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