*
*                                                                            *
******************************************************************************
*                                                                            *
*    THE SUBROUTINES THAT ARE USED TO FIND WHETHER THE LINE SEGMENT HAS HIT  *
*  ANY OF THE SURFACES AND LOST OR NOT. Check Buckley's Thesis               *
******************************************************************************

      SUBROUTINE CHECKHIT(HIT,NSURF)

      IMPLICIT NONE
      INTEGER I,J,HIT,MAXSURF,NSURF,NTY,NVERT
      PARAMETER (MAXSURF=69)
      REAL*8 A,B,C,COEFF(MAXSURF,4),D,VERT(10,3)
      logical done
      COMMON /COEFF/COEFF,/VERT/VERT,NVERT,/NTY/NTY

      NSURF=1
      HIT=0
      done=.false.
      DO WHILE((NSURF.LE.NTY).AND.(.not.done))
       CALL FINDVERT(VERT,NVERT,NSURF)
D       WRITE(6,*) 'NO. OF VERTICES ARE: ',NVERT
D       DO I=1,NVERT
D        WRITE(6,*) (VERT(I,J),J=1,3)
D       END DO
       A=COEFF(NSURF,1)
       B=COEFF(NSURF,2)
       C=COEFF(NSURF,3)
       D=COEFF(NSURF,4)
D      WRITE(6,*)(COEFF(NSURF,1),COEFF(NSURF,2),COEFF(NSURF,3))

      CALL SCANPOLY(done,A,B,C,D,HIT,NSURF)	
       NSURF=NSURF+1
      END DO
      RETURN
      END
C------------------------------------------------------------------
      SUBROUTINE FINDVERT(VERT,NVERT,NSURF)

      INTEGER I1,I2,J,MAXSURF,MAXCOO,NVERT
      PARAMETER (MAXSURF=69,MAXCOO=69)
      REAL*8 CHAN(MAXSURF,MAXCOO),VERT(10,3)  
      COMMON /CHAN/CHAN
  

      NVERT=IDINT(CHAN(NSURF,1))
      I2=2
      DO I1=1,NVERT
       DO J=1,3
        VERT(I1,J)=CHAN(NSURF,I2)
        I2=I2+1
       END DO
      END DO
      RETURN
      END

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

      SUBROUTINE SCANPOLY(done,A,B,C,D,HIT,I)

 
      IMPLICIT NONE
      REAL*8 A,B,C,D,SECT(3),TLN(2,3),VERT(10,3),c1(4)
      INTEGER FLAG,HIT,I,NCOUNT,NGAMAP,NLOWAP,NERR,NVERT,NOPAP,
     & scatter,wscatter,ape 
      logical done
	common /scatter/scatter,/wscatter/wscatter
	common /sect/sect,/c1/c1
      PARAMETER (NOPAP=70,NGAMAP=70,NLOWAP=70,ape=7)
C HIT : WHOSE VALUE INDICATES WHETHER TRAJECTORY CALCULATION IS TO BE CONTINUED
C      OR NOT.
C     = 0  PARTICLE IS NOT LOST-SO CONTINUE WITH THE TRAJECTORY CALCULATION.
C     = 1  PARTICLE IS LOST-DO NOT CONTINUE WITH THE TRAJECTORY CALCULATION.
C     = 2  PARTICLE ESCAPES FROM THE SENSOR ASSEMBLY WITHOUT HITTING 
C          ANY OF THE OTHER SURFACES. NO NEED TO CONTINUE WITH THE 
C          TRAJECTORY CALCULATION - SO STOP CALCULATION.
          

	c1(1)=A
	c1(2)=B
	c1(3)=C
	c1(4)=D

      	CALL INTERSECT(SECT,A,B,C,D,NERR)

C                   NERR = 0 :NO INTERSECTION BETWEEN THE LINE & THE PLANE

      	IF (NERR .EQ. 0) THEN
       	    HIT=0
D           WRITE(6,*) 'LINE SEGMENT PARALLEL TO THE PLANE'
       	RETURN
      	END IF
   

C        SINCE AT THS POINT INTERSECTION IS POSSIBLE,FIND WHETHER THE POINT 
C     		BELONGS TO THE LINE SEGMENT OR NOT.

D      	WRITE(6,*) 'INTERSECTION POINT :',(SECT(NCOUNT),NCOUNT=1,3)
D	WRITE(6,*) I
D      	WRITE(6,*) 'TEST WHETHER IT BELONGS TO THE LINE-SEGMENT OR NOT'
      

	CALL BETWEEN(SECT,NERR)

