program xfdh
!
! Test double precision piecewise minimax rational approximations
! of Fermi-Dirac integrals of half integer orders k=-9/2(1)21/2
! and of integer orders k=0(1)10
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
implicit integer (i-n)
implicit real*8 (a-h,o-z)
!
! approximate functions for half integer orders, k=-9/2(1)21/2
!
real*8 fdm9h,fdm7h,fdm5h,fdm3h,fdm1h
real*8 fd1h,fd3h,fd5h,fd7h,fd9h
real*8 fd11h,fd13h,fd15h,fd17h,fd19h
real*8 fd21h
!
! approximate functions for integer orders, k=0(1)10
!
real*8 fd0h,fd2h,fd4h,fd6h,fd8h
real*8 fd10h,fd12h,fd14h,fd16h,fd18h
real*8 fd20h
!
real*8 eta(8)	! test arguments
!
eta(1)=-3.d0
eta(2)=-1.d0
eta(3)=1.d0
eta(4)=3.5d0
eta(5)=7.5d0
eta(6)=15d0
eta(7)=30d0
eta(8)=60d0
!
do n=-9,21  ! n=2k
    if(n.EQ.-8.or.n.EQ.-6.or.n.EQ.-4.or.n.EQ.-2) cycle
    do j=1,8
        x=eta(j)
        if(n.EQ.-9) then
            fd=fdm9h(x)
        elseif(n.EQ.-7) then
            fd=fdm7h(x)
        elseif(n.EQ.-5) then
            fd=fdm5h(x)
        elseif(n.EQ.-3) then
            fd=fdm3h(x)
        elseif(n.EQ.-1) then
            fd=fdm1h(x)
        elseif(n.EQ.0) then
            fd=fd0h(x)
        elseif(n.EQ.1) then
            fd=fd1h(x)
        elseif(n.EQ.2) then
            fd=fd2h(x)
        elseif(n.EQ.3) then
            fd=fd3h(x)
        elseif(n.EQ.4) then
            fd=fd4h(x)
        elseif(n.EQ.5) then
            fd=fd5h(x)
        elseif(n.EQ.6) then
            fd=fd6h(x)
        elseif(n.EQ.7) then
            fd=fd7h(x)
        elseif(n.EQ.8) then
            fd=fd8h(x)
        elseif(n.EQ.9) then
            fd=fd9h(x)
        elseif(n.EQ.10) then
            fd=fd10h(x)
        elseif(n.EQ.11) then
            fd=fd11h(x)
        elseif(n.EQ.12) then
            fd=fd12h(x)
        elseif(n.EQ.13) then
            fd=fd13h(x)
        elseif(n.EQ.14) then
            fd=fd14h(x)
        elseif(n.EQ.15) then
            fd=fd15h(x)
        elseif(n.EQ.16) then
            fd=fd16h(x)
        elseif(n.EQ.17) then
            fd=fd17h(x)
        elseif(n.EQ.18) then
            fd=fd18h(x)
        elseif(n.EQ.19) then
            fd=fd19h(x)
        elseif(n.EQ.20) then
            fd=fd20h(x)
        elseif(n.EQ.21) then
            fd=fd21h(x)
        endif
        write (*,"(a1,0p2f5.1,1pe25.15)") "!",dble(n)*0.5d0,x,fd
    enddo
    write (*,"(a1)") "!"
