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