*		         PROGRAM I.1			  	  *

******************************************************************************
*                               TRAJ6PT                                      *
* THIS IS THE MAIN TRAJECTORY PROGRAM.  ALL VARIABLES AND PARAMETERS ARE     *
* DEFINED IN THE COMMENTS FOUND ON PAGE 65 OF SHODHAN'S THESIS.              *
******************************************************************************
 

      PROGRAM TRAJ1PHI
      IMPLICIT NONE
      CHARACTER*72 FNAME,FNAME1,FNAME2,FNAME3,FNAME4
      INTEGER I,I1,IHLF,J,MAXE,NCOUNT,NDIM,NHIT,NPHI,NTHETA,NU
      INTEGER NXK,NYK,NZK,MAXSURF
      PARAMETER (MAXE=100)
      REAL*8 C,CON,ERRWT
      PARAMETER (CON=2.540005D0,NDIM=6)
      REAL*8 STEPHI,STEPTHETA,EK,ENEK,PHIMIN,EVELO,
     & PHI,QMC,QMCP,
     & THETA,THETAMAX,THETAMID,V,VX,VY,VZ
      REAL*8 AUX(16,NDIM),DERY(NDIM),PRMT(5),X0,Y(NDIM),Y0,Z0
      REAL CPUTIME,TIMER,ZTIM0

      INCLUDE 'PASS5.CMN'
      COMMON /QMC/QMC
      COMMON /NCOUNT/NCOUNT,/NU/NU,/NHIT/NHIT
C     COMMON /FNAME/FNAME
      COMMON /FNAME4/FNAME4
      COMMON /FNAME1/FNAME1
      COMMON /FNAME2/FNAME2,FNAME3
      DATA C/2.998D0/,QMCP/0.175602D0/
      EXTERNAL FCT,OUTP
      ZTIM0=TIMER()    
      CALL GEOM    
    

C     OPEN(UNIT=7,FILE='TRAJSH.DAT',ACCESS='SEQUENTIAL',STATUS='NEW')

      WRITE(6,*) 'ENTER THE INITIAL COORDINATES OF THE e IN INCHES'
      READ(5,*) X0,Y0,Z0

C     TO ENTER QUANTITIES FOR PRMT
      WRITE(6,*) 'ENTER LOWER BOUND ON TIME (t=0)'
      READ(5,*) PRMT(1)
      WRITE(6,*) 'ENTER THE UPPER BOUND ON TIME'
      READ(5,*) PRMT(2)
      WRITE(6,*) 'ENTER THE TIME STEP'
      READ(5,*) PRMT(3)
      WRITE(6,*) 'ENTER THE ERROR BOUND'
      READ(5,*) PRMT(4)
   

      WRITE(6,*) 'ENTER THE ENERGY OF THE PARTICLE IN MEV'
      READ(5,*) EK

      ENEK=IDINT(1000.0D0*EK)  
      NXK=IDINT(-10000.0D0*X0)
      NYK=IDINT(-10000.0D0*Y0)
      NZK=IDINT(10000.0D0*Z0)

      ENCODE (17,60,FNAME1)NXK,NYK,NZK,IDINT(ENEK)
      OPEN(UNIT=8,STATUS='NEW',ACCESS='SEQUENTIAL',FILE=FNAME1)

      ENCODE (19,110,FNAME2)NXK,NYK,NZK,IDINT(ENEK)
      OPEN(UNIT=1,STATUS='NEW',ACCESS='SEQUENTIAL',FILE=FNAME2)

      ENCODE (20,120,FNAME3)NXK,NYK,NZK,IDINT(ENEK)
      OPEN(UNIT=2,STATUS='NEW',ACCESS='SEQUENTIAL',FILE=FNAME3)

      ENCODE (23,125,FNAME4)NXK,NYK,NZK,IDINT(ENEK)
      OPEN(UNIT=3,STATUS='NEW',ACCESS='SEQUENTIAL',FILE=FNAME4)

      V = EVELO(EK)          !COMPUTE V FOR THIS ENERGY OF e
      QMC = -QMCP*DSQRT(1 - (V/C)**2)  

      WRITE (6,*) 'ENTER MID. THETA,THETAMAX,DEL THETA'
      READ (5,*) THETAMID,THETAMAX,STEPTHETA

      WRITE (6,*) 'ENTER MIN. PHI,NO. OF STEPS,DEL PHI'
      READ (5,*) PHIMIN,NPHI,STEPHI

      WRITE (8,70) EK,V
      WRITE (8,80) THETAMID,THETAMAX,STEPTHETA
      WRITE (8,90) PHIMIN,NPHI,STEPHI
      WRITE (8,20) PRMT(1),PRMT(2),PRMT(3),PRMT(4)
      WRITE (8,*)
      WRITE (8,95) X0,Y0,Z0

      WRITE (1,130)
	WRITE (1,*) (1) 
      WRITE (2,140)
	WRITE (2,*) (1)
      WRITE (3,150)

