# Calculating the Earth Geopotential effect

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

path="C:/xampp/htdocs/sat_prop/";
JGM_coeff_file=string(path,"coeff.txt");

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= [
dU(y,y,y),
dU(y,y,y),
dU(y,y,y)
]
end

function potential()

@vars x y z myu

longitude=atan( y/x );
r=(x^2+y^2+z^2)^(1/2)

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

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

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

P_ij=(((1-myu^2)^(j/2))/(factorial(i)*2^i))*P_i

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

U=GMe*(U_sum)/r

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])
end

dU=potential();
dy=zero(y)
@time harmonics(dy,y,dU,date)
@time harmonics(dy,y,dU,date)

• In what frame is the X, Y, and Z coordinates? They need to be in ECEF. Jul 29, 2018 at 18:59
• +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!
– uhoh
Jul 29, 2018 at 23:17
• @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. Jul 30, 2018 at 0:03
• @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. Jul 30, 2018 at 18:05
• @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?
– uhoh
Aug 1, 2018 at 8:20

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

• 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. Aug 3, 2018 at 8:52