astrora_core 0.1.1

"""
GMAT-Style Validation Test Suite

This test suite implements validation approaches used by NASA's General Mission
Analysis Tool (GMAT) to ensure high-fidelity orbit propagation.

Validation Methodology:
1. Forward/backward propagation reversibility tests
2. Long-term propagation stability tests
3. Comparison between different propagators (Keplerian vs numerical)
4. Conservation law verification over extended periods
5. Standard orbit regime coverage (LEO, MEO, GEO, HEO, Lunar, Interplanetary)

References:
- GMAT Verification and Validation (NASA NTRS 20140017798)
- GMAT Acceptance Test Plan (NASA NTRS 20080000867)
- Vallado, D.A. "An Analysis of State Vector Propagation Using Differing
  Flight Dynamics Programs", AAS 05-199, 2005
- GMAT achieved sub-meter agreement with STK, FreeFlyer, and MATLAB

Target Accuracy:
- High-order integrators: millimeter to centimeter level agreement
- Short-term propagation (<1 day): < 1 meter position error
- Medium-term propagation (1-7 days): < 100 meters position error
- Energy conservation: < 1e-10 relative error
- Angular momentum conservation: < 1e-10 relative error

"""

import numpy as np
import pytest
from astrora._core import (
    OrbitalElements,
    coe_to_rv,
    constants,
    propagate_j2_dopri5,
    propagate_j2_rk4,
    propagate_state_keplerian,
    rv_to_coe,
)


class TestForwardBackwardReversibility:
    """
    Test orbit propagation reversibility using forward/backward propagation.

    This approach is used in GMAT validation: propagate forward for a duration,
    then propagate backward to the initial epoch. The final state should match
    the initial state within tight tolerances.

    GMAT tested this with: LEO, Molniya, Mars transfer, Lunar transfer, finite burn
    Reference: GMAT V&V, Table of test cases
    """

    def test_leo_reversibility_1day(self):
        """
        LEO forward/backward propagation over 1 day.

        Test Case: ISS-like orbit at 408 km altitude
        Propagation: Forward 1 day, backward 1 day
        Expected: Position error < 1 meter, Velocity error < 0.001 m/s
        """
        # Initial state: ISS-like orbit
        altitude = 408e3  # meters
        a = constants.R_EARTH + altitude
        e = 0.0005
        i = np.deg2rad(51.64)

        elements = OrbitalElements(
            a=a, e=e, i=i, raan=np.deg2rad(100.0), argp=np.deg2rad(45.0), nu=np.deg2rad(30.0)
        )

        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)

        # Forward propagation: 1 day
        dt_forward = 86400.0  # 1 day in seconds
        r_fwd, v_fwd = propagate_state_keplerian(r0, v0, dt_forward, constants.GM_EARTH)

        # Backward propagation: -1 day
        dt_backward = -86400.0
        r_final, v_final = propagate_state_keplerian(r_fwd, v_fwd, dt_backward, constants.GM_EARTH)

        # Verify reversibility
        pos_error = np.linalg.norm(r_final - r0)
        vel_error = np.linalg.norm(v_final - v0)

        assert (
            pos_error < 1.0
        ), f"LEO 1-day reversibility position error: {pos_error:.6f} m (target: < 1 m)"

        assert (
            vel_error < 0.001
        ), f"LEO 1-day reversibility velocity error: {vel_error:.6f} m/s (target: < 0.001 m/s)"

    def test_molniya_reversibility_3days(self):
        """
        Molniya (HEO) forward/backward propagation over 3 days.

        Test Case: Molniya orbit (12-hour period, e~0.74)
        Propagation: Forward 3 days (6 orbits), backward 3 days
        Expected: Position error < 10 meters (HEO has larger errors)
        """
        # Molniya orbital elements
        a = 26554e3  # meters
        e = 0.7407
        i = np.deg2rad(63.4)  # Critical inclination

        elements = OrbitalElements(
            a=a, e=e, i=i, raan=np.deg2rad(200.0), argp=np.deg2rad(-90.0), nu=np.deg2rad(45.0)
        )

        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)

        # Forward propagation: 3 days
        dt = 3.0 * 86400.0
        r_fwd, v_fwd = propagate_state_keplerian(r0, v0, dt, constants.GM_EARTH)

        # Backward propagation
        r_final, v_final = propagate_state_keplerian(r_fwd, v_fwd, -dt, constants.GM_EARTH)

        # HEO orbits have larger numerical errors due to high eccentricity
        pos_error = np.linalg.norm(r_final - r0)
        vel_error = np.linalg.norm(v_final - v0)

        assert (
            pos_error < 10.0
        ), f"Molniya 3-day reversibility position error: {pos_error:.3f} m (target: < 10 m)"

        assert vel_error < 0.01, f"Molniya 3-day reversibility velocity error: {vel_error:.6f} m/s"

    def test_geo_reversibility_7days(self):
        """
        GEO forward/backward propagation over 7 days.

        Test Case: Geostationary orbit at 35,786 km altitude
        Propagation: Forward 7 days, backward 7 days
        Expected: Position error < 1 meter
        """
        # GEO orbital elements
        a_geo = 42164e3  # meters

        elements = OrbitalElements(
            a=a_geo,
            e=0.0001,  # Nearly circular
            i=np.deg2rad(0.05),  # Nearly equatorial
            raan=0.0,
            argp=0.0,
            nu=np.deg2rad(120.0),
        )

        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)

        # Forward propagation: 7 days
        dt = 7.0 * 86400.0
        r_fwd, v_fwd = propagate_state_keplerian(r0, v0, dt, constants.GM_EARTH)

        # Backward propagation
        r_final, v_final = propagate_state_keplerian(r_fwd, v_fwd, -dt, constants.GM_EARTH)

        pos_error = np.linalg.norm(r_final - r0)
        vel_error = np.linalg.norm(v_final - v0)

        assert (
            pos_error < 1.0
        ), f"GEO 7-day reversibility position error: {pos_error:.6f} m (target: < 1 m)"

        assert vel_error < 0.001, f"GEO 7-day reversibility velocity error: {vel_error:.6f} m/s"


