function [IpQ,WkQ,WpQ,VpQ,FtQ,Ft0Q,OpQ] = ...
  meqintQ(L,F0,F1, rBt , ag, Fx, Opy, sIp)
%MEQINTQ Computes integral quantities on pQ grid.
% [IpQ,WkQ,WpQ,VpQ,FtQ,Ft0Q,OpQ] = meqintQ(L,FA,FB, rBt , ag, Fx, Opy, sIp)
% see MEQINT
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

FQ = L.pQ.^2;
nB = find(F0~=F1,1,'last');

[IpQ,WkQ,WpQ,VpQ,FtQ,Ft0Q,OpQ] = deal(zeros(L.nQ,L.nD));

if all(ag==0), return; end
for jj = 2:L.nQ
  % Extract the flux inside FQ(jj) flux surface
  Opynew = Opy;
  for iD = 1:nB
    Opynew(Opy == iD & (sIp*(Fx(2:end-1,2:end-1)-F0(iD)-(F1(iD)-F0(iD))*FQ(jj)) <= 0)) = 0;
  end
  % mask TpDg,ITpDg by new Opy
  [~,TpDg,ITpDg] = L.bfct(1,L.bfp, Fx, F0, F1,Opynew,L.ry,L.iry);
  
  IpQ(jj,:) =  TpDg *ag ;
  x         = ITpDg.*ag';
  WkQ(jj,:) = x*L.fPg*1.5;
  WtQ       = x*L.fTg    ;
  FtQ(jj,:) = 2e-7*WtQ/rBt;
  for iD = 1:nB
    [WpQ(jj,iD),Ft0Q(jj,iD),VpQ(jj,iD)] = vizrmex(Fx,Opynew==iD,L.ry,L.iry,rBt,L.drx,L.dzx);
    OpQ(jj,iD) = sum(Opynew(:) == iD);
  end
end

end
