      PROGRAM centralforce
!Numerically solves for 3 particles interacting via a central force ~ 1/r^2 
!via the Runge-Kutta 4 method
!
!The length are in astronomical units (AU) and the time is in years/(2*pi) units.
!Measures the Perihelion precession for Mercury
      IMPLICIT NONE
      INTEGER           :: i,j,k
      INTEGER           :: imax
      REAL*16,PARAMETER :: pi = 4.q0*atan(1.q0)
      REAL*16           :: e,a!excentricity and major semi-axis
      REAL*16           :: m2,m3
      REAL*16           :: x0,y0!initial position
      REAL*16           :: vx0,vy0!initial velocity 
      REAL*16           :: dt,dtm,dt2!time discretization
      REAL*16           :: fi = -2!interaction exponent
      REAL*16           :: aux,bux!auxiliar variables
      REAL*16           :: r(3,2)!position of the 3 particles
      REAL*16           :: v(3,2)!velocity of the 3 particles
      REAL*16           :: aphelion(3,1000)!time, aphelion angle, aphelion distance
      INTEGER           :: krev!counts the number of aphelions measured during the simulation
      REAL*16           :: MSdistance(3)!Mercury-Sun distance
      REAL*16           :: MSposition(2)!Mercury position relative to the Sun
      REAL*16           :: F(3,4,4)!auxiliar variables
      REAL*16           :: x1,x2,x3,y1,y2,y3!auxiliar variables
      OPEN(200,file='precession-rate-Mercury-Jupiter.dat',status='unknown')
      write(200,*)'program 3body.f95: implements RK4 for the Sun-Mercury-Jupiter system'
      dt   = 1.q-8
      dtm  = dt / 2.q0
      dt2  = dt*dt
      aux  = 2.q0!time simulated
      imax = int(2.q0*pi*aux/dt)!total no. of iterations
      
      write(200,*)'dt =',dt
      write(200,*)'time simulated =',imax*dt/(2*pi)
      write(200,*)''
      
      m2 = 5.5q-2 / 3.33q5!Mercury mass in terms of solar mass (Earth mass = 1)
      m3 = 3.18q2 / 3.33q5!Jupiter mass in terms of solar mass (Earth mass = 1)
      
      e   = 0.206q0!Mercury excentricity
      a   = 0.39q0!Mercury major semiaxis
      x0  = (1.q0 + e) * a!Initial position from the Sun: the aphelion
      y0  = 0.q0!Initial position from the Sun: the aphelion
      vx0 = 0.q0!Initial velocity at aphelion
      vy0 = sqrt( (1.q0 - e) / x0 )!Initial velocity at aphelion
      
      write(200,*)'Mercury major semiaxis =',a
      write(200,*)'Mercury excentricity =',e
      write(200,*)'Initial x position =',x0
      write(200,*)'Initial y position =',y0
      write(200,*)'Initial x velocity =',vx0
      write(200,*)'Initial y velocity =',vy0
      write(200,*)''
      
      r(2,1) =  x0!Mercury initial conditions (from the Sun)
      r(2,2) =  y0
      v(2,1) =  vx0
      v(2,2) =  vy0
      
      e   = 0.049q0!Jupiter excentricity
      a   = 5.2q0!Jupiter major semiaxis
      x0  = (1.q0 + e) * a!Initial position from the Sun: the aphelion
!      x0  = -(1.q0 - e) * a!Initial position from the Sun: the perihelion
      y0  = 0.q0
      vx0 = 0.q0
      vy0 = sqrt( (1.q0 - e) / x0 )!Initial velocity at aphelion
