function ag=meqfitprof(L,FN,FA,FB,PpQ,TTpQ)
%function ag = meqfitprof(L,FN,FA,FB,PpQ,TTpQ)
% fit PpQ,TTpQ profiles with ag coefficients
% FN: normalized flux radial coordinate
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

nt = size(FA,2);

if ~any(PpQ(:)) && ~any(TTpQ(:))
  ag = zeros(L.ng,nt);
  return
end

if all(size(L.bfp)>1)
  assert(~any(any(diff(L.bfp,[],2))),'time-varying bfp non supported');
end

% find ag that fits PpQ,TTpQ profiles
ag  = ones(L.ng,nt);
ago = ones(L.ng, 1); % ones just to get scaling
for kt=1:nt
  [F0,F1] = meqpdomlim(FA(:,kt),FB(:,kt),L.nD);
  [gN]       = L.bfct(2,L.bfp,FN,F0,F1);
  % get scaling between ag=1 and p',TT' equation scaling
  [aPp,aTTp] = L.bfct(3,L.bfp,ago,F0,F1,L.fPg,L.fTg,L.idsx);
  % fit ag coefficients to desired p',TT'
  ag(:,kt) = [gN.*aPp';gN.*aTTp']\[PpQ(:,kt);TTpQ(:,kt)];
end