/**
 * Mercury retrograde windows, computed rather than copied.
 *
 * The site's rule is that a published number is measured or computed from a
 * stated model. Retrograde dates are easy to find on the web and impossible to
 * verify by looking at them, so they are derived here from orbital elements and
 * the derivation is testable: if the arithmetic were wrong, the intervals
 * between successive retrogrades would not land on Mercury's known synodic
 * period of 115.88 days, and the tests would say so.
 *
 * Mercury is not going backwards. It appears to because Earth overtakes it on
 * the inside; the apparent reversal is a projection effect of two bodies moving
 * at different speeds around the same star. That is what makes the audit
 * interesting: the windows are perfectly real and perfectly predictable, and
 * nothing about them touches a flight schedule.
 *
 * Elements: JPL/Standish approximate mean elements, J2000 with linear rates.
 * Accurate to well under a degree over 1800–2050 — far tighter than needed to
 * find a station to the day.
 */

type Elements = {
  a: number; // semi-major axis, AU
  e: number; // eccentricity
  i: number; // inclination, degrees
  L: number; // mean longitude, degrees
  peri: number; // longitude of perihelion, degrees
  node: number; // longitude of ascending node, degrees
};

const MERCURY = {
  base: { a: 0.38709927, e: 0.20563593, i: 7.00497902, L: 252.2503235, peri: 77.45779628, node: 48.33076593 },
  rate: { a: 0.00000037, e: 0.00001906, i: -0.00594749, L: 149472.67411175, peri: 0.16047689, node: -0.12534081 },
};

/** Earth–Moon barycentre; the Moon's wobble is irrelevant at this precision. */
const EARTH = {
  base: { a: 1.00000261, e: 0.01671123, i: -0.00001531, L: 100.46457166, peri: 102.93768193, node: 0 },
  rate: { a: 0.00000562, e: -0.00004392, i: -0.01294668, L: 35999.37244981, peri: 0.32327364, node: 0 },
};

const DEG = Math.PI / 180;

/** Julian centuries since J2000.0 (2000-01-01T12:00:00Z). */
export function julianCenturies(date: Date): number {
  const jd = date.getTime() / 86_400_000 + 2440587.5;
  return (jd - 2451545.0) / 36525;
}

function elementsAt(body: typeof MERCURY, t: number): Elements {
  return {
    a: body.base.a + body.rate.a * t,
    e: body.base.e + body.rate.e * t,
    i: body.base.i + body.rate.i * t,
    L: body.base.L + body.rate.L * t,
    peri: body.base.peri + body.rate.peri * t,
    node: body.base.node + body.rate.node * t,
  };
}

/** Solve M = E − e·sin E for E. Newton converges in a few steps at these e. */
function eccentricAnomaly(meanAnomalyDeg: number, e: number): number {
  const M = meanAnomalyDeg * DEG;
  let E = M + e * Math.sin(M);
  for (let i = 0; i < 12; i++) {
    const dE = (E - e * Math.sin(E) - M) / (1 - e * Math.cos(E));
    E -= dE;
    if (Math.abs(dE) < 1e-12) break;
  }
  return E;
}

/** Heliocentric ecliptic rectangular coordinates, AU. */
function heliocentric(body: typeof MERCURY, t: number): { x: number; y: number; z: number } {
  const el = elementsAt(body, t);

  // Mean anomaly, wrapped to ±180° before solving.
  const M = ((el.L - el.peri) % 360 + 540) % 360 - 180;
  const E = eccentricAnomaly(M, el.e);

  // Position in the orbital plane.
  const xOrb = el.a * (Math.cos(E) - el.e);
  const yOrb = el.a * Math.sqrt(1 - el.e * el.e) * Math.sin(E);

  const omega = (el.peri - el.node) * DEG; // argument of perihelion
  const node = el.node * DEG;
  const inc = el.i * DEG;

  const cosW = Math.cos(omega), sinW = Math.sin(omega);
  const cosN = Math.cos(node), sinN = Math.sin(node);
  const cosI = Math.cos(inc), sinI = Math.sin(inc);

  return {
    x: (cosW * cosN - sinW * sinN * cosI) * xOrb + (-sinW * cosN - cosW * sinN * cosI) * yOrb,
    y: (cosW * sinN + sinW * cosN * cosI) * xOrb + (-sinW * sinN + cosW * cosN * cosI) * yOrb,
    z: sinW * sinI * xOrb + cosW * sinI * yOrb,
  };
}

