I am trying to plot multiple data sets on one one diagram. When I use the plot function, I am not getting any data plots.

I am a novice to Matlab. My goal for this code is to calculate exit Mach numbers and total pressure ratios using a set of equations. The following code is what I generated based on my knowledge (I am sure there is a better way) but for some reason, I am not getting any data plots. Code:

clc

clear all

M1 = [1:0.1:2.5]; %%Free Stream Mach

a = [-2.5:0.5:2.5]; %%Angles of Attack

gamma = 1.4;

Beta1 = [5 8 12]; %%Design 1

Beta2 = [7 8 10];%%Design 2

Beta3 = [5 10 10];%%Design 3

Beta4 = [8.33 8.33 8.33];%%Design 4

Beta5 = [3 5 17];%%Design 5

%%normal shock relation

for i=1:length(M1)

M2 = sqrt((1+((gamma-1)./2)*(M1(i).^2))./((gamma.*(M1(i).^2))-((gamma-1)./2)))

p2p1 = (1+((2.*gamma)./(gamma+1))*((M1(i).^2)-1))

p02p2 = (1+((gamma-1)./2).*M2.^2).^(gamma./(gamma-1))

p1p01 = (1+((gamma-1)./2).*M1(i).^2).^-(gamma./(gamma-1))

p02p01 = (p2p1).*(p02p2).*(p1p01)

hold on

end

%%Oblique shock relation

for i=1:length(M1)

for j=1:length(Beta1)

theta1 = atand(2.*cotd(Beta1(j)).*((((M1(i).^2).*(sin(Beta1(j))).^2)-1)./((M1(i).^2).*(gamma+cos(2.*Beta1(j)))+2)))

Mn1 = M1(i).*sind(Beta1(j));

Mn2 = sqrt((Mn1.^2)+(2./(gamma-1))./(((2.*gamma)./(gamma-1)).*(Mn1.^2)-1));

M2 = Mn2./(sind(Beta1(j)-theta1));

p2p1 = (1+((2.*gamma)./(gamma+1)).*((Mn1.^2)-1))

p02p2 = (1+((gamma-1)./2).*Mn2.^2).^(gamma./(gamma-1))

p1p01 = (1+((gamma-1)./2).*Mn1.^2).^-(gamma./(gamma-1))

p02p01 = (p2p1).*(p02p2).*(p1p01)

hold on

end

for j=1:length(Beta2)

theta2 = atand(2.*cotd(Beta2(j)).*((((M1(i).^2).*(sin(Beta2(j))).^2)-1)./((M1(i).^2).*(gamma+cos(2.*Beta2(j)))+2)))

Mn1 = M1(i).*sind(Beta2(j));

Mn2 = sqrt((Mn1.^2)+(2./(gamma-1))./(((2.*gamma)./(gamma-1)).*(Mn1.^2)-1));

M2 = Mn2./(sind(Beta2(j)-theta2));

p2p1 = (1+((2.*gamma)./(gamma+1)).*((Mn1.^2)-1))

p02p2 = (1+((gamma-1)./2).*Mn2.^2).^(gamma./(gamma-1))

p1p01 = (1+((gamma-1)./2).*Mn1.^2).^-(gamma./(gamma-1))

p02p01 = (p2p1).*(p02p2).*(p1p01)

hold on

end

for j=1:length(Beta3)

theta3 = atand(2.*cotd(Beta3(j)).*((((M1(i).^2).*(sin(Beta3(j))).^2)-1)./((M1(i).^2).*(gamma+cos(2.*Beta3(j)))+2)))

Mn1 = M1(i).*sind(Beta3(j));

Mn2 = sqrt((Mn1.^2)+(2./(gamma-1))./(((2.*gamma)./(gamma-1)).*(Mn1.^2)-1));

M2 = Mn2./(sind(Beta3(j)-theta3));

p2p1 = (1+((2.*gamma)./(gamma+1)).*((Mn1.^2)-1))

p02p2 = (1+((gamma-1)./2).*Mn2.^2).^(gamma./(gamma-1))

p1p01 = (1+((gamma-1)./2).*Mn1.^2).^-(gamma./(gamma-1))

p02p01 = (p2p1).*(p02p2).*(p1p01)

hold on

end

for j=1:length(Beta4)

theta4 = atand(2.*cotd(Beta4(j)).*((((M1(i).^2).*(sin(Beta4(j))).^2)-1)./((M1(i).^2).*(gamma+cos(2.*Beta4(j)))+2)))

Mn1 = M1(i).*sind(Beta4(j));

Mn2 = sqrt((Mn1.^2)+(2./(gamma-1))./(((2.*gamma)./(gamma-1)).*(Mn1.^2)-1));

M2 = Mn2./(sind(Beta4(j)-theta4));

p2p1 = (1+((2.*gamma)./(gamma+1)).*((Mn1.^2)-1))

p02p2 = (1+((gamma-1)./2).*Mn2.^2).^(gamma./(gamma-1))

