Swapi proton moments

docs/swapi/figure_src/build_validation_table.py

@@ -0,0 +1,170 @@
+#!/usr/bin/env python3
+"""Generate the integrator-validation table in `docs/swapi/solar-wind-moments.md`."""
+
+import sys
+import types
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[3]))
+
+import numpy as np
+import pandas as pd
+
+from imap_l3_processing.constants import (
+    EV_TO_KELVIN,
+    PROTON_MASS_KG,
+    PROTON_MASS_PER_CHARGE_M_P_PER_E,
+)
+from imap_l3_processing.swapi.l3a.science.solar_wind.forward_model import (
+    calculate_integral,
+)
+from imap_l3_processing.swapi.l3a.science.solar_wind.params import SolarWindParams
+
+from figure_utils import (
+    REPO_ROOT,
+    bulk_velocity_rtn_from_swapi_angles,
+    load_swapi_response,
+    peak_esa_voltage_for_proton_bulk_speed,
+    run_parallel_map,
+)
+
+_worker_state: types.SimpleNamespace | None = None
+
+_REFERENCE_INTEGRALS_PATH = (
+    REPO_ROOT / "tests" / "test_data" / "swapi" / "reference_integrals.csv"
+)
+_DOC_PATH = REPO_ROOT / "docs" / "swapi" / "solar-wind-moments.md"
+_TABLE_BEGIN = "<!-- BEGIN: validation_table"
+_TABLE_END = "<!-- END: validation_table -->"
+
+_BANDS = [
+    (r"$\lt 0.1$", 0.0, 0.1),
+    ("$0.1$ – $1$", 0.1, 1.0),
+    ("$1$ – $10$", 1.0, 10.0),
+    ("$10$ – $10^2$", 10.0, 100.0),
+    ("$10^2$ – $10^3$", 100.0, 1000.0),
+    ("$10^3$ – $10^4$", 1000.0, 1e4),
+    ("$10^4$ – $10^5$", 1e4, 1e5),
+    ("$\\geq 10^5$", 1e5, np.inf),
+]
+
+
+def _build_table(rel: np.ndarray, refs_v: np.ndarray) -> str:
+    lines = [
+        "| Reference (Hz)  |     N |  Median |     95% |     99% |     Max |",
+        "|-----------------|-------|---------|---------|---------|---------|",
+    ]
+    for label, lo, hi in _BANDS:
+        mask = (refs_v >= lo) & (refs_v < hi)
+        n = int(mask.sum())
+        if n == 0:
+            lines.append(
+                f"| {label:<15s} | {n:>5d} |       — |       — |       — |       — |"
+            )
+            continue
+        band = rel[mask]
+        med = np.median(band) * 100
+        p95 = np.percentile(band, 95) * 100
+        p99 = np.percentile(band, 99) * 100
+        mx = band.max() * 100
+        lines.append(
+            f"| {label:<15s} | {n:>5d} | {med:>6.2f}% | {p95:>6.2f}% | {p99:>6.2f}% | {mx:>6.2f}% |"
+        )
+    return "\n".join(lines)
+
+
+def _update_doc(table_md: str) -> None:
+    text = _DOC_PATH.read_text()
+    begin_idx = text.find(_TABLE_BEGIN)
+    end_idx = text.find(_TABLE_END)
+    if begin_idx < 0 or end_idx < 0 or end_idx <= begin_idx:
+        raise RuntimeError(
+            f"Could not find '{_TABLE_BEGIN}' / '{_TABLE_END}' markers in {_DOC_PATH}"
+        )
+    begin_line_end = text.find("\n", begin_idx) + 1
+    new_text = text[:begin_line_end] + table_md + "\n" + text[end_idx:]
+    _DOC_PATH.write_text(new_text)
+
+
+def _initialize_worker_state(
+    rows: list[tuple[float, float, float, float, float]],
+    peak_voltages: np.ndarray,
+    swapi_response,
+) -> None:
+    global _worker_state
+    _worker_state = types.SimpleNamespace(
+        rows=rows,
+        peak_voltages=peak_voltages,
+        swapi_response=swapi_response,
+        rotation_matrix=np.eye(3),
+    )
+
+
+def _process_one(i: int) -> float:
+    state = _worker_state
+    bulk_speed, azimuth, elevation, density, temperature_ev = state.rows[i]
+    sw = SolarWindParams(
+        density=density,
+        bulk_velocity_rtn=bulk_velocity_rtn_from_swapi_angles(
+            bulk_speed, azimuth, elevation
+        ),
+        temperature=temperature_ev * EV_TO_KELVIN,
+        mass=PROTON_MASS_KG,
+    )
+    response_grid = state.swapi_response.get_response_grid(
+        state.peak_voltages[i], PROTON_MASS_PER_CHARGE_M_P_PER_E
+    )
+    return calculate_integral(sw, response_grid, state.rotation_matrix)[0]
+
+
+def main():
+    print("Loading calibration data...")
+    swapi_response = load_swapi_response()
+    df = pd.read_csv(_REFERENCE_INTEGRALS_PATH)
+    references = df["integral"].to_numpy()
+    rows = [
+        (
+            float(row.bulk_speed),
+            float(row.bulk_azimuth),
+            float(row.bulk_elevation),
+            float(row.density),
+            float(row.temperature_ev),
+        )
+        for row in df.itertuples(index=False)
+    ]
+
+    print(f"Warming passband cache for {len(df)} rows...")
+    peak_voltages = np.array(
+        [peak_esa_voltage_for_proton_bulk_speed(r[0]) for r in rows]
+    )
+    swapi_response.warm_cache(peak_voltages)
+    for unique_voltage in np.unique(peak_voltages):
+        swapi_response.get_response_grid(
+            float(unique_voltage), PROTON_MASS_PER_CHARGE_M_P_PER_E
+        )
+
+    _initialize_worker_state(rows, peak_voltages, swapi_response)
+    optimized = np.array(
+        run_parallel_map(_process_one, len(df), desc="integrals", chunksize=64)
+    )
+
+    valid = references > 0
+    refs_v = references[valid]
+    rel = np.abs(optimized[valid] / refs_v - 1.0)
+    # Cases where both integrators round to ~0 give meaningless ratios; clamp.
+    rel = np.minimum(rel, 1.0)
+
+    print(f"\nMax |ratio - 1|:    {rel.max():.2%}")
+    print(f"Median |ratio - 1|: {np.median(rel):.2%}")
+    print(f"95th pct:           {np.percentile(rel, 95):.2%}")
+    print(f"99th pct:           {np.percentile(rel, 99):.2%}")
+
+    table_md = _build_table(rel, refs_v)
+    print()
+    print(table_md)
+    _update_doc(table_md)
+    print(f"\nUpdated {_DOC_PATH.relative_to(REPO_ROOT)}")
+
+
+if __name__ == "__main__":
+    main()

