classdef meqpdom_test < meq_test
  % Tests for meq plasma domain finding module
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbose = 0;
    writegif = 1;
    tmpfilename;
    L;
    Ip;
  end

  properties(ClassSetupParameter)
    ilim   = struct('standard',    1,'refined',   3);
    icsint = struct('standard',false,'refined',true);
  end
  
  properties(MethodSetupParameter)
    sIp = {-1,1};  % Ip signs to test
  end
  
  properties(TestParameter)
    type = {
      'Lim','Lim-X','SND','DND',...
      'Double-Snowflake-Minus','Double-Snowflake-Plus',...
      'Doublet','Droplet',...
      'Doublet-div','Doublet-div-nomantle',...
      'Triplet','Triplet-madness'};
  end
  
  methods(TestClassSetup, ParameterCombination = 'sequential')
    function getL(testCase,ilim,icsint)
      testCase.L = fgs('ana',0,0,'cappav',2,'nc',16,'ac',0.6,...
        'debug',testCase.verbose,'nz',32,'nr',32,'infct',@qintmex,...
        'ilim',ilim,'icsint',icsint);
    end
    
    function setupGif(testCase)
      if testCase.writegif
        tmpfilename_ = sprintf('%s_animated_tmp.gif',mfilename);
        filename     = sprintf('%s_animated.gif',mfilename);
        testCase.addTeardown(@() movefile(tmpfilename_,filename));
        testCase.tmpfilename = tmpfilename_;
      end
    end
  end
  
  methods(TestMethodSetup)
    function setIp(testCase,sIp)
      testCase.Ip = sIp;
    end
  end
  
  methods(Test, TestTags = {'Unit'})
    function test_configurations(testCase,type)
      L = testCase.L; %#ok<*PROPLC>
      sIp = sign(testCase.Ip);
      isaddl = 1;
      %% configurations
      Ia = zeros(L.G.na,1);
      
      S = testCase.generate_flux_map(type,sIp);
      
      Fx = zeros(size(L.rrx));
      Fxp = S.Fxh(L.rrx,L.zzx);
      Fx(:) = L.G.Mxa*Ia + Fxp(:);
      
      [rA,zA,~,~,~,~,rX,zX,~,~,~,~,~,~,...
        rB,zB,~,~,lX,Opy] = meqpdom(Fx,sIp,isaddl,L);
      
      % Try meqpdom on this
      if testCase.verbose || testCase.writegif
        subplot(121); cla;
        contourf(L.G.rx,L.G.zx,Fx,21); hold on;
        plot(L.G.rl,L.G.zl,'k','linewidth',2);
        plot(rB,zB,'r.','markersize',10);
        axis equal tight
        title(sprintf('%s,sIp=%+1.0f',type,sIp));
        subplot(122); cla
        drawnow;
      end
      
      testCase.verifyEqual(numel(rA),S.nA,'Mismatch in number of numel(rA)');
      testCase.verifyEqual(numel(zA),S.nA,'Mismatch in number of numel(zA)');
      testCase.verifyEqual(numel(rX),S.nX,'Mismatch in number of numel(rX)');
      testCase.verifyEqual(numel(zX),S.nX,'Mismatch in number of numel(zX)');
      testCase.verifyEqual(numel(unique(Opy))-1,S.ndom,'Mismatch in number of domains');
      testCase.verifyEqual(lX,S.lX,'Mismatch in diverted state');
            
      if testCase.writegif
        if testCase.verbose, fprintf('Writing gif for %s\n',type);end
        % Capture the plot as an image
        im = frame2im(getframe(gcf));
        [imind,cm] = rgb2ind(im,256);
        % Write to the GIF File
        if ~exist(testCase.tmpfilename,'file')
          imwrite(imind,cm,testCase.tmpfilename,'gif', 'Loopcount',inf);
        else
          imwrite(imind,cm,testCase.tmpfilename,'gif','WriteMode','append');
        end
      end
    end
  end
end
