classdef nfdb_test < meq_test
  % tests of nfdb (normal difference on flux boundary)
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    tol = 1e-2;
    L;
  end
  
  properties (MethodSetupParameter)
    nz = {32,64};
    nr = {27,54};
  end
  
  methods (TestMethodSetup, ParameterCombination='sequential')
    function set_L(testCase,nr,nz)
      testCase.L = liu('ana',1,0,'nr',nr,'nz',nz);
    end
  end
  
  methods (Test, TestTags={'Unit'})
    function test_Lackner_trick(testCase)
      L = testCase.L; %#ok<*PROP>
      
      % Square Iy distribution
      rA = 0.9;
      zA = 0;
      A = 0.1;
      
      Iy = (double(abs(L.zy-zA)<A) * double( abs(L.ry - rA) < A)');
      % contourf(L.ry,L.zy,Iy)
      
      %% Lackner trick: contour integral using normal flux at boundary
      % Flux for this Iy
      Fb = zeros(2*(L.nzx+L.nrx-2),1);
      Fx = gszrmex(Fb,Iy,L.cx,L.cq,L.cr,L.cs,L.ci,L.co,0);
      nf = nfdbmex(Fx); % nf just computes flux normal derivative
      Fb = L.Tbc*nf; % Tbc does the contour integration
      
      % This is the exact same operation done by
      Fbfbp = meqfbp(Iy,L);
      testCase.assertEqual(Fb,Fbfbp,'Fb not equal to meqfbp output')
      
      %% Slow method: direct mutual evaluation (integral over grid)
      [rry,zzy] = meshgrid(L.ry,L.zy);
      Mby = greenem('mut',L.rb  ,L.zb  ,rry , zzy ); % mutuals y->b
      Fb2 = Mby*reshape(Iy,L.ny,1);
      
      testCase.assertLessThan(norm(Fb-Fb2)/norm(Fb),testCase.tol,...
        'Fb calculation using Lackner trick not accurate');
      
    end
  end
end
