      PROGRAM centralforce
!Numerically solves for a particle in a central force ~ r^(fi)
!Implements the Euler-Crommer, Verlet, and Runge-Kutta 4 methods
!
!For the Kepler problem (fi = -2), 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
      INTEGER           :: framestep!step for discarding useless frames
      INTEGER           :: imax
      INTEGER,PARAMETER :: nmax = 1000000!max. no. of iterations
      REAL*16,PARAMETER  :: pi   = 4.q0*atan(1.q0)
      REAL*16            :: e,a!excentricity and major semi-axis
      REAL*16            :: x0,y0!initial position
      REAL*16            :: vx0,vy0!initial velocity 
      REAL*16            :: dt,dtm,dt2!time discretization
      REAL*16            :: fi!interaction exponent
      REAL*16            :: xmin,ymin,xmax,ymax
      REAL*16            :: t!time
      REAL*16            :: alpha!relativistic correction
      REAL*16            :: xRK(nmax),yRK(nmax)!Runge-Kutta
      REAL*16            :: vxRK(nmax),vyRK(nmax)!Runge-Kutta
      REAL*16            :: xEC(nmax),yEC(nmax)!Euler-Cromer
      REAL*16            :: vxEC(nmax),vyEC(nmax)!Euler-Cromer
      REAL*16            :: xV(nmax),yV(nmax)!Verlet
      REAL*16            :: vxV(nmax),vyV(nmax)!Verlet
      REAL*16            :: aux!auxiliar variables
      REAL*16            :: F11,F12,F13,F14!auxiliar variables
      REAL*16            :: F21,F22,F23,F24!auxiliar variables
      REAL*16            :: F31,F32,F33,F34!auxiliar variables
      REAL*16            :: F41,F42,F43,F44!auxiliar variables
      fi = -2.q0
      e   = 0.206q0
      e   = 0.6q0
      a   = 0.39q0
      x0  = (1.q0 + e) * a
      y0  = 0.q0
      vx0 = 0.q0
      vy0 = sqrt( (1.q0 - e) / x0 )
      
!      x0  = 9.q0
!      y0  = 0.q0
!      vx0 = 0.q0
!      vy0 = 1.q0/3.q0
      
      xRK(1) = x0
      yRK(1) = y0
      xEC(1) = x0
      yEC(1) = y0
      xV(1)  = x0
      yV(1)  = y0
      
      vxRK(1) = vx0
      vyRK(1) = vy0
      vxEC(1) = vx0
      vyEC(1) = vy0
      vxV(1)  = vx0
      vyV(1)  = vy0
      
!      dt = 1.q-4 * (x0**1.5q0)
      dt = 1.q-3
      imax = 10000
      if(imax.gt.nmax)then
         write(6,*)'increase nmax to',imax
         stop
      endif
      
      dtm = dt/2.q0
      do i=1,imax-1!loop for the RK4 method
         j = i + 1!auxiliar variable representing t+dt
       
         F11 = vxRK(i)
         F21 = vyRK(i)
         CALL accel( F31,F41,fi,xRK(i),vxRK(i),yRK(i),vyRK(i) )
         
         F12 = vxRK(i) + dtm*F31
         F22 = vyRK(i) + dtm*F41
         CALL accel( F32,F42,fi,xRK(i)+dtm*F11,vxRK(i)+dtm*F31,yRK(i)+dtm*F21,vyRK(i)+dtm*F41 )
         
         F13 = vxRK(i) + dtm*F32
         F23 = vyRK(i) + dtm*F42
         CALL accel( F33,F43,fi,xRK(i)+dtm*F12,vxRK(i)+dtm*F32,yRK(i)+dtm*F22,vyRK(i)+dtm*F42 )
         
         F14 = vxRK(i) + dt*F33
         F24 = vyRK(i) + dt*F43
         CALL accel( F34,F44,fi,xRK(i)+dt*F13,vxRK(i)+dt*F33,yRK(i)+dt*F23,vyRK(i)+dt*F43 )
         
         xRK(j)  = xRK(i)  + dt*(F11 + 2.q0*(F12 + F13) + F14)/6.q0
         yRK(j)  = yRK(i)  + dt*(F21 + 2.q0*(F22 + F23) + F24)/6.q0
         vxRK(j) = vxRK(i) + dt*(F31 + 2.q0*(F32 + F33) + F34)/6.q0
         vyRK(j) = vyRK(i) + dt*(F41 + 2.q0*(F42 + F43) + F44)/6.q0
      enddo
      
      do i=1,nmax-1!loop for the Euler-Cromer method
         j = i + 1
         
         CALL accel( F31,F41,fi,xEC(i),vxEC(i),yEC(i),vyEC(i) )
         vxEC(j) = vxEC(i) + dt*F31
         vyEC(j) = vyEC(i) + dt*F41
         
         xEC(j) = xEC(i) + dt*vxEC(j)
         yEC(j) = yEC(i) + dt*vyEC(j)
      enddo
      
      xV(2)  = xEC(2)
      yV(2)  = yEC(2)
