PROGRAM three_body_problem_in_rotating_frame_by_2_rungekutta

!--変数の宣言をする---------------------------------------------------
   implicit none
   integer :: p,q,i,j
   real(8) :: x,y,vx,vy,fx,fy
   real(8) :: xh,yh,vxh,vyh,fxh,fyh
   real(8) :: xn,yn,vxn,vyn,fxn,fyn
   real(8) :: x1,y1,x2,y2,r1,r2,xx,yy
   real(8) :: m1,m2,mu,mu1,mu2
   real(8) :: w,t,dt,maxrevo,Cj,Cj0,U,forbid,PI
   real(8) :: kx(4),ky(4),lx(4),ly(4)
   real(8) :: w1,w2,Tk1,Tk2

!--書き込むfileの準備-------------------------------------------------
   open(unit=100, file="2orbit.dat", status="unknown")
   open(unit=200, file="locate.dat", status="unknown" )
   open(unit=300, file="2jacobi.dat", status="unknown" )
   open(unit=400, file="forbidden.dat", status="unknown")
   write(400,*) " " !-forbidden.dat の中身のクリア 

!--定数の値を決める---------------------------------------------------

   w = 1.0d0
   PI = 3.1415926535898d0

  !--中心星と惑星の質量比を求める-----------------------------------
   m1 = 1.9891d33
   m2 = 1.8986d30

   mu = m2/(m1+m2)
   mu1 = 1.0d0-mu
   mu2 = mu

  !--中心星と惑星の位置を決める-------------------------------------
   x1 = -mu2
   y1 = 0.0d0

   x2 = mu1
   y2 = 0.0d0

   write(*,*) "factor"
   write(*,*) "mu1=",mu1,"mu2=",mu2
   write(*,*) " x1=",x1," y1=",y1
   write(*,*) " x2=",x2," y2=",y2
   write(*,*)
   write(200,*) x1,y1
   write(200,*) x2,y2

!--初期条件-----------------------------------------------------------
  !--粒子の初期位置および初速度を与える-----------------------------
   x = 1.4960d0 / 5.2026d0
   y = 0.0d0

   vx = 0.0d0
   vy = 1.57d0

   write(*,*) "initial condition"
   write(*,*) " x =",x," y =",y
   write(*,*) " vx=",vx," vy=",vy
   write(*,*)

!--jacobi定数から定まる粒子の運動の禁止領域を求める-------------------

  !--jacobi定数の計算-----------------------------------------------
   r1 = sqrt( (x-x1)**2.0d0 + y**2.0d0 )
   r2 = sqrt( (x-x2)**2.0d0 + y**2.0d0 )
   Cj = w**2.0d0/2.0d0*(x**2.0d0+y**2.0d0) + mu1/r1 + mu2/r2 - &
                (vx**2.0d0+vy**2.0d0)/2.0d0
   Cj0 = Cj

   write(*,*) "Cj=",Cj

  !--xとyを-2から2までroopさせる------------------------------------
   do i = 0,200
      do j = 0,200

        xx = -2.0d0+dble(i)*2.0d-2
        yy = -2.0d0+dble(j)*2.0d-2

        r1 = sqrt( (xx-x1)**2.0d0 + yy**2.0d0 )
        r2 = sqrt( (xx-x2)**2.0d0 + yy**2.0d0 )

       !--ある点(x,y)でのポテンシャル U とjacobi定数の差をとる----
        U = w**2.0d0/2.0d0*(xx**2.0d0+yy**2.0d0) + mu1/r1 + mu2/r2
        forbid = U - Cj

       !--差が0以下である(粒子の速度が0となる)場合、その点を書き出す-
        if (forbid .le. 0.0d0) then
          write(400,*) xx,yy,forbid
        end if

      end do !-j
   end do !-i

