Initial Orbit Determination

This example recovers an orbit from a handful of position fixes using three classic IOD methods: Lambert (two positions), Gibbs (three positions), and Herrick-Gibbs (three closely spaced positions). Each solution is then fitted to a TLE so the results can be compared side by side and fed to SGP4 consumers.

import {
  DEG2RAD, Degrees, EpochUTC, J2000, Kilometers, KilometersPerSecond, RAE, Radians, Tle, Vector3D, calcGmst,
  lla2eci, LambertIOD, GibbsIOD, HerrickGibbsIOD,
} from 'ootk';

Run it

npm run build
npx tsx ./examples/iod.ts

Setup Observations

Three radar-style range/azimuth/elevation (RAE) measurements taken 10 seconds apart from a fixed ground sensor, plus the sensor's geodetic location (latitude and longitude in radians, altitude in km).

const lambert = new LambertIOD();

const rae1 = {
  t: EpochUTC.fromDateTime(new Date(1704628462000)),
  rng: 1599.89 as Kilometers,
  az: 174 as Degrees,
  el: 13.6 as Degrees,
};
const rae2 = {
  t: EpochUTC.fromDateTime(new Date(1704628462000 + 10 * 1000)),
  rng: 1568.76 as Kilometers,
  az: 171 as Degrees,
  el: 14.2 as Degrees,
};
const rae3 = {
  t: EpochUTC.fromDateTime(new Date(1704628462000 + 20 * 1000)),
  rng: 1540.09 as Kilometers,
  az: 169 as Degrees,
  el: 14.7 as Degrees,
};
const sensor = {
  lat: (41.754785 * DEG2RAD) as Radians,
  lon: (-70.539151 * DEG2RAD) as Radians,
  alt: 0.085 as Kilometers,
};

RAE to Position

IOD methods work on inertial positions, so each RAE measurement is converted: calcGmst and lla2eci produce the sensor's ECI position, and RAE.fromDegrees(...).toStateVector(site) turns the measurement into a J2000 state relative to that site. The site velocity is left at zero here, which is acceptable because only the positions are used downstream.

const gmst = calcGmst(rae1.t.toDateTime());
const sensorEci = lla2eci(sensor, gmst.gmst);

const p1 = RAE.fromDegrees(rae1.t, rae1.rng, rae1.az, rae1.el).toStateVector(
  new J2000(rae1.t,
    new Vector3D(sensorEci.x, sensorEci.y, sensorEci.z), Vector3D.origin as Vector3D<KilometersPerSecond>,
  ),
);
const p2 = RAE.fromDegrees(rae2.t, rae2.rng, rae2.az, rae2.el).toStateVector(
  new J2000(rae2.t,
    new Vector3D(sensorEci.x, sensorEci.y, sensorEci.z), Vector3D.origin as Vector3D<KilometersPerSecond>,
  ),
);
const p3 = RAE.fromDegrees(rae3.t, rae3.rng, rae3.az, rae3.el).toStateVector(
  new J2000(rae3.t,
    new Vector3D(sensorEci.x, sensorEci.y, sensorEci.z), Vector3D.origin as Vector3D<KilometersPerSecond>,
  ),
);

ECI Positions

A second data set of known ECI positions for the same object (recorded from a propagator) at 10 and 30 second spacings. These bypass the sensor geometry entirely and represent the cleanest possible IOD input.

const eci1 = {
  x: -4901.84521484375 as Kilometers,
  y: -3592.527587890625 as Kilometers,
  z: 3322.875732421875 as Kilometers,
};
const eci2 = {
  x: -4847.90185546875 as Kilometers,
  y: -3631.424560546875 as Kilometers,
  z: 3359.44482421875 as Kilometers,
};
const eci3 = {
  x: -4793.376953125 as Kilometers,
  y: -3669.885986328125 as Kilometers,
  z: 3395.60986328125 as Kilometers,
};
const p1b = new J2000(rae1.t, new Vector3D(eci1.x, eci1.y, eci1.z), Vector3D.origin as Vector3D<KilometersPerSecond>);
const p2b = new J2000(rae2.t, new Vector3D(eci2.x, eci2.y, eci2.z), Vector3D.origin as Vector3D<KilometersPerSecond>);
const p3b = new J2000(rae3.t, new Vector3D(eci3.x, eci3.y, eci3.z), Vector3D.origin as Vector3D<KilometersPerSecond>);