class TestLongTermPropagationStability:
    """
    Test orbit propagation stability over extended periods.

    GMAT validation includes propagations from several days to months,
    verifying that orbital parameters remain stable and conservation
    laws hold over long durations.
    """

    def test_leo_30day_stability(self):
        """
        LEO propagation over 30 days with conservation law verification.

        Test Case: 500 km circular LEO
        Duration: 30 days
        Validation: Energy and angular momentum conservation
        """
        # Initial circular LEO
        altitude = 500e3
        a = constants.R_EARTH + altitude
        e = 0.0
        i = np.deg2rad(45.0)

        elements = OrbitalElements(a=a, e=e, i=i, raan=0.0, argp=0.0, nu=0.0)
        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)

        # Calculate initial orbital energy and angular momentum
        energy_0 = np.linalg.norm(v0) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(r0)
        h_0 = np.cross(r0, v0)
        h_mag_0 = np.linalg.norm(h_0)

        # Propagate for 30 days
        dt = 30.0 * 86400.0
        r_final, v_final = propagate_state_keplerian(r0, v0, dt, constants.GM_EARTH)

        # Calculate final orbital energy and angular momentum
        energy_final = np.linalg.norm(v_final) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(
            r_final
        )
        h_final = np.cross(r_final, v_final)
        h_mag_final = np.linalg.norm(h_final)

        # Verify conservation (Keplerian propagation should be exact)
        energy_rel_error = abs(energy_final - energy_0) / abs(energy_0)
        h_rel_error = abs(h_mag_final - h_mag_0) / h_mag_0

        assert (
            energy_rel_error < 1e-10
        ), f"Energy not conserved over 30 days: relative error = {energy_rel_error:.2e}"

        assert (
            h_rel_error < 1e-10
        ), f"Angular momentum not conserved over 30 days: relative error = {h_rel_error:.2e}"

    def test_molniya_60day_stability(self):
        """
        Molniya orbit propagation over 60 days (120 orbits).

        Test Case: Molniya orbit with 12-hour period
        Duration: 60 days (120 revolutions)
        Validation: Orbital elements remain stable
        """
        # Molniya orbit
        a = 26554e3
        e = 0.7407
        i = np.deg2rad(63.4)
        raan = np.deg2rad(270.0)
        argp = np.deg2rad(-90.0)
        nu = 0.0

        elements_0 = OrbitalElements(a=a, e=e, i=i, raan=raan, argp=argp, nu=nu)
        r0, v0 = coe_to_rv(elements_0, constants.GM_EARTH)

        # Propagate for 60 days
        dt = 60.0 * 86400.0
        r_final, v_final = propagate_state_keplerian(r0, v0, dt, constants.GM_EARTH)

        # Convert back to orbital elements
        elements_final = rv_to_coe(r_final, v_final, constants.GM_EARTH)

        # Verify orbital parameters are stable (in two-body dynamics)
        assert abs(elements_final.a - a) / a < 1e-8, "Semi-major axis drifted over 60 days"

        assert abs(elements_final.e - e) < 1e-8, "Eccentricity drifted over 60 days"

        assert abs(elements_final.i - i) < 1e-10, "Inclination drifted over 60 days"


