qpms/amos/zbknu.f

569 lines
17 KiB
FortranFixed
Raw Permalink Normal View History

SUBROUTINE ZBKNU(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM,
* ALIM)
C***BEGIN PROLOGUE ZBKNU
C***REFER TO ZBESI,ZBESK,ZAIRY,ZBESH
C
C ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE.
C
C***ROUTINES CALLED DGAMLN,I1MACH,D1MACH,ZKSCL,ZSHCH,ZUCHK,AZABS,ZDIV,
C AZEXP,AZLOG,ZMLT,AZSQRT
C***END PROLOGUE ZBKNU
C
DOUBLE PRECISION AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ,
* CBI, CBR, CC, CCHI, CCHR, CKI, CKR, COEFI, COEFR, CONEI, CONER,
* CRSCR, CSCLR, CSHI, CSHR, CSI, CSR, CSRR, CSSR, CTWOR,
* CZEROI, CZEROR, CZI, CZR, DNU, DNU2, DPI, ELIM, ETEST, FC, FHS,
* FI, FK, FKS, FMUI, FMUR, FNU, FPI, FR, G1, G2, HPI, PI, PR, PTI,
* PTR, P1I, P1R, P2I, P2M, P2R, QI, QR, RAK, RCAZ, RTHPI, RZI,
* RZR, R1, S, SMUI, SMUR, SPI, STI, STR, S1I, S1R, S2I, S2R, TM,
* TOL, TTH, T1, T2, YI, YR, ZI, ZR, DGAMLN, D1MACH, AZABS, ELM,
* CELMR, ZDR, ZDI, AS, ALAS, HELIM, CYR, CYI
INTEGER I, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N, NZ,
* IDUM, I1MACH, J, IC, INUB, NW
DIMENSION YR(N), YI(N), CC(8), CSSR(3), CSRR(3), BRY(3), CYR(2),
* CYI(2)
C COMPLEX Z,Y,A,B,RZ,SMU,FU,FMU,F,FLRZ,CZ,S1,S2,CSH,CCH
C COMPLEX CK,P,Q,COEF,P1,P2,CBK,PT,CZERO,CONE,CTWO,ST,EZ,CS,DK
C
DATA KMAX / 30 /
DATA CZEROR,CZEROI,CONER,CONEI,CTWOR,R1/
1 0.0D0 , 0.0D0 , 1.0D0 , 0.0D0 , 2.0D0 , 2.0D0 /
DATA DPI, RTHPI, SPI ,HPI, FPI, TTH /
1 3.14159265358979324D0, 1.25331413731550025D0,
2 1.90985931710274403D0, 1.57079632679489662D0,
3 1.89769999331517738D0, 6.66666666666666666D-01/
DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)/
1 5.77215664901532861D-01, -4.20026350340952355D-02,
2 -4.21977345555443367D-02, 7.21894324666309954D-03,
3 -2.15241674114950973D-04, -2.01348547807882387D-05,
4 1.13302723198169588D-06, 6.11609510448141582D-09/
C
CAZ = AZABS(ZR,ZI)
CSCLR = 1.0D0/TOL
CRSCR = TOL
CSSR(1) = CSCLR
CSSR(2) = 1.0D0
CSSR(3) = CRSCR
CSRR(1) = CRSCR
CSRR(2) = 1.0D0
CSRR(3) = CSCLR
BRY(1) = 1.0D+3*D1MACH(1)/TOL
BRY(2) = 1.0D0/BRY(1)
BRY(3) = D1MACH(2)
NZ = 0
IFLAG = 0
KODED = KODE
RCAZ = 1.0D0/CAZ
STR = ZR*RCAZ
STI = -ZI*RCAZ
RZR = (STR+STR)*RCAZ
RZI = (STI+STI)*RCAZ
INU = INT(SNGL(FNU+0.5D0))
DNU = FNU - DBLE(FLOAT(INU))
IF (DABS(DNU).EQ.0.5D0) GO TO 110
DNU2 = 0.0D0
IF (DABS(DNU).GT.TOL) DNU2 = DNU*DNU
IF (CAZ.GT.R1) GO TO 110
C-----------------------------------------------------------------------
C SERIES FOR CABS(Z).LE.R1
C-----------------------------------------------------------------------
FC = 1.0D0
CALL AZLOG(RZR, RZI, SMUR, SMUI, IDUM)
FMUR = SMUR*DNU
FMUI = SMUI*DNU
CALL ZSHCH(FMUR, FMUI, CSHR, CSHI, CCHR, CCHI)
IF (DNU.EQ.0.0D0) GO TO 10
FC = DNU*DPI
FC = FC/DSIN(FC)
SMUR = CSHR/DNU
SMUI = CSHI/DNU
10 CONTINUE
A2 = 1.0D0 + DNU
C-----------------------------------------------------------------------
C GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
C-----------------------------------------------------------------------
T2 = DEXP(-DGAMLN(A2,IDUM))
T1 = 1.0D0/(T2*FC)
IF (DABS(DNU).GT.0.1D0) GO TO 40
C-----------------------------------------------------------------------
C SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
C-----------------------------------------------------------------------
AK = 1.0D0
S = CC(1)
DO 20 K=2,8
AK = AK*DNU2
TM = CC(K)*AK
S = S + TM
IF (DABS(TM).LT.TOL) GO TO 30
20 CONTINUE
30 G1 = -S
GO TO 50
40 CONTINUE
G1 = (T1-T2)/(DNU+DNU)
50 CONTINUE
G2 = (T1+T2)*0.5D0
FR = FC*(CCHR*G1+SMUR*G2)
FI = FC*(CCHI*G1+SMUI*G2)
CALL AZEXP(FMUR, FMUI, STR, STI)
PR = 0.5D0*STR/T2
PI = 0.5D0*STI/T2
CALL ZDIV(0.5D0, 0.0D0, STR, STI, PTR, PTI)
QR = PTR/T1
QI = PTI/T1
S1R = FR
S1I = FI
S2R = PR
S2I = PI
AK = 1.0D0
A1 = 1.0D0
CKR = CONER
CKI = CONEI
BK = 1.0D0 - DNU2
IF (INU.GT.0 .OR. N.GT.1) GO TO 80
C-----------------------------------------------------------------------
C GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1
C-----------------------------------------------------------------------
IF (CAZ.LT.TOL) GO TO 70
CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI)
CZR = 0.25D0*CZR
CZI = 0.25D0*CZI
T1 = 0.25D0*CAZ*CAZ
60 CONTINUE
FR = (FR*AK+PR+QR)/BK
FI = (FI*AK+PI+QI)/BK
STR = 1.0D0/(AK-DNU)
PR = PR*STR
PI = PI*STR
STR = 1.0D0/(AK+DNU)
QR = QR*STR
QI = QI*STR
STR = CKR*CZR - CKI*CZI
RAK = 1.0D0/AK
CKI = (CKR*CZI+CKI*CZR)*RAK
CKR = STR*RAK
S1R = CKR*FR - CKI*FI + S1R
S1I = CKR*FI + CKI*FR + S1I
A1 = A1*T1*RAK
BK = BK + AK + AK + 1.0D0
AK = AK + 1.0D0
IF (A1.GT.TOL) GO TO 60
70 CONTINUE
YR(1) = S1R
YI(1) = S1I
IF (KODED.EQ.1) RETURN
CALL AZEXP(ZR, ZI, STR, STI)
CALL ZMLT(S1R, S1I, STR, STI, YR(1), YI(1))
RETURN
C-----------------------------------------------------------------------
C GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
C-----------------------------------------------------------------------
80 CONTINUE
IF (CAZ.LT.TOL) GO TO 100
CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI)
CZR = 0.25D0*CZR
CZI = 0.25D0*CZI
T1 = 0.25D0*CAZ*CAZ
90 CONTINUE
FR = (FR*AK+PR+QR)/BK
FI = (FI*AK+PI+QI)/BK
STR = 1.0D0/(AK-DNU)
PR = PR*STR
PI = PI*STR
STR = 1.0D0/(AK+DNU)
QR = QR*STR
QI = QI*STR
STR = CKR*CZR - CKI*CZI
RAK = 1.0D0/AK
CKI = (CKR*CZI+CKI*CZR)*RAK
CKR = STR*RAK
S1R = CKR*FR - CKI*FI + S1R
S1I = CKR*FI + CKI*FR + S1I
STR = PR - FR*AK
STI = PI - FI*AK
S2R = CKR*STR - CKI*STI + S2R
S2I = CKR*STI + CKI*STR + S2I
A1 = A1*T1*RAK
BK = BK + AK + AK + 1.0D0
AK = AK + 1.0D0
IF (A1.GT.TOL) GO TO 90
100 CONTINUE
KFLAG = 2
A1 = FNU + 1.0D0
AK = A1*DABS(SMUR)
IF (AK.GT.ALIM) KFLAG = 3
STR = CSSR(KFLAG)
P2R = S2R*STR
P2I = S2I*STR
CALL ZMLT(P2R, P2I, RZR, RZI, S2R, S2I)
S1R = S1R*STR
S1I = S1I*STR
IF (KODED.EQ.1) GO TO 210
CALL AZEXP(ZR, ZI, FR, FI)
CALL ZMLT(S1R, S1I, FR, FI, S1R, S1I)
CALL ZMLT(S2R, S2I, FR, FI, S2R, S2I)
GO TO 210
C-----------------------------------------------------------------------
C IFLAG=0 MEANS NO UNDERFLOW OCCURRED
C IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
C RECURSION
C-----------------------------------------------------------------------
110 CONTINUE
CALL AZSQRT(ZR, ZI, STR, STI)
CALL ZDIV(RTHPI, CZEROI, STR, STI, COEFR, COEFI)
KFLAG = 2
IF (KODED.EQ.2) GO TO 120
IF (ZR.GT.ALIM) GO TO 290
C BLANK LINE
STR = DEXP(-ZR)*CSSR(KFLAG)
STI = -STR*DSIN(ZI)
STR = STR*DCOS(ZI)
CALL ZMLT(COEFR, COEFI, STR, STI, COEFR, COEFI)
120 CONTINUE
IF (DABS(DNU).EQ.0.5D0) GO TO 300
C-----------------------------------------------------------------------
C MILLER ALGORITHM FOR CABS(Z).GT.R1
C-----------------------------------------------------------------------
AK = DCOS(DPI*DNU)
AK = DABS(AK)
IF (AK.EQ.CZEROR) GO TO 300
FHS = DABS(0.25D0-DNU2)
IF (FHS.EQ.CZEROR) GO TO 300
C-----------------------------------------------------------------------
C COMPUTE R2=F(E). IF CABS(Z).GE.R2, USE FORWARD RECURRENCE TO
C DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
C 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(14))=
C TOL WHERE B IS THE BASE OF THE ARITHMETIC.
C-----------------------------------------------------------------------
T1 = DBLE(FLOAT(I1MACH(14)-1))
T1 = T1*D1MACH(5)*3.321928094D0
T1 = DMAX1(T1,12.0D0)
T1 = DMIN1(T1,60.0D0)
T2 = TTH*T1 - 6.0D0
IF (ZR.NE.0.0D0) GO TO 130
T1 = HPI
GO TO 140
130 CONTINUE
T1 = DATAN(ZI/ZR)
T1 = DABS(T1)
140 CONTINUE
IF (T2.GT.CAZ) GO TO 170
C-----------------------------------------------------------------------
C FORWARD RECURRENCE LOOP WHEN CABS(Z).GE.R2
C-----------------------------------------------------------------------
ETEST = AK/(DPI*CAZ*TOL)
FK = CONER
IF (ETEST.LT.CONER) GO TO 180
FKS = CTWOR
CKR = CAZ + CAZ + CTWOR
P1R = CZEROR
P2R = CONER
DO 150 I=1,KMAX
AK = FHS/FKS
CBR = CKR/(FK+CONER)
PTR = P2R
P2R = CBR*P2R - P1R*AK
P1R = PTR
CKR = CKR + CTWOR
FKS = FKS + FK + FK + CTWOR
FHS = FHS + FK + FK
FK = FK + CONER
STR = DABS(P2R)*FK
IF (ETEST.LT.STR) GO TO 160
150 CONTINUE
GO TO 310
160 CONTINUE
FK = FK + SPI*T1*DSQRT(T2/CAZ)
FHS = DABS(0.25D0-DNU2)
GO TO 180
170 CONTINUE
C-----------------------------------------------------------------------
C COMPUTE BACKWARD INDEX K FOR CABS(Z).LT.R2
C-----------------------------------------------------------------------
A2 = DSQRT(CAZ)
AK = FPI*AK/(TOL*DSQRT(A2))
AA = 3.0D0*T1/(1.0D0+CAZ)
BB = 14.7D0*T1/(28.0D0+CAZ)
AK = (DLOG(AK)+CAZ*DCOS(AA)/(1.0D0+0.008D0*CAZ))/DCOS(BB)
FK = 0.12125D0*AK*AK/CAZ + 1.5D0
180 CONTINUE
C-----------------------------------------------------------------------
C BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
C-----------------------------------------------------------------------
K = INT(SNGL(FK))
FK = DBLE(FLOAT(K))
FKS = FK*FK
P1R = CZEROR
P1I = CZEROI
P2R = TOL
P2I = CZEROI
CSR = P2R
CSI = P2I
DO 190 I=1,K
A1 = FKS - FK
AK = (FKS+FK)/(A1+FHS)
RAK = 2.0D0/(FK+CONER)
CBR = (FK+ZR)*RAK
CBI = ZI*RAK
PTR = P2R
PTI = P2I
P2R = (PTR*CBR-PTI*CBI-P1R)*AK
P2I = (PTI*CBR+PTR*CBI-P1I)*AK
P1R = PTR
P1I = PTI
CSR = CSR + P2R
CSI = CSI + P2I
FKS = A1 - FK + CONER
FK = FK - CONER
190 CONTINUE
C-----------------------------------------------------------------------
C COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER
C SCALING
C-----------------------------------------------------------------------
TM = AZABS(CSR,CSI)
PTR = 1.0D0/TM
S1R = P2R*PTR
S1I = P2I*PTR
CSR = CSR*PTR
CSI = -CSI*PTR
CALL ZMLT(COEFR, COEFI, S1R, S1I, STR, STI)
CALL ZMLT(STR, STI, CSR, CSI, S1R, S1I)
IF (INU.GT.0 .OR. N.GT.1) GO TO 200
ZDR = ZR
ZDI = ZI
IF(IFLAG.EQ.1) GO TO 270
GO TO 240
200 CONTINUE
C-----------------------------------------------------------------------
C COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING
C-----------------------------------------------------------------------
TM = AZABS(P2R,P2I)
PTR = 1.0D0/TM
P1R = P1R*PTR
P1I = P1I*PTR
P2R = P2R*PTR
P2I = -P2I*PTR
CALL ZMLT(P1R, P1I, P2R, P2I, PTR, PTI)
STR = DNU + 0.5D0 - PTR
STI = -PTI
CALL ZDIV(STR, STI, ZR, ZI, STR, STI)
STR = STR + 1.0D0
CALL ZMLT(STR, STI, S1R, S1I, S2R, S2I)
C-----------------------------------------------------------------------
C FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH
C SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
C-----------------------------------------------------------------------
210 CONTINUE
STR = DNU + 1.0D0
CKR = STR*RZR
CKI = STR*RZI
IF (N.EQ.1) INU = INU - 1
IF (INU.GT.0) GO TO 220
IF (N.GT.1) GO TO 215
S1R = S2R
S1I = S2I
215 CONTINUE
ZDR = ZR
ZDI = ZI
IF(IFLAG.EQ.1) GO TO 270
GO TO 240
220 CONTINUE
INUB = 1
IF(IFLAG.EQ.1) GO TO 261
225 CONTINUE
P1R = CSRR(KFLAG)
ASCLE = BRY(KFLAG)
DO 230 I=INUB,INU
STR = S2R
STI = S2I
S2R = CKR*STR - CKI*STI + S1R
S2I = CKR*STI + CKI*STR + S1I
S1R = STR
S1I = STI
CKR = CKR + RZR
CKI = CKI + RZI
IF (KFLAG.GE.3) GO TO 230
P2R = S2R*P1R
P2I = S2I*P1R
STR = DABS(P2R)
STI = DABS(P2I)
P2M = DMAX1(STR,STI)
IF (P2M.LE.ASCLE) GO TO 230
KFLAG = KFLAG + 1
ASCLE = BRY(KFLAG)
S1R = S1R*P1R
S1I = S1I*P1R
S2R = P2R
S2I = P2I
STR = CSSR(KFLAG)
S1R = S1R*STR
S1I = S1I*STR
S2R = S2R*STR
S2I = S2I*STR
P1R = CSRR(KFLAG)
230 CONTINUE
IF (N.NE.1) GO TO 240
S1R = S2R
S1I = S2I
240 CONTINUE
STR = CSRR(KFLAG)
YR(1) = S1R*STR
YI(1) = S1I*STR
IF (N.EQ.1) RETURN
YR(2) = S2R*STR
YI(2) = S2I*STR
IF (N.EQ.2) RETURN
KK = 2
250 CONTINUE
KK = KK + 1
IF (KK.GT.N) RETURN
P1R = CSRR(KFLAG)
ASCLE = BRY(KFLAG)
DO 260 I=KK,N
P2R = S2R
P2I = S2I
S2R = CKR*P2R - CKI*P2I + S1R
S2I = CKI*P2R + CKR*P2I + S1I
S1R = P2R
S1I = P2I
CKR = CKR + RZR
CKI = CKI + RZI
P2R = S2R*P1R
P2I = S2I*P1R
YR(I) = P2R
YI(I) = P2I
IF (KFLAG.GE.3) GO TO 260
STR = DABS(P2R)
STI = DABS(P2I)
P2M = DMAX1(STR,STI)
IF (P2M.LE.ASCLE) GO TO 260
KFLAG = KFLAG + 1
ASCLE = BRY(KFLAG)
S1R = S1R*P1R
S1I = S1I*P1R
S2R = P2R
S2I = P2I
STR = CSSR(KFLAG)
S1R = S1R*STR
S1I = S1I*STR
S2R = S2R*STR
S2I = S2I*STR
P1R = CSRR(KFLAG)
260 CONTINUE
RETURN
C-----------------------------------------------------------------------
C IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
C-----------------------------------------------------------------------
261 CONTINUE
HELIM = 0.5D0*ELIM
ELM = DEXP(-ELIM)
CELMR = ELM
ASCLE = BRY(1)
ZDR = ZR
ZDI = ZI
IC = -1
J = 2
DO 262 I=1,INU
STR = S2R
STI = S2I
S2R = STR*CKR-STI*CKI+S1R
S2I = STI*CKR+STR*CKI+S1I
S1R = STR
S1I = STI
CKR = CKR+RZR
CKI = CKI+RZI
AS = AZABS(S2R,S2I)
ALAS = DLOG(AS)
P2R = -ZDR+ALAS
IF(P2R.LT.(-ELIM)) GO TO 263
CALL AZLOG(S2R,S2I,STR,STI,IDUM)
P2R = -ZDR+STR
P2I = -ZDI+STI
P2M = DEXP(P2R)/TOL
P1R = P2M*DCOS(P2I)
P1I = P2M*DSIN(P2I)
CALL ZUCHK(P1R,P1I,NW,ASCLE,TOL)
IF(NW.NE.0) GO TO 263
J = 3 - J
CYR(J) = P1R
CYI(J) = P1I
IF(IC.EQ.(I-1)) GO TO 264
IC = I
GO TO 262
263 CONTINUE
IF(ALAS.LT.HELIM) GO TO 262
ZDR = ZDR-ELIM
S1R = S1R*CELMR
S1I = S1I*CELMR
S2R = S2R*CELMR
S2I = S2I*CELMR
262 CONTINUE
IF(N.NE.1) GO TO 270
S1R = S2R
S1I = S2I
GO TO 270
264 CONTINUE
KFLAG = 1
INUB = I+1
S2R = CYR(J)
S2I = CYI(J)
J = 3 - J
S1R = CYR(J)
S1I = CYI(J)
IF(INUB.LE.INU) GO TO 225
IF(N.NE.1) GO TO 240
S1R = S2R
S1I = S2I
GO TO 240
270 CONTINUE
YR(1) = S1R
YI(1) = S1I
IF(N.EQ.1) GO TO 280
YR(2) = S2R
YI(2) = S2I
280 CONTINUE
ASCLE = BRY(1)
CALL ZKSCL(ZDR,ZDI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM)
INU = N - NZ
IF (INU.LE.0) RETURN
KK = NZ + 1
S1R = YR(KK)
S1I = YI(KK)
YR(KK) = S1R*CSRR(1)
YI(KK) = S1I*CSRR(1)
IF (INU.EQ.1) RETURN
KK = NZ + 2
S2R = YR(KK)
S2I = YI(KK)
YR(KK) = S2R*CSRR(1)
YI(KK) = S2I*CSRR(1)
IF (INU.EQ.2) RETURN
T2 = FNU + DBLE(FLOAT(KK-1))
CKR = T2*RZR
CKI = T2*RZI
KFLAG = 1
GO TO 250
290 CONTINUE
C-----------------------------------------------------------------------
C SCALE BY DEXP(Z), IFLAG = 1 CASES
C-----------------------------------------------------------------------
KODED = 2
IFLAG = 1
KFLAG = 2
GO TO 120
C-----------------------------------------------------------------------
C FNU=HALF ODD INTEGER CASE, DNU=-0.5
C-----------------------------------------------------------------------
300 CONTINUE
S1R = COEFR
S1I = COEFI
S2R = COEFR
S2I = COEFI
GO TO 210
C
C
310 CONTINUE
NZ=-2
RETURN
END