classdef fbt_test<meq_test
  % Tests of FBT  
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbosity = 0;
    tokamak = 'ana';
    tol = sqrt(eps);
  end
  
  properties (TestParameter)
    shot = struct('circular',1,'diverted',2,'doublenull',4);
    bp = {0.1 1 2};
    rBt = {1}; % adding -1 makes convergence more difficult, to be investigated
    qA = {0.8 2};
    Ip = {50e3,-150e3};
    bfct = {'bffbt','bf3pmex'};
    margin = {'default','margin','margin_single'};
  end
  
  methods(Test,TestTags={'fbt'})
    function test_scan(testCase,shot,rBt,Ip,bp,qA,bfct)
      
      tok = testCase.tokamak;
      qA = sign(Ip)*sign(rBt)*qA;
      
      [L,LX] = fbt(tok,shot,[],...
        'agfitfct',@meqfit3,...
        'bfct',str2func(bfct),'bfp',[],...
        'debug',0,...
        'iterq',50,'izgrid',testCase.verbosity>0,...
        'agfitp',[true,true,true],...
        'pql',0.2,...
        'ihole',false); % Necessary at high betap for FBT initial iterations in these configurations
     
      LX.IpD = Ip;
      LX.bpD = bp;
      LX.qA  = qA;
      LX.rBt = rBt;
      % Solve FBT
      LY = fbtt(L,LX);
      
      testCase.assertNotEmpty(LY,'LY returned empty');
      testCase.verifyTrue(LY.isconverged,'fbt did not converge')
      testCase.verifyEqual(LY.bp,bp,'RelTol',testCase.tol);
      testCase.verifyEqual(LY.qA,qA,'RelTol',testCase.tol);
      testCase.verifyEqual(LY.Ip,Ip,'RelTol',testCase.tol);
      
      if testCase.verbosity
        subplot(211)
        meqplotfancy(L,LY);
        title(sprintf('bp=%2.1f,qA=%2.1f,Ip=%2.2fkA,rBt=%2.1fTm',LY.bp,LY.qA,LY.Ip/1e3,LY.rBt))
        subplot(223)
        plot(L.pQ.^2,LY.PQ/1e3);
        title('P [kPa]'); xlabel('\psi_N');
        ylim([0 20]);

        subplot(224)
        plot(L.pQ.^2,1./LY.iqQ);
        title('q'); xlabel('\psi_N');
        if sign(Ip)*sign(rBt)>0, ylim([0 7]); else, ylim([-7 0]); end
        drawnow;
      end
    end
        
    function test_fbtxdisp(testCase,shot)
      [L,LX] = fbt('ana',shot);
      
      % Call without output
      testCase.assertWarningFree(@()fbtxdisp(L,LX));
      
      % Call with output
      s = fbtxdisp(L,LX);
      testCase.assertClass(s,'cell');
    end
    
    function test_limm(testCase,shot,margin)
      % Test margin factor parameter limm for coil current combinations
      tok = testCase.tokamak;
      
      L = fbt(tok,shot); na = L.G.na;
      switch margin
        case 'default', limm=1;
        case 'margin', limm=0.5;
        case 'margin_single', limm=ones(na,1); limm(4)=0.1;
      end
      
      [L,~,LY] = fbt(tok,shot,[],'limu', 5e5*ones(na,1),...
                                 'liml',-5e5*ones(na,1),...
                                 'limm',limm,...
                                 'debug',    testCase.verbosity,...
                                 'izgrid',testCase.verbosity);
      
      % Verify FBT converged and solution is within prescribed limits
      testCase.assertNotEmpty(LY,'LY returned empty');
      testCase.verifyTrue(LY.isconverged,'fbt did not converge')
      testCase.verifyLessThan( L.P.limc*LY.Ia, L.P.limm.*L.P.limu + testCase.tol);
      testCase.verifyLessThan(-L.P.limc*LY.Ia,-L.P.limm.*L.P.liml + testCase.tol);
      
      if testCase.verbosity
        clf;
        subplot(121); meqplotfancy(L,LY);
        title(sprintf('limm=%s',margin),'interpreter','none');
        
        subplot(222); bar( L.P.limc*LY.Ia/1e3); hold on; plot(1:na,L.P.limm.*L.P.limu/1e3,'x');
        title('upper limit [kA]')

        subplot(224); bar(-L.P.limc*LY.Ia/1e3); hold on; plot(1:na,-L.P.limm.*L.P.liml/1e3,'x'); 
        title('lower limit [kA]')

        drawnow;
      end                           
    end
    
    function test_novessel(testCase)
      % fbt should work without a vessel definition at all
      
      % Default run should give empty Iu and zero Iv
      myshot = 1; t=0; tok='ana';
      [L,~,LY] = fbt(tok,myshot,t);
      testCase.assertEqual(LY.Iv,zeros(L.G.nv,numel(t)))
      testCase.assertEqual(LY.Iu,zeros(     0,numel(t)))

      % Run without vessel should give empty Iv
      [~,~,LY] = fbt(tok,myshot,t,'selu','n','nv',0);
      testCase.assertEqual(LY.Iv,zeros(0,numel(t)));
      testCase.assertEqual(LY.Iu,zeros(0,numel(t)));
    end
    
    function test_fbtx(testCase)
      tol_  = 1e-6;
      eqtol = 1e-11;
      %% Test for reordering
      [L,LX] = fbt('ana',1,0,'tol',tol_);
      LX = fbtgp(LX,1,1,0,[],[],[],NaN,NaN,NaN,0,[],[],[],[],[],[],[],[]);
      LX = fbtgp(LX,1,2,0,[],[],[],NaN,NaN,NaN,0,[],[],[],[],[],[],[],[]);
      LX = fbtgp(LX,1,1,0,[],[],[],NaN,NaN,NaN,0,[],[],[],[],[],[],[],[]);
      LY = fbtt(L,LX); % baseline
      %LX = fbtx(L,LX);
      % check reordering
      testCase.verifyEqual(LY,fbtt(L,LX),'AbsTol',1e-12,...
        'fbt solutions with/without reordering are not the same');
      
      %% Test for eliminating redundant constraints
      
      % If separate Br=0 Bz=0, merge and eliminate the second one
      [L,LX] = fbt('ana',1,0,'tol',tol_); % br bz ba be
      LX = fbtgp(LX,1,1,0,[],[],[],        0  ,NaN,NaN,0,[],[],[],[],[],[],[],[]);
      LX = fbtgp(LX,1,1,0,[],[],[],        NaN,0  ,NaN,0,[],[],[],[],[],[],[],[]);
      LX.gpbd = 1;
      LY = fbtt(L,LX); % baseline
      LX = fbtx(L,LX);
      % expect
      testCase.assertEqual(structcmp(LY,fbtt(L,LX),eqtol,{'chie','res'}),true,'fbt solutions with/without grouping are not the same');
      
      % If redundant Bz=0, eliminate the second one
      [L,LX] = fbt('ana',1,0,'tol',tol_); % br bz  ba be
      LX = fbtgp(LX,1,1,0,[],[],[],        0  ,0  ,NaN,0,[],[],[],[],[],[],[],[]);
      LX = fbtgp(LX,1,1,0,[],[],[],        NaN,0  ,NaN,0,[],[],[],[],[],[],[],[]);
      LX.gpbd = 1;
      LY  = fbtt(L,LX); % baseline
      LX = fbtx(L,LX);
      % expect
      testCase.assertEqual(structcmp(LY,fbtt(L,LX),eqtol,{'chie','res'}),true,'fbt solutions with/without grouping are not the same');
      
      % If not a constraint, don't eliminate it
      [L,LX] = fbt('ana',1,0,'tol',tol_); % br bz   ba  be
      LX = fbtgp(LX,1,1,1,[],[],[],        0  ,NaN ,NaN,1,[],[],[],[],[],[],[],[]);
      LX = fbtgp(LX,1,1,1,[],[],[],        NaN,0   ,NaN,0,[],[],[],[],[],[],[],[]);
      LX.gpbd = 1;
      % expect
      testCase.verifyEqual(LX,fbtx(L,LX),'LX changed while it should not have');
      
      % if not at same R,Z, don't eliminate it
      [L,LX] = fbt('ana',1,0,'tol',tol_); % br  bz  ba be
      LX = fbtgp(LX,2,1,0,[],[],[],        0  ,NaN,NaN,0,[],[],[],[],[],[],[],[]);
      LX = fbtgp(LX,1,1,0,[],[],[],        NaN,0  ,NaN,0,[],[],[],[],[],[],[],[]);
      LX.gpbd = 1;
      testCase.verifyEqual(LX,fbtx(L,LX),'LX changed while it should not have');
      
    end
    
  end
end
