GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zbesy.f
Go to the documentation of this file.
1  SUBROUTINE zbesy(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, CWRKR, CWRKI,
2  * IERR)
3 C***BEGIN PROLOGUE ZBESY
4 C***DATE WRITTEN 830501 (YYMMDD)
5 C***REVISION DATE 890801 (YYMMDD)
6 C***CATEGORY NO. B5K
7 C***KEYWORDS Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT,
8 C BESSEL FUNCTION OF SECOND KIND
9 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
10 C***PURPOSE TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT
11 C***DESCRIPTION
12 C
13 C ***A DOUBLE PRECISION ROUTINE***
14 C
15 C ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
16 C BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE
17 C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
18 C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESY RETURNS THE SCALED
19 C FUNCTIONS
20 C
21 C CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z)
22 C
23 C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
24 C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
25 C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
26 C (REF. 1).
27 C
28 C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION
29 C ZR,ZI - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0),
30 C -PI.LT.ARG(Z).LE.PI
31 C FNU - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0D0
32 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
33 C KODE= 1 RETURNS
34 C CY(I)=Y(FNU+I-1,Z), I=1,...,N
35 C = 2 RETURNS
36 C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N
37 C WHERE Y=AIMAG(Z)
38 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
39 C CWRKR, - DOUBLE PRECISION WORK VECTORS OF DIMENSION AT
40 C CWRKI AT LEAST N
41 C
42 C OUTPUT CYR,CYI ARE DOUBLE PRECISION
43 C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
44 C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
45 C CY(I)=Y(FNU+I-1,Z) OR
46 C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N
47 C DEPENDING ON KODE.
48 C NZ - NZ=0 , A NORMAL RETURN
49 C NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO
50 C UNDERFLOW (GENERALLY ON KODE=2)
51 C IERR - ERROR FLAG
52 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
53 C IERR=1, INPUT ERROR - NO COMPUTATION
54 C IERR=2, OVERFLOW - NO COMPUTATION, FNU IS
55 C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
56 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
57 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
58 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE
59 C ACCURACY
60 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
61 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
62 C CANCE BY ARGUMENT REDUCTION
63 C IERR=5, ERROR - NO COMPUTATION,
64 C ALGORITHM TERMINATION CONDITION NOT MET
65 C
66 C***LONG DESCRIPTION
67 C
68 C THE COMPUTATION IS CARRIED OUT BY THE FORMULA
69 C
70 C Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I
71 C
72 C WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z)
73 C AND H(2,FNU,Z) ARE CALCULATED IN CBESH.
74 C
75 C FOR NEGATIVE ORDERS,THE FORMULA
76 C
77 C Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
78 C
79 C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD
80 C INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE
81 C POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)*
82 C SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS
83 C NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A
84 C LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM
85 C CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS,
86 C WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF
87 C ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z).
88 C
89 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
90 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
91 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
92 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
93 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
94 C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
95 C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
96 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
97 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
98 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
99 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
100 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
101 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
102 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
103 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
104 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
105 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
106 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
107 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
108 C
109 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
110 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
111 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
112 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
113 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
114 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
115 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
116 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
117 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
118 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
119 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
120 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
121 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
122 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
123 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
124 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
125 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
126 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
127 C OR -PI/2+P.
128 C
129 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
130 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
131 C COMMERCE, 1955.
132 C
133 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
134 C BY D. E. AMOS, SAND83-0083, MAY, 1983.
135 C
136 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
137 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
138 C
139 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
140 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
141 C 1018, MAY, 1985
142 C
143 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
144 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
145 C MATH. SOFTWARE, 1986
146 C
147 C***ROUTINES CALLED ZBESH,I1MACH,D1MACH
148 C***END PROLOGUE ZBESY
149 C
150 C COMPLEX CWRK,CY,C1,C2,EX,HCI,Z,ZU,ZV
151  DOUBLE PRECISION CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2R,
152  * elim, exi, exr, ey, fnu, hcii, sti, str, tay, zi, zr, dexp,
153  * d1mach, ascle, rtol, atol, aa, bb, tol
154  INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH
155  dimension cyr(n), cyi(n), cwrkr(n), cwrki(n)
156 C***FIRST EXECUTABLE STATEMENT ZBESY
157  ierr = 0
158  nz=0
159  IF (zr.EQ.0.0d0 .AND. zi.EQ.0.0d0) ierr=1
160  IF (fnu.LT.0.0d0) ierr=1
161  IF (kode.LT.1 .OR. kode.GT.2) ierr=1
162  IF (n.LT.1) ierr=1
163  IF (ierr.NE.0) RETURN
164  hcii = 0.5d0
165  CALL zbesh(zr, zi, fnu, kode, 1, n, cyr, cyi, nz1, ierr)
166  IF (ierr.NE.0.AND.ierr.NE.3) GO TO 170
167  CALL zbesh(zr, zi, fnu, kode, 2, n, cwrkr, cwrki, nz2, ierr)
168  IF (ierr.NE.0.AND.ierr.NE.3) GO TO 170
169  nz = min0(nz1,nz2)
170  IF (kode.EQ.2) GO TO 60
171  DO 50 i=1,n
172  str = cwrkr(i) - cyr(i)
173  sti = cwrki(i) - cyi(i)
174  cyr(i) = -sti*hcii
175  cyi(i) = str*hcii
176  50 CONTINUE
177  RETURN
178  60 CONTINUE
179  tol = dmax1(d1mach(4),1.0d-18)
180  k1 = i1mach(15)
181  k2 = i1mach(16)
182  k = min0(iabs(k1),iabs(k2))
183  r1m5 = d1mach(5)
184 C-----------------------------------------------------------------------
185 C ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
186 C-----------------------------------------------------------------------
187  elim = 2.303d0*(dble(float(k))*r1m5-3.0d0)
188  exr = dcos(zr)
189  exi = dsin(zr)
190  ey = 0.0d0
191  tay = dabs(zi+zi)
192  IF (tay.LT.elim) ey = dexp(-tay)
193  IF (zi.LT.0.0d0) GO TO 90
194  c1r = exr*ey
195  c1i = exi*ey
196  c2r = exr
197  c2i = -exi
198  70 CONTINUE
199  nz = 0
200  rtol = 1.0d0/tol
201  ascle = d1mach(1)*rtol*1.0d+3
202  DO 80 i=1,n
203 C STR = C1R*CYR(I) - C1I*CYI(I)
204 C STI = C1R*CYI(I) + C1I*CYR(I)
205 C STR = -STR + C2R*CWRKR(I) - C2I*CWRKI(I)
206 C STI = -STI + C2R*CWRKI(I) + C2I*CWRKR(I)
207 C CYR(I) = -STI*HCII
208 C CYI(I) = STR*HCII
209  aa = cwrkr(i)
210  bb = cwrki(i)
211  atol = 1.0d0
212  IF (dmax1(dabs(aa),dabs(bb)).GT.ascle) GO TO 75
213  aa = aa*rtol
214  bb = bb*rtol
215  atol = tol
216  75 CONTINUE
217  str = (aa*c2r - bb*c2i)*atol
218  sti = (aa*c2i + bb*c2r)*atol
219  aa = cyr(i)
220  bb = cyi(i)
221  atol = 1.0d0
222  IF (dmax1(dabs(aa),dabs(bb)).GT.ascle) GO TO 85
223  aa = aa*rtol
224  bb = bb*rtol
225  atol = tol
226  85 CONTINUE
227  str = str - (aa*c1r - bb*c1i)*atol
228  sti = sti - (aa*c1i + bb*c1r)*atol
229  cyr(i) = -sti*hcii
230  cyi(i) = str*hcii
231  IF (str.EQ.0.0d0 .AND. sti.EQ.0.0d0 .AND. ey.EQ.0.0d0) nz = nz
232  * + 1
233  80 CONTINUE
234  RETURN
235  90 CONTINUE
236  c1r = exr
237  c1i = exi
238  c2r = exr*ey
239  c2i = -exi*ey
240  GO TO 70
241  170 CONTINUE
242  nz = 0
243  RETURN
244  END
double precision function d1mach(i)
Definition: d1mach.f:20
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)