1      SUBROUTINE ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
2     * TOL, ELIM, ALIM)
3C***BEGIN PROLOGUE  ZUNI2
4C***REFER TO  ZBESI,ZBESK
5C
6C     ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF
7C     UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I
8C     OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.
9C
10C     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
11C     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
12C     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
13C     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
14C     Y(I)=CZERO FOR I=NLAST+1,N
15C
16C***ROUTINES CALLED  ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,XZABS
17C***END PROLOGUE  ZUNI2
18C     COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS,
19C    *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN
20      DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI,
21     * ARGR, ASCLE, ASUMI, ASUMR, BRY, BSUMI, BSUMR, CIDI, CIPI, CIPR,
22     * CONER, CRSC, CSCL, CSRR, CSSR, C1R, C2I, C2M, C2R, DAII,
23     * DAIR, ELIM, FN, FNU, FNUL, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI,
24     * RZR, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZBI, ZBR, ZEROI,
25     * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, ZNI, ZNR, ZR, CYR,
26     * CYI, D1MACH, XZABS, CAR, SAR
27      INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
28     * NN, NUF, NW, NZ, IDUM
29      DIMENSION BRY(3), YR(N), YI(N), CIPR(4), CIPI(4), CSSR(3),
30     * CSRR(3), CYR(2), CYI(2)
31      DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
32      DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4),
33     * CIPI(4)/ 1.0D0,0.0D0, 0.0D0,1.0D0, -1.0D0,0.0D0, 0.0D0,-1.0D0/
34      DATA HPI, AIC  /
35     1      1.57079632679489662D+00,     1.265512123484645396D+00/
36C
37      NZ = 0
38      ND = N
39      NLAST = 0
40C-----------------------------------------------------------------------
41C     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
42C     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
43C     EXP(ALIM)=EXP(ELIM)*TOL
44C-----------------------------------------------------------------------
45      CSCL = 1.0D0/TOL
46      CRSC = TOL
47      CSSR(1) = CSCL
48      CSSR(2) = CONER
49      CSSR(3) = CRSC
50      CSRR(1) = CRSC
51      CSRR(2) = CONER
52      CSRR(3) = CSCL
53      BRY(1) = 1.0D+3*D1MACH(1)/TOL
54C-----------------------------------------------------------------------
55C     ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
56C-----------------------------------------------------------------------
57      ZNR = ZI
58      ZNI = -ZR
59      ZBR = ZR
60      ZBI = ZI
61      CIDI = -CONER
62      INU = INT(SNGL(FNU))
63      ANG = HPI*(FNU-DBLE(FLOAT(INU)))
64      C2R = DCOS(ANG)
65      C2I = DSIN(ANG)
66      CAR = C2R
67      SAR = C2I
68      IN = INU + N - 1
69      IN = MOD(IN,4) + 1
70      STR = C2R*CIPR(IN) - C2I*CIPI(IN)
71      C2I = C2R*CIPI(IN) + C2I*CIPR(IN)
72      C2R = STR
73      IF (ZI.GT.0.0D0) GO TO 10
74      ZNR = -ZNR
75      ZBI = -ZBI
76      CIDI = -CIDI
77      C2I = -C2I
78   10 CONTINUE
79C-----------------------------------------------------------------------
80C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
81C-----------------------------------------------------------------------
82      FN = DMAX1(FNU,1.0D0)
83      CALL ZUNHJ(ZNR, ZNI, FN, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R,
84     * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
85      IF (KODE.EQ.1) GO TO 20
86      STR = ZBR + ZETA2R
87      STI = ZBI + ZETA2I
88      RAST = FN/XZABS(STR,STI)
89      STR = STR*RAST*RAST
90      STI = -STI*RAST*RAST
91      S1R = -ZETA1R + STR
92      S1I = -ZETA1I + STI
93      GO TO 30
94   20 CONTINUE
95      S1R = -ZETA1R + ZETA2R
96      S1I = -ZETA1I + ZETA2I
97   30 CONTINUE
98      RS1 = S1R
99      IF (DABS(RS1).GT.ELIM) GO TO 150
100   40 CONTINUE
101      NN = MIN0(2,ND)
102      DO 90 I=1,NN
103        FN = FNU + DBLE(FLOAT(ND-I))
104        CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR, PHII, ARGR, ARGI,
105     *   ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
106        IF (KODE.EQ.1) GO TO 50
107        STR = ZBR + ZETA2R
108        STI = ZBI + ZETA2I
109        RAST = FN/XZABS(STR,STI)
110        STR = STR*RAST*RAST
111        STI = -STI*RAST*RAST
112        S1R = -ZETA1R + STR
113        S1I = -ZETA1I + STI + DABS(ZI)
114        GO TO 60
115   50   CONTINUE
116        S1R = -ZETA1R + ZETA2R
117        S1I = -ZETA1I + ZETA2I
118   60   CONTINUE
119C-----------------------------------------------------------------------
120C     TEST FOR UNDERFLOW AND OVERFLOW
121C-----------------------------------------------------------------------
122        RS1 = S1R
123        IF (DABS(RS1).GT.ELIM) GO TO 120
124        IF (I.EQ.1) IFLAG = 2
125        IF (DABS(RS1).LT.ALIM) GO TO 70
126C-----------------------------------------------------------------------
127C     REFINE  TEST AND SCALE
128C-----------------------------------------------------------------------
129C-----------------------------------------------------------------------
130        APHI = XZABS(PHIR,PHII)
131        AARG = XZABS(ARGR,ARGI)
132        RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC
133        IF (DABS(RS1).GT.ELIM) GO TO 120
134        IF (I.EQ.1) IFLAG = 1
135        IF (RS1.LT.0.0D0) GO TO 70
136        IF (I.EQ.1) IFLAG = 3
137   70   CONTINUE
138C-----------------------------------------------------------------------
139C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
140C     EXPONENT EXTREMES
141C-----------------------------------------------------------------------
142        CALL ZAIRY(ARGR, ARGI, 0, 2, AIR, AII, NAI, IDUM)
143        CALL ZAIRY(ARGR, ARGI, 1, 2, DAIR, DAII, NDAI, IDUM)
144        STR = DAIR*BSUMR - DAII*BSUMI
145        STI = DAIR*BSUMI + DAII*BSUMR
146        STR = STR + (AIR*ASUMR-AII*ASUMI)
147        STI = STI + (AIR*ASUMI+AII*ASUMR)
148        S2R = PHIR*STR - PHII*STI
149        S2I = PHIR*STI + PHII*STR
150        STR = DEXP(S1R)*CSSR(IFLAG)
151        S1R = STR*DCOS(S1I)
152        S1I = STR*DSIN(S1I)
153        STR = S2R*S1R - S2I*S1I
154        S2I = S2R*S1I + S2I*S1R
155        S2R = STR
156        IF (IFLAG.NE.1) GO TO 80
157        CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
158        IF (NW.NE.0) GO TO 120
159   80   CONTINUE
160        IF (ZI.LE.0.0D0) S2I = -S2I
161        STR = S2R*C2R - S2I*C2I
162        S2I = S2R*C2I + S2I*C2R
163        S2R = STR
164        CYR(I) = S2R
165        CYI(I) = S2I
166        J = ND - I + 1
167        YR(J) = S2R*CSRR(IFLAG)
168        YI(J) = S2I*CSRR(IFLAG)
169        STR = -C2I*CIDI
170        C2I = C2R*CIDI
171        C2R = STR
172   90 CONTINUE
173      IF (ND.LE.2) GO TO 110
174      RAZ = 1.0D0/XZABS(ZR,ZI)
175      STR = ZR*RAZ
176      STI = -ZI*RAZ
177      RZR = (STR+STR)*RAZ
178      RZI = (STI+STI)*RAZ
179      BRY(2) = 1.0D0/BRY(1)
180      BRY(3) = D1MACH(2)
181      S1R = CYR(1)
182      S1I = CYI(1)
183      S2R = CYR(2)
184      S2I = CYI(2)
185      C1R = CSRR(IFLAG)
186      ASCLE = BRY(IFLAG)
187      K = ND - 2
188      FN = DBLE(FLOAT(K))
189      DO 100 I=3,ND
190        C2R = S2R
191        C2I = S2I
192        S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I)
193        S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R)
194        S1R = C2R
195        S1I = C2I
196        C2R = S2R*C1R
197        C2I = S2I*C1R
198        YR(K) = C2R
199        YI(K) = C2I
200        K = K - 1
201        FN = FN - 1.0D0
202        IF (IFLAG.GE.3) GO TO 100
203        STR = DABS(C2R)
204        STI = DABS(C2I)
205        C2M = DMAX1(STR,STI)
206        IF (C2M.LE.ASCLE) GO TO 100
207        IFLAG = IFLAG + 1
208        ASCLE = BRY(IFLAG)
209        S1R = S1R*C1R
210        S1I = S1I*C1R
211        S2R = C2R
212        S2I = C2I
213        S1R = S1R*CSSR(IFLAG)
214        S1I = S1I*CSSR(IFLAG)
215        S2R = S2R*CSSR(IFLAG)
216        S2I = S2I*CSSR(IFLAG)
217        C1R = CSRR(IFLAG)
218  100 CONTINUE
219  110 CONTINUE
220      RETURN
221  120 CONTINUE
222      IF (RS1.GT.0.0D0) GO TO 140
223C-----------------------------------------------------------------------
224C     SET UNDERFLOW AND UPDATE PARAMETERS
225C-----------------------------------------------------------------------
226      YR(ND) = ZEROR
227      YI(ND) = ZEROI
228      NZ = NZ + 1
229      ND = ND - 1
230      IF (ND.EQ.0) GO TO 110
231      CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM)
232      IF (NUF.LT.0) GO TO 140
233      ND = ND - NUF
234      NZ = NZ + NUF
235      IF (ND.EQ.0) GO TO 110
236      FN = FNU + DBLE(FLOAT(ND-1))
237      IF (FN.LT.FNUL) GO TO 130
238C      FN = CIDI
239C      J = NUF + 1
240C      K = MOD(J,4) + 1
241C      S1R = CIPR(K)
242C      S1I = CIPI(K)
243C      IF (FN.LT.0.0D0) S1I = -S1I
244C      STR = C2R*S1R - C2I*S1I
245C      C2I = C2R*S1I + C2I*S1R
246C      C2R = STR
247      IN = INU + ND - 1
248      IN = MOD(IN,4) + 1
249      C2R = CAR*CIPR(IN) - SAR*CIPI(IN)
250      C2I = CAR*CIPI(IN) + SAR*CIPR(IN)
251      IF (ZI.LE.0.0D0) C2I = -C2I
252      GO TO 40
253  130 CONTINUE
254      NLAST = ND
255      RETURN
256  140 CONTINUE
257      NZ = -1
258      RETURN
259  150 CONTINUE
260      IF (RS1.GT.0.0D0) GO TO 140
261      NZ = N
262      DO 160 I=1,N
263        YR(I) = ZEROR
264        YI(I) = ZEROI
265  160 CONTINUE
266      RETURN
267      END
268