1 SUBROUTINE zunk1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
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 /
45 bry(1) = 1.0d+3*
d1mach(1)/tol
50 IF (zr.GE.0.0d0)
GO TO 10
60 fn = fnu + dble(float(i-1))
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
68 rast = fn/xzabs(str,sti)
75 s1r = zeta1r(j) - zeta2r(j)
76 s1i = zeta1i(j) - zeta2i(j)
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
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
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)
104 str = s2r*s1r - s2i*s1i
105 s2i = s1r*s2i + s2r*s1i
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
113 yr(i) = s2r*csrr(kflag)
114 yi(i) = s2i*csrr(kflag)
115 IF (kdflg.EQ.2)
GO TO 75
119 IF (rs1.GT.0.0d0)
GO TO 300
123 IF (zr.LT.0.0d0)
GO TO 300
129 IF ((yr(i-1).EQ.zeror).AND.(yi(i-1).EQ.zeroi))
GO TO 70
136 razr = 1.0d0/xzabs(zrr,zri)
144 IF (n.LT.ib)
GO TO 160
149 fn = fnu + dble(float(n-1))
151 IF (mr.NE.0) ipard = 0
153 CALL zunik(zrr, zri, fn, 2, ipard, tol, initd, phidr, phidi,
154 * zet1dr, zet1di, zet2dr, zet2di, sumdr, sumdi, cwrkr(1,3),
156 IF (kode.EQ.1)
GO TO 80
159 rast = fn/xzabs(str,sti)
166 s1r = zet1dr - zet2dr
167 s1i = zet1di - zet2di
170 IF (dabs(rs1).GT.elim)
GO TO 95
171 IF (dabs(rs1).LT.alim)
GO TO 100
175 aphi = xzabs(phidr,phidi)
177 IF (dabs(rs1).LT.elim)
GO TO 100
179 IF (dabs(rs1).GT.0.0d0)
GO TO 300
183 IF (zr.LT.0.0d0)
GO TO 300
203 s2r = ckr*c2r - cki*c2i + s1r
204 s2i = ckr*c2i + cki*c2r + s1i
213 IF (kflag.GE.3)
GO TO 120
217 IF (c2m.LE.ascle)
GO TO 120
224 s1r = s1r*cssr(kflag)
225 s1i = s1i*cssr(kflag)
226 s2r = s2r*cssr(kflag)
227 s2i = s2i*cssr(kflag)
236 fmr = dble(float(mr))
243 fnf = fnu - dble(float(inu))
248 IF (
mod(ifn,2).EQ.0)
GO TO 170
259 fn = fnu + dble(float(kk-1))
265 IF (n.GT.2)
GO TO 175
280 IF ((kk.EQ.n).AND.(ib.LT.n))
GO TO 180
281 IF ((kk.EQ.ib).OR.(kk.EQ.ic))
GO TO 172
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
290 rast = fn/xzabs(str,sti)
297 s1r = -zet1dr + zet2dr
298 s1i = -zet1di + zet2di
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
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
317 str = phidr*sumdr - phidi*sumdi
318 sti = phidr*sumdi + phidi*sumdr
321 str = dexp(s1r)*cssr(iflag)
324 str = s2r*s1r - s2i*s1i
325 s2i = s2r*s1i + s2i*s1r
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
337 s2r = s2r*csrr(iflag)
338 s2i = s2i*csrr(iflag)
344 IF (kode.EQ.1)
GO TO 250
345 CALL zs1s2(zrr, zri, s1r, s1i, s2r, s2i, nw, asc, alim, iuf)
348 yr(kk) = s1r*cspnr - s1i*cspni + s2r
349 yi(kk) = cspnr*s1i + cspni*s1r + s2i
353 IF (c2r.NE.0.0d0 .OR. c2i.NE.0.0d0)
GO TO 255
357 IF (kdflg.EQ.2)
GO TO 275
361 IF (rs1.GT.0.0d0)
GO TO 300
381 fn = dble(float(inu+il))
385 s2r = s1r + (fn+fnf)*(rzr*c2r-rzi*c2i)
386 s2i = s1i + (fn+fnf)*(rzr*c2i+rzi*c2r)
396 IF (kode.EQ.1)
GO TO 280
397 CALL zs1s2(zrr, zri, c1r, c1i, c2r, c2i, nw, asc, alim, iuf)
400 yr(kk) = c1r*cspnr - c1i*cspni + c2r
401 yi(kk) = c1r*cspni + c1i*cspnr + c2i
405 IF (iflag.GE.3)
GO TO 290
409 IF (c2m.LE.ascle)
GO TO 290
416 s1r = s1r*cssr(iflag)
417 s1i = s1i*cssr(iflag)
418 s2r = s2r*cssr(iflag)
419 s2i = s2i*cssr(iflag)
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
double precision function d1mach(i)
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)