function [integr] = int_SjdotB(quantity,g,model,varargin)
%function [integral] = int_SjdotB(quantity,g,model,indices)
% integrate flux-sufrace-averaged current densities over plasma surface
%
% Ipl = F/2pi int V'/F^2 <jdotB> drho
% quantity is a <jdotB> current denstity in (TA/m^2)
%
% quantity: quantity to be integrated on gauss grid
% 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

Vpgauss = model.geom.LamGgauss*g(model.geom.ind_Vp);
Fgauss  = model.geom.LamGgauss*g(model.geom.ind_F);
F       = model.geom.LamG*g(model.geom.ind_F);

if numel(indices) == 1 && (indices == model.rgrid.nrho)
    integr = F(end)/(2*pi)*sum(model.rgrid.wgauss.*(Vpgauss./Fgauss.^2).*quantity);
else
    cumsumall = cumsum(model.rgrid.wgauss.*(Vpgauss./Fgauss.^2).*quantity);
    integrall = 1/(2*pi)* F.*[0;cumsumall(model.rgrid.ngauss:model.rgrid.ngauss:end)];
    integr = integrall(indices);
end

return