!      vy0 = -sqrt( -(1.q0 + e) / x0 )!Initial velocity at perihelion
      
      write(200,*)'Jupiter major semiaxis =',a
      write(200,*)'jupiter excentricity =',e
      write(200,*)'Initial x position =',x0
      write(200,*)'Initial y position =',y0
      write(200,*)'Initial x velocity =',vx0
      write(200,*)'Initial y velocity =',vy0
      write(200,*)''
      
      r(3,1) = x0!Jupiter initial conditions (from the Sun)
      r(3,2) = y0
      v(3,1) = vx0
      v(3,2) = vy0
      
      r(1,1) = 0.q0!Sun at origin
      r(1,2) = 0.q0
      v(1,1) = 0.q0!Sun at rest
      v(1,2) = 0.q0
      
      aux = ( r(1,1) + m2*r(2,1) + m3*r(3,1) ) / (1.q0 + m2 + m3)!Center of mass x-coordinate
      do i=1,3
         r(i,1) = r(i,1) - aux!make the center of mass at the origin
      enddo
      aux = ( v(1,2) + m2*v(2,2) + m3*v(3,2) ) / (1.q0 + m2 + m3)!Center of mass y-velocity
      do i=1,3
         v(i,2) = v(i,2) - aux!make the center of mass at rest
      enddo
      
      MSdistance(1) = -1.q0!initialize with a negative value for convinience
      MSdistance(2) = (r(2,1) - r(1,1))**2.q0!This will be the 1st aphelion^2
      MSposition(1) = x0!store the x-distance between Mercury and Sun
      MSposition(2) = y0!store the y-distance between Mercury and Sun
      krev          = 0
      do i=2,imax!loop for the RK4 method
         
         do k=1,3!particle label
            F(k,1,1) = v(k,1)
            F(k,2,1) = v(k,2)
         enddo
         k = 1!1st term in RK4
         x1 = r(1,1)
         y1 = r(1,2)
         x2 = r(2,1)
         y2 = r(2,2)
         x3 = r(3,1)
         y3 = r(3,2)
         CALL accel(m2,m3,k,F,x1,y1,x2,y2,x3,y3)
         
         do k=1,3!particle label
            F(k,1,2) = v(k,1) + dtm*F(k,3,1)
            F(k,2,2) = v(k,2) + dtm*F(k,4,1)
         enddo
         k = 2!2nd term in RK4
         x1 = r(1,1) + dtm*F(1,1,k-1)
         y1 = r(1,2) + dtm*F(1,2,k-1)
         x2 = r(2,1) + dtm*F(2,1,k-1)
         y2 = r(2,2) + dtm*F(2,2,k-1)
         x3 = r(3,1) + dtm*F(3,1,k-1)
         y3 = r(3,2) + dtm*F(3,2,k-1)
         CALL accel(m2,m3,k,F,x1,y1,x2,y2,x3,y3)
         
         do k=1,3!particle label
            F(k,1,3) = v(k,1) + dtm*F(k,3,2)
            F(k,2,3) = v(k,2) + dtm*F(k,4,2)
         enddo
         k = 3!3rd term in RK4
         x1 = r(1,1) + dtm*F(1,1,k-1)
         y1 = r(1,2) + dtm*F(1,2,k-1)
         x2 = r(2,1) + dtm*F(2,1,k-1)
         y2 = r(2,2) + dtm*F(2,2,k-1)
         x3 = r(3,1) + dtm*F(3,1,k-1)
         y3 = r(3,2) + dtm*F(3,2,k-1)
         CALL accel(m2,m3,k,F,x1,y1,x2,y2,x3,y3)
         
         do k=1,3!particle label
            F(k,1,4) = v(k,1) + dt*F(k,3,3)
            F(k,2,4) = v(k,2) + dt*F(k,4,3)
         enddo
         k = 4!4th term in RK4
         x1 = r(1,1) + dt*F(1,1,k-1)
         y1 = r(1,2) + dt*F(1,2,k-1)
         x2 = r(2,1) + dt*F(2,1,k-1)
         y2 = r(2,2) + dt*F(2,2,k-1)
         x3 = r(3,1) + dt*F(3,1,k-1)
         y3 = r(3,2) + dt*F(3,2,k-1)
         CALL accel(m2,m3,k,F,x1,y1,x2,y2,x3,y3)
         
         do k=1,3
            r(k,1) = r(k,1) + dt*(F(k,1,1)+2.q0*(F(k,1,2)+F(k,1,3))+F(k,1,4))/6.q0
            r(k,2) = r(k,2) + dt*(F(k,2,1)+2.q0*(F(k,2,2)+F(k,2,3))+F(k,2,4))/6.q0
            v(k,1) = v(k,1) + dt*(F(k,3,1)+2.q0*(F(k,3,2)+F(k,3,3))+F(k,3,4))/6.q0
            v(k,2) = v(k,2) + dt*(F(k,4,1)+2.q0*(F(k,4,2)+F(k,4,3))+F(k,4,4))/6.q0
         enddo
         MSdistance(3) = (r(2,1) - r(1,1))**2.q0 + (r(2,2) - r(1,2))**2.q0!new distance^2, update the distance at time step i
         if(MSdistance(2).gt.MSdistance(1).and.MSdistance(2).gt.MSdistance(3))then
            krev             = krev + 1
            aphelion(1,krev) = real(i-2,8)*dt / (2.q0*pi)!time at the kth aphelion, this happens at the previous time step: i-1
            aphelion(2,krev) = atan(MSposition(2),MSposition(1))!aphelion angle
            aphelion(3,krev) = sqrt(MSposition(1)**2.q0 + MSposition(2)**2.q0)!aphelion distance
         endif
         MSdistance(1) = MSdistance(2)!update the distance at time step i-2
         MSdistance(2) = MSdistance(3)!update the distance at time step i-1
         MSposition(1) = r(2,1) - r(1,1)!store the x-distance between Mercury and Sun
         MSposition(2) = r(2,2) - r(1,2)!store the y-distance between Mercury and Sun
      enddo!end of the loop for the RK4 method
      
      OPEN(unit=100,file='aphelion.dat',status='unknown')
      do i=1,krev
         write(100,*)(aphelion(j,i),j=1,3)
      enddo
      CLOSE(100)
      
      aux = 0.q0!computing the precession rate by the linear fit of precession angle = aux*time
      bux = 0.q0
      do i=1,krev!loop for the linear regression aphelionangle = precessionrate * t
         aux = aux + aphelion(2,i)*aphelion(1,i)
         bux = bux + aphelion(1,i)*aphelion(1,i)
      enddo
      aux = aux/bux!this is the precession rate
      write(6,*)'precession rate =',aux,'rad/year'
      write(200,*)'precession rate =',aux,'rad/year'
      aux = aux * 180.q0 / pi !transform radians -> degrees
      aux = aux * 3600 !transform degrees -> seconds
      aux = aux * 100.q0 !transform years -> centuries
      write(6,*)'precession rate =',aux,'arcseconds/century'
      write(200,*)'precession rate =',aux,'arcseconds/century'
      CLOSE(200)
      STOP
      END PROGRAM centralforce
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!ENDO OF PROGRAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE accel(m2,m3,k,F,x1,y1,x2,y2,x3,y3)
      IMPLICIT NONE
      INTEGER,INTENT(in)    :: k
      REAL*16,INTENT(in)    :: x1,y1,x2,y2,x3,y3,m2,m3
      REAL*16,INTENT(inout) :: F(3,4,4)
      REAL*16               :: d,ax,ay
      
      d  = ( (x1-x2)**2.q0 + (y1-y2)**2.q0 )**1.5q0
      ax = - (x1-x2)/d
      ay = - (y1-y2)/d
      
      F(1,3,k) =  m2*ax
      F(1,4,k) =  m2*ay
      F(2,3,k) = -ax!recall m1=1
      F(2,4,k) = -ay!recall m1=1
      
      d  = ( (x1-x3)**2.q0 + (y1-y3)**2.q0 )**1.5q0
      ax = - (x1-x3)/d
      ay = - (y1-y3)/d
      
      F(1,3,k) =  F(1,3,k) + m3*ax
      F(1,4,k) =  F(1,4,k) + m3*ay
      F(3,3,k) = -ax!recal m1=1
      F(3,4,k) = -ay!recal m1=1
      
      d  = ( (x2-x3)**2.q0 + (y2-y3)**2.q0 )**1.5q0
      ax = - (x2-x3)/d
      ay = - (y2-y3)/d
      
      F(2,3,k) = F(2,3,k) + m3*ax
      F(2,4,k) = F(2,4,k) + m3*ay
      F(3,3,k) = F(3,3,k) - m2*ax
      F(3,4,k) = F(3,4,k) - m2*ay
      RETURN
      END SUBROUTINE accel
