function [xint, shiftx, shiftix, shifty, shiftiy] = shiftinterp(tx,x, ty,y, nlagx, defaultshift)
% Shift high-frequency x(tx) to times of lower frequency ty
% by shifts <= dt*nlag to maximize autocorrelation of x and y.
% xint is the x quantity shifted and interpolated.
% shifty, shiftx are time shift (in days) on the x and y grid, respectively.
%
% Simon de Szoeke, PISTON

y=y(:); x=x(:);
nx = length(x);
ny = length(y);
dtx = mean( diff(tx) );
dty = mean( diff(ty) );

% interpolate y to the higher frequency x times tx
yint = interp1(ty(:),y(:), tx(:), 'linear');

weightx = [1:nlagx+1 nlagx:-1:1]';
% index  = 1:(2*nlag+1);
nwin = round(nx/nlagx);
ctrwinx = round( linspace(1,nx,nwin) )'; % center of window index
stawinx = max(1, ctrwinx - nlagx     );  % start  of window
endwinx = min(   ctrwinx + nlagx, nx );  % end
% ctrwiny = interp1( ty, 1:ny, tx(ctrwinx), 'nearest','extrap');
% stawiny = interp1( ty, 1:ny, tx(stawinx), 'nearest','extrap');
% endwiny = interp1( ty, 1:ny, tx(endwinx), 'nearest','extrap');

[maxcor, shift] = deal( zeros(nwin,1) );
for iw = 1:nwin
    iseriesx = (ctrwinx(iw) - nlagx):(ctrwinx(iw) + nlagx);
    indexx = iseriesx > 0 & iseriesx <= nx;
%     iseriesy = (stawinx(iw):endwinx(iw) ) - stawinx(iw) + 1;
%     indexy = iseriesy > 0 & iseriesy <= ny;
    y4cov = yint( stawinx(iw):endwinx(iw) ) .* weightx(indexx);
    x4cov =    x( stawinx(iw):endwinx(iw) ) .* weightx(indexx);
    y4cov = y4cov - nanmean( y4cov );
    x4cov = x4cov - nanmean( x4cov );
    y4cov( isnan(y4cov) ) = 0;
    x4cov( isnan(x4cov) ) = 0;
    XYcor = xcov( x4cov(:), y4cov(:), nlagx, 'coeff');
    [maxcor(iw), i] = max(XYcor);
    shift(iw) = i - (nlagx + 1); % subscript of x shift for each window
end

twin = tx( ctrwinx(:) ); % time at center of each window
mc = maxcor>0.5;
if sum(mc) > 1
    shiftx = interp1(twin(mc), dtx*shift(mc), tx, 'linear','extrap'); % time shifts on the x coordinate
    shifty = interp1(twin(mc), dtx*shift(mc), ty, 'linear','extrap'); % time shifts on the y coordinate
elseif sum(mc) == 1
    shiftx = dtx*shift(mc) + 0*tx;
    shifty = dtx*shift(mc) + 0*ty;
else
    shiftx = defaultshift  + 0*tx;
    shifty = defaultshift  + 0*ty;
end
xint = interp1(tx - shiftx, x, ty, 'linear', 'extrap');
xint(isnan(y)) = NaN;
shiftix = shiftx/dtx;
shiftiy = shifty/dty;