C       Between is a subroutine to test whether segment of tln(1,i) to tln(2,i)
C	    has the intercept with any planes
C                  NERR = 0 :POINT DOES NOT BELONG TO THE LINE-SEGMENT
C                            SO IT HAS NOT REACHED THE PLANE YET-HANGING THERE

      	IF (NERR .EQ. 0) THEN
       	    HIT=0
D       WRITE(6,*) 'POINT DOES NOT BELONG TO THE LINE SEGMENT'
       	

	RETURN
      	

	END IF      

C
C        POINT BELONGS TO THE LINE-SEGMENT-TEST WHETHER THIS 
C        INTERSECTION POINT LIES ON THE EDGES OF THE POLYGON OR NOT
C

D      	WRITE(6,*) 'TEST WHETHER IT LIES ON THE EDGES OF THE 
D    &  POLYGON OR NOT'
      	

	CALL HITEDGE(SECT,NERR)

C                 NERR = 0 :IT DOES NOT LIE ON ANY OF THE 
C                           EDGES OF THE POLYGON
C                      = 1 :PARTICLE LOST-SINCE IT LIES ON ONE OF 
C                           THE EDGES OF THE POLYGON.
                        

      	IF (NERR .EQ. 1) THEN
            if ((wscatter.eq.0) .or.(I.eq.nopap)) then
		hit=1
		scatter=1
	    else           
	    	scatter=1
	    	hit=0
	    end if
            done=.true.

D       WRITE(6,*) 'IT LIES ON THE EDGES OF THE POLYGON'
        

	RETURN
      	END IF

C     	TEST WHETHER THE PARTICLE LIES INSIDE THE POLYGON OR NOT
D      	WRITE(6,*) 'TEST WHETHER IT LIES INSIDE THE POLYGON OR NOT'

      	CALL INOUT(SECT,FLAG)

C                 FLAG = 0 :PARTCLE LIES OUTSIDE THE POLYGON
C                      = 1 :IT LIES INSIDE THE POLYGON
 

      	IF (FLAG .EQ. 1) THEN 
	   		!done with the loop of checkhit
       	    IF ((I .EQ. NLOWAP) .OR. (I .EQ. NGAMAP))THEN 

                      !LIES INSIDE THE LOWER APERTURE OR GAMMA OPENING
c		hit=1
c		scatter=1
            	HIT=0
D           	WRITE(6,*) 'IT LIES INSIDE THE LOWER POLYGON'
	

	    RETURN
       	

	    ELSE
            	IF (I .le.ape) THEN      !INSIDE THE OUTER APERTURE
		    done=.true.
         	    HIT=2
C	            write(6,*) wscatter,scatter,sect(1),sect(2),sect(3)
D          	    WRITE(6,*) 'IT LIES INSIDE THE OUTER OPENING APERUTRE'
                RETURN
                

		END IF
       	

	    END IF
            done=.true.
	    if ((wscatter.eq.0).or.(i.eq.nopap)) then
	    	hit=1
		scatter=1
	    else
	    	scatter=1
	    	hit=0
	    end if
                         !INSIDE THE OTHER POLYGONS  

D        WRITE(6,*) 'IT LIES INSIDE THE OTHER POLYGONS'
       	RETURN
      	END IF
   

      	IF (FLAG .NE. 0) THEN

D           WRITE(6,*) 'ERROR IN INOUT ROUTINE-CHECK'
            HIT=10
	    done=.true.       	
	RETURN      	
	END IF

      	IF ((I .EQ. NLOWAP).OR.(I .EQ. NOPAP).OR.(I.EQ.NGAMAP)) THEN 

                         !LIES OUTSIDE THE APERTURE SURFACES
	    done=.true.
       	    HIT=1    
            scatter=1                                 

D           WRITE(6,*) 'LIES OUTSIDE THE APERTURE SURFACES'

      	ELSE
            HIT=0                          !OUTSIDE THE OTHER PLANE     
D           WRITE(6,*) 'LIES OUTSIDE OTHER PLANES-SO NOT LOST YET'

      	END IF                             !POLYGONS-MAY NOT BE LOST YET
        

      	RETURN
      	END
C ***********************************************************************
      SUBROUTINE CHECKHIT1(HIT,NSURF)

      IMPLICIT NONE
      INTEGER I,J,HIT,MAXSURF,NSURF,NTY,NVERT
      PARAMETER (MAXSURF=69)
      REAL*8 A,B,C,COEFF(MAXSURF,4),D,VERT(10,3)
      COMMON /COEFF/COEFF,/VERT/VERT,NVERT,/NTY/NTY

      NSURF=1
      HIT=0
      DO WHILE((NSURF.LE.NTY).AND.(HIT.EQ.0))
       CALL FINDVERT(VERT,NVERT,NSURF)
