Optimize over phase angles

.gitignore

@@ -22,3 +22,4 @@ dist-ssr
 *.njsproj
 *.sln
 *.sw?
+*.sh

package-lock.json

@@ -7,11 +7,22 @@
     "": {
       "name": "kerbal-calculations",
       "version": "0.0.0",
+      "dependencies": {
+        "three": "^0.183.2"
+      },
       "devDependencies": {
+        "@types/three": "^0.183.1",
         "typescript": "^6.0.2",
         "vite": "^7.3.1"
       }
     },
+    "node_modules/@dimforge/rapier3d-compat": {
+      "version": "0.12.0",
+      "resolved": "https://registry.npmjs.org/@dimforge/rapier3d-compat/-/rapier3d-compat-0.12.0.tgz",
+      "integrity": "sha512-uekIGetywIgopfD97oDL5PfeezkFpNhwlzlaEYNOA0N6ghdsOvh/HYjSMek5Q2O1PYvRSDFcqFVJl4r4ZBwOow==",
+      "dev": true,
+      "license": "Apache-2.0"
+    },
     "node_modules/@esbuild/aix-ppc64": {
       "version": "0.27.3",
       "resolved": "https://registry.npmjs.org/@esbuild/aix-ppc64/-/aix-ppc64-0.27.3.tgz",
@@ -804,6 +815,13 @@
         "win32"
       ]
     },
+    "node_modules/@tweenjs/tween.js": {
+      "version": "23.1.3",
+      "resolved": "https://registry.npmjs.org/@tweenjs/tween.js/-/tween.js-23.1.3.tgz",
+      "integrity": "sha512-vJmvvwFxYuGnF2axRtPYocag6Clbb5YS7kLL+SO/TeVFzHqDIWrNKYtcsPMibjDx9O+bu+psAy9NKfWklassUA==",
+      "dev": true,
+      "license": "MIT"
+    },
     "node_modules/@types/estree": {
       "version": "1.0.8",
       "resolved": "https://registry.npmjs.org/@types/estree/-/estree-1.0.8.tgz",
@@ -811,6 +829,43 @@
       "dev": true,
       "license": "MIT"
     },
+    "node_modules/@types/stats.js": {
+      "version": "0.17.4",
+      "resolved": "https://registry.npmjs.org/@types/stats.js/-/stats.js-0.17.4.tgz",
+      "integrity": "sha512-jIBvWWShCvlBqBNIZt0KAshWpvSjhkwkEu4ZUcASoAvhmrgAUI2t1dXrjSL4xXVLB4FznPrIsX3nKXFl/Dt4vA==",
+      "dev": true,
+      "license": "MIT"
+    },
+    "node_modules/@types/three": {
+      "version": "0.183.1",
+      "resolved": "https://registry.npmjs.org/@types/three/-/three-0.183.1.tgz",
+      "integrity": "sha512-f2Pu5Hrepfgavttdye3PsH5RWyY/AvdZQwIVhrc4uNtvF7nOWJacQKcoVJn0S4f0yYbmAE6AR+ve7xDcuYtMGw==",
+      "dev": true,
+      "license": "MIT",
+      "dependencies": {
+        "@dimforge/rapier3d-compat": "~0.12.0",
+        "@tweenjs/tween.js": "~23.1.3",
+        "@types/stats.js": "*",
+        "@types/webxr": ">=0.5.17",
+        "@webgpu/types": "*",
+        "fflate": "~0.8.2",
+        "meshoptimizer": "~1.0.1"
+      }
+    },
+    "node_modules/@types/webxr": {
+      "version": "0.5.24",
+      "resolved": "https://registry.npmjs.org/@types/webxr/-/webxr-0.5.24.tgz",
+      "integrity": "sha512-h8fgEd/DpoS9CBrjEQXR+dIDraopAEfu4wYVNY2tEPwk60stPWhvZMf4Foo5FakuQ7HFZoa8WceaWFervK2Ovg==",
+      "dev": true,
+      "license": "MIT"
+    },
+    "node_modules/@webgpu/types": {
+      "version": "0.1.69",
+      "resolved": "https://registry.npmjs.org/@webgpu/types/-/types-0.1.69.tgz",
+      "integrity": "sha512-RPmm6kgRbI8e98zSD3RVACvnuktIja5+yLgDAkTmxLr90BEwdTXRQWNLF3ETTTyH/8mKhznZuN5AveXYFEsMGQ==",
+      "dev": true,
+      "license": "BSD-3-Clause"
+    },
     "node_modules/esbuild": {
       "version": "0.27.3",
       "resolved": "https://registry.npmjs.org/esbuild/-/esbuild-0.27.3.tgz",
@@ -871,6 +926,13 @@
         }
       }
     },
+    "node_modules/fflate": {
+      "version": "0.8.2",
+      "resolved": "https://registry.npmjs.org/fflate/-/fflate-0.8.2.tgz",
+      "integrity": "sha512-cPJU47OaAoCbg0pBvzsgpTPhmhqI5eJjh/JIu8tPj5q+T7iLvW/JAYUqmE7KOB4R1ZyEhzBaIQpQpardBF5z8A==",
+      "dev": true,
+      "license": "MIT"
+    },
     "node_modules/fsevents": {
       "version": "2.3.3",
       "resolved": "https://registry.npmjs.org/fsevents/-/fsevents-2.3.3.tgz",
@@ -886,6 +948,13 @@
         "node": "^8.16.0 || ^10.6.0 || >=11.0.0"
       }
     },
+    "node_modules/meshoptimizer": {
+      "version": "1.0.1",
+      "resolved": "https://registry.npmjs.org/meshoptimizer/-/meshoptimizer-1.0.1.tgz",
+      "integrity": "sha512-Vix+QlA1YYT3FwmBBZ+49cE5y/b+pRrcXKqGpS5ouh33d3lSp2PoTpCw19E0cKDFWalembrHnIaZetf27a+W2g==",
+      "dev": true,
+      "license": "MIT"
+    },
     "node_modules/nanoid": {
       "version": "3.3.11",
       "resolved": "https://registry.npmjs.org/nanoid/-/nanoid-3.3.11.tgz",
@@ -1009,6 +1078,12 @@
         "node": ">=0.10.0"
       }
     },
+    "node_modules/three": {
+      "version": "0.183.2",
+      "resolved": "https://registry.npmjs.org/three/-/three-0.183.2.tgz",
+      "integrity": "sha512-di3BsL2FEQ1PA7Hcvn4fyJOlxRRgFYBpMTcyOgkwJIaDOdJMebEFPA+t98EvjuljDx4hNulAGwF6KIjtwI5jgQ==",
+      "license": "MIT"
+    },
     "node_modules/tinyglobby": {
       "version": "0.2.15",
       "resolved": "https://registry.npmjs.org/tinyglobby/-/tinyglobby-0.2.15.tgz",

package.json

@@ -9,7 +9,11 @@
     "preview": "vite preview"
   },
   "devDependencies": {
+    "@types/three": "^0.183.1",
     "typescript": "^6.0.2",
     "vite": "^7.3.1"
+  },
+  "dependencies": {
+    "three": "^0.183.2"
   }
 }

