Discrepancies between J2 Propagator and STK J2Perturbations

[xinlongliu]

Hi, I am using this package to simulate a J2 propagator. However, there are discrepancies between the results and STK J2Perturbations. The discrepancy can be around 1.5 km after 5 days. Did anyone run into this issue? Are there any investigations on why there are such discrepancies?
(I have tried TwoBody Propagator, and the results are the same as STK TwoBody. So I assume there are no mistakes about how I use the package.)
Thank you very much in advance!

[ronisbr]

Hi @xinlongliu !

Can you provide a detail example of the STK configuration and SatelliteToolbox.jl configuration?

[xinlongliu]

Thank you for your quick reply!
Sure, below is the command that I used to call J2 Propagator:
SatelliteToolbox.j2!(SatelliteToolbox.j2_init(2459945.5, 8000000.0, 0.015, 28.5pi/180, 100pi/180, 200pi/180, 45pi/180), 345600.0)
The screenshot below is the STK J2 configuration:

Thank you!

[ronisbr]

Hi @xinlongliu !

Can you please provide a CSV file with the results from STK?

[xinlongliu]

Sure, attached please find the CSV file with the results from STK J2 Perturbation. Thanks!
Satellite1_Inertial_Position_Velocity.csv

[ronisbr]

Thanks! I will analyze the file. Is there a way in STK to obtain the values of the constants it is using to propagate the orbit? Like the gravitational parameter of the Earth and the J2 coefficient?

[xinlongliu]

