classdef (SharedTestFixtures={mexm_fixture}) ...
 asxy_other_test < meq_test
  % Use asxy_test for generic asxy tests across multiple configurations
  % Use asxy_other_test (this file) for specific issues necessiting a precise configuration
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbosity = 0;
  end
  
  properties (TestParameter)
    method = {'asxymex','asxymexm'}
    R0     = struct('symetric',1,'asymetric',1.001)
    lasy   = {true,false}
  end
  
  methods (Test, TestTags={'Unit'})
    function duplicate_value_near_axis(testCase,method,R0)
      L  = fgs('ana',1);
      
      % Fake parabolic Fx
      Z0 = 0;
      FA = 1; FB = 0;
      rBt = 1;
      [L,LY] = testCase.getCircularEquilibrium(L,R0,Z0,FA,FB,rBt);
      
      % Localise extrema and saddle points of the Fx map
      zA = feval(method,LY.Fx,L.G.zx,L.G.rx,L.P.dasm,L.dzx,L.drx,L.idzx,L.idrx,L.Oasx);
      testCase.verifyEqual(numel(zA),1,'Simple flux map has a single axis');
    end
    
    function flat_gradient_null(testCase,method)
      
      % Without clamping this map produces an X-point at [-23,15]
      %   (this was randomly generated until the above result was obtained)
      Fx =    [0.235077636686170   0.729397698256153   0.428442475951073;...
               0.461602080891622   0.724916747219119   0.864686144763534;...
               0.198591503053926   0.330391746361703   0.657065424622890;...
               0.025476235778820   0.843494126063655   0.581035657418090];
             
      rx = 0:2; drx = 1;
      zx = 0:3; dzx = 1;
      dasm = 0;
      
      
      [zA,rA,~,~,~,~,zX,rX] = feval(method,Fx,zx,rx,dasm,dzx,drx,1./dzx,1./drx,true(size(Fx)));
      
      if ~isempty(zA)
        testCase.verifyLessThanOrEqual(rA,max(rx));
        testCase.verifyLessThanOrEqual(zA,max(zx));
        testCase.verifyGreaterThanOrEqual(rA,min(rx));
        testCase.verifyGreaterThanOrEqual(zA,min(zx));
      end
      if ~isempty(zX)
        testCase.verifyLessThanOrEqual(rX,max(rx));
        testCase.verifyLessThanOrEqual(zX,max(zx));
        testCase.verifyGreaterThanOrEqual(rX,min(rx));
        testCase.verifyGreaterThanOrEqual(zX,min(zx));
      end   
    end
    
    function lasy_test(testCase,lasy)
      
      % Run FBT case with X-point outside limiter
      [L,LX,LY] = fbt('ana',21,[],'lasy',lasy);
      
      % Check whether X-point is identified
      dtol = max(L.drx,L.dzx)/10; % tolerance: 1/10th grid cell
      iX = (LX.gpbr==0 & LX.gpbz==0);
      check = any(sqrt((LY.rX - LX.gpr(iX)).^2 + (LY.zX - LX.gpz(iX)).^2) < dtol);
      
      testCase.verifyEqual(check,lasy,sprintf('Error in out of limiter X-point identification, lasy=%d, check(is X-point in LY)=%d',lasy,check))
      
    end

  end
  
end