p1p01 = (1+((gamma-1)./2).*Mn1.^2).^-(gamma./(gamma-1))

p02p01 = (p2p1).*(p02p2).*(p1p01)

hold on

end

for j=1:length(Beta5)

theta5 = atand(2.*cotd(Beta5(j))*((((M1(i).^2).*(sin(Beta5(j))).^2)-1)./((M1(i).^2).*(gamma+cos(2.*Beta5(j)))+2)))

Mn1 = M1(i).*sind(Beta5(j));

Mn2 = sqrt((Mn1.^2)+(2./(gamma-1))./(((2.*gamma)./(gamma-1)).*(Mn1.^2)-1));

M2 = Mn2./(sind(Beta5(j)-theta5));

p2p1 = (1+((2.*gamma)./(gamma+1)).*((Mn1.^2)-1))

p02p2 = (1+((gamma-1)./2).*Mn2.^2).^(gamma./(gamma-1))

p1p01 = (1+((gamma-1)./2).*Mn1.^2).^-(gamma./(gamma-1))

p02p01 = (p2p1).*(p02p2).*(p1p01)

hold on

end

end

plot(M1,p02p01, 'LineWidth', 4);

Answer 1

[DGM]

p02p01 is a scalar, so there's not really anything to plot. I imagine that's supposed to be a vector, but you're just overwriting the same value in the loops. If that's the case, you'll have to decide which of the intermediate calculations are to be saved and then save them to a vector.

[Jeffrey Lewis]

I figured that was the case but couldn't figure out how to identify the issue. So, I should erase the subsequent for loops and just have one for loop that iterates over the range of Machs and Betas to get my pressure ratios?

I guess my biggest question about the logic is, how do I get the loop to consider each Beta case?

Thanks for your input thus far!

[DGM]

Just what you want to store and how you want to store it. Instead of writing to the same scalar every time, write it to an array of some sort.

M1 = [1:0.1:2.5]; %%Free Stream Mach 
a = [-2.5:0.5:2.5]; %%Angles of Attack
gamma = 1.4;
Beta1 = [5 8 12]; %%Design 1
Beta2 = [7 8 10];%%Design 2
Beta3 = [5 10 10];%%Design 3
Beta4 = [8.33 8.33 8.33];%%Design 4
Beta5 = [3 5 17];%%Design 5

%%normal shock relation
p02p01 = zeros(size(M1));
for i=1:length(M1)
    M2 = sqrt((1+((gamma-1)./2)*(M1(i).^2))./((gamma.*(M1(i).^2))-((gamma-1)./2)));
    p2p1 = (1+((2.*gamma)./(gamma+1))*((M1(i).^2)-1));
    p02p2 = (1+((gamma-1)./2).*M2.^2).^(gamma./(gamma-1));
    p1p01 = (1+((gamma-1)./2).*M1(i).^2).^-(gamma./(gamma-1));
    p02p01(i) = (p2p1).*(p02p2).*(p1p01);
end

plot(M1,p02p01,'LineWidth',2);

