Calcolare e disegnare un'orbita a partire dal vettore di stato con Javascript e ThreeJS

In un precedente post ho raccontato come valutare l'evoluzione della posizione di un satellite in orbita propagandone lo stato (posizione/velocità) tenendo in considerazione l'effetto di attrazione gravitazionale del corpo centrale.

Nella nostra simulazione ci troviamo in un caso ideale: un corpo orbitante di massa trascurabile rispetto all'attrattore, e che risente solo dell'influenza dell'unico corpo attrattore, senza anomalie varie. Questo problema è noto in fisica come "Problema dei due corpi" ed è risolvibile in forma analitica. Conoscendo lo stato del corpo (vettori velocità e posizione) e la massa del corpo centrale possiamo ricavare tutti i parametri necessari per disegnare la sua orbita, senza bisogno di iterare nel tempo.

In questa prima parte descriverò i calcoli necessari a ricavare i parametri orbitali, con qualche cenno teorico. In un post successivo mostrerò come rappresentare l'orbita in Javascript e ThreeJS.

Le regole del gioco

Un'orbita è sempre un'ellisse, con il corpo centrale in uno dei due fuochi. A meno che non superiamo la velocità di fuga, nel qual caso diventa un'iperbole. Un'orbita può essere circolare, poiché il cerchio è un caso particolare di ellisse in cui i due fuochi coincidono.
Un'orbita è per sempre: finché il satellite non subisce altre forze, continuerà in eterno a seguire la stessa traiettoria chiusa.
Un'orbita ha due punti particolarmente interessanti: il periapside/pericentro (o perigeo se siamo attorno alla Terra), che è il punto di minor distanza tra il satellite e il corpo attrattore, e l'apoapside/apocentro (o apogeo), che è il punto in cui la distanza è maggiore. Se l'orbita è circolare sono indefiniti (non esistono punti più distanti o meno distanti).
La velocità del corpo orbitante varia di continuo, e raggiunge un minimo all'apoapside e un massimo al periapside (più ci si avvicina all'attrattore più la velocità aumenta).
Un'orbita è una figura geometrica bidimensionale, quindi giace su un piano.

Come si descrive un'orbita

Esistono diversi modi tutti diversi e tutti validi per definire un'orbita (la geometria è una brutta bestia), ma i più tradizionali sono i sei valori definiti "parametri kepleriani" (in realtà per disegnare la traiettoria ce ne servono solo 5, poichè il sesto, l'anomalia vera, rappresenta la posizione del satellite all'interno dell'orbita).

I primi due parametri servono a definirne la forma e dimensione e sono:

Semiasse maggiore: la metà della lunghezza dell'ellisse (distanza tra gli apsidi). Non è indicata nel disegno sopra.
Eccentricità: indicazione di quanto l'ellisse è "allungata", ha un valore compreso tra zero (orbita perfettamente circolare) a uno non incluso (ellisse molto allungata). Non indicata nel disegno sopra.

Gli altri tre descrivono la posizione dell'orbita nello spazio tridimensionale:

Inclinazione: inclinazione del piano orbitale rispetto al piano di riferimento (per la Terra il riferimento è tipicamente l'eclittica, ovvero il piano su cui giace l'orbita che il pianeta segue attorno al Sole).
Longitudine del nodo ascendente. L'orbita attraversa in due punti distinti il piano di riferimento: questi punti sono detti "nodi". In uno di questi punti il satellite si sta muovendo verso l'"alto" ("alto" secondo una convenzione, nel caso della Terra è la direzione del polo nord): questo è il nodo ascendente. La longitudine del nodo ascendente è la distanza angolare tra questo nodo e un riferimento celeste, il punto vernale. Nota: non ha nulla a che vedere con la longitudine in quanto coordinata geografica!
Argomento del pericentro: l'angolo, lungo il piano orbitale, tra il periapside e il nodo ascendente.

Per finire il sesto, che descrive la posizione corrente del satellite all'interno della propria orbita: l'anomalia vera, ovvero l'angolo in cui il satellite si trova rispetto al proprio periapside (se il satellite è nel periapside, l'anomalia vera è zero; all'apoapside è 180 gradi). Gli angoli si considerano originare dal fuoco dell'ellisse, ovvero dal centro del corpo centrale.

Esiste un altro tipo di anomalia, detta anomalia eccentrica, che descrive un parametro equivalente ma considerando i vettori con origine il centro geometrico dell'ellisse. E, volendo esiste ancora un altro tipo di anomalia, l'anomalia media, che descrive l'angolo del satellite rispetto al periapside (centrando i vettori sul fuoco, come l'anomalia vera) ma fingendo che il satellite abbia una velocità angolare costante. Fingiamo quindi che, in un tempo pari a un decimo del periodo orbitale, il satellite percorra un decimo dell'angolo giro, quindi 36 gradi. Naturalmente non è così, poiché la velocità angolare varia sempre nel corso dell'orbita e quindi l'angolo vero sarà diverso, ma questo parametro è molto facile da calcolare semplicemente valutando il tempo percorso dall'ultimo passaggio dal periapside ed è quindi il parametro che viene riportato normalmente nelle descrizioni ufficiali delle orbite dei satelliti, come le TLE.

