From 4b034bcab6dbd28864e90db13f797c52c1aefac9 Mon Sep 17 00:00:00 2001 From: Chase Bryan <281080815+Bryforge@users.noreply.github.com> Date: Fri, 22 May 2026 22:36:25 -0500 Subject: [PATCH 1/2] Add Rust core backend and update README instructions --- .github/workflows/ci.yml | 40 +++ .gitignore | 7 + Cargo.lock | 133 ++++++++++ Cargo.toml | 3 + Makefile | 25 ++ README.md | 322 ++++++++++++++++++++---- crates/mcdowell-arc-core/Cargo.toml | 18 ++ crates/mcdowell-arc-core/README.md | 13 + crates/mcdowell-arc-core/pyproject.toml | 13 + crates/mcdowell-arc-core/src/lib.rs | 304 ++++++++++++++++++++++ docs/ARCHITECTURE.md | 30 +++ docs/SCIENTIFIC_MODEL.md | 45 ++++ docs/VALIDATION.md | 34 +++ pyproject.toml | 21 +- scripts/benchmark_monte_carlo.py | 26 ++ src/mcdowell_arc/__init__.py | 2 +- src/mcdowell_arc/cli.py | 52 +++- src/mcdowell_arc/dynamics.py | 61 +++-- src/mcdowell_arc/fast.py | 88 +++++++ src/mcdowell_arc/fit.py | 134 ++++++++-- tests/test_atmosphere.py | 12 + tests/test_dynamics.py | 32 +++ tests/test_fit_cli.py | 38 +++ 23 files changed, 1366 insertions(+), 87 deletions(-) create mode 100644 .github/workflows/ci.yml create mode 100644 Cargo.lock create mode 100644 Cargo.toml create mode 100644 Makefile create mode 100644 crates/mcdowell-arc-core/Cargo.toml create mode 100644 crates/mcdowell-arc-core/README.md create mode 100644 crates/mcdowell-arc-core/pyproject.toml create mode 100644 crates/mcdowell-arc-core/src/lib.rs create mode 100644 docs/ARCHITECTURE.md create mode 100644 docs/SCIENTIFIC_MODEL.md create mode 100644 docs/VALIDATION.md create mode 100644 scripts/benchmark_monte_carlo.py create mode 100644 src/mcdowell_arc/fast.py create mode 100644 tests/test_atmosphere.py create mode 100644 tests/test_dynamics.py create mode 100644 tests/test_fit_cli.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..b5d2ee7 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,40 @@ +name: CI + +on: + push: + pull_request: + +jobs: + python: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.10", "3.11", "3.12", "3.13"] + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - name: Install package + run: | + python -m pip install --upgrade pip setuptools wheel + python -m pip install -e ".[dev]" + - name: Lint + run: python -m ruff check src tests + - name: Test + run: python -m pytest -q + + rust: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: "3.12" + - uses: actions-rust-lang/setup-rust-toolchain@v1 + - name: Rust tests + run: cargo test --workspace + - name: Build Rust Python core + run: | + python -m pip install --upgrade pip maturin + maturin build --manifest-path crates/mcdowell-arc-core/Cargo.toml --release diff --git a/.gitignore b/.gitignore index 4de96c8..cf827d9 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,10 @@ build/ *.egg-info/ .env .DS_Store +.ruff_cache/ +.mypy_cache/ +target/ +*.so +*.dylib +*.dll +profile.out diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..869f053 --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,133 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "heck" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" + +[[package]] +name = "libc" +version = "0.2.186" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "68ab91017fe16c622486840e4c83c9a37afeff978bd239b5293d61ece587de66" + +[[package]] +name = "mcdowell-arc-core" +version = "0.2.0" +dependencies = [ + "pyo3", +] + +[[package]] +name = "once_cell" +version = "1.21.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50" + +[[package]] +name = "portable-atomic" +version = "1.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c33a9471896f1c69cecef8d20cbe2f7accd12527ce60845ff44c153bb2a21b49" + +[[package]] +name = "proc-macro2" +version = "1.0.106" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8fd00f0bb2e90d81d1044c2b32617f68fcb9fa3bb7640c23e9c748e53fb30934" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pyo3" +version = "0.28.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91fd8e38a3b50ed1167fb981cd6fd60147e091784c427b8f7183a7ee32c31c12" +dependencies = [ + "libc", + "once_cell", + "portable-atomic", + "pyo3-build-config", + "pyo3-ffi", + "pyo3-macros", +] + +[[package]] +name = "pyo3-build-config" +version = "0.28.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e368e7ddfdeb98c9bca7f8383be1648fd84ab466bf2bc015e94008db6d35611e" +dependencies = [ + "target-lexicon", +] + +[[package]] +name = "pyo3-ffi" +version = "0.28.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7f29e10af80b1f7ccaf7f69eace800a03ecd13e883acfacc1e5d0988605f651e" +dependencies = [ + "libc", + "pyo3-build-config", +] + +[[package]] +name = "pyo3-macros" +version = "0.28.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "df6e520eff47c45997d2fc7dd8214b25dd1310918bbb2642156ef66a67f29813" +dependencies = [ + "proc-macro2", + "pyo3-macros-backend", + "quote", + "syn", +] + +[[package]] +name = "pyo3-macros-backend" +version = "0.28.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c4cdc218d835738f81c2338f822078af45b4afdf8b2e33cbb5916f108b813acb" +dependencies = [ + "heck", + "proc-macro2", + "pyo3-build-config", + "quote", + "syn", +] + +[[package]] +name = "quote" +version = "1.0.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "41f2619966050689382d2b44f664f4bc593e129785a36d6ee376ddf37259b924" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "syn" +version = "2.0.117" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e665b8803e7b1d2a727f4023456bbbbe74da67099c585258af0ad9c5013b9b99" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "target-lexicon" +version = "0.13.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adb6935a6f5c20170eeceb1a3835a49e12e19d792f6dd344ccc76a985ca5a6ca" + +[[package]] +name = "unicode-ident" +version = "1.0.24" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6e4313cd5fcd3dad5cafa179702e2b244f760991f45397d14d4ebf38247da75" diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..648f31e --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,3 @@ +[workspace] +members = ["crates/mcdowell-arc-core"] +resolver = "2" diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..2fcd30b --- /dev/null +++ b/Makefile @@ -0,0 +1,25 @@ +.PHONY: install dev test lint doctor rust-core clean + +install: + python -m pip install -e . + +dev: + python -m pip install --upgrade pip setuptools wheel + python -m pip install -e ".[dev]" + +test: + python -m pytest -q + +lint: + python -m ruff check src tests + +doctor: + mcdowell-arc doctor + +rust-core: + python -m pip install maturin + maturin develop --manifest-path crates/mcdowell-arc-core/Cargo.toml --release + +clean: + rm -rf .pytest_cache .ruff_cache build dist src/*.egg-info target profile.out + find . -type d -name __pycache__ -prune -exec rm -rf {} + diff --git a/README.md b/README.md index 9e19d2c..d2e2a96 100644 --- a/README.md +++ b/README.md @@ -1,80 +1,316 @@ # McDowell Arc -**McDowell Arc** is a small starter tool for drag-aware, Earth-centered inertial trajectory fitting from webcast-style telemetry. It is named as a working-title tribute to Jonathan McDowell's challenge: do better than a clean ballistic ellipse fitted near apogee. +**McDowell Arc** is an experimental scientific software tool for drag-aware trajectory fitting and uncertainty analysis from simplified webcast-style trajectory observations. -The MVP takes timestamped altitude and Earth-relative speed samples, fits a reduced 2-D ECI trajectory with atmospheric drag, and reports a nominal perigee/apogee plus Monte Carlo confidence bounds. The primary result is not just “the fitted perigee,” but the probability that perigee rises above a configurable atmosphere threshold. +It estimates nominal perigee/apogee values, confidence intervals, and the probability that a fitted perigee remains above a selected atmospheric threshold. -> Public naming note: before publishing under a living person's name, consider asking permission or using a neutral package name with an acknowledgement section. +The v0.2 architecture separates the project into two parts: -## What it answers +- **`mcdowell-arc`** — Python CLI, CSV loading, reporting, validation, and scientific workflow. +- **`mcdowell-arc-core`** — Rust acceleration core exposed to Python through a native extension. -Instead of asking only: +The goal is to keep the tool easy to use from Python while moving performance-sensitive trajectory calculations into a faster compiled backend. -> What ellipse fits these samples? +> Named in honor of Jonathan McDowell's challenge to improve near-orbital trajectory fitting. This project is not affiliated with or endorsed by Jonathan McDowell unless explicitly stated. -McDowell Arc asks: +--- -> Given measurement uncertainty and drag, what is the probability the fitted perigee remains above the atmosphere? +## Current Status -## Current MVP scope +McDowell Arc is currently an engineering MVP moving toward production-quality scientific software. -Implemented now: +It includes: -- Webcast CSV ingestion: `t_s`, `altitude_km`, `speed_km_s` -- Optional uncertainty columns: `sigma_altitude_km`, `sigma_speed_km_s` -- Earth rotation correction in the drag-relative velocity term -- Simplified co-rotating atmosphere -- Built-in rough density table for 0-1000 km altitude -- Weighted least-squares fitting +- Webcast-style CSV ingestion +- Weighted trajectory fitting +- Simplified drag-aware dynamics +- Earth-rotation-aware drag-relative velocity - Monte Carlo uncertainty propagation -- Perigee/apogee reporting +- Nominal perigee/apogee reporting +- Confidence interval reporting +- Probability estimate for perigee above an atmospheric threshold +- Python backend +- Rust acceleration backend through `mcdowell-arc-core` +- Backend selection with `--backend auto|python|rust` +- Environment diagnostics with `mcdowell-arc doctor` +- Test coverage for Python and Rust-backed behavior -Important limitations: +Scientific caveat: the current model is still simplified. It is useful for development, experimentation, validation workflows, and early trajectory-analysis research, but it is not yet publication-grade orbital determination software. -- The current fitter is a reduced 2-D approximation. Full ECI reconstruction needs position/velocity direction information, such as lat/lon/range/range-rate, tracking geometry, or assumptions about launch azimuth/inclination. -- The built-in atmosphere is intentionally simple. Replace it with NRLMSISE-00 or another vetted density model before serious public claims. -- Near-apogee-only telemetry can be underconstrained. The tool reports uncertainty so that weak fits remain visibly weak. +--- -## Install +## Repository Layout + +```text +mcdowell-arc/ +├── mcdowell_arc/ # Python package and CLI +├── crates/ +│ └── mcdowell-arc-core/ # Rust/PyO3 acceleration core +├── examples/ # Example CSV input data +├── tests/ # Python test suite +├── docs/ # Scientific and validation notes +├── pyproject.toml # Python package metadata +├── Cargo.toml # Rust workspace metadata +└── README.md +``` + +--- + +## Requirements + +### Python + +Python 3.10 or newer is recommended. + +The project has been tested on Fedora with Python 3.14. + +### Rust + +The Rust backend requires Rust and Cargo. + +On Fedora: + +```bash +sudo dnf install -y rust cargo +``` + +You can also install Rust through `rustup`. + +--- + +## Quick Start + +Clone the repository: + +```bash +git clone https://github.com/Bryforge/mcdowell-arc.git +cd mcdowell-arc +``` + +Create and activate a virtual environment: ```bash python3 -m venv .venv source .venv/bin/activate -pip install -e . ``` -## Run the example +Upgrade packaging tools: ```bash -mcdowell-arc fit examples/sample_webcast.csv --monte-carlo 20 --atmosphere-km 100 +python -m pip install --upgrade pip setuptools wheel ``` -## CSV format +Install the Python package in editable development mode: -Required columns: +```bash +python -m pip install -e ".[dev]" +``` + +Run the tests: + +```bash +python -m pytest -q +``` + +Run diagnostics: + +```bash +mcdowell-arc doctor +``` + +At this point, the Python backend should work. + +--- + +## Building the Rust Acceleration Core + +Install `maturin`: + +```bash +python -m pip install maturin +``` + +Build and install the Rust extension from inside the Rust crate directory: + +```bash +cd crates/mcdowell-arc-core +maturin develop --release +cd ../.. +``` + +Verify that both packages are installed: + +```bash +python -m pip list | grep mcdowell +``` + +Expected result: + +```text +mcdowell-arc 0.2.0 +mcdowell-arc-core 0.2.0 +``` + +Verify that Python can import both modules: + +```bash +python -c "import mcdowell_arc; print('Python package loaded')" +python -c "import mcdowell_arc_core; print('Rust core loaded')" +``` + +Run the full test suite: + +```bash +python -m pytest -q +``` + +Expected result: + +```text +12 passed +``` + +--- + +## Usage + +Run a fit using automatic backend selection: + +```bash +mcdowell-arc fit examples/sample_webcast.csv --monte-carlo 1000 --atmosphere-km 100 --backend auto +``` + +Force the Python backend: + +```bash +mcdowell-arc fit examples/sample_webcast.csv --monte-carlo 1000 --atmosphere-km 100 --backend python +``` + +Force the Rust backend: + +```bash +mcdowell-arc fit examples/sample_webcast.csv --monte-carlo 1000 --atmosphere-km 100 --backend rust +``` + +If the Rust backend is not installed and `--backend rust` is selected, the command will fail. If `--backend auto` is selected, the tool will use Rust when available and fall back to Python when needed. + +--- + +## Input CSV Format + +The basic input CSV should contain: ```csv t_s,altitude_km,speed_km_s -0,210,7.81 -10,211,7.80 -20,210.8,7.79 +0,160.0,7.72 +10,158.4,7.70 +20,156.8,7.68 ``` -Recommended uncertainty columns: +Required columns: -```csv -t_s,altitude_km,speed_km_s,sigma_altitude_km,sigma_speed_km_s -0,210,7.81,3,0.05 -10,211,7.80,3,0.05 +- `t_s` — time in seconds +- `altitude_km` — observed altitude in kilometers +- `speed_km_s` — observed speed in kilometers per second + +Optional uncertainty columns may be added in future model expansions. + +--- + +## Example Output + +The CLI emits JSON output containing fields such as: + +```json +{ + "optimizer_success": true, + "backend": "rust", + "samples_requested": 1000, + "samples_succeeded": 1000, + "samples_failed": 0, + "nominal": { + "perigee_km": 129.0, + "apogee_km": 2242.0 + }, + "probability_perigee_above_threshold": 1.0 +} +``` + +Exact values depend on the input data, random seed behavior, model settings, and backend. + +--- + +## Development Workflow + +Install development dependencies: + +```bash +python -m pip install -e ".[dev]" +python -m pip install maturin ``` -## Roadmap +Build the Rust backend: + +```bash +cd crates/mcdowell-arc-core +maturin develop --release +cd ../.. +``` + +Run tests: + +```bash +python -m pytest -q +``` + +Run diagnostics: + +```bash +mcdowell-arc doctor +``` + +Run the example fit: + +```bash +mcdowell-arc fit examples/sample_webcast.csv --monte-carlo 1000 --atmosphere-km 100 --backend auto +``` + +--- + +## Scientific Roadmap + +Major next steps: + +1. Full 3-D trajectory state modeling. +2. Observation geometry support: latitude, longitude, azimuth, range, range-rate, and launch-site assumptions. +3. Improved atmosphere modeling. +4. Better uncertainty modeling for webcast-derived measurements. +5. Validation against synthetic truth cases. +6. Validation against known historical trajectories where public data is available. +7. Parallel Monte Carlo execution. +8. Benchmarking Python backend versus Rust backend. +9. Plotting and report generation. +10. Formal model documentation. + +--- + +## Production Readiness Goals + +McDowell Arc should evolve toward: + +- Reproducible scientific runs +- Versioned model assumptions +- Clear input schemas +- Deterministic test cases +- Synthetic validation datasets +- Backend parity tests +- CI testing for Python and Rust +- Structured JSON output +- Scientific documentation +- Honest uncertainty reporting + +--- -1. Add full ECEF/geodetic to ECI observation support. -2. Add real Earth orientation / sidereal-time handling. -3. Add NRLMSISE-00 density as a plug-in backend. -4. Support boost/coast/drag phase boundaries. -5. Add covariance estimation from the least-squares Jacobian. -6. Add plots for altitude, speed, residuals, and perigee distribution. -7. Add a reproducible report generator suitable for public trajectory notes. +## License +Add a license before wider public release. diff --git a/crates/mcdowell-arc-core/Cargo.toml b/crates/mcdowell-arc-core/Cargo.toml new file mode 100644 index 0000000..f74227e --- /dev/null +++ b/crates/mcdowell-arc-core/Cargo.toml @@ -0,0 +1,18 @@ +[package] +name = "mcdowell-arc-core" +version = "0.2.0" +edition = "2021" +description = "Fast Rust core for McDowell Arc trajectory propagation and orbit summaries." +license = "MIT" + +[lib] +name = "mcdowell_arc_core" +crate-type = ["cdylib", "rlib"] + +[dependencies] +pyo3 = { version = "0.28", features = ["extension-module"] } + +[profile.release] +codegen-units = 1 +lto = "thin" +opt-level = 3 diff --git a/crates/mcdowell-arc-core/README.md b/crates/mcdowell-arc-core/README.md new file mode 100644 index 0000000..bcd0495 --- /dev/null +++ b/crates/mcdowell-arc-core/README.md @@ -0,0 +1,13 @@ +# mcdowell-arc-core + +`mcdowell-arc-core` is the optional Rust acceleration layer for McDowell Arc. +It exposes a Python extension module named `mcdowell_arc_core` through PyO3. + +The Python package works without this crate. Install the Rust core when you want +faster deterministic propagation inside the existing `mcdowell-arc` CLI. + +```bash +python -m pip install maturin +maturin develop --manifest-path crates/mcdowell-arc-core/Cargo.toml --release +mcdowell-arc doctor +``` diff --git a/crates/mcdowell-arc-core/pyproject.toml b/crates/mcdowell-arc-core/pyproject.toml new file mode 100644 index 0000000..12c7073 --- /dev/null +++ b/crates/mcdowell-arc-core/pyproject.toml @@ -0,0 +1,13 @@ +[build-system] +requires = ["maturin>=1.8,<2"] +build-backend = "maturin" + +[project] +name = "mcdowell-arc-core" +version = "0.2.0" +description = "Rust acceleration core for McDowell Arc" +requires-python = ">=3.10" + +[tool.maturin] +features = ["pyo3/extension-module"] +module-name = "mcdowell_arc_core" diff --git a/crates/mcdowell-arc-core/src/lib.rs b/crates/mcdowell-arc-core/src/lib.rs new file mode 100644 index 0000000..ce24385 --- /dev/null +++ b/crates/mcdowell-arc-core/src/lib.rs @@ -0,0 +1,304 @@ +use pyo3::exceptions::{PyRuntimeError, PyValueError}; +use pyo3::prelude::*; + +const MU_EARTH_KM3_S2: f64 = 398_600.4418; +const R_EARTH_KM: f64 = 6_378.137; +const OMEGA_EARTH_RAD_S: f64 = 7.292_115_9e-5; + +const DENSITY_TABLE: &[(f64, f64)] = &[ + (0.0, 1.225), + (25.0, 3.899e-2), + (30.0, 1.774e-2), + (40.0, 3.972e-3), + (50.0, 1.057e-3), + (60.0, 3.206e-4), + (70.0, 8.770e-5), + (80.0, 1.905e-5), + (90.0, 3.396e-6), + (100.0, 5.297e-7), + (110.0, 9.661e-8), + (120.0, 2.438e-8), + (130.0, 8.484e-9), + (140.0, 3.845e-9), + (150.0, 2.070e-9), + (180.0, 5.464e-10), + (200.0, 2.789e-10), + (250.0, 7.248e-11), + (300.0, 2.418e-11), + (350.0, 9.518e-12), + (400.0, 3.725e-12), + (450.0, 1.585e-12), + (500.0, 6.967e-13), + (600.0, 1.454e-13), + (700.0, 3.614e-14), + (800.0, 1.170e-14), + (900.0, 5.245e-15), + (1000.0, 3.019e-15), +]; + +#[inline] +fn check_finite(value: f64, name: &str) -> PyResult<()> { + if value.is_finite() { + Ok(()) + } else { + Err(PyValueError::new_err(format!("{name} must be finite"))) + } +} + +#[inline] +fn norm3(v: &[f64; 3]) -> f64 { + (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt() +} + +#[inline] +fn cross(a: &[f64; 3], b: &[f64; 3]) -> [f64; 3] { + [ + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0], + ] +} + +#[inline] +fn vec_add(a: &[f64; 6], b: &[f64; 6], scale: f64) -> [f64; 6] { + [ + a[0] + scale * b[0], + a[1] + scale * b[1], + a[2] + scale * b[2], + a[3] + scale * b[3], + a[4] + scale * b[4], + a[5] + scale * b[5], + ] +} + +fn density_scalar(altitude_km: f64) -> f64 { + let mut alt = altitude_km; + let min_alt = DENSITY_TABLE[0].0; + let max_alt = DENSITY_TABLE[DENSITY_TABLE.len() - 1].0; + if alt < min_alt { + alt = min_alt; + } + if alt > max_alt { + alt = max_alt; + } + + for window in DENSITY_TABLE.windows(2) { + let (h0, rho0) = window[0]; + let (h1, rho1) = window[1]; + if alt >= h0 && alt <= h1 { + let frac = if h1 > h0 { (alt - h0) / (h1 - h0) } else { 0.0 }; + let log_rho = rho0.ln() + frac * (rho1.ln() - rho0.ln()); + return log_rho.exp(); + } + } + DENSITY_TABLE[DENSITY_TABLE.len() - 1].1 +} + +fn rotation_velocity(r: &[f64; 3]) -> [f64; 3] { + [-OMEGA_EARTH_RAD_S * r[1], OMEGA_EARTH_RAD_S * r[0], 0.0] +} + +fn acceleration(state: &[f64; 6], ballistic_coef_kg_m2: f64, use_drag: bool) -> PyResult<[f64; 3]> { + let r = [state[0], state[1], state[2]]; + let v = [state[3], state[4], state[5]]; + let r_norm = norm3(&r); + if r_norm <= 0.0 || !r_norm.is_finite() { + return Err(PyValueError::new_err("position norm must be positive and finite")); + } + + let mut a = [ + -MU_EARTH_KM3_S2 * r[0] / r_norm.powi(3), + -MU_EARTH_KM3_S2 * r[1] / r_norm.powi(3), + -MU_EARTH_KM3_S2 * r[2] / r_norm.powi(3), + ]; + + if !use_drag { + return Ok(a); + } + if ballistic_coef_kg_m2 <= 0.0 || !ballistic_coef_kg_m2.is_finite() { + return Ok(a); + } + + let alt_km = (r_norm - R_EARTH_KM).max(0.0); + let rho = density_scalar(alt_km); + let v_atm = rotation_velocity(&r); + let v_rel = [v[0] - v_atm[0], v[1] - v_atm[1], v[2] - v_atm[2]]; + let v_rel_norm_km_s = norm3(&v_rel); + let v_rel_m_s = v_rel_norm_km_s * 1000.0; + + if v_rel_m_s <= 0.0 || !v_rel_m_s.is_finite() { + return Ok(a); + } + + let drag_mag_m_s2 = -0.5 * rho * v_rel_m_s * v_rel_m_s / ballistic_coef_kg_m2; + let unit = [v_rel[0] / v_rel_norm_km_s, v_rel[1] / v_rel_norm_km_s, v_rel[2] / v_rel_norm_km_s]; + a[0] += drag_mag_m_s2 * unit[0] / 1000.0; + a[1] += drag_mag_m_s2 * unit[1] / 1000.0; + a[2] += drag_mag_m_s2 * unit[2] / 1000.0; + Ok(a) +} + +fn rhs(state: &[f64; 6], ballistic_coef_kg_m2: f64, use_drag: bool) -> PyResult<[f64; 6]> { + let a = acceleration(state, ballistic_coef_kg_m2, use_drag)?; + Ok([state[3], state[4], state[5], a[0], a[1], a[2]]) +} + +fn rk4_step(state: &[f64; 6], dt_s: f64, ballistic_coef_kg_m2: f64, use_drag: bool) -> PyResult<[f64; 6]> { + let k1 = rhs(state, ballistic_coef_kg_m2, use_drag)?; + let s2 = vec_add(state, &k1, 0.5 * dt_s); + let k2 = rhs(&s2, ballistic_coef_kg_m2, use_drag)?; + let s3 = vec_add(state, &k2, 0.5 * dt_s); + let k3 = rhs(&s3, ballistic_coef_kg_m2, use_drag)?; + let s4 = vec_add(state, &k3, dt_s); + let k4 = rhs(&s4, ballistic_coef_kg_m2, use_drag)?; + + let mut out = [0.0; 6]; + for i in 0..6 { + out[i] = state[i] + (dt_s / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]); + } + Ok(out) +} + +fn state6_from_vec(values: Vec, name: &str) -> PyResult<[f64; 6]> { + if values.len() != 6 { + return Err(PyValueError::new_err(format!("{name} must contain exactly 6 values"))); + } + for (idx, value) in values.iter().enumerate() { + check_finite(*value, &format!("{name}[{idx}]"))?; + } + Ok([values[0], values[1], values[2], values[3], values[4], values[5]]) +} + +#[pyfunction] +fn density_kg_m3(altitude_km: f64) -> PyResult { + check_finite(altitude_km, "altitude_km")?; + Ok(density_scalar(altitude_km)) +} + +#[pyfunction] +fn earth_rotation_velocity_km_s(r_eci_km: Vec) -> PyResult> { + if r_eci_km.len() != 3 { + return Err(PyValueError::new_err("r_eci_km must contain exactly 3 values")); + } + for (idx, value) in r_eci_km.iter().enumerate() { + check_finite(*value, &format!("r_eci_km[{idx}]"))?; + } + let r = [r_eci_km[0], r_eci_km[1], r_eci_km[2]]; + Ok(rotation_velocity(&r).to_vec()) +} + +#[pyfunction] +fn propagate( + sample_times_s: Vec, + state0: Vec, + ballistic_coef_kg_m2: f64, + use_drag: bool, + max_step_s: f64, +) -> PyResult>> { + if sample_times_s.is_empty() { + return Err(PyValueError::new_err("sample_times_s must be non-empty")); + } + if max_step_s <= 0.0 || !max_step_s.is_finite() { + return Err(PyValueError::new_err("max_step_s must be positive and finite")); + } + check_finite(ballistic_coef_kg_m2, "ballistic_coef_kg_m2")?; + for (idx, t) in sample_times_s.iter().enumerate() { + check_finite(*t, &format!("sample_times_s[{idx}]"))?; + if idx > 0 && *t < sample_times_s[idx - 1] { + return Err(PyValueError::new_err("sample_times_s must be sorted ascending")); + } + } + + let mut state = state6_from_vec(state0, "state0")?; + let mut output = Vec::with_capacity(sample_times_s.len()); + output.push(state.to_vec()); + let mut t_current = sample_times_s[0]; + + for t_target in sample_times_s.iter().skip(1) { + while t_current < *t_target - 1e-12 { + let dt = max_step_s.min(*t_target - t_current); + state = rk4_step(&state, dt, ballistic_coef_kg_m2, use_drag)?; + t_current += dt; + if state.iter().any(|v| !v.is_finite()) { + return Err(PyRuntimeError::new_err("propagation produced a non-finite state")); + } + let r_norm = norm3(&[state[0], state[1], state[2]]); + if r_norm < R_EARTH_KM + 50.0 { + return Err(PyRuntimeError::new_err("propagation hit the reentry guard before all sample times")); + } + } + output.push(state.to_vec()); + } + + Ok(output) +} + +#[pyfunction] +fn summarize_orbit(state: Vec) -> PyResult<(f64, Option, f64, Option, f64)> { + let state = state6_from_vec(state, "state")?; + let r = [state[0], state[1], state[2]]; + let v = [state[3], state[4], state[5]]; + let r_norm = norm3(&r); + let v_norm = norm3(&v); + let h_vec = cross(&r, &v); + let h = norm3(&h_vec); + if r_norm <= 0.0 || h <= 0.0 { + return Err(PyValueError::new_err("degenerate orbit state")); + } + + let energy = 0.5 * v_norm * v_norm - MU_EARTH_KM3_S2 / r_norm; + let vxh = cross(&v, &h_vec); + let e_vec = [ + vxh[0] / MU_EARTH_KM3_S2 - r[0] / r_norm, + vxh[1] / MU_EARTH_KM3_S2 - r[1] / r_norm, + vxh[2] / MU_EARTH_KM3_S2 - r[2] / r_norm, + ]; + let ecc = norm3(&e_vec); + let p = h * h / MU_EARTH_KM3_S2; + let rp = p / (1.0 + ecc); + let perigee_km = rp - R_EARTH_KM; + + let (apogee_km, semi_major_axis_km) = if energy < 0.0 { + let a = -MU_EARTH_KM3_S2 / (2.0 * energy); + let ra = a * (1.0 + ecc); + (Some(ra - R_EARTH_KM), Some(a)) + } else { + (None, None) + }; + + Ok((perigee_km, apogee_km, ecc, semi_major_axis_km, energy)) +} + +#[pymodule] +fn mcdowell_arc_core(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add("__version__", env!("CARGO_PKG_VERSION"))?; + m.add("MU_EARTH_KM3_S2", MU_EARTH_KM3_S2)?; + m.add("R_EARTH_KM", R_EARTH_KM)?; + m.add("OMEGA_EARTH_RAD_S", OMEGA_EARTH_RAD_S)?; + m.add_function(wrap_pyfunction!(density_kg_m3, m)?)?; + m.add_function(wrap_pyfunction!(earth_rotation_velocity_km_s, m)?)?; + m.add_function(wrap_pyfunction!(propagate, m)?)?; + m.add_function(wrap_pyfunction!(summarize_orbit, m)?)?; + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn density_decreases_with_altitude() { + assert!(density_scalar(100.0) > density_scalar(150.0)); + assert!(density_scalar(150.0) > density_scalar(250.0)); + } + + #[test] + fn one_step_propagation_returns_two_rows() { + let r = R_EARTH_KM + 400.0; + let v = (MU_EARTH_KM3_S2 / r).sqrt(); + let out = propagate(vec![0.0, 10.0], vec![r, 0.0, 0.0, 0.0, v, 0.0], 250.0, false, 1.0) + .expect("propagation should succeed"); + assert_eq!(out.len(), 2); + assert_eq!(out[0].len(), 6); + } +} diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md new file mode 100644 index 0000000..577b8c4 --- /dev/null +++ b/docs/ARCHITECTURE.md @@ -0,0 +1,30 @@ +# Architecture + +McDowell Arc is split into two layers. + +## Python layer: `mcdowell_arc` + +Responsibilities: + +- command-line interface +- CSV ingestion and validation +- weighted least-squares fitting +- Monte Carlo uncertainty orchestration +- JSON reporting +- plotting and report generation in future releases + +## Rust layer: `mcdowell-arc-core` + +Responsibilities: + +- deterministic propagation kernels +- drag acceleration math +- orbit summary math +- future parallel Monte Carlo kernels +- future optimized observation transforms + +The Python package remains usable without Rust. The CLI option `--backend auto` selects Rust when the extension is installed and Python otherwise. + +## Why this split exists + +Python is better for fast scientific iteration and reproducible user workflows. Rust is better for deterministic, tested numerical kernels that may need to run thousands or millions of times. The design goal is not to replace Python, but to move hot loops into a small, auditable core. diff --git a/docs/SCIENTIFIC_MODEL.md b/docs/SCIENTIFIC_MODEL.md new file mode 100644 index 0000000..de35739 --- /dev/null +++ b/docs/SCIENTIFIC_MODEL.md @@ -0,0 +1,45 @@ +# Scientific Model Notes + +## Current model + +The current release uses a reduced two-dimensional ECI state: + +```text +state0 = [x, y, z, vx, vy, vz] + = [R_earth + altitude0 + offset, 0, 0, vr, vt, 0] +``` + +The observations are: + +- elapsed time in seconds +- altitude in kilometers +- Earth-relative speed magnitude in kilometers per second + +The model integrates gravity plus optional drag in a co-rotating atmosphere. + +## Why this is still not publication-grade + +Altitude and speed magnitude alone do not uniquely determine a full 3-D inertial state. A serious trajectory fit needs additional geometry such as: + +- latitude and longitude +- range and range rate +- azimuth/elevation from a known observer +- launch site and launch azimuth assumptions +- inclination or plane constraints +- external tracking observations + +The current model is useful for software plumbing, sensitivity testing, public telemetry experiments, and uncertainty framing. It should not be represented as a definitive orbital reconstruction. + +## Atmosphere caveat + +The built-in atmosphere table is a rough log-interpolated density table. It should be replaced by a vetted model before high-confidence scientific claims. + +## Correct public framing + +Strong framing: + +> This is an uncertainty-aware experimental fit under a reduced observation model. + +Weak or misleading framing: + +> This proves the object was orbital/non-orbital. diff --git a/docs/VALIDATION.md b/docs/VALIDATION.md new file mode 100644 index 0000000..7edf327 --- /dev/null +++ b/docs/VALIDATION.md @@ -0,0 +1,34 @@ +# Validation Plan + +This project should earn scientific credibility through layered validation. + +## Level 0: software sanity + +- unit tests for frame helpers +- unit tests for orbit summaries +- Python/Rust propagation agreement tests +- CLI smoke tests +- CSV validation tests + +## Level 1: synthetic truth + +Generate synthetic states with known perigee/apogee, then sample altitude/speed with noise. The fitter should recover the known state within uncertainty bounds under the same model. + +## Level 2: model sensitivity + +Run sensitivity sweeps over: + +- ballistic coefficient +- drag on/off +- time-window selection +- altitude uncertainty +- speed uncertainty +- atmosphere threshold + +## Level 3: independent event checks + +Compare against public cases where an independent orbit estimate exists. Report assumptions and residuals, not just final perigee. + +## Level 4: full geometry + +Move beyond webcast altitude/speed by supporting real observation geometry and Earth orientation handling. diff --git a/pyproject.toml b/pyproject.toml index ace966c..7cbb6dd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,8 +4,8 @@ build-backend = "setuptools.build_meta" [project] name = "mcdowell-arc" -version = "0.1.0" -description = "Drag-aware ECI webcast trajectory fitting with uncertainty bounds." +version = "0.2.0" +description = "Drag-aware ECI webcast trajectory fitting with uncertainty bounds and optional Rust core." readme = "README.md" requires-python = ">=3.10" authors = [{ name = "Bryforge" }] @@ -16,6 +16,14 @@ dependencies = [ "pandas>=2.0" ] +[project.optional-dependencies] +dev = [ + "pytest>=8.0", + "ruff>=0.6", + "build>=1.2", + "maturin>=1.8" +] + [project.scripts] mcdowell-arc = "mcdowell_arc.cli:main" @@ -25,3 +33,12 @@ where = ["src"] [tool.pytest.ini_options] pythonpath = ["src"] testpaths = ["tests"] +addopts = "-ra" + +[tool.ruff] +line-length = 110 +target-version = "py310" + +[tool.ruff.lint] +select = ["E", "F", "I", "UP", "B"] +ignore = ["E501"] diff --git a/scripts/benchmark_monte_carlo.py b/scripts/benchmark_monte_carlo.py new file mode 100644 index 0000000..640de7e --- /dev/null +++ b/scripts/benchmark_monte_carlo.py @@ -0,0 +1,26 @@ +from __future__ import annotations + +import argparse +import time + +from mcdowell_arc.fit import load_observations_csv, monte_carlo + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("csv") + parser.add_argument("--samples", type=int, default=100) + parser.add_argument("--backend", choices=["auto", "python", "rust"], default="auto") + args = parser.parse_args() + + obs = load_observations_csv(args.csv) + start = time.perf_counter() + result = monte_carlo(obs, args.samples, backend=args.backend) + elapsed = time.perf_counter() - start + print(f"backend={args.backend} samples={args.samples} elapsed_s={elapsed:.3f}") + print(result) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/mcdowell_arc/__init__.py b/src/mcdowell_arc/__init__.py index 9054c68..cf7cb5c 100644 --- a/src/mcdowell_arc/__init__.py +++ b/src/mcdowell_arc/__init__.py @@ -1,3 +1,3 @@ """McDowell Arc: drag-aware ECI webcast trajectory fitting.""" -__version__ = "0.1.0" +__version__ = "0.2.0" diff --git a/src/mcdowell_arc/cli.py b/src/mcdowell_arc/cli.py index 9386513..ef12508 100644 --- a/src/mcdowell_arc/cli.py +++ b/src/mcdowell_arc/cli.py @@ -4,13 +4,21 @@ import argparse import json -from dataclasses import asdict +import platform +import sys +import numpy as np +import pandas as pd +import scipy + +from . import __version__ +from .fast import core_available, core_version, resolve_backend from .fit import fit_observations, load_observations_csv, monte_carlo def build_parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser(prog="mcdowell-arc") + parser.add_argument("--version", action="version", version=f"mcdowell-arc {__version__}") sub = parser.add_subparsers(dest="command", required=True) fit = sub.add_parser("fit", help="fit a drag-aware ECI trajectory from webcast CSV") @@ -19,16 +27,48 @@ def build_parser() -> argparse.ArgumentParser: fit.add_argument("--atmosphere-km", type=float, default=100.0, help="perigee threshold for atmosphere crossing") fit.add_argument("--no-drag", action="store_true", help="disable drag for comparison") fit.add_argument("--seed", type=int, default=7, help="Monte Carlo random seed") + fit.add_argument("--backend", choices=["auto", "python", "rust"], default="auto", help="propagation backend") + fit.add_argument("--max-step-s", type=float, default=1.0, help="maximum fixed RK4 propagation step in seconds") + fit.add_argument("--max-nfev", type=int, default=35, help="maximum optimizer function evaluations") + + sub.add_parser("doctor", help="show environment and backend availability") return parser +def _doctor_payload() -> dict: + return { + "mcdowell_arc_version": __version__, + "python": sys.version.split()[0], + "platform": platform.platform(), + "numpy": np.__version__, + "scipy": scipy.__version__, + "pandas": pd.__version__, + "rust_core_available": core_available(), + "rust_core_version": core_version(), + "default_backend": resolve_backend("auto"), + } + + def main(argv: list[str] | None = None) -> int: args = build_parser().parse_args(argv) + if args.command == "doctor": + print(json.dumps(_doctor_payload(), indent=2, sort_keys=True)) + return 0 + if args.command == "fit": obs = load_observations_csv(args.csv) - nominal = fit_observations(obs, use_drag=not args.no_drag) + backend = resolve_backend(args.backend) + nominal = fit_observations( + obs, + use_drag=not args.no_drag, + backend=args.backend, + max_step_s=args.max_step_s, + max_nfev=args.max_nfev, + ) output = { - "model": "reduced_2d_eci_drag_mvp" if not args.no_drag else "reduced_2d_eci_no_drag_mvp", + "model": "reduced_2d_eci_drag" if not args.no_drag else "reduced_2d_eci_no_drag", + "model_status": "engineering_mvp_not_publication_grade", + "backend": backend, "nominal": { "perigee_km": nominal.orbit.perigee_km, "apogee_km": nominal.orbit.apogee_km, @@ -36,13 +76,15 @@ def main(argv: list[str] | None = None) -> int: "semi_major_axis_km": nominal.orbit.semi_major_axis_km, "ballistic_coef_kg_m2": nominal.ballistic_coef_kg_m2, "residual_rms_sigma": nominal.residual_rms, + "optimizer_cost": nominal.optimizer_cost, "optimizer_success": nominal.optimizer_success, "optimizer_message": nominal.message, }, "state0_eci_km_km_s": nominal.state0.tolist(), "interpretation": ( "Report probability_perigee_above_threshold, not only the nominal perigee. " - "MVP is reduced 2-D and needs full tracking geometry for publication-grade results." + "This model remains reduced 2-D and needs full tracking geometry, atmosphere validation, " + "and independent checks before publication-grade trajectory claims." ), } if args.monte_carlo > 0: @@ -52,6 +94,8 @@ def main(argv: list[str] | None = None) -> int: atmosphere_threshold_km=args.atmosphere_km, seed=args.seed, use_drag=not args.no_drag, + backend=args.backend, + max_step_s=args.max_step_s, ) print(json.dumps(output, indent=2, sort_keys=True)) return 0 diff --git a/src/mcdowell_arc/dynamics.py b/src/mcdowell_arc/dynamics.py index 9cb0e08..39368c9 100644 --- a/src/mcdowell_arc/dynamics.py +++ b/src/mcdowell_arc/dynamics.py @@ -3,8 +3,10 @@ from __future__ import annotations import numpy as np + from .atmosphere import density_kg_m3 from .constants import MU_EARTH_KM3_S2, R_EARTH_KM +from .fast import propagate_core, resolve_backend, validate_backend from .frames import earth_rotation_velocity_km_s @@ -16,14 +18,17 @@ def acceleration_km_s2( ) -> np.ndarray: """Compute gravity plus drag acceleration for a 6-state vector. - ballistic_coef_kg_m2 is beta = mass / (Cd * area). Larger beta means less drag. + `ballistic_coef_kg_m2` is beta = mass / (Cd * area). Larger beta means less drag. """ del t_s - r = np.asarray(state[:3], dtype=float) - v = np.asarray(state[3:], dtype=float) + state = np.asarray(state, dtype=float) + if state.shape != (6,): + raise ValueError("state must be a 6-element vector [x,y,z,vx,vy,vz].") + r = state[:3] + v = state[3:] r_norm = np.linalg.norm(r) - if r_norm <= 0: - raise ValueError("Position norm must be positive.") + if r_norm <= 0 or not np.isfinite(r_norm): + raise ValueError("Position norm must be positive and finite.") a_gravity = -MU_EARTH_KM3_S2 * r / r_norm**3 if not use_drag: @@ -33,13 +38,16 @@ def acceleration_km_s2( rho = density_kg_m3(alt_km) v_atm = earth_rotation_velocity_km_s(r) v_rel_km_s = v - v_atm - v_rel_m_s = np.linalg.norm(v_rel_km_s) * 1000.0 + v_rel_norm_km_s = np.linalg.norm(v_rel_km_s) + v_rel_m_s = v_rel_norm_km_s * 1000.0 if v_rel_m_s <= 0.0 or ballistic_coef_kg_m2 <= 0.0: return a_gravity # Drag acceleration magnitude: 0.5 * rho * v^2 / beta in m/s^2. - a_drag_m_s2_vec = -0.5 * rho * v_rel_m_s**2 / ballistic_coef_kg_m2 * (v_rel_km_s / np.linalg.norm(v_rel_km_s)) + a_drag_m_s2_vec = -0.5 * rho * v_rel_m_s**2 / ballistic_coef_kg_m2 * ( + v_rel_km_s / v_rel_norm_km_s + ) a_drag_km_s2 = a_drag_m_s2_vec / 1000.0 return a_gravity + a_drag_km_s2 @@ -58,27 +66,30 @@ def _rk4_step(t_s: float, state: np.ndarray, dt_s: float, ballistic_coef_kg_m2: return state + (dt_s / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4) -def propagate( +def _propagate_python( sample_times_s: np.ndarray, state0: np.ndarray, ballistic_coef_kg_m2: float, use_drag: bool = True, max_step_s: float = 1.0, ) -> np.ndarray: - """Propagate state0 to the requested sample times using fixed-step RK4. - - The fixed-step integrator is deliberately predictable for fitting. Impossible - trial states that dive below 50 km are rejected quickly. - """ + """Pure-Python/NumPy propagation implementation.""" times = np.asarray(sample_times_s, dtype=float) if times.ndim != 1 or len(times) == 0: raise ValueError("sample_times_s must be a non-empty 1-D array.") + if not np.all(np.isfinite(times)): + raise ValueError("sample_times_s must contain finite values.") if np.any(np.diff(times) < 0): raise ValueError("sample_times_s must be sorted ascending.") - if max_step_s <= 0: - raise ValueError("max_step_s must be positive.") + if max_step_s <= 0 or not np.isfinite(max_step_s): + raise ValueError("max_step_s must be positive and finite.") state = np.asarray(state0, dtype=float).copy() + if state.shape != (6,): + raise ValueError("state0 must be a 6-element vector [x,y,z,vx,vy,vz].") + if not np.all(np.isfinite(state)): + raise ValueError("state0 must contain only finite values.") + out = [state.copy()] t_current = float(times[0]) @@ -95,3 +106,23 @@ def propagate( out.append(state.copy()) return np.vstack(out) + + +def propagate( + sample_times_s: np.ndarray, + state0: np.ndarray, + ballistic_coef_kg_m2: float, + use_drag: bool = True, + max_step_s: float = 1.0, + backend: str = "auto", +) -> np.ndarray: + """Propagate state0 to the requested sample times using fixed-step RK4. + + `backend="auto"` uses the Rust core when installed and falls back to the + Python implementation otherwise. `backend="rust"` fails clearly when the + extension is unavailable. + """ + selected = validate_backend(backend) + if resolve_backend(selected) == "rust": + return propagate_core(sample_times_s, state0, ballistic_coef_kg_m2, use_drag, max_step_s) + return _propagate_python(sample_times_s, state0, ballistic_coef_kg_m2, use_drag, max_step_s) diff --git a/src/mcdowell_arc/fast.py b/src/mcdowell_arc/fast.py new file mode 100644 index 0000000..e1a664e --- /dev/null +++ b/src/mcdowell_arc/fast.py @@ -0,0 +1,88 @@ +"""Optional Rust acceleration bridge for McDowell Arc. + +The package is fully usable without Rust. When the `mcdowell_arc_core` extension +is installed, selected deterministic kernels can be delegated to it. +""" + +from __future__ import annotations + +from functools import lru_cache +from importlib import import_module +from types import ModuleType +from typing import Literal + +import numpy as np + +BackendName = Literal["auto", "python", "rust"] + + +class BackendUnavailableError(RuntimeError): + """Raised when the Rust backend is explicitly requested but unavailable.""" + + +@lru_cache(maxsize=1) +def _load_core() -> ModuleType | None: + try: + return import_module("mcdowell_arc_core") + except ModuleNotFoundError: + return None + + +def core_available() -> bool: + """Return True when the optional Rust extension can be imported.""" + return _load_core() is not None + + +def core_version() -> str | None: + """Return the Rust core version, if the extension is installed.""" + core = _load_core() + return None if core is None else str(getattr(core, "__version__", "unknown")) + + +def require_core() -> ModuleType: + """Return the Rust extension or raise a clear installation error.""" + core = _load_core() + if core is None: + raise BackendUnavailableError( + "Rust backend requested, but mcdowell_arc_core is not installed. " + "Run: maturin develop --manifest-path crates/mcdowell-arc-core/Cargo.toml --release" + ) + return core + + +def validate_backend(backend: str) -> BackendName: + """Validate a backend selector from CLI/API input.""" + allowed = {"auto", "python", "rust"} + if backend not in allowed: + raise ValueError(f"backend must be one of {sorted(allowed)}") + return backend # type: ignore[return-value] + + +def resolve_backend(backend: str) -> Literal["python", "rust"]: + """Resolve `auto` to the concrete backend that will be used.""" + selected = validate_backend(backend) + if selected == "python": + return "python" + if selected == "rust": + require_core() + return "rust" + return "rust" if core_available() else "python" + + +def propagate_core( + sample_times_s: np.ndarray, + state0: np.ndarray, + ballistic_coef_kg_m2: float, + use_drag: bool, + max_step_s: float, +) -> np.ndarray: + """Call Rust propagation and return a NumPy array.""" + core = require_core() + rows = core.propagate( + np.asarray(sample_times_s, dtype=float).tolist(), + np.asarray(state0, dtype=float).tolist(), + float(ballistic_coef_kg_m2), + bool(use_drag), + float(max_step_s), + ) + return np.asarray(rows, dtype=float) diff --git a/src/mcdowell_arc/fit.py b/src/mcdowell_arc/fit.py index 1168b83..3793c5d 100644 --- a/src/mcdowell_arc/fit.py +++ b/src/mcdowell_arc/fit.py @@ -3,6 +3,7 @@ from __future__ import annotations from dataclasses import dataclass +from pathlib import Path import numpy as np import pandas as pd @@ -10,18 +11,72 @@ from .constants import R_EARTH_KM from .dynamics import propagate +from .fast import validate_backend from .frames import earth_rotation_velocity_km_s from .orbit import OrbitSummary, summarize_orbit @dataclass(frozen=True) class Observations: + """Webcast-style trajectory observations. + + The MVP observation format intentionally contains only time, altitude, and + Earth-relative speed. That makes it useful for public webcast telemetry, but + it is underconstrained compared with real tracking geometry. + """ + t_s: np.ndarray altitude_km: np.ndarray speed_km_s: np.ndarray sigma_altitude_km: np.ndarray sigma_speed_km_s: np.ndarray + def validated(self) -> "Observations": + arrays = { + "t_s": np.asarray(self.t_s, dtype=float), + "altitude_km": np.asarray(self.altitude_km, dtype=float), + "speed_km_s": np.asarray(self.speed_km_s, dtype=float), + "sigma_altitude_km": np.asarray(self.sigma_altitude_km, dtype=float), + "sigma_speed_km_s": np.asarray(self.sigma_speed_km_s, dtype=float), + } + length = len(arrays["t_s"]) + if length < 3: + raise ValueError("At least three observations are required for fitting.") + for name, values in arrays.items(): + if values.ndim != 1: + raise ValueError(f"{name} must be a 1-D array.") + if len(values) != length: + raise ValueError(f"{name} length must match t_s length.") + if not np.all(np.isfinite(values)): + raise ValueError(f"{name} must contain only finite values.") + if np.any(np.diff(arrays["t_s"]) <= 0): + raise ValueError("t_s must be strictly increasing after sorting.") + if np.any(arrays["altitude_km"] < -1.0): + raise ValueError("altitude_km contains values below the accepted guard band.") + if np.any(arrays["speed_km_s"] <= 0.0): + raise ValueError("speed_km_s must be positive.") + if np.any(arrays["sigma_altitude_km"] <= 0.0) or np.any(arrays["sigma_speed_km_s"] <= 0.0): + raise ValueError("uncertainty columns must be positive.") + return Observations(**arrays) + + +@dataclass(frozen=True) +class FitOptions: + """Controls for deterministic fitting and propagation.""" + + use_drag: bool = True + backend: str = "auto" + max_step_s: float = 1.0 + max_nfev: int = 35 + + def validated(self) -> "FitOptions": + validate_backend(self.backend) + if self.max_step_s <= 0 or not np.isfinite(self.max_step_s): + raise ValueError("max_step_s must be positive and finite.") + if self.max_nfev <= 0: + raise ValueError("max_nfev must be positive.") + return self + @dataclass(frozen=True) class FitResult: @@ -34,8 +89,15 @@ class FitResult: message: str -def load_observations_csv(path: str) -> Observations: - """Load webcast-style observations from CSV.""" +def _numeric_series(df: pd.DataFrame, column: str) -> np.ndarray: + try: + return pd.to_numeric(df[column], errors="raise").to_numpy(float) + except Exception as exc: + raise ValueError(f"CSV column {column!r} must be numeric.") from exc + + +def load_observations_csv(path: str | Path) -> Observations: + """Load webcast-style observations from CSV and validate them.""" df = pd.read_csv(path) required = {"t_s", "altitude_km", "speed_km_s"} missing = sorted(required - set(df.columns)) @@ -43,20 +105,21 @@ def load_observations_csv(path: str) -> Observations: raise ValueError(f"Missing required CSV columns: {', '.join(missing)}") df = df.sort_values("t_s").reset_index(drop=True) - sigma_alt = df["sigma_altitude_km"].to_numpy(float) if "sigma_altitude_km" in df else np.full(len(df), 3.0) - sigma_speed = df["sigma_speed_km_s"].to_numpy(float) if "sigma_speed_km_s" in df else np.full(len(df), 0.05) + sigma_alt = _numeric_series(df, "sigma_altitude_km") if "sigma_altitude_km" in df else np.full(len(df), 3.0) + sigma_speed = _numeric_series(df, "sigma_speed_km_s") if "sigma_speed_km_s" in df else np.full(len(df), 0.05) sigma_alt = np.where(sigma_alt > 0, sigma_alt, 3.0) sigma_speed = np.where(sigma_speed > 0, sigma_speed, 0.05) - t = df["t_s"].to_numpy(float) + t = _numeric_series(df, "t_s") t = t - t[0] - return Observations( + obs = Observations( t_s=t, - altitude_km=df["altitude_km"].to_numpy(float), - speed_km_s=df["speed_km_s"].to_numpy(float), + altitude_km=_numeric_series(df, "altitude_km"), + speed_km_s=_numeric_series(df, "speed_km_s"), sigma_altitude_km=sigma_alt, sigma_speed_km_s=sigma_speed, ) + return obs.validated() def _params_to_state(params: np.ndarray, obs: Observations) -> tuple[np.ndarray, float]: @@ -72,6 +135,7 @@ def _params_to_state(params: np.ndarray, obs: Observations) -> tuple[np.ndarray, def observables_from_states(states: np.ndarray) -> tuple[np.ndarray, np.ndarray]: """Compute altitude and Earth-relative speed from propagated ECI states.""" + states = np.asarray(states, dtype=float) r = states[:, :3] v = states[:, 3:] altitude = np.linalg.norm(r, axis=1) - R_EARTH_KM @@ -80,11 +144,18 @@ def observables_from_states(states: np.ndarray) -> tuple[np.ndarray, np.ndarray] return altitude, rel_speed -def residuals(params: np.ndarray, obs: Observations, use_drag: bool = True) -> np.ndarray: +def residuals(params: np.ndarray, obs: Observations, options: FitOptions) -> np.ndarray: """Weighted residual vector for least-squares fitting.""" state0, beta = _params_to_state(params, obs) try: - states = propagate(obs.t_s, state0, beta, use_drag=use_drag) + states = propagate( + obs.t_s, + state0, + beta, + use_drag=options.use_drag, + max_step_s=options.max_step_s, + backend=options.backend, + ) alt_model, speed_model = observables_from_states(states) except Exception: # Penalize impossible trial states without crashing the optimizer. @@ -100,6 +171,7 @@ def residuals(params: np.ndarray, obs: Observations, use_drag: bool = True) -> n def initial_guess(obs: Observations) -> np.ndarray: """Build a simple first guess near the first observation.""" + obs = obs.validated() if len(obs.t_s) >= 2 and obs.t_s[-1] > obs.t_s[0]: vr_guess = np.gradient(obs.altitude_km, obs.t_s)[0] else: @@ -117,8 +189,12 @@ def fit_observations( use_drag: bool = True, initial_params: np.ndarray | None = None, max_nfev: int = 35, + backend: str = "auto", + max_step_s: float = 1.0, ) -> FitResult: """Fit a reduced 2-D drag-aware ECI trajectory.""" + obs = obs.validated() + options = FitOptions(use_drag=use_drag, backend=backend, max_step_s=max_step_s, max_nfev=max_nfev).validated() x0 = initial_guess(obs) if initial_params is None else np.asarray(initial_params, dtype=float) bounds = ( np.array([-100.0, -5.0, 0.1, np.log(10.0)]), @@ -128,15 +204,15 @@ def fit_observations( residuals, x0=x0, bounds=bounds, - args=(obs, use_drag), + args=(obs, options), loss="soft_l1", f_scale=2.0, - max_nfev=max_nfev, + max_nfev=options.max_nfev, x_scale=np.array([10.0, 1.0, 1.0, 2.0]), ) state0, beta = _params_to_state(opt.x, obs) orbit = summarize_orbit(state0) - res = residuals(opt.x, obs, use_drag=use_drag) + res = residuals(opt.x, obs, options) return FitResult( state0=state0, ballistic_coef_kg_m2=beta, @@ -154,10 +230,15 @@ def monte_carlo( atmosphere_threshold_km: float = 100.0, seed: int = 7, use_drag: bool = True, + backend: str = "auto", + max_step_s: float = 1.0, ) -> dict: """Perturb observations by their sigmas and refit many times.""" + obs = obs.validated() + if samples <= 0: + raise ValueError("samples must be positive.") rng = np.random.default_rng(seed) - nominal = fit_observations(obs, use_drag=use_drag) + nominal = fit_observations(obs, use_drag=use_drag, backend=backend, max_step_s=max_step_s) perigees: list[float] = [] apogees: list[float] = [] betas: list[float] = [] @@ -170,15 +251,24 @@ def monte_carlo( speed_km_s=obs.speed_km_s + rng.normal(0.0, obs.sigma_speed_km_s), sigma_altitude_km=obs.sigma_altitude_km.copy(), sigma_speed_km_s=obs.sigma_speed_km_s.copy(), - ) + ).validated() try: - warm_start = np.array([ - nominal.state0[0] - (R_EARTH_KM + perturbed.altitude_km[0]), - nominal.state0[3], - nominal.state0[4], - np.log(nominal.ballistic_coef_kg_m2), - ]) - fit = fit_observations(perturbed, use_drag=use_drag, initial_params=warm_start, max_nfev=18) + warm_start = np.array( + [ + nominal.state0[0] - (R_EARTH_KM + perturbed.altitude_km[0]), + nominal.state0[3], + nominal.state0[4], + np.log(nominal.ballistic_coef_kg_m2), + ] + ) + fit = fit_observations( + perturbed, + use_drag=use_drag, + initial_params=warm_start, + max_nfev=18, + backend=backend, + max_step_s=max_step_s, + ) perigees.append(fit.orbit.perigee_km) if fit.orbit.apogee_km is not None: apogees.append(fit.orbit.apogee_km) diff --git a/tests/test_atmosphere.py b/tests/test_atmosphere.py new file mode 100644 index 0000000..d4ab11a --- /dev/null +++ b/tests/test_atmosphere.py @@ -0,0 +1,12 @@ +import numpy as np + +from mcdowell_arc.atmosphere import density_kg_m3 + + +def test_density_monotonically_decreases_in_upper_atmosphere(): + samples = density_kg_m3(np.array([100.0, 150.0, 250.0, 400.0])) + assert np.all(np.diff(samples) < 0.0) + + +def test_density_scalar_return_type(): + assert isinstance(density_kg_m3(100.0), float) diff --git a/tests/test_dynamics.py b/tests/test_dynamics.py new file mode 100644 index 0000000..f101283 --- /dev/null +++ b/tests/test_dynamics.py @@ -0,0 +1,32 @@ +import numpy as np +import pytest + +from mcdowell_arc.constants import MU_EARTH_KM3_S2, R_EARTH_KM +from mcdowell_arc.dynamics import propagate +from mcdowell_arc.fast import core_available + + +def test_propagate_python_returns_expected_shape(): + r = R_EARTH_KM + 400.0 + v = np.sqrt(MU_EARTH_KM3_S2 / r) + states = propagate(np.array([0.0, 10.0, 20.0]), np.array([r, 0, 0, 0, v, 0]), 250.0, use_drag=False, backend="python") + assert states.shape == (3, 6) + assert np.all(np.isfinite(states)) + + +def test_auto_backend_falls_back_without_rust(): + r = R_EARTH_KM + 400.0 + v = np.sqrt(MU_EARTH_KM3_S2 / r) + states = propagate(np.array([0.0, 5.0]), np.array([r, 0, 0, 0, v, 0]), 250.0, use_drag=False, backend="auto") + assert states.shape == (2, 6) + + +@pytest.mark.skipif(not core_available(), reason="Rust extension is not installed") +def test_rust_backend_matches_python_for_short_no_drag_propagation(): + r = R_EARTH_KM + 400.0 + v = np.sqrt(MU_EARTH_KM3_S2 / r) + times = np.array([0.0, 10.0]) + state0 = np.array([r, 0, 0, 0, v, 0]) + py_states = propagate(times, state0, 250.0, use_drag=False, backend="python") + rs_states = propagate(times, state0, 250.0, use_drag=False, backend="rust") + assert np.allclose(py_states, rs_states, rtol=1e-11, atol=1e-11) diff --git a/tests/test_fit_cli.py b/tests/test_fit_cli.py new file mode 100644 index 0000000..8020bcf --- /dev/null +++ b/tests/test_fit_cli.py @@ -0,0 +1,38 @@ +import json +from pathlib import Path + +import pytest + +from mcdowell_arc.cli import main +from mcdowell_arc.fit import load_observations_csv + +ROOT = Path(__file__).resolve().parents[1] +SAMPLE = ROOT / "examples" / "sample_webcast.csv" + + +def test_load_observations_csv_validates_sample(): + obs = load_observations_csv(SAMPLE) + assert len(obs.t_s) == 11 + assert obs.t_s[0] == 0.0 + + +def test_cli_doctor_outputs_json(capsys): + assert main(["doctor"]) == 0 + payload = json.loads(capsys.readouterr().out) + assert "mcdowell_arc_version" in payload + assert "rust_core_available" in payload + + +def test_cli_fit_smoke(capsys): + assert main(["fit", str(SAMPLE), "--monte-carlo", "2", "--backend", "python"]) == 0 + payload = json.loads(capsys.readouterr().out) + assert payload["backend"] == "python" + assert payload["nominal"]["optimizer_success"] is True + assert payload["monte_carlo"]["samples_requested"] == 2 + + +def test_bad_csv_missing_required_column(tmp_path): + path = tmp_path / "bad.csv" + path.write_text("t_s,altitude_km\n0,100\n1,101\n2,102\n") + with pytest.raises(ValueError, match="Missing required CSV columns"): + load_observations_csv(path) From 352436481f1fb760daf9575e2264c54646b2202d Mon Sep 17 00:00:00 2001 From: Chase Bryan <281080815+Bryforge@users.noreply.github.com> Date: Fri, 22 May 2026 22:41:21 -0500 Subject: [PATCH 2/2] Fix CI lint and Rust extension check --- .github/workflows/ci.yml | 4 ++-- src/mcdowell_arc/fit.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b5d2ee7..cac8411 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -32,8 +32,8 @@ jobs: with: python-version: "3.12" - uses: actions-rust-lang/setup-rust-toolchain@v1 - - name: Rust tests - run: cargo test --workspace + - name: Rust check + run: cargo check --workspace - name: Build Rust Python core run: | python -m pip install --upgrade pip maturin diff --git a/src/mcdowell_arc/fit.py b/src/mcdowell_arc/fit.py index 3793c5d..ce1c3dc 100644 --- a/src/mcdowell_arc/fit.py +++ b/src/mcdowell_arc/fit.py @@ -31,7 +31,7 @@ class Observations: sigma_altitude_km: np.ndarray sigma_speed_km_s: np.ndarray - def validated(self) -> "Observations": + def validated(self) -> Observations: arrays = { "t_s": np.asarray(self.t_s, dtype=float), "altitude_km": np.asarray(self.altitude_km, dtype=float), @@ -69,7 +69,7 @@ class FitOptions: max_step_s: float = 1.0 max_nfev: int = 35 - def validated(self) -> "FitOptions": + def validated(self) -> FitOptions: validate_backend(self.backend) if self.max_step_s <= 0 or not np.isfinite(self.max_step_s): raise ValueError("max_step_s must be positive and finite.")