Cálculo del efecto geopotencial de la Tierra

Estoy calculando la aceleración debida a los armónicos en el marco ECEF. Aquí se muestra el potencial de gravedad en armónicos esféricos (acabo de eliminar "$1+$", considerando sólo el efecto de los armónicos).

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

donde$d$y$o$- grado y orden,$\phi$y$\lambda$- latitud y longitud respectivamente.

Comparo los resultados con GMAT. Para grado 2 y orden 0 (J2) el error de propagación fue de 5m. ¡Pero para grado/orden = 8, el error es de 350 km!

Los pasos:

1. en un bucle$i\in[2,2]$,$j\in[0,0]$

   • Calcula el$P_i=\frac{d^{i+j}}{d\mu^{i+j}}(\mu^2-1)^i$
   • Calcular el polinomio de Legendre$P_{ij}=\frac{(1-\mu^2)^{\frac{j}{2}}}{i!*2^i}P_i$
   • Calcula el$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. Calcular el potencial$U_{har}=\mu\frac{sum}{r}$
3. Calcula el$a_x$:$f=\frac{dU_{har}}{dx}$
4. Calcular el valor en el punto$f(r_x,r_y,r_z)$

Como se ve, la latitud es$asin(\frac{z}{r})$y la longitud es$atan(\frac{y}{x})$

Los coeficientes (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

He escrito un código en lenguaje Julia, que construye la expresión (según el grado y el orden).

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)