solargod/js/ephem.js

// SOLARGOD ephemeris engine — JPL SSD "Approximate Positions of the Planets"
// (Standish). Tables extracted 2026-07-15 from
// https://ssd.jpl.nasa.gov/planets/approx_pos.html — authoritative, verbatim.
// Angles in degrees, rates per Julian century; `e` is dimensionless (the page's
// "rad" header for e is its own typo). No n-body integration — this IS the spec.

import { centuriesSinceJ2000, normDegPM180, DEG } from './lib.js';

// Row layout: el = [a(au), e, I(deg), L(deg), ϖ long.peri(deg), Ω long.node(deg)]
//             rate = same order, per century.
const TABLE1 = {
  mercury: { el: [0.38709927, 0.20563593, 7.00497902, 252.25032350, 77.45779628, 48.33076593],
             rate: [0.00000037, 0.00001906, -0.00594749, 149472.67411175, 0.16047689, -0.12534081] },
  venus:   { el: [0.72333566, 0.00677672, 3.39467605, 181.97909950, 131.60246718, 76.67984255],
             rate: [0.00000390, -0.00004107, -0.00078890, 58517.81538729, 0.00268329, -0.27769418] },
  earth:   { el: [1.00000261, 0.01671123, -0.00001531, 100.46457166, 102.93768193, 0.0],
             rate: [0.00000562, -0.00004392, -0.01294668, 35999.37244981, 0.32327364, 0.0] },
  mars:    { el: [1.52371034, 0.09339410, 1.84969142, -4.55343205, -23.94362959, 49.55953891],
             rate: [0.00001847, 0.00007882, -0.00813131, 19140.30268499, 0.44441088, -0.29257343] },
  jupiter: { el: [5.20288700, 0.04838624, 1.30439695, 34.39644051, 14.72847983, 100.47390909],
             rate: [-0.00011607, -0.00013253, -0.00183714, 3034.74612775, 0.21252668, 0.20469106] },
  saturn:  { el: [9.53667594, 0.05386179, 2.48599187, 49.95424423, 92.59887831, 113.66242448],
             rate: [-0.00125060, -0.00050991, 0.00193609, 1222.49362201, -0.41897216, -0.28867794] },
  uranus:  { el: [19.18916464, 0.04725744, 0.77263783, 313.23810451, 170.95427630, 74.01692503],
             rate: [-0.00196176, -0.00004397, -0.00242939, 428.48202785, 0.40805281, 0.04240589] },
  neptune: { el: [30.06992276, 0.00859048, 1.77004347, -55.12002969, 44.96476227, 131.78422574],
             rate: [0.00026291, 0.00005105, 0.00035372, 218.45945325, -0.32241464, -0.00508664] },
};

// Table 2a — valid 3000 BC–3000 AD. Shipped for the extended-range config swap
// (CONFIG.useExtendedRange); unused by v1 UI.
const TABLE2A = {
  mercury: { el: [0.38709843, 0.20563661, 7.00559432, 252.25166724, 77.45771895, 48.33961819],
             rate: [0.00000000, 0.00002123, -0.00590158, 149472.67486623, 0.15940013, -0.12214182] },
  venus:   { el: [0.72332102, 0.00676399, 3.39777545, 181.97970850, 131.76755713, 76.67261496],
             rate: [-0.00000026, -0.00005107, 0.00043494, 58517.81560260, 0.05679648, -0.27274174] },
  earth:   { el: [1.00000018, 0.01673163, -0.00054346, 100.46691572, 102.93005885, -5.11260389],
             rate: [-0.00000003, -0.00003661, -0.01337178, 35999.37306329, 0.31795260, -0.24123856] },
  mars:    { el: [1.52371243, 0.09336511, 1.85181869, -4.56813164, -23.91744784, 49.71320984],
             rate: [0.00000097, 0.00009149, -0.00724757, 19140.29934243, 0.45223625, -0.26852431] },
  jupiter: { el: [5.20248019, 0.04853590, 1.29861416, 34.33479152, 14.27495244, 100.29282654],
             rate: [-0.00002864, 0.00018026, -0.00322699, 3034.90371757, 0.18199196, 0.13024619] },
  saturn:  { el: [9.54149883, 0.05550825, 2.49424102, 50.07571329, 92.86136063, 113.63998702],
             rate: [-0.00003065, -0.00032044, 0.00451969, 1222.11494724, 0.54179478, -0.25015002] },
  uranus:  { el: [19.18797948, 0.04685740, 0.77298127, 314.20276625, 172.43404441, 73.96250215],
             rate: [-0.00020455, -0.00001550, -0.00180155, 428.49512595, 0.09266985, 0.05739699] },
  neptune: { el: [30.06952752, 0.00895439, 1.77005520, 304.22289287, 46.68158724, 131.78635853],
             rate: [0.00006447, 0.00000818, 0.00022400, 218.46515314, 0.01009938, -0.00606302] },
};