const eci4 = { x: -4738.27734375 as Kilometers, y: -3707.9072265625 as Kilometers, z: 3431.36669921875 as Kilometers };
const p4b = new J2000(
  EpochUTC.fromDateTime(new Date(1704628492000)),
  new Vector3D(eci4.x, eci4.y, eci4.z),
  Vector3D.origin as Vector3D<KilometersPerSecond>,
);

const eci5 = {
  x: -4569.595703125 as Kilometers,
  y: -3819.285400390625 as Kilometers,
  z: 3536.144287109375 as Kilometers,
};
const p5b = new J2000(
  EpochUTC.fromDateTime(new Date(1704628522000)),
  new Vector3D(eci5.x, eci5.y, eci5.z),
  Vector3D.origin as Vector3D<KilometersPerSecond>,
);

Solve IOD

Four estimates are produced: LambertIOD.estimate from two positions and their epochs, HerrickGibbsIOD.solve from the three RAE-derived positions, GibbsIOD.solve from the three ECI positions, and HerrickGibbsIOD.solve again from ECI positions 30 seconds apart. Each returns a J2000 state (position and velocity) at one of the observation epochs. Herrick-Gibbs is preferred over Gibbs when the positions are closely spaced (small angular separation).

const estimate = lambert.estimate(p1.position, p2.position, p1.epoch, p2.epoch);
const estimate2 = new HerrickGibbsIOD().solve(p1.position, p1.epoch, p2.position, p2.epoch, p3.position, p3.epoch);
const estimate3 = new GibbsIOD().solve(p1b.position, p2b.position, p3b.position, p2b.epoch, p3b.epoch);
const estimate4 = new HerrickGibbsIOD().solve(
  p1b.position,
  p1b.epoch,
  p4b.position,
  p4b.epoch,
  p5b.position,
  p5b.epoch,
);

Fit TLEs

toClassicalElements() converts each solved state to Keplerian elements, and Tle.fromClassicalElements() fits a TLE. Note how the two ECI-based solutions (Gibbs and Herrick-Gibbs) agree closely, while the Lambert solution from only two noisy RAE fixes is far less constrained (its eccentricity is clearly wrong), showing how sensitive IOD is to input quality and observation count.

const tleLambert = Tle.fromClassicalElements(estimate.toClassicalElements());
const tleHerrickGibbsRae = Tle.fromClassicalElements(estimate2.toClassicalElements());
const tle = Tle.fromClassicalElements(estimate3.toClassicalElements());
const tle2 = Tle.fromClassicalElements(estimate4.toClassicalElements());

console.log('Lambert IOD (two RAE-derived positions, 10 s apart) fitted to a TLE:');
console.log(`  ${tleLambert.line1}`);
console.log(`  ${tleLambert.line2}`);

console.log('\nHerrick-Gibbs IOD (three RAE-derived positions, 10 s spacing) fitted to a TLE:');
console.log(`  ${tleHerrickGibbsRae.line1}`);
console.log(`  ${tleHerrickGibbsRae.line2}`);

console.log('\nGibbs IOD (three ECI positions, 10 s spacing) fitted to a TLE:');
console.log(`  ${tle.line1}`);
console.log(`  ${tle.line2}`);

console.log('\nHerrick-Gibbs IOD (three ECI positions, 30 s spacing) fitted to a TLE:');
console.log(`  ${tle2.line1}`);
console.log(`  ${tle2.line2}`);

Output

Lambert IOD (two RAE-derived positions, 10 s apart) fitted to a TLE:
  1 00001U 58001A   24007.49608796  .00000000  00000+0  00000+0 0  9995
  2 00001  40.7526 176.3707 3432410  46.3020   0.6275 08.01238333    01

Herrick-Gibbs IOD (three RAE-derived positions, 10 s spacing) fitted to a TLE:
  1 00001U 58001A   24007.49620370  .00000000  00000+0  00000+0 0  9997
  2 00001  43.5386 180.7830 0388414 224.2656 180.9250 15.93538194    07

Gibbs IOD (three ECI positions, 10 s spacing) fitted to a TLE:
  1 00001U 58001A   24007.49620370  .00000000  00000+0  00000+0 0  9997
  2 00001  42.9994 180.3394 0006341 105.1309 300.2630 15.05329933    05

Herrick-Gibbs IOD (three ECI positions, 30 s spacing) fitted to a TLE:
  1 00001U 58001A   24007.49643519  .00000000  00000+0  00000+0 0  9997
  2 00001  42.9985 180.3381 0010714  76.6282 330.0209 15.03982493    09