program Twobody

   implicit none

!---- 惑星の位置の時間進化と速度の時間進化を次々とつなげていって
!---- 惑星が長い時間をかけてどのような軌道をたどるか調べてみよう
!---- 今までに作ったプログラムの時間発展を計算する部分くり返し計算を
!---- するように改造する。

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

   integer step
   real(8) time, dt
   real(8) xp, yp, rp
   real(8) vx, vy
   real(8) fx, fy, fr 
   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

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

   do step = 1, 20000

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

!---- 太陽惑星間の距離は

      rp = Sqrt(??????????)

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

      fr = ??????????

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

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

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

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

!---- 時刻が time+dt での粒子の位置は?
!---- 新しい位置 = 昔の位置 + 速度経過時間

      time = time + dt
      xp = ??????????
      yp = ??????????
      vx = ??????????
      vy = ??????????

   end do

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

end program Twobody

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