"""Orbital validation against known objects (ISS fixture) and classification rules."""
from __future__ import annotations

import json
import math
from datetime import UTC, datetime, timedelta
from pathlib import Path

import numpy as np
import pytest
from sgp4 import omm
from sgp4.api import Satrec

from satelliteindex.orbital.elements import classify_orbit, derived_geometry
from satelliteindex.orbital.propagate import BatchPropagator, Elements, ground_track, jd_fr, propagate_one, teme_to_geodetic

FIX = Path(__file__).parent / "fixtures"


def iss_elements() -> tuple[Elements, dict]:
    d = json.loads((FIX / "celestrak_stations.json").read_text())[0]
    assert d["NORAD_CAT_ID"] == 25544
    el = Elements(satellite_id="sat_test", norad_id=25544, epoch=datetime.fromisoformat(d["EPOCH"]).replace(tzinfo=UTC), mean_motion=d["MEAN_MOTION"],
                  eccentricity=d["ECCENTRICITY"], inclination=d["INCLINATION"], raan=d["RA_OF_ASC_NODE"], arg_of_perigee=d["ARG_OF_PERICENTER"],
                  mean_anomaly=d["MEAN_ANOMALY"], bstar=d["BSTAR"], mean_motion_dot=d["MEAN_MOTION_DOT"], mean_motion_ddot=d["MEAN_MOTION_DDOT"])
    return el, d


def test_iss_geometry_and_class():
    el, _ = iss_elements()
    g = derived_geometry(el.mean_motion, el.eccentricity)
    assert 400 < g["perigee_km"] < 440 and 400 < g["apogee_km"] < 440
    assert 92 < g["period_minutes"] < 94
    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"


def test_iss_position_matches_reference_sgp4():
    el, d = iss_elements()
    t = el.epoch + timedelta(minutes=37)
    mine = propagate_one(el, t)
    ref = Satrec()
    omm.initialize(ref, {k: str(v) for k, v in d.items()})
    jd, fr = jd_fr(t)
    e, r, v = ref.sgp4(jd, fr)
    assert e == 0 and mine["error"] is None
    assert np.linalg.norm(np.array(r) - np.array(mine["position_teme_km"])) < 0.5  # km
    assert 400 < mine["altitude_km"] < 440
    assert 7.5 < mine["velocity_km_s"] < 7.8
    assert -52 < mine["lat"] < 52  # ISS inclination bound


def test_ground_track_shape_and_wrap():
    el, _ = iss_elements()
    pts = ground_track(el, el.epoch, minutes_before=45, minutes_after=90, step_s=60)
    assert len(pts) == 136
    assert all(-180 <= p["lon"] <= 180 for p in pts)
    assert sum(1 for p in pts if p["future"]) == 91


def test_batch_matches_single():
    el, _ = iss_elements()
    t = el.epoch + timedelta(hours=1)
    b = BatchPropagator([el, el])
    p = b.positions(t)
    s = propagate_one(el, t)
    assert p["ok"].all()
    assert abs(p["lat"][0] - s["lat"]) < 1e-6 and abs(p["alt"][1] - s["altitude_km"]) < 1e-6


def test_geodetic_conversion_equator_and_pole():
    jd = np.array(2460000.5)
    r_eq = np.array([6378.137 + 500.0, 0.0, 0.0])
    lat, lon, alt = teme_to_geodetic(r_eq, jd)
    assert abs(lat) < 1e-6 and abs(alt - 500.0) < 1e-3
    r_pole = np.array([0.0, 0.0, 6356.7523 + 800.0])
    lat, lon, alt = teme_to_geodetic(r_pole, jd)
    assert abs(lat - 90) < 1e-4 and abs(alt - 800.0) < 0.01


@pytest.mark.parametrize("mm,e,i,expected", [
    (1.00271, 0.0002, 0.05, "GEO"),      # geostationary
    (2.00565, 0.01, 55.0, "MEO"),        # GPS
    (15.5, 0.0005, 51.6, "LEO"),         # ISS
    (2.0, 0.74, 63.4, "HEO"),            # Molniya
    (1.00271, 0.0002, 7.0, "GEO"),       # inclined geosynchronous
])
def test_classification_table(mm, e, i, expected):
    g = derived_geometry(mm, e)
    assert classify_orbit(period_minutes=g["period_minutes"], eccentricity=e, inclination_deg=i, apogee_km=g["apogee_km"], perigee_km=g["perigee_km"]) == expected


def test_semi_major_axis_gps():
    g = derived_geometry(2.00565, 0.0)
    assert math.isclose(g["semi_major_axis_km"], 26560, rel_tol=0.002)