bessel_olver.c

多个常用的特殊函数的例子

  -1, 1,
  29
};


/* Chebyshev fit for f(x) = z(x) B_2(z(x)), z(x) = 2(x+1)/5 */
static double B2_lt1_data[40] = {
  0.00073681565841337130021924199490,
  0.00033803599647571227535304316937,
 -0.00008251723219239754024210552679,
 -0.00003390879948656432545900779710,
  0.00001961398056848881816694014889,
 -2.35593745904151401624656805567e-06,
 -1.79055017080406086541563835433e-06,
  1.33129571185610681090725934031e-06,
 -5.38879444715436544130673956170e-07,
  1.49603056041381416881299945557e-07,
 -1.83377228267274327911131293091e-08,
 -1.33191430762944336526965187651e-08,
  1.60642096463700438411396889489e-08,
 -1.28932576330421806740136816643e-08,
  9.6169275086179165484403221944e-09,
 -7.1818502280703532276832887290e-09,
  5.4744009217215145730697754561e-09,
 -4.2680446690508456935030086136e-09,
  3.3941665009266174865683284781e-09,
 -2.7440714072221673882163135170e-09,
  2.2488361522108255229193038962e-09,
 -1.8638240716608748862087923337e-09,
  1.5592350940805373500866440401e-09,
 -1.3145743937732330609242633070e-09,
  1.1153716777215047842790244968e-09,
 -9.5117576805266622854647303110e-10,
  8.1428799553234876296804561100e-10,
 -6.9893770813548773664326279169e-10,
  6.0073113636087448745018831981e-10,
 -5.1627434258513453901420776514e-10,
  4.4290993195074905891788459756e-10,
 -3.7852978599966867611179315200e-10,
  3.2143959338863177145307610452e-10,
 -2.7025926680620777594992221143e-10,
  2.2384857772457918539228234321e-10,
 -1.8125071664276678046551271701e-10,
  1.4164870008713668767293008546e-10,
 -1.0433101857132782485813325981e-10,
  6.8663910168392483929411418190e-11,
 -3.4068313177952244040559740439e-11
};
static cheb_series B2_lt1_cs = {
  B2_lt1_data,
  39,
  -1, 1,
  39
};


/* Chebyshev fit for f(x) = B_3(2(x+1)/5) */
static double B3_lt1_data[40] = {
 -0.00137160820526992057354001614451,
 -0.00025474937951101049982680561302,
  0.00024762975547895881652073467771,
  0.00005229657281480196749313930265,
 -0.00007488354272621512385016593760,
  0.00001416880012891046449980449746,
  0.00001528986060172183690742576230,
 -0.00001668672297078590514293325326,
  0.00001061765189536459018739585094,
 -5.8220577442406209989680801335e-06,
  3.3322423743855900506302033234e-06,
 -2.23292405803003860894449897815e-06,
  1.74816651036678291794777245325e-06,
 -1.49581306041395051804547535093e-06,
  1.32759146107893129050610165582e-06,
 -1.19376077392564467408373553343e-06,
  1.07878303863211630544654040875e-06,
 -9.7743335011819134006676476250e-07,
  8.8729318903693324226127054792e-07,
 -8.0671146292125665050876015280e-07,
  7.3432860378667354971042255937e-07,
 -6.6897926072697370325310483359e-07,
  6.0966619703735610352576581485e-07,
 -5.5554095284507959561958605420e-07,
  5.0588335673197236002812826526e-07,
 -4.6008146297767601862670079590e-07,
  4.1761348515688145911438168306e-07,
 -3.7803230006989446874174476515e-07,
  3.4095248501364300041684648230e-07,
 -3.0603959751354749520615015472e-07,
  2.7300134179365690589640458993e-07,
 -2.4158028250762304756044254231e-07,
  2.1154781038298751985689113868e-07,
 -1.8269911328756771201465223313e-07,
  1.5484895085808513749026173074e-07,
 -1.2782806851555809369226440495e-07,
  1.0148011725394892565174207341e-07,
 -7.5658969771439627809239950461e-08,
  5.0226342286491286957075289622e-08,
 -2.5049645660259882970547555831e-08
};
static cheb_series B3_lt1_cs = {
  B3_lt1_data,
  39,
  -1, 1,
  39
};


