From 8d8c65bdb2d1f41282c2dfcf311e552a22d73493 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 11:19:45 +0000 Subject: [PATCH] Propagate exoplanets with their real orbital period MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Every exoplanet was propagated with gmForParent(undefined) — the Sun's gravitational parameter — so the whole catalogue orbited as though each host were exactly one solar mass. Most hosts are red dwarfs far lighter than that, and a heavier central mass pulls harder and shortens the period, so their planets were whirling round much too fast: TRAPPIST-1 is 0.09 solar masses, and its planets were completing an orbit in roughly a third of the true time. pl_orbper was already in the TAP query and was being discarded on the way into the record. It is now kept, along with st_mass. A period and a semi-major axis together pin the host's gravitational parameter exactly, via GM = n^2 a^3 — no stellar model, no assumption, just the inverse of the orbitalPeriodDays helper that was already there. resolveGravitationalParameter picks the best available source: the measured period, else the published host mass, else one solar mass as before. A derived value implying something outside 0.01-150 solar masses is rejected and falls through, since a period and axis taken from disagreeing solutions would otherwise send a planet spinning at a visibly absurd rate. Note the direction of the error, which is the opposite of what it looks like: assuming a *heavier* host than reality makes a planet orbit *faster*. A test pins it, and caught me stating it backwards first. The NASA Exoplanet Archive is unreachable from this environment (egress policy returns 403 on CONNECT), so exoplanets.json cannot be regenerated here and still carries no periods. Behaviour is therefore unchanged until `npm run etl` is run somewhere with archive access, at which point every planet with a published period starts moving correctly with no further code changes. build.ts reports how many records gained a period, and rejects non-positive ones. Tests: 171 passing, up from 151, including a new end-to-end check that TRAPPIST-1 b with its real period completes exactly one orbit in 1.51088 days and sits a full diameter away at half that. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01WaySiNst4HhDXBHnMy8p5G --- .../system-orbits-renderer.spec.ts | 100 ++++++++++++++++++ .../galaxy-system/system-orbits-renderer.ts | 11 +- src/app/shared/astro/kepler.spec.ts | 89 ++++++++++++++++ src/app/shared/astro/kepler.ts | 68 +++++++++++- src/app/shared/models/exoplanet.model.ts | 8 ++ tools/etl/build.ts | 14 +++ tools/etl/fetchExoplanets.ts | 6 ++ 7 files changed, 293 insertions(+), 3 deletions(-) create mode 100644 src/app/features/galaxy-system/system-orbits-renderer.spec.ts diff --git a/src/app/features/galaxy-system/system-orbits-renderer.spec.ts b/src/app/features/galaxy-system/system-orbits-renderer.spec.ts new file mode 100644 index 0000000..1a257d9 --- /dev/null +++ b/src/app/features/galaxy-system/system-orbits-renderer.spec.ts @@ -0,0 +1,100 @@ +import * as THREE from 'three/webgpu'; +import { describe, expect, it } from 'vitest'; + +import { DEFAULT_EPOCH_JD } from '../../shared/astro/constants'; +import { ExoplanetRecord } from '../../shared/models/exoplanet.model'; +import { SystemOrbitsRenderer } from './system-orbits-renderer'; + +/** 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; + +function exoplanet(overrides: Partial = {}): ExoplanetRecord { + return { + id: 'TRAPPIST-1 b', + hostStarId: 1, + hostStarName: 'TRAPPIST-1', + name: 'TRAPPIST-1 b', + orbit: { semiMajorAxisAu: TRAPPIST_1B_SEMI_MAJOR_AXIS_AU, eccentricity: 0 }, + ...overrides + }; +} + +/** Marker position for the system's single exoplanet at a given Julian date. */ +function positionAt(renderer: SystemOrbitsRenderer, epochJd: number): THREE.Vector3 { + renderer.update(epochJd); + return renderer.members[0].marker.position.clone(); +} + +describe('SystemOrbitsRenderer exoplanet propagation', () => { + it('completes exactly one orbit over the measured period', () => { + // The end-to-end check that the period actually reaches the propagator: after one full + // published period the planet must be back where it started. + const renderer = new SystemOrbitsRenderer([], [exoplanet({ periodDays: TRAPPIST_1B_PERIOD_DAYS })]); + + const start = positionAt(renderer, DEFAULT_EPOCH_JD); + const afterOnePeriod = positionAt(renderer, DEFAULT_EPOCH_JD + TRAPPIST_1B_PERIOD_DAYS); + const afterHalfPeriod = positionAt(renderer, DEFAULT_EPOCH_JD + TRAPPIST_1B_PERIOD_DAYS / 2); + + expect(afterOnePeriod.distanceTo(start)).toBeLessThan(1e-6); + // Half an orbit of a circle is the far side, a full diameter away. + expect(afterHalfPeriod.distanceTo(start)).toBeCloseTo(2 * TRAPPIST_1B_SEMI_MAJOR_AXIS_AU, 6); + renderer.dispose(); + }); + + it('moves a red dwarf planet more slowly than the old solar-mass assumption did', () => { + // Assuming a solar-mass host made TRAPPIST-1's planets orbit about 3.3x too fast, so the + // corrected planet must have travelled less far after the same elapsed time. + const corrected = new SystemOrbitsRenderer([], [exoplanet({ periodDays: TRAPPIST_1B_PERIOD_DAYS })]); + const assumingSolar = new SystemOrbitsRenderer([], [exoplanet()]); + + const elapsed = TRAPPIST_1B_PERIOD_DAYS / 8; + const correctedTravel = positionAt(corrected, DEFAULT_EPOCH_JD).distanceTo(positionAt(corrected, DEFAULT_EPOCH_JD + elapsed)); + const solarTravel = positionAt(assumingSolar, DEFAULT_EPOCH_JD).distanceTo( + positionAt(assumingSolar, DEFAULT_EPOCH_JD + elapsed) + ); + + expect(correctedTravel).toBeLessThan(solarTravel); + corrected.dispose(); + assumingSolar.dispose(); + }); + + it('uses the host star mass when no period is published', () => { + const fromMass = new SystemOrbitsRenderer([], [exoplanet({ hostStarMassSolar: 0.0898 })]); + const fromPeriod = new SystemOrbitsRenderer([], [exoplanet({ periodDays: TRAPPIST_1B_PERIOD_DAYS })]); + + const elapsed = 0.3; + const massTravel = positionAt(fromMass, DEFAULT_EPOCH_JD).distanceTo(positionAt(fromMass, DEFAULT_EPOCH_JD + elapsed)); + const periodTravel = positionAt(fromPeriod, DEFAULT_EPOCH_JD).distanceTo(positionAt(fromPeriod, DEFAULT_EPOCH_JD + elapsed)); + + // The published mass and the period-derived mass agree, so the two must nearly coincide. + expect(massTravel).toBeCloseTo(periodTravel, 4); + fromMass.dispose(); + fromPeriod.dispose(); + }); + + it('still renders an exoplanet that has neither a period nor a host mass', () => { + const renderer = new SystemOrbitsRenderer([], [exoplanet()]); + + expect(renderer.members).toHaveLength(1); + expect(positionAt(renderer, DEFAULT_EPOCH_JD).length()).toBeCloseTo(TRAPPIST_1B_SEMI_MAJOR_AXIS_AU, 6); + renderer.dispose(); + }); + + it('skips an exoplanet with no semi-major axis rather than crashing', () => { + const renderer = new SystemOrbitsRenderer([], [exoplanet({ orbit: { eccentricity: 0 } })]); + + expect(renderer.members).toHaveLength(0); + renderer.dispose(); + }); + + it('keeps every propagated position finite', () => { + const renderer = new SystemOrbitsRenderer([], [exoplanet({ periodDays: TRAPPIST_1B_PERIOD_DAYS, orbit: { semiMajorAxisAu: TRAPPIST_1B_SEMI_MAJOR_AXIS_AU, eccentricity: 0.62 } })]); + + for (const offset of [0, 0.1, 1, 10, 1000]) { + const { x, y, z } = positionAt(renderer, DEFAULT_EPOCH_JD + offset); + expect([x, y, z].every(Number.isFinite)).toBe(true); + } + renderer.dispose(); + }); +}); diff --git a/src/app/features/galaxy-system/system-orbits-renderer.ts b/src/app/features/galaxy-system/system-orbits-renderer.ts index 20b9b01..f6c3623 100644 --- a/src/app/features/galaxy-system/system-orbits-renderer.ts +++ b/src/app/features/galaxy-system/system-orbits-renderer.ts @@ -1,7 +1,7 @@ import * as THREE from 'three/webgpu'; import { gmForParent } from '../../shared/astro/constants'; -import { orbitEllipsePoints, propagateOrbit, resolveOrbitalElements } from '../../shared/astro/kepler'; +import { orbitEllipsePoints, propagateOrbit, resolveGravitationalParameter, resolveOrbitalElements } from '../../shared/astro/kepler'; import { BodyRecord, OrbitalElements } from '../../shared/models/body.model'; import { ExoplanetRecord } from '../../shared/models/exoplanet.model'; @@ -160,7 +160,14 @@ export class SystemOrbitsRenderer { epochJd: exoplanet.orbit.epochJd }); const radiusKm = exoplanet.radiusEarth ? exoplanet.radiusEarth * EARTH_RADIUS_KM : undefined; - const tracked = this.addTopLevelBody(exoplanet.id, 'exoplanet', elements, gmForParent(undefined), radiusKm); + // 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. + const gm = resolveGravitationalParameter({ + semiMajorAxisAu: exoplanet.orbit.semiMajorAxisAu, + periodDays: exoplanet.periodDays, + hostStarMassSolar: exoplanet.hostStarMassSolar + }); + const tracked = this.addTopLevelBody(exoplanet.id, 'exoplanet', elements, gm, radiusKm); members.push({ id: exoplanet.id, kind: 'exoplanet', marker: tracked.marker }); } diff --git a/src/app/shared/astro/kepler.spec.ts b/src/app/shared/astro/kepler.spec.ts index 99a1d6f..6b29761 100644 --- a/src/app/shared/astro/kepler.spec.ts +++ b/src/app/shared/astro/kepler.spec.ts @@ -2,11 +2,13 @@ import { describe, expect, it } from 'vitest'; import { GM_SUN_AU3_PER_DAY2, DEFAULT_EPOCH_JD } from './constants'; import { + gravitationalParameterFromPeriod, meanMotionRadPerDay, orbitEllipsePoints, orbitalPeriodDays, positionAtTrueAnomaly, propagateOrbit, + resolveGravitationalParameter, resolveOrbitalElements, solveEccentricAnomaly, trueAnomalyFromEccentricAnomaly @@ -163,3 +165,90 @@ describe('resolveOrbitalElements', () => { expect(resolved.argumentOfPeriapsisDeg).toBe(50); }); }); + +describe('gravitationalParameterFromPeriod', () => { + it('round-trips with orbitalPeriodDays', () => { + const derived = gravitationalParameterFromPeriod(1, 365.256); + expect(orbitalPeriodDays(1, derived)).toBeCloseTo(365.256, 9); + }); + + it('recovers the Sun from Earth\'s orbit', () => { + // 1 AU in one sidereal year is the definition of the solar gravitational parameter. + const derived = gravitationalParameterFromPeriod(1, 365.256363); + expect(derived / GM_SUN_AU3_PER_DAY2).toBeCloseTo(1, 4); + }); + + it('recovers a red dwarf from a real short-period orbit', () => { + // TRAPPIST-1 b: 0.01154 AU in 1.51088 days around a 0.0898 solar-mass star. + const derived = gravitationalParameterFromPeriod(0.01154, 1.51088); + expect(derived / GM_SUN_AU3_PER_DAY2).toBeCloseTo(0.09, 2); + }); + + it('scales as a^3 at fixed period', () => { + const single = gravitationalParameterFromPeriod(1, 100); + const doubled = gravitationalParameterFromPeriod(2, 100); + expect(doubled / single).toBeCloseTo(8, 9); + }); +}); + +describe('resolveGravitationalParameter', () => { + it('prefers the measured period over everything else', () => { + // The period says 0.09 solar masses; the (deliberately wrong) host mass says 5. + const gm = resolveGravitationalParameter({ semiMajorAxisAu: 0.01154, periodDays: 1.51088, hostStarMassSolar: 5 }); + expect(gm / GM_SUN_AU3_PER_DAY2).toBeCloseTo(0.09, 2); + }); + + it('corrects a red dwarf planet that the solar-mass assumption spun too fast', () => { + const withPeriod = resolveGravitationalParameter({ semiMajorAxisAu: 0.01154, periodDays: 1.51088 }); + const assumingSolar = resolveGravitationalParameter({ semiMajorAxisAu: 0.01154 }); + + // A heavier central mass pulls harder, so it shortens the period: T scales as 1/sqrt(GM). + // Assuming the Sun for TRAPPIST-1's 0.09 solar masses therefore made its planets orbit + // sqrt(0.09) = 0.3x the true period — about 3.3x too fast, not too slow. + expect(orbitalPeriodDays(0.01154, withPeriod)).toBeCloseTo(1.51088, 4); + const ratio = orbitalPeriodDays(0.01154, assumingSolar) / orbitalPeriodDays(0.01154, withPeriod); + expect(ratio).toBeCloseTo(Math.sqrt(0.09), 2); + }); + + it('falls back to the host star mass when no period is published', () => { + const gm = resolveGravitationalParameter({ semiMajorAxisAu: 0.5, hostStarMassSolar: 0.31 }); + expect(gm).toBeCloseTo(GM_SUN_AU3_PER_DAY2 * 0.31, 12); + }); + + it('falls back to one solar mass when nothing is known', () => { + expect(resolveGravitationalParameter({ semiMajorAxisAu: 1 })).toBe(GM_SUN_AU3_PER_DAY2); + }); + + it('ignores a period that is missing, zero, negative or not a number', () => { + for (const periodDays of [undefined, 0, -5, Number.NaN, Number.POSITIVE_INFINITY]) { + expect(resolveGravitationalParameter({ semiMajorAxisAu: 1, periodDays })).toBe(GM_SUN_AU3_PER_DAY2); + } + }); + + it('ignores a host mass that is missing, zero or negative', () => { + for (const hostStarMassSolar of [undefined, 0, -1, Number.NaN]) { + expect(resolveGravitationalParameter({ semiMajorAxisAu: 1, hostStarMassSolar })).toBe(GM_SUN_AU3_PER_DAY2); + } + }); + + it('rejects a period implying something that cannot be a star, and falls through', () => { + // 1 AU in a single day implies thousands of solar masses. + const gm = resolveGravitationalParameter({ semiMajorAxisAu: 1, periodDays: 1, hostStarMassSolar: 0.5 }); + expect(gm).toBeCloseTo(GM_SUN_AU3_PER_DAY2 * 0.5, 12); + }); + + it('rejects a period implying far too little mass', () => { + // 1 AU taking a million days implies a mass far below any star. + expect(resolveGravitationalParameter({ semiMajorAxisAu: 1, periodDays: 1e6 })).toBe(GM_SUN_AU3_PER_DAY2); + }); + + it('rejects an implausible host mass too', () => { + expect(resolveGravitationalParameter({ semiMajorAxisAu: 1, hostStarMassSolar: 5000 })).toBe(GM_SUN_AU3_PER_DAY2); + }); + + it('keeps a real short-period hot Jupiter around a sun-like star', () => { + // 51 Pegasi b: 0.0527 AU in 4.23 days around a ~1.1 solar-mass star. + const gm = resolveGravitationalParameter({ semiMajorAxisAu: 0.0527, periodDays: 4.230785 }); + expect(gm / GM_SUN_AU3_PER_DAY2).toBeCloseTo(1.1, 1); + }); +}); diff --git a/src/app/shared/astro/kepler.ts b/src/app/shared/astro/kepler.ts index a7a8783..d6dda4c 100644 --- a/src/app/shared/astro/kepler.ts +++ b/src/app/shared/astro/kepler.ts @@ -1,5 +1,5 @@ import { CartesianCoordinates } from './coordinates'; -import { DEFAULT_EPOCH_JD } from './constants'; +import { DEFAULT_EPOCH_JD, GM_SUN_AU3_PER_DAY2 } from './constants'; import { OrbitalElements } from '../models/body.model'; const DEG_TO_RAD = Math.PI / 180; @@ -34,6 +34,72 @@ export function orbitalPeriodDays(semiMajorAxisAu: number, gmAu3PerDay2: number) return TWO_PI / meanMotionRadPerDay(semiMajorAxisAu, gmAu3PerDay2); } +/** + * Gravitational parameter implied by a measured orbital period — the inverse of + * {@link orbitalPeriodDays}: `GM = n^2 * a^3`, with `n = 2*pi / T`. + * + * This is how an exoplanet's host star gets its mass into the propagator. Nothing about the + * star needs to be known or guessed: the period and the semi-major axis between them pin the + * gravitational parameter exactly. + */ +export function gravitationalParameterFromPeriod(semiMajorAxisAu: number, periodDays: number): number { + const meanMotion = TWO_PI / periodDays; + return meanMotion * meanMotion * semiMajorAxisAu * semiMajorAxisAu * semiMajorAxisAu; +} + +/** + * Plausible range for a host star's mass, in solar masses — from below the hydrogen-burning + * limit to beyond the heaviest known stars. Used only to reject a derived value that cannot be + * a star, which would otherwise send a planet spinning at a visibly absurd rate. + */ +const MIN_PLAUSIBLE_STELLAR_MASS_SOLAR = 0.01; +const MAX_PLAUSIBLE_STELLAR_MASS_SOLAR = 150; + +function isPositiveFinite(value: number | undefined): value is number { + return value !== undefined && Number.isFinite(value) && value > 0; +} + +/** + * The gravitational parameter to propagate a planet with, in AU^3/day^2, best source first: + * + * 1. **Its measured orbital period.** Exact, and independent of any stellar model. + * 2. **Its host star's measured mass.** + * 3. **One solar mass**, as a last resort. + * + * Falling back to the Sun is a real approximation, not a neutral default. Most exoplanet hosts + * are red dwarfs far lighter than the Sun, and a heavier central mass pulls harder and shortens + * the period, so assuming solar mass makes their planets whirl round far too fast. TRAPPIST-1 + * is 0.09 solar masses; its planets were completing an orbit in roughly a third of the true + * time. + */ +export function resolveGravitationalParameter(input: { + semiMajorAxisAu: number; + periodDays?: number; + hostStarMassSolar?: number; +}): number { + const { semiMajorAxisAu, periodDays, hostStarMassSolar } = input; + + if (isPositiveFinite(periodDays) && isPositiveFinite(semiMajorAxisAu)) { + const derived = gravitationalParameterFromPeriod(semiMajorAxisAu, periodDays); + const impliedMassSolar = derived / GM_SUN_AU3_PER_DAY2; + // A period and axis drawn from disagreeing solutions can imply something that is not a + // star; prefer a known-approximate answer over a confidently wrong one. + if (impliedMassSolar >= MIN_PLAUSIBLE_STELLAR_MASS_SOLAR && impliedMassSolar <= MAX_PLAUSIBLE_STELLAR_MASS_SOLAR) { + return derived; + } + } + + if ( + isPositiveFinite(hostStarMassSolar) && + hostStarMassSolar >= MIN_PLAUSIBLE_STELLAR_MASS_SOLAR && + hostStarMassSolar <= MAX_PLAUSIBLE_STELLAR_MASS_SOLAR + ) { + return GM_SUN_AU3_PER_DAY2 * hostStarMassSolar; + } + + return GM_SUN_AU3_PER_DAY2; +} + /** Normalizes an angle (radians) into [0, 2*pi). */ function normalizeAngle(angleRad: number): number { const wrapped = angleRad % TWO_PI; diff --git a/src/app/shared/models/exoplanet.model.ts b/src/app/shared/models/exoplanet.model.ts index 3ff4ea2..7260017 100644 --- a/src/app/shared/models/exoplanet.model.ts +++ b/src/app/shared/models/exoplanet.model.ts @@ -12,5 +12,13 @@ export interface ExoplanetRecord { radiusEarth?: number; massEarth?: number; discoveryYear?: number; + /** + * Measured orbital period in days (`pl_orbper`). Together with the semi-major axis this + * pins the host star's gravitational parameter exactly, so the planet can be propagated at + * its real rate instead of as though it orbited the Sun — see `resolveGravitationalParameter`. + */ + periodDays?: number; + /** Host star mass in solar masses (`st_mass`); the fallback when no period is published. */ + hostStarMassSolar?: number; orbit: Partial; } diff --git a/tools/etl/build.ts b/tools/etl/build.ts index 05636b9..9eea63a 100644 --- a/tools/etl/build.ts +++ b/tools/etl/build.ts @@ -70,9 +70,23 @@ function validateExoplanets(exoplanets: ExoplanetRecord[], starIds: Set) ); crossReferenced++; } + + assertCondition( + exoplanet.periodDays === undefined || exoplanet.periodDays > 0, + `Exoplanet ${exoplanet.id} has a non-positive orbital period.` + ); + assertCondition( + exoplanet.hostStarMassSolar === undefined || exoplanet.hostStarMassSolar > 0, + `Exoplanet ${exoplanet.id} has a non-positive host star mass.` + ); } console.log(` ${crossReferenced}/${exoplanets.length} exoplanets cross-referenced to a HYG host star.`); + + // How many can be propagated at their real rate rather than as if the host were the Sun. + const withPeriod = exoplanets.filter((exoplanet) => exoplanet.periodDays !== undefined).length; + const withHostMass = exoplanets.filter((exoplanet) => exoplanet.hostStarMassSolar !== undefined).length; + console.log(` ${withPeriod}/${exoplanets.length} have a measured period, ${withHostMass} a host star mass.`); } const UNIT_VECTOR_TOLERANCE = 1e-6; diff --git a/tools/etl/fetchExoplanets.ts b/tools/etl/fetchExoplanets.ts index c579d93..55d6335 100644 --- a/tools/etl/fetchExoplanets.ts +++ b/tools/etl/fetchExoplanets.ts @@ -22,6 +22,7 @@ const TAP_COLUMNS = [ 'pl_orbper', 'pl_rade', 'pl_bmasse', + 'st_mass', 'disc_year' ].join(','); const TAP_QUERY = `select+${TAP_COLUMNS}+from+ps+where+default_flag=1&format=csv`; @@ -69,6 +70,11 @@ export async function fetchExoplanets(stars?: StarRecord[]): Promise