/* this program figures out which stars in the extinction star file
	   to observe (at a given time) in order to make some measurement
	   of the extinction.  It tries to find stars that range between airmass
	   of 1.0 and 3.0, or 5.0 at most.  Once it has decided on a set of stars,
	   it generates a request file that will take all their pictures in
	   all filters UBVRI, adjusting the exposure time accordingly. */

	/* usage:
	            genextinct [help] [filters=] [ut=] [num=] [outfile=] [print]
	   where
	            filters=BVR        which filters to use (default is UBRVI)
	            ut=4:56:43         UT of observations (default is current time)
	            num=3              number of shots per filter per star (def 1)
	            outfile=ext1.rqs   name of output file (def is 'extinct.rqs')
				print              print out list of stars
	 */
	

	/* fixed bug that improperly printed out DEC lines with values between
	   -1 and 0.   MWR 6/19/92 */

	/* increased exposure times for all bands by factor of about 5 (except
	   for U band) ; also added "exit(0)" statement to fix non-zero 
	   exit codes. MWR 6/20/92 */

	/* added 'sunrise' and 'sunset' flags - they prevent stars from the
	   sunlit side of the sky from being chosen.  If neither is given,
	   stars on either side of the meridian may be chosen.   MWR 7/2/92 */

	/* added "GUIDEMOD= 'tryguide'" to each stanza of the request file  
	   MWR 7/3/92 */

	/* changed format of input file slightly, and decreased exposure times
	   for all filters somewhat.  All this in preparation for using a new
	   set of extinction stars, much brighter than before.  MWR 7/24/92 */

	/* added a check for valid Hour Angle to stars.  Note that if a single
	   observation is scheduled at an invalid HA, the entire request 
	   (from that time on) is cancelled!  MWR 7/25/92 */

	/* modified exposure times slightly, and added MIN_EXPTIME.  MWR 7/25/92 */

	/* modified exposure times again, increasing them by factors 5-10. 
	      --  MWR 8/1/92 */

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "bait.h"
#include "pf.h"
	      
	/* define whichever of the following machines the code is running on */
#undef BAIT	
#define  PCYGNI

#undef DEBUG

#define STRLEN    256
#define MAXSTAR   500
#ifdef PCYGNI
#define DATAFILE  "/home/bait/genextinct.dat"
#else
#define DATAFILE  "/u/richmond/scope/genextinct.dat"
#endif
#define EPOCH     1989.5			/* epoch of data coords */
#define FILTERS   "UBVRI"			/* default filters */
#define OUTFILE   "extinct.rqs"		/* default output file name */
#define NUM_EXP   1					/* def. number of shots per filter */
#define NOSUN     0					/* don't know if sunrise/set */
#define SUNRISE   1					/* indicates this is near sunrise */
#define SUNSET    2					/* indicates this is near sunset */
#define DEF_LIMEAST   -6			/* default East Hour Angle limit (hours) */
#define DEF_LIMWEST    6			/* default West Hour Angle limit (hours) */
#define MIN_EXPTIME  0.1		/* minimum exposure time, seconds */


static void init();
static double find_lst( /* char *utstr */ );
static void mark_stars( /* double lst */ );
static void pick_stars();
static void print_stars();
static void create_file( /* char *filename, *utstr, *filters, int num_exp */ );
static void do_stanza( /* struct s_star *star, char filter, FILE *fp,
                          int num_exp */ );
static double calc_exp( /* struct s_star *star, char filter */ );
static int filter_index( /* char filter */ );
static void find_utday( /* char *str */ );

	/* this is the appropriate exposure time for sixth-mag star in 
	   each band, in seconds */

#define NUM_BANDS   5
static char bands[NUM_BANDS]        = { 'U',   'B',  'V', 'R',  'I' };
static double sixth_time[NUM_BANDS] = { 25.0, 50.0, 15.0, 7.0, 10.0 };
 
struct s_star {
	char name[STRLEN];
	double ra, dec, epoch;
	double mag[NUM_BANDS];
	double airmass;
	double ha;
	int use_flag;
} stars[MAXSTAR];
	
int nstar;
char progname[] = "genextinct";
char usage[] = "genextinct [filters=] [ut=] [num=] [outfile=] [sunset|sunrise] [print] ";


