classdef fl4pinterp_test < meq_test
  % Tests for fl4pinterp
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbosity = 0;
    L,LY;
  end
  
  properties(ClassSetupParameter)
    sIp = {-1,1};
    shot = {1,82}
  end
  
  properties(TestParameter)
    spike = {true,false};
    dr = num2cell([-0.04,0.04]);
    dz = num2cell([-0.04,0.04]);
  end
  
  methods(TestClassSetup)
    function setup_eq(testCase,sIp,shot)
      [L,LX] = fbt('ana',shot,[],'nl',16,...
        'ilim',2); % ilim=2 enables fl4pinterp
      LX.IpD = sIp*LX.IpD;
      LX.qA  = sIp*LX.qA;
      LY = fbtt(L,LX);
      if testCase.verbosity
        clf; meqplott(L,LY)
      end
      testCase.L=L; testCase.LY=LY; %#ok<*PROP>
    end
  end
  
  methods(Test, TestTags={'Unit'})
    function test_limiter_point(testCase,spike,dr,dz)
      %% Test that limiter points ends up where expected
      % Keep flux distribution fixed, just recompute limiter point
      
      L = testCase.L; LY = testCase.LY; %#ok<*PROPLC>
      
      rA = LY.rA; zA = LY.zA;
      Fx = LY.Fx; % flip with Ip sign
      sIp = sign(LY.Ip); Ip = LY.Ip;
      
      % Move limiter around
      rl = 0.9*(L.G.rl - L.P.r0)+ dr + L.P.r0;
      zl = 0.9*L.G.zl + dz;
      
      if spike
        % stick piece of limiter into the plasma, this should be the
        % limiter point;
        rl = [rl;rl(1)-0.8*(rl(1)-rA(1))];
        zl = [zl;zl(1)-0.8*(zl(1)-zA(1))];
      end
      G = L.G; P = L.P;
      G.rl = rl; G.zl=zl;        % hack limiter
      G = meqg(G,P);             % recompute geometry
      L = meqc(P,G,'rrx','zzx'); % consolidate
      
      % Find boundary point in this situation
      [~,~,~,~,~,~,~,~,~,~,~,~,~,~,...
       rB,zB,FB] = meqpdom(Fx,Ip,L.P.isaddl,L);
      
      % Check spike specifically
      if spike
        testCase.verifyTrue(any(rB==rl(end)) && any(zB==zl(end)),...
          'Boundary should be on spike since it sticks into the plasma')
      end
      
      % check against high-precision limiter
      rrl = [L.G.rl;L.G.rl(1)]; zzl = [L.G.zl;L.G.zl(1)];
      [rrll,zzll] = deal([]); % interpolated limiter points
      nint = 30; % interpolation factor
      for ii=1:numel(rrl)-1
        rrll = [rrll;rrl(ii) + linspace(0,1,nint)'*(rrl(ii+1)-rrl(ii))]; %#ok<AGROW>
        zzll = [zzll;zzl(ii) + linspace(0,1,nint)'*(zzl(ii+1)-zzl(ii))]; %#ok<AGROW>
      end
      FFll = interp2(L.rx,L.zx,Fx,rrll,zzll);
      
      [~,imax] = max(sIp*FFll); [~,iFBmax] = max(sIp*FB);
      testCase.verifyEqual(rrll(imax),rB(iFBmax),'absTol',2e-2)
      testCase.verifyEqual(zzll(imax),zB(iFBmax),'absTol',2e-2)
      testCase.verifyEqual(FFll(imax),FB(iFBmax),'absTol',2e-3)
      
      if testCase.verbosity
        %%
        clf; dg = [0 0.6 0]; % dark green
        contour(L.rx,L.zx,LY.Fx-LY.FB(1),21); hold on;
        contour(L.rx,L.zx,LY.Fx,unique(LY.FB.*[1,1]),'linewidth',2,'color','k');
        plot(rrl,zzl,'k-',L.G.rl,L.G.zl,'o-k'); % original limiter
        plot(rB,zB,'o','markersize',16,'markerfacecolor',dg); % winner points
        plot(rrll(imax),zzll(imax),'.r','markersize',15); % point on high-density interpolated limiter
        plot(rl,zl,'go','markersize',10); % interpolated limiter candidates
        plot(L.rrx(L.kxl),L.zzx(L.kxl),'sb')
        plot(L.rrx(L.kxl+[0 L.nzx 1 L.nzx+1]),L.zzx(L.kxl+[0 L.nzx 1 L.nzx+1]),'sb')
        contour(L.rx,L.zx,Fx,unique(FB.*[1,1]),'r','linewidth',2)
        
        axis tight equal;
        drawnow;
      end
    end
  end
end
