code: def hugoniot_t_single(rho, rho0, c0, s, gamma0, q, theta0, n, mass,
                        three_r=3. * constants.R, t_ref=300., c_v=0.):
      
      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
docs: 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