%LIUGCREATE CREATE(ITER) LIUQE geometry
% LIUGCREATE(EQUIL,P) returns a structure with geometry for ITER. See also LIUG.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function G = liugcreate(equil,P)

%% Load equilibrium file
file = createf(equil);
if isempty(file)
  warning('Cannot find file for equilibrium #%05d',equil)
  return
end

D = load(file{1});

%% Poloidal field coils
rc = D.coils_data.R(:);
zc = D.coils_data.Z(:);
wc = D.coils_data.DR(:);
hc = D.coils_data.DZ(:);
nc = D.coils_data.turns(:);
G.dimc = extractAfter(D.coils_data.names,'coil ');
G.wc = wc; G.hc = hc;
G.rc = rc; G.zc = zc;

% Discretise
d = .05; % typical filament size
nr = round(wc/d);
nz = round(hc/d);
nt = nr .* nz;
dr = wc ./ nr;
dz = hc ./ nz;
r1 = rc-0.5*wc+0.5*dr;
z1 = zc-0.5*hc+0.5*dz;
G.Twc = zeros(sum(nt),length(nc));
[G.rw,G.zw,G.ww,G.hw] = deal(zeros(size(G.Twc,1),1));
l = 0;
for kc = 1:length(nc)
  k = (1:nt(kc))+l;
  l = k(end);
  [rr,zz] = meshgrid(r1(kc)+(0:nr(kc)-1)*dr(kc),...
    z1(kc)+(0:nz(kc)-1)*dz(kc));
  G.rw(k) = rr(:);
  G.zw(k) = zz(:);
  G.ww(k) = dr(kc);
  G.hw(k) = dz(kc);
  G.Twc(k,kc) = nc(kc)/nt(kc);
end
G.dw = nan(size(G.rw));

Tca = eye(numel(nc)); dima = G.dimc; % init
if not(isfield(D.coils_data,'connection'))
  % connect CS1U+CS1L and VS3U-VS3L (default for ITER)
  [Tca,dima] = connect_coils(dima,Tca,'CS1U','CS1L','CS1',+1);
  [Tca,dima] = connect_coils(dima,Tca,'VS3U','VS3L','VS3',-1);
else
  dima = G.dimc;
  for i = 1 : size(D.coils_data.connection,1)
    % different connection specified
    coil1  = D.coils_data.connection{i,1};
    coil2  = D.coils_data.connection{i,2};
    circ   = D.coils_data.connection{i,3};
    series = D.coils_data.connection{i,4};
    [Tca,dima] = connect_coils(dima,Tca,coil1,coil2,circ,series);
  end
end
G.dima = dima;
G.Twa = G.Twc*Tca;
G.Tca = Tca;

% P*->PF* for naming consistency
for ic=1:numel(G.dimc)
  if ~isempty(regexp(G.dimc{ic},'P[\d]', 'once'))
    G.dimc{ic} = ['PF',G.dimc{ic}(2:end)];
  end
end
for ia=1:numel(G.dima)
  if ~isempty(regexp(G.dima{ia},'P[\d]', 'once'))
    G.dima{ia} = ['PF',G.dima{ia}(2:end)];
  end
end
na = size(G.Tca,2);

% Resistance
G.Ra = zeros(na,1); % Supercondicting coils
G.Ra(contains(G.dima,'VS3')) = 12e-3; % Non-superconducting VS3 internal coil, resistance value from GRT47

% limits
G.Ialim = zeros(na,1);
G.Ialim(startsWith(G.dima,'PF')) = 48e3; % from B.Lim  et al, IEEE trans. Supra 2011 (doi:10.1109/TASC.2010.2092732)
G.Ialim(startsWith(G.dima,'VS')) = 60e3; % from I.Song et al, IEEE trans. Supra 2013 10.1109/TASC.2013.2292671
G.Ialim(startsWith(G.dima,'CS')) = 45e3; % from USDA Central Solenoid Fact Sheet

%% Flux loops
ksel = find(~cellfun(@isempty,regexp([D.saddle.name],P.selfm)));
G.rf   = [D.saddle(ksel).r   ]';
G.zf   = [D.saddle(ksel).z   ]';
G.dimf = [D.saddle(ksel).name]';
tf = reshape(([D.saddle(ksel).TorAngEnd] - [D.saddle(ksel).TorAngStart])/(2*pi),[],1);
G.Tff = zeros(length(ksel),length(G.rf));
l = 0;
for kf = 1:size(G.Tff,1)
  k = (1:length(D.saddle(ksel(kf)).r))+l; l = k(end);
  G.Tff(kf,k) = tf(k);
end

%% Magnetic probes
ksel = find(~cellfun(@isempty,regexp(D.field_sensors.names,P.selfm)));
G.rm   = D.field_sensors.r(ksel);
G.zm   = D.field_sensors.z(ksel);
G.am   = D.field_sensors.theta(ksel);
G.dimm = D.field_sensors.names(ksel);

%% Limiter
G.rl = D.limiter.r(:);
G.zl = D.limiter.z(:);

%% Vessel
if any(P.vesselmode==[1,2])
  assert(isfield(D,'conductors'),'need ''conductors'' field to be specified for vesselmode>0');
  C = D.conductors;
end

