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