src/calculations/orbit-calculations.ts

@@ -1,6 +1,6 @@
 import type { InterpolationParameters } from "../gui/interpolate";
 import type { Body } from "./constants";
-import { addMatrix, addVector, getVectorMagnitude, invertTwoByTwoMatrix, matrixMultiply, multiplyMatrixWithScalar, normalizeVector, subtractVector, vectorCrossProduct, vectorDotProduct } from "./mathematics";
+import { addVector, getVectorMagnitude, matrixMultiply, multiplyMatrixWithScalar, normalizeVector, subtractVector, vectorCrossProduct, vectorDotProduct } from "./mathematics";
 
 export interface Orbit {
 	semiLatusRectum: number,
@@ -192,13 +192,8 @@ export class LambertSolutions {
 			transferGoalTrueAnomaly -= 2 * Math.PI;
 		}
 
-		if (transferStartTrueAnomaly < 2 * Math.PI && transferGoalTrueAnomaly > 2 * Math.PI) {
-			transferStartTrueAnomaly -= 2 * Math.PI;
-			transferGoalTrueAnomaly -= 2 * Math.PI;
-		}
-
 		let closestPoint;
-		if (transferStartTrueAnomaly < 0 && transferGoalTrueAnomaly >= 0) {
+		if ((transferStartTrueAnomaly < 0 && transferGoalTrueAnomaly >= 0) || (transferStartTrueAnomaly < 2 * Math.PI && transferGoalTrueAnomaly >= 2 * Math.PI)) {
 			closestPoint = semiLatusRectum / (1 + eccentricity);
 		} else {
 			closestPoint = Math.min(this.positionOneMagnitude, this.positionTwoMagnitude);
@@ -570,7 +565,7 @@ export function getOrbitalPeriod(semiMajor: number, body: Body) {
 	return Math.sqrt(semiMajor**3 / body.gravitationalParameter);
 }
 
-export function perform2dGradientDescent(functionToMinimize: ((variableOne: number, variableTwo: number) => number | null), initialGuessVariableOne: number, initialGuessVaribaleTwo: number, maxDifference?: number, maxTries?: number): [number, number] | null {
+export function perform2dGradientDescent(functionToMinimize: ((variableOne: number, variableTwo: number) => number | null), initialGuessVariableOne: number, initialGuessVaribaleTwo: number, maxDifference?: number, maxTries?: number, variableOneBounds?: [number, number], variableTwoBounds?: [number, number]): [number, number] | null {
 	let variableOne = initialGuessVariableOne;
 	let variableTwo = initialGuessVaribaleTwo;
 	let bestValue = functionToMinimize(variableOne, variableTwo);
@@ -649,6 +644,18 @@ export function perform2dGradientDescent(functionToMinimize: ((variableOne: numb
 			break;
 		}
 
+		if (variableOneBounds) {
+			if (variableOneUpdate < variableOneBounds[0] || variableOneUpdate > variableOneBounds[1]) {
+				break;
+			}
+		}
+
+		if (variableTwoBounds) {
+			if (variableTwoUpdate < variableTwoBounds[0] || variableTwoUpdate > variableTwoBounds[1]) {
+				break;
+			}
+		}
+
 		variableOne = variableOneUpdate;
 		variableTwo = variableTwoUpdate;
 		bestValue = bestValueImprovement;
@@ -992,17 +999,26 @@ export function findCheapestIntercept(startingSituation: OrbitalCoordinates, tar
 			
 		}
 
-		let extraStartOrbits = Math.floor(startOrbitMeanAnomaly / (2 * Math.PI));
+		let startTrueAnomaly: number = 0;
-		let extraEndOrbits = Math.floor(endOrbitMeanAnomaly / (2 * Math.PI));
+		let endTrueAnomaly: number = 0;
 
-		startOrbitMeanAnomaly = startOrbitMeanAnomaly % (2 * Math.PI);
+		["start", "end"].forEach(orbit => {
-		endOrbitMeanAnomaly = endOrbitMeanAnomaly % (2 * Math.PI);
+			let meanAnomaly = orbit == "start" ? startOrbitMeanAnomaly : endOrbitMeanAnomaly;
+			let eccentricity = orbit == "start" ? startingSituation.orbit.eccentricity : targetSituation.orbit.eccentricity;
+			let extraOrbits = 0;
+			if (eccentricity < 1) {
+				extraOrbits = Math.floor(meanAnomaly / (2 * Math.PI));
+				meanAnomaly = meanAnomaly % (2 * Math.PI);
+			}
+			let [_, trueAnomaly] = getEccentricAndTrueAnomalyFromMeanAnomaly(meanAnomaly, eccentricity);
+			trueAnomaly += extraOrbits * 2 * Math.PI;
 
-		let [_, startTrueAnomaly] = getEccentricAndTrueAnomalyFromMeanAnomaly(startOrbitMeanAnomaly, startingSituation.orbit.eccentricity);
+			if (orbit == "start") {
-		let [__, endTrueAnomaly] = getEccentricAndTrueAnomalyFromMeanAnomaly(endOrbitMeanAnomaly, targetSituation.orbit.eccentricity);
+				startTrueAnomaly = trueAnomaly;
-
+			} else if (orbit == "end") {
-		startTrueAnomaly += extraStartOrbits * 2 * Math.PI;
+				endTrueAnomaly = trueAnomaly;
-		endTrueAnomaly += extraEndOrbits * 2 * Math.PI;
+			}
+		});
 
 		let targetTransferTime = endTime - startTime;
 		let lambertSolutions = new LambertSolutions(startingSituation.orbit, startTrueAnomaly, targetSituation.orbit, endTrueAnomaly + extraTrueAnomaly, body, backwards);
@@ -1024,7 +1040,12 @@ export function findCheapestIntercept(startingSituation: OrbitalCoordinates, tar
 				return null;
 			}
 
-			if (transfer.farthestPointDistance > body.sphereOfInfluence || transfer.closestPointDistance < body.closestSafeDistance) {
+			if (transfer.farthestPointDistance >= body.sphereOfInfluence || transfer.closestPointDistance <= body.closestSafeDistance) {
+				return null;
+			}
+
+			// Sometimes, we get eccentric orbits that are impossible to do
+			if (transfer.transferOrbitTrueAnomalyAtManoeuvreOne < Math.PI && transfer.transferOrbitTrueanomalyAtManoeuvreTwo > Math.PI && transfer.transferOrbit.eccentricity - 1 >= 0) {
 				return null;
 			}
 
@@ -1066,7 +1087,7 @@ export function findCheapestIntercept(startingSituation: OrbitalCoordinates, tar
 				return;
 			}
 
-			let result = perform2dGradientDescent(interceptFunction, foundStartTime, foundEndTime, 0, 30);
+			let result = perform2dGradientDescent(interceptFunction, foundStartTime, foundEndTime, 0, 30, [0, 1e99]);
 			if (result) {
 				let foundTransfer = getTransfer(result[0], result[1], backwards);
 				if (foundTransfer) {
@@ -1337,7 +1358,7 @@ export function findCheapestLanding(startingCoordinates: OrbitalCoordinates, sta
 	return bestLanding;
 }
 
-export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates, interpolationParameters: InterpolationParameters, body: Body): [OrbitalCoordinates, number] | null {
+export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates, interpolationParameters: InterpolationParameters, body: Body): [OrbitalCoordinates, number][] {
 	let targetPeriapsis = interpolationParameters.targetPeriapsis;
 
 	let targetEccentricity;
@@ -1390,6 +1411,14 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 	const targetFirstTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
 	const targetSecondTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
 
+	// We can use the true anomalies and such to calculate the distance between the two target points
+	let angleDifference = targetSecondTrueAnomaly - targetFirstTrueAnomaly;
+	let expectedDistance = Math.sqrt(
+		interpolationParameters.firstTargetAltitude**2 
+		+ interpolationParameters.secondTargetAltitude**2 
+		- 2 * interpolationParameters.firstTargetAltitude * interpolationParameters.secondTargetAltitude * Math.cos(angleDifference)
+	);
+
 	const getVectors = (trueAnomaly: number, orbit: Orbit) => {
 		const radius = orbit.semiLatusRectum / (1 + orbit.eccentricity * Math.cos(trueAnomaly));
 		const localX = radius * Math.cos(trueAnomaly);
@@ -1425,285 +1454,194 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 	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);
+	let results: [OrbitalCoordinates, number][] = [];
-
-	// 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)
-			)
-		);
-
-		let p2 = addVector(
-			c2,
-			addVector(
-				multiplyMatrixWithScalar(Math.cos(angleTwo), n21),
-				multiplyMatrixWithScalar(Math.sin(angleTwo), n22)
-			)
-		);
-
-		return getVectorMagnitude(addVector(p2, multiplyMatrixWithScalar(-1, p1))) - expectedDistance;
-	}
-
-	const distancePartialDerivativeAngleOne = (angleOne: number, angleTwo: number) => {
-		let distance = distanceFunction(angleOne, angleTwo) + expectedDistance;
-		return (
-			Math.sin(angleOne) * vectorDotProduct(n11, c2)
-			+ Math.sin(angleOne) * Math.cos(angleTwo) * vectorDotProduct(n21, n11)
-			- Math.sin(angleOne) * Math.cos(angleOne) * vectorDotProduct(n11, n11)
-			- Math.cos(angleOne) * vectorDotProduct(n12, c2)
-			- Math.cos(angleOne) * Math.sin(angleTwo) * vectorDotProduct(n12, n22)
-			+ Math.sin(angleOne) * Math.sin(angleOne) * vectorDotProduct(n12, n12)
-		) / distance;
-	}
-
-	const distancePartialDerivativeAngleTwo = (angleOne: number, angleTwo: number) => {
-		let distance = distanceFunction(angleOne, angleTwo) + expectedDistance;
-		return (
-			Math.cos(angleTwo) * Math.sin(angleTwo) * vectorDotProduct(n22, n22)
-			- Math.cos(angleTwo) * vectorDotProduct(n22, c1)
-			- 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 horizontalTwo = addVector(
-			c2,
-			multiplyMatrixWithScalar(Math.cos(angleTwo), n21)
-		);
-
-		return vectorDotProduct(horizontalOne, horizontalTwo) / (getVectorMagnitude(horizontalOne) * getVectorMagnitude(horizontalTwo)) - Math.cos(phaseAngleDifference);
-	}
-
-	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;
-
-			while (Math.abs(trialFunctionVector[0][0]) >= Math.abs(functionVector[0][0])) {
-				counter += 1;
-				if (counter >= 16) {
-					deadEnd = true;
-					break;
-				}
-
-				update = multiplyMatrixWithScalar(0.5, update);
-				trialFunctionVector = getFunctionVector(addVector(anglesVector, multiplyMatrixWithScalar(-1, update)));
-			}
-
-			if (deadEnd) {
-				break;
-			}*/
-			anglesVector = addVector(anglesVector, multiplyMatrixWithScalar(-1, update));
-		};
-		return anglesVector;
-	}
-
-	const getAnglesFromPhaseAngles = (firstPhaseAngle: number, secondPhaseAngle: number) => {
-		let angleOne;
-		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);
-		}
-
-		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);
-		}
-
-		return [angleOne, angleTwo];
-	}
-
-	let bestAngles = [[0], [0]];
-	let bestDistance: number | null = null;
-
-	for (var i = -100; i < 101; i++) {
-		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(
-			multiplyMatrixWithScalar(Math.cos(bestAngles[0][0]) * firstDistancePerpendicularToDirection, firstVectors[1]),
-			multiplyMatrixWithScalar(Math.sin(bestAngles[0][0]) * firstDistancePerpendicularToDirection, firstVectors[2])
-		)
-	);
-
-	let targetPositionTwo = addVector(multiplyMatrixWithScalar(secondDistanceAlongDirection, secondVectors[0]),
-		addVector(
-			multiplyMatrixWithScalar(Math.cos(bestAngles[1][0]) * secondDistancePerpendicularToDirection, secondVectors[1]),
-			multiplyMatrixWithScalar(Math.sin(bestAngles[1][0]) * secondDistancePerpendicularToDirection, secondVectors[2])
-		)
-	);
-
-	let normalVector = normalizeVector(vectorCrossProduct(targetPositionTwo, targetPositionOne));
-	
-	// Rotate the position vector about this normal vector to get the direction of the periapsis
-	let ux = normalVector[0][0];
-	let uy = normalVector[1][0];
-	let uz = normalVector[2][0];
-
-	let rotationMatrix = [
-		[
-			ux*ux*(1 - Math.cos(-targetFirstTrueAnomaly)) + Math.cos(-targetFirstTrueAnomaly),
-			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 periapsisVector = normalizeVector(matrixMultiply(rotationMatrix, targetPositionOne));
-	let targetOrbit: Orbit = {
-		semiLatusRectum: targetSemiLatusRectum,
-		eccentricity: targetEccentricity,
-		coordinateAxes: [
-			periapsisVector,
-			normalizeVector(vectorCrossProduct(normalVector, periapsisVector)),
-			normalVector
-		]
-	};
-
-	let targetCoordinates: OrbitalCoordinates = getOrbitalCoordinatesFromAltitude(
-		interpolationParameters.secondTargetAltitude,
-		targetHeadedInwards,
-		targetOrbit,
-		body
-	);
-
 	let epochTrueAnomaly = ownCoordinates.trueAnomaly;
 	while (epochTrueAnomaly + 2 * Math.PI < ownFirstTrueAnomaly) {
 		epochTrueAnomaly += 2 * Math.PI;
 	}
-
 	const timeElapsed = getTimeBetweenTrueAnomalies(epochTrueAnomaly, ownSecondTrueAnomaly, ownCoordinates.orbit, body);
-	return [
+
-		targetCoordinates,
+	// Try all four possible arrangements of angles
-		timeElapsed
+	[-1, 1].forEach(angleOneMultiplier => {
-	]
+		[-1, 1].forEach(angleTwoMultiplier => {
+
+			const phaseAnglesToDistanceFunctionWithDerivatives = (phaseAngleOne: number, phaseAngleTwo: number): [number, number, number] => {
+				let angleOne = angleOneMultiplier * Math.acos(firstDistanceAlongDirection*Math.tan(phaseAngleOne) / firstDistancePerpendicularToDirection);
+				let angleTwo = angleTwoMultiplier * Math.acos(secondDistanceAlongDirection*Math.tan(phaseAngleTwo) / secondDistancePerpendicularToDirection);
+
+				let positionOne = addVector(multiplyMatrixWithScalar(firstDistanceAlongDirection, firstVectors[0]),
+					addVector(
+						multiplyMatrixWithScalar(firstDistancePerpendicularToDirection * Math.cos(angleOne), firstVectors[1]),
+						multiplyMatrixWithScalar(firstDistancePerpendicularToDirection * Math.sin(angleOne), firstVectors[2])
+					)
+				);
+
+				let positionTwo = addVector(multiplyMatrixWithScalar(secondDistanceAlongDirection, secondVectors[0]),
+					addVector(
+						multiplyMatrixWithScalar(secondDistancePerpendicularToDirection * Math.cos(angleTwo), secondVectors[1]),
+						multiplyMatrixWithScalar(secondDistancePerpendicularToDirection * Math.sin(angleTwo), secondVectors[2])						
+					)
+				);
+
+				let difference = addVector(positionTwo, multiplyMatrixWithScalar(-1, positionOne));
+				let distance = getVectorMagnitude(difference);
+
+				// Find derivative with respect to phase angles
+				let dAngleOneByPhaseAngleOne = -firstDistanceAlongDirection / (firstDistancePerpendicularToDirection * Math.cos(phaseAngleOne)**2 * Math.sin(angleOne));
+				let dAngleTwoByPhaseAngleTwo = -secondDistanceAlongDirection / (secondDistancePerpendicularToDirection * Math.cos(phaseAngleTwo)**2 * Math.sin(angleTwo));
+
+				let dPositionOneByPhaseAngleOne = addVector(
+					multiplyMatrixWithScalar(-Math.sin(angleOne)*dAngleOneByPhaseAngleOne*firstDistancePerpendicularToDirection, firstVectors[1]),
+					multiplyMatrixWithScalar(Math.cos(angleOne)*dAngleOneByPhaseAngleOne*firstDistancePerpendicularToDirection, firstVectors[2])
+				);
+
+				let dPositionTwoByPhaseAngleTwo = addVector(
+					multiplyMatrixWithScalar(-Math.sin(angleTwo)*dAngleTwoByPhaseAngleTwo*secondDistancePerpendicularToDirection, secondVectors[1]),
+					multiplyMatrixWithScalar(Math.cos(angleTwo)*dAngleTwoByPhaseAngleTwo*secondDistancePerpendicularToDirection, secondVectors[2])
+				);
+
+				let dDistanceByPhaseAngleOne = -vectorDotProduct(difference, dPositionOneByPhaseAngleOne) / distance;
+				let dDistanceByPhaseAngleTwo = vectorDotProduct(difference, dPositionTwoByPhaseAngleTwo) / distance;
+
+				return [distance - expectedDistance, dDistanceByPhaseAngleOne, dDistanceByPhaseAngleTwo];
+			}
+
+			// Try to find the phase angles that minimise the function above
+			let phaseAngleOne = interpolationParameters.firstPhaseAngle;
+			let phaseAngleTwo = interpolationParameters.secondPhaseAngle;
+
+			for (var i = 0; i < 1000; i++) {
+				// Minimise phase angles one by one
+				let [value, dfd1, dfd2] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne, phaseAngleTwo);
+				
+				let updateOne = -value / dfd1;
+				let learningRate = 1;
+
+				while(learningRate*updateOne > Math.PI / 3600 && Math.abs(phaseAngleOne + learningRate * updateOne - interpolationParameters.firstPhaseAngle) > Math.PI / 3600) {
+					learningRate *= 0.5;
+				}
+
+				// We'll also check if the update skips over a minimum of the function
+				let deadEndAngleOne = false;
+				let [testValue, _, __] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne + learningRate * updateOne, phaseAngleTwo);
+				let counter = 0;
+				while (Math.abs(testValue) > Math.abs(value)) {
+					counter++;
+					if (counter >= 32) {
+						deadEndAngleOne = true;
+						break;
+					}
+
+					learningRate *= 0.5;
+					[testValue, _, __] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne + learningRate * updateOne, phaseAngleTwo);
+				}
+
+				if (!deadEndAngleOne) {
+					phaseAngleOne += learningRate * updateOne;
+				}
+
+				// Now, do the same for angle two
+				[value, dfd1, dfd2] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne, phaseAngleTwo);
+				let updateTwo = -value / dfd2;
+				learningRate = 1;
+
+				while (learningRate * updateTwo > Math.PI / 36000 && Math.abs(phaseAngleTwo + learningRate * updateTwo - interpolationParameters.secondPhaseAngle) > Math.PI / 3600) {
+					learningRate *= 0.5;
+				}
+
+				let deadEndAngleTwo = false;
+				[testValue, _, __] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne, phaseAngleTwo + learningRate * updateTwo);
+				counter = 0;
+				while (Math.abs(testValue) > Math.abs(value)) {
+					counter++;
+					if (counter >= 32) {
+						deadEndAngleTwo = true;
+						break;
+					}
+
+					learningRate *= 0.5;
+					[testValue, _, __] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne, phaseAngleTwo + learningRate * updateTwo);
+				}
+
+				if (!deadEndAngleTwo) {
+					phaseAngleTwo += learningRate * updateTwo;
+				}
+
+				if (deadEndAngleOne && deadEndAngleTwo) {
+					break;
+				}
+			};
+
+			// If we're still far away from a good solution, stop here
+			let [value, _, __] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne, phaseAngleTwo);
+			if (Math.abs(value) > 1000 || isNaN(value)) {
+				return;
+			}
+
+			let angleOne = angleOneMultiplier * Math.acos(firstDistanceAlongDirection * Math.tan(phaseAngleOne) / firstDistancePerpendicularToDirection)
+			let angleTwo = angleTwoMultiplier * Math.acos(secondDistanceAlongDirection * Math.tan(phaseAngleTwo) / secondDistancePerpendicularToDirection);
+
+
+			let targetPositionOne = addVector(multiplyMatrixWithScalar(firstDistanceAlongDirection, firstVectors[0]),
+				addVector(
+					multiplyMatrixWithScalar(Math.cos(angleOne) * firstDistancePerpendicularToDirection, firstVectors[1]),
+					multiplyMatrixWithScalar(Math.sin(angleOne) * firstDistancePerpendicularToDirection, firstVectors[2])
+				)
+			);
+
+			let targetPositionTwo = addVector(multiplyMatrixWithScalar(secondDistanceAlongDirection, secondVectors[0]),
+				addVector(
+					multiplyMatrixWithScalar(Math.cos(angleTwo) * secondDistancePerpendicularToDirection, secondVectors[1]),
+					multiplyMatrixWithScalar(Math.sin(angleTwo) * secondDistancePerpendicularToDirection, secondVectors[2])
+				)
+			);
+
+			let normalVector = normalizeVector(vectorCrossProduct(targetPositionOne, targetPositionTwo));
+			
+			// Rotate the position vector about this normal vector to get the direction of the periapsis
+			let ux = normalVector[0][0];
+			let uy = normalVector[1][0];
+			let uz = normalVector[2][0];
+
+			let rotationMatrix = [
+				[
+					ux*ux*(1 - Math.cos(-targetFirstTrueAnomaly)) + Math.cos(-targetFirstTrueAnomaly),
+					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 periapsisVector = normalizeVector(matrixMultiply(rotationMatrix, targetPositionOne));
+			let targetOrbit: Orbit = {
+				semiLatusRectum: targetSemiLatusRectum,
+				eccentricity: targetEccentricity,
+				coordinateAxes: [
+					periapsisVector,
+					normalizeVector(vectorCrossProduct(normalVector, periapsisVector)),
+					normalVector
+				]
+			};
+
+			let targetCoordinates: OrbitalCoordinates = getOrbitalCoordinatesFromAltitude(
+				interpolationParameters.secondTargetAltitude,
+				targetHeadedInwards,
+				targetOrbit,
+				body
+			);
+
+			results.push([targetCoordinates, timeElapsed]);
+		});
+	});
+
+	return results;
 }

