Move the solar system on JPL's mean elements, so it stays right as the clock runs

Every body carried one set of osculating elements from Horizons at 2025-01-01, run forward by Kepler with a GM from a table of mass ratios. That set is exact at its instant and drifts from then on, and the clock now runs a month a second: the Moon, with Earth's mass ratio lacking its own and the osculating axis, went round in 27.70 days instead of 27.32, 66 degrees out after a year, and its locked face was spun at the same wrong rate. Planets and Pluto now take Standish's Table 2a/2b ("Keplerian Elements for Approximate Positions of the Major Planets"): elements against the J2000 ecliptic, their rates per century, and the b, c, s, f terms of Jupiter to Pluto, fit for 3000 BC to AD 3000. Table 1 is closer near the present (Saturn 0.23 degrees at worst 1950-2100, against 0.32 here) but is only fit for 1800-2050, and by AD 3000 has Saturn 4.3 degrees out where Table 2 holds every planet within 0.3. The moons take JPL SSD's satellite mean elements: sidereal mean motion n to ten figures, the periods of their node and periapsis, and each one's local Laplace plane by its pole. They propagate with n itself, never a GM: gmForParent and its mass table are gone. Horizons still gives size, spin and obliquity. Both tables are read from the Internet Archive's copy of JPL's pages, pinned to one capture: the live approx_pos page has dropped Pluto, and the live sats/elem page has dropped n and rounds the period to four or five figures (Phobos 0.3187 d, a revolution out within a decade). What the tables leave implicit, measured against Horizons before it was accepted: - The precession periods are magnitudes. A node regresses on a prograde orbit and advances on a retrograde one; a periapsis advances except where a resonance forces the eccentricity. Io's and Europa's follow their conjunction line backwards at 2 n(Europa) - n(Io) = 0.74 degrees a day, which is exactly the 1.625- and 1.394-year periods in the table. Read as advancing, Io was 0.9 degrees out and Europa 2.1. - On a retrograde orbit the node's turning is added back to the mean anomaly. Taken off, Triton drifted a degree a year, 105 degrees by 2100. - The Laplace frame's x axis is where the plane rises through the ICRF equator, RA of the pole plus 90. Read against the ecliptic, Io was 2.8 degrees out, Phobos 54 and Titan 127. Orbit lines are now drawn in their own plane and turned by a quaternion each tick, so a turning node carries the line with the body: fixed at one date, the Moon's line would be up to 69 000 km off it nine years on. The Earth row is the Earth-Moon barycentre, 4 700 km from Earth, 0.002 degrees from the Sun. A tidally locked moon's day is now 360 / n, its sidereal period (the Moon 27.321662 d), so it stays locked to the orbit it is drawn on. Angular error against Horizons VECTORS (ICRF, TDB; heliocentric for planets, planet-centred for moons), degrees, read from the live renderer's markers in the running app: body 1950-01-01 1975-01-01 1987-07-23 2000-01-01 2025-01-01 2037-03-06 2050-01-01 2075-01-01 2100-01-01 max mercury 0.004 0.002 0.003 0.002 0.002 0.001 0.000 0.002 0.000 0.004 venus 0.003 0.007 0.003 0.004 0.004 0.004 0.003 0.004 0.004 0.007 earth 0.003 0.008 0.002 0.005 0.004 0.009 0.003 0.002 0.003 0.009 mars 0.009 0.010 0.008 0.024 0.009 0.012 0.009 0.011 0.028 0.028 jupiter 0.063 0.030 0.171 0.135 0.013 0.020 0.056 0.041 0.075 0.171 saturn 0.080 0.064 0.018 0.320 0.066 0.114 0.044 0.164 0.177 0.320 uranus 0.018 0.169 0.068 0.050 0.101 0.015 0.141 0.017 0.114 0.169 neptune 0.070 0.028 0.004 0.021 0.036 0.037 0.013 0.029 0.072 0.072 pluto 0.045 0.054 0.041 0.033 0.019 0.020 0.023 0.027 0.026 0.054 moon 0.486 1.928 0.127 0.631 1.407 1.086 0.720 0.339 1.180 1.928 phobos 2.068 0.294 0.881 1.113 0.313 0.636 2.089 5.862 11.099 11.099 deimos 0.077 0.043 0.310 0.066 0.164 0.068 0.034 0.468 0.044 0.468 io 0.021 0.015 0.010 0.019 0.009 0.035 0.006 0.011 0.022 0.035 europa 0.036 0.039 0.053 0.064 0.078 0.032 0.006 0.034 0.044 0.078 ganymede 0.132 0.103 0.018 0.007 0.023 0.054 0.091 0.118 0.044 0.132 callisto 0.040 0.019 0.023 0.019 0.038 0.008 0.060 0.119 0.056 0.119 titan 0.003 0.019 0.023 0.023 0.027 0.028 0.048 0.008 0.014 0.048 triton 0.051 0.029 0.009 0.021 0.052 0.048 0.063 0.089 0.137 0.137 Three miss what was hoped for, and why: - Jupiter 0.17, Saturn 0.32, Uranus 0.17 against the 0.1 hoped for: short-period perturbations of the giants by one another, which no Keplerian fit carries. Standish states his own Table 2 errors as 600, 1 000 and 2 000 arcseconds (0.17, 0.28, 0.56 degrees). Out to AD 3000, measured at 1800, 2200, 2400, 2600 and 3000, every planet stays within 0.3. - The Moon, 1.9: evection (1.27) and variation (0.66), which a mean ellipse leaves out. - Phobos, 2.1 until 2050, then 5.9 in 2075 and 11.1 in 2100, growing as the square of the time: its tidal acceleration, which the table has no column for. Its elements are MAR080's, epoch 1950. The map's dates are also UTC where the elements are TDB, 69 s today, which is 0.9 degrees of Phobos and nothing for anything else. Held in place by: - build.ts: each body's mean elements against Horizons' own osculating elements on the ETL's 2025-01-01, at most 0.25 degrees for a planet and 2.5 for a moon (measured: Uranus 0.101, the Moon 1.407; a regressing Triton node reads 10.24 and fails), and every moon's day equal to its sidereal period (a 1% error fails). - Unit tests freezing nine Horizons vectors (Earth 2100, Jupiter 1950, Saturn 2075, Pluto 1975, the Moon 2050, Io and Europa 1950, Titan and Triton 2100) through SystemOrbitsRenderer, the Moon kept on its own turning line, the retrograde rule, the Standish terms, the Laplace frame, and both table parsers. Nine mutants each fail the test named for them, and the two validators each refuse a mutated build of the real catalogue. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
2026-09-24 20:18:57 +02:00
parent c38a42cbcb
commit 48319c3fe2
19 changed files with 1062 additions and 244 deletions
@@ -5,12 +5,13 @@ import { ExoplanetRecord } from '../models/exoplanet.model';
 import { appearanceForBody, appearanceForExoplanet, heliocentricDistanceAu } from './body-appearance';
 import { DEFAULT_EPOCH_JD } from './constants';

