classdef fge_vacuum_tests < meq_test
  % Tests of vacuum running case of fge
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbosity = 0;
    tol = 1e-10;
  end
  
  methods (Test,TestTags={'fgelin'})
    
    function test_vs_ODE(testCase)
      
      % Tests for FGE vacuum solver
      
      Ts = 1e-3;
      t = (-0.3:Ts:-0.2);
      [L,LX,LY] = fge('tcv',66998,t,'insrc','mag',...
        'selu','e','nu',30);
      
      %% Solve using direct ODE
      A = (L.Mee + Ts*diag(L.Re))\(L.Mee);
      B = (L.Mee + Ts*diag(L.Re))\(Ts*eye(L.ne,L.G.na));
      
      nt = numel(LX.t);
      ia = 1:L.G.na; iu = L.G.na+(1:L.G.nu);
      X = [LX.Ia(:,1);LX.Iu(:,1)]; u = LX.Va;
      
      for kt = 2:nt
        X(:,kt) = A*X(:,kt-1) + B*u(:,kt);
      end
      
      if testCase.verbosity
        subplot(211)
        plot(t,X(ia,:),'b',t,LY.Ia,'r--')
        subplot(212)
        plot(t,X(ia,:)-LY.Ia)
      end
      
      testCase.verifyEqual(X(ia,:),LY.Ia,'AbsTol',testCase.tol);
      testCase.verifyEqual(X(iu,:),LY.Iu,'AbsTol',testCase.tol);
    end
  end
end
