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