| |
Fundamental TechnologiesVoyager LECP Pages |
by Sheela Shodhan
************************************************************************* * 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)
* PURPOSE : TO SCAN THROUGH ALL THE POLYGONS, TO DETERMINE THE FATE OF * * PARTICLE * *************************************************************************
IMPLICIT NONE
INTEGER I,J,HIT,MAXSURF,NSURF,NTY,NVERT
PARAMETER (MAXSURF=50)
REAL*8 A,B,C,COEFF(MAXSURF,4),D,VERT(10,3)
COMMON /COEFF/COEFF,/VERT/VERT,NVERT,/NTY/NTY
D WRITE(6,*) ' ---------- '
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)
CALL SCANPOLY(A,B,C,D,HIT,NSURF)
NSURF=NSURF+1
END DO
D WRITE(6,*) 'HIT: ',HIT,'SURFACE NO.: ',NSURF-1
RETURN
END
C---------------------------------------------------------------------------
C---------------------------------------------------------------------------
SUBROUTINE FINDVERT(VERT,NVERT,NSURF)
* PURPOSE : INITIALIZES THE ARRAY VERT WITH THE VERTICES FOR A GIVEN * * SURFACE NSURF * ****************************************************************************
INTEGER I1,I2,J,MAXSURF,MAXCOO,NVERT
PARAMETER (MAXSURF=50,MAXCOO=31)
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----------------------------------------------------------------------------
C----------------------------------------------------------------------------
SUBROUTINE SCANPOLY(A,B,C,D,HIT,I)
* PURPOSE : TO DETERMINE WHETHER THE PARTICLE HAS HIT THE PLANE SPECIFIED * * BY THE COEFFICIENTS A,B,C & D AND LOST OR NOT. * * VARIABLES: * * HIT : WHOSE VALUE INDICATES WHETHER TRAJECTORY CALCULATION IS TO BE * * CONTINUED OR NOT. * * = 0 PARTICLE IS NOT LOST-SO CONTINUE WITH THE TRAJECTORY CALCULATION.* * = 1 PARTICLE IS LOST-DO NOT CONTINUE WITH THE TRAJECTORY CALCULATION.* * = 2 PARTICLE ESCAPES FROM THE SENSOR ASSEMBLY WITHOUT HITTING * * ANY OF THE OTHER SURFACES. NO NEED TO CONTINUE WITH THE * * TRAJECTORY CALCULATION - SO STOP CALCULATION. * * NGAMAP,NLOWAP,NOPAP : APERTURES OF THE HOUSINGS OF THE DETECTORS,DEFLECTIO* * N SYSTEM AND SENSOR SUBSYSTEM RESPECTIVELY * * SECT : COORDINATES OF THE INTERSECTION POINT * *****************************************************************************
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 (NGAMAP=14,NLOWAP=19,NOPAP=26)
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,*) '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 .EQ. NOPAP) 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----------------------------------------------------------------------------
C----------------------------------------------------------------------------
SUBROUTINE INTERSECT(SECT,A,B,C,D,NERR)
*
* PURPOSE : TO FIND THE INTERSECTION POINT OF THE PLANE AND THE LINE FORMED *
* BY THE LINE-SEGMENT *
*****************************************************************************
IMPLICIT NONE
REAL*8 A,B,C,D,DX,DY,DZ,DET,RATIO
REAL*8 TLN(2,3),SECT(3),NUM
INTEGER NERR
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)
DET=A*DX+B*DY+C*DZ
NUM=-(A*TLN(1,1)+B*TLN(1,2)+C*TLN(1,3)+D)
IF (DET.EQ. 0.0) THEN
NERR=0
RETURN
END IF
RATIO=NUM/DET
C
SECT(1)=DX*RATIO+TLN(1,1)
SECT(2)=DY*RATIO+TLN(1,2)
SECT(3)=DZ*RATIO+TLN(1,3)
C
RETURN
END
C----------------------------------------------------------------------------
C----------------------------------------------------------------------------
SUBROUTINE BETWEEN(SECT,NERR)
*
* PURPOSE : DETERMINE WHETHER THE INTERSECTION POINT LIES BETWEEN THE TWO *
* END POINTS OF THE SEGMENT OR NOT *
*****************************************************************************
IMPLICIT NONE
REAL*8 DD,DT1,DT2,DT3,DIS(2),SECT(3),TLN(2,3)
INTEGER I,NERR
COMMON /TLN/TLN
C
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)
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)
*
* PURPOSE : TO DETERMINE WHETHER THE INTERSECTION POINT LIES ON THE EDGE OF *
* THE POLYGON OR NOT *
*****************************************************************************
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)
* PURPOSE : TO DETERMINE WHETHER THE INTERSECTION POINT LIES INSIDE THE * * BOUNDED POLYGON OR NOT * ******************************************************************************
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 CORSS 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----------------------------------------------------------------------------
SUBROUTINE CRSPRCTSIGN(CRSPRCT,SINAL)
* PURPOSE : IT DETERMINES THE SIGN OF THE CROSS-PRODUCT. * *****************************************************************************
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---------------------------------------------------------------------------- C----------------------------------------------------------------------------
Return to thesis table of contents.
Return to Voyager LECP Data Analysis Handbook Table of
Contents.
Return to Fundamental Technologies Home Page.
Last modified 12/9/02, Tizby Hunt-Ward
tizby@ftecs.com