#  --------------------------------------------------------------
#
   FDp12_CT := proc(x::numeric)
#
#  --------------------------------------------------------------
#
#  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 nearly 8 significant decimal
#    digits of accuracy using Cody and Thacher's pecewise
#    rational approximation scheme. (See reference below.)
#
#    The magnitude of the relative error is at most 0.7x10^(-8)
#    on the x domain (-inf,inf). 
#
#  USAGE
#
#    y = FDp12_CT(x)
#
#  FILE NAME 
#
#    FDp12_CT.txt
#
#  DATE
#
#    October 3, 2000
#
#  Author
#
#    Dr. F.G. Lether
#    Dept. Math. 
#    Univ. Georgia
#    Athens, GA. 30602 USA
#
#    email: fglether@math.uga.edu
#
#  REFERENCE
#
#    Cody,W.J. and Thacher,H.C.(1967),Rational Chebyshev
#    Approximations for Fermi-Dirac Integrals of Orders -1/2,
#    1/2 and 3/2, Math. Comp. 21, 30-40.
#
#  --------------------------------------------------------------
#
   local const,p1,q1,p2,q2,p3,q3,R,t,FDp12;
   Digits := 13;
#
#  Initialize constant 1/GAMMA(3/2) = 2/sqrt(Pi)
#	
   const := 1.128379167096; 
#   
#  Region 1 : Left tail (-inf,1] Cody-Thacher coefficients
#
   p1 := array(0..4 , [ -3.133285305570e-1 ,
                        -4.161873852293e-1 ,
                        -1.502208400588e-1 ,
                        -1.339579375173e-2 ,
                        -1.513350700138e-5 ] );
   q1 := array(0..4 , [  1.000000000000    ,
                         1.872608675902    ,
                         1.145204446578    ,
                         2.570225587573e-1 ,
                         1.639902543568e-2 ] );
#
#  Region 2 : Middle region [1,4] Cody-Thacher coefficients
#
   p2 := array(0..4 , [ 6.78176626660e-1 ,
                        6.33124017910e-1 ,
                        2.94479651772e-1 ,
                        8.01320711419e-2 ,
                        1.33918212940e-2 ] );
   q2 := array(0..4 , [ 1.000000000000   ,
                        1.43740400397e-1 ,
                        7.08662148450e-2 ,
                        2.34579494735e-3 ,
                       -1.29449928835e-5 ] );
#
#  Region 3 : Right tail [4,inf) Cody-Thacher coefficients
#
   p3 := array(0..4 , [ 8.2244997626e-1 ,
                        2.0046303393e1  ,
                        1.8268093446e3  ,
                        1.2226530374e4  ,
                        1.4040750092e5  ] );
   q3 := array(0..4 , [ 1.000000000000   ,
                        2.3486207659e1 ,
                        2.2013483743e3 ,
                        1.1442673596e4 ,
                        1.6584715900e5 ] );
#
#  Explicit function to evaluate Horner nested form of rational
#  function Rn,n(z) = Pn(x)/Qn(z), with n = 4
#
   R := (P,Q,z) ->
        ( P[0] + (P[1] + (P[2] + (P[3] + P[4]*z)*z)*z)*z )
      / ( Q[0] + (Q[1] + (Q[2] + (Q[3] + Q[4]*z)*z)*z)*z );
#
#  Compute Cody-Thacher approxiamtion depending on x in
#  (-infinity,1), [1,4) or [4,infinity). (Use the Maple
#  hardware floating point evaluation instruction, evalhf,
#  to gain speed.)  
#
   if x < 1 then
     t := evalhf(exp(x));
     FDp12 := evalhf( t * ( 1 + const * t * R(p1,q1,t) ) );
   elif x < 4 then
     FDp12 := evalhf( const * R(p2,q2,x) );
   else
     t := evalhf(1/(x^2));
     FDp12 := evalhf( const * x * evalf(sqrt(x))
                            * ( 2/3 + t * R(p3,q3,t) ) );
   fi;
#
# Round floating point result to Digits = 13 significant digits
# with nearly 8 of these significant digits being correct
#
   FDp12 := evalf(FDp12,Digits);
   end: