classdef (SharedTestFixtures={mexm_fixture}) qint_test < meq_test
  % Tests of N-point quadratic interpolation
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties (ClassSetupParameter)
    dr = {0.2 1 10}; % radial grid size
    fz = {0.5 1}; % fraction dz/dr
    n = {6,7,8,9}; % number of interpolation points
  end

  properties (TestParameter)
    method = {'qintmex','qintmexm'};
  end
  
  properties
    verbose = 0;
    nint = 11; % number of interpolation points
    rx, zx;
    nr, nz;
    rri, zzi;
    inM; % parameter interpolation matrix
    tolF = 5e-3; % lower tolerances for F,Br,Bz since Fx not a truly quadratic function
    tolB = 5e-3;
    toldB = 5e-2; % higher more inaccurate
  end
  
  methods (TestClassSetup)
    function generateGrid(testCase,dr,fz,n)
      % Evaluate on non-square grid
      nnr = 13; nnz = 5;
      dz = dr*fz;
      
      testCase.rx = (10+(0:nnr-1))*dr; % R grid for data
      testCase.zx = (-nnz:nnz)*dz;     % Z grid for data
      
      testCase.nr = nnr;
      testCase.nz = 2*nnz+1;
      
      intp.n = n;
      testCase.inM = qintc(intp,dr,dz); % contains interpolation coefficient matrix
      
      % points on which to perform interpolation
      % avoidgrid boundary - these are understandably inaccurate
      ri = linspace(testCase.rx(2),testCase.rx(end-1),testCase.nint);
      zi = linspace(testCase.zx(2),testCase.zx(end-1),testCase.nint);
      [testCase.rri,testCase.zzi] = ...
        meshgrid(ri,zi);
      
    end
  end
  
  methods (Test,TestTags={'Unit'})
    function test_value(testCase, method)
      % test value interpolation vs analytical fit
      
      % A quadratic test function
      H = 0.1*[-1 .1;.1 1]; G = 0.1*[-1;1]; c = 0;
      q = @(x) sum(x.*(H*x.').',2) + x*G + c; % size(q)=[n,1] for size(x) = [n,2] 
      
      % Expressions for ficticious Br,Bz for this quadratic Fx
      qBz = @(x)  (2*(H(1,:)*x.').' + G(1))./x(:,1)/2/pi; % Bz =  1/2piR * dF/dR
      qBr = @(x) -(2*(H(2,:)*x.').' + G(2))./x(:,1)/2/pi; % Br = -1/2piR * dF/dZ
      
      % Expressions for ficticious Brr,Brz,Bzr,Bzz for this quadratic Fx
      qBzz = @(x)  (2*H(1,2))./x(:,1)/2/pi; % Bzz =  1/2piR * dF2/dRdZ
      qBrz = @(x) -(2*H(2,2))./x(:,1)/2/pi; % Brz = -1/2piR * dF2/dZdZ
      qBzr = @(x)  (2*H(1,1))./x(:,1)/2/pi - qBz(x)./x(:,1); % Bzr =  1/2piR * dF2/dRdR - Bz/R
      qBrr = @(x) -(2*H(2,1))./x(:,1)/2/pi - qBr(x)./x(:,1); % Brr = -1/2piR * dF2/dZdR - Br/R
      
      [rrx,zzx] = meshgrid(testCase.rx,testCase.zx);
      Yx = reshape(q([rrx(:),zzx(:)]),size(rrx));
      
      % test interpolation
      % grid of points where to interpolate
      ri = testCase.rri(:);  zi = testCase.zzi(:);
      
      if testCase.verbose
        clf;
        pcolor(testCase.rx,testCase.zx,Yx);
        colorbar;hold on;
        plot(testCase.rri,testCase.zzi,'ok');
        drawnow;
      end
      
      [Ye,Bze,Bre,Brre,Brze,Bzre,Bzze] = deal(zeros(size(testCase.rri)));
      [Yi,Bzi,Bri,Brri,Brzi,Bzri,Bzzi] = deal(zeros(size(testCase.rri)));
      % interpolate
      [Yi(:),Bri(:),Bzi(:),Brri(:),Brzi(:),Bzri(:),Bzzi(:)] = ...
        feval(method,testCase.rx,testCase.zx,Yx,ri,zi,testCase.inM);
      
      % directly evaluate
      Ye(:) = q([ri,zi]); % evaluate
      Bre(:) = qBr([ri,zi]); % evaluate
      Bze(:) = qBz([ri,zi]); % evaluate
      
      Brre(:) = qBrr([ri,zi]); % evaluate
      Brze(:) = qBrz([ri,zi]); % evaluate
      Bzre(:) = qBzr([ri,zi]); % evaluate
      Bzze(:) = qBzz([ri,zi]); % evaluate
      
      testCase.verifyEqual(Yi,Ye,'AbsTol',sqrt(eps));
      testCase.verifyEqual(Bri,Bre,'AbsTol',sqrt(eps));
      testCase.verifyEqual(Bzi,Bze,'AbsTol',sqrt(eps));
      testCase.verifyEqual(Brri,Brre,'AbsTol',sqrt(eps));
      testCase.verifyEqual(Bzri,Bzre,'AbsTol',sqrt(eps));
      testCase.verifyEqual(Brzi,Brze,'AbsTol',sqrt(eps));
      testCase.verifyEqual(Bzzi,Bzze,'AbsTol',sqrt(eps));
      
    end
    
    function test_fields(testCase, method)
      % Test actual flux, field calculation vs greenem
      
      % grid of points on which to generate Fx
      [rrx,zzx] = meshgrid(testCase.rx,testCase.zx);
      
      % a fictictious winding outside the grid
      rw = max(rrx(:))*1.5;  zw = mean(zzx(:))+0.1*max(zzx(:));
      % Flux on this grid using greenem
      R0 = mean(testCase.rx); % scale current by R0 to keep same fields
      Iw = R0;
      Fx = reshape(greenem('mut',rrx,zzx,rw,zw)*Iw,testCase.nz,testCase.nr);
      % Br,Bz by interpolation
      ri = testCase.rri(:);  zi = testCase.zzi(:);
      
      % init
      [Fi,Bri,Bzi,Brri,Brzi,Bzri,Bzzi] = deal(zeros(size(testCase.rri)));
      [Fe,Bre,Bze,Brre,Brze,Bzre,Bzze] = deal(zeros(size(testCase.rri)));

      % Fluxes, fields, field deriatives by evaluation
      Fe(:) = greenem('mut',ri,zi,rw,zw)*Iw; 
      Bre(:) = greenem('br',ri,zi,rw,zw)*Iw;
      Bze(:) = greenem('bz',ri,zi,rw,zw)*Iw;
      Brre(:) = greenem('dr1br',ri,zi,rw,zw)*Iw;
      Bzre(:) = greenem('dr1bz',ri,zi,rw,zw)*Iw;
      Brze(:) = greenem('dz1br',ri,zi,rw,zw)*Iw;
      Bzze(:) = greenem('dz1bz',ri,zi,rw,zw)*Iw;
      
      [Fi(:),Bri(:),Bzi(:),Brri(:),Brzi(:),Bzri(:),Bzzi(:)] = ...
        feval(method,testCase.rx,testCase.zx,Fx,ri,zi,testCase.inM);     % interpolate
      
      if testCase.verbose
        clf;
        subplot(231)
        pcolor(testCase.rx,testCase.zx,Fx);
        colorbar;hold on;
        plot(testCase.rri,testCase.zzi,'ok');
        title('original');
        subplot(234);
        pcolor(testCase.rri,testCase.zzi,Fe); colorbar
        title('interpolated')
        
        subplot(232)
        pcolor(testCase.rri,testCase.zzi,Bre); colorbar
        colorbar;
        title('Bre original');
        subplot(235);
        pcolor(testCase.rri,testCase.zzi,Bri); colorbar
        title('Bre interpolated')
        
        subplot(233)
        pcolor(testCase.rri,testCase.zzi,Bze); colorbar
        colorbar;
        title('Bze original');
        subplot(236);
        pcolor(testCase.rri,testCase.zzi,Bzi); colorbar
        title('Bze interpolated')
        
        drawnow;
      end
      
      Bpe = sqrt(Bre.^2+Bze.^2); % local Bp for normalization
      
      switch size(testCase.inM,2)
        case {6,7,8}
          fact = 5; % more lax test tolerance for these fits
        case {9}
          fact = 1;
        otherwise
          error('invalid inM size');
      end
      tolF = testCase.tolF*fact; %#ok<*PROP>
      tolB = testCase.tolB*fact;
      toldB = testCase.toldB*fact;
      
      testCase.verifyEqual(Fi,Fe,'RelTol',tolF); % lower tolerance since Fx is not a quadratic function!
      testCase.verifyLessThan(abs((Bri-Bre)./Bpe),tolB);
      testCase.verifyLessThan(abs((Bzi-Bze)./Bpe),tolB);
      
      testCase.verifyLessThan(abs((Brri-Brre)./Bpe),toldB);
      testCase.verifyLessThan(abs((Brzi-Brze)./Bpe),toldB);
      testCase.verifyLessThan(abs((Bzri-Bzre)./Bpe),toldB);
      testCase.verifyLessThan(abs((Bzzi-Bzze)./Bpe),toldB);
    end
  end
end

