1 SUBROUTINE zunk2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
20 DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGDI,
21 * argdr, argi, argr, asc, ascle, asumdi, asumdr, asumi, asumr,
22 * bry, bsumdi, bsumdr, bsumi, bsumr, car, cipi, cipr, cki, ckr,
23 * coner, crsc, cr1i, cr1r, cr2i, cr2r, cscl, csgni, csi,
24 * cspni, cspnr, csr, csrr, cssr, cyi, cyr, c1i, c1r, c2i, c2m,
25 * c2r, daii, dair, elim, fmr, fn, fnf, fnu, hpi, phidi, phidr,
26 * phii, phir, pi, pti, ptr, rast, razr, rs1, rzi, rzr, sar, sgn,
27 * sti, str, s1i, s1r, s2i, s2r, tol, yi, yr, yy, zbi, zbr, zeroi,
28 * zeror, zeta1i, zeta1r, zeta2i, zeta2r, zet1di, zet1dr, zet2di,
30 INTEGER I, IB, IFLAG, IFN, IL, IN, INU, IUF, K, KDFLG, KFLAG, KK,
31 * kode, mr, n, nai, ndai, nw, nz, idum, j, ipard, ic
32 dimension bry(3), yr(n), yi(n), asumr(2), asumi(2), bsumr(2),
33 * bsumi(2), phir(2), phii(2), argr(2), argi(2), zeta1r(2),
34 * zeta1i(2), zeta2r(2), zeta2i(2), cyr(2), cyi(2), cipr(4),
35 * cipi(4), cssr(3), csrr(3)
36 DATA zeror,zeroi,coner,cr1r,cr1i,cr2r,cr2i /
37 1 0.0d0, 0.0d0, 1.0d0,
38 1 1.0d0,1.73205080756887729d0 , -0.5d0,-8.66025403784438647d-01 /
40 1 1.57079632679489662d+00, 3.14159265358979324d+00,
41 1 1.26551212348464539d+00/
42 DATA cipr(1),cipi(1),cipr(2),cipi(2),cipr(3),cipi(3),cipr(4),
44 1 1.0d0,0.0d0 , 0.0d0,-1.0d0 , -1.0d0,0.0d0 , 0.0d0,1.0d0 /
60 bry(1) = 1.0d+3*
d1mach(1)/tol
65 IF (zr.GE.0.0d0)
GO TO 10
75 fnf = fnu - dble(float(inu))
82 str = c2r*cipr(kk) - c2i*cipi(kk)
83 sti = c2r*cipi(kk) + c2i*cipr(kk)
84 csr = cr1r*str - cr1i*sti
85 csi = cr1r*sti + cr1i*str
86 IF (yy.GT.0.0d0)
GO TO 20
101 fn = fnu + dble(float(i-1))
102 CALL zunhj(znr, zni, fn, 0, tol, phir(j), phii(j), argr(j),
103 * argi(j), zeta1r(j), zeta1i(j), zeta2r(j), zeta2i(j), asumr(j),
104 * asumi(j), bsumr(j), bsumi(j))
105 IF (kode.EQ.1)
GO TO 30
106 str = zbr + zeta2r(j)
107 sti = zbi + zeta2i(j)
108 rast = fn/
xzabs(str,sti)
111 s1r = zeta1r(j) - str
112 s1i = zeta1i(j) - sti
115 s1r = zeta1r(j) - zeta2r(j)
116 s1i = zeta1i(j) - zeta2i(j)
122 IF (dabs(rs1).GT.elim)
GO TO 70
123 IF (kdflg.EQ.1) kflag = 2
124 IF (dabs(rs1).LT.alim)
GO TO 50
128 aphi =
xzabs(phir(j),phii(j))
129 aarg =
xzabs(argr(j),argi(j))
130 rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
131 IF (dabs(rs1).GT.elim)
GO TO 70
132 IF (kdflg.EQ.1) kflag = 1
133 IF (rs1.LT.0.0d0)
GO TO 50
134 IF (kdflg.EQ.1) kflag = 3
140 c2r = argr(j)*cr2r - argi(j)*cr2i
141 c2i = argr(j)*cr2i + argi(j)*cr2r
142 CALL zairy(c2r, c2i, 0, 2, air, aii, nai, idum)
143 CALL zairy(c2r, c2i, 1, 2, dair, daii, ndai, idum)
144 str = dair*bsumr(j) - daii*bsumi(j)
145 sti = dair*bsumi(j) + daii*bsumr(j)
146 ptr = str*cr2r - sti*cr2i
147 pti = str*cr2i + sti*cr2r
148 str = ptr + (air*asumr(j)-aii*asumi(j))
149 sti = pti + (air*asumi(j)+aii*asumr(j))
150 ptr = str*phir(j) - sti*phii(j)
151 pti = str*phii(j) + sti*phir(j)
152 s2r = ptr*csr - pti*csi
153 s2i = ptr*csi + pti*csr
154 str = dexp(s1r)*cssr(kflag)
157 str = s2r*s1r - s2i*s1i
158 s2i = s1r*s2i + s2r*s1i
160 IF (kflag.NE.1)
GO TO 60
161 CALL zuchk(s2r, s2i, nw, bry(1), tol)
162 IF (nw.NE.0)
GO TO 70
164 IF (yy.LE.0.0d0) s2i = -s2i
167 yr(i) = s2r*csrr(kflag)
168 yi(i) = s2i*csrr(kflag)
172 IF (kdflg.EQ.2)
GO TO 85
176 IF (rs1.GT.0.0d0)
GO TO 320
180 IF (zr.LT.0.0d0)
GO TO 320
189 IF ((yr(i-1).EQ.zeror).AND.(yi(i-1).EQ.zeroi))
GO TO 80
196 razr = 1.0d0/
xzabs(zrr,zri)
204 IF (n.LT.ib)
GO TO 180
209 fn = fnu + dble(float(n-1))
211 IF (mr.NE.0) ipard = 0
212 CALL zunhj(znr, zni, fn, ipard, tol, phidr, phidi, argdr, argdi,
213 * zet1dr, zet1di, zet2dr, zet2di, asumdr, asumdi, bsumdr, bsumdi)
214 IF (kode.EQ.1)
GO TO 90
217 rast = fn/
xzabs(str,sti)
224 s1r = zet1dr - zet2dr
225 s1i = zet1di - zet2di
228 IF (dabs(rs1).GT.elim)
GO TO 105
229 IF (dabs(rs1).LT.alim)
GO TO 120
233 aphi =
xzabs(phidr,phidi)
235 IF (dabs(rs1).LT.elim)
GO TO 120
237 IF (rs1.GT.0.0d0)
GO TO 320
241 IF (zr.LT.0.0d0)
GO TO 320
258 s2r = ckr*c2r - cki*c2i + s1r
259 s2i = ckr*c2i + cki*c2r + s1i
268 IF (kflag.GE.3)
GO TO 130
272 IF (c2m.LE.ascle)
GO TO 130
279 s1r = s1r*cssr(kflag)
280 s1i = s1i*cssr(kflag)
281 s2r = s2r*cssr(kflag)
282 s2i = s2i*cssr(kflag)
291 fmr = dble(float(mr))
297 IF (yy.LE.0.0d0) csgni = -csgni
302 IF (
mod(ifn,2).EQ.0)
GO TO 190
317 str = csr*c2r + csi*c2i
318 csi = -csr*c2i + csi*c2r
327 fn = fnu + dble(float(kk-1))
332 IF (n.GT.2)
GO TO 175
349 IF ((kk.EQ.n).AND.(ib.LT.n))
GO TO 210
350 IF ((kk.EQ.ib).OR.(kk.EQ.ic))
GO TO 172
351 CALL zunhj(znr, zni, fn, 0, tol, phidr, phidi, argdr,
352 * argdi, zet1dr, zet1di, zet2dr, zet2di, asumdr,
353 * asumdi, bsumdr, bsumdi)
355 IF (kode.EQ.1)
GO TO 220
358 rast = fn/
xzabs(str,sti)
365 s1r = -zet1dr + zet2dr
366 s1i = -zet1di + zet2di
372 IF (dabs(rs1).GT.elim)
GO TO 280
373 IF (kdflg.EQ.1) iflag = 2
374 IF (dabs(rs1).LT.alim)
GO TO 240
378 aphi =
xzabs(phidr,phidi)
379 aarg =
xzabs(argdr,argdi)
380 rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
381 IF (dabs(rs1).GT.elim)
GO TO 280
382 IF (kdflg.EQ.1) iflag = 1
383 IF (rs1.LT.0.0d0)
GO TO 240
384 IF (kdflg.EQ.1) iflag = 3
386 CALL zairy(argdr, argdi, 0, 2, air, aii, nai, idum)
387 CALL zairy(argdr, argdi, 1, 2, dair, daii, ndai, idum)
388 str = dair*bsumdr - daii*bsumdi
389 sti = dair*bsumdi + daii*bsumdr
390 str = str + (air*asumdr-aii*asumdi)
391 sti = sti + (air*asumdi+aii*asumdr)
392 ptr = str*phidr - sti*phidi
393 pti = str*phidi + sti*phidr
394 s2r = ptr*csr - pti*csi
395 s2i = ptr*csi + pti*csr
396 str = dexp(s1r)*cssr(iflag)
399 str = s2r*s1r - s2i*s1i
400 s2i = s2r*s1i + s2i*s1r
402 IF (iflag.NE.1)
GO TO 250
403 CALL zuchk(s2r, s2i, nw, bry(1), tol)
404 IF (nw.EQ.0)
GO TO 250
408 IF (yy.LE.0.0d0) s2i = -s2i
413 s2r = s2r*csrr(iflag)
414 s2i = s2i*csrr(iflag)
420 IF (kode.EQ.1)
GO TO 270
421 CALL zs1s2(zrr, zri, s1r, s1i, s2r, s2i, nw, asc, alim, iuf)
424 yr(kk) = s1r*cspnr - s1i*cspni + s2r
425 yi(kk) = s1r*cspni + s1i*cspnr + s2i
432 IF (c2r.NE.0.0d0 .OR. c2i.NE.0.0d0)
GO TO 255
436 IF (kdflg.EQ.2)
GO TO 295
440 IF (rs1.GT.0.0d0)
GO TO 320
460 fn = dble(float(inu+il))
464 s2r = s1r + (fn+fnf)*(rzr*c2r-rzi*c2i)
465 s2i = s1i + (fn+fnf)*(rzr*c2i+rzi*c2r)
475 IF (kode.EQ.1)
GO TO 300
476 CALL zs1s2(zrr, zri, c1r, c1i, c2r, c2i, nw, asc, alim, iuf)
479 yr(kk) = c1r*cspnr - c1i*cspni + c2r
480 yi(kk) = c1r*cspni + c1i*cspnr + c2i
484 IF (iflag.GE.3)
GO TO 310
488 IF (c2m.LE.ascle)
GO TO 310
495 s1r = s1r*cssr(iflag)
496 s1i = s1i*cssr(iflag)
497 s2r = s2r*cssr(iflag)
498 s2i = s2i*cssr(iflag)
double precision function d1mach(i)
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
double precision function xzabs(ZR, ZI)
subroutine zairy(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR)
subroutine zs1s2(ZRR, ZRI, S1R, S1I, S2R, S2I, NZ, ASCLE, ALIM, IUF)
subroutine zuchk(YR, YI, NZ, ASCLE, TOL)
subroutine zunhj(ZR, ZI, FNU, IPMTR, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
subroutine zunk2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, ALIM)