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