function [bfp,tbfp] = bfpfbt(P)
% [bfp,tbfp] = bfpfbt(P)
% Computes basis function parameters for BFFBT
%
% Parameters: [TO BE POPULATED]
% .ppfac
% .ppal
% .ell1
% .ttfac
% .omega
% .emm1 
%
% Returns: bfp = [ppfac,ppal,ell1,P.ttfac1.*rh.^2,omega,emm1]
%             where rh radius is the center of the box bounding all points on the
%             plasma LCFS, as determined by P.gpr and P.gpb. In the absence
%             of P.gpr,P.gpb, rh=P.r0;
%          tbfp: corresponding equilibrium time
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%% Compute center of bounding box
if isfield(P,'gpr') && ~isempty(P.gpr)
  rH = P.gpr;
  rH(isnan(P.gpr) | isnan(P.gpz) | ~isfinite(P.gpb) | P.gpb == 0) = NaN; % Eliminate points outside of the plasma domain
  rh = 0.5*(min(rH,[],1) + max(rH,[],1));
else
  rh = P.r0; % in case no gpr specified
end

% defaults
D.ppfac  = 1;
D.ppal   = 1;
D.ell1   = 1;
D.ttfac1 = 1;
D.omega  = 0;
D.emm1   = 1.8;
D.t      = 0;

for field=fieldnames(D)'
  if ~isfield(P,field{1})
    P.(field{1}) = D.(field{1});
  end
end

%% Base function parameters
bfp = [P.ppfac;
       P.ppal;
       P.ell1;
       P.ttfac1.*rh.^2;
       P.omega;
       P.emm1];

tbfp = P.t;
end