martin/KerbalCalculations
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Simplified interpolation

Browse Source
This commit is contained in:
Martin Asprusten 2026-04-04 11:17:30 +02:00
parent dd63c64443
commit b3366eba01
4 changed files with 348 additions and 208 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
index.html
View File

@ -1,5 +1,5 @@
<!doctype html> <!doctype html>
<html lang="en"> <html lang="en-US">
<head> <head>
<meta charset="UTF-8" /> <meta charset="UTF-8" />
<link rel="icon" type="image/png" href="/satellite.png" /> <link rel="icon" type="image/png" href="/satellite.png" />

43
src/calculations/mathematics.ts
View File

@ -133,6 +133,20 @@ export function addVector(vectorOne: number[][], vectorTwo: number[][]): number[
return result; return result;
} }
export function addVectors(vectors: number[][][]) {
let resultVector = [
[0],
[0],
[0]
];
vectors.forEach(vector => {
resultVector = addVector(resultVector, vector);
});
return resultVector;
}
export function addMatrix(matrixOne: number[][], matrixTwo: number[][]): number[][] { export function addMatrix(matrixOne: number[][], matrixTwo: number[][]): number[][] {
if (!checkIfValidMatrix(matrixOne) || !checkIfValidMatrix(matrixTwo)) { if (!checkIfValidMatrix(matrixOne) || !checkIfValidMatrix(matrixTwo)) {
throw new TypeError("Two valid matrices are required"); throw new TypeError("Two valid matrices are required");
@ -227,3 +241,32 @@ export function multiplyMatrixWithScalar(scalar: number, matrix: number[][]): nu
return result; return result;
} }
export function invertTwoByTwoMatrix(matrix: number[][]): number[][] {
if (!checkIfValidMatrix(matrix)) {
throw new TypeError("A valid matrix is required");
}
if (matrix.length != 2) {
throw new TypeError("Can only invert 2x2 matrices");
}
if (matrix[0].length != 2) {
throw new TypeError("Can only invert 2x2 matrices");
}
let divisor = matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
let multiplier = 1 / divisor;
return [
[
matrix[1][1] * multiplier,
-matrix[0][1] * multiplier
],
[
-matrix[1][0] * multiplier,
matrix[0][0] * multiplier
]
];
}

496
src/calculations/orbit-calculations.ts
View File

@ -1,6 +1,6 @@
import type { InterpolationParameters } from "../gui/interpolate"; import type { InterpolationParameters } from "../gui/interpolate";
import type { Body } from "./constants"; import type { Body } from "./constants";
import { addMatrix, addVector, getVectorMagnitude, matrixMultiply, multiplyMatrixWithScalar, normalizeVector, subtractVector, vectorCrossProduct, vectorDotProduct } from "./mathematics"; import { addMatrix, addVector, getVectorMagnitude, invertTwoByTwoMatrix, matrixMultiply, multiplyMatrixWithScalar, normalizeVector, subtractVector, vectorCrossProduct, vectorDotProduct } from "./mathematics";
export interface Orbit { export interface Orbit {
semiLatusRectum: number, semiLatusRectum: number,
@ -362,7 +362,14 @@ export function getEccentricAndTrueAnomalyFromMeanAnomaly(meanAnomaly: number, e
} }
export function getOrbitalCoordinates(timeToPeriapsis: number, orbit: Orbit, planet: Body): OrbitalCoordinates { export function getOrbitalCoordinates(timeToPeriapsis: number, orbit: Orbit, planet: Body): OrbitalCoordinates {
const meanAngularMotion = Math.sqrt(planet.gravitationalParameter * Math.abs(1 - orbit.eccentricity**2)**3 / orbit.semiLatusRectum**3); let meanAngularMotion;
if (Math.abs(orbit.eccentricity - 1) < 0.0001) {
meanAngularMotion = Math.sqrt(8 * planet.gravitationalParameter / orbit.semiLatusRectum**3);
} else {
meanAngularMotion = Math.sqrt(planet.gravitationalParameter * Math.abs(1 - orbit.eccentricity**2)**3 / orbit.semiLatusRectum**3);
}
const meanAnomaly = meanAngularMotion * -timeToPeriapsis; const meanAnomaly = meanAngularMotion * -timeToPeriapsis;
@ -412,9 +419,11 @@ export function getOrbitalCoordinatesFromAltitude(altitude: number, headingInwar
let eccentricAnomaly; let eccentricAnomaly;
let meanAnomaly; let meanAnomaly;
let meanAngularMotion: number | null = null;
if (Math.abs(orbit.eccentricity - 1) < 0.0001) { if (Math.abs(orbit.eccentricity - 1) < 0.0001) {
eccentricAnomaly = trueAnomaly; eccentricAnomaly = trueAnomaly;
meanAnomaly = eccentricAnomaly; meanAnomaly = Math.sqrt(2) * (Math.tan(trueAnomaly / 2) + Math.tan(trueAnomaly / 2)**3 / 3);
meanAngularMotion = Math.sqrt(8 * planet.gravitationalParameter / orbit.semiLatusRectum**3);
} else if (orbit.eccentricity < 1) { } else if (orbit.eccentricity < 1) {
eccentricAnomaly = Math.atan2(Math.sqrt(1 - orbit.eccentricity**2)*Math.sin(trueAnomaly), orbit.eccentricity + Math.cos(trueAnomaly)); eccentricAnomaly = Math.atan2(Math.sqrt(1 - orbit.eccentricity**2)*Math.sin(trueAnomaly), orbit.eccentricity + Math.cos(trueAnomaly));
meanAnomaly = eccentricAnomaly - orbit.eccentricity * Math.sin(eccentricAnomaly); meanAnomaly = eccentricAnomaly - orbit.eccentricity * Math.sin(eccentricAnomaly);
@ -423,7 +432,9 @@ export function getOrbitalCoordinatesFromAltitude(altitude: number, headingInwar
meanAnomaly = orbit.eccentricity * Math.sinh(eccentricAnomaly) - eccentricAnomaly; meanAnomaly = orbit.eccentricity * Math.sinh(eccentricAnomaly) - eccentricAnomaly;
} }
const meanAngularMotion = Math.sqrt(planet.gravitationalParameter * Math.abs(1 - orbit.eccentricity**2)**3 / orbit.semiLatusRectum**3); if (meanAngularMotion == null) {
meanAngularMotion = Math.sqrt(planet.gravitationalParameter * Math.abs(1 - orbit.eccentricity**2)**3 / orbit.semiLatusRectum**3);
}
return { return {
orbit: orbit, orbit: orbit,
@ -977,6 +988,10 @@ export function findCheapestIntercept(startingSituation: OrbitalCoordinates, tar
let startOrbitMeanAnomaly = startingSituation.meanAnomaly + startTime * startingSituation.meanAngularMotion; let startOrbitMeanAnomaly = startingSituation.meanAnomaly + startTime * startingSituation.meanAngularMotion;
let endOrbitMeanAnomaly = targetSituation.meanAnomaly + endTime * targetSituation.meanAngularMotion; let endOrbitMeanAnomaly = targetSituation.meanAnomaly + endTime * targetSituation.meanAngularMotion;
if (Math.abs(startingSituation.orbit.eccentricity - 1) < 0.0001) {
}
let extraStartOrbits = Math.floor(startOrbitMeanAnomaly / (2 * Math.PI)); let extraStartOrbits = Math.floor(startOrbitMeanAnomaly / (2 * Math.PI));
let extraEndOrbits = Math.floor(endOrbitMeanAnomaly / (2 * Math.PI)); let extraEndOrbits = Math.floor(endOrbitMeanAnomaly / (2 * Math.PI));
@ -1141,7 +1156,7 @@ export function findCheapestLanding(startingCoordinates: OrbitalCoordinates, sta
let position = startingPosition; let position = startingPosition;
let mass = ship.mass; let mass = ship.mass;
let massUse = ship.thrust / ship.specificImpulse; let massUse = ship.thrust / (ship.specificImpulse * 9.81);
let timeSteps = 0; let timeSteps = 0;
let lengthOfStep = 0.01; let lengthOfStep = 0.01;
@ -1369,261 +1384,326 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
return trueAnomaly; return trueAnomaly;
} }
const ownFirstTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards);
const ownSecondTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards);
const getPossibleTargetVectors = (ownOrbit: Orbit, ownAltitude: number, targetAltitude: number, distanceToTarget: number): [number[][], number[][], number[][], number, number] => { const targetFirstTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
let ownTrueAnomaly = getTrueAnomalyFromAltitude(ownAltitude, ownOrbit.semiLatusRectum, ownOrbit.eccentricity, ownShipHeadedInwards); const targetSecondTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
let ownLocalX = ownAltitude * Math.cos(ownTrueAnomaly); const getVectors = (trueAnomaly: number, orbit: Orbit) => {
let ownLocalY = ownAltitude * Math.sin(ownTrueAnomaly); const radius = orbit.semiLatusRectum / (1 + orbit.eccentricity * Math.cos(trueAnomaly));
const localX = radius * Math.cos(trueAnomaly);
const localY = radius * Math.sin(trueAnomaly);
let ownDirection = normalizeVector( const vectorPointingTowardsShip = normalizeVector(
addVector( addVector(
multiplyMatrixWithScalar(ownLocalX, ownOrbit.coordinateAxes[0]), multiplyMatrixWithScalar(localX, orbit.coordinateAxes[0]),
multiplyMatrixWithScalar(ownLocalY, ownOrbit.coordinateAxes[1]) multiplyMatrixWithScalar(localY, orbit.coordinateAxes[1])
) )
); );
let perpendicularInPlane = normalizeVector(vectorCrossProduct(ownOrbit.coordinateAxes[2], ownDirection)); const perpendicularVectorInPlane = normalizeVector(
let perpendicularOutOfPlane = normalizeVector(vectorCrossProduct(ownDirection, perpendicularInPlane)); vectorCrossProduct(orbit.coordinateAxes[2], vectorPointingTowardsShip)
// Find phase angle to target
let phaseAngle = Math.acos((ownAltitude**2 + targetAltitude**2 - distanceToTarget**2) / (2 * ownAltitude * targetAltitude));
let distanceAlongDirection = targetAltitude * Math.cos(phaseAngle);
let distanceAlongNormal = targetAltitude * Math.sin(phaseAngle);
let root = multiplyMatrixWithScalar(distanceAlongDirection, ownDirection);
let vectorOne = multiplyMatrixWithScalar(distanceAlongNormal, perpendicularInPlane);
let vectorTwo = multiplyMatrixWithScalar(distanceAlongNormal, perpendicularOutOfPlane);
return [root, vectorOne, vectorTwo, distanceAlongDirection, distanceAlongNormal];
}
let firstPositionPossibleVectors = getPossibleTargetVectors(ownCoordinates.orbit, interpolationParameters.firstOwnAltitude, interpolationParameters.firstTargetAltitude, interpolationParameters.firstDistance);
let secondPositionPossibleVectors = getPossibleTargetVectors(ownCoordinates.orbit, interpolationParameters.secondOwnAltitude, interpolationParameters.secondTargetAltitude, interpolationParameters.secondDistance);
let targetFirstTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
let targetSecondTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
let c1 = firstPositionPossibleVectors[0];
let p11 = firstPositionPossibleVectors[1];
let p12 = firstPositionPossibleVectors[2];
let c2 = secondPositionPossibleVectors[0];
let p21 = secondPositionPossibleVectors[1];
let p22 = secondPositionPossibleVectors[2];
// Scale everything to avoid numerical explosions
let scaleFactor = 1 / getVectorMagnitude(c2);
c1 = multiplyMatrixWithScalar(scaleFactor, c1);
p11 = multiplyMatrixWithScalar(scaleFactor, p11);
p12 = multiplyMatrixWithScalar(scaleFactor, p12);
c2 = multiplyMatrixWithScalar(scaleFactor, c2);
p21 = multiplyMatrixWithScalar(scaleFactor, p21);
p22 = multiplyMatrixWithScalar(scaleFactor, p22);
let scaledAltitudeOne = interpolationParameters.firstTargetAltitude * scaleFactor;
let scaledAltitudeTwo = interpolationParameters.secondTargetAltitude * scaleFactor;
const constantPart = scaledAltitudeOne * scaledAltitudeTwo * Math.cos(targetSecondTrueAnomaly - targetFirstTrueAnomaly);
const functionToMinimize = (angleOne: number, angleTwo: number) => {
let d1 = addVector(c1, addVector(
multiplyMatrixWithScalar(Math.cos(angleOne), p11),
multiplyMatrixWithScalar(Math.sin(angleOne), p12)
));
let d2 = addVector(c2, addVector(
multiplyMatrixWithScalar(Math.cos(angleTwo), p21),
multiplyMatrixWithScalar(Math.sin(angleTwo), p22)
));
return vectorDotProduct(d1, d2) - constantPart;
}
const partialDerivativeAngleOne = (angleOne: number, angleTwo: number) => {
let d2 = addVector(c2, addVector(
multiplyMatrixWithScalar(Math.cos(angleTwo), p21),
multiplyMatrixWithScalar(Math.sin(angleTwo), p22)
));
let d1Dx = addVector(
multiplyMatrixWithScalar(-Math.sin(angleOne), p11),
multiplyMatrixWithScalar(Math.cos(angleOne), p12)
); );
return vectorDotProduct(d1Dx, d2); const perpendicularVectorOutOfPlane = orbit.coordinateAxes[2];
return [vectorPointingTowardsShip, perpendicularVectorInPlane, perpendicularVectorOutOfPlane];
} }
const partialDerivativeAngleTwo = (angleOne: number, angleTwo: number) => { const firstVectors = getVectors(ownFirstTrueAnomaly, ownCoordinates.orbit);
let d1 = addVector(c1, addVector( const secondVectors = getVectors(ownSecondTrueAnomaly, ownCoordinates.orbit);
multiplyMatrixWithScalar(Math.cos(angleOne), p11),
multiplyMatrixWithScalar(Math.sin(angleOne), p12)
));
let d2Dy = addVector( const getDistances = (ownAltitude: number, targetAltitude: number, targetDistance: number) => {
multiplyMatrixWithScalar(-Math.sin(angleTwo), p21), const angleBetweenPlanetAndTarget = Math.acos((ownAltitude**2 + targetAltitude**2 - targetDistance**2) / (2*ownAltitude*targetAltitude));
multiplyMatrixWithScalar(Math.cos(angleTwo), p22) const distanceAlongDirection = Math.cos(angleBetweenPlanetAndTarget) * targetAltitude;
const distanceAlongPerpendicular = Math.sin(angleBetweenPlanetAndTarget) * targetAltitude;
return [distanceAlongDirection, distanceAlongPerpendicular];
}
const [firstDistanceAlongDirection, firstDistancePerpendicularToDirection] = getDistances(interpolationParameters.firstOwnAltitude, interpolationParameters.firstTargetAltitude, interpolationParameters.firstDistance);
const [secondDistanceAlongDirection, secondDistancePerpendicularToDirection] = getDistances(interpolationParameters.secondOwnAltitude, interpolationParameters.secondTargetAltitude, interpolationParameters.secondDistance);
const phaseAngleDifference = (ownSecondTrueAnomaly + interpolationParameters.secondPhaseAngle) - (ownFirstTrueAnomaly + interpolationParameters.firstPhaseAngle);
// Now, try some Newton's method to estimate the two position's the target has gone through
// To avoid the maths exploding, we'll scale everything down a bit
const scaleFactor = 1;
const c1 = multiplyMatrixWithScalar(firstDistanceAlongDirection * scaleFactor, firstVectors[0]);
const n11 = multiplyMatrixWithScalar(firstDistancePerpendicularToDirection * scaleFactor, firstVectors[1]);
const n12 = multiplyMatrixWithScalar(firstDistancePerpendicularToDirection * scaleFactor, firstVectors[2]);
const c2 = multiplyMatrixWithScalar(secondDistanceAlongDirection * scaleFactor, secondVectors[0]);
const n21 = multiplyMatrixWithScalar(secondDistancePerpendicularToDirection * scaleFactor, secondVectors[1]);
const n22 = multiplyMatrixWithScalar(secondDistancePerpendicularToDirection * scaleFactor, secondVectors[2]);
let expectedDistance = Math.sqrt(interpolationParameters.firstTargetAltitude**2 + interpolationParameters.secondTargetAltitude**2 - 2 * interpolationParameters.firstTargetAltitude * interpolationParameters.secondTargetAltitude * Math.cos(targetSecondTrueAnomaly - targetFirstTrueAnomaly));
// Scale the expected distance too
expectedDistance *= scaleFactor;
const distanceFunction = (angleOne: number, angleTwo: number) => {
let p1 = addVector(
c1,
addVector(
multiplyMatrixWithScalar(Math.cos(angleOne), n11),
multiplyMatrixWithScalar(Math.sin(angleOne), n12)
)
); );
return vectorDotProduct(d1, d2Dy); let p2 = addVector(
c2,
addVector(
multiplyMatrixWithScalar(Math.cos(angleTwo), n21),
multiplyMatrixWithScalar(Math.sin(angleTwo), n22)
)
);
return getVectorMagnitude(addVector(p2, multiplyMatrixWithScalar(-1, p1))) - expectedDistance;
} }
let firstOwnTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards); const distancePartialDerivativeAngleOne = (angleOne: number, angleTwo: number) => {
let secondOwnTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards); let distance = distanceFunction(angleOne, angleTwo) + expectedDistance;
return (
let horizontalAngleDifference = (secondOwnTrueAnomaly - interpolationParameters.secondPhaseAngle) - (firstOwnTrueAnomaly - interpolationParameters.firstPhaseAngle); Math.sin(angleOne) * vectorDotProduct(n11, c2)
+ Math.sin(angleOne) * Math.cos(angleTwo) * vectorDotProduct(n21, n11)
const phaseAngleFunction = (angleOne: number, angleTwo: number): number => { - Math.sin(angleOne) * Math.cos(angleOne) * vectorDotProduct(n11, n11)
let d1 = addVector(c1, multiplyMatrixWithScalar(Math.cos(angleOne), p11)); - Math.cos(angleOne) * vectorDotProduct(n12, c2)
let d2 = addVector(c2, multiplyMatrixWithScalar(Math.cos(angleTwo), p21)); - Math.cos(angleOne) * Math.sin(angleTwo) * vectorDotProduct(n12, n22)
return vectorDotProduct(d1, d2) - Math.cos(horizontalAngleDifference) * getVectorMagnitude(d1) * getVectorMagnitude(d2); + Math.sin(angleOne) * Math.sin(angleOne) * vectorDotProduct(n12, n12)
) / distance;
} }
const phaseAnglePartialDerivativeAngleOne = (angleOne: number, angleTwo: number) => { const distancePartialDerivativeAngleTwo = (angleOne: number, angleTwo: number) => {
let d1 = addVector(c1, multiplyMatrixWithScalar(Math.cos(angleOne), p11)); let distance = distanceFunction(angleOne, angleTwo) + expectedDistance;
let d1Dx = multiplyMatrixWithScalar(-Math.sin(angleOne), p11); return (
let d2 = addVector(c2, multiplyMatrixWithScalar(Math.cos(angleTwo), p21)); Math.cos(angleTwo) * Math.sin(angleTwo) * vectorDotProduct(n22, n22)
- Math.cos(angleTwo) * vectorDotProduct(n22, c1)
return vectorDotProduct(d1Dx, d2) - Math.cos(horizontalAngleDifference) * getVectorMagnitude(d2) * Math.sin(angleOne) * Math.cos(angleOne) * firstPositionPossibleVectors[4]**2 / getVectorMagnitude(d1); - Math.cos(angleTwo) * Math.sin(angleOne) * vectorDotProduct(n22, n12)
- Math.sin(angleTwo) * Math.cos(angleTwo) * vectorDotProduct(n21, n21)
+ Math.sin(angleTwo) * vectorDotProduct(n21, c1)
+ Math.sin(angleTwo) * Math.cos(angleOne) * vectorDotProduct(n21, n11)
) / distance;
} }
const planeAngleFunction = (angleOne: number, angleTwo: number) => {
const horizontalOne = addVector(
c1,
multiplyMatrixWithScalar(Math.cos(angleOne), n11)
);
const gradientDescent = (startingGuessAngleOne: number, startingGuessAngleTwo: number): [number, number, number] => { const horizontalTwo = addVector(
let angleOne = startingGuessAngleOne; c2,
let angleTwo = startingGuessAngleTwo; multiplyMatrixWithScalar(Math.cos(angleTwo), n21)
);
let distanceAlongDirection = getVectorMagnitude(firstPositionPossibleVectors[0]); return vectorDotProduct(horizontalOne, horizontalTwo) / (getVectorMagnitude(horizontalOne) * getVectorMagnitude(horizontalTwo)) - Math.cos(phaseAngleDifference);
let distanceAlongPerpendicular = getVectorMagnitude(firstPositionPossibleVectors[1]);
let value = minimizeFunction(angleOne, angleTwo);
// Do up to 100 minimizing steps
for (var i = 0; i < 100; i++) {
let dfDx = partialDerivativeAngleOne(angleOne, angleTwo);
let dfDy = partialDerivativeAngleTwo(angleOne, angleTwo);
let phaseAngleValue = distanceAlongPerpendicular * Math.cos(angleOne) / distanceAlongDirection - Math.tan(interpolationParameters.firstPhaseAngle);
let phaseAngleDerivative = -distanceAlongPerpendicular * Math.sin(angleOne) / distanceAlongDirection;
let changeAngleOne = dfDx * phaseAngleValue / (dfDy * phaseAngleDerivative);
let changeAngleTwo = (dfDx * phaseAngleValue - phaseAngleDerivative * value) / (dfDy * phaseAngleDerivative);
let trialAngleOne = angleOne;
let trialAngleTwo = angleTwo;
let trialValue = value;
let divisor = 1;
// We don't want things to blow up
while (Math.abs(changeAngleOne) / 2**divisor > Math.PI / 1000 || Math.abs(changeAngleTwo) / 2**divisor > Math.PI / 1000) {
divisor +=1;
} }
const planeAnglePartialDerivativeAngleOne = (angleOne: number, angleTwo: number) => {
const horizontalOne = addVector(
c1,
multiplyMatrixWithScalar(Math.cos(angleOne), n11)
);
const horizontalTwo = addVector(
c2,
multiplyMatrixWithScalar(Math.cos(angleTwo), n21)
);
return (
-Math.sin(angleOne) * vectorDotProduct(n11, c2)
-Math.sin(angleOne) * Math.cos(angleTwo) * vectorDotProduct(n11, n21)
-Math.cos(angleOne) * Math.sin(angleOne) * vectorDotProduct(n11, n11) * vectorDotProduct(horizontalOne, horizontalTwo) / getVectorMagnitude(horizontalOne)
) / (getVectorMagnitude(horizontalOne) * getVectorMagnitude(horizontalTwo));
}
const planeAnglePartialDerivativeAngleTwo = (angleOne: number, angleTwo: number) => {
const horizontalOne = addVector(
c1,
multiplyMatrixWithScalar(Math.cos(angleOne), n11)
);
const horizontalTwo = addVector(
c2,
multiplyMatrixWithScalar(Math.cos(angleTwo), n21)
);
return (
-Math.sin(angleTwo) * vectorDotProduct(n21, c1)
-Math.sin(angleTwo) * Math.cos(angleOne) * vectorDotProduct(n21, n11)
-Math.cos(angleTwo) * Math.sin(angleTwo) * vectorDotProduct(n21, n21) * vectorDotProduct(horizontalOne, horizontalTwo) / getVectorMagnitude(horizontalTwo)
) / (getVectorMagnitude(horizontalOne) * getVectorMagnitude(horizontalTwo))
}
// Try to Newton's method up in this bitch
const getFunctionVector = (anglesVector: number[][]) => {
return [
[distanceFunction(anglesVector[0][0], anglesVector[1][0])],
[planeAngleFunction(anglesVector[0][0], anglesVector[1][0])]
]
};
const getJacobianMatrix = (anglesVector: number[][]) => {
let angleOne = anglesVector[0][0];
let angleTwo = anglesVector[1][0];
return[
[distancePartialDerivativeAngleOne(angleOne, angleTwo), distancePartialDerivativeAngleTwo(angleOne, angleTwo)],
[planeAnglePartialDerivativeAngleOne(angleOne, angleTwo), planeAnglePartialDerivativeAngleTwo(angleOne, angleTwo)]
]
};
const performNewtonMethod = (initialGuess: number[][], iterations: number) => {
let anglesVector = [
[initialGuess[0][0]],
[initialGuess[1][0]]
];
for (var i = 0; i < iterations; i++) {
let functionVector = getFunctionVector(anglesVector);
let jacobianMatrix = getJacobianMatrix(anglesVector);
let update = matrixMultiply(invertTwoByTwoMatrix(jacobianMatrix), functionVector);
// Don't make updates too big
while (getVectorMagnitude(update) > Math.PI / 100) {
update = multiplyMatrixWithScalar(0.5, update);
}
/**
// Also don't update if it brings us further away from where we want to be
let trialFunctionVector = getFunctionVector(addVector(anglesVector, multiplyMatrixWithScalar(-1, update)));
let counter = 0;
let deadEnd = false; let deadEnd = false;
while (Math.abs(trialValue) >= Math.abs(value)) { while (Math.abs(trialFunctionVector[0][0]) >= Math.abs(functionVector[0][0])) {
trialAngleOne = angleOne - changeAngleOne / (2 ** divisor); counter += 1;
trialAngleTwo = angleTwo - changeAngleTwo / (2 ** divisor); if (counter >= 16) {
trialValue = minimizeFunction(trialAngleOne, trialAngleTwo);
divisor += 1;
if (divisor >= 32) {
deadEnd = true; deadEnd = true;
break; break;
} }
update = multiplyMatrixWithScalar(0.5, update);
trialFunctionVector = getFunctionVector(addVector(anglesVector, multiplyMatrixWithScalar(-1, update)));
} }
if (deadEnd) { if (deadEnd) {
return [angleOne, angleTwo, value]; break;
}*/
anglesVector = addVector(anglesVector, multiplyMatrixWithScalar(-1, update));
};
return anglesVector;
} }
value = trialValue; const getAnglesFromPhaseAngles = (firstPhaseAngle: number, secondPhaseAngle: number) => {
angleOne = trialAngleOne; let angleOne;
angleTwo = trialAngleTwo; let angleTwo;
if (Math.abs(interpolationParameters.firstPhaseAngle) < Math.PI / 2) {
angleOne = Math.acos(Math.tan(firstPhaseAngle) * firstDistanceAlongDirection / firstDistancePerpendicularToDirection);
} else {
angleOne = Math.acos(Math.tan(Math.PI - firstPhaseAngle) * -firstDistanceAlongDirection / firstDistancePerpendicularToDirection);
} }
return [angleOne, angleTwo, value]; if (Math.abs(interpolationParameters.secondPhaseAngle) < Math.PI / 2) {
angleTwo = Math.acos(Math.tan(secondPhaseAngle) * secondDistanceAlongDirection / secondDistancePerpendicularToDirection);
} else {
angleTwo = Math.acos(Math.tan(Math.PI - secondPhaseAngle) * -secondDistanceAlongDirection / secondDistancePerpendicularToDirection);
} }
let firstCenterLength = getVectorMagnitude(firstPositionPossibleVectors[0]); return [angleOne, angleTwo];
let firstPerpendicularLength = getVectorMagnitude(firstPositionPossibleVectors[1]);
let secondCenterLength = getVectorMagnitude(secondPositionPossibleVectors[0]);
let secondPerpendicularLength = getVectorMagnitude(secondPositionPossibleVectors[1]);
let angleOneGuess = Math.acos(firstCenterLength * Math.tan(interpolationParameters.firstPhaseAngle) / firstPerpendicularLength);
let angleTwoGuess = Math.acos(secondCenterLength * Math.tan(interpolationParameters.secondPhaseAngle) / secondPerpendicularLength);
let bestAngleOne = 0;
let bestAngleTwo = 0;
let bestValue: number | null = null;
[angleOneGuess, -angleOneGuess].forEach(angleOne => {
[angleTwoGuess, -angleTwoGuess].forEach(angleTwo => {
let value = minimizeFunction(angleOne, angleTwo);
if (bestValue == null || value < bestValue) {
bestValue = value;
bestAngleOne = angleOne;
bestAngleTwo = angleTwo;
} }
});
});
[bestAngleOne, bestAngleOne, bestValue] = gradientDescent(bestAngleOne, bestAngleTwo); let bestAngles = [[0], [0]];
let bestDistance: number | null = null;
let bestFirstPosition = addVector( for (var i = -100; i < 101; i++) {
firstPositionPossibleVectors[0], for (var j = -100; j < 101; j++) {
let firstPhaseAngle = interpolationParameters.firstPhaseAngle + i * Math.PI / 18000;
let secondPhaseAngle = interpolationParameters.secondPhaseAngle + j * Math.PI / 18000;
let [angleOne, angleTwo] = getAnglesFromPhaseAngles(firstPhaseAngle, secondPhaseAngle);
let distance = distanceFunction(angleOne, angleTwo);
if (bestDistance == null || Math.abs(distance) < Math.abs(bestDistance)) {
bestDistance = distance;
bestAngles = [[angleOne], [angleTwo]];
}
}
}
if (!interpolationParameters.targetAbovePlane) {
bestAngles = [[-bestAngles[0][0]], [-bestAngles[1][0]]];
}
// Do some additional Newtoning on the best angles we've found so far
//bestAngles = performNewtonMethod(bestAngles, 10000);
let targetPositionOne = addVector(multiplyMatrixWithScalar(firstDistanceAlongDirection, firstVectors[0]),
addVector( addVector(
multiplyMatrixWithScalar(Math.cos(bestAngleOne), firstPositionPossibleVectors[1]), multiplyMatrixWithScalar(Math.cos(bestAngles[0][0]) * firstDistancePerpendicularToDirection, firstVectors[1]),
multiplyMatrixWithScalar(Math.sin(bestAngleOne), firstPositionPossibleVectors[2]) multiplyMatrixWithScalar(Math.sin(bestAngles[0][0]) * firstDistancePerpendicularToDirection, firstVectors[2])
)
);
let bestSecondPosition = addVector(
secondPositionPossibleVectors[0],
addVector(
multiplyMatrixWithScalar(Math.cos(bestAngleTwo), secondPositionPossibleVectors[1]),
multiplyMatrixWithScalar(Math.sin(bestAngleTwo), secondPositionPossibleVectors[2])
) )
); );
if (bestFirstPosition == null || bestSecondPosition == null) { let targetPositionTwo = addVector(multiplyMatrixWithScalar(secondDistanceAlongDirection, secondVectors[0]),
return null; addVector(
} multiplyMatrixWithScalar(Math.cos(bestAngles[1][0]) * secondDistancePerpendicularToDirection, secondVectors[1]),
// We can now find the vector that defines the plane the target is in multiplyMatrixWithScalar(Math.sin(bestAngles[1][0]) * secondDistancePerpendicularToDirection, secondVectors[2])
let normalVector = normalizeVector(vectorCrossProduct(bestFirstPosition, bestSecondPosition)); )
let firstTargetPositionDirection = normalizeVector(bestFirstPosition); );
// Rotate the position vector about the normal vector by the true anomaly to find the vector pointing to the periapsis let normalVector = normalizeVector(vectorCrossProduct(targetPositionTwo, targetPositionOne));
let identityMatrix = [[1, 0, 0], [0, 1, 0], [0, 0, 1]];
let crossProductMatrix = [ // Rotate the position vector about this normal vector to get the direction of the periapsis
[0, -normalVector[2][0], normalVector[1][0]], let ux = normalVector[0][0];
[normalVector[2][0], 0, -normalVector[0][0]], let uy = normalVector[1][0];
[-normalVector[1][0], normalVector[0][0], 0] let uz = normalVector[2][0];
];
let outerProductMatrix = [ let rotationMatrix = [
[normalVector[0][0]**2, normalVector[0][0]*normalVector[1][0], normalVector[0][0]*normalVector[2][0]], [
[normalVector[0][0]*normalVector[1][0], normalVector[1][0]**2, normalVector[1][0]*normalVector[2][0]], ux*ux*(1 - Math.cos(-targetFirstTrueAnomaly)) + Math.cos(-targetFirstTrueAnomaly),
[normalVector[0][0]*normalVector[2][0], normalVector[1][0]*normalVector[2][0], normalVector[2][0]**2] ux*uy*(1 - Math.cos(-targetFirstTrueAnomaly)) - uz * Math.sin(-targetFirstTrueAnomaly),
ux*uz*(1 - Math.cos(-targetFirstTrueAnomaly)) + uy * Math.sin(-targetFirstTrueAnomaly)
],
[
ux*ux*(1 - Math.cos(-targetFirstTrueAnomaly)) + uz * Math.sin(-targetFirstTrueAnomaly),
uy*uy*(1 - Math.cos(-targetFirstTrueAnomaly)) + Math.cos(-targetFirstTrueAnomaly),
uy*uz*(1 - Math.cos(-targetFirstTrueAnomaly)) - ux * Math.sin(-targetFirstTrueAnomaly)
],
[
ux*uz*(1 - Math.cos(-targetFirstTrueAnomaly)) - uy * Math.sin(-targetFirstTrueAnomaly),
uy*uz*(1 - Math.cos(-targetFirstTrueAnomaly)) + ux * Math.sin(-targetFirstTrueAnomaly),
uz*uz*(1 - Math.cos(-targetFirstTrueAnomaly)) + Math.cos(-targetFirstTrueAnomaly)
]
]; ];
let rotationMatrix = addMatrix( let periapsisVector = normalizeVector(matrixMultiply(rotationMatrix, targetPositionOne));
multiplyMatrixWithScalar(Math.cos(-firstTargetTrueAnomaly), identityMatrix),
addMatrix(
multiplyMatrixWithScalar(Math.sin(-firstTargetTrueAnomaly), crossProductMatrix),
multiplyMatrixWithScalar(1 - Math.cos(-firstTargetTrueAnomaly), outerProductMatrix)
)
);
let targetOrbitXVector = matrixMultiply(rotationMatrix, firstTargetPositionDirection);
let targetOrbitYVector = vectorCrossProduct(normalVector, targetOrbitXVector);
let targetOrbit: Orbit = { let targetOrbit: Orbit = {
semiLatusRectum: targetSemiLatusRectum, semiLatusRectum: targetSemiLatusRectum,
eccentricity: targetEccentricity, eccentricity: targetEccentricity,
coordinateAxes: [targetOrbitXVector, targetOrbitYVector, normalVector] coordinateAxes: [
periapsisVector,
normalizeVector(vectorCrossProduct(normalVector, periapsisVector)),
normalVector
]
}; };
let timePassed = getTimeBetweenTrueAnomalies(ownCoordinates.trueAnomaly, getTrueAnomalyFromAltitude(interpolationParameters.secondOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards), ownCoordinates.orbit, body); let targetCoordinates: OrbitalCoordinates = getOrbitalCoordinatesFromAltitude(
interpolationParameters.secondTargetAltitude,
targetHeadedInwards,
targetOrbit,
body
);
return [getOrbitalCoordinatesFromAltitude(interpolationParameters.secondTargetAltitude, targetHeadedInwards, targetOrbit, body), timePassed]; let epochTrueAnomaly = ownCoordinates.trueAnomaly;
while (epochTrueAnomaly + 2 * Math.PI < ownFirstTrueAnomaly) {
epochTrueAnomaly += 2 * Math.PI;
}
const timeElapsed = getTimeBetweenTrueAnomalies(epochTrueAnomaly, ownSecondTrueAnomaly, ownCoordinates.orbit, body);
return [
targetCoordinates,
timeElapsed
]
} }

17
src/gui/interpolate.ts
View File

@ -7,6 +7,7 @@ export interface InterpolationParameters {
targetApoapsis: number, targetApoapsis: number,
targetSpeed: number, targetSpeed: number,
targetAltitude: number, targetAltitude: number,
targetAbovePlane: boolean,
firstOwnAltitude: number, firstOwnAltitude: number,
firstTargetAltitude: number, firstTargetAltitude: number,
firstDistance: number, firstDistance: number,
@ -23,6 +24,7 @@ export const DefaultInterpolationParameters: InterpolationParameters = {
targetApoapsis: 200000, targetApoapsis: 200000,
targetSpeed: 0, targetSpeed: 0,
targetAltitude: 0, targetAltitude: 0,
targetAbovePlane: true,
firstOwnAltitude: 0, firstOwnAltitude: 0,
firstTargetAltitude: 0, firstTargetAltitude: 0,
firstDistance: 0, firstDistance: 0,
@ -78,6 +80,13 @@ export class InterpolateOrbitGui {
let speedInput = createInputWithLabel("Target speed:", 0); let speedInput = createInputWithLabel("Target speed:", 0);
let altitudeInput = createInputWithLabel("Target altitude:", 0); let altitudeInput = createInputWithLabel("Target altitude:", 0);
let abovePlaneId = crypto.randomUUID();
let abovePlaneButton = document.createElement("input");
abovePlaneButton.setAttribute("type", "checkbox");
abovePlaneButton.setAttribute("id", abovePlaneId);
let abovePlaneLabel = createLabel(abovePlaneId, "Target is currently above orbital plane");
addToParent([abovePlaneButton, abovePlaneLabel]);
let instantOneHeader = document.createElement("h4"); let instantOneHeader = document.createElement("h4");
instantOneHeader.appendChild(document.createTextNode("Measurement one:")); instantOneHeader.appendChild(document.createTextNode("Measurement one:"));
this.parentDiv.appendChild(instantOneHeader); this.parentDiv.appendChild(instantOneHeader);
@ -102,6 +111,7 @@ export class InterpolateOrbitGui {
let apoapsis = parseFloat(apoapsisInput.value); let apoapsis = parseFloat(apoapsisInput.value);
let speed = parseFloat(speedInput.value); let speed = parseFloat(speedInput.value);
let altitude = parseFloat(altitudeInput.value); let altitude = parseFloat(altitudeInput.value);
let abovePlane = abovePlaneButton.checked;
let firstTargetAltitude = parseFloat(firstTargetAltitudeInput.value); let firstTargetAltitude = parseFloat(firstTargetAltitudeInput.value);
let firstOwnAltitude = parseFloat(firstOwnAltitudeInput.value); let firstOwnAltitude = parseFloat(firstOwnAltitudeInput.value);
let firstDistance = parseFloat(firstDistanceInput.value); let firstDistance = parseFloat(firstDistanceInput.value);
@ -122,6 +132,7 @@ export class InterpolateOrbitGui {
targetApoapsis: apoapsis, targetApoapsis: apoapsis,
targetSpeed: speed, targetSpeed: speed,
targetAltitude: altitude, targetAltitude: altitude,
targetAbovePlane: abovePlane,
firstTargetAltitude: firstTargetAltitude, firstTargetAltitude: firstTargetAltitude,
firstOwnAltitude: firstOwnAltitude, firstOwnAltitude: firstOwnAltitude,
firstDistance: firstDistance, firstDistance: firstDistance,
@ -140,6 +151,7 @@ export class InterpolateOrbitGui {
apoapsisInput, apoapsisInput,
speedInput, speedInput,
altitudeInput, altitudeInput,
abovePlaneButton,
firstOwnAltitudeInput, firstOwnAltitudeInput,
firstTargetAltitudeInput, firstTargetAltitudeInput,
firstDistanceInput, firstDistanceInput,
@ -171,6 +183,11 @@ export class InterpolateOrbitGui {
apoapsisInput.value = value.targetApoapsis.toString(); apoapsisInput.value = value.targetApoapsis.toString();
speedInput.value = value.targetSpeed.toString(); speedInput.value = value.targetSpeed.toString();
altitudeInput.value = value.targetAltitude.toString(); altitudeInput.value = value.targetAltitude.toString();
if (value.targetAbovePlane) {
abovePlaneButton.setAttribute("checked", "");
} else {
abovePlaneButton.removeAttribute("checked");
}
firstOwnAltitudeInput.value = value.firstOwnAltitude.toString(); firstOwnAltitudeInput.value = value.firstOwnAltitude.toString();
firstTargetAltitudeInput.value = value.firstTargetAltitude.toString(); firstTargetAltitudeInput.value = value.firstTargetAltitude.toString();
firstDistanceInput.value = value.firstDistance.toString(); firstDistanceInput.value = value.firstDistance.toString();