%% Script to test TCV doublet preparation programmatically
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

clear;
saveshot = 103123;

tag = 'closer-doublets';

M=mgp(103890); % starting point
t1 = 0.01; t2=0.04; nt = 4;

M.timeeq = linspace(t1,t2,nt);

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

nth = 12;

switch tag
  case 'closer-doublets'
    at = 0.2; ab=at;
    z0 = 0.40; r0=0.85;
    
    dtu=0.3; dbl=0.3;
    dtl=0.6; dbu=0.6;
    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;
  otherwise, error('unknown case')
end

zt = z0*ones(size(M.timeeq)) + dzt;
zb = z0*ones(size(M.timeeq)) + dzb;
rt = r0*ones(1,nt)+drt;
rb = r0*ones(1,nt)+drb;

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);
  
  rlia1(:,kt) = [rat;rab];
  zlia1(:,kt) = [zat;zab];
  
  cla; plot(rlia1(:,kt),zlia1(:,kt)  ,'x'); pause(0.1);
end
M.rlia1   = rlia1;
M.zlia1   = zlia1;
M.placu1  = linspace(400,400,nt);
M.ierat   = -1;
M.iscramb = -101;
M.istop   = -320;
M.nzaxel  = nt;
M.save(saveshot);

%
[Lfbt,LXfbt] = fbt('tcv',saveshot,[],'selx','X',...
  'bfct',@bfabmex,'bfp',[1 1],...
  'agfitfct',@meqfit3,'agfitp',[true;false;true],'bp',0,...
  'fancyplot',true,'icoilon',ones(16,nt),'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

LXfbt.IpD(1,:) = LXfbt.IpD(1,:).*(1+Iptscal); % scale top current optionally

LYfbt = fbtt(Lfbt,LXfbt);
% %% 
% clf;
% meqplott(Lfbt,meqxk(LYfbt,2))
%
clf;
meqmovie(Lfbt,LYfbt);
plot(LXfbt.gpr(:,end),LXfbt.gpz(:,end),'ow'); 
drawnow;
%% FGE
[L] = fge('tcv',saveshot,[],'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 ii=1:numel(LYfbt.t)
L.P.debug=1;
LX = fgeinit(L,meqxconvert(Lfbt,meqxk(LYfbt,ii),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
  fig = gcf;
  set(gcf,'paperpositionmode','auto');
  fig_pos = fig.PaperPosition;
  fig.PaperSize = [fig_pos(3) fig_pos(4)];
  print('-depsc',sprintf('FBT_doublet_access_%s_%d',tag,ii));
end

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

figure(2); clf;
plot(LX.FB(end)/LX.FX,E(1),'x');
end
%%
dosave=1;
if dosave
 
  save doublet_control_target L LX
end