1 SUBROUTINE cbknu(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
10 COMPLEX CCH, CK, COEF, CONE, CRSC, CS, CSCL, CSH, CSR, CSS, CTWO,
11 * cz, czero,
f, fmu,
p, pt, p1, p2, q, rz, smu, st, s1, s2, y, z,
13 REAL AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ, CC, DNU,
14 * dnu2, elim, etest, fc, fhs, fk, fks, fnu, fpi, g1, g2, hpi, pi,
15 * p2i, p2m, p2r, rk, rthpi, r1,
s, spi, tm, tol, tth, t1, t2, xx,
16 * yy,
gamln, r1mach, helim, elm, xd, yd, alas, as
17 INTEGER I, IDUM, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N,
18 * nz, i1mach, nw, j, ic, inub
19 dimension bry(3), cc(8), css(3), csr(3), y(n), cy(2)
23 DATA czero,cone,ctwo /(0.0e0,0.0e0),(1.0e0,0.0e0),(2.0e0,0.0e0)/
25 DATA pi, rthpi, spi ,hpi, fpi, tth /
26 1 3.14159265358979324e0, 1.25331413731550025e0,
27 2 1.90985931710274403e0, 1.57079632679489662e0,
28 3 1.89769999331517738e0, 6.66666666666666666
e-01/
30 DATA cc(1), cc(2), cc(3), cc(4), cc(5), cc(6), cc(7), cc(8)/
31 1 5.77215664901532861
e-01, -4.20026350340952355
e-02,
32 2 -4.21977345555443367
e-02, 7.21894324666309954
e-03,
33 3 -2.15241674114950973
e-04, -2.01348547807882387
e-05,
34 4 1.13302723198169588
e-06, 6.11609510448141582
e-09/
39 cscl =
cmplx(1.0e0/tol,0.0e0)
40 crsc =
cmplx(tol,0.0e0)
47 bry(1) = 1.0
e+3*r1mach(1)/tol
55 dnu = fnu - float(inu)
56 IF (
abs(dnu).EQ.0.5e0) go
to 110
58 IF (
abs(dnu).GT.tol) dnu2 = dnu*dnu
59 IF (caz.GT.r1) go
to 110
65 fmu = smu*
cmplx(dnu,0.0e0)
66 CALL
cshch(fmu, csh, cch)
67 IF (dnu.EQ.0.0e0) go
to 10
70 smu = csh*
cmplx(1.0e0/dnu,0.0e0)
78 IF (
abs(dnu).GT.0.1e0) go
to 40
88 IF (
abs(tm).LT.tol) go
to 30
93 g1 = (t1-t2)/(dnu+dnu)
99 p =
cmplx(0.5e0/t2,0.0e0)*pt
100 q =
cmplx(0.5e0/t1,0.0e0)/pt
107 IF (inu.GT.0 .OR. n.GT.1) go
to 80
111 IF (caz.LT.tol) go
to 70
112 cz = z*z*
cmplx(0.25e0,0.0e0)
116 p =
p*
cmplx(1.0e0/(ak-dnu),0.0e0)
117 q = q*
cmplx(1.0e0/(ak+dnu),0.0e0)
119 ck = ck*cz*
cmplx(rk,0.0)
122 bk = bk + ak + ak + 1.0e0
124 IF (a1.GT.tol) go
to 60
127 IF (koded.EQ.1)
RETURN
134 IF (caz.LT.tol) go
to 100
135 cz = z*z*
cmplx(0.25e0,0.0e0)
139 p =
p*
cmplx(1.0e0/(ak-dnu),0.0e0)
140 q = q*
cmplx(1.0e0/(ak+dnu),0.0e0)
142 ck = ck*cz*
cmplx(rk,0.0e0)
144 s2 = s2 + ck*(
p-
f*
cmplx(ak,0.0e0))
146 bk = bk + ak + ak + 1.0e0
148 IF (a1.GT.tol) go
to 90
154 IF (ak.GT.alim) kflag = 3
158 IF (koded.EQ.1) go
to 210
170 coef =
cmplx(rthpi,0.0e0)/csqrt(z)
172 IF (koded.EQ.2) go
to 120
173 IF (xx.GT.alim) go
to 290
175 a1 =
exp(-xx)*
REAL(css(kflag))
179 IF (
abs(dnu).EQ.0.5e0) go
to 300
185 IF (ak.EQ.0.0e0) go
to 300
186 fhs =
abs(0.25e0-dnu2)
187 IF (fhs.EQ.0.0e0) go
to 300
194 t1 = float(i1mach(11)-1)*r1mach(5)*3.321928094e0
195 t1 = amax1(t1,12.0e0)
196 t1 = amin1(t1,60.0e0)
198 IF (xx.NE.0.0e0) go
to 130
205 IF (t2.GT.caz) go
to 170
209 etest = ak/(pi*caz*tol)
211 IF (etest.LT.1.0e0) go
to 180
213 rk = caz + caz + 2.0e0
223 fks = fks + fk + fk + 2.0e0
227 IF (etest.LT.tm) go
to 160
231 fk = fk + spi*t1*
sqrt(t2/caz)
232 fhs =
abs(0.25e0-dnu2)
239 ak = fpi*ak/(tol*
sqrt(a2))
240 aa = 3.0e0*t1/(1.0e0+caz)
241 bb = 14.7e0*t1/(28.0e0+caz)
242 ak = (alog(ak)+caz*
cos(aa)/(1.0e0+0.008e0*caz))/
cos(bb)
243 fk = 0.12125e0*ak*ak/caz + 1.5e0
252 p2 =
cmplx(tol,0.0e0)
256 a2 = (fks+fk)/(a1+fhs)
257 rk = 2.0e0/(fk+1.0e0)
264 fks = a1 - fk + 1.0e0
272 pt =
cmplx(1.0e0/tm,0.0e0)
276 IF (inu.GT.0 .OR. n.GT.1) go
to 200
278 IF(iflag.EQ.1) go
to 270
285 pt =
cmplx(1.0e0/tm,0.0e0)
289 s2 = s1*(cone+(
cmplx(dnu+0.5e0,0.0e0)-pt)/z)
295 ck =
cmplx(dnu+1.0e0,0.0e0)*rz
296 IF (n.EQ.1) inu = inu - 1
297 IF (inu.GT.0) go
to 220
300 IF(iflag.EQ.1) go
to 270
304 IF (iflag.EQ.1) go
to 261
313 IF (kflag.GE.3) go
to 230
320 IF (p2m.LE.ascle) go
to 230
348 IF (kflag.GE.3) go
to 260
354 IF (p2m.LE.ascle) go
to 260
370 celm =
cmplx(elm,0.0)
385 IF(p2r.LT.(-elim)) go
to 263
391 CALL
cuchk(p1,nw,ascle,tol)
392 IF(nw.NE.0) go
to 263
395 IF(ic.EQ.(i-1)) go
to 264
399 IF(alas.LT.helim) go
to 262
413 IF(inub.LE.inu) go
to 225
418 IF (n.EQ.1) go
to 280
422 CALL
ckscl(zd, fnu, n, y, nz, rz, ascle, tol, elim)
433 t2 = fnu + float(kk-1)
434 ck =
cmplx(t2,0.0e0)*rz
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE * f
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)
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))<
subroutine ckscl(ZR, FNU, N, Y, NZ, RZ, ASCLE, TOL, ELIM)
subroutine cshch(Z, CSH, CCH)
subroutine cbknu(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
octave_value sqrt(void) const