c****************************************************************************** SUBROUTINE fermi_dirac(x,fd) c calcul des intégrales de Fermi-Dirac from MHD package c routine public du module mod_numerique c fd(1) = f (x) c -1/2 c fd(2) = f (x) c 1/2 c fd(3) = f (x) c 3/2 c fd(4) = f' (x) c -1/2 c fd(5) = f" (x) c -1/2 c version F95 : P.Morel, Département J.D. Cassini, O.C.A. c-------------------------------------------------------------------------- USE mod_kind IMPLICIT NONE REAL (kind=dp), INTENT(in) :: x REAL (kind=dp), INTENT(out) :: fd(5) REAL (kind=dp), DIMENSION(5), SAVE :: p1, p2, p3, p4, p5, p6, 1 p7, p8, p9, q1, q2, q3, q4, q5, q6, q7, q8, q9 REAL (kind=dp) :: y, p, dr, d2p, q, dq, d2q, root, xsq DATA p1 1/-1.253314128820d+0,-1.723663557701d+0,-6.559045729258d-1, 2 -6.342283197682d-2,-1.488383106116d-5/ DATA q1 1/+1.000000000000d+0,+2.191780925980d+0,+1.605812955406d+0, 2 +4.443669527481d-1,+3.624232288112d-2/ DATA p2 1/-3.133285305570d-1,-4.161873852293d-1,-1.502208400588d-1, 2 -1.339579375173d-2,-1.513350700138d-5/ DATA q2 1/+1.000000000000d+0,+1.872608675902d+0,+1.145204446578d+0, 2 +2.570225587573d-1,+1.639902543568d-2/ DATA p3 1/-2.349963985406d-1,-2.927373637547d-1,-9.883097588738d-2, 2 -8.251386379551d-3,-1.874384153223d-5/ DATA q3 1/+1.000000000000d+0,+1.608597109146d+0,+8.275289530880d-1, 2 +1.522322382850d-1,+7.695120475064d-3/ DATA p4 1/+1.073812769400d+0,+5.600330366000d+0,+3.688221127000d+0, 2 +1.174339281600d+0,+2.364193552700d-1/ DATA q4 1/+1.000000000000d+0,+4.603184066700d+0,+4.307591067400d-1, 2 +4.215113214500d-1,+1.183260160100d-2/ DATA p5 1/+6.781766266600d-1,+6.331240179100d-1,+2.944796517720d-1, 2 +8.013207114190d-2,+1.339182129400d-2/ DATA q5 1/+1.000000000000d+0,+1.437404003970d-1,+7.086621484500d-2, 2 +2.345794947350d-3,-1.294499288350d-5/ DATA p6 1/+1.153021340200d+0,+1.059155897200d+0,+4.689880309500d-1, 2 +1.188290878400d-1,+1.943875578700d-2/ DATA q6 1/+1.000000000000d+0,+3.734895384100d-2,+2.324845813700d-2, 2 -1.376677087400d-3, +4.646639278100d-5/ DATA p7 1/-8.222559330000d-1, -3.620369345000d+1,-3.015385410000d+3, 2 -7.049871579000d+4, -5.698145924000d+4/ DATA q7 1/+1.000000000000d+0, +3.935689841000d+1,+3.568756266000d+3, 2 +4.181893625000d+4, +3.385138907000d+5/ DATA p8 1/+8.224499762600d-1, +2.004630339300d+1,+1.826809344600d+3, 2 +1.222653037400d+4, +1.404075009200d+5/ DATA q8 1/+1.000000000000d+0, +2.348620765900d+1,+2.201348374300d+3, 2 +1.144267359600d+4, +1.658471590000d+5/ DATA p9 1/+2.467400236840d+0, +2.191675823680d+2,+1.238293790750d+4, 2 +2.206677249680d+5, +8.494429200340d+5/ DATA q9 1/+1.000000000000d+0, +8.911251406190d+1,+5.045756696670d+3, 2 +9.090759463040d+4, +3.899609156410d+5/ c-------------------------------------------------------------------------- IF(x < 1.d0)THEN y= dexp(x) p = y**2*( p1(1) + y*( p1(2) + y*( p1(3) + y*( p1(4) 1 + y* p1(5))))) dr = y**2*(2.*p1(1) + y*(3.*p1(2) + y*( 4.*p1(3) + y*( 5.*p1(4) 1 + y* 6.*p1(5))))) d2p = y**2*(4.*p1(1) + y*(9.*p1(2) + y*(16.*p1(3) + y*(25.*p1(4) 1 + y*36.*p1(5))))) q = q1(1) +y*(q1(2) + y*( q1(3) + y*( q1(4) + y* q1(5)))) dq = y*(q1(2) + y*(2.*q1(3) + y*(3.*q1(4) + y* 4.*q1(5)))) d2q = y*(q1(2) + y*(4.*q1(3) + y*(9.*q1(4) + y*16.*q1(5)))) fd(1) = 1.7724 53850 90552 d0*y + p/q fd(2) = y*(0.8862 26925 45276 d0 + y* 1 (p2(1) + y*(p2(2) + y*(p2(3) + y*(p2(4) + y*p2(5)))))/ 2 (q2(1) + y*(q2(2) + y*(q2(3) + y*(q2(4) + y*q2(5)))))) fd(3) = y*(1.3293 40388 17914 d0 + y* 1 (p3(1) + y*(p3(2) + y*(p3(3) + y*(p3(4) + y*p3(5)))))/ 2 (q3(1) + y*(q3(2) + y*(q3(3) + y*(q3(4) + y*q3(5)))))) fd(4) = 1.7724 53850 90552 d0*y + (dr*q - p*dq)/q**2 fd(5) = 1.7724 53850 90552 d0*y + 1 ((d2p*q - p*d2q)*q - 2.*(dr*q -p*dq)*dq)/q**3 ELSEIF(x < 4.d0)THEN p = p4(1)+ x*(p4(2) + x*( p4(3) + x*( p4(4) + x* p4(5)))) dr = p4(2) + x*(2.*p4(3) + x*(3.*p4(4) + x*4.*p4(5))) d2p = 2.*p4(3) + 6.*x*(p4(4) + x*2.*p4(5)) q = q4(1)+ x*(q4(2) + x*( q4(3) + x*( q4(4) + x* q4(5)))) dq = q4(2) + x*(2.*q4(3) + x*(3.*q4(4) + x*4.*q4(5))) d2q = 2.*q4(3) + 6.*x*(q4(4) + x*2.*q4(5)) fd(1) = p/q fd(2) = (p5(1) + x*(p5(2) + x*(p5(3) + x*(p5(4) + x*p5(5)))))/ 1 (q5(1) + x*(q5(2) + x*(q5(3) + x*(q5(4) + x*q5(5))))) fd(3) = (p6(1) + x*(p6(2) + x*(p6(3) + x*(p6(4) + x*p6(5)))))/ 1 (q6(1) + x*(q6(2) + x*(q6(3) + x*(q6(4) + x*q6(5))))) fd(4) = (dr *q - p*dq )/q**2 fd(5) = (d2p*q - p*d2q)/q**2 - 2.*fd(4)*dq/q ELSE root = sqrt(x) xsq = x**2 y = 1.0d0/xsq p = y * (p7(1) + y*( p7(2) + y*( p7(3) + y*( p7(4) 1 + y* p7(5))))) dr = (-2.*y/x)* (p7(1) + y*(2.*p7(2) + y*(3.*p7(3) + y*( 4.*p7(4) 1 + y* 5.*p7(5))))) d2p = 4.*y**2 * (p7(1) + y*(4.*p7(2) + y*(9.*p7(3) + y*(16.*p7(4) 1 + y*25.*p7(5))))) - dr/x q = q7(1) + y * (q7(2) +y*( q7(3) +y*( q7(4) +y* q7(5)))) dq = (-2.*y/x) * (q7(2) +y*(2.*q7(3) +y*(3.*q7(4) +y* 4.*q7(5)))) d2q = 4.*y**2 * (q7(2) +y*(4.*q7(3) +y*(9.*q7(4) +y*16.*q7(5)))) 1 - dq/x fd(1) = root * (2.0d0 + p/q) fd(2) = x*root*(0.66666 66666 66667 d0 + y* 1 (p8(1) + y*(p8(2) + y*(p8(3) + y*(p8(4) + y*p8(5)))))/ 2 (q8(1) + y*(q8(2) + y*(q8(3) + y*(q8(4) + y*q8(5)))))) fd(3) = xsq*root*(0.4d0 + y* 1 (p9(1) + y*(p9(2) + y*(p9(3) + y*(p9(4) + y*p9(5)))))/ 2 (q9(1) + y*(q9(2) + y*(q9(3) + y*(q9(4) + y*q9(5)))))) fd(4) = fd(1)/(2.d0*x) + root*(dr*q - p*dq)/q**2 fd(5) = ( -fd(1)/x + fd(4))/(2.d0*x) 1 + root* (dr *q - p*dq )/(2.*x*q**2) 2 + root*((d2p*q - p*d2q) - 2.*dq*(dr*q - p*dq)/q)/q**2 ENDIF RETURN END SUBROUTINE fermi_dirac