Optimize over phase angles

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 { addVector, getVectorMagnitude, matrixMultiply, multiplyMatrixWithScalar, normalizeVector, subtractVector, vectorCrossProduct, vectorDotProduct } from "./mathematics";
+import { addVector, getVectorMagnitude, invertTwoByTwoMatrix, matrixMultiply, multiplyMatrixWithScalar, normalizeVector, subtractVector, vectorCrossProduct, vectorDotProduct } from "./mathematics";
 
 export interface Orbit {
 	semiLatusRectum: number,
@@ -1358,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;
@@ -1411,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);
@@ -1446,105 +1454,145 @@ 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);
 
-	// Now, try some Newton's method to estimate the two position's the target has gone through
+	let results: [OrbitalCoordinates, number][] = [];
+	let epochTrueAnomaly = ownCoordinates.trueAnomaly;
+	while (epochTrueAnomaly + 2 * Math.PI < ownFirstTrueAnomaly) {
+		epochTrueAnomaly += 2 * Math.PI;
+	}
+	const timeElapsed = getTimeBetweenTrueAnomalies(epochTrueAnomaly, ownSecondTrueAnomaly, ownCoordinates.orbit, body);
 
-	// To avoid the maths exploding, we'll scale everything down a bit
+	// Try all four possible arrangements of angles
-	const scaleFactor = 1;
+	[-1, 1].forEach(angleOneMultiplier => {
+		[-1, 1].forEach(angleTwoMultiplier => {
 
-	const c1 = multiplyMatrixWithScalar(firstDistanceAlongDirection * scaleFactor, firstVectors[0]);
+			const phaseAnglesToDistanceFunctionWithDerivatives = (phaseAngleOne: number, phaseAngleTwo: number): [number, number, number] => {
-	const n11 = multiplyMatrixWithScalar(firstDistancePerpendicularToDirection * scaleFactor, firstVectors[1]);
+				let angleOne = angleOneMultiplier * Math.acos(firstDistanceAlongDirection*Math.tan(phaseAngleOne) / firstDistancePerpendicularToDirection);
-	const n12 = multiplyMatrixWithScalar(firstDistancePerpendicularToDirection * scaleFactor, firstVectors[2]);
+				let angleTwo = angleTwoMultiplier * Math.acos(secondDistanceAlongDirection*Math.tan(phaseAngleTwo) / secondDistancePerpendicularToDirection);
 
-	const c2 = multiplyMatrixWithScalar(secondDistanceAlongDirection * scaleFactor, secondVectors[0]);
+				let positionOne = addVector(multiplyMatrixWithScalar(firstDistanceAlongDirection, firstVectors[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(firstDistancePerpendicularToDirection * Math.cos(angleOne), firstVectors[1]),
-				multiplyMatrixWithScalar(Math.sin(angleOne), n12)
+						multiplyMatrixWithScalar(firstDistancePerpendicularToDirection * Math.sin(angleOne), firstVectors[2])
 					)
 				);
 
-		let p2 = addVector(
+				let positionTwo = addVector(multiplyMatrixWithScalar(secondDistanceAlongDirection, secondVectors[0]),
-			c2,
 					addVector(
-				multiplyMatrixWithScalar(Math.cos(angleTwo), n21),
+						multiplyMatrixWithScalar(secondDistancePerpendicularToDirection * Math.cos(angleTwo), secondVectors[1]),
-				multiplyMatrixWithScalar(Math.sin(angleTwo), n22)
+						multiplyMatrixWithScalar(secondDistancePerpendicularToDirection * Math.sin(angleTwo), secondVectors[2])						
 					)
 				);
 
-		return getVectorMagnitude(addVector(p2, multiplyMatrixWithScalar(-1, p1))) - expectedDistance;
+				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];
 			}
 
-	const getAnglesFromPhaseAngles = (firstPhaseAngle: number, secondPhaseAngle: number) => {
+			// Try to find the phase angles that minimise the function above
-		let cosAngleOne;
+			let phaseAngleOne = interpolationParameters.firstPhaseAngle;
-		let cosAngleTwo;
+			let phaseAngleTwo = interpolationParameters.secondPhaseAngle;
 
-		if (Math.abs(interpolationParameters.firstPhaseAngle) < Math.PI / 2) {
+			for (var i = 0; i < 1000; i++) {
-			cosAngleOne = Math.tan(firstPhaseAngle) * firstDistanceAlongDirection / firstDistancePerpendicularToDirection;
+				// Minimise phase angles one by one
-		} else {
+				let [value, dfd1, dfd2] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne, phaseAngleTwo);
-			cosAngleOne = Math.tan(Math.PI - firstPhaseAngle) * -firstDistanceAlongDirection / firstDistancePerpendicularToDirection;
+				
+				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;
 				}
 
-		if (Math.abs(interpolationParameters.secondPhaseAngle) < Math.PI / 2) {
+				// We'll also check if the update skips over a minimum of the function
-			cosAngleTwo = Math.tan(secondPhaseAngle) * secondDistanceAlongDirection / secondDistancePerpendicularToDirection;
+				let deadEndAngleOne = false;
-		} else {
+				let [testValue, _, __] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne + learningRate * updateOne, phaseAngleTwo);
-			cosAngleTwo = Math.tan(Math.PI - secondPhaseAngle) * -secondDistanceAlongDirection / secondDistancePerpendicularToDirection;
+				let counter = 0;
+				while (Math.abs(testValue) > Math.abs(value)) {
+					counter++;
+					if (counter >= 32) {
+						deadEndAngleOne = true;
+						break;
 					}
 
-		if (Math.abs(cosAngleOne) > 1) {
+					learningRate *= 0.5;
-			cosAngleOne = 1 * Math.sign(cosAngleOne);
+					[testValue, _, __] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne + learningRate * updateOne, phaseAngleTwo);
 				}
 
-		if (Math.abs(cosAngleTwo) > 1) {
+				if (!deadEndAngleOne) {
-			cosAngleTwo = 1 * Math.sign(cosAngleTwo);
+					phaseAngleOne += learningRate * updateOne;
 				}
 
-		return [Math.acos(cosAngleOne), Math.acos(cosAngleTwo)];
+				// 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 bestAngles = [[0], [0]];
+				let deadEndAngleTwo = false;
-	let bestDistance: number | null = null;
+				[testValue, _, __] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne, phaseAngleTwo + learningRate * updateTwo);
-
+				counter = 0;
-	for (var i = -100; i < 101; i++) {
+				while (Math.abs(testValue) > Math.abs(value)) {
-		for (var j = -100; j < 101; j++) {
+					counter++;
-			let firstPhaseAngle = interpolationParameters.firstPhaseAngle + i * Math.PI / 18000;
+					if (counter >= 32) {
-			let secondPhaseAngle = interpolationParameters.secondPhaseAngle + j * Math.PI / 18000;
+						deadEndAngleTwo = true;
-
+						break;
-			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) {
+					learningRate *= 0.5;
-		bestAngles = [[-bestAngles[0][0]], [-bestAngles[1][0]]];
+					[testValue, _, __] = phaseAnglesToDistanceFunctionWithDerivatives(phaseAngleOne, phaseAngleTwo + learningRate * updateTwo);
 				}
 
-	// Do some additional Newtoning on the best angles we've found so far
+				if (!deadEndAngleTwo) {
-	//bestAngles = performNewtonMethod(bestAngles, 10000);
+					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(bestAngles[0][0]) * firstDistancePerpendicularToDirection, firstVectors[1]),
+					multiplyMatrixWithScalar(Math.cos(angleOne) * firstDistancePerpendicularToDirection, firstVectors[1]),
-			multiplyMatrixWithScalar(Math.sin(bestAngles[0][0]) * firstDistancePerpendicularToDirection, firstVectors[2])
+					multiplyMatrixWithScalar(Math.sin(angleOne) * firstDistancePerpendicularToDirection, firstVectors[2])
 				)
 			);
 
 			let targetPositionTwo = addVector(multiplyMatrixWithScalar(secondDistanceAlongDirection, secondVectors[0]),
 				addVector(
-			multiplyMatrixWithScalar(Math.cos(bestAngles[1][0]) * secondDistancePerpendicularToDirection, secondVectors[1]),
+					multiplyMatrixWithScalar(Math.cos(angleTwo) * secondDistancePerpendicularToDirection, secondVectors[1]),
-			multiplyMatrixWithScalar(Math.sin(bestAngles[1][0]) * secondDistancePerpendicularToDirection, secondVectors[2])
+					multiplyMatrixWithScalar(Math.sin(angleTwo) * secondDistancePerpendicularToDirection, secondVectors[2])
 				)
 			);
 
@@ -1591,14 +1639,9 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 				body
 			);
 
-	let epochTrueAnomaly = ownCoordinates.trueAnomaly;
+			results.push([targetCoordinates, timeElapsed]);
-	while (epochTrueAnomaly + 2 * Math.PI < ownFirstTrueAnomaly) {
+		});
-		epochTrueAnomaly += 2 * Math.PI;
+	});
-	}
 
-	const timeElapsed = getTimeBetweenTrueAnomalies(epochTrueAnomaly, ownSecondTrueAnomaly, ownCoordinates.orbit, body);
+	return results;
-	return [
-		targetCoordinates,
-		timeElapsed
-	]
 }

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/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);
+		});
+	}
+}