program euler

!  Program in Fortran90 to compute the evolution of the two-body
!  problem with the first-order Euler method.

!  Wladimir Lyra, March 19, 2021.
!  Astro 506, New Mexico State University.  

  implicit none

!  Declaration of variables.
  
  integer, parameter :: npar = 2, ndim = 3
  real, parameter :: pi=3.14159265358979323846264338327950d0

  real, dimension(npar,ndim) :: x=0.,v=0.,dxdt=0.,dvdt=0.
  real, dimension(npar) :: mass
  
  integer :: ixp=1,iyp=2,izp=3,iplanet=2,istar=1,it,nt,i,j,k,it0
  real :: t,dt,q,e,norbits,torb,mean_motion,tmax,rr2,Gnewton=1.0
!
  real :: start, finish

! Start the timer.   

  call cpu_time(start)

!  These can be put in a input list.
  
  q = 0.
  e = 0.
  dt = 1e-3
  it0 = 10
  norbits = 10.

! Masses, positions, and velocities.
!
  call initial_conditions(q,e,mass,x,v,iplanet,istar,ixp,iyp,mean_motion,npar,ndim)

! Calculate the number of timesteps nt given the (fixed) timestep dt.

  torb = 2*pi/mean_motion
  tmax  = norbits * torb
  nt = int(tmax/dt)
!
  do it=1,nt

! Reset the derivative arrays.

     dxdt = 0
     dvdt = 0

! Calculate the evolution of particle k.     

     do k=1,npar

! dx/dt is simply the velocity.

        dxdt(k,:) = v(k,:)

! For the acceleration, sum over the gravity
! of all other massive particles in the system.        

        do j=1,npar
           if (j/=k) then
              ! Calculate the scalar separation.
              rr2=0
              do i=1,ndim
                 rr2 = rr2 + (x(k,i)-x(j,i))**2 
              enddo

! Add gravity from particle j on particle k.

              dvdt(k,:) = dvdt(k,:) - Gnewton*mass(j)*(x(k,:)-x(j,:))*rr2**(-1.5)
           endif
!
        enddo
     enddo

! Update position and velocity.

     x = x + dxdt*dt
     v = v + dvdt*dt

! Advance time. 

     t = t + dt

! Write diagnostics to standard output. 

     if (mod(it,it0)==0) then
        print*,it,t,dt,&
             x(istar,ixp),x(istar,iyp),x(iplanet,ixp),x(iplanet,iyp),&
             v(istar,ixp),v(istar,iyp),v(iplanet,ixp),v(iplanet,iyp)
     endif
  enddo

! Close the timer and print the wall time. 
  
  call cpu_time(finish)
  print*,"Wall time = ",finish-start," seconds."
 
endprogram euler
!**********************************************************************************
subroutine initial_conditions(q,e,mass,x,v,iplanet,istar,ixp,iyp,mean_motion,npar,ndim)
!
! Initial condition subroutine. 
!
!  The position of the center of mass is 
!
!     xcm = (m1*x1 + m2*x2)/(m1+m2)
!
!  Use barycentric coordinates where xcm = 0. Initialize the planet and 
!  star at the x direction. The velocity of the center of mass is 
!
!     vcm = (m1*v1 + m2*v2)/(m1+m2)
!
!  use baricentric coordinates where vcm = 0. The velocities in barycentric 
!  coordinates are 
!
!    v_planet =  m1 * vcirc  * sqrt((1-e)/(1+e))
!    v_star   = -m2 * vcirc  * sqrt((1-e)/(1+e))
! 
!  where the circular velocity is vcirc = n*a, with n the mean motion and 
!  a is the semimajor axis. If the particles positions are apocenter, then 
!  (xp-xs) = a*(1+e). Given the units, xp-xs=1. The mean motion is n=G(m1+m2)/a**1.5.
!
! 
  real, dimension(npar,ndim) :: x,v
  real, dimension(npar) :: mass
  integer :: npar,ndim,ixp,iyp,istar,iplanet
  real :: q,e,sma,vcirc,mean_motion
  
  intent(in) :: q,e,npar,ndim,ixp,iyp,istar,iplanet
  intent(out) :: x,v,mass,mean_motion
  
  mass(iplanet) = q
  mass(istar)   = 1-q

  x(istar  ,ixp) =  -q
  x(iplanet,ixp) = 1-q

  sma = 1/(1+e)
  mean_motion = sma**(-1.5)
  vcirc = mean_motion * sma
  
  v(istar  ,iyp) =   -q  * vcirc * sqrt((1-e)/(1+e))
  v(iplanet,iyp) = (1-q) * vcirc * sqrt((1-e)/(1+e))
  
endsubroutine initial_conditions
!**********************************************************************************