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