"""
For more details on the demonstrated formula below, please find the following references:

Ref 1) Extensions of the tetrahedron method for evaluating spectral properties of solids, 
Journal of Physics C 12, 2991 (1979).
Ref 2) FermiSurfer: Fermi-surface viewer providing multiple representation schemes,
Computer Physics Communications 239, 197 (2019).
"""
import numpy as np
from typing import List


def tetrahedron_1d(_freq: np.ndarray,
                   _pdos: np.ndarray,
                   k_points: List[np.ndarray],
                   w_q: np.ndarray,
                   v_q: np.ndarray,
                   num_k_points: List[int],
                   _ind_pbc: np.ndarray) -> np.ndarray:
    """
    The 1D version of linear tetrahedron method: **linear line method**.

    :param _freq: '(num_dos,) size' frequencies where DOS will be evaluated
    :type _freq: np.ndarray[float]
    :param _pdos: '(num_mode, num_dos) size' partial DOS for each phonon mode
    :type _pdos: np.ndarray[float]
    :param k_points: List of '(3,) size' k-point
    :type k_points: List[np.ndarray]
    :param w_q: '(num_k_points, num_mode) size' eigen-frequency
    :type w_q: np.ndarray[float]
    :param v_q: '(num_k_points, num_mode, num_mode) size' eigen-mode
    :type v_q: np.ndarray[float]
    :param num_k_points: '(3,) size' list for the number of automatic grid
    :type num_k_points: List[int]
    :param _ind_pbc: Indices for the lattice direction of periodic boundary
    :type _ind_pbc: np.ndarray[int]
    :return: '(num_mode, num_dos) size' partial DOS for each phonon mode
    :rtype: np.ndarray[float]
    """
    #  1st Brillouin Zone is segmented by line components connecting the sampled k-points grid.
    #  Integration is performed for each segmentation by linear interpolation.
    #
    #  line components:
    #  0-------–1
    #
    #  line 1: 0-1
    num_1st_kpt = num_k_points[_ind_pbc[0]]
    length_bz = 1.0 ** 1
    length_ith = np.sqrt(np.dot(k_points[1 % num_1st_kpt] - k_points[0], k_points[1 % num_1st_kpt] - k_points[0]))

    for ind_x, __freq in enumerate(_freq):
        for ind_freq in range(v_q.shape[1]):
            _w_q = w_q[:, ind_freq]
            _v_q = v_q[:, ind_freq, :]
            for ind_1st_kpt in range(num_1st_kpt):
                # line 1: 0-1
                # for ind_mode in range(v_q.shape[2]):
                __w0 = _w_q[ind_1st_kpt]
                __F0 = abs(_v_q[ind_1st_kpt, :]) ** 2
                __w1 = _w_q[(ind_1st_kpt + 1) % num_1st_kpt]
                __F1 = abs(_v_q[(ind_1st_kpt + 1) % num_1st_kpt, :]) ** 2

                __w0, __w1 = np.asfarray([__w0, __w1])[np.argsort(np.asfarray([__w0, __w1]))]
                __F0, __F1 = np.asfarray([__F0, __F1])[np.argsort(np.asfarray([__w0, __w1])), :]

                if __w0 <= __freq < __w1:
                    __g = 1.0 / (__w1 - __w0)

                    __I0 = (__freq - __w1) / (__w0 - __w1)
                    __I1 = (__freq - __w0) / (__w1 - __w0)

                    _pdos[:, ind_x] += length_ith / length_bz * __g \
                                                 * (__I0 * __F0 + __I1 * __F1)

                else:
                    pass

    return _pdos


