      PROGRAM centralforce
!Numerically solves for a particle in a central force ~ 1/r^2 (1 + alpha/r^2)
!via the Runge-Kutta 4 method
!
!The length are in astronomical units (AU) and the time is in years/(2*pi) units.
!The energy is in units of (2 Pi AU / year)^2 * (mass of the planet)
      IMPLICIT NONE
      INTEGER           :: i,j,k,ialpha
      INTEGER           :: imax
      INTEGER,PARAMETER :: nmax = 1000000!max. no. of iterations
      REAL*8,PARAMETER  :: pi   = 4.q0*atan(1.q0)
      REAL*8            :: e,a!excentricity and major semi-axis
      REAL*8            :: x0,y0!initial position
      REAL*8            :: vx0,vy0!initial velocity 
      REAL*8            :: dt,dtm,dt2!time discretization
      REAL*8            :: fi = -2!interaction exponent
      REAL*8            :: xmin,ymin,xmax,ymax
      REAL*8            :: t!time
      REAL*8            :: alpha!relativistic correction
      REAL*8            :: x(nmax),y(nmax)!position
      REAL*8            :: vx(nmax),vy(nmax)!velocity
      REAL*8            :: aux!auxiliar variables
      REAL*8            :: F11,F12,F13,F14!auxiliar variables
      REAL*8            :: F21,F22,F23,F24!auxiliar variables
      REAL*8            :: F31,F32,F33,F34!auxiliar variables
      REAL*8            :: F41,F42,F43,F44!auxiliar variables
      REAL*8            :: aphelion(3,1000)!aphelion angle vs time
      REAL*8            :: precessionrate(2,9)!precession rate vs alpha
      CHARACTER         :: u(1)
      dt   = 2.q-5
      dtm  = dt / 2.q0
      dt2  = dt*dt
      imax = 800000      
      if(imax.gt.nmax)then
         write(6,*)'increase nmax to',imax
         stop
      endif
      
      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
      vx0 = 0.q0!Initial velocity at aphelion
      vy0 = sqrt( (1.q0 - e) / x0 )
      
      alpha = 3.q-2!2 times the initial value of alpha
      do ialpha=1,9!loop for different values of alpha
         alpha = alpha / 3.q0!change alpha
      
         x(1) = x0!initial condition
         y(1) = y0
      
         vx(1) = vx0!initial condition
         vy(1) = vy0
         do i=1,imax-1!loop for the RK4 method
            j = i + 1!auxiliar variable representing t+dt
       
            F11 = vx(i)
            F21 = vy(i)
            CALL accel( F31,F41,fi,alpha,x(i),vx(i),y(i),vy(i) )
         
            F12 = vx(i) + dtm*F31
            F22 = vy(i) + dtm*F41
            CALL accel( F32,F42,fi,alpha,x(i)+dtm*F11,vx(i)+dtm*F31,y(i)+dtm*F21,vy(i)+dtm*F41 )
         
            F13 = vx(i) + dtm*F32
            F23 = vy(i) + dtm*F42
            CALL accel( F33,F43,fi,alpha,x(i)+dtm*F12,vx(i)+dtm*F32,y(i)+dtm*F22,vy(i)+dtm*F42 )
         
            F14 = vx(i) + dt*F33
            F24 = vy(i) + dt*F43
            CALL accel( F34,F44,fi,alpha,x(i)+dt*F13,vx(i)+dt*F33,y(i)+dt*F23,vy(i)+dt*F43 )
         
            x(j)  = x(i)  + dt*(F11 + 2.q0*(F12 + F13) + F14)/6.q0
            y(j)  = y(i)  + dt*(F21 + 2.q0*(F22 + F23) + F24)/6.q0
            vx(j) = vx(i) + dt*(F31 + 2.q0*(F32 + F33) + F34)/6.q0
            vy(j) = vy(i) + dt*(F41 + 2.q0*(F42 + F43) + F44)/6.q0
         enddo
         u(1) = char(48+ialpha)
         OPEN(unit=100,file='energy-'//u(1)//'.dat',status='unknown')
         OPEN(unit=101,file='distance-'//u(1)//'.dat',status='unknown')
         OPEN(unit=102,file='trajectory-difference-'//u(1)//'.dat',status='unknown')
         OPEN(unit=103,file='trajectory-'//u(1)//'.dat',status='unknown')  
         OPEN(unit=104,file='aphelion-'//u(1)//'.dat',status='unknown')
         do i=1,imax,100
            t = real(i-1,8)*dt/(2.q0*pi)
            CALL energy(aux,fi,alpha,x(i),vx(i),y(i),vy(i)) 
            write(100,*)t,aux
         
            write(101,*)t,sqrt((x(i))**2 + y(i)**2)/x0
         
            write(102,*)t,abs( abs(y(i))/(a*sqrt(1.q0-e*e)) - sqrt(abs(1.q0-(x(i)/a-e)**2)) )
         
            write(103,*)x(i),y(i)
         enddo
         do i=1,imax
            vx(i) = sqrt( x(i)**2 + y(i)**2 )!now, vx becomes the distance from the SUN
         enddo
         aphelion(1,1) = 0.q0!time at the 1st aphelion
         aphelion(2,1) = 0.q0!angle of the 1st aphelion
         aphelion(3,1) = x0!distance of the 1st aphelion
         k             = 1
         do i=2,imax-1
            if(vx(i).gt.vx(i-1).and.vx(i).gt.vx(i+1))then
               k             = k + 1!this is the kth aphelion
               aphelion(1,k) = real(i-1,8)*dt / (2.q0*pi)!time
               aphelion(2,k) = atan(y(i),x(i))!aphelion angle
               do while(aphelion(2,k).lt.aphelion(2,k-1))
                  aphelion(2,k) = aphelion(2,k) + 2.q0*pi
               enddo
               aphelion(3,k) = vx(i)!aphelion distance
            endif
         enddo
         do i=1,k
            write(104,*)(aphelion(j,i),j=1,3)
         enddo            
         CLOSE(100)
         CLOSE(101)
         CLOSE(102)
         CLOSE(103)
         CLOSE(104)
         precessionrate(1,ialpha) = alpha
         aux = 0.q0
         t   = 0.q0!just another auxiliar variable
         do i=1,k!loop for the linear regression aphelionangle = precessionrate * t
            aux = aux + aphelion(2,i)*aphelion(1,i)
            t   = t   + aphelion(1,i)*aphelion(1,i)
         enddo
         precessionrate(2,ialpha) = aux/t
      enddo!end of the loop for different value of alpha
      OPEN(unit=100,file='precessionrate-alpha.dat',status='unknown')
      do i=1,9
         write(100,*)(precessionrate(j,i),j=1,2)
      enddo
      CLOSE(100)
      aux = 0.q0!computing the precession rate by the linear fit of log(precession ) = log(aux) + log(alpha)
      do i=1,9
         aux = aux + log(precessionrate(2,i)/precessionrate(1,i))
      enddo
      aux = exp(aux / 9.q0) !this is precessionrate = aux * alpha
      write(6,*)'precession rate =',aux,'*alpha (in rad/year)'
      aux = aux * 1.1q-8 !using the alpha for Mercury
      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,'seconds of degrees/century'
      STOP
      END PROGRAM centralforce
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!ENDO OF PROGRAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE accel(ax,ay,fi,alpha,x,vx,y,vy)
      IMPLICIT NONE
      REAL*8,INTENT(out)  :: ax,ay
      REAL*8,INTENT(in)   :: fi
      REAL*8,INTENT(in)   :: alpha
      REAL*8,INTENT(in)   :: x,y
      REAL*8,INTENT(in)   :: vx,vy
      REAL*8              :: d
      d  = sqrt(x*x + y*y)
      d  = d**(fi - 1.q0)*( 1.q0 + alpha/(d*d) )
      
      ax = - x*d
      ay = - y*d
      RETURN
      END SUBROUTINE accel
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE energy(E,fi,alpha,x,vx,y,vy)
      IMPLICIT NONE
      REAL*8,INTENT(out)  :: E
      REAL*8,INTENT(in)   :: fi
      REAL*8,INTENT(in)   :: alpha
      REAL*8,INTENT(in)   :: x,y
      REAL*8,INTENT(in)   :: vx,vy
      REAL*8              :: d
      d = dsqrt(x*x + y*y)
      E = d**(fi + 1.q0)/(fi + 1.q0) + alpha*d**(fi - 1.q0)/(fi - 1.q0) + 0.5q0 * (vx*vx + vy*vy)
      RETURN
      END SUBROUTINE energy
