Mach Number Area Relation Several Inputs

The following code inputs a single input (ARatio) and outputs two terms (Msub) & (Msup)

I want to be able to input several ARatio values, more specifically an array from 0.1 to 10 with an increment of 0.1

and have Msub & Msup for each of these ARatio stored in an array in order to plot ARatio and Mach number.

clear;
clc;

%% INPUTS

% Define some paramters
g   = 1.4;
gm1 = g-1;
gp1 = g+1;

% Define anonymous function with two inputs (M and ARatio)
% - Will be used in the methods below
% - Pass M and ARatio as arguments to AM_EQN to get function value
% -     funVal = AM_EQN(M,ARatio)
AM_EQN = @(M,ARatio) sqrt((1/M^2)*(((2+gm1*M^2)/gp1)^...
                            (gp1/gm1)))-ARatio;

% Solve for Msub and Msup using this area ratio (A/A*)
ARatio = 1.5;

% Error tolerance
errTol = 1e-4;

% Flags for printing iterations to screen
verboseBisection   = 0;
verboseIncremental = 0;


%% MATLAB SOLVER

% Set up the solver
problem.objective = @(M) (1/M^2)*(((2+gm1*M^2)/gp1)^(gp1/gm1))-ARatio^2;    % Objective function
problem.solver    = 'fzero';                                                % Find the zero
problem.options   = optimset(@fzero);                                       % Default options

% Solve subsonic root
problem.x0 = [1e-6 1];                                                      % Subsonic solver bounds
Msub       = fzero(problem);                                                % Solve for subsonic M
problem.x0 = [1+1e-6 50];                                                   % Supersonic solver bounds
Msup       = fzero(problem);                                                % Solve for supersonic M

% Print solutions to command window
fprintf('==== MATLAB SOLVER ====\n');
fprintf('Msub: %3.4f\n',Msub);
fprintf('Msup: %3.4f\n',Msup);
fprintf('=======================\n\n');

[Jim Riggs]

I am curious to the origin of this function equation.

[Steven Castrillon]

if you set the function i gave equal to 0 and move AR to the other side of the equation you will get the equation for Area Ratio - Mach Number Relation.

[Star Strider]

Related: Area Mach Number Relation

[Jim Riggs]

Thank you. It works.

Answer 1

[David Hill]

function [Msub,Msup] = sub_super(ARatio)%ARatio needs to be >1
g   = 1.4;
gm1 = g-1;
gp1 = g+1;
Msub=zeros(size(ARatio));
Msup=zeros(size(ARatio));
for i=1:length(ARatio)
    problem.objective = @(M) (1/M^2)*(((2+gm1*M^2)/gp1)^(gp1/gm1))-ARatio(i)^2;    
    problem.solver    = 'fzero';                                              
    problem.options   = optimset(@fzero);                                      
    problem.x0 = [1e-6 1];                                                     
    Msub(i)       = fzero(problem);                                               
    problem.x0 = [1+1e-6 50];                                                   
    Msup(i)       = fzero(problem);  
end
plot(ARatio,Msub);
hold on
plot(ARatio,Msup);
end