*			PROGRAM I.11			   *

******************************************************************************
*                                  TRAJ2PT                                   *
* THIS IS THE MODIFICATION TO THE MAIN PROGRAM NEEDED TO DO TRAJECTORY PLOTS.*
* PRIMARY DIFFERENCES BETWEEN THIS AND TRAJ6PT ARE EXPLAINED ON PAGE 108 OF  *
* SHODHAN'S THESIS.                                                          *
******************************************************************************
 

      PROGRAM TRAJ1PHI

      IMPLICIT NONE
      CHARACTER*72 FNAME,FNAME1,FNAME2,FNAME3
      INTEGER I,I1,IHLF,J,MAXE,NCOUNT,NDIM,NHIT,NPAS,NPHI,NU,NTHETA
      INTEGER NXK,NYK,NZK,NTRAJ
      PARAMETER (MAXE=100)
      REAL*8 C,CON,ERRWT,QMCP
      REAL*8 STEPHI,STEPTHETA,EK,ENEK,PHIMIN,PHIMAX,EVELO,PHI(MAXE),QMC,QMCP,
     & V,VX,VY,VZ,THETA(MAXE),THETAMIN

      INCLUDE 'PASS.CMN'
      COMMON /QMC/QMC
      COMMON /NCOUNT/NCOUNT,/NU/NU,/FNAME/FNAME,/NHIT/NHIT
      COMMON /FNAME1/FNAME1
      COMMON /FNAME2/FNAME2,FNAME3
      DATA C/2.998D0/,QMCP/0.175602D0/
      PARAMETER (CON=2.540005D0,NDIM=6)
      REAL*8 AUX(16,NDIM),DERY(NDIM),PRMT(5),X0,Y(NDIM),Y0,Z0
      REAL CPUTIME,TIMER,ZTIM0

      EXTERNAL FCT,OUTP
      ZTIM0=TIMER()
      

      CALL GEOM    
    

C     OPEN(UNIT=3,FILE='TRAJ.DAT',ACCESS='SEQUENTIAL',STATUS='NEW')
C     OPEN(UNIT=4,FILE='TRAJCO.DAT',ACCESS='SEQUENTIAL',STATUS='NEW')
      OPEN(UNIT=7,FILE='TRAJSH.DAT',ACCESS='SEQUENTIAL',STATUS='NEW')
C     OPEN(UNIT=8,FILE='PASS.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*abs(X0))
      NYK=IDINT(10000.0D0*abs(Y0))
      NZK=IDINT(10000.0D0*abs(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)

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

      WRITE(6,*) 'ENTER MIN. THETA,NO. OF STEPS,DEL THETA'
      READ (5,*) THETAMIN,NTHETA,STEPTHETA

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

      WRITE (8,70) EK,V
      WRITE (8,80) THETAMIN,NTHETA,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 (2,140)

      X0 = X0 * CON
      Y0 = Y0 * CON
      Z0 = Z0 * CON

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

      DO I=1,NTHETA
	THETA(I)=THETAMIN + DFLOAT(I-1)*STEPTHETA
      END DO

      DO I=1,NPHI
       PHI(I) = PHIMIN + (I-1)*STEPHI
      END DO
 

      NPAS=0
      DO J=1,NTHETA 
       DO I=1,NPHI

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

        NCOUNT = 0        
        NU = 1
 

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

        WRITE(4,*) 
        WRITE(4,05) (I+49)
	WRITE(4,*)

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

C       WRITE(6,*) ' X: ',X0,' Y: ',Y0,' Z: ',Z0
C       WRITE(6,*) 'ENERGY (in mev) ',EK,' V ',V,' QMC ',QMC
        WRITE(6,*) ' THETA( ',J,')',THETA(J),' PHI( ',I,')',PHI(I)
C       WRITE(6,*) ' VX ',VX,' VY ',VY,' VZ ',VZ

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(J)
         PAS(NPAS,2) = PHI(I)
	END IF

        WRITE(7,40)
C       WRITE(3,40)
        WRITE(4,40)
        WRITE(6,*) 
C       WRITE(6,*) ' Y ',(Y(I1),I1=1,NDIM)
        WRITE(6,*) ' IHLF NO. OF BISECTIONS OF STEP: ',IHLF
        WRITE(6,*) 'TOTAL NO. OF POINTS IN THE TRAJECTORY: ',NU-1
       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,'INITIAL THETA',F10.2,2X,'NO. OF STEPS',
     &   I5,2X,'DEL THETA',F10.2)
 90   FORMAT (1X,'INITIAL PHI',F10.2,2X,'NO. OF STEPS',
     &   I5,2X,'DEL PHI',F10.2)
 95   FORMAT (1X,'STARTING POSITION(INCHES)',2X,'X(+02)',D14.7,
     &   2X,'Y(+02)',D14.7,2X,'Z(+02)',D14.7)
100   FORMAT (1X,'STARTING POSITION(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')
130   FORMAT (1X,'POLAR & AZIM ANGLES AT DETECTOR')
140   FORMAT (1X,'POLAR & AZIM ANGLES AT THE APERATURE')

      CLOSE(1)
      CLOSE(2)
      CLOSE(7)      
      CLOSE(4)
      CLOSE(8)
      STOP
      END

Return to the Table of Contents for Buckley's Thesis

Return to HISCALE List of Appendices

Return to Ulysses HISCALE Data Analysis Handbook Table of Contents