Carry every body's IAU rotational elements, read from NAIF's kernel of the 2015 report
bodies.json now holds, for 33 of the 38 bodies, the pole right ascension and declination and the prime meridian W of the IAU WGCCRE 2015 report (Archinal et al. 2018), with their rates and the periodic terms. They are read from NAIF's pck00011.tpc, which carries the report in a form a program can read, periodic terms and their angles included. Hyperion (chaotic), Nereid, Eris, Haumea and Makemake have no model in the report. The parser, src/app/shared/astro/rotational-elements.ts, sits beside the other source readers so the unit suite covers it. It reads data blocks only where \begindata stands alone on a line, as the kernel's own prose mentions the token mid-sentence. It reads the Fortran exponent (the Moon's -1.4D-12 d² term) and the degree-2 angles of the Mars system, where Phobos's tidal acceleration lives. NAIF numbers a small body 2 000 000 past its catalogue number, so Ceres is 2000001. Periodic terms under 0.01 degrees are left out. 0.01 degrees moves a point by 0.11 px on the largest body ever drawn (Jupiter at 641 px of radius). That drops 32 terms: - Mercury: 4 (0.0011 degrees and less) - the Moon: 8 of 13 (0.0072 and less) - Mars: 13 (0.00024 and less); its three 0.42-1.59 degree long-period terms stay - Phobos: 1 (0.0063) - Jupiter: 5 (0.0022 and less) - Europa: 1 (0.009) Kept, among others: Mimas's 44.85-degree libration, Triton's 32-degree precession, Miranda's 4.4 and Phobos's 1.14-degree libration. build.ts now checks the elements against Horizons on the real catalogue: - Every body but those five carries elements, and they do not. - The IAU day, 360 over W's rate, is within 1e-4 of Horizons' period. Measured: at most 1.8e-5 (Jupiter). Neptune gets a 0.01 ceiling: 0.89 per cent, because the report takes Karkoschka's 15.9663 h where Horizons keeps Voyager's 16.11. - The spin axis, the pole turned end for end where W runs backwards, is within 0.1 degrees of Horizons' obliquity. Measured: at most 0.058 (Venus, 177.358 against 177.3); Uranus 97.771, Pluto 119.610, Earth 23.435. Full npm run etl passes. Three mutants each fail it on the named check: - W's sign dropped: "Venus's IAU spin axis is 2.642 degrees". - Ceres looked up by catalogue number: "Body ceres has no IAU rotational elements". - W's rate read per century: "Mercury's IAU day ... 3.65e+4". Nothing is drawn from these yet. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5.5
@@ -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);
|
||||
});
|
||||
});
|
||||
@@ -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<string, number[]> {
|
||||
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<string, number[]>();
|
||||
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<RotationalElements['terms']> = [];
|
||||
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 };
|
||||
}
|
||||
@@ -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 }>;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user