E' un pò contorto, ma lo schema di prima dovrebbe aiutare. Ve l'ho detto, la geometria è una brutta bestia.

Un pò di calcoli: disegnare l'ellisse

Il nostro punto di partenza è questo:

Abbiamo un satellite di cui conosciamo posizione (espressa con il vettore r0 che rappresenta la congiungente tra il centro del corpo attrattore e il satellite stesso), velocità v0 e angolo phi (tra r0 e v0).

Per prima cosa dobbiamo calcolare l'energia del nostro satellite. L'energia è pari alla somma dell'energia cinetica e dell'energia potenziale (quest'ultima, ricordiamo, con segno negativo). Naturalmente in questo contesto usiamo i moduli dei vettori v e r. Quindi:

dove m è la massa del satellite, v la velocità, G la costante di gravitazione universale e M la massa del corpo attrattore. Siccome sappiamo che le caratteristiche di un'orbita non sono influenzate dalla massa del satellite (che abbiamo definito trascurabile), usiamo l'energia specifica, che è pari all'energia diviso la massa.

Si noti che l'energia di un'orbita non varia con il passare del tempo ed è sempre negativa (se è positiva abbiamo un'orbita aperta, quindi una iperbole).

Dall'energia specifica possiamo calcolare il semiasse maggiore:

Il passo successivo è ricordarci la formula del momento angolare (altro parametro costante per un corpo in orbita), che è uguale a:

A questo punto sfoderiamo una delle formule per il calcolo dell'eccentricità a partire dai parametri di cui sopra, e facciamo un pò di semplificazioni.

A questo punto abbiamo calcolato eccentricità e semiasse maggiore, e abbiamo tutti gli strumenti per disegnare la nostra orbita usando l'equazione polare dell'ellisse a partire da un fuoco: scegliamo una serie di angoli tra 0 e 360 e per ciascuno possiamo calcolare la distanza r. Uniamo i punti e abbiamo la nostra ellisse!

Altri numeri interessanti (anche se non ci servono)

A partire dai parametri dell'ellisse possiamo calcolare il raggio di apoapside e periapside (in questo contesto il raggio è sempre la distanza tra l'oggetto orbitante e il centro del corpo attrattore; se vogliamo calcolare l'altezza del satellite rispetto al terreno dobbiamo sottrarre il raggio medio terrestre, ovvero 6'371 km).

Inoltre, a partire dall'equazione "vis-viva" (legge dell'invarianza dell'energia orbitale) possiamo calcolare la velocità che il nostro satellite assumerà nei due punti:

Per finire, possiamo anche ricavare il periodo orbitale, ovvero il tempo che il satellite impiega per percorrere una singola orbita:

Ancora un pò di calcoli: ruotare l'ellisse

Una volta disegnata l'ellisse dobbiamo ricavare gli altri tre parametri orbitali per poterla correttamente ruotare nello spazio tridimensionale. Finora abbiamo usato i moduli della distanza e della velocità iniziali, ora dobbiamo iniziare a fare i conti con i vettori.

Per prima cosa dobbiamo calcolare tre vettori fondamentali, che ci serviranno per calcolare il resto dei parametri. Il primo è il vettore del momento angolare h, calcolato come prodotto vettoriale tra la distanza e la velocità. Per definizione, è un vettore perpendicolare al piano della nostra orbita.

Il secondo è il vettore n, il vettore dei nodi, calcolato come il prodotto vettoriale del vettore unitario Y perpendicolare al piano equatoriale (quindi verticale) e del momento angolare. Per essere perpendicolare a Y, n deve giacere nel piano equatoriale; per essere perpendicolare ad h, n deve giacere nel piano della nostra orbita. Di conseguenza, n non può che giacere all'intersezione dell'orbita con l'eclittica, dove si trovano i nodi.

Per concludere, calcoliamo il vettore eccentricità, che è un vettore che parte dal centro del corpo attrattore e punta verso il perigeo.

Poichè l'inclinazione è l'angolo tra K e h, la si può ricavare come:

hy è la componente del vettore h lungo l'asse verticale, e h il modulo del vettore. Il prossimo passo è il calcolo della longitudine del nodo ascendente. E' l'angolo tra l'asse x (il riferimento vernale, nel nostro sistema di riferimento) e il vettore dei nodi n, quindi si può calcolare così:

Successivamente calcoliamo l'argomento del perigeo, che è l'angolo tra i vettori e ed n.

E con questo, se il cielo vuole, abbiamo finito :-). Ora dobbiamo solo ruotare la nostra ellisse:

Attorno all'asse verticale usando la longitudine del nodo ascendente.
Attorno al vettore h, usando l'argomento di perigeo.
Attorno al vettore e, usando l'inclinazione.

Come informazione aggiuntiva possiamo calcolare anche l'anomalia vera con questa formula:

E con questo la teoria è finita. Appena mi riprendo dalle serate passando a scrivere formule in LaTeX e sfogliare un libro faticoso (ma molto interessante e consigliato a chi vuole approfondire), provo a mettere giù la seconda parte in cui proveremo a mettere a frutto quando faticosamente imparato qui sopra per provare a fare in Javascript qualcosa del genere:

Nel frattempo buonanotte, alla prossima!