Calculating the Earth Geopotential effect

Question

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

  CS=open(readdlm,JGM_coeff_file)

  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)

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.