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 051bd5b..2c228ec 100644 --- a/src/app/features/galaxy-system/system-orbits-renderer.spec.ts +++ b/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 = { @@ -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 { diff --git a/src/app/shared/astro/rotational-elements.spec.ts b/src/app/shared/astro/rotational-elements.spec.ts index 7a2dd3d..5900f47 100644 --- a/src/app/shared/astro/rotational-elements.spec.ts +++ b/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(); + }); +}); diff --git a/src/app/shared/astro/rotational-elements.ts b/src/app/shared/astro/rotational-elements.ts index c0da58a..54da1de 100644 --- a/src/app/shared/astro/rotational-elements.ts +++ b/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 + }; +} diff --git a/src/assets/data/bodies.json b/src/assets/data/bodies.json index 16593db..adc1bed 100644 --- a/src/assets/data/bodies.json +++ b/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, diff --git a/tools/etl/build.ts b/tools/etl/build.ts index 9eaa97d..7af145b 100644 --- a/tools/etl/build.ts +++ b/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 = { mimas: 46, hyperion: 21, iapetus: 11, nereid: 3 }; +const MOON_OFFSET_CEILINGS_DEG: Record = { 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 = { 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 = { 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 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 — diff --git a/tools/etl/fetchSolarSystem.ts b/tools/etl/fetchSolarSystem.ts index aa11574..e09f3a3 100644 --- a/tools/etl/fetchSolarSystem.ts +++ b/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,