1 SUBROUTINE ZBESH(ZR, ZI, FNU, KODE, M, N, CYR, CYI, NZ, IERR) 2C***BEGIN PROLOGUE ZBESH 3C***DATE WRITTEN 830501 (YYMMDD) 4C***REVISION DATE 890801 (YYMMDD) 5C***CATEGORY NO. B5K 6C***KEYWORDS H-BESSEL FUNCTIONS,BESSEL FUNCTIONS OF COMPLEX ARGUMENT, 7C BESSEL FUNCTIONS OF THIRD KIND,HANKEL FUNCTIONS 8C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 9C***PURPOSE TO COMPUTE THE H-BESSEL FUNCTIONS OF A COMPLEX ARGUMENT 10C***DESCRIPTION 11C 12C ***A DOUBLE PRECISION ROUTINE*** 13C ON KODE=1, ZBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 14C HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1 15C OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX 16C Z.NE.CMPLX(0.0,0.0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI. 17C ON KODE=2, ZBESH RETURNS THE SCALED HANKEL FUNCTIONS 18C 19C CY(I)=EXP(-MM*Z*I)*H(M,FNU+J-1,Z) MM=3-2*M, I**2=-1. 20C 21C WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER AND 22C LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN THE 23C NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1). 24C 25C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION 26C ZR,ZI - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0), 27C -PT.LT.ARG(Z).LE.PI 28C FNU - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0D0 29C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 30C KODE= 1 RETURNS 31C CY(J)=H(M,FNU+J-1,Z), J=1,...,N 32C = 2 RETURNS 33C CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) 34C J=1,...,N , I**2=-1 35C M - KIND OF HANKEL FUNCTION, M=1 OR 2 36C N - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1 37C 38C OUTPUT CYR,CYI ARE DOUBLE PRECISION 39C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS 40C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE 41C CY(J)=H(M,FNU+J-1,Z) OR 42C CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) J=1,...,N 43C DEPENDING ON KODE, I**2=-1. 44C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, 45C NZ= 0 , NORMAL RETURN 46C NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE 47C TO UNDERFLOW, CY(J)=CMPLX(0.0D0,0.0D0) 48C J=1,...,NZ WHEN Y.GT.0.0 AND M=1 OR 49C Y.LT.0.0 AND M=2. FOR THE COMPLMENTARY 50C HALF PLANES, NZ STATES ONLY THE NUMBER 51C OF UNDERFLOWS. 52C IERR - ERROR FLAG 53C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 54C IERR=1, INPUT ERROR - NO COMPUTATION 55C IERR=2, OVERFLOW - NO COMPUTATION, FNU TOO 56C LARGE OR CABS(Z) TOO SMALL OR BOTH 57C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE 58C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT 59C REDUCTION PRODUCE LESS THAN HALF OF MACHINE 60C ACCURACY 61C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- 62C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- 63C CANCE BY ARGUMENT REDUCTION 64C IERR=5, ERROR - NO COMPUTATION, 65C ALGORITHM TERMINATION CONDITION NOT MET 66C 67C***LONG DESCRIPTION 68C 69C THE COMPUTATION IS CARRIED OUT BY THE RELATION 70C 71C H(M,FNU,Z)=(1/MP)*EXP(-MP*FNU)*K(FNU,Z*EXP(-MP)) 72C MP=MM*HPI*I, MM=3-2*M, HPI=PI/2, I**2=-1 73C 74C FOR M=1 OR 2 WHERE THE K BESSEL FUNCTION IS COMPUTED FOR THE 75C RIGHT HALF PLANE RE(Z).GE.0.0. THE K FUNCTION IS CONTINUED 76C TO THE LEFT HALF PLANE BY THE RELATION 77C 78C K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z) 79C MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1 80C 81C WHERE I(FNU,Z) IS THE I BESSEL FUNCTION. 82C 83C EXPONENTIAL DECAY OF H(M,FNU,Z) OCCURS IN THE UPPER HALF Z 84C PLANE FOR M=1 AND THE LOWER HALF Z PLANE FOR M=2. EXPONENTIAL 85C GROWTH OCCURS IN THE COMPLEMENTARY HALF PLANES. SCALING 86C BY EXP(-MM*Z*I) REMOVES THE EXPONENTIAL BEHAVIOR IN THE 87C WHOLE Z PLANE FOR Z TO INFINITY. 88C 89C FOR NEGATIVE ORDERS,THE FORMULAE 90C 91C H(1,-FNU,Z) = H(1,FNU,Z)*CEXP( PI*FNU*I) 92C H(2,-FNU,Z) = H(2,FNU,Z)*CEXP(-PI*FNU*I) 93C I**2=-1 94C 95C CAN BE USED. 96C 97C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 98C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS 99C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. 100C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN 101C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG 102C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS 103C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. 104C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS 105C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS 106C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE 107C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS 108C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 109C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION 110C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION 111C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN 112C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT 113C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS 114C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. 115C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. 116C 117C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 118C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 119C ROUNDOFF,1.0D-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 120C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 121C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 122C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 123C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 124C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 125C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 126C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 127C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 128C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 129C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 130C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 131C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 132C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 133C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 134C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 135C OR -PI/2+P. 136C 137C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 138C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 139C COMMERCE, 1955. 140C 141C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 142C BY D. E. AMOS, SAND83-0083, MAY, 1983. 143C 144C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 145C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 146C 147C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 148C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 149C 1018, MAY, 1985 150C 151C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 152C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 153C MATH. SOFTWARE, 1986 154C 155C***ROUTINES CALLED ZACON,ZBKNU,ZBUNK,ZUOIK,ZABS,I1MACH,D1MACH 156C***END PROLOGUE ZBESH 157C 158C COMPLEX CY,Z,ZN,ZT,CSGN 159 DOUBLE PRECISION AA, ALIM, ALN, ARG, AZ, CYI, CYR, DIG, ELIM, 160 * FMM, FN, FNU, FNUL, HPI, RHPI, RL, R1M5, SGN, STR, TOL, UFL, ZI, 161 * ZNI, ZNR, ZR, ZTI, D1MACH, ZABS, BB, ASCLE, RTOL, ATOL, STI, 162 * CSGNR, CSGNI 163 INTEGER I, IERR, INU, INUH, IR, K, KODE, K1, K2, M, 164 * MM, MR, N, NN, NUF, NW, NZ, I1MACH 165 DIMENSION CYR(N), CYI(N) 166C 167 DATA HPI /1.57079632679489662D0/ 168C 169C***FIRST EXECUTABLE STATEMENT ZBESH 170 IERR = 0 171 NZ=0 172 IF (ZR.EQ.0.0D0 .AND. ZI.EQ.0.0D0) IERR=1 173 IF (FNU.LT.0.0D0) IERR=1 174 IF (M.LT.1 .OR. M.GT.2) IERR=1 175 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 176 IF (N.LT.1) IERR=1 177 IF (IERR.NE.0) RETURN 178 NN = N 179C----------------------------------------------------------------------- 180C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 181C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. 182C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 183C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 184C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 185C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 186C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 187C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 188C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU 189C----------------------------------------------------------------------- 190 TOL = DMAX1(D1MACH(4),1.0D-18) 191 K1 = I1MACH(15) 192 K2 = I1MACH(16) 193 R1M5 = D1MACH(5) 194 K = MIN0(IABS(K1),IABS(K2)) 195 ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) 196 K1 = I1MACH(14) - 1 197 AA = R1M5*DBLE(FLOAT(K1)) 198 DIG = DMIN1(AA,18.0D0) 199 AA = AA*2.303D0 200 ALIM = ELIM + DMAX1(-AA,-41.45D0) 201 FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) 202 RL = 1.2D0*DIG + 3.0D0 203 FN = FNU + DBLE(FLOAT(NN-1)) 204 MM = 3 - M - M 205 FMM = DBLE(FLOAT(MM)) 206 ZNR = FMM*ZI 207 ZNI = -FMM*ZR 208C----------------------------------------------------------------------- 209C TEST FOR PROPER RANGE 210C----------------------------------------------------------------------- 211 AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 212 AA = 0.5D0/TOL 213 BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 214 AA = DMIN1(AA,BB) 215 IF (AZ.GT.AA) GO TO 260 216 IF (FN.GT.AA) GO TO 260 217 AA = DSQRT(AA) 218 IF (AZ.GT.AA) IERR=3 219 IF (FN.GT.AA) IERR=3 220C----------------------------------------------------------------------- 221C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE 222C----------------------------------------------------------------------- 223 UFL = D1MACH(1)*1.0D+3 224 IF (AZ.LT.UFL) GO TO 230 225 IF (FNU.GT.FNUL) GO TO 90 226 IF (FN.LE.1.0D0) GO TO 70 227 IF (FN.GT.2.0D0) GO TO 60 228 IF (AZ.GT.TOL) GO TO 70 229 ARG = 0.5D0*AZ 230 ALN = -FN*DLOG(ARG) 231 IF (ALN.GT.ELIM) GO TO 230 232 GO TO 70 233 60 CONTINUE 234 CALL ZUOIK(ZNR, ZNI, FNU, KODE, 2, NN, CYR, CYI, NUF, TOL, ELIM, 235 * ALIM) 236 IF (NUF.LT.0) GO TO 230 237 NZ = NZ + NUF 238 NN = NN - NUF 239C----------------------------------------------------------------------- 240C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK 241C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I 242C----------------------------------------------------------------------- 243 IF (NN.EQ.0) GO TO 140 244 70 CONTINUE 245 IF ((ZNR.LT.0.0D0) .OR. (ZNR.EQ.0.0D0 .AND. ZNI.LT.0.0D0 .AND. 246 * M.EQ.2)) GO TO 80 247C----------------------------------------------------------------------- 248C RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR. 249C YN.GE.0. .OR. M=1) 250C----------------------------------------------------------------------- 251 CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NZ, TOL, ELIM, ALIM) 252 GO TO 110 253C----------------------------------------------------------------------- 254C LEFT HALF PLANE COMPUTATION 255C----------------------------------------------------------------------- 256 80 CONTINUE 257 MR = -MM 258 CALL ZACON(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, RL, FNUL, 259 * TOL, ELIM, ALIM) 260 IF (NW.LT.0) GO TO 240 261 NZ=NW 262 GO TO 110 263 90 CONTINUE 264C----------------------------------------------------------------------- 265C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL 266C----------------------------------------------------------------------- 267 MR = 0 268 IF ((ZNR.GE.0.0D0) .AND. (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0 .OR. 269 * M.NE.2)) GO TO 100 270 MR = -MM 271 IF (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0) GO TO 100 272 ZNR = -ZNR 273 ZNI = -ZNI 274 100 CONTINUE 275 CALL ZBUNK(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, TOL, ELIM, 276 * ALIM) 277 IF (NW.LT.0) GO TO 240 278 NZ = NZ + NW 279 110 CONTINUE 280C----------------------------------------------------------------------- 281C H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT) 282C 283C ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2 284C----------------------------------------------------------------------- 285 SGN = DSIGN(HPI,-FMM) 286C----------------------------------------------------------------------- 287C CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE 288C WHEN FNU IS LARGE 289C----------------------------------------------------------------------- 290 INU = INT(SNGL(FNU)) 291 INUH = INU/2 292 IR = INU - 2*INUH 293 ARG = (FNU-DBLE(FLOAT(INU-IR)))*SGN 294 RHPI = 1.0D0/SGN 295C ZNI = RHPI*DCOS(ARG) 296C ZNR = -RHPI*DSIN(ARG) 297 CSGNI = RHPI*DCOS(ARG) 298 CSGNR = -RHPI*DSIN(ARG) 299 IF (MOD(INUH,2).EQ.0) GO TO 120 300C ZNR = -ZNR 301C ZNI = -ZNI 302 CSGNR = -CSGNR 303 CSGNI = -CSGNI 304 120 CONTINUE 305 ZTI = -FMM 306 RTOL = 1.0D0/TOL 307 ASCLE = UFL*RTOL 308 DO 130 I=1,NN 309C STR = CYR(I)*ZNR - CYI(I)*ZNI 310C CYI(I) = CYR(I)*ZNI + CYI(I)*ZNR 311C CYR(I) = STR 312C STR = -ZNI*ZTI 313C ZNI = ZNR*ZTI 314C ZNR = STR 315 AA = CYR(I) 316 BB = CYI(I) 317 ATOL = 1.0D0 318 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 135 319 AA = AA*RTOL 320 BB = BB*RTOL 321 ATOL = TOL 322 135 CONTINUE 323 STR = AA*CSGNR - BB*CSGNI 324 STI = AA*CSGNI + BB*CSGNR 325 CYR(I) = STR*ATOL 326 CYI(I) = STI*ATOL 327 STR = -CSGNI*ZTI 328 CSGNI = CSGNR*ZTI 329 CSGNR = STR 330 130 CONTINUE 331 RETURN 332 140 CONTINUE 333 IF (ZNR.LT.0.0D0) GO TO 230 334 RETURN 335 230 CONTINUE 336 NZ=0 337 IERR=2 338 RETURN 339 240 CONTINUE 340 IF(NW.EQ.(-1)) GO TO 230 341 NZ=0 342 IERR=5 343 RETURN 344 260 CONTINUE 345 NZ=0 346 IERR=4 347 RETURN 348 END 349