src/physics/events.js

// ————————————————————————————————————————————————
// Astronomical event detection on the JD timeline.
// Phase 2 of #6.
//
// All scanners take a starting JD and a max-window (in days). They return
// either { jd, ...metadata } or `null` if no event is found inside the
// window. The scanners are pure: they do not mutate state and they call
// only into formulas.js + the planets data.
// ————————————————————————————————————————————————

import { keplerHeliocentric, meanAnomalyAtJD, elementsAtJD } from './formulas.js';
import { planetById, PLANETS } from '../data/planets.js';
import { phaseAngle, lunarPosition } from './moon.js';
import { PHASE_RATE } from '../data/moon.js';

const TWO_PI = 2 * Math.PI;
const SUN_ANGULAR_RADIUS_RAD = ((16 / 60) * Math.PI) / 180; // ~16 arcmin

// Wrap an angle to (-π, π].
function wrapPi(a) {
  let x = (a + Math.PI) % TWO_PI;
  if (x < 0) x += TWO_PI;
  return x - Math.PI;
}

// Heliocentric (ecliptic) position of a planet at JD. Returns {x, y, z}.
export function helioPositionAtJD(planet, jd) {
  if (planet.isStar) return { x: 0, y: 0, z: 0 };
  const M = meanAnomalyAtJD(planet, jd);
  const el = elementsAtJD(planet, jd);
  return keplerHeliocentric(el.a, el.e, M, el.argPerihelionRad, el.ascNodeRad, el.iRad);
}

// Geocentric position of a planet (heliocentric of planet minus heliocentric
// of Earth). Returns {x, y, z}.
export function geocentricPosition(planet, jd) {
  const earth = planetById('earth');
  const p = helioPositionAtJD(planet, jd);
  const e = helioPositionAtJD(earth, jd);
  return { x: p.x - e.x, y: p.y - e.y, z: p.z - e.z };
}

// Heliocentric ecliptic longitude (radians) — the projection angle in the
// orbital plane. Range (-π, π].
export function helioLongitude(planet, jd) {
  const p = helioPositionAtJD(planet, jd);
  return Math.atan2(p.z, p.x);
}

// Sun-Earth-planet angle (Earth's elongation from the Sun towards `planet`).
// Returns radians in [0, π].
export function sunEarthPlanetAngle(planet, jd) {
  const earth = planetById('earth');
  const eHelio = helioPositionAtJD(earth, jd);
  const sunFromEarth = { x: -eHelio.x, y: -eHelio.y, z: -eHelio.z };
  const planetFromEarth = geocentricPosition(planet, jd);
  return vectorAngle(sunFromEarth, planetFromEarth);
}

// Geocentric angular separation between two planets seen from Earth (rad).
export function geocentricAngle(planetA, planetB, jd) {
  const a = geocentricPosition(planetA, jd);
  const b = geocentricPosition(planetB, jd);
  return vectorAngle(a, b);
}

function vectorAngle(a, b) {
  const dot = a.x * b.x + a.y * b.y + a.z * b.z;
  const na = Math.sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
  const nb = Math.sqrt(b.x * b.x + b.y * b.y + b.z * b.z);
  if (na === 0 || nb === 0) return 0;
  // Clamp against floating-point excursions.
  const c = Math.max(-1, Math.min(1, dot / (na * nb)));
  return Math.acos(c);
}

// Scan for a *local minimum* of `scalarFn(jd)` below `threshold`.
// `direction` = +1 walks forward in time (returns the *first* min after
// fromJD), −1 walks backward (returns the *first* min before fromJD).
// Coarse step locates the bucket; golden-section search refines.
function scanForLocalMin(
  scalarFn,
  fromJD,
  windowDays,
  stepDays,
  threshold = Infinity,
  direction = 1
) {
  const sign = direction >= 0 ? 1 : -1;
  let bestT = null;
  let bestV = Infinity;
  let t2 = fromJD;
  let t1 = fromJD + sign * stepDays;
  let prev2 = scalarFn(t2);
  let prev1 = scalarFn(t1);
  const end = fromJD + sign * windowDays;
  for (let t0 = fromJD + sign * 2 * stepDays; sign * (t0 - end) <= 0; t0 += sign * stepDays) {
    const cur = scalarFn(t0);
    if (prev1 < prev2 && prev1 < cur) {
      // Refine via golden section in [min(t2,t0), max(t2,t0)].
      const lo = Math.min(t2, t0),
        hi = Math.max(t2, t0);
      const t = goldenSectionMin(scalarFn, lo, hi);
      const v = scalarFn(t);
      if (v < bestV && v < threshold) {
        return { jd: t, value: v };
      }
    }
    prev2 = prev1;
    prev1 = cur;
    t2 = t1;
    t1 = t0;
  }
  return bestT == null ? null : { jd: bestT, value: bestV };
}

