classdef lih_tests < meq_test
  % Tests for LIH
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    tok = 'ana';
    t = 0;
    verbosity = 0;
    itert = 1;
    L,LX; 
    Lfbt,LYfbt;
  end

  properties(ClassSetupParameter)
    config = {'limited','diverted','doublet'};
    iterh = {0,10};
  end
  
  methods(TestClassSetup)
    function setup_LLX(testCase,config,iterh)
      switch config
        case 'limited'
          shot = 1;
        case 'diverted'
          shot = 2;
        case 'droplet'
          shot = 81;
        case 'doublet'
          shot = 82;
        otherwise
          error('undefined config %s\n',config)
      end
      tok = testCase.tok; t=testCase.t; %#ok<*PROPLC>
      
      [Lfbt,~,LYfbt] = fbt(tok,shot,t);
      [L] = lih(tok,shot,t,'iterh',iterh,...
        'debug',testCase.verbosity,'itert',testCase.itert,...
        'nelem',12);
      LX = meqxconvert(Lfbt,LYfbt,L);
      
      testCase.LX = LX;
      testCase.L = L;
      testCase.LYfbt = LYfbt;
      testCase.Lfbt = Lfbt;
    end
  end
  
  methods(Test,TestTags={'Unit'})
    function test_lih(testCase)
      L   = testCase.L; LX=testCase.LX; %#ok<*PROP>
      LYh = liht(L,LX);
      Lf  = testCase.Lfbt; LYf = testCase.LYfbt;
   
      if testCase.verbosity
        clf
        subplot(121);         
        imagesc(Lf.ry,Lf.zy,LYf.Iy); axis xy; title('FBT'); hold on;
        plot(LYf.rA,LYf.zA,'+k');
        plot(LYh.rA,LYh.zA,'or');

        subplot(122); 
        imagesc(L.ry,L.zy,LYh.Iy); axis xy; title('LIH'); hold on;
        plot(LYh.rA,LYh.zA,'or');
        drawnow
      end
      
      testCase.verifyTrue(~isempty(LYh),'liht returned empty');
      testCase.verifyEqual(LYh.Ip,LX.Ip,'AbsTol',2*L.P.Iperr,'Ip does not match')
      
      testCase.verifyEqual(LYh.zA,LYf.zA,'AbsTol',0.1)
      testCase.verifyEqual(LYh.rA,LYf.rA,'AbsTol',0.1 )
    end
  end
end
