%BFGEND generic multi-domain base function set
% BFGEND assembles single-domain basis function sets into one multi-domain
% basis function set. The parameters P are defined as a cell array where
% each line is composed of the function handle for a single-domain basis
% function set, its parameters and finally a column vector indicating which
% domains it applies to (must be of size [nD,1]).
%
% Example: {@bf3pmex,0,[1;0;0];@bf3pmex,0,[0;1;0];@bf3pmex,1,[1;1;1]}
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function [y1,y2,y3,y4] = bfgenD(mode,P,F,F0,F1,q1,q2,q3)

nbf  = size(P,1);
bfct = P(:,1);
bfp  = P(:,2);
lD   = P(:,3);
if nbf, nD = size(lD{1},1);end

% Get number of basis functions
ngi = zeros(nbf,1);
for ibf = 1:nbf
  ngi(ibf) = numel(P{ibf,1}(0,P{ibf,2}));
end
ng  = sum(ngi);

% Get axis and boundary flux for each sub-bf
if nargin >= 5
  Finc = all(F1>=F0);
  [F0bf,F1bf] = deal(zeros(nbf,1));
  if Finc
    for ibf = 1:nbf
      iD = logical(lD{ibf}) & F1 > F0;
      if any(iD)
        F0bf(ibf) = min(F0(iD));
        F1bf(ibf) = max(F1(iD));
      end
    end
  else
    for ibf = 1:nbf
      iD = logical(lD{ibf}) & F1 < F0;
      if any(iD)
        F0bf(ibf) = max(F0(iD));
        F1bf(ibf) = min(F1(iD));
      end
    end
  end
end

switch mode
  case 0 % [FP,FT,TDG       ] = BF(0,PAR                    )
    [fPgi,fTgi,TDgi] = deal(cell(1,nbf));
    for ibf = 1:nbf
      [fPgi{ibf},fTgi{ibf}] = bfct{ibf}(0,bfp{ibf});
      TDgi{ibf} = repmat(lD{ibf},1,ngi(ibf));
    end
    y1 = vertcat(fPgi{:});
    y2 = vertcat(fTgi{:});
    y3 = horzcat(TDgi{:});
  case 1 % [TYG,TPDG,ITPDG  ] = BF(1,PAR, F,F0,F1,OY,RY,IRY )
    ny = numel(q1);
    Tyg = zeros(ny,ng);
    [TpDg,ITpDg] = deal(zeros(nD,ng));
    ig = 0;
    for ibf = 1:nbf
      ig = ig(end)+(1:ngi(ibf));
      for iD = find(lD{ibf} & (F1 - F0)).'
        [Tygi,TpDg(iD,ig),ITpDg(iD,ig)] = bfct{ibf}(1,bfp{ibf},F,F0bf(ibf),F1bf(ibf),int8(q1==iD),q2,q3);
        Tyg(:,ig) = Tyg(:,ig) + Tygi;
      end
    end
    y1 =   Tyg;
    y2 =  TpDg;
    y3 = ITpDg;
  case 2 % [GQDG,IGQDG      ] = BF(2,PAR,FQ,F0,F1           )
    nQ = numel(F);
    [gQDg,IgQDg] = deal(zeros(nQ,nD,ng));
    ig = 0;
    for ibf = 1:nbf
      ig = ig(end)+(1:ngi(ibf));
      for iD = find(lD{ibf} & (F1 - F0)).'
        FN = (F0(iD) + F*(F1(iD) - F0(iD)) - F0bf(ibf))/(F1bf(ibf) - F0bf(ibf));
        [gQDg(:,iD,ig),IgQDg(:,iD,ig)] = bfct{ibf}(2,bfp{ibf},FN,F0bf(ibf),F1bf(ibf));
      end
    end
    y1 =  gQDg;
    y2 = IgQDg;
    if nD == 1, y1 = reshape(y1,nQ,ng);
                y2 = reshape(y2,nQ,ng);
    end
  case 3 % [APP,ATTP,AP,AHQT] = BF(3,PAR, A,F0,F1,FP,FT,IDS )
    [aPpg,aTTpg,aPg,ahqTg] = deal(zeros(ng,1));
    ig = 0;
    for ibf = 1:nbf
      ig = ig(end)+(1:ngi(ibf));
      [aPpg(ig),aTTpg(ig),aPg(ig),ahqTg(ig)] = bfct{ibf}(3,bfp{ibf},F(ig),F0bf(ibf),F1bf(ibf),q1(ig),q2(ig),q3);
    end
    y1 =  aPpg;
    y2 = aTTpg;
    y3 =   aPg;
    y4 = ahqTg;
         % [TDG,TGY         ] = BF(4,PAR, F,F0,F1,O ,KD, FD )
  case 5 % [G0G,IG0G        ] = BF(5,PAR, ~,F0,F1           )
    ag = ones(ng,1); [fPg,fTg] = bfgenD(0,P);
    [gQDg,IgQDg           ] = bfgenD(2,P,0,F0,F1); % size(gQDg) = [1,nD,ng]
    [aPpg,aTTpg, aPg,ahqTg] = bfgenD(3,P,ag,F0,F1,fPg,fTg,1);
    y1 = reshape( gQDg,nD,ng).*(aPpg.*fPg+aTTpg.*fTg/mu0).'*2*pi;
    y2 = reshape(IgQDg,nD,ng).*( aPg.*fPg+ahqTg.*fTg/mu0).'*2*pi;
  case 6 % [QQG,XQ          ] = BF(6,PAR, ~,F0,F1,R0,IR0,IDS)
    [Qqgi,Xqi] = deal(cell(1,nbf));
    for ibf = 1:nbf
      iD = find(F0bf(ibf) == F0 & lD{ibf});
      [Qqgi{ibf},Xqi{ibf}] = bfct{ibf}(6,bfp{ibf},F,F0bf(ibf),F1bf(ibf),q1(iD),q2(iD),q3);
    end
    y1 = blkdiag(Qqgi{:});
    y2 = vertcat( Xqi{:});
  case 7 % [QCG,XC          ] = BF(7,PAR, ~,F0,F1           )
    [Qcgi,Xci] = deal(cell(1,nbf));
    for ibf = 1:nbf
      [Qcgi{ibf},Xci{ibf}] = bfct{ibf}(7,bfp{ibf},F,F0bf(ibf),F1bf(ibf));
    end
    y1 = blkdiag(Qcgi{:});
    y2 = vertcat( Xci{:});
         % [BTY             ] = BF(8,PAR, F,F0,F1,O ,  A,RBT,IDS,IRY)
  otherwise
    error('MEQ:BFGEND:CaseNotImplemented','Mode %d not implemented for %s',mode,mfilename);
end
end