Phase 0: monorepo skeleton (hub, live-map, api, packages, infra, CI)

packages/orbital-math/src/index.ts

@@ -0,0 +1,9 @@
+export { solveKepler } from './kepler.js';
+export {
+  meanMotion,
+  propagateToEpoch,
+  positionAt,
+  sampleOrbit,
+} from './propagate.js';
+export { shadowFraction } from './occultation.js';
+export { phaseAngle, hohmannDeltaV, findTransferWindows } from './transfer.js';

packages/orbital-math/src/kepler.ts

@@ -0,0 +1,28 @@
+/**
+ * Solve Kepler's equation:  M = E - e * sin(E)
+ * for the eccentric anomaly E, given mean anomaly M and eccentricity e.
+ *
+ * Uses Newton-Raphson iteration with a sane initial guess.
+ * Converges in ~5 iterations for any reasonable e (< 0.9).
+ * For near-parabolic / hyperbolic orbits, use a different solver.
+ */
+export function solveKepler(meanAnomaly: number, eccentricity: number): number {
+  // Normalize M to [-π, π] for faster convergence.
+  const TWO_PI = Math.PI * 2;
+  let M = meanAnomaly % TWO_PI;
+  if (M > Math.PI) M -= TWO_PI;
+  if (M < -Math.PI) M += TWO_PI;
+
+  // Initial guess: E₀ = M + e·sin(M) is a good first order approximation.
+  let E = M + eccentricity * Math.sin(M);
+
+  for (let i = 0; i < 30; i++) {
+    const f = E - eccentricity * Math.sin(E) - M;
+    const fp = 1 - eccentricity * Math.cos(E);
+    const dE = f / fp;
+    E -= dE;
+    if (Math.abs(dE) < 1e-12) break;
+  }
+
+  return E;
+}

packages/orbital-math/src/occultation.ts

@@ -0,0 +1,53 @@
+/**
+ * Geometric occultation: given a position (relative to a body's center)
+ * and the radii of the occluder (R1) and the body the observer is on (R2),
+ * is the observer currently in shadow?
+ *
+ * Used for both:
+ *   - "is this vessel in the planet's shadow?" (R1 = planet radius, R2 ≈ 0)
+ *   - "is this ground station blocked by the local terrain?" (R1 = planet, R2 = earth station)
+ *
+ * Returns the fraction (0..1) of the line of sight to the sun that is
+ * occluded. 0 = full sun, 1 = total eclipse.
+ *
+ * Note: the canonical way to do this is to compute the half-angle between
+ * the sun and the occluding body as seen by the observer. We treat the
+ * sun as effectively at infinity (parallel rays) which is fine for KSP
+ * since Kerbol is the system root and we're never going to need parallax
+ * precision at this scale.
+ */
+export function shadowFraction(
+  observerToSun: { x: number; y: number; z: number },
+  occluderToObserver: { x: number; y: number; z: number },
+  occluderRadius: number,
+): number {
+  // Vector from observer to sun, normalized
+  const sunDist = Math.hypot(observerToSun.x, observerToSun.y, observerToSun.z);
+  if (sunDist === 0) return 0;
+  const sx = observerToSun.x / sunDist;
+  const sy = observerToSun.y / sunDist;
+  const sz = observerToSun.z / sunDist;
+
+  // Project occluder center onto the sun-direction line
+  const proj = occluderToObserver.x * sx + occluderToObserver.y * sy + occluderToObserver.z * sz;
+  if (proj >= 0) {
+    // Occluder is behind the observer relative to the sun → no eclipse
+    return 0;
+  }
+  // Perpendicular distance from occluder center to sun ray
+  const px = occluderToObserver.x - proj * sx;
+  const py = occluderToObserver.y - proj * sy;
+  const pz = occluderToObserver.z - proj * sz;
+  const perpDist = Math.hypot(px, py, pz);
+
+  if (perpDist >= occluderRadius) return 0;
+  // Approximate chord length through the occluder disc
+  const halfChord = Math.sqrt(occluderRadius * occluderRadius - perpDist * perpDist);
+  // Approximate the angular size of the sun as seen from the occluder
+  // vs the angular size of the occluder; we use 1.0 for the sun
+  // (i.e. effectively point source) — good enough for visualization.
+  // For a "fraction in shadow" treat the occluder disc as fully shadowing
+  // when perpDist + halfChord reaches the observer; that simplifies to
+  // perpDist < occluderRadius which we already check.
+  return Math.min(1, 1 - perpDist / occluderRadius);
+}

packages/orbital-math/src/propagate.ts

