diff --git a/examples/coordinate-transforms.ts b/examples/coordinate-transforms.ts new file mode 100644 index 00000000..cfba177f --- /dev/null +++ b/examples/coordinate-transforms.ts @@ -0,0 +1,242 @@ +/* eslint-disable no-console */ +/** + * Example demonstrating coordinate transformations. + * + * This example shows: + * - Converting between different coordinate frames (ECI, ECF, LLA) + * - Working with different state vector representations (J2000, TEME, ITRF) + * - Relative coordinates (RIC, Hill) + * - Observation coordinates (RAE, SEZ) + */ + +import { + calcGmst, + Degrees, + ecf2eci, + eci2ecf, + eci2lla, + EpochUTC, + J2000, + Kilometers, + KilometersPerSecond, + lla2ecf, + lla2eci, + Radians, + Satellite, + TEME, + TleLine1, + TleLine2, + Vector3D, +} from '../dist/main.js'; + +// Example 1: ECI ↔ ECF transformations +console.log('=== Example 1: ECI ↔ ECF Transformations ===\n'); + +const date = new Date('2024-01-28T12:00:00.000Z'); +const gmst = calcGmst(date); + +// ECI position vector +const eciPos = { + x: 6778 as Kilometers, + y: 0 as Kilometers, + z: 0 as Kilometers, +}; + +console.log('ECI Position:'); +console.log(` X: ${eciPos.x.toFixed(2)} km`); +console.log(` Y: ${eciPos.y.toFixed(2)} km`); +console.log(` Z: ${eciPos.z.toFixed(2)} km`); + +// Convert to ECF +const ecfPos = eci2ecf(eciPos, gmst.gmst); + +console.log('\nECF Position:'); +console.log(` X: ${ecfPos.x.toFixed(2)} km`); +console.log(` Y: ${ecfPos.y.toFixed(2)} km`); +console.log(` Z: ${ecfPos.z.toFixed(2)} km`); + +// Convert back to ECI +const eciPos2 = ecf2eci(ecfPos, gmst.gmst); + +console.log('\nConverted back to ECI:'); +console.log(` X: ${eciPos2.x.toFixed(2)} km`); +console.log(` Y: ${eciPos2.y.toFixed(2)} km`); +console.log(` Z: ${eciPos2.z.toFixed(2)} km`); + +// Example 2: ECI ↔ LLA transformations +console.log('\n=== Example 2: ECI ↔ Geodetic (LLA) Transformations ===\n'); + +// Convert ECI to Geodetic +const lla = eci2lla(eciPos, gmst.gmst); + +console.log('Geodetic Coordinates:'); +console.log(` Latitude: ${lla.lat.toFixed(4)}°`); +console.log(` Longitude: ${lla.lon.toFixed(4)}°`); +console.log(` Altitude: ${lla.alt.toFixed(2)} km`); + +// Convert geodetic to ECI +const llaRad = { + lat: (lla.lat * (Math.PI / 180)) as Radians, + lon: (lla.lon * (Math.PI / 180)) as Radians, + alt: lla.alt, +}; + +const eciFromLla = lla2eci(llaRad, gmst.gmst); + +console.log('\nConverted back to ECI:'); +console.log(` X: ${eciFromLla.x.toFixed(2)} km`); +console.log(` Y: ${eciFromLla.y.toFixed(2)} km`); +console.log(` Z: ${eciFromLla.z.toFixed(2)} km`); + +// Example 3: Working with different state vector frames +console.log('\n=== Example 3: Different State Vector Frames (J2000, TEME, ITRF) ===\n'); + +const sat = new Satellite({ + tle1: '1 25544U 98067A 24028.54545847 .00031576 00000-0 57240-3 0 9991' as TleLine1, + tle2: '2 25544 51.6418 292.2590 0002595 167.5319 252.0460 15.49326324436741' as TleLine2, +}); + +// Get state in J2000 frame +const j2000State = sat.toJ2000(date); + +console.log('J2000 Frame (ECI):'); +console.log(` Position: [${j2000State.position.x.toFixed(2)}, ${j2000State.position.y.toFixed(2)}, ${j2000State.position.z.toFixed(2)}] km`); +console.log(` Velocity: [${j2000State.velocity.x.toFixed(6)}, ${j2000State.velocity.y.toFixed(6)}, ${j2000State.velocity.z.toFixed(6)}] km/s`); + +// Convert to TEME frame +const temeState = j2000State.toTEME(); + +console.log('\nTEME Frame:'); +console.log(` Position: [${temeState.position.x.toFixed(2)}, ${temeState.position.y.toFixed(2)}, ${temeState.position.z.toFixed(2)}] km`); +console.log(` Velocity: [${temeState.velocity.x.toFixed(6)}, ${temeState.velocity.y.toFixed(6)}, ${temeState.velocity.z.toFixed(6)}] km/s`); + +// Convert to ITRF frame (Earth-fixed) +const itrfState = j2000State.toITRF(); + +console.log('\nITRF Frame (Earth-fixed):'); +console.log(` Position: [${itrfState.position.x.toFixed(2)}, ${itrfState.position.y.toFixed(2)}, ${itrfState.position.z.toFixed(2)}] km`); +console.log(` Velocity: [${itrfState.velocity.x.toFixed(6)}, ${itrfState.velocity.y.toFixed(6)}, ${itrfState.velocity.z.toFixed(6)}] km/s`); + +// Convert ITRF to Geodetic +const geodeticFromItrf = itrfState.toGeodetic(); + +console.log('\nGeodetic from ITRF:'); +console.log(` Latitude: ${geodeticFromItrf.lat.toFixed(4)}°`); +console.log(` Longitude: ${geodeticFromItrf.lon.toFixed(4)}°`); +console.log(` Altitude: ${geodeticFromItrf.alt.toFixed(2)} km`); + +// Example 4: RIC (Radial, In-track, Cross-track) coordinates +console.log('\n=== Example 4: Relative Coordinates (RIC Frame) ===\n'); + +// Create two satellites +const sat1 = new Satellite({ + tle1: '1 25544U 98067A 24028.54545847 .00031576 00000-0 57240-3 0 9991' as TleLine1, + tle2: '2 25544 51.6418 292.2590 0002595 167.5319 252.0460 15.49326324436741' as TleLine2, +}); + +// Create a second satellite slightly ahead in the same orbit +const sat2State = sat1.toJ2000(date); +const sat2Elements = sat2State.toClassicalElements(); + +sat2Elements.trueAnomaly = (sat2Elements.trueAnomaly + 0.1) as Radians; // Slightly ahead + +const sat2J2000 = sat2Elements.toJ2000(); + +console.log('Satellite 1 (Chief) Position:'); +console.log(` [${sat2State.position.x.toFixed(2)}, ${sat2State.position.y.toFixed(2)}, ${sat2State.position.z.toFixed(2)}] km`); + +console.log('\nSatellite 2 (Deputy) Position:'); +console.log(` [${sat2J2000.position.x.toFixed(2)}, ${sat2J2000.position.y.toFixed(2)}, ${sat2J2000.position.z.toFixed(2)}] km`); + +// Convert to RIC coordinates +const ricState = sat2J2000.toRIC(sat2State); + +console.log('\nRelative Position in RIC Frame:'); +console.log(` Radial: ${ricState.position.x.toFixed(4)} km`); +console.log(` In-track: ${ricState.position.y.toFixed(4)} km`); +console.log(` Cross-track: ${ricState.position.z.toFixed(4)} km`); + +console.log('\nRelative Velocity in RIC Frame:'); +console.log(` Radial: ${ricState.velocity.x.toFixed(6)} km/s`); +console.log(` In-track: ${ricState.velocity.y.toFixed(6)} km/s`); +console.log(` Cross-track: ${ricState.velocity.z.toFixed(6)} km/s`); + +// Example 5: Creating state vectors directly +console.log('\n=== Example 5: Creating State Vectors Directly ===\n'); + +const epoch = EpochUTC.fromDateTime(date); + +// Create a J2000 state vector +const customJ2000 = new J2000( + epoch, + new Vector3D( + 6778 as Kilometers, + 0 as Kilometers, + 0 as Kilometers, + ), + new Vector3D( + 0 as KilometersPerSecond, + 7.7 as KilometersPerSecond, + 0 as KilometersPerSecond, + ), +); + +console.log('Custom J2000 State:'); +console.log(` Position: [${customJ2000.position.x.toFixed(2)}, ${customJ2000.position.y.toFixed(2)}, ${customJ2000.position.z.toFixed(2)}] km`); +console.log(` Velocity: [${customJ2000.velocity.x.toFixed(6)}, ${customJ2000.velocity.y.toFixed(6)}, ${customJ2000.velocity.z.toFixed(6)}] km/s`); + +// Convert to TEME +const customTEME = new TEME( + epoch, + new Vector3D( + 6778 as Kilometers, + 0 as Kilometers, + 0 as Kilometers, + ), + new Vector3D( + 0 as KilometersPerSecond, + 7.7 as KilometersPerSecond, + 0 as KilometersPerSecond, + ), +); + +console.log('\nCustom TEME State:'); +console.log(` Position: [${customTEME.position.x.toFixed(2)}, ${customTEME.position.y.toFixed(2)}, ${customTEME.position.z.toFixed(2)}] km`); +console.log(` Velocity: [${customTEME.velocity.x.toFixed(6)}, ${customTEME.velocity.y.toFixed(6)}, ${customTEME.velocity.z.toFixed(6)}] km/s`); + +// Convert TEME to J2000 +const temeToJ2000 = customTEME.toJ2000(); + +console.log('\nTEME converted to J2000:'); +console.log(` Position: [${temeToJ2000.position.x.toFixed(2)}, ${temeToJ2000.position.y.toFixed(2)}, ${temeToJ2000.position.z.toFixed(2)}] km`); +console.log(` Velocity: [${temeToJ2000.velocity.x.toFixed(6)}, ${temeToJ2000.velocity.y.toFixed(6)}, ${temeToJ2000.velocity.z.toFixed(6)}] km/s`); + +// Example 6: LLA to ECF transformation +console.log('\n=== Example 6: Geodetic to ECF ===\n'); + +const observerLla = { + lat: 41.754785 as Degrees, + lon: -70.539151 as Degrees, + alt: 0.060966 as Kilometers, +}; + +console.log('Observer Location:'); +console.log(` Latitude: ${observerLla.lat}°`); +console.log(` Longitude: ${observerLla.lon}°`); +console.log(` Altitude: ${observerLla.alt} km`); + +const observerEcf = lla2ecf(observerLla); + +console.log('\nECF Position:'); +console.log(` X: ${observerEcf.x.toFixed(4)} km`); +console.log(` Y: ${observerEcf.y.toFixed(4)} km`); +console.log(` Z: ${observerEcf.z.toFixed(4)} km`); + +const distance = Math.sqrt( + observerEcf.x * observerEcf.x + + observerEcf.y * observerEcf.y + + observerEcf.z * observerEcf.z, +); + +console.log(`\nDistance from Earth center: ${distance.toFixed(4)} km`); +console.log(`Earth equatorial radius: 6378.137 km`); diff --git a/examples/doppler.ts b/examples/doppler.ts new file mode 100644 index 00000000..984492a4 --- /dev/null +++ b/examples/doppler.ts @@ -0,0 +1,234 @@ +/* eslint-disable no-console */ +/** + * Example demonstrating Doppler shift calculations. + * + * This example shows: + * - Calculating Doppler shift for satellite communications + * - Computing frequency shifts due to relative motion + * - Tracking Doppler changes during a satellite pass + * - Applying Doppler corrections + */ + +import { + Degrees, + dopplerFactor, + Kilometers, + Satellite, + Sensor, + TleLine1, + TleLine2, +} from '../dist/main.js'; + +// Example 1: Basic Doppler calculation +console.log('=== Example 1: Basic Doppler Shift Calculation ===\n'); + +// Create ground station +const groundStation = new Sensor({ + lat: 41.754785 as Degrees, + lon: -70.539151 as Degrees, + alt: 0.060966 as Kilometers, + name: 'Cape Cod Ground Station', +}); + +// Create ISS satellite +const iss = new Satellite({ + tle1: '1 25544U 98067A 24028.54545847 .00031576 00000-0 57240-3 0 9991' as TleLine1, + tle2: '2 25544 51.6418 292.2590 0002595 167.5319 252.0460 15.49326324436741' as TleLine2, +}); + +const transmitFreq = 437.8e6; // 437.8 MHz (UHF amateur radio frequency) +const date = new Date('2024-01-28T12:00:00.000Z'); + +console.log(`Ground Station: ${groundStation.name}`); +console.log(`Transmit Frequency: ${(transmitFreq / 1e6).toFixed(1)} MHz`); +console.log(`Time: ${date.toISOString()}\n`); + +// Calculate Doppler shift +const doppler = iss.dopplerFactor(groundStation, date); +const receivedFreq = iss.applyDoppler(transmitFreq, groundStation, date); +const freqShift = receivedFreq - transmitFreq; + +console.log('Doppler Calculations:'); +console.log(` Doppler Factor: ${doppler.toFixed(8)}`); +console.log(` Transmit Frequency: ${(transmitFreq / 1e6).toFixed(4)} MHz`); +console.log(` Received Frequency: ${(receivedFreq / 1e6).toFixed(4)} MHz`); +console.log(` Frequency Shift: ${(freqShift / 1e3).toFixed(2)} kHz`); + +// Get look angles for context +const rae = groundStation.rae(iss, date); + +console.log(`\nSatellite Position:`); +console.log(` Azimuth: ${rae.az.toFixed(1)}°`); +console.log(` Elevation: ${rae.el.toFixed(1)}°`); +console.log(` Range: ${rae.rng.toFixed(1)} km`); + +// Example 2: Doppler during a satellite pass +console.log('\n=== Example 2: Doppler Shift During a Pass ===\n'); + +// Simulate a pass over 10 minutes +const passStart = new Date('2024-01-28T12:00:00.000Z'); +const interval = 60; // seconds + +console.log('Time El Range Doppler Freq Shift'); +console.log('──────── ─── ──────── ────────── ────────────'); + +for (let i = 0; i <= 10; i++) { + const time = new Date(passStart.getTime() + i * interval * 1000); + const timeStr = time.toISOString().substring(11, 19); + + const passRae = groundStation.rae(iss, time); + const passDoppler = iss.dopplerFactor(groundStation, time); + const passReceivedFreq = iss.applyDoppler(transmitFreq, groundStation, time); + const passFreqShift = passReceivedFreq - transmitFreq; + + const elStr = passRae.el.toFixed(1).padStart(5); + const rngStr = passRae.rng.toFixed(1).padStart(8); + const dopplerStr = passDoppler.toFixed(8); + const shiftStr = (passFreqShift / 1e3).toFixed(2).padStart(9); + + console.log(`${timeStr} ${elStr}° ${rngStr} km ${dopplerStr} ${shiftStr} kHz`); +} + +// Example 3: Different frequency bands +console.log('\n=== Example 3: Doppler Shift Across Different Bands ===\n'); + +const frequencies = [ + { band: 'VHF', freq: 145.8e6, name: '145.8 MHz' }, + { band: 'UHF', freq: 437.8e6, name: '437.8 MHz' }, + { band: 'L-band', freq: 1575.42e6, name: '1575.42 MHz (GPS L1)' }, + { band: 'S-band', freq: 2200e6, name: '2.2 GHz' }, + { band: 'X-band', freq: 8400e6, name: '8.4 GHz' }, + { band: 'Ku-band', freq: 12e9, name: '12 GHz' }, +]; + +const bandCheckTime = new Date('2024-01-28T12:05:00.000Z'); + +console.log(`Time: ${bandCheckTime.toISOString()}\n`); +console.log('Band Frequency Received Freq Shift'); +console.log('──────── ────────────── ────────────── ─────────'); + +frequencies.forEach((f) => { + const bandReceivedFreq = iss.applyDoppler(f.freq, groundStation, bandCheckTime); + const bandFreqShift = bandReceivedFreq - f.freq; + + const bandStr = f.band.padEnd(8); + const freqStr = f.name.padEnd(14); + const recvStr = formatFrequency(bandReceivedFreq).padEnd(14); + const shiftStr = formatFrequency(Math.abs(bandFreqShift)).padStart(9); + + console.log(`${bandStr} ${freqStr} ${recvStr} ${bandFreqShift >= 0 ? '+' : '-'}${shiftStr}`); +}); + +// Example 4: Maximum Doppler shift +console.log('\n=== Example 4: Maximum Doppler Shift ===\n'); + +// For a LEO satellite, maximum Doppler occurs when satellite velocity +// is directly toward or away from the observer + +// Calculate approximate maximum radial velocity +// ISS orbital velocity ~ 7.66 km/s +const orbitalVelocity = 7.66; // km/s +const maxRadialVelocity = orbitalVelocity; // km/s (worst case) + +// Calculate theoretical maximum Doppler factor +const speedOfLight = 299792.458; // km/s +const maxDopplerFactor = maxRadialVelocity / speedOfLight; + +console.log(`ISS Orbital Velocity: ~${orbitalVelocity} km/s`); +console.log(`Speed of Light: ${speedOfLight} km/s`); +console.log(`Max Doppler Factor: ±${(maxDopplerFactor * 100).toFixed(6)}%`); + +console.log('\nTheoretical Maximum Frequency Shifts:'); + +frequencies.slice(0, 4).forEach((f) => { + const maxShift = f.freq * maxDopplerFactor; + + console.log(` ${f.name}: ±${formatFrequency(maxShift)}`); +}); + +// Example 5: Doppler rate of change +console.log('\n=== Example 5: Doppler Rate of Change ===\n'); + +const rateStart = new Date('2024-01-28T12:00:00.000Z'); +const deltaTime = 10; // seconds + +console.log('Measuring Doppler rate of change over time:\n'); + +for (let i = 0; i < 5; i++) { + const t1 = new Date(rateStart.getTime() + i * 60 * 1000); + const t2 = new Date(t1.getTime() + deltaTime * 1000); + + const freq1 = iss.applyDoppler(transmitFreq, groundStation, t1); + const freq2 = iss.applyDoppler(transmitFreq, groundStation, t2); + + const freqChange = freq2 - freq1; + const rateOfChange = freqChange / deltaTime; // Hz per second + + const timeStr = t1.toISOString().substring(11, 19); + + console.log(`At ${timeStr}:`); + console.log(` Frequency: ${(freq1 / 1e6).toFixed(4)} MHz`); + console.log(` Rate of change: ${rateOfChange.toFixed(2)} Hz/s`); + console.log(''); +} + +// Example 6: Using dopplerFactor function directly +console.log('=== Example 6: Using dopplerFactor Utility Function ===\n'); + +const observerVelocity = [0, 0, 0]; // Ground station (stationary in ECEF) +const satelliteState = iss.eci(date); +const satelliteVelocity = [ + satelliteState.velocity.x, + satelliteState.velocity.y, + satelliteState.velocity.z, +]; + +const observerState = groundStation.eci(date); +const observerPosition = [ + observerState.position.x, + observerState.position.y, + observerState.position.z, +]; +const satellitePosition = [ + satelliteState.position.x, + satelliteState.position.y, + satelliteState.position.z, +]; + +// Calculate relative position vector +const relativePos = [ + satellitePosition[0] - observerPosition[0], + satellitePosition[1] - observerPosition[1], + satellitePosition[2] - observerPosition[2], +]; + +// Calculate relative velocity +const relativeVel = [ + satelliteVelocity[0] - observerVelocity[0], + satelliteVelocity[1] - observerVelocity[1], + satelliteVelocity[2] - observerVelocity[2], +]; + +console.log('Position Vector (km):'); +console.log(` [${relativePos[0].toFixed(2)}, ${relativePos[1].toFixed(2)}, ${relativePos[2].toFixed(2)}]`); + +console.log('\nVelocity Vector (km/s):'); +console.log(` [${relativeVel[0].toFixed(4)}, ${relativeVel[1].toFixed(4)}, ${relativeVel[2].toFixed(4)}]`); + +const manualDoppler = dopplerFactor(relativePos, relativeVel); + +console.log(`\nCalculated Doppler Factor: ${manualDoppler.toFixed(8)}`); +console.log(`Satellite method result: ${doppler.toFixed(8)}`); + +// Helper function to format frequencies +function formatFrequency(freq: number): string { + if (freq >= 1e9) { + return `${(freq / 1e9).toFixed(4)} GHz`; + } else if (freq >= 1e6) { + return `${(freq / 1e6).toFixed(4)} MHz`; + } else if (freq >= 1e3) { + return `${(freq / 1e3).toFixed(2)} kHz`; + } else { + return `${freq.toFixed(2)} Hz`; + } +} diff --git a/examples/maneuvers.ts b/examples/maneuvers.ts new file mode 100644 index 00000000..e99ecb99 --- /dev/null +++ b/examples/maneuvers.ts @@ -0,0 +1,186 @@ +/* eslint-disable no-console */ +/** + * Example demonstrating orbital maneuver calculations. + * + * This example shows: + * - Hohmann transfer calculations + * - Two-burn orbit transfers + * - Delta-V calculations + * - Transfer orbit parameters + */ + +import { + ClassicalElements, + Degrees, + EpochUTC, + Kilometers, + TwoBurnOrbitTransfer, +} from '../dist/main.js'; + +console.log('=== Example 1: Hohmann Transfer from LEO to GEO ===\n'); + +const date = new Date('2024-01-28T00:00:00.000Z'); + +// Create initial LEO orbit (circular) +const leoOrbit = new ClassicalElements({ + epoch: EpochUTC.fromDateTime(date), + semimajorAxis: 6778 as Kilometers, // ~400 km altitude + eccentricity: 0.001, // Nearly circular + inclination: 28.5 as Degrees, + rightAscension: 0 as Degrees, + argPerigee: 0 as Degrees, + trueAnomaly: 0 as Degrees, +}); + +// Create target GEO orbit (circular) +const geoOrbit = new ClassicalElements({ + epoch: EpochUTC.fromDateTime(date), + semimajorAxis: 42164 as Kilometers, // GEO altitude + eccentricity: 0.001, + inclination: 0 as Degrees, // Equatorial + rightAscension: 0 as Degrees, + argPerigee: 0 as Degrees, + trueAnomaly: 0 as Degrees, +}); + +console.log('Initial Orbit (LEO):'); +console.log(` Altitude: ${(leoOrbit.semimajorAxis - 6378.137).toFixed(2)} km`); +console.log(` Period: ${leoOrbit.period.toFixed(2)} minutes`); +console.log(` Velocity: ${Math.sqrt(398600.4418 / leoOrbit.semimajorAxis).toFixed(3)} km/s`); + +console.log('\nTarget Orbit (GEO):'); +console.log(` Altitude: ${(geoOrbit.semimajorAxis - 6378.137).toFixed(2)} km`); +console.log(` Period: ${geoOrbit.period.toFixed(2)} minutes`); +console.log(` Velocity: ${Math.sqrt(398600.4418 / geoOrbit.semimajorAxis).toFixed(3)} km/s`); + +// Calculate Hohmann transfer +const transfer = new TwoBurnOrbitTransfer(leoOrbit, geoOrbit); +const hohmann = transfer.hohmannTransfer(); + +console.log('\nHohmann Transfer:'); +console.log(` Transfer semi-major axis: ${hohmann.semimajorAxis.toFixed(2)} km`); +console.log(` Transfer eccentricity: ${hohmann.eccentricity.toFixed(6)}`); +console.log(` Transfer period: ${hohmann.period.toFixed(2)} minutes`); +console.log(` Transfer time: ${(hohmann.period / 2).toFixed(2)} minutes (half orbit)`); + +// Calculate transfer orbit periapsis and apoapsis +const transferPeriapsis = hohmann.semimajorAxis * (1 - hohmann.eccentricity); +const transferApoapsis = hohmann.semimajorAxis * (1 + hohmann.eccentricity); + +console.log(` Periapsis altitude: ${(transferPeriapsis - 6378.137).toFixed(2)} km`); +console.log(` Apoapsis altitude: ${(transferApoapsis - 6378.137).toFixed(2)} km`); + +// Calculate delta-V requirements +const v1 = Math.sqrt(398600.4418 / leoOrbit.semimajorAxis); // LEO velocity +const vTransferPerigee = Math.sqrt(398600.4418 * (2 / leoOrbit.semimajorAxis - 1 / hohmann.semimajorAxis)); +const deltaV1 = vTransferPerigee - v1; + +const v2 = Math.sqrt(398600.4418 / geoOrbit.semimajorAxis); // GEO velocity +const vTransferApogee = Math.sqrt(398600.4418 * (2 / geoOrbit.semimajorAxis - 1 / hohmann.semimajorAxis)); +const deltaV2 = v2 - vTransferApogee; + +const totalDeltaV = deltaV1 + deltaV2; + +console.log('\nDelta-V Requirements:'); +console.log(` First burn (at LEO): ${deltaV1.toFixed(3)} km/s`); +console.log(` Second burn (at GEO): ${deltaV2.toFixed(3)} km/s`); +console.log(` Total delta-V: ${totalDeltaV.toFixed(3)} km/s`); + +// Example 2: Transfer between elliptical orbits +console.log('\n=== Example 2: Transfer from LEO to Molniya Orbit ===\n'); + +const molniyaOrbit = new ClassicalElements({ + epoch: EpochUTC.fromDateTime(date), + semimajorAxis: 26554 as Kilometers, + eccentricity: 0.74, + inclination: 63.4 as Degrees, + rightAscension: 0 as Degrees, + argPerigee: 270 as Degrees, + trueAnomaly: 0 as Degrees, +}); + +console.log('Target Orbit (Molniya):'); +console.log(` Semi-major axis: ${molniyaOrbit.semimajorAxis.toFixed(2)} km`); +console.log(` Eccentricity: ${molniyaOrbit.eccentricity.toFixed(4)}`); +console.log(` Apogee altitude: ${((molniyaOrbit.semimajorAxis * (1 + molniyaOrbit.eccentricity)) - 6378.137).toFixed(2)} km`); +console.log(` Perigee altitude: ${((molniyaOrbit.semimajorAxis * (1 - molniyaOrbit.eccentricity)) - 6378.137).toFixed(2)} km`); +console.log(` Period: ${molniyaOrbit.period.toFixed(2)} minutes`); + +const molniyaTransfer = new TwoBurnOrbitTransfer(leoOrbit, molniyaOrbit); +const molniyaHohmann = molniyaTransfer.hohmannTransfer(); + +console.log('\nTransfer Orbit:'); +console.log(` Semi-major axis: ${molniyaHohmann.semimajorAxis.toFixed(2)} km`); +console.log(` Eccentricity: ${molniyaHohmann.eccentricity.toFixed(4)}`); +console.log(` Period: ${molniyaHohmann.period.toFixed(2)} minutes`); +console.log(` Transfer time: ${(molniyaHohmann.period / 2).toFixed(2)} minutes`); + +// Example 3: Orbit raising maneuver +console.log('\n=== Example 3: Simple Orbit Raising (LEO to MEO) ===\n'); + +const meoOrbit = new ClassicalElements({ + epoch: EpochUTC.fromDateTime(date), + semimajorAxis: 26560 as Kilometers, // GPS-like orbit + eccentricity: 0.001, + inclination: 55 as Degrees, + rightAscension: 0 as Degrees, + argPerigee: 0 as Degrees, + trueAnomaly: 0 as Degrees, +}); + +console.log('Initial Orbit (LEO):'); +console.log(` Altitude: ${(leoOrbit.semimajorAxis - 6378.137).toFixed(2)} km`); +console.log(` Period: ${leoOrbit.period.toFixed(2)} minutes`); + +console.log('\nTarget Orbit (MEO - GPS-like):'); +console.log(` Altitude: ${(meoOrbit.semimajorAxis - 6378.137).toFixed(2)} km`); +console.log(` Period: ${meoOrbit.period.toFixed(2)} minutes`); + +const meoTransfer = new TwoBurnOrbitTransfer(leoOrbit, meoOrbit); +const meoHohmann = meoTransfer.hohmannTransfer(); + +// Calculate delta-V for this transfer +const vLeo = Math.sqrt(398600.4418 / leoOrbit.semimajorAxis); +const vMeoTransferPerigee = Math.sqrt(398600.4418 * (2 / leoOrbit.semimajorAxis - 1 / meoHohmann.semimajorAxis)); +const meoDeltaV1 = vMeoTransferPerigee - vLeo; + +const vMeo = Math.sqrt(398600.4418 / meoOrbit.semimajorAxis); +const vMeoTransferApogee = Math.sqrt(398600.4418 * (2 / meoOrbit.semimajorAxis - 1 / meoHohmann.semimajorAxis)); +const meoDeltaV2 = vMeo - vMeoTransferApogee; + +console.log('\nTransfer Parameters:'); +console.log(` Transfer semi-major axis: ${meoHohmann.semimajorAxis.toFixed(2)} km`); +console.log(` Transfer time: ${(meoHohmann.period / 2).toFixed(2)} minutes`); +console.log(` First burn delta-V: ${meoDeltaV1.toFixed(3)} km/s`); +console.log(` Second burn delta-V: ${meoDeltaV2.toFixed(3)} km/s`); +console.log(` Total delta-V: ${(meoDeltaV1 + meoDeltaV2).toFixed(3)} km/s`); + +// Example 4: Comparison of different transfers +console.log('\n=== Example 4: Comparison of Common Transfers ===\n'); + +const transfers = [ + { name: 'LEO to GEO', initial: 6778, final: 42164 }, + { name: 'LEO to MEO (GPS)', initial: 6778, final: 26560 }, + { name: 'LEO to ISS-like', initial: 6678, final: 6778 }, +]; + +transfers.forEach((t) => { + const r1 = t.initial as Kilometers; + const r2 = t.final as Kilometers; + const at = (r1 + r2) / 2; // Transfer orbit semi-major axis + + const v1 = Math.sqrt(398600.4418 / r1); + const vt1 = Math.sqrt(398600.4418 * (2 / r1 - 1 / at)); + const dv1 = Math.abs(vt1 - v1); + + const v2 = Math.sqrt(398600.4418 / r2); + const vt2 = Math.sqrt(398600.4418 * (2 / r2 - 1 / at)); + const dv2 = Math.abs(v2 - vt2); + + const period = 2 * Math.PI * Math.sqrt(Math.pow(at, 3) / 398600.4418) / 60; // minutes + + console.log(`${t.name}:`); + console.log(` Total delta-V: ${(dv1 + dv2).toFixed(3)} km/s`); + console.log(` Transfer time: ${(period / 2).toFixed(2)} minutes`); + console.log(''); +}); diff --git a/examples/moon.ts b/examples/moon.ts new file mode 100644 index 00000000..65c038e9 --- /dev/null +++ b/examples/moon.ts @@ -0,0 +1,178 @@ +/* eslint-disable no-console */ +/** + * Example demonstrating lunar position calculations. + * + * This example shows: + * - Getting Moon position in different coordinate frames + * - Calculating Moon rise/set times + * - Computing lunar distance and angular diameter + * - Finding Moon phase and illumination + */ + +import { + Degrees, + eci2lla, + Kilometers, + Meters, + Moon, + calcGmst, +} from '../dist/main.js'; + +// Example 1: Get Moon position at a specific time +console.log('=== Example 1: Moon Position ===\n'); + +const date = new Date('2024-01-28T12:00:00.000Z'); + +// Get Moon position in ECI coordinates +const moonEci = Moon.eci(date); + +console.log('Moon Position (ECI):'); +console.log(` X: ${moonEci.position.x.toFixed(2)} km`); +console.log(` Y: ${moonEci.position.y.toFixed(2)} km`); +console.log(` Z: ${moonEci.position.z.toFixed(2)} km`); + +// Calculate distance from Earth +const distance = Math.sqrt( + moonEci.position.x * moonEci.position.x + + moonEci.position.y * moonEci.position.y + + moonEci.position.z * moonEci.position.z, +); + +console.log(`\nDistance from Earth: ${distance.toFixed(2)} km`); +console.log(` ${(distance / 384400).toFixed(4)} × mean lunar distance`); + +// Convert to latitude/longitude (sub-lunar point) +const gmst = calcGmst(date); +const moonLla = eci2lla(moonEci.position, gmst.gmst); + +console.log('\nSub-Lunar Point:'); +console.log(` Latitude: ${moonLla.lat.toFixed(4)}°`); +console.log(` Longitude: ${moonLla.lon.toFixed(4)}°`); + +// Example 2: Moon rise and set times +console.log('\n=== Example 2: Moon Rise/Set Times ===\n'); + +const observerLat = 41 as Degrees; +const observerLon = -71 as Degrees; +const observerAlt = 0 as Meters; + +const moonTimes = Moon.getTimes(date, observerLat, observerLon, observerAlt); + +console.log(`Observer Location: ${observerLat}° N, ${Math.abs(observerLon)}° W`); +console.log(`\nMoon Times for ${date.toDateString()}:`); + +if (moonTimes.rise) { + console.log(` Moonrise: ${moonTimes.rise.toLocaleTimeString()}`); +} else { + console.log(' Moonrise: No moonrise today'); +} + +if (moonTimes.set) { + console.log(` Moonset: ${moonTimes.set.toLocaleTimeString()}`); +} else { + console.log(' Moonset: No moonset today'); +} + +// Example 3: Moon phase and illumination +console.log('\n=== Example 3: Moon Phase and Illumination ===\n'); + +const illumination = Moon.getIllumination(date); + +console.log(`Phase: ${(illumination.phase * 100).toFixed(2)}%`); +console.log(` 0% = New Moon`); +console.log(` 25% = First Quarter`); +console.log(` 50% = Full Moon`); +console.log(` 75% = Last Quarter`); + +console.log(`\nIllumination: ${(illumination.fraction * 100).toFixed(2)}%`); +console.log(`Phase Angle: ${(illumination.phaseAngle * (180 / Math.PI)).toFixed(2)}°`); + +// Determine moon phase name +let phaseName = ''; + +if (illumination.phase < 0.03 || illumination.phase > 0.97) { + phaseName = 'New Moon'; +} else if (illumination.phase < 0.22) { + phaseName = 'Waxing Crescent'; +} else if (illumination.phase < 0.28) { + phaseName = 'First Quarter'; +} else if (illumination.phase < 0.47) { + phaseName = 'Waxing Gibbous'; +} else if (illumination.phase < 0.53) { + phaseName = 'Full Moon'; +} else if (illumination.phase < 0.72) { + phaseName = 'Waning Gibbous'; +} else if (illumination.phase < 0.78) { + phaseName = 'Last Quarter'; +} else { + phaseName = 'Waning Crescent'; +} + +console.log(`Moon Phase: ${phaseName}`); + +// Example 4: Angular diameter +console.log('\n=== Example 4: Moon Angular Size ===\n'); + +// Mean lunar radius is 1737.4 km +const lunarRadius = 1737.4 as Kilometers; + +// Calculate angular diameter in radians +const angularDiameterRad = 2 * Math.atan(lunarRadius / distance); + +// Convert to degrees and arcminutes +const angularDiameterDeg = angularDiameterRad * (180 / Math.PI); +const angularDiameterArcmin = angularDiameterDeg * 60; + +console.log(`Angular Diameter: ${angularDiameterArcmin.toFixed(2)} arcminutes`); +console.log(` ${angularDiameterDeg.toFixed(4)}°`); +console.log(`\nMean angular diameter: ~31 arcminutes`); + +// Example 5: Moon position over a day +console.log('\n=== Example 5: Moon Position Every 6 Hours ===\n'); + +for (let hour = 0; hour < 24; hour += 6) { + const timePoint = new Date(date.getTime()); + + timePoint.setUTCHours(hour, 0, 0, 0); + + const moonPos = Moon.eci(timePoint); + const dist = Math.sqrt( + moonPos.position.x * moonPos.position.x + + moonPos.position.y * moonPos.position.y + + moonPos.position.z * moonPos.position.z, + ); + + const gmstTime = calcGmst(timePoint); + const lla = eci2lla(moonPos.position, gmstTime.gmst); + + console.log(`${timePoint.toUTCString()}:`); + console.log(` Distance: ${dist.toFixed(2)} km`); + console.log(` Sub-lunar point: ${lla.lat.toFixed(2)}° N, ${lla.lon.toFixed(2)}° E`); + console.log(''); +} + +// Example 6: Find next full moon (approximate) +console.log('=== Example 6: Finding Next Full Moon (Approximate) ===\n'); + +const startDate = new Date('2024-01-01T00:00:00.000Z'); +let fullMoonDate: Date | null = null; + +// Check every day for a month +for (let day = 0; day < 60; day++) { + const checkDate = new Date(startDate.getTime() + day * 24 * 60 * 60 * 1000); + const illum = Moon.getIllumination(checkDate); + + // Full moon is when phase is closest to 0.5 + if (Math.abs(illum.phase - 0.5) < 0.02) { + fullMoonDate = checkDate; + console.log(`Approximate Full Moon: ${fullMoonDate.toDateString()}`); + console.log(` Phase: ${(illum.phase * 100).toFixed(2)}%`); + console.log(` Illumination: ${(illum.fraction * 100).toFixed(2)}%`); + + break; + } +} + +if (!fullMoonDate) { + console.log('No full moon found in search period'); +} diff --git a/examples/observations.ts b/examples/observations.ts new file mode 100644 index 00000000..195fe06f --- /dev/null +++ b/examples/observations.ts @@ -0,0 +1,237 @@ +/* eslint-disable no-console */ +/** + * Example demonstrating different observation formats. + * + * This example shows: + * - Right Ascension and Declination (RADEC) observations + * - Converting between observation formats + * - Topocentric vs Geocentric observations + * - Creating state vectors from observations + */ + +import { + Degrees, + EpochUTC, + Kilometers, + RadecGeocentric, + RadecTopocentric, + Radians, + Satellite, + Sensor, + TleLine1, + TleLine2, +} from '../dist/main.js'; + +// Example 1: Topocentric RADEC observation +console.log('=== Example 1: Topocentric RADEC Observations ===\n'); + +// Create ground station +const observatory = new Sensor({ + lat: 34.5 as Degrees, // Example: Southern California + lon: -117.9 as Degrees, + alt: 1.2 as Kilometers, + name: 'Example Observatory', +}); + +// Create satellite +const satellite = new Satellite({ + tle1: '1 25544U 98067A 24028.54545847 .00031576 00000-0 57240-3 0 9991' as TleLine1, + tle2: '2 25544 51.6418 292.2590 0002595 167.5319 252.0460 15.49326324436741' as TleLine2, +}); + +const observationTime = new Date('2024-01-28T12:00:00.000Z'); +const epoch = EpochUTC.fromDateTime(observationTime); + +console.log(`Observatory: ${observatory.name}`); +console.log(` Location: ${observatory.lat}° N, ${Math.abs(observatory.lon)}° W`); +console.log(` Altitude: ${observatory.alt} km`); +console.log(`\nObservation Time: ${observationTime.toISOString()}\n`); + +// Get satellite state +const satState = satellite.toJ2000(observationTime); + +// Create topocentric observation +const topoRadec = RadecTopocentric.fromStateVector(satState, observatory); + +console.log('Topocentric Observation (from ground station):'); +console.log(` Right Ascension: ${topoRadec.rightAscension.toFixed(6)} rad`); +console.log(` ${(topoRadec.rightAscension * (180 / Math.PI)).toFixed(4)}°`); +console.log(` ${raToHMS(topoRadec.rightAscension)}`); + +console.log(` Declination: ${topoRadec.declination.toFixed(6)} rad`); +console.log(` ${(topoRadec.declination * (180 / Math.PI)).toFixed(4)}°`); +console.log(` ${decToDMS(topoRadec.declination)}`); + +console.log(` Range: ${topoRadec.range.toFixed(2)} km`); + +// Example 2: Geocentric RADEC observation +console.log('\n=== Example 2: Geocentric RADEC Observations ===\n'); + +// Create geocentric observation (from Earth's center) +const geoRadec = RadecGeocentric.fromStateVector(satState); + +console.log('Geocentric Observation (from Earth center):'); +console.log(` Right Ascension: ${geoRadec.rightAscension.toFixed(6)} rad`); +console.log(` ${(geoRadec.rightAscension * (180 / Math.PI)).toFixed(4)}°`); +console.log(` ${raToHMS(geoRadec.rightAscension)}`); + +console.log(` Declination: ${geoRadec.declination.toFixed(6)} rad`); +console.log(` ${(geoRadec.declination * (180 / Math.PI)).toFixed(4)}°`); +console.log(` ${decToDMS(geoRadec.declination)}`); + +console.log(` Range: ${geoRadec.range.toFixed(2)} km`); + +// Example 3: Compare different observation formats +console.log('\n=== Example 3: Comparing Observation Formats ===\n'); + +// Get RAE (Range, Azimuth, Elevation) for comparison +const rae = observatory.rae(satellite, observationTime); + +console.log('Same satellite observed in different coordinate systems:\n'); + +console.log('RAE (Range-Azimuth-Elevation):'); +console.log(` Azimuth: ${rae.az.toFixed(4)}°`); +console.log(` Elevation: ${rae.el.toFixed(4)}°`); +console.log(` Range: ${rae.rng.toFixed(2)} km`); + +console.log('\nTopocentric RADEC:'); +console.log(` RA: ${raToHMS(topoRadec.rightAscension)}`); +console.log(` Dec: ${decToDMS(topoRadec.declination)}`); +console.log(` Range: ${topoRadec.range.toFixed(2)} km`); + +console.log('\nGeocentric RADEC:'); +console.log(` RA: ${raToHMS(geoRadec.rightAscension)}`); +console.log(` Dec: ${decToDMS(geoRadec.declination)}`); +console.log(` Range: ${geoRadec.range.toFixed(2)} km`); + +// Example 4: Creating observations from angles +console.log('\n=== Example 4: Creating RADEC from Angles ===\n'); + +// Create a topocentric observation manually +const manualTopoRadec = new RadecTopocentric( + epoch, + 1.5 as Radians, // Right ascension + 0.5 as Radians, // Declination + 1200 as Kilometers, // Range + observatory, +); + +console.log('Manually created topocentric observation:'); +console.log(` RA: ${raToHMS(manualTopoRadec.rightAscension)}`); +console.log(` Dec: ${decToDMS(manualTopoRadec.declination)}`); +console.log(` Range: ${manualTopoRadec.range.toFixed(2)} km`); + +// Convert to state vector +const stateFromRadec = manualTopoRadec.toJ2000(); + +console.log('\nConverted to J2000 State Vector:'); +console.log(` Position: [${stateFromRadec.position.x.toFixed(2)}, ${stateFromRadec.position.y.toFixed(2)}, ${stateFromRadec.position.z.toFixed(2)}] km`); + +// Example 5: Tracking object across the sky +console.log('\n=== Example 5: Tracking Satellite Motion in RADEC ===\n'); + +console.log('Time Right Ascension Declination Range'); +console.log('──────── ───────────────── ───────────── ─────────'); + +for (let i = 0; i < 6; i++) { + const trackTime = new Date(observationTime.getTime() + i * 5 * 60 * 1000); + const trackState = satellite.toJ2000(trackTime); + const trackRadec = RadecTopocentric.fromStateVector(trackState, observatory); + + const timeStr = trackTime.toISOString().substring(11, 19); + const raStr = raToHMS(trackRadec.rightAscension); + const decStr = decToDMS(trackRadec.declination); + const rngStr = trackRadec.range.toFixed(1).padStart(9); + + console.log(`${timeStr} ${raStr} ${decStr} ${rngStr} km`); +} + +// Example 6: Geocentric observations for different satellites +console.log('\n=== Example 6: Multiple Satellites in Geocentric RADEC ===\n'); + +const satellites = [ + { + name: 'ISS', + sat: new Satellite({ + tle1: '1 25544U 98067A 24028.54545847 .00031576 00000-0 57240-3 0 9991' as TleLine1, + tle2: '2 25544 51.6418 292.2590 0002595 167.5319 252.0460 15.49326324436741' as TleLine2, + }), + }, + { + name: 'Hubble', + sat: new Satellite({ + tle1: '1 20580U 90037B 24028.50123227 .00000825 00000-0 39644-4 0 9997' as TleLine1, + tle2: '2 20580 28.4696 273.2640 0002975 297.7865 189.2151 15.09696656316758' as TleLine2, + }), + }, +]; + +console.log(`Time: ${observationTime.toISOString()}\n`); +console.log('Satellite Right Ascension Declination Range'); +console.log('───────── ───────────────── ───────────── ─────────'); + +satellites.forEach((satInfo) => { + const satJ2000 = satInfo.sat.toJ2000(observationTime); + const satGeoRadec = RadecGeocentric.fromStateVector(satJ2000); + + const nameStr = satInfo.name.padEnd(9); + const raStr = raToHMS(satGeoRadec.rightAscension); + const decStr = decToDMS(satGeoRadec.declination); + const rngStr = satGeoRadec.range.toFixed(1).padStart(9); + + console.log(`${nameStr} ${raStr} ${decStr} ${rngStr} km`); +}); + +// Example 7: Angular separation between objects +console.log('\n=== Example 7: Angular Separation Between Objects ===\n'); + +const sat1State = satellites[0].sat.toJ2000(observationTime); +const sat2State = satellites[1].sat.toJ2000(observationTime); + +const sat1Radec = RadecGeocentric.fromStateVector(sat1State); +const sat2Radec = RadecGeocentric.fromStateVector(sat2State); + +// Calculate angular separation using spherical trigonometry +const ra1 = sat1Radec.rightAscension; +const dec1 = sat1Radec.declination; +const ra2 = sat2Radec.rightAscension; +const dec2 = sat2Radec.declination; + +const angularSep = Math.acos( + Math.sin(dec1) * Math.sin(dec2) + + Math.cos(dec1) * Math.cos(dec2) * Math.cos(ra1 - ra2), +); + +console.log(`${satellites[0].name}:`); +console.log(` RA: ${raToHMS(ra1)}`); +console.log(` Dec: ${decToDMS(dec1)}`); + +console.log(`\n${satellites[1].name}:`); +console.log(` RA: ${raToHMS(ra2)}`); +console.log(` Dec: ${decToDMS(dec2)}`); + +console.log(`\nAngular Separation:`); +console.log(` ${(angularSep * (180 / Math.PI)).toFixed(4)}°`); +console.log(` ${(angularSep * (180 / Math.PI) * 60).toFixed(2)} arcminutes`); + +// Helper function: Convert radians to Hours:Minutes:Seconds +function raToHMS(radians: number): string { + const hours = (radians * (12 / Math.PI)) % 24; + const h = Math.floor(hours); + const m = Math.floor((hours - h) * 60); + const s = ((hours - h) * 60 - m) * 60; + + return `${h.toString().padStart(2, '0')}h ${m.toString().padStart(2, '0')}m ${s.toFixed(2).padStart(5, '0')}s`; +} + +// Helper function: Convert radians to Degrees:Minutes:Seconds +function decToDMS(radians: number): string { + const degrees = radians * (180 / Math.PI); + const sign = degrees >= 0 ? '+' : '-'; + const absDegrees = Math.abs(degrees); + const d = Math.floor(absDegrees); + const m = Math.floor((absDegrees - d) * 60); + const s = ((absDegrees - d) * 60 - m) * 60; + + return `${sign}${d.toString().padStart(2, '0')}° ${m.toString().padStart(2, '0')}' ${s.toFixed(2).padStart(5, '0')}"`; +} diff --git a/examples/orbital-elements.ts b/examples/orbital-elements.ts new file mode 100644 index 00000000..82c8d6c7 --- /dev/null +++ b/examples/orbital-elements.ts @@ -0,0 +1,184 @@ +/* eslint-disable no-console */ +/** + * Example demonstrating orbital elements and conversions. + * + * This example shows: + * - Creating satellites from TLEs + * - Extracting classical orbital elements + * - Converting between state vectors and classical elements + * - Creating TLEs from classical elements + */ + +import { + ClassicalElements, + Degrees, + EpochUTC, + J2000, + Kilometers, + KilometersPerSecond, + Radians, + Satellite, + Tle, + TleLine1, + TleLine2, + Vector3D, +} from '../dist/main.js'; + +// Example 1: Extract orbital elements from a TLE +console.log('=== Example 1: Extract Orbital Elements from TLE ===\n'); + +const issTle = new Satellite({ + tle1: '1 25544U 98067A 24028.54545847 .00031576 00000-0 57240-3 0 9991' as TleLine1, + tle2: '2 25544 51.6418 292.2590 0002595 167.5319 252.0460 15.49326324436741' as TleLine2, +}); + +const date = new Date('2024-01-28T13:05:27.451Z'); +const j2000State = issTle.toJ2000(date); +const elements = j2000State.toClassicalElements(); + +console.log('ISS Orbital Elements:'); +console.log(` Semi-major axis: ${elements.semimajorAxis.toFixed(2)} km`); +console.log(` Eccentricity: ${elements.eccentricity.toFixed(6)}`); +console.log(` Inclination: ${(elements.inclination * (180 / Math.PI)).toFixed(4)}°`); +console.log(` Right Ascension: ${(elements.rightAscension * (180 / Math.PI)).toFixed(4)}°`); +console.log(` Arg of Perigee: ${(elements.argPerigee * (180 / Math.PI)).toFixed(4)}°`); +console.log(` True Anomaly: ${(elements.trueAnomaly * (180 / Math.PI)).toFixed(4)}°`); + +// Calculate derived parameters +console.log(`\nDerived Parameters:`); +console.log(` Period: ${elements.period.toFixed(2)} minutes`); +console.log(` Apogee altitude: ${((elements.semimajorAxis * (1 + elements.eccentricity)) - 6378.137).toFixed(2)} km`); +console.log(` Perigee altitude: ${((elements.semimajorAxis * (1 - elements.eccentricity)) - 6378.137).toFixed(2)} km`); + +// Example 2: Create classical elements and convert to state vector +console.log('\n=== Example 2: Create Orbital Elements and Convert to State Vector ===\n'); + +const customElements = new ClassicalElements({ + epoch: EpochUTC.fromDateTime(date), + semimajorAxis: 8000 as Kilometers, + eccentricity: 0.1, + inclination: 45 as Degrees, + rightAscension: 90 as Degrees, + argPerigee: 30 as Degrees, + trueAnomaly: 0 as Degrees, +}); + +console.log('Custom Orbit:'); +console.log(` Semi-major axis: ${customElements.semimajorAxis.toFixed(2)} km`); +console.log(` Eccentricity: ${customElements.eccentricity.toFixed(6)}`); +console.log(` Inclination: ${customElements.inclination.toFixed(4)}°`); +console.log(` Period: ${customElements.period.toFixed(2)} minutes`); + +const stateVector = customElements.toJ2000(); + +console.log(`\nState Vector:`); +console.log(` Position: [${stateVector.position.x.toFixed(2)}, ${stateVector.position.y.toFixed(2)}, ${stateVector.position.z.toFixed(2)}] km`); +console.log(` Velocity: [${stateVector.velocity.x.toFixed(6)}, ${stateVector.velocity.y.toFixed(6)}, ${stateVector.velocity.z.toFixed(6)}] km/s`); + +// Example 3: Convert state vector to classical elements +console.log('\n=== Example 3: Create State Vector and Convert to Classical Elements ===\n'); + +const customState = new J2000( + EpochUTC.fromDateTime(date), + new Vector3D( + -4040.9257 as Kilometers, + -4884.0906 as Kilometers, + 3522.9643 as Kilometers, + ), + new Vector3D( + 5.4662 as KilometersPerSecond, + -3.4425 as KilometersPerSecond, + -2.4854 as KilometersPerSecond, + ), +); + +const derivedElements = customState.toClassicalElements(); + +console.log('Derived Orbital Elements:'); +console.log(` Semi-major axis: ${derivedElements.semimajorAxis.toFixed(2)} km`); +console.log(` Eccentricity: ${derivedElements.eccentricity.toFixed(6)}`); +console.log(` Inclination: ${(derivedElements.inclination * (180 / Math.PI)).toFixed(4)}°`); +console.log(` Right Ascension: ${(derivedElements.rightAscension * (180 / Math.PI)).toFixed(4)}°`); +console.log(` Arg of Perigee: ${(derivedElements.argPerigee * (180 / Math.PI)).toFixed(4)}°`); +console.log(` True Anomaly: ${(derivedElements.trueAnomaly * (180 / Math.PI)).toFixed(4)}°`); + +// Example 4: Create TLE from classical elements +console.log('\n=== Example 4: Create TLE from Classical Elements ===\n'); + +const geoElements = new ClassicalElements({ + epoch: EpochUTC.fromDateTime(date), + semimajorAxis: 42164 as Kilometers, // GEO altitude + eccentricity: 0.0001, + inclination: 0.1 as Degrees, + rightAscension: 0 as Degrees, + argPerigee: 0 as Degrees, + trueAnomaly: 0 as Degrees, +}); + +const geoTle = Tle.fromClassicalElements(geoElements); + +console.log('Generated TLE for GEO satellite:'); +console.log(geoTle.line1); +console.log(geoTle.line2); + +// Verify by reading back +const verifyTle = new Satellite({ + tle1: geoTle.line1, + tle2: geoTle.line2, +}); + +const verifyState = verifyTle.toJ2000(date); +const verifyElements = verifyState.toClassicalElements(); + +console.log(`\nVerification - Semi-major axis: ${verifyElements.semimajorAxis.toFixed(2)} km`); +console.log(`Verification - Period: ${verifyElements.period.toFixed(2)} minutes (should be ~1436 min for GEO)`); + +// Example 5: Different orbit types +console.log('\n=== Example 5: Different Orbit Types ===\n'); + +const orbits = [ + { + name: 'LEO (Circular)', + elements: new ClassicalElements({ + epoch: EpochUTC.fromDateTime(date), + semimajorAxis: 6778 as Kilometers, + eccentricity: 0.001, + inclination: 51.6 as Degrees, + rightAscension: 0 as Degrees, + argPerigee: 0 as Degrees, + trueAnomaly: 0 as Degrees, + }), + }, + { + name: 'MEO (GPS-like)', + elements: new ClassicalElements({ + epoch: EpochUTC.fromDateTime(date), + semimajorAxis: 26560 as Kilometers, + eccentricity: 0.01, + inclination: 55 as Degrees, + rightAscension: 0 as Degrees, + argPerigee: 0 as Degrees, + trueAnomaly: 0 as Degrees, + }), + }, + { + name: 'HEO (Molniya)', + elements: new ClassicalElements({ + epoch: EpochUTC.fromDateTime(date), + semimajorAxis: 26554 as Kilometers, + eccentricity: 0.74, + inclination: 63.4 as Degrees, + rightAscension: 0 as Degrees, + argPerigee: 270 as Degrees, + trueAnomaly: 0 as Degrees, + }), + }, +]; + +orbits.forEach((orbit) => { + console.log(`${orbit.name}:`); + console.log(` Altitude at apogee: ${((orbit.elements.semimajorAxis * (1 + orbit.elements.eccentricity)) - 6378.137).toFixed(2)} km`); + console.log(` Altitude at perigee: ${((orbit.elements.semimajorAxis * (1 - orbit.elements.eccentricity)) - 6378.137).toFixed(2)} km`); + console.log(` Period: ${orbit.elements.period.toFixed(2)} minutes`); + console.log(''); +}); diff --git a/examples/satellite-passes.ts b/examples/satellite-passes.ts new file mode 100644 index 00000000..6f54ac3c --- /dev/null +++ b/examples/satellite-passes.ts @@ -0,0 +1,234 @@ +/* eslint-disable no-console */ +/** + * Example demonstrating satellite pass prediction and visibility calculations. + * + * This example shows: + * - Calculating satellite passes + * - Finding when satellites are visible from a ground station + * - Computing look angles (azimuth, elevation, range) + * - Checking field of view constraints + */ + +import { + Degrees, + DetailedSensor, + Kilometers, + Satellite, + Sensor, + SpaceObjectType, + TleLine1, + TleLine2, +} from '../dist/main.js'; + +// Example 1: Calculate passes for ISS +console.log('=== Example 1: ISS Pass Prediction ===\n'); + +// Create a ground station +const groundStation = new Sensor({ + lat: 41.754785 as Degrees, + lon: -70.539151 as Degrees, + alt: 0.060966 as Kilometers, + name: 'Cape Cod', +}); + +// Create ISS satellite +const iss = new Satellite({ + tle1: '1 25544U 98067A 24028.54545847 .00031576 00000-0 57240-3 0 9991' as TleLine1, + tle2: '2 25544 51.6418 292.2590 0002595 167.5319 252.0460 15.49326324436741' as TleLine2, +}); + +console.log(`Ground Station: ${groundStation.name}`); +console.log(` Location: ${groundStation.lat}° N, ${Math.abs(groundStation.lon)}° W`); +console.log(` Altitude: ${groundStation.alt} km`); + +console.log('\nSatellite: ISS'); +console.log(` Inclination: ${iss.inclination}°`); +console.log(` Period: ${iss.period} minutes`); + +// Calculate passes (checking every 30 seconds for next 24 hours) +const passes = groundStation.calculatePasses(30, iss); + +console.log(`\nFound ${passes.length} passes in the next 24 hours:\n`); + +passes.slice(0, 5).forEach((pass, index) => { + console.log(`Pass ${index + 1}:`); + console.log(` Rise Time: ${pass.start.toLocaleString()}`); + console.log(` Set Time: ${pass.end.toLocaleString()}`); + console.log(` Duration: ${((pass.end.getTime() - pass.start.getTime()) / 1000 / 60).toFixed(1)} minutes`); + console.log(` Max Elevation: ${pass.maxEl.toFixed(1)}°`); + + if (pass.maxEl > 45) { + console.log(` ⭐ Excellent pass!`); + } else if (pass.maxEl > 20) { + console.log(` ✓ Good pass`); + } + + console.log(''); +}); + +// Example 2: Check visibility at a specific time +console.log('=== Example 2: Satellite Visibility Check ===\n'); + +const checkTime = new Date('2024-01-28T12:00:00.000Z'); + +// Get current look angles +const rae = groundStation.rae(iss, checkTime); + +console.log(`Time: ${checkTime.toISOString()}`); +console.log(`\nLook Angles:`); +console.log(` Azimuth: ${rae.az.toFixed(2)}°`); +console.log(` Elevation: ${rae.el.toFixed(2)}°`); +console.log(` Range: ${rae.rng.toFixed(2)} km`); + +// Check if satellite is visible +const isVisible = rae.el > 0; + +console.log(`\nSatellite is ${isVisible ? 'VISIBLE' : 'BELOW HORIZON'}`); + +if (isVisible) { + console.log(`\nDirection: ${getCardinalDirection(rae.az)}`); + console.log(`Elevation: ${getElevationDescription(rae.el)}`); +} + +// Example 3: Detailed sensor with field of view constraints +console.log('\n=== Example 3: Field of View Constraints ===\n'); + +const radar = new DetailedSensor({ + lat: 41.754785 as Degrees, + lon: -70.539151 as Degrees, + alt: 0.060966 as Kilometers, + minAz: 0 as Degrees, + maxAz: 360 as Degrees, + minEl: 10 as Degrees, // Minimum elevation 10° + maxEl: 85 as Degrees, + minRng: 100 as Kilometers, // Minimum range 100 km + maxRng: 5556 as Kilometers, // Maximum range 5556 km + name: 'Cape Cod Radar', + type: SpaceObjectType.PHASED_ARRAY_RADAR, +}); + +console.log(`Sensor: ${radar.name}`); +console.log(`Field of View Constraints:`); +console.log(` Azimuth: ${radar.minAz}° - ${radar.maxAz}°`); +console.log(` Elevation: ${radar.minEl}° - ${radar.maxEl}°`); +console.log(` Range: ${radar.minRng} - ${radar.maxRng} km`); + +// Check if satellite is in FOV +const inFov = radar.isSatInFov(iss, checkTime); + +console.log(`\nAt ${checkTime.toISOString()}:`); +console.log(` Satellite in FOV: ${inFov ? 'YES ✓' : 'NO ✗'}`); + +if (inFov) { + console.log(` The satellite meets all FOV constraints`); +} else { + const raeCheck = radar.rae(iss, checkTime); + + console.log(` Constraints not met:`); + + if (raeCheck.el < radar.minEl) { + console.log(` - Elevation too low (${raeCheck.el.toFixed(1)}° < ${radar.minEl}°)`); + } + + if (raeCheck.el > radar.maxEl) { + console.log(` - Elevation too high (${raeCheck.el.toFixed(1)}° > ${radar.maxEl}°)`); + } + + if (raeCheck.rng < radar.minRng) { + console.log(` - Range too close (${raeCheck.rng.toFixed(0)} km < ${radar.minRng} km)`); + } + + if (raeCheck.rng > radar.maxRng) { + console.log(` - Range too far (${raeCheck.rng.toFixed(0)} km > ${radar.maxRng} km)`); + } +} + +// Example 4: Track satellite across the sky +console.log('\n=== Example 4: Satellite Tracking Over Time ===\n'); + +const trackStart = new Date('2024-01-28T12:00:00.000Z'); + +console.log('ISS Position every 5 minutes:\n'); +console.log('Time Az El Range Status'); +console.log('───────────────────── ────── ────── ─────── ────────'); + +for (let i = 0; i < 12; i++) { + const trackTime = new Date(trackStart.getTime() + i * 5 * 60 * 1000); + const trackRae = groundStation.rae(iss, trackTime); + + const timeStr = trackTime.toISOString().substring(11, 19); + const azStr = trackRae.az.toFixed(1).padStart(6); + const elStr = trackRae.el.toFixed(1).padStart(6); + const rngStr = trackRae.rng.toFixed(0).padStart(7); + + let status = 'Below horizon'; + + if (trackRae.el > 0) { + status = 'Visible'; + + if (radar.isSatInFov(iss, trackTime)) { + status = 'In FOV ✓'; + } + } + + console.log(`${timeStr} ${azStr}° ${elStr}° ${rngStr} km ${status}`); +} + +// Example 5: Multiple satellites +console.log('\n=== Example 5: Tracking Multiple Satellites ===\n'); + +const satellites = [ + { + name: 'ISS', + sat: new Satellite({ + tle1: '1 25544U 98067A 24028.54545847 .00031576 00000-0 57240-3 0 9991' as TleLine1, + tle2: '2 25544 51.6418 292.2590 0002595 167.5319 252.0460 15.49326324436741' as TleLine2, + }), + }, + { + name: 'HST (Hubble)', + sat: new Satellite({ + tle1: '1 20580U 90037B 24028.50123227 .00000825 00000-0 39644-4 0 9997' as TleLine1, + tle2: '2 20580 28.4696 273.2640 0002975 297.7865 189.2151 15.09696656316758' as TleLine2, + }), + }, +]; + +const multiCheckTime = new Date('2024-01-28T18:00:00.000Z'); + +console.log(`Time: ${multiCheckTime.toISOString()}\n`); +console.log('Satellite Az El Range Visible In FOV'); +console.log('──────────── ────── ────── ─────── ──────── ──────'); + +satellites.forEach((satInfo) => { + const satRae = groundStation.rae(satInfo.sat, multiCheckTime); + const satVisible = satRae.el > 0; + const satInFov = radar.isSatInFov(satInfo.sat, multiCheckTime); + + const nameStr = satInfo.name.padEnd(12); + const azStr = satRae.az.toFixed(1).padStart(6); + const elStr = satRae.el.toFixed(1).padStart(6); + const rngStr = satRae.rng.toFixed(0).padStart(7); + const visStr = (satVisible ? 'Yes' : 'No').padEnd(8); + const fovStr = satInFov ? 'Yes ✓' : 'No'; + + console.log(`${nameStr} ${azStr}° ${elStr}° ${rngStr} km ${visStr} ${fovStr}`); +}); + +// Helper functions +function getCardinalDirection(azimuth: number): string { + const directions = ['N', 'NNE', 'NE', 'ENE', 'E', 'ESE', 'SE', 'SSE', 'S', 'SSW', 'SW', 'WSW', 'W', 'WNW', 'NW', 'NNW']; + const index = Math.round(azimuth / 22.5) % 16; + + return directions[index]; +} + +function getElevationDescription(elevation: number): string { + if (elevation > 80) return 'Nearly overhead'; + if (elevation > 60) return 'Very high in the sky'; + if (elevation > 45) return 'High in the sky'; + if (elevation > 30) return 'Medium elevation'; + if (elevation > 15) return 'Low in the sky'; + + return 'Near the horizon'; +} diff --git a/examples/time-systems.ts b/examples/time-systems.ts new file mode 100644 index 00000000..6378d578 --- /dev/null +++ b/examples/time-systems.ts @@ -0,0 +1,199 @@ +/* eslint-disable no-console */ +/** + * Example demonstrating different time systems and epoch conversions. + * + * This example shows: + * - Different epoch types (UTC, TAI, TT, TDB, GPS) + * - Converting between time systems + * - Julian dates + * - GMST (Greenwich Mean Sidereal Time) + * - Time differences and offsets + */ + +import { + calcGmst, + Days, + EpochGPS, + EpochTAI, + EpochTDB, + EpochTT, + EpochUTC, + jday, +} from '../dist/main.js'; + +// Example 1: Create epochs in different time systems +console.log('=== Example 1: Different Time Systems ===\n'); + +const date = new Date('2024-01-28T12:00:00.000Z'); + +const utcEpoch = EpochUTC.fromDateTime(date); +const taiEpoch = EpochTAI.fromDateTime(date); +const ttEpoch = EpochTT.fromDateTime(date); +const tdbEpoch = EpochTDB.fromDateTime(date); +const gpsEpoch = EpochGPS.fromDateTime(date); + +console.log(`UTC (Coordinated Universal Time): ${utcEpoch.toDateTime().toISOString()}`); +console.log(`TAI (International Atomic Time): ${taiEpoch.toDateTime().toISOString()}`); +console.log(`TT (Terrestrial Time): ${ttEpoch.toDateTime().toISOString()}`); +console.log(`TDB (Barycentric Dynamical Time): ${tdbEpoch.toDateTime().toISOString()}`); +console.log(`GPS (GPS Time): ${gpsEpoch.toDateTime().toISOString()}`); + +// Example 2: Time system conversions +console.log('\n=== Example 2: Converting Between Time Systems ===\n'); + +console.log('Starting with UTC:'); +console.log(` UTC: ${utcEpoch.toDateTime().toISOString()}`); + +// Convert UTC to other systems +const utcToTai = utcEpoch.toTai(); +const utcToTt = utcEpoch.toTt(); +const utcToTdb = utcEpoch.toTdb(); +const utcToGps = utcEpoch.toGps(); + +console.log(` → TAI: ${utcToTai.toDateTime().toISOString()}`); +console.log(` → TT: ${utcToTt.toDateTime().toISOString()}`); +console.log(` → TDB: ${utcToTdb.toDateTime().toISOString()}`); +console.log(` → GPS: ${utcToGps.toDateTime().toISOString()}`); + +console.log('\nStarting with TAI:'); +console.log(` TAI: ${taiEpoch.toDateTime().toISOString()}`); + +const taiToUtc = taiEpoch.toUtc(); +const taiToTt = taiEpoch.toTt(); + +console.log(` → UTC: ${taiToUtc.toDateTime().toISOString()}`); +console.log(` → TT: ${taiToTt.toDateTime().toISOString()}`); + +// Example 3: Julian dates +console.log('\n=== Example 3: Julian Dates ===\n'); + +const jd = jday( + date.getUTCFullYear(), + date.getUTCMonth() + 1, + date.getUTCDate(), + date.getUTCHours(), + date.getUTCMinutes(), + date.getUTCSeconds(), +); + +console.log(`Date: ${date.toISOString()}`); +console.log(`Julian Date: ${jd.toFixed(6)}`); +console.log(`Modified Julian Date: ${(jd - 2400000.5).toFixed(6)}`); + +// J2000 epoch (January 1, 2000, 12:00:00 TT) +const j2000Date = new Date('2000-01-01T12:00:00.000Z'); +const j2000Jd = jday( + j2000Date.getUTCFullYear(), + j2000Date.getUTCMonth() + 1, + j2000Date.getUTCDate(), + j2000Date.getUTCHours(), + j2000Date.getUTCMinutes(), + j2000Date.getUTCSeconds(), +); + +console.log(`\nJ2000 Epoch:`); +console.log(` Date: ${j2000Date.toISOString()}`); +console.log(` Julian Date: ${j2000Jd.toFixed(6)}`); + +// Example 4: GMST (Greenwich Mean Sidereal Time) +console.log('\n=== Example 4: Greenwich Mean Sidereal Time ===\n'); + +const gmstResult = calcGmst(date); + +console.log(`Date: ${date.toISOString()}`); +console.log(`GMST: ${gmstResult.gmst.toFixed(6)} radians`); +console.log(` ${(gmstResult.gmst * (180 / Math.PI)).toFixed(6)}°`); +console.log(` ${((gmstResult.gmst * (180 / Math.PI)) / 15).toFixed(6)} hours`); + +// Convert to hour:minute:second format +const gmstHours = (gmstResult.gmst * (180 / Math.PI)) / 15; +const hours = Math.floor(gmstHours); +const minutes = Math.floor((gmstHours - hours) * 60); +const seconds = ((gmstHours - hours) * 60 - minutes) * 60; + +console.log(` ${hours}h ${minutes}m ${seconds.toFixed(2)}s`); + +// Example 5: Time differences +console.log('\n=== Example 5: Time System Offsets ===\n'); + +// TAI-UTC offset (leap seconds) +const taiUtcDiff = (utcToTai.toDateTime().getTime() - utcEpoch.toDateTime().getTime()) / 1000; + +console.log(`TAI - UTC = ${taiUtcDiff} seconds (leap seconds)`); + +// TT-TAI offset (constant 32.184 seconds) +const ttTaiDiff = (utcToTt.toDateTime().getTime() - utcToTai.toDateTime().getTime()) / 1000; + +console.log(`TT - TAI = ${ttTaiDiff.toFixed(3)} seconds (constant)`); + +// GPS-UTC offset +const gpsUtcDiff = (utcToGps.toDateTime().getTime() - utcEpoch.toDateTime().getTime()) / 1000; + +console.log(`GPS - UTC = ${gpsUtcDiff} seconds`); + +// Example 6: Working with epoch arithmetic +console.log('\n=== Example 6: Epoch Arithmetic ===\n'); + +const startEpoch = EpochUTC.fromDateTime(new Date('2024-01-01T00:00:00.000Z')); + +console.log(`Start: ${startEpoch.toDateTime().toISOString()}`); + +// Add 1 day +const oneDayLater = startEpoch.roll(1 as Days); + +console.log(`+1 day: ${oneDayLater.toDateTime().toISOString()}`); + +// Add 7 days +const oneWeekLater = startEpoch.roll(7 as Days); + +console.log(`+7 days: ${oneWeekLater.toDateTime().toISOString()}`); + +// Example 7: Multiple date conversions +console.log('\n=== Example 7: Common Date/Time Scenarios ===\n'); + +const scenarios = [ + { name: 'GPS Epoch Start', date: new Date('1980-01-06T00:00:00.000Z') }, + { name: 'J2000.0', date: new Date('2000-01-01T12:00:00.000Z') }, + { name: 'Unix Epoch', date: new Date('1970-01-01T00:00:00.000Z') }, + { name: 'Current Example', date: new Date('2024-01-28T12:00:00.000Z') }, +]; + +scenarios.forEach((scenario) => { + const jdValue = jday( + scenario.date.getUTCFullYear(), + scenario.date.getUTCMonth() + 1, + scenario.date.getUTCDate(), + scenario.date.getUTCHours(), + scenario.date.getUTCMinutes(), + scenario.date.getUTCSeconds(), + ); + + const gmstValue = calcGmst(scenario.date); + + console.log(`${scenario.name}:`); + console.log(` Date: ${scenario.date.toISOString()}`); + console.log(` JD: ${jdValue.toFixed(6)}`); + console.log(` GMST: ${(gmstValue.gmst * (180 / Math.PI) / 15).toFixed(4)} hours`); + console.log(''); +}); + +// Example 8: Precision timing +console.log('=== Example 8: High-Precision Time Comparison ===\n'); + +const preciseDate = new Date('2024-01-28T12:34:56.789Z'); +const utcPrecise = EpochUTC.fromDateTime(preciseDate); + +console.log(`Input: ${preciseDate.toISOString()}`); +console.log(`UTC epoch posix: ${utcPrecise.posix.toFixed(9)} seconds since Unix epoch`); + +// Show difference between time systems in milliseconds +const taiPrecise = utcPrecise.toTai(); +const ttPrecise = utcPrecise.toTt(); + +const utcMillis = utcPrecise.toDateTime().getTime(); +const taiMillis = taiPrecise.toDateTime().getTime(); +const ttMillis = ttPrecise.toDateTime().getTime(); + +console.log(`\nTime differences (from UTC):`); +console.log(` TAI: +${(taiMillis - utcMillis).toFixed(0)} ms`); +console.log(` TT: +${(ttMillis - utcMillis).toFixed(0)} ms`);