classdef integrals_test < meq_test
  % Tests for surface and volume integrals
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    tol = 1e-4;
    L,LY
    Wp,Ft0,Vp
  end
  
  properties (TestParameter)
    shot = struct('diverted',2,'droplet',81,'doublet',82);
    Ip = {-300e3,350e3};
  end
  
  methods(TestClassSetup)
    
    function setup_mexm(testCase)
      % add mexm path
      testCase.applyFixture(mexm_fixture);
    end
    
    function setupLLY(testCase)
      [testCase.L,~,testCase.LY] = fbt('ana',1,[],'pql',0.1);
    end
  end
  
  methods (Test,TestTags={'Unit'})
    function test_vizr_m(testCase)
      L = testCase.L; %#ok<*PROP>
      LY = testCase.LY;
      
      Opy = logical(testCase.LY.Iy);
      
      % test vs m file implementation
      
      [Wp_,Ft0_,Vp_] = vizrmex (LY.Fx,Opy,L.ry,L.iry,LY.rBt,L.drx,L.dzx);
      [Wp ,Ft0, Vp ] = vizrmexm(LY.Fx,Opy,L.ry,L.iry,LY.rBt,L.drx,L.dzx);
      
      % verify
      testCase.verifyEqual(Vp ,Vp_ ,'AbsTol',testCase.tol)
      testCase.verifyEqual(Ft0,Ft0_,'AbsTol',testCase.tol)
      testCase.verifyEqual(Wp ,Wp_ ,'RelTol',testCase.tol)
      
    end
    
    function test_meqintQ(testCase,Ip)
      %%
      L = testCase.L; %#ok<*PROPLC>
      
      %% Fake circular equilibrium
      r0 = 1.001; z0 = 0;
      FA = sign(Ip); FB = 0;
      rBt = 1; bp=1; qA = 1.5;
      [L,LY] = testCase.getCircularEquilibrium(L,r0,z0,FA,FB,rBt);
      
      % plasma current from basis function coefficients
      [~,TpDg,ITpDg] = L.bfct(1,L.bfp,LY.Fx,LY.FA,LY.FB,LY.Opy,L.ry,L.iry);
      [g0g,Ig0g    ] = L.bfct(5,L.bfp,   [],LY.FA,LY.FB);
      
      %% Coefficient scaling
      ag = meqfit3(L,LY.rA,LY.dr2FA,LY.dz2FA,LY.drzFA,...
        L.fPg,L.fTg,TpDg,ITpDg,g0g,Ig0g,rBt,Ip,bp,qA);
      
      [IpQ,WkQ,WpQ,VpQ,FtQ,Ft0Q] = ...
        meqintQ(L,FA,FB,rBt,ag, LY.Fx, LY.Opy, Ip);
      
      [Ip,Wk,Wp,~,~,~,Vp,Ft,Ft0] = ...
        meqint(L.fPg,L.fTg,TpDg,ITpDg,ag,LY.Fx,LY.Opy,L,rBt);
      
      testCase.verifyEqual(Ip, IpQ(end), 'AbsTol',1e-10, 'IpQ final value not equal to Ip');
      testCase.verifyEqual(Wk, WkQ(end), 'WkQ final value not equal to Wk');
      testCase.verifyEqual(Wp, WpQ(end), 'WpQ final value not equal to Wp');
      testCase.verifyEqual(Vp, VpQ(end), 'VpQ final value not equal to Vp');
      testCase.verifyEqual(Ft, FtQ(end), 'FtQ final value not equal to Ft');
      testCase.verifyEqual(Ft0,Ft0Q(end),'Ft0Q final value not equal to Ft0');
      
      % check expected volume dependency
      a = sqrt((LY.rA-LY.rB).^2 + (LY.zA-LY.zB).^2);
      V = L.pQ'.^2 * (2*pi^2*r0*a.^2);
      testCase.verifyEqual(V/Vp,VpQ/Vp,'AbsTol',1e-2,'Unexpected volume result for analytical case');
    end
    
    function test_meqintQ_multiD(testCase,shot)
      % Use directly FBT simulation and results from meqpost
      [~,~,LY] = fbt('ana',shot,[]);
      
      testCase.verifyEqual(LY.IpQ (end,:).', LY.IpD        , 'AbsTol',1e-10, 'IpQ  final value does not match IpD');
      testCase.verifyEqual(LY.VpQ (end,:).', LY.VpD        , 'VpQ  final value does not match VpD');
      testCase.verifyEqual(LY.FtPQ(end,:).', LY.FtD+LY.Ft0D, 'FtpQ final value does not match FtD+Ft0D');
    end
  end
end
