1 Orbite dei Pianeti
In questo tutorial vediamo come calcolare l'orbita di un pianeta attorno al sole. Le equazioni da cui partire sono quelle che descrivono le componenti della forza di attrazione gravitazionale a cui è soggetto il pianeta in esame da parte del sole, considerato fisso in uno dei due fuochi dell'orbita:
Anche in questo caso abbiamo a che fare con le equazioni differenziali ordinarie, come nel caso dell'oscillatore. La novità è che trattandosi di un sistema a due dimensioni, ne dovremo risolvere due contemporaneamente, una per componente.
Vediamo di seguito il corpo principale del programma, constituito da un ciclo for che contiene uno dei metodi di integrazione visti in precedenza e lo esegue iterativamente:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | for istep in range(Nstep):
Fx=-GM*x/r**3
Fy=-GM*y/r**3
#Eulero-Cromer
#vx+=tau*Fx
#vy+=tau*Fy
#x+=tau*vx
#y+=tau*vy
#r=sqrt(x**2+y**2)
#velocity-Verlet
x+=vx*tau+0.5*Fx*tau**2
y+=vy*tau+0.5*Fy*tau**2
r=sqrt(x**2+y**2)
F1x=-GM*x/r**3
F1y=-GM*y/r**3
vx+=0.5*tau*(Fx+F1x)
vy+=0.5*tau*(Fy+F1y)
x_arr[istep+1]=x
y_arr[istep+1]=y
vx_arr[istep+1]=vx
vy_arr[istep+1]=vy
|
- abbiamo scritto anche qui due metodi, quello di Cromer e quello di Verlet, in modo che commentando uno dei due possiamo provarne l'accuratezza e l'affidabilità;
- poichè il sistema è bidimensionale abbiamo due componenti per la posizione , la velocità , e l'accelerazione , ognuna delle quali si calcola allo stesso modo di quello previsto dallo schema di integrazione;
- calcoliamo il raggio a partire dalle coordinate appena calcolate, perchè esso è richiesto per il calcolo dell'accelerazione;
- al posto della forma var1=var0+..., con due variabili, denominate con 0 e 1, che indicavano la variabile allo tempo e , usiamo la forma var+=... che aggiunge al valore gia contenuto nella variabile il nuovo valore; questo ci permette di non dover aggiornare la variabile 0 con quella 1 dopo averla calcolata (var0=var1);
- anche qui salviamo la posizione e la velocità in appositi arrays.
Affinchè il programma sia completo ci servono le definizioni delle costanti usate dentro il ciclo e delle condizioni iniziali, vediamole di seguito:
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 | #Costanti moto
GM = 4*pi**2 # G*M in unita astronomiche (UA)
m=5.9722E+24 # massa terra (kg)
#Posizioni e velocità iniziali in UA, al perielio
x=0.98
y=0.0
r=sqrt(x**2+y**2)
vx=0.0
vy=2*pi
# Nstep*tau è il tempo totale in anni
Nstep=1000 # Numero di passi temporali
tau=0.001 # lunghezza singolo passo temporale
#arrays posizione e velocità
x_arr=zeros(Nstep+1)
y_arr=zeros(Nstep+1)
vx_arr=zeros(Nstep+1)
vy_arr=zeros(Nstep+1)
#Salvataggio delle condizioni iniziali
x_arr[0]=x
y_arr[0]=y
vx_arr[0]=vx
vy_arr[0]=vy
|
Notiamo come si siano usate le unità astronomiche: quindi le lunghezze espresse in termini di distanza media Terra-Sole e il tempo in anni. Queste unità sono più adatte al sistema che stiamo trattando e ci permettono di ridurre il numero di iterazioni da compiere, quindi il tempo di attesa e il costo computazionale.
Le condizioni iniziali proposte sono relative al pianeta Terra. Ci aspettiamo quindi di osservare una singola orbita (Nstep*tau=1 anni) quasi circolare. Facciamo il grafico al solito modo, poi eguagliamo la scala dei due assi in maniera che l'orbita non venga distorta:
1 2 | plot(x_arr,y_arr)
axis('equal')
|
Per generare le orbite di altri pianeti basterà impostare la posizione e la velocità iniziali come fatto per la Terra, regolare il numero di step, eseguire nuovamente il programma e infine fare il grafico. Potremo per comodità chiedere all'utente di inserire questi parametri, come abbiamo gia visto per l'oscillatore:
1 2 3 | x=input('Posizione iniziale perielio')
vy=input('Velocità iniziale perielio')
Nstep=input('Numero di passi')
|
Eseguendo il programma più volte, lasciando la prima finestra del grafico aperta, inseriamo i seguenti paremetri: Nstep=1000, , ; quindi lasciamo fisso il numero di passi e la posizione iniziale e cambiamo la velocità iniziale lungo y. Questo è il grafico che si ottiene:
Come osserviamo, le due orbite più grandi della blu (terrestre) non si chiudono per via del fatto che Nstep=1000 non è sufficiente. Viceversa le orbite più piccole della blu sono piu lunghe di un singolo giro.
Se volessimo calcolare altre quantità legate al moto del pianeta in orbita attorno al sole, come ad esempio il raggio, l'angolo tra raggio e asse x, velocità angolare, energia cinetica e potenziale, il metodo è sempre lo stesso: aggiungiamo un array per la grandezza scelta, impostiamo la prima componente in base alle condizioni iniziali e poi scriviamo ad ogni iterazione il nuovo valore nelle componenti successive, come di seguito mostrato per l'angolo :
1 2 3 4 5 6 7 8 9 10 | ...
theta_arr=zeros(Nstep+1)
theta_arr[0]=arctan2(y,x)
...
for istep in range(Nstep):
...
theta_arr[istep+1]=arctan2(y,x)
time_arr=linspace(0,Nstep*tau,Nstep+1)
plot(time_arr,theta_arr)
|
La funzione arctan2(y,x) ci permette di calcolare l'angolo compreso tra raggio e asse x in base alle coordinate della posizione del pianeta. Questo è il grafico che si ottinene per il caso della terra e per il caso di maggiore eccentricità dell'orbita (ottenuto con):
L'angolo cresce linearmente nel caso di orbita circolare, come nel caso approsimato della terra. I salti che si vedono son dovuti al fatto che la funzione artan2() restituisce solo angoli compresi tra . Quando l'orbita è ellitica, la velocità del pianeta non è piu costante per cui l'angolo varia piu velocemente in corrispondenza del afelio.