function [rv,zv,wv,hv,Rv,nv,dimv] = ITERvessel2filaments(vessel)
% convert ITER vessel description to filaments
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
[rv,zv,wv,Rv,dimv] = deal([]);
for ives=fieldnames(vessel)'
  [rrv,zzv,wwv,RRv] = vessel2filaments(vessel.(ives{:}));
  rv=[rv;rrv]; zv=[zv;zzv]; wv=[wv;wwv]; Rv = [Rv;RRv]; %#ok<AGROW>
  nnv = numel(rrv); dimv = [dimv;cellstr(num2str((1:nnv)',[ives{:},'_%03d']))]; %#ok<AGROW>
end
hv = wv; % squares
nv = numel(rv);
end

function [rv,zv,wv,Rv] = vessel2filaments(insrc)
[rv,zv] = filaments(insrc.r,insrc.z,insrc.thick);
wv = insrc.thick*ones(size(rv));
rho = insrc.R*1e-6; % resistivity [Ohm m]
Rv = rho * 2*pi*rv./wv.^2; % vessel filament resistance
end

function [rv,zv] = filaments(r,z,w)
% create vessel filaments spaced as much as the specified width.

% determine if closed or not
isclosed = (r(end)-r(1)).^2+(z(end)-z(1)).^2 < 1.1*mean(sum(diff([r,z]).^2,2));
isclosed = isclosed && numel(r)>2;

if isclosed
  rrr = [r;r(1)]; zzz=[z;z(1)];
else
  rrr = r; zzz=z;
end

dr = diff(rrr); dz=diff(zzz);
ds = sqrt(dr.^2+dz.^2);
s = [0;cumsum(ds)];
snew = linspace(w/2,s(end)-w/2,ceil(s(end)/(w)));
rv = interp1(s,rrr,snew)';
zv = interp1(s,zzz,snew)';

end
