function y = firmfoc(mu, z, D, p)   


q       = p.sigma^(p.sigma/p.epsi)*(1 - 1./mu).^(p.sigma/p.epsi); 


Uq      = (p.sigma - 1)/p.sigma*exp((1 - q.^(p.epsi/p.sigma))/p.epsi); 


y       = Uq.*q - D*mu.*(q./z)/(1 + p.xis);   % this is markup over cost including the subsidy
                   