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