1      SUBROUTINE ZUNK1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
2     * ALIM)
3C***BEGIN PROLOGUE  ZUNK1
4C***REFER TO  ZBESK
5C
6C     ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
7C     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
8C     UNIFORM ASYMPTOTIC EXPANSION.
9C     MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
10C     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
11C
12C***ROUTINES CALLED  ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,ZABS
13C***END PROLOGUE  ZUNK1
14C     COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO,
15C    *C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR
16      DOUBLE PRECISION ALIM, ANG, APHI, ASC, ASCLE, BRY, CKI, CKR,
17     * CONER, CRSC, CSCL, CSGNI, CSPNI, CSPNR, CSR, CSRR, CSSR,
18     * CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2M, C2R, ELIM, FMR, FN,
19     * FNF, FNU, PHIDI, PHIDR, PHII, PHIR, PI, RAST, RAZR, RS1, RZI,
20     * RZR, SGN, STI, STR, SUMDI, SUMDR, SUMI, SUMR, S1I, S1R, S2I,
21     * S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R,
22     * ZET1DI, ZET1DR, ZET2DI, ZET2DR, ZI, ZR, ZRI, ZRR, D1MACH, ZABS
23      INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG,
24     * KK, KODE, MR, N, NW, NZ, INITD, IC, IPARD, J
25      DIMENSION BRY(3), INIT(2), YR(N), YI(N), SUMR(2), SUMI(2),
26     * ZETA1R(2), ZETA1I(2), ZETA2R(2), ZETA2I(2), CYR(2), CYI(2),
27     * CWRKR(16,3), CWRKI(16,3), CSSR(3), CSRR(3), PHIR(2), PHII(2)
28      DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
29      DATA PI / 3.14159265358979324D0 /
30C
31      KDFLG = 1
32      NZ = 0
33C-----------------------------------------------------------------------
34C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
35C     THE UNDERFLOW LIMIT
36C-----------------------------------------------------------------------
37      CSCL = 1.0D0/TOL
38      CRSC = TOL
39      CSSR(1) = CSCL
40      CSSR(2) = CONER
41      CSSR(3) = CRSC
42      CSRR(1) = CRSC
43      CSRR(2) = CONER
44      CSRR(3) = CSCL
45      BRY(1) = 1.0D+3*D1MACH(1)/TOL
46      BRY(2) = 1.0D0/BRY(1)
47      BRY(3) = D1MACH(2)
48      ZRR = ZR
49      ZRI = ZI
50      IF (ZR.GE.0.0D0) GO TO 10
51      ZRR = -ZR
52      ZRI = -ZI
53   10 CONTINUE
54      J = 2
55      DO 70 I=1,N
56C-----------------------------------------------------------------------
57C     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
58C-----------------------------------------------------------------------
59        J = 3 - J
60        FN = FNU + DBLE(FLOAT(I-1))
61        INIT(J) = 0
62        CALL ZUNIK(ZRR, ZRI, FN, 2, 0, TOL, INIT(J), PHIR(J), PHII(J),
63     *   ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), SUMR(J), SUMI(J),
64     *   CWRKR(1,J), CWRKI(1,J))
65        IF (KODE.EQ.1) GO TO 20
66        STR = ZRR + ZETA2R(J)
67        STI = ZRI + ZETA2I(J)
68        RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0)))
69        STR = STR*RAST*RAST
70        STI = -STI*RAST*RAST
71        S1R = ZETA1R(J) - STR
72        S1I = ZETA1I(J) - STI
73        GO TO 30
74   20   CONTINUE
75        S1R = ZETA1R(J) - ZETA2R(J)
76        S1I = ZETA1I(J) - ZETA2I(J)
77   30   CONTINUE
78        RS1 = S1R
79C-----------------------------------------------------------------------
80C     TEST FOR UNDERFLOW AND OVERFLOW
81C-----------------------------------------------------------------------
82        IF (DABS(RS1).GT.ELIM) GO TO 60
83        IF (KDFLG.EQ.1) KFLAG = 2
84        IF (DABS(RS1).LT.ALIM) GO TO 40
85C-----------------------------------------------------------------------
86C     REFINE  TEST AND SCALE
87C-----------------------------------------------------------------------
88        APHI = ZABS(CMPLX(PHIR(J),PHII(J),kind=KIND(1.0D0)))
89        RS1 = RS1 + DLOG(APHI)
90        IF (DABS(RS1).GT.ELIM) GO TO 60
91        IF (KDFLG.EQ.1) KFLAG = 1
92        IF (RS1.LT.0.0D0) GO TO 40
93        IF (KDFLG.EQ.1) KFLAG = 3
94   40   CONTINUE
95C-----------------------------------------------------------------------
96C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
97C     EXPONENT EXTREMES
98C-----------------------------------------------------------------------
99        S2R = PHIR(J)*SUMR(J) - PHII(J)*SUMI(J)
100        S2I = PHIR(J)*SUMI(J) + PHII(J)*SUMR(J)
101        STR = DEXP(S1R)*CSSR(KFLAG)
102        S1R = STR*DCOS(S1I)
103        S1I = STR*DSIN(S1I)
104        STR = S2R*S1R - S2I*S1I
105        S2I = S1R*S2I + S2R*S1I
106        S2R = STR
107        IF (KFLAG.NE.1) GO TO 50
108        CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
109        IF (NW.NE.0) GO TO 60
110   50   CONTINUE
111        CYR(KDFLG) = S2R
112        CYI(KDFLG) = S2I
113        YR(I) = S2R*CSRR(KFLAG)
114        YI(I) = S2I*CSRR(KFLAG)
115        IF (KDFLG.EQ.2) GO TO 75
116        KDFLG = 2
117        GO TO 70
118   60   CONTINUE
119        IF (RS1.GT.0.0D0) GO TO 300
120C-----------------------------------------------------------------------
121C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
122C-----------------------------------------------------------------------
123        IF (ZR.LT.0.0D0) GO TO 300
124        KDFLG = 1
125        YR(I)=ZEROR
126        YI(I)=ZEROI
127        NZ=NZ+1
128        IF (I.EQ.1) GO TO 70
129        IF ((YR(I-1).EQ.ZEROR).AND.(YI(I-1).EQ.ZEROI)) GO TO 70
130        YR(I-1)=ZEROR
131        YI(I-1)=ZEROI
132        NZ=NZ+1
133   70 CONTINUE
134      I = N
135   75 CONTINUE
136      RAZR = 1.0D0/ZABS(CMPLX(ZRR,ZRI,kind=KIND(1.0D0)))
137      STR = ZRR*RAZR
138      STI = -ZRI*RAZR
139      RZR = (STR+STR)*RAZR
140      RZI = (STI+STI)*RAZR
141      CKR = FN*RZR
142      CKI = FN*RZI
143      IB = I + 1
144      IF (N.LT.IB) GO TO 160
145C-----------------------------------------------------------------------
146C     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
147C     ON UNDERFLOW.
148C-----------------------------------------------------------------------
149      FN = FNU + DBLE(FLOAT(N-1))
150      IPARD = 1
151      IF (MR.NE.0) IPARD = 0
152      INITD = 0
153      CALL ZUNIK(ZRR, ZRI, FN, 2, IPARD, TOL, INITD, PHIDR, PHIDI,
154     * ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI, CWRKR(1,3),
155     * CWRKI(1,3))
156      IF (KODE.EQ.1) GO TO 80
157      STR = ZRR + ZET2DR
158      STI = ZRI + ZET2DI
159      RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0)))
160      STR = STR*RAST*RAST
161      STI = -STI*RAST*RAST
162      S1R = ZET1DR - STR
163      S1I = ZET1DI - STI
164      GO TO 90
165   80 CONTINUE
166      S1R = ZET1DR - ZET2DR
167      S1I = ZET1DI - ZET2DI
168   90 CONTINUE
169      RS1 = S1R
170      IF (DABS(RS1).GT.ELIM) GO TO 95
171      IF (DABS(RS1).LT.ALIM) GO TO 100
172C----------------------------------------------------------------------------
173C     REFINE ESTIMATE AND TEST
174C-------------------------------------------------------------------------
175      APHI = ZABS(CMPLX(PHIDR,PHIDI,kind=KIND(1.0D0)))
176      RS1 = RS1+DLOG(APHI)
177      IF (DABS(RS1).LT.ELIM) GO TO 100
178   95 CONTINUE
179      IF (DABS(RS1).GT.0.0D0) GO TO 300
180C-----------------------------------------------------------------------
181C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
182C-----------------------------------------------------------------------
183      IF (ZR.LT.0.0D0) GO TO 300
184      NZ = N
185      DO 96 I=1,N
186        YR(I) = ZEROR
187        YI(I) = ZEROI
188   96 CONTINUE
189      RETURN
190C---------------------------------------------------------------------------
191C     FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
192C----------------------------------------------------------------------------
193  100 CONTINUE
194      S1R = CYR(1)
195      S1I = CYI(1)
196      S2R = CYR(2)
197      S2I = CYI(2)
198      C1R = CSRR(KFLAG)
199      ASCLE = BRY(KFLAG)
200      DO 120 I=IB,N
201        C2R = S2R
202        C2I = S2I
203        S2R = CKR*C2R - CKI*C2I + S1R
204        S2I = CKR*C2I + CKI*C2R + S1I
205        S1R = C2R
206        S1I = C2I
207        CKR = CKR + RZR
208        CKI = CKI + RZI
209        C2R = S2R*C1R
210        C2I = S2I*C1R
211        YR(I) = C2R
212        YI(I) = C2I
213        IF (KFLAG.GE.3) GO TO 120
214        STR = DABS(C2R)
215        STI = DABS(C2I)
216        C2M = DMAX1(STR,STI)
217        IF (C2M.LE.ASCLE) GO TO 120
218        KFLAG = KFLAG + 1
219        ASCLE = BRY(KFLAG)
220        S1R = S1R*C1R
221        S1I = S1I*C1R
222        S2R = C2R
223        S2I = C2I
224        S1R = S1R*CSSR(KFLAG)
225        S1I = S1I*CSSR(KFLAG)
226        S2R = S2R*CSSR(KFLAG)
227        S2I = S2I*CSSR(KFLAG)
228        C1R = CSRR(KFLAG)
229  120 CONTINUE
230  160 CONTINUE
231      IF (MR.EQ.0) RETURN
232C-----------------------------------------------------------------------
233C     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
234C-----------------------------------------------------------------------
235      NZ = 0
236      FMR = DBLE(FLOAT(MR))
237      SGN = -DSIGN(PI,FMR)
238C-----------------------------------------------------------------------
239C     CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
240C-----------------------------------------------------------------------
241      CSGNI = SGN
242      INU = INT(SNGL(FNU))
243      FNF = FNU - DBLE(FLOAT(INU))
244      IFN = INU + N - 1
245      ANG = FNF*SGN
246      CSPNR = DCOS(ANG)
247      CSPNI = DSIN(ANG)
248      IF (MOD(IFN,2).EQ.0) GO TO 170
249      CSPNR = -CSPNR
250      CSPNI = -CSPNI
251  170 CONTINUE
252      ASC = BRY(1)
253      IUF = 0
254      KK = N
255      KDFLG = 1
256      IB = IB - 1
257      IC = IB - 1
258      DO 270 K=1,N
259        FN = FNU + DBLE(FLOAT(KK-1))
260C-----------------------------------------------------------------------
261C     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
262C     FUNCTION ABOVE
263C-----------------------------------------------------------------------
264        M=3
265        IF (N.GT.2) GO TO 175
266  172   CONTINUE
267        INITD = INIT(J)
268        PHIDR = PHIR(J)
269        PHIDI = PHII(J)
270        ZET1DR = ZETA1R(J)
271        ZET1DI = ZETA1I(J)
272        ZET2DR = ZETA2R(J)
273        ZET2DI = ZETA2I(J)
274        SUMDR = SUMR(J)
275        SUMDI = SUMI(J)
276        M = J
277        J = 3 - J
278        GO TO 180
279  175   CONTINUE
280        IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 180
281        IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 172
282        INITD = 0
283  180   CONTINUE
284        CALL ZUNIK(ZRR, ZRI, FN, 1, 0, TOL, INITD, PHIDR, PHIDI,
285     *   ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI,
286     *   CWRKR(1,M), CWRKI(1,M))
287        IF (KODE.EQ.1) GO TO 200
288        STR = ZRR + ZET2DR
289        STI = ZRI + ZET2DI
290        RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0)))
291        STR = STR*RAST*RAST
292        STI = -STI*RAST*RAST
293        S1R = -ZET1DR + STR
294        S1I = -ZET1DI + STI
295        GO TO 210
296  200   CONTINUE
297        S1R = -ZET1DR + ZET2DR
298        S1I = -ZET1DI + ZET2DI
299  210   CONTINUE
300C-----------------------------------------------------------------------
301C     TEST FOR UNDERFLOW AND OVERFLOW
302C-----------------------------------------------------------------------
303        RS1 = S1R
304        IF (DABS(RS1).GT.ELIM) GO TO 260
305        IF (KDFLG.EQ.1) IFLAG = 2
306        IF (DABS(RS1).LT.ALIM) GO TO 220
307C-----------------------------------------------------------------------
308C     REFINE  TEST AND SCALE
309C-----------------------------------------------------------------------
310        APHI = ZABS(CMPLX(PHIDR,PHIDI,kind=KIND(1.0D0)))
311        RS1 = RS1 + DLOG(APHI)
312        IF (DABS(RS1).GT.ELIM) GO TO 260
313        IF (KDFLG.EQ.1) IFLAG = 1
314        IF (RS1.LT.0.0D0) GO TO 220
315        IF (KDFLG.EQ.1) IFLAG = 3
316  220   CONTINUE
317        STR = PHIDR*SUMDR - PHIDI*SUMDI
318        STI = PHIDR*SUMDI + PHIDI*SUMDR
319        S2R = -CSGNI*STI
320        S2I = CSGNI*STR
321        STR = DEXP(S1R)*CSSR(IFLAG)
322        S1R = STR*DCOS(S1I)
323        S1I = STR*DSIN(S1I)
324        STR = S2R*S1R - S2I*S1I
325        S2I = S2R*S1I + S2I*S1R
326        S2R = STR
327        IF (IFLAG.NE.1) GO TO 230
328        CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
329        IF (NW.EQ.0) GO TO 230
330        S2R = ZEROR
331        S2I = ZEROI
332  230   CONTINUE
333        CYR(KDFLG) = S2R
334        CYI(KDFLG) = S2I
335        C2R = S2R
336        C2I = S2I
337        S2R = S2R*CSRR(IFLAG)
338        S2I = S2I*CSRR(IFLAG)
339C-----------------------------------------------------------------------
340C     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
341C-----------------------------------------------------------------------
342        S1R = YR(KK)
343        S1I = YI(KK)
344        IF (KODE.EQ.1) GO TO 250
345        CALL ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NW, ASC, ALIM, IUF)
346        NZ = NZ + NW
347  250   CONTINUE
348        YR(KK) = S1R*CSPNR - S1I*CSPNI + S2R
349        YI(KK) = CSPNR*S1I + CSPNI*S1R + S2I
350        KK = KK - 1
351        CSPNR = -CSPNR
352        CSPNI = -CSPNI
353        IF (C2R.NE.0.0D0 .OR. C2I.NE.0.0D0) GO TO 255
354        KDFLG = 1
355        GO TO 270
356  255   CONTINUE
357        IF (KDFLG.EQ.2) GO TO 275
358        KDFLG = 2
359        GO TO 270
360  260   CONTINUE
361        IF (RS1.GT.0.0D0) GO TO 300
362        S2R = ZEROR
363        S2I = ZEROI
364        GO TO 230
365  270 CONTINUE
366      K = N
367  275 CONTINUE
368      IL = N - K
369      IF (IL.EQ.0) RETURN
370C-----------------------------------------------------------------------
371C     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
372C     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
373C     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
374C-----------------------------------------------------------------------
375      S1R = CYR(1)
376      S1I = CYI(1)
377      S2R = CYR(2)
378      S2I = CYI(2)
379      CSR = CSRR(IFLAG)
380      ASCLE = BRY(IFLAG)
381      FN = DBLE(FLOAT(INU+IL))
382      DO 290 I=1,IL
383        C2R = S2R
384        C2I = S2I
385        S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I)
386        S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R)
387        S1R = C2R
388        S1I = C2I
389        FN = FN - 1.0D0
390        C2R = S2R*CSR
391        C2I = S2I*CSR
392        CKR = C2R
393        CKI = C2I
394        C1R = YR(KK)
395        C1I = YI(KK)
396        IF (KODE.EQ.1) GO TO 280
397        CALL ZS1S2(ZRR, ZRI, C1R, C1I, C2R, C2I, NW, ASC, ALIM, IUF)
398        NZ = NZ + NW
399  280   CONTINUE
400        YR(KK) = C1R*CSPNR - C1I*CSPNI + C2R
401        YI(KK) = C1R*CSPNI + C1I*CSPNR + C2I
402        KK = KK - 1
403        CSPNR = -CSPNR
404        CSPNI = -CSPNI
405        IF (IFLAG.GE.3) GO TO 290
406        C2R = DABS(CKR)
407        C2I = DABS(CKI)
408        C2M = DMAX1(C2R,C2I)
409        IF (C2M.LE.ASCLE) GO TO 290
410        IFLAG = IFLAG + 1
411        ASCLE = BRY(IFLAG)
412        S1R = S1R*CSR
413        S1I = S1I*CSR
414        S2R = CKR
415        S2I = CKI
416        S1R = S1R*CSSR(IFLAG)
417        S1I = S1I*CSSR(IFLAG)
418        S2R = S2R*CSSR(IFLAG)
419        S2I = S2I*CSSR(IFLAG)
420        CSR = CSRR(IFLAG)
421  290 CONTINUE
422      RETURN
423  300 CONTINUE
424      NZ = -1
425      RETURN
426      END
427