00001 SUBROUTINE CBIRY(Z, ID, KODE, BI, IERR) 00002 C***BEGIN PROLOGUE CBIRY 00003 C***DATE WRITTEN 830501 (YYMMDD) 00004 C***REVISION DATE 890801 (YYMMDD) 00005 C***CATEGORY NO. B5K 00006 C***KEYWORDS AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD 00007 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 00008 C***PURPOSE TO COMPUTE AIRY FUNCTIONS BI(Z) AND DBI(Z) FOR COMPLEX Z 00009 C***DESCRIPTION 00010 C 00011 C ON KODE=1, CBIRY COMPUTES THE COMPLEX AIRY FUNCTION BI(Z) OR 00012 C ITS DERIVATIVE DBI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON 00013 C KODE=2, A SCALING OPTION CEXP(-AXZTA)*BI(Z) OR CEXP(-AXZTA)* 00014 C DBI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL BEHAVIOR IN 00015 C BOTH THE LEFT AND RIGHT HALF PLANES WHERE 00016 C ZTA=(2/3)*Z*CSQRT(Z)=CMPLX(XZTA,YZTA) AND AXZTA=ABS(XZTA). 00017 C DEFINITIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF 00018 C MATHEMATICAL FUNCTIONS (REF. 1). 00019 C 00020 C INPUT 00021 C Z - Z=CMPLX(X,Y) 00022 C ID - ORDER OF DERIVATIVE, ID=0 OR ID=1 00023 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 00024 C KODE= 1 RETURNS 00025 C BI=BI(Z) ON ID=0 OR 00026 C BI=DBI(Z)/DZ ON ID=1 00027 C = 2 RETURNS 00028 C BI=CEXP(-AXZTA)*BI(Z) ON ID=0 OR 00029 C BI=CEXP(-AXZTA)*DBI(Z)/DZ ON ID=1 WHERE 00030 C ZTA=(2/3)*Z*CSQRT(Z)=CMPLX(XZTA,YZTA) 00031 C AND AXZTA=ABS(XZTA) 00032 C 00033 C OUTPUT 00034 C BI - COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND 00035 C KODE 00036 C IERR - ERROR FLAG 00037 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 00038 C IERR=1, INPUT ERROR - NO COMPUTATION 00039 C IERR=2, OVERFLOW - NO COMPUTATION, REAL(Z) 00040 C TOO LARGE WITH KODE=1 00041 C IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED 00042 C LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION 00043 C PRODUCE LESS THAN HALF OF MACHINE ACCURACY 00044 C IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION 00045 C COMPLETE LOSS OF ACCURACY BY ARGUMENT 00046 C REDUCTION 00047 C IERR=5, ERROR - NO COMPUTATION, 00048 C ALGORITHM TERMINATION CONDITION NOT MET 00049 C 00050 C***LONG DESCRIPTION 00051 C 00052 C BI AND DBI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE I BESSEL 00053 C FUNCTIONS BY 00054 C 00055 C BI(Z)=C*SQRT(Z)*( I(-1/3,ZTA) + I(1/3,ZTA) ) 00056 C DBI(Z)=C * Z * ( I(-2/3,ZTA) + I(2/3,ZTA) ) 00057 C C=1.0/SQRT(3.0) 00058 C ZTA=(2/3)*Z**(3/2) 00059 C 00060 C WITH THE POWER SERIES FOR CABS(Z).LE.1.0. 00061 C 00062 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 00063 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES 00064 C OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF 00065 C THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR), 00066 C THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR 00067 C FLAG IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. 00068 C ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN 00069 C ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT 00070 C FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE 00071 C LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA 00072 C MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, 00073 C AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE 00074 C PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE 00075 C PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT- 00076 C ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG- 00077 C NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN 00078 C DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN 00079 C EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, 00080 C NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE 00081 C PRECISION ARITHMETIC. 00082 C 00083 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 00084 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 00085 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 00086 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 00087 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 00088 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 00089 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 00090 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 00091 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 00092 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 00093 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 00094 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 00095 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 00096 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 00097 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 00098 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 00099 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 00100 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 00101 C OR -PI/2+P. 00102 C 00103 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 00104 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 00105 C COMMERCE, 1955. 00106 C 00107 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00108 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 00109 C 00110 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 00111 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 00112 C 1018, MAY, 1985 00113 C 00114 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 00115 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 00116 C MATH. SOFTWARE, 1986 00117 C 00118 C***ROUTINES CALLED CBINU,I1MACH,R1MACH 00119 C***END PROLOGUE CBIRY 00120 COMPLEX BI, CONE, CSQ, CY, S1, S2, TRM1, TRM2, Z, ZTA, Z3 00121 REAL AA, AD, AK, ALIM, ATRM, AZ, AZ3, BB, BK, CK, COEF, C1, C2, 00122 * DIG, DK, D1, D2, ELIM, FID, FMR, FNU, FNUL, PI, RL, R1M5, SFAC, 00123 * TOL, TTH, ZI, ZR, Z3I, Z3R, R1MACH 00124 INTEGER ID, IERR, K, KODE, K1, K2, NZ, I1MACH 00125 DIMENSION CY(2) 00126 DATA TTH, C1, C2, COEF, PI /6.66666666666666667E-01, 00127 * 6.14926627446000736E-01,4.48288357353826359E-01, 00128 * 5.77350269189625765E-01,3.14159265358979324E+00/ 00129 DATA CONE / (1.0E0,0.0E0) / 00130 C***FIRST EXECUTABLE STATEMENT CBIRY 00131 IERR = 0 00132 NZ=0 00133 IF (ID.LT.0 .OR. ID.GT.1) IERR=1 00134 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 00135 IF (IERR.NE.0) RETURN 00136 AZ = CABS(Z) 00137 TOL = AMAX1(R1MACH(4),1.0E-18) 00138 FID = FLOAT(ID) 00139 IF (AZ.GT.1.0E0) GO TO 60 00140 C----------------------------------------------------------------------- 00141 C POWER SERIES FOR CABS(Z).LE.1. 00142 C----------------------------------------------------------------------- 00143 S1 = CONE 00144 S2 = CONE 00145 IF (AZ.LT.TOL) GO TO 110 00146 AA = AZ*AZ 00147 IF (AA.LT.TOL/AZ) GO TO 40 00148 TRM1 = CONE 00149 TRM2 = CONE 00150 ATRM = 1.0E0 00151 Z3 = Z*Z*Z 00152 AZ3 = AZ*AA 00153 AK = 2.0E0 + FID 00154 BK = 3.0E0 - FID - FID 00155 CK = 4.0E0 - FID 00156 DK = 3.0E0 + FID + FID 00157 D1 = AK*DK 00158 D2 = BK*CK 00159 AD = AMIN1(D1,D2) 00160 AK = 24.0E0 + 9.0E0*FID 00161 BK = 30.0E0 - 9.0E0*FID 00162 Z3R = REAL(Z3) 00163 Z3I = AIMAG(Z3) 00164 DO 30 K=1,25 00165 TRM1 = TRM1*CMPLX(Z3R/D1,Z3I/D1) 00166 S1 = S1 + TRM1 00167 TRM2 = TRM2*CMPLX(Z3R/D2,Z3I/D2) 00168 S2 = S2 + TRM2 00169 ATRM = ATRM*AZ3/AD 00170 D1 = D1 + AK 00171 D2 = D2 + BK 00172 AD = AMIN1(D1,D2) 00173 IF (ATRM.LT.TOL*AD) GO TO 40 00174 AK = AK + 18.0E0 00175 BK = BK + 18.0E0 00176 30 CONTINUE 00177 40 CONTINUE 00178 IF (ID.EQ.1) GO TO 50 00179 BI = S1*CMPLX(C1,0.0E0) + Z*S2*CMPLX(C2,0.0E0) 00180 IF (KODE.EQ.1) RETURN 00181 ZTA = Z*CSQRT(Z)*CMPLX(TTH,0.0E0) 00182 AA = REAL(ZTA) 00183 AA = -ABS(AA) 00184 BI = BI*CMPLX(EXP(AA),0.0E0) 00185 RETURN 00186 50 CONTINUE 00187 BI = S2*CMPLX(C2,0.0E0) 00188 IF (AZ.GT.TOL) BI = BI + Z*Z*S1*CMPLX(C1/(1.0E0+FID),0.0E0) 00189 IF (KODE.EQ.1) RETURN 00190 ZTA = Z*CSQRT(Z)*CMPLX(TTH,0.0E0) 00191 AA = REAL(ZTA) 00192 AA = -ABS(AA) 00193 BI = BI*CMPLX(EXP(AA),0.0E0) 00194 RETURN 00195 C----------------------------------------------------------------------- 00196 C CASE FOR CABS(Z).GT.1.0 00197 C----------------------------------------------------------------------- 00198 60 CONTINUE 00199 FNU = (1.0E0+FID)/3.0E0 00200 C----------------------------------------------------------------------- 00201 C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 00202 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. 00203 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 00204 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 00205 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 00206 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 00207 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 00208 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 00209 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. 00210 C----------------------------------------------------------------------- 00211 K1 = I1MACH(12) 00212 K2 = I1MACH(13) 00213 R1M5 = R1MACH(5) 00214 K = MIN0(IABS(K1),IABS(K2)) 00215 ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) 00216 K1 = I1MACH(11) - 1 00217 AA = R1M5*FLOAT(K1) 00218 DIG = AMIN1(AA,18.0E0) 00219 AA = AA*2.303E0 00220 ALIM = ELIM + AMAX1(-AA,-41.45E0) 00221 RL = 1.2E0*DIG + 3.0E0 00222 FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) 00223 C----------------------------------------------------------------------- 00224 C TEST FOR RANGE 00225 C----------------------------------------------------------------------- 00226 AA=0.5E0/TOL 00227 BB=FLOAT(I1MACH(9))*0.5E0 00228 AA=AMIN1(AA,BB) 00229 AA=AA**TTH 00230 IF (AZ.GT.AA) GO TO 190 00231 AA=SQRT(AA) 00232 IF (AZ.GT.AA) IERR=3 00233 CSQ=CSQRT(Z) 00234 ZTA=Z*CSQ*CMPLX(TTH,0.0E0) 00235 C----------------------------------------------------------------------- 00236 C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL 00237 C----------------------------------------------------------------------- 00238 SFAC = 1.0E0 00239 ZI = AIMAG(Z) 00240 ZR = REAL(Z) 00241 AK = AIMAG(ZTA) 00242 IF (ZR.GE.0.0E0) GO TO 70 00243 BK = REAL(ZTA) 00244 CK = -ABS(BK) 00245 ZTA = CMPLX(CK,AK) 00246 70 CONTINUE 00247 IF (ZI.EQ.0.0E0 .AND. ZR.LE.0.0E0) ZTA = CMPLX(0.0E0,AK) 00248 AA = REAL(ZTA) 00249 IF (KODE.EQ.2) GO TO 80 00250 C----------------------------------------------------------------------- 00251 C OVERFLOW TEST 00252 C----------------------------------------------------------------------- 00253 BB = ABS(AA) 00254 IF (BB.LT.ALIM) GO TO 80 00255 BB = BB + 0.25E0*ALOG(AZ) 00256 SFAC = TOL 00257 IF (BB.GT.ELIM) GO TO 170 00258 80 CONTINUE 00259 FMR = 0.0E0 00260 IF (AA.GE.0.0E0 .AND. ZR.GT.0.0E0) GO TO 90 00261 FMR = PI 00262 IF (ZI.LT.0.0E0) FMR = -PI 00263 ZTA = -ZTA 00264 90 CONTINUE 00265 C----------------------------------------------------------------------- 00266 C AA=FACTOR FOR ANALYTIC CONTINUATION OF I(FNU,ZTA) 00267 C KODE=2 RETURNS EXP(-ABS(XZTA))*I(FNU,ZTA) FROM CBINU 00268 C----------------------------------------------------------------------- 00269 CALL CBINU(ZTA, FNU, KODE, 1, CY, NZ, RL, FNUL, TOL, ELIM, ALIM) 00270 IF (NZ.LT.0) GO TO 180 00271 AA = FMR*FNU 00272 Z3 = CMPLX(SFAC,0.0E0) 00273 S1 = CY(1)*CMPLX(COS(AA),SIN(AA))*Z3 00274 FNU = (2.0E0-FID)/3.0E0 00275 CALL CBINU(ZTA, FNU, KODE, 2, CY, NZ, RL, FNUL, TOL, ELIM, ALIM) 00276 CY(1) = CY(1)*Z3 00277 CY(2) = CY(2)*Z3 00278 C----------------------------------------------------------------------- 00279 C BACKWARD RECUR ONE STEP FOR ORDERS -1/3 OR -2/3 00280 C----------------------------------------------------------------------- 00281 S2 = CY(1)*CMPLX(FNU+FNU,0.0E0)/ZTA + CY(2) 00282 AA = FMR*(FNU-1.0E0) 00283 S1 = (S1+S2*CMPLX(COS(AA),SIN(AA)))*CMPLX(COEF,0.0E0) 00284 IF (ID.EQ.1) GO TO 100 00285 S1 = CSQ*S1 00286 BI = S1*CMPLX(1.0E0/SFAC,0.0E0) 00287 RETURN 00288 100 CONTINUE 00289 S1 = Z*S1 00290 BI = S1*CMPLX(1.0E0/SFAC,0.0E0) 00291 RETURN 00292 110 CONTINUE 00293 AA = C1*(1.0E0-FID) + FID*C2 00294 BI = CMPLX(AA,0.0E0) 00295 RETURN 00296 170 CONTINUE 00297 NZ=0 00298 IERR=2 00299 RETURN 00300 180 CONTINUE 00301 IF(NZ.EQ.(-1)) GO TO 170 00302 NZ=0 00303 IERR=5 00304 RETURN 00305 190 CONTINUE 00306 IERR=4 00307 NZ=0 00308 RETURN 00309 END
1.7.1
