1*DECK DBESK 2 SUBROUTINE DBESK (X, FNU, KODE, N, Y, NZ,ierr) 3C***BEGIN PROLOGUE DBESK 4C***PURPOSE Implement forward recursion on the three term recursion 5C relation for a sequence of non-negative order Bessel 6C functions K/SUB(FNU+I-1)/(X), or scaled Bessel functions 7C EXP(X)*K/SUB(FNU+I-1)/(X), I=1,...,N for real, positive 8C X and non-negative orders FNU. 9C***LIBRARY SLATEC 10C***CATEGORY C10B3 11C***TYPE DOUBLE PRECISION (BESK-S, DBESK-D) 12C***KEYWORDS K BESSEL FUNCTION, SPECIAL FUNCTIONS 13C***AUTHOR Amos, D. E., (SNLA) 14C***DESCRIPTION 15C 16C Abstract **** a double precision routine **** 17C DBESK implements forward recursion on the three term 18C recursion relation for a sequence of non-negative order Bessel 19C functions K/sub(FNU+I-1)/(X), or scaled Bessel functions 20C EXP(X)*K/sub(FNU+I-1)/(X), I=1,..,N for real X .GT. 0.0D0 and 21C non-negative orders FNU. If FNU .LT. NULIM, orders FNU and 22C FNU+1 are obtained from DBSKNU to start the recursion. If 23C FNU .GE. NULIM, the uniform asymptotic expansion is used for 24C orders FNU and FNU+1 to start the recursion. NULIM is 35 or 25C 70 depending on whether N=1 or N .GE. 2. Under and overflow 26C tests are made on the leading term of the asymptotic expansion 27C before any extensive computation is done. 28C 29C The maximum number of significant digits obtainable 30C is the smaller of 14 and the number of digits carried in 31C double precision arithmetic. 32C 33C Description of Arguments 34C 35C Input X,FNU are double precision 36C X - X .GT. 0.0D0 37C FNU - order of the initial K function, FNU .GE. 0.0D0 38C KODE - a parameter to indicate the scaling option 39C KODE=1 returns Y(I)= K/sub(FNU+I-1)/(X), 40C I=1,...,N 41C KODE=2 returns Y(I)=EXP(X)*K/sub(FNU+I-1)/(X), 42C I=1,...,N 43C N - number of members in the sequence, N .GE. 1 44C 45C Output Y is double precision 46C Y - a vector whose first N components contain values 47C for the sequence 48C Y(I)= k/sub(FNU+I-1)/(X), I=1,...,N or 49C Y(I)=EXP(X)*K/sub(FNU+I-1)/(X), I=1,...,N 50C depending on KODE 51C NZ - number of components of Y set to zero due to 52C underflow with KODE=1, 53C NZ=0 , normal return, computation completed 54C NZ .NE. 0, first NZ components of Y set to zero 55C due to underflow, Y(I)=0.0D0, I=1,...,NZ 56C 57C Error Conditions 58C Improper input arguments - a fatal error 59C Overflow - a fatal error 60C Underflow with KODE=1 - a non-fatal error (NZ .NE. 0) 61C 62C***REFERENCES F. W. J. Olver, Tables of Bessel Functions of Moderate 63C or Large Orders, NPL Mathematical Tables 6, Her 64C Majesty's Stationery Office, London, 1962. 65C N. M. Temme, On the numerical evaluation of the modified 66C Bessel function of the third kind, Journal of 67C Computational Physics 19, (1975), pp. 324-337. 68C***ROUTINES CALLED D1MACH, DASYIK, DBESK0, DBESK1, DBSK0E, DBSK1E, 69C DBSKNU, I1MACH, XERMSG 70C***REVISION HISTORY (YYMMDD) 71C 790201 DATE WRITTEN 72C 890531 Changed all specific intrinsics to generic. (WRB) 73C 890911 Removed unnecessary intrinsics. (WRB) 74C 890911 REVISION DATE from Version 3.2 75C 891214 Prologue converted to Version 4.0 format. (BAB) 76C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 77C 920501 Reformatted the REFERENCES section. (WRB) 78C***END PROLOGUE DBESK 79C 80 INTEGER I, J, K, KODE, MZ, N, NB, ND, NN, NUD, NULIM, NZ 81 INTEGER I1MACH 82 DOUBLE PRECISION CN,DNU,ELIM,ETX,FLGIK,FN,FNN,FNU,GLN,GNU,RTZ, 83 1 S, S1, S2, T, TM, TRX, W, X, XLIM, Y, ZN 84 DOUBLE PRECISION DBESK0, DBESK1, DBSK1E, DBSK0E, D1MACH 85 DIMENSION W(2), NULIM(2), Y(*) 86 SAVE NULIM 87 DATA NULIM(1),NULIM(2) / 35 , 70 / 88C***FIRST EXECUTABLE STATEMENT DBESK 89 ierr=0 90 NN = -I1MACH(15) 91 ELIM = 2.303D0*(NN*D1MACH(5)-3.0D0) 92 XLIM = D1MACH(1)*1.0D+3 93 IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 280 94 IF (FNU.LT.0.0D0) GO TO 290 95 IF (X.LE.0.0D0) GO TO 300 96 IF (X.LT.XLIM) GO TO 320 97 IF (N.LT.1) GO TO 310 98 ETX = KODE - 1 99C 100C ND IS A DUMMY VARIABLE FOR N 101C GNU IS A DUMMY VARIABLE FOR FNU 102C NZ = NUMBER OF UNDERFLOWS ON KODE=1 103C 104 ND = N 105 NZ = 0 106 NUD = INT(FNU) 107 DNU = FNU - NUD 108 GNU = FNU 109 NN = MIN(2,ND) 110 FN = FNU + N - 1 111 FNN = FN 112 IF (FN.LT.2.0D0) GO TO 150 113C 114C OVERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION) 115C FOR THE LAST ORDER, FNU+N-1.GE.NULIM 116C 117 ZN = X/FN 118 IF (ZN.EQ.0.0D0) GO TO 320 119 RTZ = SQRT(1.0D0+ZN*ZN) 120 GLN = LOG((1.0D0+RTZ)/ZN) 121 T = RTZ*(1.0D0-ETX) + ETX/(ZN+RTZ) 122 CN = -FN*(T-GLN) 123 IF (CN.GT.ELIM) GO TO 320 124 IF (NUD.LT.NULIM(NN)) GO TO 30 125 IF (NN.EQ.1) GO TO 20 126 10 CONTINUE 127C 128C UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION) 129C FOR THE FIRST ORDER, FNU.GE.NULIM 130C 131 FN = GNU 132 ZN = X/FN 133 RTZ = SQRT(1.0D0+ZN*ZN) 134 GLN = LOG((1.0D0+RTZ)/ZN) 135 T = RTZ*(1.0D0-ETX) + ETX/(ZN+RTZ) 136 CN = -FN*(T-GLN) 137 20 CONTINUE 138 IF (CN.LT.-ELIM) GO TO 230 139C 140C ASYMPTOTIC EXPANSION FOR ORDERS FNU AND FNU+1.GE.NULIM 141C 142 FLGIK = -1.0D0 143 CALL DASYIK(X,GNU,KODE,FLGIK,RTZ,CN,NN,Y) 144 IF (NN.EQ.1) GO TO 240 145 TRX = 2.0D0/X 146 TM = (GNU+GNU+2.0D0)/X 147 GO TO 130 148C 149 30 CONTINUE 150 IF (KODE.EQ.2) GO TO 40 151C 152C UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION IN X) 153C FOR ORDER DNU 154C 155 IF (X.GT.ELIM) GO TO 230 156 40 CONTINUE 157 IF (DNU.NE.0.0D0) GO TO 80 158 IF (KODE.EQ.2) GO TO 50 159 S1 = DBESK0(X) 160 GO TO 60 161 50 S1 = DBSK0E(X) 162 60 CONTINUE 163 IF (NUD.EQ.0 .AND. ND.EQ.1) GO TO 120 164 IF (KODE.EQ.2) GO TO 70 165 S2 = DBESK1(X) 166 GO TO 90 167 70 S2 = DBSK1E(X) 168 GO TO 90 169 80 CONTINUE 170 NB = 2 171 IF (NUD.EQ.0 .AND. ND.EQ.1) NB = 1 172 CALL DBSKNU(X, DNU, KODE, NB, W, NZ) 173 S1 = W(1) 174 IF (NB.EQ.1) GO TO 120 175 S2 = W(2) 176 90 CONTINUE 177 TRX = 2.0D0/X 178 TM = (DNU+DNU+2.0D0)/X 179C FORWARD RECUR FROM DNU TO FNU+1 TO GET Y(1) AND Y(2) 180 IF (ND.EQ.1) NUD = NUD - 1 181 IF (NUD.GT.0) GO TO 100 182 IF (ND.GT.1) GO TO 120 183 S1 = S2 184 GO TO 120 185 100 CONTINUE 186 DO 110 I=1,NUD 187 S = S2 188 S2 = TM*S2 + S1 189 S1 = S 190 TM = TM + TRX 191 110 CONTINUE 192 IF (ND.EQ.1) S1 = S2 193 120 CONTINUE 194 Y(1) = S1 195 IF (ND.EQ.1) GO TO 240 196 Y(2) = S2 197 130 CONTINUE 198 IF (ND.EQ.2) GO TO 240 199C FORWARD RECUR FROM FNU+2 TO FNU+N-1 200 DO 140 I=3,ND 201 Y(I) = TM*Y(I-1) + Y(I-2) 202 TM = TM + TRX 203 140 CONTINUE 204 GO TO 240 205C 206 150 CONTINUE 207C UNDERFLOW TEST FOR KODE=1 208 IF (KODE.EQ.2) GO TO 160 209 IF (X.GT.ELIM) GO TO 230 210 160 CONTINUE 211C OVERFLOW TEST 212 IF (FN.LE.1.0D0) GO TO 170 213 IF (-FN*(LOG(X)-0.693D0).GT.ELIM) GO TO 320 214 170 CONTINUE 215 IF (DNU.EQ.0.0D0) GO TO 180 216 CALL DBSKNU(X, FNU, KODE, ND, Y, MZ) 217 GO TO 240 218 180 CONTINUE 219 J = NUD 220 IF (J.EQ.1) GO TO 210 221 J = J + 1 222 IF (KODE.EQ.2) GO TO 190 223 Y(J) = DBESK0(X) 224 GO TO 200 225 190 Y(J) = DBSK0E(X) 226 200 IF (ND.EQ.1) GO TO 240 227 J = J + 1 228 210 IF (KODE.EQ.2) GO TO 220 229 Y(J) = DBESK1(X) 230 GO TO 240 231 220 Y(J) = DBSK1E(X) 232 GO TO 240 233C 234C UPDATE PARAMETERS ON UNDERFLOW 235C 236 230 CONTINUE 237 NUD = NUD + 1 238 ND = ND - 1 239 IF (ND.EQ.0) GO TO 240 240 NN = MIN(2,ND) 241 GNU = GNU + 1.0D0 242 IF (FNN.LT.2.0D0) GO TO 230 243 IF (NUD.LT.NULIM(NN)) GO TO 230 244 GO TO 10 245 240 CONTINUE 246 NZ = N - ND 247 IF (NZ.EQ.0) RETURN 248 IF (ND.EQ.0) GO TO 260 249 DO 250 I=1,ND 250 J = N - I + 1 251 K = ND - I + 1 252 Y(J) = Y(K) 253 250 CONTINUE 254 260 CONTINUE 255 DO 270 I=1,NZ 256 Y(I) = 0.0D0 257 270 CONTINUE 258 RETURN 259C 260C 261C 262 280 CONTINUE 263C CALL XERMSG ('SLATEC', 'DBESK', 264C + 'SCALING OPTION, KODE, NOT 1 OR 2', 2, 1) 265 ierr=1 266 RETURN 267 290 CONTINUE 268C CALL XERMSG ('SLATEC', 'DBESK', 'ORDER, FNU, LESS THAN ZERO', 2, 269C + 1) 270 ierr=1 271 RETURN 272 300 CONTINUE 273C CALL XERMSG ('SLATEC', 'DBESK', 'X LESS THAN OR EQUAL TO ZERO', 274C + 2, 1) 275 ierr=1 276 RETURN 277 310 CONTINUE 278C CALL XERMSG ('SLATEC', 'DBESK', 'N LESS THAN ONE', 2, 1) 279 ierr=1 280 RETURN 281 320 CONTINUE 282C CALL XERMSG ('SLATEC', 'DBESK', 283C + 'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL', 6, 1) 284 ierr=2 285 RETURN 286 END 287