/* Chebyshev fit for f(x) = A_3(9(x+1)/20) */
static double A3_lt1_data[40] = {
  -0.00017982561472134418587634980117,
  -0.00036558603837525275836608884064,
  -0.00002819398055929628850294406363,
   0.00016704539863875736769812786067,
  -0.00007098969970347674307623044850,
  -8.4470843942344237748899879940e-06,
   0.0000273413090343147765148014327150,
  -0.0000199073838489821681991178018081,
   0.0000100004176278235088881096950105,
  -3.9739852013143676487867902026e-06,
   1.2265357766449574306882693267e-06,
  -1.88755584306424047416914864854e-07,
  -1.37482206060161206336523452036e-07,
   2.10326379301853336795686477738e-07,
  -2.05583778245412633433934301948e-07,
   1.82377384812654863038691147988e-07,
  -1.58130247846381041027699152436e-07,
   1.36966982725588978654041029615e-07,
  -1.19250280944620257443805710485e-07,
   1.04477169029350256435316644493e-07,
  -9.2064832489437534542041040184e-08,
   8.1523798290458784610230199344e-08,
  -7.2471794980050867512294061891e-08,
   6.4614432955971132569968860233e-08,
  -5.7724095125560946811081322985e-08,
   5.1623107567436835158110947901e-08,
  -4.6171250746798606260216486042e-08,
   4.1256621998650164023254101585e-08,
  -3.6788925543159819135102047082e-08,
   3.2694499457951844422299750661e-08,
  -2.89125899697964696586521743928e-08,
   2.53925288725374047626589488217e-08,
  -2.20915707933726481321465184207e-08,
   1.89732166352720474944407102940e-08,
  -1.60058977893259856012119939554e-08,
   1.31619294542205876946742394494e-08,
  -1.04166651771938038563454275883e-08,
   7.7478015858156185064152078434e-09,
  -5.1347942579352613057675111787e-09,
   2.5583541594586723967261504321e-09
};
static cheb_series A3_lt1_cs = {
  A3_lt1_data,
  39,
  -1, 1,
  39
};

/* chebyshev fit for f(x) = A_4(2(x+1)/5) */
static double A4_lt1_data[30] = {
  0.00009054703770051610946958226736,
  0.00033066000498098017589672988293,
  0.00019737453734363989127226073272,
 -0.00015490809725932037720034762889,
 -0.00004514948935538730085479280454,
  0.00007976881782603940889444573924,
 -0.00003314566154544740986264993251,
 -1.88212148790135672249935711657e-06,
  0.0000114788756505519986352882940648,
 -9.2263039911196207101468331210e-06,
  5.1401128250377780476084336340e-06,
 -2.38418218951722002658891397905e-06,
  1.00664292214481531598338960828e-06,
 -4.23224678096490060264249970540e-07,
  2.00132031535793489976535190025e-07,
 -1.18689501178886741400633921047e-07,
  8.7819524319114212999768013738e-08,
 -7.3964150324206644900787216386e-08,
  6.5780431507637165113885884236e-08,
 -5.9651053193022652369837650411e-08,
  5.4447762662767276209052293773e-08,
 -4.9802057381568863702541294988e-08,
  4.5571368194694340198117635845e-08,
 -4.1682117173547642845382848197e-08,
  3.8084701352766049815367147717e-08,
 -3.4740302885185237434662649907e-08,
  3.1616557064701510611273692060e-08,
 -2.8685739487689556252374879267e-08,
  2.5923752117132254429002796600e-08,
 -2.3309428552190587304662883477e-08
};
static cheb_series A4_lt1_cs = {
  A4_lt1_data,
  29,
  -1, 1,
  29
};


static double olver_B0(double z, double abs_zeta)
{
  if(z < 0.98) {
    const double t = 1.0/sqrt(1.0-z*z);
    return -5.0/(48.0*abs_zeta*abs_zeta) + t*(-3.0 + 5.0*t*t)/(24.0*sqrt(abs_zeta));
  }
  else if(z < 1.02) {
    const double a = 1.0-z;
    const double c0 =  0.0179988721413553309252458658183;
    const double c1 =  0.0111992982212877614645974276203;
    const double c2 =  0.0059404069786014304317781160605;
    const double c3 =  0.0028676724516390040844556450173;
    const double c4 =  0.0012339189052567271708525111185;
    const double c5 =  0.0004169250674535178764734660248;
    const double c6 =  0.0000330173385085949806952777365;
    const double c7 = -0.0001318076238578203009990106425;
    const double c8 = -0.0001906870370050847239813945647;
    return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*(c6 + a*(c7 + a*c8)))))));
  }
  else {
    const double t = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
    return -5.0/(48.0*abs_zeta*abs_zeta) + t*(3.0 + 5.0*t*t)/(24.0*sqrt(abs_zeta));
  }
}


