FDp12_M := proc(x)
#
#  ---------------------------------------------------------------
#  
#  PURPOSE
#
#    This Maple 6 procedure Computes the Fermi-Dirac integral
#
#       Fj(x) = 1/GAMMMA(j+1)*int(t^j/(exp(t-x)+1),t=0..infinity)
#
#    of order j = 1/2, correct to at least 15 significant decimal
#    digits of accuracy using MacLeod's piecewise Chebyshev series
#    approximation scheme. (See reference below.)
#
#    The magnitude of the relative error is at most 0.4x10^(-15)
#    on the x domain (-infinity,infinity). 
#
#  USAGE
#
#    y = FDp12_M(x)
#
#  FILE NAMES 
#
#    FDp12_M.mpl    -- MacLeod's piecewise Chebyshev series approx 
#    Clenshaw_M.mpl -- Clenshaw's evaluation of Chebyshev sum
#
#  DATE
#
#    October 13, 2000
#
#  PROGRAMMER
#
#    Dr. F.G. Lether
#    Dept. Math. 
#    Univ. Georgia
#    Athens, GA. 30602 USA
#
#    email : fglether@math.uga.edu
#
#  REFERENCE
#
#    MacLeod, Allan.J.(1998), Algorithm 778: Fermi-Dirac Functions
#    of Order -1/2,1/2,3/2,5/2, Trans. Math. Soft. (TOMS) 24,1-12
#
#  ---------------------------------------------------------------
#  
   local twoe,const,c1,c2,c3,ex,t,xsq,ChebSum,FDp12_Cheb;
#
   Digits := 19;
#
#  Initialize constants 2*e and 1/GAMMA(5/2) = 4/(3*sqrt(Pi)) 
#
   twoe  := 5.4365636569180904e0;
   const := 7.5225277806367505e-1;
#   
#  Left tail region (-inf,-1) Chebyshev coefficients
#
   c1 := array(0..13 , [ 1.8862968392734597e0 ,
                        -0.543580817644053e-1 ,
                         0.23644975439720e-2  ,
                        -0.1216929365880e-3   ,
                         0.68695130622e-5     ,
                        -0.4112076172e-6      ,
                         0.256351628e-7       ,
                        -0.16465008e-8        ,
                         0.1081948e-9         ,
                        -0.72392e-11          ,
                         0.4915e-12           ,
                        -0.338e-13            ,
                         0.23e-14             ,           
                        -0.2e-15              ] );
#
#  Middle region I [-1,2) Chebyshev coefficients
#
   c2 := array(0..23 , [ 2.6982492788170612e0 ,
                         1.2389914141133012e0 ,
                         0.2291439379816278e0 ,
                         0.90316534687279e-2  ,
                        -0.25776524691246e-2  ,
                        -0.583681605388e-4    ,
                         0.693609458725e-4    ,
                        -0.18061670265e-5     ,
                        -0.21321530005e-5     ,
                         0.1754983951e-6      ,
                         0.665325470e-7       ,
                        -0.101675977e-7       , 
                        -0.19637597e-8        ,
                         0.5075769e-9         , 
                         0.491469e-10         , 
                        -0.233737e-10         ,
                        -0.66450e-12          ,
                         0.10115e-11          ,
                        -0.313e-13            ,
                        -0.412e-13            ,
                         0.38e-14             ,
                         0.16e-14             ,
                        -0.3e-15              ,
                        -0.1e-15              ] );
#
#  Right tail region [2,infinity) Chebyshev coefficients
#
   c3 := array(0..53 , [ 2.5484384198009122e0 ,
                         0.510439408960652e-1 ,
                         0.77493527628294e-2  ,
                        -0.75041656584953e-2  ,
                        -0.77540826320296e-2  , 
                        -0.45810844539977e-2  ,
                        -0.23431641587363e-2  ,
                        -0.11788049513591e-2  ,
                        -0.5802739359702e-3   ,
                        -0.2825350700537e-3   , 
                        -0.1388136651799e-3   ,
                        -0.680695084875e-4    ,
                        -0.335356350608e-4    ,
                        -0.166533018734e-4    ,
                        -0.82714908266e-5     ,
                        -0.41425714409e-5     ,
                        -0.20805255294e-5     ,   
                        -0.10479767478e-5     , 
                        -0.5315273802e-6      ,
                        -0.2694061178e-6      ,
                        -0.1374878749e-6      , 
                        -0.702308887e-7       ,
                        -0.359543942e-7       ,
                        -0.185106126e-7       ,
                        -0.95023937e-8        ,            
                        -0.49184811e-8        , 
                        -0.25371950e-8        ,
                        -0.13151532e-8        ,
                        -0.6835168e-9         ,
                        -0.3538244e-9         ,
                        -0.1853182e-9         ,
                        -0.958983e-10         ,
                        -0.504083e-10         , 
                        -0.262238e-10         , 
                        -0.137255e-10         ,
                        -0.72340e-11          , 
                        -0.37429e-11          ,
                        -0.20059e-11          , 
                        -0.10269e-11          ,
                        -0.5551e-12           ,
                        -0.2857e-12           ,
                        -0.1520e-12           , 
                        -0.811e-13            , 
                        -0.410e-13            , 
                        -0.234e-13            ,
                        -0.110e-13            , 
                        -0.67e-14             ,
                        -0.30e-14             ,
                        -0.19e-14             ,
                        -0.9e-15              , 
                        -0.5e-15              ,
                        -0.3e-15              , 
                        -0.1e-15              , 
                        -0.1e-15              ] ); 
#
#  Compute MacLeod approximation depending on x in
#  (-infinity,-1), [1,2) or [2,infinity)
#   
   if x < -1 then
#
     ex := evalf(exp(x)); 
     t :=  twoe*ex - 1;
     ChebSum := Clenshaw_M(c1,13,t);
     FDp12_Cheb := ex*ChebSum;
#
   elif x < 2 then
#
     t := (2*x-1)/3;
     ChebSum := Clenshaw_M(c2,23,t);
     FDp12_Cheb := ChebSum;
#
   else
#
     FDp12_Cheb := const*x*evalf(sqrt(x));
     xsq := x*x;
     t := (50 - xsq)/(42 + xsq);
     ChebSum := Clenshaw_M(c3,53,t);
     FDp12_Cheb := FDp12_Cheb*(1 + ChebSum/xsq);
   fi;
#
   end: