Voyager LECP Data Analysis Handbook

Instrument Modeling Reports

An Analysis of the Performance of the Magnetic Deflection System
 in the Voyager Low Energy Charged Particle Experiment

by Sheela Shodhan

E.5 TDHPCG

*******************************************************************
* THIS PROGRAM IS THE DIFFERENTIAL EQUATION SOLVER-DHPCG WHICH    *
* USES THE HAMMING'S MODIFIED FOURTH ORDER PREDICTOR-CORRECTOR    *
* DOUBLE-PRECISION METHOD.                                        *
* THIS IS A MODIFICATION OF SUBROUTINE DHPCG FROM THE             *
* SCIENTIFIC SUBROUTINE PACKAGE (SSP) PUBLISHED BY                *
* IBM. THE ORIGINAL VERSION OF THE ROUTINE WAS WRITTEN IN         *
* FORTRAN-4 WITH NO INDENTATION AS WELL AS POORLY                 *
* STRUCTURED.                                                     *
* THE FOLLOWING VERSION IS VERY MUCH MODERNIZED TO A STRUCTURED   *
* CODE AND INDENTED IN FORTRAN-77. THE USERS ARE REFERED TO THE   *
* MANUAL OF SSP FOR FURTHER IMFORMATION ABOUT THE METHODS AND     *
* DETAILED DESCRIPTION.                                           *
*******************************************************************

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

      REAL*8 PRMT(1),X,H,Z,DELT,PV(1),DERY(1),AUX(16,1)
      INTEGER NHIT
      COMMON /NHIT/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(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-----------------------------------------------------------------------

Return to thesis table of contents. 

Return to Voyager LECP Data Analysis Handbook Table of Contents.
Return to Fundamental Technologies Home Page.

VOYAGER 1 ELAPSED TIME

VOYAGER 2 ELAPSED TIME

QUICK FACTS

RELATED LINKS