docs/swapi/figure_src/figure_utils.py

@@ -0,0 +1,163 @@
+"""Shared helpers for the SWAPI documentation figures."""
+
+import multiprocessing
+import os
+import sys
+import time
+from pathlib import Path
+from typing import Callable, TypeVar
+
+import numpy as np
+
+REPO_ROOT = Path(__file__).resolve().parents[3]
+sys.path.insert(0, str(REPO_ROOT))
+
+from imap_l3_processing.constants import (
+    METERS_PER_KILOMETER,
+    PROTON_CHARGE_COULOMBS,
+    PROTON_MASS_KG,
+)
+from imap_l3_processing.swapi.constants import SWAPI_K_FACTOR
+from imap_l3_processing.swapi.response.swapi_response import SwapiResponse
+
+FIGURES_DIR = REPO_ROOT / "docs" / "swapi" / "figures"
+
+_INSTRUMENT_DATA_DIR = REPO_ROOT / "instrument_team_data" / "swapi"
+
+
+def load_swapi_response() -> SwapiResponse:
+    return SwapiResponse.from_files(
+        _INSTRUMENT_DATA_DIR / "imap_swapi_azimuthal-transmission_20260425_v001.csv",
+        _INSTRUMENT_DATA_DIR / "imap_swapi_central-effective-area_20260425_v001.csv",
+        _INSTRUMENT_DATA_DIR / "imap_swapi_passband-fit-coefficients_20260425_v001.csv",
+    )
+
+
+SWEEP_DURATION_S = 12.0
+BINS_PER_SWEEP = 72
+BIN_PERIOD_S = SWEEP_DURATION_S / BINS_PER_SWEEP
+SPIN_PERIOD_S = 15.13
+
+COARSE_BIN_INDICES_IN_SWEEP = np.arange(1, 63)
+
+COARSE_SWEEP_VOLTAGES_MEAN_V = np.array(
+    [
+        9895.52, 9088.69, 8348.80, 7667.55, 7042.16, 6469.31, 5941.77, 5457.31,
+        5013.22, 4603.65, 4230.77, 3886.92, 3569.16, 3278.72, 3011.13, 2766.25,
+        2539.54, 2333.83, 2144.24, 1969.31, 1808.74, 1660.86, 1525.75, 1401.82,
+        1287.58, 1182.24, 1085.15, 995.55, 914.31, 839.94, 771.70, 709.46,
+        651.59, 598.47, 549.91, 505.12, 463.89, 425.92, 391.18, 359.35,
+        329.94, 303.02, 278.25, 255.55, 234.77, 215.61, 197.95, 181.82,
+        167.04, 153.46, 140.91, 129.50, 118.91, 109.20, 100.30, 92.11,
+        84.61, 77.73, 71.40, 65.59, 60.23, 55.34,
+    ]
+)
+
+# SWAPI→RTN at the first-sweep midpoint of a real SPICE attitude near 2026-01-01;
+# spin axis (+Y column) sits ~4° off −R̂_RTN. Per-bin matrices are built by
+# spinning this about its own +Y at the nominal IMAP spin period.
+ANCHOR_ROTATION_SWAPI_TO_RTN = np.array(
+    [
+        [+0.0705, +0.9157, +0.3955],
+        [-0.9968, +0.0792, -0.0057],
+        [-0.0365, -0.3939, +0.9184],
+    ]
+).T
+ANCHOR_TIME_S = 0.5 * SWEEP_DURATION_S
+# Sign reproduces independent SPICE-derived midpoints across a 5-sweep cycle.
+SPIN_OMEGA_RAD_S = -2.0 * np.pi / SPIN_PERIOD_S
+
+
+def compute_per_bin_rotation_matrices(
+    n_sweeps: int,
+    bin_indices_in_sweep: np.ndarray = COARSE_BIN_INDICES_IN_SWEEP,
+) -> np.ndarray:
+    """Synthetic per-bin SWAPI→RTN matrices: anchor spun about its spin axis.
+
+    Returns shape (n_sweeps · n_bins, 3, 3) in sweep-major, bin-minor order.
+    """
+    sweep_index = np.repeat(np.arange(n_sweeps), len(bin_indices_in_sweep))
+    bin_index = np.tile(bin_indices_in_sweep, n_sweeps)
+    sample_times_s = sweep_index * SWEEP_DURATION_S + bin_index * BIN_PERIOD_S
+
+    spin_axis = ANCHOR_ROTATION_SWAPI_TO_RTN[:, 1] / np.linalg.norm(
+        ANCHOR_ROTATION_SWAPI_TO_RTN[:, 1]
+    )
+    delta_phi = SPIN_OMEGA_RAD_S * (sample_times_s - ANCHOR_TIME_S)
+
+    ax, ay, az = spin_axis
+    K = np.array([[0, -az, ay], [az, 0, -ax], [-ay, ax, 0]])
+    sin_dp = np.sin(delta_phi)[:, None, None]
+    one_minus_cos = (1.0 - np.cos(delta_phi))[:, None, None]
+    rot = np.eye(3) + sin_dp * K + one_minus_cos * (K @ K)
+    return rot @ ANCHOR_ROTATION_SWAPI_TO_RTN
+
+
+def peak_esa_voltage_for_proton_bulk_speed(bulk_speed_km_s: float) -> float:
+    return float(
+        PROTON_MASS_KG
+        * (bulk_speed_km_s * METERS_PER_KILOMETER) ** 2
+        / (2 * SWAPI_K_FACTOR * PROTON_CHARGE_COULOMBS)
+    )
+
+
+_T = TypeVar("_T")
+
+
+def run_parallel_map(
+    process_one: Callable[[int], _T],
+    n_items: int,
+    *,
+    desc: str,
+    chunksize: int = 10,
+) -> list[_T]:
+    """Run `process_one(i)` for i in [0, n_items) across forked workers.
+
+    Module-level `_worker_state` is inherited by children via fork — callers
+    should populate it in the parent before invoking this helper.
+    """
+    if multiprocessing.get_start_method(allow_none=True) != "fork":
+        multiprocessing.set_start_method("fork", force=True)
+    n_workers = os.cpu_count() or 1
+    print(f"Running {n_items} {desc} across {n_workers} processes...")
+    start = time.perf_counter()
+    results: list[_T | None] = [None] * n_items
+    report_every = max(1, n_items // 20)
+    with multiprocessing.get_context("fork").Pool(processes=n_workers) as pool:
+        completed = 0
+        for index, value in pool.imap_unordered(
+            _IndexedCall(process_one), range(n_items), chunksize=chunksize
+        ):
+            results[index] = value
+            completed += 1
+            if completed % report_every == 0 or completed == n_items:
+                elapsed = time.perf_counter() - start
+                print(f"  {desc}: {completed}/{n_items} ({elapsed:.1f}s)")
+    print(f"  {desc} done in {time.perf_counter() - start:.1f}s.")
+    return results  # type: ignore[return-value]
+
+
+class _IndexedCall:
+    """Picklable wrapper that returns (index, func(index)) for ordered reassembly."""
+
+    def __init__(self, func: Callable[[int], _T]):
+        self._func = func
+
+    def __call__(self, index: int):
+        return index, self._func(index)
+
+
+def bulk_velocity_rtn_from_swapi_angles(
+    bulk_speed_km_s: float, azimuth_deg: float, elevation_deg: float
+) -> np.ndarray:
+    """Inverse of `bulk_angles_in_instrument_frame` for R = identity (RTN ≡ SWAPI)."""
+    az = np.radians(azimuth_deg)
+    el = np.radians(elevation_deg)
+    horizontal = bulk_speed_km_s * np.cos(el)
+    return np.array(
+        [
+            -horizontal * np.sin(az),
+            -horizontal * np.cos(az),
+            -bulk_speed_km_s * np.sin(el),
+        ]
+    )