1      SUBROUTINE CUNK1(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
2C***BEGIN PROLOGUE  CUNK1
3C***REFER TO  CBESK
4C
5C     CUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
6C     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
7C     UNIFORM ASYMPTOTIC EXPANSION.
8C     MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
9C     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
10C
11C***ROUTINES CALLED  CS1S2,CUCHK,CUNIK,R1MACH
12C***END PROLOGUE  CUNK1
13      COMPLEX CFN, CK, CONE, CRSC, CS, CSCL, CSGN, CSPN, CSR, CSS,
14     * CWRK, CY, CZERO, C1, C2, PHI,  RZ, SUM,  S1, S2, Y, Z,
15     * ZETA1,  ZETA2,  ZR, PHID, ZETA1D, ZETA2D, SUMD
16      REAL ALIM, ANG, APHI, ASC, ASCLE, BRY, CPN, C2I, C2M, C2R, ELIM,
17     * FMR, FN, FNF, FNU, PI, RS1, SGN, SPN, TOL, X, R1MACH
18      INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG,
19     * KK, KODE, MR, N, NW, NZ, J, IPARD, INITD, IC
20      DIMENSION BRY(3), INIT(2), Y(N), SUM(2), PHI(2), ZETA1(2),
21     * ZETA2(2), CY(2), CWRK(16,3), CSS(3), CSR(3)
22      DATA CZERO, CONE / (0.0E0,0.0E0) , (1.0E0,0.0E0) /
23      DATA PI / 3.14159265358979324E0 /
24C
25      KDFLG = 1
26      NZ = 0
27C-----------------------------------------------------------------------
28C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
29C     THE UNDERFLOW LIMIT
30C-----------------------------------------------------------------------
31      CSCL = CMPLX(1.0E0/TOL,0.0E0)
32      CRSC = CMPLX(TOL,0.0E0)
33      CSS(1) = CSCL
34      CSS(2) = CONE
35      CSS(3) = CRSC
36      CSR(1) = CRSC
37      CSR(2) = CONE
38      CSR(3) = CSCL
39      BRY(1) = 1.0E+3*R1MACH(1)/TOL
40      BRY(2) = 1.0E0/BRY(1)
41      BRY(3) = R1MACH(2)
42      X = REAL(Z)
43      ZR = Z
44      IF (X.LT.0.0E0) ZR = -Z
45      J=2
46      DO 70 I=1,N
47C-----------------------------------------------------------------------
48C     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
49C-----------------------------------------------------------------------
50        J = 3 - J
51        FN = FNU + FLOAT(I-1)
52        INIT(J) = 0
53        CALL CUNIK(ZR, FN, 2, 0, TOL, INIT(J), PHI(J), ZETA1(J),
54     *   ZETA2(J), SUM(J), CWRK(1,J))
55        IF (KODE.EQ.1) GO TO 20
56        CFN = CMPLX(FN,0.0E0)
57        S1 = ZETA1(J) - CFN*(CFN/(ZR+ZETA2(J)))
58        GO TO 30
59   20   CONTINUE
60        S1 = ZETA1(J) - ZETA2(J)
61   30   CONTINUE
62C-----------------------------------------------------------------------
63C     TEST FOR UNDERFLOW AND OVERFLOW
64C-----------------------------------------------------------------------
65        RS1 = REAL(S1)
66        IF (ABS(RS1).GT.ELIM) GO TO 60
67        IF (KDFLG.EQ.1) KFLAG = 2
68        IF (ABS(RS1).LT.ALIM) GO TO 40
69C-----------------------------------------------------------------------
70C     REFINE  TEST AND SCALE
71C-----------------------------------------------------------------------
72        APHI = CABS(PHI(J))
73        RS1 = RS1 + ALOG(APHI)
74        IF (ABS(RS1).GT.ELIM) GO TO 60
75        IF (KDFLG.EQ.1) KFLAG = 1
76        IF (RS1.LT.0.0E0) GO TO 40
77        IF (KDFLG.EQ.1) KFLAG = 3
78   40   CONTINUE
79C-----------------------------------------------------------------------
80C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
81C     EXPONENT EXTREMES
82C-----------------------------------------------------------------------
83        S2 = PHI(J)*SUM(J)
84        C2R = REAL(S1)
85        C2I = AIMAG(S1)
86        C2M = EXP(C2R)*REAL(CSS(KFLAG))
87        S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I))
88        S2 = S2*S1
89        IF (KFLAG.NE.1) GO TO 50
90        CALL CUCHK(S2, NW, BRY(1), TOL)
91        IF (NW.NE.0) GO TO 60
92   50   CONTINUE
93        CY(KDFLG) = S2
94        Y(I) = S2*CSR(KFLAG)
95        IF (KDFLG.EQ.2) GO TO 75
96        KDFLG = 2
97        GO TO 70
98   60   CONTINUE
99        IF (RS1.GT.0.0E0) GO TO 290
100C-----------------------------------------------------------------------
101C     FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
102C-----------------------------------------------------------------------
103        IF (X.LT.0.0E0) GO TO 290
104        KDFLG = 1
105        Y(I) = CZERO
106        NZ=NZ+1
107        IF (I.EQ.1) GO TO 70
108        IF (Y(I-1).EQ.CZERO) GO TO 70
109        Y(I-1) = CZERO
110        NZ=NZ+1
111   70 CONTINUE
112      I=N
113   75 CONTINUE
114      RZ = CMPLX(2.0E0,0.0E0)/ZR
115      CK = CMPLX(FN,0.0E0)*RZ
116      IB = I+1
117      IF (N.LT.IB) GO TO 160
118C-----------------------------------------------------------------------
119C     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO
120C     ON UNDERFLOW
121C-----------------------------------------------------------------------
122      FN = FNU+FLOAT(N-1)
123      IPARD = 1
124      IF (MR.NE.0) IPARD = 0
125      INITD = 0
126      CALL CUNIK(ZR,FN,2,IPARD,TOL,INITD,PHID,ZETA1D,ZETA2D,SUMD,
127     *CWRK(1,3))
128      IF (KODE.EQ.1) GO TO 80
129      CFN=CMPLX(FN,0.0E0)
130      S1=ZETA1D-CFN*(CFN/(ZR+ZETA2D))
131      GO TO 90
132   80 CONTINUE
133      S1=ZETA1D-ZETA2D
134   90 CONTINUE
135      RS1=REAL(S1)
136      IF (ABS(RS1).GT.ELIM) GO TO 95
137      IF (ABS(RS1).LT.ALIM) GO TO 100
138C-----------------------------------------------------------------------
139C     REFINE ESTIMATE AND TEST
140C-----------------------------------------------------------------------
141      APHI=CABS(PHID)
142      RS1=RS1+ALOG(APHI)
143      IF (ABS(RS1).LT.ELIM) GO TO 100
144   95 CONTINUE
145      IF (RS1.GT.0.0E0) GO TO 290
146C-----------------------------------------------------------------------
147C     FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
148C-----------------------------------------------------------------------
149      IF (X.LT.0.0E0) GO TO 290
150      NZ=N
151      DO 96 I=1,N
152        Y(I) = CZERO
153   96 CONTINUE
154      RETURN
155  100 CONTINUE
156C-----------------------------------------------------------------------
157C     RECUR FORWARD FOR REMAINDER OF THE SEQUENCE
158C-----------------------------------------------------------------------
159      S1 = CY(1)
160      S2 = CY(2)
161      C1 = CSR(KFLAG)
162      ASCLE = BRY(KFLAG)
163      DO 120 I=IB,N
164        C2 = S2
165        S2 = CK*S2 + S1
166        S1 = C2
167        CK = CK + RZ
168        C2 = S2*C1
169        Y(I) = C2
170        IF (KFLAG.GE.3) GO TO 120
171        C2R = REAL(C2)
172        C2I = AIMAG(C2)
173        C2R = ABS(C2R)
174        C2I = ABS(C2I)
175        C2M = AMAX1(C2R,C2I)
176        IF (C2M.LE.ASCLE) GO TO 120
177        KFLAG = KFLAG + 1
178        ASCLE = BRY(KFLAG)
179        S1 = S1*C1
180        S2 = C2
181        S1 = S1*CSS(KFLAG)
182        S2 = S2*CSS(KFLAG)
183        C1 = CSR(KFLAG)
184  120 CONTINUE
185  160 CONTINUE
186      IF (MR.EQ.0) RETURN
187C-----------------------------------------------------------------------
188C     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0E0
189C-----------------------------------------------------------------------
190      NZ = 0
191      FMR = FLOAT(MR)
192      SGN = -SIGN(PI,FMR)
193C-----------------------------------------------------------------------
194C     CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
195C-----------------------------------------------------------------------
196      CSGN = CMPLX(0.0E0,SGN)
197      INU = INT(FNU)
198      FNF = FNU - FLOAT(INU)
199      IFN = INU + N - 1
200      ANG = FNF*SGN
201      CPN = COS(ANG)
202      SPN = SIN(ANG)
203      CSPN = CMPLX(CPN,SPN)
204      IF (MOD(IFN,2).EQ.1) CSPN = -CSPN
205      ASC = BRY(1)
206      KK = N
207      IUF = 0
208      KDFLG = 1
209      IB = IB-1
210      IC = IB-1
211      DO 260 K=1,N
212        FN = FNU + FLOAT(KK-1)
213C-----------------------------------------------------------------------
214C     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
215C     FUNCTION ABOVE
216C-----------------------------------------------------------------------
217        M=3
218        IF (N.GT.2) GO TO 175
219  170   CONTINUE
220        INITD = INIT(J)
221        PHID = PHI(J)
222        ZETA1D = ZETA1(J)
223        ZETA2D = ZETA2(J)
224        SUMD = SUM(J)
225        M = J
226        J = 3 - J
227        GO TO 180
228  175   CONTINUE
229        IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 180
230        IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 170
231        INITD = 0
232  180   CONTINUE
233        CALL CUNIK(ZR, FN, 1, 0, TOL, INITD, PHID, ZETA1D,
234     *   ZETA2D, SUMD, CWRK(1,M))
235        IF (KODE.EQ.1) GO TO 190
236        CFN = CMPLX(FN,0.0E0)
237        S1 = -ZETA1D + CFN*(CFN/(ZR+ZETA2D))
238        GO TO 200
239  190   CONTINUE
240        S1 = -ZETA1D + ZETA2D
241  200   CONTINUE
242C-----------------------------------------------------------------------
243C     TEST FOR UNDERFLOW AND OVERFLOW
244C-----------------------------------------------------------------------
245        RS1 = REAL(S1)
246        IF (ABS(RS1).GT.ELIM) GO TO 250
247        IF (KDFLG.EQ.1) IFLAG = 2
248        IF (ABS(RS1).LT.ALIM) GO TO 210
249C-----------------------------------------------------------------------
250C     REFINE  TEST AND SCALE
251C-----------------------------------------------------------------------
252        APHI = CABS(PHID)
253        RS1 = RS1 + ALOG(APHI)
254        IF (ABS(RS1).GT.ELIM) GO TO 250
255        IF (KDFLG.EQ.1) IFLAG = 1
256        IF (RS1.LT.0.0E0) GO TO 210
257        IF (KDFLG.EQ.1) IFLAG = 3
258  210   CONTINUE
259        S2 = CSGN*PHID*SUMD
260        C2R = REAL(S1)
261        C2I = AIMAG(S1)
262        C2M = EXP(C2R)*REAL(CSS(IFLAG))
263        S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I))
264        S2 = S2*S1
265        IF (IFLAG.NE.1) GO TO 220
266        CALL CUCHK(S2, NW, BRY(1), TOL)
267        IF (NW.NE.0) S2 = CMPLX(0.0E0,0.0E0)
268  220   CONTINUE
269        CY(KDFLG) = S2
270        C2 = S2
271        S2 = S2*CSR(IFLAG)
272C-----------------------------------------------------------------------
273C     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
274C-----------------------------------------------------------------------
275        S1 = Y(KK)
276        IF (KODE.EQ.1) GO TO 240
277        CALL CS1S2(ZR, S1, S2, NW, ASC, ALIM, IUF)
278        NZ = NZ + NW
279  240   CONTINUE
280        Y(KK) = S1*CSPN + S2
281        KK = KK - 1
282        CSPN = -CSPN
283        IF (C2.NE.CZERO) GO TO 245
284        KDFLG = 1
285        GO TO 260
286  245   CONTINUE
287        IF (KDFLG.EQ.2) GO TO 265
288        KDFLG = 2
289        GO TO 260
290  250   CONTINUE
291        IF (RS1.GT.0.0E0) GO TO 290
292        S2 = CZERO
293        GO TO 220
294  260 CONTINUE
295      K = N
296  265 CONTINUE
297      IL = N - K
298      IF (IL.EQ.0) RETURN
299C-----------------------------------------------------------------------
300C     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
301C     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
302C     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
303C-----------------------------------------------------------------------
304      S1 = CY(1)
305      S2 = CY(2)
306      CS = CSR(IFLAG)
307      ASCLE = BRY(IFLAG)
308      FN = FLOAT(INU+IL)
309      DO 280 I=1,IL
310        C2 = S2
311        S2 = S1 + CMPLX(FN+FNF,0.0E0)*RZ*S2
312        S1 = C2
313        FN = FN - 1.0E0
314        C2 = S2*CS
315        CK = C2
316        C1 = Y(KK)
317        IF (KODE.EQ.1) GO TO 270
318        CALL CS1S2(ZR, C1, C2, NW, ASC, ALIM, IUF)
319        NZ = NZ + NW
320  270   CONTINUE
321        Y(KK) = C1*CSPN + C2
322        KK = KK - 1
323        CSPN = -CSPN
324        IF (IFLAG.GE.3) GO TO 280
325        C2R = REAL(CK)
326        C2I = AIMAG(CK)
327        C2R = ABS(C2R)
328        C2I = ABS(C2I)
329        C2M = AMAX1(C2R,C2I)
330        IF (C2M.LE.ASCLE) GO TO 280
331        IFLAG = IFLAG + 1
332        ASCLE = BRY(IFLAG)
333        S1 = S1*CS
334        S2 = CK
335        S1 = S1*CSS(IFLAG)
336        S2 = S2*CSS(IFLAG)
337        CS = CSR(IFLAG)
338  280 CONTINUE
339      RETURN
340  290 CONTINUE
341      NZ = -1
342      RETURN
343      END
344