Расчет эффекта геопотенциала Земли

Я рассчитываю ускорение из-за гармоник в кадре ECEF. Здесь показан гравитационный потенциал в сферических гармониках (я просто убрал "$1+$", так как учитывается только эффект гармоник).

$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})]$,

куда$d$и$o$- степень и порядок,$\phi$и$\lambda$- широта и долгота соответственно.

Я сравниваю результаты с GMAT. Для степени 2 и порядка 0 (J2) ошибка распространения составила 5m. Но для степени/порядка=8 ошибка составляет 350 км!

Шаги:

1. В петле$i\in[2,2]$,$j\in[0,0]$

   • Рассчитать$P_i=\frac{d^{i+j}}{d\mu^{i+j}}(\mu^2-1)^i$
   • Вычислить полином Лежандра$P_{ij}=\frac{(1-\mu^2)^{\frac{j}{2}}}{i!*2^i}P_i$
   • Рассчитать$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. Рассчитать потенциал$U_{har}=\mu\frac{sum}{r}$
3. Рассчитать$a_x$:$f=\frac{dU_{har}}{dx}$
4. Рассчитать стоимость в точке$f(r_x,r_y,r_z)$

Как видно, широта$asin(\frac{z}{r})$а долгота$atan(\frac{y}{x})$

Коэффициенты (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

Я написал код на языке Julia, который строит выражение (в зависимости от степени и порядка).

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)