      function finv(f)       
*
* Compute Fermi-Dirac integral F1/2 to accuracy of 10^-4 using 
* rational function approximation of Anita (1993, ApJS, 84, 101).
*
      implicit double precision(a-h,o-z)
      dimension a1(3),b1(3),a2(3),b2(3)
*
* low precision formulae
*
      data an,m1,k1,m2,k2/0.5d0,2,2,2,2/
      data a1/4.4593646d+01,+1.1288764d+01,1.0000000d+00/
      data b1/3.9519346d+01,-5.7517464d+00,2.6594291d-01/
      data a2/3.4873722d+01,-2.6922515d+01,1.0000000d+00/
      data b2/2.6612832d+01,-2.0452930d+01,1.1808945d+01/
*
* small eta rational function       
*
      if(f .lt. 4.d0)then
        rn = f + a1(m1)
        do 100 i=m1-1,1,-1
          rn = rn*f + a1(i)
100     continue
        den = b1(k1+1)
        do 200 i=k1,1,-1
          den = den*f + b1(i)
200     continue
        finv = dlog(f * rn / den)
       else
*
* large eta rational function       
*
        ff = 1.d0 / f**(1.d0/(1.d0+an))
        rn = ff + a2(m2)
        do 300 i=m2-1,1,-1
          rn = rn*ff + a2(i)
300     continue
        den = b2(k2+1)
        do 400 i=k2,1,-1
          den = den*ff + b2(i)
400     continue
        finv =  rn / (den * ff)
      endif

      return      
      end 

************************************************************************

