Add Lambert solver, transfer matrix, Dijkstra routing, and route planner UI

- Lambert's problem solver using universal variable method with bisection
  (handles elliptic, parabolic, hyperbolic transfers + anti-podal cases)
- Transfer matrix: precompute pairwise station transfers over time window
  using Lambert solver with configurable launch velocity
- Dijkstra routing on time-expanded graph (station × week nodes, transfer
  + wait edges) to find minimum-time routes
- Route Planner UI: from/to station dropdowns, search window selector
  (1-10 years), "Find Optimal Route" button with results card
- Route visualization: orange dashed trajectory lines with arrow heads
  and leg numbers on the 2D canvas
- Tested: Mercury L1 → Jupiter L1 computes 6-month direct transfer at
  30 km/s — physically reasonable

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

crates/orbital-mechanics/src/lambert.rs

@@ -0,0 +1,300 @@
+//! Lambert's problem solver using universal variable method.
+//!
+//! Given two position vectors r1 and r2, and a time of flight tof,
+//! find the orbit connecting them. Returns departure and arrival velocities.
+//!
+//! Uses the algorithm from Curtis, "Orbital Mechanics for Engineering Students",
+//! with bisection for robust convergence.
+
+use nalgebra::Vector3;
+
+/// Solution to Lambert's problem.
+#[derive(Debug, Clone)]
+pub struct LambertSolution {
+    /// Departure velocity (km/s)
+    pub v1: Vector3<f64>,
+    /// Arrival velocity (km/s)
+    pub v2: Vector3<f64>,
+}
+
+/// Solve Lambert's problem.
+///
+/// # Arguments
+/// * `r1` - Initial position vector (km)
+/// * `r2` - Final position vector (km)
+/// * `tof` - Time of flight (seconds)
+/// * `mu` - Gravitational parameter (km³/s²)
+/// * `prograde` - If true, use prograde (short-way) transfer
+///
+/// # Returns
+/// `Some(LambertSolution)` if a solution exists, `None` otherwise.
+pub fn solve_lambert(
+    r1: Vector3<f64>,
+    r2: Vector3<f64>,
+    tof: f64,
+    mu: f64,
+    prograde: bool,
+) -> Option<LambertSolution> {
+    if tof <= 0.0 || mu <= 0.0 {
+        return None;
+    }
+
+    let r1_norm = r1.norm();
+    let r2_norm = r2.norm();
+
+    if r1_norm < 1e-6 || r2_norm < 1e-6 {
+        return None;
+    }
+
+    // Handle near-180° transfers by adding a tiny out-of-plane perturbation
+    let cos_angle = r1.dot(&r2) / (r1_norm * r2_norm);
+    let r2 = if cos_angle < -0.99999 {
+        // Near anti-podal: perturb slightly to break degeneracy
+        Vector3::new(r2.x, r2.y, r2.z + r2_norm * 1e-6)
+    } else {
+        r2
+    };
+    let r2_norm = r2.norm();
+
+    // Change in true anomaly
+    let cos_dnu = r1.dot(&r2) / (r1_norm * r2_norm);
+    let cos_dnu = cos_dnu.clamp(-1.0, 1.0);
+
+    let cross = r1.cross(&r2);
+
+    let sin_dnu = if prograde {
+        if cross.z >= 0.0 { (1.0 - cos_dnu * cos_dnu).sqrt() }
+        else { -(1.0 - cos_dnu * cos_dnu).sqrt() }
+    } else {
+        if cross.z < 0.0 { (1.0 - cos_dnu * cos_dnu).sqrt() }
+        else { -(1.0 - cos_dnu * cos_dnu).sqrt() }
+    };
+
+    // Parameter A (from Curtis Eq. 5.35)
+    let a = sin_dnu * (r1_norm * r2_norm / (1.0 - cos_dnu)).sqrt();
+
+    if a.abs() < 1e-10 {
+        return None;
+    }
+
+    // Function to compute TOF for a given z using Stumpff functions
+    let tof_from_z = |z: f64| -> Option<f64> {
+        let (c2, c3) = stumpff(z);
+        if c2.abs() < 1e-20 {
+            return None;
+        }
+
+        let y = r1_norm + r2_norm + a * (z * c3 - 1.0) / c2.sqrt();
+        if y < 0.0 {
+            return None;
+        }
+
+        let chi = (y / c2).sqrt();
+        let t = (chi * chi * chi * c3 + a * y.sqrt()) / mu.sqrt();
+        Some(t)
+    };
+
+    // Find z using bisection + Newton hybrid
+    // For elliptic orbits: z > 0
+    // For hyperbolic orbits: z < 0
+
+    // Find a bracket [z_low, z_high] containing the solution.
+    // z_low should give TOF < target, z_high should give TOF > target.
+    let mut z_low = -4.0 * std::f64::consts::PI * std::f64::consts::PI;
+    let mut z_high = 4.0 * std::f64::consts::PI * std::f64::consts::PI;
+
+    // Ensure z_high gives TOF >= target (keep doubling if needed)
+    for _ in 0..50 {
+        match tof_from_z(z_high) {
+            Some(t) if t >= tof => break,
+            _ => z_high *= 2.0,
+        }
+        if z_high > 1e10 { return None; }
+    }
+
+    // Ensure z_low gives TOF < target (keep decreasing if needed)
+    for _ in 0..50 {
+        match tof_from_z(z_low) {
+            Some(t) if t < tof => break,
+            None => break, // y < 0 means TOF undefined → z is too low, which is fine
+            _ => z_low -= 10.0,
+        }
+        if z_low < -1e10 { return None; }
+    }
+
+    // Bisection iteration
+    // Higher z → larger semi-major axis → longer TOF
+    let mut z = 0.0;
+    for _ in 0..200 {
+        z = (z_low + z_high) / 2.0;
+
+        let t = match tof_from_z(z) {
+            Some(t) => t,
+            None => {
+                // y < 0 means z is too low; increase lower bound
+                z_low = z;
+                continue;
+            }
+        };
+
+        if (t - tof).abs() < 1e-6 {
+            break; // Converged
+        }
+
+        if t < tof {
+            z_low = z; // Need more time → increase z
+        } else {
+            z_high = z; // Need less time → decrease z
+        }
+
+        if (z_high - z_low).abs() < 1e-12 {
+            break;
+        }
+    }
+
+    // Compute velocities from the converged z
+    let (c2, c3) = stumpff(z);
+    if c2.abs() < 1e-20 {
+        return None;
+    }
+
+    let y = r1_norm + r2_norm + a * (z * c3 - 1.0) / c2.sqrt();
+    if y < 0.0 {
+        return None;
+    }
+
+    // Lagrange coefficients
+    let f = 1.0 - y / r1_norm;
+    let g = a * (y / mu).sqrt();
+    let g_dot = 1.0 - y / r2_norm;
+
+    if g.abs() < 1e-20 {
+        return None;
+    }
+
+    let v1 = (r2 - f * r1) / g;
+    let v2 = (g_dot * r2 - r1) / g;
+
+    Some(LambertSolution { v1, v2 })
+}
+
+/// Stumpff functions c2(z) and c3(z).
+fn stumpff(z: f64) -> (f64, f64) {
+    if z > 1e-6 {
+        let sz = z.sqrt();
+        let c2 = (1.0 - sz.cos()) / z;
+        let c3 = (sz - sz.sin()) / (z * sz);
+        (c2, c3)
+    } else if z < -1e-6 {
+        let sz = (-z).sqrt();
+        let c2 = (sz.cosh() - 1.0) / (-z);
+        let c3 = (sz.sinh() - sz) / ((-z) * sz);
+        (c2, c3)
+    } else {
+        // Taylor series near z = 0
+        (1.0 / 2.0 - z / 24.0 + z * z / 720.0,
+         1.0 / 6.0 - z / 120.0 + z * z / 5040.0)
+    }
+}
+
+/// Compute the delta-v required for a Lambert transfer.
+/// Returns (dv_departure, dv_arrival) in km/s, where dv is relative to
+/// the velocity at each position.
+pub fn transfer_delta_v(
+    r1: Vector3<f64>,
+    v1_body: Vector3<f64>,
+    r2: Vector3<f64>,
+    v2_body: Vector3<f64>,
+    tof: f64,
+    mu: f64,
+) -> Option<(f64, f64)> {
+    let sol = solve_lambert(r1, r2, tof, mu, true)?;
+    let dv1 = (sol.v1 - v1_body).norm();
+    let dv2 = (sol.v2 - v2_body).norm();
+    Some((dv1, dv2))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::constants::*;
+
+    #[test]
+    fn test_hohmann_earth_mars() {
+        // Earth at 1 AU, Mars at 1.524 AU (opposite side)
+        let r1 = Vector3::new(AU_KM, 0.0, 0.0);
+        let r2 = Vector3::new(-1.524 * AU_KM, 0.0, 0.0);
+
+        // Hohmann transfer time ≈ 259 days
+        let tof = 259.0 * SECONDS_PER_DAY;
+
+        let result = solve_lambert(r1, r2, tof, MU_SUN, true);
+        assert!(result.is_some(), "Lambert solver should find a solution");
+
+        let sol = result.unwrap();
+        let v1_mag = sol.v1.norm();
+        // Departure velocity should be ~30-35 km/s (Earth orbital velocity + transfer excess)
+        assert!(
+            v1_mag > 20.0 && v1_mag < 45.0,
+            "Departure velocity {} km/s should be reasonable",
+            v1_mag
+        );
+    }
+
+    #[test]
+    fn test_short_transfer() {
+        // 30-degree transfer at 1 AU
+        let r1 = Vector3::new(AU_KM, 0.0, 0.0);
+        let angle = std::f64::consts::FRAC_PI_6;
+        let r2 = Vector3::new(AU_KM * angle.cos(), AU_KM * angle.sin(), 0.0);
+
+        let tof = 50.0 * SECONDS_PER_DAY;
+
+        let result = solve_lambert(r1, r2, tof, MU_SUN, true);
+        assert!(result.is_some(), "Should solve short transfer");
+
+        let sol = result.unwrap();
+        let v1_mag = sol.v1.norm();
+        assert!(
+            v1_mag > 10.0 && v1_mag < 60.0,
+            "Velocity {} km/s should be reasonable",
+            v1_mag
+        );
+    }
+
+    #[test]
+    fn test_circular_orbit_consistency() {
+        // If we solve Lambert for two points on a circular orbit with
+        // the correct TOF, the departure velocity should match orbital velocity
+        let r = AU_KM; // 1 AU
+        let v_circ = (MU_SUN / r).sqrt(); // ~29.8 km/s
+
+        let angle = std::f64::consts::FRAC_PI_4; // 45 degrees
+        let r1 = Vector3::new(r, 0.0, 0.0);
+        let r2 = Vector3::new(r * angle.cos(), r * angle.sin(), 0.0);
+
+        // Time for 45° of circular orbit
+        let period = 2.0 * std::f64::consts::PI * (r.powi(3) / MU_SUN).sqrt();
+        let tof = period * 45.0 / 360.0;
+
+        let result = solve_lambert(r1, r2, tof, MU_SUN, true);
+        assert!(result.is_some(), "Should solve circular orbit transfer");
+
+        let sol = result.unwrap();
+        let v1_mag = sol.v1.norm();
+        // Should be close to circular velocity
+        assert!(
+            (v1_mag - v_circ).abs() / v_circ < 0.05,
+            "Departure velocity {} km/s should match circular velocity {} km/s",
+            v1_mag, v_circ
+        );
+    }
+
+    #[test]
+    fn test_degenerate_zero_tof() {
+        let r1 = Vector3::new(AU_KM, 0.0, 0.0);
+        let r2 = Vector3::new(0.0, AU_KM, 0.0);
+        let result = solve_lambert(r1, r2, 0.0, MU_SUN, true);
+        assert!(result.is_none(), "Zero TOF should return None");
+    }
+}

crates/orbital-mechanics/src/lib.rs

@@ -2,6 +2,7 @@ pub mod bodies;
 pub mod constants;
 pub mod kepler;
 pub mod lagrange;
+pub mod lambert;
 pub mod orbits;
 
 pub use nalgebra::Vector3;