C Takes 3 passes through the data....
C Changed to 5 passes 7/4/05
C Sometimes it starts off with a bum first value though...
	subroutine best3(xrefs,yrefs,arefs,nf,
     &                   xnews,ynews,anews,nn,
     &       avexshift,aveyshift,aveashift,std,nk)
C maximum shift in x and y
	parameter (maxshift=4000)
	parameter (maxmags=200)
	parameter (maxhits=10000)
	dimension xrefs(*),yrefs(*),arefs(*)
	dimension xnews(*),ynews(*),anews(*)
	dimension ixshift(maxshift),iyshift(maxshift),iashift(maxmags)
	dimension xtol(5),ytol(5),atol(5)
	dimension xshift(maxhits),yshift(maxhits),ashift(maxhits),
     &   index(maxhits)
	data xtol/4000.,500.,100.,50.,10/
	data ytol/4000.,500.,100.,50.,10/
	data atol/2.,1.,1.,0.5,0.5/
	ipass=0
	avexshift=0.
	aveyshift=0.
	aveashift=0.
C shifts are defined as follows:
C	xref=xnew+xshift
C	yref=ynew+yshift
C	aref=anew+ashift
C  so, shift=ref-new

C First, zero the integer bins:
C
900	continue
C First, zero the bins:
	do 18 k=1,maxshift
	ixshift(k)=0
	iyshift(k)=0
18	continue
	do 19 k=1,maxmags
	iashift(k)=0
19	continue
C
C
C New pass.... 
	nk=0
	z2=0
	ipass=ipass+1
	do 20 i=1,nf
	do 20 j=1,nn
	dmag=arefs(i)-anews(j)
	if(abs(dmag-aveashift).gt.atol(ipass)) go to 20
	dx=xrefs(i)-xnews(j)
	if(abs(dx-avexshift).gt.xtol(ipass)) go to 20
	dy=yrefs(i)-ynews(j)
	if(abs(dy-aveyshift).gt.ytol(ipass)) go to 20
C Now turn things into integers for determining the mode:
	idx=dx+0.5+maxshift/2.
	idy=dy+0.5+maxshift/2.
	iam=dmag*10.+maxmags/2.+0.5
	if(idx.lt.1) go to 20
	if(idx.gt.maxshift) go to 20
	if(idy.lt.1) go to 20
	if(idy.gt.maxshift) go to 20
	if(iam.lt.1) go to 20
	if(iam.gt.maxmags) go to 20
C OK, it passed all these tests.  Now fill up the hoppers:
	nk=nk+1
	if(nk.gt.maxhits) stop 'too many hits...'
	iashift(iam)=iashift(iam)+1
	ixshift(idx)=ixshift(idx)+1
	iyshift(idy)=iyshift(idy)+1
	z2=z2+(dx-avexshift)**2+(dy-aveyshift)**2
	xshift(nk)=dx
	yshift(nk)=dy
	ashift(nk)=dmag
20	continue
	std=sqrt(z2/(nk-1.))
C OK, now how well did we do?
	maxy=0
	maxx=0
	maxa=0
	bestx=9999.
	besty=9999.
	bestmag=9999.
	do 21 i=2,maxshift-1
C Smooth by 3:
	valx=ixshift(i-1)+ixshift(i)+ixshift(i+1)
	if(valx.gt.maxx) then
	maxx=valx
	bestx=((i-1)*ixshift(i-1)+i*ixshift(i)+
     &                               (i+1)*ixshift(i+1))/valx
	end if
	valy=iyshift(i-1)+iyshift(i)+iyshift(i+1)
	if(valy.gt.maxy) then
	maxy=valy
	besty=((i-1)*iyshift(i-1)+i*iyshift(i)+
     &                               (i+1)*iyshift(i+1))/valy
	end if
	write(15,*) i,valx,valy
21	continue
	do 22 i=2,maxmags-1
	valmag=iashift(i-1)+iashift(i)+iashift(i+1)
	if(valmag.gt.maxa) then
	maxa=valmag
	besta=((i-1)*iashift(i-1)+i*iashift(i)+
     &                               (i+1)*iashift(i+1))/valmag
	end if
22	continue
C OK, what are the averageshifts?
	avexshift=bestx-maxshift/2.
	aveyshift=besty-maxshift/2.
	aveashift=(besta-maxmags/2.)/10.
Cwrite(6,134) ipass,avexshift,aveyshift,aveashift,std,nk
134	format(i6,4f10.3,i6)
	if (ipass.lt.5) go to 900
	if(nk.gt.5) then
	call sort(nk,xshift,index)
Cwrite(6,*) avexshift,xshift(index(nk/2))
	avexshift=xshift(index(nk/2))
	call sort(nk,yshift,index)
Cwrite(6,*) aveyshift,yshift(index(nk/2))
	aveyshift=yshift(index(nk/2))
	call sort(nk,ashift,index)
Cwrite(6,*) aveashift,ashift(index(nk/2))
	aveashift=ashift(index(nk/2))
	else
	std=1000.
	end if
	return
	end