1      SUBROUTINE CUNI2(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM,
2     * ALIM)
3C***BEGIN PROLOGUE  CUNI2
4C***REFER TO  CBESI,CBESK
5C
6C     CUNI2 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  CAIRY,CUCHK,CUNHJ,CUOIK,R1MACH
17C***END PROLOGUE  CUNI2
18      COMPLEX AI, ARG, ASUM, BSUM, CFN, CI, CID, CIP, CONE, CRSC, CSCL,
19     * CSR, CSS, CY, CZERO, C1, C2, DAI, PHI, RZ, S1, S2, Y, Z, ZB,
20     * ZETA1, ZETA2, ZN, ZAR
21      REAL AARG, AIC, ALIM, ANG, APHI, ASCLE, AY, BRY, CAR, C2I, C2M,
22     * C2R, ELIM, FN, FNU, FNUL, HPI, RS1, SAR, TOL, YY, R1MACH
23      INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
24     * NN, NUF, NW, NZ, IDUM
25      DIMENSION BRY(3), Y(N), CIP(4), CSS(3), CSR(3), CY(2)
26      DATA CZERO,CONE,CI/(0.0E0,0.0E0),(1.0E0,0.0E0),(0.0E0,1.0E0)/
27      DATA CIP(1),CIP(2),CIP(3),CIP(4)/
28     1 (1.0E0,0.0E0), (0.0E0,1.0E0), (-1.0E0,0.0E0), (0.0E0,-1.0E0)/
29      DATA HPI, AIC  /
30     1      1.57079632679489662E+00,     1.265512123484645396E+00/
31C
32      NZ = 0
33      ND = N
34      NLAST = 0
35C-----------------------------------------------------------------------
36C     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
37C     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
38C     EXP(ALIM)=EXP(ELIM)*TOL
39C-----------------------------------------------------------------------
40      CSCL = CMPLX(1.0E0/TOL,0.0E0)
41      CRSC = CMPLX(TOL,0.0E0)
42      CSS(1) = CSCL
43      CSS(2) = CONE
44      CSS(3) = CRSC
45      CSR(1) = CRSC
46      CSR(2) = CONE
47      CSR(3) = CSCL
48      BRY(1) = 1.0E+3*R1MACH(1)/TOL
49      YY = AIMAG(Z)
50C-----------------------------------------------------------------------
51C     ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
52C-----------------------------------------------------------------------
53      ZN = -Z*CI
54      ZB = Z
55      CID = -CI
56      INU = INT(FNU)
57      ANG = HPI*(FNU-FLOAT(INU))
58      CAR = COS(ANG)
59      SAR = SIN(ANG)
60      C2 = CMPLX(CAR,SAR)
61      ZAR = C2
62      IN = INU + N - 1
63      IN = MOD(IN,4)
64      C2 = C2*CIP(IN+1)
65      IF (YY.GT.0.0E0) GO TO 10
66      ZN = CONJG(-ZN)
67      ZB = CONJG(ZB)
68      CID = -CID
69      C2 = CONJG(C2)
70   10 CONTINUE
71C-----------------------------------------------------------------------
72C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
73C-----------------------------------------------------------------------
74      FN = AMAX1(FNU,1.0E0)
75      CALL CUNHJ(ZN, FN, 1, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM)
76      IF (KODE.EQ.1) GO TO 20
77      CFN = CMPLX(FNU,0.0E0)
78      S1 = -ZETA1 + CFN*(CFN/(ZB+ZETA2))
79      GO TO 30
80   20 CONTINUE
81      S1 = -ZETA1 + ZETA2
82   30 CONTINUE
83      RS1 = REAL(S1)
84      IF (ABS(RS1).GT.ELIM) GO TO 150
85   40 CONTINUE
86      NN = MIN0(2,ND)
87      DO 90 I=1,NN
88        FN = FNU + FLOAT(ND-I)
89        CALL CUNHJ(ZN, FN, 0, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM)
90        IF (KODE.EQ.1) GO TO 50
91        CFN = CMPLX(FN,0.0E0)
92        AY = ABS(YY)
93        S1 = -ZETA1 + CFN*(CFN/(ZB+ZETA2)) + CMPLX(0.0E0,AY)
94        GO TO 60
95   50   CONTINUE
96        S1 = -ZETA1 + ZETA2
97   60   CONTINUE
98C-----------------------------------------------------------------------
99C     TEST FOR UNDERFLOW AND OVERFLOW
100C-----------------------------------------------------------------------
101        RS1 = REAL(S1)
102        IF (ABS(RS1).GT.ELIM) GO TO 120
103        IF (I.EQ.1) IFLAG = 2
104        IF (ABS(RS1).LT.ALIM) GO TO 70
105C-----------------------------------------------------------------------
106C     REFINE  TEST AND SCALE
107C-----------------------------------------------------------------------
108C-----------------------------------------------------------------------
109        APHI = CABS(PHI)
110        AARG = CABS(ARG)
111        RS1 = RS1 + ALOG(APHI) - 0.25E0*ALOG(AARG) - AIC
112        IF (ABS(RS1).GT.ELIM) GO TO 120
113        IF (I.EQ.1) IFLAG = 1
114        IF (RS1.LT.0.0E0) GO TO 70
115        IF (I.EQ.1) IFLAG = 3
116   70   CONTINUE
117C-----------------------------------------------------------------------
118C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
119C     EXPONENT EXTREMES
120C-----------------------------------------------------------------------
121        CALL CAIRY(ARG, 0, 2, AI, NAI, IDUM)
122        CALL CAIRY(ARG, 1, 2, DAI, NDAI, IDUM)
123        S2 = PHI*(AI*ASUM+DAI*BSUM)
124        C2R = REAL(S1)
125        C2I = AIMAG(S1)
126        C2M = EXP(C2R)*REAL(CSS(IFLAG))
127        S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I))
128        S2 = S2*S1
129        IF (IFLAG.NE.1) GO TO 80
130        CALL CUCHK(S2, NW, BRY(1), TOL)
131        IF (NW.NE.0) GO TO 120
132   80   CONTINUE
133        IF (YY.LE.0.0E0) S2 = CONJG(S2)
134        J = ND - I + 1
135        S2 = S2*C2
136        CY(I) = S2
137        Y(J) = S2*CSR(IFLAG)
138        C2 = C2*CID
139   90 CONTINUE
140      IF (ND.LE.2) GO TO 110
141      RZ = CMPLX(2.0E0,0.0E0)/Z
142      BRY(2) = 1.0E0/BRY(1)
143      BRY(3) = R1MACH(2)
144      S1 = CY(1)
145      S2 = CY(2)
146      C1 = CSR(IFLAG)
147      ASCLE = BRY(IFLAG)
148      K = ND - 2
149      FN = FLOAT(K)
150      DO 100 I=3,ND
151        C2 = S2
152        S2 = S1 + CMPLX(FNU+FN,0.0E0)*RZ*S2
153        S1 = C2
154        C2 = S2*C1
155        Y(K) = C2
156        K = K - 1
157        FN = FN - 1.0E0
158        IF (IFLAG.GE.3) GO TO 100
159        C2R = REAL(C2)
160        C2I = AIMAG(C2)
161        C2R = ABS(C2R)
162        C2I = ABS(C2I)
163        C2M = AMAX1(C2R,C2I)
164        IF (C2M.LE.ASCLE) GO TO 100
165        IFLAG = IFLAG + 1
166        ASCLE = BRY(IFLAG)
167        S1 = S1*C1
168        S2 = C2
169        S1 = S1*CSS(IFLAG)
170        S2 = S2*CSS(IFLAG)
171        C1 = CSR(IFLAG)
172  100 CONTINUE
173  110 CONTINUE
174      RETURN
175  120 CONTINUE
176      IF (RS1.GT.0.0E0) GO TO 140
177C-----------------------------------------------------------------------
178C     SET UNDERFLOW AND UPDATE PARAMETERS
179C-----------------------------------------------------------------------
180      Y(ND) = CZERO
181      NZ = NZ + 1
182      ND = ND - 1
183      IF (ND.EQ.0) GO TO 110
184      CALL CUOIK(Z, FNU, KODE, 1, ND, Y, NUF, TOL, ELIM, ALIM)
185      IF (NUF.LT.0) GO TO 140
186      ND = ND - NUF
187      NZ = NZ + NUF
188      IF (ND.EQ.0) GO TO 110
189      FN = FNU + FLOAT(ND-1)
190      IF (FN.LT.FNUL) GO TO 130
191C      FN = AIMAG(CID)
192C      J = NUF + 1
193C      K = MOD(J,4) + 1
194C      S1 = CIP(K)
195C      IF (FN.LT.0.0E0) S1 = CONJG(S1)
196C      C2 = C2*S1
197      IN = INU + ND - 1
198      IN = MOD(IN,4) + 1
199      C2 = ZAR*CIP(IN)
200      IF (YY.LE.0.0E0)C2=CONJG(C2)
201      GO TO 40
202  130 CONTINUE
203      NLAST = ND
204      RETURN
205  140 CONTINUE
206      NZ = -1
207      RETURN
208  150 CONTINUE
209      IF (RS1.GT.0.0E0) GO TO 140
210      NZ = N
211      DO 160 I=1,N
212        Y(I) = CZERO
213  160 CONTINUE
214      RETURN
215      END
216