classdef vessel_tests < meq_test
  % Test for vessel representations and calculations
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties(TestParameter)
    ne = {256,30}
    selu = {'v','s','e','n'};
    tok  = {'tcv','ana'};
    code = {'liu','fbt','fgs','fge'}
  end
  
  properties
    L,LX
    solvetol = 1e-10;
    shot = 61400;
    insrc = 'liu';
    t = 0.5;
    nz = 16;
    nr = 16;
    verbosity = 0;
    excludedFields = {'Iu','res','niter','mkryl','nfeval','rese','resy','resg','resFx'}; % fields to exclude in checks
  end
  
  methods(TestClassSetup)
    function define_base_case(testCase)
      [testCase.L,testCase.LX] = ...
        fgs('tcv',testCase.shot,testCase.t,...
        'insrc',testCase.insrc,...
        'tolF',testCase.solvetol,...
        'debug',testCase.verbosity);
    end
  end
    
  methods (Test,TestTags={'fgs'})    
    function test_eigenmodes(testCase,ne)
      %% Project on smaller eigenmode space and solve
      [Le,LXe] = ...
        fgs('tcv',testCase.shot,testCase.t,...
        'insrc',testCase.insrc,...
        'tolF',testCase.solvetol,...
        'selu','e','nu',ne,'debug',testCase.verbosity);
      
      Iu2 = Le.G.Tvu*(Le.G.Tvu\testCase.LX.Iu); % project to eigenmode currents
      
      if testCase.verbosity>1
        plot([testCase.LX.Iu,Iu2]); drawnow
      end
      
      testCase.LX.Iu = Iu2;
      LY  = fgst(testCase.L,testCase.LX); % rerun full case with projected current
      LYe = fgst(Le,LXe); % solve directly with less eigenmodes

      testCase.verifyEqual(size(LYe.Iu,1),ne,'size of Iu does not match expectation')
      % compare excluding Iu (which changed) and res (which is very small)
      for ifield = fieldnames(LY)'
        field=ifield{:};
        switch field
          case testCase.excludedFields,  continue;
          otherwise
            if isempty(LY.(field))
              testCase.assertEmpty(LYe.(field));
            else
              testCase.assertEqual(LY.(field),LYe.(field),'AbsTol',sqrt(testCase.solvetol)*max(abs(LY.(field)(:))),...
                sprintf('field %s is not equal',field));
            end
        end
      end
    end
    
    function test_segments(testCase)
      % solve FGS using segment currents directly or projecting segment
      % currents onto vessel currents. Should give the same result.
      [Ls,LXs] = ...
        fgs('tcv',testCase.shot,testCase.t,...
        'insrc',testCase.insrc,...
        'tolF',testCase.solvetol,...
        'selu','s','debug',testCase.verbosity);
      LYs = fgst(Ls,LXs);      
      
      % project to segment currents
      Iv = Ls.G.Tvu*(Ls.G.Tvu\testCase.LX.Iu); % project to segment currents
      
      if testCase.verbosity>1
        plot([testCase.LX.Iu,Iv]); drawnow
      end
      
      LXm = testCase.LX; LXm.Iu = Iv; % modify Iv
      LY  = fgst(testCase.L,LXm); % run full case
      
      for ifield = fieldnames(LY)'
        field=ifield{:};
        switch field
          case testCase.excludedFields,  continue;
          otherwise
            if isempty(LY.(field))
              testCase.assertEmpty(LYs.(field));
            else
              testCase.assertEqual(LY.(field),LYs.(field),'AbsTol',sqrt(testCase.solvetol)*max(abs(LY.(field)(:))),...
                sprintf('failed check for %s',field));
            end
        end
      end
    end
  end
    
  methods(Test,TestTags={'Integration'})
    function test_inits(testCase,code,tok,selu)
      % test initialization for all codes for all vessel types for TCV and
      % anamak
      
      testCase.assumeFalse(isequal(tok,'tcv') && isempty(which('mdsopen')),...
        'skipping TCV test since mdsopen() not found');
      
      switch tok
        case 'tcv', shot = 61400; %#ok<*PROPLC>
        case 'ana', shot = 0;
      end
      fcode = str2func(code);
      try
        L = fcode(tok,shot,[],'selu',selu);
        testCase.verifyEqual(L.P.selu,selu);
      catch ME
        % expected error when running FGE with selu='s'
        if isequal(ME.identifier,'FGE:INVALIDVESSELMODE')
          return
        else
          rethrow(ME)
        end
      end
    end
  end
end