def tetrahedron_2d(_freq: np.ndarray,
                   _pdos: np.ndarray,
                   k_points: List[np.ndarray],
                   w_q: np.ndarray,
                   v_q: np.ndarray,
                   num_k_points: List[int],
                   _ind_pbc: np.ndarray) -> np.ndarray:
    """
    The 2D version of linear tetrahedron method: **linear triangle method**.

    :param _freq: '(num_dos,) size' frequencies where DOS will be evaluated
    :type _freq: np.ndarray[float]
    :param _pdos: '(num_mode, num_dos) size' partial DOS for each phonon mode
    :type _pdos: np.ndarray[float]
    :param k_points: List of '(3,) size' k-point
    :type k_points: List[np.ndarray]
    :param w_q: '(num_k_points, num_mode) size' eigen-frequency
    :type w_q: np.ndarray[float]
    :param v_q: '(num_k_points, num_mode, num_mode) size' eigen-mode
    :type v_q: np.ndarray[float]
    :param num_k_points: '(3,) size' list for the number of automatic grid
    :type num_k_points: List[int]
    :param _ind_pbc: Indices for the lattice direction of periodic boundary
    :type _ind_pbc: np.ndarray[int]
    :return: '(num_mode, num_dos) size' partial DOS for each phonon mode
    :rtype: np.ndarray[float]
    """
    #  1st Brillouin Zone is segmented by triangle components connecting the sampled k-points grid.
    #  Integration is performed for each segmentation by linear interpolation.
    #
    #  triangle components:
    #  0--------1
    #  |        |
    #  |        |
    #  |        |
    #  2-------–3
    #
    # triangle 1: 0-1-3
    # triangle 2: 0-2-3
    num_1st_kpt = num_k_points[_ind_pbc[0]]
    num_2nd_kpt = num_k_points[_ind_pbc[1]]

    __k01 = k_points[1 % num_2nd_kpt] - k_points[0]
    __k02 = k_points[num_2nd_kpt] - k_points[0]

    __tmp = np.ones((3,))
    __tmp[_ind_pbc[0]] = 0.0
    __tmp[_ind_pbc[1]] = 0.0

    area_bz = 1.0 ** 2
    area_ith = 1.0 / 2.0 * abs(np.dot(__tmp, np.cross(__k01, __k02)))

    for ind_x, __freq in enumerate(_freq):
        for ind_freq in range(v_q.shape[1]):
            _w_q = w_q[:, ind_freq]
            _v_q = v_q[:, ind_freq, :]
            for ind_1st_kpt in range(num_1st_kpt):
                for ind_2nd_kpt in range(num_2nd_kpt):
                    # triangle 1: 0-1-3
                    # for ind_mode in range(v_q.shape[2]):
                    __w0 = _w_q[ind_1st_kpt * num_2nd_kpt +
                                ind_2nd_kpt]
                    __F0 = abs(_v_q[ind_1st_kpt * num_2nd_kpt +
                                    ind_2nd_kpt, :]) ** 2
                    __w1 = _w_q[ind_1st_kpt * num_2nd_kpt +
                                (ind_2nd_kpt + 1) % num_2nd_kpt]
                    __F1 = abs(_v_q[ind_1st_kpt * num_2nd_kpt +
                                    (ind_2nd_kpt + 1) % num_2nd_kpt, :]) ** 2
                    __w2 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt +
                                (ind_2nd_kpt + 1) % num_2nd_kpt]
                    __F2 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt +
                                    (ind_2nd_kpt + 1) % num_2nd_kpt, :]) ** 2

                    __w0, __w1, __w2 = np.asfarray([__w0, __w1, __w2])[np.argsort(np.asfarray([__w0, __w1, __w2]))]
                    __F0, __F1, __F2 = np.asfarray([__F0, __F1, __F2])[np.argsort(np.asfarray([__w0, __w1, __w2])), :]

                    if __w0 <= __freq < __w1:
                        __g = 2 * (__freq - __w0) / (__w1 - __w0) / (__w2 - __w0)

                        __I0 = 1.0 / 2.0 * ((__freq - __w1) / (__w0 - __w1) + (__freq - __w2) / (__w0 - __w2))
                        __I1 = 1.0 / 2.0 * (__freq - __w0) / (__w1 - __w0)
                        __I2 = 1.0 / 2.0 * (__freq - __w0) / (__w2 - __w0)

                        _pdos[:, ind_x] += area_ith / area_bz * __g \
                                                     * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2)

                    elif __w1 <= __freq < __w2:
                        __g = 2 * (__w2 - __freq) / (__w2 - __w1) / (__w2 - __w0)

                        __I0 = 1.0 / 2.0 * (__freq - __w2) / (__w0 - __w2)
                        __I1 = 1.0 / 2.0 * (__freq - __w2) / (__w1 - __w2)
                        __I2 = 1.0 / 2.0 * ((__freq - __w0) / (__w2 - __w0) + (__freq - __w1) / (__w2 - __w1))

                        _pdos[:, ind_x] += area_ith / area_bz * __g \
                                                     * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2)

                    else:
                        pass

                    # triangle 2: 0-2-3
                    # for ind_mode in range(v_q.shape[2]):
                    __w0 = _w_q[ind_1st_kpt * num_2nd_kpt +
                                ind_2nd_kpt]
                    __F0 = abs(_v_q[ind_1st_kpt * num_2nd_kpt +
                                    ind_2nd_kpt, :]) ** 2
                    __w1 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt +
                                ind_2nd_kpt]
                    __F1 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt +
                                    ind_2nd_kpt, :]) ** 2
                    __w2 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt +
                                (ind_2nd_kpt + 1) % num_2nd_kpt]
                    __F2 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt +
                                    (ind_2nd_kpt + 1) % num_2nd_kpt, :]) ** 2

                    __w0, __w1, __w2 = np.asfarray([__w0, __w1, __w2])[np.argsort(np.asfarray([__w0, __w1, __w2]))]
                    __F0, __F1, __F2 = np.asfarray([__F0, __F1, __F2])[np.argsort(np.asfarray([__w0, __w1, __w2])), :]

                    if __w0 <= __freq < __w1:
                        __g = 2 * (__freq - __w0) / (__w1 - __w0) / (__w2 - __w0)

                        __I0 = 1.0 / 2.0 * ((__freq - __w1) / (__w0 - __w1) + (__freq - __w2) / (__w0 - __w2))
                        __I1 = 1.0 / 2.0 * (__freq - __w0) / (__w1 - __w0)
                        __I2 = 1.0 / 2.0 * (__freq - __w0) / (__w2 - __w0)

                        _pdos[:, ind_x] += area_ith / area_bz * __g \
                                                     * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2)

                    elif __w1 <= __freq < __w2:
                        __g = 2 * (__w2 - __freq) / (__w2 - __w1) / (__w2 - __w0)

                        __I0 = 1.0 / 2.0 * (__freq - __w2) / (__w0 - __w2)
                        __I1 = 1.0 / 2.0 * (__freq - __w2) / (__w1 - __w2)
                        __I2 = 1.0 / 2.0 * ((__freq - __w0) / (__w2 - __w0) + (__freq - __w1) / (__w2 - __w1))

                        _pdos[:, ind_x] += area_ith / area_bz * __g \
                                                     * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2)

                    else:
                        pass

    return _pdos


