classdef greenem_test < genlib_test
  % Tests for greenem
  %
  % [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties (TestParameter)
    mode = {'s','m','br','bz',...
      'dr1m','dz1m','dr2m','dz2m',...
      'dr1r1m','dr1r2m','dr1z1m','dr1z2m',...
      'dr2r1m','dr2r2m','dr2z1m','dr2z2m',...
      'dz1r1m','dz1r2m','dz1z1m','dz1z2m',...
      'dz2r1m','dz2r2m','dz2z1m','dz2z2m',...
      'dbrdr' 'dbrdz','dbrdr2','dbrdz2',...
      'dbzdr','dbzdr1','dbzdz','dbzdz1',...
      'dbzdr2','dbzdz2','zero'}
    mode2 = {'m','dr2m','dz2m'}
  end
  
  properties
    verbosity = false;
  end
  
  methods (Test)
    function test_greenem(testCase,mode)
      % Verify all modes are returning correctly for simple inputs
      r1 = 1;
      z1 = 0;
      r2 = 1;
      z2 = 1;
      y = greenem(mode,r1,z1,r2,z2);
      testCase.verifyFalse(any(isnan(y)),'greenem returned NaN');
    end
    
    function test_greenem_close(testCase,mode)
      % Verify that greenem still returns finite values even for close pair of loops.
      d  = logspace(-16,-1,16);
      t  = 0.3;
      r1 = 1;
      z1 = 0;
      r2 = r1 + d*cos(t);
      z2 = z1 + d*sin(t);
      y = greenem(mode,r1,z1,r2,z2);
      testCase.assertFalse(any(isnan(y)),'greenem returned NaN for close pair of loops');
    end
    
    function test_greenem_ellipke(testCase,mode2)
      % Compare mutual values from greenem with computation using ellipke
      dr = 1e-6;
      r1 = 1;
      z1 = 0;
      dd = 0.01;
      r2min =  0.3; r2max =  4.0;
      z2min = -4.0; z2max =  4.0;
      [r2,z2] = meshgrid(r2min:dd:r2max,z2min:dd:z2max);
      r2_ = reshape(r2,1,[]);
      z2_ = reshape(z2,1,[]);
      
      mask = (r2-r1).^2 + (z2-z1).^2 > 0;
      
      switch mode2
        case 'm'
          M1 = greenem('m',    r1,z1,r2_,z2_);
          M2 = greenem_ellipke(r1,z1,r2_,z2_);
      
          testCase.verifyEqual(M1(mask),M2(mask),'RelTol',1e-5,'Failed check for mut');
          
        case 'dr2m'
          M1 = greenem('dr2m',    r1,z1,r2_   ,z2_);
          M2 = (  greenem_ellipke(r1,z1,r2_+dr,z2_) ...
                - greenem_ellipke(r1,z1,r2_-dr,z2_))/2/dr;
          
          testCase.verifyEqual(M1(mask),M2(mask),'AbsTol',1e-5,'Failed check for dr2m');
     
        case 'dz2m'
          M1 = greenem('dz2m',    r1,z1,r2_,z2_   );
          M2 = (  greenem_ellipke(r1,z1,r2_,z2_+dr) ...
                - greenem_ellipke(r1,z1,r2_,z2_-dr))/2/dr;
          
          testCase.verifyEqual(M1(mask),M2(mask),'AbsTol',1e-5,'Failed check for dz2m');
          
      end
      
      if(testCase.verbosity)
        figure;
        subplot(1,2,1);
        contourf(r2,z2,reshape(log10(abs(M1./M2-1)),size(r2)),41,'Linestyle','none');
        title(mode2);
        colorbar;
        caxis([-8,-3]);
        axis equal;
        subplot(1,2,2);
        contourf(r2,z2,reshape(log10(abs(M1)),size(r2)),41,'Linestyle','none');
        title(mode2);
        colorbar;
        axis equal;
      end
    end

    function test_greenem_class(testCase)
      % Verify greenem correctly handles different float types
      r1 = 1;
      z1 = 0;
      r2 = 1;
      z2 = 1;
      yd = greenem('m',double(r1),double(z1),double(r2),double(z2));
      ys = greenem('m',single(r1),single(z1),single(r2),single(z2));
      testCase.verifyEqual(yd,ys,'RelTol',eps('single'),'Double and single precision results do not match');
    end
  end
end

function [M,m] = greenem_ellipke(r1,z1,r2,z2)
   mu0 = 4e-7*pi;
   
   r   = r1+r2;
   h   = z1-z2;
   s   = r.*r+h.*h;
   is  = 1./s;
   m   = 4*r1.*r2.*is;
   [K,E] = ellipke(m);
   v0  = (1.0-0.5*m).*K-E;
   u0  = sqrt(s);
   
   M = mu0*u0.*v0;

end
