function [betaN,dbetaN_dx] = eval_betaN(x,g,v,model)
% #codegen

if nargout>1
    der = true; % compute derivative
    dbetaN_dx = zeros(1,numel(x)); % init
else
    der = false;
end

% out.pav = 2/3*out.Wth ./ out.Volume(end,:); % volume averaged pressure
% out.beta = out.pav/(model.equi.B0.^2 / (2*4e-7*pi)); % percent
% out.betaN = 100*1e6*out.beta./(out.Ip(end,:)./(model.equi.B0*model.equi.epsilon*model.equi.R0));

geop = geometry_profiles(g,0*g,model);
V=int_Vtot(ones(size(model.rgrid.rhogauss)),geop,model);

[Ip,dIp_dx] = eval_Ipl(x,g,v,model,false);
[Wth,dWth_dx] = eval_Wth(x,g,v,model,false);
betaN_constant=100*2/3*1e6*(model.equi.epsilon*model.equi.R0*2*4e-7*pi)/(V*model.equi.B0);

betaN = betaN_constant*Wth/Ip(end);

if der
    dbetaN_dx = bsxfun(@times,betaN_constant/Ip(end),dWth_dx)-bsxfun(@times,betaN,dIp_dx(end,:)/1e6);
end

return
