Intermediate stage

src/calculations/orbit-calculations.ts

@@ -1361,36 +1361,6 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 	let ownShipHeadedInwards = interpolationParameters.secondOwnAltitude < interpolationParameters.firstOwnAltitude;
 	let targetHeadedInwards = interpolationParameters.secondTargetAltitude < interpolationParameters.firstTargetAltitude;
 
-	// Define some helper functions for later
-	const getRandomPerpendicularVector = (vector: number[][]) => {
-		let perpendicularX = Math.random();
-		let perpendicularY = Math.random();
-		let perpendicularZ = Math.random();
-
-		let normalVector: number[][];
-		if (Math.abs(vector[0][0]) > 0.0001) {
-			normalVector = [
-				[(-vector[1][0] * perpendicularY - vector[2][0]*perpendicularZ) / vector[0][0]],
-				[perpendicularY],
-				[perpendicularZ]
-			];
-		} else if (Math.abs(vector[1][0]) > 0.0001) {
-			normalVector = [
-				[perpendicularX],
-				[(-vector[0][0] / perpendicularX - vector[2][0] / perpendicularZ) / vector[1][0]],
-				[perpendicularZ]
-			];
-		} else {
-			normalVector = [
-				[perpendicularX],
-				[perpendicularY],
-				[(-vector[0][0] * perpendicularX - vector[1][0] * perpendicularY) / vector[2][0]]
-			]
-		};
-
-		return normalizeVector(normalVector);
-	};
-
 	const getTrueAnomalyFromAltitude = (altitude: number, semiLatusRectum: number, eccentricity: number, headingInwards: boolean) => {
 		let trueAnomaly = Math.acos((semiLatusRectum - altitude) / (altitude * eccentricity));
 		if (headingInwards) {
@@ -1400,7 +1370,7 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 	}
 
 
-	const getPossibleTargetVectors = (ownOrbit: Orbit, ownAltitude: number, targetAltitude: number, distanceToTarget: number): [number[][], number[][], number[][]] => {
+	const getPossibleTargetVectors = (ownOrbit: Orbit, ownAltitude: number, targetAltitude: number, distanceToTarget: number): [number[][], number[][], number[][], number, number] => {
 		let ownTrueAnomaly = getTrueAnomalyFromAltitude(ownAltitude, ownOrbit.semiLatusRectum, ownOrbit.eccentricity, ownShipHeadedInwards);
 
 		let ownLocalX = ownAltitude * Math.cos(ownTrueAnomaly);
@@ -1413,7 +1383,7 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 			)
 		);
 
-		let perpendicularInPlane = normalizeVector(vectorCrossProduct(ownDirection, ownOrbit.coordinateAxes[2]));
+		let perpendicularInPlane = normalizeVector(vectorCrossProduct(ownOrbit.coordinateAxes[2], ownDirection));
 		let perpendicularOutOfPlane = normalizeVector(vectorCrossProduct(ownDirection, perpendicularInPlane));
 
 		// Find phase angle to target
@@ -1426,114 +1396,134 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 		let vectorOne = multiplyMatrixWithScalar(distanceAlongNormal, perpendicularInPlane);
 		let vectorTwo = multiplyMatrixWithScalar(distanceAlongNormal, perpendicularOutOfPlane);
 
-		return [root, vectorOne, vectorTwo];
+		return [root, vectorOne, vectorTwo, distanceAlongDirection, distanceAlongNormal];
-	}
-
-	const getFunctionAndDerivatives = (possiblePositionOne: [number[][], number[][], number[][]], possiblePositionTwo: [number[][], number[][], number[][]], altitudeOne: number, altitudeTwo: number, anomalyDifference: number): [(angleOne: number, angleTwo: number) => number, (angleOne: number, angleTwo: number) => number, (angleOne: number, angleTwo: number) => number] => {
-		let c1 = possiblePositionOne[0];
-		let p11 = possiblePositionOne[1];
-		let p12 = possiblePositionOne[2];
-
-		let c2 = possiblePositionTwo[0];
-		let p21 = possiblePositionTwo[1];
-		let p22 = possiblePositionTwo[2];
-
-		// Scale everything to avoid numerical explosions
-		let scaleFactor = 1 / getVectorMagnitude(c2);
-
-		c1 = multiplyMatrixWithScalar(scaleFactor, c1);
-		p11 = multiplyMatrixWithScalar(scaleFactor, p11);
-		p12 = multiplyMatrixWithScalar(scaleFactor, p12);
-		c2 = multiplyMatrixWithScalar(scaleFactor, c2);
-		p21 = multiplyMatrixWithScalar(scaleFactor, p21);
-		p22 = multiplyMatrixWithScalar(scaleFactor, p22);
-
-		let scaledAltitudeOne = altitudeOne * scaleFactor;
-		let scaledAltitudeTwo = altitudeTwo * scaleFactor;
-
-		const constantPart = scaledAltitudeOne * scaledAltitudeTwo * Math.cos(anomalyDifference);
-
-		const functionToMinimize = (angleOne: number, angleTwo: number) => {
-			let d1 = addVector(c1, addVector(
-				multiplyMatrixWithScalar(Math.cos(angleOne), p11),
-				multiplyMatrixWithScalar(Math.sin(angleOne), p12)
-			));
-
-			let d2 = addVector(c2, addVector(
-				multiplyMatrixWithScalar(Math.cos(angleTwo), p21),
-				multiplyMatrixWithScalar(Math.sin(angleTwo), p22)
-			));
-
-			return vectorDotProduct(d1, d2) - constantPart;
-		}
-
-		const partialDerivativeAngleOne = (angleOne: number, angleTwo: number) => {
-			let d2 = addVector(c2, addVector(
-				multiplyMatrixWithScalar(Math.cos(angleTwo), p21),
-				multiplyMatrixWithScalar(Math.sin(angleTwo), p22)
-			));
-			
-			let d1Dx = addVector(
-				multiplyMatrixWithScalar(-Math.sin(angleOne), p11),
-				multiplyMatrixWithScalar(Math.cos(angleOne), p12)
-			);
-
-			return vectorDotProduct(d1Dx, d2);
-		}
-
-		const partialDerivativeAngleTwo = (angleOne: number, angleTwo: number) => {
-			let d1 = addVector(c1, addVector(
-				multiplyMatrixWithScalar(Math.cos(angleOne), p11),
-				multiplyMatrixWithScalar(Math.sin(angleOne), p12)
-			));
-
-			let d2Dy = addVector(
-				multiplyMatrixWithScalar(-Math.sin(angleTwo), p21),
-				multiplyMatrixWithScalar(Math.cos(angleTwo), p22)
-			);
-
-			return vectorDotProduct(d1, d2Dy);
-		}
-
-		return [functionToMinimize, partialDerivativeAngleOne, partialDerivativeAngleTwo];
 	}
 
 	let firstPositionPossibleVectors = getPossibleTargetVectors(ownCoordinates.orbit, interpolationParameters.firstOwnAltitude, interpolationParameters.firstTargetAltitude, interpolationParameters.firstDistance);
 	let secondPositionPossibleVectors = getPossibleTargetVectors(ownCoordinates.orbit, interpolationParameters.secondOwnAltitude, interpolationParameters.secondTargetAltitude, interpolationParameters.secondDistance);
 
