function [integr] = int_Wpol(quantity,g,model,varargin)
%function [integral] = int_Wpol(quantity,g,model,indices)
% integrate poloidal field energy over plasma
%
% WI = 1/(8*pi^2*mu0)* int g2/V' *(dpsi/drho).^2 drho
%
% quantity: quantity to be integrated on gauss grid (usually dpsi/drho for
% this case)
% g: geometry
% model: RAPTOR model
% indices (optional): indices of integral to return (on rho grid)

% #codegen

if numel(varargin)>0
    indices = varargin{1};
    assert(all(indices>0),'index must be >1')
    assert(all(indices<=model.rgrid.nrho),'indices can not be larger than nrho')
else
    indices = 1:model.rgrid.nrho; % return all
end

g2gauss = eval_g2([],g,[],model,true);
Vpgauss = eval_Vp([],g,[],model,true);
mu0 = 4e-7*pi;

if numel(indices) == 1 && (indices == model.rgrid.nrho)
    integr = 1/(8*pi^2*mu0)*sum(bsxfun(@times,model.rgrid.wgauss.*(g2gauss./Vpgauss),quantity));
else
    cumsumall = cumsum(bsxfun(@times,model.rgrid.wgauss.*(g2gauss./Vpgauss),quantity));
    integrall = 1/(8*pi^2*mu0)*...
        [zeros(1,size(cumsumall,2));...
        cumsumall(model.rgrid.ngauss:model.rgrid.ngauss:end,:)];
    integr = integrall(indices,:);
end

return
