Move Mimas, Tethys and Phobos along their orbits by the terms their IAU W already carried

The IAU W of a locked moon follows its mean longitude, so a term of W that is the moon running
ahead of and behind its mean motion is its orbit's too. Two were in bodies.json's W and in no
orbit: the 71-year libration of the Mimas-Tethys 4:2 resonance, -44.85 degrees on Mimas and
+2.23 on Tethys on the angle S5 = 316.45 + 506.2 T of pck00011.tpc, and Phobos's tidal
quadratic, 9.536e-9 degrees a day squared about J2000. JPL's satellite table has a column for
neither. orbitalTermsOfPrimeMeridian now turns each into the row's meanAnomalyTerms about the
row's own epoch (Phobos's 1950 row gets the quadratic re-centred, which adds to its mean motion
and mean anomaly at the epoch), and the ETL takes them for the three moons named in their specs.

Against Horizons: Mimas on 2026 May 27, near the libration's extreme, 2.24 degrees instead of
43.3; Tethys the same day 0.18 instead of 2.05; Phobos in 2100 1.25 instead of 11.1. On the
ETL's 2025-01-01 check Mimas is 1.56 degrees, so its named 46-degree ceiling is gone. The renderer
spec freezes the Mimas and Phobos vectors.

The day-equals-orbit check checked a number that turns no locked moon: since cdf474b every one is
turned by its IAU W. build.ts now checks what is drawn instead: the east longitude of the planet
on the moon's IAU body-fixed frame, from where the mean elements put it, every 135 days from 1950
to 2100. Measured: at most 6.70 degrees (the Moon's own eccentricity swing), with Mimas 10.15
(its physical libration, which Horizons shows too), Iapetus 18.33 (its row sits 9.4 degrees
behind Horizons) and Proteus 8.18 under their own ceilings. The lock holds only over that span:
Proteus's W turns 6.3e-7 of its rate slower than its orbit, which drifts its face 74 degrees by
AD 3000, and Mimas's W runs 6e-5 degrees a day ahead of the table's mean motion. Without the new
terms the check fails: Phobos's face turned 13.78 degrees from Mars by 2100, and Mimas's 54.5.

The comments that promised the lock "however long the clock runs" now say what the check covers.
The ceilings comment in build.ts gives Hyperion's and Nereid's worst offsets sampled daily (22.2
and 11.2 degrees, where twelve New Year's Days gave 20.2 and 2.6).

Controls: dropping the libration's cosine term, the quadratic, the quadratic's re-centring or
swapping sine and cosine each fails its named test; the ETL without Mimas's term fails at 44.68
degrees and without Phobos's at 13.78.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>

src/app/features/galaxy-system/system-orbits-renderer.spec.ts

@@ -506,6 +506,10 @@ describe('solar-system bodies against Horizons', () => {
     charon: {kind: 'moon', orbit: {semiMajorAxisAu: 0.00013095774631236113, eccentricity: 0.0002, inclinationDeg: 0.08, longitudeOfAscendingNodeDeg: 26.928, argumentOfPeriapsisDeg: 146.106, meanAnomalyAtEpochDeg: 131.07, epochJd: 2451545}, rates: {meanMotionDegPerDay: 56.362521, longitudeOfAscendingNodeDegPerDay: -0.00010926638529337138, argumentOfPeriapsisDegPerDay: 0.00009683851540842405}, laplacePole: {raDeg: 132.993, decDeg: -6.163}, parentBodyId: 'pluto', massRatio: 0.1220485755631374},
     venus: {kind: 'planet', orbit: {semiMajorAxisAu: 0.72332102, eccentricity: 0.00676399, inclinationDeg: 3.39777545, longitudeOfAscendingNodeDeg: 76.67261496, argumentOfPeriapsisDeg: 55.094942169999996, meanAnomalyAtEpochDeg: 50.21215136999999, epochJd: 2451545}, rates: {meanMotionDegPerDay: 1.6021304750882956, longitudeOfAscendingNodeDegPerDay: -0.000007467261875427789, argumentOfPeriapsisDegPerDay: 0.000009022264750171116, semiMajorAxisAuPerDay: -7.118412046543463e-12, eccentricityPerDay: -1.398220396988364e-9, inclinationDegPerDay: 1.1908008213552361e-8}},
     mars: {kind: 'planet', orbit: {semiMajorAxisAu: 1.52371243, eccentricity: 0.09336511, inclinationDeg: 1.85181869, longitudeOfAscendingNodeDeg: 49.71320984, argumentOfPeriapsisDeg: -73.63065768, meanAnomalyAtEpochDeg: 19.3493162, epochJd: 2451545}, rates: {meanMotionDegPerDay: 0.5240328362061601, longitudeOfAscendingNodeDegPerDay: -0.000007351794934976044, argumentOfPeriapsisDegPerDay: 0.00001973334866529774, semiMajorAxisAuPerDay: 2.6557152635181385e-11, eccentricityPerDay: 2.5048596851471597e-9, inclinationDegPerDay: -1.9842765229295004e-7}},
+    // Mimas and Phobos carry the terms of their IAU W that are motion along the orbit: the Mimas-Tethys
+    // libration and Phobos's tidal acceleration (see `orbitalTermsOfPrimeMeridian`).
+    mimas: {kind: 'moon', orbit: {semiMajorAxisAu: 0.0012402516100785653, eccentricity: 0.0196, inclinationDeg: 1.574, longitudeOfAscendingNodeDeg: 173.027, argumentOfPeriapsisDeg: 332.499, meanAnomalyAtEpochDeg: 14.848, epochJd: 2451545}, rates: {meanMotionDegPerDay: 381.9944948, longitudeOfAscendingNodeDegPerDay: -0.9996209770462032, argumentOfPeriapsisDegPerDay: 1.9992419540924065, meanAnomalyTerms: {b: 0, c: 30.901081394463684, s: -32.50608664008527, f: 506.2}}, laplacePole: {raDeg: 40.589, decDeg: 83.536}, parentBodyId: 'saturn'},
+    phobos: {kind: 'moon', orbit: {semiMajorAxisAu: 0.00006267468885838895, eccentricity: 0.0151, inclinationDeg: 1.075, longitudeOfAscendingNodeDeg: 207.784, argumentOfPeriapsisDeg: 150.057, meanAnomalyAtEpochDeg: 94.2394819925, epochJd: 2433282.5}, rates: {meanMotionDegPerDay: 1128.8444085925948, longitudeOfAscendingNodeDegPerDay: -0.43579001784832494, argumentOfPeriapsisDegPerDay: 0.871002371303956, meanAnomalyTerms: {b: 12.721927969999998, c: 0, s: 0, f: 0}}, laplacePole: {raDeg: 317.671, decDeg: 52.893}, parentBodyId: 'mars'},
   };
   // The IAU WGCCRE 2015 rotational elements bodies.json carries for them, from pck00011.tpc.
   const ROTATION: Record<string, RotationalElements> = {
@@ -521,7 +525,9 @@ describe('solar-system bodies against Horizons', () => {
   // Each ceiling sits just above what these elements measure on that date: Earth 0.003 degrees,
   // Jupiter 0.063, Saturn 0.164, Pluto 0.054, the Moon 0.72 (no mean ellipse has its evection or
   // variation), Io 0.021, Europa 0.036, Titan 0.014, Triton 0.137, Titania 0.62 (against Uranus's
-  // equator, 120 years from its 1980 epoch), Charon 0.37.
+  // equator, 120 years from its 1980 epoch), Charon 0.37, Mimas 2.24 on 2026 May 27, when its libration has it
+  // 44 degrees ahead of its mean motion (43.3 without the term), and Phobos 1.25 in 2100 (11.1 without its
+  // tidal acceleration).
   const HORIZONS: Array<[id: string, jd: number, x: number, y: number, z: number, maxDeg: number]> = [
     ['earth', 2488069.5, -0.1574071329883954, 0.890666220858489, 0.3859132211165683, 0.02],
     ['jupiter', 2433282.5, 3.406605247558555, -3.425997624196318, -1.551719750032203, 0.1],
@@ -534,6 +540,8 @@ describe('solar-system bodies against Horizons', () => {
     ['triton', 2488069.5, -0.001421151845853369, -0.0001894510477241482, 0.001888790702926415, 0.2],
     ['titania', 2488069.5, -0.00151919968294745, -0.0003387914082135071, 0.002465657830788125, 0.75],
     ['charon', 2488069.5, -0.00003046411046017432, -0.000009404114448552256, 0.0001270457155789907, 0.5],
+    ['mimas', 2461187.5, -3.840685116962088e-4, 1.182766232573268e-3, -2.006928678577268e-5, 3],
+    ['phobos', 2488069.5, 4.269855297288105e-5, -2.759158950396543e-5, -3.816269137755616e-5, 1.5],
   ];
 
   function record(id: string): BodyRecord {

src/app/shared/astro/rotational-elements.spec.ts

@@ -1,6 +1,8 @@
 import { describe, expect, it } from 'vitest';
 
-import { orientationAt, parsePckRotationalElements } from './rotational-elements';
+import { meanElementsAt } from './kepler';
+import { orbitalTermsOfPrimeMeridian, orientationAt, parsePckRotationalElements } from './rotational-elements';
+import { MeanElementRates, OrbitalElements, RotationalElements } from '../models/body.model';
 
 // Excerpts of pck00011.tpc as NAIF publishes it: prose, then data blocks.
 const KERNEL = String.raw`KPL/PCK
@@ -131,3 +133,68 @@ describe('orientationAt', () => {
     expect(orientationAt(phobos, J2000 + days).primeMeridianDeg).toBeCloseTo(((expected % 360) + 360) % 360, 5);
   });
 });
+
+describe('orbitalTermsOfPrimeMeridian', () => {
+  const J2000 = 2451545.0;
+  // Mimas's and Phobos's rows and W as pck00011.tpc and JPL's table give them, with each mean
+  // motion set to W's rate, so that only the terms can part the two.
+  const MIMAS: RotationalElements = {
+    poleRaDeg: [40.66, -0.036],
+    poleDecDeg: [83.52, -0.004],
+    primeMeridianDeg: [333.46, 381.994555, 0],
+    terms: [
+      { angleDeg: [177.4, -36505.5], ra: 13.56, dec: -1.53, pm: -13.48 },
+      { angleDeg: [316.45, 506.2], ra: 0, dec: 0, pm: -44.85 }
+    ]
+  };
+  const PHOBOS: RotationalElements = {
+    poleRaDeg: [317.67071657, -0.10844326, 0],
+    poleDecDeg: [52.88627266, -0.06134706, 0],
+    primeMeridianDeg: [35.1877444, 1128.84475928, 9.536137031212154e-9]
+  };
+  const orbit = (epochJd: number): OrbitalElements => ({
+    semiMajorAxisAu: 0.001,
+    eccentricity: 0.02,
+    inclinationDeg: 1.5,
+    longitudeOfAscendingNodeDeg: 170,
+    argumentOfPeriapsisDeg: 60,
+    meanAnomalyAtEpochDeg: 10,
+    epochJd
+  });
+
+  /** The moon's mean longitude less W, which a locked moon holds still whatever the date. */
+  function lead(elements: RotationalElements, epochJd: number, angleRate: number | undefined, jd: number): number {
+    const { meanAnomalyTerms, meanMotionDegPerDay, meanAnomalyDeg } = orbitalTermsOfPrimeMeridian(elements, epochJd, angleRate);
+    const start = orbit(epochJd);
+    const rates: MeanElementRates = {
+      meanMotionDegPerDay: elements.primeMeridianDeg[1] + meanMotionDegPerDay,
+      longitudeOfAscendingNodeDegPerDay: -1,
+      argumentOfPeriapsisDegPerDay: 2,
+      meanAnomalyTerms
+    };
+    const moved = meanElementsAt({ ...start, meanAnomalyAtEpochDeg: start.meanAnomalyAtEpochDeg + meanAnomalyDeg }, rates, jd);
+    const longitude = moved.longitudeOfAscendingNodeDeg + moved.argumentOfPeriapsisDeg + moved.meanAnomalyAtEpochDeg;
+    // W with the pole's nodding term left out: that one is the pole's, not the orbit's.
+    const w = orientationAt({ ...elements, terms: elements.terms?.filter((term) => term.angleDeg[1] === angleRate) }, jd).primeMeridianDeg;
+    return (((longitude - w) % 360) + 540) % 360 - 180;
+  }
+
+  it('moves Mimas along its orbit by the libration its W carries, over the 71 years it takes', () => {
+    const atEpoch = lead(MIMAS, J2000, 506.2, J2000);
+    for (const years of [-130, -40, 17.8, 35.5, 100]) {
+      expect(lead(MIMAS, J2000, 506.2, J2000 + years * 365.25)).toBeCloseTo(atEpoch, 8);
+    }
+  });
+
+  it('speeds Phobos up by the tidal quadratic its W carries about J2000, from a row whose epoch is 1950', () => {
+    const epoch = 2433282.5;
+    const atEpoch = lead(PHOBOS, epoch, undefined, epoch);
+    for (const years of [-150, 0, 50, 100, 150, 400]) {
+      expect(lead(PHOBOS, epoch, undefined, J2000 + years * 365.25)).toBeCloseTo(atEpoch, 6);
+    }
+  });
+
+  it('refuses a term W does not carry', () => {
+    expect(() => orbitalTermsOfPrimeMeridian(PHOBOS, J2000, 506.2)).toThrow();
+  });
+});

src/app/shared/astro/rotational-elements.ts

@@ -121,3 +121,44 @@ export function orientationAt(elements: RotationalElements, jdTdb: number): { po
   }
   return { poleRaDeg, poleDecDeg, primeMeridianDeg: ((primeMeridianDeg % 360) + 360) % 360 };
 }
+
+/**
+ * What a locked moon's W says of its going round that a row of mean elements leaves out, as terms
+ * of that row. W follows the moon's mean longitude, so a term of W that is the moon running ahead
+ * of and behind its mean motion, rather than its pole nodding, is its orbit's too. Two are here,
+ * and JPL's satellite table has a column for neither: Mimas's -44.85 degrees and Tethys's +2.23 on
+ * the angle that turns 506.2 degrees a century, the 71-year libration of their 4:2 resonance, and
+ * Phobos's quadratic, 12.72 degrees per century squared about J2000, the tidal acceleration
+ * drawing it in. Carried by W and not by the orbit, they left the drawn Mimas up to 45 degrees from
+ * where Horizons has it and its face as far from Saturn, and Phobos 11 degrees out by 2100.
+ *
+ * `angleRateDegPerCentury` names the term by its angle's rate; W's quadratic, where it has one, is
+ * always taken. Both come back about `epochJd`, which is where `meanAnomalyTerms` counts T from:
+ * the sine as its `c` and `s`, and the quadratic re-centred from J2000 onto that epoch, as `b` plus
+ * what the re-centring adds to the mean motion and to the mean anomaly at the epoch.
+ */
+export function orbitalTermsOfPrimeMeridian(
+  elements: RotationalElements,
+  epochJd: number,
+  angleRateDegPerCentury?: number
+): { meanAnomalyTerms: { b: number; c: number; s: number; f: number }; meanMotionDegPerDay: number; meanAnomalyDeg: number } {
+  const epochCenturies = (epochJd - J2000_JD) / DAYS_PER_JULIAN_CENTURY;
+  // W turns clockwise about the pole the IAU names where its rate is negative; the orbit does not.
+  const sense = Math.sign(elements.primeMeridianDeg[1]);
+  const quadratic = sense * (elements.primeMeridianDeg[2] ?? 0) * DAYS_PER_JULIAN_CENTURY * DAYS_PER_JULIAN_CENTURY;
+  let sine = { c: 0, s: 0, f: 0 };
+  if (angleRateDegPerCentury !== undefined) {
+    const term = elements.terms?.find((candidate) => candidate.angleDeg[1] === angleRateDegPerCentury);
+    if (!term || (term.angleDeg[2] ?? 0) !== 0) {
+      throw new Error(`No term of W turns linearly at ${angleRateDegPerCentury} degrees a century.`);
+    }
+    const phase = (term.angleDeg[0] + term.angleDeg[1] * epochCenturies) * DEG_TO_RAD;
+    sine = { c: sense * term.pm * Math.sin(phase), s: sense * term.pm * Math.cos(phase), f: term.angleDeg[1] };
+  }
+  // q (T + T0)², T from the epoch and T0 the epoch from J2000, is q T² + 2 q T0 T + q T0².
+  return {
+    meanAnomalyTerms: { b: quadratic, ...sine },
+    meanMotionDegPerDay: (2 * quadratic * epochCenturies) / DAYS_PER_JULIAN_CENTURY,
+    meanAnomalyDeg: quadratic * epochCenturies * epochCenturies
+  };
+}

src/assets/data/bodies.json

@@ -689,13 +689,19 @@
       "inclinationDeg": 1.075,
       "longitudeOfAscendingNodeDeg": 207.784,
       "argumentOfPeriapsisDeg": 150.057,
-      "meanAnomalyAtEpochDeg": 91.059,
+      "meanAnomalyAtEpochDeg": 94.2394819925,
       "epochJd": 2433282.5
     },
     "rates": {
-      "meanMotionDegPerDay": 1128.8447569,
+      "meanMotionDegPerDay": 1128.8444085925948,
       "longitudeOfAscendingNodeDegPerDay": -0.43579001784832494,
-      "argumentOfPeriapsisDegPerDay": 0.871002371303956
+      "argumentOfPeriapsisDegPerDay": 0.871002371303956,
+      "meanAnomalyTerms": {
+        "b": 12.721927969999998,
+        "c": 0,
+        "s": 0,
+        "f": 0
+      }
     },
     "laplacePole": {
       "raDeg": 317.671,
@@ -703,7 +709,7 @@
     },
     "orbitSource": "JPL SSD satellite mean elements, epoch 1950 Jan 1",
     "parentBodyId": "mars",
-    "rotationPeriodHours": 7.653842521027348,
+    "rotationPeriodHours": 7.653844882637157,
     "rotationalElements": {
       "poleRaDeg": [
         317.67071657,
@@ -1166,7 +1172,13 @@
     "rates": {
       "meanMotionDegPerDay": 381.9944948,
       "longitudeOfAscendingNodeDegPerDay": -0.9996209770462032,
-      "argumentOfPeriapsisDegPerDay": 1.9992419540924065
+      "argumentOfPeriapsisDegPerDay": 1.9992419540924065,
+      "meanAnomalyTerms": {
+        "b": 0,
+        "c": 30.901081394463684,
+        "s": -32.50608664008527,
+        "f": 506.2
+      }
     },
     "laplacePole": {
       "raDeg": 40.589,
@@ -1276,7 +1288,13 @@
     "rates": {
       "meanMotionDegPerDay": 190.6979109,
       "longitudeOfAscendingNodeDegPerDay": -0.1978374715711675,
-      "argumentOfPeriapsisDegPerDay": 0.3958338487419905
+      "argumentOfPeriapsisDegPerDay": 0.3958338487419905,
+      "meanAnomalyTerms": {
+        "b": 0,
+        "c": -1.5364417281974139,
+        "s": 1.6162446646017872,
+        "f": 506.2
+      }
     },
     "laplacePole": {
       "raDeg": 40.578,

tools/etl/build.ts

@@ -150,20 +150,27 @@ const KM_PER_AU = 149597870.7;
 const DEG_TO_RAD = Math.PI / 180;
 
 /**
- * The moons whose table row cannot come within that, each for a reason no mean ellipse carries,
- * with a ceiling just above its worst offset from Horizons at twelve dates from 1980 to 2100:
+ * The moons whose table row cannot come within that on this one date, each for a reason no mean
+ * ellipse carries, with a ceiling just above its offset here. What each reaches elsewhere is
+ * larger, and neither this check nor twelve New Year's Days sampled from 1980 to 2100 see it:
  *
- * - Mimas, 44.7 degrees: its resonance with Tethys swings its mean longitude 44 degrees either
- *   way over 70.8 years, and the table has no column for it (Tethys, on the other end, swings 2).
- * - Hyperion, 20.2: held in a 4:3 resonance by Titan; the row's eccentricity, 0.0232, is less than
- *   a quarter of the 0.105 JPL's current table gives.
- * - Iapetus, 10.1: the row sits 9.4 degrees behind Horizons at its own epoch, 2000 Jan 1.5, and
- *   keeps that offset; its plane agrees with Horizons' to 0.07 degrees and its period to 0.001 per
- *   cent, so the fault is in the row's longitude, which this has no second source to correct.
- * - Nereid, 2.6: an eccentricity of 0.75, the largest here, which a mean ellipse follows least
- *   well: under 0.9 degrees in every year measured but 2025 and 2030 (2.6 and 2.3) and 2100 (1.7).
+ * - Hyperion, 9.4 degrees here and 22.2 at worst, sampled every other day from 1980 to 2100: held
+ *   in a 4:3 resonance by Titan; the row's eccentricity, 0.0232, is less than a quarter of the
+ *   0.105 JPL's current table gives.
+ * - Iapetus, 9.6 here and 10.1 at worst: the row sits 9.4 degrees behind Horizons at its own epoch,
+ *   2000 Jan 1.5, and keeps that offset; its plane agrees with Horizons' to 0.07 degrees and its
+ *   period to 0.001 per cent, so the fault is in the row's longitude, which this has no second
+ *   source to correct.
+ * - Nereid, 2.6 here and 11.2 at worst, sampled daily: an eccentricity of 0.75, the largest here,
+ *   which a mean ellipse follows least well near periapsis, where the true anomaly runs ten times
+ *   faster than the mean. Its 360-day year put all twelve New Year's Days, under 2.6, far from
+ *   periapsis; sampled daily it is past 2.6 on 92 days in 2010-2020 and on 295 in 2040-2050, each
+ *   near a periapsis, and under 0.4 on most days.
+ *
+ * Mimas's swing of 44 degrees either way, the libration of its resonance with Tethys, which also
+ * needed a ceiling here once, is now in its orbit; see `orbitFromW` in `fetchSolarSystem.ts`.
  */
-const MOON_OFFSET_CEILINGS_DEG: Record<string, number> = { mimas: 46, hyperion: 21, iapetus: 11, nereid: 3 };
+const MOON_OFFSET_CEILINGS_DEG: Record<string, number> = { hyperion: 21, iapetus: 11, nereid: 3 };
 
 /**
  * The bodies the IAU WGCCRE 2015 report gives no rotational elements for: Hyperion tumbles, and
@@ -192,6 +199,36 @@ const DAY_OFFSET_CEILINGS: Record<string, number> = { neptune: 0.01 };
  */
 const MAX_OBLIQUITY_OFFSET_DEG = 0.1;
 
+/**
+ * How far from its planet a locked moon's drawn face may turn: the east longitude, on the IAU's
+ * body-fixed frame, of the direction to the planet from where the mean elements put the moon,
+ * sampled every 135 days from 1950 to 2100, where both the tables and the IAU's elements hold.
+ *
+ * Measured on this catalogue: at most 6.70 degrees (the Moon, whose longitude swings 6.3 either
+ * way with its eccentricity; Horizons has the same). Three need their own: Mimas 10.15, whose
+ * physical libration W carries and Horizons shows as 5 to 9 degrees at its true place; Iapetus
+ * 18.33, whose row sits 9.4 degrees behind Horizons; and Proteus 8.18, whose W turns 6.3e-7 of
+ * its rate slower than its orbit, a drift of 74 degrees by AD 3000. What this catches is an orbit
+ * and a W that go round at different rates: the tidal acceleration W carried and the orbit did not
+ * turned Phobos 13.8 degrees from Mars by 2100, and the Mimas-Tethys libration Mimas 54.5.
+ */
+const MAX_SUB_PLANET_LONGITUDE_DEG = 7;
+const SUB_PLANET_CEILINGS_DEG: Record<string, number> = { mimas: 11, iapetus: 19, proteus: 9 };
+const LOCK_DATES_JD = Array.from({ length: 407 }, (_, index) => 2433282.5 + index * 135);
+
+/** The planet's east longitude on a moon's IAU body-fixed frame, from the moon's mean place, at a TDB date. */
+function subPlanetLongitudeDeg(body: BodyRecord, jd: number): number {
+  const own = positionAtEpoch(meanElementsAt(body.orbit, body.rates, jd));
+  const place = body.laplacePole ? laplacePlaneToEquatorial(own, body.laplacePole) : eclipticToEquatorial(own);
+  const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(body.rotationalElements!, jd);
+  const pole = { raDeg: poleRaDeg, decDeg: poleDecDeg };
+  const w = primeMeridianDeg * DEG_TO_RAD;
+  const meridian = laplacePlaneToEquatorial({ x: Math.cos(w), y: Math.sin(w), z: 0 }, pole);
+  const east = laplacePlaneToEquatorial({ x: -Math.sin(w), y: Math.cos(w), z: 0 }, pole);
+  const along = (axis: { x: number; y: number; z: number }) => -(place.x * axis.x + place.y * axis.y + place.z * axis.z);
+  return Math.atan2(along(east), along(meridian)) / DEG_TO_RAD;
+}
+
 function angleBetweenDeg(a: { x: number; y: number; z: number }, b: { x: number; y: number; z: number }): number {
   const cosine = (a.x * b.x + a.y * b.y + a.z * b.z) / (Math.hypot(a.x, a.y, a.z) * Math.hypot(b.x, b.y, b.z));
   return (Math.acos(Math.min(1, Math.max(-1, cosine))) * 180) / Math.PI;
@@ -271,13 +308,16 @@ function validateBodies(bodies: BodyRecord[], horizonsOrbits: Map<string, Orbita
           `Moon ${body.id} does not keep one face to its planet, yet turns once in ${body.rotationPeriodHours} hours against an orbit of ${orbitHours}.`
         );
       } else {
-        // Every other moon here is tidally locked: its day is its orbit, from the same mean motion
-        // that carries it round, or its face turns away from its planet: the Kepler period of the
-        // osculating orbit this used to take would turn the Moon's five degrees an orbit.
+        // Every other moon here is tidally locked, and drawn by its orbit and its IAU W: the two
+        // have to agree, or its face turns away from its planet.
+        assertCondition(rotation !== undefined, `Moon ${body.id} is locked but has no W to keep its face to its planet by.`);
+        const ceiling = SUB_PLANET_CEILINGS_DEG[body.id] ?? MAX_SUB_PLANET_LONGITUDE_DEG;
+        const worst = Math.max(...LOCK_DATES_JD.map((jd) => Math.abs(subPlanetLongitudeDeg(body, jd))));
         assertCondition(
-          body.rotationPeriodHours !== undefined && Math.abs(body.rotationPeriodHours - orbitHours) <= orbitHours * 1e-9,
-          `Moon ${body.id} turns once in ${body.rotationPeriodHours} hours but goes round in ${orbitHours} — it will not keep one face to its planet.`
+          worst <= ceiling,
+          `Moon ${body.id} turns its face up to ${worst.toFixed(2)} degrees from its planet between 1950 and 2100 (at most ${ceiling} expected) — its orbit and its W disagree.`
         );
+        spins.push(`${body.id} faces ${worst.toFixed(2)}`);
       }
       if (body.massRatio !== undefined) {
         // The pair's barycentre, which the planet's elements place, must lie outside the planet —

tools/etl/fetchSolarSystem.ts

@@ -5,7 +5,7 @@ import { SUN_STAR_ID } from '../../src/app/shared/models/star.model';
 import { fetchHorizonsBody } from './lib/horizons';
 import { MeanOrbit, parsePlanetMeanElements, parseSatelliteMeanElements, parseSmallBodyElements } from '../../src/app/shared/astro/mean-elements';
 import { fetchPlanetMeanElementsText, fetchSatelliteMeanElementsHtml, fetchSmallBodyAnswer } from './lib/mean-elements';
-import { MIN_PERIODIC_TERM_DEG, parsePckRotationalElements } from '../../src/app/shared/astro/rotational-elements';
+import { MIN_PERIODIC_TERM_DEG, orbitalTermsOfPrimeMeridian, parsePckRotationalElements } from '../../src/app/shared/astro/rotational-elements';
 import { fetchPckText } from './lib/pck';
 import { dataPath, ensureDataDir } from './lib/paths';
 
@@ -42,8 +42,17 @@ interface BodySpec {
   nodeOffsetDeg?: number;
   epochJd?: number;
   periodDays?: number;
+  /**
+   * The terms of the IAU's W that are this locked moon's motion along its orbit, which its row has
+   * no column for: W's quadratic, and the term whose angle turns at `angleRateDegPerCentury`, if
+   * given. See `orbitalTermsOfPrimeMeridian`.
+   */
+  orbitFromW?: { angleRateDegPerCentury?: number };
 }
 
+/** S5 in pck00011.tpc, 316.45 + 506.2 T: the libration of Mimas and Tethys in their 4:2 resonance. */
+const MIMAS_TETHYS_LIBRATION = { angleRateDegPerCentury: 506.2 };
+
 /**
  * The poles of the equators JPL refers Uranus's and Pluto's moons to, from the IAU WGCCRE 2015
  * report, each taken at the end the table's inclinations are measured from (Titania 0.079
@@ -82,15 +91,15 @@ const BODY_SPECS: BodySpec[] = [
   { id: 'haumea', name: 'Haumea', kind: 'dwarf', horizonsCommand: '136108;', center: '500@10', sbdb: 'Haumea', radiusKm: 797.6 },
   { id: 'makemake', name: 'Makemake', kind: 'dwarf', horizonsCommand: '136472;', center: '500@10', sbdb: 'Makemake', radiusKm: 715 },
   { id: 'moon', name: 'Moon', kind: 'moon', horizonsCommand: '301', center: '500@399', parentBodyId: 'earth' },
-  { id: 'phobos', name: 'Phobos', kind: 'moon', horizonsCommand: '401', center: '500@499', parentBodyId: 'mars' },
+  { id: 'phobos', name: 'Phobos', kind: 'moon', horizonsCommand: '401', center: '500@499', parentBodyId: 'mars', orbitFromW: {} },
   { id: 'deimos', name: 'Deimos', kind: 'moon', horizonsCommand: '402', center: '500@499', parentBodyId: 'mars' },
   { id: 'io', name: 'Io', kind: 'moon', horizonsCommand: '501', center: '500@599', parentBodyId: 'jupiter', apsidesRegress: true },
   { id: 'europa', name: 'Europa', kind: 'moon', horizonsCommand: '502', center: '500@599', parentBodyId: 'jupiter', apsidesRegress: true },
   { id: 'ganymede', name: 'Ganymede', kind: 'moon', horizonsCommand: '503', center: '500@599', parentBodyId: 'jupiter' },
   { id: 'callisto', name: 'Callisto', kind: 'moon', horizonsCommand: '504', center: '500@599', parentBodyId: 'jupiter' },
-  { id: 'mimas', name: 'Mimas', kind: 'moon', horizonsCommand: '601', center: '500@699', parentBodyId: 'saturn' },
+  { id: 'mimas', name: 'Mimas', kind: 'moon', horizonsCommand: '601', center: '500@699', parentBodyId: 'saturn', orbitFromW: MIMAS_TETHYS_LIBRATION },
   { id: 'enceladus', name: 'Enceladus', kind: 'moon', horizonsCommand: '602', center: '500@699', parentBodyId: 'saturn' },
-  { id: 'tethys', name: 'Tethys', kind: 'moon', horizonsCommand: '603', center: '500@699', parentBodyId: 'saturn' },
+  { id: 'tethys', name: 'Tethys', kind: 'moon', horizonsCommand: '603', center: '500@699', parentBodyId: 'saturn', orbitFromW: MIMAS_TETHYS_LIBRATION },
   { id: 'dione', name: 'Dione', kind: 'moon', horizonsCommand: '604', center: '500@699', parentBodyId: 'saturn' },
   { id: 'rhea', name: 'Rhea', kind: 'moon', horizonsCommand: '605', center: '500@699', parentBodyId: 'saturn' },
   { id: 'titan', name: 'Titan', kind: 'moon', horizonsCommand: '606', center: '500@699', parentBodyId: 'saturn' },
@@ -148,6 +157,15 @@ export async function fetchSolarSystem(): Promise<{ bodies: BodyRecord[]; horizo
     horizonsOrbits.set(spec.id, result.orbit);
     gmById.set(spec.id, result.gmKm3PerS2);
 
+    // NAIF numbers a small body 2 000 000 past its catalogue number: Ceres, "1;" to Horizons, is 2000001.
+    const naifId = spec.horizonsCommand.endsWith(';') ? 2_000_000 + Number.parseInt(spec.horizonsCommand, 10) : Number(spec.horizonsCommand);
+    const rotation = parsePckRotationalElements(pck, naifId);
+    if (!rotation) {
+      console.warn(`  no IAU rotational elements for ${spec.name}; its pole and meridian are not known.`);
+    } else if (rotation.skippedDeg.length > 0) {
+      console.log(`  ${spec.name}: ${rotation.skippedDeg.length} periodic terms under ${MIN_PERIODIC_TERM_DEG} degrees left out, the largest ${Math.max(...rotation.skippedDeg)}.`);
+    }
+
     const parentName = BODY_SPECS.find((candidate) => candidate.id === spec.parentBodyId)?.name;
     const smallBody = spec.sbdb ? parseSmallBodyElements(await fetchSmallBodyAnswer(spec.sbdb, `sbdb-${spec.id}.json`)) : undefined;
     const read: MeanOrbit =
@@ -155,7 +173,7 @@ export async function fetchSolarSystem(): Promise<{ bodies: BodyRecord[]; horizo
       (parentName
         ? parseSatelliteMeanElements(satelliteElements, parentName, spec.name, spec.apsidesRegress ?? false, spec.equatorPole)
         : parsePlanetMeanElements(planetElements, spec.id));
-    const mean: MeanOrbit = {
+    const corrected: MeanOrbit = {
       ...read,
       orbit: {
         ...read.orbit,
@@ -164,17 +182,28 @@ export async function fetchSolarSystem(): Promise<{ bodies: BodyRecord[]; horizo
       },
       rates: spec.periodDays ? { ...read.rates, meanMotionDegPerDay: 360 / spec.periodDays } : read.rates
     };
+    if (spec.orbitFromW && !rotation) {
+      throw new Error(`${spec.name}'s orbit takes terms from a W the kernel does not give.`);
+    }
+    const fromW = spec.orbitFromW && orbitalTermsOfPrimeMeridian(rotation!.elements, corrected.orbit.epochJd, spec.orbitFromW.angleRateDegPerCentury);
+    const mean: MeanOrbit = fromW
+      ? {
+          ...corrected,
+          orbit: { ...corrected.orbit, meanAnomalyAtEpochDeg: corrected.orbit.meanAnomalyAtEpochDeg + fromW.meanAnomalyDeg },
+          rates: { ...corrected.rates, meanMotionDegPerDay: corrected.rates.meanMotionDegPerDay + fromW.meanMotionDegPerDay, meanAnomalyTerms: fromW.meanAnomalyTerms }
+        }
+      : corrected;
     const radiusKm = smallBody?.radiusKm ?? spec.radiusKm ?? result.radiusKm;
     if (radiusKm === undefined) {
       console.warn(`  no physical radius found for ${spec.name}; defaulting to 0.`);
     }
 
     // A moon listed here is tidally locked unless its spec says otherwise, so its day is its
-    // orbit: the sidereal period from the same mean motion that carries it round, which keeps one
-    // face towards the parent however long the clock runs. Not every page says so — the Moon's
-    // gives a rate, Titan's and Proteus's nothing. The Kepler period of the osculating orbit this
-    // used to take, 27.70 days for the Moon, would now turn its face five degrees an orbit away
-    // from the orbit it is drawn on.
+    // orbit: the sidereal period from the same mean motion that carries it round. Not every page
+    // says so — the Moon's gives a rate, Titan's and Proteus's nothing. Every locked moon here is
+    // turned by its IAU W rather than by this day, and `build.ts` checks that W and the orbit keep
+    // its face to its planet from 1950 to 2100; this day is what the renderer would turn a moon
+    // without W by.
     const rotationPeriodHours = result.tidallyLocked || (spec.kind === 'moon' && !spec.spinsFreely)
       ? (360 / mean.rates.meanMotionDegPerDay) * HOURS_PER_DAY
       : smallBody
@@ -187,14 +216,6 @@ export async function fetchSolarSystem(): Promise<{ bodies: BodyRecord[]; horizo
     if (rotationPeriodHours === undefined) {
       console.warn(`  no rotation period found for ${spec.name}; it will not turn.`);
     }
-    // NAIF numbers a small body 2 000 000 past its catalogue number: Ceres, "1;" to Horizons, is 2000001.
-    const naifId = spec.horizonsCommand.endsWith(';') ? 2_000_000 + Number.parseInt(spec.horizonsCommand, 10) : Number(spec.horizonsCommand);
-    const rotation = parsePckRotationalElements(pck, naifId);
-    if (!rotation) {
-      console.warn(`  no IAU rotational elements for ${spec.name}; its pole and meridian are not known.`);
-    } else if (rotation.skippedDeg.length > 0) {
-      console.log(`  ${spec.name}: ${rotation.skippedDeg.length} periodic terms under ${MIN_PERIODIC_TERM_DEG} degrees left out, the largest ${Math.max(...rotation.skippedDeg)}.`);
-    }
 
     bodies.push({
       id: spec.id,