GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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  35 CONTINUE
223  ufl = r1mach(1)*1.0e+3
224  IF (az.LT.ufl) GO TO 220
225  IF (fnu.GT.fnul) GO TO 90
226  IF (fn.LE.1.0e0) GO TO 70
227  IF (fn.GT.2.0e0) GO TO 60
228  IF (az.GT.tol) GO TO 70
229  arg = 0.5e0*az
230  aln = -fn*alog(arg)
231  IF (aln.GT.elim) GO TO 220
232  GO TO 70
233  60 CONTINUE
234  CALL cuoik(zn, fnu, kode, 2, nn, cy, nuf, tol, elim, alim)
235  IF (nuf.LT.0) GO TO 220
236  nz = nz + nuf
237  nn = nn - nuf
238 C-----------------------------------------------------------------------
239 C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
240 C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
241 C-----------------------------------------------------------------------
242  IF (nn.EQ.0) GO TO 130
243  70 CONTINUE
244  IF ((xn.LT.0.0e0) .OR. (xn.EQ.0.0e0 .AND. yn.LT.0.0e0 .AND.
245  * m.EQ.2)) GO TO 80
246 C-----------------------------------------------------------------------
247 C RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR.
248 C YN.GE.0. .OR. M=1)
249 C-----------------------------------------------------------------------
250  CALL cbknu(zn, fnu, kode, nn, cy, nz, tol, elim, alim)
251  GO TO 110
252 C-----------------------------------------------------------------------
253 C LEFT HALF PLANE COMPUTATION
254 C-----------------------------------------------------------------------
255  80 CONTINUE
256  mr = -mm
257  CALL cacon(zn, fnu, kode, mr, nn, cy, nw, rl, fnul, tol, elim,
258  * alim)
259  IF (nw.LT.0) GO TO 230
260  nz=nw
261  GO TO 110
262  90 CONTINUE
263 C-----------------------------------------------------------------------
264 C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
265 C-----------------------------------------------------------------------
266  mr = 0
267  IF ((xn.GE.0.0e0) .AND. (xn.NE.0.0e0 .OR. yn.GE.0.0e0 .OR.
268  * m.NE.2)) GO TO 100
269  mr = -mm
270  IF (xn.EQ.0.0e0 .AND. yn.LT.0.0e0) zn = -zn
271  100 CONTINUE
272  CALL cbunk(zn, fnu, kode, mr, nn, cy, nw, tol, elim, alim)
273  IF (nw.LT.0) GO TO 230
274  nz = nz + nw
275  110 CONTINUE
276 C-----------------------------------------------------------------------
277 C H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT)
278 C
279 C ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2
280 C-----------------------------------------------------------------------
281  sgn = sign(hpi,-fmm)
282 C-----------------------------------------------------------------------
283 C CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
284 C WHEN FNU IS LARGE
285 C-----------------------------------------------------------------------
286  inu = int(fnu)
287  inuh = inu/2
288  ir = inu - 2*inuh
289  arg = (fnu-float(inu-ir))*sgn
290  rhpi = 1.0e0/sgn
291  cpn = rhpi*cos(arg)
292  spn = rhpi*sin(arg)
293 C ZN = CMPLX(-SPN,CPN)
294  csgn = cmplx(-spn,cpn)
295 C IF (MOD(INUH,2).EQ.1) ZN = -ZN
296  IF (mod(inuh,2).EQ.1) csgn = -csgn
297  zt = cmplx(0.0e0,-fmm)
298  rtol = 1.0e0/tol
299  ascle = ufl*rtol
300  DO 120 i=1,nn
301 C CY(I) = CY(I)*ZN
302 C ZN = ZN*ZT
303  zn=cy(i)
304  aa=REAL(ZN)
305  bb=aimag(zn)
306  atol=1.0e0
307  IF (amax1(abs(aa),abs(bb)).GT.ascle) GO TO 125
308  zn = zn*cmplx(rtol,0.0e0)
309  atol = tol
310  125 CONTINUE
311  zn = zn*csgn
312  cy(i) = zn*cmplx(atol,0.0e0)
313  csgn = csgn*zt
314  120 CONTINUE
315  RETURN
316  130 CONTINUE
317  IF (xn.LT.0.0e0) GO TO 220
318  RETURN
319  220 CONTINUE
320  ierr=2
321  nz=0
322  RETURN
323  230 CONTINUE
324  IF(nw.EQ.(-1)) GO TO 220
325  nz=0
326  ierr=5
327  RETURN
328  240 CONTINUE
329  ierr=4
330  GO TO 35
331  END
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:860
static T abs(T x)
Definition: pr-output.cc:1696
double complex cmplx
Definition: Faddeeva.cc:217
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)