function [res,AA,bb] = meqagcon(L,LX,LY,F0,F1,rA,dr2FA,dz2FA,drzFA,IpD,ag,Fx,Opy,fPg,fTg,TpDg,ITpDg)
% function res = meqagcon(L,LX,LY,F0,F1,rA,gA,ag,Fx,Opy,fPg,fTg,TpDg,ITpDg)
% Compute residuals for basis function coefficient constraints in FGE
% Also return AA,bb such that res = AA*ag - bb;
% Evaluates residual vector based on constraint specifications in L.agconc
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

nB = numel(F1);
irr = 0; % equation row index
AA=zeros(L.ng,L.ng); bb=zeros(L.ng,1);

if all(IpD==0)
  res = zeros(L.ng,1); return;
end

for iC = 1:L.nC % for all constraints
  agconc = L.agconc{iC,1}; iD = L.agconc{iC,2};
  ii = L.agconc{iC,4}; % index for this constraint
  if iD<=nB
    % call constraint function
    [A,b] = agconc(L,LX,LY,...
      F0,F1,rA,dr2FA,dz2FA,drzFA,...
      IpD,ag,Fx,Opy,fPg,fTg,TpDg,ITpDg,ii);
  else
    % for non-active constraints, keep previous ag
    A = zeros(1,L.ng);
    A(:,ii) = -1; b = LY.ag(ii);
  end
  
  if L.P.debug>4
    fprintf('evaluated constraint %d: %s, size=%d, ires=%d, iD=%d, ii = %d, res=%+5.3e\n',...
      iC,func2str(agconc),size(A,1),irr(end),iD,ii,A*ag-b);
  end
  
  % indices
  irr  = irr(end) + (1:(size(A,1)));
  
  % build linear system
  AA(irr,:) = A;
  bb(irr,1) = b;
end
assert(irr(end) == (L.ng),...
  'Number of residuals (%d) does not match total number of basis functions (%d)',irr(end),L.ng)

% final calculation of residual
res = AA*ag - bb;
end