main(argc, argv)
int argc;
char *argv[];
{
	int i, num_exp, sun_flag, print_flag;
	char filters[STRLEN];
	char utstr[STRLEN];
	char outfile[STRLEN];
	double lst;

	strcpy(filters, FILTERS);
	strcpy(utstr, "");
	strcpy(outfile, OUTFILE);
	num_exp = NUM_EXP;
	sun_flag = NOSUN;
	print_flag = 0;

	for (i = 1; i < argc; i++) {
		if (strncmp(argv[i], "help", 4) == 0) {
			printf("%s \n   create requests for extinction measurements\n", usage);
			exit(0);
		}
		if (strncmp(argv[i], "filters=", 8) == 0) {
			if (strlen(argv[i]) > 8) {
				strcpy(filters, argv[i] + 8);
			}
			else {
				fprintf("filters string too short: ..%s..\n", argv[i]);
				exit(-1);
			}
			continue;
		}
		if (strncmp(argv[i], "ut=", 3) == 0) {
			if (strlen(argv[i]) > 3) {
				strcpy(utstr, argv[i] + 3);
			}
			else {	
				fprintf("UT string too short: ..%s..\n", argv[i]);
				exit(-1);
			}
			continue;
		}
		if (strncmp(argv[i], "num=", 4) == 0) {
			if (strlen(argv[i]) > 4) {
				if (sscanf(argv[i] + 4, "%d", &num_exp) != 1) {
					fprintf(stderr, "bad value for num_exp ..%s..\n", argv[i]);
					exit(-1);
				}
			}
			else {
				fprintf("num_exp string too short: ..%s..\n", argv[i]);
				exit(-1);
			}
		}
		if (strncmp(argv[i], "outfile=", 8) == 0) {
			if (strlen(argv[i]) > 8) {
				strcpy(outfile, argv[i] + 8);
			}
			else {	
				fprintf("outfile string too short: ..%s..\n", argv[i]);
				exit(-1);
			}
			continue;
		}
		if (strcmp(argv[i], "sunrise") == 0) {
			sun_flag = SUNRISE;
			continue;
		}
		if (strcmp(argv[i], "sunset") == 0) {
			sun_flag = SUNSET;
			continue;
		}
		if (strcmp(argv[i], "print") == 0) {
			print_flag = 1;
			continue;
		}
	}

	init();

	/* determine the LST of the desired moment */
	lst = find_lst(utstr);

	mark_stars(lst);
	pick_stars(sun_flag);

	if (print_flag == 1) {
		print_stars();
	}

	create_file(outfile, utstr, filters, num_exp);

	exit(0);
}

	/* initialize the stars array */

static void
init()
{
	int i, j;

	for (i = 0; i < MAXSTAR; i++) {
		stars[i].ra = -1.0;
		stars[i].dec = -1.0;
		stars[i].epoch = -1.0;
		for (j = 0; j < NUM_BANDS; j++) 
			stars[i].mag[j] = 99.0;
		stars[i].use_flag = 0;
	}
	nstar = 0;
}

	/* determine the LST of the desired time, and return it as fractional
	   hours.  If the passed UT string is empty, it is filled with the
	   current UT.  */

static double 
find_lst(utstr)
char *utstr;
{
	int day, month, year, h, m;
	double ut, jd, lst, s;

	get_today(&day, &month, &year);

	if (strcmp(utstr, "") == 0) {
		get_ut(&ut);
		hhhtohms(ut, &h, &m, &s);
		sprintf(utstr, "%02d:%02d:%02.0lf", h, m, s);	
	}
	else {
		if (sscanf(utstr, "%d:%d:%lf", &h, &m, &s) != 3) {
			fprintf(stderr, "error in UT string ..%s..\n", utstr);
			exit(-1);
		}
		hmstohhh(h, m, s, &ut);
	}

	get_jd(day, month, year, ut, &jd);
	get_lst(jd, &lst);

#ifdef DEBUG
	printf("LST is %6.3lf\n", lst);
#endif

	return(lst);
}


	/* read in all stars from the star database file.  
	   Calculate the airmass for all available stars for the given LST and mark
	   them for use.   */

