%% Droplet-double transition scan for TCV
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

clear;
shot = 103890;
tag = 'symmetric-stretchdoublet';

t1 = 0.01; t2=0.04; nt = 11;

z0 = 0.47; r0=0.87; 
k=1.2; nth = 11;

Ip0 = 200e3;
switch tag
  case 'closer-doublets'
    at = 0.2; ab=at;
    z0 = 0.40; r0=0.85;
    
    dtu=0.3; dbl=0.3;
    dtl=0.4; dbu=0.4;
    dzt = 0; dzb = 0;
    su = 0.0; sl = 0;
    
    drt = 0; drb = 0;
    Iptscal = 0;
    [ktt,kbt,kbb,kbt] = deal(k);
    ktb = ktt + linspace(0.0,0.5,nt);
    kbt = ktb;
    
    nopptile = false;
    
  case 'asymmetric-porttiles'
    at=0.24; ab=0.2;
    
    dtu=0.3; dbl=0.3;
    dtl=0.6; dbu=0.6;
    dzt = - [0 1 2 3]/100;
    dzb = 0;
    su = 0.0; sl = -0.2;
    
    drt = +0.01; drb = -0.02;
    Iptscal = linspace(0,0.5,nt);
    [ktt,kbt,kbb,kbt] = deal(k);
    ktb = ktt + linspace(0.4,0.8,nt);
    kbt = ktb + linspace(0.4,0.8,nt);
    
    nopptile = false;

  case 'symmetric-porttiles'
    at=0.21; ab=0.21;
    
    r0 = 0.85; z0 = 0.47;
    dtu=0.3; dbl=0.3;
    dtl=0.6; dbu=0.6;
    su = 0.0; sl = 0;
    
    [drt,drb,dzt,dzb] = deal(0);

    Iptscal = 0;
    [ktt,kbt,kbb,kbt] = deal(k);
    ktb = ktt + linspace(0.4,1,nt);
    kbt = ktb;
    nopptile = false;

  case 'symmetric-noporttiles'
    % symm
    at=0.23; ab=0.23;
    dtu=0.2; dbl=0.2;
    dtl=0.6; dbu=0.6;
    dzt = 0; dzb = 0;

    % sym
    su = 0.0; sl = 0;
    drt = 0; drb = 0;
    
    Iptscal = 0 ;
    [ktt,kbt,kbb,kbt] = deal(k);
    ktb = ktt + linspace(0.3,0.8,nt);
    kbt = kbt + linspace(0.3,0.8,nt);
    nopptile = true;
  case 'symmetric-stretchdoublet'
    at=0.21; ab=0.21;
    
    r0 = 0.85; z0 = 0.47;
    dtu=0.3; dbl=0.3;
    dtl=0.6; dbu=0.6;
    su = 0.0; sl = 0;
    
    [drt,drb,dzt,dzb] = deal(0);
    drt = linspace(-2 ,0,nt)/100;
    drb = linspace(-2 ,0,nt)/100;
    dzt = linspace(-10,0,nt)/100;
    dzb = linspace(-10,0,nt)/100;

    Iptscal = 0;
    [ktt,kbt,kbb,kbt] = deal(k);
    ktb = ktt + linspace(0,1,nt);
    kbt = ktb;
    nopptile = false;
  otherwise, error('unknown case')
end

zt = z0*ones(1,nt) + dzt;
zb = z0*ones(1,nt) + dzb;
rt = r0*ones(1,nt)+drt;
rb = r0*ones(1,nt)+drb;
Ip  = linspace(Ip0,Ip0,nt);


%% Prep FBT
%
[Lfbt,LXfbt] = fbt('tcv',shot,0.05,'selx','X',...
  'bfct',@bfabmex,'bfp',[1 1],...
  'agfitfct',@meqfit3,'agfitp',[true;false;true],'bp',0,...
  'fancyplot',true,'icoilon',ones(16,1),'dissi',2e-15,'debug',0,'iterq',50);

