internal function to calculate pressure along Hugoniot

    :param rho: density in g/cm^3
    :param rho0: density at 1 bar in g/cm^3
    :param c0: velocity at 1 bar in km/s
    :param s: slope of the velocity change
    :param gamma0: Gruneisen parameter at 1 bar
    :param q: logarithmic derivative of Gruneisen parameter
    :param theta0: Debye temperature in K
    :param n: number of elements in a chemical formula
    :param mass: molar mass in gram
    :param three_r: 3 times gas constant.
        Jamieson modified this value to compensate for mismatches
    :param t_ref: reference temperature, 300 K
    :param c_v: heat capacity, see Jamieson 1983 for detail
    :return: temperature along hugoniot

def hugoniot_t_single(rho, rho0, c0, s, gamma0, q, theta0, n, mass,
                      three_r=3. * constants.R, t_ref=300., c_v=0.):
    """
    internal function to calculate pressure along Hugoniot

    :param rho: density in g/cm^3
    :param rho0: density at 1 bar in g/cm^3
    :param c0: velocity at 1 bar in km/s
    :param s: slope of the velocity change
    :param gamma0: Gruneisen parameter at 1 bar
    :param q: logarithmic derivative of Gruneisen parameter
    :param theta0: Debye temperature in K
    :param n: number of elements in a chemical formula
    :param mass: molar mass in gram
    :param three_r: 3 times gas constant.
        Jamieson modified this value to compensate for mismatches
    :param t_ref: reference temperature, 300 K
    :param c_v: heat capacity, see Jamieson 1983 for detail
    :return: temperature along hugoniot
    """
    eta = 1. - rho0 / rho
    if eta == 0.0:
        return 300.
    threenk = three_r / mass * n  # [J/mol/K] / [g/mol] = [J/g/K]
    k = [rho0, c0, s, gamma0, q, theta0 / 1.e3]
    t_h = odeint(_dT_h_delta, t_ref / 1.e3, [0., eta],
                 args=(k, threenk, c_v), full_output=1)
    temp_h = np.squeeze(t_h[0][1])
    return temp_h * 1.e3