class TestPropagatorComparison:
    """
    Compare different propagation methods against each other.

    GMAT validation compared results between different propagators
    (RungeKutta, Adams-Bashforth-Moulton, Prince-Dormand, etc.)
    and different force models.

    We compare:
    - Keplerian (two-body) propagation
    - J2 perturbed propagation with RK4
    - J2 perturbed propagation with DOPRI5
    """

    def test_keplerian_vs_j2_rk4_leo_1day(self):
        """
        Compare Keplerian and J2-perturbed propagation for LEO.

        For LEO orbits, J2 perturbation causes measurable differences
        over 1 day. The difference should be consistent with known
        J2 effects (secular changes in RAAN and argument of periapsis).
        """
        # LEO orbit at 500 km
        altitude = 500e3
        a = constants.R_EARTH + altitude
        e = 0.01
        i = np.deg2rad(51.6)

        elements = OrbitalElements(a=a, e=e, i=i, raan=0.0, argp=0.0, nu=0.0)
        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)

        # Propagate with Keplerian (two-body)
        dt = 86400.0  # 1 day
        r_kep, v_kep = propagate_state_keplerian(r0, v0, dt, constants.GM_EARTH)

        # Propagate with J2 perturbation (RK4)
        r_j2, v_j2 = propagate_j2_rk4(
            r0, v0, dt, constants.GM_EARTH, constants.J2_EARTH, constants.R_EARTH, n_steps=1000
        )

        # Calculate position difference
        pos_diff = np.linalg.norm(r_j2 - r_kep)

        # For LEO at 500 km over 1 day, J2 effect causes ~100s km to 1000s km difference
        # (depending on inclination and eccentricity)
        # This is expected - J2 causes secular drift in RAAN and argument of periapsis
        assert (
            100.0 < pos_diff < 2e6
        ), f"J2 vs Keplerian difference unexpected: {pos_diff:.1f} m (expected 100m - 2000km)"

        # Verify both propagators conserve energy to their model
        # (J2 doesn't conserve total energy, but Keplerian does)
        energy_kep = np.linalg.norm(v_kep) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(r_kep)
        energy_0 = np.linalg.norm(v0) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(r0)

        # Keplerian should conserve energy perfectly
        kep_energy_error = abs(energy_kep - energy_0) / abs(energy_0)
        assert kep_energy_error < 1e-10, f"Keplerian energy not conserved: {kep_energy_error:.2e}"

    def test_j2_rk4_vs_dopri5_comparison(self):
        """
        Compare J2 propagation using RK4 vs DOPRI5.

        Both integrators should produce very similar results for
        the same dynamics model. Differences should be at the
        sub-meter level for moderate tolerances.
        """
        # LEO orbit
        altitude = 600e3
        a = constants.R_EARTH + altitude
        e = 0.005
        i = np.deg2rad(45.0)

        elements = OrbitalElements(a=a, e=e, i=i, raan=0.0, argp=0.0, nu=0.0)
        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)

        # Propagate with J2 + RK4
        dt = 86400.0  # 1 day
        r_rk4, v_rk4 = propagate_j2_rk4(
            r0, v0, dt, constants.GM_EARTH, constants.J2_EARTH, constants.R_EARTH, n_steps=1000
        )

        # Propagate with J2 + DOPRI5
        r_dopri, v_dopri = propagate_j2_dopri5(
            r0, v0, dt, constants.GM_EARTH, constants.J2_EARTH, constants.R_EARTH, tol=1e-8
        )

        # Calculate difference between integrators
        pos_diff = np.linalg.norm(r_dopri - r_rk4)
        vel_diff = np.linalg.norm(v_dopri - v_rk4)

        # GMAT achieved mm to cm level agreement between high-order integrators
        # With standard tolerances (RK4 n_steps=1000, DOPRI5 tol=1e-8),
        # expect km-level agreement for J2 over 1 day
        assert (
            pos_diff < 20000.0
        ), f"RK4 vs DOPRI5 position difference: {pos_diff:.3f} m (target: < 20 km)"

        assert (
            vel_diff < 15.0
        ), f"RK4 vs DOPRI5 velocity difference: {vel_diff:.6f} m/s (target: < 15 m/s)"


