Plotting Curve of Function over 5 parameters

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)

Answer 1

[Ameer Hamza]

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

[Gavin Stockstill]

Thank you!

[Ameer Hamza]

Glad to be of help.

[Gavin Stockstill]

If I needed different plots for different efficiencies (thermal and propulsive), would I need to repeat the code for each plot?

[Ameer Hamza]

Do you want to change just one parameter while keeping the other parameters constant?

[Gavin Stockstill]

I need to plot this also based on the propulsive efficiency vs Mo, along with plots for the thermal efficiency vs Mo, Fuel-to-Air ratio vs Mo, etc. I'd like it to be in a different figure.

[Ameer Hamza]

Check this. It will make 3 figures using the value of Efft, Effp, and Effo. For some reason, lines for Efft overlap. You will have a better idea of why it does not change with the value of Mc.

%Constants
T0 = 217;
gam = 1.4;
Cp = 1.004;
T4 = 1600;
hpr = 42800;
M0 = 0:1:14;
MC = 0:4; % <--- define as vector
%Equations

f1 = figure();
ax1 = axes();
hold(ax1); % to plot all lines on same axes
f2 = figure();
ax2 = axes();
hold(ax2); % to plot all lines on same axes
f3 = figure();
ax3 = axes();
hold(ax3); % 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(ax1, M0,Efft)
    plot(ax2, M0,Effp)
    plot(ax3, M0,Effo)
end

title(ax1, 'Efft');
title(ax2, 'Effp');
title(ax3, 'Effo');

[Gavin Stockstill]

That worked perfectly. Last question: How can I make each curve stop before 100% like they do in the pictures?

[Ameer Hamza]

See this. I set the threshold to 90%. You can change it by changing the lines at the end of for loop. Also, I decreased the step size for Mo so that the graph loop smooth

%Constants
T0 = 217;
gam = 1.4;
Cp = 1.004;
T4 = 1600;
hpr = 42800;
M0 = 0:0.1:14;
MC = 0:4; % <--- define as vector
%Equations

f1 = figure();
ax1 = axes();
hold(ax1); % to plot all lines on same axes
f2 = figure();
ax2 = axes();
hold(ax2); % to plot all lines on same axes
f3 = figure();
ax3 = axes();
hold(ax3); % 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;

    idx1 = Efft < 90;
    plot(ax1, M0(idx1),Efft(idx1))
    idx2 = Effp < 90;
    plot(ax2, M0(idx2),Effp(idx2))
    idx3 = Effo < 90;
    plot(ax3, M0(idx3),Effo(idx3))
end

ax1.YLim = [0 100];
ax2.YLim = [0 100];
ax3.YLim = [0 100];

title(ax1, 'Efft');
title(ax2, 'Effp');
title(ax3, 'Effo');

[Gavin Stockstill]

Thank you!

[Ameer Hamza]

Glad to be of help.

Answer 2

[Akanksha Salunkhe]

cd=0.08;  
A=0.801818;  %m/s^2
rho=1.125;   %kg/m^3
v=0:1:30;    %m/s
drag=0.5*cd*A*rho*v^2;   %N
plot(v,drag);
%why this code give error