function [y1,y2,y3,y4] = bffbt(mode,P,F,F0,F1,q1,q2,q3)
%BFFBT  FBT base functions for p',TT'
% BFFBT produces base functions pf*(pl+(1-pl)*FN)*(1-FN)^pe for p' and
% tf*(1-FN)^te and tf*(tl*(FN-1)-1)*FN*(1-FN)^te for TT'. 
% P = [pf pl pe tf tl te]. 
%
% See also BFHELP and BFPFBT
%
% For details see [MEQ-redbook]
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if nargin>4, assert(numel(F0) <= 1 && numel(F1) <= 1,'bffbt:multiD',...
    'To use bffbt with multiple domains, please use bfgenD');end

if nargin>4 && ~isempty(F1)
  if ~any(mode == [2,3]), F  = F-F1;end
  F0 = F0-F1;
end

% Check shape and size of parameters
assert(isvector(P) && numel(P) == 6,'Parameters for bffbt must be in the form of a 6 element vector');
pf=P(1);pl=P(2);pe=P(3);tf=P(4);tl=P(5);te=P(6); 


fP = [1;0;0];
fT = [0;1;1];
switch mode
 case 0 % [FP,FT,TDG       ] = BF(0,PAR                    )
  y1 = fP;
  y2 = fT;
  y3 = ones(1,3);
 case 1 % [TYG,TPDG,ITPDG  ] = BF(1,PAR, F,F0,F1,OY,RY,IRY )
  FN = 1-1/F0*F(2:end-1,2:end-1);
  FN(~q1) = 1;
  [y1,y2] = bffbt(2,P,FN(:));
  y1 = [reshape(reshape(y1(:,  1),[size(F)-2 1]).*reshape(    q2,1,[]),[],1) ...
        reshape(reshape(y1(:,2:3),[size(F)-2 2]).*reshape(    q3,1,[]),[],2)];
  y2 = [reshape(reshape(y2(:,  1),[size(F)-2 1]).*reshape(-F0*q2,1,[]),[],1) ...
        reshape(reshape(y2(:,2:3),[size(F)-2 2]).*reshape(-F0*q3,1,[]),[],2)];
  y3 = sum(y2,1);
  y2 = sum(y1,1);
 case 2 % [GDQG,IGDQG      ] = BF(2,PAR,FN,F0,F1           )
  F  = reshape(F,[],1);
  f  = 1-F;
  x1 =      f.^pe;
  x2 = tf*f.^te;
  y1 = [(pf*pl+pf*(1-pl)*F).*x1 x2 (tl*F-(1+tl)).*F.*x2];
  x1 = x1.*f;
  x2 = x2.*f;
  y2 = [((-pf*(pl+1/(pe+1))/(pe+2))+pf*((pl-1)/(pe+2))*F).*x1 ...
        (-1/(te+1))*x2 ...
        ((tl/((te+2)*(te+3))+1/((te+1)*(te+2)))+((tl*(te+1)/((te+2)*(te+3))+1/(te+2))-(tl/(te+3))*F).*F).*x2];        
 case 3 % [APP,ATTP,AP,AHQT] = BF(3,PAR, A,F0,F1,FP,FT,IDS )
  cP = q3/(2*pi);
  cT = q3* 2e-7 ;
  y1 = cP*fP(:).*F;
  y2 = cT*fT(:).*F;
  y3 = -F0*y1;
  y4 = -F0*y2;
        % [TDG,TGY         ] = BF(4,PAR, F,F0,F1,O ,KD, FD )
 case 5 % [G0G,IG0G        ] = BF(5,PAR, ~,F0,F1           )
  y1 = [pf*pl tf 0];
  y2 = F0*[pf*(pl+1/(pe+1))/(pe+2) tf/(1+te) -tf*(tl/((te+2)*(te+3))+1/((te+1)*(te+2)))];
 case 6 % [QQG,XQ          ] = BF(6,PAR, ~,F0,F1,R0,IR0,IDS)
  y1 = zeros(0,3);
  y2 = zeros(0,1);
 case 7 % [QCG,XC          ] = BF(7,PAR, ~,F0,F1           )
  y1 = zeros(0,3);
  y2 = zeros(0,1);
otherwise
  error('MEQ:BFFBT:CaseNotImplemented','Mode %d not implemented for %s',mode,mfilename);
end
end