C Determines the MODAL fwhm produced by fitpsf
C Modified April 21/22 to recompute fwhm if mean is significantly
C higher than the mode (--->CRs are dominating)
CModified to allow for bad seeing SMARTS images...
C Changed to MEDIAN! 8/26/2005
	program radmode
	implicit double precision(a-h,o-z)
	parameter(maxstars=300000)
	dimension data(maxstars),sig1(maxstars)
	open(unit=1,file='allrad',status='old')
	open(unit=2,file='tempfwhm',status='new')
	ix=0
C start off not counting anything less than 2.3 pixels fwhm
	alimit=1.0
	iflag=0
1	read(1,*,end=99) rad
C exclude silly values:
	if((rad.ge.alimit).and.(rad.lt.6.5)) then
	ix=ix+1
	if(ix.ge.maxstars) go to 99
C factor of 2.355 is to go from sigma to fwhm
	data(ix)=rad*2.355
	sig1(ix)=1.0
	end if
	go to 1
99	continue
Cwrite(6,921) ix,data(1)
921	format('About to call median with',i5,1x,f9.2)
	call median(data,ix,amed)
	call modal(data,ix,amode)
	call xfit(data,sig1,ix,0,xmean,sigmam,sigma2)
	write(6,1101) ix,amed,amode,xmean,iflag,alimit
1101	format('Radbig: ix,amed,amode,mean,iflag,alimit=',i10,3f6.2,i2,
     &   f6.2)
	if(amode.eq.-99.) stop "HELP!"
Cif((xmean.gt.(amode+0.25)).and.(iflag.eq.0)) then
C I think there's a problem!
Ciflag=1
Calimit=xmean/2.355
Cix=0
Crewind(1)
Cgo to 1
Cend if
C Changed to median!
	write(2,101) amed
101	format(f8.2)
	close(unit=2)
	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(300000)
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)
	parameter(maxstars=300000)
	dimension ibin(200)
	dimension data(maxstars)
	do 1 i=1,200
	ibin(i)=0
1	continue
	if(ix.gt.maxstars) then
	amode=-99.
	return
	end if
	do 100 i=1,ix
	val=data(i)
	ival=10.*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=20,199
	write(4,*) i/10.,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/10.
	return
	end