function a=n_h_fun(n_h0,alpha,w_r,g,gamma,sigma,phi)
% anonymous function
 myfun = @(a) n_h-alpha*(w_r* ((n_h-1)/(g*(1-gamma)))^(-sigma))^(1/phi);
a = fzero(myfun, n_h0)
end
