from __future__ import annotations

import math

import pytest

from moira.coordinates import mat_mul, nutation_matrix_equatorial, precession_matrix_equatorial
from moira.julian import (
    apparent_sidereal_time_at,
    calendar_from_jd,
    earth_rotation_angle,
    greenwich_mean_sidereal_time,
    julian_day,
    ut_to_tt,
)
from moira.obliquity import nutation, true_obliquity
from moira.precession import general_precession_in_longitude, mean_obliquity_p03, precession_matrix

erfa = pytest.importorskip("erfa")

ARCSEC = 3600.0
PASS_THRESHOLD_ARCSEC = 0.001
_BCE_EPOCHS = [
    ("500 BCE (astronomical -499)", -499, 1538803.5),
    ("200 BCE (astronomical -199)", -199, 1648376.5),
]
TEST_EPOCHS = [
    *((label, jd, jd) for label, _astronomical_year, jd in _BCE_EPOCHS),
    ("J0000 (1 CE)", 1721045.5, 1721045.5),
    ("J1000.0", 2086308.0, 2086308.0),
    ("J1500.0", 2268923.5, 2268923.5),
    ("J1800.0", 2378496.5, 2378496.5),
    ("J1900.0", 2415020.5, 2415020.5),
    ("J2000.0", 2451545.0, 2451545.0),
    ("J2010.0", 2455196.5, 2455196.5),
    ("J2024.0", 2460310.5, 2460310.5),
    ("J2050.0", 2469807.5, 2469807.5),
    ("J2100.0", 2488069.5, 2488069.5),
]


def _split(jd: float) -> tuple[float, float]:
    return float(int(jd)), float(jd - int(jd))


def _wrapped_arcsec(a_deg: float, b_deg: float) -> float:
    delta_arcsec = abs(a_deg - b_deg) * ARCSEC
    if delta_arcsec > 180.0 * ARCSEC:
        delta_arcsec = 360.0 * ARCSEC - delta_arcsec
    return delta_arcsec


def _matrix_max_diff_arcsec(erfa_matrix, moira_matrix) -> float:
    max_diff = 0.0
    for i in range(3):
        for j in range(3):
            max_diff = max(max_diff, abs(float(erfa_matrix[i][j]) - moira_matrix[i][j]))
    return max_diff * (180.0 / math.pi) * ARCSEC


@pytest.mark.integration
@pytest.mark.parametrize(
    ("label", "astronomical_year", "expected_jd"),
    _BCE_EPOCHS,
    ids=[label.split(" (")[0] for label, _, _ in _BCE_EPOCHS],
)
def test_bce_epoch_calendar_identity(
    label: str,
    astronomical_year: int,
    expected_jd: float,
) -> None:
    erfa_d1, erfa_d2 = erfa.cal2jd(astronomical_year, 1, 1)

    assert julian_day(astronomical_year, 1, 1) == expected_jd
    assert float(erfa_d1 + erfa_d2) == expected_jd
    assert calendar_from_jd(expected_jd) == (astronomical_year, 1, 1, 0.0)
    assert (label, expected_jd, expected_jd) in TEST_EPOCHS


@pytest.mark.integration
@pytest.mark.parametrize(("label", "jd_tt", "jd_ut"), TEST_EPOCHS, ids=[label for label, _, _ in TEST_EPOCHS])
def test_gmst_matches_erfa(label: str, jd_tt: float, jd_ut: float) -> None:
    d1u, d2u = _split(jd_ut)
    d1t, d2t = _split(jd_tt)
    erfa_gmst = math.degrees(erfa.gmst06(d1u, d2u, d1t, d2t)) % 360.0
    moira_gmst = greenwich_mean_sidereal_time(jd_ut) % 360.0

    assert _wrapped_arcsec(erfa_gmst, moira_gmst) < PASS_THRESHOLD_ARCSEC