def tetrahedron_3d(_freq: np.ndarray,
                   _pdos: np.ndarray,
                   k_points: List[np.ndarray],
                   w_q: np.ndarray,
                   v_q: np.ndarray,
                   num_k_points: List[int],
                   _ind_pbc: np.ndarray) -> np.ndarray:
    """
    The 3D version of linear tetrahedron method: **linear tetrahedron method**.

    :param _freq: '(num_dos,) size' frequencies where DOS will be evaluated
    :type _freq: np.ndarray[float]
    :param _pdos: '(num_mode, num_dos) size' partial DOS for each phonon mode
    :type _pdos: np.ndarray[float]
    :param k_points: List of '(3,) size' k-point
    :type k_points: List[np.ndarray]
    :param w_q: '(num_k_points, num_mode) size' eigen-frequency
    :type w_q: np.ndarray[float]
    :param v_q: '(num_k_points, num_mode, num_mode) size' eigen-mode
    :type v_q: np.ndarray[float]
    :param num_k_points: '(3,) size' list for the number of automatic grid
    :type num_k_points: List[int]
    :param _ind_pbc: Indices for the lattice direction of periodic boundary
    :type _ind_pbc: np.ndarray[int]
    :return: '(num_mode, num_dos) size' partial DOS for each phonon mode
    :rtype: np.ndarray[float]
    """
    #  1st Brillouin Zone is segmented by tetrahedron components connecting the sampled k-points grid.
    #  Integration is performed for each segmentation by linear interpolation.
    #
    #  tetrahedron components:
    #     4--------5
    #    /|       /|
    #   / |      / |
    #  6--------7  |
    #  |  0-----|--1
    #  | /      | /
    #  |/       |/
    #  2--------3
    #
    # tetrahedron 1: 0-1-5-7
    # tetrahedron 2: 0-1-3-7
    # tetrahedron 3: 0-2-3-7
    # tetrahedron 4: 0-4-5-7
    # tetrahedron 5: 0-4-6-7
    # tetrahedron 6: 0-2-6-7
    num_1st_kpt = num_k_points[_ind_pbc[0]]
    num_2nd_kpt = num_k_points[_ind_pbc[1]]
    num_3rd_kpt = num_k_points[_ind_pbc[2]]

    __k01 = k_points[1 % num_2nd_kpt] - k_points[0]
    __k02 = k_points[num_2nd_kpt] - k_points[0]
    __k03 = k_points[num_2nd_kpt * num_3rd_kpt] - k_points[0]

    volume_bz = 1.0 ** 3
    volume_ith = 1.0 / 6.0 * abs(np.dot(__k01, np.cross(__k02, __k03)))

    for ind_x, __freq in enumerate(_freq):
        for ind_freq in range(v_q.shape[1]):
            _w_q = w_q[:, ind_freq]
            _v_q = v_q[:, ind_freq, :]
            for ind_1st_kpt in range(num_1st_kpt):
                for ind_2nd_kpt in range(num_2nd_kpt):
                    for ind_3rd_kpt in range(num_3rd_kpt):
                        # tetrahedron 1: 0-1-5-7
                        # for ind_mode in range(v_q.shape[2]):
                        __w0 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F0 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w1 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F1 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2
                        __w2 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F2 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2
                        __w3 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F3 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2

                        __w0, __w1, __w2, __w3 = np.asfarray([__w0,
                                                              __w1,
                                                              __w2,
                                                              __w3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3]))]
                        __F0, __F1, __F2, __F3 = np.asfarray([__F0,
                                                              __F1,
                                                              __F2,
                                                              __F3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3])), :]

                        if __w0 <= __freq < __w1:
                            __g = 3 * (__freq - __w0) ** 2 / (__w1 - __w0) / (__w2 - __w0) / (__w3 - __w0)

                            __I0 = 1.0 / 3.0 * ((__freq - __w1) / (__w0 - __w1)
                                                + (__freq - __w2) / (__w0 - __w2)
                                                + (__freq - __w3) / (__w0 - __w3))
                            __I1 = 1.0 / 3.0 * (__freq - __w0) / (__w1 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w0) / (__w2 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                                         * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w1 <= __freq < __w2:
                            __g = 3 / (__w3 - __w0) \
                                  * ((__freq - __w1) * (__freq - __w3) / ((__w2 - __w1) * (__w1 - __w3))
                                     + (__freq - __w0) * (__freq - __w2) / ((__w2 - __w0) * (__w1 - __w2)))

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3) \
                                   + (__freq - __w2) / (__w0 - __w2) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I1 = 1.0 / 3.0 * (__freq - __w2) / (__w1 - __w2) \
                                   + (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w1) / (__w2 - __w1) \
                                   + (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0) \
                                   + (__freq - __w1) / (__w3 - __w1) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                                         * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w2 <= __freq < __w3:
                            __g = 3 * (__w3 - __freq) ** 2 / (__w3 - __w0) / (__w3 - __w1) / (__w3 - __w2)

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3)
                            __I1 = 1.0 / 3.0 * (__freq - __w3) / (__w1 - __w3)
                            __I2 = 1.0 / 3.0 * (__freq - __w3) / (__w2 - __w3)
                            __I3 = 1.0 / 3.0 * ((__freq - __w0) / (__w3 - __w0)
                                                + (__freq - __w1) / (__w3 - __w1)
                                                + (__freq - __w2) / (__w3 - __w2))

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                                         * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        else:
                            pass

                        # tetrahedron 2: 0-1-3-7
                        # for ind_mode in range(v_q.shape[2]):
                        __w0 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F0 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w1 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F1 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2
                        __w2 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F2 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2
                        __w3 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F3 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2

                        __w0, __w1, __w2, __w3 = np.asfarray([__w0,
                                                              __w1,
                                                              __w2,
                                                              __w3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3]))]
                        __F0, __F1, __F2, __F3 = np.asfarray([__F0,
                                                              __F1,
                                                              __F2,
                                                              __F3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3])), :]

                        if __w0 <= __freq < __w1:
                            __g = 3 * (__freq - __w0) ** 2 / (__w1 - __w0) / (__w2 - __w0) / (__w3 - __w0)

                            __I0 = 1.0 / 3.0 * ((__freq - __w1) / (__w0 - __w1)
                                                + (__freq - __w2) / (__w0 - __w2)
                                                + (__freq - __w3) / (__w0 - __w3))
                            __I1 = 1.0 / 3.0 * (__freq - __w0) / (__w1 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w0) / (__w2 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w1 <= __freq < __w2:
                            __g = 3 / (__w3 - __w0) \
                                  * ((__freq - __w1) * (__freq - __w3) / ((__w2 - __w1) * (__w1 - __w3))
                                     + (__freq - __w0) * (__freq - __w2) / ((__w2 - __w0) * (__w1 - __w2)))

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3) \
                                   + (__freq - __w2) / (__w0 - __w2) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I1 = 1.0 / 3.0 * (__freq - __w2) / (__w1 - __w2) \
                                   + (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w1) / (__w2 - __w1) \
                                   + (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0) \
                                   + (__freq - __w1) / (__w3 - __w1) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w2 <= __freq < __w3:
                            __g = 3 * (__w3 - __freq) ** 2 / (__w3 - __w0) / (__w3 - __w1) / (__w3 - __w2)

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3)
                            __I1 = 1.0 / 3.0 * (__freq - __w3) / (__w1 - __w3)
                            __I2 = 1.0 / 3.0 * (__freq - __w3) / (__w2 - __w3)
                            __I3 = 1.0 / 3.0 * ((__freq - __w0) / (__w3 - __w0)
                                                + (__freq - __w1) / (__w3 - __w1)
                                                + (__freq - __w2) / (__w3 - __w2))

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        else:
                            pass

                        # tetrahedron 3: 0-2-3-7
                        # for ind_mode in range(v_q.shape[2]):
                        __w0 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F0 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w1 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F1 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w2 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F2 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2
                        __w3 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F3 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2

                        __w0, __w1, __w2, __w3 = np.asfarray([__w0,
                                                              __w1,
                                                              __w2,
                                                              __w3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3]))]
                        __F0, __F1, __F2, __F3 = np.asfarray([__F0,
                                                              __F1,
                                                              __F2,
                                                              __F3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3])), :]

                        if __w0 <= __freq < __w1:
                            __g = 3 * (__freq - __w0) ** 2 / (__w1 - __w0) / (__w2 - __w0) / (__w3 - __w0)

                            __I0 = 1.0 / 3.0 * ((__freq - __w1) / (__w0 - __w1)
                                                + (__freq - __w2) / (__w0 - __w2)
                                                + (__freq - __w3) / (__w0 - __w3))
                            __I1 = 1.0 / 3.0 * (__freq - __w0) / (__w1 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w0) / (__w2 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w1 <= __freq < __w2:
                            __g = 3 / (__w3 - __w0) \
                                  * ((__freq - __w1) * (__freq - __w3) / ((__w2 - __w1) * (__w1 - __w3))
                                     + (__freq - __w0) * (__freq - __w2) / ((__w2 - __w0) * (__w1 - __w2)))

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3) \
                                   + (__freq - __w2) / (__w0 - __w2) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I1 = 1.0 / 3.0 * (__freq - __w2) / (__w1 - __w2) \
                                   + (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w1) / (__w2 - __w1) \
                                   + (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0) \
                                   + (__freq - __w1) / (__w3 - __w1) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w2 <= __freq < __w3:
                            __g = 3 * (__w3 - __freq) ** 2 / (__w3 - __w0) / (__w3 - __w1) / (__w3 - __w2)

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3)
                            __I1 = 1.0 / 3.0 * (__freq - __w3) / (__w1 - __w3)
                            __I2 = 1.0 / 3.0 * (__freq - __w3) / (__w2 - __w3)
                            __I3 = 1.0 / 3.0 * ((__freq - __w0) / (__w3 - __w0)
                                                + (__freq - __w1) / (__w3 - __w1)
                                                + (__freq - __w2) / (__w3 - __w2))

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        else:
                            pass

                        # tetrahedron 4: 0-4-5-7
                        # for ind_mode in range(v_q.shape[2]):
                        __w0 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F0 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w1 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F1 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w2 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F2 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2
                        __w3 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F3 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2

                        __w0, __w1, __w2, __w3 = np.asfarray([__w0,
                                                              __w1,
                                                              __w2,
                                                              __w3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3]))]
                        __F0, __F1, __F2, __F3 = np.asfarray([__F0,
                                                              __F1,
                                                              __F2,
                                                              __F3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3])), :]

                        if __w0 <= __freq < __w1:
                            __g = 3 * (__freq - __w0) ** 2 / (__w1 - __w0) / (__w2 - __w0) / (__w3 - __w0)

                            __I0 = 1.0 / 3.0 * ((__freq - __w1) / (__w0 - __w1)
                                                + (__freq - __w2) / (__w0 - __w2)
                                                + (__freq - __w3) / (__w0 - __w3))
                            __I1 = 1.0 / 3.0 * (__freq - __w0) / (__w1 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w0) / (__w2 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w1 <= __freq < __w2:
                            __g = 3 / (__w3 - __w0) \
                                  * ((__freq - __w1) * (__freq - __w3) / ((__w2 - __w1) * (__w1 - __w3))
                                     + (__freq - __w0) * (__freq - __w2) / ((__w2 - __w0) * (__w1 - __w2)))

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3) \
                                   + (__freq - __w2) / (__w0 - __w2) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I1 = 1.0 / 3.0 * (__freq - __w2) / (__w1 - __w2) \
                                   + (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w1) / (__w2 - __w1) \
                                   + (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0) \
                                   + (__freq - __w1) / (__w3 - __w1) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w2 <= __freq < __w3:
                            __g = 3 * (__w3 - __freq) ** 2 / (__w3 - __w0) / (__w3 - __w1) / (__w3 - __w2)

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3)
                            __I1 = 1.0 / 3.0 * (__freq - __w3) / (__w1 - __w3)
                            __I2 = 1.0 / 3.0 * (__freq - __w3) / (__w2 - __w3)
                            __I3 = 1.0 / 3.0 * ((__freq - __w0) / (__w3 - __w0)
                                                + (__freq - __w1) / (__w3 - __w1)
                                                + (__freq - __w2) / (__w3 - __w2))

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        else:
                            pass

                        # tetrahedron 5: 0-4-6-7
                        # for ind_mode in range(v_q.shape[2]):
                        __w0 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F0 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w1 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F1 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w2 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F2 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w3 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F3 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2

                        __w0, __w1, __w2, __w3 = np.asfarray([__w0,
                                                              __w1,
                                                              __w2,
                                                              __w3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3]))]
                        __F0, __F1, __F2, __F3 = np.asfarray([__F0,
                                                              __F1,
                                                              __F2,
                                                              __F3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3])), :]

                        if __w0 <= __freq < __w1:
                            __g = 3 * (__freq - __w0) ** 2 / (__w1 - __w0) / (__w2 - __w0) / (__w3 - __w0)

                            __I0 = 1.0 / 3.0 * ((__freq - __w1) / (__w0 - __w1)
                                                + (__freq - __w2) / (__w0 - __w2)
                                                + (__freq - __w3) / (__w0 - __w3))
                            __I1 = 1.0 / 3.0 * (__freq - __w0) / (__w1 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w0) / (__w2 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w1 <= __freq < __w2:
                            __g = 3 / (__w3 - __w0) \
                                  * ((__freq - __w1) * (__freq - __w3) / ((__w2 - __w1) * (__w1 - __w3))
                                     + (__freq - __w0) * (__freq - __w2) / ((__w2 - __w0) * (__w1 - __w2)))

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3) \
                                   + (__freq - __w2) / (__w0 - __w2) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I1 = 1.0 / 3.0 * (__freq - __w2) / (__w1 - __w2) \
                                   + (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w1) / (__w2 - __w1) \
                                   + (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0) \
                                   + (__freq - __w1) / (__w3 - __w1) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w2 <= __freq < __w3:
                            __g = 3 * (__w3 - __freq) ** 2 / (__w3 - __w0) / (__w3 - __w1) / (__w3 - __w2)

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3)
                            __I1 = 1.0 / 3.0 * (__freq - __w3) / (__w1 - __w3)
                            __I2 = 1.0 / 3.0 * (__freq - __w3) / (__w2 - __w3)
                            __I3 = 1.0 / 3.0 * ((__freq - __w0) / (__w3 - __w0)
                                                + (__freq - __w1) / (__w3 - __w1)
                                                + (__freq - __w2) / (__w3 - __w2))

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        else:
                            pass

                        # tetrahedron 6: 0-2-6-7
                        # for ind_mode in range(v_q.shape[2]):
                        __w0 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ind_2nd_kpt * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F0 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ind_2nd_kpt * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w1 = _w_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F1 = abs(_v_q[ind_1st_kpt * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w2 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    ind_3rd_kpt]
                        __F2 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        ind_3rd_kpt, :]) ** 2
                        __w3 = _w_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                    ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                    (ind_3rd_kpt + 1) % num_3rd_kpt]
                        __F3 = abs(_v_q[((ind_1st_kpt + 1) % num_1st_kpt) * num_2nd_kpt * num_3rd_kpt +
                                        ((ind_2nd_kpt + 1) % num_2nd_kpt) * num_3rd_kpt +
                                        (ind_3rd_kpt + 1) % num_3rd_kpt, :]) ** 2

                        __w0, __w1, __w2, __w3 = np.asfarray([__w0,
                                                              __w1,
                                                              __w2,
                                                              __w3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3]))]
                        __F0, __F1, __F2, __F3 = np.asfarray([__F0,
                                                              __F1,
                                                              __F2,
                                                              __F3])[np.argsort(np.asfarray([__w0,
                                                                                             __w1,
                                                                                             __w2,
                                                                                             __w3])), :]

                        if __w0 <= __freq < __w1:
                            __g = 3 * (__freq - __w0) ** 2 / (__w1 - __w0) / (__w2 - __w0) / (__w3 - __w0)

                            __I0 = 1.0 / 3.0 * ((__freq - __w1) / (__w0 - __w1)
                                                + (__freq - __w2) / (__w0 - __w2)
                                                + (__freq - __w3) / (__w0 - __w3))
                            __I1 = 1.0 / 3.0 * (__freq - __w0) / (__w1 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w0) / (__w2 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w1 <= __freq < __w2:
                            __g = 3 / (__w3 - __w0) \
                                  * ((__freq - __w1) * (__freq - __w3) / ((__w2 - __w1) * (__w1 - __w3))
                                     + (__freq - __w0) * (__freq - __w2) / ((__w2 - __w0) * (__w1 - __w2)))

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3) \
                                   + (__freq - __w2) / (__w0 - __w2) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I1 = 1.0 / 3.0 * (__freq - __w2) / (__w1 - __w2) \
                                   + (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)
                            __I2 = 1.0 / 3.0 * (__freq - __w1) / (__w2 - __w1) \
                                   + (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w0) / (__w2 - __w0) \
                                   * (__freq - __w2) / (__w1 - __w2) / __g / (__w3 - __w0)
                            __I3 = 1.0 / 3.0 * (__freq - __w0) / (__w3 - __w0) \
                                   + (__freq - __w1) / (__w3 - __w1) \
                                   * (__freq - __w3) / (__w1 - __w3) \
                                   * (__freq - __w1) / (__w2 - __w1) / __g / (__w3 - __w0)

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        elif __w2 <= __freq < __w3:
                            __g = 3 * (__w3 - __freq) ** 2 / (__w3 - __w0) / (__w3 - __w1) / (__w3 - __w2)

                            __I0 = 1.0 / 3.0 * (__freq - __w3) / (__w0 - __w3)
                            __I1 = 1.0 / 3.0 * (__freq - __w3) / (__w1 - __w3)
                            __I2 = 1.0 / 3.0 * (__freq - __w3) / (__w2 - __w3)
                            __I3 = 1.0 / 3.0 * ((__freq - __w0) / (__w3 - __w0)
                                                + (__freq - __w1) / (__w3 - __w1)
                                                + (__freq - __w2) / (__w3 - __w2))

                            _pdos[:, ind_x] += volume_ith / volume_bz * __g \
                                               * (__I0 * __F0 + __I1 * __F1 + __I2 * __F2 + __I3 * __F3)

                        else:
                            pass

    return _pdos
