1 SUBROUTINE cunk1(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
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 /
31 cscl =
cmplx(1.0e0/tol,0.0e0)
32 crsc =
cmplx(tol,0.0e0)
39 bry(1) = 1.0
e+3*r1mach(1)/tol
44 IF (x.LT.0.0e0) zr = -z
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
57 s1 = zeta1(j) - cfn*(cfn/(zr+zeta2(j)))
60 s1 = zeta1(j) - zeta2(j)
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
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
86 c2m =
exp(c2r)*
REAL(css(kflag))
89 IF (kflag.NE.1) go
to 50
90 CALL
cuchk(s2, nw, bry(1), tol)
95 IF (kdflg.EQ.2) go
to 75
99 IF (rs1.GT.0.0e0) go
to 290
103 IF (x.LT.0.0e0) go
to 290
108 IF (y(i-1).EQ.czero) go
to 70
114 rz =
cmplx(2.0e0,0.0e0)/zr
115 ck =
cmplx(fn,0.0e0)*rz
117 IF (n.LT.ib) go
to 160
124 IF (mr.NE.0) ipard = 0
126 CALL
cunik(zr,fn,2,ipard,tol,initd,phid,zeta1d,zeta2d,sumd,
128 IF (kode.EQ.1) go
to 80
130 s1=zeta1d-cfn*(cfn/(zr+zeta2d))
136 IF (
abs(rs1).GT.elim) go
to 95
137 IF (
abs(rs1).LT.alim) go
to 100
143 IF (
abs(rs1).LT.elim) go
to 100
145 IF (rs1.GT.0.0e0) go
to 290
149 IF (x.LT.0.0e0) go
to 290
170 IF (kflag.GE.3) go
to 120
176 IF (c2m.LE.ascle) go
to 120
196 csgn =
cmplx(0.0e0,sgn)
198 fnf = fnu - float(inu)
203 cspn =
cmplx(cpn,spn)
204 IF (
mod(ifn,2).EQ.1) cspn = -cspn
212 fn = fnu + float(kk-1)
218 IF (n.GT.2) go
to 175
229 IF ((kk.EQ.n).AND.(ib.LT.n)) go
to 180
230 IF ((kk.EQ.ib).OR.(kk.EQ.ic)) go
to 170
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))
240 s1 = -zeta1d + zeta2d
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
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
262 c2m =
exp(c2r)*
REAL(css(iflag))
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)
276 IF (kode.EQ.1) go
to 240
277 CALL
cs1s2(zr, s1, s2, nw, asc, alim, iuf)
283 IF (c2.NE.czero) go
to 245
287 IF (kdflg.EQ.2) go
to 265
291 IF (rs1.GT.0.0e0) go
to 290
311 s2 = s1 +
cmplx(fn+fnf,0.0e0)*rz*s2
317 IF (kode.EQ.1) go
to 270
318 CALL
cs1s2(zr, c1, c2, nw, asc, alim, iuf)
324 IF (iflag.GE.3) go
to 280
330 IF (c2m.LE.ascle) go
to 280
subroutine cunk1(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
subroutine cuchk(Y, NZ, ASCLE, TOL)
octave_value sin(void) const
subroutine cs1s2(ZR, S1, S2, NZ, ASCLE, ALIM, IUF)
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
subroutine cunik(ZR, FNU, IKFLG, IPMTR, TOL, INIT, PHI, ZETA1, ZETA2, SUM, CWRK)
octave_value cos(void) const
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<