function [Wpol,dWpol_dx] = eval_Wpol(x,g,~,model,gaussgrid)
% function [Wpol,dWpol_dx] = eval_Wpol(x,g,v,model)
% #codegen

if gaussgrid
    error('Wpol is an integral quantity, can not be evaluated on gauss grid')
end

if nargout>1
    der = true; % compute derivative
    dWpol_dx = zeros(model.rgrid.nrho,numel(x)); % init
else
    der = false;
end

bpol = (model.psi.Lampgauss*x(model.psi.xind,:));
Wpol = int_Wpol(bpol.^2,g,model); % Wi (magnetic energy)

if der
    dbpol_dpsi = model.psi.Lampgauss;
    dWpol_dpsi = int_Wpol(2*(bsxfun(@times,bpol,dbpol_dpsi)),g,model);
    dWpol_dx(:,model.psi.xind) = dWpol_dpsi;
end

return
