classdef mds2meq_test < meq_test
% tests of MDS retrieval of MEQ structures
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

  properties(TestParameter)
    source = {'LIUQE.M', 'FBTE','RTLIUQE'};
    shots  = {'70474','70477','70563','73306'};
  end
  
  methods(Test,TestTags={'TCV'})
    function mds2meq_data_retrieval(testCase,source,shots)
      shot = str2double(shots);
      if strcmp(source,'RTLIUQE')
        testCase.assumeFalse(ismember(shot,[70563,73306]),...
          sprintf('skip shot %d for RTLIUQE since SCD did not run',shot))
      end
      
      skipfields = {'shot'};
      [L,LY] = mds2meq(shot,source);
      nt = numel(LY.t);
      for field=fieldnames(LY)'
        myfield = field{:};
        if contains(skipfields,myfield), continue; end
        tsize = meqsize(L,myfield);
        mysize = size(LY.(myfield));
        if tsize(2)==1
          expectedsize = [tsize(1) nt];
        else
          expectedsize = [tsize,nt];
          if nt==1
            mysize  = [mysize,nt]; %#ok<AGROW>
          end
        end
        testCase.verifyEqual(mysize,expectedsize,...
          sprintf('invalid size of retrieved data for %s',myfield))
      end
    end
    
    function mds2meq_units(testCase,source,shots)
      shot = str2double(shots);
      if strcmp(source,'RTLIUQE')
        testCase.assumeFalse(ismember(shot,[70563,73306]),...
          sprintf('skip shot %d for RTLIUQE since SCD did not run',shot))
      end
      
      [L,LY] = mds2meq(shot,source);
        
      testCase.assumeTrue(L.nD==1,sprintf('Skipping mds2meq Iy test for shot %d: multiple domains',shot));
      
      % Test slice with largest current
      [~,kt] = max(abs(LY.Ip));

      Fx = LY.Fx(:,:,kt);
      if ~isfield(LY,'Opy')
        if strcmp(source,'FBTE'), isaddl = L.P.isaddl(kt);
        else,                     isaddl = L.P.isaddl;    end
        
        [~,~,FA,~,~,~,~,~,~,~,~,~,~,~,...
          ~,~,FB,lB,~,Opy] = meqpdom(Fx,LY.Ip(kt),isaddl,L);
        
        testCase.assumeTrue(lB && numel(FA)==1,sprintf('Skipping mds2meq Iy test for shot %d: meqpdom failed',shot));
      else
        FA  = LY.FA (    kt);
        FB  = LY.FB (    kt);
        Opy = LY.Opy(:,:,kt);
      end
        
      FNy = (Fx(2:end-1,2:end-1)-FA)/(FB-FA);
      Ppy  = interp1(L.pQ.^2,LY.PpQ(:,kt) ,FNy,'spline'); Ppy (~Opy) = 0;
      TTpy = interp1(L.pQ.^2,LY.TTpQ(:,kt),FNy,'spline'); TTpy(~Opy) = 0;
      
      Iy = 2*pi*(L.ry.'.*Ppy + TTpy.*L.iry.'/mu0)*L.dsx;
      
      % Test with low tolerance, allowing for differences in models
      testCase.verifyEqual(Iy,LY.Iy(:,:,kt),'AbsTol',max(max(abs(LY.Iy(:,:,kt))))*1e-2,...
        'Iy from PpQ/TTpQ does not match MDS storage within required tolerance');
      
    end
  end
  
  
end
