src/physics/formulas.js

// ————————————————————————————————————————————————
// Physics formulas used across the views.
// ————————————————————————————————————————————————

import {
  C_KMS,
  KM_PER_AU,
  KM_PER_LY,
  SECONDS_PER_YEAR,
  EARTH_SIM_PERIOD_SEC,
} from '../data/constants.js';

// ————————————————————————————————————————————————
// Orbital mechanics
// ————————————————————————————————————————————————

/**
 * Kepler's third law (normalised to Earth).
 *
 * For a planet with orbital period `periodDays`, returns its orbital angular
 * velocity in the simulation's "sim seconds" — i.e. how many radians it
 * sweeps per second of real animation time at 1× speed.
 *
 *   ω_sim = 2π / (EARTH_SIM_PERIOD_SEC × periodYears)
 *
 * This keeps the RELATIVE speed between planets exactly correct: Mercury
 * orbits 4.15× faster than Earth because its real period is 87.97 vs 365.25
 * days. Only the absolute time is rescaled to make the animation watchable.
 */
export function angularVelocityAtOneX(periodDays) {
  const periodYears = periodDays / 365.25;
  return (2 * Math.PI) / (EARTH_SIM_PERIOD_SEC * periodYears);
}

/**
 * Keplerian radius on an elliptical orbit as a function of true anomaly θ.
 *
 *   r(θ) = a(1 - e²) / (1 + e·cos θ)
 *
 * Used by the "Planet X" view (now removed) to animate an eccentric orbit.
 */
export function keplerRadius(a, e, theta) {
  return (a * (1 - e * e)) / (1 + e * Math.cos(theta));
}

// ————————————————————————————————————————————————
// Galactic rotation
// ————————————————————————————————————————————————

/**
 * Flat rotation curve approximation (valid for the Milky Way outside the bulge).
 *
 * The Milky Way's circular velocity v_c is approximately CONSTANT with radius,
 * not falling off with 1/√r as Kepler would predict (this is the key evidence
 * for dark matter). So:
 *
 *   ω(r) = v_c / r
 *
 * Inner stars whip around faster (angular), outer stars lag. At the Sun's
 * radius (~26 000 ly), v_c ≈ 220 km/s gives a full orbit in ~225 million years.
 *
 * This is what the galaxy view's vertex shader computes per-star.
 */
export function galacticAngularVelocity(vCircular, radius) {
  return vCircular / radius;
}

// ————————————————————————————————————————————————
// Special relativity
// ————————————————————————————————————————————————

/**
 * Lorentz factor γ at velocity v (in units of c).
 *
 *   γ = 1 / √(1 - β²)     where β = v/c
 *
 * Used for relativistic-ship time dilation: a journey that takes T seconds
 * in Earth's frame is experienced as T/γ seconds by the ship's crew.
 */
export function lorentzFactor(betaC) {
  if (betaC >= 1) return Infinity;
  return 1 / Math.sqrt(1 - betaC * betaC);
}

/**
 * Travel time at a given velocity (km/s) over a given distance (km).
 * Returns both Earth-frame time and ship proper time, in years.
 */
export function relativisticTravel(distanceKm, velocityKms) {
  const earthSec = distanceKm / velocityKms;
  const beta = velocityKms / C_KMS;
  const gamma = lorentzFactor(beta);
  return {
    earthYears: earthSec / SECONDS_PER_YEAR,
    shipYears: earthSec / SECONDS_PER_YEAR / gamma,
    gamma,
  };
}

// ————————————————————————————————————————————————
// Unit conversions
// ————————————————————————————————————————————————

// ————————————————————————————————————————————————
// Earth satellite orbital mechanics
// ————————————————————————————————————————————————

// Earth's standard gravitational parameter μ = GM, km³/s².
// From IAU 2015 nominal values: GM_E = 3.986 004 × 10⁵ km³/s².
const GM_EARTH_KM3_S2 = 398600.4418;

// Earth's equatorial radius, km (WGS-84).
export const EARTH_RADIUS_KM = 6378.137;

/**
 * Orbital period of a circular Earth satellite at a given altitude (km).
 *
 *   T = 2π √((R_E + h)³ / μ)
 *
 * Returns period in minutes.
 */
export function orbitalPeriodFromAltitude(altitudeKm) {
  const a = EARTH_RADIUS_KM + altitudeKm; // semi-major axis in km
  return (2 * Math.PI * Math.sqrt((a * a * a) / GM_EARTH_KM3_S2)) / 60;
}

