function ag = meqfit3(L,rA,dr2FA,dz2FA,drzFA,fPg,fTg,Tpg,ITpg,gAg,IgAg,rBt,Ip,bp,qA)
% function ag = meqfit3(L,rA,dr2FA,dz2FA,drzFA,fPg,fTg,Tpg,ITpg,gAg,IgAg,rBt,Ip,bp,qA)
% find ag such that the profiles match Ip,qA,bp
% Parameter L.P.agfitp = true(3,1) (default) uses all 3 constraints
% set one to false to disable it (for less BFs)
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if isfield(L.P,'agfitp') && ~isempty(L.P.agfitp), agfitp = L.P.agfitp;
else, agfitp=true(3,1); end
ng = numel(fPg);
assert(sum(agfitp)==ng,'must have same number of basis functions as constraints for meqfit3');
assert(~agfitp(3) || isfinite(bp),'must have finite bp if used for meqfit3');

if agfitp(2)
  if ~isempty(rA) % regular call for axis domain
    assert(~agfitp(2) || isfinite(qA),'must have finite qA if used for meqfit3');
    
    IgAg = 2e-7*L.idsx*IgAg;
    
    % qA equation [(88) in 2015 paper]: jphi*rA^2*mu0*qA*gamma = TA
    % with gamma = sqrt(det(H))/tr(H) (H hessian of axis)
    % Multiply rBt on both sides and use T = rBt + F^2/2/rBt (small diamag approximation, with F^2/2 primitive of TT')
    % gives: caA*TyAg*ag = rBt.^2 + (IgAg*fTg.')*ag
    % NOTE: for a coarse grid, this can cause a difference between the
    % target value and the achieved value of qA.
    
    TyAg = gAg.*(rA*fPg+fTg/rA).'; % TyAg: Axis current (jphi.dR.dZ) contribution from each basis function
    cqA = qA*4e-7*pi*rA^2*rBt*sqrt(dr2FA*dz2FA-drzFA^2)/(abs(dr2FA+dz2FA)*L.dsx);
    AqA = cqA*TyAg - IgAg.*fTg.'; bqA = rBt^2;
  else
    % call for mantle domain - replace qA equation by constraint making all
    % fTg coefficients except the first one equal to 0.
    AqA = fTg'; i1=find(fTg>0,1,'first'); AqA(i1) = 0; bqA = 0;
  end
else
  AqA = zeros(0,ng); bqA = [];
end

if agfitp(1)
  % Ip equation Tpg*ag = Ip;
  AIp = 1e-6*Tpg; bIp = 1e-6*Ip;
else
  AIp = zeros(0,ng); bIp=[];
end

if agfitp(3)
  % bp equation betap = 2/3*Wk/WN with Wk = (ITpg.*fPg.')*ag
  WN  = 1e-7*pi*L.P.r0*Ip^2;
  Abp = 1e-3*ITpg.*fPg'; bbp = 1e-3*bp*WN;
else
  Abp = zeros(0,ng); bbp=[];
end

A = [ AIp; AqA; Abp];
b = [ bIp; bqA; bbp];
% A*ag=b with Row 1: Ip constraint. Row 2: qzero constraint. Row 3: 3/2*betap*WN=Wk.

assert(rank(A)==size(A,2),'Coefficient fitting matrix is rank deficient, insufficient freedom in basis functions?')
ag = A \ b; % fit selecting only active constraints
end