if nopptile
  Go = fbtgtcv(55331,Lfbt.P);
  Lfbt.G.rl = Go.rl; Lfbt.G.zl = Go.zl;
  Lfbt = fbtc(Lfbt.P,Lfbt.G);
end
%% Loop over equils

for kt=1:nt
  [rat,zat] = bouD(nth,rt(kt),+zt(kt),at,ktt,ktb(kt),dtu,dtl,su,sl,0,0);
  [rab,zab] = bouD(nth,rb(kt),-zb(kt),ab,kbt(kt),kbb,dbu,dbl,sl,su,0,0);

  k=~isnan(LXfbt.gpr);
  LXfbt.IpD(1:2) = Ip(kt); % assign custom/different Ip
  LXfbt.IpD(1,:) = LXfbt.IpD(1,:).*(1+Iptscal); % scale top current optionally

  ra = [rat;rab]; za = [zat;zab];
  LXfbt.gpr(k) = ra; % assign new r positions of control points
  LXfbt.gpz(k) = za; % assign new z positions of control points


  LYfbtt(kt) = fbtt(Lfbt,LXfbt);
  cla;
  meqplotfancy(Lfbt,LYfbtt(kt));
  plot(ra,za  ,'ow'); pause(0.1);

end
LYfbt = meqlpack(LYfbtt); % pack

%% FGE
[L] = fge('tcv',shot,[],'bfct',@bfabmex,'bfp',[1 1],...
  'agcon',{{'Ip','bp'},{'Ip','bp'},{'ag','ag'}},...
  'idoublet',1,'iterq',50,'selu','e','nu',50,'fancyplot',true);

if nopptile
Go = fgegtcv(55331,L.P);
L.G.rl = Go.rl; L.G.zl = Go.zl;
L = fgec(L.P,L.G);
end

for kt=1:numel(LYfbt.t)
L.P.debug=1;
LX = fgeinit(L,meqxconvert(Lfbt,meqxk(LYfbt,kt),L));
L = fgel(L,LX);

%%
figure(1); clf;
subplot(131)
meqplotfancy(L,LX); 
title(sprintf('t=%01.3f\nI_{p,t|b}=%3.0f|%3.0fkA\nb_{p,1|2}=%2.2f|%2.2f',...
 LX.t,LX.IpD(1)/1e3,LX.IpD(2)/1e3,LX.bpD(1),LX.bpD(2)));
subplot(132);
meqploteig(L,'imodes',1);
subplot(133);
meqploteig(L,'imodes',2);

doplot = true;
if doplot
  set(gcf,'paperpositionmode','auto');
  print('-depsc',sprintf('FBT_doublet_access_%s_%d',tag,kt));
end

%%
sys = fgess(L);
[V,D] = eig(sys.A);
[E,isort] = esort(diag(D));

Eigs(:,kt) = E;
LYfge(kt) = LX;
end

%%
figure(2); clf;
LY = meqlpack(LYfge);
semilogy(nanmin(sign(LY.Ip).*LY.FB)./LY.FX,Eigs(1,:),'x');
xlabel('F_{LCFS}/F_X'); ylabel('Largest growth rate [1/s]')
ha = annotation('textbox'); 
ha.Position = [0.3 0.5 0.2 0.1]; ha.String='Doublet'; ha.EdgeColor='none';
ha = annotation('textbox'); 
ha.Position = [0.7 0.5 0.2 0.1]; ha.String='Droplet'; ha.EdgeColor='none';
vline(1,'k--')
title(sprintf('growth rates: %s',tag))

if doplot
  set(gcf,'paperpositionmode','auto');
  print('-depsc',sprintf('FBT_doublet_access_%s_growthrates_flux',tag));
end

%%
figure(3); clf;
LY = meqlpack(LYfge);
semilogy(-diff(LY.zA),Eigs(1,:),'x');
xlabel('z_{A,top}-z_{A,bot}'); ylabel('Largest growth rate [1/s]')
title(sprintf('growth rates: %s',tag))

if doplot
  set(gcf,'paperpositionmode','auto');
  print('-depsc',sprintf('FBT_doublet_access_%s_growthrates_dz',tag));
end

