```/*
* 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);
}
```