Ulysses HISCALE

ULYSSES

Ulysses HISCALE Data Analysis Handbook

Appendix 9 Geometric Factor Study for the Deflected and Unscattered Electrons of HISCALE (Buckley MS Thesis)

*			  PROGRAM I.13			    *

*******************************************************************************
*                                  TRAJRT1BDET                                *
* THIS SUBROUTINE SHOULD BE LINKED WITH TRAJ2PT.  IT NOT ONLY INCLUDES        *
* TRAJRT2DET, BUT ALSO THE DIFFERENTIAL EQUATION SOLVER TDHPCG.               *
*******************************************************************************

      SUBROUTINE DHPCG(PRMT,PV,DERY,NDIM,IHLF,FCT,OUTP,AUX)

      REAL*8 PRMT(5),X,H,Z,DELT,PV(6),DERY(6),AUX(16,6)
      INTEGER NHIT

      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)  
      END DO
C     ERROR RETURNS
      IF((H*PRMT(2)).LT.X*H) THEN
         IHLF=13
      ELSE IF((H*PRMT(2)).EQ.X*H) THEN
         IHLF=12
      END IF
      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
            END IF
            DO I=1,NDIM  
               AUX(8,I)=DERY(I) 
            END DO
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)  
               END DO
            END IF
            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)  
            END DO
            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
      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=0.02540005D0)
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,NPAS,NDIM,NSURF,NU,NTRAJ
      REAL*8 DERY(NDIM),PRMT(5),X,Y(NDIM),TLN(2,3),TOTV,BX,BY,BZ,B,
     &   X00,Y00,Z00,FAC
      PARAMETER (FAC=0.02540005D0)
      COMMON /NCOUNT/NCOUNT,/FNAME/FNAME
C     COMMON /BX/BX,BY,BZ,/X00/X00,Y00,Z00
      COMMON /TLN/TLN,/NU/NU,/NPAS/NPAS
      DATA MULC/40/,MULC1/5/

C     OPEN(UNIT=3,FILE='TRAJ.DAT',ACCESS='SEQUENTIAL',STATUS='OLD')
      OPEN(UNIT=4,FILE=FNAME,ACCESS='SEQUENTIAL',STATUS='OLD')
      OPEN(UNIT=7,FILE='TRAJSH.DAT',ACCESS='SEQUENTIAL',STATUS='OLD')
C     OPEN(UNIT=8,FILE='PASS.DAT',ACCESS='SEQUENTIAL',STATUS='OLD')

      TOTV=DSQRT(Y(4)*Y(4) + Y(5)*Y(5) + Y(6)*Y(6))

      IF (NCOUNT .EQ. 0) THEN
       WRITE(7,15) X,(Y(I),I=1,NDIM),TOTV
       WRITE(4,20) (Y(I)/FAC,I=1,3)
      END IF

      NCOUNT = NCOUNT + 1
      IF (MOD(NCOUNT,MULC) .EQ. 0) THEN
        WRITE(7,15) X,(Y(I),I=1,NDIM),TOTV
      END IF

      IF (MOD(NCOUNT,MULC1) .EQ. 0) THEN
        WRITE(4,20) (Y(I)/FAC,I=1,3)
      END IF

D     B = DSQRT(BX*BX + BY*BY + BZ*BZ)

      IF (NU .EQ. 1) THEN
       DO I=1,3
        TLN(2,I)=Y(I)/FAC
       END DO
      ELSE
       DO I=1,3
        TLN(1,I)=TLN(2,I)
        TLN(2,I)=Y(I)/FAC
       END DO

       CALL CHECKHIT(HIT,NSURF)
                        !HIT=0 :CONTINUE WITH TRAJECTORY CALCULATION
                        !   =1 :HIT THE NSURF(NTH SURFACE), LOST
                        !   =2 :SUCCESSFULLY ESCAPED FROM THE SENSOR 
       IF (HIT .EQ. 1) THEN
         WRITE(4,20) (Y(I)/FAC,I=1,3)
D        WRITE(6,*) 'HIT',HIT,'SURFACE NO.',NSURF-1 
D        WRITE(4,*) 'HIT',HIT,'SURFACE NO.',NSURF-1
       ELSE 
         IF (HIT .EQ. 2) THEN
          WRITE(4,20) (Y(I)/FAC,I=1,3)
D         WRITE(6,*) 'ESCAPED'
D         WRITE(4,*) 'ESCAPED'
          NPAS=NPAS+1
         END IF
       END IF

      END IF 
      NU=NU+1

C     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,F9.6))
      RETURN
      END

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

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

Appendix 9 Geometric Factor Study for the Deflected and Unscattered Electrons of HISCALE (Buckley MS Thesis)

QUICK FACTS

RELATED LINKS