static double olver_B1(double z, double abs_zeta)
{
  if(z < 0.88) {
    const double t   = 1.0/sqrt(1.0-z*z);
    const double t2  = t*t;
    const double rz  = sqrt(abs_zeta);
    const double z32 = rz*rz*rz;
    const double z92 = z32*z32*z32;
    const double term1 = t*t*t * (30375.0 - 369603.0*t2 + 765765.0*t2*t2 - 425425.0*t2*t2*t2)/414720.0;
    const double term2 = 85085.0/(663552.0*z92);
    const double term3 = 385.0/110592.*t*(3.0-5.0*t2)/(abs_zeta*abs_zeta*abs_zeta);
    const double term4 = 5.0/55296.0*t2*(81.0 - 462.0*t2 + 385.0*t2*t2)/z32;
    return -(term1 + term2 + term3 + term4)/rz;
  }
  else if(z < 1.12) {
    const double a = 1.0-z;
    const double c0  = -0.00149282953213429172050073403334;
    const double c1  = -0.00175640941909277865678308358128;
    const double c2  = -0.00113346148874174912576929663517;
    const double c3  = -0.00034691090981382974689396961817;
    const double c4  =  0.00022752516104839243675693256916;
    const double c5  =  0.00051764145724244846447294636552;
    const double c6  =  0.00058906174858194233998714243010;
    const double c7  =  0.00053485514521888073087240392846;
    const double c8  =  0.00042891792986220150647633418796;
    const double c9  =  0.00031639765900613633260381972850;
    const double c10 =  0.00021908147678699592975840749194;
    return c0+a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10)))))))));
  }
  else {
    const double t   = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
    const double t2  = t*t;
    const double rz  = sqrt(abs_zeta);
    const double z32 = rz*rz*rz;
    const double z92 = z32*z32*z32;
    const double term1 = -t2*t * (30375.0 + 369603.0*t2 + 765765.0*t2*t2 + 425425.0*t2*t2*t2)/414720.0;
    const double term2 = 85085.0/(663552.0*z92);
    const double term3 = -385.0/110592.0*t*(3.0+5.0*t2)/(abs_zeta*abs_zeta*abs_zeta);
    const double term4 = 5.0/55296.0*t2*(81.0 + 462.0*t2 + 385.0*t2*t2)/z32;
    return (term1 + term2 + term3 + term4)/rz;
  }
}


static double olver_B2(double z, double abs_zeta)
{
  if(z < 0.8) {
    const double x = 5.0*z/2.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&B2_lt1_cs, x, &c);
    return  c.val / z;
  }
  else if(z <= 1.2) {
    const double a = 1.0-z;
    const double c0 = 0.00055221307672129279005986982501;
    const double c1 = 0.00089586516310476929281129228969;
    const double c2 = 0.00067015003441569770883539158863;
    const double c3 = 0.00010166263361949045682945811828;
    const double c4 = -0.00044086345133806887291336488582;
    const double c5 = -0.00073963081508788743392883072523;
    const double c6 = -0.00076745494377839561259903887331;
    const double c7 = -0.00060829038106040362291568012663;
    const double c8 = -0.00037128707528893496121336168683;
    const double c9 = -0.00014116325105702609866850307176;
    return c0+a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*c9))))))));
  }
  else {
    const double zi = 1.0/z;
    const double x  = 12.0/5.0 * zi - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&B2_gt1_cs, x, &c);
    return c.val * zi*zi*zi;
  }
}


static double olver_B3(double z, double abs_zeta)
{
  if(z < 0.8) {
    const double x = 5.0*z/2.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&B3_lt1_cs, x, &c);
    return c.val;
  }
  else if(z < 1.2) {
    const double a = 1.0-z;
    const double c0 = -0.00047461779655995980754441833105;
    const double c1 = -0.00095572913429464297452176811898;
    const double c2 = -0.00080369634512082892655558133973;
    const double c3 = -0.00000727921669154784138080600339;
    const double c4 =  0.00093162500331581345235746518994;
    const double c5 =  0.00149848796913751497227188612403;
    const double c6 =  0.00148406039675949727870390426462;
    return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*c6)))));
  }
  else {
    const double x   = 12.0/(5.0*z) - 1.0;
    const double zi2 = 1.0/(z*z);
    gsl_sf_result c;
    cheb_eval_e(&B3_gt1_cs, x, &c);
    return  c.val * zi2*zi2*zi2;
  }
}