*		        PROGRAM I.2				   *

*******************************************************************************
*                            CHECHITSEN1                                      *
* THIS FILE CONTAINS ALL THE SUBROUTINES THAT ARE USED TO FIND WHETHER        *
* THE LINE SEGMENT HAS HIT ANY OF THE SURFACES AND LOST OR NOT.               *
*******************************************************************************

      SUBROUTINE CHECKHIT(HIT,NSURF)

      IMPLICIT NONE
      INTEGER I,J,HIT,MAXSURF,NSURF,NTY,NVERT
      PARAMETER (MAXSURF=55)
      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 SCANPOLY(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=55,MAXCOO=55)
      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(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 
      PARAMETER (NOPAP=40,NGAMAP=56,NLOWAP=57)
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
        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
         WRITE(6,*) 'IT LIES INSIDE THE LOWER POLYGON'
        RETURN
       ELSE
        IF (I .EQ. NOPAP) THEN      	!INSIDE THE OUTER APERTURE
         HIT=2
          WRITE(6,*) 'IT LIES INSIDE THE OUTER OPENING APERUTRE'
         RETURN
        END IF
       END IF
       HIT=1     			!INSIDE THE OTHER POLYGONS        
        WRITE(6,*) 'IT LIES INSIDE THE OTHER POLYGONS'
       RETURN
      END IF
   

      IF (FLAG .NE. 0) THEN
       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                                     
        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 ***********************************************************************
C
      SUBROUTINE INTERSECT(SECT,A,B,C,D,NERR)
C
      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
		write(6,*) h
		dt1=dt1/(10000000000000000.0d0)
		dt2=dt2/(10000000000000000.0d0)
		dt3=dt3/(10000000000000000.0d0)
		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         
C ************************************************************************

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