classdef (SharedTestFixtures={mexm_fixture}) ...
    shap_test < meq_test
  % test of shape parameter calculations
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    R0 = 1;
    nth = 8;
    verbosity = 0;
    tol = 100*eps;
  end
  
  properties(TestParameter)
    cappa = {1.0,1.2,2.0};
    a = {0.2};
    Z0 = {0};
    nrho = {1,2,7};
    xpt = {[0;0],[1;-1],[1;1],[-1;1]}; % x point
    shot = struct('Snowflake',40000,'SNL',61400);
    t= struct('Snowflake',1,'SNL',1);
  end
  
  methods (TestClassSetup)
    %  function compile(testCase) %#ok<MANU>
    %    %% Compile - OSX case kept for debugging
    %    [s,w]=system('gcc -c libmeq/shap.c -o obj/dshap.o -I. -fPIC -Wno-extra-tokens -I -O3 -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk');
    %    disp(w);
    %    if s;error('compilation failed\n%s',w); end
    %    mex CFLAGS='$CFLAGS -fPIC -I. -I$CPATH' -output shapmex.mexmaci64 mexc/shapmex.c obj/dshap.o -ldl /opt/intel/compilers_and_libraries_2019.5.281/mac/mkl/lib/libmkl_intel_ilp64.a /opt/intel/compilers_and_libraries_2019.5.281/mac/mkl/lib/libmkl_intel_thread.a /opt/intel/compilers_and_libraries_2019.5.281/mac/mkl/lib/libmkl_core.a -liomp5 -lpthread -lm
    %  end
   
  end
  
  methods (Static)
    function printresults(m)
      fprintf('\n rgeom, zgeom, amino, epsil, kappa, delta, deltl, deltu\n')
      for ii=1:size(m.rgeom,2)
        fprintf('%6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\n',...
          m.rgeom(ii),m.zgeom(ii),m.amino(ii),m.epsil(ii),m.kappa(ii),m.delta(ii),m.deltl(ii),m.deltu(ii))
      end
    end
  end
  
  methods (Test,TestTags={'Unit'})
    
    function test_shap_analytical(testCase,Z0,a,nrho,cappa,xpt)
      
      %%
      tgr = 2*pi*(1:testCase.nth)'/testCase.nth;
      agr  = a*(1:nrho)/nrho;
      
      rq = testCase.R0 + cos(tgr)*agr;
      zq = Z0 + cappa*sin(tgr)*agr;
      rA = testCase.R0;
      zA = Z0;
      dr2FA = 1;
      dz2FA = 1;
      dr = 0; % minimum grid size distance from boundary point to reject contour point
      dz = 0;
      
      if all(xpt==0)
        rB = [];
        zB = [];
      else
        % add x point at various locations
        rB = testCase.R0+xpt(1)*1.1*a;
        zB = Z0+xpt(2)*1.1*cappa*a;
      end
      
      %% meqshape (matlab version)
      [rmat.rgeom,rmat.zgeom,rmat.amino,rmat.epsil,rmat.kappa,rmat.delta,rmat.deltl,rmat.deltu,...
        rmat.rrmax,rmat.zrmax,rmat.rrmin,rmat.zrmin,rmat.rzmax,rmat.zzmax,rmat.rzmin,rmat.zzmin] ...
        = shapmexm(rq,zq,rB,zB,rA,zA,dr2FA,dz2FA,dr,dz,testCase.verbosity);
      
      if testCase.verbosity
        shap_test.printresults(rmat)
        clf;
        
        plot(rq,zq,'+k'); hold on;
        plot(rB,zB,'xk','markersize',12);
        plot(rmat.rrmin,rmat.zrmin,'rs',rmat.rrmax,rmat.zrmax,'rs','markersize',15);
        plot(rmat.rzmin,rmat.zzmin,'ro',rmat.rzmax,rmat.zzmax,'ro','markersize',15);
        plot(rmat.rgeom,rmat.zgeom,'co');
        axis equal;
      end
      %% Run mex version
      [rmex.rgeom,rmex.zgeom,rmex.amino,rmex.epsil,rmex.kappa,rmex.delta,rmex.deltl,rmex.deltu,...
        rmex.rrmax,rmex.zrmax,rmex.rrmin,rmex.zrmin,rmex.rzmax,rmex.zzmax,rmex.rzmin,rmex.zzmin] = ...
        shapmex(rq,zq,rB,zB,rA,zA,dr2FA,dz2FA,0,0);
      
      if testCase.verbosity
        shap_test.printresults(rmex)
        
        plot(rmex.rrmin,rmex.zrmin,'bs',rmex.rrmax,rmex.zrmax,'bs','markersize',12);
        plot(rmex.rzmin,rmex.zzmin,'bo',rmex.rzmax,rmex.zzmax,'bo','markersize',12);
        plot(rmex.rgeom,rmex.zgeom,'co');
        drawnow
      end
      
      %% Compare
      fields = {'rgeom','zgeom','amino','epsil','kappa','delta','deltl','deltu',...
        'rrmax','zrmax','rrmin','zrmin','rzmax','zzmax','rzmin','zzmin'};
      
      for myfield=fields
        testCase.verifyEqual(rmat.(myfield{1}),rmex.(myfield{1}),'AbsTol',testCase.tol,sprintf('Fields %s do not match',myfield{1}));
      end
    end
  end
  
  methods (Test, TestTags={'TCV'}, ParameterCombination='sequential')
    function test_shape_TCV(testCase,shot,t)
      [L,~,LY] = liuqe(shot,t,'iterq',50); % run liuqe and force contour processing
      %%
      tic
      N=100;
      for ii=1:N
        [s.rgeom,s.zgeom,s.amino,s.epsil,s.kappa,s.delta,s.deltl,s.deltu,...
          s.rrmax,s.zrmax,s.rrmin,s.zrmin,s.rzmax,s.zzmax,s.rzmin,s.zzmin] = ...
          shapmex(LY.rq,LY.zq,LY.rB,LY.zB,LY.rA,LY.zA,LY.dr2FA,LY.dz2FA,L.drx,L.dzx);
      end
      fprintf('shapmex in %3.3fus\n',toc/N*1e6);
      
      if testCase.verbosity
        clf;
        contourf(L.rx,L.zx,LY.Fx-LY.FB,51); hold on;
        contour(L.rx,L.zx,LY.Fx-LY.FB,[0 0],'k','linewidth',2); hold on;
        axis equal tight;
        plot(LY.rq,LY.zq,'.'); hold on;
        plot(s.rrmin,s.zrmin,'ks',s.rrmax,s.zrmax,'ks','markersize',15);
        plot(s.rzmin,s.zzmin,'ko',s.rzmax,s.zzmax,'ko','markersize',15);
        plot(s.rgeom,s.zgeom,'co');
        title(sprintf('TCV %d//%3.3f',shot,t));
        
        shap_test.printresults(s)
      end
      %%
      
      tic
      for ii=1:N
        [m.rgeom,m.zgeom,m.amino,m.epsil,m.kappa,m.delta,m.deltl,m.deltu,...
          m.rrmax,m.zrmax,m.rrmin,m.zrmin,m.rzmax,m.zzmax,m.rzmin,m.zzmin] = ...
          shapmexm(LY.rq,LY.zq,LY.rB,LY.zB,LY.rA,LY.zA,LY.dr2FA,LY.dz2FA,L.drx,L.dzx,testCase.verbosity);
      end
      fprintf('shapmexm in %3.3fus\n',toc/N*1e6);
      
      if testCase.verbosity
        plot(m.rrmin,m.zrmin,'ws',m.rrmax,m.zrmax,'ws','markersize',12);
        plot(m.rzmin,m.zzmin,'wo',m.rzmax,m.zzmax,'wo','markersize',12);
        plot(m.rgeom,m.zgeom,'co');
        
        shap_test.printresults(m)
      end
    end
  end
end


