classdef bf_multiD_tests < meq_test
  % Tests for multidimensional basis functions
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    L,LY,bfct,bfp
  end
  
  methods (TestClassSetup)
    function setupLLY(testCase)
      [L,~,LY] = fbt('ana',82,[]); %#ok<*PROP>
      testCase.L  = L;
      testCase.LY = LY;
      
      testCase.bfct = @bfgenD;
      testCase.bfp  = {@bf3pmex,false,[1;0;0];@bf3pmex,false,[0;1;0];@bf3pmex,false,[1;1;1]};
    end
  end
      
  
  methods (Test, TestTags = {'Unit','bf'})
    %% Verify that each mode runs and yields expected size of outputs
    function check_mode_0(testCase)
      [~,~,bfct,bfp] = deal(testCase.L,testCase.LY,testCase.bfct,testCase.bfp);
      
      [fPg,fTg,TDg] = bfct(0,bfp);
      
      testCase.verifyEqual(size(fPg),[9,1]);
      testCase.verifyEqual(size(fTg),[9,1]);
      testCase.verifyEqual(size(TDg),[3,9]);
    end
    
    function check_mode_1(testCase)
      [L,LY,bfct,bfp] = deal(testCase.L,testCase.LY,testCase.bfct,testCase.bfp);
      
      [F0,F1] = meqpdomlim(LY.FA,LY.FB,L.nD);
      
      [Tyg,TpDg,ITpDg] = bfct(1,bfp,LY.Fx,F0,F1,LY.Opy,L.ry,L.iry);
      
      testCase.verifyEqual(size(  Tyg),[L.ny,9]);
      testCase.verifyEqual(size( TpDg),[L.nD,9]);
      testCase.verifyEqual(size(ITpDg),[L.nD,9]);
    end
    
    function check_mode_2(testCase)
      [L,LY,bfct,bfp] = deal(testCase.L,testCase.LY,testCase.bfct,testCase.bfp);
      
      [F0,F1] = meqpdomlim(LY.FA,LY.FB,L.nD);
      
      [gQDg,IgQDg] = bfct(2,bfp,L.pQ.^2,F0,F1);
      
      if L.nD > 1, size_ = [L.nQ,L.nD,9];
      else,        size_ = [L.nQ,9];
      end
      testCase.verifyEqual(size( gQDg),size_);
      testCase.verifyEqual(size(IgQDg),size_);
    end
    
    function check_mode_3(testCase)
      [L,LY,bfct,bfp] = deal(testCase.L,testCase.LY,testCase.bfct,testCase.bfp);
      
      [F0,F1] = meqpdomlim(LY.FA,LY.FB,L.nD);
      [fPg,fTg,~] = bfct(0,bfp);
      
      [aPpg,aTTpg,aPg,ahqTg] = bfct(3,bfp,ones(9,1),F0,F1,fPg,fTg,L.idsx);
      
      testCase.verifyEqual(size( aPpg),[9,1]);
      testCase.verifyEqual(size(aTTpg),[9,1]);
      testCase.verifyEqual(size(  aPg),[9,1]);
      testCase.verifyEqual(size(ahqTg),[9,1]);
    end
    
    function check_mode_5(testCase)
      [L,LY,bfct,bfp] = deal(testCase.L,testCase.LY,testCase.bfct,testCase.bfp);
      
      [F0,F1] = meqpdomlim(LY.FA,LY.FB,L.nD);
      
      [g0g,Ig0g] = bfct(5,bfp,[],F0,F1);
      
      testCase.verifyEqual(size( g0g),[L.nD,9]);
      testCase.verifyEqual(size(Ig0g),[L.nD,9]);
    end
    
    function check_mode_6(testCase)
      [L,LY,bfct,bfp] = deal(testCase.L,testCase.LY,testCase.bfct,testCase.bfp);
      
      [F0,F1] = meqpdomlim(LY.FA,LY.FB,L.nD);
      r0  = [   LY.rA;zeros(L.nD-LY.nA,1)];
      ir0 = [1./LY.rA;zeros(L.nD-LY.nA,1)];
      
      [Qqg,Xq] = bfct(6,bfp,[],F0,F1,r0,ir0,L.idsx);
      
      testCase.verifyEqual(size(Qqg,2),9);
      testCase.verifyEqual(size(Xq,1),size(Qqg,1));
    end
    
    function check_mode_7(testCase)
      [L,LY,bfct,bfp] = deal(testCase.L,testCase.LY,testCase.bfct,testCase.bfp);
      
      [F0,F1] = meqpdomlim(LY.FA,LY.FB,L.nD);
      
      [Qcg,Xc] = bfct(7,bfp,[],F0,F1);
      
      testCase.verifyEqual(size(Qcg,2),9);
      testCase.verifyEqual(size(Xc,1),size(Qcg,1));
    end
    
    %% check continuity across separatrix of basis functions
    function bf_continuity_test(testCase)
      [L,LY,bfct,bfp] = deal(testCase.L,testCase.LY,testCase.bfct,testCase.bfp);
      
      [F0,F1] = meqpdomlim(LY.FA,LY.FB,L.nD);
      
      [gQDg,IgQDg] = bfct(2,bfp,L.pQ.^2,F0,F1);
      
      testCase.verifyEqual( gQDg(end,1,:), gQDg(1,3,:),'bf functions are discontinuous across separatrix 1->3');
      testCase.verifyEqual(IgQDg(end,1,:),IgQDg(1,3,:),'bf integrals are discontinuous across separatrix 1->3');
      testCase.verifyEqual( gQDg(end,2,:), gQDg(1,3,:),'bf functions are discontinuous across separatrix 2->3');
      testCase.verifyEqual(IgQDg(end,2,:),IgQDg(1,3,:),'bf integrals are discontinuous across separatrix 2->3');
      
    end
  end
  
end