classdef fbttcv_test < meq_test
  % Tests of FBT for TCV
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    % Should be enough to account for differences in models
    tolI = 5e-2;
    tolF = 5e-2;
    tolq = 5e-2;
    verbosity = 0; % optional plots
  end

  properties (TestParameter)
    shot = {64988};
    sip = struct('positive',1,'negative',-1);
    sb0 = struct('positive',1,'negative',-1); % Check qA constraint works with all b0 signs
    
    % For vacuum field constraints
    vrz = struct('none',0,'pos',0.1,'neg',-0.05); % value for dBrdz
    ve  = struct('none',NaN','pos',0.001,'neg',0.001); % weight for dBrdz constraint
    
    % Comparison with stored runs
    mds_shot={70153,65414,68225}; % Only problematic runs so far
  end
  
  methods(TestClassSetup)
    function connect_to_tcvdata(testCase)
      server = 'tcvdata.epfl.ch';
      status = mdsconnect(server);
      testCase.assertTrue(isodd(status),'mdsconnect to tcvdata.epfl.ch failed');
    end
  end
  
  methods(Test, TestTags = {'TCV'})
    function test_fbt_shot(testCase,shot,sip,sb0)
      %% Tests for FBT
      [L,LX,LY] = fbt('TCV',shot,[],'iohfb',sip,'if36fb',sb0);
      
      % Check fbtx displaying option for all time slices (output to void)
      for it=1:numel(LX.t)
        [~] = fbtxdisp(L,meqxk(LX,it));
      end
      
      % Check qA constraint is verified
      qA = 2*pi./LY.rA./sign(LY.FA)./sqrt(LY.dr2FA.*LY.dz2FA-LY.drzFA.^2) .* LY.TQ(1,:);
      testCase.verifyEqual(qA,L.P.qA,'relTol',testCase.tolq);

      % compare with results of FORTRAN version for existing TCV shots
      mdsopen('pcs',shot);
      sip0 = mdsvalue('\mgams.data:iohfb');
      sb00 = mdsvalue('\mgams.data:if36fb');
      % "RAMP" always has as many times as LY.t
      Ia  = sip*sip0*mdsvalue('tcv_eq("i_pol"   ,"ramp")'); % FBT  nodes, with correct signs through tcv_eq
      FAB = sip*sip0*mdsvalue('tcv_eq("psi_axis","ramp")'); % FBTE nodes, has correct signs
      Ia  = Ia(1:16,:);
      FAB = reshape(FAB,1,[]);
      
      testCase.verifyEqual(LY.Ia        , Ia  ,'absTol',testCase.tolI*repmat(max(abs(Ia),[],1),size(Ia,1),1));
      testCase.verifyEqual(LY.FA - LY.FB, FAB ,'relTol',testCase.tolF);
    end
    
    function test_fbt_grid(testCase,shot)
      %% Tests for different grids in FBTE
      [~,~,LY1] = fbt('TCV',shot);
      [~,~,LY2] = fbt('TCV',shot,[],'selx','X');
      
      testCase.verifyEqual(LY2.Ia, LY1.Ia, 'absTol', testCase.tolI*repmat(max(abs(LY1.Ia),[],1),size(LY1.Ia,1),1));
    end

    function test_ia_constraint(testCase,shot)
      %% Tests that Ia constraint using gpia is working
      P = mgp(shot);
      P.gpia(16,20) = -6e3;
      P.gpie(:,20) = [Inf*ones(15,1);0];
      P.gpid(20) = 1;
      [~,~,LY] = fbt('TCV',shot,[],'gpia',P.gpia,'gpie',P.gpie,'gpid',P.gpid);
      % Verify that exact constraint is enforced
      testCase.verifyEqual(P.iohfb*P.gpia(16,20),LY.Ia(16,20),'absTol',1);
    end
    
    function test_ia_conflict(testCase,shot)      
      %% Tests that conflicts between gpia and icoilon are identified
      P = mgp(shot);
      P.gpia(16,end) = -6e3;
      P.gpie(:,end) = [Inf*ones(15,1);0];
      P.gpid(end) = 1;
      % There is no way to specify the number of return arguments for the call in R2017a (fine for this test)
      testCase.verifyError(@() fbt('TCV',shot,[],'gpia',P.gpia,'gpie',P.gpie,'gpid',P.gpid),'fbtptcv:IaConflict');
    end
  end
  
  methods(Test, TestTags = {'TCV'},ParameterCombination='sequential')

    function test_vrz(testCase,vrz,ve)
      %% Tests enforcing vacuum field curvature
      P = mgp(100000);
      
      gpr = P.gpr;
      gpz = P.gpz;
      gpb = P.gpb;
      [gpvrr,gpvrz,gpvzz,gpve] = deal(gpr);
      gpvd = P.gpbd;
      gpr(1) = P.rmajo1(1);
      gpz(1) = P.zmajo1(1);
      gpb(1) = 0;
      gpvrr(1) = NaN;
      gpvrz(1) = vrz; % Maximum value where convergence is obtained
      gpvzz(1) = NaN;
      gpve(1) = ve; % NaN to disable, 0 for exact constraint, 0.001 is ok, 0.1 too high to get desired value
      gpvd(1) = 1;
      [L,~,LY] = fbt('tcv',100000,0.01,'gpr',gpr,'gpz',gpz,'gpb',gpb,'gpvrr',gpvrr,'gpvrz',gpvrz,'gpvzz',gpvzz,'gpvd',gpvd,'gpve',gpve);
      testCase.verifyTrue(~isempty(LY) && logical(LY.isconverged),'FBT has not converged');
      
      % Compute vacuum field derivatives using interpolation
      Fxa = reshape(L.G.Mxa*LY.Ia,L.nzx,L.nrx);
      inM = qintc([],L.drx,L.dzx);
      [~,~,~,~,Brzi,~,~] = qintmex(L.G.rx,L.G.zx,Fxa,gpr(1),gpz(1),inM);
      
      if ~isnan(ve)
        testCase.verifyEqual(P.iohfb*gpvrz(1),Brzi,'relTol',0.1); % Cannot ask much better agreement with qint.
      end
      
      if testCase.verbosity
        %%
        clf;
        contour(L.rx,L.zx,LY.Fx);
        hold on
        contour(L.rx,L.zx,Fxa,'k--'); 
        axis equal;
        title(sprintf('Brz=%2.2f',Brzi));
        drawnow;
        %%
      end
    end
    
    function compare_fbt_vs_mds(testCase,mds_shot)
      %% Tests for FBTE
      [~,~,LY] = fbt('TCV',mds_shot,[]);
      testCase.verifyTrue(~isempty(LY) && all(logical(LY.isconverged)),'FBT has not converged');
      
      % compare with results of FORTRAN version for existing TCV shots
      mdsopen('pcs',mds_shot);
      % "RAMP" always has as many times as LY.t
      Ia  = mdsvalue('tcv_eq("i_pol"   ,"ramp")');
      FAB = mdsvalue('tcv_eq("psi_axis","ramp")');
      Ia  = Ia(1:16,:);
      FAB = reshape(FAB,1,[]);
      
      testCase.verifyEqual(LY.Ia        , Ia  ,'absTol',testCase.tolI*repmat(max(abs(Ia),[],1),size(Ia,1),1));
      testCase.verifyEqual(LY.FA - LY.FB, FAB ,'relTol',testCase.tolF);     
    end
    
    function compare_merge_nomerge(testCase,mds_shot)
      %% Tests for FBTE
      tol = 1e-6;
      [L,LX] = fbt('TCV',mds_shot,[],'mergex',false,'debug',0);
      LX.tol(:) = tol;
      LY = fbtt(L,LX);
      % reprocess with merge
      L.P.mergex = true;
      LX2 = fbtx(L,LX );
      LY2 = fbtt(L,LX2);
      
      testCase.verifyTrue(structcmp(LY,LY2,1e-12,{'chie'}),'LX with merged constraints yields different solution');
    end
  end
end
