Take the orbits at TDB as the spins already were, so a locked moon faces the planet it is drawn round

Every element set here runs on TDB: Standish's T_eph, the SSD satellite and SBDB epochs, the
IAU's d and T. bodyOrientation already took the clock's UTC to TDB, but SystemOrbitsRenderer.update
and the body page's heliocentricPosition fed the UTC date straight to meanElementsAt, so in one
frame each body's place was 69.184 s behind its spin. That is n x 69 s of orbit: Phobos 0.90
degrees, Mimas 0.31, Deimos 0.23, Enceladus 0.21, Miranda 0.20, Io 0.16, Tethys 0.15, Europa
0.08, the Moon 0.011. 48319c3's table measured the app at a UTC date against Horizons at the same
number read as TDB, which hid it, and its "nothing for anything else" was wrong: Io's 0.16 is four
to five times Io's worst model error there (0.035).

tdbFromUtc, in constants.ts, is now the one conversion, and positions and spins both go through
it. In the running app, clock pinned to 2025-06-01 12:00 UTC, Io's face towards Jupiter is at
0.024 E, latitude -0.009, where Horizons (observer quantity 14 from Jupiter's centre) has 0.036 E
and -0.003: 0.012 degrees apart, where it was 0.175. The renderer spec checks that point, and now
hands its frozen Horizons vectors, which are TDB, to update() as the UTC dates that name them,
69.184 s earlier; the same frozen rows fed at the UTC date fail for Io and Europa. A body-page test
checks that the Sun lights the point it stands over at the same TDB instant the body is turned for.

Controls: taking the renderer's orbits at the clock's UTC fails "faces jupiter and the Sun with the
points Horizons gives on io"; taking the page's Sun there fails "takes the Sun where it stands at
the same TDB instant".

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
This commit is contained in:
2026-09-29 18:45:26 +02:00
co-authored by Claude Opus 5.5
parent f8ee3ac3ab
commit d4808788ec
5 changed files with 56 additions and 20 deletions
@@ -1,7 +1,7 @@
import * as THREE from 'three/webgpu'; import * as THREE from 'three/webgpu';
import { describe, expect, it } from 'vitest'; import { describe, expect, it } from 'vitest';
import { DEFAULT_EPOCH_JD, GM_SUN_AU3_PER_DAY2 } from '../../shared/astro/constants'; import { DEFAULT_EPOCH_JD, GM_SUN_AU3_PER_DAY2, TT_MINUS_UTC_DAYS } from '../../shared/astro/constants';
import { keplerRates } from '../../shared/astro/kepler'; import { keplerRates } from '../../shared/astro/kepler';
import { eclipticToEquatorial, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; import { eclipticToEquatorial, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates';
import { BodyRecord, RotationalElements } from '../../shared/models/body.model'; import { BodyRecord, RotationalElements } from '../../shared/models/body.model';
@@ -203,7 +203,8 @@ describe('SystemOrbitsRenderer exoplanet propagation', () => {
// A body at ecliptic longitude 0 sits on the +X axis in both frames, so it must not move. // A body at ecliptic longitude 0 sits on the +X axis in both frames, so it must not move.
const atEquinox: BodyRecord = { ...EARTH, orbit: { ...EARTH.orbit, eccentricity: 0 } }; const atEquinox: BodyRecord = { ...EARTH, orbit: { ...EARTH.orbit, eccentricity: 0 } };
const renderer = new SystemOrbitsRenderer([atEquinox], []); const renderer = new SystemOrbitsRenderer([atEquinox], []);
renderer.update(DEFAULT_EPOCH_JD); // The clock's UTC date whose TDB is the elements' epoch.
renderer.update(DEFAULT_EPOCH_JD - TT_MINUS_UTC_DAYS);
const p = renderer.members[0].marker.position; const p = renderer.members[0].marker.position;
expect(p.x).toBeCloseTo(1, 6); expect(p.x).toBeCloseTo(1, 6);
@@ -489,7 +490,8 @@ describe('solar-system bodies against Horizons', () => {
// Real records from bodies.json, and Horizons' own positions for them (ICRF, AU; heliocentric // 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 — // 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 // mean elements, their rates, the Laplace planes and the scene's frame — is checked against
// JPL's ephemeris rather than against itself. // JPL's ephemeris rather than against itself. Horizons' dates are TDB and the renderer's are the
// clock's UTC, so each is handed over 69.184 s earlier.
const RECORDS: Record<string, Pick<BodyRecord, 'kind' | 'orbit' | 'rates' | 'laplacePole' | 'parentBodyId' | 'massRatio'>> = { const RECORDS: Record<string, Pick<BodyRecord, 'kind' | 'orbit' | 'rates' | 'laplacePole' | 'parentBodyId' | 'massRatio'>> = {
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}}, 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}}}, 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}}},
@@ -517,6 +519,7 @@ describe('solar-system bodies against Horizons', () => {
earth: {poleRaDeg: [0, -0.641, 0], poleDecDeg: [90, -0.557, 0], primeMeridianDeg: [190.147, 360.9856235, 0]}, earth: {poleRaDeg: [0, -0.641, 0], poleDecDeg: [90, -0.557, 0], primeMeridianDeg: [190.147, 360.9856235, 0]},
mars: {poleRaDeg: [317.269202, -0.10927547, 0], poleDecDeg: [54.432516, -0.05827105, 0], primeMeridianDeg: [176.049863, 350.891982443297, 0], terms: [{angleDeg: [79.398797, 0.5042615, 0], ra: 0.419057, dec: 0, pm: 0}, {angleDeg: [166.325722, 0.5042615, 0], ra: 0, dec: 1.591274, pm: 0}, {angleDeg: [95.391654, 0.5042615, 0], ra: 0, dec: 0, pm: 0.584542}]}, mars: {poleRaDeg: [317.269202, -0.10927547, 0], poleDecDeg: [54.432516, -0.05827105, 0], primeMeridianDeg: [176.049863, 350.891982443297, 0], terms: [{angleDeg: [79.398797, 0.5042615, 0], ra: 0.419057, dec: 0, pm: 0}, {angleDeg: [166.325722, 0.5042615, 0], ra: 0, dec: 1.591274, pm: 0}, {angleDeg: [95.391654, 0.5042615, 0], ra: 0, dec: 0, pm: 0.584542}]},
jupiter: {poleRaDeg: [268.056595, -0.006499, 0], poleDecDeg: [64.495303, 0.002413, 0], primeMeridianDeg: [284.95, 870.536, 0]}, jupiter: {poleRaDeg: [268.056595, -0.006499, 0], poleDecDeg: [64.495303, 0.002413, 0], primeMeridianDeg: [284.95, 870.536, 0]},
io: {poleRaDeg: [268.05, -0.009, 0], poleDecDeg: [64.5, 0.003, 0], primeMeridianDeg: [200.39, 203.4889538, 0], terms: [{angleDeg: [283.9, 4850.7], ra: 0.094, dec: 0.04, pm: -0.085}, {angleDeg: [355.8, 1191.3], ra: 0.024, dec: 0.011, pm: -0.022}]},
saturn: {poleRaDeg: [40.589, -0.036, 0], poleDecDeg: [83.537, -0.004, 0], primeMeridianDeg: [38.9, 810.7939024, 0]}, saturn: {poleRaDeg: [40.589, -0.036, 0], poleDecDeg: [83.537, -0.004, 0], primeMeridianDeg: [38.9, 810.7939024, 0]},
uranus: {poleRaDeg: [257.311, 0, 0], poleDecDeg: [-15.175, 0, 0], primeMeridianDeg: [203.81, -501.1600928, 0]}, uranus: {poleRaDeg: [257.311, 0, 0], poleDecDeg: [-15.175, 0, 0], primeMeridianDeg: [203.81, -501.1600928, 0]},
pluto: {poleRaDeg: [132.993, 0, 0], poleDecDeg: [-6.163, 0, 0], primeMeridianDeg: [302.695, 56.3625225, 0]}, pluto: {poleRaDeg: [132.993, 0, 0], poleDecDeg: [-6.163, 0, 0], primeMeridianDeg: [302.695, 56.3625225, 0]},
@@ -552,7 +555,7 @@ describe('solar-system bodies against Horizons', () => {
for (const [id, jd, x, y, z, maxDeg] of HORIZONS) { for (const [id, jd, x, y, z, maxDeg] of HORIZONS) {
it(`puts ${id} within ${maxDeg} degrees of Horizons on JD ${jd}`, () => { it(`puts ${id} within ${maxDeg} degrees of Horizons on JD ${jd}`, () => {
renderer.update(jd); renderer.update(jd - TT_MINUS_UTC_DAYS);
const drawn = renderer.members.find((member) => member.id === id)!.marker.position; const drawn = renderer.members.find((member) => member.id === id)!.marker.position;
const angleDeg = (drawn.angleTo(new THREE.Vector3(x, y, z)) * 180) / Math.PI; const angleDeg = (drawn.angleTo(new THREE.Vector3(x, y, z)) * 180) / Math.PI;
expect(angleDeg).toBeLessThan(maxDeg); expect(angleDeg).toBeLessThan(maxDeg);
@@ -560,9 +563,9 @@ describe('solar-system bodies against Horizons', () => {
} }
it('puts Pluto where Horizons has it round its barycentre with Charon, 2 131 km out and opposite Charon', () => { it('puts Pluto where Horizons has it round its barycentre with Charon, 2 131 km out and opposite Charon', () => {
// Horizons, Pluto (999) from the Pluto-system barycentre (9), on JD 2488069.5 (2100). // Horizons, Pluto (999) from the Pluto-system barycentre (9), on JD 2488069.5 TDB (2100).
const horizons = new THREE.Vector3(0.000003313612032581019, 0.000001023040948538272, -0.00001381793390079716); const horizons = new THREE.Vector3(0.000003313612032581019, 0.000001023040948538272, -0.00001381793390079716);
renderer.update(2488069.5); renderer.update(2488069.5 - TT_MINUS_UTC_DAYS);
const charon = renderer.members.find((member) => member.id === 'charon')!.marker; const charon = renderer.members.find((member) => member.id === 'charon')!.marker;
const barycentre = charon.parent!.position; const barycentre = charon.parent!.position;
const pluto = renderer.members.find((member) => member.id === 'pluto')!.marker.position.clone().sub(barycentre); const pluto = renderer.members.find((member) => member.id === 'pluto')!.marker.position.clone().sub(barycentre);
@@ -712,12 +715,14 @@ describe('solar-system bodies against Horizons', () => {
// //
// Measured: every longitude within 0.09 degrees and every latitude within 0.03, but for the // Measured: every longitude within 0.09 degrees and every latitude within 0.03, but for the
// Moon's face towards Earth, 0.70 and 0.09 out because its mean orbit is (its evection alone is // Moon's face towards Earth, 0.70 and 0.09 out because its mean orbit is (its evection alone is
// 1.27 degrees); its face towards the Sun is within 0.002. // 1.27 degrees); its face towards the Sun is within 0.002. Io's face towards Jupiter is 0.012 out:
// with its orbit taken at the clock's UTC and its spin at TDB it was 0.175, the 69 s between them.
const SUB_POINTS: Array<[id: string, observer: string | undefined, lightMinutes: number, west: boolean, flattening: number, observerLon: number, observerLat: number, sunLon: number, sunLat: number, maxObserverDeg: number]> = [ const SUB_POINTS: Array<[id: string, observer: string | undefined, lightMinutes: number, west: boolean, flattening: number, observerLon: number, observerLat: number, sunLon: number, sunLat: number, maxObserverDeg: number]> = [
['earth', undefined, 8.43351424, false, 1 / 298.257, 1.5855, 22.261204, 1.579501, 22.260426, 0.1], ['earth', undefined, 8.43351424, false, 1 / 298.257, 1.5855, 22.261204, 1.579501, 22.260426, 0.1],
['mars', 'earth', 14.13295841, true, 1 - 3376.2 / 3396.19, 307.365389, 21.27653, 269.287887, 25.451264, 0.1], ['mars', 'earth', 14.13295841, true, 1 - 3376.2 / 3396.19, 307.365389, 21.27653, 269.287887, 25.451264, 0.1],
['moon', 'earth', 0.02150549, false, 0, 7.256763, -3.462104, 116.285934, 1.503004, 0.8], ['moon', 'earth', 0.02150549, false, 0, 7.256763, -3.462104, 116.285934, 1.503004, 0.8],
['jupiter', 'earth', 50.70337676, true, 1 - 66854 / 71492, 251.139846, 2.58787, 247.855871, 2.572658, 0.1] ['jupiter', 'earth', 50.70337676, true, 1 - 66854 / 71492, 251.139846, 2.58787, 247.855871, 2.572658, 0.1],
['io', 'jupiter', 0.02340584, true, 0, 359.964094, -0.002537, 355.673108, 2.26528, 0.05]
]; ];
for (const [id, observer, lightMinutes, west, flattening, observerLon, observerLat, sunLon, sunLat, maxObserverDeg] of SUB_POINTS) { for (const [id, observer, lightMinutes, west, flattening, observerLon, observerLat, sunLon, sunLat, maxObserverDeg] of SUB_POINTS) {
@@ -6,6 +6,7 @@ import { planetTexture } from '../../shared/rendering/procedural-planet-texture'
import { bodyTexturePath, loadCachedTexture, saturnRing } from '../../shared/rendering/texture-catalog'; import { bodyTexturePath, loadCachedTexture, saturnRing } from '../../shared/rendering/texture-catalog';
import { isPropagatableOrbit, keplerRates, meanElementsAt, orbitEllipsePoints, positionAtEpoch, resolveGravitationalParameter, resolveOrbitalElements } from '../../shared/astro/kepler'; import { isPropagatableOrbit, keplerRates, meanElementsAt, orbitEllipsePoints, positionAtEpoch, resolveGravitationalParameter, resolveOrbitalElements } from '../../shared/astro/kepler';
import { CartesianCoordinates, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; import { CartesianCoordinates, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates';
import { tdbFromUtc } from '../../shared/astro/constants';
import { BodyRecord, MeanElementRates, OrbitalElements, RotationalElements } from '../../shared/models/body.model'; import { BodyRecord, MeanElementRates, OrbitalElements, RotationalElements } from '../../shared/models/body.model';
import { bodyOrientation, poleFrame } from '../../shared/rendering/body-orientation'; import { bodyOrientation, poleFrame } from '../../shared/rendering/body-orientation';
import { bodyMarkerRadiusAu, systemGridRingsAu } from './system-framing'; import { bodyMarkerRadiusAu, systemGridRingsAu } from './system-framing';
@@ -466,10 +467,14 @@ export class SystemOrbitsRenderer {
this.object.add(starLight()); this.object.add(starLight());
} }
/** Recomputes every marker's position for the given Julian date. Call once per tick. */ /**
* Recomputes every marker's position for the given Julian date, UTC as the map's clock gives it:
* the orbits are taken at its TDB, as the spins are. Call once per tick.
*/
update(epochJd: number): void { update(epochJd: number): void {
const jdTdb = tdbFromUtc(epochJd);
for (const body of this.topLevelBodies) { for (const body of this.topLevelBodies) {
const current = meanElementsAt(body.elements, body.rates, epochJd); const current = meanElementsAt(body.elements, body.rates, jdTdb);
const orbital = positionAtEpoch(current); const orbital = positionAtEpoch(current);
body.position.set(orbital.x, orbital.y, orbital.z).applyQuaternion(body.frame); body.position.set(orbital.x, orbital.y, orbital.z).applyQuaternion(body.frame);
body.marker.position.copy(body.position); body.marker.position.copy(body.position);
@@ -477,7 +482,7 @@ export class SystemOrbitsRenderer {
if (body.rotationalElements) { if (body.rotationalElements) {
bodyOrientation(body.rotationalElements, epochJd, body.marker.quaternion); bodyOrientation(body.rotationalElements, epochJd, body.marker.quaternion);
} else if (body.rotationPeriodHours) { } else if (body.rotationPeriodHours) {
body.marker.quaternion.copy(spinFor(current, body.frame, body.rotationPeriodHours, epochJd - body.elements.epochJd)); body.marker.quaternion.copy(spinFor(current, body.frame, body.rotationPeriodHours, jdTdb - body.elements.epochJd));
} }
} }
@@ -487,7 +492,7 @@ export class SystemOrbitsRenderer {
continue; continue;
} }
moon.pivot.position.copy(parent.position); moon.pivot.position.copy(parent.position);
const current = meanElementsAt(moon.elements, moon.rates, epochJd); const current = meanElementsAt(moon.elements, moon.rates, jdTdb);
const orbital = positionAtEpoch(current); const orbital = positionAtEpoch(current);
moon.marker.position.set(orbital.x, orbital.y, orbital.z).applyQuaternion(moon.frame); moon.marker.position.set(orbital.x, orbital.y, orbital.z).applyQuaternion(moon.frame);
orientOrbit(moon.orbitLine.quaternion, current, moon.frame); orientOrbit(moon.orbitLine.quaternion, current, moon.frame);
@@ -502,7 +507,7 @@ export class SystemOrbitsRenderer {
if (moon.rotationalElements) { if (moon.rotationalElements) {
bodyOrientation(moon.rotationalElements, epochJd, moon.marker.quaternion); bodyOrientation(moon.rotationalElements, epochJd, moon.marker.quaternion);
} else if (moon.rotationPeriodHours) { } else if (moon.rotationPeriodHours) {
moon.marker.quaternion.copy(spinFor(current, moon.frame, moon.rotationPeriodHours, epochJd - moon.elements.epochJd)); moon.marker.quaternion.copy(spinFor(current, moon.frame, moon.rotationPeriodHours, jdTdb - moon.elements.epochJd));
} }
} }
+11
View File
@@ -27,6 +27,17 @@ export const GM_SUN_AU3_PER_DAY2 = 0.01720209895 * 0.01720209895;
*/ */
export const TT_MINUS_UTC_DAYS = 69.184 / 86400; export const TT_MINUS_UTC_DAYS = 69.184 / 86400;
/**
* The TDB date every element set here is evaluated at, for a date on the map's clock, which is
* UTC: Standish's T_eph, the SSD satellite and SBDB epochs and the IAU's d and T all run on TDB.
* Positions and spins both go through this, so a locked moon's face and the orbit it is drawn on
* are taken at the same instant; taken at the clock's date, the orbits ran 69 s behind the spins,
* which is 0.9 degrees of Phobos's orbit and 0.16 of Io's.
*/
export function tdbFromUtc(jdUtc: number): number {
return jdUtc + TT_MINUS_UTC_DAYS;
}
/** Converts a JS `Date` into a Julian date (days), for driving the Kepler propagator "now". */ /** Converts a JS `Date` into a Julian date (days), for driving the Kepler propagator "now". */
export function dateToJulianDate(date: Date = new Date()): number { export function dateToJulianDate(date: Date = new Date()): number {
return date.getTime() / 86400000 + 2440587.5; return date.getTime() / 86400000 + 2440587.5;
@@ -1,8 +1,11 @@
import * as THREE from 'three/webgpu'; import * as THREE from 'three/webgpu';
import { describe, expect, it } from 'vitest'; import { describe, expect, it } from 'vitest';
import { TT_MINUS_UTC_DAYS } from '../astro/constants';
import { eclipticToEquatorial } from '../astro/coordinates';
import { meanElementsAt, positionAtEpoch } from '../astro/kepler';
import { BodyRecord } from '../models/body.model'; import { BodyRecord } from '../models/body.model';
import { bodyPageView } from './body-orientation'; import { bodyOrientation, bodyPageView } from './body-orientation';
// Earth (the Earth-Moon barycentre's mean elements) and the Moon as bodies.json carries them. // Earth (the Earth-Moon barycentre's mean elements) and the Moon as bodies.json carries them.
const EARTH: BodyRecord = { const EARTH: BodyRecord = {
@@ -51,6 +54,17 @@ describe('bodyPageView', () => {
expect(Math.abs(moon.latDeg - 1.503004)).toBeLessThan(0.05); expect(Math.abs(moon.latDeg - 1.503004)).toBeLessThan(0.05);
}); });
it('takes the Sun where it stands at the same TDB instant the body is turned for', () => {
// Earth's own sphere, turned as the system view turns it, and the Sun seen from Earth's mean
// place at the clock's date taken to TDB: the page must light that same point of its map.
const planet = new THREE.Quaternion();
const sun = new THREE.Vector3();
bodyPageView(EARTH, BODIES, JUNE_1_2025_NOON_UTC, SUN_AZIMUTH, planet, sun);
const place = eclipticToEquatorial(positionAtEpoch(meanElementsAt(EARTH.orbit, EARTH.rates, JUNE_1_2025_NOON_UTC + TT_MINUS_UTC_DAYS)));
const expected = new THREE.Vector3(-place.x, -place.y, -place.z).normalize().applyQuaternion(bodyOrientation(EARTH.rotationalElements!, JUNE_1_2025_NOON_UTC).invert());
expect(sun.clone().applyQuaternion(planet.clone().invert()).angleTo(expected)).toBeLessThan(1e-9);
});
it('keeps the pole up and the Sun where the page’s light stands, turning the body under it', () => { it('keeps the pole up and the Sun where the page’s light stands, turning the body under it', () => {
const planet = new THREE.Quaternion(); const planet = new THREE.Quaternion();
const sun = new THREE.Vector3(); const sun = new THREE.Vector3();
+7 -6
View File
@@ -1,6 +1,6 @@
import * as THREE from 'three/webgpu'; import * as THREE from 'three/webgpu';
import { TT_MINUS_UTC_DAYS } from '../astro/constants'; import { tdbFromUtc } from '../astro/constants';
import { CartesianCoordinates, eclipticToEquatorial, laplacePlaneToEquatorial } from '../astro/coordinates'; import { CartesianCoordinates, eclipticToEquatorial, laplacePlaneToEquatorial } from '../astro/coordinates';
import { meanElementsAt, positionAtEpoch } from '../astro/kepler'; import { meanElementsAt, positionAtEpoch } from '../astro/kepler';
import { orientationAt } from '../astro/rotational-elements'; import { orientationAt } from '../astro/rotational-elements';
@@ -51,10 +51,10 @@ export function poleFrame(pole: { raDeg: number; decDeg: number }, target = new
* and on to where the pole points. * and on to where the pole points.
* *
* The clock is UTC and the IAU's elements run on TDB, 69.184 s ahead; in that time Earth turns * The clock is UTC and the IAU's elements run on TDB, 69.184 s ahead; in that time Earth turns
* 0.29 degrees, Jupiter 0.70 and Phobos 0.90, so the difference is added here. * 0.29 degrees, Jupiter 0.70 and Phobos 0.90, so the date is taken to TDB here (see `tdbFromUtc`).
*/ */
export function bodyOrientation(elements: RotationalElements, jdUtc: number, target = new THREE.Quaternion()): THREE.Quaternion { export function bodyOrientation(elements: RotationalElements, jdUtc: number, target = new THREE.Quaternion()): THREE.Quaternion {
const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(elements, jdUtc + TT_MINUS_UTC_DAYS); const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(elements, tdbFromUtc(jdUtc));
return poleFrame({ raDeg: poleRaDeg, decDeg: poleDecDeg }, target) return poleFrame({ raDeg: poleRaDeg, decDeg: poleDecDeg }, target)
.multiply(scratchTurn.setFromAxisAngle(Z_AXIS, primeMeridianDeg * DEG_TO_RAD)) .multiply(scratchTurn.setFromAxisAngle(Z_AXIS, primeMeridianDeg * DEG_TO_RAD))
.multiply(MAP_TO_BODY); .multiply(MAP_TO_BODY);
@@ -62,13 +62,14 @@ export function bodyOrientation(elements: RotationalElements, jdUtc: number, tar
/** Where a body is from the Sun at a date, in the ICRF, AU: a moon's planet's place plus its own. */ /** Where a body is from the Sun at a date, in the ICRF, AU: a moon's planet's place plus its own. */
function heliocentricPosition(body: BodyRecord, bodies: readonly BodyRecord[], jdUtc: number): CartesianCoordinates { function heliocentricPosition(body: BodyRecord, bodies: readonly BodyRecord[], jdUtc: number): CartesianCoordinates {
const own = positionAtEpoch(meanElementsAt(body.orbit, body.rates, jdUtc)); const jdTdb = tdbFromUtc(jdUtc);
const own = positionAtEpoch(meanElementsAt(body.orbit, body.rates, jdTdb));
const parent = body.parentBodyId ? bodies.find((candidate) => candidate.id === body.parentBodyId) : undefined; const parent = body.parentBodyId ? bodies.find((candidate) => candidate.id === body.parentBodyId) : undefined;
if (!parent) { if (!parent) {
return eclipticToEquatorial(own); return eclipticToEquatorial(own);
} }
const offset = body.laplacePole ? laplacePlaneToEquatorial(own, body.laplacePole) : eclipticToEquatorial(own); const offset = body.laplacePole ? laplacePlaneToEquatorial(own, body.laplacePole) : eclipticToEquatorial(own);
const centre = eclipticToEquatorial(positionAtEpoch(meanElementsAt(parent.orbit, parent.rates, jdUtc))); const centre = eclipticToEquatorial(positionAtEpoch(meanElementsAt(parent.orbit, parent.rates, jdTdb)));
return { x: centre.x + offset.x, y: centre.y + offset.y, z: centre.z + offset.z }; return { x: centre.x + offset.x, y: centre.y + offset.y, z: centre.z + offset.z };
} }
@@ -92,7 +93,7 @@ export function bodyPageView(body: BodyRecord, bodies: readonly BodyRecord[], jd
if (!elements) { if (!elements) {
return false; return false;
} }
const { poleRaDeg, poleDecDeg } = orientationAt(elements, jdUtc + TT_MINUS_UTC_DAYS); const { poleRaDeg, poleDecDeg } = orientationAt(elements, tdbFromUtc(jdUtc));
// From the ICRF into the body's frame with its pole on +Y, before the turn about that pole. // From the ICRF into the body's frame with its pole on +Y, before the turn about that pole.
const toPage = poleFrame({ raDeg: poleRaDeg, decDeg: poleDecDeg }, scratchPage).multiply(MAP_TO_BODY).invert(); const toPage = poleFrame({ raDeg: poleRaDeg, decDeg: poleDecDeg }, scratchPage).multiply(MAP_TO_BODY).invert();
const position = heliocentricPosition(body, bodies, jdUtc); const position = heliocentricPosition(body, bodies, jdUtc);