1 SUBROUTINE cuni2(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM,
18 COMPLEX AI, ARG, ASUM, BSUM, CFN, CI, CID, CIP, CONE, CRSC, CSCL,
19 * csr, css, cy, czero, c1, c2, dai, phi, rz, s1, s2, y, z, zb,
20 * zeta1, zeta2, zn, zar
21 REAL AARG, AIC, ALIM, ANG, APHI, ASCLE, AY, BRY, CAR, C2I, C2M,
22 * c2r, elim, fn, fnu, fnul, hpi, rs1, sar, tol, yy, r1mach
23 INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
24 * nn, nuf, nw, nz, idum
25 dimension bry(3), y(n), cip(4), css(3), csr(3), cy(2)
26 DATA czero,cone,ci/(0.0e0,0.0e0),(1.0e0,0.0e0),(0.0e0,1.0e0)/
27 DATA cip(1),cip(2),cip(3),cip(4)/
28 1 (1.0e0,0.0e0), (0.0e0,1.0e0), (-1.0e0,0.0e0), (0.0e0,-1.0e0)/
30 1 1.57079632679489662
e+00, 1.265512123484645396
e+00/
40 cscl =
cmplx(1.0e0/tol,0.0e0)
41 crsc =
cmplx(tol,0.0e0)
48 bry(1) = 1.0
e+3*r1mach(1)/tol
57 ang = hpi*(fnu-float(inu))
65 IF (yy.GT.0.0e0) go
to 10
75 CALL
cunhj(zn, fn, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
76 IF (kode.EQ.1) go
to 20
77 cfn =
cmplx(fnu,0.0e0)
78 s1 = -zeta1 + cfn*(cfn/(zb+zeta2))
84 IF (
abs(rs1).GT.elim) go
to 150
88 fn = fnu + float(nd-i)
89 CALL
cunhj(zn, fn, 0, tol, phi, arg, zeta1, zeta2, asum, bsum)
90 IF (kode.EQ.1) go
to 50
93 s1 = -zeta1 + cfn*(cfn/(zb+zeta2)) +
cmplx(0.0e0,ay)
102 IF (
abs(rs1).GT.elim) go
to 120
103 IF (i.EQ.1) iflag = 2
104 IF (
abs(rs1).LT.alim) go
to 70
111 rs1 = rs1 + alog(aphi) - 0.25e0*alog(aarg) - aic
112 IF (
abs(rs1).GT.elim) go
to 120
113 IF (i.EQ.1) iflag = 1
114 IF (rs1.LT.0.0e0) go
to 70
115 IF (i.EQ.1) iflag = 3
121 CALL
cairy(arg, 0, 2, ai, nai, idum)
122 CALL
cairy(arg, 1, 2, dai, ndai, idum)
123 s2 = phi*(ai*asum+dai*bsum)
126 c2m =
exp(c2r)*
REAL(css(iflag))
129 IF (iflag.NE.1) go
to 80
130 CALL
cuchk(s2, nw, bry(1), tol)
131 IF (nw.NE.0) go
to 120
133 IF (yy.LE.0.0e0) s2 = conjg(s2)
140 IF (nd.LE.2) go
to 110
141 rz =
cmplx(2.0e0,0.0e0)/z
142 bry(2) = 1.0e0/bry(1)
152 s2 = s1 +
cmplx(fnu+fn,0.0e0)*rz*s2
158 IF (iflag.GE.3) go
to 100
164 IF (c2m.LE.ascle) go
to 100
176 IF (rs1.GT.0.0e0) go
to 140
183 IF (nd.EQ.0) go
to 110
184 CALL
cuoik(z, fnu, kode, 1, nd, y, nuf, tol, elim, alim)
185 IF (nuf.LT.0) go
to 140
188 IF (nd.EQ.0) go
to 110
189 fn = fnu + float(nd-1)
190 IF (fn.LT.fnul) go
to 130
200 IF (yy.LE.0.0e0)c2=conjg(c2)
209 IF (rs1.GT.0.0e0) go
to 140
subroutine cairy(Z, ID, KODE, AI, NZ, IERR)
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
subroutine cuni2(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM, ALIM)
subroutine cuchk(Y, NZ, ASCLE, TOL)
octave_value sin(void) const
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 cuoik(Z, FNU, KODE, IKFLG, N, Y, NUF, TOL, ELIM, ALIM)
subroutine cunhj(Z, FNU, IPMTR, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM)
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))<