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