function goldenSectionMin(f, a, b, tol = 1e-3, maxIter = 60) {
  const phi = (Math.sqrt(5) - 1) / 2;
  let c = b - phi * (b - a);
  let d = a + phi * (b - a);
  let fc = f(c),
    fd = f(d);
  for (let i = 0; i < maxIter && b - a > tol; i++) {
    if (fc < fd) {
      b = d;
      d = c;
      fd = fc;
      c = b - phi * (b - a);
      fc = f(c);
    } else {
      a = c;
      c = d;
      fc = fd;
      d = a + phi * (b - a);
      fd = f(d);
    }
  }
  return 0.5 * (a + b);
}

// ————————————————————————————————————————————————
// Public scanners
// ————————————————————————————————————————————————

// Δλ-cosine. cos(Δλ) = +1 when planet and Earth share heliocentric
// longitude (opposition for outer / inferior conjunction for inner);
// cos(Δλ) = −1 at superior conjunction. Avoids the (-π, π] wrap that
// trips a sign-change scanner near superior conjunction.
function cosDeltaLambda(planet, jd) {
  const lp = helioLongitude(planet, jd);
  const le = helioLongitude(planetById('earth'), jd);
  return Math.cos(lp - le);
}

// Next (or previous) opposition of an *outer* planet: local maximum of
// cos(Δλ) (i.e. local minimum of −cos). `direction` = ±1 toggles search
// direction in time. Returns { jd, planet } or null.
export function nextOpposition(planet, fromJD, windowDays = 4 * 365, { direction = 1 } = {}) {
  const earth = planetById('earth');
  if (planet.id === 'earth' || planet.isStar) return null;
  if (planet.orbitKm < earth.orbitKm) return null; // outer only
  const f = (jd) => -cosDeltaLambda(planet, jd);
  const hit = scanForLocalMin(
    f,
    fromJD,
    windowDays,
    5,
    -Math.cos((0.5 * Math.PI) / 180),
    direction
  );
  if (hit == null) return null;
  return { jd: hit.jd, planet };
}

// Next (or previous) inferior conjunction of an *inner* planet
// (Mercury / Venus). Local maximum of cos(Δλ); planet sits between Sun
// and Earth.
export function nextInferiorConjunction(
  planet,
  fromJD,
  windowDays = 4 * 365,
  { direction = 1 } = {}
) {
  const earth = planetById('earth');
  if (planet.id === 'earth' || planet.orbitKm >= earth.orbitKm) return null;
  const f = (jd) => -cosDeltaLambda(planet, jd);
  const hit = scanForLocalMin(
    f,
    fromJD,
    windowDays,
    5,
    -Math.cos((0.5 * Math.PI) / 180),
    direction
  );
  if (hit == null) return null;
  return { jd: hit.jd, planet, kind: 'inferior' };
}

// Next (or previous) superior conjunction of any non-Earth planet. Local
// minimum of cos(Δλ); planet on the far side of the Sun from Earth.
export function nextSuperiorConjunction(
  planet,
  fromJD,
  windowDays = 4 * 365,
  { direction = 1 } = {}
) {
  if (planet.id === 'earth' || planet.isStar) return null;
  const f = (jd) => cosDeltaLambda(planet, jd);
  const hit = scanForLocalMin(
    f,
    fromJD,
    windowDays,
    5,
    Math.cos((179.5 * Math.PI) / 180),
    direction
  );
  if (hit == null) return null;
  return { jd: hit.jd, planet, kind: 'superior' };
}

