3982 int zbknu(
double *zr,
double *zi,
double *fnu,
3983 long *kode,
long *
n,
double *yr,
double *yi,
3984 long *nz,
double *tol,
double *elim,
3996 double dpi = 3.14159265358979324;
3997 double rthpi = 1.25331413731550025;
3998 double spi = 1.90985931710274403;
3999 double hpi = 1.57079632679489662;
4000 double fpi = 1.89769999331517738;
4001 double tth = .666666666666666666;
4002 double cc[8] = { .577215664901532861, -.0420026350340952355,
4003 -.0421977345555443367, .00721894324666309954,
4004 -2.15241674114950973e-4, -2.01348547807882387e-5,
4005 1.13302723198169588e-6, 6.11609510448141582e-9
4013 double sin(
double), exp(
double), cos(
double), atan(
double)
4014 , sqrt(
double), log(
double);
4018 double s, a1, a2, g1, g2, t1, t2, aa, bb, fc, ak, bk;
4022 double fr, pi, qi, tm, pr, qr;
4024 double p1i, p2i, s1i, s2i, p2m, p1r, p2r, s1r, s2r, cbi, cbr,
4025 cki, caz, csi, ckr, fhs, fks, rak, czi, dnu, csr, elm, zdi, bry[3]
4026 , pti, czr, sti, zdr, cyr[2], rzi, ptr, cyi[2];
4028 double str, rzr, dnu2, cchi, cchr, alas, cshi;
4030 double cshr, fmui, rcaz, csrr[3], cssr[3], fmur;
4031 extern int zdiv(
double *,
double *,
double *
4032 ,
double *,
double *,
double *);
4034 extern int zmlt(
double *,
double *,
double *
4035 ,
double *,
double *,
double *);
4039 double ascle, coefr, helim;
4040 extern double azabs(
double *,
double *);
4041 double celmr, csclr, crscr;
4042 extern int azlog(
double *,
double *,
double
4043 *,
double *,
long *),
zshch(
double *,
4051 extern int zuchk(
double *,
double *,
long *,
4052 double *,
double *),
azexp(
double *,
4056 zkscl(
double *,
double *,
double *,
long *,
double *,
double *,
4057 long *,
double *,
double *,
double *,
double *,
double *);
4058 extern double d1mach(
long *);
4059 extern long i1mach(
long *);
4060 extern double dgamln(
double *,
long *);
4061 extern int azsqrt(
double *,
double *,
4062 double *,
double *);
4082 caz =
azabs(zr, zi);
4091 bry[0] =
d1mach(&c__1) * 1e3 / *tol;
4092 bry[1] = 1. / bry[0];
4099 sti = -(*zi) * rcaz;
4100 rzr = (str + str) * rcaz;
4101 rzi = (sti + sti) * rcaz;
4102 inu = (long) ((
float) (*fnu + .5));
4103 dnu = *fnu - (double) ((
float) inu);
4104 if (abs(dnu) == .5) {
4108 if (abs(dnu) > *tol) {
4118 azlog(&rzr, &rzi, &smur, &smui, &idum);
4121 zshch(&fmur, &fmui, &cshr, &cshi, &cchr, &cchi);
4134 t2 = exp(-
dgamln(&a2, &idum));
4135 t1 = 1. / (t2 * fc);
4136 if (abs(dnu) > .1) {
4144 for (k = 2; k <= 8; ++k) {
4146 tm = cc[k - 1] * ak;
4148 if (abs(tm) < *tol) {
4157 g1 = (t1 - t2) / (dnu + dnu);
4159 g2 = (t1 + t2) * .5;
4160 fr = fc * (cchr * g1 + smur * g2);
4161 fi = fc * (cchi * g1 + smui * g2);
4162 azexp(&fmur, &fmui, &str, &sti);
4165 zdiv(&c_b219, &c_b220, &str, &sti, &ptr, &pti);
4177 if (inu > 0 || *
n > 1) {
4186 zmlt(zr, zi, zr, zi, &czr, &czi);
4189 t1 = caz * .25 * caz;
4191 fr = (fr * ak + pr + qr) / bk;
4192 fi = (fi * ak + pi + qi) / bk;
4193 str = 1. / (ak - dnu);
4196 str = 1. / (ak + dnu);
4199 str = ckr * czr - cki * czi;
4201 cki = (ckr * czi + cki * czr) * rak;
4203 s1r = ckr * fr - cki * fi + s1r;
4204 s1i = ckr * fi + cki * fr + s1i;
4206 bk = bk + ak + ak + 1.;
4217 azexp(zr, zi, &str, &sti);
4218 zmlt(&s1r, &s1i, &str, &sti, &yr[1], &yi[1]);
4227 zmlt(zr, zi, zr, zi, &czr, &czi);
4230 t1 = caz * .25 * caz;
4232 fr = (fr * ak + pr + qr) / bk;
4233 fi = (fi * ak + pi + qi) / bk;
4234 str = 1. / (ak - dnu);
4237 str = 1. / (ak + dnu);
4240 str = ckr * czr - cki * czi;
4242 cki = (ckr * czi + cki * czr) * rak;
4244 s1r = ckr * fr - cki * fi + s1r;
4245 s1i = ckr * fi + cki * fr + s1i;
4248 s2r = ckr * str - cki * sti + s2r;
4249 s2i = ckr * sti + cki * str + s2i;
4251 bk = bk + ak + ak + 1.;
4259 ak = a1 * abs(smur);
4263 str = cssr[kflag - 1];
4266 zmlt(&p2r, &p2i, &rzr, &rzi, &s2r, &s2i);
4272 azexp(zr, zi, &fr, &fi);
4273 zmlt(&s1r, &s1i, &fr, &fi, &s1r, &s1i);
4274 zmlt(&s2r, &s2i, &fr, &fi, &s2r, &s2i);
4283 azsqrt(zr, zi, &str, &sti);
4284 zdiv(&rthpi, &czeroi, &str, &sti, &coefr, &coefi);
4293 str = exp(-(*zr)) * cssr[kflag - 1];
4294 sti = -str * sin(*zi);
4296 zmlt(&coefr, &coefi, &str, &sti, &coefr, &coefi);
4298 if (abs(dnu) == .5) {
4304 ak = cos(dpi * dnu);
4309 fhs = (d__1 = .25 - dnu2, abs(d__1));
4310 if (fhs == czeror) {
4319 t1 = (double) ((
float) (
i1mach(&c__14) - 1));
4320 t1 = t1 *
d1mach(&c__5) * 3.321928094;
4330 t1 = atan(*zi / *zr);
4339 etest = ak / (dpi * caz * *tol);
4341 if (etest < coner) {
4345 ckr = caz + caz + ctwor;
4349 for (i__ = 1; i__ <= i__1; ++i__) {
4351 cbr = ckr / (fk + coner);
4353 p2r = cbr * p2r - p1r * ak;
4356 fks = fks + fk + fk + ctwor;
4357 fhs = fhs + fk + fk;
4359 str = abs(p2r) * fk;
4367 fk += spi * t1 * sqrt(t2 / caz);
4368 fhs = (d__1 = .25 - dnu2, abs(d__1));
4375 ak = fpi * ak / (*tol * sqrt(a2));
4376 aa = t1 * 3. / (caz + 1.);
4377 bb = t1 * 14.7 / (caz + 28.);
4378 ak = (log(ak) + caz * cos(aa) / (caz * .008 + 1.)) / cos(bb);
4379 fk = ak * .12125 * ak / caz + 1.5;
4384 k = (long) ((
float) fk);
4385 fk = (double) ((
float) k);
4394 for (i__ = 1; i__ <= i__1; ++i__) {
4396 ak = (fks + fk) / (a1 + fhs);
4397 rak = 2. / (fk + coner);
4398 cbr = (fk + *zr) * rak;
4402 p2r = (ptr * cbr - pti * cbi - p1r) * ak;
4403 p2i = (pti * cbr + ptr * cbi - p1i) * ak;
4408 fks = a1 - fk + coner;
4416 tm =
azabs(&csr, &csi);
4422 zmlt(&coefr, &coefi, &s1r, &s1i, &str, &sti);
4423 zmlt(&str, &sti, &csr, &csi, &s1r, &s1i);
4424 if (inu > 0 || *
n > 1) {
4437 tm =
azabs(&p2r, &p2i);
4443 zmlt(&p1r, &p1i, &p2r, &p2i, &ptr, &pti);
4444 str = dnu + .5 - ptr;
4446 zdiv(&str, &sti, zr, zi, &str, &sti);
4448 zmlt(&str, &sti, &s1r, &s1i, &s2r, &s2i);
4481 p1r = csrr[kflag - 1];
4482 ascle = bry[kflag - 1];
4484 for (i__ = inub; i__ <= i__1; ++i__) {
4487 s2r = ckr * str - cki * sti + s1r;
4488 s2i = ckr * sti + cki * str + s1i;
4500 p2m =
max(str, sti);
4505 ascle = bry[kflag - 1];
4510 str = cssr[kflag - 1];
4515 p1r = csrr[kflag - 1];
4525 str = csrr[kflag - 1];
4542 p1r = csrr[kflag - 1];
4543 ascle = bry[kflag - 1];
4545 for (i__ = kk; i__ <= i__1; ++i__) {
4548 s2r = ckr * p2r - cki * p2i + s1r;
4549 s2i = cki * p2r + ckr * p2i + s1i;
4563 p2m =
max(str, sti);
4568 ascle = bry[kflag - 1];
4573 str = cssr[kflag - 1];
4578 p1r = csrr[kflag - 1];
4588 elm = exp(-(*elim));
4596 for (i__ = 1; i__ <= i__1; ++i__) {
4599 s2r = str * ckr - sti * cki + s1r;
4600 s2i = sti * ckr + str * cki + s1i;
4605 as =
azabs(&s2r, &s2i);
4608 if (p2r < -(*elim)) {
4611 azlog(&s2r, &s2i, &str, &sti, &idum);
4614 p2m = exp(p2r) / *tol;
4615 p1r = p2m * cos(p2i);
4616 p1i = p2m * sin(p2i);
4617 zuchk(&p1r, &p1i, &nw, &ascle, tol);
4624 if (ic == i__ - 1) {
4674 zkscl(&zdr, &zdi, fnu,
n, &yr[1], &yi[1], nz, &rzr, &rzi, &ascle, tol,
4683 yr[kk] = s1r * csrr[0];
4684 yi[kk] = s1i * csrr[0];
4691 yr[kk] = s2r * csrr[0];
4692 yi[kk] = s2i * csrr[0];
4696 t2 = *fnu + (double) ((
float) (kk - 1));
6300 int zunhj(
double *zr,
double *zi,
double *fnu,
6301 long *ipmtr,
double *tol,
double *phir,
6302 double *phii,
double *argr,
double *argi,
6303 double *zeta1r,
double *zeta1i,
double *zeta2r,
6304 double *zeta2i,
double *asumr,
double *asumi,
6305 double *bsumr,
double *bsumi)
6309 double ar[14] = { 1., .104166666666666667, .0835503472222222222,
6310 .12822657455632716, .291849026464140464, .881627267443757652,
6311 3.32140828186276754, 14.9957629868625547, 78.9230130115865181,
6312 474.451538868264323, 3207.49009089066193, 24086.5496408740049,
6313 198923.119169509794, 1791902.00777534383
6315 double br[14] = { 1., -.145833333333333333,
6316 -.0987413194444444444, -.143312053915895062, -.317227202678413548,
6317 -.942429147957120249, -3.51120304082635426, -15.7272636203680451,
6318 -82.2814390971859444, -492.355370523670524, -3316.21856854797251,
6319 -24827.6742452085896, -204526.587315129788, -1838444.9170682099
6321 double c__[105] = { 1., -.208333333333333333, .125,
6322 .334201388888888889, -.401041666666666667, .0703125,
6323 -1.02581259645061728, 1.84646267361111111, -.8912109375,
6325 4.66958442342624743, -11.2070026162229938, 8.78912353515625,
6326 -2.3640869140625, .112152099609375, -28.2120725582002449,
6327 84.6362176746007346, -91.8182415432400174, 42.5349987453884549,
6328 -7.3687943594796317, .227108001708984375, 212.570130039217123,
6329 -765.252468141181642, 1059.99045252799988, -699.579627376132541,
6330 218.19051174421159, -26.4914304869515555, .572501420974731445,
6331 -1919.457662318407, 8061.72218173730938, -13586.5500064341374,
6332 11655.3933368645332, -5305.64697861340311, 1200.90291321635246,
6333 -108.090919788394656, 1.7277275025844574, 20204.2913309661486,
6334 -96980.5983886375135, 192547.001232531532, -203400.177280415534,
6335 122200.46498301746, -41192.6549688975513, 7109.51430248936372,
6336 -493.915304773088012, 6.07404200127348304, -242919.187900551333,
6337 1311763.6146629772, -2998015.91853810675, 3763271.297656404,
6338 -2813563.22658653411, 1268365.27332162478, -331645.172484563578,
6339 45218.7689813627263, -2499.83048181120962, 24.3805296995560639,
6340 3284469.85307203782, -19706819.1184322269, 50952602.4926646422,
6341 -74105148.2115326577, 66344512.2747290267, -37567176.6607633513,
6342 13288767.1664218183, -2785618.12808645469, 308186.404612662398,
6343 -13886.0897537170405, 110.017140269246738, -49329253.664509962,
6344 325573074.185765749, -939462359.681578403, 1553596899.57058006,
6345 -1621080552.10833708, 1106842816.82301447, -495889784.275030309,
6346 142062907.797533095, -24474062.7257387285, 2243768.17792244943,
6347 -84005.4336030240853, 551.335896122020586, 814789096.118312115,
6348 -5866481492.05184723, 18688207509.2958249, -34632043388.1587779,
6349 41280185579.753974, -33026599749.8007231, 17954213731.1556001,
6350 -6563293792.61928433, 1559279864.87925751, -225105661.889415278,
6351 17395107.5539781645, -549842.327572288687, 3038.09051092238427,
6352 -14679261247.6956167, 114498237732.02581, -399096175224.466498,
6353 819218669548.577329, -1098375156081.22331, 1008158106865.38209,
6354 -645364869245.376503, 287900649906.150589, -87867072178.0232657,
6355 17634730606.8349694, -2167164983.22379509, 143157876.718888981,
6356 -3871833.44257261262, 18257.7554742931747
6358 double alfa[180] = { -.00444444444444444444,
6359 -9.22077922077922078e-4, -8.84892884892884893e-5,
6360 1.65927687832449737e-4, 2.4669137274179291e-4,
6361 2.6599558934625478e-4, 2.61824297061500945e-4,
6362 2.48730437344655609e-4, 2.32721040083232098e-4,
6363 2.16362485712365082e-4, 2.00738858762752355e-4,
6364 1.86267636637545172e-4, 1.73060775917876493e-4,
6365 1.61091705929015752e-4, 1.50274774160908134e-4,
6366 1.40503497391269794e-4, 1.31668816545922806e-4,
6367 1.23667445598253261e-4, 1.16405271474737902e-4,
6368 1.09798298372713369e-4, 1.03772410422992823e-4,
6369 9.82626078369363448e-5, 9.32120517249503256e-5,
6370 8.85710852478711718e-5, 8.42963105715700223e-5,
6371 8.03497548407791151e-5, 7.66981345359207388e-5,
6372 7.33122157481777809e-5, 7.01662625163141333e-5,
6373 6.72375633790160292e-5, 6.93735541354588974e-4,
6374 2.32241745182921654e-4, -1.41986273556691197e-5,
6375 -1.1644493167204864e-4, -1.50803558053048762e-4,
6376 -1.55121924918096223e-4, -1.46809756646465549e-4,
6377 -1.33815503867491367e-4, -1.19744975684254051e-4,
6378 -1.0618431920797402e-4, -9.37699549891194492e-5,
6379 -8.26923045588193274e-5, -7.29374348155221211e-5,
6380 -6.44042357721016283e-5, -5.69611566009369048e-5,
6381 -5.04731044303561628e-5, -4.48134868008882786e-5,
6382 -3.98688727717598864e-5, -3.55400532972042498e-5,
6383 -3.1741425660902248e-5, -2.83996793904174811e-5,
6384 -2.54522720634870566e-5, -2.28459297164724555e-5,
6385 -2.05352753106480604e-5, -1.84816217627666085e-5,
6386 -1.66519330021393806e-5, -1.50179412980119482e-5,
6387 -1.35554031379040526e-5, -1.22434746473858131e-5,
6388 -1.10641884811308169e-5, -3.54211971457743841e-4,
6389 -1.56161263945159416e-4, 3.0446550359493641e-5,
6390 1.30198655773242693e-4, 1.67471106699712269e-4,
6391 1.70222587683592569e-4, 1.56501427608594704e-4,
6392 1.3633917097744512e-4, 1.14886692029825128e-4,
6393 9.45869093034688111e-5, 7.64498419250898258e-5,
6394 6.07570334965197354e-5, 4.74394299290508799e-5,
6395 3.62757512005344297e-5, 2.69939714979224901e-5,
6396 1.93210938247939253e-5, 1.30056674793963203e-5,
6397 7.82620866744496661e-6, 3.59257485819351583e-6,
6398 1.44040049814251817e-7, -2.65396769697939116e-6,
6399 -4.9134686709848591e-6, -6.72739296091248287e-6,
6400 -8.17269379678657923e-6, -9.31304715093561232e-6,
6401 -1.02011418798016441e-5, -1.0880596251059288e-5,
6402 -1.13875481509603555e-5, -1.17519675674556414e-5,
6403 -1.19987364870944141e-5, 3.78194199201772914e-4,
6404 2.02471952761816167e-4, -6.37938506318862408e-5,
6405 -2.38598230603005903e-4, -3.10916256027361568e-4,
6406 -3.13680115247576316e-4, -2.78950273791323387e-4,
6407 -2.28564082619141374e-4, -1.75245280340846749e-4,
6408 -1.25544063060690348e-4, -8.22982872820208365e-5,
6409 -4.62860730588116458e-5, -1.72334302366962267e-5,
6410 5.60690482304602267e-6, 2.313954431482868e-5,
6411 3.62642745856793957e-5, 4.58006124490188752e-5,
6412 5.2459529495911405e-5, 5.68396208545815266e-5,
6413 5.94349820393104052e-5, 6.06478527578421742e-5,
6414 6.08023907788436497e-5, 6.01577894539460388e-5,
6415 5.891996573446985e-5, 5.72515823777593053e-5,
6416 5.52804375585852577e-5, 5.3106377380288017e-5,
6417 5.08069302012325706e-5, 4.84418647620094842e-5,
6418 4.6056858160747537e-5, -6.91141397288294174e-4,
6419 -4.29976633058871912e-4, 1.83067735980039018e-4,
6420 6.60088147542014144e-4, 8.75964969951185931e-4,
6421 8.77335235958235514e-4, 7.49369585378990637e-4,
6422 5.63832329756980918e-4, 3.68059319971443156e-4,
6423 1.88464535514455599e-4, 3.70663057664904149e-5,
6424 -8.28520220232137023e-5, -1.72751952869172998e-4,
6425 -2.36314873605872983e-4, -2.77966150694906658e-4,
6426 -3.02079514155456919e-4, -3.12594712643820127e-4,
6427 -3.12872558758067163e-4, -3.05678038466324377e-4,
6428 -2.93226470614557331e-4, -2.77255655582934777e-4,
6429 -2.59103928467031709e-4, -2.39784014396480342e-4,
6430 -2.20048260045422848e-4, -2.00443911094971498e-4,
6431 -1.81358692210970687e-4, -1.63057674478657464e-4,
6432 -1.45712672175205844e-4, -1.29425421983924587e-4,
6433 -1.14245691942445952e-4, .00192821964248775885,
6434 .00135592576302022234, -7.17858090421302995e-4,
6435 -.00258084802575270346, -.00349271130826168475,
6436 -.00346986299340960628, -.00282285233351310182,
6437 -.00188103076404891354, -8.895317183839476e-4,
6438 3.87912102631035228e-6, 7.28688540119691412e-4,
6439 .00126566373053457758, .00162518158372674427,
6440 .00183203153216373172,
6441 .00191588388990527909, .00190588846755546138,
6442 .00182798982421825727,
6443 .0017038950642112153, .00155097127171097686, .00138261421852276159,
6444 .00120881424230064774, .00103676532638344962,
6445 8.71437918068619115e-4, 7.16080155297701002e-4,
6446 5.72637002558129372e-4, 4.42089819465802277e-4,
6447 3.24724948503090564e-4, 2.20342042730246599e-4,
6448 1.28412898401353882e-4, 4.82005924552095464e-5
6450 double beta[210] = { .0179988721413553309,
6451 .00559964911064388073, .00288501402231132779,
6452 .00180096606761053941,
6453 .00124753110589199202, 9.22878876572938311e-4,
6454 7.14430421727287357e-4, 5.71787281789704872e-4,
6455 4.69431007606481533e-4, 3.93232835462916638e-4,
6456 3.34818889318297664e-4, 2.88952148495751517e-4,
6457 2.52211615549573284e-4, 2.22280580798883327e-4,
6458 1.97541838033062524e-4, 1.76836855019718004e-4,
6459 1.59316899661821081e-4, 1.44347930197333986e-4,
6460 1.31448068119965379e-4, 1.20245444949302884e-4,
6461 1.10449144504599392e-4, 1.01828770740567258e-4,
6462 9.41998224204237509e-5, 8.74130545753834437e-5,
6463 8.13466262162801467e-5, 7.59002269646219339e-5,
6464 7.09906300634153481e-5, 6.65482874842468183e-5,
6465 6.25146958969275078e-5, 5.88403394426251749e-5,
6466 -.00149282953213429172, -8.78204709546389328e-4,
6467 -5.02916549572034614e-4, -2.94822138512746025e-4,
6468 -1.75463996970782828e-4, -1.04008550460816434e-4,
6469 -5.96141953046457895e-5, -3.1203892907609834e-5,
6470 -1.26089735980230047e-5, -2.42892608575730389e-7,
6471 8.05996165414273571e-6, 1.36507009262147391e-5,
6472 1.73964125472926261e-5, 1.9867297884213378e-5,
6473 2.14463263790822639e-5, 2.23954659232456514e-5,
6474 2.28967783814712629e-5, 2.30785389811177817e-5,
6475 2.30321976080909144e-5, 2.28236073720348722e-5,
6476 2.25005881105292418e-5, 2.20981015361991429e-5,
6477 2.16418427448103905e-5, 2.11507649256220843e-5,
6478 2.06388749782170737e-5, 2.01165241997081666e-5,
6479 1.95913450141179244e-5, 1.9068936791043674e-5,
6480 1.85533719641636667e-5, 1.80475722259674218e-5,
6481 5.5221307672129279e-4, 4.47932581552384646e-4,
6482 2.79520653992020589e-4, 1.52468156198446602e-4,
6483 6.93271105657043598e-5, 1.76258683069991397e-5,
6484 -1.35744996343269136e-5, -3.17972413350427135e-5,
6485 -4.18861861696693365e-5, -4.69004889379141029e-5,
6486 -4.87665447413787352e-5, -4.87010031186735069e-5,
6487 -4.74755620890086638e-5, -4.55813058138628452e-5,
6488 -4.33309644511266036e-5, -4.09230193157750364e-5,
6489 -3.84822638603221274e-5, -3.60857167535410501e-5,
6490 -3.37793306123367417e-5, -3.15888560772109621e-5,
6491 -2.95269561750807315e-5, -2.75978914828335759e-5,
6492 -2.58006174666883713e-5, -2.413083567612802e-5,
6493 -2.25823509518346033e-5, -2.11479656768912971e-5,
6494 -1.98200638885294927e-5, -1.85909870801065077e-5,
6495 -1.74532699844210224e-5, -1.63997823854497997e-5,
6496 -4.74617796559959808e-4, -4.77864567147321487e-4,
6497 -3.20390228067037603e-4, -1.61105016119962282e-4,
6498 -4.25778101285435204e-5, 3.44571294294967503e-5,
6499 7.97092684075674924e-5, 1.031382367082722e-4,
6500 1.12466775262204158e-4, 1.13103642108481389e-4,
6501 1.08651634848774268e-4, 1.01437951597661973e-4,
6502 9.29298396593363896e-5, 8.40293133016089978e-5,
6503 7.52727991349134062e-5, 6.69632521975730872e-5,
6504 5.92564547323194704e-5, 5.22169308826975567e-5,
6505 4.58539485165360646e-5, 4.01445513891486808e-5,
6506 3.50481730031328081e-5, 3.05157995034346659e-5,
6507 2.64956119950516039e-5, 2.29363633690998152e-5,
6508 1.97893056664021636e-5, 1.70091984636412623e-5,
6509 1.45547428261524004e-5, 1.23886640995878413e-5,
6510 1.04775876076583236e-5, 8.79179954978479373e-6,
6511 7.36465810572578444e-4, 8.72790805146193976e-4,
6512 6.22614862573135066e-4, 2.85998154194304147e-4,
6513 3.84737672879366102e-6, -1.87906003636971558e-4,
6514 -2.97603646594554535e-4, -3.45998126832656348e-4,
6515 -3.53382470916037712e-4, -3.35715635775048757e-4,
6516 -3.04321124789039809e-4, -2.66722723047612821e-4,
6517 -2.27654214122819527e-4, -1.89922611854562356e-4,
6518 -1.5505891859909387e-4, -1.2377824076187363e-4,
6519 -9.62926147717644187e-5, -7.25178327714425337e-5,
6520 -5.22070028895633801e-5, -3.50347750511900522e-5,
6521 -2.06489761035551757e-5, -8.70106096849767054e-6,
6522 1.1369868667510029e-6, 9.16426474122778849e-6,
6523 1.5647778542887262e-5, 2.08223629482466847e-5,
6524 2.48923381004595156e-5, 2.80340509574146325e-5,
6525 3.03987774629861915e-5, 3.21156731406700616e-5,
6526 -.00180182191963885708, -.00243402962938042533,
6527 -.00183422663549856802, -7.62204596354009765e-4,
6528 2.39079475256927218e-4, 9.49266117176881141e-4,
6529 .00134467449701540359, .00148457495259449178,
6530 .00144732339830617591,
6531 .00130268261285657186, .00110351597375642682,
6532 8.86047440419791759e-4, 6.73073208165665473e-4,
6533 4.77603872856582378e-4, 3.05991926358789362e-4,
6534 1.6031569459472163e-4, 4.00749555270613286e-5,
6535 -5.66607461635251611e-5, -1.32506186772982638e-4,
6536 -1.90296187989614057e-4, -2.32811450376937408e-4,
6537 -2.62628811464668841e-4, -2.82050469867598672e-4,
6538 -2.93081563192861167e-4, -2.97435962176316616e-4,
6539 -2.96557334239348078e-4, -2.91647363312090861e-4,
6540 -2.83696203837734166e-4, -2.73512317095673346e-4,
6541 -2.6175015580676858e-4, .00638585891212050914,
6542 .00962374215806377941, .00761878061207001043,
6543 .00283219055545628054,
6544 -.0020984135201272009, -.00573826764216626498,
6545 -.0077080424449541462, -.00821011692264844401,
6546 -.00765824520346905413, -.00647209729391045177,
6547 -.00499132412004966473, -.0034561228971313328,
6548 -.00201785580014170775, -7.59430686781961401e-4,
6549 2.84173631523859138e-4, .00110891667586337403,
6550 .00172901493872728771, .00216812590802684701,
6551 .00245357710494539735,
6552 .00261281821058334862, .00267141039656276912, .0026520307339598043,
6553 .00257411652877287315, .00245389126236094427,
6554 .00230460058071795494,
6555 .00213684837686712662, .00195896528478870911,
6556 .00177737008679454412,
6557 .00159690280765839059, .00142111975664438546
6559 double gama[30] = { .629960524947436582, .251984209978974633,
6560 .154790300415655846, .110713062416159013, .0857309395527394825,
6561 .0697161316958684292, .0586085671893713576, .0504698873536310685,
6562 .0442600580689154809, .0393720661543509966, .0354283195924455368,
6563 .0321818857502098231, .0294646240791157679, .0271581677112934479,
6564 .0251768272973861779, .0234570755306078891, .0219508390134907203,
6565 .020621082823564624, .0194388240897880846, .0183810633800683158,
6566 .0174293213231963172, .0165685837786612353, .0157865285987918445,
6567 .0150729501494095594, .0144193250839954639, .0138184805735341786,
6568 .0132643378994276568, .0127517121970498651, .0122761545318762767,
6569 .0118338262398482403
6571 double ex1 = .333333333333333333;
6572 double ex2 = .666666666666666667;
6573 double hpi = 1.57079632679489662;
6574 double gpi = 3.14159265358979324;
6575 double thpi = 4.71238898038468986;
6586 double log(
double), pow_dd(
double *,
double *), atan(
double),
6587 cos(
double), sin(
double), sqrt(
double);
6590 long j, k, l, m, l1, l2;
6591 double ac, ap[30], pi[30];
6592 long is, jr, ks, ju;
6593 double pp, wi, pr[30];
6597 double t2i, w2i, t2r, w2r, ang, fn13, fn23;
6599 double cri[14], dri[14];
6601 double zai, zbi, zci, crr[14], drr[14], raw, zar, upi[14], sti,
6602 zbr, zcr, upr[14], str, raw2;
6606 double atol, btol, tfni;
6608 double azth, tzai, tfnr, rfnu;
6609 extern int zdiv(
double *,
double *,
double *
6610 ,
double *,
double *,
double *);
6611 double zthi, test, tzar, zthr, rfnu2;
6612 extern double azabs(
double *,
double *);
6614 extern int azlog(
double *,
double *,
double
6615 *,
double *,
long *);
6616 double ptfni, sumai, sumbi, zetar, ptfnr, razth, sumar, sumbr,
6618 extern double d1mach(
long *);
6619 double rzthr, rtzti, rtztr, przthi;
6620 extern int azsqrt(
double *,
double *,
6621 double *,
double *);
6663 test =
d1mach(&c__1) * 1e3;
6665 if (abs(*zr) > ac || abs(*zi) > ac) {
6668 *zeta1r = (d__1 = log(test), abs(d__1)) * 2. + *fnu;
6680 rfnu2 = rfnu * rfnu;
6684 fn13 = pow_dd(fnu, &ex1);
6687 w2r = coner - zbr * zbr + zbi * zbi;
6688 w2i = conei - zbr * zbi - zbr * zbi;
6689 aw2 =
azabs(&w2r, &w2i);
6705 for (k = 2; k <= 30; ++k) {
6706 pr[k - 1] = pr[k - 2] * w2r - pi[k - 2] * w2i;
6707 pi[k - 1] = pr[k - 2] * w2i + pi[k - 2] * w2r;
6708 sumar += pr[k - 1] * gama[k - 1];
6709 sumai += pi[k - 1] * gama[k - 1];
6710 ap[k - 1] = ap[k - 2] * aw2;
6711 if (ap[k - 1] < *tol) {
6719 zetar = w2r * sumar - w2i * sumai;
6720 zetai = w2r * sumai + w2i * sumar;
6721 *argr = zetar * fn23;
6722 *argi = zetai * fn23;
6723 azsqrt(&sumar, &sumai, &zar, &zai);
6724 azsqrt(&w2r, &w2i, &str, &sti);
6725 *zeta2r = str * *fnu;
6726 *zeta2i = sti * *fnu;
6727 str = coner + ex2 * (zetar * zar - zetai * zai);
6728 sti = conei + ex2 * (zetar * zai + zetai * zar);
6729 *zeta1r = str * *zeta2r - sti * *zeta2i;
6730 *zeta1i = str * *zeta2i + sti * *zeta2r;
6733 azsqrt(&zar, &zai, &str, &sti);
6734 *phir = str * rfn13;
6735 *phii = sti * rfn13;
6745 for (k = 1; k <= i__1; ++k) {
6746 sumbr += pr[k - 1] * beta[k - 1];
6747 sumbi += pi[k - 1] * beta[k - 1];
6756 btol = *tol * (abs(*bsumr) + abs(*bsumi));
6764 for (is = 2; is <= 7; ++is) {
6773 for (k = 1; k <= i__1; ++k) {
6775 sumar += pr[k - 1] * alfa[m - 1];
6776 sumai += pi[k - 1] * alfa[m - 1];
6777 if (ap[k - 1] < atol) {
6783 *asumr += sumar * pp;
6784 *asumi += sumai * pp;
6795 for (k = 1; k <= i__1; ++k) {
6797 sumbr += pr[k - 1] * beta[m - 1];
6798 sumbi += pi[k - 1] * beta[m - 1];
6799 if (ap[k - 1] < atol) {
6805 *bsumr += sumbr * pp;
6806 *bsumi += sumbi * pp;
6811 if (ias == 1 && ibs == 1) {
6829 azsqrt(&w2r, &w2i, &wr, &wi);
6838 zdiv(&str, &sti, &zbr, &zbi, &zar, &zai);
6839 azlog(&zar, &zai, &zcr, &zci, &idum);
6849 zthr = (zcr - wr) * 1.5;
6850 zthi = (zci - wi) * 1.5;
6851 *zeta1r = zcr * *fnu;
6852 *zeta1i = zci * *fnu;
6853 *zeta2r = wr * *fnu;
6854 *zeta2i = wi * *fnu;
6855 azth =
azabs(&zthr, &zthi);
6857 if (zthr >= 0. && zthi < 0.) {
6864 ang = atan(zthi / zthr);
6869 pp = pow_dd(&azth, &ex2);
6871 zetar = pp * cos(ang);
6872 zetai = pp * sin(ang);
6876 *argr = zetar * fn23;
6877 *argi = zetai * fn23;
6878 zdiv(&zthr, &zthi, &zetar, &zetai, &rtztr, &rtzti);
6879 zdiv(&rtztr, &rtzti, &wr, &wi, &zar, &zai);
6882 azsqrt(&tzar, &tzai, &str, &sti);
6883 *phir = str * rfn13;
6884 *phii = sti * rfn13;
6888 raw = 1. / sqrt(aw2);
6891 tfnr = str * rfnu * raw;
6892 tfni = sti * rfnu * raw;
6895 sti = -zthi * razth;
6896 rzthr = str * razth * rfnu;
6897 rzthi = sti * razth * rfnu;
6898 zcr = rzthr * ar[1];
6899 zci = rzthi * ar[1];
6905 str = t2r * c__[1] + c__[2];
6907 upr[1] = str * tfnr - sti * tfni;
6908 upi[1] = str * tfni + sti * tfnr;
6909 *bsumr = upr[1] + zcr;
6910 *bsumi = upi[1] + zci;
6923 btol = *tol * (abs(*bsumr) + abs(*bsumi));
6929 for (lr = 2; lr <= 12; lr += 2) {
6936 for (k = lr; k <= i__1; ++k) {
6943 for (j = 2; j <= i__2; ++j) {
6945 str = zar * t2r - t2i * zai + c__[l - 1];
6946 zai = zar * t2i + zai * t2r;
6950 str = ptfnr * tfnr - ptfni * tfni;
6951 ptfni = ptfnr * tfni + ptfni * tfnr;
6953 upr[kp1 - 1] = ptfnr * zar - ptfni * zai;
6954 upi[kp1 - 1] = ptfni * zar + ptfnr * zai;
6955 crr[ks - 1] = przthr * br[ks];
6956 cri[ks - 1] = przthi * br[ks];
6957 str = przthr * rzthr - przthi * rzthi;
6958 przthi = przthr * rzthi + przthi * rzthr;
6960 drr[ks - 1] = przthr * ar[ks + 1];
6961 dri[ks - 1] = przthi * ar[ks + 1];
6968 sumar = upr[lrp1 - 1];
6969 sumai = upi[lrp1 - 1];
6972 for (jr = 1; jr <= i__1; ++jr) {
6975 sumar + crr[jr - 1] * upr[ju - 1] - cri[jr - 1] * upi[ju -
6978 sumai + crr[jr - 1] * upi[ju - 1] + cri[jr - 1] * upr[ju -
6984 test = abs(sumar) + abs(sumai);
6985 if (pp < *tol && test < *tol) {
6992 sumbr = upr[lr + 1] + upr[lrp1 - 1] * zcr - upi[lrp1 - 1] * zci;
6993 sumbi = upi[lr + 1] + upr[lrp1 - 1] * zci + upi[lrp1 - 1] * zcr;
6996 for (jr = 1; jr <= i__1; ++jr) {
6999 sumbr + drr[jr - 1] * upr[ju - 1] - dri[jr - 1] * upi[ju -
7002 sumbi + drr[jr - 1] * upi[ju - 1] + dri[jr - 1] * upr[ju -
7008 test = abs(sumbr) + abs(sumbi);
7009 if (pp < btol && test < btol) {
7013 if (ias == 1 && ibs == 1) {
7020 str = -(*bsumr) * rfn13;
7021 sti = -(*bsumi) * rfn13;
7022 zdiv(&str, &sti, &rtztr, &rtzti, bsumr, bsumi);
7696 int zunik(
double *zrr,
double *zri,
double *fnu,
7697 long *ikflg,
long *ipmtr,
double *tol,
7698 long *init,
double *phir,
double *phii,
7699 double *zeta1r,
double *zeta1i,
double *zeta2r,
7700 double *zeta2i,
double *sumr,
double *sumi,
7701 double *cwrkr,
double *cwrki)
7709 double con[2] = { .398942280401432678, 1.25331413731550025 };
7710 double c__[120] = { 1., -.208333333333333333, .125,
7711 .334201388888888889, -.401041666666666667, .0703125,
7712 -1.02581259645061728, 1.84646267361111111, -.8912109375,
7714 4.66958442342624743, -11.2070026162229938, 8.78912353515625,
7715 -2.3640869140625, .112152099609375, -28.2120725582002449,
7716 84.6362176746007346, -91.8182415432400174, 42.5349987453884549,
7717 -7.3687943594796317, .227108001708984375, 212.570130039217123,
7718 -765.252468141181642, 1059.99045252799988, -699.579627376132541,
7719 218.19051174421159, -26.4914304869515555, .572501420974731445,
7720 -1919.457662318407, 8061.72218173730938, -13586.5500064341374,
7721 11655.3933368645332, -5305.64697861340311, 1200.90291321635246,
7722 -108.090919788394656, 1.7277275025844574, 20204.2913309661486,
7723 -96980.5983886375135, 192547.001232531532, -203400.177280415534,
7724 122200.46498301746, -41192.6549688975513, 7109.51430248936372,
7725 -493.915304773088012, 6.07404200127348304, -242919.187900551333,
7726 1311763.6146629772, -2998015.91853810675, 3763271.297656404,
7727 -2813563.22658653411, 1268365.27332162478, -331645.172484563578,
7728 45218.7689813627263, -2499.83048181120962, 24.3805296995560639,
7729 3284469.85307203782, -19706819.1184322269, 50952602.4926646422,
7730 -74105148.2115326577, 66344512.2747290267, -37567176.6607633513,
7731 13288767.1664218183, -2785618.12808645469, 308186.404612662398,
7732 -13886.0897537170405, 110.017140269246738, -49329253.664509962,
7733 325573074.185765749, -939462359.681578403, 1553596899.57058006,
7734 -1621080552.10833708, 1106842816.82301447, -495889784.275030309,
7735 142062907.797533095, -24474062.7257387285, 2243768.17792244943,
7736 -84005.4336030240853, 551.335896122020586, 814789096.118312115,
7737 -5866481492.05184723, 18688207509.2958249, -34632043388.1587779,
7738 41280185579.753974, -33026599749.8007231, 17954213731.1556001,
7739 -6563293792.61928433, 1559279864.87925751, -225105661.889415278,
7740 17395107.5539781645, -549842.327572288687, 3038.09051092238427,
7741 -14679261247.6956167, 114498237732.02581, -399096175224.466498,
7742 819218669548.577329, -1098375156081.22331, 1008158106865.38209,
7743 -645364869245.376503, 287900649906.150589, -87867072178.0232657,
7744 17634730606.8349694, -2167164983.22379509, 143157876.718888981,
7745 -3871833.44257261262, 18257.7554742931747, 286464035717.679043,
7746 -2406297900028.50396, 9109341185239.89896, -20516899410934.4374,
7747 30565125519935.3206, -31667088584785.1584, 23348364044581.8409,
7748 -12320491305598.2872, 4612725780849.13197, -1196552880196.1816,
7749 205914503232.410016, -21822927757.5292237, 1247009293.51271032,
7750 -29188388.1222208134, 118838.426256783253
7762 double ac, si, ti, sr, tr, t2i, t2r, rfn, sri, sti, zni, srr,
7765 extern int zdiv(
double *,
double *,
double *
7766 ,
double *,
double *,
double *);
7767 double test, crfni, crfnr;
7768 extern int azlog(
double *,
double *,
double
7769 *,
double *,
long *);
7770 extern double d1mach(
long *);
7771 extern int azsqrt(
double *,
double *,
7772 double *,
double *);
7812 test =
d1mach(&c__1) * 1e3;
7814 if (abs(*zrr) > ac || abs(*zri) > ac) {
7817 *zeta1r = (d__1 = log(test), abs(d__1)) * 2. + *fnu;
7827 sr = coner + (tr * tr - ti * ti);
7828 si = conei + (tr * ti + ti * tr);
7829 azsqrt(&sr, &si, &srr, &sri);
7832 zdiv(&str, &sti, &tr, &ti, &znr, &zni);
7833 azlog(&znr, &zni, &str, &sti, &idum);
7834 *zeta1r = *fnu * str;
7835 *zeta1i = *fnu * sti;
7836 *zeta2r = *fnu * srr;
7837 *zeta2i = *fnu * sri;
7838 zdiv(&coner, &conei, &srr, &sri, &tr, &ti);
7841 azsqrt(&srr, &sri, &cwrkr[16], &cwrki[16]);
7842 *phir = cwrkr[16] * con[*ikflg - 1];
7843 *phii = cwrki[16] * con[*ikflg - 1];
7847 zdiv(&coner, &conei, &sr, &si, &t2r, &t2i);
7854 for (k = 2; k <= 15; ++k) {
7858 for (j = 1; j <= i__1; ++j) {
7860 str = sr * t2r - si * t2i + c__[l - 1];
7861 si = sr * t2i + si * t2r;
7865 str = crfnr * srr - crfni * sri;
7866 crfni = crfnr * sri + crfni * srr;
7868 cwrkr[k] = crfnr * sr - crfni * si;
7869 cwrki[k] = crfnr * si + crfni * sr;
7871 test = (d__1 = cwrkr[k], abs(d__1)) + (d__2 = cwrki[k], abs(d__2));
7872 if (ac < *tol && test < *tol) {
7890 for (i__ = 1; i__ <= i__1; ++i__) {
7897 *phir = cwrkr[16] * con[0];
7898 *phii = cwrki[16] * con[0];
7908 for (i__ = 1; i__ <= i__1; ++i__) {
7909 sr += tr * cwrkr[i__];
7910 si += tr * cwrki[i__];
7916 *phir = cwrkr[16] * con[1];
7917 *phii = cwrki[16] * con[1];
8511 int zunk2(
double *zr,
double *zi,
double *fnu,
8512 long *kode,
long *mr,
long *
n,
double *yr,
8513 double *yi,
long *nz,
double *tol,
double *elim,
8522 double cr1i = 1.73205080756887729;
8524 double cr2i = -.866025403784438647;
8525 double hpi = 1.57079632679489662;
8526 double pi = 3.14159265358979324;
8527 double aic = 1.26551212348464539;
8528 double cipr[4] = { 1., 0., -1., 0. };
8529 double cipi[4] = { 0., -1., 0., 1. };
8535 double cos(
double), sin(
double), log(
double), exp(
double),
8536 d_sign(
double *,
double *);
8539 long i__, j, k, ib, ic;
8541 long il, kk, in, nw;
8542 double yy, c1i, c2i, c2m, c1r, c2r, s1i, s2i, rs1, s1r, s2r,
8543 aii, ang, asc, car, cki, fnf;
8549 double cyi[2], fmr, sar, csr, sgn, zbi;
8551 double bry[3], cyr[2], pti, sti, zbr, zni, rzi, ptr, zri, str,
8552 znr, rzr, zrr, daii, aarg;
8554 double dair, aphi, argi[2], cscl, phii[2], crsc, argr[2];
8556 double phir[2], csrr[3], cssr[3], rast, razr;
8557 extern int zs1s2(
double *,
double *,
double
8558 *,
double *,
double *,
double *,
8559 long *,
double *,
double *,
long *);
8561 double argdi, ascle;
8563 double phidi, argdr;
8564 extern double azabs(
double *,
double *);
8566 double csgni, phidr, cspni, asumi[2], bsumi[2];
8567 extern int zuchk(
double *,
double *,
long *,
8568 double *,
double *);
8569 double cspnr, asumr[2], bsumr[2];
8570 extern double d1mach(
long *);
8571 extern int zunhj(
double *,
double *,
double
8572 *,
long *,
double *,
double *,
8573 double *,
double *,
double *,
8574 double *,
double *,
double *,
8575 double *,
double *,
double *,
8576 double *,
double *),
zairy(
double *,
8584 double zeta1i[2], zeta2i[2], zet1di, zet2di, zeta1r[2],
8585 zeta2r[2], zet1dr, zet2dr, asumdi, bsumdi, asumdr, bsumdr;
8624 bry[0] =
d1mach(&c__1) * 1e3 / *tol;
8625 bry[1] = 1. / bry[0];
8640 inu = (long) ((
float) (*fnu));
8641 fnf = *fnu - (double) ((
float) inu);
8648 str = c2r * cipr[kk - 1] - c2i * cipi[kk - 1];
8649 sti = c2r * cipi[kk - 1] + c2i * cipr[kk - 1];
8650 csr = cr1r * str - cr1i * sti;
8651 csi = cr1r * sti + cr1i * str;
8665 for (i__ = 1; i__ <= i__1; ++i__) {
8670 fn = *fnu + (double) ((
float) (i__ - 1));
8671 zunhj(&znr, &zni, &fn, &c__0, tol, &phir[j - 1], &phii[j - 1],
8672 &argr[j - 1], &argi[j - 1], &zeta1r[j - 1], &zeta1i[j - 1],
8673 &zeta2r[j - 1], &zeta2i[j - 1], &asumr[j - 1], &asumi[j - 1],
8674 &bsumr[j - 1], &bsumi[j - 1]);
8678 str = zbr + zeta2r[j - 1];
8679 sti = zbi + zeta2i[j - 1];
8680 rast = fn /
azabs(&str, &sti);
8681 str = str * rast * rast;
8682 sti = -sti * rast * rast;
8683 s1r = zeta1r[j - 1] - str;
8684 s1i = zeta1i[j - 1] - sti;
8687 s1r = zeta1r[j - 1] - zeta2r[j - 1];
8688 s1i = zeta1i[j - 1] - zeta2i[j - 1];
8694 if (abs(rs1) > *elim) {
8700 if (abs(rs1) < *alim) {
8706 aphi =
azabs(&phir[j - 1], &phii[j - 1]);
8707 aarg =
azabs(&argr[j - 1], &argi[j - 1]);
8708 rs1 = rs1 + log(aphi) - log(aarg) * .25 - aic;
8709 if (abs(rs1) > *elim) {
8726 c2r = argr[j - 1] * cr2r - argi[j - 1] * cr2i;
8727 c2i = argr[j - 1] * cr2i + argi[j - 1] * cr2r;
8728 zairy(&c2r, &c2i, &c__0, &c__2, &air, &aii, &nai, &idum);
8729 zairy(&c2r, &c2i, &c__1, &c__2, &dair, &daii, &ndai, &idum);
8730 str = dair * bsumr[j - 1] - daii * bsumi[j - 1];
8731 sti = dair * bsumi[j - 1] + daii * bsumr[j - 1];
8732 ptr = str * cr2r - sti * cr2i;
8733 pti = str * cr2i + sti * cr2r;
8734 str = ptr + (air * asumr[j - 1] - aii * asumi[j - 1]);
8735 sti = pti + (air * asumi[j - 1] + aii * asumr[j - 1]);
8736 ptr = str * phir[j - 1] - sti * phii[j - 1];
8737 pti = str * phii[j - 1] + sti * phir[j - 1];
8738 s2r = ptr * csr - pti * csi;
8739 s2i = ptr * csi + pti * csr;
8740 str = exp(s1r) * cssr[kflag - 1];
8741 s1r = str * cos(s1i);
8742 s1i = str * sin(s1i);
8743 str = s2r * s1r - s2i * s1i;
8744 s2i = s1r * s2i + s2r * s1i;
8749 zuchk(&s2r, &s2i, &nw, bry, tol);
8757 cyr[kdflg - 1] = s2r;
8758 cyi[kdflg - 1] = s2i;
8759 yr[i__] = s2r * csrr[kflag - 1];
8760 yi[i__] = s2i * csrr[kflag - 1];
8789 if (yr[i__ - 1] == zeror && yi[i__ - 1] == zeroi) {
8792 yr[i__ - 1] = zeror;
8793 yi[i__ - 1] = zeroi;
8800 razr = 1. /
azabs(&zrr, &zri);
8803 rzr = (str + str) * razr;
8804 rzi = (sti + sti) * razr;
8815 fn = *fnu + (double) ((
float) (*
n - 1));
8820 zunhj(&znr, &zni, &fn, &ipard, tol, &phidr, &phidi, &argdr, &argdi,
8821 &zet1dr, &zet1di, &zet2dr, &zet2di, &asumdr, &asumdi, &bsumdr,
8828 rast = fn /
azabs(&str, &sti);
8829 str = str * rast * rast;
8830 sti = -sti * rast * rast;
8835 s1r = zet1dr - zet2dr;
8836 s1i = zet1di - zet2di;
8839 if (abs(rs1) > *elim) {
8842 if (abs(rs1) < *alim) {
8848 aphi =
azabs(&phidr, &phidi);
8850 if (abs(rs1) < *elim) {
8865 for (i__ = 1; i__ <= i__1; ++i__) {
8876 c1r = csrr[kflag - 1];
8877 ascle = bry[kflag - 1];
8879 for (i__ = ib; i__ <= i__1; ++i__) {
8882 s2r = ckr * c2r - cki * c2i + s1r;
8883 s2i = ckr * c2i + cki * c2r + s1i;
8897 c2m =
max(str, sti);
8902 ascle = bry[kflag - 1];
8907 s1r *= cssr[kflag - 1];
8908 s1i *= cssr[kflag - 1];
8909 s2r *= cssr[kflag - 1];
8910 s2i *= cssr[kflag - 1];
8911 c1r = csrr[kflag - 1];
8923 fmr = (double) ((
float) (*mr));
8924 sgn = -
d_sign(&pi, &fmr);
8953 str = csr * c2r + csi * c2i;
8954 csi = -csr * c2i + csi * c2r;
8963 for (k = 1; k <= i__1; ++k) {
8964 fn = *fnu + (double) ((
float) (kk - 1));
8973 phidr = phir[j - 1];
8974 phidi = phii[j - 1];
8975 argdr = argr[j - 1];
8976 argdi = argi[j - 1];
8977 zet1dr = zeta1r[j - 1];
8978 zet1di = zeta1i[j - 1];
8979 zet2dr = zeta2r[j - 1];
8980 zet2di = zeta2i[j - 1];
8981 asumdr = asumr[j - 1];
8982 asumdi = asumi[j - 1];
8983 bsumdr = bsumr[j - 1];
8984 bsumdi = bsumi[j - 1];
8988 if (kk == *
n && ib < *
n) {
8991 if (kk == ib || kk == ic) {
8994 zunhj(&znr, &zni, &fn, &c__0, tol, &phidr, &phidi, &argdr, &argdi,
8995 &zet1dr, &zet1di, &zet2dr, &zet2di, &asumdr, &asumdi,
9003 rast = fn /
azabs(&str, &sti);
9004 str = str * rast * rast;
9005 sti = -sti * rast * rast;
9006 s1r = -zet1dr + str;
9007 s1i = -zet1di + sti;
9010 s1r = -zet1dr + zet2dr;
9011 s1i = -zet1di + zet2di;
9017 if (abs(rs1) > *elim) {
9023 if (abs(rs1) < *alim) {
9029 aphi =
azabs(&phidr, &phidi);
9030 aarg =
azabs(&argdr, &argdi);
9031 rs1 = rs1 + log(aphi) - log(aarg) * .25 - aic;
9032 if (abs(rs1) > *elim) {
9045 zairy(&argdr, &argdi, &c__0, &c__2, &air, &aii, &nai, &idum);
9046 zairy(&argdr, &argdi, &c__1, &c__2, &dair, &daii, &ndai, &idum);
9047 str = dair * bsumdr - daii * bsumdi;
9048 sti = dair * bsumdi + daii * bsumdr;
9049 str += air * asumdr - aii * asumdi;
9050 sti += air * asumdi + aii * asumdr;
9051 ptr = str * phidr - sti * phidi;
9052 pti = str * phidi + sti * phidr;
9053 s2r = ptr * csr - pti * csi;
9054 s2i = ptr * csi + pti * csr;
9055 str = exp(s1r) * cssr[iflag - 1];
9056 s1r = str * cos(s1i);
9057 s1i = str * sin(s1i);
9058 str = s2r * s1r - s2i * s1i;
9059 s2i = s2r * s1i + s2i * s1r;
9064 zuchk(&s2r, &s2i, &nw, bry, tol);
9074 cyr[kdflg - 1] = s2r;
9075 cyi[kdflg - 1] = s2i;
9078 s2r *= csrr[iflag - 1];
9079 s2i *= csrr[iflag - 1];
9088 zs1s2(&zrr, &zri, &s1r, &s1i, &s2r, &s2i, &nw, &asc, alim, &iuf);
9091 yr[kk] = s1r * cspnr - s1i * cspni + s2r;
9092 yi[kk] = s1r * cspni + s1i * cspnr + s2i;
9099 if (c2r != 0. || c2i != 0.) {
9135 csr = csrr[iflag - 1];
9136 ascle = bry[iflag - 1];
9137 fn = (double) ((
float) (inu + il));
9139 for (i__ = 1; i__ <= i__1; ++i__) {
9142 s2r = s1r + (fn + fnf) * (rzr * c2r - rzi * c2i);
9143 s2i = s1i + (fn + fnf) * (rzr * c2i + rzi * c2r);
9156 zs1s2(&zrr, &zri, &c1r, &c1i, &c2r, &c2i, &nw, &asc, alim, &iuf);
9159 yr[kk] = c1r * cspnr - c1i * cspni + c2r;
9160 yi[kk] = c1r * cspni + c1i * cspnr + c2i;
9169 c2m =
max(c2r, c2i);
9174 ascle = bry[iflag - 1];
9179 s1r *= cssr[iflag - 1];
9180 s1i *= cssr[iflag - 1];
9181 s2r *= cssr[iflag - 1];
9182 s2i *= cssr[iflag - 1];
9183 csr = csrr[iflag - 1];