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