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 db780fe..298e5a7 100644 --- a/src/app/features/galaxy-system/system-orbits-renderer.spec.ts +++ b/src/app/features/galaxy-system/system-orbits-renderer.spec.ts @@ -1,13 +1,17 @@ import * as THREE from 'three/webgpu'; import { describe, expect, it } from 'vitest'; -import { DEFAULT_EPOCH_JD, GM_SUN_AU3_PER_DAY2, TT_MINUS_UTC_DAYS } from '../../shared/astro/constants'; +import { DEFAULT_EPOCH_JD, GM_SUN_AU3_PER_DAY2, TT_MINUS_UTC_DAYS, ttMinusUtSeconds } from '../../shared/astro/constants'; import { keplerRates } from '../../shared/astro/kepler'; -import { eclipticToEquatorial, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; +import { eclipticToEquatorial, laplacePlaneToEquatorial, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; +import { orientationAt } from '../../shared/astro/rotational-elements'; import { BodyRecord, RotationalElements } from '../../shared/models/body.model'; import { ExoplanetRecord } from '../../shared/models/exoplanet.model'; import { SystemOrbitsRenderer } from './system-orbits-renderer'; +/** The clock's UT date that names a TDB one: TT - UT, which moves by under a second a year, earlier. */ +const utOf = (jdTdb: number): number => jdTdb - ttMinusUtSeconds(jdTdb) / 86400; + /** TRAPPIST-1 b: a real short-period planet around a 0.09 solar-mass red dwarf. */ const TRAPPIST_1B_SEMI_MAJOR_AXIS_AU = 0.01154; const TRAPPIST_1B_PERIOD_DAYS = 1.51088; @@ -203,8 +207,8 @@ describe('SystemOrbitsRenderer exoplanet propagation', () => { // 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 renderer = new SystemOrbitsRenderer([atEquinox], []); - // The clock's UTC date whose TDB is the elements' epoch. - renderer.update(DEFAULT_EPOCH_JD - TT_MINUS_UTC_DAYS); + // The clock's UT date whose TDB is the elements' epoch. + renderer.update(utOf(DEFAULT_EPOCH_JD)); const p = renderer.members[0].marker.position; expect(p.x).toBeCloseTo(1, 6); @@ -491,7 +495,7 @@ describe('solar-system bodies against Horizons', () => { // 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. Horizons' dates are TDB and the renderer's are the - // clock's UTC, so each is handed over 69.184 s earlier. + // clock's UT, so each is handed over TT - UT earlier: 69.184 s today, 29 in 1950. 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}}}, @@ -555,7 +559,7 @@ describe('solar-system bodies against Horizons', () => { for (const [id, jd, x, y, z, maxDeg] of HORIZONS) { it(`puts ${id} within ${maxDeg} degrees of Horizons on JD ${jd}`, () => { - renderer.update(jd - TT_MINUS_UTC_DAYS); + renderer.update(utOf(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); @@ -565,7 +569,7 @@ 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', () => { // 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); - renderer.update(2488069.5 - TT_MINUS_UTC_DAYS); + renderer.update(utOf(2488069.5)); const charon = renderer.members.find((member) => member.id === 'charon')!.marker; const barycentre = charon.parent!.position; const pluto = renderer.members.find((member) => member.id === 'pluto')!.marker.position.clone().sub(barycentre); @@ -656,6 +660,33 @@ describe('solar-system bodies against Horizons', () => { /** Degrees between two longitudes, the short way round. */ const apart = (a: number, b: number): number => Math.abs(((((a - b) % 360) + 540) % 360) - 180); + /** Where the IAU puts a body's prime meridian at a TDB date, in the scene. */ + function iauPrimeMeridian(id: string, jdTdb: number): THREE.Vector3 { + const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(ROTATION[id], jdTdb); + const w = (primeMeridianDeg * Math.PI) / 180; + const meridian = laplacePlaneToEquatorial({ x: Math.cos(w), y: Math.sin(w), z: 0 }, { raDeg: poleRaDeg, decDeg: poleDecDeg }); + return new THREE.Vector3(meridian.x, meridian.y, meridian.z); + } + + /** The drawn sphere's longitude 0 on its equator: +X of the sphere as `SphereGeometry` wraps its map. */ + function drawnPrimeMeridian(id: string): THREE.Vector3 { + return new THREE.Vector3(1, 0, 0).applyQuaternion(renderer.members.find((member) => member.id === id)!.marker.quaternion); + } + + it('turns Jupiter at AD 1000 by its W at that date’s TT, 1 574 s after the UT the clock names', () => { + // Espenak and Meeus's ΔT for JD 2086307.5, 1 January 1000 in the Julian calendar, where TT - UT was 23 times what it is today: held + // at today's 69 s, Jupiter was drawn 15 degrees short of its W. + const jdUt = 2086307.5; + renderer.update(jdUt); + expect((drawnPrimeMeridian('jupiter').angleTo(iauPrimeMeridian('jupiter', jdUt + 1574.1 / 86400)) * 180) / Math.PI).toBeLessThan(0.01); + }); + + it('turns Earth by the UT the clock names, which is its turning: at AD 1000 its W is not moved on by ΔT', () => { + const jdUt = 2086307.5; + renderer.update(jdUt); + expect((drawnPrimeMeridian('earth').angleTo(iauPrimeMeridian('earth', jdUt + TT_MINUS_UTC_DAYS)) * 180) / Math.PI).toBeLessThan(0.01); + }); + it('lights Earth where the Sun really stands: within 4 degrees of Greenwich at noon UTC', () => { // The equation of time is all that separates them: on 1 June 2025 it puts the Sun over 0.53 W, // and the drawn sphere has it over 0.43 W. diff --git a/src/app/features/galaxy-system/system-orbits-renderer.ts b/src/app/features/galaxy-system/system-orbits-renderer.ts index 5950d2d..fdadb0f 100644 --- a/src/app/features/galaxy-system/system-orbits-renderer.ts +++ b/src/app/features/galaxy-system/system-orbits-renderer.ts @@ -480,7 +480,7 @@ export class SystemOrbitsRenderer { body.marker.position.copy(body.position); orientOrbit(body.orbitLine.quaternion, current, body.frame); if (body.rotationalElements) { - bodyOrientation(body.rotationalElements, epochJd, body.marker.quaternion); + bodyOrientation(body.rotationalElements, epochJd, body.marker.quaternion, body.id === 'earth'); } else if (body.rotationPeriodHours) { body.marker.quaternion.copy(spinFor(current, body.frame, body.rotationPeriodHours, jdTdb - body.elements.epochJd)); } diff --git a/src/app/shared/astro/constants.spec.ts b/src/app/shared/astro/constants.spec.ts new file mode 100644 index 0000000..40a597b --- /dev/null +++ b/src/app/shared/astro/constants.spec.ts @@ -0,0 +1,36 @@ +import { describe, expect, it } from 'vitest'; + +import { tdbFromUtc, ttMinusUtSeconds } from './constants'; + +const jd = (year: number, month = 1, day = 1): number => Date.UTC(year, month - 1, day) / 86400000 + 2440587.5; + +describe('ttMinusUtSeconds', () => { + it('follows the historical record before 1972: within 1 per cent of Horizons at AD 1, 6 per cent at AD 1000, 0.2 s in 1950', () => { + // Horizons' TDB - UT (observer quantity 30) on JD 1721600, 2086455 and 2433282.5. + expect(Math.abs(ttMinusUtSeconds(1721600) - 10465.73)).toBeLessThan(105); + expect(Math.abs(ttMinusUtSeconds(2086455) - 1658.0)).toBeLessThan(100); + expect(Math.abs(ttMinusUtSeconds(2433282.5) - 28.93)).toBeLessThan(0.2); + }); + + it('counts the leap seconds from 1972, and holds the last from 2017 on', () => { + expect(ttMinusUtSeconds(jd(1972, 6, 30))).toBe(42.184); + expect(ttMinusUtSeconds(jd(1972, 7, 1))).toBe(43.184); + expect(ttMinusUtSeconds(jd(2016, 12, 31))).toBe(68.184); + expect(ttMinusUtSeconds(jd(2017, 1, 1))).toBe(69.184); + expect(ttMinusUtSeconds(jd(2999, 1, 1))).toBe(69.184); + }); + + it('joins its pieces without a jump of more than a second', () => { + for (const year of [500, 1600, 1700, 1800, 1860, 1900, 1920, 1941, 1961, 1972]) { + const at = 2451544.5 + (year - 2000) * 365.2425; + expect(Math.abs(ttMinusUtSeconds(at + 0.01) - ttMinusUtSeconds(at - 0.01))).toBeLessThan(1); + } + }); +}); + +describe('tdbFromUtc', () => { + it('puts the clock’s date that far on', () => { + expect((tdbFromUtc(jd(2025)) - jd(2025)) * 86400).toBeCloseTo(69.184, 3); + expect((tdbFromUtc(2086455) - 2086455) * 86400).toBeCloseTo(ttMinusUtSeconds(2086455), 3); + }); +}); diff --git a/src/app/shared/astro/constants.ts b/src/app/shared/astro/constants.ts index b54dae3..5146705 100644 --- a/src/app/shared/astro/constants.ts +++ b/src/app/shared/astro/constants.ts @@ -20,22 +20,86 @@ export const DEFAULT_EPOCH_JD = 2451545.0; export const GM_SUN_AU3_PER_DAY2 = 0.01720209895 * 0.01720209895; /** - * TT - UTC, in days: 32.184 s plus the 37 leap seconds UTC has taken since 1972, the last at the - * end of 2016. TDB, which ephemerides run on, stays within 2 ms of TT. Held constant, as Horizons - * holds it for dates past the last announced leap second; before 2017 it was smaller, about 29 s - * in 1950. + * TT - UTC today, in days: 32.184 s plus the 37 leap seconds UTC has taken since 1972, the last at + * the end of 2016. TDB, which ephemerides run on, stays within 2 ms of TT. See {@link ttMinusUtSeconds} + * for other dates. */ export const TT_MINUS_UTC_DAYS = 69.184 / 86400; +/** The first day of each month UTC took a leap second at the start of, from its 10 s of 1972. */ +const LEAP_SECONDS_FROM = [ + [1972, 7], [1973, 1], [1974, 1], [1975, 1], [1976, 1], [1977, 1], [1978, 1], [1979, 1], [1980, 1], [1981, 7], + [1982, 7], [1983, 7], [1985, 7], [1988, 1], [1990, 1], [1991, 1], [1992, 7], [1993, 7], [1994, 7], [1996, 1], + [1997, 7], [1999, 1], [2006, 1], [2009, 1], [2012, 7], [2015, 7], [2017, 1] +].map(([year, month]) => Date.UTC(year, month - 1, 1) / 86400000 + 2440587.5); +const JD_1972 = Date.UTC(1972, 0, 1) / 86400000 + 2440587.5; + +/** + * TT - UT, in seconds, at a date on the map's clock: how far Earth's turning, which UT counts, + * has fallen behind the uniform time the ephemerides run on. + * + * From 1972 the clock is UTC, held to within 0.9 s of UT by leap seconds, and TT - UTC is exact: + * 32.184 s plus the 10 to 37 of them. After the last, at the start of 2017, it is held at 69.184 s, + * as Horizons holds it: no one knows the leap seconds to come. Before 1972 it is ΔT from the + * Espenak-Meeus polynomials (NASA's Five Millennium Canon, 2006), which fit the historical record + * of eclipses and occultations: 10 570 s at AD 1, 1 574 at AD 1000, 29 in 1950. Held at 69 s there, + * as it was, every spin but Earth's was a turn of (ΔT - 69 s) times its rate out, 15 degrees for + * Jupiter at AD 1000 and 106 at AD 1, and the Moon 0.22 and 1.43 degrees along its orbit. + */ +export function ttMinusUtSeconds(jdUt: number): number { + if (jdUt >= JD_1972) { + return 32.184 + 10 + LEAP_SECONDS_FROM.filter((from) => jdUt >= from).length; + } + const y = 2000 + (jdUt - 2451544.5) / 365.2425; + if (y < 500) { + const u = y / 100; + return 10583.6 - 1014.41 * u + 33.78311 * u ** 2 - 5.952053 * u ** 3 - 0.1798452 * u ** 4 + 0.022174192 * u ** 5 + 0.0090316521 * u ** 6; + } + if (y < 1600) { + const u = (y - 1000) / 100; + return 1574.2 - 556.01 * u + 71.23472 * u ** 2 + 0.319781 * u ** 3 - 0.8503463 * u ** 4 - 0.005050998 * u ** 5 + 0.0083572073 * u ** 6; + } + if (y < 1700) { + const t = y - 1600; + return 120 - 0.9808 * t - 0.01532 * t ** 2 + t ** 3 / 7129; + } + if (y < 1800) { + const t = y - 1700; + return 8.83 + 0.1603 * t - 0.0059285 * t ** 2 + 0.00013336 * t ** 3 - t ** 4 / 1174000; + } + if (y < 1860) { + const t = y - 1800; + return 13.72 - 0.332447 * t + 0.0068612 * t ** 2 + 0.0041116 * t ** 3 - 0.00037436 * t ** 4 + 0.0000121272 * t ** 5 - 0.0000001699 * t ** 6 + 0.000000000875 * t ** 7; + } + if (y < 1900) { + const t = y - 1860; + return 7.62 + 0.5737 * t - 0.251754 * t ** 2 + 0.01680668 * t ** 3 - 0.0004473624 * t ** 4 + t ** 5 / 233174; + } + if (y < 1920) { + const t = y - 1900; + return -2.79 + 1.494119 * t - 0.0598939 * t ** 2 + 0.0061966 * t ** 3 - 0.000197 * t ** 4; + } + if (y < 1941) { + const t = y - 1920; + return 21.2 + 0.84493 * t - 0.0761 * t ** 2 + 0.0020936 * t ** 3; + } + if (y < 1961) { + const t = y - 1950; + return 29.07 + 0.407 * t - t ** 2 / 233 + t ** 3 / 2547; + } + const t = y - 1975; + return 45.45 + 1.067 * t - t ** 2 / 260 - t ** 3 / 718; +} + /** * 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. + * UT: 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; + return jdUtc + ttMinusUtSeconds(jdUtc) / 86400; } /** Converts a JS `Date` into a Julian date (days), for driving the Kepler propagator "now". */ diff --git a/src/app/shared/rendering/body-orientation.spec.ts b/src/app/shared/rendering/body-orientation.spec.ts index 82bef8a..5519052 100644 --- a/src/app/shared/rendering/body-orientation.spec.ts +++ b/src/app/shared/rendering/body-orientation.spec.ts @@ -1,7 +1,7 @@ import * as THREE from 'three/webgpu'; import { describe, expect, it } from 'vitest'; -import { TT_MINUS_UTC_DAYS } from '../astro/constants'; +import { tdbFromUtc } from '../astro/constants'; import { eclipticToEquatorial } from '../astro/coordinates'; import { meanElementsAt, positionAtEpoch } from '../astro/kepler'; import { BodyRecord } from '../models/body.model'; @@ -54,15 +54,18 @@ describe('bodyPageView', () => { 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('takes the Sun where it stands at the same TDB instant the body is turned for, and Earth turned as the system view turns it', () => { + // Earth's own sphere, turned as the system view turns it (by UT, see `bodyOrientation`), 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, today and at AD 1000, when TT was 1 574 s past UT. + for (const jdUt of [JUNE_1_2025_NOON_UTC, 2086307.5]) { + const planet = new THREE.Quaternion(); + const sun = new THREE.Vector3(); + bodyPageView(EARTH, BODIES, jdUt, SUN_AZIMUTH, planet, sun); + const place = eclipticToEquatorial(positionAtEpoch(meanElementsAt(EARTH.orbit, EARTH.rates, tdbFromUtc(jdUt)))); + const expected = new THREE.Vector3(-place.x, -place.y, -place.z).normalize().applyQuaternion(bodyOrientation(EARTH.rotationalElements!, jdUt, undefined, true).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', () => { diff --git a/src/app/shared/rendering/body-orientation.ts b/src/app/shared/rendering/body-orientation.ts index 3755e84..ccf270e 100644 --- a/src/app/shared/rendering/body-orientation.ts +++ b/src/app/shared/rendering/body-orientation.ts @@ -1,6 +1,6 @@ import * as THREE from 'three/webgpu'; -import { tdbFromUtc } from '../astro/constants'; +import { tdbFromUtc, TT_MINUS_UTC_DAYS } from '../astro/constants'; import { CartesianCoordinates, eclipticToEquatorial, laplacePlaneToEquatorial } from '../astro/coordinates'; import { meanElementsAt, positionAtEpoch } from '../astro/kepler'; import { orientationAt } from '../astro/rotational-elements'; @@ -50,11 +50,14 @@ export function poleFrame(pole: { raDeg: number; decDeg: number }, target = new * into the ICRF at the map's own clock: the map onto the body's frame, turned by W about the pole, * 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 - * 0.29 degrees, Jupiter 0.70 and Phobos 0.90, so the date is taken to TDB here (see `tdbFromUtc`). + * The clock is UT and the IAU's elements run on TDB, 69.184 s ahead today and 1 574 s at AD 1000; + * in 69 s Earth turns 0.29 degrees, Jupiter 0.70 and Phobos 0.90, so the date is taken to TDB here + * (see `tdbFromUtc`). Earth, `followsUt`, is the one exception: its turning is what UT counts, + * so the clock's date already says how far it has turned, and its W, fitted to today, is taken at + * that date plus today's TT - UTC. Taken at TDB, it would turn ΔT further: 44 degrees at AD 1. */ -export function bodyOrientation(elements: RotationalElements, jdUtc: number, target = new THREE.Quaternion()): THREE.Quaternion { - const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(elements, tdbFromUtc(jdUtc)); +export function bodyOrientation(elements: RotationalElements, jdUtc: number, target = new THREE.Quaternion(), followsUt = false): THREE.Quaternion { + const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(elements, followsUt ? jdUtc + TT_MINUS_UTC_DAYS : tdbFromUtc(jdUtc)); return poleFrame({ raDeg: poleRaDeg, decDeg: poleDecDeg }, target) .multiply(scratchTurn.setFromAxisAngle(Z_AXIS, primeMeridianDeg * DEG_TO_RAD)) .multiply(MAP_TO_BODY); @@ -100,6 +103,6 @@ export function bodyPageView(body: BodyRecord, bodies: readonly BodyRecord[], jd sun.set(-position.x, -position.y, -position.z).normalize().applyQuaternion(toPage); const turn = scratchPageTurn.setFromAxisAngle(Y_AXIS, sunAzimuthRad - Math.atan2(sun.x, sun.z)); sun.applyQuaternion(turn); - planet.copy(turn).multiply(toPage).multiply(bodyOrientation(elements, jdUtc, scratchBody)); + planet.copy(turn).multiply(toPage).multiply(bodyOrientation(elements, jdUtc, scratchBody, body.id === 'earth')); return true; }