// Table 2b — extra terms added to M (deg) for Jupiter–Neptune, Table 2a only:
// M += b·T² + c·cos(f·T) + s·sin(f·T).  [b, c, s, f]
const TABLE2B = {
  jupiter: [-0.00012452, 0.06064060, -0.35635438, 38.35125000],
  saturn:  [0.00025899, -0.13434469, 0.87320147, 38.35125000],
  uranus:  [0.00058331, -0.97731848, 0.17689245, 7.67025000],
  neptune: [-0.00041348, 0.68346318, -0.10162547, 7.67025000],
};

export const PLANET_IDS = ['mercury', 'venus', 'earth', 'mars', 'jupiter', 'saturn', 'uranus', 'neptune'];

let useExtended = false;
export function setExtendedRange(on) { useExtended = !!on; }

// Solve Kepler's equation E − e·sinE = M (radians) by Newton–Raphson.
// Seed E₀ = M + e·sinM; converges in a few iters even for e→0.97 (verify gate).
export function solveKepler(M, e, tol = 1e-9, maxIter = 20) {
  let E = M + e * Math.sin(M);
  for (let i = 0; i < maxIter; i++) {
    const dE = (E - e * Math.sin(E) - M) / (1 - e * Math.cos(E));
    E -= dE;
    if (Math.abs(dE) < tol) break;
  }
  return E;
}

// The reusable core: resolved classical elements → heliocentric J2000 ecliptic
// {x,y,z} in AU. a[au], e, and iDeg/OmDeg(Ω)/wDeg(ω arg.peri)/Mdeg all in degrees.
// Shared by planets (element+rate) and small bodies (epoch elements) — brief §3.
export function eclFromElements(a, e, iDeg, OmDeg, wDeg, Mdeg, out = {}) {
  const M = normDegPM180(Mdeg) * DEG;
  const E = solveKepler(M, e);
  // Orbital-plane coords (perifocal), AU.
  const xp = a * (Math.cos(E) - e);
  const yp = a * Math.sqrt(1 - e * e) * Math.sin(E);
  const w = wDeg * DEG, Om = OmDeg * DEG, I = iDeg * DEG;
  const cw = Math.cos(w), sw = Math.sin(w);
  const cO = Math.cos(Om), sO = Math.sin(Om);
  const cI = Math.cos(I), sI = Math.sin(I);
  // Perifocal → ecliptic rotation (brief §3 step 5 = standard R_z(Ω)R_x(I)R_z(ω)).
  out.x = (cw * cO - sw * sO * cI) * xp + (-sw * cO - cw * sO * cI) * yp;
  out.y = (cw * sO + sw * cO * cI) * xp + (-sw * sO + cw * cO * cI) * yp;
  out.z = (sw * sI) * xp + (cw * sI) * yp;
  return out;
}

