classdef anamak_all_code_test < meq_test
  % Test all MEQ suite codes on anamak
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbose = 0;
    tokamak = 'ana'
    iterq = 50;
    tol = 3e-5;
  end
  
  properties (TestParameter)
    shot = struct('circular',1,'diverted',2,'elongated',11);
    usecs = struct('True',true,'False',false);
  end
  
  methods(Test,TestTags={'Integration'})
    function test_all(testCase,shot,usecs)
      dt = 2.5e-5;
      t = 0:dt:1e-3;
      
      tok = testCase.tokamak;
      fprintf('\n** Find initial FBT solution **\n')
      
      if usecs, csargs = {'icsint',true,'ilim',3};
      else,     csargs = {'icsint',false};end

      [Lfbt,LXfbt] = fbt(tok,shot,[],...
        'debug',testCase.verbose,'iterq',testCase.iterq,'tol',testCase.tol,...
        csargs{:});
      
      if testCase.verbose
        meqgplot(Lfbt.G); drawnow; % plot geometry
      end
      
      % Solve FBT
      LYfbt = fbtt(Lfbt,LXfbt);
      testCase.assertNotEmpty(LYfbt,'FBT returned empty structure');
      
      %% Use this as starting point for FGS
      fprintf('** FGS forward equilibrium from FBT **\n')
      
      [Lfgs] = fgs(tok,shot,t(1),'debug',testCase.verbose,...
        'iterq',testCase.iterq,csargs{:});
      LXfgs = meqxconvert(Lfbt,LYfbt,Lfgs);
      LYfgs = fgst(Lfgs,LXfgs);
      
      % Compare
      if testCase.verbose
        clf; meqplotQ(Lfgs,LYfbt,LYfgs); drawnow;
      end
      reltol = 1e-4; % Same tolerance as FGS solver
      testCase.verifyEqual(LYfbt.FB,LYfgs.FB,'AbsTol',reltol,'inaccurate FB')
      testCase.verifyEqual(LYfbt.rA,LYfgs.rA,'AbsTol',reltol,'inaccurate rA')
      testCase.verifyEqual(LYfbt.zA,LYfgs.zA,'AbsTol',reltol*LYfbt.rA,'inaccurate rA')
      
      %% And then evolve using FGE...
      fprintf('** Evolve using FGE forward equilibrium evolution **\n')
         
      [Lfge,LXfge] = fge(tok,shot,t,'usepreconditioner',1,...
        'debug',testCase.verbose,'izgrid',(testCase.verbose>0),...
        'iterq',testCase.iterq,csargs{:});
      
      % run time evolution without control
      LYfge = fget(Lfge,LXfge);
      
      if testCase.verbose
        clf;
        meqmovie(Lfge,LYfge,'decimate',5);
      end
      
      if shot==2
        % run time evolution with control
        [~,~,ctrlpar] = meqctrl('setup',[],[],Lfge,LXfge,'dt',dt);
        Lfge.P.ctrlfct  = @meqctrl;
        Lfge.P.ctrlpar = ctrlpar;
        Lfge.P.ctrlpar.debug = 1;
        % slightly perturb one coil current to displace initial equilibrium
        LXfge.Ia(2) = LXfge.Ia(2)*1.01;
        [~] = fget(Lfge,LXfge,'debug',1);
      end
      %% Do time-varying equil. reconstruction using LIUQE for uncontrolled case
      fprintf('** LIUQE reconstruction **\n')
      
      Lliu = liu(tok,shot,t,'pQ',Lfbt.pQ,'Fterr',1e-4,...
        'debug',testCase.verbose,'psichco',1e-6,csargs{:});
      it =1:10:numel(LYfge.t); % select time slices where to do reconstruction
      LYliut_ref = meqxk(LYfge,it);
      LXliut = meqxconvert(Lfge,LYliut_ref,Lliu,true);
      
      % Deselect some sensors
      wBm = ones(Lliu.G.nm,1);wBm(5) = 0;
      wFf = ones(Lliu.G.nf,1);wFf(11) = 0;
      
      LYliut = liut(Lliu,LXliut,'debug',testCase.verbose,...
        'wBm',wBm,'wFf',wFf);
      testCase.assertNotEmpty(LYliut,'LIUQE did not find any equilibrium')
      nt = numel(LYliut.t);
      
      if testCase.verbose
        meqplotQ(Lliu,meqxk(LYliut,nt),meqxk(LYfge,nt))
      end
      
      reltol = 2e-3; % tight tolerances: should be close to equal
      testCase.verifyEqual(LYliut.FB,LYliut_ref.FB,'RelTol',reltol,'inaccurate FB')
      testCase.verifyEqual(LYliut.rA,LYliut_ref.rA,'RelTol',reltol,'inaccurate rA')
      testCase.verifyEqual(LYliut.zA,LYliut_ref.zA,'AbsTol',reltol*LYfbt.rA,'inaccurate zA')
      
      testCase.verifyEqual(LYliut.Bm,LXliut.Bm,'AbsTol',5*reltol*max(abs(LXliut.Bm(:))),'Verification of Bm failed')
      testCase.verifyEqual(LYliut.Ff,LXliut.Ff,'AbsTol',5*reltol*max(abs(LXliut.Ff(:))),'Verification of Ff failed')
    end
  end
end
