function [beta,dbeta_dx] = eval_beta(x,g,v,model)
% beta = 2/3* Wth / (B0^2/2*mu0) * V))
% #codegen

if nargout>1
    der = true; % compute derivative
    dbeta_dx = zeros(1,numel(x)); % init
else
    der = false;
end
geop = geometry_profiles(g,0*g,model); % cheating
V = int_Vtot(ones(size(model.rgrid.rhogauss)),geop,model); % integral (dV)
[Wth,dWth_dx] = eval_Wth(x,g,v,model,false);

beta_constant = 2/3 ./ (model.equi.B0.^2 / (2*4e-7*pi) .* V);
beta = beta_constant.*Wth ; % percent

if der
    dbeta_dx = bsxfun(@times,beta_constant,dWth_dx);
end

return
