function error = LagError(k, data, t, a, b);

% error = LagError(k, data, t, a, b);
%
%   This function computes the error for (data(t,a) - k(1)*data(t-k(2),b))^2
%   For temperature and conductivity, we have:
%
%     t = m_present_time
%     a = water_temp
%     b = water_cond
%     k = [12.5 -0.35]  (best first guess)
%
%   This function does NOT assume the parameters have been sampled on a
%   uniform time interval.  It resamples the parameters onto a uniform
%   time-base using MATLABs INTERP1Q function.  This makes the function
%   a little slow, but hopefully more accurate.
%
%   This function high-pass filters the data before looking for the error
%   since we are looking for a short lag time.  It uses the fft function
%   for the filtering, by simply truncating the low frequency component.
%
%   J. G. Bellingham
%   December 27, 1998
%

% Determine bounds of data.  Don't use first or last 45 seconds...
lastn = length(data(:, t));
t1 = data(5*45, t);
t2 = data(lastn-5*45, t);

% Establish new time-base with uniform interval
time = (t1:0.2:t2)';
n = length(time);

% Resample parameters on new (uniform) time-base
ok = find(~isnan(data(:,a)));
x1 = interp1q(data(ok,t), data(ok,a), time);

% plot(data(:, t), data(:,a), 'b', time, x1, '+r');

% Add lag...
ok = find(~isnan(data(:,b)));
x2 = interp1q(data(ok,t), data(ok,b), (time-k(2)));

% High pass filter each of the parameters

center = ((1:n)*0)';
n4 = round(n/4);
center(n4:(n-n4)) = 1;
x1 = ifft(fft(x1) .* center);
x2 = ifft(fft(x2) .* center);

% Computer error between parameters...
error = std(x1 - k(1)*x2);

return;

