/* go through all the files on the command line (which should be
	   standard-star .pht files) and read in information on airmass,
	   time, filter, etc.  then, from the part of the file after
	   the header, figure out which star was the standard and get its
	   instrumental magnitudes.

	   after that, somehow go through all the data for each star and 	
	   figure out how well all the instrumental mags match each other.
	   perhaps figure out extinction (?), but certainly give some indication
	   of fluctuations in the transparency, and the scatter among all
	   measurements. 

	   modified 3/30/92 to ignore files with a REQID that doesn't begin
	   with "st_"  MWR 

	   modified 5/30/92 to merely set WINDSPEED and HUMIDITY to -1 
	   if we can't get their values from the header, rather than
	   terminating execution. MWR */

#include <stdio.h>
#include <math.h>
#include "pf.h"

#define MAX_STAR  5000		/* max number of stars (= #.pht + .coo files) */
#define STRLEN    40		/* max length of star names and filters */
#define MAX_APER  5			/* max number of apertures a star may have */
#define LINELEN   256		/* max length of data lines */
#define NODATA    99.0		/* value assigned to a variable w/ no data */
#define PHT_POS   38		/* position in a .pht data line at which */
                    		/*    mag/err pairs begin */
#define PHT_LEN   15		/* length of a mag/err pair in chars */
#define AVG_TIME  0.10		/* fraction of an hour within which observations */
                      		/*    may be considered "simultaneous" */
#define PHT       1			/* denotes a .pht star */
#define COO       2			/* denotes a .coo star */
#define NOT_STANDARD 1		/* this star is not a standard star */

	/* flags for functionality (OR'ed together) */
#define PRINT     1			/* print out the data */
#define SCATTER   2			/* do scatter analysis within groups; DEFAULT */
#define EXTINCT   4			/* do extinction analysis */

	/* various flags for sorting */
#define BY_TIME   1			/* sort observations by time */
#define BY_NAME   2			/* sort by star primarily, by filter secondarily */
                   			/*    and by time tertially */

static struct u_coo {
	int peak;
	double fwhm;
	double round;
	double sharp;
};

static struct u_pht {
	int n_aper;						/* number of aperture mags measured */
	double inst_mag[MAX_APER];		/* the observed instrumental magnitude */
	double inst_err[MAX_APER];		/* ditto uncertainties */
};

static struct s_star {
	char name[STRLEN];
	char filter[STRLEN];
	double exptime;
	double airmass;
	double windspeed;
	double humidity;
	double teltemp;
	double telfocus;
	double real_mag;		/* the mag of this star in standard system */
	double ut;
	int type;				/* either PHT or COO */
	union {	
		struct u_coo c;
		struct u_pht p;
	} dat;
};

static struct s_coo {
	int id;
	double xc;
	double yc;
	double peak;
	double fwhm;
	double round;
	double sharp;
};

static struct s_pht {
	int id;
	double xc;
	double yc;
	double sky;
	double skysig;
	double mag[MAX_APER];
	double err[MAX_APER];
};


	/* definitions of static functions in this file */
static int do_pht( /* char *file */ );
static int do_coo( /* char *file */ );
static int slurp_header( /* char *file, struct s_star *star */ );
static int slurp_pht_data( /* char *file, struct s_star *star */ );
static int slurp_coo_data( /* char *file, struct s_star *star */ );
static int read_coo_line( /* char *line, s_coo *coo */ );
static int read_pht_line( /* char *line, s_pht *pht, int *naper */ );
static void print_star( /* struct s_star *s */ );
static void proc_pht( /* int flags */ );
static void sortem( /* int flag */ );
static void do_print( /* */ );
static void do_scatter( /* */ );
static void do_extinct( /* */ );
static int comp_name( /* struct s_star *s1, *s2 */ );
static int comp_time( /* struct s_star *s1, *s2 */ );
static void replace_spaces( /* char *str */ );


