function [integr] = int_V(quantity,g,model,varargin)
%function [integral] = int_V(quantity,g,model,indices)
% integrate quantity over plasma volume
% int(quantity dV/drho) drho
% quantity: quantity on rhogauss
% integr: integral on rhotor
%
% indices: optional choice of indices, if none given, all are returned

% #codegen

if numel(varargin)>0
    indices = varargin{1};
else
    indices = 1:model.rgrid.nrho; % return all
end
%%
Vpgauss = eval_Vp([],g,[],model,true);
if numel(indices) == 1 && (indices == model.rgrid.nrho)
    integr = sum(bsxfun(@times,model.rgrid.wgauss.*Vpgauss,quantity));
else
    cumsumall = cumsum(bsxfun(@times,model.rgrid.wgauss.*Vpgauss,quantity));
    integrall = [0;cumsumall(model.rgrid.ngauss:model.rgrid.ngauss:end)];
    integr = integrall(indices);
end
%%
return
