%BFHELP  Calling scheme for base functions for p',TT'
% [FP,FT           ] = BF(0,PAR                            )
% [TYG,TPDG,ITPDG  ] = BF(1,PAR, F,F0,F1,O ,RY,IRY         )
% [GQDG,IGQDG      ] = BF(2,PAR,FN,F0,F1                   )
% [APP,ATTP,AP,AHQT] = BF(3,PAR, A,F0,F1,FP,FT,IDS         )
% [TDG,TGY         ] = BF(4,PAR, F,F0,F1,O ,KD, FD         )
% [G0G,IG0G        ] = BF(5,PAR, ~,F0,F1                   )
% [QQG,XQ          ] = BF(6,PAR, ~,F0,F1,R0,IR0,IDS        )
% [QCG,XC          ] = BF(7,PAR, ~,F0,F1                   )
% [BTY             ] = BF(8,PAR, F,F0,F1,O ,  A,RBT,IDS,IRY)
% 
%
% BF(0,...) returns the factors FP[#g] and FT[#g] weighting p' and TT'
%
% BF(1,...) for a given flux map F(z,r) returns TYG[(#z-2)*(#r-2),#g], the
% #g base functions {g'(F(z,r))} evaluated on the inner grid of F,
% multiplied by FP*r+FT/r. TPDG=SUM_Y(TYG), ITPDG=SUM_Y(ITYG) with ITYG
% similarly built from the corresponding primitives {g}.
%
% BF(2,...) returns GN(nN,ng)={g'(FN)} and IGN={g(FN)} evaluated for #N
% normalised flux values FN.
%
% BF(3,...) for coefficients A such that the plasma current distribution
% IY=TYG*A, returns coefficients such that p'=GN*APP, TT'=GN*ATTP,
% p=IGN*AP, T^2/2=IGN*AHQT.
%
% BF(4,...) is used for diagnostic with vertical chord integrals and
% returns TDG={SUM_Z(g'(F(z,r(KD)))*FD(1,:)+g'(F(z,r(KD+1)))*FD(2,:))}.
%
% BF(5,...) returns the value on magnetic axis of G[#g]={g'} and IG={g}.
%
% BF(6,...) returns QQG[#q,#g] and XQ[#q,1] for the #q regularisation
% least square equations XQ=QQG*AG.
%
% BF(7,...) returns QCG[#c,#g] and XC[#c,1] for the #c constraint
% equations XC=QCG*AG.
%
% BF(8,...) returns the toroidal field for the coefficients A and vacuum
% field RBT
%
% F0,F1 are the flux on axis and at the boundary of each domain; the int8 
% array O(z,r)) flags the points of the inner grid in the plasma domain,
% IRY=1/r(2:#r-1), IDS=1/dr/dz. BF is called using the function handle
% P.bfct and PAR=P.bfp, P as returned by LIUP.
%
% For details see: [MEQ-redbook] 
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function bfhelp

 %% Check
 n = 41;
 FN = linspace(0,1,n)';
 Fx = zeros(n+2,3); Fx(2:end-1,2) = FN;
 bf = {@bffbt                    @bfabmex @bf3pmex};
 p  = {[1.1 .25 1.5 2.2 .75 2.5] [3 0]    false};
 FA = [-1                        0      	 0];
 FB = [0                         1        1];
 for k = 1:length(bf)
  [fP,fT   ] = bf{k}(0,p{k});
  [~,~,~   ] = bf{k}(1,p{k},Fx*(FB(k)-FA(k))+FA(k),FA(k),FB(k),true(n,1),ones(n,1),ones(n,1));
  [GN ,IGN ] = bf{k}(2,p{k},FN                    ,FA(k),FB(k)                              );
  [~,~,~,~ ] = bf{k}(3,p{k},ones(length(fP),1)    ,FA(k),FB(k),fP,fT,1.                     );
  [~,~     ] = bf{k}(5,p{k},[]                    ,FA(k),FB(k)                              );
  [~,~     ] = bf{k}(6,p{k},[]                    ,FA(k),FB(k),1,1,1                        );
  IGN2 = bfprmex(GN);
  subplot(2,2,k)
  plot(FN,GN,'-k',FN,IGN,'-xk',FN,IGN2,'-+r')
  title(func2str(bf{k}))
 end
 
 %% Maths for analytical base functions

 %% BFAB, BF3P
 disp('BFAB,BF3P')
 syms F FA FB
 syms FN positive
  g  = [ F - FB;
        (F - FB)*(F - FA);
        (F - FB)*(F - FA)*((F-FB)+(F-FA))];
 Ig  = int(g,F);
 Ig  = Ig - subs(Ig,F,FB); % 0 on LCFS
  gN = subs(g,F,(FB-FA)*FN+FA);
 IgN = subs(Ig,F,(FB-FA)*FN+FA);
 dsp1(g,Ig,gN,IgN)

 %% BFFBT
 disp('BFFBT')
 syms F FA FB a1 b1 a2 b2
 syms FN f c1 c2 positive
  gN = [a1*(b1+(1-b1)*FN)*(1-FN)^c1;
        a2*(1-FN)^c2;
        a2*(b2*FN-1-b2)*FN*(1-FN)^c2];
  gN = [a1*(b1+(1-b1)*(1-f))*f^c1;
        a2*f^c2;
        a2*(f-1)*(1+b2*f)*f^c2];
  gN = subs(gN,f,1-FN);
 IgN = int(gN,FN);
 IgN = IgN - subs(IgN,FN,1); % 0 on LCFS
 dsp1(sym([]),sym([]),gN,IgN)
  gf = subs( gN,FN,1-f);
 Igf = subs(IgN,FN,1-f);
 dsp2( gf,f,1)
 dsp2(Igf,f,1)

end

function dsp1(g,Ig,gN,IgN)
 syms F FN FA
 dsp2( g ,sym('F' ),sym('FA'))
 dsp2(Ig ,sym('F' ),sym('FA'))
 dsp2( gN,sym('FN'),0        )
 dsp2(IgN,sym('FN'),0        )
end

function dsp2(y,x,a)
 fprintf('%s\n',inputname(1)); pretty(y)
 for k = 1:length(y)
  f{k} = factor(y(k),x);
  fprintf(1,'factors for %s%d\n',inputname(1),k); pretty(f{k})
  for l = 1:length(f{k})
   fprintf(1,'coeffs for %s%d, factor %d, increasing power\n',inputname(1),k,l); pretty(coeffs(f{k}(l),x))
  end
 end
 fprintf('value on axis %s\n',inputname(1));
 pretty(subs(y,x,a))
end
