program TwobodyRK

   implicit none

!---- 惑星の位置の時間進化と速度の時間進化を次々とつなげていって惑星が
!---- 長い時間をかけてどのような軌道をたどるか調べてみよう。前に作った
!---- プログラムは惑星の楕円軌道があまり正確におえていなかった。
!---- そこで、惑星の位置や速度を計算する時に使う速度や力を時間が dt だけ
!---- 進んだ位置の値をつかってdt進む間の平均的な速度、力の代表値だと思って
!---- 計算するプログラムに書き換える。

!---- 各種変数を定義する

   integer step
   real(8) time, dt, dth
   real(8) xp, xph, yp, yph, rp, rph
   real(8) vx, vxh, vy, vyh
   real(8) fx, fxh, fy, fyh, fr, frh 
   real(8) PI

!---- Main

   PI = 3.1415926535898d0

!---- 時刻を 0 に設定する

   time = 0.0d0

!---- 初期に惑星を(1,0)に配置する

   xp = 1.0d0
   yp = 0.0d0

!---- おいた惑星を ypの正の方向に初速を与えて放り出す

   vx = 0.0d0
   vy = 1.1d0

!---- 与えた初期値を画面に表示する

   write(*,"(A,F10.4,A,F10.4)") "xp=", xp ," yp=", yp
   write(*,"(A,F10.4,A,F10.4)") "vx=", vx ," vy=", vy

!---- 進める時間 dt を決める

   dt = 2.0d0*PI*1.0d-4

   write(13,'(5(1X,E15.7))') time, xp, yp, vx, vy

!---- くり返し計算をするために Do ループで囲む

   do step = 1, 20000

!---- 惑星にかかる力を計算する(万有引力)
!---- 今簡単のため万有引力定数を1, 太陽 + 惑星の質量を 1 とする
!---- 太陽惑星間の距離は

      rp = Sqrt(xp*xp + yp*yp)

!---- 太陽が惑星を引っ張る力は万有引力の法則から

      fr = -1.0d0/(rp*rp)

!---- 上で求めた力は太陽から地球に向かう方向に働く力
!---- これを x方向の力と y方向の力に分解する

      fx = fr*xp/rp
      fy = fr*yp/rp

!---- 時刻を一気に dt 進めるのでなく 半分だけ進めてみる

      dth = dt*0.5d0
      xph = xp + vx*dth
      yph = yp + vy*dth
      vxh = vx + fx*dth
      vyh = vy + fy*dth

!---- 時刻が半分進んだ位置での惑星にかかる力を計算する

      rph = Sqrt(xph*xph+yph*yph)
      frh = -1.0d0/(rph*rph)
      fxh = frh*xph/rph
      fyh = frh*yph/rph

!---- 時刻が半分進んだ位置での 速度と力を使って 元の位置から
!---- の時間発展を計算する

      time = time + dt
      xp   = xp + vxh*dt
      yp   = yp + vyh*dt
      vx   = vx + fxh*dt
      vy   = vy + fyh*dt

!---- 位置と速度をファイルに書き出す
!---- 書き出し方はおまじないみたいなものをつかう
!---- これを使うとファイル名は fort.13 になる
!---- 時刻 位置xp 位置yp 考えている粒子の速度vx, 速度vy の順で書き出す

      write(13,'(5(1X,E15.7))') time, xp, yp, vx, vy

   end do

!---- プログラム終了

end program TwobodyRK

! --- このプログラムは2006 - 2008年度惑星科学実習用として
! --- 今枝佑輔氏が作成されたプログラムを基にして作成されました。