@@ -0,0 +1,141 @@
+import type { KeplerianElements, CelestialBody } from '@kerbal-rt/shared-types';
+import { solveKepler } from './kepler.js';
+
+const TWO_PI = Math.PI * 2;
+
+/**
+ * Compute the mean motion n (rad/s) from semi-major axis and the
+ * gravitational parameter μ of the central body.
+ *
+ *   n = sqrt(μ / a^3)
+ */
+export function meanMotion(semiMajorAxis: number, mu: number): number {
+  return Math.sqrt(mu / Math.pow(semiMajorAxis, 3));
+}
+
+/**
+ * Propagate Keplerian elements forward in time to a new epoch.
+ *
+ * Only `meanAnomalyAtEpoch` and `epoch` change. (For a fully accurate
+ * J2-perturbed propagation you'd also adjust RAAN and argPe due to
+ * precession, but for visualization purposes this is plenty.)
+ */
+export function propagateToEpoch(
+  elements: KeplerianElements,
+  mu: number,
+  newEpoch: number,
+): KeplerianElements {
+  const dt = newEpoch - elements.epoch;
+  const n = meanMotion(elements.semiMajorAxis, mu);
+  return {
+    ...elements,
+    epoch: newEpoch,
+    meanAnomalyAtEpoch: elements.meanAnomalyAtEpoch + n * dt,
+  };
+}
+
+/**
+ * Propagate an orbit to a target UT and return the position in the
+ * parent body's inertial frame.
+ *
+ * For an ellipse (e < 1):
+ *   1. Compute mean anomaly at target time
+ *   2. Solve Kepler for eccentric anomaly E
+ *   3. Position in perifocal frame: (a(cosE - e), a√(1-e²) sinE, 0)
+ *   4. Rotate by ω, i, Ω to get the inertial frame
+ *
+ * Returns the cartesian position vector (m).
+ */
+export function positionAt(
+  elements: KeplerianElements,
+  mu: number,
+  ut: number,
+): { x: number; y: number; z: number } {
+  const e = elements.eccentricity;
+  const a = elements.semiMajorAxis;
+  const i = elements.inclination;
+  const O = elements.longitudeOfAscendingNode; // RAAN
+  const w = elements.argumentOfPeriapsis;
+
+  // Mean anomaly at the requested time
+  const n = meanMotion(a, mu);
+  const M = elements.meanAnomalyAtEpoch + n * (ut - elements.epoch);
+
+  // Eccentric anomaly
+  const E = solveKepler(M, e);
+
+  // True anomaly
+  //   cos ν = (cos E - e) / (1 - e cos E)
+  //   sin ν = (√(1-e²) sin E) / (1 - e cos E)
+  const cosE = Math.cos(E);
+  const sinE = Math.sin(E);
+  const denom = 1 - e * cosE;
+  const cosNu = (cosE - e) / denom;
+  const sinNu = (Math.sqrt(1 - e * e) * sinE) / denom;
+
+  // Distance
+  const r = a * (1 - e * cosE);
+
+  // Position in perifocal frame (P, Q, W)
+  const xP = r * cosNu;
+  const yQ = r * sinNu;
+
+  // Rotation matrix from perifocal to inertial: Rz(-Ω) · Rx(-i) · Rz(-ω)
+  // Applied to (xP, yQ, 0):
+  const cosO = Math.cos(O);
+  const sinO = Math.sin(O);
+  const cosi = Math.cos(i);
+  const sini = Math.sin(i);
+  const cosw = Math.cos(w);
+  const sinw = Math.sin(w);
+
+  const x =
+    (cosO * cosw - sinO * sinw * cosi) * xP + (-cosO * sinw - sinO * cosw * cosi) * yQ;
+  const y =
+    (sinO * cosw + cosO * sinw * cosi) * xP + (-sinO * sinw + cosO * cosw * cosi) * yQ;
+  const z = sinw * sini * xP + cosw * sini * yQ;
+
+  return { x, y, z };
+}
+
+/**
+ * Sample a full orbit as a list of cartesian points. Useful for drawing
+ * the orbit line. Returns `steps` evenly-spaced points in true anomaly
+ * around the conic.
+ */
+export function sampleOrbit(
+  elements: KeplerianElements,
+  mu: number,
+  steps: number = 128,
+): { x: number; y: number; z: number }[] {
+  const points: { x: number; y: number; z: number }[] = [];
+  const e = elements.eccentricity;
+  const a = elements.semiMajorAxis;
+  const i = elements.inclination;
+  const O = elements.longitudeOfAscendingNode;
+  const w = elements.argumentOfPeriapsis;
+
+  const cosO = Math.cos(O);
+  const sinO = Math.sin(O);
+  const cosi = Math.cos(i);
+  const sini = Math.sin(i);
+  const cosw = Math.cos(w);
+  const sinw = Math.sin(w);
+
+  for (let k = 0; k < steps; k++) {
+    // True anomaly uniformly from 0 to 2π
+    const nu = (k / steps) * TWO_PI;
+    const cosNu = Math.cos(nu);
+    const sinNu = Math.sin(nu);
+    const r = (a * (1 - e * e)) / (1 + e * cosNu);
+    const xP = r * cosNu;
+    const yQ = r * sinNu;
+    const x =
+      (cosO * cosw - sinO * sinw * cosi) * xP + (-cosO * sinw - sinO * cosw * cosi) * yQ;
+    const y =
+      (sinO * cosw + cosO * sinw * cosi) * xP + (-sinO * sinw + cosO * cosw * cosi) * yQ;
+    const z = sinw * sini * xP + cosw * sini * yQ;
+    points.push({ x, y, z });
+  }
+  return points;
+}

