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,XZABS,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, XZABS, 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 = XZABS(ZR,ZI) 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 35 CONTINUE 224 UFL = D1MACH(1)*1.0D+3 225 IF (AZ.LT.UFL) GO TO 230 226 IF (FNU.GT.FNUL) GO TO 90 227 IF (FN.LE.1.0D0) GO TO 70 228 IF (FN.GT.2.0D0) GO TO 60 229 IF (AZ.GT.TOL) GO TO 70 230 ARG = 0.5D0*AZ 231 ALN = -FN*DLOG(ARG) 232 IF (ALN.GT.ELIM) GO TO 230 233 GO TO 70 234 60 CONTINUE 235 CALL ZUOIK(ZNR, ZNI, FNU, KODE, 2, NN, CYR, CYI, NUF, TOL, ELIM, 236 * ALIM) 237 IF (NUF.LT.0) GO TO 230 238 NZ = NZ + NUF 239 NN = NN - NUF 240C----------------------------------------------------------------------- 241C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK 242C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I 243C----------------------------------------------------------------------- 244 IF (NN.EQ.0) GO TO 140 245 70 CONTINUE 246 IF ((ZNR.LT.0.0D0) .OR. (ZNR.EQ.0.0D0 .AND. ZNI.LT.0.0D0 .AND. 247 * M.EQ.2)) GO TO 80 248C----------------------------------------------------------------------- 249C RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR. 250C YN.GE.0. .OR. M=1) 251C----------------------------------------------------------------------- 252 CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NZ, TOL, ELIM, ALIM) 253 GO TO 110 254C----------------------------------------------------------------------- 255C LEFT HALF PLANE COMPUTATION 256C----------------------------------------------------------------------- 257 80 CONTINUE 258 MR = -MM 259 CALL ZACON(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, RL, FNUL, 260 * TOL, ELIM, ALIM) 261 IF (NW.LT.0) GO TO 240 262 NZ=NW 263 GO TO 110 264 90 CONTINUE 265C----------------------------------------------------------------------- 266C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL 267C----------------------------------------------------------------------- 268 MR = 0 269 IF ((ZNR.GE.0.0D0) .AND. (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0 .OR. 270 * M.NE.2)) GO TO 100 271 MR = -MM 272 IF (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0) GO TO 100 273 ZNR = -ZNR 274 ZNI = -ZNI 275 100 CONTINUE 276 CALL ZBUNK(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, TOL, ELIM, 277 * ALIM) 278 IF (NW.LT.0) GO TO 240 279 NZ = NZ + NW 280 110 CONTINUE 281C----------------------------------------------------------------------- 282C H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT) 283C 284C ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2 285C----------------------------------------------------------------------- 286 SGN = DSIGN(HPI,-FMM) 287C----------------------------------------------------------------------- 288C CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE 289C WHEN FNU IS LARGE 290C----------------------------------------------------------------------- 291 INU = INT(SNGL(FNU)) 292 INUH = INU/2 293 IR = INU - 2*INUH 294 ARG = (FNU-DBLE(FLOAT(INU-IR)))*SGN 295 RHPI = 1.0D0/SGN 296C ZNI = RHPI*DCOS(ARG) 297C ZNR = -RHPI*DSIN(ARG) 298 CSGNI = RHPI*DCOS(ARG) 299 CSGNR = -RHPI*DSIN(ARG) 300 IF (MOD(INUH,2).EQ.0) GO TO 120 301C ZNR = -ZNR 302C ZNI = -ZNI 303 CSGNR = -CSGNR 304 CSGNI = -CSGNI 305 120 CONTINUE 306 ZTI = -FMM 307 RTOL = 1.0D0/TOL 308 ASCLE = UFL*RTOL 309 DO 130 I=1,NN 310C STR = CYR(I)*ZNR - CYI(I)*ZNI 311C CYI(I) = CYR(I)*ZNI + CYI(I)*ZNR 312C CYR(I) = STR 313C STR = -ZNI*ZTI 314C ZNI = ZNR*ZTI 315C ZNR = STR 316 AA = CYR(I) 317 BB = CYI(I) 318 ATOL = 1.0D0 319 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 135 320 AA = AA*RTOL 321 BB = BB*RTOL 322 ATOL = TOL 323 135 CONTINUE 324 STR = AA*CSGNR - BB*CSGNI 325 STI = AA*CSGNI + BB*CSGNR 326 CYR(I) = STR*ATOL 327 CYI(I) = STI*ATOL 328 STR = -CSGNI*ZTI 329 CSGNI = CSGNR*ZTI 330 CSGNR = STR 331 130 CONTINUE 332 RETURN 333 140 CONTINUE 334 IF (ZNR.LT.0.0D0) GO TO 230 335 RETURN 336 230 CONTINUE 337 NZ=0 338 IERR=2 339 RETURN 340 240 CONTINUE 341 IF(NW.EQ.(-1)) GO TO 230 342 NZ=0 343 IERR=5 344 RETURN 345 260 CONTINUE 346 IERR=4 347 GO TO 35 348 END 349