diff --git a/README.md b/README.md index 0890633..b44be02 100644 --- a/README.md +++ b/README.md @@ -56,7 +56,10 @@ in it is measured and what is not. cutting to a new scene. The Sun gets the real solar-system bodies, moving on JPL's mean orbital elements — Standish's for the planets, JPL SSD's satellite table for the moons, the Small-Body Database for Ceres, Eris, Haumea and Makemake — and turned by the IAU's rotational elements, Earth -by the IERS Earth Rotation Angle; other +by the IERS Earth Rotation Angle; a tidally locked moon's prime meridian (but the Moon's and +Phobos's) turns at its JPL mean motion and its pole's terms on its node at its JPL node rate, both +re-phased to the IAU's values on 2025-01-01, and Iapetus's pole follows its orbit normal +(`lockedToOrbit`), so each keeps its face to its planet from AD 1 to 3000; other stars get their confirmed exoplanets. Orbits are drawn as ellipses and bodies are propagated along them by a Kepler solver to the date on the map's clock. Under them, a dashed grid marks out round distances in AU — 5 AU rings for the solar system, 0.01 AU rings for TRAPPIST-1 — with a @@ -330,7 +333,7 @@ plugin's own files are kept so it can be listed from a marketplace of its own la ## Data credits Star catalogue: [HYG database](https://github.com/astronexus/HYG-Database) (Hipparcos, Yale -Bright Star, Gliese) — 68 388 stars within 250 pc. Solar-system orbits: JPL approximate planetary mean elements (Standish), JPL SSD satellite mean elements and the JPL Small-Body Database; rotation: the IAU WGCCRE 2015 report via NAIF's pck00011, and for Earth the IERS Conventions 2010; physical data, and the positions the orbits are checked against: NASA/JPL Horizons. Exoplanets: NASA Exoplanet +Bright Star, Gliese) — 68 388 stars within 250 pc. Solar-system orbits: JPL approximate planetary mean elements (Standish), JPL SSD satellite mean elements and the JPL Small-Body Database; rotation: the IAU WGCCRE 2015 report via NAIF's pck00011, with a locked moon's W and node terms re-rated to its JPL mean elements and Iapetus's pole carried round its orbit normal, and for Earth the IERS Conventions 2010; physical data, and the positions the orbits are checked against: NASA/JPL Horizons. Exoplanets: NASA Exoplanet Archive. Deep-sky objects: [OpenNGC](https://github.com/mattiaverga/OpenNGC). Body and skybox imagery: NASA/JPL/USGS public domain and Solar System Scope (CC BY 4.0) — per-file provenance is recorded in `src/assets/textures/README.md`. 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 333be1e..5dc6f5a 100644 --- a/src/app/features/galaxy-system/system-orbits-renderer.spec.ts +++ b/src/app/features/galaxy-system/system-orbits-renderer.spec.ts @@ -923,7 +923,7 @@ describe('locked moons across the clock’s window', () => { } 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 + // Measured: Proteus 2.6 degrees at most over AD 1-3000, Miranda 2.4, 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]) { @@ -934,4 +934,21 @@ describe('locked moons across the clock’s window', () => { expect(Math.abs(facingPlanet('iapetus'))).toBeLessThan(16.5); } }); + + it('keeps the axes of Miranda, Mimas and Iapetus on their drawn orbits’ normals, as a Cassini state holds them, at AD 1, today and AD 3000', () => { + // Measured over AD 1-3000: Miranda 0.42 degrees at most, Mimas 0.47, Iapetus 0.74. With their + // poles going round at the IAU's node rates, Miranda was 7.6 off at AD 1 and Mimas 2.6; with + // Iapetus's pole on its Laplace pole, 8.3 off at every date. + for (const jd of [1721425.5, 2460676.5, 2816787.4]) { + renderer.update(jd); + for (const id of ['miranda', 'mimas', 'iapetus']) { + const moon = renderer.members.find((member) => member.id === id)!.marker; + const line = moon.parent!.children.find((child) => child.name === 'orbit-line')!; + const axis = new THREE.Vector3(0, 1, 0).applyQuaternion(moon.quaternion); + const normal = new THREE.Vector3(0, 0, 1).applyQuaternion(line.quaternion); + // A line, not a direction: Miranda turns backwards against the IAU's pole. + expect((Math.acos(Math.min(1, Math.abs(axis.dot(normal)))) * 180) / Math.PI).toBeLessThan(1); + } + } + }); }); diff --git a/src/app/shared/models/body.model.ts b/src/app/shared/models/body.model.ts index aa717c2..50611a6 100644 --- a/src/app/shared/models/body.model.ts +++ b/src/app/shared/models/body.model.ts @@ -105,8 +105,9 @@ export interface BodyRecord { /** * 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, 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 + * moon's W, and its pole's terms on its node, turn at its drawn orbit's rates and Iapetus's pole + * goes round with its orbit's, so they keep their faces to their planets, and their poles round + * their orbits', 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. diff --git a/src/assets/data/bodies.json b/src/assets/data/bodies.json index 978b3d1..b6a7305 100644 --- a/src/assets/data/bodies.json +++ b/src/assets/data/bodies.json @@ -821,9 +821,8 @@ "terms": [ { "angleDeg": [ - 121.46893664, - 660.22803474, - 0 + 121.49945940952146, + 660.1059470044942 ], "ra": 3.09217726, "dec": 1.83936004, @@ -831,9 +830,8 @@ }, { "angleDeg": [ - 231.05028581, - 660.9912354, - 0 + 231.2716139683453, + 660.1059470044942 ], "ra": 0.22980637, "dec": 0.1432532, @@ -841,9 +839,8 @@ }, { "angleDeg": [ - 251.37314025, - 1320.50145245, - 0 + 251.44553184217244, + 1320.2118940089883 ], "ra": 0.06418655, "dec": 0.01911409, @@ -918,8 +915,8 @@ "terms": [ { "angleDeg": [ - 283.9, - 4850.7 + 283.6369874084692, + 4851.752021563342 ], "ra": 0.094, "dec": 0.04, @@ -983,8 +980,8 @@ "terms": [ { "angleDeg": [ - 355.8, - 1191.3 + 355.4537739825443, + 1192.6848661542538 ], "ra": 1.086, "dec": 0.468, @@ -1066,8 +1063,8 @@ }, { "angleDeg": [ - 119.9, - 262.1 + 117.57926275587673, + 271.38269483015966 ], "ra": 0.431, "dec": 0.186, @@ -1211,8 +1208,8 @@ "terms": [ { "angleDeg": [ - 177.4, - -36505.5 + 178.81408536763, + -36511.15618661257 ], "ra": 13.56, "dec": -1.53, @@ -1327,8 +1324,8 @@ "terms": [ { "angleDeg": [ - 300, - -7225.9 + 300.02841306210905, + -7226.013649136892 ], "ra": 9.66, "dec": -1.09, @@ -1437,8 +1434,8 @@ "terms": [ { "angleDeg": [ - 345.2, - -1016.3 + 342.2970571615749, + -1004.6885465505693 ], "ra": 3.1, "dec": -0.35, @@ -1703,8 +1700,8 @@ "terms": [ { "angleDeg": [ - 102.23, - -2024.22 + 103.87516350424971, + -2030.8004738534437 ], "ra": 4.41, "dec": 4.25, @@ -1721,8 +1718,8 @@ }, { "angleDeg": [ - 204.46, - -4048.44 + 207.75032700849943, + -4061.6009477068874 ], "ra": -0.04, "dec": -0.02, @@ -2028,8 +2025,8 @@ "terms": [ { "angleDeg": [ - 177.85, - 52.316 + 177.83706224270992, + 52.367749612333185 ], "ra": -32.35, "dec": 22.55, @@ -2037,8 +2034,8 @@ }, { "angleDeg": [ - 355.7, - 104.632 + 355.67412448541984, + 104.73549922466637 ], "ra": -6.28, "dec": 2.1, @@ -2046,8 +2043,8 @@ }, { "angleDeg": [ - 533.55, - 156.948 + 533.5111867281297, + 157.10324883699957 ], "ra": -2.08, "dec": 0.55, @@ -2055,8 +2052,8 @@ }, { "angleDeg": [ - 711.4, - 209.264 + 711.3482489708397, + 209.47099844933274 ], "ra": -0.74, "dec": 0.16, @@ -2064,8 +2061,8 @@ }, { "angleDeg": [ - 889.25, - 261.58 + 889.1853112135495, + 261.8387480616659 ], "ra": -0.28, "dec": 0.05, @@ -2073,8 +2070,8 @@ }, { "angleDeg": [ - 1067.1, - 313.896 + 1067.0223734562594, + 314.20649767399914 ], "ra": -0.11, "dec": 0.02, @@ -2082,8 +2079,8 @@ }, { "angleDeg": [ - 1244.95, - 366.212 + 1244.8594356989695, + 366.5742472863323 ], "ra": -0.07, "dec": 0.01, @@ -2091,8 +2088,8 @@ }, { "angleDeg": [ - 1422.8, - 418.528 + 1422.6964979416794, + 418.9419968986655 ], "ra": -0.02, "dec": 0, @@ -2100,8 +2097,8 @@ }, { "angleDeg": [ - 1600.65, - 470.844 + 1600.5335601843892, + 471.30974651099865 ], "ra": -0.01, "dec": 0, @@ -2193,8 +2190,8 @@ }, { "angleDeg": [ - 142.61, - 2824.6 + 141.60196187337888, + 2828.6320421151877 ], "ra": -0.05, "dec": -0.04, diff --git a/tools/etl/build.ts b/tools/etl/build.ts index 5cf738c..e6f20ba 100644 --- a/tools/etl/build.ts +++ b/tools/etl/build.ts @@ -1,6 +1,6 @@ import { statSync } from 'node:fs'; -import { BodyRecord, OrbitalElements } from '../../src/app/shared/models/body.model'; +import { BodyRecord, OrbitalElements, RotationalElements } from '../../src/app/shared/models/body.model'; import { eclipticToEquatorial, laplacePlaneToEquatorial, raDecToUnitVector } from '../../src/app/shared/astro/coordinates'; import { meanElementsAt, positionAtEpoch } from '../../src/app/shared/astro/kepler'; import { orientationAt } from '../../src/app/shared/astro/rotational-elements'; @@ -257,7 +257,7 @@ const MAX_OBLIQUITY_OFFSET_DEG = 0.1; * * 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 + * tidal slowing its orbit here does not carry, adds 0.75 by then. Mimas 8.89: 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 @@ -269,6 +269,22 @@ const MAX_OBLIQUITY_OFFSET_DEG = 0.1; */ const MAX_SUB_PLANET_LONGITUDE_DEG = 7; const SUB_PLANET_CEILINGS_DEG: Record = { moon: 8, mimas: 9.5, iapetus: 16.5 }; +/** + * How far a locked moon's spin axis may lean from the normal of the orbit it is drawn going round, + * over the same dates. A locked moon sits in a Cassini state, its axis on its orbit normal as the + * node carries both round the Laplace pole, and the IAU's pole goes round on a term of the node's + * angle; at the rate the IAU's source had for it and not the drawn orbit's, Miranda's axis was 7.89 + * degrees off at AD 9 and Mimas's 2.63, and an Iapetus pole left on the Laplace pole is 8.30 off at + * every date (see `lockedToOrbit`). + * + * Measured on this catalogue: at most 0.97 degrees (Tethys, whose IAU pole sits 0.69 from its orbit + * normal today; Titan 0.94, whose pole the IAU holds still while its node turns in 705 years) but + * for four. The Moon 6.98, its real 6.7-degree tilt to its orbit. Phobos 1.81 and Deimos 1.74, and + * Proteus 1.09: their IAU poles nod with Mars's and Neptune's precessing poles, the Laplace poles + * their orbits are drawn round are fixed. + */ +const MAX_AXIS_FROM_ORBIT_DEG = 1; +const AXIS_FROM_ORBIT_CEILINGS_DEG: Record = { moon: 7.1, phobos: 2, deimos: 2, proteus: 1.2 }; /** 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; @@ -279,6 +295,19 @@ function angleBetweenDeg(a: { x: number; y: number; z: number }, b: { x: number; return (Math.acos(Math.min(1, Math.max(-1, cosine))) * 180) / Math.PI; } +/** Degrees between a body's spin axis — its IAU pole, turned over where W runs backwards — and the normal of the orbit it is drawn going round, at a TDB date. */ +function axisFromOrbitDeg(body: BodyRecord, rotation: RotationalElements, jd: number): number { + const pole = orientationAt(rotation, jd); + const pointing = raDecToUnitVector(pole.poleRaDeg / 15, pole.poleDecDeg); + const sense = Math.sign(rotation.primeMeridianDeg[1]); + const axis = { x: sense * pointing.x, y: sense * pointing.y, z: sense * pointing.z }; + const { inclinationDeg, longitudeOfAscendingNodeDeg } = meanElementsAt(body.orbit, body.rates, jd); + const tilt = inclinationDeg * DEG_TO_RAD; + const node = longitudeOfAscendingNodeDeg * DEG_TO_RAD; + const normal = { x: Math.sin(tilt) * Math.sin(node), y: -Math.sin(tilt) * Math.cos(node), z: Math.cos(tilt) }; + return angleBetweenDeg(axis, body.laplacePole ? laplacePlaneToEquatorial(normal, body.laplacePole) : eclipticToEquatorial(normal)); +} + function validateBodies(bodies: BodyRecord[], horizonsOrbits: Map, horizonsTracks: Map): void { assertCondition(bodies.length > 0, 'No solar-system bodies were produced.'); @@ -364,14 +393,7 @@ function validateBodies(bodies: BodyRecord[], horizonsOrbits: Map axisFromOrbitDeg(body, rotation!, jd))); + assertCondition( + worstAxis <= axisCeiling, + `Moon ${body.id}'s spin axis leans up to ${worstAxis.toFixed(2)} degrees from its orbit's normal between AD 1 and 3000 (at most ${axisCeiling} expected) — its pole does not go round with its node.` + ); + spins.push(`${body.id} axis ${worstAxis.toFixed(2)}`); } if (body.massRatio !== undefined) { // The pair's barycentre, which the planet's elements place, must lie outside the planet — diff --git a/tools/etl/lib/locked-spin.ts b/tools/etl/lib/locked-spin.ts index 2f85ac4..82de187 100644 --- a/tools/etl/lib/locked-spin.ts +++ b/tools/etl/lib/locked-spin.ts @@ -21,6 +21,16 @@ 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; +/** + * How far a periodic term's angle may turn from a multiple of the node's rate, as a fraction of it, + * and still be taken for the node's angle as the IAU's source had it. Measured: at most 3.4e-2 + * (Ganymede's J5), then Rhea's R4 1.2e-2 and Miranda's U11 3.2e-3; the nearest that is not a node is + * a term of Miranda's W alone, 6.0e-2 from three times it. Multiples go up to the ninth, the most the + * report takes (Triton's N7); past that, Umbriel's W has a term 1.0e-2 from ten times its node's. + */ +const MAX_NODE_RATE_OFFSET = 0.05; +const MAX_NODE_HARMONIC = 9; + /** 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)); @@ -41,15 +51,25 @@ export function subPlanetLongitudeDeg(body: Pick 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] }; + const nodeRate = mean.rates.longitudeOfAscendingNodeDegPerDay * DAYS_PER_JULIAN_CENTURY; + const present = (PRESENT_JD - J2000_JD) / DAYS_PER_JULIAN_CENTURY; + const terms = elements.terms?.map((term) => { + const [constant, rate, quadratic = 0] = term.angleDeg; + const k = Math.round(rate / nodeRate); + if (quadratic !== 0 || k === 0 || Math.abs(k) > MAX_NODE_HARMONIC || Math.abs(rate / (k * nodeRate) - 1) > MAX_NODE_RATE_OFFSET) { + return term; + } + return { ...term, angleDeg: [constant + (rate - k * nodeRate) * present, k * nodeRate] }; + }); + const locked: RotationalElements = { ...elements, primeMeridianDeg: [w0 + (w1 - n) * (PRESENT_JD - J2000_JD), n, 0], ...(terms ? { terms } : {}) }; return poleFollowsOrbit ? poleRoundOrbit(locked, mean) : locked; }