src/gui/intercept.ts

@@ -2,9 +2,9 @@ import { createLabel, createNumberInput, createRadioButton, getCoordinatesFromPa
 import { OrbitalParametersGui } from "./orbit";
 import { type Body } from "../calculations/constants";
 import type { FindBestInterceptMessage, ProgressMessage } from "./worker";
-import type { ChangingStorageValue } from "../storage";
+import { ChangingStorageValue } from "../storage";
 import { ManoeuvresGui } from "./manoeuvres";
-import { extrapolateTrajectory, findOrbitThroughInterpolation, type OrbitalCoordinates } from "../calculations/orbit-calculations";
+import { extrapolateTrajectory, type OrbitalCoordinates } from "../calculations/orbit-calculations";
 import { InterpolateOrbitGui, type InterpolationParameters } from "./interpolate";
 
 export type TargetOrbitChoice = "KnownOrbit" | "InterpolateOrbit";
@@ -24,6 +24,8 @@ export class InterceptTargetGui {
 	additionalTrueAnomaly: ChangingStorageValue<number>;
 	targetOrbitChoice: ChangingStorageValue<TargetOrbitChoice>;
 	interpolationParameters: ChangingStorageValue<InterpolationParameters>;
+	interpolatedOrbitCoordinates: ChangingStorageValue<OrbitalCoordinates | null>;
+	interpolatedOrbitTime: ChangingStorageValue<number | null>;
 
 	parentDiv: HTMLDivElement;
 	orbitGui: OrbitalParametersGui;
@@ -42,9 +44,11 @@ export class InterceptTargetGui {
 		this.additionalTrueAnomaly = additionalTrueAnomaly;
 		this.targetOrbitChoice = targetOrbitChoice;
 		this.interpolationParameters = interpolationParameters;
+		this.interpolatedOrbitCoordinates = new ChangingStorageValue<OrbitalCoordinates | null>(null);
+		this.interpolatedOrbitTime = new ChangingStorageValue<number | null>(null);
 
 		this.orbitGui = new OrbitalParametersGui(targetOrbitalParameters, "orbitWithPosition");
-		this.interpolateGui = new InterpolateOrbitGui(this.interpolationParameters);
+		this.interpolateGui = new InterpolateOrbitGui(this.interpolationParameters, this.startingOrbitalParameters, this.interpolatedOrbitCoordinates, this.interpolatedOrbitTime, this.body);
 		this.manoeuvresGui = new ManoeuvresGui(true);
 
 		this.worker = null;
@@ -114,24 +118,14 @@ export class InterceptTargetGui {
 				this.worker = null;
 				searchButton.innerHTML = "Search for cheapest intercept";
 			} else {
-				searchButton.innerHTML = "Cancel search";
-
 				let currentTime = this.currentTime.getCurrentValue();
 				let body = this.body.getCurrentValue();
 				let startingOrbitalParameters = this.startingOrbitalParameters.getCurrentValue();
 				let targetOrbitalParameters = this.targetOrbitalParameters.getCurrentValue();
 				let additionalTrueAnomaly = this.additionalTrueAnomaly.getCurrentValue();
 				let targetOrbitChoice = this.targetOrbitChoice.getCurrentValue();
-				let interpolationParameters = this.interpolationParameters.getCurrentValue();
+				let interpolatedOrbitCoordinates = this.interpolatedOrbitCoordinates.getCurrentValue();
-
+				let interpolatedOrbitTime = this.interpolatedOrbitTime.getCurrentValue();
-				interpolationParameters = structuredClone(interpolationParameters);
-				interpolationParameters.firstOwnAltitude += body.radius;
-				interpolationParameters.firstTargetAltitude += body.radius;
-				interpolationParameters.secondOwnAltitude += body.radius;
-				interpolationParameters.secondTargetAltitude += body.radius;
-				interpolationParameters.targetAltitude += body.radius;
-				interpolationParameters.targetPeriapsis += body.radius;
-				interpolationParameters.targetApoapsis += body.radius;
 
 
 				let startingCoordinates: OrbitalCoordinates | null = getCoordinatesFromParameters(startingOrbitalParameters, body);
@@ -141,10 +135,9 @@ export class InterceptTargetGui {
 
 				let targetCoordinates: OrbitalCoordinates | null = null;
 				if (targetOrbitChoice == "InterpolateOrbit") {
-					let result = findOrbitThroughInterpolation(startingCoordinates, interpolationParameters, body);
+					if (interpolatedOrbitCoordinates && interpolatedOrbitTime) {
-					if (result) {
+						targetCoordinates = interpolatedOrbitCoordinates;
-						targetCoordinates = result[0];
+						targetStartTime = interpolatedOrbitTime;
-						targetStartTime = startingOrbitalParameters.currentTimeAtReading + result[1];
 					}
 				} else if (targetOrbitChoice = "KnownOrbit") {
 					targetCoordinates = getCoordinatesFromParameters(targetOrbitalParameters, body);
@@ -167,6 +160,7 @@ export class InterceptTargetGui {
 
 				if (startingCoordinates && targetCoordinates && startTime) {
 					this.worker = new Worker(new URL('./worker.ts', import.meta.url), {type: 'module'});
+					searchButton.innerHTML = "Cancel search";
 					this.worker.addEventListener("message", event => {
 						let transferResponse = event.data as ProgressMessage;
 						if (transferResponse) {

src/gui/interpolate.ts

@@ -1,5 +1,8 @@
+import type { Body } from "../calculations/constants";
+import { extrapolateTrajectory, findOrbitThroughInterpolation, type OrbitalCoordinates } from "../calculations/orbit-calculations";
 import type { ChangingStorageValue } from "../storage";
-import { createLabel, createNumberInput, createRadioButton } from "./common";
+import { createLabel, createNumberInput, createRadioButton, getCoordinatesFromParameters, type OrbitalParameters } from "./common";
+import { Renderer } from "./renderer";
 
 export interface InterpolationParameters {
 	targetPeriapsis: number,
@@ -7,7 +10,6 @@ export interface InterpolationParameters {
 	targetApoapsis: number,
 	targetSpeed: number,
 	targetAltitude: number,
-	targetAbovePlane: boolean,
 	firstOwnAltitude: number,
 	firstTargetAltitude: number,
 	firstDistance: number,
@@ -24,7 +26,6 @@ export const DefaultInterpolationParameters: InterpolationParameters = {
 	targetApoapsis: 200000,
 	targetSpeed: 0,
 	targetAltitude: 0,
-	targetAbovePlane: true,
 	firstOwnAltitude: 0,
 	firstTargetAltitude: 0,
 	firstDistance: 0,
@@ -37,12 +38,22 @@ export const DefaultInterpolationParameters: InterpolationParameters = {
 
 export class InterpolateOrbitGui {
 	parentDiv: HTMLDivElement;
-	interpolationParameters: ChangingStorageValue<InterpolationParameters>
+	startingOrbitalParameters: ChangingStorageValue<OrbitalParameters>;
+	body: ChangingStorageValue<Body>;
+	interpolationParameters: ChangingStorageValue<InterpolationParameters>;
+	interpolatedOrbitCoordinates: ChangingStorageValue<OrbitalCoordinates | null>;
+	interpolatedOrbitTime: ChangingStorageValue<number | null>;
 	sourceId: string;
+	renderDiv: HTMLDivElement;
 
-	constructor(interpolationParameters: ChangingStorageValue<InterpolationParameters>) {
+	constructor(interpolationParameters: ChangingStorageValue<InterpolationParameters>, startingOrbitParameters: ChangingStorageValue<OrbitalParameters>, interpolatedOrbitCoordinates: ChangingStorageValue<OrbitalCoordinates | null>, interpolatedOrbitTime: ChangingStorageValue<number | null>, body: ChangingStorageValue<Body>) {
 		this.parentDiv = document.createElement("div");
+		this.startingOrbitalParameters = startingOrbitParameters;
 		this.interpolationParameters = interpolationParameters;
+		this.interpolatedOrbitCoordinates = interpolatedOrbitCoordinates;
+		this.interpolatedOrbitTime = interpolatedOrbitTime;
+		this.renderDiv = document.createElement("div");
+		this.body = body;
 		
 		let targetOrbitHeader = document.createElement("h4");
 		targetOrbitHeader.appendChild(document.createTextNode("Describe target orbit:"));
@@ -80,13 +91,6 @@ export class InterpolateOrbitGui {
 		let speedInput = createInputWithLabel("Target speed:", 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");
 		instantOneHeader.appendChild(document.createTextNode("Measurement one:"));
 		this.parentDiv.appendChild(instantOneHeader);
@@ -105,13 +109,17 @@ export class InterpolateOrbitGui {
 		let secondDistanceInput = createInputWithLabel("Distance to target:", 0);
 		let secondPhaseAngleInput = createInputWithLabel("Phase angle:", -180, 180);
 
+		let interpolateOrbitButton = document.createElement("button");
+		interpolateOrbitButton.appendChild(document.createTextNode("Find interpolated orbit"));
+		this.parentDiv.appendChild(interpolateOrbitButton);
+		this.parentDiv.appendChild(this.renderDiv);
+
 		this.sourceId = crypto.randomUUID();
 		const onChange = () => {
 			let periapsis = parseFloat(periapsisInput.value);
 			let apoapsis = parseFloat(apoapsisInput.value);
 			let speed = parseFloat(speedInput.value);
 			let altitude = parseFloat(altitudeInput.value);
-			let abovePlane = abovePlaneButton.checked;
 			let firstTargetAltitude = parseFloat(firstTargetAltitudeInput.value);
 			let firstOwnAltitude = parseFloat(firstOwnAltitudeInput.value);
 			let firstDistance = parseFloat(firstDistanceInput.value);
@@ -132,7 +140,6 @@ export class InterpolateOrbitGui {
 				targetApoapsis: apoapsis,
 				targetSpeed: speed,
 				targetAltitude: altitude,
-				targetAbovePlane: abovePlane,
 				firstTargetAltitude: firstTargetAltitude,
 				firstOwnAltitude: firstOwnAltitude,
 				firstDistance: firstDistance,
@@ -151,7 +158,6 @@ export class InterpolateOrbitGui {
 			apoapsisInput,
 			speedInput,
 			altitudeInput,
-			abovePlaneButton,
 			firstOwnAltitudeInput,
 			firstTargetAltitudeInput,
 			firstDistanceInput,
@@ -183,11 +189,6 @@ export class InterpolateOrbitGui {
 			apoapsisInput.value = value.targetApoapsis.toString();
 			speedInput.value = value.targetSpeed.toString();
 			altitudeInput.value = value.targetAltitude.toString();
-			if (value.targetAbovePlane) {
-				abovePlaneButton.setAttribute("checked", "");
-			} else {
-				abovePlaneButton.removeAttribute("checked");
-			}
 			firstOwnAltitudeInput.value = value.firstOwnAltitude.toString();
 			firstTargetAltitudeInput.value = value.firstTargetAltitude.toString();
 			firstDistanceInput.value = value.firstDistance.toString();
@@ -198,5 +199,61 @@ export class InterpolateOrbitGui {
 			secondDistanceInput.value = value.secondDistance.toString();
 			secondPhaseAngleInput.value = (value.secondPhaseAngle * 180 / Math.PI).toString();
 		});
+
+
+		// Interpolation action
+		interpolateOrbitButton.addEventListener("click", _ => {
+			this.renderDiv.innerHTML = "";
+
+			let startingOrbitalParameters = this.startingOrbitalParameters.getCurrentValue();
+			let interpolationParameters = structuredClone(this.interpolationParameters.getCurrentValue());
+			let body = this.body.getCurrentValue();
+
+			// All interpolation altitude parameters need to get the planet radius added to them
+			interpolationParameters.firstTargetAltitude += body.radius;
+			interpolationParameters.firstOwnAltitude += body.radius;
+			interpolationParameters.secondTargetAltitude += body.radius;
+			interpolationParameters.secondOwnAltitude += body.radius;
+			interpolationParameters.targetAltitude += body.radius;
+			interpolationParameters.targetApoapsis += body.radius;
+			interpolationParameters.targetPeriapsis += body.radius;
+			
+			let ownCoordinates = getCoordinatesFromParameters(startingOrbitalParameters, body);
+			let results = findOrbitThroughInterpolation(ownCoordinates, interpolationParameters, body);
+
+			if (results.length == 1) {
+				this.interpolatedOrbitCoordinates.set(results[0][0], "null");
+				this.interpolatedOrbitTime.set(results[0][1], "null");
+			} else {
+				let chooseHeader = document.createElement("h3");
+				chooseHeader.appendChild(document.createTextNode("Choose interpolated orbit:"));
+				this.renderDiv.appendChild(chooseHeader);
+
+				let interpolationChoiceId = crypto.randomUUID();
+				results.forEach(([coordinates, time]) => {
+					let buttonId = crypto.randomUUID();
+					let button = createRadioButton(interpolationChoiceId, buttonId);
+					let label = createLabel(buttonId, "The orbit below fits best");
+					let renderer = new Renderer(body);
+					let extrapolatedOwnCoordinates = extrapolateTrajectory(time, ownCoordinates, body);
+					let usedCoordinates = ownCoordinates;
+					if (extrapolatedOwnCoordinates) {
+						usedCoordinates = extrapolatedOwnCoordinates;
+					}
+					renderer.addCoordinates(usedCoordinates, 0x00ffff);
+					renderer.addCoordinates(coordinates, 0xff00ff);
+
+					this.renderDiv.appendChild(button);
+					this.renderDiv.appendChild(label);
+					this.renderDiv.appendChild(renderer.parentDiv);
+					this.renderDiv.appendChild(document.createElement("br"));
+
+					button.addEventListener("change", _ => {
+						this.interpolatedOrbitCoordinates.set(coordinates, "null");
+						this.interpolatedOrbitTime.set(startingOrbitalParameters.currentTimeAtReading + time, "null");
+					});
+				});
+			}
+		});
 	}
 }

src/gui/landing.ts

@@ -1,7 +1,7 @@
 import type { Body } from "../calculations/constants";
-import { extrapolateTrajectory, findCheapestLanding, getOrbitalCoordinates, type OrbitalCoordinates, type ShipParameters } from "../calculations/orbit-calculations";
+import { extrapolateTrajectory, type OrbitalCoordinates, type ShipParameters } from "../calculations/orbit-calculations";
 import type { ChangingStorageValue } from "../storage";
-import { createDisabledInput, createLabel, getCoordinatesFromParameters, getOrbitFromParameters, type OrbitalParameters } from "./common";
+import { createDisabledInput, createLabel, getCoordinatesFromParameters, type OrbitalParameters } from "./common";
 import { type LandingParameters } from "../calculations/orbit-calculations";
 import { LandingParametersGui } from "./landingparameters";
 import { ShipGui } from "./ship";

src/gui/renderer.ts

@@ -0,0 +1,186 @@
+import * as THREE from 'three';
+import { OrbitControls } from 'three/examples/jsm/Addons.js';
+import type { Body } from '../calculations/constants';
+import type { Orbit, OrbitalCoordinates } from '../calculations/orbit-calculations';
+import { addVector, multiplyMatrixWithScalar, normalizeVector } from '../calculations/mathematics';
+
+function getPosition(trueAnomaly: number, orbit: Orbit): number[][] {
+	var radius = orbit.semiLatusRectum / (1 + orbit.eccentricity * Math.cos(trueAnomaly));
+	var localX = radius * Math.cos(trueAnomaly);
+	var localY = radius * Math.sin(trueAnomaly);
+
+	let position = addVector(
+		multiplyMatrixWithScalar(localX, orbit.coordinateAxes[0]),
+		multiplyMatrixWithScalar(localY, orbit.coordinateAxes[1])
+	);
+
+	return position;
+}
+
+function getScalingFunction(mesh: THREE.Mesh) {
+	return function (_: number, cameraPosition: THREE.Vector3, cameraDirection: THREE.Vector3) {
+		let displacement = new THREE.Vector3();
+		displacement.subVectors(cameraPosition, mesh.position);
+		let distanceFromCamera = displacement.dot(cameraDirection);
+		mesh.scale.setScalar(distanceFromCamera / 5000000);
+	}
+}
+
+export class Renderer {
+	parentDiv: HTMLDivElement;
+
+	scene: THREE.Scene
+	camera: THREE.PerspectiveCamera;
+	renderer: THREE.WebGLRenderer;
+	controls: OrbitControls;
+
+	body: Body;
+
+	animationFunctions: ((time: number, cameraPosition: THREE.Vector3, cameraDirection: THREE.Vector3) => void)[];
+
+	constructor(body: Body) {
+		this.parentDiv = document.createElement("div");
+		this.body = body;
+		this.animationFunctions = [];
+
+		// Rendering
+		this.scene = new THREE.Scene();
+		this.scene.background = new THREE.Color(0x888888);
+		this.camera = new THREE.PerspectiveCamera(75, 400 / 400, 0.1, 1e99);
+		this.camera.position.set(0, 0, 5000000);
+		this.renderer = new THREE.WebGLRenderer();
+
+		this.renderer.setSize(400, 400);
+		this.parentDiv.appendChild(this.renderer.domElement);
+
+		this.controls = new OrbitControls(this.camera, this.renderer.domElement);	
+		this.controls.enableDamping = true;
+		this.controls.dampingFactor = 0.05;
+		this.controls.screenSpacePanning = false;
+		this.controls.minDistance = this.body.radius * 2;
+		this.controls.maxDistance = 1e99;
+		this.controls.cursorStyle = 'grab';
+		this.controls.maxPolarAngle = Math.PI;
+
+		this.renderer.setAnimationLoop(this.animate.bind(this));
+
+		let bodySphere = new THREE.SphereGeometry(this.body.radius);
+		let material = new THREE.MeshBasicMaterial({color: 0x000000});
+		let mesh = new THREE.Mesh(bodySphere, material);
+		this.scene.add(mesh);
+	}
+
+	animate ( time: number) {
+		this.controls.update();
+
+		let cameraPosition = new THREE.Vector3();
+		this.camera.getWorldPosition(cameraPosition);
+		let cameraDirection = new THREE.Vector3();
+		this.camera.getWorldDirection(cameraDirection);
+		this.animationFunctions.forEach(f => f(time, cameraPosition, cameraDirection));
+		this.renderer.render(this.scene, this.camera);
+	}
+
+	addOrbit(orbit: Orbit) {
+		let stableOrbit = false;
+		if (orbit.eccentricity < 1) {
+			let maxDistance = orbit.semiLatusRectum / (1 - orbit.eccentricity);
+			if (maxDistance < this.body.sphereOfInfluence) {
+				stableOrbit = true;
+			}
+		}
+
+		let minimumTrueAnomaly = -Math.PI;
+		let maximumTrueAnomaly = Math.PI
+
+		if (!stableOrbit) {
+			maximumTrueAnomaly = Math.acos((orbit.semiLatusRectum - this.body.sphereOfInfluence) / (this.body.sphereOfInfluence * orbit.eccentricity));
+			minimumTrueAnomaly = -maximumTrueAnomaly;
+		}
+
+		let angles = [];
+		for (var i = 0; i <= 360; i++) {
+			let angle = minimumTrueAnomaly + (maximumTrueAnomaly - minimumTrueAnomaly) * i / 360;
+			angles.push(angle);
+		}
+
+		let points: THREE.Vector3[] = [];
+		angles.forEach(angle => {
+			let position = getPosition(angle, orbit);
+
+			points.push(new THREE.Vector3(
+				position[0][0],
+				position[2][0],
+				-position[1][0]
+			));
+		});
+
+		const geometry = new THREE.BufferGeometry().setFromPoints(points);
+		const material = new THREE.LineBasicMaterial({color: 0xffffff});
+
+		const line = new THREE.Line(geometry, material);
+
+		this.scene.add(line);
+
+		if (orbit.eccentricity > 0.001) {
+			// Add the periapsis
+			let periapsisPoint = new THREE.SphereGeometry(this.body.radius / 20);
+			let periapsisMaterial = new THREE.MeshBasicMaterial({color: 0x00ff00});
+			let periapsisMesh = new THREE.Mesh(periapsisPoint, periapsisMaterial);
+
+			let periapsisPosition = getPosition(0, orbit);
+			periapsisMesh.position.x = periapsisPosition[0][0];
+			periapsisMesh.position.y = periapsisPosition[2][0];
+			periapsisMesh.position.z = -periapsisPosition[1][0];
+
+			this.scene.add(periapsisMesh);
+			this.animationFunctions.push(getScalingFunction(periapsisMesh));
+
+			if (stableOrbit) {
+				let apoapsisPoint = new THREE.SphereGeometry(this.body.radius / 20);
+				let apoapsisMaterial = new THREE.MeshBasicMaterial({color: 0xff0000});
+				let apoapsisMesh = new THREE.Mesh(apoapsisPoint, apoapsisMaterial);
+
+				let apoapsisPosition = getPosition(Math.PI, orbit);
+				apoapsisMesh.position.x = apoapsisPosition[0][0];
+				apoapsisMesh.position.y = apoapsisPosition[2][0];
+				apoapsisMesh.position.z = -apoapsisPosition[1][0];
+
+				this.scene.add(apoapsisMesh);
+				this.animationFunctions.push(getScalingFunction(apoapsisMesh));
+			}
+		}
+	}
+
+	addCoordinates(coordinates: OrbitalCoordinates, color: THREE.ColorRepresentation) {
+		this.addOrbit(coordinates.orbit);
+
+		let pointSphere = new THREE.SphereGeometry(this.body.radius / 5);
+		let material = new THREE.MeshBasicMaterial({color: color});
+		let mesh = new THREE.Mesh(pointSphere, material);
+
+		let position = getPosition(coordinates.trueAnomaly, coordinates.orbit);
+
+		mesh.position.x = position[0][0];
+		mesh.position.y = position[2][0];
+		mesh.position.z = -position[1][0];
+		
+		this.scene.add(mesh);
+
+		// Create an arrow that points in the correct direction
+		let nextPosition = getPosition(coordinates.trueAnomaly + 0.0001, coordinates.orbit);
+		let direction = normalizeVector(addVector(nextPosition, multiplyMatrixWithScalar(-1, position)));
+		let renderDirection = new THREE.Vector3(direction[0][0], direction[2][0], -direction[1][0]);
+
+		let arrowHelper = new THREE.ArrowHelper(renderDirection, mesh.position, this.body.radius / 1.666, color, this.body.radius / 5, this.body.radius/5);
+		this.scene.add(arrowHelper);
+
+		this.animationFunctions.push((_, cameraPosition) => {
+			let displacement = new THREE.Vector3();
+			displacement.subVectors(cameraPosition, mesh.position);
+			let distance = displacement.length();
+			mesh.scale.setScalar(distance / 5000000);
+			arrowHelper.scale.setScalar(distance / 5000000);
+		});
+	}
+}

src/gui/worker.ts

@@ -1,6 +1,6 @@
 import type { Body } from "../calculations/constants";
 import { findCheapestIntercept, findCheapestLanding, findCheapestTransfer, type Orbit, type OrbitalCoordinates, type ShipParameters, type Transfer } from "../calculations/orbit-calculations";
-import type { Landing, LandingParameters, LandingProgressCallbackFunction, Manoeuvre } from "../calculations/orbit-calculations";
+import type { Landing, LandingParameters, LandingProgressCallbackFunction } from "../calculations/orbit-calculations";
 
 const ctx: Worker = self as any;