1*DECK DBSYNU
2      SUBROUTINE DBSYNU (X, FNU, N, Y)
3C***BEGIN PROLOGUE  DBSYNU
4C***SUBSIDIARY
5C***PURPOSE  Subsidiary to DBESY
6C***LIBRARY   SLATEC
7C***TYPE      DOUBLE PRECISION (BESYNU-S, DBSYNU-D)
8C***AUTHOR  Amos, D. E., (SNLA)
9C***DESCRIPTION
10C
11C     Abstract  **** A DOUBLE PRECISION routine ****
12C         DBSYNU computes N member sequences of Y Bessel functions
13C         Y/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 Y/SUB(DNU)/(X) and Y/SUB(DNU+1)/(X).
16C         Forward recursion with the three term recursion relation
17C         generates higher orders FNU+I-1, I=1,...,N.
18C
19C         To start the recursion FNU is normalized to the interval
20C         -0.5.LE.DNU.LT.0.5. A special form of the power series is
21C         implemented on 0.LT.X.LE.X1 while the Miller algorithm for the
22C         K Bessel function in terms of the confluent hypergeometric
23C         function U(FNU+0.5,2*FNU+1,I*X) is implemented on X1.LT.X.LE.X
24C         Here I is the complex number SQRT(-1.).
25C         For X.GT.X2, the asymptotic expansion for large X is used.
26C         When FNU is a half odd integer, a special formula for
27C         DNU=-0.5 and DNU+1.0=0.5 is used to start the recursion.
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         DBSYNU assumes that a significant digit SINH function is
34C         available.
35C
36C     Description of Arguments
37C
38C         INPUT
39C           X      - X.GT.0.0D0
40C           FNU    - Order of initial Y function, FNU.GE.0.0D0
41C           N      - Number of members of the sequence, N.GE.1
42C
43C         OUTPUT
44C           Y      - A vector whose first N components contain values
45C                    for the sequence Y(I)=Y/SUB(FNU+I-1), I=1,N.
46C
47C     Error Conditions
48C         Improper input arguments - a fatal error
49C         Overflow - a fatal error
50C
51C***SEE ALSO  DBESY
52C***REFERENCES  N. M. Temme, On the numerical evaluation of the ordinary
53C                 Bessel function of the second kind, Journal of
54C                 Computational Physics 21, (1976), pp. 343-350.
55C               N. M. Temme, On the numerical evaluation of the modified
56C                 Bessel function of the third kind, Journal of
57C                 Computational Physics 19, (1975), pp. 324-337.
58C***ROUTINES CALLED  D1MACH, DGAMMA, XERMSG
59C***REVISION HISTORY  (YYMMDD)
60C   800501  DATE WRITTEN
61C   890531  Changed all specific intrinsics to generic.  (WRB)
62C   890911  Removed unnecessary intrinsics.  (WRB)
63C   891214  Prologue converted to Version 4.0 format.  (BAB)
64C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
65C   900326  Removed duplicate information from DESCRIPTION section.
66C           (WRB)
67C   900328  Added TYPE section.  (WRB)
68C   900727  Added EXTERNAL statement.  (WRB)
69C   910408  Updated the AUTHOR and REFERENCES sections.  (WRB)
70C   920501  Reformatted the REFERENCES section.  (WRB)
71C***END PROLOGUE  DBSYNU
72C
73      INTEGER I, INU, J, K, KK, N, NN
74      DOUBLE PRECISION A,AK,ARG,A1,A2,BK,CB,CBK,CC,CCK,CK,COEF,CPT,
75     1 CP1, CP2, CS, CS1, CS2, CX, DNU, DNU2, ETEST, ETX, F, FC, FHS,
76     2 FK, FKS, FLRX, FMU, FN, FNU, FX, G, G1, G2, HPI, P, PI, PT, Q,
77     3 RB, RBK, RCK, RELB, RPT, RP1, RP2, RS, RS1, RS2, RTHPI, RX, S,
78     4 SA, SB, SMU, SS, ST, S1, S2, TB, TM, TOL, T1, T2, X, X1, X2, Y
79      DIMENSION A(120), RB(120), CB(120), Y(*), CC(8)
80      DOUBLE PRECISION DGAMMA, D1MACH
81      EXTERNAL DGAMMA
82      SAVE X1, X2,PI, RTHPI, HPI, CC
83      DATA X1, X2 / 3.0D0, 20.0D0 /
84      DATA PI,RTHPI        / 3.14159265358979D+00, 7.97884560802865D-01/
85      DATA HPI             / 1.57079632679490D+00/
86      DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)
87     1                     / 5.77215664901533D-01,-4.20026350340952D-02,
88     2-4.21977345555443D-02, 7.21894324666300D-03,-2.15241674114900D-04,
89     3-2.01348547807000D-05, 1.13302723200000D-06, 6.11609500000000D-09/
90C***FIRST EXECUTABLE STATEMENT  DBSYNU
91      AK = D1MACH(3)
92      TOL = MAX(AK,1.0D-15)
93      IF (X.LE.0.0D0) GO TO 270
94      IF (FNU.LT.0.0D0) GO TO 280
95      IF (N.LT.1) GO TO 290
96      RX = 2.0D0/X
97      INU = INT(FNU+0.5D0)
98      DNU = FNU - INU
99      IF (ABS(DNU).EQ.0.5D0) GO TO 260
100      DNU2 = 0.0D0
101      IF (ABS(DNU).LT.TOL) GO TO 10
102      DNU2 = DNU*DNU
103   10 CONTINUE
104      IF (X.GT.X1) GO TO 120
105C
106C     SERIES FOR X.LE.X1
107C
108      A1 = 1.0D0 - DNU
109      A2 = 1.0D0 + DNU
110      T1 = 1.0D0/DGAMMA(A1)
111      T2 = 1.0D0/DGAMMA(A2)
112      IF (ABS(DNU).GT.0.1D0) GO TO 40
113C     SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
114      S = CC(1)
115      AK = 1.0D0
116      DO 20 K=2,8
117        AK = AK*DNU2
118        TM = CC(K)*AK
119        S = S + TM
120        IF (ABS(TM).LT.TOL) GO TO 30
121   20 CONTINUE
122   30 G1 = -(S+S)
123      GO TO 50
124   40 CONTINUE
125      G1 = (T1-T2)/DNU
126   50 CONTINUE
127      G2 = T1 + T2
128      SMU = 1.0D0
129      FC = 1.0D0/PI
130      FLRX = LOG(RX)
131      FMU = DNU*FLRX
132      TM = 0.0D0
133      IF (DNU.EQ.0.0D0) GO TO 60
134      TM = SIN(DNU*HPI)/DNU
135      TM = (DNU+DNU)*TM*TM
136      FC = DNU/SIN(DNU*PI)
137      IF (FMU.NE.0.0D0) SMU = SINH(FMU)/FMU
138   60 CONTINUE
139      F = FC*(G1*COSH(FMU)+G2*FLRX*SMU)
140      FX = EXP(FMU)
141      P = FC*T1*FX
142      Q = FC*T2/FX
143      G = F + TM*Q
144      AK = 1.0D0
145      CK = 1.0D0
146      BK = 1.0D0
147      S1 = G
148      S2 = P
149      IF (INU.GT.0 .OR. N.GT.1) GO TO 90
150      IF (X.LT.TOL) GO TO 80
151      CX = X*X*0.25D0
152   70 CONTINUE
153      F = (AK*F+P+Q)/(BK-DNU2)
154      P = P/(AK-DNU)
155      Q = Q/(AK+DNU)
156      G = F + TM*Q
157      CK = -CK*CX/AK
158      T1 = CK*G
159      S1 = S1 + T1
160      BK = BK + AK + AK + 1.0D0
161      AK = AK + 1.0D0
162      S = ABS(T1)/(1.0D0+ABS(S1))
163      IF (S.GT.TOL) GO TO 70
164   80 CONTINUE
165      Y(1) = -S1
166      RETURN
167   90 CONTINUE
168      IF (X.LT.TOL) GO TO 110
169      CX = X*X*0.25D0
170  100 CONTINUE
171      F = (AK*F+P+Q)/(BK-DNU2)
172      P = P/(AK-DNU)
173      Q = Q/(AK+DNU)
174      G = F + TM*Q
175      CK = -CK*CX/AK
176      T1 = CK*G
177      S1 = S1 + T1
178      T2 = CK*(P-AK*G)
179      S2 = S2 + T2
180      BK = BK + AK + AK + 1.0D0
181      AK = AK + 1.0D0
182      S = ABS(T1)/(1.0D0+ABS(S1)) + ABS(T2)/(1.0D0+ABS(S2))
183      IF (S.GT.TOL) GO TO 100
184  110 CONTINUE
185      S2 = -S2*RX
186      S1 = -S1
187      GO TO 160
188  120 CONTINUE
189      COEF = RTHPI/SQRT(X)
190      IF (X.GT.X2) GO TO 210
191C
192C     MILLER ALGORITHM FOR X1.LT.X.LE.X2
193C
194      ETEST = COS(PI*DNU)/(PI*X*TOL)
195      FKS = 1.0D0
196      FHS = 0.25D0
197      FK = 0.0D0
198      RCK = 2.0D0
199      CCK = X + X
200      RP1 = 0.0D0
201      CP1 = 0.0D0
202      RP2 = 1.0D0
203      CP2 = 0.0D0
204      K = 0
205  130 CONTINUE
206      K = K + 1
207      FK = FK + 1.0D0
208      AK = (FHS-DNU2)/(FKS+FK)
209      PT = FK + 1.0D0
210      RBK = RCK/PT
211      CBK = CCK/PT
212      RPT = RP2
213      CPT = CP2
214      RP2 = RBK*RPT - CBK*CPT - AK*RP1
215      CP2 = CBK*RPT + RBK*CPT - AK*CP1
216      RP1 = RPT
217      CP1 = CPT
218      RB(K) = RBK
219      CB(K) = CBK
220      A(K) = AK
221      RCK = RCK + 2.0D0
222      FKS = FKS + FK + FK + 1.0D0
223      FHS = FHS + FK + FK
224      PT = MAX(ABS(RP1),ABS(CP1))
225      FC = (RP1/PT)**2 + (CP1/PT)**2
226      PT = PT*SQRT(FC)*FK
227      IF (ETEST.GT.PT) GO TO 130
228      KK = K
229      RS = 1.0D0
230      CS = 0.0D0
231      RP1 = 0.0D0
232      CP1 = 0.0D0
233      RP2 = 1.0D0
234      CP2 = 0.0D0
235      DO 140 I=1,K
236        RPT = RP2
237        CPT = CP2
238        RP2 = (RB(KK)*RPT-CB(KK)*CPT-RP1)/A(KK)
239        CP2 = (CB(KK)*RPT+RB(KK)*CPT-CP1)/A(KK)
240        RP1 = RPT
241        CP1 = CPT
242        RS = RS + RP2
243        CS = CS + CP2
244        KK = KK - 1
245  140 CONTINUE
246      PT = MAX(ABS(RS),ABS(CS))
247      FC = (RS/PT)**2 + (CS/PT)**2
248      PT = PT*SQRT(FC)
249      RS1 = (RP2*(RS/PT)+CP2*(CS/PT))/PT
250      CS1 = (CP2*(RS/PT)-RP2*(CS/PT))/PT
251      FC = HPI*(DNU-0.5D0) - X
252      P = COS(FC)
253      Q = SIN(FC)
254      S1 = (CS1*Q-RS1*P)*COEF
255      IF (INU.GT.0 .OR. N.GT.1) GO TO 150
256      Y(1) = S1
257      RETURN
258  150 CONTINUE
259      PT = MAX(ABS(RP2),ABS(CP2))
260      FC = (RP2/PT)**2 + (CP2/PT)**2
261      PT = PT*SQRT(FC)
262      RPT = DNU + 0.5D0 - (RP1*(RP2/PT)+CP1*(CP2/PT))/PT
263      CPT = X - (CP1*(RP2/PT)-RP1*(CP2/PT))/PT
264      CS2 = CS1*CPT - RS1*RPT
265      RS2 = RPT*CS1 + RS1*CPT
266      S2 = (RS2*Q+CS2*P)*COEF/X
267C
268C     FORWARD RECURSION ON THE THREE TERM RECURSION RELATION
269C
270  160 CONTINUE
271      CK = (DNU+DNU+2.0D0)/X
272      IF (N.EQ.1) INU = INU - 1
273      IF (INU.GT.0) GO TO 170
274      IF (N.GT.1) GO TO 190
275      S1 = S2
276      GO TO 190
277  170 CONTINUE
278      DO 180 I=1,INU
279        ST = S2
280        S2 = CK*S2 - S1
281        S1 = ST
282        CK = CK + RX
283  180 CONTINUE
284      IF (N.EQ.1) S1 = S2
285  190 CONTINUE
286      Y(1) = S1
287      IF (N.EQ.1) RETURN
288      Y(2) = S2
289      IF (N.EQ.2) RETURN
290      DO 200 I=3,N
291        Y(I) = CK*Y(I-1) - Y(I-2)
292        CK = CK + RX
293  200 CONTINUE
294      RETURN
295C
296C     ASYMPTOTIC EXPANSION FOR LARGE X, X.GT.X2
297C
298  210 CONTINUE
299      NN = 2
300      IF (INU.EQ.0 .AND. N.EQ.1) NN = 1
301      DNU2 = DNU + DNU
302      FMU = 0.0D0
303      IF (ABS(DNU2).LT.TOL) GO TO 220
304      FMU = DNU2*DNU2
305  220 CONTINUE
306      ARG = X - HPI*(DNU+0.5D0)
307      SA = SIN(ARG)
308      SB = COS(ARG)
309      ETX = 8.0D0*X
310      DO 250 K=1,NN
311        S1 = S2
312        T2 = (FMU-1.0D0)/ETX
313        SS = T2
314        RELB = TOL*ABS(T2)
315        T1 = ETX
316        S = 1.0D0
317        FN = 1.0D0
318        AK = 0.0D0
319        DO 230 J=1,13
320          T1 = T1 + ETX
321          AK = AK + 8.0D0
322          FN = FN + AK
323          T2 = -T2*(FMU-FN)/T1
324          S = S + T2
325          T1 = T1 + ETX
326          AK = AK + 8.0D0
327          FN = FN + AK
328          T2 = T2*(FMU-FN)/T1
329          SS = SS + T2
330          IF (ABS(T2).LE.RELB) GO TO 240
331  230   CONTINUE
332  240   S2 = COEF*(S*SA+SS*SB)
333        FMU = FMU + 8.0D0*DNU + 4.0D0
334        TB = SA
335        SA = -SB
336        SB = TB
337  250 CONTINUE
338      IF (NN.GT.1) GO TO 160
339      S1 = S2
340      GO TO 190
341C
342C     FNU=HALF ODD INTEGER CASE
343C
344  260 CONTINUE
345      COEF = RTHPI/SQRT(X)
346      S1 = COEF*SIN(X)
347      S2 = -COEF*COS(X)
348      GO TO 160
349C
350C
351  270 CALL XERMSG ('SLATEC', 'DBSYNU', 'X NOT GREATER THAN ZERO', 2, 1)
352      RETURN
353  280 CALL XERMSG ('SLATEC', 'DBSYNU', 'FNU NOT ZERO OR POSITIVE', 2,
354     +   1)
355      RETURN
356  290 CALL XERMSG ('SLATEC', 'DBSYNU', 'N NOT GREATER THAN 0', 2, 1)
357      RETURN
358      END
359