function [betap,dbetap_dx] = eval_betap(x,g,v,model)
% function [betap,dbetap_dx] = eval_betap(x,g,v,model)
% #codegen

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

[Ip,dIp_dx] = eval_Ipl(x,g,v,model,false);
[Wth,dWth_dx] = eval_Wth(x,g,v,model,false);

betap = 1e7/(pi*model.equi.R0)*2/3*Wth(end)./(Ip(end).^2);
if der
    dbetap_dx = bsxfun(@times,1e7/(pi*model.equi.R0)*2/3./(Ip(end)^2),dWth_dx(end,:)) + ...
        -2*bsxfun(@times,betap./Ip(end),dIp_dx(end,:));
end

return