/**
 * Geocentric apparent ecliptic longitude of Mercury, degrees.
 *
 * Geometric, not corrected for light-time or aberration. Both shift the answer
 * by far less than the day-scale resolution this is used at.
 */
export function mercuryLongitude(date: Date): number {
  const t = julianCenturies(date);
  const m = heliocentric(MERCURY, t);
  const e = heliocentric(EARTH, t);
  const lon = Math.atan2(m.y - e.y, m.x - e.x) / DEG;
  return (lon % 360 + 360) % 360;
}

/**
 * Geocentric ecliptic longitude of the Sun, degrees — which is simply Earth's
 * heliocentric longitude plus 180°.
 *
 * Used to state Mercury's elongation. Mercury is an inferior planet, so it
 * never appears far from the Sun: that bound is what makes retrograde a
 * near-Sun event and is worth being able to check.
 */
export function sunLongitude(date: Date): number {
  const e = heliocentric(EARTH, julianCenturies(date));
  const lon = Math.atan2(-e.y, -e.x) / DEG;
  return (lon % 360 + 360) % 360;
}

/** Angular distance from the Sun, 0–180°. */
export function elongation(date: Date): number {
  const d = Math.abs(mercuryLongitude(date) - sunLongitude(date)) % 360;
  return d > 180 ? 360 - d : d;
}

/** Signed change in longitude, handling the 360°→0° wrap. */
function longitudeDelta(a: number, b: number): number {
  let d = b - a;
  if (d > 180) d -= 360;
  if (d < -180) d += 360;
  return d;
}

export type RetrogradeWindow = { start: Date; end: Date; days: number };

const DAY = 86_400_000;

/**
 * Every interval in [from, to] during which Mercury's geocentric longitude is
 * decreasing — the retrograde windows.
 *
 * Stations are found by scanning daily for a sign change, then bisecting to the
 * hour. Daily sampling is safe because retrograde lasts about three weeks;
 * nothing shorter can hide between samples.
 */
export function retrogradeWindows(from: Date, to: Date): RetrogradeWindow[] {
  const windows: RetrogradeWindow[] = [];

  const isRetro = (time: number) =>
    longitudeDelta(mercuryLongitude(new Date(time - DAY / 2)), mercuryLongitude(new Date(time + DAY / 2))) < 0;

  /** Bisect for the moment the motion changes direction. */
  const station = (before: number, after: number): Date => {
    let lo = before;
    let hi = after;
    const target = isRetro(after);
    for (let i = 0; i < 24; i++) {
      const mid = (lo + hi) / 2;
      if (isRetro(mid) === target) hi = mid;
      else lo = mid;
    }
    return new Date(hi);
  };

  // If the scan begins mid-retrograde, that window's real station lies before
  // `from` and is unknowable from here. Reporting it would publish the edge of
  // the scan as though it were an astronomical event, so it is skipped and only
  // complete station-to-station windows are returned. Same at the far end: a
  // window still open at `to` is never pushed.
  let previous = isRetro(from.getTime());
  let openedAt: Date | null = null;

  for (let time = from.getTime() + DAY; time <= to.getTime(); time += DAY) {
    const current = isRetro(time);
    if (current === previous) continue;

    if (current) {
      openedAt = station(time - DAY, time);
    } else if (openedAt) {
      const end = station(time - DAY, time);
      windows.push({
        start: openedAt,
        end,
        days: Math.round((end.getTime() - openedAt.getTime()) / DAY),
      });
      openedAt = null;
    }
    previous = current;
  }

  return windows;
}

/** The window containing `now`, if any. */
export function currentWindow(now: Date, windows: RetrogradeWindow[]): RetrogradeWindow | null {
  return windows.find((w) => now >= w.start && now <= w.end) ?? null;
}

/** Mercury's synodic period — the interval between successive inferior conjunctions. */
export const SYNODIC_PERIOD_DAYS = 115.88;