classdef test_mga_voltage_strays < matlab.unittest.TestCase
  
  properties
    refshot = 66300;
    testshot;
    verbose = false;
  end
  
  properties(TestParameter)
    nofeedback = {true,false};
    inores = {true,false};
    iscramb = {-2,-1};
  end
  
  methods(TestClassSetup)
    function getmytestshot(testCase)
      testCase.testshot = mgu.mytestshot();
    end

    function working_folder_fixture(testCase)
      testCase.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture);
    end
  end
  
  methods(Test, TestTags={'single_instance'})
    function test_preparation(testCase,nofeedback,inores,iscramb)
      refshot  = testCase.refshot; %#ok<*PROPLC>
      testshot = testCase.testshot;
      loadshot = testshot;
      verbose  = testCase.verbose;
      MP = mgp(refshot);
      
      MP.iscramb = iscramb + sign(iscramb)*(inores*100); % encoding
      MP.inova   =  3;
      MP.mvloop  =  7;
      MP.nfast   =  0;
      
      if nofeedback
        MP.emgain = 0;
      end
      
      fname = fullfile(pwd,'voltage_stray_traces.mat');
      MP.debug_p3 = sprintf('''fname'',''%s''',fname);

      MP.save(testshot)
      
      %% FBT
      [~,~,~] = fbte(testshot);

      %%
      [Ia,Va,tout] = voltage_stray_traces('dosave',true,'tC',{[1 2 3],[]},'dt',1e-3,'fname',fname);
      
      %% modify
      mga(testshot);
      fbte18(testshot,108,'rtdtacq');

      %%
      mdsopen('pcs',loadshot);
 
      ipol = [1:16,22]; % IE, IF, IOH1
      ipolin =  77:93 ;
      AA = mdsvalue('\pcs::phys_mat_a');
      assert(norm(AA(ipol,ipolin)-eye(numel(ipol)))==0,...
        'Coil control directions should be identity for IE, IF, IOH')
      idoh = 17;
      assert(sum(logical(AA(idoh,:)))==1,'Error for dOH observer, more than one element is non-zero');
      
      G1 = mdsvalue('\pcs::phys_mat_g1');
      G2 = mdsvalue('\pcs::phys_mat_g2');
      G3 = mdsvalue('\pcs::phys_mat_g3');
      
      ioh = 17:18;
      % assert(norm(G1(ioh,:))==0,'G1 rows for IOH must be empty');
      assert(norm(G2(1:16,1:16) - diag(diag(G2(1:16,1:16))))==0,...
        'G2 matrix should be diagonal for E,F');
      assert(norm(G2(ioh,[1:16,18:21,23:24]))==0,...
        'G2 matrix should only have non-zero elements on columns 17 and 22 for OH');
      assert(~any(G3(:)),'G3 must be empty');
      
      %% check M matrix resistive part
      MM = mdsvalue('\pcs::phys_mat_m');
      jIPOL = 23:40;
      kIPOL = 1:18;
      MMIPOL = MM(kIPOL,jIPOL);
      if inores
        testCase.verifyTrue(~any(any(MMIPOL)),'Res part of M matrix must be absent for this case');
      else
        testCase.verifyTrue(all(diag(MMIPOL)),'Res part of M matrix not fully populated')
        testCase.verifyTrue(~any(any(diag(diag(MMIPOL))-MMIPOL)),'Res part of M matrix not diagonal');
      end
      
      % check gas is excluded
      inG = [1:18,20]; % exclude gas
      if nofeedback
        assert(all(all(MM(inG,inG)))==0);
      end
      %%
        
      rr = [1:8,9:11,20,19,14:16]; % refs ordering
      for ii = 1:8
        V_E{ii} = tdi(sprintf('\\DRAW_FEEDFOR_E_U:ALIM_%03d', ii))'; %#ok<*AGROW,*SAGROW>
        V_F{ii} = tdi(sprintf('\\DRAW_FEEDFOR_F_U:ALIM_%03d', ii))';
        I_E{ii} = tdi(sprintf('\\DRAW_REFS:REF_%03d',rr(ii)   ));
        I_F{ii} = tdi(sprintf('\\DRAW_REFS:REF_%03d',rr(ii+8) ));
      end
      
      for ii=1:2
        V_OH{ii} = tdi(sprintf('\\DRAW_FEEDFOR_OH_U:ALIM_%03d',ii))';
      end
      I_OH001 = tdi('\DRAW_REFS:REF_022');
      I_DOH   = tdi('\DRAW_REFS:REF_017');
      

      %% Plot and check Vs
      tol = 1e-3;
      scale = 1e4;
      if verbose, clf; end
      for ii=1:2
        [~,ia] = intersect(round(tout*scale),round(V_OH{ii}.dim{1}*scale));
        testCase.verifyEqual(V_OH{ii}.data,Va(16+ii,ia)','AbsTol',tol);
        if verbose
          cla;
          plot(V_OH{ii},'o-'); hold on; plot(tout,Va(16+ii,:),'k--');
          title(sprintf('OH%03d',ii))
          drawnow
        end
      end
      
      for ii=1:8
        [~,ia] = intersect(round(tout*scale),round(V_E{ii}.dim{1}*scale));
        testCase.verifyEqual(V_E{ii}.data,Va(0 +ii,ia)','AbsTol',tol);
        if verbose
          cla;
          plot(V_E{ii},'o-'); hold on; plot(tout,Va(0+ii,:),'k.--');
          title(sprintf('E%03d',ii))
          drawnow;
        end
      end
      
      for ii=1:8
        [~,ia] = intersect(round(tout*scale),round(V_F{ii}.dim{1}*scale));
        testCase.verifyEqual(V_F{ii}.data,Va(8+ii,ia)','AbsTol',tol);
        if verbose
          cla;
          plot(V_F{ii},'o-'); hold on; plot(tout,Va(8+ii,:),'k.--');
          title(sprintf('F%03d',ii))
          drawnow
        end
      end
      
      %% E
      for ii=1:8
        [~,ia] = intersect(round(tout*scale),round(I_E{ii}.dim{1}*scale));
        testCase.verifyEqual(I_E{ii}.data,Ia(0+ii,ia)','AbsTol',tol);
        if verbose
          cla;
          plot(I_E{ii},'.-'); hold on; plot(tout,Ia(0+ii,:),'k.--');
          title(sprintf('E%03d',ii))
          drawnow
        end
      end
      
      % F
      for ii=1:8
        [~,ia] = intersect(round(tout*scale),round(I_F{ii}.dim{1}*scale));
        testCase.verifyEqual(I_F{ii}.data,Ia(8+ii,ia)','AbsTol',tol);
        if verbose
          cla;
          plot(I_F{ii},'.-'); hold on; plot(tout,Ia(8+ii,:),'k.--');
          title(sprintf('F%03d',ii))
          drawnow
        end
      end
      
      % oh
      [~,ia] = intersect(round(tout*scale),round(I_OH001.dim{1}*scale));
      testCase.verifyEqual(I_OH001.data,Ia(17,ia)','AbsTol',tol);
      [~,ia] = intersect(round(tout*scale),round(I_DOH.dim{1}*scale));
      assert(norm(I_DOH.data  + diff(Ia(17:18,ia))')<tol);
      
      if verbose
        cla;
        plot(I_OH001,'o-'); hold on; plot(tout,Ia(16+1,:),'k.--');
        drawnow;
        cla
        plot(I_DOH,'o-'); hold on; plot(tout,-diff(Ia(17:18,:)),'k.--');
      end
    end
  end
end
