Unformatted text preview: 9/22/2010 Problem: Aircraft are typically equipped with a static port and a Pitot tube. The static port measures the static pressure (P) of the airflow past the aircraft and the Pitot tube measures the stagnation pressure (Po) of this flow. If an aircraft's speed (expressed as a Mach number) is known, the ratio of the two pressures (Po/P) can be computed: for M <= 1: for M >= 1: PR( M ) 1 0.2M 2 PR( M ) 1.2 M 2 3.5 3.5 (7 / 6)M 2 1 / 6 2.5 where M is the aircraft's Mach number Typically the situation is the other way around. The pressure ratio is known (from the static port and Pitot tube measurements) and the Mach number must be determined. This is another root finding problem. We need to find M such that PR(M) – observed pressure ratio = 0 Function PR requires a “function m‐file”. function [ ratio ] = PR( M ) %PR Given a Mach number, returns the corresponding pressure ratio % Inputs: M = Mach number % Outputs: ratio = pressure ratio (stagnation pressure / static pressure) if M <= 1 ratio = (1 + 0.2 * M^2)^3.5; else ratio = (1.2 * M^2)^3.5 * ((7/6)* M^2 ‐ 1/6)^‐2.5; end end Such functions are created in the editor window and each function is stored in a file called function.m, where function is the name of the function (e.g. for this function, the file is PR.m). 1 9/22/2010 Two kinds of m‐files: • script files : contain “canned” instructions • function files: contain a function (file name = function name) Two kinds of functions (at least for now): • “anonymous” functions: defined by a single command • file functions: generally more complex, defined within a file File functions are comparable to functions in C++ (and other languages) • self‐contained units • do NOT share variables with the command window • accept inputs and produce outputs The initial comments document the function: >> lookfor mach TargetsComms_Machine ‐ class to hold information about a machine PR ‐ Given a Mach number, returns the corresponding pressure ratio PRv ‐ Given a Mach number, returns the corresponding pressure ratio sfopen ‐ Opens a new machine. sfsave ‐ Saves a machine. wmachdep ‐ Machine dependent values. >> help PR PR Given a Mach number, returns the corresponding pressure ratio Inputs: M = mach number Outputs: ratio = pressure ratio (stagnation pressure / static pressure) The lookfor command looks for functions having the specified keyword in their initial comment line. The help command outputs all of a function’s initial comments. 2 9/22/2010 The Matlab “if” is logically equivalent to a C++ “if”. Only the details differ. if M <= 1 % round brackets around expression allowed but not required ratio = (1 + 0.2 * M^2)^3.5; else ratio = (1.2 * M^2)^3.5 * ((7/6)* M^2 ‐ 1/6)^‐2.5; end Relational and logical operators: == ~= < > <= >= ~ && || is equal to is not equal to ** different from C++ ** is less than is greater than is less than or equal to is greater than or equal to NOT ** different from C++ ** AND note: & is array version OR note: | is array version The script file below plots the pressure ratio for Mach numbers from 0 to 7 (world record: the X‐15 rocket plane achieved Mach 6.72), reads in a pressure ratio, and outputs the corresponding Mach number. figure (1); fplot (@PR, [0, 7]); % Note the [email protected] xlabel ('Mach Number'); ylabel ('Pressure Ratio'); grid on; ratio = input ('Enter pressure ratio: '); % define function for root finding f = @(M) PR(M) ‐ ratio; M = fzero (f, [0 7]); fprintf ('The mach number is %f\n', M); 3 9/22/2010 Functions like “fplot” expect to be given a function handle. “Anonymous” functions: The function itself is anonymous What looks like a function name is actually the name of a “function handle” >> f = @(x) x + 4; % “f” is a function handle >>fplot (f, [6 20]); % which is exactly what “fplot” expects File functions: The name relates to the function itself When a function handle is required, one must be created by using [email protected] >>fplot (@PR, [0, 7]); % [email protected] creates a function handle Built‐in functions: These are also named >>fplot (@sin, [0, 7]); % a function handle must be created Function PR is not vector friendly. Unlike most Matlab functions, it will not work properly if it is given a vector of input values. In order to convert a vector of inputs into a vector of outputs it is necessary to use a loop: M = linspace (0, 7, 100); % 100 points between 0 and 7 for k = 1 : length(M) % from 1 to length(M) in steps of 1 ratio(k) = PR(M(k)); end figure (1) plot (M, ratio); xlabel ('Mach Number'); ylabel ('Pressure Ratio'); grid on; This code takes advantage of the fact that vectors grow automatically when a value is assigned to a previously non‐existent value. 4 9/22/2010 Function “length” returns the length of a row or column vector. Function “size” returns a two element row vector containing the dimensions of a variable (number of rows and numbers of columns). >> A = [1 3 5 7]; >> length(A) ans = 4 % A is 4 elements long >> size(A) ans = 1 4 % A has 1 row and 4 columns >> x = 4; >> size(x) ans = 1 1 % Matlab treats scalars as 1x1 arrays Having vectors grow one element at a time is not very efficient. It is much better create a vector of the required size and then fill in the values. This is called “preallocation”. M = linspace (0, 7, 100); % 100 points between 0 and 7 ratio = ones(size(M)); % preallocate vector for efficiency for k = 1 : length(M) % from 1 to length(M) in steps of 1 ratio(k) = PR(M(k)); end Function “ones” creates an array full of ones. It accepts a two element row vector containing the dimensions (# of rows, # of columns) of the array. >> ones([2 3]) ans = 1 1 1 1 1 1 Function “zeros” (note that Americans don’t spell properly) is analogous but creates an array of zeroes. 5 9/22/2010 The code necessary to handle vectors can be incorporated in the function. Function “PRv” does the same job as “PR” did but works properly when the input value is a vector. function [ ratio ] = PRv( M ) %PRv Given a Mach number, returns the corresponding pressure ratio % Vector friendly version of PR % Inputs: M = Mach number (may be a vector of Mach numbers) % Outputs: ratio = pressure ratio (stagnation pressure / static pressure) ratio = ones(size(M)); % preallocate for efficiency for k = 1:length(M) if M(k) <= 1 ratio(k) = (1 + 0.2 * M(k)^2)^3.5; else ratio(k) = (1.2 * M(k)^2)^3.5 * ((7/6)* M(k)^2 ‐ 1/6)^‐2.5; end end end The Matlab “for” loop is actually for var = array ... end The loop is executed once for each value in the array. “val1:val2 “ creates an array of values from “val1” to “val2” in steps of 1: >> A = 1:4 A = 1 2 3 4 A step size may be specified: >> A = 1:3:10 A = 1 4 7 10 The series stops at the last value that does not exceed the upper limit: >> A = 1:3:9 A = 1 4 7 6 9/22/2010 Summary of ways of creating vectors: >> v1 = [2 3 4 5 6] % the values may be separated with commas if desired v1 = 2 3 4 5 6 >> v2 = linspace (2, 6, 5) % 5 evenly spaced values between 2 and 6 v2 = 2 3 4 5 6 >> v3 = 2 : 6 % from 2 to 6 in steps of 1 (the default) v3 = 2 3 4 5 6 >> v4 = 1 : 2 : 7 % from 1 to 7 in steps of 2 v4 = 1 3 5 7 >> v5 = ones([1 5]) % row vector of ones v5 = 1 1 1 1 1 >> v6 = zeros([1 5]) % row vector of zeroes v6 = 0 0 0 0 0 >> v7 = ones(size(v5)) % size specified using function “size” v7 = 1 1 1 1 1 Matlab also has while loops (just like C++). The script below uses an infinite loop to process pressure ratios until a ratio that is less than 1 is entered. while true % true is 1, false is 0 ratio = input ('Enter pressure ratio: '); if ratio < 1; break; end % could be return instead of break % define function for root finding f = @(M) PR(M) ‐ ratio; M = fzero (f, [0 7]); fprintf ('The mach number is %f\n', M); end 7 9/22/2010 Perhaps what we would really like is a function that converts pressure ratios to Mach numbers… function [ M ] = PR2Mach( ratio ) %PR2Mach Given a pressure ratio, returns the corresponding Mach number % Inputs: ratio = pressure ratio (stagnation pressure / static % pressure) % Outputs: M = Mach number if ratio < PR(0) || ratio > PR(7) error 'Pressure ratio is not reasonable.'; % abort function execution end % define function for root finding f = @(M) PR(M) ‐ ratio; % find root M = fzero (f, [0 7]); end 8 ...
View Full Document
This note was uploaded on 10/18/2010 for the course ECOR2606 2606 taught by Professor Aaaa during the Spring '10 term at Carleton CA.
- Spring '10