日月歳差の影響を計算する: セッピーナの趣味の天文計算

2013年8月20日 (火)

日月歳差の影響を計算する

天体の赤経・赤緯に対する(日月)歳差の影響を計算する方法です。

一般歳差の場合も座標の回転を三回するという点は異ならないのでExcelのシートは流用可能です (^^;;

--------

これまでの経緯（計算方法だけ必要な方は読み飛ばしてください）

太陽の自転軸を表すパラメータB0は「太陽の自転軸の方向」と「地心から見た太陽の方向」のなす角から得られます。

「太陽の自転軸の方向」と「地心から見た太陽の方向」は国立天文台の暦から得られるのですが「太陽の自転軸の方向」は有効数字が不足しており必要とする精度でB0を計算することができません。

NASA JPLの暦から「太陽の自転軸の方向」の一桁精度の高い値が得られますのでこれを使って計算することにします。

ただ国立天文台の暦から得られる「地心から見た太陽の方向」は真春分点に対するものであり一方NASA JPLの暦から得られる「太陽の自転軸の方向」はJ2000.0に対するものなのでこれらを使って計算してしまうと正しい結果がえられません。

そこでNASA JPLの「太陽の自転軸の方向」に歳差の影響を加味すると正しい結果が得られました。

　　「太陽の自転軸のパラメーターP,B0,L0を「計算で」求める (5)」　編集

ただこの記事では具体的な計算方法を示さなかったので今回はExcelによる具体的な方法を記します。

注

NASA JPLの暦からは「太陽の自転軸の方向」も「地心から見た太陽の方向」もJ2000.0に対するものが得られますのでそれらを使えばいいだけの話ですし、もっと言えばB0の値は国立天文台の暦にあります。
この記事は歳差の影響をどう計算するかの例として読んでいただければと思います。

-------

今回はまず日月歳差のみを考慮した場合です。

これは比較的簡単で

0.赤経・赤緯を赤道座標（方向余弦）に変換する
1.赤道座標を黄道座標に変換する
　（x軸=春分点方向を中心に地球の自転軸の分だけ回転する）
2.z軸（黄道面に垂直な方向）を中心に歳差の分だけ座標を回転する
3.黄道座標を赤道座標に変換する
　（x軸=春分点方向を中心に地球の自転軸の分だけ回転する。1.とは逆向き）
4.赤道座標（方向余弦）を赤経・赤緯に変換する。

という操作で歳差の影響を計算することができます。

　　「座標の回転の記法」

に書いた表し方だと

1. R1(ε)　　ε: 赤道面の黄道に対する傾斜角
2. R3(-P）　　P: 歳差による回転角（春分点は天頂から見て時計回りに回転します）
3. R1(-ε)

の三つの操作を行えばいいということで三つ合わせると

　　R1(-ε)R3(-P)R1(ε) = R131(ε,-P,-ε)

の操作ということになります。

それから P≒0 であることを使って式を簡略化できるのですが、ここではそのままにしておきます。いくら計算が面倒でも計算するのは私じゃなくてExcelですから (^^:;
簡略化した式は極に近い天体で誤差が大きくなるという問題もあります。

------

まず R131(ε,-P,-ε) を求めるまでの計算式です。
R1(ε)にR3(-P）をかけてR13(ε,-P）を求めさらにそれにR1(-ε)をかけてR131(ε,-P,-ε) を求めます。

=mmult(.....) というのがたくさんありますがこれはマトリックスの（あるいはベクトル）の積を求めるための関数です。入力方法にコツ(?)があります。

例えばA20のセルのR13(ε,-P)のところですが....

ドラッグしてB20からD22の9個のセルを選択します。
そして「=MMULT(B14:D16,B17:D19)」を入力してCtrl+Shift+Enterを押下します。
「=MMULT(B14:D16,B17:D19)」は「=MMULT(」まで入力しB14:D16を選択し「,」、B17:D19を選択、「)」でもいいです。

次に上で求めたR131(ε,-P,-ε)を使って赤経・赤緯の変換（補正？）を行います。
赤経・赤緯から方向余弦を求めこれにR131(ε,-P,-ε)をかけて歳差の影響を考慮した方向余弦を求めます。そしてそれから赤経・赤緯を求めます。