function xc = gaps(x,dx)
% function xc = gaps(x,dx)
%
% find gaps in data, i.e. timeseries, delineated by unusually large
% differences in the value, i.e. exceeding (by 10%) the median difference
% or "dx" if specified
%
% x:    input data
% dx:   optional. maximum data difference
%
% xc:   data arranged into cell vector, delineated by gaps

% T.Hilmer, UH
% 2011.01.8  version 3
%   fixed bug for no gaps
% 2010.03.25  version 2
%   Drastically revised and simplified

if sum(size(x)~=1) > 1, error('haven''t coded for N-D yet'), end

Dx = abs(diff(x));

if nargin < 2
    dx = 1.10*nanmedian(Dx);
end

% "j" is the starting index for each section:
j = 1 + find(Dx > dx);

if isempty(j) % no gaps found.  trivial solution
    xc = {x};
    return
end

if j(1) ~= 1, j = [1; j]; end
% except for the last value, which ensures the tail section is handled
% correctly even if there is a gap between the last two values:
j(end+1) = length(x)+1; % that simple

for n = 1:length(j)-1
    xc{n} = x(j(n):j(n+1)-1);
end


return

% This was the old function. Preserve for posterity.

% function to break up gappy timeseries into subsets.  Identifies gaps
% based on a nifty sliding median difference between 't' values, or a
% specified dt.  Optional factor 'p' will disregard   gaps < median_gap*p
%
% t:    time vector in any units
% dt:   specified time difference, i.e. sampling rate
% p:    optional gap allowance factor. defaults to 1
%
% x:    indices into time vector of continuous subsets
% dt:   detected time difference: median difference between 't' values'
%
% there is a fundamental error to this method.  It identifies the median of
% the time intervals.  If you input a time vector composed of two sampling
% frequencies, e.g. 1 s and 5 min data, it is hard to say what will happen.
%
% I suggest using a p-value of ~0.9.  1 is too restrictive

if nargin < 4, showme = 0; end

if nargin < 3 || isempty(p), p = 1; end

td = diff(t);

if nargin < 2 || isempty(dt),
    % sliding median filter. For-loop for the each point in the vector.
    % Slow but accurate.
    Nd = length(td);
    dt = nan(Nd,1);
    for n = 1:Nd
        disp(sprintf('sliding median filter completion: %.3f',n/Nd))
        ind = n + (-10:10); ind(ind < 1 | ind > Nd) = [];
        dt(n) = median(td(ind));
    end    
end

f = [find(abs(td-dt) > (1-p)*dt); length(t)];

n = 1; m = 1;
while true
    
    x{m} = n : f(m);

    n = f(m) + 1;
    m = m + 1; 
    
    if n > length(t), break, end

end

if length([x{:}]) ~= length(t), error('function ''gaps'' is faulty'), end

if showme
    figure
    c = colormap(lines);
    h.a1 = subplot(211);
    plot(t,1:length(t),'ko')
    h.a2 = subplot(212);
    plot(t(1:end-1),diff(t)*86400,'ko')
    set([h.a1 h.a2],'NextPlot','add')

    m = 0;% m is a cyclical color index
    for n = 1:length(x)
        m = m + 1; if m > length(c), m = 1; end
        plot(h.a1,t(x{n}),x{n},'Linestyle','none','Marker','.','Color',c(m,:))
        plot(h.a2,t(x{n}(1:end-1)),diff(t(x{n}))*86400,'Linestyle','none','Marker','.','Color',c(m,:))
    end
    datetick(h.a1), datetick(h.a2)
end