Add animated package along trajectory with Kepler propagation

- Trajectory sampling: propagate Lambert transfer orbits using universal
  variable Kepler propagation, sample 60 points per leg
- Smooth curved trajectory lines on 2D canvas (real orbital arcs, not
  straight lines)
- Animated green package dot with radial glow traveling along trajectory
- Route progress slider (0-100%) with play/pause animation toggle
- Auto-animating ~20-second cycle through the full route
- WASM API: sample_route_trajectory() with per-leg NaN separators

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

crates/mass-driver-wasm/src/api.rs

@@ -102,6 +102,76 @@ pub fn get_station_names(stations_json: &str) -> String {
     serde_json::to_string(&names).unwrap()
 }
 
+/// Sample trajectory points along a route for visualization.
+/// Takes the stations JSON, route JSON (from compute_route), and current JD.
+/// Returns flat Float64Array of [x,y,z, ...] in AU for all legs concatenated,
+/// with a separator of NaN,NaN,NaN between legs.
+#[wasm_bindgen]
+pub fn sample_route_trajectory(
+    stations_json: &str,
+    route_json: &str,
+    jd: f64,
+    samples_per_leg: usize,
+) -> Vec<f64> {
+    use mass_driver_core::route::RouteResult;
+    use orbital_mechanics::lambert;
+
+    let stations: Vec<station::Station> = match serde_json::from_str(stations_json) {
+        Ok(s) => s,
+        Err(_) => return vec![],
+    };
+    let route: RouteResult = match serde_json::from_str(route_json) {
+        Ok(r) => r,
+        Err(_) => return vec![],
+    };
+
+    let all_bodies = bodies::all_bodies();
+    let mut result = Vec::new();
+
+    for (i, leg) in route.legs.iter().enumerate() {
+        if i > 0 {
+            // Separator between legs
+            result.push(f64::NAN);
+            result.push(f64::NAN);
+            result.push(f64::NAN);
+        }
+
+        // Get station positions at departure/arrival times
+        let dep_jd = jd + (leg.departure_week as f64) * 7.0;
+        let arr_jd = jd + (leg.arrival_week as f64) * 7.0;
+        let tof_seconds = (arr_jd - dep_jd) * orbital_mechanics::constants::SECONDS_PER_DAY;
+
+        let r1_au = station::station_position(
+            &stations[leg.from_station], &all_bodies, dep_jd,
+        );
+        let r2_au = station::station_position(
+            &stations[leg.to_station], &all_bodies, arr_jd,
+        );
+
+        // Convert to km for Lambert solver
+        let au_km = orbital_mechanics::constants::AU_KM;
+        let r1_km = r1_au * au_km;
+        let r2_km = r2_au * au_km;
+
+        let points = lambert::sample_trajectory(
+            r1_km, r2_km, tof_seconds,
+            orbital_mechanics::constants::MU_SUN,
+            samples_per_leg,
+        );
+
+        // Convert back to AU
+        for chunk in points.chunks(3) {
+            if chunk.len() == 3 {
+                result.push(chunk[0] / au_km);
+                result.push(chunk[1] / au_km);
+                result.push(chunk[2] / au_km);
+            }
+        }
+    }
+
+    result
+}
+
 /// Compute transfer matrix and find optimal route between two stations.
 ///
 /// This is the main computation function. It:

crates/orbital-mechanics/src/lagrange.rs

