GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zunk1.f
Go to the documentation of this file.
1  SUBROUTINE zunk1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
2  * ALIM)
3 C***BEGIN PROLOGUE ZUNK1
4 C***REFER TO ZBESK
5 C
6 C ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
7 C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
8 C UNIFORM ASYMPTOTIC EXPANSION.
9 C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
10 C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
11 C
12 C***ROUTINES CALLED ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,XZABS
13 C***END PROLOGUE ZUNK1
14 C COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO,
15 C *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, xzabs
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 /
30 C
31  kdflg = 1
32  nz = 0
33 C-----------------------------------------------------------------------
34 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
35 C THE UNDERFLOW LIMIT
36 C-----------------------------------------------------------------------
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
56 C-----------------------------------------------------------------------
57 C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
58 C-----------------------------------------------------------------------
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/xzabs(str,sti)
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
79 C-----------------------------------------------------------------------
80 C TEST FOR UNDERFLOW AND OVERFLOW
81 C-----------------------------------------------------------------------
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
85 C-----------------------------------------------------------------------
86 C REFINE TEST AND SCALE
87 C-----------------------------------------------------------------------
88  aphi = xzabs(phir(j),phii(j))
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
95 C-----------------------------------------------------------------------
96 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
97 C EXPONENT EXTREMES
98 C-----------------------------------------------------------------------
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
120 C-----------------------------------------------------------------------
121 C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
122 C-----------------------------------------------------------------------
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/xzabs(zrr,zri)
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
145 C-----------------------------------------------------------------------
146 C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
147 C ON UNDERFLOW.
148 C-----------------------------------------------------------------------
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/xzabs(str,sti)
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
172 C----------------------------------------------------------------------------
173 C REFINE ESTIMATE AND TEST
174 C-------------------------------------------------------------------------
175  aphi = xzabs(phidr,phidi)
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
180 C-----------------------------------------------------------------------
181 C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
182 C-----------------------------------------------------------------------
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
190 C---------------------------------------------------------------------------
191 C FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
192 C----------------------------------------------------------------------------
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
232 C-----------------------------------------------------------------------
233 C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
234 C-----------------------------------------------------------------------
235  nz = 0
236  fmr = dble(float(mr))
237  sgn = -dsign(pi,fmr)
238 C-----------------------------------------------------------------------
239 C CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
240 C-----------------------------------------------------------------------
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))
260 C-----------------------------------------------------------------------
261 C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
262 C FUNCTION ABOVE
263 C-----------------------------------------------------------------------
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/xzabs(str,sti)
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
300 C-----------------------------------------------------------------------
301 C TEST FOR UNDERFLOW AND OVERFLOW
302 C-----------------------------------------------------------------------
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
307 C-----------------------------------------------------------------------
308 C REFINE TEST AND SCALE
309 C-----------------------------------------------------------------------
310  aphi = xzabs(phidr,phidi)
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)
339 C-----------------------------------------------------------------------
340 C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
341 C-----------------------------------------------------------------------
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
370 C-----------------------------------------------------------------------
371 C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
372 C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
373 C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
374 C-----------------------------------------------------------------------
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
mod
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:860
d1mach
double precision function d1mach(i)
Definition: d1mach.f:20
dimension
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)