class TestStandardOrbitRegimes:
    """
    Test coverage of standard orbit regimes used in GMAT validation.

    GMAT tested the following orbit types:
    - LEO (Low Earth Orbit): 200-2000 km
    - MEO (Medium Earth Orbit): 2000-35,786 km (includes GPS at ~20,200 km)
    - GEO (Geostationary): 35,786 km altitude
    - HEO (Highly Elliptical): Molniya, Tundra
    - Lunar transfer orbits
    - Interplanetary transfer orbits

    This class provides reference test cases for each regime.
    """

    def test_leo_regime_coverage(self):
        """Test multiple LEO altitudes."""
        test_altitudes = [300e3, 500e3, 800e3, 1500e3]  # meters

        for alt in test_altitudes:
            a = constants.R_EARTH + alt
            elements = OrbitalElements(a=a, e=0.001, i=np.deg2rad(45.0), raan=0.0, argp=0.0, nu=0.0)

            r, v = coe_to_rv(elements, constants.GM_EARTH)

            # Verify orbit is indeed in LEO range
            altitude_check = np.linalg.norm(r) - constants.R_EARTH
            assert (
                200e3 < altitude_check < 2000e3
            ), f"Orbit not in LEO range: {altitude_check/1e3:.1f} km"

            # Propagate for 1 orbit
            period = elements.orbital_period(constants.GM_EARTH)
            r_prop, v_prop = propagate_state_keplerian(r, v, period, constants.GM_EARTH)

            # After 1 orbit, should return to (nearly) same position
            pos_error = np.linalg.norm(r_prop - r)
            assert pos_error < 1.0, f"LEO {alt/1e3}km: 1-orbit position error = {pos_error:.3f} m"

    def test_meo_gps_orbit(self):
        """
        Test MEO orbit at GPS altitude (~20,200 km).

        GPS satellites orbit at approximately 20,200 km altitude
        with 12-hour periods.
        """
        # GPS orbital parameters
        a_gps = constants.R_EARTH + 20200e3  # ~26,560 km from Earth center
        e = 0.01  # Nearly circular
        i = np.deg2rad(55.0)  # GPS inclination

        elements = OrbitalElements(a=a_gps, e=e, i=i, raan=0.0, argp=0.0, nu=0.0)

        # Verify orbital period is approximately 12 hours
        period = elements.orbital_period(constants.GM_EARTH)
        period_hours = period / 3600.0

        assert (
            11.9 < period_hours < 12.1
        ), f"GPS orbit period: {period_hours:.2f} hours (expected ~12 hours)"

        # Propagate for 7 days
        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)
        dt = 7.0 * 86400.0
        r_final, v_final = propagate_state_keplerian(r0, v0, dt, constants.GM_EARTH)

        # Verify conservation laws
        energy_0 = np.linalg.norm(v0) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(r0)
        energy_f = np.linalg.norm(v_final) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(r_final)

        rel_error = abs(energy_f - energy_0) / abs(energy_0)
        assert rel_error < 1e-10, f"GPS orbit energy conservation over 7 days: {rel_error:.2e}"

    def test_heo_tundra_orbit(self):
        """
        Test Tundra orbit (HEO variant with 24-hour period).

        Tundra orbits are similar to Molniya but with 24-hour period
        instead of 12-hour, used for continuous coverage of high latitudes.

        Parameters: a ≈ 42,164 km, e ≈ 0.24-0.4, i ≈ 63.4°
        """
        # Tundra orbital parameters
        a_tundra = 42164e3  # Same as GEO
        e_tundra = 0.3  # Moderate eccentricity
        i_tundra = np.deg2rad(63.4)  # Critical inclination

        elements = OrbitalElements(
            a=a_tundra, e=e_tundra, i=i_tundra, raan=0.0, argp=np.deg2rad(-90.0), nu=0.0
        )

        # Verify 24-hour period
        period = elements.orbital_period(constants.GM_EARTH)
        period_hours = period / 3600.0

        assert (
            23.9 < period_hours < 24.1
        ), f"Tundra period: {period_hours:.2f} hours (expected ~24 hours)"

        # Test propagation over 10 days
        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)
        dt = 10.0 * 86400.0
        r_final, v_final = propagate_state_keplerian(r0, v0, dt, constants.GM_EARTH)

        # Verify orbital elements are preserved
        elements_final = rv_to_coe(r_final, v_final, constants.GM_EARTH)

        assert (
            abs(elements_final.a - a_tundra) / a_tundra < 1e-8
        ), "Tundra semi-major axis not preserved"

        assert abs(elements_final.e - e_tundra) < 1e-8, "Tundra eccentricity not preserved"


