I'm calculating the acceleration due to the harmonics in ECEF frame. Gravity potential in spherical harmonics is shown here (I just removed "$1+$", as considering only harmonics effect).

$U_{har}=\frac{\mu}{r}[\sum_{i=2}^d\sum_{j=0}^o (\frac{R_{eq}}{r})^iP_{ij}(\sin\phi)(S_{ij}\sin{j\lambda}+C_{ij}\cos{j\lambda})]$,

where $d$ and $o$- degree and order, $\phi$ and $\lambda$- latitude and longitude respectively.

I compare the results with GMAT. For degree 2 and order 0 (J2) the error in propagation was 5m. But for degree/order=8 the error is 350km!

The steps:

  1. In a loop $i\in[2,2]$, $j\in[0,0]$

    • Calculate the $P_i=\frac{d^{i+j}}{d\mu^{i+j}}(\mu^2-1)^i$
    • Calculate the Legendre polynom $P_{ij}=\frac{(1-\mu^2)^{\frac{j}{2}}}{i!*2^i}P_i$
    • Calculate the $sum+=P_{ij}(\frac{z}{r})*(C\cos(j*atan(\frac{y}{x}))+S\sin(j*atan(\frac{y}{x})))(\frac{R_{eq}}{r})^i$
  2. Calculate the potential $U_{har}=\mu\frac{sum}{r}$
  3. Calculate the $a_x$: $f=\frac{dU_{har}}{dx}$
  4. Calculate the value in the point $f(r_x,r_y,r_z)$

As it seen, the latitude is $asin(\frac{z}{r})$ and longitude is $atan(\frac{y}{x})$

The coefficients (JGM-3):

2 0 -0.10826360229840e-02 0.0
2 1 -0.24140000522221e-09 0.15430999737844e-08
2 2 0.15745360427672e-05 -0.90386807301869e-06

I have written a code on Julia language, which builds the expression (depending on the degree and order).

using SatelliteToolbox
using SymEngine


const date  = DatetoJD(2100,01,01,0,0,0)
const degree = 8

y = [-4617E+03, 1709E+03, -5040E+03]

const Req = 6378136.3
const GMe = 398600.4415E+9

function harmonics(dy,y,dU,date)
  dy= [

function potential()

  @vars x y z myu


  longitude=atan( y/x );

  for i=2:degree
      for j=0:degree

          index=1+j; for ll=2:i-1 index+=ll+1; end

          for k=1:i+j P_i=diff(P_i,myu) end


            U_sum+= P_ij(z/r)*(CS[index,3]*cos(j*longitude)+CS[index,4]*sin(j*longitude))*(Req/r)^i


  return lambdify(expand(diff(U,x)),[x,y,z]),lambdify(expand(diff(U,y)),[x,y,z]),lambdify(expand(diff(U,z)),[x,y,z])

@time harmonics(dy,y,dU,date)
@time harmonics(dy,y,dU,date)
  • $\begingroup$ In what frame is the X, Y, and Z coordinates? They need to be in ECEF. $\endgroup$
    – ChrisR
    Jul 29, 2018 at 18:59
  • 1
    $\begingroup$ +1 Great editing, this looks much better, very nice! I'll try to take a careful look today, thanks for taking the time to dig in with MathJax! $\endgroup$
    – uhoh
    Jul 29, 2018 at 23:17
  • 1
    $\begingroup$ @Leeloo, okay, thanks for the edits. To be honest, I'm not sure how to answer this question sadly, as I've only used the Pines equations to compute the harmonics. $\endgroup$
    – ChrisR
    Jul 30, 2018 at 0:03
  • $\begingroup$ @Leeloo the equations are several pages long. A good reference is Chapter 2 of the PhD dissertation of Brandon A. Jones 2010. It's available freely on the University of Colorado website. The filename should be called "bajones2010" I think. He also mentions two other formulations which should speed up calculations (Fantino 2008), but they seemed significantly more complex to me. I've implemented Pines here, in Rust, as part of a toolkit I'm writing to automate most of the analysis I do, & with the fidelity of GMAT. $\endgroup$
    – ChrisR
    Jul 30, 2018 at 18:05
  • $\begingroup$ @Leeloo I'll take a look, that looks very promising! Possibly helpful here as well: How can I verify my reconstructed gravity field of Ceres from spherical harmonics? $\endgroup$
    – uhoh
    Aug 1, 2018 at 8:20

1 Answer 1


Be careful with arc tangent

I am unfamiliar with Julia, but in most languages $\arctan(y/x)$ returns a value between $[-\pi/2, \pi/2]$. It does this because it doesn't know if $y$ or $x$ is positive or negative. So when you calculate longitude as:

$$ \lambda = \arctan(\frac{y}{x}) = \arctan(\frac{4617\times10^3}{-1709\times 10^3}) = -69.7^\circ$$

That is in the 4th quadrant, which is $180^\circ$ off from where the longitude actually is (the 2nd quadrant when $y$ is positive and $x$ is negative). This doesn't matter for the $j=0$ terms like J2, but it inverts the tesseral terms, since $\sin( \theta + 180 ) = - \sin \theta$ and $\cos( \theta + 180) = -\cos ( \theta )$. That would explain the discrepancy in agreement when adding more terms.

Use a arctan2

Most languages have a special function (usually called "atan2(x,y)" or some variant) that takes two parameters: an x and y parameter. If you don't have one, you will have to use some if statements to consider which quadrant the longitude is actually in.

  • $\begingroup$ You're absolutely right. However, there're some troubles with atan2 in Julia, so I used longitude=atan(y/x)+pi*sign(y)*(1-sign(x))/2. Is it OK? Assumed, that x and y are never 0. $\endgroup$
    – Leeloo
    Aug 3, 2018 at 8:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.