function [lag,i1,i2,co]=bestlag(a1,a2,nlag,npoint)
% function [lag,i1,i2,co]=bestlag(a1,a2,nlag,npoint)
%
% function to find the best shift for two vectors
% uses median difference
% 
%  input:  a1: first vector
%          a2: second vector
%          nlag = number of point to try shifiting , or vector
%
% M. Visbeck LDEO August-2002

if nargin<3, nlag=fix((length(a2)+length(a1))/8); end
if nargin<4, npoint=fix((length(a2)+length(a1))/2); end 
a1=a1(:);
a2=a2(:);

l1=length(a1);
l2=length(a2);
lmin=min([l1, l2]);

% set up lag index 
if length(nlag)>1
 ilag=nlag;
else
 ilag=-nlag:nlag; 
end

% set up target index 
lag1=min(ilag);
lag2=max(ilag);
ic=(max([1,-lag1])+1):(lmin-max([0,lag2])-1);

% reference matrix
c1=meshgrid(a1(ic),ilag)';

% create lag matrix 
for i=1:length(ilag)
   icc=ic+ilag(i);
   c2(:,i)=a2(icc);
end

% remove NaN and other bad data
igood=find(isfinite(sum(c1')+sum(c2')));

if length(igood)<10
 disp(' not enough data to determine lag')
 lag=0;
 co=1;
else
 % check only npoint points

 % use median as estimator
 cm=medianan(abs(gradient(c2(igood,:)')'-gradient(c1(igood,:)')'),ceil(length(igood)/6));
 [m1,im1]=min(cm);

 % find best match by fitting parabola to result
 il=floor(length(ilag)/4);
 ifit=[-il:il]+im1;
 ii=find(ifit<1 | ifit>=length(ilag));
 ifit(ii)=[];
 if length(ifit>10);
  c=polyfit(ilag(ifit),cm(ifit),2);
  cmb=cm+mean(cm)*2;
  cmf=cm+nan;
  cmf(ifit)=polyval(c,ilag(ifit));
  cmb(ifit)=cmf(ifit)/6+5/6*cm(ifit);
  [m,im]=min(cmb);
 else
  m=m1;
  im=im1;
 end

% check correlation of best lag 
 c3=[c1(igood,1),c2(igood,im)];

 % remove bad spikes that degrade correlation
 nspike=0;
 for i=1:4
  c3d=diff(c3');
  ii=find(abs(c3d)>(2*(1+i/8)*std(c3d)));
  if (length(c3)-length(ii))<(0.8*length(igood)), break, end
  nspike=nspike+length(ii);
  c3(ii,:)=[];
 end
 if nspike>10
   disp([' bestlag removed ',int2str(nspike),' spikes'])
 end
 
 % compute correlation
 co=corrcoef(c3);
 co=co(1,2);
 lag=ilag(im);
 
 %clf
 %plot(ilag,cm,ilag,cmf,'-g',ilag,cmb,'-r',ilag(im1),m1,'bp',ilag(im),m,'rp')
 %pause

end


% prepare joint index vector
i1=1:l1;
i2=i1+lag;
ii=find(i2>0 & i2<=l2);
i1=i1(ii); 
i2=i2(ii);