-	// Since we know a lot about the target orbit, we can figure out how far it is supposed to be between the two points
+	let targetFirstTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
-	let firstTargetTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
+	let targetSecondTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
-	let secondTargetTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
+
+	let c1 = firstPositionPossibleVectors[0];
+	let p11 = firstPositionPossibleVectors[1];
+	let p12 = firstPositionPossibleVectors[2];
+
+	let c2 = secondPositionPossibleVectors[0];
+	let p21 = secondPositionPossibleVectors[1];
+	let p22 = secondPositionPossibleVectors[2];
+
+	// Scale everything to avoid numerical explosions
+	let scaleFactor = 1 / getVectorMagnitude(c2);
+
+	c1 = multiplyMatrixWithScalar(scaleFactor, c1);
+	p11 = multiplyMatrixWithScalar(scaleFactor, p11);
+	p12 = multiplyMatrixWithScalar(scaleFactor, p12);
+	c2 = multiplyMatrixWithScalar(scaleFactor, c2);
+	p21 = multiplyMatrixWithScalar(scaleFactor, p21);
+	p22 = multiplyMatrixWithScalar(scaleFactor, p22);
+
+	let scaledAltitudeOne = interpolationParameters.firstTargetAltitude * scaleFactor;
+	let scaledAltitudeTwo = interpolationParameters.secondTargetAltitude * scaleFactor;
+
+	const constantPart = scaledAltitudeOne * scaledAltitudeTwo * Math.cos(targetSecondTrueAnomaly - targetFirstTrueAnomaly);
+
+	const functionToMinimize = (angleOne: number, angleTwo: number) => {
+		let d1 = addVector(c1, addVector(
+			multiplyMatrixWithScalar(Math.cos(angleOne), p11),
+			multiplyMatrixWithScalar(Math.sin(angleOne), p12)
+		));
+
+		let d2 = addVector(c2, addVector(
+			multiplyMatrixWithScalar(Math.cos(angleTwo), p21),
+			multiplyMatrixWithScalar(Math.sin(angleTwo), p22)
+		));
+
+		return vectorDotProduct(d1, d2) - constantPart;
+	}
+
+	const partialDerivativeAngleOne = (angleOne: number, angleTwo: number) => {
+		let d2 = addVector(c2, addVector(
+			multiplyMatrixWithScalar(Math.cos(angleTwo), p21),
+			multiplyMatrixWithScalar(Math.sin(angleTwo), p22)
+		));
+		
+		let d1Dx = addVector(
+			multiplyMatrixWithScalar(-Math.sin(angleOne), p11),
+			multiplyMatrixWithScalar(Math.cos(angleOne), p12)
+		);
+
+		return vectorDotProduct(d1Dx, d2);
+	}
+
+	const partialDerivativeAngleTwo = (angleOne: number, angleTwo: number) => {
+		let d1 = addVector(c1, addVector(
+			multiplyMatrixWithScalar(Math.cos(angleOne), p11),
+			multiplyMatrixWithScalar(Math.sin(angleOne), p12)
+		));
+
+		let d2Dy = addVector(
+			multiplyMatrixWithScalar(-Math.sin(angleTwo), p21),
+			multiplyMatrixWithScalar(Math.cos(angleTwo), p22)
+		);
+
+		return vectorDotProduct(d1, d2Dy);
+	}
+
+	let firstOwnTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards);
+	let secondOwnTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards);
+
+	let horizontalAngleDifference = (secondOwnTrueAnomaly - interpolationParameters.secondPhaseAngle) - (firstOwnTrueAnomaly - interpolationParameters.firstPhaseAngle);
+
+	const phaseAngleFunction = (angleOne: number, angleTwo: number): number => {
+		let d1 = addVector(c1, multiplyMatrixWithScalar(Math.cos(angleOne), p11));
+		let d2 = addVector(c2, multiplyMatrixWithScalar(Math.cos(angleTwo), p21));
+		return vectorDotProduct(d1, d2) - Math.cos(horizontalAngleDifference) * getVectorMagnitude(d1) * getVectorMagnitude(d2);
+	}
+
+	const phaseAnglePartialDerivativeAngleOne = (angleOne: number, angleTwo: number) => {
+		let d1 = addVector(c1, multiplyMatrixWithScalar(Math.cos(angleOne), p11));
+		let d1Dx = multiplyMatrixWithScalar(-Math.sin(angleOne), p11);
+		let d2 = addVector(c2, multiplyMatrixWithScalar(Math.cos(angleTwo), p21));
+
+		return vectorDotProduct(d1Dx, d2) - Math.cos(horizontalAngleDifference) * getVectorMagnitude(d2) * Math.sin(angleOne) * Math.cos(angleOne) * firstPositionPossibleVectors[4]**2 / getVectorMagnitude(d1);
+	}
 	
