From 3b4fd1af6cfa2c36080b9fd1004a32faf593cfc0 Mon Sep 17 00:00:00 2001 From: Senrokai Date: Wed, 30 Sep 2026 15:57:54 +0200 Subject: [PATCH] Turn each locked moon at its orbit's rate and Iapetus's pole round its orbit, so they face their planets at every date the clock reaches The IAU gives a locked moon's W the mean motion of whichever orbit its authors had, and JPL's table has another. Near the present the difference is nothing; over the clock's AD 1 to 3000 it turned Proteus's far side to Neptune at AD 1 (146 degrees), Mimas 52 degrees from Saturn and Miranda 23. Iapetus was worse for another reason: its IAU pole is a straight line, 3.9 degrees a century in right ascension, through its orbit normal's 3 439-year circle round the Laplace pole, which by AD 1 has run past the celestial pole (Dec 97.9), 11 degrees off the orbit, with the face 87 degrees from Saturn. Mimas and Iapetus carry mission maps, so a wrong hemisphere was drawn facing Saturn. The ETL's lock check sampled only 1950-2100, so none of it failed. tools/etl/lib/locked-spin.ts, lockedToOrbit, called for every locked moon: - W's rate becomes the orbit's own mean motion, its constant moved so W is unchanged on 2025-01-01; the pole and every periodic term stay the IAU's. A W with a quadratic is left (Phobos's orbit already takes it; the Moon's is its tidal slowing, 0.75 degrees at AD 1). The kernel's rate must be within 1e-5 of the orbit's first (at most 3.4e-6, Iapetus). - Iapetus (poleFollowsOrbit): the pole follows its orbit normal, as a moon in a Cassini state does, in the IAU's own form: sines of the node's angle and four harmonics on right ascension, cosines on declination, fitted to the normal's circle and pinned to the IAU pole at the present; W takes sines of the same angles, fitted to hold the face where it is today. build.ts samples the lock over AD 1 to 3000 (8 114 dates, every 135 days) instead of 1950-2100, and subPlanetLongitudeDeg moved to the lib, shared by both. Measured on the real catalogue: at most 5.36 degrees (Titan) but the Moon 7.62 (its eccentricity, and W's quadratic at AD 1: a new named ceiling of 8), Mimas 8.94 (ceiling 11 -> 9.5) and Iapetus 15.95 (19 -> 16.5, 9.4 of it its row's lag); Proteus's own ceiling of 9 is gone, at 2.66. Iapetus's axis stays within 0.74 degrees of its orbit normal (11.06 before) and its pole is the IAU's at the present to 1e-4 degrees. Live on :4301, the longitude facing the planet at AD 1 / 1000 / 2025 / 2999: Proteus 2.6 / 2.6 / 2.6 / 2.6 (was -146.5 / -72.9 / 2.6 / 74.5), Iapetus -15.9 / -15.4 / -15.3 / -9.3 (-87.0 / -50.9 / -15.3 / 22.8), Mimas 4.1 / 6.8 / 5.7 / 8.4 (48.6 / 29.3 / 5.7 / -13.0), Miranda -0.1 / 2.2 / 0.0 / 1.4 (-23.4 / -9.6 / 0.0 / 12.6); the present is unchanged. The renderer spec now takes the rotational elements from bodies.json too, so no hand copy is left, and a new test turns Proteus, Miranda, Mimas and Iapetus to their planets at AD 1 and AD 3000. Guarded mutants, each run through the solar ETL and then the full suite on what it wrote: lockedToOrbit bypassed (validator: Mimas 52.30, ceiling 9.5; the new test fails), Iapetus on the IAU's straight pole (98.48), and its pole round the orbit without W's terms (73.13); each fails the new test and only it. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../system-orbits-renderer.spec.ts | 64 +++++--- src/app/shared/models/body.model.ts | 10 +- src/assets/data/bodies.json | 135 +++++++++++------ tools/etl/build.ts | 53 +++---- tools/etl/fetchSolarSystem.ts | 16 +- tools/etl/lib/locked-spin.ts | 141 ++++++++++++++++++ 6 files changed, 312 insertions(+), 107 deletions(-) create mode 100644 tools/etl/lib/locked-spin.ts diff --git a/src/app/features/galaxy-system/system-orbits-renderer.spec.ts b/src/app/features/galaxy-system/system-orbits-renderer.spec.ts index 38daa82..0c01ec2 100644 --- a/src/app/features/galaxy-system/system-orbits-renderer.spec.ts +++ b/src/app/features/galaxy-system/system-orbits-renderer.spec.ts @@ -8,7 +8,7 @@ import { DEFAULT_EPOCH_JD, GM_SUN_AU3_PER_DAY2, ttMinusUtSeconds } from '../../s import { keplerRates } from '../../shared/astro/kepler'; import { eclipticToEquatorial, laplacePlaneToEquatorial, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; import { orientationAt } from '../../shared/astro/rotational-elements'; -import { BodyRecord, RotationalElements } from '../../shared/models/body.model'; +import { BodyRecord } from '../../shared/models/body.model'; import { ExoplanetRecord } from '../../shared/models/exoplanet.model'; import { SystemOrbitsRenderer } from './system-orbits-renderer'; import { bodyTexturePath, loadCachedTexture } from '../../shared/rendering/texture-catalog'; @@ -593,29 +593,18 @@ describe('exoplanet size without a measured radius', () => { }); describe('solar-system bodies against Horizons', () => { - // The records the app ships, read from bodies.json, and Horizons' own positions for them (ICRF, - // AU; heliocentric for the planets, planet-centred for the moons) at dates across 1950-2100, so the - // whole path — the ETL's reading of the mean elements, their rates, the Laplace planes and the - // scene's frame — is checked against JPL's ephemeris rather than against itself. A hand copy of the - // records stood here, and an ETL that dropped Standish's a, e and i rates or Io's and Europa's - // backward periapses passed the whole suite on the data it wrote. Horizons' dates are TDB and the - // renderer's are the clock's UT, so each is handed over TT - UT earlier: 69.184 s today, 29 in 1950. + // The records the app ships, read from bodies.json with their IAU rotational elements, and + // Horizons' own positions for them (ICRF, AU; heliocentric for the planets, planet-centred for the + // moons) at dates across 1950-2100, so the whole path — the ETL's reading of the mean elements, + // their rates, the Laplace planes and the scene's frame — is checked against JPL's ephemeris rather + // than against itself. A hand copy of the records stood here, and an ETL that dropped Standish's a, + // e and i rates or Io's and Europa's backward periapses passed the whole suite on the data it + // wrote. Horizons' dates are TDB and the renderer's are the clock's UT, so each is handed over + // TT - UT earlier: 69.184 s today, 29 in 1950. const SHIPPED: BodyRecord[] = JSON.parse(readFileSync(`${process.cwd()}/src/assets/data/bodies.json`, 'utf8')); // Mimas and Phobos among them for the terms of their IAU W that are motion along the orbit: the // Mimas-Tethys libration and Phobos's tidal acceleration (see `orbitalTermsOfPrimeMeridian`). const IDS = ['earth', 'jupiter', 'saturn', 'neptune', 'pluto', 'moon', 'io', 'europa', 'titan', 'triton', 'uranus', 'titania', 'charon', 'venus', 'mars', 'mimas', 'phobos']; - // The IAU WGCCRE 2015 rotational elements bodies.json carries for them, from pck00011.tpc. - const ROTATION: Record = { - venus: {poleRaDeg: [272.76, 0, 0], poleDecDeg: [67.16, 0, 0], primeMeridianDeg: [160.2, -1.4813688, 0]}, - earth: {poleRaDeg: [0, -0.641, 0], poleDecDeg: [90, -0.557, 0], primeMeridianDeg: [190.147, 360.9856235, 0]}, - mars: {poleRaDeg: [317.269202, -0.10927547, 0], poleDecDeg: [54.432516, -0.05827105, 0], primeMeridianDeg: [176.049863, 350.891982443297, 0], terms: [{angleDeg: [79.398797, 0.5042615, 0], ra: 0.419057, dec: 0, pm: 0}, {angleDeg: [166.325722, 0.5042615, 0], ra: 0, dec: 1.591274, pm: 0}, {angleDeg: [95.391654, 0.5042615, 0], ra: 0, dec: 0, pm: 0.584542}]}, - jupiter: {poleRaDeg: [268.056595, -0.006499, 0], poleDecDeg: [64.495303, 0.002413, 0], primeMeridianDeg: [284.95, 870.536, 0]}, - io: {poleRaDeg: [268.05, -0.009, 0], poleDecDeg: [64.5, 0.003, 0], primeMeridianDeg: [200.39, 203.4889538, 0], terms: [{angleDeg: [283.9, 4850.7], ra: 0.094, dec: 0.04, pm: -0.085}, {angleDeg: [355.8, 1191.3], ra: 0.024, dec: 0.011, pm: -0.022}]}, - saturn: {poleRaDeg: [40.589, -0.036, 0], poleDecDeg: [83.537, -0.004, 0], primeMeridianDeg: [38.9, 810.7939024, 0]}, - uranus: {poleRaDeg: [257.311, 0, 0], poleDecDeg: [-15.175, 0, 0], primeMeridianDeg: [203.81, -501.1600928, 0]}, - pluto: {poleRaDeg: [132.993, 0, 0], poleDecDeg: [-6.163, 0, 0], primeMeridianDeg: [302.695, 56.3625225, 0]}, - moon: {poleRaDeg: [269.9949, 0.0031, 0], poleDecDeg: [66.5392, 0.013, 0], primeMeridianDeg: [38.3213, 13.17635815, -1.4e-12], terms: [{angleDeg: [125.045, -1935.5364525], ra: -3.8787, dec: 1.5419, pm: 3.561}, {angleDeg: [250.089, -3871.072905], ra: -0.1204, dec: 0.0239, pm: 0.1208}, {angleDeg: [260.008, 475263.3328725], ra: 0.07, dec: -0.0278, pm: -0.0642}, {angleDeg: [176.625, 487269.629985], ra: -0.0172, dec: 0.0068, pm: 0.0158}, {angleDeg: [357.529, 35999.0509575], ra: 0, dec: 0, pm: 0.0252}]}, - }; // 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 @@ -644,8 +633,8 @@ describe('solar-system bodies against Horizons', () => { ]; function record(id: string): BodyRecord { - const { kind, orbit, rates, laplacePole, parentBodyId, massRatio } = SHIPPED.find((body) => body.id === id)!; - return { id, systemStarId: 0, name: id, radiusKm: 1000, orbitSource: 'test', kind, orbit, rates, laplacePole, parentBodyId, massRatio, rotationalElements: ROTATION[id] }; + const { kind, orbit, rates, laplacePole, parentBodyId, massRatio, rotationalElements } = SHIPPED.find((body) => body.id === id)!; + return { id, systemStarId: 0, name: id, radiusKm: 1000, orbitSource: 'test', kind, orbit, rates, laplacePole, parentBodyId, massRatio, rotationalElements }; } const renderer = new SystemOrbitsRenderer(IDS.map(record), []); @@ -796,7 +785,7 @@ describe('solar-system bodies against Horizons', () => { /** Where the IAU puts a body's prime meridian at a TDB date, in the scene. */ function iauPrimeMeridian(id: string, jdTdb: number): THREE.Vector3 { - const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(ROTATION[id], jdTdb); + const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(SHIPPED.find((body) => body.id === id)!.rotationalElements!, jdTdb); const w = (primeMeridianDeg * Math.PI) / 180; const meridian = laplacePlaneToEquatorial({ x: Math.cos(w), y: Math.sin(w), z: 0 }, { raDeg: poleRaDeg, decDeg: poleDecDeg }); return new THREE.Vector3(meridian.x, meridian.y, meridian.z); @@ -907,3 +896,32 @@ describe('solar-system bodies against Horizons', () => { }); } }); + +describe('locked moons across the clock’s window', () => { + // As shipped, pole, W and all: the IAU gives each a W fitted near the present, and its rate is + // not quite its orbit's, nor Iapetus's pole a line for twenty centuries. + const shipped: BodyRecord[] = JSON.parse(readFileSync(`${process.cwd()}/src/assets/data/bodies.json`, 'utf8')); + const renderer = new SystemOrbitsRenderer(shipped.filter((body) => ['saturn', 'uranus', 'neptune', 'mimas', 'iapetus', 'miranda', 'proteus'].includes(body.id)), []); + + /** East longitude, on its map, of the point on a moon's drawn sphere that faces its planet. */ + function facingPlanet(id: string): number { + const moon = renderer.members.find((member) => member.id === id)!.marker; + // Its position is from the planet, which is its pivot; SphereGeometry wraps u = atan2(z, -x) / 2 pi. + const toPlanet = moon.position.clone().negate().applyQuaternion(moon.quaternion.clone().invert()); + const u = Math.atan2(toPlanet.z, -toPlanet.x) / (2 * Math.PI); + return ((((u - 0.5) * 360) % 360) + 540) % 360 - 180; + } + + it('keeps Proteus, Miranda, Mimas and Iapetus facing their planets at AD 1 and AD 3000', () => { + // Measured: Proteus 2.6 degrees at most over AD 1-3000, Miranda 2.8, Mimas 8.9, Iapetus 16 (9.4 + // of it the lag of the row its orbit is drawn from). On the IAU's own W and Iapetus's straight + // pole they were 146, 23, 49 and 87 degrees at AD 1. + for (const jd of [1721425.5, 2816787.4]) { + renderer.update(jd); + expect(Math.abs(facingPlanet('proteus'))).toBeLessThan(3); + expect(Math.abs(facingPlanet('miranda'))).toBeLessThan(3); + expect(Math.abs(facingPlanet('mimas'))).toBeLessThan(9.5); + expect(Math.abs(facingPlanet('iapetus'))).toBeLessThan(16.5); + } + }); +}); diff --git a/src/app/shared/models/body.model.ts b/src/app/shared/models/body.model.ts index aa607ea..25fe0f8 100644 --- a/src/app/shared/models/body.model.ts +++ b/src/app/shared/models/body.model.ts @@ -95,10 +95,12 @@ export interface BodyRecord { obliquityDeg?: number; /** * Where the body's pole points and which way its prime meridian faces at any date, from the IAU - * WGCCRE 2015 report (Archinal et al. 2018) as NAIF's `pck00011.tpc` carries it. Where present - * it alone sets how the body is drawn, and the ETL checks the period and obliquity above against - * it. Absent where the report gives none: Hyperion tumbles, and Nereid, Eris, Haumea and - * Makemake have no model. + * WGCCRE 2015 report (Archinal et al. 2018) as NAIF's `pck00011.tpc` carries it, but that a locked + * moon's W turns at its drawn orbit's rate and Iapetus's pole goes round with its orbit's, so they + * keep their faces to their planets over the clock's AD 1 to 3000 (see `lockedToOrbit` in the + * ETL). Where present it alone sets how the body is drawn, and the ETL checks the period and + * obliquity above against it. Absent where the report gives none: Hyperion tumbles, and Nereid, + * Eris, Haumea and Makemake have no model. */ rotationalElements?: RotationalElements; } diff --git a/src/assets/data/bodies.json b/src/assets/data/bodies.json index 55b0051..2fde121 100644 --- a/src/assets/data/bodies.json +++ b/src/assets/data/bodies.json @@ -809,8 +809,8 @@ 0 ], "primeMeridianDeg": [ - 79.39932954, - 285.16188899, + 79.49055322521608, + 285.161879, 0 ], "terms": [ @@ -906,8 +906,8 @@ 0 ], "primeMeridianDeg": [ - 200.39, - 203.4889538, + 200.34890824984421, + 203.4889583, 0 ], "terms": [ @@ -971,8 +971,8 @@ 0 ], "primeMeridianDeg": [ - 36.022, - 101.3747235, + 36.015607949990184, + 101.3747242, 0 ], "terms": [ @@ -1045,8 +1045,8 @@ 0 ], "primeMeridianDeg": [ - 44.064, - 50.3176081, + 44.07221835003116, + 50.3176072, 0 ], "terms": [ @@ -1119,8 +1119,8 @@ 0 ], "primeMeridianDeg": [ - 259.51, - 21.5710715, + 259.498129049991, + 21.5710728, 0 ], "terms": [ @@ -1199,8 +1199,8 @@ 0 ], "primeMeridianDeg": [ - 333.46, - 381.994555, + 334.00971630006546, + 381.9944948, 0 ], "terms": [ @@ -1264,8 +1264,8 @@ 0 ], "primeMeridianDeg": [ - 6.32, - 262.7318996, + 6.336436700062315, + 262.7318978, 0 ] } @@ -1315,8 +1315,8 @@ 0 ], "primeMeridianDeg": [ - 8.95, - 190.6979085, + 8.928084400003424, + 190.6979109, 0 ], "terms": [ @@ -1380,8 +1380,8 @@ 0 ], "primeMeridianDeg": [ - 357.6, - 131.5349316, + 357.60821835003117, + 131.5349307, 0 ] } @@ -1425,8 +1425,8 @@ 0 ], "primeMeridianDeg": [ - 235.16, - 79.6900478, + 235.1773498500081, + 79.6900459, 0 ], "terms": [ @@ -1481,8 +1481,8 @@ 0 ], "primeMeridianDeg": [ - 186.5855, - 22.5769768, + 186.5964577999983, + 22.5769756, 0 ] } @@ -1544,19 +1544,66 @@ "rotationPeriodHours": 1903.9469348834284, "rotationalElements": { "poleRaDeg": [ - 318.16, - -3.949, + 283.4205056238858, + 0, 0 ], "poleDecDeg": [ - 75.03, - -1.143, + 77.15116551386498, + 0, 0 ], "primeMeridianDeg": [ - 355.2, - 4.5379572, + 388.27970766227594, + 4.5379416, 0 + ], + "terms": [ + { + "angleDeg": [ + 261.10499999998126, + -10.46898128087986 + ], + "ra": -42.60490881707297, + "dec": 7.692796952777751, + "pm": 41.78612687698857 + }, + { + "angleDeg": [ + 522.2099999999625, + -20.93796256175972 + ], + "ra": -15.535995393995726, + "dec": 1.2845933786338362, + "pm": 15.53892376069157 + }, + { + "angleDeg": [ + 783.3149999999438, + -31.406943842639578 + ], + "ra": -7.628387941010053, + "dec": 0.45238994454766307, + "pm": 7.62837456853055 + }, + { + "angleDeg": [ + 1044.419999999925, + -41.87592512351944 + ], + "ra": -4.2134447101107595, + "dec": 0.20197385694504744, + "pm": 4.213441809239813 + }, + { + "angleDeg": [ + 1305.5249999999064, + -52.344906404399296 + ], + "ra": -2.482395388436937, + "dec": 0.10190190070186696, + "pm": 2.4609277085117798 + } ] } }, @@ -1644,8 +1691,8 @@ 0 ], "primeMeridianDeg": [ - 30.7, - -254.6906892, + 30.41144460000184, + -254.6906576, 0 ], "terms": [ @@ -1727,8 +1774,8 @@ 0 ], "primeMeridianDeg": [ - 156.22, - -142.8356681, + 156.12685870007942, + -142.8356579, 0 ], "terms": [ @@ -1792,8 +1839,8 @@ 0 ], "primeMeridianDeg": [ - 108.05, - -86.8688923, + 108.00982140004953, + -86.8688879, 0 ], "terms": [ @@ -1857,8 +1904,8 @@ 0 ], "primeMeridianDeg": [ - 77.74, - -41.3514316, + 77.67607950003162, + -41.3514246, 0 ], "terms": [ @@ -1913,8 +1960,8 @@ 0 ], "primeMeridianDeg": [ - 6.77, - -26.7394932, + 6.729821400017093, + -26.7394888, 0 ], "terms": [ @@ -1969,8 +2016,8 @@ 0 ], "primeMeridianDeg": [ - 296.53, - -61.2572637, + 296.5309131499458, + -61.2572638, 0 ], "terms": [ @@ -2124,8 +2171,8 @@ 0 ], "primeMeridianDeg": [ - 93.38, - 320.7654228, + 91.53817645008237, + 320.7656245, 0 ], "terms": [ @@ -2190,8 +2237,8 @@ 0 ], "primeMeridianDeg": [ - 122.695, - 56.3625225, + 122.70869724996541, + 56.362521, 0 ] } diff --git a/tools/etl/build.ts b/tools/etl/build.ts index bea1968..f549c0f 100644 --- a/tools/etl/build.ts +++ b/tools/etl/build.ts @@ -11,6 +11,7 @@ import { fetchDeepSky } from './fetchDeepSky'; import { fetchExoplanets } from './fetchExoplanets'; import { fetchSolarSystem, FREELY_SPINNING_MOONS, offsetFromTrackDeg } from './fetchSolarSystem'; import { TrackPoint } from './lib/horizons'; +import { subPlanetLongitudeDeg } from './lib/locked-spin'; import { BYTES_PER_STAR_META, BYTES_PER_STAR_POSITION, decodeStarCatalog, encodeStarCatalog } from '../../src/app/shared/models/star-catalog'; import { fetchStars } from './fetchStars'; import { describeSources } from './sources/registry'; @@ -223,37 +224,29 @@ 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. + * sampled every 135 days over the clock's AD 1 to 3000. Every locked moon's W turns at its orbit's + * own rate (see `lockedToOrbit`); at the IAU's own rates, and sampled only from 1950 to 2100, this + * let Proteus turn its far side to Neptune at AD 1 (146 degrees), Iapetus 87 degrees, Mimas 52 and + * Miranda 23, on dates the clock offers. * - * 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: its drawn - * face runs from 2.5 to 10.15 degrees, about 6.3 off on average because the IAU's W and JPL's mean - * longitude disagree, drifting 3.3 over the span because W turns 6.0e-5 degrees a day faster than - * the row's n, and swung 2.3 either way (2e) by its eccentricity. None of that is Mimas: its - * measured physical libration is 0.84 degrees (Tajeddine et al. 2014, Science 346, 322), and W - * carries none; Horizons, on the same W against its integrated orbit, runs from -2.7 to 12.7 - * degrees over 1950-2100 with the 71-year S5 term the orbit here cancels. 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. + * Measured on this catalogue: at most 5.36 degrees (Titan) but for three. The Moon 7.62, at AD 1: + * its longitude swings 6.3 either way with its eccentricity, Horizons' too, and W's quadratic, the + * tidal slowing its orbit here does not carry, adds 0.75 by then. Mimas 8.94: about 6.3 off on + * average because the IAU's W and JPL's mean longitude disagree, and swung 2.3 either way (2e) by + * its eccentricity. None of that is Mimas: its measured physical libration is 0.84 degrees + * (Tajeddine et al. 2014, Science 346, 322), and W carries none; Horizons, on the same W against its + * integrated orbit, runs from -2.7 to 12.7 degrees over 1950-2100 with the 71-year S5 term the + * orbit here cancels. Iapetus 15.95, whose row sits 9.4 degrees behind Horizons. 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 = { 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; -} +const SUB_PLANET_CEILINGS_DEG: Record = { moon: 8, mimas: 9.5, iapetus: 16.5 }; +/** The clock's window, AD 1 to 3000 (`CLOCK_WINDOW` in `time.store.ts`), as Julian dates. */ +const CLOCK_START_JD = Date.parse('0001-01-01T00:00Z') / 86400000 + 2440587.5; +const CLOCK_END_JD = Date.parse('3000-01-01T00:00Z') / 86400000 + 2440587.5; +const LOCK_DATES_JD = Array.from({ length: Math.floor((CLOCK_END_JD - CLOCK_START_JD) / 135) + 1 }, (_, index) => CLOCK_START_JD + index * 135); 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)); @@ -368,10 +361,10 @@ function validateBodies(bodies: BodyRecord[], horizonsOrbits: Map Math.abs(subPlanetLongitudeDeg(body, jd)))); + const worst = Math.max(...LOCK_DATES_JD.map((jd) => Math.abs(subPlanetLongitudeDeg(body, rotation!, jd)))); assertCondition( 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.` + `Moon ${body.id} turns its face up to ${worst.toFixed(2)} degrees from its planet between AD 1 and 3000 (at most ${ceiling} expected) — its orbit and its W disagree.` ); spins.push(`${body.id} faces ${worst.toFixed(2)}`); } diff --git a/tools/etl/fetchSolarSystem.ts b/tools/etl/fetchSolarSystem.ts index aceec1e..cfe4239 100644 --- a/tools/etl/fetchSolarSystem.ts +++ b/tools/etl/fetchSolarSystem.ts @@ -9,6 +9,7 @@ import { MeanOrbit, parsePlanetMeanElements, parseSatelliteMeanElements, parseSm import { fetchPlanetMeanElementsText, fetchSatelliteMeanElementsHtml, fetchSmallBodyAnswer } from './lib/mean-elements'; import { MIN_PERIODIC_TERM_DEG, orbitalTermsOfPrimeMeridian, parsePckRotationalElements, SUN_ROTATIONAL_ELEMENTS } from '../../src/app/shared/astro/rotational-elements'; import { fetchPckText } from './lib/pck'; +import { lockedToOrbit } from './lib/locked-spin'; import { dataPath, ensureDataDir } from './lib/paths'; const HOURS_PER_DAY = 24; @@ -54,6 +55,8 @@ interface BodySpec { * given. See `orbitalTermsOfPrimeMeridian`. */ orbitFromW?: { angleRateDegPerCentury?: number }; + /** A locked moon whose pole is carried round with its orbit's, as Iapetus's; see `lockedToOrbit`. */ + poleFollowsOrbit?: boolean; /** * Days between the Horizons positions the orbit is checked against from 1950 to 2100; 2 unless * the error changes faster than that. Nereid, at an eccentricity of 0.75, sweeps through its @@ -131,7 +134,7 @@ const BODY_SPECS: BodySpec[] = [ // Hyperion tumbles ("Rotational period = Chaotic") and Phoebe, captured, turns in 9.27 hours. // Hyperion's eccentricity is 0.105 in JPL's current table (ssd.jpl.nasa.gov/sats/elem, SAT441). { id: 'hyperion', name: 'Hyperion', kind: 'moon', horizonsCommand: '607', center: '500@699', parentBodyId: 'saturn', spinsFreely: true, measuredEccentricity: 0.105, trackStepDays: 1 }, - { id: 'iapetus', name: 'Iapetus', kind: 'moon', horizonsCommand: '608', center: '500@699', parentBodyId: 'saturn' }, + { id: 'iapetus', name: 'Iapetus', kind: 'moon', horizonsCommand: '608', center: '500@699', parentBodyId: 'saturn', poleFollowsOrbit: true }, // Phoebe's row gives a mean motion of 0.6569114 degrees a day, a 548.02-day year, where its // Horizons page and JPL's current table (SAT441) give 550.30: the table's own note warns that // its source misstated the mean motions of retrograde moons. On the row's figure Phoebe was @@ -245,10 +248,11 @@ export async function fetchSolarSystem(): Promise<{ bodies: BodyRecord[]; horizo // 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. 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) + // turned by its IAU W, at this same rate (see `lockedToOrbit`), and `build.ts` checks that W and + // the orbit keep its face to its planet from AD 1 to 3000; this day is what the renderer would + // turn a moon without W by. + const locked = spec.kind === 'moon' && !spec.spinsFreely; + const rotationPeriodHours = result.tidallyLocked || locked ? (360 / mean.rates.meanMotionDegPerDay) * HOURS_PER_DAY : (spec.rotationPeriodHours ?? (smallBody ? smallBody.rotationPeriodHours : result.rotationPeriodHours)); const parentGm = spec.barycentric && spec.parentBodyId ? gmById.get(spec.parentBodyId) : undefined; @@ -274,7 +278,7 @@ export async function fetchSolarSystem(): Promise<{ bodies: BodyRecord[]; horizo ...(parentGm !== undefined ? { massRatio: result.gmKm3PerS2! / parentGm } : {}), ...(rotationPeriodHours !== undefined ? { rotationPeriodHours } : {}), ...((result.obliquityDeg ?? spec.obliquityDeg) !== undefined ? { obliquityDeg: result.obliquityDeg ?? spec.obliquityDeg } : {}), - ...(rotation ? { rotationalElements: rotation.elements } : {}) + ...(rotation ? { rotationalElements: locked ? lockedToOrbit(rotation.elements, mean, spec.name, spec.poleFollowsOrbit) : rotation.elements } : {}) }); } diff --git a/tools/etl/lib/locked-spin.ts b/tools/etl/lib/locked-spin.ts new file mode 100644 index 0000000..2f85ac4 --- /dev/null +++ b/tools/etl/lib/locked-spin.ts @@ -0,0 +1,141 @@ +import { eclipticToEquatorial, laplacePlaneToEquatorial } from '../../../src/app/shared/astro/coordinates'; +import { meanElementsAt, positionAtEpoch } from '../../../src/app/shared/astro/kepler'; +import { MeanOrbit } from '../../../src/app/shared/astro/mean-elements'; +import { orientationAt } from '../../../src/app/shared/astro/rotational-elements'; +import { BodyRecord, RotationalElements } from '../../../src/app/shared/models/body.model'; + +const J2000_JD = 2451545; +const DAYS_PER_JULIAN_CENTURY = 36525; +const DEG_TO_RAD = Math.PI / 180; + +/** 2025-01-01, the date the ETL asks Horizons about: where a locked moon's W is left as the IAU has it. */ +export const PRESENT_JD = 2460676.5; + +/** + * How far the IAU's W rate for a locked moon may be from the mean motion its orbit is drawn at, as + * a fraction of it, before it is taken for the orbit's. Measured: at most 3.4e-6 (Iapetus; Proteus + * 6.3e-7). What this catches is a rate read for the wrong body: Oberon's for Titania's is 55 per cent out. + */ +const MAX_LOCKED_RATE_OFFSET = 1e-5; + +/** Harmonics of the node's angle that carry a pole round its orbit's; see {@link lockedToOrbit}. */ +const POLE_HARMONICS = 5; + +/** The planet's east longitude on a moon's IAU body-fixed frame, from the moon's mean place, at a TDB date. */ +export function subPlanetLongitudeDeg(body: Pick, elements: RotationalElements, 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(elements, 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; +} + +/** + * A locked moon's IAU elements, turned at the rate its orbit is drawn at, so it keeps its face to + * its planet over the clock's AD 1 to 3000 and not only near the present its W was fitted to. + * + * The report gives a locked moon's W the mean motion of whichever orbit its authors had, and JPL's + * table has another: Proteus's W turns 6.3e-7 of its rate slower than its row, which turned its far + * side to Neptune at AD 1 (146 degrees), Mimas's 1.6e-7 faster (52 at AD 1) and Miranda's (23). + * W's rate is set to the orbit's here, its constant moved so W is unchanged at {@link PRESENT_JD}, + * and the pole and every periodic term are the IAU's. Measured over AD 1-3000: Proteus 2.7 degrees, + * Mimas 8.9, Miranda 2.8, Ariel 1.0. A W with a quadratic is left: Phobos's orbit already takes the + * quadratic from W (see `orbitalTermsOfPrimeMeridian`), and the Moon's, its tidal slowing, is 0.75 + * degrees at AD 1. + * + * `poleFollowsOrbit` is for Iapetus, whose IAU pole moves 3.9 degrees a century in right ascension + * and 1.1 in declination: a straight line through its orbit normal's 3 439-year circle round the + * Laplace pole, 8.3 degrees across, which by AD 1 has run past the celestial pole (Dec 97.9) and 11 + * degrees off the orbit, and turned its face 87 degrees from Saturn. Its axis sits on its orbit normal + * (0.04 degrees apart today), as a moon in a Cassini state keeps it, so its pole is given the circle: the + * normal's right ascension as sines and declination as cosines of the node's angle and its first + * {@link POLE_HARMONICS} harmonics, the IAU's own form for a precessing pole, with its constants + * set so the pole is the IAU's at the present. W counts from where the equator crosses the ICRF + * equator, which swings as the pole goes round, so W takes sines of the same angles, fitted to hold + * the face where it is today. Measured over AD 1-3000: the axis within 0.74 degrees of the orbit + * normal (the IAU's line, 11.06), and the face within 16 of Saturn (87), which is what the row's own + * 9.4-degree lag and its eccentricity make it from 1950 to 2100 as well (15.9). + */ +export function lockedToOrbit(elements: RotationalElements, mean: Pick, name: string, poleFollowsOrbit = false): RotationalElements { + const [w0, w1, w2 = 0] = elements.primeMeridianDeg; + if (w2 !== 0) { + return elements; + } + const n = Math.sign(w1) * mean.rates.meanMotionDegPerDay; + if (Math.abs(w1 / n - 1) > MAX_LOCKED_RATE_OFFSET) { + throw new Error(`${name}'s IAU W turns at ${w1} degrees a day, ${Math.abs(w1 / n - 1).toExponential(2)} of its orbit's ${n}: not the rate of the orbit it keeps its face to.`); + } + const locked: RotationalElements = { ...elements, primeMeridianDeg: [w0 + (w1 - n) * (PRESENT_JD - J2000_JD), n, 0] }; + return poleFollowsOrbit ? poleRoundOrbit(locked, mean) : locked; +} + +function poleRoundOrbit(elements: RotationalElements, mean: Pick): RotationalElements { + if (!mean.laplacePole) { + throw new Error('A pole that follows its orbit is carried round the orbit\'s Laplace pole, and this orbit has none.'); + } + const laplacePole = mean.laplacePole; + const normalAt = (jd: number) => { + const { inclinationDeg, longitudeOfAscendingNodeDeg } = meanElementsAt(mean.orbit, mean.rates, jd); + const tilt = inclinationDeg * DEG_TO_RAD; + const node = longitudeOfAscendingNodeDeg * DEG_TO_RAD; + const normal = laplacePlaneToEquatorial({ x: Math.sin(tilt) * Math.sin(node), y: -Math.sin(tilt) * Math.cos(node), z: Math.cos(tilt) }, laplacePole); + return { raDeg: Math.atan2(normal.y, normal.x) / DEG_TO_RAD, decDeg: Math.asin(normal.z) / DEG_TO_RAD }; + }; + // The node's angle, T in centuries, turned so that 0 is where the normal is furthest north: the + // circle is then even in declination and odd in right ascension about it, as the form requires. + const nodeRate = mean.rates.longitudeOfAscendingNodeDegPerDay * DAYS_PER_JULIAN_CENTURY; + const nodeAtJ2000 = meanElementsAt(mean.orbit, mean.rates, J2000_JD).longitudeOfAscendingNodeDeg; + const jdAtAngle = (angleDeg: number, phaseDeg: number) => J2000_JD + ((angleDeg - phaseDeg - nodeAtJ2000) / nodeRate) * DAYS_PER_JULIAN_CENTURY; + let phase = 0; + let northmost = -Infinity; + for (let candidate = 0; candidate < 360; candidate += 0.01) { + const dec = normalAt(jdAtAngle(0, candidate)).decDeg; + if (dec > northmost) { + northmost = dec; + phase = candidate; + } + } + const centre = laplacePole; + const samples = 3600; + const ra = new Array(POLE_HARMONICS + 1).fill(0); + const dec = new Array(POLE_HARMONICS + 1).fill(0); + for (let sample = 0; sample < samples; sample++) { + const angle = (sample / samples) * 360; + const normal = normalAt(jdAtAngle(angle, phase)); + const raOffset = ((((normal.raDeg - centre.raDeg) % 360) + 540) % 360) - 180; + for (let k = 0; k <= POLE_HARMONICS; k++) { + ra[k] += (2 / samples) * raOffset * Math.sin(k * angle * DEG_TO_RAD); + dec[k] += ((k === 0 ? 1 : 2) / samples) * (normal.decDeg - centre.decDeg) * Math.cos(k * angle * DEG_TO_RAD); + } + } + const terms = Array.from({ length: POLE_HARMONICS }, (_, index) => { + const k = index + 1; + return { angleDeg: [k * (nodeAtJ2000 + phase), k * nodeRate], ra: ra[k], dec: dec[k], pm: 0 }; + }); + const round: RotationalElements = { ...elements, poleRaDeg: [centre.raDeg, 0, 0], poleDecDeg: [centre.decDeg + dec[0], 0, 0], terms: [...(elements.terms ?? []), ...terms] }; + // The IAU's pole at the present, exactly: the fitted circle's constants moved onto it. + const iau = orientationAt(elements, PRESENT_JD); + const fitted = orientationAt(round, PRESENT_JD); + round.poleRaDeg = [round.poleRaDeg[0] + iau.poleRaDeg - fitted.poleRaDeg, 0, 0]; + round.poleDecDeg = [round.poleDecDeg[0] + iau.poleDecDeg - fitted.poleDecDeg, 0, 0]; + // W's sines on the same angles, fitted over a turn of the node to what the face drifts by. + const pm = new Array(POLE_HARMONICS + 1).fill(0); + const wSamples = 36000; + for (let sample = 0; sample < wSamples; sample++) { + const angle = (sample / wSamples) * 360; + const drift = subPlanetLongitudeDeg(mean, round, jdAtAngle(angle, phase)); + for (let k = 1; k <= POLE_HARMONICS; k++) { + pm[k] += (2 / wSamples) * drift * Math.sin(k * angle * DEG_TO_RAD); + } + } + const ownTerms = elements.terms?.length ?? 0; + const turned: RotationalElements = { ...round, terms: round.terms!.map((term, index) => (index < ownTerms ? term : { ...term, pm: pm[index - ownTerms + 1] })) }; + // And W the IAU's at the present. + const shift = orientationAt(turned, PRESENT_JD).primeMeridianDeg - orientationAt(elements, PRESENT_JD).primeMeridianDeg; + turned.primeMeridianDeg = [turned.primeMeridianDeg[0] - shift, turned.primeMeridianDeg[1], 0]; + return turned; +}