AbstractPropType (singleton base)
├── CowellPropagator         # Cartesian state
├── GaussVEPropagator        # Gauss Variational Equations
├── EDromoPropagator          # Regularized EDromo elements
├── KustaanheimoStiefelPropagator  # KS transform
├── MilankovichPropagator     # Milankovich elements
├── StiefelScheifelePropagator     # Stiefel-Scheifele
├── USM7Propagator            # Unified State Model (7-element)
├── USM6Propagator            # Unified State Model (6-element)
└── USMEMPropagator           # Unified State Model (exponential map)

src/
  AstroPropagators.jl       # Module entry point
  api.jl                    # propagate() public API
  auxiliary/util.jl          # skew_sym, quaternions2DCM, RTN_frame
  propagators/               # One file per propagator (EOM + EOM!)
  events/impulsive_maneuvers.jl  # Impulsive burn callbacks

# Out-of-place (returns SVector, zero-allocation)
function Cowell_EOM(u, p, t, model_list)
    # ... compute accelerations ...
    return SVector{6}(du1, du2, du3, du4, du5, du6)
end

# In-place (calls out-of-place, assigns with .=)
function Cowell_EOM!(du, u, p, t, model_list)
    du .= Cowell_EOM(u, p, t, model_list)
    return nothing
end

sol = propagate(
    u0, p, models, tspan;
    prop_type=CowellPropagator(),
    config=nothing,                    # RegularizedCoordinateConfig for EDromo/KS/SS
    ODE_solver=VCABM(),
    abstol=1E-13, reltol=1E-13,
)

EOM!(du, u, p, t) = Cowell_EOM!(du, u, p, t, model_list)
prob = ODEProblem(EOM!, u0, tspan, p)
sol = solve(prob, VCABM(); abstol=1e-13, reltol=1e-13)

grav = GravityHarmonicsAstroModel(; gravity_model=..., eop_data=..., order=36, degree=36)
sun = ThirdBodyModel(; body=SunBody(), eop_data=eop_data)
moon = ThirdBodyModel(; body=MoonBody(), eop_data=eop_data)
model_list = CentralBodyDynamicsModel(grav, (sun, moon, srp, drag))

end_EDromo_integration(stop_time, config)       # ContinuousCallback for termination
EDromo_burn(burn_time, ΔV, config)              # ContinuousCallback for impulsive burn