	program mode3pick
C written by Phil Massey, Lowell Observatory, March 2002
C Selects psf stars
C revised June 2003, made things a little looser (.02-->.05) 
C since we are now using the fwhm as a filter
C CHANGED TO THE MEDIAN 6/16/03
C MADE MORE FLEXIBLE 6/16/03
C Added support for the second aperture 6/16/03
	implicit double precision(a-h,o-z)
	dimension data1(1000),data2(1000)
C this file contains all the ap=3.5,5,15 stuff:
	open(unit=1,file='tempmode',status='old')
C this file contains all the uncrowded stars
	open(unit=9,file='pospst',status='old')
	open(unit=2,file='temppst3',status='new')
	open(unit=10,file='temppst4',status='new')
	write(2,2)
2	format('#N ID    XCENTER   YCENTER   MAG         '/
     &         '#U ##    pixels    pixels    magnitudes'/
     &         '#F %9d   %10.3f    %10.3f    %12.3f  '/
     &         '#')
	ix=0
1	read(1,*,end=99,err=1) id,xc,yc,amag1,amag2,amag3,err1,err2,err3
	if(err1.gt.0.01) go to 1
	if(ix.eq.1000) go to 99
	ix=ix+1
	data1(ix)=amag1-amag2
	data2(ix)=amag1-amag3
	write(3,3) amag1-amag2,amag1-amag3,xc,yc,amag1,amag2,amag3,
     &  err1,err2,err3
3	format(10f10.3)
	go to 1
99	continue
	call median(data1,ix,amed1)
	call median(data2,ix,amed2)
	call modal(data1,ix,amode1)
	call modal(data2,ix,amode2)
	write(10,1011) amed1,amed2
	close(unit=10)
1011	format(2f10.3)
	write(6,100) ix,amed1,amode1,amed2,amode2
100	format('For',i5,' stars the med, mode aper cor was',4f9.3)
C OK, now find at least 7 stars.
C First, be conservative:
	errlim=0.01
	close=0.03
	fact=1.
1010    rewind(1)
	ik=0
11	read(1,*,end=999) id,xc,yc,amag1,amag2,amag3,err1,err2,err3
	if(err1.gt.fact*errlim) go to 11
	dif1=amag1-amag2
	dif2=amag1-amag3
C made this more generous
	if(abs(dif1-amed1).gt.fact*close) go to 11
	if(abs(dif2-amed2).gt.fact*close) go to 11
C OK, we've made it this far.  Now lets find its real id:
	rewind(9)
101	read(9,*,end=11) idl,xcl,ycl
	dist2=(xcl-xc)**2+(ycl-yc)**2
	if(dist2.lt.1.) go to 102
	go to 101
102	ik=ik+1
C Assumes that the stars have been sorted by mag
	if(ik.gt.100) stop
C Note that we are using the id system for the original mag file
	write(2,12) idl,xcl,ycl,amag1
12	format(i9,2f10.3,f12.3)
	go to 11
999	continue
	if(ik.lt.7) then
	fact=fact+0.75
	if(fact.gt.4.0) then
	write(6,1014) 
1014	format('fact is too large!  I quit!')
	fact=fact-0.75
	go to 10155
	end if
	rewind(2)
	write(2,2)
	go to 1010
	end if
	close(unit=2)
10155	write(6,1015) ik,fact,fact*errlim,fact*close
1015	format('Settled on',i3,' stars with factor',f5.1,' errlim,close=',
     &  2f8.3)
	stop
	end
	subroutine xfit(x,sigmax,npts,mode,xmean,sigmam,sigma)
	implicit double precision (a-h,o-z)
	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 median(data,n,amed)
	implicit double precision(a-h,o-z)
	dimension data(*),indx(1000)
Cwrite(6,71) n,data(5)
71	format('Calling sort',i5,f8.3)
	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)
	implicit double precision(a-h,o-z)
        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
	subroutine modal(data,ix,amode)
	implicit double precision (a-h, o-z)
	dimension ibin(200)
	dimension data(1000)
	do 1 i=1,200
	ibin(i)=0
1	continue
	if(ix.gt.1000) then
	amode=-99.
	return
	end if
	do 100 i=1,ix
	val=data(i)
	ival=100.*val+.5
	if(ival.lt.1) ival=1
	if(ival.gt.200) ival=200
	ibin(ival)=ibin(ival)+1
100	continue
	ibig=0
	ibest=0
	do 200 i=4,197
	write(4,*) i/100.,ibin(i)
	if(ibin(i).gt.ibig) then
	ibig=ibin(i)
	ibest=i
	end if
200	continue
	try2=(ibest*(ibin(ibest))+(ibest-1)*ibin(ibest-1)+
     &    (ibest+1)*ibin(ibest+1))/float((ibin(ibest)+ibin(ibest-1)+
     &    ibin(ibest+1)))
	amode=try2/100.
	return
	end
