ULYSSES
Ulysses HISCALE Data Analysis Handbook
Appendix 10. Effect of Backscattered Electrons on the Geometric Factors of the LEMS30 Telescope (Hong MS Thesis)
****************************************************************************** * * * 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 Appendix 10.11 List of Programs
Return to the Table of Contents for Hong's MS Thesis
Return to HISCALE List of Appendices
Return to Ulysses HISCALE Data Analysis Handbook Table of Contents
Updated 8/8/19, Cameron Crane
QUICK FACTS
Manufacturer:
ESA provided the Ulysses spacecraft, NASA provided the power
supply, and various others provided its instruments.
Mission End Date: June 30, 2009
Destination: The inner heliosphere of the sun away from the ecliptic plane
Orbit: Elliptical orbit transversing the polar regions of the sun outside of the ecliptic plane
Mission End Date: June 30, 2009
Destination: The inner heliosphere of the sun away from the ecliptic plane
Orbit: Elliptical orbit transversing the polar regions of the sun outside of the ecliptic plane