GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
zbesh.f
Go to the documentation of this file.
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,XZABS,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, 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)
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 = 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
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