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 ************************************************************************