D       WRITE(6,*) 'NO. OF VERTICES ARE: ',NVERT
D       DO I=1,NVERT
D        WRITE(6,*) (VERT(I,J),J=1,3)
D       END DO
       A=COEFF(NSURF,1)
       B=COEFF(NSURF,2)
       C=COEFF(NSURF,3)
       D=COEFF(NSURF,4)
D      WRITE(6,*)(COEFF(NSURF,1),COEFF(NSURF,2),COEFF(NSURF,3))
       CALL SCANPOLY1(A,B,C,D,HIT,NSURF)
       NSURF=NSURF+1
      END DO
      RETURN
      END
C------------------------------------------------------------------
      SUBROUTINE SCANPOLY1(A,B,C,D,HIT,I)
 

      IMPLICIT NONE
      REAL*8 A,B,C,D,SECT(3),TLN(2,3),VERT(10,3)
      INTEGER FLAG,HIT,I,NCOUNT,NGAMAP,NLOWAP,NERR,NVERT,NOPAP,ape 
      PARAMETER (NOPAP=70,NGAMAP=70,NLOWAP=70,ape=7)
C HIT : WHOSE VALUE INDICATES WHETHER TRAJECTORY CALCULATION IS TO BE CONTINUED
C      OR NOT.
C     = 0  PARTICLE IS NOT LOST-SO CONTINUE WITH THE TRAJECTORY CALCULATION.
C     = 1  PARTICLE IS LOST-DO NOT CONTINUE WITH THE TRAJECTORY CALCULATION.
C     = 2  PARTICLE ESCAPES FROM THE SENSOR ASSEMBLY WITHOUT HITTING 
C          ANY OF THE OTHER SURFACES. NO NEED TO CONTINUE WITH THE 
C          TRAJECTORY CALCULATION - SO STOP CALCULATION.
          

      CALL INTERSECT(SECT,A,B,C,D,NERR)
C                   NERR = 0 :NO INTERSECTION BETWEEN THE LINE & THE PLANE
      IF (NERR .EQ. 0) THEN
       HIT=0
D       WRITE(6,*) 'LINE SEGMENT PARALLEL TO THE PLANE'
       RETURN
      END IF
   

C        SINCE AT THS POINT INTERSECTION IS POSSIBLE,FIND WHETHER THE POINT 
C     BELONGS TO THE LINE SEGMENT OR NOT.
D      WRITE(6,*) 'INTERSECTION POINT :',(SECT(NCOUNT),NCOUNT=1,3)
D	WRITE(6,*) I
D      WRITE(6,*) 'TEST WHETHER IT BELONGS TO THE LINE-SEGMENT OR NOT'
      CALL BETWEEN(SECT,NERR)
C                  NERR = 0 :POINT DOES NOT BELONG TO THE LINE-SEGMENT
C                            SO IT HAS NOT REACHED THE PLANE YET-HANGING THERE

      IF (NERR .EQ. 0) THEN
       HIT=0
D       WRITE(6,*) 'POINT DOES NOT BELONG TO THE LINE SEGMENT'
       RETURN
      END IF      

C        POINT BELONGS TO THE LINE-SEGMENT-TEST WHETHER THIS 
C        INTERSECTION POINT 
C     LIES ON THE EDGES OF THE POLYGON OR NOT

D      WRITE(6,*) 'TEST WHETHER IT LIES ON THE EDGES OF THE 
D    &  POLYGON OR NOT'
      CALL HITEDGE(SECT,NERR)
C                 NERR = 0 :IT DOES NOT LIE ON ANY OF THE 
C                           EDGES OF THE POLYGON
C                      = 1 :PARTICLE LOST-SINCE IT LIES ON ONE OF 
C                           THE EDGES OF THE POLYGON.
                        

      IF (NERR .EQ. 1) THEN
        HIT=1
D        WRITE(6,*) 'IT LIES ON THE EDGES OF THE POLYGON'
        RETURN
      END IF

C     TEST WHETHER THE PARTICLE LIES INSIDE THE POLYGON OR NOT

D      WRITE(6,*) 'TEST WHETHER IT LIES INSIDE THE POLYGON OR NOT'
      CALL INOUT(SECT,FLAG)
C                 FLAG = 0 :PARTCLE LIES OUTSIDE THE POLYGON
C                      = 1 :IT LIES INSIDE THE POLYGON
 

      IF (FLAG .EQ. 1) THEN
       IF ((I .EQ. NLOWAP) .OR. (I .EQ. NGAMAP))THEN 
                      !LIES INSIDE THE LOWER APERTURE OR GAMMA OPENING
        HIT=0
D         WRITE(6,*) 'IT LIES INSIDE THE LOWER POLYGON'
        RETURN
       ELSE
        IF (I .le. ape) THEN      !INSIDE THE OUTER APERTURE
         HIT=2
D          WRITE(6,*) 'IT LIES INSIDE THE OUTER OPENING APERUTRE'
         RETURN
        END IF
       END IF
       HIT=1                     !INSIDE THE OTHER POLYGONS         
D        WRITE(6,*) 'IT LIES INSIDE THE OTHER POLYGONS'
       RETURN
      END IF
   

      IF (FLAG .NE. 0) THEN
D       WRITE(6,*) 'ERROR IN INOUT ROUTINE-CHECK'
       HIT=10
       RETURN
      END IF

      IF ((I .EQ. NLOWAP).OR.(I .EQ. NOPAP).OR.(I.EQ.NGAMAP)) THEN 
                         !LIES OUTSIDE THE APERTURE SURFACES
       HIT=1                                     
D        WRITE(6,*) 'LIES OUTSIDE THE APERTURE SURFACES'
      ELSE
       HIT=0                                      !OUTSIDE THE OTHER PLANE     
D        WRITE(6,*) 'LIES OUTSIDE OTHER PLANES-SO NOT LOST YET'
      END IF                                      !POLYGONS-MAY NOT BE LOST YET
   

      
      RETURN
      END
C ***********************************************************************

      SUBROUTINE INTERSECT(SECT,A,B,C,D,NERR)

      IMPLICIT NONE
      REAL*8 A,B,C,D,DX,DY,DZ,DET,RATIO
      REAL*8 TLN(2,3),SECT(3),NUM
      INTEGER NERR,i
      COMMON /TLN/TLN     
      NERR = 1
      DX=TLN(2,1)-TLN(1,1)
      DY=TLN(2,2)-TLN(1,2)
      DZ=TLN(2,3)-TLN(1,3)
D	WRITE(6,*) A,B,C,D
D	WRITE(6,*) (DX,DY,DZ)
      DET=A*DX+B*DY+C*DZ
      NUM=-(A*TLN(1,1)+B*TLN(1,2)+C*TLN(1,3)+D)

D	WRITE(6,*) (DET,NUM)
      IF (DET.EQ. 0.0) THEN
        NERR=0
        RETURN
      END IF
      RATIO=NUM/DET
D	WRITE(6,*) RATIO
C
      SECT(1)=DX*RATIO+TLN(1,1)
      SECT(2)=DY*RATIO+TLN(1,2)
      SECT(3)=DZ*RATIO+TLN(1,3)
      RETURN
      END

C *********************************************************************
C 
      SUBROUTINE BETWEEN(SECT,NERR)
      IMPLICIT NONE
      REAL*8 DD,DT1,DT2,DT3,DIS(2),SECT(3),TLN(2,3)
      INTEGER I,NERR,h
      COMMON /TLN/TLN
      NERR = 1
      DT1=TLN(1,1)-TLN(2,1)
      DT2=TLN(1,2)-TLN(2,2)
      DT3=TLN(1,3)-TLN(2,3)
      DD=DSQRT(DT1**2+DT2**2+DT3**2)
      DO I=1,2
         DT1=TLN(I,1)-SECT(1)
         DT2=TLN(I,2)-SECT(2)
         DT3=TLN(I,3)-SECT(3)
d 	 write(6,*) dt1,dt2,dt3
	 if (abs(dt1).gt.(100000000)) then
		h = h + 1
C		write(6,*) h
		dt1=dt1/(10000000000000000.0d0)
		dt2=dt2/(10000000000000000.0d0)
		dt3=dt3/(10000000000000000.0d0)
C		write(6,*) dt1,dt2,dt3
	 end if
         DIS(I)=DSQRT(DT1**2+DT2**2+DT3**2)
      END DO
      IF (DIS(1).GT.DD.OR.DIS(2).GT.DD) THEN
         NERR=0
      END IF

C
      RETURN
      END 
C *********************************************************************
C
      SUBROUTINE HITEDGE(SECT,NERR)
C
C
      IMPLICIT NONE
      REAL*8 VERT(10,3),SECT(3),TOL
      REAL*8 DS,DS1,DS2,DIF
      PARAMETER (TOL=5.D-05)
      INTEGER I,K,NERR,NVERT
      COMMON /VERT/VERT,NVERT
