internal function for calculation of temperature along a Hugoniot

    :param T_in_kK: temperature in kK scale, see Jamieson for detail
    :param eta: = 1 - rho0/rho
    :param k: = [rho0, c0, s, gamma0, q, theta0]
    :param threenk: see the definition in Jamieson 1983,
        it is a correction term mostly for Jamieson gold scale
    :param c_v: manual input of Cv value,
        if 0 calculated through Debye function
    :return: eta derivative of temperature

def _dT_h_delta(T_in_kK, eta, k, threenk, c_v):
    """
    internal function for calculation of temperature along a Hugoniot

    :param T_in_kK: temperature in kK scale, see Jamieson for detail
    :param eta: = 1 - rho0/rho
    :param k: = [rho0, c0, s, gamma0, q, theta0]
    :param threenk: see the definition in Jamieson 1983,
        it is a correction term mostly for Jamieson gold scale
    :param c_v: manual input of Cv value,
        if 0 calculated through Debye function
    :return: eta derivative of temperature
    """
    rho0 = k[0]  # g/m^3
    gamma0 = k[3]  # no unit
    q = k[4]  # no unit
    theta0_in_kK = k[5]  # K, see Jamieson 1983 for detail
    rho = rho0 / (1. - eta)
    c0 = k[1]  # km/s
    s = k[2]  # no unit
    dPhdelta_H = rho0 * c0 * c0 * (1. + s * eta) / \
        np.power((1. - s * eta), 3.)
    # [g/cm^3][km/s]^2 = 1e9[kg m^2/s^2] = [GPa]
    Ph = hugoniot_p(rho, rho0, c0, s)  # in [GPa]
    # calculate Cv
    gamma = gamma0 * np.power((1. - eta), q)
    theta_in_kK = theta0_in_kK * np.exp((gamma0 - gamma) / q)
    x = theta_in_kK / T_in_kK
    debye3 = debye_E(x)
    if c_v == 0.:
        c_v = threenk * (4. * debye3 - 3. * x / (np.exp(x) - 1.))  # [J/g/K]
    # calculate dYdX
    dYdX = (gamma / (1. - eta) * T_in_kK) + (dPhdelta_H * eta - Ph) / \
        (2. * c_v * rho0)
    # print('dYdX', dYdX)
    return dYdX