|
3217
|
1 SUBROUTINE ZSERI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, |
|
|
2 * ALIM) |
|
|
3 C***BEGIN PROLOGUE ZSERI |
|
|
4 C***REFER TO ZBESI,ZBESK |
|
|
5 C |
|
|
6 C ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY |
|
|
7 C MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE |
|
|
8 C REGION CABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN. |
|
|
9 C NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO |
|
|
10 C DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE |
|
|
11 C CONDITION CABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE |
|
|
12 C COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ). |
|
|
13 C |
|
|
14 C***ROUTINES CALLED DGAMLN,D1MACH,ZUCHK,XZABS,ZDIV,XZLOG,ZMLT |
|
|
15 C***END PROLOGUE ZSERI |
|
|
16 C COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z |
|
|
17 DOUBLE PRECISION AA, ACZ, AK, AK1I, AK1R, ALIM, ARM, ASCLE, ATOL, |
|
|
18 * AZ, CKI, CKR, COEFI, COEFR, CONEI, CONER, CRSCR, CZI, CZR, DFNU, |
|
|
19 * ELIM, FNU, FNUP, HZI, HZR, RAZ, RS, RTR1, RZI, RZR, S, SS, STI, |
|
|
20 * STR, S1I, S1R, S2I, S2R, TOL, YI, YR, WI, WR, ZEROI, ZEROR, ZI, |
|
|
21 * ZR, DGAMLN, D1MACH, XZABS |
|
|
22 INTEGER I, IB, IDUM, IFLAG, IL, K, KODE, L, M, N, NN, NZ, NW |
|
|
23 DIMENSION YR(N), YI(N), WR(2), WI(2) |
|
|
24 DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / |
|
|
25 C |
|
|
26 NZ = 0 |
|
|
27 AZ = XZABS(ZR,ZI) |
|
|
28 IF (AZ.EQ.0.0D0) GO TO 160 |
|
|
29 ARM = 1.0D+3*D1MACH(1) |
|
|
30 RTR1 = DSQRT(ARM) |
|
|
31 CRSCR = 1.0D0 |
|
|
32 IFLAG = 0 |
|
|
33 IF (AZ.LT.ARM) GO TO 150 |
|
|
34 HZR = 0.5D0*ZR |
|
|
35 HZI = 0.5D0*ZI |
|
|
36 CZR = ZEROR |
|
|
37 CZI = ZEROI |
|
|
38 IF (AZ.LE.RTR1) GO TO 10 |
|
|
39 CALL ZMLT(HZR, HZI, HZR, HZI, CZR, CZI) |
|
|
40 10 CONTINUE |
|
|
41 ACZ = XZABS(CZR,CZI) |
|
|
42 NN = N |
|
|
43 CALL XZLOG(HZR, HZI, CKR, CKI, IDUM) |
|
|
44 20 CONTINUE |
|
|
45 DFNU = FNU + DBLE(FLOAT(NN-1)) |
|
|
46 FNUP = DFNU + 1.0D0 |
|
|
47 C----------------------------------------------------------------------- |
|
|
48 C UNDERFLOW TEST |
|
|
49 C----------------------------------------------------------------------- |
|
|
50 AK1R = CKR*DFNU |
|
|
51 AK1I = CKI*DFNU |
|
|
52 AK = DGAMLN(FNUP,IDUM) |
|
|
53 AK1R = AK1R - AK |
|
|
54 IF (KODE.EQ.2) AK1R = AK1R - ZR |
|
|
55 IF (AK1R.GT.(-ELIM)) GO TO 40 |
|
|
56 30 CONTINUE |
|
|
57 NZ = NZ + 1 |
|
|
58 YR(NN) = ZEROR |
|
|
59 YI(NN) = ZEROI |
|
|
60 IF (ACZ.GT.DFNU) GO TO 190 |
|
|
61 NN = NN - 1 |
|
|
62 IF (NN.EQ.0) RETURN |
|
|
63 GO TO 20 |
|
|
64 40 CONTINUE |
|
|
65 IF (AK1R.GT.(-ALIM)) GO TO 50 |
|
|
66 IFLAG = 1 |
|
|
67 SS = 1.0D0/TOL |
|
|
68 CRSCR = TOL |
|
|
69 ASCLE = ARM*SS |
|
|
70 50 CONTINUE |
|
|
71 AA = DEXP(AK1R) |
|
|
72 IF (IFLAG.EQ.1) AA = AA*SS |
|
|
73 COEFR = AA*DCOS(AK1I) |
|
|
74 COEFI = AA*DSIN(AK1I) |
|
|
75 ATOL = TOL*ACZ/FNUP |
|
|
76 IL = MIN0(2,NN) |
|
|
77 DO 90 I=1,IL |
|
|
78 DFNU = FNU + DBLE(FLOAT(NN-I)) |
|
|
79 FNUP = DFNU + 1.0D0 |
|
|
80 S1R = CONER |
|
|
81 S1I = CONEI |
|
|
82 IF (ACZ.LT.TOL*FNUP) GO TO 70 |
|
|
83 AK1R = CONER |
|
|
84 AK1I = CONEI |
|
|
85 AK = FNUP + 2.0D0 |
|
|
86 S = FNUP |
|
|
87 AA = 2.0D0 |
|
|
88 60 CONTINUE |
|
|
89 RS = 1.0D0/S |
|
|
90 STR = AK1R*CZR - AK1I*CZI |
|
|
91 STI = AK1R*CZI + AK1I*CZR |
|
|
92 AK1R = STR*RS |
|
|
93 AK1I = STI*RS |
|
|
94 S1R = S1R + AK1R |
|
|
95 S1I = S1I + AK1I |
|
|
96 S = S + AK |
|
|
97 AK = AK + 2.0D0 |
|
|
98 AA = AA*ACZ*RS |
|
|
99 IF (AA.GT.ATOL) GO TO 60 |
|
|
100 70 CONTINUE |
|
|
101 S2R = S1R*COEFR - S1I*COEFI |
|
|
102 S2I = S1R*COEFI + S1I*COEFR |
|
|
103 WR(I) = S2R |
|
|
104 WI(I) = S2I |
|
|
105 IF (IFLAG.EQ.0) GO TO 80 |
|
|
106 CALL ZUCHK(S2R, S2I, NW, ASCLE, TOL) |
|
|
107 IF (NW.NE.0) GO TO 30 |
|
|
108 80 CONTINUE |
|
|
109 M = NN - I + 1 |
|
|
110 YR(M) = S2R*CRSCR |
|
|
111 YI(M) = S2I*CRSCR |
|
|
112 IF (I.EQ.IL) GO TO 90 |
|
|
113 CALL ZDIV(COEFR, COEFI, HZR, HZI, STR, STI) |
|
|
114 COEFR = STR*DFNU |
|
|
115 COEFI = STI*DFNU |
|
|
116 90 CONTINUE |
|
|
117 IF (NN.LE.2) RETURN |
|
|
118 K = NN - 2 |
|
|
119 AK = DBLE(FLOAT(K)) |
|
|
120 RAZ = 1.0D0/AZ |
|
|
121 STR = ZR*RAZ |
|
|
122 STI = -ZI*RAZ |
|
|
123 RZR = (STR+STR)*RAZ |
|
|
124 RZI = (STI+STI)*RAZ |
|
|
125 IF (IFLAG.EQ.1) GO TO 120 |
|
|
126 IB = 3 |
|
|
127 100 CONTINUE |
|
|
128 DO 110 I=IB,NN |
|
|
129 YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2) |
|
|
130 YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2) |
|
|
131 AK = AK - 1.0D0 |
|
|
132 K = K - 1 |
|
|
133 110 CONTINUE |
|
|
134 RETURN |
|
|
135 C----------------------------------------------------------------------- |
|
|
136 C RECUR BACKWARD WITH SCALED VALUES |
|
|
137 C----------------------------------------------------------------------- |
|
|
138 120 CONTINUE |
|
|
139 C----------------------------------------------------------------------- |
|
|
140 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE |
|
|
141 C UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1.0D+3 |
|
|
142 C----------------------------------------------------------------------- |
|
|
143 S1R = WR(1) |
|
|
144 S1I = WI(1) |
|
|
145 S2R = WR(2) |
|
|
146 S2I = WI(2) |
|
|
147 DO 130 L=3,NN |
|
|
148 CKR = S2R |
|
|
149 CKI = S2I |
|
|
150 S2R = S1R + (AK+FNU)*(RZR*CKR-RZI*CKI) |
|
|
151 S2I = S1I + (AK+FNU)*(RZR*CKI+RZI*CKR) |
|
|
152 S1R = CKR |
|
|
153 S1I = CKI |
|
|
154 CKR = S2R*CRSCR |
|
|
155 CKI = S2I*CRSCR |
|
|
156 YR(K) = CKR |
|
|
157 YI(K) = CKI |
|
|
158 AK = AK - 1.0D0 |
|
|
159 K = K - 1 |
|
|
160 IF (XZABS(CKR,CKI).GT.ASCLE) GO TO 140 |
|
|
161 130 CONTINUE |
|
|
162 RETURN |
|
|
163 140 CONTINUE |
|
|
164 IB = L + 1 |
|
|
165 IF (IB.GT.NN) RETURN |
|
|
166 GO TO 100 |
|
|
167 150 CONTINUE |
|
|
168 NZ = N |
|
|
169 IF (FNU.EQ.0.0D0) NZ = NZ - 1 |
|
|
170 160 CONTINUE |
|
|
171 YR(1) = ZEROR |
|
|
172 YI(1) = ZEROI |
|
|
173 IF (FNU.NE.0.0D0) GO TO 170 |
|
|
174 YR(1) = CONER |
|
|
175 YI(1) = CONEI |
|
|
176 170 CONTINUE |
|
|
177 IF (N.EQ.1) RETURN |
|
|
178 DO 180 I=2,N |
|
|
179 YR(I) = ZEROR |
|
|
180 YI(I) = ZEROI |
|
|
181 180 CONTINUE |
|
|
182 RETURN |
|
|
183 C----------------------------------------------------------------------- |
|
|
184 C RETURN WITH NZ.LT.0 IF CABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE |
|
|
185 C THE CALCULATION IN CBINU WITH N=N-IABS(NZ) |
|
|
186 C----------------------------------------------------------------------- |
|
|
187 190 CONTINUE |
|
|
188 NZ = -NZ |
|
|
189 RETURN |
|
|
190 END |
