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