function [g,g1,g2] = consfn(cf,cm,gpar);

g     = cf+cm;
g1    = ones(length(cf),1);
g2    = g1;

return
etaf  = gpar(1);
etam  = gpar(2);
rho   = gpar(3);
g     = (etaf*cf.^rho+etam*cm.^rho).^(1/rho);
g1    = etaf*cf.^(rho-1).*(etaf*cf.^rho+etam*cm.^rho).^(1/rho-1);
g2    = etam*cm.^(rho-1).*(etaf*cf.^rho+etam*cm.^rho).^(1/rho-1);



