program threebody

   implicit none

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

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

   integer step, nstep, count
   integer ip, jp

   Real(8) time, maxrevo, 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.0d8
   maxrevo = 5.0d4
   maxtime = maxrevo*2.0d0*PI
   
!---- 時刻と出力回数を 0 に設定する

   time = 0.0d0
   count = 0

!---- 各惑星の質量を決める

   factor = 1.0d0
   mp(1)  = 1.0d0
   mp(2)  = 9.548d-4*factor
   mp(3)  = 2.859d-4*factor

   write(*, *) mp(2), mp(3)

!---- 初期に惑星を次のように配置する

   xp(1) = 0.0d0
   yp(1) = 0.0d0

   xp(2) = 1.0d0
   yp(2) = 0.0d0

   xp(3) = 1.837d0
   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年度惑星科学実習用として、
! --- 今枝佑輔氏が作成されたプログラムを基にして作成されました。