function tetrahedron_fermidirac(F, E_ref, E2, t, hw, Ef, type, tet_cutoff, avoid_deg)
!=============================================================================
! Calculates contribution from a single tetrahedron, for the Kubo formula. !
! "integral d3k f_{nk} * (...)
!=============================================================================
use w90_constants, only: dp
implicit none
real(kind=dp) :: tetrahedron_fermidirac
real(kind=dp), dimension(4), intent(in) :: F, E_ref, E2
real(kind=dp), dimension(3, 3), intent(in) :: t
real(kind=dp), intent(in) :: hw, Ef
integer, intent(in) :: type
real(kind=dp), intent(in) :: tet_cutoff
real(kind=dp), intent(in) :: avoid_deg
!intermediate variables, s:sorted
real(kind=dp) :: t_s(3, 3), t_small(3, 3), x(3), y
real(kind=dp), dimension(4) :: D, F_s, E1_s, E2_s, F_small, D_small
!result
real(kind=dp) :: Ans
integer :: i
Ans = 0.0_dp
!sorting vertices according to E_ref
F_s = F; E1_s = E_ref; E2_s = E2; t_s = t !_s: sorted
call tetrahedron_sort(E1_s, E2_s, F_s, t_s)
D = E1_s - E2_s
!case 1,2,3,4,5
if (Ef < E1_s(1)) then ! case 1: zero
Ans = Ans + 0.0_dp
else if (Ef < E1_s(2)) then ! case 2: a small tet.
x(1) = (Ef - E1_s(1))/(E1_s(2) - E1_s(1))
x(2) = (Ef - E1_s(1))/(E1_s(3) - E1_s(1))
x(3) = (Ef - E1_s(1))/(E1_s(4) - E1_s(1))
F_small(1) = F_s(1); D_small(1) = D(1)
do i = 1, 3
F_small(i + 1) = F_s(1) + (F_s(i + 1) - F_s(1))*x(i)
D_small(i + 1) = D(1) + (D(i + 1) - D(1))*x(i)
t_small(:, i) = t_s(:, i)*x(i)
end do
Ans = Ans + tetrahedron_integral(F_small, D_small, t_small, hw, type, tet_cutoff, avoid_deg)
else if (Ef < E1_s(3)) then ! case 3: two tet.'s with cases 2 and 4
x(1) = (Ef - E1_s(4))/(E1_s(2) - E1_s(4))
x(2) = (Ef - E1_s(1))/(E1_s(3) - E1_s(1))
x(3) = (Ef - E1_s(1))/(E1_s(4) - E1_s(1))
y = (Ef - E1_s(3))/(E1_s(2) - E1_s(3))
F_small = F_s; D_small = D; t_small = t_s
F_small(4) = F_s(1) + (F_s(4) - F_s(1))*x(3)
D_small(4) = D(1) + (D(4) - D(1))*x(3)
t_small(:, 3) = t_s(:, 3)*x(3)
Ans = Ans + tetrahedron_integral(F_small, D_small, t_small, hw, type, tet_cutoff, avoid_deg)
F_small(1) = F_s(1) + (F_s(3) - F_s(1))*x(2); F_small(2) = F_s(3) + (F_s(2) - F_s(3))*y
D_small(1) = D(1) + (D(3) - D(1))*x(2); D_small(2) = D(3) + (D(2) - D(3))*y
t_small(:, 1) = t_s(:, 1)*y + t_s(:, 2)*(1 - y - x(2))
t_small(:, 2) = t_s(:, 2)*(1 - x(2))
t_small(:, 3) = t_s(:, 3)*x(3) - t_s(:, 2)*x(2)
Ans = Ans - tetrahedron_integral(F_small, D_small, t_small, hw, type, tet_cutoff, avoid_deg)
F_small(1) = F_s(4) + (F_s(2) - F_s(4))*x(1)
F_small(3) = F_small(2); F_small(2) = F_s(2)
D_small(1) = D(4) + (D(2) - D(4))*x(1)
D_small(3) = D_small(2); D_small(2) = D(2)
t_small(:, 1) = (t_s(:, 1) - t_s(:, 3))*(1 - x(1))
t_small(:, 2) = t_s(:, 1)*(y - x(1)) + t_s(:, 2)*(1 - y) + t_s(:, 3)*(x(1) - 1)
t_small(:, 3) = -t_s(:, 1)*x(1) + t_s(:, 3)*(x(1) + x(3) - 1)
Ans = Ans + tetrahedron_integral(F_small, D_small, t_small, hw, type, tet_cutoff, avoid_deg)
else if (Ef < E1_s(4)) then ! case 4: a large tet. - a small tet.
x(1) = (Ef - E1_s(4))/(E1_s(2) - E1_s(4))
x(2) = (Ef - E1_s(4))/(E1_s(3) - E1_s(4))
x(3) = (Ef - E1_s(1))/(E1_s(4) - E1_s(1))
F_small(1) = F_s(4); F_small(2) = F_s(1) + (F_s(4) - F_s(1))*x(3)
F_small(3) = F_s(4) + (F_s(2) - F_s(4))*x(1)
F_small(4) = F_s(4) + (F_s(3) - F_s(4))*x(2)
D_small(1) = D(4); D_small(2) = D(1) + (D(4) - D(1))*x(3)
D_small(3) = D(4) + (D(2) - D(4))*x(1)
D_small(4) = D(4) + (D(3) - D(4))*x(2)
t_small(:, 1) = -t_s(:, 3)*(1 - x(3))
t_small(:, 2) = (t_s(:, 1) - t_s(:, 3))*x(1)
t_small(:, 3) = (t_s(:, 2) - t_s(:, 3))*x(2)
Ans = Ans + tetrahedron_integral(F_s, D, t_s, hw, type, tet_cutoff, avoid_deg) &
- tetrahedron_integral(F_small, D_small, t_small, hw, type, tet_cutoff, avoid_deg)
else ! case 5: a large tet.
Ans = Ans + tetrahedron_integral(F_s, D, t_s, hw, type, tet_cutoff, avoid_deg)
end if
tetrahedron_fermidirac = Ans
end function tetrahedron_fermidirac
