Parametric_plot gives half-broken picture for transcendental function

edit

• Goal: Solve 2 functions about Ma1, delta, beta, sigma, and give a contour plot using Ma1 & sigma as axes. (an Aerodynamic problem)
• Problem: One of them is transcendental. I tried:

1. For every give delta & Ma1, create a loop to iterate beta, and numerically solve sigma for every point. But find_root gives Runtime Error, and interrupted the script.
2. Solve the equations for Ma1 & sigma under given delta, and use parametric_plot to give the picture by iterating beta. This works half. The other half of these lines are missed.
3. Substitute delta and use implicit_plot. This works on plotting eq1 using Ma1 & beta as axes. But I have no idea of introducing 3 variables in order to give a 2D plot.

I think it's not a numerically unsolvable problem since my classmates got a picture using MATLAB. But I couldn't figure out what's wrong with sagemath. Maybe it's because of the order of solving equations? Or, maybe it's the limit of sagemath? Please give me a little help, thanks.

Script of attempt 2:

from sage.plot.colors import rainbow
clr=rainbow(8);
#Shock Wave 4
delta,Ma1,Ma2,gamma,beta,pp,sigma=SR.var('delta Ma1 Ma2 gamma beta pp sigma')
eq1=tan(delta)==(Ma1^2*sin(beta)^2-1)/(tan(beta)*(Ma1^2*((gamma+1)/2-sin(beta)^2)+1))
eq1=eq1.subs(gamma=1.4)
eq4=sigma==(2*gamma/(gamma+1)*Ma1^2*sin(beta)^2-(gamma-1)/(gamma+1))^(-1/(gamma-1))*((gamma+1)*Ma1^2*sin(beta)^2/((gamma-1)*Ma1^2*sin(beta)^2+2))^(gamma/(gamma-1))
eq4=eq4.subs(gamma=1.4)
sol4=solve([eq1,eq4],[Ma1,sigma],solution_dict=True)
sud=0;bb=0.05
fg4=parametric_plot((sol4[1][Ma1].subs(delta=sud),sol4[1][sigma].subs(delta=sud)),(beta,0.1,pi/4-bb),color=clr[0],
              axes_labels=(r'$Ma_1$',r'$\sigma$'),title=r'$\sigma\ vs\ Ma_1$')
fg4+=text(r'$\delta=0^\circ$',(1.1,0.9),color=clr[0])
for i in range(1,8):
    sud=i*pi/36
    fg4+=parametric_plot((sol4[1][Ma1].subs(delta=sud),sol4[1][sigma].subs(delta=sud)),(beta,0.1,pi/4-bb),color=clr[i])
    fg4+=parametric_plot((sol4[0][Ma1].subs(delta=sud),sol4[0][sigma].subs(delta=sud)),(beta,0.1,pi/4-bb),color=clr[i])
    tx=sol4[1][Ma1].subs(delta=sud,beta=pi/2.4)+0.05;
    ty=sol4[1][sigma].subs(delta=sud,beta=pi/2.4)+0.02;
    delt='$'+str(i*5)+'^\\circ$'
    fg4+=text(r'%s'%delt,(tx,ty),color=clr[i])
end
fg4.show(xmin=1,xmax=4,ymin=0,ymax=1,aspect_ratio=2)