function [spos,vpos,vstdpos,sneg,vneg,vstdneg,vcl] = split(s,v,vstd,Ri)
% split  Separate a velocity profile taken across a pipe into two profiles
%
% Synopsis:  [spos,vpos,vstdpos,sneg,vneg,vstdneg,vcl] = split(s,v,vstd,Ri)
%
% Input: s = vector of positions across the duct
%        v = velocity values at each s(i)
%        vstd = standar deviation of velocity readings at each v(i)
%        Ri = inner radius of the pipe
%
% Output: spos, sneg = value of s relative to the centerline of the pipe
%                      spos is the set of all positive s values.  sneg
%                      is the set of all negative values.  Both sets
%                      include s=0.  spos includes spos = Ri, and sneg
%                      includes sneg = -Ri.
%         vpos, vneg = velocity values corresponding to spos and sneg.
%                      The velocity value corresponding to s=0 is obtained
%                      by interpolation.  Values of vpos =0 and vneg = 0
%                      are included to account for spos = Ri and sneg = -Ri
%         vstdpos, vstdneg = standard deviation of velocity values corresponding
%                      to spos and sneg.  See vpos and vneg.
%         vcl = value of velocity at the centerline (s=0).  vcl is obtained
%               by linear interpolation in the original data set.

verbose = false;

% --- Sort the data in order of increasing s
[s,is] = sort(s);   v = v(is);  vstd = vstd(is);

% --- Interpolate to find velocity and standard deviation of velocity
%     on the centerline (s=0)
vcl = interp1(s,v,0);
vstdcl = interp1(s,vstd,0);

% --- Find index of point with negative s that is closest to zero
imid = max( find(s<0) );

% --- Copy data into two sets:  one for s<0 and one for s>0
%     Add centerline values to each set
spos = [0; s(imid+1:end)];
vpos = [vcl; v(imid+1:end)];
vstdpos = [vstdcl; vstd(imid+1:end)];
sneg = [s(1:imid); 0];
vneg = [v(1:imid); vcl];
vstdneg = [vstd(1:imid); vstdcl];

% --- Trim data by eliminating any points with |spos|>Ri or |sneg|>Ri
%     and add data point corresonding to the no slip condition.
imax = min( find(spos>Ri) );
if ~isempty(imax)
  spos(imax:end) = [];   vpos(imax:end) = [];  vstdpos(imax:end) = [];
  if verbose
    fprintf('\nspos values out of bounds:\n');
    disp(spos(imax+1:end));
  end
end
spos = [spos; Ri];    vpos = [vpos; 0];     vstdpos = [vstdpos; 0];

imin = max( find(abs(sneg)>Ri) );
if ~isempty(imin)
  sneg(1:imin) = [];   vneg(1:imin) = [];  vstdneg(1:imin) = [];
  if verbose
    fprintf('\nsneg values out of bounds:\n');
    disp(sneg(1:imin));
  end
end
sneg = [-Ri; sneg];    vneg = [0; vneg];     vstdneg = [0; vstdneg];

if verbose
  imax = min( find(abs(sneg)>Ri) )
  fprintf('\nsneg values out of bounds:\n');
  disp(sneg(imax+1:end))
end
