1 SUBROUTINE ZAIRY(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR) 2C***BEGIN PROLOGUE ZAIRY 3C***DATE WRITTEN 830501 (YYMMDD) 4C***REVISION DATE 890801 (YYMMDD) 5C***CATEGORY NO. B5K 6C***KEYWORDS AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD 7C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 8C***PURPOSE TO COMPUTE AIRY FUNCTIONS AI(Z) AND DAI(Z) FOR COMPLEX Z 9C***DESCRIPTION 10C 11C ***A DOUBLE PRECISION ROUTINE*** 12C ON KODE=1, ZAIRY COMPUTES THE COMPLEX AIRY FUNCTION AI(Z) OR 13C ITS DERIVATIVE DAI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON 14C KODE=2, A SCALING OPTION CEXP(ZTA)*AI(Z) OR CEXP(ZTA)* 15C DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN 16C -PI/3.LT.ARG(Z).LT.PI/3 AND THE EXPONENTIAL GROWTH IN 17C PI/3.LT.ABS(ARG(Z)).LT.PI WHERE ZTA=(2/3)*Z*CSQRT(Z). 18C 19C WHILE THE AIRY FUNCTIONS AI(Z) AND DAI(Z)/DZ ARE ANALYTIC IN 20C THE WHOLE Z PLANE, THE CORRESPONDING SCALED FUNCTIONS DEFINED 21C FOR KODE=2 HAVE A CUT ALONG THE NEGATIVE REAL AXIS. 22C DEFINTIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF 23C MATHEMATICAL FUNCTIONS (REF. 1). 24C 25C INPUT ZR,ZI ARE DOUBLE PRECISION 26C ZR,ZI - Z=CMPLX(ZR,ZI) 27C ID - ORDER OF DERIVATIVE, ID=0 OR ID=1 28C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 29C KODE= 1 RETURNS 30C AI=AI(Z) ON ID=0 OR 31C AI=DAI(Z)/DZ ON ID=1 32C = 2 RETURNS 33C AI=CEXP(ZTA)*AI(Z) ON ID=0 OR 34C AI=CEXP(ZTA)*DAI(Z)/DZ ON ID=1 WHERE 35C ZTA=(2/3)*Z*CSQRT(Z) 36C 37C OUTPUT AIR,AII ARE DOUBLE PRECISION 38C AIR,AII- COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND 39C KODE 40C NZ - UNDERFLOW INDICATOR 41C NZ= 0 , NORMAL RETURN 42C NZ= 1 , AI=CMPLX(0.0D0,0.0D0) DUE TO UNDERFLOW IN 43C -PI/3.LT.ARG(Z).LT.PI/3 ON KODE=1 44C IERR - ERROR FLAG 45C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 46C IERR=1, INPUT ERROR - NO COMPUTATION 47C IERR=2, OVERFLOW - NO COMPUTATION, REAL(ZTA) 48C TOO LARGE ON KODE=1 49C IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED 50C LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION 51C PRODUCE LESS THAN HALF OF MACHINE ACCURACY 52C IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION 53C COMPLETE LOSS OF ACCURACY BY ARGUMENT 54C REDUCTION 55C IERR=5, ERROR - NO COMPUTATION, 56C ALGORITHM TERMINATION CONDITION NOT MET 57C 58C***LONG DESCRIPTION 59C 60C AI AND DAI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE K BESSEL 61C FUNCTIONS BY 62C 63C AI(Z)=C*SQRT(Z)*K(1/3,ZTA) , DAI(Z)=-C*Z*K(2/3,ZTA) 64C C=1.0/(PI*SQRT(3.0)) 65C ZTA=(2/3)*Z**(3/2) 66C 67C WITH THE POWER SERIES FOR CABS(Z).LE.1.0. 68C 69C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 70C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES 71C OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF 72C THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR), 73C THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR 74C FLAG IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS 75C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. 76C ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN 77C ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT 78C FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE 79C LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA 80C MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, 81C AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE 82C PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE 83C PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT- 84C ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG- 85C NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN 86C DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN 87C EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, 88C NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE 89C PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER 90C MACHINES. 91C 92C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 93C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 94C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 95C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 96C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 97C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 98C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 99C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 100C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 101C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 102C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 103C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 104C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 105C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 106C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 107C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 108C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 109C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 110C OR -PI/2+P. 111C 112C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 113C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 114C COMMERCE, 1955. 115C 116C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 117C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 118C 119C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 120C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 121C 1018, MAY, 1985 122C 123C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 124C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 125C MATH. SOFTWARE, 1986 126C 127C***ROUTINES CALLED ZACAI,ZBKNU,ZEXP,ZSQRT,I1MACH,D1MACH 128C***END PROLOGUE ZAIRY 129C COMPLEX AI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3 130 DOUBLE PRECISION AA, AD, AII, AIR, AK, ALIM, ATRM, AZ, AZ3, BK, 131 * CC, CK, COEF, CONEI, CONER, CSQI, CSQR, CYI, CYR, C1, C2, DIG, 132 * DK, D1, D2, ELIM, FID, FNU, PTR, RL, R1M5, SFAC, STI, STR, 133 * S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, TRM2R, TTH, ZEROI, 134 * ZEROR, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, ZABS, ALAZ, BB 135 INTEGER ID, IERR, IFLAG, K, KODE, K1, K2, MR, NN, NZ, I1MACH 136 DIMENSION CYR(1), CYI(1) 137 DATA TTH, C1, C2, COEF /6.66666666666666667D-01, 138 * 3.55028053887817240D-01,2.58819403792806799D-01, 139 * 1.83776298473930683D-01/ 140 DATA ZEROR, ZEROI, CONER, CONEI /0.0D0,0.0D0,1.0D0,0.0D0/ 141C***FIRST EXECUTABLE STATEMENT ZAIRY 142 IERR = 0 143 NZ=0 144 IF (ID.LT.0 .OR. ID.GT.1) IERR=1 145 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 146 IF (IERR.NE.0) RETURN 147 AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 148 TOL = DMAX1(D1MACH(4),1.0D-18) 149 FID = DBLE(FLOAT(ID)) 150 IF (AZ.GT.1.0D0) GO TO 70 151C----------------------------------------------------------------------- 152C POWER SERIES FOR CABS(Z).LE.1. 153C----------------------------------------------------------------------- 154 S1R = CONER 155 S1I = CONEI 156 S2R = CONER 157 S2I = CONEI 158 IF (AZ.LT.TOL) GO TO 170 159 AA = AZ*AZ 160 IF (AA.LT.TOL/AZ) GO TO 40 161 TRM1R = CONER 162 TRM1I = CONEI 163 TRM2R = CONER 164 TRM2I = CONEI 165 ATRM = 1.0D0 166 STR = ZR*ZR - ZI*ZI 167 STI = ZR*ZI + ZI*ZR 168 Z3R = STR*ZR - STI*ZI 169 Z3I = STR*ZI + STI*ZR 170 AZ3 = AZ*AA 171 AK = 2.0D0 + FID 172 BK = 3.0D0 - FID - FID 173 CK = 4.0D0 - FID 174 DK = 3.0D0 + FID + FID 175 D1 = AK*DK 176 D2 = BK*CK 177 AD = DMIN1(D1,D2) 178 AK = 24.0D0 + 9.0D0*FID 179 BK = 30.0D0 - 9.0D0*FID 180 DO 30 K=1,25 181 STR = (TRM1R*Z3R-TRM1I*Z3I)/D1 182 TRM1I = (TRM1R*Z3I+TRM1I*Z3R)/D1 183 TRM1R = STR 184 S1R = S1R + TRM1R 185 S1I = S1I + TRM1I 186 STR = (TRM2R*Z3R-TRM2I*Z3I)/D2 187 TRM2I = (TRM2R*Z3I+TRM2I*Z3R)/D2 188 TRM2R = STR 189 S2R = S2R + TRM2R 190 S2I = S2I + TRM2I 191 ATRM = ATRM*AZ3/AD 192 D1 = D1 + AK 193 D2 = D2 + BK 194 AD = DMIN1(D1,D2) 195 IF (ATRM.LT.TOL*AD) GO TO 40 196 AK = AK + 18.0D0 197 BK = BK + 18.0D0 198 30 CONTINUE 199 40 CONTINUE 200 IF (ID.EQ.1) THEN 201 GO TO 50 202 END IF 203 AIR = S1R*C1 - C2*(ZR*S2R-ZI*S2I) 204 AII = S1I*C1 - C2*(ZR*S2I+ZI*S2R) 205 IF (KODE.EQ.1) RETURN 206 CALL ZSQRT(ZR, ZI, STR, STI) 207 ZTAR = TTH*(ZR*STR-ZI*STI) 208 ZTAI = TTH*(ZR*STI+ZI*STR) 209 CALL ZEXP(ZTAR, ZTAI, STR, STI) 210 PTR = AIR*STR - AII*STI 211 AII = AIR*STI + AII*STR 212 AIR = PTR 213 RETURN 214 50 CONTINUE 215 AIR = -S2R*C2 216 AII = -S2I*C2 217 IF (AZ.LE.TOL) GO TO 60 218 STR = ZR*S1R - ZI*S1I 219 STI = ZR*S1I + ZI*S1R 220 CC = C1/(1.0D0+FID) 221 AIR = AIR + CC*(STR*ZR-STI*ZI) 222 AII = AII + CC*(STR*ZI+STI*ZR) 223 60 CONTINUE 224 IF (KODE.EQ.1) RETURN 225 CALL ZSQRT(ZR, ZI, STR, STI) 226 ZTAR = TTH*(ZR*STR-ZI*STI) 227 ZTAI = TTH*(ZR*STI+ZI*STR) 228 CALL ZEXP(ZTAR, ZTAI, STR, STI) 229 PTR = STR*AIR - STI*AII 230 AII = STR*AII + STI*AIR 231 AIR = PTR 232 RETURN 233C----------------------------------------------------------------------- 234C CASE FOR CABS(Z).GT.1.0 235C----------------------------------------------------------------------- 236 70 CONTINUE 237 FNU = (1.0D0+FID)/3.0D0 238C----------------------------------------------------------------------- 239C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 240C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0D-18. 241C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 242C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 243C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 244C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 245C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 246C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 247C----------------------------------------------------------------------- 248 K1 = I1MACH(15) 249 K2 = I1MACH(16) 250 R1M5 = D1MACH(5) 251 K = MIN0(IABS(K1),IABS(K2)) 252 ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) 253 K1 = I1MACH(14) - 1 254 AA = R1M5*DBLE(FLOAT(K1)) 255 DIG = DMIN1(AA,18.0D0) 256 AA = AA*2.303D0 257 ALIM = ELIM + DMAX1(-AA,-41.45D0) 258 RL = 1.2D0*DIG + 3.0D0 259 ALAZ = DLOG(AZ) 260C-------------------------------------------------------------------------- 261C TEST FOR PROPER RANGE 262C----------------------------------------------------------------------- 263 AA=0.5D0/TOL 264 BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 265 AA=DMIN1(AA,BB) 266 AA=AA**TTH 267 IF (AZ.GT.AA) GO TO 260 268 AA=DSQRT(AA) 269 IF (AZ.GT.AA) IERR=3 270 CALL ZSQRT(ZR, ZI, CSQR, CSQI) 271 ZTAR = TTH*(ZR*CSQR-ZI*CSQI) 272 ZTAI = TTH*(ZR*CSQI+ZI*CSQR) 273C----------------------------------------------------------------------- 274C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL 275C----------------------------------------------------------------------- 276 IFLAG = 0 277 SFAC = 1.0D0 278 AK = ZTAI 279 IF (ZR.GE.0.0D0) GO TO 80 280 BK = ZTAR 281 CK = -DABS(BK) 282 ZTAR = CK 283 ZTAI = AK 284 80 CONTINUE 285 IF (ZI.NE.0.0D0) GO TO 90 286 IF (ZR.GT.0.0D0) GO TO 90 287 ZTAR = 0.0D0 288 ZTAI = AK 289 90 CONTINUE 290 AA = ZTAR 291 IF (AA.GE.0.0D0 .AND. ZR.GT.0.0D0) GO TO 110 292 IF (KODE.EQ.2) GO TO 100 293C----------------------------------------------------------------------- 294C OVERFLOW TEST 295C----------------------------------------------------------------------- 296 IF (AA.GT.(-ALIM)) GO TO 100 297 AA = -AA + 0.25D0*ALAZ 298 IFLAG = 1 299 SFAC = TOL 300 IF (AA.GT.ELIM) GO TO 270 301 100 CONTINUE 302C----------------------------------------------------------------------- 303C CBKNU AND CACON RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2 304C----------------------------------------------------------------------- 305 MR = 1 306 IF (ZI.LT.0.0D0) MR = -1 307 CALL ZACAI(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, NN, RL, TOL, 308 * ELIM, ALIM) 309 IF (NN.LT.0) GO TO 280 310 NZ = NZ + NN 311 GO TO 130 312 110 CONTINUE 313 IF (KODE.EQ.2) GO TO 120 314C----------------------------------------------------------------------- 315C UNDERFLOW TEST 316C----------------------------------------------------------------------- 317 IF (AA.LT.ALIM) GO TO 120 318 AA = -AA - 0.25D0*ALAZ 319 IFLAG = 2 320 SFAC = 1.0D0/TOL 321 IF (AA.LT.(-ELIM)) GO TO 210 322 120 CONTINUE 323 CALL ZBKNU(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM, 324 * ALIM) 325 130 CONTINUE 326 S1R = CYR(1)*COEF 327 S1I = CYI(1)*COEF 328 IF (IFLAG.NE.0) GO TO 150 329 IF (ID.EQ.1) GO TO 140 330 AIR = CSQR*S1R - CSQI*S1I 331 AII = CSQR*S1I + CSQI*S1R 332 RETURN 333 140 CONTINUE 334 AIR = -(ZR*S1R-ZI*S1I) 335 AII = -(ZR*S1I+ZI*S1R) 336 RETURN 337 150 CONTINUE 338 S1R = S1R*SFAC 339 S1I = S1I*SFAC 340 IF (ID.EQ.1) GO TO 160 341 STR = S1R*CSQR - S1I*CSQI 342 S1I = S1R*CSQI + S1I*CSQR 343 S1R = STR 344 AIR = S1R/SFAC 345 AII = S1I/SFAC 346 RETURN 347 160 CONTINUE 348 STR = -(S1R*ZR-S1I*ZI) 349 S1I = -(S1R*ZI+S1I*ZR) 350 S1R = STR 351 AIR = S1R/SFAC 352 AII = S1I/SFAC 353 RETURN 354 170 CONTINUE 355 AA = 1.0D+3*D1MACH(1) 356 S1R = ZEROR 357 S1I = ZEROI 358 IF (ID.EQ.1) GO TO 190 359 IF (AZ.LE.AA) GO TO 180 360 S1R = C2*ZR 361 S1I = C2*ZI 362 180 CONTINUE 363 AIR = C1 - S1R 364 AII = -S1I 365 RETURN 366 190 CONTINUE 367 AIR = -C2 368 AII = 0.0D0 369 AA = DSQRT(AA) 370 IF (AZ.LE.AA) GO TO 200 371 S1R = 0.5D0*(ZR*ZR-ZI*ZI) 372 S1I = ZR*ZI 373 200 CONTINUE 374 AIR = AIR + C1*S1R 375 AII = AII + C1*S1I 376 RETURN 377 210 CONTINUE 378 NZ = 1 379 AIR = ZEROR 380 AII = ZEROI 381 RETURN 382 270 CONTINUE 383 NZ = 0 384 IERR=2 385 RETURN 386 280 CONTINUE 387 IF(NN.EQ.(-1)) GO TO 270 388 NZ=0 389 IERR=5 390 RETURN 391 260 CONTINUE 392 IERR=4 393 NZ=0 394 RETURN 395 END 396