00001 SUBROUTINE ZBESH(ZR, ZI, FNU, KODE, M, N, CYR, CYI, NZ, IERR) 00002 C***BEGIN PROLOGUE ZBESH 00003 C***DATE WRITTEN 830501 (YYMMDD) 00004 C***REVISION DATE 890801 (YYMMDD) 00005 C***CATEGORY NO. B5K 00006 C***KEYWORDS H-BESSEL FUNCTIONS,BESSEL FUNCTIONS OF COMPLEX ARGUMENT, 00007 C BESSEL FUNCTIONS OF THIRD KIND,HANKEL FUNCTIONS 00008 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 00009 C***PURPOSE TO COMPUTE THE H-BESSEL FUNCTIONS OF A COMPLEX ARGUMENT 00010 C***DESCRIPTION 00011 C 00012 C ***A DOUBLE PRECISION ROUTINE*** 00013 C ON KODE=1, ZBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 00014 C HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1 00015 C OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX 00016 C Z.NE.CMPLX(0.0,0.0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI. 00017 C ON KODE=2, ZBESH RETURNS THE SCALED HANKEL FUNCTIONS 00018 C 00019 C CY(I)=EXP(-MM*Z*I)*H(M,FNU+J-1,Z) MM=3-2*M, I**2=-1. 00020 C 00021 C WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER AND 00022 C LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN THE 00023 C NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1). 00024 C 00025 C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION 00026 C ZR,ZI - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0), 00027 C -PT.LT.ARG(Z).LE.PI 00028 C FNU - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0D0 00029 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 00030 C KODE= 1 RETURNS 00031 C CY(J)=H(M,FNU+J-1,Z), J=1,...,N 00032 C = 2 RETURNS 00033 C CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) 00034 C J=1,...,N , I**2=-1 00035 C M - KIND OF HANKEL FUNCTION, M=1 OR 2 00036 C N - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1 00037 C 00038 C OUTPUT CYR,CYI ARE DOUBLE PRECISION 00039 C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS 00040 C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE 00041 C CY(J)=H(M,FNU+J-1,Z) OR 00042 C CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) J=1,...,N 00043 C DEPENDING ON KODE, I**2=-1. 00044 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, 00045 C NZ= 0 , NORMAL RETURN 00046 C NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE 00047 C TO UNDERFLOW, CY(J)=CMPLX(0.0D0,0.0D0) 00048 C J=1,...,NZ WHEN Y.GT.0.0 AND M=1 OR 00049 C Y.LT.0.0 AND M=2. FOR THE COMPLMENTARY 00050 C HALF PLANES, NZ STATES ONLY THE NUMBER 00051 C OF UNDERFLOWS. 00052 C IERR - ERROR FLAG 00053 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 00054 C IERR=1, INPUT ERROR - NO COMPUTATION 00055 C IERR=2, OVERFLOW - NO COMPUTATION, FNU TOO 00056 C LARGE OR CABS(Z) TOO SMALL OR BOTH 00057 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE 00058 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT 00059 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE 00060 C ACCURACY 00061 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- 00062 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- 00063 C CANCE BY ARGUMENT REDUCTION 00064 C IERR=5, ERROR - NO COMPUTATION, 00065 C ALGORITHM TERMINATION CONDITION NOT MET 00066 C 00067 C***LONG DESCRIPTION 00068 C 00069 C THE COMPUTATION IS CARRIED OUT BY THE RELATION 00070 C 00071 C H(M,FNU,Z)=(1/MP)*EXP(-MP*FNU)*K(FNU,Z*EXP(-MP)) 00072 C MP=MM*HPI*I, MM=3-2*M, HPI=PI/2, I**2=-1 00073 C 00074 C FOR M=1 OR 2 WHERE THE K BESSEL FUNCTION IS COMPUTED FOR THE 00075 C RIGHT HALF PLANE RE(Z).GE.0.0. THE K FUNCTION IS CONTINUED 00076 C TO THE LEFT HALF PLANE BY THE RELATION 00077 C 00078 C K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z) 00079 C MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1 00080 C 00081 C WHERE I(FNU,Z) IS THE I BESSEL FUNCTION. 00082 C 00083 C EXPONENTIAL DECAY OF H(M,FNU,Z) OCCURS IN THE UPPER HALF Z 00084 C PLANE FOR M=1 AND THE LOWER HALF Z PLANE FOR M=2. EXPONENTIAL 00085 C GROWTH OCCURS IN THE COMPLEMENTARY HALF PLANES. SCALING 00086 C BY EXP(-MM*Z*I) REMOVES THE EXPONENTIAL BEHAVIOR IN THE 00087 C WHOLE Z PLANE FOR Z TO INFINITY. 00088 C 00089 C FOR NEGATIVE ORDERS,THE FORMULAE 00090 C 00091 C H(1,-FNU,Z) = H(1,FNU,Z)*CEXP( PI*FNU*I) 00092 C H(2,-FNU,Z) = H(2,FNU,Z)*CEXP(-PI*FNU*I) 00093 C I**2=-1 00094 C 00095 C CAN BE USED. 00096 C 00097 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 00098 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS 00099 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. 00100 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN 00101 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG 00102 C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS 00103 C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. 00104 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS 00105 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS 00106 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE 00107 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS 00108 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 00109 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION 00110 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION 00111 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN 00112 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT 00113 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS 00114 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. 00115 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. 00116 C 00117 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 00118 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 00119 C ROUNDOFF,1.0D-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 00120 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 00121 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 00122 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 00123 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 00124 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 00125 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 00126 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 00127 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 00128 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 00129 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 00130 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 00131 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 00132 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 00133 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 00134 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 00135 C OR -PI/2+P. 00136 C 00137 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 00138 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 00139 C COMMERCE, 1955. 00140 C 00141 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00142 C BY D. E. AMOS, SAND83-0083, MAY, 1983. 00143 C 00144 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 00145 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 00146 C 00147 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 00148 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 00149 C 1018, MAY, 1985 00150 C 00151 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 00152 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 00153 C MATH. SOFTWARE, 1986 00154 C 00155 C***ROUTINES CALLED ZACON,ZBKNU,ZBUNK,ZUOIK,XZABS,I1MACH,D1MACH 00156 C***END PROLOGUE ZBESH 00157 C 00158 C COMPLEX CY,Z,ZN,ZT,CSGN 00159 DOUBLE PRECISION AA, ALIM, ALN, ARG, AZ, CYI, CYR, DIG, ELIM, 00160 * FMM, FN, FNU, FNUL, HPI, RHPI, RL, R1M5, SGN, STR, TOL, UFL, ZI, 00161 * ZNI, ZNR, ZR, ZTI, D1MACH, XZABS, BB, ASCLE, RTOL, ATOL, STI, 00162 * CSGNR, CSGNI 00163 INTEGER I, IERR, INU, INUH, IR, K, KODE, K1, K2, M, 00164 * MM, MR, N, NN, NUF, NW, NZ, I1MACH 00165 DIMENSION CYR(N), CYI(N) 00166 C 00167 DATA HPI /1.57079632679489662D0/ 00168 C 00169 C***FIRST EXECUTABLE STATEMENT ZBESH 00170 IERR = 0 00171 NZ=0 00172 IF (ZR.EQ.0.0D0 .AND. ZI.EQ.0.0D0) IERR=1 00173 IF (FNU.LT.0.0D0) IERR=1 00174 IF (M.LT.1 .OR. M.GT.2) IERR=1 00175 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 00176 IF (N.LT.1) IERR=1 00177 IF (IERR.NE.0) RETURN 00178 NN = N 00179 C----------------------------------------------------------------------- 00180 C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 00181 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. 00182 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 00183 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 00184 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 00185 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 00186 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 00187 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 00188 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU 00189 C----------------------------------------------------------------------- 00190 TOL = DMAX1(D1MACH(4),1.0D-18) 00191 K1 = I1MACH(15) 00192 K2 = I1MACH(16) 00193 R1M5 = D1MACH(5) 00194 K = MIN0(IABS(K1),IABS(K2)) 00195 ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) 00196 K1 = I1MACH(14) - 1 00197 AA = R1M5*DBLE(FLOAT(K1)) 00198 DIG = DMIN1(AA,18.0D0) 00199 AA = AA*2.303D0 00200 ALIM = ELIM + DMAX1(-AA,-41.45D0) 00201 FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) 00202 RL = 1.2D0*DIG + 3.0D0 00203 FN = FNU + DBLE(FLOAT(NN-1)) 00204 MM = 3 - M - M 00205 FMM = DBLE(FLOAT(MM)) 00206 ZNR = FMM*ZI 00207 ZNI = -FMM*ZR 00208 C----------------------------------------------------------------------- 00209 C TEST FOR PROPER RANGE 00210 C----------------------------------------------------------------------- 00211 AZ = XZABS(ZR,ZI) 00212 AA = 0.5D0/TOL 00213 BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 00214 AA = DMIN1(AA,BB) 00215 IF (AZ.GT.AA) GO TO 260 00216 IF (FN.GT.AA) GO TO 260 00217 AA = DSQRT(AA) 00218 IF (AZ.GT.AA) IERR=3 00219 IF (FN.GT.AA) IERR=3 00220 C----------------------------------------------------------------------- 00221 C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE 00222 C----------------------------------------------------------------------- 00223 UFL = D1MACH(1)*1.0D+3 00224 IF (AZ.LT.UFL) GO TO 230 00225 IF (FNU.GT.FNUL) GO TO 90 00226 IF (FN.LE.1.0D0) GO TO 70 00227 IF (FN.GT.2.0D0) GO TO 60 00228 IF (AZ.GT.TOL) GO TO 70 00229 ARG = 0.5D0*AZ 00230 ALN = -FN*DLOG(ARG) 00231 IF (ALN.GT.ELIM) GO TO 230 00232 GO TO 70 00233 60 CONTINUE 00234 CALL ZUOIK(ZNR, ZNI, FNU, KODE, 2, NN, CYR, CYI, NUF, TOL, ELIM, 00235 * ALIM) 00236 IF (NUF.LT.0) GO TO 230 00237 NZ = NZ + NUF 00238 NN = NN - NUF 00239 C----------------------------------------------------------------------- 00240 C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK 00241 C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I 00242 C----------------------------------------------------------------------- 00243 IF (NN.EQ.0) GO TO 140 00244 70 CONTINUE 00245 IF ((ZNR.LT.0.0D0) .OR. (ZNR.EQ.0.0D0 .AND. ZNI.LT.0.0D0 .AND. 00246 * M.EQ.2)) GO TO 80 00247 C----------------------------------------------------------------------- 00248 C RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR. 00249 C YN.GE.0. .OR. M=1) 00250 C----------------------------------------------------------------------- 00251 CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NZ, TOL, ELIM, ALIM) 00252 GO TO 110 00253 C----------------------------------------------------------------------- 00254 C LEFT HALF PLANE COMPUTATION 00255 C----------------------------------------------------------------------- 00256 80 CONTINUE 00257 MR = -MM 00258 CALL ZACON(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, RL, FNUL, 00259 * TOL, ELIM, ALIM) 00260 IF (NW.LT.0) GO TO 240 00261 NZ=NW 00262 GO TO 110 00263 90 CONTINUE 00264 C----------------------------------------------------------------------- 00265 C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL 00266 C----------------------------------------------------------------------- 00267 MR = 0 00268 IF ((ZNR.GE.0.0D0) .AND. (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0 .OR. 00269 * M.NE.2)) GO TO 100 00270 MR = -MM 00271 IF (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0) GO TO 100 00272 ZNR = -ZNR 00273 ZNI = -ZNI 00274 100 CONTINUE 00275 CALL ZBUNK(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, TOL, ELIM, 00276 * ALIM) 00277 IF (NW.LT.0) GO TO 240 00278 NZ = NZ + NW 00279 110 CONTINUE 00280 C----------------------------------------------------------------------- 00281 C H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT) 00282 C 00283 C ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2 00284 C----------------------------------------------------------------------- 00285 SGN = DSIGN(HPI,-FMM) 00286 C----------------------------------------------------------------------- 00287 C CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE 00288 C WHEN FNU IS LARGE 00289 C----------------------------------------------------------------------- 00290 INU = INT(SNGL(FNU)) 00291 INUH = INU/2 00292 IR = INU - 2*INUH 00293 ARG = (FNU-DBLE(FLOAT(INU-IR)))*SGN 00294 RHPI = 1.0D0/SGN 00295 C ZNI = RHPI*DCOS(ARG) 00296 C ZNR = -RHPI*DSIN(ARG) 00297 CSGNI = RHPI*DCOS(ARG) 00298 CSGNR = -RHPI*DSIN(ARG) 00299 IF (MOD(INUH,2).EQ.0) GO TO 120 00300 C ZNR = -ZNR 00301 C ZNI = -ZNI 00302 CSGNR = -CSGNR 00303 CSGNI = -CSGNI 00304 120 CONTINUE 00305 ZTI = -FMM 00306 RTOL = 1.0D0/TOL 00307 ASCLE = UFL*RTOL 00308 DO 130 I=1,NN 00309 C STR = CYR(I)*ZNR - CYI(I)*ZNI 00310 C CYI(I) = CYR(I)*ZNI + CYI(I)*ZNR 00311 C CYR(I) = STR 00312 C STR = -ZNI*ZTI 00313 C ZNI = ZNR*ZTI 00314 C ZNR = STR 00315 AA = CYR(I) 00316 BB = CYI(I) 00317 ATOL = 1.0D0 00318 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 135 00319 AA = AA*RTOL 00320 BB = BB*RTOL 00321 ATOL = TOL 00322 135 CONTINUE 00323 STR = AA*CSGNR - BB*CSGNI 00324 STI = AA*CSGNI + BB*CSGNR 00325 CYR(I) = STR*ATOL 00326 CYI(I) = STI*ATOL 00327 STR = -CSGNI*ZTI 00328 CSGNI = CSGNR*ZTI 00329 CSGNR = STR 00330 130 CONTINUE 00331 RETURN 00332 140 CONTINUE 00333 IF (ZNR.LT.0.0D0) GO TO 230 00334 RETURN 00335 230 CONTINUE 00336 NZ=0 00337 IERR=2 00338 RETURN 00339 240 CONTINUE 00340 IF(NW.EQ.(-1)) GO TO 230 00341 NZ=0 00342 IERR=5 00343 RETURN 00344 260 CONTINUE 00345 NZ=0 00346 IERR=4 00347 RETURN 00348 END
1.7.1
