function [t2,theta2,rho2] = interval_avg_vectors(t1, theta1, rho1, tbins)
%{
computes vector average for high-rate angle/magnitude data over
intervals tbins.  Converts angle in (theta)/magnitude(rho) to uv components
before averaging, then converts result back to polar.

Theta input/output is in degrees and is assumed to be compass
angle since the sd to uv functions are used to convert.

All inputs must be 1-D, tbins is interval bin edges.

Times are assumed to mark the start point of each interval.
HOWEVER, THE LAST VALUE OF tbins IS ASSUMED TO BE THE STOP POINT AND
NO AVG WILL BE COMPUTED FOR THIS VALUE.  THE RESULT ARRAY WILL
THEREFORE HAVE length(tbins)-1 ELEMENTS OR ROWS.

Uniform delta-t is not required in either t1 or tbins, although in most
applications delta-tbins will be uniform.  Ranges of t1 and tbins must overlap,
but the start and end points for each time series may be different.
Intervals in tbins outside the range of t1 are set to NaN.

Returns tbins(1:end-1) as timestamp variable for output
%}

% check inputs for validity
test_t1 = (size(t1,1)~=1 && size(t1,2)~=1);
test_t2 = (size(tbins,1)~=1 && size(tbins,2)~=1);
if (test_t1 || test_t2)
    error('interval_avg_vectors: time arrays must be 1-D');
end
test_x1 = (size(theta1,1)~=1 && size(theta1,2)~=1);
test_x2 = (size(rho1,1)~=1 && size(rho1,2)~=1);
if (test_x1 || test_x2)
    error('interval_avg_vectors: data arrays must be 1-D');
end
theta1 = theta1(:);
rho1 = rho1(:);
if length(t1) ~= size(theta1,1)
    error('interval_avg_vectors: t1 and x1 not equal length');
end
if size(theta1) ~= size(rho1)
    error('interval_avg_vectors: x1 and x2 not the same size');
end
if mean(diff(tbins)) <= mean(diff(t1))
    error('interval_avg_vectors: delta tbins should be > delta t1');
end
if tbins(1) > t1(end) || tbins(end) < t1(1)
    error('interval_avg_vectors: intervals do not overlap');
end

rows = length(tbins)-1;
u2 = NaN(rows,1);
v2 = NaN(rows,1);

[u1,v1] = sd_to_uv(rho1,theta1);

% We operate over intervals in tbins, but begin and end times
% of tbins may differ from begin and end of t1.
% Therefore, find the range of tbins over which we will operate.
% Start point in tbins:
ii = find(tbins <= t1(1),1,'last');   % find last t2  <= t1(1)
if isempty(ii)
    startPt = 1;     % t2 starts later than t1, so start at t2(1)
else
    startPt = ii;    % beginning of first interval is tt2[ii[-1]]
end
    
% End point in tbins will be the end of the last interval
ii = find(tbins >= t1(end),1,'first'); % find first value in tbins >= t1(end)
if isempty(ii)
    endPt = length(tbins);   % tbins ends before t1, so end at tbins(end)
else
    endPt = ii;        % end of last interval is tbins(ii(1))
end

intvl = startPt:endPt;    % indices in tbins for valid bin edges
zz = length(intvl)-1;     % number of bin aves to compute
% compute bin covariances
for jj=1:zz
    [~,ii] = min(abs(t1-tbins(intvl(jj))));  % last t1 <= start of intvl
    kk = find(t1<tbins(intvl(jj)+1),1,'last'); % last t1 < next intvl
    if ii==kk  % no data for this bin
        continue
    end
    u2(intvl(jj),:) = nanmean1(u1(ii:kk,:));
    v2(intvl(jj),:) = nanmean1(v1(ii:kk,:));
end

[rho2,theta2] = uv_to_sd(u2,v2);
t2 = tbins(1:end-1);

end
