function [PpQ,TTpQ,PQ,TQ,iTQ,PpQg,TTpQg] = meqprof(fPg,fTg,ag,FN,FA,FB,rBt,bfct,bfp,idsx,smalldia)
%MEQPROF Computes profiles from basis functions and fitting parameters
% function [PpQ,TTpQ,PQ,TQ,iTQ,PpQg,TTpQg] = meqprof(fPg,fTg,ag,FN,FA,FB,rBt,bfct,bfp,idsx,smalldia)
% -
% INPUTS
% fPg, fTg : basis function assignment flags (see BFHELP)
% ag       : basis function coefficients
% FN       : Normalized flux on Q grid = L.pQ^2;
% FA, FB   : Axis, boundary flux
% rBt      : Vacuum R*Bt
% bfct, bfp: basis function handle and parameters (see BFHELP)
% idsx     : 1/(drx*dzx) = 1/(grid element surface)
% smalldia : use small diamagnetism approximation  for TQ computation
% - 
% OUTPUTS
% PpQ,TTpQ  :  pprime and TTprime. 
% PQ,TQ     :  Their primitives, assuming P(1)=0, T(1)=rBt
% iTQ       :  1/TQ,
% PpQg,TTpQg:  Contributions of each basis function to PpQ and TTpQ.
% The grid Q has associated L.pQ: the value of rhopsi=sqrt((Fx-FA)/(FB-FA)).
% See also BFHELP. For details, see [MEQ-redbook] documentation.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if nargin<11
  smalldia=false; % use small diamagnetic approximation optionally
end
nQ = numel(FN);
nD = numel(FA);

% Computes profiles from basis functions
[ gQg, IgQg          ] = bfct(2,bfp,FN,FA,FB);
[aPpg,aTTpg,aPg,ahqTg] = bfct(3,bfp,ag,FA,FB,fPg,fTg,idsx);
PpQg  = reshape(reshape(gQg,nQ*nD,[]).* aPpg.',nQ,nD,[]);
TTpQg = reshape(reshape(gQg,nQ*nD,[]).*aTTpg.',nQ,nD,[]);
PpQ   = sum( PpQg,3);
TTpQ  = sum(TTpQg,3);
PQ    = reshape(reshape(IgQg,nQ*nD,[])*aPg  ,nQ,nD);
hqTQ  = reshape(reshape(IgQg,nQ*nD,[])*ahqTg,nQ,nD);

if nargout>3
  irBt = 1./rBt;
  if smalldia
    TQ =   rBt +      irBt*hqTQ; % small diamagnetic approximation
    iTQ = irBt - (irBt*irBt*irBt)*hqTQ;
  else
    TQ = sign(rBt)*sqrt(2*hqTQ+rBt*rBt);
    iTQ = 1./TQ;
  end
end
end
