gitea-mirror/KSP-MGA-Planner
1
0
Fork 0
mirror of https://github.com/Krafpy/KSP-MGA-Planner.git synced 2026-01-09 03:42:04 -08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5 from Krafpy/new-physics3d

Browse Source 
New physics3d
This commit is contained in:
Krafpy
2021-12-24 17:40:17 +01:00
committed by GitHub
parent 0210eb9390 01a733fbd1
commit bbd0792a6d
31 changed files with 2560 additions and 2257 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
data/config.yml
View File

@@ -41,11 +41,12 @@ trajectorySearch:
splitLimit: 1000 # maximum input chunk size per worker in the worker pool, exceeded if all workers are already used
crossoverProba: 0.9 # crossover probability (CR) of the DE algorithm
diffWeight: 0.8 # differential weight (F) of the DE algorithm
depDVScaleMin: 1.05 # The minimum ejection velocity, in terms of scale of the minimum velocity required to escape the body
depDVScaleMin: 1.01 # The minimum ejection velocity, in terms of scale of the minimum velocity required to escape the body
depDVScaleMax: 3 # The maximum ejection velocity
dsmOffsetMin: 0.01 # The minimum offset of a DSM on an interplanetary leg
dsmOffsetMax: 0.99 # The maximum offset of a DSM
minLegDuration: 21600 # The minimum duration of a leg (s)
fbRadiusMaxScale: 4
popSizeDimScale: 500 # The population size is equal to this value times the dimension of the search space (number of compnents agent vector)
maxGenerations: 200 # Maximum number of evolution iterations

56
dist/dedicated-workers/libs/common.js vendored
View File

