@
Add star-map Angular app, ETL pipeline, and caveman plugin Angular 3D star map (galaxy/system/body views, Three.js rendering, navigation store) plus the NASA ETL tooling that builds the star, exoplanet and solar-system datasets, Playwright e2e suite, and the cs:caveman Claude Code plugin (command, agent, skill). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> @
This commit is contained in:
@@ -0,0 +1,52 @@
|
||||
/** Astronomical unit conversion and gravitational constants shared by the astro math modules. */
|
||||
|
||||
/** Number of astronomical units in one parsec (IAU exact definition). */
|
||||
export const AU_PER_PARSEC = 206264.80624709636;
|
||||
|
||||
/**
|
||||
* Reference epoch (Julian date, J2000.0) used when orbital data lacks an explicit epoch —
|
||||
* e.g. exoplanets from the NASA Exoplanet Archive only report a handful of elements
|
||||
* (semi-major axis, eccentricity, sometimes argument of periapsis), not a mean-anomaly/epoch
|
||||
* pair. Defaulting the missing epoch to J2000 still lets the body's real orbital period
|
||||
* carry it around a plausible (if not phase-accurate) orbit over time.
|
||||
*/
|
||||
export const DEFAULT_EPOCH_JD = 2451545.0;
|
||||
|
||||
/**
|
||||
* Heliocentric gravitational parameter (GM of the Sun), in AU^3/day^2 — the square of the
|
||||
* Gaussian gravitational constant `k = 0.01720209895 rad/day`. Used to derive a body's mean
|
||||
* motion from its semi-major axis via Kepler's third law.
|
||||
*/
|
||||
export const GM_SUN_AU3_PER_DAY2 = 0.01720209895 * 0.01720209895;
|
||||
|
||||
/**
|
||||
* Approximate planet/Sun mass ratios for the major planets that host moons in `bodies.json`.
|
||||
* Used to derive each planet's gravitational parameter (for propagating its moons) as
|
||||
* `GM_SUN_AU3_PER_DAY2 * massRatio`. Precise enough for visualization; not JPL-grade.
|
||||
*/
|
||||
const PLANET_TO_SUN_MASS_RATIO: Record<string, number> = {
|
||||
earth: 3.003e-6,
|
||||
mars: 3.227e-7,
|
||||
jupiter: 9.545e-4,
|
||||
saturn: 2.857e-4,
|
||||
uranus: 4.365e-5,
|
||||
neptune: 5.151e-5
|
||||
};
|
||||
|
||||
/**
|
||||
* Gravitational parameter (AU^3/day^2) to use when propagating a body's orbit: the Sun's
|
||||
* for planets/dwarfs/exoplanets, or the host planet's (derived from its Sun mass ratio) for
|
||||
* moons. Falls back to the Sun's GM if `parentBodyId` isn't a known planet.
|
||||
*/
|
||||
export function gmForParent(parentBodyId: string | undefined): number {
|
||||
if (!parentBodyId) {
|
||||
return GM_SUN_AU3_PER_DAY2;
|
||||
}
|
||||
const massRatio = PLANET_TO_SUN_MASS_RATIO[parentBodyId];
|
||||
return massRatio ? GM_SUN_AU3_PER_DAY2 * massRatio : GM_SUN_AU3_PER_DAY2;
|
||||
}
|
||||
|
||||
/** Converts a JS `Date` into a Julian date (days), for driving the Kepler propagator "now". */
|
||||
export function dateToJulianDate(date: Date = new Date()): number {
|
||||
return date.getTime() / 86400000 + 2440587.5;
|
||||
}
|
||||
@@ -0,0 +1,67 @@
|
||||
import { describe, expect, it } from 'vitest';
|
||||
|
||||
import { distanceBetween, parallaxMasToParsecs, raDecDistanceToXyz, raDegDecDistanceToXyz } from './coordinates';
|
||||
|
||||
// Reference values taken directly from the HYG v4.1 database (RA/Dec/dist and its own
|
||||
// precomputed x/y/z, which uses the same equatorial-Cartesian convention we implement).
|
||||
describe('raDecDistanceToXyz', () => {
|
||||
it('matches the HYG reference position for Sirius', () => {
|
||||
const result = raDecDistanceToXyz(6.752481, -16.716116, 2.6371);
|
||||
|
||||
expect(result.x).toBeCloseTo(-0.494323, 3);
|
||||
expect(result.y).toBeCloseTo(2.476731, 3);
|
||||
expect(result.z).toBeCloseTo(-0.758485, 3);
|
||||
});
|
||||
|
||||
it('matches the HYG reference position for Proxima Centauri', () => {
|
||||
const result = raDecDistanceToXyz(14.495985, -62.679485, 1.2959);
|
||||
|
||||
expect(result.x).toBeCloseTo(-0.472264, 3);
|
||||
expect(result.y).toBeCloseTo(-0.361451, 3);
|
||||
expect(result.z).toBeCloseTo(-1.151219, 3);
|
||||
});
|
||||
|
||||
it('places a star on RA 6h / Dec 0 entirely on the +Y axis', () => {
|
||||
const result = raDecDistanceToXyz(6, 0, 10);
|
||||
|
||||
expect(result.x).toBeCloseTo(0, 9);
|
||||
expect(result.y).toBeCloseTo(10, 9);
|
||||
expect(result.z).toBeCloseTo(0, 9);
|
||||
});
|
||||
|
||||
it('places the vernal equinox direction entirely on the +X axis', () => {
|
||||
const result = raDecDistanceToXyz(0, 0, 10);
|
||||
|
||||
expect(result.x).toBeCloseTo(10, 9);
|
||||
expect(result.y).toBeCloseTo(0, 9);
|
||||
expect(result.z).toBeCloseTo(0, 9);
|
||||
});
|
||||
});
|
||||
|
||||
describe('raDegDecDistanceToXyz', () => {
|
||||
it('is equivalent to raDecDistanceToXyz with RA converted from degrees to hours', () => {
|
||||
const fromHours = raDecDistanceToXyz(6.752481, -16.716116, 2.6371);
|
||||
const fromDegrees = raDegDecDistanceToXyz(6.752481 * 15, -16.716116, 2.6371);
|
||||
|
||||
expect(fromDegrees.x).toBeCloseTo(fromHours.x, 9);
|
||||
expect(fromDegrees.y).toBeCloseTo(fromHours.y, 9);
|
||||
expect(fromDegrees.z).toBeCloseTo(fromHours.z, 9);
|
||||
});
|
||||
});
|
||||
|
||||
describe('parallaxMasToParsecs', () => {
|
||||
it('converts a positive parallax to the expected distance', () => {
|
||||
expect(parallaxMasToParsecs(769.33)).toBeCloseTo(1.3, 2); // Proxima Centauri
|
||||
});
|
||||
|
||||
it('returns Infinity for zero or negative parallax', () => {
|
||||
expect(parallaxMasToParsecs(0)).toBe(Infinity);
|
||||
expect(parallaxMasToParsecs(-5)).toBe(Infinity);
|
||||
});
|
||||
});
|
||||
|
||||
describe('distanceBetween', () => {
|
||||
it('computes the Euclidean distance between two points', () => {
|
||||
expect(distanceBetween({ x: 0, y: 0, z: 0 }, { x: 3, y: 4, z: 0 })).toBeCloseTo(5, 9);
|
||||
});
|
||||
});
|
||||
@@ -0,0 +1,53 @@
|
||||
const DEG_TO_RAD = Math.PI / 180;
|
||||
const HOURS_TO_DEG = 15;
|
||||
|
||||
export interface CartesianCoordinates {
|
||||
x: number;
|
||||
y: number;
|
||||
z: number;
|
||||
}
|
||||
|
||||
/**
|
||||
* Converts right ascension (hours), declination (degrees) and distance (parsecs) into
|
||||
* equatorial Cartesian coordinates (parsecs). Matches the HYG database convention:
|
||||
* +X toward the vernal equinox (epoch 2000), +Z toward the north celestial pole,
|
||||
* +Y toward RA 6h / Dec 0.
|
||||
*/
|
||||
export function raDecDistanceToXyz(raHours: number, decDeg: number, distancePc: number): CartesianCoordinates {
|
||||
const raRad = raHours * HOURS_TO_DEG * DEG_TO_RAD;
|
||||
const decRad = decDeg * DEG_TO_RAD;
|
||||
const cosDec = Math.cos(decRad);
|
||||
|
||||
return {
|
||||
x: distancePc * cosDec * Math.cos(raRad),
|
||||
y: distancePc * cosDec * Math.sin(raRad),
|
||||
z: distancePc * Math.sin(decRad)
|
||||
};
|
||||
}
|
||||
|
||||
/**
|
||||
* Same conversion as {@link raDecDistanceToXyz}, but for sources (e.g. exoplanet host
|
||||
* stars, deep-sky catalogs) that report right ascension in degrees rather than hours.
|
||||
*/
|
||||
export function raDegDecDistanceToXyz(raDeg: number, decDeg: number, distancePc: number): CartesianCoordinates {
|
||||
return raDecDistanceToXyz(raDeg / HOURS_TO_DEG, decDeg, distancePc);
|
||||
}
|
||||
|
||||
/**
|
||||
* Converts a parallax (milliarcseconds) into a distance in parsecs.
|
||||
* Returns `Infinity` for non-positive parallax (unmeasured/negative parallax).
|
||||
*/
|
||||
export function parallaxMasToParsecs(parallaxMas: number): number {
|
||||
return parallaxMas > 0 ? 1000 / parallaxMas : Infinity;
|
||||
}
|
||||
|
||||
/**
|
||||
* Euclidean distance (parsecs) between two Cartesian points, e.g. for nearest-neighbour
|
||||
* star matching during exoplanet host-star cross-referencing.
|
||||
*/
|
||||
export function distanceBetween(a: CartesianCoordinates, b: CartesianCoordinates): number {
|
||||
const dx = a.x - b.x;
|
||||
const dy = a.y - b.y;
|
||||
const dz = a.z - b.z;
|
||||
return Math.sqrt(dx * dx + dy * dy + dz * dz);
|
||||
}
|
||||
@@ -0,0 +1,64 @@
|
||||
import { describe, expect, it } from 'vitest';
|
||||
|
||||
import { buildStarNameIndex, normalizeStarName, resolveHostStarId } from './host-star-matching';
|
||||
import { StarRecord } from '../models/star.model';
|
||||
|
||||
// A small fixture standing in for a slice of the HYG star index, used to exercise the
|
||||
// exoplanet host-star cross-referencing logic without hitting any real API.
|
||||
const FIXTURE_STARS: StarRecord[] = [
|
||||
{ id: 1, name: 'Proxima Centauri', x: -0.472264, y: -0.361451, z: -1.151219, magnitude: 11.01, spectralType: 'M5Ve', colorIndex: 1.807 },
|
||||
{ id: 2, name: 'Sirius', x: -0.494323, y: 2.476731, z: -0.758485, magnitude: -1.44, spectralType: 'A0m...', colorIndex: 0.009 },
|
||||
{ id: 3, name: 'GJ 3512', x: 3.0, y: 4.0, z: 5.0, magnitude: 11.0, spectralType: 'M5.5', colorIndex: 1.6 }
|
||||
];
|
||||
|
||||
describe('normalizeStarName', () => {
|
||||
it('lowercases and strips non-alphanumeric characters', () => {
|
||||
expect(normalizeStarName('GJ 3512')).toBe('gj3512');
|
||||
expect(normalizeStarName('Proxima Centauri')).toBe('proximacentauri');
|
||||
});
|
||||
});
|
||||
|
||||
describe('resolveHostStarId', () => {
|
||||
it('matches by exact (normalized) host star name', () => {
|
||||
const id = resolveHostStarId({ hostname: 'Proxima Centauri', raDeg: NaN, decDeg: NaN, distancePc: NaN }, FIXTURE_STARS, 0.5);
|
||||
|
||||
expect(id).toBe(1);
|
||||
});
|
||||
|
||||
it('matches by name regardless of case/spacing differences', () => {
|
||||
const id = resolveHostStarId({ hostname: 'gj3512', raDeg: NaN, decDeg: NaN, distancePc: NaN }, FIXTURE_STARS, 0.5);
|
||||
|
||||
expect(id).toBe(3);
|
||||
});
|
||||
|
||||
it('falls back to nearest-neighbour position matching when the name is unknown', () => {
|
||||
// Slightly off from Sirius's exact position, within tolerance.
|
||||
const id = resolveHostStarId({ hostname: 'Sirius A', raDeg: 101.29, decDeg: -16.72, distancePc: 2.64 }, FIXTURE_STARS, 0.5);
|
||||
|
||||
expect(id).toBe(2);
|
||||
});
|
||||
|
||||
it('returns null when no name match and no star is within tolerance', () => {
|
||||
const id = resolveHostStarId({ hostname: 'Unknown Star XYZ', raDeg: 0, decDeg: 0, distancePc: 100 }, FIXTURE_STARS, 0.5);
|
||||
|
||||
expect(id).toBeNull();
|
||||
});
|
||||
|
||||
it('returns null when there is no name match and no position is available', () => {
|
||||
const id = resolveHostStarId({ hostname: 'Unknown Star XYZ', raDeg: NaN, decDeg: NaN, distancePc: NaN }, FIXTURE_STARS, 0.5);
|
||||
|
||||
expect(id).toBeNull();
|
||||
});
|
||||
|
||||
it('picks the closest star when more than one falls within tolerance', () => {
|
||||
const stars: StarRecord[] = [
|
||||
{ id: 10, name: 'Near', x: 0, y: 0, z: 0, magnitude: 5, spectralType: 'G', colorIndex: 0.5 },
|
||||
{ id: 11, name: 'Far', x: 0.4, y: 0, z: 0, magnitude: 5, spectralType: 'G', colorIndex: 0.5 }
|
||||
];
|
||||
const nameIndex = buildStarNameIndex(stars);
|
||||
|
||||
const id = resolveHostStarId({ hostname: 'Unmatched', raDeg: 0, decDeg: 0, distancePc: 0.2 }, stars, 0.5, nameIndex);
|
||||
|
||||
expect(id).toBe(10);
|
||||
});
|
||||
});
|
||||
@@ -0,0 +1,59 @@
|
||||
import { CartesianCoordinates, distanceBetween, raDegDecDistanceToXyz } from './coordinates';
|
||||
import { StarRecord } from '../models/star.model';
|
||||
|
||||
/** Normalizes a star name for comparison: lowercase, alphanumeric characters only. */
|
||||
export function normalizeStarName(name: string): string {
|
||||
return name.toLowerCase().replace(/[^a-z0-9]/g, '');
|
||||
}
|
||||
|
||||
export interface HostStarQuery {
|
||||
hostname: string;
|
||||
raDeg: number;
|
||||
decDeg: number;
|
||||
distancePc: number;
|
||||
}
|
||||
|
||||
/** Builds a lookup of normalized star name -> star, for fast repeated name matching. */
|
||||
export function buildStarNameIndex(stars: readonly StarRecord[]): Map<string, StarRecord> {
|
||||
return new Map(stars.map((star) => [normalizeStarName(star.name), star]));
|
||||
}
|
||||
|
||||
/**
|
||||
* Cross-references an exoplanet host star to the HYG star index: first by (normalized)
|
||||
* name, then by nearest-neighbour position matching within `toleranceInPc`. Returns `null`
|
||||
* when neither approach finds a confident match, rather than guessing.
|
||||
*
|
||||
* `nameIndex` should be built once (via {@link buildStarNameIndex}) and reused across calls
|
||||
* when resolving many queries against the same star list.
|
||||
*/
|
||||
export function resolveHostStarId(
|
||||
query: HostStarQuery,
|
||||
stars: readonly StarRecord[],
|
||||
toleranceInPc: number,
|
||||
nameIndex: Map<string, StarRecord> = buildStarNameIndex(stars)
|
||||
): number | null {
|
||||
const byName = nameIndex.get(normalizeStarName(query.hostname));
|
||||
if (byName) {
|
||||
return byName.id;
|
||||
}
|
||||
|
||||
if (![query.raDeg, query.decDeg, query.distancePc].every(Number.isFinite)) {
|
||||
return null;
|
||||
}
|
||||
|
||||
const hostPosition = raDegDecDistanceToXyz(query.raDeg, query.decDeg, query.distancePc);
|
||||
return findNearestStarWithin(hostPosition, stars, toleranceInPc);
|
||||
}
|
||||
|
||||
function findNearestStarWithin(position: CartesianCoordinates, stars: readonly StarRecord[], toleranceInPc: number): number | null {
|
||||
let closest: { id: number; distance: number } | null = null;
|
||||
|
||||
for (const star of stars) {
|
||||
const distance = distanceBetween(position, star);
|
||||
if (distance <= toleranceInPc && (!closest || distance < closest.distance)) {
|
||||
closest = { id: star.id, distance };
|
||||
}
|
||||
}
|
||||
|
||||
return closest ? closest.id : null;
|
||||
}
|
||||
@@ -0,0 +1,165 @@
|
||||
import { describe, expect, it } from 'vitest';
|
||||
|
||||
import { GM_SUN_AU3_PER_DAY2, DEFAULT_EPOCH_JD } from './constants';
|
||||
import {
|
||||
meanMotionRadPerDay,
|
||||
orbitEllipsePoints,
|
||||
orbitalPeriodDays,
|
||||
positionAtTrueAnomaly,
|
||||
propagateOrbit,
|
||||
resolveOrbitalElements,
|
||||
solveEccentricAnomaly,
|
||||
trueAnomalyFromEccentricAnomaly
|
||||
} from './kepler';
|
||||
|
||||
// Earth's actual orbital elements (osculating, ~J2000), used as a real-world reference case.
|
||||
const EARTH_ELEMENTS = {
|
||||
semiMajorAxisAu: 1.00000011,
|
||||
eccentricity: 0.01671022,
|
||||
inclinationDeg: 0.00005,
|
||||
longitudeOfAscendingNodeDeg: -11.26064,
|
||||
argumentOfPeriapsisDeg: 102.94719,
|
||||
meanAnomalyAtEpochDeg: 100.46435,
|
||||
epochJd: DEFAULT_EPOCH_JD
|
||||
};
|
||||
|
||||
describe('solveEccentricAnomaly', () => {
|
||||
it('satisfies Keplers equation for a range of eccentricities', () => {
|
||||
for (const eccentricity of [0, 0.0167, 0.3, 0.6, 0.9]) {
|
||||
for (const meanAnomalyRad of [0, 0.5, 1.5, 3.0, 5.5]) {
|
||||
const e = solveEccentricAnomaly(meanAnomalyRad, eccentricity);
|
||||
const residual = e - eccentricity * Math.sin(e) - meanAnomalyRad;
|
||||
// residual is computed against the (possibly un-normalized) input, but the solver
|
||||
// normalizes internally, so compare against the normalized mean anomaly instead.
|
||||
const normalizedMeanAnomaly = ((meanAnomalyRad % (2 * Math.PI)) + 2 * Math.PI) % (2 * Math.PI);
|
||||
expect(e - eccentricity * Math.sin(e)).toBeCloseTo(normalizedMeanAnomaly, 6);
|
||||
expect(Number.isFinite(residual)).toBe(true);
|
||||
}
|
||||
}
|
||||
});
|
||||
});
|
||||
|
||||
describe('trueAnomalyFromEccentricAnomaly', () => {
|
||||
it('returns 0 at periapsis and pi at apoapsis', () => {
|
||||
expect(trueAnomalyFromEccentricAnomaly(0, 0.3)).toBeCloseTo(0, 9);
|
||||
expect(trueAnomalyFromEccentricAnomaly(Math.PI, 0.3)).toBeCloseTo(Math.PI, 9);
|
||||
});
|
||||
|
||||
it('matches the eccentric anomaly exactly for a circular orbit', () => {
|
||||
expect(trueAnomalyFromEccentricAnomaly(1.234, 0)).toBeCloseTo(1.234, 9);
|
||||
});
|
||||
});
|
||||
|
||||
describe('positionAtTrueAnomaly', () => {
|
||||
it('places a circular, unrotated orbit at radius = semiMajorAxisAu for every true anomaly', () => {
|
||||
const circular = resolveOrbitalElements({ semiMajorAxisAu: 2.5, eccentricity: 0 });
|
||||
|
||||
for (const trueAnomalyRad of [0, Math.PI / 2, Math.PI, (3 * Math.PI) / 2]) {
|
||||
const { x, y, z } = positionAtTrueAnomaly(circular, trueAnomalyRad);
|
||||
expect(Math.sqrt(x * x + y * y + z * z)).toBeCloseTo(2.5, 9);
|
||||
}
|
||||
});
|
||||
|
||||
it('reaches periapsis distance a*(1-e) and apoapsis distance a*(1+e)', () => {
|
||||
const elements = resolveOrbitalElements({ semiMajorAxisAu: 10, eccentricity: 0.2 });
|
||||
|
||||
const periapsis = positionAtTrueAnomaly(elements, 0);
|
||||
const apoapsis = positionAtTrueAnomaly(elements, Math.PI);
|
||||
|
||||
expect(Math.hypot(periapsis.x, periapsis.y, periapsis.z)).toBeCloseTo(8, 9);
|
||||
expect(Math.hypot(apoapsis.x, apoapsis.y, apoapsis.z)).toBeCloseTo(12, 9);
|
||||
});
|
||||
|
||||
it('tilts a 90-degree-inclined orbit entirely onto the z axis at true anomaly 90 degrees', () => {
|
||||
const elements = resolveOrbitalElements({ semiMajorAxisAu: 1, eccentricity: 0, inclinationDeg: 90 });
|
||||
|
||||
const { x, y, z } = positionAtTrueAnomaly(elements, Math.PI / 2);
|
||||
|
||||
expect(x).toBeCloseTo(0, 9);
|
||||
expect(y).toBeCloseTo(0, 9);
|
||||
expect(z).toBeCloseTo(1, 9);
|
||||
});
|
||||
});
|
||||
|
||||
describe('meanMotionRadPerDay / orbitalPeriodDays', () => {
|
||||
it('reproduces Earths ~365.25-day year from its semi-major axis', () => {
|
||||
const period = orbitalPeriodDays(EARTH_ELEMENTS.semiMajorAxisAu, GM_SUN_AU3_PER_DAY2);
|
||||
expect(period).toBeCloseTo(365.25, 0);
|
||||
});
|
||||
|
||||
it('is the inverse of orbitalPeriodDays', () => {
|
||||
const n = meanMotionRadPerDay(1, GM_SUN_AU3_PER_DAY2);
|
||||
const period = orbitalPeriodDays(1, GM_SUN_AU3_PER_DAY2);
|
||||
expect(n * period).toBeCloseTo(2 * Math.PI, 9);
|
||||
});
|
||||
});
|
||||
|
||||
describe('propagateOrbit', () => {
|
||||
it('reduces to the instantaneous position at the elements own epoch', () => {
|
||||
const eccentricAnomalyRad = solveEccentricAnomaly((EARTH_ELEMENTS.meanAnomalyAtEpochDeg * Math.PI) / 180, EARTH_ELEMENTS.eccentricity);
|
||||
const trueAnomalyRad = trueAnomalyFromEccentricAnomaly(eccentricAnomalyRad, EARTH_ELEMENTS.eccentricity);
|
||||
const expected = positionAtTrueAnomaly(EARTH_ELEMENTS, trueAnomalyRad);
|
||||
|
||||
const actual = propagateOrbit(EARTH_ELEMENTS, GM_SUN_AU3_PER_DAY2, EARTH_ELEMENTS.epochJd);
|
||||
|
||||
expect(actual.x).toBeCloseTo(expected.x, 9);
|
||||
expect(actual.y).toBeCloseTo(expected.y, 9);
|
||||
expect(actual.z).toBeCloseTo(expected.z, 9);
|
||||
});
|
||||
|
||||
it('stays within the periapsis/apoapsis distance bounds after propagating forward a year', () => {
|
||||
const period = orbitalPeriodDays(EARTH_ELEMENTS.semiMajorAxisAu, GM_SUN_AU3_PER_DAY2);
|
||||
const { x, y, z } = propagateOrbit(EARTH_ELEMENTS, GM_SUN_AU3_PER_DAY2, EARTH_ELEMENTS.epochJd + period * 0.37);
|
||||
const distance = Math.hypot(x, y, z);
|
||||
|
||||
const { semiMajorAxisAu: a, eccentricity: e } = EARTH_ELEMENTS;
|
||||
expect(distance).toBeGreaterThanOrEqual(a * (1 - e) - 1e-6);
|
||||
expect(distance).toBeLessThanOrEqual(a * (1 + e) + 1e-6);
|
||||
});
|
||||
|
||||
it('returns to (very nearly) the same position after exactly one full orbital period', () => {
|
||||
const period = orbitalPeriodDays(EARTH_ELEMENTS.semiMajorAxisAu, GM_SUN_AU3_PER_DAY2);
|
||||
const start = propagateOrbit(EARTH_ELEMENTS, GM_SUN_AU3_PER_DAY2, EARTH_ELEMENTS.epochJd + 12.3);
|
||||
const afterOneOrbit = propagateOrbit(EARTH_ELEMENTS, GM_SUN_AU3_PER_DAY2, EARTH_ELEMENTS.epochJd + 12.3 + period);
|
||||
|
||||
expect(afterOneOrbit.x).toBeCloseTo(start.x, 6);
|
||||
expect(afterOneOrbit.y).toBeCloseTo(start.y, 6);
|
||||
expect(afterOneOrbit.z).toBeCloseTo(start.z, 6);
|
||||
});
|
||||
});
|
||||
|
||||
describe('orbitEllipsePoints', () => {
|
||||
it('samples a closed loop whose distances stay within the periapsis/apoapsis bounds', () => {
|
||||
const elements = resolveOrbitalElements({ semiMajorAxisAu: 5, eccentricity: 0.4 });
|
||||
const points = orbitEllipsePoints(elements, 64);
|
||||
|
||||
expect(points).toHaveLength(65);
|
||||
for (const { x, y, z } of points) {
|
||||
const distance = Math.hypot(x, y, z);
|
||||
expect(distance).toBeGreaterThanOrEqual(5 * (1 - 0.4) - 1e-9);
|
||||
expect(distance).toBeLessThanOrEqual(5 * (1 + 0.4) + 1e-9);
|
||||
}
|
||||
|
||||
// The first and last sampled points (true anomaly 0 and 2*pi) should coincide.
|
||||
expect(points[0].x).toBeCloseTo(points[64].x, 9);
|
||||
expect(points[0].y).toBeCloseTo(points[64].y, 9);
|
||||
expect(points[0].z).toBeCloseTo(points[64].z, 9);
|
||||
});
|
||||
});
|
||||
|
||||
describe('resolveOrbitalElements', () => {
|
||||
it('defaults missing angles to 0 and the missing epoch to J2000', () => {
|
||||
const resolved = resolveOrbitalElements({ semiMajorAxisAu: 1.5, eccentricity: 0.1 });
|
||||
|
||||
expect(resolved.inclinationDeg).toBe(0);
|
||||
expect(resolved.longitudeOfAscendingNodeDeg).toBe(0);
|
||||
expect(resolved.argumentOfPeriapsisDeg).toBe(0);
|
||||
expect(resolved.meanAnomalyAtEpochDeg).toBe(0);
|
||||
expect(resolved.epochJd).toBe(DEFAULT_EPOCH_JD);
|
||||
});
|
||||
|
||||
it('preserves explicitly provided fields', () => {
|
||||
const resolved = resolveOrbitalElements({ semiMajorAxisAu: 1.5, eccentricity: 0.1, argumentOfPeriapsisDeg: 50 });
|
||||
expect(resolved.argumentOfPeriapsisDeg).toBe(50);
|
||||
});
|
||||
});
|
||||
@@ -0,0 +1,136 @@
|
||||
import { CartesianCoordinates } from './coordinates';
|
||||
import { DEFAULT_EPOCH_JD } from './constants';
|
||||
import { OrbitalElements } from '../models/body.model';
|
||||
|
||||
const DEG_TO_RAD = Math.PI / 180;
|
||||
const TWO_PI = Math.PI * 2;
|
||||
|
||||
/**
|
||||
* Fills in the elements the Kepler propagator needs but that some sources (e.g. exoplanets,
|
||||
* see `ExoplanetRecord.orbit: Partial<OrbitalElements>`) don't report: inclination, longitude
|
||||
* of ascending node, mean anomaly at epoch, and the epoch itself. Missing angles default to
|
||||
* zero (a face-on, unrotated ellipse) and the missing epoch defaults to J2000 — enough to draw
|
||||
* a plausible, period-correct orbit even without full data.
|
||||
*/
|
||||
export function resolveOrbitalElements(partial: Partial<OrbitalElements> & Pick<OrbitalElements, 'semiMajorAxisAu' | 'eccentricity'>): OrbitalElements {
|
||||
return {
|
||||
semiMajorAxisAu: partial.semiMajorAxisAu,
|
||||
eccentricity: partial.eccentricity,
|
||||
inclinationDeg: partial.inclinationDeg ?? 0,
|
||||
longitudeOfAscendingNodeDeg: partial.longitudeOfAscendingNodeDeg ?? 0,
|
||||
argumentOfPeriapsisDeg: partial.argumentOfPeriapsisDeg ?? 0,
|
||||
meanAnomalyAtEpochDeg: partial.meanAnomalyAtEpochDeg ?? 0,
|
||||
epochJd: partial.epochJd ?? DEFAULT_EPOCH_JD
|
||||
};
|
||||
}
|
||||
|
||||
/** Mean motion (rad/day) of a body via Kepler's third law: n = sqrt(GM / a^3). */
|
||||
export function meanMotionRadPerDay(semiMajorAxisAu: number, gmAu3PerDay2: number): number {
|
||||
return Math.sqrt(gmAu3PerDay2 / (semiMajorAxisAu * semiMajorAxisAu * semiMajorAxisAu));
|
||||
}
|
||||
|
||||
/** Orbital period (days) of a body via Kepler's third law: T = 2*pi / n. */
|
||||
export function orbitalPeriodDays(semiMajorAxisAu: number, gmAu3PerDay2: number): number {
|
||||
return TWO_PI / meanMotionRadPerDay(semiMajorAxisAu, gmAu3PerDay2);
|
||||
}
|
||||
|
||||
/** Normalizes an angle (radians) into [0, 2*pi). */
|
||||
function normalizeAngle(angleRad: number): number {
|
||||
const wrapped = angleRad % TWO_PI;
|
||||
return wrapped < 0 ? wrapped + TWO_PI : wrapped;
|
||||
}
|
||||
|
||||
/**
|
||||
* Solves Kepler's equation `M = E - e*sin(E)` for the eccentric anomaly `E` (radians) via
|
||||
* Newton-Raphson iteration.
|
||||
*/
|
||||
export function solveEccentricAnomaly(meanAnomalyRad: number, eccentricity: number, tolerance = 1e-8, maxIterations = 30): number {
|
||||
const m = normalizeAngle(meanAnomalyRad);
|
||||
let e = eccentricity < 0.8 ? m : Math.PI;
|
||||
|
||||
for (let i = 0; i < maxIterations; i++) {
|
||||
const delta = (e - eccentricity * Math.sin(e) - m) / (1 - eccentricity * Math.cos(e));
|
||||
e -= delta;
|
||||
if (Math.abs(delta) < tolerance) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
return e;
|
||||
}
|
||||
|
||||
/** Converts an eccentric anomaly (radians) into the true anomaly (radians). */
|
||||
export function trueAnomalyFromEccentricAnomaly(eccentricAnomalyRad: number, eccentricity: number): number {
|
||||
const cosE = Math.cos(eccentricAnomalyRad);
|
||||
const sinE = Math.sin(eccentricAnomalyRad);
|
||||
return Math.atan2(Math.sqrt(1 - eccentricity * eccentricity) * sinE, cosE - eccentricity);
|
||||
}
|
||||
|
||||
/**
|
||||
* Places a point at the given true anomaly (radians) along the orbit described by
|
||||
* `elements`, in AU, relative to the central body (the Sun for planets/dwarfs, the host
|
||||
* planet for moons — see `BodyRecord.parentBodyId`). Standard perifocal-to-reference-frame
|
||||
* rotation: argument of periapsis, then inclination, then longitude of ascending node.
|
||||
*/
|
||||
export function positionAtTrueAnomaly(elements: OrbitalElements, trueAnomalyRad: number): CartesianCoordinates {
|
||||
const { semiMajorAxisAu: a, eccentricity: e } = elements;
|
||||
const semiLatusRectum = a * (1 - e * e);
|
||||
const radius = semiLatusRectum / (1 + e * Math.cos(trueAnomalyRad));
|
||||
|
||||
// Position in the perifocal (orbital-plane) frame: +x toward periapsis.
|
||||
const xPerifocal = radius * Math.cos(trueAnomalyRad);
|
||||
const yPerifocal = radius * Math.sin(trueAnomalyRad);
|
||||
|
||||
const omega = elements.argumentOfPeriapsisDeg * DEG_TO_RAD; // argument of periapsis
|
||||
const inclination = elements.inclinationDeg * DEG_TO_RAD;
|
||||
const raan = elements.longitudeOfAscendingNodeDeg * DEG_TO_RAD; // right ascension of ascending node
|
||||
|
||||
const cosOmega = Math.cos(omega);
|
||||
const sinOmega = Math.sin(omega);
|
||||
const cosInclination = Math.cos(inclination);
|
||||
const sinInclination = Math.sin(inclination);
|
||||
const cosRaan = Math.cos(raan);
|
||||
const sinRaan = Math.sin(raan);
|
||||
|
||||
// Rotate by argument of periapsis within the orbital plane first.
|
||||
const xOrbitPlane = xPerifocal * cosOmega - yPerifocal * sinOmega;
|
||||
const yOrbitPlane = xPerifocal * sinOmega + yPerifocal * cosOmega;
|
||||
|
||||
// Tilt by inclination, then rotate by the longitude of the ascending node.
|
||||
const xTilted = xOrbitPlane;
|
||||
const yTilted = yOrbitPlane * cosInclination;
|
||||
const zTilted = yOrbitPlane * sinInclination;
|
||||
|
||||
return {
|
||||
x: xTilted * cosRaan - yTilted * sinRaan,
|
||||
y: xTilted * sinRaan + yTilted * cosRaan,
|
||||
z: zTilted
|
||||
};
|
||||
}
|
||||
|
||||
/**
|
||||
* Propagates `elements` to Julian date `epochJdEval`, returning the body's position (AU)
|
||||
* relative to its central body. This is the app's "current epoch" evaluation used for live
|
||||
* (and future time-scrubbable) positions, as opposed to {@link orbitEllipsePoints} which
|
||||
* samples the fixed orbit shape independent of time.
|
||||
*/
|
||||
export function propagateOrbit(elements: OrbitalElements, gmAu3PerDay2: number, epochJdEval: number): CartesianCoordinates {
|
||||
const meanMotion = meanMotionRadPerDay(elements.semiMajorAxisAu, gmAu3PerDay2);
|
||||
const meanAnomalyRad = elements.meanAnomalyAtEpochDeg * DEG_TO_RAD + meanMotion * (epochJdEval - elements.epochJd);
|
||||
const eccentricAnomalyRad = solveEccentricAnomaly(meanAnomalyRad, elements.eccentricity);
|
||||
const trueAnomalyRad = trueAnomalyFromEccentricAnomaly(eccentricAnomalyRad, elements.eccentricity);
|
||||
return positionAtTrueAnomaly(elements, trueAnomalyRad);
|
||||
}
|
||||
|
||||
/**
|
||||
* Samples `segments` points around the fixed shape of the orbit (AU, relative to the central
|
||||
* body), for drawing the orbit ellipse. Independent of epoch/time — unlike {@link propagateOrbit}.
|
||||
*/
|
||||
export function orbitEllipsePoints(elements: OrbitalElements, segments = 128): CartesianCoordinates[] {
|
||||
const points: CartesianCoordinates[] = [];
|
||||
for (let i = 0; i <= segments; i++) {
|
||||
const trueAnomalyRad = (i / segments) * TWO_PI;
|
||||
points.push(positionAtTrueAnomaly(elements, trueAnomalyRad));
|
||||
}
|
||||
return points;
|
||||
}
|
||||
Reference in New Issue
Block a user