!      vxV(2) = vxEC(2)
!      vyV(2) = vyEC(2)
      dt2    = dt*dt
      do i=2,nmax-1!loop for the Verlet method
         j = i + 1
         
         CALL accel( F31,F41,fi,xV(i),vxV(i),yV(i),vyV(i) )
         
         xV(j) = 2.q0*xV(i) - xV(i-1) + dt2*F31
         yV(j) = 2.q0*yV(i) - yV(i-1) + dt2*F41
!         vxV(j) = ( xV(j) - xV(i) ) / dt
!         vyV(j) = ( yV(j) - yV(i) ) / dt
      enddo
      dt2 = 2.q0*dt
      do i=2,nmax-1!loop for computing the velocities in the Verlet method
         j = i + 1
         k = i - 1
         
         vxV(i) = ( xV(j) - xV(k) ) / dt2!3 point symmetric derivative
         vyV(i) = ( yV(j) - yV(k) ) / dt2!3 point symmetric derivative
      enddo
      vxV(nmax) = ( xV(nmax) - xV(nmax-1) ) / dt!2 point backward derivative
      vyV(nmax) = ( yV(nmax) - yV(nmax-1) ) / dt!2 point backward derivative
      
      
      OPEN(unit=100,file='energy-RK.dat',status='unknown')
      OPEN(unit=200,file='energy-EC.dat',status='unknown')
      OPEN(unit=300,file='energy-V.dat',status='unknown')
      OPEN(unit=101,file='distance-RK.dat',status='unknown')
      OPEN(unit=201,file='distance-EC.dat',status='unknown')
      OPEN(unit=301,file='distance-V.dat',status='unknown')
      OPEN(unit=102,file='trajectory-difference-RK.dat',status='unknown')
      OPEN(unit=202,file='trajectory-difference-EC.dat',status='unknown')
      OPEN(unit=302,file='trajectory-difference-V.dat',status='unknown')
      OPEN(unit=103,file='trajectory-RK.dat',status='unknown')
      OPEN(unit=203,file='trajectory-EC.dat',status='unknown')
      OPEN(unit=303,file='trajectory-V.dat',status='unknown')
      do i=1,imax
         t = real(i-1,16)*dt/(2.q0*pi)
         CALL energy(aux,fi,xRK(i),vxRK(i),yRK(i),vyRK(i)) 
         write(100,*)t,aux
         CALL energy(aux,fi,xEC(i),vxEC(i),yEC(i),vyEC(i))
         write(200,*)t,aux
         CALL energy(aux,fi,xV(i),vxV(i),yV(i),vyV(i))
         write(300,*)t,aux
         
         write(101,*)t,sqrt((xRK(i)-e*a)**2 + yRK(i)**2)/x0
         write(201,*)t,sqrt((xEC(i)-e*a)**2 + yEC(i)**2)/x0
         write(301,*)t,sqrt((xV(i)-e*a)**2  + yV(i)**2)/x0
         
         write(102,*)t,abs( abs(yRK(i))/(a*sqrt(1.q0-e*e)) - sqrt(abs(1.q0-(xRK(i)/a-e)**2)) )
         write(202,*)t,abs( abs(yEC(i))/(a*sqrt(1.q0-e*e)) - sqrt(abs(1.q0-(xEC(i)/a-e)**2)) )
         write(302,*)t,abs( abs(yV(i))/(a*sqrt(1.q0-e*e))  - sqrt(abs(1.q0-(xV(i)/a-e)**2)) )
         
         write(103,*)xRK(i)-e*a,yRK(i)
         write(203,*)xEC(i)-e*a,yEC(i)
         write(303,*)xV(i)-e*a,yV(i)
      enddo
      CLOSE(100)
      CLOSE(200)
      CLOSE(300)
      CLOSE(101)
      CLOSE(201)
      CLOSE(301)
      CLOSE(102)
      CLOSE(202)
      CLOSE(302)
      CLOSE(103)
      CLOSE(203)
      CLOSE(303)
      STOP
      END PROGRAM centralforce
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!ENDO OF PROGRAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE accel(ax,ay,fi,x,vx,y,vy)
      IMPLICIT NONE
      REAL*16,INTENT(out)  :: ax,ay
      REAL*16,INTENT(in)   :: fi
      REAL*16,INTENT(in)   :: x,y
      REAL*16,INTENT(in)   :: vx,vy
      REAL*16              :: d
      d  = sqrt(x*x + y*y)
      d  = d**(fi - 1.q0)
      
      ax = - x*d
      ay = - y*d
      RETURN
      END SUBROUTINE accel
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE energy(E,fi,x,vx,y,vy)
      IMPLICIT NONE
      REAL*16,INTENT(out)  :: E
      REAL*16,INTENT(in)   :: fi
      REAL*16,INTENT(in)   :: x,y
      REAL*16,INTENT(in)   :: vx,vy
      REAL*16              :: d
      d = sqrt(x*x + y*y)
      E = d**(fi + 1.q0)/(fi + 1.q0) + 0.5q0 * (vx*vx + vy*vy)
      RETURN
      END SUBROUTINE energy
