%FBTXANA  ANA FBT equilibrium specification
% LX = FBXPANA(t,L) returns a structure LX with cost function/
% constraint data for the ANA tokamak.
% Several shots are pre-defined:
%
%  1: limiter
%  2: single null
%  3: single null with secondary x point
%  4: double null
%  5: limited kappa<1
%  6: negative triangularity limited
% 11: Elongated
% 21: Circular with X point outside limiter
% 81: Droplets
% 82: Doublet
% 83: Diverted droplet with second droplet in private flux region
% 84: Diverted droplet with second droplet in public flux region
% 85: Two diverted droplets
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function LX = fbtxana(t,L)

P  = L.P;
LX = struct();
% flux points
npts = 8;
thp  = linspace(0,2*pi,npts+1)'; thp=thp(1:end-1);
a0 = 0.40;
z0 = 0;

% defaults
if isempty(t)
  LX.t = 0;
elseif numel(t)>1
  LX.t = t(1); % take first (overwrite at the end)
else
  LX.t = t;
end

LX.qA       = 1;
LX.rBt      = 1*P.r0;
LX.bpD      = 1;
LX.IpD      = 200e3 ;
LX.isaddl   = 1;
LX.niter    = 100; 
LX.idoublet = false;
LX.shot     = P.shot;

% shot dependent overrides
switch floor(LX.shot/10)
  case 8
    % doublets
    LX.idoublet = true;
    LX.IpD = [100;100;0]*1e3;
    LX.bpD = [0.25;0.25;0];
    LX.qA  = [1;1];
end


