find_gaps.m

function [gap, comment, errors] = find_gaps(GPS,seq,stat)

% Usage: [gap, comment, errors] = find_gaps(GPS,seq,stat)
%
% Scans unwrapped NIMSread sequence string for gaps.
% Outputs indices of the last block before each gap
% and the number of missing scans.
% For efficiency, can also output the time series start time.
%
% Anna Kelbert, 24 Aug 2007
%
% Last modified 11 Aug 2020 to identify spikes if they are not already
% removed by the removeSpike() subroutine - hopefully they will all be
% identified and fixed on the input to Matlab so this condition will no
% longer be encountered (double, consecutive gaps that sum to the length
% of 255 and are associated with a spike in the data). More importantly, 
% no longer issue an error if the GPS string is available on both sides
% of a gap, even if the sequence numbers don't match. This will remove
% false errors, allowing us to only focus on those that matter.

errors = '';

j = 1;
gap.ind = [];
gap.len = [];
gap.err = [];
gap.spike = [];
gap.start = {};
gap.end = {};
for i = 2:length(seq)
    diff = seq(i) - seq(i-1);
    if diff ~= 1
        gap.ind(j) = i-1;
        gap.len(j) = diff - 1;
        gap.err(j) = 1;
        gap.spike(j) = 0;
        j = j+1;
    end
end

% In some special cases, e.g. when a duplicate data point
% has been deleted, the gap is zero length; still worth
% checking that everything is consistent!
gap.len = rem(gap.len,256);

if nargin < 3
    comment = 'Not enough information to compute gaps. ';
    disp(comment);
    return
elseif isempty(gap.ind)
    comment = 'No data gaps. ';
    disp(comment);
    return
else
    comment = ['Found data gaps (' int2str(j-1) '). '];
    errors = [errors comment];
    disp(comment);    
end

ndupl = length(find(gap.len==0));
if ndupl > 0
    comment = ['Of these, ' num2str(ndupl) ' duplicate blocks [deleted]. '];
    errors = [errors comment];
end
    
gap.ind(j) = length(seq);
gap.len(j) = 0;
nsecs = [];
num2sec = 24*3600;
sec2num = 1/(3600*24);
section = {};

% Compute the start and end times for the first section of the time series
ind(1) = 1;
ind(2) = gap.ind(1);
GPSj = GPS(ind(1):ind(2));
seqj = seq(ind(1):ind(2)) - seq(ind(1));
statj = stat(ind(1):ind(2));
disp('Searching for start time of data section 1 ...');
t = startTimeFromGPS(GPSj,seqj,statj);

if ~isempty(t)
    disp(['Start time of section 1 is ' datestr(t)]);
    section(1).start = datenum(t);
    section(1).nscans = gap.ind(1);
    section(1).end = section(1).start + sec2num*section(1).nscans;
else
    comment = ['Unable to compute start time because of a gap ' ...
        'in sequence numbers before the first valid GPS message.' ...
        'Please delete all data up to and including index ' int2str(gap.ind(1)) '. '];
    errors = [errors comment];
    disp(comment)
end 

bad = [];
good = [];
for j = 1:length(gap.ind)-1
    ind(1) = gap.ind(j) + 1; % index of section start
    ind(2) = gap.ind(j+1);   % index of section end
    GPSj = GPS(ind(1):ind(2));
    seqj = seq(ind(1):ind(2)) - seq(ind(1));
    statj = stat(ind(1):ind(2));
    tprev = t;
    disp(['Searching for start time of data section ' int2str(j+1) ' ...']);
    t = startTimeFromGPS(GPSj,seqj,statj);
    if ~ isempty(t)
        disp(['Start time of section ' int2str(j+1) ' is ' datestr(t)]);
        section(j+1).start = datenum(t);
        section(j+1).nscans = gap.ind(j+1) - gap.ind(j);
        section(j+1).end = section(j+1).start + sec2num*section(j+1).nscans;
        if ~ isempty(tprev)
            nsecs(j) = int32(num2sec*(section(j+1).start - section(j).end));
            % Check whether the time difference is consistent with the sequence numbers
            if gap.len(j) == rem(nsecs(j),256)
                gap.len(j) = nsecs(j);
                gap.err(j) = 0;
                gap.end{j} = datestr(t);
                gap.start{j} = datestr(datenum(t)-nsecs(j)/24/3600);
                if gap.len(j) > 0 % Don't count duplicate data blocks
                    good = [good j];
                end
            else
                comment = ['The difference in sequence numbers (nscans = ' int2str(gap.len(j)) ') ' ...
                    'for gap # ' int2str(j) ' does not match the GPS time difference ' ...
                    '(nsecs = ' int2str(nsecs(j)) '). Using the GPS gap length. '];
                disp(comment)
                % trust the GPS gap info in this case - usually not a huge
                % problem as long as the GPS is good, so no longer an error!
                gap.len(j) = nsecs(j);
                gap.err(j) = 0;
                gap.end{j} = datestr(t);
                gap.start{j} = datestr(datenum(t)-nsecs(j)/24/3600);
                good = [good j]; % this is still good...
            end
        elseif j>1
            if (gap.ind(j) - gap.ind(j-1) == 1) && (gap.len(j) + gap.len(j-1) == 255) && (stat(gap.ind(j)) ~= 129)
                comment = ['False double gap identified at ' int2str(gap.ind(j)) '. Using the GPS time to count back.'];
                disp(comment)
                % Identify spikes associated with false double gaps, and remove the gaps
                % safely since the GPS string is available!
                gap.len(j-1) = 0;
                gap.len(j) = 1;
                gap.err(j) = 0; % note that gap.err(j-1) = 1 still
                gap.spike(j) = 1;
                gap.end{j} = datestr(t);
                gap.start{j} = datestr(datenum(t)-1/24/3600);
                good = [good j]; % even this is still good!
            end
        else
            comment = ['Unknown length: gap ' int2str(j) ' [' int2str(gap.ind(j)) '] - no start time. '];
            bad = [bad j];
            disp(comment)      
        end
    elseif j>1
        if (gap.ind(j) - gap.ind(j-1) == 1) && (gap.len(j) + gap.len(j-1) == 255) && (stat(gap.ind(j)) ~= 129)
            comment = ['False double gap identified at ' int2str(gap.ind(j)) '. Removed.'];
            disp(comment)
            % Identify spikes associated with false double gaps, and remove the gaps
            % but issue an error - no GPS to check against!
            gap.len(j-1) = 0;
            gap.len(j) = 1;
            gap.err(j-1) = 1;
            gap.err(j) = 1;
            gap.spike(j) = 1;
            bad = [bad j]; % no longer good...
        end
    else
        comment = ['Unknown length: gap ' int2str(j) ' [' int2str(gap.ind(j)) ']. '];
        bad = [bad j];
        disp(comment)
    end
end

gap.ind = gap.ind(1:end-1);
gap.len = gap.len(1:end-1);

% bad = find(gap.err==1);
if ~isempty(bad)
    bad_gaps = num2str(gap.ind(bad),'%d; ');
    bad_gaps = bad_gaps(1:end-1);
    comment = ['Gaps of unknown length: ' int2str(length(bad)) ' [' bad_gaps ']. '];
    errors = [errors comment];
end

if ~isempty(good)
    good_gaps = num2str(gap.ind(good),'%d; ');
    good_gaps = good_gaps(1:end-1);
    comment = ['GPS used to compute length of ' int2str(length(good)) ' gaps [' good_gaps ']. '];
    errors = [errors comment];
end