static void
mark_stars(lst)
double lst;
{
	FILE *fp;
	int h, m, ngood;
	char line[STRLEN], name[STRLEN], rastr[STRLEN], decstr[STRLEN];
	char class[STRLEN];
	double v, bv, ub, vr, ri, vi, s, x, ha, alt, az, airmass;

	ngood = 0;

	if ((fp = fopen(DATAFILE, "r")) == NULL) {
		fprintf(stderr, "can't open data file %s\n", DATAFILE);
		exit(-1);
	}
	while (fgets(line, STRLEN, fp) != NULL) {
		if (line[0] == '#') 
			continue;
		if (sscanf(line, "%s %s %s %lf %lf %lf %lf %lf %s", 
				stars[nstar].name, rastr, decstr, 
				&v, &ub, &bv, &vr, &vi, class) != 9) {
			fprintf(stderr, "format error in line ..%s..\n", line);
			exit(-1);
		}
		if (sscanf(rastr, "%d:%d:%lf", &h, &m, &s) != 3) {
			fprintf(stderr, "error in RA string ..%s..\n", rastr);
			exit(-1);
		}
		ratodeg(h, m, s, &(stars[nstar].ra));
		if (sscanf(decstr, "%d:%d:%lf", &h, &m, &s) != 3) {
			fprintf(stderr, "error in RA string ..%s..\n", rastr);
			exit(-1);
		}
		hmstohhh(h, m, s, &(stars[nstar].dec));
		/* check to see if there was a negative sign in the Dec string */
		if (strchr(decstr, '-') != NULL) {
			if (stars[nstar].dec >= 0)
				stars[nstar].dec *= -1.0;
		}

		stars[nstar].epoch = EPOCH;
		stars[nstar].mag[2] = v;			/* V mag */
		stars[nstar].mag[1] = v + bv;		/* B mag */
		stars[nstar].mag[0] = stars[nstar].mag[1] + ub;	/* U mag */
		stars[nstar].mag[3] = v - vr;		/* R mag */
		stars[nstar].mag[4] = v - vi;		/* I mag */
		
		/* now figure out what its airmass is at the given LST */
		/* ignore precession for now */
		hhhtodeg(lst, &x);
		ha = x - stars[nstar].ra;
		degtohhh(ha, &x);
		/* make sure that the Hour Angle is in the range -12 < HA < 12 */
		if ((stars[nstar].ha = x) > 12)
			stars[nstar].ha -= 24;
		else if (stars[nstar].ha < -12)
			stars[nstar].ha += 24;
		to_altaz(x, stars[nstar].dec, &alt, &az, &airmass);
		stars[nstar].airmass = airmass;

		nstar++;
	}

}

	/* pick out only the best of the available stars, trying to get 
	        one star in range   1.0 < x < 1.5
	        one                 1.5 < x < 2.0
	        one                 2.0 < x < 2.5
            one                 2.5 < x < 3.0
            one                 3.0 < x < 3.5
	   if it is not possible to fill any one range, don't worry. */

	/* added 7/2/92: if 'sun_flag == SUNSET', then try to pick stars that
	                      are on the eastern side of the sky (HA < 0)
	                 if 'sun_flag == SUNRISE', try to pick stars that are 
	                      on the western side of the sky (HA > 0)
	                 in other words, try to avoid areas in which the
	                 sky is very bright */

	/* added 7/25/92: make sure that the Hour Angle of each star is
	   within the valid limits, as described in the "ait.config" file  MWR */

