# 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[1](y[1],y[2],y[3]),
dU[2](y[1],y[2],y[3]),
dU[3](y[1],y[2],y[3])
]
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. – ChrisR Jul 29 '18 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 '18 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. – ChrisR Jul 30 '18 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. – ChrisR Jul 30 '18 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 '18 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. – Leeloo Aug 3 '18 at 8:52