@pytest.mark.integration
@pytest.mark.parametrize(("label", "jd_tt", "jd_ut"), TEST_EPOCHS, ids=[label for label, _, _ in TEST_EPOCHS])
def test_era_matches_erfa(label: str, jd_tt: float, jd_ut: float) -> None:
    d1, d2 = _split(jd_ut)
    erfa_era = math.degrees(erfa.era00(d1, d2)) % 360.0
    moira_era = earth_rotation_angle(jd_ut) % 360.0

    assert _wrapped_arcsec(erfa_era, moira_era) < PASS_THRESHOLD_ARCSEC


@pytest.mark.integration
@pytest.mark.parametrize(("label", "jd_tt", "_jd_ut"), TEST_EPOCHS, ids=[label for label, _, _ in TEST_EPOCHS])
def test_mean_obliquity_matches_erfa(label: str, jd_tt: float, _jd_ut: float) -> None:
    d1, d2 = _split(jd_tt)
    erfa_eps = math.degrees(erfa.obl06(d1, d2))
    moira_eps = mean_obliquity_p03(jd_tt)

    assert abs(erfa_eps - moira_eps) * ARCSEC < PASS_THRESHOLD_ARCSEC


@pytest.mark.integration
@pytest.mark.parametrize("jd_tt", [
    *(jd_tt for _, jd_tt, _ in TEST_EPOCHS),
    2435553.5,  # Historical Lahiri origin epoch, already in TT.
    2461206.808439629,  # June 2026 ingress neighbourhood, already in TT.
    2451545.0 - 50 * 36525.0,
    2451545.0 + 50 * 36525.0,
])
def test_general_precession_matches_erfa_p06e(jd_tt: float) -> None:
    # eraP06e output 12 is p_A, general precession; output 15 is the
    # distinct Fukushima-Williams rotation angle psi. Authority:
    # https://github.com/liberfa/erfa/blob/v2.0.1/src/p06e.c
    # Compare a signed, unwrapped scalar over the stated +/-50-century
    # polynomial interval. The existing 0.001-arcsecond gate is unchanged.
    d1, d2 = _split(jd_tt)
    expected = math.degrees(erfa.p06e(d1, d2)[12])
    actual = general_precession_in_longitude(jd_tt)
    assert abs(actual - expected) * ARCSEC < PASS_THRESHOLD_ARCSEC


@pytest.mark.integration
@pytest.mark.parametrize(("label", "jd_tt", "_jd_ut"), TEST_EPOCHS, ids=[label for label, _, _ in TEST_EPOCHS])
def test_nutation_matches_erfa(label: str, jd_tt: float, _jd_ut: float) -> None:
    d1, d2 = _split(jd_tt)
    erfa_dpsi, erfa_deps = erfa.nut06a(d1, d2)
    moira_dpsi_deg, moira_deps_deg = nutation(jd_tt)

    dpsi_error_arcsec = abs(math.degrees(erfa_dpsi) - moira_dpsi_deg) * ARCSEC
    deps_error_arcsec = abs(math.degrees(erfa_deps) - moira_deps_deg) * ARCSEC

    assert dpsi_error_arcsec < PASS_THRESHOLD_ARCSEC
    assert deps_error_arcsec < PASS_THRESHOLD_ARCSEC


@pytest.mark.integration
@pytest.mark.parametrize(("label", "jd_tt", "_jd_ut"), TEST_EPOCHS, ids=[label for label, _, _ in TEST_EPOCHS])
def test_precession_matrix_matches_erfa(label: str, jd_tt: float, _jd_ut: float) -> None:
    d1, d2 = _split(jd_tt)
    erfa_matrix = erfa.pmat06(d1, d2)
    moira_matrix = precession_matrix(jd_tt)

    assert _matrix_max_diff_arcsec(erfa_matrix, moira_matrix) < PASS_THRESHOLD_ARCSEC


