C===============================================================
C--- CPS500 Homework#4
C--- Higher order equations by RUNGE-KUTTA method
C---   Programmed by Sung-Hoon Ko (shko@npac.syr.edu)
C--------------------------------------------------------------
C                          ..          ..
C--- Solve the equations m*X1=-F and m*X2=+F
C---   where F=(-6*b)/r**7+(12*c)/r**13
C---   where r=ABS(X1-X2)
C--- Initial Conditions;
C---   X1=0, X2=R, R=4
C---   m=b=c=1, V=0.1, DT=0.001
C--- Carry the integration up to the point where
C--- the particles have exactly the same positions
C--- but opposite velocities to initial conditions.

      PROGRAM HW4_RUNGE_KUTTA
      implicit none
      INTEGER N
      REAL DT
      parameter (DT=0.001)  !Delta t
      parameter (N=100000)
      REAL    X1(N),X2(N),V1(N),V2(N),DT
      INTEGER i               ! Array index
      EXTERNAL dp_initialize  ! External routines to initialize
      EXTERNAL dp_draw_graph  !  and draw 

C--- Initial conditions
      X1(1)=0.     
      X2(1)=4.
      V1(1)=0.1
      V2(1)=-V1(1)
      
      i=1
  10  CALL RUNGE(X1,X2,V1,V2,DT,i)
      i=i+1
C--- Check to terminate
      IF( (X1(i).LE.0.).AND.(X2(i).GE.4.)
     *    .AND.(V1(i).LE.-0.1) ) then
          goto 999
      ELSE 
          goto 10
      ENDIF    

  999 write(6,*) 'XX1=',X1(i),'X2=',X2(i),'V1=',V1(i),'V2=',V2(i)
      write(6,*) 'I=',i 
     
      CALL dp_initialize      
      CALL dp_draw_graph(X1,V1,X2,V2,i)

      PAUSE
      STOP
      END

C-------------------------------------------------------
      REAL FUNCTION FN(x1,x2)
      REAL x1,x2
      REAL r

      r=(ABS(x1-x2))**7
      FN=-6/r+12/(r**6)
      RETURN
      END


C---------------------------------------------------------
C---  DX1=V1, DX2=V2
C---  DV1=-F, DV2=+F

      SUBROUTINE  RUNGE(x1,x2,v1,v2,dt,j)
      REAL x1(*),x2(*),v1(*),v2(*),dt
      integer j
      
      REAL k1,k2,k3,k4
      REAL m1,m2,m3,m4
                   
      k1=v1(j)
      m1=FN(x1(j),x2(j))
      k2=v1(j)+m1*dt/2.
      m2=FN(x1(j)+k1*dt/2.,x2(j)+k1*dt/2.)
      k3=v1(j)+m2*dt/2.
      m3=FN(x1(j)+k2*dt/2.,x2(j)+k2*dt/2.)
      k4=v1(j)+m3*dt
      m4=FN(x1(j)+k3*dt,x2(j)+k3*dt)

C--- Calculate next x1 and v1 by Runge-Kutta method.
      x1(j+1)=x1(j)+dt*(k1+2*k2+2*k3+k4)/6.
      x2(j+1)=4.-x1(j)

      v1(j+1)=v1(j)+dt*(m1+2*m2+2*m3+m4)/6.
      v2(j+1)=-v1(j)

      RETURN
      END
