qpms/amos/zmlri.f

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,ZABS,ZEXP,ZLOG,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, ZABS
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 = ZABS(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 = ZABS(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 = ZABS(P2R,P2I)
IF (AP.LT.TST) GO TO 30
IF (ITIME.EQ.2) GO TO 40
ACK = ZABS(CKR,CKI)
FLAM = ACK + DSQRT(ACK*ACK-1.0D0)
FKAP = AP/ZABS(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 ZLOG(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 = ZABS(P2R,P2I)
P1R = 1.0D0/AP
CALL ZEXP(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