// Resolve a planet's elements at jd (Table 1, 1800–2050) → {a,e,I,Om,w,M} deg.
function resolvePlanet(bodyId, jd) {
  const T = centuriesSinceJ2000(jd);
  const row = TABLE1[bodyId];
  if (!row) throw new Error(`ephem: unknown planet ${bodyId}`);
  const a  = row.el[0] + row.rate[0] * T;
  const e  = row.el[1] + row.rate[1] * T;
  const I  = row.el[2] + row.rate[2] * T;
  const L  = row.el[3] + row.rate[3] * T;
  const wb = row.el[4] + row.rate[4] * T; // ϖ  longitude of perihelion
  const Om = row.el[5] + row.rate[5] * T; // Ω  longitude of ascending node
  const w = wb - Om;                       // ω  argument of perihelion
  const M = L - wb;                        // mean anomaly (normalize AFTER subtract)
  return { a, e, I, Om, w, M };
}
// Extended-range resolver (Table 2a + the 2b M-correction for Jupiter–Neptune,
// valid 3000 BC–3000 AD). f·T uses T in centuries; f is deg-per-century so f·T
// is degrees → convert to rad for the trig. Only used when useExtended is set.
function resolvePlanetExt(bodyId, jd) {
  const T = centuriesSinceJ2000(jd);
  const row = TABLE2A[bodyId];
  const a = row.el[0] + row.rate[0] * T;
  const e = row.el[1] + row.rate[1] * T;
  const I = row.el[2] + row.rate[2] * T;
  const L = row.el[3] + row.rate[3] * T;
  const wb = row.el[4] + row.rate[4] * T;
  const Om = row.el[5] + row.rate[5] * T;
  const w = wb - Om;
  let M = L - wb;
  const b2 = TABLE2B[bodyId];
  if (b2) {
    const [b, c, s, f] = b2;
    M += b * T * T + c * Math.cos(f * T * DEG) + s * Math.sin(f * T * DEG);
  }
  return { a, e, I, Om, w, M };
}

// Public: heliocentric J2000 ecliptic position of a planet at jd, AU.
export function helioEcl(bodyId, jd, out = {}) {
  const p = useExtended ? resolvePlanetExt(bodyId, jd) : resolvePlanet(bodyId, jd);
  return eclFromElements(p.a, p.e, p.I, p.Om, p.w, p.M, out);
}

// Mean motion [deg/day] from semi-major axis via Kepler's third law (Gaussian).
export function meanMotionDegPerDay(a) { return 0.9856076686 / Math.pow(a, 1.5); }

// Resolved semi-major axis [AU] of a planet at jd (for vis-viva speed, year length).
export function semiMajor(bodyId, jd) {
  const p = useExtended ? resolvePlanetExt(bodyId, jd) : resolvePlanet(bodyId, jd);
  return p.a;
}

// Small-body position (asteroids/comets/Pluto): epoch elements propagated by
// M = M0 + n·(jd − epoch). el = {a,e,i,om(Ω),w(ω),ma(M0),epoch}.
export function smallBodyEcl(el, jd, out = {}) {
  const n = meanMotionDegPerDay(el.a);
  const M = el.ma + n * (jd - el.epoch);
  return eclFromElements(el.a, el.e, el.i, el.om, el.w, M, out);
}

// One full orbit sampled in mean anomaly, as a Float64Array of xyz triples
// (nSamples+1 points; last == first so a LineLoop/LineSegments closes cleanly —
// verify gate "orbitPath endpoints join"). Elements frozen at jd (they drift
// slowly; caller rebuilds only when |ΔT| > 0.1 Cy, brief §11).
export function orbitPath(bodyId, jd, nSamples = 256) {
  const p = useExtended ? resolvePlanetExt(bodyId, jd) : resolvePlanet(bodyId, jd);
  return samplePath(p.a, p.e, p.I, p.Om, p.w, nSamples);
}
export function orbitPathFromElements(el, nSamples = 256) {
  return samplePath(el.a, el.e, el.i, el.om, el.w, nSamples);
}
function samplePath(a, e, I, Om, w, nSamples) {
  const out = new Float64Array((nSamples + 1) * 3);
  const tmp = {};
  for (let k = 0; k <= nSamples; k++) {
    const M = (360 * k) / nSamples; // 0..360; k=nSamples reproduces k=0
    eclFromElements(a, e, I, Om, w, M, tmp);
    out[k * 3] = tmp.x; out[k * 3 + 1] = tmp.y; out[k * 3 + 2] = tmp.z;
  }
  return out;
}