/* given two filenames  (all pertaining to SN 1994D)

	               1. one with a set of raw diff mags in b
	               2. one with a set of raw diff mags in v

	   with format
	          date      magSN-A   magB-A  magC-A ...

	   calculate the color-extinction correction factor for the 
	   raw diff b mag.

	   spit out the CORRECTED diff b mags to stdout.

	 */

    /* 11/20/93: relax restriction on simultaneity after MJD 9250.
                 allow measurements from up to +/- 3 nights in order
                 to measure color terms */


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


#define LINELEN 100
#define NMAX 1000
#define KBV   0.033       /* coeff of color-extinction term k''_bv */

double date[4][NMAX], magsn[4][NMAX], magb[4][NMAX];

char progname[] = "kbv94d";

main(argc, argv)
int argc;
char *argv[];
{
	char line[LINELEN + 1];
	int i, j, k, nj, nk, num[4], plusflag;
	double avg, nmag, sum, coeff, nx, snsum, bsum;
	double diff, min, color, x, airmass, air93(), zz;
	FILE *fp;

	if (argc != 3) {
		fprintf(stderr, "usage: kbv93 raw.b raw.v\n");
		exit(1);
	}

	/* read in data from all files */
	for (i = 0; i < 2; i++) {
		if ((fp = fopen(argv[i + 1], "r")) == NULL) {
			fprintf(stderr, "can't open file %s\n", argv[i]);
			exit(1);
		}
		num[i] = 0;
		while (fgets(line, LINELEN, fp) != NULL) {
			if (line[0] == '#')
				continue;
			if (sscanf(line, "%lf %lf %lf", &(date[i][num[i]]), 
						&(magsn[i][num[i]]), &(magb[i][num[i]])) != 3) {
				fprintf(stderr, "bad format line ..%s..\n", line);
				exit(1);
			}
			if (num[i]++ >= NMAX - 1) {
				fprintf(stderr, "too many lines!\n");
				exit(1);
			}
		}
		fclose(fp);
	}

	/* use the dates in FIRST file as fiducials for output */
	for (i = 0; i < num[0]; i++) {
		airmass = air93(date[0][i]);
		plusflag = 0;
		min = 1e6;
		for (j = 0; j < num[1]; j++) {
			diff = fabs(date[0][i] - date[1][j]);
			if (diff < min) {
				min = diff;
			}
		}
		if (min > 0.5) {
#if 0
            if (((date[0][i] <= 9250) && (min < 1.5)) ||
                ((date[0][i] > 9250) && (min < 3.5))) {
#else
            if (min < 1.5) {
#endif

				fprintf(stderr, "using diff %5.2lf days for raw mag2 %7.2lf in file %s\n",
						min, date[0][i], argv[1]);
				plusflag = 1;
			}
			else {
				fprintf(stderr, "no date in raw mag2 file %s for %7.2lf .. skipping\n",
						argv[1], date[0][i]);
				continue;
			}
		}

		/* now average all v-band measurements with the matching time */
		snsum = 0.0;
		bsum = 0.0;
		nx = 0.0;
		for (j = 0; j < num[1]; j++) {
			diff = fabs(date[0][i] - date[1][j]);
			if (diff == min) {
				snsum += magsn[1][j];
				bsum += magb[1][j];
				nx++;
			}
		}
		if (nx < 1) {
			fprintf(stderr, "no v-band values to average for date %lf?! skipping\n",
					date[0][i]);
			continue;
		}
		snsum /= nx;
		bsum /= nx;

#ifdef DEBUG
		zz = (magsn[0][i] - snsum) - (magb[0][i] - bsum);
printf("%7.2lf  SN B-V=%6.2lf  B B-V=%6.2lf  SN-B=%6.2lf\n", 
date[0][i], magsn[0][i] - snsum, magb[0][i]-bsum, zz);
#endif

		/* now, calculate the correction for the b-band values in this
		   line, using the AVERAGE of all v-band values at the matching time */

		x = KBV*airmass*(magsn[0][i] - snsum);
		magsn[0][i]  -= x;
		x = KBV*airmass*(magb[0][i] - bsum);
		magb[0][i]  -= x;


		printf("%7.2lf %6.2lf %6.2lf \n", date[0][i], magsn[0][i], magb[0][i]);
	}

}