@pytest.mark.integration
@pytest.mark.parametrize(("label", "jd_tt", "_jd_ut"), TEST_EPOCHS, ids=[label for label, _, _ in TEST_EPOCHS])
def test_precession_nutation_matrix_matches_erfa(label: str, jd_tt: float, _jd_ut: float) -> None:
    d1, d2 = _split(jd_tt)
    erfa_matrix = erfa.pnm06a(d1, d2)
    moira_matrix = mat_mul(nutation_matrix_equatorial(jd_tt), precession_matrix_equatorial(jd_tt))

    assert _matrix_max_diff_arcsec(erfa_matrix, moira_matrix) < PASS_THRESHOLD_ARCSEC


@pytest.mark.integration
@pytest.mark.parametrize(("label", "jd_tt", "_jd_ut"), TEST_EPOCHS, ids=[label for label, _, _ in TEST_EPOCHS])
def test_true_obliquity_matches_erfa(label: str, jd_tt: float, _jd_ut: float) -> None:
    d1, d2 = _split(jd_tt)
    _dpsi, erfa_deps = erfa.nut06a(d1, d2)
    erfa_true = math.degrees(erfa.obl06(d1, d2) + erfa_deps)
    moira_true = true_obliquity(jd_tt)

    assert abs(erfa_true - moira_true) * ARCSEC < PASS_THRESHOLD_ARCSEC


@pytest.mark.integration
@pytest.mark.parametrize(("label", "jd_tt", "jd_ut"), TEST_EPOCHS, ids=[label for label, _, _ in TEST_EPOCHS])
def test_gast_approximation_matches_erfa(label: str, jd_tt: float, jd_ut: float) -> None:
    d1u, d2u = _split(jd_ut)
    d1t, d2t = _split(jd_tt)
    erfa_dpsi, erfa_deps = erfa.nut06a(d1t, d2t)
    erfa_eps = erfa.obl06(d1t, d2t)
    erfa_eqeq_rad = erfa_dpsi * math.cos(erfa_eps + erfa_deps)
    erfa_gast_deg = math.degrees(erfa.gmst06(d1u, d2u, d1t, d2t) + erfa_eqeq_rad) % 360.0

    moira_gmst = greenwich_mean_sidereal_time(jd_ut)
    moira_dpsi_deg, moira_deps_deg = nutation(jd_tt)
    moira_eps_rad = math.radians(mean_obliquity_p03(jd_tt) + moira_deps_deg)
    moira_gast_deg = (moira_gmst + moira_dpsi_deg * math.cos(moira_eps_rad)) % 360.0

    assert _wrapped_arcsec(erfa_gast_deg, moira_gast_deg) < PASS_THRESHOLD_ARCSEC


# Modern epochs only (J1500–J2100): full GAST via erfa.gst06a.
# Ancient epochs (pre-J1000) diverge up to ~1.1" due to a model-basis
# difference between Moira's equation-of-equinoxes approach and ERFA's
# equation-of-origins approach; see VALIDATION_ASTRONOMY.md §3.7.1.
_MODERN_EPOCHS = [
    ("J1500.0",  2268923.5),
    ("J1800.0",  2378496.5),
    ("J1900.0",  2415020.5),
    ("J2000.0",  2451545.0),
    ("J2010.0",  2455196.5),
    ("J2024.0",  2460310.5),
    ("J2050.0",  2469807.5),
    ("J2100.0",  2488069.5),
]


@pytest.mark.integration
@pytest.mark.parametrize(("label", "jd_ut"), _MODERN_EPOCHS, ids=[l for l, _ in _MODERN_EPOCHS])
def test_full_gast_matches_erfa_gst06a(label: str, jd_ut: float) -> None:
    jd_tt = ut_to_tt(jd_ut)
    d1u, d2u = _split(jd_ut)
    d1t, d2t = _split(jd_tt)
    erfa_gast = math.degrees(erfa.gst06a(d1u, d2u, d1t, d2t)) % 360.0
    moira_gast = apparent_sidereal_time_at(jd_ut) % 360.0

    assert _wrapped_arcsec(erfa_gast, moira_gast) < PASS_THRESHOLD_ARCSEC