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
|
*DECK ZMLRI
SUBROUTINE ZMLRI (ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)
C***BEGIN PROLOGUE ZMLRI
C***SUBSIDIARY
C***PURPOSE Subsidiary to ZBESI and ZBESK
C***LIBRARY SLATEC
C***TYPE ALL (CMLRI-A, ZMLRI-A)
C***AUTHOR Amos, D. E., (SNL)
C***DESCRIPTION
C
C ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
C MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
C
C***SEE ALSO ZBESI, ZBESK
C***ROUTINES CALLED D1MACH, DGAMLN, ZABS, ZEXP, ZLOG, ZMLT
C***REVISION HISTORY (YYMMDD)
C 830501 DATE WRITTEN
C 910415 Prologue converted to Version 4.0 format. (BAB)
C 930122 Added ZEXP and ZLOG to EXTERNAL statement. (RWC)
C***END PROLOGUE ZMLRI
C COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z
DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI,
* CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, FNU, PTI, PTR, P1I,
* P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI,
* SUMR, TFNF, TOL, TST, YI, YR, ZEROI, ZEROR, ZI, ZR, DGAMLN,
* D1MACH, ZABS
INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ
DIMENSION YR(N), YI(N)
EXTERNAL ZABS, ZEXP, ZLOG
DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
C***FIRST EXECUTABLE STATEMENT ZMLRI
SCLE = D1MACH(1)/TOL
NZ=0
AZ = ZABS(ZR,ZI)
IAZ = AZ
IFNU = FNU
INU = IFNU + N - 1
AT = IAZ + 1.0D0
RAZ = 1.0D0/AZ
STR = ZR*RAZ
STI = -ZI*RAZ
CKR = STR*AT*RAZ
CKI = STI*AT*RAZ
RZR = (STR+STR)*RAZ
RZI = (STI+STI)*RAZ
P1R = ZEROR
P1I = ZEROI
P2R = CONER
P2I = CONEI
ACK = (AT+1.0D0)*RAZ
RHO = ACK + SQRT(ACK*ACK-1.0D0)
RHO2 = RHO*RHO
TST = (RHO2+RHO2)/((RHO2-1.0D0)*(RHO-1.0D0))
TST = TST/TOL
C-----------------------------------------------------------------------
C COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
C-----------------------------------------------------------------------
AK = AT
DO 10 I=1,80
PTR = P2R
PTI = P2I
P2R = P1R - (CKR*PTR-CKI*PTI)
P2I = P1I - (CKI*PTR+CKR*PTI)
P1R = PTR
P1I = PTI
CKR = CKR + RZR
CKI = CKI + RZI
AP = ZABS(P2R,P2I)
IF (AP.GT.TST*AK*AK) GO TO 20
AK = AK + 1.0D0
10 CONTINUE
GO TO 110
20 CONTINUE
I = I + 1
K = 0
IF (INU.LT.IAZ) GO TO 40
C-----------------------------------------------------------------------
C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
C-----------------------------------------------------------------------
P1R = ZEROR
P1I = ZEROI
P2R = CONER
P2I = CONEI
AT = INU + 1.0D0
STR = ZR*RAZ
STI = -ZI*RAZ
CKR = STR*AT*RAZ
CKI = STI*AT*RAZ
ACK = AT*RAZ
TST = SQRT(ACK/TOL)
ITIME = 1
DO 30 K=1,80
PTR = P2R
PTI = P2I
P2R = P1R - (CKR*PTR-CKI*PTI)
P2I = P1I - (CKR*PTI+CKI*PTR)
P1R = PTR
P1I = PTI
CKR = CKR + RZR
CKI = CKI + RZI
AP = ZABS(P2R,P2I)
IF (AP.LT.TST) GO TO 30
IF (ITIME.EQ.2) GO TO 40
ACK = ZABS(CKR,CKI)
FLAM = ACK + SQRT(ACK*ACK-1.0D0)
FKAP = AP/ZABS(P1R,P1I)
RHO = MIN(FLAM,FKAP)
TST = TST*SQRT(RHO/(RHO*RHO-1.0D0))
ITIME = 2
30 CONTINUE
GO TO 110
40 CONTINUE
C-----------------------------------------------------------------------
C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
C-----------------------------------------------------------------------
K = K + 1
KK = MAX(I+IAZ,K+INU)
FKK = KK
P1R = ZEROR
P1I = ZEROI
C-----------------------------------------------------------------------
C SCALE P2 AND SUM BY SCLE
C-----------------------------------------------------------------------
P2R = SCLE
P2I = ZEROI
FNF = FNU - IFNU
TFNF = FNF + FNF
BK = DGAMLN(FKK+TFNF+1.0D0,IDUM) - DGAMLN(FKK+1.0D0,IDUM) -
* DGAMLN(TFNF+1.0D0,IDUM)
BK = EXP(BK)
SUMR = ZEROR
SUMI = ZEROI
KM = KK - INU
DO 50 I=1,KM
PTR = P2R
PTI = P2I
P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
P1R = PTR
P1I = PTI
AK = 1.0D0 - TFNF/(FKK+TFNF)
ACK = BK*AK
SUMR = SUMR + (ACK+BK)*P1R
SUMI = SUMI + (ACK+BK)*P1I
BK = ACK
FKK = FKK - 1.0D0
50 CONTINUE
YR(N) = P2R
YI(N) = P2I
IF (N.EQ.1) GO TO 70
DO 60 I=2,N
PTR = P2R
PTI = P2I
P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
P1R = PTR
P1I = PTI
AK = 1.0D0 - TFNF/(FKK+TFNF)
ACK = BK*AK
SUMR = SUMR + (ACK+BK)*P1R
SUMI = SUMI + (ACK+BK)*P1I
BK = ACK
FKK = FKK - 1.0D0
M = N - I + 1
YR(M) = P2R
YI(M) = P2I
60 CONTINUE
70 CONTINUE
IF (IFNU.LE.0) GO TO 90
DO 80 I=1,IFNU
PTR = P2R
PTI = P2I
P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR)
P1R = PTR
P1I = PTI
AK = 1.0D0 - TFNF/(FKK+TFNF)
ACK = BK*AK
SUMR = SUMR + (ACK+BK)*P1R
SUMI = SUMI + (ACK+BK)*P1I
BK = ACK
FKK = FKK - 1.0D0
80 CONTINUE
90 CONTINUE
PTR = ZR
PTI = ZI
IF (KODE.EQ.2) PTR = ZEROR
CALL ZLOG(RZR, RZI, STR, STI, IDUM)
P1R = -FNF*STR + PTR
P1I = -FNF*STI + PTI
AP = DGAMLN(1.0D0+FNF,IDUM)
PTR = P1R - AP
PTI = P1I
C-----------------------------------------------------------------------
C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
C-----------------------------------------------------------------------
P2R = P2R + SUMR
P2I = P2I + SUMI
AP = ZABS(P2R,P2I)
P1R = 1.0D0/AP
CALL ZEXP(PTR, PTI, STR, STI)
CKR = STR*P1R
CKI = STI*P1R
PTR = P2R*P1R
PTI = -P2I*P1R
CALL ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI)
DO 100 I=1,N
STR = YR(I)*CNORMR - YI(I)*CNORMI
YI(I) = YR(I)*CNORMI + YI(I)*CNORMR
YR(I) = STR
100 CONTINUE
RETURN
110 CONTINUE
NZ=-2
RETURN
END
|
