function [rv,zv,wv,hv,Tvs,Rv,dims,dimv] = create2ves(C,doplot)
% [rv,zv,wv,hv,Tvs,Rv,dims,dimv] = create2ves(C,doplot)
% define vessel directly from CREATE segments
% C: CREATE vessel description
% doplot = true to plot
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if nargin==1, doplot=false; end

nelem = numel(C.PassiveConductor);
iv = 0;
for ielem = 1:nelem
  cond = C.PassiveConductor(ielem);
  rr = cond.Rbound; rr=[rr(end-1),rr]; %#ok<AGROW>
  zz = cond.Zbound; zz=[zz(end-1),zz]; %#ok<AGROW>
  
  % eliminate colinear points, keep only extrema
  col = nan(numel(cond.Rbound,1));
  for ipts=2:(numel(rr)-1)
    v1 = [rr(ipts)-rr(ipts-1);zz(ipts)-zz(ipts-1)];
    v2 = [rr(ipts+1)-rr(ipts);zz(ipts+1)-zz(ipts)];
    % colinearity indicator: 1 means they are 'in a row'
    col(ipts-1) = (v1'*v2)/(norm(v1)*norm(v2));
  end
  ielim = col>0.9;
  % rr,zz are corners of elements to keep
  rr =  cond.Rbound(~ielim); zz= cond.Zbound(~ielim);
  
  if doplot
    plot(rr,zz,'ok'); hold on;
    plot(cond.Rcenter,cond.Zcenter,'*k')
  end
  
  % use SVD to find principal direction
  XX = [rr;zz]-[cond.Rcenter;cond.Zcenter];
  [U,S,~] = svd(XX);
  aspect = S(1,1)/S(2,2); % aspect ratio (length/width) of segment
  
  % Store reults for this segment
  if aspect>2 % long segment, take a mean width
    npts = round(aspect); % number of points
    ivv=iv+(1:npts);
    xx = [linspace(-1,1,npts)*(npts)/(2*(npts+1));zeros(1,npts)];
    a=U*S(1:2,1:2)*xx + [cond.Rcenter;cond.Zcenter];
    rv(ivv,1) = a(1,:); zv(ivv,1) = a(2,:); %#ok<AGROW>
    wv(ivv,1) = sqrt(mean(sum([diff(rv(ivv)),diff(zv(ivv))].^2,2))); %#ok<AGROW>
    hv(ivv,1) = wv(ivv,1); %#ok<AGROW> % square filaments   
  else % very short segment, use single filament
    npts = 1; % override npts
    ivv=iv+(1:npts);
    rv(ivv,1) = cond.Rcenter; zv(ivv,1) = cond.Zcenter; %#ok<AGROW>
    wv(ivv,1) = max(abs(diff(XX(1,:))),[],2); %#ok<AGROW>
    hv(ivv,1) = max(abs(diff(XX(2,:))),[],2); %#ok<AGROW>
  end
  Tvs(ivv,ielem) = 1/npts; %#ok<AGROW> % assumes segment current spread evenly across segment %#ok<AGROW>
  iv = iv+npts;
  
  if doplot
    for ivi = 1:numel(ivv)
      iii = ivv(ivi);
      xx = rv(iii)+wv(iii)*[-1 -1 1 1]/2;
      yy = zv(iii)+hv(iii)*[-1 1 1 -1]/2;
      patch(xx,yy,'r');
      plot(rv(iii),zv(iii),'.b')
    end
  end
end

%% Reverse engineer Resistances Rv from element Resistances
% Rss = Tvs'*Rvv*Tvs (JMM toolbox formulae)
% assuming constant Rv for each filament in segment:
% Rss(is) = Tvs(:,1)'*I*Tvs(:,1)*Rv;
Rs = C.R(contains(C.x_type,'Passive')); % segment resistances
nv = size(Tvs,1);
Rv = zeros(nv,1);
for is = 1:numel(Rs)
  Rv((Tvs(:,is)~=0)) = Rs(is) / (Tvs(:,is)'*eye(nv)*Tvs(:,is));
end

%% Vessel element names
% Element subdivision per component
iin = 1:50; iout = 51:101; itria = 101:105; irail = 106:110;
dims = cell(nelem,1);
dimv = cell(numel(rv),1);
segcat = {'in',iin,'out',iout,'triangular_support',itria,'div_inboard_rail',irail};
for icat=1:2:numel(segcat)
  segstr = segcat{icat};
  iis = segcat{icat+1}; nns = numel(iis);
  dims(iis) = cellstr(num2str((1:nns)',[segstr,'_s%03d']));
  
  % for all segments of this type, assign vessel label
  ivv = any((Tvs(:,iis)~=0),2); nnv = sum(ivv);
  dimv(ivv) = cellstr(num2str((1:nnv)',[segstr,'_%03d']));
end
  
end