static void
pick_stars(sun_flag)
int sun_flag;
{
	int i, j, found;
	double upper, lower, east_ha, west_ha;
	char string[PF_LEN + 1], *config_file;
	pf_handle pf_config;

	/* get the limits on Hour Angle */
	if ((config_file = getenv("CONFIG_FILE")) == NULL) {
		fprintf(stderr, "can't get environment var CONFIG_FILE - using defaults\n");
		east_ha = DEF_LIMEAST;		/* East HA limit, in hours */
		west_ha = DEF_LIMWEST;		/* West HA limit, in hours */
	}
	else {
		if ((pf_config = pf_open(config_file, 1, PF_EXIST)) < 0) {
			fprintf(stderr, "can't get environment var CONFIG_FILE - using defaults\n");
			east_ha = DEF_LIMEAST;		/* East HA limit, in hours */
			west_ha = DEF_LIMWEST;		/* West HA limit, in hours */
		}
		else {
			if (pf_getval(pf_config, "TELLIME", string) == NULL) {
				fprintf(stderr, "can't pf_getval for TELLIME - using default\n");
				east_ha = DEF_LIMEAST;
			}
			else if (pf_valtodouble(string, &east_ha) == PF_ERROR) {
				fprintf(stderr, "can't pf_valtodouble for TELLIME - using default\n");
				east_ha = DEF_LIMEAST; 
			}
			else {
				east_ha /= 15.0;		/* convert degrees to hours */
			}

			if (pf_getval(pf_config, "TELLIMW", string) == NULL) {
				fprintf(stderr, "can't pf_getval for TELLIMW - using default\n");
				west_ha = DEF_LIMWEST;
			}
			else if (pf_valtodouble(string, &west_ha) == PF_ERROR) {
				fprintf(stderr, "can't pf_valtodouble for TELLIMW - using default\n");
				west_ha = DEF_LIMWEST; 
			}
			else {
				west_ha /= 15.0;		/* convert degrees to hours */
			}
			pf_close(pf_config);
		}
	}
	
	/* now, give ourselves a little slack on the HA limits, since the
	   entire set of stars is likely to take around 30 minutes to run,
	   and if we do try to point too far, the entire request is 
	   cancelled */
	east_ha += 1.0;			/* give ourselves one whole extra hour */
	west_ha -= 1.0;			/* ditto */
	
	for (i = 0; i < 5; i++) {
		upper = 1.0 + i*0.5;
		lower = upper + 0.5;
		found = 0;
		for (j = 0; j < nstar; j++) {
			if ((stars[j].airmass > upper) && (stars[j].airmass < lower)) {
				if ((sun_flag == SUNRISE) && (stars[j].ha < 0))
					continue;
				if ((sun_flag == SUNSET) && (stars[j].ha > 0))
					continue;
				if ((stars[j].ha < east_ha) || (stars[j].ha > west_ha))
					continue;
				stars[j].use_flag = 1;
				found = 1;
				break;
			}
		}
#ifdef DEBUG
		if (found == 0)
			printf("no star for range %4.2lf - %4.2lf \n", upper, lower);
#endif
	}
}

	/* print out all the stars chosen for observations */

static void
print_stars()
{
	int i;

	for (i = 0; i < nstar; i++) {
		if (stars[i].use_flag > 0) {
			printf("%8.3lf %8.3lf %5.2lf %4.2lf %7.2lf %d\n", stars[i].ra,
					stars[i].dec, stars[i].mag[2], stars[i].airmass, 
					stars[i].ha, stars[i].use_flag);
		}
	}
}

	/* create the output file, a multi-stanza request file with 
	   entries for each star in each filter.  The first stanza contains
	   a little extra information, including the UT-START of the
	   observation. */

static void 
create_file(filename, utstr, filters, num_exp)
char *filename, *utstr, *filters;
int num_exp;
{
	int i, j, max_prior;
	FILE *fp;
	char string[PF_LEN + 1];
	pf_handle pf_config;

	if ((fp = fopen(filename, "w")) == NULL) {
		fprintf(stderr, "can't open output file %s\n", filename);
		exit(-1);
	}

	if ((pf_config = pf_open(CONFIG_FILE, 1, PF_EXIST)) < 0) {
		fprintf(stderr, "process_header: can't open config file ..%s..\n", 
				CONFIG_FILE);
		exit(-1);
	}
	if (pf_getval(pf_config, "MAXPRIOR", string) < 0) {
		printf("can't get value for MAXPRIOR from config file \n");
		exit(-1);
	}
	if (pf_valtoint(string, &max_prior) < 0) {
		printf("invalid value for MAXPRIOR from config file ..%s..\n", 
				string);
		exit(-1);
	}
	pf_close(pf_config);

	/* information for the first stanza only */
	find_utday(string);
	fprintf(fp, "REQID   = 'extinct' \n");
	fprintf(fp, "DAYSTART= '%s' \n", string);
	fprintf(fp, "DAYEND  = '%s' \n", string);
	fprintf(fp, "UT-START= '%s' \n", utstr);
	fprintf(fp, "NUM-OBS = 1 \n");
	fprintf(fp, "PRIORITY= %d \n", max_prior);

	for (i = 0; i < nstar; i++) {
		if (stars[i].use_flag > 0) {
			for (j = 0; j < strlen(filters); j++) {
				do_stanza(&(stars[i]), filters[j], fp, num_exp);
			}
		}
	}

	fclose(fp);
}
	

	/* print out a request file stanza for the given star in the given
	   band.  Take the number of images in each filter given by the
	   'num_exp' parameter. */

