function [weq, beq] = congctrl_equil(N, rho, c);

% Compute out the equilibrium values
alpha = (1/rho)^2 / (2 * N^2);
brho = (-2*3^(1/3)*alpha + 2^(1/3) * ...
  (9*alpha.^2 + sqrt(3)*sqrt(alpha.^3.*(4 + 27*alpha))).^(2/3)) ./ ...
  (6^(2/3)*alpha.*(9*alpha.^2 + sqrt(3).*sqrt(alpha.^3.*(4 + 27*alpha))).^(1/3));
beq = brho / rho;
weq = beq/N;
