classdef doublet_fbt_tests < meq_test
  % Tests involving doublet cases of FBT
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    doprint = false;
    verbosity = 1;
    shot = 55331;
    t = 0.2;
    tol = 1e-10;
  end
  
  methods(Test,TestTags={'fbt'})
    
    function fbt_doublet_meqint_test(testCase)
      % Verify domain integrals
      [~,~,LY] = fbt('tcv',testCase.shot);
      
      tol = 1e-4;
      diag = 'Sum of per-domain values do not match total value for ''%s''';
      testCase.verifyEqual(sum(LY.FtD) ,LY.Ft ,'RelTol',tol,sprintf(diag,'Ft'));
      testCase.verifyEqual(sum(LY.Ft0D),LY.Ft0,'RelTol',tol,sprintf(diag,'Ft0'));
      testCase.verifyEqual(sum(LY.IpD) ,LY.Ip ,'RelTol',tol,sprintf(diag,'Ip'));
      testCase.verifyEqual(sum(LY.VpD) ,LY.Vp ,'RelTol',tol,sprintf(diag,'Vp'));
      testCase.verifyEqual(sum(LY.WkD) ,LY.Wk ,'RelTol',tol,sprintf(diag,'Wk'));
      testCase.verifyEqual(sum(LY.WpD) ,LY.Wp ,'RelTol',tol,sprintf(diag,'Wp'));
      testCase.verifyEqual(sum(LY.WtD) ,LY.Wt ,'RelTol',tol,sprintf(diag,'Wt'));
      testCase.verifyEqual(sum(LY.Wt0D),LY.Wt0,'RelTol',tol,sprintf(diag,'Wt0'));
    end
    
    function fbt_doublet_scan(testCase)
      % scan some FBT equilibria transitioning between doublet and droplet
      
      [L,LX] = fbt('tcv',testCase.shot,0.02,'izgrid',true,...
        'agfitfct',@meqfit3,...
        'bfct',@bf3pmex,'bfp',false,...
        'iterq',10,...
        'ifield',true,...
        'ivacuum',true,...
        'selx','X');
      LX.bpD(1:2) = 0.25;
      
      %% Three FBT Droplet/Doublet solutions
      fac = [0,0.12,0.03];
      titstr={'Droplet','Doublet-Droplet','Doublet'};
      for ii=1:3
        k = ~isnan(LX.gpz); % valid points
        gpz = LX.gpz(k);
        gpz = gpz-sign(gpz).*fac(ii);
        LX.gpz(k) = gpz;
        LY = fbtt(L,LX);
        
        % check constraints are satisfied
        testCase.verifyEqual(LX.bpD,LY.bpD,'AbsTol',testCase.tol)
        testCase.verifyEqual(LX.IpD,LY.IpD,'AbsTol',testCase.tol)
        testCase.verifyEqual(LX.qA,LY.qA,'AbsTol',testCase.tol)
       
        %% plot
        if testCase.verbosity
          if ii==1, clf; end
          subplot(1,3,ii)
          meqplotfancy(L,LY,'vacuum',true,'openFieldlineStyle','--','openFieldlineColor','k');
          title(titstr{ii})
        end
      end
      
      if testCase.verbosity && testCase.doprint
        set(gcf,'position',[50 50 800 400]); shg;
        print('-depsc','doublet_droplet_fbt')
      end
      
    end
    
    function fbt_asymmetric_scan(testCase)
      % scan some FBT equilibria transitioning between doublet and droplet
      
      [L,LX] = fbt('tcv',testCase.shot,0.02,...
        'izgrid',true,'ivacuum',true,...
        'agfitfct',@meqfit3,...
        'bfct',@bf3pmex,'bfp',false,...
        'issym',false,'iterq',10);
      
      LX.bpD = [0.25;0.25;0];
      %% Three droplets at different Ip and zA
      scalgrid = [-0.3 0 0.4]; rscalgrid = [-0.1,0,0.1];
      Ip0 = LX.IpD(1); gpr0 = LX.gpr;
      for jj = 1:2
        for ii=1:numel(scalgrid)
          switch jj
            case 1
              LX.IpD(1:2) = Ip0*(1+[1;-1]*scalgrid(ii));
            case 2
              LX.gpr(LX.gpz>0) = gpr0(LX.gpz>0)*(1-rscalgrid(ii));
          end
          LY = fbtt(L,LX);
          
          % check constraints are satisfied
          testCase.verifyEqual(LX.bpD,LY.bpD,'AbsTol',testCase.tol)
          testCase.verifyEqual(LX.IpD,LY.IpD,'AbsTol',testCase.tol)
          testCase.verifyEqual(LX.qA ,LY.qA ,'AbsTol',testCase.tol)
          
          if testCase.verbosity
            if ii==1 && jj==1, clf; end
            subplot(2,3,ii+(jj-1)*3)
            meqplotfancy(L,LY); hold on;
            title(sprintf('IpD=[%2.0f,%2.0f]kA\nrA=[%2.0f,%2.0f]cm',LY.IpD(1:2)/1e3,LY.rA(1:2)*100),...
              'FontSize',10)
          end
        end
      end
      
      if testCase.verbosity && testCase.doprint
        set(gcf,'position',[50 50 800 400]); shg;
        print('-depsc','asymmetric_doublet')
      end
    end
    
    function fbt_stable_doublet(testCase)
      % construct stable droplet equilibria imposing vacuum dBr/dz values
      % at axis
      shot = 103888; %#ok<*PROP>
      P = mgp(shot);
      
      vrz = 0.03; % dBrdz desired
      
      [gpr,gpz,gpb,gpvrr,gpvrz,gpvzz,gpve] = deal(nan(2,numel(P.timeeq)));
      gpvd = nan(1,numel(P.timeeq));
      gpr(1:2,:) = [1; 1]*P.rmajo1;
      gpz(1:2,:) = [1;-1]*P.zmajo1;
      gpb(1:2,:) = 0;
      gpvrr(1:2,:) = NaN;
      gpvrz(1:2,:) = vrz;
      gpvzz(1:2,:) = NaN;
      gpve(1:2,:) = 0; % NaN to disable, 0 for exact constraint, 0.001 is ok, 0.1 too high to get desired value
      gpvd(:) = 1;
      
      %%
      [L,LX] = fbt('tcv',shot,[],...
        'bfct',@bf3pmex,'bfp',false,...
        'agfitfct',@meqfit3,...
        'gpr',gpr,'gpz',gpz,'gpb',gpb,...
        'gpvrr',gpvrr,'gpvrz',gpvrz,'gpvzz',gpvzz,'gpvd',gpvd,'gpve',gpve,...
        'izgrid',true,'ivacuum',true);

      LX.bpD = repmat([0;0;0],1,numel(LX.t));
      %%
      LY = fbtt(L,LX);
      kt = 1;
      LYt = meqxk(LY,kt);
      %% Compute vacuum field derivatives using interpolation
      Fxa = reshape(L.G.Mxa*LYt.Ia,L.nzx,L.nrx);
      inM = qintc([],L.drx,L.dzx);
      [~,~,~,~,Brzi,~,~] = qintmex(L.G.rx,L.G.zx,Fxa,gpr(end-1:end,kt),gpz(end-1:end,kt),inM);
      target = L.P.iohfb*vrz*[1;1];
      testCase.verifyEqual(Brzi,target,'RelTol',0.15); % can not expect more for qint
      %%
      if testCase.verbosity
        clf;
        meqplotfancy(L,LYt); hold on;
        contour(L.rx,L.zx,reshape(L.G.Mxa*LYt.Ia,L.nzx,L.nrx),'k--');
        LXt = meqxk(LX,kt); plot(LXt.gpr,LXt.gpz,'ok')
      end
    end
  end
end
