1*DECK DBESK
2      SUBROUTINE DBESK (X, FNU, KODE, N, Y, NZ)
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      NN = -I1MACH(15)
90      ELIM = 2.303D0*(NN*D1MACH(5)-3.0D0)
91      XLIM = D1MACH(1)*1.0D+3
92      IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 280
93      IF (FNU.LT.0.0D0) GO TO 290
94      IF (X.LE.0.0D0) GO TO 300
95      IF (X.LT.XLIM) GO TO 320
96      IF (N.LT.1) GO TO 310
97      ETX = KODE - 1
98C
99C     ND IS A DUMMY VARIABLE FOR N
100C     GNU IS A DUMMY VARIABLE FOR FNU
101C     NZ = NUMBER OF UNDERFLOWS ON KODE=1
102C
103      ND = N
104      NZ = 0
105      NUD = INT(FNU)
106      DNU = FNU - NUD
107      GNU = FNU
108      NN = MIN(2,ND)
109      FN = FNU + N - 1
110      FNN = FN
111      IF (FN.LT.2.0D0) GO TO 150
112C
113C     OVERFLOW TEST  (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)
114C     FOR THE LAST ORDER, FNU+N-1.GE.NULIM
115C
116      ZN = X/FN
117      IF (ZN.EQ.0.0D0) GO TO 320
118      RTZ = SQRT(1.0D0+ZN*ZN)
119      GLN = LOG((1.0D0+RTZ)/ZN)
120      T = RTZ*(1.0D0-ETX) + ETX/(ZN+RTZ)
121      CN = -FN*(T-GLN)
122      IF (CN.GT.ELIM) GO TO 320
123      IF (NUD.LT.NULIM(NN)) GO TO 30
124      IF (NN.EQ.1) GO TO 20
125   10 CONTINUE
126C
127C     UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)
128C     FOR THE FIRST ORDER, FNU.GE.NULIM
129C
130      FN = GNU
131      ZN = X/FN
132      RTZ = SQRT(1.0D0+ZN*ZN)
133      GLN = LOG((1.0D0+RTZ)/ZN)
134      T = RTZ*(1.0D0-ETX) + ETX/(ZN+RTZ)
135      CN = -FN*(T-GLN)
136   20 CONTINUE
137      IF (CN.LT.-ELIM) GO TO 230
138C
139C     ASYMPTOTIC EXPANSION FOR ORDERS FNU AND FNU+1.GE.NULIM
140C
141      FLGIK = -1.0D0
142      CALL DASYIK(X,GNU,KODE,FLGIK,RTZ,CN,NN,Y)
143      IF (NN.EQ.1) GO TO 240
144      TRX = 2.0D0/X
145      TM = (GNU+GNU+2.0D0)/X
146      GO TO 130
147C
148   30 CONTINUE
149      IF (KODE.EQ.2) GO TO 40
150C
151C     UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION IN X)
152C     FOR ORDER DNU
153C
154      IF (X.GT.ELIM) GO TO 230
155   40 CONTINUE
156      IF (DNU.NE.0.0D0) GO TO 80
157      IF (KODE.EQ.2) GO TO 50
158      S1 = DBESK0(X)
159      GO TO 60
160   50 S1 = DBSK0E(X)
161   60 CONTINUE
162      IF (NUD.EQ.0 .AND. ND.EQ.1) GO TO 120
163      IF (KODE.EQ.2) GO TO 70
164      S2 = DBESK1(X)
165      GO TO 90
166   70 S2 = DBSK1E(X)
167      GO TO 90
168   80 CONTINUE
169      NB = 2
170      IF (NUD.EQ.0 .AND. ND.EQ.1) NB = 1
171      CALL DBSKNU(X, DNU, KODE, NB, W, NZ)
172      S1 = W(1)
173      IF (NB.EQ.1) GO TO 120
174      S2 = W(2)
175   90 CONTINUE
176      TRX = 2.0D0/X
177      TM = (DNU+DNU+2.0D0)/X
178C     FORWARD RECUR FROM DNU TO FNU+1 TO GET Y(1) AND Y(2)
179      IF (ND.EQ.1) NUD = NUD - 1
180      IF (NUD.GT.0) GO TO 100
181      IF (ND.GT.1) GO TO 120
182      S1 = S2
183      GO TO 120
184  100 CONTINUE
185      DO 110 I=1,NUD
186        S = S2
187        S2 = TM*S2 + S1
188        S1 = S
189        TM = TM + TRX
190  110 CONTINUE
191      IF (ND.EQ.1) S1 = S2
192  120 CONTINUE
193      Y(1) = S1
194      IF (ND.EQ.1) GO TO 240
195      Y(2) = S2
196  130 CONTINUE
197      IF (ND.EQ.2) GO TO 240
198C     FORWARD RECUR FROM FNU+2 TO FNU+N-1
199      DO 140 I=3,ND
200        Y(I) = TM*Y(I-1) + Y(I-2)
201        TM = TM + TRX
202  140 CONTINUE
203      GO TO 240
204C
205  150 CONTINUE
206C     UNDERFLOW TEST FOR KODE=1
207      IF (KODE.EQ.2) GO TO 160
208      IF (X.GT.ELIM) GO TO 230
209  160 CONTINUE
210C     OVERFLOW TEST
211      IF (FN.LE.1.0D0) GO TO 170
212      IF (-FN*(LOG(X)-0.693D0).GT.ELIM) GO TO 320
213  170 CONTINUE
214      IF (DNU.EQ.0.0D0) GO TO 180
215      CALL DBSKNU(X, FNU, KODE, ND, Y, MZ)
216      GO TO 240
217  180 CONTINUE
218      J = NUD
219      IF (J.EQ.1) GO TO 210
220      J = J + 1
221      IF (KODE.EQ.2) GO TO 190
222      Y(J) = DBESK0(X)
223      GO TO 200
224  190 Y(J) = DBSK0E(X)
225  200 IF (ND.EQ.1) GO TO 240
226      J = J + 1
227  210 IF (KODE.EQ.2) GO TO 220
228      Y(J) = DBESK1(X)
229      GO TO 240
230  220 Y(J) = DBSK1E(X)
231      GO TO 240
232C
233C     UPDATE PARAMETERS ON UNDERFLOW
234C
235  230 CONTINUE
236      NUD = NUD + 1
237      ND = ND - 1
238      IF (ND.EQ.0) GO TO 240
239      NN = MIN(2,ND)
240      GNU = GNU + 1.0D0
241      IF (FNN.LT.2.0D0) GO TO 230
242      IF (NUD.LT.NULIM(NN)) GO TO 230
243      GO TO 10
244  240 CONTINUE
245      NZ = N - ND
246      IF (NZ.EQ.0) RETURN
247      IF (ND.EQ.0) GO TO 260
248      DO 250 I=1,ND
249        J = N - I + 1
250        K = ND - I + 1
251        Y(J) = Y(K)
252  250 CONTINUE
253  260 CONTINUE
254      DO 270 I=1,NZ
255        Y(I) = 0.0D0
256  270 CONTINUE
257      RETURN
258C
259C
260C
261  280 CONTINUE
262      CALL XERMSG ('SLATEC', 'DBESK',
263     +   'SCALING OPTION, KODE, NOT 1 OR 2', 2, 1)
264      RETURN
265  290 CONTINUE
266      CALL XERMSG ('SLATEC', 'DBESK', 'ORDER, FNU, LESS THAN ZERO', 2,
267     +   1)
268      RETURN
269  300 CONTINUE
270      CALL XERMSG ('SLATEC', 'DBESK', 'X LESS THAN OR EQUAL TO ZERO',
271     +   2, 1)
272      RETURN
273  310 CONTINUE
274      CALL XERMSG ('SLATEC', 'DBESK', 'N LESS THAN ONE', 2, 1)
275      RETURN
276  320 CONTINUE
277      CALL XERMSG ('SLATEC', 'DBESK',
278     +   'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL', 6, 1)
279      RETURN
280      END
281