classdef analytical_map_test < meq_test
  % Basic tests involving analytical flux maps
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    x,y,xx,yy,dx,dy,xl,yl,fl,xs,ys,fs,dx2fs,dy2fs,dxyfs,xa,ya,fa,dx2fa,dy2fa,dxyfa,F;
    tolxy    = 5e-3;
    tolf    = 1e-2;
    told2f = 5e-2;
    verbosity = 0;
  end
  
  properties (MethodSetupParameter)
    NY = {2^6+1};
    NX = {2^5-1}; % This needs to be 2^k-1 for asxy to work
    limiter = {'circular','box'}
  end
  
  methods(TestMethodSetup)
    function analyticalTestSetup(testCase,NX,NY,limiter)
      
      % Use a 3rd order polyinomal to generate a synthetic map
      a = 1; b = 1; c = 0; d = 1; e = 0; f = 0;
      f = @(x,y) -(a/3*y.^3 + b/2*x.^2 + c*x.*y + d/2*y.^2 + e*x + f*y);
      f2x = @(x,y) -(b);
      f2y = @(x,y) -(2*a*y + d);
      fxy = @(x,y) -(c);
      
      % Define grid
      testCase.x = linspace(-1,1,NX);
      testCase.y = linspace(-1.5,1,NY-1);
      testCase.dx = testCase.x(2)-testCase.x(1);
      testCase.dy = testCase.y(2)-testCase.y(1);
      [testCase.xx,testCase.yy] = meshgrid(testCase.x,testCase.y);
      
      % Axis and X-point locations
      testCase.xs = 0; 
      testCase.ys = -1; 
      testCase.fs = f(testCase.xs,testCase.ys);
      testCase.xa = 0; 
      testCase.ya =  0;  
      testCase.fa = f(testCase.xa,testCase.ya);
      testCase.dx2fs = f2x(testCase.xs,testCase.ys);
      testCase.dy2fs = f2y(testCase.xs,testCase.ys);
      testCase.dxyfs = fxy(testCase.xs,testCase.ys);
      testCase.dx2fa = f2x(testCase.xa,testCase.ya);
      testCase.dy2fa = f2y(testCase.xa,testCase.ya);
      testCase.dxyfa = fxy(testCase.xa,testCase.ya);
      
      % Compute synthetic map over the grid
      testCase.F = testCase.xx*0;
      for ix = 1 : length(testCase.x)
        for iy = 1 : length(testCase.y)
          testCase.F(iy,ix) = f(testCase.xx(iy,ix),testCase.yy(iy,ix));
        end
      end
      
      % Find a fake limiter
      switch limiter
        case 'box'
          NL = 20;
          testCase.xl = [linspace(-0.8,0.8,NL) 0.8*ones(1,NL) linspace(0.8,-0.8,NL) -0.8*ones(1,NL)];
          testCase.yl = [-1.2*ones(1,NL) linspace(-1.2,0.6,NL) 0.6*ones(1,NL) linspace(0.6,-1.2,NL)];
        case 'circular'
          NL = 100;
          r = abs(testCase.ys-testCase.ya)*0.75;
          t = linspace(0,2*pi,NL+1); t = t(1:end-1);
          testCase.xl = testCase.xa + r*cos(t);
          testCase.yl = testCase.ya + r*sin(t);
      end
      testCase.fl = f(testCase.xl,testCase.yl);
    end
  end
  
  
  methods(Test, TestTags = {'Unit'})
    
    function test_asxy(testCase)
      
      % test asxy
      [XA,YA,FA,DX2FA,DY2FA,DXYFA,XS,YS,FS,DX2FS,DY2FS,DXYFS] = asxymex(testCase.F',testCase.x,testCase.y,2,...
        testCase.dx,testCase.dy,1/testCase.dx,1/testCase.dy,true(size(testCase.F')));
      
      testCase.verifyLessThan(max(abs(XA-testCase.xa)),testCase.tolxy,'XA asxy result failure')
      testCase.verifyLessThan(max(abs(YA-testCase.ya)),testCase.tolxy,'YA asxy result failure')
      testCase.verifyLessThan(max(abs(FA-testCase.fa)),testCase.tolf,'FA asxy result failure')
      
      testCase.verifyLessThan(max(abs(XS-testCase.xs)),testCase.tolxy,'XS asxy result failure')
      testCase.verifyLessThan(max(abs(YS-testCase.ys)),testCase.tolxy,'YS asxy result failure')
      testCase.verifyLessThan(max(abs(FS-testCase.fs)),testCase.tolf,'FS asxy result failure')
      
      testCase.verifyLessThan(max(abs(DX2FA-testCase.dx2fa)),testCase.told2f,'DX2FA asxy result failure')
      testCase.verifyLessThan(max(abs(DX2FS-testCase.dx2fs)),testCase.told2f,'DX2FS asxy result failure')
      testCase.verifyLessThan(max(abs(DY2FA-testCase.dy2fa)),testCase.told2f,'DY2FA asxy result failure')
      testCase.verifyLessThan(max(abs(DY2FS-testCase.dy2fs)),testCase.told2f,'DY2FS asxy result failure')
      testCase.verifyLessThan(max(abs(DXYFA-testCase.dxyfa)),testCase.told2f,'DXYFA asxy result failure')
      testCase.verifyLessThan(max(abs(DXYFS-testCase.dxyfs)),testCase.told2f,'DXYFS asxy result failure')
      
      if testCase.verbosity
        figure
        subplot(121)
        contour(testCase.xx,testCase.yy,testCase.F,100)
        hold on, grid on, axis equal, orient tall
        plot(testCase.xl,testCase.yl,'-k','linewidth',2)
        plot(testCase.xa,testCase.ya,'ok',testCase.xs,testCase.ys,'ok')
        plot(XA,YA,'xr',XS,YS,'xr','markersize',10)
      end
    
      % test bavx/baax
      OXY = bavxmex(XS,YS,XA-XS,YA-YS,0,testCase.xl,testCase.yl);
      
      % If x-point is inside the limiter, the X-point polygon should be different from the limiter polygon
      if inpolygon(XS,YS,testCase.xl,testCase.yl)
        testCase.verifyNotEqual(testCase.xl,testCase.xl(OXY),'bavx test result failure')
      else
        testCase.verifyEqual(testCase.xl,testCase.xl(OXY),'bavx test result failure')
      end
          
      % test fbnd
      [FB,XB,YB,LB,LX] = fbndmex(testCase.fl(OXY),testCase.xl(OXY),testCase.yl(OXY),FS,XS,YS,-1);
      
      % Return false if no boundary defining point is found
      testCase.assertTrue(LB, 'Boundary defining point search failure');
      
      if testCase.verbosity
        contour(testCase.xx,testCase.yy,testCase.F,[FB FB],'linewidth',2,'color','k')
        plot(XB,YB,'^g')
        if LX
          plot([XB testCase.xl(OXY) XB],[YB testCase.yl(OXY) YB],'.-r')
        else
          plot([testCase.xl(OXY) testCase.xl(OXY(1))],[testCase.yl(OXY) testCase.yl(OXY(1))],'.r')
        end
      end

      %% test pddom
      X  = testCase.xx(2:end-1,2:end-1);
      Y  = testCase.yy(2:end-1,2:end-1);
      FF = testCase.F(2:end-1,2:end-1);
      
      OPY = pdommex(testCase.F,FB,1,bavxmex(XS,YS,XA-XS,YA-YS,0,X,Y));
      
      % Return false if the flux inside the plasma boundary is smaller than the flux at the boundary
      testCase.verifyGreaterThanOrEqual(FF(OPY),FB,'pdom test failure');
      
      if testCase.verbosity
        plot(X(OPY),Y(OPY),'.b','markersize',0.8)
      end
    end
    
    function test_bbox(testCase)
      
      %% test bbox
      O1 = inpolygon(testCase.xx,testCase.yy,testCase.xl,testCase.yl);
      B  = bboxmex(O1',testCase.x,testCase.y);
      O2 = inpolygon(testCase.xx,testCase.yy,B([1 3 3 1]),B([2 2 4 4]));
      
      % Verify that all the points in the limiter are also inside the box and that no points outside the limiter are outside the box
      testCase.verifyEqual(O1,O1&O2,'bbox test failure');
      testCase.verifyEqual(O2,O1|O2,'bbox test failure');
      
      if testCase.verbosity
        subplot(122)
        plot(testCase.xx,testCase.yy,'.k','markersize',0.8)
        hold on, plot(testCase.xl,testCase.yl,'.r',B([1 3 3 1 1]),B([2 2 4 4 2]),'-sb')
        axis equal, orient tall, grid on
        xlim([min(testCase.x),max(testCase.x)])
        ylim([min(testCase.y),max(testCase.y)])
        title('bbox test')
        shg
      end
    end
  end
end