+const RATES = { meanMotionDegPerDay: 1, longitudeOfAscendingNodeDegPerDay: 0, argumentOfPeriapsisDegPerDay: 0 };
 const ORBIT = { eccentricity: 0, inclinationDeg: 0, longitudeOfAscendingNodeDeg: 0, argumentOfPeriapsisDeg: 0, meanAnomalyAtEpochDeg: 0, epochJd: DEFAULT_EPOCH_JD };

-const JUPITER: BodyRecord = { id: 'jupiter', systemStarId: 0, name: 'Jupiter', kind: 'planet', radiusKm: 69911, orbit: { ...ORBIT, semiMajorAxisAu: 5.204 } };
+const JUPITER: BodyRecord = { id: 'jupiter', systemStarId: 0, name: 'Jupiter', kind: 'planet', radiusKm: 69911, orbit: { ...ORBIT, semiMajorAxisAu: 5.204 }, rates: RATES, orbitSource: 'test' };
 /** Europa's own orbit is around Jupiter: 671,000 km, which is 0.00449 AU. */
-const EUROPA: BodyRecord = { id: 'europa', systemStarId: 0, name: 'Europa', kind: 'moon', radiusKm: 1560, parentBodyId: 'jupiter', orbit: { ...ORBIT, semiMajorAxisAu: 0.00449 } };
-const EARTH: BodyRecord = { id: 'earth', systemStarId: 0, name: 'Earth', kind: 'planet', radiusKm: 6371, orbit: { ...ORBIT, semiMajorAxisAu: 1 } };
+const EUROPA: BodyRecord = { id: 'europa', systemStarId: 0, name: 'Europa', kind: 'moon', radiusKm: 1560, parentBodyId: 'jupiter', orbit: { ...ORBIT, semiMajorAxisAu: 0.00449 }, rates: RATES, orbitSource: 'test' };
+const EARTH: BodyRecord = { id: 'earth', systemStarId: 0, name: 'Earth', kind: 'planet', radiusKm: 6371, orbit: { ...ORBIT, semiMajorAxisAu: 1 }, rates: RATES, orbitSource: 'test' };
 const ORPHAN: BodyRecord = { ...EUROPA, id: 'orphan', parentBodyId: 'nowhere' };

 const BODIES = [JUPITER, EUROPA, EARTH, ORPHAN];
@@ -19,33 +19,6 @@ export const DEFAULT_EPOCH_JD = 2451545.0;
 */
 export const GM_SUN_AU3_PER_DAY2 = 0.01720209895 * 0.01720209895;

-/**
- * Approximate planet/Sun mass ratios for the major planets that host moons in `bodies.json`.
- * Used to derive each planet's gravitational parameter (for propagating its moons) as
- * `GM_SUN_AU3_PER_DAY2 * massRatio`. Precise enough for visualization; not JPL-grade.
- */
-const PLANET_TO_SUN_MASS_RATIO: Record<string, number> = {
-  earth: 3.003e-6,
-  mars: 3.227e-7,
-  jupiter: 9.545e-4,
-  saturn: 2.857e-4,
-  uranus: 4.365e-5,
-  neptune: 5.151e-5
-};
-
-/**
- * Gravitational parameter (AU^3/day^2) to use when propagating a body's orbit: the Sun's
- * for planets/dwarfs/exoplanets, or the host planet's (derived from its Sun mass ratio) for
- * moons. Falls back to the Sun's GM if `parentBodyId` isn't a known planet.
- */
-export function gmForParent(parentBodyId: string | undefined): number {
-  if (!parentBodyId) {
-    return GM_SUN_AU3_PER_DAY2;
-  }
-  const massRatio = PLANET_TO_SUN_MASS_RATIO[parentBodyId];
-  return massRatio ? GM_SUN_AU3_PER_DAY2 * massRatio : GM_SUN_AU3_PER_DAY2;
-}
-
 /** Converts a JS `Date` into a Julian date (days), for driving the Kepler propagator "now". */
 export function dateToJulianDate(date: Date = new Date()): number {
  return date.getTime() / 86400000 + 2440587.5;
@@ -4,6 +4,7 @@ import {
  distanceBetween,
  eclipticToEquatorial,
  equatorialToEcliptic,
+  laplacePlaneToEquatorial,
  OBLIQUITY_J2000_DEG,
  parallaxMasToParsecs,
  parseSexagesimal,
@@ -206,6 +207,26 @@ describe('eclipticToEquatorial', () => {
  });
 });