// Next (or previous) geocentric conjunction between two planets (angular
// separation drops below `thresholdRad`). Local minimum of geocentricAngle.
export function nextGeocentricConjunction(
  planetA,
  planetB,
  fromJD,
  {
    windowDays = 4 * 365,
    thresholdRad = (1 * Math.PI) / 180, // 1° default
    direction = 1,
  } = {}
) {
  if (planetA.id === planetB.id) return null;
  const f = (jd) => geocentricAngle(planetA, planetB, jd);
  const hit = scanForLocalMin(f, fromJD, windowDays, 5, thresholdRad, direction);
  if (hit == null) return null;
  return { jd: hit.jd, planetA, planetB, separationRad: hit.value };
}

// Next (or previous) heliocentric multi-planet alignment: local minimum
// of the std-dev of the chosen planets' heliocentric longitudes.
export function nextHeliocentricAlignment(
  planetIds,
  fromJD,
  {
    windowDays = 50 * 365,
    thresholdRad = (5 * Math.PI) / 180, // 5° default
    direction = 1,
  } = {}
) {
  const planets = planetIds.map(planetById).filter(Boolean);
  if (planets.length < 2) return null;
  const f = (jd) => longitudeStdDev(planets.map((p) => helioLongitude(p, jd)));
  const hit = scanForLocalMin(f, fromJD, windowDays, 30, thresholdRad, direction);
  if (hit == null) return null;
  return { jd: hit.jd, planetIds, stdRad: hit.value };
}

// Std-dev of a set of angles, computed on the unwrapped values around the
// circular mean. Good enough for "are these things clustered" — for
// multi-decade spreads we want a different metric, but multi-planet
// alignments cluster within a small arc by definition.
function longitudeStdDev(longitudes) {
  if (longitudes.length === 0) return 0;
  // Circular mean.
  let sx = 0,
    sy = 0;
  for (const l of longitudes) {
    sx += Math.cos(l);
    sy += Math.sin(l);
  }
  const mean = Math.atan2(sy / longitudes.length, sx / longitudes.length);
  let s = 0;
  for (const l of longitudes) {
    const d = wrapPi(l - mean);
    s += d * d;
  }
  return Math.sqrt(s / longitudes.length);
}

// Next (or previous) transit (Mercury / Venus across the Sun's disc, seen
// from Earth). Walks through inferior conjunctions; at each one, checks
// whether the geocentric angular distance from the Sun is inside the
// Sun's disc.
export function nextTransit(planet, fromJD, windowDays = 50 * 365, { direction = 1 } = {}) {
  const earth = planetById('earth');
  if (planet.orbitKm >= earth.orbitKm) return null;
  const sign = direction >= 0 ? 1 : -1;
  let cursor = fromJD;
  let remaining = windowDays;
  while (remaining > 0) {
    const ic = nextInferiorConjunction(planet, cursor, remaining, { direction });
    if (!ic) return null;
    const eHelio = helioPositionAtJD(earth, ic.jd);
    const sunFromEarth = { x: -eHelio.x, y: -eHelio.y, z: -eHelio.z };
    const planetFromEarth = geocentricPosition(planet, ic.jd);
    const sep = vectorAngle(sunFromEarth, planetFromEarth);
    if (sep < SUN_ANGULAR_RADIUS_RAD) {
      return { jd: ic.jd, planet, separationRad: sep };
    }
    // Step past this conjunction in the current scan direction.
    const newCursor = ic.jd + sign * 1;
    remaining -= sign * (newCursor - cursor);
    cursor = newCursor;
  }
  return null;
}

// Convenience: tag every orbital planet as 'inner' / 'outer' relative to Earth.
export function isOuterPlanet(planet) {
  const earth = planetById('earth');
  return planet.orbitKm > earth.orbitKm;
}
export function isInnerPlanet(planet) {
  const earth = planetById('earth');
  return planet.orbitKm > 0 && planet.orbitKm < earth.orbitKm;
}

// Re-export the static catalog of orbital planets for UI dropdowns.
export const ORBITAL_PLANETS = PLANETS.filter((p) => !p.isStar && p.id !== 'earth');