C

      I=1
      NERR = 0
      DO WHILE(I.LE.NVERT.AND.NERR.EQ.0) 
         DS=0.D0
         DS1=0.D0
         DS2=0.D0
         IF (I.NE.NVERT) THEN
            DO K=1,3
               DS=DS+(VERT(I,K)-VERT(I+1,K))**2
               DS1=DS1+(VERT(I,K)-SECT(K))**2
               DS2=DS2+(VERT(I+1,K)-SECT(K))**2
            END DO
         ELSE
            DO K=1,3
               DS=DS+(VERT(I,K)-VERT(1,K))**2
               DS1=DS1+(VERT(I,K)-SECT(K))**2
               DS2=DS2+(VERT(1,K)-SECT(K))**2
            END DO
         END IF
         DS=DSQRT(DS)
         DS1=DSQRT(DS1)
         DS2=DSQRT(DS2)
         DIF=DABS(DS-DS1-DS2)
C
         IF (DIF.LE.TOL) THEN
            NERR=1
         END IF
         I=I+1
      END DO       
      RETURN
      END
C ************************************************************************
C
      SUBROUTINE INOUT(SECT,FLAG)
C
      IMPLICIT NONE
      REAL*8   VECT1(3),VECT2(3),CRSPRCT(3)
      REAL*8   VERT(10,3),SECT(3)
      INTEGER I,J,SINAL(3),NSINAL(3),FLAG,NVERT
      COMMON /VERT/VERT,NVERT
C
C
C     COMPUTE THE CROSS PRODUCT OF THE FIRST TWO VECTORS FROM THE FIRST TWO
C     VERTICES AND INTERSECTION POINT. THE SINAL OF THIS CURL VECTOR WILL 
C     BE USED AS A REFERENCE TO TEST THE REST OF THE CROSS PRODUCT VECTORS.
C
      FLAG=1
      DO J=1,3
         VECT1(J)=VERT(1,J)-SECT(J)
         VECT2(J)=VERT(2,J)-SECT(J)
      END DO
C 
C     COMPUTE THE COMPONENTS OF THE COMPONENTS OF THE FIRST CROSS PRODUCT
C     
      CRSPRCT(1)=VECT1(2)*VECT2(3)-VECT1(3)*VECT2(2)
      CRSPRCT(2)=VECT1(3)*VECT2(1)-VECT1(1)*VECT2(3)
      CRSPRCT(3)=VECT1(1)*VECT2(2)-VECT1(2)*VECT2(1)
      CALL CRSPRCTSIGN(CRSPRCT,SINAL)
C
C     TEST IF THE OTHER CROSS PRODUCT VECTORS HAVE THE SAME DIRECTIONS 
C     AS THE FIRST ONE. IF YES ,THEN "IN". IF NO,THEN "OUT". 
C
      I=2
      DO WHILE (FLAG.EQ.1.AND.I.LE.NVERT)
         DO J=1,3
            VECT1(J)=VERT(I,J)-SECT(J)
            IF (I.NE.NVERT) THEN
               VECT2(J)=VERT(I+1,J)-SECT(J)
            ELSE
               VECT2(J)=VERT(1,J)-SECT(J)     
            END IF
         END DO
C
C        COMPUTE THE COMPONENTS OF THE CROSS PRODUCT VECTORS OF THE REST
C        OF THE VERTICES.
C
         CRSPRCT(1)=VECT1(2)*VECT2(3)-VECT1(3)*VECT2(2)
         CRSPRCT(2)=VECT1(3)*VECT2(1)-VECT1(1)*VECT2(3)
         CRSPRCT(3)=VECT1(1)*VECT2(2)-VECT1(2)*VECT2(1)
C
         CALL CRSPRCTSIGN(CRSPRCT,NSINAL)
         IF (NSINAL(1).EQ.SINAL(1).AND.NSINAL(2).EQ.SINAL(2).
     &                           AND.NSINAL(3).EQ.SINAL(3)) THEN
            FLAG=1
         ELSE
            FLAG=0
         END IF
         I=I+1
      END DO
      RETURN
      END
C
C ********************************************************************
C
      SUBROUTINE CRSPRCTSIGN(CRSPRCT,SINAL)

      IMPLICIT NONE
      REAL*8 TOL
      PARAMETER (TOL=5.D-05)
      REAL*8 CRSPRCT(3)
      INTEGER I,SINAL(3)
C
      DO I=1,3
         IF (DABS(CRSPRCT(I)).LT.TOL) THEN
            SINAL(I)=0
         ELSE
            IF (CRSPRCT(I).LT.0.0D0) THEN
               SINAL(I)=-1
            ELSE
               SINAL(I)=1
            END IF
         END IF
      END DO
      RETURN
      END         
*************************************************************************

 

Return to Appendix 10.11 List of Programs