I wrote a code to replicate the image to the right but for some reason the M3 and the M2 are on opposite sides and the M1 is completely wrong. Please help me fix this

clc
clear all
%Given:
M1 = 3;
G = 1.4;
tot = 180;
Theta= zeros(1,tot);
p2_p1= zeros(1,tot);


n=1;


for i = 1 : 0.01 : tot

    Theta(n)  = atand(2*cotd(i)*((M1*(sind(i)))^2-1)/(2+M1^2*(G+cosd(2*i))));
    p2_p1(n) = 1+(2*G)/((G+1)*(M1*sind(i))^2-1);

    n = n+1;
end

error = 1;

i = 0;

while error > 0.001
    i=i+1;
    error = 20 - Theta(i);

end

d=i;

b2 = 38;
Mn1 = M1*sind(b2);
Mn2 = sqrt((Mn1^(2)+(2/(G-1)))/(2*G/(G-1)*Mn1^(2)-1));    
M2 = Mn2/(sind(b2-20));

n = 1;

for i = 30 : 0.01 : 149

    theta2(n)  = 20 + atand(2*cotd(i)*((M2*(sind(i)))^2-1)/(2+M2^2*(G+cosd(2*i))));
    p3_p1(n) = p2_p1(d)*(1+(2*G)/(G+1)*((M2*(sind(i)))^2-1));

    n = n+1;
end

error = 1;
i = 0;
while error > 0.001
    i = i+1;
    error = abs(-15-Theta(i));

end

a=i;

b3 = 32;
Mn1 = M1*sind(b3);
Mn3 = sqrt((Mn1^(2)+(2/(G-1)))/(2*G/(G-1)*Mn1^(2)-1));    
M3 = Mn3/sind(b3-15);

s = 1;

for i = 26 : 0.01 : 153

    theta3(s)  = -15 + atand(2*cotd(i)*((M3*(sind(i)))^2-1)/(2+M3^2*(G+cosd(2*i))));
    p4_p1(s) = p2_p1(a)*(1+(2*G)/(G+1)*((M3*(sind(i)))^2-1));

    s = s+1;
end

figure(1)
plot(Theta,p2_p1)
hold on
plot(theta2,p3_p1)
hold on
plot(theta3,p4_p1)

title('PD-Diagram')
xlabel('Theta (\theta)')
ylabel('Pressure Ratio')
axis([-45 45 1 17])
legend('M1','M2','M3')

Answer 1

[Milan Bansal]

Hi Aijalon Marsh,

It is my understanding that you are trying to plot the graphs M1, M2 and M3 present in the attached image but the M1 curve is not plotted as expected and M2 and M3 are on the opposite side.

For correcting M1 curve, please modify the "for loop" for M1 in the code as showm in the code snippet below.

% iterate from 20 to 160 instead of 1 to 180
for i = 20 : 0.01 : 160
    Theta(n)  = atand(2*cotd(i)*((M1*(sind(i)))^2-1)/(2+M1^2*(G+cosd(2*i))));
    p2_p1(n) = 1+(2*G)/(G+1)*((M1*(sind(i)))^2-1); % modified equation
               
    n = n+1;
end

The output firgure after the modification is shown below.

The M2 and M3 plots are at the correct positions and M1 plot is in the correct form. Please modify the "M1", "M2" and "M3" variables to change the size of plots and the threshold values to set the positions to replicate the plots given in the image.

Hope it helps!