classdef (SharedTestFixtures={mexm_fixture}) fsgi_test < meq_test
  % Tests for flux surface geometric integrals calculations
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    L; LY;
    rel_tol = 1e-2; % Error is dominated by pq resolution at the edge.
    verbosity = 0;
  end
  
  properties (MethodSetupParameter)
    npq = {32  64};
  end
  
  properties (TestParameter)
    method = {'fsgimex','fsgimexm'};
  end
    
  methods(TestMethodSetup)
    function fsgiTestSetup(testCase,npq)
      L = liu('ana',1,1,'iterq',1,'pq',linspace(0,1,npq).^0.5);

      % calculate axis location and radius of circle a0
      r0 = (L.rx(end)+L.rx(1))*0.5; 
      z0 = 0;
      FA = 0;
      FB = 1;
      rBt = 1.4;
      
      [testCase.L,testCase.LY] = testCase.getCircularEquilibrium(L,r0,z0,FA,FB,rBt);
    end
  end

  methods(Test, TestTags = {'Unit'})
    function test_fsgi(testCase,method)
      
      %%
      L = testCase.L; %#ok<*PROPLC> % for brevity
      LY = testCase.LY;
      
      sIp = sign(LY.Ip);
      
      qaq = LY.aq.^2;
      M1q = qaq*L.M1q; % -da^2/qdq
      M2q = ((L.M2q*qaq).^2./qaq+qaq)./M1q; % ((da^2/2do)^2/a^2+a^2)/(-da^2/qdq)
  
      % contour integrals
      [Q0Q,Q1Q,Q2Q,Q3Q,Q4Q,iqQ] = feval(method,...
        M1q,M2q,LY.rq,LY.irq,LY.rA,LY.FA,LY.FB,...
        sIp*sqrt(LY.dr2FA*LY.dz2FA-LY.drzFA.^2),...
        LY.lX,LY.rB,LY.iTQ,1/L.doq);
      
      tol = testCase.rel_tol;
      testCase.verifyEqual(Q0Q,LY.Q0Q,'RelTol',tol,'Q0Q value from fsgimex does not match analytical value within tolerance');
      testCase.verifyEqual(Q1Q,LY.Q1Q,'RelTol',tol,'Q1Q value from fsgimex does not match analytical value within tolerance');
      testCase.verifyEqual(Q2Q,LY.Q2Q,'RelTol',tol,'Q2Q value from fsgimex does not match analytical value within tolerance');
      testCase.verifyEqual(Q3Q,LY.Q3Q,'RelTol',tol,'Q3Q value from fsgimex does not match analytical value within tolerance');
      testCase.verifyEqual(Q4Q,LY.Q4Q,'RelTol',tol,'Q4Q value from fsgimex does not match analytical value within tolerance');
      testCase.verifyEqual(iqQ,LY.iqQ,'RelTol',tol,'iqQ value from fsgimex does not match analytical value within tolerance');
      
      if testCase.verbosity
        %%
        clf;
        subplot(231)
        plot(L.pQ,   Q0Q);
        plot(L.pQ,LY.Q0Q,'--','Linewidth',2);
        axis([0 1 0.9 1.1]);
        title('Q0Q')
        
        subplot(232)
        plot(L.pQ,   Q1Q);
        plot(L.pQ,LY.Q1Q,'--','Linewidth',2);
        axis([0 1 -1.1 -0.9]);
        title('Q1Q')
        
        subplot(233)
        plot(L.pQ,   Q2Q);
        plot(L.pQ,LY.Q2Q,'--','Linewidth',2);
        axis([0 1 0.9 1.1]);
        title('Q2Q')
        
        subplot(234)
        plot(L.pQ,   Q3Q);
        plot(L.pQ,LY.Q3Q,'--','Linewidth',2);
        axis([0 1 0 100]);
        title('Q3Q')
        
        subplot(235)
        plot(L.pQ,   Q4Q);
        plot(L.pQ,LY.Q4Q,'--','Linewidth',2);
        axis([0 1 0 100]);
        title('Q4Q')
        
        subplot(236)
        plot(L.pQ,   iqQ);
        plot(L.pQ,LY.iqQ,'--','Linewidth',2);
        axis([0 1 -10 -5]);
        title('iqQ')
        drawnow;
      end
      
    end
    
  end
  
  
end
