Extraction d'un maximum

[HSbF6]

Bonjour,

j'ai programmé un simulateur d'orbite qui permet d'afficher l'orbite d'un satellite autour d'un astre quelconque (constante de l'astre dans le programme). Ce simulateur prend en compte les frottements atmosphériques (uniquement dans le cas d'une trajectoire autour de la Terre). Le simulateur marche bien. Seulement, dans le cas d'une trajectoire qui rentre dans l'atmosphère, j'aimerais pouvoir afficher le maximum de décélération (noté amax ligne 100) ressenti par le satellite sphérique (masse et diamètre en constante). Or amax, faisant appel à un tableau numpy et un à une fonction auxiliaire, il m'est impossible d'utiliser la fonction max() pour extraire amax. Comment faire svp ?


Le code:

from scipy.integrate import odeint
import numpy as np
from matplotlib import pyplot as plt
from math import pi
from math import sqrt
from math import exp

t = np.linspace(0, 50000, 200000)  

G = 6.6743e-11
M = 5.9737e24
R = 6371
A = 120
Cx = 0.5
S = pi*4
masse = 4000


def density(altitude):
    if altitude <= 11000:
        rho = 1.225 * (1-(0.0065*altitude)/288.15)**4.5
    elif altitude <= 12000:
        rho = -(6e-5)*altitude+1
    elif altitude <= 14000:
        rho = -(4.5e-5)*altitude+0.81
    elif altitude <= 16000:
        rho = -(3.5e-5)*altitude+0.68
    elif altitude <= 20000:
        rho = -(1.65e-5)*altitude+0.38
    elif altitude <= 50000:
        rho = 0.14388 * 1.225 * exp(-altitude/16800)-0.004
    elif altitude <= 120000:  
        rho = 0.0001388 * 1.225 * exp(-altitude/30300)
    else:
        rho = 0
    return rho

def F(f, t):
    [r, rp, theta] = f
    altitude = r - R*1000
    rho = density(altitude) 
    thetap = C /(r**2)
    drag = -0.5 * Cx * rho * (thetap*r)**2 * S
    return [rp, r*thetap**2 - G*M/(r**2) + (drag/masse), thetap]

n = input("\n" + "Vitesse orthoradiale initiale en m/s: ") 
p = input("\n" + "Vitesse radiale initiale en m/s: ") 
m = input("\n" + "Altitude initiale par rapport à la surface en km: ")

r0 = float(m)*1000 + float(R)*1000
thetap0 = float(n) / r0
C = thetap0 * (r0 ** 2) 
f0 = [r0, p, 0]
f = odeint(F, f0, t)
r = f[:, 0]
theta = f[:, 2]
# =============================================================================
# 
# ============================================================================= 

Rm = R*1000
r[r<Rm] = None


thetap = C /(r**2)
vmax = round((max(thetap*r)), 2)
VMI®®®n = round((min(thetap*r)), 2)


x = (r * np.cos(theta + pi/2))/1000
y = (r * np.sin(theta + pi/2))/1000


circle1 = plt.Circle((0, 0), R+A , color='b')
circle2 = plt.Circle((0, 0), R , color='sienna')
fig, ax = plt.subplots()
ax.set_aspect(1)
ax.add_artist(circle1)
ax.add_artist(circle2)

plt.plot(x, y, "steelblue")
plt.axis("equal") 
plt.show()
    
    
apo = round(((max(r)/1000)-R), 2)
           
peri = round(((min(r)/1000)-R), 2)

d = ((apo+peri+2*R)/2)*1000

P = sqrt((4*(pi**2)*(d**3))/(G*M))

k = round((P/60), 2)




if np.any(r < (R+A)*1000):
    print("\n" + """Votre orbite pénètre l'atmosphère !""" + "\n" + "\n" + "Décélération maximale = " + str(amax) + " g")
