classdef fields_test<meq_test
  % Test involving magnetic fields and their calculation
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbosity = 0;
  end
  
  properties(TestParameter)
    bp = {0,2}; % high and low bp
    doublet = {false,true};
  end
  
  methods(Test,TestTags = {'Integration'})
    function field_calculation_on_custom_grid(testCase)
      % Call calculation of Br,Bz,Bt fields on custom 'z' grid
      [L,~,LY] = fbt('ana',1,[],...
        'nr',64,'nz',32,...
        'izgrid',true,'rz',linspace(0.01,2,61),'zz',linspace(-1.5,1.5,61),...  % specify extra grid
        'ifield',true,... % compute fields
        'ivacuum',true); % also compute vacuum fields
      
      testCase.assertEqual(LY.Iu,zeros(L.G.nu,1),'Vessel currents must be zero for subsequent tests to work')
      
      % Verify current matches ampere's law estimate: integrate around square outer contour
      % First order method: very coarse estimate
      % x grid
      I_expected = LY.Ip;
      I_ampere_x = ampere_integrate(LY.Brx,LY.Bzx,L.drx,L.dzx);
      testCase.verifyEqual(I_ampere_x,I_expected,'RelTol',0.05,'Current from poloidal field integral does not match equilibrium current')
      
      % Same for z grid: encompasses all conductor currents
      I_expected = LY.Ip+sum(LY.Ia);
      I_ampere_z = ampere_integrate(LY.Brz,LY.Bzz,L.drz,L.dzz);
      testCase.verifyEqual(I_ampere_z,I_expected,'RelTol',0.05,'Current from poloidal field integral does not match equilibrium current')
      
      % Also use direct method using meqg Green's functions for vacuum part
      L.G = meqg(L.G,L.P,'Brza','Brzu','Bzza','Bzzu');
      Brz = reshape(L.G.Brza*LY.Ia + L.G.Brzu*LY.Iu,L.nzz,L.nrz);
      Bzz = reshape(L.G.Bzza*LY.Ia + L.G.Bzzu*LY.Iu,L.nzz,L.nrz);
      
      I_ampere_z = ampere_integrate(Brz,Bzz,L.drz,L.dzz);
      I_expected = sum(LY.Ia); % only PF currents - no plasma
      testCase.verifyEqual(I_ampere_z,I_expected,'RelTol',0.05,'Current from poloidal field integral does not match equilibrium current')
      
      % vacuum calculation should give Brz,Bzz directly, allowing tolerance
      % for different calculation method.
      
      I_expected  = 0; % expect no current if integrating vacuum field internally to the PF coils
      I_ampere_x0 = ampere_integrate(LY.Br0x,LY.Bz0x,L.drx,L.dzx); %#ok<NASGU>
      testCase.verifyEqual(I_expected,I_ampere_x0,'AbsTol',1e3);
      
      I_expected  = sum(LY.Ia);
      I_ampere_z0 = ampere_integrate(LY.Br0z,LY.Bz0z,L.drz,L.dzz);
      testCase.verifyEqual(I_ampere_z0,I_expected,'RelTol',1e-3,'Current from vacuum field integral does not match sum of coil currents'); % these should be really equal
      
      if testCase.verbosity
        clf; 
        subplot(121)
        meqplotfancy(L,LY); hold on; 
        meqplotfield(L,LY,'izgrid',true);
        subplot(122)
        meqplott(L,LY); hold on; 
        meqplotfield(L,LY);
      end
    end
    
    function toroidal_field_should_match_TQ_interpolation_result(testCase,bp,doublet)
      
      if doublet, shot = 82; bpD = [bp/4;bp/4;0];
      else,       shot = 1;  bpD = bp;
      end
      
      % Call calculation of Br,Bz,Bt fields on custom 'z' grid
      [L,LX] = fbt('ana',shot,[],'ifield',true,...  % compute fields
        'bfct',@bf3pmex,'bfp',[],... % Direct evaluation of Btx from basis function only for this bf
        'itert',1); % this forces small diamagnetism approximation for TQ
      LX.bpD = bpD;
      LY = fbtt(L,LX);
      
      % Compute toroidal field using TQ interpolation
      Btx  = LY.rBt*repmat(1./L.rx',L.nzx,1);
      Bty = Btx(L.lxy);  
      [F0,F1] = meqpdomlim(LY.FA,LY.FB,L.nD); % order for multi-domain
      nB = numel(F1);
      for iB = 1:nB
        Opyi = (LY.Opy==iB); % mask for this domain
        FyN = (LY.Fx(L.lxy)-F0(iB))/(F1(iB)-F0(iB)); % rhopol.^2
        Bty(Opyi) = interp1(L.pQ.^2,LY.TQ(:,iB)',FyN(Opyi),'linear')./L.rry(Opyi);
      end
      Btx(L.lxy) = Bty;

      % bfct(8..) method used in meqpost() 
      % for evaluating Btx uses direct TQ evaluation
      % while we use linear interpolation. This causes a small error.
      testCase.verifyEqual(Btx,LY.Btx,...
        'RelTol',1e-4,'Plasma toroidal field calculation is incorrect')
      
      if testCase.verbosity
        clf;
        imagesc(LY.Btx-Btx0)
      end
    end
  end
end

function I_ampere = ampere_integrate(Br,Bz,dr,dz)
% auxiliary function for integral along computational domain boundary
I_ampere = 1/(4e-7*pi)*(sum(Bz(1:end-1,1)*dz) - sum(Br(1,2:end)*dr) - sum(Bz(2:end,end)*dz) + sum(Br(end,1:end-1)*dr));
end

  
