classdef (SharedTestFixtures={mexm_fixture}) fl4p_test < meq_test
  % Test for evaluation of fl4p 4-point flux interpolater
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbosity = 0;
  end
  
  methods(Test,TestTags = {'Unit'})
    function test_fl4p_c_and_m(testCase)
      
      % Generate test data
      L  = fgs('ana',1);
      
      R0 = 1.;
      Z0 = 0.;
      FA = 1; FB = 0;
      rBt = 1;
      [L,LY] = testCase.getCircularEquilibrium(L,R0,Z0,FA,FB,rBt);
      % End of data generation
      
      [zA,rA,~,~,~,~,~,~,~,~,~,~,stat] = ...
        asxymex(LY.Fx,L.G.zx,L.G.rx,L.P.dasm,L.dzx,L.drx,L.idzx,L.idrx,L.Oasx);
      testCase.verifyTrue(logical(stat));

      FN = L.FN*sign(LY.Ip);
      Fl = fl4pmex(LY.Fx,int32(L.kxl-1),L.clx,FN,L.G.rl,L.G.zl,rA,zA); % interpolate and keep only compatible sign of extrema

      Flm = fl4pmexm(LY.Fx,int32(L.kxl-1),L.clx,FN,L.G.rl,L.G.zl,rA,zA);
      %%
      if testCase.verbosity
        Flc = (Fl~=FN);
        contourf(L.G.rx,L.G.zx,LY.Fx,41); hold on;
        plot(L.G.rl,L.G.zl,'c'); axis equal
        plot(L.G.rl(Flc),L.G.zl(Flc),'ko'); axis equal
        title(sprintf('%s',mfilename));
      end
      %%
      testCase.verifyTrue(any(Fl~=FN),'No candidate limiter point found');
      testCase.verifyEqual(Fl,Flm,'RelTol',sqrt(eps));
    end
  end
  
end
