J4 Osculating Analytical Orbit Propagator

This algorithm uses the J4 propagator to obtain the secular effects of the Keplerian elements caused by the terms $J_2$, $J_2^2$, and $J_4$of the geopotential field. Afterward, it adds short-term perturbations using only the $J_2$ term. This model is useful when fitting an orbit to a set of mean elements for the J4 orbit propagator.

Algorithm

The algorithm implemented here is based on [1].

Initialization

We can initialize the J4 osculating analytical orbit propagator with the following function:

Propagators.init(Val(:J4osc), orb₀::KeplerianElements; kwargs...) -> OrbitPropagatorJ4Osculating

which creates a J4 osculating propagator structure OrbitPropagatorJ4Osculating with the mean Keplerian elements orb₀.

The following keyword selects the gravitational constants for the propagation algorithm:

• j4c::J4PropagatorConstants: J4 orbit propagator constants (see J4PropagatorConstants). (Default: J4C_EGM2008)

This package contains some pre-built propagation constants 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(:J4osc), orb)OrbitPropagatorJ4Osculating{Float64, Float64} (J4 Osculating 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.22855879 rev/day
  │    RAAN Rate              : 0.955513678 °/day
  │    Arg. of Periapsis Rate : -2.921463798 °/day
  ├─ Constants
  │    R₀ : 6378.137 km
  │    μm : 0.001239447462 rad/s
  │    J₂ : 0.001082626174
  │    J₄ : -1.6198976e-6
  └─ Propagation
       Last Instant : 0.0 s

Fitting Mean Elements

We can use the function:

Propagators.fit_mean_elements(::Val{:J4osc}, 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 J4 osculating 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:

• j4c::J4PropagatorConstants{T}: J4 orbit propagator constants (see J4PropagatorConstants), whose number type T is used in the fitting. (Default: J4C_EGM2008)
• 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(:J4osc), vjd, vr_i, vv_i)ACTION:   Fitting the mean elements for the J4 osculating propagator.
           Iteration        Position RMSE        Velocity RMSE           Total RMSE       RMSE Variation
                                     [km]             [km / s]                  [ ]


PROGRESS:          1              6.24852           0.00466097              6248.52                  ---

PROGRESS:          2          0.000816956            0.0025964              2.72189             -99.9564 %

PROGRESS:          3          1.68789e-05           0.00259611              2.59616              -4.6192 %

PROGRESS:          4          1.68763e-05           0.00259611              2.59616          3.00631e-10 %

(KeplerianElements{MeanAnomaly, Float64, Float64}: Epoch = 2.46003e6 (2023-03-24T16:33:40.388), [0.9999427680963674 -0.004901179136806686 … -0.00012021347596864274 -0.00011550159599902191; -0.00490117915398981 1.0007325998974053 … 0.0032661543929103955 3.8582051397170914e-5; … ; -0.00012021347602330139 0.003266154392909064 … 2.20482279968684e-5 5.398238168448781e-8; -0.00011550159599370338 3.858205139704896e-5 … 5.3982381681982574e-8 2.1657607939733534e-5], (converged = true, iterations = 4, position_rmse = 0.016876287593669836, velocity_rmse = 2.5961099245046393, total_rmse = 2.5961647769728384))
julia> orbKeplerianElements{MeanAnomaly, Float64, Float64}:
  Epoch             : 2.46003e6 (2023-03-24T16:33:40.388)
  Semi-Major Axis   : 7135.792461 km
  Eccentricity      : 0.001353135365
  Inclination       : 98.4304116°
  RA of Asc. Node   : 162.1131631°
  Arg. of Periapsis : 64.96868276°
  Mean Anomaly      : 313.1552992°
julia> stats(converged = true, iterations = 4, position_rmse = 0.016876287593669836, velocity_rmse = 2.5961099245046393, total_rmse = 2.5961647769728384)

References

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