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