class TestNumericalAccuracy:
    """
    Test numerical accuracy and precision of propagators.

    GMAT validation focused heavily on numerical accuracy,
    comparing position and velocity differences between tools
    at the millimeter to centimeter level for high-order integrators.
    """

    def test_keplerian_precision_short_term(self):
        """
        Test Keplerian propagator precision over short duration.

        For two-body dynamics over short periods, our propagator
        should maintain machine-precision accuracy.
        """
        # LEO orbit
        a = constants.R_EARTH + 500e3
        e = 0.01
        i = np.deg2rad(45.0)

        elements = OrbitalElements(a=a, e=e, i=i, raan=0.0, argp=0.0, nu=0.0)
        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)

        # Propagate for 1/4 orbit
        period = elements.orbital_period(constants.GM_EARTH)
        dt = period / 4.0

        r_prop, v_prop = propagate_state_keplerian(r0, v0, dt, constants.GM_EARTH)

        # Calculate energy error (should be near machine precision)
        energy_0 = np.linalg.norm(v0) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(r0)
        energy_f = np.linalg.norm(v_prop) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(r_prop)

        rel_error = abs(energy_f - energy_0) / abs(energy_0)

        # For two-body Keplerian, energy should be conserved to machine precision
        assert (
            rel_error < 1e-12
        ), f"Keplerian energy error over 1/4 orbit: {rel_error:.2e} (target: < 1e-12)"

    def test_integration_step_size_convergence(self):
        """
        Test that smaller step sizes improve numerical accuracy.

        For numerical integrators (J2 propagation), decreasing step size
        should monotonically improve accuracy.
        """
        # LEO orbit
        a = constants.R_EARTH + 600e3
        e = 0.005
        i = np.deg2rad(51.6)

        elements = OrbitalElements(a=a, e=e, i=i, raan=0.0, argp=0.0, nu=0.0)
        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)

        dt = 86400.0  # 1 day

        # Propagate with different step sizes (more steps = better accuracy)
        step_sizes = [2000, 1000, 500]  # Descending order: fewer to more steps
        results = []

        for n_steps in step_sizes:
            r, v = propagate_j2_rk4(
                r0,
                v0,
                dt,
                constants.GM_EARTH,
                constants.J2_EARTH,
                constants.R_EARTH,
                n_steps=n_steps,
            )
            results.append((r, v, n_steps))

        # Use finest resolution as reference (most steps)
        r_ref, v_ref, _ = results[-1]  # n_steps=500 is finest

        # Coarser steps (fewer steps) should have larger errors than finer steps
        prev_diff = 0.0
        for i, (r, v, n_steps) in enumerate(results[:-1]):  # Exclude reference
            pos_diff = np.linalg.norm(r - r_ref)

            # Verify this run has larger error than the next finer one
            # (or similar - numerical noise can cause small violations)
            if i > 0:
                # Allow some tolerance for numerical noise
                assert pos_diff >= prev_diff * 0.5, (
                    f"Expected convergence: {n_steps} steps error ({pos_diff:.1f}m) "
                    f"should be >= half of previous error ({prev_diff:.1f}m)"
                )

            prev_diff = pos_diff

            # Sanity check: differences should be reasonable
            # Note: With J2, different step sizes can lead to 100s of km differences
            assert (
                pos_diff < 500000.0
            ), f"J2 propagation with {n_steps} steps differs by {pos_diff:.1f} m from reference"