@@ -72,7 +72,7 @@ pub fn lagrange_point_position(
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::bodies::all_bodies;
+    use crate::bodies::{self, all_bodies};
 
     #[test]
     fn test_l1_between_sun_and_planet() {

crates/orbital-mechanics/src/lambert.rs

@@ -214,6 +214,114 @@ pub fn transfer_delta_v(
     Some((dv1, dv2))
 }
 
+/// Sample points along a Lambert transfer trajectory.
+/// Returns positions at evenly spaced times between departure and arrival.
+///
+/// # Arguments
+/// * `r1`, `r2` - Start/end positions (km)
+/// * `tof` - Time of flight (seconds)
+/// * `mu` - Gravitational parameter (km³/s²)
+/// * `num_samples` - Number of sample points (including start and end)
+///
+/// Returns flat [x,y,z, x,y,z, ...] in km, or empty if Lambert fails.
+pub fn sample_trajectory(
+    r1: Vector3<f64>,
+    r2: Vector3<f64>,
+    tof: f64,
+    mu: f64,
+    num_samples: usize,
+) -> Vec<f64> {
+    let sol = match solve_lambert(r1, r2, tof, mu, true) {
+        Some(s) => s,
+        None => return vec![],
+    };
+
+    if num_samples < 2 {
+        return vec![r1.x, r1.y, r1.z];
+    }
+
+    let mut result = Vec::with_capacity(num_samples * 3);
+
+    for i in 0..num_samples {
+        let t_frac = i as f64 / (num_samples - 1) as f64;
+        let dt = t_frac * tof;
+
+        if i == 0 {
+            result.push(r1.x);
+            result.push(r1.y);
+            result.push(r1.z);
+            continue;
+        }
+        if i == num_samples - 1 {
+            result.push(r2.x);
+            result.push(r2.y);
+            result.push(r2.z);
+            continue;
+        }
+
+        // Propagate from r1 with v1 using the universal variable
+        let pos = propagate_kepler(r1, sol.v1, dt, mu);
+        result.push(pos.x);
+        result.push(pos.y);
+        result.push(pos.z);
+    }
+
+    result
+}
+
+/// Propagate a Keplerian orbit forward by time dt from given state (r, v).
+fn propagate_kepler(r0: Vector3<f64>, v0: Vector3<f64>, dt: f64, mu: f64) -> Vector3<f64> {
+    let r0_mag = r0.norm();
+    let v0_mag = v0.norm();
+    let vr0 = r0.dot(&v0) / r0_mag; // radial velocity
+    let alpha = 2.0 / r0_mag - v0_mag * v0_mag / mu; // 1/a
+
+    // Solve for universal variable chi using Newton-Raphson
+    let sqrt_mu = mu.sqrt();
+    // Initial guess for universal variable
+    let mut chi = if alpha > 1e-10 {
+        sqrt_mu * dt * alpha
+    } else {
+        sqrt_mu * dt / r0_mag
+    };
+
+    for _ in 0..50 {
+        let chi2 = chi * chi;
+        let z = alpha * chi2;
+        let (c2, c3) = super::lambert::stumpff(z);
+
+        let _r = chi2 * c2 + vr0 / sqrt_mu * chi * (1.0 - z * c3) + r0_mag * (1.0 - z * c2);
+
+        // f(chi) = r0 * chi * (1 - z*c3) + vr0/sqrt_mu * chi^2 * c2 + r0_mag * chi^3 * c3 - sqrt_mu * dt
+        let f_chi = vr0 / sqrt_mu * chi2 * c2
+            + (1.0 - r0_mag * alpha) * chi * chi2 * c3
+            + r0_mag * chi
+            - sqrt_mu * dt;
+
+        let f_prime = chi2 * c2 + vr0 / sqrt_mu * chi * (1.0 - z * c3) + r0_mag * (1.0 - z * c2);
+
+        if f_prime.abs() < 1e-20 {
+            break;
+        }
+
+        chi -= f_chi / f_prime;
+
+        if (f_chi / f_prime).abs() < 1e-10 {
+            break;
+        }
+    }
+
+    let chi2 = chi * chi;
+    let z = alpha * chi2;
+    let (c2, c3) = stumpff(z);
+
+    // Lagrange coefficients
+    let f = 1.0 - chi2 / r0_mag * c2;
+    let g = dt - chi * chi2 / sqrt_mu * c3;
+
+    f * r0 + g * v0
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;