Two Body Analytical Propagator

The two-body analytical orbit propagator considers the Earth a perfect sphere with uniform density. Hence, it propagates the orbit using the solution considering Newtonian gravity. It has an extremely low precision but with a minimal computational burden. Thus, it is helpful in some analysis that requires propagating the orbit many times for short periods.

Algorithm

The algorithm implemented here is based on [1].

Since we are considering a spherical Earth with uniform density, gravity points towards the center of Earth. Thus, we propagate the orbit by updating the satellite mean anomaly since all other Keplerian elements do not change. The equation to correct the mean anomaly is:

\[M(t) = M_0 + \sqrt{\frac{\mu}{a_0^3}} \cdot \left(t - t_0\right)\]

where $t_0$ is the initial mean elements' epoch, $a_0$ is the mean semi-major axis, and $\mu$ the Earth's standard gravitational parameter.

Initialization

We can initialize the two-body analytical propagator with the following function:

Propagators.init(Val(:TwoBody), orb₀::KeplerianElements; kwargs...) -> OrbitPropagatorTwoBody

which creates a two-body propagator structure OrbitPropagatorTwoBody with the mean Keplerian elements orb₀. The following keyword selects the standard gravitational parameter for the propagation algorithm:

• m0::T: Standard gravitational parameter of the central body [m³/s²]. (Default: TBC_M0)

This package contains some pre-built gravitational parameters of the Earth for this propagator:

julia> orb = KeplerianElements(
                  date_to_jd(2023, 1, 1, 0, 0, 0),
                  7190.982e3,
                  0.001111,
                  98.405 |> deg2rad,
                  100    |> deg2rad,
                  90     |> deg2rad,
                  19     |> deg2rad
              )KeplerianElements{TrueAnomaly, Float64, Float64}:
  Epoch             : 2.45995e6 (2023-01-01T00:00:00)
  Semi-Major Axis   : 7190.982 km
  Eccentricity      : 0.001111
  Inclination       : 98.405°
  RA of Asc. Node   : 100.0°
  Arg. of Periapsis : 90.0°
  True Anomaly      : 19.0°
julia> orbp = Propagators.init(Val(:TwoBody), orb)OrbitPropagatorTwoBody{Float64, Float64} (Two-Body Orbit Propagator):
  ├─ Mean Elements
  │    Epoch             : 2.45995e6 (2023-01-01T00:00:00)
  │    Semi-Major Axis   : 7190.982 km
  │    Eccentricity      : 0.001111
  │    Inclination       : 98.405°
  │    RA of Asc. Node   : 100.0°
  │    Arg. of Periapsis : 90.0°
  │    Mean Anomaly      : 18.95858415°
  ├─ Secular Rates
  │    Mean Motion : 14.23706472 rev/day
  ├─ Constants
  │    μ : 3.986004415e14 m³/s²
  └─ Propagation
       Last Instant : 0.0 s

Fitting Mean Elements

We can use the function:

Propagators.fit_mean_elements(::Val{:TwoBody}, vjd::AbstractVector{Tjd}, vr_i::AbstractVector{Tv}, vv_i::AbstractVector{Tv}; kwargs...) -> KeplerianElements{MeanAnomaly, Float64, T}, SMatrix{6, 6, T}, NamedTuple

to fit a set of mean Keplerian elements for the two-body orbit propagator using the osculating elements represented by a set of position vectors vr_i [m] and a set of velocity vectors vv_i [m / s] represented in an inertial reference frame at instants in the array vjd [Julian Day].

It returns the fitted Keplerian elements, the final covariance matrix of the least-square algorithm, and a NamedTuple with its statistics: converged, iterations, position_rmse [m], velocity_rmse [m / s], and total_rmse.

The following keywords are available to configure the fitting process:

