--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bestlag.m	Tue Oct 20 16:23:49 2009 -0400
@@ -0,0 +1,108 @@
+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); 
+