diff --git a/src/calculations/orbit-calculations.ts b/src/calculations/orbit-calculations.ts index 8ff68ba..4576098 100644 --- a/src/calculations/orbit-calculations.ts +++ b/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]; - } - - 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]; + 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); - // 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); + let targetFirstTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.firstTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards); + let targetSecondTrueAnomaly = getTrueAnomalyFromAltitude(interpolationParameters.secondTargetAltitude, targetSemiLatusRectum, targetEccentricity, targetHeadedInwards); - let [minimizeFunction, partialDerivativeAngleOne, partialDerivativeAngleTwo] = getFunctionAndDerivatives(firstPositionPossibleVectors, secondPositionPossibleVectors, interpolationParameters.firstTargetAltitude, interpolationParameters.secondTargetAltitude, secondTargetTrueAnomaly - firstTargetTrueAnomaly); + 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); + } + 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 directionX = -dfDy / scaling; - let directionY = -dfDx / scaling; + let phaseAngleValue = distanceAlongPerpendicular * Math.cos(angleOne) / distanceAlongDirection - Math.tan(interpolationParameters.firstPhaseAngle); + let phaseAngleDerivative = -distanceAlongPerpendicular * Math.sin(angleOne) / distanceAlongDirection; - 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); - let stepAngleTwo = - value * directionY / (changeAlongDirection * 2 ** divisor); - trialAngleOne = angleOne + stepAngleOne; - trialAngleTwo = angleTwo + stepAngleTwo; + trialAngleOne = angleOne - changeAngleOne / (2 ** divisor); + trialAngleTwo = angleTwo - changeAngleTwo / (2 ** divisor); 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 - for (var i = 0; i < 500; i++) { - for (var j = 0; j < 500; j++) { - let startingAngleOne = i * 2 * Math.PI / 500.0; - let startingAngleTwo = j * 2 * Math.PI / 500.0 - let startingValue = minimizeFunction(startingAngleOne, startingAngleTwo); - - if (bestValue == null || Math.abs(startingValue) < Math.abs(bestValue)) { - bestAngleOne = startingAngleOne; - bestAngleTwo = startingAngleTwo; - bestValue = startingValue; + [angleOneGuess, -angleOneGuess].forEach(angleOne => { + [angleTwoGuess, -angleTwoGuess].forEach(angleTwo => { + let value = minimizeFunction(angleOne, angleTwo); + if (bestValue == null || value < bestValue) { + bestValue = value; + bestAngleOne = angleOne; + bestAngleTwo = angleTwo; } - } - } + }); + }); [bestAngleOne, bestAngleOne, bestValue] = gradientDescent(bestAngleOne, bestAngleTwo); diff --git a/src/gui/interpolate.ts b/src/gui/interpolate.ts index aa0eb4f..2fdd231 100644 --- a/src/gui/interpolate.ts +++ b/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(); }); } } \ No newline at end of file