#  --------------------------------------------------------------
#
   FDp32_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= 3/2, correct to at least 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.26x10^(-8)
#    on the x domain (-infinity,infinity). 
#
#  USAGE
#
#    y = FDp32_CT(x)
#
#  FILE NAME 
#
#    FDp32_CT.txt
#
#  DATE
#
#    October 3, 2000
#
#  Author
#
#    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,FDp32;
   Digits := 13;
#
#  Initialize constant 1/GAMMA(5/2) = 4/(3*sqrt(Pi))
#	
   const := .7522527780637; 
#   
#  Region 1 : Left tail (-infinity,1] Cody-Thacher coefficients
#
   p1 := array(0..4 , [ -2.349963985406e-1 ,
                        -2.927373637547e-1 ,
                        -9.883097588738e-2 ,
                        -8.251386379551e-3 ,
                        -1.874384153223e-5 ] );
   q1 := array(0..4 , [  1.000000000000    ,
                         1.608597109146    ,
                         8.275289530880e-1 ,
                         1.522322382850e-1 ,
                         7.695120475064e-3 ] );
#
#  Region 2 : Middle region [1,4] Cody-Thacher coefficients
#
   p2 := array(0..4 , [ 1.1530213402    ,
                        1.0591558972    ,
                        4.6898803095e-1 ,
                        1.1882908784e-1 ,
                        1.9438755787e-2 ] );
   q2 := array(0..4 , [ 1.0000000000    ,
                        3.7348953841e-2 ,
                        2.3248458137e-2 ,
                       -1.3766770874e-3 ,
                        4.6466392781e-5 ] );
#
#  Region 3 : Right tail [4,infinity) Cody-Thacher coefficients
#
   p3 := array(0..4 , [ 2.46740023684   ,
                        2.19167582368e2 ,
                        1.23829379075e4 ,
                        2.20667724968e5  ,
                        8.49442920034e5  ] );
   q3 := array(0..4 , [ 1.000000000000   ,
                        8.91125140619e1 ,
                        5.04575669667e3 ,
                        9.09075946304e4 ,
                        3.89960915641e5 ] );
#
#  Explicit function to evaluate Horner nested form of a 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
#  (-inf,1), [1,4) or [4,inf). (Use Maple hardware floating
#  point evaluation instruction, evalhf, to gain speed.) 
#
   if x < 1 then
     t := evalf(exp(x));
     FDp32 := evalhf( t * ( 1 + const * t * R(p1,q1,t) ) );
   elif x < 4 then
     FDp32 := evalhf( const * R(p2,q2,x) );
   else
     t := evalhf(1/(x^2));
     FDp32 := evalhf( const * x*x*evalf(sqrt(x))
                      * ( 2/5 + t * R(p3,q3,t) ) );
   fi;
#
#  Round floating point result to Digits = 13 significant digits,
#  with at least 8 of these significant digits being correct
#
   FDp32 := evalf(FDp32,Digits);
#   
   end: