日月惑星の位置の求め方 - 見上げれば、空

いくつかあります。

理科年表の値を用いて任意時刻の値を補間
- 視赤経・視赤緯
国立天文台暦計算室暦象年表の値を用いて任意時刻の値を補間（2009〜2014年；執筆時現在）
- 視赤経・視赤緯、視黄経・視黄緯、日心黄経・日心黄緯
海上保安庁海洋情報部天測暦付録の計算式を用いて任意時刻の値を計算（2008〜2013年；執筆時現在）
- 視赤経・視赤緯
JPL DE/LE を用いて任意時刻の値を計算（JPL DE405 の場合1600〜2200年）
- ICRS 直交座標
The Astronomical Almanac の計算式を用いて任意時刻の値を計算（2001〜2013年；執筆時現在）
- 視赤経・視赤緯；月のみ
The Astronomical Almanac の軌道要素を用いて任意時刻の値を計算（2008〜2012年；執筆時現在）
- 瞬時の平均黄道と春分点における日心の真位置；地球月重心を含め惑星のみ
JPL's HORIZONS の値を用いて任意時刻の値を補間（1600〜2200年）
- オプションにより様々な形式

このうち暦象年表と HORIZONS のウェブ版は非常に細かく時間の指定ができるため、秒レベルであれば値を取得し、補間計算を行う必要はありません。

暦象年表の値

で、暦象年表です。

データブックとして理科年表を使用している人も多いと思います。理科年表のように太陽・月・各惑星の位置を調べたいときに、便利なのが国立天文台の暦計算室サイト内にある暦象年表のコーナー。

秒単位で指定でき、また最大400サイクルまで出力できるので大量に使用する際に重宝します。十分に細かいステップでデータが得られるので、一次補間や三次補間など単純な補間法で秒以下の時刻についても近い値が得られます。ただし世界時・日本標準時は地球時からのオフセットが整数秒に固定されているため、地球時で取得し適宜補間するとよいでしょう。

ですが、ここで得られるデータは今現在2009年から2014年までに限られています。

そこで暦象年表サービスで出力される値に合致するように太陽・月・惑星の位置について計算してみます。同じ方法を用いることで、任意の時刻における値も得ることができます。

暦象年表で使われている理論

2003年度および2009年度に暦象年表は改訂が行われています。2003年度のものはトピックス「暦の改訂について」、2009年度のものは「暦象年表の改訂について」で簡単に、国立天文台報 Volume 11 p.57-67 の「暦象年表の改訂について」[PDF]では詳細に説明されています。

おおまかに抜粋すると2003年度から JPL DE405/LE405 の使用、2009年度から IAU 2006 の歳差章動理論の採用が挙げられます。

国立天文台報の「暦象年表の改訂について」[PDF]に沿って説明していきます。数式などはこの文書に載っているため併読しているものとして、式などを省きます。

JPL DE405/LE405

JPL DE405/LE405 の使い方については過去の記事を参照してください。

時刻系は JPL DE405/LE405 における t_eph と同一とみなしてもよい TDB を使います。

座標変換

JPL DE405/LE405 で得られる座標は ICRS 座標系です。これは原点が太陽系重心で、ある方向を軸にしたものです。地球の自転軸によって定まる赤道座標系や地球の公転軌道によって定まる黄道座標系とは異なり、太陽系においては常に固定されている座標系といえます。

また、 GCRS は原点が地球重心が原点の座標系で、 ICRS とは相対論的な差がありますが、ここでは単純に原点を移動させた座標系として考えます。日心の座標系においても同様です。

$X_{\mathrm{equat}} = R_1(-\epsilon_A - \Delta \epsilon) R_3(-\bar{\psi} - \Delta \psi) R_1(\bar{\phi}) R_3(\bar{\gamma}) X_{\mathrm{GCRS}}$
$X_{\mathrm{eclip}} = \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; R_3(-\bar{\psi} - \Delta \psi) R_1(\bar{\phi}) R_3(\bar{\gamma}) X_{\mathrm{GCRS}}$

$X = \left( \begin{array}{c} x \\ y \\ z \end{array} \right)$
$R_1(\theta) = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos \theta & \sin \theta \\ 0 & -\sin \theta & \cos \theta \end{array} \right)$
$R_3(\theta) = \left( \begin{array}{ccc} \cos \theta & \sin \theta & 0 \\ -\sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{array} \right)$

歳差における角度 $\bar{\gamma},\; \bar{\phi},\; \bar{\psi},\; \epsilon_A$ と章動における角度 $\Delta \psi,\; \Delta \epsilon$ を用いて座標軸を回転させて座標系を変換します。

右3つの R₃R₁R₃ の回転で黄道座標系に、更に左の R₁ の回転で赤道座標系に変換されます。章動の $\Delta \psi$ があれば真春分点、なければ平均春分点を x 原点にします。また $\Delta \epsilon$ があれば真赤道、なければ平均赤道を xy 平面にします。

地球の自転軸は約25,800年かけゆっくりと回転する歳差と、短い周期で細かく振動する章動があります。この章動項を考えない場合の赤道を平均赤道と言い、平均赤道と黄道とで定まる春分点を平均春分点と言います。章動を考えた場合は真赤道、真春分点です。

歳差・章動

歳差に関しては J2000.0 = 2000/01/01 12:00:00 TT からの経過ユリウス世紀である T の5次式が載っていますので、それを用います。ここでの T は TDB ではなく TT です。

章動は IAU 2000A の MHB2000 章動理論を修正したものを用います。その章動の計算は係数も多く、面倒です。実装として IAU SOFA の nut00a.c を紹介しておきます。もっとも、他の関数も必要ですが。プログラム群すべてのダウンロードページもあります。

JPL DE405/LE405 にも章動の項があるため計算で取得できますが、これは IAU 1980 の章動理論によるもののようで若干値が異なります。

地心座標

太陽・月・惑星の地心座標を求めます。各用語は「こよみ用語解説」に載っています。

視赤経・視赤緯、視黄経・視黄緯

地心座標なので、求める時刻 t₂ における地心 r₂(t₂) が原点にくるように平行移動させます。

「光行差を含みます」とあるので、まず対象の天体が光を発した時刻 t₁ を求めます。対象天体および重力源天体Jは時刻 t₁ の位置を用いるため、直接は求まりません。何度か繰り返し漸近させていきます。

次に光の曲がりの量 $\hat{\mathbf{d}}_{12}$ を求め、補正し $\hat{\mathbf{d}}'_{12}$ を出します。書かれていませんが（あるいは誤植ですが）、ここで $\mathbf{v}_2 = \frac{\mathbf{v}_{\mathrm{earth}}(t_2)}{c}$ であり、光速度 c で割られています。 v₂ も同様です。 JPL DE405/LE405 で求まる速度は1日あたりの変化量なので、86,400で割って秒速にしておくことも必要です。

視赤経・視赤緯、視黄経・視黄緯は「瞬時の真春分点と真赤道にもとづく座標系における値」とあるので章動を含めた回転処理を行います。 $\hat{\mathbf{d}}'_{12}$ は地心から見て光がやってきた方角なので、すなわち天体の視位置です。 $\hat{\mathbf{d}}'_{12}$ から経度・緯度に変換すると視赤経・視赤緯、視黄経・視黄緯です。

距離

距離は「地心距離は天文単位（au，1au = 149,597,870.7 km）で表した地球からの真距離です」とあるので、時刻 t₂ における真距離 |r₁(t₂) - r₂(t₂)| を求めます。

ただ、見ているとどうも IAU で定義された天文単位ではなく、 JPL DE405/LE405 の計算内で用いている 149,597,870.691 km で割っているように見えます。おそらくバグでしょうね。

視差・視半径

視差・視半径は書かれているとおりです。天体の赤道半径として使用している定数はおそらく水星 2439.7 km 、金星 6051.8 km 、地球 6378.137 km 、月 1737.4 km 、火星 3396.19 km 、木星 71492.0 km 、土星 60268.0 km だと思われます。

等級

面倒なので出してません。

日心座標

惑星の日心座標も求めてみます。

日心黄経・日心黄緯は「光行差は含みません」とあるので、真位置を用います。

「瞬時の平均春分点と黄道にもとづく座標系における値」とあるので、章動を含めず歳差のみの回転を施します。

動径は「太陽からの真距離です」とあるので同様に真距離を求めます。

軌道経度、日変化

軌道経度は $\Omega + \omega$ だと思いますが、ぴったり値が一致するような計算ができませんでした。春分点や昇交点がいつにおけるものなのかがポイントだと思うのですが……。

日変化もどこから求めた値なのかわかりませんでした。

準惑星・小惑星

JPL DE405/LE405 には含まれていません。

計算可能なエフェメリスはあるのでしょうか……？

グリニッジ恒星時・均時差

未稿。

プログラム

整形してそのうち出すかも。