%BFEF  Constant and polynomial base functions for p',TT'
% BFEF uses P=[nP nT] and produces base functions for p',TT' as the first
% nP<=4,nT<=4 functions in {1, F-F1, (F-F1)*(F-F0),
% (F-F1)*(F-F0)*(2*F-F1-F0)}. See also BFHELP.
%
% For details see: [MEQ-redbook] 
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function [y1,y2,y3,y4] = bfef(mode,P,F,F0,F1,q1,q2,q3)
 if nargin>4, assert(numel(F0) <= 1 && numel(F1) <= 1,'bfef:multiD',...
     'To use bfef with multiple domains, please use bfgenD');end
 nP = P(1); nT = P(2); ng = nP + nT;
 kg = true(1,ng);
 if nP > 0, kg(1)    = false; end
 if nT > 0, kg(1+nP) = false; end
 switch mode
  case 0 % [FP,FT,TDG       ] = BF(0,PAR                    )
   y1 = [ ones(nP,1) ; zeros(nT,1)];
   y2 = [zeros(nP,1) ;  ones(nT,1)];
   y3 = ones(1,nP+nT);
  case 1 % [TYG,TPDG,ITPDG  ] = BF(1,PAR, F,F0,F1,OY,RY,IRY )
   y1      =      zeros(numel(q1),ng);
   [y2,y3] = deal(zeros(1        ,ng));
   if nP > 0
    y1(:,1) = reshape(logical(q1).*q2(:)',[],1); % Tyg(Opy,1) = ry
    y2(1,1) = sum(y1(:,1));
    y3(1,1) = sum(reshape(F(2:end-1,2:end-1)-F1,[],1).*y1(:,1));
   end
   if nT > 0
    y1(:,1+nP) = reshape(logical(q1).*q3(:)',[],1); % Tyg(Opy,1) = 1/ry
    y2(1,1+nP) = sum(y1(:,1+nP));
    y3(1,1+nP) = sum(reshape(F(2:end-1,2:end-1)-F1,[],1).*y1(:,1+nP));
   end
   [y1(:,kg),y2(:,kg),y3(:,kg)] = bfabmex(mode,max(P-1,0),F,F0,F1,q1,q2,q3);
  case 2 % [GQDG,IGQDG      ] = BF(2,PAR,FN,F0,F1           )
   [y1,y2] = deal(ones(numel(F),ng));
   if nP > 0, y2(:,1   ) = (F-1); end
   if nT > 0, y2(:,1+nP) = (F-1); end
   [y1(:,kg),y2(:,kg)] = bfabmex(mode,max(P-1,0),F,F0,F1);
  case 3 % [APP,ATTP,AP,AHQT] = BF(3,PAR, A,F0,F1,FP,FT,IDS )
   [y1,y2,y3,y4] = deal(zeros(ng,1));
   y1(~kg) = 0.5/pi        *q3 * F(~kg).*q1(~kg);
   y2(~kg) = 2e-7          *q3 * F(~kg).*q2(~kg);
   y3(~kg) = 0.5/pi*(F1-F0)*q3 * F(~kg).*q1(~kg);
   y4(~kg) = 2e-7  *(F1-F0)*q3 * F(~kg).*q2(~kg);
   [y1(kg),y2(kg),y3(kg),y4(kg)] = bfabmex(mode,max(P-1,0),F(kg),F0,F1,q1(kg),q2(kg),q3);
         % [TDG,TGY         ] = BF(4,PAR, F,F0,F1,O ,KD, FD )
  case 5 % [G0G,IG0G        ] = BF(5,PAR, ~,F0,F1           )
   [y1,y2] = deal(zeros(1,ng));
   y1(~kg) = 1;
   y2(~kg) = F0-F1;
   [y1(:,kg),y2(:,kg)] = bfabmex(mode,max(P-1,0),[],F0,F1);
  case 6 % [QQG,XQ          ] = BF(6,PAR, ~,F0,F1,R0,IR0,IDS)
   [y1(:,kg),y2] = bfabmex(mode,max(P-1,0),F,F0,F1,q1,q2,q3);
   % Adding constraint to minimise jphi at the boundary, while minimising
   % the difference between the p',TT' contributions therein.
   nq = size(y1,1)+2;
   if nP > 0, y1(nq+(-1:0),   1) = q1*q3*[ 1; 1]; y2(nq+(-1:0),1) = 0; end % rA
   if nT > 0, y1(nq+(-1:0),nP+1) = q2*q3*[ 1;-1]; y2(nq+(-1:0),1) = 0; end % 1/rA
  case 7 % [QCG,XC          ] = BF(7,PAR, ~,F0,F1           )
   y1 = zeros(0,ng);
   y2 = zeros(0, 1);
  otherwise
    error('MEQ:BFEF:CaseNotImplemented','Mode %d not implemented for %s',mode,mfilename);
 end
end

