FDp32_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 = 3/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 = FDp32_M(x)
#
#  FILE NAMES 
#
#    FDp32_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,FDp32_Cheb;
#
   Digits := 19;
#
#  Initialize constants 2*e and 1/GAMMA(7/2) = 8/(15*sqrt(Pi)) 
#
   twoe  := 5.4365636569180904e0;
   const := 3.0090111122547002e-1;
#   
#  Left tail region (-inf,-1) Chebyshev coefficients
#
   c1 := array(0..12 , [ 1.9406549210378650e0 ,
                        -0.287867475518043e-1 ,
                         0.8509157952313e-3   ,
                        -0.332784525669e-4    ,
                         0.15171202058e-5     ,
                        -0.762200874e-7       ,
                         0.40955489e-8        ,
                        -0.2311964e-9         ,
                         0.135537e-10         ,
                        -0.8187e-12           ,
                         0.507e-13            ,
                        -0.32e-14             ,
                         0.2e-15              ] );
#
#  Middle region I [-1,2) Chebyshev coefficients
#
   c2 := array(0..22 , [ 3.5862251615634306e0 ,
                         1.8518290056265751e0 ,
                         0.4612349102417150e0 ,
                         0.579303976126881e-1 ,
                         0.17043790554875e-2  ,
                        -0.3970520122496e-3   ,
                        -0.70702491890e-5     ,
                         0.76599748792e-5     , 
                        -0.1857811333e-6      ,
                        -0.1832237956e-6      , 
                         0.139249495e-7       ,
                         0.46702027e-8        , 
                        -0.6671984e-9         ,
                        -0.1161292e-9         ,
                         0.284438e-10         ,
                         0.24906e-11          ,
                        -0.11431e-11          , 
                        -0.279e-13            ,
                         0.439e-13            , 
                        -0.14e-14             , 
                        -0.16e-14             ,
                         0.1e-15              ,
                         0.1e-15              ] );
#
#  Right tail region [2,infinity) Chebyshev coefficients
#
   c3 := array(0..55 , [ 12.1307581736884627e0 ,
                         -0.1547501111287255e0 ,
                         -0.739007388850999e-1 , 
                         -0.307235377959258e-1 ,
                         -0.114548579330328e-1 ,
                         -0.40567636809539e-2  ,
                         -0.13980158373227e-2  ,
                         -0.4454901810153e-3   ,
                         -0.1173946112704e-3   ,
                         -0.148408980093e-4    ,
                          0.118895154223e-4    ,
                          0.146476958178e-4    ,
                          0.113228741730e-4    , 
                          0.75762292948e-5     , 
                          0.47120400466e-5     , 
                          0.28132628202e-5     ,
                          0.16370517341e-5     ,
                          0.9351076272e-6      ,
                          0.5278689210e-6      ,
                          0.2951079870e-6      , 
                          0.1638600190e-6      ,
                          0.905205409e-7       ,
                          0.497756975e-7       ,
                          0.272955863e-7       ,
                          0.149214585e-7       ,
                          0.81420359e-8        ,
                          0.44349200e-8        ,
                          0.24116032e-8        ,
                          0.13105018e-8        ,
                          0.7109736e-9         ,
                          0.3856721e-9         , 
                          0.2089529e-9         ,
                          0.1131735e-9         ,
                          0.612785e-10         ,
                          0.331448e-10         , 
                          0.179419e-10         ,
                          0.96953e-11          ,
                          0.52463e-11          ,
                          0.28343e-11          ,
                          0.15323e-11          ,
                          0.8284e-12           ,
                          0.4472e-12           , 
                          0.2421e-12           ,
                          0.1304e-12           ,
                          0.707e-13            , 
                          0.381e-13            ,
                          0.206e-13            , 
                          0.111e-13            ,
                          0.60e-14             ,
                          0.33e-14             ,
                          0.17e-14             ,
                          0.11e-14             ,
                          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,12,t);
     FDp32_Cheb := ex*ChebSum;
#
   elif x < 2 then
#
     t := (2*x-1)/3;
     ChebSum := Clenshaw_M(c2,22,t);
     FDp32_Cheb := ChebSum;
#
   else
#
     xsq := x*x;
     FDp32_Cheb := const*xsq*evalf(sqrt(x));
     t := (50 - xsq)/(42 + xsq);
     ChebSum := Clenshaw_M(c3,55,t);
     FDp32_Cheb := FDp32_Cheb*(1 + ChebSum/xsq);
   fi;
#
   end: