*******************************************************************************
*                                                                             *
*			PROGRAM I.1			    *
*                                                                             *
*******************************************************************************
*                                                                             *
*                               TRAJ4P.FOR                                    *
* 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,wscatter,SCATTER,nSCATTER,startnsurf
      	PARAMETER (MAXE=180)
      	REAL*8 C,CON,ERRWT
      	PARAMETER (CON=1.0D0,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
        logical firstcheck,secondcheck
	Integer SEED1,SEED2,SEED3,Seed4
	COMMON /SEED1/SEED1,/SEED2/SEED2,/SEED3/SEED3,/Seed4/Seed4

      	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
	common /wscatter/wscatter,/scatter/scatter,/Nscatter/Nscatter
	common /theta/theta,/phi/phi,/EK/EK,/startnsurf/startnsurf
	common /firstcheck/firstcheck,/secondcheck/secondcheck
	data scatter/0/
	DATA SEED1/1234441519/
	DATA SEED2/278161611/
	DATA SEED3/467321899/
	Data Seed4/598516711/
      	DATA C/2.998D0/,QMCP/0.175602D0/
      	EXTERNAL FCT,OUTP
C    	ZTIM0=TIMER()

	write(6,*) 'ENTER THE CHOICE OF MODEL OF THE SCATTER'
	write(6,*) 'NON SCATTER----------->0'
	write(6,*) 'SPECULAR SCATTER------>1'
	write(6,*) 'DIFFUSE SCATTER------->2'
	read(5,*) wscatter    

	if (wscatter.eq.2) then
	    write(6,*) 'ENTER NSCATTER'
	    read(5,*) Nscatter
	end if

      	CALL GEOM    
    

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

      	WRITE(6,*) 'ENTER THE INITIAL COORDINATES OF THE e IN CM/100s'
      	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*dabs(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,*)
C      	WRITE (8,95) X0,Y0,Z0

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

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

	write(6,*)X0,Y0,Z0
      	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)
    

C            	    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
		    startnsurf=1
		    firstcheck=.false.
		    secondcheck=.false.
C	            write(6,*)' THETA ',THETA,' PHI ',PHI
	    	    CALL DHPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)

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

C            	    WRITE(6,*) 
C            	    WRITE(6,*) ' IHLF NO. OF BISECTIONS OF STEP: ',IHLF
C            	    WRITE(6,*) 'TOTAL NO. OF POINTS IN THE TRAJECTORY: ',NU-1	    
	    	    IF ((THETA.LT.Thetamax).AND.(THETA.GT.60)) NHIT = 2
	    		THETA = THETA+DFLOAT(J)*STEPTHETA
          	END DO
            END DO 
      	END DO

	WRITE (1,*) NPAS
	WRITE (2,*) NPAS	

      	CALL PASSOUTPUT

      	WRITE(6,*) 'NPAS:',NPAS
C     	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)     
      	close(8)
      	stop
      	END           

 

 

Return to Appendix 10.11 List of Programs