function f = rho(z,plist)

   A = z(1);
   P = z(2);
   r = plist(:,1);
   c = plist(:,2);
   f = r - (2*A*P)./(P*(1-c) + A*(1+c));

  return
end

