function [interpData] = smart_interp_pres(rawTime, rawData, interpTime, maxGap);
%
%   [interpData] = smart_interp_pres(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.  Removes the
%   mean of the pressure data

%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) = [];

%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+1) = [];
rawData(theRepeats+1) = [];

%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); continue; end

    inData = interp1(rawTime(theStarts(indexSet):theEnds(indexSet)), ...
        rawData(theStarts(indexSet):theEnds(indexSet)), ...
        interpTime(startInterp:endInterp));
    
    interpData(startInterp:endInterp) = inData - nanmean(inData);
end