|
| |
※「セキュリティ保護のため...」というメッセージが出る方・日本語が入力できない方へ
天体計算の基礎(8)
第８章 実験的な近似式による惑星位置の計算
【２体問題の修正】
惑星の位置を２体問題の発展として捉えた場合、天体の位置を「乱す」他の惑星や
衛星・彗星などから受ける引力のことを摂動と呼びます。この摂動の主なものを拾い
上げ、実測値に基づいて天体の位置を推算する計算式を多くの人が開発しています。
なお、ここまで、似たような言葉が幾つか出てきたので少し整理しておきましょう。
歳差 precession 春分点が黄道に沿って毎年前進して行く現象。地球は球体ではない
ので、太陽や月から受ける引力の中心が地球の自転軸上にない。このため自転軸
（＝赤道面に垂直な直線）は黄道の極の回りを約25,800年を周期に１回転する。
結果、春分点も同じ周期で黄道上を１回転することになる。これは70年で1度
くらいだから結構大きな変動量である。
これは太陽と月による歳差（日月歳差）であるが、当然他の惑星からの歳差も
ある。日月歳差に惑星歳差を加えたものを一般歳差という。
章動 nutation 同じく太陽や月の引力によって自転軸の方向が周期的に揺れる現象
をいう。別の言い方をすれば、太陽や月の引力による春分点の変動の内、最も
大きな周期の変動を歳差と呼び、それ以下の変動を章動というのである！！
章動の中で最も大きな成分は周期が約18.6年、振幅10秒程度のものである。
摂動 perturbation 惑星の軌道が太陽との２体問題と考えたときからずれる現象。
これは惑星が他の惑星からの引力も受けることにより生じるものである。
衛星の軌道が母惑星との２体問題と考えたときからずれる現象も同様に摂動
と呼ばれる。
【摂動について】
摂動については少しだけ触れるに留める。
３次元空間の２体問題を考えた場合、引力の場のポテンシャルを座標の各成分で
偏微分すると、その場にある質点にかかる力が求まることを想起しておく。
今、空間上に２つの質点と＋摂動を与える質点１個があった場合、１つの質点の
位置ベクトルを r とすると、引力の法則により、dr/dt = -G(M+m)r/abs(r)^3 + F
と書ける。最後の + F が摂動であり、余分な１点と従来の２点との間に働く重力に
基づく力である。（d は本当は偏微分の記号「∂」だが旧ＪＩＳのマシンでは表示さ
れないので d で代用する） ↑ミエルカナァ? ソノ 旧JISマシン デ ウッタノデ...
この時、場所に依存するある関数 R が存在して、ポテンシャルと同様に座標による
偏微分によって、この摂動力 F を得る事が出来る。
つまり、
dr/dt + G(M+m)r/abs(r)^3 = dR/dr
であるような R が存在する。この引力のポテンシャルに似た性質を持つ、摂動力を
与える関数を、摂動関数(disturbing function)と呼ぶ。
この関数は、摂動を与える点の位置ベクトルを j、質量を N とすると、摂動関数は
G・N(1/abs(r-j) - j・r/abs(j)^3)
の形で表現できる。ここに j・r は内積である。
この式の第１項は摂動を与える点（他の惑星）が問題とする点（惑星）に与える引力
を表しているので直接摂動(direct perturbation)という。第２項は摂動を与える点
（別の惑星）が中心とする点（太陽）に与える引力の影響であり間接摂動(indirect
perturbation)という。
この式は単純な形はしているが、このままでは実際の処理には困るので、実際の摂動
を計算するにあたっては、この関数を多項式や三角関数のべき級数の形に変形する
ことが必要である。この過程は作業としては単純な作業だが、作業量としては膨大な
ものになる。
長沢工氏によれば、それは「２次項まで展開するのは相当に粘り強い努力を長いこと
続けてやっとできるといった程度であり、３次項の展開となると通常の努力では
とうてい成し遂げることができないほど大変な量の計算をしなければならない」と
いうことである。
次項で出所を説明する、本プログラムの中核となる近似式は、恐らくこの展開を実行
した結果、得られた物なのではないかと想像する。そういうものを公開して頂けた
というのは、大変な感激である。
【本プログラムの近似計算式について】
このプログラムでとりあげている各惑星の位置の概算式は、
暦計算研究会編「新こよみ便利帳」恒星社厚生閣
に掲載されていたものを入力したものです。著作者である暦計算研究会殿からは
当フォーラムＦＦＯＲＴＵＮＥにＵＰするということと、商業的な利用・漁業への
転用を禁止するとうたう、と述べて、ソースを含めた公開の許可を得ています。
【この近似計算式の精度】
この計算式の精度は次の「程度はある」と推定される。範囲をせばめると、あと
１桁は精度がいいのではないかとも想像するのだが、精度の高い位置表をまだ入手
していないので、何とも言えない。（値段が張りそうなので、お金に余裕が出たら
入手に向けて動き始めたい。^^;; ）
太陽 1850-2050で0.01度以内
月 1850-2050で0.03度以内
水星 1850-2050で0.02度以内
金星 1850-2050で0.01度以内
火星 1850-2050で0.01度以内
木星 1850-2050で0.07度以内 1900-2050では0.03度以内
土星 1920-2050で0.03度以内
天王星 1945-2050で0.02度以内
海王星 1935-2050で0.02度以内
冥王星 1940-2050で0.07度以内 1945-2005では0.02度以内
【計算式の補正について】
今回のプログラムでは、1880〜1895年付近を処理のターゲットにした為、上記の計算
式では木星〜冥王星の位置の計算に困ることとなった。そこで、これは各々について、
計算式の補正・或は別の方法による計算を試みることにした。
なお、今回は超手抜きで、地心黄経のレベルでの補正を行った。もし日心黄経で処理
していたら、もう少し精度の良いものが出来たかも知れないが、（後で気付いたので、
少々疲れたこともあり実行していない）
（１）木星
1850〜1935について、誤差を調査したところ、概ね指数関数的な
誤差の拡大が見られたので、まず指数関数により補正し、その後
必要なら他の補正を試みることにした。この補正に必要な定数は
Microsoft-Excelで少しずつ数値を変えながらシミュレーション
して、決定した。結局指数関数の補正だけで、間にあった。
（２）土星
最初1850〜1920でやろうとしたが、それではかなり沢山補正項を
加える必要がありそうだったので見送り、1880〜1895について、
処理することにした。誤差を調査したところ、概ね直線的に誤差
の拡大が見られたので、まずフィットする直線の成分を決定して
補正すると、今度は周期的な変動が残った。そこでこれを三角関数
で補正したが、まだ揺れがあったので、更にもう１個三角関数を
加えた。このシミュレーションの過程はやはりMicrosoft-Excelで
処理した。
（３）海王星
1850〜1935について、誤差を調査したところ、指数関数的な伸びと
わりと大きな周期的な揺れが見られたので、まず指数関数の成分を
大雑把に除去し、その後は土星と同様２つの三角関数で補正する
ことにした。処理はやはりMicrosoft-Excelである。
（４）冥王星
冥王星はターゲットとする期間での計算結果が余りにも実際の推測
位置と掛け離れていたので、この近似式による計算を諦め、軌道
要素による２体問題で解いてみることにした。今回はこの計算は
殆ど目算^^;;で行った。軌道要素の内、最も影響の強いのは元期平均
近点離角と近日点引数なので、これ以外の要素を 長沢工「天体の位置
計算・増補版」に転載されていた軌道要素の計算式で計算した上で
この２つの数値の推測値を少しずつ手作業で変更しながら、うまく
フィットする値を見つけ出した。（よく見つかったものだ ^^;;）
（５）天王星
天王星もターゲット期間での計算結果が、かなり実際の推測位置と
離れていたので、冥王星と同様に軌道要素からの２体問題で処理
しようとした。しかし今度はなかなか手作業ではうまい値が見つか
らなかった為、ターゲット期間の代表点数カ所を取った最小自乗法
により最も誤差の少ない軌道要素の値を見つけることにした。この
場合も元期平均近点離角と近日点引数以外は同様の固定推測値である。
（或はこの推測値が悪かったのかも知れない）
これを実行するプログラムは、方程式を数値的に解く処理と同様に
両側から解の値を挟み込んで行くやり方で最も質の良いものを自動
的に求めるようにした。このプログラムは恐らく組み方が悪かった
のであろうが、２時間程動き続けて、やっと答えを出してくれた。
しかし、こうして求めた値を使って実際に軌道計算を行ってみると
かなり近い値は出るようになったが、まだまだ誤差が含まれていた。
そこで、非常に不本意ながら、これに更に地心黄経に対する補正を
行うことにし、海王星などと同様にsin関数で補正した。この部分の
処理は、またもやMicrosoft Excelである。
【誤差吸収用のｓｉｎ関数の決め方】
これは次の手順で決めた。こういう方法が正しい方法であるかどうかは
定かではない。(1)揺れ幅を見て、振幅（sinに掛ける定数）を決める。
(2)揺れの底から天までの時間及び、次の底までの時間等を総合的に判断
して周波数（sin( a・(x-t) ) の a ）を決める。(3)揺れの底から天に昇
って行く中点付近に左記の t の値を求め、この a(x-t) を展開し ax - at
として、 -at にて位相を決める。今回はこういった数値は目分量で決め
ている。（本当はちゃんと出し方があるのであろうが。。。）
※ 今回比較に使用した実際の位置の推測値としては、★秋津★氏のＳtagazer
の表示値（これは数値積分に基づく値である）と長沢工「天体の位置計算・
増補版」に転載されていた Seidelmann-Doggett-DeLuccia の方法による
計算結果を使用した。
実際の作業にあたっては、後者だとプログラム内で計算できる為、だいたい
の評価をするのに後者を使用し、最終的な値の合わせ付けを行う段階では
Ｓtargazerの値を尊重した。Seidelmann-Doggett-DeLuccia の方法を使用
する部分については今回の公開版からは外している。
（アメリカまで許可を求める手紙を出す気力が沸かない。。。。。）