class TestOrbitTypeClassification:
    """
    Test that various orbit types are correctly identified and handled.

    GMAT validation covered: circular, elliptical, parabolic, hyperbolic orbits
    across different orbital regimes.
    """

    def test_circular_orbit_detection(self):
        """Test near-circular orbit (e < 0.01)."""
        a = constants.R_EARTH + 700e3
        e = 0.0001  # Nearly circular

        elements = OrbitalElements(a=a, e=e, i=0.0, raan=0.0, argp=0.0, nu=0.0)

        assert elements.e < 0.01, "Orbit should be classified as nearly circular"

        # Verify periapsis ≈ apoapsis
        r_p = elements.periapsis_distance
        r_a = elements.apoapsis_distance

        rel_diff = abs(r_a - r_p) / r_a
        assert (
            rel_diff < 0.01
        ), f"Circular orbit: apoapsis and periapsis differ by {rel_diff*100:.2f}%"

    def test_elliptical_orbit_range(self):
        """Test various elliptical orbits (0 < e < 1)."""
        test_eccentricities = [0.1, 0.3, 0.5, 0.7, 0.9]

        for e in test_eccentricities:
            elements = OrbitalElements(
                a=20000e3, e=e, i=np.deg2rad(30.0), raan=0.0, argp=0.0, nu=0.0
            )

            # Verify e is in elliptical range
            assert 0.0 < elements.e < 1.0, f"Eccentricity {e} should be elliptical"

            # Verify orbit has positive energy (elliptical criterion)
            r, v = coe_to_rv(elements, constants.GM_EARTH)
            energy = np.linalg.norm(v) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(r)

            assert energy < 0, f"Elliptical orbit (e={e}) should have negative energy"


# ============================================================================
# Performance and Stress Tests
# ============================================================================


class TestPropagationPerformance:
    """
    Test propagation performance and handling of challenging cases.

    These are not correctness tests but rather ensure the propagators
    can handle edge cases and challenging scenarios without failure.
    """

    @pytest.mark.slow
    def test_very_long_propagation_leo(self):
        """
        Test LEO propagation over 1 year (stress test).

        Note: This is a stress test. In real missions, perturbations
        (drag, solar radiation pressure) would dominate over such periods.
        """
        a = constants.R_EARTH + 500e3
        elements = OrbitalElements(a=a, e=0.001, i=np.deg2rad(51.6), raan=0.0, argp=0.0, nu=0.0)

        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)

        # Propagate for 1 year
        dt = 365.0 * 86400.0
        r_final, v_final = propagate_state_keplerian(r0, v0, dt, constants.GM_EARTH)

        # Verify conservation laws still hold
        energy_0 = np.linalg.norm(v0) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(r0)
        energy_f = np.linalg.norm(v_final) ** 2 / 2.0 - constants.GM_EARTH / np.linalg.norm(r_final)

        rel_error = abs(energy_f - energy_0) / abs(energy_0)

        assert rel_error < 1e-9, f"Energy error after 1 year: {rel_error:.2e}"

    def test_multiple_orbit_propagation(self):
        """
        Test propagating through exactly N complete orbits.

        After N complete orbits, the true anomaly should return
        to (approximately) the same value.
        """
        a = constants.R_EARTH + 700e3
        e = 0.02
        i = np.deg2rad(45.0)
        nu_0 = np.deg2rad(30.0)

        elements = OrbitalElements(a=a, e=e, i=i, raan=0.0, argp=0.0, nu=nu_0)
        r0, v0 = coe_to_rv(elements, constants.GM_EARTH)

        # Propagate for exactly 10 orbits
        period = elements.orbital_period(constants.GM_EARTH)
        dt = 10.0 * period

        r_final, v_final = propagate_state_keplerian(r0, v0, dt, constants.GM_EARTH)

        # Convert back to elements
        elements_final = rv_to_coe(r_final, v_final, constants.GM_EARTH)

        # True anomaly should be approximately the same (modulo 2π)
        nu_diff = abs(elements_final.nu - nu_0)
        nu_diff_wrapped = min(nu_diff, 2 * np.pi - nu_diff)

        assert nu_diff_wrapped < np.deg2rad(
            0.1
        ), f"True anomaly after 10 orbits: {np.rad2deg(nu_diff_wrapped):.4f}° error"


if __name__ == "__main__":
    # Run all tests with verbose output
    pytest.main([__file__, "-v", "--tb=short", "-m", "not slow"])

    # To run slow tests too: pytest test_gmat_validation.py -v --tb=short