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