function [RL, RV] = runlength(D,zrhos)
%
%function [RL, RV] = runlength(D,zrhos)
%
% Overturn Toolbox
% David S. Ciochetto
% Dalhousie University Physical Oceanography
% 2005 Mar 21 (Last Modified)
%
% This function is to calculate the statistical property known as runlength
% given in Galbraith and Kelley (1996).
%
% "The standard statistical property known as run length is defined as
% follows. The points in a time series are examined sequentially and
% adjacent values of one sign are grouped into sets called 'runs'. For a
% random, un-correlated time series with equal numbers of positive and
% negative values, the probability density function (PDF) of run length is
% P(n) = 2^-n; where P(n) is the probability of a run of length n."
%
% INPUTS:
%   D  =  A single vector of data (row or column)
%   zrhos (optional) = a string ['y'/'n'] indicating the users desire to
%         insert the number of zeros found in D back into the center of the
%         result. Default value = 'n'.
%
% OUTPUTS:
%   RL = A vector of data with the run lengths for negative values listed
%        first, any zero values listed in the middle and the run lengths
%        for the positive values listed last.
%   RV = A vector of the same length as D with the elements = the value of
%        the runlength, if any, matching the indices. Zeros are included
%        where rho(z) = rho^(z).
%
% NOTES:
%   1. The theoretical PDF works for random noise not Gaussian distribted
%   noise.
%
%   2. This is meant for primary use on Thorpe fluctuation profiles which
%   are rho(z) - rho^(z) where rho(z) is a density profile and rho^(z) is
%   that profile sorted such that it is gravitationally stable. Other data
%   sent should be mean zero data.
%
% REFERENCES:
%   Galbraith, P. S. and Kelley, D. E., 1996 "Identifying overturns in CTD
%   profiles," Journal of Atmospheric	and Oceanic Technology, 13, 688-702.
%
% USES:
%
% SEE ALSO:
%	  test_runlength.m, ot_cutoff.m, overturns.m, ot_example.m,
%

if ~exist('zrhos')
	zrhos = 'n';
end

D=D(:);	% Convert to a column if not.

r1 = [0; find(D>=0); length(D)+1];
r2 = diff(r1);
r3 = r2(r2>1)-1;
% Gives (-) run lengths

r4 = [0; find(D<=0); length(D)+1];
r5 = diff(r4);
r6 = r5(r5>1)-1;
% Gives (+) run lengths

% For code development and change only. Check if it was done correctly.
% by plotting the original fluctuation signal, the positive fluctuations
% and the negative fluctuations. Note that they get flipped when taken.
% I=[1:length(D)];plot(I,D,'.-k');line(get(gca,'xlim'),[0 0],'color','r');hold on;
% plot(I(r1(2:end-1)),D(r1(2:end-1)),'c+:',I(r4(2:end-1)),D(r4(2:end-1)),'m+:');
% hold off;zoom on;clear I rng

% Check for zeros
z = sum(D==0);

% When we add zeros and calculate the histogram later, the zeros weight the
% total number of runs and skew the percentages so they will not be
% returned here in the result. Will not need to know the order of the runs
% for overturn analysis, just the distribution.

if lower(zrhos)=='y'	% Add zeros into the middle of the result
	% The Result [(-) 0 0 (+)]
	RL = [r3; zeros(z,1); r6];
else
	% The Result [(-) (+)]
	RL = [r3; r6];
end

% Check the result for consistant length
if sum(r3)+sum(r6)+z ~= length(D)
	disp('Problem; the lengths in runlength.m do not match')
end

% Make a vector with the same length as D with the elements = the value of the runlength,
% if any, matching the indices. Zeros are included where rho(z) = rho^(z).
RV = zeros(size(D));
inegRL = setdiff([1:length(D)]',r1); % Gives the indicies for the (-) values
iposRL = setdiff([1:length(D)]',r4); % Gives the indicies for the (+) values
% Make a vector of the same length as inegRL and fill it with the proper RL values to match
nRV=[]; clear k
for k=1:length(r3)
  nRV=[nRV; r3(k)*ones(r3(k),1)];
end
% Make a vector of the same length as iposRL and fill it with the proper RL values to match
pRV=[]; clear k
for k=1:length(r6)
  pRV=[pRV; r6(k)*ones(r6(k),1)];
end
RV(inegRL) = nRV;
RV(iposRL) = pRV;