%%Oblique shock relation
% assuming all Beta vectors are the same length
allP = zeros(numel(Beta1),numel(M1),5); 
for i=1:length(M1)
    for j=1:length(Beta1)
        theta1 = atand(2.*cotd(Beta1(j)).*((((M1(i).^2).*(sin(Beta1(j))).^2)-1)./((M1(i).^2).*(gamma+cos(2.*Beta1(j)))+2)));
        Mn1 = M1(i).*sind(Beta1(j));
        Mn2 = sqrt((Mn1.^2)+(2./(gamma-1))./(((2.*gamma)./(gamma-1)).*(Mn1.^2)-1));
        M2 = Mn2./(sind(Beta1(j)-theta1));
        p2p1 = (1+((2.*gamma)./(gamma+1)).*((Mn1.^2)-1));
        p02p2 = (1+((gamma-1)./2).*Mn2.^2).^(gamma./(gamma-1));
        p1p01 = (1+((gamma-1)./2).*Mn1.^2).^-(gamma./(gamma-1));
        allP(j,i,1) = (p2p1).*(p02p2).*(p1p01);
    end
    for j=1:length(Beta2)
        theta2 = atand(2.*cotd(Beta2(j)).*((((M1(i).^2).*(sin(Beta2(j))).^2)-1)./((M1(i).^2).*(gamma+cos(2.*Beta2(j)))+2)));
        Mn1 = M1(i).*sind(Beta2(j));
        Mn2 = sqrt((Mn1.^2)+(2./(gamma-1))./(((2.*gamma)./(gamma-1)).*(Mn1.^2)-1));
        M2 = Mn2./(sind(Beta2(j)-theta2));
        p2p1 = (1+((2.*gamma)./(gamma+1)).*((Mn1.^2)-1));
        p02p2 = (1+((gamma-1)./2).*Mn2.^2).^(gamma./(gamma-1));
        p1p01 = (1+((gamma-1)./2).*Mn1.^2).^-(gamma./(gamma-1));
        allP(j,i,2) = (p2p1).*(p02p2).*(p1p01);
    end
    for j=1:length(Beta3)
        theta3 = atand(2.*cotd(Beta3(j)).*((((M1(i).^2).*(sin(Beta3(j))).^2)-1)./((M1(i).^2).*(gamma+cos(2.*Beta3(j)))+2)));
        Mn1 = M1(i).*sind(Beta3(j));
        Mn2 = sqrt((Mn1.^2)+(2./(gamma-1))./(((2.*gamma)./(gamma-1)).*(Mn1.^2)-1));
        M2 = Mn2./(sind(Beta3(j)-theta3));
        p2p1 = (1+((2.*gamma)./(gamma+1)).*((Mn1.^2)-1));
        p02p2 = (1+((gamma-1)./2).*Mn2.^2).^(gamma./(gamma-1));
        p1p01 = (1+((gamma-1)./2).*Mn1.^2).^-(gamma./(gamma-1));
        allP(j,i,3) = (p2p1).*(p02p2).*(p1p01);
    end
    for j=1:length(Beta4)
        theta4 = atand(2.*cotd(Beta4(j)).*((((M1(i).^2).*(sin(Beta4(j))).^2)-1)./((M1(i).^2).*(gamma+cos(2.*Beta4(j)))+2)));
        Mn1 = M1(i).*sind(Beta4(j));
        Mn2 = sqrt((Mn1.^2)+(2./(gamma-1))./(((2.*gamma)./(gamma-1)).*(Mn1.^2)-1));
        M2 = Mn2./(sind(Beta4(j)-theta4));
        p2p1 = (1+((2.*gamma)./(gamma+1)).*((Mn1.^2)-1));
        p02p2 = (1+((gamma-1)./2).*Mn2.^2).^(gamma./(gamma-1));
        p1p01 = (1+((gamma-1)./2).*Mn1.^2).^-(gamma./(gamma-1));
        allP(j,i,4) = (p2p1).*(p02p2).*(p1p01);
    end
    for j=1:length(Beta5)
        theta5 = atand(2.*cotd(Beta5(j))*((((M1(i).^2).*(sin(Beta5(j))).^2)-1)./((M1(i).^2).*(gamma+cos(2.*Beta5(j)))+2)));
        Mn1 = M1(i).*sind(Beta5(j));
        Mn2 = sqrt((Mn1.^2)+(2./(gamma-1))./(((2.*gamma)./(gamma-1)).*(Mn1.^2)-1));
        M2 = Mn2./(sind(Beta5(j)-theta5));
        p2p1 = (1+((2.*gamma)./(gamma+1)).*((Mn1.^2)-1));
        p02p2 = (1+((gamma-1)./2).*Mn2.^2).^(gamma./(gamma-1));
        p1p01 = (1+((gamma-1)./2).*Mn1.^2).^-(gamma./(gamma-1));
        allP(j,i,5) = (p2p1).*(p02p2).*(p1p01);
    end
end

figure
subplot(5,1,1)
plot(M1,allP(:,:,1),'LineWidth',1);

Warning: Imaginary parts of complex X and/or Y arguments ignored.

subplot(5,1,2)
plot(M1,allP(:,:,2),'LineWidth',1);

Warning: Imaginary parts of complex X and/or Y arguments ignored.

subplot(5,1,3)
plot(M1,allP(:,:,3),'LineWidth',1);

Warning: Imaginary parts of complex X and/or Y arguments ignored.

subplot(5,1,4)
plot(M1,allP(:,:,4),'LineWidth',1);

Warning: Imaginary parts of complex X and/or Y arguments ignored.

subplot(5,1,5)
plot(M1,allP(:,:,5),'LineWidth',1);

Warning: Imaginary parts of complex X and/or Y arguments ignored.

Note that the results in allP are complex-valued. I don't know if that's expected.

[Jeffrey Lewis]

This is great help! This makes more sense. The only other question I had is what is the syntax needed to have all data sets on the same plot. I.e. if I wanted the same trend lines from the oblique shocks to compare to the Normal shock and, say, have those lines in different colors?

Could I do a 'hold on' for each data set and do my plot function at the end? Sorry if this is super novice questions but, that's exactly what my level is.

Also, I have double checked my equations, I imagine that the complex numbers are based on the betas chosen so I will run with it and talk to any errors that are questioned.

Thanks again for the help.

[DGM]

Yes, you should be able to work the first plot into each of the later plots if you want. You might need to play with the line properties, etc. to make them distinguishable.

figure
subplot(5,1,1)
% plot p02p01 in a heavy black line
plot(M1,p02p01,'k','LineWidth',2); hold on
% plot others in a lighter weight line, picking colors 
% from the axes 'colororder' property
plot(M1,allP(:,:,1),'LineWidth',1);

subplot(5,1,2)
% ... and so on