*DECK ZBKNU SUBROUTINE ZBKNU (ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, + ALIM) C***BEGIN PROLOGUE ZBKNU C***SUBSIDIARY C***PURPOSE Subsidiary to ZAIRY, ZBESH, ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CBKNU-A, ZBKNU-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE. C C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESK C***ROUTINES CALLED D1MACH, DGAMLN, I1MACH, ZABS, ZDIV, ZEXP, ZKSCL, C ZLOG, ZMLT, ZSHCH, ZSQRT, ZUCHK C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C 930122 Added ZEXP, ZLOG and ZSQRT to EXTERNAL statement. (RWC) 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, ZABS, 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) EXTERNAL ZABS, ZEXP, ZLOG, ZSQRT 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***FIRST EXECUTABLE STATEMENT ZBKNU CAZ = ZABS(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 = FNU+0.5D0 DNU = FNU - INU IF (ABS(DNU).EQ.0.5D0) GO TO 110 DNU2 = 0.0D0 IF (ABS(DNU).GT.TOL) DNU2 = DNU*DNU IF (CAZ.GT.R1) GO TO 110 C----------------------------------------------------------------------- C SERIES FOR ABS(Z).LE.R1 C----------------------------------------------------------------------- FC = 1.0D0 CALL ZLOG(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/SIN(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 = EXP(-DGAMLN(A2,IDUM)) T1 = 1.0D0/(T2*FC) IF (ABS(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 (ABS(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 ZEXP(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 ZEXP(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*ABS(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 ZEXP(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 ZSQRT(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 = EXP(-ZR)*CSSR(KFLAG) STI = -STR*SIN(ZI) STR = STR*COS(ZI) CALL ZMLT(COEFR, COEFI, STR, STI, COEFR, COEFI) 120 CONTINUE IF (ABS(DNU).EQ.0.5D0) GO TO 300 C----------------------------------------------------------------------- C MILLER ALGORITHM FOR ABS(Z).GT.R1 C----------------------------------------------------------------------- AK = COS(DPI*DNU) AK = ABS(AK) IF (AK.EQ.CZEROR) GO TO 300 FHS = ABS(0.25D0-DNU2) IF (FHS.EQ.CZEROR) GO TO 300 C----------------------------------------------------------------------- C COMPUTE R2=F(E). IF ABS(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 = I1MACH(14)-1 T1 = T1*D1MACH(5)*3.321928094D0 T1 = MAX(T1,12.0D0) T1 = MIN(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 = ABS(T1) 140 CONTINUE IF (T2.GT.CAZ) GO TO 170 C----------------------------------------------------------------------- C FORWARD RECURRENCE LOOP WHEN ABS(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 = ABS(P2R)*FK IF (ETEST.LT.STR) GO TO 160 150 CONTINUE GO TO 310 160 CONTINUE FK = FK + SPI*T1*SQRT(T2/CAZ) FHS = ABS(0.25D0-DNU2) GO TO 180 170 CONTINUE C----------------------------------------------------------------------- C COMPUTE BACKWARD INDEX K FOR ABS(Z).LT.R2 C----------------------------------------------------------------------- A2 = SQRT(CAZ) AK = FPI*AK/(TOL*SQRT(A2)) AA = 3.0D0*T1/(1.0D0+CAZ) BB = 14.7D0*T1/(28.0D0+CAZ) AK = (LOG(AK)+CAZ*COS(AA)/(1.0D0+0.008D0*CAZ))/COS(BB) FK = 0.12125D0*AK*AK/CAZ + 1.5D0 180 CONTINUE C----------------------------------------------------------------------- C BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM C----------------------------------------------------------------------- K = FK FK = 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/ABS(CS))*(CONJG(CS)/ABS(CS)) FOR BETTER C SCALING C----------------------------------------------------------------------- TM = ZABS(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/ABS(P2)*CONJG(P2)/ABS(P2) FOR SCALING C----------------------------------------------------------------------- TM = ZABS(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 = ABS(P2R) STI = ABS(P2I) P2M = MAX(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 = ABS(P2R) STI = ABS(P2I) P2M = MAX(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 = EXP(-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 = ZABS(S2R,S2I) ALAS = LOG(AS) P2R = -ZDR+ALAS IF(P2R.LT.(-ELIM)) GO TO 263 CALL ZLOG(S2R,S2I,STR,STI,IDUM) P2R = -ZDR+STR P2I = -ZDI+STI P2M = EXP(P2R)/TOL P1R = P2M*COS(P2I) P1I = P2M*SIN(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 + (KK-1) CKR = T2*RZR CKI = T2*RZI KFLAG = 1 GO TO 250 290 CONTINUE C----------------------------------------------------------------------- C SCALE BY EXP(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