C Written, more or less, by Erin Darnell
C Added a new criterion 5/9/06
	character*80 infiles, inline
        character*18 name,namet
        parameter(maxobs=1000)
	integer alast, fileNo, lineCount, i,k,ik
	double precision hjd
	real magmean, errmean, sd, interr
	real eoveri, momaoveri, mimaoveri
	real minv, maxv, mode, fudge, ratiog
	real v(maxobs), verr(maxobs),vn(maxobs)
	fileNo = 0
	fudge=0.020
C starts here
	open(unit=3,file='mcheck3.out',status='new')
	open(unit=1,file='allstars',status='old')
	write(3, 16) 
16      format ('#file,  meanV,   interr,  sd,  e/i,  moma/i,',
     &  '  mima/i,  ratiog')
1	read(1,2,end=999) infiles
	fileNo = fileNo + 1
2	format(a80)
	alast=lnblnk(infiles)
	open(unit=2,file=infiles(1:alast),status='old')
	lineCount = 0
3	read(2,2,end=88) inline
        if(inline(1:1).eq.'#') then
           if(inline(2:2).eq.'F') go to 3
           name=inline(54:71)
	   read(inline(33:38),*) bmv
           go to 3
           END IF
        read(inline,*) hjd,vtemp,verrtemp
        if(vtemp.gt.90.) go to 3 
	lineCount = lineCount + 1
	if(lineCount.gt.maxobs) stop 'too many stars'
        read(inline,*) hjd, v(lineCount), verr(lineCount)
	go to 3

88	continue
	call clean(v,lineCount,vn,nvn,sign)
C "sign" is now the sigma excluding outlyers
	magmean = calcmagmean(v, lineCount) / lineCount
	errmean = calcerrmean(verr, lineCount) / lineCount
	maxv = calcmaxv(v, lineCount)
	minv = calcminv(v, lineCount)
	mode = calcmode(v, lineCount)
	sum2=0.
	do 99 i = 1, lineCount
	sum2=sum2+(magmean-v(i))**2
99	sd = sqrt(sum2/(lineCount-1.))
	interr=sqrt(errmean**2+fudge**2)
        eoveri=sd/interr
	eovere=sd/sign
	mimaoveri=(maxv-minv)/interr
	momaoveri=(mode-magmean)/interr	
	ratiog=calcg3s(v, verr, magmean, fudge, lineCount)
        ik=0
        namet=''
        do k=1,18
        if(name(k:k).ne.' ')  then
        ik=ik+1
        namet(ik:ik)=name(k:k)
        end if
        end do
        if(ik.eq.0) namet='nocrossid'
	write(6, 17) infiles(1:alast),magmean,interr,sd,eoveri,momaoveri,mimaoveri,ratiog,namet,bmv
	write(3, 17) infiles(1:alast),magmean,interr,sd,eoveri,momaoveri,mimaoveri,ratiog,namet,bmv
17      format(a,f6.2,f5.2,f5.2,f7.1,2f7.2,f7.1,1x,a,f6.2)


	close(unit=2)
	goto 1
999	print *, fileNo, ' stars analyzed'
	close(unit=3)
	stop
	end

	real function calcmagmean(v, lineCount)
	integer lc, lineCount
	real v(*)
	calcmagmean = 0.0
	do 20 lc = 1, lineCount
		calcmagmean = calcmagmean + v(lc)
20	calcmagmean = calcmagmean
	return
	end

	real function calcerrmean(verr, lineCount)
	integer lc, lineCount
	real verr(*)
	calcerrmean = 0.0
	do 30 lc = 1, lineCount
		calcerrmean = calcerrmean + verr(lc)
30	calcerrmean = calcerrmean
	return
	end

	real function calcsdmean(sd, lineCount)
	integer lc, lineCount
	real sd(*)
	calcsdmean = 0.0
	do 40 lc = 1, lineCount
		calcsdmean = calcsdmean + sd(lc)
40	calcsdmean = calcsdmean
	return
	end

	real function calcmaxv(v, lineCount)
	integer lc, lineCount
	real v(*)
	calcmaxv = 0.0
	do 50 lc = 1, lineCount
		if (v(lc) .GT. calcmaxv) then 
			calcmaxv = v(lc)
		endif
