/* given an .rqs file, calculate the probability that should
	   be assigned to it TONIGHT; the scheme used is

	      first, let
	                           # of obs wanted, total
	             interval = --------------------------------
	                           # of days, total


	      if (days-since-last-observed < interval)

	                                 interval - days since observed 
	             probabil = exp(-6.0*  -----------------------------)
	                                         interval

	      else
	             probabil = 1.0
	

	   if there is no information on when the object was last 
	   observed, then use 1/interval as the probability.

	   fixed 2 bugs 5/21/92 MWR:
	   1. now write INTERVAL into the request file, so that we can read it
	      in and use it as the interval from now on.  So we only calculate
	      this once, instead of doing it every day (which gave 
	      ever-shrinking probabilities, due to ever-shrinking NUM-OBS). 
	   2. fixed exponental calculation - I had been doing it wrong! */

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

char usage[] = "probrqs file.rqs ...";
char progname[] = "probrqs";

static void find_today( /* int *today */ );

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

main(argc, argv)
int argc;
char *argv[];
{
	int i, today;
	char file[100];

	if ((argc == 1) || (strcmp(argv[1], "help")) == 0) {
		fprintf(stderr, "usage: %s - assign probability for request file(s)\n",
				 usage);
		exit(0);
	}

	find_today(&today);

	for (i = 1; i < argc; i++) {
		strcpy(file, argv[i]);
		calc_prob(file, today);
	}
	
	exit(0);
}


calc_prob(file, today)
char *file;
int today;
{
	char value[PF_LEN + 1];
	int day, month, year, dstart, dend, dlast, numobs, interval_flag;
	pf_handle pfp;
	double interval, prob, x;

	interval_flag = 0;

	if ((pfp = pf_open(file, 1, PF_EXIST)) < 0) {
		fprintf(stderr, "%s.calc_prob: can't pf_open file %s\n", progname, file);
		exit(-1);
	}
	if (pf_getval(pfp, "DAYSTART", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read DAYSTART from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (read_date(value, &day, &month, &year) < 0) {
		fprintf(stderr, "%s.slurp_header: bad format for DAYSTART ..%s.. in file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
	dstart = since_1950(day, month, year);

	if (pf_getval(pfp, "DAYEND", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read DAYEND from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (read_date(value, &day, &month, &year) < 0) {
		fprintf(stderr, "%s.slurp_header: bad format for DAYEND ..%s.. in file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
	dend = since_1950(day, month, year);

	if (pf_getval(pfp, "DAYLAST", value) == NULL) {
		dlast = -1;
	}
	else if (read_date(value, &day, &month, &year) < 0) {
		fprintf(stderr, "%s.slurp_header: bad format for DAYLAST ..%s.. in file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}
	else
		dlast = since_1950(day, month, year);

	if (pf_getval(pfp, "NUM-OBS", value) == NULL) {
		fprintf(stderr, "%s.slurp_header: can't read NUM-OBS from file %s\n",
					progname, file);
		pf_close(pfp);
		return(-1);
	}
	if (pf_valtoint(value, &numobs) != 0) {
		fprintf(stderr, "%s.slurp_header: can't valtoint NUM-OBS ..%s.. from file %s\n",
					progname, value, file);
		pf_close(pfp);
		return(-1);
	}

	if (pf_getval(pfp, "INTERVAL", value) != NULL) {
		if (pf_valtodouble(value, &interval) != 0) {
			fprintf(stderr, "%s.slurp_header: can't valtoint NUM-OBS ..%s.. from file %s\n",
						progname, value, file);
			pf_close(pfp);
			return(-1);
		}
		interval_flag = 1;
	}
	else
		interval = -1.0;

	if (numobs == 0) {
		prob = 0.0;
	}
	else {
		if (interval == -1.0)
			interval = (double)(dend - dstart)/(double)numobs;
		if ((dstart == dend) || (dlast == -1)) {
			prob = 1.0;
		}
		else {
			if ((today - dlast) >= interval)
				prob = 1.0;
			else {
				x = ((interval - ((double)today - (double)dlast))/interval);
				prob = exp(-6.0*x);
			}
		}
	}

	if (prob > 1.0)
		prob = 1.0;

	/* now put the value of PROBABIL into the header of the file */
	pf_doubletoval(prob, value);
	if (pf_putval(pfp, "PROBABIL", value) == NULL) {
		fprintf(stderr, "%s.calc_prob: can't pf_putval PROBABIL = ..%s.. in file %s\n",
				progname, value, file);
		pf_close(pfp);
		return(-1);
	}
	/* and, if it wasn't there already, put the INTERVAL value into the
	   header of the file */
	pf_doubletoval(interval, value);
	if (pf_putval(pfp, "INTERVAL", value) == NULL) {
		fprintf(stderr, "%s.calc_prob: can't pf_putval INTERVAL = ..%s.. in file %s\n",
				progname, value, file);
		pf_close(pfp);
		return(-1);
	}
	
	return(0);
}


	/* figure out the UT date for the current, or NEXT, night, and put the
	   number of days from 1950 to that date in the passed parameter. */

static void
find_today(today)
int *today;
{
	int day, month, year;
	long t, time();
	double ut, jd, jds, jde;
	struct tm *tm;

	/* 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++;
			}
		}
	}

	*today = since_1950(day, month, year);
}