今回は3番目の要素変化法による軌道要素の最適化です。
@地心太陽座標の計算
A軌道要素の計算
B要素変化法による軌道要素の最適化
計算例は引き続きサイディング・スプリング(C/2013A1)彗星です。
計算内容の詳細は以下の本を参照下さい。
・「天体の軌道計算」中野主一著、誠文堂新光社(1992)
・「マイコン天文学 I」中野主一、恒星社厚生閣(1983)
・「天体軌道論」長谷川一郎、恒星社厚生閣(1983)
今回は「天体の軌道計算」第7章に記載のプログラムルーチン(Fortran版)の必要部分をC言語に移植したルーチンを使用しています。(光度関係は割愛しています。)なお「天体の軌道計算」第7章では観測値を配列で保持していますが添付のプログラムではファイルで保持するようにしています。
@ベクトルの計算
次の7通りの軌道要素によるベクトルの計算を行います。
1)オリジナルの軌道6要素
2)T のみ0.001変化させた軌道要素
3)ωのみ0.00001変化させた軌道要素
4)Ωのみ0.00001変化させた軌道要素
5)iのみ0.00001変化させた軌道要素
6)qのみ0.0000001変化させた軌道要素
7)eのみ0.0000001変化させた軌道要素
A残差の計算
@で計算したそれぞれのベクトル要素をもとに、観測時刻における座標を計算し、観測座標との残差をもとめ残差の配列に代入。
B変化分の計算
1)Aで計算した残差の配列から残差を最小とするような軌道要素の変化分の解を求める。
2)オリジナルの軌道要素に変化分を加える。
3)改良された軌道要素をファイルに書き込む。(上書きではなく追記しています。)
4)変化分が大きければ@へ帰り繰り返す。
(添付ソフトでは、ループ5回 or Tの変化が10の-10乗以下でループを抜けています。)
B計算結果
上記手順で計算したときのディスプレイ表示結果の後半がresults31.txtです。 今回の計算例では残差の大きい観測値も見受けられますが、次に摂動を考慮した改良を行いますのでそのままとしています。
軌道計算(1) 地心太陽座標の計算のC言語ソース:
comet_v3.0.zip
注:windows上のCコンパイラ(VC等)でコンパイルください。
JPL/NASAのDE405のファイル(lnxp1600p2200.405)がほかに必要です。(ftp://ssd.jpl.nasa.gov/pub/eph/planets/Linux/de405/ 参照)
"jpl_eph.c"の217行目の天体暦ファイルのディレクトリ/ファイル名は環境に合わせて修整ください。
|