C     TO CONVERT FROM INCHES TO CMS. 
      X0 = X0 * CON
      Y0 = Y0 * CON
      Z0 = Z0 * CON

      WRITE (8,100) X0,Y0,Z0
      WRITE (8,*)
 

      NPAS=0 

      DO I=1,NPHI
 

        PHI=PHIMIN + DFLOAT(I-1)*STEPHI
       

        DO J=-1,1,2
          IF (J .EQ. -1) THEN
           THETA=THETAMID
          ELSE
           THETA=THETAMID + DFLOAT(J)*STEPTHETA
          END IF
          NHIT=2
          DO WHILE (NHIT .EQ. 2 .AND. THETA .LT. THETAMAX)

C            ENCODE (23,50,FNAME)NXK,NYK,NZK,IDINT(ENEK),
C     &                          IDINT(THETA),IDINT(PHI)
C            OPEN(UNIT=4,STATUS='NEW',ACCESS='SEQUENTIAL',FILE=FNAME)

            NCOUNT = 0        
            NU = 1

C           WRITE(7,05) I
C           WRITE(7,10) EK,V,THETA(J),PHI(I)
C           WRITE(7,20) PRMT(1),PRMT(2),PRMT(3),PRMT(4)
C           WRITE(7,30)
C           WRITE(6,*) J

C           TO COMPUTE THE VELOCITY PROJECTIONS Vx,Vy,Vz
            CALL VELOPROJ(V,VX,VY,VZ,THETA,PHI)
    

            WRITE(6,*) ' THETA ',THETA,' PHI ',PHI

C           TO INITIALISE THE INITIAL VALUES 
            Y(1) = X0
            Y(2) = Y0
            Y(3) = Z0
            Y(4) = VX
            Y(5) = VY
            Y(6) = VZ

            ERRWT = 1.0D0/6.0D0
            DO I1=1,NDIM
             DERY(I1) = ERRWT
            END DO
            CALL DHPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)

            IF (NHIT .EQ. 2) THEN
             PAS(NPAS,1) = THETA
             PAS(NPAS,2) = PHI
            END IF

            WRITE(6,*) 
            WRITE(6,*) ' IHLF NO. OF BISECTIONS OF STEP: ',IHLF
            WRITE(6,*) 'TOTAL NO. OF POINTS IN THE TRAJECTORY: ',NU-1
	    

	    IF ((THETA.LT.140).AND.(THETA.GT.40)) NHIT = 2
	    THETA = THETA+DFLOAT(J)*STEPTHETA

          END DO
        END DO 
      END DO

      CALL PASSOUTPUT

      WRITE(6,*) 'NPAS:',NPAS
      CPUTIME=TIMER()-ZTIM0
      WRITE(6,*) 'C.P.U. TIME: ',CPUTIME
      WRITE(1,40)
      WRITE(2,40)

 05   FORMAT(1X,I3)
 10   FORMAT(1X,'ENERGY(in mev)',D10.3,2X,'VELOCITY(*10+10)',F12.5,2X,
     & 'THETA(in deg.)',F10.3,2X,'PHI(in deg.)',F10.3)
 20   FORMAT(1X,'INITIAL TIME(*10-08)',F12.6,2X,'FINAL TIME(*10-08)',
     & F12.6,2X,'INITIAL STEP(*10-08)',F14.8,1X,'ERROR BOUND',F19.12)
 30   FORMAT(4X,'T(-08)',10X,'X(+02)',10X,'Y(+02)',10X,'Z(+02)',10X,
     & 'VX(+10)',10X,'VY(+10)',10X,'VZ(+10)',10X,'(V+10)')
 40   FORMAT(X,'---*---*---*--- END OF ENERGY ---*---*---*---')
 50   FORMAT (I3,I3,I4,'.E',I5,I3,I3)
 60   FORMAT (I3,I3,I4,'.E',I5)
 70   FORMAT (1X,'ENERGY(in mev)',D10.3,2X,'VELOCITY(*10+10)',F12.5)
 80   FORMAT (1X,'MIDDLE THETA(in deg.)',F10.2,2X,'MAX.THETA',
     & F10.2,2X,'DEL THETA',F10.2)
 90   FORMAT (1X,'INITIAL PHI  (in deg.)',F10.2,2X,'NO. OF STEPS',
     & I5,2X,'DEL PHI  ',F10.2)
 95   FORMAT (1X,'STARTING POSITION(in in.)',2X,'X(+02)',D14.7,
     & 2X,'Y(+02)',D14.7,2X,'Z(+02)',D14.7)
100   FORMAT (1X,'STARTING POSITION(in cm.)',2X,'X(+02)',D14.7,
     & 2X,'Y(+02)',D14.7,2X,'Z(+02)',D14.7)
110   FORMAT (I3,I3,I4,'.E',I5,'IN')
120   FORMAT (I3,I3,I4,'.E',I5,'OUT')
125   FORMAT (I3,I3,I4,'.E',I5,'OUTVEL')
130   FORMAT(1X,'POLAR & AZIMUTHAL ANGLES AT THE DETECTOR')
140   FORMAT(1X,'POLAR & AZIMUTHAL ANGLES AT THE APERTURE')
150   FORMAT(1X,'POLAR & AZIMUTHAL ANGLES AT THE APERTURE FROM VEL.')
      CLOSE(1)
      CLOSE(2)
      CLOSE(3)
      END

C-----------------------------------------------------------------------------
C-----------------------------------------------------------------------------