classdef (SharedTestFixtures={meq_fixture}) ...
    meq_test < matlab.unittest.TestCase
  % Superclass for meq tests
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
 
  properties
    meqpath = fullfile(fileparts(mfilename('fullpath')),'..');
  end
  
  methods(TestMethodSetup)
    function rng_seed_setup(testCase)
      s = rng(1);
      testCase.addTeardown(@() rng(s));
    end
  end
  
  % Generall useful static methods
   methods (Static)
    function [L,LY] = getCircularEquilibrium(L,R0,Z0,FA,FB,rBt)
      % Flux map and plasma domain
      Fx = ((L.rx'-R0).^2 + (L.zx-Z0).^2)*(FB-FA) + FA; % Dummy flux map - initial guess
      [qdl,kl] = min((L.G.rl-R0).^2 + (L.G.zl-Z0).^2);
      rB = L.G.rl(kl); zB = L.G.zl(kl);
      Fx = FA + 1/qdl*(Fx-FA); % Rescale to get correct FB
      a0 = sqrt(qdl);          % Plasma minor radius
      
      Fy = Fx(2:end-1,2:end-1); % on inner grid
      FyN = (Fy-FA)/(FB-FA);
      Opy = int8(FyN>=0 & FyN<=1); % logical where there is plasma
      
      % Current density
      Imax = 1e4*sign(FA-FB);
      Iy = Imax*(1-FyN).*double(Opy>0);
      
      % Pack all quantities
      LY.Ip = sum(Iy(:));
      LY.Fx = Fx;
      LY.Iy = Iy;
      LY.FA = FA;
      LY.FB = FB;
      LY.Opy = Opy;
      LY.rA = R0;
      LY.zA = Z0;
      LY.dr2FA = 2/a0.^2*(FB-FA);
      LY.dz2FA = 2/a0.^2*(FB-FA);
      LY.drzFA = 0;
      LY.rB = rB;
      LY.zB = zB;
      LY.lX = false;
      LY.rX = zeros(0,1);
      LY.zX = zeros(0,1);
      LY.rBt = rBt;
      
      % Flux surfaces positions
      aq = L.pq./sqrt(((cos(L.oq)/a0).^2+(sin(L.oq)/a0).^2));
      rq = aq.*L.crq + R0; zq = aq.*L.czq + Z0;
      irq = 1./rq;
        
      % Toroidal field profile
      iTQ = 1/rBt*ones(L.npq+1,1);
        
      % Flux surface integrals
      FQ = FA + (FB-FA)*(L.pQ.').^2;
      Q0Q = 1/R0 *ones(L.npq+1,1);
      Q1Q = -(FB-FA)/(2*pi*pi*R0*a0*a0)*ones(L.npq+1,1);
      Q2Q = 1./(R0^2*sqrt(1-(a0/R0)^2*FQ));
      Q3Q = Q2Q*4.*FQ/a0^2;
      Q4Q = 4*FQ/a0^2/R0;
      iqQ = 2*pi*Q1Q./Q2Q.*iTQ;
      
      LY.aq = aq;
      LY.rq = rq;
      LY.zq = zq;
      LY.irq = irq;
      LY.Q0Q = Q0Q;
      LY.Q1Q = Q1Q;
      LY.Q2Q = Q2Q;
      LY.Q3Q = Q3Q;
      LY.Q4Q = Q4Q;
      LY.iqQ = iqQ;
      LY.iTQ = iTQ;
      
      % Gaps
      co = cos(L.G.oW)/a0;
      so = sin(L.G.oW)/a0;
      dr = (L.G.rW - R0)/a0;
      dz = (L.G.zW - Z0)/a0;
      
      % Looking for solution of a*x^2 + 2*b*x + c = 0
      a = co.^2  + so.^2;
      b = dz.*so - dr.*co;
      c = dr.^2  + dz.^2 - FB;
      aW1 = (-b+sqrt(b.^2-a.*c))./a;aW1(aW1 < 0) = Inf;
      aW2 = (-b-sqrt(b.^2-a.*c))./a;aW2(aW2 < 0) = Inf;
      aW = min(aW1,aW2);
      
      LY.aW = aW;
    end
    
    function [S] = generate_flux_map(type,sIp)
      % defaults
      ro = 1; zo = 0.4; wo = 0.2; % gaussian flux distribution parameters
      Fo = 0.01; % up/down flux level
      switch type
        case 'SND'
          Fo = Fo*[1,1]; ro = [ro,0.6]; zo = [0,-0.8];
          wo = wo*[1,0.4];
          ndom = 1; nA = 1; nX = 1; lX = true;
        case 'DND'
          Fo = Fo*[1,1,1]; ro = [ro,0.6,0.6]; zo = [0,-0.8,0.8];
          wo = wo*[1,0.4,0.4];
          ndom = 1; nA = 1; nX = 2; lX = true;
        case 'Lim'
          zo = 0;
          ndom = 1; nA = 1; nX = 0;  lX = false;
        case 'Lim-X'           % limited with a non-limiting X-point
          Fo = Fo*[1,1]; ro = [ro,0.6]; zo = [zo-0.2,-0.8];
          ndom = 1; nA = 1; nX = 1; lX = true;
        case 'Double-Snowflake-Minus'
          ndom = 1; nA = 1; nX = 4; lX = true;
          Fo = Fo*[1.6,1,1,1,1];
          ro = [ro,0.6,0.6,1.4,1.4];
          zo = [0,-0.8,0.8,-0.8,0.8];
        case 'Double-Snowflake-Plus'
          ndom = 1; nA = 1; nX = 4; lX = true;
          Fo = Fo*[1.1,1,1,1,1];
          ro = [ro,0.6,0.6,1.4,1.4];
          zo = [0,-0.8,0.8,-0.8,0.8];
        case 'Doublet'
          ro = ro*[1,1]; zo = zo*[-1,1]; wo = wo*[1,1];  Fo = Fo*[1,1];
          ndom = 3; nA = 2; nX = 1;  lX = [true;true;false];
        case 'Droplet'
          zo = (zo+0.2)*[-1,1];
          ndom = 2; nA = 2; nX = 1; lX = [false;false];
        case 'Doublet-div-nomantle'
          ro = ro*[1,1,1]; zo = [-zo,zo+0.1,-1.1]; wo = wo*[1,1,.3];  Fo = Fo*[1,1.1,2];
          ndom = 2; nA = 2; nX = 2; lX = [true;true];
        case 'Doublet-div'
          ro = ro*[1,1,1]; zo = [-zo,zo,-1]; wo = wo*[1,1,.2];  Fo = Fo*[1,1,1];
          ndom = 3; nA = 2; nX = 2; lX = [true;true;true];
        case 'Triplet'
          ro = ro*[1,1,1]; zo = 1.2*[-zo,0,zo]; wo = 0.2*wo*[1,1,1];  Fo = Fo*[2,1,1];
          ndom = 5; nA = 3; nX = 2; lX = [true;true;true;true;false];
        case 'Triplet-madness'
          ro = ro*[1,1,1,1]; zo = [1.2*[-zo-0.1,0,zo],-1]; wo = 0.2*wo*[1,1,1,1];  Fo = Fo*[1,1,1.3,1];
          ndom = 4; nA = 3; nX = 3; lX = [true;true;true;true];
        case 'Boundary-X'
          ro = [ro,0.55*2-ro+0.01]; zo = zo*[-1,1]; wo = wo*[1,1]; Fo = Fo*[1,1];
          ndom = 1; nA = 1; nX = 1; lX = true;
        otherwise
          error('unknown type %s',type)
      end
      
      S.Fxh = @(rrx,zzx) reshape(sum(sIp*(Fo.*exp( -((rrx(:)-ro).^2 + (zzx(:)-zo).^2)./wo)),2),size(rrx));
      S.drFxh = @(rrx,zzx) reshape(sum(sIp*(Fo.*(-2*(rrx(:)-ro)./wo).*exp( -((rrx(:)-ro).^2 + (zzx(:)-zo).^2)./wo)),2),size(rrx));
      S.dzFxh = @(rrx,zzx) reshape(sum(sIp*(Fo.*(-2*(zzx(:)-zo)./wo).*exp( -((rrx(:)-ro).^2 + (zzx(:)-zo).^2)./wo)),2),size(rrx));
      S.dr2Fxh = @(rrx,zzx) reshape(sum(sIp*(Fo.*(4*(rrx(:)-ro).^2./wo.^2 -2./wo).*exp( -((rrx(:)-ro).^2 + (zzx(:)-zo).^2)./wo)),2),size(rrx));
      S.dz2Fxh = @(rrx,zzx) reshape(sum(sIp*(Fo.*(4*(zzx(:)-zo).^2./wo.^2 -2./wo).*exp( -((rrx(:)-ro).^2 + (zzx(:)-zo).^2)./wo)),2),size(rrx));
      S.drzFxh = @(rrx,zzx) reshape(sum(sIp*(Fo.*(4*(rrx(:)-ro).*(zzx(:)-zo)./wo.^2).*exp( -((rrx(:)-ro).^2 + (zzx(:)-zo).^2)./wo)),2),size(rrx));
      S.ndom = ndom;
      S.nA = nA;
      S.nX = nX;
      S.lX = lX;
    end

    function [ok,msg] = check_tok(tok,shot)
      ok = true; msg = '';
      if strcmpi(tok, 'CREATE')
        % check existence of CREATE file and filter otherwise
        createfile = shot{1};
        ok = exist(createfile,'file');
        msg = sprintf('Can not find CREATE file %s. Filtering test for now...',createfile);
      end

      if strcmpi(tok,'TCV')
        hostname = getenv('HOSTNAME');
        ok = ~isempty(regexp(hostname,'(lac\d*\.epfl\.ch|spc-ci)','once'));
        msg = 'Ignoring TCV case since hostname is not lac*.epfl.ch or spc-ci';
      end
    end
  end
  
end
