#include    "focas1.h"


/* FINDSTARS -- Find stars based on the first radial moments. */

findstars (input, output, nstars, sigmamin, nstart, lsigma, hsigma)
char	*input;
char	*output;
int	*nstars;
float	*sigmamin;
int	*nstart;
float	*lsigma;
float	*hsigma;
{
    struct objrec ob, *obs;
    int		i, j, n, nmax, nlast;
    double	ir1, prob, mean, sigma, ratio, resid, residmax, sum1, sum2;	
    char	*malloc();

    /* Get the objects with the largest "prob" values */
    cfd = catopen (input, 0);

    nmax = *nstart;
    obs = (struct objrec *) malloc (nmax * sizeof (struct objrec));
    sum1 = 0.;
    sum2 = 0.;
    n = 0;
    i = 0;
    for (;;) {
	if (rdcatob (cfd, 0L, &ob))
	    break;
	if (ffilter (&ob) != 1)
	    continue;
	if ((ob.eflgs & EVAL) == 0)
	    continue;
	if ((ob.eflgs & (EDGE | PEAK | DARK)) != 0)
	    continue;

	prob = ob.prob;
	ir1 = ob.ir1;
	if (n < nmax) {
	    obs[n] = ob;
	    if (prob < obs[i].prob)
		i = n;
	    sum1 += ir1;
	    sum2 += ir1 * ir1;
	    n++;
	} else {
	    if (prob < obs[i].prob)
		continue;
	    sum1 += ir1 - obs[i].ir1;
	    sum2 += ir1 * ir1 - obs[i].ir1 * obs[i].ir1;
	    obs[i] = ob;
	    ir1 = obs[0].prob;
	    for (i=0, n=1; n<nmax; n++) {
		prob = obs[n].prob;
		if (prob < ir1) {
		    ir1 = prob;
		    i = n;
		}
	    }
	}
    }
    fclose (cfd);

    if (n > 1) {
	mean = sum1 / n;
	sigma = (n * sum2 - sum1 * sum1) / (n * (n - 1));
	sigma = sigma > 0. ? sigma : 0.;
	sigma = sqrt (sigma);
	ratio = sigma / mean;
    }

    nmax = n;
    nlast = 0;
    while (n > *nstars && ratio > *sigmamin && n != nlast) {
	nlast = n;
	j = -1;
	residmax = 0.;
	for (i = 0; i < nmax; i++) {
	    if (obs[i].ir1 == 0.)
		continue;
	    resid = (obs[i].ir1 - mean) / sigma;
	    if (resid > *hsigma || resid < -(*lsigma)) {
		if (abs (resid) > residmax) {
		    j = i;
		    residmax = abs (resid);
		}
	    }
	}
	if (j > 0) {
	    ir1 = obs[j].ir1;
	    sum1 -= ir1;
	    sum2 -= ir1 * ir1;
	    n--;
	    if (n > 1) {
		mean = sum1 / n;
		sigma = (n * sum2 - sum1 * sum1) / (n * (n - 1));
		sigma = sigma > 0. ? sigma : 0.;
		sigma = sqrt (sigma);
		ratio = sigma / mean;
	    }
	    obs[j].ir1 = 0.;
	}
    }

    if (n == 0)
	focaserr (1, "findstars:  No objects found - ", input);

    fprintf (stderr, "Final statistics:\n");
    fprintf (stderr,
       "Objects = %d, <ir1> = %.6g, sigma = %.6g, sigma/<ir1> = %.6g\n",
       n, mean, sigma, ratio); 

    if (strcmp (output, "STDOUT") == 0) {
	cfd = stdout;
	wtcathdr (cfd, 0);
    } else
	cfd = catopen (output, 1);

    for (i = 0; i < nmax; i++) {
	if (obs[i].ir1 == 0.)
	    continue;
	wtcat (&obs[i], cfd);
    }
    fclose (cfd);

    free (obs);
}