brahe 1.6.2

Brahe is a modern satellite dynamics library for research and engineering applications designed to be easy-to-learn, high-performance, and quick-to-deploy. The north-star of the development is enabling users to solve meaningful problems and answer questions quickly, easily, and correctly.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

"""GMAT benchmark implementations for state frame conversions.

Uses CoordinateConverter with per-iteration CoordinateSystem construction
(EarthMJ2000Eq <-> EarthFixed). The gmat.Clear() call inside each iteration
wipes the named-object registry, so CoordinateSystem objects must be rebuilt
each time (matches the Basilisk per-iteration SimBaseClass pattern).

Epoch conversion: the task provides Julian Dates (UTC). GMAT's TimeSystemConverter
uses a modified Julian Date referenced to JD 2430000.0 (not the standard MJD
reference of JD 2400000.5). The correct mapping is:
    GMAT UTCMJD = JD - 2430000.0
    A1MJD = tcv.Convert(gmat_utcmjd, tcv.UTCMJD, tcv.A1MJD)

CoordinateConverter.Convert signature (5-arg, output is filled in-place):
    cc.Convert(a1mjd: float, in_state: Rvector6, in_cs: CoordinateSystem,
               out_state: Rvector6, out_cs: CoordinateSystem)

Unit conventions: task states are in meters / m/s; GMAT uses km / km/s.

Accuracy note: GMAT's EarthFixed uses FK5 + IAU 1980 nutation whereas
brahe/Orekit use IERS 2010 ITRF. Expect bounded errors ≤ ~10 m at LEO.
"""

from benchmarks.comparative.implementations.gmat.base import (
    build_task_result,
    gmat_clear,
    km_to_m_state,
    m_to_km_state,
    time_iterations,
)

# GMAT UTCMJD uses JD - 2430000.0 as its reference (not the standard JD - 2400000.5).
_GMAT_JD_REF = 2430000.0


def _jd_to_a1mjd(jd: float) -> float:
    """JD (UTC) -> GMAT A1MJD using TimeSystemConverter.

    GMAT's UTCMJD reference is JD 2430000.0, not the standard MJD reference
    of JD 2400000.5. Passing (jd - 2400000.5) to UTCMJD yields epoch errors
    of ~30,000 days, producing completely wrong rotation angles.
    """
    import gmatpy as gmat
    tcv = gmat.TimeSystemConverter.Instance()
    gmat_utcmjd = jd - _GMAT_JD_REF
    return float(tcv.Convert(gmat_utcmjd, tcv.UTCMJD, tcv.A1MJD))


def _build_systems():
    """Construct and initialize ECI and ECEF CoordinateSystems for one iteration.

    Must be called after gmat_clear() each iteration because Clear() wipes the
    named-object registry. SetSolarSystem + Initialize are required for the
    CoordinateConverter to resolve the body-fixed rotation.
    """
    import gmatpy as gmat
    ss = gmat.GetSolarSystem()
    cs_eci = gmat.Construct("CoordinateSystem", "EciCs", "Earth", "MJ2000Eq")
    cs_eci.SetSolarSystem(ss)
    cs_eci.Initialize()
    cs_ecef = gmat.Construct("CoordinateSystem", "EcefCs", "Earth", "BodyFixed")
    cs_ecef.SetSolarSystem(ss)
    cs_ecef.Initialize()
    return cs_eci, cs_ecef


def _rvec6_to_list(rv6) -> list[float]:
    """Extract a Python list[float] from a GMAT Rvector6."""
    return [float(rv6.GetElement(i)) for i in range(6)]


def state_eci_to_ecef(params: dict, iterations: int):
    """Convert ECI states (EarthMJ2000Eq) to ECEF (EarthFixed).

    Input and output states are in meters / m/s. Internal GMAT calls use km / km/s.
    Task params format: params["cases"] is a list of {"jd": float, "state": [6 floats]}.
    """
    import gmatpy as gmat

    cases = params["cases"]
    # Pre-compute A1MJD epochs and convert state units OUTSIDE the timed region.
    prepared = [
        (_jd_to_a1mjd(c["jd"]), m_to_km_state(c["state"]))
        for c in cases
    ]

    def run():
        out = []
        for a1mjd, state_km in prepared:
            gmat_clear()
            cs_eci, cs_ecef = _build_systems()
            rv_in = gmat.Rvector6(*state_km)
            rv_out = gmat.Rvector6()
            gmat.CoordinateConverter().Convert(a1mjd, rv_in, cs_eci, rv_out, cs_ecef)
            out.append(km_to_m_state(_rvec6_to_list(rv_out)))
        return out

    times, results = time_iterations(run, iterations)
    return build_task_result(
        "frames.state_eci_to_ecef",
        iterations,
        times,
        results,
        extra_metadata={"earth_fixed_frame": "EarthFixed (FK5/IAU 1980)"},
    )


def state_ecef_to_eci(params: dict, iterations: int):
    """Convert ECEF states (EarthFixed) to ECI (EarthMJ2000Eq).

    Input and output states are in meters / m/s. Internal GMAT calls use km / km/s.
    Task params format: params["cases"] is a list of {"jd": float, "state": [6 floats]}.
    """
    import gmatpy as gmat

    cases = params["cases"]
    prepared = [
        (_jd_to_a1mjd(c["jd"]), m_to_km_state(c["state"]))
        for c in cases
    ]

    def run():
        out = []
        for a1mjd, state_km in prepared:
            gmat_clear()
            cs_eci, cs_ecef = _build_systems()
            rv_in = gmat.Rvector6(*state_km)
            rv_out = gmat.Rvector6()
            # Swap source/dest to convert from ECEF to ECI.
            gmat.CoordinateConverter().Convert(a1mjd, rv_in, cs_ecef, rv_out, cs_eci)
            out.append(km_to_m_state(_rvec6_to_list(rv_out)))
        return out

    times, results = time_iterations(run, iterations)
    return build_task_result(
        "frames.state_ecef_to_eci",
        iterations,
        times,
        results,
        extra_metadata={"earth_fixed_frame": "EarthFixed (FK5/IAU 1980)"},
    )