これはFortranによる 3次元3体問題のプログラムです。
Fortran で計算し、gnuplot で描画（アニメーション）させることができます。


使い方


f90 プログラムソース
gnuplot 用テキストファイル
  1. どこかに新しいフォルダを作る。
  2. そのフォルダ内に、上の f90 ファイルを拡張子を f90 として保存する。
    同じフォルダに、上のテキストファイルを 3objs_move.txt として保存する。
    Windows の場合、そのフォルダに gnuplot も入れておく。
  3. f90 ファイルをコンパイルし、実行ファイルを実行する。
  4. gnuplot を起動し、そのコマンドラインで

    > load '3objs_move.txt'

    と打つと動画が見られる。

コメント：

解説


f90 を使える人であれば、解説がなくともソースを見ればわかるとは思いますが、一応少しだけ書いておきます。

f90 ファイルをコンパイル・実行すると、計算されたデータが 3objs.dat というファイルに書き出されます。
このファイルはプログラムを実行した後は不要ですので消しても問題はありません。
（約 54 MB：HDD の容量が気になる場合は消してください）
書き出された後、すぐに読み出されます。このとき、間引きされます。
そして、一定時間ごとのデータが数字名のファイルに書き出されます。
よって、数字名のファイル（拡張子なし）のファイルが 1990 個できるはずです。
このとき、便宜上 number というファイルも作られますが、すぐに使用されるので 3objs.dat とともに不要のファイルになります。
一定時間ごとのデータと書きましたが、実際は1ファイル当たり 50 時点のデータを書き出しています。
これにより、gnuplot で描画するとき、処理速度を遅くし、かつ残像（軌跡）が見やすくなります。
１つのデータに1点分のデータしか入れないと、これらのファイルがあっという間に処理されてしまい、動きがよくわかりません。

3objs_move.txt というファイルには、gnuplot でファイルを次々に呼び出すコマンドが書かれています。
上のをコピーして使えば必要ないですが、このファイルの内容を書き出すプログラムはこれ（f90）です。

数値計算手段として、4次の Runge-Kutta 法を使用しています。


実験・補足


プログラムをいじる場合、以下のようにしてください。
（！なお、いじった場合、コンパイル・実行するのを忘れないでください）
まず、星がいくつかある場合（上のプログラムでは3つ）、それぞれの星の情報は配列の番号で表されています。
星 1 の質量は m（1）、x座標、y座標、z座標はそれぞれ a（1）、b（1）、c（1）、同様に速度はvx（1）、vy（1）、vz（1）です。

星の数を変えたい場合、3行目の star= についている数字を変えてください。（整数）
そして、それにともなって、質量・初期位置・初速度の配列を増やし（減らし）てください。
例えば、星を4つにしたい場合は以下のようにします。（4箇所）

character(len=4)::num(0:1990)
integer,parameter::n= 400000,s=50000
integer,parameter::star=4                          ! 星の数
real,parameter:: h = 0.003                         ! 時間間隔
real,parameter:: pi=3.1415926536
real :: m(star),vx(star),vy(star),vz(star),k1,k2,k3,k4,l1,l2,l3,l4
real :: a(star),b(star),c(star),next_a(star),next_b(star),next_c(star),at(star,s),bt(star,s),ct(star,s)

! 質量
m(1)=1000.0
m(2)=1000.0
m(3)=1000.0
m(4)=700.0

! 初期位置
a(1)=50.0 ;b(1)=0 ;c(1)=0
a(2)=-25.0 ;b(2)=25.0*sqrt(3.0) ;c(2)=0
a(3)=-25.0 ;b(3)=-25.0*sqrt(3.0) ;c(3)=0
a(4)=0 ;b(4)=0 ;c(4)=100.0

! 初期速度
vx(1)=0 ;vy(1)=2.0 ;vz(1)=1.0
vx(2)=-2.0*sqrt(3.0)/2.0 ;vy(2)=-2.0/2.0 ;vz(2)=-1.0
vx(3)=2.0*sqrt(3.0)/2.0 ;vy(3)=-2.0/2.0 ;vz(3)=0
vx(4)=0 ;vy(4)=0 ;vz(4)=0

open(21,file='3objs.dat')

do j = 0, n-1
   do k=1,star
      k1 = vx(k)
・・・・・・・・

配列番号を4にするのを忘れないでください。よく忘れます。
あと、初期位置がダブらないように気をつけてください。
コピー・貼り付けで書くとたまに忘れたりします。

5個以上のときも同様にしたらできます。

また、時間間隔（h）を上げれば、より先の結果を見ることができます。（その分描画時の動きは速くなる）
ただし、あまり上げすぎると計算が不安定になります。

gnuplot での定義域（見ている範囲）を変えたい場合、3objs_move.txt を開き、
下の太字の数字を変えてください。（[下限：上限]）

 set xrange [-600:600]
 set yrange [-600:600]
 set zrange [-400:400]
splot '0'
 pause 0.8
splot '1'
 pause 0.8
splot '2'
 pause 0.8
splot '3'
・・・・・・・

x, y, z は比率が同じくらいになるようにした方がいいです。