!--軌道計算-----------------------------------------------------------

  !--時間が maxrevo となるまで計算する------------------------------
   t = 0.0d0
   maxrevo = 2.0d0*PI * 1.0d2

   do while (t .lt. maxrevo)

    !--進める時間 dt を決める-------------------------------------
     !--現在の位置でそれぞれの天体の回りを円運動するときの周期の-
     !--1/1000と1.0d-3を比較して、最も小さいものを dt とする-----
      dt = 1.0d-3

      r1 = sqrt( (x-x1)**2.0d0 + y**2.0d0 )
      r2 = sqrt( (x-x2)**2.0d0 + y**2.0d0 )

      w1 = sqrt( mu1 / r1**3.0d0 )
      Tk1 = 2.0d0*PI/w1

      w2 = sqrt( mu2 / r2**3.0d0 )
      Tk2 = 2.0d0*PI/w2
      dt = min(dt,Tk1/1.0d3,Tk2/1.0d3)

    !--2 次のrunge-kutta法----------------------------------------
     !--現在の位置(x,y,t)での傾き(速度) kx(1),ky(1) を求める-----
     !--力についても同様に(x,y,vx,vy,t)での lx(1),ly(1) を求める-
     !--力 fx,fy は最後の Subroutine で計算している--------------
      call f(x,y,vx,vy,fx,fy)

      kx(1) = vx
      ky(1) = vy

      lx(1) = fx
      ly(1) = fy

     !--(x,y,t)から傾き(1)で時間を dt/2.0 だけ進める-------------
     !--進んだ先(x+kx(1)*dt/2.0, y+ky(1)*dt/2.0, t+dt/2.0)での傾き(2)を求める-
     !--力についても同様に求める---------------------------------
      xh = x + kx(1)*dt/2.0
      yh = y + ky(1)*dt/2.0
      vxh = vx + lx(1)*dt/2.0
      vyh = vy + ly(1)*dt/2.0

      call f(xh,yh,vxh,vyh,fxh,fyh)

      kx(2) = vxh
      ky(2) = vyh

      lx(2) = fxh
      ly(2) = fyh

     !--上で求めた傾き(2)で時間 dt 進める-------------------------
      x = x + kx(2)*dt
      y = y + ky(2)*dt
      vx = vx + lx(2)*dt
      vy = vy + ly(2)*dt

      write(100,*) x,y

    !--jacobi定数の計算--------------------------------------------
      r1 = sqrt( (x-x1)**2.0d0 + y**2.0d0 )
      r2 = sqrt( (x-x2)**2.0d0 + y**2.0d0 )
      Cj = (w**2.0d0)/2.0d0*(x**2.0d0+y**2.0d0) + mu1/r1 + mu2/r2 - &
                (vx**2.0d0+vy**2.0d0)/2.0d0

      write(300,*) t,Cj,abs(Cj-Cj0)

      t = t + dt
 
   end do !-while

!---------------------------------------------------------------------

 close(unit=100,status="keep")
 close(unit=200,status="keep")
 close(unit=300,status="keep")
 close(unit=400,status="keep")

 CONTAINS

!-力を求めるSubroutine------------------------------------------------
 !--本文中で call f(1,2,3,4,answer1,answer2) と書かれたとき、-------
 !--1,2,3,4 を X,Y,Vx,Vy へ代入して Fx,Fy を計算し、それらを--------
 !--answer1,answer2 として返す--------------------------------------
  SUBROUTINE f(X,Y,Vx,Vy,Fx,Fy)
    implicit none
    real(8) :: X,Y,Vx,Vy,Fx,Fy

      r1 = sqrt( (X-x1)**2.0d0 + Y**2.0d0 )
      r2 = sqrt( (X-x2)**2.0d0 + Y**2.0d0 )
      Fx = -mu1*(X+mu2)/(r1**3.0d0) - mu2*(X-mu1)/(r2**3.0d0) + &
                      2.0d0*w*Vy + w**2.0d0*X
      Fy = -mu1*Y/(r1**3.0d0) - mu2*Y/(r2**3.0d0) - &
                      2.0d0*w*Vx + w**2.0d0*Y

  END SUBROUTINE f
!---------------------------------------------------------------------

END PROGRAM