const PHASE_TARGETS = {
  new: 0,
  firstQuarter: Math.PI / 2,
  full: Math.PI,
  lastQuarter: (3 * Math.PI) / 2,
};

// Newton iteration on the Meeus-precision phaseAngle. The first step is the
// synodic-rate estimate from fromJD's current phase to the next target in the
// chosen direction (≤ ~30 days). Subsequent steps clamp to ±3 days so a
// correction can't accidentally jump across a lunation. Converges to ~1e-7
// rad in 3–4 iterations.

const NEWTON_TOL_RAD = 1e-9;
const NEWTON_MAX_ITER = 8;
const NEWTON_MAX_STEP_DAYS = 3;

export function nextLunarPhase(kind, fromJD, windowDays = 365, { direction = 1 } = {}) {
  const target = PHASE_TARGETS[kind];
  if (target === undefined) return null;
  const current = phaseAngle(fromJD);
  // Normalize delta to (-π, π] so the EPS rollover behaves correctly
  // regardless of which side of the target Newton converged to last call.
  // Without this, a previous Newton pass that landed at phase = 2π − tiny
  // (i.e. just shy of target=0 in unwrapped form) yields delta = −2π + tiny,
  // which < EPS triggers a +2π rollover that collapses to a tiny advance —
  // and Newton re-converges to the same lunation (#87).
  let delta = target - current;
  delta = ((((delta + Math.PI) % (2 * Math.PI)) + 2 * Math.PI) % (2 * Math.PI)) - Math.PI;
  // Force advance to the *next* target in `direction`. EPS << real phase
  // advance (~0.21 rad/day) but >> Newton tolerance, so a fromJD already
  // at the target rolls over to the next lunation instead of returning ~0.
  const EPS = 1e-6;
  if (direction >= 0) {
    if (delta < EPS) delta += 2 * Math.PI;
  } else {
    if (delta > -EPS) delta -= 2 * Math.PI;
  }
  let jd = fromJD + delta / PHASE_RATE;
  for (let i = 0; i < NEWTON_MAX_ITER; i++) {
    const phase = phaseAngle(jd);
    let d = target - phase;
    if (d > Math.PI) d -= 2 * Math.PI;
    else if (d < -Math.PI) d += 2 * Math.PI;
    if (Math.abs(d) < NEWTON_TOL_RAD) break;
    let step = d / PHASE_RATE;
    if (step > NEWTON_MAX_STEP_DAYS) step = NEWTON_MAX_STEP_DAYS;
    else if (step < -NEWTON_MAX_STEP_DAYS) step = -NEWTON_MAX_STEP_DAYS;
    jd += step;
  }
  if (Math.abs(jd - fromJD) > windowDays) return null;
  return { jd, kind };
}

// ─────────────────────────────────────────────────────────────────────────
// Solstice / equinox scanner (#82). The four cardinal solar events
// correspond to the Sun's geocentric ecliptic longitude crossing 0,
// π/2, π, 3π/2. Sun's geocentric longitude = Earth's heliocentric
// longitude + π.
//
// Newton iteration on the residual; same direction/window contract as
// nextLunarPhase.
// ─────────────────────────────────────────────────────────────────────────

const SOLAR_EVENT_TARGETS = {
  marchEquinox: 0,
  juneSolstice: Math.PI / 2,
  septemberEquinox: Math.PI,
  decemberSolstice: (3 * Math.PI) / 2,
};

// Sun's geocentric ecliptic longitude (rad), wrapped to [0, 2π).
function sunGeocentricLongitude(jd) {
  const earth = planetById('earth');
  let lon = helioLongitude(earth, jd) + Math.PI;
  // helioLongitude returns (-π, π]; +π pushes us into (0, 2π]. Wrap to [0, 2π).
  lon = ((lon % (2 * Math.PI)) + 2 * Math.PI) % (2 * Math.PI);
  return lon;
}

// Earth's mean longitude rate, expressed in rad/day (the data file stores
// rad per Julian century).
const EARTH_MEAN_MOTION_RAD_PER_DAY = (() => {
  const earth = planetById('earth');
  return earth.meanLongitudeRadRate / 36525;
})();

const SOLAR_NEWTON_TOL_RAD = 1e-9;
const SOLAR_NEWTON_MAX_ITER = 8;
const SOLAR_NEWTON_MAX_STEP_DAYS = 5;

