Berechnung des Erd-Geopotential-Effekts

Ich berechne die Beschleunigung aufgrund der Harmonischen im ECEF-Rahmen. Das Gravitationspotential in sphärischen Harmonischen wird hier gezeigt (ich habe gerade "$1+$", da nur der Oberwelleneffekt berücksichtigt wird).

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

wo$d$und$o$- Grad und Ordnung,$\phi$und$\lambda$- Breitengrad bzw. Längengrad.

Ich vergleiche die Ergebnisse mit GMAT. Für Grad 2 und Ordnung 0 (J2) betrug der Ausbreitungsfehler 5 m. Aber für Grad/Ordnung = 8 beträgt der Fehler 350 km!

Die Schritte:

1. In einer Schleife$i\in[2,2]$,$j\in[0,0]$

   • Berechne das$P_i=\frac{d^{i+j}}{d\mu^{i+j}}(\mu^2-1)^i$
   • Berechnen Sie das Legendre-Polynom$P_{ij}=\frac{(1-\mu^2)^{\frac{j}{2}}}{i!*2^i}P_i$
   • Berechne das$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. Berechnen Sie das Potenzial$U_{har}=\mu\frac{sum}{r}$
3. Berechne das$a_x$:$f=\frac{dU_{har}}{dx}$
4. Berechnen Sie den Wert im Punkt$f(r_x,r_y,r_z)$

Wie es gesehen wird, ist der Breitengrad$asin(\frac{z}{r})$und Längengrad ist$atan(\frac{y}{x})$

Die Koeffizienten (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

Ich habe einen Code auf Julia-Sprache geschrieben, der den Ausdruck (je nach Grad und Reihenfolge) aufbaut.

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)