|
3217
|
1 SUBROUTINE ZBESI(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR) |
|
|
2 C***BEGIN PROLOGUE ZBESI |
|
|
3 C***DATE WRITTEN 830501 (YYMMDD) |
|
|
4 C***REVISION DATE 890801 (YYMMDD) |
|
|
5 C***CATEGORY NO. B5K |
|
|
6 C***KEYWORDS I-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION, |
|
|
7 C MODIFIED BESSEL FUNCTION OF THE FIRST KIND |
|
|
8 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES |
|
|
9 C***PURPOSE TO COMPUTE I-BESSEL FUNCTIONS OF COMPLEX ARGUMENT |
|
|
10 C***DESCRIPTION |
|
|
11 C |
|
|
12 C ***A DOUBLE PRECISION ROUTINE*** |
|
|
13 C ON KODE=1, ZBESI COMPUTES AN N MEMBER SEQUENCE OF COMPLEX |
|
|
14 C BESSEL FUNCTIONS CY(J)=I(FNU+J-1,Z) FOR REAL, NONNEGATIVE |
|
|
15 C ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z IN THE CUT PLANE |
|
|
16 C -PI.LT.ARG(Z).LE.PI. ON KODE=2, ZBESI RETURNS THE SCALED |
|
|
17 C FUNCTIONS |
|
|
18 C |
|
|
19 C CY(J)=EXP(-ABS(X))*I(FNU+J-1,Z) J = 1,...,N , X=REAL(Z) |
|
|
20 C |
|
|
21 C WITH THE EXPONENTIAL GROWTH REMOVED IN BOTH THE LEFT AND |
|
|
22 C RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION |
|
|
23 C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS |
|
|
24 C (REF. 1). |
|
|
25 C |
|
|
26 C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION |
|
|
27 C ZR,ZI - Z=CMPLX(ZR,ZI), -PI.LT.ARG(Z).LE.PI |
|
|
28 C FNU - ORDER OF INITIAL I FUNCTION, FNU.GE.0.0D0 |
|
|
29 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION |
|
|
30 C KODE= 1 RETURNS |
|
|
31 C CY(J)=I(FNU+J-1,Z), J=1,...,N |
|
|
32 C = 2 RETURNS |
|
|
33 C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)), J=1,...,N |
|
|
34 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 |
|
|
35 C |
|
|
36 C OUTPUT CYR,CYI ARE DOUBLE PRECISION |
|
|
37 C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS |
|
|
38 C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE |
|
|
39 C CY(J)=I(FNU+J-1,Z) OR |
|
|
40 C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)) J=1,...,N |
|
|
41 C DEPENDING ON KODE, X=REAL(Z) |
|
|
42 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, |
|
|
43 C NZ= 0 , NORMAL RETURN |
|
|
44 C NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO |
|
|
45 C TO UNDERFLOW, CY(J)=CMPLX(0.0D0,0.0D0) |
|
|
46 C J = N-NZ+1,...,N |
|
|
47 C IERR - ERROR FLAG |
|
|
48 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED |
|
|
49 C IERR=1, INPUT ERROR - NO COMPUTATION |
|
|
50 C IERR=2, OVERFLOW - NO COMPUTATION, REAL(Z) TOO |
|
|
51 C LARGE ON KODE=1 |
|
|
52 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE |
|
|
53 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT |
|
|
54 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE |
|
|
55 C ACCURACY |
|
|
56 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- |
|
|
57 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- |
|
|
58 C CANCE BY ARGUMENT REDUCTION |
|
|
59 C IERR=5, ERROR - NO COMPUTATION, |
|
|
60 C ALGORITHM TERMINATION CONDITION NOT MET |
|
|
61 C |
|
|
62 C***LONG DESCRIPTION |
|
|
63 C |
|
|
64 C THE COMPUTATION IS CARRIED OUT BY THE POWER SERIES FOR |
|
|
65 C SMALL CABS(Z), THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z), |
|
|
66 C THE MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN AND A |
|
|
67 C NEUMANN SERIES FOR IMTERMEDIATE MAGNITUDES, AND THE |
|
|
68 C UNIFORM ASYMPTOTIC EXPANSIONS FOR I(FNU,Z) AND J(FNU,Z) |
|
|
69 C FOR LARGE ORDERS. BACKWARD RECURRENCE IS USED TO GENERATE |
|
|
70 C SEQUENCES OR REDUCE ORDERS WHEN NECESSARY. |
|
|
71 C |
|
|
72 C THE CALCULATIONS ABOVE ARE DONE IN THE RIGHT HALF PLANE AND |
|
|
73 C CONTINUED INTO THE LEFT HALF PLANE BY THE FORMULA |
|
|
74 C |
|
|
75 C I(FNU,Z*EXP(M*PI)) = EXP(M*PI*FNU)*I(FNU,Z) REAL(Z).GT.0.0 |
|
|
76 C M = +I OR -I, I**2=-1 |
|
|
77 C |
|
|
78 C FOR NEGATIVE ORDERS,THE FORMULA |
|
|
79 C |
|
|
80 C I(-FNU,Z) = I(FNU,Z) + (2/PI)*SIN(PI*FNU)*K(FNU,Z) |
|
|
81 C |
|
|
82 C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE |
|
|
83 C THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE |
|
|
84 C INTEGER,THE MAGNITUDE OF I(-FNU,Z)=I(FNU,Z) IS A LARGE |
|
|
85 C NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER, |
|
|
86 C K(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF |
|
|
87 C TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY |
|
|
88 C UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN |
|
|
89 C OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE, |
|
|
90 C LARGE MEANS FNU.GT.CABS(Z). |
|
|
91 C |
|
|
92 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- |
|
|
93 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS |
|
|
94 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. |
|
|
95 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN |
|
|
96 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG |
|
|
97 C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS |
|
|
98 C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. |
|
|
99 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS |
|
|
100 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS |
|
|
101 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE |
|
|
102 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS |
|
|
103 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 |
|
|
104 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION |
|
|
105 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION |
|
|
106 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN |
|
|
107 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT |
|
|
108 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS |
|
|
109 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. |
|
|
110 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. |
|
|
111 C |
|
|
112 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX |
|
|
113 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT |
|
|
114 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- |
|
|
115 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE |
|
|
116 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), |
|
|
117 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF |
|
|
118 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY |
|
|
119 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN |
|
|
120 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY |
|
|
121 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER |
|
|
122 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, |
|
|
123 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS |
|
|
124 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER |
|
|
125 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY |
|
|
126 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER |
|
|
127 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE |
|
|
128 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, |
|
|
129 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, |
|
|
130 C OR -PI/2+P. |
|
|
131 C |
|
|
132 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ |
|
|
133 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF |
|
|
134 C COMMERCE, 1955. |
|
|
135 C |
|
|
136 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT |
|
|
137 C BY D. E. AMOS, SAND83-0083, MAY, 1983. |
|
|
138 C |
|
|
139 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT |
|
|
140 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 |
|
|
141 C |
|
|
142 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX |
|
|
143 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- |
|
|
144 C 1018, MAY, 1985 |
|
|
145 C |
|
|
146 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX |
|
|
147 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. |
|
|
148 C MATH. SOFTWARE, 1986 |
|
|
149 C |
|
|
150 C***ROUTINES CALLED ZBINU,I1MACH,D1MACH |
|
|
151 C***END PROLOGUE ZBESI |
|
|
152 C COMPLEX CONE,CSGN,CW,CY,CZERO,Z,ZN |
|
|
153 DOUBLE PRECISION AA, ALIM, ARG, CONEI, CONER, CSGNI, CSGNR, CYI, |
|
|
154 * CYR, DIG, ELIM, FNU, FNUL, PI, RL, R1M5, STR, TOL, ZI, ZNI, ZNR, |
|
|
155 * ZR, D1MACH, AZ, BB, FN, XZABS, ASCLE, RTOL, ATOL, STI |
|
|
156 INTEGER I, IERR, INU, K, KODE, K1,K2,N,NZ,NN, I1MACH |
|
|
157 DIMENSION CYR(N), CYI(N) |
|
|
158 DATA PI /3.14159265358979324D0/ |
|
|
159 DATA CONER, CONEI /1.0D0,0.0D0/ |
|
|
160 C |
|
|
161 C***FIRST EXECUTABLE STATEMENT ZBESI |
|
|
162 IERR = 0 |
|
|
163 NZ=0 |
|
|
164 IF (FNU.LT.0.0D0) IERR=1 |
|
|
165 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 |
|
|
166 IF (N.LT.1) IERR=1 |
|
|
167 IF (IERR.NE.0) RETURN |
|
|
168 C----------------------------------------------------------------------- |
|
|
169 C SET PARAMETERS RELATED TO MACHINE CONSTANTS. |
|
|
170 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. |
|
|
171 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. |
|
|
172 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND |
|
|
173 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR |
|
|
174 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. |
|
|
175 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. |
|
|
176 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). |
|
|
177 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. |
|
|
178 C----------------------------------------------------------------------- |
|
|
179 TOL = DMAX1(D1MACH(4),1.0D-18) |
|
|
180 K1 = I1MACH(15) |
|
|
181 K2 = I1MACH(16) |
|
|
182 R1M5 = D1MACH(5) |
|
|
183 K = MIN0(IABS(K1),IABS(K2)) |
|
|
184 ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) |
|
|
185 K1 = I1MACH(14) - 1 |
|
|
186 AA = R1M5*DBLE(FLOAT(K1)) |
|
|
187 DIG = DMIN1(AA,18.0D0) |
|
|
188 AA = AA*2.303D0 |
|
|
189 ALIM = ELIM + DMAX1(-AA,-41.45D0) |
|
|
190 RL = 1.2D0*DIG + 3.0D0 |
|
|
191 FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) |
|
|
192 C----------------------------------------------------------------------------- |
|
|
193 C TEST FOR PROPER RANGE |
|
|
194 C----------------------------------------------------------------------- |
|
|
195 AZ = XZABS(ZR,ZI) |
|
|
196 FN = FNU+DBLE(FLOAT(N-1)) |
|
|
197 AA = 0.5D0/TOL |
|
|
198 BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 |
|
|
199 AA = DMIN1(AA,BB) |
|
|
200 IF (AZ.GT.AA) GO TO 260 |
|
|
201 IF (FN.GT.AA) GO TO 260 |
|
|
202 AA = DSQRT(AA) |
|
|
203 IF (AZ.GT.AA) IERR=3 |
|
|
204 IF (FN.GT.AA) IERR=3 |
|
|
205 ZNR = ZR |
|
|
206 ZNI = ZI |
|
|
207 CSGNR = CONER |
|
|
208 CSGNI = CONEI |
|
|
209 IF (ZR.GE.0.0D0) GO TO 40 |
|
|
210 ZNR = -ZR |
|
|
211 ZNI = -ZI |
|
|
212 C----------------------------------------------------------------------- |
|
|
213 C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE |
|
|
214 C WHEN FNU IS LARGE |
|
|
215 C----------------------------------------------------------------------- |
|
|
216 INU = INT(SNGL(FNU)) |
|
|
217 ARG = (FNU-DBLE(FLOAT(INU)))*PI |
|
|
218 IF (ZI.LT.0.0D0) ARG = -ARG |
|
|
219 CSGNR = DCOS(ARG) |
|
|
220 CSGNI = DSIN(ARG) |
|
|
221 IF (MOD(INU,2).EQ.0) GO TO 40 |
|
|
222 CSGNR = -CSGNR |
|
|
223 CSGNI = -CSGNI |
|
|
224 40 CONTINUE |
|
|
225 CALL ZBINU(ZNR, ZNI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, |
|
|
226 * ELIM, ALIM) |
|
|
227 IF (NZ.LT.0) GO TO 120 |
|
|
228 IF (ZR.GE.0.0D0) RETURN |
|
|
229 C----------------------------------------------------------------------- |
|
|
230 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE |
|
|
231 C----------------------------------------------------------------------- |
|
|
232 NN = N - NZ |
|
|
233 IF (NN.EQ.0) RETURN |
|
|
234 RTOL = 1.0D0/TOL |
|
|
235 ASCLE = D1MACH(1)*RTOL*1.0D+3 |
|
|
236 DO 50 I=1,NN |
|
|
237 C STR = CYR(I)*CSGNR - CYI(I)*CSGNI |
|
|
238 C CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR |
|
|
239 C CYR(I) = STR |
|
|
240 AA = CYR(I) |
|
|
241 BB = CYI(I) |
|
|
242 ATOL = 1.0D0 |
|
|
243 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 55 |
|
|
244 AA = AA*RTOL |
|
|
245 BB = BB*RTOL |
|
|
246 ATOL = TOL |
|
|
247 55 CONTINUE |
|
|
248 STR = AA*CSGNR - BB*CSGNI |
|
|
249 STI = AA*CSGNI + BB*CSGNR |
|
|
250 CYR(I) = STR*ATOL |
|
|
251 CYI(I) = STI*ATOL |
|
|
252 CSGNR = -CSGNR |
|
|
253 CSGNI = -CSGNI |
|
|
254 50 CONTINUE |
|
|
255 RETURN |
|
|
256 120 CONTINUE |
|
|
257 IF(NZ.EQ.(-2)) GO TO 130 |
|
|
258 NZ = 0 |
|
|
259 IERR=2 |
|
|
260 RETURN |
|
|
261 130 CONTINUE |
|
|
262 NZ=0 |
|
|
263 IERR=5 |
|
|
264 RETURN |
|
|
265 260 CONTINUE |
|
|
266 NZ=0 |
|
|
267 IERR=4 |
|
|
268 RETURN |
|
|
269 END |
