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
cbesh.f
Go to the documentation of this file.
1  SUBROUTINE cbesh(Z, FNU, KODE, M, N, CY, NZ, IERR)
2 C***BEGIN PROLOGUE CBESH
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 ON KODE=1, CBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
13 C HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1
14 C OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX
15 C Z.NE.CMPLX(0.0E0,0.0E0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI.
16 C ON KODE=2, CBESH COMPUTES THE SCALED HANKEL FUNCTIONS
17 C
18 C CY(I)=H(M,FNU+J-1,Z)*EXP(-MM*Z*I) MM=3-2M, I**2=-1.
19 C
20 C WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER
21 C AND LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN
22 C THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1).
23 C
24 C INPUT
25 C Z - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI
26 C FNU - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0E0
27 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
28 C KODE= 1 RETURNS
29 C CY(J)=H(M,FNU+J-1,Z), J=1,...,N
30 C = 2 RETURNS
31 C CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M))
32 C J=1,...,N , I**2=-1
33 C M - KIND OF HANKEL FUNCTION, M=1 OR 2
34 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
35 C
36 C OUTPUT
37 C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
38 C VALUES FOR THE SEQUENCE
39 C CY(J)=H(M,FNU+J-1,Z) OR
40 C CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) J=1,...,N
41 C DEPENDING ON KODE, I**2=-1.
42 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
43 C NZ= 0 , NORMAL RETURN
44 C NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO
45 C DUE TO UNDERFLOW, CY(J)=CMPLX(0.0,0.0)
46 C J=1,...,NZ WHEN Y.GT.0.0 AND M=1 OR
47 C Y.LT.0.0 AND M=2. FOR THE COMPLMENTARY
48 C HALF PLANES, NZ STATES ONLY THE NUMBER
49 C OF UNDERFLOWS.
50 C IERR -ERROR FLAG
51 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
52 C IERR=1, INPUT ERROR - NO COMPUTATION
53 C IERR=2, OVERFLOW - NO COMPUTATION, FNU+N-1 TOO
54 C LARGE OR CABS(Z) TOO SMALL OR BOTH
55 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
56 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
57 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE
58 C ACCURACY
59 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
60 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
61 C CANCE BY ARGUMENT REDUCTION
62 C IERR=5, ERROR - NO COMPUTATION,
63 C ALGORITHM TERMINATION CONDITION NOT MET
64 C
65 C***LONG DESCRIPTION
66 C
67 C THE COMPUTATION IS CARRIED OUT BY THE RELATION
68 C
69 C H(M,FNU,Z)=(1/MP)*EXP(-MP*FNU)*K(FNU,Z*EXP(-MP))
70 C MP=MM*HPI*I, MM=3-2*M, HPI=PI/2, I**2=-1
71 C
72 C FOR M=1 OR 2 WHERE THE K BESSEL FUNCTION IS COMPUTED FOR THE
73 C RIGHT HALF PLANE RE(Z).GE.0.0. THE K FUNCTION IS CONTINUED
74 C TO THE LEFT HALF PLANE BY THE RELATION
75 C
76 C K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z)
77 C MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1
78 C
79 C WHERE I(FNU,Z) IS THE I BESSEL FUNCTION.
80 C
81 C EXPONENTIAL DECAY OF H(M,FNU,Z) OCCURS IN THE UPPER HALF Z
82 C PLANE FOR M=1 AND THE LOWER HALF Z PLANE FOR M=2. EXPONENTIAL
83 C GROWTH OCCURS IN THE COMPLEMENTARY HALF PLANES. SCALING
84 C BY EXP(-MM*Z*I) REMOVES THE EXPONENTIAL BEHAVIOR IN THE
85 C WHOLE Z PLANE FOR Z TO INFINITY.
86 C
87 C FOR NEGATIVE ORDERS,THE FORMULAE
88 C
89 C H(1,-FNU,Z) = H(1,FNU,Z)*CEXP( PI*FNU*I)
90 C H(2,-FNU,Z) = H(2,FNU,Z)*CEXP(-PI*FNU*I)
91 C I**2=-1
92 C
93 C CAN BE USED.
94 C
95 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
96 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
97 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
98 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
99 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
100 C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO
101 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
102 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
103 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
104 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
105 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
106 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
107 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
108 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
109 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
110 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
111 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
112 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
113 C
114 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
115 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
116 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
117 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
118 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
119 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
120 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
121 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
122 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
123 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
124 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
125 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
126 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
127 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
128 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
129 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
130 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
131 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
132 C OR -PI/2+P.
133 C
134 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
135 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
136 C COMMERCE, 1955.
137 C
138 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
139 C BY D. E. AMOS, SAND83-0083, MAY, 1983.
140 C
141 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
142 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
143 C
144 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
145 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
146 C 1018, MAY, 1985
147 C
148 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
149 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
150 C MATH. SOFTWARE, 1986
151 C
152 C***ROUTINES CALLED CACON,CBKNU,CBUNK,CUOIK,I1MACH,R1MACH
153 C***END PROLOGUE CBESH
154 C
155  COMPLEX cy, z, zn, zt, csgn
156  REAL aa, alim, aln, arg, az, cpn, dig, elim, fmm, fn, fnu, fnul,
157  * hpi, rhpi, rl, r1m5, sgn, spn, tol, ufl, xn, xx, yn, yy, r1mach,
158  * bb, ascle, rtol, atol
159  INTEGER i, ierr, inu, inuh, ir, k, kode, k1, k2, m,
160  * mm, mr, n, nn, nuf, nw, nz, i1mach
161  dimension cy(n)
162 C
163  DATA hpi /1.57079632679489662e0/
164 C
165 C***FIRST EXECUTABLE STATEMENT CBESH
166  nz=0
167  xx = REAL(z)
168  yy = aimag(z)
169  ierr = 0
170  IF (xx.EQ.0.0e0 .AND. yy.EQ.0.0e0) ierr=1
171  IF (fnu.LT.0.0e0) ierr=1
172  IF (m.LT.1 .OR. m.GT.2) ierr=1
173  IF (kode.LT.1 .OR. kode.GT.2) ierr=1
174  IF (n.LT.1) ierr=1
175  IF (ierr.NE.0) RETURN
176  nn = n
177 C-----------------------------------------------------------------------
178 C SET PARAMETERS RELATED TO MACHINE CONSTANTS.
179 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
180 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
181 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
182 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
183 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
184 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
185 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
186 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
187 C-----------------------------------------------------------------------
188  tol = amax1(r1mach(4),1.0e-18)
189  k1 = i1mach(12)
190  k2 = i1mach(13)
191  r1m5 = r1mach(5)
192  k = min0(iabs(k1),iabs(k2))
193  elim = 2.303e0*(float(k)*r1m5-3.0e0)
194  k1 = i1mach(11) - 1
195  aa = r1m5*float(k1)
196  dig = amin1(aa,18.0e0)
197  aa = aa*2.303e0
198  alim = elim + amax1(-aa,-41.45e0)
199  fnul = 10.0e0 + 6.0e0*(dig-3.0e0)
200  rl = 1.2e0*dig + 3.0e0
201  fn = fnu + float(nn-1)
202  mm = 3 - m - m
203  fmm = float(mm)
204  zn = z*cmplx(0.0e0,-fmm)
205  xn = REAL(zn)
206  yn = aimag(zn)
207  az = cabs(z)
208 C-----------------------------------------------------------------------
209 C TEST FOR RANGE
210 C-----------------------------------------------------------------------
211  aa = 0.5e0/tol
212  bb=float(i1mach(9))*0.5e0
213  aa=amin1(aa,bb)
214  IF(az.GT.aa) go to 240
215  IF(fn.GT.aa) go to 240
216  aa=sqrt(aa)
217  IF(az.GT.aa) ierr=3
218  IF(fn.GT.aa) ierr=3
219 C-----------------------------------------------------------------------
220 C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
221 C-----------------------------------------------------------------------
222  ufl = r1mach(1)*1.0e+3
223  IF (az.LT.ufl) go to 220
224  IF (fnu.GT.fnul) go to 90
225  IF (fn.LE.1.0e0) go to 70
226  IF (fn.GT.2.0e0) go to 60
227  IF (az.GT.tol) go to 70
228  arg = 0.5e0*az
229  aln = -fn*alog(arg)
230  IF (aln.GT.elim) go to 220
231  go to 70
232  60 CONTINUE
233  CALL cuoik(zn, fnu, kode, 2, nn, cy, nuf, tol, elim, alim)
234  IF (nuf.LT.0) go to 220
235  nz = nz + nuf
236  nn = nn - nuf
237 C-----------------------------------------------------------------------
238 C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
239 C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
240 C-----------------------------------------------------------------------
241  IF (nn.EQ.0) go to 130
242  70 CONTINUE
243  IF ((xn.LT.0.0e0) .OR. (xn.EQ.0.0e0 .AND. yn.LT.0.0e0 .AND.
244  * m.EQ.2)) go to 80
245 C-----------------------------------------------------------------------
246 C RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR.
247 C YN.GE.0. .OR. M=1)
248 C-----------------------------------------------------------------------
249  CALL cbknu(zn, fnu, kode, nn, cy, nz, tol, elim, alim)
250  go to 110
251 C-----------------------------------------------------------------------
252 C LEFT HALF PLANE COMPUTATION
253 C-----------------------------------------------------------------------
254  80 CONTINUE
255  mr = -mm
256  CALL cacon(zn, fnu, kode, mr, nn, cy, nw, rl, fnul, tol, elim,
257  * alim)
258  IF (nw.LT.0) go to 230
259  nz=nw
260  go to 110
261  90 CONTINUE
262 C-----------------------------------------------------------------------
263 C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
264 C-----------------------------------------------------------------------
265  mr = 0
266  IF ((xn.GE.0.0e0) .AND. (xn.NE.0.0e0 .OR. yn.GE.0.0e0 .OR.
267  * m.NE.2)) go to 100
268  mr = -mm
269  IF (xn.EQ.0.0e0 .AND. yn.LT.0.0e0) zn = -zn
270  100 CONTINUE
271  CALL cbunk(zn, fnu, kode, mr, nn, cy, nw, tol, elim, alim)
272  IF (nw.LT.0) go to 230
273  nz = nz + nw
274  110 CONTINUE
275 C-----------------------------------------------------------------------
276 C H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT)
277 C
278 C ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2
279 C-----------------------------------------------------------------------
280  sgn = sign(hpi,-fmm)
281 C-----------------------------------------------------------------------
282 C CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
283 C WHEN FNU IS LARGE
284 C-----------------------------------------------------------------------
285  inu = int(fnu)
286  inuh = inu/2
287  ir = inu - 2*inuh
288  arg = (fnu-float(inu-ir))*sgn
289  rhpi = 1.0e0/sgn
290  cpn = rhpi*cos(arg)
291  spn = rhpi*sin(arg)
292 C ZN = CMPLX(-SPN,CPN)
293  csgn = cmplx(-spn,cpn)
294 C IF (MOD(INUH,2).EQ.1) ZN = -ZN
295  IF (mod(inuh,2).EQ.1) csgn = -csgn
296  zt = cmplx(0.0e0,-fmm)
297  rtol = 1.0e0/tol
298  ascle = ufl*rtol
299  DO 120 i=1,nn
300 C CY(I) = CY(I)*ZN
301 C ZN = ZN*ZT
302  zn=cy(i)
303  aa=REAL(zn)
304  bb=aimag(zn)
305  atol=1.0e0
306  IF (amax1(abs(aa),abs(bb)).GT.ascle) go to 125
307  zn = zn*cmplx(rtol,0.0e0)
308  atol = tol
309  125 CONTINUE
310  zn = zn*csgn
311  cy(i) = zn*cmplx(atol,0.0e0)
312  csgn = csgn*zt
313  120 CONTINUE
314  RETURN
315  130 CONTINUE
316  IF (xn.LT.0.0e0) go to 220
317  RETURN
318  220 CONTINUE
319  ierr=2
320  nz=0
321  RETURN
322  230 CONTINUE
323  IF(nw.EQ.(-1)) go to 220
324  nz=0
325  ierr=5
326  RETURN
327  240 CONTINUE
328  nz=0
329  ierr=4
330  RETURN
331  END
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2