/**
 * Angular velocity (rad / real second) for a satellite with the given
 * orbital period (minutes), normalised to the same sim-time convention as
 * angularVelocityAtOneX: at 1× speed Earth orbits the Sun in
 * EARTH_SIM_PERIOD_SEC real seconds.
 *
 * Satellites are fast — the ISS (period 92 min) laps Earth ~15.5 times per
 * day. At 1× sim speed that's extremely fast on screen; callers should
 * typically run the launch view at much higher sim speeds than the solar-
 * system views.
 *
 * At 1× speed: ω_sat = 2π / (periodMin × 60) rad/s_real
 * This is NOT normalised to EARTH_SIM_PERIOD_SEC — satellite orbits need
 * wall-clock seconds so their motion is physically correct relative to the
 * simulation's own clock (which the launch view advances in real seconds).
 */
export function satelliteAngularVelocity(periodMin) {
  return (2 * Math.PI) / (periodMin * 60); // rad / real_second at 1× speed
}

export const kmToAU = (km) => km / KM_PER_AU;
export const auToKm = (au) => au * KM_PER_AU;
export const kmToLY = (km) => km / KM_PER_LY;
export const lyToKm = (ly) => ly * KM_PER_LY;
export const lightTimeS = (km) => km / C_KMS;

// ————————————————————————————————————————————————
// Kepler orbit (planar, J2000.0 elements)
// ————————————————————————————————————————————————
//
// Phase 1a of #6: real eccentricity and orientation, but inclination is
// kept at zero so everything stays in the orbital plane (XZ in our scenes).
// Real-date anchoring and inclination are Phase 1b.

// Solve Keplers equation M = E - e·sin(E) for the eccentric anomaly E.
// Newton iteration; converges in 3–6 iterations for e < 0.5. Returns E in
// radians; M may be any real (not modulo-reduced).
export function solveKepler(meanAnomaly, eccentricity, { tol = 1e-12, maxIter = 12 } = {}) {
  let E = meanAnomaly;
  for (let i = 0; i < maxIter; i++) {
    const f = E - eccentricity * Math.sin(E) - meanAnomaly;
    const fp = 1 - eccentricity * Math.cos(E);
    const dE = f / fp;
    E -= dE;
    if (Math.abs(dE) < tol) break;
  }
  return E;
}

// Mean anomaly → true anomaly ν, via the eccentric anomaly. Returns ν in
// (-π, π]. e = 0 yields ν ≡ M (modulo 2π); identity-on-circles.
export function meanToTrueAnomaly(meanAnomaly, eccentricity) {
  const E = solveKepler(meanAnomaly, eccentricity);
  return (
    2 *
    Math.atan2(
      Math.sqrt(1 + eccentricity) * Math.sin(E / 2),
      Math.sqrt(1 - eccentricity) * Math.cos(E / 2)
    )
  );
}

// Heliocentric XZ position on a Kepler ellipse (i = 0). a in any length
// unit, M in radians. ascNode (Ω) and argPerihelion (ω) orient the orbit
// in the plane. Returns { x, z } in the same length unit as a.
export function keplerXZ(a, e, M, argPerihelion, ascNode) {
  const E = solveKepler(M, e);
  const r = a * (1 - e * Math.cos(E));
  const nu = 2 * Math.atan2(Math.sqrt(1 + e) * Math.sin(E / 2), Math.sqrt(1 - e) * Math.cos(E / 2));
  const ang = nu + argPerihelion + ascNode;
  return { x: r * Math.cos(ang), z: r * Math.sin(ang) };
}

// Sample N+1 points along the ellipse for orbit-ring rendering. Uniform in
// eccentric anomaly E gives smoother visual sampling than uniform M (which
// bunches near aphelion). Returns an array of { x, z } pairs, closing the
// loop (last point == first).
export function ellipsePoints(a, e, argPerihelion, ascNode, segments = 192) {
  const b = a * Math.sqrt(1 - e * e);
  const ang = argPerihelion + ascNode;
  const ca = Math.cos(ang),
    sa = Math.sin(ang);
  const pts = new Array(segments + 1);
  for (let i = 0; i <= segments; i++) {
    const E = (i / segments) * 2 * Math.PI;
    const xo = a * (Math.cos(E) - e); // ellipse focus at origin
    const yo = b * Math.sin(E);
    pts[i] = { x: xo * ca - yo * sa, z: xo * sa + yo * ca };
  }
  return pts;
}

// ————————————————————————————————————————————————
// Real-date anchor: Julian Date helpers and full 3D Kepler elements.
// Phase 1b of #6.
// ————————————————————————————————————————————————