switch LX.shot
  case 0 % no plasma
    LX = fbtgp(LX,[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]);
  case 1 % limited, circular
    r0 = P.r0 + 0.05;
    ra = r0 + a0*cos(thp);
    za = z0 + a0*sin(thp);
    %     fbtgp(X,r ,z ,b,fa,fb,fe,br,bz,ba,be,cr,cz,ca,ce,vrr,vrz,vzz,ve)
    LX = fbtgp(LX,ra,za,1 ,0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
    % Exact boundary point
    LX.gpfe(2) = 0;
  case 2 % diverted
    thp = linspace(0,2*pi,10)'; thp=thp(1:end-1);
    r0 = P.r0 + 0.02;
    a0 = 0.4;
    ra = r0 + a0*cos(thp); ra(1)=0.95*ra(1);
    za = z0 + a0*sin(thp);
    
    %   fbtgp(LX,r ,z ,b,fa,fb,fe,br,bz,ba,be,cr,cz,ca,ce,vrr,vrz,vzz,ve)
    LX = fbtgp(LX,ra,za,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
    % Exact boundary point
    LX.gpfe(1) = 0;
    % Specify an X-point
    LX.gpfe(7) = 0; % Lower one is an exact boundary point
    LX.gpbr(7) = 0; % Br == 0
    LX.gpbz(7) = 0; % Bz == 0
    LX.gpbe(7) = 0; % B field constraints are exact
    LX.gpr(7)  =  0.85;
    LX.gpz(7)  = -0.4;
    
  case 3 % diverted with secondary x-point
    thp = linspace(0,2*pi,4)'; thp=thp(1:end-1);
    r0 = P.r0;
    a0 = 0.45;
    ra = r0 + a0*cos(thp); ra(1)=0.95*ra(1);
    za = z0 + a0*sin(thp);
    
    %    fbtgp(P ,r ,z ,b,fa,fb,fe,br,bz,ba,be,cr,cz,ca,ce,vrr,vrz,vzz,ve)
    LX = fbtgp(LX,ra,za,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
    % Specify an X-point
    LX.gpfe( 2   )   = 0; % Lower one is an exact boundary point
    LX.gpbr([1,2,3]) = 0; % Br == 0
    LX.gpbz([1,2,3]) = 0; % Bz == 0
    LX.gpbe([2,3]) = 0; % B field constraints are exact
  case {4,94} % double null diverted
    r0 = P.r0 + 0.02;
    a0 = 0.4;
    % X points
    rX = r0 - [1; 1]*a0*0.5;
    zX = z0 + [1;-1]*a0*1.0;
    % boundary and strike points
    rr = [1; 1]*(r0+[0.25 0.0 -0.23 -0.19 -0.28]);
    zz = [1;-1]*[0.1 0.32 0.2 0.45 0.40];
    r = [rX;rr(:)]; z = [zX;zz(:)];
    
    LX = fbtgp(LX,r,z,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
    % Specify an X-point
    LX.gpfe(1:2)   = 0; % Make X-points exact constraints
    LX.gpbr([1,2]) = 0; % Br == 0
    LX.gpbz([1,2]) = 0; % Bz == 0
    LX.gpbe([1,2]) = 0; % B field constraints are exact
    % Remove strike points from initial guess
    LX.gpb(end-3:end) = 0;
  case 5
    % squashed (vertically stable)
    r0 = P.r0 + 0.05;
    ra = r0 + a0*cos(thp);
    za = z0 + 0.9*a0*sin(thp);
    %   fbtgp(P,r ,z ,b,fa,fb,fe,br,bz,ba,be,cr,cz,ca,ce,vrr,vrz,vzz,ve)
    LX = fbtgp(LX,ra,za,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
    % Exact boundary point
    LX.gpfe(2) = 0;
  case 6 % NT limited with secondary x points
    r0 = P.r0 + 0.02;
    a0 = 0.4;
    % X points
    rX = r0 + [1; 1]*a0*0.5;
    zX = z0 + [1;-1]*a0*1.0;
    % boundary points
    rr = [1; 1]*(r0-0.15+[0.25 0.0 -0.23]);
    zz = [1;-1]*[0.2 0.32 0.2];
    r = [rX;rr(:)]; z = [zX;zz(:)];
    
    LX = fbtgp(LX,r,z,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
    % Specify an X-point
    LX.gpbr([1,2]) = 0; % Br == 0
    LX.gpbz([1,2]) = 0; % Bz == 0
    LX.gpbe([1,2]) = 0; % B field constraints are exact
  case 81 % two droplets
    z0 = 0.55; k = 1;
    [ra1,za1] = bouD(31,P.r0, z0,0.25,k,k,0,0,0,0,0,0);
    [ra2,za2] = bouD(31,P.r0,-z0,0.25,k,k,0,0,0,0,0,0);
    ra = [ra1;ra2]; za = [za1;za2];
    LX = fbtgp(LX,ra,za,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
  case 82 % doublet!
    [ra,za] = bouD(31,P.r0,0,0.7,1.2,1.2,0,0,0,0,0.75,0.75);
    LX = fbtgp(LX,ra,za,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
  case 83 % diverted droplet with second droplet in private flux region
    z0 = 0.55; k = 1;
    [ra1,za1] = bouD(31,P.r0, z0+0.1,0.25,k,k,0,0,0,0,0,0);
    [ra2,za2] = bouD(31,P.r0,-z0+0.1,0.25,k,k,0,0,0,0,0,0);
    ra = [ra1;ra2]; za = [za1;za2];
    LX = fbtgp(LX,ra,za,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
  case 84 % diverted droplet with second droplet in public flux region
    z0 = 0.53; k = 1; npt = 31;
    [ra1,za1] = bouD(npt,P.r0, z0,0.35,k,k,0,0,0,0,0,0);
    [ra2,za2] = bouD(npt,P.r0,-z0,0.35,k,k,0,0,0,0,0,0);
    ra = [ra1;ra2]; za = [za1;za2];
    LX = fbtgp(LX,ra,za,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
    ix = (npt+1)/4+1;
    % Specify an X-point
    LX.gpfe(ix) = 0; % Lower one is an exact boundary point
    LX.gpbr(ix) = 0; % Br == 0
    LX.gpbz(ix) = 0; % Bz == 0
    LX.gpbe(ix) = 0; % B field constraints are exact
  case 85 % two diverted droplets
    z0 = 0.53; k = 1; npt = 31;
    [ra1,za1] = bouD(npt,P.r0, z0,0.35,k,k,0,0,0,0,0,0);
    [ra2,za2] = bouD(npt,P.r0,-z0,0.35,k,k,0,0,0,0,0,0);
    ra = [ra1;ra2]; za = [za1;za2];
    LX = fbtgp(LX,ra,za,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
    ix = (npt+1)/4+1;
    % Specify an X-point
    LX.gpfe(ix) = 0; % Lower one is an exact boundary point
    LX.gpbr(ix) = 0; % Br == 0
    LX.gpbz(ix) = 0; % Bz == 0
    LX.gpbe(ix) = 0; % B field constraints are exact
    
    % Specify an X-point
    ix = (7*npt+3)/4;
    LX.gpfe(ix) = 0; % Lower one is an exact boundary point
    LX.gpbr(ix) = 0; % Br == 0
    LX.gpbz(ix) = 0; % Bz == 0
    LX.gpbe(ix) = 0; % B field constraints are exact
  case 11 % high elongation with triangularity
    z0 = 0; k = 2.5; d = 0.5; a = 0.32;
    [ra,za] = bouD(31,P.r0+0.1, z0,a,k,k,d,d,0,0,0,0);
    LX = fbtgp(LX,ra,za,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
    LX.Ip = 500e3;
  case 21 % circular with X-point outside limiter
    r0 = P.r0 + 0.05;
    ra = r0 + a0*cos(thp);
    za = z0 + a0*sin(thp);
    %    fbtgp(LX,r ,z ,b,fa,fb,fe,br,bz,ba,be,cr,cz,ca,ce,vrr,vrz,vzz,ve)
    LX = fbtgp(LX,ra,za,1, 0, 1, 1,[],[],[],[],[],[],[],[],[] ,[] ,[] ,[]);
    rx = r0 + .9*a0;
    zx = z0 + .9*a0;
    LX = fbtgp(LX,rx,zx,0, 0, 0,NaN,0,0,NaN,0,[],[],[],[],[] ,[] ,[] ,[]);
  otherwise
    error('Shot %d not available for anamak',LX.shot);
end

na = P.nc; % number of coils

% Cost function weights
LX.gpfd = 1e-3/size(LX.gpr,1);
LX.gpbd = 1e-3/size(LX.gpr,1);
LX.gpcd = 1e-3/size(LX.gpr,1);
LX.gpid = LX.gpfd./(na*2*pi*sqrt(2*P.dissi));
LX.gpdd = LX.gpfd./(   2*pi*sqrt(  P.dipol));

% repeat for multiple time slices
if numel(t)>1
  LX = repmat(LX,numel(t),1);
  LX = meqlpack(LX);
  LX.t = t;
end

% PFC currents
nt = numel(LX.t);
LX.gpia = 0*ones(na,nt);     % Minimize total current in the coils
LX.gpie = 1*ones(na,nt);     % id.
LX.gpdw = ones(na,nt);       % Simple current dipoles
LX.gpde = NaN*ones(na-1,nt); % no weight for coil dipoles

end
