classdef (SharedTestFixtures={mexm_fixture}) rtci_test < meq_test
  % Tests for RTCI (real-time contour integrals)
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    L; LY;
    rel_tol = 1;
    verbosity = 0;
  end
  
  properties (MethodSetupParameter)
    nz = {32  64};
    icsint = {false, true};
  end
  
  properties (TestParameter)
    method = {'rtcimex','rtcimexm'};
  end
    
  methods(TestMethodSetup)
    function rtciTestSetup(testCase,nz,icsint)
      % Several surfaces to track
      fq = [0.5 0.7 0.95 1];
      
      L = liu('ana',1,1,'iterq',1,'nz',nz,'pq',sqrt(fq),'icsint',icsint);

      % calculate axis location and radius of circle a0
      r0 = 1.05; 
      z0 = 0.05;
      FA = 0;
      FB = 1;
      rBt = 1.4;
      
      [testCase.L,testCase.LY] = testCase.getCircularEquilibrium(L,r0,z0,FA,FB,rBt);
    end
  end

  methods(Test, TestTags = {'Unit'})
    function test_rtci_contour(testCase,method)
      
      %%
      L = testCase.L; %#ok<*PROPLC> % for brevity
      LY = testCase.LY;
      
      % Main iteration to compute level curve for Fx = 1 (this call was adapted from meqpostq.m)
      aq = repmat(min([LY.rA-L.rx(1) L.rx(end)-LY.rA LY.zA-L.zx(1) L.zx(end)-LY.zA])*...
        L.pinit*L.pq,L.noq,1); % initial a
      
      niter = 200; % iterations to converge contour calculation
      F = LY.FB-L.fq*(LY.FB-LY.FA);
  
      if L.P.icsint
        aq = rtcics(LY.Fx,F,LY.rA,LY.zA,L.crq,L.czq,L);
        daq = 0;
      else
        ero=(LY.rA-L.rx(1))*L.idrx;
        ezo=(LY.zA-L.zx(1))*L.idzx;
        for ki = 1:niter
          [aq,daq] = feval(method,aq,ero,ezo,LY.Fx,F,L.cdrq,L.cdzq,LY.Opy,LY.FA,int8(1),L.drx);
        end
      end
      % final call to also get r,z of final contour
      rq = LY.rA + L.crq.*aq;
      zq = LY.zA + L.czq.*aq;

      tol = max(L.drx,L.dzx)*testCase.rel_tol;
      testCase.verifyLessThan(max(abs(LY.aq-aq)),tol,'.mex result failure')
      testCase.verifyLessThan(daq,tol,'rtci did not converge')
            
      if testCase.verbosity
        %%

        txt = sprintf('%dx%d',L.nzx,L.nrx);
        subplot(121)
        plot(LY.rA,LY.zA,'ok')
        hold on
        contour(L.rx,L.zx,LY.Fx,linspace(0,3,31),'-r');
        contour(L.rx,L.zx,LY.Fx,LY.FB*[1 1],'-k','linewidth',2)
        plot(LY.rq,LY.zq,'xm',rq,zq,'ob')
        hold off, axis equal
        title(sprintf('RTCI test: %s, nz=%d',L.P.tokamak,L.nzx))
        
        subplot(122)
        semilogy(L.oq,abs(LY.aq-aq))
        xlabel('\theta'), xlim([0 2*pi]), 
        title('LCFS position error')
        legend(txt,'location','best')
        orient tall
        drawnow; shg;
      end
      
    end
    
    function test_rtci_gaps(testCase,method)
      
      %%
      L = testCase.L;
      LY = testCase.LY;
      
      % Main iteration to compute level curve for Fx = 1 (this call was adapted from liut.m)
      aW = zeros(L.G.nW,L.P.nFW); % initial a
      
      % Values at origin points
      FoW = bintmex(LY.Fx,int32(L.kxW-1),L.cWx);
      
      niter = 200; % iterations to converge contour calculation
      FW = LY.FB;
  
      if L.P.icsint
        aW = rtcics(LY.Fx,FW,L.G.rW,L.G.zW,-cos(L.G.oW),sin(L.G.oW),L);
        daW = 0;
      else
        for ki = 1:niter
          [aW,daW] = feval(method,aW,L.erW,L.ezW,LY.Fx,FW,L.cdrW,L.cdzW,LY.Opy,FoW,int8(0),L.drx);
        end
      end

      tol = max(L.drx,L.dzx)*testCase.rel_tol;
      testCase.verifyLessThan(max(abs(LY.aW-aW)),tol,'.mex result failure')
      testCase.verifyLessThan(daW,tol,'rtci did not converge')
            
      if testCase.verbosity
        %%
        % Compute coordinates of the level curves considered points
        r0 = L.G.rW-LY.aW.*cos(L.G.oW); z0 = L.G.zW+LY.aW.*sin(L.G.oW); % analytical
        r1 = L.G.rW-   aW.*cos(L.G.oW); z1 = L.G.zW+   aW.*sin(L.G.oW); % answer

        clf;
        plot(LY.rA,LY.zA,'ok')
        hold on
        contour(L.rx,L.zx,LY.Fx,linspace(0,3,31),'-r');
        contour(L.rx,L.zx,LY.Fx,LY.FB*[1 1],'-k','linewidth',2)
        plot([L.G.rW,r0].',[L.G.zW,z0].','-xm')
        plot([L.G.rW,r1].',[L.G.zW,z1].','-ob')
        plot(L.G.rl,L.G.zl,'-k')
        hold off, axis equal
        title(sprintf('RTCI test: %s, nz=%d',L.P.tokamak,L.nzx))
      end
      
    end
  end
  
  
end
