classdef (SharedTestFixtures={meq_simulink_fixture}) lihsim_tests < meq_test
  % Tests for Simulink version of LIH
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbosity = 0;
  end
  
  properties (TestParameter)

    tok  = {'TCV'    ,'ana', 'TCV'    ,'ana'};
    t1   = {2.40     , 0   , 2.40     , 0   };
    t2   = {2.43     , 0.1 , 2.43     , 0.1 };
    dt   = {1e-3     , 1e-3, 1e-3     , 1e-3};
    shot = {68898    , 2   , 68898    , 2   };
    iterh= {0        , 0   , 10       , 10  };
 idoublet= {1        , 0   , 1        ,0    };

  end
  
  methods (Test, TestTags={'Simulink'},ParameterCombination='sequential')
    
    function slx_vs_m(testCase,tok,shot,t1,t2,dt,iterh,idoublet)
          
      [ok,msg] = meq_test.check_tok(tok);
      testCase.assumeTrue(ok,msg);
      
      time = t1:dt:t2; nt = numel(time);
      
      if idoublet, ndom=2; else, ndom=1; end
      %% Get L,LX
      [L,LX] = lih(tok,shot,time,'iterh',iterh,'slx',true,'iterq',0,...
        'itert',1,'iqR',1/2,'raS',0.5,'idoublet',idoublet,...
        'bfct',@bf3pmex,'bfp',false,'ndom',ndom);

      % run Matlab
      LYm = liht(L,LX,'slx',false);
      
      % run Simulink
      LYs = lihtsim(L,LX);      
      
      %% Check outputs
      fields=fieldnames(LYs); nfields = numel(fields);

      % Errors on rA,zA can reach high values when IpD nears 0.
      ignoredFields = {'nA','lB','Fx','Iy','rA','zA','shot','rX','zX','FX'};

      for ifield = 1:nfields
        myfield = fields{ifield};
        if contains(myfield,ignoredFields), continue; end % ignore this field
        if ~isfield(LYm,myfield), continue; end % skip non-existent fields
        
        LYsval = LYs.(myfield);
        LYmval = LYm.(myfield);
        
        dimssame = isequal(size(LYsval),size(LYmval));
        testCase.verifyEqual(size(LYsval),size(LYmval),sprintf('dimensions not equal: %s\n',myfield));
       
        if dimssame && ismatrix(LYsval)
          % check numerical error
          ds = reshape(LYsval,[],nt);
          dm = reshape(LYmval,[],nt);
          absErr = norm(double(ds-dm));
          relErr = absErr/norm(double(ds));
          if testCase.verbosity
            %% plot
            clf;
            subplot(2,1,1)
            plot(LYs.t,ds,'b',LYm.t,dm,'r--')
            title(sprintf('%s: abs;%4.3e rel:%4.3e: .slx(b), .m(r--) \n',myfield,absErr,relErr));
            
            subplot(2,1,2)
            plot(LYs.t,sum(ds-dm,1)/size(ds,1));
            title('error');
            pause(0.1);
          end
        
          %% compare
          ok = relErr < 1e-2 || absErr < 1e-3;
          testCase.verifyTrue(ok,sprintf('%s values did not match within tolerance',myfield))
        end
      end
    end
  end
end
