1      SUBROUTINE CMLRI(Z, FNU, KODE, N, Y, NZ, TOL)
2C***BEGIN PROLOGUE  CMLRI
3C***REFER TO  CBESI,CBESK
4C
5C     CMLRI 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  GAMLN,R1MACH
9C***END PROLOGUE  CMLRI
10      COMPLEX CK, CNORM, CONE, CTWO, CZERO, PT, P1, P2, RZ, SUM, Y, Z
11      REAL ACK, AK, AP, AT, AZ, BK, FKAP, FKK, FLAM, FNF, FNU, RHO,
12     * RHO2, SCLE, TFNF, TOL, TST, X, GAMLN, R1MACH
13      INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N
14      DIMENSION Y(N)
15      DATA CZERO,CONE,CTWO /(0.0E0,0.0E0),(1.0E0,0.0E0),(2.0E0,0.0E0)/
16      SCLE = 1.0E+3*R1MACH(1)/TOL
17      NZ=0
18      AZ = CABS(Z)
19      X = REAL(Z)
20      IAZ = INT(AZ)
21      IFNU = INT(FNU)
22      INU = IFNU + N - 1
23      AT = FLOAT(IAZ) + 1.0E0
24      CK = CMPLX(AT,0.0E0)/Z
25      RZ = CTWO/Z
26      P1 = CZERO
27      P2 = CONE
28      ACK = (AT+1.0E0)/AZ
29      RHO = ACK + SQRT(ACK*ACK-1.0E0)
30      RHO2 = RHO*RHO
31      TST = (RHO2+RHO2)/((RHO2-1.0E0)*(RHO-1.0E0))
32      TST = TST/TOL
33C-----------------------------------------------------------------------
34C     COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
35C-----------------------------------------------------------------------
36      AK = AT
37      DO 10 I=1,80
38        PT = P2
39        P2 = P1 - CK*P2
40        P1 = PT
41        CK = CK + RZ
42        AP = CABS(P2)
43        IF (AP.GT.TST*AK*AK) GO TO 20
44        AK = AK + 1.0E0
45   10 CONTINUE
46      GO TO 110
47   20 CONTINUE
48      I = I + 1
49      K = 0
50      IF (INU.LT.IAZ) GO TO 40
51C-----------------------------------------------------------------------
52C     COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
53C-----------------------------------------------------------------------
54      P1 = CZERO
55      P2 = CONE
56      AT = FLOAT(INU) + 1.0E0
57      CK = CMPLX(AT,0.0E0)/Z
58      ACK = AT/AZ
59      TST = SQRT(ACK/TOL)
60      ITIME = 1
61      DO 30 K=1,80
62        PT = P2
63        P2 = P1 - CK*P2
64        P1 = PT
65        CK = CK + RZ
66        AP = CABS(P2)
67        IF (AP.LT.TST) GO TO 30
68        IF (ITIME.EQ.2) GO TO 40
69        ACK = CABS(CK)
70        FLAM = ACK + SQRT(ACK*ACK-1.0E0)
71        FKAP = AP/CABS(P1)
72        RHO = AMIN1(FLAM,FKAP)
73        TST = TST*SQRT(RHO/(RHO*RHO-1.0E0))
74        ITIME = 2
75   30 CONTINUE
76      GO TO 110
77   40 CONTINUE
78C-----------------------------------------------------------------------
79C     BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
80C-----------------------------------------------------------------------
81      K = K + 1
82      KK = MAX0(I+IAZ,K+INU)
83      FKK = FLOAT(KK)
84      P1 = CZERO
85C-----------------------------------------------------------------------
86C     SCALE P2 AND SUM BY SCLE
87C-----------------------------------------------------------------------
88      P2 = CMPLX(SCLE,0.0E0)
89      FNF = FNU - FLOAT(IFNU)
90      TFNF = FNF + FNF
91      BK = GAMLN(FKK+TFNF+1.0E0,IDUM) - GAMLN(FKK+1.0E0,IDUM)
92     *     -GAMLN(TFNF+1.0E0,IDUM)
93      BK = EXP(BK)
94      SUM = CZERO
95      KM = KK - INU
96      DO 50 I=1,KM
97        PT = P2
98        P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2
99        P1 = PT
100        AK = 1.0E0 - TFNF/(FKK+TFNF)
101        ACK = BK*AK
102        SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1
103        BK = ACK
104        FKK = FKK - 1.0E0
105   50 CONTINUE
106      Y(N) = P2
107      IF (N.EQ.1) GO TO 70
108      DO 60 I=2,N
109        PT = P2
110        P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2
111        P1 = PT
112        AK = 1.0E0 - TFNF/(FKK+TFNF)
113        ACK = BK*AK
114        SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1
115        BK = ACK
116        FKK = FKK - 1.0E0
117        M = N - I + 1
118        Y(M) = P2
119   60 CONTINUE
120   70 CONTINUE
121      IF (IFNU.LE.0) GO TO 90
122      DO 80 I=1,IFNU
123        PT = P2
124        P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2
125        P1 = PT
126        AK = 1.0E0 - TFNF/(FKK+TFNF)
127        ACK = BK*AK
128        SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1
129        BK = ACK
130        FKK = FKK - 1.0E0
131   80 CONTINUE
132   90 CONTINUE
133      PT = Z
134      IF (KODE.EQ.2) PT = PT - CMPLX(X,0.0E0)
135      P1 = -CMPLX(FNF,0.0E0)*CLOG(RZ) + PT
136      AP = GAMLN(1.0E0+FNF,IDUM)
137      PT = P1 - CMPLX(AP,0.0E0)
138C-----------------------------------------------------------------------
139C     THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
140C     IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
141C-----------------------------------------------------------------------
142      P2 = P2 + SUM
143      AP = CABS(P2)
144      P1 = CMPLX(1.0E0/AP,0.0E0)
145      CK = CEXP(PT)*P1
146      PT = CONJG(P2)*P1
147      CNORM = CK*PT
148      DO 100 I=1,N
149        Y(I) = Y(I)*CNORM
150  100 CONTINUE
151      RETURN
152  110 CONTINUE
153      NZ=-2
154      RETURN
155      END
156