Copyright 1989- Jun Makino
tar xvzf pccode.tgzカレントディレクトリに直接ファイルを展開するので、あらかじめこれ用のディ レクトリを作って移動しておくこと。
make pccodeケプラー問題(1体問題)用のプログラムは
make testbed3
Enter order, hermit_flag, nrot, dtmax:とでる。order は積分公式の次数。多分16くらいまでは動く。hermit_flag は エルミート型か古典的なアーセス型かのフラグ。1ならエルミート。 nrot はケプラー問題を何回転解くかの回数。dtmax は可変時間刻みの上限。
これは(多分)軌道長半径1、中心天体の質量1、重力定数1のケプラー問題を解 く。離心率だけが初期条件のパラメータでこれはあとでいれる。これらを入力すると
Enter dttype(1-aarseth, 2-abs, lmode, 3-abs, gmode
4-rel, lmode,5-rel, gmode
とでる。ケプラー問題の場合は、1 のアーセス型でもまあ使えるのでこれを使
う。入力は実際には dttype, tol1, tol2, exfact の4つが必要で、aarseth
step の時には tol2 は使われないが指定は必要。例えば
1 0.01 0 1.2とかいれる。これは、 aarseth step のいわゆる eta (これの平方根がタイム ステップ条件の係数になる)が 0.01、タイムステップを大きくする時の上限係 数が 1.2 である。 次に
Enter eccentricity:とでるので、計算したい軌道の離心率をいれる。 最後に Enter epsinit: とでる。これは、このコードで実際されているアルゴリズムは可変次数なので、 初期には低次の公式から出発する。そのため、上の時間刻み条件とは別に、十 分小さいステップにする必要がある。これは、単純に
dt = epsinit/|加速度の導関数の絶対値|としているので、それで大丈夫な程度に小さくする。0.00001 とかにしておけ ば多分大丈夫。以下、実行結果の例
% testbed3
Enter order, hermit_flag, nrot, dtmax:
8 1 10 0.01
Order, Hermite, Nrot, dt : 8 1 10 0.10000E-01
Enter dttype(1-aarseth, 2-abs, lmode, 3-abs, gmode
4-rel, lmode,5-rel, gmode
1 0.01 0 1.2
dt expand limit = 1.200000
Aarseth dt, eta= 0.1000000E-01
Enter eccentricity:
0.9
eccentricity= 0.90000000000000002
Enter epsinit:
0.00001
xp 0 1.0000048595286 0.23535842029826E-05 0.0000000000000
vp 0 0.89999460053377 0.43588989434771 0.0000000000000
steps, tnow : 50 0.39707100043171E-01
x 0 1.0349662175239 0.17303607260498E-01 0.0000000000000
v 0 0.86164932739810 0.43556932417576 0.0000000000000
E, Eerr, L, Lerr : -0.5000000 2.7755576E-16 0.4358899 -7.2164497E-16
steps, tnow : 100 0.53716634517873
x 0 1.3712808368650 0.22766757989514 0.0000000000000
v 0 0.52425479684155 0.40491028553187 0.0000000000000
(途中一杯でるが省略)
E, Eerr, L, Lerr : -0.5000000 -4.2428283E-12 0.4358899 1.8818280E-14
steps, tnow : 7150 62.583964355194
x 0 0.74028398469568 -0.10639616606190 0.0000000000000
v 0 1.2263706048950 0.41255622185008 0.0000000000000
E, Eerr, L, Lerr : -0.5000000 -5.4761751E-12 0.4358899 -5.7509553E-14
steps, tnow : 7175 62.833964355194
E, Eerr, L, Lerr : -0.5000000 -5.4751759E-12 0.4358899 -5.8508753E-14
EerrMAX, LerrMAX = 5.5981886E-12 8.2489571E-14
対話的に入力する代わりに
8 1 10 0.1 1 0.01 0 1.2 0.9 0.0001と書いてあるファイル(これを testinという名前でセーブしたとして)を作っ ておけば
# testbed3 < testinで実行できる。位置、速度は現在のところ上のような感じに画面にでる。そのへんは適宜ソースファイルをいじるか、この出力を加工すること。もっと違う初期条件にするならそれもファイルから読むようにするとか改造すること。