@@ -1,4 +1,7 @@
"use strict";
function postMessageSafe(msg) {
postMessage(msg);
}
function sendProgress(progress, data) {
postMessageSafe({ label: "progress", progress, data });
}
@@ -8,38 +11,35 @@ function debug(...data) {
function sendResult(result) {
postMessageSafe({ label: "complete", result: result });
}
function postMessageSafe(msg) {
postMessage(msg);
}
class WorkerEnvironment {
static init(Env) {
const env = new Env();
onmessage = ({ data }) => {
switch (data.label) {
case "initialize":
env.onWorkerInitialize(data.config);
postMessageSafe({ label: "initialized" });
break;
case "run":
env.onWorkerRun(data.input);
break;
case "continue":
env.onWorkerContinue();
break;
case "stop":
env.onWorkerStop();
postMessageSafe({ label: "stopped" });
break;
case "pass":
env.onWorkerDataPass(data.data);
postMessageSafe({ label: "received" });
break;
}
};
}
onWorkerInitialize(data) { }
onWorkerRun(input) { }
onWorkerContinue(input) { }
onWorkerStop() { }
onWorkerDataPass(data) { }
}
function initWorker(Env) {
const env = new Env();
onmessage = ({ data }) => {
switch (data.label) {
case "initialize":
env.onWorkerInitialize(data.config);
postMessageSafe({ label: "initialized" });
break;
case "run":
env.onWorkerRun(data.input);
break;
case "continue":
env.onWorkerContinue();
break;
case "stop":
env.onWorkerStop();
postMessageSafe({ label: "stopped" });
break;
case "pass":
env.onWorkerDataPass(data.data);
postMessageSafe({ label: "received" });
break;
}
};
}

138
dist/dedicated-workers/libs/evolution.js vendored
View File

@@ -1,64 +1,94 @@
"use strict";
function populationChunkFitness(popChunk, fitness) {
const fitChunk = Array(popChunk.length).fill(0);
for (let i = 0; i < popChunk.length; i++) {
fitChunk[i] = fitness(popChunk[i]);
class ChunkedEvolver {
constructor(chunkStart, chunkEnd, agentDim, fitness, cr, f) {
this.chunkStart = chunkStart;
this.chunkEnd = chunkEnd;
this.agentDim = agentDim;
this.fitness = fitness;
this.cr = cr;
this.f = f;
const ce = chunkEnd;
const cs = chunkStart;
this.chunkSize = ce - cs + 1;
}
return fitChunk;
}
function evolvePopulationChunk(population, fitnesses, chunkStart, chunkEnd, cr, f, fitness) {
const chunkSize = chunkEnd - chunkStart + 1;
const dim = population[0].length;
const nextPopChunk = Array(chunkSize).fill(null);
const nextFitChunk = Array(chunkSize).fill(0);
for (let j = chunkStart; j <= chunkEnd; j++) {
const x = population[j];
const fx = fitnesses[j];
const [a, b, c] = pick3(population, j);
const ri = randint(0, dim - 1);
const y = Array(dim).fill(0);
for (let i = 0; i < dim; i++) {
if (Math.random() < cr || i == ri) {
y[i] = a[i] + f * (b[i] - c[i]);
createRandomAgent() {
const agent = Array(this.agentDim);
this.randomizeAgent(agent);
return agent;
}
randomizeAgent(agent) {
for (let i = 0; i < this.agentDim; i++) {
agent[i] = Math.random();
}
}
createRandomPopulationChunk() {
const popChunk = [];
for (let i = 0; i < this.chunkSize; i++) {
popChunk.push(this.createRandomAgent());
}
return popChunk;
}
evaluateChunkFitness(popChunk) {
const fitChunk = Array(popChunk.length);
for (let i = 0; i < popChunk.length; i++) {
fitChunk[i] = this.fitness(popChunk[i]);
}
return fitChunk;
}
evolvePopulationChunk(population, fitnesses) {
const dim = population[0].length;
const nextPopChunk = Array(this.chunkSize);
const nextFitChunk = Array(this.chunkSize);
for (let j = this.chunkStart; j <= this.chunkEnd; j++) {
const x = population[j];
const fx = fitnesses[j];
const [a, b, c] = this._pick3(population, j);
const ri = randint(0, dim - 1);
const y = Array(dim).fill(0);
for (let i = 0; i < dim; i++) {
if (Math.random() < this.cr || i == ri) {
y[i] = a[i] + this.f * (b[i] - c[i]);
}
else {
y[i] = x[i];
}
y[i] = clamp(y[i], 0, 1);
}
const fy = this.fitness(y);
const index = j - this.chunkStart;
if (fy < fx) {
nextPopChunk[index] = y;
nextFitChunk[index] = fy;
}
else {
y[i] = x[i];
nextPopChunk[index] = x;
nextFitChunk[index] = fx;
}
}
const fy = fitness(y);
const index = j - chunkStart;
if (fy < fx) {
nextPopChunk[index] = y;
nextFitChunk[index] = fy;
return {
popChunk: nextPopChunk,
fitChunk: nextFitChunk
};
}
_pick3(population, parentIndex) {
const swap = (arr, i, j) => {
const t = arr[j];
arr[j] = arr[i];
arr[i] = t;
};
swap(population, parentIndex, 0);
const picked = [null, null, null];
const pickedIndices = [0, 0, 0];
for (let i = 0; i <= 2; i++) {
const ri = randint(1 + i, population.length - 1);
picked[i] = population[ri];
pickedIndices[i] = ri;
swap(population, ri, i);
}
else {
nextPopChunk[index] = x;
nextFitChunk[index] = fx;
for (let i = 2; i >= 0; i--) {
swap(population, pickedIndices[i], i);
}
swap(population, parentIndex, 0);
return picked;
}
return {
popChunk: nextPopChunk,
fitChunk: nextFitChunk
};
}
function pick3(population, parentIndex) {
swap(population, parentIndex, 0);
const picked = [null, null, null];
const pickedIndices = [0, 0, 0];
for (let i = 0; i <= 2; i++) {
const ri = randint(1 + i, population.length - 1);
picked[i] = population[ri];
pickedIndices[i] = ri;
swap(population, ri, i);
}
for (let i = 2; i >= 0; i--) {
swap(population, pickedIndices[i], i);
}
swap(population, parentIndex, 0);
return picked;
}
function swap(arr, i, j) {
const t = arr[j];
arr[j] = arr[i];
arr[i] = t;
}

294
dist/dedicated-workers/libs/lambert.js vendored
View File

@@ -1,151 +1,155 @@
"use strict";
function solveLambert(r1vec, r2vec, tof, attractor) {
const mu = attractor.stdGravParam;
const r1 = mag3(r1vec);
const r2 = mag3(r2vec);
const c = mag3(sub3(r2vec, r1vec));
const s = 0.5 * (r1 + r2 + c);
const ir1 = div3(r1vec, r1);
const ir2 = div3(r2vec, r2);
const ih = normalize3(cross(ir1, ir2));
const lambda2 = 1 - c / s;
let lambda = Math.sqrt(lambda2);
let it1, it2;
if (ih.y < 0) {
lambda = -lambda;
it1 = cross(ir1, ih);
it2 = cross(ir2, ih);
}
else {
it1 = cross(ih, ir1);
it2 = cross(ih, ir2);
}
it1 = normalize3(it1);
it2 = normalize3(it2);
const lambda3 = lambda * lambda2;
const T = Math.sqrt(2 * mu / (s * s * s)) * tof;
const T0 = Math.acos(lambda) + lambda * Math.sqrt(1 - lambda2);
const T1 = 2 / 3 * (1 - lambda3);
let x0;
if (T >= T0) {
x0 = -(T - T0) / (T - T0 + 4);
}
else if (T <= T1) {
x0 = T1 * (T1 - T) / (2 / 5 * (1 - lambda2 * lambda3) * T) + 1;
}
else {
x0 = Math.pow(T / T0, 0.69314718055994529 / Math.log(T1 / T0)) - 1;
}
const x = householderIterations(T, x0, 1e-15, lambda, 15);
const gamma = Math.sqrt(mu * s / 2.0);
const rho = (r1 - r2) / c;
const sigma = Math.sqrt(1 - rho * rho);
const y = Math.sqrt(1.0 - lambda2 + lambda2 * x * x);
const vr1 = gamma * ((lambda * y - x) - rho * (lambda * y + x)) / r1;
const vr2 = -gamma * ((lambda * y - x) + rho * (lambda * y + x)) / r2;
const vt = gamma * sigma * (y + lambda * x);
const vt1 = vt / r1;
const vt2 = vt / r2;
const v1 = add3(mult3(ir1, vr1), mult3(it1, vt1));
const v2 = add3(mult3(ir2, vr2), mult3(it2, vt2));
return { v1, v2 };
}
function householderIterations(T, x0, eps, lambda, maxIters) {
let err = 1;
let xnew = 0;
let tof = 0;
let delta = 0;
let DT = 0, DDT = 0, DDDT = 0;
for (let it = 0; err > eps && it < maxIters; it++) {
tof = x2tof(x0, lambda);
const DTs = dTdx(DT, DDT, DDDT, x0, tof, lambda);
DT = DTs.DT;
DDT = DTs.DDT;
DDDT = DTs.DDDT;
delta = tof - T;
const DT2 = DT * DT;
xnew = x0 - delta * (DT2 - delta * DDT / 2) / (DT * (DT2 - delta * DDT) + DDDT * delta * delta / 6);
err = Math.abs(x0 - xnew);
x0 = xnew;
}
return x0;
}
function dTdx(DT, DDT, DDDT, x, T, lambda) {
const l2 = lambda * lambda;
const l3 = l2 * lambda;
const umx2 = 1.0 - x * x;
const y = Math.sqrt(1.0 - l2 * umx2);
const y2 = y * y;
const y3 = y2 * y;
DT = 1.0 / umx2 * (3.0 * T * x - 2.0 + 2.0 * l3 * x / y);
DDT = 1.0 / umx2 * (3.0 * T + 5.0 * x * DT + 2.0 * (1.0 - l2) * l3 / y3);
DDDT = 1.0 / umx2 * (7.0 * x * DDT + 8.0 * DT - 6.0 * (1.0 - l2) * l2 * l3 * x / y3 / y2);
return { DT, DDT, DDDT };
}
function x2tof2(x, lambda) {
const a = 1.0 / (1.0 - x * x);
if (a > 0) {
let alfa = 2.0 * Math.acos(x);
let beta = 2.0 * Math.asin(Math.sqrt(lambda * lambda / a));
if (lambda < 0.0)
beta = -beta;
return ((a * Math.sqrt(a) * ((alfa - Math.sin(alfa)) - (beta - Math.sin(beta)))) / 2.0);
}
else {
let alfa = 2.0 * Math.acosh(x);
let beta = 2.0 * Math.asinh(Math.sqrt(-lambda * lambda / a));
if (lambda < 0.0)
beta = -beta;
return (-a * Math.sqrt(-a) * ((beta - Math.sinh(beta)) - (alfa - Math.sinh(alfa))) / 2.0);
}
}
function x2tof(x, lambda) {
const battin = 0.01;
const lagrange = 0.2;
const dist = Math.abs(x - 1);
if (dist < lagrange && dist > battin) {
return x2tof2(x, lambda);
}
const K = lambda * lambda;
const E = x * x - 1.0;
const rho = Math.abs(E);
const z = Math.sqrt(1 + K * E);
if (dist < battin) {
const eta = z - lambda * x;
const S1 = 0.5 * (1.0 - lambda - x * eta);
let Q = hypergeometricF(S1, 1e-11);
Q = 4.0 / 3.0 * Q;
return (eta * eta * eta * Q + 4 * lambda * eta) / 2.0;
}
else {
const y = Math.sqrt(rho);
const g = x * z - lambda * E;
let d = 0.0;
if (E < 0) {
const l = Math.acos(g);
d = l;
var Lambert;
(function (Lambert) {
function solve(r1vec, r2vec, tof, attractor) {
const mu = attractor.stdGravParam;
const r1 = mag3(r1vec);
const r2 = mag3(r2vec);
const c = mag3(sub3(r2vec, r1vec));
const s = 0.5 * (r1 + r2 + c);
const ir1 = div3(r1vec, r1);
const ir2 = div3(r2vec, r2);
const ih = normalize3(cross(ir1, ir2));
const lambda2 = 1 - c / s;
let lambda = Math.sqrt(lambda2);
let it1, it2;
if (ih.y < 0) {
lambda = -lambda;
it1 = cross(ir1, ih);
it2 = cross(ir2, ih);
}
else {
const f = y * (z - lambda * x);
d = Math.log(f + g);
it1 = cross(ih, ir1);
it2 = cross(ih, ir2);
}
return (x - lambda * z - d / y) / E;
it1 = normalize3(it1);
it2 = normalize3(it2);
const lambda3 = lambda * lambda2;
const T = Math.sqrt(2 * mu / (s * s * s)) * tof;
const T0 = Math.acos(lambda) + lambda * Math.sqrt(1 - lambda2);
const T1 = 2 / 3 * (1 - lambda3);
let x0;
if (T >= T0) {
x0 = -(T - T0) / (T - T0 + 4);
}
else if (T <= T1) {
x0 = T1 * (T1 - T) / (2 / 5 * (1 - lambda2 * lambda3) * T) + 1;
}
else {
x0 = Math.pow(T / T0, 0.69314718055994529 / Math.log(T1 / T0)) - 1;
}
const x = householderIterations(T, x0, 1e-15, lambda, 15);
const gamma = Math.sqrt(mu * s / 2.0);
const rho = (r1 - r2) / c;
const sigma = Math.sqrt(1 - rho * rho);
const y = Math.sqrt(1.0 - lambda2 + lambda2 * x * x);
const vr1 = gamma * ((lambda * y - x) - rho * (lambda * y + x)) / r1;
const vr2 = -gamma * ((lambda * y - x) + rho * (lambda * y + x)) / r2;
const vt = gamma * sigma * (y + lambda * x);
const vt1 = vt / r1;
const vt2 = vt / r2;
const v1 = add3(mult3(ir1, vr1), mult3(it1, vt1));
const v2 = add3(mult3(ir2, vr2), mult3(it2, vt2));
return { v1, v2 };
}
}
function hypergeometricF(z, tol) {
let Sj = 1.0;
let Cj = 1.0;
let err = 1.0;
let Cj1 = 0.0;
let Sj1 = 0.0;
let j = 0;
while (err > tol) {
Cj1 = Cj * (3.0 + j) * (1.0 + j) / (2.5 + j) * z / (j + 1);
Sj1 = Sj + Cj1;
err = Math.abs(Cj1);
Sj = Sj1;
Cj = Cj1;
j++;
Lambert.solve = solve;
function householderIterations(T, x0, eps, lambda, maxIters) {
let err = 1;
let xnew = 0;
let tof = 0;
let delta = 0;
let DT = 0, DDT = 0, DDDT = 0;
for (let it = 0; err > eps && it < maxIters; it++) {
tof = x2tof(x0, lambda);
const DTs = dTdx(DT, DDT, DDDT, x0, tof, lambda);
DT = DTs.DT;
DDT = DTs.DDT;
DDDT = DTs.DDDT;
delta = tof - T;
const DT2 = DT * DT;
xnew = x0 - delta * (DT2 - delta * DDT / 2) / (DT * (DT2 - delta * DDT) + DDDT * delta * delta / 6);
err = Math.abs(x0 - xnew);
x0 = xnew;
}
return x0;
}
return Sj;
}
function dTdx(DT, DDT, DDDT, x, T, lambda) {
const l2 = lambda * lambda;
const l3 = l2 * lambda;
const umx2 = 1.0 - x * x;
const y = Math.sqrt(1.0 - l2 * umx2);
const y2 = y * y;
const y3 = y2 * y;
DT = 1.0 / umx2 * (3.0 * T * x - 2.0 + 2.0 * l3 * x / y);
DDT = 1.0 / umx2 * (3.0 * T + 5.0 * x * DT + 2.0 * (1.0 - l2) * l3 / y3);
DDDT = 1.0 / umx2 * (7.0 * x * DDT + 8.0 * DT - 6.0 * (1.0 - l2) * l2 * l3 * x / y3 / y2);
return { DT, DDT, DDDT };
}
function x2tof2(x, lambda) {
const a = 1.0 / (1.0 - x * x);
if (a > 0) {
let alfa = 2.0 * Math.acos(x);
let beta = 2.0 * Math.asin(Math.sqrt(lambda * lambda / a));
if (lambda < 0.0)
beta = -beta;
return ((a * Math.sqrt(a) * ((alfa - Math.sin(alfa)) - (beta - Math.sin(beta)))) / 2.0);
}
else {
let alfa = 2.0 * Math.acosh(x);
let beta = 2.0 * Math.asinh(Math.sqrt(-lambda * lambda / a));
if (lambda < 0.0)
beta = -beta;
return (-a * Math.sqrt(-a) * ((beta - Math.sinh(beta)) - (alfa - Math.sinh(alfa))) / 2.0);
}
}
function x2tof(x, lambda) {
const battin = 0.01;
const lagrange = 0.2;
const dist = Math.abs(x - 1);
if (dist < lagrange && dist > battin) {
return x2tof2(x, lambda);
}
const K = lambda * lambda;
const E = x * x - 1.0;
const rho = Math.abs(E);
const z = Math.sqrt(1 + K * E);
if (dist < battin) {
const eta = z - lambda * x;
const S1 = 0.5 * (1.0 - lambda - x * eta);
let Q = hypergeometricF(S1, 1e-11);
Q = 4.0 / 3.0 * Q;
return (eta * eta * eta * Q + 4 * lambda * eta) / 2.0;
}
else {
const y = Math.sqrt(rho);
const g = x * z - lambda * E;
let d = 0.0;
if (E < 0) {
const l = Math.acos(g);
d = l;
}
else {
const f = y * (z - lambda * x);
d = Math.log(f + g);
}
return (x - lambda * z - d / y) / E;
}
}
function hypergeometricF(z, tol) {
let Sj = 1.0;
let Cj = 1.0;
let err = 1.0;
let Cj1 = 0.0;
let Sj1 = 0.0;
let j = 0;
while (err > tol) {
Cj1 = Cj * (3.0 + j) * (1.0 + j) / (2.5 + j) * z / (j + 1);
Sj1 = Sj + Cj1;
err = Math.abs(Cj1);
Sj = Sj1;
Cj = Cj1;
j++;
}
return Sj;
}
})(Lambert || (Lambert = {}));

22
dist/dedicated-workers/libs/math.js vendored
View File

@@ -1,5 +1,6 @@
"use strict";
const TWO_PI = 2 * Math.PI;
const HALF_PI = 0.5 * Math.PI;
function clamp(x, min, max) {
return x > max ? max : x < min ? min : x;
}
@@ -17,27 +18,6 @@ function newtonRootSolve(f, df, x0, eps, maxIters = 1000) {
}
return x;
}
function randomPointOnSphereRing(thetaMin, thetaMax) {
return {
theta: randomInInterval(thetaMin, thetaMax),
phi: randomInInterval(0, TWO_PI)
};
}
function pointOnSphereRing(dir, theta, phi) {
let axis = normalize3(cross(dir, vec3(1, 0, 0)));
let point = rotate3(dir, axis, theta);
point = rotate3(point, dir, phi);
return point;
}
function minMaxRingAngle(r, rmin, rmax) {
return {
thetaMin: Math.asin(rmin / r),
thetaMax: Math.asin(rmax / r)
};
}
function randomInInterval(a, b) {
return a + Math.random() * (b - a);
}
function randint(a, b) {
return Math.floor(a + Math.random() * (b - a + 1));
}

8
dist/dedicated-workers/libs/physics-2d.js vendored
View File

@@ -18,11 +18,11 @@ var Physics2D;
const a = -mu * r0 / (r0 * v0 * v0 - 2 * mu);
const p = a * (1 - e * e);
const rp = a * (1 - e);
const pdir = mult2(div2(evec, e), rp);
const pvec = mult2(div2(evec, e), rp);
const clk = det(pos, vel) < 0;
return {
eccentricity: e,
periapsisDir: pdir,
periapsisVec: pvec,
semiMajorAxis: a,
orbitalParam: p,
clockwise: clk
@@ -65,7 +65,7 @@ var Physics2D;
const orbitElts = stateToOrbitElements(attractor, state);
const e = orbitElts.eccentricity;
const p = orbitElts.orbitalParam;
const pdir = orbitElts.periapsisDir;
const pvec = orbitElts.periapsisVec;
const pos0 = state.pos;
const vel0 = state.vel;
const tgRadius = target.orbit.semiMajorAxis;
@@ -87,7 +87,7 @@ var Physics2D;
x: -sinNu * vmag * vdir,
y: (e + cosNu) * vmag * vdir
};
const pangle = Math.atan2(pdir.y, pdir.x);
const pangle = Math.atan2(pvec.y, pvec.x);
return {
pos: rotate2(pos, pangle),
vel: rotate2(vel, pangle)

553
dist/dedicated-workers/libs/physics-3d.js vendored
View File

@@ -1,275 +1,324 @@
"use strict";
function orbitalElementsFromOrbitData(orbit) {
const right = vec3(1, 0, 0), up = vec3(0, 1, 0);
const ascNodeDir = rotate3(right, up, orbit.ascNodeLongitude);
const normal = rotate3(up, ascNodeDir, orbit.inclination);
const periapsisDir = rotate3(ascNodeDir, normal, orbit.argOfPeriapsis);
return {
semiMajorAxis: orbit.semiMajorAxis,
eccentricity: orbit.eccentricity,
periapsiDir: periapsisDir,
inclination: orbit.inclination,
argOfPeriapsis: orbit.argOfPeriapsis,
ascNodeLongitude: orbit.ascNodeLongitude,
ascNodeDir: ascNodeDir,
orbitalParam: orbit.orbitalParam
};
}
function equatorialCircularOrbit(radius) {
return {
eccentricity: 0,
periapsisDir: vec3(1, 0, 0),
semiMajorAxis: radius,
inclination: 0,
argOfPeriapsis: 0,
ascNodeLongitude: 0,
ascNodeDir: vec3(1, 0, 0),
orbitalParam: radius
};
}
function getBodyStateAtDate(body, date, mu) {
const { orbit } = body;
const n = orbit.meanMotion;
const M = body.meanAnomaly0 + n * date;
const nu = solveTrueAnomalyFromMeanAnomaly(orbit.eccentricity, M);
const orbitElts = orbitalElementsFromOrbitData(orbit);
return stateFromOrbitElements(orbitElts, mu, nu);
}
function soiExitTrueAnomaly(orbit, soi) {
const p = orbit.orbitalParam;
const e = orbit.eccentricity;
return Math.acos((p / soi - 1) / e);
}
function periapsisRadius(orbit) {
return orbit.semiMajorAxis * (1 - orbit.eccentricity);
}
function tofBetweenAnomalies(orbit, attractor, nu1, nu2) {
const mu = attractor.stdGravParam;
const a = orbit.semiMajorAxis;
const e = orbit.eccentricity;
let tof = 0;
if (e < 1) {
let E1 = eccentricAnomalyFromTrueAnomaly(nu1, e);
let E2 = eccentricAnomalyFromTrueAnomaly(nu2, e);
const invN = Math.sqrt(a * a * a / mu);
if (E2 < E1) {
E2 += TWO_PI;
var Physics3D;
(function (Physics3D) {
function orbitElementsFromOrbitData(orbit) {
const right = vec3(1, 0, 0), up = vec3(0, 1, 0);
const ascNodeDir = rotate3(right, up, orbit.ascNodeLongitude);
const normal = rotate3(up, ascNodeDir, orbit.inclination);
const periapsisDir = rotate3(ascNodeDir, normal, orbit.argOfPeriapsis);
return {
semiMajorAxis: orbit.semiMajorAxis,
eccentricity: orbit.eccentricity,
periapsiDir: periapsisDir,
inclination: orbit.inclination,
argOfPeriapsis: orbit.argOfPeriapsis,
ascNodeLongitude: orbit.ascNodeLongitude,
ascNodeDir: ascNodeDir,
orbitalParam: orbit.orbitalParam
};
}
Physics3D.orbitElementsFromOrbitData = orbitElementsFromOrbitData;
function equatorialCircularOrbit(radius) {
return {
eccentricity: 0,
periapsisDir: vec3(1, 0, 0),
semiMajorAxis: radius,
inclination: 0,
argOfPeriapsis: 0,
ascNodeLongitude: 0,
ascNodeDir: vec3(1, 0, 0),
orbitalParam: radius
};
}
Physics3D.equatorialCircularOrbit = equatorialCircularOrbit;
function bodyStateAtDate(body, orbit, attractor, date) {
const n = body.orbit.meanMotion;
const M = body.meanAnomaly0 + n * date;
const nu = trueAnomalyFromMeanAnomaly(orbit.eccentricity, M);
return orbitElementsToState(orbit, attractor, nu);
}
Physics3D.bodyStateAtDate = bodyStateAtDate;
function circularVelocity(attractor, radius) {
return Math.sqrt(attractor.stdGravParam / radius);
}
Physics3D.circularVelocity = circularVelocity;
function ejectionVelocity(attractor, radius) {
return Math.SQRT2 * circularVelocity(attractor, radius);
}
Physics3D.ejectionVelocity = ejectionVelocity;
function periapsisRadius(orbit) {
return orbit.semiMajorAxis * (1 - orbit.eccentricity);
}
Physics3D.periapsisRadius = periapsisRadius;
function trueAnomalyAtRadius(orbit, radius) {
const p = orbit.orbitalParam;
const e = orbit.eccentricity;
const cosNu = (p / radius - 1) / e;
if (Math.abs(cosNu) > 1)
throw new Error("The orbit never reaches that radius.");
return Math.acos(cosNu);
}
Physics3D.trueAnomalyAtRadius = trueAnomalyAtRadius;
function velocityAtRadius(orbit, attractor, radius) {
const mu = attractor.stdGravParam;
const a = orbit.semiMajorAxis;
const v2 = mu * (2 / radius - 1 / a);
if (v2 < 0)
throw new Error("Invalid radius.");
return Math.sqrt(v2);
}
Physics3D.velocityAtRadius = velocityAtRadius;
function velocityToReachAltitude(body, r0, r) {
const mu = body.stdGravParam;
const a = 0.5 * (r0 + r);
const v2 = mu * (2 / r0 - 1 / a);
if (v2 < 0)
throw new Error("Invalid radius.");
const vej2 = 2 * mu / r0;
return Math.sqrt(Math.min(vej2, v2));
}
Physics3D.velocityToReachAltitude = velocityToReachAltitude;
function orbitPeriod(attractor, a) {
return TWO_PI * Math.sqrt(a * a * a / attractor.stdGravParam);
}
Physics3D.orbitPeriod = orbitPeriod;
function specificEnergy(attractor, r, v) {
return 0.5 * v * v - attractor.stdGravParam / r;
}
Physics3D.specificEnergy = specificEnergy;
function deduceVelocityAtRadius(attractor, r0, v0, r) {
const se = specificEnergy(attractor, r0, v0);
return Math.sqrt(2 * (se + attractor.stdGravParam / r));
}
Physics3D.deduceVelocityAtRadius = deduceVelocityAtRadius;
function stateToOrbitElements(state, attractor) {
const mu = attractor.stdGravParam;
const pos = state.pos;
const vel = state.vel;
const r = mag3(pos);
const v2 = magSq3(vel);
const nullEps = 1e-10;
const h = cross(pos, vel);
let evec = sub3(div3(cross(vel, h), mu), div3(pos, r));
let e = mag3(evec);
if (e <= nullEps) {
evec = vec3(0, 0, 0);
e = 0;
}
tof = invN * (E2 - E1 + e * (Math.sin(E1) - Math.sin(E2)));
}
else {
let H1 = eccentricAnomalyFromTrueAnomaly(nu1, e);
let H2 = eccentricAnomalyFromTrueAnomaly(nu2, e);
if (H2 < H1) {
const t = H2;
H2 = H1;
H1 = t;
let nvec = cross(vec3(0, 1, 0), h);
const n = mag3(nvec);
let i = Math.acos(h.y / mag3(h));
let inXZPlane = false;
if (Math.abs(i) < nullEps) {
i = 0;
inXZPlane = true;
}
const invN = Math.sqrt(-a * a * a / mu);
tof = -invN * (H2 - H1 + e * (Math.sinh(H1) - Math.sinh(H2)));
else if (Math.abs(i - Math.PI) < nullEps) {
i = Math.PI;
inXZPlane = true;
}
else if (Math.abs(i + Math.PI) < nullEps) {
i = -Math.PI;
inXZPlane = true;
}
const cosEps = 1e-5;
let cos_lan = nvec.x / n;
if (Math.abs(cos_lan) < 1 + cosEps) {
cos_lan = clamp(cos_lan, -1, 1);
}
const t_lan = Math.acos(cos_lan);
let lan = nvec.z <= 0 ? t_lan : TWO_PI - t_lan;
let cos_arg = dot3(nvec, evec) / (e * n);
if (Math.abs(cos_arg) < 1 + cosEps) {
cos_arg = clamp(cos_arg, -1, 1);
}
const t_arg = Math.acos(cos_arg);
let arg = evec.y >= 0 ? t_arg : TWO_PI - t_arg;
const a = 1 / (2 / r - v2 / mu);
const p = a * (1 - e * e);
if (inXZPlane && e != 0) {
const t_plong = Math.acos(dot3(vec3(1, 0, 0), evec) / e);
const plong = evec.z < 0 ? t_plong : TWO_PI - t_plong;
lan = 0;
nvec = vec3(1, 0, 0);
arg = plong;
}
else if (!inXZPlane && e == 0) {
arg = 0;
evec = clone3(nvec);
}
else if (inXZPlane && e == 0) {
lan = 0;
nvec = vec3(1, 0, 0);
arg = 0;
evec = vec3(1, 0, 0);
}
return {
eccentricity: e,
periapsisDir: normalize3(evec),
semiMajorAxis: a,
inclination: i,
argOfPeriapsis: arg,
ascNodeLongitude: lan,
ascNodeDir: normalize3(nvec),
orbitalParam: p
};
}
return tof;
}
function circularVelocity(attractor, radius) {
return Math.sqrt(attractor.stdGravParam / radius);
}
function ejectionVelocity(attractor, radius) {
return Math.SQRT2 * circularVelocity(attractor, radius);
}
function velocityAtRadius(orbit, attractor, radius) {
return Math.sqrt(attractor.stdGravParam * (2 / radius - 1 / orbit.semiMajorAxis));
}
function idealEjectionDirection(bodyPos, retrograde) {
const tangent = normalize3(vec3(bodyPos.z, 0, -bodyPos.x));
let nu = Math.acos(bodyPos.x / mag3(bodyPos));
nu = bodyPos.z >= 0 ? nu : TWO_PI - nu;
if (retrograde) {
tangent.x *= -1;
tangent.z *= -1;
nu = (nu + Math.PI) % TWO_PI;
}
return { tangent, nu };
}
function hyperbolicEjectionOffsetAngle(perigeeVel, perigeeRadius, attractor) {
const e = perigeeRadius * perigeeVel * perigeeVel / attractor.stdGravParam - 1;
return Math.PI - Math.acos(-1 / e);
}
function stateToOrbitElements(state, attractor) {
const mu = attractor.stdGravParam;
const pos = state.pos;
const vel = state.vel;
const r = mag3(pos);
const v2 = magSq3(vel);
const h = cross(pos, vel);
let evec = sub3(div3(cross(vel, h), mu), div3(pos, r));
let e = mag3(evec);
if (e <= 1e-15) {
evec = vec3(0, 0, 0);
e = 0;
}
let nvec = cross(vec3(0, 1, 0), h);
const n = mag3(nvec);
const i = Math.acos(h.y / mag3(h));
const t_lan = Math.acos(nvec.x / n);
let lan = nvec.z <= 0 ? t_lan : TWO_PI - t_lan;
const t_arg = Math.acos(dot3(nvec, evec) / (e * n));
let arg = evec.y >= 0 ? t_arg : TWO_PI - t_arg;
const a = 1 / (2 / r - v2 / mu);
const p = a * (1 - e * e);
if (i == 0 && e != 0) {
const t_plong = Math.acos(dot3(vec3(1, 0, 0), evec) / e);
const plong = evec.z < 0 ? t_plong : TWO_PI - t_plong;
lan = 0;
nvec = vec3(1, 0, 0);
arg = plong;
}
else if (i != 0 && e == 0) {
arg = 0;
evec = clone3(nvec);
}
else if (i == 0 && e == 0) {
lan = 0;
nvec = vec3(1, 0, 0);
arg = 0;
evec = vec3(1, 0, 0);
}
return {
eccentricity: e,
periapsisDir: normalize3(evec),
semiMajorAxis: a,
inclination: i,
argOfPeriapsis: arg,
ascNodeLongitude: lan,
ascNodeDir: normalize3(nvec),
orbitalParam: p
};
}
function trueAnomalyFromOrbitalState(orbitElts, state) {
const pos = state.pos;
const vel = state.vel;
const r = mag3(pos);
const e = orbitElts.eccentricity;
const i = orbitElts.inclination;
let nu;
if (e != 0) {
const evec = orbitElts.periapsisDir;
const t_nu = Math.acos(dot3(evec, pos) / r);
const d = dot3(pos, vel);
Physics3D.stateToOrbitElements = stateToOrbitElements;
function orbitElementsToState(orbit, attractor, nu) {
const mu = attractor.stdGravParam;
const a = orbit.semiMajorAxis;
const e = orbit.eccentricity;
const p = orbit.orbitalParam;
nu *= -1;
const r = p / (1 + e * Math.cos(nu));
let pos = vec3(r * Math.cos(nu), 0, r * Math.sin(nu));
let vel;
if (e < 1) {
nu = d >= 0 ? t_nu : TWO_PI - t_nu;
const v = Math.sqrt(mu * a) / r;
const E = eccentricAnomalyFromTrueAnomaly(nu, e);
vel = vec3(v * Math.sin(E), 0, -v * Math.sqrt(1 - e * e) * Math.cos(E));
}
else {
nu = d >= 0 ? t_nu : -t_nu;
const v = Math.sqrt(-mu * a) / r;
const H = eccentricAnomalyFromTrueAnomaly(nu, e);
vel = vec3(v * Math.sinh(H), 0, -v * Math.sqrt(e * e - 1) * Math.cosh(H));
}
const right = vec3(1, 0, 0), up = vec3(0, 1, 0);
const ascNodeDir = rotate3(right, up, orbit.ascNodeLongitude);
pos = rotate3(pos, up, orbit.ascNodeLongitude);
pos = rotate3(pos, up, orbit.argOfPeriapsis);
pos = rotate3(pos, ascNodeDir, orbit.inclination);
vel = rotate3(vel, up, orbit.ascNodeLongitude);
vel = rotate3(vel, up, orbit.argOfPeriapsis);
vel = rotate3(vel, ascNodeDir, orbit.inclination);
return { pos, vel };
}
else if (i != 0) {
const nvec = orbitElts.ascNodeDir;
const t_u = Math.acos(dot3(nvec, pos) / r);
let u;
if (e < 1) {
u = pos.y >= 0 ? t_u : TWO_PI - t_u;
Physics3D.orbitElementsToState = orbitElementsToState;
function trueAnomalyFromOrbitalState(orbit, state) {
const pos = state.pos;
const vel = state.vel;
const r = mag3(pos);
const e = orbit.eccentricity;
const i = orbit.inclination;
let nu;
if (e != 0) {
const evec = orbit.periapsisDir;
const t_nu = Math.acos(dot3(evec, pos) / r);
const d = dot3(pos, vel);
if (e < 1) {
nu = d >= 0 ? t_nu : TWO_PI - t_nu;
}
else {
nu = d >= 0 ? t_nu : -t_nu;
}
}
else if (i != 0) {
const nvec = orbit.ascNodeDir;
const t_u = Math.acos(dot3(nvec, pos) / r);
let u;
if (e < 1) {
u = pos.y >= 0 ? t_u : TWO_PI - t_u;
}
else {
u = pos.y >= 0 ? t_u : -t_u;
}
nu = u;
}
else {
u = pos.y >= 0 ? t_u : -t_u;
const t_l = Math.acos(pos.x / r);
let l;
if (e < 1) {
l = vel.x <= 0 ? t_l : TWO_PI - t_l;
}
else {
l = vel.x <= 0 ? t_l : -t_l;
}
nu = l;
}
nu = u;
return nu;
}
else {
const t_l = Math.acos(pos.x / r);
let l;
Physics3D.trueAnomalyFromOrbitalState = trueAnomalyFromOrbitalState;
function eccentricAnomalyFromTrueAnomaly(nu, e) {
if (e < 1) {
l = vel.x <= 0 ? t_l : TWO_PI - t_l;
return 2 * Math.atan(Math.tan(nu * 0.5) * Math.sqrt((1 - e) / (1 + e)));
}
else {
l = vel.x <= 0 ? t_l : -t_l;
return 2 * Math.atanh(Math.tan(nu * 0.5) * Math.sqrt((e - 1) / (e + 1)));
}
nu = l;
}
return nu;
}
function eccentricAnomalyFromTrueAnomaly(nu, e) {
if (e < 1) {
return 2 * Math.atan(Math.tan(nu * 0.5) * Math.sqrt((1 - e) / (1 + e)));
Physics3D.eccentricAnomalyFromTrueAnomaly = eccentricAnomalyFromTrueAnomaly;
function meanAnomalyFromEccentricAnomaly(EH, e) {
if (e < 1) {
return EH - e * Math.sin(EH);
}
else {
return e * Math.sinh(EH) - EH;
}
}
else {
return 2 * Math.atanh(Math.tan(nu * 0.5) * Math.sqrt((e - 1) / (e + 1)));
Physics3D.meanAnomalyFromEccentricAnomaly = meanAnomalyFromEccentricAnomaly;
function meanAnomalyFromTrueAnomaly(nu, e) {
if (e < 1) {
const E = 2 * Math.atan(Math.tan(nu * 0.5) * Math.sqrt((1 - e) / (1 + e)));
const M = E - e * Math.sin(E);
return M;
}
else {
const H = 2 * Math.atanh(Math.tan(nu * 0.5) * Math.sqrt((e - 1) / (e + 1)));
const M = e * Math.sinh(H) - H;
return M;
}
}
}
function meanAnomalyFromEccentricAnomaly(EH, e) {
if (e < 1) {
return EH - e * Math.sin(EH);
Physics3D.meanAnomalyFromTrueAnomaly = meanAnomalyFromTrueAnomaly;
function trueAnomalyFromMeanAnomaly(e, M) {
if (e < 1) {
const E = newtonRootSolve(x => x - e * Math.sin(x) - M, x => 1 - e * Math.cos(x), M, 1e-15);
return 2 * Math.atan(Math.sqrt((1 + e) / (1 - e)) * Math.tan(E * 0.5));
}
else {
const H = newtonRootSolve(x => e * Math.sinh(x) - x - M, x => e * Math.cosh(x) - 1, M, 1e-15);
return 2 * Math.atan(Math.sqrt((e + 1) / (e - 1)) * Math.tanh(H * 0.5));
}
}
else {
return e * Math.sinh(EH) - EH;
Physics3D.trueAnomalyFromMeanAnomaly = trueAnomalyFromMeanAnomaly;
function propagateStateFromTrueAnomaly(orbit, attractor, nu0, deltaTime) {
const M0 = meanAnomalyFromTrueAnomaly(nu0, orbit.eccentricity);
const a = Math.abs(orbit.semiMajorAxis);
const mu = attractor.stdGravParam;
const n = Math.sqrt(mu / (a * a * a));
const MNext = M0 + n * deltaTime;
const nuNext = trueAnomalyFromMeanAnomaly(orbit.eccentricity, MNext);
return orbitElementsToState(orbit, attractor, nuNext);
}
}
function meanAnomalyFromTrueAnomaly(nu, e) {
if (e < 1) {
const E = 2 * Math.atan(Math.tan(nu * 0.5) * Math.sqrt((1 - e) / (1 + e)));
const M = E - e * Math.sin(E);
return M;
Physics3D.propagateStateFromTrueAnomaly = propagateStateFromTrueAnomaly;
function tofBetweenAnomalies(orbit, attractor, nu1, nu2) {
const mu = attractor.stdGravParam;
const a = orbit.semiMajorAxis;
const e = orbit.eccentricity;
let tof = 0;
if (e < 1) {
let E1 = eccentricAnomalyFromTrueAnomaly(nu1, e);
let E2 = eccentricAnomalyFromTrueAnomaly(nu2, e);
const invN = Math.sqrt(a * a * a / mu);
if (E2 < E1) {
E2 += TWO_PI;
}
tof = invN * (E2 - E1 + e * (Math.sin(E1) - Math.sin(E2)));
}
else {
let H1 = eccentricAnomalyFromTrueAnomaly(nu1, e);
let H2 = eccentricAnomalyFromTrueAnomaly(nu2, e);
if (H2 < H1) {
const t = H2;
H2 = H1;
H1 = t;
}
const invN = Math.sqrt(-a * a * a / mu);
tof = -invN * (H2 - H1 + e * (Math.sinh(H1) - Math.sinh(H2)));
}
if (tof < 0)
throw new Error("Negative TOF calculated.");
return tof;
}
else {
const H = 2 * Math.atanh(Math.tan(nu * 0.5) * Math.sqrt((e - 1) / (e + 1)));
const M = e * Math.sinh(H) - H;
return M;
}
}
function solveTrueAnomalyFromMeanAnomaly(e, M) {
if (e < 1) {
const E = newtonRootSolve(x => x - e * Math.sin(x) - M, x => 1 - e * Math.cos(x), M, 1e-15);
return 2 * Math.atan(Math.sqrt((1 + e) / (1 - e)) * Math.tan(E * 0.5));
}
else {
const H = newtonRootSolve(x => e * Math.sinh(x) - x - M, x => e * Math.cosh(x) - 1, M, 1e-15);
return 2 * Math.atan(Math.sqrt((e + 1) / (e - 1)) * Math.tanh(H * 0.5));
}
}
function stateFromOrbitElements(orbit, mu, nu) {
const a = orbit.semiMajorAxis;
const e = orbit.eccentricity;
const p = orbit.orbitalParam;
nu *= -1;
const r = p / (1 + e * Math.cos(nu));
let pos = vec3(r * Math.cos(nu), 0, r * Math.sin(nu));
let vel;
if (e < 1) {
const v = Math.sqrt(mu * a) / r;
const E = eccentricAnomalyFromTrueAnomaly(nu, e);
vel = vec3(v * Math.sin(E), 0, -v * Math.sqrt(1 - e * e) * Math.cos(E));
}
else {
const v = Math.sqrt(-mu * a) / r;
const H = eccentricAnomalyFromTrueAnomaly(nu, e);
vel = vec3(v * Math.sinh(H), 0, -v * Math.sqrt(e * e - 1) * Math.cosh(H));
}
const right = vec3(1, 0, 0), up = vec3(0, 1, 0);
const ascNodeDir = rotate3(right, up, orbit.ascNodeLongitude);
pos = rotate3(pos, up, orbit.ascNodeLongitude);
pos = rotate3(pos, up, orbit.argOfPeriapsis);
pos = rotate3(pos, ascNodeDir, orbit.inclination);
vel = rotate3(vel, up, orbit.ascNodeLongitude);
vel = rotate3(vel, up, orbit.argOfPeriapsis);
vel = rotate3(vel, ascNodeDir, orbit.inclination);
return { pos, vel };
}
function propagateState(orbit, nu0, deltaTime, attractor) {
const M0 = meanAnomalyFromTrueAnomaly(nu0, orbit.eccentricity);
const a = Math.abs(orbit.semiMajorAxis);
const mu = attractor.stdGravParam;
const n = Math.sqrt(mu / (a * a * a));
const MNext = M0 + n * deltaTime;
const nuNext = solveTrueAnomalyFromMeanAnomaly(orbit.eccentricity, MNext);
return stateFromOrbitElements(orbit, mu, nuNext);
}
function calculateOrbitPeriod(a, mu) {
return TWO_PI * Math.sqrt(a * a * a / mu);
}
function getHohmannPeriod(body1, body2, attractor) {
const orbit1 = body1.orbit;
const orbit2 = body2.orbit;
const { stdGravParam } = attractor;
const a = (orbit1.semiMajorAxis + orbit2.semiMajorAxis) * 0.5;
const period = calculateOrbitPeriod(a, stdGravParam);
return period;
}
Physics3D.tofBetweenAnomalies = tofBetweenAnomalies;
})(Physics3D || (Physics3D = {}));

510
dist/dedicated-workers/libs/trajectory-calculator.js vendored
View File

@@ -6,299 +6,339 @@ class TrajectoryCalculator {
this.sequence = sequence;
this.steps = [];
this._legs = [];
this.totalDeltaV = 0;
this._missionTime = 0;
this._sequenceIndex = 0;
this._departureDate = 0;
this._departureDVScale = 0;
this._flybys = [];
this._secondArcsData = [];
const attractorId = this._departureBody.orbiting;
this._mainAttractor = this.system[attractorId];
}
setParameters(depAltitude, startDateMin, startDateMax, params) {
this.depAltitude = depAltitude;
this.startDateMin = startDateMin;
this.startDateMax = startDateMax;
this.params = params;
this._getDepartureInfo();
this._clampDepartureInfos();
this._getLegsInfo();
}
compute() {
this.noError = false;
this._missionTime = this._departureDate;
this._calculateParkingOrbit();
this._calculateEjectionOrbit();
this._missionTime += this._lastStep.duration;
for (let i = 1; i < this.sequence.length; i++) {
this._sequenceIndex = i;
const leg = this._legs[i - 1];
if (!this._calculateLegFirstArc()) {
return;
}
this._calculateLegSecondArc();
this._updateLegInfos(i - 1);
this._missionTime += leg.duration;
this._calculateFlybyTrajectory();
this._missionTime += this._lastStep.duration;
if (i == this.sequence.length - 1) {
this._calculateArrivalCircularization();
}
}
this.noError = true;
}
get _destinationBody() {
const id = this.sequence[this.sequence.length - 1];
return this.system[id];
}
get _departureBody() {
return this.system[this.sequence[0]];
}
get _attractorMu() {
return this._mainAttractor.stdGravParam;
get _destinationBody() {
const last = this.sequence.length - 1;
return this.system[this.sequence[last]];
}
get _lastStep() {
return this.steps[this.steps.length - 1];
}
get _lastStepDateOfEnd() {
const { dateOfStart, duration } = this._lastStep;
get _lastStepEndDate() {
const last = this._lastStep;
const { dateOfStart, duration } = last;
return dateOfStart + duration;
}
_getDepartureInfo() {
this._departureDate = this.params[0];
this._departureDVScale = this.params[1];
addPrecomputedOrbits(bodiesOrbits) {
this._bodiesOrbits = bodiesOrbits;
}
_clampDepartureInfos() {
this._departureDate = clamp(this._departureDate, this.startDateMin, this.startDateMax);
const { depDVScaleMin, depDVScaleMax } = this.config;
this._departureDVScale = clamp(this._departureDVScale, depDVScaleMin, depDVScaleMax);
this.params[0] = this._departureDate;
this.params[1] = this._departureDVScale;
setParameters(depAltitude, startDateMin, startDateMax, params) {
this._depAltitude = depAltitude;
this._startDateMin = startDateMin;
this._startDateMax = startDateMax;
this._params = params;
this._getDepartureSettings();
this._getLegsSettings();
this._getFlybySettings();
}
_getLegsInfo() {
for (let i = 2; i < this.params.length; i += 4) {
this._legs.push({
duration: this.params[i],
dsmOffset: this.params[i + 1],
theta: this.params[i + 2],
phi: this.params[i + 3]
});
reset() {
this._secondArcsData = [];
this._legs = [];
this._flybys = [];
this.steps = [];
}
_getDepartureSettings() {
this._departureInfos = {
dateParam: this._params[0],
phaseParam: this._params[1],
ejVelParam: this._params[2]
};
}
_getLegsSettings() {
const { dsmOffsetMin, dsmOffsetMax } = this.config;
for (let i = 1; i < this.sequence.length; i++) {
const idx = 3 + (i - 1) * 4;
const infos = {
exitedBodyId: this.sequence[i - 1],
targetBodyId: this.sequence[i],
durationParam: this._params[idx],
duration: 0,
dsmParam: lerp(dsmOffsetMin, dsmOffsetMax, this._params[idx + 1])
};
this._computeLegDuration(infos);
this._legs.push(infos);
}
}
_updateLegInfos(legIndex) {
const leg = this._legs[legIndex];
this.params[2 + legIndex * 4] = leg.duration;
this.params[2 + legIndex * 4 + 1] = leg.dsmOffset;
this.params[2 + legIndex * 4 + 2] = leg.theta;
this.params[2 + legIndex * 4 + 3] = leg.phi;
_getFlybySettings() {
for (let i = 1; i < this.sequence.length - 1; i++) {
const idx = 3 + (i - 1) * 4 + 2;
const infos = {
flybyBodyId: this.sequence[i],
normAngleParam: this._params[idx],
periRadiParam: this._params[idx + 1]
};
this._flybys.push(infos);
}
}
_calculateArrivalCircularization() {
const body = this._destinationBody;
const flybyOrbit = this._lastStep.orbitElts;
const periapsisState = stateFromOrbitElements(flybyOrbit, body.stdGravParam, 0);
const progradeDir = normalize3(periapsisState.vel);
const circularVel = circularVelocity(body, mag3(periapsisState.pos));
const deltaVMag = Math.abs(circularVel - mag3(periapsisState.vel));
const deltaV = mult3(progradeDir, -deltaVMag);
this.totalDeltaV += deltaVMag;
const circularState = {
pos: periapsisState.pos,
vel: mult3(progradeDir, circularVel)
};
const circularOrbit = stateToOrbitElements(circularState, body);
const maneuvre = {
deltaVToPrevStep: deltaV,
progradeDir: progradeDir,
manoeuvrePosition: periapsisState.pos,
context: {
type: "circularization",
compute() {
this._appendDepartureParkingOrbit();
this._computeDepartureEjectionOrbit();
const numOfLegs = this._legs.length;
for (let i = 0; i < numOfLegs - 1; i++) {
const leg = this._legs[i];
const flyby = this._flybys[i];
this._computeFirstLegArc(leg);
this._computeLegSecondArcSimple(leg);
this._computeFlyby(flyby);
}
const last = this._legs[numOfLegs - 1];
this._computeFirstLegArc(last);
this._computeLegSecondArcSimple(last);
}
get totalDeltaV() {
let total = 0;
for (const step of this.steps) {
if (step.maneuvre) {
const { maneuvre } = step;
const { deltaVToPrevStep } = maneuvre;
const dv = mag3(deltaVToPrevStep);
total += dv;
}
};
this.steps.push({
orbitElts: circularOrbit,
attractorId: body.id,
beginAngle: 0,
endAngle: TWO_PI,
duration: 0,
dateOfStart: this._lastStep.dateOfStart + this._lastStep.duration,
maneuvre: maneuvre
});
}
return total;
}
_calculateFlybyTrajectory() {
const body = this.system[this.sequence[this._sequenceIndex]];
const localIncomingState = {
pos: sub3(this._flybyIncomingState.pos, this._flybyBodyState.pos),
vel: sub3(this._flybyIncomingState.vel, this._flybyBodyState.vel)
_computeLegDuration(infos) {
const exitedBody = this.system[infos.exitedBodyId];
const { durationParam } = infos;
const attr = this._mainAttractor;
if (infos.exitedBodyId != infos.targetBodyId) {
const targetBody = this.system[infos.targetBodyId];
const exBodyA = exitedBody.orbit.semiMajorAxis;
const tgBodyA = targetBody.orbit.semiMajorAxis;
const a = 0.5 * (exBodyA + tgBodyA);
const period = Physics3D.orbitPeriod(attr, a);
infos.duration = lerp(0.35 * period, 0.75 * period, durationParam);
}
else {
const a = exitedBody.orbit.semiMajorAxis;
const period = Physics3D.orbitPeriod(attr, a);
infos.duration = lerp(1 * period, 2 * period, durationParam);
}
const { minLegDuration } = this.config;
infos.duration = Math.max(minLegDuration, infos.duration);
}
recomputeLegsSecondArcs() {
for (let i = 0; i < this._secondArcsData.length; i++) {
const data = this._secondArcsData[i];
const step = this.steps[3 + i * 3];
this._recomputeSecondArc(step, data);
}
}
_recomputeSecondArc(lambertStep, arcData) {
const fbBody = this.system[arcData.fbBodyId];
const { flybyOrbit, soiEnterAngle } = arcData;
const localEnterState = Physics3D.orbitElementsToState(flybyOrbit, fbBody, soiEnterAngle);
const { fbBodyState } = arcData;
const targetEnterPos = add3(fbBodyState.pos, localEnterState.pos);
const { preDSMState } = arcData;
const { v1, v2 } = Lambert.solve(preDSMState.pos, targetEnterPos, lambertStep.duration, this._mainAttractor);
const postDSMState = { pos: preDSMState.pos, vel: v1 };
const encounterState = { pos: targetEnterPos, vel: v2 };
const arcOrbit = Physics3D.stateToOrbitElements(postDSMState, this._mainAttractor);
const angles = {
begin: Physics3D.trueAnomalyFromOrbitalState(arcOrbit, postDSMState),
end: Physics3D.trueAnomalyFromOrbitalState(arcOrbit, encounterState)
};
const flybyOrbit = stateToOrbitElements(localIncomingState, body);
const isDestinationBody = body.id == this._destinationBody.id;
let nu1 = trueAnomalyFromOrbitalState(flybyOrbit, localIncomingState);
let nu2 = isDestinationBody ? 0 : -nu1;
const tof = tofBetweenAnomalies(flybyOrbit, body, nu1, nu2);
const drawAngles = {
begin: angles.begin,
end: angles.end
};
if (drawAngles.begin > drawAngles.end) {
drawAngles.end += TWO_PI;
}
lambertStep.orbitElts = arcOrbit;
lambertStep.angles = angles;
lambertStep.drawAngles = drawAngles;
const maneuvre = lambertStep.maneuvre;
const dsmDV = sub3(postDSMState.vel, preDSMState.vel);
maneuvre.deltaVToPrevStep = dsmDV;
}
_computeFlyby(flybyInfo) {
const body = this.system[flybyInfo.flybyBodyId];
const bodyVel = this._fbBodyState.vel;
const globalIncomingVel = this._vesselState.vel;
const relativeIncomingVel = sub3(globalIncomingVel, bodyVel);
const incomingVelMag = mag3(relativeIncomingVel);
const incomingVelDir = div3(relativeIncomingVel, incomingVelMag);
const normAngle = lerp(0, TWO_PI, flybyInfo.normAngleParam);
const t_normal = normalize3(cross(relativeIncomingVel, bodyVel));
const normal = rotate3(t_normal, incomingVelDir, normAngle);
const t_periDir = incomingVelDir;
const { fbRadiusMaxScale } = this.config;
const periRadius = lerp(body.radius, fbRadiusMaxScale * body.radius, flybyInfo.periRadiParam);
const t_periPos = mult3(t_periDir, periRadius);
const periVelMag = Physics3D.deduceVelocityAtRadius(body, body.soi, incomingVelMag, periRadius);
const t_periVelDir = rotate3(t_periDir, normal, HALF_PI);
const t_periVel = mult3(t_periVelDir, periVelMag);
const t_periapsisState = { pos: t_periPos, vel: t_periVel };
const t_flybyOrbit = Physics3D.stateToOrbitElements(t_periapsisState, body);
const exitAngle = Physics3D.trueAnomalyAtRadius(t_flybyOrbit, body.soi);
const angles = { begin: -exitAngle, end: exitAngle };
const drawAngles = { begin: angles.begin, end: angles.end };
const t_incomingVel = Physics3D.orbitElementsToState(t_flybyOrbit, body, angles.begin).vel;
const t_incomingVelDir = normalize3(t_incomingVel);
const dotCross = dot3(cross(incomingVelDir, t_incomingVelDir), normal);
const dotVel = dot3(incomingVelDir, t_incomingVelDir);
const offsetAngle = Math.atan2(dotCross, dotVel);
const periPos = rotate3(t_periPos, normal, -offsetAngle);
const periVel = rotate3(t_periVel, normal, -offsetAngle);
const periapsisState = { pos: periPos, vel: periVel };
const flybyOrbit = Physics3D.stateToOrbitElements(periapsisState, body);
const tof = Physics3D.tofBetweenAnomalies(flybyOrbit, body, angles.begin, angles.end);
this.steps.push({
orbitElts: flybyOrbit,
attractorId: body.id,
beginAngle: nu1,
endAngle: nu2,
angles: angles,
drawAngles: drawAngles,
duration: tof,
dateOfStart: this._missionTime
dateOfStart: this._lastStepEndDate
});
const exitState = Physics3D.orbitElementsToState(flybyOrbit, body, exitAngle);
this._vesselState = exitState;
this._secondArcsData.push({
preDSMState: this._preDSMState,
fbBodyId: body.id,
fbBodyState: this._fbBodyState,
flybyOrbit: flybyOrbit,
soiEnterAngle: angles.begin
});
}
_calculateLegSecondArc() {
const leg = this._legs[this._sequenceIndex - 1];
const preDSMState = this._preDSMState;
const origin = this.system[this.sequence[this._sequenceIndex - 1]];
const target = this.system[this.sequence[this._sequenceIndex]];
const startPos = preDSMState.pos;
const encounterDate = this._missionTime + leg.duration;
const targetState = getBodyStateAtDate(target, encounterDate, this._attractorMu);
const arcDuration = (1 - leg.dsmOffset) * leg.duration;
const incomingVel = solveLambert(startPos, targetState.pos, arcDuration, this._mainAttractor).v2;
const relativeVel = sub3(incomingVel, targetState.vel);
const incomingDir = normalize3(mult3(relativeVel, -1));
const { thetaMin, thetaMax } = minMaxRingAngle(target.soi, target.radius, target.radius * 2);
const pointIsValid = leg.theta >= thetaMin && leg.theta <= thetaMax;
if (!pointIsValid) {
const { theta, phi } = randomPointOnSphereRing(thetaMin, thetaMax);
leg.theta = theta;
leg.phi = phi;
_computeLegSecondArcSimple(legInfo) {
const attr = this._mainAttractor;
const lastStep = this._lastStep;
const preDSMState = this._vesselState;
const encounterDate = lastStep.dateOfStart + legInfo.duration;
const targetBody = this.system[legInfo.targetBodyId];
const tgBodyOrbit = this._bodiesOrbits[targetBody.id];
const tgBodyState = Physics3D.bodyStateAtDate(targetBody, tgBodyOrbit, attr, encounterDate);
const arcDuration = legInfo.duration * (1 - legInfo.dsmParam);
const { v1, v2 } = Lambert.solve(preDSMState.pos, tgBodyState.pos, arcDuration, attr);
const postDSMState = { pos: preDSMState.pos, vel: v1 };
const encounterState = { pos: tgBodyState.pos, vel: v2 };
const arcOrbit = Physics3D.stateToOrbitElements(postDSMState, attr);
const angles = {
begin: Physics3D.trueAnomalyFromOrbitalState(arcOrbit, postDSMState),
end: Physics3D.trueAnomalyFromOrbitalState(arcOrbit, encounterState)
};
const drawAngles = {
begin: angles.begin,
end: angles.end
};
if (drawAngles.begin > drawAngles.end) {
drawAngles.end += TWO_PI;
}
const soiPointLocal = mult3(pointOnSphereRing(incomingDir, leg.theta, leg.phi), target.soi);
const soiPoint = add3(targetState.pos, soiPointLocal);
const { v1, v2 } = solveLambert(startPos, soiPoint, arcDuration, this._mainAttractor);
const secondLegArc = stateToOrbitElements({ pos: startPos, vel: v1 }, this._mainAttractor);
let nu1 = trueAnomalyFromOrbitalState(secondLegArc, { pos: startPos, vel: v1 });
let nu2 = trueAnomalyFromOrbitalState(secondLegArc, { pos: soiPoint, vel: v2 });
if (nu1 > nu2) {
nu2 = TWO_PI + nu2;
}
this.totalDeltaV += mag3(sub3(v1, preDSMState.vel));
const progradeDir = normalize3(preDSMState.vel);
const dsmDV = sub3(postDSMState.vel, preDSMState.vel);
const maneuvre = {
deltaVToPrevStep: sub3(v1, preDSMState.vel),
progradeDir: normalize3(preDSMState.vel),
manoeuvrePosition: preDSMState.pos,
position: preDSMState.pos,
deltaVToPrevStep: dsmDV,
progradeDir: progradeDir,
context: {
type: "dsm",
originId: origin.id,
targetId: target.id
originId: legInfo.exitedBodyId,
targetId: legInfo.targetBodyId
}
};
this.steps.push({
orbitElts: secondLegArc,
attractorId: this._mainAttractor.id,
beginAngle: nu1,
endAngle: nu2,
orbitElts: arcOrbit,
attractorId: attr.id,
angles: angles,
drawAngles: drawAngles,
duration: arcDuration,
dateOfStart: this._lastStep.dateOfStart + this._lastStep.duration,
dateOfStart: this._lastStepEndDate,
maneuvre: maneuvre
});
this._flybyBodyState = targetState;
this._flybyIncomingState = { pos: soiPoint, vel: v2 };
this._vesselState = encounterState;
this._fbBodyState = tgBodyState;
}
_calculateLegFirstArc() {
const seqIndex = this._sequenceIndex;
const leg = this._legs[seqIndex - 1];
const isResonant = this.sequence[seqIndex - 1] == this.sequence[seqIndex];
const lastStep = this._lastStep;
const exitedBody = this.system[lastStep.attractorId];
const localState = stateFromOrbitElements(lastStep.orbitElts, exitedBody.stdGravParam, lastStep.endAngle);
const exitedBodyState = getBodyStateAtDate(exitedBody, this._missionTime, this._attractorMu);
_computeFirstLegArc(legInfo) {
const localExitState = this._vesselState;
const exitedBody = this.system[this._lastStep.attractorId];
const bodyOrbit = this._bodiesOrbits[exitedBody.id];
const exitDate = this._lastStepEndDate;
const exitedBodyState = Physics3D.bodyStateAtDate(exitedBody, bodyOrbit, this._mainAttractor, exitDate);
const exitState = {
pos: add3(localState.pos, exitedBodyState.pos),
vel: add3(localState.vel, exitedBodyState.vel)
pos: add3(exitedBodyState.pos, localExitState.pos),
vel: add3(exitedBodyState.vel, localExitState.vel),
};
const firstLegArc = stateToOrbitElements(exitState, this._mainAttractor);
if (firstLegArc.semiMajorAxis >= this._mainAttractor.soi) {
return false;
const arcOrbit = Physics3D.stateToOrbitElements(exitState, this._mainAttractor);
const beginAngle = Physics3D.trueAnomalyFromOrbitalState(arcOrbit, exitState);
const arcDuration = legInfo.dsmParam * legInfo.duration;
const preDSMState = Physics3D.propagateStateFromTrueAnomaly(arcOrbit, this._mainAttractor, beginAngle, arcDuration);
let endAngle = Physics3D.trueAnomalyFromOrbitalState(arcOrbit, preDSMState);
const angles = { begin: beginAngle, end: endAngle };
const drawAngles = { begin: angles.begin, end: angles.end };
const period = Physics3D.orbitPeriod(this._mainAttractor, arcOrbit.semiMajorAxis);
if (arcDuration > period) {
drawAngles.begin = 0;
drawAngles.end = TWO_PI;
}
const { dsmOffsetMax, dsmOffsetMin, minLegDuration } = this.config;
if (firstLegArc.eccentricity < 1) {
const orbitPeriod = calculateOrbitPeriod(firstLegArc.semiMajorAxis, this._attractorMu);
let minScale = 0.5;
let maxScale = 1;
if (isResonant) {
minScale = 1;
maxScale = 3;
}
const minDuration = minScale * orbitPeriod;
const maxDuration = maxScale * orbitPeriod;
leg.duration = clamp(leg.duration, minDuration, maxDuration);
const revolutions = leg.duration / orbitPeriod;
const minDSMOffset = Math.max((revolutions - 1) / revolutions, dsmOffsetMin);
leg.dsmOffset = clamp(leg.dsmOffset, minDSMOffset, dsmOffsetMax);
}
leg.duration = clamp(leg.duration, minLegDuration, Infinity);
leg.dsmOffset = clamp(leg.dsmOffset, dsmOffsetMin, dsmOffsetMax);
const arcDuration = leg.dsmOffset * leg.duration;
let legStartNu = trueAnomalyFromOrbitalState(firstLegArc, exitState);
const preDSMState = propagateState(firstLegArc, legStartNu, arcDuration, this._mainAttractor);
let dsmNu = trueAnomalyFromOrbitalState(firstLegArc, preDSMState);
if (legStartNu > dsmNu) {
dsmNu += TWO_PI;
}
if (isNaN(legStartNu) || isNaN(dsmNu)) {
return false;
else if (beginAngle > endAngle) {
drawAngles.end += TWO_PI;
}
this.steps.push({
orbitElts: firstLegArc,
orbitElts: arcOrbit,
attractorId: this._mainAttractor.id,
beginAngle: legStartNu,
endAngle: dsmNu,
angles: angles,
drawAngles: drawAngles,
duration: arcDuration,
dateOfStart: this._lastStep.dateOfStart + this._lastStep.duration
dateOfStart: exitDate,
});
this._vesselState = preDSMState;
this._preDSMState = preDSMState;
return true;
}
_calculateEjectionOrbit() {
const startBody = this.system[this.sequence[0]];
const nextBody = this.system[this.sequence[1]];
const parkingRadius = startBody.radius + this.depAltitude;
const parkingVel = circularVelocity(startBody, parkingRadius);
const ejectionVel = this._departureDVScale * ejectionVelocity(startBody, parkingRadius);
const ejectionDV = ejectionVel - parkingVel;
this.totalDeltaV += ejectionDV;
const bodyPos = getBodyStateAtDate(startBody, this._departureDate, this._attractorMu).pos;
const r1 = startBody.orbit.semiMajorAxis;
const r2 = nextBody.orbit.semiMajorAxis;
let { tangent, nu } = idealEjectionDirection(bodyPos, r2 < r1);
const ejectionState = {
pos: vec3(parkingRadius * Math.cos(nu), 0, parkingRadius * Math.sin(nu)),
vel: mult3(tangent, ejectionVel),
};
const offsetAngle = hyperbolicEjectionOffsetAngle(ejectionVel, parkingRadius, startBody);
const up = vec3(0, 1, 0);
ejectionState.pos = rotate3(ejectionState.pos, up, -offsetAngle);
ejectionState.vel = rotate3(ejectionState.vel, up, -offsetAngle);
const ejectionOrbit = stateToOrbitElements(ejectionState, startBody);
const exitAnomaly = soiExitTrueAnomaly(ejectionOrbit, startBody.soi);
const tof = tofBetweenAnomalies(ejectionOrbit, startBody, 0, exitAnomaly);
const progradeDir = normalize3(ejectionState.vel);
_computeDepartureEjectionOrbit() {
const curOrbit = this._lastStep.orbitElts;
const depBody = this._departureBody;
const { phaseParam, ejVelParam } = this._departureInfos;
const phase = lerp(0, TWO_PI, phaseParam);
const ejVelMag0 = Physics3D.velocityToReachAltitude(depBody, curOrbit.semiMajorAxis, depBody.soi);
const { depDVScaleMin, depDVScaleMax } = this.config;
const ejVelMag = lerp(depDVScaleMin, depDVScaleMax, ejVelParam) * ejVelMag0;
const { vel, pos } = Physics3D.orbitElementsToState(curOrbit, depBody, phase);
const progradeDir = normalize3(vel);
const ejVel = mult3(progradeDir, ejVelMag);
const ejOrbit = Physics3D.stateToOrbitElements({ pos, vel: ejVel }, depBody);
const soiExitAngle = Physics3D.trueAnomalyAtRadius(ejOrbit, depBody.soi);
const ejDV = ejVelMag - mag3(vel);
const maneuvre = {
deltaVToPrevStep: mult3(progradeDir, ejectionDV),
position: pos,
deltaVToPrevStep: mult3(progradeDir, ejDV),
progradeDir: progradeDir,
manoeuvrePosition: ejectionState.pos,
context: { type: "ejection" }
};
const tof = Physics3D.tofBetweenAnomalies(ejOrbit, depBody, 0, soiExitAngle);
this.steps.push({
orbitElts: ejectionOrbit,
attractorId: startBody.id,
beginAngle: 0,
endAngle: exitAnomaly,
orbitElts: ejOrbit,
attractorId: depBody.id,
angles: { begin: 0, end: soiExitAngle },
drawAngles: { begin: 0, end: soiExitAngle },
duration: tof,
dateOfStart: this._lastStepDateOfEnd,
dateOfStart: this._lastStepEndDate,
maneuvre: maneuvre
});
const exitState = Physics3D.orbitElementsToState(ejOrbit, depBody, soiExitAngle);
this._vesselState = exitState;
}
_calculateParkingOrbit() {
const startBody = this.system[this.sequence[0]];
const orbitRadius = startBody.radius + this.depAltitude;
_appendDepartureParkingOrbit() {
const { dateParam } = this._departureInfos;
const radius = this._departureBody.radius + this._depAltitude;
const dateMin = this._startDateMin;
const dateMax = this._startDateMax;
this.steps.push({
orbitElts: equatorialCircularOrbit(orbitRadius),
attractorId: startBody.id,
beginAngle: 0,
endAngle: TWO_PI,
orbitElts: Physics3D.equatorialCircularOrbit(radius),
attractorId: this._departureBody.id,
angles: { begin: 0, end: 0 },
drawAngles: { begin: 0, end: TWO_PI },
duration: 0,
dateOfStart: this._departureDate
dateOfStart: lerp(dateMin, dateMax, dateParam)
});
}
}

2
dist/dedicated-workers/sequence-evaluator.js vendored
View File

@@ -99,4 +99,4 @@ class SequenceEvaluator extends WorkerEnvironment {
}
}
}
initWorker(SequenceEvaluator);
WorkerEnvironment.init(SequenceEvaluator);

2
dist/dedicated-workers/sequence-generator.js vendored
View File

@@ -85,4 +85,4 @@ class SequenceGenerator extends WorkerEnvironment {
info.resonant++;
}
}
initWorker(SequenceGenerator);
WorkerEnvironment.init(SequenceGenerator);

148
dist/dedicated-workers/trajectory-optimizer.js vendored
View File

@@ -4,6 +4,12 @@ class TrajectoryOptimizer extends WorkerEnvironment {
onWorkerInitialize(data) {
this._config = data.config;
this._system = data.system;
this._bodiesOrbits = [null];
for (let i = 1; i < this._system.length; i++) {
const data = this._system[i].orbit;
const orbit = Physics3D.orbitElementsFromOrbitData(data);
this._bodiesOrbits.push(orbit);
}
}
onWorkerDataPass(data) {
this._depAltitude = data.depAltitude;
@@ -12,101 +18,87 @@ class TrajectoryOptimizer extends WorkerEnvironment {
this._startDateMax = data.startDateMax;
}
onWorkerRun(input) {
const evaluate = (params) => this._evaluate(params);
if (input.start) {
const numLegs = this._sequence.length - 1;
const agentDim = 3 + numLegs * 4 - 2;
const fitness = (agent) => {
const trajectory = this._computeTrajectory(agent);
if (trajectory.totalDeltaV < this._bestDeltaV) {
this._bestDeltaV = trajectory.totalDeltaV;
this._bestTrajectory = trajectory;
}
return trajectory.totalDeltaV;
};
const trajConfig = this._config.trajectorySearch;
const { crossoverProba, diffWeight } = trajConfig;
const { chunkStart, chunkEnd } = input;
this._evolver = new ChunkedEvolver(chunkStart, chunkEnd, agentDim, fitness, crossoverProba, diffWeight);
this._bestDeltaV = Infinity;
const chunkSize = input.chunkEnd - input.chunkStart + 1;
const popChunk = this._createRandomMGAPopulationChunk(chunkSize);
const dvChunk = populationChunkFitness(popChunk, evaluate);
const popChunk = this._evolver.createRandomPopulationChunk();
const dvChunk = this._evolver.evaluateChunkFitness(popChunk);
sendResult({
bestSteps: this._bestSteps,
bestSteps: this._bestTrajectory.steps,
bestDeltaV: this._bestDeltaV,
fitChunk: dvChunk,
popChunk: popChunk,
});
}
else {
const { crossoverProba, diffWeight } = this._config.trajectorySearch;
const results = evolvePopulationChunk(input.population, input.deltaVs, input.chunkStart, input.chunkEnd, crossoverProba, diffWeight, evaluate);
const { population, deltaVs } = input;
const { popChunk, fitChunk } = this._evolver.evolvePopulationChunk(population, deltaVs);
sendResult({
...results,
bestSteps: this._bestSteps,
popChunk, fitChunk,
bestSteps: this._bestTrajectory.steps,
bestDeltaV: this._bestDeltaV
});
}
}
_evaluate(params) {
const trajectory = this._computeTrajectory(params);
if (trajectory.totalDeltaV < this._bestDeltaV) {
this._bestDeltaV = trajectory.totalDeltaV;
this._bestSteps = trajectory.steps;
_computeTrajectory(agent, maxAttempts = 1000) {
const trajConfig = this._config.trajectorySearch;
const trajectory = new TrajectoryCalculator(this._system, trajConfig, this._sequence);
trajectory.addPrecomputedOrbits(this._bodiesOrbits);
let attempts = 0;
while (attempts < maxAttempts) {
trajectory.setParameters(this._depAltitude, this._startDateMin, this._startDateMax, agent);
let failed = false;
try {
trajectory.compute();
trajectory.recomputeLegsSecondArcs();
}
catch {
failed = true;
}
if (failed || this._hasNaNValuesInSteps(trajectory)) {
this._evolver.randomizeAgent(agent);
trajectory.reset();
}
else {
return trajectory;
}
attempts++;
}
return trajectory.totalDeltaV;
throw new Error("Impossible to compute the trajectory.");
}
;
_computeTrajectory(params) {
const calculate = () => {
const trajectory = new TrajectoryCalculator(this._system, this._config.trajectorySearch, this._sequence);
trajectory.setParameters(this._depAltitude, this._startDateMin, this._startDateMax, params);
trajectory.compute();
return trajectory;
_hasNaNValuesInSteps(trajectory) {
const hasNaN = obj => {
for (const value of Object.values(obj)) {
if (typeof value == "object") {
if (hasNaN(value))
return true;
}
else if (typeof value == "number") {
if (isNaN(value))
return true;
}
}
return false;
};
let trajectory = calculate();
while (!trajectory.noError) {
this._randomizeExistingAgent(params);
trajectory = calculate();
}
return trajectory;
}
_createRandomMGAPopulationChunk(chunkSize) {
const popChunk = [];
for (let i = 0; i < chunkSize; i++) {
popChunk.push(this._createRandomMGAAgent());
}
return popChunk;
}
_createRandomMGAAgent() {
const dim = 4 * (this._sequence.length - 1) + 2;
const solution = Array(dim).fill(0);
solution[0] = randomInInterval(this._startDateMin, this._startDateMax);
solution[1] = randomInInterval(1, 3);
for (let i = 1; i < this._sequence.length; i++) {
const j = 2 + (i - 1) * 4;
const { legDuration, dsmOffset } = this._legDurationAndDSM(i);
solution[j] = legDuration;
solution[j + 1] = dsmOffset;
const { theta, phi } = randomPointOnSphereRing(0, Math.PI / 2);
solution[j + 2] = theta;
solution[j + 3] = phi;
}
return solution;
}
_randomizeExistingAgent(agent) {
const newAgent = this._createRandomMGAAgent();
for (let i = 0; i < agent.length; i++) {
agent[i] = newAgent[i];
}
}
_legDurationAndDSM(index) {
const isResonant = this._sequence[index - 1] == this._sequence[index];
const { dsmOffsetMin, dsmOffsetMax } = this._config.trajectorySearch;
if (isResonant) {
const sideralPeriod = this._system[this._sequence[index]].orbit.sideralPeriod;
const revs = randomInInterval(1, 4);
const legDuration = revs * sideralPeriod;
const minOffset = clamp((revs - 1) / revs, dsmOffsetMin, dsmOffsetMax);
const dsmOffset = randomInInterval(minOffset, dsmOffsetMax);
return { legDuration, dsmOffset };
}
else {
const body1 = this._system[this._sequence[index - 1]];
const body2 = this._system[this._sequence[index]];
const attractor = this._system[body1.orbiting];
const period = getHohmannPeriod(body1, body2, attractor);
const legDuration = randomInInterval(0.1, 1) * period;
const dsmOffset = randomInInterval(dsmOffsetMin, dsmOffsetMax);
return { legDuration, dsmOffset };
const { steps } = trajectory;
for (let i = steps.length - 1; i >= 0; i--) {
if (hasNaN(steps[i]))
return true;
}
return false;
}
}
initWorker(TrajectoryOptimizer);
WorkerEnvironment.init(TrajectoryOptimizer);

2
dist/main/editor/editor.js vendored
View File

@@ -124,7 +124,7 @@ export function initEditor(controls, system, config, canvas) {
const displayFoundTrajectory = () => {
trajectory = new Trajectory(solver.bestTrajectorySteps, system, config);
trajectory.draw(canvas);
trajectory.fillResultControls(maneuvreSelector, resultSpans, stepSlider, systemTime);
trajectory.fillResultControls(maneuvreSelector, resultSpans, stepSlider, systemTime, controls);
maneuvreSelector.select(0);
maneuvreSelector.enable();
stepSlider.enable();

103
dist/main/solvers/trajectory-solver.js vendored
View File

@@ -17,49 +17,41 @@ export class TrajectorySolver {
this._workerPool = new WorkerPool("dist/dedicated-workers/trajectory-optimizer.js", this.config);
this._workerPool.initialize({ system: this.system.data, config: this.config });
}
_initPlot() {
this.plot.clearPlot();
}
_updatePlot(iteration) {
const { best, mean } = this._bestMeanDeltaV;
let mean = 0;
for (const dv of this._deltaVs)
mean += dv;
mean /= this.popSize;
const best = this.bestDeltaV;
this.plot.addIterationData(iteration, mean, best);
}
cancel() {
if (this._running)
this._cancelled = true;
}
_checkCancellation() {
if (this._cancelled) {
this._cancelled = false;
this._running = false;
throw new Error("TRAJECTORY FINDER CANCELLED");
}
}
async searchOptimalTrajectory(sequence, startDateMin, startDateMax, depAltitude) {
this._running = true;
this._initPlot();
this.plot.clearPlot();
this._calculatePopulationSize(sequence);
this._calculatePopulationChunks();
await this._passSettingsData(sequence.ids, startDateMin, startDateMax, depAltitude);
await this._createStartPopulation();
this._updatePlot(0);
this._checkCancellation();
const { maxGenerations } = this.config.trajectorySearch;
for (let i = 0; i < maxGenerations; i++) {
if (this._cancelled) {
this._cancelled = false;
this._running = false;
throw new Error("TRAJECTORY FINDER CANCELLED");
}
await this._generateNextPopulation();
this._updatePlot(1 + i);
this._checkCancellation();
}
this._running = false;
}
async _passSettingsData(sequence, startDateMin, startDateMax, depAltitude) {
return this._workerPool.passData({ depAltitude, sequence, startDateMin, startDateMax });
}
async _createStartPopulation() {
const inputs = this._firstGenerationInputs();
const results = await this._workerPool.runPool(inputs);
this._mergeResultsChunks(results);
}
_calculatePopulationChunks() {
const { splitLimit } = this.config.trajectorySearch;
const numChunks = this._workerPool.optimizeUsedWorkersCount(this.popSize, splitLimit);
@@ -74,15 +66,35 @@ export class TrajectorySolver {
this._numChunks = numChunks;
this._chunkIndices = chunkIndices;
}
_calculatePopulationSize(sequence) {
const { popSizeDimScale } = this.config.trajectorySearch;
this.popSize = popSizeDimScale * (4 * (sequence.length - 1) + 2);
}
async _generateNextPopulation() {
const inputs = this._nextGenerationInputs();
async _createStartPopulation() {
const inputs = [];
for (let i = 0; i < this._numChunks; i++) {
const start = this._chunkIndices[i * 2];
const end = this._chunkIndices[i * 2 + 1];
inputs.push({
start: true,
chunkStart: start,
chunkEnd: end
});
}
const results = await this._workerPool.runPool(inputs);
this._mergeResultsChunks(results);
}
async _generateNextPopulation() {
const inputs = [];
for (let i = 0; i < this._numChunks; i++) {
inputs[i] = {
population: this._population,
deltaVs: this._deltaVs,
};
}
const results = await this._workerPool.runPool(inputs);
this._mergeResultsChunks(results);
}
_calculatePopulationSize(sequence) {
const { popSizeDimScale } = this.config.trajectorySearch;
this.popSize = popSizeDimScale * sequence.length;
}
_mergeResultsChunks(results) {
const popChunks = [];
const dVChunks = [];
@@ -102,45 +114,4 @@ export class TrajectorySolver {
this.bestTrajectorySteps = bestSteps;
this.bestDeltaV = bestDeltaV;
}
_firstGenerationInputs() {
const inputs = [];
for (let i = 0; i < this._numChunks; i++) {
const { start, end } = this._chunkStartEnd(i);
inputs.push({
start: true,
chunkStart: start,
chunkEnd: end
});
}
return inputs;
}
_nextGenerationInputs() {
const inputs = [];
for (let i = 0; i < this._numChunks; i++) {
const { start, end } = this._chunkStartEnd(i);
inputs[i] = {
population: this._population,
deltaVs: this._deltaVs,
chunkStart: start,
chunkEnd: end
};
}
return inputs;
}
_chunkStartEnd(index) {
return {
start: this._chunkIndices[index * 2],
end: this._chunkIndices[index * 2 + 1]
};
}
get _bestMeanDeltaV() {
let mean = 0;
for (const dv of this._deltaVs)
mean += dv;
mean /= this.popSize;
return {
mean: mean,
best: this.bestDeltaV
};
}
}

20
dist/main/solvers/trajectory.js vendored
View File

@@ -35,8 +35,8 @@ export class Trajectory {
const { scale } = this.config.rendering;
for (let i = 0; i < this.orbits.length; i++) {
const orbit = this.orbits[i];
const { beginAngle, endAngle } = this.steps[i];
const orbitPoints = createOrbitPoints(orbit, samplePoints, scale, beginAngle, endAngle);
const { begin, end } = this.steps[i].drawAngles;
const orbitPoints = createOrbitPoints(orbit, samplePoints, scale, begin, end);
const color = new THREE.Color(`hsl(${i * 35 % 360}, 100%, 85%)`);
const orbitLine = createLine(orbitPoints, resolution, {
color: color.getHex(),
@@ -54,7 +54,7 @@ export class Trajectory {
if (step.maneuvre) {
const group = this.system.objectsOfBody(step.attractorId);
const sprite = createSprite(Trajectory.arrowMaterial, 0xFFFFFF, false, maneuvreArrowSize);
const { x, y, z } = step.maneuvre.manoeuvrePosition;
const { x, y, z } = step.maneuvre.position;
sprite.position.set(x, y, z);
sprite.position.multiplyScalar(scale);
group.add(sprite);
@@ -84,14 +84,15 @@ export class Trajectory {
}
}
}
fillResultControls(maneuvreSelector, resultSpans, stepSlider, systemTime) {
fillResultControls(maneuvreSelector, resultSpans, stepSlider, systemTime, controls) {
const depDate = new KSPTime(this.steps[0].dateOfStart, this.config.time);
resultSpans.totalDVSpan.innerHTML = this._totalDeltaV.toFixed(1);
resultSpans.depDateSpan.innerHTML = depDate.stringYDHMS("hms", "ut");
resultSpans.depDateSpan.onclick = () => {
this.system.date = depDate.dateSeconds;
controls.centerOnTarget();
systemTime.time.dateSeconds = depDate.dateSeconds;
systemTime.update();
systemTime.onChange();
};
stepSlider.setMinMax(0, this.steps.length - 1);
stepSlider.input((index) => this._displayStepsUpTo(index));
@@ -128,11 +129,16 @@ export class Trajectory {
resultSpans.radialDVSpan.innerHTML = details.radialDV.toFixed(1);
resultSpans.maneuvreNumber.innerHTML = (index + 1).toString();
resultSpans.dateSpan.onclick = () => {
systemTime.time.dateSeconds = depDate.dateSeconds + dateEMT.dateSeconds;
const date = depDate.dateSeconds + dateEMT.dateSeconds;
this.system.date = date;
controls.centerOnTarget();
systemTime.time.dateSeconds = date;
systemTime.update();
systemTime.onChange();
};
});
for (const step of this.steps) {
console.log(step);
}
}
_displayStepsUpTo(index) {
for (let i = 0; i < this.steps.length; i++) {

BIN
imgs/jool_system_trajectory.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 22 KiB

29
index.html
View File

@@ -280,7 +280,7 @@
<h3>Planetary sequence</h3>
<p>
The flyby sequence is the sequence of bodies encountered during the mission's
flight time. It starts with the departure body and ends with the arrival body.
flight time. It starts with the departure body and ends with the destination body.
All intermediate bodies of the sequence are used for gravity assist. <br>
Sequences are represented using the first two characters of each body's name. For example, a mission
going from Kerbin to Jool, with a first flyby of Eve followed by a flyby of Duna would be written:
@@ -299,12 +299,19 @@
<h3>Sequence generation</h3>
<p>
The first step of the trajectory planning consists of generating a <em>possible</em> planetary sequence.
The first step of the trajectory planning consists of generating a planetary sequence.
It will list a bunch of <em>possible</em> sequences considering a very simplified model of the solar system
and swing-bys, without any phasing.
Therefore there is <strong>no guarantee</strong> to have the most optimal sequence generated. But it tries to
provide feasible sequences. <br>
The following settings must be defined:
</p>
<p>
Alternatively, it is possible to set a custom sequence if no interesting sequences is proposed by the generator,
or to define custom routes. The input text must be a valid sequence representation.
The trajectory calculation step will choose the custom sequence first if it is defined.
</p>
<p>
For the sequence generator, the following settings must be defined:
</p>
<ul>
<li>The <strong>origin</strong> and <strong>destination</strong> body of the trajectory.</li>
@@ -332,11 +339,6 @@
this back leg has a spacing of 1.
</li>
</ul>
<p>
Alternatively, it is possible to set a custom sequence if no interesting sequences is proposed by the generator,
or to define custom routes. The input text must be a valid sequence representation.
The trajectory calculation step will choose the custom sequence first if defined.
</p>
<h3>Trajectory calculation</h3>
<p>
@@ -381,13 +383,14 @@
and then choose the best one.
</p>
<p>
Currently, although the tool allows it, searching trajectories between moons (e.g. in Jool's system)
may give absurd trajectories due to implementation details.
The computed trajectory itself may not be feasible in game because
the current implementation of the application may not consider some in game particularities.
Therefore the maneuvers' details are more indicative and fine tunings are necessary.
</p>
<p>
The computed trajectory itself may not be feasible in game because
it follows the current implementation of the application and may not consider some in game particularities.
Therefore the maneuvers' details are more indicative and fine tunings are necessary.
For the Jool's system, since the SOIs of the satellites are large compared to the radius of their
orbits, it is possible that a produced trajectory is not feasible because it intersects a satellites's SOI
before reaching the final state of an arc.
</p>
</article>
<!-- End of Issues section -->

118
src/dedicated-workers/libs/common.ts
View File

@@ -1,49 +1,101 @@
// Common file for all worker scripts. It adds useful communication function
// with the main thread to be automatically handled by the ComputeWorker class.
/**
* Type enforced postMessage.
* @param msg A message
*/
function postMessageSafe(msg: MessageFromWorker) {
postMessage(msg);
}
/**
* Sends the progress status back to the main thread.
* @param progress The progress
* @param data Data to pass along
*/
function sendProgress(progress: number, data?: any){
postMessageSafe({label: "progress", progress, data});
}
/**
* Sends data to be shown in the console. Use it as a
* console.log
* @param data Data to be debuged
*/
function debug(...data: any[]){
postMessageSafe({label: "debug", data: data});
}
/**
* Sends back a result object to the main thread,
* indicating an end of run.
* @param result The results of the ComputeWorker to send back
*/
function sendResult(result: any) {
postMessageSafe({label: "complete", result: result});
}
function postMessageSafe(msg: MessageFromWorker) {
postMessage(msg);
}
/**
* Encapsulates the necessary functions for all ComputeWorkers to use.
*/
class WorkerEnvironment {
public onWorkerInitialize(data: any) {}
public onWorkerRun(input?: any) {}
public onWorkerContinue(input?: any) {}
public onWorkerStop() {}
public onWorkerDataPass(data: any) {}
}
function initWorker(Env: typeof WorkerEnvironment){
const env = new Env();
onmessage = ({data}: MessageEvent<MessageToWorker>) => {
switch(data.label) {
case "initialize":
env.onWorkerInitialize(data.config);
postMessageSafe({label: "initialized"});
break;
case "run":
env.onWorkerRun(data.input);
break;
case "continue":
env.onWorkerContinue();
break;
case "stop":
env.onWorkerStop();
postMessageSafe({label: "stopped"});
break;
case "pass":
env.onWorkerDataPass(data.data);
postMessageSafe({label: "received"});
break;
/**
* Initializes the events for the worker to listen and react to
* @param Env The environment type, inheriting from `WorkerEnvironment`
*/
public static init(Env: typeof WorkerEnvironment) {
const env = new Env();
onmessage = ({data}: MessageEvent<MessageToWorker>) => {
switch(data.label) {
case "initialize":
env.onWorkerInitialize(data.config);
postMessageSafe({label: "initialized"});
break;
case "run":
env.onWorkerRun(data.input);
break;
case "continue":
env.onWorkerContinue();
break;
case "stop":
env.onWorkerStop();
postMessageSafe({label: "stopped"});
break;
case "pass":
env.onWorkerDataPass(data.data);
postMessageSafe({label: "received"});
break;
}
}
}
/**
* Called on an initialize message.
* @param data The received initalization data
*/
public onWorkerInitialize(data: any) {}
/**
* Called on a run message.
* @param input The input provided for this run instance.
*/
public onWorkerRun(input?: any) {}
/**
* Called on a continue message.
* @param input
*/
public onWorkerContinue(input?: any) {}
/**
* Called on a stop message.
*/
public onWorkerStop() {}
/**
* Called on a data pass message.
* @param data The data passed to the worker
*/
public onWorkerDataPass(data: any) {}
}

229
src/dedicated-workers/libs/evolution.ts
View File

@@ -1,103 +1,156 @@
/**
* @param popChunk The population chunk
* @param fitness The fitness function
* @returns The fitness value for each agent in the chunk.
* Encapsulates a set of functions and variables to evolve a particular population chunk.
*/
function populationChunkFitness(popChunk: Agent[], fitness: (x: Agent) => number){
const fitChunk: number[] = Array(popChunk.length).fill(0);
for(let i = 0; i < popChunk.length; i++) {
fitChunk[i] = fitness(popChunk[i]);
class ChunkedEvolver {
public readonly chunkSize!: number;
/**
* Instantiates an evolver object, used to evolve a specified chunk of a whole population.
* @param chunkStart The start index of the chunk in the population
* @param chunkEnd The end index (inclusive) of the chunk in the population
* @param agentDim The dimension size of agents
* @param fitness The fitness function
* @param cr The CR coefficient of the DE
* @param f The F coefficient of the DE
*/
constructor(
public readonly chunkStart: number,
public readonly chunkEnd: number,
public readonly agentDim: number,
public readonly fitness: (x: Agent) => number,
public readonly cr: number,
public readonly f: number,
){
const ce = chunkEnd;
const cs = chunkStart;
this.chunkSize = ce - cs + 1;
}
return fitChunk;
}
/**
* Evolves a chunk of population from generation G to generation G + 1.
* @param population The whole population at the generation G
* @param fitnesses The fitnesses of all agents in the population at generation G
* @param chunkStart The first index of the chunk to evolve
* @param chunkEnd The last index of the chunk to evolve
* @param cr The CR coefficient of DE
* @param f The F coefficient of DE
* @param fitness The fitness function
* @returns The new population chunk with it's fitness at generation G + 1
*/
function evolvePopulationChunk(
population: Agent[],
fitnesses: number[],
chunkStart: number,
chunkEnd: number,
cr: number,
f: number,
fitness: (x: Agent) => number
){
const chunkSize = chunkEnd - chunkStart + 1;
const dim = population[0].length;
/**
* Creates a random agent with its parameters all set to random values
* between 0 and 1.
* @returns A random agent
*/
public createRandomAgent(){
const agent: Agent = Array<number>(this.agentDim);
this.randomizeAgent(agent);
return agent;
}
const nextPopChunk: (Agent | null)[] = Array(chunkSize).fill(null);
const nextFitChunk: number[] = Array(chunkSize).fill(0);
/**
* Randomizes an existing agent by setting its parameters to random values
* between 0 and 1.
* @param agent The agent to randomize
*/
public randomizeAgent(agent: Agent){
for(let i = 0; i < this.agentDim; i++){
agent[i] = Math.random();
}
}
for(let j = chunkStart; j <= chunkEnd; j++) {
const x = population[j];
const fx = fitnesses[j];
const [a, b, c] = pick3(population, j);
const ri = randint(0, dim - 1);
const y: Agent = Array(dim).fill(0);
for(let i = 0; i < dim; i++){
if(Math.random() < cr || i == ri) {
y[i] = a[i] + f * (b[i] - c[i]);
/**
* Returns a population chunk with agents whose parameters have been set to a random
* number between 0 and 1.
* @returns A random population chunk
*/
public createRandomPopulationChunk(){
const popChunk: Agent[] = [];
for(let i = 0; i < this.chunkSize; i++){
popChunk.push(this.createRandomAgent());
}
return popChunk;
}
/**
* Returns the fitness of each agent in the chunk.
* @param popChunk The population chunk
* @param fitness The fitness function
* @returns The fitness value for each agent in the chunk.
*/
public evaluateChunkFitness(popChunk: Agent[]){
const fitChunk: number[] = Array(popChunk.length);
for(let i = 0; i < popChunk.length; i++) {
fitChunk[i] = this.fitness(popChunk[i]);
}
return fitChunk;
}
/**
* Evolves a chunk of population from generation G to generation G + 1.
* @param population The whole population at the generation G
* @param fitnesses The fitnesses of all agents in the population at generation G
* @returns The new population chunk with it's fitness at generation G + 1
*/
public evolvePopulationChunk(population: Agent[], fitnesses: number[]){
const dim = population[0].length;
const nextPopChunk: Agent[] = Array(this.chunkSize);
const nextFitChunk: number[] = Array(this.chunkSize);
for(let j = this.chunkStart; j <= this.chunkEnd; j++) {
const x = population[j];
const fx = fitnesses[j];
const [a, b, c] = this._pick3(population, j);
const ri = randint(0, dim - 1);
const y: Agent = Array(dim).fill(0);
for(let i = 0; i < dim; i++){
if(Math.random() < this.cr || i == ri) {
y[i] = a[i] + this.f * (b[i] - c[i]);
} else {
y[i] = x[i];
}
y[i] = clamp(y[i], 0, 1);
}
const fy = this.fitness(y);
const index = j - this.chunkStart;
if(fy < fx) {
nextPopChunk[index] = y;
nextFitChunk[index] = fy;
} else {
y[i] = x[i];
nextPopChunk[index] = x;
nextFitChunk[index] = fx;
}
}
const fy = fitness(y);
const index = j - chunkStart;
if(fy < fx) {
nextPopChunk[index] = y;
nextFitChunk[index] = fy;
} else {
nextPopChunk[index] = x;
nextFitChunk[index] = fx;
return {
popChunk: nextPopChunk,
fitChunk: nextFitChunk
}
}
return {
popChunk: nextPopChunk as Agent[],
fitChunk: nextFitChunk
/**
* Selects three random agents in the population that are distinct from each other
* and distinct from the parent agent.
* @param population The whole population
* @param parentIndex The index of the parent, index not to pick
* @returns Three random distincy agents from the population.
*/
private _pick3(population: Agent[], parentIndex: number) {
const swap = (arr: any[], i: number, j: number) => {
const t = arr[j];
arr[j] = arr[i];
arr[i] = t;
}
swap(population, parentIndex, 0);
const picked: (Agent | null)[] = [null, null, null];
const pickedIndices = [0, 0, 0];
for(let i = 0; i <= 2; i++) {
const ri = randint(1 + i, population.length - 1);
picked[i] = population[ri];
pickedIndices[i] = ri;
swap(population, ri, i);
}
for(let i = 2; i >= 0; i--){
swap(population, pickedIndices[i], i);
}
swap(population, parentIndex, 0);
return picked as Agent[];
}
}
/**
* Selects three random agents in the population that are distinct from each other
* and distinct from the parent agent, in a O(1) time complexity.
* @param population The whole population
* @param parentIndex The index of the parent, index not to pick
* @returns Three random distincy agents from the population.
*/
function pick3(population: Agent[], parentIndex: number) {
swap(population, parentIndex, 0);
const picked: (Agent | null)[] = [null, null, null];
const pickedIndices = [0, 0, 0];
for(let i = 0; i <= 2; i++) {
const ri = randint(1 + i, population.length - 1);
picked[i] = population[ri];
pickedIndices[i] = ri;
swap(population, ri, i);
}
for(let i = 2; i >= 0; i--){
swap(population, pickedIndices[i], i);
}
swap(population, parentIndex, 0);
return picked as Agent[];
}
function swap(arr: any[], i: number, j: number){
const t = arr[j];
arr[j] = arr[i];
arr[i] = t;
}

349
src/dedicated-workers/libs/lambert.ts
View File

@@ -26,187 +26,192 @@
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
*****************************************************************************/
/**
* Solves the Lambert's problem considering 0 revolutions and prograde direction only.
* @param r1vec The start position in space
* @param r2vec The end position in space
* @param tof The time of flight between the two positions
* @param attractor The attractor body
* @returns The velocities at each points
*/
function solveLambert(r1vec: Vector3, r2vec: Vector3, tof: number, attractor: ICelestialBody) {
const mu = attractor.stdGravParam;
// Calculating lambda and T
const r1 = mag3(r1vec);
const r2 = mag3(r2vec);
const c = mag3(sub3(r2vec, r1vec));
namespace Lambert
{
/**
* Solves the Lambert's problem considering 0 revolutions and prograde direction only.
* @param r1vec The start position in space
* @param r2vec The end position in space
* @param tof The time of flight between the two positions
* @param attractor The attractor body
* @returns The velocities at each points
*/
export function solve(r1vec: Vector3, r2vec: Vector3, tof: number, attractor: ICelestialBody) {
const mu = attractor.stdGravParam;
const s = 0.5 * (r1 + r2 + c);
// Calculating lambda and T
const r1 = mag3(r1vec);
const r2 = mag3(r2vec);
const c = mag3(sub3(r2vec, r1vec));
const ir1 = div3(r1vec, r1);
const ir2 = div3(r2vec, r2);
const s = 0.5 * (r1 + r2 + c);
const ih = normalize3(cross(ir1, ir2));
const ir1 = div3(r1vec, r1);
const ir2 = div3(r2vec, r2);
const lambda2 = 1 - c / s;
let lambda = Math.sqrt(lambda2);
let it1: Vector3, it2: Vector3;
const ih = normalize3(cross(ir1, ir2));
if(ih.y < 0) { // transfer angle is larger than 180° (viewed from upper y axis)
lambda = -lambda;
it1 = cross(ir1, ih);
it2 = cross(ir2, ih);
} else {
it1 = cross(ih, ir1);
it2 = cross(ih, ir2);
}
it1 = normalize3(it1);
it2 = normalize3(it2);
const lambda3 = lambda * lambda2;
const T = Math.sqrt(2 * mu / (s*s*s)) * tof;
// There are 0 revolutions, we can directly calculate the unique x
const T0 = Math.acos(lambda) + lambda * Math.sqrt(1 - lambda2);
const T1 = 2 / 3 * (1 - lambda3);
// Initial guess
let x0: number;
if(T >= T0) {
x0 = -(T - T0) / (T - T0 + 4);
} else if(T <= T1) {
x0 = T1 * (T1 - T) / (2 / 5 * (1 - lambda2 * lambda3) * T) + 1;
} else {
x0 = Math.pow(T / T0, 0.69314718055994529 / Math.log(T1 / T0)) - 1;
}
// Householder iterations
const x = householderIterations(T, x0, 1e-15, lambda, 15);
// Reconstruct the terminal velocities from the calculated x
const gamma = Math.sqrt(mu * s / 2.0);
const rho = (r1 - r2) / c;
const sigma = Math.sqrt(1 - rho * rho);
const y = Math.sqrt(1.0 - lambda2 + lambda2 * x * x);
const vr1 = gamma * ((lambda * y - x) - rho * (lambda * y + x)) / r1;
const vr2 = -gamma * ((lambda * y - x) + rho * (lambda * y + x)) / r2;
const vt = gamma * sigma * (y + lambda * x);
const vt1 = vt / r1;
const vt2 = vt / r2;
const v1 = add3(mult3(ir1, vr1), mult3(it1, vt1));
const v2 = add3(mult3(ir2, vr2), mult3(it2, vt2));
return {v1, v2};
}
function householderIterations(T: number, x0: number, eps: number, lambda: number, maxIters: number) {
let err = 1;
let xnew = 0;
let tof = 0;
let delta = 0;
let DT = 0, DDT = 0, DDDT = 0;
for(let it = 0; err > eps && it < maxIters; it++) {
tof = x2tof(x0, lambda);
const lambda2 = 1 - c / s;
let lambda = Math.sqrt(lambda2);
const DTs = dTdx(DT, DDT, DDDT, x0, tof, lambda);
DT = DTs.DT;
DDT = DTs.DDT;
DDDT = DTs.DDDT;
let it1: Vector3, it2: Vector3;
delta = tof - T;
const DT2 = DT * DT;
xnew = x0 - delta * (DT2 - delta * DDT / 2) / (DT * (DT2 - delta * DDT) + DDDT * delta * delta / 6);
err = Math.abs(x0 - xnew);
x0 = xnew;
}
return x0;
}
function dTdx(DT: number, DDT: number, DDDT: number, x: number, T: number, lambda: number) {
const l2 = lambda * lambda;
const l3 = l2 * lambda;
const umx2 = 1.0 - x * x;
const y = Math.sqrt(1.0 - l2 * umx2);
const y2 = y * y;
const y3 = y2 * y;
DT = 1.0 / umx2 * (3.0 * T * x - 2.0 + 2.0 * l3 * x / y);
DDT = 1.0 / umx2 * (3.0 * T + 5.0 * x * DT + 2.0 * (1.0 - l2) * l3 / y3);
DDDT = 1.0 / umx2 * (7.0 * x * DDT + 8.0 * DT - 6.0 * (1.0 - l2) * l2 * l3 * x / y3 / y2);
return {DT, DDT, DDDT};
}
function x2tof2(x: number, lambda: number)
{
const a = 1.0 / (1.0 - x * x);
if (a > 0) // ellipse
{
let alfa = 2.0 * Math.acos(x);
let beta = 2.0 * Math.asin(Math.sqrt(lambda * lambda / a));
if (lambda < 0.0) beta = -beta;
return ((a * Math.sqrt(a) * ((alfa - Math.sin(alfa)) - (beta - Math.sin(beta)))) / 2.0);
} else {
let alfa = 2.0 * Math.acosh(x);
let beta = 2.0 *Math.asinh(Math.sqrt(-lambda * lambda / a));
if (lambda < 0.0) beta = -beta;
return (-a * Math.sqrt(-a) * ((beta - Math.sinh(beta)) - (alfa - Math.sinh(alfa))) / 2.0);
}
}
function x2tof(x: number, lambda: number)
{
const battin = 0.01;
const lagrange = 0.2;
const dist = Math.abs(x - 1);
if (dist < lagrange && dist > battin) { // We use Lagrange tof expression
return x2tof2(x, lambda);
}
const K = lambda * lambda;
const E = x * x - 1.0;
const rho = Math.abs(E);
const z = Math.sqrt(1 + K * E);
if (dist < battin) { // We use Battin series tof expression
const eta = z - lambda * x;
const S1 = 0.5 * (1.0 - lambda - x * eta);
let Q = hypergeometricF(S1, 1e-11);
Q = 4.0 / 3.0 * Q;
return (eta * eta * eta * Q + 4 * lambda * eta) / 2.0;
} else { // We use Lancaster tof expresion
const y = Math.sqrt(rho);
const g = x * z - lambda * E;
let d = 0.0;
if (E < 0) {
const l = Math.acos(g);
d = l;
if(ih.y < 0) { // transfer angle is larger than 180° (viewed from upper y axis)
lambda = -lambda;
it1 = cross(ir1, ih);
it2 = cross(ir2, ih);
} else {
const f = y * (z - lambda * x);
d = Math.log(f + g);
it1 = cross(ih, ir1);
it2 = cross(ih, ir2);
}
return (x - lambda * z - d / y) / E;
}
}
it1 = normalize3(it1);
it2 = normalize3(it2);
function hypergeometricF(z: number, tol: number)
{
let Sj = 1.0;
let Cj = 1.0;
let err = 1.0;
let Cj1 = 0.0;
let Sj1 = 0.0;
let j = 0;
while (err > tol) {
Cj1 = Cj * (3.0 + j) * (1.0 + j) / (2.5 + j) * z / (j + 1);
Sj1 = Sj + Cj1;
err = Math.abs(Cj1);
Sj = Sj1;
Cj = Cj1;
j++;
const lambda3 = lambda * lambda2;
const T = Math.sqrt(2 * mu / (s*s*s)) * tof;
// There are 0 revolutions, we can directly calculate the unique x
const T0 = Math.acos(lambda) + lambda * Math.sqrt(1 - lambda2);
const T1 = 2 / 3 * (1 - lambda3);
// Initial guess
let x0: number;
if(T >= T0) {
x0 = -(T - T0) / (T - T0 + 4);
} else if(T <= T1) {
x0 = T1 * (T1 - T) / (2 / 5 * (1 - lambda2 * lambda3) * T) + 1;
} else {
x0 = Math.pow(T / T0, 0.69314718055994529 / Math.log(T1 / T0)) - 1;
}
// Householder iterations
const x = householderIterations(T, x0, 1e-15, lambda, 15);
// Reconstruct the terminal velocities from the calculated x
const gamma = Math.sqrt(mu * s / 2.0);
const rho = (r1 - r2) / c;
const sigma = Math.sqrt(1 - rho * rho);
const y = Math.sqrt(1.0 - lambda2 + lambda2 * x * x);
const vr1 = gamma * ((lambda * y - x) - rho * (lambda * y + x)) / r1;
const vr2 = -gamma * ((lambda * y - x) + rho * (lambda * y + x)) / r2;
const vt = gamma * sigma * (y + lambda * x);
const vt1 = vt / r1;
const vt2 = vt / r2;
const v1 = add3(mult3(ir1, vr1), mult3(it1, vt1));
const v2 = add3(mult3(ir2, vr2), mult3(it2, vt2));
return {v1, v2};
}
return Sj;
}
function householderIterations(T: number, x0: number, eps: number, lambda: number, maxIters: number) {
let err = 1;
let xnew = 0;
let tof = 0;
let delta = 0;
let DT = 0, DDT = 0, DDDT = 0;
for(let it = 0; err > eps && it < maxIters; it++) {
tof = x2tof(x0, lambda);
const DTs = dTdx(DT, DDT, DDDT, x0, tof, lambda);
DT = DTs.DT;
DDT = DTs.DDT;
DDDT = DTs.DDDT;
delta = tof - T;
const DT2 = DT * DT;
xnew = x0 - delta * (DT2 - delta * DDT / 2) / (DT * (DT2 - delta * DDT) + DDDT * delta * delta / 6);
err = Math.abs(x0 - xnew);
x0 = xnew;
}
return x0;
}
function dTdx(DT: number, DDT: number, DDDT: number, x: number, T: number, lambda: number) {
const l2 = lambda * lambda;
const l3 = l2 * lambda;
const umx2 = 1.0 - x * x;
const y = Math.sqrt(1.0 - l2 * umx2);
const y2 = y * y;
const y3 = y2 * y;
DT = 1.0 / umx2 * (3.0 * T * x - 2.0 + 2.0 * l3 * x / y);
DDT = 1.0 / umx2 * (3.0 * T + 5.0 * x * DT + 2.0 * (1.0 - l2) * l3 / y3);
DDDT = 1.0 / umx2 * (7.0 * x * DDT + 8.0 * DT - 6.0 * (1.0 - l2) * l2 * l3 * x / y3 / y2);
return {DT, DDT, DDDT};
}
function x2tof2(x: number, lambda: number)
{
const a = 1.0 / (1.0 - x * x);
if (a > 0) // ellipse
{
let alfa = 2.0 * Math.acos(x);
let beta = 2.0 * Math.asin(Math.sqrt(lambda * lambda / a));
if (lambda < 0.0) beta = -beta;
return ((a * Math.sqrt(a) * ((alfa - Math.sin(alfa)) - (beta - Math.sin(beta)))) / 2.0);
} else {
let alfa = 2.0 * Math.acosh(x);
let beta = 2.0 * Math.asinh(Math.sqrt(-lambda * lambda / a));
if (lambda < 0.0) beta = -beta;
return (-a * Math.sqrt(-a) * ((beta - Math.sinh(beta)) - (alfa - Math.sinh(alfa))) / 2.0);
}
}
function x2tof(x: number, lambda: number)
{
const battin = 0.01;
const lagrange = 0.2;
const dist = Math.abs(x - 1);
if (dist < lagrange && dist > battin) { // We use Lagrange tof expression
return x2tof2(x, lambda);
}
const K = lambda * lambda;
const E = x * x - 1.0;
const rho = Math.abs(E);
const z = Math.sqrt(1 + K * E);
if (dist < battin) { // We use Battin series tof expression
const eta = z - lambda * x;
const S1 = 0.5 * (1.0 - lambda - x * eta);
let Q = hypergeometricF(S1, 1e-11);
Q = 4.0 / 3.0 * Q;
return (eta * eta * eta * Q + 4 * lambda * eta) / 2.0;
} else { // We use Lancaster tof expresion
const y = Math.sqrt(rho);
const g = x * z - lambda * E;
let d = 0.0;
if (E < 0) {
const l = Math.acos(g);
d = l;
} else {
const f = y * (z - lambda * x);
d = Math.log(f + g);
}
return (x - lambda * z - d / y) / E;
}
}
function hypergeometricF(z: number, tol: number)
{
let Sj = 1.0;
let Cj = 1.0;
let err = 1.0;
let Cj1 = 0.0;
let Sj1 = 0.0;
let j = 0;
while (err > tol) {
Cj1 = Cj * (3.0 + j) * (1.0 + j) / (2.5 + j) * z / (j + 1);
Sj1 = Sj + Cj1;
err = Math.abs(Cj1);
Sj = Sj1;
Cj = Cj1;
j++;
}
return Sj;
}
}

26
src/dedicated-workers/libs/math.ts
View File

@@ -1,6 +1,7 @@
/* Useful little functions and constants */
const TWO_PI = 2 * Math.PI;
const HALF_PI = 0.5 * Math.PI;
function clamp(x: number, min: number, max: number) {
return x > max ? max : x < min ? min : x;
@@ -28,31 +29,6 @@ function newtonRootSolve(
return x;
}
function randomPointOnSphereRing(thetaMin: number, thetaMax: number){
return {
theta: randomInInterval(thetaMin, thetaMax),
phi: randomInInterval(0, TWO_PI)
}
}
function pointOnSphereRing(dir: Vector3, theta: number, phi: number){
let axis = normalize3(cross(dir, vec3(1, 0, 0)));
let point = rotate3(dir, axis, theta);
point = rotate3(point, dir, phi);
return point;
}
function minMaxRingAngle(r: number, rmin: number, rmax: number){
return {
thetaMin: Math.asin(rmin/r),
thetaMax: Math.asin(rmax/r)
};
}
function randomInInterval(a: number, b: number) {
return a + Math.random() * (b - a);
}
function randint(a: number, b: number){
return Math.floor(a + Math.random() * (b - a + 1));
}

8
src/dedicated-workers/libs/physics-2d.ts
View File

@@ -35,14 +35,14 @@ namespace Physics2D
const p = a * (1 - e*e);
// Periapsis radius and direction
const rp = a*(1 - e);
const pdir = mult2(div2(evec, e), rp);
const pvec = mult2(div2(evec, e), rp);
// Check if orbit is clockwise or counter-clockwise
const clk = det(pos, vel) < 0;
return {
eccentricity: e,
periapsisDir: pdir,
periapsisVec: pvec,
semiMajorAxis: a,
orbitalParam: p,
clockwise: clk
@@ -155,7 +155,7 @@ namespace Physics2D
const orbitElts = stateToOrbitElements(attractor, state);
const e = orbitElts.eccentricity;
const p = orbitElts.orbitalParam;
const pdir = orbitElts.periapsisDir;
const pvec = orbitElts.periapsisVec;
const pos0 = state.pos;
const vel0 = state.vel;
@@ -189,7 +189,7 @@ namespace Physics2D
y: (e + cosNu) * vmag * vdir
};
const pangle = Math.atan2(pdir.y,pdir.x);
const pangle = Math.atan2(pvec.y,pvec.x);
return {
pos: rotate2(pos, pangle),

936
src/dedicated-workers/libs/physics-3d.ts
View File

@@ -1,446 +1,526 @@
function orbitalElementsFromOrbitData(orbit: IOrbit) : OrbitalElements{
// Calculate main direction vectors
// with reference direction (1, 0, 0)
const right = vec3(1, 0, 0), up = vec3(0, 1, 0);
const ascNodeDir = rotate3(right, up, orbit.ascNodeLongitude);
const normal = rotate3(up, ascNodeDir, orbit.inclination);
const periapsisDir = rotate3(ascNodeDir, normal, orbit.argOfPeriapsis);
return {
semiMajorAxis: orbit.semiMajorAxis,
eccentricity: orbit.eccentricity,
periapsiDir: periapsisDir,
inclination: orbit.inclination,
argOfPeriapsis: orbit.argOfPeriapsis,
ascNodeLongitude: orbit.ascNodeLongitude,
ascNodeDir: ascNodeDir,
// @ts-ignore
orbitalParam: orbit.orbitalParam
};
}
function equatorialCircularOrbit(radius: number) : OrbitalElements{
return {
eccentricity: 0,
periapsisDir: vec3(1, 0, 0),
semiMajorAxis: radius,
inclination: 0,
argOfPeriapsis: 0,
ascNodeLongitude: 0,
ascNodeDir: vec3(1, 0, 0),
orbitalParam: radius
};
}
/**
* Calculates the state vectors of an orbiting body at a given date
* @param body The body to calculate the state of
* @param date The date (in s)
* @param mu The standard gravitational parameter of the the attractor body around which `body` is orbiting
* @returns The state vectors of the body at the specified date
*/
function getBodyStateAtDate(body: IOrbitingBody, date: number, mu: number) {
const {orbit} = body;
const n = <number>orbit.meanMotion;
const M = body.meanAnomaly0 + n * date;
const nu = solveTrueAnomalyFromMeanAnomaly(orbit.eccentricity, M);
const orbitElts = orbitalElementsFromOrbitData(orbit);
return stateFromOrbitElements(orbitElts, mu, nu);
}
/**
* @param orbit An orbit
* @param soi The sphere of influence radius of the orbit's attractor body
* @returns The true anomaly of the point on that orbit leaving the SOI.
*/
function soiExitTrueAnomaly(orbit: OrbitalElements, soi: number){
const p = orbit.orbitalParam;
const e = orbit.eccentricity;
return Math.acos((p/soi - 1) / e);
}
/**
* @param orbit An orbit
* @returns The radius of the periapsis of that orbit.
*/
function periapsisRadius(orbit: OrbitalElements){
return orbit.semiMajorAxis * (1 - orbit.eccentricity);
}
/**
* Calculates the time of flight between two anomalies on a given orbit.
* @param orbit The orbit
* @param attractor The orbit's attractor
* @param nu1 The start true anomaly
* @param nu2 The end true anomaly
* @returns The time of flight from the first to the second anomaly provided.
*/
function tofBetweenAnomalies(orbit: OrbitalElements, attractor: ICelestialBody, nu1: number, nu2: number){
const mu = attractor.stdGravParam;
const a = orbit.semiMajorAxis;
const e = orbit.eccentricity;
/*let EH1 = eccentricAnomalyFromTrueAnomaly(nu1, e);
let EH2 = eccentricAnomalyFromTrueAnomaly(nu2, e);*/
let tof = 0;
if(e < 1) {
let E1 = eccentricAnomalyFromTrueAnomaly(nu1, e);
let E2 = eccentricAnomalyFromTrueAnomaly(nu2, e);
const invN = Math.sqrt(a*a*a/mu);
if(E2 < E1) {
E2 += TWO_PI;
}
tof = invN * ( E2 - E1 + e*(Math.sin(E1) - Math.sin(E2)) );
} else {
let H1 = eccentricAnomalyFromTrueAnomaly(nu1, e);
let H2 = eccentricAnomalyFromTrueAnomaly(nu2, e);
if(H2 < H1) {
const t = H2;
H2 = H1;
H1 = t;
}
const invN = Math.sqrt(-a*a*a/mu);
tof = -invN * ( H2 - H1 + e*(Math.sinh(H1) - Math.sinh(H2)) );
}
/*if(tof < 0){
debug("ERROR");
}*/
return tof;
}
/**
*
* @param attractor The attractor body
* @param radius The orbit radius
* @returns The velocity for a circular orbit around the specified body at the specified radius.
*/
function circularVelocity(attractor: ICelestialBody, radius: number) {
return Math.sqrt(attractor.stdGravParam / radius);
}
/**
* @param attractor The attractor body
* @param radius The departure radius
* @returns The minimum velocity required to reach escape the body's gravtational attraction.
*/
function ejectionVelocity(attractor: ICelestialBody, radius: number){
return Math.SQRT2 * circularVelocity(attractor, radius);
}
/**
* Calculates the magnitude of the velocity for a given radius on an orbit.
* @param orbit The orbit
* @param attractor The attractor body of the orbit
* @param radius The radius from the attractor (must exist on that orbit)
* @returns The velocity at the specified radius of the specified orbit
*/
function velocityAtRadius(orbit: OrbitalElements, attractor: ICelestialBody, radius: number){
return Math.sqrt(attractor.stdGravParam * (2/radius - 1/orbit.semiMajorAxis));
}
/**
* @param bodyPos The position the body relative to its attractor
* @param retrograde Whether the direction should be prograde or retrograde to the body's movement
* @returns The tangent vector pointing prograde or retrograde in the body's velocity direction
* and the angle of departure from the reference direction
*/
function idealEjectionDirection(bodyPos: Vector3, retrograde: boolean) {
const tangent = normalize3(vec3(bodyPos.z, 0, -bodyPos.x));
let nu = Math.acos(bodyPos.x / mag3(bodyPos));
nu = bodyPos.z >= 0 ? nu : TWO_PI - nu;
if(retrograde) {
tangent.x *= -1;
tangent.z *= -1;
nu = (nu + Math.PI) % TWO_PI;
}
return {tangent, nu};
}
/**
* Calculates the offset angle of the ejection manoeuvre to perform to get an ejection velocity parallel to
* The body's movement direction
* @param perigeeVel The velocity at the periapsis of the orbit
* @param perigeeRadius The radius of the periapsis
* @param attractor The attractor body of the orbit
* @returns The angle needed to have an asymptotic ejection velocity parallel to the body's direction
*/
function hyperbolicEjectionOffsetAngle(perigeeVel: number, perigeeRadius: number, attractor: ICelestialBody){
const e = perigeeRadius*perigeeVel*perigeeVel/attractor.stdGravParam - 1;
return Math.PI - Math.acos(-1/e);
}
/**
* Calculates the orbital elements of an orbit from cartesian state vectors, based on:
* https://downloads.rene-schwarz.com/download/M002-Cartesian_State_Vectors_to_Keplerian_Orbit_Elements.pdf
* @param state The orbital state containing the cartesian state vectors, relative to the attractor
* @param attractor The attractor body of the orbit
* @returns The corresponding orbital elements.
*/
function stateToOrbitElements(state: {pos: Vector3, vel: Vector3}, attractor: ICelestialBody) : OrbitalElements {
const mu = attractor.stdGravParam;
const pos = state.pos;
const vel = state.vel;
const r = mag3(pos);
const v2 = magSq3(vel);
namespace Physics3D
{
/**
* Converts orbit data to the corresponding orbital elements
* @param orbit The orbit data
* @returns The corresponding orbital elements
*/
export function orbitElementsFromOrbitData(orbit: IOrbit) : OrbitalElements3D {
// Calculate main direction vectors
// with reference direction (1, 0, 0)
const right = vec3(1, 0, 0), up = vec3(0, 1, 0);
const ascNodeDir = rotate3(right, up, orbit.ascNodeLongitude);
const normal = rotate3(up, ascNodeDir, orbit.inclination);
const periapsisDir = rotate3(ascNodeDir, normal, orbit.argOfPeriapsis);
// Momentum
const h = cross(pos, vel);
// Eccentricity vector
let evec = sub3(div3(cross(vel, h), mu), div3(pos, r));
let e = mag3(evec);
if(e <= 1e-15) {
evec = vec3(0, 0, 0);
e = 0;
return {
semiMajorAxis: orbit.semiMajorAxis,
eccentricity: orbit.eccentricity,
periapsiDir: periapsisDir,
inclination: orbit.inclination,
argOfPeriapsis: orbit.argOfPeriapsis,
ascNodeLongitude: orbit.ascNodeLongitude,
ascNodeDir: ascNodeDir,
// @ts-ignore
orbitalParam: orbit.orbitalParam
};
}
// Vector pointing towards the ascending node
let nvec = cross(vec3(0, 1, 0), h);
const n = mag3(nvec);
// Orbital inclination
const i = Math.acos(h.y / mag3(h));
// Longitude of the ascending node
const t_lan = Math.acos(nvec.x / n);
let lan = nvec.z <= 0 ? t_lan : TWO_PI - t_lan;
// Argument of the periapsis
const t_arg = Math.acos(dot3(nvec, evec) / (e * n));
let arg = evec.y >= 0 ? t_arg : TWO_PI - t_arg;
// Semi major axis
const a = 1 / ( 2/r - v2/mu );
// Calculate the orbital parameter
const p = a * (1 - e*e);
// Use alternate orbital elements in special cases
if(i == 0 && e != 0) {
// Longitude of perigee
const t_plong = Math.acos(dot3(vec3(1, 0, 0), evec) / e);
const plong = evec.z < 0 ? t_plong : TWO_PI - t_plong;
lan = 0;
nvec = vec3(1, 0, 0);
arg = plong;
} else if(i != 0 && e == 0) {
arg = 0;
evec = clone3(nvec);
} else if(i == 0 && e == 0) {
lan = 0;
nvec = vec3(1, 0, 0);
arg = 0;
evec = vec3(1, 0, 0);
/**
* Returns the orbital elements corresponding to an equatorial circular orbit
* with the specified radius.
* @param radius The radius of the orbit
* @returns The orbital elements
*/
export function equatorialCircularOrbit(radius: number) : OrbitalElements3D {
return {
eccentricity: 0,
periapsisDir: vec3(1, 0, 0),
semiMajorAxis: radius,
inclination: 0,
argOfPeriapsis: 0,
ascNodeLongitude: 0,
ascNodeDir: vec3(1, 0, 0),
orbitalParam: radius
};
}
return {
eccentricity: e,
periapsisDir: normalize3(evec),
semiMajorAxis: a,
inclination: i,
argOfPeriapsis: arg,
ascNodeLongitude: lan,
ascNodeDir: normalize3(nvec),
orbitalParam: p
};
}
/**
* Calculates the state vectors of an orbiting body at a given date
* @param body The body to calculate the state of
* @param orbit The orbital elements of the body (built from `body.orbit`)
* @param attractor The attractor around which the body is orbiting
* @param date The date (in s)
* @returns The state vectors of the body at the specified date
*/
export function bodyStateAtDate(body: IOrbitingBody, orbit: OrbitalElements3D, attractor: ICelestialBody, date: number) {
const n = body.orbit.meanMotion as number;
const M = body.meanAnomaly0 + n * date;
const nu = trueAnomalyFromMeanAnomaly(orbit.eccentricity, M);
return orbitElementsToState(orbit, attractor, nu);
}
/**
* Calculates the true anomaly of state on an orbit, considering extreme cases with alternate
* orbital elements described here:
* https://www.faa.gov/about/office_org/headquarters_offices/avs/offices/aam/cami/library/online_libraries/aerospace_medicine/tutorial/media/III.4.1.4_Describing_Orbits.pdf
* @param orbitElts The orbital elements describing the orbit
* @param state A state on that orbit
* @returns The true anomaly of that state on the specified orbit
*/
function trueAnomalyFromOrbitalState(orbitElts: OrbitalElements, state: {pos: Vector3, vel: Vector3}) {
const pos = state.pos;
const vel = state.vel;
const r = mag3(pos);
/**
* Returns the velocity on a circular orbit of radius `radius` around the specified
* attractor.
* @param attractor The attractor body
* @param radius The orbit radius
* @returns The velocity for a circular orbit around the specified body at the specified radius.
*/
export function circularVelocity(attractor: ICelestialBody, radius: number) {
return Math.sqrt(attractor.stdGravParam / radius);
}
const e = orbitElts.eccentricity;
const i = orbitElts.inclination;
/**
* Returns the minimum velocity required to reach a hypberolic ejection orbit from a circular
* orbit with radius `radius` around the specified attractor.
* @param attractor The attractor body
* @param radius The departure radius
* @returns The minimum velocity required to reach escape the body's gravtational attraction.
*/
export function ejectionVelocity(attractor: ICelestialBody, radius: number){
return Math.SQRT2 * circularVelocity(attractor, radius);
}
let nu: number;
/**
* @param orbit An orbit
* @returns The radius of the periapsis of that orbit.
*/
export function periapsisRadius(orbit: OrbitalElements3D){
return orbit.semiMajorAxis * (1 - orbit.eccentricity);
}
if(e != 0) {
// Elliptic orbit : we directly compute the true anomaly
const evec = orbitElts.periapsisDir;
const t_nu = Math.acos(dot3(evec, pos) / r);
const d = dot3(pos, vel);
/**
* Returns the positive true anomaly of the point in the orbit reaching the
* specified radius. Throws an error if the orbit never reaches the specified radius.
* @param orbit The orbital elements of the orbit
* @param radius The radius at which calculate the true anomaly
* @returns The positive true anomaly
*/
export function trueAnomalyAtRadius(orbit: OrbitalElements3D, radius: number){
const p = orbit.orbitalParam;
const e = orbit.eccentricity;
const cosNu = (p/radius - 1) / e;
if(Math.abs(cosNu) > 1)
throw new Error("The orbit never reaches that radius.");
return Math.acos(cosNu);
}
/**
* Returns the velocity at the specified radius on the specified orbit.
* @param orbit The orbit
* @param attractor The attracor body of the orbit
* @param radius The radius at which compute the velocity
* @returns The velocity at the specified radius
*/
export function velocityAtRadius(orbit: OrbitalElements3D, attractor: ICelestialBody, radius: number){
const mu = attractor.stdGravParam;
const a = orbit.semiMajorAxis;
const v2 = mu*(2/radius - 1/a);
if(v2 < 0)
throw new Error("Invalid radius.");
return Math.sqrt(v2);
}
/**
* Computes the required velocity that a vessel must have at radius `r0` around the specified
* body to reach the altitude `r`, considering a Hohmann like transfer. If the computed velocity
* is greater than the ejection velocity at radius `r0`, the ejection velocity returned.
* @param body The orbited body
* @param r0 The radius from which we calculate the velocity
* @param r The altitude we want to reach
* @returns The required velocity
*/
export function velocityToReachAltitude(body: ICelestialBody, r0: number, r: number){
const mu = body.stdGravParam;
const a = 0.5 * (r0 + r);
const v2 = mu*(2/r0 - 1/a);
if(v2 < 0)
throw new Error("Invalid radius.");
const vej2 = 2 * mu / r0;
return Math.sqrt(Math.min(vej2, v2));
}
/**
* Returns the orbital period of the orbit with the given attractor and semi major axis.
* @param attractor The attractor body of the orbit
* @param a The semi major axis of the orbit
* @returns The orbital period
*/
export function orbitPeriod(attractor: ICelestialBody, a: number){
return TWO_PI * Math.sqrt(a*a*a/attractor.stdGravParam);
}
/**
* Calculates the specific energy of an orbit from the radius and velocity magnitude
* at one of its point.
* @param attractor The attractor body of the orbit
* @param r A radius reached by the orbit
* @param v The velocity at the specified radius
* @returns The specific energy
*/
export function specificEnergy(attractor: ICelestialBody, r: number, v: number){
return 0.5*v*v - attractor.stdGravParam/r;
}
/**
* Computes the magnitude of the velocity at some radius orbit, knowing the radius and velocity
* at some other point on the orbit.
* @param attractor The attractor body of the orbit
* @param r0 A known radius reached by the orbit
* @param v0 The known velocity at the known radius
* @param r A (reached) radius at which we want to know the velocity magnitude
* @returns The velocity at radius `r`
*/
export function deduceVelocityAtRadius(attractor: ICelestialBody, r0: number, v0: number, r: number){
const se = specificEnergy(attractor, r0, v0);
return Math.sqrt(2*(se + attractor.stdGravParam/r));
}
/**
* Calculates the orbital elements of an orbit from cartesian state vectors, based on:
* https://downloads.rene-schwarz.com/download/M002-Cartesian_State_Vectors_to_Keplerian_Orbit_Elements.pdf
* @param state The orbital state containing the cartesian state vectors, relative to the attractor
* @param attractor The attractor body of the orbit
* @returns The corresponding orbital elements.
*/
export function stateToOrbitElements(state: OrbitalState3D, attractor: ICelestialBody) : OrbitalElements3D {
const mu = attractor.stdGravParam;
const pos = state.pos;
const vel = state.vel;
const r = mag3(pos);
const v2 = magSq3(vel);
const nullEps = 1e-10;
// Momentum
const h = cross(pos, vel);
// Eccentricity vector
let evec = sub3(div3(cross(vel, h), mu), div3(pos, r));
let e = mag3(evec);
if(e <= nullEps) {
evec = vec3(0, 0, 0);
e = 0;
}
// Vector pointing towards the ascending node
let nvec = cross(vec3(0, 1, 0), h);
const n = mag3(nvec);
// Orbital inclination
let i = Math.acos(h.y / mag3(h));
let inXZPlane = false;
if(Math.abs(i) < nullEps) {
i = 0;
inXZPlane = true;
} else if(Math.abs(i - Math.PI) < nullEps) {
i = Math.PI;
inXZPlane = true;
} else if(Math.abs(i + Math.PI) < nullEps) {
i = -Math.PI;
inXZPlane = true;
}
const cosEps = 1e-5;
// Longitude of the ascending node
let cos_lan = nvec.x / n;
if(Math.abs(cos_lan) < 1 + cosEps) {
cos_lan = clamp(cos_lan, -1, 1);
}
const t_lan = Math.acos(cos_lan);
let lan = nvec.z <= 0 ? t_lan : TWO_PI - t_lan;
// Argument of the periapsis
let cos_arg = dot3(nvec, evec) / (e * n);
if(Math.abs(cos_arg) < 1 + cosEps){
cos_arg = clamp(cos_arg, -1, 1);
}
const t_arg = Math.acos(cos_arg);
let arg = evec.y >= 0 ? t_arg : TWO_PI - t_arg;
// Semi major axis
const a = 1 / ( 2/r - v2/mu );
// Calculate the orbital parameter
const p = a * (1 - e*e);
// Use alternate orbital elements in special cases
if(inXZPlane && e != 0) {
// Longitude of perigee
const t_plong = Math.acos(dot3(vec3(1, 0, 0), evec) / e);
const plong = evec.z < 0 ? t_plong : TWO_PI - t_plong;
lan = 0;
nvec = vec3(1, 0, 0);
arg = plong;
} else if(!inXZPlane && e == 0) {
arg = 0;
evec = clone3(nvec);
} else if(inXZPlane && e == 0) {
lan = 0;
nvec = vec3(1, 0, 0);
arg = 0;
evec = vec3(1, 0, 0);
}
return {
eccentricity: e,
periapsisDir: normalize3(evec),
semiMajorAxis: a,
inclination: i,
argOfPeriapsis: arg,
ascNodeLongitude: lan,
ascNodeDir: normalize3(nvec),
orbitalParam: p
};
}
/**
* Calculates the state vectors at the specified true anomaly for a given orbit.
* @param orbit The orbital elements describing the orbit
* @param attractor The attractor body
* @param nu The true anomaly at which get the state
* @returns The state vectors
*/
export function orbitElementsToState(orbit: OrbitalElements3D, attractor: ICelestialBody, nu: number) : OrbitalState3D {
const mu = attractor.stdGravParam;
const a = orbit.semiMajorAxis;
const e = orbit.eccentricity;
const p = orbit.orbitalParam;
nu *= -1; // Without this the resulting position is on the opposite side
const r = p / (1 + e * Math.cos(nu)); // radius of the orbit
// 2D position with periapsis pointing towards (1, 0, 0) (reference direction)
let pos = vec3(r * Math.cos(nu), 0, r * Math.sin(nu));
// the corresponding 2D velocity
let vel: Vector3;
if(e < 1) { // Elliptical orbit
const v = Math.sqrt(mu * a) / r;
const E = eccentricAnomalyFromTrueAnomaly(nu, e);
vel = vec3(v * Math.sin(E), 0, -v * Math.sqrt(1 - e*e) * Math.cos(E));
} else { // Hyperbolic orbit (the case e = 1, parabolic orbit, will never be reached)
const v = Math.sqrt(-mu * a) / r;
const H = eccentricAnomalyFromTrueAnomaly(nu, e);
vel = vec3(v * Math.sinh(H), 0, -v * Math.sqrt(e*e - 1) * Math.cosh(H));
}
// Rotating the plane vectors towards their real 3D position
const right = vec3(1, 0, 0), up = vec3(0, 1, 0);
const ascNodeDir = rotate3(right, up, orbit.ascNodeLongitude);
pos = rotate3(pos, up, orbit.ascNodeLongitude);
pos = rotate3(pos, up, orbit.argOfPeriapsis);
pos = rotate3(pos, ascNodeDir, orbit.inclination);
vel = rotate3(vel, up, orbit.ascNodeLongitude);
vel = rotate3(vel, up, orbit.argOfPeriapsis);
vel = rotate3(vel, ascNodeDir, orbit.inclination);
return {pos, vel};
}
/**
* Calculates the true anomaly of state on an orbit, considering extreme cases with alternate
* orbital elements described here:
* https://www.faa.gov/about/office_org/headquarters_offices/avs/offices/aam/cami/library/online_libraries/aerospace_medicine/tutorial/media/III.4.1.4_Describing_Orbits.pdf
* @param orbit The orbital elements describing the orbit
* @param state A state on that orbit
* @returns The true anomaly of that state on the specified orbit
*/
export function trueAnomalyFromOrbitalState(orbit: OrbitalElements3D, state: OrbitalState3D) {
const pos = state.pos;
const vel = state.vel;
const r = mag3(pos);
const e = orbit.eccentricity;
const i = orbit.inclination;
let nu: number;
if(e != 0) {
// Elliptic orbit : we directly compute the true anomaly
const evec = orbit.periapsisDir;
const t_nu = Math.acos(dot3(evec, pos) / r);
const d = dot3(pos, vel);
if(e < 1) {
nu = d >= 0 ? t_nu : TWO_PI - t_nu;
} else {
nu = d >= 0 ? t_nu : -t_nu;
}
} else if(i != 0) { // && e == 0
// Circular inclined orbit: the true anomaly is the argument of latitude
const nvec = orbit.ascNodeDir;
const t_u = Math.acos(dot3(nvec, pos) / r);
let u: number;
if(e < 1) {
u = pos.y >= 0 ? t_u : TWO_PI - t_u;
} else {
u = pos.y >= 0 ? t_u : -t_u;
}
nu = u;
} else { // i == 0 && e == 0
// Circular and equatorial orbit : the true anomaly is the true longitude
const t_l = Math.acos(pos.x / r);
let l: number;
if(e < 1) {
l = vel.x <= 0 ? t_l : TWO_PI - t_l;
} else {
l = vel.x <= 0 ? t_l : -t_l;
}
nu = l;
}
return nu;
}
/**
* Computes the eccentric (or hyperbolic) anomaly for a given eccentricity and true anomaly
* @param nu The true anomaly
* @param e The eccentricity
* @returns The eccentric anomaly
*/
export function eccentricAnomalyFromTrueAnomaly(nu: number, e: number) {
if(e < 1) { // Elliptic orbit
// Eccentric anomaly
return 2 * Math.atan( Math.tan(nu*0.5) * Math.sqrt((1-e)/(1+e)) );
} else { // Hyperbolic orbit, the case e = 1 (parabolic orbit) will never be reached
// Hyperbolic eccentric anomaly
return 2 * Math.atanh( Math.tan(nu*0.5) * Math.sqrt((e-1)/(e+1)) );
}
}
/**
* Computes the mean anomaly from the eccentric (or hyperbolic) anomaly for a given eccentricity
* @param EH The eccentric or hyperbolic anomaly
* @param e The eccentricity
* @returns The mean anomaly
*/
export function meanAnomalyFromEccentricAnomaly(EH: number, e: number) {
if(e < 1) { // Elliptic orbit
return EH - e * Math.sin(EH);
} else { // Hyperbolic orbit, the case e = 1 (parabolic orbit) will never be reached
return e * Math.sinh(EH) - EH;
}
}
/**
* Calculates the mean anomaly corresponding to the specified true anomaly of an orbit
* with the specified eccentricity.
* @param nu The true anomaly
* @param e The eccentricty of the orbit
* @returns The corresponding mean anomaly
*/
export function meanAnomalyFromTrueAnomaly(nu: number, e: number) {
if(e < 1) { // Elliptic orbit
// Eccentric anomaly
const E = 2 * Math.atan( Math.tan(nu*0.5) * Math.sqrt((1-e)/(1+e)) );
// Mean anomaly
const M = E - e * Math.sin(E);
return M;
} else { // Hyperbolic orbit, the case e = 1 (parabolic orbit) will never be reached
// Hyperbolic eccentric anomaly
const H = 2 * Math.atanh( Math.tan(nu*0.5) * Math.sqrt((e-1)/(e+1)) );
// Mean anomaly
const M = e * Math.sinh(H) - H;
return M;
}
}
/**
* Computes the true anomaly from a mean anomaly for an orbit with the specifed eccentricity.
* @param e The eccentric anomaly of the orbit
* @param M The mean anomaly at which compute the true anomaly
* @returns The true anomaly
*/
export function trueAnomalyFromMeanAnomaly(e: number, M: number){
// Solving Kepler's equation for eccentric anomaly with Newton's method.
if(e < 1) {
nu = d >= 0 ? t_nu : TWO_PI - t_nu;
const E = newtonRootSolve(
x => x - e * Math.sin(x) - M,
x => 1 - e * Math.cos(x),
M,
1e-15
);
return 2 * Math.atan(Math.sqrt((1 + e)/(1 - e)) * Math.tan(E * 0.5));
} else {
nu = d >= 0 ? t_nu : -t_nu;
const H = newtonRootSolve(
x => e * Math.sinh(x) - x - M,
x => e * Math.cosh(x) - 1,
M,
1e-15
);
return 2 * Math.atan(Math.sqrt((e + 1)/(e - 1)) * Math.tanh(H * 0.5));
}
}
} else if(i != 0) { // && e == 0
// Circular inclined orbit: the true anomaly is the argument of latitude
const nvec = orbitElts.ascNodeDir;
const t_u = Math.acos(dot3(nvec, pos) / r);
let u: number;
/**
* Calculates the orbital state reached after the specified time of flight from the
* starting point defined by the specified initial true anomaly.
* @param orbit The orbital elements describing the orbit
* @param attractor The attractor body of the orbit
* @param nu0 The initial true anomaly
* @param deltaTime The duration of the flight between the initial point
* @returns The orbital state vectors.
*/
export function propagateStateFromTrueAnomaly(orbit: OrbitalElements3D, attractor: ICelestialBody, nu0: number, deltaTime: number){
const M0 = meanAnomalyFromTrueAnomaly(nu0, orbit.eccentricity);
const a = Math.abs(orbit.semiMajorAxis);
const mu = attractor.stdGravParam;
const n = Math.sqrt(mu / (a*a*a));
const MNext = M0 + n * deltaTime;
const nuNext = trueAnomalyFromMeanAnomaly(orbit.eccentricity, MNext);
return orbitElementsToState(orbit, attractor, nuNext);
}
/**
* Calculates the time of flight between two true anomalies on a given orbit.
* @param orbit The orbit
* @param attractor The orbit's attractor
* @param nu1 The start true anomaly
* @param nu2 The end true anomaly
* @returns The time of flight from the start to the end true anomaly.
*/
export function tofBetweenAnomalies(orbit: OrbitalElements3D, attractor: ICelestialBody, nu1: number, nu2: number){
const mu = attractor.stdGravParam;
const a = orbit.semiMajorAxis;
const e = orbit.eccentricity;
let tof = 0;
if(e < 1) {
u = pos.y >= 0 ? t_u : TWO_PI - t_u;
let E1 = eccentricAnomalyFromTrueAnomaly(nu1, e);
let E2 = eccentricAnomalyFromTrueAnomaly(nu2, e);
const invN = Math.sqrt(a*a*a/mu);
if(E2 < E1) {
E2 += TWO_PI;
}
tof = invN * ( E2 - E1 + e*(Math.sin(E1) - Math.sin(E2)) );
} else {
u = pos.y >= 0 ? t_u : -t_u;
let H1 = eccentricAnomalyFromTrueAnomaly(nu1, e);
let H2 = eccentricAnomalyFromTrueAnomaly(nu2, e);
if(H2 < H1) {
const t = H2;
H2 = H1;
H1 = t;
}
const invN = Math.sqrt(-a*a*a/mu);
tof = -invN * ( H2 - H1 + e*(Math.sinh(H1) - Math.sinh(H2)) );
}
nu = u;
} else { // i == 0 && e == 0
// Circular and equatorial orbit : the true anomaly is the true longitude
const t_l = Math.acos(pos.x / r);
let l: number;
if(e < 1) {
l = vel.x <= 0 ? t_l : TWO_PI - t_l;
} else {
l = vel.x <= 0 ? t_l : -t_l;
}
nu = l;
if(tof < 0)
throw new Error("Negative TOF calculated.");
return tof;
}
return nu;
}
function eccentricAnomalyFromTrueAnomaly(nu: number, e: number) {
if(e < 1) { // Elliptic orbit
// Eccentric anomaly
return 2 * Math.atan( Math.tan(nu*0.5) * Math.sqrt((1-e)/(1+e)) );
} else { // Hyperbolic orbit, the case e = 1 (parabolic orbit) will never be reached
// Hyperbolic eccentric anomaly
return 2 * Math.atanh( Math.tan(nu*0.5) * Math.sqrt((e-1)/(e+1)) );
}
}
function meanAnomalyFromEccentricAnomaly(EH: number, e: number) {
if(e < 1) { // Elliptic orbit
return EH - e * Math.sin(EH);
} else { // Hyperbolic orbit, the case e = 1 (parabolic orbit) will never be reached
return e * Math.sinh(EH) - EH;
}
}
/**
* Calculates the mean anomaly corresponding to the specified true anomaly of an orbit
* with the specified eccentricity.
* @param nu The true anomaly
* @param e The eccentricty of the orbit
* @returns The corresponding mean anomaly
*/
function meanAnomalyFromTrueAnomaly(nu: number, e: number) {
if(e < 1) { // Elliptic orbit
// Eccentric anomaly
const E = 2 * Math.atan( Math.tan(nu*0.5) * Math.sqrt((1-e)/(1+e)) );
// Mean anomaly
const M = E - e * Math.sin(E);
return M;
} else { // Hyperbolic orbit, the case e = 1 (parabolic orbit) will never be reached
// Hyperbolic eccentric anomaly
const H = 2 * Math.atanh( Math.tan(nu*0.5) * Math.sqrt((e-1)/(e+1)) );
// Mean anomaly
const M = e * Math.sinh(H) - H;
return M;
}
}
function solveTrueAnomalyFromMeanAnomaly(e: number, M: number){
// Solving Kepler's equation for eccentric anomaly with Newton's method.
if(e < 1) {
const E = newtonRootSolve(
x => x - e * Math.sin(x) - M,
x => 1 - e * Math.cos(x),
M,
1e-15
);
return 2 * Math.atan(Math.sqrt((1 + e)/(1 - e)) * Math.tan(E * 0.5));
} else {
const H = newtonRootSolve(
x => e * Math.sinh(x) - x - M,
x => e * Math.cosh(x) - 1,
M,
1e-15
);
return 2 * Math.atan(Math.sqrt((e + 1)/(e - 1)) * Math.tanh(H * 0.5));
}
}
/**
* Calculates the state vectors at the specified true anomaly for a given orbit.
* @param orbit The orbital elements describing the orbit
* @param mu The attractor's body standard gravitational parameter
* @param nu The true anomaly at which get the state
* @returns The state vectors
*/
function stateFromOrbitElements(orbit: OrbitalElements, mu: number, nu: number){
const a = orbit.semiMajorAxis;
const e = orbit.eccentricity;
const p = orbit.orbitalParam;
nu *= -1; // Without this the resulting position is on the opposite side
const r = p / (1 + e * Math.cos(nu)); // radius of the orbit
// 2D position with periapsis pointing towards (1, 0, 0) (reference direction)
let pos = vec3(r * Math.cos(nu), 0, r * Math.sin(nu));
// the corresponding 2D velocity
let vel: Vector3;
if(e < 1) { // Elliptical orbit
const v = Math.sqrt(mu * a) / r;
const E = eccentricAnomalyFromTrueAnomaly(nu, e);
vel = vec3(v * Math.sin(E), 0, -v * Math.sqrt(1 - e*e) * Math.cos(E));
} else { // Hyperbolic orbit (the case e = 1, parabolic orbit, will never be reached)
const v = Math.sqrt(-mu * a) / r;
const H = eccentricAnomalyFromTrueAnomaly(nu, e);
vel = vec3(v * Math.sinh(H), 0, -v * Math.sqrt(e*e - 1) * Math.cosh(H));
}
// Rotating the plane vectors towards their real 3D position
const right = vec3(1, 0, 0), up = vec3(0, 1, 0);
const ascNodeDir = rotate3(right, up, orbit.ascNodeLongitude);
pos = rotate3(pos, up, orbit.ascNodeLongitude);
pos = rotate3(pos, up, orbit.argOfPeriapsis);
pos = rotate3(pos, ascNodeDir, orbit.inclination);
vel = rotate3(vel, up, orbit.ascNodeLongitude);
vel = rotate3(vel, up, orbit.argOfPeriapsis);
vel = rotate3(vel, ascNodeDir, orbit.inclination);
return {pos, vel};
}
/**
* Calculates the orbital state reached after the specified time of flight from the
* starting point defined by the specified initial true anomaly.
* @param orbit The orbital elements describing the orbit
* @param nu0 The initial true anomaly
* @param deltaTime The duration of the flight between the initial point
* @param mu The attractor body's standard gravitational parameter
* @returns The orbital state vectors.
*/
function propagateState(orbit: OrbitalElements, nu0: number, deltaTime: number, attractor: ICelestialBody){
const M0 = meanAnomalyFromTrueAnomaly(nu0, orbit.eccentricity);
const a = Math.abs(orbit.semiMajorAxis);
const mu = attractor.stdGravParam;
const n = Math.sqrt(mu / (a*a*a));
const MNext = M0 + n * deltaTime;
const nuNext = solveTrueAnomalyFromMeanAnomaly(orbit.eccentricity, MNext);
return stateFromOrbitElements(orbit, mu, nuNext);
}
function calculateOrbitPeriod(a: number, mu: number){
return TWO_PI * Math.sqrt(a*a*a/mu);
}
function getHohmannPeriod(body1: IOrbitingBody, body2: IOrbitingBody, attractor: ICelestialBody){
const orbit1 = body1.orbit;
const orbit2 = body2.orbit;
const {stdGravParam} = attractor;
const a = (orbit1.semiMajorAxis + orbit2.semiMajorAxis) * 0.5;
const period = calculateOrbitPeriod(a, stdGravParam);
return period;
}

829
src/dedicated-workers/libs/trajectory-calculator.ts
View File

@@ -1,444 +1,535 @@
class TrajectoryCalculator {
public readonly steps: TrajectoryStep[] = [];
private readonly _mainAttractor!: ICelestialBody;
private readonly _legs: LegInfo[] = [];
public totalDeltaV: number = 0;
public noError!: boolean;
public steps: TrajectoryStep[] = [];
public depAltitude!: number;
public startDateMin!: number;
public startDateMax!: number;
public params!: Agent;
private _vesselState!: OrbitalState3D;
private _fbBodyState!: OrbitalState3D;
private _preDSMState!: OrbitalState3D;
private _missionTime: number = 0;
private _sequenceIndex: number = 0;
private _preDSMState!: OrbitalState3D;
private _flybyIncomingState!: OrbitalState3D;
private _flybyBodyState!: OrbitalState3D;
private _depAltitude!: number;
private _startDateMin!: number;
private _startDateMax!: number;
private _params!: Agent;
private _departureDate: number = 0;
private _departureDVScale: number = 0;
private _bodiesOrbits!: OrbitalElements3D[];
constructor(
public readonly system: IOrbitingBody[],
public readonly config: TrajectorySearchSettings,
public readonly sequence: number[],
) {
private _legs: LegInfo[] = [];
private _flybys: FlybyInfo[] = [];
private _departureInfos!: DepartureInfo;
private _secondArcsData: SecondArcData[] = [];
constructor(public readonly system: IOrbitingBody[], public readonly config: TrajectorySearchSettings, public readonly sequence: number[]){
const attractorId = this._departureBody.orbiting;
this._mainAttractor = this.system[attractorId];
}
public setParameters(
depAltitude: number,
startDateMin: number,
startDateMax: number,
params: Agent,
){
this.depAltitude = depAltitude;
this.startDateMin = startDateMin;
this.startDateMax = startDateMax;
this.params = params;
this._getDepartureInfo();
this._clampDepartureInfos();
this._getLegsInfo();
}
public compute(){
this.noError = false;
this._missionTime = this._departureDate;
this._calculateParkingOrbit();
this._calculateEjectionOrbit();
/*if(this._lastStep.duration < 0){
return;
}*/
this._missionTime += this._lastStep.duration;
for(let i = 1; i < this.sequence.length; i++) {
this._sequenceIndex = i;
const leg = this._legs[i-1];
if(!this._calculateLegFirstArc()){
return;
}
this._calculateLegSecondArc();
this._updateLegInfos(i-1);
this._missionTime += leg.duration;
this._calculateFlybyTrajectory();
/*if(this._lastStep.duration < 0){
return;
}*/
this._missionTime += this._lastStep.duration;
if(i == this.sequence.length - 1){
this._calculateArrivalCircularization();
}
}
this.noError = true;
}
private get _destinationBody() {
const id = this.sequence[this.sequence.length - 1];
return this.system[id];
}
private get _departureBody(){
return this.system[this.sequence[0]];
}
private get _attractorMu(){
return this._mainAttractor.stdGravParam;
private get _destinationBody(){
const last = this.sequence.length - 1;
return this.system[this.sequence[last]];
}
private get _lastStep() {
private get _lastStep(){
return this.steps[this.steps.length - 1];
}
private get _lastStepDateOfEnd(){
const {dateOfStart, duration} = this._lastStep;
private get _lastStepEndDate(){
const last = this._lastStep;
const {dateOfStart, duration} = last;
return dateOfStart + duration;
}
/**
* Gets the departure information to be used in a more convenient way.
*/
private _getDepartureInfo(){
this._departureDate = this.params[0];
this._departureDVScale = this.params[1];
public addPrecomputedOrbits(bodiesOrbits: OrbitalElements3D[]){
this._bodiesOrbits = bodiesOrbits;
}
/**
* Apply constraints on the trajectory parameters. The DE may make the agents escape the initial constraints,
* thus the need to apply the constraints on each trajectory evaluation.
public setParameters(depAltitude: number, startDateMin: number, startDateMax: number, params: Agent){
this._depAltitude = depAltitude;
this._startDateMin = startDateMin;
this._startDateMax = startDateMax;
this._params = params;
this._getDepartureSettings();
this._getLegsSettings();
this._getFlybySettings();
}
public reset(){
this._secondArcsData = [];
this._legs = [];
this._flybys = [];
this.steps = [];
}
/**
* Extracts the departures settings from the parameters to be used in a more convenient way.
*/
private _clampDepartureInfos(){
this._departureDate = clamp(this._departureDate, this.startDateMin, this.startDateMax);
const {depDVScaleMin, depDVScaleMax} = this.config;
this._departureDVScale = clamp(this._departureDVScale, depDVScaleMin, depDVScaleMax);
this.params[0] = this._departureDate;
this.params[1] = this._departureDVScale;
private _getDepartureSettings(){
this._departureInfos = {
dateParam: this._params[0],
phaseParam: this._params[1],
ejVelParam: this._params[2]
};
}
/**
* Retrieves the leg informations for the parameters to be used in a more convenient way.
*/
private _getLegsInfo(){
for(let i = 2; i < this.params.length; i += 4) {
this._legs.push({
duration: this.params[i],
dsmOffset: this.params[i+1],
theta: this.params[i+2],
phi: this.params[i+3]
});
private _getLegsSettings(){
const {dsmOffsetMin, dsmOffsetMax} = this.config;
for(let i = 1; i < this.sequence.length; i++) {
const idx = 3 + (i-1)*4;
const infos = {
exitedBodyId: this.sequence[i-1],
targetBodyId: this.sequence[i],
durationParam: this._params[idx],
duration: 0,
dsmParam: lerp(dsmOffsetMin, dsmOffsetMax, this._params[idx+1])
};
this._computeLegDuration(infos);
this._legs.push(infos);
}
}
/** Retrieves the flyby information from the parameters to be used in a more convenient way.
*/
private _getFlybySettings(){
for(let i = 1; i < this.sequence.length - 1; i++) {
const idx = 3 + (i-1)*4 + 2;
const infos = {
flybyBodyId: this.sequence[i],
normAngleParam: this._params[idx],
periRadiParam: this._params[idx+1]
};
this._flybys.push(infos);
}
}
/**
* Updates the parameters of the trajectory (direct agent access) to apply eventual modification
* relative to constraints, which can only be calculated when the trajectory is calculated.
* @param legIndex The index of the leg being updated.
* Computes the whole trajectory with its steps and maneuvers, and stores
* the total delta V.
*/
private _updateLegInfos(legIndex: number){
const leg = this._legs[legIndex];
this.params[2 + legIndex * 4] = leg.duration;
this.params[2 + legIndex * 4 + 1] = leg.dsmOffset;
this.params[2 + legIndex * 4 + 2] = leg.theta;
this.params[2 + legIndex * 4 + 3] = leg.phi;
}
public compute(){
// Compute the departure circular orbit and ejection orbit
this._appendDepartureParkingOrbit();
this._computeDepartureEjectionOrbit();
/**
* Calculates the circularization orbit and maneuver at the arrival body
*/
private _calculateArrivalCircularization(){
const body = this._destinationBody;
const flybyOrbit = this._lastStep.orbitElts;
const periapsisState = stateFromOrbitElements(flybyOrbit, body.stdGravParam, 0);
const progradeDir = normalize3(periapsisState.vel);
const circularVel = circularVelocity(body, mag3(periapsisState.pos));
const numOfLegs = this._legs.length;
// Compute all the interplanetary arcs and flybys
for(let i = 0; i < numOfLegs - 1; i++){
const leg = this._legs[i];
const flyby = this._flybys[i];
this._computeFirstLegArc(leg);
this._computeLegSecondArcSimple(leg);
this._computeFlyby(flyby);
}
const deltaVMag = Math.abs(circularVel - mag3(periapsisState.vel));
const deltaV = mult3(progradeDir, -deltaVMag);
this.totalDeltaV += deltaVMag;
const circularState = {
pos: periapsisState.pos,
vel: mult3(progradeDir, circularVel)
}
const circularOrbit = stateToOrbitElements(circularState, body);
const maneuvre: ManeuvreInfo = {
deltaVToPrevStep: deltaV,
progradeDir: progradeDir,
manoeuvrePosition: periapsisState.pos,
context: {
type: "circularization",
}
};
this.steps.push({
orbitElts: circularOrbit,
attractorId: body.id,
beginAngle: 0,
endAngle: TWO_PI,
duration: 0,
dateOfStart: this._lastStep.dateOfStart + this._lastStep.duration,
maneuvre: maneuvre
});
// Compute the last leg reaching the destination body
// with no flyby
const last = this._legs[numOfLegs-1];
this._computeFirstLegArc(last);
this._computeLegSecondArcSimple(last);
}
/**
* Calculates the swing-by orbit from the external incoming state.
* If the body is the destination body of the sequence, a deltaV is calculated as the difference
* in velocity needed to circularize around the body at the arrival radius.
*/
private _calculateFlybyTrajectory(){
const body = this.system[this.sequence[this._sequenceIndex]];
const localIncomingState = {
pos: sub3(this._flybyIncomingState.pos, this._flybyBodyState.pos),
vel: sub3(this._flybyIncomingState.vel, this._flybyBodyState.vel)
};
const flybyOrbit = stateToOrbitElements(localIncomingState, body);
const isDestinationBody = body.id == this._destinationBody.id;
let nu1 = trueAnomalyFromOrbitalState(flybyOrbit, localIncomingState);
let nu2 = isDestinationBody ? 0 : -nu1;
const tof = tofBetweenAnomalies(flybyOrbit, body, nu1, nu2);
this.steps.push({
orbitElts: flybyOrbit,
attractorId: body.id,
beginAngle: nu1,
endAngle: nu2,
duration: tof,
dateOfStart: this._missionTime
});
}
/**
* Calculates the second arc of a leg : solves the Lambert problem to join the pre-DSM state
* and the targeted body from the sequence. The parameters relative to the leg are clamped
* to ensure a feasible trajectory and respect of the constraints.
*/
private _calculateLegSecondArc(){
const leg = this._legs[this._sequenceIndex-1];
const preDSMState = this._preDSMState;
const origin = this.system[this.sequence[this._sequenceIndex-1]];
const target = this.system[this.sequence[this._sequenceIndex]];
const startPos = preDSMState.pos;
const encounterDate = this._missionTime + leg.duration;
const targetState = getBodyStateAtDate(target, encounterDate, this._attractorMu);
const arcDuration = (1 - leg.dsmOffset) * leg.duration;
const incomingVel = solveLambert(startPos, targetState.pos, arcDuration, this._mainAttractor).v2;
const relativeVel = sub3(incomingVel, targetState.vel);
const incomingDir = normalize3(mult3(relativeVel, -1));
// Calculate the position of the entrance point on the SOI
const {thetaMin, thetaMax} = minMaxRingAngle(target.soi, target.radius, target.radius * 2);
const pointIsValid = leg.theta >= thetaMin && leg.theta <= thetaMax;
if(!pointIsValid) {
const {theta, phi} = randomPointOnSphereRing(thetaMin, thetaMax);
leg.theta = theta;
leg.phi = phi;
}
const soiPointLocal = mult3(pointOnSphereRing(incomingDir, leg.theta, leg.phi), target.soi);
const soiPoint = add3(targetState.pos, soiPointLocal);
// Solve lambert problem to reach this point
const {v1, v2} = solveLambert(startPos, soiPoint, arcDuration, this._mainAttractor);
const secondLegArc = stateToOrbitElements({pos: startPos, vel: v1}, this._mainAttractor);
let nu1 = trueAnomalyFromOrbitalState(secondLegArc, {pos: startPos, vel: v1});
let nu2 = trueAnomalyFromOrbitalState(secondLegArc, {pos: soiPoint, vel: v2});
if(nu1 > nu2) {
nu2 = TWO_PI + nu2;
}
this.totalDeltaV += mag3(sub3(v1, preDSMState.vel));
// Calculate the DSM maneuvre
const maneuvre: ManeuvreInfo = {
deltaVToPrevStep: sub3(v1, preDSMState.vel),
progradeDir: normalize3(preDSMState.vel),
manoeuvrePosition: preDSMState.pos,
context: {
type: "dsm",
originId: origin.id,
targetId: target.id
public get totalDeltaV(){
let total = 0;
for(const step of this.steps){
if(step.maneuvre) {
const {maneuvre} = step;
const {deltaVToPrevStep} = maneuvre;
const dv = mag3(deltaVToPrevStep);
total += dv;
}
}
this.steps.push({
orbitElts: secondLegArc,
attractorId: this._mainAttractor.id,
beginAngle: nu1,
endAngle: nu2,
duration: arcDuration,
dateOfStart: this._lastStep.dateOfStart + this._lastStep.duration,
maneuvre: maneuvre
});
this._flybyBodyState = targetState;
this._flybyIncomingState = {pos: soiPoint, vel: v2};
return total;
}
/**
* Propagates the orbit from the last body's SOI exit point to the date of the DSM.
* The parameters relative to the leg are clamped to ensure a feasible trajectory and respect
* of the constraints.
* @returns Returns whether the propagation returns a valid result
* Completes the provided leg infos by calculating the leg duration from the already given
* parameters
* @param infos The leg infos to complete
*/
private _calculateLegFirstArc() {
const seqIndex = this._sequenceIndex;
const leg = this._legs[seqIndex-1];
const isResonant = this.sequence[seqIndex-1] == this.sequence[seqIndex];
private _computeLegDuration(infos: LegInfo){
const exitedBody = this.system[infos.exitedBodyId];
const {durationParam} = infos;
const attr = this._mainAttractor;
if(infos.exitedBodyId != infos.targetBodyId) { // if the leg is not resonant
// The flight duration of the leg is between 35% and 75% of an ideal elliptic transfer
// orbit from the exited body to the target body
const targetBody = this.system[infos.targetBodyId];
const exBodyA = exitedBody.orbit.semiMajorAxis;
const tgBodyA = targetBody.orbit.semiMajorAxis;
const a = 0.5 * (exBodyA + tgBodyA);
const period = Physics3D.orbitPeriod(attr, a);
infos.duration = lerp(0.35*period, 0.75*period, durationParam);
} else { // if the leg is resonant
// The flight duration of the leg is between 100% and 200% the orbital period of the
// exited body
const a = exitedBody.orbit.semiMajorAxis;
const period = Physics3D.orbitPeriod(attr, a);
infos.duration = lerp(1*period, 2*period, durationParam);
}
const {minLegDuration} = this.config;
infos.duration = Math.max(minLegDuration, infos.duration);
}
/**
* Recomputes the second arc of a all legs accounting this time for SOI enter points,
* thus providing a more precise visualization of the orbit.
*/
public recomputeLegsSecondArcs(){
for(let i = 0; i < this._secondArcsData.length; i++){
const data = this._secondArcsData[i];
const step = this.steps[3 + i*3];
this._recomputeSecondArc(step, data);
}
}
/**
* Recomputes the second arc of a leg accounting this time for SOI enter points.
* @param lambertStep The step when the lambert arc to this flyby is described
* @param arcData The data related to the lambert arc and flyby stored during the flyby calculation
*/
private _recomputeSecondArc(lambertStep: TrajectoryStep, arcData: SecondArcData){
const fbBody = this.system[arcData.fbBodyId];
const {flybyOrbit, soiEnterAngle} = arcData;
const localEnterState = Physics3D.orbitElementsToState(flybyOrbit, fbBody, soiEnterAngle);
const {fbBodyState} = arcData;
const targetEnterPos = add3(fbBodyState.pos, localEnterState.pos);
const {preDSMState} = arcData;
const {v1, v2} = Lambert.solve(preDSMState.pos, targetEnterPos, lambertStep.duration, this._mainAttractor);
const postDSMState = {pos: preDSMState.pos, vel: v1};
const encounterState = {pos: targetEnterPos, vel: v2};
// Compute the orbit associated to this Lambert arc
const arcOrbit = Physics3D.stateToOrbitElements(postDSMState, this._mainAttractor);
// Compute the begin and end angles of the arc
const angles = {
begin: Physics3D.trueAnomalyFromOrbitalState(arcOrbit, postDSMState),
end: Physics3D.trueAnomalyFromOrbitalState(arcOrbit, encounterState)
};
const drawAngles = {
begin: angles.begin,
end: angles.end
};
if(drawAngles.begin > drawAngles.end){
drawAngles.end += TWO_PI;
}
lambertStep.orbitElts = arcOrbit;
lambertStep.angles = angles;
lambertStep.drawAngles = drawAngles;
const maneuvre = lambertStep.maneuvre as ManeuvreInfo;
const dsmDV = sub3(postDSMState.vel, preDSMState.vel);
maneuvre.deltaVToPrevStep = dsmDV;
}
/**
* Computes the flyby orbit with the provided parameters.
* @param flybyInfo The flyby parameters
*/
private _computeFlyby(flybyInfo: FlybyInfo){
const body = this.system[flybyInfo.flybyBodyId];
// Compute the relative velocity to the body when entering its SOI
const bodyVel = this._fbBodyState.vel;
const globalIncomingVel = this._vesselState.vel;
const relativeIncomingVel = sub3(globalIncomingVel, bodyVel);
const incomingVelMag = mag3(relativeIncomingVel);
const incomingVelDir = div3(relativeIncomingVel, incomingVelMag);
// Compute the normal vector of the orbit : we compute the vector perpendicular to the incoming
// relative velocity and the body velocity (relative to the main attractor), and rotate it
// around the incoming relative velocity direction by angle angle defined in the parameters.
const normAngle = lerp(0, TWO_PI, flybyInfo.normAngleParam);
const t_normal = normalize3(cross(relativeIncomingVel, bodyVel));
const normal = rotate3(t_normal, incomingVelDir, normAngle);
// Compute the temporary periapsis position vector, it is colinear to the incoming velcity
// vector
const t_periDir = incomingVelDir;
const {fbRadiusMaxScale} = this.config;
const periRadius = lerp(body.radius, fbRadiusMaxScale * body.radius, flybyInfo.periRadiParam);
const t_periPos = mult3(t_periDir, periRadius);
// Compute the temporary periapsis velocity velocity, at this point on the orbit,
// the velocity vector is perpendicular to the periapsis vector.
const periVelMag = Physics3D.deduceVelocityAtRadius(body, body.soi, incomingVelMag, periRadius);
const t_periVelDir = rotate3(t_periDir, normal, HALF_PI);
const t_periVel = mult3(t_periVelDir, periVelMag);
// Compute the temporary flyby orbit
const t_periapsisState = {pos: t_periPos, vel: t_periVel};
const t_flybyOrbit = Physics3D.stateToOrbitElements(t_periapsisState, body);
// Compute the SOI enter and exit angles
const exitAngle = Physics3D.trueAnomalyAtRadius(t_flybyOrbit, body.soi);
const angles = {begin: -exitAngle, end: exitAngle};
const drawAngles = {begin: angles.begin, end: angles.end};
// Compute the angle between the incoming velocity vector and the velocity vector
// on the calculated orbit at the SOI enter point
const t_incomingVel = Physics3D.orbitElementsToState(t_flybyOrbit, body, angles.begin).vel;
const t_incomingVelDir = normalize3(t_incomingVel);
const dotCross = dot3(cross(incomingVelDir, t_incomingVelDir), normal);
const dotVel = dot3(incomingVelDir, t_incomingVelDir);
const offsetAngle = Math.atan2(dotCross, dotVel);
// Recompute the flyby orbit by rotating the periapsis so that the incoming velocity calculated
// on the orbit matches the incoming velocity computed at the beginning.
// FIX: When `normAngle` is exactly 0 or TWO_PI, we must rotate by offsetAngle
// and not -offsetAngle... Why this happens ? I don't know, probably some floating point error
// stuff...
const periPos = rotate3(t_periPos, normal, -offsetAngle);
const periVel = rotate3(t_periVel, normal, -offsetAngle);
const periapsisState = {pos: periPos, vel: periVel};
const flybyOrbit = Physics3D.stateToOrbitElements(periapsisState, body);
// Compute the duration of the flyby
const tof = Physics3D.tofBetweenAnomalies(flybyOrbit, body, angles.begin, angles.end);
// Append the flyby orbit
this.steps.push({
orbitElts: flybyOrbit,
attractorId: body.id,
angles: angles,
drawAngles: drawAngles,
duration: tof,
dateOfStart: this._lastStepEndDate
});
// Compute the vessel state relative to the body when exiting its SOI
const exitState = Physics3D.orbitElementsToState(flybyOrbit, body, exitAngle);
this._vesselState = exitState;
// Store the required information for eventual further recalculations of the
// legs' lambert arcs
this._secondArcsData.push({
preDSMState: this._preDSMState,
fbBodyId: body.id,
fbBodyState: this._fbBodyState,
flybyOrbit: flybyOrbit,
soiEnterAngle: angles.begin
});
}
/**
* Computes the second leg arc using a Lambert solver to aim at the body position
* supposing a null radius SOI.
* @param legInfo The leg parameters
*/
private _computeLegSecondArcSimple(legInfo: LegInfo){
const attr = this._mainAttractor;
const lastStep = this._lastStep;
const exitedBody = this.system[lastStep.attractorId];
const localState = stateFromOrbitElements(lastStep.orbitElts, exitedBody.stdGravParam, lastStep.endAngle);
const exitedBodyState = getBodyStateAtDate(exitedBody, this._missionTime, this._attractorMu);
const exitState = {
pos: add3(localState.pos, exitedBodyState.pos),
vel: add3(localState.vel, exitedBodyState.vel)
const preDSMState = this._vesselState;
// Compute the state of the body at the encounter date
const encounterDate = lastStep.dateOfStart + legInfo.duration;
const targetBody = this.system[legInfo.targetBodyId];
const tgBodyOrbit = this._bodiesOrbits[targetBody.id];
const tgBodyState = Physics3D.bodyStateAtDate(targetBody, tgBodyOrbit, attr, encounterDate);
// Solve the Lambert problem to reach the body.
// We suppose for now that the body has a null radius SOI, and so we target directly the body.
const arcDuration = legInfo.duration * (1 - legInfo.dsmParam);
const {v1, v2} = Lambert.solve(preDSMState.pos, tgBodyState.pos, arcDuration, attr);
const postDSMState = {pos: preDSMState.pos, vel: v1};
const encounterState = {pos: tgBodyState.pos, vel: v2};
// Compute the orbit associated to this Lambert arc
const arcOrbit = Physics3D.stateToOrbitElements(postDSMState, attr);
// Compute the begin and end angles of the arc
const angles = {
begin: Physics3D.trueAnomalyFromOrbitalState(arcOrbit, postDSMState),
end: Physics3D.trueAnomalyFromOrbitalState(arcOrbit, encounterState)
};
const firstLegArc = stateToOrbitElements(exitState, this._mainAttractor);
if(firstLegArc.semiMajorAxis >= this._mainAttractor.soi){
return false;
const drawAngles = {
begin: angles.begin,
end: angles.end
};
if(drawAngles.begin > drawAngles.end){
drawAngles.end += TWO_PI;
}
const {dsmOffsetMax, dsmOffsetMin, minLegDuration} = this.config;
if(firstLegArc.eccentricity < 1) {
const orbitPeriod = calculateOrbitPeriod(firstLegArc.semiMajorAxis, this._attractorMu);
let minScale = 0.5;
let maxScale = 1;
if(isResonant){
minScale = 1;
maxScale = 3;
// Compute DSM maneuver
const progradeDir = normalize3(preDSMState.vel);
const dsmDV = sub3(postDSMState.vel, preDSMState.vel);
const maneuvre: ManeuvreInfo = {
position: preDSMState.pos,
deltaVToPrevStep: dsmDV,
progradeDir: progradeDir,
context: {
type: "dsm",
originId: legInfo.exitedBodyId,
targetId: legInfo.targetBodyId
}
const minDuration = minScale * orbitPeriod;
const maxDuration = maxScale * orbitPeriod;
leg.duration = clamp(leg.duration, minDuration, maxDuration);
const revolutions = leg.duration / orbitPeriod;
const minDSMOffset = Math.max((revolutions - 1) / revolutions, dsmOffsetMin);
leg.dsmOffset = clamp(leg.dsmOffset, minDSMOffset, dsmOffsetMax);
}
leg.duration = clamp(leg.duration, minLegDuration, Infinity);
leg.dsmOffset = clamp(leg.dsmOffset, dsmOffsetMin, dsmOffsetMax);
const arcDuration = leg.dsmOffset * leg.duration;
let legStartNu = trueAnomalyFromOrbitalState(firstLegArc, exitState);
const preDSMState = propagateState(firstLegArc, legStartNu, arcDuration, this._mainAttractor);
let dsmNu = trueAnomalyFromOrbitalState(firstLegArc, preDSMState);
if(legStartNu > dsmNu) {
dsmNu += TWO_PI;
}
if(isNaN(legStartNu) || isNaN(dsmNu)){
return false;
}
};
// Append the second arc orbit
this.steps.push({
orbitElts: firstLegArc,
attractorId: this._mainAttractor.id,
beginAngle: legStartNu,
endAngle: dsmNu,
duration: arcDuration,
dateOfStart: this._lastStep.dateOfStart + this._lastStep.duration
orbitElts: arcOrbit,
attractorId: attr.id,
angles: angles,
drawAngles: drawAngles,
duration: arcDuration,
dateOfStart: this._lastStepEndDate,
maneuvre: maneuvre
});
this._preDSMState = preDSMState;
return true;
// Store the incoming state of the vessel (relative to the main attractor) at encounter
// and the flyby body for the flyby calculation
this._vesselState = encounterState;
this._fbBodyState = tgBodyState;
}
/**
* Calculates the departure body's ejection orbit : hyperbolic orbit from the circular
* parking orbit to SOI exit.
* Computes the next leg first arc with the provided parameters.
* @param legInfo The parameters of the leg
*/
private _calculateEjectionOrbit(){
const startBody = this.system[this.sequence[0]];
const nextBody = this.system[this.sequence[1]];
private _computeFirstLegArc(legInfo: LegInfo){
// Compute the state of the vessel relatived to the exited body after exiting its SOI
/*const lastStep = this._lastStep;
const exitedBody = this.system[lastStep.attractorId];
const localExitState = Physics3D.orbitElementsToState(lastStep.orbitElts, exitedBody, lastStep.angles.end);*/
const localExitState = this._vesselState;
const exitedBody = this.system[this._lastStep.attractorId];
const parkingRadius = startBody.radius + this.depAltitude;
const parkingVel = circularVelocity(startBody, parkingRadius);
const ejectionVel = this._departureDVScale * ejectionVelocity(startBody, parkingRadius);
const ejectionDV = ejectionVel - parkingVel;
this.totalDeltaV += ejectionDV;
const bodyPos = getBodyStateAtDate(startBody, this._departureDate, this._attractorMu).pos;
const r1 = startBody.orbit.semiMajorAxis;
const r2 = nextBody.orbit.semiMajorAxis;
let {tangent, nu} = idealEjectionDirection(bodyPos, r2 < r1);
const ejectionState = {
pos: vec3(parkingRadius * Math.cos(nu), 0, parkingRadius * Math.sin(nu)),
vel: mult3(tangent, ejectionVel),
};
const offsetAngle = hyperbolicEjectionOffsetAngle(ejectionVel, parkingRadius, startBody);
const up = vec3(0, 1, 0);
ejectionState.pos = rotate3(ejectionState.pos, up, -offsetAngle);
ejectionState.vel = rotate3(ejectionState.vel, up, -offsetAngle);
const ejectionOrbit = stateToOrbitElements(ejectionState, startBody);
const exitAnomaly = soiExitTrueAnomaly(ejectionOrbit, startBody.soi);
const tof = tofBetweenAnomalies(ejectionOrbit, startBody, 0, exitAnomaly);
// Compute the exited body state relative to its attractor at the moment the vessel exits its SOI
const bodyOrbit = this._bodiesOrbits[exitedBody.id];
const exitDate = this._lastStepEndDate;
const exitedBodyState = Physics3D.bodyStateAtDate(exitedBody, bodyOrbit, this._mainAttractor, exitDate);
// Calculate the maneuvre informations
const progradeDir = normalize3(ejectionState.vel);
const maneuvre: ManeuvreInfo = {
deltaVToPrevStep: mult3(progradeDir, ejectionDV),
progradeDir: progradeDir,
manoeuvrePosition: ejectionState.pos,
context: {type: "ejection"}
// Compute the vessel SOI-exit state relative to the main attractor (exited body's attractor)
const exitState = {
pos: add3(exitedBodyState.pos, localExitState.pos),
vel: add3(exitedBodyState.vel, localExitState.vel),
};
// Compute the orbit of the vessel around the main attractor after exiting the exited body SOI.
const arcOrbit = Physics3D.stateToOrbitElements(exitState, this._mainAttractor);
const beginAngle = Physics3D.trueAnomalyFromOrbitalState(arcOrbit, exitState);
// Propagate the vessel on the arc for the arc duration
const arcDuration = legInfo.dsmParam * legInfo.duration;
const preDSMState = Physics3D.propagateStateFromTrueAnomaly(arcOrbit, this._mainAttractor, beginAngle, arcDuration);
let endAngle = Physics3D.trueAnomalyFromOrbitalState(arcOrbit, preDSMState);
const angles = {begin: beginAngle, end: endAngle};
// Compute the ends of the arc to show in the 3D view if the trajectory happens to be shown
const drawAngles = {begin: angles.begin, end: angles.end};
// It is possible that the duration of the arc is longer than the period of the arc orbit.
// In that case we must display all the orbit but keep the last angle meaningful.
const period = Physics3D.orbitPeriod(this._mainAttractor, arcOrbit.semiMajorAxis);
if(arcDuration > period){
drawAngles.begin = 0;
drawAngles.end = TWO_PI;
} else if(beginAngle > endAngle) {
drawAngles.end += TWO_PI;
}
// Append the arc to the steps
this.steps.push({
orbitElts: ejectionOrbit,
attractorId: startBody.id,
beginAngle: 0,
endAngle: exitAnomaly,
duration: tof,
dateOfStart: this._lastStepDateOfEnd,
maneuvre: maneuvre
orbitElts: arcOrbit,
attractorId: this._mainAttractor.id,
angles: angles,
drawAngles: drawAngles,
duration: arcDuration,
dateOfStart: exitDate,
});
// Store the pre-DSM state to use it in the second arc maneuver calculation
this._vesselState = preDSMState;
this._preDSMState = preDSMState;
}
/**
* Computes the ejection arc from the departure body parking orbit to exiting its SOI.
*/
private _computeDepartureEjectionOrbit(){
const curOrbit = this._lastStep.orbitElts;
const depBody = this._departureBody;
const {phaseParam, ejVelParam} = this._departureInfos;
// Compute the phase angle, measured from the vector pointing to the body's
// attractor
const phase = lerp(0, TWO_PI, phaseParam);
// Compute the ejection velocity magnitude, i.e. the velocity at the parking
// orbit altitude which leads to exit the SOI with the specified
// scale given in parameters.
const ejVelMag0 = Physics3D.velocityToReachAltitude(depBody, curOrbit.semiMajorAxis, depBody.soi);
const {depDVScaleMin, depDVScaleMax} = this.config;
const ejVelMag = lerp(depDVScaleMin, depDVScaleMax, ejVelParam) * ejVelMag0;
// Compute the ejection velocity vector; rotated by the phase angle relative to the
// reference vector (1, 0, 0)
const {vel, pos} = Physics3D.orbitElementsToState(curOrbit, depBody, phase);
const progradeDir = normalize3(vel);
const ejVel = mult3(progradeDir, ejVelMag);
// Compute the ejection orbit
const ejOrbit = Physics3D.stateToOrbitElements({pos, vel: ejVel}, depBody);
// Compute the angle at which the SOI is exited
const soiExitAngle = Physics3D.trueAnomalyAtRadius(ejOrbit, depBody.soi);
// Compute the ejection maneuvre
const ejDV = ejVelMag - mag3(vel);
const maneuvre: ManeuvreInfo = {
position: pos,
deltaVToPrevStep: mult3(progradeDir, ejDV),
progradeDir: progradeDir,
context: {type: "ejection"}
};
// Compute the duration of the flight from the parking orbit to SOI exit.
const tof = Physics3D.tofBetweenAnomalies(ejOrbit, depBody, 0, soiExitAngle);
// Append the ejection orbit
this.steps.push({
orbitElts: ejOrbit,
attractorId: depBody.id,
angles: {begin: 0, end: soiExitAngle},
drawAngles: {begin: 0, end: soiExitAngle},
duration: tof,
dateOfStart: this._lastStepEndDate,
maneuvre: maneuvre
});
// Compute the state of the vessel relative to the exited body when exiting its SOI
const exitState = Physics3D.orbitElementsToState(ejOrbit, depBody, soiExitAngle);
this._vesselState = exitState;
}
/**
* Calculates the orbital elements of the circular parking orbit
* around the departure body.
* around the departure body and appends it to the steps.
*/
private _calculateParkingOrbit(){
const startBody = this.system[this.sequence[0]];
const orbitRadius = startBody.radius + this.depAltitude;
private _appendDepartureParkingOrbit(){
const {dateParam} = this._departureInfos;
const radius = this._departureBody.radius + this._depAltitude;
const dateMin = this._startDateMin;
const dateMax = this._startDateMax;
this.steps.push({
orbitElts: equatorialCircularOrbit(orbitRadius),
attractorId: startBody.id,
beginAngle: 0,
endAngle: TWO_PI,
duration: 0,
dateOfStart: this._departureDate
orbitElts: Physics3D.equatorialCircularOrbit(radius),
attractorId: this._departureBody.id,
angles: {begin: 0, end: 0},
drawAngles: {begin: 0, end: TWO_PI},
duration: 0,
dateOfStart: lerp(dateMin, dateMax, dateParam)
});
}
}

2
src/dedicated-workers/sequence-evaluator.ts
View File

@@ -154,4 +154,4 @@ class SequenceEvaluator extends WorkerEnvironment {
}
}
initWorker(SequenceEvaluator);
WorkerEnvironment.init(SequenceEvaluator);

2
src/dedicated-workers/sequence-generator.ts
View File

@@ -124,4 +124,4 @@ class SequenceGenerator extends WorkerEnvironment {
}
}
initWorker(SequenceGenerator);
WorkerEnvironment.init(SequenceGenerator);

229
src/dedicated-workers/trajectory-optimizer.ts
View File

@@ -3,18 +3,30 @@ importScripts("libs/common.js", "libs/evolution.js", "libs/math.js", "libs/physi
class TrajectoryOptimizer extends WorkerEnvironment {
private _config!: Config;
private _system!: IOrbitingBody[];
private _bodiesOrbits!: OrbitalElements3D[];
private _depAltitude!: number;
private _sequence!: number[];
private _startDateMin!: number;
private _startDateMax!: number;
private _bestDeltaV!: number;
private _bestSteps!: TrajectoryStep[];
private _bestTrajectory!: TrajectoryCalculator;
private _bestDeltaV!: number;
private _evolver!: ChunkedEvolver;
override onWorkerInitialize(data: any){
this._config = data.config;
this._system = data.system;
// Precompute bodies' orbital elements
//@ts-ignore
this._bodiesOrbits = [null];
for(let i = 1; i < this._system.length; i++) {
const data = this._system[i].orbit;
const orbit = Physics3D.orbitElementsFromOrbitData(data);
this._bodiesOrbits.push(orbit);
}
}
override onWorkerDataPass(data: any){
@@ -25,158 +37,119 @@ class TrajectoryOptimizer extends WorkerEnvironment {
}
override onWorkerRun(input: any){
const evaluate = (params: Agent) => this._evaluate(params);
if(input.start){
// If it's the first generation, we generate configure the evolver and generate
// a new random population.
const numLegs = this._sequence.length - 1;
const agentDim = 3 + numLegs*4 - 2;
const fitness = (agent: Agent) => {
// The fitness function calculates the trajectory represented
// by the agent and update the best trajectory found yet.
const trajectory = this._computeTrajectory(agent);
if(trajectory.totalDeltaV < this._bestDeltaV){
this._bestDeltaV = trajectory.totalDeltaV;
this._bestTrajectory = trajectory;
}
return trajectory.totalDeltaV;
};
const trajConfig = this._config.trajectorySearch;
const {crossoverProba, diffWeight} = trajConfig;
const {chunkStart, chunkEnd} = input;
this._evolver = new ChunkedEvolver(
chunkStart, chunkEnd,
agentDim, fitness, crossoverProba, diffWeight
);
if(input.start) {
this._bestDeltaV = Infinity;
const chunkSize = input.chunkEnd - input.chunkStart + 1;
const popChunk = this._createRandomMGAPopulationChunk(chunkSize);
const dvChunk = populationChunkFitness(popChunk, evaluate);
// Create the first generation and evaluate it
const popChunk = this._evolver.createRandomPopulationChunk();
const dvChunk = this._evolver.evaluateChunkFitness(popChunk);
sendResult({
bestSteps: this._bestSteps,
bestSteps: this._bestTrajectory.steps,
bestDeltaV: this._bestDeltaV,
fitChunk: dvChunk,
popChunk: popChunk,
});
} else {
const {crossoverProba, diffWeight} = this._config.trajectorySearch;
const results = evolvePopulationChunk(
input.population,
input.deltaVs,
input.chunkStart,
input.chunkEnd,
crossoverProba, diffWeight, evaluate
);
// If not the first generation, then evolve the current population
const {population, deltaVs} = input;
const {popChunk, fitChunk} = this._evolver.evolvePopulationChunk(population, deltaVs);
sendResult({
...results,
bestSteps: this._bestSteps,
popChunk, fitChunk,
bestSteps: this._bestTrajectory.steps,
bestDeltaV: this._bestDeltaV
});
}
}
/**
* Evaluates a trajectory cost and stores it if it's the best found so far.
* @param params A trajectory paremeters (agent)
* @returns The fitness (total deltaV) of the trajectory, used for the DE algorithm.
* Computes the trajectory represented by an agent. Throws an error if the trajectory fails at
* being computed.
* @param agent An agent representing a trajectory
* @param maxAttempts The maximum number of attempts to compute the trajectory before throwing an error
* @returns The computed trajectory
*/
private _evaluate(params: Agent) {
const trajectory = this._computeTrajectory(params);
if(trajectory.totalDeltaV < this._bestDeltaV){
this._bestDeltaV = trajectory.totalDeltaV;
this._bestSteps = trajectory.steps;
}
return trajectory.totalDeltaV;
};
private _computeTrajectory(agent: Agent, maxAttempts: number = 1000){
const trajConfig = this._config.trajectorySearch;
const trajectory = new TrajectoryCalculator(this._system, trajConfig, this._sequence);
trajectory.addPrecomputedOrbits(this._bodiesOrbits);
/**
* Calculates a trajectory descibed by the given parameters.
* Edits the parameters if needed to satisfy the constraints.
* @param params A trajectory paremeters (agent)
* @returns The calculated trajectory
*/
private _computeTrajectory(params: Agent){
const calculate = () => {
const trajectory = new TrajectoryCalculator(this._system, this._config.trajectorySearch, this._sequence);
trajectory.setParameters(this._depAltitude, this._startDateMin, this._startDateMax, params);
trajectory.compute();
return trajectory;
}
let trajectory = calculate();
while(!trajectory.noError) {
this._randomizeExistingAgent(params);
trajectory = calculate();
}
return trajectory;
}
/**
* Generates a random population chunk of agent representing MGA trajectories based
* on the given body sequence.
* @param chunkSize The size of the population chunk to generate
* @returns A random population chunk.
*/
private _createRandomMGAPopulationChunk(chunkSize: number, ){
const popChunk: Agent[] = [];
for(let i = 0; i < chunkSize; i++) {
popChunk.push(this._createRandomMGAAgent());
}
return popChunk;
}
/**
* Generates a random agent representing an MGA trajectory based on the given sequence.
* @returns A random trajectory agent
*/
private _createRandomMGAAgent() {
const dim = 4 * (this._sequence.length - 1) + 2;
const solution: Agent = Array<number>(dim).fill(0);
// Initial ejection condition on departure body
solution[0] = randomInInterval(this._startDateMin, this._startDateMax); // departure date
solution[1] = randomInInterval(1, 3); // ejection deltaV scale
// Interplanetary legs' parameters
for(let i = 1; i < this._sequence.length; i++){
const j = 2 + (i - 1) * 4;
// There can be errors happening during the calculation (like NaN values),
// therefore we try to calculate the trajectory until the maximum number of attempts is reached.
// If an error occurs, the agent is randomized.
let attempts = 0;
while(attempts < maxAttempts){
trajectory.setParameters(this._depAltitude, this._startDateMin, this._startDateMax, agent);
let failed = false;
// FIX: "This radius is never reached" error thrown... why ?
try {
trajectory.compute();
trajectory.recomputeLegsSecondArcs();
} catch {
failed = true;
}
const {legDuration, dsmOffset} = this._legDurationAndDSM(i);
solution[j] = legDuration;
solution[j+1] = dsmOffset;
if(failed || this._hasNaNValuesInSteps(trajectory)) {
this._evolver.randomizeAgent(agent);
trajectory.reset();
} else {
return trajectory;
}
const {theta, phi} = randomPointOnSphereRing(0, Math.PI/2);
solution[j+2] = theta;
solution[j+3] = phi;
attempts++;
}
return solution;
throw new Error("Impossible to compute the trajectory.");
}
/**
* Randomizes an existing agent instance.
* @param agent The angent to randomize
* Checks if there is a NaN value in the computed steps (caused by a math error)
* @param trajectory The trajectory we want to check its steps for NaN values.
* @returns true if there is a NaN value in the computed steps, false otherwise.
*/
private _randomizeExistingAgent(agent: Agent){
const newAgent = this._createRandomMGAAgent();
for(let i = 0; i < agent.length; i++){
agent[i] = newAgent[i];
}
}
/**
* @param index The index on the planetary sequence.
* @returns A random leg duration and DSM offset
*/
private _legDurationAndDSM(index: number){
const isResonant = this._sequence[index-1] == this._sequence[index];
const {dsmOffsetMin, dsmOffsetMax} = this._config.trajectorySearch;
if(isResonant) {
const sideralPeriod = <number>this._system[this._sequence[index]].orbit.sideralPeriod;
const revs = randomInInterval(1, 4);
const legDuration = revs * sideralPeriod;
const minOffset = clamp((revs - 1) / revs, dsmOffsetMin, dsmOffsetMax);
const dsmOffset = randomInInterval(minOffset, dsmOffsetMax);
return {legDuration, dsmOffset};
} else {
const body1 = this._system[this._sequence[index-1]];
const body2 = this._system[this._sequence[index]];
const attractor = this._system[body1.orbiting];
const period = getHohmannPeriod(body1, body2, attractor);
const legDuration = randomInInterval(0.1, 1) * period;
const dsmOffset = randomInInterval(dsmOffsetMin, dsmOffsetMax);
return {legDuration, dsmOffset};
private _hasNaNValuesInSteps(trajectory: TrajectoryCalculator){
const hasNaN: (obj: Object) => boolean = obj => {
for(const value of Object.values(obj)){
if(typeof value == "object"){
if(hasNaN(value))
return true;
} else if(typeof value == "number") {
if(isNaN(value))
return true;
}
}
return false;
};
const {steps} = trajectory;
for(let i = steps.length - 1; i >= 0; i--){
if(hasNaN(steps[i]))
return true;
}
return false;
}
}
initWorker(TrajectoryOptimizer);
WorkerEnvironment.init(TrajectoryOptimizer);

2
src/main/editor/editor.ts
View File

@@ -163,7 +163,7 @@ export function initEditor(controls: CameraController, system: SolarSystem, conf
const displayFoundTrajectory = () => {
trajectory = new Trajectory(solver.bestTrajectorySteps, system, config);
trajectory.draw(canvas);
trajectory.fillResultControls(maneuvreSelector, resultSpans, stepSlider, systemTime);
trajectory.fillResultControls(maneuvreSelector, resultSpans, stepSlider, systemTime, controls);
maneuvreSelector.select(0);
maneuvreSelector.enable();

2
src/main/objects/orbit.ts
View File

@@ -63,7 +63,7 @@ export class Orbit implements IOrbit {
};
}
public static fromOrbitalElements(elements: OrbitalElements, attractor: ICelestialBody, config: OrbitSettings){
public static fromOrbitalElements(elements: OrbitalElements3D, attractor: ICelestialBody, config: OrbitSettings){
const data = {
eccentricity: elements.eccentricity,
semiMajorAxis: elements.semiMajorAxis,

111
src/main/solvers/trajectory-solver.ts
View File

@@ -23,13 +23,13 @@ export class TrajectorySolver {
this._workerPool = new WorkerPool("dist/dedicated-workers/trajectory-optimizer.js", this.config);
this._workerPool.initialize({system: this.system.data, config: this.config});
}
private _initPlot(){
this.plot.clearPlot();
}
private _updatePlot(iteration: number){
const {best, mean} = this._bestMeanDeltaV;
let mean = 0;
for(const dv of this._deltaVs)
mean += dv;
mean /= this.popSize;
const best = this.bestDeltaV;
this.plot.addIterationData(iteration, mean, best);
}
@@ -37,18 +37,10 @@ export class TrajectorySolver {
if(this._running) this._cancelled = true;
}
private _checkCancellation(){
if(this._cancelled){
this._cancelled = false;
this._running = false;
throw new Error("TRAJECTORY FINDER CANCELLED");
}
}
public async searchOptimalTrajectory(sequence: FlybySequence, startDateMin: number, startDateMax: number, depAltitude: number){
this._running = true;
this._initPlot();
this.plot.clearPlot()
this._calculatePopulationSize(sequence);
this._calculatePopulationChunks();
@@ -57,15 +49,18 @@ export class TrajectorySolver {
await this._createStartPopulation();
this._updatePlot(0);
this._checkCancellation();
const {maxGenerations} = this.config.trajectorySearch;
for(let i = 0; i < maxGenerations; i++) {
if(this._cancelled){
this._cancelled = false;
this._running = false;
throw new Error("TRAJECTORY FINDER CANCELLED");
}
await this._generateNextPopulation();
this._updatePlot(1 + i);
this._checkCancellation();
}
this._running = false;
}
@@ -73,12 +68,6 @@ export class TrajectorySolver {
return this._workerPool.passData({depAltitude, sequence, startDateMin, startDateMax});
}
private async _createStartPopulation(){
const inputs = this._firstGenerationInputs();
const results = await this._workerPool.runPool<GenerationResult>(inputs);
this._mergeResultsChunks(results);
}
private _calculatePopulationChunks(){
const {splitLimit} = this.config.trajectorySearch;
const numChunks = this._workerPool.optimizeUsedWorkersCount(this.popSize, splitLimit);
@@ -97,17 +86,38 @@ export class TrajectorySolver {
this._chunkIndices = chunkIndices;
}
private _calculatePopulationSize(sequence: FlybySequence){
const {popSizeDimScale} = this.config.trajectorySearch;
this.popSize = popSizeDimScale * (4 * (sequence.length - 1) + 2);
private async _createStartPopulation(){
const inputs = [];
for(let i = 0; i < this._numChunks; i++) {
const start = this._chunkIndices[i * 2];
const end = this._chunkIndices[i * 2 + 1];
inputs.push({
start: true,
chunkStart: start,
chunkEnd: end
});
}
const results = await this._workerPool.runPool<GenerationResult>(inputs);
this._mergeResultsChunks(results);
}
private async _generateNextPopulation(){
const inputs = this._nextGenerationInputs();
const inputs = [];
for(let i = 0; i < this._numChunks; i++) {
inputs[i] = {
population: this._population,
deltaVs: this._deltaVs,
};
}
const results = await this._workerPool.runPool<GenerationResult>(inputs);
this._mergeResultsChunks(results);
}
private _calculatePopulationSize(sequence: FlybySequence){
const {popSizeDimScale} = this.config.trajectorySearch;
this.popSize = popSizeDimScale * sequence.length;
}
private _mergeResultsChunks(results: GenerationResult[]){
const popChunks: Agent[][] = [];
const dVChunks: number[][] = [];
@@ -127,49 +137,4 @@ export class TrajectorySolver {
this.bestTrajectorySteps = bestSteps;
this.bestDeltaV = bestDeltaV;
}
private _firstGenerationInputs() {
const inputs = [];
for(let i = 0; i < this._numChunks; i++) {
const {start, end} = this._chunkStartEnd(i);
inputs.push({
start: true,
chunkStart: start,
chunkEnd: end
});
}
return inputs;
}
private _nextGenerationInputs() {
const inputs = [];
for(let i = 0; i < this._numChunks; i++) {
const {start, end} = this._chunkStartEnd(i);
inputs[i] = {
population: this._population,
deltaVs: this._deltaVs,
chunkStart: start,
chunkEnd: end
};
}
return inputs;
}
private _chunkStartEnd(index: number){
return {
start: this._chunkIndices[index * 2],
end: this._chunkIndices[index * 2 + 1]
};
}
private get _bestMeanDeltaV() {
let mean = 0;
for(const dv of this._deltaVs)
mean += dv;
mean /= this.popSize;
return {
mean: mean,
best: this.bestDeltaV
};
}
}

22
src/main/solvers/trajectory.ts
View File

@@ -5,6 +5,7 @@ import { createOrbitPoints, createLine, createSprite } from "../utilities/geomet
import { Orbit } from "../objects/orbit.js";
import { SolarSystem } from "../objects/system.js";
import { KSPTime } from "../utilities/time.js";
import { CameraController } from "../objects/camera.js";
export class Trajectory {
public readonly orbits: Orbit[] = [];
@@ -45,8 +46,8 @@ export class Trajectory {
for(let i = 0; i < this.orbits.length; i++) {
const orbit = this.orbits[i];
const {beginAngle, endAngle} = this.steps[i];
const orbitPoints = createOrbitPoints(orbit, samplePoints, scale, beginAngle, endAngle);
const {begin, end} = this.steps[i].drawAngles;
const orbitPoints = createOrbitPoints(orbit, samplePoints, scale, begin, end);
const color = new THREE.Color(`hsl(${i*35 % 360}, 100%, 85%)`);
const orbitLine = createLine(orbitPoints, resolution, {
color: color.getHex(),
@@ -65,7 +66,7 @@ export class Trajectory {
if(step.maneuvre){
const group = this.system.objectsOfBody(step.attractorId);
const sprite = createSprite(Trajectory.arrowMaterial, 0xFFFFFF, false, maneuvreArrowSize);
const {x, y, z} = step.maneuvre.manoeuvrePosition;
const {x, y, z} = step.maneuvre.position;
sprite.position.set(x, y, z);
sprite.position.multiplyScalar(scale);
group.add(sprite);
@@ -109,16 +110,17 @@ export class Trajectory {
}
}
public fillResultControls(maneuvreSelector: Selector, resultSpans: ResultPanelSpans, stepSlider: DiscreteRange, systemTime: TimeSelector){
public fillResultControls(maneuvreSelector: Selector, resultSpans: ResultPanelSpans, stepSlider: DiscreteRange, systemTime: TimeSelector, controls: CameraController){
const depDate = new KSPTime(this.steps[0].dateOfStart, this.config.time);
resultSpans.totalDVSpan.innerHTML = this._totalDeltaV.toFixed(1);
resultSpans.depDateSpan.innerHTML = depDate.stringYDHMS("hms", "ut");
resultSpans.depDateSpan.onclick = () => {
this.system.date = depDate.dateSeconds;
controls.centerOnTarget();
systemTime.time.dateSeconds = depDate.dateSeconds;
systemTime.update();
systemTime.onChange();
};
stepSlider.setMinMax(0, this.steps.length - 1);
@@ -158,11 +160,17 @@ export class Trajectory {
resultSpans.maneuvreNumber.innerHTML = (index + 1).toString();
resultSpans.dateSpan.onclick = () => {
systemTime.time.dateSeconds = depDate.dateSeconds + dateEMT.dateSeconds;
const date = depDate.dateSeconds + dateEMT.dateSeconds;
this.system.date = date;
controls.centerOnTarget();
systemTime.time.dateSeconds = date;
systemTime.update();
systemTime.onChange();
};
});
for(const step of this.steps){
console.log(step);
}
}
private _displayStepsUpTo(index: number){

62
src/types.d.ts vendored
View File

@@ -97,6 +97,7 @@ interface TrajectorySearchSettings {
readonly minLegDuration: number;
readonly popSizeDimScale: number;
readonly maxGenerations: number;
readonly fbRadiusMaxScale: number;
}
interface Config {
@@ -171,7 +172,15 @@ type OrbitalState3D = {
vel: Vector3,
};
type OrbitalElements = {
type OrbitalElements2D = {
readonly eccentricity: number,
readonly periapsisVec: Vector2,
readonly semiMajorAxis: number,
readonly orbitalParam: number,
readonly clockwise: boolean
}
type OrbitalElements3D = {
readonly eccentricity: number,
readonly periapsisDir: Vector3,
readonly semiMajorAxis: number,
@@ -182,29 +191,23 @@ type OrbitalElements = {
readonly orbitalParam: number
};
type OrbitalElements2D = {
readonly eccentricity: number,
readonly periapsisDir: Vector2,
readonly semiMajorAxis: number,
readonly orbitalParam: number,
readonly clockwise: boolean
}
type ArcEndsAngles = {begin: number, end: number};
type TrajectoryStep = {
orbitElts: OrbitalElements,
orbitElts: OrbitalElements3D,
attractorId: number,
beginAngle: number,
endAngle: number,
angles: ArcEndsAngles,
drawAngles: ArcEndsAngles,
dateOfStart: number,
duration: number,
maneuvre?: ManeuvreInfo
};
type ManeuvreInfo = {
deltaVToPrevStep: Vector3,
progradeDir: Vector3,
manoeuvrePosition: Vector3,
context: ManeuvreContext
position: Vector3,
deltaVToPrevStep: Vector3,
progradeDir: Vector3,
context: ManeuvreContext
}
type ManeuvreContext =
@@ -213,13 +216,34 @@ type ManeuvreContext =
| {type: "circularization"}
;
type DepartureInfo = {
phaseParam: number,
ejVelParam: number,
dateParam: number,
};
type LegInfo = {
duration: number,
dsmOffset: number,
theta: number,
phi: number
exitedBodyId: number,
targetBodyId: number,
durationParam: number,
duration: number,
dsmParam: number,
}
type FlybyInfo = {
flybyBodyId: number,
normAngleParam: number,
periRadiParam: number
};
type SecondArcData = {
preDSMState: OrbitalState3D,
fbBodyId: number,
fbBodyState: OrbitalState3D,
flybyOrbit: OrbitalElements3D,
soiEnterAngle: number,
};
type GenerationResult = {
bestSteps: TrajectoryStep[],
bestDeltaV: number,