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