J4 Analytical Orbit Propagator

The J4 analytical orbit propagator considers the perturbation caused by the terms $J_2$, $J_2^2$, and $J_4$. It is a slightly more precise version of the J2 analytical orbit propagator useful for mission design and analysis.

Algorithm

The algorithm implemented here is based on [1, 2]. This version takes into account only the secular effects in the orbit elements.

Considering only the terms $J_2$, $J_2^2$, and $J_4$, we obtain that only the right accession of the ascending node (RAAN), the argument of perigee, and mean motion are perturbed by the geopotential field with secular effects. Hence, the algorithm propagates the mean elements using the following equations:

\[\begin{aligned} & p_0 = a_0 \cdot (1 - e_0^2) \\ & \beta = \sqrt{1 - e_0^2} \\ & \bar{n} = n_0 \cdot \left[ \right. \\ & \begin{aligned} \quad 1 + \\ \quad \frac{ 3}{ 4} & \cdot J_2 \cdot \left(\frac{R_0}{p_0}\right)^2 \cdot \beta \cdot (2 - 3 \sin^2 i_0) + \\ \quad \frac{ 3}{128} & \cdot J_2^2 \cdot \left(\frac{R_0}{p_0}\right)^4 \cdot \beta \cdot \left(120 + 64\beta - 40\beta^2 + (-240 - 192\beta + 40\beta^2) \sin{i_0}^2 + (105 + 144\beta + 25\beta^2) \sin{i_0}^4\right) -\\ \quad \frac{45}{128} & \cdot J_4 \cdot \left(\frac{R_0}{p_0}\right)^4 \cdot \beta \cdot e_0^2 \cdot \left(-8 + 40\sin{i_0}^2 - 35\sin{i_0}^4\right) \end{aligned} \\ & \left. \right] \\ & \Omega(t) = \Omega_0 + \left[ \right. \\ & \begin{aligned} \quad -\frac{ 3}{ 2} & \cdot \bar{n} \cdot J_2 \left(\frac{R_0}{p_0}\right)^2 \cdot \cos i_0 + \\ \quad \frac{ 3}{32} & \cdot \bar{n} \cdot J_2^2 \left(\frac{R_0}{p_0}\right)^4 \cdot \cos i_0 \cdot \left(-36 - 4e_0^2 + 48\beta + (40 - 5e_0^2 - 72\beta) \sin{i_0}^2\right) + \\ \quad \frac{15}{32} & \cdot n_0 \cdot J_4 \left(\frac{R_0}{p_0}\right)^4 \cdot \cos i_0 \cdot \left(8 + 12e_0^2 - (14 + 21e_0^2) \sin{i_0}^2 \right) \end{aligned} \\ & \left. \right] \cdot \left(t - t_0\right) \\ & \omega(t) = \omega_0 + \left[ \right. \\ & \begin{aligned} \quad \frac{ 3}{ 4} & \cdot \bar{n} \cdot J_2 \left(\frac{R_0}{p_0}\right)^2 \cdot (4 - 5\sin^2 i_0) + \\ \quad \frac{ 3}{128} & \cdot \bar{n} \cdot J_2^2 \left(\frac{R_0}{p_0}\right)^4 \cdot \left(384 + 96e_0^2 - 384\beta + (-824 - 116e_0^2 + 1056\beta) \sin{i_0}^2 + (430 - 5e_0^2 - 720\beta) \sin{i_0}^4\right) - \\ \quad \frac{15}{ 16} & \cdot n_0 \cdot J_2^2 \left(\frac{R_0}{p_0}\right)^4 \cdot e_0^2 \cdot \cos{i_0}^4 -\\ \quad \frac{15}{128} & \cdot n_0 \cdot J_4 \left(\frac{R_0}{p_0}\right)^4 \cdot (64 + 72e_0^2 - (248 + 252e_0^2) \sin{i_0}^2 + (196 + 189e_0^2) \sin{i_0}^4) \\ \end{aligned} \\ & \left. \right] \cdot \left(t - t_0\right) \\ & M(t) = M_0 + \bar{n} \cdot \left(t - t_0\right) \\ \end{aligned}\]

where the subscript $_0$ indicates the initial element, $t_0$ is the initial mean elements' epoch, $a$ is the mean semi-major axis, $n_0$ is mean motion, $\bar{n}$ is the perturbed mean motion, $e$ is the eccentricity, $p$ is the semi-latus rectum, $R_0$ is the Earth's equatorial radius, $\Omega$ is the RAAN, and $M$ is the mean anomaly.

Warning

This propagator provides the same results as the ANSYS STK® for all Keplerian elements except for the right ascension of the ascending node (RAAN). In this case, we would need to flip the sign of the $J_4$ perturbation term to obtain similar values. However, this modification does not seem right given the references [1] and [2], and the version implemented here leads to lower errors when fitting orbits. Hence, further investigation is required. For more information, see this issue.

Initialization

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

Propagators.init(Val(:J4), orb₀::KeplerianElements; kwargs...) -> OrbitPropagatorJ4

which creates a J4 propagator structure OrbitPropagatorJ4 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:

