%BFSP  B-spline basis functions
%
% Basis-functions are constructed using order 5 B-splines (degree 4) for
% the integral functions and their exact derivatives for p' and TT'. The
% knot distribution is so far the same for basis functions associated to p'
% or TT'. The parameters P must be obtained by calling the BFPSP function.
%
% So far only mode 0-3,5 have been fully implemented. In particular no
% regularisation has been implemented yet.
%
% See also: BFPSP
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function [y1,y2,y3,y4] = bfsp(mode,P,F,F0,F1,q1,q2,q3)
 if nargin>4, assert(numel(F0) <= 1 && numel(F1) <= 1,'bfsp:multiD',...
     'To use bfsp with multiple domains, please use bfgenD');end
 nP = size(P.M,2); nT = nP;
 ng = nP + nT;
 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,ng);
  case 1 % [TYG,TPDG,ITPDG  ] = BF(1,PAR, F,F0,F1,OY,RY,IRY )
   q1 = logical(q1); % Convert Opy to logical
   FN = reshape((F(2:end-1,2:end-1)-F0)/(F1-F0),[],1);
   y1 = bspsum(P.taug,P.M,FN,1);
   y3 = bspsum(P.taug,P.M,FN,0)*(F1-F0);
   % Combine with r or 1/r
   y1 = [y1.*reshape(q1.*q2(:).',[],1),...
         y1.*reshape(q1.*q3(:).',[],1)];
   y3 = [y3.*reshape(q1.*q2(:).',[],1),...
         y3.*reshape(q1.*q3(:).',[],1)];
   % Compute the integrals   
   y2 = sum(y1,1);
   y3 = sum(y3,1);
  case 2 % [GDQG,IGDQG      ] = BF(2,PAR,FN,F0,F1           )
   y1 = bspsum(P.taug,P.M,F,1);
   y2 = bspsum(P.taug,P.M,F,0);
   y1 = [y1,y1];
   y2 = [y2,y2];
  case 3 % [APP,ATTP,AP,AHQT] = BF(3,PAR, A,F0,F1,FP,FT,IDS )
   y1 = 0.5/pi        *q3 * F.*q1;
   y2 = 2e-7          *q3 * F.*q2;
   y3 = 0.5/pi*(F1-F0)*q3 * F.*q1;
   y4 = 2e-7  *(F1-F0)*q3 * F.*q2;
         % [TDG,TGY         ] = BF(4,PAR, F,F0,F1,O ,KD, FD )
  case 5 % [G0G,IG0G        ] = BF(5,PAR, ~,F0,F1           )
   y1 = bspsum(P.taug,P.M,0,1);
   y2 = bspsum(P.taug,P.M,0,0)*(F1-F0);
   y1 = [y1,y1];
   y2 = [y2,y2];
  case 6 % [QQG,XQ          ] = BF(6,PAR, ~,F0,F1,R0,IR0,IDS)
   y1 = zeros(0,ng);
   y2 = zeros(0, 1);
  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

