program threebody implicit none !---- 惑星の位置の時間進化と速度の時間進化を次々とつなげていって !---- 惑星が長い時間をかけてどのような軌道をたどるか調べてみよう !---- 今までに作ったプログラムの時間発展を計算する部分くり返し計算を !---- するように改造する !---- 各種変数を定義する integer step, nstep, count integer ip, jp Real(8) time, maxtime Real(8) dt, dth, dtmin Real(8) mp(3) Real(8) xp(3), xph(3), yp(3), yph(3) Real(8) vx(3), vxh(3), vy(3), vyh(3) Real(8) xg, yg Real(8) vxg, vyg Real(8) rp(3,3), rph(3,3) Real(8) vp, fr, frh Real(8) fx(3,3), fxh(3,3), fy(3,3), fyh(3,3) Real(8) fxtotal, fytotal, mtotal Real(8) factor Real(8) PI, EPSILON !---- Main PI = 3.1415926535898d0 EPSILON = 1.0d-30 nstep = 1.0d9 maxtime = 3.0d5 !---- 時刻を 0に設定する time = 0.0d0 !---- 各惑星の質量を決める factor = 1.0d0 mp(1) = 1.0d0 mp(2) = 9.5479d-4*factor mp(3) = 2.859d-4 *factor write(*, *) mp(2), mp(3) !---- 初期に惑星を(1,0)に配置する xp(1) = 0.0d0 yp(1) = 0.0d0 xp(2) = 1.0d0 yp(2) = 0.0d0 xp(3) = 9.5549d0/5.2026d0 yp(3) = 0.0d0 !---- おいた惑星を ypの正の方向に初速を与えて放り出す vx(1) = 0.0d0 vy(1) = 0.0d0 vx(2) = 0.0d0 vy(2) = sqrt(mp(1)+mp(2))/sqrt(Abs(xp(2))) vx(3) = 0.0d0 vy(3) = sqrt(mp(1)+mp(3))/sqrt(Abs(xp(3))) !---- 与えた初期値を画面に表示する write(*,'(e12.3)') factor write(*,*) write(*,"(A,F10.4,A,F10.4)") "xp1=",xp(1)," yp1=",yp(1) write(*,"(A,F10.4,A,F10.4)") "vx1=",vx(1)," vy1=",vy(1) write(*,*) write(*,"(A,F10.4,A,F10.4)") "xp2=",xp(2)," yp2=",yp(2) write(*,"(A,F10.4,A,F10.4)") "vx2=",vx(2)," vy2=",vy(2) write(*,*) write(*,"(A,F10.4,A,F10.4)") "xp3=",xp(3)," yp3=",yp(3) write(*,"(A,F10.4,A,F10.4)") "vx3=",vx(3)," vy3=",vy(3) write(*,*) !---- 重心を求める xg = 0.0d0 yg = 0.0d0 vxg = 0.0d0 vyg = 0.0d0 mtotal = 0.0d0 do ip = 1, 3 xg = xg + mp(ip)*xp(ip) yg = yg + mp(ip)*yp(ip) vxg = vxg + mp(ip)*vx(ip) vyg = vyg + mp(ip)*vy(ip) mtotal = mtotal + mp(ip) end do xg = xg/mtotal yg = yg/mtotal vxg = vxg/mtotal vyg = vyg/mtotal !---- 重心を原点にとるように座標を移動 do ip = 1, 3 xp(ip) = xp(ip) - xg yp(ip) = yp(ip) - yg vx(ip) = vx(ip) - vxg vy(ip) = vy(ip) - vyg end do !---- 初期設定終わり !---- 進める時間 dt を決める(最初は0とする) dt = 0.0d0 !---- くり返し計算をするために Do ループで囲む do step = 1, nstep dtmin = 1.0d11 !---- 惑星にかかる力を計算する(万有引力) !---- 今簡単のため万有引力定数を1 太陽の質量を1とする !---- 併せて、1ステップに進める微小時間を決める do ip = 1, 3 do jp = 1, 3 if (jp.ne.ip) then rp(jp,ip) = ?????????? fr = ?????????? fx(jp,ip) = fr*(xp(jp)-xp(ip))/rp(jp,ip) fy(jp,ip) = fr*(yp(jp)-yp(ip))/rp(jp,ip) vp = Sqrt((vx(jp)-vx(ip))**2.0d0 + & (vy(jp)-vy(ip))**2.0d0) dtmin = Min(dtmin, rp(jp,ip)/vp) else rp(jp,ip) = 0.0d0 fx(jp,ip) = 0.0d0 fy(jp,ip) = 0.0d0 endif end do end do !---- dt は 惑星が 1 公転するのを 1001 ステップで刻む幅にとる dt = 2.0d0*PI*dtmin*1.001d-3 !---- 時刻を一気に dt 進めるのでなく 半分だけ進めてみる dth = dt*0.5d0 do ip = 1, 3 fxtotal = 0.0d0 fytotal = 0.0d0 do jp = 1,3 fxtotal = fxtotal + fx(ip,jp) fytotal = fytotal + fy(ip,jp) end do xph(ip) = ?????????? yph(ip) = ?????????? vxh(ip) = ?????????? vyh(ip) = ?????????? end do !---- 時刻が半分進んだ位置での惑星にかかる力を計算する do ip = 1, 3 do jp = 1, 3 if (ip.ne.jp) then rph(jp,ip) = sqrt( (xph(jp)-xph(ip))**2.0d0 + & (yph(jp)-yph(ip))**2.0d0 ) frh = - mp(ip)/(rph(jp,ip)**2.0d0) fxh(jp,ip) = frh*(xph(jp)-xph(ip))/rph(jp,ip) fyh(jp,ip) = frh*(yph(jp)-yph(ip))/rph(jp,ip) else rph(jp,ip) = 0.0d0 fxh(jp,ip) = 0.0d0 fyh(jp,ip) = 0.0d0 endif end do end do !---- 時刻が半分進んだ位置での 速度と力を使って 元の位置から !---- の時間発展を計算する do ip = 1, 3 fxtotal = 0.0d0 fytotal = 0.0d0 do jp = 1, 3 fxtotal = fxtotal + fxh(ip,jp) fytotal = fytotal + fyh(ip,jp) end do xp(ip) = ?????????? yp(ip) = ?????????? vx(ip) = ?????????? vy(ip) = ?????????? end do !---- 時間をdt進めた時刻での惑星の位置と速度が求まったので !---- 時刻もdtだけ進める time = time + dt !---- 位置と速度をファイルに書き出す !---- 書き出し方はおまじないみたいなものをつかう !---- これを使うとファイル名は fort.14 になる !---- 時刻 位置xp 位置yp 考えている粒子の速度vx, 速度vy の順で書き出す if(mod(step, 50000).eq.1) then write(*, *) step write(14,'(13(1X,E15.7))') time, xp(1), yp(1), xp(2), & yp(2), xp(3), yp(3) !c$$$ & ,vx(1),vy(1),vx(2),vy(2),vx(3),vy(3) count = count +1 if (count.eq.5000) Goto 2001 endif !c$$$c---- 重心の位置が時間が進むとともにどの程度変化するかを調べる !c$$$c---- 正しくは重心の位置は時間とともに変化しない !c$$$c---- 元の位置からほとんど動かなければ、正しい計算ができている証拠 !c$$$ xg = 0.0d0 !c$$$ yg = 0.0d0 !c$$$ vxg = 0.0d0 !c$$$ vyg = 0.0d0 !c$$$ mtotal = 0.0d0 !c$$$ !c$$$ do ip = 1,3 !c$$$ !c$$$ xg = xg + mp(ip)*xp(ip) !c$$$ yg = yg + mp(ip)*yp(ip) !c$$$ vxg = vxg + mp(ip)*vx(ip) !c$$$ vyg = vyg + mp(ip)*vy(ip) !c$$$ mtotal = mtotal + mp(ip) !c$$$ !c$$$ end do !c$$$ !c$$$ xg = xg/mtotal !c$$$ yg = yg/mtotal !c$$$ vxg = vxg/mtotal !c$$$ vyg = vyg/mtotal !c$$$ !c$$$ write(15,'(5(1X,E15.7))') time, xg, yg, vxg, vyg if(time.gt.maxtime) Goto 1001 end do write(*,*) "Step >",nstep Stop 1001 Continue write(*,*) "time >",maxtime Stop 2001 Continue write(*,*) "File is full" Stop !---- プログラム終了 end program threebody ! --- このプログラムは2006 - 2008年度惑星科学実習用として、 ! --- 今枝佑輔氏が作成されたプログラムを基にして作成されました。