diff --git a/src/app/shared/astro/rotational-elements.spec.ts b/src/app/shared/astro/rotational-elements.spec.ts new file mode 100644 index 0000000..7a2dd3d --- /dev/null +++ b/src/app/shared/astro/rotational-elements.spec.ts @@ -0,0 +1,133 @@ +import { describe, expect, it } from 'vitest'; + +import { orientationAt, parsePckRotationalElements } from './rotational-elements'; + +// Excerpts of pck00011.tpc as NAIF publishes it: prose, then data blocks. +const KERNEL = String.raw`KPL/PCK + + The portion of the file preceding the first data block is treated + as a comment. + +\begindata + + BODY399_POLE_RA = ( 0. -0.641 0. ) + BODY399_POLE_DEC = ( 90. -0.557 0. ) + BODY399_PM = ( 190.147 360.9856235 0. ) + +\begintext + + A data block starts with the \begindata token only when that token + sits on a line by itself, so BODY399_PM = ( 1 2 3 ) here is prose. + +\begindata + + BODY301_POLE_RA = ( 269.9949 0.0031 0. ) + BODY301_POLE_DEC = ( 66.5392 0.0130 0. ) + BODY301_PM = ( 38.3213 13.17635815 -1.4D-12 ) + + BODY899_POLE_RA = ( 299.36 0. 0. ) + BODY899_POLE_DEC = ( 43.46 0. 0. ) + BODY899_PM = ( 249.978 541.1397757 0. ) + BODY899_NUT_PREC_RA = ( 0.70 0. 0. 0. 0. 0. 0. 0. ) + BODY899_NUT_PREC_DEC = ( -0.51 0. 0. 0. 0. 0. 0. 0. ) + BODY899_NUT_PREC_PM = ( -0.48 0. 0. 0. 0. 0. 0. 0. ) + BODY8_NUT_PREC_ANGLES = ( 357.85 52.316 + 323.92 62606.6 ) + + BODY199_POLE_RA = ( 281.0103 -0.0328 0. ) + BODY199_POLE_DEC = ( 61.4155 -0.0049 0. ) + BODY199_PM = ( 329.5988 6.1385108 0. ) + BODY199_NUT_PREC_RA = ( 0. 0. ) + BODY199_NUT_PREC_DEC = ( 0. 0. ) + BODY199_NUT_PREC_PM = ( 0.01067257 + -0.00112309 ) + BODY1_NUT_PREC_ANGLES = ( 174.7910857 0.14947253587500003E+06 + 349.5821714 0.29894507175000006E+06 ) + + BODY401_POLE_RA = ( 317.67071657 -0.10844326 0. ) + BODY401_POLE_DEC = ( 52.88627266 -0.06134706 0. ) + BODY401_PM = ( 35.18774440 1128.84475928 + 9.536137031212154e-09 ) + BODY401_NUT_PREC_RA = ( -1.78428399 ) + BODY401_NUT_PREC_DEC = ( -1.07516537 ) + BODY401_NUT_PREC_PM = ( 1.42421769 + -1.143 ) + BODY4_MAX_PHASE_DEGREE = 2 + BODY4_NUT_PREC_ANGLES = ( + 190.72646643 15917.10818695 0 + 189.63271560 41215158.18420050 12.711923222 ) + +\begintext +`; + +describe('parsePckRotationalElements', () => { + it('reads a pole and a prime meridian from the data blocks, not from the prose around them', () => { + expect(parsePckRotationalElements(KERNEL, 399)).toEqual({ + elements: { poleRaDeg: [0, -0.641, 0], poleDecDeg: [90, -0.557, 0], primeMeridianDeg: [190.147, 360.9856235, 0] }, + skippedDeg: [] + }); + }); + + it('reads the exponent the Fortran way, as the Moon’s quadratic is written', () => { + expect(parsePckRotationalElements(KERNEL, 301)!.elements.primeMeridianDeg).toEqual([38.3213, 13.17635815, -1.4e-12]); + }); + + it('pairs each periodic term with its system’s angle', () => { + expect(parsePckRotationalElements(KERNEL, 899)!.elements.terms).toEqual([{ angleDeg: [357.85, 52.316], ra: 0.7, dec: -0.51, pm: -0.48 }]); + }); + + it('reads angles to the degree the system states, Phobos’s quadratic among them', () => { + expect(parsePckRotationalElements(KERNEL, 401)!.elements.terms).toEqual([ + { angleDeg: [190.72646643, 15917.10818695, 0], ra: -1.78428399, dec: -1.07516537, pm: 1.42421769 }, + { angleDeg: [189.6327156, 41215158.1842005, 12.711923222], ra: 0, dec: 0, pm: -1.143 } + ]); + }); + + it('leaves out a term under a hundredth of a degree, and says how large it was', () => { + const mercury = parsePckRotationalElements(KERNEL, 199)!; + + expect(mercury.elements.terms).toEqual([{ angleDeg: [174.7910857, 149472.53587500003], ra: 0, dec: 0, pm: 0.01067257 }]); + expect(mercury.skippedDeg).toEqual([0.00112309]); + }); + + it('gives nothing for a body the kernel has no model for', () => { + expect(parsePckRotationalElements(KERNEL, 802)).toBeUndefined(); + }); +}); + +describe('orientationAt', () => { + const J2000 = 2451545.0; + + it('turns the prime meridian at its rate per day and moves the pole at its rate per century', () => { + const earth = parsePckRotationalElements(KERNEL, 399)!.elements; + + const epoch = orientationAt(earth, J2000); + expect([epoch.poleRaDeg, epoch.poleDecDeg]).toEqual([0, 90]); + expect(epoch.primeMeridianDeg).toBeCloseTo(190.147, 9); + const century = orientationAt(earth, J2000 + 36525); + expect(century.poleRaDeg).toBeCloseTo(-0.641, 12); + expect(century.poleDecDeg).toBeCloseTo(90 - 0.557, 12); + expect(orientationAt(earth, J2000 + 1).primeMeridianDeg).toBeCloseTo(190.147 + 360.9856235 - 360, 9); + }); + + it('adds a term as a sine to the right ascension and the meridian and a cosine to the declination', () => { + const neptune = parsePckRotationalElements(KERNEL, 899)!.elements; + const days = 9000; + const angle = ((357.85 + (52.316 * days) / 36525) * Math.PI) / 180; + const drawn = orientationAt(neptune, J2000 + days); + + expect(drawn.poleRaDeg).toBeCloseTo(299.36 + 0.7 * Math.sin(angle), 12); + expect(drawn.poleDecDeg).toBeCloseTo(43.46 - 0.51 * Math.cos(angle), 12); + expect(drawn.primeMeridianDeg).toBeCloseTo((249.978 + 541.1397757 * days - 0.48 * Math.sin(angle)) % 360, 6); + }); + + it('carries the quadratic in the meridian and in the angle, which is how Phobos falls inward', () => { + const phobos = parsePckRotationalElements(KERNEL, 401)!.elements; + const days = 36525; + const first = (190.72646643 + 15917.10818695) * (Math.PI / 180); + const second = (189.6327156 + 41215158.1842005 + 12.711923222) * (Math.PI / 180); + const expected = 35.1877444 + 1128.84475928 * days + 9.536137031212154e-9 * days * days + 1.42421769 * Math.sin(first) - 1.143 * Math.sin(second); + + expect(orientationAt(phobos, J2000 + days).primeMeridianDeg).toBeCloseTo(((expected % 360) + 360) % 360, 5); + }); +}); diff --git a/src/app/shared/astro/rotational-elements.ts b/src/app/shared/astro/rotational-elements.ts new file mode 100644 index 0000000..c0da58a --- /dev/null +++ b/src/app/shared/astro/rotational-elements.ts @@ -0,0 +1,123 @@ +import { RotationalElements } from '../models/body.model'; + +/** + * Reads the IAU WGCCRE 2015 rotational elements from NAIF's text kernel `pck00011.tpc`, which the + * ETL fetches (see `tools/etl/lib/pck.ts`), and evaluates them at a date. + */ + +const J2000_JD = 2451545.0; +const DAYS_PER_JULIAN_CENTURY = 36525; +const DEG_TO_RAD = Math.PI / 180; + +/** + * The smallest periodic term kept, in degrees. A term turns the body, or tips its pole, by at most + * its amplitude, and the largest a body is ever drawn is Jupiter filling the screen at 641 px of + * radius, where 0.01 degrees moves a point on its surface by 0.11 px. In `pck00011.tpc` this + * leaves out 32 terms: Mercury's four smaller librations (0.0011 degrees and less), eight of the + * Moon's thirteen (0.0072 and less), the thirteen short-period terms of Mars's pole and meridian + * (0.00024 and less; its three 0.42-1.59 degree long-period ones stay), one of Phobos's (0.0063), + * Jupiter's five (0.0022 and less) and one of Europa's (0.009). Mimas's 44.85-degree libration, + * Triton's 32-degree precession and Miranda's 4.4 are kept, down to Triton's 0.01. + */ +export const MIN_PERIODIC_TERM_DEG = 0.01; + +/** + * Every `NAME = ( values )` assignment in the kernel's data blocks. A data block runs from a line + * holding only `\begindata` to one holding only `\begintext`; the kernel's own prose mentions both + * tokens mid-sentence, which is why they are only read alone on a line. Exponents are written + * the Fortran way, `-1.4D-12`. + */ +function pckVariables(text: string): Map { + const data = text + .split(/^\s*\\begindata\s*$/m) + .slice(1) + .map((block) => block.split(/^\s*\\begintext\s*$/m)[0]) + .join('\n'); + const variables = new Map(); + for (const [, name, value] of data.matchAll(/(\w+)\s*=\s*(\([^)]*\)|\S+)/g)) { + variables.set( + name, + value + .replace(/[()]/g, ' ') + .trim() + .split(/[\s,]+/) + .filter(Boolean) + .map((token) => Number(token.replace(/d/i, 'e'))) + ); + } + return variables; +} + +/** + * One body's elements, by its NAIF id: 399 for Earth, 301 for the Moon, 2000001 for Ceres. + * Undefined where the kernel has none. + * + * The periodic terms' angles belong to the planet's whole system, `BODY5_NUT_PREC_ANGLES` for + * Jupiter and its moons, each a polynomial in T whose degree `BODYn_MAX_PHASE_DEGREE` gives: 1 + * unless stated, 2 for Mars, where Phobos's angle carries the tidal acceleration that is drawing + * it in. A term is kept if any of its three amplitudes reaches {@link MIN_PERIODIC_TERM_DEG}; + * the largest amplitude of each term left out comes back in `skippedDeg`, for the ETL to say so. + */ +export function parsePckRotationalElements(text: string, naifId: number): { elements: RotationalElements; skippedDeg: number[] } | undefined { + const variables = pckVariables(text); + const poleRaDeg = variables.get(`BODY${naifId}_POLE_RA`); + const poleDecDeg = variables.get(`BODY${naifId}_POLE_DEC`); + const primeMeridianDeg = variables.get(`BODY${naifId}_PM`); + if (!poleRaDeg || !poleDecDeg || !primeMeridianDeg) { + return undefined; + } + if (![...poleRaDeg, ...poleDecDeg, ...primeMeridianDeg].every(Number.isFinite)) { + throw new Error(`Body ${naifId}'s pole or prime meridian did not parse.`); + } + + const ra = variables.get(`BODY${naifId}_NUT_PREC_RA`) ?? []; + const dec = variables.get(`BODY${naifId}_NUT_PREC_DEC`) ?? []; + const pm = variables.get(`BODY${naifId}_NUT_PREC_PM`) ?? []; + const system = naifId < 1000 ? Math.floor(naifId / 100) : undefined; + const angles = system === undefined ? [] : (variables.get(`BODY${system}_NUT_PREC_ANGLES`) ?? []); + const coefficients = (variables.get(`BODY${system}_MAX_PHASE_DEGREE`)?.[0] ?? 1) + 1; + + const terms: NonNullable = []; + const skippedDeg: number[] = []; + for (let index = 0; index < Math.max(ra.length, dec.length, pm.length); index++) { + const term = { ra: ra[index] ?? 0, dec: dec[index] ?? 0, pm: pm[index] ?? 0 }; + const largest = Math.max(Math.abs(term.ra), Math.abs(term.dec), Math.abs(term.pm)); + if (largest === 0) { + continue; + } + if (largest < MIN_PERIODIC_TERM_DEG) { + skippedDeg.push(largest); + continue; + } + const angleDeg = angles.slice(index * coefficients, (index + 1) * coefficients); + if (angleDeg.length !== coefficients || !angleDeg.every(Number.isFinite)) { + throw new Error(`Body ${naifId}'s periodic term ${index + 1} has no angle among BODY${system}_NUT_PREC_ANGLES.`); + } + terms.push({ angleDeg, ...term }); + } + + return { + elements: { poleRaDeg, poleDecDeg, primeMeridianDeg, ...(terms.length > 0 ? { terms } : {}) }, + skippedDeg + }; +} + +function polynomial(coefficients: readonly number[], x: number): number { + return (coefficients[0] ?? 0) + (coefficients[1] ?? 0) * x + (coefficients[2] ?? 0) * x * x; +} + +/** The pole's right ascension and declination and the prime meridian W, in degrees, at a TDB Julian date. */ +export function orientationAt(elements: RotationalElements, jdTdb: number): { poleRaDeg: number; poleDecDeg: number; primeMeridianDeg: number } { + const days = jdTdb - J2000_JD; + const centuries = days / DAYS_PER_JULIAN_CENTURY; + let poleRaDeg = polynomial(elements.poleRaDeg, centuries); + let poleDecDeg = polynomial(elements.poleDecDeg, centuries); + let primeMeridianDeg = polynomial(elements.primeMeridianDeg, days); + for (const term of elements.terms ?? []) { + const angle = polynomial(term.angleDeg, centuries) * DEG_TO_RAD; + poleRaDeg += term.ra * Math.sin(angle); + poleDecDeg += term.dec * Math.cos(angle); + primeMeridianDeg += term.pm * Math.sin(angle); + } + return { poleRaDeg, poleDecDeg, primeMeridianDeg: ((primeMeridianDeg % 360) + 360) % 360 }; +} diff --git a/src/app/shared/models/body.model.ts b/src/app/shared/models/body.model.ts index 2d554c7..d1a53ae 100644 --- a/src/app/shared/models/body.model.ts +++ b/src/app/shared/models/body.model.ts @@ -86,4 +86,34 @@ export interface BodyRecord { */ rotationPeriodHours?: number; 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. + */ + rotationalElements?: RotationalElements; +} + +/** + * The IAU's rotational elements for one body: polynomials in time, plus periodic terms. + * + * The pole's right ascension and declination are in degrees in the ICRF, `[c0, c1, c2]` for + * `c0 + c1 T + c2 T²`, T in Julian centuries from J2000.0 TDB. The prime meridian W is the angle + * along the body's equator, anticlockwise seen from above that pole, from where the equator rises + * through the ICRF equator to the body's longitude 0, `c0 + c1 d + c2 d²` with d in days. A + * negative rate turns the body clockwise about the pole the IAU names: Venus, Uranus and its + * moons, Triton. + */ +export interface RotationalElements { + poleRaDeg: number[]; + poleDecDeg: number[]; + primeMeridianDeg: number[]; + /** + * Each adds `ra sin θ` to the right ascension, `dec cos θ` to the declination and `pm sin θ` to + * W, θ being `angleDeg[0] + angleDeg[1] T + angleDeg[2] T²`. The smallest are left out; see + * `parsePckRotationalElements`. + */ + terms?: Array<{ angleDeg: number[]; ra: number; dec: number; pm: number }>; } diff --git a/src/assets/data/bodies.json b/src/assets/data/bodies.json index e929fd3..16593db 100644 --- a/src/assets/data/bodies.json +++ b/src/assets/data/bodies.json @@ -24,7 +24,35 @@ }, "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 1407.512239412851, - "obliquityDeg": 0.035166666666666666 + "obliquityDeg": 0.035166666666666666, + "rotationalElements": { + "poleRaDeg": [ + 281.0103, + -0.0328, + 0 + ], + "poleDecDeg": [ + 61.4155, + -0.0049, + 0 + ], + "primeMeridianDeg": [ + 329.5988, + 6.1385108, + 0 + ], + "terms": [ + { + "angleDeg": [ + 174.7910857, + 149472.53587500003 + ], + "ra": 0, + "dec": 0, + "pm": 0.01067257 + } + ] + } }, { "id": "venus", @@ -51,7 +79,24 @@ }, "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": -5832.539941165383, - "obliquityDeg": 177.3 + "obliquityDeg": 177.3, + "rotationalElements": { + "poleRaDeg": [ + 272.76, + 0, + 0 + ], + "poleDecDeg": [ + 67.16, + 0, + 0 + ], + "primeMeridianDeg": [ + 160.2, + -1.4813688, + 0 + ] + } }, { "id": "earth", @@ -78,7 +123,24 @@ }, "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 23.934472399219285, - "obliquityDeg": 23.4392911 + "obliquityDeg": 23.4392911, + "rotationalElements": { + "poleRaDeg": [ + 0, + -0.641, + 0 + ], + "poleDecDeg": [ + 90, + -0.557, + 0 + ], + "primeMeridianDeg": [ + 190.147, + 360.9856235, + 0 + ] + } }, { "id": "mars", @@ -105,7 +167,56 @@ }, "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 24.622955438662025, - "obliquityDeg": 25.19 + "obliquityDeg": 25.19, + "rotationalElements": { + "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 + } + ] + } }, { "id": "jupiter", @@ -138,7 +249,24 @@ }, "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 9.925102371306965, - "obliquityDeg": 3.13 + "obliquityDeg": 3.13, + "rotationalElements": { + "poleRaDeg": [ + 268.056595, + -0.006499, + 0 + ], + "poleDecDeg": [ + 64.495303, + 0.002413, + 0 + ], + "primeMeridianDeg": [ + 284.95, + 870.536, + 0 + ] + } }, { "id": "saturn", @@ -171,7 +299,24 @@ }, "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 10.656221583138441, - "obliquityDeg": 26.73 + "obliquityDeg": 26.73, + "rotationalElements": { + "poleRaDeg": [ + 40.589, + -0.036, + 0 + ], + "poleDecDeg": [ + 83.537, + -0.004, + 0 + ], + "primeMeridianDeg": [ + 38.9, + 810.7939024, + 0 + ] + } }, { "id": "uranus", @@ -204,7 +349,24 @@ }, "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": -17.24003330792427, - "obliquityDeg": 97.77 + "obliquityDeg": 97.77, + "rotationalElements": { + "poleRaDeg": [ + 257.311, + 0, + 0 + ], + "poleDecDeg": [ + -15.175, + 0, + 0 + ], + "primeMeridianDeg": [ + 203.81, + -501.1600928, + 0 + ] + } }, { "id": "neptune", @@ -237,7 +399,35 @@ }, "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 16.110037586020876, - "obliquityDeg": 28.32 + "obliquityDeg": 28.32, + "rotationalElements": { + "poleRaDeg": [ + 299.36, + 0, + 0 + ], + "poleDecDeg": [ + 43.46, + 0, + 0 + ], + "primeMeridianDeg": [ + 249.978, + 541.1397757, + 0 + ], + "terms": [ + { + "angleDeg": [ + 357.85, + 52.316 + ], + "ra": 0.7, + "dec": -0.51, + "pm": -0.48 + } + ] + } }, { "id": "pluto", @@ -270,7 +460,24 @@ }, "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 153.29335198, - "obliquityDeg": 119.6 + "obliquityDeg": 119.6, + "rotationalElements": { + "poleRaDeg": [ + 132.993, + 0, + 0 + ], + "poleDecDeg": [ + -6.163, + 0, + 0 + ], + "primeMeridianDeg": [ + 302.695, + 56.3625225, + 0 + ] + } }, { "id": "ceres", @@ -293,7 +500,24 @@ "argumentOfPeriapsisDegPerDay": 0 }, "orbitSource": "JPL SBDB osculating elements, epoch 2026 Jun 9", - "rotationPeriodHours": 9.07417 + "rotationPeriodHours": 9.07417, + "rotationalElements": { + "poleRaDeg": [ + 291.418, + 0, + 0 + ], + "poleDecDeg": [ + 66.764, + 0, + 0 + ], + "primeMeridianDeg": [ + 170.65, + 952.1532, + 0 + ] + } }, { "id": "eris", @@ -387,7 +611,71 @@ "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "earth", "rotationPeriodHours": 655.7198886065481, - "obliquityDeg": 6.67 + "obliquityDeg": 6.67, + "rotationalElements": { + "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 + } + ] + } }, { "id": "phobos", @@ -415,7 +703,66 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1950 Jan 1", "parentBodyId": "mars", - "rotationPeriodHours": 7.653842521027348 + "rotationPeriodHours": 7.653842521027348, + "rotationalElements": { + "poleRaDeg": [ + 317.67071657, + -0.10844326, + 0 + ], + "poleDecDeg": [ + 52.88627266, + -0.06134706, + 0 + ], + "primeMeridianDeg": [ + 35.1877444, + 1128.84475928, + 9.536137031212154e-9 + ], + "terms": [ + { + "angleDeg": [ + 190.72646643, + 15917.10818695, + 0 + ], + "ra": -1.78428399, + "dec": -1.07516537, + "pm": 1.42421769 + }, + { + "angleDeg": [ + 21.4689247, + 31834.27934054, + 0 + ], + "ra": 0.02212824, + "dec": 0.00668626, + "pm": -0.02273783 + }, + { + "angleDeg": [ + 332.86082793, + 19139.89694742, + 0 + ], + "ra": -0.01028251, + "dec": -0.0064874, + "pm": 0.00410711 + }, + { + "angleDeg": [ + 189.6327156, + 41215158.1842005, + 12.711923222 + ], + "ra": 0, + "dec": 0, + "pm": -1.143 + } + ] + } }, { "id": "deimos", @@ -443,7 +790,76 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1950 Jan 1", "parentBodyId": "mars", - "rotationPeriodHours": 30.29857998656265 + "rotationPeriodHours": 30.29857998656265, + "rotationalElements": { + "poleRaDeg": [ + 316.65705808, + -0.10518014, + 0 + ], + "poleDecDeg": [ + 53.50992033, + -0.05979094, + 0 + ], + "primeMeridianDeg": [ + 79.39932954, + 285.16188899, + 0 + ], + "terms": [ + { + "angleDeg": [ + 121.46893664, + 660.22803474, + 0 + ], + "ra": 3.09217726, + "dec": 1.83936004, + "pm": -2.73954829 + }, + { + "angleDeg": [ + 231.05028581, + 660.9912354, + 0 + ], + "ra": 0.22980637, + "dec": 0.1432532, + "pm": -0.39968606 + }, + { + "angleDeg": [ + 251.37314025, + 1320.50145245, + 0 + ], + "ra": 0.06418655, + "dec": 0.01911409, + "pm": -0.06563259 + }, + { + "angleDeg": [ + 217.98635955, + 38279.9612555, + 0 + ], + "ra": 0.02533537, + "dec": -0.0148259, + "pm": -0.0291294 + }, + { + "angleDeg": [ + 196.19729402, + 19139.83628608, + 0 + ], + "ra": 0.00778695, + "dec": 0.0019243, + "pm": 0.0169916 + } + ] + } }, { "id": "io", @@ -471,7 +887,44 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1997 Jan 16", "parentBodyId": "jupiter", - "rotationPeriodHours": 42.45930625514436 + "rotationPeriodHours": 42.45930625514436, + "rotationalElements": { + "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 + } + ] + } }, { "id": "europa", @@ -499,7 +952,53 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1997 Jan 16", "parentBodyId": "jupiter", - "rotationPeriodHours": 85.22834531173994 + "rotationPeriodHours": 85.22834531173994, + "rotationalElements": { + "poleRaDeg": [ + 268.08, + -0.009, + 0 + ], + "poleDecDeg": [ + 64.51, + 0.003, + 0 + ], + "primeMeridianDeg": [ + 36.022, + 101.3747235, + 0 + ], + "terms": [ + { + "angleDeg": [ + 355.8, + 1191.3 + ], + "ra": 1.086, + "dec": 0.468, + "pm": -0.98 + }, + { + "angleDeg": [ + 119.9, + 262.1 + ], + "ra": 0.06, + "dec": 0.026, + "pm": -0.054 + }, + { + "angleDeg": [ + 229.8, + 64.3 + ], + "ra": 0.015, + "dec": 0.007, + "pm": -0.014 + } + ] + } }, { "id": "ganymede", @@ -527,7 +1026,53 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1997 Jan 16", "parentBodyId": "jupiter", - "rotationPeriodHours": 171.70927794038664 + "rotationPeriodHours": 171.70927794038664, + "rotationalElements": { + "poleRaDeg": [ + 268.2, + -0.009, + 0 + ], + "poleDecDeg": [ + 64.57, + 0.003, + 0 + ], + "primeMeridianDeg": [ + 44.064, + 50.3176081, + 0 + ], + "terms": [ + { + "angleDeg": [ + 355.8, + 1191.3 + ], + "ra": -0.037, + "dec": -0.016, + "pm": 0.033 + }, + { + "angleDeg": [ + 119.9, + 262.1 + ], + "ra": 0.431, + "dec": 0.186, + "pm": -0.389 + }, + { + "angleDeg": [ + 229.8, + 64.3 + ], + "ra": 0.091, + "dec": 0.039, + "pm": -0.082 + } + ] + } }, { "id": "callisto", @@ -555,7 +1100,53 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1997 Jan 16", "parentBodyId": "jupiter", - "rotationPeriodHours": 400.5364072574082 + "rotationPeriodHours": 400.5364072574082, + "rotationalElements": { + "poleRaDeg": [ + 268.72, + -0.009, + 0 + ], + "poleDecDeg": [ + 64.83, + 0.003, + 0 + ], + "primeMeridianDeg": [ + 259.51, + 21.5710715, + 0 + ], + "terms": [ + { + "angleDeg": [ + 119.9, + 262.1 + ], + "ra": -0.068, + "dec": -0.029, + "pm": 0.061 + }, + { + "angleDeg": [ + 229.8, + 64.3 + ], + "ra": 0.59, + "dec": 0.254, + "pm": -0.533 + }, + { + "angleDeg": [ + 113.35, + 6070 + ], + "ra": 0.01, + "dec": -0.004, + "pm": -0.009 + } + ] + } }, { "id": "mimas", @@ -583,7 +1174,44 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "saturn", - "rotationPeriodHours": 22.618127008672275 + "rotationPeriodHours": 22.618127008672275, + "rotationalElements": { + "poleRaDeg": [ + 40.66, + -0.036, + 0 + ], + "poleDecDeg": [ + 83.52, + -0.004, + 0 + ], + "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 + } + ] + } }, { "id": "enceladus", @@ -611,7 +1239,24 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "saturn", - "rotationPeriodHours": 32.88523423439451 + "rotationPeriodHours": 32.88523423439451, + "rotationalElements": { + "poleRaDeg": [ + 40.66, + -0.036, + 0 + ], + "poleDecDeg": [ + 83.52, + -0.004, + 0 + ], + "primeMeridianDeg": [ + 6.32, + 262.7318996, + 0 + ] + } }, { "id": "tethys", @@ -639,7 +1284,44 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "saturn", - "rotationPeriodHours": 45.30726088829954 + "rotationPeriodHours": 45.30726088829954, + "rotationalElements": { + "poleRaDeg": [ + 40.66, + -0.036, + 0 + ], + "poleDecDeg": [ + 83.52, + -0.004, + 0 + ], + "primeMeridianDeg": [ + 8.95, + 190.6979085, + 0 + ], + "terms": [ + { + "angleDeg": [ + 300, + -7225.9 + ], + "ra": 9.66, + "dec": -1.09, + "pm": -9.6 + }, + { + "angleDeg": [ + 316.45, + 506.2 + ], + "ra": 0, + "dec": 0, + "pm": 2.23 + } + ] + } }, { "id": "dione", @@ -667,7 +1349,24 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "saturn", - "rotationPeriodHours": 65.68597371070803 + "rotationPeriodHours": 65.68597371070803, + "rotationalElements": { + "poleRaDeg": [ + 40.66, + -0.036, + 0 + ], + "poleDecDeg": [ + 83.52, + -0.004, + 0 + ], + "primeMeridianDeg": [ + 357.6, + 131.5349316, + 0 + ] + } }, { "id": "rhea", @@ -695,7 +1394,35 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "saturn", - "rotationPeriodHours": 108.42006554798584 + "rotationPeriodHours": 108.42006554798584, + "rotationalElements": { + "poleRaDeg": [ + 40.38, + -0.036, + 0 + ], + "poleDecDeg": [ + 83.55, + -0.004, + 0 + ], + "primeMeridianDeg": [ + 235.16, + 79.6900478, + 0 + ], + "terms": [ + { + "angleDeg": [ + 345.2, + -1016.3 + ], + "ra": 3.1, + "dec": -0.35, + "pm": -3.08 + } + ] + } }, { "id": "titan", @@ -723,7 +1450,24 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "saturn", - "rotationPeriodHours": 382.69076217631203 + "rotationPeriodHours": 382.69076217631203, + "rotationalElements": { + "poleRaDeg": [ + 39.4827, + 0, + 0 + ], + "poleDecDeg": [ + 83.4279, + 0, + 0 + ], + "primeMeridianDeg": [ + 186.5855, + 22.5769768, + 0 + ] + } }, { "id": "hyperion", @@ -778,7 +1522,24 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "saturn", - "rotationPeriodHours": 1903.9469348834284 + "rotationPeriodHours": 1903.9469348834284, + "rotationalElements": { + "poleRaDeg": [ + 318.16, + -3.949, + 0 + ], + "poleDecDeg": [ + 75.03, + -1.143, + 0 + ], + "primeMeridianDeg": [ + 355.2, + 4.5379572, + 0 + ] + } }, { "id": "phoebe", @@ -806,7 +1567,24 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "saturn", - "rotationPeriodHours": 9.273966666666666 + "rotationPeriodHours": 9.273966666666666, + "rotationalElements": { + "poleRaDeg": [ + 356.9, + 0, + 0 + ], + "poleDecDeg": [ + 77.8, + 0, + 0 + ], + "primeMeridianDeg": [ + 178.58, + 931.639, + 0 + ] + } }, { "id": "miranda", @@ -834,7 +1612,62 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1980 Jan 1", "parentBodyId": "uranus", - "rotationPeriodHours": 33.9235057988244 + "rotationPeriodHours": 33.9235057988244, + "rotationalElements": { + "poleRaDeg": [ + 257.43, + 0, + 0 + ], + "poleDecDeg": [ + -15.08, + 0, + 0 + ], + "primeMeridianDeg": [ + 30.7, + -254.6906892, + 0 + ], + "terms": [ + { + "angleDeg": [ + 102.23, + -2024.22 + ], + "ra": 4.41, + "dec": 4.25, + "pm": 1.15 + }, + { + "angleDeg": [ + 316.41, + 2863.96 + ], + "ra": 0, + "dec": 0, + "pm": -1.27 + }, + { + "angleDeg": [ + 204.46, + -4048.44 + ], + "ra": -0.04, + "dec": -0.02, + "pm": -0.09 + }, + { + "angleDeg": [ + 632.82, + 5727.92 + ], + "ra": 0, + "dec": 0, + "pm": 0.15 + } + ] + } }, { "id": "ariel", @@ -862,7 +1695,44 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1980 Jan 1", "parentBodyId": "uranus", - "rotationPeriodHours": 60.48909723963263 + "rotationPeriodHours": 60.48909723963263, + "rotationalElements": { + "poleRaDeg": [ + 257.43, + 0, + 0 + ], + "poleDecDeg": [ + -15.1, + 0, + 0 + ], + "primeMeridianDeg": [ + 156.22, + -142.8356681, + 0 + ], + "terms": [ + { + "angleDeg": [ + 316.41, + 2863.96 + ], + "ra": 0, + "dec": 0, + "pm": 0.05 + }, + { + "angleDeg": [ + 304.01, + -51.94 + ], + "ra": 0.29, + "dec": 0.28, + "pm": 0.08 + } + ] + } }, { "id": "umbriel", @@ -890,7 +1760,44 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1980 Jan 1", "parentBodyId": "uranus", - "rotationPeriodHours": 99.46023494563465 + "rotationPeriodHours": 99.46023494563465, + "rotationalElements": { + "poleRaDeg": [ + 257.43, + 0, + 0 + ], + "poleDecDeg": [ + -15.1, + 0, + 0 + ], + "primeMeridianDeg": [ + 108.05, + -86.8688923, + 0 + ], + "terms": [ + { + "angleDeg": [ + 316.41, + 2863.96 + ], + "ra": 0, + "dec": 0, + "pm": -0.09 + }, + { + "angleDeg": [ + 308.71, + -93.17 + ], + "ra": 0.21, + "dec": 0.2, + "pm": 0.06 + } + ] + } }, { "id": "titania", @@ -918,7 +1825,35 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1980 Jan 1", "parentBodyId": "uranus", - "rotationPeriodHours": 208.94080635857947 + "rotationPeriodHours": 208.94080635857947, + "rotationalElements": { + "poleRaDeg": [ + 257.43, + 0, + 0 + ], + "poleDecDeg": [ + -15.1, + 0, + 0 + ], + "primeMeridianDeg": [ + 77.74, + -41.3514316, + 0 + ], + "terms": [ + { + "angleDeg": [ + 340.82, + -75.32 + ], + "ra": 0.29, + "dec": 0.28, + "pm": 0.08 + } + ] + } }, { "id": "oberon", @@ -946,7 +1881,35 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 1980 Jan 1", "parentBodyId": "uranus", - "rotationPeriodHours": 323.1176207078424 + "rotationPeriodHours": 323.1176207078424, + "rotationalElements": { + "poleRaDeg": [ + 257.43, + 0, + 0 + ], + "poleDecDeg": [ + -15.1, + 0, + 0 + ], + "primeMeridianDeg": [ + 6.77, + -26.7394932, + 0 + ], + "terms": [ + { + "angleDeg": [ + 259.14, + -504.81 + ], + "ra": 0.16, + "dec": 0.16, + "pm": 0.04 + } + ] + } }, { "id": "triton", @@ -974,7 +1937,107 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "neptune", - "rotationPeriodHours": 141.0444976486201 + "rotationPeriodHours": 141.0444976486201, + "rotationalElements": { + "poleRaDeg": [ + 299.36, + 0, + 0 + ], + "poleDecDeg": [ + 41.17, + 0, + 0 + ], + "primeMeridianDeg": [ + 296.53, + -61.2572637, + 0 + ], + "terms": [ + { + "angleDeg": [ + 177.85, + 52.316 + ], + "ra": -32.35, + "dec": 22.55, + "pm": 22.25 + }, + { + "angleDeg": [ + 355.7, + 104.632 + ], + "ra": -6.28, + "dec": 2.1, + "pm": 6.73 + }, + { + "angleDeg": [ + 533.55, + 156.948 + ], + "ra": -2.08, + "dec": 0.55, + "pm": 2.05 + }, + { + "angleDeg": [ + 711.4, + 209.264 + ], + "ra": -0.74, + "dec": 0.16, + "pm": 0.74 + }, + { + "angleDeg": [ + 889.25, + 261.58 + ], + "ra": -0.28, + "dec": 0.05, + "pm": 0.28 + }, + { + "angleDeg": [ + 1067.1, + 313.896 + ], + "ra": -0.11, + "dec": 0.02, + "pm": 0.11 + }, + { + "angleDeg": [ + 1244.95, + 366.212 + ], + "ra": -0.07, + "dec": 0.01, + "pm": 0.05 + }, + { + "angleDeg": [ + 1422.8, + 418.528 + ], + "ra": -0.02, + "dec": 0, + "pm": 0.02 + }, + { + "angleDeg": [ + 1600.65, + 470.844 + ], + "ra": -0.01, + "dec": 0, + "pm": 0.01 + } + ] + } }, { "id": "nereid", @@ -1029,7 +2092,44 @@ }, "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "neptune", - "rotationPeriodHours": 26.93555462331033 + "rotationPeriodHours": 26.93555462331033, + "rotationalElements": { + "poleRaDeg": [ + 299.27, + 0, + 0 + ], + "poleDecDeg": [ + 42.91, + 0, + 0 + ], + "primeMeridianDeg": [ + 93.38, + 320.7654228, + 0 + ], + "terms": [ + { + "angleDeg": [ + 357.85, + 52.316 + ], + "ra": 0.7, + "dec": -0.51, + "pm": -0.48 + }, + { + "angleDeg": [ + 142.61, + 2824.6 + ], + "ra": -0.05, + "dec": -0.04, + "pm": 0.04 + } + ] + } }, { "id": "charon", @@ -1058,6 +2158,23 @@ "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "pluto", "massRatio": 0.1220485755631374, - "rotationPeriodHours": 153.29335605836368 + "rotationPeriodHours": 153.29335605836368, + "rotationalElements": { + "poleRaDeg": [ + 132.993, + 0, + 0 + ], + "poleDecDeg": [ + -6.163, + 0, + 0 + ], + "primeMeridianDeg": [ + 122.695, + 56.3625225, + 0 + ] + } } ] \ No newline at end of file diff --git a/tools/etl/build.ts b/tools/etl/build.ts index f95df3a..9eaa97d 100644 --- a/tools/etl/build.ts +++ b/tools/etl/build.ts @@ -1,8 +1,9 @@ import { statSync } from 'node:fs'; import { BodyRecord, OrbitalElements } from '../../src/app/shared/models/body.model'; -import { eclipticToEquatorial, laplacePlaneToEquatorial } from '../../src/app/shared/astro/coordinates'; +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'; import { DeepSkyRecord } from '../../src/app/shared/models/deepsky.model'; import { ExoplanetRecord } from '../../src/app/shared/models/exoplanet.model'; import { StarRecord, SUN_STAR_ID } from '../../src/app/shared/models/star.model'; @@ -146,6 +147,7 @@ function validateMerge(stars: StarRecord[]): void { const MAX_PLANET_OFFSET_DEG = 0.25; const MAX_MOON_OFFSET_DEG = 2.5; 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, @@ -163,6 +165,33 @@ const KM_PER_AU = 149597870.7; */ const MOON_OFFSET_CEILINGS_DEG: Record = { mimas: 46, hyperion: 21, iapetus: 11, nereid: 3 }; +/** + * The bodies the IAU WGCCRE 2015 report gives no rotational elements for: Hyperion tumbles, and + * Nereid, Eris, Haumea and Makemake have no model. Every other body must carry them, or the + * kernel was read wrongly and the body would be drawn on an invented pole. + */ +const WITHOUT_ROTATIONAL_ELEMENTS = new Set(['hyperion', 'nereid', 'eris', 'haumea', 'makemake']); + +/** + * How far the IAU's day, 360 degrees over W's rate, may be from the one Horizons states, as a + * fraction of it. Measured on this catalogue: at most 1.8e-5 (Jupiter's System III, 9.92492 hours + * against 9.92510). Neptune is 0.89 per cent out, because the report takes 15.9663 hours from the + * cloud features Karkoschka (2011) tracked, where Horizons keeps Voyager's radio period, 16.11. What + * this catches is a rate read in the wrong unit or for the wrong body: Oberon's day for Titania's is + * 55 per cent out. + */ +const MAX_DAY_OFFSET = 1e-4; +const DAY_OFFSET_CEILINGS: Record = { neptune: 0.01 }; + +/** + * How far the tilt of the IAU's spin axis from the orbit may be from the obliquity Horizons + * states. The axis is the IAU's pole, turned end for end where W runs backwards: the report names + * a planet's north pole by the side of the solar system it lies on, whichever way the planet turns. + * Measured on this catalogue: at most 0.058 degrees (Venus, 177.358 against 177.3). Taken as the + * pole alone, Venus comes out at 2.6 degrees and Uranus at 82.2, which is what this catches. + */ +const MAX_OBLIQUITY_OFFSET_DEG = 0.1; + 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; @@ -175,6 +204,7 @@ function validateBodies(bodies: BodyRecord[], horizonsOrbits: Map 0, `Body ${body.id} has no radius; its page states it in a form the ETL does not read.`); + const rotation = body.rotationalElements; + assertCondition( + (rotation === undefined) === WITHOUT_ROTATIONAL_ELEMENTS.has(body.id), + `Body ${body.id} ${rotation ? 'has' : 'has no'} IAU rotational elements, which the report ${rotation ? 'does not give' : 'gives'} for it.` + ); + if (rotation) { + const rate = rotation.primeMeridianDeg[1]; + if (body.rotationPeriodHours !== undefined) { + const dayOffset = Math.abs(((360 / Math.abs(rate)) * 24) / Math.abs(body.rotationPeriodHours) - 1); + const dayCeiling = DAY_OFFSET_CEILINGS[body.id] ?? MAX_DAY_OFFSET; + assertCondition( + dayOffset <= dayCeiling, + `${body.name}'s IAU day, ${((360 / Math.abs(rate)) * 24).toFixed(5)} hours, is ${dayOffset.toExponential(2)} of its length from Horizons' ${Math.abs(body.rotationPeriodHours).toFixed(5)} (at most ${dayCeiling} expected).` + ); + spins.push(`${body.id} day ${dayOffset.toExponential(1)}`); + } + if (body.obliquityDeg !== undefined) { + const pole = orientationAt(rotation, horizons!.epochJd); + const pointing = raDecToUnitVector(pole.poleRaDeg / 15, pole.poleDecDeg); + const axis = { x: Math.sign(rate) * pointing.x, y: Math.sign(rate) * pointing.y, z: Math.sign(rate) * pointing.z }; + const { inclinationDeg, longitudeOfAscendingNodeDeg } = meanElementsAt(body.orbit, body.rates, horizons!.epochJd); + 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) }; + const obliquity = angleBetweenDeg(axis, body.laplacePole ? laplacePlaneToEquatorial(normal, body.laplacePole) : eclipticToEquatorial(normal)); + assertCondition( + Math.abs(obliquity - body.obliquityDeg) <= MAX_OBLIQUITY_OFFSET_DEG, + `${body.name}'s IAU spin axis is ${obliquity.toFixed(3)} degrees from its orbit's pole, where Horizons gives an obliquity of ${body.obliquityDeg} (at most ${MAX_OBLIQUITY_OFFSET_DEG} apart expected) — the pole or the sense of W was read wrongly.` + ); + spins.push(`${body.id} tilt ${obliquity.toFixed(3)}`); + } + } + if (body.kind === 'moon') { const parent = bodies.find((candidate) => candidate.id === body.parentBodyId); assertCondition(parent !== undefined, `Moon ${body.id} has no valid parentBodyId.`); @@ -233,6 +296,7 @@ function validateBodies(bodies: BodyRecord[], horizonsOrbits: Map body.kind === 'dwarf').length; assertCondition(dwarfCount === 5, `Expected the IAU's 5 dwarf planets, found ${dwarfCount}.`); console.log(` mean elements against Horizons, degrees: ${offsets.join(', ')}.`); + console.log(` IAU rotation against Horizons (day as a fraction of it, tilt in degrees): ${spins.join(', ')}.`); } function validateExoplanets(exoplanets: ExoplanetRecord[], starIds: Set): void { diff --git a/tools/etl/fetchSolarSystem.ts b/tools/etl/fetchSolarSystem.ts index bbfb8b1..aa11574 100644 --- a/tools/etl/fetchSolarSystem.ts +++ b/tools/etl/fetchSolarSystem.ts @@ -5,6 +5,8 @@ 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 { fetchPckText } from './lib/pck'; import { dataPath, ensureDataDir } from './lib/paths'; const HOURS_PER_DAY = 24; @@ -123,15 +125,17 @@ export const FREELY_SPINNING_MOONS = new Set(BODY_SPECS.filter((spec) => spec.sp * Writes `bodies.json` for the major planets, the five dwarf planets, and every moon in JPL's * mean-element table more than 100 km in mean radius — Phoebe, at 106.6, the smallest: JPL's * mean orbital elements for where they go, or the SBDB's osculating ones where there are none, - * and JPL Horizons for their size and spin. Horizons' osculating elements for the same date come - * back alongside, for `build.ts` to check the mean ones against. + * JPL Horizons for their size and spin, and the IAU's rotational elements for where their poles + * point and which face is where. Horizons' osculating elements for the same date come back + * alongside, for `build.ts` to check the mean ones against. */ export async function fetchSolarSystem(): Promise<{ bodies: BodyRecord[]; horizonsOrbits: Map }> { - console.log(`Fetching ${BODY_SPECS.length} solar-system bodies from JPL (mean elements, Horizons)...`); + console.log(`Fetching ${BODY_SPECS.length} solar-system bodies from JPL (mean elements, Horizons, NAIF's PCK)...`); const bodies: BodyRecord[] = []; const horizonsOrbits = new Map(); const planetElements = await fetchPlanetMeanElementsText(); const satelliteElements = await fetchSatelliteMeanElementsHtml(); + const pck = await fetchPckText(); const gmById = new Map(); for (const spec of BODY_SPECS) { @@ -183,6 +187,14 @@ 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, @@ -197,7 +209,8 @@ export async function fetchSolarSystem(): Promise<{ bodies: BodyRecord[]; horizo ...(spec.parentBodyId ? { parentBodyId: spec.parentBodyId } : {}), ...(parentGm !== undefined ? { massRatio: result.gmKm3PerS2! / parentGm } : {}), ...(rotationPeriodHours !== undefined ? { rotationPeriodHours } : {}), - ...((result.obliquityDeg ?? spec.obliquityDeg) !== undefined ? { obliquityDeg: result.obliquityDeg ?? spec.obliquityDeg } : {}) + ...((result.obliquityDeg ?? spec.obliquityDeg) !== undefined ? { obliquityDeg: result.obliquityDeg ?? spec.obliquityDeg } : {}), + ...(rotation ? { rotationalElements: rotation.elements } : {}) }); } diff --git a/tools/etl/lib/pck.ts b/tools/etl/lib/pck.ts new file mode 100644 index 0000000..780d1c9 --- /dev/null +++ b/tools/etl/lib/pck.ts @@ -0,0 +1,13 @@ +import { fetchTextCached } from './http'; + +/** + * NAIF's generic text PCK, which carries the IAU WGCCRE 2015 report's rotational elements + * (Archinal et al. 2018, Celest Mech Dyn Astr 130:22) for every body here that has them, periodic + * terms included, in a form a program can read rather than a table typeset in a paper. A released + * kernel is never edited, only superseded under a new name, so the URL pins the numbers. + */ +const PCK_URL = 'https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00011.tpc'; + +export async function fetchPckText(): Promise { + return fetchTextCached(PCK_URL, 'naif-pck00011.tpc'); +}