205 lines
5.9 KiB
Fortran
205 lines
5.9 KiB
Fortran
SUBROUTINE ZMLRI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)
|
|
C***BEGIN PROLOGUE ZMLRI
|
|
C***REFER TO ZBESI,ZBESK
|
|
C
|
|
C ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
|
|
C MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
|
|
C
|
|
C***ROUTINES CALLED DGAMLN,D1MACH,AZABS,AZEXP,AZLOG,ZMLT
|
|
C***END PROLOGUE ZMLRI
|
|
C COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z
|
|
DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI,
|
|
* CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, FNU, PTI, PTR, P1I,
|
|
* P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI,
|
|
* SUMR, TFNF, TOL, TST, YI, YR, ZEROI, ZEROR, ZI, ZR, DGAMLN,
|
|
* D1MACH, AZABS
|
|
INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ
|
|
DIMENSION YR(N), YI(N)
|
|
DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
|
|
SCLE = D1MACH(1)/TOL
|
|
NZ=0
|
|
AZ = AZABS(ZR,ZI)
|
|
IAZ = INT(SNGL(AZ))
|
|
IFNU = INT(SNGL(FNU))
|
|
INU = IFNU + N - 1
|
|
AT = DBLE(FLOAT(IAZ)) + 1.0D0
|
|
RAZ = 1.0D0/AZ
|
|
STR = ZR*RAZ
|
|
STI = -ZI*RAZ
|
|
CKR = STR*AT*RAZ
|
|
CKI = STI*AT*RAZ
|
|
RZR = (STR+STR)*RAZ
|
|
RZI = (STI+STI)*RAZ
|
|
P1R = ZEROR
|
|
P1I = ZEROI
|
|
P2R = CONER
|
|
P2I = CONEI
|
|
ACK = (AT+1.0D0)*RAZ
|
|
RHO = ACK + DSQRT(ACK*ACK-1.0D0)
|
|
RHO2 = RHO*RHO
|
|
TST = (RHO2+RHO2)/((RHO2-1.0D0)*(RHO-1.0D0))
|
|
TST = TST/TOL
|
|
C-----------------------------------------------------------------------
|
|
C COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
|
|
C-----------------------------------------------------------------------
|
|
AK = AT
|
|
DO 10 I=1,80
|
|
PTR = P2R
|
|
PTI = P2I
|
|
P2R = P1R - (CKR*PTR-CKI*PTI)
|
|
P2I = P1I - (CKI*PTR+CKR*PTI)
|
|
P1R = PTR
|
|
P1I = PTI
|
|
CKR = CKR + RZR
|
|
CKI = CKI + RZI
|
|
AP = AZABS(P2R,P2I)
|
|
IF (AP.GT.TST*AK*AK) GO TO 20
|
|
AK = AK + 1.0D0
|
|
10 CONTINUE
|
|
GO TO 110
|
|
20 CONTINUE
|
|
I = I + 1
|
|
K = 0
|
|
IF (INU.LT.IAZ) GO TO 40
|
|
C-----------------------------------------------------------------------
|
|
C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
|
|
C-----------------------------------------------------------------------
|
|
P1R = ZEROR
|
|
P1I = ZEROI
|
|
P2R = CONER
|
|
P2I = CONEI
|
|
AT = DBLE(FLOAT(INU)) + 1.0D0
|
|
STR = ZR*RAZ
|
|
STI = -ZI*RAZ
|
|
CKR = STR*AT*RAZ
|
|
CKI = STI*AT*RAZ
|
|
ACK = AT*RAZ
|
|
TST = DSQRT(ACK/TOL)
|
|
ITIME = 1
|
|
DO 30 K=1,80
|
|
PTR = P2R
|
|
PTI = P2I
|
|
P2R = P1R - (CKR*PTR-CKI*PTI)
|
|
P2I = P1I - (CKR*PTI+CKI*PTR)
|
|
P1R = PTR
|
|
P1I = PTI
|
|
CKR = CKR + RZR
|
|
CKI = CKI + RZI
|
|
AP = AZABS(P2R,P2I)
|
|
IF (AP.LT.TST) GO TO 30
|
|
IF (ITIME.EQ.2) GO TO 40
|
|
ACK = AZABS(CKR,CKI)
|
|
FLAM = ACK + DSQRT(ACK*ACK-1.0D0)
|
|
FKAP = AP/AZABS(P1R,P1I)
|
|
RHO = DMIN1(FLAM,FKAP)
|
|
TST = TST*DSQRT(RHO/(RHO*RHO-1.0D0))
|
|
ITIME = 2
|
|
30 CONTINUE
|
|
GO TO 110
|
|
40 CONTINUE
|
|
C-----------------------------------------------------------------------
|
|
C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
|
|
C-----------------------------------------------------------------------
|
|
K = K + 1
|
|
KK = MAX0(I+IAZ,K+INU)
|
|
FKK = DBLE(FLOAT(KK))
|
|
P1R = ZEROR
|
|
P1I = ZEROI
|
|
C-----------------------------------------------------------------------
|
|
C SCALE P2 AND SUM BY SCLE
|
|
C-----------------------------------------------------------------------
|
|
P2R = SCLE
|
|
P2I = ZEROI
|
|
FNF = FNU - DBLE(FLOAT(IFNU))
|
|
TFNF = FNF + FNF
|
|
BK = DGAMLN(FKK+TFNF+1.0D0,IDUM) - DGAMLN(FKK+1.0D0,IDUM) -
|
|
* DGAMLN(TFNF+1.0D0,IDUM)
|
|
BK = DEXP(BK)
|
|
SUMR = ZEROR
|
|
SUMI = ZEROI
|
|
KM = KK - INU
|
|
DO 50 I=1,KM
|
|
PTR = P2R
|
|
PTI = P2I
|
|
P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
|
|
P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
|
|
P1R = PTR
|
|
P1I = PTI
|
|
AK = 1.0D0 - TFNF/(FKK+TFNF)
|
|
ACK = BK*AK
|
|
SUMR = SUMR + (ACK+BK)*P1R
|
|
SUMI = SUMI + (ACK+BK)*P1I
|
|
BK = ACK
|
|
FKK = FKK - 1.0D0
|
|
50 CONTINUE
|
|
YR(N) = P2R
|
|
YI(N) = P2I
|
|
IF (N.EQ.1) GO TO 70
|
|
DO 60 I=2,N
|
|
PTR = P2R
|
|
PTI = P2I
|
|
P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
|
|
P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
|
|
P1R = PTR
|
|
P1I = PTI
|
|
AK = 1.0D0 - TFNF/(FKK+TFNF)
|
|
ACK = BK*AK
|
|
SUMR = SUMR + (ACK+BK)*P1R
|
|
SUMI = SUMI + (ACK+BK)*P1I
|
|
BK = ACK
|
|
FKK = FKK - 1.0D0
|
|
M = N - I + 1
|
|
YR(M) = P2R
|
|
YI(M) = P2I
|
|
60 CONTINUE
|
|
70 CONTINUE
|
|
IF (IFNU.LE.0) GO TO 90
|
|
DO 80 I=1,IFNU
|
|
PTR = P2R
|
|
PTI = P2I
|
|
P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
|
|
P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR)
|
|
P1R = PTR
|
|
P1I = PTI
|
|
AK = 1.0D0 - TFNF/(FKK+TFNF)
|
|
ACK = BK*AK
|
|
SUMR = SUMR + (ACK+BK)*P1R
|
|
SUMI = SUMI + (ACK+BK)*P1I
|
|
BK = ACK
|
|
FKK = FKK - 1.0D0
|
|
80 CONTINUE
|
|
90 CONTINUE
|
|
PTR = ZR
|
|
PTI = ZI
|
|
IF (KODE.EQ.2) PTR = ZEROR
|
|
CALL AZLOG(RZR, RZI, STR, STI, IDUM)
|
|
P1R = -FNF*STR + PTR
|
|
P1I = -FNF*STI + PTI
|
|
AP = DGAMLN(1.0D0+FNF,IDUM)
|
|
PTR = P1R - AP
|
|
PTI = P1I
|
|
C-----------------------------------------------------------------------
|
|
C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
|
|
C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
|
|
C-----------------------------------------------------------------------
|
|
P2R = P2R + SUMR
|
|
P2I = P2I + SUMI
|
|
AP = AZABS(P2R,P2I)
|
|
P1R = 1.0D0/AP
|
|
CALL AZEXP(PTR, PTI, STR, STI)
|
|
CKR = STR*P1R
|
|
CKI = STI*P1R
|
|
PTR = P2R*P1R
|
|
PTI = -P2I*P1R
|
|
CALL ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI)
|
|
DO 100 I=1,N
|
|
STR = YR(I)*CNORMR - YI(I)*CNORMI
|
|
YI(I) = YR(I)*CNORMI + YI(I)*CNORMR
|
|
YR(I) = STR
|
|
100 CONTINUE
|
|
RETURN
|
|
110 CONTINUE
|
|
NZ=-2
|
|
RETURN
|
|
END
|