50	calcmaxv = calcmaxv
	return
	end


	real function calcminv(v, lineCount)
	integer lc, lineCount
	real v(*)
	calcminv = 100
	do 60 lc = 1, lineCount
	if (v(lc) .LT. calcminv) calcminv=v(lc)
60	calcminv = calcminv
	return
	end


	real function calcmode(v, lineCount)
	integer lc, lineCount, imhist(400), ix, i, ibig
	real v(*)
	ix=0
	ibig=0.
        do 50 i=1,300
        imhist(i)=0
 50	continue
	do 70 lc = 1, lineCount
	ix=((v(lc)-5.)/.05)
	if (ix .GT. 300) go to 70
	if (ix .LT. 1) go to 70
	imhist(ix)=imhist(ix)+ 1
	if (imhist(ix) .GT. ibig) then
        calcmode= ix
        ibig=imhist(ix)
        end if
70	continue
        calcmode=calcmode/20.+5.+0.025
	return
	end


	real function calcg3s(v, verr, magmean, fudge, lineCount)
	integer lc, lineCount, n
	real v(*), verr(*)
	real expR, magmean, fudge, err, ratio, nexpfound
	n=0
        do 80 lc = 1, lineCount
		expR=v(lc)-magmean
		err=sqrt(verr(lc)**2+fudge**2)
	if (expR .GT. 3*err) n=n+1
 80	CONTINUE
	nexpfound= .0025*lineCount
	ratio=float(n)/nexpfound
	calcg3s = ratio
	return
	end
C Here's the good stuff
        subroutine clean(x,n,xnew,nnew,sign)
        dimension x(*),xnew(*)
        sig=2.5
C Fill up the new array with everything:
        do 2 i=1,n
        xnew(i)=x(i)
2       continue
        nnew=n
        do 10 j=1,10
C First, find the median:
        call median(xnew,nnew,amed)
C Next find the average and sigma around the median:
        call findsigma(xnew,nnew,amed,sigma)
Cwrite(6,*) amed,sigma,nnew,n
C Reload array:
        nnew=0
C Now eliminate points sig points away
        do 1 i=1,n
        dif=x(i)-amed
Cwrite(6,*) i,x(i),amed,dif,sig*sigma
        if(abs(dif).gt.sig*sigma) go to 1
        nnew=nnew+1
        xnew(nnew)=x(i)
1       continue
10      continue
        sign=sigma
        return
        end
C
        subroutine findsigma(x,n,amed,sigma)
        dimension x(*)
        sum=0.0
        do 1 i=1,n
        sum=sum+(x(i)-amed)**2
1       continue
        sigma=sqrt(sum/(n-1.))
        return
        end
c subroutine xfit.f
c
c source
c   Bevington, page 76.
c
c purpose
c   calculate the mean and estimated errors for a set of data points
c
c usage
c   call xfit (x, sigmax, npts, mode, xmean, sigmam, sigma)
c
c description of parameters
c   x      - array of data points
c   sigmax - array of standard deviations for data points
c   npts   - number of data points
c   mode   - determines method of weighting
c            +1 (instrumental) weight(i) = 1./sigmax(i)**2
c             0 (no weighting) weight(i) = 1.
c            -1 (statistical)  weight(i) = 1.
c   xmean  - weighted mean
c   sigmam - standard deviation of mean
c   sigma  - standard deviation of data
c
c subroutines and function subprograms required
c   none
c
	subroutine xfit (x,sigmax,npts,mode,xmean,sigmam,sigma)
	double precision sum,sumx,weight,free
	dimension x(1),sigmax(1)
c
c accumulate weighted sums
c
11	sum=0.
	sumx=0.
	sigma=0.
	sigmam=0.
20	do 32 i=1,npts
21	if (mode) 22,22,24
22	weight=1.
	goto 31
24	weight=1./sigmax(i)**2
31	sum=sum+weight
32	sumx=sumx+weight*x(i)
c
c evaluate mean and standard deviations
c
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)
	goto 70
64	sigmam=sigma/dsqrt(sum)
	goto 70
66	sigmam=dsqrt(1./sum)
70	return
	end
	subroutine median(data,n,amed)
	dimension data(*),indx(300000)
	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