1*DECK ZBIRY 2 SUBROUTINE ZBIRY (ZR, ZI, ID, KODE, BIR, BII, IERR) 3C***BEGIN PROLOGUE ZBIRY 4C***PURPOSE Compute the Airy function Bi(z) or its derivative dBi/dz 5C for complex argument z. A scaling option is available 6C to help avoid overflow. 7C***LIBRARY SLATEC 8C***CATEGORY C10D 9C***TYPE COMPLEX (CBIRY-C, ZBIRY-C) 10C***KEYWORDS AIRY FUNCTION, BESSEL FUNCTION OF ORDER ONE THIRD, 11C BESSEL FUNCTION OF ORDER TWO THIRDS 12C***AUTHOR Amos, D. E., (SNL) 13C***DESCRIPTION 14C 15C ***A DOUBLE PRECISION ROUTINE*** 16C On KODE=1, ZBIRY computes the complex Airy function Bi(z) 17C or its derivative dBi/dz on ID=0 or ID=1 respectively. 18C On KODE=2, a scaling option exp(abs(Re(zeta)))*Bi(z) or 19C exp(abs(Re(zeta)))*dBi/dz is provided to remove the 20C exponential behavior in both the left and right half planes 21C where zeta=(2/3)*z**(3/2). 22C 23C The Airy functions Bi(z) and dBi/dz are analytic in the 24C whole z-plane, and the scaling option does not destroy this 25C property. 26C 27C Input 28C ZR - DOUBLE PRECISION real part of argument Z 29C ZI - DOUBLE PRECISION imag part of argument Z 30C ID - Order of derivative, ID=0 or ID=1 31C KODE - A parameter to indicate the scaling option 32C KODE=1 returns 33C BI=Bi(z) on ID=0 34C BI=dBi/dz on ID=1 35C at z=Z 36C =2 returns 37C BI=exp(abs(Re(zeta)))*Bi(z) on ID=0 38C BI=exp(abs(Re(zeta)))*dBi/dz on ID=1 39C at z=Z where zeta=(2/3)*z**(3/2) 40C 41C Output 42C BIR - DOUBLE PRECISION real part of result 43C BII - DOUBLE PRECISION imag part of result 44C IERR - Error flag 45C IERR=0 Normal return - COMPUTATION COMPLETED 46C IERR=1 Input error - NO COMPUTATION 47C IERR=2 Overflow - NO COMPUTATION 48C (Re(Z) too large with KODE=1) 49C IERR=3 Precision warning - COMPUTATION COMPLETED 50C (Result has less than half precision) 51C IERR=4 Precision error - NO COMPUTATION 52C (Result has no precision) 53C IERR=5 Algorithmic error - NO COMPUTATION 54C (Termination condition not met) 55C 56C *Long Description: 57C 58C Bi(z) and dBi/dz are computed from I Bessel functions by 59C 60C Bi(z) = c*sqrt(z)*( I(-1/3,zeta) + I(1/3,zeta) ) 61C dBi/dz = c* z *( I(-2/3,zeta) + I(2/3,zeta) ) 62C c = 1/sqrt(3) 63C zeta = (2/3)*z**(3/2) 64C 65C when abs(z)>1 and from power series when abs(z)<=1. 66C 67C In most complex variable computation, one must evaluate ele- 68C mentary functions. When the magnitude of Z is large, losses 69C of significance by argument reduction occur. Consequently, if 70C the magnitude of ZETA=(2/3)*Z**(3/2) exceeds U1=SQRT(0.5/UR), 71C then losses exceeding half precision are likely and an error 72C flag IERR=3 is triggered where UR=MAX(D1MACH(4),1.0D-18) is 73C double precision unit roundoff limited to 18 digits precision. 74C Also, if the magnitude of ZETA is larger than U2=0.5/UR, then 75C all significance is lost and IERR=4. In order to use the INT 76C function, ZETA must be further restricted not to exceed 77C U3=I1MACH(9)=LARGEST INTEGER. Thus, the magnitude of ZETA 78C must be restricted by MIN(U2,U3). In IEEE arithmetic, U1,U2, 79C and U3 are approximately 2.0E+3, 4.2E+6, 2.1E+9 in single 80C precision and 4.7E+7, 2.3E+15, 2.1E+9 in double precision. 81C This makes U2 limiting is single precision and U3 limiting 82C in double precision. This means that the magnitude of Z 83C cannot exceed approximately 3.4E+4 in single precision and 84C 2.1E+6 in double precision. This also means that one can 85C expect to retain, in the worst cases on 32-bit machines, 86C no digits in single precision and only 6 digits in double 87C precision. 88C 89C The approximate relative error in the magnitude of a complex 90C Bessel function can be expressed as P*10**S where P=MAX(UNIT 91C ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre- 92C sents the increase in error due to argument reduction in the 93C elementary functions. Here, S=MAX(1,ABS(LOG10(ABS(Z))), 94C ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF 95C ABS(Z),ABS(EXPONENT OF FNU)) ). However, the phase angle may 96C have only absolute accuracy. This is most likely to occur 97C when one component (in magnitude) is larger than the other by 98C several orders of magnitude. If one component is 10**K larger 99C than the other, then one can expect only MAX(ABS(LOG10(P))-K, 100C 0) significant digits; or, stated another way, when K exceeds 101C the exponent of P, no significant digits remain in the smaller 102C component. However, the phase angle retains absolute accuracy 103C because, in complex arithmetic with precision P, the smaller 104C component will not (as a rule) decrease below P times the 105C magnitude of the larger component. In these extreme cases, 106C the principal phase angle is on the order of +P, -P, PI/2-P, 107C or -PI/2+P. 108C 109C***REFERENCES 1. M. Abramowitz and I. A. Stegun, Handbook of Mathe- 110C matical Functions, National Bureau of Standards 111C Applied Mathematics Series 55, U. S. Department 112C of Commerce, Tenth Printing (1972) or later. 113C 2. D. E. Amos, Computation of Bessel Functions of 114C Complex Argument and Large Order, Report SAND83-0643, 115C Sandia National Laboratories, Albuquerque, NM, May 116C 1983. 117C 3. D. E. Amos, A Subroutine Package for Bessel Functions 118C of a Complex Argument and Nonnegative Order, Report 119C SAND85-1018, Sandia National Laboratory, Albuquerque, 120C NM, May 1985. 121C 4. D. E. Amos, A portable package for Bessel functions 122C of a complex argument and nonnegative order, ACM 123C Transactions on Mathematical Software, 12 (September 124C 1986), pp. 265-273. 125C 126C***ROUTINES CALLED D1MACH, I1MACH, ZABS, ZBINU, ZDIV, ZSQRT 127C***REVISION HISTORY (YYMMDD) 128C 830501 DATE WRITTEN 129C 890801 REVISION DATE from Version 3.2 130C 910415 Prologue converted to Version 4.0 format. (BAB) 131C 920128 Category corrected. (WRB) 132C 920811 Prologue revised. (DWL) 133C 930122 Added ZSQRT to EXTERNAL statement. (RWC) 134C***END PROLOGUE ZBIRY 135C COMPLEX BI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3 136 DOUBLE PRECISION AA, AD, AK, ALIM, ATRM, AZ, AZ3, BB, BII, BIR, 137 * BK, CC, CK, COEF, CONEI, CONER, CSQI, CSQR, CYI, CYR, C1, C2, 138 * DIG, DK, D1, D2, EAA, ELIM, FID, FMR, FNU, FNUL, PI, RL, R1M5, 139 * SFAC, STI, STR, S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, 140 * TRM2R, TTH, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, ZABS 141 INTEGER ID, IERR, K, KODE, K1, K2, NZ, I1MACH 142 DIMENSION CYR(2), CYI(2) 143 EXTERNAL ZABS, ZSQRT 144 DATA TTH, C1, C2, COEF, PI /6.66666666666666667D-01, 145 * 6.14926627446000736D-01,4.48288357353826359D-01, 146 * 5.77350269189625765D-01,3.14159265358979324D+00/ 147 DATA CONER, CONEI /1.0D0,0.0D0/ 148C***FIRST EXECUTABLE STATEMENT ZBIRY 149 IERR = 0 150 NZ=0 151 IF (ID.LT.0 .OR. ID.GT.1) IERR=1 152 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 153 IF (IERR.NE.0) RETURN 154 AZ = ZABS(ZR,ZI) 155 TOL = MAX(D1MACH(4),1.0D-18) 156 FID = ID 157 IF (AZ.GT.1.0E0) GO TO 70 158C----------------------------------------------------------------------- 159C POWER SERIES FOR ABS(Z).LE.1. 160C----------------------------------------------------------------------- 161 S1R = CONER 162 S1I = CONEI 163 S2R = CONER 164 S2I = CONEI 165 IF (AZ.LT.TOL) GO TO 130 166 AA = AZ*AZ 167 IF (AA.LT.TOL/AZ) GO TO 40 168 TRM1R = CONER 169 TRM1I = CONEI 170 TRM2R = CONER 171 TRM2I = CONEI 172 ATRM = 1.0D0 173 STR = ZR*ZR - ZI*ZI 174 STI = ZR*ZI + ZI*ZR 175 Z3R = STR*ZR - STI*ZI 176 Z3I = STR*ZI + STI*ZR 177 AZ3 = AZ*AA 178 AK = 2.0D0 + FID 179 BK = 3.0D0 - FID - FID 180 CK = 4.0D0 - FID 181 DK = 3.0D0 + FID + FID 182 D1 = AK*DK 183 D2 = BK*CK 184 AD = MIN(D1,D2) 185 AK = 24.0D0 + 9.0D0*FID 186 BK = 30.0D0 - 9.0D0*FID 187 DO 30 K=1,25 188 STR = (TRM1R*Z3R-TRM1I*Z3I)/D1 189 TRM1I = (TRM1R*Z3I+TRM1I*Z3R)/D1 190 TRM1R = STR 191 S1R = S1R + TRM1R 192 S1I = S1I + TRM1I 193 STR = (TRM2R*Z3R-TRM2I*Z3I)/D2 194 TRM2I = (TRM2R*Z3I+TRM2I*Z3R)/D2 195 TRM2R = STR 196 S2R = S2R + TRM2R 197 S2I = S2I + TRM2I 198 ATRM = ATRM*AZ3/AD 199 D1 = D1 + AK 200 D2 = D2 + BK 201 AD = MIN(D1,D2) 202 IF (ATRM.LT.TOL*AD) GO TO 40 203 AK = AK + 18.0D0 204 BK = BK + 18.0D0 205 30 CONTINUE 206 40 CONTINUE 207 IF (ID.EQ.1) GO TO 50 208 BIR = C1*S1R + C2*(ZR*S2R-ZI*S2I) 209 BII = C1*S1I + C2*(ZR*S2I+ZI*S2R) 210 IF (KODE.EQ.1) RETURN 211 CALL ZSQRT(ZR, ZI, STR, STI) 212 ZTAR = TTH*(ZR*STR-ZI*STI) 213 ZTAI = TTH*(ZR*STI+ZI*STR) 214 AA = ZTAR 215 AA = -ABS(AA) 216 EAA = EXP(AA) 217 BIR = BIR*EAA 218 BII = BII*EAA 219 RETURN 220 50 CONTINUE 221 BIR = S2R*C2 222 BII = S2I*C2 223 IF (AZ.LE.TOL) GO TO 60 224 CC = C1/(1.0D0+FID) 225 STR = S1R*ZR - S1I*ZI 226 STI = S1R*ZI + S1I*ZR 227 BIR = BIR + CC*(STR*ZR-STI*ZI) 228 BII = BII + CC*(STR*ZI+STI*ZR) 229 60 CONTINUE 230 IF (KODE.EQ.1) RETURN 231 CALL ZSQRT(ZR, ZI, STR, STI) 232 ZTAR = TTH*(ZR*STR-ZI*STI) 233 ZTAI = TTH*(ZR*STI+ZI*STR) 234 AA = ZTAR 235 AA = -ABS(AA) 236 EAA = EXP(AA) 237 BIR = BIR*EAA 238 BII = BII*EAA 239 RETURN 240C----------------------------------------------------------------------- 241C CASE FOR ABS(Z).GT.1.0 242C----------------------------------------------------------------------- 243 70 CONTINUE 244 FNU = (1.0D0+FID)/3.0D0 245C----------------------------------------------------------------------- 246C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 247C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. 248C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 249C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 250C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 251C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 252C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 253C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 254C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. 255C----------------------------------------------------------------------- 256 K1 = I1MACH(15) 257 K2 = I1MACH(16) 258 R1M5 = D1MACH(5) 259 K = MIN(ABS(K1),ABS(K2)) 260 ELIM = 2.303D0*(K*R1M5-3.0D0) 261 K1 = I1MACH(14) - 1 262 AA = R1M5*K1 263 DIG = MIN(AA,18.0D0) 264 AA = AA*2.303D0 265 ALIM = ELIM + MAX(-AA,-41.45D0) 266 RL = 1.2D0*DIG + 3.0D0 267 FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) 268C----------------------------------------------------------------------- 269C TEST FOR RANGE 270C----------------------------------------------------------------------- 271 AA=0.5D0/TOL 272 BB=I1MACH(9)*0.5D0 273 AA=MIN(AA,BB) 274 AA=AA**TTH 275 IF (AZ.GT.AA) GO TO 260 276 AA=SQRT(AA) 277 IF (AZ.GT.AA) IERR=3 278 CALL ZSQRT(ZR, ZI, CSQR, CSQI) 279 ZTAR = TTH*(ZR*CSQR-ZI*CSQI) 280 ZTAI = TTH*(ZR*CSQI+ZI*CSQR) 281C----------------------------------------------------------------------- 282C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL 283C----------------------------------------------------------------------- 284 SFAC = 1.0D0 285 AK = ZTAI 286 IF (ZR.GE.0.0D0) GO TO 80 287 BK = ZTAR 288 CK = -ABS(BK) 289 ZTAR = CK 290 ZTAI = AK 291 80 CONTINUE 292 IF (ZI.NE.0.0D0 .OR. ZR.GT.0.0D0) GO TO 90 293 ZTAR = 0.0D0 294 ZTAI = AK 295 90 CONTINUE 296 AA = ZTAR 297 IF (KODE.EQ.2) GO TO 100 298C----------------------------------------------------------------------- 299C OVERFLOW TEST 300C----------------------------------------------------------------------- 301 BB = ABS(AA) 302 IF (BB.LT.ALIM) GO TO 100 303 BB = BB + 0.25D0*LOG(AZ) 304 SFAC = TOL 305 IF (BB.GT.ELIM) GO TO 190 306 100 CONTINUE 307 FMR = 0.0D0 308 IF (AA.GE.0.0D0 .AND. ZR.GT.0.0D0) GO TO 110 309 FMR = PI 310 IF (ZI.LT.0.0D0) FMR = -PI 311 ZTAR = -ZTAR 312 ZTAI = -ZTAI 313 110 CONTINUE 314C----------------------------------------------------------------------- 315C AA=FACTOR FOR ANALYTIC CONTINUATION OF I(FNU,ZTA) 316C KODE=2 RETURNS EXP(-ABS(XZTA))*I(FNU,ZTA) FROM CBESI 317C----------------------------------------------------------------------- 318 CALL ZBINU(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, RL, FNUL, TOL, 319 * ELIM, ALIM) 320 IF (NZ.LT.0) GO TO 200 321 AA = FMR*FNU 322 Z3R = SFAC 323 STR = COS(AA) 324 STI = SIN(AA) 325 S1R = (STR*CYR(1)-STI*CYI(1))*Z3R 326 S1I = (STR*CYI(1)+STI*CYR(1))*Z3R 327 FNU = (2.0D0-FID)/3.0D0 328 CALL ZBINU(ZTAR, ZTAI, FNU, KODE, 2, CYR, CYI, NZ, RL, FNUL, TOL, 329 * ELIM, ALIM) 330 CYR(1) = CYR(1)*Z3R 331 CYI(1) = CYI(1)*Z3R 332 CYR(2) = CYR(2)*Z3R 333 CYI(2) = CYI(2)*Z3R 334C----------------------------------------------------------------------- 335C BACKWARD RECUR ONE STEP FOR ORDERS -1/3 OR -2/3 336C----------------------------------------------------------------------- 337 CALL ZDIV(CYR(1), CYI(1), ZTAR, ZTAI, STR, STI) 338 S2R = (FNU+FNU)*STR + CYR(2) 339 S2I = (FNU+FNU)*STI + CYI(2) 340 AA = FMR*(FNU-1.0D0) 341 STR = COS(AA) 342 STI = SIN(AA) 343 S1R = COEF*(S1R+S2R*STR-S2I*STI) 344 S1I = COEF*(S1I+S2R*STI+S2I*STR) 345 IF (ID.EQ.1) GO TO 120 346 STR = CSQR*S1R - CSQI*S1I 347 S1I = CSQR*S1I + CSQI*S1R 348 S1R = STR 349 BIR = S1R/SFAC 350 BII = S1I/SFAC 351 RETURN 352 120 CONTINUE 353 STR = ZR*S1R - ZI*S1I 354 S1I = ZR*S1I + ZI*S1R 355 S1R = STR 356 BIR = S1R/SFAC 357 BII = S1I/SFAC 358 RETURN 359 130 CONTINUE 360 AA = C1*(1.0D0-FID) + FID*C2 361 BIR = AA 362 BII = 0.0D0 363 RETURN 364 190 CONTINUE 365 IERR=2 366 NZ=0 367 RETURN 368 200 CONTINUE 369 IF(NZ.EQ.(-1)) GO TO 190 370 NZ=0 371 IERR=5 372 RETURN 373 260 CONTINUE 374 IERR=4 375 NZ=0 376 RETURN 377 END 378