-	let [minimizeFunction, partialDerivativeAngleOne, partialDerivativeAngleTwo] = getFunctionAndDerivatives(firstPositionPossibleVectors, secondPositionPossibleVectors, interpolationParameters.firstTargetAltitude, interpolationParameters.secondTargetAltitude, secondTargetTrueAnomaly - firstTargetTrueAnomaly);
 
 	const gradientDescent = (startingGuessAngleOne: number, startingGuessAngleTwo: number): [number, number, number] => {
 		let angleOne = startingGuessAngleOne;
 		let angleTwo = startingGuessAngleTwo;
 
+		let distanceAlongDirection = getVectorMagnitude(firstPositionPossibleVectors[0]);
+		let distanceAlongPerpendicular = getVectorMagnitude(firstPositionPossibleVectors[1]);
+
 		let value = minimizeFunction(angleOne, angleTwo);
 		// Do up to 100 minimizing steps
 		for (var i = 0; i < 100; i++) {
 			let dfDx = partialDerivativeAngleOne(angleOne, angleTwo);
 			let dfDy = partialDerivativeAngleTwo(angleOne, angleTwo);
 
-			let scaling = Math.sqrt(dfDx**2 + dfDy**2);
+			let phaseAngleValue = distanceAlongPerpendicular * Math.cos(angleOne) / distanceAlongDirection - Math.tan(interpolationParameters.firstPhaseAngle);
-			let directionX = -dfDy / scaling;
+			let phaseAngleDerivative = -distanceAlongPerpendicular * Math.sin(angleOne) / distanceAlongDirection;
-			let directionY = -dfDx / scaling;
 
-			let changeAlongDirection = directionX * dfDx + directionY * dfDy;
+			let changeAngleOne = dfDx * phaseAngleValue / (dfDy * phaseAngleDerivative);
+			let changeAngleTwo = (dfDx * phaseAngleValue - phaseAngleDerivative * value) / (dfDy * phaseAngleDerivative);
 
-			let deadEnd = false;
-			let trialValue = value;
 			let trialAngleOne = angleOne;
 			let trialAngleTwo = angleTwo;
+			let trialValue = value;
 			let divisor = 1;
 
+			// We don't want things to blow up
+			while (Math.abs(changeAngleOne) / 2**divisor > Math.PI / 1000 || Math.abs(changeAngleTwo) / 2**divisor > Math.PI / 1000) {
+				divisor +=1;
+			}
+
+			let deadEnd = false;
+
 			while (Math.abs(trialValue) >= Math.abs(value)) {
-				let stepAngleOne = - value * directionX / (changeAlongDirection * 2 ** divisor);
+				trialAngleOne = angleOne - changeAngleOne / (2 ** divisor);
-				let stepAngleTwo = - value * directionY / (changeAlongDirection * 2 ** divisor);
+				trialAngleTwo = angleTwo - changeAngleTwo / (2 ** divisor);
-				trialAngleOne = angleOne + stepAngleOne;
-				trialAngleTwo = angleTwo + stepAngleTwo;
 				trialValue = minimizeFunction(trialAngleOne, trialAngleTwo);
 				divisor += 1;
 
@@ -1555,24 +1545,29 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 		return [angleOne, angleTwo, value];
 	}
 
+	let firstCenterLength = getVectorMagnitude(firstPositionPossibleVectors[0]);
+	let firstPerpendicularLength = getVectorMagnitude(firstPositionPossibleVectors[1]);
+
+	let secondCenterLength = getVectorMagnitude(secondPositionPossibleVectors[0]);
+	let secondPerpendicularLength = getVectorMagnitude(secondPositionPossibleVectors[1]);
+
+	let angleOneGuess = Math.acos(firstCenterLength * Math.tan(interpolationParameters.firstPhaseAngle) / firstPerpendicularLength);
+	let angleTwoGuess = Math.acos(secondCenterLength * Math.tan(interpolationParameters.secondPhaseAngle) / secondPerpendicularLength);
+
 	let bestAngleOne = 0;
 	let bestAngleTwo = 0;
-	let bestValue = null;
+	let bestValue: number | null = null;
 
-	// Try to start at a bunch of different positions and do gradient descent from there
+	[angleOneGuess, -angleOneGuess].forEach(angleOne => {
-	for (var i = 0; i < 500; i++) {
+		[angleTwoGuess, -angleTwoGuess].forEach(angleTwo => {
-		for (var j = 0; j < 500; j++) {
+			let value = minimizeFunction(angleOne, angleTwo);
-			let startingAngleOne = i * 2 * Math.PI / 500.0;
+			if (bestValue == null || value < bestValue) {
-			let startingAngleTwo = j * 2 * Math.PI / 500.0
+				bestValue = value;
-			let startingValue = minimizeFunction(startingAngleOne, startingAngleTwo);
+				bestAngleOne = angleOne;
-
+				bestAngleTwo = angleTwo;
-			if (bestValue == null || Math.abs(startingValue) < Math.abs(bestValue)) {
-				bestAngleOne = startingAngleOne;
-				bestAngleTwo = startingAngleTwo;
-				bestValue = startingValue;
 			}
-		}
+		});
-	}
+	});
 
 	[bestAngleOne, bestAngleOne, bestValue] = gradientDescent(bestAngleOne, bestAngleTwo);
 

src/gui/interpolate.ts

@@ -10,9 +10,11 @@ export interface InterpolationParameters {
 	firstOwnAltitude: number,
 	firstTargetAltitude: number,
 	firstDistance: number,
+	firstPhaseAngle: number,
 	secondOwnAltitude: number,
 	secondTargetAltitude: number,
 	secondDistance: number,
+	secondPhaseAngle: number,
 }
 
 export const DefaultInterpolationParameters: InterpolationParameters = {
@@ -24,9 +26,11 @@ export const DefaultInterpolationParameters: InterpolationParameters = {
 	firstOwnAltitude: 0,
 	firstTargetAltitude: 0,
 	firstDistance: 0,
+	firstPhaseAngle: 0,
 	secondOwnAltitude: 0,
 	secondTargetAltitude: 0,
-	secondDistance: 0
+	secondDistance: 0,
+	secondPhaseAngle: 0
 }
 
 export class InterpolateOrbitGui {
@@ -81,6 +85,7 @@ export class InterpolateOrbitGui {
 		let firstOwnAltitudeInput = createInputWithLabel("Own altitude:", 0);
 		let firstTargetAltitudeInput = createInputWithLabel("Target altitude:", 0);
 		let firstDistanceInput = createInputWithLabel("Distance to target:", 0);
+		let firstPhaseAngleInput = createInputWithLabel("Phase angle:", -180, 180);
 
 		let instantTwoHeader = document.createElement("h4");
 		instantTwoHeader.appendChild(document.createTextNode("Measurement two:"));
@@ -89,6 +94,7 @@ export class InterpolateOrbitGui {
 		let secondOwnAltitudeInput = createInputWithLabel("Own altitude:", 0);
 		let secondTargetAltitudeInput = createInputWithLabel("Target altitude:", 0);
 		let secondDistanceInput = createInputWithLabel("Distance to target:", 0);
+		let secondPhaseAngleInput = createInputWithLabel("Phase angle:", -180, 180);
 
 		this.sourceId = crypto.randomUUID();
 		const onChange = () => {
@@ -99,9 +105,11 @@ export class InterpolateOrbitGui {
 			let firstTargetAltitude = parseFloat(firstTargetAltitudeInput.value);
 			let firstOwnAltitude = parseFloat(firstOwnAltitudeInput.value);
 			let firstDistance = parseFloat(firstDistanceInput.value);
+			let firstPhaseAngle = parseFloat(firstPhaseAngleInput.value) * Math.PI / 180.0;
 			let secondTargetAltitude = parseFloat(secondTargetAltitudeInput.value);
 			let secondOwnAltitude = parseFloat(secondOwnAltitudeInput.value);
 			let secondDistance = parseFloat(secondDistanceInput.value);
+			let secondPhaseAngle = parseFloat(secondPhaseAngleInput.value) * Math.PI / 180.0;
 
 			let orbitChoice: "apoapsis" | "speedAndAltitude" = "apoapsis";
 			if (speedAndAltitudeButton.checked) {
@@ -117,9 +125,11 @@ export class InterpolateOrbitGui {
 				firstTargetAltitude: firstTargetAltitude,
 				firstOwnAltitude: firstOwnAltitude,
 				firstDistance: firstDistance,
+				firstPhaseAngle: firstPhaseAngle,
 				secondTargetAltitude: secondTargetAltitude,
 				secondOwnAltitude: secondOwnAltitude,
 				secondDistance: secondDistance,
+				secondPhaseAngle: secondPhaseAngle
 			}, this.sourceId)
 		};
 
@@ -133,9 +143,11 @@ export class InterpolateOrbitGui {
 			firstOwnAltitudeInput,
 			firstTargetAltitudeInput,
 			firstDistanceInput,
+			firstPhaseAngleInput,
 			secondOwnAltitudeInput,
 			secondTargetAltitudeInput,
-			secondDistanceInput
+			secondDistanceInput,
+			secondPhaseAngleInput
 		].forEach(input => input.addEventListener("change", onChange));
 
 		interpolationParameters.listenToValue((source, value) => {
@@ -162,10 +174,12 @@ export class InterpolateOrbitGui {
 			firstOwnAltitudeInput.value = value.firstOwnAltitude.toString();
 			firstTargetAltitudeInput.value = value.firstTargetAltitude.toString();
 			firstDistanceInput.value = value.firstDistance.toString();
+			firstPhaseAngleInput.value = (value.firstPhaseAngle * 180 / Math.PI).toString();
 
 			secondOwnAltitudeInput.value = value.secondOwnAltitude.toString();
 			secondTargetAltitudeInput.value = value.secondTargetAltitude.toString();
 			secondDistanceInput.value = value.secondDistance.toString();
+			secondPhaseAngleInput.value = (value.secondPhaseAngle * 180 / Math.PI).toString();
 		});
 	}
 }