function [interpData] = smart_interp(rawTime, rawData, interpTime, maxGap);
%
% [interpData] = smart_interp(rawTime, rawData, interpTime, maxGap);
%
% Function which interpolates an original data set (rawTime, rawData,
% with rawTime in julian days) to a new time base (interpTime); maxGap is
% the maximum time (IN HOURS) that can exist in the raw data set (as
% NaN's or missing data) over which smart_interp will interpolate. The
% program splits the data into groups if there are gaps.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is self described.
% Program written in Matlab v7.1.0 SP3
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty,
% expressed or implied, is made by the USGS or the United States
% Government as to the accuracy and functioning of the program
% and related program material nor shall the fact of distribution
% constitute any such warranty, and no responsibility is assumed
% by the USGS in connection therewith."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Check in the input interpTime
checkTime = diff(interpTime);
if (checkTime > (mean(checkTime) + std(checkTime))) | ...
(checkTime < (mean(checkTime) - std(checkTime))) | checkTime < 0
error('The time base to interpolate to must be monotonically increasing!');
return
end
%Set up the empty data set
interpData = repmat(NaN, size(interpTime));
%Eliminates NaN's from the data
theNans = find(isnan(rawData));
rawData(theNans) = [];
rawTime(theNans) = [];
if isempty(rawData)
interpData = repmat(NaN, size(interpTime));
return
end
%Insure the data is in order
[rawTime, I] = sort(rawTime);
rawData = rawData(I);
%Delete any repeats
diffTime = diff(rawTime);
theRepeats = find(diffTime == 0);
rawTime(theRepeats) = [];
rawData(theRepeats) = [];
%Convert maxGap to Julian Days
maxGap = maxGap/24;
%Find the large gaps in the data set
timeDiff = diff(rawTime);
largeGaps = find(timeDiff > maxGap);
largeGaps = largeGaps(:);
theStarts = [1];
theEnds = length(rawTime);
%Group the data into sets
theStarts = [theStarts(:); (largeGaps+1)];
theEnds = [largeGaps; theEnds(:)];
%Eliminate sets of one
theBads = find((theEnds - theStarts) == 1);
theStarts(theBads) = [];
theEnds(theBads) = [];
%Go through and do the interpolations
for indexSet = 1:length(theStarts)
setStart = rawTime(theStarts(indexSet));
setEnd = rawTime(theEnds(indexSet));
startInterp = min(find(interpTime >= setStart));
endInterp = max(find(interpTime <= setEnd));
%Outside interpolate time
if isempty(startInterp) | isempty(endInterp) | ...
(theStarts(indexSet) == theEnds(indexSet))
continue
end
interpData(startInterp:endInterp) = ...
interp1(rawTime(theStarts(indexSet):theEnds(indexSet)), ...
rawData(theStarts(indexSet):theEnds(indexSet)), ...
interpTime(startInterp:endInterp));
end