Cassini MIMI Investigation at Fundamental Technologies

Numerical Computation of Energy-Dependent Geometric Factors of E and F Electron Detectors of CASSINI/MIMI/LEMMS

Technical Report by Xiaodong Hong and Thomas P. Armstrong, May 10, 1997

*
*                                                                            *
******************************************************************************
*                                                                            *
*    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,scatter,ape,startnsurf
      PARAMETER (MAXSURF=62)
      REAL*8 A,B,C,COEFF(MAXSURF,4),D,VERT(10,3)
      logical done
      COMMON /COEFF/COEFF,/VERT/VERT,NVERT,/NTY/NTY,/scatter/scatter
      common /startnsurf/startnsurf, /ape/ape

      done=.false.
C      edgehit = .false.	
      nsurf = startnsurf
      if (hit.eq.1) then
	done = .true.
      end if	
      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------------------------------------------------------------------
C------------------------------------------------------------------

      SUBROUTINE FINDVERT(VERT,NVERT,NSURF)

      INTEGER I1,I2,J,MAXSURF,MAXCOO,NVERT
      PARAMETER (MAXSURF=62,MAXCOO=62)
      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,nsurf)

 
      IMPLICIT NONE
      REAL*8 A,B,C,D,SECT(3),TLN(2,3),VERT(10,3),c1(4)
      INTEGER FLAG,HIT,NCOUNT,NGAMAP,NLOWAP,NERR,NVERT,NOPAP,
     & scatter,wscatter,ape, startnsurf, nty, nsurf 
      logical done
	common /scatter/scatter,/wscatter/wscatter
	common /sect/sect,/c1/c1,/ape/ape, /nty/nty, /startnsurf/startnsurf
      PARAMETER (NOPAP=63,NGAMAP=63,NLOWAP=63)
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,*) nsurf
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.(nsurf.eq.nopap)) then
		hit = 1
		scatter = 1
c		edgehit = .true.
	    else           
	    	scatter=1
	    	hit=0
	    end if

            done=.true.

D	write(6,*) wscatter,scatter,sect(1),sect(2),sect(3)
D       WRITE(6,*) 'IT LIES ON THE EDGES OF THE POLYGON',nsurf
	  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 ( nsurf .eq.ape ) THEN     ! it escapes
                HIT=2
		done = .true.
D	write(6,*) wscatter,scatter,sect(1),sect(2),sect(3)
D        WRITE(6,*) 'IT LIES INSIDE THE OUTER OPENING APERUTRE',nsurf
D	Write(6,*) 'It ESCAPES!!!!!!!!!!!'

                RETURN
	    END IF

	    if (nsurf.le.5) then
		hit = 1
		done = .true.      !hit walls other than apertures
		

		return	
            end if

	    if (nsurf .eq. 6) then  !potential to escape the apertures
		hit = 0             !need furtur check
D		WRITE(6,*) 'It hit polygon #6'
		return
	    end if
 	

	    if ((nsurf.ge.7) .or. (nsurf.le.13)) then
		hit = 0
D		WRITE(6,*) 'It is on the way from ', nsurf
		nsurf = (nsurf-7)*7 + 13      ! it would add 1 later  
	        startnsurf = nsurf + 1        ! set the new start polygon number
	 	ape = startnsurf + 6          ! set the new escape polygon number
		nty = ape
		return
	    end if	

D            WRITE(6,*) 'It hits one wall and lost'
       	    done = .true.                     ! trajectory hit one of walls in a channel

	    return
	else	

	    RETURN                            ! trajectory has not lost yet

      	END IF
   

      	IF (FLAG .NE. 0) THEN

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

            
	HIT=0                          !OUTSIDE THE OTHER PLANE     

c       WRITE(6,*) 'LIES OUTSIDE OTHER PLANES-SO NOT LOST YET'
                          !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 D 

Return to Technical Report Table of Contents

Return to Cassini MIMI table of contents page.
Return to Fundamental Technologies Home Page.

QUICK FACTS

RELATED LINKS