/*
* fermi.i - $Id$
* Fermi-Dirac integrals and inverses of orders -1/2, 1/2, 3/2, 5/2
*/
local fermi ;
/* DOCUMENT #include "fermi.i"
* Fermi-Dirac integrals and inverses of orders -1/2, 1/2, 3/2, 5/2
*
* Antia, H. M., Aph.J. 84, p.101-108 (1993)
*
* SEE ALSO: fdm12, fd12, fd32, fd52, ifdm12, ifd12, ifd32, ifd52
*/
func fdm12 (x)
/* DOCUMENT fdm12(x)
* return Fermi-Dirac integral of order -1/2,
* fdm12(x) = integral[0 to inf]{ dt * t^-0.5 / (exp(t-x)+1) }
* accurate to about 1e-12
* SEE ALSO: fd12, fd32, fd52, ifdm12, ifd12, ifd32, ifd52
*/
{
mask = x < 2.;
list = where(mask);
if (numberof(list)) {
a = exp(x(list));
a = a*poly(a, 1.71446374704454e7, 3.88148302324068e7,
3.16743385304962e7, 1.14587609192151e7,
1.83696370756153e6, 1.14980998186874e5,
1.98276889924768e3, 1.0) /
poly(a, 9.67282587452899e6, 2.87386436731785e7,
3.26070130734158e7, 1.77657027846367e7,
4.81648022267831e6, 6.13709569333207e5,
3.13595854332114e4, 4.35061725080755e2);
}
list = where(!mask);
if (numberof(list)) {
x = x(list);
b = 1./(x*x);
b = sqrt(x)*poly(b, -4.46620341924942e-15, -1.58654991146236e-12,
-4.44467627042232e-10, -6.84738791621745e-8,
-6.64932238528105e-6, -3.69976170193942e-4,
-1.12295393687006e-2, -1.60926102124442e-1,
-8.52408612877447e-1, -7.45519953763928e-1,
2.98435207466372e0, 1.0) /
poly(b, -2.23310170962369e-15, -7.94193282071464e-13,
-2.22564376956228e-10, -3.43299431079845e-8,
-3.33919612678907e-6, -1.86432212187088e-4,
-5.69764436880529e-3, -8.34904593067194e-2,
-4.78770844009440e-1, -4.99759250374148e-1,
1.86795964993052e0, 4.16485970495288e-1);
}
return merge(a, b, mask);
}
func fd12 (x)
/* DOCUMENT fd12(x)
* return Fermi-Dirac integral of order 1/2,
* fd12(x) = integral[0 to inf]{ dt * t^0.5 / (exp(t-x)+1) }
* accurate to about 1e-12
* SEE ALSO: fdm12, fd32, fd52, ifdm12, ifd12, ifd32, ifd52
*/
{
mask = x < 2.;
list = where(mask);
if (numberof(list)) {
a = exp(x(list));
a = a*poly(a, 5.75834152995465e6, 1.30964880355883e7,
1.07608632249013e7, 3.93536421893014e6,
6.42493233715640e5, 4.16031909245777e4,
7.77238678539648e2, 1.0) /
poly(a, 6.49759261942269e6, 1.70750501625775e7,
1.69288134856160e7, 7.95192647756086e6,
1.83167424554505e6, 1.95155948326832e5,
8.17922106644547e3, 9.02129136642157e1);
}
list = where(!mask);
if (numberof(list)) {
x = x(list);
b = 1./(x*x);
b = x*sqrt(x)*poly(b, 4.85378381173415e-14, 1.64429113030738e-11,
3.76794942277806e-9, 4.69233883900644e-7,
3.40679845803144e-5, 1.32212995937796e-3,
2.60768398973913e-2, 2.48653216266227e-1,
1.08037861921488e0, 1.91247528779676e0, 1.0) /
poly(b, 7.28067571760518e-14, 2.45745452167585e-11,
5.62152894375277e-9, 6.96888634549649e-7,
5.02360015186394e-5, 1.92040136756592e-3,
3.66887808002874e-2, 3.24095226486468e-1,
1.16434871200131e0, 1.34981244060549e0,
2.01311836975930e-1, -2.14562434782759e-2);
}
return merge(a, b, mask);
}
func fd32 (x)
/* DOCUMENT fd32(x)
* return Fermi-Dirac integral of order 3/2,
* fd32(x) = integral[0 to inf]{ dt * t^1.5 / (exp(t-x)+1) }
* accurate to about 1e-12
* SEE ALSO: fdm12, fd12, fd52, ifdm12, ifd12, ifd32, ifd52
*/
{
mask = x < 2.;
list = where(mask);
if (numberof(list)) {
a = exp(x(list));
a = a*poly(a, 4.32326386604283e4, 8.55472308218786e4,
5.95275291210962e4, 1.77294861572005e4,
2.21876607796460e3, 9.90562948053193e1, 1.0) /
poly(a, 3.25218725353467e4, 7.01022511904373e4,
5.50859144223638e4, 1.95942074576400e4,
3.20803912586318e3, 2.20853967067789e2,
5.05580641737527e0, 1.99507945223266e-2);
}
list = where(!mask);
if (numberof(list)) {
x = x(list);
b = 1./(x*x);
b = sqrt(x)*poly(b, 2.80452693148553e-13, 8.60096863656367e-11,
1.62974620742993e-8, 1.63598843752050e-6,
9.12915407846722e-5, 2.62988766922117e-3,
3.85682997219346e-2, 2.78383256609605e-1,
9.02250179334496e-1, 1.0) /
(b*poly(b, 7.01131732871184e-13, 2.10699282897576e-10,
3.94452010378723e-8, 3.84703231868724e-6,
2.04569943213216e-4, 5.31999109566385e-3,
6.39899717779153e-2, 3.14236143831882e-1,
4.70252591891375e-1, -2.15540156936373e-2,
2.34829436438087e-3));
}
return merge(a, b, mask);
}
func fd52 (x)
/* DOCUMENT fd52(x)
* return Fermi-Dirac integral of order 5/2,
* fd52(x) = integral[0 to inf]{ dt * t^2.5 / (exp(t-x)+1) }
* accurate to about 1e-12
* SEE ALSO: fdm12, fd12, fd32, ifdm12, ifd12, ifd32, ifd52
*/
{
mask = x < 2.;
list = where(mask);
if (numberof(list)) {
a = exp(x(list));
a = a*poly(a, 6.61606300631656e4, 1.20132462801652e5,
7.67255995316812e4, 2.10427138842443e4,
2.44325236813275e3, 1.02589947781696e2, 1.0) /
poly(a, 1.99078071053871e4, 3.79076097261066e4,
2.60117136841197e4, 7.97584657659364e3,
1.10886130159658e3, 6.35483623268093e1,
1.16951072617142e0, 3.31482978240026e-3);
}
list = where(!mask);
if (numberof(list)) {
x = x(list);
b = 1./(x*x);
b = x*sqrt(x)*poly(b, 8.42667076131315e-12, 2.31618876821567e-9,
3.54323824923987e-7, 2.77981736000034e-5,
1.14008027400645e-3, 2.32779790773633e-2,
2.39564845938301e-1, 1.24415366126179e0,
3.18831203950106e0, 3.42040216997894e0, 1.0) /
(b*poly(b, 2.94933476646033e-11, 7.68215783076936e-9,
1.12919616415947e-6, 8.09451165406274e-5,
2.81111224925648e-3, 3.99937801931919e-2,
2.27132567866839e-1, 5.31886045222680e-1,
3.70866321410385e-1, 2.27326643192516e-2));
}
return merge(a, b, mask);
}
func ifdm12 (x)
/* DOCUMENT ifdm12(y)
* return x = inverse of Fermi-Dirac integral of order -1/2,
* y = integral[0 to inf]{ dt * t^-0.5 / (exp(t-x)+1) }
* accurate to about 1e-8
* SEE ALSO: ifd12, ifd32, ifd52, fdm12, fd12, fd32, fd52
*/
{
mask = x < 4.;
list = where(mask);
if (numberof(list)) {
a = x(list);
a = log(a*poly(a, -1.570044577033e4, 1.001958278442e4,
-2.805343454951e3, 4.121170498099e2,
-3.174780572961e1, 1.0) /
poly(a, -2.782831558471e4, 2.886114034012e4,
-1.274243093149e4, 3.063252215963e3,
-4.225615045074e2, 3.168918168284e1,
-1.008561571363e0));
}
list = where(!mask);
if (numberof(list)) {
x = x(list);
b = 1./(x*x);
b = poly(b, 2.206779160034e-8, -1.437701234283e-6,
6.103116850636e-5, -1.169411057416e-3,
1.814141021608e-2, -9.588603457639e-2, 1.0) /
(b*poly(b, 8.827116613576e-8, -5.750804196059e-6,
2.429627688357e-4, -4.601959491394e-3,
6.932122275919e-2, -3.217372489776e-1,
3.124344749296e0));
}
return merge(a, b, mask);
}
func ifd12 (x)
/* DOCUMENT ifd12(y)
* return x = inverse of Fermi-Dirac integral of order 1/2,
* y = integral[0 to inf]{ dt * t^0.5 / (exp(t-x)+1) }
* accurate to about 1e-8
* SEE ALSO: ifdm12, ifd32, ifd52, fdm12, fd12, fd32, fd52
*/
{
mask = x < 4.;
list = where(mask);
if (numberof(list)) {
a = x(list);
a = log(a*poly(a, 1.999266880833e4, 5.702479099336e3,
6.610132843877e2, 3.818838129486e1, 1.0) /
poly(a, 1.771804140488e4, -2.014785161019e3,
9.130355392717e1, -1.670718177489e0));
}
list = where(!mask);
if (numberof(list)) {
x = x(list);
b = x^(-2./3.);
b = poly(b, -1.277060388085e-2, 7.187946804945e-2,
-4.262314235106e-1, 4.997559426872e-1,
-1.285579118012e0, -3.930805454272e-1, 1.0) /
(b*poly(b, -9.745794806288e-3, 5.485432756838e-2,
-3.299466243260e-1, 4.077841975923e-1,
-1.145531476975e0, -6.067091689181e-2));
}
return merge(a, b, mask);
}
func ifd32 (x)
/* DOCUMENT ifd32(y)
* return x = inverse of Fermi-Dirac integral of order 3/2,
* y = integral[0 to inf]{ dt * t^1.5 / (exp(t-x)+1) }
* accurate to about 1e-8
* SEE ALSO: ifdm12, ifd12, ifd52, fdm12, fd12, fd32, fd52
*/
{
mask = x < 4.;
list = where(mask);
if (numberof(list)) {
a = x(list);
a = log(a*poly(a, 1.715627994191e2, 1.125926232897e2,
2.056296753055e1, 1.0) /
poly(a, 2.280653583157e2, 1.193456203021e2,
1.167743113540e1, -3.226808804038e-1,
3.519268762788e-3));
}
list = where(!mask);
if (numberof(list)) {
x = x(list);
b = x^(-0.4);
b = poly(b, -6.321828169799e-3, -2.183147266896e-2,
-1.057562799320e-1, -4.657944387545e-1,
-5.951932864088e-1, 3.684471177100e-1, 1.0) /
(b*poly(b, -4.381942605018e-3, -1.513236504100e-2,
-7.850001283886e-2, -3.407561772612e-1,
-5.074812565486e-1, -1.387107009074e-1));
}
return merge(a, b, mask);
}
func ifd52 (x)
/* DOCUMENT ifd52(y)
* return x = inverse of Fermi-Dirac integral of order 5/2,
* y = integral[0 to inf]{ dt * t^2.5 / (exp(t-x)+1) }
* accurate to about 1e-8
* SEE ALSO: ifdm12, ifd12, ifd32, fdm12, fd12, fd32, fd52
*/
{
mask = x < 4.;
list = where(mask);
if (numberof(list)) {
a = x(list);
a = log(a*poly(a, 2.138969250409e2, 3.539903493971e1, 1.0) /
poly(a, 7.108545512710e2, 9.873746988121e1,
1.067755522895e0, -1.182798726503e-2));
}
list = where(!mask);
if (numberof(list)) {
x = x(list);
b = x^(-2./7.);
b = poly(b, -3.312041011227e-2, 1.315763372315e-1,
-4.820942898296e-1, 5.099038074944e-1,
5.495613498630e-1, -1.498867562255e0, 1.0) /
(b*poly(b, -2.315515517515e-2, 9.198776585252e-2,
-3.835879295548e-1, 5.415026856351e-1,
-3.847241692193e-1, 3.739781456585e-2,
-3.008504449098e-2));
}
return merge(a, b, mask);
}
