function theResult = timterp(src, dst)

% timterp -- Interpolation of time data.
%  timterp(src, dst) interpolates the "src" data
%   at the times corresponding to the "dst" data.
%   The "src" and "dst" are NetCDF variables, both
%   of whose leftmost dimension represents the same
%   kind of time, possibly EPIC time with "time"
%   (Julian Day) augmented by "time2" (ms).
%   The result is a Matlab array of interpolated values,
%   oriented in the same fashion as the "dst" data.  If
%   no output argument is given, the result is silently
%   assigned to "ans" in the caller's workspace.


%%% START USGS BOILERPLATE -------------% There is no published documentation for these programs- use the 
%  internal comments and help function.
% Programs written in Matlab- as early as 7.2.0.232 (R2006a) and 
% most recent modifications in  v7.8.0 (R2009a)
% Programs ran on PC with Windows XP Professional OS, and RHEL4 linux.
%
% "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."
%%% END USGS BOILERPLATE --------------

 
% Fran Hotchkiss and Chuck Denham, USGS.
% Version of 22-Dec-1997 14:31:10.

if nargin < 2, help(mfilename), return, end

SEC_PER_DAY = 86400;   % 24 * 60 * 60.
MSEC_PER_DAY = 86400000;   % 24 * 60 * 60 * 1000.

% Time for "src".

fsrc = parent(src);
d = dim(src);
tsrc = fsrc{'time'};
t2src = fsrc{'time2'};
if ~isempty(t2src)   % EPIC.
    tsrc = EP_DateNum([tsrc(:) t2src(:)]);   % Decimal days.
elseif ~isempty(tsrc.base_date)   % Seconds + base_date.
    b = tsrc.base_date(:);
    t2src = datenum(b(1), b(2), b(3), b(4), b(5), b(6));
    tsrc = t2src + tsrc(:) / SEC_PER_DAY;
else
    tsrc = tsrc(:);   % Unknown units.
end

% Time for "dst".

fdst = parent(dst);
d = dim(dst);
tdst = fdst{'time'};
t2dst = fdst{'time2'};
if ~isempty(t2dst)   % EPIC.
    tdst = EP_DateNum([tdst(:) t2dst(:)]);   % Decimal days.
elseif ~isempty(tdst.base_date)   % Seconds + base_date.
    b = tdst.base_date(:);
    t2dst = datenum(b(1), b(2), b(3), b(4), b(5), b(6));
    tdst = t2dst + tdst(:) / SEC_PER_DAY;
else
    tdst = tdst(:);   % Unknown units.
end

% Set-up for Matlab "interp1" function.

x = tsrc;
y = src(:);
xi = tdst;

yi = interp1(x, y, xi, 'linear');

% Output.

if nargout > 0
    theResult = yi;
else
    assignin('caller', 'ans', yi)
end