enddo
!
end program xfdh
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fdm9h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=-9/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=-2.d0/7.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fdm9h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(0.270088205852269109d0 &
    -ex*(98457.373700674016d0 &
    +t*(7895.1716828987114d0 &
    +t*(1170.39154599810719d0 &
    +t*(84.347808600628057d0 &
    +t*(1.94779771326323869d0 &
    +t*0.00062804476066590963d0 &
    )))))/(32220.8981607867840d0 &
    +t*(20608.4545396913128d0 &
    +t*(5234.8569367171782d0 &
    +t*(658.75071261067250d0 &
    +t*(40.9380138615149699d0 &
    +t))))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    fd=(-1.8432713155273516d0 &
    +s*(-7.2420131297087391d0 &
    +s*(-2.2155513134917048d0 &
    +s*(0.04304761419364448d0 &
    +s*(-0.16668060134499780d0 &
    +s*(0.10403025949642735d0 &
    +s*(0.016505424797869711d0 &
    +s*(-0.0063406328104926592d0 &
    +s*(-0.0009078668865214013d0 &
    +s*0.00024692958093140298d0 &
    )))))))))/(129.8291336259148d0 &
    +s*(-1.9134679779378d0 &
    +s*(249.6123960853833d0 &
    +s*(-3.3454980273582d0 &
    +s*(190.7566926415911d0 &
    +s*(-2.2190711173347d0 &
    +s*(72.29929273955658d0 &
    +s*(-0.66382343556531d0 &
    +s*(13.551723146833901d0 &
    +s*(-0.0757657170023758d0 &
    +s))))))))))*(x+1.82715356570827096d0)	! care for zero point
elseif(x.lt.2.d0) then
    t=0.5d0*x
    fd=(27.9376225634323064d0 &
    +t*(17.6954702922794847d0 &
    +t*(-18.2160790482833106d0 &
    +t*(-2.57045608672749918d0 &
    +t*(2.20866272949862340d0 &
    +t*(-0.104202929779896185d0 &
    +t*(-0.177712759450918142d0 &
    +t*(0.055743070525419128d0 &
    -t*0.0059642251814009821d0 &
    ))))))))/(599.97601432112870d0 &
    +t*(-64.538219224509587d0 &
    +t*(1048.07497435809354d0 &
    +t*(-104.535376308139931d0 &
    +t*(694.31701140352634d0 &
    +t*(-61.402988847817226d0 &
    +t*(207.271703307329581d0 &
    +t*(-14.7241003912299059d0 &
    +t*(23.5639226214135345d0 &
    -t)))))))))*(x-0.557104412737456300d0)	! care for zero point
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=-(0.105736863491337914d0 &
    +t*(0.182686468618144972d0 &
    +t*(-0.0830333880807117100d0 &
    +t*(0.0123471104343183151d0 &
    +t*(0.0213438566082710205d0 &
    +t*(-0.0431206240183745280d0 &
    +t*(0.0311190249582737067d0 &
    +t*(-0.0115670162945002128d0 &
    +t*(0.00222720333535487902d0 &
    -t*0.000179769912955733772d0 &
    )))))))))/(11.4369029896970731d0 &
    +t*(29.5790675947434036d0 &
    +t*(55.2958932613908644d0 &
    +t*(59.7192644653455683d0 &
    +t*(51.8646589384828875d0 &
    +t*(26.5061101661226671d0 &
    +t*(11.1515306816966566d0 &
    +t*(-0.140562900980389991d0 &
    +t*(-0.566848363630388431d0 &
    -t)))))))))*(x-3.79898975457237574d0)	! care for zero point
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=-(0.0131615741211877996d0 &
    +t*(0.0630819955070757629d0 &
    +t*(0.0347025335431948214d0 &
    +t*(0.00659613983600166143d0 &
    +t*(0.0239129779624141424d0 &
    +t*(0.000528952248906432005d0 &
    +t*(0.00411353218448862483d0 &
    +t*(-0.00119448333536781374d0 &
    +t*(0.000497581590580634773d0 &
    +t*(-0.0000934288845342437109d0 &
    +t*7.29489294899966855d-6 &
    ))))))))))/(12.2828253912160684d0 &
    +t*(63.2143680817101004d0 &
    +t*(155.880993448869246d0 &
    +t*(239.067686638304706d0 &
    +t*(254.656949674257383d0 &
    +t*(200.134181826630681d0 &
    +t*(120.784221710133443d0 &
    +t*(56.5037249583191013d0 &
    +t*(20.2196452284175417d0 &
    +t*(5.22412387722958575d0 &
    +t))))))))))
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=-(0.000243271785585373458d0 &
    +t*(0.000534699823019827761d0 &
    +t*(0.000824373965146324914d0 &
    +t*(0.00106321349769688454d0 &
    +t*(0.000488890808530810197d0 &
    +t*(0.000510404811590977338d0 &
    +t*(0.0000543711407291589868d0 &
    +t*(0.0000597066095966549458d0 &
    +t*(-9.52409541937957252d-6 &
    +t*(1.26499979908255436d-6 &
    -t*8.23879904655758674d-8 &
    ))))))))))/(1.89001039755418574d0 &
    +t*(12.0291520111708943d0 &
    +t*(36.3016148236874403d0 &
    +t*(69.0194326307265861d0 &
    +t*(92.5327210822279397d0 &
    +t*(92.3443444983825834d0 &
    +t*(70.3437819943334403d0 &
    +t*(41.0546070181362401d0 &
    +t*(17.9813587502434782d0 &
    +t*(5.52467389311001165d0 &
    +t))))))))))
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=-(5.76414987722574846d-6 &
    +t*(0.0000175838621123388022d0 &
    +t*(0.0000306497936833891724d0 &
    +t*(0.0000245465851262614617d0 &
    +t*(9.68986921084327441d-6 &
    +t*(-9.60429091051474115d-6 &
    +t*(-2.71258276002580026d-6 &
    +t*(2.14306002093744082d-7 &
    +t*(-2.02878733208210684d-8 &
    +t*1.11972866705178940d-9 &
    )))))))))/(0.672666709110930353d0 &
    +t*(4.51114986457223474d0 &
    +t*(14.2130577502093159d0 &
    +t*(27.0924755033205922d0 &
    +t*(33.3161322251216687d0 &
    +t*(25.3582658290231325d0 &
    +t*(8.61719049364763904d0 &
    +t*(-2.98613950390747585d0 &
    +t*(-3.72974082078125757d0 &
    -t)))))))))
else
    w=1.d0/(x*x)
    s=1.d0-1600.d0*w
    fd=factor/(sqrt(x)*x*x*x)*(1.d0 &
    +w*(20360.4093941608469d0 &
    +s*(48690.0157659639679d0 &
    +s*(56148.8553220593957d0 &
    +s*(16506.5714937251887d0 &
    +s*(1066.30280063583769d0 &
    +s*3.73715435019239458d0 &
    )))))/(765.319028032818290d0 &
    +s*(1851.12566867527470d0 &
    +s*(2160.21729234460541d0 &
    +s*(677.261952266814024d0 &
    +s*(56.0179124057448609d0 &
    +s))))))
endif
!write(*,"(a20,1p2e15.7)") "(fdm9h) t,fd=",t,fd
fdm9h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fdm7h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=-7/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=-2.d0/5.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fdm7h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(-0.945308720482941881d0 &
    +ex*(22521.6817696480830d0 &
    +t*(3006.50943115266211d0 &
    +t*(223.777970320362539d0 &
    +t*(7.1096231598119406d0 &
    -t*0.00186000064211079826d0 &
    ))))/(4211.64894546505314d0 &
    +t*(2132.92301652550982d0 &
    +t*(400.935813755588393d0 &
    +t*(33.0304914363496116d0 &
    +t)))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    fd=(5.49736264373955718d0 &
    +s*(3.93868182790099018d0 &
    +s*(0.00680535603918657089d0 &
    +s*(-0.0610218211998865281d0 &
    +s*(-0.195616572333524455d0 &
    +s*(0.0172573301744899007d0 &
    +s*(0.0148839062617757320d0 &
    +s*(-0.00259857880173559174d0 &
    +s*0.0000170523346905109458d0 &
    ))))))))/(48.4238197836838959d0 &
    +s*(8.18206322265668686d0 &
    +s*(73.9370055146725832d0 &
    +s*(10.7170918787834126d0 &
    +s*(42.0907078084017798d0 &
    +s*(4.72976153773643018d0 &
    +s*(10.5965223842451046d0 &
    +s*(0.704336166377684061d0 &
    +s))))))))*(x+0.731435761340666539d0)	! care for zero point
elseif(x.lt.2.d0) then
    t=0.5d0*x
    fd=(143.837264704139569d0 &
    +t*(322.683997310438763d0 &
    +t*(-209.374575032955500d0 &
    +t*(-54.681932464048843d0 &
    +t*(-30.3617985514093593d0 &
    +t*(-5.7605448223972412d0 &
    +t*(3.65745825030136209d0 &
    +t*(0.145440718467921945d0 &
    +t*(-0.340001125998911237d0 &
    +t*0.0446577756904410773d0 &
    )))))))))/(1732.20737492411446d0 &
    +t*(97.952930820143076d0 &
    +t*(2660.50150777614871d0 &
    +t*(115.730278615026183d0 &
    +t*(1524.22026694958607d0 &
    +t*(40.4116469657132500d0 &
    +t*(385.907052068138736d0 &
    +t*(1.45068964456455148d0 &
    +t*(36.4544527882983039d0 &
    -t)))))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=-(66.1873851488861175d0 &
    +t*(57.9197865216898287d0 &
    +t*(-10.6658833177171402d0 &
    +t*(10.2066946495410176d0 &
    +t*(1.19595124198914370d0 &
    +t*(-1.81286879962085966d0 &
    +t*(1.52800227451442796d0 &
    +t*(-0.546498161708704000d0 &
    +t*(0.0988077919447287808d0 &
    -t*0.00745279173673929794d0 &
    )))))))))/(1525.13818882033488d0 &
    +t*(3966.07893650424194d0 &
    +t*(6837.27734004268945d0 &
    +t*(7086.45653543153192d0 &
    +t*(5612.71274301484627d0 &
    +t*(2928.56588821761611d0 &
    +t*(1256.17847516062581d0 &
    +t*(310.558929999645079d0 &
    +t*(89.3091775022809646d0 &
    -t)))))))))*(x-2.59355543171509487d0)	! care for zero point
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=-(0.184179377827614791d0 &
    +t*(0.516542442597467614d0 &
    +t*(0.292799259677321999d0 &
    +t*(0.166303092691678597d0 &
    +t*(0.258615034442190300d0 &
    +t*(-0.0303157249245734253d0 &
    +t*(0.0373447714795897992d0 &
    +t*(-0.00332926057033243781d0 &
    +t*(0.000563408857595472728d0 &
    -t*0.0000334037183021934960d0 &
    )))))))))/(18.2615126816087432d0 &
    +t*(85.1686926876657081d0 &
    +t*(181.378990080239764d0 &
    +t*(232.528454347497368d0 &
    +t*(200.564064511137621d0 &
    +t*(125.206874911765040d0 &
    +t*(58.6883968744415564d0 &
    +t*(20.6768210567814309d0 &
    +t*(5.02620698086580395d0 &
    +t)))))))))
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=-(0.00521027824483942692d0 &
    +t*(0.0159223834333424812d0 &
    +t*(0.0290786393625134859d0 &
    +t*(0.0350554456604475612d0 &
    +t*(0.0276022017688728517d0 &
    +t*(0.0190056432579866130d0 &
    +t*(0.00725681379911524042d0 &
    +t*(0.00278845414163336780d0 &
    +t*(0.000285395695659601419d0 &
    +t*(-9.95541819385254380d-6 &
    +t*2.96666654254966665d-7 &
    ))))))))))/(3.40047012308625195d0 &
    +t*(20.3896764592442369d0 &
    +t*(58.0983301256082527d0 &
    +t*(104.466101980605120d0 &
    +t*(132.574327309248600d0 &
    +t*(125.279491711147983d0 &
    +t*(90.3637198347191251d0 &
    +t*(49.8942487159252830d0 &
    +t*(20.6309250519127594d0 &
    +t*(5.96507377885482873d0 &
    +t))))))))))
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=-(0.000222755850099679261d0 &
    +t*(0.00115329637370955141d0 &
    +t*(0.00290802444182764526d0 &
    +t*(0.00445062601336644735d0 &
    +t*(0.00436605816073874863d0 &
    +t*(0.00260703930450196953d0 &
    +t*(0.000708204675930321953d0 &
    +t*(0.0000387281377930932522d0 &
    +t*(-1.07635110532534644d-6 &
    +t*3.45876932703497152d-8 &
    )))))))))/(0.958823449511661332d0 &
    +t*(7.43981964725756497d0 &
    +t*(27.2010614567095770d0 &
    +t*(61.4593095996041459d0 &
    +t*(94.2239246717949674d0 &
    +t*(100.959076669111257d0 &
    +t*(74.4305466149867624d0 &
    +t*(35.3635458846636560d0 &
    +t*(9.43423871492993288d0 &
    +t)))))))))
else
    w=1.d0/(x*x)
    s=1.d0-1600.d0*w
    fd=factor/(sqrt(x)*x*x)*(1.d0 &
    +w*(114389.966467649306d0 &
    +s*(109888.157408212290d0 &
    +s*(21209.9082486576965d0 &
    +s*(896.625771647944984d0 &
    +s*2.27005555539139462d0 &
    ))))/(7803.86522209355981d0 &
    +s*(7642.62657004470466d0 &
    +s*(1584.72834559231955d0 &
    +s*(86.0988387732258338d0 &
    +s)))))
endif
!write(*,"(a20,1p2e15.7)") "(fdm7h) t,fd=",t,fd
fdm7h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fdm5h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=-5/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=-2.d0/3.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fdm5h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(2.36327180120735470d0 &
    -ex*(55891.891151428516d0 &
    +t*(10900.4507673238614d0 &
    +t*(731.61871282829311d0 &
    +t*(17.3576202287841595d0 &
    +t*0.00111269690244199535d0 &
    ))))/(8361.6144419623308d0 &
    +t*(3709.66537451151168d0 &
    +t*(598.60597628606205d0 &
    +t*(41.1772384480908825d0 &
    +t)))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    fd=(46.7273884990861909d0 &
    +s*(65.1460094720913271d0 &
    +s*(-23.3560330148569366d0 &
    +s*(16.3832748811699383d0 &
    +s*(-6.19159773777140396d0 &
    +s*(0.799041627235852088d0 &
    +s*(-0.136702158227027886d0 &
    +s*(0.0592594484852497570d0 &
    -s*0.00813537888028213710d0 &
    ))))))))/(166.600811683081505d0 &
    +s*(-14.3475443594245037d0 &
    +s*(207.623210767457661d0 &
    +s*(-16.0651771000195284d0 &
    +s*(91.9265825199149767d0 &
    +s*(-5.91543082309762067d0 &
    +s*(16.6946542439476255d0 &
    +s*(-0.714551880206870314d0 &
    +s))))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    fd=-(55.096201117757232d0 &
    +t*(14.3294333176185961d0 &
    +t*(8.7057675255700556d0 &
    +t*(2.85717540287481278d0 &
    +t*(0.153440611811510950d0 &
    +t*(0.111339313377987928d0 &
    +t*(0.052394928476702755d0 &
    -t*0.00308189450941126788d0 &
    )))))))/(217.840649625384406d0 &
    +t*(-13.7539547994891024d0 &
    +t*(264.475360271691825d0 &
    +t*(-16.3936568402444332d0 &
    +t*(112.575929318729487d0 &
    +t*(-6.5735021347402422d0 &
    +t*(19.1324946771515682d0 &
    +t*(-0.89572129073702419d0 &
    +t))))))))*(x-1.10894923342229868d0)	! care for zero point
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=-(4.91753613607304828d0 &
    +t*(20.4599278677959969d0 &
    +t*(14.7299603843029675d0 &
    +t*(6.07034924563989291d0 &
    +t*(2.84749597355562421d0 &
    +t*(1.42067667370006629d0 &
    +t*(0.194025010869198784d0 &
    +t*(0.0296942990922311325d0 &
    +t*0.0000687311336894144097d0 &
    ))))))))/(39.1942895728295075d0 &
    +t*(102.720679402185990d0 &
    +t*(163.814228574815100d0 &
    +t*(160.042363559709804d0 &
    +t*(112.736106158227044d0 &
    +t*(53.5881251269350506d0 &
    +t*(19.3013585080264391d0 &
    +t*(4.60099930165203599d0 &
    +t))))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=-(2.27928384600199307d0 &
    +t*(4.20924981693908136d0 &
    +t*(3.71772880855517432d0 &
    +t*(3.45602273475008183d0 &
    +t*(1.35886950179549082d0 &
    +t*(0.536408369715153184d0 &
    +t*(0.133773439887674841d0 &
    +t*(0.0125358162552656152d0 &
    -t*0.000258659892241001720d0 &
    ))))))))/(29.5538409791709387d0 &
    +t*(102.889070328418668d0 &
    +t*(171.483234327128671d0 &
    +t*(174.069526985568413d0 &
    +t*(121.638517061339171d0 &
    +t*(61.2471517251153565d0 &
    +t*(22.6513730231817499d0 &
    +t*(5.63585485331701445d0 &
    +t))))))))
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=-(0.138865047109769882d0 &
    +t*(0.495879942014454703d0 &
    +t*(0.975841030964797014d0 &
    +t*(1.19756546203286062d0 &
    +t*(1.04262591465861665d0 &
    +t*(0.636815947333856586d0 &
    +t*(0.267431185397286764d0 &
    +t*(0.0680541575349993884d0 &
    +t*(0.00339987875007404643d0 &
    -t*0.0000596473600392971721d0 &
    )))))))))/(6.10307238045677639d0 &
    +t*(32.0683844393622288d0 &
    +t*(81.7709446404731603d0 &
    +t*(131.907213444156949d0 &
    +t*(148.807381678792024d0 &
    +t*(122.126271732385001d0 &
    +t*(73.4713378422568434d0 &
    +t*(31.3211204388750577d0 &
    +t*(8.63324324695581677d0 &
    +t)))))))))
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=-(0.294313222065093082d0 &
    +t*(1.31885924825320721d0 &
    +t*(2.84488318326892294d0 &
    +t*(3.47529649317584998d0 &
    +t*(2.44599317338319296d0 &
    +t*(0.711594188260176831d0 &
    +t*(0.0568475228400688404d0 &
    +t*(0.000310240841471278615d0 &
    +t*7.25941094961112814d-6 &
    ))))))))/(38.8567648628867585d0 &
    +t*(233.714122704911173d0 &
    +t*(657.094399788013655d0 &
    +t*(1097.58958431940912d0 &
    +t*(1157.52830939122072d0 &
    +t*(742.054932624983009d0 &
    +t*(251.716945386093114d0 &
    +t*(34.4183681886106367d0 &
    +t))))))))
else
    w=1.d0/(x*x)
    s=1.d0-1600.d0*w
    fd=factor/(sqrt(x)*x)*(1.d0 &
    +w*(563249.531577994933d0 &
    +s*(348717.634501234702d0 &
    +s*(40765.2886083808789d0 &
    +s*865.059113888543852d0 &
    )))/(90262.6676538587293d0 &
    +s*(56944.6907876652433d0 &
    +s*(7175.67113380361324d0 &
    +s*(207.376531486457753d0 &
    +s)))))
endif
!write(*,"(a20,1p2e15.7)") "(fdm5h) t,fd=",t,fd
fdm5h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fdm3h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=-3/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=-2.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fdm3h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(-3.54490770181103205d0 &
    +ex*(82737.595643818605d0 &
    +t*(18481.5553495836940d0 &
    +t*(1272.73919064487495d0 &
    +t*(26.3420403338352574d0 &
    -t*0.00110648970639283347d0 &
    ))))/(16503.7625405383183d0 &
    +t*(6422.0552658801394d0 &
    +t*(890.85389683932154d0 &
    +t*(51.251447078851450d0 &
    +t)))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=-(946.638483706348559d0 &
    +t*(76.3328330396778450d0 &
    +t*(62.7809183134124193d0 &
    +t*(83.8442376534073219d0 &
    +t*(23.2285755924515097d0 &
    +t*(3.21516808559640925d0 &
    +t*(1.58754232369392539d0 &
    +t*(0.687397326417193593d0 &
    +t*0.111510355441975495d0 &
    ))))))))/(889.4123665319664d0 &
    +s*(126.7054690302768d0 &
    +s*(881.4713137175090d0 &
    +s*(108.2557767973694d0 &
    +s*(289.38131234794585d0 &
    +s*(27.75902071820822d0 &
    +s*(34.252606975067480d0 &
    +s*(1.9592981990370705d0 &
    +s))))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    fd=-(754.61690882095729d0 &
    +t*(565.56180911009650d0 &
    +t*(494.901267018948095d0 &
    +t*(267.922900418996927d0 &
    +t*(110.418683240337860d0 &
    +t*(39.4050164908951420d0 &
    +t*(10.8654460206463482d0 &
    +t*(2.11194887477009033d0 &
    +t*0.246843599687496060d0 &
    ))))))))/(560.03894899770103d0 &
    +t*(70.007586553114572d0 &
    +t*(582.42052644718871d0 &
    +t*(56.181678606544951d0 &
    +t*(205.248662395572799d0 &
    +t*(12.5169006932790528d0 &
    +t*(27.2916998671096202d0 &
    +t*(0.53299717876883183d0 &
    +t))))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=-(526.022770226139287d0 &
    +t*(631.116211478274904d0 &
    +t*(516.367876532501329d0 &
    +t*(267.894697896892166d0 &
    +t*(91.3331816844847913d0 &
    +t*(17.5723541971644845d0 &
    +t*(1.46434478819185576d0 &
    +t*(1.29615441010250662d0 &
    +t*0.223495452221465265d0 &
    ))))))))/(354.867400305615304d0 &
    +t*(560.931137013002977d0 &
    +t*(666.070260050472570d0 &
    +t*(363.745894096653220d0 &
    +t*(172.272943258816724d0 &
    +t*(23.7751062504377332d0 &
    +t*(12.5916012142616255d0 &
    +t*(-0.888604976123420661d0 &
    +t))))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=-(18.0110784494455205d0 &
    +t*(36.1225408181257913d0 &
    +t*(38.4464752521373310d0 &
    +t*(24.1477896166966673d0 &
    +t*(9.27772356782901602d0 &
    +t*(2.49074754470533706d0 &
    +t*(0.163824586249464178d0 &
    -t*0.00329391807590771789d0 &
    )))))))/(18.8976860386360201d0 &
    +t*(49.3696375710309920d0 &
    +t*(60.9273314194720251d0 &
    +t*(43.6334649971575003d0 &
    +t*(20.6568810936423065d0 &
    +t*(6.11094689399482273d0 &
    +t))))))
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=-(4.10698092142661427d0 &
    +t*(17.1412152818912658d0 &
    +t*(32.6347877674122945d0 &
    +t*(36.6653101837618939d0 &
    +t*(25.9424894559624544d0 &
    +t*(11.2179995003884922d0 &
    +t*(2.30099511642112478d0 &
    +t*(0.0928307248942099967d0 &
    -t*0.00146397877054988411d0 &
    ))))))))/(6.40341731836622598d0 &
    +t*(30.1333068545276116d0 &
    +t*(64.0494725642004179d0 &
    +t*(80.5635003792282196d0 &
    +t*(64.9297873014508805d0 &
    +t*(33.3013900893183129d0 &
    +t*(9.61549304470339929d0 &
    +t)))))))
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=-(95.2141371910496454d0 &
    +t*(420.050572604265456d0 &
    +t*(797.778374374075796d0 &
    +t*(750.378359146985564d0 &
    +t*(324.818150247463736d0 &
    +t*(50.3115388695905757d0 &
    +t*(0.372431961605507103d0 &
    +t*(-0.103162211894757911d0 &
    +t*0.00191752611445211151d0 &
    ))))))))/(212.232981736099697d0 &
    +t*(1043.79079070035083d0 &
    +t*(2224.50099218470684d0 &
    +t*(2464.84669868672670d0 &
    +t*(1392.55318009810070d0 &
    +t*(346.597189642259199d0 &
    +t*(22.7314613168652593d0 &
    -t)))))))
else
    w=1.d0/(x*x)
    s=1.d0-1600.d0*w
    fd=factor/sqrt(x)*(1.d0 &
    +w*(12264.3569103180524d0 &
    +s*(3204.34872454052352d0 &
    +s*(140.119604748253961d0 &
    +s*0.523918919699235590d0 &
    )))/(9877.87829948067200d0 &
    +s*(2644.71979353906092d0 &
    +s*(128.863768007644572d0 &
    +s))))
endif
!write(*,"(a20,1p2e15.7)") "(fdm3h) t,fd=",t,fd
fdm3h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fdm1h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=-1/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fdm1h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(1.77245385090551603d0 &
    -ex*(40641.4537510284430d0 &
    +t*(9395.7080940846442d0 &
    +t*(649.96168315267301d0 &
    +t*(12.7972295804758967d0 &
    +t*0.00153864350767585460d0 &
    ))))/(32427.1884765292940d0 &
    +t*(11079.9205661274782d0 &
    +t*(1322.96627001478859d0 &
    +t*(63.738361029333467d0 &
    +t)))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(272.770092131932696d0 &
    +t*(30.8845653844682850d0 &
    +t*(-6.43537632380366113d0 &
    +t*(14.8747473098217879d0 &
    +t*(4.86928862842142635d0 &
    +t*(-1.53265834550673654d0 &
    +t*(-1.02698898315597491d0 &
    +t*(-0.177686820928605932d0 &
    -t*0.00377141325509246441d0 &
    ))))))))/(293.075378187667857d0 &
    +s*(305.818162686270816d0 &
    +s*(299.962395449297620d0 &
    +s*(207.640834087494249d0 &
    +s*(92.0384803181851755d0 &
    +s*(37.0164914112791209d0 &
    +s*(7.88500950271420583d0 &
    +s)))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    fd=(3531.50360568243046d0 &
    +t*(6077.5339658420037d0 &
    +t*(6199.7700433981326d0 &
    +t*(4412.78701919567594d0 &
    +t*(2252.27343092810898d0 &
    +t*(811.84098649224085d0 &
    +t*(191.836401053637121d0 &
    +t*23.2881838959183802d0 &
    )))))))/(3293.83702584796268d0 &
    +t*(1528.97474029789098d0 &
    +t*(2568.48562814986046d0 &
    +t*(925.64264653555825d0 &
    +t*(574.23248354035988d0 &
    +t*(132.803859320667262d0 &
    +t*(29.8447166552102115d0 &
    +t)))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=(4060.70753404118265d0 &
    +t*(10812.7291333052766d0 &
    +t*(13897.5649482242583d0 &
    +t*(10628.4749852740029d0 &
    +t*(5107.70670190679021d0 &
    +t*(1540.84330126003381d0 &
    +t*(284.452720112970331d0 &
    +t*29.5214417358484151d0 &
    )))))))/(1564.58195612633534d0 &
    +t*(2825.75172277850406d0 &
    +t*(3189.16066169981562d0 &
    +t*(1955.03979069032571d0 &
    +t*(828.000333691814748d0 &
    +t*(181.498111089518376d0 &
    +t*(32.0352857794803750d0 &
    +t)))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=(1198.41719029557508d0 &
    +t*(3263.51454554908654d0 &
    +t*(3874.97588471376487d0 &
    +t*(2623.13060317199813d0 &
    +t*(1100.41355637121217d0 &
    +t*(267.469532490503605d0 &
    +t*(25.4207671812718340d0 &
    +t*0.389887754234555773d0 &
    )))))))/(273.407957792556998d0 &
    +t*(595.918318952058643d0 &
    +t*(605.202452261660849d0 &
    +t*(343.183302735619981d0 &
    +t*(122.187622015695729d0 &
    +t*(20.9016359079855933d0 &
    +t))))))
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=(9446.00169435237637d0 &
    +t*(36843.4448474028632d0 &
    +t*(63710.1115419926191d0 &
    +t*(62985.2197361074768d0 &
    +t*(37634.5231395700921d0 &
    +t*(12810.9898627807754d0 &
    +t*(1981.56896138920963d0 &
    +t*81.4930171897667580d0 &
    )))))))/(1500.04697810133666d0 &
    +t*(5086.91381052794059d0 &
    +t*(7730.01593747621895d0 &
    +t*(6640.83376239360596d0 &
    +t*(3338.99590300826393d0 &
    +t*(860.499043886802984d0 &
    +t*(78.8565824186926692d0 &
    +t)))))))
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=(22977.9657855367223d0 &
    +t*(123416.616813887781d0 &
    +t*(261153.765172355107d0 &
    +t*(274618.894514095795d0 &
    +t*(149710.718389924860d0 &
    +t*(40129.3371700184546d0 &
    +t*(4470.46495881415076d0 &
    +t*132.684346831002976d0 &
    )))))))/(2571.68842525335676d0 &
    +t*(12521.4982290775358d0 &
    +t*(23268.1574325055341d0 &
    +t*(20477.2320119758141d0 &
    +t*(8726.52577962268114d0 &
    +t*(1647.42896896769909d0 &
    +t*(106.475275142076623d0 &
    +t)))))))
else
    w=1.d0/(x*x)
    t=1600.d0*w
    fd=sqrt(x)*factor*(1.d0 &
    -w*(0.411233516712009968d0 &
    +t*(0.00110980410034088951d0 &
    +t*(0.0000113689298990173683d0 &
    +t*(2.56931790679436797d-7 &
    +t*(9.97897786755446178d-9 &
    +t*8.67667698791108582d-10))))))
endif
!write(*,"(a20,1p2e15.7)") "(fdm1h) t,fd=",t,fd
fdm1h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd0h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=0
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,w,s,fd
!
!write(*,"(a20,1pe15.7)") "(fd0h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(1.d0 &
    -ex*(22696.2126132366633d0 &
    +t*(5222.0667923565138d0 &
    +t*(357.623326425354522d0 &
    +t*(6.9167792879948140d0 &
    +t*0.00200096064827815813d0 &
    ))))/(45392.4252264733267d0 &
    +t*(14539.5980679273792d0 &
    +t*(1611.36476693109675d0 &
    +t*(71.072178562726798d0 &
    +t)))))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(159.601717762460980d0 &
    +t*(23.7193942338278703d0 &
    +t*(0.377783268730614356d0 &
    +t*(10.5181677709577503d0 &
    +t*(3.78181326142271599d0 &
    +t*(-0.441998676614933572d0 &
    +t*(-0.450072959113928254d0 &
    +t*(-0.0734798777056723512d0 &
    +t*0.000915454570009894267d0 &
    ))))))))/(284.26032127745967d0 &
    +s*(315.2592651624449d0 &
    +s*(310.2713981221035d0 &
    +s*(206.21640678892182d0 &
    +s*(96.77898293084927d0 &
    +s*(35.456591489081173d0 &
    +s*(8.1762315442738975d0 &
    +s)))))))
endif
if(y.gt.0.d0) then
    fd=fd+y
endif
!write(*,"(a20,1p2e15.7)") "(fd0h) y,fd=",y,fd
fd0h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd1h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=1/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/3.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd1h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(0.886226925452758014d0 &
    -ex*(19894.4553386951666d0 &
    +t*(4509.64329955948557d0 &
    +t*(303.461789035142376d0 &
    +t*(5.7574879114754736d0 &
    +t*0.00275088986849762610d0 &
    ))))/(63493.915041308052d0 &
    +t*(19070.1178243603945d0 &
    +t*(1962.19362141235102d0 &
    +t*(79.250704958640158d0 &
    +t)))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(149.462587768865243d0 &
    +t*(22.8125889885050154d0 &
    +t*(-0.629256395534285422d0 &
    +t*(9.08120441515995244d0 &
    +t*(3.35357478401835299d0 &
    +t*(-0.473677696915555805d0 &
    +t*(-0.467190913556185953d0 &
    +t*(-0.0880610317272330793d0 &
    -t*0.00262208080491572673d0 &
    ))))))))/(269.94660938022644d0 &
    +s*(343.6419926336247d0 &
    +s*(323.9049470901941d0 &
    +s*(218.89170769294024d0 &
    +s*(102.31331350098315d0 &
    +s*(36.319337289702664d0 &
    +s*(8.3317401231389461d0 &
    +s)))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    fd=(71652.717119215557d0 &
    +t*(134954.734070223743d0 &
    +t*(153693.833350315645d0 &
    +t*(123247.280745703400d0 &
    +t*(72886.293647930726d0 &
    +t*(32081.2499422362952d0 &
    +t*(10210.9967337762918d0 &
    +t*(2152.71110381320778d0 &
    +t*232.906588165205042d0 &
    ))))))))/(105667.839854298798d0 &
    +t*(31946.0752989314444d0 &
    +t*(71158.788776422211d0 &
    +t*(15650.8990138187414d0 &
    +t*(13521.8033657783433d0 &
    +t*(1646.98258283527892d0 &
    +t*(618.90691969249409d0 &
    +t*(-3.36319591755394735d0 &
    +t))))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=(23744.8706993314289d0 &
    +t*(68257.8589855623002d0 &
    +t*(89327.4467683334597d0 &
    +t*(62766.3415600442563d0 &
    +t*(20093.6622609901994d0 &
    +t*(-2213.89084119777949d0 &
    +t*(-3901.66057267577389d0 &
    -t*948.642895944858861d0 &
    )))))))/(9488.61972919565851d0 &
    +t*(12514.8125526953073d0 &
    +t*(9903.44088207450946d0 &
    +t*(2138.15420910334305d0 &
    +t*(-528.394863730838233d0 &
    +t*(-661.033633995449691d0 &
    +t*(-51.4481470250962337d0 &
    +t)))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=(311337.452661582536d0 &
    +t*(1.11267074416648198d6 &
    +t*(1.75638628895671735d6 &
    +t*(1.59630855803772449d6 &
    +t*(910818.935456183774d0 &
    +t*(326492.733550701245d0 &
    +t*(65507.2624972852908d0 &
    +t*4809.45649527286889d0 &
    )))))))/(39721.6641625089685d0 &
    +t*(86424.7529107662431d0 &
    +t*(88163.7255252151780d0 &
    +t*(50615.7363511157353d0 &
    +t*(17334.9774805008209d0 &
    +t*(2712.13170809042550d0 &
    +t*(82.2205828354629102d0 &
    -t)))))))*0.999999999999999877d0
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=(7.26870063003059784d6 &
    +t*(2.79049734854776025d7 &
    +t*(4.42791767759742390d7 &
    +t*(3.63735017512363365d7 &
    +t*(1.55766342463679795d7 &
    +t*(2.97469357085299505d6 &
    +t*154516.447031598403d0 &
    ))))))/(340542.544360209743d0 &
    +t*(805021.468647620047d0 &
    +t*(759088.235455002605d0 &
    +t*(304686.671371640343d0 &
    +t*(39289.4061400542309d0 &
    +t*(582.426138126398363d0 &
    +t*(11.2728194581586028d0 &
    -t)))))))
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=(4.81449797541963104d6 &
    +t*(1.85162850713127602d7 &
    +t*(2.77630967522574435d7 &
    +t*(2.03275937688070624d7 &
    +t*(7.41578871589369361d6 &
    +t*(1.21193113596189034d6 &
    +t*63211.9545144644852d0 &
    ))))))/(80492.7765975237449d0 &
    +t*(189328.678152654840d0 &
    +t*(151155.890651482570d0 &
    +t*(48146.3242253837259d0 &
    +t*(5407.08878394180588d0 &
    +t*(112.195044410775577d0 &
    -t))))))
else
    w=1.d0/(x*x)
    s=1.d0-1600.d0*w
    fd=x*sqrt(x)*0.666666666666666667d0*(1.d0+w &
    *(8109.79390744477921d0 &
    +s*(342.069867454704106d0 &
    +s*1.07141702293504595d0)) &
    /(6569.98472532829094d0 &
    +s*(280.706465851683809d0 &
    +s)))
endif
!write(*,"(a20,1p2e15.7)") "(fd1h) t,fd=",t,fd
fd1h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd2h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=1
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,w,s,fd
!
!write(*,"(a20,1pe15.7)") "(fd2h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(1.d0 &
    -ex*(22189.1070807945062d0 &
    +t*(4915.92700908746777d0 &
    +t*(322.901386168881348d0 &
    +t*(5.9897442965804548d0 &
    +t*0.00397641173774375092d0 &
    ))))/(88756.428323178025d0 &
    +t*(25002.3197546553836d0 &
    +t*(2389.06277237306633d0 &
    +t*(88.376214553692756d0 &
    +t)))))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(145.488167182330098d0 &
    +t*(251.392824471576922d0 &
    +t*(56.6537141912783024d0 &
    +t*(17.9918985363509694d0 &
    +t*(20.1369115558099802d0 &
    +t*(7.09659390228556164d0 &
    +t*(0.199701180197912643d0 &
    +t*(-0.403173132925886253d0 &
    -t*0.0792966701498222697d0 &
    ))))))))/(606.0757707716040d0 &
    +s*(374.1806357435014d0 &
    +s*(252.1367051536344d0 &
    +s*(27.2746245830016d0 &
    +s*(-61.57766112137513d0 &
    +s*(-53.72117554363975d0 &
    +s*(-25.678454878692950d0 &
    +s*(-7.1995819520154718d0 &
    -s))))))))
endif
if(y.gt.0.d0) then
    fd=-fd+1.64493406684822644d0 &
    +0.5d0*y*y
endif
!write(*,"(a20,1p2e15.7)") "(fd2h) y,fd=",y,fd
fd2h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd3h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=3/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/5.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd3h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(1.32934038817913702d0 &
    -ex*(1346.14119046566636d0 &
    +t*(199.946876779712712d0 &
    +t*(6.5210149677288048d0 &
    +t*0.0108588591982722183d0 &
    )))/(5728.3481201778541d0 &
    +t*(1132.17837281710987d0 &
    +t*(64.805243148002602d0 &
    +t))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(631.667081787115831d0 &
    +t*(504.131655805666135d0 &
    +t*(113.449065431934917d0 &
    +t*(56.0939647772947784d0 &
    +t*(43.3374223200846752d0 &
    +t*(12.8047010597109577d0 &
    +t*(0.219164386586949410d0 &
    +t*(-0.678659552658390139d0 &
    -t*0.126533769309899232d0 &
    ))))))))/(1180.5112183558028d0 &
    +s*(1101.0159189871135d0 &
    +s*(864.4448234404281d0 &
    +s*(392.2018227840790d0 &
    +s*(89.58093202779063d0 &
    +s*(-9.95066218572899d0 &
    +s*(-17.312068771997626d0 &
    +s*(-6.4116162917822773d0 &
    -s))))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    fd=(90122.488639370400d0 &
    +t*(157095.208147064037d0 &
    +t*(166879.962599668589d0 &
    +t*(125708.597728045460d0 &
    +t*(69968.278213181390d0 &
    +t*(29035.3292989055404d0 &
    +t*(8736.4439472398517d0 &
    +t*(1747.16784760309227d0 &
    +t*180.132410666734053d0 &
    ))))))))/(78176.777123671727d0 &
    +t*(-1681.44633240543085d0 &
    +t*(38665.7913035496031d0 &
    +t*(-2527.29685826087874d0 &
    +t*(5062.6683078100048d0 &
    +t*(-553.21165462054589d0 &
    +t*(165.395637981775430d0 &
    +t*(- 18.0295465153725544d0 &
    +t))))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=(912944.432058014054d0 &
    +t*(3.28217091334054338d6 &
    +t*(5.59250227196369585d6 &
    +t*(5.76136129685687470d6 &
    +t*(3.84331519034749983d6 &
    +t*(1.65284168824947710d6 &
    +t*(423452.676670436605d0 &
    +t*49835.4127241373113d0 &
    )))))))/(164873.145721762182d0 &
    +t*(257442.511191094986d0 &
    +t*(225604.160532840884d0 &
    +t*(99932.1955662320024d0 &
    +t*(24761.0878784286761d0 &
    +t*(1398.26392212830777d0 &
    +t*(-36.4450237523474167d0 &
    +t)))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=(1.88412548327216052d6 &
    +t*(8.08838896259910792d6 &
    +t*(1.56869793001790529d7 &
    +t*(1.79109792599373447d7 &
    +t*(1.31345142328147214d7 &
    +t*(6.29500412046744325d6 &
    +t*(1.89326213154091054d6 &
    +t*(312372.643127575407d0 &
    +t*18814.7420442630170d0 &
    ))))))))/(67768.3347951202583d0 &
    +t*(147635.914444221358d0 &
    +t*(151908.303165069423d0 &
    +t*(86671.1222110642970d0 &
    +t*(27855.9481608626219d0 &
    +t*(3833.22697473114940d0 &
    +t*(98.3384567064269554d0 &
    -t)))))))*0.999999999999999876d0	! correction to remove bias
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=(1.59656593348660977d9 &
    +t*(7.32769737561517060d9 &
    +t*(1.42662658588280191d10 &
    +t*(1.51238422045169918d10 &
    +t*(9.27233604548095476d9 &
    +t*(3.18834513406577423d9 &
    +t*(5.36061988605886123d8 &
    +t*3.03619219668246382d7 &
    )))))))/(1.18906980815759995d7 &
    +t*(2.62209219322122975d7 &
    +t*(2.28143701746618729d7 &
    +t*(8.57156701742181874d6 &
    +t*(1.13860063870524239d6 &
    +t*(27091.7884208687379d0 &
    +t*(-275.664733379090447d0 &
    +t)))))))*0.999999999999999829d0	! correction to remove bias
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=(2.60437581212904589d8 &
    +t*(1.08771546307370080d9 &
    +t*(1.81531350939088943d9 &
    +t*(1.52833764636304939d9 &
    +t*(6.70684451492750149d8 &
    +t*(1.40870639531414149d8 &
    +t*1.04957900377463854d7 &
    ))))))/(358448.871166784200d0 &
    +t*(611808.419702466190d0 &
    +t*(326307.561591723775d0 &
    +t*(58407.9904827573816d0 &
    +t*(2049.50040323021794d0 &
    +t*(-39.8767861209088081d0 &
    +t))))))*0.999999999999999828d0
else
    w=1.d0/(x*x)
    s=1.d0-1600.d0*w
    fd=x*x*sqrt(x)*factor*(1.d0 &
    +w*(6.16739021212286242d0 &
    +s*(0.00111530123694574981d0 &
    +s*(-2.79156524536560815d-6 &
    +s*(2.95571462110856359d-8 &
    -s*6.70917556862133933d-10)))))
endif
!write(*,"(a20,1p2e15.7)") "(fd3h) t,fd=",t,fd
fd3h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd4h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,w,s,fd
!
!write(*,"(a20,1pe15.7)") "(fd4h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(2.d0 &
    -ex*(1914.06748184935743d0 &
    +t*(273.085756700981399d0 &
    +t*(8.5861610217850095d0 &
    +t*0.0161890243763741414d0 &
    )))/(7656.2699273974454d0 &
    +t*(1399.35442210906621d0 &
    +t*(72.929152915475392d0 &
    +t))))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(2711.49678259128843d0 &
    +t*(1299.85460914884154d0 &
    +t*(222.606134197895041d0 &
    +t*(172.881855215582924d0 &
    +t*(112.951038040682055d0 &
    +t*(24.0376482128898634d0 &
    +t*(-2.68393549333878715d0 &
    +t*(-2.14077421411719935d0 &
    -t*0.326188299771397236d0 &
    ))))))))/(2517.1726659917047d0 &
    +s*(3038.7689794575778d0 &
    +s*(2541.7823512372631d0 &
    +s*(1428.0589853413436d0 &
    +s*(531.62378035996132d0 &
    +s*(122.54595216479181d0 &
    +s*(8.395768655115050d0 &
    +s*(-3.9142702096919080d0 &
    -s))))))))
endif
if(y.gt.0.d0) then
    fd=fd+y*(3.28986813369645287d0 &
    +0.333333333333333333d0*y*y)
endif
!write(*,"(a20,1p2e15.7)") "(fd4h) y,fd=",y,fd
fd4h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd5h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=5/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/7.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd5h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(3.32335097044784255d0 &
    -ex*(3004.41138112148735d0 &
    +t*(409.582975848055860d0 &
    +t*(12.3422465543050559d0 &
    +t*0.0252600128832787650d0 &
    )))/(10227.9400759349389d0 &
    +t*(1729.21912667445507d0 &
    +t*(82.075224736463695d0 &
    +t))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(3273.17630052701392d0 &
    +t*(1698.67676090446982d0 &
    +t*(498.075924016264615d0 &
    +t*(286.399893838802205d0 &
    +t*(153.449411299292925d0 &
    +t*(46.3172989461474687d0 &
    +t*(7.47269262560340991d0 &
    +t*0.512665223813025153d0 &
    )))))))/(1934.7654167995427d0 &
    +s*(2174.8582387970533d0 &
    +s*(1803.7716559637348d0 &
    +s*(970.4325125469826d0 &
    +s*(368.29988963760778d0 &
    +s*(94.85059993048974d0 &
    +s*(14.845914179579222d0 &
    +s)))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    s=1.d0-t
    fd=(46746.0667679244140d0 &
    +t*(67556.152887188131d0 &
    +t*(58097.414724032516d0 &
    +t*(35613.1718933770642d0 &
    +t*(15859.9235610267359d0 &
    +t*(5040.9166211327297d0 &
    +t*(1064.44778209849372d0 &
    +t*117.741363279815688d0 &
    )))))))/(13126.266942469915d0 &
    +s*(-1313.6945396119670d0 &
    +s*(2636.3767383046264d0 &
    +s*(453.2433336267639d0 &
    +s*(201.68980019364685d0 &
    +s*(50.35944043583614d0 &
    +s*(9.319958361334161d0 &
    +s)))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=(2.96446396463084245d6 &
    +t*(1.13985063437515144d7 &
    +t*(2.03551842990364919d7 &
    +t*(2.18469168036723457d7 &
    +t*(1.52194204967854489d7 &
    +t*(6.89564335059367232d6 &
    +t*(1.88985082880325956d6 &
    +t*243639.263893338434d0 &
    )))))))/(169113.644493386066d0 &
    +t*(249598.556346488591d0 &
    +t*(162470.264413303189d0 &
    +t*(48830.6689634566512d0 &
    +t*(4080.41491690869679d0 &
    +t*(-227.369223251313984d0 &
    +t*(19.1547319985914212d0 &
    -t)))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=(5.98759698946959779d7 &
    +t*(2.87159410448355925d8 &
    +t*(5.99247230055005978d8 &
    +t*(7.05261579500244084d8 &
    +t*(5.04521860843522470d8 &
    +t*(2.18416589668256432d8 &
    +t*(5.25509373690209805d7 &
    +t*5.33490422476777307d6 &
    )))))))/(469653.962099017877d0 &
    +t*(972160.483958338969d0 &
    +t*(696829.387839904810d0 &
    +t*(200021.997093042021d0 &
    +t*(12604.5784254924892d0 &
    +t*(-464.256227099766735d0 &
    +t*(26.4866799762414562d0 &
    -t)))))))*0.999999999999999786d0	! correction to remove bias
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=(1.11471644837527339d9 &
    +t*(5.74225610798369390d9 &
    +t*(1.26551657157853213d10 &
    +t*(1.53960861052359422d10 &
    +t*(1.10946777927365321d10 &
    +t*(4.68865671930554474d9 &
    +t*(1.05763202432846973d9 &
    +t*9.49124183370767767d7 &
    )))))))/(1.07734938834750844d6 &
    +t*(2.05459707311873616d6 &
    +t*(1.39824422108531691d6 &
    +t*(347716.997197363113d0 &
    +t*(17894.5194245484999d0 &
    +t*(-553.162195184268593d0 &
    +t*(28.1090136251865326d0 &
    -t)))))))*0.999999999999999759d0	! correction to remove bias
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=(8.02670881104191218d9 &
    +t*(3.67306963017003546d10 &
    +t*(6.83085243661572356d10 &
    +t*(6.55107149365528043d10 &
    +t*(3.37406319128317264d10 &
    +t*(8.67814222875818408d9 &
    +t*8.43844503352450216d8 &
    ))))))/(757905.984443885971d0 &
    +t*(868424.806294231233d0 &
    +t*(260670.917865642513d0 &
    +t*(14550.0472712579662d0 &
    +t*(-476.120164041067762d0 &
    +t*(25.8288614974100332d0 &
    -t))))))*0.999999999999999820d0	! correction to remove bias
else
    w=1.d0/(x*x)
    t=1600.d0*w
    fd=x*x*x*sqrt(x)*factor*(1.d0 &
    +w*(14.3931730849220041d0 &
    +t*(0.00776862867834285253d0 &
    +t*(3.78966458769690333d-6 &
    +t*(2.09248095155530095d-8 &
    +t*3.52097438532254351d-10)))))
endif
!write(*,"(a20,1p2e15.7)") "(fd5h) t,fd=",t,fd
fd5h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd6h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=3
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,w,s,fd,y2
!
!write(*,"(a20,1pe15.7)") "(fd6h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(6.d0 &
    -ex*(5121.6401850302408d0 &
    +t*(664.28706260743472d0 &
    +t*(19.0856927562699544d0 &
    +t*0.0410982603688952131d0 &
    )))/(13657.7071600806539d0 &
    +t*(2136.54222460571183d0 &
    +t*(92.376788603062645d0 &
    +t))))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(7881.24597452900838d0 &
    +t*(4323.07526636309661d0 &
    +t*(1260.13125873282465d0 &
    +t*(653.359212389160499d0 &
    +t*(354.630774329461644d0 &
    +t*(113.373708671587772d0 &
    +t*(19.9559488532742796d0 &
    +t*1.59407954898394322d0 &
    )))))))/(2570.7250703533430d0 &
    +s*(2972.7443644211129d0 &
    +s*(2393.9995533270879d0 &
    +s*(1259.0724833462608d0 &
    +s*(459.86413596901097d0 &
    +s*(112.60906419590854d0 &
    +s*(16.468882811659000d0 &
    +s)))))))
endif
if(y.gt.0.d0) then
    y2=y*y
    fd=-fd+11.3643939539669510d0 &
    +y2*(4.93480220054467931d0 &
    +y2*0.25d0)
endif
!write(*,"(a20,1p2e15.7)") "(fd6h) y,fd=",y,fd
fd6h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd7h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=7/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/9.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd7h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(11.6317283965674489d0 &
    -ex*(9371.8868848378853d0 &
    +t*(1152.33581795824468d0 &
    +t*(31.4225568934398918d0 &
    +t*0.069508514261022902d0 &
    )))/(18231.3053891121107d0 &
    +t*(2639.60073595942887d0 &
    +t*(103.985056337236794d0 &
    +t))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(12412.9547163940264d0 &
    +t*(6557.39898920194282d0 &
    +t*(1745.01251050433134d0 &
    +t*(883.903030170203178d0 &
    +t*(485.680209936338180d0 &
    +t*(150.395267497504100d0 &
    +t*(23.6841796284116177d0 &
    +t*(1.07158463967401877d0 &
    -t*0.157435982722068683d0 &
    ))))))))/(1990.3886647182248d0 &
    +s*(2450.5836860178091d0 &
    +s*(1959.3508463399357d0 &
    +s*(1039.6469110978007d0 &
    +s*(381.70374149314543d0 &
    +s*(94.69566180593596d0 &
    +s*(14.401931782580504d0 &
    +s)))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    s=1.d0-t
    fd=(138724.961458493789d0 &
    +t*(169546.564556791848d0 &
    +t*(125336.240767414531d0 &
    +t*(68469.193938175092d0 &
    +t*(27692.3811560029092d0 &
    +t*(8054.8076265209314d0 &
    +t*(1569.98144698357220d0 &
    +t*162.669853437143155d0 &
    )))))))/(7624.675802050973d0 &
    +s*(2096.5891309081106d0 &
    +s*(1735.7563869500611d0 &
    +s*(662.0787128606726d0 &
    +s*(218.91333929478294d0 &
    +s*(55.27329667089387d0 &
    +s*(9.904579892966869d0 &
    +s)))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=(2.20638302461716245d7 &
    +t*(7.19733683266519997d7 &
    +t*(1.11705692645922359d8 &
    +t*(1.06328515554623708d8 &
    +t*(6.68969329369869135d7 &
    +t*(2.78374130823093789d7 &
    +t*(7.11915156250133512d6 &
    +t*870011.051674696596d0 &
    )))))))/(311792.108500755594d0 &
    +t*(206109.613235058227d0 &
    +t*(81814.2443130826819d0 &
    +t*(3513.24062956494342d0 &
    +t*(368.155702298381570d0 &
    +t*(-106.091905490453452d0 &
    +t*(15.0961925962965432d0 &
    -t)))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=(3.0709767203984838d9 &
    +t*(1.3145330402583295d10 &
    +t*(2.5004908399964440d10 &
    +t*(2.7299220224056906d10 &
    +t*(1.8396311263825673d10 &
    +t*(7.6098692716640943d9 &
    +t*(1.7746200575457532d9 &
    +t*1.7740961200347283d8 &
    )))))))/(4.46811391083410494d6 &
    +t*(4.62189039302323844d6 &
    +t*(1.60932873023944955d6 &
    +t*(113458.965071032695d0 &
    +t*(-4751.36769214919833d0 &
    +t*(351.424157693641750d0 &
    +t*(-24.0392624870048634d0 &
    +t)))))))*0.999999999999999779d0	! correction to remove bias
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=(4.91615192873221833d9 &
    +t*(2.30941143669878184d10 &
    +t*(4.52380364776557202d10 &
    +t*(4.68413409944225048d10 &
    +t*(2.65134543816397678d10 &
    +t*(7.35961621586797944d9 &
    +t*(5.52546586238663078d8 &
    -t*7.67154234216335832d7 &
    )))))))/(550749.688306094769d0 &
    +t*(352810.964448809157d0 &
    +t*(12162.0407442978353d0 &
    +t*(-4414.91208383702609d0 &
    +t*(288.821012224214490d0 &
    +t*(-21.7191096434783647d0 &
    +t))))))*0.999999999999999647d0	! correction to remove bias
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=(1.23356074447405045d12 &
    +t*(6.8241092859415704d12 &
    +t*(1.5939237639946117d13 &
    +t*(2.0277830436278028d13 &
    +t*(1.5074920317497838d13 &
    +t*(6.4878414656043323d12 &
    +t*(1.4761927839220043d12 &
    +t*1.3410778049868313d11 &
    )))))))/(7.2811112110966361d6 &
    +t*(8.4187326609915276d6 &
    +t*(2.59806064027658564d6 &
    +t*(153974.958146639672d0 &
    +t*(-5666.5977917092119d0 &
    +t*(387.205694265311079d0 &
    +t*(-25.1500172248070365d0 &
    +t)))))))*0.999999999999999750d0	! correction to remove bias
else
    w=1.d0/(x*x)
    t=1600.d0*w
    fd=x*x*x*x*sqrt(x)*factor*(1.d0 &
    +w*(25.9077115528595532d0 &
    +t*(0.0699176581158797143d0 &
    +t*(-0.0000113689722343055157d0 &
    +t*(-2.69205155161558844d-8 &
    -t*2.81838487282327867d-10)))))
endif
!write(*,"(a20,1p2e15.7)") "(fd7h) t,fd=",t,fd
fd7h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd8h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=4
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,w,s,fd,y2
!
!write(*,"(a20,1pe15.7)") "(fd8h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(24.d0 &
    -ex*(18247.2542465629199d0 &
    +t*(2120.56302902849207d0 &
    +t*(54.659507299984584d0 &
    +t*0.121876197098273914d0 &
    )))/(24329.6723287505678d0 &
    +t*(3261.01909656925521d0 &
    +t*(117.071576489684022d0 &
    +t))))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(33093.9102025608137d0 &
    +t*(19031.4783357798975d0 &
    +t*(5431.08357152245897d0 &
    +t*(2393.94262931609398d0 &
    +t*(1268.40017257978070d0 &
    +t*(418.662172475927160d0 &
    +t*(77.4108960876508190d0 &
    +t*6.67374311450268357d0 &
    )))))))/(2645.4885670047153d0 &
    +s*(3236.2237702166948d0 &
    +s*(2500.5977847175497d0 &
    +s*(1278.0109577275445d0 &
    +s*(448.99020896813485d0 &
    +s*(105.86020755838874d0 &
    +s*(15.216887271751039d0 &
    +s)))))))
endif
if(y.gt.0.d0) then
    y2=y*y
    fd=fd+y*(45.4575758158678040d0 &
    +y2*(6.57973626739290575d0 &
    +y2*0.2d0))
endif
!write(*,"(a20,1p2e15.7)") "(fd8h) y,fd=",y,fd
fd8h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd9h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=9/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/11.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd9h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(52.3427777845535202d0 &
    -ex*(37544.8291647986266d0 &
    +t*(4113.48125371067712d0 &
    +t*(99.866490020209337d0 &
    +t*0.221035794126711913d0 &
    )))/(32460.7344732907336d0 &
    +t*(4028.81223496037193d0 &
    +t*(131.831032502212893d0 &
    +t))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(69576.0137991453185d0 &
    +t*(40945.8214269928341d0 &
    +t*(11636.1644381151941d0 &
    +t*(4627.86250549447612d0 &
    +t*(2350.13392862571985d0 &
    +t*(772.917515779084181d0 &
    +t*(143.282768049011730d0 &
    +t*12.4161720203166775d0 &
    )))))))/(2535.8419026162019d0 &
    +s*(3176.5677917330904d0 &
    +s*(2403.1400307314004d0 &
    +s*(1205.9132705211921d0 &
    +s*(416.75320914741449d0 &
    +s*(97.44802030226960d0 &
    +s*(14.192502951043942d0 &
    +s)))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    s=1.d0-t
    fd=(514822.53904173208d0 &
    +t*(522785.81416122782d0 &
    +t*(325144.232801546845d0 &
    +t*(158667.557186950475d0 &
    +t*(59144.454873935605d0 &
    +t*(15870.4677825005999d0 &
    +t*(2843.37449801975925d0 &
    +t*272.600692278110437d0 &
    )))))))/(4614.934746888508d0 &
    +s*(2755.8568288489921d0 &
    +s*(1644.4943294235255d0 &
    +s*(712.2181577760027d0 &
    +s*(238.45583147179604d0 &
    +s*(59.96344339125020d0 &
    +s*(10.461476374827578d0 &
    +s)))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=(3.42730803487596668d7 &
    +t*(9.75806658972896811d7 &
    +t*(1.34864192616165924d8 &
    +t*(1.16686466196565430d8 &
    +t*(6.80744250678140765d7 &
    +t*(2.67765233037235380d7 &
    +t*(6.59736845425507999d6 &
    +t*792820.676193592762d0 &
    )))))))/(98882.7648048876841d0 &
    +t*(8989.62788188007021d0 &
    +t*(9880.17356654342362d0 &
    +t*(-2706.34081341857522d0 &
    +t*(689.741654284866509d0 &
    +t*(-130.873077409052406d0 &
    +t*(16.2961910923971688d0 &
    -t)))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=(4.89775289337907328d9 &
    +t*(2.18308560190839931d10 &
    +t*(4.37929196886803725d10 &
    +t*(5.12376841913497741d10 &
    +t*(3.77548677711700977d10 &
    +t*(1.75249365659518690d10 &
    +t*(4.74974437725451705d9 &
    +t*5.81264585063986499d8 &
    )))))))/(1.16603582458148748d6 &
    +t*(904387.524377664561d0 &
    +t*(128592.404626994618d0 &
    +t*(-10482.2282242040897d0 &
    +t*(1464.66129358886087d0 &
    +t*(-193.730497094764856d0 &
    +t*(19.0689920747890676d0 &
    -t)))))))*0.999999999999999823d0	! correction to remove bias
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=(2.2676209606423816d13 &
    +t*(1.24113302112340516d14 &
    +t*(2.9428117821474437d14 &
    +t*(3.9045588584587502d14 &
    +t*(3.1161725148614816d14 &
    +t*(1.4855808500445768d14 &
    +t*(3.8726601091729722d13 &
    +t*4.1673615587969372d12 &
    )))))))/(2.72286381952640596d8 &
    +t*(1.76997321474048891d8 &
    +t*(1.58695998922133285d7 &
    +t*(-757839.84310180434d0 &
    +t*(61065.315196529561d0 &
    +t*(-4351.64276565798757d0 &
    +t*(174.373773100542411d0 &
    +t)))))))*0.999999999999999697d0	! correction to remove bias
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=(1.65627579413599954d13 &
    +t*(9.8868092265530109d13 &
    +t*(2.51632614382698224d14 &
    +t*(3.53003554017750602d14 &
    +t*(2.93707971643967604d14 &
    +t*(1.44188847117098389d14 &
    +t*(3.83794461491062424d13 &
    +t*4.22247796386900843d12 &
    )))))))/(5.7640009494074862d6 &
    +t*(3.82118971474499980d6 &
    +t*(375185.320814949919d0 &
    +t*(-21189.7277758964890d0 &
    +t*(2225.13154291024219d0 &
    +t*(-240.494104779578652d0 &
    +t*(20.7062411414387229d0 &
    -t)))))))*0.999999999999999697d0	! correction to remove bias
else
    w=1.d0/(x*x)
    t=1600.d0*w
    fd=x*x*x*x*x*sqrt(x)*factor*(1.d0 &
    +w*(40.7121181544936151d0 &
    +t*(0.256364746421960192d0 &
    +t*(0.000125058641885852090d0 &
    +t*(5.92447825215879480d-8 &
    +t*3.38237202703194112d-10)))))
endif
!write(*,"(a20,1p2e15.7)") "(fd9h) t,fd=",t,fd
fd9h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd10h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=5
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,s,fd,y2
!
!write(*,"(a20,1pe15.7)") "(fd10h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(120.d0 &
    -ex*(81190.938912603315d0 &
    +t*(8368.4990332049831d0 &
    +t*(190.753830813273698d0 &
    +t*0.413800735538960261d0 &
    )))/(43301.8340867217726d0 &
    +t*(4977.68099709243407d0 &
    +t*(148.484432990224213d0 &
    +t))))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(159651.547840031111d0 &
    +t*(96307.2005742063042d0 &
    +t*(26923.5693307648389d0 &
    +t*(9274.54751848368696d0 &
    +t*(4445.76333033698006d0 &
    +t*(1461.45267097859337d0 &
    +t*(272.164427980501432d0 &
    +t*23.6526046298891619d0 &
    )))))))/(2522.7839609396783d0 &
    +s*(3244.5527999477567d0 &
    +s*(2403.0532924519756d0 &
    +s*(1176.7202478443275d0 &
    +s*(397.7596246691212d0 &
    +s*(91.84661231161838d0 &
    +s*(13.491911254479298d0 &
    +s)))))))
endif
if(y.gt.0.d0) then
    y2=y*y
    fd=-fd+236.532261911384425d0 &
    +y2*(113.643939539669510d0 &
    +y2*(8.22467033424113218d0 &
    +y2*0.166666666666666667d0))
endif
!write(*,"(a20,1p2e15.7)") "(fd10h) y,fd=",y,fd
fd10h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd11h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=11/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/13.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd11h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(287.885277815044361d0 &
    -ex*(183706.438585548927d0 &
    +t*(17781.2198813711389d0 &
    +t*(379.468701417596835d0 &
    +t*0.79823396490661890d0 &
    )))/(57756.370490854840d0 &
    +t*(6150.6214356632731d0 &
    +t*(167.282750039076672d0 &
    +t))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(399729.728975269086d0 &
    +t*(245048.290430557943d0 &
    +t*(64937.4257386015642d0 &
    +t*(17956.2886271385592d0 &
    +t*(8102.86054746075370d0 &
    +t*(2739.56619011305740d0 &
    +t*(522.197272134368319d0 &
    +t*45.9656786040236010d0 &
    )))))))/(2594.1469579416305d0 &
    +s*(3458.1847294224010d0 &
    +s*(2522.2748339106920d0 &
    +s*(1201.9024264852580d0 &
    +s*(394.94834682669296d0 &
    +s*(89.38046041124415d0 &
    +s*(13.110503750759739d0 &
    +s)))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    s=1.d0-t
    fd=(2.40644550206125935d6 &
    +t*(2.04903536564989627d6 &
    +t*(1.06804765545666824d6 &
    +t*(478219.980682203232d0 &
    +t*(171969.059741023943d0 &
    +t*(43927.0876916260199d0 &
    +t*(7256.0914773211089d0 &
    +t*626.34457422853533d0 &
    )))))))/(3109.1225672233219d0 &
    +s*(2636.2757909353112d0 &
    +s*(1626.1798939585643d0 &
    +s*(743.1451899485297d0 &
    +s*(255.41213934754826d0 &
    +s*(64.38943570079975d0 &
    +s*(10.994511861529934d0 &
    +s)))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    fd=(7.70852968624936390d8 &
    +t*(2.17991233625412353d9 &
    +t*(3.06662111607913460d9 &
    +t*(2.78373598102826860d9 &
    +t*(1.77557231061500545d9 &
    +t*(8.12988256861068453d8 &
    +t*(2.60886497846886753d8 &
    +t*(5.37757509332685362d7 &
    +t*5.49029147083874521d6 &
    ))))))))/(384975.653810726401d0 &
    +t*(-10859.8352547471390d0 &
    +t*(47234.6496298591730d0 &
    +t*(-15689.8880335197160d0 &
    +t*(4580.68161072037214d0 &
    +t*(-1041.20834072320464d0 &
    +t*(172.979919201040627d0 &
    +t*(-18.7593019718622089d0 &
    +t))))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=(1.71681959469725307d10 &
    +t*(7.11557910478070659d10 &
    +t*(1.35103963154256158d11 &
    +t*(1.52010553996923224d11 &
    +t*(1.09288906897533869d11 &
    +t*(5.01817283608162539d10 &
    +t*(1.36404590384355901d10 &
    +t*1.69961340800732270d9 &
    )))))))/(600169.964027912926d0 &
    +t*(63989.8050690151317d0 &
    +t*(3306.49827047230707d0 &
    +t*(-2248.31637427985521d0 &
    +t*(649.075943625483938d0 &
    +t*(-126.914632653258205d0 &
    +t*(16.0289311039308149d0 &
    -t)))))))*0.999999999999999735d0
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=(1.05700178515928121d12 &
    +t*(6.12117465680051512d12 &
    +t*(1.55582880442596747d13 &
    +t*(2.24736526624247277d13 &
    +t*(1.99069164118763497d13 &
    +t*(1.08056362404911750d13 &
    +t*(3.32622572257609455d12 &
    +t*4.47682895372084249d11 &
    )))))))/(1.27888301811042482d6 &
    +t*(318603.811469994946d0 &
    +t*(-33834.0645458632405d0 &
    +t*(6000.96391500182524d0 &
    +t*(-1070.68701259288375d0 &
    +t*(161.612080255393227d0 &
    +t*(-17.5331213409665917d0 &
    +t)))))))*0.999999999999999653d0
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=(1.06038260567473087d14 &
    +t*(6.8238792169646891d14 &
    +t*(1.88917697801412952d15 &
    +t*(2.91401892625408501d15 &
    +t*(2.70150682400302739d15 &
    +t*(1.50286942831009197d15 &
    +t*(4.63446459437813772d14 &
    +t*6.0885766731916815d13 &
    )))))))/(2.08723140369481691d6 &
    +t*(445823.690247689681d0 &
    +t*(-42239.1125174643989d0 &
    +t*(6891.8137509933955d0 &
    +t*(-1158.79785029794968d0 &
    +t*(168.017431473136672d0 &
    +t*(-17.7743748531881933d0 &
    +t)))))))*0.999999999999999637d0
else
    w=1.d0/(x*x)
    t=1600.d0*w
    fd=x*x*x*x*x*x*sqrt(x)*factor*(1.d0 &
    +w*(58.8063928898240827d0 &
    +t*(0.666548340699108862d0 &
    +t*(0.00162576228417182612d0 &
    +t*(-2.56764926438640722d-7 &
    -t*6.18964665583859548d-10)))))
endif
!write(*,"(a20,1p2e15.7)") "(fd11h) t,fd=",t,fd
fd11h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd12h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=6
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,s,fd,y2
!
!write(*,"(a20,1pe15.7)") "(fd12h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(720.d0 &
    -ex*(433290.557356514403d0 &
    +t*(39323.1678283287030d0 &
    +t*(783.71684376947655d0 &
    +t*1.58412947146158337d0 &
    )))/(77029.432418935897d0 &
    +t*(7600.9245809611507d0 &
    +t*(188.511069473956679d0 &
    +t))))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(1.06809887479312876d6 &
    +t*(651074.246191348755d0 &
    +t*(152197.899924352192d0 &
    +t*(27269.0203707062592d0 &
    +t*(11937.9088600476726d0 &
    +t*(4659.76125900467198d0 &
    +t*(964.100791156809939d0 &
    +t*88.5841245838029230d0 &
    )))))))/(2681.3731718905701d0 &
    +s*(3764.9048490469408d0 &
    +s*(2739.9504946219358d0 &
    +s*(1276.0933863294022d0 &
    +s*(406.93880737411049d0 &
    +s*(89.70143752313466d0 &
    +s*(12.997281214703279d0 &
    +s)))))))
endif
if(y.gt.0.d0) then
    y2=y*y
    fd=fd+y*(1419.19357146830655d0 &
    +y2*(227.287879079339020d0 &
    +y2*(9.86960440108935862d0 &
    +y2*0.142857142857142857d0)))
endif
!write(*,"(a20,1p2e15.7)") "(fd12h) y,fd=",y,fd
fd12h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd13h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=13/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/15.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd13h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    s=1.d0-t
    fd=ex*(1871.25430579778835d0 &
    -ex*(10.2714824358192133d0 &
    +s*(0.0648393897767320738d0 &
    +s*(0.000971410002224865156d0 &
    +s*(0.0000232470542763426993d0 &
    +s*(7.41058598210549775d-7 &
    +s*3.37424008840868627d-8))))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(2.94504307654357499d6 &
    +t*(1.72142318914178000d6 &
    +t*(302664.670131641659d0 &
    +t*(-4392.91160671565775d0 &
    +t*(1886.14265293074511d0 &
    +t*(5292.33536016064607d0 &
    +t*(1532.39748567462704d0 &
    +t*160.488909859013319d0 &
    )))))))/(2672.0262186872157d0 &
    +s*(4050.0534640335863d0 &
    +s*(2998.5673483066873d0 &
    +s*(1381.1929813352472d0 &
    +s*(429.69401904358381d0 &
    +s*(92.16501861997658d0 &
    +s*(13.091741877025778d0 &
    +s)))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    s=1.d0-t
    fd=(1.45543624525939873d7 &
    +t*(1.09901378809591560d7 &
    +t*(5.0278309600981774d6 &
    +t*(2.17920053476107569d6 &
    +t*(801203.84933034017d0 &
    +t*(204535.562039585858d0 &
    +t*(32251.6713113453623d0 &
    +t*2520.49966644549104d0 &
    )))))))/(2525.1970865635211d0 &
    +s*(2519.0633694107075d0 &
    +s*(1640.0689226392724d0 &
    +s*(779.5069098915619d0 &
    +s*(273.74080946976357d0 &
    +s*(69.07956623994342d0 &
    +s*(11.541634429054885d0 &
    +s)))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    s=1.d0-t
    fd=(1.93509493140212933d9 &
    +t*(4.73031979570464235d9 &
    +t*(5.87850060979944364d9 &
    +t*(4.82135265081185025d9 &
    +t*(2.83819099819062720d9 &
    +t*(1.22325945046779759d9 &
    +t*(3.76487745749965788d8 &
    +t*(7.58303164598034022d7 &
    +t*7.71973955404921694d6 &
    ))))))))/(101717.15224358499d0 &
    +s*(25406.290690742596d0 &
    +s*(11612.352340539603d0 &
    +s*(4174.0221825716411d0 &
    +s*(1281.7558184421595d0 &
    +s*(331.85386698947176d0 &
    +s*(69.634528227703809d0 &
    +s*(10.868655326785347d0 &
    +s))))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    fd=(5.2085062280477356d11 &
    +t*(2.2833424052547072d12 &
    +t*(4.6942573239060036d12 &
    +t*(5.8958834991587034d12 &
    +t*(4.9360061465395771d12 &
    +t*(2.8153080727145281d12 &
    +t*(1.06730054566227826d12 &
    +t*(2.4600629203104023d11 &
    +t*2.6466259504577480d10 &
    ))))))))/(2.42061644583385800d6 &
    +t*(153096.878206723536d0 &
    +t*(38844.6431417054012d0 &
    +t*(-16885.0152281400560d0 &
    +t*(4865.05644899569952d0 &
    +t*(-1067.44742924509046d0 &
    +t*(173.008767319757249d0 &
    +t*(-18.6099093504480426d0 &
    +t))))))))*0.999999999999999770d0	! correction to remove bias
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    fd=(7.7494086475291441d13 &
    +t*(5.0061227151520991d14 &
    +t*(1.45146838868293917d15 &
    +t*(2.46379854838342831d15 &
    +t*(2.67523099025981612d15 &
    +t*(1.90126580903196655d15 &
    +t*(8.6325468139568610d14 &
    +t*(2.28901218799706877d14 &
    +t*2.71416555283189345d13 &
    ))))))))/(8.9593891719102305d6 &
    +t*(2.23010123362691961d6 &
    +t*(-239561.734722929070d0 &
    +t*(43938.8978835860601d0 &
    +t*(-8408.3951285804076d0 &
    +t*(1448.44898008864648d0 &
    +t*(-201.658017368512975d0 &
    +t*(19.7129633131788616d0 &
    -t))))))))*0.999999999999999586d0	! correction to remove bias
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    fd=(1.22298919425291664d16 &
    +t*(8.9726625177108898d16 &
    +t*(2.89429751278928300d17 &
    +t*(5.3577901105140076d17 &
    +t*(6.2206348573306568d17 &
    +t*(4.63426239396617756d17 &
    +t*(2.16065376561259960d17 &
    +t*(5.7540345485220279d16 &
    +t*6.6835199866465357d15 &
    ))))))))/(1.31871866498468464d7 &
    +t*(2.83885323472743888d6 &
    +t*(-275298.111415499841d0 &
    +t*(47078.0598194149218d0 &
    +t*(-8612.3789225919783d0 &
    +t*(1446.29041011108831d0 &
    +t*(-199.315359842233309d0 &
    +t*(19.5132779151984627d0 &
    -t))))))))*0.999999999999999580d0	! correction to remove bias
else
    w=1.d0/(x*x)
    t=1600.d0*w
    fd=x*x*x*x*x*x*x*sqrt(x)*factor*(1.d0 &
    +w*(80.1905357588510589d0 &
    +t*(1.42831787292433216d0 &
    +t*(0.00812881145350198890d0 &
    +t*(3.85165680520411590d-6 &
    +t*1.83598087001386478d-9)))))
endif
!write(*,"(a20,1p2e15.7)") "(fd13h) t,fd=",t,fd
fd13h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd14h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=7
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,s,fd,y2
!
!write(*,"(a20,1pe15.7)") "(fd14h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    s=1.d0-t
    fd=ex*(5040.d0 &
    -ex*(19.5849162780581217d0 &
    +s*(0.101236264600248646d0 &
    +s*(0.00131827373096852460d0 &
    +s*(0.0000283211236200358235d0 &
    +s*(8.27819847443616991d-7 &
    +s*3.46664379401620231d-8))))))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(8.05214238081846197d6 &
    +t*(4.23358827341165564d6 &
    +t*(262627.912342619479d0 &
    +t*(-323968.614686001584d0 &
    +t*(-98175.8116505823446d0 &
    +t*(-7562.84276896647246d0 &
    +t*(1252.57372771721279d0 &
    +t*240.247242901336449d0 &
    )))))))/(2413.8835947192718d0 &
    +s*(4139.8252477104615d0 &
    +s*(3204.5076382872745d0 &
    +s*(1486.4720567970787d0 &
    +s*(456.50435752749115d0 &
    +s*(95.80622002968487d0 &
    +s*(13.317952574698873d0 &
    +s)))))))
endif
if(y.gt.0.d0) then
    y2=y*y
    fd=-fd+10042.0286586746908d0 &
    +y2*(4967.17750013907292d0 &
    +y2*(397.753788388843285d0 &
    +y2*(11.5145384679375851d0 &
    +y2*0.125d0)))
endif
!write(*,"(a20,1p2e15.7)") "(fd14h) y,fd=",y,fd
fd14h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd15h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=15/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/17.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd15h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(14034.4072934834126d0 &
    -ex*(171704.258234536602d0 &
    +t*(6049.5894628442878d0 &
    +t*18.1763510883857458d0 &
    ))/(4429.36992795387362d0 &
    +t*(175.155804618692908d0 &
    +t)))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(2.06790092766575034d7 &
    +t*(8.26618273142399596d6 &
    +t*(-2.08115753677480168d6 &
    +t*(-2.17301177649999531d6 &
    +t*(-641552.116750411401d0 &
    +t*(-91807.9577537648588d0 &
    +t*(-5144.12042361256270d0 &
    +t*130.857943513312398d0 &
    )))))))/(1711.2966803109211d0 &
    +s*(3804.4823677048040d0 &
    +s*(3232.4717528599742d0 &
    +s*(1550.7278916841484d0 &
    +s*(478.50436506151174d0 &
    +s*(99.39616883636118d0 &
    +s*(13.585756243903013d0 &
    +s)))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    s=1.d0-t
    fd=(1.14403769376860621d8 &
    +t*(8.4972567501965974d7 &
    +t*(3.81753609315202610d7 &
    +t*(1.65881819811148351d7 &
    +t*(6.1705051839981204d6 &
    +t*(1.58079386467916017d6 &
    +t*(246507.752966115175d0 &
    +t*18602.2895319461278d0 &
    )))))))/(2572.2942565752562d0 &
    +s*(2659.5730458609857d0 &
    +s*(1732.1057976186518d0 &
    +s*(829.4655494427400d0 &
    +s*(293.34263153462293d0 &
    +s*(73.73539965269387d0 &
    +s*(12.058936149293461d0 &
    +s)))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    s=1.d0-t
    fd=(6.26769725573910255d9 &
    +t*(1.32711699403202671d10 &
    +t*(1.45560142717062095d10 &
    +t*(1.07860795716164403d10 &
    +t*(5.86969120046064229d9 &
    +t*(2.38796511791514625d9 &
    +t*(7.07324240427174132d8 &
    +t*(1.39791587840460999d8 &
    +t*1.42694865767192645d7 &
    ))))))))/(30361.411562018280d0 &
    +s*(17707.813587539272d0 &
    +s*(8523.204548822773d0 &
    +s*(3390.1996842891790d0 &
    +s*(1127.1366883800605d0 &
    +s*(309.69425190845478d0 &
    +s*(67.751604120864452d0 &
    +s*(10.834010773474702d0 &
    +s))))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    s=1.d0-t
    fd=(5.69341524462926732d11 &
    +t*(2.36489781301030086d12 &
    +t*(4.67723648136547218d12 &
    +t*(5.73677938363737470d12 &
    +t*(4.75966348844286492d12 &
    +t*(2.72977815388509051d12 &
    +t*(1.05619383963810221d12 &
    +t*(2.52484808268258635d11 &
    +t*2.87129101561538399d10 &
    ))))))))/(233403.35536095456d0 &
    +s*(60565.447915936497d0 &
    +s*(18770.529358157892d0 &
    +s*(5494.0347508271098d0 &
    +s*(1458.2356321178129d0 &
    +s*(341.00432855822186d0 &
    +s*(67.360611688552025d0 &
    +s*(10.355894270482783d0 &
    +s))))))))*0.999999999999999687d0	! correction to remove bias
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    s=1.d0-t
    fd=(6.79824678667324436d13 &
    +t*(4.47889271605157194d14 &
    +t*(1.33539302658610566d15 &
    +t*(2.35081455707368313d15 &
    +t*(2.67105097393954020d15 &
    +t*(2.00614187037046236d15 &
    +t*(9.73640078908216993d14 &
    +t*(2.79804648854039015d14 &
    +t*3.66107518187916541d13 &
    ))))))))/(585167.34458490542d0 &
    +s*(99101.77544747319d0 &
    +s*(23659.896891961921d0 &
    +s*(5847.0608231871316d0 &
    +s*(1394.1981365390388d0 &
    +s*(307.36846956466131d0 &
    +s*(59.734384735038341d0 &
    +s*(9.438593558054937d0 &
    +s))))))))*0.999999999999999620d0	! correction to remove bias
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    s=1.d0-t
    fd=(1.99221037305471970d16 &
    +t*(1.55301826780028291d17 &
    +t*(5.3562790972794281d17 &
    +t*(1.06787028331863552d18 &
    +t*(1.34677794609411928d18 &
    +t*(1.10118289362700600d18 &
    +t*(5.7076444557408842d17 &
    +t*(1.71793041083229002d17 &
    +t*2.30600690809359773d16 &
    ))))))))/(959281.46439522639d0 &
    +s*(146009.54484661068d0 &
    +s*(31788.065669901676d0 &
    +s*(7257.6052635887384d0 &
    +s*(1618.0190118885021d0 &
    +s*(337.33224463190506d0 &
    +s*(62.698857808593658d0 &
    +s*(9.5893450572160768d0 &
    +s))))))))*0.999999999999999649d0	! correction to remove bias
else
    w=1.d0/(x*x)
    t=1600.d0*w
    fd=x*x*x*x*x*x*x*x*sqrt(x)*factor*(1.d0 &
    +w*(104.864546761574387d0 &
    +t*(2.69793375997248149d0 &
    +t*(0.0276379588953803004d0 &
    +t*(0.0000654786006189425666d0 &
    -t*1.03209499898071117d-8)))))
endif
!write(*,"(a20,1p2e15.7)") "(fd15h) t,fd=",t,fd
fd15h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd16h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=8
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,s,fd,y2
!
!write(*,"(a20,1pe15.7)") "(fd16h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(40320.d0 &
    -ex*(438381.668209835602d0 &
    +t*(14167.5566380822417d0 &
    +t*39.1239843501671193d0 &
    ))/(5566.7513423478282d0 &
    +t*(199.502568378336498d0 &
    +t)))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(2.14460848478648315d8 &
    +t*(1.66513461191585055d8 &
    +t*(5.17491656951190442d7 &
    +t*(6.60393533940475286d6 &
    +t*(-301119.943009454383d0 &
    +t*(-202842.867158865294d0 &
    -t*20669.3789063593406d0 &
    ))))))/(10903.784789014007d0 &
    +s*(14623.305826998018d0 &
    +s*(9179.893497142397d0 &
    +s*(3522.2940212081378d0 &
    +s*(905.1498999681050d0 &
    +s*(159.07578981158691d0 &
    +s*(18.015010064955669d0 &
    +s)))))))
endif
if(y.gt.0.d0) then
    y2=y*y
    fd=fd+y*(80336.2292693975266d0 &
    +y2*(13245.8066670375278d0 &
    +y2*(636.406061422149257d0 &
    +y2*(13.1594725347858115d0 &
    +y2*0.111111111111111111d0))))
endif
!write(*,"(a20,1p2e15.7)") "(fd16h) y,fd=",y,fd
fd16h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd17h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=17/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/19.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd17h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    s=1.d0-t
    fd=ex*(119292.461994609007d0 &
    -ex*(164.281538211493590d0 &
    +s*(0.465418918711191171d0 &
    +s*(0.00397356161994318188d0 &
    +s*(0.0000615968461595874662d0 &
    +s*1.50191715167058393d-6)))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(1.43233835112073902d8 &
    +t*(1.39644553704144873d8 &
    +t*(6.13766014358067700d7 &
    +t*(1.56246791599180488d7 &
    +t*(2.43295219093432312d6 &
    +t*(218176.915876999296d0 &
    +t*8461.46979813306275d0 &
    ))))))/(3043.1981753546660d0 &
    +s*(3390.9687009373892d0 &
    +s*(1761.0786885568746d0 &
    +s*(550.68360665579008d0 &
    +s*(112.08528557398394d0 &
    +s*(14.640373761706403d0 &
    +s))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    s=1.d0-t
    fd=(1.21473360069923149d9 &
    +t*(1.01176690309922349d9 &
    +t*(4.95296056344158641d8 &
    +t*(2.06442737728860908d8 &
    +t*(7.0192333337028334d7 &
    +t*(1.68572449951275985d7 &
    +t*(2.53144559406161612d6 &
    +t*186718.775025914870d0 &
    )))))))/(3455.142635040908d0 &
    +s*(3415.073282258220d0 &
    +s*(2023.9784836523565d0 &
    +s*(903.6597935624893d0 &
    +s*(308.93682588995655d0 &
    +s*(76.49008105409234d0 &
    +s*(12.338920623517048d0 &
    +s)))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    s=1.d0-t
    fd=(2.67440099736518804d10 &
    +t*(4.91370104189578919d10 &
    +t*(4.73732500104890339d10 &
    +t*(3.16238224173646457d10 &
    +t*(1.59135818804699047d10 &
    +t*(6.12061461496369111d9 &
    +t*(1.74615616554871192d9 &
    +t*(3.38446352294444219d8 &
    +t*3.46053509281020592d7 &
    ))))))))/(11121.047017722056d0 &
    +s*(9829.804681621007d0 &
    +s*(5748.503365098324d0 &
    +s*(2606.1477102230902d0 &
    +s*(953.9183383047316d0 &
    +s*(281.64461270648236d0 &
    +s*(64.892044676666959d0 &
    +s*(10.720667663673541d0 &
    +s))))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    s=1.d0-t
    fd=(1.21689593180741239d12 &
    +t*(4.79662223812667505d12 &
    +t*(9.12838410654446290d12 &
    +t*(1.09342887556403655d13 &
    +t*(8.99572172197880880d12 &
    +t*(5.19592600765981107d12 &
    +t*(2.05751238857523311d12 &
    +t*(5.12310527144160372d11 &
    +t*6.19902241659841813d10 &
    ))))))))/(39255.392055313674d0 &
    +s*(22199.471524635690d0 &
    +s*(9415.210433779872d0 &
    +s*(3354.6940213500329d0 &
    +s*(1027.9861625160893d0 &
    +s*(269.15608801102034d0 &
    +s*(58.326091875891505d0 &
    +s*(9.690839898241749d0 &
    +s))))))))*0.999999999999999611d0
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    s=1.d0-t
    fd=(8.36136698862497927d13 &
    +t*(5.63471803594835697d14 &
    +t*(1.73344511753855139d15 &
    +t*(3.17713365223451161d15 &
    +t*(3.79488820022826435d15 &
    +t*(3.02846836744419910d15 &
    +t*(1.58129984938571763d15 &
    +t*(4.96500350365744464d14 &
    +t*7.24731492178964020d13 &
    ))))))))/(43633.188864890011d0 &
    +s*(21099.956516665469d0 &
    +s*(8040.798889514530d0 &
    +s*(2678.4837772160305d0 &
    +s*(794.3651890403060d0 &
    +s*(208.01106751201808d0 &
    +s*(46.655407061024955d0 &
    +s*(8.3567416627775934d0 &
    +s))))))))*0.999999999999999524d0
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    s=1.d0-t
    fd=(2.85846778936881832d16 &
    +t*(2.36061789976170730d17 &
    +t*(8.67582358035872224d17 &
    +t*(1.85563671344211925d18 &
    +t*(2.53068942615802160d18 &
    +t*(2.25901769990104319d18 &
    +t*(1.29359182644194722d18 &
    +t*(4.36848670280351452d17 &
    +t*6.72189933519645798d16 &
    ))))))))/(50522.192743560920d0 &
    +s*(23134.306793549278d0 &
    +s*(8428.128974587599d0 &
    +s*(2708.9200442813342d0 &
    +s*(782.3736655720094d0 &
    +s*(201.43688280680891d0 &
    +s*(44.899007032888684d0 &
    +s*(8.0956175513436493d0 &
    +s))))))))*0.999999999999999478d0
else
    w=1.d0/(x*x)
    t=1600.d0*w
    fd=x*x*x*x*x*x*x*x*x*sqrt(x)*factor*(1.d0 &
    +w*(132.828425897994463d0 &
    +t*(4.66006740357572083d0 &
    +t*(0.0750173171258174240d0 &
    +t*(0.000414697262104018271d0 &
    +t*1.95167395079296505d-7)))))
endif
!write(*,"(a20,1p2e15.7)") "(fd17h) t,fd=",t,fd
fd17h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd18h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=9
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,s,fd,y2
!
!write(*,"(a20,1pe15.7)") "(fd18h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    fd=ex*(362880.d0 &
    -ex*(3.11619728815890900d6 &
    +t*(84459.263086739990d0 &
    +t*194.729051789880448d0 &
    ))/(8793.5020477151129d0 &
    +t*(258.970670923338654d0 &
    +t)))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(5.34284039564434465d8 &
    +t*(5.33773366525970850d8 &
    +t*(2.42092992302655038d8 &
    +t*(6.42814979745526752d7 &
    +t*(1.06272445991848042d7 &
    +t*(1.04623859092301352d6 &
    +t*48187.9003961411735d0 &
    ))))))/(3823.5401164205447d0 &
    +s*(4167.5033491759870d0 &
    +s*(2107.9832247353650d0 &
    +s*(639.46059745678970d0 &
    +s*(125.53444672868253d0 &
    +s*(15.642276104269218d0 &
    +s))))))
endif
if(y.gt.0.d0) then
    y2=y*y
    fd=-fd+725062.913034521571d0 &
    +y2*(361513.031712288870d0 &
    +y2*(29803.0650008344375d0 &
    +y2*(954.609092133223885d0 &
    +y2*(14.8044066016340379d0 &
    +y2*0.1d0))))
endif
!write(*,"(a20,1p2e15.7)") "(fd18h) y,fd=",y,fd
fd18h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd19h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=19/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/21.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd19h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    s=1.d0-t
    fd=ex*(1.13327838894878557d6 &
    -ex*(781.077395263087427d0 &
    +s*(1.48017751152694790d0 &
    +s*(0.00952643972760576450d0 &
    +s*(0.000118887473010527875d0 &
    +s*2.40377778104657177d-6)))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(1.98943681340803829d9 &
    +t*(2.02504701261695886d9 &
    +t*(9.40112997876683365d8 &
    +t*(2.57205491041485272d8 &
    +t*(4.42589158347859287d7 &
    +t*(4.61202797242277802d6 &
    +t*232199.034580285755d0 &
    ))))))/(4645.3648503361484d0 &
    +s*(4976.8851201391641d0 &
    +s*(2464.6700288281232d0 &
    +s*(729.19053696862632d0 &
    +s*(138.78343942961328d0 &
    +s*(16.586263194372411d0 &
    +s))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    s=1.d0-t
    fd=(2.0904667585793652d10 &
    +t*(1.9423314417935927d10 &
    +t*(9.4466078397530459d9 &
    +t*(3.3380904172376415d9 &
    +t*(9.2799283366292994d8 &
    +t*(1.9279930343506925d8 &
    +t*(2.6584305328117208d7 &
    +t*1.8721659680393183d6 &
    )))))))/(6510.453486834102d0 &
    +s*(6583.316420184059d0 &
    +s*(3538.873468101783d0 &
    +s*(1337.9108450581962d0 &
    +s*(388.7430265102922d0 &
    +s*(85.60461680432659d0 &
    +s*(12.861326654789701d0 &
    +s)))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    s=1.d0-t
    fd=(1.46863114809614622d11 &
    +t*(2.34916998394298104d11 &
    +t*(1.98120796400802870d11 &
    +t*(1.18887551091630769d11 &
    +t*(5.55928633602927965d10 &
    +t*(2.03490255735994097d10 &
    +t*(5.60239395267172792d9 &
    +t*(1.05980329750305191d9 &
    +t*1.07474670751003016d8 &
    ))))))))/(4912.718774356158d0 &
    +s*(5658.595947000664d0 &
    +s*(3911.882548894156d0 &
    +s*(2003.2886387850459d0 &
    +s*(805.2407096838510d0 &
    +s*(255.52601531171211d0 &
    +s*(62.070351460815245d0 &
    +s*(10.603837914443282d0 &
    +s))))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    s=1.d0-t
    fd=(4.12413857105017004d12 &
    +t*(1.53777683873017139d13 &
    +t*(2.80178767823520436d13 &
    +t*(3.25874492182871964d13 &
    +t*(2.64367160554342534d13 &
    +t*(1.53023378021225309d13 &
    +t*(6.17571294140739606d12 &
    +t*(1.59629111145863694d12 &
    +t*2.05057950138673222d11 &
    ))))))))/(9870.066781117703d0 &
    +s*(8442.670237353289d0 &
    +s*(4641.969549915204d0 &
    +s*(1995.0747902033748d0 &
    +s*(707.3070308400987d0 &
    +s*(208.47240484486439d0 &
    +s*(49.869236685972587d0 &
    +s*(9.0113347734666514d0 &
    +s))))))))*0.999999999999999624d0	! correction to remove bias
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    s=1.d0-t
    fd=(4.83390017425457148d15 &
    +t*(3.49375860214157535d16 &
    +t*(1.17174577652359819d17 &
    +t*(2.38956570804562992d17 &
    +t*(3.26219845844040472d17 &
    +t*(3.09099843298003147d17 &
    +t*(2.03426110846076766d17 &
    +t*(8.98530419590583886d16 &
    +t*(2.42666818156258567d16 &
    +t*3.07539397786247273d15 &
    )))))))))/(207639.00554015383d0 &
    +s*(101594.67784669785d0 &
    +s*(39318.472360728631d0 &
    +s*(13385.302295791795d0 &
    +s*(4099.2485074048775d0 &
    +s*(1128.0104124272086d0 &
    +s*(274.41287517597127d0 &
    +s*(56.961169613859215d0 &
    +s*(9.354997067031336d0 &
    +s)))))))))*0.999999999999999475d0	! correction to remove bias
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    s=1.d0-t
    fd=(2.81147322924282096d18 &
    +t*(2.55872925207166842d19 &
    +t*(1.05126940075992923d20 &
    +t*(2.56086191531691937d20 &
    +t*(4.08006931722902169d20 &
    +t*(4.41541112498720876d20 &
    +t*(3.25219940350974164d20 &
    +t*(1.57672402355387744d20 &
    +t*(4.58547965512713670d19 &
    +t*6.13678661161381742d18 &
    )))))))))/(253648.48647211156d0 &
    +s*(116226.83914995012d0 &
    +s*(42576.299330336023d0 &
    +s*(13857.868987762036d0 &
    +s*(4097.2448088308021d0 &
    +s*(1099.1091735270184d0 &
    +s*(263.31471897218326d0 &
    +s*(54.436947056837050d0 &
    +s*(9.0281927555123252d0 &
    +s)))))))))*0.999999999999999480d0	! correction to remove bias
else
    w=1.d0/(x*x)
    t=1600.d0*w
    fd=x*x*x*x*x*x*x*x*x*x*sqrt(x)*factor*(1.d0 &
    +w*(164.082173168110588d0 &
    +t*(7.52780119040598426d0 &
    +t*(0.175040406459615661d0 &
    +t*(0.00174172934550196671d0 &
    +t*(4.09368331231164635d-6 &
    -t*6.44677689509362889d-10))))))
endif
!write(*,"(a20,1p2e15.7)") "(fd19h) t,fd=",t,fd
fd19h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd20h(y)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=10
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 y,x,ex,t,s,fd,y2
!
!write(*,"(a20,1pe15.7)") "(fd20h) y=",y
x=-abs(y)
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    s=1.d0-t
    fd=ex*(3628800.d0 &
    -ex*(1769.11836499240874d0 &
    +s*(2.74114748349296929d0 &
    +s*(0.0153129509121134183d0 &
    +s*(0.000171415624076050339d0 &
    +s*3.15741345799313491d-6)))))
elseif(x.le.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(7.23923596166699396d9 &
    +t*(7.46011115766110324d9 &
    +t*(3.51560076972821513d9 &
    +t*(9.79963332403905796d8 &
    +t*(1.72722667850513422d8 &
    +t*(1.85872779338803730d7 &
    +t*980119.027043354796d0 &
    ))))))/(5345.1738990257193d0 &
    +s*(5661.6190861926977d0 &
    +s*(2763.5410988949565d0 &
    +s*(803.48749805350521d0 &
    +s*(149.54574471403198d0 &
    +s*(17.324182954806273d0 &
    +s))))))
endif
if(y.gt.0.d0) then
    y2=y*y
    fd=fd+y*(7.25062913034521571d6 &
    +y2*(1.20504343904096290d6 &
    +y2*(59606.1300016688751d0 &
    +y2*(1363.72727447603412d0 &
    +y2*(16.4493406684822644d0 &
    +y2*0.0909090909090909091d0)))))
endif
!write(*,"(a20,1p2e15.7)") "(fd20h) y,fd=",y,fd
fd20h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real*8 function fd21h(x)
!
! double precision rational minimax approximation of Fermi-Dirac integral of order k=21/2
!
! Reference: Fukushima, T. (2014, submitted to App. Math. Comp.)
!
! Author: Fukushima, T. <Toshio.Fukushima@nao.ac.jp>
!
real*8 x,ex,t,w,s,fd,factor
parameter (factor=2.d0/23.d0)    ! = 1/(k+1)
!
!write(*,"(a20,1pe15.7)") "(fd21h) x=",x
if(x.lt.-2.d0) then
    ex=exp(x)
    t=ex*7.38905609893065023d0
    s=1.d0-t
    fd=ex*(1.18994230839622485d7 &
    -ex*(4103.25503466860405d0 &
    +s*(5.19741906649655513d0 &
    +s*(0.0251967194122553158d0 &
    +s*(0.000252954565705636841d0 &
    +s*4.24573915094571073d-6)))))
elseif(x.lt.0.d0) then
    s=-0.5d0*x
    t=1.d0-s
    fd=(2.57000236613575175d10 &
    +t*(2.66775768017824432d10 &
    +t*(1.26823522594150685d10 &
    +t*(3.57329482988205525d9 &
    +t*(6.38366360660921950d8 &
    +t*(6.99165387034866263d7 &
    +t*3.77747237079641640d6 &
    ))))))/(5829.6137460125218d0 &
    +s*(6133.1261260721688d0 &
    +s*(2967.5333196500928d0 &
    +s*(853.66205862182849d0 &
    +s*(156.68183527916256d0 &
    +s*(17.793529538010685d0 &
    +s))))))
elseif(x.lt.2.d0) then
    t=0.5d0*x
    s=1.d0-t
    fd=(3.9256463321926022d11 &
    +t*(3.5272008321450459d11 &
    +t*(1.4960338986194452d11 &
    +t*(4.0840964995951679d10 &
    +t*(8.3525752792093004d9 &
    +t*(1.3652900099806473d9 &
    +t*(1.6653931710965140d8 &
    +t*1.13722223824443501d7 &
    )))))))/(10780.778916243366d0 &
    +s*(12211.902429240048d0 &
    +s*(6810.538241268853d0 &
    +s*(2444.9651837399508d0 &
    +s*(622.4594160475152d0 &
    +s*(115.21708456233468d0 &
    +s*(14.653186311439082d0 &
    +s)))))))
elseif(x.lt.5.d0) then
    t=0.3333333333333333333d0*(x-2.d0)
    s=1.d0-t
    fd=(1.02038438555518024d12 &
    +t*(1.44086845004445038d12 &
    +t*(1.06951376053758491d12 &
    +t*(5.83349119667798955d11 &
    +t*(2.59184513055987834d11 &
    +t*(9.22344989320630316d10 &
    +t*(2.46767885435893916d10 &
    +t*(4.50298910194390990d9 &
    +t*4.40533249962065063d8 &
    ))))))))/(2627.6481435132529d0 &
    +s*(3598.485403241250d0 &
    +s*(2816.3184982451496d0 &
    +s*(1591.2052134075520d0 &
    +s*(693.3250333689612d0 &
    +s*(234.74251456309597d0 &
    +s*(59.835610870372227d0 &
    +s*(10.529060212037859d0 &
    +s))))))))
elseif(x.lt.10.d0) then
    t=0.2d0*x-1.d0
    s=1.d0-t
    fd=(1.92539440493876721d13 &
    +t*(6.76427905453399646d13 &
    +t*(1.17204818714740525d14 &
    +t*(1.31312575057612394d14 &
    +t*(1.04191636638567130d14 &
    +t*(5.99784170385635521d13 &
    +t*(2.45015915838649510d13 &
    +t*(6.53409875052899427d12 &
    +t*8.86443190814042498d11 &
    ))))))))/(3208.0283085434934d0 &
    +s*(3638.3347710577402d0 &
    +s*(2454.8591095140763d0 &
    +s*(1237.5705806010665d0 &
    +s*(499.5532233797705d0 &
    +s*(164.10632643205828d0 &
    +s*(43.045792336177945d0 &
    +s*(8.4174363747925276d0 &
    +s))))))))*0.999999999999999649d0	! correction to remove bias
elseif(x.lt.20.d0) then
    t=0.1d0*x-1.d0
    s=1.d0-t
    fd=(1.35032839928658828d16 &
    +t*(9.84341732357224541d16 &
    +t*(3.35732954478237964d17 &
    +t*(7.02136022561204639d17 &
    +t*(9.91445873700748081d17 &
    +t*(9.80394164297009879d17 &
    +t*(6.79928848238746956d17 &
    +t*(3.19982319837818266d17 &
    +t*(9.33078432378309346d16 &
    +t*1.29960135147133318d16 &
    )))))))))/(32616.320007579440d0 &
    +s*(25996.868562405952d0 &
    +s*(13806.569048984017d0 &
    +s*(5938.718527590630d0 &
    +s*(2189.8437434853080d0 &
    +s*(703.2974320942947d0 &
    +s*(195.48717008127781d0 &
    +s*(45.691859712028360d0 &
    +s*(8.3721801701244348d0 &
    +s)))))))))*0.999999999999999529d0	! correction to remove bias
elseif(x.lt.40.d0) then
    t=0.05d0*x-1.d0
    s=1.d0-t
    fd=(1.01458237235662575d19 &
    +t*(9.67883827639700444d19 &
    +t*(4.18744967361252695d20 &
    +t*(1.07961610867931522d21 &
    +t*(1.83106271883557796d21 &
    +t*(2.12354501642149972d21 &
    +t*(1.68956127221959607d21 &
    +t*(8.93514699913180709d20 &
    +t*(2.87036545075886579d20 &
    +t*4.31795348604572483d19 &
    )))))))))/(32461.378904270757d0 &
    +s*(24848.502722479329d0 &
    +s*(12785.219096902698d0 &
    +s*(5373.990729362032d0 &
    +s*(1953.3664643688677d0 &
    +s*(624.0842567143388d0 &
    +s*(174.30734991710066d0 &
    +s*(41.427358351772922d0 &
    +s*(7.8380857354891532d0 &
    +s)))))))))*0.999999999999999445d0	! correction to remove bias
else
    w=1.d0/(x*x)
    t=1600.d0*w
    fd=x*x*x*x*x*x*x*x*x*x*x*sqrt(x)*factor*(1.d0 &
    +w*(198.625788571923339d0 &
    +t*(11.5426284919561671d0 &
    +t*(0.365993577138676996d0 &
    +t*(0.00572282501335319477d0 &
    +t*(0.0000313848422097767888d0 &
    +t*1.47431097320479878d-8))))))
endif
!write(*,"(a20,1p2e15.7)") "(fd21h) t,fd=",t,fd
fd21h=fd
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! -4.5 -3.0    7.240015921650838E-03
! -4.5 -1.0   -2.442145618002432E-02
! -4.5  1.0    1.647586588274429E-02
! -4.5  3.5    1.029804636049462E-03
! -4.5  7.5   -4.060975578007697E-04
! -4.5 15.0   -2.508232176525567E-05
! -4.5 30.0   -1.990373968491516E-06
! -4.5 60.0   -1.720097159754510E-07
!
! -3.5 -3.0   -3.545676388861523E-02
! -3.5 -1.0   -2.657136470551858E-02
! -3.5  1.0    9.517508954351293E-02
! -3.5  3.5   -1.294450181808404E-02
! -3.5  7.5   -3.613853399968076E-03
! -3.5 15.0   -4.935123501719669E-04
! -3.5 30.0   -8.248542579808003E-05
! -3.5 60.0   -1.440219364000877E-05
!
! -2.5 -3.0    1.024984312340880E-01
! -2.5 -1.0    3.492797177355019E-01
! -2.5  1.0    2.502420253443888E-02
! -2.5  3.5   -1.233961959341110E-01
! -2.5  7.5   -3.760442897358426E-02
! -2.5 15.0   -1.182183957963569E-02
! -2.5 30.0   -4.085597149537218E-03
! -2.5 60.0   -1.436908659724343E-03
!
! -1.5 -3.0   -1.647804107730068E-01
! -1.5 -1.0   -8.394844152932343E-01
! -1.5  1.0   -1.580053773447138E+00
! -1.5  3.5   -1.176280028545422E+00
! -1.5  7.5   -7.505784586840675E-01
! -1.5 15.0   -5.193758746338315E-01
! -1.5 30.0   -3.656546825930124E-01
! -1.5 60.0   -2.582876225442183E-01
!
! -0.5 -3.0    8.525970123326887E-02
! -0.5 -1.0    5.211503831079912E-01
! -0.5  1.0    1.820411357146962E+00
! -0.5  3.5    3.591334824651103E+00
! -0.5  7.5    5.432873628864860E+00
! -0.5 15.0    7.731513057173728E+00
! -0.5 30.0    1.094942130440661E+01
! -0.5 60.0    1.549016158518217E+01
!
!  0.0 -3.0    4.858735157374205E-02
!  0.0 -1.0    3.132616875182229E-01
!  0.0  1.0    1.313261687518223E+00
!  0.0  3.5    3.529750418272620E+00
!  0.0  7.5    7.500552931475361E+00
!  0.0 15.0    1.500000030590227E+01
!  0.0 30.0    3.000000000000009E+01
!  0.0 60.0    6.000000000000000E+01
!
!  0.5 -3.0    4.336636755041557E-02
!  0.5 -1.0    2.905008961699176E-01
!  0.5  1.0    1.396375280666564E+00
!  0.5  3.5    4.837065897622567E+00
!  0.5  7.5    1.399909743357599E+01
!  0.5 15.0    3.894304660093270E+01
!  0.5 30.0    1.096948183372665E+02
!  0.5 60.0    3.099448732700438E+02
!
!  1.0 -3.0    4.918072033882423E-02
!  1.0 -1.0    3.386479964034522E-01
!  1.0  1.0    1.806286070444774E+00
!  1.0  3.5    7.739961645298564E+00
!  1.0  7.5    2.976938105893487E+01
!  1.0 15.0    1.141449337609459E+02
!  1.0 30.0    4.516449340668481E+02
!  1.0 60.0    1.801644934066848E+03
!
!  1.5 -3.0    6.561173880637544E-02
!  1.5 -1.0    4.608488062901017E-01
!  1.5  1.0    2.661682624732004E+00
!  1.5  3.5    1.365420168610915E+01
!  1.5  7.5    6.833812856132876E+01
!  1.5 15.0    3.581122477085265E+02
!  1.5 30.0    1.985311377746038E+03
!  1.5 60.0    1.117330291388515E+04
!
!  2.0 -3.0    9.896340290959225E-02
!  2.0 -1.0    7.051297585956156E-01
!  2.0  1.0    4.328331225625401E+00
!  2.0  3.5    2.586637394510408E+01
!  2.0  7.5    1.653001170950006E+02
!  2.0 15.0    1.174348022617251E+03
!  2.0 30.0    9.098696044010894E+03
!  2.0 60.0    7.219739208802179E+04
!
!  2.5 -3.0    1.647403937322051E-01
!  2.5 -1.0    1.185968175443467E+00
!  2.5  1.0    7.626535355005596E+00
!  2.5  3.5    5.186981146923186E+01
!  2.5  7.5    4.158852838077063E+02
!  2.5 15.0    3.974487881941656E+03
!  2.5 30.0    4.292925758509993E+04
!  2.5 60.0    4.799485008429281E+05
!
!  3.0 -3.0    2.978018784709205E-01
!  3.0 -1.0    2.159839661017986E+00
!  3.0  1.0    1.438935649349364E+01
!  3.0  3.5    1.091505015453568E+02
!  3.0  7.5    1.079959324343085E+03
!  3.0 15.0    1.377794488724111E+04
!  3.0 30.0    2.069526863744442E+05
!  3.0 60.0    3.257776652315915E+06
!
!  3.5 -3.0    5.778455375087482E-01
!  3.5 -1.0    4.213264071926359E+00
!  3.5  1.0    2.883131841599375E+01
!  3.5  3.5    2.396666058159899E+02
!  3.5  7.5    2.880213529254126E+03
!  3.5 15.0    4.868423690758540E+04
!  3.5 30.0    1.014417201742514E+06
!  3.5 60.0    2.246912083856067E+07
!
!  4.0 -3.0    1.193042623617263E+00
!  4.0 -1.0    8.732138288905944E+00
!  4.0  1.0    6.096945037216665E+01
!  4.0  3.5    5.469755138110205E+02
!  4.0  7.5    7.862865080220864E+03
!  4.0 15.0    1.747634735470307E+05
!  4.0 30.0    5.039016606494085E+06
!  4.0 60.0    1.569439504883058E+08
!
!  4.5 -3.0    2.603141667351258E+00
!  4.5 -1.0    1.910506806522883E+01
!  4.5  1.0    1.354192084601931E+02
!  4.5  3.5    1.293868145709697E+03
!  4.5  7.5    2.192168139830375E+04
!  4.5 15.0    6.358400491901155E+05
!  4.5 30.0    2.530631914465860E+07
!  4.5 60.0    1.107558362743446E+09
!
!  5.0 -3.0    5.969820680488160E+00
!  5.0 -1.0    4.389948426765305E+01
!  5.0  1.0    3.146680541843087E+02
!  5.0  3.5    3.165640736730491E+03
!  5.0  7.5    6.231539440840174E+04
!  5.0 15.0    2.340617854292586E+06
!  5.0 30.0    1.282644990485829E+08
!  5.0 60.0    7.883001082246370E+09
!
!  5.5 -3.0    1.432510773316834E+01
!  5.5 -1.0    1.054873701421355E+02
!  5.5  1.0    7.624071301487685E+02
!  5.5  3.5    7.997774066508814E+03
!  5.5  7.5    1.804014784983563E+05
!  5.5 15.0    8.706355279385276E+06
!  5.5 30.0    6.552423229390368E+08
!  5.5 60.0    5.651215902520975E+10
!
!  6.0 -3.0    3.583278660538172E+01
!  6.0 -1.0    2.641275790125708E+02
!  6.0  1.0    1.920621491104163E+03
!  6.0  3.5    2.083671641045758E+04
!  6.0  7.5    5.314330295476191E+05
!  6.0 15.0    3.269159748061942E+07
!  6.0 30.0    3.370296449774471E+09
!  6.0 60.0    4.076323551443539E+11
!
!  6.5 -3.0    9.313870302507033E+01
!  6.5 -1.0    6.870205993243595E+02
!  6.5  1.0    5.017993045210247E+03
!  6.5  3.5    5.591696834104130E+04
!  6.5  7.5    1.592092071043047E+06
!  6.5 15.0    1.238219532435030E+08
!  6.5 30.0    1.744018597946511E+10
!  6.5 60.0    2.956078846984141E+12
!
!  7.0 -3.0    2.508781184723398E+02
!  7.0 -1.0    1.851484886980967E+03
!  7.0  1.0    1.356711459868958E+04
!  7.0  3.5    1.544073285162487E+05
!  7.0  7.5    4.848712724195143E+06
!  7.0 15.0    4.727830603631850E+08
!  7.0 30.0    9.073325961350024E+10
!  7.0 60.0    2.153759508773864E+13
!
!  7.5 -3.0    6.986360585012507E+02
!  7.5 -1.0    5.157783274159222E+03
!  7.5  1.0    3.788356128997377E+04
!  7.5  3.5    4.383205315095304E+05
!  7.5  7.5    1.500759095117367E+07
!  7.5 15.0    1.818973214398967E+09
!  7.5 30.0    4.743320667852189E+11
!  7.5 60.0    1.575715160529580E+14
!
!  8.0 -3.0    2.007219646720646E+03
!  8.0 -1.0    1.482234071460382E+04
!  8.0  1.0    1.090540532961069E+05
!  8.0  3.5    1.277985060309017E+06
!  8.0  7.5    4.720135119144394E+07
!  8.0 15.0    7.049084121468240E+09
!  8.0 30.0    2.490622378495212E+12
!  8.0 60.0    1.157079836302754E+15
!
!  8.5 -3.0    5.938814014650166E+03
!  8.5 -1.0    4.386311784955336E+04
!  8.5  1.0    3.231138561347001E+05
!  8.5  3.5    3.823781896168859E+06
!  8.5  7.5    1.508420732678749E+08
!  8.5 15.0    2.750798130249417E+10
!  8.5 30.0    1.313056860664467E+13
!  8.5 60.0    8.524941900028644E+15
!
!  9.0 -3.0    1.806585371781600E+04
!  9.0 -1.0    1.334484320310655E+05
!  9.0  1.0    9.839000912153142E+05
!  9.0  3.5    1.173074963158710E+07
!  9.0  7.5    4.897833845586907E+08
!  9.0 15.0    1.080715538012005E+11
!  9.0 30.0    6.948254776893934E+13
!  9.0 60.0    6.299767361155796E+16
!
!  9.5 -3.0    5.642067020436533E+04
!  9.5 -1.0    4.168044535669152E+05
!  9.5  1.0    3.074986313749142E+06
!  9.5  3.5    3.686910477586340E+07
!  9.5  7.5    1.615856749616024E+09
!  9.5 15.0    4.273879330919564E+11
!  9.5 30.0    3.689532972232714E+14
!  9.5 60.0    4.668155109009914E+17
!
! 10.0 -3.0    1.806629241770091E+05
! 10.0 -1.0    1.334722123421517E+06
! 10.0  1.0    9.851381090333600E+06
! 10.0  3.5    1.186176377953597E+08
! 10.0  7.5    5.416646367586654E+09
! 10.0 15.0    1.701154132363905E+12
! 10.0 30.0    1.965505543248569E+15
! 10.0 60.0    3.467790683237043E+18
!
! 10.5 -3.0    5.924272114982210E+05
! 10.5 -1.0    4.376998997938984E+06
! 10.5  1.0    3.231634166275790E+07
! 10.5  3.5    3.903364896044762E+08
! 10.5  7.5    1.845016484905844E+10
! 10.5 15.0    6.814687278457835E+12
! 10.5 30.0    1.050272666569349E+16
! 10.5 60.0    2.582022264434903E+19
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