char progname[] = "quality";
static char usage[] =    "quality [debug] [print] [scatter] [extinct] *.pht *.coo ";
static int num_stars = 0;
static int debug_flag = 0;				/* if set to 1, print stuff a lot */
static struct s_star stars[MAX_STAR];


	/* usage: quality [help] [debug] [print] [scatter] [extinct] *.pht *.coo */ 

main(argc, argv)
int argc;
char *argv[];
{
	int i, flags;
	char file[100];
	struct s_star *s;

	if (argc == 1) {
		fprintf(stderr, "usage: %s\n", usage);
		exit(0);
	}
	
	/* read in data from all the files */
	for (i = 1; i < argc; i++) {
		if (strcmp(argv[i], "help") == 0) {
			fprintf(stderr, "usage: %s\n", usage);
			exit(0);
		}
		if (strcmp(argv[i], "debug") == 0) {
			debug_flag = 1;
			continue;
		}
		if (strcmp(argv[i], "print") == 0) {
			flags |= PRINT;
			continue;
		}
		if (strcmp(argv[i], "scatter") == 0) {
			flags |= SCATTER;
			continue;
		}
		if (strcmp(argv[i], "extinct") == 0) {
			flags |= EXTINCT;
			continue;
		}
		strcpy(file, argv[i]);
		if (strncmp(file + (strlen(file) - 3), "pht", 3) == 0) {
			do_pht(file);
		}
		else if (strncmp(file + (strlen(file) - 3), "coo", 3) == 0) {
			do_coo(file);
		}
		else {
			fprintf(stderr, "file %s neither .coo nor .pht \n", file);
			exit(1);
		}
	}

	if (flags == 0)
		flags = SCATTER;

	/* now that we've read in all the data, time to process it */
	proc_pht(flags);


	exit(0);
}

	/* do all the necessary things for a .pht file, adding the new star
	   into the array if possible. return 0 if all went well, or -1 if
	   there's a problem. */

static int
do_pht(file)
char *file;
{
	struct s_star *s;
	int ret;

	s = &(stars[num_stars]);
	s->type = PHT;
	if ((ret = slurp_header(file, s)) == -1) {
		fprintf(stderr, "%s.do_pht: slurp_header returns -1 for file %s\n",
					progname, file);
		return(-1);
	}
	else if (ret == NOT_STANDARD) {
		return(0);
	}
	if (slurp_pht_data(file, s) != 0) {
		fprintf(stderr, "%s.do_pht: slurp_pht_data returns non-zero for file %s\n",
					progname, file);
		return(-1);
	}
	num_stars++;

	return(0);
}


	/* do all the necessary things for a .coo file */

static int
do_coo(file)
char *file;
{
	struct s_star *s;
	int ret;


	s = &(stars[num_stars]);
	s->type = COO;
	if ((ret = slurp_header(file, s)) == -1) {
		fprintf(stderr, "%s.do_coo: slurp_header returns -1 for file %s\n",
					progname, file);
		return(-1);
	}
	else if (ret == NOT_STANDARD) {
		return(0);
	}
	if (slurp_coo_data(file, s) != 0) {
		fprintf(stderr, "%s.do_coo: slurp_coo_data returns non-zero for file %s\n",
					progname, file);
		return(-1);
	}
	num_stars++;

	return(0);
}


	/* read in the appropriate quantities from a pseudo-FITS header.
	   note that we copy only the portion of the star's name before
	   a parenthesis, thus converting the name
	                     SA 101 311 (V)
	   to
	                     SA 101 311
	   we read in filter information later, so we don't need it in the
	   name.

	   return 0 if all went well, NOT_STANDARD if the file's REQID does
	   not begin with "st_", or -1 otherwise. */

static int
slurp_header(file, star)
char *file;
struct s_star *star;
{
	int i;
	char value[PF_LEN + 1], str[PF_LEN + 1];
	pf_handle pfp;

