|
ここでは数値積分による彗星の軌道計算のまとめとして、ハレー彗星の数値積分をやってみたいと思います。ハレー彗星については過去から検討されていますが、次回の回帰が５０年余り先とのこともあり、WEBで探してもあまり新しいまとまった研究は無いようです。なのでここは前回の回帰(1986年)の前に発表されているD.K.Yeomans/T.Kiang著「The long-term motion of comet Halley」(1981)をもとに検討してみました。
この論文を読むと、D.K.Yeomans/T.Kiang両氏は古代についてはAD1607年からBC1404まで一貫して計算しているとのことですが、837/04/11の地球への接近(0.04AU)を調整するために、837年の近日点通過時刻を-0.88日変更し、また軌道離心率を-7.2x10-6だけ変更しているとのこと。またこの時代なので惑星暦(DE97の初期値を使用)は4日間隔で400年毎にテープに収められていたということなので、時間もかかっていたのではないかと思われます。数値積分としてはアダムス法を使い、地球と月は重心位置で考え、他に8惑星も考慮していたとのこと。
彗星の計算になると非重力効果を考慮するのが重要となるので、「天体軌道論」(p329-330)を参考に非重力効果のルーチンをプログラムに組み込みました。なおA3及びB1/B2/B3は無視(0)しています。論文では期間中(AD1607年からBC1404)A2=0.2767/A2=0.015を使用したとされていますが、ここでは「Catalogue of cometary orbits 1997(12th)」にある切り上げた数値A1=0.28/A2=0.015をベースに考えています。
以上のような状況を考慮しここでは以下を目標に計算してみました。
①1986年の回帰から古代まで一貫して計算する。（非重力定数以外は変更しない）
②「The long-term motion of comet Halley」(1981)の近日点通過の計算結果になるべく近づける。
（中国の古代の観測を基にした観測推定値とは近日点通過で最大数日の差がある。）
③論文ではBC1404まで計算しているが、記録とされるものは-239年までしかない。
したがい、あまり遡っても意味が無いので、-239年までを近づけることを目標とした。
④AD1600年以降については観測値を使い最適化する必要があるがここでは周期のみに着目。
⑤積分間隔は論文と同じ0.5日間隔とする。（計算時間はAD1986-BC1000で約10分間）
⑥近日点通過日は正確には軌道要素を計算して出さないといけないが、ここでは簡易的に
最接近の日付とした。（刻みは0.5日）
数値積分を行ううえでのハレー彗星の初期値は以下の軌道要素から計算した。計算開始は1986年の近日点通過の2日後を初日としそこから過去へ遡る計算を行う。
Comet Halley Orbits from catalogue of cometary orbits 1997
ORB[1]= 2446470.95890// T(TPASS) :近日点通過時刻 (1986/ 2/10)
ORB[2]= 111.8657000 // ω(PERI/Somg) :近日点引数
ORB[3]= 58.86010000 // Ω(NODE/Omg) :昇交点黄経
ORB[4]= 162.2422000 // i (INC) :軌道傾斜角
ORB[5]= 0.58710400 // q (QDIS) :近日点距離
ORB[6]= 0.96727700 // e (ECC) :軌道離心率
EPOCH = 2446472.500 // EPOCH :元期
以上の軌道要素からプログラムで使用した初期値（位置/速度）
CD[I][16][1]= 0.292488509596; CD[I][16][2]=-0.508161688283; CD[I][16][3]=-0.045506195962;
VL[I][16][1]=-0.025373343081; VL[I][16][2]=-0.015110455835; VL[I][16][3]=-0.010850178120;
AMASS[16]=0.0;// 重量
この初期値をもとに計算した結果が以下。以下の図は計算結果と論文にある近日点通過日からのずれの日数をしめしたもので、青色の線が非重力効果を考慮しない場合。赤色の線が「catalogue of cometary orbits 1997」にある非重力効果の定数をそのまま使った場合の線です。一回の周期の違いが累計して効いてくるので、非重力効果を考慮しない場合、紀元前後では３年以上違ってきてしまっています。また考慮した赤い線の場合でも-239年あたりでは一年の違いがでてきてます。
そこで非重力効果の定数A1を調整し、それぞれの近日点通過の日時のずれが１日前後になるようにラフに調整したのが黄色の線です。これからいくと論文のように離心率等を調整しなくても、非重力効果の定数のわずかな変化で説明できる範囲のようです。
地球への最接近としては論文にあるように837/4/10に0.033AUまで近づいていました。(論文では837/04/11に0.04AU)
図１ D.K.Yeomans/T.Kiang論文の近日点通過日からのずれた日数
表１ 近日点通過日からのずれた日数 (計算の刻みは0.5日）
(2014/10/29追加)
またそれぞれの場合の回帰周期をまとめたものが図２です。こえを見ると現代からAD374年ぐらいまでは、周期に目立った違いは見られず、これから非重力効果自体は年単位の変化を与えるものではないことが分かります。しかし、少しの違いの累積が惑星との位置関係に違いを与え、惑星から受ける摂動が大きくなったり、逆に小さくったりするため、年単位で周期が変わる原因になっているようです。
図２ ハレー彗星の回帰周期(年)の推移
以上の内容にて作成したソフトで彗星の軌道計算にも問題は無いようです。
(念のため、0.125日間隔でもやってみましたが-390年の回帰まではほぼ同じ結果でした。)
ハレー彗星の計算ソフトのC言語ソース(2015/06/26修正版UP)：
n_berg3_2.0.zip
注：windows上のCコンパイラ(VC等)でコンパイルください。
JPL/NASAのDE431のファイル(lnxm13000p17000.431)がほかに必要です。
|