github.com/gopherd/gonum@v0.0.4/mathext/internal/amos/amoslib/zbesh.f (about) 1 SUBROUTINE ZBESH(ZR, ZI, FNU, KODE, M, N, CYR, CYI, NZ, IERR) 2 C***BEGIN PROLOGUE ZBESH 3 C***DATE WRITTEN 830501 (YYMMDD) 4 C***REVISION DATE 890801 (YYMMDD) 5 C***CATEGORY NO. B5K 6 C***KEYWORDS H-BESSEL FUNCTIONS,BESSEL FUNCTIONS OF COMPLEX ARGUMENT, 7 C BESSEL FUNCTIONS OF THIRD KIND,HANKEL FUNCTIONS 8 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 9 C***PURPOSE TO COMPUTE THE H-BESSEL FUNCTIONS OF A COMPLEX ARGUMENT 10 C***DESCRIPTION 11 C 12 C ***A DOUBLE PRECISION ROUTINE*** 13 C ON KODE=1, ZBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 14 C HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1 15 C OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX 16 C Z.NE.CMPLX(0.0,0.0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI. 17 C ON KODE=2, ZBESH RETURNS THE SCALED HANKEL FUNCTIONS 18 C 19 C CY(I)=EXP(-MM*Z*I)*H(M,FNU+J-1,Z) MM=3-2*M, I**2=-1. 20 C 21 C WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER AND 22 C LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN THE 23 C NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1). 24 C 25 C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION 26 C ZR,ZI - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0), 27 C -PT.LT.ARG(Z).LE.PI 28 C FNU - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0D0 29 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 30 C KODE= 1 RETURNS 31 C CY(J)=H(M,FNU+J-1,Z), J=1,...,N 32 C = 2 RETURNS 33 C CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) 34 C J=1,...,N , I**2=-1 35 C M - KIND OF HANKEL FUNCTION, M=1 OR 2 36 C N - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1 37 C 38 C OUTPUT CYR,CYI ARE DOUBLE PRECISION 39 C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS 40 C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE 41 C CY(J)=H(M,FNU+J-1,Z) OR 42 C CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) J=1,...,N 43 C DEPENDING ON KODE, I**2=-1. 44 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, 45 C NZ= 0 , NORMAL RETURN 46 C NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE 47 C TO UNDERFLOW, CY(J)=CMPLX(0.0D0,0.0D0) 48 C J=1,...,NZ WHEN Y.GT.0.0 AND M=1 OR 49 C Y.LT.0.0 AND M=2. FOR THE COMPLMENTARY 50 C HALF PLANES, NZ STATES ONLY THE NUMBER 51 C OF UNDERFLOWS. 52 C IERR - ERROR FLAG 53 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 54 C IERR=1, INPUT ERROR - NO COMPUTATION 55 C IERR=2, OVERFLOW - NO COMPUTATION, FNU TOO 56 C LARGE OR CABS(Z) TOO SMALL OR BOTH 57 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE 58 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT 59 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE 60 C ACCURACY 61 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- 62 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- 63 C CANCE BY ARGUMENT REDUCTION 64 C IERR=5, ERROR - NO COMPUTATION, 65 C ALGORITHM TERMINATION CONDITION NOT MET 66 C 67 C***LONG DESCRIPTION 68 C 69 C THE COMPUTATION IS CARRIED OUT BY THE RELATION 70 C 71 C H(M,FNU,Z)=(1/MP)*EXP(-MP*FNU)*K(FNU,Z*EXP(-MP)) 72 C MP=MM*HPI*I, MM=3-2*M, HPI=PI/2, I**2=-1 73 C 74 C FOR M=1 OR 2 WHERE THE K BESSEL FUNCTION IS COMPUTED FOR THE 75 C RIGHT HALF PLANE RE(Z).GE.0.0. THE K FUNCTION IS CONTINUED 76 C TO THE LEFT HALF PLANE BY THE RELATION 77 C 78 C K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z) 79 C MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1 80 C 81 C WHERE I(FNU,Z) IS THE I BESSEL FUNCTION. 82 C 83 C EXPONENTIAL DECAY OF H(M,FNU,Z) OCCURS IN THE UPPER HALF Z 84 C PLANE FOR M=1 AND THE LOWER HALF Z PLANE FOR M=2. EXPONENTIAL 85 C GROWTH OCCURS IN THE COMPLEMENTARY HALF PLANES. SCALING 86 C BY EXP(-MM*Z*I) REMOVES THE EXPONENTIAL BEHAVIOR IN THE 87 C WHOLE Z PLANE FOR Z TO INFINITY. 88 C 89 C FOR NEGATIVE ORDERS,THE FORMULAE 90 C 91 C H(1,-FNU,Z) = H(1,FNU,Z)*CEXP( PI*FNU*I) 92 C H(2,-FNU,Z) = H(2,FNU,Z)*CEXP(-PI*FNU*I) 93 C I**2=-1 94 C 95 C CAN BE USED. 96 C 97 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 98 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS 99 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. 100 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN 101 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG 102 C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS 103 C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. 104 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS 105 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS 106 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE 107 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS 108 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 109 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION 110 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION 111 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN 112 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT 113 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS 114 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. 115 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. 116 C 117 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 118 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 119 C ROUNDOFF,1.0D-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 120 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 121 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 122 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 123 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 124 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 125 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 126 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 127 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 128 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 129 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 130 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 131 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 132 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 133 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 134 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 135 C OR -PI/2+P. 136 C 137 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 138 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 139 C COMMERCE, 1955. 140 C 141 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 142 C BY D. E. AMOS, SAND83-0083, MAY, 1983. 143 C 144 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 145 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 146 C 147 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 148 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 149 C 1018, MAY, 1985 150 C 151 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 152 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 153 C MATH. SOFTWARE, 1986 154 C 155 C***ROUTINES CALLED ZACON,ZBKNU,ZBUNK,ZUOIK,ZABS,I1MACH,D1MACH 156 C***END PROLOGUE ZBESH 157 C 158 C 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) 166 C 167 DATA HPI /1.57079632679489662D0/ 168 C 169 C***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 179 C----------------------------------------------------------------------- 180 C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 181 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. 182 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 183 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 184 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 185 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 186 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 187 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 188 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU 189 C----------------------------------------------------------------------- 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 208 C----------------------------------------------------------------------- 209 C TEST FOR PROPER RANGE 210 C----------------------------------------------------------------------- 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 220 C----------------------------------------------------------------------- 221 C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE 222 C----------------------------------------------------------------------- 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 239 C----------------------------------------------------------------------- 240 C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK 241 C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I 242 C----------------------------------------------------------------------- 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 247 C----------------------------------------------------------------------- 248 C RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR. 249 C YN.GE.0. .OR. M=1) 250 C----------------------------------------------------------------------- 251 CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NZ, TOL, ELIM, ALIM) 252 GO TO 110 253 C----------------------------------------------------------------------- 254 C LEFT HALF PLANE COMPUTATION 255 C----------------------------------------------------------------------- 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 264 C----------------------------------------------------------------------- 265 C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL 266 C----------------------------------------------------------------------- 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 280 C----------------------------------------------------------------------- 281 C H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT) 282 C 283 C ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2 284 C----------------------------------------------------------------------- 285 SGN = DSIGN(HPI,-FMM) 286 C----------------------------------------------------------------------- 287 C CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE 288 C WHEN FNU IS LARGE 289 C----------------------------------------------------------------------- 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 295 C ZNI = RHPI*DCOS(ARG) 296 C ZNR = -RHPI*DSIN(ARG) 297 CSGNI = RHPI*DCOS(ARG) 298 CSGNR = -RHPI*DSIN(ARG) 299 IF (MOD(INUH,2).EQ.0) GO TO 120 300 C ZNR = -ZNR 301 C 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 309 C STR = CYR(I)*ZNR - CYI(I)*ZNI 310 C CYI(I) = CYR(I)*ZNI + CYI(I)*ZNR 311 C CYR(I) = STR 312 C STR = -ZNI*ZTI 313 C ZNI = ZNR*ZTI 314 C 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