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!
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.