Propagazione orbitale con Javascript
Nei giorni scorsi ho voluto provare a mettere assieme un esperimento in Javascript che potesse mettere a frutto la mia passione per la dinamica orbitale (si, lo so, non è una frase che si legge tutti i giorni).
L'animazione qui sopra si può zoomare con la rotella del mouse, ruotare trascinando il mouse, e si può variare il parametro in alto a destra per cambiare la velocità orbitale. NB: in questo articolo è spiegato il calcolo, non come realizzare questa animazione bellissimissima che sarà oggetto di post futuri.
Nella sua forma più grezza e semplificata, il comportamento di un corpo in orbita deriva dall'equilibrio tra due effetti: la velocità tangenziale che tende a lanciarlo lontano dal corpo centrale e l'attrazione gravitazionale che ne cambia la traiettoria. Quando questi due effetti raggiungono un equilibrio abbiamo un'orbita. Secondo le leggi di Newton un'orbita può essere ellittica/circolare o parabolica, ed è eterna in assenza di forze esterne.
Proviamo a simulare il tutto con una manciata di righe in javascript, e vediamo che succede. Proveremo a modellizzare un satellite in una orbita circolare con dei parametri paragonabili a quelli della Stazione Spaziale Internazionale (circa 400km sopra la superficie terrestre). La simulazione avverrà attraverso una propagazione numerica: ogni tot secondi ricalcoliamo la posizione e la velocità del satellite prendendo in considerazione il vettore gravitazionale corrispondente in quell'istante. Il risultato saranno una serie di punti che potrebbero poi essere disegnati per mostrare l'orbita.
Per prima cosa ci occorre un oggetto che rappresenti le tre dimensioni di una forza, una posizione o una velocità, ovverossia un vettore. Questo oggetto deve contenere le tre variabili di cui sopra e prevedere alcune operazioni standard: somme, calcolo del modulo, inversione, normalizzazione. Possiamo usare una classe preesistente (ce ne sono diverse in giro) o ce lo possiamo fare noi:
export default class Vector{
constructor(x = 0, y = 0, z = 0){
this.x = x;
this.y = y;
this.z = z;
}
/**
* Sums another vector to this. Returns a new Vector.
* @param {Vector} other Other vector.
* @returns {Vector} This.
*/
add(other){
return new Vector(
this.x + other.x,
this.y + other.y,
this.z + other.z
);
}
/**
* Returns a new Vector as difference between this and the
* other.
* @param {Vector} other Other vector.
* @returns {Vector} This.
*/
diff(other){
return new Vector(
this.x - other.x,
this.y - other.y,
this.z - other.z
);
}
/**
* Returns the length of this Vector.
*/
module(){
return Math.sqrt(this.x*this.x + this.y*this.y + this.z*this.z);
}
/**
* Normalize. Returns a new Vector with the same direction
* but unitary length.
*/
norm(){
const module = this.module();
return new Vector(
this.x / module,
this.y / module,
this.z / module
)
}
/**
* Returns a scaled version of this Vector.
* @param {number} value
*/
scale(value){
return new Vector(
this.x * value,
this.y * value,
this.z * value
)
}
/**
* Returns an inverted version of this Vector.
*/
minus(){
return new Vector(
-this.x,
-this.y,
-this.z
)
}
/**
* Returns a copy of this Vector.
*/
clone(){
return new Vector(
this.x,
this.y,
this.z
)
}
}A questo punto dobbiamo modellizzare dei corpi celesti, di cui ci interessano velocità, posizione, massa.
import Vector from './vector.js'
export default class Body{
constructor(name){
this.name = name;
this.mass = 0;
this.radius = 0;
this.position = new Vector();
this.velocity = new Vector();
}
}Nel file principale possiamo inizializzare un modello della Terra e della nostra astronave. Nella semplificazione che stiamo operando la massa dell'astronave è considerata ininfluente e quindi la possiamo lasciare a zero. In compenso ne dobbiamo impostare posizione (400 km sopra la superficie, grossomodo la quota della stazione spaziale internazionale) e velocità (tangente alla superficie, di circa 7600 m/s che sappiamo essere la velocità orbitale della stazione spaziale internazionale). Questi parametri ci dovrebbero garantire un'orbita più o meno circolare. Ulteriori dettagli qui per chi non ha paura della matematica.
var earth = new Body('Earth');
earth.mass = 5.972e24;
earth.position = new Vector(0,0,0);
var ship = new Body('Ship');
ship.mass = 0;
ship.position = earth.position.add(new Vector(-6371e3-400e3, 0, 0));
ship.velocity = new Vector(0, 0, -7650)A questo punto dobbiamo modellizzare le forze agenti sulla nostra astronave a causa dell'attrazione gravitazionale terrestre. Per generalizzare costruiamo una funzione che a partire dall'elenco degli attrattori e dalla posizione del corpo attratto costruisce un vettore rappresentante la forza complessiva. L'accelerazione di gravità è calcolata come un vettore diretto verso il centro del corpo attrattore con una intensità pari a:
(dove M è la massa dell'attrattore in Kg, r è la distanza tra il satellite e il centro dell'attrattore in metri, G è la costante di gravitazione universale). Nel caso di più attrattori questi effetti si sommano (come somma tra vettori).
// Universal gravitation constant
export const G = 6.67e-11;
/**
* Acceleration function of a body.
* @param {Body[]} attractors
* @param {Vector} position
* @return {Vector} Acceleration.
*/
export function acceleration(attractors, position){
return attractors
.map(attr => {
let distVector = position.diff(attr.position);
return distVector.norm().minus().scale(
G * attr.mass / Math.pow(distVector.module(), 2)
);
})
.reduce((prev, curr) => prev.add(curr), new Vector())
}Per effettuare la propagazione potremmo fare così:
- Partiamo da posizione e velocità iniziali all'istante 0.
- Dopo TOT secondi, moltiplichiamo la velocità iniziale per TOT (per calcolare la distanza percorsa nel TOT di tempo) e aggiungiamo il risultato alla posizione iniziale. Così abbiamo la posizione nell'istante 1. Inoltre ricalcoliamo la nuova velocità aggiungendo alla velocità iniziale la velocità causata dall'accelerazione gravitazionale moltiplicata anch'essa per il TOT tempo.
- E così via, per tutti i punti successivi della simulazione.
Il problema è che questo metodo di propagazione è estremamente grossolano (non tiene conto in alcun modo delle variazioni di velocità che avvengono continuamente anche all'interno dei nostri TOT di tempo), e può rapidamente causare errori nella traiettoria simulata. Per ridurre questi problemi esistono varie metodologie di integrazione numerica, tra le quali il più noto è il metodo Runge-Kutta. Lascio i dettagli alla spiegazione esaustiva su Wikipedia. Siccome stasera mi sento particolarmente esoso voglio implementare un integratore con metodo Runge-Kutta di quarto ordine.
/**
* Returns final (position, velocity) array after time dt has passed.
* @param {Body[]} attractors Attractors influencing this object.
* @param {Vector} initialPosition Initial position.
* @param {Vector} initialVelocity Initial velocity.
* @param {Function} accFunction Acceleration function (attractors, position, velocity, deltaTime).
* @param {Number} dt Time step (seconds).
* @returns {Vector[]} Status vector [position, velocity]
*/
export function rk4(attractors, initialPosition, initialVelocity, accFunction, dt){
let position1 = initialPosition.clone();
let velocity1 = initialVelocity.clone();
let acceleration1 = accFunction(attractors, position1, velocity1, 0);
let position2 = initialPosition.add(velocity1.scale(0.5*dt));
let velocity2 = initialVelocity.add(acceleration1.scale(0.5*dt));
let acceleration2 = accFunction(attractors, position2, velocity2, dt/2);
let position3 = initialPosition.add(velocity2.scale(0.5*dt));
let velocity3 = initialVelocity.add(acceleration2.scale(0.5*dt));
let acceleration3 = accFunction(attractors, position3, velocity3, dt/2);
let position4 = initialPosition.add(velocity3.scale(dt));
let velocity4 = initialVelocity.add(acceleration3.scale(dt));
let acceleration4 = accFunction(attractors, position4, velocity4, dt);
let finalPosition = position1.add(
(velocity1
.add(velocity2.scale(2))
.add(velocity3.scale(2))
.add(velocity4)
).scale(dt/6)
);
let finalVelocity = velocity1.add(
(acceleration1
.add(acceleration2.scale(2))
.add(acceleration3.scale(2))
.add(acceleration4)
).scale(dt/6)
);
return [finalPosition, finalVelocity];
}Per effettuare uno step di propagazione basterà dare in pasto a questo algoritmo posizione e velocità attuale del nostro satellite, la lista di tutti gli attrattori e il tempo su cui calcolare la propagazione, oltre che la funzione di accelerazione che abbiamo scritto in precedenza.
/**
* Propagates the state vector of timestep seconds.
* @param {Vector} actualPosition Position of the orbiter.
* @param {Vector} actualVelocity Velocity of the orbiter.
* @param {Body[]} attractors Attractors inlfuencing this orbiter.
* @param {number} timestep Time step (in seconds).
* @return {Vector[]} New state [position, velocity]
*/
function propagate(actualPosition, actualVelocity, attractors, timestep){
return rk4(attractors, actualPosition, actualVelocity, acceleration, timestep);
}Chiamando iterativamente la funzione di propagazione e salvando i vettori risultanti possiamo costruire una simulazione dell'evoluzione dell'orbita del nostro satellite. Mettendo assieme quanto visto nel mio precedente post su Three.js ho voluto provare a disegnare l'orbita risultante. Il risultato è rappresentato dall'animazione in cima alla pagina. I dettagli delle implementazioni in Three.js li scriverò in qualche post futuro. La simulazione disegnata sopra è composta da 10k step di 10o secondi l'uno. La velocità può essere simulata a partire dal valore necessario ad ottenere un'orbita circolare fino a +5000 m/s (ben oltre la velocità di fuga dall'attrazione terrestre).
Grazie per la lettura e alla prossima!