• m0::T: Standard gravitational parameter of the central body [m³ / s²], whose number type T is used in the fitting. (Default: TBC_M0)
• atol::Number: Tolerance for the residual absolute value. If the residual is lower than atol at any iteration, the computation loop stops. (Default: 2e-4)
• rtol::Number: Tolerance for the relative difference between the residuals. If the relative difference between the residuals in two consecutive iterations is lower than rtol, the computation loop stops. (Default: 2e-4)
• initial_guess::Union{Nothing, KeplerianElements}: Initial guess for the mean elements fitting process. If it is nothing, the algorithm will obtain an initial estimate from the osculating elements in vr_i and vv_i. (Default: nothing)
• jacobian_method::AbstractJacobianMethod: Method used to compute the Jacobian matrix. It can be FiniteDiffJacobian() for finite differences or ForwardDiffJacobian() for ForwardDiff.jl automatic differentiation. (Default: FiniteDiffJacobian())
• jacobian_perturbation::Number: Initial state perturbation to compute the finite-difference when calculating the Jacobian matrix. Only used with FiniteDiffJacobian(). (Default: 1e-3)
• jacobian_perturbation_tol::Number: Tolerance to accept the perturbation when calculating the Jacobian matrix. If the computed perturbation is lower than jacobian_perturbation_tol, we increase it until its absolute value is higher than jacobian_perturbation_tol. Only used with FiniteDiffJacobian(). (Default: 1e-7)
• max_iterations::Int: Maximum number of iterations allowed for the least-square fitting. (Default: 50)
• mean_elements_epoch::Union{Number, DateTime}: Epoch of the fitted mean elements, represented by a Julian Day [UTC] or a DateTime [UTC]. (Default: vjd[end])
• verbose::Bool: If true, the algorithm prints debugging information to stdout. (Default: true)
• weight_vector::AbstractVector: Vector with the measurements weights for the least-square algorithm. We assemble the weight matrix W as a diagonal matrix with the elements in weight_vector at its diagonal. (Default: @SVector(ones(Bool, 6)))

julia> vr_i = [
           [-6792.402703741442, 2192.6458461287293, 0.18851758695295118] .* 1000,
           [-6357.88873265975, 2391.9476768911686, 2181.838771262736] .* 1000
       ];
julia> vv_i = [
           [0.3445760107690598, 1.0395135806993514, 7.393686131436984] .* 1000,
           [2.5285015912807003, 0.27812476784300005, 7.030323100703928] .* 1000
       ];
julia> vjd = [
           2.46002818657856e6,
           2.460028190050782e6
       ];
julia> orb, P, stats = Propagators.fit_mean_elements(Val(:TwoBody), vjd, vr_i, vv_i)ACTION:   Fitting the mean elements for the two-body propagator.
           Iteration        Position RMSE        Velocity RMSE           Total RMSE       RMSE Variation
                                     [km]             [km / s]                  [ ]


PROGRESS:          1             0.306117           0.00205737              306.124                  ---

PROGRESS:          2          2.60894e-06           0.00142531              1.42531             -99.5344 %

PROGRESS:          3          1.67398e-06           0.00142531              1.42531         -4.87677e-05 %

(KeplerianElements{MeanAnomaly, Float64, Float64}: Epoch = 2.46003e6 (2023-03-24T16:33:40.388), [0.9999772153874656 3.4362845724651176e-7 … -9.849257846270436e-5 -6.85147872850814e-5; 3.4362850838986486e-7 0.9999780011263167 … 0.003258577044732411 2.5060982376383854e-5; … ; -9.84925784624993e-5 0.0032585770447324676 … 2.1999356161871043e-5 1.1527516890200296e-7; -6.851478728036799e-5 2.5060982371481454e-5 … 1.1527516888462523e-7 2.196482928505477e-5], (converged = true, iterations = 3, position_rmse = 0.0016739781438326025, velocity_rmse = 1.4253086768039351, total_rmse = 1.4253096598197217))
julia> orbKeplerianElements{MeanAnomaly, Float64, Float64}:
  Epoch             : 2.46003e6 (2023-03-24T16:33:40.388)
  Semi-Major Axis   : 7139.282635 km
  Eccentricity      : 0.001327104531
  Inclination       : 98.42538016°
  RA of Asc. Node   : 162.1096997°
  Arg. of Periapsis : 79.45743611°
  Mean Anomaly      : 298.683562°
julia> stats(converged = true, iterations = 3, position_rmse = 0.0016739781438326025, velocity_rmse = 1.4253086768039351, total_rmse = 1.4253096598197217)

References

• [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications. 4th ed. Microcosm Press, Hawthorn, CA, USA.