J4 Propagator Constant	Description	Type
`J4C_EGM2008`	EGM-2008 gravitational constants	`Float64`
`J4C_EGM2008_F32`	EGM-2008 gravitational constants	`Float32`
`J4C_EGM1996`	EGM-1996 gravitational constants	`Float64`
`J4C_EGM1996_F32`	EGM-1996 gravitational constants	`Float32`
`J4C_JGM02`	JGM-02 gravitational constants	`Float64`
`J4C_JGM02_F32`	JGM-02 gravitational constants	`Float32`
`J4C_JGM03`	JGM-03 gravitational constants	`Float64`
`J4C_JGM03_F32`	JGM-03 gravitational constants	`Float32`

Note

The type used in the propagation will be the same as used to define the constants in the structure j4c.

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(:J4), orb; j4c = J4C_JGM03)OrbitPropagatorJ4{Float64, Float64} (J4 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.22855872 rev/day
  │    RAAN Rate              : 0.9555228851 °/day
  │    Arg. of Periapsis Rate : -2.921491706 °/day
  ├─ Constants
  │    R₀ : 6378.1363 km
  │    μm : 0.001239447666 rad/s
  │    J₂ : 0.001082636023
  │    J₄ : -1.619331205e-6
  └─ Propagation
       Last Instant : 0.0 s

Fitting Mean Elements

We can use the function:

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

Note

This algorithm version will allocate a new J4 propagator with the constants j4c. If the allocation must be avoided, use the function Propagators.fit_mean_elements! instead.

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,
           [-1781.214419290065, 1619.7795321872854, 6707.771633846665]    .* 1000,
           [ 5693.643675547716, -1192.342828671633, 4123.976025977494]    .* 1000,
           [ 5291.613719530499, -2354.5417593130833, -4175.561367156414]  .* 1000,
           [-2416.3705905186903, -268.74923235392623, -6715.411357310478] .* 1000,
           [-6795.043410709359, 2184.4414321930635, -0.4327055325971031]  .* 1000,
       ];
julia> vv_i = [
           [0.3445760107690598, 1.0395135806993514, 7.393686131436984]    .* 1000,
           [6.875680282038698, -1.864319399615942, 2.270603214569518]     .* 1000,
           [3.8964090757666496, -2.1887896252945875, -5.9960180359219075] .* 1000,
           [-4.470258022565413, 0.5119576359985208, -5.9608372367141635]  .* 1000,
           [-6.647358060413909, 2.495415251255861, 2.292118747543002]     .* 1000,
           [0.3427096905434428, 1.040125572862349, 7.3936887585116855]    .* 1000,
       ];
julia> vjd = [
           2.46002818657856e6
           2.460028200467449e6
           2.460028214356338e6
           2.4600282282452267e6
           2.4600282421341157e6
           2.4600282560230047e6
       ];
julia> orb, P, stats = Propagators.fit_mean_elements(Val(:J4), vjd, vr_i, vv_i)ACTION:   Fitting the mean elements for the J4 propagator.
           Iteration        Position RMSE        Velocity RMSE           Total RMSE       RMSE Variation
                                     [km]             [km / s]                  [ ]


PROGRESS:          1              53.2399            0.0547615              53239.9                  ---

PROGRESS:          2              4.34243           0.00537588              4342.43             -91.8436 %

PROGRESS:          3              4.33861            0.0053998              4338.61           -0.0880663 %

PROGRESS:          4              4.33863           0.00539961              4338.63          0.000476165 %

(KeplerianElements{MeanAnomaly, Float64, Float64}: Epoch = 2.46003e6 (2023-03-24T18:08:40.388), [0.16604866268470456 0.06643592989410553 … -3.8553207101755296e-5 0.00012403204407738838; 0.06643592989492662 0.26633435566583386 … -1.7942562826191596e-5 -1.9568111153594842e-5; … ; -3.855320709217859e-5 -1.7942562825277312e-5 … 4.397201329008156e-7 -8.092682611441456e-8; 0.00012403204407903322 -1.9568111154084217e-5 … -8.092682612285625e-8 1.2451901446128845e-7], (converged = true, iterations = 4, position_rmse = 4338.627383011516, velocity_rmse = 5.39961217416161, total_rmse = 4338.630743037368))
julia> orbKeplerianElements{MeanAnomaly, Float64, Float64}:
  Epoch             : 2.46003e6 (2023-03-24T18:08:40.388)
  Semi-Major Axis   : 7131.636674 km
  Eccentricity      : 0.001142978092
  Inclination       : 98.43659578°
  RA of Asc. Node   : 162.1773538°
  Arg. of Periapsis : 101.2819928°
  Mean Anomaly      : 258.8210223°
julia> stats(converged = true, iterations = 4, position_rmse = 4338.627383011516, velocity_rmse = 5.39961217416161, total_rmse = 4338.630743037368)

References

[1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications. 4th ed. Microcosm Press, Hawthorn, CA, USA.
[2] Kozai, Y (1959). The Motion of a Close Earth Satellite. The Astronomical Journal, v. 64, no. 1274, pp. 367 – 377.