classdef doublet_tests < meq_test
  % Tests involving doublets
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties (TestParameter)
   % Deselect shot 68710@1.41 temp. since FBT eq is the same as 68466@2.43
   tok    = struct('TCV_droplet','tcv', 'ana_droplet','ana','ana_doublet','ana'); % ,'TCV_doublet','tcv'
   shot   = struct('TCV_droplet',68466, 'ana_droplet',   81,'ana_doublet',   82); % ,'TCV_doublet',68710
   t      = struct('TCV_droplet', 2.43, 'ana_droplet',    0,'ana_doublet',    0); % ,'TCV_doublet', 1.41
  end
  
  properties
    verbosity = 0;
    bfp = [1 1]; 
    iterq = 10;
    agcon = {'Ip','bp'};
  end
  
  methods(Test,TestTags={'fge'},ParameterCombination='sequential')
    function fbt_liu_fgs_fge(testCase,tok,shot,t)
          
      [ok,msg] = meq_test.check_tok(tok,shot);
      testCase.assumeTrue(ok,msg);
      
      % ----------------- fbt starting point -------------------------------
      testCase.assumeTrue(~isequal(tok,'tcv') || ~isempty(which('mgp')),'skipping test because mgp is not available')
      [Lfbt,LXfbt] = fbt(tok,shot,t,...
        'izgrid',true,'agfitfct',@meqfit3,'agfitp',[true;false;true],...
        'bfct',@bfabmex,'bfp',testCase.bfp,...
        'selx','X',...
        'debug',testCase.verbosity);
      LXfbt.bpD = [0.25;0.25;0];
      LYfbt = fbtt(Lfbt,LXfbt);
      
      if testCase.verbosity
        clf; meqplotfancy(Lfbt,LYfbt);
      end
      testCase.assertTrue(~isempty(LYfbt) && LYfbt.isconverged,'FBT did not converge')
     
      %% LIUQE from fbt
      nD = 3; % number of domain considered 
      ng = sum(testCase.bfp);
      wag = zeros(ng*nD,1);
      wag((nD-1)*ng+1:end) = Inf; % Exact constraints on third domains ag = 0;
      for method = {'jfnk','picard'}
        inPar = liup_doublet('idoublet',1,...
           'izgrid',true,'bfct',@bfabmex,'bfp',testCase.bfp,...
           'selx','X','debug',testCase.verbosity, 'algoNL', method{1}, 'wag', wag);
        if strcmp(method{1},'picard'), inPar = [inPar,{'stabz',1,'wdz',[0;0;Inf]}];end %#ok<AGROW>
        Lliu = liu(tok,shot,t,inPar{:});
        
        minimal=true; % return minimal LXliu
        LXliu =  meqxconvert(Lfbt, LYfbt, Lliu,minimal);
        switch method{1}
          case 'jfnk',   LYliu = liutnwt(Lliu, LXliu);
          case 'picard', LYliu = liut(Lliu, LXliu);
        end
        
        testCase.assertTrue(~isempty(LYliu) ,['liu from fbt did not converge with method ',method{1}])
        testCase.assertEqual(LYliu.ag , LYfbt.ag             ,'AbsTol', 3e-2*norm(LYfbt.ag ,Inf), ['liu does not recover ag  fbt solution for noise free data with method ',method{1}]);
        testCase.assertEqual(LYliu.Bm , LYfbt.Bm             ,'AbsTol', 3e-2*norm(LYfbt.Bm ,Inf), ['liu does not recover Bm  fbt solution for noise free data with method ',method{1}]);
        testCase.assertEqual(LYliu.IpD, LYfbt.IpD(1:Lliu.nD) ,'AbsTol', 1e-3*norm(LYfbt.IpD,Inf), ['liu does not recover IpD fbt solution for noise free data with method ',method{1}]);
        
        if testCase.verbosity
          clf; meqplotliu(Lliu, LXliu,LYliu);
        end
      end
      
      %% FGS from fbt
      [Lfgs] = fgs(tok,shot,t,'insrc','fbt',...
        'bfct',@bfabmex,'bfp',testCase.bfp,...
        'agcon',testCase.agcon,...
        'izgrid',true,'iterq',testCase.iterq,'idoublet',true);
      LXfgs = meqxconvert(Lfbt,LYfbt,Lfgs);
      LYfgs = fgst(Lfgs,LXfgs,'debug',testCase.verbosity,'mkryl',150,'tolF',1e-6);
      testCase.assertTrue(~isempty(LYfgs) && LYfgs.isconverged,'FGS did not converge')
   
      %% LIUQE from FGS TODO
  
      %% FGE with fixed Ip from FGS 
      t = t + (0:1e-4:5e-3); %#ok<*PROP>
      [Lfge] = fge(tok,shot,t,...
        'bfct',@bfabmex,'bfp',testCase.bfp,...
        'agcon',testCase.agcon,...
        'iterq',testCase.iterq,...
        'idoublet',true,...
        'izgrid',true);
      LXfge = fgex(testCase.tok,t,Lfge,meqxconvert(Lfgs,LYfgs,Lfge));
      LXfge.Va(1,:) = LXfge.Va(1,:) + randn(size(LXfge.t));

      Lfge = fgel(Lfge,meqxk(LXfge,1));
      
      LYfge = fget(Lfge,LXfge,'debug',testCase.verbosity);
      testCase.assertTrue(~isempty(LYfge) && all(LYfge.isconverged),'FGE did not converge')
      
      if testCase.verbosity
        meqplotevo(Lfge,LYfge); drawnow;
      end
      
      %% FGE with CDE from FGS
      [LfgeCDE] = fge(tok,shot,t,...
       'bfct',@bfabmex,'bfp',testCase.bfp,...
       'agcon',testCase.agcon,...
       'iterq',testCase.iterq,...
       'cde','cde0Dss',...
       'ssinit',true,'idoublet',true,... % start with stationary initial condition
       'izgrid',true);
      
      LXfgeCDE = fgex(testCase.tok,t,LfgeCDE,meqxconvert(Lfgs,LYfgs,LfgeCDE));
      Lfge = fgel(Lfge,meqxk(LXfge,1));
      LYfge = fget(Lfge,LXfge,'debug',testCase.verbosity);
      testCase.assertTrue(~isempty(LYfge) && all(LYfge.isconverged),'FGE did not converge')
      
      if testCase.verbosity
       meqplotevo(Lfge,LYfge); drawnow;
      end
      
      if testCase.verbosity
        clf;
        subplot(121);
        meqplotfancy(Lfge ,meqxk(LXfge ,1),'vacuum',true); title('eq. field - no Iv');
        subplot(122);
        meqplotfancy(LfgeCDE,meqxk(LXfgeCDE,1),'vacuum',true); title('eq field - w Iv');
        drawnow;
      end

    end
  end
end
