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 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
|
import numpy as np
import matplotlib.pyplot as plt # importations des modules numpy et matplotlib.pyplot
#données initiales
xM=0 # variable position x de Mars
yM=0 # variable position y de Mars
xT=0
yT=0
Ms= 2*10**30
G= 6.6743* 10**(-11) # constante de gravitation
Rm=3389.5 * 10**(3) # rayon de Mars
Rt=6371* 10**3 # Rayon de la Terre
Mt=5.972* 10**(24) # Masse de la Terre
Mm=6.42* 10**(23) # Masse de Mars
dSM=207* 10**9 # distance Soleil-Mars
dST=150* 10**9 # distance Soleil-Terre
phi=0 # variable angle de déphasage de Mars
def rk4(x, dx, y, deriv):
# /* d1, d2, d3, d4 = estimations des derivees
# yp = estimations intermediaires des fonctions */
ddx = dx/2. # /* demi-pas */
d1 = deriv(x,y) # /* 1ere estimation */
yp = y + d1*ddx
d2 = deriv(x+ddx,yp) #/* 2eme estimat. (1/2 pas) */
yp = y + d2*ddx
d3 = deriv(x+ddx,yp) #/* 3eme estimat. (1/2 pas) */
yp = y + d3*dx
d4 = deriv(x+dx,yp) # /* 4eme estimat. (1 pas) */
#/* estimation de y pour le pas suivant en utilisant
# une moyenne ponderee des derivees en remarquant
# que : 1/6 + 1/3 + 1/3 + 1/6 = 1 */
return y + dx*( d1 + 2*d2 + 2*d3 + d4 )/6
def integrationEDO(start, end, step, v_ini, derivee):
'''
Application de la méthode de RK4 pour chaque temps
'''
n = v_ini.size # nombre d'equa-diff du 1er ordre
# Création du tableau temps
interval = end - start # Intervalle
num_points = int(interval / step) + 1 # Nombre d'éléments
t = np.linspace(start, end, num_points) # Tableau temps t
print('Pas de temps = ',t[1]-t[0],' s') # affiche intervalle de temps
# Initialisation du tableau v
v = np.empty((n, num_points))
# tableau y qui contient les valeurs actuelles de v au temps t[i]
y = np.empty(n)
# Condition initiale
v[:, 0] = v_ini
y[:] = v_ini
distVm=np.empty(num_points)
# Boucle for
for i in range(num_points - 1):
if (y[0]-xT)**2 + (y[1]-yT)**2 < Rt**2:
print("Le vaisseau atterit sur la Terre") #vérification que le vaisseau n'a pas crashé sur Terre
break
if (y[0]-xM)**2 + (y[1]-yM)**2 < Rm**2:
print("Le vaisseau atterit sur Mars") #vérification que le vaisseau n'a pas crashé sur Mars
break
y = rk4(t[i], step, y, derivee) # modifie le tableau y
v[:, i + 1] = y
distVm[i]= ((y[0]-xM)**2 + (y[1]-yM)**2)**(1/2)
# Argument de sortie
return t[:i], v[:,:i], i, distVm[:i]
def derivee_u (t, u):
global xM
global yM
global xT
global yT
global phi
#position de la Terre
xT= dST * np.cos(2*np.pi*t/(365.25*3600*24))
yT= dST * np.sin(2*np.pi*t/(365.25*3600*24))
#positions de Mars
xM= dSM * np.cos(2*np.pi*t/(686.7*3600*24) + np.radians(phi))
yM= dSM * np.sin(2*np.pi*t/(686.7*3600*24) + np.radians(phi))
# Initialisation de la dérivée
du = np.empty(u.size)
# Déclaration des équations
dxt= u[0]-xT
dyt= u[1]-yT
dxm= u[0]-xM
dym= u[1]-yM
du[0] = u[2]
du[1]= u[3]
du[2] = -(G *Mt * dxt)/((dxt**2 + dyt**2)**(3/2)) - (G *Mm * dxm)/((dxm**2 + dym**2)**(3/2))
du[3]= -(G *Mt * dyt)/((dxt**2 + dyt**2)**(3/2)) - (G *Mm * dym)/((dxm**2 + dym**2)**(3/2))
return du
def decoupe(secondes):
""" convertit les secondes en jours/heures/minutes/secondes"""
minutes = secondes // 60
secondes %= 60
heures = minutes // 60
minutes %= 60
jours = heures // 24
heures %= 24
return(f"{jours}jour(s), {heures}heure(s), {minutes}minute(s), {secondes}seconde(s)")
def fonction(tmax, dt, v0, angle, phi_local):
""" affiche la trajectoire du vaisseau (les positions x et y) en fonction du temps ,
trouve les valeurs du temps et des coordonnées à l'aide des condition initales imposés et de la fonction integrationEDO"""
global phi
phi=phi_local # pour que la variable phi soit utilisable dans les autres fonctions
# Conditions initiales
u_ini = np.array([dST, Rt, v0*np.cos(np.radians(angle)) - dST * (2*np.pi/(365.25*3600*24))*np.sin(2*np.pi/(365.25*3600*24)) , v0*np.sin(np.radians(angle))+ dST * (2*np.pi/(365.25*3600*24))*np.cos(2*np.pi/(365.25*3600*24))])
# Méthode d'Euler
t, u, i, distVm = integrationEDO(0, tmax, dt, u_ini, derivee_u)
d = min(list(distVm)) # valeur minimum de de la distance vaisseau-Mars
D = list(distVm) # Liste des distances vaisseau-Mars
#position de la Terre
xT= dST * np.cos(2*np.pi*t/(365.25*3600*24))
yT= dST * np.sin(2*np.pi*t/(365.25*3600*24))
#positions de Mars
xm= dSM * np.cos(2*np.pi*t/(686.7*3600*24) + np.radians(phi))
ym= dSM * np.sin(2*np.pi*t/(686.7*3600*24) + np.radians(phi))
#dimensions du cadre
plt.figure(figsize=(7,7))
#traçage trajectoire du vaisseau
plt.plot(u[0,:],u[1,:],'r',label='Trajectoire du vaisseau')
#traçage Terre
plt.plot(xT, yT,'k',label='Tracé de la Terre')
#traçage trajectoire de la Mars
plt.plot(xm,ym,'b',label='Trajectoire de Mars')
#légendes
plt.title('Trajectoire de retour libre')
plt.xlabel('Position en x')
plt.ylabel('Position en y')
plt.legend()
plt.grid()
plt.figure()
#traçage distance vaisseau Mars dans le temps
plt.plot(t,distVm)
plt.title('Distance Vaisseau Mars en fonction du temps')
plt.xlabel('Temps(s)')
plt.ylabel('Distance Vaisseau Mars (m)')
#plt.ylim(1500* 10**3,8000* 10**3)
#traçage force attraction mars sur vaisseau
#plt.plot(t, )
plt.grid()
plt.figure()
#traçage vitesse vaisseau dans le temps
plt.plot(t, (u[2,:]**2+u[3,:]**2)**(1/2) ,'b')
plt.title('Vitesse du vaisseau en fonction du temps')
plt.ylabel('Vitesse Vaisseau (m.s-1)')
plt.xlabel('Temps(s)')
plt.grid()
plt.figure()
print("distance minimale = ",(distVm.min()-Rm)*10**(-3), 'km', 'à', decoupe(D.index(d)*dt))
plt.show()
print('durée voyage=',decoupe(i*dt)) |
Partager