2
SUBROUTINE ZBESH (ZR, ZI, FNU, KODE, M, N, CYR, CYI, NZ, IERR)
3
C***BEGIN PROLOGUE ZBESH
4
C***PURPOSE Compute a sequence of the Hankel functions H(m,a,z)
5
C for superscript m=1 or 2, real nonnegative orders a=b,
6
C b+1,... where b>0, and nonzero complex argument z. A
7
C scaling option is available to help avoid overflow.
10
C***TYPE COMPLEX (CBESH-C, ZBESH-C)
11
C***KEYWORDS BESSEL FUNCTIONS OF COMPLEX ARGUMENT,
12
C BESSEL FUNCTIONS OF THE THIRD KIND, H BESSEL FUNCTIONS,
14
C***AUTHOR Amos, D. E., (SNL)
17
C ***A DOUBLE PRECISION ROUTINE***
18
C On KODE=1, ZBESH computes an N member sequence of complex
19
C Hankel (Bessel) functions CY(L)=H(M,FNU+L-1,Z) for super-
20
C script M=1 or 2, real nonnegative orders FNU+L-1, L=1,...,
21
C N, and complex nonzero Z in the cut plane -pi<arg(Z)<=pi
22
C where Z=ZR+i*ZI. On KODE=2, CBESH returns the scaled
25
C CY(L) = H(M,FNU+L-1,Z)*exp(-(3-2*M)*Z*i), i**2=-1
27
C which removes the exponential behavior in both the upper
28
C and lower half planes. Definitions and notation are found
29
C in the NBS Handbook of Mathematical Functions (Ref. 1).
32
C ZR - DOUBLE PRECISION real part of nonzero argument Z
33
C ZI - DOUBLE PRECISION imag part of nonzero argument Z
34
C FNU - DOUBLE PRECISION initial order, FNU>=0
35
C KODE - A parameter to indicate the scaling option
37
C CY(L)=H(M,FNU+L-1,Z), L=1,...,N
39
C CY(L)=H(M,FNU+L-1,Z)*exp(-(3-2M)*Z*i),
41
C M - Superscript of Hankel function, M=1 or 2
42
C N - Number of terms in the sequence, N>=1
45
C CYR - DOUBLE PRECISION real part of result vector
46
C CYI - DOUBLE PRECISION imag part of result vector
47
C NZ - Number of underflows set to zero
49
C NZ>0 CY(L)=0 for NZ values of L (if M=1 and
50
C Im(Z)>0 or if M=2 and Im(Z)<0, then
51
C CY(L)=0 for L=1,...,NZ; in the com-
52
C plementary half planes, the underflows
53
C may not be in an uninterrupted sequence)
55
C IERR=0 Normal return - COMPUTATION COMPLETED
56
C IERR=1 Input error - NO COMPUTATION
57
C IERR=2 Overflow - NO COMPUTATION
58
C (abs(Z) too small and/or FNU+N-1
60
C IERR=3 Precision warning - COMPUTATION COMPLETED
61
C (Result has half precision or less
62
C because abs(Z) or FNU+N-1 is large)
63
C IERR=4 Precision error - NO COMPUTATION
64
C (Result has no precision because
65
C abs(Z) or FNU+N-1 is too large)
66
C IERR=5 Algorithmic error - NO COMPUTATION
67
C (Termination condition not met)
71
C The computation is carried out by the formula
73
C H(m,a,z) = (1/t)*exp(-a*t)*K(a,z*exp(-t))
76
C where the K Bessel function is computed as described in the
79
C Exponential decay of H(m,a,z) occurs in the upper half z
80
C plane for m=1 and the lower half z plane for m=2. Exponential
81
C growth occurs in the complementary half planes. Scaling
82
C by exp(-(3-2*m)*z*i) removes the exponential behavior in the
83
C whole z plane as z goes to infinity.
85
C For negative orders, the formula
87
C H(m,-a,z) = H(m,a,z)*exp((3-2*m)*a*pi*i)
91
C In most complex variable computation, one must evaluate ele-
92
C mentary functions. When the magnitude of Z or FNU+N-1 is
93
C large, losses of significance by argument reduction occur.
94
C Consequently, if either one exceeds U1=SQRT(0.5/UR), then
95
C losses exceeding half precision are likely and an error flag
96
C IERR=3 is triggered where UR=MAX(D1MACH(4),1.0D-18) is double
97
C precision unit roundoff limited to 18 digits precision. Also,
98
C if either is larger than U2=0.5/UR, then all significance is
99
C lost and IERR=4. In order to use the INT function, arguments
100
C must be further restricted not to exceed the largest machine
101
C integer, U3=I1MACH(9). Thus, the magnitude of Z and FNU+N-1
102
C is restricted by MIN(U2,U3). In IEEE arithmetic, U1,U2, and
103
C U3 approximate 2.0E+3, 4.2E+6, 2.1E+9 in single precision
104
C and 4.7E+7, 2.3E+15 and 2.1E+9 in double precision. This
105
C makes U2 limiting in single precision and U3 limiting in
106
C double precision. This means that one can expect to retain,
107
C in the worst cases on IEEE machines, no digits in single pre-
108
C cision and only 6 digits in double precision. Similar con-
109
C siderations hold for other machines.
111
C The approximate relative error in the magnitude of a complex
112
C Bessel function can be expressed as P*10**S where P=MAX(UNIT
113
C ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre-
114
C sents the increase in error due to argument reduction in the
115
C elementary functions. Here, S=MAX(1,ABS(LOG10(ABS(Z))),
116
C ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF
117
C ABS(Z),ABS(EXPONENT OF FNU)) ). However, the phase angle may
118
C have only absolute accuracy. This is most likely to occur
119
C when one component (in magnitude) is larger than the other by
120
C several orders of magnitude. If one component is 10**K larger
121
C than the other, then one can expect only MAX(ABS(LOG10(P))-K,
122
C 0) significant digits; or, stated another way, when K exceeds
123
C the exponent of P, no significant digits remain in the smaller
124
C component. However, the phase angle retains absolute accuracy
125
C because, in complex arithmetic with precision P, the smaller
126
C component will not (as a rule) decrease below P times the
127
C magnitude of the larger component. In these extreme cases,
128
C the principal phase angle is on the order of +P, -P, PI/2-P,
131
C***REFERENCES 1. M. Abramowitz and I. A. Stegun, Handbook of Mathe-
132
C matical Functions, National Bureau of Standards
133
C Applied Mathematics Series 55, U. S. Department
134
C of Commerce, Tenth Printing (1972) or later.
135
C 2. D. E. Amos, Computation of Bessel Functions of
136
C Complex Argument, Report SAND83-0086, Sandia National
137
C Laboratories, Albuquerque, NM, May 1983.
138
C 3. D. E. Amos, Computation of Bessel Functions of
139
C Complex Argument and Large Order, Report SAND83-0643,
140
C Sandia National Laboratories, Albuquerque, NM, May
142
C 4. D. E. Amos, A Subroutine Package for Bessel Functions
143
C of a Complex Argument and Nonnegative Order, Report
144
C SAND85-1018, Sandia National Laboratory, Albuquerque,
146
C 5. D. E. Amos, A portable package for Bessel functions
147
C of a complex argument and nonnegative order, ACM
148
C Transactions on Mathematical Software, 12 (September
149
C 1986), pp. 265-273.
151
C***ROUTINES CALLED D1MACH, I1MACH, ZABS, ZACON, ZBKNU, ZBUNK, ZUOIK
152
C***REVISION HISTORY (YYMMDD)
153
C 830501 DATE WRITTEN
154
C 890801 REVISION DATE from Version 3.2
155
C 910415 Prologue converted to Version 4.0 format. (BAB)
156
C 920128 Category corrected. (WRB)
157
C 920811 Prologue revised. (DWL)
158
C***END PROLOGUE ZBESH
160
C COMPLEX CY,Z,ZN,ZT,CSGN
161
DOUBLE PRECISION AA, ALIM, ALN, ARG, AZ, CYI, CYR, DIG, ELIM,
162
* FMM, FN, FNU, FNUL, HPI, RHPI, RL, R1M5, SGN, STR, TOL, UFL, ZI,
163
* ZNI, ZNR, ZR, ZTI, D1MACH, ZABS, BB, ASCLE, RTOL, ATOL, STI,
165
INTEGER I, IERR, INU, INUH, IR, K, KODE, K1, K2, M,
166
* MM, MR, N, NN, NUF, NW, NZ, I1MACH
167
DIMENSION CYR(N), CYI(N)
170
DATA HPI /1.57079632679489662D0/
172
C***FIRST EXECUTABLE STATEMENT ZBESH
175
IF (ZR.EQ.0.0D0 .AND. ZI.EQ.0.0D0) IERR=1
176
IF (FNU.LT.0.0D0) IERR=1
177
IF (M.LT.1 .OR. M.GT.2) IERR=1
178
IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
180
IF (IERR.NE.0) RETURN
182
C-----------------------------------------------------------------------
183
C SET PARAMETERS RELATED TO MACHINE CONSTANTS.
184
C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
185
C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
186
C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
187
C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
188
C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
189
C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
190
C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
191
C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
192
C-----------------------------------------------------------------------
193
TOL = MAX(D1MACH(4),1.0D-18)
197
K = MIN(ABS(K1),ABS(K2))
198
ELIM = 2.303D0*(K*R1M5-3.0D0)
203
ALIM = ELIM + MAX(-AA,-41.45D0)
204
FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
205
RL = 1.2D0*DIG + 3.0D0
211
C-----------------------------------------------------------------------
212
C TEST FOR PROPER RANGE
213
C-----------------------------------------------------------------------
218
IF (AZ.GT.AA) GO TO 260
219
IF (FN.GT.AA) GO TO 260
223
C-----------------------------------------------------------------------
224
C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
225
C-----------------------------------------------------------------------
226
UFL = D1MACH(1)*1.0D+3
227
IF (AZ.LT.UFL) GO TO 230
228
IF (FNU.GT.FNUL) GO TO 90
229
IF (FN.LE.1.0D0) GO TO 70
230
IF (FN.GT.2.0D0) GO TO 60
231
IF (AZ.GT.TOL) GO TO 70
234
IF (ALN.GT.ELIM) GO TO 230
237
CALL ZUOIK(ZNR, ZNI, FNU, KODE, 2, NN, CYR, CYI, NUF, TOL, ELIM,
239
IF (NUF.LT.0) GO TO 230
242
C-----------------------------------------------------------------------
243
C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
244
C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
245
C-----------------------------------------------------------------------
246
IF (NN.EQ.0) GO TO 140
248
IF ((ZNR.LT.0.0D0) .OR. (ZNR.EQ.0.0D0 .AND. ZNI.LT.0.0D0 .AND.
250
C-----------------------------------------------------------------------
251
C RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR.
253
C-----------------------------------------------------------------------
254
CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NZ, TOL, ELIM, ALIM)
256
C-----------------------------------------------------------------------
257
C LEFT HALF PLANE COMPUTATION
258
C-----------------------------------------------------------------------
261
CALL ZACON(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, RL, FNUL,
263
IF (NW.LT.0) GO TO 240
267
C-----------------------------------------------------------------------
268
C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
269
C-----------------------------------------------------------------------
271
IF ((ZNR.GE.0.0D0) .AND. (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0 .OR.
274
IF (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0) GO TO 100
278
CALL ZBUNK(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, TOL, ELIM,
280
IF (NW.LT.0) GO TO 240
283
C-----------------------------------------------------------------------
284
C H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT)
286
C ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2
287
C-----------------------------------------------------------------------
288
SGN = DSIGN(HPI,-FMM)
289
C-----------------------------------------------------------------------
290
C CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
292
C-----------------------------------------------------------------------
296
ARG = (FNU-(INU-IR))*SGN
298
C ZNI = RHPI*COS(ARG)
299
C ZNR = -RHPI*SIN(ARG)
300
CSGNI = RHPI*COS(ARG)
301
CSGNR = -RHPI*SIN(ARG)
302
IF (MOD(INUH,2).EQ.0) GO TO 120
312
C STR = CYR(I)*ZNR - CYI(I)*ZNI
313
C CYI(I) = CYR(I)*ZNI + CYI(I)*ZNR
321
IF (MAX(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 135
326
STR = AA*CSGNR - BB*CSGNI
327
STI = AA*CSGNI + BB*CSGNR
336
IF (ZNR.LT.0.0D0) GO TO 230
343
IF(NW.EQ.(-1)) GO TO 230
added
removed
deps/openlibm/slatec/zbesh.f
