code: def hugoniot_rho_single(p, rho0, c0, s, min_strain=0.01):
      
      if p <= 1.e-5:
          return rho0
  
      def f_diff(rho):
          return hugoniot_p(rho, rho0, c0, s) - p
      rho = brenth(f_diff, rho0, rho0 / min_strain)
      return rho
docs: calculate density in g/cm^3 from a hugoniot curve
  
      :param p: pressure in GPa
      :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 min_strain: defining minimum v/v0 value to search volume for
      :return: density in g/cm^3