docs: def hugoniot_p(rho, rho0, c0, s):
      """
      calculate pressure along a 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
      :return: pressure in GPa
      """
      eta = 1. - (rho0 / rho)
      Ph = rho0 * c0 * c0 * eta / np.power((1. - s * eta), 2.)
      return Ph