diff --git a/src/app/features/body-detail/body-view-model.spec.ts b/src/app/features/body-detail/body-view-model.spec.ts index a76f661..4792435 100644 --- a/src/app/features/body-detail/body-view-model.spec.ts +++ b/src/app/features/body-detail/body-view-model.spec.ts @@ -34,6 +34,9 @@ const earth: BodyRecord = { kind: 'planet', radiusKm: 6371, orbit: orbit(), + // Standish's mean longitude rate, 35 999.373 degrees a century. + rates: { meanMotionDegPerDay: 35999.37306329 / 36525, longitudeOfAscendingNodeDegPerDay: 0, argumentOfPeriapsisDegPerDay: 0 }, + orbitSource: 'JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000', }; const luna: BodyRecord = { id: 'luna', @@ -43,6 +46,9 @@ const luna: BodyRecord = { radiusKm: 1737, parentBodyId: 'earth', orbit: orbit({ semiMajorAxisAu: 0.00257 }), + // JPL SSD's sidereal mean motion for the Moon. + rates: { meanMotionDegPerDay: 13.176358, longitudeOfAscendingNodeDegPerDay: -0.05299, argumentOfPeriapsisDegPerDay: 0.16435 }, + orbitSource: 'JPL SSD satellite mean elements, epoch 2000 Jan 1', }; describe('heliocentricPeriodDays', () => { diff --git a/src/app/features/galaxy-system/galaxy-system-scene.component.spec.ts b/src/app/features/galaxy-system/galaxy-system-scene.component.spec.ts index e470776..89a5a7e 100644 --- a/src/app/features/galaxy-system/galaxy-system-scene.component.spec.ts +++ b/src/app/features/galaxy-system/galaxy-system-scene.component.spec.ts @@ -6,6 +6,8 @@ import { afterEach, beforeEach, describe, expect, it, MockInstance, vi } from 'v import { DataLoaderService, StarField } from '../../core/data/data-loader.service'; import { EngineService, EngineTickCallback } from '../../core/engine/engine.service'; import { BodyRecord } from '../../shared/models/body.model'; +import { GM_SUN_AU3_PER_DAY2 } from '../../shared/astro/constants'; +import { keplerRates } from '../../shared/astro/kepler'; import { DeepSkyRecord } from '../../shared/models/deepsky.model'; import { ExoplanetRecord } from '../../shared/models/exoplanet.model'; import { StarRecord } from '../../shared/models/star.model'; @@ -63,7 +65,8 @@ const EARTH: BodyRecord = { argumentOfPeriapsisDeg: 0, meanAnomalyAtEpochDeg: 0, epochJd: 2451545.0 - } + }, + rates: keplerRates(1, GM_SUN_AU3_PER_DAY2), orbitSource: 'test' }; /** Minimal stand-in for `EngineService` that skips real WebGPU/WebGL initialization entirely, 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 daa00af..0ee9c88 100644 --- a/src/app/features/galaxy-system/system-orbits-renderer.spec.ts +++ b/src/app/features/galaxy-system/system-orbits-renderer.spec.ts @@ -1,7 +1,8 @@ import * as THREE from 'three/webgpu'; import { describe, expect, it } from 'vitest'; -import { DEFAULT_EPOCH_JD } from '../../shared/astro/constants'; +import { DEFAULT_EPOCH_JD, GM_SUN_AU3_PER_DAY2 } from '../../shared/astro/constants'; +import { keplerRates } from '../../shared/astro/kepler'; import { eclipticToEquatorial, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; import { BodyRecord } from '../../shared/models/body.model'; import { ExoplanetRecord } from '../../shared/models/exoplanet.model'; @@ -164,7 +165,8 @@ describe('SystemOrbitsRenderer exoplanet propagation', () => { argumentOfPeriapsisDeg: 0, meanAnomalyAtEpochDeg: 0, epochJd: DEFAULT_EPOCH_JD - } + }, + rates: keplerRates(1, GM_SUN_AU3_PER_DAY2), orbitSource: 'test' }; it('places an ecliptic orbit in the ecliptic plane of the equatorial scene', () => { @@ -234,7 +236,9 @@ describe('SystemOrbitsRenderer exoplanet propagation', () => { name: 'Jupiter', kind: 'planet', radiusKm: 69911, - orbit: { semiMajorAxisAu: 5.2, eccentricity: 0.048, inclinationDeg: 1.3, longitudeOfAscendingNodeDeg: 100, argumentOfPeriapsisDeg: 275, meanAnomalyAtEpochDeg: 20, epochJd: DEFAULT_EPOCH_JD } + orbit: { semiMajorAxisAu: 5.2, eccentricity: 0.048, inclinationDeg: 1.3, longitudeOfAscendingNodeDeg: 100, argumentOfPeriapsisDeg: 275, meanAnomalyAtEpochDeg: 20, epochJd: DEFAULT_EPOCH_JD }, + rates: keplerRates(5.2, GM_SUN_AU3_PER_DAY2), + orbitSource: 'test' }; /** The grid and the tethers are the only line objects the renderer adds outside a pivot. */ @@ -383,6 +387,7 @@ describe('rotation', () => { kind: 'planet', radiusKm: 6371, orbit: { semiMajorAxisAu: 1, eccentricity: 0.0167, inclinationDeg: 0, longitudeOfAscendingNodeDeg: 0, argumentOfPeriapsisDeg: 0, meanAnomalyAtEpochDeg: 0, epochJd: DEFAULT_EPOCH_JD }, + rates: keplerRates(1, GM_SUN_AU3_PER_DAY2), orbitSource: 'test', rotationPeriodHours: 23.934, obliquityDeg: 23.4392911, ...overrides @@ -464,3 +469,63 @@ describe('exoplanet size without a measured radius', () => { expect(radiusOf({ radiusEarth: 1.88, massEarth: 2829 }) / EARTH_AU).toBeCloseTo(1.88, 2); }); }); + +describe('solar-system bodies against Horizons', () => { + // Real records from bodies.json, and Horizons' own positions for them (ICRF, AU; heliocentric + // for the planets, planet-centred for the moons) at dates across 1950-2100, so the whole path — + // mean elements, their rates, the Laplace planes and the scene's frame — is checked against + // JPL's ephemeris rather than against itself. + const RECORDS: Record> = { + earth: {kind: 'planet', orbit: {semiMajorAxisAu: 1.00000018, eccentricity: 0.01673163, inclinationDeg: -0.00054346, longitudeOfAscendingNodeDeg: -5.11260389, argumentOfPeriapsisDeg: 108.04266274, meanAnomalyAtEpochDeg: -2.4631431299999917, epochJd: 2451545}, rates: {meanMotionDegPerDay: 0.9856091187759068, longitudeOfAscendingNodeDegPerDay: -0.000006604751813826146, argumentOfPeriapsisDegPerDay: 0.000015309819575633124, semiMajorAxisAuPerDay: -8.213552361396303e-13, eccentricityPerDay: -1.002327173169062e-9, inclinationDegPerDay: -3.6609938398357287e-7}}, + jupiter: {kind: 'planet', orbit: {semiMajorAxisAu: 5.20248019, eccentricity: 0.0485359, inclinationDeg: 1.29861416, longitudeOfAscendingNodeDeg: 100.29282654, argumentOfPeriapsisDeg: -86.0178741, meanAnomalyAtEpochDeg: 20.059839080000003, epochJd: 2451545}, rates: {meanMotionDegPerDay: 0.08309113532019165, longitudeOfAscendingNodeDegPerDay: 0.0000035659463381245725, argumentOfPeriapsisDegPerDay: 0.0000014167219712525667, semiMajorAxisAuPerDay: -7.841204654346339e-10, eccentricityPerDay: 4.935249828884326e-9, inclinationDegPerDay: -8.83501711156742e-8, meanAnomalyTerms: {b: -0.00012452, c: 0.0606406, s: -0.35635438, f: 38.35125}}}, + saturn: {kind: 'planet', orbit: {semiMajorAxisAu: 9.54149883, eccentricity: 0.05550825, inclinationDeg: 2.49424102, longitudeOfAscendingNodeDeg: 113.63998702, argumentOfPeriapsisDeg: -20.778626390000014, meanAnomalyAtEpochDeg: -42.78564733999999, epochJd: 2451545}, rates: {meanMotionDegPerDay: 0.033459683702669406, longitudeOfAscendingNodeDegPerDay: -0.000006848734291581108, argumentOfPeriapsisDegPerDay: 0.000021682266940451745, semiMajorAxisAuPerDay: -8.391512662559891e-10, eccentricityPerDay: -8.773169062286106e-9, inclinationDegPerDay: 1.2374236824093085e-7, meanAnomalyTerms: {b: 0.00025899, c: -0.13434469, s: 0.87320147, f: 38.35125}}}, + neptune: {kind: 'planet', orbit: {semiMajorAxisAu: 30.06952752, eccentricity: 0.00895439, inclinationDeg: 1.7700552, longitudeOfAscendingNodeDeg: 131.78635853, argumentOfPeriapsisDeg: -85.10477129, meanAnomalyAtEpochDeg: 257.54130563, epochJd: 2451545}, rates: {meanMotionDegPerDay: 0.005981249914852841, longitudeOfAscendingNodeDegPerDay: -1.6599644079397672e-7, argumentOfPeriapsisDegPerDay: 4.4250239561943875e-7, semiMajorAxisAuPerDay: 1.7650924024640657e-9, eccentricityPerDay: 2.2395619438740589e-10, inclinationDegPerDay: 6.132785763175907e-9, meanAnomalyTerms: {b: -0.00041348, c: 0.68346318, s: -0.10162547, f: 7.67025}}}, + pluto: {kind: 'dwarf', orbit: {semiMajorAxisAu: 39.48686035, eccentricity: 0.24885238, inclinationDeg: 17.1410426, longitudeOfAscendingNodeDeg: 110.30167986, argumentOfPeriapsisDeg: 113.79534612000002, meanAnomalyAtEpochDeg: 14.86832412999999, epochJd: 2451545}, rates: {meanMotionDegPerDay: 0.003974823518959616, longitudeOfAscendingNodeDegPerDay: -2.2176071184120468e-7, argumentOfPeriapsisDegPerDay: -4.3489664613278575e-8, semiMajorAxisAuPerDay: 1.2313511293634495e-7, eccentricityPerDay: 1.6470910335386722e-9, inclinationDegPerDay: 1.3716632443531827e-10, meanAnomalyTerms: {b: -0.01262724, c: 0, s: 0, f: 0}}}, + moon: {kind: 'moon', orbit: {semiMajorAxisAu: 0.0025695552897999907, eccentricity: 0.0554, inclinationDeg: 5.16, longitudeOfAscendingNodeDeg: 125.08, argumentOfPeriapsisDeg: 318.15, meanAnomalyAtEpochDeg: 135.27, epochJd: 2451545}, rates: {meanMotionDegPerDay: 13.176358, longitudeOfAscendingNodeDegPerDay: -0.052990660396105185, argumentOfPeriapsisDegPerDay: 0.164353223839846}, parentBodyId: 'earth'}, + io: {kind: 'moon', orbit: {semiMajorAxisAu: 0.0028195588481728304, eccentricity: 0.0041, inclinationDeg: 0.036, longitudeOfAscendingNodeDeg: 43.977, argumentOfPeriapsisDeg: 84.129, meanAnomalyAtEpochDeg: 342.021, epochJd: 2450464.5}, rates: {meanMotionDegPerDay: 203.4889583, longitudeOfAscendingNodeDegPerDay: -0.1328337309120696, argumentOfPeriapsisDegPerDay: -0.6065392513031117}, laplacePole: {raDeg: 268.057, decDeg: 64.495}, parentBodyId: 'jupiter'}, + europa: {kind: 'moon', orbit: {semiMajorAxisAu: 0.004486026417754354, eccentricity: 0.0094, inclinationDeg: 0.466, longitudeOfAscendingNodeDeg: 219.106, argumentOfPeriapsisDeg: 88.97, meanAnomalyAtEpochDeg: 171.016, epochJd: 2450464.5}, rates: {meanMotionDegPerDay: 101.3747242, longitudeOfAscendingNodeDegPerDay: -0.03265393199600969, argumentOfPeriapsisDegPerDay: -0.7070489837643877}, laplacePole: {raDeg: 268.084, decDeg: 64.506}, parentBodyId: 'jupiter'}, + titan: {kind: 'moon', orbit: {semiMajorAxisAu: 0.008167663044150534, eccentricity: 0.0288, inclinationDeg: 0.306, longitudeOfAscendingNodeDeg: 28.06, argumentOfPeriapsisDeg: 180.532, meanAnomalyAtEpochDeg: 163.31, epochJd: 2451545}, rates: {meanMotionDegPerDay: 22.5769756, longitudeOfAscendingNodeDegPerDay: -0.001398845136769169, argumentOfPeriapsisDegPerDay: 0.002799120423059061}, laplacePole: {raDeg: 36.214, decDeg: 83.949}, parentBodyId: 'saturn'}, + triton: {kind: 'moon', orbit: {semiMajorAxisAu: 0.002371417442908832, eccentricity: 0, inclinationDeg: 156.865, longitudeOfAscendingNodeDeg: 177.608, argumentOfPeriapsisDeg: 66.142, meanAnomalyAtEpochDeg: 352.257, epochJd: 2451545}, rates: {meanMotionDegPerDay: 61.2572638, longitudeOfAscendingNodeDegPerDay: 0.001433750844964632, argumentOfPeriapsisDegPerDay: 0.0025509841146658433}, laplacePole: {raDeg: 299.456, decDeg: 43.414}, parentBodyId: 'neptune'}, + }; + // 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. + 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], + ['saturn', 2478938.5, -3.51309768447752, -8.723317933082274, -3.452662390556131, 0.25], + ['pluto', 2442413.5, -29.2488165026956, -7.1421817246801, 6.58403957591589, 0.1], + ['moon', 2469807.5, 0.00240364781322315, 0.0006554283236619424, 0.0004472719300783614, 2], + ['io', 2433282.5, 0.0004488349204269952, 0.002519633434577752, 0.00120678715190893, 0.05], + ['europa', 2433282.5, 0.004084372287322533, -0.001665375585011311, -0.0007673072324795899, 0.1], + ['titan', 2488069.5, 0.007800850235156121, -0.001556932380983438, -0.0006078959246502567, 0.05], + ['triton', 2488069.5, -0.001421151845853369, -0.0001894510477241482, 0.001888790702926415, 0.2], + ]; + + function record(id: string): BodyRecord { + return { id, systemStarId: 0, name: id, radiusKm: 1000, orbitSource: 'test', ...RECORDS[id] }; + } + + const renderer = new SystemOrbitsRenderer(Object.keys(RECORDS).map(record), []); + + for (const [id, jd, x, y, z, maxDeg] of HORIZONS) { + it(`puts ${id} within ${maxDeg} degrees of Horizons on JD ${jd}`, () => { + renderer.update(jd); + const drawn = renderer.members.find((member) => member.id === id)!.marker.position; + const angleDeg = (drawn.angleTo(new THREE.Vector3(x, y, z)) * 180) / Math.PI; + expect(angleDeg).toBeLessThan(maxDeg); + }); + } + + it('turns the Moon’s drawn orbit with its node, so the Moon stays on its own line', () => { + // Half the node's 18.6-year turn on, the ellipse drawn at the epoch has the Moon 10 degrees off + // its plane at the worst. + const moon = renderer.members.find((member) => member.id === 'moon')!.marker; + const line = moon.parent!.children.find((child) => child.name === 'orbit-line')!; + for (const days of [0, 1700, 3397, 3400]) { + renderer.update(DEFAULT_EPOCH_JD + days); + const normal = new THREE.Vector3(0, 0, 1).applyQuaternion(line.quaternion); + expect(Math.abs(moon.position.clone().normalize().dot(normal))).toBeLessThan(1e-9); + } + }); +}); diff --git a/src/app/features/galaxy-system/system-orbits-renderer.ts b/src/app/features/galaxy-system/system-orbits-renderer.ts index 37f5076..0922353 100644 --- a/src/app/features/galaxy-system/system-orbits-renderer.ts +++ b/src/app/features/galaxy-system/system-orbits-renderer.ts @@ -1,13 +1,12 @@ import * as THREE from 'three/webgpu'; import { appearanceForBody, appearanceForExoplanet } from '../../shared/astro/body-appearance'; -import { gmForParent } from '../../shared/astro/constants'; import { PlanetAppearance } from '../../shared/astro/planet-appearance'; import { planetTexture } from '../../shared/rendering/procedural-planet-texture'; import { bodyTexturePath, loadCachedTexture } from '../../shared/rendering/texture-catalog'; -import { isPropagatableOrbit, orbitEllipsePoints, propagateOrbit, resolveGravitationalParameter, resolveOrbitalElements } from '../../shared/astro/kepler'; -import { CartesianCoordinates, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; -import { BodyRecord, OrbitalElements } from '../../shared/models/body.model'; +import { isPropagatableOrbit, keplerRates, meanElementsAt, orbitEllipsePoints, positionAtEpoch, resolveGravitationalParameter, resolveOrbitalElements } from '../../shared/astro/kepler'; +import { CartesianCoordinates, laplacePlaneToEquatorial, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; +import { BodyRecord, MeanElementRates, OrbitalElements } from '../../shared/models/body.model'; import { bodyMarkerRadiusAu, systemGridRingsAu } from './system-framing'; import { PolarGridPlane, TetherField } from './grid-plane'; import { ExoplanetRecord } from '../../shared/models/exoplanet.model'; @@ -45,11 +44,45 @@ const SYSTEM_TETHER_OPACITY = 0.3; /** * Rotation carrying the **ecliptic** frame into the scene's equatorial one — a turn of the - * obliquity about the shared vernal-equinox axis. Solar-system elements come from Horizons - * against the ecliptic, so this is their frame. + * obliquity about the shared vernal-equinox axis. The planets' and the Moon's mean elements are + * given against the J2000 ecliptic, so this is their frame. */ const ECLIPTIC_FRAME = new THREE.Quaternion().setFromAxisAngle(new THREE.Vector3(1, 0, 0), OBLIQUITY_J2000_DEG * DEG_TO_RAD); +/** + * Rotation carrying a moon's element frame into the scene: its local Laplace plane where JPL + * gives one, the ecliptic otherwise. Built from the three axes {@link laplacePlaneToEquatorial} + * sends, so the scene and the ETL's check against Horizons share the one conversion. + */ +function moonFrame(body: BodyRecord): THREE.Quaternion { + const pole = body.laplacePole; + if (!pole) { + return ECLIPTIC_FRAME.clone(); + } + const axis = (x: number, y: number, z: number): THREE.Vector3 => { + const turned = laplacePlaneToEquatorial({ x, y, z }, pole); + return new THREE.Vector3(turned.x, turned.y, turned.z); + }; + return new THREE.Quaternion().setFromRotationMatrix(new THREE.Matrix4().makeBasis(axis(1, 0, 0), axis(0, 1, 0), axis(0, 0, 1))); +} + +const X_AXIS = new THREE.Vector3(1, 0, 0); +const Z_AXIS = new THREE.Vector3(0, 0, 1); +const scratchTurn = new THREE.Quaternion(); + +/** + * Sets `target` to the rotation carrying an orbit's own plane, periapsis along +X, into the + * scene: the argument of periapsis, then the inclination, then the node, as + * `positionAtTrueAnomaly` turns a point, and then the frame the elements are measured in. + */ +function orientOrbit(target: THREE.Quaternion, elements: OrbitalElements, frame: THREE.Quaternion): THREE.Quaternion { + return target + .copy(frame) + .multiply(scratchTurn.setFromAxisAngle(Z_AXIS, elements.longitudeOfAscendingNodeDeg * DEG_TO_RAD)) + .multiply(scratchTurn.setFromAxisAngle(X_AXIS, elements.inclinationDeg * DEG_TO_RAD)) + .multiply(scratchTurn.setFromAxisAngle(Z_AXIS, elements.argumentOfPeriapsisDeg * DEG_TO_RAD)); +} + /** * Rotation carrying the frame an **exoplanet's** elements are measured in into the scene. * @@ -94,17 +127,23 @@ function colorForKind(kind: SystemMemberKind): THREE.Color { /** Marks orbit lines so the whole layer can be toggled without touching the bodies. */ const ORBIT_LINE_NAME = 'orbit-line'; +/** + * The orbit's ellipse, drawn in its own plane and turned into place by the line's quaternion (see + * {@link orientOrbit}), which `update` sets again each tick: a node and a periapsis that turn cost + * a quaternion rather than a new geometry. The Moon's node goes right round in 18.6 years, so an + * ellipse fixed at one date has the Moon up to 2 sin 5.16° of its distance, 69 000 km, off its own + * line nine years on. + * + * The shape is the epoch's. The planets' axes and eccentricities do drift, but Pluto's axis, the + * fastest, moves 0.0045 AU a century against 39.5, which no drawn line shows. + */ function buildOrbitLine(elements: OrbitalElements, kind: SystemMemberKind, frame: THREE.Quaternion): THREE.Line { - const points = orbitEllipsePoints(elements); + const points = orbitEllipsePoints({ ...elements, inclinationDeg: 0, longitudeOfAscendingNodeDeg: 0, argumentOfPeriapsisDeg: 0 }); const positions = new Float32Array(points.length * 3); - const scratch = new THREE.Vector3(); points.forEach((point, index) => { - // Elements are measured against their source's own reference plane; `frame` rotates that - // plane into the scene's equatorial one. - const { x, y, z } = scratch.set(point.x, point.y, point.z).applyQuaternion(frame); - positions[index * 3] = x; - positions[index * 3 + 1] = y; - positions[index * 3 + 2] = z; + positions[index * 3] = point.x; + positions[index * 3 + 1] = point.y; + positions[index * 3 + 2] = point.z; }); const geometry = new THREE.BufferGeometry(); @@ -118,6 +157,7 @@ function buildOrbitLine(elements: OrbitalElements, kind: SystemMemberKind, frame const line = new THREE.Line(geometry, material); line.name = ORBIT_LINE_NAME; + orientOrbit(line.quaternion, elements, frame); return line; } @@ -201,8 +241,9 @@ const HOURS_PER_DAY = 24; * The obliquity fixes how far the pole leans from the orbit normal, and nothing more: which way * it leans needs the pole's right ascension, which the Horizons pages this reads do not carry. The * lean is taken about the orbit's ascending node because that is the one line the elements name, - * not because the data says so — so the tilt is real and its azimuth is not. Likewise the phase: - * each body starts at its elements' epoch (2025-01-01 here) in an arbitrary orientation, the + * not because the data says so — so the tilt is real and its azimuth is not. It is the node on the + * date drawn, so the lean follows the orbit as the node turns. Likewise the phase: each body + * starts at its elements' epoch (J2000 for the planets) in an arbitrary orientation, the * shortest rotation of +Y onto its axis, and turns from there. The rate and the sense are real; * the face towards the camera is not. * @@ -212,7 +253,7 @@ const HOURS_PER_DAY = 24; * obliquity is given it carries the sense, and the period is taken as a magnitude; the sign of the * period is only read for a body with no obliquity at all. */ -function spinFor(elements: OrbitalElements, frame: THREE.Quaternion, rotationPeriodHours: number, obliquityDeg: number | undefined, epochJd: number): THREE.Quaternion { +function spinFor(elements: OrbitalElements, frame: THREE.Quaternion, rotationPeriodHours: number, obliquityDeg: number | undefined, daysSinceEpoch: number): THREE.Quaternion { const node = elements.longitudeOfAscendingNodeDeg * DEG_TO_RAD; const inclination = elements.inclinationDeg * DEG_TO_RAD; const nodeDirection = new THREE.Vector3(Math.cos(node), Math.sin(node), 0); @@ -220,7 +261,7 @@ function spinFor(elements: OrbitalElements, frame: THREE.Quaternion, rotationPer .applyAxisAngle(nodeDirection, (obliquityDeg ?? 0) * DEG_TO_RAD) .applyQuaternion(frame); const period = obliquityDeg === undefined ? rotationPeriodHours : Math.abs(rotationPeriodHours); - const turns = ((epochJd - elements.epochJd) * HOURS_PER_DAY) / period; + const turns = (daysSinceEpoch * HOURS_PER_DAY) / period; return new THREE.Quaternion() .setFromUnitVectors(SPIN_AXIS, axis) .multiply(new THREE.Quaternion().setFromAxisAngle(SPIN_AXIS, turns * 2 * Math.PI)); @@ -230,8 +271,9 @@ interface TrackedTopLevelBody { id: string; kind: SystemMemberKind; elements: OrbitalElements; - gmAu3PerDay2: number; + rates: MeanElementRates; marker: THREE.Mesh; + orbitLine: THREE.Line; /** Rotation from this body's own element frame into the scene's equatorial one. */ frame: THREE.Quaternion; /** AU position last computed for this body; moons read their parent's here. */ @@ -244,8 +286,9 @@ interface TrackedTopLevelBody { interface TrackedMoon { id: string; elements: OrbitalElements; - gmAu3PerDay2: number; + rates: MeanElementRates; marker: THREE.Mesh; + orbitLine: THREE.Line; frame: THREE.Quaternion; pivot: THREE.Group; parentId: string; @@ -322,7 +365,7 @@ export class SystemOrbitsRenderer { } // A body reaches here only when it has no parentBodyId, so `kind` is 'planet' or 'dwarf'. const kind: SystemMemberKind = body.kind; - const tracked = this.addTopLevelBody(body.id, kind, body.orbit, gmForParent(undefined), body.radiusKm, ECLIPTIC_FRAME, appearanceForBody(body, bodies, hostLuminositySolar), { periodHours: body.rotationPeriodHours, obliquityDeg: body.obliquityDeg }); + const tracked = this.addTopLevelBody(body.id, kind, body.orbit, body.rates, body.radiusKm, ECLIPTIC_FRAME, appearanceForBody(body, bodies, hostLuminositySolar), { periodHours: body.rotationPeriodHours, obliquityDeg: body.obliquityDeg }); members.push({ id: body.id, kind, marker: tracked.marker }); } @@ -335,7 +378,7 @@ export class SystemOrbitsRenderer { if (!parentTracked) { continue; // orphaned moon reference; skip rather than crash. } - const moon = this.addMoon(body.id, body.orbit, gmForParent(body.parentBodyId), body.radiusKm, parentTracked, ECLIPTIC_FRAME, appearanceForBody(body, bodies, hostLuminositySolar), { periodHours: body.rotationPeriodHours, obliquityDeg: body.obliquityDeg }); + const moon = this.addMoon(body.id, body.orbit, body.rates, body.radiusKm, parentTracked, moonFrame(body), appearanceForBody(body, bodies, hostLuminositySolar), { periodHours: body.rotationPeriodHours, obliquityDeg: body.obliquityDeg }); members.push({ id: body.id, kind: 'moon', marker: moon.marker, parentId: parent.id }); } @@ -353,21 +396,21 @@ export class SystemOrbitsRenderer { const elements = resolveOrbitalElements(exoplanet.orbit); const radiusEarth = exoplanet.radiusEarth ?? radiusFromMassEarth(exoplanet.massEarth); const radiusKm = radiusEarth ? radiusEarth * EARTH_RADIUS_KM : undefined; - // Not `gmForParent(undefined)`: that assumes a solar-mass host for every system, and - // most exoplanet hosts are red dwarfs a fraction of the Sun's mass. + // Not the Sun's: that assumes a solar-mass host for every system, and most exoplanet hosts + // are red dwarfs a fraction of the Sun's mass. const gm = resolveGravitationalParameter({ semiMajorAxisAu: exoplanet.orbit.semiMajorAxisAu, periodDays: exoplanet.periodDays, hostStarMassSolar: exoplanet.hostStarMassSolar }); - const tracked = this.addTopLevelBody(exoplanet.id, 'exoplanet', elements, gm, radiusKm, exoplanetFrame, appearanceForExoplanet(exoplanet, hostLuminositySolar)); + const tracked = this.addTopLevelBody(exoplanet.id, 'exoplanet', elements, keplerRates(elements.semiMajorAxisAu, gm), radiusKm, exoplanetFrame, appearanceForExoplanet(exoplanet, hostLuminositySolar)); members.push({ id: exoplanet.id, kind: 'exoplanet', marker: tracked.marker }); } this.members = members; // Which plane the system is read against follows from where its elements came from. Only the - // Sun has Horizons bodies and no system has both, so this is a choice between the two rather + // Sun has JPL bodies and no system has both, so this is a choice between the two rather // than a compromise: the ecliptic if there are solar-system bodies, the sky plane otherwise. this.referenceFrame = bodies.some((body) => !body.parentBodyId) ? ECLIPTIC_FRAME.clone() : exoplanetFrame; @@ -403,11 +446,13 @@ export class SystemOrbitsRenderer { /** Recomputes every marker's position for the given Julian date. Call once per tick. */ update(epochJd: number): void { for (const body of this.topLevelBodies) { - const orbital = propagateOrbit(body.elements, body.gmAu3PerDay2, epochJd); + const current = meanElementsAt(body.elements, body.rates, epochJd); + const orbital = positionAtEpoch(current); body.position.set(orbital.x, orbital.y, orbital.z).applyQuaternion(body.frame); body.marker.position.copy(body.position); + orientOrbit(body.orbitLine.quaternion, current, body.frame); if (body.rotationPeriodHours) { - body.marker.quaternion.copy(spinFor(body.elements, body.frame, body.rotationPeriodHours, body.obliquityDeg, epochJd)); + body.marker.quaternion.copy(spinFor(current, body.frame, body.rotationPeriodHours, body.obliquityDeg, epochJd - body.elements.epochJd)); } } @@ -417,10 +462,12 @@ export class SystemOrbitsRenderer { continue; } moon.pivot.position.copy(parent.position); - const orbital = propagateOrbit(moon.elements, moon.gmAu3PerDay2, epochJd); + const current = meanElementsAt(moon.elements, moon.rates, epochJd); + const orbital = positionAtEpoch(current); moon.marker.position.set(orbital.x, orbital.y, orbital.z).applyQuaternion(moon.frame); + orientOrbit(moon.orbitLine.quaternion, current, moon.frame); if (moon.rotationPeriodHours) { - moon.marker.quaternion.copy(spinFor(moon.elements, moon.frame, moon.rotationPeriodHours, moon.obliquityDeg, epochJd)); + moon.marker.quaternion.copy(spinFor(current, moon.frame, moon.rotationPeriodHours, moon.obliquityDeg, epochJd - moon.elements.epochJd)); } } @@ -473,7 +520,7 @@ export class SystemOrbitsRenderer { id: string, kind: SystemMemberKind, elements: OrbitalElements, - gmAu3PerDay2: number, + rates: MeanElementRates, radiusKm: number | undefined, frame: THREE.Quaternion, appearance?: PlanetAppearance, @@ -485,7 +532,7 @@ export class SystemOrbitsRenderer { this.trackDisposable(orbitLine.geometry, orbitLine.material as THREE.Material); this.trackDisposable(marker.geometry, marker.material as THREE.Material); - const tracked: TrackedTopLevelBody = { id, kind, elements, gmAu3PerDay2, marker, frame, position: new THREE.Vector3(), rotationPeriodHours: rotation?.periodHours, obliquityDeg: rotation?.obliquityDeg }; + const tracked: TrackedTopLevelBody = { id, kind, elements, rates, marker, orbitLine, frame, position: new THREE.Vector3(), rotationPeriodHours: rotation?.periodHours, obliquityDeg: rotation?.obliquityDeg }; this.topLevelBodies.push(tracked); return tracked; } @@ -493,7 +540,7 @@ export class SystemOrbitsRenderer { private addMoon( id: string, elements: OrbitalElements, - gmAu3PerDay2: number, + rates: MeanElementRates, radiusKm: number | undefined, parent: TrackedTopLevelBody, frame: THREE.Quaternion, @@ -508,7 +555,7 @@ export class SystemOrbitsRenderer { this.trackDisposable(orbitLine.geometry, orbitLine.material as THREE.Material); this.trackDisposable(marker.geometry, marker.material as THREE.Material); - const moon: TrackedMoon = { id, elements, gmAu3PerDay2, marker, frame, pivot, parentId: parent.id, rotationPeriodHours: rotation?.periodHours, obliquityDeg: rotation?.obliquityDeg }; + const moon: TrackedMoon = { id, elements, rates, marker, orbitLine, frame, pivot, parentId: parent.id, rotationPeriodHours: rotation?.periodHours, obliquityDeg: rotation?.obliquityDeg }; this.moons.push(moon); return moon; } diff --git a/src/app/features/search/search.component.spec.ts b/src/app/features/search/search.component.spec.ts index f5eeb7f..b2fa00c 100644 --- a/src/app/features/search/search.component.spec.ts +++ b/src/app/features/search/search.component.spec.ts @@ -34,7 +34,9 @@ const IO: BodyRecord = { argumentOfPeriapsisDeg: 0, meanAnomalyAtEpochDeg: 0, epochJd: 2451545.0 - } + }, + rates: { meanMotionDegPerDay: 203.4889583, longitudeOfAscendingNodeDegPerDay: 0, argumentOfPeriapsisDegPerDay: 0 }, + orbitSource: 'test' }; const PROXIMA_B: ExoplanetRecord = { diff --git a/src/app/shared/astro/body-appearance.spec.ts b/src/app/shared/astro/body-appearance.spec.ts index 4d49c4f..80ec84f 100644 --- a/src/app/shared/astro/body-appearance.spec.ts +++ b/src/app/shared/astro/body-appearance.spec.ts @@ -5,12 +5,13 @@ import { ExoplanetRecord } from '../models/exoplanet.model'; import { appearanceForBody, appearanceForExoplanet, heliocentricDistanceAu } from './body-appearance'; import { DEFAULT_EPOCH_JD } from './constants'; +const RATES = { meanMotionDegPerDay: 1, longitudeOfAscendingNodeDegPerDay: 0, argumentOfPeriapsisDegPerDay: 0 }; const ORBIT = { eccentricity: 0, inclinationDeg: 0, longitudeOfAscendingNodeDeg: 0, argumentOfPeriapsisDeg: 0, meanAnomalyAtEpochDeg: 0, epochJd: DEFAULT_EPOCH_JD }; -const JUPITER: BodyRecord = { id: 'jupiter', systemStarId: 0, name: 'Jupiter', kind: 'planet', radiusKm: 69911, orbit: { ...ORBIT, semiMajorAxisAu: 5.204 } }; +const JUPITER: BodyRecord = { id: 'jupiter', systemStarId: 0, name: 'Jupiter', kind: 'planet', radiusKm: 69911, orbit: { ...ORBIT, semiMajorAxisAu: 5.204 }, rates: RATES, orbitSource: 'test' }; /** Europa's own orbit is around Jupiter: 671,000 km, which is 0.00449 AU. */ -const EUROPA: BodyRecord = { id: 'europa', systemStarId: 0, name: 'Europa', kind: 'moon', radiusKm: 1560, parentBodyId: 'jupiter', orbit: { ...ORBIT, semiMajorAxisAu: 0.00449 } }; -const EARTH: BodyRecord = { id: 'earth', systemStarId: 0, name: 'Earth', kind: 'planet', radiusKm: 6371, orbit: { ...ORBIT, semiMajorAxisAu: 1 } }; +const EUROPA: BodyRecord = { id: 'europa', systemStarId: 0, name: 'Europa', kind: 'moon', radiusKm: 1560, parentBodyId: 'jupiter', orbit: { ...ORBIT, semiMajorAxisAu: 0.00449 }, rates: RATES, orbitSource: 'test' }; +const EARTH: BodyRecord = { id: 'earth', systemStarId: 0, name: 'Earth', kind: 'planet', radiusKm: 6371, orbit: { ...ORBIT, semiMajorAxisAu: 1 }, rates: RATES, orbitSource: 'test' }; const ORPHAN: BodyRecord = { ...EUROPA, id: 'orphan', parentBodyId: 'nowhere' }; const BODIES = [JUPITER, EUROPA, EARTH, ORPHAN]; diff --git a/src/app/shared/astro/constants.ts b/src/app/shared/astro/constants.ts index 4fb1b8f..9ccb835 100644 --- a/src/app/shared/astro/constants.ts +++ b/src/app/shared/astro/constants.ts @@ -19,33 +19,6 @@ export const DEFAULT_EPOCH_JD = 2451545.0; */ export const GM_SUN_AU3_PER_DAY2 = 0.01720209895 * 0.01720209895; -/** - * Approximate planet/Sun mass ratios for the major planets that host moons in `bodies.json`. - * Used to derive each planet's gravitational parameter (for propagating its moons) as - * `GM_SUN_AU3_PER_DAY2 * massRatio`. Precise enough for visualization; not JPL-grade. - */ -const PLANET_TO_SUN_MASS_RATIO: Record = { - earth: 3.003e-6, - mars: 3.227e-7, - jupiter: 9.545e-4, - saturn: 2.857e-4, - uranus: 4.365e-5, - neptune: 5.151e-5 -}; - -/** - * Gravitational parameter (AU^3/day^2) to use when propagating a body's orbit: the Sun's - * for planets/dwarfs/exoplanets, or the host planet's (derived from its Sun mass ratio) for - * moons. Falls back to the Sun's GM if `parentBodyId` isn't a known planet. - */ -export function gmForParent(parentBodyId: string | undefined): number { - if (!parentBodyId) { - return GM_SUN_AU3_PER_DAY2; - } - const massRatio = PLANET_TO_SUN_MASS_RATIO[parentBodyId]; - return massRatio ? GM_SUN_AU3_PER_DAY2 * massRatio : GM_SUN_AU3_PER_DAY2; -} - /** Converts a JS `Date` into a Julian date (days), for driving the Kepler propagator "now". */ export function dateToJulianDate(date: Date = new Date()): number { return date.getTime() / 86400000 + 2440587.5; diff --git a/src/app/shared/astro/coordinates.spec.ts b/src/app/shared/astro/coordinates.spec.ts index 6812d4c..4c4d6f6 100644 --- a/src/app/shared/astro/coordinates.spec.ts +++ b/src/app/shared/astro/coordinates.spec.ts @@ -4,6 +4,7 @@ import { distanceBetween, eclipticToEquatorial, equatorialToEcliptic, + laplacePlaneToEquatorial, OBLIQUITY_J2000_DEG, parallaxMasToParsecs, parseSexagesimal, @@ -206,6 +207,26 @@ describe('eclipticToEquatorial', () => { }); }); +describe('laplacePlaneToEquatorial', () => { + const RAD = Math.PI / 180; + /** Jupiter's moons' Laplace pole, as JPL gives it for Io. */ + const POLE = { raDeg: 268.057, decDeg: 64.495 }; + + it('sends the plane’s own pole to the right ascension and declination it is named by', () => { + const pole = laplacePlaneToEquatorial({ x: 0, y: 0, z: 1 }, POLE); + expect(Math.asin(pole.z) / RAD).toBeCloseTo(POLE.decDeg, 9); + expect(((Math.atan2(pole.y, pole.x) / RAD) + 360) % 360).toBeCloseTo(POLE.raDeg, 9); + }); + + it('counts the node from where the plane rises through the equator, 90 degrees past the pole', () => { + const node = laplacePlaneToEquatorial({ x: 1, y: 0, z: 0 }, POLE); + expect(node.z).toBeCloseTo(0, 12); + expect(((Math.atan2(node.y, node.x) / RAD) + 360) % 360).toBeCloseTo((POLE.raDeg + 90) % 360, 9); + // Rising: a quarter-turn on along the plane is north of the equator. + expect(laplacePlaneToEquatorial({ x: 0, y: 1, z: 0 }, POLE).z).toBeGreaterThan(0); + }); +}); + describe('equatorialToEcliptic', () => { it('is the exact inverse of eclipticToEquatorial', () => { for (const point of [ diff --git a/src/app/shared/astro/coordinates.ts b/src/app/shared/astro/coordinates.ts index be08da2..532447c 100644 --- a/src/app/shared/astro/coordinates.ts +++ b/src/app/shared/astro/coordinates.ts @@ -59,8 +59,8 @@ export const OBLIQUITY_J2000_DEG = 23.4392911; * * The app has to span both because its two sources disagree. Star positions come from HYG as * equatorial coordinates, which `raDecDistanceToXyz` produces and which the galaxy view renders - * directly. Orbital elements come from JPL Horizons, whose default reference plane for element - * output is the ecliptic — the ETL never overrides it. The two are tilted + * directly. The planets' and the Moon's orbital elements are JPL mean elements against the J2000 + * ecliptic. The two are tilted * {@link OBLIQUITY_J2000_DEG} apart about the shared vernal-equinox axis, so orbits have to be * rotated before they can share a scene with the stars. */ @@ -78,6 +78,29 @@ export function eclipticToEquatorial(position: CartesianCoordinates): CartesianC }; } +/** + * Rotates a vector from a moon's local **Laplace plane** frame into the equatorial one. + * + * JPL gives the giant planets' moons against the plane their orbits precess about, which lies + * between the planet's equator and its orbit, and names it by its pole. The frame's x axis is + * where that plane rises through the ICRF equator, at right ascension 90 degrees past the pole's, + * which is what the node is counted from; its z axis is the pole, 90 degrees less its declination + * away from the celestial one. Read against the ecliptic instead, Io was up to 2.8 degrees from + * where Horizons has it between 1950 and 2100, Phobos 54 and Titan 127: their nodes are counted + * from a different line altogether. + */ +export function laplacePlaneToEquatorial(position: CartesianCoordinates, pole: { raDeg: number; decDeg: number }): CartesianCoordinates { + const tilt = (90 - pole.decDeg) * DEG_TO_RAD; + const node = (pole.raDeg + 90) * DEG_TO_RAD; + const y = position.y * Math.cos(tilt) - position.z * Math.sin(tilt); + const z = position.y * Math.sin(tilt) + position.z * Math.cos(tilt); + return { + x: position.x * Math.cos(node) - y * Math.sin(node), + y: position.x * Math.sin(node) + y * Math.cos(node), + z + }; +} + /** Inverse of {@link eclipticToEquatorial}. */ export function equatorialToEcliptic(position: CartesianCoordinates): CartesianCoordinates { const obliquity = OBLIQUITY_J2000_DEG * DEG_TO_RAD; diff --git a/src/app/shared/astro/kepler.spec.ts b/src/app/shared/astro/kepler.spec.ts index eaa1a51..62f918d 100644 --- a/src/app/shared/astro/kepler.spec.ts +++ b/src/app/shared/astro/kepler.spec.ts @@ -4,9 +4,11 @@ import { GM_SUN_AU3_PER_DAY2, DEFAULT_EPOCH_JD } from './constants'; import { gravitationalParameterFromPeriod, isPropagatableOrbit, + meanElementsAt, meanMotionRadPerDay, orbitEllipsePoints, orbitalPeriodDays, + positionAtEpoch, positionAtTrueAnomaly, propagateOrbit, resolveGravitationalParameter, @@ -131,6 +133,46 @@ describe('propagateOrbit', () => { }); }); +describe('meanElementsAt', () => { + /** A circle in the reference plane, prograde (0) or retrograde (180), whose node turns. */ + function circle(inclinationDeg: number) { + return { semiMajorAxisAu: 1, eccentricity: 0, inclinationDeg, longitudeOfAscendingNodeDeg: 0, argumentOfPeriapsisDeg: 0, meanAnomalyAtEpochDeg: 0, epochJd: DEFAULT_EPOCH_JD }; + } + const RATES = { meanMotionDegPerDay: 10, longitudeOfAscendingNodeDegPerDay: 0.5, argumentOfPeriapsisDegPerDay: 0.2 }; + + /** Longitude in the reference plane a day on, in degrees, signed. */ + function longitudeAfterOneDay(inclinationDeg: number): number { + const { x, y } = positionAtEpoch(meanElementsAt(circle(inclinationDeg), RATES, DEFAULT_EPOCH_JD + 1)); + return (Math.atan2(y, x) * 180) / Math.PI; + } + + it('goes round at its mean motion however its node and periapsis turn', () => { + expect(longitudeAfterOneDay(0)).toBeCloseTo(10, 9); + }); + + it('goes round a retrograde orbit backwards at the same rate, the node’s turning added back', () => { + // Taking the node off as for a prograde orbit made this 9 degrees, and Triton drifted a + // degree a year from where Horizons has it. + expect(longitudeAfterOneDay(180)).toBeCloseTo(-10, 9); + }); + + it('turns the node and periapsis at their own rates, and dates the result', () => { + const later = meanElementsAt(circle(0), RATES, DEFAULT_EPOCH_JD + 4); + expect(later.longitudeOfAscendingNodeDeg).toBeCloseTo(2, 12); + expect(later.argumentOfPeriapsisDeg).toBeCloseTo(0.8, 12); + expect(later.epochJd).toBe(DEFAULT_EPOCH_JD + 4); + }); + + it('adds Standish’s b T² + c cos(fT) + s sin(fT) to the mean anomaly', () => { + const terms = { b: -0.00012452, c: 0.0606406, s: -0.35635438, f: 38.35125 }; + const T = 0.7; + const withTerms = meanElementsAt(circle(0), { ...RATES, meanAnomalyTerms: terms }, DEFAULT_EPOCH_JD + T * 36525); + const without = meanElementsAt(circle(0), RATES, DEFAULT_EPOCH_JD + T * 36525); + const f = (terms.f * T * Math.PI) / 180; + expect(withTerms.meanAnomalyAtEpochDeg - without.meanAnomalyAtEpochDeg).toBeCloseTo(terms.b * T * T + terms.c * Math.cos(f) + terms.s * Math.sin(f), 9); + }); +}); + describe('orbitEllipsePoints', () => { it('samples a closed loop whose distances stay within the periapsis/apoapsis bounds', () => { const elements = resolveOrbitalElements({ semiMajorAxisAu: 5, eccentricity: 0.4 }); diff --git a/src/app/shared/astro/kepler.ts b/src/app/shared/astro/kepler.ts index 05dff65..7354934 100644 --- a/src/app/shared/astro/kepler.ts +++ b/src/app/shared/astro/kepler.ts @@ -1,9 +1,10 @@ import { CartesianCoordinates } from './coordinates'; import { DEFAULT_EPOCH_JD, GM_SUN_AU3_PER_DAY2 } from './constants'; -import { OrbitalElements } from '../models/body.model'; +import { MeanElementRates, OrbitalElements } from '../models/body.model'; const DEG_TO_RAD = Math.PI / 180; const TWO_PI = Math.PI * 2; +const DAYS_PER_JULIAN_CENTURY = 36525; /** * Fills in the elements the Kepler propagator needs but that some sources (e.g. exoplanets, @@ -208,17 +209,63 @@ export function positionAtTrueAnomaly(elements: OrbitalElements, trueAnomalyRad: } /** - * Propagates `elements` to Julian date `epochJdEval`, returning the body's position (AU) - * relative to its central body. This is the app's "current epoch" evaluation used for live - * (and future time-scrubbable) positions, as opposed to {@link orbitEllipsePoints} which + * The rates of an orbit that only goes round: Kepler's mean motion from the central mass, with + * nothing turning. What an exoplanet has, since the archive publishes no precession. + */ +export function keplerRates(semiMajorAxisAu: number, gmAu3PerDay2: number): MeanElementRates { + return { + meanMotionDegPerDay: meanMotionRadPerDay(semiMajorAxisAu, gmAu3PerDay2) / DEG_TO_RAD, + longitudeOfAscendingNodeDegPerDay: 0, + argumentOfPeriapsisDegPerDay: 0 + }; +} + +/** + * The elements at `epochJdEval`, each moved from its epoch at its own rate, and returned with that + * date as their epoch — so {@link positionAtEpoch} places the body, and the node and periapsis + * say where to draw the orbit it is on. + * + * The mean anomaly is what is left of the body's motion once the node and periapsis have turned: + * `meanMotionDegPerDay` is how fast it goes round in space, and a periapsis that has moved on is + * that much further to reach. On a retrograde orbit, past 90 degrees, the body runs against the + * direction the node is counted in, so the node's turning is added back rather than taken off. + * Taken off, Triton — whose node turns half a degree a year — drifted a degree a year from where + * Horizons has it, 105 degrees by 2100. + */ +export function meanElementsAt(elements: OrbitalElements, rates: MeanElementRates, epochJdEval: number): OrbitalElements { + const days = epochJdEval - elements.epochJd; + const node = rates.longitudeOfAscendingNodeDegPerDay * days; + const periapsis = rates.argumentOfPeriapsisDegPerDay * days; + const nodeAlongOrbit = elements.inclinationDeg > 90 ? -node : node; + const terms = rates.meanAnomalyTerms; + const centuries = days / DAYS_PER_JULIAN_CENTURY; + const extra = terms + ? terms.b * centuries * centuries + terms.c * Math.cos(terms.f * centuries * DEG_TO_RAD) + terms.s * Math.sin(terms.f * centuries * DEG_TO_RAD) + : 0; + return { + semiMajorAxisAu: elements.semiMajorAxisAu + (rates.semiMajorAxisAuPerDay ?? 0) * days, + eccentricity: elements.eccentricity + (rates.eccentricityPerDay ?? 0) * days, + inclinationDeg: elements.inclinationDeg + (rates.inclinationDegPerDay ?? 0) * days, + longitudeOfAscendingNodeDeg: elements.longitudeOfAscendingNodeDeg + node, + argumentOfPeriapsisDeg: elements.argumentOfPeriapsisDeg + periapsis, + meanAnomalyAtEpochDeg: elements.meanAnomalyAtEpochDeg + rates.meanMotionDegPerDay * days - periapsis - nodeAlongOrbit + extra, + epochJd: epochJdEval + }; +} + +/** Where `elements` put the body at their own epoch (AU, relative to the central body). */ +export function positionAtEpoch(elements: OrbitalElements): CartesianCoordinates { + const eccentricAnomalyRad = solveEccentricAnomaly(elements.meanAnomalyAtEpochDeg * DEG_TO_RAD, elements.eccentricity); + return positionAtTrueAnomaly(elements, trueAnomalyFromEccentricAnomaly(eccentricAnomalyRad, elements.eccentricity)); +} + +/** + * Propagates `elements` to Julian date `epochJdEval` around a central mass, returning the body's + * position (AU) relative to its central body, as opposed to {@link orbitEllipsePoints} which * samples the fixed orbit shape independent of time. */ export function propagateOrbit(elements: OrbitalElements, gmAu3PerDay2: number, epochJdEval: number): CartesianCoordinates { - const meanMotion = meanMotionRadPerDay(elements.semiMajorAxisAu, gmAu3PerDay2); - const meanAnomalyRad = elements.meanAnomalyAtEpochDeg * DEG_TO_RAD + meanMotion * (epochJdEval - elements.epochJd); - const eccentricAnomalyRad = solveEccentricAnomaly(meanAnomalyRad, elements.eccentricity); - const trueAnomalyRad = trueAnomalyFromEccentricAnomaly(eccentricAnomalyRad, elements.eccentricity); - return positionAtTrueAnomaly(elements, trueAnomalyRad); + return positionAtEpoch(meanElementsAt(elements, keplerRates(elements.semiMajorAxisAu, gmAu3PerDay2), epochJdEval)); } /** diff --git a/src/app/shared/astro/mean-elements.spec.ts b/src/app/shared/astro/mean-elements.spec.ts new file mode 100644 index 0000000..32b61a7 --- /dev/null +++ b/src/app/shared/astro/mean-elements.spec.ts @@ -0,0 +1,110 @@ +import { describe, expect, it } from 'vitest'; + +import { parsePlanetMeanElements, parseSatelliteMeanElements } from './mean-elements'; + +/** Standish's p_elem_t2.txt, cut to the lines that matter here, as JPL published them. */ +const TABLE_2 = `Keplerian elements and their rates, with respect to the mean ecliptic and equinox of J2000, +valid for the time-interval 3000 BC -- 3000 AD. NOTE: the computation of M for Jupiter through +Pluto *must* be augmented by the additional terms given in Table 2b (below). + +EM Bary 1.00000018 0.01673163 -0.00054346 100.46691572 102.93005885 -5.11260389 + -0.00000003 -0.00003661 -0.01337178 35999.37306329 0.31795260 -0.24123856 +Jupiter 5.20248019 0.04853590 1.29861416 34.33479152 14.27495244 100.29282654 + -0.00002864 0.00018026 -0.00322699 3034.90371757 0.18199196 0.13024619 +Pluto 39.48686035 0.24885238 17.14104260 238.96535011 224.09702598 110.30167986 + 0.00449751 0.00006016 0.00000501 145.18042903 -0.00968827 -0.00809981 + +Table 2b. +Jupiter -0.00012452 0.06064060 -0.35635438 38.35125000 +Pluto -0.01262724 +`; + +/** The satellite page's markup around three rows, as the 2021 page served it. */ +const SATELLITES = ` +Satellites of Earth +jump to: Earth, Mars +

Mean ecliptic orbital elements

+Epoch 2000 Jan. 1.50 TT
+Moon +384400.0.0554318.15135.275.16125.08 +13.17635827.3225.99718.600 +1 +Satellites of Jupiter +jump to: Earth, Mars +

Mean orbital elements referred to the local Laplace planes

+Epoch 1997 Jan. 16.00 TT
+Io421800.0.0041 +84.129342.0210.03643.977203.4889583 +1.7691.6257.420268.05764.495 +0.000 +11 +Satellites of Neptune +jump to: Earth, Mars +

Mean orbital elements referred to the local Laplace planes

+Epoch 2000 Jan. 1.50 TT
+Triton354759.0.0000 +66.142352.257156.865177.608 +61.25726385.877386.371687.446 +299.45643.4140.010 +54 +`; + +describe('parsePlanetMeanElements', () => { + it('turns Standish’s longitudes into the argument of periapsis and mean anomaly', () => { + const { orbit } = parsePlanetMeanElements(TABLE_2, 'jupiter'); + expect(orbit.argumentOfPeriapsisDeg).toBeCloseTo(14.27495244 - 100.29282654, 8); + expect(orbit.meanAnomalyAtEpochDeg).toBeCloseTo(34.33479152 - 14.27495244, 8); + expect(orbit.epochJd).toBe(2451545); + }); + + it('gives the rates per day, the mean motion being the mean longitude’s', () => { + const { rates } = parsePlanetMeanElements(TABLE_2, 'earth'); + // 35 999.373 degrees a century is the sidereal year. + expect(360 / rates.meanMotionDegPerDay).toBeCloseTo(365.2564, 4); + expect(rates.argumentOfPeriapsisDegPerDay * 36525).toBeCloseTo(0.3179526 + 0.24123856, 8); + }); + + it('carries Table 2b’s terms for Jupiter and beyond, and none for the inner planets', () => { + expect(parsePlanetMeanElements(TABLE_2, 'jupiter').rates.meanAnomalyTerms).toEqual({ b: -0.00012452, c: 0.0606406, s: -0.35635438, f: 38.35125 }); + expect(parsePlanetMeanElements(TABLE_2, 'earth').rates.meanAnomalyTerms).toBeUndefined(); + }); + + it('reads Pluto’s row, not the note above the table that starts a line with its name', () => { + const pluto = parsePlanetMeanElements(TABLE_2, 'pluto'); + expect(pluto.orbit.semiMajorAxisAu).toBe(39.48686035); + expect(pluto.rates.meanAnomalyTerms).toEqual({ b: -0.01262724, c: 0, s: 0, f: 0 }); + }); +}); + +describe('parseSatelliteMeanElements', () => { + it('reads a Laplace-plane row with its pole and its section’s epoch', () => { + const io = parseSatelliteMeanElements(SATELLITES, 'Jupiter', 'Io', true); + expect(io.laplacePole).toEqual({ raDeg: 268.057, decDeg: 64.495 }); + expect(io.orbit.epochJd).toBe(2450464.5); + expect(io.rates.meanMotionDegPerDay).toBe(203.4889583); + expect(io.orbitSource).toBe('JPL SSD satellite mean elements, epoch 1997 Jan 16'); + }); + + it('reads the Moon against the ecliptic, with no pole', () => { + const moon = parseSatelliteMeanElements(SATELLITES, 'Earth', 'Moon', false); + expect(moon.laplacePole).toBeUndefined(); + expect(moon.orbit.epochJd).toBe(2451545); + expect(moon.orbit.semiMajorAxisAu * 149597870.7).toBeCloseTo(384400, 3); + }); + + it('regresses a prograde node and advances a periapsis, as the planet’s oblateness turns them', () => { + const { rates } = parseSatelliteMeanElements(SATELLITES, 'Earth', 'Moon', false); + expect(rates.longitudeOfAscendingNodeDegPerDay).toBeCloseTo(-360 / (18.6 * 365.25), 9); + expect(rates.argumentOfPeriapsisDegPerDay).toBeCloseTo(360 / (5.997 * 365.25), 9); + }); + + it('advances the node of a retrograde orbit', () => { + const { rates } = parseSatelliteMeanElements(SATELLITES, 'Neptune', 'Triton', false); + expect(rates.longitudeOfAscendingNodeDegPerDay).toBeCloseTo(360 / (687.446 * 365.25), 9); + }); + + it('turns the periapsis backwards where a resonance holds it', () => { + const { rates } = parseSatelliteMeanElements(SATELLITES, 'Jupiter', 'Io', true); + expect(rates.argumentOfPeriapsisDegPerDay).toBeCloseTo(-360 / (1.625 * 365.25), 9); + }); +}); diff --git a/src/app/shared/astro/mean-elements.ts b/src/app/shared/astro/mean-elements.ts new file mode 100644 index 0000000..bbd8251 --- /dev/null +++ b/src/app/shared/astro/mean-elements.ts @@ -0,0 +1,151 @@ +import { MeanElementRates, OrbitalElements } from '../models/body.model'; + +/** + * Reads JPL's two tables of mean orbital elements, which the ETL fetches (see + * `tools/etl/lib/mean-elements.ts`), into the elements and rates `bodies.json` carries. + */ + +const KM_PER_AU = 149597870.7; +const J2000_JD = 2451545.0; +const DAYS_PER_JULIAN_CENTURY = 36525; +const DAYS_PER_JULIAN_YEAR = 365.25; + +const PLANET_ORBIT_SOURCE = 'JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000'; + +export interface MeanOrbit { + orbit: OrbitalElements; + rates: MeanElementRates; + laplacePole?: { raDeg: number; decDeg: number }; + orbitSource: string; +} + +/** Table 2a's name for each planet; Earth's row is the Earth-Moon barycentre, 4 700 km off Earth. */ +const PLANET_ROW_NAMES: Record = { + mercury: 'Mercury', + venus: 'Venus', + earth: 'EM Bary', + mars: 'Mars', + jupiter: 'Jupiter', + saturn: 'Saturn', + uranus: 'Uranus', + neptune: 'Neptune', + pluto: 'Pluto' +}; + +function numbers(text: string): number[] { + return text.trim().split(/\s+/).map(Number); +} + +/** + * Reads one planet's row pair from Table 2a, and its Table 2b terms where it has them. The + * elements are Standish's own — a, e, I, mean longitude L, longitude of perihelion ϖ, node Ω — + * turned into the argument of periapsis ϖ - Ω and mean anomaly L - ϖ the propagator takes. + */ +export function parsePlanetMeanElements(text: string, bodyId: string): MeanOrbit { + const name = PLANET_ROW_NAMES[bodyId]; + const lines = text.split(/\r?\n/); + // A row is the name followed by a number: the notes above the table start a line with "Pluto" too. + const rows = lines.flatMap((line, index) => (name && new RegExp(`^${name}\\s+-?[\\d.]`).test(line) ? [index] : [])); + if (rows.length === 0) { + throw new Error(`No row for ${bodyId} in Standish's Table 2a.`); + } + const values = [...numbers(lines[rows[0]].slice(name.length)), ...numbers(lines[rows[0] + 1])]; + if (values.length !== 12 || !values.every(Number.isFinite)) { + throw new Error(`Standish's Table 2a rows for ${name} did not parse: ${values.join(' ')}`); + } + const [a, e, inclination, meanLongitude, perihelion, node, aRate, eRate, inclinationRate, meanLongitudeRate, perihelionRate, nodeRate] = values; + // Table 2b repeats the name further down, with b, c, s, f (Pluto has b alone). + const extra = rows[1] === undefined ? undefined : numbers(lines[rows[1]].slice(name.length)); + const [b = 0, c = 0, s = 0, f = 0] = extra ?? []; + + return { + orbit: { + semiMajorAxisAu: a, + eccentricity: e, + inclinationDeg: inclination, + longitudeOfAscendingNodeDeg: node, + argumentOfPeriapsisDeg: perihelion - node, + meanAnomalyAtEpochDeg: meanLongitude - perihelion, + epochJd: J2000_JD + }, + rates: { + meanMotionDegPerDay: meanLongitudeRate / DAYS_PER_JULIAN_CENTURY, + longitudeOfAscendingNodeDegPerDay: nodeRate / DAYS_PER_JULIAN_CENTURY, + argumentOfPeriapsisDegPerDay: (perihelionRate - nodeRate) / DAYS_PER_JULIAN_CENTURY, + semiMajorAxisAuPerDay: aRate / DAYS_PER_JULIAN_CENTURY, + eccentricityPerDay: eRate / DAYS_PER_JULIAN_CENTURY, + inclinationDegPerDay: inclinationRate / DAYS_PER_JULIAN_CENTURY, + ...(extra ? { meanAnomalyTerms: { b, c, s, f } } : {}) + }, + orbitSource: PLANET_ORBIT_SOURCE + }; +} + +const MONTHS = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']; + +/** `1997 Jan. 16.00` as a Julian date. TT and TDB differ by under two milliseconds. */ +function julianDate(year: number, month: string, day: number): number { + const monthIndex = MONTHS.indexOf(month); + if (monthIndex < 0) { + throw new Error(`Unknown month ${month}.`); + } + return Date.UTC(year, monthIndex, 1) / 86400000 + 2440587.5 + day - 1; +} + +/** + * Reads one moon's row from the satellite page: `a e w M i node n P Pw Pnode`, then the Laplace + * pole `RA Dec Tilt` where the section is referred to one, then a reference number. + * + * The page gives the two precession periods as magnitudes, so their sense is supplied here. A + * node driven by the planet's oblateness regresses on a prograde orbit and advances on a + * retrograde one, and the orbit's inclination says which. A periapsis advances — except where a + * resonance forces the eccentricity, which `apsidesRegress` names: Io's and Europa's are held to + * the line of their conjunctions, which turns backwards at 2 n(Europa) - n(Io) = 0.74 degrees a + * day, and that is exactly the 1.625- and 1.394-year periods the table gives for them. Read as + * advancing, Io was 0.9 degrees out and Europa 2.1. + */ +export function parseSatelliteMeanElements(html: string, planetName: string, moonName: string, apsidesRegress: boolean): MeanOrbit { + const text = html.replace(/<[^>]+>/g, ' ').replace(/ /g, ' ').replace(/\s+/g, ' '); + const section = text.indexOf(`Satellites of ${planetName} jump to`); + if (section < 0) { + throw new Error(`No section for the satellites of ${planetName}.`); + } + const row = text.slice(section).match(new RegExp(` ${moonName} ((?:-?[\\d.]+ )+)`)); + if (!row || row.index === undefined) { + throw new Error(`No row for ${moonName} among the satellites of ${planetName}.`); + } + const before = text.slice(section, section + row.index); + const epoch = [...before.matchAll(/Epoch (\d{4}) (\w{3})\. ([\d.]+) T/g)].at(-1); + if (!epoch) { + throw new Error(`No epoch above ${moonName}'s row.`); + } + const laplace = before.lastIndexOf('Laplace plane') > before.lastIndexOf('Mean ecliptic'); + const values = numbers(row[1]); + const expected = laplace ? 14 : 11; + if (values.length !== expected || !values.every(Number.isFinite)) { + throw new Error(`${moonName}'s row has ${values.length} numbers, ${expected} expected: ${row[1]}`); + } + const [aKm, e, periapsis, meanAnomaly, inclination, node, meanMotion, , periapsisPeriodYears, nodePeriodYears, raDeg, decDeg] = values; + const nodeSense = inclination > 90 ? 1 : -1; + const periapsisSense = apsidesRegress ? -1 : 1; + const perDay = (periodYears: number): number => (periodYears > 0 ? 360 / (periodYears * DAYS_PER_JULIAN_YEAR) : 0); + + return { + orbit: { + semiMajorAxisAu: aKm / KM_PER_AU, + eccentricity: e, + inclinationDeg: inclination, + longitudeOfAscendingNodeDeg: node, + argumentOfPeriapsisDeg: periapsis, + meanAnomalyAtEpochDeg: meanAnomaly, + epochJd: julianDate(Number(epoch[1]), epoch[2], Number(epoch[3])) + }, + rates: { + meanMotionDegPerDay: meanMotion, + longitudeOfAscendingNodeDegPerDay: nodeSense * perDay(nodePeriodYears), + argumentOfPeriapsisDegPerDay: periapsisSense * perDay(periapsisPeriodYears) + }, + ...(laplace ? { laplacePole: { raDeg, decDeg } } : {}), + orbitSource: `JPL SSD satellite mean elements, epoch ${epoch[1]} ${epoch[2]} ${Math.floor(Number(epoch[3]))}` + }; +} diff --git a/src/app/shared/models/body.model.ts b/src/app/shared/models/body.model.ts index 5bf97a5..d33798b 100644 --- a/src/app/shared/models/body.model.ts +++ b/src/app/shared/models/body.model.ts @@ -1,7 +1,7 @@ /** - * Osculating Keplerian orbital elements at a reference epoch. Positions are derived - * client-side by propagating these elements forward/backward from `epochJd` (see - * `shared/astro/kepler.ts`), rather than fetching per-frame positions. + * Keplerian orbital elements at a reference epoch. Positions are derived client-side by + * propagating these elements forward/backward from `epochJd` (see `shared/astro/kepler.ts`), + * rather than fetching per-frame positions. */ export interface OrbitalElements { semiMajorAxisAu: number; @@ -14,8 +14,36 @@ export interface OrbitalElements { } /** - * A solar-system planet, moon, or dwarf planet, sourced from JPL Horizons/SSD orbital - * elements. `systemStarId` links back to the HYG star index (the Sun, see `SUN_STAR_ID`). + * How a body's mean elements move away from their epoch, per day. + * + * Mean elements rather than one osculating set, because the map's clock runs decades in minutes. + * An osculating orbit is exact at its instant and drifts from then on: fed to Kepler with a mass + * ratio, the Moon's went round in 27.70 days instead of 27.32 and was 66 degrees out after a year. + * A mean set carries its own measured motion, and the slow turning of its node and periapsis, so + * it holds for as long as its source was fit over. + */ +export interface MeanElementRates { + /** + * How fast the body goes round in space, in degrees per day: the rate of its mean longitude. + * 360 over this is its sidereal period. + */ + meanMotionDegPerDay: number; + longitudeOfAscendingNodeDegPerDay: number; + argumentOfPeriapsisDegPerDay: number; + semiMajorAxisAuPerDay?: number; + eccentricityPerDay?: number; + inclinationDegPerDay?: number; + /** + * Standish's extra terms in the mean anomaly of Jupiter and beyond, `b T² + c cos(f T) + + * s sin(f T)` degrees, with T in Julian centuries from the epoch and f in degrees per century: + * the great-inequality wobble his 3000 BC to AD 3000 fit needs on top of its linear rates. + */ + meanAnomalyTerms?: { b: number; c: number; s: number; f: number }; +} + +/** + * A solar-system planet, moon, or dwarf planet: JPL mean orbital elements, and JPL Horizons + * physical data. `systemStarId` links back to the HYG star index (the Sun, see `SUN_STAR_ID`). */ export interface BodyRecord { id: string; @@ -23,7 +51,18 @@ export interface BodyRecord { name: string; kind: 'planet' | 'moon' | 'dwarf'; radiusKm: number; + /** Mean elements at `orbit.epochJd`, moving at `rates`. */ orbit: OrbitalElements; + rates: MeanElementRates; + /** + * The pole of the plane a moon's elements are measured against, where that is its local + * Laplace plane: right ascension and declination in the ICRF. The node is then counted from + * where that plane crosses the ICRF equator. Absent means the J2000 ecliptic, as for the + * planets and the Moon. + */ + laplacePole?: { raDeg: number; decDeg: number }; + /** Where the elements come from and the span they hold over, as the card prints it. */ + orbitSource: string; /** * For `kind: 'moon'`, the `id` of the planet it orbits — its `orbit` is expressed * relative to that planet, not heliocentrically. Undefined for planets/dwarfs. diff --git a/src/assets/data/bodies.json b/src/assets/data/bodies.json index b8ccae8..ea8f8e6 100644 --- a/src/assets/data/bodies.json +++ b/src/assets/data/bodies.json @@ -6,14 +6,23 @@ "kind": "planet", "radiusKm": 2439.4, "orbit": { - "semiMajorAxisAu": 0.38709857292289906, - "eccentricity": 0.2056388859154261, - "inclinationDeg": 7.003502017708398, - "longitudeOfAscendingNodeDeg": 48.29980532881731, - "argumentOfPeriapsisDeg": 29.19561108948743, - "meanAnomalyAtEpochDeg": 103.9465145977117, - "epochJd": 2460676.5 + "semiMajorAxisAu": 0.38709843, + "eccentricity": 0.20563661, + "inclinationDeg": 7.00559432, + "longitudeOfAscendingNodeDeg": 48.33961819, + "argumentOfPeriapsisDeg": 29.118100759999997, + "meanAnomalyAtEpochDeg": 174.79394829, + "epochJd": 2451545 }, + "rates": { + "meanMotionDegPerDay": 4.092338805372484, + "longitudeOfAscendingNodeDegPerDay": -0.0000033440607802874744, + "argumentOfPeriapsisDegPerDay": 0.000007708198494182067, + "semiMajorAxisAuPerDay": 0, + "eccentricityPerDay": 5.812457221081451e-10, + "inclinationDegPerDay": -1.6157645448323065e-7 + }, + "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 1407.512239412851, "obliquityDeg": 2.11 }, @@ -24,14 +33,23 @@ "kind": "planet", "radiusKm": 6051.84, "orbit": { - "semiMajorAxisAu": 0.7233281694940689, - "eccentricity": 0.006746576710187382, - "inclinationDeg": 3.394393253252075, - "longitudeOfAscendingNodeDeg": 76.61185393039594, - "argumentOfPeriapsisDeg": 55.15075425343929, - "meanAnomalyAtEpochDeg": 280.0749102629981, - "epochJd": 2460676.5 + "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 + }, + "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": -5832.539941165383, "obliquityDeg": 177.3 }, @@ -42,14 +60,23 @@ "kind": "planet", "radiusKm": 6371.01, "orbit": { - "semiMajorAxisAu": 1.0009125494279245, - "eccentricity": 0.01756190256786036, - "inclinationDeg": 0.002977082685807642, - "longitudeOfAscendingNodeDeg": 190.1716375775361, - "argumentOfPeriapsisDeg": 272.9783142442708, - "meanAnomalyAtEpochDeg": 357.4122246804211, - "epochJd": 2460676.5 + "semiMajorAxisAu": 1.00000018, + "eccentricity": 0.01673163, + "inclinationDeg": -0.00054346, + "longitudeOfAscendingNodeDeg": -5.11260389, + "argumentOfPeriapsisDeg": 108.04266274, + "meanAnomalyAtEpochDeg": -2.4631431299999917, + "epochJd": 2451545 }, + "rates": { + "meanMotionDegPerDay": 0.9856091187759068, + "longitudeOfAscendingNodeDegPerDay": -0.000006604751813826146, + "argumentOfPeriapsisDegPerDay": 0.000015309819575633124, + "semiMajorAxisAuPerDay": -8.213552361396303e-13, + "eccentricityPerDay": -1.002327173169062e-9, + "inclinationDegPerDay": -3.6609938398357287e-7 + }, + "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 23.934472399219285, "obliquityDeg": 23.4392911 }, @@ -60,14 +87,23 @@ "kind": "planet", "radiusKm": 3389.92, "orbit": { - "semiMajorAxisAu": 1.5237367794976557, - "eccentricity": 0.09343026241683372, - "inclinationDeg": 1.847583389630714, - "longitudeOfAscendingNodeDeg": 49.48673257018439, - "argumentOfPeriapsisDeg": 286.7114828288203, - "meanAnomalyAtEpochDeg": 124.444888195349, - "epochJd": 2460676.5 + "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 + }, + "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 24.622955438662025, "obliquityDeg": 25.19 }, @@ -78,14 +114,29 @@ "kind": "planet", "radiusKm": 69911, "orbit": { - "semiMajorAxisAu": 5.20282710471751, - "eccentricity": 0.04830624138918495, - "inclinationDeg": 1.303459692430419, - "longitudeOfAscendingNodeDeg": 100.5202095061913, - "argumentOfPeriapsisDeg": 273.6090683600047, - "meanAnomalyAtEpochDeg": 58.98282567286461, - "epochJd": 2460676.5 + "semiMajorAxisAu": 5.20248019, + "eccentricity": 0.0485359, + "inclinationDeg": 1.29861416, + "longitudeOfAscendingNodeDeg": 100.29282654, + "argumentOfPeriapsisDeg": -86.0178741, + "meanAnomalyAtEpochDeg": 20.059839080000003, + "epochJd": 2451545 }, + "rates": { + "meanMotionDegPerDay": 0.08309113532019165, + "longitudeOfAscendingNodeDegPerDay": 0.0000035659463381245725, + "argumentOfPeriapsisDegPerDay": 0.0000014167219712525667, + "semiMajorAxisAuPerDay": -7.841204654346339e-10, + "eccentricityPerDay": 4.935249828884326e-9, + "inclinationDegPerDay": -8.83501711156742e-8, + "meanAnomalyTerms": { + "b": -0.00012452, + "c": 0.0606406, + "s": -0.35635438, + "f": 38.35125 + } + }, + "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 9.925102371306965, "obliquityDeg": 3.13 }, @@ -96,14 +147,29 @@ "kind": "planet", "radiusKm": 58232, "orbit": { - "semiMajorAxisAu": 9.555678383881443, - "eccentricity": 0.05522302318620472, - "inclinationDeg": 2.485819253218824, - "longitudeOfAscendingNodeDeg": 113.5595556634874, - "argumentOfPeriapsisDeg": 337.1663598259115, - "meanAnomalyAtEpochDeg": 264.9877650588944, - "epochJd": 2460676.5 + "semiMajorAxisAu": 9.54149883, + "eccentricity": 0.05550825, + "inclinationDeg": 2.49424102, + "longitudeOfAscendingNodeDeg": 113.63998702, + "argumentOfPeriapsisDeg": -20.778626390000014, + "meanAnomalyAtEpochDeg": -42.78564733999999, + "epochJd": 2451545 }, + "rates": { + "meanMotionDegPerDay": 0.033459683702669406, + "longitudeOfAscendingNodeDegPerDay": -0.000006848734291581108, + "argumentOfPeriapsisDegPerDay": 0.000021682266940451745, + "semiMajorAxisAuPerDay": -8.391512662559891e-10, + "eccentricityPerDay": -8.773169062286106e-9, + "inclinationDegPerDay": 1.2374236824093085e-7, + "meanAnomalyTerms": { + "b": 0.00025899, + "c": -0.13434469, + "s": 0.87320147, + "f": 38.35125 + } + }, + "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 10.656221583138441, "obliquityDeg": 26.73 }, @@ -114,14 +180,29 @@ "kind": "planet", "radiusKm": 25362, "orbit": { - "semiMajorAxisAu": 19.30135052232751, - "eccentricity": 0.04562415515293685, - "inclinationDeg": 0.7728969735590447, - "longitudeOfAscendingNodeDeg": 74.0126541901479, - "argumentOfPeriapsisDeg": 90.49593456147204, - "meanAnomalyAtEpochDeg": 255.8822481742851, - "epochJd": 2460676.5 + "semiMajorAxisAu": 19.18797948, + "eccentricity": 0.0468574, + "inclinationDeg": 0.77298127, + "longitudeOfAscendingNodeDeg": 73.96250215, + "argumentOfPeriapsisDeg": 98.47154226, + "meanAnomalyAtEpochDeg": 141.76872184, + "epochJd": 2451545 }, + "rates": { + "meanMotionDegPerDay": 0.011731557178644764, + "longitudeOfAscendingNodeDegPerDay": 0.0000015714439425051334, + "argumentOfPeriapsisDegPerDay": 9.65718275154004e-7, + "semiMajorAxisAuPerDay": -5.600273785078713e-9, + "eccentricityPerDay": -4.2436687200547574e-10, + "inclinationDegPerDay": -4.932375085557837e-8, + "meanAnomalyTerms": { + "b": 0.00058331, + "c": -0.97731848, + "s": 0.17689245, + "f": 7.67025 + } + }, + "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": -17.24003330792427, "obliquityDeg": 97.77 }, @@ -132,14 +213,29 @@ "kind": "planet", "radiusKm": 24624, "orbit": { - "semiMajorAxisAu": 30.183632709237763, - "eccentricity": 0.01266098870528033, - "inclinationDeg": 1.774849832507214, - "longitudeOfAscendingNodeDeg": 131.9473970418126, - "argumentOfPeriapsisDeg": 268.1146665565159, - "meanAnomalyAtEpochDeg": 319.6858384317641, - "epochJd": 2460676.5 + "semiMajorAxisAu": 30.06952752, + "eccentricity": 0.00895439, + "inclinationDeg": 1.7700552, + "longitudeOfAscendingNodeDeg": 131.78635853, + "argumentOfPeriapsisDeg": -85.10477129, + "meanAnomalyAtEpochDeg": 257.54130563, + "epochJd": 2451545 }, + "rates": { + "meanMotionDegPerDay": 0.005981249914852841, + "longitudeOfAscendingNodeDegPerDay": -1.6599644079397672e-7, + "argumentOfPeriapsisDegPerDay": 4.4250239561943875e-7, + "semiMajorAxisAuPerDay": 1.7650924024640657e-9, + "eccentricityPerDay": 2.2395619438740589e-10, + "inclinationDegPerDay": 6.132785763175907e-9, + "meanAnomalyTerms": { + "b": -0.00041348, + "c": 0.68346318, + "s": -0.10162547, + "f": 7.67025 + } + }, + "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 16.110037586020876, "obliquityDeg": 28.32 }, @@ -150,14 +246,29 @@ "kind": "dwarf", "radiusKm": 1188.3, "orbit": { - "semiMajorAxisAu": 39.28778257358678, - "eccentricity": 0.2438605399689669, - "inclinationDeg": 16.93906659887321, - "longitudeOfAscendingNodeDeg": 110.1714547539489, - "argumentOfPeriapsisDeg": 113.5754868232679, - "meanAnomalyAtEpochDeg": 51.93655343727463, - "epochJd": 2460676.5 + "semiMajorAxisAu": 39.48686035, + "eccentricity": 0.24885238, + "inclinationDeg": 17.1410426, + "longitudeOfAscendingNodeDeg": 110.30167986, + "argumentOfPeriapsisDeg": 113.79534612000002, + "meanAnomalyAtEpochDeg": 14.86832412999999, + "epochJd": 2451545 }, + "rates": { + "meanMotionDegPerDay": 0.003974823518959616, + "longitudeOfAscendingNodeDegPerDay": -2.2176071184120468e-7, + "argumentOfPeriapsisDegPerDay": -4.3489664613278575e-8, + "semiMajorAxisAuPerDay": 1.2313511293634495e-7, + "eccentricityPerDay": 1.6470910335386722e-9, + "inclinationDegPerDay": 1.3716632443531827e-10, + "meanAnomalyTerms": { + "b": -0.01262724, + "c": 0, + "s": 0, + "f": 0 + } + }, + "orbitSource": "JPL approximate mean elements (Standish), fit for 3000 BC to AD 3000", "rotationPeriodHours": 153.29335198, "obliquityDeg": 119.6 }, @@ -168,16 +279,22 @@ "kind": "moon", "radiusKm": 1737.53, "orbit": { - "semiMajorAxisAu": 0.0025850155421324465, - "eccentricity": 0.04034704696474715, - "inclinationDeg": 5.004154175941119, - "longitudeOfAscendingNodeDeg": 0.531732989166472, - "argumentOfPeriapsisDeg": 6.557682814330136, - "meanAnomalyAtEpochDeg": 290.7825171697369, - "epochJd": 2460676.5 + "semiMajorAxisAu": 0.0025695552897999907, + "eccentricity": 0.0554, + "inclinationDeg": 5.16, + "longitudeOfAscendingNodeDeg": 125.08, + "argumentOfPeriapsisDeg": 318.15, + "meanAnomalyAtEpochDeg": 135.27, + "epochJd": 2451545 }, + "rates": { + "meanMotionDegPerDay": 13.176358, + "longitudeOfAscendingNodeDegPerDay": -0.052990660396105185, + "argumentOfPeriapsisDegPerDay": 0.164353223839846 + }, + "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "earth", - "rotationPeriodHours": 664.8546956215548, + "rotationPeriodHours": 655.7198886065481, "obliquityDeg": 6.67 }, { @@ -187,16 +304,26 @@ "kind": "moon", "radiusKm": 13.1, "orbit": { - "semiMajorAxisAu": 0.00006269319601715431, - "eccentricity": 0.01559367482014246, - "inclinationDeg": 26.37766999967838, - "longitudeOfAscendingNodeDeg": 85.13091230996636, - "argumentOfPeriapsisDeg": 356.4060881660871, - "meanAnomalyAtEpochDeg": 342.6005509941174, - "epochJd": 2460676.5 + "semiMajorAxisAu": 0.00006267468885838895, + "eccentricity": 0.0151, + "inclinationDeg": 1.075, + "longitudeOfAscendingNodeDeg": 207.784, + "argumentOfPeriapsisDeg": 150.057, + "meanAnomalyAtEpochDeg": 91.059, + "epochJd": 2433282.5 }, + "rates": { + "meanMotionDegPerDay": 1128.8447569, + "longitudeOfAscendingNodeDegPerDay": -0.43579001784832494, + "argumentOfPeriapsisDegPerDay": 0.871002371303956 + }, + "laplacePole": { + "raDeg": 317.671, + "decDeg": 52.893 + }, + "orbitSource": "JPL SSD satellite mean elements, epoch 1950 Jan 1", "parentBodyId": "mars", - "rotationPeriodHours": 7.660212212136228 + "rotationPeriodHours": 7.653842521027348 }, { "id": "deimos", @@ -205,16 +332,26 @@ "kind": "moon", "radiusKm": 7.8, "orbit": { - "semiMajorAxisAu": 0.00015681797919999554, - "eccentricity": 0.0002637395029404219, - "inclinationDeg": 24.26157855280465, - "longitudeOfAscendingNodeDeg": 80.74605390828611, - "argumentOfPeriapsisDeg": 21.65246201037153, - "meanAnomalyAtEpochDeg": 273.8943716897566, - "epochJd": 2460676.5 + "semiMajorAxisAu": 0.00015680704471417322, + "eccentricity": 0.0002, + "inclinationDeg": 1.788, + "longitudeOfAscendingNodeDeg": 24.525, + "argumentOfPeriapsisDeg": 260.729, + "meanAnomalyAtEpochDeg": 325.329, + "epochJd": 2433282.5 }, + "rates": { + "meanMotionDegPerDay": 285.161879, + "longitudeOfAscendingNodeDegPerDay": -0.018072715865968356, + "argumentOfPeriapsisDegPerDay": 0.03601079576648982 + }, + "laplacePole": { + "raDeg": 316.657, + "decDeg": 53.529 + }, + "orbitSource": "JPL SSD satellite mean elements, epoch 1950 Jan 1", "parentBodyId": "mars", - "rotationPeriodHours": 30.304279685850094 + "rotationPeriodHours": 30.29857998656265 }, { "id": "io", @@ -223,16 +360,26 @@ "kind": "moon", "radiusKm": 1821.49, "orbit": { - "semiMajorAxisAu": 0.0028210852715482176, - "eccentricity": 0.004228419931613188, - "inclinationDeg": 2.184262830697117, - "longitudeOfAscendingNodeDeg": 338.066968919613, - "argumentOfPeriapsisDeg": 164.6533514126175, - "meanAnomalyAtEpochDeg": 74.88524962049125, - "epochJd": 2460676.5 + "semiMajorAxisAu": 0.0028195588481728304, + "eccentricity": 0.0041, + "inclinationDeg": 0.036, + "longitudeOfAscendingNodeDeg": 43.977, + "argumentOfPeriapsisDeg": 84.129, + "meanAnomalyAtEpochDeg": 342.021, + "epochJd": 2450464.5 }, + "rates": { + "meanMotionDegPerDay": 203.4889583, + "longitudeOfAscendingNodeDegPerDay": -0.1328337309120696, + "argumentOfPeriapsisDegPerDay": -0.6065392513031117 + }, + "laplacePole": { + "raDeg": 268.057, + "decDeg": 64.495 + }, + "orbitSource": "JPL SSD satellite mean elements, epoch 1997 Jan 16", "parentBodyId": "jupiter", - "rotationPeriodHours": 42.51537583211252 + "rotationPeriodHours": 42.45930625514436 }, { "id": "europa", @@ -241,16 +388,26 @@ "kind": "moon", "radiusKm": 1560.8, "orbit": { - "semiMajorAxisAu": 0.004486887602977481, - "eccentricity": 0.009740025569562362, - "inclinationDeg": 2.245018421035425, - "longitudeOfAscendingNodeDeg": 326.0277596344292, - "argumentOfPeriapsisDeg": 349.2045643684611, - "meanAnomalyAtEpochDeg": 40.72575250295771, - "epochJd": 2460676.5 + "semiMajorAxisAu": 0.004486026417754354, + "eccentricity": 0.0094, + "inclinationDeg": 0.466, + "longitudeOfAscendingNodeDeg": 219.106, + "argumentOfPeriapsisDeg": 88.97, + "meanAnomalyAtEpochDeg": 171.016, + "epochJd": 2450464.5 }, + "rates": { + "meanMotionDegPerDay": 101.3747242, + "longitudeOfAscendingNodeDegPerDay": -0.03265393199600969, + "argumentOfPeriapsisDegPerDay": -0.7070489837643877 + }, + "laplacePole": { + "raDeg": 268.084, + "decDeg": 64.506 + }, + "orbitSource": "JPL SSD satellite mean elements, epoch 1997 Jan 16", "parentBodyId": "jupiter", - "rotationPeriodHours": 85.27848729079142 + "rotationPeriodHours": 85.22834531173994 }, { "id": "ganymede", @@ -259,16 +416,26 @@ "kind": "moon", "radiusKm": 2631.2, "orbit": { - "semiMajorAxisAu": 0.007156574251479677, - "eccentricity": 0.001723569500229484, - "inclinationDeg": 2.332721144446764, - "longitudeOfAscendingNodeDeg": 339.4869144536639, - "argumentOfPeriapsisDeg": 0.02754303249978196, - "meanAnomalyAtEpochDeg": 355.7249344187845, - "epochJd": 2460676.5 + "semiMajorAxisAu": 0.007155182055676145, + "eccentricity": 0.0013, + "inclinationDeg": 0.177, + "longitudeOfAscendingNodeDeg": 63.552, + "argumentOfPeriapsisDeg": 192.417, + "meanAnomalyAtEpochDeg": 317.54, + "epochJd": 2450464.5 }, + "rates": { + "meanMotionDegPerDay": 50.3176072, + "longitudeOfAscendingNodeDegPerDay": -0.007430053246547835, + "argumentOfPeriapsisDegPerDay": 0.015509705634511267 + }, + "laplacePole": { + "raDeg": 268.168, + "decDeg": 64.543 + }, + "orbitSource": "JPL SSD satellite mean elements, epoch 1997 Jan 16", "parentBodyId": "jupiter", - "rotationPeriodHours": 171.78271980469088 + "rotationPeriodHours": 171.70927794038664 }, { "id": "callisto", @@ -277,16 +444,26 @@ "kind": "moon", "radiusKm": 2410.3, "orbit": { - "semiMajorAxisAu": 0.012583604816889162, - "eccentricity": 0.007245216892383587, - "inclinationDeg": 1.94976938450231, - "longitudeOfAscendingNodeDeg": 336.7558983297979, - "argumentOfPeriapsisDeg": 32.65151937817475, - "meanAnomalyAtEpochDeg": 126.0198540648681, - "epochJd": 2460676.5 + "semiMajorAxisAu": 0.012585072175094802, + "eccentricity": 0.0074, + "inclinationDeg": 0.192, + "longitudeOfAscendingNodeDeg": 298.848, + "argumentOfPeriapsisDeg": 52.643, + "meanAnomalyAtEpochDeg": 181.408, + "epochJd": 2450464.5 }, + "rates": { + "meanMotionDegPerDay": 21.5710728, + "longitudeOfAscendingNodeDegPerDay": -0.002908996763377476, + "argumentOfPeriapsisDegPerDay": 0.004790407209562851 + }, + "laplacePole": { + "raDeg": 268.639, + "decDeg": 64.749 + }, + "orbitSource": "JPL SSD satellite mean elements, epoch 1997 Jan 16", "parentBodyId": "jupiter", - "rotationPeriodHours": 400.52470451330635 + "rotationPeriodHours": 400.5364072574082 }, { "id": "titan", @@ -295,16 +472,26 @@ "kind": "moon", "radiusKm": 2575.5, "orbit": { - "semiMajorAxisAu": 0.00816836581932601, - "eccentricity": 0.0288317905724529, - "inclinationDeg": 27.71117736818323, - "longitudeOfAscendingNodeDeg": 169.0716057233001, - "argumentOfPeriapsisDeg": 177.4798658320667, - "meanAnomalyAtEpochDeg": 32.18862839469676, - "epochJd": 2460676.5 + "semiMajorAxisAu": 0.008167663044150534, + "eccentricity": 0.0288, + "inclinationDeg": 0.306, + "longitudeOfAscendingNodeDeg": 28.06, + "argumentOfPeriapsisDeg": 180.532, + "meanAnomalyAtEpochDeg": 163.31, + "epochJd": 2451545 }, + "rates": { + "meanMotionDegPerDay": 22.5769756, + "longitudeOfAscendingNodeDegPerDay": -0.001398845136769169, + "argumentOfPeriapsisDegPerDay": 0.002799120423059061 + }, + "laplacePole": { + "raDeg": 36.214, + "decDeg": 83.949 + }, + "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "saturn", - "rotationPeriodHours": 382.87527117206236 + "rotationPeriodHours": 382.69076217631203 }, { "id": "triton", @@ -313,15 +500,25 @@ "kind": "moon", "radiusKm": 1352.6, "orbit": { - "semiMajorAxisAu": 0.002371477351073512, - "eccentricity": 0.0001253161638938937, - "inclinationDeg": 129.1766133444893, - "longitudeOfAscendingNodeDeg": 222.3992799388247, - "argumentOfPeriapsisDeg": 97.422342803062, - "meanAnomalyAtEpochDeg": 273.1157438944629, - "epochJd": 2460676.5 + "semiMajorAxisAu": 0.002371417442908832, + "eccentricity": 0, + "inclinationDeg": 156.865, + "longitudeOfAscendingNodeDeg": 177.608, + "argumentOfPeriapsisDeg": 66.142, + "meanAnomalyAtEpochDeg": 352.257, + "epochJd": 2451545 }, + "rates": { + "meanMotionDegPerDay": 61.2572638, + "longitudeOfAscendingNodeDegPerDay": 0.001433750844964632, + "argumentOfPeriapsisDegPerDay": 0.0025509841146658433 + }, + "laplacePole": { + "raDeg": 299.456, + "decDeg": 43.414 + }, + "orbitSource": "JPL SSD satellite mean elements, epoch 2000 Jan 1", "parentBodyId": "neptune", - "rotationPeriodHours": 141.05626027628523 + "rotationPeriodHours": 141.0444976486201 } ] \ No newline at end of file diff --git a/tools/etl/build.ts b/tools/etl/build.ts index 5ca1709..791e384 100644 --- a/tools/etl/build.ts +++ b/tools/etl/build.ts @@ -1,6 +1,8 @@ import { statSync } from 'node:fs'; -import { BodyRecord } from '../../src/app/shared/models/body.model'; +import { BodyRecord, OrbitalElements } from '../../src/app/shared/models/body.model'; +import { eclipticToEquatorial, laplacePlaneToEquatorial } from '../../src/app/shared/astro/coordinates'; +import { meanElementsAt, positionAtEpoch } from '../../src/app/shared/astro/kepler'; 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'; @@ -130,23 +132,66 @@ function validateMerge(stars: StarRecord[]): void { console.log(` ${survivors} HYG stars have no Gaia counterpart; ${twins} unmerged cross-catalogue pairs within an arcsecond.`); } -function validateBodies(bodies: BodyRecord[]): void { +/** + * How far a body's mean elements may put it from where Horizons has it, on the one date the ETL + * asks Horizons about (2025-01-01), seen from the Sun for a planet and from its planet for a moon. + * + * Measured on this catalogue: the planets at most 0.10 degrees (Uranus; Standish's own stated + * error for his fit is 2 000 arcseconds, 0.56 degrees), the moons at most 1.41 (the Moon, whose + * evection and variation, 1.27 and 0.66 degrees, no mean ellipse has). What this catches is a + * table read wrongly: a moon read against the ecliptic instead of its Laplace plane, a precession + * run the wrong way, or a column taken for its neighbour, which put Triton 26 degrees out and Io + * 0.9. + */ +const MAX_PLANET_OFFSET_DEG = 0.25; +const MAX_MOON_OFFSET_DEG = 2.5; + +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; +} + +function validateBodies(bodies: BodyRecord[], horizonsOrbits: Map): void { assertCondition(bodies.length > 0, 'No solar-system bodies were produced.'); const ids = new Set(bodies.map((body) => body.id)); assertCondition(ids.size === bodies.length, 'Duplicate body ids were found.'); + const offsets: string[] = []; for (const body of bodies) { const orbitValues = Object.values(body.orbit); assertCondition(orbitValues.every(Number.isFinite), `Body ${body.id} has non-finite orbital elements.`); + assertCondition(body.rates.meanMotionDegPerDay > 0, `Body ${body.id} has no mean motion.`); + + // Horizons' elements are osculating, exact at their own epoch; both sets are placed there. + const horizons = horizonsOrbits.get(body.id); + assertCondition(horizons !== undefined, `Body ${body.id} has no Horizons elements to be checked against.`); + const truth = eclipticToEquatorial(positionAtEpoch(horizons!)); + const mean = positionAtEpoch(meanElementsAt(body.orbit, body.rates, horizons!.epochJd)); + const offset = angleBetweenDeg(body.laplacePole ? laplacePlaneToEquatorial(mean, body.laplacePole) : eclipticToEquatorial(mean), truth); + const ceiling = body.kind === 'moon' ? MAX_MOON_OFFSET_DEG : MAX_PLANET_OFFSET_DEG; + assertCondition( + offset <= ceiling, + `${body.name}'s mean elements put it ${offset.toFixed(2)} degrees from where Horizons has it (at most ${ceiling} expected) — the elements were read wrongly.` + ); + offsets.push(`${body.id} ${offset.toFixed(3)}`); if (body.kind === 'moon') { assertCondition(!!body.parentBodyId && ids.has(body.parentBodyId), `Moon ${body.id} has no valid parentBodyId.`); + // Every moon here is tidally locked: its day is its orbit, from the same mean motion that + // carries it round, or its face turns away from its planet: the Kepler period of the + // osculating orbit this used to take would turn the Moon's five degrees an orbit. + const orbitHours = (360 / body.rates.meanMotionDegPerDay) * 24; + 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.` + ); } } const planetCount = bodies.filter((body) => body.kind === 'planet').length; assertCondition(planetCount === 8, `Expected 8 planets, found ${planetCount}.`); + console.log(` mean elements against Horizons, degrees: ${offsets.join(', ')}.`); } function validateExoplanets(exoplanets: ExoplanetRecord[], starIds: Set): void { @@ -234,7 +279,7 @@ async function build(): Promise { const stars = await fetchStars(); console.log(); - const bodies = await fetchSolarSystem(); + const { bodies, horizonsOrbits } = await fetchSolarSystem(); console.log(); const exoplanets = await fetchExoplanets(stars); console.log(); @@ -244,7 +289,7 @@ async function build(): Promise { console.log('Validating output...'); validateStars(stars); validateMerge(stars); - validateBodies(bodies); + validateBodies(bodies, horizonsOrbits); validateExoplanets(exoplanets, new Set(stars.map((star) => star.id))); validateDeepSky(deepSky); diff --git a/tools/etl/fetchSolarSystem.ts b/tools/etl/fetchSolarSystem.ts index 59aabab..516e839 100644 --- a/tools/etl/fetchSolarSystem.ts +++ b/tools/etl/fetchSolarSystem.ts @@ -1,10 +1,10 @@ import { writeFileSync } from 'node:fs'; -import { BodyRecord } from '../../src/app/shared/models/body.model'; -import { gmForParent } from '../../src/app/shared/astro/constants'; -import { orbitalPeriodDays } from '../../src/app/shared/astro/kepler'; +import { BodyRecord, OrbitalElements } from '../../src/app/shared/models/body.model'; import { SUN_STAR_ID } from '../../src/app/shared/models/star.model'; import { fetchHorizonsBody } from './lib/horizons'; +import { parsePlanetMeanElements, parseSatelliteMeanElements } from '../../src/app/shared/astro/mean-elements'; +import { fetchPlanetMeanElementsText, fetchSatelliteMeanElementsHtml } from './lib/mean-elements'; import { dataPath, ensureDataDir } from './lib/paths'; const HOURS_PER_DAY = 24; @@ -21,6 +21,8 @@ interface BodySpec { * WGCCRE 2015 pole (RA 132.99, Dec -6.16), 119.6 degrees: past 90, so it turns retrograde. */ obliquityDeg?: number; + /** The periapsis turns backwards; see `parseSatelliteMeanElements`. */ + apsidesRegress?: boolean; } // Sun-centered planets/dwarf, then their major moons (planetocentric elements). @@ -37,8 +39,8 @@ const BODY_SPECS: BodySpec[] = [ { 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: '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' }, - { id: 'europa', name: 'Europa', kind: 'moon', horizonsCommand: '502', center: '500@599', parentBodyId: 'jupiter' }, + { 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: 'titan', name: 'Titan', kind: 'moon', horizonsCommand: '606', center: '500@699', parentBodyId: 'saturn' }, @@ -46,13 +48,17 @@ const BODY_SPECS: BodySpec[] = [ ]; /** - * Queries JPL Horizons for the osculating orbital elements (and mean radius, where - * reported) of the major planets, Pluto, and a curated set of major moons, and writes - * `bodies.json`. + * Writes `bodies.json` for the major planets, Pluto, and a curated set of major moons: JPL's + * mean orbital elements for where they go, 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. */ -export async function fetchSolarSystem(): Promise { - console.log(`Fetching ${BODY_SPECS.length} solar-system bodies from JPL Horizons...`); +export async function fetchSolarSystem(): Promise<{ bodies: BodyRecord[]; horizonsOrbits: Map }> { + console.log(`Fetching ${BODY_SPECS.length} solar-system bodies from JPL (mean elements, Horizons)...`); const bodies: BodyRecord[] = []; + const horizonsOrbits = new Map(); + const planetElements = await fetchPlanetMeanElementsText(); + const satelliteElements = await fetchSatelliteMeanElementsHtml(); for (const spec of BODY_SPECS) { const result = await fetchHorizonsBody({ @@ -61,17 +67,23 @@ export async function fetchSolarSystem(): Promise { cacheKey: `horizons-${spec.id}.txt` }); + horizonsOrbits.set(spec.id, result.orbit); if (result.radiusKm === undefined) { console.warn(` no physical radius found for ${spec.name}; defaulting to 0.`); } - // Every moon listed here is tidally locked, so its day is its orbit — as drawn, from these - // elements and the parent's mass by Kepler. Not every page says so: the Moon's gives a rate, - // the true sidereal month, 1.4% off the orbit these elements trace, so its face drifted five - // degrees an orbit; Titan's gives nothing, so it did not turn. Taking the orbit keeps one face - // towards the parent, which is what synchronous means. + const parentName = BODY_SPECS.find((candidate) => candidate.id === spec.parentBodyId)?.name; + const mean = parentName + ? parseSatelliteMeanElements(satelliteElements, parentName, spec.name, spec.apsidesRegress ?? false) + : parsePlanetMeanElements(planetElements, spec.id); + + // Every moon listed here is tidally locked, 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 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. const rotationPeriodHours = result.tidallyLocked || spec.kind === 'moon' - ? orbitalPeriodDays(result.orbit.semiMajorAxisAu, gmForParent(spec.parentBodyId)) * HOURS_PER_DAY + ? (360 / mean.rates.meanMotionDegPerDay) * HOURS_PER_DAY : result.rotationPeriodHours; if (rotationPeriodHours === undefined) { console.warn(` no rotation period found for ${spec.name}; it will not turn.`); @@ -83,7 +95,10 @@ export async function fetchSolarSystem(): Promise { name: spec.name, kind: spec.kind, radiusKm: result.radiusKm ?? 0, - orbit: result.orbit, + orbit: mean.orbit, + rates: mean.rates, + ...(mean.laplacePole ? { laplacePole: mean.laplacePole } : {}), + orbitSource: mean.orbitSource, ...(spec.parentBodyId ? { parentBodyId: spec.parentBodyId } : {}), ...(rotationPeriodHours !== undefined ? { rotationPeriodHours } : {}), ...((result.obliquityDeg ?? spec.obliquityDeg) !== undefined ? { obliquityDeg: result.obliquityDeg ?? spec.obliquityDeg } : {}) @@ -93,7 +108,7 @@ export async function fetchSolarSystem(): Promise { ensureDataDir(); writeFileSync(dataPath('bodies.json'), JSON.stringify(bodies, null, 2)); console.log(` wrote ${bodies.length} bodies.`); - return bodies; + return { bodies, horizonsOrbits }; } if (require.main === module) { diff --git a/tools/etl/lib/horizons.ts b/tools/etl/lib/horizons.ts index 067b2a5..aad7f45 100644 --- a/tools/etl/lib/horizons.ts +++ b/tools/etl/lib/horizons.ts @@ -64,8 +64,8 @@ const SECONDS_PER_HOUR = 3600; /** * True where the page gives no number because the body keeps one face to its parent, so its day - * is its orbit. The period itself is then Kepler's, which the caller - * works out from the elements above and the parent's mass. + * is its orbit. The period itself is then the orbit's, which the caller takes from the body's + * mean motion. */ export function isTidallyLocked(text: string): boolean { return SYNCHRONOUS_PATTERN.test(text); diff --git a/tools/etl/lib/mean-elements.ts b/tools/etl/lib/mean-elements.ts new file mode 100644 index 0000000..25e66be --- /dev/null +++ b/tools/etl/lib/mean-elements.ts @@ -0,0 +1,31 @@ +import { fetchTextCached } from './http'; + +/** + * Standish's "Keplerian Elements for Approximate Positions of the Major Planets", Table 2a/2b: + * elements against the J2000 ecliptic and their rates per century, fit to the JPL ephemeris for + * 3000 BC to AD 3000. Table 1 is closer near the present — Saturn within 0.23 degrees of Horizons + * from 1950 to 2100 against this table's 0.32 — but it is only fit for 1800-2050, which the clock + * leaves in minutes, and by AD 3000 it has Saturn 4.3 degrees out where this one is within 0.3 + * of every planet. The page at ssd.jpl.nasa.gov/planets/approx_pos.html carries the same numbers but has + * dropped Pluto, so this reads the plain-text file as JPL last published it, from the Internet + * Archive's copy — pinned to one capture, so the numbers cannot move under the cache. + */ +const PLANET_ELEMENTS_URL = 'https://web.archive.org/web/20210420020242id_/https://ssd.jpl.nasa.gov/txt/p_elem_t2.txt'; + +/** + * JPL SSD's planetary satellite mean elements, as the page stood until 2021: each moon's elements, + * its sidereal mean motion to ten figures, and how fast its node and periapsis turn, against its + * local Laplace plane (the Moon against the ecliptic). The current page, ssd.jpl.nasa.gov/sats/elem, + * has dropped the mean motion and rounds the period to four or five figures — 0.3187 days for + * Phobos, which is a revolution out within a decade — so a period from it would not hold. Pinned + * to one Internet Archive capture for the same reason as the planets. + */ +const SATELLITE_ELEMENTS_URL = 'https://web.archive.org/web/20210203000649id_/https://ssd.jpl.nasa.gov/?sat_elem'; + +export async function fetchPlanetMeanElementsText(): Promise { + return fetchTextCached(PLANET_ELEMENTS_URL, 'jpl-planet-mean-elements-t2.txt'); +} + +export async function fetchSatelliteMeanElementsHtml(): Promise { + return fetchTextCached(SATELLITE_ELEMENTS_URL, 'jpl-satellite-mean-elements.html'); +}