else:
    print("\n" + "Apogée: " + "Altitude = " + str(apo) + " km, " + "Vitesse = " + str(VMI®®®n) + " m/s ")
               
    print("\n" + "Périgée: " + "Altitude = " + str(peri) + " km, " + "Vitesse = " + str(vmax) + " m/s ")

    print("\n" + "Période de révolution = " + str(k) + " min" + " (" + str(round((k/1440),2)) + " j" + ")" )

-----

[HSbF6]

PS: La décélération apparente a s'obtient en faisant a = drag/masse. Drag contient des tableaux (r et thetap) et une fonction annexe (density)

[umfred]

il est où le calcul de amax ? sinon il y a plusieurs fonctions numpy pour récupérer un maximum (np.amax(), np.maximum(), ...)

[HSbF6]

Justement il n’y a pas encore de calcul de amax, c’est ce que je demande. Jai juste préparer le amax a la fin. Et j’avais essayer de fzire une boucle for a part pour calculer amax, avec np.amax etc, mais ca parcours les 120000 valeurs donc c’était tres long et en plus ca me mettait des messages d’erreur. Il n’y a pas un moyen plus rapide d’extraire amax directement de la fonction F(f, t) ? Comme je l’ai fait pour vmin et vmax qui dépendent de r.

[A voir en vidéo sur Futura]

[umfred]

les calculs matriciels ça existent ... donc calculer les décélérations subies aux divers moments et ensuite en extraire le maximum

[HSbF6]

C'est à dire ? Des listes de listes ? Peut-être mais je pensais qu'on pouvait extraire amax aussi facilement que j'ai extrait vmin et vmax, vu que vmin et vmax dépendent aussi de la variable r dans F(f, t).

[umfred]

PS: La décélération apparente a s'obtient en faisant a = drag/masse. Drag contient des tableaux (r et thetap) et une fonction annexe (density)

Du coup c'est quoi la formule de calcul avec les éléments que tu as ?

[HSbF6]

Ca serait amax = max(drag/(masse*9,81)) pour avoir le résultat en équivalent g. Mais du coup j’ai mis max par ca marche pas comme fonction.

[umfred]

et sans le max ça donne quoi comme résultat ? un vecteur ? tu as essayé les fonction dont j'ai parlé plus haut (np.amax, ...)

[HSbF6]

Bah j’ai essayé les fonctions max pour les tableau numpy. Mais le problème c’est qu’il faut que je redéfinisse la variable rho (vu qu’elle est que dans la fonction density). Donc ca m’oblige à faire une boucle sur l’altitude qui allonge énormément le temps d’execution. Le faite de passer la variable rho en variable globale (avec global rho) ne change rien. Alors que pour vmin et vmax, le programme extrait directement le maximum du tableau, sans boucle rallongeant le programme.

[umfred]

on ne sait toujours pas ce que ta fonction drag.*
density recalcule rho selon l'altitude (si j'en crois sa définition)
pour trouver le max de quelque chose, il faut bien calculer toutes les valeurs du quelque chose.

*ah si je viens de le retrouver, il est calculé dans F(), il faut peut être le faire sortir de F aussi
A tester:

def F(f, t):
    [r, rp, theta,drag] = f
    altitude = r - R*1000
    rho = density(altitude) 
    thetap = C /(r**2)
    drag = -0.5 * Cx * rho * (thetap*r)**2 * S
    return [rp, r*thetap**2 - G*M/(r**2) + (drag/masse), thetap,drag]
...
f0 = [r0, p, 0,0]
f = odeint(F, f0, t)
r = f[:, 0]
theta = f[:, 2]
drag=f[:,3]
..
amax = max(drag/(masse*9,81))

[HSbF6]

Ah oui ce qui me manquait c’était de savoir comment mettre en forme drag ans l’equa diff (avec la ligne drag = f[:,3]. Je vais essayer

[HSbF6]

La forme marche mais ca donne amax nul ou enorme (de l’ordre de 10^303). Peut-être le fait que la trajectoire soit saccadée en intervalle de temps donne une décélération énorme a chaque intervalle.