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);
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];
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;
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);
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,
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];
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];
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);
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);
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];