Much guessing

src/calculations/orbit-calculations.ts

@@ -555,6 +555,97 @@ export function getSpeed(trueAnomaly: number, orbit: Orbit, planet: Body): numbe
 	return Math.sqrt(planet.gravitationalParameter * (1 + 2 * orbit.eccentricity * Math.cos(trueAnomaly) + orbit.eccentricity**2) / orbit.semiLatusRectum);
 }
 
+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 {
+	let variableOne = initialGuessVariableOne;
+	let variableTwo = initialGuessVaribaleTwo;
+	let bestValue = functionToMinimize(variableOne, variableTwo);
+
+	if (bestValue == null) {
+		return null;
+	}
+
+	maxDifference = maxDifference ?? 0.5;
+	maxTries = maxTries ?? 1000;
+
+	// Do some gradient descent on the best values found so far
+	let tries = 0;
+	while (bestValue != null && bestValue > maxDifference && tries < maxTries) {
+		tries++;
+
+		let dfx = functionToMinimize(variableOne+0.0001, variableTwo);
+		let dfy = functionToMinimize(variableOne, variableTwo+0.0001);
+
+		let estimateDfDx = null;
+		let estimateDfDy = null;
+		
+		if (dfx) {
+			estimateDfDx = (dfx - bestValue) / 0.0001;
+		}
+
+		if (dfy) {
+			estimateDfDy = (functionToMinimize(variableOne, variableTwo+0.0001) ?? bestValue - bestValue) / 0.0001;
+		}
+
+		let direction: number[][] | null = null;
+
+		if (estimateDfDx && estimateDfDy) {
+			let slopeNormal = vectorCrossProduct(
+				[[1], [0], [estimateDfDx]],
+				[[0], [1], [estimateDfDy]]
+			);	
+
+			direction = normalizeVector([[slopeNormal[0][0]], [slopeNormal[1][0]], [0]]);
+		} else if (estimateDfDx) {
+			direction = [[1], [0], [0]];
+		} else if (estimateDfDy) {
+			direction = [[0], [1], [0]];
+		}
+
+		if (!direction) {
+			break;
+		}
+
+		let df = functionToMinimize(variableOne + direction[0][0]*0.0001, variableTwo + direction[1][0]*0.0001);
+		if (!df) {
+			break;
+		}
+
+		let estimateDfDv = (df - bestValue) / 0.0001;
+
+		let bestValueImprovement: number = bestValue;
+		let variableOneUpdate = variableOne;
+		let variableTwoUpdate = variableTwo;
+		let divisor = 0;
+		let deadEnd = false;
+
+		while (bestValueImprovement >= bestValue) {
+			variableOneUpdate = variableOne - bestValue * direction[0][0] / (estimateDfDv * 2**divisor);
+			variableTwoUpdate = variableTwo - bestValue * direction[1][0] / (estimateDfDv * 2**divisor);
+			bestValueImprovement = functionToMinimize(variableOneUpdate, variableTwoUpdate) ?? bestValue;
+			divisor += 1;
+
+			if (divisor >= 32) {
+				deadEnd = true;
+				break;
+			}
+		}
+
+		if (deadEnd) {
+			break;
+		}
+
+		variableOne = variableOneUpdate;
+		variableTwo = variableTwoUpdate;
+		bestValue = bestValueImprovement;
+	}
+
+	return [variableOne, variableTwo];
+}
+
 export function calculateSimplePlaneChange(coordinates: OrbitalCoordinates, planet: Body, targetInclination: number, targetLongitudeOfAscendingNode: number, circularizeOrbit: boolean): SimplePlaneChange {
 	const otherPlaneNormal = [
 		[Math.sin(targetLongitudeOfAscendingNode)*Math.sin(targetInclination)],
@@ -807,8 +898,8 @@ export function findLambertSolutionsWithCorrectTime(lambertSolutions: LambertSol
 export function findCheapestIntercept(startingSituation: OrbitalCoordinates, targetSituation: OrbitalCoordinates, body: Body, extraTrueAnomaly: number, progressCallback?: ProgressCallbackFunction): Transfer | null {
 	// Find acceptable intercept times
 	// First number is intercept time, second number is true anomaly
-	let starts: [number, number][] = [];
+	let startTimes: number[] = [];
-	let intercepts: [number, number][] = [];
+	let endTimes: number[] = [];
 
 	// Check if the starting orbit is stable0
 	let startingOrbitStable = false;
@@ -828,95 +919,147 @@ export function findCheapestIntercept(startingSituation: OrbitalCoordinates, tar
 		}
 	}
 
-	let maxStartTime = 1e99;
+	let maxEndTime: number | null = null;
-	let maxStartTrueAnomaly = 1e99;
+	let maxStartTime: number | null = null;
-	let maxEndTime = 1e99;
-	let maxEndTrueAnomaly = 1e99;
 
-	if (interceptOrbitStable && startingOrbitStable) {
+	if (!startingOrbitStable) {
-		// If both orbits are stable, we'll check for three whole orbits of the largest orbit
+		let maxStartTrueAnomaly = Math.abs(Math.acos((startingSituation.orbit.semiLatusRectum - body.sphereOfInfluence) / (body.sphereOfInfluence * startingSituation.orbit.eccentricity)));
-		let startingOrbitPeriod = getTimeBetweenTrueAnomalies(0, 2*Math.PI, startingSituation.orbit, body);
+		maxStartTime = getTimeBetweenTrueAnomalies(startingSituation.trueAnomaly, maxStartTrueAnomaly, startingSituation.orbit, body);
-		let targetOrbitPeriod = getTimeBetweenTrueAnomalies(0, 2*Math.PI, targetSituation.orbit, body);
+		// If we're leaving the galaxy, set the max end time to the time it would take to orbit at the edge of the sphere of influence
-
+		maxEndTime = getOrbitalPeriod(body.sphereOfInfluence, body);
-		maxStartTime = 3 * Math.max(startingOrbitPeriod, targetOrbitPeriod);
-		maxEndTime = 4 * Math.max(startingOrbitPeriod, targetOrbitPeriod);
-	} else {
-		maxStartTime = 1e99;
-		maxStartTrueAnomaly = 1e99;
-		maxEndTime = 1e99;
-		maxEndTrueAnomaly = 1e99;
-
-		if (!startingOrbitStable) {
-			maxStartTrueAnomaly = Math.abs(Math.acos((startingSituation.orbit.semiLatusRectum - body.sphereOfInfluence) / (body.sphereOfInfluence * startingSituation.orbit.eccentricity)));
-			maxStartTime = getTimeBetweenTrueAnomalies(startingSituation.trueAnomaly, maxStartTrueAnomaly, startingSituation.orbit, body);
-		}
-
-		if (!interceptOrbitStable) {
-			maxEndTrueAnomaly = Math.abs(Math.acos((targetSituation.orbit.semiLatusRectum - body.sphereOfInfluence) / (body.sphereOfInfluence * targetSituation.orbit.eccentricity)));
-			maxEndTime = getTimeBetweenTrueAnomalies(targetSituation.trueAnomaly, maxEndTrueAnomaly, targetSituation.orbit, body);
-		} else {
-			// Set the max end time to the max start time plus two full orbits of the target
-			maxEndTime = maxStartTime + 4 * Math.PI * Math.sqrt((targetSituation.orbit.semiLatusRectum / (1 - targetSituation.orbit.eccentricity**2))**3 / body.gravitationalParameter);
-		}
 	}
 
-	let anomalyCounter = 1;
+	if (!interceptOrbitStable) {
-	while (true) {
+		let maxEndTrueAnomaly = Math.abs(Math.acos((targetSituation.orbit.semiLatusRectum - body.sphereOfInfluence) / (body.sphereOfInfluence * targetSituation.orbit.eccentricity)));
-		let possibleStart = startingSituation.trueAnomaly + anomalyCounter * 2 * Math.PI / 100;
+		maxEndTime = Math.min(maxEndTime ?? 1e99, getTimeBetweenTrueAnomalies(targetSituation.trueAnomaly, maxEndTrueAnomaly, targetSituation.orbit, body));
-		if (possibleStart > maxStartTrueAnomaly) {
-			break;
-		}
-		let possibleStartTime = getTimeBetweenTrueAnomalies(startingSituation.trueAnomaly, possibleStart, startingSituation.orbit, body);
-		if (possibleStartTime > maxStartTime) {
-			break;
-		}
-
-		starts.push([possibleStartTime, possibleStart]);
-		anomalyCounter += 1;
 	}
 
-	anomalyCounter = 0;
+	if (startingOrbitStable && interceptOrbitStable) {
-	while (true) {
+		// If both orbits are stable, try for four of the longest orbit
-		let possibleEnd = targetSituation.trueAnomaly + anomalyCounter * 2 * Math.PI / 100.0;
+		let startingSemiMajor = startingSituation.orbit.semiLatusRectum / (1 - startingSituation.orbit.eccentricity**2);
-		if (possibleEnd > maxEndTrueAnomaly) {
+		let targetSemiMajor = targetSituation.orbit.semiLatusRectum / (1 - targetSituation.orbit.eccentricity**2);
-			break;
-		}
-		let possibleEndTime = getTimeBetweenTrueAnomalies(targetSituation.trueAnomaly, possibleEnd, targetSituation.orbit, body);
-		possibleEnd += extraTrueAnomaly;
 
-		if (possibleEndTime > maxEndTime) {
+		let startingOrbitalPeriod = getOrbitalPeriod(startingSemiMajor, body);
-			break;
+		let targetOrbitalPeriod = getOrbitalPeriod(targetSemiMajor, body);
-		}
 
-		intercepts.push([possibleEndTime, possibleEnd]);
+		maxStartTime = 4 * Math.max(startingOrbitalPeriod, targetOrbitalPeriod);
-		anomalyCounter += 1;
+		maxEndTime = 5 * Math.max(startingOrbitalPeriod, targetOrbitalPeriod);
 	}
 
-	let pairs: [number, number, number, number][] = [];
+	// Max end time should have been set by now, but Typescript doesn't seem to know that
-	starts.forEach(([startingTime, startingTrueAnomaly]) => {
+	maxEndTime = maxEndTime ?? 0;
-		intercepts.forEach(([endingTime, endingTrueAnomaly]) => {
+
-			if (endingTime > startingTime) {
+	// No point in having start times after the end times
-				pairs.push([startingTime, startingTrueAnomaly, endingTime, endingTrueAnomaly]);
+	if (!maxStartTime || maxStartTime > maxEndTime) {
+		maxStartTime = maxEndTime;
+	}
+
+	// We'll try starting at 100 different start and end times
+	let startTimeStep = maxStartTime / 100;
+	let endTimeStep = maxEndTime / 100;
+
+	for (var i = 0; i < 100; i++) {
+		startTimes.push(i * startTimeStep);
+		endTimes.push(i * endTimeStep);
+	}
+
+	// Create pairs to test
+	let pairs: [number, number][] = [];
+	startTimes.forEach(startTime => {
+		endTimes.forEach(endTime => {
+			if (endTime > startTime) {
+				pairs.push([startTime, endTime]);
 			}
-		});
+		})
 	});
 
-	let bestDeltaV: number | null = null;
+	const getTransfer = (startTime: number, endTime: number, backwards: boolean): Transfer | null => {
+		let startOrbitMeanAnomaly = startingSituation.meanAnomaly + startTime * startingSituation.meanAngularMotion;
+		let endOrbitMeanAnomaly = targetSituation.meanAnomaly + endTime * targetSituation.meanAngularMotion;
+
+		let extraStartOrbits = Math.floor(startOrbitMeanAnomaly / (2 * Math.PI));
+		let extraEndOrbits = Math.floor(endOrbitMeanAnomaly / (2 * Math.PI));
+
+		startOrbitMeanAnomaly = startOrbitMeanAnomaly % (2 * Math.PI);
+		endOrbitMeanAnomaly = endOrbitMeanAnomaly % (2 * Math.PI);
+
+		let [_, startTrueAnomaly] = getEccentricAndTrueAnomalyFromMeanAnomaly(startOrbitMeanAnomaly, startingSituation.orbit.eccentricity);
+		let [__, endTrueAnomaly] = getEccentricAndTrueAnomalyFromMeanAnomaly(endOrbitMeanAnomaly, targetSituation.orbit.eccentricity);
+
+		startTrueAnomaly += extraStartOrbits * 2 * Math.PI;
+		endTrueAnomaly += extraEndOrbits * 2 * Math.PI;
+
+		let targetTransferTime = endTime - startTime;
+		let lambertSolutions = new LambertSolutions(startingSituation.orbit, startTrueAnomaly, targetSituation.orbit, endTrueAnomaly + extraTrueAnomaly, body, backwards);
+
+		let transfer = findLambertSolutionsWithCorrectTime(lambertSolutions, targetTransferTime);
+		if (transfer) {
+			transfer.firstManoeuvre.time += startTime;
+			transfer.secondManoeuvre.time += startTime;
+		}
+
+		return transfer;
+	}
+
+	// Create the function we want to gradient descend
+	const getInterceptFunction = (backwards: boolean) => {
+		return function(startTime: number, endTime: number): number | null {
+			let transfer = getTransfer(startTime, endTime, backwards);
+			if (!transfer) {
+				return null;
+			}
+
+			if (transfer.farthestPointDistance > body.sphereOfInfluence || transfer.closestPointDistance < body.closestSafeDistance) {
+				return null;
+			}
+
+			return transfer.firstManoeuvre.totalDeltaV + transfer.secondManoeuvre.totalDeltaV;
+		}
+	};
+
+	// Start at each pair, and do a little gradient descending
+	let bestDeltaV: number | null = null;;
 	let bestTransfer: Transfer | null = null;
 
-	pairs.forEach(([startingTime, startingTrueAnomaly, endingTime, endingTrueAnomaly], index) => {
+	pairs.forEach(([startTime, endTime], index) => {
-		let transferTime = endingTime - startingTime;
+		// Try both forwards and backwards
-		[true, false].forEach(goBackwards => {
+		[true, false].forEach(backwards => {
-			let lambertSolutions = new LambertSolutions(startingSituation.orbit, startingTrueAnomaly, targetSituation.orbit, endingTrueAnomaly, body, goBackwards);
+			let interceptFunction = getInterceptFunction(backwards);
-			let transfer = findLambertSolutionsWithCorrectTime(lambertSolutions, transferTime);
-			if (transfer !== null && transfer.closestPointDistance > body.closestSafeDistance && transfer.farthestPointDistance < body.sphereOfInfluence) {
-				let totalDeltaV = transfer.firstManoeuvre.totalDeltaV + transfer.secondManoeuvre.totalDeltaV;
-				if (bestDeltaV == null || totalDeltaV < bestDeltaV) {
-					bestDeltaV = totalDeltaV;
-					bestTransfer = transfer;
 
-					bestTransfer.firstManoeuvre.time += startingTime;
+			let foundStartTime = startTime;
-					bestTransfer.secondManoeuvre.time += startingTime;
+			let foundEndTime = endTime;
+
+			let foundWorkingTransfer = false;
+			let trialCounter = 0;
+			while (!foundWorkingTransfer && trialCounter < 20) {
+				trialCounter++;
+
+				let testTransfer = interceptFunction(foundStartTime, foundEndTime);
+				if (testTransfer) {
+					foundWorkingTransfer = true;
+				} else {
+					foundStartTime = Math.max(0,startTime + Math.random() * 0.5 * startTimeStep);
+					foundEndTime = Math.max(0,endTime + Math.random() * 0.5 * endTimeStep);
+
+					if (foundEndTime < foundStartTime) {
+						foundEndTime += startTimeStep;
+					}
+				}
+			}
+
+			if (!foundWorkingTransfer) {
+				return;
+			}
+
+			let result = perform2dGradientDescent(interceptFunction, foundStartTime, foundEndTime, 0, 30);
+			if (result) {
+				let foundTransfer = getTransfer(result[0], result[1], backwards);
+				if (foundTransfer) {
+					let transferDeltaV = foundTransfer.firstManoeuvre.totalDeltaV + foundTransfer.secondManoeuvre.totalDeltaV;
+					if (bestDeltaV == null || transferDeltaV < bestDeltaV) {
+						bestDeltaV = transferDeltaV;
+						bestTransfer = foundTransfer;
+					}
 				}
 			}
 		});
@@ -1257,11 +1400,11 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 	}
 
 
-	const getPossibleTargetLocationFunction = (ownOrbit: Orbit, ownAltitude: number, targetAltitude: number, distanceToTarget: number): ((angle: number) => number[][]) => {
+	const getPossibleTargetVectors = (ownOrbit: Orbit, ownAltitude: number, targetAltitude: number, distanceToTarget: number): [number[][], number[][], number[][]] => {
 		let ownTrueAnomaly = getTrueAnomalyFromAltitude(ownAltitude, ownOrbit.semiLatusRectum, ownOrbit.eccentricity, ownShipHeadedInwards);
 
-		let ownLocalX = interpolationParameters.firstOwnAltitude * Math.cos(ownTrueAnomaly);
+		let ownLocalX = ownAltitude * Math.cos(ownTrueAnomaly);
-		let ownLocalY = interpolationParameters.firstOwnAltitude * Math.sin(ownTrueAnomaly);
+		let ownLocalY = ownAltitude * Math.sin(ownTrueAnomaly);
 
 		let ownDirection = normalizeVector(
 			addVector(
@@ -1270,8 +1413,8 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 			)
 		);
 
-		let perpendicularOne = getRandomPerpendicularVector(ownDirection);
+		let perpendicularInPlane = normalizeVector(vectorCrossProduct(ownDirection, ownOrbit.coordinateAxes[2]));
-		let perpendicularTwo = normalizeVector(vectorCrossProduct(ownDirection, perpendicularOne));
+		let perpendicularOutOfPlane = normalizeVector(vectorCrossProduct(ownDirection, perpendicularInPlane));
 
 		// Find phase angle to target
 		let phaseAngle = Math.acos((ownAltitude**2 + targetAltitude**2 - distanceToTarget**2) / (2 * ownAltitude * targetAltitude));
@@ -1280,95 +1423,173 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates
 		let distanceAlongNormal = targetAltitude * Math.sin(phaseAngle);
 
 		let root = multiplyMatrixWithScalar(distanceAlongDirection, ownDirection);
+		let vectorOne = multiplyMatrixWithScalar(distanceAlongNormal, perpendicularInPlane);
+		let vectorTwo = multiplyMatrixWithScalar(distanceAlongNormal, perpendicularOutOfPlane);
 
-		return function(angle: number) {
+		return [root, vectorOne, vectorTwo];
-			let addition = addVector(
-				multiplyMatrixWithScalar(distanceAlongNormal * Math.cos(angle), perpendicularOne),
-				multiplyMatrixWithScalar(distanceAlongNormal * Math.sin(angle), perpendicularTwo)
-			);
-
-			return addVector(root, addition);
-		}
 	}
 
-	let firstTargetFunction = getPossibleTargetLocationFunction(ownCoordinates.orbit, interpolationParameters.firstOwnAltitude, interpolationParameters.firstTargetAltitude, interpolationParameters.firstDistance);
+	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 secondTargetFunction = getPossibleTargetLocationFunction(ownCoordinates.orbit, interpolationParameters.secondOwnAltitude, interpolationParameters.secondTargetAltitude, interpolationParameters.secondDistance);
+		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 firstTargetTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
 	let secondTargetTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards);
 
-	// Use cosine rule
+	let [minimizeFunction, partialDerivativeAngleOne, partialDerivativeAngleTwo] = getFunctionAndDerivatives(firstPositionPossibleVectors, secondPositionPossibleVectors, interpolationParameters.firstTargetAltitude, interpolationParameters.secondTargetAltitude, secondTargetTrueAnomaly - firstTargetTrueAnomaly);
-	let targetMovement = Math.sqrt(interpolationParameters.firstTargetAltitude**2 + interpolationParameters.secondTargetAltitude**2 - 2 * interpolationParameters.firstTargetAltitude * interpolationParameters.secondTargetAltitude * Math.cos(secondTargetTrueAnomaly - firstTargetTrueAnomaly));
 
-	// Now, we need to do multivariate optimization to find the two angles to give to the functions that make the distance between
+	const gradientDescent = (startingGuessAngleOne: number, startingGuessAngleTwo: number): [number, number, number] => {
-	// two target positions as equal to the target movement as possible
+		let angleOne = startingGuessAngleOne;
-	const functionToOptimize = (angleOne: number, angleTwo: number) => {
+		let angleTwo = startingGuessAngleTwo;
-		let firstTargetPosition = firstTargetFunction(angleOne);
-		let secondTargetPosition = secondTargetFunction(angleTwo);
 
-		return Math.abs(getVectorMagnitude(subtractVector(secondTargetPosition, firstTargetPosition)) - targetMovement);
+		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 directionX = -dfDy / scaling;
+			let directionY = -dfDx / scaling;
+
+			let changeAlongDirection = directionX * dfDx + directionY * dfDy;
+
+			let deadEnd = false;
+			let trialValue = value;
+			let trialAngleOne = angleOne;
+			let trialAngleTwo = angleTwo;
+			let divisor = 1;
+
+			while (Math.abs(trialValue) >= Math.abs(value)) {
+				let stepAngleOne = - value * directionX / (changeAlongDirection * 2 ** divisor);
+				let stepAngleTwo = - value * directionY / (changeAlongDirection * 2 ** divisor);
+				trialAngleOne = angleOne + stepAngleOne;
+				trialAngleTwo = angleTwo + stepAngleTwo;
+				trialValue = minimizeFunction(trialAngleOne, trialAngleTwo);
+				divisor += 1;
+
+				if (divisor >= 32) {
+					deadEnd = true;
+					break;
+				}
+			}
+
+			if (deadEnd) {
+				return [angleOne, angleTwo, value];
+			}
+
+			value = trialValue;
+			angleOne = trialAngleOne;
+			angleTwo = trialAngleTwo;
+		}
+
+		return [angleOne, angleTwo, value];
 	}
 
 	let bestAngleOne = 0;
 	let bestAngleTwo = 0;
-	let bestValue = 1e99;
+	let bestValue = null;
 
-	for (var i = 0; i < 100; i++) {
+	// Try to start at a bunch of different positions and do gradient descent from there
-		for (var j = 0; j < 100; j++) {
+	for (var i = 0; i < 500; i++) {
-			let angleOne = 2 * Math.PI * i / 100;
+		for (var j = 0; j < 500; j++) {
-			let angleTwo = 2 * Math.PI * j / 100;
+			let startingAngleOne = i * 2 * Math.PI / 500.0;
+			let startingAngleTwo = j * 2 * Math.PI / 500.0
+			let startingValue = minimizeFunction(startingAngleOne, startingAngleTwo);
 
-			let value = functionToOptimize(angleOne, angleTwo);
+			if (bestValue == null || Math.abs(startingValue) < Math.abs(bestValue)) {
-			if (value < bestValue) {
+				bestAngleOne = startingAngleOne;
-				bestAngleOne = angleOne;
+				bestAngleTwo = startingAngleTwo;
-				bestAngleTwo = angleTwo;
+				bestValue = startingValue;
-				bestValue = value;
 			}
 		}
 	}
 
-	// Do some gradient descent on the best values found so far
+	[bestAngleOne, bestAngleOne, bestValue] = gradientDescent(bestAngleOne, bestAngleTwo);
-	while (bestValue > 0.5) {
-		let estimateDfDx = (functionToOptimize(bestAngleOne+0.0000001, bestAngleTwo) - bestValue) / 0.0000001;
-		let estimateDfDy = (functionToOptimize(bestAngleOne, bestAngleTwo+0.0000001) - bestValue) / 0.0000001;
 
-		let slopeNormal = vectorCrossProduct(
+	let bestFirstPosition = addVector(
-			[[1], [0], [estimateDfDx]],
+		firstPositionPossibleVectors[0],
-			[[0], [1], [estimateDfDy]]
+		addVector(
-		);
+			multiplyMatrixWithScalar(Math.cos(bestAngleOne), firstPositionPossibleVectors[1]),
-
+			multiplyMatrixWithScalar(Math.sin(bestAngleOne), firstPositionPossibleVectors[2])
-		let direction = normalizeVector([[slopeNormal[0][0]], [slopeNormal[1][0]], [0]]);
+		)
-		let estimateDfDv = (functionToOptimize(bestAngleOne + direction[0][0]*0.0001, bestAngleTwo + direction[1][0]*0.0001) - bestValue) / 0.0001;
+	);
-
+	let bestSecondPosition = addVector(
-		let bestValueImprovement = bestValue;
+		secondPositionPossibleVectors[0],
-		let bestAngleOneUpdate = bestAngleOne;
+		addVector(
-		let bestAngleTwoUpdate = bestAngleTwo;
+			multiplyMatrixWithScalar(Math.cos(bestAngleTwo), secondPositionPossibleVectors[1]),
-		let divisor = 0;
+			multiplyMatrixWithScalar(Math.sin(bestAngleTwo), secondPositionPossibleVectors[2])
-		let deadEnd = false;
+		)
-
+	);
-		while (bestValueImprovement >= bestValue) {
-			bestAngleOneUpdate = bestAngleOne - bestValue * direction[0][0] / (estimateDfDv * 2**divisor);
-			bestAngleTwoUpdate = bestAngleTwo - bestValue * direction[1][0] / (estimateDfDv * 2**divisor);
-			bestValueImprovement = functionToOptimize(bestAngleOneUpdate, bestAngleTwoUpdate);
-			divisor += 1;
-
-			if (divisor >= 32) {
-				deadEnd = true;
-				break;
-			}
-		}
-
-		if (deadEnd) {
-			break;
-		}
-		bestAngleOne = bestAngleOneUpdate;
-		bestAngleTwo = bestAngleTwoUpdate;
-		bestValue = bestValueImprovement;
-	}
-
-	let bestFirstPosition = firstTargetFunction(bestAngleOne);
-	let bestSecondPosition = secondTargetFunction(bestAngleTwo);
 
 	if (bestFirstPosition == null || bestSecondPosition == null) {
 		return null;

src/gui/worker.ts

@@ -46,7 +46,7 @@ export interface LandingProgressMessage {
 
 ctx.addEventListener("message", event => {
 	const progressCallback = (numberChecked: number, totalNumber: number, bestDeltaV: number | null, bestTransfer: Transfer | null) => {
-		if (numberChecked % 100 == 0) {
+		if (numberChecked % 1 == 0) {
 			let percentDone = numberChecked * 100 / totalNumber;
 			let message: ProgressMessage = {
 				type: "ProgressMessage",