Take TT - UT from the historical record before 1972, so the far dates the clock reaches turn every body by the right amount
The clock reaches AD 1, but TT - UT was held at today's 69.184 s. At AD 1000 it was 1 574 s and at AD 1 about 10 570 (Espenak and Meeus, NASA's Five Millennium Canon; Horizons' TDB - UT gives 1 658 and 10 466 on JD 2086455 and 1721600). So every spin but Earth's was (ΔT - 69 s) times its rate out, Jupiter 15.2 degrees at AD 1000 and 106 at AD 1, Mars 6 and 43, and every orbit that much behind: the Moon about 0.2 and 1.4 degrees. ttMinusUtSeconds gives TT - UT for a date on the clock: the Espenak-Meeus polynomials before 1972, 32.184 s plus UTC's leap seconds from 1972 to the last one, at the start of 2017, and 69.184 s held after it, as Horizons holds it. Its pieces join within 0.1 s. tdbFromUtc, which positions and spins already share, now adds it. Within 0.2 s of Horizons in 1950, 105 s at AD 1 and 86 s at AD 1000, where the historical record itself is that uncertain. Earth is the exception: its turning is what UT counts, so the clock's date already says how far it has turned, and ΔT would turn it again, 44 degrees at AD 1. Its W, fitted to today, keeps today's 69.184 s (bodyOrientation's followsUt, set for Earth in the system view and on its page). In the running app at 1000-01-01 00:00 UT, Jupiter's drawn prime meridian sits 0.000 degrees from its IAU W at TT and 15.164 from where the held offset put it; Earth's sits on its W at UT + 69.184 s, 6.288 degrees short of what TT would have turned it to. The renderer spec now hands its frozen Horizons vectors over as the UT dates that name them through the same TT - UT, and checks Jupiter's and Earth's prime meridians at AD 1000. Controls: the leap-second rule used before 1972 fails "follows the historical record before 1972"; TT - UT held at 69 s fails "turns Jupiter at AD 1000 by its W"; Earth turned at TDB, or the renderer or the page not keeping it on UT, fails the Earth tests. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
This commit is contained in:
@@ -1,13 +1,17 @@
|
|||||||
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, 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 { 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 { BodyRecord, RotationalElements } from '../../shared/models/body.model';
|
||||||
import { ExoplanetRecord } from '../../shared/models/exoplanet.model';
|
import { ExoplanetRecord } from '../../shared/models/exoplanet.model';
|
||||||
import { SystemOrbitsRenderer } from './system-orbits-renderer';
|
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. */
|
/** 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_SEMI_MAJOR_AXIS_AU = 0.01154;
|
||||||
const TRAPPIST_1B_PERIOD_DAYS = 1.51088;
|
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.
|
// 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], []);
|
||||||
// The clock's UTC date whose TDB is the elements' epoch.
|
// The clock's UT date whose TDB is the elements' epoch.
|
||||||
renderer.update(DEFAULT_EPOCH_JD - TT_MINUS_UTC_DAYS);
|
renderer.update(utOf(DEFAULT_EPOCH_JD));
|
||||||
|
|
||||||
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);
|
||||||
@@ -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 —
|
// 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. Horizons' dates are TDB and the renderer's are the
|
// 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<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}}},
|
||||||
@@ -555,7 +559,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 - TT_MINUS_UTC_DAYS);
|
renderer.update(utOf(jd));
|
||||||
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);
|
||||||
@@ -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', () => {
|
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).
|
// 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 - TT_MINUS_UTC_DAYS);
|
renderer.update(utOf(2488069.5));
|
||||||
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);
|
||||||
@@ -656,6 +660,33 @@ describe('solar-system bodies against Horizons', () => {
|
|||||||
/** Degrees between two longitudes, the short way round. */
|
/** Degrees between two longitudes, the short way round. */
|
||||||
const apart = (a: number, b: number): number => Math.abs(((((a - b) % 360) + 540) % 360) - 180);
|
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', () => {
|
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,
|
// 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.
|
// and the drawn sphere has it over 0.43 W.
|
||||||
|
|||||||
@@ -480,7 +480,7 @@ export class SystemOrbitsRenderer {
|
|||||||
body.marker.position.copy(body.position);
|
body.marker.position.copy(body.position);
|
||||||
orientOrbit(body.orbitLine.quaternion, current, body.frame);
|
orientOrbit(body.orbitLine.quaternion, current, body.frame);
|
||||||
if (body.rotationalElements) {
|
if (body.rotationalElements) {
|
||||||
bodyOrientation(body.rotationalElements, epochJd, body.marker.quaternion);
|
bodyOrientation(body.rotationalElements, epochJd, body.marker.quaternion, body.id === 'earth');
|
||||||
} else if (body.rotationPeriodHours) {
|
} else if (body.rotationPeriodHours) {
|
||||||
body.marker.quaternion.copy(spinFor(current, body.frame, body.rotationPeriodHours, jdTdb - body.elements.epochJd));
|
body.marker.quaternion.copy(spinFor(current, body.frame, body.rotationPeriodHours, jdTdb - body.elements.epochJd));
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -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);
|
||||||
|
});
|
||||||
|
});
|
||||||
@@ -20,22 +20,86 @@ export const DEFAULT_EPOCH_JD = 2451545.0;
|
|||||||
export const GM_SUN_AU3_PER_DAY2 = 0.01720209895 * 0.01720209895;
|
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
|
* TT - UTC today, in days: 32.184 s plus the 37 leap seconds UTC has taken since 1972, the last at
|
||||||
* end of 2016. TDB, which ephemerides run on, stays within 2 ms of TT. Held constant, as Horizons
|
* the end of 2016. TDB, which ephemerides run on, stays within 2 ms of TT. See {@link ttMinusUtSeconds}
|
||||||
* holds it for dates past the last announced leap second; before 2017 it was smaller, about 29 s
|
* for other dates.
|
||||||
* in 1950.
|
|
||||||
*/
|
*/
|
||||||
export const TT_MINUS_UTC_DAYS = 69.184 / 86400;
|
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
|
* 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
|
* 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,
|
* 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.
|
* which is 0.9 degrees of Phobos's orbit and 0.16 of Io's.
|
||||||
*/
|
*/
|
||||||
export function tdbFromUtc(jdUtc: number): number {
|
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". */
|
/** Converts a JS `Date` into a Julian date (days), for driving the Kepler propagator "now". */
|
||||||
|
|||||||
@@ -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 { TT_MINUS_UTC_DAYS } from '../astro/constants';
|
import { tdbFromUtc } from '../astro/constants';
|
||||||
import { eclipticToEquatorial } from '../astro/coordinates';
|
import { eclipticToEquatorial } from '../astro/coordinates';
|
||||||
import { meanElementsAt, positionAtEpoch } from '../astro/kepler';
|
import { meanElementsAt, positionAtEpoch } from '../astro/kepler';
|
||||||
import { BodyRecord } from '../models/body.model';
|
import { BodyRecord } from '../models/body.model';
|
||||||
@@ -54,15 +54,18 @@ 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', () => {
|
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, and the Sun seen from Earth's mean
|
// Earth's own sphere, turned as the system view turns it (by UT, see `bodyOrientation`), and
|
||||||
// place at the clock's date taken to TDB: the page must light that same point of its map.
|
// the Sun seen from Earth's mean place at the clock's date taken to TDB: the page must light that
|
||||||
const planet = new THREE.Quaternion();
|
// same point of its map, today and at AD 1000, when TT was 1 574 s past UT.
|
||||||
const sun = new THREE.Vector3();
|
for (const jdUt of [JUNE_1_2025_NOON_UTC, 2086307.5]) {
|
||||||
bodyPageView(EARTH, BODIES, JUNE_1_2025_NOON_UTC, SUN_AZIMUTH, planet, sun);
|
const planet = new THREE.Quaternion();
|
||||||
const place = eclipticToEquatorial(positionAtEpoch(meanElementsAt(EARTH.orbit, EARTH.rates, JUNE_1_2025_NOON_UTC + TT_MINUS_UTC_DAYS)));
|
const sun = new THREE.Vector3();
|
||||||
const expected = new THREE.Vector3(-place.x, -place.y, -place.z).normalize().applyQuaternion(bodyOrientation(EARTH.rotationalElements!, JUNE_1_2025_NOON_UTC).invert());
|
bodyPageView(EARTH, BODIES, jdUt, SUN_AZIMUTH, planet, sun);
|
||||||
expect(sun.clone().applyQuaternion(planet.clone().invert()).angleTo(expected)).toBeLessThan(1e-9);
|
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', () => {
|
it('keeps the pole up and the Sun where the page’s light stands, turning the body under it', () => {
|
||||||
|
|||||||
@@ -1,6 +1,6 @@
|
|||||||
import * as THREE from 'three/webgpu';
|
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 { 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';
|
||||||
@@ -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,
|
* 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.
|
* 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 UT and the IAU's elements run on TDB, 69.184 s ahead today and 1 574 s at AD 1000;
|
||||||
* 0.29 degrees, Jupiter 0.70 and Phobos 0.90, so the date is taken to TDB here (see `tdbFromUtc`).
|
* 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 {
|
export function bodyOrientation(elements: RotationalElements, jdUtc: number, target = new THREE.Quaternion(), followsUt = false): THREE.Quaternion {
|
||||||
const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(elements, tdbFromUtc(jdUtc));
|
const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(elements, followsUt ? jdUtc + TT_MINUS_UTC_DAYS : 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);
|
||||||
@@ -100,6 +103,6 @@ export function bodyPageView(body: BodyRecord, bodies: readonly BodyRecord[], jd
|
|||||||
sun.set(-position.x, -position.y, -position.z).normalize().applyQuaternion(toPage);
|
sun.set(-position.x, -position.y, -position.z).normalize().applyQuaternion(toPage);
|
||||||
const turn = scratchPageTurn.setFromAxisAngle(Y_AXIS, sunAzimuthRad - Math.atan2(sun.x, sun.z));
|
const turn = scratchPageTurn.setFromAxisAngle(Y_AXIS, sunAzimuthRad - Math.atan2(sun.x, sun.z));
|
||||||
sun.applyQuaternion(turn);
|
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;
|
return true;
|
||||||
}
|
}
|
||||||
|
|||||||
Reference in New Issue
Block a user