function ag = fbtfit(L,rA,dr2FA,dz2FA,drzFA,fPg,fTg,Tpg,~,gAg,IgAg,rBt,Ip,~,qA)
% function ag = fbtfit(L,rA,dr2FA,dz2FA,drzFA,fPg,fTg,Tpg,ITpg,gAg,IgAg,rBt,Ip,bp,qA)
% FBT legacy fit
% Using 3 basis functions (1 p', 2 TT'), where the first 2 have the same coefficient
%   Solving for 2 constraints, Ip and qA
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
assert(~isempty(rA),'rA can not be empty for fbtfit, use something else for doublets');

IgAg = 2e-7*L.idsx*IgAg;
TyAg = gAg.*(rA*fPg+fTg/rA).'; % TyAg: Axis current IA=(jphi.dR.dZ) contribution from each basis function

% qA constraints starts from the expression for qA, (eq88 of 2015 paper), which
% depends on TA and IA. Both these depend in turn on ag.
% NOTE: for a coarse grid, this can cause a difference between the
% target value and the achieved value of qA.
cqA = qA*4e-7*pi*rA^2*rBt*sqrt(dr2FA*dz2FA-drzFA^2)/(abs(dr2FA+dz2FA)*L.dsx);

AA=[Tpg ; cqA*TyAg - IgAg.*fTg.'; 1 -1 0 ];
assert(rank(AA)==3,'Coefficient fitting matrix is rank deficient, insufficient freedom in basis functions?')

ag = AA \ [Ip ; rBt^2 ; 0]; % Row 1: Ip constraint. Row 2: qA constraint. Row 3: ag(1) = ag(2).
end
