# Calculating J2 perturbations in cartesian coordinates?

I am an engineering student and we have taken upon a project that involves modelling the motion of earth satellites. So far we have a model that uses numerical symplectic integration to solve the Newtonian gravitational equations in cartesian coordinates, giving us the motion of a satellite. We wish to include the pertubation due to the oblateness of the Earth. I am aware that it makes more sense to do this using keplerian orbital elements, but as the majority of our code thus far is based on cartesian coordinates, I am wondering whether it is possible at all to include pertubations using cartesian coordinates, or if it is easier for us to just back track and redesign our model to use keplerian elements. Clearly we have no idea what we are doing, so go easy on us!

• It's a great question! Please check out this answer where I've shown how propagating a Sun-synchronous orbit for 91 days at n=14 orbits per day nicely precesses by 90 degrees, just using simple Python and default SciPy integrator, which I'm assuming only uses symplectic methods. If you are feeling more adventurous, you can look at this and this answer as well. Leave a comment or update to your question if there is something more you need.
– uhoh
Mar 25, 2018 at 19:40
• You can also poke around on this site using search for other Q&A, as well as Astronomy SE and Physics SE as well. Definitely stick with cartesian coordinates, this is how the "pros" have been doing it ever since enough computing power looked like this or this. For more insight on J2 see @DavidHammen's nice table: physics.stackexchange.com/a/141981/83380
– uhoh
Mar 25, 2018 at 19:42
• ...as well as iag-aig.org/attach/e354a3264d1e420ea0a9920fe762f2a0/… and also en.wikipedia.org/wiki/…
– uhoh
Mar 25, 2018 at 19:46
• I haven't tried it with odeint yet, but @jit worked nicely for me as described here and it should probably be just as effective in this case as well, and for the same reasons.
– uhoh
Mar 25, 2018 at 19:58
• update: No, SciPy's odeint does not use symplectic methods! See the excellent answers to What does “symplectic” mean in reference to numerical integrators, and does SciPy's odeint use them?
– uhoh
Mar 26, 2018 at 15:47

Great question.

The perturbations due to non spherical nature of earth is accounted using spherical harmonics which are the general solution of laplace equations. Celestial bodies such as Earth, venus, moon and mars have their geopotential models defined by zonals and tesserals terms, measured by NASA with their probes. Geopotential Model in this wiki page explains the math behind this.

So, most of the time 30 x 30 model suffice for orbit propogations. Greater the complexity, the closer you are to the exact gravity( Well on average though)

So, GeographicLib is an excellent library, which has models for earth gravity inbuilt, which can then be used with appropriate integrator( Runge-Kutta suffices) to get your orbit.

If you want to do this propogations for other bodies than Earth, then you must download the Snm and Cnm variables( Zonals and Tesserals) from appropriate website and then use the library spherical harmonics calculator class to calculate gravity.

Edit: As pointed out by uhoh, how does one gets gravity from gravity scalar potential?

The geopotential model has two coefficients ( zonals and tesserals ) It is intuitively like 3D fourier transform. Now once you have function that gives you scalar potential, it is then differtiated in spherical coordinates to get force function. This has already a known analytic standard form. If you use that library, you dont have to do any differentiation, the library has no function for calculating scalar potential( no use of it, isnt it?) But directly the force function.

• +1 the GeographicLib looks really useful, great link! But one can't use a gravitational potential directly in an integrator. Can you add to your answer something explaining how to determine an acceleration vector from this scalar potential? See my comments and links below the question.
– uhoh
Mar 26, 2018 at 7:41
• Well, geopotential model gives you gravity directly, it is basically differentiated spherical harmonics. Mar 26, 2018 at 8:45
• In order to better answer the OP's question as asked, it would be helpful to describe how one gets a gravitational acceleration vector from a gravitational scalar potential. I'm not sure what "gives you gravity" means.
– uhoh
Mar 26, 2018 at 9:10
• I added bit of explanation, the math is lightly bit involved to be put here. Interested reader can refer wiki and further, references in wiki for rigorous derivations. Mar 26, 2018 at 10:43
• Consider adding the word gradient somewhere! $\mathbf{a_G}=-\nabla \phi$. You need to take the gradient of the potential described by the geopotential model in order to obtain the acceleration vector. Since it is expressed in spherical harmonics, that should be fairly easy to do analytically.
– uhoh
Mar 26, 2018 at 15:55