export function nextSolarEvent(kind, fromJD, windowDays = 365, { direction = 1 } = {}) {
  const target = SOLAR_EVENT_TARGETS[kind];
  if (target === undefined) return null;
  const current = sunGeocentricLongitude(fromJD);
  // Normalize delta to (-π, π] before the rollover — see nextLunarPhase for
  // the same wrap-around stall (#87).
  let delta = target - current;
  delta = ((((delta + Math.PI) % (2 * Math.PI)) + 2 * Math.PI) % (2 * Math.PI)) - Math.PI;
  // Force advance to the *next* target in `direction` (same EPS rollover
  // logic as nextLunarPhase).
  const EPS = 1e-6;
  if (direction >= 0) {
    if (delta < EPS) delta += 2 * Math.PI;
  } else {
    if (delta > -EPS) delta -= 2 * Math.PI;
  }
  let jd = fromJD + delta / EARTH_MEAN_MOTION_RAD_PER_DAY;
  for (let i = 0; i < SOLAR_NEWTON_MAX_ITER; i++) {
    const lon = sunGeocentricLongitude(jd);
    let d = target - lon;
    if (d > Math.PI) d -= 2 * Math.PI;
    else if (d < -Math.PI) d += 2 * Math.PI;
    if (Math.abs(d) < SOLAR_NEWTON_TOL_RAD) break;
    let step = d / EARTH_MEAN_MOTION_RAD_PER_DAY;
    if (step > SOLAR_NEWTON_MAX_STEP_DAYS) step = SOLAR_NEWTON_MAX_STEP_DAYS;
    else if (step < -SOLAR_NEWTON_MAX_STEP_DAYS) step = -SOLAR_NEWTON_MAX_STEP_DAYS;
    jd += step;
  }
  if (Math.abs(jd - fromJD) > windowDays) return null;
  return { jd, kind };
}

// ─────────────────────────────────────────────────────────────────────────
// Eclipse scanners (#17). A lunar eclipse happens at full moon when the
// Moon is close enough to the ecliptic that it enters Earth's shadow
// cone; a solar eclipse happens at new moon when the Moon crosses the
// Sun-Earth line within the Moon-disc + Sun-disc + parallax window.
//
// We piggyback on nextLunarPhase to step through new/full instants, then
// reject lunations whose geocentric ecliptic latitude exceeds the
// eclipse-possibility threshold. ~1.5° covers any-contact for both
// types (lunar penumbra + Moon disc; solar disc + Moon disc + parallax).
// ─────────────────────────────────────────────────────────────────────────

// 1.5° in radians — generous "any eclipse visible somewhere on Earth"
// threshold. Sources: Meeus §54 (penumbra geometry), §51 (solar parallax).
const ECLIPSE_LAT_THRESHOLD_RAD = (1.5 * Math.PI) / 180;

function nextEclipse(kind, fromJD, windowDays, direction) {
  const phaseKind = kind === 'lunar' ? 'full' : 'new';
  let cursor = fromJD;
  // Walk lunations until we find one that qualifies or the window runs out.
  // ~13 lunations per year covers any plausible windowDays argument.
  const MAX_STEPS = Math.ceil(windowDays / 29 + 5);
  for (let i = 0; i < MAX_STEPS; i++) {
    const seed = nextLunarPhase(phaseKind, cursor, windowDays, { direction });
    if (!seed) return null;
    const { latRad } = lunarPosition(seed.jd);
    if (Math.abs(latRad) < ECLIPSE_LAT_THRESHOLD_RAD) {
      return { jd: seed.jd, kind, latRad };
    }
    // Step past this lunation so the next nextLunarPhase moves on.
    cursor = seed.jd + (direction >= 0 ? 1 : -1);
  }
  return null;
}

export function nextLunarEclipse(fromJD, windowDays = 3650, { direction = 1 } = {}) {
  return nextEclipse('lunar', fromJD, windowDays, direction);
}

export function nextSolarEclipse(fromJD, windowDays = 3650, { direction = 1 } = {}) {
  return nextEclipse('solar', fromJD, windowDays, direction);
}