1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
| alt,vit=[0],[0]
z2,y,x = 0,0,0.6371009*(10**(-19))
w=2*np.pi/(23*3600+59*60+4) #rad/s#vitesse de rotation de la terre(465,1 m/s )
for i in range (l-3):
if z2<100000:
k1 = -g + (pousse[i] - 0.5*Cx*maitre_couple(tps[i])*masse_volumique(z2)*v**2)/masse[i] + x*w**2
else:
k1 = -g + pousse[i]/masse[i] + x*w**2
if z2< 100000:
k2 = -g + ((pousse[i]+pousse[i+1])/2 - 0.5*Cx*(np.pi*(5.425/2)**2 + 2*np.pi*(3.15/2)**2)*masse_volumique(z2)*((y+(pas[i]*k1/2))**2))/(masse[i]+masse[i+1])/2 + x*w**2
else:
k2 = -g +(pousse[i]+pousse[i+1])/2 + x*w**2
if z2< 100000:
k3 = -g + ((pousse[i]+pousse[i+1])/2 -0.5*Cx*(np.pi*(5.425/2)**2 + 2*np.pi*(3.15/2)**2)*masse_volumique(z2)*(y+1/2*pas[i]*k2)**2)/(masse[i]+masse[i+1])/2 + x*w**2
else:
k3 = -g + (pousse[i]+pousse[i+1])/2 + x*w**2
if z2<100000:
k4 = -g + ((pousse[i+1]+(y+pas[i]*k3)*debitm[i+1])/masse[i+1]) + x*w**2
else:
k4 = -g + (poussee[i+1]/masse[i+1]) + x*w**2
y+=pas[i]*(k1/6 + k2/3 + k3/3 + k4/6)
z2+=y*pas[i]
x+=y*pas[i]
alt.append(z2)
vit.append(y)
#plt.plot(tps[0:683],alt, label='altirkf')
plt.plot(tps[0:684], vit)
plt.legend
plt.show() |
Partager