function [integr] = int_SjdotB_edge(quantity,geop,model)
%function [integral] = int_SjdotB(quantity,geop,model)
% integrate flux-sufrace-averaged current densities over entire 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
% geop: geometry

% #codegen
integr = geop.Fb/(2*pi)*sum(model.rgrid.wgauss.*(geop.Vp./geop.F.^2).*quantity);

return
