| |
Fundamental TechnologiesCassini MIMI Pages |
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
****************************************************************************** * * * PROGRAM I.3 * * * ****************************************************************************** * * * TRAJRT4BDET * * THIS SUBROUTINE SHOULD BE LINKED WITH TRAJ4PT. IT NOT ONLY INCLUDES * * OUTP, BUT ALSO THE DIFFERENTIAL EQUATION SOLVER TDHPCG. * * * ******************************************************************************
SUBROUTINE DHPCG(PRMT,PV,DERY,NDIM,IHLF,FCT,OUTP,AUX)
implicit none
REAL*8 PRMT(5),X,H,Z,DELT,PV(6),DERY(6),AUX(16,6)
INTEGER NHIT,wscatter,scatter,i,n,istep,imod,ihlf,ihlf1,isw,ndim
common /nhit/nhit,/wscatter/wscatter,/scatter/scatter
integer startnsurf,ape
logical firstcheck,secondcheck
common /startnsurf/startnsurf,/ape/ape
common /firstcheck/firstcheck,/secondcheck/secondcheck
scatter=0
NHIT=0
N=1
IHLF1=0
IHLF=0
X=PRMT(1)
H=PRMT(3)
PRMT(5)=0.D0
DO I=1,NDIM
AUX(16,I)=0.D0
AUX(15,I)=DERY(I)
AUX(1,I)=PV(I)
ENDDO
C ERROR RETURNS
IF((H*PRMT(2)).LT.X*H) THEN
IHLF=13
ELSE IF((H*PRMT(2)).EQ.X*H) THEN
IHLF=12
ENDIF
ISW=0
C ****************************************************************
C
C COMPUTATION OF DERY FOR STARTING VALUES
C BLOCK 1 VS OLD LABELS 4_20
C
C ***************************************************************
DO WHILE (ISW.NE.4)
IF(ISW.EQ.0) THEN
IHLF1=0
CALL FCT(X,PV,DERY)
C RECORDING OF STARTING VALUES
CALL OUTP(X,PV,DERY,IHLF,NDIM,PRMT,NHIT)
IF (NHIT.NE.0) RETURN
IF(PRMT(5).NE.0.) THEN
RETURN
ENDIF
DO I=1,NDIM
AUX(8,I)=DERY(I)
ENDDO
C COMPUTATION OF AUX(2,I)
ISW=1
ELSE IF(ISW.EQ.1) THEN
IF(IHLF1.NE.1) THEN
X=X+H
DO I=1,NDIM
AUX(2,I)=PV(I)
ENDDO
ENDIF
IHLF1=0
C INCREMENT IS TESTED BY MEANS OF BISECTION
IHLF=IHLF+1
X=X-H
DO I=1,NDIM
AUX(4,I)=AUX(2,I)
ENDDO
H=.5D0*H
N=1
ISW=2
ELSE IF(ISW.EQ.2) THEN
X=X+H
CALL FCT(X,PV,DERY)
N=2
DO I=1,NDIM
AUX(2,I)=PV(I)
AUX(9,I)=DERY(I)
ENDDO
ISW=3
ELSE IF(ISW.EQ.3) THEN
C COMPUTATION OF TEST VALUE DELT
DELT=0.D0
DO I=1,NDIM
DELT=DELT+AUX(15,I)*DABS(PV(I)-AUX(4,I))
ENDDO
DELT=.66666666667D-1*DELT
if (delt.le.prmt(4)) then
c IF(DELE.LE.PRMT(4)) THEN
ISW=5
ELSE
IF(IHLF.GE.10) THEN
IHLF=11
X=X+H
ISW=0
IHLF1=2
ELSE
C NO SATISFACTORY ACCURACY AFTER 10 BISECTIONS. ERROR MESSAGE.
IHLF1=1
ISW=1
ENDIF
ENDIF
ELSE IF(ISW.EQ.5) THEN
C THERE IS SATISFACTORY ACCURACY AFTER LESS THAN 11 BISECTIONS.
X=X+H
CALL FCT(X,PV,DERY)
DO I=1,NDIM
AUX(3,I)=PV(I)
AUX(10,I)=DERY(I)
ENDDO
N=3
ISW=4
ENDIF
C *******************************************************************
C
C THE FOLLOWING PART OF SUBROUTINE HPCG COMPUTES BY MEANS OF
C RUNGE-KUTTA METHOD STARTING VALUES FOR THE NOT SELF-STARTING
C PREDICTOR-CORRECTOR METHOD.
C BLOCK 2 VS OLD LABELS 100_104
C
C *******************************************************************
IF(ISW.NE.5.AND.IHLF1.EQ.0) THEN
DO I=1,NDIM
Z=H*AUX(N+7,I)
AUX(5,I)=Z
PV(I)=AUX(N,I)+.4D0*Z
ENDDO
C Z IS AN AUXILIARY STORAGE LOCATION
Z=X+.4D0*H
CALL FCT(X,PV,DERY)
DO I=1,NDIM
Z=H*DERY(I)
AUX(6,I)=Z
PV(I)=AUX(N,I)+.29697760925D0*AUX(5,I)+.15875964497D0*Z
ENDDO
Z=X+.45573725422D0*H
CALL FCT(X,PV,DERY)
DO I=1,NDIM
Z=H*DERY(I)
AUX(7,I)=Z
PV(I)=AUX(N,I)+0.21810038823D0*AUX(5,I)-3.0509651487D0*AUX(6,I)
* + 3.8328647605D0*Z
ENDDO
Z=X+H
CALL FCT(X,PV,DERY)
DO I=1,NDIM
PV(I)=AUX(N,I)+.17476028223D0*AUX(5,I)-.55148066288D0*AUX(6,I)
* +1.2055355994D0*AUX(7,I)+.17118478122D0*H*DERY(I)
ENDDO
ENDIF
ENDDO
C *******************************************************************
C
C BLOCK 3 VS OLD LABELS 21_22
C
C *******************************************************************
N=1
X=X+H
CALL FCT(X,PV,DERY)
X=PRMT(1)
DO I=1,NDIM
AUX(11,I)=DERY(I)
PV(I)=AUX(1,I)+H*(.375D0*AUX(8,I)+.79166666667D0*AUX(9,I)
* -.20833333333D0*AUX(10,I)+.41666666667D-1*DERY(I))
ENDDO
C *********************************************************************
C
C BLOCK 4 VS OLD LABELS 23_30
C
C *********************************************************************
DO WHILE (N.LT.4)
X=X+H
N=N+1
CALL FCT(X,PV,DERY)
CALL OUTP(X,PV,DERY,IHLF,NDIM,PRMT,NHIT)
IF (NHIT.NE.0) RETURN
IF(PRMT(5).NE.0) THEN
RETURN
ELSE IF((N-4).LT.0) THEN
DO I=1,NDIM
AUX(N,I)=PV(I)
AUX(N+7,I)=DERY(I)
ENDDO
IF(N.LT.3) THEN
DO I=1,NDIM
DELT=AUX(9,I)+AUX(9,I)
DELT=DELT+DELT
PV(I)=AUX(1,I)+.33333333333D0*H*(AUX(8,I)+DELT+AUX(10,I))
ENDDO
ELSE
DO I=1,NDIM
DELT=AUX(9,I)+AUX(10,I)
DELT=DELT+DELT+DELT
PV(I)=AUX(1,I)+.375D0*H*(AUX(8,I)+DELT+AUX(11,I))
ENDDO
ENDIF
ENDIF
ENDDO
ISW=6
C ******************************************************************
C
C POSSIBLE BREAK POINT FOR LINKAGE
C STARTING VALUES ARE COMPUTED
C NOW START HAMMINGS MODIFIED PREDICTOR-CORRECTOR METHOD.
C BLOCK 5 VS OLD LABELS 200_226
C
C ******************************************************************
ISTEP=3
DO WHILE (ISW.GE.6)
IF(ISW.EQ.6) THEN
IF(N.EQ.8) THEN
C N=8 CAUSES THE ROWS OF AUX TO CHANGE THEIR STORAGE LOCATIONS
DO N=2,7
DO I=1,NDIM
AUX(N-1,I)=AUX(N,I)
AUX(N+6,I)=AUX(N+7,I)
ENDDO
ENDDO
N=7
ENDIF
C N LESS THAN 8 CAUSES N+1 TO GET N
N=N+1
C COMPUTATION OF NEXT VECTOR PV
DO I=1,NDIM
AUX(N-1,I)=PV(I)
AUX(N+6,I)=DERY(I)
ENDDO
X=X+H
ISW=7
ELSE IF(ISW.EQ.7) THEN
ISTEP=ISTEP+1
DO I=1,NDIM
DELT=AUX(N-4,I)+1.3333333333D0*H*(AUX(N+6,I)+AUX(N+6,I)
* -AUX(N+5,I)+AUX(N+4,I)+AUX(N+4,I))
PV(I)=DELT-.92561983471D0*AUX(16,I)
AUX(16,I)=DELT
ENDDO
C PREDICTOR IS NOW GENERATED IN ROW 16 OF AUX, MODIFIED PREDICTOR
C IS GENERATED I PV. DELT MEANS AN AUXILIARY STORAGE.
CALL FCT(X,PV,DERY)
C DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY
DO I=1,NDIM
DELT=.125D0*(9.D0*AUX(N-1,I)-AUX(N-3,I)+3.D0*H*(DERY(I)
* +AUX(N+6,I)+AUX(N+6,I)-AUX(N+5,I)))
AUX(16,I)=AUX(16,I)-DELT
PV(I)=DELT+.7438016529D-1*AUX(16,I)
ENDDO
C TEST WHETHER H MUST BE HALVED OR DOUBLED
DELT=0.D0
DO I=1,NDIM
DELT=DELT+AUX(15,I)*DABS(AUX(16,I))
ENDDO
ISW=8
ELSE IF(ISW.EQ.8) THEN
IF(DELT.LT.PRMT(4).OR.IHLF1.EQ.5) THEN
IHLF1=0
C H MUST NOT BE HALVED. THAT MEANS PV(I) ARE GOOD.
CALL FCT(X,PV,DERY)
CALL OUTP(X,PV,DERY,IHLF,NDIM,PRMT,NHIT)
IF (NHIT.NE.0) RETURN
IF(PRMT(5).NE.0.OR.IHLF.GE.11.OR.(H*(X-PRMT(2))).GE.0
* .OR.DABS(X-PRMT(2)).LT.0.1D0*DABS(H)) THEN
1150 FORMAT(I5)
RETURN
ELSE
IF(DELT.LE.(.2D-1*PRMT(4)).AND.IHLF.GT.0.AND.N.GE.7
* .AND.ISTEP.GE.4) THEN
C H COULD BE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE AVAILABLE
IMOD=ISTEP/2
IF((ISTEP-IMOD-IMOD).EQ.0) THEN
H=H+H
IHLF=IHLF-1
ISTEP=0
DO I=1,NDIM
AUX(N-1,I)=AUX(N-2,I)
AUX(N-2,I)=AUX(N-4,I)
AUX(N-3,I)=AUX(N-6,I)
AUX(N+6,I)=AUX(N+5,I)
AUX(N+5,I)=AUX(N+3,I)
AUX(N+4,I)=AUX(N+1,I)
DELT=AUX(N+6,I)+AUX(N+5,I)
DELT=DELT+DELT+DELT
AUX(16,I)=8.962962963D0*(PV(I)-AUX(N-3,I))-3.361111111D0*H*
* (DERY(I)+DELT+AUX(N+4,I))
ENDDO
ISW=6
ELSE
ISW=6
ENDIF
ELSE
ISW=6
ENDIF
ENDIF
ISW=6
ELSE
ISW=9
ENDIF
ELSE IF(ISW.EQ.9) THEN
C H MUST BE HALVED
IHLF=IHLF+1
IF(IHLF.LE.10) THEN
ISW=10
ELSE
ISW=8
IHLF1=5
ENDIF
ELSE IF(ISW.EQ.10) THEN
H=.5D0*H
ISTEP=0
DO I=1,NDIM
PV(I)=.390625D-2*(8.D1*AUX(N-1,I)+135.D0*AUX(N-2,I)
* + 4.D1*AUX(N-1,I)+AUX(N-4,I))-.1171875D0*(AUX(N+6,I)
* - 6.D0*AUX(N+5,I)-AUX(N+4,I))*H
AUX(N-4,I)=.390625D-2*(12.D0*AUX(N-1,I)+135.D0*AUX(N-2,I)+
* 08.D0*AUX(N-3,I)+AUX(N-4,I))-.234375D-1*(AUX(N+6,I)+
* 8.D0*AUX(N+5,I)-9.D0*AUX(N+4,I))*H
AUX(N-3,I)=AUX(N-2,I)
AUX(N+4,I)=AUX(N+5,I)
ENDDO
X=X-H
DELT=X-(H+H)
CALL FCT(X,PV,DERY)
DO I=1,NDIM
AUX(N-2,I)=PV(I)
AUX(N+5,I)=DERY(I)
PV(I)=AUX(N-4,I)
ENDDO
DELT=DELT-(H+H)
CALL FCT(X,PV,DERY)
DO I=1,NDIM
DELT=AUX(N+5,I)+AUX(N+4,I)
DELT=DELT+DELT+DELT
AUX(16,I)=8.9629696296D0*(AUX(N-1,I)-PV(I))
* - 3.3611111111D0*H*(AUX(N+6,I)+DELT+DERY(I))
AUX(N+3,I)=DERY(I)
ENDDO
ISW=7
ENDIF
ENDDO
c print*, nhit
c pause
return
END C----------------------------------------------------------------------- C----------------------------------------------------------------------- C THIS IS A FILE OF ROUTINES THAT HAVE BEEN DECLARED EXTERNAL IN THE C MAIN PROGRAM & WHICH ARE PASSED INTO THE ROUTINE DHPCG THAT SOLVES THE C SYSTEM OF DIFFERENTIAL EQUATIONS.
SUBROUTINE FCT(X,Y,DERY)
IMPLICIT NONE
INTEGER NDIM
PARAMETER (NDIM=6)
REAL*8 BX,BY,BZ,DERY(NDIM),FAC,QMC,X,X00,Y(NDIM),Y00,Z00
COMMON /QMC/QMC
C COMMON /BX/BX,BY,BZ,/X00/X00,Y00,Z00
PARAMETER (FAC=1.0D-2)
C SUBROUTINE TO COMPUTE THE RIGHT HAND SIDE DERY OF THE SYSTEM AT THE
C GIVEN VALUES OF X & Y.
X00 = Y(1)/FAC
Y00 = Y(2)/FAC
Z00 = Y(3)/FAC
CALL FDMOD(X00,Y00,Z00,BX,BY,BZ)
D PRINT *,' QMC ',QMC
D PRINT *,' X00 ',X00,' Y00 ',Y00,' Z00 ',Z00
D PRINT *,' BX ',BX,' BY ',BY,' BZ ',BZ
DERY(1) = Y(4)
DERY(2) = Y(5)
DERY(3) = Y(6)
DERY(4) = QMC * (Y(5)*BZ - Y(6)*BY)
DERY(5) = QMC * (Y(6)*BX - Y(4)*BZ)
DERY(6) = QMC * (Y(4)*BY - Y(5)*BX)
RETURN
END
C--------------------------------------------------------------------
C---------------------------------------------------------------------
SUBROUTINE OUTP(X,Y,DERY,IHLF,NDIM,PRMT,HIT)
C PURPOSE : THIS IS A ROUTINE THAT IS DECLARED EXTERNAL IN THE MAIN C ROUTINE. IT PRINTS THE OUTPUT VALUES OBTAINED FROM THE ROUTINE C DHPCG - WHICH SOLVES THE SYSTEM OF DIFF. EQUATIONS.
IMPLICIT NONE
CHARACTER*72 FNAME
INTEGER HIT,I,IHLF,MULC,MULC1,NCOUNT,NDIM,NSURF,NU,ek
integer scatter,wscatter,change,Nscatter,nsurf1
REAL*8 DERY(NDIM),PRMT(5),X,Y(NDIM),TLN(2,3),TOTV,BX,BY,BZ,B,
& X00,Y00,Z00,FAC,PASVEL(10000,2)
real*8 x01(3),x02(3),aux(16,6),sect(3),c1(4),theta,phi
REAL*8 FPHI,DIST,X1,Y1,Z1
integer startnsurf,ape
logical firstcheck,secondcheck
common /startnsurf/startnsurf,/ape/ape
common /firstcheck/firstcheck,/secondcheck/secondcheck
external fct
PARAMETER (FAC=1.0D-2)
COMMON /NCOUNT/NCOUNT,/ek/ek
COMMON /TLN/TLN,/NU/NU
common /wscatter/wscatter,/scatter/scatter
common /sect/sect,/nsurf1/nsurf1,/c1/c1,/Nscatter/Nscatter
common /theta/theta,/phi/phi
C COMMON /PASVEL/PASVEL
INCLUDE 'PASS5.CMN'
DATA MULC/40/,MULC1/10/
TOTV=DSQRT(Y(4)*Y(4) + Y(5)*Y(5) + Y(6)*Y(6))
NCOUNT = NCOUNT + 1
D B = DSQRT(BX*BX + BY*BY + BZ*BZ)
IF (NU .EQ. 1) THEN
DO I=1,3
TLN(2,I)=Y(I)/FAC
ENDDO
ELSE
DO I=1,3
TLN(1,I)=TLN(2,I)
TLN(2,I)=Y(I)/FAC
ENDDO
do 31 i=1,3
x01(i)=tln(1,i)
x02(i)=tln(2,i)
31 enddo
c if (secondcheck) then
c write(6,*)'-----',startnsurf,ape
c endif
CALL CHECKHIT(HIT,NSURF)
!HIT=0 :CONTINUE WITH TRAJECTORY CALCULATION
! =1 :HIT THE NSURF(NTH SURFACE), LOST
! =2 :SUCCESSFULLY ESCAPED FROM THE SENSOR
if ( ((x01(1).gt.7.0).or.(x01(1).lt.0.0).or.(x01(2).gt.4)
& .or.(x01(2).lt.1.0).or.(x01(3).gt.1.0).or.(x01(3).lt.-1.0))
& .and. (hit.eq.0) ) then
c write(6,*)'Special case, trajectory has run out of assembly'
c & ,x02(1),x02(2),x02(3)
hit = 1
endif
IF (HIT .EQ. 2) THEN
NPAS = NPAS + 1
X1=TLN(2,1)-TLN(1,1)
Y1=TLN(2,2)-TLN(1,2)
Z1=TLN(2,3)-TLN(1,3)
DIST=DSQRT(X1*X1 + Y1*Y1 + Z1*Z1)
pas(npas,1)=theta
pas(npas,2)=phi
PAS(NPAS,3)=DACOS(Z1/DIST)
PAS(NPAS,4)=FPHI(X1,Y1)
pas(npas,11)=1
C PASVEL(NPAS,1)=DACOS(Y(6)/TOTV)
C PASVEL(NPAS,2)=FPHI(Y(4),Y(5))
DO I=1,6
PAS(NPAS,I+4)=Y(I)
ENDDO
ENDIF
C if ((hit.eq.0).and.(scatter.eq.1)) then
C if (wscatter.eq.1) then
C call SCATTERSUB(HIT,sect,X01,X02,C1,y)
C nsurf1 = nsurf - 1
C if (hit.eq.0) then
C call fct(x,y,dery)
C call DHPCG1(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP1,AUX)
C endif
C endif
C if (wscatter.eq.2) then
C do i=1,nscatter
C call SCATTERSUB(HIT,sect,X01,X02,C1,y)
C nsurf1 = nsurf - 1
C if (hit.eq.0) then
C call fct(x,y,dery)
C call DHPCG1(PRMT,Y,DERY,NDIM,IHLF,
C & FCT,OUTP1,AUX)
C------------------------------------------------------------------------
C Here is a modification to restrict only one diffuse scattering
C occurs during the process, Be careful with this modification
C need more testing for this.
C
C if (hit.eq.2) then
C nu=nu+1
C return
C endif
C------------------------------------------------------------------------
C endif
C enddo
C endif
C hit=1
C endif
ENDIF
NU=NU+1
D WRITE(3,15) X,(Y(I),I=1,NDIM),TOTV C15 FORMAT (1X,D14.7,6(1X,D14.7),1X,D14.7,1X,F10.3) 15 FORMAT (1X,D14.7,6(1X,D15.8),1X,D21.14) 20 FORMAT(3(X,F13.6))
RETURN
END **************************************************************************
Return to Technical Report
table of contents.
Return to Cassini MIMI table of contents
page.
Return to Fundamental Technologies Home Page.
Updated 6/6/02, T. Hunt-Ward
tizby@ftecs.com