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