1*DECK DBSKNU 2 SUBROUTINE DBSKNU (X, FNU, KODE, N, Y, NZ) 3C***BEGIN PROLOGUE DBSKNU 4C***SUBSIDIARY 5C***PURPOSE Subsidiary to DBESK 6C***LIBRARY SLATEC 7C***TYPE DOUBLE PRECISION (BESKNU-S, DBSKNU-D) 8C***AUTHOR Amos, D. E., (SNLA) 9C***DESCRIPTION 10C 11C Abstract **** A DOUBLE PRECISION routine **** 12C DBSKNU computes N member sequences of K Bessel functions 13C K/SUB(FNU+I-1)/(X), I=1,N for non-negative orders FNU and 14C positive X. Equations of the references are implemented on 15C small orders DNU for K/SUB(DNU)/(X) and K/SUB(DNU+1)/(X). 16C Forward recursion with the three term recursion relation 17C generates higher orders FNU+I-1, I=1,...,N. The parameter 18C KODE permits K/SUB(FNU+I-1)/(X) values or scaled values 19C EXP(X)*K/SUB(FNU+I-1)/(X), I=1,N to be returned. 20C 21C To start the recursion FNU is normalized to the interval 22C -0.5.LE.DNU.LT.0.5. A special form of the power series is 23C implemented on 0.LT.X.LE.X1 while the Miller algorithm for the 24C K Bessel function in terms of the confluent hypergeometric 25C function U(FNU+0.5,2*FNU+1,X) is implemented on X1.LT.X.LE.X2. 26C For X.GT.X2, the asymptotic expansion for large X is used. 27C When FNU is a half odd integer, a special formula for 28C DNU=-0.5 and DNU+1.0=0.5 is used to start the recursion. 29C 30C The maximum number of significant digits obtainable 31C is the smaller of 14 and the number of digits carried in 32C DOUBLE PRECISION arithmetic. 33C 34C DBSKNU assumes that a significant digit SINH function is 35C available. 36C 37C Description of Arguments 38C 39C INPUT X,FNU are DOUBLE PRECISION 40C X - X.GT.0.0D0 41C FNU - Order of initial K function, FNU.GE.0.0D0 42C N - Number of members of the sequence, N.GE.1 43C KODE - A parameter to indicate the scaling option 44C KODE= 1 returns 45C Y(I)= K/SUB(FNU+I-1)/(X) 46C I=1,...,N 47C = 2 returns 48C Y(I)=EXP(X)*K/SUB(FNU+I-1)/(X) 49C I=1,...,N 50C 51C OUTPUT Y is DOUBLE PRECISION 52C Y - A vector whose first N components contain values 53C for the sequence 54C Y(I)= K/SUB(FNU+I-1)/(X), I=1,...,N or 55C Y(I)=EXP(X)*K/SUB(FNU+I-1)/(X), I=1,...,N 56C depending on KODE 57C NZ - Number of components set to zero due to 58C underflow, 59C NZ= 0 , normal return 60C NZ.NE.0 , first NZ components of Y set to zero 61C due to underflow, Y(I)=0.0D0,I=1,...,NZ 62C 63C Error Conditions 64C Improper input arguments - a fatal error 65C Overflow - a fatal error 66C Underflow with KODE=1 - a non-fatal error (NZ.NE.0) 67C 68C***SEE ALSO DBESK 69C***REFERENCES N. M. Temme, On the numerical evaluation of the modified 70C Bessel function of the third kind, Journal of 71C Computational Physics 19, (1975), pp. 324-337. 72C***ROUTINES CALLED D1MACH, DGAMMA, I1MACH, XERMSG 73C***REVISION HISTORY (YYMMDD) 74C 790201 DATE WRITTEN 75C 890531 Changed all specific intrinsics to generic. (WRB) 76C 890911 Removed unnecessary intrinsics. (WRB) 77C 891214 Prologue converted to Version 4.0 format. (BAB) 78C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 79C 900326 Removed duplicate information from DESCRIPTION section. 80C (WRB) 81C 900328 Added TYPE section. (WRB) 82C 900727 Added EXTERNAL statement. (WRB) 83C 910408 Updated the AUTHOR and REFERENCES sections. (WRB) 84C 920501 Reformatted the REFERENCES section. (WRB) 85C***END PROLOGUE DBSKNU 86C 87 INTEGER I, IFLAG, INU, J, K, KK, KODE, KODED, N, NN, NZ 88 INTEGER I1MACH 89 DOUBLE PRECISION A,AK,A1,A2,B,BK,CC,CK,COEF,CX,DK,DNU,DNU2,ELIM, 90 1 ETEST, EX, F, FC, FHS, FK, FKS, FLRX, FMU, FNU, G1, G2, P, PI, 91 2 PT, P1, P2, Q, RTHPI, RX, S, SMU, SQK, ST, S1, S2, TM, TOL, T1, 92 3 T2, X, X1, X2, Y 93 DIMENSION A(160), B(160), Y(*), CC(8) 94 DOUBLE PRECISION DGAMMA, D1MACH 95 EXTERNAL DGAMMA 96 SAVE X1, X2, PI, RTHPI, CC 97 DATA X1, X2 / 2.0D0, 17.0D0 / 98 DATA PI,RTHPI / 3.14159265358979D+00, 1.25331413731550D+00/ 99 DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8) 100 1 / 5.77215664901533D-01,-4.20026350340952D-02, 101 2-4.21977345555443D-02, 7.21894324666300D-03,-2.15241674114900D-04, 102 3-2.01348547807000D-05, 1.13302723200000D-06, 6.11609500000000D-09/ 103C***FIRST EXECUTABLE STATEMENT DBSKNU 104 KK = -I1MACH(15) 105 ELIM = 2.303D0*(KK*D1MACH(5)-3.0D0) 106 AK = D1MACH(3) 107 TOL = MAX(AK,1.0D-15) 108 IF (X.LE.0.0D0) GO TO 350 109 IF (FNU.LT.0.0D0) GO TO 360 110 IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 370 111 IF (N.LT.1) GO TO 380 112 NZ = 0 113 IFLAG = 0 114 KODED = KODE 115 RX = 2.0D0/X 116 INU = INT(FNU+0.5D0) 117 DNU = FNU - INU 118 IF (ABS(DNU).EQ.0.5D0) GO TO 120 119 DNU2 = 0.0D0 120 IF (ABS(DNU).LT.TOL) GO TO 10 121 DNU2 = DNU*DNU 122 10 CONTINUE 123 IF (X.GT.X1) GO TO 120 124C 125C SERIES FOR X.LE.X1 126C 127 A1 = 1.0D0 - DNU 128 A2 = 1.0D0 + DNU 129 T1 = 1.0D0/DGAMMA(A1) 130 T2 = 1.0D0/DGAMMA(A2) 131 IF (ABS(DNU).GT.0.1D0) GO TO 40 132C SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU) 133 S = CC(1) 134 AK = 1.0D0 135 DO 20 K=2,8 136 AK = AK*DNU2 137 TM = CC(K)*AK 138 S = S + TM 139 IF (ABS(TM).LT.TOL) GO TO 30 140 20 CONTINUE 141 30 G1 = -S 142 GO TO 50 143 40 CONTINUE 144 G1 = (T1-T2)/(DNU+DNU) 145 50 CONTINUE 146 G2 = (T1+T2)*0.5D0 147 SMU = 1.0D0 148 FC = 1.0D0 149 FLRX = LOG(RX) 150 FMU = DNU*FLRX 151 IF (DNU.EQ.0.0D0) GO TO 60 152 FC = DNU*PI 153 FC = FC/SIN(FC) 154 IF (FMU.NE.0.0D0) SMU = SINH(FMU)/FMU 155 60 CONTINUE 156 F = FC*(G1*COSH(FMU)+G2*FLRX*SMU) 157 FC = EXP(FMU) 158 P = 0.5D0*FC/T2 159 Q = 0.5D0/(FC*T1) 160 AK = 1.0D0 161 CK = 1.0D0 162 BK = 1.0D0 163 S1 = F 164 S2 = P 165 IF (INU.GT.0 .OR. N.GT.1) GO TO 90 166 IF (X.LT.TOL) GO TO 80 167 CX = X*X*0.25D0 168 70 CONTINUE 169 F = (AK*F+P+Q)/(BK-DNU2) 170 P = P/(AK-DNU) 171 Q = Q/(AK+DNU) 172 CK = CK*CX/AK 173 T1 = CK*F 174 S1 = S1 + T1 175 BK = BK + AK + AK + 1.0D0 176 AK = AK + 1.0D0 177 S = ABS(T1)/(1.0D0+ABS(S1)) 178 IF (S.GT.TOL) GO TO 70 179 80 CONTINUE 180 Y(1) = S1 181 IF (KODED.EQ.1) RETURN 182 Y(1) = S1*EXP(X) 183 RETURN 184 90 CONTINUE 185 IF (X.LT.TOL) GO TO 110 186 CX = X*X*0.25D0 187 100 CONTINUE 188 F = (AK*F+P+Q)/(BK-DNU2) 189 P = P/(AK-DNU) 190 Q = Q/(AK+DNU) 191 CK = CK*CX/AK 192 T1 = CK*F 193 S1 = S1 + T1 194 T2 = CK*(P-AK*F) 195 S2 = S2 + T2 196 BK = BK + AK + AK + 1.0D0 197 AK = AK + 1.0D0 198 S = ABS(T1)/(1.0D0+ABS(S1)) + ABS(T2)/(1.0D0+ABS(S2)) 199 IF (S.GT.TOL) GO TO 100 200 110 CONTINUE 201 S2 = S2*RX 202 IF (KODED.EQ.1) GO TO 170 203 F = EXP(X) 204 S1 = S1*F 205 S2 = S2*F 206 GO TO 170 207 120 CONTINUE 208 COEF = RTHPI/SQRT(X) 209 IF (KODED.EQ.2) GO TO 130 210 IF (X.GT.ELIM) GO TO 330 211 COEF = COEF*EXP(-X) 212 130 CONTINUE 213 IF (ABS(DNU).EQ.0.5D0) GO TO 340 214 IF (X.GT.X2) GO TO 280 215C 216C MILLER ALGORITHM FOR X1.LT.X.LE.X2 217C 218 ETEST = COS(PI*DNU)/(PI*X*TOL) 219 FKS = 1.0D0 220 FHS = 0.25D0 221 FK = 0.0D0 222 CK = X + X + 2.0D0 223 P1 = 0.0D0 224 P2 = 1.0D0 225 K = 0 226 140 CONTINUE 227 K = K + 1 228 FK = FK + 1.0D0 229 AK = (FHS-DNU2)/(FKS+FK) 230 BK = CK/(FK+1.0D0) 231 PT = P2 232 P2 = BK*P2 - AK*P1 233 P1 = PT 234 A(K) = AK 235 B(K) = BK 236 CK = CK + 2.0D0 237 FKS = FKS + FK + FK + 1.0D0 238 FHS = FHS + FK + FK 239 IF (ETEST.GT.FK*P1) GO TO 140 240 KK = K 241 S = 1.0D0 242 P1 = 0.0D0 243 P2 = 1.0D0 244 DO 150 I=1,K 245 PT = P2 246 P2 = (B(KK)*P2-P1)/A(KK) 247 P1 = PT 248 S = S + P2 249 KK = KK - 1 250 150 CONTINUE 251 S1 = COEF*(P2/S) 252 IF (INU.GT.0 .OR. N.GT.1) GO TO 160 253 GO TO 200 254 160 CONTINUE 255 S2 = S1*(X+DNU+0.5D0-P1/P2)/X 256C 257C FORWARD RECURSION ON THE THREE TERM RECURSION RELATION 258C 259 170 CONTINUE 260 CK = (DNU+DNU+2.0D0)/X 261 IF (N.EQ.1) INU = INU - 1 262 IF (INU.GT.0) GO TO 180 263 IF (N.GT.1) GO TO 200 264 S1 = S2 265 GO TO 200 266 180 CONTINUE 267 DO 190 I=1,INU 268 ST = S2 269 S2 = CK*S2 + S1 270 S1 = ST 271 CK = CK + RX 272 190 CONTINUE 273 IF (N.EQ.1) S1 = S2 274 200 CONTINUE 275 IF (IFLAG.EQ.1) GO TO 220 276 Y(1) = S1 277 IF (N.EQ.1) RETURN 278 Y(2) = S2 279 IF (N.EQ.2) RETURN 280 DO 210 I=3,N 281 Y(I) = CK*Y(I-1) + Y(I-2) 282 CK = CK + RX 283 210 CONTINUE 284 RETURN 285C IFLAG=1 CASES 286 220 CONTINUE 287 S = -X + LOG(S1) 288 Y(1) = 0.0D0 289 NZ = 1 290 IF (S.LT.-ELIM) GO TO 230 291 Y(1) = EXP(S) 292 NZ = 0 293 230 CONTINUE 294 IF (N.EQ.1) RETURN 295 S = -X + LOG(S2) 296 Y(2) = 0.0D0 297 NZ = NZ + 1 298 IF (S.LT.-ELIM) GO TO 240 299 NZ = NZ - 1 300 Y(2) = EXP(S) 301 240 CONTINUE 302 IF (N.EQ.2) RETURN 303 KK = 2 304 IF (NZ.LT.2) GO TO 260 305 DO 250 I=3,N 306 KK = I 307 ST = S2 308 S2 = CK*S2 + S1 309 S1 = ST 310 CK = CK + RX 311 S = -X + LOG(S2) 312 NZ = NZ + 1 313 Y(I) = 0.0D0 314 IF (S.LT.-ELIM) GO TO 250 315 Y(I) = EXP(S) 316 NZ = NZ - 1 317 GO TO 260 318 250 CONTINUE 319 RETURN 320 260 CONTINUE 321 IF (KK.EQ.N) RETURN 322 S2 = S2*CK + S1 323 CK = CK + RX 324 KK = KK + 1 325 Y(KK) = EXP(-X+LOG(S2)) 326 IF (KK.EQ.N) RETURN 327 KK = KK + 1 328 DO 270 I=KK,N 329 Y(I) = CK*Y(I-1) + Y(I-2) 330 CK = CK + RX 331 270 CONTINUE 332 RETURN 333C 334C ASYMPTOTIC EXPANSION FOR LARGE X, X.GT.X2 335C 336C IFLAG=0 MEANS NO UNDERFLOW OCCURRED 337C IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH 338C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD 339C RECURSION 340 280 CONTINUE 341 NN = 2 342 IF (INU.EQ.0 .AND. N.EQ.1) NN = 1 343 DNU2 = DNU + DNU 344 FMU = 0.0D0 345 IF (ABS(DNU2).LT.TOL) GO TO 290 346 FMU = DNU2*DNU2 347 290 CONTINUE 348 EX = X*8.0D0 349 S2 = 0.0D0 350 DO 320 K=1,NN 351 S1 = S2 352 S = 1.0D0 353 AK = 0.0D0 354 CK = 1.0D0 355 SQK = 1.0D0 356 DK = EX 357 DO 300 J=1,30 358 CK = CK*(FMU-SQK)/DK 359 S = S + CK 360 DK = DK + EX 361 AK = AK + 8.0D0 362 SQK = SQK + AK 363 IF (ABS(CK).LT.TOL) GO TO 310 364 300 CONTINUE 365 310 S2 = S*COEF 366 FMU = FMU + 8.0D0*DNU + 4.0D0 367 320 CONTINUE 368 IF (NN.GT.1) GO TO 170 369 S1 = S2 370 GO TO 200 371 330 CONTINUE 372 KODED = 2 373 IFLAG = 1 374 GO TO 120 375C 376C FNU=HALF ODD INTEGER CASE 377C 378 340 CONTINUE 379 S1 = COEF 380 S2 = COEF 381 GO TO 170 382C 383C 384 350 CALL XERMSG ('SLATEC', 'DBSKNU', 'X NOT GREATER THAN ZERO', 2, 1) 385 RETURN 386 360 CALL XERMSG ('SLATEC', 'DBSKNU', 'FNU NOT ZERO OR POSITIVE', 2, 387 + 1) 388 RETURN 389 370 CALL XERMSG ('SLATEC', 'DBSKNU', 'KODE NOT 1 OR 2', 2, 1) 390 RETURN 391 380 CALL XERMSG ('SLATEC', 'DBSKNU', 'N NOT GREATER THAN 0', 2, 1) 392 RETURN 393 END 394