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