code: def hugoniot_p_nlin(rho, rho0, a, b, c):
      
      eta = 1. - (rho0 / rho)
      Up = np.zeros_like(eta)
      if isuncertainties([rho, rho0, a, b, c]):
          Up[eta != 0.] = ((b * eta - 1.) + unp.sqrt(
              np.power((1. - b * eta), 2.) - 4. * np.power(eta, 2.) * a * c)) /\
              (-2. * eta * c)
      else:
          Up[eta != 0.] = ((b * eta - 1.) + np.sqrt(
              np.power((1. - b * eta), 2.) - 4. * np.power(eta, 2.) * a * c)) /\
              (-2. * eta * c)
      Us = a + Up * b + Up * Up * c
      Ph = rho0 * Up * Us
      return Ph
docs: calculate pressure along a Hugoniot throug nonlinear equations
      presented in Jameison 1982
  
      :param rho: density in g/cm^3
      :param rho0: density at 1 bar in g/cm^3
      :param a: prefactor for nonlinear fit of Hugoniot data
      :param b: prefactor for nonlinear fit of Hugoniot data
      :param c: prefactor for nonlinear fit of Hugoniot data
      :return: pressure along Hugoniot in GPa