GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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  35 CONTINUE
224  ufl = d1mach(1)*1.0d+3
225  IF (az.LT.ufl) GO TO 230
226  IF (fnu.GT.fnul) GO TO 90
227  IF (fn.LE.1.0d0) GO TO 70
228  IF (fn.GT.2.0d0) GO TO 60
229  IF (az.GT.tol) GO TO 70
230  arg = 0.5d0*az
231  aln = -fn*dlog(arg)
232  IF (aln.GT.elim) GO TO 230
233  GO TO 70
234  60 CONTINUE
235  CALL zuoik(znr, zni, fnu, kode, 2, nn, cyr, cyi, nuf, tol, elim,
236  * alim)
237  IF (nuf.LT.0) GO TO 230
238  nz = nz + nuf
239  nn = nn - nuf
240 C-----------------------------------------------------------------------
241 C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
242 C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
243 C-----------------------------------------------------------------------
244  IF (nn.EQ.0) GO TO 140
245  70 CONTINUE
246  IF ((znr.LT.0.0d0) .OR. (znr.EQ.0.0d0 .AND. zni.LT.0.0d0 .AND.
247  * m.EQ.2)) GO TO 80
248 C-----------------------------------------------------------------------
249 C RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR.
250 C YN.GE.0. .OR. M=1)
251 C-----------------------------------------------------------------------
252  CALL zbknu(znr, zni, fnu, kode, nn, cyr, cyi, nz, tol, elim, alim)
253  GO TO 110
254 C-----------------------------------------------------------------------
255 C LEFT HALF PLANE COMPUTATION
256 C-----------------------------------------------------------------------
257  80 CONTINUE
258  mr = -mm
259  CALL zacon(znr, zni, fnu, kode, mr, nn, cyr, cyi, nw, rl, fnul,
260  * tol, elim, alim)
261  IF (nw.LT.0) GO TO 240
262  nz=nw
263  GO TO 110
264  90 CONTINUE
265 C-----------------------------------------------------------------------
266 C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
267 C-----------------------------------------------------------------------
268  mr = 0
269  IF ((znr.GE.0.0d0) .AND. (znr.NE.0.0d0 .OR. zni.GE.0.0d0 .OR.
270  * m.NE.2)) GO TO 100
271  mr = -mm
272  IF (znr.NE.0.0d0 .OR. zni.GE.0.0d0) GO TO 100
273  znr = -znr
274  zni = -zni
275  100 CONTINUE
276  CALL zbunk(znr, zni, fnu, kode, mr, nn, cyr, cyi, nw, tol, elim,
277  * alim)
278  IF (nw.LT.0) GO TO 240
279  nz = nz + nw
280  110 CONTINUE
281 C-----------------------------------------------------------------------
282 C H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT)
283 C
284 C ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2
285 C-----------------------------------------------------------------------
286  sgn = dsign(hpi,-fmm)
287 C-----------------------------------------------------------------------
288 C CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
289 C WHEN FNU IS LARGE
290 C-----------------------------------------------------------------------
291  inu = int(sngl(fnu))
292  inuh = inu/2
293  ir = inu - 2*inuh
294  arg = (fnu-dble(float(inu-ir)))*sgn
295  rhpi = 1.0d0/sgn
296 C ZNI = RHPI*DCOS(ARG)
297 C ZNR = -RHPI*DSIN(ARG)
298  csgni = rhpi*dcos(arg)
299  csgnr = -rhpi*dsin(arg)
300  IF (mod(inuh,2).EQ.0) GO TO 120
301 C ZNR = -ZNR
302 C ZNI = -ZNI
303  csgnr = -csgnr
304  csgni = -csgni
305  120 CONTINUE
306  zti = -fmm
307  rtol = 1.0d0/tol
308  ascle = ufl*rtol
309  DO 130 i=1,nn
310 C STR = CYR(I)*ZNR - CYI(I)*ZNI
311 C CYI(I) = CYR(I)*ZNI + CYI(I)*ZNR
312 C CYR(I) = STR
313 C STR = -ZNI*ZTI
314 C ZNI = ZNR*ZTI
315 C ZNR = STR
316  aa = cyr(i)
317  bb = cyi(i)
318  atol = 1.0d0
319  IF (dmax1(dabs(aa),dabs(bb)).GT.ascle) GO TO 135
320  aa = aa*rtol
321  bb = bb*rtol
322  atol = tol
323  135 CONTINUE
324  str = aa*csgnr - bb*csgni
325  sti = aa*csgni + bb*csgnr
326  cyr(i) = str*atol
327  cyi(i) = sti*atol
328  str = -csgni*zti
329  csgni = csgnr*zti
330  csgnr = str
331  130 CONTINUE
332  RETURN
333  140 CONTINUE
334  IF (znr.LT.0.0d0) GO TO 230
335  RETURN
336  230 CONTINUE
337  nz=0
338  ierr=2
339  RETURN
340  240 CONTINUE
341  IF(nw.EQ.(-1)) GO TO 230
342  nz=0
343  ierr=5
344  RETURN
345  260 CONTINUE
346  ierr=4
347  GO TO 35
348  END
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:932
subroutine zacon(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: zacon.f:3
subroutine zbesh(ZR, ZI, FNU, KODE, M, N, CYR, CYI, NZ, IERR)
Definition: zbesh.f:2
subroutine zbknu(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM)
Definition: zbknu.f:3
subroutine zbunk(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, ALIM)
Definition: zbunk.f:3
subroutine zuoik(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL, ELIM, ALIM)
Definition: zuoik.f:3