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
|
SUBROUTINE ZSERI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM,
* ALIM)
C***BEGIN PROLOGUE ZSERI
C***REFER TO ZBESI,ZBESK
C
C ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
C MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE
C REGION CABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN.
C NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO
C DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE
C CONDITION CABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE
C COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
C
C***ROUTINES CALLED DGAMLN,D1MACH,ZUCHK,AZABS,ZDIV,AZLOG,ZMLT
C***END PROLOGUE ZSERI
C COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z
DOUBLE PRECISION AA, ACZ, AK, AK1I, AK1R, ALIM, ARM, ASCLE, ATOL,
* AZ, CKI, CKR, COEFI, COEFR, CONEI, CONER, CRSCR, CZI, CZR, DFNU,
* ELIM, FNU, FNUP, HZI, HZR, RAZ, RS, RTR1, RZI, RZR, S, SS, STI,
* STR, S1I, S1R, S2I, S2R, TOL, YI, YR, WI, WR, ZEROI, ZEROR, ZI,
* ZR, DGAMLN, D1MACH, AZABS
INTEGER I, IB, IDUM, IFLAG, IL, K, KODE, L, M, N, NN, NZ, NW
DIMENSION YR(N), YI(N), WR(2), WI(2)
DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
C
NZ = 0
AZ = AZABS(ZR,ZI)
IF (AZ.EQ.0.0D0) GO TO 160
ARM = 1.0D+3*D1MACH(1)
RTR1 = DSQRT(ARM)
CRSCR = 1.0D0
IFLAG = 0
IF (AZ.LT.ARM) GO TO 150
HZR = 0.5D0*ZR
HZI = 0.5D0*ZI
CZR = ZEROR
CZI = ZEROI
IF (AZ.LE.RTR1) GO TO 10
CALL ZMLT(HZR, HZI, HZR, HZI, CZR, CZI)
10 CONTINUE
ACZ = AZABS(CZR,CZI)
NN = N
CALL AZLOG(HZR, HZI, CKR, CKI, IDUM)
20 CONTINUE
DFNU = FNU + DBLE(FLOAT(NN-1))
FNUP = DFNU + 1.0D0
C-----------------------------------------------------------------------
C UNDERFLOW TEST
C-----------------------------------------------------------------------
AK1R = CKR*DFNU
AK1I = CKI*DFNU
AK = DGAMLN(FNUP,IDUM)
AK1R = AK1R - AK
IF (KODE.EQ.2) AK1R = AK1R - ZR
IF (AK1R.GT.(-ELIM)) GO TO 40
30 CONTINUE
NZ = NZ + 1
YR(NN) = ZEROR
YI(NN) = ZEROI
IF (ACZ.GT.DFNU) GO TO 190
NN = NN - 1
IF (NN.EQ.0) RETURN
GO TO 20
40 CONTINUE
IF (AK1R.GT.(-ALIM)) GO TO 50
IFLAG = 1
SS = 1.0D0/TOL
CRSCR = TOL
ASCLE = ARM*SS
50 CONTINUE
AA = DEXP(AK1R)
IF (IFLAG.EQ.1) AA = AA*SS
COEFR = AA*DCOS(AK1I)
COEFI = AA*DSIN(AK1I)
ATOL = TOL*ACZ/FNUP
IL = MIN0(2,NN)
DO 90 I=1,IL
DFNU = FNU + DBLE(FLOAT(NN-I))
FNUP = DFNU + 1.0D0
S1R = CONER
S1I = CONEI
IF (ACZ.LT.TOL*FNUP) GO TO 70
AK1R = CONER
AK1I = CONEI
AK = FNUP + 2.0D0
S = FNUP
AA = 2.0D0
60 CONTINUE
RS = 1.0D0/S
STR = AK1R*CZR - AK1I*CZI
STI = AK1R*CZI + AK1I*CZR
AK1R = STR*RS
AK1I = STI*RS
S1R = S1R + AK1R
S1I = S1I + AK1I
S = S + AK
AK = AK + 2.0D0
AA = AA*ACZ*RS
IF (AA.GT.ATOL) GO TO 60
70 CONTINUE
S2R = S1R*COEFR - S1I*COEFI
S2I = S1R*COEFI + S1I*COEFR
WR(I) = S2R
WI(I) = S2I
IF (IFLAG.EQ.0) GO TO 80
CALL ZUCHK(S2R, S2I, NW, ASCLE, TOL)
IF (NW.NE.0) GO TO 30
80 CONTINUE
M = NN - I + 1
YR(M) = S2R*CRSCR
YI(M) = S2I*CRSCR
IF (I.EQ.IL) GO TO 90
CALL ZDIV(COEFR, COEFI, HZR, HZI, STR, STI)
COEFR = STR*DFNU
COEFI = STI*DFNU
90 CONTINUE
IF (NN.LE.2) RETURN
K = NN - 2
AK = DBLE(FLOAT(K))
RAZ = 1.0D0/AZ
STR = ZR*RAZ
STI = -ZI*RAZ
RZR = (STR+STR)*RAZ
RZI = (STI+STI)*RAZ
IF (IFLAG.EQ.1) GO TO 120
IB = 3
100 CONTINUE
DO 110 I=IB,NN
YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2)
YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2)
AK = AK - 1.0D0
K = K - 1
110 CONTINUE
RETURN
C-----------------------------------------------------------------------
C RECUR BACKWARD WITH SCALED VALUES
C-----------------------------------------------------------------------
120 CONTINUE
C-----------------------------------------------------------------------
C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
C UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1.0D+3
C-----------------------------------------------------------------------
S1R = WR(1)
S1I = WI(1)
S2R = WR(2)
S2I = WI(2)
DO 130 L=3,NN
CKR = S2R
CKI = S2I
S2R = S1R + (AK+FNU)*(RZR*CKR-RZI*CKI)
S2I = S1I + (AK+FNU)*(RZR*CKI+RZI*CKR)
S1R = CKR
S1I = CKI
CKR = S2R*CRSCR
CKI = S2I*CRSCR
YR(K) = CKR
YI(K) = CKI
AK = AK - 1.0D0
K = K - 1
IF (AZABS(CKR,CKI).GT.ASCLE) GO TO 140
130 CONTINUE
RETURN
140 CONTINUE
IB = L + 1
IF (IB.GT.NN) RETURN
GO TO 100
150 CONTINUE
NZ = N
IF (FNU.EQ.0.0D0) NZ = NZ - 1
160 CONTINUE
YR(1) = ZEROR
YI(1) = ZEROI
IF (FNU.NE.0.0D0) GO TO 170
YR(1) = CONER
YI(1) = CONEI
170 CONTINUE
IF (N.EQ.1) RETURN
DO 180 I=2,N
YR(I) = ZEROR
YI(I) = ZEROI
180 CONTINUE
RETURN
C-----------------------------------------------------------------------
C RETURN WITH NZ.LT.0 IF CABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE
C THE CALCULATION IN CBINU WITH N=N-IABS(NZ)
C-----------------------------------------------------------------------
190 CONTINUE
NZ = -NZ
RETURN
END
|
