1      SUBROUTINE ZMLRI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)
2C***BEGIN PROLOGUE  ZMLRI
3C***REFER TO  ZBESI,ZBESK
4C
5C     ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
6C     MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
7C
8C***ROUTINES CALLED  DGAMLN,D1MACH,ZABS,ZEXP,ZLOG,ZMLT
9C***END PROLOGUE  ZMLRI
10C     COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z
11      DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI,
12     * CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, FNU, PTI, PTR, P1I,
13     * P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI,
14     * SUMR, TFNF, TOL, TST, YI, YR, ZEROI, ZEROR, ZI, ZR, DGAMLN,
15     * D1MACH, ZABS
16      INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ
17      DIMENSION YR(N), YI(N)
18      DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
19      SCLE = D1MACH(1)/TOL
20      NZ=0
21      AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0)))
22      IAZ = INT(SNGL(AZ))
23      IFNU = INT(SNGL(FNU))
24      INU = IFNU + N - 1
25      AT = DBLE(FLOAT(IAZ)) + 1.0D0
26      RAZ = 1.0D0/AZ
27      STR = ZR*RAZ
28      STI = -ZI*RAZ
29      CKR = STR*AT*RAZ
30      CKI = STI*AT*RAZ
31      RZR = (STR+STR)*RAZ
32      RZI = (STI+STI)*RAZ
33      P1R = ZEROR
34      P1I = ZEROI
35      P2R = CONER
36      P2I = CONEI
37      ACK = (AT+1.0D0)*RAZ
38      RHO = ACK + DSQRT(ACK*ACK-1.0D0)
39      RHO2 = RHO*RHO
40      TST = (RHO2+RHO2)/((RHO2-1.0D0)*(RHO-1.0D0))
41      TST = TST/TOL
42C-----------------------------------------------------------------------
43C     COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
44C-----------------------------------------------------------------------
45      AK = AT
46      DO 10 I=1,80
47        PTR = P2R
48        PTI = P2I
49        P2R = P1R - (CKR*PTR-CKI*PTI)
50        P2I = P1I - (CKI*PTR+CKR*PTI)
51        P1R = PTR
52        P1I = PTI
53        CKR = CKR + RZR
54        CKI = CKI + RZI
55        AP = ZABS(CMPLX(P2R,P2I,kind=KIND(1.0D0)))
56        IF (AP.GT.TST*AK*AK) THEN
57          GO TO 20
58        END IF
59        AK = AK + 1.0D0
60   10 CONTINUE
61      GO TO 110
62   20 CONTINUE
63      I = I + 1
64      K = 0
65      IF (INU.LT.IAZ) GO TO 40
66C-----------------------------------------------------------------------
67C     COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
68C-----------------------------------------------------------------------
69      P1R = ZEROR
70      P1I = ZEROI
71      P2R = CONER
72      P2I = CONEI
73      AT = DBLE(FLOAT(INU)) + 1.0D0
74      STR = ZR*RAZ
75      STI = -ZI*RAZ
76      CKR = STR*AT*RAZ
77      CKI = STI*AT*RAZ
78      ACK = AT*RAZ
79      TST = DSQRT(ACK/TOL)
80      ITIME = 1
81      DO 30 K=1,80
82        PTR = P2R
83        PTI = P2I
84        P2R = P1R - (CKR*PTR-CKI*PTI)
85        P2I = P1I - (CKR*PTI+CKI*PTR)
86        P1R = PTR
87        P1I = PTI
88        CKR = CKR + RZR
89        CKI = CKI + RZI
90        AP = ZABS(CMPLX(P2R,P2I,kind=KIND(1.0D0)))
91        IF (AP.LT.TST) GO TO 30
92        IF (ITIME.EQ.2) GO TO 40
93        ACK = ZABS(CMPLX(CKR,CKI,kind=KIND(1.0D0)))
94        FLAM = ACK + DSQRT(ACK*ACK-1.0D0)
95        FKAP = AP/ZABS(CMPLX(P1R,P1I,kind=KIND(1.0D0)))
96        RHO = DMIN1(FLAM,FKAP)
97        TST = TST*DSQRT(RHO/(RHO*RHO-1.0D0))
98        ITIME = 2
99   30 CONTINUE
100      GO TO 110
101   40 CONTINUE
102C-----------------------------------------------------------------------
103C     BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
104C-----------------------------------------------------------------------
105      K = K + 1
106      KK = MAX0(I+IAZ,K+INU)
107      FKK = DBLE(FLOAT(KK))
108      P1R = ZEROR
109      P1I = ZEROI
110C-----------------------------------------------------------------------
111C     SCALE P2 AND SUM BY SCLE
112C-----------------------------------------------------------------------
113      P2R = SCLE
114      P2I = ZEROI
115      FNF = FNU - DBLE(FLOAT(IFNU))
116      TFNF = FNF + FNF
117      BK = DGAMLN(FKK+TFNF+1.0D0,IDUM) - DGAMLN(FKK+1.0D0,IDUM) -
118     * DGAMLN(TFNF+1.0D0,IDUM)
119      BK = DEXP(BK)
120      SUMR = ZEROR
121      SUMI = ZEROI
122      KM = KK - INU
123      DO 50 I=1,KM
124        PTR = P2R
125        PTI = P2I
126        P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
127        P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
128        P1R = PTR
129        P1I = PTI
130        AK = 1.0D0 - TFNF/(FKK+TFNF)
131        ACK = BK*AK
132        SUMR = SUMR + (ACK+BK)*P1R
133        SUMI = SUMI + (ACK+BK)*P1I
134        BK = ACK
135        FKK = FKK - 1.0D0
136   50 CONTINUE
137      YR(N) = P2R
138      YI(N) = P2I
139      IF (N.EQ.1) GO TO 70
140      DO 60 I=2,N
141        PTR = P2R
142        PTI = P2I
143        P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
144        P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
145        P1R = PTR
146        P1I = PTI
147        AK = 1.0D0 - TFNF/(FKK+TFNF)
148        ACK = BK*AK
149        SUMR = SUMR + (ACK+BK)*P1R
150        SUMI = SUMI + (ACK+BK)*P1I
151        BK = ACK
152        FKK = FKK - 1.0D0
153        M = N - I + 1
154        YR(M) = P2R
155        YI(M) = P2I
156   60 CONTINUE
157   70 CONTINUE
158      IF (IFNU.LE.0) GO TO 90
159      DO 80 I=1,IFNU
160        PTR = P2R
161        PTI = P2I
162        P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
163        P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR)
164        P1R = PTR
165        P1I = PTI
166        AK = 1.0D0 - TFNF/(FKK+TFNF)
167        ACK = BK*AK
168        SUMR = SUMR + (ACK+BK)*P1R
169        SUMI = SUMI + (ACK+BK)*P1I
170        BK = ACK
171        FKK = FKK - 1.0D0
172   80 CONTINUE
173   90 CONTINUE
174      PTR = ZR
175      PTI = ZI
176      IF (KODE.EQ.2) PTR = ZEROR
177      CALL ZLOG(RZR, RZI, STR, STI, IDUM)
178      P1R = -FNF*STR + PTR
179      P1I = -FNF*STI + PTI
180      AP = DGAMLN(1.0D0+FNF,IDUM)
181      PTR = P1R - AP
182      PTI = P1I
183C-----------------------------------------------------------------------
184C     THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
185C     IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
186C-----------------------------------------------------------------------
187      P2R = P2R + SUMR
188      P2I = P2I + SUMI
189      AP = ZABS(CMPLX(P2R,P2I,kind=KIND(1.0D0)))
190      P1R = 1.0D0/AP
191      CALL ZEXP(PTR, PTI, STR, STI)
192      CKR = STR*P1R
193      CKI = STI*P1R
194      PTR = P2R*P1R
195      PTI = -P2I*P1R
196      CALL ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI)
197      DO 100 I=1,N
198        STR = YR(I)*CNORMR - YI(I)*CNORMI
199        YI(I) = YR(I)*CNORMI + YI(I)*CNORMR
200        YR(I) = STR
201  100 CONTINUE
202      RETURN
203  110 CONTINUE
204      NZ=-2
205      RETURN
206      END
207