classdef (SharedTestFixtures={mexm_fixture}) locR_tests < meq_test
  % tests of locRmex
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  properties
    verbosity = 0;
  end
  
  properties(TestParameter)
    method = {'locRmex','locRmexm'};
    sq     = struct('positive',-1,'negative',1);
  end
  
  methods(Test,TestTags={'Unit'})
    
    function test_basic(testCase,method,sq)
      x = [0;0.5;1]; y = [0;0.5;1]; yy=[0;0.5;1];
      xx = feval(method,x,sq.*y,sq.*yy);
      testCase.verifyEqual(yy,xx,'basic straight line interpolation failed')
    end
    
    function test_specials(testCase,method,sq)
      % check repeated x points
      x = [0;0;1];
      xx = feval(method,x,sq.*[0 1 2],sq.*[0 1 2],2);
      testCase.verifyEqual([x,nan(size(x))],xx,'failed for repeated x points')

      % check assigned default
      xx = feval(method,[0;1],sq.*[0 1],sq.*10,1,NaN);
      testCase.verifyTrue(all(isnan(xx)));      
    end
    
    function test_v_interp1(testCase,method,sq)
      %% A simple q profile
      aQ = 1 + linspace(0,1,21)'.^1.5; % unequally spaced grid
      qQ = 1 + 20*abs((aQ-1.4))'.^3; % some nonmonotonic profile
      
      
      % requested q values - some do not exist
      qR = [1,1.5,4,10];
      
      % amount of points to find
      
      % run locRmex()
      aR = feval(method,aQ,sq./qQ,sq./qR);
      
      %% check vs interp1 for monotonic part
      imono = find(diff(qQ)>0)+1;
      aR1 = interp1(1./qQ(imono),aQ(imono),1./qR,'linear',NaN)';
      
      testCase.verifyEqual(aR,aR1,'AbsTol',mean(diff(aQ))*1e-2)
      
      %% check that correct number of cases are found
      naR=10; % allow many
      % run locRmex()
      aR = feval(method,aQ,sq./qQ,sq./qR,naR);
      nR = sum(~isnan(aR),2);
      testCase.verifyEqual(nR,[0;2;1;0]); % hard coded number of points expected
      
      % optional plots
      if testCase.verbosity
        clf;
        plot(aQ,qQ);
        hold on
        for ii=1:naR
          iok = (aR(ii,:)~=-1);
          plot(aR(ii,iok),qR(iok),'ok')
        end
      end
      
      %% test empty
      aR = feval(method,aQ,sq./qQ,[],naR);
      testCase.assertEmpty(aR)
      
    end
  end
end