ssmichael1 · ssmichael1 · Mar 21, 2026 · Mar 21, 2026 · Mar 21, 2026 · Mar 21, 2026
diff --git a/.gitattributes b/.gitattributes
@@ -0,0 +1,2 @@
+*.ipynb filter=nbstripout
+*.ipynb diff=ipynb
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
@@ -46,16 +46,23 @@ jobs:
           path: astro-data
           key: astro-data-v1
 
+      - name: Install system dependencies for cartopy
+        run: sudo apt-get update && sudo apt-get install -y libgeos-dev libproj-dev proj-data
+
       - name: Install satkit
         run: pip install -e .
 
       - name: Install docs dependencies
         run: pip install -r docs/requirements.txt
 
+      - name: Generate static plots
+        env:
+          SATKIT_DATA: astro-data
+        run: python docs/examples/gen_plots.py
+
       - name: Build docs
         env:
           SATKIT_DATA: astro-data
-          PLOTLY_RENDERER: notebook
         run: mkdocs build
 
       - name: Upload pages artifact

diff --git a/.gitignore b/.gitignore
@@ -30,3 +30,6 @@ wheel
 
 # MkDocs build output
 site
+
+# Generated doc plots (built by docs/examples/gen_plots.py)
+docs/images/*.svg
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,6 +1,108 @@
 # Changelog
 
 
+## 0.14.1 - 2026-03-21
+
+### Lambert Targeting
+
+- Add Lambert's problem solver using Izzo's algorithm (2015) with Householder 4th-order iteration
+- Handles all orbit types (elliptic, parabolic, hyperbolic) and 180-degree transfers
+- Multi-revolution solution support
+- Rust API: `satkit::lambert::lambert()` in new `lambert` module
+- Python API: `satkit.lambert()` with full type stub and docstring
+- 9 Rust tests, 6 Python tests
+
+### Documentation Overhaul
+
+- Switch all tutorial plots from Plotly to matplotlib with SciencePlots
+- STIX serif fonts, SVG output, colorblind-friendly palette matching numeris docs
+- Match numeris site theme: blue grey header with navigation tabs
+- Muted steel blue link color (#4a7c96)
+- Static SVG plots (density, forces) generated at build time via `docs/examples/gen_plots.py`
+- Add Lambert Targeting tutorial with delta-v analysis, pork-chop plot, and orbit visualization
+- Add Lambert solver to User Guide and API Reference
+- Fix docstring formatting across all Python stubs: bullet style, Returns sections, clickable URLs, Notes admonitions
+- Add concrete defs for all `@typing.overload` functions so mkdocstrings renders them
+- Add `@typing.overload` signatures to `propresult.interp` for all call patterns
+- CI/CD: add cartopy/certifi system deps, generate plots before build, remove plotly
+
+### Internal
+
+- 114 Rust tests pass (9 new Lambert)
+- 58 Python tests pass (6 new Lambert)
+
+## 0.14.0 - 2026-03-20
+
+### Breaking: Replace nalgebra with numeris
+
+The `nalgebra` dependency has been replaced with `numeris` 0.5.6 for all linear algebra. The built-in ODE solver module (`src/ode/`) has been removed in favor of the ODE solvers provided by `numeris`. This is a **breaking change** for Rust API consumers; the Python API is unchanged.
+
+### Rust API Changes
+
+- **Math types** (`satkit::mathtypes`): All type aliases now point to `numeris` types instead of `nalgebra`. `Vector<N>`, `Matrix<M,N>`, `Quaternion`, and `DMatrix<T>` remain available with the same names.
+- **Vector construction**: `Vector3::new(x, y, z)` is replaced by `numeris::vector![x, y, z]`
+- **Matrix construction**: `Matrix3::new(a,b,c,d,e,f,g,h,i)` is replaced by `Matrix3::new([[a,b,c],[d,e,f],[g,h,i]])`
+- **Identity matrix**: `Matrix::identity()` is replaced by `Matrix::eye()`
+- **Quaternion axis rotations**: `Quaternion::from_axis_angle(&Vector3::z_axis(), θ)` is replaced by `Quaternion::rotz(θ)` (also `rotx`, `roty`)
+- **Quaternion vector rotation**: `q.transform_vector(&v)` is replaced by `q * v`
+- **Quaternion storage order**: Changed from nalgebra's `[x,y,z,w]` to numeris `[w,x,y,z]` (scalar-first). Component access uses `.w`, `.x`, `.y`, `.z` fields.
+- **Block extraction**: `m.fixed_view::<R,C>(i,j)` is replaced by `m.block::<R,C>(i,j)` (returns owned copy)
+- **Block insertion**: `m.fixed_view_mut::<R,C>(i,j).copy_from(&src)` is replaced by `m.set_block(i, j, &src)`
+- **Matrix inverse**: `m.try_inverse()` (returning `Option`) is replaced by `m.inverse()` (returning `Result`)
+- **Cholesky**: `m.cholesky()` now returns `Result<CholeskyDecomposition, LinalgError>` instead of `Option`
+
+### ODE Module Removed
+
+The `satkit::ode` module (7 adaptive solvers, Rosenbrock, ODEState trait, ~2,500 lines) has been removed. Orbit propagation now uses `numeris::ode` solvers directly. The same solver algorithms are available (RKF45, RKTS54, RKV65, RKV87, RKV98, RKV98NoInterp, RODAS4). The `PropSettings` API and `Integrator` enum are unchanged.
+
+### nalgebra Interoperability
+
+If you need nalgebra types for interoperability with other crates, enable the `nalgebra` feature on `numeris`:
+
+```toml
+numeris = { version = "0.5.6", features = ["nalgebra"] }
+```
+
+This provides zero-cost `From`/`Into` conversions between numeris and nalgebra matrix, vector, and dynamic matrix types. Both libraries use identical column-major storage, so conversions are a `memcpy`.
+
+### Python Bindings
+
+- No breaking changes to the Python API
+- All internal conversions updated for numeris types
+- Quaternion component order handling updated internally (transparent to Python users)
+
+### Code Simplification
+
+- Use `numeris::vector!` macro throughout instead of `Vector3::from_array([...])`
+- Use `Quaternion::rotation_between()` from numeris instead of hand-rolled helper
+- Remove `qrot_xcoord`/`qrot_ycoord`/`qrot_zcoord` wrappers; use `Quaternion::rotx`/`roty`/`rotz` directly
+- Remove `satkit::filters` module (UKF); use `numeris::estimate::Ukf` instead
+- Enable `estimate` feature on numeris
+
+### Dependency Cleanup
+
+- Remove `nalgebra` dependency entirely
+- Remove `ndarray` dependency (unused; `numpy` crate re-exports it for Python bindings)
+- Remove `cty` dependency; use `std::ffi::{c_double, c_int}` instead
+- Remove `once_cell` dependency; use `std::sync::OnceLock` (stable since Rust 1.70)
+- Remove redundant reference-based `Add`/`Sub` operator impls from `ITRFCoord` (both types are `Copy`)
+
+### Internal
+
+- Net reduction of ~11,000 lines of code (removed ODE module and filters module)
+- 141 Rust tests pass (105 lib + 36 doc)
+- 52 Python tests pass
+
+
+## 0.13.0 - 2026-03-15
+
+### Integrator and Gravity Model Selection
+
+- Add `Integrator` enum to `PropSettings` for selecting ODE solver (RKV98, RKV87, RKV65, RKTS54, RODAS4)
+- Add `GravityModel` enum for selecting Earth gravity model (JGM3, JGM2, EGM96, ITU GRACE16)
+- Optimize ODE solver hot paths
+
+
 ## 0.12.0 - 2026-03-02
 
 ### API Improvements

diff --git a/Cargo.toml b/Cargo.toml
@@ -3,7 +3,7 @@ members = [".", "python"]
 
 [package]
 name = "satkit"
-version = "0.13.0"
+version = "0.14.1"
 edition = "2021"
 description = "Satellite Toolkit"
 readme = "README.md"
@@ -22,12 +22,9 @@ name = "satkit"
 
 
 [dependencies]
-nalgebra = { version = "0.34.0", features = ["serde-serialize"] }
-ndarray = "0.17.1"
-cty = "0.2.2"
+numeris = { version = "0.5.6", features = ["serde", "estimate", "ode"] }
 num-traits = "0.2.19"
 thiserror = "2.0"
-once_cell = "1.21"
 ureq = "3.1.2"
 json = "0.12.4"
 process_path = "0.1.4"
@@ -41,7 +38,6 @@ rmpfit = "0.3.0"
 chrono = { version = "0.4", optional = true }
 
 [build-dependencies]
-cc = { version = "1.2", features = ["parallel"] }
 chrono = "0.4"
 
 [dev-dependencies]

diff --git a/README.md b/README.md
@@ -141,6 +141,14 @@ Seamless conversion between UTC, TAI, TT, TDB, UT1, and GPS time scales with ful
 - Fast analytical Sun/Moon models for lower-precision work
 - Sunrise/sunset and Moon phase calculations
 
+### Linear Algebra
+
+SatKit uses [numeris](https://crates.io/crates/numeris) for all linear algebra (vectors, matrices, quaternions, ODE integration). If you also use nalgebra in your project, enable the `nalgebra` feature on numeris for zero-cost `From`/`Into` conversions between types:
+
+```toml
+numeris = { version = "0.5.5", features = ["nalgebra"] }
+```
+
 ### Cargo Features
 
 | Feature | Default | Description |
@@ -165,7 +173,7 @@ The library is validated against:
 - **ICGEM** reference values for gravity field calculations
 - **GPS SP3** precise ephemerides for multi-day numerical propagation
 
-106 unit tests and 36 doc-tests run on every commit across Linux, macOS, and Windows.
+142 tests (106 unit + 36 doc-tests) run on every commit across Linux, macOS, and Windows.
 
 ## Documentation
 

diff --git a/build.rs b/build.rs
@@ -1,10 +1,7 @@
 use std::process::Command;
 
 fn main() {
-    cc::Build::new()
-        .file("extern/nrlmsise/nrlmsise-00.c")
-        .file("extern/nrlmsise/nrlmsise-00_data.c")
-        .compile("nrlmsise");
+    // NRLMSISE-00 is now pure Rust (src/nrlmsise.rs), no C compilation needed
 
     // Record git hash to compile-time environment variable
     let output = Command::new("git")

diff --git a/docs/api/lambert.md b/docs/api/lambert.md
@@ -0,0 +1,3 @@
+# Lambert Solver
+
+::: satkit.lambert
diff --git a/docs/examples/gen_plots.py b/docs/examples/gen_plots.py
@@ -0,0 +1,168 @@
+#!/usr/bin/env python3
+"""Generate SVG plots for MkDocs documentation using matplotlib + SciencePlots.
+
+    python docs/examples/gen_plots.py
+
+Writes SVG files to docs/images/.
+"""
+
+import os
+from pathlib import Path
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+
+import scienceplots  # noqa: F401
+
+plt.style.use(["science", "no-latex"])
+
+plt.rcParams.update({
+    "mathtext.fontset": "stix",
+    "font.family": "serif",
+    "font.serif": ["STIX Two Text", "STIXGeneral", "DejaVu Serif"],
+    "font.size": 13,
+    "axes.labelsize": 14,
+    "axes.titlesize": 15,
+    "xtick.labelsize": 12,
+    "ytick.labelsize": 12,
+    "legend.fontsize": 12,
+    "axes.formatter.use_mathtext": True,
+    "svg.fonttype": "none",
+    "axes.grid": True,
+    "grid.alpha": 0.3,
+    "grid.linestyle": "--",
+    "axes.prop_cycle": plt.cycler(color=[
+        "#0077BB", "#EE7733", "#009988", "#CC3311",
+        "#33BBEE", "#EE3377", "#BBBBBB",
+    ]),
+})
+
+IMAGES = Path(__file__).resolve().parent.parent / "images"
+
+COLORS = [
+    "#0077BB",  # blue
+    "#EE7733",  # orange
+    "#009988",  # teal
+    "#CC3311",  # red
+    "#33BBEE",  # cyan
+    "#EE3377",  # magenta
+    "#BBBBBB",  # grey
+]
+
+
+def savefig(fig, name):
+    path = IMAGES / f"{name}.svg"
+    fig.savefig(path, format="svg", bbox_inches="tight")
+    plt.close(fig)
+    print(f"  -> {path}")
+
+
+# ── Air Density vs Solar Cycle ─────────────────────────────────────────────
+
+def make_density_plot():
+    import satkit as sk
+
+    start = sk.time(1995, 1, 1)
+    end = sk.time(2022, 12, 31)
+    duration = end - start
+    timearray = [start + sk.duration(days=x)
+                 for x in np.linspace(0, duration.days, 4000)]
+
+    rho_400 = [sk.density.nrlmsise(400e3, 0, 0, x)[0] for x in timearray]
+    rho_500 = [sk.density.nrlmsise(500e3, 0, 0, x)[0] for x in timearray]
+
+    fig, ax = plt.subplots(figsize=(8, 4.5))
+    dates = [t.as_datetime() for t in timearray]
+    ax.semilogy(dates, rho_400, color=COLORS[0], linewidth=1,
+                label="Altitude = 400 km")
+    ax.semilogy(dates, rho_500, color=COLORS[1], linewidth=1,
+                linestyle="--", label="Altitude = 500 km")
+    ax.set_xlabel("Year")
+    ax.set_ylabel(r"Density [kg/m$^3$]")
+    ax.set_title("Air Density Changes with Solar Cycle")
+    ax.legend()
+    fig.autofmt_xdate()
+    savefig(fig, "density_vs_solar_cycle")
+
+
+# ── Satellite Forces vs Altitude ───────────────────────────────────────────
+
+def make_force_plot():
+    import satkit as sk
+    import math as m
+
+    N = 1000
+    range_arr = np.logspace(m.log10(6378.2e3 + 100e3), m.log10(50e6), N)
+
+    grav1v = np.array([sk.gravity(np.array([a, 0, 0]), order=1)
+                       for a in range_arr])
+    grav1 = np.linalg.norm(grav1v, axis=1)
+
+    grav2v = np.array([sk.gravity(np.array([a, 0, 0]), order=2)
+                       for a in range_arr])
+    grav2 = np.linalg.norm(grav2v - grav1v, axis=1)
+
+    grav6v = np.array([sk.gravity(np.array([a, 0, 0]), order=6)
+                       for a in range_arr])
+    grav5v = np.array([sk.gravity(np.array([a, 0, 0]), order=5)
+                       for a in range_arr])
+    grav6 = np.linalg.norm(grav6v - grav5v, axis=1)
+
+    aoverm = 0.01
+    Cd = 2.2
+    Cr = 1.0
+
+    didx = np.argwhere(range_arr - sk.consts.earth_radius < 800e3).flatten()
+    tm_max = sk.time(2001, 12, 1)
+    tm_min = sk.time(1996, 12, 1)
+    rho_max = np.array([sk.density.nrlmsise(a - sk.consts.earth_radius, 0, 0, tm_max)[0]
+                        for a in range_arr[didx]])
+    rho_min = np.array([sk.density.nrlmsise(a - sk.consts.earth_radius, 0, 0, tm_min)[0]
+                        for a in range_arr[didx]])
+    varr = np.sqrt(sk.consts.mu_earth / (range_arr + sk.consts.earth_radius))
+    drag_max = 0.5 * rho_max * varr[didx]**2 * Cd * aoverm
+    drag_min = 0.5 * rho_min * varr[didx]**2 * Cd * aoverm
+
+    moon_range = np.linalg.norm(
+        sk.jplephem.geocentric_pos(sk.solarsystem.Moon, sk.time(2023, 1, 1)))
+    moon = sk.consts.mu_moon * (
+        (moon_range - range_arr)**(-2) - moon_range**(-2))
+    sun = sk.consts.mu_sun * (
+        (sk.consts.au - range_arr)**(-2) - sk.consts.au**(-2))
+
+    a_radiation = 4.56e-6 * 0.5 * Cr * aoverm * np.ones(range_arr.shape)
+
+    def add_line(ax, x, y, text, frac=0.5, dx=-20, dy=-20):
+        ax.loglog(x, y, "k-", linewidth=1.5)
+        idx = int(len(x) * frac)
+        ax.annotate(text, xy=(x[idx], y[idx]), fontsize=10,
+                    xytext=(dx, dy), textcoords="offset points",
+                    arrowprops=dict(arrowstyle="->", color="gray"))
+
+    fig, ax = plt.subplots(figsize=(8, 8))
+    add_line(ax, range_arr / 1e3, grav1 / 1e3, "Gravity")
+    add_line(ax, range_arr / 1e3, grav2 / 1e3, "J2", 0.2, 0, -15)
+    add_line(ax, range_arr / 1e3, grav6 / 1e3, "J6", 0.8, 0, -15)
+    add_line(ax, range_arr[didx] / 1e3, drag_max / 1e3,
+             "Drag Max", 0.7, 30, 0)
+    add_line(ax, range_arr[didx] / 1e3, drag_min / 1e3,
+             "Drag Min", 0.8, 10, 30)
+    add_line(ax, range_arr / 1e3, moon / 1e3, "Moon", 0.8, -10, -15)
+    add_line(ax, range_arr / 1e3, sun / 1e3, "Sun", 0.7, -10, 15)
+    add_line(ax, range_arr / 1e3, a_radiation / 1e3,
+             "Radiation\nPressure", 0.3, -10, 15)
+    ax.set_xlabel("Distance from Earth Origin [km]")
+    ax.set_ylabel(r"Acceleration [km/s$^2$]")
+    ax.set_title("Satellite Forces vs Altitude")
+    ax.set_xlim(6378.1, 50e3)
+    savefig(fig, "force_vs_altitude")
+
+
+if __name__ == "__main__":
+    os.makedirs(IMAGES, exist_ok=True)
+    print("Generating SVG plots...")
+    make_density_plot()
+    make_force_plot()
+    print("Done.")