docs: 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