Calcolare e disegnare un'orbita a partire dal vettore di stato con Javascript e ThreeJS - seconda parte
Come promesso nel precedente post, voglio descrivere sommariamente come si può calcolare un'orbita a partire da un vettore di stato e disegnarla con ThreeJS.
La spiegazione non sarà particolarmente discorsiva poichè è già stato detto tutto nel precedente post, ma ho pensato che potesse essere utile condividere il codice javascript che ho usato per mettere in pratica quanto detto in precedenza. Il codice è semplice e (decentemente) commentato.
La prima cosa da fare è definire una classe Vector che descriva un vettore, e che servirà a contenere tutte le informazioni geometriche. Inoltre in questa classe saranno implementate tutte le operazioni che riguardano i vettori: somme, scalature, prodotto scalare e prodotto vettoriale, modulo, normalizzazione. La classe avrà questo aspetto:
import * as THREE from '../build/three.module.js'
export default class Vector{
constructor(x = 0, y = 0, z = 0){
this.x = x;
this.y = y;
this.z = z;
}
static fromPolar(radius, inclination, argument){
return new Vector(
radius * Math.sin(inclination) * Math.cos(argument),
radius * Math.sin(inclination) * Math.cos(inclination),
radius * Math.cos(inclination)
)
}
/**
* 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
)
}
print(){
return({
x: this.x.toExponential(2),
y: this.y.toExponential(2),
z: this.z.toExponential(2),
module: this.module().toExponential(2)
})
}
/**
* Returns a copy of this Vector.
*/
clone(){
return new Vector(
this.x,
this.y,
this.z
)
}
/**
* Dot product of this Vector vs another.
* @param {Vector} other Other Vector
* @returns Dot product
*/
dot(other){
return (this.x*other.x + this.y*other.y + this.z*other.z)
}
/**
* Cross product of this Vector vs another.
* @param {Vector} other Other Vector
* @returns Cross product
*/
cross(other){
let x = this.y * other.z - this.z * other.y
let y = this.z * other.x - this.x * other.z
let z = this.x * other.y - this.y * other.x
return new Vector(x, y, z)
}
/**
* Angle between this and another Vector
* @param {Vector} other Other Vector.
* @returns Angle (rad)
*/
angle(other){
return Math.acos(this.dot(other) / (this.module() * other.module()))
}
/**
* Returns this vector as THREE.Vector3
* @returns
*/
toTHREEVector3(){
return new THREE.Vector3(
this.x,
this.y,
this.z
)
}
}Il passo successivo è costruire una classe destinata a contenere le informazioni relative a un singolo corpo celeste, in modo che siano raggruppate in maniera logica: positione, velocità, massa... Questa classe contiene anche i riferimenti alle mesh threejs che useremo per il disegno.
import Vector from './vector.js'
export default class Body{
constructor(name, mass, radius, SOIRadius){
// Nome del corpo celeste
this.name = name;
// Massa [Kg]
this.mass = mass;
// Raggio medio [m]
this.radius = radius;
// Raggio sfera influenza [m]
this.SOIRadius = SOIRadius;
// Posizione [m, m, m]
this.position = new Vector();
// Velocità [m^2, m^2, m^2]
this.velocity = new Vector();
this.mesh = null;
this.speedMesh = null;
this.lineMesh = null;
}
/**
* Calcola la posizione di un satellite.
* @param {*} height Altezza sopra la superficie [m]
* @param {*} inclination Inclinazione [°]
* @param {*} azimuth Azimuth [°]
* @returns
*/
calcSatellitePosition(height, inclination, azimuth){
let r = height + this.radius;
let theta = 2 * Math.PI * inclination / 360;
let phi = 2 * Math.PI * azimuth / 360;
let x = r * Math.cos(phi) * Math.sin(theta);
let y = r * Math.sin(phi) * Math.sin(theta);
let z = r * Math.cos(theta);
return this.position.add(new Vector(x, y, z));
}
clone(){
let newBody = new Body(this.name);
newBody.name = this.name;
newBody.mass = this.mass;
newBody.radius = this.radius;
newBody.SOIRadius = this.SOIRadius;
newBody.position = this.position.clone();
newBody.velocity = this.velocity.clone();
newBody.mesh = null;
return newBody;
}
}A questo punto dobbiamo creare un metodo che effettui tutti i calcoli. Non entro nel dettaglio dell'implementazione poichè segue passo passo le formule descritte nel precedente post.
/**
* Calculater orbital parameters of orbiting body.
* @param {Body} orbitingBody
* @param {Body} centreBody
* @returns
*/
export function orbitalCalcBody(orbitingBody, centreBody){
let velocity = orbitingBody.velocity
let position = orbitingBody.position
let attractorPosition = centreBody.position
// Constants
let M = centreBody.mass
// Initial state vector
let relPosition = position.diff(attractorPosition)
let v_0 = velocity.module()
let r_0 = relPosition.module()
let angle_0 = relPosition.angle(velocity)
// Specific energy (Earth orbit)
let specificEnergy = Math.pow(v_0, 2) / 2.0 - G * M / r_0
// Semimajor axis
let semiMajorAxis = - G * M / (2 * specificEnergy)
// Orbit eccentricity
let eccentricity = Math.sqrt(1+(
(2 * Math.pow(v_0, 2) * Math.pow(r_0, 2) * Math.pow(Math.sin(angle_0), 2) * specificEnergy)
/(Math.pow(G, 2) * Math.pow(M, 2))
))
// Orbit shape
let rApoapsis = (semiMajorAxis * (1 + eccentricity))
let rPeriapsis = (semiMajorAxis * (1 - eccentricity))
// Orbital period
let period = 2 * Math.PI * Math.sqrt(Math.pow(semiMajorAxis, 3) / (G * M))
// Specific angular moment vector
let h = relPosition.cross(velocity)
// Inclination
let inclination = Math.acos(h.y / h.module())
// Nodes vector
let n = new Vector(0, 1, 0).cross(h)
// Longitude of ascending node
let longAscNode = Math.acos(n.x / n.module())
if (n.z < 0) longAscNode = 2*Math.PI - longAscNode
if (isNaN(longAscNode)) longAscNode = 0
// Eccentricity vector
let eccVector = velocity
.cross(h)
.scale(1/(G*M))
.diff(relPosition.scale(1/relPosition.module()))
// Argument of periapsis
let argPeriapsis = Math.acos(n.dot(eccVector)/(n.module() * eccVector.module()))
if (eccVector.y < 0) argPeriapsis = 2*Math.PI - argPeriapsis
if (isNaN(argPeriapsis)) argPeriapsis = 0
// Velocità
let vApoapsis = Math.sqrt(G * M * ((2/rApoapsis)-(1/semiMajorAxis)))
let vPeriapsis = Math.sqrt(G * M * ((2/rPeriapsis)-(1/semiMajorAxis)))
// True anomaly
let trueAnomaly = Math.acos(eccVector.dot(relPosition) / (eccVector.module() * relPosition.module()))
if (relPosition.dot(velocity) < 0){
trueAnomaly = 2 * Math.PI - trueAnomaly
}
return({
v_0,
r_0,
angle_0,
specificEnergy,
eccentricity,
semiMajorAxis,
rApoapsis,
rPeriapsis,
period,
inclination,
argPeriapsis,
longAscNode,
vApoapsis,
vPeriapsis,
eccVector,
trueAnomaly,
orbitingBody,
centreBody
})
}Finalmente possiamo iniziare a disegnare qualcosa. Partiamo da un progetto ThreeJS (come ad esempio quanto descritto in questo post). Per prima cosa definiamo un metodo setMeshPosition, che ci servirà a posizionare le nostre mesh al punto giusto.
// Fattore di scala globale (la Terra misura 10 unità)
const scaleFactor = 10 / 6371000;
/**
* Set the position of the mesh.
* @param {Body} body Corpo da posizionare.
*/
export function setMeshPosition(body){
body.mesh.position.x = body.position.x * scaleFactor;
body.mesh.position.y = body.position.y * scaleFactor;
body.mesh.position.z = body.position.z * scaleFactor;
}Definiamo anche un metodo che ci servirà per inizializzare la mesh che useremo per disegnare l'orbita:
/**
* Inizializza la mesh lineare per disegnare la propagazione di una orbita.
* @param {*} scene
* @param {*} simSize Numero di step.
* @param {*} color Colore.
* @param {boolean} dashed
* @returns
*/
export function buildLineMesh(simSize, color = 'red', dashed = false){
// Create mesh with fake points, because the BufferGeometry has to be
// initialized with the right size.
const newSimPoints = []
for (let i = 0; i < simSize; i++){
newSimPoints.push(new THREE.Vector3(0,0,0));
}
const simGeometry = new THREE.BufferGeometry().setFromPoints(newSimPoints);
let simMaterial;
if(!dashed){
simMaterial = new THREE.LineBasicMaterial({
color : color
});
}
else {
simMaterial = new THREE.LineDashedMaterial({
linewidth: 1,
color: color,
dashSize: .5,
gapSize: .1
})
}
const mesh = new THREE.Line( simGeometry, simMaterial );
mesh.visible = false;
return mesh;
}Già che ci siamo, definiamo il metodo che riceverà i risultati dei nostri calcoli e li userà per piazzare la mesh lungo il percorso dell'orbita.
L'orbita viene disegnata a partire dal fuoco (la posizione del corpo attraente) dividendo l'angolo giro in ANGLESTEPS angoli e calcolando per ciascuno il raggio, usando la formula dell'ellisse. Non è il metodo ideale (l'angolo zero corrisponde al periapside, quindi le posizioni vicine a 180 gradi saranno più lontani tra loro e quindi l'ellisse risulterà meno fluida), ma per ora ci accontentiamo.
La mesh viene poi ruotata tre volte usando i valori di inclinazione, aegomento di periapside e longitudine del nodo ascendente.
export function orbitDraw(calcOrbit, orbitSim, angleSteps){
let eccentricity = calcOrbit.eccentricity
let semiMajorAxis = calcOrbit.semiMajorAxis
let inclination = calcOrbit.inclination
let argPeriapsis = calcOrbit.argPeriapsis
let longAscNode = calcOrbit.longAscNode
let eccVector = calcOrbit.eccVector
// Simulate orbit from parameters
for (let angleStep = 0; angleStep < angleSteps; angleStep++){
let theta = (angleStep * 360 / angleSteps) * Math.PI / 180;
// Ellisse
let r_theta = semiMajorAxis * (1 - Math.pow(eccentricity,2)) / (1 + eccentricity * Math.cos(theta))
let x = r_theta * Math.cos(theta) * scaleFactor;
let y = 0;
let z = r_theta * Math.sin(theta) * scaleFactor;
let posArray = orbitSim.geometry.getAttribute('position').array;
posArray[angleStep*3] = x;
posArray[angleStep*3+1] = y;
posArray[angleStep*3+2] = z;
}
orbitSim.geometry.setDrawRange(0, angleSteps)
orbitSim.geometry.attributes.position.needsUpdate = true;
orbitSim.visible = true;
orbitSim.computeLineDistances();
// Rotazione su parametri orbitali
// Reset rotation
orbitSim.rotation.x = 0
orbitSim.rotation.y = 0
orbitSim.rotation.z = 0
// Apply longitude of ascending node
orbitSim.rotateOnAxis(new THREE.Vector3(0,1,0).normalize(), -longAscNode)
// Apply argument of periapsis
let eccVectorPerp = calcOrbit.orbitingBody.position.diff(calcOrbit.centreBody.position).cross(calcOrbit.orbitingBody.velocity).norm()
orbitSim.rotateOnWorldAxis(eccVectorPerp.toTHREEVector3(), argPeriapsis)
// Apply inclination
orbitSim.rotateOnWorldAxis(eccVector.norm().toTHREEVector3(), inclination)
}Descriviamo la Terra:
var earth = new Body('Earth', 5.972e24, 6371e3, 0.929e9);
// Costruzione mesh
let earthGeometry = new THREE.SphereGeometry(
earth.radius * scaleFactor,
50,
50
);
let earthMaterial = new THREE.MeshPhongMaterial({
map: new THREE.ImageUtils.loadTexture("./images/2k_earth_daymap.jpg"),
color: 0xaaaaaa,
specular: 0x333333,
shininess: 5
});
earth.mesh = new THREE.Mesh(earthGeometry, earthMaterial);
scene.add(earth.mesh);
setMeshPosition(earth);..e il nostro satellite:
var ship = new Body('Ship', 100000, 200000, 0);
// Create ship mesh
var shipGeometry = new THREE.SphereGeometry(200000 * scaleFactor, 50, 50 );
var shipMaterial = new THREE.MeshPhongMaterial({
color: 'lightgreen',
transparent: true,
opacity: .8
});
ship.mesh = new THREE.Mesh(shipGeometry, shipMaterial);
scene.add(ship.mesh);
setMeshPosition(ship);
// Impostiamo posizione e velocità iniziali
ship.position: earth.calcSatellitePosition(400e3, 80, -40),
ship.velocity: new Vector(0, 1000, -7670-2500) Inizializziamo la mesh che useremo per disegnare l'orbita del satellite. Il numero di step definisce quanti punti saranno usati per disegnare la traiettoria.
let angleSteps = 36*2;
let orbitSim = buildLineMesh(angleSteps, 'yellow', true)
scene.add(orbitSim)A questo punto abbiamo tutti i tasselli, dobbiamo solo richiamare i metodi di calcolo e di disegno:
// Calcola e disegna orbita simulata
let calcOrbit = orbitalCalcBody(ship, earth)
orbitDraw(calcOrbit, orbitSim, angleSteps, scaleFactor)Voilà!
Alla prossima!