+describe('laplacePlaneToEquatorial', () => {
+  const RAD = Math.PI / 180;
+  /** Jupiter's moons' Laplace pole, as JPL gives it for Io. */
+  const POLE = { raDeg: 268.057, decDeg: 64.495 };
+
+  it('sends the plane’s own pole to the right ascension and declination it is named by', () => {
+    const pole = laplacePlaneToEquatorial({ x: 0, y: 0, z: 1 }, POLE);
+    expect(Math.asin(pole.z) / RAD).toBeCloseTo(POLE.decDeg, 9);
+    expect(((Math.atan2(pole.y, pole.x) / RAD) + 360) % 360).toBeCloseTo(POLE.raDeg, 9);
+  });
+
+  it('counts the node from where the plane rises through the equator, 90 degrees past the pole', () => {
+    const node = laplacePlaneToEquatorial({ x: 1, y: 0, z: 0 }, POLE);
+    expect(node.z).toBeCloseTo(0, 12);
+    expect(((Math.atan2(node.y, node.x) / RAD) + 360) % 360).toBeCloseTo((POLE.raDeg + 90) % 360, 9);
+    // Rising: a quarter-turn on along the plane is north of the equator.
+    expect(laplacePlaneToEquatorial({ x: 0, y: 1, z: 0 }, POLE).z).toBeGreaterThan(0);
+  });
+});
+
 describe('equatorialToEcliptic', () => {
  it('is the exact inverse of eclipticToEquatorial', () => {
    for (const point of [
@@ -59,8 +59,8 @@ export const OBLIQUITY_J2000_DEG = 23.4392911;
 *
 * The app has to span both because its two sources disagree. Star positions come from HYG as
 * equatorial coordinates, which `raDecDistanceToXyz` produces and which the galaxy view renders
- * directly. Orbital elements come from JPL Horizons, whose default reference plane for element
- * output is the ecliptic — the ETL never overrides it. The two are tilted
+ * directly. The planets' and the Moon's orbital elements are JPL mean elements against the J2000
+ * ecliptic. The two are tilted
 * {@link OBLIQUITY_J2000_DEG} apart about the shared vernal-equinox axis, so orbits have to be
 * rotated before they can share a scene with the stars.
 */
@@ -78,6 +78,29 @@ export function eclipticToEquatorial(position: CartesianCoordinates): CartesianC
  };
 }

+/**
+ * Rotates a vector from a moon's local **Laplace plane** frame into the equatorial one.
+ *
+ * JPL gives the giant planets' moons against the plane their orbits precess about, which lies
+ * between the planet's equator and its orbit, and names it by its pole. The frame's x axis is
+ * where that plane rises through the ICRF equator, at right ascension 90 degrees past the pole's,
+ * which is what the node is counted from; its z axis is the pole, 90 degrees less its declination
+ * away from the celestial one. Read against the ecliptic instead, Io was up to 2.8 degrees from
+ * where Horizons has it between 1950 and 2100, Phobos 54 and Titan 127: their nodes are counted
+ * from a different line altogether.
+ */
+export function laplacePlaneToEquatorial(position: CartesianCoordinates, pole: { raDeg: number; decDeg: number }): CartesianCoordinates {
+  const tilt = (90 - pole.decDeg) * DEG_TO_RAD;
+  const node = (pole.raDeg + 90) * DEG_TO_RAD;
+  const y = position.y * Math.cos(tilt) - position.z * Math.sin(tilt);
+  const z = position.y * Math.sin(tilt) + position.z * Math.cos(tilt);
+  return {
+    x: position.x * Math.cos(node) - y * Math.sin(node),
+    y: position.x * Math.sin(node) + y * Math.cos(node),
+    z
+  };
+}
+
 /** Inverse of {@link eclipticToEquatorial}. */
 export function equatorialToEcliptic(position: CartesianCoordinates): CartesianCoordinates {
  const obliquity = OBLIQUITY_J2000_DEG * DEG_TO_RAD;
@@ -4,9 +4,11 @@ import { GM_SUN_AU3_PER_DAY2, DEFAULT_EPOCH_JD } from './constants';
 import {
  gravitationalParameterFromPeriod,
  isPropagatableOrbit,
+  meanElementsAt,
  meanMotionRadPerDay,
  orbitEllipsePoints,
  orbitalPeriodDays,
+  positionAtEpoch,
  positionAtTrueAnomaly,
  propagateOrbit,
  resolveGravitationalParameter,
@@ -131,6 +133,46 @@ describe('propagateOrbit', () => {
  });
 });

+describe('meanElementsAt', () => {
+  /** A circle in the reference plane, prograde (0) or retrograde (180), whose node turns. */
+  function circle(inclinationDeg: number) {
+    return { semiMajorAxisAu: 1, eccentricity: 0, inclinationDeg, longitudeOfAscendingNodeDeg: 0, argumentOfPeriapsisDeg: 0, meanAnomalyAtEpochDeg: 0, epochJd: DEFAULT_EPOCH_JD };
+  }
+  const RATES = { meanMotionDegPerDay: 10, longitudeOfAscendingNodeDegPerDay: 0.5, argumentOfPeriapsisDegPerDay: 0.2 };
+
+  /** Longitude in the reference plane a day on, in degrees, signed. */
+  function longitudeAfterOneDay(inclinationDeg: number): number {
+    const { x, y } = positionAtEpoch(meanElementsAt(circle(inclinationDeg), RATES, DEFAULT_EPOCH_JD + 1));
+    return (Math.atan2(y, x) * 180) / Math.PI;
+  }
+
+  it('goes round at its mean motion however its node and periapsis turn', () => {
+    expect(longitudeAfterOneDay(0)).toBeCloseTo(10, 9);
+  });
+
+  it('goes round a retrograde orbit backwards at the same rate, the node’s turning added back', () => {
+    // Taking the node off as for a prograde orbit made this 9 degrees, and Triton drifted a
+    // degree a year from where Horizons has it.
+    expect(longitudeAfterOneDay(180)).toBeCloseTo(-10, 9);
+  });
+
+  it('turns the node and periapsis at their own rates, and dates the result', () => {
+    const later = meanElementsAt(circle(0), RATES, DEFAULT_EPOCH_JD + 4);
+    expect(later.longitudeOfAscendingNodeDeg).toBeCloseTo(2, 12);
+    expect(later.argumentOfPeriapsisDeg).toBeCloseTo(0.8, 12);
+    expect(later.epochJd).toBe(DEFAULT_EPOCH_JD + 4);
+  });
+
+  it('adds Standish’s b T² + c cos(fT) + s sin(fT) to the mean anomaly', () => {
+    const terms = { b: -0.00012452, c: 0.0606406, s: -0.35635438, f: 38.35125 };
+    const T = 0.7;
+    const withTerms = meanElementsAt(circle(0), { ...RATES, meanAnomalyTerms: terms }, DEFAULT_EPOCH_JD + T * 36525);
+    const without = meanElementsAt(circle(0), RATES, DEFAULT_EPOCH_JD + T * 36525);
+    const f = (terms.f * T * Math.PI) / 180;
+    expect(withTerms.meanAnomalyAtEpochDeg - without.meanAnomalyAtEpochDeg).toBeCloseTo(terms.b * T * T + terms.c * Math.cos(f) + terms.s * Math.sin(f), 9);
+  });
+});
+
 describe('orbitEllipsePoints', () => {
  it('samples a closed loop whose distances stay within the periapsis/apoapsis bounds', () => {
    const elements = resolveOrbitalElements({ semiMajorAxisAu: 5, eccentricity: 0.4 });
@@ -1,9 +1,10 @@
 import { CartesianCoordinates } from './coordinates';
 import { DEFAULT_EPOCH_JD, GM_SUN_AU3_PER_DAY2 } from './constants';
-import { OrbitalElements } from '../models/body.model';
+import { MeanElementRates, OrbitalElements } from '../models/body.model';

 const DEG_TO_RAD = Math.PI / 180;
 const TWO_PI = Math.PI * 2;
+const DAYS_PER_JULIAN_CENTURY = 36525;

 /**
 * Fills in the elements the Kepler propagator needs but that some sources (e.g. exoplanets,
@@ -208,17 +209,63 @@ export function positionAtTrueAnomaly(elements: OrbitalElements, trueAnomalyRad:
 }

 /**
- * Propagates `elements` to Julian date `epochJdEval`, returning the body's position (AU)
- * relative to its central body. This is the app's "current epoch" evaluation used for live
- * (and future time-scrubbable) positions, as opposed to {@link orbitEllipsePoints} which
+ * The rates of an orbit that only goes round: Kepler's mean motion from the central mass, with
+ * nothing turning. What an exoplanet has, since the archive publishes no precession.
+ */
+export function keplerRates(semiMajorAxisAu: number, gmAu3PerDay2: number): MeanElementRates {
+  return {
+    meanMotionDegPerDay: meanMotionRadPerDay(semiMajorAxisAu, gmAu3PerDay2) / DEG_TO_RAD,
+    longitudeOfAscendingNodeDegPerDay: 0,
+    argumentOfPeriapsisDegPerDay: 0
+  };
+}
+
+/**
+ * The elements at `epochJdEval`, each moved from its epoch at its own rate, and returned with that
+ * date as their epoch — so {@link positionAtEpoch} places the body, and the node and periapsis
+ * say where to draw the orbit it is on.
+ *
+ * The mean anomaly is what is left of the body's motion once the node and periapsis have turned:
+ * `meanMotionDegPerDay` is how fast it goes round in space, and a periapsis that has moved on is
+ * that much further to reach. On a retrograde orbit, past 90 degrees, the body runs against the
+ * direction the node is counted in, so the node's turning is added back rather than taken off.
+ * Taken off, Triton — whose node turns half a degree a year — drifted a degree a year from where
+ * Horizons has it, 105 degrees by 2100.
+ */
+export function meanElementsAt(elements: OrbitalElements, rates: MeanElementRates, epochJdEval: number): OrbitalElements {
+  const days = epochJdEval - elements.epochJd;
+  const node = rates.longitudeOfAscendingNodeDegPerDay * days;
+  const periapsis = rates.argumentOfPeriapsisDegPerDay * days;
+  const nodeAlongOrbit = elements.inclinationDeg > 90 ? -node : node;
+  const terms = rates.meanAnomalyTerms;
+  const centuries = days / DAYS_PER_JULIAN_CENTURY;
+  const extra = terms
+    ? terms.b * centuries * centuries + terms.c * Math.cos(terms.f * centuries * DEG_TO_RAD) + terms.s * Math.sin(terms.f * centuries * DEG_TO_RAD)
+    : 0;
+  return {
+    semiMajorAxisAu: elements.semiMajorAxisAu + (rates.semiMajorAxisAuPerDay ?? 0) * days,
+    eccentricity: elements.eccentricity + (rates.eccentricityPerDay ?? 0) * days,
+    inclinationDeg: elements.inclinationDeg + (rates.inclinationDegPerDay ?? 0) * days,
+    longitudeOfAscendingNodeDeg: elements.longitudeOfAscendingNodeDeg + node,
+    argumentOfPeriapsisDeg: elements.argumentOfPeriapsisDeg + periapsis,
+    meanAnomalyAtEpochDeg: elements.meanAnomalyAtEpochDeg + rates.meanMotionDegPerDay * days - periapsis - nodeAlongOrbit + extra,
+    epochJd: epochJdEval
+  };
+}
+
+/** Where `elements` put the body at their own epoch (AU, relative to the central body). */
+export function positionAtEpoch(elements: OrbitalElements): CartesianCoordinates {
+  const eccentricAnomalyRad = solveEccentricAnomaly(elements.meanAnomalyAtEpochDeg * DEG_TO_RAD, elements.eccentricity);
+  return positionAtTrueAnomaly(elements, trueAnomalyFromEccentricAnomaly(eccentricAnomalyRad, elements.eccentricity));
+}
+
+/**
+ * Propagates `elements` to Julian date `epochJdEval` around a central mass, returning the body's
+ * position (AU) relative to its central body, as opposed to {@link orbitEllipsePoints} which
 * samples the fixed orbit shape independent of time.
 */
 export function propagateOrbit(elements: OrbitalElements, gmAu3PerDay2: number, epochJdEval: number): CartesianCoordinates {
-  const meanMotion = meanMotionRadPerDay(elements.semiMajorAxisAu, gmAu3PerDay2);
-  const meanAnomalyRad = elements.meanAnomalyAtEpochDeg * DEG_TO_RAD + meanMotion * (epochJdEval - elements.epochJd);
-  const eccentricAnomalyRad = solveEccentricAnomaly(meanAnomalyRad, elements.eccentricity);
-  const trueAnomalyRad = trueAnomalyFromEccentricAnomaly(eccentricAnomalyRad, elements.eccentricity);
-  return positionAtTrueAnomaly(elements, trueAnomalyRad);
+  return positionAtEpoch(meanElementsAt(elements, keplerRates(elements.semiMajorAxisAu, gmAu3PerDay2), epochJdEval));
 }

 /**
@@ -0,0 +1,110 @@
+import { describe, expect, it } from 'vitest';
+
+import { parsePlanetMeanElements, parseSatelliteMeanElements } from './mean-elements';
+
+/** Standish's p_elem_t2.txt, cut to the lines that matter here, as JPL published them. */
+const TABLE_2 = `Keplerian elements and their rates, with respect to the mean ecliptic and equinox of J2000,
+valid for the time-interval 3000 BC -- 3000 AD.  NOTE: the computation of M for Jupiter through
+Pluto *must* be augmented by the additional terms given in Table 2b (below).
+
+EM Bary   1.00000018      0.01673163     -0.00054346      100.46691572    102.93005885     -5.11260389
+         -0.00000003     -0.00003661     -0.01337178    35999.37306329      0.31795260     -0.24123856
+Jupiter   5.20248019      0.04853590      1.29861416       34.33479152     14.27495244    100.29282654
+         -0.00002864      0.00018026     -0.00322699     3034.90371757      0.18199196      0.13024619
+Pluto    39.48686035      0.24885238     17.14104260      238.96535011    224.09702598    110.30167986
+          0.00449751      0.00006016      0.00000501      145.18042903     -0.00968827     -0.00809981
+
+Table 2b.
+Jupiter   -0.00012452    0.06064060   -0.35635438   38.35125000
+Pluto     -0.01262724
+`;
+
+/** The satellite page's markup around three rows, as the 2021 page served it. */
+const SATELLITES = `
+<td align="left" nowrap><b>Satellites of Earth</b></td>
+<td align="right" nowrap><b>jump to:</b> <a href="#earth">Earth</a>, <a href="#mars">Mars</a></td>
+<H3>Mean <a href="?glossary&term=ecliptic">ecliptic</a> orbital elements</H3>
+Epoch 2000 Jan. 1.50 TT<BR>
+<TR ALIGN=right><TD ALIGN=left>Moon</TD>
+<TD>384400.</TD><TD>0.0554</TD><TD>318.15</TD><TD>135.27</TD><TD>5.16</TD><TD>125.08</TD>
+<TD>13.176358</TD><TD>27.322</TD><TD>5.997</TD><TD>18.600</TD>
+<TD ALIGN=right><A HREF="#ref1">1</A></TD></TR>
+<td align="left" nowrap><b>Satellites of Jupiter</b></td>
+<td align="right" nowrap><b>jump to:</b> <a href="#earth">Earth</a>, <a href="#mars">Mars</a></td>
+<H3>Mean orbital elements referred to the local <a href="?glossary&term=lp">Laplace planes</a></H3>
+Epoch 1997 Jan. 16.00 TT<BR>
+<TR ALIGN=right><TD ALIGN=left>Io</TD><TD>421800.</TD><TD>0.0041</TD>
+<TD>84.129</TD><TD>342.021</TD><TD>0.036</TD><TD>43.977</TD><TD>203.4889583</TD>
+<TD>1.769</TD><TD>1.625</TD><TD>7.420</TD><TD>268.057</TD><TD>64.495</TD>
+<TD>0.000</TD>
+<TD ALIGN=right><A HREF="#ref11">11</A></TD></TR>
+<td align="left" nowrap><b>Satellites of Neptune</b></td>
+<td align="right" nowrap><b>jump to:</b> <a href="#earth">Earth</a>, <a href="#mars">Mars</a></td>
+<H3>Mean orbital elements referred to the local <a href="?glossary&term=lp">Laplace planes</a></H3>
+Epoch 2000 Jan.  1.50 TT<BR>
+<TR ALIGN=right><TD ALIGN=left>Triton</TD><TD>354759.</TD><TD>0.0000</TD>
+<TD>66.142</TD><TD>352.257</TD><TD>156.865</TD><TD>177.608</TD>
+<TD>61.2572638</TD><TD>5.877</TD><TD>386.371</TD><TD>687.446</TD>
+<TD>299.456</TD><TD>43.414</TD><TD>0.010</TD>
+<TD ALIGN=right><A HREF="#ref54">54</A></TD></TR>
+`;
+
+describe('parsePlanetMeanElements', () => {
+  it('turns Standish’s longitudes into the argument of periapsis and mean anomaly', () => {
+    const { orbit } = parsePlanetMeanElements(TABLE_2, 'jupiter');
+    expect(orbit.argumentOfPeriapsisDeg).toBeCloseTo(14.27495244 - 100.29282654, 8);
+    expect(orbit.meanAnomalyAtEpochDeg).toBeCloseTo(34.33479152 - 14.27495244, 8);
+    expect(orbit.epochJd).toBe(2451545);
+  });
+
+  it('gives the rates per day, the mean motion being the mean longitude’s', () => {
+    const { rates } = parsePlanetMeanElements(TABLE_2, 'earth');
+    // 35 999.373 degrees a century is the sidereal year.
+    expect(360 / rates.meanMotionDegPerDay).toBeCloseTo(365.2564, 4);
+    expect(rates.argumentOfPeriapsisDegPerDay * 36525).toBeCloseTo(0.3179526 + 0.24123856, 8);
+  });
+
+  it('carries Table 2b’s terms for Jupiter and beyond, and none for the inner planets', () => {
+    expect(parsePlanetMeanElements(TABLE_2, 'jupiter').rates.meanAnomalyTerms).toEqual({ b: -0.00012452, c: 0.0606406, s: -0.35635438, f: 38.35125 });
+    expect(parsePlanetMeanElements(TABLE_2, 'earth').rates.meanAnomalyTerms).toBeUndefined();
+  });
+
+  it('reads Pluto’s row, not the note above the table that starts a line with its name', () => {
+    const pluto = parsePlanetMeanElements(TABLE_2, 'pluto');
+    expect(pluto.orbit.semiMajorAxisAu).toBe(39.48686035);
+    expect(pluto.rates.meanAnomalyTerms).toEqual({ b: -0.01262724, c: 0, s: 0, f: 0 });
+  });
+});
+
+describe('parseSatelliteMeanElements', () => {
+  it('reads a Laplace-plane row with its pole and its section’s epoch', () => {
+    const io = parseSatelliteMeanElements(SATELLITES, 'Jupiter', 'Io', true);
+    expect(io.laplacePole).toEqual({ raDeg: 268.057, decDeg: 64.495 });
+    expect(io.orbit.epochJd).toBe(2450464.5);
+    expect(io.rates.meanMotionDegPerDay).toBe(203.4889583);
+    expect(io.orbitSource).toBe('JPL SSD satellite mean elements, epoch 1997 Jan 16');
+  });
+
+  it('reads the Moon against the ecliptic, with no pole', () => {
+    const moon = parseSatelliteMeanElements(SATELLITES, 'Earth', 'Moon', false);
+    expect(moon.laplacePole).toBeUndefined();
+    expect(moon.orbit.epochJd).toBe(2451545);
+    expect(moon.orbit.semiMajorAxisAu * 149597870.7).toBeCloseTo(384400, 3);
+  });
+
+  it('regresses a prograde node and advances a periapsis, as the planet’s oblateness turns them', () => {
+    const { rates } = parseSatelliteMeanElements(SATELLITES, 'Earth', 'Moon', false);
+    expect(rates.longitudeOfAscendingNodeDegPerDay).toBeCloseTo(-360 / (18.6 * 365.25), 9);
+    expect(rates.argumentOfPeriapsisDegPerDay).toBeCloseTo(360 / (5.997 * 365.25), 9);
+  });
+
+  it('advances the node of a retrograde orbit', () => {
+    const { rates } = parseSatelliteMeanElements(SATELLITES, 'Neptune', 'Triton', false);
+    expect(rates.longitudeOfAscendingNodeDegPerDay).toBeCloseTo(360 / (687.446 * 365.25), 9);
+  });
+
+  it('turns the periapsis backwards where a resonance holds it', () => {
+    const { rates } = parseSatelliteMeanElements(SATELLITES, 'Jupiter', 'Io', true);
+    expect(rates.argumentOfPeriapsisDegPerDay).toBeCloseTo(-360 / (1.625 * 365.25), 9);
+  });
+});
@@ -0,0 +1,151 @@
+import { MeanElementRates, OrbitalElements } from '../models/body.model';
+
+/**
+ * Reads JPL's two tables of mean orbital elements, which the ETL fetches (see
+ * `tools/etl/lib/mean-elements.ts`), into the elements and rates `bodies.json` carries.
+ */
+
+const KM_PER_AU = 149597870.7;
+const J2000_JD = 2451545.0;
+const DAYS_PER_JULIAN_CENTURY = 36525;
+const DAYS_PER_JULIAN_YEAR = 365.25;
+
+const PLANET_ORBIT_SOURCE = 'JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000';
+
+export interface MeanOrbit {
+  orbit: OrbitalElements;
+  rates: MeanElementRates;
+  laplacePole?: { raDeg: number; decDeg: number };
+  orbitSource: string;
+}
+
+/** Table 2a's name for each planet; Earth's row is the Earth-Moon barycentre, 4 700 km off Earth. */
+const PLANET_ROW_NAMES: Record<string, string> = {
+  mercury: 'Mercury',
+  venus: 'Venus',
+  earth: 'EM Bary',
+  mars: 'Mars',
+  jupiter: 'Jupiter',
+  saturn: 'Saturn',
+  uranus: 'Uranus',
+  neptune: 'Neptune',
+  pluto: 'Pluto'
+};
+
+function numbers(text: string): number[] {
+  return text.trim().split(/\s+/).map(Number);
+}
+
+/**
+ * Reads one planet's row pair from Table 2a, and its Table 2b terms where it has them. The
+ * elements are Standish's own — a, e, I, mean longitude L, longitude of perihelion ϖ, node Ω —
+ * turned into the argument of periapsis ϖ - Ω and mean anomaly L - ϖ the propagator takes.
+ */
+export function parsePlanetMeanElements(text: string, bodyId: string): MeanOrbit {
+  const name = PLANET_ROW_NAMES[bodyId];
+  const lines = text.split(/\r?\n/);
+  // A row is the name followed by a number: the notes above the table start a line with "Pluto" too.
+  const rows = lines.flatMap((line, index) => (name && new RegExp(`^${name}\\s+-?[\\d.]`).test(line) ? [index] : []));
+  if (rows.length === 0) {
+    throw new Error(`No row for ${bodyId} in Standish's Table 2a.`);
+  }
+  const values = [...numbers(lines[rows[0]].slice(name.length)), ...numbers(lines[rows[0] + 1])];
+  if (values.length !== 12 || !values.every(Number.isFinite)) {
+    throw new Error(`Standish's Table 2a rows for ${name} did not parse: ${values.join(' ')}`);
+  }
+  const [a, e, inclination, meanLongitude, perihelion, node, aRate, eRate, inclinationRate, meanLongitudeRate, perihelionRate, nodeRate] = values;
+  // Table 2b repeats the name further down, with b, c, s, f (Pluto has b alone).
+  const extra = rows[1] === undefined ? undefined : numbers(lines[rows[1]].slice(name.length));
+  const [b = 0, c = 0, s = 0, f = 0] = extra ?? [];
+
+  return {
+    orbit: {
+      semiMajorAxisAu: a,
+      eccentricity: e,
+      inclinationDeg: inclination,
+      longitudeOfAscendingNodeDeg: node,
+      argumentOfPeriapsisDeg: perihelion - node,
+      meanAnomalyAtEpochDeg: meanLongitude - perihelion,
+      epochJd: J2000_JD
+    },
+    rates: {
+      meanMotionDegPerDay: meanLongitudeRate / DAYS_PER_JULIAN_CENTURY,
+      longitudeOfAscendingNodeDegPerDay: nodeRate / DAYS_PER_JULIAN_CENTURY,
+      argumentOfPeriapsisDegPerDay: (perihelionRate - nodeRate) / DAYS_PER_JULIAN_CENTURY,
+      semiMajorAxisAuPerDay: aRate / DAYS_PER_JULIAN_CENTURY,
+      eccentricityPerDay: eRate / DAYS_PER_JULIAN_CENTURY,
+      inclinationDegPerDay: inclinationRate / DAYS_PER_JULIAN_CENTURY,
+      ...(extra ? { meanAnomalyTerms: { b, c, s, f } } : {})
+    },
+    orbitSource: PLANET_ORBIT_SOURCE
+  };
+}
+
+const MONTHS = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'];
+
+/** `1997 Jan. 16.00` as a Julian date. TT and TDB differ by under two milliseconds. */
+function julianDate(year: number, month: string, day: number): number {
+  const monthIndex = MONTHS.indexOf(month);
+  if (monthIndex < 0) {
+    throw new Error(`Unknown month ${month}.`);
+  }
+  return Date.UTC(year, monthIndex, 1) / 86400000 + 2440587.5 + day - 1;
+}
+
+/**
+ * Reads one moon's row from the satellite page: `a e w M i node n P Pw Pnode`, then the Laplace
+ * pole `RA Dec Tilt` where the section is referred to one, then a reference number.
+ *
+ * The page gives the two precession periods as magnitudes, so their sense is supplied here. A
+ * node driven by the planet's oblateness regresses on a prograde orbit and advances on a
+ * retrograde one, and the orbit's inclination says which. A periapsis advances — except where a
+ * resonance forces the eccentricity, which `apsidesRegress` names: Io's and Europa's are held to
+ * the line of their conjunctions, which turns backwards at 2 n(Europa) - n(Io) = 0.74 degrees a
+ * day, and that is exactly the 1.625- and 1.394-year periods the table gives for them. Read as
+ * advancing, Io was 0.9 degrees out and Europa 2.1.
+ */
+export function parseSatelliteMeanElements(html: string, planetName: string, moonName: string, apsidesRegress: boolean): MeanOrbit {
+  const text = html.replace(/<[^>]+>/g, ' ').replace(/&nbsp;/g, ' ').replace(/\s+/g, ' ');
+  const section = text.indexOf(`Satellites of ${planetName} jump to`);
+  if (section < 0) {
+    throw new Error(`No section for the satellites of ${planetName}.`);
+  }
+  const row = text.slice(section).match(new RegExp(` ${moonName} ((?:-?[\\d.]+ )+)`));
+  if (!row || row.index === undefined) {
+    throw new Error(`No row for ${moonName} among the satellites of ${planetName}.`);
+  }
+  const before = text.slice(section, section + row.index);
+  const epoch = [...before.matchAll(/Epoch (\d{4}) (\w{3})\. ([\d.]+) T/g)].at(-1);
+  if (!epoch) {
+    throw new Error(`No epoch above ${moonName}'s row.`);
+  }
+  const laplace = before.lastIndexOf('Laplace plane') > before.lastIndexOf('Mean ecliptic');
+  const values = numbers(row[1]);
+  const expected = laplace ? 14 : 11;
+  if (values.length !== expected || !values.every(Number.isFinite)) {
+    throw new Error(`${moonName}'s row has ${values.length} numbers, ${expected} expected: ${row[1]}`);
+  }
+  const [aKm, e, periapsis, meanAnomaly, inclination, node, meanMotion, , periapsisPeriodYears, nodePeriodYears, raDeg, decDeg] = values;
+  const nodeSense = inclination > 90 ? 1 : -1;
+  const periapsisSense = apsidesRegress ? -1 : 1;
+  const perDay = (periodYears: number): number => (periodYears > 0 ? 360 / (periodYears * DAYS_PER_JULIAN_YEAR) : 0);
+
+  return {
+    orbit: {
+      semiMajorAxisAu: aKm / KM_PER_AU,
+      eccentricity: e,
+      inclinationDeg: inclination,
+      longitudeOfAscendingNodeDeg: node,
+      argumentOfPeriapsisDeg: periapsis,
+      meanAnomalyAtEpochDeg: meanAnomaly,
+      epochJd: julianDate(Number(epoch[1]), epoch[2], Number(epoch[3]))
+    },
+    rates: {
+      meanMotionDegPerDay: meanMotion,
+      longitudeOfAscendingNodeDegPerDay: nodeSense * perDay(nodePeriodYears),
+      argumentOfPeriapsisDegPerDay: periapsisSense * perDay(periapsisPeriodYears)
+    },
+    ...(laplace ? { laplacePole: { raDeg, decDeg } } : {}),
+    orbitSource: `JPL SSD satellite mean elements, epoch ${epoch[1]} ${epoch[2]} ${Math.floor(Number(epoch[3]))}`
+  };
+}
@@ -1,7 +1,7 @@
 /**
- * Osculating Keplerian orbital elements at a reference epoch. Positions are derived
- * client-side by propagating these elements forward/backward from `epochJd` (see
- * `shared/astro/kepler.ts`), rather than fetching per-frame positions.
+ * Keplerian orbital elements at a reference epoch. Positions are derived client-side by
+ * propagating these elements forward/backward from `epochJd` (see `shared/astro/kepler.ts`),
+ * rather than fetching per-frame positions.
 */
 export interface OrbitalElements {
  semiMajorAxisAu: number;
@@ -14,8 +14,36 @@ export interface OrbitalElements {
 }

 /**
- * A solar-system planet, moon, or dwarf planet, sourced from JPL Horizons/SSD orbital
- * elements. `systemStarId` links back to the HYG star index (the Sun, see `SUN_STAR_ID`).
+ * How a body's mean elements move away from their epoch, per day.
+ *
+ * Mean elements rather than one osculating set, because the map's clock runs decades in minutes.
+ * An osculating orbit is exact at its instant and drifts from then on: fed to Kepler with a mass
+ * ratio, the Moon's went round in 27.70 days instead of 27.32 and was 66 degrees out after a year.
+ * A mean set carries its own measured motion, and the slow turning of its node and periapsis, so
+ * it holds for as long as its source was fit over.
+ */
+export interface MeanElementRates {
+  /**
+   * How fast the body goes round in space, in degrees per day: the rate of its mean longitude.
+   * 360 over this is its sidereal period.
+   */
+  meanMotionDegPerDay: number;
+  longitudeOfAscendingNodeDegPerDay: number;
+  argumentOfPeriapsisDegPerDay: number;
+  semiMajorAxisAuPerDay?: number;
+  eccentricityPerDay?: number;
+  inclinationDegPerDay?: number;
+  /**
+   * Standish's extra terms in the mean anomaly of Jupiter and beyond, `b T² + c cos(f T) +
+   * s sin(f T)` degrees, with T in Julian centuries from the epoch and f in degrees per century:
+   * the great-inequality wobble his 3000 BC to AD 3000 fit needs on top of its linear rates.
+   */
+  meanAnomalyTerms?: { b: number; c: number; s: number; f: number };
+}
+
+/**
+ * A solar-system planet, moon, or dwarf planet: JPL mean orbital elements, and JPL Horizons
+ * physical data. `systemStarId` links back to the HYG star index (the Sun, see `SUN_STAR_ID`).
 */
 export interface BodyRecord {
  id: string;
@@ -23,7 +51,18 @@ export interface BodyRecord {
  name: string;
  kind: 'planet' | 'moon' | 'dwarf';
  radiusKm: number;
+  /** Mean elements at `orbit.epochJd`, moving at `rates`. */
  orbit: OrbitalElements;
+  rates: MeanElementRates;
+  /**
+   * The pole of the plane a moon's elements are measured against, where that is its local
+   * Laplace plane: right ascension and declination in the ICRF. The node is then counted from
+   * where that plane crosses the ICRF equator. Absent means the J2000 ecliptic, as for the
+   * planets and the Moon.
+   */
+  laplacePole?: { raDeg: number; decDeg: number };
+  /** Where the elements come from and the span they hold over, as the card prints it. */
+  orbitSource: string;
  /**
   * For `kind: 'moon'`, the `id` of the planet it orbits — its `orbit` is expressed
   * relative to that planet, not heliocentrically. Undefined for planets/dwarfs.