1 SUBROUTINE ZBKNU
(ZR
, ZI
, FNU
, KODE
, N
, YR
, YI
, NZ
, TOL
, ELIM
,
3 C***BEGIN PROLOGUE ZBKNU
4 C***REFER TO ZBESI,ZBESK,ZAIRY,ZBESH
6 C ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE.
8 C***ROUTINES CALLED DGAMLN,I1MACH,D1MACH,ZKSCL,ZSHCH,ZUCHK,AZABS,ZDIV,
9 C AZEXP,AZLOG,ZMLT,AZSQRT
10 C***END PROLOGUE ZBKNU
12 DOUBLE PRECISION AA
, AK
, ALIM
, ASCLE
, A1
, A2
, BB
, BK
, BRY
, CAZ
,
13 * CBI
, CBR
, CC
, CCHI
, CCHR
, CKI
, CKR
, COEFI
, COEFR
, CONEI
, CONER
,
14 * CRSCR
, CSCLR
, CSHI
, CSHR
, CSI
, CSR
, CSRR
, CSSR
, CTWOR
,
15 * CZEROI
, CZEROR
, CZI
, CZR
, DNU
, DNU2
, DPI
, ELIM
, ETEST
, FC
, FHS
,
16 * FI
, FK
, FKS
, FMUI
, FMUR
, FNU
, FPI
, FR
, G1
, G2
, HPI
, PI
, PR
, PTI
,
17 * PTR
, P1I
, P1R
, P2I
, P2M
, P2R
, QI
, QR
, RAK
, RCAZ
, RTHPI
, RZI
,
18 * RZR
, R1
, S
, SMUI
, SMUR
, SPI
, STI
, STR
, S1I
, S1R
, S2I
, S2R
, TM
,
19 * TOL
, TTH
, T1
, T2
, YI
, YR
, ZI
, ZR
, DGAMLN
, D1MACH
, AZABS
, ELM
,
20 * CELMR
, ZDR
, ZDI
, AS
, ALAS
, HELIM
, CYR
, CYI
21 INTEGER I
, IFLAG
, INU
, K
, KFLAG
, KK
, KMAX
, KODE
, KODED
, N
, NZ
,
22 * IDUM
, I1MACH
, J
, IC
, INUB
, NW
23 DIMENSION YR
(N
), YI
(N
), CC
(8), CSSR
(3), CSRR
(3), BRY
(3), CYR
(2),
25 C COMPLEX Z,Y,A,B,RZ,SMU,FU,FMU,F,FLRZ,CZ,S1,S2,CSH,CCH
26 C COMPLEX CK,P,Q,COEF,P1,P2,CBK,PT,CZERO,CONE,CTWO,ST,EZ,CS,DK
29 DATA CZEROR
,CZEROI
,CONER
,CONEI
,CTWOR
,R1
/
30 1 0.0D0
, 0.0D0
, 1.0D0
, 0.0D0
, 2.0D0
, 2.0D0
/
31 DATA DPI
, RTHPI
, SPI
,HPI
, FPI
, TTH
/
32 1 3.14159265358979324D0
, 1.25331413731550025D0
,
33 2 1.90985931710274403D0
, 1.57079632679489662D0
,
34 3 1.89769999331517738D0
, 6.66666666666666666D
-01/
35 DATA CC
(1), CC
(2), CC
(3), CC
(4), CC
(5), CC
(6), CC
(7), CC
(8)/
36 1 5.77215664901532861D
-01, -4.20026350340952355D
-02,
37 2 -4.21977345555443367D
-02, 7.21894324666309954D
-03,
38 3 -2.15241674114950973D
-04, -2.01348547807882387D
-05,
39 4 1.13302723198169588D
-06, 6.11609510448141582D
-09/
50 BRY
(1) = 1.0D
+3*D1MACH
(1)/TOL
61 INU
= INT
(SNGL
(FNU
+0.5D0
))
62 DNU
= FNU
- DBLE
(FLOAT
(INU
))
63 IF (DABS
(DNU
).EQ
.0.5D0
) GO TO 110
65 IF (DABS
(DNU
).GT
.TOL
) DNU2
= DNU*DNU
66 IF (CAZ
.GT
.R1
) GO TO 110
67 C-----------------------------------------------------------------------
68 C SERIES FOR CABS(Z).LE.R1
69 C-----------------------------------------------------------------------
71 CALL AZLOG
(RZR
, RZI
, SMUR
, SMUI
, IDUM
)
74 CALL ZSHCH
(FMUR
, FMUI
, CSHR
, CSHI
, CCHR
, CCHI
)
75 IF (DNU
.EQ
.0.0D0
) GO TO 10
82 C-----------------------------------------------------------------------
83 C GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
84 C-----------------------------------------------------------------------
85 T2
= DEXP
(-DGAMLN
(A2
,IDUM
))
87 IF (DABS
(DNU
).GT
.0.1D0
) GO TO 40
88 C-----------------------------------------------------------------------
89 C SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
90 C-----------------------------------------------------------------------
97 IF (DABS
(TM
).LT
.TOL
) GO TO 30
102 G1
= (T1
-T2
)/(DNU
+DNU
)
105 FR
= FC*
(CCHR*G1
+SMUR*G2
)
106 FI
= FC*
(CCHI*G1
+SMUI*G2
)
107 CALL AZEXP
(FMUR
, FMUI
, STR
, STI
)
110 CALL ZDIV
(0.5D0
, 0.0D0
, STR
, STI
, PTR
, PTI
)
122 IF (INU
.GT
.0 .OR
. N
.GT
.1) GO TO 80
123 C-----------------------------------------------------------------------
124 C GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1
125 C-----------------------------------------------------------------------
126 IF (CAZ
.LT
.TOL
) GO TO 70
127 CALL ZMLT
(ZR
, ZI
, ZR
, ZI
, CZR
, CZI
)
132 FR
= (FR*AK
+PR
+QR
)/BK
133 FI
= (FI*AK
+PI
+QI
)/BK
140 STR
= CKR*CZR
- CKI*CZI
142 CKI
= (CKR*CZI
+CKI*CZR
)*RAK
144 S1R
= CKR*FR
- CKI*FI
+ S1R
145 S1I
= CKR*FI
+ CKI*FR
+ S1I
147 BK
= BK
+ AK
+ AK
+ 1.0D0
149 IF (A1
.GT
.TOL
) GO TO 60
153 IF (KODED
.EQ
.1) RETURN
154 CALL AZEXP
(ZR
, ZI
, STR
, STI
)
155 CALL ZMLT
(S1R
, S1I
, STR
, STI
, YR
(1), YI
(1))
157 C-----------------------------------------------------------------------
158 C GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
159 C-----------------------------------------------------------------------
161 IF (CAZ
.LT
.TOL
) GO TO 100
162 CALL ZMLT
(ZR
, ZI
, ZR
, ZI
, CZR
, CZI
)
167 FR
= (FR*AK
+PR
+QR
)/BK
168 FI
= (FI*AK
+PI
+QI
)/BK
175 STR
= CKR*CZR
- CKI*CZI
177 CKI
= (CKR*CZI
+CKI*CZR
)*RAK
179 S1R
= CKR*FR
- CKI*FI
+ S1R
180 S1I
= CKR*FI
+ CKI*FR
+ S1I
183 S2R
= CKR*STR
- CKI*STI
+ S2R
184 S2I
= CKR*STI
+ CKI*STR
+ S2I
186 BK
= BK
+ AK
+ AK
+ 1.0D0
188 IF (A1
.GT
.TOL
) GO TO 90
193 IF (AK
.GT
.ALIM
) KFLAG
= 3
197 CALL ZMLT
(P2R
, P2I
, RZR
, RZI
, S2R
, S2I
)
200 IF (KODED
.EQ
.1) GO TO 210
201 CALL AZEXP
(ZR
, ZI
, FR
, FI
)
202 CALL ZMLT
(S1R
, S1I
, FR
, FI
, S1R
, S1I
)
203 CALL ZMLT
(S2R
, S2I
, FR
, FI
, S2R
, S2I
)
205 C-----------------------------------------------------------------------
206 C IFLAG=0 MEANS NO UNDERFLOW OCCURRED
207 C IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
208 C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
210 C-----------------------------------------------------------------------
212 CALL AZSQRT
(ZR
, ZI
, STR
, STI
)
213 CALL ZDIV
(RTHPI
, CZEROI
, STR
, STI
, COEFR
, COEFI
)
215 IF (KODED
.EQ
.2) GO TO 120
216 IF (ZR
.GT
.ALIM
) GO TO 290
218 STR
= DEXP
(-ZR
)*CSSR
(KFLAG
)
221 CALL ZMLT
(COEFR
, COEFI
, STR
, STI
, COEFR
, COEFI
)
223 IF (DABS
(DNU
).EQ
.0.5D0
) GO TO 300
224 C-----------------------------------------------------------------------
225 C MILLER ALGORITHM FOR CABS(Z).GT.R1
226 C-----------------------------------------------------------------------
229 IF (AK
.EQ
.CZEROR
) GO TO 300
230 FHS
= DABS
(0.25D0
-DNU2
)
231 IF (FHS
.EQ
.CZEROR
) GO TO 300
232 C-----------------------------------------------------------------------
233 C COMPUTE R2=F(E). IF CABS(Z).GE.R2, USE FORWARD RECURRENCE TO
234 C DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
235 C 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(14))=
236 C TOL WHERE B IS THE BASE OF THE ARITHMETIC.
237 C-----------------------------------------------------------------------
238 T1
= DBLE
(FLOAT
(I1MACH
(14)-1))
239 T1
= T1*D1MACH
(5)*3.321928094D0
240 T1
= DMAX1
(T1
,12.0D0
)
241 T1
= DMIN1
(T1
,60.0D0
)
243 IF (ZR
.NE
.0.0D0
) GO TO 130
250 IF (T2
.GT
.CAZ
) GO TO 170
251 C-----------------------------------------------------------------------
252 C FORWARD RECURRENCE LOOP WHEN CABS(Z).GE.R2
253 C-----------------------------------------------------------------------
254 ETEST
= AK
/(DPI*CAZ*TOL
)
256 IF (ETEST
.LT
.CONER
) GO TO 180
258 CKR
= CAZ
+ CAZ
+ CTWOR
265 P2R
= CBR*P2R
- P1R*AK
268 FKS
= FKS
+ FK
+ FK
+ CTWOR
272 IF (ETEST
.LT
.STR
) GO TO 160
276 FK
= FK
+ SPI*T1*DSQRT
(T2
/CAZ
)
277 FHS
= DABS
(0.25D0
-DNU2
)
280 C-----------------------------------------------------------------------
281 C COMPUTE BACKWARD INDEX K FOR CABS(Z).LT.R2
282 C-----------------------------------------------------------------------
284 AK
= FPI*AK
/(TOL*DSQRT
(A2
))
285 AA
= 3.0D0*T1
/(1.0D0
+CAZ
)
286 BB
= 14.7D0*T1
/(28.0D0
+CAZ
)
287 AK
= (DLOG
(AK
)+CAZ*DCOS
(AA
)/(1.0D0
+0.008D0*CAZ
))/DCOS
(BB
)
288 FK
= 0.12125D0*AK*AK
/CAZ
+ 1.5D0
290 C-----------------------------------------------------------------------
291 C BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
292 C-----------------------------------------------------------------------
304 AK
= (FKS
+FK
)/(A1
+FHS
)
305 RAK
= 2.0D0
/(FK
+CONER
)
310 P2R
= (PTR*CBR
-PTI*CBI
-P1R
)*AK
311 P2I
= (PTI*CBR
+PTR*CBI
-P1I
)*AK
316 FKS
= A1
- FK
+ CONER
319 C-----------------------------------------------------------------------
320 C COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER
322 C-----------------------------------------------------------------------
329 CALL ZMLT
(COEFR
, COEFI
, S1R
, S1I
, STR
, STI
)
330 CALL ZMLT
(STR
, STI
, CSR
, CSI
, S1R
, S1I
)
331 IF (INU
.GT
.0 .OR
. N
.GT
.1) GO TO 200
334 IF(IFLAG
.EQ
.1) GO TO 270
337 C-----------------------------------------------------------------------
338 C COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING
339 C-----------------------------------------------------------------------
346 CALL ZMLT
(P1R
, P1I
, P2R
, P2I
, PTR
, PTI
)
347 STR
= DNU
+ 0.5D0
- PTR
349 CALL ZDIV
(STR
, STI
, ZR
, ZI
, STR
, STI
)
351 CALL ZMLT
(STR
, STI
, S1R
, S1I
, S2R
, S2I
)
352 C-----------------------------------------------------------------------
353 C FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH
354 C SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
355 C-----------------------------------------------------------------------
360 IF (N
.EQ
.1) INU
= INU
- 1
361 IF (INU
.GT
.0) GO TO 220
362 IF (N
.GT
.1) GO TO 215
368 IF(IFLAG
.EQ
.1) GO TO 270
372 IF(IFLAG
.EQ
.1) GO TO 261
379 S2R
= CKR*STR
- CKI*STI
+ S1R
380 S2I
= CKR*STI
+ CKI*STR
+ S1I
385 IF (KFLAG
.GE
.3) GO TO 230
391 IF (P2M
.LE
.ASCLE
) GO TO 230
405 IF (N
.NE
.1) GO TO 240
425 S2R
= CKR*P2R
- CKI*P2I
+ S1R
426 S2I
= CKI*P2R
+ CKR*P2I
+ S1I
435 IF (KFLAG
.GE
.3) GO TO 260
439 IF (P2M
.LE
.ASCLE
) GO TO 260
454 C-----------------------------------------------------------------------
455 C IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
456 C-----------------------------------------------------------------------
469 S2R
= STR*CKR
-STI*CKI
+S1R
470 S2I
= STI*CKR
+STR*CKI
+S1I
478 IF(P2R
.LT
.(-ELIM
)) GO TO 263
479 CALL AZLOG
(S2R
,S2I
,STR
,STI
,IDUM
)
485 CALL ZUCHK
(P1R
,P1I
,NW
,ASCLE
,TOL
)
486 IF(NW
.NE
.0) GO TO 263
490 IF(IC
.EQ
.(I
-1)) GO TO 264
494 IF(ALAS
.LT
.HELIM
) GO TO 262
513 IF(INUB
.LE
.INU
) GO TO 225
526 CALL ZKSCL
(ZDR
,ZDI
,FNU
,N
,YR
,YI
,NZ
,RZR
,RZI
,ASCLE
,TOL
,ELIM
)
541 T2
= FNU
+ DBLE
(FLOAT
(KK
-1))
547 C-----------------------------------------------------------------------
548 C SCALE BY DEXP(Z), IFLAG = 1 CASES
549 C-----------------------------------------------------------------------
554 C-----------------------------------------------------------------------
555 C FNU=HALF ODD INTEGER CASE, DNU=-0.5
556 C-----------------------------------------------------------------------