// J2000.0 epoch in Julian Date.
export const J2000_JD = 2451545.0;

// Julian Date for "right now". 1970-01-01T00:00:00Z = JD 2440587.5.
export function jdNow() {
  return jdFromMs(Date.now());
}

// Julian Date for a Date / epoch milliseconds value.
export function jdFromMs(ms) {
  return ms / 86400000 + 2440587.5;
}

// Inverse of jdFromMs.
export function msFromJD(jd) {
  return (jd - 2440587.5) * 86400000;
}

// Linearly extrapolate the six Kepler elements from their J2000.0 values
// + per-Julian-century rates to the given JD. Returns a fresh object
// holding the six numbers; the input planet entry isn't mutated.
//
// The JPL "approximate-elements" fit is valid AD 1800–2050 to about an
// arcminute, multi-century to a few; rates are linear in time.
export function elementsAtJD(planet, jd) {
  const T = (jd - J2000_JD) / 36525; // Julian centuries since J2000
  return {
    a: planet.orbitKm, // a not stored as a rate
    e: planet.e + planet.eRate * T,
    iRad: planet.iRad + planet.iRadRate * T,
    ascNodeRad: planet.ascNodeRad + planet.ascNodeRadRate * T,
    argPerihelionRad: planet.argPerihelionRad + planet.argPerihelionRadRate * T,
    meanLongitudeRad: planet.meanLongitudeJ2000Rad + planet.meanLongitudeRadRate * T,
  };
}

// Mean anomaly at JD = L − ϖ where ϖ = Ω + ω. Wrapped to (-π, π] so
// solveKepler stays in its well-conditioned range.
export function meanAnomalyAtJD(planet, jd) {
  const el = elementsAtJD(planet, jd);
  let M = el.meanLongitudeRad - (el.ascNodeRad + el.argPerihelionRad);
  M = ((((M + Math.PI) % (2 * Math.PI)) + 2 * Math.PI) % (2 * Math.PI)) - Math.PI;
  return M;
}

// Heliocentric ecliptic position via the standard 3-2-3 rotation
// (perifocal → heliocentric ecliptic): rotate by ω in-plane, tilt by i,
// rotate by Ω. Returns { x, y, z } in the same length unit as `a`.
//
// Reduces to keplerXZ when iRad = 0, with y = 0.
export function keplerHeliocentric(a, e, M, argPerihelion, ascNode, iRad) {
  const E = solveKepler(M, e);
  const cE = Math.cos(E),
    sE = Math.sin(E);
  // Perifocal frame (focus at origin, +x toward perihelion).
  const xp = a * (cE - e);
  const yp = a * Math.sqrt(1 - e * e) * sE;
  // Rotate by ω about z (perifocal z = orbit normal).
  const co = Math.cos(argPerihelion),
    so = Math.sin(argPerihelion);
  const x1 = co * xp - so * yp;
  const y1 = so * xp + co * yp;
  // Tilt by i about x.
  const ci = Math.cos(iRad),
    si = Math.sin(iRad);
  const x2 = x1;
  const y2 = ci * y1;
  const z2 = si * y1;
  // Rotate by Ω about z (ecliptic normal).
  const cO = Math.cos(ascNode),
    sO = Math.sin(ascNode);
  return {
    x: cO * x2 - sO * y2,
    y: z2, // we use Y-up; orbital z → world y
    z: sO * x2 + cO * y2,
  };
}

// Sample N+1 points along the inclined ellipse for orbit-ring rendering.
// Uniform in eccentric anomaly E. Returns { x, y, z } in the same length
// unit as a, with the orbit's plane tilted by iRad about the line of
// nodes (ascNode).
export function inclinedEllipsePoints(a, e, argPerihelion, ascNode, iRad, segments = 192) {
  const b = a * Math.sqrt(1 - e * e);
  const co = Math.cos(argPerihelion),
    so = Math.sin(argPerihelion);
  const ci = Math.cos(iRad),
    si = Math.sin(iRad);
  const cO = Math.cos(ascNode),
    sO = Math.sin(ascNode);
  const pts = new Array(segments + 1);
  for (let i = 0; i <= segments; i++) {
    const E = (i / segments) * 2 * Math.PI;
    const xp = a * (Math.cos(E) - e);
    const yp = b * Math.sin(E);
    const x1 = co * xp - so * yp;
    const y1 = so * xp + co * yp;
    const y2 = ci * y1;
    const z2 = si * y1;
    pts[i] = {
      x: cO * x1 - sO * y2,
      y: z2,
      z: sO * x1 + cO * y2,
    };
  }
  return pts;
}