|
3217
|
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 |
