Il y a aussi quelques améliorations qu'on peut y apporter
* initialiser directement les listes avec r+1 zéros et les remplir dans la boucle sera plus efficace que de faire à chaque fois un append et réallouer toute la mémoire.
* précalculer G*M_T et R_T**2
* écrire plutôt x, y = x + vx*h, y + vy*h est susceptible de faire comprendre à l'interpréteur qu'il peut paralléliser ces instructions. Efficacité marginale ici, mais ça peut servir ailleurs, je précise au cas où
* remplacer le t = tf par un break évite une assignation, mais en fait c'est inutile si tu remplaces ton while par un for avec le bon nombre d'itérations, puisque t = i*h où i est l'itération.
Tu auras automatiquement t>=tf dès que i*h >= tf ie i >= tf/h, ie i >= ceil(tf/h). Donc la bonne borne est effectivement r = ceil(tf/h) (for i in range(r): ...)
* aussi, on remplit VX et VY mais on ne les renvoie pas et on n'en fait rien, ce qui est du gâchis de mémoire. V2X n'apparaît même pas dans la suite du code, il faut le dégager aussi
* remplacer les *10**... par une notation scientifique qui est faite pour ça et donne directement la mantisse et l'exposant à Python. L'avantage d'utiliser 10** est que la masse et le rayon restent des entiers ce qui peut être utile si tu fais du calcul exact (fractions) mais ici on se contente d'un float
Code : Tout sélectionner
import numpy as np
import math
G,M_T,R_T = 6.67e-11, 5.972e24, 6.4e6
K, RT2 = G*M_T, R_T**2
def sat_in_geocentric_ref(x0, y0, vx0, vy0, n, tf):
""" Le calcul est mené dans le cadre de la mécanique Newtonienne dans le référientiel géocentrique.
x0,y0 correspondent aux coordonnées de départ
vx0, vy0 correspondent aux cordonnées du vecteur vitesse du satellite
"""
h = tf / float(n)
r = int(math.ceil(tf/h))
T = [i*h for i in range(r)]
X = np.zeros(r+1) # remplaçable par [0]*(r+1) ici parce que 0 est une constante. Sinon gare au piège de la référence
Y = np.zeros(r+1)
#VX = np.zeros(r+1)
#VY = np.zeros(r+1)
#V2X = np.zeros(r+1)
x, y, vx, vy = x0, y0, vx0, vy0
hK = h*K
for i in range(r): # xrange si tu es encore sur Python 2
d2 = x**2+y**2
if d2 < RT2:
print("CRASH")
break # ou mieux : return X[:i],Y[:i],T[:i]
d = d2**.5
x_n, y_n = x/d, y/d # vecteur normalisé
vx, vy = vx - (hK*x_n)/d2, vy - (hK*y_n)/d2
x, y = x + h*vx, y + h*vy
X[i], Y[i] = x, y
#VX[i], VY[i] = vx, vy
return X, Y, T
X, Y, T = sat_in_geocentric_ref(36000e3 + R_T, 0, 0, 3e3, 86400, 86400)
for u in (T,X,Y): print(u)
* enfin dernière remarque, le break dans le for peut te faire retourner un tuple de listes ayant de nombreux zéros. Tu peux le remplacer par return X[:i],Y[:i],T[:i]