From b3366eba010549e1eb075a70aa49599de1c5d921 Mon Sep 17 00:00:00 2001 From: Martin Asprusten Date: Sat, 4 Apr 2026 11:17:30 +0200 Subject: [PATCH] Simplified interpolation --- index.html | 2 +- src/calculations/mathematics.ts | 43 +++ src/calculations/orbit-calculations.ts | 494 ++++++++++++++----------- src/gui/interpolate.ts | 17 + 4 files changed, 348 insertions(+), 208 deletions(-) diff --git a/index.html b/index.html index abc4a5d..2781fca 100644 --- a/index.html +++ b/index.html @@ -1,5 +1,5 @@ - + diff --git a/src/calculations/mathematics.ts b/src/calculations/mathematics.ts index 4dd7751..ed21150 100644 --- a/src/calculations/mathematics.ts +++ b/src/calculations/mathematics.ts @@ -133,6 +133,20 @@ export function addVector(vectorOne: number[][], vectorTwo: number[][]): number[ return result; } +export function addVectors(vectors: number[][][]) { + let resultVector = [ + [0], + [0], + [0] + ]; + + vectors.forEach(vector => { + resultVector = addVector(resultVector, vector); + }); + + return resultVector; +} + export function addMatrix(matrixOne: number[][], matrixTwo: number[][]): number[][] { if (!checkIfValidMatrix(matrixOne) || !checkIfValidMatrix(matrixTwo)) { throw new TypeError("Two valid matrices are required"); @@ -226,4 +240,33 @@ export function multiplyMatrixWithScalar(scalar: number, matrix: number[][]): nu } return result; +} + +export function invertTwoByTwoMatrix(matrix: number[][]): number[][] { + if (!checkIfValidMatrix(matrix)) { + throw new TypeError("A valid matrix is required"); + } + + if (matrix.length != 2) { + throw new TypeError("Can only invert 2x2 matrices"); + } + + if (matrix[0].length != 2) { + throw new TypeError("Can only invert 2x2 matrices"); + } + + let divisor = matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]; + + let multiplier = 1 / divisor; + return [ + [ + matrix[1][1] * multiplier, + -matrix[0][1] * multiplier + ], + [ + -matrix[1][0] * multiplier, + matrix[0][0] * multiplier + ] + + ]; } \ No newline at end of file diff --git a/src/calculations/orbit-calculations.ts b/src/calculations/orbit-calculations.ts index 4576098..d64315f 100644 --- a/src/calculations/orbit-calculations.ts +++ b/src/calculations/orbit-calculations.ts @@ -1,6 +1,6 @@ import type { InterpolationParameters } from "../gui/interpolate"; import type { Body } from "./constants"; -import { addMatrix, addVector, getVectorMagnitude, matrixMultiply, multiplyMatrixWithScalar, normalizeVector, subtractVector, vectorCrossProduct, vectorDotProduct } from "./mathematics"; +import { addMatrix, addVector, getVectorMagnitude, invertTwoByTwoMatrix, matrixMultiply, multiplyMatrixWithScalar, normalizeVector, subtractVector, vectorCrossProduct, vectorDotProduct } from "./mathematics"; export interface Orbit { semiLatusRectum: number, @@ -362,7 +362,14 @@ export function getEccentricAndTrueAnomalyFromMeanAnomaly(meanAnomaly: number, e } export function getOrbitalCoordinates(timeToPeriapsis: number, orbit: Orbit, planet: Body): OrbitalCoordinates { - const meanAngularMotion = Math.sqrt(planet.gravitationalParameter * Math.abs(1 - orbit.eccentricity**2)**3 / orbit.semiLatusRectum**3); + let meanAngularMotion; + + if (Math.abs(orbit.eccentricity - 1) < 0.0001) { + meanAngularMotion = Math.sqrt(8 * planet.gravitationalParameter / orbit.semiLatusRectum**3); + } else { + meanAngularMotion = Math.sqrt(planet.gravitationalParameter * Math.abs(1 - orbit.eccentricity**2)**3 / orbit.semiLatusRectum**3); + } + const meanAnomaly = meanAngularMotion * -timeToPeriapsis; @@ -412,9 +419,11 @@ export function getOrbitalCoordinatesFromAltitude(altitude: number, headingInwar let eccentricAnomaly; let meanAnomaly; + let meanAngularMotion: number | null = null; if (Math.abs(orbit.eccentricity - 1) < 0.0001) { eccentricAnomaly = trueAnomaly; - meanAnomaly = eccentricAnomaly; + meanAnomaly = Math.sqrt(2) * (Math.tan(trueAnomaly / 2) + Math.tan(trueAnomaly / 2)**3 / 3); + meanAngularMotion = Math.sqrt(8 * planet.gravitationalParameter / orbit.semiLatusRectum**3); } else if (orbit.eccentricity < 1) { eccentricAnomaly = Math.atan2(Math.sqrt(1 - orbit.eccentricity**2)*Math.sin(trueAnomaly), orbit.eccentricity + Math.cos(trueAnomaly)); meanAnomaly = eccentricAnomaly - orbit.eccentricity * Math.sin(eccentricAnomaly); @@ -423,7 +432,9 @@ export function getOrbitalCoordinatesFromAltitude(altitude: number, headingInwar meanAnomaly = orbit.eccentricity * Math.sinh(eccentricAnomaly) - eccentricAnomaly; } - const meanAngularMotion = Math.sqrt(planet.gravitationalParameter * Math.abs(1 - orbit.eccentricity**2)**3 / orbit.semiLatusRectum**3); + if (meanAngularMotion == null) { + meanAngularMotion = Math.sqrt(planet.gravitationalParameter * Math.abs(1 - orbit.eccentricity**2)**3 / orbit.semiLatusRectum**3); + } return { orbit: orbit, @@ -977,6 +988,10 @@ export function findCheapestIntercept(startingSituation: OrbitalCoordinates, tar let startOrbitMeanAnomaly = startingSituation.meanAnomaly + startTime * startingSituation.meanAngularMotion; let endOrbitMeanAnomaly = targetSituation.meanAnomaly + endTime * targetSituation.meanAngularMotion; + if (Math.abs(startingSituation.orbit.eccentricity - 1) < 0.0001) { + + } + let extraStartOrbits = Math.floor(startOrbitMeanAnomaly / (2 * Math.PI)); let extraEndOrbits = Math.floor(endOrbitMeanAnomaly / (2 * Math.PI)); @@ -1141,7 +1156,7 @@ export function findCheapestLanding(startingCoordinates: OrbitalCoordinates, sta let position = startingPosition; let mass = ship.mass; - let massUse = ship.thrust / ship.specificImpulse; + let massUse = ship.thrust / (ship.specificImpulse * 9.81); let timeSteps = 0; let lengthOfStep = 0.01; @@ -1368,262 +1383,327 @@ export function findOrbitThroughInterpolation(ownCoordinates: OrbitalCoordinates } return trueAnomaly; } + + const ownFirstTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards); + const ownSecondTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards); + const targetFirstTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards); + const targetSecondTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards); - const getPossibleTargetVectors = (ownOrbit: Orbit, ownAltitude: number, targetAltitude: number, distanceToTarget: number): [number[][], number[][], number[][], number, number] => { - let ownTrueAnomaly = getTrueAnomalyFromAltitude(ownAltitude, ownOrbit.semiLatusRectum, ownOrbit.eccentricity, ownShipHeadedInwards); + const getVectors = (trueAnomaly: number, orbit: Orbit) => { + const radius = orbit.semiLatusRectum / (1 + orbit.eccentricity * Math.cos(trueAnomaly)); + const localX = radius * Math.cos(trueAnomaly); + const localY = radius * Math.sin(trueAnomaly); - let ownLocalX = ownAltitude * Math.cos(ownTrueAnomaly); - let ownLocalY = ownAltitude * Math.sin(ownTrueAnomaly); - - let ownDirection = normalizeVector( + const vectorPointingTowardsShip = normalizeVector( addVector( - multiplyMatrixWithScalar(ownLocalX, ownOrbit.coordinateAxes[0]), - multiplyMatrixWithScalar(ownLocalY, ownOrbit.coordinateAxes[1]) + multiplyMatrixWithScalar(localX, orbit.coordinateAxes[0]), + multiplyMatrixWithScalar(localY, orbit.coordinateAxes[1]) ) ); - let perpendicularInPlane = normalizeVector(vectorCrossProduct(ownOrbit.coordinateAxes[2], ownDirection)); - let perpendicularOutOfPlane = normalizeVector(vectorCrossProduct(ownDirection, perpendicularInPlane)); - - // Find phase angle to target - let phaseAngle = Math.acos((ownAltitude**2 + targetAltitude**2 - distanceToTarget**2) / (2 * ownAltitude * targetAltitude)); - - let distanceAlongDirection = targetAltitude * Math.cos(phaseAngle); - let distanceAlongNormal = targetAltitude * Math.sin(phaseAngle); - - let root = multiplyMatrixWithScalar(distanceAlongDirection, ownDirection); - let vectorOne = multiplyMatrixWithScalar(distanceAlongNormal, perpendicularInPlane); - let vectorTwo = multiplyMatrixWithScalar(distanceAlongNormal, perpendicularOutOfPlane); - - return [root, vectorOne, vectorTwo, distanceAlongDirection, distanceAlongNormal]; - } - - let firstPositionPossibleVectors = getPossibleTargetVectors(ownCoordinates.orbit, interpolationParameters.firstOwnAltitude, interpolationParameters.firstTargetAltitude, interpolationParameters.firstDistance); - let secondPositionPossibleVectors = getPossibleTargetVectors(ownCoordinates.orbit, interpolationParameters.secondOwnAltitude, interpolationParameters.secondTargetAltitude, interpolationParameters.secondDistance); - - let targetFirstTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards); - let targetSecondTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards); - - let c1 = firstPositionPossibleVectors[0]; - let p11 = firstPositionPossibleVectors[1]; - let p12 = firstPositionPossibleVectors[2]; - - let c2 = secondPositionPossibleVectors[0]; - let p21 = secondPositionPossibleVectors[1]; - let p22 = secondPositionPossibleVectors[2]; - - // Scale everything to avoid numerical explosions - let scaleFactor = 1 / getVectorMagnitude(c2); - - c1 = multiplyMatrixWithScalar(scaleFactor, c1); - p11 = multiplyMatrixWithScalar(scaleFactor, p11); - p12 = multiplyMatrixWithScalar(scaleFactor, p12); - c2 = multiplyMatrixWithScalar(scaleFactor, c2); - p21 = multiplyMatrixWithScalar(scaleFactor, p21); - p22 = multiplyMatrixWithScalar(scaleFactor, p22); - - let scaledAltitudeOne = interpolationParameters.firstTargetAltitude * scaleFactor; - let scaledAltitudeTwo = interpolationParameters.secondTargetAltitude * scaleFactor; - - const constantPart = scaledAltitudeOne * scaledAltitudeTwo * Math.cos(targetSecondTrueAnomaly - targetFirstTrueAnomaly); - - const functionToMinimize = (angleOne: number, angleTwo: number) => { - let d1 = addVector(c1, addVector( - multiplyMatrixWithScalar(Math.cos(angleOne), p11), - multiplyMatrixWithScalar(Math.sin(angleOne), p12) - )); - - let d2 = addVector(c2, addVector( - multiplyMatrixWithScalar(Math.cos(angleTwo), p21), - multiplyMatrixWithScalar(Math.sin(angleTwo), p22) - )); - - return vectorDotProduct(d1, d2) - constantPart; - } - - const partialDerivativeAngleOne = (angleOne: number, angleTwo: number) => { - let d2 = addVector(c2, addVector( - multiplyMatrixWithScalar(Math.cos(angleTwo), p21), - multiplyMatrixWithScalar(Math.sin(angleTwo), p22) - )); - - let d1Dx = addVector( - multiplyMatrixWithScalar(-Math.sin(angleOne), p11), - multiplyMatrixWithScalar(Math.cos(angleOne), p12) + const perpendicularVectorInPlane = normalizeVector( + vectorCrossProduct(orbit.coordinateAxes[2], vectorPointingTowardsShip) ); - return vectorDotProduct(d1Dx, d2); + const perpendicularVectorOutOfPlane = orbit.coordinateAxes[2]; + + return [vectorPointingTowardsShip, perpendicularVectorInPlane, perpendicularVectorOutOfPlane]; } - const partialDerivativeAngleTwo = (angleOne: number, angleTwo: number) => { - let d1 = addVector(c1, addVector( - multiplyMatrixWithScalar(Math.cos(angleOne), p11), - multiplyMatrixWithScalar(Math.sin(angleOne), p12) - )); + const firstVectors = getVectors(ownFirstTrueAnomaly, ownCoordinates.orbit); + const secondVectors = getVectors(ownSecondTrueAnomaly, ownCoordinates.orbit); - let d2Dy = addVector( - multiplyMatrixWithScalar(-Math.sin(angleTwo), p21), - multiplyMatrixWithScalar(Math.cos(angleTwo), p22) + const getDistances = (ownAltitude: number, targetAltitude: number, targetDistance: number) => { + const angleBetweenPlanetAndTarget = Math.acos((ownAltitude**2 + targetAltitude**2 - targetDistance**2) / (2*ownAltitude*targetAltitude)); + const distanceAlongDirection = Math.cos(angleBetweenPlanetAndTarget) * targetAltitude; + const distanceAlongPerpendicular = Math.sin(angleBetweenPlanetAndTarget) * targetAltitude; + + return [distanceAlongDirection, distanceAlongPerpendicular]; + } + + const [firstDistanceAlongDirection, firstDistancePerpendicularToDirection] = getDistances(interpolationParameters.firstOwnAltitude, interpolationParameters.firstTargetAltitude, interpolationParameters.firstDistance); + const [secondDistanceAlongDirection, secondDistancePerpendicularToDirection] = getDistances(interpolationParameters.secondOwnAltitude, interpolationParameters.secondTargetAltitude, interpolationParameters.secondDistance); + + const phaseAngleDifference = (ownSecondTrueAnomaly + interpolationParameters.secondPhaseAngle) - (ownFirstTrueAnomaly + interpolationParameters.firstPhaseAngle); + + // Now, try some Newton's method to estimate the two position's the target has gone through + + // To avoid the maths exploding, we'll scale everything down a bit + const scaleFactor = 1; + + const c1 = multiplyMatrixWithScalar(firstDistanceAlongDirection * scaleFactor, firstVectors[0]); + const n11 = multiplyMatrixWithScalar(firstDistancePerpendicularToDirection * scaleFactor, firstVectors[1]); + const n12 = multiplyMatrixWithScalar(firstDistancePerpendicularToDirection * scaleFactor, firstVectors[2]); + + const c2 = multiplyMatrixWithScalar(secondDistanceAlongDirection * scaleFactor, secondVectors[0]); + const n21 = multiplyMatrixWithScalar(secondDistancePerpendicularToDirection * scaleFactor, secondVectors[1]); + const n22 = multiplyMatrixWithScalar(secondDistancePerpendicularToDirection * scaleFactor, secondVectors[2]); + + let expectedDistance = Math.sqrt(interpolationParameters.firstTargetAltitude**2 + interpolationParameters.secondTargetAltitude**2 - 2 * interpolationParameters.firstTargetAltitude * interpolationParameters.secondTargetAltitude * Math.cos(targetSecondTrueAnomaly - targetFirstTrueAnomaly)); + // Scale the expected distance too + expectedDistance *= scaleFactor; + + const distanceFunction = (angleOne: number, angleTwo: number) => { + let p1 = addVector( + c1, + addVector( + multiplyMatrixWithScalar(Math.cos(angleOne), n11), + multiplyMatrixWithScalar(Math.sin(angleOne), n12) + ) ); - return vectorDotProduct(d1, d2Dy); + let p2 = addVector( + c2, + addVector( + multiplyMatrixWithScalar(Math.cos(angleTwo), n21), + multiplyMatrixWithScalar(Math.sin(angleTwo), n22) + ) + ); + + return getVectorMagnitude(addVector(p2, multiplyMatrixWithScalar(-1, p1))) - expectedDistance; } - let firstOwnTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards); - 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 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 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); + 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 gradientDescent = (startingGuessAngleOne: number, startingGuessAngleTwo: number): [number, number, number] => { - let angleOne = startingGuessAngleOne; - let angleTwo = startingGuessAngleTwo; + const planeAngleFunction = (angleOne: number, angleTwo: number) => { + const horizontalOne = addVector( + c1, + multiplyMatrixWithScalar(Math.cos(angleOne), n11) + ); - let distanceAlongDirection = getVectorMagnitude(firstPositionPossibleVectors[0]); - let distanceAlongPerpendicular = getVectorMagnitude(firstPositionPossibleVectors[1]); + const horizontalTwo = addVector( + c2, + multiplyMatrixWithScalar(Math.cos(angleTwo), n21) + ); - 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); + return vectorDotProduct(horizontalOne, horizontalTwo) / (getVectorMagnitude(horizontalOne) * getVectorMagnitude(horizontalTwo)) - Math.cos(phaseAngleDifference); + } - let phaseAngleValue = distanceAlongPerpendicular * Math.cos(angleOne) / distanceAlongDirection - Math.tan(interpolationParameters.firstPhaseAngle); - let phaseAngleDerivative = -distanceAlongPerpendicular * Math.sin(angleOne) / distanceAlongDirection; + const planeAnglePartialDerivativeAngleOne = (angleOne: number, angleTwo: number) => { + const horizontalOne = addVector( + c1, + multiplyMatrixWithScalar(Math.cos(angleOne), n11) + ); - let changeAngleOne = dfDx * phaseAngleValue / (dfDy * phaseAngleDerivative); - let changeAngleTwo = (dfDx * phaseAngleValue - phaseAngleDerivative * value) / (dfDy * phaseAngleDerivative); + const horizontalTwo = addVector( + c2, + multiplyMatrixWithScalar(Math.cos(angleTwo), n21) + ); - let trialAngleOne = angleOne; - let trialAngleTwo = angleTwo; - let trialValue = value; - let divisor = 1; + 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)); + } - // We don't want things to blow up - while (Math.abs(changeAngleOne) / 2**divisor > Math.PI / 1000 || Math.abs(changeAngleTwo) / 2**divisor > Math.PI / 1000) { - divisor +=1; + const 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(trialValue) >= Math.abs(value)) { - trialAngleOne = angleOne - changeAngleOne / (2 ** divisor); - trialAngleTwo = angleTwo - changeAngleTwo / (2 ** divisor); - trialValue = minimizeFunction(trialAngleOne, trialAngleTwo); - divisor += 1; - - if (divisor >= 32) { + 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) { - return [angleOne, angleTwo, value]; - } + break; + }*/ + anglesVector = addVector(anglesVector, multiplyMatrixWithScalar(-1, update)); + }; + return anglesVector; + } - value = trialValue; - angleOne = trialAngleOne; - angleTwo = trialAngleTwo; + 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); } - return [angleOne, angleTwo, value]; + if (Math.abs(interpolationParameters.secondPhaseAngle) < Math.PI / 2) { + angleTwo = Math.acos(Math.tan(secondPhaseAngle) * secondDistanceAlongDirection / secondDistancePerpendicularToDirection); + } else { + angleTwo = Math.acos(Math.tan(Math.PI - secondPhaseAngle) * -secondDistanceAlongDirection / secondDistancePerpendicularToDirection); + } + + return [angleOne, angleTwo]; } - let firstCenterLength = getVectorMagnitude(firstPositionPossibleVectors[0]); - let firstPerpendicularLength = getVectorMagnitude(firstPositionPossibleVectors[1]); + let bestAngles = [[0], [0]]; + let bestDistance: number | null = null; - let secondCenterLength = getVectorMagnitude(secondPositionPossibleVectors[0]); - let secondPerpendicularLength = getVectorMagnitude(secondPositionPossibleVectors[1]); + 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 angleOneGuess = Math.acos(firstCenterLength * Math.tan(interpolationParameters.firstPhaseAngle) / firstPerpendicularLength); - let angleTwoGuess = Math.acos(secondCenterLength * Math.tan(interpolationParameters.secondPhaseAngle) / secondPerpendicularLength); - - let bestAngleOne = 0; - let bestAngleTwo = 0; - let bestValue: number | null = null; - - [angleOneGuess, -angleOneGuess].forEach(angleOne => { - [angleTwoGuess, -angleTwoGuess].forEach(angleTwo => { - let value = minimizeFunction(angleOne, angleTwo); - if (bestValue == null || value < bestValue) { - bestValue = value; - bestAngleOne = angleOne; - bestAngleTwo = angleTwo; + 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]]; } - }); - }); - - [bestAngleOne, bestAngleOne, bestValue] = gradientDescent(bestAngleOne, bestAngleTwo); - - let bestFirstPosition = addVector( - firstPositionPossibleVectors[0], - addVector( - multiplyMatrixWithScalar(Math.cos(bestAngleOne), firstPositionPossibleVectors[1]), - multiplyMatrixWithScalar(Math.sin(bestAngleOne), firstPositionPossibleVectors[2]) - ) - ); - let bestSecondPosition = addVector( - secondPositionPossibleVectors[0], - addVector( - multiplyMatrixWithScalar(Math.cos(bestAngleTwo), secondPositionPossibleVectors[1]), - multiplyMatrixWithScalar(Math.sin(bestAngleTwo), secondPositionPossibleVectors[2]) - ) - ); - - if (bestFirstPosition == null || bestSecondPosition == null) { - return null; + } } - // We can now find the vector that defines the plane the target is in - let normalVector = normalizeVector(vectorCrossProduct(bestFirstPosition, bestSecondPosition)); - let firstTargetPositionDirection = normalizeVector(bestFirstPosition); - // Rotate the position vector about the normal vector by the true anomaly to find the vector pointing to the periapsis - let identityMatrix = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]; - let crossProductMatrix = [ - [0, -normalVector[2][0], normalVector[1][0]], - [normalVector[2][0], 0, -normalVector[0][0]], - [-normalVector[1][0], normalVector[0][0], 0] - ]; - let outerProductMatrix = [ - [normalVector[0][0]**2, normalVector[0][0]*normalVector[1][0], normalVector[0][0]*normalVector[2][0]], - [normalVector[0][0]*normalVector[1][0], normalVector[1][0]**2, normalVector[1][0]*normalVector[2][0]], - [normalVector[0][0]*normalVector[2][0], normalVector[1][0]*normalVector[2][0], normalVector[2][0]**2] - ]; + if (!interpolationParameters.targetAbovePlane) { + bestAngles = [[-bestAngles[0][0]], [-bestAngles[1][0]]]; + } - let rotationMatrix = addMatrix( - multiplyMatrixWithScalar(Math.cos(-firstTargetTrueAnomaly), identityMatrix), - addMatrix( - multiplyMatrixWithScalar(Math.sin(-firstTargetTrueAnomaly), crossProductMatrix), - multiplyMatrixWithScalar(1 - Math.cos(-firstTargetTrueAnomaly), outerProductMatrix) + // 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 targetOrbitXVector = matrixMultiply(rotationMatrix, firstTargetPositionDirection); - let targetOrbitYVector = vectorCrossProduct(normalVector, targetOrbitXVector); + 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: [targetOrbitXVector, targetOrbitYVector, normalVector] + coordinateAxes: [ + periapsisVector, + normalizeVector(vectorCrossProduct(normalVector, periapsisVector)), + normalVector + ] }; - let timePassed = getTimeBetweenTrueAnomalies(ownCoordinates.trueAnomaly, getTrueAnomalyFromAltitude(interpolationParameters.secondOwnAltitude, ownCoordinates.orbit.semiLatusRectum, ownCoordinates.orbit.eccentricity, ownShipHeadedInwards), ownCoordinates.orbit, body); + let targetCoordinates: OrbitalCoordinates = getOrbitalCoordinatesFromAltitude( + interpolationParameters.secondTargetAltitude, + targetHeadedInwards, + targetOrbit, + body + ); - return [getOrbitalCoordinatesFromAltitude(interpolationParameters.secondTargetAltitude, targetHeadedInwards, targetOrbit, body), timePassed]; + let epochTrueAnomaly = ownCoordinates.trueAnomaly; + while (epochTrueAnomaly + 2 * Math.PI < ownFirstTrueAnomaly) { + epochTrueAnomaly += 2 * Math.PI; + } + + const timeElapsed = getTimeBetweenTrueAnomalies(epochTrueAnomaly, ownSecondTrueAnomaly, ownCoordinates.orbit, body); + return [ + targetCoordinates, + timeElapsed + ] } \ No newline at end of file diff --git a/src/gui/interpolate.ts b/src/gui/interpolate.ts index 2fdd231..810b124 100644 --- a/src/gui/interpolate.ts +++ b/src/gui/interpolate.ts @@ -7,6 +7,7 @@ export interface InterpolationParameters { targetApoapsis: number, targetSpeed: number, targetAltitude: number, + targetAbovePlane: boolean, firstOwnAltitude: number, firstTargetAltitude: number, firstDistance: number, @@ -23,6 +24,7 @@ export const DefaultInterpolationParameters: InterpolationParameters = { targetApoapsis: 200000, targetSpeed: 0, targetAltitude: 0, + targetAbovePlane: true, firstOwnAltitude: 0, firstTargetAltitude: 0, firstDistance: 0, @@ -78,6 +80,13 @@ 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); @@ -102,6 +111,7 @@ export class InterpolateOrbitGui { 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); @@ -122,6 +132,7 @@ export class InterpolateOrbitGui { targetApoapsis: apoapsis, targetSpeed: speed, targetAltitude: altitude, + targetAbovePlane: abovePlane, firstTargetAltitude: firstTargetAltitude, firstOwnAltitude: firstOwnAltitude, firstDistance: firstDistance, @@ -140,6 +151,7 @@ export class InterpolateOrbitGui { apoapsisInput, speedInput, altitudeInput, + abovePlaneButton, firstOwnAltitudeInput, firstTargetAltitudeInput, firstDistanceInput, @@ -171,6 +183,11 @@ 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();