switch P.vesselmode
  case 0
    % use only ITER vessel information as filaments (filaments)
    [G.rv,G.zv,G.wv,G.hv,G.Rv,G.nv,G.dimv] = ITERvessel2filaments(ITER_vessel_data());
  case 1
    % map CREATE segments to filaments based on ITER vessel
    [G.rv,G.zv,G.wv,G.hv,G.Rv,G.nv,G.dimv] = ITERvessel2filaments(ITER_vessel_data());
    
    % overwrite resistance data to match CREATE GRT47 reverse-engineered values
    rho = zeros(G.nv,1);
    for iv=1:G.nv
      if any(startsWith(G.dimv{iv},{'in','out'}))
        rho(iv) = 0.8e-6;
      elseif startsWith(G.dimv{iv},{'div'})
        rho(iv) = 0.9e-6;
      elseif startsWith(G.dimv{iv},{'tria'})
        rho(iv) = 0.3224e-6;
      else
        error('unrecognized dimv %s',G.dimv{iv})
      end
    end
    G.Rv = 2*pi*rho.*G.rv./(G.wv.*G.hv);
    
    % assign filaments to a segment based on whether they are inside one
    G.ns = numel(C.PassiveConductor); % number of CREATE Finite Elements
    G.Tvs = zeros(G.nv,G.ns); % init
    G.dims = cell(G.ns,1);
    
    for iseg=1:G.ns
      mycond = C.PassiveConductor(iseg);
      inpoly = inpolygon(G.rv,G.zv,mycond.Rbound,mycond.Zbound);
      G.Tvs(:,iseg) = inpoly/sum(inpoly);
    
      % dims label
      dims_prefix = G.dimv{find(inpoly,1,'first')}(1:end-4);
      if iseg==1 || ~isequal(dims_prefix,prev_prefix)
        iii=1; % reset counter
      else
        iii=iii+1; % increment dims counter
      end
      prev_prefix = dims_prefix; % store prefix
      G.dims{iseg} = sprintf('%s_s%03d',dims_prefix,iii);
    end
    G.Tivs = double(G.Tvs~=0)'; % matrix with ones/zeros such that `Is = G.Tivs*Iv` (see JMM toolbox documentation)
  case 2
    % directly use vessel filament description derived from CREATE Finite Elements
    [G.rv,G.zv,G.wv,G.hv,G.Tvs,G.Rv,G.dims,G.dimv] = create2ves(C);
    G.Tivs = double(G.Tvs~=0)'; % matrix with ones/zeros such that `Is = G.Tivs*Iv` (see JMM toolbox documentation)
  otherwise
    error('unsupported vessel mode')
end

%% Gaps [r1 r2 z1 z2]
if isfield(P,'confxml') % take gaps from xml
  XMLConfig = xml_read(P.confxml);
  Gaps = XMLConfig.OutputParameters.Gaps;
  rzW = [];
  k = 0;
  for i = 1 : numel(Gaps.Gap)
    if Gaps.Gap(i).ATTRIBUTE.compute % also accepts values different from 1
      k = k+1;
      G.dimW{k,1} = Gaps.Gap(i).CONTENT;
      rzW(k,:) = [Gaps.Gap(i).ATTRIBUTE.r1 Gaps.Gap(i).ATTRIBUTE.r2 Gaps.Gap(i).ATTRIBUTE.z1 Gaps.Gap(i).ATTRIBUTE.z2];
    end
  end
else
  rzW = [ % Old case for tests fallback
    4.4155    4.7101   -3.4043   -3.8083
    5.2620    5.6389   -3.9852   -3.8511
    4.0970    6.0970   -1.5000   -1.5000
    4.0670    6.0670   -0.4836   -0.4836
    4.0460    6.0460    0.5328    0.5328
    4.0460    6.0460    1.5492    1.5492
    4.0760    6.0760    2.5656    2.5656
    4.1254    7.0813    3.5697    3.0570
    4.3159    6.6717    4.3272    2.4697
    4.9131    5.4667    4.7101    1.7616
    5.7602    4.5070    4.5265    1.8008
    5.7629    5.7629    4.5244    1.5244
    6.5925    4.6620    3.8891    1.5928
    7.4708    5.9591    3.0785    1.7691
    7.9349    6.1980    2.3978    1.4063
    8.2703    6.3467    1.6788    1.1314
    8.3950    6.3953    0.6328    0.5984
    8.3062    6.3685   -0.4191    0.0760
    7.9008    6.1501   -1.3379   -0.3710
    7.2865    5.8251   -2.2518   -0.8864
    6.2723    5.0458   -3.0419   -1.4621];
  G.dimW = [cellstr(num2str((1:2)','cg%d')) ; cellstr(num2str((1:9)','g%d')) ; cellstr(num2str((10:19)','g%d'))];
  G.dimW([17 14 12 5]) = cellstr(num2str((3:6)','cg%d'));
end
G.rW = rzW(:,1);
G.zW = rzW(:,3);
G.oW = atan2(rzW(:,4)-rzW(:,3)     ,rzW(:,1)-rzW(:,2)    );
G.aW = sqrt((rzW(:,2)-rzW(:,1)).^2+(rzW(:,4)-rzW(:,3)).^2);

end

function [Taa,dima] = connect_coils(dima,Taa,coil1,coil2,circname,series)
% [Taa,dima] = connect_coils(dima,Taa,coil1,coil2,circname,series)
%   Connect coils [coil1] and [coil2] in series/antiseries.
%   dima     = coil names
%   Taa      = coils connection matrix
%   coil1    = first coil name
%   coil2    = second coil name
%   circname = series circuit name
%   series   = 1 for series, -1 for antiseries

k = find(strcmp(dima,coil1));
l = find(strcmp(dima,coil2));
Taa(:,k) = Taa(:,k) + series*Taa(:,l); Taa(:,l) = [];
dima{k} = circname; dima(l) = [];
end