packages/orbital-math/src/transfer.ts

@@ -0,0 +1,95 @@
+import type { KeplerianElements, CelestialBody } from '@kerbal-rt/shared-types';
+import { meanMotion } from './propagate.js';
+
+const TWO_PI = Math.PI * 2;
+
+/**
+ * Compute the phase angle between two orbiting bodies at a given UT.
+ * Phase angle is the angle at the central body between the two bodies,
+ * measured in the direction of motion.
+ */
+export function phaseAngle(
+  a: KeplerianElements,
+  aMu: number,
+  b: KeplerianElements,
+  bMu: number,
+  ut: number,
+): number {
+  const nA = meanMotion(a.semiMajorAxis, aMu);
+  const nB = meanMotion(b.semiMajorAxis, bMu);
+  const MA = a.meanAnomalyAtEpoch + nA * (ut - a.epoch);
+  const MB = b.meanAnomalyAtEpoch + nB * (ut - b.epoch);
+  let phase = MB - MA;
+  // Normalize to [0, 2π)
+  phase = phase % TWO_PI;
+  if (phase < 0) phase += TWO_PI;
+  return phase;
+}
+
+/**
+ * Very rough Hohmann transfer Δv estimate between two coplanar circular
+ * orbits. Good enough for a "transfer window" calculator; for precise
+ * numbers you'd solve Lambert's problem.
+ *
+ * Returns the total Δv (m/s) for the two-burn transfer.
+ */
+export function hohmannDeltaV(
+  r1: number,
+  r2: number,
+  mu: number,
+): { dv1: number; dv2: number; total: number; transferTime: number } {
+  // Semi-major axis of transfer ellipse
+  const aTransfer = (r1 + r2) / 2;
+  // Circular velocities
+  const v1 = Math.sqrt(mu / r1);
+  const v2 = Math.sqrt(mu / r2);
+  // Velocities on transfer ellipse at periapsis (r1) and apoapsis (r2)
+  const vt1 = Math.sqrt(mu * (2 / r1 - 1 / aTransfer));
+  const vt2 = Math.sqrt(mu * (2 / r2 - 1 / aTransfer));
+  const dv1 = Math.abs(vt1 - v1);
+  const dv2 = Math.abs(v2 - vt2);
+  const transferTime = Math.PI * Math.sqrt(Math.pow(aTransfer, 3) / mu);
+  return { dv1, dv2, total: dv1 + dv2, transferTime };
+}
+
+/**
+ * Find the next N transfer windows from `from` body to `to` body,
+ * defined as times when the phase angle is within ±tolerance of the
+ * ideal Hohmann transfer phase.
+ *
+ * Ideal phase = π · (1 - (1/2) · ((r2/r1)^(3/2) + 1)^(-2/3))
+ * (this is the classic approximation for the case r2 > r1)
+ *
+ * We just scan forward from `ut` until we've found `n` windows.
+ */
+export function findTransferWindows(
+  from: { elements: KeplerianElements; mu: number },
+  to: { elements: KeplerianElements; mu: number },
+  ut: number,
+  options: { count?: number; toleranceRad?: number; maxSearchTime?: number } = {},
+): { ut: number; phaseAngle: number }[] {
+  const { count = 3, toleranceRad = 0.15, maxSearchTime = 5 * 365 * 24 * 3600 } = options;
+  const r1 = from.elements.semiMajorAxis;
+  const r2 = to.elements.semiMajorAxis;
+  const ratio = Math.pow(r2 / r1, 1.5);
+  const idealPhase = Math.PI * (1 - Math.pow(1 + ratio, -2 / 3) * 0.5);
+  // For inner→outer transfers; for outer→inner the phase wraps differently.
+  // For simplicity, accept either ±(2π - idealPhase) too.
+  const phases = [idealPhase, TWO_PI - idealPhase];
+
+  const results: { ut: number; phaseAngle: number }[] = [];
+  const stepSec = 3600; // 1h steps for the coarse scan
+  const nA = meanMotion(from.elements.semiMajorAxis, from.mu);
+
+  for (let t = ut; t < ut + maxSearchTime && results.length < count; t += stepSec) {
+    const phase = phaseAngle(from.elements, from.mu, to.elements, to.mu, t);
+    for (const target of phases) {
+      const diff = Math.abs(((phase - target + Math.PI) % TWO_PI) - Math.PI);
+      if (diff < toleranceRad) {
+        results.push({ ut: t, phaseAngle: phase });
+        break;
+      }
+    }
+  }
+  return results;
+}