C Written Phil Massey, Lowell Obs. Nov 13, 2001 
C Takes "kfiles" worth of stellar photometry, and determines the pixel
C shifts among stars found in common.
	parameter(kfiles=12)
	parameter(maxstars=200000)
CC
CC
CC
	dimension x(maxstars,kfiles),y(maxstars,kfiles),
     &           dx(maxstars,kfiles),dy(maxstars,kfiles),
     &           number(kfiles)
	dimension inum(8,8),ax(8,8),ay(8,8),xbin(8,8,maxstars),
     &    ybin(8,8,maxstars)
	dimension xval(maxstars),yval(maxstars)

	character*25 newfile
	character*25 namelist
C
	open(unit=1,file='temp',status='old')
	read(1,2) namelist
	close(unit=1)
	open(unit=1,file=namelist,status='old')
	nfiles=0
C read everything into memory.
1	read(1,2,end=20) newfile
Cwrite(6,2) newfile
2	format(a25)
	alast=lnblnk(newfile)
	nfiles=nfiles+1
	if(nfiles.gt.kfiles) stop 'too many files'
	open(unit=2,file=newfile(1:alast)//'comxy',status='old')
	istar=0
9	read(2,*,end=15)   xk,yk
     	istar=istar+1
	x(istar,nfiles)=xk
	y(istar,nfiles)=yk
        go to 9                                  
15	number(nfiles)=istar
	close(unit=2)
	write(6,16) nfiles,istar
16	format('File=',i3,' Number=',i5)
	if(nfiles.lt.kfiles) then
	go to 1
	else
	write(6,17) kfiles
17	format('Number of images exceeds',i5)
	stop
	end if
CCCCCCCCCCCCCCCCCCCCCCCCC
20	close(unit=1)
Cwrite(6,21) nfiles,(number(j),j=1,nfiles)
21	format(i5,30i5)
	do 22 j=2,nfiles
	if(number(j).ne.number(1)) stop 'not all the same number of stars'
22	continue
	num=number(1)
C Now go through and see what the deltas are:
	do 100 j=1,nfiles
	do 100 i=1,num
	dx(i,j)=x(i,j)-x(i,1)
	dy(i,j)=y(i,j)-y(i,1)
Cwrite(6,*) i,j,x(i,j),x(i,1),dx(i,j),dy(i,j)
100	continue
C Now do some statistics.  Divide things up into 1000x1000 sections. There
C will be 64 of them.
	do 301 j=1,nfiles
	write(6,300) j
300	format('FILE=',i5)
C Now for each section first zero the counter
	do 201 k=1,8
	do 201 l=1,8
	inum(k,l)=0
201	continue
	do 250  i=1,num
	ix=x(i,j)/1000+1
	iy=y(i,j)/1000+1
	if(ix.lt.1) ix=1
	if(ix.gt.8) ix=8
	if(iy.lt.1) iy=1
	if(iy.gt.8) iy=8
	inum(ix,iy)=inum(ix,iy)+1
	xbin(ix,iy,inum(ix,iy))=dx(i,j)
	ybin(ix,iy,inum(ix,iy))=dy(i,j)
250	continue
	do 350 k=1,8
	do 350 l=1,8
Cwrite(6,650) k,l, inum(k,l)
650	format('Number for',2i5,'=',i5)
	if(inum(k,l).lt.9) then
	ax(k,l)=-99.
	ay(k,l)=-99.
	else
	do 251 i=1,inum(k,l)
	xval(i)=xbin(k,l,i)
Cwrite(6,254) xval(i)
254	format('Loading xval',f8.3)
	yval(i)=ybin(k,l,i)
251	continue
	call amedian(xval,inum(k,l),amed)
Cwrite(6,252) xval(3),inum(k,l),amed
252	format('Here boy!',f8.3,i5,f8.3)
	ax(k,l)=amed
	call amedian(yval,inum(k,l),amed)
	ay(k,l)=amed
	lfirstx=(k-1)*1000+1
	lastx=lfirstx+998
	lfirsty=(l-1)*1000+1
	lasty=lfirsty+998
	write(6,203) lfirstx,lastx,lfirsty,lasty, inum(k,l),ax(k,l),
     &   ay(k,l)
203	format("X:",i4.4,"-",i4.4," Y:",i4.4,"-",i4.4," N:",i4," Med:",
     &   2f8.2)
	end if
350	continue
301     continue
	stop
	end

	subroutine xfit(x,sigmax,npts,mode,xmean,sigmam,sigma)
	double precision sum,sumx,weight,free
	dimension x(*),sigmax(*)
11	sum=0.
	sumx=0.
	sigma=0.
	sigmam=0.
	if(npts.eq.0) return
	if(npts.eq.1) then
	xmean=x(1)
	sigma=sigmax(1)
	sigmam=sigmax(1)
	return
	end if
20	do 32 i=1,npts
21	if(mode) 22,22,24
22	weight=1.
	go to 31
24	if(sigmax(i).lt.0.007) sigmax(i)=0.007
	weight=1./sigmax(i)**2
C small little cheat removed here, damn it!
31	sum=sum+weight
32	sumx=sumx+weight*x(i)
c find mean and s.d.
41	xmean=sumx/sum
51	do 52 i=1,npts
52	sigma=sigma+(x(i)-xmean)**2
	free=npts-1
54	sigma=dsqrt(sigma/free)
61	if(mode) 62,64,66
62	sigmam=dsqrt(xmean/sum)
	go to 70
64	sigmam=sigma/dsqrt(sum)
	go to 70
66	sigmam=dsqrt(1./sum)
70	return
	end
	subroutine amedian(data,n,amed)
	dimension data(*),indx(200000)
	call sort(n,data,indx)
	if( ((n/2)*2).eq.n) then
	amed=(data(indx(n/2))+data(indx(n/2+1)))/2.
	else
	amed=data(indx(n/2+1))
	end if
	return
	end
        subroutine sort(n,arrin,indx)
        dimension arrin(*),indx(*)
        do 11 j=1,n
        indx(j)=j
11      continue
        if(n.eq.1) return
        l=n/2+1
        ir=n
10      continue
                if(l.gt.1) then
                        l=l-1
                        indxt=indx(l)
                        q=arrin(indxt)
                ELSE
                        indxt=indx(ir)
                        q=arrin(indxt)
                        indx(ir)=indx(1)
                        ir=ir-1
                        if(ir.eq.1) then
                                indx(1)=indxt
                                return
                        endif
                endif
        i=l
        j=l+l
20      if(j.le.ir) then
                if(j.lt.ir) then
                        if(arrin(indx(j)).lt.arrin(indx(j+1)))j=j+1
                endif
                if(q.lt.arrin(indx(j)))then
                        indx(i)=indx(j)
                        i=j
                        j=j+j
                else
                        j=ir+1
                endif
        go to 20
        endif
        indx(i)=indxt
        go to 10
        end

