MATLAB: Plotting Curve of Function over 5 parameters

MATLABplotting

I am trying to find out how to make a curve like the one in the attached picture where n is plotted as a function of Mo with Mc as a parameter. Basically, I am trying to put all 5 curves for Mc (0 to 4) on the plot for Mo vs n. Here is the code I have so far. I currently have Mc = 4 and the curve is working, I just am not sure how to plot all 5 Mc curves on one graph.
%Constants
T0 = 217;
gam = 1.4;
Cp = 1.004;
T4 = 1600;
hpr = 42800;
M0 = 0:1:14;
Mc = 4;
%Equations
% Gas Constant
R = ((gam-1)/ gam)*Cp;
%Velocity of Sound
a0 = (1000*gam*R*T0).^(1/2);
%Total to Static Temperature Ratio
Taur = 1 + (((gam-1)/2)*(M0.^2));
%Burner Exit Total Temp, T'max
Tpmax = T4*(1+((gam-1)/2)*Mc.^2);
%Burner Exit Enthalpy to Ambient Enthalpy Ratio
Taugam = Tpmax/T0;
%Velocity Ratio
VR = M0.*((Taugam./Taur).^(1/2));
%Exit Velocity
v9 = a0*VR;
%Temperature at Nozzle Exit
T9 = (((v9/a0).^2)/M0.^2)*T0;
%Burner Exit Temp to Burner Inlet Temp
Taub = T9/T0;
%Temp at station t2
Tt2 = Tpmax/Taub;
%Specific Thrust
ST = (a0.*M0).*((sqrt(Taugam./Taur)-1));
%Fuel to Air Ratio
f = ((Cp*Tt2)/hpr)*((Tpmax/Tt2)-1);
%Thrust Specific Fuel Consumption
TSFC = (Cp*T0*(Taugam-Taur))./(a0.*M0.*hpr.*(sqrt(Taugam./Taur)-1));
%Thermal Efficiency
Efft = (1-(1./Taur))*100;
%Propulsive Efficiency
Effp = (2./(((Taugam./Taur).^(1/2))+1))*100;
%Overall Efficiency
Effo = (Efft.*Effp)/100;
%Total to Static Temperature Ratio for Maximum Specific Thrust
Taurmax = Taugam.^(1/3);
%Mach Number of Flight for Maximum Specific Thrust
M0max = ((2/(gam-1))*((T4/T0)*(1+((gam-1)/2)*Mc.^2)).^1/3).^1/2;
%Temperature at Station Tt0
Tt0 = Taur*T0;
plot(M0,Effo)

Best Answer

See this code
%Constants
T0 = 217;
gam = 1.4;
Cp = 1.004;
T4 = 1600;
hpr = 42800;
M0 = 0:1:14;
MC = 0:4; % <--- define as vector
%Equations
f = figure();
ax = gca;
hold(ax); % to plot all lines on same axes
for i=1:numel(MC)
    Mc = MC(i); % select one value of Mc
    % Gas Constant
    R = ((gam-1)/ gam)*Cp;
    %Velocity of Sound
    a0 = (1000*gam*R*T0).^(1/2);
    %Total to Static Temperature Ratio
    Taur = 1 + (((gam-1)/2)*(M0.^2));
    %Burner Exit Total Temp, T'max
    Tpmax = T4*(1+((gam-1)/2)*Mc.^2);
    %Burner Exit Enthalpy to Ambient Enthalpy Ratio
    Taugam = Tpmax/T0;
    %Velocity Ratio
    VR = M0.*((Taugam./Taur).^(1/2));
    %Exit Velocity
    v9 = a0*VR;
    %Temperature at Nozzle Exit
    T9 = (((v9/a0).^2)/M0.^2)*T0;
    %Burner Exit Temp to Burner Inlet Temp
    Taub = T9/T0;
    %Temp at station t2
    Tt2 = Tpmax/Taub;
    %Specific Thrust
    ST = (a0.*M0).*((sqrt(Taugam./Taur)-1));
    %Fuel to Air Ratio
    f = ((Cp*Tt2)/hpr)*((Tpmax/Tt2)-1);
    %Thrust Specific Fuel Consumption
    TSFC = (Cp*T0*(Taugam-Taur))./(a0.*M0.*hpr.*(sqrt(Taugam./Taur)-1));
    %Thermal Efficiency
    Efft = (1-(1./Taur))*100;
    %Propulsive Efficiency
    Effp = (2./(((Taugam./Taur).^(1/2))+1))*100;
    %Overall Efficiency
    Effo = (Efft.*Effp)/100;
    %Total to Static Temperature Ratio for Maximum Specific Thrust
    Taurmax = Taugam.^(1/3);
    %Mach Number of Flight for Maximum Specific Thrust
    M0max = ((2/(gam-1))*((T4/T0)*(1+((gam-1)/2)*Mc.^2)).^1/3).^1/2;
    %Temperature at Station Tt0
    Tt0 = Taur*T0;
    plot(M0,Effo)
end
Related Question