KSP-MissionControl/apps/live-map/src/calculators/eclipse.ts

/**
 * Eclipse calculator — finds the next N times when the sun (Kerbol)
 * is occluded from an observer's point of view by another body.
 *
 * Algorithm: scan forward in UT at fixed steps, refining around
 * the transition points with bisection. Treats the sun as a point
 * source (parallel rays) — accurate for KSP scale.
 */
import type { CelestialBody } from '@kerbal-rt/shared-types';
import { bodyPositionAt } from '../scene/layout.js';
import { shadowFraction } from '@kerbal-rt/orbital-math';

export interface EclipseWindow {
  /** UT at which the eclipse begins (sun starts being occluded). */
  utStart: number;
  /** UT at which the eclipse ends. */
  utEnd: number;
  /** Maximum shadow fraction during the window (0..1). */
  maxFraction: number;
  /** UT at which the max occurs. */
  utPeak: number;
}

export interface EclipseOptions {
  observerId: string;
  eclipserId: string;
  /** Scan starts at this UT. */
  startUt: number;
  /** How many windows to find. */
  count?: number;
  /** Max time to scan forward (default ~1 KSP year). */
  maxSearchTime?: number;
  /** Coarse scan step (default 600 = 10 KSP minutes). */
  stepSec?: number;
  /** Threshold for "eclipse started". */
  threshold?: number;
}

/**
 * Find the next `count` eclipse windows where `eclipser` occludes
 * the sun as seen from `observer`.
 *
 * Returns an empty array if observer and eclipser are the same body
 * or if either is Kerbol (the sun can't eclipse itself, and
 * eclipsers behind the sun don't make sense).
 */
export function findEclipseWindows(
  bodies: CelestialBody[],
  options: EclipseOptions,
): EclipseWindow[] {
  const {
    observerId,
    eclipserId,
    startUt,
    count = 3,
    maxSearchTime = 426 * 6 * 3600, // 1 KSP year
    stepSec = 600,
    threshold = 0.05,
  } = options;

  if (observerId === eclipserId) return [];
  const observer = bodies.find((b) => b.id === observerId);
  const eclipser = bodies.find((b) => b.id === eclipserId);
  if (!observer || !eclipser) return [];

  const windows: EclipseWindow[] = [];
  const end = startUt + maxSearchTime;

  // Coarse scan: find the start of an eclipse (transition from below
  // threshold to above threshold). Then refine to find utStart, utPeak,
  // utEnd.
  let inEclipse = false;
  let current: Partial<EclipseWindow> = {};

  // If we start in the middle of an eclipse, bisect backwards to find utStart.
  const startF = computeShadowFraction(bodies, observerId, eclipserId, startUt);
  if (startF > threshold) {
    inEclipse = true;
    const utStart = refine(
      bodies,
      observerId,
      eclipserId,
      Math.max(0, startUt - stepSec),
      startUt,
      threshold,
      'down',
    );
    current = { utStart };
  }

  let t = startUt + stepSec;
  while (t < end && windows.length < count) {
    const f = computeShadowFraction(bodies, observerId, eclipserId, t);
    if (!inEclipse && f > threshold) {
      // Eclipse just started — bisect to find exact start
      const utStart = refine(bodies, observerId, eclipserId, t - stepSec, t, threshold, 'up');
      inEclipse = true;
      current = { utStart };
    } else if (inEclipse && f < threshold) {
      // Eclipse just ended
      const utEnd = refine(bodies, observerId, eclipserId, t - stepSec, t, threshold, 'down');
      const utPeak = refine(
        bodies,
        observerId,
        eclipserId,
        current.utStart!,
        utEnd,
        threshold,
        'max',
      );
      const maxFraction = computeShadowFraction(bodies, observerId, eclipserId, utPeak);
      windows.push({
        utStart: current.utStart!,
        utPeak,
        utEnd,
        maxFraction,
      });
      inEclipse = false;
      current = {};
    }
    t += stepSec;
  }

  // If we ended while still in eclipse, close it.
  if (inEclipse) {
    const utEnd = refine(
      bodies,
      observerId,
      eclipserId,
      t - stepSec,
      t + stepSec,
      threshold,
      'down',
    );
    const utPeak = refine(
      bodies,
      observerId,
      eclipserId,
      current.utStart!,
      utEnd,
      threshold,
      'max',
    );
    const maxFraction = computeShadowFraction(bodies, observerId, eclipserId, utPeak);
    windows.push({
      utStart: current.utStart!,
      utPeak,
      utEnd,
      maxFraction,
    });
  }

  return windows;
}

