1*DECK DBESY 2 SUBROUTINE DBESY (X, FNU, N, Y) 3C***BEGIN PROLOGUE DBESY 4C***PURPOSE Implement forward recursion on the three term recursion 5C relation for a sequence of non-negative order Bessel 6C functions Y/SUB(FNU+I-1)/(X), I=1,...,N for real, positive 7C X and non-negative orders FNU. 8C***LIBRARY SLATEC 9C***CATEGORY C10A3 10C***TYPE DOUBLE PRECISION (BESY-S, DBESY-D) 11C***KEYWORDS SPECIAL FUNCTIONS, Y BESSEL FUNCTION 12C***AUTHOR Amos, D. E., (SNLA) 13C***DESCRIPTION 14C 15C Abstract **** a double precision routine **** 16C DBESY implements forward recursion on the three term 17C recursion relation for a sequence of non-negative order Bessel 18C functions Y/sub(FNU+I-1)/(X), I=1,N for real X .GT. 0.0D0 and 19C non-negative orders FNU. If FNU .LT. NULIM, orders FNU and 20C FNU+1 are obtained from DBSYNU which computes by a power 21C series for X .LE. 2, the K Bessel function of an imaginary 22C argument for 2 .LT. X .LE. 20 and the asymptotic expansion for 23C X .GT. 20. 24C 25C If FNU .GE. NULIM, the uniform asymptotic expansion is coded 26C in DASYJY for orders FNU and FNU+1 to start the recursion. 27C NULIM is 70 or 100 depending on whether N=1 or N .GE. 2. An 28C overflow test is made on the leading term of the asymptotic 29C expansion before any extensive computation is done. 30C 31C The maximum number of significant digits obtainable 32C is the smaller of 14 and the number of digits carried in 33C double precision arithmetic. 34C 35C Description of Arguments 36C 37C Input 38C X - X .GT. 0.0D0 39C FNU - order of the initial Y function, FNU .GE. 0.0D0 40C N - number of members in the sequence, N .GE. 1 41C 42C Output 43C Y - a vector whose first N components contain values 44C for the sequence Y(I)=Y/sub(FNU+I-1)/(X), I=1,N. 45C 46C Error Conditions 47C Improper input arguments - a fatal error 48C Overflow - a fatal error 49C 50C***REFERENCES F. W. J. Olver, Tables of Bessel Functions of Moderate 51C or Large Orders, NPL Mathematical Tables 6, Her 52C Majesty's Stationery Office, London, 1962. 53C N. M. Temme, On the numerical evaluation of the modified 54C Bessel function of the third kind, Journal of 55C Computational Physics 19, (1975), pp. 324-337. 56C N. M. Temme, On the numerical evaluation of the ordinary 57C Bessel function of the second kind, Journal of 58C Computational Physics 21, (1976), pp. 343-350. 59C***ROUTINES CALLED D1MACH, DASYJY, DBESY0, DBESY1, DBSYNU, DYAIRY, 60C I1MACH, XERMSG 61C***REVISION HISTORY (YYMMDD) 62C 800501 DATE WRITTEN 63C 890531 Changed all specific intrinsics to generic. (WRB) 64C 890911 Removed unnecessary intrinsics. (WRB) 65C 890911 REVISION DATE from Version 3.2 66C 891214 Prologue converted to Version 4.0 format. (BAB) 67C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 68C 920501 Reformatted the REFERENCES section. (WRB) 69C***END PROLOGUE DBESY 70C 71 EXTERNAL DYAIRY 72 INTEGER I, IFLW, J, N, NB, ND, NN, NUD, NULIM 73 INTEGER I1MACH 74 DOUBLE PRECISION AZN,CN,DNU,ELIM,FLGJY,FN,FNU,RAN,S,S1,S2,TM,TRX, 75 1 W,WK,W2N,X,XLIM,XXN,Y 76 DOUBLE PRECISION DBESY0, DBESY1, D1MACH 77 DIMENSION W(2), NULIM(2), Y(*), WK(7) 78 SAVE NULIM 79 DATA NULIM(1),NULIM(2) / 70 , 100 / 80C***FIRST EXECUTABLE STATEMENT DBESY 81 NN = -I1MACH(15) 82 ELIM = 2.303D0*(NN*D1MACH(5)-3.0D0) 83 XLIM = D1MACH(1)*1.0D+3 84 IF (FNU.LT.0.0D0) GO TO 140 85 IF (X.LE.0.0D0) GO TO 150 86 IF (X.LT.XLIM) GO TO 170 87 IF (N.LT.1) GO TO 160 88C 89C ND IS A DUMMY VARIABLE FOR N 90C 91 ND = N 92 NUD = INT(FNU) 93 DNU = FNU - NUD 94 NN = MIN(2,ND) 95 FN = FNU + N - 1 96 IF (FN.LT.2.0D0) GO TO 100 97C 98C OVERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION) 99C FOR THE LAST ORDER, FNU+N-1.GE.NULIM 100C 101 XXN = X/FN 102 W2N = 1.0D0-XXN*XXN 103 IF(W2N.LE.0.0D0) GO TO 10 104 RAN = SQRT(W2N) 105 AZN = LOG((1.0D0+RAN)/XXN) - RAN 106 CN = FN*AZN 107 IF(CN.GT.ELIM) GO TO 170 108 10 CONTINUE 109 IF (NUD.LT.NULIM(NN)) GO TO 20 110C 111C ASYMPTOTIC EXPANSION FOR ORDERS FNU AND FNU+1.GE.NULIM 112C 113 FLGJY = -1.0D0 114 CALL DASYJY(DYAIRY,X,FNU,FLGJY,NN,Y,WK,IFLW) 115 IF(IFLW.NE.0) GO TO 170 116 IF (NN.EQ.1) RETURN 117 TRX = 2.0D0/X 118 TM = (FNU+FNU+2.0D0)/X 119 GO TO 80 120C 121 20 CONTINUE 122 IF (DNU.NE.0.0D0) GO TO 30 123 S1 = DBESY0(X) 124 IF (NUD.EQ.0 .AND. ND.EQ.1) GO TO 70 125 S2 = DBESY1(X) 126 GO TO 40 127 30 CONTINUE 128 NB = 2 129 IF (NUD.EQ.0 .AND. ND.EQ.1) NB = 1 130 CALL DBSYNU(X, DNU, NB, W) 131 S1 = W(1) 132 IF (NB.EQ.1) GO TO 70 133 S2 = W(2) 134 40 CONTINUE 135 TRX = 2.0D0/X 136 TM = (DNU+DNU+2.0D0)/X 137C FORWARD RECUR FROM DNU TO FNU+1 TO GET Y(1) AND Y(2) 138 IF (ND.EQ.1) NUD = NUD - 1 139 IF (NUD.GT.0) GO TO 50 140 IF (ND.GT.1) GO TO 70 141 S1 = S2 142 GO TO 70 143 50 CONTINUE 144 DO 60 I=1,NUD 145 S = S2 146 S2 = TM*S2 - S1 147 S1 = S 148 TM = TM + TRX 149 60 CONTINUE 150 IF (ND.EQ.1) S1 = S2 151 70 CONTINUE 152 Y(1) = S1 153 IF (ND.EQ.1) RETURN 154 Y(2) = S2 155 80 CONTINUE 156 IF (ND.EQ.2) RETURN 157C FORWARD RECUR FROM FNU+2 TO FNU+N-1 158 DO 90 I=3,ND 159 Y(I) = TM*Y(I-1) - Y(I-2) 160 TM = TM + TRX 161 90 CONTINUE 162 RETURN 163C 164 100 CONTINUE 165C OVERFLOW TEST 166 IF (FN.LE.1.0D0) GO TO 110 167 IF (-FN*(LOG(X)-0.693D0).GT.ELIM) GO TO 170 168 110 CONTINUE 169 IF (DNU.EQ.0.0D0) GO TO 120 170 CALL DBSYNU(X, FNU, ND, Y) 171 RETURN 172 120 CONTINUE 173 J = NUD 174 IF (J.EQ.1) GO TO 130 175 J = J + 1 176 Y(J) = DBESY0(X) 177 IF (ND.EQ.1) RETURN 178 J = J + 1 179 130 CONTINUE 180 Y(J) = DBESY1(X) 181 IF (ND.EQ.1) RETURN 182 TRX = 2.0D0/X 183 TM = TRX 184 GO TO 80 185C 186C 187C 188 140 CONTINUE 189 CALL XERMSG ('SLATEC', 'DBESY', 'ORDER, FNU, LESS THAN ZERO', 2, 190 + 1) 191 RETURN 192 150 CONTINUE 193 CALL XERMSG ('SLATEC', 'DBESY', 'X LESS THAN OR EQUAL TO ZERO', 194 + 2, 1) 195 RETURN 196 160 CONTINUE 197 CALL XERMSG ('SLATEC', 'DBESY', 'N LESS THAN ONE', 2, 1) 198 RETURN 199 170 CONTINUE 200 CALL XERMSG ('SLATEC', 'DBESY', 201 + 'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL', 6, 1) 202 RETURN 203 END 204