1. 惑星と小惑星のN体数値積分




 彗星の数値積分をおこなうためには、まず一般的な計算ルーチンが必要ということで、まずは惑星(太陽・地球・月を含む)と 小惑星のプログラムを作ってみることにした。

 最初に「天体の軌道計算」にあるMCDというFortranで書かれたソフトをCに移植して計算してみた。MCDはDE28の パラメータがベースになっており、修正しないそのままのソフトで一番古い1804年(jd=238000.5)日の計算結果でも小数点5〜6桁の誤差。しかし初期値が1950年分点で、決算結果は2000年分点に変更しているため使い勝手が悪かった。
 それなら最初からJ2000のDE430/431の初期値を与えればそのまま使えると考えたのだが、こんどはどうしたことか結果がオリジナルの計算値にくらべ2桁ぐらい悪くこれ以上遡った計算には問題があると思われた。ここで1,2週間いろんなことを試してみたがこの差は埋まらなかった。

 そこでDE200やDE430/431の文献をみてみると、ニュートン力学の摂動力(=の次の項)だけでなく、相対性理論からくる力(1/C2がかかっている項)が働いていることが判明。具体的には以下の(27)式になる。(実際にはこれ以外にも潮汐力とか他の力が働いているようだが、彗星の計算にはとりあえず効いてこなそうなの感じなので略した。*1) この式はベクトルの内積とかが表現されているようだが一見してわからず、昔NIFTYでダウンロードしていたNIPEのソースを見てようやく計算方法が理解できた。
*1 最終的には彗星の数値積分には太陽/惑星/月の位置/速度はDE431のデータを使用するので。

「The Planetary and Lunar Ephemerides DE430 and DE431」William M. Folkner,* James G. Williams, Dale H. Boggs, Ryan S. Park, and Petr Kuchynka(2014/02/15)より(なお、βとγは1です)


 この式を計算してみた結果が以下の2つの表にまとめてある。上の表がそのまま計算したもの。下の表がニュートン力学(最初の項)のみの結果。これでようやく使える精度となった。(小惑星については比較するものがないので表には載せていません。)
計算しているのは太陽・9大惑星(冥王星を含む)および5個の小惑星である。なお数値積分は誤差の見積もりもあるのでルンゲクッタ法の一種であるフェールベルク法(6段5次)を使用し積分間隔は0.125日。 (計算時間は、40日間隔出力で1分程度)


 誤差が使用できる領域となったので、そのまま過去3000年を計算して、その一日毎の誤差の平均(X,Y,Zの誤差の2乗の和のルート)を計算したのが以下のグラフ。惑星としては悪い水星や地球でマイナス4乗(AU)程度。


 最後に参考までにこのソフトのソースを添付します。昔の古いCのテンプレートを使用して「MS Visual studio 2008」でコンパイルしています。 最初からVisual Studioでやる場合には.cppに拡張子を変える等c++環境にする必要があるかもしれません。
またgccでもコンパイルできることを確認済みです。gccではいかのような感じでコンパイル/実行してみてください。(-mlは数値演算のオプションです。)

>gcc -o n_bergu n_bergu.c -ml
>./n_bergu

N体惑星数値積分ソフト(N_Bergu):
n_bergu_2.0.zip


2014/08/30 Up
2015/05/16 修正版Up

Copyright(C) 2014 Shinobu Takesako
All rights reserved