C234567890
C     Program to calculate the motion of an atom
C     that 'scatters' off a surface of atoms, interacting
C     by the Lennard-Jones potential
C
C     The units are chosen so that energy units E=1, mass units
C     M=1, and distance units s=1.
C     Inputs:
C     X0, Y0: intial starting position
C     VX0, VY0: initial velocity
C     TS: time step scaling factor. Time units are sqrt[M s**2 E]. 
C         The time step size DT is calculated by taking the distance
C         to the scattering atoms, say R, and setting
C         DT = (R/V0) * 1/TS, where V0 is the initial velocity.
C         
C     Outputs:
C      A file is created, scatter.dat, which contains the
C      x and y coordinates of the atom along its trajectory.
C      Some diagnostics are printed on the console during the program.
C
C234567
C      Declarations
      DOUBLE PRECISION X0, Y0, VX0, VY0, V0, DT, TS, R, RR, M, S, E
      DOUBLE PRECISION XN, YN, VXN, VYN, XN1, YN1, VXN1, VYN1
      INTEGER*4 II,JJ,KK,LL
C
C     INPUT
      WRITE(*,*) 'Input X0, Y0 on one line'
      READ(*,*) X0, Y0
      WRITE(*,*) 'Input VX0, VY0 on one line'
      READ(*,*) VX0, VY0
      WRITE(*,*) 'Input time scaling factor, TS'
      READ(*,*) TS
C
C     SET UP INITIAL CONDITIONS
      M=1
      R=SQRT(X0*X0+Y0*Y0)
      V0=SQRT(VX0*VX0+VY0*VY0)
      E=0.5*M*V0*V0
      DT=(R/V0)*TS
      VXN=VX0
      VYN=VY0
      XN=X0
      YN=Y0
C    
      WRITE(*,*) 'R, V0, E, DT'
      WRITE(*,*) R,V0,E,DT
      WRITE(*,*) 'XN, YN, VXN, VYN'
      WRITE(*,*) XN,YN,VXN,VYN
      PAUSE 'Hit CR to continue'
C
      OPEN(UNIT=11,FILE='scatter.dat',STATUS='NEW')
C
C     Calculate the trajectory, step by step, until
C     the maximum number of steps is exceeded or the
C     the maximum range is exceeded.
C     KK is maximum number of steps
      KK=1000
      WRITE(11,*) XN,YN
      DO 100 II=1,KK
C     REPORT ON HOW IT'S GOING
      WRITE(*,15) 'LOOP NO. ',II
15    FORMAT(X,A10,I8)
      CALL NEXTSTEP(XN,YN,VXN,VYN,XN1,YN1,VXN1,VYN1,DT)
      WRITE(11,20) XN1,YN1,VXN1,VYN1
20    FORMAT(' ',F15.6,F15.6,F15.6,F15.6)
C    UPDATE POSITIONS
      XN=XN1
      YN=YN1
      VXN=VXN1
      VYN=VYN1
      RR=SQRT(XN1*XN1 + YN1*YN1)
C   THIS IS FORTRAN 'SPAGHETTI' CODE! WILL GENERATE A WARNING!
      IF(RR.GE.10*R) II=KK
100   CONTINUE
C     CLEAN UP THINGS
      CLOSE(11)
      STOP
      END
C
C  SUBROUTINES AND FUNCTIONS
C
      SUBROUTINE NEXTSTEP(XN,YN,VXN,VYN,XN1,YN1,VXN1,VYN1,DT)
C
      DOUBLE PRECISION XN,YN,VXN,VYN,XN1,YN1,VXN1,VYN1,DT
      DOUBLE PRECISION ACCELX, ACCELY, AVACCELX, AVACCELY
      DOUBLE PRECISION ACCELXN, ACCELYN
C     VERLET ALGORITHM:
      ACCELXN=ACCELX(XN,YN)
      ACCELYN=ACCELY(XN,YN)
      XN1=XN + VXN*DT + 0.5 * ACCELXN * DT**2
      YN1=YN + VYN*DT + 0.5 * ACCELYN * DT**2
      AVACCELX = 0.5 * (ACCELXN + ACCELX(XN1,YN1))
      AVACCELY = 0.5 * (ACCELYN + ACCELY(XN1,YN1))
      VXN1=VXN + AVACCELX*DT
      VYN1=VYN + AVACCELY*DT
C
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION ACCELX(X,Y)
C
      DOUBLE PRECISION X, Y
C     USE FOR DEBUGGING PURPOSES:
      ACCELX = 0.0
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION ACCELY(X,Y)
C
      DOUBLE PRECISION X, Y
      ACCELY=0.0
      RETURN
      END

