JPL DE/LE シリーズ
利用し方など、基本的なことは過去の記事を。
いくつものバージョンを試したことで気付いた差異とかをメモ。
DE/LE の種類
start final step ksize ncoeff consts ------------------------------------------------------------- 102 1206160.5 2817872.5 64.0 1546 773 152 200 2305424.5 2513392.5 32.0 1652 826 200 202 2414992.5 2469808.5 32.0 1652 826 32 403 2305200.5 2524400.5 32.0 2036 1018 144 403 2305200.5 2524400.5 32.0 2036 1018 144 405 2305424.5 2525008.5 32.0 2036 1018 156 406 625360.5 2816912.5 64.0 1456 728 156 410 2415056.5 2458832.5 32.0 2036 1018 186 413 2414992.5 2469872.5 32.0 2036 1018 235 414 2305424.5 2525008.5 32.0 2036 1018 259 418 2414992.5 2470192.5 32.0 2036 1018 228 421 2414992.5 2524624.5 32.0 2036 1018 228 422 625648.5 2816816.5 32.0 2036 1018 222 423 2378480.5 2524624.5 32.0 2036 1018 222 424 625296.5 2816816.5 32.0 2036 1018 222 430 2287184.5 2688976.5 32.0 2036 1018 572 (229) 431 -3100015.5 8000016.5 32.0 2036 1018 572 (229)
この start と final は JPL/DE としての対応値であって、実際に頒布されているアスキー／バイナリファイルはこれらよりも範囲の小さいものであることもある。ので、チェック大事。
DE430/431 は定数が572個記述されたヘッダファイルの他、229個のものも同時に頒布されている。定数が400個までしか対応していないソフトウェアの場合は229個版を用いる必要がある。頒布されてるバイナリファイルは572個版。
拡張されたバイナリ形式
DE430/431のヘッダファイルには400個以上の定数が記述されている。また、DE431では12個の ipt テーブル。1個の lpt （秤動）テーブルの他、 rpt テーブルと tpt テーブルがそれぞれ記述されている。（拡張されたといってもDE431の場合は rpt と tpt の係数の数が0個なので、係数としては含まれていない）
rpt テーブルは月のオイラー角に関するデータ（自転軸？ lunar euler angle rates と書かれている）。 tpt テーブルは TT-TDB 。
バイナリにおける各データは lpt のあとのパディング領域に記述され、400を超える場合の定数名、 rpt 、 tpt の順となり、残りはパディングされる。
2ブロック目の定数値はそのまま連続して記述される。
エポックの位置
テスト用かは知らないが、定数表の中にはあるエポックについて、ある地点からの位置が含まれている。が、微妙にバージョン間に差異がある。
まず、 X1, Y1, Z1, XD1, YD1, ZD1 は水星における位置（au）と速度（au/day）。以下 2, B, 4, 5, 6, 7, 8, 9 は金星、地球月重心、火星、……。 M は地心の月（CENTER に依らず）。
ここまではいいが、 S と C が微妙。
JDEPOC CTR XS YS ZS XC YC ZC ------------------------------------------------------------------------------- 102 1240400.5 0 0.008349 0.003809 0.001338 200 2440400.5 11 0.002954 0.005187 0.001658 202 403 2440400.5 11 0.004503 0.000767 0.000266 405 2440400.5 0 0.004503 0.000767 0.000266 406 2440400.5 0 0.004503 0.000767 0.000266 410 2447952.5 11 -0.000342 0.000235 0.000066 413 2447952.5 11 0.0 0.0 0.0 414 2440400.5 11 0.0 0.0 0.0 418 2440400.5 11 0.004503 0.000767 0.000266 0.0 0.0 0.0 421 2440400.5 0 0.004503 0.000767 0.000266 0.0 0.0 0.0 422 2440400.5 11 -0.004503 -0.000767 -0.000266 423 2440400.5 11 -0.004503 -0.000767 -0.000266 424 2440400.5 11 -0.004503 -0.000767 -0.000266 430 2440400.5 0 0.004503 0.000767 0.000266 0.0 0.0 0.0 431 2440400.5 0 0.004503 0.000767 0.000266 0.0 0.0 0.0
CENTER の0は原点、11は太陽を意味している。
基本的に、 CENTER に関わらず S は原点から太陽の位置を示しており、 C は原点を示している。
が、DE413/414は太陽中心として太陽の位置をみている（つまりゼロ）。一方で、DE422/423/424の C は太陽中心から原点をみている（つまり正負逆）。
なぜ微妙に仕様を変えるのか。解せぬ。
あと、DE200の S の値が何を表すか不明。
よくよくみてると、DE410/418は CENTER=11 にも関わらず、原点中心の値ですね。
まとめると
CTR 12B456789 M S C
---------------------------------------------------------------
102 0 中心から惑星 地心から月 原点から太陽
200 11 中心から惑星 地心から月 不明
202
403 11 中心から惑星 地心から月 原点から太陽
405 0 中心から惑星 地心から月 原点から太陽
406 0 中心から惑星 地心から月 原点から太陽
410 11 原点から惑星 地心から月 原点から太陽
413 11 中心から惑星 地心から月 中心から太陽
414 11 中心から惑星 地心から月 中心から太陽
418 11 原点から惑星 地心から月 原点から太陽 原点から原点
421 0 中心から惑星 地心から月 原点から太陽 原点から原点
422 11 中心から惑星 地心から月 中心から原点
423 11 中心から惑星 地心から月 中心から原点
424 11 中心から惑星 地心から月 中心から原点
430 0 中心から惑星 地心から月 原点から太陽 原点から原点
431 0 中心から惑星 地心から月 原点から太陽 原点から原点中心とは CENTER = 0, 11 の天体から
章動
章動は（比較的単純な）数式で求められる歳差とともに、 GCRS から赤道・黄道座標へと変換する際に必要となります。
値の導出のもととなる惑星の運行データに JPL/DE で計算した値を使っているのでしょう。各 DE でごくごく微量の差がありますが、最新の DE430/431 に至るまで IAU 1980 章動モデルと思われる章動値が含まれています（ DE102/406 の2種のバージョンには入っていない）。
章動モデルには IAU 2000A とそれを少し補正した IAU 2006A というモデルがあり、これらを使うには別途計算する必要があります。なお、国立天文台暦計算室では IAU 2006A 章動モデルを用いているそうです。
Fortran プログラム
使い方は userguide.txt に載っているが。
asc2eph.f
アスキーファイルをバイナリ化するもの。バイナリファイルの配布が行われていないDE424などで使用。
- 74行目か75行目どちらかのコメントを削除（一般的には74行目）
- $ fort77 -o asc2eph asc2eph.f でコンパイル
- $ cat header.424 ascp1600.424 ascp1700.424 | ./asc2eph
- JPLEPH というファイルが生成される（ただし JPLEPH が既に存在しているとエラー）
testeph.f
計算値のチェックをするもの。
- 220行目、285行目、357行目の NRECL に値を代入する（とりあえず 4 ？）
- 386行目の KSIZE に値を代入する（目的のバージョンにあわせて）
- 908行目から910行目までのどれかのコメントを削除（上記の KSIZE の設定が適切であれば910行目で構わない）
- 97行目など、いくつか DATA 行のあとに INTEGER などの宣言文がきているとコンパイラに怒られる（そうでないコンパイラもあるのかな？）ので、適宜、行を入れ替える
- $ fort77 -o testeph testeph.f
- あらかじめ、目的のバージョンのバイナリファイルを JPLEPH という名前で用意する
- $ cat testpo.405 | ./testeph
testpo ファイルについて
testpo には、秤動の3つ目の角 psi が、時間経過とともに精度が落ちるという理由で補正されて記述されています。そのため、 JPL DE/LE で計算した値との差を意図的に小さくしています。
上記の testeph.f では角度のみ補正し、角速度の補正が入っていないため、DE414の testpo.414 で許容値以上の差となり、その値が出力されます。
しかし、正しく計算できたかを確認するためのファイルに対して、 testpo の値を弄っているのは……