static void
do_stanza(star, filter, fp, num_exp)
struct s_star *star;
char filter;
FILE *fp;
int num_exp;
{
	int i, h, m;
	double hh, s, t;
	char buf[STRLEN];
	
	fprintf(fp, "OBSERVAT= 'Leuschner' \n");
	fprintf(fp, "INSTRUME= 'B.A.I.T.'  \n");
	fprintf(fp, "PROCEDUR= 'photo_all' \n");
	fprintf(fp, "GUIDEMOD= 'tryguide' \n");

	fprintf(fp, "OBJECT  = '%s' \n", star->name);
	degtohhh(star->ra, &hh);
	hhhtohms(hh, &h, &m, &s);
	fprintf(fp, "RA      = '%02d:%02d:%04.2lf' \n", h, m, s);
	hhhtohms(star->dec, &h, &m, &s);
	format_hms(h, m, s, buf);
	fprintf(fp, "DEC     = '%s' \n", buf);
	fprintf(fp, "EPOCH   = %6.1lf \n", star->epoch);
	fprintf(fp, "FILTERS = '%c' \n", filter);
	t = calc_exp(star, filter);
	fprintf(fp, "EXPTIME = %6.2lf \n", t);
	i = filter_index(filter);
	fprintf(fp, "MAGNITUD= %5.3lf \n", star->mag[i]);
	fprintf(fp, "END \n");

	for (i = 1; i < num_exp; i++) {
		fprintf(fp, "FILTERS = '%c' \n", filter);
		fprintf(fp, "END \n");
	}
}
	

	
	/* calculate the exposure time for a star, using its magnitude and 
	   the filter as variables.  Base the exposure time on the 
	   following formula:
                                   b  
                      time = k * ( - ) 
                                   B 
	   where
	              b = brightness of this star in the given band
	              B = brightness of sixth-mag star in given band
	              k = appropriate exposure time for sixth-mag star in
	                     the given band

	   Since all of our stars are very bright, the signal rises basically
	   as the exposure time, and we just want to make sure that we don't
	   saturate the chip.

	   If the calculated time is less than MIN_EXPTIME, then return MIN_EXPTIME.
	 */

static double 
calc_exp(star, filter)
struct s_star *star;
char filter;
{
	int i, found;
	double sixth, x, z, time;
		
	found = 0;
	for (i = 0; i < NUM_BANDS; i++) {
		if (filter == bands[i]) {
			found = 1;
			break;
		}
	}
	if (found == 0) {
		fprintf(stderr, "unknown filter band ..%c..\n", filter);
		exit(-1);
	}

	sixth = pow((double) 10.0, (double) (-6.0*0.4));
	x = pow((double) 10.0, (star->mag[i]*(-0.4)));
	z = (sixth/x);
	time = z*sixth_time[i];
	if (time < MIN_EXPTIME)
		return(MIN_EXPTIME);

	return(time);
}	

	/* return the index of the given filter.  Exit with error message
	   if the filter is not found. */

static int
filter_index(filter)
char filter;
{
	int i;

	for (i = 0; i < NUM_BANDS; i++) {
		if (filter == bands[i]) {
			return(i);
		}
	}
	fprintf(stderr, "unknown filter band ..%c..\n", filter);
	exit(-1);
}

static int days[] = { 0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };

	/* print into the given string the UT date of the next night */

static void
find_utday(str)
char *str;
{
	int day, month, year;
	double ut, jd, jds, jde;

	/* first, figure out the current JD */
	get_today(&day, &month, &year);
	get_ut(&ut);
	get_jd(day, month, year, ut, &jd);

	/* now determine the JD limits of the current or next night */
	night_limits(jd, -TWI_DEGREES, &jds, &jde);

	/* use the UT date corresponding to the JD of the START of night 
	   as the returned date */
	if ((floor(jd + 0.5) - 0.5) < (floor(jds + 0.5) - 0.5)) {
		if ((year % 4 == 0) && (year % 400 != 0))
			days[1] = 29;			/* February has 29 days in leap year */
		if (++day > days[month]) {
			day = 1;
			if (++month > 12) {
				month = 1;
				year++;
			}
		}
	}
	sprintf(str, "%02d/%02d/%04d", day, month, year);
}