11 #include "tbci/specfun/prototypes.h"
12 #include "tbci/specfun/prototypes2.h"
30 int zunhj_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *ipmtr,
doublereal *tol,
doublereal *phir,
doublereal *phii,
doublereal *argr,
doublereal *argi,
doublereal *zeta1r,
doublereal *zeta1i,
doublereal *zeta2r,
doublereal *zeta2i,
doublereal *asumr,
doublereal *asumi,
doublereal *bsumr,
doublereal *bsumi);
31 int zuoik_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *ikflg,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nuf,
doublereal *tol,
doublereal *elim,
doublereal *alim);
35 integer*,
double*,
double*,
double*,
39 integer*,
double*,
double*,
double*);
42 integer*,
double*,
double*,
double*);
45 integer*,
double*,
double*,
double*);
55 integer*,
double*,
double*,
double*),
79 integer*,
double*,
double*,
double*);
83 extern int zkscl_(
doublereal *zrr,
doublereal *zri,
doublereal *fnu,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
doublereal *rzr,
doublereal *rzi,
doublereal *
ascle,
doublereal *tol,
doublereal *elim);
92 extern int zunhj_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *ipmtr,
doublereal *tol,
doublereal *phir,
doublereal *phii,
doublereal *argr,
doublereal *argi,
doublereal *zeta1r,
doublereal *zeta1i,
doublereal *zeta2r,
doublereal *zeta2i,
doublereal *asumr,
doublereal *asumi,
doublereal *bsumr,
doublereal *bsumi);
94 extern int zbuni_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
integer *nui,
integer *nlast,
doublereal *fnul,
doublereal *tol,
doublereal *elim,
doublereal *alim),
107 zuoik_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *ikflg,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nuf,
doublereal *tol,
doublereal *elim,
doublereal *alim),
108 zwrsk_(
doublereal *zrr,
doublereal *zri,
doublereal *fnu,
integer *kode,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
doublereal *cwr,
doublereal *cwi,
doublereal *tol,
doublereal *elim,
doublereal *alim);
110 extern int zuni1_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
integer *nlast,
doublereal *fnul,
doublereal *tol,
doublereal *elim,
doublereal *alim),
zuni2_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
integer *nlast,
doublereal *fnul,
doublereal *tol,
doublereal *elim,
doublereal *alim);
114 extern int zunk1_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *mr,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
doublereal *tol,
doublereal *elim,
doublereal *alim),
zunk2_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *mr,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
doublereal *tol,
doublereal *elim,
doublereal *alim);
121 extern int zunik_(
doublereal *zrr,
doublereal *zri,
doublereal *fnu,
integer *ikflg,
integer *ipmtr,
doublereal *tol,
integer *init,
doublereal *phir,
doublereal *phii,
doublereal *zeta1r,
doublereal *zeta1i,
doublereal *zeta2r,
doublereal *zeta2i,
doublereal *sumr,
doublereal *sumi,
doublereal *cwrkr,
doublereal *cwrki),
zuoik_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *ikflg,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nuf,
doublereal *tol,
doublereal *elim,
doublereal *alim);
164 static doublereal tol, sti, zni, zti, str, znr;
333 if (*zr == 0. && *zi == 0.) {
339 if (*m < 1 || *m > 2) {
342 if (*kode < 1 || *kode > 2) {
365 tol =
max(d__1,1e-18);
370 i__1 =
abs(k1), i__2 =
abs(k2);
379 alim = elim +
max(d__1,-41.45);
380 fnul = (dig - 3.) * 6. + 10.;
427 aln = -fn *
log(arg);
433 zuoik_(&znr, &zni, fnu, kode, &
c__2, &nn, &cyr[1], &cyi[1], &nuf, &tol, &
448 if (znr < 0. || (znr == 0. && zni < 0. && *m == 2)) {
455 zbknu_(&znr, &zni, fnu, kode, &nn, &cyr[1], &cyi[1], nz, &tol, &elim, &
463 zacon_(&znr, &zni, fnu, kode, &mr, &nn, &cyr[1], &cyi[1], &nw, &rl, &fnul,
475 if (znr >= 0. && (znr != 0. || zni >= 0. || *m != 2)) {
479 if (znr != 0. || zni >= 0.) {
485 zbunk_(&znr, &zni, fnu, kode, &mr, &nn, &cyr[1], &cyi[1], &nw, &tol, &
498 sgn =
d_sign(&hpi, &d__1);
505 ir = inu - (inuh << 1);
510 csgni = rhpi *
cos(arg);
511 csgnr = -rhpi *
sin(arg);
524 for (i__ = 1; i__ <= i__1; ++i__) {
535 d__1 =
abs(aa), d__2 =
abs(bb);
536 if (
max(d__1,d__2) > ascle) {
543 str = aa * csgnr - bb *
csgni;
544 sti = aa * csgni + bb *
csgnr;
545 cyr[i__] = str * atol;
546 cyi[i__] = sti * atol;
584 } equiv_4 = { {2002288515, 1050897, 1487780761, 2146426097,
585 -1209488034, 1017118298, -1209488034, 1018166874, 1352628735,
586 1070810131}, {0}, 0. };
593 #define log10 ((integer *)&equiv_4 + 8)
594 #define dmach ((doublereal *)&equiv_4)
595 #define large ((integer *)&equiv_4 + 2)
596 #define small ((integer *)&equiv_4)
597 #define diver ((integer *)&equiv_4 + 6)
598 #define right ((integer *)&equiv_4 + 4)
1001 if (*i__ < 1 || *i__ > 5) {
1006 ret_val =
dmach[*i__ - 1];
1023 static doublereal gln[100] = { 0.,0.,.693147180559945309,
1024 1.791759469228055,3.17805383034794562,4.78749174278204599,
1025 6.579251212010101,8.5251613610654143,10.6046029027452502,
1026 12.8018274800814696,15.1044125730755153,17.5023078458738858,
1027 19.9872144956618861,22.5521638531234229,25.1912211827386815,
1028 27.8992713838408916,30.6718601060806728,33.5050734501368889,
1029 36.3954452080330536,39.339884187199494,42.335616460753485,
1030 45.380138898476908,48.4711813518352239,51.6066755677643736,
1031 54.7847293981123192,58.0036052229805199,61.261701761002002,
1032 64.5575386270063311,67.889743137181535,71.257038967168009,
1033 74.6582363488301644,78.0922235533153106,81.5579594561150372,
1034 85.0544670175815174,88.5808275421976788,92.1361756036870925,
1035 95.7196945421432025,99.3306124547874269,102.968198614513813,
1036 106.631760260643459,110.320639714757395,114.034211781461703,
1037 117.771881399745072,121.533081515438634,125.317271149356895,
1038 129.123933639127215,132.95257503561631,136.802722637326368,
1039 140.673923648234259,144.565743946344886,148.477766951773032,
1040 152.409592584497358,156.360836303078785,160.331128216630907,
1041 164.320112263195181,168.327445448427652,172.352797139162802,
1042 176.395848406997352,180.456291417543771,184.533828861449491,
1043 188.628173423671591,192.739047287844902,196.866181672889994,
1044 201.009316399281527,205.168199482641199,209.342586752536836,
1045 213.532241494563261,217.736934113954227,221.956441819130334,
1046 226.190548323727593,230.439043565776952,234.701723442818268,
1047 238.978389561834323,243.268849002982714,247.572914096186884,
1048 251.890402209723194,256.221135550009525,260.564940971863209,
1049 264.921649798552801,269.291097651019823,273.673124285693704,
1050 278.067573440366143,282.474292687630396,286.893133295426994,
1051 291.323950094270308,295.766601350760624,300.220948647014132,
1052 304.686856765668715,309.164193580146922,313.652829949879062,
1053 318.152639620209327,322.663499126726177,327.185287703775217,
1054 331.717887196928473,336.261181979198477,340.815058870799018,
1055 345.379407062266854,349.954118040770237,354.539085519440809,
1056 359.134205369575399 };
1057 static doublereal cf[22] = { .0833333333333333333,-.00277777777777777778,
1058 7.93650793650793651e-4,-5.95238095238095238e-4,
1059 8.41750841750841751e-4,-.00191752691752691753,
1060 .00641025641025641026,-.0295506535947712418,.179644372368830573,
1061 -1.39243221690590112,13.402864044168392,-156.848284626002017,
1062 2193.10333333333333,-36108.7712537249894,691472.268851313067,
1063 -15238221.5394074162,382900751.391414141,-10882266035.7843911,
1064 347320283765.002252,-12369602142269.2745,488788064793079.335,
1065 -21320333960919373.9 };
1082 static doublereal fln, tlg, rln, trm, tst, zsq;
1135 fz = *z__ - (
real) nz;
1142 ret_val = gln[nz - 1];
1146 wdtol =
max(wdtol,5e-19);
1152 zm = fln * .3875 + 1.8;
1160 zinc = zmin - (
real) nz;
1171 for (k = 2; k <= 22; ++k) {
1173 trm = cf[k - 1] * zp;
1174 if (
abs(trm) < tst) {
1185 ret_val = *z__ * (tlg - 1.) + (con - tlg) * .5 + s;
1191 for (i__ = 1; i__ <= i__1; ++i__) {
1192 zp *= *z__ + (
real) (i__ - 1);
1196 ret_val = zdmy * (tlg - 1.) -
log(zp) + (con - tlg) * .5 + s;
1212 } equiv_0 = { {5, 6, 0, 0, 32, 4, 2, 31, 2147483647, 2, 24, -125, 127,
1217 static char fmt_9000[] =
"(\0021ERROR 1 IN I1MACH - I OUT OF BOUND\
1226 #define imach ((integer *)&equiv_0)
1227 #define output ((integer *)&equiv_0 + 3)
1230 static cilist io___72 = { 0, 0, 0, fmt_9000, 0 };
1871 if (*i__ < 1 || *i__ > 16) {
1875 ret_val =
imach[*i__ - 1];
1899 static char fmt_900[] =
"(/)";
1907 static integer kmin, i__, k, nn, nr;
1910 static cilist io___76 = { 0, 6, 0, fmt_900, 0 };
1911 static cilist io___79 = { 0, 6, 0, 0, 0 };
1912 static cilist io___80 = { 0, 6, 0, fmt_900, 0 };
1922 nr = *nmess - nn * 70;
1930 for (i__ = 1; i__ <= i__1; ++i__) {
1933 kmin =
min(i__2,*nmess);
1978 ret_val = v *
sqrt(q * q + 1.);
1982 ret_val = u *
sqrt(q * q + 1.);
2045 if (az * az * .25 > dfnu + 1.) {
2052 zseri_(&znr, &zni, fnu, kode, &nn, &yr[1], &yi[1], &nw, tol, elim, alim);
2061 zasyi_(&znr, &zni, fnu, kode, &nn, &yr[1], &yi[1], &nw, rl, tol, elim,
2071 zmlri_(&znr, &zni, fnu, kode, &nn, &yr[1], &yi[1], &nw, tol);
2079 zbknu_(&znr, &zni, fnu, kode, &
c__1, cyr, cyi, &nw, tol, elim, alim);
2084 sgn = -
d_sign(&pi, &fmr);
2117 zs1s2_(&znr, &zni, &c1r, &c1i, &c2r, &c2i, &nw, &
ascle, alim, &iuf);
2131 int zacon_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *mr,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
doublereal *rl,
doublereal *fnul,
doublereal *tol,
doublereal *elim,
doublereal *alim)
2145 static doublereal cscl, cscr, csrr[3], cssr[3], razn;
2151 static doublereal as2, c1r, c2r, s1i, s2i, s1r, s2r, cki,
arg, ckr, cpn;
2153 static doublereal cyi[2], fmr, csr, azn, sgn;
2155 static doublereal bry[3], cyr[2], pti, spn, sti, zni, rzi, ptr, str, znr,
2156 rzr, sc1i, sc2i, sc1r, sc2r;
2182 zbinu_(&znr, &zni, fnu, kode, &nn, &yr[1], &yi[1], &nw, rl, fnul, tol,
2191 zbknu_(&znr, &zni, fnu, kode, &nn, cyr, cyi, &nw, tol, elim, alim);
2198 sgn = -
d_sign(&pi, &fmr);
2207 zmlt_(&csgnr, &csgni, &cpn, &spn, &csgnr, &csgni);
2234 zs1s2_(&znr, &zni, &c1r, &c1i, &c2r, &c2i, &nw, &ascle, alim, &iuf);
2239 zmlt_(&cspnr, &cspni, &c1r, &c1i, &str, &sti);
2240 zmlt_(&csgnr, &csgni, &c2r, &c2i, &ptr, &pti);
2257 zs1s2_(&znr, &zni, &c1r, &c1i, &c2r, &c2i, &nw, &ascle, alim, &iuf);
2262 zmlt_(&cspnr, &cspni, &c1r, &c1i, &str, &sti);
2263 zmlt_(&csgnr, &csgni, &c2r, &c2i, &ptr, &pti);
2275 rzr = (str + str) * razn;
2276 rzi = (sti + sti) * razn;
2292 bry[1] = 1. /
ascle;
2307 bscle = bry[kflag - 1];
2308 s1r *= cssr[kflag - 1];
2309 s1i *= cssr[kflag - 1];
2310 s2r *= cssr[kflag - 1];
2311 s2i *= cssr[kflag - 1];
2312 csr = csrr[kflag - 1];
2314 for (i__ = 3; i__ <= i__1; ++i__) {
2317 s2r = ckr * str - cki * sti + s1r;
2318 s2i = ckr * sti + cki * str + s1i;
2333 zs1s2_(&znr, &zni, &c1r, &c1i, &c2r, &c2i, &nw, &ascle, alim, &iuf);
2343 s1r = sc1r * cssr[kflag - 1];
2344 s1i = sc1i * cssr[kflag - 1];
2345 s2r = sc2r * cssr[kflag - 1];
2346 s2i = sc2i * cssr[kflag - 1];
2350 ptr = cspnr * c1r - cspni * c1i;
2351 pti = cspnr * c1i + cspni * c1r;
2352 yr[i__] = ptr + csgnr * c2r - csgni * c2i;
2353 yi[i__] = pti + csgnr * c2i + csgni * c2r;
2368 bscle = bry[kflag - 1];
2373 s1r *= cssr[kflag - 1];
2374 s1i *= cssr[kflag - 1];
2375 s2r *= cssr[kflag - 1];
2376 s2i *= cssr[kflag - 1];
2377 csr = csrr[kflag - 1];
2400 static doublereal rtpi = .159154943091895336;
2423 static doublereal p1i, s2i, p1r, s2r, cki, dki, fdn,
arg, aez, arm, ckr,
2426 static doublereal raz, czr, ezr, sqk, sti, rzi, tzi, str, rzr, tzr, ak1i,
2427 ak1r, cs1i, cs2i, cs1r, cs2r, dnu2, rtr1;
2458 ak1r = rtpi * str * raz;
2459 ak1i = rtpi * sti * raz;
2460 zsqrt_(&ak1r, &ak1i, &ak1r, &ak1i);
2469 if (
abs(czr) > *elim) {
2474 if (
abs(czr) > *alim && *n > 2) {
2478 zexp_(&czr, &czi, &str, &sti);
2479 zmlt_(&ak1r, &ak1i, &str, &sti, &ak1r, &ak1i);
2506 inu = inu + *n - il;
2521 for (k = 1; k <= i__1; ++k) {
2523 atol = s *
abs(sqk);
2537 for (j = 1; j <= i__2; ++j) {
2538 zdiv_(&ckr, &cki, &dkr, &dki, &str, &sti);
2548 aa = aa *
abs(sqk) / bb;
2561 if (*zr + *zr >= *elim) {
2568 zexp_(&d__1, &d__2, &str, &sti);
2569 zmlt_(&str, &sti, &p1r, &p1i, &str, &sti);
2570 zmlt_(&str, &sti, &cs2r, &cs2i, &str, &sti);
2574 fdn = fdn + dfnu * 8. + 4.;
2578 yr[m] = s2r * ak1r - s2i * ak1i;
2579 yi[m] = s2r * ak1i + s2i * ak1r;
2590 rzr = (str + str) * raz;
2591 rzi = (sti + sti) * raz;
2594 for (i__ = ib; i__ <= i__1; ++i__) {
2595 yr[k] = (ak + *fnu) * (rzr * yr[k + 1] - rzi * yi[k + 1]) + yr[k + 2];
2596 yi[k] = (ak + *fnu) * (rzr * yi[k + 1] + rzi * yr[k + 1]) + yi[k + 2];
2604 zexp_(&czr, &czi, &ckr, &cki);
2606 for (i__ = 1; i__ <= i__1; ++i__) {
2607 str = yr[i__] * ckr - yi[i__] * cki;
2608 yi[i__] = yr[i__] * cki + yi[i__] * ckr;
2621 int zbinu_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *n,
doublereal *cyr,
doublereal *cyi,
integer *nz,
doublereal *rl,
doublereal *fnul,
doublereal *tol,
doublereal *elim,
doublereal *alim)
2659 if (az * az * .25 > dfnu + 1.) {
2666 zseri_(zr, zi, fnu, kode, &nn, &cyr[1], &cyi[1], &nw, tol, elim, alim);
2684 if (az + az < dfnu * dfnu) {
2691 zasyi_(zr, zi, fnu, kode, &nn, &cyr[1], &cyi[1], &nw, rl, tol, elim, alim)
2705 zuoik_(zr, zi, fnu, kode, &
c__1, &nn, &cyr[1], &cyi[1], &nw, tol, elim,
2730 zmlri_(zr, zi, fnu, kode, &nn, &cyr[1], &cyi[1], &nw, tol);
2742 zuoik_(zr, zi, fnu, kode, &
c__2, &
c__2, cwr, cwi, &nw, tol, elim, alim);
2748 for (i__ = 1; i__ <= i__1; ++i__) {
2758 zwrsk_(zr, zi, fnu, kode, &nn, &cyr[1], &cyi[1], &nw, cwr, cwi, tol, elim,
2770 zbuni_(zr, zi, fnu, kode, &nn, &cyr[1], &cyi[1], &nw, &nui, &nlast, fnul,
2808 static doublereal rthpi = 1.25331413731550025;
2813 static doublereal cc[8] = { .577215664901532861,-.0420026350340952355,
2814 -.0421977345555443367,.00721894324666309954,
2815 -2.15241674114950973e-4,-2.01348547807882387e-5,
2816 1.13302723198169588e-6,6.11609510448141582e-9 };
2827 static doublereal cshr, fmui, rcaz, csrr[3], cssr[3], fmur;
2830 static integer i__, j, k, iflag;
2845 static doublereal p1i, p2i, s1i, s2i, p2m, p1r, p2r, s1r, s2r, cbi, cbr,
2846 cki, caz, csi, ckr, fhs, fks, rak, czi, dnu, csr, elm, zdi, bry[3]
2847 , pti, czr, sti, zdr, cyr[2], rzi, ptr, cyi[2];
2879 bry[1] = 1. / bry[0];
2886 sti = -(*zi) * rcaz;
2887 rzr = (str + str) * rcaz;
2888 rzi = (sti + sti) * rcaz;
2891 if (
abs(dnu) == .5) {
2895 if (
abs(dnu) > *tol) {
2905 zlog_(&rzr, &rzi, &smur, &smui, &idum);
2908 zshch_(&fmur, &fmui, &cshr, &cshi, &cchr, &cchi);
2922 t1 = 1. / (t2 * fc);
2923 if (
abs(dnu) > .1) {
2931 for (k = 2; k <= 8; ++k) {
2933 tm = cc[k - 1] * ak;
2935 if (
abs(tm) < *tol) {
2944 g1 = (t1 - t2) / (dnu + dnu);
2946 g2 = (t1 + t2) * .5;
2947 fr = fc * (cchr * g1 + smur * g2);
2948 fi = fc * (cchi * g1 + smui * g2);
2949 zexp_(&fmur, &fmui, &str, &sti);
2964 if (inu > 0 || *n > 1) {
2973 zmlt_(zr, zi, zr, zi, &czr, &czi);
2976 t1 = caz * .25 * caz;
2978 fr = (fr * ak + pr + qr) / bk;
2979 fi = (fi * ak + pi + qi) / bk;
2980 str = 1. / (ak - dnu);
2983 str = 1. / (ak + dnu);
2986 str = ckr * czr - cki * czi;
2988 cki = (ckr * czi + cki * czr) * rak;
2990 s1r = ckr * fr - cki * fi + s1r;
2991 s1i = ckr * fi + cki * fr + s1i;
2993 bk = bk + ak + ak + 1.;
3004 zexp_(zr, zi, &str, &sti);
3005 zmlt_(&s1r, &s1i, &str, &sti, &yr[1], &yi[1]);
3014 zmlt_(zr, zi, zr, zi, &czr, &czi);
3017 t1 = caz * .25 * caz;
3019 fr = (fr * ak + pr + qr) / bk;
3020 fi = (fi * ak + pi + qi) / bk;
3021 str = 1. / (ak - dnu);
3024 str = 1. / (ak + dnu);
3027 str = ckr * czr - cki * czi;
3029 cki = (ckr * czi + cki * czr) * rak;
3031 s1r = ckr * fr - cki * fi + s1r;
3032 s1i = ckr * fi + cki * fr + s1i;
3035 s2r = ckr * str - cki * sti + s2r;
3036 s2i = ckr * sti + cki * str + s2i;
3038 bk = bk + ak + ak + 1.;
3046 ak = a1 *
abs(smur);
3050 str = cssr[kflag - 1];
3053 zmlt_(&p2r, &p2i, &rzr, &rzi, &s2r, &s2i);
3059 zexp_(zr, zi, &fr, &fi);
3060 zmlt_(&s1r, &s1i, &fr, &fi, &s1r, &s1i);
3061 zmlt_(&s2r, &s2i, &fr, &fi, &s2r, &s2i);
3070 zsqrt_(zr, zi, &str, &sti);
3071 zdiv_(&rthpi, &czeroi, &str, &sti, &coefr, &coefi);
3080 str =
exp(-(*zr)) * cssr[kflag - 1];
3081 sti = -str *
sin(*zi);
3083 zmlt_(&coefr, &coefi, &str, &sti, &coefr, &coefi);
3085 if (
abs(dnu) == .5) {
3091 ak =
cos(dpi * dnu);
3096 fhs = (d__1 = .25 - dnu2,
abs(d__1));
3097 if (fhs == czeror) {
3117 t1 =
atan(*zi / *zr);
3126 etest = ak / (dpi * caz * *tol);
3128 if (etest < coner) {
3132 ckr = caz + caz + ctwor;
3136 for (i__ = 1; i__ <= i__1; ++i__) {
3138 cbr = ckr / (fk + coner);
3140 p2r = cbr * p2r - p1r * ak;
3143 fks = fks + fk + fk + ctwor;
3144 fhs = fhs + fk + fk;
3146 str =
abs(p2r) * fk;
3154 fk += spi * t1 *
sqrt(t2 / caz);
3155 fhs = (d__1 = .25 - dnu2,
abs(d__1));
3162 ak = fpi * ak / (*tol *
sqrt(a2));
3163 aa = t1 * 3. / (caz + 1.);
3164 bb = t1 * 14.7 / (caz + 28.);
3165 ak = (
log(ak) + caz *
cos(aa) / (caz * .008 + 1.)) /
cos(bb);
3166 fk = ak * .12125 * ak / caz + 1.5;
3181 for (i__ = 1; i__ <= i__1; ++i__) {
3183 ak = (fks + fk) / (a1 + fhs);
3184 rak = 2. / (fk + coner);
3185 cbr = (fk + *zr) * rak;
3189 p2r = (ptr * cbr - pti * cbi - p1r) * ak;
3190 p2i = (pti * cbr + ptr * cbi - p1i) * ak;
3195 fks = a1 - fk + coner;
3209 zmlt_(&coefr, &coefi, &s1r, &s1i, &str, &sti);
3210 zmlt_(&str, &sti, &csr, &csi, &s1r, &s1i);
3211 if (inu > 0 || *n > 1) {
3230 zmlt_(&p1r, &p1i, &p2r, &p2i, &ptr, &pti);
3231 str = dnu + .5 - ptr;
3233 zdiv_(&str, &sti, zr, zi, &str, &sti);
3235 zmlt_(&str, &sti, &s1r, &s1i, &s2r, &s2i);
3268 p1r = csrr[kflag - 1];
3269 ascle = bry[kflag - 1];
3271 for (i__ = inub; i__ <= i__1; ++i__) {
3274 s2r = ckr * str - cki * sti + s1r;
3275 s2i = ckr * sti + cki * str + s1i;
3292 ascle = bry[kflag - 1];
3297 str = cssr[kflag - 1];
3302 p1r = csrr[kflag - 1];
3312 str = csrr[kflag - 1];
3329 p1r = csrr[kflag - 1];
3330 ascle = bry[kflag - 1];
3332 for (i__ = kk; i__ <= i__1; ++i__) {
3335 s2r = ckr * p2r - cki * p2i + s1r;
3336 s2i = cki * p2r + ckr * p2i + s1i;
3355 ascle = bry[kflag - 1];
3360 str = cssr[kflag - 1];
3365 p1r = csrr[kflag - 1];
3375 elm =
exp(-(*elim));
3383 for (i__ = 1; i__ <= i__1; ++i__) {
3386 s2r = str * ckr - sti * cki + s1r;
3387 s2i = sti * ckr + str * cki + s1i;
3395 if (p2r < -(*elim)) {
3398 zlog_(&s2r, &s2i, &str, &sti, &idum);
3401 p2m =
exp(p2r) / *tol;
3402 p1r = p2m *
cos(p2i);
3403 p1i = p2m *
sin(p2i);
3404 zuchk_(&p1r, &p1i, &nw, &ascle, tol);
3411 if (ic == i__ - 1) {
3461 zkscl_(&zdr, &zdi, fnu, n, &yr[1], &yi[1], nz, &rzr, &rzi, &ascle, tol,
3470 yr[kk] = s1r * csrr[0];
3471 yi[kk] = s1i * csrr[0];
3478 yr[kk] = s2r * csrr[0];
3479 yi[kk] = s2i * csrr[0];
3512 int zbuni_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
integer *nui,
integer *nlast,
doublereal *fnul,
doublereal *tol,
doublereal *elim,
doublereal *alim)
3525 static doublereal c1r, s1i, s2i, s1r, s2r, cyi[2], gnu, raz, cyr[2], sti,
3526 bry[3], rzi, str, rzr;
3546 ax =
abs(*zr) * 1.7321;
3565 zuni1_(zr, zi, &gnu, kode, &
c__2, cyr, cyi, &nw, nlast, fnul, tol, elim,
3574 zuni2_(zr, zi, &gnu, kode, &
c__2, cyr, cyi, &nw, nlast, fnul, tol, elim,
3588 bry[1] = 1. / bry[0];
3609 s1r = cyr[1] * csclr;
3610 s1i = cyi[1] * csclr;
3611 s2r = cyr[0] * csclr;
3612 s2i = cyi[0] * csclr;
3616 rzr = (str + str) * raz;
3617 rzi = (sti + sti) * raz;
3619 for (i__ = 1; i__ <= i__1; ++i__) {
3622 s2r = (dfnu + fnui) * (rzr * str - rzi * sti) + s1r;
3623 s2i = (dfnu + fnui) * (rzr * sti + rzi * str) + s1i;
3639 ascle = bry[iflag - 1];
3653 yr[*n] = s2r * cscrr;
3654 yi[*n] = s2i * cscrr;
3662 for (i__ = 1; i__ <= i__1; ++i__) {
3665 s2r = (*fnu + fnui) * (rzr * str - rzi * sti) + s1r;
3666 s2i = (*fnu + fnui) * (rzr * sti + rzi * str) + s1i;
3685 ascle = bry[iflag - 1];
3714 zuni1_(zr, zi, fnu, kode, n, &yr[1], &yi[1], &nw, nlast, fnul, tol, elim,
3723 zuni2_(zr, zi, fnu, kode, n, &yr[1], &yi[1], &nw, nlast, fnul, tol, elim,
3736 int zbunk_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *mr,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
doublereal *tol,
doublereal *elim,
doublereal *alim)
3756 ax =
abs(*zr) * 1.7321;
3765 zunk1_(zr, zi, fnu, kode, mr, n, &yr[1], &yi[1], nz, tol, elim, alim);
3773 zunk2_(zr, zi, fnu, kode, mr, n, &yr[1], &yi[1], nz, tol, elim, alim);
3792 ca = (*ar * cc + *ai * cd) * bm;
3793 cb = (*ai * cc - *ar * cd) * bm;
3821 int zkscl_(
doublereal *zrr,
doublereal *zri,
doublereal *fnu,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
doublereal *rzr,
doublereal *rzi,
doublereal *
ascle,
doublereal *tol,
doublereal *elim)
3841 static doublereal s1i, s2i, s1r, s2r, acs, cki, elm, csi, ckr, cyi[2],
3842 zdi, csr, cyr[2], zdr, str;
3864 for (i__ = 1; i__ <= i__1; ++i__) {
3870 acs = -(*zrr) +
log(as);
3874 if (acs < -(*elim)) {
3877 zlog_(&s1r, &s1i, &csr, &csi, &idum);
3880 str =
exp(csr) / *tol;
3881 csr = str *
cos(csi);
3882 csi = str *
sin(csi);
3883 zuchk_(&csr, &csi, &nw, ascle, tol);
3918 elm =
exp(-(*elim));
3927 for (i__ = 3; i__ <= i__1; ++i__) {
3931 s2r = ckr * csr - cki * csi + s1r;
3932 s2i = cki * csr + ckr * csi + s1i;
3943 if (acs < -(*elim)) {
3946 zlog_(&s2r, &s2i, &csr, &csi, &idum);
3949 str =
exp(csr) / *tol;
3950 csr = str *
cos(csi);
3951 csi = str *
sin(csi);
3952 zuchk_(&csr, &csi, &nw, ascle, tol);
3985 for (i__ = 1; i__ <= i__1; ++i__) {
3997 static doublereal dpi = 3.141592653589793238462643383;
3998 static doublereal dhpi = 1.570796326794896619231321696;
4020 dtheta =
atan(*ai / *ar);
4085 static integer i__, k, m, itime;
4090 static doublereal p1i, p2i, p1r, p2r, ack, cki, fnf, fkk, ckr;
4094 static doublereal pti, raz, sti, rzi, ptr, str, tst, rzr, rho2;
4115 inu = ifnu + *n - 1;
4120 ckr = str * at * raz;
4121 cki = sti * at * raz;
4122 rzr = (str + str) * raz;
4123 rzi = (sti + sti) * raz;
4128 ack = (at + 1.) * raz;
4129 rho = ack +
sqrt(ack * ack - 1.);
4131 tst = (rho2 + rho2) / ((rho2 - 1.) * (rho - 1.));
4137 for (i__ = 1; i__ <= 80; ++i__) {
4140 p2r = p1r - (ckr * ptr - cki * pti);
4141 p2i = p1i - (cki * ptr + ckr * pti);
4147 if (ap > tst * ak * ak) {
4170 ckr = str * at * raz;
4171 cki = sti * at * raz;
4173 tst =
sqrt(ack / *tol);
4175 for (k = 1; k <= 80; ++k) {
4178 p2r = p1r - (ckr * ptr - cki * pti);
4179 p2i = p1i - (ckr * pti + cki * ptr);
4192 flam = ack +
sqrt(ack * ack - 1.);
4193 fkap = ap /
myzabs_(&p1r, &p1i);
4194 rho =
min(flam,fkap);
4195 tst *=
sqrt(rho / (rho * rho - 1.));
4207 i__1 = i__ + iaz, i__2 = k + inu;
4208 kk =
max(i__1,i__2);
4219 d__1 = fkk + tfnf + 1.;
4229 for (i__ = 1; i__ <= i__1; ++i__) {
4232 p2r = p1r + (fkk + fnf) * (rzr * ptr - rzi * pti);
4233 p2i = p1i + (fkk + fnf) * (rzi * ptr + rzr * pti);
4236 ak = 1. - tfnf / (fkk + tfnf);
4238 sumr += (ack + bk) * p1r;
4239 sumi += (ack + bk) * p1i;
4250 for (i__ = 2; i__ <= i__1; ++i__) {
4253 p2r = p1r + (fkk + fnf) * (rzr * ptr - rzi * pti);
4254 p2i = p1i + (fkk + fnf) * (rzi * ptr + rzr * pti);
4257 ak = 1. - tfnf / (fkk + tfnf);
4259 sumr += (ack + bk) * p1r;
4260 sumi += (ack + bk) * p1i;
4273 for (i__ = 1; i__ <= i__1; ++i__) {
4276 p2r = p1r + (fkk + fnf) * (rzr * ptr - rzi * pti);
4277 p2i = p1i + (fkk + fnf) * (rzr * pti + rzi * ptr);
4280 ak = 1. - tfnf / (fkk + tfnf);
4282 sumr += (ack + bk) * p1r;
4283 sumi += (ack + bk) * p1i;
4294 zlog_(&rzr, &rzi, &str, &sti, &idum);
4295 p1r = -fnf * str + ptr;
4296 p1i = -fnf * sti + pti;
4309 zexp_(&ptr, &pti, &str, &sti);
4314 zmlt_(&ckr, &cki, &ptr, &pti, &cnormr, &cnormi);
4316 for (i__ = 1; i__ <= i__1; ++i__) {
4317 str = yr[i__] * cnormr - yi[i__] * cnormi;
4318 yi[i__] = yr[i__] * cnormi + yi[i__] * cnormr;
4339 ca = *ar * *br - *ai * *bi;
4340 cb = *ar * *bi + *ai * *br;
4372 static doublereal az, cdfnui, cdfnur, ap1, ap2;
4373 static doublereal p1i, p2i, t1i, p1r, p2r, t1r,
arg, rak, rho;
4375 static doublereal pti, tti, rzi, ptr, ttr, rzr, rap1;
4397 idnu = inu + *n - 1;
4401 fnup =
max(amagz,fdnu);
4402 id = idnu - magz - 1;
4406 rzr = ptr * (*zr + *zr) * ptr;
4407 rzi = -ptr * (*zi + *zi) * ptr;
4427 arg = (ap2 + ap2) / (ap1 * *tol);
4441 p2r = p1r - (t1r * ptr - t1i * pti);
4442 p2i = p1i - (t1r * pti + t1i * ptr);
4454 ak =
myzabs_(&t1r, &t1i) * .5;
4455 flam = ak +
sqrt(ak * ak - 1.);
4458 rho =
min(d__1,flam);
4459 test = test1 *
sqrt(rho / (rho * rho - 1.));
4473 for (i__ = 1; i__ <= i__1; ++i__) {
4479 p1r = ptr * ttr - pti * tti + p2r;
4480 p1i = ptr * tti + pti * ttr + p2i;
4486 if (p1r != czeror || p1i != czeroi) {
4492 zdiv_(&p2r, &p2i, &p1r, &p1i, &cyr[*n], &cyi[*n]);
4500 cdfnur = *fnu * rzr;
4501 cdfnui = *fnu * rzi;
4503 for (i__ = 2; i__ <= i__1; ++i__) {
4504 ptr = cdfnur + (t1r * rzr - t1i * rzi) + cyr[k + 1];
4505 pti = cdfnui + (t1r * rzi + t1i * rzr) + cyi[k + 1];
4515 cyr[k] = rak * ptr * rak;
4516 cyi[k] = -rak * pti * rak;
4555 if (*s1r == 0. && *s1i == 0.) {
4561 aln = -(*zrr) - *zrr +
log(as1);
4567 if (aln < -(*alim)) {
4570 zlog_(&s1dr, &s1di, &c1r, &c1i, &idum);
4571 c1r = c1r - *zrr - *zrr;
4572 c1i = c1i - *zri - *zri;
4573 zexp_(&c1r, &c1i, s1r, s1i);
4612 static integer i__, k, l, m, iflag;
4624 static doublereal s1i, s2i, s1r, s2r, cki, acz, arm, ckr, czi, hzi, raz,
4625 czr, sti, hzr, rzi, str, rzr, ak1i, ak1r, rtr1;
4666 zmlt_(&hzr, &hzi, &hzr, &hzi, &czr, &czi);
4670 zlog_(&hzr, &hzi, &ckr, &cki, &idum);
4684 if (ak1r > -(*elim)) {
4700 if (ak1r > -(*alim)) {
4712 coefr = aa *
cos(ak1i);
4713 coefi = aa *
sin(ak1i);
4714 atol = *tol * acz / fnup;
4717 for (i__ = 1; i__ <= i__1; ++i__) {
4722 if (acz < *tol * fnup) {
4732 str = ak1r * czr - ak1i * czi;
4733 sti = ak1r * czi + ak1i * czr;
4745 s2r = s1r * coefr - s1i * coefi;
4746 s2i = s1r * coefi + s1i * coefr;
4752 zuchk_(&s2r, &s2i, &nw, &ascle, tol);
4758 yr[m] = s2r * crscr;
4759 yi[m] = s2i * crscr;
4763 zdiv_(&coefr, &coefi, &hzr, &hzi, &str, &sti);
4777 rzr = (str + str) * raz;
4778 rzi = (sti + sti) * raz;
4785 for (i__ = ib; i__ <= i__1; ++i__) {
4786 yr[k] = (ak + *fnu) * (rzr * yr[k + 1] - rzi * yi[k + 1]) + yr[k + 2];
4787 yi[k] = (ak + *fnu) * (rzr * yi[k + 1] + rzi * yr[k + 1]) + yi[k + 2];
4806 for (l = 3; l <= i__1; ++l) {
4809 s2r = s1r + (ak + *fnu) * (rzr * ckr - rzi * cki);
4810 s2i = s1i + (ak + *fnu) * (rzr * cki + rzi * ckr);
4819 if (
myzabs_(&ckr, &cki) > ascle) {
4849 for (i__ = 2; i__ <= i__1; ++i__) {
4895 static doublereal drt = .7071067811865475244008443621;
4896 static doublereal dpi = 3.141592653589793238462643383;
4917 dtheta =
atan(*ai / *ar);
4952 *br = zm *
cos(dtheta);
4953 *bi = zm *
sin(dtheta);
4999 int zunhj_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *ipmtr,
doublereal *tol,
doublereal *phir,
doublereal *phii,
doublereal *argr,
doublereal *argi,
doublereal *zeta1r,
doublereal *zeta1i,
doublereal *zeta2r,
doublereal *zeta2i,
doublereal *asumr,
doublereal *asumi,
doublereal *bsumr,
doublereal *bsumi)
5003 static doublereal ar[14] = { 1.,.104166666666666667,.0835503472222222222,
5004 .12822657455632716,.291849026464140464,.881627267443757652,
5005 3.32140828186276754,14.9957629868625547,78.9230130115865181,
5006 474.451538868264323,3207.49009089066193,24086.5496408740049,
5007 198923.119169509794,1791902.00777534383 };
5008 static doublereal br[14] = { 1.,-.145833333333333333,
5009 -.0987413194444444444,-.143312053915895062,-.317227202678413548,
5010 -.942429147957120249,-3.51120304082635426,-15.7272636203680451,
5011 -82.2814390971859444,-492.355370523670524,-3316.21856854797251,
5012 -24827.6742452085896,-204526.587315129788,-1838444.9170682099 };
5013 static doublereal c__[105] = { 1.,-.208333333333333333,.125,
5014 .334201388888888889,-.401041666666666667,.0703125,
5015 -1.02581259645061728,1.84646267361111111,-.8912109375,.0732421875,
5016 4.66958442342624743,-11.2070026162229938,8.78912353515625,
5017 -2.3640869140625,.112152099609375,-28.2120725582002449,
5018 84.6362176746007346,-91.8182415432400174,42.5349987453884549,
5019 -7.3687943594796317,.227108001708984375,212.570130039217123,
5020 -765.252468141181642,1059.99045252799988,-699.579627376132541,
5021 218.19051174421159,-26.4914304869515555,.572501420974731445,
5022 -1919.457662318407,8061.72218173730938,-13586.5500064341374,
5023 11655.3933368645332,-5305.64697861340311,1200.90291321635246,
5024 -108.090919788394656,1.7277275025844574,20204.2913309661486,
5025 -96980.5983886375135,192547.001232531532,-203400.177280415534,
5026 122200.46498301746,-41192.6549688975513,7109.51430248936372,
5027 -493.915304773088012,6.07404200127348304,-242919.187900551333,
5028 1311763.6146629772,-2998015.91853810675,3763271.297656404,
5029 -2813563.22658653411,1268365.27332162478,-331645.172484563578,
5030 45218.7689813627263,-2499.83048181120962,24.3805296995560639,
5031 3284469.85307203782,-19706819.1184322269,50952602.4926646422,
5032 -74105148.2115326577,66344512.2747290267,-37567176.6607633513,
5033 13288767.1664218183,-2785618.12808645469,308186.404612662398,
5034 -13886.0897537170405,110.017140269246738,-49329253.664509962,
5035 325573074.185765749,-939462359.681578403,1553596899.57058006,
5036 -1621080552.10833708,1106842816.82301447,-495889784.275030309,
5037 142062907.797533095,-24474062.7257387285,2243768.17792244943,
5038 -84005.4336030240853,551.335896122020586,814789096.118312115,
5039 -5866481492.05184723,18688207509.2958249,-34632043388.1587779,
5040 41280185579.753974,-33026599749.8007231,17954213731.1556001,
5041 -6563293792.61928433,1559279864.87925751,-225105661.889415278,
5042 17395107.5539781645,-549842.327572288687,3038.09051092238427,
5043 -14679261247.6956167,114498237732.02581,-399096175224.466498,
5044 819218669548.577329,-1098375156081.22331,1008158106865.38209,
5045 -645364869245.376503,287900649906.150589,-87867072178.0232657,
5046 17634730606.8349694,-2167164983.22379509,143157876.718888981,
5047 -3871833.44257261262,18257.7554742931747 };
5048 static doublereal alfa[180] = { -.00444444444444444444,
5049 -9.22077922077922078e-4,-8.84892884892884893e-5,
5050 1.65927687832449737e-4,2.4669137274179291e-4,
5051 2.6599558934625478e-4,2.61824297061500945e-4,
5052 2.48730437344655609e-4,2.32721040083232098e-4,
5053 2.16362485712365082e-4,2.00738858762752355e-4,
5054 1.86267636637545172e-4,1.73060775917876493e-4,
5055 1.61091705929015752e-4,1.50274774160908134e-4,
5056 1.40503497391269794e-4,1.31668816545922806e-4,
5057 1.23667445598253261e-4,1.16405271474737902e-4,
5058 1.09798298372713369e-4,1.03772410422992823e-4,
5059 9.82626078369363448e-5,9.32120517249503256e-5,
5060 8.85710852478711718e-5,8.42963105715700223e-5,
5061 8.03497548407791151e-5,7.66981345359207388e-5,
5062 7.33122157481777809e-5,7.01662625163141333e-5,
5063 6.72375633790160292e-5,6.93735541354588974e-4,
5064 2.32241745182921654e-4,-1.41986273556691197e-5,
5065 -1.1644493167204864e-4,-1.50803558053048762e-4,
5066 -1.55121924918096223e-4,-1.46809756646465549e-4,
5067 -1.33815503867491367e-4,-1.19744975684254051e-4,
5068 -1.0618431920797402e-4,-9.37699549891194492e-5,
5069 -8.26923045588193274e-5,-7.29374348155221211e-5,
5070 -6.44042357721016283e-5,-5.69611566009369048e-5,
5071 -5.04731044303561628e-5,-4.48134868008882786e-5,
5072 -3.98688727717598864e-5,-3.55400532972042498e-5,
5073 -3.1741425660902248e-5,-2.83996793904174811e-5,
5074 -2.54522720634870566e-5,-2.28459297164724555e-5,
5075 -2.05352753106480604e-5,-1.84816217627666085e-5,
5076 -1.66519330021393806e-5,-1.50179412980119482e-5,
5077 -1.35554031379040526e-5,-1.22434746473858131e-5,
5078 -1.10641884811308169e-5,-3.54211971457743841e-4,
5079 -1.56161263945159416e-4,3.0446550359493641e-5,
5080 1.30198655773242693e-4,1.67471106699712269e-4,
5081 1.70222587683592569e-4,1.56501427608594704e-4,
5082 1.3633917097744512e-4,1.14886692029825128e-4,
5083 9.45869093034688111e-5,7.64498419250898258e-5,
5084 6.07570334965197354e-5,4.74394299290508799e-5,
5085 3.62757512005344297e-5,2.69939714979224901e-5,
5086 1.93210938247939253e-5,1.30056674793963203e-5,
5087 7.82620866744496661e-6,3.59257485819351583e-6,
5088 1.44040049814251817e-7,-2.65396769697939116e-6,
5089 -4.9134686709848591e-6,-6.72739296091248287e-6,
5090 -8.17269379678657923e-6,-9.31304715093561232e-6,
5091 -1.02011418798016441e-5,-1.0880596251059288e-5,
5092 -1.13875481509603555e-5,-1.17519675674556414e-5,
5093 -1.19987364870944141e-5,3.78194199201772914e-4,
5094 2.02471952761816167e-4,-6.37938506318862408e-5,
5095 -2.38598230603005903e-4,-3.10916256027361568e-4,
5096 -3.13680115247576316e-4,-2.78950273791323387e-4,
5097 -2.28564082619141374e-4,-1.75245280340846749e-4,
5098 -1.25544063060690348e-4,-8.22982872820208365e-5,
5099 -4.62860730588116458e-5,-1.72334302366962267e-5,
5100 5.60690482304602267e-6,2.313954431482868e-5,
5101 3.62642745856793957e-5,4.58006124490188752e-5,
5102 5.2459529495911405e-5,5.68396208545815266e-5,
5103 5.94349820393104052e-5,6.06478527578421742e-5,
5104 6.08023907788436497e-5,6.01577894539460388e-5,
5105 5.891996573446985e-5,5.72515823777593053e-5,
5106 5.52804375585852577e-5,5.3106377380288017e-5,
5107 5.08069302012325706e-5,4.84418647620094842e-5,
5108 4.6056858160747537e-5,-6.91141397288294174e-4,
5109 -4.29976633058871912e-4,1.83067735980039018e-4,
5110 6.60088147542014144e-4,8.75964969951185931e-4,
5111 8.77335235958235514e-4,7.49369585378990637e-4,
5112 5.63832329756980918e-4,3.68059319971443156e-4,
5113 1.88464535514455599e-4,3.70663057664904149e-5,
5114 -8.28520220232137023e-5,-1.72751952869172998e-4,
5115 -2.36314873605872983e-4,-2.77966150694906658e-4,
5116 -3.02079514155456919e-4,-3.12594712643820127e-4,
5117 -3.12872558758067163e-4,-3.05678038466324377e-4,
5118 -2.93226470614557331e-4,-2.77255655582934777e-4,
5119 -2.59103928467031709e-4,-2.39784014396480342e-4,
5120 -2.20048260045422848e-4,-2.00443911094971498e-4,
5121 -1.81358692210970687e-4,-1.63057674478657464e-4,
5122 -1.45712672175205844e-4,-1.29425421983924587e-4,
5123 -1.14245691942445952e-4,.00192821964248775885,
5124 .00135592576302022234,-7.17858090421302995e-4,
5125 -.00258084802575270346,-.00349271130826168475,
5126 -.00346986299340960628,-.00282285233351310182,
5127 -.00188103076404891354,-8.895317183839476e-4,
5128 3.87912102631035228e-6,7.28688540119691412e-4,
5129 .00126566373053457758,.00162518158372674427,.00183203153216373172,
5130 .00191588388990527909,.00190588846755546138,.00182798982421825727,
5131 .0017038950642112153,.00155097127171097686,.00138261421852276159,
5132 .00120881424230064774,.00103676532638344962,
5133 8.71437918068619115e-4,7.16080155297701002e-4,
5134 5.72637002558129372e-4,4.42089819465802277e-4,
5135 3.24724948503090564e-4,2.20342042730246599e-4,
5136 1.28412898401353882e-4,4.82005924552095464e-5 };
5137 static doublereal beta[210] = { .0179988721413553309,
5138 .00559964911064388073,.00288501402231132779,.00180096606761053941,
5139 .00124753110589199202,9.22878876572938311e-4,
5140 7.14430421727287357e-4,5.71787281789704872e-4,
5141 4.69431007606481533e-4,3.93232835462916638e-4,
5142 3.34818889318297664e-4,2.88952148495751517e-4,
5143 2.52211615549573284e-4,2.22280580798883327e-4,
5144 1.97541838033062524e-4,1.76836855019718004e-4,
5145 1.59316899661821081e-4,1.44347930197333986e-4,
5146 1.31448068119965379e-4,1.20245444949302884e-4,
5147 1.10449144504599392e-4,1.01828770740567258e-4,
5148 9.41998224204237509e-5,8.74130545753834437e-5,
5149 8.13466262162801467e-5,7.59002269646219339e-5,
5150 7.09906300634153481e-5,6.65482874842468183e-5,
5151 6.25146958969275078e-5,5.88403394426251749e-5,
5152 -.00149282953213429172,-8.78204709546389328e-4,
5153 -5.02916549572034614e-4,-2.94822138512746025e-4,
5154 -1.75463996970782828e-4,-1.04008550460816434e-4,
5155 -5.96141953046457895e-5,-3.1203892907609834e-5,
5156 -1.26089735980230047e-5,-2.42892608575730389e-7,
5157 8.05996165414273571e-6,1.36507009262147391e-5,
5158 1.73964125472926261e-5,1.9867297884213378e-5,
5159 2.14463263790822639e-5,2.23954659232456514e-5,
5160 2.28967783814712629e-5,2.30785389811177817e-5,
5161 2.30321976080909144e-5,2.28236073720348722e-5,
5162 2.25005881105292418e-5,2.20981015361991429e-5,
5163 2.16418427448103905e-5,2.11507649256220843e-5,
5164 2.06388749782170737e-5,2.01165241997081666e-5,
5165 1.95913450141179244e-5,1.9068936791043674e-5,
5166 1.85533719641636667e-5,1.80475722259674218e-5,
5167 5.5221307672129279e-4,4.47932581552384646e-4,
5168 2.79520653992020589e-4,1.52468156198446602e-4,
5169 6.93271105657043598e-5,1.76258683069991397e-5,
5170 -1.35744996343269136e-5,-3.17972413350427135e-5,
5171 -4.18861861696693365e-5,-4.69004889379141029e-5,
5172 -4.87665447413787352e-5,-4.87010031186735069e-5,
5173 -4.74755620890086638e-5,-4.55813058138628452e-5,
5174 -4.33309644511266036e-5,-4.09230193157750364e-5,
5175 -3.84822638603221274e-5,-3.60857167535410501e-5,
5176 -3.37793306123367417e-5,-3.15888560772109621e-5,
5177 -2.95269561750807315e-5,-2.75978914828335759e-5,
5178 -2.58006174666883713e-5,-2.413083567612802e-5,
5179 -2.25823509518346033e-5,-2.11479656768912971e-5,
5180 -1.98200638885294927e-5,-1.85909870801065077e-5,
5181 -1.74532699844210224e-5,-1.63997823854497997e-5,
5182 -4.74617796559959808e-4,-4.77864567147321487e-4,
5183 -3.20390228067037603e-4,-1.61105016119962282e-4,
5184 -4.25778101285435204e-5,3.44571294294967503e-5,
5185 7.97092684075674924e-5,1.031382367082722e-4,
5186 1.12466775262204158e-4,1.13103642108481389e-4,
5187 1.08651634848774268e-4,1.01437951597661973e-4,
5188 9.29298396593363896e-5,8.40293133016089978e-5,
5189 7.52727991349134062e-5,6.69632521975730872e-5,
5190 5.92564547323194704e-5,5.22169308826975567e-5,
5191 4.58539485165360646e-5,4.01445513891486808e-5,
5192 3.50481730031328081e-5,3.05157995034346659e-5,
5193 2.64956119950516039e-5,2.29363633690998152e-5,
5194 1.97893056664021636e-5,1.70091984636412623e-5,
5195 1.45547428261524004e-5,1.23886640995878413e-5,
5196 1.04775876076583236e-5,8.79179954978479373e-6,
5197 7.36465810572578444e-4,8.72790805146193976e-4,
5198 6.22614862573135066e-4,2.85998154194304147e-4,
5199 3.84737672879366102e-6,-1.87906003636971558e-4,
5200 -2.97603646594554535e-4,-3.45998126832656348e-4,
5201 -3.53382470916037712e-4,-3.35715635775048757e-4,
5202 -3.04321124789039809e-4,-2.66722723047612821e-4,
5203 -2.27654214122819527e-4,-1.89922611854562356e-4,
5204 -1.5505891859909387e-4,-1.2377824076187363e-4,
5205 -9.62926147717644187e-5,-7.25178327714425337e-5,
5206 -5.22070028895633801e-5,-3.50347750511900522e-5,
5207 -2.06489761035551757e-5,-8.70106096849767054e-6,
5208 1.1369868667510029e-6,9.16426474122778849e-6,
5209 1.5647778542887262e-5,2.08223629482466847e-5,
5210 2.48923381004595156e-5,2.80340509574146325e-5,
5211 3.03987774629861915e-5,3.21156731406700616e-5,
5212 -.00180182191963885708,-.00243402962938042533,
5213 -.00183422663549856802,-7.62204596354009765e-4,
5214 2.39079475256927218e-4,9.49266117176881141e-4,
5215 .00134467449701540359,.00148457495259449178,.00144732339830617591,
5216 .00130268261285657186,.00110351597375642682,
5217 8.86047440419791759e-4,6.73073208165665473e-4,
5218 4.77603872856582378e-4,3.05991926358789362e-4,
5219 1.6031569459472163e-4,4.00749555270613286e-5,
5220 -5.66607461635251611e-5,-1.32506186772982638e-4,
5221 -1.90296187989614057e-4,-2.32811450376937408e-4,
5222 -2.62628811464668841e-4,-2.82050469867598672e-4,
5223 -2.93081563192861167e-4,-2.97435962176316616e-4,
5224 -2.96557334239348078e-4,-2.91647363312090861e-4,
5225 -2.83696203837734166e-4,-2.73512317095673346e-4,
5226 -2.6175015580676858e-4,.00638585891212050914,
5227 .00962374215806377941,.00761878061207001043,.00283219055545628054,
5228 -.0020984135201272009,-.00573826764216626498,
5229 -.0077080424449541462,-.00821011692264844401,
5230 -.00765824520346905413,-.00647209729391045177,
5231 -.00499132412004966473,-.0034561228971313328,
5232 -.00201785580014170775,-7.59430686781961401e-4,
5233 2.84173631523859138e-4,.00110891667586337403,
5234 .00172901493872728771,.00216812590802684701,.00245357710494539735,
5235 .00261281821058334862,.00267141039656276912,.0026520307339598043,
5236 .00257411652877287315,.00245389126236094427,.00230460058071795494,
5237 .00213684837686712662,.00195896528478870911,.00177737008679454412,
5238 .00159690280765839059,.00142111975664438546 };
5239 static doublereal gama[30] = { .629960524947436582,.251984209978974633,
5240 .154790300415655846,.110713062416159013,.0857309395527394825,
5241 .0697161316958684292,.0586085671893713576,.0504698873536310685,
5242 .0442600580689154809,.0393720661543509966,.0354283195924455368,
5243 .0321818857502098231,.0294646240791157679,.0271581677112934479,
5244 .0251768272973861779,.0234570755306078891,.0219508390134907203,
5245 .020621082823564624,.0194388240897880846,.0183810633800683158,
5246 .0174293213231963172,.0165685837786612353,.0157865285987918445,
5247 .0150729501494095594,.0144193250839954639,.0138184805735341786,
5248 .0132643378994276568,.0127517121970498651,.0122761545318762767,
5249 .0118338262398482403 };
5254 static doublereal thpi = 4.71238898038468986;
5272 static doublereal zthi, test, tzar, zthr, rfnu2;
5274 static doublereal zetai, ptfni, sumai, sumbi, zetar, ptfnr, razth, sumar,
5279 static integer is, jr, ks, ju;
5285 static doublereal przthi, t2i, w2i, t2r, przthr, w2r, ang, fn13, fn23;
5289 static doublereal zai, zbi, zci, crr[14], drr[14], raw, zar, upi[14], sti,
5290 zbr, zcr, upr[14], str, raw2;
5334 if (
abs(*zr) > ac ||
abs(*zi) > ac) {
5337 *zeta1r = (d__1 =
log(test),
abs(d__1)) * 2. + *fnu;
5349 rfnu2 = rfnu * rfnu;
5353 fn13 =
pow_dd(fnu, &ex1);
5356 w2r = coner - zbr * zbr + zbi * zbi;
5357 w2i = conei - zbr * zbi - zbr * zbi;
5374 for (k = 2; k <= 30; ++k) {
5375 pr[k - 1] = pr[k - 2] * w2r - pi[k - 2] * w2i;
5376 pi[k - 1] = pr[k - 2] * w2i + pi[k - 2] * w2r;
5377 sumar += pr[k - 1] * gama[k - 1];
5378 sumai += pi[k - 1] * gama[k - 1];
5379 ap[k - 1] = ap[k - 2] * aw2;
5380 if (ap[k - 1] < *tol) {
5388 zetar = w2r * sumar - w2i * sumai;
5389 zetai = w2r * sumai + w2i * sumar;
5390 *argr = zetar * fn23;
5391 *argi = zetai * fn23;
5392 zsqrt_(&sumar, &sumai, &zar, &zai);
5393 zsqrt_(&w2r, &w2i, &str, &sti);
5394 *zeta2r = str * *fnu;
5395 *zeta2i = sti * *fnu;
5396 str = coner + ex2 * (zetar * zar - zetai * zai);
5397 sti = conei + ex2 * (zetar * zai + zetai * zar);
5398 *zeta1r = str * *zeta2r - sti * *zeta2i;
5399 *zeta1i = str * *zeta2i + sti * *zeta2r;
5402 zsqrt_(&zar, &zai, &str, &sti);
5403 *phir = str * rfn13;
5404 *phii = sti * rfn13;
5414 for (k = 1; k <= i__1; ++k) {
5415 sumbr += pr[k - 1] * beta[k - 1];
5416 sumbi += pi[k - 1] * beta[k - 1];
5425 btol = *tol * (
abs(*bsumr) +
abs(*bsumi));
5433 for (is = 2; is <= 7; ++is) {
5442 for (k = 1; k <= i__1; ++k) {
5444 sumar += pr[k - 1] * alfa[m - 1];
5445 sumai += pi[k - 1] * alfa[m - 1];
5446 if (ap[k - 1] < atol) {
5452 *asumr += sumar * pp;
5453 *asumi += sumai * pp;
5464 for (k = 1; k <= i__1; ++k) {
5466 sumbr += pr[k - 1] * beta[m - 1];
5467 sumbi += pi[k - 1] * beta[m - 1];
5468 if (ap[k - 1] < atol) {
5474 *bsumr += sumbr * pp;
5475 *bsumi += sumbi * pp;
5480 if (ias == 1 && ibs == 1) {
5498 zsqrt_(&w2r, &w2i, &wr, &wi);
5507 zdiv_(&str, &sti, &zbr, &zbi, &zar, &zai);
5508 zlog_(&zar, &zai, &zcr, &zci, &idum);
5518 zthr = (zcr - wr) * 1.5;
5519 zthi = (zci - wi) * 1.5;
5520 *zeta1r = zcr * *fnu;
5521 *zeta1i = zci * *fnu;
5522 *zeta2r = wr * *fnu;
5523 *zeta2i = wi * *fnu;
5526 if (zthr >= 0. && zthi < 0.) {
5533 ang =
atan(zthi / zthr);
5538 pp =
pow_dd(&azth, &ex2);
5540 zetar = pp *
cos(ang);
5541 zetai = pp *
sin(ang);
5545 *argr = zetar * fn23;
5546 *argi = zetai * fn23;
5547 zdiv_(&zthr, &zthi, &zetar, &zetai, &rtztr, &rtzti);
5548 zdiv_(&rtztr, &rtzti, &wr, &wi, &zar, &zai);
5551 zsqrt_(&tzar, &tzai, &str, &sti);
5552 *phir = str * rfn13;
5553 *phii = sti * rfn13;
5557 raw = 1. /
sqrt(aw2);
5560 tfnr = str * rfnu * raw;
5561 tfni = sti * rfnu * raw;
5564 sti = -zthi * razth;
5565 rzthr = str * razth * rfnu;
5566 rzthi = sti * razth * rfnu;
5567 zcr = rzthr * ar[1];
5568 zci = rzthi * ar[1];
5574 str = t2r * c__[1] + c__[2];
5576 upr[1] = str * tfnr - sti * tfni;
5577 upi[1] = str * tfni + sti * tfnr;
5578 *bsumr = upr[1] + zcr;
5579 *bsumi = upi[1] + zci;
5592 btol = *tol * (
abs(*bsumr) +
abs(*bsumi));
5598 for (lr = 2; lr <= 12; lr += 2) {
5605 for (k = lr; k <= i__1; ++k) {
5612 for (j = 2; j <= i__2; ++j) {
5614 str = zar * t2r - t2i * zai + c__[l - 1];
5615 zai = zar * t2i + zai * t2r;
5619 str = ptfnr * tfnr - ptfni * tfni;
5620 ptfni = ptfnr * tfni + ptfni * tfnr;
5622 upr[kp1 - 1] = ptfnr * zar - ptfni * zai;
5623 upi[kp1 - 1] = ptfni * zar + ptfnr * zai;
5624 crr[ks - 1] = przthr * br[ks];
5625 cri[ks - 1] = przthi * br[ks];
5626 str = przthr * rzthr - przthi * rzthi;
5627 przthi = przthr * rzthi + przthi * rzthr;
5629 drr[ks - 1] = przthr * ar[ks + 1];
5630 dri[ks - 1] = przthi * ar[ks + 1];
5637 sumar = upr[lrp1 - 1];
5638 sumai = upi[lrp1 - 1];
5641 for (jr = 1; jr <= i__1; ++jr) {
5643 sumar = sumar + crr[jr - 1] * upr[ju - 1] - cri[jr - 1] * upi[ju
5645 sumai = sumai + crr[jr - 1] * upi[ju - 1] + cri[jr - 1] * upr[ju
5651 test =
abs(sumar) +
abs(sumai);
5652 if (pp < *tol && test < *tol) {
5659 sumbr = upr[lr + 1] + upr[lrp1 - 1] * zcr - upi[lrp1 - 1] * zci;
5660 sumbi = upi[lr + 1] + upr[lrp1 - 1] * zci + upi[lrp1 - 1] * zcr;
5663 for (jr = 1; jr <= i__1; ++jr) {
5665 sumbr = sumbr + drr[jr - 1] * upr[ju - 1] - dri[jr - 1] * upi[ju
5667 sumbi = sumbi + drr[jr - 1] * upi[ju - 1] + dri[jr - 1] * upr[ju
5673 test =
abs(sumbr) +
abs(sumbi);
5674 if (pp < btol && test < btol) {
5678 if (ias == 1 && ibs == 1) {
5685 str = -(*bsumr) * rfn13;
5686 sti = -(*bsumi) * rfn13;
5687 zdiv_(&str, &sti, &rtztr, &rtzti, bsumr, bsumi);
5691 int zuni1_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
integer *nlast,
doublereal *fnul,
doublereal *tol,
doublereal *elim,
doublereal *alim)
5705 static doublereal aphi, cscl, phii, crsc, phir;
5707 static doublereal csrr[3], cssr[3], rast, sumi, sumr;
5708 static integer i__, k, m, iflag;
5711 static doublereal zeta1i, zeta2i, zeta1r, zeta2r;
5716 static doublereal c1r, c2r, s1i, s2i, rs1, s1r, s2r, cyi[2];
5718 static doublereal bry[3], cyr[2], sti, rzi, str, rzr;
5764 zunik_(zr, zi, &fn, &
c__1, &
c__1, tol, &init, &phir, &phii, &zeta1r, &
5765 zeta1i, &zeta2r, &zeta2i, &sumr, &sumi, cwrkr, cwrki);
5771 rast = fn /
myzabs_(&str, &sti);
5772 str = str * rast * rast;
5773 sti = -sti * rast * rast;
5774 s1r = -zeta1r + str;
5775 s1i = -zeta1i + sti;
5778 s1r = -zeta1r + zeta2r;
5779 s1i = -zeta1i + zeta2i;
5782 if (
abs(rs1) > *elim) {
5788 for (i__ = 1; i__ <= i__1; ++i__) {
5791 zunik_(zr, zi, &fn, &
c__1, &
c__0, tol, &init, &phir, &phii, &zeta1r, &
5792 zeta1i, &zeta2r, &zeta2i, &sumr, &sumi, cwrkr, cwrki);
5798 rast = fn /
myzabs_(&str, &sti);
5799 str = str * rast * rast;
5800 sti = -sti * rast * rast;
5801 s1r = -zeta1r + str;
5802 s1i = -zeta1i + sti + *zi;
5805 s1r = -zeta1r + zeta2r;
5806 s1i = -zeta1i + zeta2i;
5812 if (
abs(rs1) > *elim) {
5818 if (
abs(rs1) < *alim) {
5826 if (
abs(rs1) > *elim) {
5842 s2r = phir * sumr - phii * sumi;
5843 s2i = phir * sumi + phii * sumr;
5844 str =
exp(s1r) * cssr[iflag - 1];
5845 s1r = str *
cos(s1i);
5846 s1i = str *
sin(s1i);
5847 str = s2r * s1r - s2i * s1i;
5848 s2i = s2r * s1i + s2i * s1r;
5853 zuchk_(&s2r, &s2i, &nw, bry, tol);
5861 yr[m] = s2r * csrr[iflag - 1];
5862 yi[m] = s2i * csrr[iflag - 1];
5870 sti = -(*zi) * rast;
5871 rzr = (str + str) * rast;
5872 rzi = (sti + sti) * rast;
5873 bry[1] = 1. / bry[0];
5879 c1r = csrr[iflag - 1];
5880 ascle = bry[iflag - 1];
5884 for (i__ = 3; i__ <= i__1; ++i__) {
5887 s2r = s1r + (*fnu + fn) * (rzr * c2r - rzi * c2i);
5888 s2i = s1i + (*fnu + fn) * (rzr * c2i + rzi * c2r);
5907 ascle = bry[iflag - 1];
5912 s1r *= cssr[iflag - 1];
5913 s1i *= cssr[iflag - 1];
5914 s2r *= cssr[iflag - 1];
5915 s2i *= cssr[iflag - 1];
5916 c1r = csrr[iflag - 1];
5936 zuoik_(zr, zi, fnu, kode, &
c__1, &nd, &yr[1], &yi[1], &nuf, tol, elim,
5961 for (i__ = 1; i__ <= i__1; ++i__) {
5969 int zuni2_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
integer *nlast,
doublereal *fnul,
doublereal *tol,
doublereal *elim,
doublereal *alim)
5976 static doublereal cipr[4] = { 1.,0.,-1.,0. };
5977 static doublereal cipi[4] = { 0.,1.,0.,-1. };
5979 static doublereal aic = 1.265512123484645396;
5989 static doublereal dair, aphi, argi, cscl, phii, crsc, argr;
5991 static doublereal phir, csrr[3], cssr[3], rast;
5992 static integer i__, j, k, iflag;
5995 static doublereal zeta1i, zeta2i, zeta1r, zeta2r;
6000 static doublereal c1r, c2r, s1i, s2i, rs1, s1r, s2r, aii, ang, car;
6004 static doublereal bry[3], raz, sti, zbr, zni, cyr[2], rzi, str, znr, rzr;
6062 str = c2r * cipr[in - 1] - c2i * cipi[in - 1];
6063 c2i = c2r * cipi[in - 1] + c2i * cipr[in - 1];
6077 zunhj_(&znr, &zni, &fn, &
c__1, tol, &phir, &phii, &argr, &argi, &zeta1r, &
6078 zeta1i, &zeta2r, &zeta2i, &asumr, &asumi, &bsumr, &bsumi);
6084 rast = fn /
myzabs_(&str, &sti);
6085 str = str * rast * rast;
6086 sti = -sti * rast * rast;
6087 s1r = -zeta1r + str;
6088 s1i = -zeta1i + sti;
6091 s1r = -zeta1r + zeta2r;
6092 s1i = -zeta1i + zeta2i;
6095 if (
abs(rs1) > *elim) {
6101 for (i__ = 1; i__ <= i__1; ++i__) {
6103 zunhj_(&znr, &zni, &fn, &
c__0, tol, &phir, &phii, &argr, &argi, &
6104 zeta1r, &zeta1i, &zeta2r, &zeta2i, &asumr, &asumi, &bsumr, &
6111 rast = fn /
myzabs_(&str, &sti);
6112 str = str * rast * rast;
6113 sti = -sti * rast * rast;
6114 s1r = -zeta1r + str;
6115 s1i = -zeta1i + sti +
abs(*zi);
6118 s1r = -zeta1r + zeta2r;
6119 s1i = -zeta1i + zeta2i;
6125 if (
abs(rs1) > *elim) {
6131 if (
abs(rs1) < *alim) {
6140 rs1 = rs1 +
log(aphi) -
log(aarg) * .25 - aic;
6141 if (
abs(rs1) > *elim) {
6160 str = dair * bsumr - daii * bsumi;
6161 sti = dair * bsumi + daii * bsumr;
6162 str += air * asumr - aii * asumi;
6163 sti += air * asumi + aii * asumr;
6164 s2r = phir * str - phii * sti;
6165 s2i = phir * sti + phii * str;
6166 str =
exp(s1r) * cssr[iflag - 1];
6167 s1r = str *
cos(s1i);
6168 s1i = str *
sin(s1i);
6169 str = s2r * s1r - s2i * s1i;
6170 s2i = s2r * s1i + s2i * s1r;
6175 zuchk_(&s2r, &s2i, &nw, bry, tol);
6183 str = s2r * c2r - s2i * c2i;
6184 s2i = s2r * c2i + s2i * c2r;
6189 yr[j] = s2r * csrr[iflag - 1];
6190 yi[j] = s2i * csrr[iflag - 1];
6202 rzr = (str + str) * raz;
6203 rzi = (sti + sti) * raz;
6204 bry[1] = 1. / bry[0];
6210 c1r = csrr[iflag - 1];
6211 ascle = bry[iflag - 1];
6215 for (i__ = 3; i__ <= i__1; ++i__) {
6218 s2r = s1r + (*fnu + fn) * (rzr * c2r - rzi * c2i);
6219 s2i = s1i + (*fnu + fn) * (rzr * c2i + rzi * c2r);
6238 ascle = bry[iflag - 1];
6243 s1r *= cssr[iflag - 1];
6244 s1i *= cssr[iflag - 1];
6245 s2r *= cssr[iflag - 1];
6246 s2i *= cssr[iflag - 1];
6247 c1r = csrr[iflag - 1];
6267 zuoik_(zr, zi, fnu, kode, &
c__1, &nd, &yr[1], &yi[1], &nuf, tol, elim,
6292 c2r = car * cipr[in - 1] - sar * cipi[in - 1];
6293 c2i = car * cipi[in - 1] + sar * cipr[in - 1];
6310 for (i__ = 1; i__ <= i__1; ++i__) {
6318 int zunik_(
doublereal *zrr,
doublereal *zri,
doublereal *fnu,
integer *ikflg,
integer *ipmtr,
doublereal *tol,
integer *init,
doublereal *phir,
doublereal *phii,
doublereal *zeta1r,
doublereal *zeta1i,
doublereal *zeta2r,
doublereal *zeta2i,
doublereal *sumr,
doublereal *sumi,
doublereal *cwrkr,
doublereal *cwrki)
6326 static doublereal con[2] = { .398942280401432678,1.25331413731550025 };
6327 static doublereal c__[120] = { 1.,-.208333333333333333,.125,
6328 .334201388888888889,-.401041666666666667,.0703125,
6329 -1.02581259645061728,1.84646267361111111,-.8912109375,.0732421875,
6330 4.66958442342624743,-11.2070026162229938,8.78912353515625,
6331 -2.3640869140625,.112152099609375,-28.2120725582002449,
6332 84.6362176746007346,-91.8182415432400174,42.5349987453884549,
6333 -7.3687943594796317,.227108001708984375,212.570130039217123,
6334 -765.252468141181642,1059.99045252799988,-699.579627376132541,
6335 218.19051174421159,-26.4914304869515555,.572501420974731445,
6336 -1919.457662318407,8061.72218173730938,-13586.5500064341374,
6337 11655.3933368645332,-5305.64697861340311,1200.90291321635246,
6338 -108.090919788394656,1.7277275025844574,20204.2913309661486,
6339 -96980.5983886375135,192547.001232531532,-203400.177280415534,
6340 122200.46498301746,-41192.6549688975513,7109.51430248936372,
6341 -493.915304773088012,6.07404200127348304,-242919.187900551333,
6342 1311763.6146629772,-2998015.91853810675,3763271.297656404,
6343 -2813563.22658653411,1268365.27332162478,-331645.172484563578,
6344 45218.7689813627263,-2499.83048181120962,24.3805296995560639,
6345 3284469.85307203782,-19706819.1184322269,50952602.4926646422,
6346 -74105148.2115326577,66344512.2747290267,-37567176.6607633513,
6347 13288767.1664218183,-2785618.12808645469,308186.404612662398,
6348 -13886.0897537170405,110.017140269246738,-49329253.664509962,
6349 325573074.185765749,-939462359.681578403,1553596899.57058006,
6350 -1621080552.10833708,1106842816.82301447,-495889784.275030309,
6351 142062907.797533095,-24474062.7257387285,2243768.17792244943,
6352 -84005.4336030240853,551.335896122020586,814789096.118312115,
6353 -5866481492.05184723,18688207509.2958249,-34632043388.1587779,
6354 41280185579.753974,-33026599749.8007231,17954213731.1556001,
6355 -6563293792.61928433,1559279864.87925751,-225105661.889415278,
6356 17395107.5539781645,-549842.327572288687,3038.09051092238427,
6357 -14679261247.6956167,114498237732.02581,-399096175224.466498,
6358 819218669548.577329,-1098375156081.22331,1008158106865.38209,
6359 -645364869245.376503,287900649906.150589,-87867072178.0232657,
6360 17634730606.8349694,-2167164983.22379509,143157876.718888981,
6361 -3871833.44257261262,18257.7554742931747,286464035717.679043,
6362 -2406297900028.50396,9109341185239.89896,-20516899410934.4374,
6363 30565125519935.3206,-31667088584785.1584,23348364044581.8409,
6364 -12320491305598.2872,4612725780849.13197,-1196552880196.1816,
6365 205914503232.410016,-21822927757.5292237,1247009293.51271032,
6366 -29188388.1222208134,118838.426256783253 };
6379 static doublereal ac, si, ti, sr, tr, t2i, t2r, rfn, sri, sti, zni, srr,
6422 if (
abs(*zrr) > ac ||
abs(*zri) > ac) {
6425 *zeta1r = (d__1 =
log(test),
abs(d__1)) * 2. + *fnu;
6435 sr = coner + (tr * tr - ti * ti);
6436 si = conei + (tr * ti + ti * tr);
6437 zsqrt_(&sr, &si, &srr, &sri);
6440 zdiv_(&str, &sti, &tr, &ti, &znr, &zni);
6441 zlog_(&znr, &zni, &str, &sti, &idum);
6442 *zeta1r = *fnu * str;
6443 *zeta1i = *fnu * sti;
6444 *zeta2r = *fnu * srr;
6445 *zeta2i = *fnu * sri;
6446 zdiv_(&coner, &conei, &srr, &sri, &tr, &ti);
6449 zsqrt_(&srr, &sri, &cwrkr[16], &cwrki[16]);
6450 *phir = cwrkr[16] * con[*ikflg - 1];
6451 *phii = cwrki[16] * con[*ikflg - 1];
6455 zdiv_(&coner, &conei, &sr, &si, &t2r, &t2i);
6462 for (k = 2; k <= 15; ++k) {
6466 for (j = 1; j <= i__1; ++j) {
6468 str = sr * t2r - si * t2i + c__[l - 1];
6469 si = sr * t2i + si * t2r;
6473 str = crfnr * srr - crfni * sri;
6474 crfni = crfnr * sri + crfni * srr;
6476 cwrkr[k] = crfnr * sr - crfni * si;
6477 cwrki[k] = crfnr * si + crfni * sr;
6479 test = (d__1 = cwrkr[k],
abs(d__1)) + (d__2 = cwrki[k],
abs(d__2));
6480 if (ac < *tol && test < *tol) {
6498 for (i__ = 1; i__ <= i__1; ++i__) {
6505 *phir = cwrkr[16] * con[0];
6506 *phii = cwrki[16] * con[0];
6516 for (i__ = 1; i__ <= i__1; ++i__) {
6517 sr += tr * cwrkr[i__];
6518 si += tr * cwrki[i__];
6524 *phir = cwrkr[16] * con[1];
6525 *phii = cwrki[16] * con[1];
6529 int zunk1_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *mr,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
doublereal *tol,
doublereal *elim,
doublereal *alim)
6544 static doublereal aphi, cscl, phii[2], crsc, phir[2];
6546 static doublereal csrr[3], cssr[3], rast, sumi[2], razr;
6548 static integer i__, j, k, m, iflag, kflag;
6557 static doublereal zeta1i[2], zeta2i[2], zet1di, zet2di, zeta1r[2], zeta2r[
6563 static doublereal c1r, c2r, s1i, s2i, rs1, s1r, s2r, ang, asc, cki, fnf;
6569 static doublereal bry[3], cyr[2], sti, rzi, zri, str, rzr, zrr;
6605 bry[1] = 1. / bry[0];
6617 for (i__ = 1; i__ <= i__1; ++i__) {
6624 zunik_(&zrr, &zri, &fn, &
c__2, &
c__0, tol, &init[j - 1], &phir[j - 1],
6625 &phii[j - 1], &zeta1r[j - 1], &zeta1i[j - 1], &zeta2r[j - 1],
6626 &zeta2i[j - 1], &sumr[j - 1], &sumi[j - 1], &cwrkr[(j << 4)
6627 - 16], &cwrki[(j << 4) - 16]);
6631 str = zrr + zeta2r[j - 1];
6632 sti = zri + zeta2i[j - 1];
6633 rast = fn /
myzabs_(&str, &sti);
6634 str = str * rast * rast;
6635 sti = -sti * rast * rast;
6636 s1r = zeta1r[j - 1] - str;
6637 s1i = zeta1i[j - 1] - sti;
6640 s1r = zeta1r[j - 1] - zeta2r[j - 1];
6641 s1i = zeta1i[j - 1] - zeta2i[j - 1];
6647 if (
abs(rs1) > *elim) {
6653 if (
abs(rs1) < *alim) {
6659 aphi =
myzabs_(&phir[j - 1], &phii[j - 1]);
6661 if (
abs(rs1) > *elim) {
6678 s2r = phir[j - 1] * sumr[j - 1] - phii[j - 1] * sumi[j - 1];
6679 s2i = phir[j - 1] * sumi[j - 1] + phii[j - 1] * sumr[j - 1];
6680 str =
exp(s1r) * cssr[kflag - 1];
6681 s1r = str *
cos(s1i);
6682 s1i = str *
sin(s1i);
6683 str = s2r * s1r - s2i * s1i;
6684 s2i = s1r * s2i + s2r * s1i;
6689 zuchk_(&s2r, &s2i, &nw, bry, tol);
6694 cyr[kdflg - 1] = s2r;
6695 cyi[kdflg - 1] = s2i;
6696 yr[i__] = s2r * csrr[kflag - 1];
6697 yi[i__] = s2i * csrr[kflag - 1];
6720 if (yr[i__ - 1] == zeror && yi[i__ - 1] == zeroi) {
6723 yr[i__ - 1] = zeror;
6724 yi[i__ - 1] = zeroi;
6731 razr = 1. /
myzabs_(&zrr, &zri);
6734 rzr = (str + str) * razr;
6735 rzi = (sti + sti) * razr;
6752 zunik_(&zrr, &zri, &fn, &
c__2, &ipard, tol, &initd, &phidr, &phidi, &
6753 zet1dr, &zet1di, &zet2dr, &zet2di, &sumdr, &sumdi, &cwrkr[32], &
6760 rast = fn /
myzabs_(&str, &sti);
6761 str = str * rast * rast;
6762 sti = -sti * rast * rast;
6767 s1r = zet1dr - zet2dr;
6768 s1i = zet1di - zet2di;
6771 if (
abs(rs1) > *elim) {
6774 if (
abs(rs1) < *alim) {
6780 aphi =
myzabs_(&phidr, &phidi);
6782 if (
abs(rs1) < *elim) {
6786 if (
abs(rs1) > 0.) {
6797 for (i__ = 1; i__ <= i__1; ++i__) {
6811 c1r = csrr[kflag - 1];
6812 ascle = bry[kflag - 1];
6814 for (i__ = ib; i__ <= i__1; ++i__) {
6817 s2r = ckr * c2r - cki * c2i + s1r;
6818 s2i = ckr * c2i + cki * c2r + s1i;
6837 ascle = bry[kflag - 1];
6842 s1r *= cssr[kflag - 1];
6843 s1i *= cssr[kflag - 1];
6844 s2r *= cssr[kflag - 1];
6845 s2i *= cssr[kflag - 1];
6846 c1r = csrr[kflag - 1];
6859 sgn = -
d_sign(&pi, &fmr);
6883 for (k = 1; k <= i__1; ++k) {
6894 initd = init[j - 1];
6895 phidr = phir[j - 1];
6896 phidi = phii[j - 1];
6897 zet1dr = zeta1r[j - 1];
6898 zet1di = zeta1i[j - 1];
6899 zet2dr = zeta2r[j - 1];
6900 zet2di = zeta2i[j - 1];
6901 sumdr = sumr[j - 1];
6902 sumdi = sumi[j - 1];
6907 if (kk == *n && ib < *n) {
6910 if (kk == ib || kk == ic) {
6915 zunik_(&zrr, &zri, &fn, &
c__1, &
c__0, tol, &initd, &phidr, &phidi, &
6916 zet1dr, &zet1di, &zet2dr, &zet2di, &sumdr, &sumdi, &cwrkr[(m
6917 << 4) - 16], &cwrki[(m << 4) - 16]);
6923 rast = fn /
myzabs_(&str, &sti);
6924 str = str * rast * rast;
6925 sti = -sti * rast * rast;
6926 s1r = -zet1dr + str;
6927 s1i = -zet1di + sti;
6930 s1r = -zet1dr + zet2dr;
6931 s1i = -zet1di + zet2di;
6937 if (
abs(rs1) > *elim) {
6943 if (
abs(rs1) < *alim) {
6949 aphi =
myzabs_(&phidr, &phidi);
6951 if (
abs(rs1) > *elim) {
6964 str = phidr * sumdr - phidi * sumdi;
6965 sti = phidr * sumdi + phidi * sumdr;
6968 str =
exp(s1r) * cssr[iflag - 1];
6969 s1r = str *
cos(s1i);
6970 s1i = str *
sin(s1i);
6971 str = s2r * s1r - s2i * s1i;
6972 s2i = s2r * s1i + s2i * s1r;
6977 zuchk_(&s2r, &s2i, &nw, bry, tol);
6984 cyr[kdflg - 1] = s2r;
6985 cyi[kdflg - 1] = s2i;
6988 s2r *= csrr[iflag - 1];
6989 s2i *= csrr[iflag - 1];
6998 zs1s2_(&zrr, &zri, &s1r, &s1i, &s2r, &s2i, &nw, &asc, alim, &iuf);
7001 yr[kk] = s1r * cspnr - s1i * cspni + s2r;
7002 yi[kk] = cspnr * s1i + cspni * s1r + s2i;
7006 if (c2r != 0. || c2i != 0.) {
7042 csr = csrr[iflag - 1];
7043 ascle = bry[iflag - 1];
7046 for (i__ = 1; i__ <= i__1; ++i__) {
7049 s2r = s1r + (fn + fnf) * (rzr * c2r - rzi * c2i);
7050 s2i = s1i + (fn + fnf) * (rzr * c2i + rzi * c2r);
7063 zs1s2_(&zrr, &zri, &c1r, &c1i, &c2r, &c2i, &nw, &asc, alim, &iuf);
7066 yr[kk] = c1r * cspnr - c1i * cspni + c2r;
7067 yi[kk] = c1r * cspni + c1i * cspnr + c2i;
7081 ascle = bry[iflag - 1];
7086 s1r *= cssr[iflag - 1];
7087 s1i *= cssr[iflag - 1];
7088 s2r *= cssr[iflag - 1];
7089 s2i *= cssr[iflag - 1];
7090 csr = csrr[iflag - 1];
7100 int zunk2_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *mr,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
doublereal *tol,
doublereal *elim,
doublereal *alim)
7108 static doublereal cr1i = 1.73205080756887729;
7110 static doublereal cr2i = -.866025403784438647;
7114 static doublereal cipr[4] = { 1.,0.,-1.,0. };
7115 static doublereal cipi[4] = { 0.,-1.,0.,1. };
7125 static doublereal dair, aphi, argi[2], cscl, phii[2], crsc, argr[2];
7127 static doublereal phir[2], csrr[3], cssr[3], rast, razr;
7128 static integer i__, k, j, iflag, kflag;
7135 static doublereal zeta1i[2], zeta2i[2], zet1di, zet2di, zeta1r[2], zeta2r[
7139 static integer il, kk, in, nw;
7140 static doublereal asumdi, bsumdi, yy, asumdr, bsumdr, c1i, c2i, c2m;
7141 static doublereal c1r, c2r, s1i, s2i, rs1, s1r, s2r, aii, ang, asc, car,
7148 static doublereal cyi[2], fmr, sar, csr, sgn, zbi;
7150 static doublereal bry[3], cyr[2], pti, sti, zbr, zni, rzi, ptr, zri, str,
7191 bry[1] = 1. / bry[0];
7214 str = c2r * cipr[kk - 1] - c2i * cipi[kk - 1];
7215 sti = c2r * cipi[kk - 1] + c2i * cipr[kk - 1];
7216 csr = cr1r * str - cr1i * sti;
7217 csi = cr1r * sti + cr1i * str;
7231 for (i__ = 1; i__ <= i__1; ++i__) {
7237 zunhj_(&znr, &zni, &fn, &
c__0, tol, &phir[j - 1], &phii[j - 1], &argr[
7238 j - 1], &argi[j - 1], &zeta1r[j - 1], &zeta1i[j - 1], &zeta2r[
7239 j - 1], &zeta2i[j - 1], &asumr[j - 1], &asumi[j - 1], &bsumr[
7240 j - 1], &bsumi[j - 1]);
7244 str = zbr + zeta2r[j - 1];
7245 sti = zbi + zeta2i[j - 1];
7246 rast = fn /
myzabs_(&str, &sti);
7247 str = str * rast * rast;
7248 sti = -sti * rast * rast;
7249 s1r = zeta1r[j - 1] - str;
7250 s1i = zeta1i[j - 1] - sti;
7253 s1r = zeta1r[j - 1] - zeta2r[j - 1];
7254 s1i = zeta1i[j - 1] - zeta2i[j - 1];
7260 if (
abs(rs1) > *elim) {
7266 if (
abs(rs1) < *alim) {
7272 aphi =
myzabs_(&phir[j - 1], &phii[j - 1]);
7273 aarg =
myzabs_(&argr[j - 1], &argi[j - 1]);
7274 rs1 = rs1 +
log(aphi) -
log(aarg) * .25 - aic;
7275 if (
abs(rs1) > *elim) {
7292 c2r = argr[j - 1] * cr2r - argi[j - 1] * cr2i;
7293 c2i = argr[j - 1] * cr2i + argi[j - 1] * cr2r;
7296 str = dair * bsumr[j - 1] - daii * bsumi[j - 1];
7297 sti = dair * bsumi[j - 1] + daii * bsumr[j - 1];
7298 ptr = str * cr2r - sti * cr2i;
7299 pti = str * cr2i + sti * cr2r;
7300 str = ptr + (air * asumr[j - 1] - aii * asumi[j - 1]);
7301 sti = pti + (air * asumi[j - 1] + aii * asumr[j - 1]);
7302 ptr = str * phir[j - 1] - sti * phii[j - 1];
7303 pti = str * phii[j - 1] + sti * phir[j - 1];
7304 s2r = ptr * csr - pti * csi;
7305 s2i = ptr * csi + pti * csr;
7306 str =
exp(s1r) * cssr[kflag - 1];
7307 s1r = str *
cos(s1i);
7308 s1i = str *
sin(s1i);
7309 str = s2r * s1r - s2i * s1i;
7310 s2i = s1r * s2i + s2r * s1i;
7315 zuchk_(&s2r, &s2i, &nw, bry, tol);
7323 cyr[kdflg - 1] = s2r;
7324 cyi[kdflg - 1] = s2i;
7325 yr[i__] = s2r * csrr[kflag - 1];
7326 yi[i__] = s2i * csrr[kflag - 1];
7355 if (yr[i__ - 1] == zeror && yi[i__ - 1] == zeroi) {
7358 yr[i__ - 1] = zeror;
7359 yi[i__ - 1] = zeroi;
7366 razr = 1. /
myzabs_(&zrr, &zri);
7369 rzr = (str + str) * razr;
7370 rzi = (sti + sti) * razr;
7386 zunhj_(&znr, &zni, &fn, &ipard, tol, &phidr, &phidi, &argdr, &argdi, &
7387 zet1dr, &zet1di, &zet2dr, &zet2di, &asumdr, &asumdi, &bsumdr, &
7394 rast = fn /
myzabs_(&str, &sti);
7395 str = str * rast * rast;
7396 sti = -sti * rast * rast;
7401 s1r = zet1dr - zet2dr;
7402 s1i = zet1di - zet2di;
7405 if (
abs(rs1) > *elim) {
7408 if (
abs(rs1) < *alim) {
7414 aphi =
myzabs_(&phidr, &phidi);
7416 if (
abs(rs1) < *elim) {
7431 for (i__ = 1; i__ <= i__1; ++i__) {
7442 c1r = csrr[kflag - 1];
7443 ascle = bry[kflag - 1];
7445 for (i__ = ib; i__ <= i__1; ++i__) {
7448 s2r = ckr * c2r - cki * c2i + s1r;
7449 s2i = ckr * c2i + cki * c2r + s1i;
7468 ascle = bry[kflag - 1];
7473 s1r *= cssr[kflag - 1];
7474 s1i *= cssr[kflag - 1];
7475 s2r *= cssr[kflag - 1];
7476 s2i *= cssr[kflag - 1];
7477 c1r = csrr[kflag - 1];
7490 sgn = -
d_sign(&pi, &fmr);
7519 str = csr * c2r + csi * c2i;
7520 csi = -csr * c2i + csi * c2r;
7529 for (k = 1; k <= i__1; ++k) {
7539 phidr = phir[j - 1];
7540 phidi = phii[j - 1];
7541 argdr = argr[j - 1];
7542 argdi = argi[j - 1];
7543 zet1dr = zeta1r[j - 1];
7544 zet1di = zeta1i[j - 1];
7545 zet2dr = zeta2r[j - 1];
7546 zet2di = zeta2i[j - 1];
7547 asumdr = asumr[j - 1];
7548 asumdi = asumi[j - 1];
7549 bsumdr = bsumr[j - 1];
7550 bsumdi = bsumi[j - 1];
7554 if (kk == *n && ib < *n) {
7557 if (kk == ib || kk == ic) {
7560 zunhj_(&znr, &zni, &fn, &
c__0, tol, &phidr, &phidi, &argdr, &argdi, &
7561 zet1dr, &zet1di, &zet2dr, &zet2di, &asumdr, &asumdi, &bsumdr,
7569 rast = fn /
myzabs_(&str, &sti);
7570 str = str * rast * rast;
7571 sti = -sti * rast * rast;
7572 s1r = -zet1dr + str;
7573 s1i = -zet1di + sti;
7576 s1r = -zet1dr + zet2dr;
7577 s1i = -zet1di + zet2di;
7583 if (
abs(rs1) > *elim) {
7589 if (
abs(rs1) < *alim) {
7595 aphi =
myzabs_(&phidr, &phidi);
7596 aarg =
myzabs_(&argdr, &argdi);
7597 rs1 = rs1 +
log(aphi) -
log(aarg) * .25 - aic;
7598 if (
abs(rs1) > *elim) {
7612 zairy_(&argdr, &argdi, &
c__1, &
c__2, &dair, &daii, &ndai, &idum);
7613 str = dair * bsumdr - daii * bsumdi;
7614 sti = dair * bsumdi + daii * bsumdr;
7615 str += air * asumdr - aii * asumdi;
7616 sti += air * asumdi + aii * asumdr;
7617 ptr = str * phidr - sti * phidi;
7618 pti = str * phidi + sti * phidr;
7619 s2r = ptr * csr - pti * csi;
7620 s2i = ptr * csi + pti * csr;
7621 str =
exp(s1r) * cssr[iflag - 1];
7622 s1r = str *
cos(s1i);
7623 s1i = str *
sin(s1i);
7624 str = s2r * s1r - s2i * s1i;
7625 s2i = s2r * s1i + s2i * s1r;
7630 zuchk_(&s2r, &s2i, &nw, bry, tol);
7640 cyr[kdflg - 1] = s2r;
7641 cyi[kdflg - 1] = s2i;
7644 s2r *= csrr[iflag - 1];
7645 s2i *= csrr[iflag - 1];
7654 zs1s2_(&zrr, &zri, &s1r, &s1i, &s2r, &s2i, &nw, &asc, alim, &iuf);
7657 yr[kk] = s1r * cspnr - s1i * cspni + s2r;
7658 yi[kk] = s1r * cspni + s1i * cspnr + s2i;
7665 if (c2r != 0. || c2i != 0.) {
7701 csr = csrr[iflag - 1];
7702 ascle = bry[iflag - 1];
7705 for (i__ = 1; i__ <= i__1; ++i__) {
7708 s2r = s1r + (fn + fnf) * (rzr * c2r - rzi * c2i);
7709 s2i = s1i + (fn + fnf) * (rzr * c2i + rzi * c2r);
7722 zs1s2_(&zrr, &zri, &c1r, &c1i, &c2r, &c2i, &nw, &asc, alim, &iuf);
7725 yr[kk] = c1r * cspnr - c1i * cspni + c2r;
7726 yi[kk] = c1r * cspni + c1i * cspnr + c2i;
7740 ascle = bry[iflag - 1];
7745 s1r *= cssr[iflag - 1];
7746 s1i *= cssr[iflag - 1];
7747 s2r *= cssr[iflag - 1];
7748 s2i *= cssr[iflag - 1];
7749 csr = csrr[iflag - 1];
7759 int zuoik_(
doublereal *zr,
doublereal *zi,
doublereal *fnu,
integer *kode,
integer *ikflg,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nuf,
doublereal *tol,
doublereal *elim,
doublereal *alim)
7765 static doublereal aic = 1.265512123484645396;
7773 static doublereal aarg, aphi, argi, phii, argr;
7783 static doublereal zeta1i, zeta2i, zeta1r, zeta2r, ax, ay;
7785 static doublereal fnn, gnn, zbi, czi, gnu, zbr, czr, rcz, sti, zni, zri,
7832 ax =
abs(*zr) * 1.7321;
7843 gnn = *fnu + fnn - 1.;
7855 zunik_(&zrr, &zri, &gnu, ikflg, &
c__1, tol, &init, &phir, &phii, &zeta1r,
7856 &zeta1i, &zeta2r, &zeta2i, &sumr, &sumi, cwrkr, cwrki);
7857 czr = -zeta1r + zeta2r;
7858 czi = -zeta1i + zeta2i;
7868 zunhj_(&znr, &zni, &gnu, &
c__1, tol, &phir, &phii, &argr, &argi, &zeta1r,
7869 &zeta1i, &zeta2r, &zeta2i, &asumr, &asumi, &bsumr, &bsumi);
7870 czr = -zeta1r + zeta2r;
7871 czi = -zeta1i + zeta2i;
7899 rcz = rcz -
log(aarg) * .25 - aic;
7909 if (rcz < -(*elim)) {
7912 if (rcz > -(*alim)) {
7917 rcz = rcz -
log(aarg) * .25 - aic;
7919 if (rcz > -(*elim)) {
7924 for (i__ = 1; i__ <= i__1; ++i__) {
7933 zlog_(&phir, &phii, &str, &sti, &idum);
7939 zlog_(&argr, &argi, &str, &sti, &idum);
7940 czr = czr - str * .25 - aic;
7943 ax =
exp(rcz) / *tol;
7947 zuchk_(&czr, &czi, &nw, &ascle, tol);
7967 zunik_(&zrr, &zri, &gnu, ikflg, &
c__1, tol, &init, &phir, &phii, &zeta1r,
7968 &zeta1i, &zeta2r, &zeta2i, &sumr, &sumi, cwrkr, cwrki);
7969 czr = -zeta1r + zeta2r;
7970 czi = -zeta1i + zeta2i;
7973 zunhj_(&znr, &zni, &gnu, &
c__1, tol, &phir, &phii, &argr, &argi, &zeta1r,
7974 &zeta1i, &zeta2r, &zeta2i, &asumr, &asumi, &bsumr, &bsumi);
7975 czr = -zeta1r + zeta2r;
7976 czi = -zeta1i + zeta2i;
7987 if (rcz < -(*elim)) {
7990 if (rcz > -(*alim)) {
7995 rcz = rcz -
log(aarg) * .25 - aic;
7997 if (rcz > -(*elim)) {
8011 zlog_(&phir, &phii, &str, &sti, &idum);
8017 zlog_(&argr, &argi, &str, &sti, &idum);
8018 czr = czr - str * .25 - aic;
8021 ax =
exp(rcz) / *tol;
8025 zuchk_(&czr, &czi, &nw, &ascle, tol);
8035 int zwrsk_(
doublereal *zrr,
doublereal *zri,
doublereal *fnu,
integer *kode,
integer *n,
doublereal *yr,
doublereal *yi,
integer *nz,
doublereal *cwr,
doublereal *cwi,
doublereal *tol,
doublereal *elim,
doublereal *alim)
8048 static doublereal c1r, c2r, act, acw, cti, ctr, pti, sti, ptr, str;
8072 zbknu_(zrr, zri, fnu, kode, &
c__2, &cwr[1], &cwi[1], &nw, tol, elim, alim)
8077 zrati_(zrr, zri, fnu, n, &yr[1], &yi[1], tol);
8096 acw =
myzabs_(&cwr[2], &cwi[2]);
8111 c1r = cwr[1] * csclr;
8112 c1i = cwi[1] * csclr;
8113 c2r = cwr[2] * csclr;
8114 c2i = cwi[2] * csclr;
8121 ptr = str * c1r - sti * c1i;
8122 pti = str * c1i + sti * c1r;
8125 ctr = *zrr * ptr - *zri * pti;
8126 cti = *zrr * pti + *zri * ptr;
8133 cinur = ptr * ctr - pti * cti;
8134 cinui = ptr * cti + pti * ctr;
8135 yr[1] = cinur * csclr;
8136 yi[1] = cinui * csclr;
8141 for (i__ = 2; i__ <= i__1; ++i__) {
8142 ptr = str * cinur - sti * cinui;
8143 cinui = str * cinui + sti * cinur;
8147 yr[i__] = cinur * csclr;
8148 yi[i__] = cinui * csclr;
double d_sign(const doublereal *, const doublereal *)
int zdiv_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci)
int zmlt_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci)
int zlog_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, integer *ierr)
integer do_lio(integer *, integer *, char *, ftnlen)
int zrati_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *n, doublereal *cyr, doublereal *cyi, doublereal *tol)
TBCI__ cplx< T > log(const TBCI__ cplx< T > &z)
int zasyi_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *rl, doublereal *tol, doublereal *elim, doublereal *alim)
TBCI__ cplx< T > cos(const TBCI__ cplx< T > &z)
int zacai_(double *, double *, double *, const integer *, integer *, integer *, double *, double *, integer *, double *, double *, double *, double *)
int zbknu_(double *, double *, double *, const integer *, integer *, double *, double *, integer *, double *, double *, double *)
int zbesh_(double *zr, double *zi, double *fnu, int *kode, int *m, int *n, double *cyr, double *cyi, int *nz, int *ierr)
int zsqrt_(const double *, const double *, double *, double *)
int zs1s2_(doublereal *zrr, doublereal *zri, doublereal *s1r, doublereal *s1i, doublereal *s2r, doublereal *s2i, integer *nz, doublereal *ascle, doublereal *alim, integer *iuf)
TBCI__ cplx< T > atan(const TBCI__ cplx< T > &z)
int zairy_(const doublereal *zr, const doublereal *zi, const integer *id, const integer *kode, doublereal *air, doublereal *aii, integer *nz, integer *ierr)
int zunhj_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *ipmtr, doublereal *tol, doublereal *phir, doublereal *phii, doublereal *argr, doublereal *argi, doublereal *zeta1r, doublereal *zeta1i, doublereal *zeta2r, doublereal *zeta2i, doublereal *asumr, doublereal *asumi, doublereal *bsumr, doublereal *bsumi)
int zunk1_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *mr, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *tol, doublereal *elim, doublereal *alim)
int zwrsk_(doublereal *zrr, doublereal *zri, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *cwr, doublereal *cwi, doublereal *tol, doublereal *elim, doublereal *alim)
T max(const FS_Vector< dims, T > &fv)
int zunik_(doublereal *zrr, doublereal *zri, doublereal *fnu, integer *ikflg, integer *ipmtr, doublereal *tol, integer *init, doublereal *phir, doublereal *phii, doublereal *zeta1r, doublereal *zeta1i, doublereal *zeta2r, doublereal *zeta2i, doublereal *sumr, doublereal *sumi, doublereal *cwrkr, doublereal *cwrki)
T arg(const TBCI__ cplx< T > &c)
doublereal myzabs_(const double *, const double *)
int zmlri_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *tol)
TBCI__ cplx< T > sin(const TBCI__ cplx< T > &z)
int zacon_(double *, double *, double *, integer *, integer *, integer *, double *, double *, integer *, double *, double *, double *, double *, double *)
int s_stop(const char *, ftnlen)
int zseri_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *tol, doublereal *elim, doublereal *alim)
TBCI__ cplx< T > exp(const TBCI__ cplx< T > &z)
int integer
barf [ba:rf] 2.
int zexp_(const double *, const double *, double *, double *)
int zuoik_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *ikflg, integer *n, doublereal *yr, doublereal *yi, integer *nuf, doublereal *tol, doublereal *elim, doublereal *alim)
int zuni2_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, integer *nlast, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim)
int zbinu_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *cyr, doublereal *cyi, integer *nz, doublereal *rl, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim)
int zuni1_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, integer *nlast, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim)
double pow_dd(const doublereal *, const doublereal *)
int zkscl_(doublereal *zrr, doublereal *zri, doublereal *fnu, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *rzr, doublereal *rzi, doublereal *ascle, doublereal *tol, doublereal *elim)
doublereal dgamln_(doublereal *z__, integer *ierr)
int zuchk_(doublereal *yr, doublereal *yi, integer *nz, doublereal *ascle, doublereal *tol)
int zunk2_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *mr, integer *n, doublereal *yr, doublereal *yi, integer *nz, doublereal *tol, doublereal *elim, doublereal *alim)
int xerror_(char *mess, integer *nmess, integer *l1, integer *l2, ftnlen mess_len)
doublereal d1mach_(const integer *)
int zshch_(doublereal *zr, doublereal *zi, doublereal *cshr, doublereal *cshi, doublereal *cchr, doublereal *cchi)
int zbunk_(double *, double *, double *, integer *, integer *, integer *, double *, double *, integer *, double *, double *, double *)
TBCI__ cplx< T > sinh(const TBCI__ cplx< T > &z)
int zbuni_(doublereal *zr, doublereal *zi, doublereal *fnu, integer *kode, integer *n, doublereal *yr, doublereal *yi, integer *nz, integer *nui, integer *nlast, doublereal *fnul, doublereal *tol, doublereal *elim, doublereal *alim)
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > long int int char v
< find minimun of func on grid with resolution res
TBCI__ cplx< T > cosh(const TBCI__ cplx< T > &z)
integer i1mach_(const integer *)