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!
1 Answer
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.
-
$\begingroup$ +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. $\endgroup$– uhohMar 26, 2018 at 7:41
-
$\begingroup$ Well, geopotential model gives you gravity directly, it is basically differentiated spherical harmonics. $\endgroup$ Mar 26, 2018 at 8:45
-
$\begingroup$ 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. $\endgroup$– uhohMar 26, 2018 at 9:10
-
$\begingroup$ 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. $\endgroup$ Mar 26, 2018 at 10:43
-
$\begingroup$ 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. $\endgroup$– uhohMar 26, 2018 at 15:55
odeintyet, 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. $\endgroup$odeintdoes 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? $\endgroup$