小惑星の軌道の数値積分プログラム

訳者注:このページのオリジナルはこちら

最終更新2004年5月3日

  • このプログラムの用途
  • このプログラムの使い方
  • 積分の C/C++ ソースコード
  • 惑星位置に JPL 推算暦を使う
  • このプログラムで使われる方法
  • 将来/検討している改良

    このプログラムの用途: ほとんど一般に使われる小惑星の軌道要素データセットは 小惑星センターの MPCORB と、Lowell 天文台の ASTORB ファイルです。 これら両方の要素は特別の標準元期で提供されます。その元期は、現在の日から100日以内です。

    これは通常、大きな問題にはなりませんが、これにより小惑星の "まさに現在の" 位置を2体問題で計算することになり、 100日までの摂動を無視することになります。メインベルト天体については、このことは通常 1秒角程度の誤差を与えます。 地球接近型天体で最接近をした後については、数分角から数度の誤差となります。また小惑星が、例えば50年前(古い乾板に それらを見出そうとして)にどこにいるかを知ろうとした時、不幸なことになるでしょう。現在の軌道要素は役に立ちません。

    一つの解決策は、要求に応じて軌道を数値積分するソフトを使うことです。私はこれを私のソフトに取り入れることを 検討し、結局使うことにしました・・・しかしその結果は速度の面と使い方の容易さで問題があり、大きな要素セットを 一度に積分する方が便利であり、そしてそのようにしました。

    私の商売上のソフトでは、Lowell データセットを使って、積分ルーチンの代りに、数十年の期間をカバーする 200日毎の 要素を生成しました。任意の日時の小惑星の軌道要素は、プログラムがどの要素セットがその日時に "最も近い" かを計算します。 これは良い解決策でしたが、莫大なファイル(バイナリーフォームを使っているにも関わらず)を必要とします。 現在登録されている小惑星では、CD に 200Mバイトの圧縮した小惑星軌道データを載せることを我慢しなければなりませんでした。 小惑星の数が増え続けるにつれて、私はタイムスパンを短くするか、何か別のデータを CD から削除しなければならないでしょう。 そして人々は CD が製造される時点で利用可能な小惑星データのみで我慢することになります。

    そこで私は、新しいオープンソースのフリーウェア・ユーティリティを提供することにします。 INTEGRATMPCORB.DATフォーマットのファイルを読み込み、新しい、ユーザーが明示した元期で作成します。

    ソフトの使い方: まず、このソフト(約152kバイト)をダウンロードします。 ディレクトリを作り、そこに解凍して、DOS シェルで INTEGRAT.EXE を走らせます。コマンドライン・パラメータ無しで これを行うと、以下のようなメッセージが表示されます:

    INTEGRAT takes as command-line arguments the name of an input file
    of the MPCORB.DAT type;  the name of the output file that is to be
    created;  and the epoch of that output file,  as either a JD or
    in YYYYMMDD form. For example,  either:
    
    integrat mpcorb.dat 2452600.mpc 2452600.5
    integrat mpcorb.dat 2452600.mpc 20021122
    
    would read in the 'mpcorb.dat' file,  and create a new file
    updated to the epoch JD 2452600.5 = 22 Nov 2002.  One can also
    specify the date as "today",  or "today+7" (i.e.,  this day next
    week) or "today-30" (a month ago). 

    こうして新しい 2452600.mpc を、MPCORB.DAT を使う多くのプラネタリウム・プログラム(私のも含めて)に 置き換えることで、そのプログラムは改訂された元期の要素で走ることになります。

    注意:この積分を行うには少し時間がかかります。計算にかかる残り時間を見積もっておくと良いでしょう。 最新の MPCORB を使い、(100日かそれ以下の)日付を進めると、500MHz マシンで約 15分かかります。 この時間はタイムスパンに比例し、思うに、数十年戻って積分すると、一晩は走らせる用意が必要でしょう。

    他の、例えば nea.dat といった、mpcorb ライクのファイルの 1つを積分すると、この問題はなくなります。 理由は単純で、このファイルは mpcorb.dat の軌道のごく一部の断片しか含んでいないからです。

    動作時間はインプットされるファイルの小惑星の数に比例します。数十個を残して全てを削除すると、 (例えば nea.dat で起こったように)長いタイムスパンでも、完全にリーズナブルな時間で行うことができます。 (このような場合、時間のほとんどは惑星位置の最初の計算の使われます。JPL 推算暦を使うと 良くスピードアップします)。

    C/C++ ソースコード: integrat.cppintegrat.mak が上の ZIP ファイルに含ままれています。 basic astronomical functions library と、 JPL ephemeris source code を使って、全てをコンパイルし、実行させる ことができます。

    惑星位置に JPL 推算暦を使う: デフォルトでは、INTEGRAT は惑星と月の位置に VSOP87法の省略版を使います。これは少しの誤差を与えますが、通常は 重要ではありません。これの最大の欠点は、軌道を長期間に渡って前方や後方に積分する場合で、プログラムをスタートさせてから、 タイムスパンをカバーする全ての位置を計算するのに長い沈黙となります。

    それらの理由のため、JPL 推算暦を使うようにするのが助けになります。 これを行うには、JPL DE ファイルを入手する必要があります。それらはインターネットからダウンロードするか、Guide 8.0 の 2枚目の CD-ROM からコピーしたり、Willmann-Bellから購入できます。これらのオプションの全ては それぞれのリンクに詳しくあります。

    DE ファイルを手に入れたら、INTEGRAT に '-f' スイッチでそれを使わせるようにします。例として:

    integrat mpcorb.dat 2452600.mpc 2452600.5 -fd:\unix.405 

    は、INTEGRAT に推算暦のソースとして d:\unix.405 を使うようにしています。

    このプログラムで使われる方法: ほとのど全てのケースで INTEGRAT は、固定ステップのルンゲ・クッタ・フェールベルグ法を使います。 固定ステップを使うことは、惑星位置を毎回計算するのに時間を費やさなければならないということがなく、 惑星の位置のテーブルを一度にセットすることができるという長所があります。 また、Encke の方法も使われています。これは、軌道を "直接" 積分する代わりに、実際の運動と 楕円軌道の2体運動との差分のみを積分する方法です。このことはプログラムの複雑さを大きなものにしましたが、 より大きなステップを扱えるようになりました。

    RKF を使うことは、大雑把な誤差の見積もりが毎ステップにあるということを意味します。誤差の見積もりが大きいと (通常、小惑星が惑星に近づいた場合)、プログラムはより小さいステップサイズで計算し、試みます。 (時々、小さいステップであるにも関わらず極端な誤差を導くことがあり、プログラムはまたより小さいステップで 計算を試みます。)最終的には、"問題のあるステップ" をカバーします。運がよければ、次のステップは不確かにならず、 プログラムは既に計算した惑星位置テーブルを使うことに戻ることができます。どのようなケースでも、"問題のあるステップ" が しばしば起こることがないように、許容誤差とステップサイズがセットされます。(それらが起こったとき、 それらは正確に扱われます。それらは "通常の" ステップほどには早くありません。)

    上記の RKF 法に加えて、プログラムがマルチステップ法(Gauss-Jackson 法のような)を持つと良いでしょう。 仮に4つ前までのステップからのデータを要求するような(プログラムが)実現できたと仮定すると、そのようなケースでは、 RKF法を 4つのステップを生成するのに使い、その後それを止め、Gauss-Jackson を走らせます。その時々に、Gauss-Jackson は 過大な誤差を生じます。そのようなケースでは、ステップは RKF が代りに行います。このようにすることで、Gauss-Jackson を ほとんどの時間に使うことができます(それによりかなり高速になります)。しかしながら、もし小惑星が惑星の近くを 通った時に誤差の蓄積があることが分かれば私達は精度を保つ方法に切り替えることができる、ということ知っておくことが 依然として安全です。

    変更: (2004年 5月 3日)先ごろ、MPC と そのミラーサイトにある mpcorb フォーマットで提供された NEA 天体に問題が生じました。特に、neatod.dat、これは 今日の元期の NEA 天体を提供することになっているものですが、それにいくつかの最近の天体が抜けていました。

    これを逃れるために、nea.dat をダウンロードした人々は、"今日" から100日前の元期でも、 今日の元期まで積分して、独自に neatod.dat を作成することができます。例として:

    integrat nea.dat neatod.dat 20040503 

    私はこの処理を単純化し、integrat の積分を一つのコマンドライン・パラメータである "today" にすると、 "nea.dat を今日の元期まで積分し、結果を neato.dat に保存する" という意味にしました。 これは、単純にこのように走らせます。

    integrat today

    すると、希望通りの結果が得られます。また、また、元期のユリウス日や年月日(YYYYMMDD)の代わりに "today" を使うことも できますし、 "today-7"(日)や "today-30" という表現にして、1週間前や 1月前の元期のファイルを得ることもできます。

    そして最後に、integrat は今回、出力ファイルのヘッダーに数行を追加して、ファイルがこれこれのスタートの 日付から現在の元期まで積分されたことを、そして integrat のバージョン日付を述べています。 たいしたことではありませんが、これは、どこかから入手した mpcorb フォーマットのファイルであるかや、 どれほどの量の処理をしたのかを忘れたときに、ヘッダーで確認できるようにしています。(私が、これを必要とする 唯一の人かもしれません。たぶん!)

    (2004年 2月 6日)Cristovao Jacques 氏の質問に対しての返答としてアップデートを登録しました。 新しく登録したバージョンは、小惑星の摂動を扱う処理に良いものです。以前は、integratmpcorb.dat にある 1番目、2番目、4番目の小惑星を、(1) Ceres、(2) Pallas、そして (4) Vesta と想定し、小惑星の軌道を積分する時に これらによる摂動を考慮していました。"通常の" mpcorb.dat ファイルを使っている場合はそれで正しいのですが、 MPC が提供する NEOs や他の天体のファイルでは、 Find_Orb's Monte Carlo ルーチンで作り出される mpcorb.dat-フォーマットファイルと同様ですが、 これらの 3天体を持っていませんintegrat は今回、正確な天体名を探し、ある小惑星がファイルのどこに あるかを基礎として、誤った結論を作りません。

    将来/検討している改良: このプログラムの主な欠点は、ユーザーフレンドリィでないということです。私はまず、他の事をする前に、この点に集中する でしょう。

    一度、ユーザーフレンドリィに、よくできた Windows インターフェイスにすれば、 マルチステップ積分法を追加することにします。これは、プログラムの操作やその精度に違いを与えるものではなく、 (思うに)かなり速くなることでしょう。

  • 訳者追加: integrat2.htm今年日本で観測に成功した小惑星による恒星食を、Guide 8.0 の build-in データと、integrat.exe で積分したデータとで食の経路を計算した結果。