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