/** Shadow fraction at a single instant. */
export function computeShadowFraction(
  bodies: CelestialBody[],
  observerId: string,
  eclipserId: string,
  ut: number,
): number {
  const observer = bodies.find((b) => b.id === observerId);
  const eclipser = bodies.find((b) => b.id === eclipserId);
  if (!observer || !eclipser) return 0;

  // Sun is Kerbol (the system root). For this to work the catalog
  // must have a single body with parentId === null — the star.
  const sun = bodies.find((b) => b.parentId === null);
  if (!sun) return 0;

  // The reference frame for shadowFraction is: vectors from the
  // observer toward the sun, and from the observer toward the
  // eclipser. We compute them in heliocentric coordinates then
  // shift.
  const sunPos = bodyPositionAt(bodies, sun.id, ut);
  const eclipserPos = bodyPositionAt(bodies, eclipserId, ut);
  const observerPos = bodyPositionAt(bodies, observerId, ut);

  // For solar eclipses, "eclipser" sits between observer and sun.
  // We treat the eclipser as the occluder and the sun as the light
  // source. shadowFraction() expects: vector from observer to sun,
  // and vector from observer to eclipser.
  // Observer-to-sun = sunPos - observerPos
  // Observer-to-eclipser = eclipserPos - observerPos
  // BUT shadowFraction actually wants: vector from observer to sun,
  // and the eclipser center RELATIVE to the observer-to-sun line.
  // Re-reading the function: it computes perpDist of occluder
  // center to sun-direction line, and uses the sun-direction
  // projection. So the right call is:
  return shadowFraction(
    { x: sunPos.x - observerPos.x, y: sunPos.y - observerPos.y, z: sunPos.z - observerPos.z },
    {
      x: eclipserPos.x - observerPos.x,
      y: eclipserPos.y - observerPos.y,
      z: eclipserPos.z - observerPos.z,
    },
    eclipser.radius,
  );
}

/**
 * Bisection to find when the shadow fraction crosses the threshold.
 * direction = 'up'  → find t such that f(t) ≈ threshold going upward
 * direction = 'down' → find t such that f(t) ≈ threshold going downward
 * direction = 'max' → find t that maximizes f in the range
 */
function refine(
  bodies: CelestialBody[],
  observerId: string,
  eclipserId: string,
  tLo: number,
  tHi: number,
  threshold: number,
  direction: 'up' | 'down' | 'max',
): number {
  if (direction === 'max') {
    // Golden-section or simple ternary search on a small interval.
    // For 600s windows this is plenty.
    let lo = tLo;
    let hi = tHi;
    for (let i = 0; i < 32; i++) {
      const m1 = lo + (hi - lo) / 3;
      const m2 = hi - (hi - lo) / 3;
      const f1 = computeShadowFraction(bodies, observerId, eclipserId, m1);
      const f2 = computeShadowFraction(bodies, observerId, eclipserId, m2);
      if (f1 < f2) lo = m1;
      else hi = m2;
    }
    return (lo + hi) / 2;
  }

  // Bisection: 30 iterations gets us to 1-second precision on a 600s range.
  let lo = tLo;
  let hi = tHi;
  for (let i = 0; i < 30; i++) {
    const mid = (lo + hi) / 2;
    const f = computeShadowFraction(bodies, observerId, eclipserId, mid);
    if (direction === 'up' ? f < threshold : f > threshold) lo = mid;
    else hi = mid;
  }
  return (lo + hi) / 2;
}