地球ジオポテンシャル効果の計算

私はECEFフレームの高調波による加速を計算しています。球面調和関数の重力ポテンシャルは、ここに示されています。 (私は “$ 1 +
$”を削除しました。これは高調波効果のみを考慮したものです)。

$ U_ {har} = frac { mu} {r} { sum_ {i = 2} ^ d sum_ {j =
0} ^ o( frac {R_ {eq}} { ij}( sin phi)(S_ {ij} sin {j
lambda} + C_ {ij} cos {j lambda})] $、

$ d $と$ o $ – degreeとorder、$ phi $と$ lambda $ –
それぞれ緯度と経度です。

結果を GMAT
と比較します。次数2および次数0(J2)の場合、伝播の誤差は5mでした。しかし、degree/order =
8の誤差は350kmです!

ステップ:

  1. In a loop $iin[2,2]$, $jin[0,0]$

    • Calculate the $P_i=frac{d^{i+j}}{dmu^{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})*(Ccos(j*atan(frac{y}{x}))+Ssin(j*atan(frac{y}{x})))(frac{R_{eq}}{r})^i$
  2. Calculate the potential $U_{har}=mufrac{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)$

見たように、緯​​度は$ 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)
ベストアンサー

アークタンジェントに注意してください

私はJuliaに慣れていませんが、ほとんどの言語で$ arctan(y/x)$は$ [ – pi/2、 pi/2]
$の間の値を返します。これは$ y $か$ x $が正か負かを知らないためです。したがって、経度を次のように計算すると、

$ lambda = arctan( frac {y} {x})= arctan( frac {4617
times10 ^ 3} { – 1709 times 10 ^ 3})= -69.7 ^ circ $$

これは第4象限にあります。これは、実際に経度がどこにあるか$ 180 ^ circ $ offです($ y $が正で$ x
$が負の場合は第2象限)。これはJ2のような$ j = 0 $の項には関係ありませんが、$ sin( theta + 180)=
– sin theta $と$ cos( theta + 180) = – cos(
theta)$です。それは、より多くの用語を追加するときの合意の不一致を説明するでしょう。

arctan2を使用する

ほとんどの言語には、xとyの2つのパラメータをとる特殊な関数(通常は
“atan2(x、y)”と呼ばれる)があります。あなたが持っていない場合は、いくつかのifステートメントを使用して、経度が実際にどの象限にあるかを検討する必要があります。

返信を残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です