J2 Analytical Orbit Propagator

The J2 analytical orbit propagator considers only the $J_2$ perturbation term of the geopotential field to obtain a simple analytic solution for propagating the mean orbital elements. It has low precision but is extremely useful during mission design. For example, when designing a Sun-Synchronous orbit for a remote sensing mission, the satellite usually has a propulsion subsystem to maintain its orbit. Hence, the J2 propagator solution can be used as the nominal orbit through the mission lifetime to perform initial studies.

Algorithm

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

After neglecting all terms with higher order than J2, 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) \\ \bar{n} &= n_0 \cdot \left[1 + \frac{3}{4} \cdot J_2 \cdot \left(\frac{R_0}{p_0}\right)^2 \cdot \sqrt{1 - e_0^2} \cdot \left(2 - 3\sin{i_0}^2\right)\right] \Omega(t) &= \Omega_0 - \frac{2}{3} \cdot J_2 \cdot \left(\frac{R_0}{p_0}\right)^2 \cdot \bar{n} \cdot \cos{i_0} \cdot \left(t - t_0\right) \\ \omega(t) &= \omega_0 + \frac{3}{4} \cdot J_2 \cdot \left(\frac{R_0}{p_0}\right)^2 \cdot \bar{n} \cdot \left(4 - 5\sin{i_0}^2\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.

Initialization

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

Propagators.init(Val(:J2), orb₀::KeplerianElements; kwargs...) -> OrbitPropagatorJ2

which creates a J2 propagator structure OrbitPropagatorJ2 with the mean Keplerian elements orb₀.

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

• j2c::J2PropagatorConstants: J2 orbit propagator constants (see J2PropagatorConstants). (Default = j2c_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{Float64, Float64}:
           Epoch :    2.45995e6 (2023-01-01T00:00:00)
 Semi-major axis : 7190.98     km
    Eccentricity :    0.001111
     Inclination :   98.405    °
            RAAN :  100.0      °
 Arg. of Perigee :   90.0      °
    True Anomaly :   19.0      °
julia> orbp = Propagators.init(Val(:J2), orb; j2c = j2c_jgm03)OrbitPropagatorJ2{Float64, Float64}:
   Propagator name : J2 Orbit Propagator
  Propagator epoch : 2023-01-01T00:00:00
  Last propagation : 2023-01-01T00:00:00

Fitting Mean Elements

We can use the function:

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

to fit a set of mean Keplerian elements for the J2 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 and the final covariance matrix of the least-square algorithm.

The following keywords are available to configure the fitting process:

• atol::Number: Tolerance for the residue absolute value. If the residue is lower than atol at any iteration, the computation loop stops. (Default = 2e-4)
• rtol::Number: Tolerance for the relative difference between the residues. If the relative difference between the residues 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_perturbation::Number: Initial state perturbation to compute the finite-difference when calculating the Jacobian matrix. (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 it absolute value is higher than jacobian_perturbation_tol. (Default = 1e-7)
• max_iterations::Int: Maximum number of iterations allowed for the least-square fitting. (Default = 50)
• mean_elements_epoch::Number: Epoch for the fitted mean elements. (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 = Propagators.fit_mean_elements(Val(:J2), vjd, vr_i, vv_i)ACTION:   Fitting the mean elements for the J2 propagator.
           Iteration        Position RMSE        Velocity RMSE           Total RMSE       RMSE Variation
                                     [km]             [km / s]                  [ ]                     


PROGRESS:          1              53.2569            0.0547793              53256.9                  ---

PROGRESS:          2              6.09001           0.00629975              6090.01             -88.5648 %

PROGRESS:          3              6.28203           0.00643173              6282.03              3.15303 %

PROGRESS:          4              6.21666           0.00639426              6216.66             -1.04058 %

PROGRESS:          5              6.14541           0.00635329              6145.41             -1.14615 %

PROGRESS:          6              6.06795           0.00630898              6067.96             -1.26036 %

PROGRESS:          7              5.98394           0.00626121              5983.94             -1.38464 %

PROGRESS:          8                5.893           0.00620986              5893.01             -1.51959 %

PROGRESS:          9              5.79487           0.00615486              5794.87              -1.6653 %

PROGRESS:         10              5.68932           0.00609619              5689.32             -1.82137 %

PROGRESS:         11               5.5763           0.00603395               5576.3             -1.98663 %

PROGRESS:         12              5.45592           0.00596836              5455.92             -2.15869 %

PROGRESS:         13              5.32862           0.00589979              5328.62             -2.33329 %

PROGRESS:         14              5.19523           0.00582885              5195.23             -2.50332 %

PROGRESS:         15              5.05717           0.00575645              5057.17              -2.6574 %

PROGRESS:         16              4.91669           0.00568391              4916.69             -2.77785 %

PROGRESS:         17              4.77715           0.00561303              4777.16             -2.83799 %

PROGRESS:         18              4.64344           0.00554626              4643.45              -2.7989 %

PROGRESS:         19              4.52243           0.00548689              4522.43             -2.60623 %

PROGRESS:         20              4.42342           0.00543916              4423.42             -2.18922 %

PROGRESS:         21              4.35859           0.00540852              4358.59             -1.46571 %

PROGRESS:         22               4.3413           0.00540076              4341.31            -0.396498 %

PROGRESS:         23               4.3413           0.00540076              4341.31         -2.90569e-06 %

(KeplerianElements{Float64, Float64}: Epoch = 2.46003e6 (2023-03-24T18:08:40.388), [0.16604846210230226 0.06643574127302715 … -3.855413642387002e-5 0.0001240322540313337; 0.0664357412698984 0.2663344825561057 … -1.7943611854657147e-5 -1.956795689447386e-5; … ; -3.85541364223557e-5 -1.7943611852843348e-5 … 4.3971984252736047e-7 -8.092704667554433e-8; 0.0001240322540325942 -1.956795689122344e-5 … -8.092704667640543e-8 1.2451922443375673e-7])
julia> orbKeplerianElements{Float64, Float64}:
           Epoch :    2.46003e6 (2023-03-24T18:08:40.388)
 Semi-major axis : 7131.63       km
    Eccentricity :    0.00114299
     Inclination :   98.4366     °
            RAAN :  162.177      °
 Arg. of Perigee :  101.286      °
    True Anomaly :  258.689      °

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.