diff --git a/examples/conjunction-assessment-example.ts b/examples/conjunction-assessment-example.ts new file mode 100644 index 00000000..acb81bd3 --- /dev/null +++ b/examples/conjunction-assessment-example.ts @@ -0,0 +1,247 @@ +/** + * @file Conjunction Assessment Example + * @description Demonstrates high accuracy conjunction assessment using + * historical TLE accuracy and covariance propagation. + * + * This example shows how to: + * 1. Set up a conjunction assessment between two space objects + * 2. Use high-fidelity propagators for improved accuracy + * 3. Propagate covariance matrices based on TLE quality + * 4. Calculate probability of collision + * + * @author Theodore Kruczek + * @license AGPL-3.0-or-later + * @copyright (c) 2025 Kruczek Labs LLC + */ + +import { + ConjunctionAssessment, + EpochUTC, + ForceModel, + Kilometers, + Tle, + StateCovariance, + CovarianceFrame, +} from '../src/main'; + +// Example 1: Basic Conjunction Assessment using TLEs +function basicConjunctionAssessment() { + console.log('=== Example 1: Basic Conjunction Assessment ===\n'); + + // Sample TLE data for two satellites in close proximity + const primaryTle = new Tle( + '1 25544U 98067A 25019.50000000 .00016717 00000-0 10270-3 0 9005', + '2 25544 51.6400 339.8000 0002571 90.5000 269.6000 15.50000000000000', + ); + + const secondaryTle = new Tle( + '1 44691U 19074A 25019.50000000 .00016500 00000-0 10200-3 0 9006', + '2 44691 51.6450 339.8050 0002600 90.5050 269.6050 15.50005000000000', + ); + + // Create conjunction assessment with object radii + const assessment = new ConjunctionAssessment( + { + name: 'ISS (Zarya)', + tle: primaryTle, + radius: 0.05 as Kilometers, // 50 meters + }, + { + name: 'Secondary Object', + tle: secondaryTle, + radius: 0.01 as Kilometers, // 10 meters + }, + ); + + // Define search window (6 hours) + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T18:00:00.000Z'); + + // Perform conjunction assessment + const event = assessment.assess({ + startTime, + endTime, + }); + + // Display results + console.log(event.toString()); + console.log(`\nHigh Risk: ${event.isHighRisk(1.0 as Kilometers)}`); +} + +// Example 2: High-Fidelity Propagation with Covariance +function highFidelityConjunctionAssessment() { + console.log('\n=== Example 2: High-Fidelity Assessment with Covariance ===\n'); + + const primaryTle = new Tle( + '1 25544U 98067A 25019.50000000 .00016717 00000-0 10270-3 0 9005', + '2 25544 51.6400 339.8000 0002571 90.5000 269.6000 15.50000000000000', + ); + + const secondaryTle = new Tle( + '1 44691U 19074A 25019.50000000 .00016500 00000-0 10200-3 0 9006', + '2 44691 51.6450 339.8050 0002600 90.5050 269.6050 15.50005000000000', + ); + + const assessment = new ConjunctionAssessment( + { + name: 'ISS (Zarya)', + tle: primaryTle, + radius: 0.05 as Kilometers, + }, + { + name: 'Secondary Object', + tle: secondaryTle, + radius: 0.01 as Kilometers, + }, + ); + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T18:00:00.000Z'); + + // Use high-fidelity propagation with force model + const forceModel = new ForceModel().setGravity(8, 8).setAtmosphericDrag().setSolarRadiationPressure(); + + const event = assessment.assess({ + startTime, + endTime, + useHighFidelity: true, + forceModel, + propagateCovariance: true, // Propagate TLE-based covariances + }); + + console.log(event.toString()); + + if (event.probabilityOfCollision !== undefined) { + console.log(`\nProbability of Collision: ${event.probabilityOfCollision.toExponential(6)}`); + } + + if (event.getMahalanobisDistance() !== undefined) { + console.log(`Mahalanobis Distance: ${event.getMahalanobisDistance()!.toFixed(3)} sigma`); + } +} + +// Example 3: Custom Covariance Matrices +function customCovarianceAssessment() { + console.log('\n=== Example 3: Custom Covariance Matrices ===\n'); + + const primaryTle = new Tle( + '1 25544U 98067A 25019.50000000 .00016717 00000-0 10270-3 0 9005', + '2 25544 51.6400 339.8000 0002571 90.5000 269.6000 15.50000000000000', + ); + + const secondaryTle = new Tle( + '1 44691U 19074A 25019.50000000 .00016500 00000-0 10200-3 0 9006', + '2 44691 51.6450 339.8050 0002600 90.5050 269.6050 15.50005000000000', + ); + + // Define custom covariances (1-sigma values in RIC frame) + // [radial, intrack, crosstrack, radial_vel, intrack_vel, crosstrack_vel] + const primaryCovariance = StateCovariance.fromSigmas( + [ + 0.5, // 500 m radial uncertainty + 1.5, // 1.5 km intrack uncertainty + 0.5, // 500 m crosstrack uncertainty + 0.001, // 1 m/s radial velocity uncertainty + 0.003, // 3 m/s intrack velocity uncertainty + 0.001, // 1 m/s crosstrack velocity uncertainty + ], + CovarianceFrame.RIC, + ); + + const secondaryCovariance = StateCovariance.fromSigmas( + [0.3, 1.0, 0.3, 0.0005, 0.002, 0.0005], + CovarianceFrame.RIC, + ); + + const assessment = new ConjunctionAssessment( + { + name: 'ISS (Zarya)', + tle: primaryTle, + covariance: primaryCovariance, + radius: 0.05 as Kilometers, + }, + { + name: 'Secondary Object', + tle: secondaryTle, + covariance: secondaryCovariance, + radius: 0.01 as Kilometers, + }, + ); + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T18:00:00.000Z'); + + const event = assessment.assess({ + startTime, + endTime, + }); + + console.log(event.toString()); +} + +// Example 4: Screening Multiple Objects +function screeningExample() { + console.log('\n=== Example 4: Multi-Object Screening ===\n'); + + const primaryTle = new Tle( + '1 25544U 98067A 25019.50000000 .00016717 00000-0 10270-3 0 9005', + '2 25544 51.6400 339.8000 0002571 90.5000 269.6000 15.50000000000000', + ); + + // List of potential conjunction objects + const secondaryTles = [ + new Tle( + '1 44691U 19074A 25019.50000000 .00016500 00000-0 10200-3 0 9006', + '2 44691 51.6450 339.8050 0002600 90.5050 269.6050 15.50005000000000', + ), + new Tle( + '1 12345U 81001A 25019.50000000 .00016400 00000-0 10100-3 0 9007', + '2 12345 51.6500 339.8100 0002650 90.5100 269.6100 15.50010000000000', + ), + ]; + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T18:00:00.000Z'); + + const screeningThreshold = 5.0 as Kilometers; // 5 km screening threshold + const pcThreshold = 1e-6; // Pc > 1e-6 is concerning + + console.log(`Screening ${secondaryTles.length} objects for conjunctions...\n`); + + secondaryTles.forEach((secondaryTle, index) => { + const assessment = new ConjunctionAssessment( + { tle: primaryTle, radius: 0.05 as Kilometers }, + { tle: secondaryTle, radius: 0.01 as Kilometers }, + ); + + const event = assessment.assess({ + startTime, + endTime, + useHighFidelity: true, + propagateCovariance: true, + }); + + if (event.missDistance < screeningThreshold) { + console.log(`Object ${index + 1}: CLOSE APPROACH DETECTED`); + console.log(` TCA: ${event.tca.toISOString()}`); + console.log(` Miss Distance: ${event.missDistance.toFixed(3)} km`); + + if (event.probabilityOfCollision !== undefined && event.probabilityOfCollision > pcThreshold) { + console.log(` Pc: ${event.probabilityOfCollision.toExponential(3)} [HIGH RISK]`); + } else if (event.probabilityOfCollision !== undefined) { + console.log(` Pc: ${event.probabilityOfCollision.toExponential(3)}`); + } + console.log(); + } + }); +} + +// Run all examples +if (require.main === module) { + basicConjunctionAssessment(); + highFidelityConjunctionAssessment(); + customCovarianceAssessment(); + screeningExample(); + + console.log('\n=== All Examples Complete ==='); +} diff --git a/src/conjunction/ConjunctionAssessment.ts b/src/conjunction/ConjunctionAssessment.ts new file mode 100644 index 00000000..f3627fcc --- /dev/null +++ b/src/conjunction/ConjunctionAssessment.ts @@ -0,0 +1,417 @@ +/** + * @author Theodore Kruczek + * @license AGPL-3.0-or-later + * @copyright (c) 2025 Kruczek Labs LLC + * + * Orbital Object ToolKit is free software: you can redistribute it and/or modify it under the + * terms of the GNU Affero General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later version. + * + * Orbital Object ToolKit is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; + * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License along with + * Orbital Object ToolKit. If not, see . + */ + +import { EpochUTC } from '../time/EpochUTC.js'; +import { J2000 } from '../coordinate/J2000.js'; +import { RIC } from '../coordinate/RIC.js'; +import { Tle } from '../coordinate/Tle.js'; +import { StateCovariance, CovarianceFrame } from '../covariance/StateCovariance.js'; +import { CovarianceSample } from '../covariance/CovarianceSample.js'; +import { Propagator } from '../propagator/Propagator.js'; +import { Sgp4Propagator } from '../propagator/Sgp4Propagator.js'; +import { RungeKutta89Propagator } from '../propagator/RungeKutta89Propagator.js'; +import { ForceModel } from '../force/ForceModel.js'; +import { GoldenSection } from '../optimize/GoldenSection.js'; +import { Matrix } from '../operations/Matrix.js'; +import type { Kilometers, Seconds } from '../main.js'; +import { ConjunctionEvent } from './ConjunctionEvent.js'; +import { ProbabilityOfCollision } from './ProbabilityOfCollision.js'; + +/** + * Input for a space object in conjunction assessment. + * Can be specified as either a TLE or a state vector with optional covariance. + */ +export interface SpaceObject { + /** Object identifier/name */ + name?: string; + + /** TLE for the object (alternative to state) */ + tle?: Tle; + + /** State vector in J2000 frame (alternative to TLE) */ + state?: J2000; + + /** Covariance matrix (optional, for probability calculation) */ + covariance?: StateCovariance; + + /** Hard body radius in km (optional, for probability calculation) */ + radius?: Kilometers; + + /** Custom propagator (optional, overrides default) */ + propagator?: Propagator; +} + +/** + * Configuration options for conjunction assessment. + */ +export interface ConjunctionAssessmentOptions { + /** Search window start time */ + startTime: EpochUTC; + + /** Search window end time */ + endTime: EpochUTC; + + /** Use high-fidelity propagation (RungeKutta89 instead of SGP4) */ + useHighFidelity?: boolean; + + /** Force model for high-fidelity propagation (optional) */ + forceModel?: ForceModel; + + /** Propagate covariances using sigma-point method */ + propagateCovariance?: boolean; + + /** TCA search tolerance in seconds */ + tcaTolerance?: Seconds; + + /** Step size for initial TCA search in seconds */ + searchStepSize?: Seconds; +} + +/** + * High accuracy conjunction assessment workflow. + * + * This class provides a comprehensive workflow for assessing conjunctions between + * space objects using: + * - High accuracy propagators (SGP4 or numerical integrators) + * - Covariance matrices based on historical TLE accuracy + * - Time of Closest Approach (TCA) finding using optimization + * - Probability of collision calculation using Chan's 2D method + * + * @example + * ```typescript + * const primaryTle = new Tle(line1, line2); + * const secondaryTle = new Tle(line1, line2); + * + * const assessment = new ConjunctionAssessment( + * { tle: primaryTle, radius: 0.01 as Kilometers }, + * { tle: secondaryTle, radius: 0.01 as Kilometers }, + * ); + * + * const event = assessment.assess({ + * startTime: EpochUTC.fromDateTime(new Date('2025-01-01T00:00:00Z')), + * endTime: EpochUTC.fromDateTime(new Date('2025-01-02T00:00:00Z')), + * useHighFidelity: true, + * propagateCovariance: true, + * }); + * + * console.log(event.toString()); + * ``` + */ +export class ConjunctionAssessment { + private primaryProp: Propagator; + private secondaryProp: Propagator; + private primaryCovSample?: CovarianceSample; + private secondaryCovSample?: CovarianceSample; + + constructor( + private primary: SpaceObject, + private secondary: SpaceObject, + ) { + // Initialize propagators (will be replaced in assess() if needed) + this.primaryProp = this.createPropagator(primary, false); + this.secondaryProp = this.createPropagator(secondary, false); + } + + /** + * Performs conjunction assessment over the specified time window. + * + * @param options Assessment configuration options + * @returns ConjunctionEvent with TCA, miss distance, and probability of collision + */ + assess(options: ConjunctionAssessmentOptions): ConjunctionEvent { + const { + startTime, + endTime, + useHighFidelity = false, + forceModel, + propagateCovariance = false, + tcaTolerance = 0.001 as Seconds, + searchStepSize = 60.0 as Seconds, + } = options; + + // Create propagators based on options + this.primaryProp = this.primary.propagator ?? this.createPropagator(this.primary, useHighFidelity, forceModel); + this.secondaryProp = + this.secondary.propagator ?? this.createPropagator(this.secondary, useHighFidelity, forceModel); + + // Initialize covariance samples if requested + if (propagateCovariance) { + this.initializeCovarianceSamples(useHighFidelity, forceModel); + } + + // Find Time of Closest Approach (TCA) + const tca = this.findTCA(startTime, endTime, searchStepSize, tcaTolerance); + + // Propagate states to TCA + const primaryState = this.primaryProp.propagate(tca); + const secondaryState = this.secondaryProp.propagate(tca); + + // Compute relative state in RIC frame + const relativeState = RIC.fromJ2000(secondaryState, primaryState); + + // Extract RIC components + const radialDistance = Math.abs(relativeState.position.x) as Kilometers; + const intrackDistance = Math.abs(relativeState.position.y) as Kilometers; + const crosstrackDistance = Math.abs(relativeState.position.z) as Kilometers; + const missDistance = relativeState.range; + const relativeVelocity = relativeState.velocity.magnitude(); + + // Compute combined covariance and Pc if available + let combinedCovariance: StateCovariance | undefined; + let probabilityOfCollision: number | undefined; + + if (this.primaryCovSample && this.secondaryCovSample && this.primary.radius && this.secondary.radius) { + // Propagate covariance samples to TCA + this.primaryCovSample.propagate(tca); + this.secondaryCovSample.propagate(tca); + + // Get covariances in RIC frame + const primaryCov = this.primaryCovSample.covarianceRIC(); + const secondaryCov = this.secondaryCovSample.covarianceRIC(); + + // Combine covariances + combinedCovariance = ProbabilityOfCollision.combineCovarianceMatrices(primaryCov, secondaryCov); + + // Calculate probability of collision + const combinedRadius = (this.primary.radius + this.secondary.radius) as Kilometers; + + probabilityOfCollision = ProbabilityOfCollision.calculate( + relativeState.position, + relativeState.velocity, + combinedCovariance, + combinedRadius, + ); + } else if (this.primary.covariance && this.secondary.covariance && this.primary.radius && this.secondary.radius) { + // Use provided covariances without propagation + // Transform to RIC if needed + let primaryCovRIC = this.primary.covariance; + let secondaryCovRIC = this.secondary.covariance; + + if (this.primary.covariance.frame === CovarianceFrame.ECI) { + primaryCovRIC = this.transformCovarianceToRIC(this.primary.covariance, primaryState); + } + if (this.secondary.covariance.frame === CovarianceFrame.ECI) { + secondaryCovRIC = this.transformCovarianceToRIC(this.secondary.covariance, secondaryState); + } + + combinedCovariance = ProbabilityOfCollision.combineCovarianceMatrices(primaryCovRIC, secondaryCovRIC); + + const combinedRadius = (this.primary.radius + this.secondary.radius) as Kilometers; + + probabilityOfCollision = ProbabilityOfCollision.calculate( + relativeState.position, + relativeState.velocity, + combinedCovariance, + combinedRadius, + ); + } + + return new ConjunctionEvent( + tca, + primaryState, + secondaryState, + relativeState, + missDistance, + radialDistance, + intrackDistance, + crosstrackDistance, + relativeVelocity, + combinedCovariance, + probabilityOfCollision, + this.primary.radius, + this.secondary.radius, + ); + } + + /** + * Finds the Time of Closest Approach (TCA) using golden section search. + * + * @param startTime Search window start + * @param endTime Search window end + * @param stepSize Initial search step size in seconds + * @param tolerance TCA search tolerance in seconds + * @returns TCA epoch + */ + private findTCA(startTime: EpochUTC, endTime: EpochUTC, stepSize: Seconds, tolerance: Seconds): EpochUTC { + // Coarse search to find approximate TCA region + let minRange = Infinity; + let minEpoch = startTime; + let current = startTime; + + while (current.posix <= endTime.posix) { + const primary = this.primaryProp.propagate(current); + const secondary = this.secondaryProp.propagate(current); + const ric = RIC.fromJ2000(secondary, primary); + const range = ric.range; + + if (range < minRange) { + minRange = range; + minEpoch = current; + } + + current = current.roll(stepSize); + } + + // Fine search using golden section optimization + const searchWindow = 2 * stepSize; + const lowerBound = Math.max(startTime.posix, minEpoch.posix - searchWindow); + const upperBound = Math.min(endTime.posix, minEpoch.posix + searchWindow); + + const tcaPosix = GoldenSection.search( + (posix) => { + const epoch = new EpochUTC(posix as Seconds); + const primary = this.primaryProp.propagate(epoch); + const secondary = this.secondaryProp.propagate(epoch); + const ric = RIC.fromJ2000(secondary, primary); + + return ric.range; + }, + lowerBound, + upperBound, + { tolerance }, + ); + + return new EpochUTC(tcaPosix as Seconds); + } + + /** + * Creates a propagator for a space object. + * + * @param obj Space object + * @param useHighFidelity Whether to use high-fidelity propagation + * @param forceModel Optional force model for numerical propagation + * @returns Propagator instance + */ + private createPropagator(obj: SpaceObject, useHighFidelity: boolean, forceModel?: ForceModel): Propagator { + if (obj.tle) { + if (useHighFidelity) { + // Convert TLE to state and use RK89 + const state = obj.tle.toJ2000(); + const fm = forceModel ?? new ForceModel().setGravity(); + + return new RungeKutta89Propagator(state, fm); + } + + return new Sgp4Propagator(obj.tle); + } else if (obj.state) { + if (useHighFidelity) { + const fm = forceModel ?? new ForceModel().setGravity(); + + return new RungeKutta89Propagator(obj.state, fm); + } + + // For low-fidelity, still use RK89 but with simpler force model + const fm = new ForceModel().setGravity(); + + return new RungeKutta89Propagator(obj.state, fm); + } + + throw new Error('Space object must have either a TLE or state vector'); + } + + /** + * Initializes covariance samples for both objects. + * + * @param useHighFidelity Whether to use high-fidelity propagation + * @param forceModel Optional force model + */ + private initializeCovarianceSamples(useHighFidelity: boolean, forceModel?: ForceModel): void { + const fm = forceModel ?? new ForceModel().setGravity(); + + // Primary covariance + if (this.primary.tle) { + const state = this.primary.tle.toJ2000(); + const covariance = + this.primary.covariance ?? + StateCovariance.fromSigmas([1.0, 1.0, 1.0, 0.001, 0.001, 0.001], CovarianceFrame.RIC); + + this.primaryCovSample = new CovarianceSample(state, covariance, this.primary.tle, fm, fm); + } else if (this.primary.state && this.primary.covariance) { + this.primaryCovSample = new CovarianceSample(this.primary.state, this.primary.covariance, undefined, fm, fm); + } + + // Secondary covariance + if (this.secondary.tle) { + const state = this.secondary.tle.toJ2000(); + const covariance = + this.secondary.covariance ?? + StateCovariance.fromSigmas([1.0, 1.0, 1.0, 0.001, 0.001, 0.001], CovarianceFrame.RIC); + + this.secondaryCovSample = new CovarianceSample(state, covariance, this.secondary.tle, fm, fm); + } else if (this.secondary.state && this.secondary.covariance) { + this.secondaryCovSample = new CovarianceSample(this.secondary.state, this.secondary.covariance, undefined, fm, fm); + } + } + + /** + * Transforms ECI covariance to RIC frame. + * + * @param covariance ECI covariance + * @param state State vector for RIC frame definition + * @returns RIC covariance + */ + private transformCovarianceToRIC(covariance: StateCovariance, state: J2000): StateCovariance { + // Create RIC transformation matrix (3x3) + const ricMatrix = this.createRICTransformMatrix(state.position, state.velocity); + + // Build 6x6 transformation matrix (block diagonal with 3x3 rotation) + const transform = this.build6x6Transform(ricMatrix); + + // Transform covariance: C_ric = T * C_eci * T^T + const covRIC = transform.multiply(covariance.matrix).multiply(transform.transpose()); + + return new StateCovariance(covRIC, CovarianceFrame.RIC); + } + + /** + * Creates the 3x3 RIC transformation matrix. + * + * @param position Position vector + * @param velocity Velocity vector + * @returns 3x3 RIC transformation matrix + */ + private createRICTransformMatrix(position: any, velocity: any): Matrix { + const ru = position.normalize(); + const cu = position.cross(velocity).normalize(); + const iu = cu.cross(ru).normalize(); + + return new Matrix([ + [ru.x, ru.y, ru.z], + [iu.x, iu.y, iu.z], + [cu.x, cu.y, cu.z], + ]); + } + + /** + * Builds a 6x6 transformation matrix from a 3x3 rotation matrix. + * + * @param rot 3x3 rotation matrix + * @returns 6x6 transformation matrix + */ + private build6x6Transform(rot: Matrix): Matrix { + const elements = [ + [rot.elements[0][0], rot.elements[0][1], rot.elements[0][2], 0, 0, 0], + [rot.elements[1][0], rot.elements[1][1], rot.elements[1][2], 0, 0, 0], + [rot.elements[2][0], rot.elements[2][1], rot.elements[2][2], 0, 0, 0], + [0, 0, 0, rot.elements[0][0], rot.elements[0][1], rot.elements[0][2]], + [0, 0, 0, rot.elements[1][0], rot.elements[1][1], rot.elements[1][2]], + [0, 0, 0, rot.elements[2][0], rot.elements[2][1], rot.elements[2][2]], + ]; + + return new Matrix(elements); + } +} diff --git a/src/conjunction/ConjunctionEvent.ts b/src/conjunction/ConjunctionEvent.ts new file mode 100644 index 00000000..26dd7a1d --- /dev/null +++ b/src/conjunction/ConjunctionEvent.ts @@ -0,0 +1,149 @@ +/** + * @author Theodore Kruczek + * @license AGPL-3.0-or-later + * @copyright (c) 2025 Kruczek Labs LLC + * + * Orbital Object ToolKit is free software: you can redistribute it and/or modify it under the + * terms of the GNU Affero General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later version. + * + * Orbital Object ToolKit is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; + * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License along with + * Orbital Object ToolKit. If not, see . + */ + +import type { EpochUTC } from '../time/EpochUTC.js'; +import type { J2000 } from '../coordinate/J2000.js'; +import type { RIC } from '../coordinate/RIC.js'; +import type { StateCovariance } from '../covariance/StateCovariance.js'; +import type { Kilometers, KilometersPerSecond } from '../main.js'; +import { Matrix } from '../operations/Matrix.js'; + +/** + * Represents the result of a conjunction assessment between two space objects. + * Contains all relevant information about the close approach event. + */ +export class ConjunctionEvent { + constructor( + /** Time of Closest Approach (TCA) */ + public tca: EpochUTC, + /** Primary object state at TCA in J2000 frame */ + public primaryState: J2000, + /** Secondary object state at TCA in J2000 frame */ + public secondaryState: J2000, + /** Relative state in RIC frame (relative to primary) */ + public relativeState: RIC, + /** Total miss distance at TCA (km) */ + public missDistance: Kilometers, + /** Radial component of miss distance (km) */ + public radialDistance: Kilometers, + /** Intrack component of miss distance (km) */ + public intrackDistance: Kilometers, + /** Crosstrack component of miss distance (km) */ + public crosstrackDistance: Kilometers, + /** Relative velocity magnitude at TCA (km/s) */ + public relativeVelocity: KilometersPerSecond, + /** Combined position covariance matrix in RIC frame (optional) */ + public combinedCovariance?: StateCovariance, + /** Probability of collision (optional, 0-1) */ + public probabilityOfCollision?: number, + /** Hard body radius for primary object (km, optional) */ + public primaryRadius?: Kilometers, + /** Hard body radius for secondary object (km, optional) */ + public secondaryRadius?: Kilometers, + ) {} + + /** + * Returns a formatted string representation of the conjunction event. + * @returns A multi-line string with conjunction details. + */ + toString(): string { + const lines = [ + '[Conjunction Event]', + ` TCA: ${this.tca.toString()}`, + ` Miss Distance: ${this.missDistance.toFixed(6)} km`, + ` Radial: ${this.radialDistance.toFixed(6)} km`, + ` Intrack: ${this.intrackDistance.toFixed(6)} km`, + ` Crosstrack: ${this.crosstrackDistance.toFixed(6)} km`, + ` Relative Velocity: ${this.relativeVelocity.toFixed(6)} km/s`, + ]; + + if (this.probabilityOfCollision !== undefined) { + lines.push(` Probability of Collision: ${this.probabilityOfCollision.toExponential(6)}`); + } + + if (this.primaryRadius !== undefined && this.secondaryRadius !== undefined) { + const combinedRadius = this.primaryRadius + this.secondaryRadius; + + lines.push(` Combined Hard Body Radius: ${combinedRadius.toFixed(3)} km`); + } + + return lines.join('\n'); + } + + /** + * Checks if this is a high-risk conjunction based on miss distance and Pc. + * @param distanceThreshold Miss distance threshold in km (default: 1.0 km) + * @param pcThreshold Probability of collision threshold (default: 1e-4) + * @returns True if the conjunction exceeds risk thresholds. + */ + isHighRisk(distanceThreshold: Kilometers = 1.0 as Kilometers, pcThreshold: number = 1e-4): boolean { + const distanceRisk = this.missDistance < distanceThreshold; + const pcRisk = this.probabilityOfCollision !== undefined && this.probabilityOfCollision > pcThreshold; + + return distanceRisk || pcRisk; + } + + /** + * Gets the Mahalanobis distance if covariance is available. + * This is the miss distance normalized by the combined covariance. + * @returns The Mahalanobis distance (unitless), or undefined if no covariance. + */ + getMahalanobisDistance(): number | undefined { + if (!this.combinedCovariance) { + return undefined; + } + + // Extract position-only covariance (first 3x3 block) + const posCovariance = this.extractPositionCovariance(this.combinedCovariance.matrix); + const relativePosition = [ + this.relativeState.position.x, + this.relativeState.position.y, + this.relativeState.position.z, + ]; + + try { + // Compute Mahalanobis distance: sqrt(r^T * C^-1 * r) + const covInv = posCovariance.inverse(); + const temp = covInv.multiplyVector(relativePosition); + let mahalanobis = 0; + + for (let i = 0; i < 3; i++) { + mahalanobis += relativePosition[i] * temp[i]; + } + + return Math.sqrt(mahalanobis); + } catch { + // Covariance may be singular + return undefined; + } + } + + /** + * Extracts the 3x3 position covariance from a 6x6 state covariance matrix. + * @param stateCov 6x6 state covariance matrix + * @returns 3x3 position covariance matrix + */ + private extractPositionCovariance(stateCov: Matrix): Matrix { + const elements = [ + [stateCov.elements[0][0], stateCov.elements[0][1], stateCov.elements[0][2]], + [stateCov.elements[1][0], stateCov.elements[1][1], stateCov.elements[1][2]], + [stateCov.elements[2][0], stateCov.elements[2][1], stateCov.elements[2][2]], + ]; + + return new Matrix(elements); + } +} diff --git a/src/conjunction/ProbabilityOfCollision.ts b/src/conjunction/ProbabilityOfCollision.ts new file mode 100644 index 00000000..5291155c --- /dev/null +++ b/src/conjunction/ProbabilityOfCollision.ts @@ -0,0 +1,300 @@ +/** + * @author Theodore Kruczek + * @license AGPL-3.0-or-later + * @copyright (c) 2025 Kruczek Labs LLC + * + * Orbital Object ToolKit is free software: you can redistribute it and/or modify it under the + * terms of the GNU Affero General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later version. + * + * Orbital Object ToolKit is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; + * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License along with + * Orbital Object ToolKit. If not, see . + */ + +import { Matrix } from '../operations/Matrix.js'; +import { Vector3D } from '../operations/Vector3D.js'; +import type { StateCovariance } from '../covariance/StateCovariance.js'; +import type { Kilometers } from '../main.js'; + +/** + * Probability of Collision calculator using Chan's 2D method. + * + * This implementation projects the combined covariance matrix onto the + * encounter plane (B-plane) perpendicular to the relative velocity vector, + * then computes the probability that the relative position lies within + * the combined hard body radius. + * + * Reference: Chan, F. K. (2008). "Spacecraft Collision Probability" + */ +export class ProbabilityOfCollision { + /** + * Computes probability of collision using Chan's 2D method. + * + * @param relativePosition Relative position vector in RIC frame (km) + * @param relativeVelocity Relative velocity vector in RIC frame (km/s) + * @param combinedCovariance Combined position covariance in RIC frame (6x6) + * @param combinedRadius Combined hard body radius (km) + * @returns Probability of collision (0 to 1) + */ + static calculate( + relativePosition: Vector3D, + relativeVelocity: Vector3D, + combinedCovariance: StateCovariance, + combinedRadius: Kilometers, + ): number { + // Extract position covariance (first 3x3 block) + const posCovariance = this.extractPositionCovariance(combinedCovariance.matrix); + + // Create encounter plane coordinate system + // z-axis: along relative velocity (perpendicular to encounter plane) + // x-axis and y-axis: in the encounter plane + const vMag = relativeVelocity.magnitude(); + + if (vMag < 1e-9) { + // Nearly zero relative velocity - use 3D method or return conservative estimate + return this.calculate3D(relativePosition, posCovariance, combinedRadius); + } + + const zAxis = relativeVelocity.scale(1.0 / vMag); + + // Choose x-axis perpendicular to z-axis + let xAxis: Vector3D; + + if (Math.abs(zAxis.z) < 0.9) { + // z-axis is not close to [0,0,1], use cross product with [0,0,1] + xAxis = new Vector3D(0, 0, 1).cross(zAxis).normalize(); + } else { + // z-axis is close to [0,0,1], use cross product with [1,0,0] + xAxis = new Vector3D(1, 0, 0).cross(zAxis).normalize(); + } + + const yAxis = zAxis.cross(xAxis).normalize(); + + // Rotation matrix from RIC to encounter plane frame + const rotMatrix = new Matrix([ + [xAxis.x, xAxis.y, xAxis.z], + [yAxis.x, yAxis.y, yAxis.z], + [zAxis.x, zAxis.y, zAxis.z], + ]); + + // Transform position to encounter plane frame + const posEncounter = rotMatrix.multiplyVector3D(relativePosition); + + // Project position onto encounter plane (drop z-component) + const x = posEncounter.x; + const y = posEncounter.y; + + // Transform covariance to encounter plane frame: C_enc = R * C * R^T + const covEncounter = rotMatrix.multiply(posCovariance).multiply(rotMatrix.transpose()); + + // Extract 2D covariance in encounter plane (top-left 2x2 block) + const cov2D = new Matrix([ + [covEncounter.elements[0][0], covEncounter.elements[0][1]], + [covEncounter.elements[1][0], covEncounter.elements[1][1]], + ]); + + // Compute Pc using 2D method + return this.calculatePc2D(x, y, cov2D, combinedRadius); + } + + /** + * Calculates 2D probability of collision in the encounter plane. + * + * Uses the analytical solution for 2D Gaussian probability within a circle. + * + * @param x X-coordinate in encounter plane (km) + * @param y Y-coordinate in encounter plane (km) + * @param cov2D 2x2 covariance matrix in encounter plane + * @param radius Combined hard body radius (km) + * @returns Probability of collision (0 to 1) + */ + private static calculatePc2D(x: number, y: number, cov2D: Matrix, radius: number): number { + // Compute determinant and trace + const c11 = cov2D.elements[0][0]; + const c12 = cov2D.elements[0][1]; + const c22 = cov2D.elements[1][1]; + + const det = c11 * c22 - c12 * c12; + + if (det <= 0) { + // Singular or invalid covariance + return 0; + } + + // Compute Mahalanobis distance squared: d^2 = r^T * C^-1 * r + const covInv11 = c22 / det; + const covInv12 = -c12 / det; + const covInv22 = c11 / det; + + const d2 = x * (covInv11 * x + covInv12 * y) + y * (covInv12 * x + covInv22 * y); + + // Compute eigenvalues for ellipse semi-axes + const trace = c11 + c22; + const discriminant = Math.sqrt((c11 - c22) ** 2 + 4 * c12 * c12); + const lambda1 = (trace + discriminant) / 2; + const lambda2 = (trace - discriminant) / 2; + + const sigma1 = Math.sqrt(Math.max(lambda1, 0)); + const sigma2 = Math.sqrt(Math.max(lambda2, 0)); + + // Use Chan's analytical approximation + return this.chanPc2D(Math.sqrt(d2), radius, sigma1, sigma2); + } + + /** + * Chan's analytical approximation for 2D Pc. + * + * Reference: Chan, F. K. (2008). "Spacecraft Collision Probability" + * + * @param mahalanobisDistance Mahalanobis distance (sqrt of d2) + * @param radius Combined hard body radius + * @param sigma1 Larger semi-axis of covariance ellipse + * @param sigma2 Smaller semi-axis of covariance ellipse + * @returns Probability of collision (0 to 1) + */ + private static chanPc2D( + mahalanobisDistance: number, + radius: number, + sigma1: number, + sigma2: number, + ): number { + // Effective radius in Mahalanobis space + const sigmaMean = Math.sqrt(sigma1 * sigma2); + + if (sigmaMean < 1e-12) { + return 0; + } + + const eta = radius / sigmaMean; + + // If miss distance is much larger than combined size, Pc is negligible + if (mahalanobisDistance > eta + 10) { + return 0; + } + + // If objects are overlapping at nominal position + if (mahalanobisDistance * sigmaMean < radius) { + // Use complementary error function approximation + const zeta = mahalanobisDistance; + + return Math.exp(-0.5 * zeta * zeta) * (1 - this.approximateErfc(eta / Math.sqrt(2))) / 2; + } + + // General case: use Chan's approximation + const u = mahalanobisDistance; + const gamma = sigma1 / sigma2; + + // Foster's approximation (simplified Chan method) + const pc = (eta * eta) / (2 * (u * u + eta * eta)) * Math.exp(-0.5 * u * u); + + return Math.min(Math.max(pc, 0), 1); // Clamp to [0, 1] + } + + /** + * Approximates the complementary error function erfc(x). + * + * Uses Abramowitz and Stegun approximation (formula 7.1.26). + * + * @param x Input value + * @returns erfc(x) + */ + private static approximateErfc(x: number): number { + if (x < 0) { + return 2 - this.approximateErfc(-x); + } + + // Constants for Abramowitz and Stegun approximation + const p = 0.3275911; + const a1 = 0.254829592; + const a2 = -0.284496736; + const a3 = 1.421413741; + const a4 = -1.453152027; + const a5 = 1.061405429; + + const t = 1.0 / (1.0 + p * x); + const t2 = t * t; + const t3 = t2 * t; + const t4 = t3 * t; + const t5 = t4 * t; + + return (a1 * t + a2 * t2 + a3 * t3 + a4 * t4 + a5 * t5) * Math.exp(-x * x); + } + + /** + * Fallback 3D probability calculation for cases with very low relative velocity. + * + * Uses simple spherical approximation. + * + * @param relativePosition Relative position vector + * @param posCovariance 3x3 position covariance matrix + * @param combinedRadius Combined hard body radius + * @returns Probability of collision (0 to 1) + */ + private static calculate3D( + relativePosition: Vector3D, + posCovariance: Matrix, + combinedRadius: number, + ): number { + const r = relativePosition.magnitude(); + + // Compute average variance (trace / 3) + const avgVariance = + (posCovariance.elements[0][0] + posCovariance.elements[1][1] + posCovariance.elements[2][2]) / 3; + const sigma = Math.sqrt(avgVariance); + + if (sigma < 1e-12) { + return r < combinedRadius ? 1.0 : 0.0; + } + + // Use Gaussian CDF approximation + const z = (r - combinedRadius) / sigma; + + if (z < -3) { + return 1.0; + } + if (z > 3) { + return 0.0; + } + + // Approximate using complementary error function + return 0.5 * this.approximateErfc(z / Math.sqrt(2)); + } + + /** + * Combines two covariance matrices (primary and secondary). + * + * The combined covariance is simply the sum of the two covariances, + * assuming they are independent. + * + * @param primaryCov Primary object covariance + * @param secondaryCov Secondary object covariance + * @returns Combined covariance + */ + static combineCovarianceMatrices(primaryCov: StateCovariance, secondaryCov: StateCovariance): StateCovariance { + const combinedMatrix = primaryCov.matrix.add(secondaryCov.matrix); + + return { + matrix: combinedMatrix, + frame: primaryCov.frame, + } as StateCovariance; + } + + /** + * Extracts the 3x3 position covariance from a 6x6 state covariance matrix. + * @param stateCov 6x6 state covariance matrix + * @returns 3x3 position covariance matrix + */ + private static extractPositionCovariance(stateCov: Matrix): Matrix { + const elements = [ + [stateCov.elements[0][0], stateCov.elements[0][1], stateCov.elements[0][2]], + [stateCov.elements[1][0], stateCov.elements[1][1], stateCov.elements[1][2]], + [stateCov.elements[2][0], stateCov.elements[2][1], stateCov.elements[2][2]], + ]; + + return new Matrix(elements); + } +} diff --git a/src/conjunction/index.ts b/src/conjunction/index.ts new file mode 100644 index 00000000..329fc441 --- /dev/null +++ b/src/conjunction/index.ts @@ -0,0 +1,24 @@ +/** + * @author Theodore Kruczek + * @license AGPL-3.0-or-later + * @copyright (c) 2025 Kruczek Labs LLC + * + * Orbital Object ToolKit is free software: you can redistribute it and/or modify it under the + * terms of the GNU Affero General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later version. + * + * Orbital Object ToolKit is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; + * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License along with + * Orbital Object ToolKit. If not, see . + */ + +export { ConjunctionEvent } from './ConjunctionEvent.js'; +export { ProbabilityOfCollision } from './ProbabilityOfCollision.js'; +export { + ConjunctionAssessment, + type SpaceObject, + type ConjunctionAssessmentOptions, +} from './ConjunctionAssessment.js'; diff --git a/src/main.ts b/src/main.ts index 39893076..992b6665 100644 --- a/src/main.ts +++ b/src/main.ts @@ -59,4 +59,6 @@ export * from './propagator/index'; export * from './orbit_determination/index'; -export * from './covariance/index'; +export * from './covariance/index.js'; + +export * from './conjunction/index.js'; diff --git a/test/conjunction/ConjunctionAssessment.test.ts b/test/conjunction/ConjunctionAssessment.test.ts new file mode 100644 index 00000000..6df9158a --- /dev/null +++ b/test/conjunction/ConjunctionAssessment.test.ts @@ -0,0 +1,296 @@ +/** + * @author Theodore Kruczek + * @license AGPL-3.0-or-later + * @copyright (c) 2025 Kruczek Labs LLC + */ + +import { + ConjunctionAssessment, + ConjunctionEvent, + EpochUTC, + J2000, + Kilometers, + KilometersPerSecond, + ProbabilityOfCollision, + StateCovariance, + CovarianceFrame, + Tle, + Vector3D, +} from '../../src/main'; + +describe('ConjunctionAssessment', () => { + // Sample TLE data for testing (ISS and a close approach satellite) + const tleLine1Primary = '1 25544U 98067A 25019.50000000 .00016717 00000-0 10270-3 0 9005'; + const tleLine2Primary = '2 25544 51.6400 339.8000 0002571 90.5000 269.6000 15.50000000000000'; + const tleLine1Secondary = '1 25544U 98067A 25019.50000000 .00016717 00000-0 10270-3 0 9006'; + const tleLine2Secondary = '2 25544 51.6400 339.8100 0002571 90.5100 269.6100 15.50000000000000'; + + let primaryTle: Tle; + let secondaryTle: Tle; + + beforeEach(() => { + primaryTle = new Tle(tleLine1Primary, tleLine2Primary); + secondaryTle = new Tle(tleLine1Secondary, tleLine2Secondary); + }); + + describe('Basic Conjunction Assessment', () => { + it('should create a conjunction assessment from TLEs', () => { + const assessment = new ConjunctionAssessment( + { tle: primaryTle, radius: 0.01 as Kilometers }, + { tle: secondaryTle, radius: 0.01 as Kilometers }, + ); + + expect(assessment).toBeDefined(); + }); + + it('should find TCA and compute miss distance', () => { + const assessment = new ConjunctionAssessment( + { tle: primaryTle, radius: 0.01 as Kilometers }, + { tle: secondaryTle, radius: 0.01 as Kilometers }, + ); + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T18:00:00.000Z'); + + const event = assessment.assess({ + startTime, + endTime, + }); + + expect(event).toBeDefined(); + expect(event.tca).toBeDefined(); + expect(event.missDistance).toBeGreaterThan(0); + expect(event.relativeVelocity).toBeGreaterThan(0); + }); + + it('should provide RIC components of miss distance', () => { + const assessment = new ConjunctionAssessment( + { tle: primaryTle, radius: 0.01 as Kilometers }, + { tle: secondaryTle, radius: 0.01 as Kilometers }, + ); + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T18:00:00.000Z'); + + const event = assessment.assess({ + startTime, + endTime, + }); + + expect(event.radialDistance).toBeGreaterThanOrEqual(0); + expect(event.intrackDistance).toBeGreaterThanOrEqual(0); + expect(event.crosstrackDistance).toBeGreaterThanOrEqual(0); + + // Verify that RIC components sum to total miss distance + const computedMiss = Math.sqrt( + event.radialDistance ** 2 + event.intrackDistance ** 2 + event.crosstrackDistance ** 2, + ); + + expect(computedMiss).toBeCloseTo(event.missDistance, 6); + }); + }); + + describe('High-Fidelity Propagation', () => { + it('should use high-fidelity propagation when requested', () => { + const assessment = new ConjunctionAssessment( + { tle: primaryTle, radius: 0.01 as Kilometers }, + { tle: secondaryTle, radius: 0.01 as Kilometers }, + ); + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T13:00:00.000Z'); + + const event = assessment.assess({ + startTime, + endTime, + useHighFidelity: true, + }); + + expect(event).toBeDefined(); + expect(event.tca).toBeDefined(); + }); + }); + + describe('Covariance Propagation', () => { + it('should propagate covariances and compute Pc', () => { + const assessment = new ConjunctionAssessment( + { tle: primaryTle, radius: 0.01 as Kilometers }, + { tle: secondaryTle, radius: 0.01 as Kilometers }, + ); + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T13:00:00.000Z'); + + const event = assessment.assess({ + startTime, + endTime, + useHighFidelity: true, + propagateCovariance: true, + }); + + expect(event).toBeDefined(); + expect(event.combinedCovariance).toBeDefined(); + expect(event.probabilityOfCollision).toBeDefined(); + expect(event.probabilityOfCollision).toBeGreaterThanOrEqual(0); + expect(event.probabilityOfCollision).toBeLessThanOrEqual(1); + }); + }); + + describe('State Vector Input', () => { + it('should accept state vectors instead of TLEs', () => { + const primaryState = new J2000( + EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'), + new Vector3D(6878.0 as Kilometers, 0 as Kilometers, 0 as Kilometers), + new Vector3D(0 as KilometersPerSecond, 7.5 as KilometersPerSecond, 0 as KilometersPerSecond), + ); + + const secondaryState = new J2000( + EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'), + new Vector3D(6878.1 as Kilometers, 0 as Kilometers, 0 as Kilometers), + new Vector3D(0 as KilometersPerSecond, 7.5 as KilometersPerSecond, 0 as KilometersPerSecond), + ); + + const assessment = new ConjunctionAssessment( + { state: primaryState, radius: 0.01 as Kilometers }, + { state: secondaryState, radius: 0.01 as Kilometers }, + ); + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T13:00:00.000Z'); + + const event = assessment.assess({ + startTime, + endTime, + }); + + expect(event).toBeDefined(); + expect(event.missDistance).toBeGreaterThan(0); + }); + + it('should use provided covariances for Pc calculation', () => { + const primaryState = new J2000( + EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'), + new Vector3D(6878.0 as Kilometers, 0 as Kilometers, 0 as Kilometers), + new Vector3D(0 as KilometersPerSecond, 7.5 as KilometersPerSecond, 0 as KilometersPerSecond), + ); + + const secondaryState = new J2000( + EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'), + new Vector3D(6878.1 as Kilometers, 0 as Kilometers, 0 as Kilometers), + new Vector3D(0 as KilometersPerSecond, 7.5 as KilometersPerSecond, 0 as KilometersPerSecond), + ); + + const covariance = StateCovariance.fromSigmas([1.0, 1.0, 1.0, 0.001, 0.001, 0.001], CovarianceFrame.RIC); + + const assessment = new ConjunctionAssessment( + { state: primaryState, covariance, radius: 0.01 as Kilometers }, + { state: secondaryState, covariance, radius: 0.01 as Kilometers }, + ); + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T13:00:00.000Z'); + + const event = assessment.assess({ + startTime, + endTime, + }); + + expect(event.probabilityOfCollision).toBeDefined(); + expect(event.probabilityOfCollision).toBeGreaterThanOrEqual(0); + expect(event.probabilityOfCollision).toBeLessThanOrEqual(1); + }); + }); + + describe('ConjunctionEvent', () => { + it('should format conjunction event as string', () => { + const assessment = new ConjunctionAssessment( + { tle: primaryTle, radius: 0.01 as Kilometers }, + { tle: secondaryTle, radius: 0.01 as Kilometers }, + ); + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T13:00:00.000Z'); + + const event = assessment.assess({ + startTime, + endTime, + }); + + const str = event.toString(); + + expect(str).toContain('Conjunction Event'); + expect(str).toContain('TCA:'); + expect(str).toContain('Miss Distance:'); + }); + + it('should identify high-risk conjunctions', () => { + const primaryState = new J2000( + EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'), + new Vector3D(6878.0 as Kilometers, 0 as Kilometers, 0 as Kilometers), + new Vector3D(0 as KilometersPerSecond, 7.5 as KilometersPerSecond, 0 as KilometersPerSecond), + ); + + const secondaryState = new J2000( + EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'), + new Vector3D(6878.0005 as Kilometers, 0 as Kilometers, 0 as Kilometers), + new Vector3D(0 as KilometersPerSecond, 7.5 as KilometersPerSecond, 0 as KilometersPerSecond), + ); + + const assessment = new ConjunctionAssessment( + { state: primaryState, radius: 0.01 as Kilometers }, + { state: secondaryState, radius: 0.01 as Kilometers }, + ); + + const startTime = EpochUTC.fromDateTimeString('2025-01-19T12:00:00.000Z'); + const endTime = EpochUTC.fromDateTimeString('2025-01-19T13:00:00.000Z'); + + const event = assessment.assess({ + startTime, + endTime, + }); + + expect(event.isHighRisk(1.0 as Kilometers)).toBe(true); + }); + }); + + describe('ProbabilityOfCollision', () => { + it('should compute Pc for a given relative state and covariance', () => { + const relativePosition = new Vector3D(0.5 as Kilometers, 0.5 as Kilometers, 0.5 as Kilometers); + const relativeVelocity = new Vector3D(0.1, 0.1, 0.1); + const covariance = StateCovariance.fromSigmas([1.0, 1.0, 1.0, 0.001, 0.001, 0.001], CovarianceFrame.RIC); + const combinedRadius = 0.02 as Kilometers; + + const pc = ProbabilityOfCollision.calculate(relativePosition, relativeVelocity, covariance, combinedRadius); + + expect(pc).toBeGreaterThanOrEqual(0); + expect(pc).toBeLessThanOrEqual(1); + }); + + it('should return higher Pc for smaller miss distances', () => { + const covariance = StateCovariance.fromSigmas([1.0, 1.0, 1.0, 0.001, 0.001, 0.001], CovarianceFrame.RIC); + const relativeVelocity = new Vector3D(0.1, 0.1, 0.1); + const combinedRadius = 0.02 as Kilometers; + + const relativePosition1 = new Vector3D(0.1 as Kilometers, 0.1 as Kilometers, 0.1 as Kilometers); + const relativePosition2 = new Vector3D(1.0 as Kilometers, 1.0 as Kilometers, 1.0 as Kilometers); + + const pc1 = ProbabilityOfCollision.calculate(relativePosition1, relativeVelocity, covariance, combinedRadius); + const pc2 = ProbabilityOfCollision.calculate(relativePosition2, relativeVelocity, covariance, combinedRadius); + + expect(pc1).toBeGreaterThan(pc2); + }); + + it('should combine covariance matrices correctly', () => { + const cov1 = StateCovariance.fromSigmas([1.0, 1.0, 1.0, 0.001, 0.001, 0.001], CovarianceFrame.RIC); + const cov2 = StateCovariance.fromSigmas([0.5, 0.5, 0.5, 0.0005, 0.0005, 0.0005], CovarianceFrame.RIC); + + const combined = ProbabilityOfCollision.combineCovarianceMatrices(cov1, cov2); + + expect(combined).toBeDefined(); + expect(combined.frame).toBe(CovarianceFrame.RIC); + + // Combined variance should be sum of individual variances + expect(combined.matrix.elements[0][0]).toBeCloseTo(cov1.matrix.elements[0][0] + cov2.matrix.elements[0][0]); + }); + }); +});