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