classdef gaps_test < meq_test
  % Tests for gaps calculations
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    tok = 'ana';
    t = 0;
    verbose = 0;
    iterq = 100; % Some targets are quite far from the origin 
    tol = 1e-3;
    L,LX;
  end

  properties(ClassSetupParameter)
    config = {'limited','diverted','secondaryX','limited-secondaryX'};
    sIp = {+1,-1}; % Ip sign
  end
  
  properties(TestParameter)
    nFW = {0,1,2,3};
    icsint = {false,true};
  end
  
  methods(TestClassSetup)
    function setup_LLX(testCase,config,sIp)
      switch config
        case 'limited'
          shot = 1;
        case 'diverted'
          shot = 2;
        case 'secondaryX'
          shot = 3;
        case 'limited-secondaryX'
          shot = 6;
        otherwise
          error('undefined config %s\n',config)
      end
      tok = testCase.tok; t=testCase.t;
      
      %FBT
      [Lfbt,LXfbt] = fbt(tok,shot,t,...
        'debug',testCase.verbose);
      LXfbt.IpD = LXfbt.IpD*sIp; % possible Ip sign flip
      LXfbt.qA  = LXfbt.qA *sIp; % flip qA for consistency
      LYfbt = fbtt(Lfbt,LXfbt);

      testCase.assertTrue(~isempty(LYfbt) && LYfbt.isconverged,'FBT did not converge')
      
      L = liu(tok,shot,t,...
        'iterq',testCase.iterq,...
        'npq',11,...
        'nW',16);
      
      testCase.LX = meqxconvert(Lfbt,LYfbt,L);
      testCase.L = L;
      
      if testCase.verbose
        clf;
        meqplott(Lfbt,LYfbt);
        drawnow;
      end
    end
  end
  
  methods(Test,TestTags={'Unit'})
    function test_gaps(testCase,nFW,icsint)
      
      %% LIUQE
      L = testCase.L; %#ok<*PROPLC>
      LX = testCase.LX;
      % Test with given amount of fluxes to track
      L.P.nFW = nFW;
      L.P.icsint = icsint;
      L = liuc(L.P,L.G);
      LY = liut(L,LX,'debug',testCase.verbose);
      
       %% Optional plotting
      if testCase.verbose
        np = ceil(numel(testCase.nFW)/2);
        subplot(np,np,nFW+1)
        meqplott(L,LY);
        [~,~,~,~,~,~,~,~,FX,~,~,~,~,~,...
          ~,~,FB,~,~,~] = meqpdom(LY.Fx,LY.Ip,L.P.isaddl,L);
        hold on;
        if numel(FX)==1, F=FX*[1 1]; else, F=FX; end
        contour(L.rx,L.zx,LY.Fx,F,'w:','linewidth',1.5)
        contour(L.rx,L.zx,LY.Fx,[FB FB],'k');
        title(sprintf('%d surfaces',nFW));
        drawnow
      end
      
      % Verifications
      testCase.assertNotEmpty(LY,'Returned empty LY');

      testCase.verifySize(LY.aW,[L.G.nW,nFW],sprintf('aW has unexpected size [%d,%d]',size(LY.aW)))
      
      if nFW>0 
        
        r = L.G.rW-LY.aW.*cos(L.G.oW);
        z = L.G.zW+LY.aW.*sin(L.G.oW);
        
        % interpolate fluxes on these points
        [rrx,zzx] = meshgrid(L.rx,L.zx);
        if icsint
          % Compute spline interpolation
          npar = L.csnthreads;
          rk = L.taur;
          zk = L.tauz;
          c = L.Mr*(L.Mz*LY.Fx).';
          Fi = reshape(bsp2deval(rk,zk,c,r(:),z(:),0,npar),size(r));
        else
          % Compute linear interpolation
          Fi = interp2(rrx,zzx,LY.Fx,r,z);
        end
        testCase.verifyEqual(Fi(:,1),repmat(LY.FW(1),L.G.nW,1),'AbsTol',1e-8*abs(LY.FA-LY.FB)); % Actual separatrix
        testCase.verifyEqual(Fi     ,repmat(LY.FW   ,L.G.nW,1),'AbsTol',1e-8*abs(LY.FA-LY.FB),...
          sprintf('Secondary separatrices were not reached for all gap specs\n\tsuccess fraction: [%s]/%d',...
          num2str(sum(abs(Fi - LY.FW )<1e-8*abs(LY.FA-LY.FB),1),'%d '),L.G.nW));
      end
    end
    
    function test_strikes(testCase)
      
      L = testCase.L; LX = testCase.LX; %#ok<*PROP>
      testCase.assumeTrue(L.P.shot == 2,'ignore test except for ''diverted'' config (shot=2)')
      
      %%
      PP=[fieldnames(L.P),struct2cell(L.P)]';
      L=liu(testCase.tok,L.P.shot,[],PP{:},...
        'aW',[  0.3  ; 1    ],...
        'rW',[  0.73 ; 1.24 ],...
        'zW',[ -0.4  ;-0.42 ],...
        'oW',[pi;0]);
      
      LY = liut(L,LX,'debug',testCase.verbose);
      if testCase.verbose
        clf;
        meqplott(L,LY);
        drawnow;
      end
      
      r = L.G.rW-LY.aW.*cos(L.G.oW);
      z = L.G.zW+LY.aW.*sin(L.G.oW);
      
      % interpolate fluxes on these points
      [rrx,zzx] = meshgrid(L.rx,L.zx);
      [Fi] = interp2(rrx,zzx,LY.Fx,r,z);
      testCase.verifyEqual(Fi,repmat(LY.FW,L.G.nW,1),'AbsTol',1e-8*abs(LY.FA-LY.FB));
      
    end
  end
end
