src/physics/seasons.js

// ————————————————————————————————————————————————
// Solar declination and seasons (#82).
//
// Pure functions — solar declination as a closed-form expression of
// Earth's heliocentric longitude, season name from the four cardinal
// solar events bracketing fromJD.
// ————————————————————————————————————————————————

import { planetById } from '../data/planets.js';
import { helioLongitude, nextSolarEvent } from './events.js';

// Sun's geocentric ecliptic longitude (rad), wrapped to [0, 2π).
// (Local copy — events.js doesn't export this helper.)
function sunGeocentricLongitude(jd) {
  const earth = planetById('earth');
  let lon = helioLongitude(earth, jd) + Math.PI;
  lon = ((lon % (2 * Math.PI)) + 2 * Math.PI) % (2 * Math.PI);
  return lon;
}

// Earth's obliquity in radians. Read once at module load — the JPL fit
// holds this constant within the simulation's accuracy budget.
const EARTH_TILT_RAD = (() => {
  const earth = planetById('earth');
  return earth.axialTiltRad;
})();

/**
 * Solar declination δ at JD: the latitude on Earth where the Sun is
 * directly overhead at local noon.
 *
 *   δ = arcsin(sin(ε) · sin(L_sun))
 *
 * where ε is Earth's obliquity (~23.44°) and L_sun is the Sun's
 * geocentric ecliptic longitude.
 *
 * Returns radians in [-ε, +ε]. Accuracy: better than ~0.01° for any JD
 * inside the J2000-fit window (1800–2050), since this rides on the same
 * JPL pipeline as the rest of the simulation.
 */
export function solarDeclinationAtJD(jd) {
  const lonSun = sunGeocentricLongitude(jd);
  return Math.asin(Math.sin(EARTH_TILT_RAD) * Math.sin(lonSun));
}

const SEASON_FROM_KIND = {
  marchEquinox: 'spring',
  juneSolstice: 'summer',
  septemberEquinox: 'autumn',
  decemberSolstice: 'winter',
};

/**
 * Northern-hemisphere season at JD. Returns 'spring' | 'summer' |
 * 'autumn' | 'winter'.
 *
 * Implementation: the season at any JD is the one that *started* at
 * the most recent cardinal event. We scan backward for each of the
 * four kinds and take the one whose JD is largest (≤ fromJD).
 */
export function seasonAtJD(jd) {
  let bestJD = -Infinity;
  let bestKind = null;
  for (const kind of Object.keys(SEASON_FROM_KIND)) {
    const r = nextSolarEvent(kind, jd, 400, { direction: -1 });
    if (r && r.jd > bestJD) {
      bestJD = r.jd;
      bestKind = kind;
    }
  }
  return bestKind ? SEASON_FROM_KIND[bestKind] : null;
}

/**
 * JD of the next sunset at a geographic location, scanning forward from
 * fromJD. Returns null if the location is in polar night or polar day
 * (no sunset within the next two UTC days).
 *
 * Closed-form: hour angle H from solar noon satisfies
 *   cos(H) = -tan(φ) · tan(δ)
 * where φ = latitude, δ = solar declination at local noon. The sunset
 * UTC clock time is then 12 − λ/15 + H/15° hours, where λ = longitude
 * (east positive).
 *
 * Approximation: ignores equation of time and atmospheric refraction.
 * Accuracy ≈ ±15 minutes — sufficient for narrative framing where the
 * goal is "the sun is on the horizon", not navigation.
 */
export function nextSunsetAtLocation({ latDeg, lonDeg }, fromJD) {
  const latRad = (latDeg * Math.PI) / 180;
  for (let dayOffset = 0; dayOffset < 2; dayOffset++) {
    const midnightJD = Math.floor(fromJD - 0.5) + 0.5 + dayOffset;
    const declRad = solarDeclinationAtJD(midnightJD + 0.5); // sample at local noon
    const cosH = -Math.tan(latRad) * Math.tan(declRad);
    if (cosH > 1 || cosH < -1) continue; // polar night or midnight sun
    const HHours = (Math.acos(cosH) * 12) / Math.PI;
    const sunsetUTCHours = 12 - lonDeg / 15 + HHours;
    const sunsetJD = midnightJD + sunsetUTCHours / 24;
    if (sunsetJD > fromJD) return sunsetJD;
  }
  return null;
}