feat: integrate PR #451 pass timing with additive catalog voice queries

src/data/satellitePass.js

@@ -0,0 +1,421 @@
+// src/data/satellitePass.js
+/**
+ * Satellite pass prediction with naked-eye visibility.
+ * Look angles come from SGP4; rise/set are bisected to sub-second precision,
+ * the peak is refined with a parabola, and visibility uses the USNO
+ * low-precision Sun position with a cylindrical Earth shadow. Scan pattern from
+ * skylight (MIT) shared/src/celestial.ts nextISSPass.
+ */
+import { propagate, gstime, eciToEcf, ecfToLookAngles } from 'satellite.js';
+
+const R2D = 180 / Math.PI;
+const D2R = Math.PI / 180;
+const EARTH_RADIUS_KM = 6371.0;
+const AU_KM = 149597870.7;
+
+/**
+ * Observer look angles at an instant, or null when propagation fails.
+ * @param {Object} satrec SGP4 satellite record
+ * @param {number} dateMs UTC epoch timestamp in milliseconds
+ * @param {number} latDeg Observer latitude in degrees [-90, 90]
+ * @param {number} lonDeg Observer longitude in degrees [-180, 180]
+ * @returns {{ elevDeg: number, azDeg: number, satECI: {x:number, y:number, z:number} } | null}
+ */
+export function lookAnglesAt(satrec, dateMs, latDeg, lonDeg) {
+  const date = new Date(dateMs);
+  let pv;
+  try { pv = propagate(satrec, date); } catch { return null; }
+  const pos = pv && pv.position;
+  if (!pos || typeof pos === 'boolean' || ![pos.x, pos.y, pos.z].every(Number.isFinite)) return null;
+  const ecf = eciToEcf(pos, gstime(date));
+  const look = ecfToLookAngles(
+    { latitude: latDeg * D2R, longitude: lonDeg * D2R, height: 0 },
+    ecf,
+  );
+  return {
+    elevDeg: look.elevation * R2D,
+    azDeg: (((look.azimuth * R2D) % 360) + 360) % 360,
+    satECI: pos,
+  };
+}
+
+/**
+ * Compute Sun position in geocentric ECI coordinates using USNO low-precision solar coordinates.
+ * Accurate to within ~1 arcminute (~0.017°) within centuries of J2000.0.
+ * @param {number} dateMs UTC epoch timestamp in milliseconds
+ * @returns {{ x: number, y: number, z: number, raDeg: number, decDeg: number }}
+ */
+export function solarPositionECI(dateMs) {
+  const JD = dateMs / 86400000 + 2440587.5;
+  const D = JD - 2451545.0; // days since J2000.0
+  const g = ((((357.529 + 0.98560028 * D) % 360) + 360) % 360) * D2R;
+  const q = (((280.459 + 0.98564736 * D) % 360) + 360) % 360;
+  const L =
+    ((((q + 1.915 * Math.sin(g) + 0.02 * Math.sin(2 * g)) % 360) + 360) % 360) *
+    D2R;
+  const e = (23.439 - 0.00000036 * D) * D2R;
+  const ra = Math.atan2(Math.cos(e) * Math.sin(L), Math.cos(L));
+  const dec = Math.asin(Math.sin(e) * Math.sin(L));
+  return {
+    x: AU_KM * Math.cos(dec) * Math.cos(ra),
+    y: AU_KM * Math.cos(dec) * Math.sin(ra),
+    z: AU_KM * Math.sin(dec),
+    raDeg: (((ra * R2D) % 360) + 360) % 360,
+    decDeg: dec * R2D,
+  };
+}
+
+/**
+ * Check whether a satellite in ECI is illuminated by the Sun (outside Earth shadow).
+ * Uses a cylindrical Earth shadow model with mean radius 6,371 km.
+ * @param {{ x: number, y: number, z: number }} satECI Satellite position in km
+ * @param {number} dateMs UTC epoch timestamp in milliseconds
+ * @returns {boolean} True if sunlit, false if in eclipse (Earth shadow)
+ */
+export function isSatelliteSunlit(satECI, dateMs) {
+  if (!satECI || typeof satECI === 'boolean') return false;
+  const sun = solarPositionECI(dateMs);
+  const sunLen = Math.hypot(sun.x, sun.y, sun.z);
+  const sx = sun.x / sunLen;
+  const sy = sun.y / sunLen;
+  const sz = sun.z / sunLen;
+  const proj = satECI.x * sx + satECI.y * sy + satECI.z * sz;
+  if (proj > 0) return true; // Sunward side of Earth
+  const dx = satECI.x - proj * sx;
+  const dy = satECI.y - proj * sy;
+  const dz = satECI.z - proj * sz;
+  return Math.hypot(dx, dy, dz) > EARTH_RADIUS_KM;
+}
+
+/**
+ * Compute the observer local solar elevation angle in degrees.
+ * @param {number} latDeg Observer latitude in degrees
+ * @param {number} lonDeg Observer longitude in degrees
+ * @param {number} dateMs UTC epoch timestamp in milliseconds
+ * @returns {number} Solar elevation in degrees [-90, 90]
+ */
+export function observerSolarElevation(latDeg, lonDeg, dateMs) {
+  const sun = solarPositionECI(dateMs);
+  const gmst = gstime(new Date(dateMs));
+  const gmstDeg = (((gmst * R2D) % 360) + 360) % 360;
+  const localHourAngleDeg =
+    (((gmstDeg + lonDeg - sun.raDeg) % 360) + 360) % 360;
+  const H = localHourAngleDeg * D2R;
+  const phi = latDeg * D2R;
+  const dec = sun.decDeg * D2R;
+  const sinElev =
+    Math.sin(phi) * Math.sin(dec) + Math.cos(phi) * Math.cos(dec) * Math.cos(H);
+  return Math.asin(Math.max(-1, Math.min(1, sinElev))) * R2D;
+}
+
+/**
+ * Check whether the observer sky is dark (civil, nautical, or astronomical twilight / night).
+ * Default threshold is -6° (civil twilight).
+ * @param {number} latDeg Observer latitude
+ * @param {number} lonDeg Observer longitude
+ * @param {number} dateMs UTC timestamp
+ * @param {number} [maxSunElevDeg=-6] Sun elevation ceiling (default -6° for civil twilight)
+ * @returns {boolean} True if sun is at or below the threshold
+ */
+export function isObserverDark(latDeg, lonDeg, dateMs, maxSunElevDeg = -6) {
+  return observerSolarElevation(latDeg, lonDeg, dateMs) <= maxSunElevDeg;
+}
+
+/**
+ * Bisect the elevation threshold boundary f(t) = elev(t) - minElevDeg = 0.
+ * Assumes one threshold crossing inside the bracket. 7 iterations resolve 20s to ~0.16s.
+ * @private
+ */
+function _bisectBoundary(
+  satrec,
+  latDeg,
+  lonDeg,
+  tLow,
+  tHigh,
+  minElevDeg,
+  isRising,
+  steps = 10,
+) {
+  let low = tLow;
+  let high = tHigh;
+  for (let i = 0; i < steps; i++) {
+    const mid = (low + high) / 2;
+    const look = lookAnglesAt(satrec, mid, latDeg, lonDeg);
+    const elev = look ? look.elevDeg : -90;
+    if (isRising) {
+      if (elev >= minElevDeg) {
+        high = mid;
+      } else {
+        low = mid;
+      }
+    } else {
+      if (elev >= minElevDeg) {
+        low = mid;
+      } else {
+        high = mid;
+      }
+    }
+  }
+  return isRising ? Math.ceil(high) : Math.floor(low);
+}
+
+/**
+ * Find the next pass of a satellite over an observer location.
+ * Uses 20s coarse scanning, bisection for sub-second rise/set times,
+ * parabolic interpolation for the peak, and optional naked-eye visibility filtering.
+ *
+ * @param {Object} options
+ * @param {Object} options.satrec SGP4 satellite record
+ * @param {number} options.latDeg Observer latitude [-90, 90]
+ * @param {number} options.lonDeg Observer longitude [-180, 180]
+ * @param {number} options.fromMs UTC start time in milliseconds
+ * @param {number} [options.minElevDeg=10] Elevation that defines rise and set (default 10°)
+ * @param {number} [options.horizonHours=24] Maximum search horizon in hours (default 24h)
+ * @param {number} [options.coarseStepSec=20] Coarse scan step in seconds (default 20s)
+ * @param {number} [options.fineStepSec=5] Fine transit step in seconds (default 5s)
+ * @param {boolean} [options.requireVisible=false] Skip passes that are not naked-eye visible
+ * @param {number} [options.maxSunElevDeg=-6] Ceiling for observer darkness (default -6° for civil twilight)
+ * @returns {{
+ *   riseMs: number,
+ *   setMs: number,
+ *   maxElevDeg: number,
+ *   maxElevMs: number,
+ *   riseAzDeg: number,
+ *   visible: boolean,
+ *   sunlit: boolean,
+ *   observerDark: boolean
+ * } | null} `visible` is true when any sampled part of the pass is sunlit under a
+ *   dark sky; `sunlit` and `observerDark` describe the peak only.
+ */
+export function findNextSatellitePass({
+  satrec,
+  latDeg,
+  lonDeg,
+  fromMs,
+  minElevDeg = 10,
+  horizonHours = 24,
+  coarseStepSec = 20,
+  fineStepSec = 5,
+  requireVisible = false,
+  maxSunElevDeg = -6,
+}) {
+  if (!satrec || ![latDeg, lonDeg, fromMs, minElevDeg, horizonHours, coarseStepSec, fineStepSec, maxSunElevDeg].every(Number.isFinite)
+    || Math.abs(latDeg) > 90 || Math.abs(lonDeg) > 180 || Math.abs(fromMs) > 8.63e15
+    || minElevDeg < 0 || minElevDeg > 90 || horizonHours <= 0 || horizonHours > 72
+    || coarseStepSec < 1 || coarseStepSec > 120 || fineStepSec < 0.1 || fineStepSec > 60) return null;
+  const horizonMs = fromMs + horizonHours * 3600_000;
+  const coarseMs = coarseStepSec * 1000;
+  const fineMs = Math.max(1000, (fineStepSec || 5) * 1000);
+
+  let searchCursor = fromMs;
+
+  while (searchCursor <= horizonMs) {
+    let tPrev = searchCursor;
+    let tHit = null;
+    let elevPrev = lookAnglesAt(satrec, tPrev, latDeg, lonDeg)?.elevDeg ?? -90;
+
+    // Check if we start already inside a pass
+    if (elevPrev >= minElevDeg) {
+      tHit = searchCursor;
+    } else {
+      // Coarse forward scan
+      for (let t = searchCursor + coarseMs; t <= horizonMs; t += coarseMs) {
+        const look = lookAnglesAt(satrec, t, latDeg, lonDeg);
+        const elev = look ? look.elevDeg : -90;
+        if (elev >= minElevDeg) {
+          tHit = t;
+          break;
+        }
+        tPrev = t;
+        elevPrev = elev;
+      }
+    }
+
+    if (tHit == null) return null; // No pass found in remaining horizon
+
+    // Resolve rise
+    let riseMs;
+    if (tHit === searchCursor && elevPrev >= minElevDeg) {
+      riseMs = searchCursor;
+    } else {
+      riseMs = _bisectBoundary(
+        satrec,
+        latDeg,
+        lonDeg,
+        tPrev,
+        tHit,
+        minElevDeg,
+        true,
+      );
+    }
+
+    // Transit the pass to track peak and bracket set between tStepPrev and tCurr
+    let tCurr = riseMs;
+    let maxElevDeg = -90;
+    let maxElevMs = riseMs;
+    let tStepPrev = riseMs;
+    let setMs = null;
+    let hasVisibleSegment = false;
+    let elevPrevSample = null;
+    let peakLeftElev = null;
+    let peakRightElev = null;
+    let peakStepMs = fineMs;
+
+    // Visibility at an instant inside [riseMs, setMs], so no elevation gate
+    const checkVisibility = (tMs) => {
+      if (hasVisibleSegment) return;
+      const look = lookAnglesAt(satrec, tMs, latDeg, lonDeg);
+      if (!look) return;
+      const sunlit = isSatelliteSunlit(look.satECI, tMs);
+      const dark = isObserverDark(latDeg, lonDeg, tMs, maxSunElevDeg);
+      if (sunlit && dark) {
+        hasVisibleSegment = true;
+      }
+    };
+
+    // Check visibility at rise
+    checkVisibility(riseMs);
+
+    // Transit forward in fine steps
+    while (tCurr <= horizonMs) {
+      const look = lookAnglesAt(satrec, tCurr, latDeg, lonDeg);
+      const elev = look ? look.elevDeg : -90;
+
+      if (look) {
+        const sunlit = isSatelliteSunlit(look.satECI, tCurr);
+        const dark = isObserverDark(latDeg, lonDeg, tCurr, maxSunElevDeg);
+        if (sunlit && dark && elev >= minElevDeg) {
+          hasVisibleSegment = true;
+        }
+      }
+
+      if (elev < minElevDeg && tCurr > riseMs) {
+        if (peakRightElev == null && tCurr > maxElevMs) {
+          peakRightElev = elev;
+          peakStepMs = tCurr - maxElevMs;
+        }
+        // Below threshold: bisect set between tStepPrev and tCurr
+        setMs = _bisectBoundary(
+          satrec,
+          latDeg,
+          lonDeg,
+          tStepPrev,
+          tCurr,
+          minElevDeg,
+          false,
+        );
+        break;
+      }
+
+      if (elev > maxElevDeg) {
+        maxElevDeg = elev;
+        maxElevMs = tCurr;
+        peakLeftElev = elevPrevSample;
+        peakRightElev = null;
+      } else if (peakRightElev == null && tCurr > maxElevMs) {
+        peakRightElev = elev;
+        peakStepMs = tCurr - maxElevMs;
+      }
+
+      elevPrevSample = elev;
+      tStepPrev = tCurr;
+      tCurr += fineMs;
+    }
+
+    if (setMs == null) {
+      setMs = Math.min(horizonMs, tCurr);
+    }
+
+    // Refine peak culmination via 3-point parabolic interpolation around maxElevMs
+    // using adjacent transit samples (zero additional SGP4 calls)
+    let refinedPeakMs = maxElevMs;
+    let refinedPeakElevDeg = maxElevDeg;
+    if (
+      Number.isFinite(peakLeftElev) &&
+      Number.isFinite(peakRightElev) &&
+      peakLeftElev - 2 * maxElevDeg + peakRightElev < 0 // concave down
+    ) {
+      const denom = 2 * (peakLeftElev - 2 * maxElevDeg + peakRightElev);
+      const shift = ((peakLeftElev - peakRightElev) / denom) * peakStepMs;
+      if (Math.abs(shift) < peakStepMs) {
+        refinedPeakMs = Math.max(riseMs, Math.min(setMs, Math.round(maxElevMs + shift)));
+        refinedPeakElevDeg =
+          Math.min(90, maxElevDeg -
+          Math.pow(peakLeftElev - peakRightElev, 2) /
+            (8 * (peakLeftElev - 2 * maxElevDeg + peakRightElev)));
+      }
+    }
+
+    // Sun and sky state at the peak
+    const peakLook = lookAnglesAt(satrec, refinedPeakMs, latDeg, lonDeg);
+    const peakSunlit = peakLook
+      ? isSatelliteSunlit(peakLook.satECI, refinedPeakMs)
+      : false;
+    const peakDark = isObserverDark(
+      latDeg,
+      lonDeg,
+      refinedPeakMs,
+      maxSunElevDeg,
+    );
+
+    // Check visibility at the refined peak and at set
+    checkVisibility(refinedPeakMs);
+    checkVisibility(setMs);
+
+    if (!hasVisibleSegment) {
+      hasVisibleSegment = hasVisibleInterval(riseMs, setMs, (time) => {
+        const look = lookAnglesAt(satrec, time, latDeg, lonDeg);
+        return [Boolean(look && isSatelliteSunlit(look.satECI, time)), isObserverDark(latDeg, lonDeg, time, maxSunElevDeg)];
+      });
+    }
+
+    // Skip passes with no visible part when requireVisible is set
+    if (requireVisible && !hasVisibleSegment) {
+      searchCursor = Math.max(setMs + 1000, searchCursor + coarseMs);
+      continue;
+    }
+
+    const riseLook = lookAnglesAt(satrec, riseMs, latDeg, lonDeg);
+    return {
+      riseMs,
+      setMs,
+      maxElevDeg: refinedPeakElevDeg,
+      maxElevMs: refinedPeakMs,
+      riseAzDeg: riseLook ? riseLook.azDeg : 0,
+      visible: hasVisibleSegment,
+      sunlit: peakSunlit,
+      observerDark: peakDark,
+    };
+  }
+
+  return null;
+}
+
+
+/** Resolve shadow/twilight crossings to 20ms so short overlap is not lost between samples. */
+export function hasVisibleInterval(start, end, conditionsAt) {
+  const visible = (states) => states[0] && states[1];
+  for (let left = start; left < end; left += 5000) {
+    const right = Math.min(end, left + 5000);
+    const a = conditionsAt(left);
+    const b = conditionsAt(right);
+    if (visible(a) || visible(b)) return true;
+    const cuts = [left, right];
+    for (let condition = 0; condition < 2; condition++) {
+      if (a[condition] === b[condition]) continue;
+      let low = left;
+      let high = right;
+      while (high - low > 20) {
+        const mid = (low + high) / 2;
+        if (conditionsAt(mid)[condition] === a[condition]) low = mid;
+        else high = mid;
+      }
+      cuts.push((low + high) / 2);
+    }
+    cuts.sort((a, b) => a - b);
+    for (let i = 1; i < cuts.length; i++) {
+      if (visible(conditionsAt((cuts[i - 1] + cuts[i]) / 2))) return true;
+    }
+  }
+  return false;
+}