JPL DE/LE の計算
JPL DE/LE のファイル構造 - 見上げれば、空 からのつづき。
時刻 t の正規化
ヘッダで読み取った(エフェメリス全体の)開始日時を start_epoch 、刻み日数を step とし、求めたい日時を jd とするとき
record_index = ((jd - start_epoch) / step).floor
と、開始日時からの日数をステップで割るとレコードのインデックスがでてきます。ただし、係数が記述されたテキストファイル/バイナリファイルは分割されていることもあるため、そこのところもおさえつつレコードを探してアクセスします。
# レコードの開始日時 start と終了日時の差が step なので # t は [0.0, 1.0) の値になる t = (jd - start) / step # 対象の天体の副区間数をたとえば 4 としたとき # [0.00, 0.25) は 0 つ目の副区間 # [0.25, 0.50) は 1 つ目の副区間 # ... とし、その区間において [0.0, 1.0) とする # 副区間数を n_subi とおくとき t -= (t * n_subi).floor # t を今度は [-1.0, 1.0) に正規化する t = 2.0 * t - 1.0
チェビシェフの多項式
T[0] = 1 T[1] = t T[2] = 2 * t^2 - 1 T[3] = 4 * t^3 - 3 * t : : T[n] = 2 * t * T[n-1] - T[n-2]
上記の漸化式で導かれる多項式があります。
pos = c[0] * T[0] + c[1] * T[1] + ... + c[n-1] * T[n-1]
で近似できるように係数 c を決定しています。この計算を行うことにより位置を求めることができます。単位は天体は km 、章動秤動に関しては rad です。
位置を時間で微分すると速度を求めることができます。時刻 t が関係するのはチェビシェフの多項式の方なので、多項式を微分すればよいのです。
T'[0] = 0 T'[1] = 1 T'[n] = 2 * (T[n-1] + t * T'[n-1]) - T'[n-2]
となります。ただし、 t の 1 の幅は step / n_subi / 2 日であるため、その日数における速度になってしまいます。そこで、 2 * n_subi / step をかけることで km/day ないし rad/day の値が得られます。
赤道座標、黄道座標に
日月惑星の位置の求め方を参照してください。
testpo
計算した値が正しいか確認するために、 JPL が用意した testpo.XXX と値をくらべてみます。
- testpo.XXX
- ascii/deXXX/testpo.XXX
- Linux/deXXX/testpo.XXX
# testpo.405 de# -- date -- -- jed -- t# c# x# -- coordinate --- EOT 405 2012.01.01 2455927.5 7 8 3 6.4557310425563
各行に、DE のバージョン、日付、ユリウス日、対象の天体番号、基準の天体番号、位置ないし速度のインデックス、値が記されています。これに書かれているユリウス日における天体の位置と速度を求め、それと記載値とを比較することで計算が正しいか、そしてどのくらいの精度で計算できているかがわかります。
天体番号は、
1 => 水星 2 => 金星 3 => 地球 4 => 火星 5 => 木星 6 => 土星 7 => 天王星 8 => 海王星 9 => 冥王星 10 => 月 11 => 太陽 12 => 太陽系重心 13 => 地球-月重心 14 => 章動 15 => 秤動
対象天体番号が 14, 15 のときの中心天体は 0 がセットされます。引きようがないので、単純に章動、秤動の値を採用します(もしくは 0 を太陽系重心とみなします)。
位置の単位は au で速度の単位は au/day です。ただし、章動と秤動においては rad, rad/day です。なお計算で求めた月は地心でしたが、ここでの月は原点中心です。
インデックスは 1 から x, y, z, vx, vy, vz です。章動、秤動も座標値ではありませんが、同様に処理します。ただし章動に関しては要素数が 2 のため、速度とあわせて 4 つしかありません。対応するインデックスは 1, 2, 4, 5 ではなく 1, 2, 3, 4 となります。
天文単位 AU および地球と月の質量比 EMRAT はヘッダの定数を用います。計算にあうように微調整された定数ですので、理科年表等から索くよりもこちらを用いた方がよいでしょう。
時間と座標
メモがてら。よくわかってないので間違いがあるはず。
JPL DE405 は TDB 相当の時刻系、 IRCS の座標系を採用。
国立天文台報 Volume 11 暦象年表の改訂について (PDF) がまとまっている。
時刻系
- 国際原子時 (TAI; International Atomic Time)
- 世界時 (UT1; Universal Time 1)
- 地球の自転速度を観測することによって決まる時刻系。 1 秒は平均太陽日の 1 / 86400 に要する時間に相当し、不規則
- 協定世界時 (UTC; Coordinated Universal Time)
- 地球時 (TT; Terrestrial Time)
- 地球のジオイド表面における座標時。 1 秒は SI 秒。以前は地球力学時(TDT; Terrestrial Dynamical Time)や力学時(TD; Dynamical Time)と呼ばれてたもので、今でもこちらで書かれることがある。 TDT の前身である暦表時との連続性をもたせるために 32.184 秒のずれをもつ
- TT = TAI + 32.184秒
- 太陽系力学時 (TDB; Barycentric Dynamical Time)
- 地球重心座標時 (TCG; Geocentric Coordinate Time)
- 地球重心を中心として運動する天体を考えるときの時刻系
- TT と比べると、1年あたり22ミリ秒ほど速く進む。周期項はない
- dTT / dTCG = 1 - L_G
- L_G = 6.969290134E-10
- 太陽系重心座標時 (TCB; Barycentric Coordinate Time)
- 太陽系重心を中心として運動する天体を考えるときの時刻系
- TT と比べると、1年あたり0.49秒ほど速く進む。また小さな周期項が存在する
- dTCG/dATCB = 1 - L_C
- L_C = 1.4808268674E-8
- TCB = ATCB + 周期項(1.6ミリ秒程度)
- ATCB = Average of TCB
まとめると
- TT は TAI と同じ速さを持ち、オフセットがあるだけ
- TDB は均すと TAI と同じ速さを持つ
- TCG は TAI より1年あたり22ミリ秒ほど速く進む
- TCB は均すと TAI より1年あたり0.49秒ほど速く進む
- TDB と TCB は最大振幅1.6ミリ秒ほどの周期項が存在する
- Barycentric な時刻系は周期項を持ち、 Coordinate な時刻系は TAI より速く進む
式にすると(たぶん)次のようになる
- 1977-01-01 00:00:00.0 を両軸の原点にとり、 TAI における時刻を t とする
- TAI = t
- TT = t + 32.184
- TDB = t + 32.184 + TDB_0 + 周期項
- TCG = t / (1-L_G) + 32.184
- TCB = t / (1-L_B) + 32.184 + 周期項
次に UTC や TT などの関係。
座標系
- 赤道座標系
- (略)
- 黄道座標系
- (略)
- 地球重力天文座標系 (GCRS; Geocentric Celestial Reference System)
- 地球重心まわりの運動を考えるとき。 TCG と対応。
- 太陽系重心天文座標系 (BCRS; Barycentric Celestial Reference System)
- 太陽系重心まわりの運動を考えるとき。 TCB と対応。
- 国際天文基準座標系 (ICRS; International Celestial Reference System)
- 太陽系重心を中心とし、超遠方の電波天体を用いて座標軸を定義したもの。
- IAU により勧告された
これにより、 ICRS という固定の座標系を用い、赤道や黄道自体が歳差運動をしているとみなす。
ICRS は今現在は赤道座標系の軸とほとんど差異がなく、原点の移動を行えばおおまかな用途にはそのまま用いることができる(ようだ)。正確に求める場合や、任意の時刻で計算する場合は座標変換が必要となる。変換方法はいくつかあるようだが、国立天文台では A New Precession Formula - Abstract - The Astronomical Journal - IOPscience が用いられている。