Menu
(C)JAXA
Akebono satellite
Since 1989
************************************************************************
SUBROUTINE ADAMS
************************************************************************
IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
*** CONSTANT
COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY
*** INPUT
PARAMETER (NSPEC=4)
COMMON /QPARAM/RREF,RBTM,RREFP,
& PF0(NSPEC),GF0(NSPEC),
& XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP,
& RRLI,HWLI,
& XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2,
& RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT,
& NPATH,NLPMAX,INDOUT,INDMSH,IEOF,
& IYEAR,IMNTH,IDATE
COMMON /QINIT /R0,T0,P0,E0,D0,F0
COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY
*** OUTPUT
COMMON /QCOND /L,LL,NN,NLOOP,A
*** Variables below are only used in the ADAMS method.
COMMON /QADAM /Y(8),P(8),F(8),
& H,HOLD,WT(8),PHI(8,14),PSIOLD(12),
& PSI(12),BETA(12),SIGMA(13),G(13),
& ERKM1,ERK,K,KOLD,
& KNEW,KS,IPHASE,IFAIL
REAL*8 SUML(8)
CALL OUT0
*** Initialization of Y
Y(1) = R0
Y(2) = T0
Y(3) = P0
EPS = E0
DEL = D0
FKC = F0
Y(4) = COS(EPS*RAD)*COS(DEL*RAD)
Y(5) = COS(EPS*RAD)*SIN(DEL*RAD)
Y(6) = SIN(EPS*RAD)
Y(7) = 0.0D0
Y(8) = 0.0D0
L = 0
LL = 0
NN = 0
NLOOP = 0
ERR0 = MAX(ERRR,ERRA)
ERRR = ERRR/ERR0
ERRA = ERRA/ERR0
* Derivative for Y
CALL FUNCT(Y,F,1)
yp20 = F(2)
do L=1,8
WT(L)=ERRR*dABS(Y(L)+ERRA)
SUML(L) = (F(L)/WT(L))**2
enddo
SUM = 0.0D0
do L=1,8
SUM = SUM+SUML(L)
enddo
SUM = DSQRT(SUM)
H = 0.00001D0*dSQRT(ERR0/SUM)
HOLD = 0.0D0
K = 1
KOLD = 0
IPHASE = 1
IFAIL = 0
do L=1,8
PHI(L,1) = F(L)
PHI(L,2) = 0.0D0
enddo
CALL OUT1(Y)
DO NLOOP=0,NLPMAX-1
* Adams-Bashforth Predictor
CALL ADAM_B(NLOOP)
* Adams-Moulton Corrector
CALL ADAM_M(ERR0,NLOOP)
*VOCL LOOP,SCALAR
IF(LL.EQ.0) THEN
IF(IFAIL.NE.0) THEN
L=11
ELSE
IF(L.GE.21) THEN
LL=L
ELSE IF(H.LT.HMIN) THEN
L =23
LL=23
ELSE IF(Y(1).LT.RRMIN) THEN
L =24
LL=24
ELSE IF(Y(1).GT.RRMAX) THEN
L =25
LL=25
C ELSE IF((T0-90.0)*F(2).GT.0.0D0) THEN
C L =26
C LL=26
C ELSE IF(yp20*yp(2).lt.0.0D0) THEN
C L =27
C LL=27
ELSE
L =0
NN=NN+1
ENDIF
ENDIF
ENDIF
CALL OUT1(Y)
IF (LL.NE.0) THEN
write(*,*) "Stop Condition: ", LL
GOTO 999
ENDIF
do I=1,7
WT(I)=ERRR*ABS(Y(I)+ERRA)
enddo
ENDDO
999 CONTINUE
write(*,*) 'Number of LOOP=',NLOOP,', Number of calculation=',A
write(*,*) 'Last interval=',H,', IFAIL=',IFAIL
write(21,*)
write(25,*)
RETURN
END
************************************************************************
SUBROUTINE ADAM_B(NLOOP)
* Calculation of Adams-Bashforth Predictor
************************************************************************
IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
COMMON /QADAM /Y(8),P(8),F(8),
& H,HOLD,WT(8),PHI(8,14),PSIOLD(12),
& PSI(12),BETA(12),SIGMA(13),G(13),
& ERKM1,ERK,K,KOLD,
& KNEW,KS,IPHASE,IFAIL
real*8 ALPHA(12),V(12),W(12)
save ALPHA,V,W
G(1) = 1.0D0
G(2) = 0.5D0
SIGMA(1) = 1.0D0
* KS is the maximum order
* In case step size is changed, KS->0.
IF (H .NE.HOLD) KS = 0
IF (KS.LE.KOLD) KS = KS + 1
IF (K.GE.KS) THEN
*** PSI
do I=KS,K-1
PSIOLD(I) = PSI(I)
enddo
PSI(KS) = H*KS
do I=KS+1,K
PSI(I) = PSIOLD(I-1)+H
enddo
*** ALPHA
ALPHA(KS) = 1.0D0/DBLE(KS)
do I=KS+1,K
ALPHA(I) = H/(PSIOLD(I-1)+H)
enddo
*** BETA
BETA(KS) = 1.0D0
do I=KS+1,K
BETA(I) = BETA(I-1)*PSI(I-1)/PSIOLD(I-1)
enddo
*** SIGMA
SIGMA(KS+1) = 1.0D0
do I=KS+1,K
SIGMA(I+1) = I*ALPHA(I)*SIGMA(I)
enddo
*** G (Integral of C)
if (KS.EQ.0) write(0,*) 'Error! KS=0'
if (KS.EQ.1) then
do I=1,K
V(I) = 1.0D0/DBLE(I*(I+1))
W(I) = V(I)
enddo
else
if (K.gt.KOLD) then
V(K) = 1.0D0/DBLE(K*(K+1))
do J=1,KS-2
V(K-J) = V(K-J)-ALPHA(J+1)*V(K-J+1)
enddo
endif
do IQ=1,K+1-KS
V(IQ) = V(IQ)-ALPHA(KS)*V(IQ+1)
W(IQ) = V(IQ)
enddo
G(KS+1) = W(1)
endif
do I=KS+1,K
do IQ=1,K+1-I
W(IQ) = W(IQ)-ALPHA(I)*W(IQ+1)
enddo
G(I+1) = W(1)
enddo
ENDIF
*** PREDICTOR
do L=1,8
do I=KS+1,K
PHI(L,I) = BETA(I)*PHI(L,I)
enddo
PHI(L,K+2) = PHI(L,K+1)
PHI(L,K+1) = 0.0D0
SUM = 0.0D0
do I=K,1, -1
SUM = SUM + G(I)*PHI(L,I)
PHI(L,I) = PHI(L,I) + PHI(L,I+1)
enddo
P(L) = Y(L) + H * SUM
enddo
*** Derivative
CALL FUNCT(P,F,0)
RETURN
END
************************************************************************
SUBROUTINE ADAM_M(ERR0,NLOOP)
* Calculation of Adams-Moulton Corrector
************************************************************************
IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
COMMON /QADAM /Y(8),P(8),F(8),
& H,HOLD,WT(8),PHI(8,14),PSIOLD(12),
& PSI(12),BETA(12),SIGMA(13),G(13),
& ERKM1,ERK,K,KOLD,
& KNEW,KS,IPHASE,IFAIL
DIMENSION GSTR(13)
DATA GSTR/0.500D0 ,0.0833D0 ,0.0417D0 ,0.0264D0 ,0.0188D0 ,
& 0.0143D0 ,0.0114D0 ,0.00936D0,0.00789D0,0.00679D0,
& 0.00592D0,0.00524D0,0.00468D0/
KNEW = K
*** ERK, ERKM1, ERKM2
ERK = 0.0D0
ERKM1 = 0.0D0
ERKM2 = 0.0D0
do L=1,8
TEMP1 = 1.0D0/WT(L)
TEMP2 = F(L)-PHI(L,1)
ERK = ERK+(TEMP2 *TEMP1)**2
IF (K.EQ.2) THEN
ERKM1 = ERKM1+((PHI(L,K )+TEMP2)*TEMP1)**2
ELSE IF (K.GT.2) THEN
ERKM1 = ERKM1+((PHI(L,K )+TEMP2)*TEMP1)**2
ERKM2 = ERKM2+((PHI(L,K-1)+TEMP2)*TEMP1)**2
ENDIF
enddo
TEMP = DABS(H)*DSQRT(ERK)
ERK = TEMP*SIGMA(K+1)*GSTR(K)
IF (k.EQ.2) THEN
ERKM1 = DABS(H)*SIGMA(K )*GSTR(K-1)*DSQRT(ERKM1)
IF(ERKM1.LE.0.5D0*ERK) KNEW = K-1
ELSE IF (K.GT.2) THEN
ERKM1 = DABS(H)*SIGMA(K )*GSTR(K-1)*DSQRT(ERKM1)
ERKM2 = DABS(H)*SIGMA(K-1)*GSTR(K-2)*DSQRT(ERKM2)
IF(DMAX1(ERKM1,ERKM2).LE.ERK) KNEW = K-1
ENDIF
*** Validity check
ERR = TEMP*(G(K)-G(K+1))
IF(ERR.LE.ERR0) THEN
CALL ADAM_UPDATE(ERR0)
ELSE
CALL ADAM_RETRY(ERR0)
ENDIF
RETURN
END
************************************************************************
SUBROUTINE ADAM_UPDATE(ERR0)
************************************************************************
IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
COMMON /QADAM /Y(8),P(8),F(8),
& H,HOLD,WT(8),PHI(8,14),PSIOLD(12),
& PSI(12),BETA(12),SIGMA(13),G(13),
& ERKM1,ERK,K,KOLD,
& KNEW,KS,IPHASE,IFAIL
DIMENSION GSTR(13),TWO(13)
DATA GSTR/0.500D0 ,0.0833D0 ,0.0417D0 ,0.0264D0 ,0.0188D0 ,
& 0.0143D0 ,0.0114D0 ,0.00936D0,0.00789D0,0.00679D0,
& 0.00592D0,0.00524D0,0.00468D0/
DATA TWO/ 2.0D0, 4.0D0, 8.0D0, 16.0D0, 32.0D0,
& 64.0D0, 128.0D0, 256.0D0, 512.0D0,1024.0D0,
& 2048.0D0,4096.0D0,8192.0D0/
do L=1,8
Y(L) = P(L) + H * G(K+1) * (F(L)-PHI(L,1))
enddo
IFAIL = 0
KOLD = K
HOLD = H
*** Derivative
CALL FUNCT(Y,F,1)
do L=1,8
PHI(L,K+1) = F(L)-PHI(L,1)
PHI(L,K+2) = PHI(L,K+1)-PHI(L,K+2)
enddo
do I=1,K
do L=1,8
PHI(L,I) = PHI(L,I)+PHI(L,K+1)
enddo
enddo
ERKP1 = 0.0D0
IF(KNEW.EQ.K-1.OR.K.EQ.12) IPHASE=0
IF(IPHASE.EQ.1) THEN
K = K+1
ERK = ERKP1
GOTO 500
ENDIF
IF(KNEW.EQ.K-1) THEN
K = K-1
ERK = ERKM1
GOTO 500
ENDIF
IF(K+1.LE.KS) THEN
do L=1,8
ERKP1 = ERKP1+(PHI(L,K+2)/WT(L))**2
enddo
ERKP1 = DABS(H)*GSTR(K+1)*DSQRT(ERKP1)
IF(K.LE.1) THEN
IF(ERKP1.LT.0.5D0*ERK) THEN
K = K+1
ERK = ERKP1
ENDIF
GOTO 500
ENDIF
IF(ERKM1.LE.DMIN1(ERK,ERKP1)) THEN
K = K-1
ERK = ERKM1
ELSE IF(ERKP1.LT.ERK .AND. K.NE.12) THEN
K = K+1
ERK = ERKP1
ENDIF
ENDIF
500 CONTINUE
*** H ->>
IF (ERK*TWO(K+1).LT.0.5D0*ERR0 .OR. IPHASE.EQ.1) THEN
HNEW = 1.01D0*H
GOTO 1000
ENDIF
*** <<- H
IF (ERK.GT.0.5D0*ERR0) THEN
R = (0.5D0*ERR0/ERK)**(1.0D0/DBLE(K+1))
HNEW = DMIN1(0.9D0,DMAX1(0.5D0,R)) * H
GOTO 1000
ENDIF
*** -> H <-
HNEW = H
1000 CONTINUE
H = HNEW
WRITE(25,800) Y(1)-6370.0, HNEW, ERK
800 FORMAT(F18.2, F18.6, F18.8)
RETURN
END
************************************************************************
SUBROUTINE ADAM_RETRY(ERR0)
************************************************************************
IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
COMMON /QADAM /Y(8),P(8),F(8),
& H,HOLD,WT(8),PHI(8,14),PSIOLD(12),
& PSI(12),BETA(12),SIGMA(13),G(13),
& ERKM1,ERK,K,KOLD,
& KNEW,KS,IPHASE,IFAIL
IPHASE = 0
do I=1,K
do L=1,8
PHI(L,I) = (PHI(L,I)-PHI(L,I+1))/BETA(I)
enddo
enddo
do I = 2,K
PSI(I-1) = PSI(I)-H
enddo
IFAIL = IFAIL+1
* Euler method
IF (IFAIL.GE.3) KNEW = 1
K = KNEW
IF (IFAIL.GT.3 .AND. 0.5D0*ERR0.LT.0.25D0*ERK) THEN
H = DSQRT(0.5D0*ERR0/ERK) * H
ELSE
H = 0.5D0 * H
ENDIF
*** Derivative
CALL FUNCT(Y,F,1)
RETURN
END