function [rv,zv,wv,hv,Tivs,Rs,Rv,dims,dimv] = vertices2meq(V,doplot)
% [rv,zv,wv,hv,Tvs,Rs,Rv,dims,dimv] = vertices2meq(V,doplot)
% define square vessel filaments from vertices of quadrilaterals.
%
% V: Array of parallellogram vessel element description structures
% V(i).Rvertex   R coordinate of vertex [m] (4x1)
% V(i).Zvertex   Z coordinate of vertex [m] (4x1)
% V(i).Rcenter   (optional) R coordinate of element center
% V(i).Zcenter   (optional) Z coordinate of element center
% V(i).Res       Resistance of element [Ohm] [optional - can also specify Resy]
% V(i).Resy      Resistivity of element [Ohm*m] [optional - can also specify Res]
%                NB: If both Res and Resy are specified, Res is used.
%
% set doplot = true to plot
%
% Outputs: documented in meqg.m
% 
% Method:
%   Find the longest edge, define the aspect ratio (long height / short height);
%   The number of filaments will be floor(aspectratio).
%   The width of the square filaments are chosen to maintain the total area
%   The filament centers are placed on the line joining the midpoints of the
%   two shortest edges, equally spaced.
%   
%   If total element resistance is given, the resistance of each filament
%   is chosen so that a constant voltage yields the same total current
%   in the filaments as in the original single element.
%  
%   If the resistivity is given, the resistance of each filament is
%   directly calculated as Res=resistivity*l/A
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if nargin==1, doplot=false; end

nelem = numel(V);
 % preallocate
nvmax = 2000; % max reasonable vessel elements
Tivs = zeros( nelem, nvmax );
Tvs  = zeros( nvmax, nelem );
Rv   = zeros( nvmax, 1     );
[rv,zv,wv,hv] = deal(zeros(nvmax,1));

iv = 0;
for ielem = 1:nelem
  cond = V(ielem);
  % rr,zz are corners of elements to keep
  rr =  cond.Rvertex; zz= cond.Zvertex;
  if iscolumn(rr), rr=rr'; end
  if iscolumn(zz), zz=zz'; end
  assert(~isempty(which('polyshape')),'this requires polyshape function. available >2017b')
  poly = polyshape(rr,zz);
  assert(size(poly.Vertices,1)==4,'only works for quadrilaterals')
  if ~isfield(cond,'Rcenter')
    % get centroid
    [cond.Rcenter,cond.Zcenter] = poly.centroid;
  end
  
  if doplot
    plot(rr,zz,'ok'); hold on;
    plot(poly)
    plot(cond.Rcenter,cond.Zcenter,'*k'); axis equal;
  end
  
  % find the longest edge
  [base,ibase] = max(sqrt(sum(diff(poly.Vertices([1:end,1],:)).^2,2)));
  % find aspect ratio as if it was a rectangle
  aspect = base.^2/poly.area;
  
  % number of filaments to use
  nfil = ceil(aspect-100*eps); % round upward to avoid gaps, avoiding roundoff errors
  ivv=iv+(1:nfil);             % update index
  N = 4; wrapN = @(x) (1 + mod(x-1, N)); % aux function for wrapping indices
  
  % take vertices of the two edges adjacent to the longest one
  Xshort1 = poly.Vertices(wrapN([ibase-1,ibase]),:);
  Xshort2 = poly.Vertices(wrapN([ibase+1,ibase+2]),:);
  
  % vectors of these same edges
  short1_midpoint = mean(Xshort1);
  short2_midpoint = mean(Xshort2);
  
  fracint = (1:2:(2*nfil))/(2*nfil); % fractional spacing of centroids along chord
  % vector between short edge midpoints, oriented towards center from short1_midpoint
  vv = short2_midpoint - short1_midpoint;
  vv = vv*sign(dot(short1_midpoint-[cond.Rcenter,cond.Zcenter],vv));
  
  % place rectangle midpoints
  rv(ivv) = short1_midpoint(1) - fracint'*vv(1);
  zv(ivv) = short1_midpoint(2) - fracint'*vv(2);
  
  % filament width to maintain total surface
  wv(ivv) = sqrt(poly.area/nfil);
  hv(ivv) = wv(ivv); % square filaments
  if ivv(end)>nvmax, warning('overran preallocation of %d vessel elements'); end
  
  % connections
  T = ones(1,nfil); % connection matrix
  % Ohm's law
  %Yv = 1 ./ Rv; % 1/Resistivity
  %Iv(1) = Yv(1)/sum(Yv(1:n)) Is(1) ,... => Iv = Yvv*T'*inv(T*Yvv*T')*Is
  %                                         Iv = Tvs                 *Is
  % with Tvs=Yvv*T'*inv(T*Yvv*T');
  
  % assuming constant resistivity, Rv~R(radius) so Yv~1/Radius
  Yvv = diag(1./rv(ivv)); % constant scaling does not matter
  myTvs = Yvv*T'/(T*Yvv*T');
  myTsv = myTvs';

  % Resistance per vessel element
  if isfield(cond,'Res') && ~isempty(cond.Res)
    % resistance given
    k = myTsv*diag(1./rv(ivv))*myTvs; % scaling
    Rv(ivv) = cond.Res/k * 1./rv(ivv);
  elseif isfield(cond,'Resy') && ~isempty(cond.Resy)
    % resistivity given: R=resy*L/A;
    Rv(ivv)   = cond.Resy * 2*pi*rv(ivv) ./ (wv(ivv).*hv(ivv));
  else
    error('must specify either resistance or resistivity')
  end
  Tivs(ielem,ivv) = T;
  Tvs(ivv,ielem) = myTvs;
  
  iv = iv+nfil;
  
  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','FaceAlpha',0.3);
      plot(rv(iii),zv(iii),'.b')
    end
    axis equal;
    drawnow;
  end
end
% clip to final dimension
nv = ivv(end);
Tvs = Tvs(1:nv,:);
Tivs = Tivs(:,1:nv);
Rv = Rv(1:nv);
rv = rv(1:nv); zv = zv(1:nv);
wv = wv(1:nv); hv = hv(1:nv);

Rs = diag(Tvs'*diag(Rv)*Tvs);

%% Vessel element names
% Element subdivision per component
dims = cell(nelem,1);
dimv = cell(numel(rv),1);
for idims = 1:numel(dims)
  % string for segment
  segstr = sprintf('Segment%03d',idims);
  dims{idims} = segstr;
  % for all segments of this type, assign vessel label
  ivv = any((Tvs(:,idims)~=0),2); nnv = sum(ivv);
  dimv(ivv) = cellstr(num2str((1:nnv)',[segstr,'_%03d']));
end
  
end