1 SUBROUTINE zbknu(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM,
12 DOUBLE PRECISION AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ,
13 * cbi, cbr, cc, cchi, cchr, cki, ckr, coefi, coefr, conei, coner,
14 * crscr, csclr, cshi, cshr, csi, csr, csrr, cssr, ctwor,
15 * czeroi, czeror, czi, czr, dnu, dnu2, dpi, elim, etest, fc, fhs,
16 *
fi, fk, fks, fmui, fmur, fnu, fpi, fr, g1, g2, hpi, pi, pr, pti,
17 * ptr, p1i, p1r, p2i, p2m, p2r, qi, qr, rak, rcaz, rthpi, rzi,
18 * rzr, r1, s, smui, smur, spi, sti,
str, s1i, s1r, s2i, s2r, tm,
19 * tol, tth, t1, t2, yi, yr, zi, zr,
dgamln,
d1mach, xzabs, elm,
20 * celmr, zdr, zdi, as, alas, helim, cyr, cyi
21 INTEGER I, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N, NZ,
22 * idum, i1mach, j, ic, inub, nw
23 dimension yr(n), yi(n), cc(8), cssr(3), csrr(3), bry(3), cyr(2),
29 DATA czeror,czeroi,coner,conei,ctwor,r1/
30 1 0.0d0 , 0.0d0 , 1.0d0 , 0.0d0 , 2.0d0 , 2.0d0 /
31 DATA dpi, rthpi, spi ,hpi, fpi, tth /
32 1 3.14159265358979324d0, 1.25331413731550025d0,
33 2 1.90985931710274403d0, 1.57079632679489662d0,
34 3 1.89769999331517738d0, 6.66666666666666666
d-01/
35 DATA cc(1), cc(2), cc(3), cc(4), cc(5), cc(6), cc(7), cc(8)/
36 1 5.77215664901532861
d-01, -4.20026350340952355
d-02,
37 2 -4.21977345555443367
d-02, 7.21894324666309954
d-03,
38 3 -2.15241674114950973
d-04, -2.01348547807882387
d-05,
39 4 1.13302723198169588
d-06, 6.11609510448141582
d-09/
61 inu = int(sngl(fnu+0.5d0))
62 dnu = fnu - dble(float(inu))
63 IF (dabs(dnu).EQ.0.5d0) go to 110
65 IF (dabs(dnu).GT.tol) dnu2 = dnu*dnu
66 IF (caz.GT.r1) go to 110
71 CALL
xzlog(rzr, rzi, smur, smui, idum)
74 CALL
zshch(fmur, fmui, cshr, cshi, cchr, cchi)
75 IF (dnu.EQ.0.0d0) go to 10
85 t2 = dexp(-
dgamln(a2,idum))
87 IF (dabs(dnu).GT.0.1d0) go to 40
97 IF (dabs(tm).LT.tol) go to 30
102 g1 = (t1-t2)/(dnu+dnu)
105 fr = fc*(cchr*g1+smur*g2)
106 fi = fc*(cchi*g1+smui*g2)
110 CALL
zdiv(0.5d0, 0.0d0,
str, sti, ptr, pti)
122 IF (inu.GT.0 .OR. n.GT.1) go to 80
126 IF (caz.LT.tol) go to 70
127 CALL
zmlt(zr, zi, zr, zi, czr, czi)
132 fr = (fr*ak+pr+qr)/bk
133 fi = (
fi*ak+pi+qi)/bk
140 str = ckr*czr - cki*czi
142 cki = (ckr*czi+cki*czr)*rak
144 s1r = ckr*fr - cki*
fi + s1r
145 s1i = ckr*
fi + cki*fr + s1i
147 bk = bk + ak + ak + 1.0d0
149 IF (a1.GT.tol) go to 60
153 IF (koded.EQ.1)
RETURN
155 CALL
zmlt(s1r, s1i,
str, sti, yr(1), yi(1))
161 IF (caz.LT.tol) go to 100
162 CALL
zmlt(zr, zi, zr, zi, czr, czi)
167 fr = (fr*ak+pr+qr)/bk
168 fi = (
fi*ak+pi+qi)/bk
175 str = ckr*czr - cki*czi
177 cki = (ckr*czi+cki*czr)*rak
179 s1r = ckr*fr - cki*
fi + s1r
180 s1i = ckr*
fi + cki*fr + s1i
183 s2r = ckr*
str - cki*sti + s2r
184 s2i = ckr*sti + cki*
str + s2i
186 bk = bk + ak + ak + 1.0d0
188 IF (a1.GT.tol) go to 90
193 IF (ak.GT.alim) kflag = 3
197 CALL
zmlt(p2r, p2i, rzr, rzi, s2r, s2i)
200 IF (koded.EQ.1) go to 210
202 CALL
zmlt(s1r, s1i, fr,
fi, s1r, s1i)
203 CALL
zmlt(s2r, s2i, fr,
fi, s2r, s2i)
213 CALL
zdiv(rthpi, czeroi,
str, sti, coefr, coefi)
215 IF (koded.EQ.2) go to 120
216 IF (zr.GT.alim) go to 290
218 str = dexp(-zr)*cssr(kflag)
221 CALL
zmlt(coefr, coefi,
str, sti, coefr, coefi)
223 IF (dabs(dnu).EQ.0.5d0) go to 300
229 IF (ak.EQ.czeror) go to 300
230 fhs = dabs(0.25d0-dnu2)
231 IF (fhs.EQ.czeror) go to 300
238 t1 = dble(float(i1mach(14)-1))
239 t1 = t1*
d1mach(5)*3.321928094d0
240 t1 = dmax1(t1,12.0d0)
241 t1 = dmin1(t1,60.0d0)
243 IF (zr.NE.0.0d0) go to 130
250 IF (t2.GT.caz) go to 170
254 etest = ak/(dpi*caz*tol)
256 IF (etest.LT.coner) go to 180
258 ckr = caz + caz + ctwor
265 p2r = cbr*p2r - p1r*ak
268 fks = fks + fk + fk + ctwor
272 IF (etest.LT.
str) go to 160
276 fk = fk + spi*t1*dsqrt(t2/caz)
277 fhs = dabs(0.25d0-dnu2)
284 ak = fpi*ak/(tol*dsqrt(a2))
285 aa = 3.0d0*t1/(1.0d0+caz)
286 bb = 14.7d0*t1/(28.0d0+caz)
287 ak = (dlog(ak)+caz*dcos(aa)/(1.0d0+0.008d0*caz))/dcos(bb)
288 fk = 0.12125d0*ak*ak/caz + 1.5d0
304 ak = (fks+fk)/(a1+fhs)
305 rak = 2.0d0/(fk+coner)
310 p2r = (ptr*cbr-pti*cbi-p1r)*ak
311 p2i = (pti*cbr+ptr*cbi-p1i)*ak
316 fks = a1 - fk + coner
329 CALL
zmlt(coefr, coefi, s1r, s1i,
str, sti)
330 CALL
zmlt(
str, sti, csr, csi, s1r, s1i)
331 IF (inu.GT.0 .OR. n.GT.1) go to 200
334 IF(iflag.EQ.1) go to 270
346 CALL
zmlt(p1r, p1i, p2r, p2i, ptr, pti)
347 str = dnu + 0.5d0 - ptr
351 CALL
zmlt(
str, sti, s1r, s1i, s2r, s2i)
360 IF (n.EQ.1) inu = inu - 1
361 IF (inu.GT.0) go to 220
362 IF (n.GT.1) go to 215
368 IF(iflag.EQ.1) go to 270
372 IF(iflag.EQ.1) go to 261
379 s2r = ckr*
str - cki*sti + s1r
380 s2i = ckr*sti + cki*
str + s1i
385 IF (kflag.GE.3) go to 230
391 IF (p2m.LE.ascle) go to 230
405 IF (n.NE.1) go to 240
425 s2r = ckr*p2r - cki*p2i + s1r
426 s2i = cki*p2r + ckr*p2i + s1i
435 IF (kflag.GE.3) go to 260
439 IF (p2m.LE.ascle) go to 260
469 s2r =
str*ckr-sti*cki+s1r
470 s2i = sti*ckr+
str*cki+s1i
478 IF(p2r.LT.(-elim)) go to 263
485 CALL
zuchk(p1r,p1i,nw,ascle,tol)
486 IF(nw.NE.0) go to 263
490 IF(ic.EQ.(i-1)) go to 264
494 IF(alas.LT.helim) go to 262
513 IF(inub.LE.inu) go to 225
526 CALL
zkscl(zdr,zdi,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim)
541 t2 = fnu + dble(float(kk-1))
std::string str(char sep= 'x') const
subroutine xzlog(AR, AI, BR, BI, IERR)
subroutine xzexp(AR, AI, BR, BI)
std::string dimension(void) const
double precision function d1mach(i)
subroutine zkscl(ZRR, ZRI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM)
F77_RET_T const double const double double * d
subroutine zdiv(AR, AI, BR, BI, CR, CI)
subroutine xzsqrt(AR, AI, BR, BI)
double precision function dgamln(Z, IERR)
subroutine zshch(ZR, ZI, CSHR, CSHI, CCHR, CCHI)
subroutine zmlt(AR, AI, BR, BI, CR, CI)
subroutine zuchk(YR, YI, NZ, ASCLE, TOL)
subroutine zbknu(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM)