      PROGRAM centralforce
!Numerically solves the following planetary problem via the Runge-Kutta 4 method:
!
!Solves the Mercury perihelion precession with the Sun fixed at the origin
!and the outer planets are replaced by a ring centered at the origin with
!the radius being the orbit of Jupiter and the total mass is the Jupiter mass.
!
!The simulation are in astronomical units with time in years/(2*pi) and length in AU.
      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           :: mJ,RJ
      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(2)!Mercurty position
      REAL*16           :: v(2)!Mercury velocity
      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(4,4)!auxiliar variables
      REAL*16           :: x,y!auxiliar variables
      OPEN(200,file='precession-rate-Mercury.dat',status='unknown')
      write(200,*)'program approx-Mercury-precession.f95: implements RK4 for the Sun-Mercury-Jupiter system'
      dt   = 1.q-5
      dtm  = dt / 2.q0
      dt2  = dt*dt
      aux  = 10.q0!intended time simulated in years
      imax = int(2.q0*pi*aux/dt)!total no. of iterations
      
      write(200,*)'dt =',dt
      write(200,*)'time simulated =',imax*dt/(2.q0*pi)
      write(200,*)''
      
      mJ = 50.q0 * 3.18q2 / 3.33q5!Jupiter mass in terms of solar mass (Earth mass = 1)
      RJ = 5.2q0!Radius of the ring representing the outer planets (mainly Jupiter)
      
      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(1) = x0!Mercury initial conditions
      r(2) = y0
      v(1) = vx0
      v(2) = vy0
      
      MSdistance(1) = -1.q0!initialize with a negative value for convinience
      MSdistance(2) = r(1)**2.q0!This will be the 1st aphelion
      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. (i=1 is simply the initial condition)
         F(1,1) = v(1)
         F(2,1) = v(2)
         k  = 1!1st term in RK4
         x  = r(1)
         y  = r(2)
         CALL accel(mJ,k,F,x,y,RJ)
         
         F(1,2) = v(1) + dtm*F(3,1)
         F(2,2) = v(2) + dtm*F(4,1)
         
         k = 2!2nd term in RK4
         x = r(1) + dtm*F(1,k-1)
         y = r(2) + dtm*F(2,k-1)
         CALL accel(mJ,k,F,x,y,RJ)
         
         F(1,3) = v(1) + dtm*F(3,2)
         F(2,3) = v(2) + dtm*F(4,2)
         
         k = 3!3rd term in RK4
         x = r(1) + dtm*F(1,k-1)
         y = r(2) + dtm*F(2,k-1)
         CALL accel(mJ,k,F,x,y,RJ)
         
         F(1,4) = v(1) + dt*F(3,3)
         F(2,4) = v(2) + dt*F(4,3)
         k = 4!4th term in RK4
         x = r(1) + dt*F(1,k-1)
         y = r(2) + dt*F(2,k-1)
         CALL accel(mJ,k,F,x,y,RJ)
         
         r(1) = r(1) + dt*(F(1,1) + 2.q0*(F(1,2)+F(1,3)) + F(1,4))/6.q0
         r(2) = r(2) + dt*(F(2,1) + 2.q0*(F(2,2)+F(2,3)) + F(2,4))/6.q0
         v(1) = v(1) + dt*(F(3,1) + 2.q0*(F(3,2)+F(3,3)) + F(3,4))/6.q0
         v(2) = v(2) + dt*(F(4,1) + 2.q0*(F(4,2)+F(4,3)) + F(4,4))/6.q0
         MSdistance(3) = r(1)**2.q0 + r(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. (i=1 corresponds to t=0)
            aphelion(2,krev) = atan2(MSposition(2),MSposition(1))!aphelion angle
            
            write(6,*) atan(MSposition(2),MSposition(1)),MSposition(2),MSposition(1)
!            write(6,*)krev,aphelion(2,krev),MSposition(1),MSposition(2)
            
            aphelion(3,krev) = sqrt(MSposition(1)**2 + MSposition(2)**2)!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(1)!store the x-distance between Mercury and Sun
         MSposition(2) = r(2)!store the y-distance between Mercury and Sun
      enddo!end of the loop for the RK4 method
      
      OPEN(unit=100,file='aphelion-approx.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(mJ,k,F,x,y,RJ)
      IMPLICIT NONE
      INTEGER,INTENT(in)    :: k
      REAL*16,INTENT(in)    :: x,y,mJ,RJ
      REAL*16,INTENT(inout) :: F(4,4)
      REAL*16               :: d,ax,ay,aux
      
      aux = sqrt(x*x+y*y)
      d   = aux**3.q0!this is d^3
      
      F(3,k) = -x/d!Mercury acceleration due to the Sun
      F(4,k) = -y/d!Mercury acceleration due to the Sun
      
      d  = aux / RJ
      d  = d + 9.q0*(d**3.q0)/8.q0 + 75.q0*(d**5.q0)/64.q0 + 1225.q0*(d**7.q0)/1024.q0
      d  = d/(2.q0*RJ*RJ)
      ax = d*x/aux
      ay = d*y/aux
      
      F(3,k) = F(3,k) + mJ*ax!correction due to the ring-Jupiter
      F(4,k) = F(4,k) + mJ*ay!correction due to the ring-Jupiter
      
      RETURN
      END SUBROUTINE accel
