spb/satelliteindex

1"""Orbital validation against known objects (ISS fixture) and classification rules."""2from __future__ import annotations34import json5import math6from datetime import UTC, datetime, timedelta7from pathlib import Path89import numpy as np10import pytest11from sgp4 import omm12from sgp4.api import Satrec1314from satelliteindex.orbital.elements import classify_orbit, derived_geometry15from satelliteindex.orbital.propagate import BatchPropagator, Elements, ground_track, jd_fr, propagate_one, teme_to_geodetic1617FIX = Path(__file__).parent / "fixtures"181920def iss_elements() -> tuple[Elements, dict]:21    d = json.loads((FIX / "celestrak_stations.json").read_text())[0]22    assert d["NORAD_CAT_ID"] == 2554423    el = Elements(satellite_id="sat_test", norad_id=25544, epoch=datetime.fromisoformat(d["EPOCH"]).replace(tzinfo=UTC), mean_motion=d["MEAN_MOTION"],24                  eccentricity=d["ECCENTRICITY"], inclination=d["INCLINATION"], raan=d["RA_OF_ASC_NODE"], arg_of_perigee=d["ARG_OF_PERICENTER"],25                  mean_anomaly=d["MEAN_ANOMALY"], bstar=d["BSTAR"], mean_motion_dot=d["MEAN_MOTION_DOT"], mean_motion_ddot=d["MEAN_MOTION_DDOT"])26    return el, d272829def test_iss_geometry_and_class():30    el, _ = iss_elements()31    g = derived_geometry(el.mean_motion, el.eccentricity)32    assert 400 < g["perigee_km"] < 440 and 400 < g["apogee_km"] < 44033    assert 92 < g["period_minutes"] < 9434    assert classify_orbit(period_minutes=g["period_minutes"], eccentricity=el.eccentricity, inclination_deg=el.inclination, apogee_km=g["apogee_km"], perigee_km=g["perigee_km"]) == "LEO"353637def test_iss_position_matches_reference_sgp4():38    el, d = iss_elements()39    t = el.epoch + timedelta(minutes=37)40    mine = propagate_one(el, t)41    ref = Satrec()42    omm.initialize(ref, {k: str(v) for k, v in d.items()})43    jd, fr = jd_fr(t)44    e, r, v = ref.sgp4(jd, fr)45    assert e == 0 and mine["error"] is None46    assert np.linalg.norm(np.array(r) - np.array(mine["position_teme_km"])) < 0.5  # km47    assert 400 < mine["altitude_km"] < 44048    assert 7.5 < mine["velocity_km_s"] < 7.849    assert -52 < mine["lat"] < 52  # ISS inclination bound505152def test_ground_track_shape_and_wrap():53    el, _ = iss_elements()54    pts = ground_track(el, el.epoch, minutes_before=45, minutes_after=90, step_s=60)55    assert len(pts) == 13656    assert all(-180 <= p["lon"] <= 180 for p in pts)57    assert sum(1 for p in pts if p["future"]) == 91585960def test_batch_matches_single():61    el, _ = iss_elements()62    t = el.epoch + timedelta(hours=1)63    b = BatchPropagator([el, el])64    p = b.positions(t)65    s = propagate_one(el, t)66    assert p["ok"].all()67    assert abs(p["lat"][0] - s["lat"]) < 1e-6 and abs(p["alt"][1] - s["altitude_km"]) < 1e-6686970def test_geodetic_conversion_equator_and_pole():71    jd = np.array(2460000.5)72    r_eq = np.array([6378.137 + 500.0, 0.0, 0.0])73    lat, lon, alt = teme_to_geodetic(r_eq, jd)74    assert abs(lat) < 1e-6 and abs(alt - 500.0) < 1e-375    r_pole = np.array([0.0, 0.0, 6356.7523 + 800.0])76    lat, lon, alt = teme_to_geodetic(r_pole, jd)77    assert abs(lat - 90) < 1e-4 and abs(alt - 800.0) < 0.01787980@pytest.mark.parametrize("mm,e,i,expected", [81    (1.00271, 0.0002, 0.05, "GEO"),      # geostationary82    (2.00565, 0.01, 55.0, "MEO"),        # GPS83    (15.5, 0.0005, 51.6, "LEO"),         # ISS84    (2.0, 0.74, 63.4, "HEO"),            # Molniya85    (1.00271, 0.0002, 7.0, "GEO"),       # inclined geosynchronous86])87def test_classification_table(mm, e, i, expected):88    g = derived_geometry(mm, e)89    assert classify_orbit(period_minutes=g["period_minutes"], eccentricity=e, inclination_deg=i, apogee_km=g["apogee_km"], perigee_km=g["perigee_km"]) == expected909192def test_semi_major_axis_gps():93    g = derived_geometry(2.00565, 0.0)94    assert math.isclose(g["semi_major_axis_km"], 26560, rel_tol=0.002)95