	if ((pfp = pf_open(file, 1, PF_EXIST)) != 0) {
		fprintf(stderr, "%s.slurp_header: can't pf_open file %s\n", 	
					progname, file);
		return(-1);
	}
	if (pf_getval(pfp, "REQID", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read REQID from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_valtostr(value, str) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtostr REQID ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
#if 0
	if (strncmp(str, "st_", 3) != 0)  {
		pf_close(pfp);
		return(NOT_STANDARD);
	}
#endif

	if (pf_getval(pfp, "OBJECT", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read OBJECT from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_valtostr(value, str) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtostr OBJECT ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
#if 0
	for (i = 0; ((str[i] != '\0') && (str[i] != '(') && (i < STRLEN - 1)); i++)
#else
	for (i = 0; ((str[i] != '\0') && (i < STRLEN - 1)); i++)
#endif
		star->name[i] = str[i];
	star->name[i] = '\0';
	replace_spaces(star->name);


	if (pf_getval(pfp, "FILTER", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read FILTER from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_valtostr(value, star->filter) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtostr FILTER ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}

	if (pf_getval(pfp, "EXPTIME", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read EXPTIME from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_valtodouble(value, &(star->exptime)) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtodouble EXPTIME ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}

	if (pf_getval(pfp, "AIRMASS", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read AIRMASS from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_valtodouble(value, &(star->airmass)) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtodouble AIRMASS ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}

	/* note the slightly different way of handling errors for weather
	   information - if we can't get values from the environment, just
	   set them to -1 and don't print any error messages.  I do this
	   because sometimes we don't have access to the weather station
	   and it's bad to terminate execution when we can certainly 
	   calculate the "seeing" and "twinkle" parameters without this
	   information.    MWR 5/30/92 */

	if (pf_getval(pfp, "WINDSPEE", value) == NULL) {
		star->windspeed = -1;
	}
	else if (pf_valtodouble(value, &(star->windspeed)) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtodouble WINDSPEE ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_getval(pfp, "HUMIDITY", value) == NULL) {
		star->humidity = -1;
	}
	else if (pf_valtodouble(value, &(star->humidity)) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtodouble HUMIDITY ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_getval(pfp, "TELFOCUS", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read TELFOCUS from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_valtodouble(value, &(star->telfocus)) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtodouble TELFOCUS ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_getval(pfp, "TELTEMP", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read TELTEMP from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_valtodouble(value, &(star->teltemp)) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtodouble TELTEMP ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
#if 0
	if (pf_getval(pfp, "MAGNITUD", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read MAGNITUD from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_valtodouble(value, &(star->real_mag)) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtodouble MAGNITUD ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
#else
	star->real_mag = -99.0;
#endif

	if (pf_getval(pfp, "UT", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read UT from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_valtostr(value, str) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtostr UT ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
	if (read_timestr(str, &(star->ut)) != 0) {
		fprintf(stderr, "%s.slurp_header: can't read_timestr UT ..%s.. from file %s\n",
					progname, str, file);
		pf_close(pfp);
		return(-1);
	}

	pf_close(pfp);
	return(0);
}
	

	/* read in the data from a .coo file, sticking the peak, FWHM, etc.
	   into the given star structure. return 0 if all OK, or -1 if there's
	   a problem. */

static int 
slurp_coo_data(file, star)
char *file;
struct s_star *star;
{
	int i, found, naper;
	char line[LINELEN];
	FILE *fp;
	struct s_coo coo;

	/* first, skip past all lines in the file until we come to the
	   one that consists of 'END' - leave us just past that line */
	if ((fp = fopen(file, "r")) == NULL) {
		fprintf(stderr, "%s.slurp_coo_data: can't open file %s\n",
					progname, file);
		return(-1);
	}
	found = 0;
	while (fgets(line, LINELEN, fp) != NULL) {
		if (strcmp(line, "END\n") == 0) {
			found = 1;
			break;
		}
	}
	if (found == 0) {
		fclose(fp);
		fprintf(stderr, "%s.slurp_coo_data: file %s has no END\n",
					progname, file);
		return(-1);
	}

	/* now that we're here, read each line in, and if it has brighter
	   values than any other (based on the peak value), 
	   replace the brightest values with its. 
	   Thus, we ASSUME that the brightest star must be the desired
	   standard */
	star->dat.c.peak = 0;
	while (fgets(line, LINELEN, fp) != NULL) {
		if (read_coo_line(line, &coo) != 0) {
			fprintf(stderr, "%s.slurp_coo_data: bad data line\n", progname);
		}
		if (coo.peak > star->dat.c.peak) {
			star->dat.c.peak = coo.peak;
			star->dat.c.fwhm = coo.fwhm;
			star->dat.c.round = coo.round;
			star->dat.c.sharp = coo.sharp;
		}
	}
	fclose(fp);

	return(0);
}
		
	/* read in the data from a .pht file, sticking the instrumental
	   magnitudes (together with their estimated uncertainties) into
	   the given star structure. return 0 if all OK, or -1 if there's
	   a problem. */

static int 
slurp_pht_data(file, star)
char *file;
struct s_star *star;
{
	int i, found, naper;
	char line[LINELEN];
	FILE *fp;
	struct s_pht pht;

	/* first, skip past all lines in the file until we come to the
	   one that consists of 'END' - leave us just past that line */
	if ((fp = fopen(file, "r")) == NULL) {
		fprintf(stderr, "%s.slurp_pht_data: can't open file %s\n",
					progname, file);
		return(-1);
	}
	found = 0;
	while (fgets(line, LINELEN, fp) != NULL) {
		if (strcmp(line, "END\n") == 0) {
			found = 1;
			break;
		}
	}
	if (found == 0) {
		fclose(fp);
		fprintf(stderr, "%s.slurp_pht_data: file %s has no END\n",
					progname, file);
		return(-1);
	}

	/* now that we're here, read each line in, and if it has brighter
	   values than any other (based on the largest aperture available), 
	   replace the brightest values with its. 
	   Thus, we ASSUME that the brightest star must be the desired
	   standard */
	for (i = 0; i < MAX_APER; i++) {
		star->dat.p.inst_mag[i] = NODATA;
		star->dat.p.inst_err[i] = NODATA;
	}
	while (fgets(line, LINELEN, fp) != NULL) {
		if (read_pht_line(line, &pht, &naper) != 0) {
			fprintf(stderr, "%s.slurp_pht_data: bad data line\n", progname);
		}
		if (pht.mag[naper - 1] < star->dat.p.inst_mag[naper - 1]) {
			for (i = 0; i < MAX_APER; i++) {
				star->dat.p.inst_mag[i] = pht.mag[i];
				star->dat.p.inst_err[i] = pht.err[i];
			}
			star->dat.p.n_aper = naper;
		}
	}
	fclose(fp);

	return(0);
}
		

	/* read in one line's worth of .coo data. 
	   return 0 if all is well, or -1 if there's a problem.  */
	
static int 
read_coo_line(line, coo)
char *line;
struct s_coo *coo;
{
	if (sscanf(line, "%d %lf %lf %lf %lf %lf %lf ", &(coo->id), &(coo->xc),
			&(coo->yc), &(coo->peak), &(coo->fwhm), &(coo->round), 
			&(coo->sharp)) != 7) {
		fprintf(stderr, "%s.read_coo_line: bad line ..%s..\n",
				progname, line);
		return(-1);
	}
	return(0);
}

	/* read in one line's worth of .pht data. place the	number of 
	   apertures found into the 'naper' parameter, and print a warning
	   message to stdout if there are more than MAX_APER. 

	   return 0 if all is well, or -1 if there's a problem.  */
	
static int 
read_pht_line(line, pht, naper)
char *line;
struct s_pht *pht;
int *naper;
{
	int i, pos;
	double foo, bar;

	if (sscanf(line, "%d %lf %lf %lf %lf ", &(pht->id), &(pht->xc),
			&(pht->yc), &(pht->sky), &(pht->skysig)) != 5) {
		fprintf(stderr, "%s.read_pht_line: bad line ..%s..\n",
				progname, line);
		return(-1);
	}
	
	pos = PHT_POS;
	for (i = 0; i < MAX_APER; i++) {
		if (sscanf(line + pos, "%lf %lf ", (&pht->mag[i]), &(pht->err[i])) != 2) {
			if (i == 0) {
				fprintf(stderr, "%s.read_pht_line: no mag data in line ..%s..\n",
						progname, line);
				return(-1);
			}
			else {
				break;
			}
		}
		pos += PHT_LEN;
	}
	*naper = i;

	/* check to see if there's more data - issue a warning message if there
	   is so that the user knows we're skipping some of it */
	if (sscanf(line + pos, "%lf %lf ", &foo, &bar) == 2) {
		fprintf(stderr, "%s.read_pht_line: mag/err pairs after number %d ignored\n",
				progname, MAX_APER);
	}
	return(0);
}

	/* print out the info in a star structure */

static void
print_star(s)
struct s_star *s;
{
	int i;

	printf("%-15s %-7s %6.2lf %4.2lf %2.0lf %2.0lf %5.1lf %4.1lf %5.2lf %5.2lf",
			s->name, s->filter, s->exptime, s->airmass, s->windspeed, 
			s->humidity, s->telfocus, s->teltemp, s->real_mag, s->ut);
	if (s->type == PHT) {
		printf("%2d ", s->dat.p.n_aper);
		for (i = 0; i < s->dat.p.n_aper; i++) {
			printf("%5.2lf %5.3lf  ", s->dat.p.inst_mag[i], s->dat.p.inst_err[i]);
		}
	}
	else {
		printf("%5d %6.3lf %6.3lf %6.3lf", s->dat.c.peak, s->dat.c.fwhm, 
				s->dat.c.round, s->dat.c.sharp);
	}

	printf("\n");

}

	/* process all the .pht data */

static void 
proc_pht(flags)
int flags;
{
	int i;

	if (flags & PRINT)
		do_print();

	if (flags & SCATTER)
		do_scatter();

	if (flags & EXTINCT)
		do_extinct();

}

	/* print out the data on all the stars in the list, using the 
	   NODATA value if there's no good data */

static void
do_print()
{
	int i;

	sortem(BY_NAME);

	for (i = 0; i < num_stars; i++)
		print_star(&(stars[i]));
}


	/* go through the entire star list: for each set of observations of a 
	   star through the same filter, figure out the extinction by fitting
	   a line of form
	                       mag = a + b(airmass)
	   to the data.

	   print out lines of the form

	           real_mag  filter  num  a[n] b[0] ... a[n] b[n]
	   where
	        real_mag       is the real mag of star through the given filter
	        filter         is the filter
	        num            is the number of measurements compared together
	        a[n]           is the intercept of line for n'th aperture
	        b[n]           is the slope of line for n'th aperture
	 */


static void
do_extinct()
{
	static int first = 0;
	int i, j, k, current, found, num;
	double a, b;
	double e_x[MAX_STAR], *e_y[MAX_APER];
	struct s_star *s, *si;

	if (first == 0) {
		first = 1;
		for (i = 0; (i < num_stars) && (stars[i].type != PHT); i++)
			;
		if (i == num_stars)
			return;
		for (j = 0; j < stars[i].dat.p.n_aper; j++) {
			if ((e_y[j] = (double *) calloc(MAX_STAR, sizeof(double))) == NULL) {
				fprintf(stderr, "%s.do_extinct: can't calloc for %d doubles\n",
							progname, MAX_STAR);
				exit(-1);
			}
		}
	}

	sortem(BY_NAME);
		
	for (i = 0; (i < num_stars) && (stars[i].type != PHT); i++)
		;
	if (i == num_stars)
		return;

	s = &(stars[i]);
	e_x[0] = s->airmass;
	for (j = 0; j < s->dat.p.n_aper; j++)
		e_y[j][0] = s->dat.p.inst_mag[j];	
	num = 0;
	current = 0;
	if (debug_flag)
		print_star(s);
	for (i = 1; i < num_stars; i++) {
		si = &(stars[i]);
		if (si->type != PHT)
			continue;
		if (debug_flag)
			print_star(si);
		found = 0;
		if (strcmp(s->name, si->name) == 0) {
			if (strcmp(s->filter, si->filter) == 0) {
				found = 1;
				current = 1;
				for (j = 0; j < s->dat.p.n_aper; j++) {
					if (si->dat.p.n_aper <= j)
						continue;
					if (si->dat.p.inst_mag[j] == NODATA)
						continue;
					num++;
					e_x[num] = si->airmass;
					e_y[j][num] = si->dat.p.inst_mag[j];
				}
			}
		}
		if ((current == 1) && (found == 0)) {
			printf("%5.2lf  %s  %3d ", s->real_mag, s->filter, num + 1);
			for (j = 0; j < s->dat.p.n_aper; j++) {
				if (fitline(e_x, e_y[j], num + 1, &a, &b) == 0) {
					printf("%6.3lf %6.3lf  ", a, b);
				}
				else {
					printf("%6.3lf %6.3lf  ", NODATA, NODATA);
				}
			}
			printf("\n");
			current = 0;
			num = 0;
			e_x[num] = si->airmass;
			for (j = 0; j < si->dat.p.n_aper; j++)
				e_y[num][j] = si->dat.p.inst_mag[j];
			s = si;
		}
		else if ((current == 1) && (found == 1)) {
			;		/* do nothing */
		}	
		else if ((current == 0) && (found == 0)) {
			num = 0;
			s = si;
			e_x[num] = s->airmass;
			for (j = 0; j < s->dat.p.n_aper; j++) {
				e_y[j][num] = s->dat.p.inst_mag[j];
			}
		}
	}
	if (current == 1) {
		printf("%5.2lf  %s  %3d ", s->real_mag, s->filter, num + 1);
		for (j = 0; j < s->dat.p.n_aper; j++) {
			if (fitline(e_x, e_y[j], num + 1, &a, &b) == 0) {
				printf("%6.3lf %6.3lf  ", a, b);
			}
			else {
				printf("%6.3lf %6.3lf  ", NODATA, NODATA);
			}
		}
		printf("\n");
	}
}

	/* go through the entire star list: for each set of observations of a 
	   star through the same filter which were made within AVG_TIME of
	   each other (i.e. for each set of three short exposures made as a 
	   single group), calculate the scatter among the measurements
	   in each aperture.  print out lines of the form

	           real_mag  time  filter  num  scatter[0] ... scatter[n]
	   where
	        real_mag       is the real mag of star through the given filter
	        time           is the ending UT of the group of observations
	        filter         is the filter
	        num            is the number of measurements compared together
	        scatter[n]     is the scatter of the observations, in mags, through
	                          the n'th aperture
	 */

static void
do_scatter()
{
	int i, j, current, found;
	double x, sum[MAX_APER], sumsq[MAX_APER], num[MAX_APER];
	struct s_star *s, *si;

	sortem(BY_NAME);

	for (i = 0; (i < num_stars) && (stars[i].type != PHT); i++)
		;
	if (i == num_stars)
		return;
	s = &(stars[i]);
	current = 0;
	for (j = 0; j < MAX_APER; j++) {
		sum[j] = s->dat.p.inst_mag[j];
		sumsq[j] = sum[j]*sum[j];
		num[j] = 1.0;
	}
	if (debug_flag)
		print_star(s);
	for (i = 1; i < num_stars; i++) {
		si = &(stars[i]);
		if (si->type != PHT)
			continue;
		if (debug_flag)
			print_star(si);
		found = 0;
		if (strcmp(s->name, si->name) == 0) {
			if (strcmp(s->filter, si->filter) == 0) {
				if (fabs(s->ut - si->ut) < AVG_TIME) {
					found = 1;
					current = 1;
					for (j = 0; j < s->dat.p.n_aper; j++) {
						if (si->dat.p.n_aper <= j)
							continue;
						if (si->dat.p.inst_mag[j] == NODATA)
							continue;
						x = si->dat.p.inst_mag[j];
						sum[j] += x;
						sumsq[j] += x*x;
						num[j] += 1;
					}
				}
			}
		}
		if ((current == 1) && (found == 0)) {
			printf("%5.2lf  %5.2lf  %s  %2d ", s->real_mag, s->ut, s->filter,
						s->dat.p.n_aper);
			for (j = 0; j < s->dat.p.n_aper; j++) {
				x = sum[j]/num[j];
				x = sqrt((sumsq[j] - num[j]*x*x)/(num[j] - 1.0));
				printf("%5.3lf ", x);
			}
			for (j = 0; j < si->dat.p.n_aper; j++) {
				sum[j] = si->dat.p.inst_mag[j];
				sumsq[j] = sum[j]*sum[j];
				num[j] = 1.0;
			}
			printf("\n");
			current = 0;
			s = si;
		}
		else if ((current == 1) && (found == 1)) {
			;		/* do nothing */
		}	
		else if ((current == 0) && (found == 0)) {
			s = si;
			for (j = 0; j < si->dat.p.n_aper; j++) {
				sum[j] = si->dat.p.inst_mag[j];
				sumsq[j] = sum[j]*sum[j];
				num[j] = 1.0;
			}
		}
	}
	if (current == 1) {
		printf("%5.2lf  %5.2lf  %s  %2d ", s->real_mag, s->ut, s->filter,
					s->dat.p.n_aper);
		for (j = 0; j < s->dat.p.n_aper; j++) {
			x = sum[j]/num[j];
			x = sqrt((sumsq[j] - num[j]*x*x)/(num[j] - 1.0));
			printf("%5.3lf ", x);
		}
		printf("\n");
	}
		
}

	
	/* sort stars by some key (or keys) */

static void
sortem(flag)
int flag;
{
	if (flag == BY_TIME) {
		qsort((char *) &(stars[0]), num_stars, sizeof(struct s_star), 
				comp_time);
	}
	else if (flag == BY_NAME) {
		qsort((char *) &(stars[0]), num_stars, sizeof(struct s_star), 
				comp_name);
	}
	else {
		fprintf(stderr, "%s.sortem: unknown sort flag %d\n", progname, flag);
	}
}
		

	/* compare two stars by name and secondarily, by filter, and 
	   tertially, by time. */

static int
comp_name(s1, s2)
struct s_star *s1, *s2;
{
	int ret;

	if ((ret = strcmp(s1->name, s2->name)) != 0)
		return(ret);
	else if ((ret = strcmp(s1->filter, s2->filter)) != 0)
		return(ret);
	else if (s1->ut < s2->ut)
		return(-1);
	else if (s1->ut == s2->ut)
		return(0);
	else
		return(1);
}

	/* compare two stars by UT time */

static int
comp_time(s1, s2)
struct s_star *s1, *s2;
{
	if (s1->ut < s2->ut)
		return(-1);
	else if (s1->ut == s2->ut)
		return(0);
	else
		return(1);
}


	/* replace all occurences of spaces ' ' by underscores '_' in the given
	   string */

static void
replace_spaces(str)
char *str;
{
	int i;

	for (i = 0; (i < strlen(str)) && (str[i] != '\0'); i++)
		if (str[i] == ' ')
			str[i] = '_';
}