Thanks! I only know that STK uses WGS84 (https://agiweb.secure.force.com/faqs/articles/Keyword/Can-the-reference-ellipsiod-used-by-STK-be-changed), but I cannot find the exact values for those constants. I have tried several combinations of the popularly-used values, but none of them got the results same as STK.

[ronisbr]

Thanks @xinlongliu! It turns out it was a bug when computing the "mean" mean motion to obtain the perturbations. Now, the difference between the STK version and the propagator here is less than 5 meters per component after 4 days. This difference is probably caused by the constants used in STK together with numerical errors.

Can you please test again master?

[xinlongliu]

Hi @ronisbr, I have tested it, and it works very well with less than 5 meters after 4 days. This is amazing. Thank you so much!
Also, is it possible to change for the J4 propagator as well? It is also different from STK. Maybe it is because of the same issue?
Attached please find the results from STK J4 Perturbation with the same input parameters.
Satellite1_Inertial_Position_Velocity_J4.csv
Thank you again!

[ronisbr]

Hi @xinlongliu !

I did the same modifications, but it is still different. I am researching to check what it is implemented in STK. I am using the algorithms in Vallado’s book. Do you have any ideas?

[xinlongliu]

I see, thanks! I do not have any ideas so far. I am new in this field, and I am still learning. I will do some research as well, and will let you know if I have any findings. Thank you!

[xinlongliu]

At the end of p693 and at the beginning of p694 in [1] Vallado, D. A (2013), it is mentioned that J2 squared term is different due to the method of averaging used. Is this information useful?

[ronisbr]

Hi @xinlongliu ! Thanks! I implemented both modifications (Kozai and Brouwer), but it is still very wrong. My current approach is to understand Kozai's paper because I think I am using wrong "mean" values in the perturbations for the semi-major axis and mean motion.

[ronisbr]

Hi @xinlongliu ! The STK report allows to generate a table with the classical orbit elements. Can you please send them to me using the J4 propagator but with 6 decimal digits?

[xinlongliu]

Sure, please find the attached csv file. Thanks!
Satellite1_Classical_Orbit_Elements_J4.csv

[ronisbr]

Thanks! I see that the Mean Anomaly is being compute just as STK. However, RAAN and the argument of perigee are very wrong. I am still investigating.

[ronisbr]

Hi @xinlongliu !

I need one more thing :) Sorry for the number of requests. I am almost concluding my analysis. I might have found a bug somewhere (in STK or Vallado's book). When I implement the equations using Kozai's method, I get very close values for the mean anomaly, and argument of perigee. However, for the RAAN, I need to flip the sign of the J4 perturbation term.

To verify if I am right, can you please provide me those two files (classical elements and RV) for a completely different orbit?

[xinlongliu]

Sure, great to hear that! Please see the attached files and screenshot of the orbital parameters. Thanks!

Satellite1_Classical_Orbit_Elements_J4_2.csv
Satellite1_Inertial_Position_Velocity_J4_2.csv

[ronisbr]

Hi @xinlongliu !

Thanks! I am starting to think that there is a bug in STK J4 propagation. Let me discuss it here.

If I implement exactly the same equations in Vallado's book using Kozai's method with only the secular perturbations, the result is equal to the STK's up to the third decimal digit for the mean anomaly and the argument of perigee. However, for the RAAN, the result is very different. If I flip the sign of the J4 perturbation term in the RAAN, then the result is equal to the STK's up to the third digit.

Here is a table comparing the orbital elements (only those that change) between my implementation and STK's for the first orbit:

	STK	Vallado's book	With J4 pert. sign flipped
RAAN [°]	84.1588	84.1134	84.1588
Arg. of perigee [°]	225.864	225.862	255.862
Mean anomaly [°]	247.19	247.193	247.193

And here are the results for the second orbit:

	STK	Vallado's book	With J4 pert. sign flipped
RAAN [°]	29.6424	29.6266	29.6424
Arg. of perigee [°]	171.583	171.589	171.589
Mean anomaly [°]	127.991	127.989	127.989

Since it worked for those two different orbits, it might not be just one coincidence.

The first thing I thought was that there is a typo in Vallado's book. However, the equation for the J4 perturbation term in RAAN is:

The SGP4 equation of the secular perturbation for RAAN is:

where θ = cos(i_0), β = sqrt(1 - e^2), and k₄ = -3 / 8 * J4.

If we ignore the terms with e^2 in Vallado's equation, we see that both algorithms provide exactly the same values. Hence, I think Vallado's equations have the correct sign and STK might have a bug.

I tried to find another analytical J4 propagator implementation but I could not.

I have no idea what to do here. I prefer to do not flip the sign since we do not have any theoretical argument, only the STK's results. We need more help now :)

Hi @helgee !

Sorry to ping you in a random thread. We are having some problems related with the orbit propagator. Can you please give us some help here?

[helgee]

No problem, I'll have a look ASAP.

[ronisbr]

Hi @xinlongliu and @helgee ,

I am pretty convinced that the version equal to the Vallado's book is the correct one and STK is providing wrong results for the J4 orbit propagator. Here is a summary of how I achieve those conclusions.

First, I implemented a numerical orbit propagator using the EGM2008 gravity model up to degree 4 and order 0. The code can be seen here: https://github.com/JuliaSpace/PropagatorTests/blob/main/J4OsculatingPropagator/j4osc_tests.jl

If we select the following mean elements:

• Epoch: 2023-01-01T00:00:00
• Semi-major axis: 7190.982 km
• Eccentricity: 0.001111
• Inclination: 98.405°
• RAAN: 90.0°
• Arg. of Perigee: 200.0°
• True Anomaly: 45.0°

we can obtain the initial osculating state vector by propagating the orbit using the J4 osculating propagator to the instant 0:

jd₀ = DateTime("2023-01-01") |> datetime2julian
r₀  = @SVector [-952882.6807035431, -3.03843825141778e6, -6.444903144699051e6]
v₀  = @SVector [-460.11262511481317, 6745.426195633091, -3115.9662885215403]

Notice that the RAAN perturbation we are verifying here does not have influence in this computation.

After propagating this state vector using a numerical propagator considering only the same perturbation terms in the J4 osculating propagator, we obtain the following result:

Numerical Propagation Results
============================================================================================

Gravity model         : EGM-2008 (Degree 2, Order 0)
Integration algorithm : AutoVern7(Rodas5())

 Time [s] │ Pos. X (TOD)  Pos. Y (TOD)  Pos. Z (TOD)  Vel. X (TOD)  Vel. Y (TOD)  Vel. Z (TOD)
          │           km            km            km        km / s        km / s        km / s
──────────┼────────────────────────────────────────────────────────────────────────────────────
      0.0 │     -952.883     -3038.438     -6444.903        -0.460         6.745        -3.116
    600.0 │    -1034.273      1320.077     -6994.342         0.197         7.315         1.342
   1200.0 │     -731.133      5188.075     -4936.908         0.781         5.164         5.295
   1800.0 │     -156.149      7126.689     -1039.323         1.074         1.090         7.277
   2400.0 │      476.940      6413.690      3245.955         0.968        -3.389         6.546
   3000.0 │      932.849      3316.832      6322.403         0.503        -6.601         3.380
   3600.0 │     1042.660     -1010.885      7047.392        -0.149        -7.361        -1.041
   4200.0 │      765.300     -4962.470      5149.030        -0.747        -5.385        -5.085
   4800.0 │      202.786     -7063.260      1327.389        -1.068        -1.388        -7.242
   5400.0 │     -435.506     -6521.694     -2991.484        -0.991         3.135        -6.686
   6000.0 │     -910.997     -3540.646     -6189.305        -0.543         6.478        -3.629

Now, if we run the J4 osculating propagator considering the equation in Vallado's book, the result for 6000s after the epoch is:

julia> using SatelliteToolboxPropagators, Dates

julia> orb = KeplerianElements(
           DateTime("2023-01-01") |> datetime2julian,
           7190.982e3,
           0.001111,
           98.405 |> deg2rad,
           90     |> deg2rad,
           200    |> deg2rad,
           45     |> 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 :   90.0      °
 Arg. of Perigee :  200.0      °
    True Anomaly :   45.0      °

julia> orbp = Propagators.init(Val(:J4osc), orb)
OrbitPropagatorJ4Osculating{Float64, Float64}:
   Propagator name : J4 Osculating Orbit Propagator
  Propagator epoch : 2023-01-01T00:00:00
  Last propagation : 2023-01-01T00:00:00

julia> r, v = Propagators.propagate!(orbp, 6000)
([-910993.7353742868, -3.5406569550705e6, -6.189314375143491e6], [-543.3929780681474, 6478.262059750446, -3629.151466375647])

julia> norm(r - [-910.997,     -3540.646,     -6189.305] .* 1000)
14.783932703269324

Hence we have an error of less than 15m that is caused probably due to the long-term perturbations that are not considered in the algorithm.

Now, if we flip the sign in RAAN to match STK result, we achieve:

julia> using SatelliteToolboxPropagators, Dates

julia> orb = KeplerianElements(
           DateTime("2023-01-01") |> datetime2julian,
           7190.982e3,
           0.001111,
           98.405 |> deg2rad,
           90     |> deg2rad,
           200    |> deg2rad,
           45     |> 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 :   90.0      °
 Arg. of Perigee :  200.0      °
    True Anomaly :   45.0      °

julia> orbp = Propagators.init(Val(:J4osc), orb)
OrbitPropagatorJ4Osculating{Float64, Float64}:
   Propagator name : J4 Osculating Orbit Propagator
  Propagator epoch : 2023-01-01T00:00:00
  Last propagation : 2023-01-01T00:00:00

julia> r, v = Propagators.propagate!(orbp, 6000)
([-910976.5188533106, -3.540661384752601e6, -6.189314375143491e6], [-543.4244787313135, 6478.25941741541, -3629.151466375647])

julia> norm(r - [-910.997,     -3540.646,     -6189.305] .* 1000)
27.277487013571008

Which has almost 2x the error of the other version. Thus, it seems the version as in Vallado's book is much more precise than the results from STK.