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