Earth geopotential effect on orbit propagation

I'm trying to take into account the Earth geopotential effect on the orbit. The only term I could consider, is J2.

1. Are the terms (J2, etc.) change with time?
2. I found other terms coefficients here. How to add them also to differential equation, taken from here?

3.Why does the shown equation give 1 km error in 1 day in comparison with GMAT (configuration is below)? The coefficients for Eq radius and J2 are the same. The integrator is 9th order Runge-Kutta with tol=10e-13.

• Are you sure there is a $J_0$ and a $J_0$? Have you ever seen them? In place of $J_0$ you are probably (I certainly hope) already using the monopole term or $GM_{Earth}$. Something like a $J_1$ or "dipole" term better not exist, as the Geocenter better be the center of mass of the Earth already, thereby zeroing out any dipole term. However, there is a small $J_{22}$ term. I would recommend you search within this site, you may find some parts to answers to your questions here already. – uhoh Jul 17 '18 at 13:20
• ilrs.gsfc.nasa.gov/docs/2014/196C.pdf starting at around page 10 is probably unhelpful, but I linked it just in case. – barrycarter Jul 17 '18 at 15:09
• While J2 is cylindrically symmetric most of the higher order terms move with the Earth, so you have to keep track of that as well. I'm still reading this answer about how to use them. Well okay, I should be reading it, I've gotten behind on some projects these days, but you many find some more helpful information there as well! – uhoh Jul 18 '18 at 10:28
• Re your third question, Why does the shown equation give 1 km error in 1 day in comparison with GMAT? This raises some questions: How did you integration the equations of motion that used $J_2$? What values did you use for $GM_\text{earth}$, $r_\text{eq}$, and $J_2$? (We're taking your word for it that you used the same values; this may be one source of the discrepancy.) Did you turn off polar motion, nutation, and precession in GMAT, if that can be done? – David Hammen Jul 18 '18 at 10:43
• I'm not seeing the values of GM, Re, or J2, or the earth radius in the question. To be consistent with JGM-2 you should be using GM=3.986004415e+14 m^3/s^2, Re=6.3781363e+06 m, and J2=0.0010826360555. However, that you used a different integrator in a different programming language might well be the answer. Not turning off polar motion, nutation, and precession might also be the answer. Note well: A one km/day discrepancy isn't bad. It's not great, but it's not bad. – David Hammen Jul 18 '18 at 10:54

1. Do the terms (J2, etc.) change with time?

Yes, they do. The Earth is still rebounding from the end of the last glaciation, and the Earth's rotation rate is decreasing due to transfer of angular momentum to the Moon's orbit. The end result is that $J2$ is decreasing by about 3 parts per million per century. You don't need to model that, however. The gravity models that GMAT uses are static models. You can enable tides, making the gravity model somewhat dynamic, but the dynamics are purely cyclical. There are no secular terms.

1. I found other terms coefficients here. How to add them also to differential equation?

Find a software package that uses spherical harmonics to model gravitation. Do not roll your own. And if you do use such a package, make sure the coefficients you use are consistent with the model. The coefficients you found are denormalized. Most modern spherical harmonics models expect fully normalized coefficients.

A much trickier issue lies in the $\bar C_{20}$ term. Some of the Earth's equatorial bulge results from tidal interactions with the Moon and the Sun. The spherical harmonics gravitational coefficients are typically computed as if the Moon and Sun were not present. These tide-free models omit the contribution of the Moon and Sun to the equatorial bulge, and hence to $J_2$. (Technically, the frequency coefficients used to model the Earth tides have a zero frequency term that is nonzero.) The Earth spherical harmonics coefficients used by GMAT are tide-free. If you want to model how satellites precess due to the tidal bulge, but don't want to use a full blown Earth tide gravity model, it's better to add a tiny bit to the the $\bar C_{20}$ term. See Section 6, Gravitation of IERS Technical Note 36 for details.

1. Why does the shown equation give 1 km error in 1 day in comparison with GMAT (configuration is below)?

The GMAT JGM-2 gravity model is in the file JGM2.cof of the GMAT source code tree. The first few non-comment lines in this file are

POTFIELD 70 70  1 3.98600441500000e+14 6.37813630000000e+06 1.00000000000000e+00
RECOEF    2  0   -4.84165390000000e-04
RECOEF    2  1   -1.86987640000000e-10 1.19528010000000e-09
RECOEF    2  2    2.43908370000000e-06-1.40010930000000e-06


You need to use compatible values to have your Julia integrator be consistent with GMAT's implementation of JGM-2. The first value, $3.986004415\times10^{14}$, is the TT-compatible value of the Earth's gravitational coefficient in $\text{m}^3/\text{s}^2$. You should be using this rather than the WGS-84 value of $3.986004418\times10^{14}$. The WGS-84 model is aimed at GPS; it is relativistically correct and hence uses a TDB-compatible value of $GM_\oplus$.

The second value, $6.3781363\times10^{6}$, is the equatorial radius of the Earth, in meters. You should be using this rather than the WGS-84 value if you want to try to match the results from GMAT.

The next three lines contain the normalized values of the cosine and sine coefficients of the JGM-2 gravitational potential model for the Earth. The first, -4.84165390000000e-04, is the tide-free value of the fully normalized $\bar C_{20}$ coefficient. The $\bar C_{20}$ is directly related to $J_2$ via $J_2 = -\sqrt{5}\bar C_{20}$.

From poking at the GMAT source code, there is a correction for the $\bar C_{20}$ term gravity coefficient if that term includes the permanent tide and if the user enables Earth tides. There does not appear to be a correction for a tide-free value that adds the permanent tide effect if the user disables Earth tides. The widely published value of $J_2$, 0.0010826359, includes the permanent tides. You probably should be using a value consistent with the tide-free value of $\bar C_{20}$, or $\sqrt{5}\,4.8416539\times10^{-4} \approx 0.00108262672$ instead of that widely published value.

The next two lines contain $\bar C_{21}$, $\bar S_{21}$, $\bar C_{22}$, and $\bar S_{22}$. Note that the $\bar C_{21}$ and $\bar S_{21}$ are nonzero. This means that the Earth's axis of rotation is not quite in line with the way we measure latitude and longitude. The Earth's rotation axis undergoes a small polar motion; you should turn this off in GMAT (if you can) to have a better chance of having your integrated state agreeing with the state computed by GMAT. If you can, you should also turn off Earth nutation and precession in GMAT.

Finally, you're using a different language, and possibly a different numerical integrator. You shouldn't be surprised on seeing differences with exactly the same code when one changes compiler optimization options or migrates to a different computer. You shouldn't be surprised at all when using a different language.

• The problem is not polar motion so much as precession and nutation. You need to be able to disable nutation and precession. If you can't, try setting the start time to 2000 Jan 01 00:00:00 TT and the end time to 2000 Jan 02 00:00:00 TT (i.e., 24 hours centered about noon TT on 1 Jan 2000, which is the J2000 epoch time). The J2000 equatorial frame is the time-smoothed earth orientation at the J2000 epoch time. If GMAT can't handle TT, UTC is close enough. You should see a much smaller error over this time interval compared to 07 Oct 2064 if it's precession and nutation. – David Hammen Jul 18 '18 at 20:08
• What's happening is that for a gravitating object with a non-spherical gravity model, the gravitational acceleration must in general be computed in the planet-fixed frame. GMAT does this by converting the inertial state to the planet fixed frame, computing the acceleration due to gravity, and converting the computed acceleration to the inertial frame. The equations in the question are valid only if the planet fixed frame doesn't undergo precession, nutation, and polar motion. (A rotating planet is okay; the zonal harmonics (e.g. $C_{20}$ or $J_2$) don't depend on longitude.) – David Hammen Jul 18 '18 at 21:17
• The primary source of the error you are still seeing is because the Earth was not aligned with the ICRF frame at noon TT, 01 Jan 2000. The ICRF frame instead represents the alignment of a fictitious mean Earth at noon TT, 01 Jan 2000. – David Hammen Jul 18 '18 at 21:19
• The equations assume a frame in which the z axis is parallel to a planet's rotation axis. Such a frame is not inertial if that planet is undergoing nutation and precession. Another way to look at: The equations give the gravitational acceleration, but expressed in a frame fixed with respect to the planet's rotation axis. You can (in fact, must) calculation the acceleration in that frame, but if the planet is nutating / precessing, integrating the equations of motion in that frame will result in errors. – David Hammen Jul 18 '18 at 21:54
• @Leeloo - That is exactly how you would calculate gravitational acceleration in ECEF, but if you would need to add accelerations due to fictitious forces if you wanted to integrate in ECEF. Fictitious forces vanish in an inertial frame, so transforming the ECEF-based gravitational acceleration to inertial means you can integrate there. What every system that deals with non-point mass gravitating objects does is (1) transform relative position to planet fixed, (2) compute acceleration, (3) transform acceleration to inertial, (4) deal with other forces / other planets, and (5) integrate. – David Hammen Jul 19 '18 at 11:42