1 subroutine DBESYN(X, ALPHA, NUM, BY) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5C>> 2009-10-28 DBESYN Krogh/Snyder Moved BGAMSQ = BGAM**2 inside else. 6C>> 2000-01-05 DBESYN Krogh Changed c$ comments to c. (Cray f90 prob.) 7C>> 1995-11-13 DBESYN Krogh Converted from SFTRAN to Fortran 8C>> 1994-10-19 DBESYN Krogh Changes to use M77CON 9C>> 1994-04-19 DBESYN CLL Edited to make DP & SP files similar. 10C>> 1991-01-14 DBESYN CLL Removed duplicate data statement. 11C>> 1989-08-09 DBESYN CLL More accurate constants for Cray 12C>> 1986-03-18 DBESYN Lawson Initial code. 13c--D replaces "?": ?BESYN, ?BESPQ, ?ERV1, ?GAMMA, ?LGAMA 14C 15c This subr computes the Y Bessel functions of X for 16c NUM orders from ALPHA through ALPHA + NUM - 1. 17c The results will be stored in BY(I), I = 1,...,NUM. 18C 19c Require X .gt. 0., ALPHA .ge. 0., NUM .ge. 1. 20C 21c The original subroutines SBYNU and BESJ/BESY were 22c designed and programmed by E. W. Ng and W. V. Snyder, JPL, 23C in 1973. Modified by Ng and S. Singletary, JPL, 1974. 24c 25c 1984 April 16, JPL, C. Lawson and S. Chan. Original subrs 26c combined into single subroutine. Converted to SFTRAN3 27c and Fortran 77 for the JPL MATH77 library. 28C 29c ------------------------------------------------------------------ 30c 31c As the machine precision, EPS, increases so does XPQ, 32c and thus so does the requirement for the storage dimension, 33c LDIMA. Here are some values of -LOG10(EPS), XPQ, and LDIMA. 34c 35c -LOG10(EPS)= 5 10 15 20 25 30 36c XPQ = 5.74 11.38 17.03 22.68 28.32 33.97 37c LDIMA = 17 33 49 63 79 95 38c 39c Since LDIMA cannot be adjusted at run-time, we are 40c setting it to 95 to handle the case of CDC double precision 41c which is probably the largest precision system likely 42c to be encountered. Note that the relative precision of results 43c from this subr cannot exceed about 16 or 17 decimal places 44c because of the limited accuracy of the polynomial coeffs used to 45c compute G0, and also the limited precision of the 46c subprograms referenced for gamma and loggamma. 47c ------------------------------------------------------------------ 48c Subprograms used: TANH, LOG10, LOG, EXP, COS, SIN, SQRT 49 external D1MACH, DBESPQ, DERV1, ERMSG, ERMOR, IERV1 50 external DGAMMA, DLGAMA 51c ---------- 52 integer I, II, K, LDIMA, M, MU, MUSAVE, ND2, NMAX, NMIN, NUM 53 parameter( LDIMA = 95, ND2 = 10) 54 double precision D1MACH, DGAMMA, DLGAMA 55 double precision AJ(LDIMA), ALPHA, ARG 56 double precision BGAM, BGAMSQ, BIG, BIGLOG, BY(NUM) 57 double precision C11293, C16, C1P5, C2BYPI, C59, CHI, CP6, CUTLOW 58 double precision D1, D2, D2SER(0:ND2), D2VAL, D3, DR 59 double precision EC, EM, EM1, EMU, EN1, EN2, EN3, EPS, ETA 60 double precision FAC, FK, FKM1, FKP1, FOUR, FV, FVM1, FVP1 61 double precision G, G0, G1, GMAX, GNU, HALF, HALFPI, HICUT 62 double precision LOGPI, LOGTWO, ONE, P, PI, PIV2, PSI1, PSIZ 63 double precision Q, Q2DXPV, SCALE, SMALL, SUM 64 double precision TEMP, TEST, THREE, THSJ, THVDX, TWO, TWODX 65 double precision V, V2, VPMU, X, XLOG, XPQ 66 double precision YLOG, YV, YVM1, YVP1, Z, ZERO 67 logical FLAG, J1SMAL 68 save EPS, HICUT, SMALL, XPQ, BIG, BIGLOG 69c 70 parameter( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0) 71 parameter( HALF = 0.5D0, THREE = 3.0D0) 72 parameter( C16 = 16.D0, FOUR = 4.0D0) 73 parameter(PI = 3.14159265358979323846264338327950288D0) 74C 75 parameter(HALFPI = 1.57079632679489661923132169163975144D0) 76c LOGPI = ln(pi) 77 parameter(LOGPI = 1.14472988584940017414342735135305869D0) 78c C2BYPI = 2/pi 79 parameter(C2BYPI = 0.63661977236758134307553505349005744D0) 80c EC = Euler's constant. 81 parameter(EC = 0.57721566490153286060651209008240243D0) 82c LOGTWO = Ln(2) 83 parameter(LOGTWO = 0.69314718055994530941723212145817657D0) 84 parameter( C1P5 = 1.5D0 ) 85 parameter( C11293 = 1.1293D0) 86 parameter( CP6 = 0.60206D0) 87 parameter( C59 = 0.59D0) 88 parameter( CUTLOW = 0.012D0) 89 parameter( THSJ = 0.12D0) 90 data EPS / ZERO / 91C 92 data D2SER / +.3674669051966159615185D+00, 93 * -.1782903980807269842231D+01, 94 * +.9411644685122855908427D+00, 95 * -.1958865250248747878077D+01, 96 * +.1557306621108283294475D+01, 97 * -.2521051413546812096437D+01, 98 * +.2184502685635110914145D+01, 99 * -.3140067153452674402872D+01, 100 * +.2817401038921461361158D+01, 101 * -.3772198775799670081858D+01, 102 * +.3452780604492584575060D+01/ 103C 104c ------------------------------------------------------------------ 105C 106c Set environmental parameters. 107c 108 if ( EPS .eq. ZERO) then 109 EPS = D1MACH(3) 110 HICUT = ONE / (EPS*C16) 111 SMALL = C16 * D1MACH(1) / EPS 112 XPQ = C11293 * (CP6 - LOG10( EPS )) - C59 113 BIG = D1MACH(2) / TWO 114 BIGLOG = LOG(BIG) 115 end if 116c ------------------------------------------------------------------ 117c 118c Compute V, NMIN, and NMAX. 119c 120 NMIN = INT(ALPHA) 121 V = ALPHA - DBLE(NMIN) 122 NMAX = NMIN + NUM -1 123c ------------------------------------------------------------------ 124c 125c Test validity of given argument values. 126C 127 if ( X .lt. ZERO .or. 128 * ALPHA .lt. ZERO .or. NUM .lt. 1) then 129C Error 1. 130 call ERMSG('DBESYN',1,0, 131 * 'Require X .gt. 0, ALPHA .ge. 0, NUM .ge. 1',',') 132 go to 400 133c 134c ------------------------------------------------------------------ 135C 136c Branch on size of X. 137c 138 else if (X .eq. ZERO ) then 139c Error 6. 140 do 20 I = 1, NUM 141 BY(I) = -BIG 142 20 continue 143 call ERMSG('DBESYN',6,0, 144 * 'When X = 0., function value is -INFINITY.', ',') 145 go to 400 146 else if ( X .lt. EPS ) then 147c 148c ********************* Code for very small X case. ******************** 149c 150c Use a single term expression for Y, valid for X very close to zero. 151c Ref NBS AMS 55 Eqs 9.1.8 & 9.1.9. 152c For GNU = 0, Y = (2/pi) * (EC + Ln(X/2)), {EC = Euler's const.} 153c For GNU .gt. 0 Y = -(1/pi) * Gamma(GNU) * (X/2)**(-GNU) 154c 155 XLOG = log( X ) 156 GNU = ALPHA 157c 158 do 140 I = 1, NUM 159 if ( GNU .eq. ZERO ) then 160 BY(I) = C2BYPI * (EC + XLOG - LOGTWO) 161 else 162 YLOG = DLGAMA(GNU) - GNU * (XLOG-LOGTWO) - LOGPI 163 if (YLOG .lt. BIGLOG) then 164 BY(I) = -EXP(YLOG) 165 else 166c Error 5. 167 do 120 II = I,NUM 168 BY(II) = -BIG 169 120 continue 170 call ERMSG('DBESYN',5,0, 171 * 'Results exceed overflow limit from BY(I) on.', ',') 172 call IERV1('I', I, ',') 173 go to 400 174 end if 175 end if 176 GNU = GNU + ONE 177 140 continue 178 return 179 else 180 TWODX = TWO / X 181 if ( X .le. XPQ ) then 182c 183c ********************* Code for the middle X case. ******************** 184c 185C 186C J-TYPE BESSEL FUNCTIONS FOLLOW THE RECURRENCE RELATION 187C F(V-1,X)=(2*V/X)*F(V,X)-F(V+1,X). 188C 189 MU = INT(X) + 1 190 DR = TWODX * (V+DBLE(MU)) 191 FKP1 = ONE 192 FK = ZERO 193C 194C RECUR FORWARD UNTIL FKP1 IS GREATER THAN PRECISION OF ARITHMETIC. 195C 196 210 if (EPS * ABS(FKP1) .LE. ONE) then 197 MU = MU + 1 198 DR = DR + TWODX 199 FKM1 = FK 200 FK = FKP1 201 FKP1 = DR * FK - FKM1 202 go to 210 203 end if 204C 205C WE ARE NOW ASSURED THAT BACKWARD RECURRENCE FROM MU WILL YIELD 206C ACCURATE RESULTS. 207C 208C GUARANTEE EVEN MU 209 if (MOD(MU,2) .NE. 0) MU = MU + 1 210 MUSAVE = MU 211c 212c Test for Error 3. 213c This error should never happen. Large MU would be due to 214c large X. But X is not larger than XPQ here. 215c See explanation at the beginning of this subroutine 216c of the relation of XPQ and LDIMA to the machine EPS. 217c 218 if ( MU + 1 .gt. LDIMA) then 219 call ERMSG('DBESYN', 3, 0, 220 * 'Need larger dimension, LDIMA, to process given X.', ',') 221 call ERMOR('Require LDIMA .ge. MU + 1', ',') 222 call IERV1('MU', MU, ',') 223 call IERV1('LDIMA', LDIMA, ',') 224 go to 400 225 end if 226c 227 FVM1 = SMALL 228 AJ(MU+1) = FVM1 229 FV = ZERO 230 ETA = ONE 231 SUM = FVM1 232 M = MU / 2 233 EM = DBLE(M) 234 EMU = DBLE(MU) 235 FAC = (V + EMU) * TWODX 236C 237c Set TEST = largest value that can be multiplied by 238c FAC without risking overflow. The present value of 239c FAC is the largest that will occur during the recursion. 240c TEST will be used to protect against overflow during 241c the recursion. 242c 243 TEST = BIG / MAX(ONE, FAC) 244C 245C Loop while MU .gt. ZERO 246C 247 230 continue 248 FVP1 = FV 249 FV = FVM1 250 if ( ABS(FV) .gt. TEST ) then 251c Rescale 252 FV = FV / SUM 253 FVP1 = FVP1 / SUM 254 do 240 II = MU+1, MUSAVE 255 AJ(II) = AJ(II) / SUM 256 240 continue 257 SUM = ONE 258 end if 259 FVM1 = FAC * FV - FVP1 260 MU = MU -1 261 EMU = EMU - ONE 262 FAC = (V + EMU) * TWODX 263 AJ(MU+1) = FVM1 264 if (MOD(MU,2) .eq. 0) then 265 if (V .eq. ZERO) then 266 SUM = SUM + FVM1 267 if (MU .eq. 0) then 268 SCALE = ONE / SUM 269 go to 260 270 end if 271 SUM = SUM + FVM1 272 else 273 if (MU .NE. 0) then 274 VPMU = V + EMU 275 ETA = ETA * (EM/(V+(EM-ONE)))*(VPMU/(VPMU+TWO)) 276 SUM = SUM + FVM1 * ETA 277 EM = EM - ONE 278 else 279c 280c Here MU = 0 and EM = 0NE. Thus the expression for 281c updating ETA reduces to the following simpler 282c expression. 283c 284 ETA = ETA / (V + TWO) 285 SUM = SUM + FVM1 * ETA 286 BGAM = DGAMMA(V+ONE) 287 Q2DXPV = TWODX ** V 288 SCALE = ( BGAM / ETA ) * SUM * Q2DXPV 289 SCALE = ONE / SCALE 290 go to 260 291 end if 292 end if 293 end if 294 go to 230 295 260 continue 296C 297C NORMALIZE AJ() TO GET VALUES OF J-BESSEL FUNCTION. 298C 299 do 270 I = 1, MUSAVE+1 300 AJ(I) = AJ(I) * SCALE 301 270 continue 302 MU = MUSAVE 303c 304c Compute Y Bessel functions, making use of previously computed J 305c Bessel functions and other previously computed values, MU, BGAM, 306c Q2DXPV, TWODX, V. 307c 308c Here V is in the range [0.,1.). The quantities G0 and G1 depend 309c on X and V and are unbounded as V approaches 1. Therefore we 310c make the the change of variables 311c V2 = V if V .le. 0.5 and 312c V2 = 1 - V if V .gt. 0.5 313c Then G0 and G1 are computed as functions of X and V2 with V2 in 314c the range (-0.5, 0.5]. 315C 316c Compute G0 and G1. 317c 318 V2 = V 319 if ( V .eq. ZERO ) then 320 Z = EC - log(TWODX) 321 G0 = Z / HALFPI 322 G1 = TWO / HALFPI 323 BGAMSQ = ONE 324 Q2DXPV = ONE 325 else 326 BGAMSQ = BGAM**2 327 if (V .gt. HALF) then 328c 329c Use the transformed variable, V2 = V - 1. 330c Make corresponding transformation of Q2DXPV & BGAMSQ. 331c 332 V2 = (V - HALF) - HALF 333 Q2DXPV = Q2DXPV / TWODX 334 BGAMSQ = BGAMSQ / V**2 335 end if 336 PIV2 = PI * V2 337c 338c Here V2 is in [-.5, .5]. Now test against CUTLOW = 0.012 339c 340 if ( ABS(V2) .lt. CUTLOW ) then 341C 342c Here we compute 343c G0 = (ONE / TAN(PIV2)) - Q2DXPV**2 * BGAMSQ / PIV2 344c by a formulation that retains accuracy for 345c V2 close to zero. 346c > The no. of coeffs from D2SER() used to compute D2VAL 347c could be fewer on lower precision computers, however this 348c computation is only done about 2.4% of the time so the 349c potential time saving would probably not be noticeable. 350c 351c This method was derived by C. Lawson and W. V. Snyder, 352c JPL, 1984 Apr 15. 353c 354c First compute EM1 = (2/X)**(2*V2) - 1 355c = exp(2 * V2 * log(2/X)) - 1 356c 357 ARG = TWO * V2 * log( TWODX ) 358 if ( ABS( ARG ) .lt. LOGTWO ) then 359 TEMP = TANH( HALF * ARG ) 360 EM1 = TWO * TEMP / (ONE - TEMP) 361 else 362 EM1 = EXP( ARG ) - ONE 363 end if 364c 365c Evaluate taylor series for 366c D2VAL = (PIV2 * cotan(PIV2) - BGAMSQ) / PIV2 367c 368 D2VAL = D2SER(ND2) 369 do 280 I = ND2-1, 0, -1 370 D2VAL = D2SER(I) + V2 * D2VAL 371 280 continue 372c 373 G0 = D2VAL - BGAMSQ * (EM1 / PIV2) 374 G1 = (Q2DXPV**2/HALFPI) * BGAMSQ * (TWO+V2) / (ONE-V2) 375 else 376 G0 = (ONE / TAN(PIV2)) - Q2DXPV**2 * BGAMSQ / PIV2 377 G1 = (Q2DXPV**2/HALFPI) * BGAMSQ * (TWO+V2) / (ONE-V2) 378 end if 379 end if 380c ---------------------------------- 381c 382C COMPUTE YO FROM SUM(J'S) FORM 383c 384 EN3 = V2 + ONE 385 EN2 = V2 + EN3 386 EN1 = V2 + FOUR 387 D1 = TWO 388 D2 = D1 - V2 389 D3 = D1 + V2 390 FLAG = .FALSE. 391c THSJ = 0.12 392 J1SMAL = ABS(AJ(1)) .lt. THSJ 393 if ( J1SMAL .or. V2 .lt. ZERO) then 394 FLAG = .TRUE. 395C 396C Y(V2+1,X) MUST ALSO BE COMPUTED BY A SUM 397C 398 THVDX = THREE * V2 / X 399 PSIZ = -BGAMSQ * Q2DXPV**2 / (HALFPI*X) 400 PSI1 = G0 - HALF * G1 401 end if 402c 403 if (V2 .ge. ZERO) then 404 M = 3 405 YV = G0 * AJ(1) 406 if ( J1SMAL ) then 407 YVP1 = PSIZ * AJ(1) + PSI1 * AJ(2) 408 end if 409 else 410 Z = TWODX * V * AJ(1)-AJ(2) 411 YV = G0 * Z 412 M = 2 413 YVP1 = PSIZ * Z + PSI1 * AJ(1) 414 end if 415c 416 do 290 I = M,MU,2 417 YV = G1 * AJ(I) + YV 418 G = G1 419 G1 = -G1 * (EN1/D1) * (EN2/D2) * (EN3/D3) 420 EN1 = EN1 + TWO 421 EN2 = EN2 + ONE 422 EN3 = EN3 + ONE 423 D1 = D1 + ONE 424 D2 = ONE + D2 425 D3 = D3 + TWO 426 if ( FLAG ) then 427 YVP1 = YVP1 + THVDX*G*AJ(I) + HALF*(G-G1)*AJ(I+1) 428 end if 429 290 continue 430c 431 if (V2 .lt. ZERO) then 432 Z = YVP1 433 YVP1 = V * Z * TWODX - YV 434 YV = Z 435 else if ( .NOT. J1SMAL ) then 436C 437C NOW COMPUTE Y(V+1) 438C WRONSKIAN PROVIDED NOT NEAR A ZERO OF J 439C 440 YVP1 = (YV*AJ(2)-ONE/(X*HALFPI)) / AJ(1) 441 end if 442 go to 350 443 else if ( X .LE. HICUT ) then 444c 445c ********************* Code for the large X case. ********************* 446c 447C 448c > Here we have X .ge. XPQ, and V in [0.,1.). 449c The asymptotic series for 450c the auxiliary functions P and Q can be used. 451c From these we will compute Y(V,X) and Y(V+1,X) and 452c then recur forward. 453c Reference: NBS AMS 55 Eqs 9.2.5 & 9.2.6 454c 455 call DBESPQ (X,V, P,Q) 456 CHI = X - (V + HALF) * HALFPI 457 YV = sqrt(ONE / (HALFPI*X)) * (P*SIN(CHI) + Q*COS(CHI)) 458C 459 if ( NMAX .gt. 0 ) then 460 call DBESPQ (X,V+ONE, P,Q) 461 CHI = X - (V + C1P5) * HALFPI 462 YVP1 = sqrt(ONE / (HALFPI*X)) * (P*SIN(CHI) + Q*COS(CHI)) 463 end if 464 go to 350 465 else 466c Error 2. 467 call ERMSG('DBESYN', 2, 0, 468 * 'Cannot obtain any accuracy when X exceeds HICUT.', ',') 469 call DERV1('HICUT', HICUT, ',') 470 go to 400 471 end if 472 end if 473c 474 350 continue 475c Do forward recursion 476c Given YV = Y(V,X), YVP1 = Y(V+1,X), TWODX = 2/X, NMIN, NUM, NMAX = 477c NMIN + NUM -1, X, ALPHA, and BIG. Recur forward and store 478c Y(NMIN+V) thru Y(NMAX+V) in BY(1) thru BY(NUM). 479c 480 if ( NMIN .eq. 0 ) then 481 BY(1) = YV 482 if ( NMAX .gt. 0 ) then 483 BY(2) = YVP1 484 end if 485 else if ( NMIN .eq. 1 ) then 486 BY(1) = YVP1 487 end if 488c 489 if ( NMAX .gt. 1 ) then 490 G = V * TWODX 491 GMAX = G + TWODX * DBLE(NMAX-1) 492 TEST = BIG / MAX(ONE, GMAX) 493c 494c Note: In the following statement, 3-NMIN can be nonpositive. 495c 496 do 370 K = 3-NMIN, NUM 497 YVM1 = YV 498 YV = YVP1 499 if (ABS(YV) .gt. TEST) then 500c 501c The recursion has reached the overflow limit. 502c Set remaining elts of BY() to a large negative value 503c and issue error message. 504c 505 do 360 II = MAX(K, 1),NUM 506 BY(II) = -BIG 507 360 continue 508c Error 4. 509 call ERMSG('DBESYN',4,0, 510 * 'Results exceed overflow limit from BY(I) on.', ',') 511 call IERV1('I', MAX(K,1), ',') 512 go to 400 513 end if 514c 515 G = G + TWODX 516 YVP1 = G * YV - YVM1 517 if ( K .ge. 1) BY(K) = YVP1 518 370 continue 519 end if 520 return 521c Error return 522 400 continue 523 call DERV1('X',X,',') 524 call DERV1('ALPHA',ALPHA,',') 525 call IERV1('NUM',NUM,'.') 526 return 527 end 528