/* this file includes routines to calculate rising, setting times, as 
	   well as similar routines that determine when an object is more than
	   some distance above the horizon. See Duffett-Smith, "Practical 
	   Astronomy with your Pocket Calculator", and Green, "Spherical
	   Astronomy", for more details. */

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


int obj_riseset( /* double jd, ra, dec, *ut_rise, *ut_set */ );
int obj_limits( /* double jd, ra, dec, alt, *ut_rise, *ut_set, 
					int limits_flag */ );
int to_altaz( /* double ha, dec, *alt, *az, *airmass */ );
void night_limits( /* double jd, degrees, *jds, *jde */ );
int obj_jdlimits( /* double jd, ra, dec, alt, int *flags, double *jds1, *jde1, 
                     *jds2, *jde2 */ );
void calc_jd( /* double jdnow, ut, *jdthen */ );

#ifdef VERBOSE
extern char *progname;
#endif

#ifndef PI
#define PI 3.14159265359
#endif  PI
#define DEGTORAD(x)    ((x)*(PI/180.0))
#define RADTODEG(x)    ((x)*(180.0/PI))

	/* these are used by the to_altaz() routine */
#define TINY          1.0e-6		/* check for division by zero */
#define MIN_ALT       5.0			/* minimum altitude in degrees */
                         			/*    for which sec z calculated */

	/* given a JD and object's (RA, Dec), find the time of object's
	   rise and set in UT. Note that circumpolar objects, and those which
	   never come above the horizon, are indicated by positive return
	   codes.
	   Input RA, Dec are in decimal degrees (J2000.0) and output ut_rise
	   and ut_set are in fractional hours from 0h UT.

	   Note that the rise time might come AFTER the set time - in this case,
	   the object was above the horizon at the given JD, set, and then rose
	   again before 24 hours had passed.

	   The observer's latitude is read from the CONFIG_FILE. It is assumed
	   that refraction lifts the object by 34 arcminutes at the horizon.

	   the routine returns 
	         0          if         object rise and set times have been found
	         ALT_ALWAYS if         object is ALWAYS above horizon 
	         ALT_NEVER  if         object is NEVER above horizon
	        -1          if         an error
	 */

int
obj_riseset(jd, ra, dec, ut_rise, ut_set)
double jd, ra, dec;
double *ut_rise, *ut_set;
{
	double latitude, phi, delta, a, ar, as, h, lstr, lsts;
	double utr, uts, psi, x, da, y, dt;

	if (read_latitude(&latitude) < 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.obj_riseset: can't read latitude\n", progname);
#endif
		return(-1);
	}
	phi = DEGTORAD(latitude);
	delta = DEGTORAD(dec);
	a = ra/15.0;					/* convert ra from degrees to hours */

	ar = sin(delta)/cos(phi);
	if (ar > 1.0)
		return(ALT_ALWAYS);
	else if (ar < -1.0)
		return(ALT_NEVER);
	ar = RADTODEG(acos(ar));		/* rise azimuth */
	as = 360.0 - ar;				/* set azimuth */
	
	h = tan(phi)*tan(delta);
	if (h > 1.0)
		return(ALT_ALWAYS);
	else if (h < -1.0)
		return(ALT_NEVER);
	h = -h;
	h = RADTODEG(acos(h))/15.0;
	lstr = 24.0 + a - h;
	if (lstr > 24.0)
		lstr -= 24.0;
	lsts = a + h;
	if (lsts > 24.0)
		lsts -= 24.0;

	/* now calculate corrections due to the atmosphere */
	psi = acos(sin(phi)/cos(delta));
	x = DEGTORAD(0.566667);		/* assumes 34 minutes of arc refraction */
	da = asin(tan(x)/tan(psi));
	da = RADTODEG(da);			/* this is the change in azimuth, in degrees */
	y = asin(sin(x)/sin(psi));
	y = RADTODEG(y);
	dt = (240.0*y)/(cos(delta));
	dt /= 3600.0;
	lstr -= dt;
	lsts += dt;
	
	/* finally, convert the LST values to UT */
	if (lst2ut(jd, lstr, &utr) < 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.obj_riseset: can't get UT from lst2ut\n", progname);
#endif
		return(-1);
	}
	if (lst2ut(jd, lsts, &uts) < 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.obj_riseset: can't get UT from lst2ut\n", progname);
#endif
		return(-1);
	}

	*ut_rise = utr;
	*ut_set = uts;

	return(0);
}

	/* calculate the UT at which an object rises above the given altitude
	   and at which it sets below the given altitude over the horizon, for
	   the given Julian Date.  Note that the object may rise AFTER it has
	   set, simply indicating that it was already in the sky at the beginning
	   of the day. 
	
	   If limits_flag is set to 1, then look at the Hour Angle limits in the 
	   CONFIG_FILE and figure those into the calculations; do not rely upon 
	   altitude alone.  If limits_flag is NOT set, then go ahead and calculate
	   times regardless of the altitude (this mode is used for getting sunrise-
	   and sunset-type times, when the altitude negative).

	   input is Julian Date, object RA and Dec in decimal degrees, altitude
	   above the horizon in decimal degrees. Output is rise and set times
	   in fractional hours UT. the function returns
	         0          if         object rise and set times have been found
	         ALT_ALWAYS if         object is ALWAYS above horizon 
	         ALT_NEVER  if         object is NEVER above horizon
	        -1          if         an error

	   no correction for refraction or other atmospheric effects is made.

	*/

int
obj_limits(jd, ra, dec, alt, ut_rise, ut_set, limits_flag)
double jd, ra, dec, alt;
double *ut_rise, *ut_set;
int limits_flag;
{
	double latitude, phi, delta, lstr, lsts, z, h, east_lim, west_lim;
	double utr, uts;

	if (read_latitude(&latitude) < 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.obj_limits: can't read latitude\n", progname);
#endif
		return(-1);
	}
	if (read_halimits(&east_lim, &west_lim) < 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.obj_limits: can't read HA limits\n", progname);
#endif
		return(-1);
	}
	phi = DEGTORAD(latitude);
	delta = DEGTORAD(dec);
	z = DEGTORAD(90.0 - alt);

	h = (cos(z) - sin(delta)*sin(phi))/(cos(delta)*cos(phi));
	if (h < -1.0)
		return(ALT_ALWAYS);
	if (h > 1.0)
		return(ALT_NEVER);
	h = RADTODEG(acos(h));
	if ((limits_flag == 1) && (h > fabs(east_lim)))
		lstr = (ra + east_lim)/15.0;		/* remember that east_lim < 0 */
	else
		lstr = (ra - h)/15.0;
	if (lstr < 0.0)
		lstr += 24.0;
	if ((limits_flag == 1) && (h > fabs(west_lim)))
		lsts = (ra + west_lim)/15.0;
	else
		lsts = (ra + h)/15.0;
	if (lsts > 24.0)
		lsts -= 24.0;

	/* finally, convert the LST values to UT */
	if (lst2ut(jd, lstr, &utr) < 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.obj_limits: can't get UT from lst2ut\n", progname);
#endif
		return(-1);
	}
	if (lst2ut(jd, lsts, &uts) < 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.obj_limits: can't get UT from lst2ut\n", progname);
#endif
		return(-1);
	}

	*ut_rise = utr;
	*ut_set = uts;

	return(0);
}

	/* convert the given Hour Angle (in hours) and Dec (decimal degrees)
	   to altitude, azimuth (in decimal degrees), and put an approximate 
	   value for the airmass of the object into the
	   'airmass' parameter (if the object is visible - otherwise just
	   use some big number).  The secant z rule is used to derive the
	   airmass. if the position is at the zenith, return an azimuth of
	   180 degrees (due south).

	   this code stolen from one of Dick Treffers' programs.

	   return 0 if all OK, or -1 if there is an error determining the latitude
	   via read_latitude(). */

int
to_altaz(ha, dec, alt, az, airmass)
double ha, dec, *alt, *az, *airmass;
{
	double lat, a, a2;
	double sl, cl, sh, ch, sd, cd, y;
	
	if (read_latitude(&lat) != 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.to_altaz: can't read latitude \n", progname);
#endif
		return(-1);
	}

	sl = sin(DEGTORAD(lat));
	cl = cos(DEGTORAD(lat));
	sh = sin(DEGTORAD(ha*15.0));
	ch = cos(DEGTORAD(ha*15.0));
	sd = sin(DEGTORAD(dec));
	cd = cos(DEGTORAD(dec));

	a = sl*sd + cl*cd*ch;
	*alt = 90.0 - RADTODEG(acos(a));
	if (*alt < MIN_ALT)
		*airmass = 1.0/cos(DEGTORAD(90.0 - MIN_ALT));
	else
		*airmass = 1.0/a;

	y = ch*sl - (sd*cl)/cd;
	if ((fabs(sh) + fabs(y)) < TINY)
		*az = 180.0;
	else
		*az = 180.0 + RADTODEG(atan2(sh, y));

	return(0);
}

	/* given a fractional Julian Date, as well as the number of degrees
	   the sun must be below the horizon for normal observing to begin
	   (as opposed to twi flats), return in the passed parameters the
	   fractional Julian Dates when the night begins and ends.  The rule is

	       1. if it's currently night, report JD when THIS night began and
	          will end.

	       2. if it's not currently night, report JD when the NEXT night
	          will begin and end. 

	   start of night is always before the end of night, of course. */

void 
night_limits(jd, degrees, jds, jde)
double jd, degrees, *jds, *jde;
{
	int h, m, up;
	double date_jd;
	double ra, dec, ha, alt, az, lst, ra2, dec2, dummy;
	double ut, utr, uts, s;

	date_jd = floor(jd + 0.5) - 0.5;
	ut = (jd - date_jd)*24.0;
	get_lst(jd, &lst);
#ifdef DEBUG
	hhhtohms(lst, &h, &m, &s);
	printf("  for JD = %.2lf, the current LST = %02d:%02d:%02lf\n", jd,
			h, m, s);
#endif

	sun_radec(date_jd, &ra, &dec);
#ifdef DEBUG
	hhhtohms(ra/15.0, &h, &m, &s);
	printf("  JD = %.4lf, the sun is at RA = %02d:%02d:%02lf ", 
				date_jd, h, m, s);
	hhhtohms(dec, &h, &m, &s);	
	printf(" Dec = %+02d:%02d:%02lf\n", h, m, s);
#endif

	ha = lst - ra/15.0;
	if (ha < -12.0)
		ha += 24.0;
	if (ha > 12.0)
		ha -= 24.0;
#ifdef DEBUG
	hhhtohms(ha, &h, &m, &s);
	printf("  the sun's Hour Angle is %02d:%02d:%02lf; it is ", h, m, s);
#endif
	to_altaz(ha, dec, &alt, &az, &dummy);	
	if (alt < -TWI_DEGREES) {
		up = 0;
#ifdef DEBUG
		printf("astronomical night\n");
#endif
	}
	else {
#ifdef DEBUG
		printf("astronomical day\n");
#endif
		up = 1;
	}

	/* now figure out the time at which the current or next "night" begins 
	   (astronomical twilight) and the time at which "night" ends 
	   (astronomical dawn), expressed in JD. */

	obj_limits(date_jd, ra, dec, (double) degrees, &utr, &uts, 0);

	/* first, figure out the JD of start of "night" */
	if (uts < ut) {	
		if (up == 0) {
			*jds = date_jd + uts/24.0;
		}
		else {
			sun_radec(date_jd + 1.0, &ra2, &dec2);
			obj_limits(date_jd + 1.0, ra2, dec2, (double) degrees,
					&dummy, &uts, 0);
#ifdef DEBUG
			hhhtohms(ra2/15.0, &h, &m, &s);
			printf("    using JD = %.2lf, solar RA = %02d:%02d:%02lf ",
					date_jd + 1.0, h, m, s);
			hhhtohms(dec2, &h, &m, &s);
			printf(" Dec = %02d:%02d:%02lf\n", h, m, s);
#endif
			*jds = date_jd + 1.0 + uts/24.0;
		}
	}
	else {
		if (up == 0) {
			sun_radec(date_jd - 1.0, &ra2, &dec2);
			obj_limits(date_jd - 1.0, ra2, dec2, (double) degrees,
					&dummy, &uts, 0);
#ifdef DEBUG
			hhhtohms(ra2/15.0, &h, &m, &s);
			printf("    using JD = %.2lf, solar RA = %02d:%02d:%02lf ",
					date_jd - 1.0, h, m, s);
			hhhtohms(dec2, &h, &m, &s);
			printf(" Dec = %02d:%02d:%02lf\n", h, m, s);
#endif
			*jds = date_jd - 1.0 + uts/24.0;
		}
		else {
			*jds = date_jd + uts/24.0;
		}
	}
#ifdef DEBUG
	printf("  start of astronomical twilight is at JD = %.2lf \n", *jds);
#endif

	/* now figure out end of night, by calculating time of astronomical dawn */
	if (utr < ut) {	
		sun_radec(date_jd + 1.0, &ra2, &dec2);
		obj_limits(date_jd + 1.0, ra2, dec2, (double) degrees, 
				&utr, &dummy, 0);
#ifdef DEBUG
		hhhtohms(ra2/15.0, &h, &m, &s);
		printf("    using JD = %.2lf, solar RA = %02d:%02d:%02lf ",
				date_jd + 1.0, h, m, s);
		hhhtohms(dec2, &h, &m, &s);
		printf(" Dec = %02d:%02d:%02lf\n", h, m, s);
#endif
		*jde = date_jd + 1.0 + utr/24.0;
	}
	else {
		*jde = date_jd + utr/24.0;
	}
#ifdef DEBUG
	printf("   end of astronomical twilight is at JD = %.2lf \n", *jde);
#endif

}
	

	/* for an object with the given RA and Dec, and for a given JD date, 
	   figure out the start and end times of the observing window(s) for
	   the object, in which it is always 'alt' degrees above horizon. 

	   'flags' is an integer variable which is set to:
	         0          object not visible tonight
	         1          object visible in one window tonight
	         2          object visible in two windows tonight

	   'jds1' and 'jde1' are set the JD start and end of the first observing
	   window, if there is one; if not, they are undefined. 'jds2' and 'jde2'
	   are likewise set to the termini of the second observing window.

	   the function returns 0 if all is OK, or -1 if there is a problem */

int 
obj_jdlimits(jd, ra, dec, alt, flags, jds1, jde1, jds2, jde2)
double jd, ra, dec, alt, *jds1, *jde1, *jds2, *jde2;
int *flags;
{
	int err_flag, alt_code;
	double night_start, night_end, night_mid, jd_base, utr, uts, jdr, jds;
	static double alt_limit = 0.0;
	pf_handle pf_config;
	char value[PF_LEN + 1];

	/* find the altitude limit of the telescope (only once) */
	if (alt_limit == 0.0) {
		err_flag = 0;
		if ((pf_config = pf_open(CONFIG_FILE, 1, PF_EXIST)) < 0) {
#ifdef VERBOSE
			fprintf(stderr, "%s.obj_jdlimits: can't open config file %s\n",
					progname, CONFIG_FILE);
#endif
			return(-1);
		}
		if (pf_getval(pf_config, "TELLIMA", value) == NULL) {
#ifdef VERBOSE
			fprintf(stderr, "%s.obj_jdlimits: can't read TELLIMA \n",
					progname);
#endif
			err_flag = 1;
		} 
		else if (pf_valtodouble(value, &alt_limit) < 0) {
#ifdef VERBOSE
			fprintf(stderr, "%s.obj_jdlimits: bad value for TELLIMA ..%s..\n",
						progname, value);
#endif
			err_flag = 1;
		}
		if (pf_close(pf_config) < 0) {
#ifdef VERBOSE
			fprintf(stderr, "%s.obj_jdlimits: can't close config file %s\n",
						progname, CONFIG_FILE);
#endif
			return(-1);
		}
		if (err_flag != 0)
			return(-1);
	}
	if (alt < alt_limit)
		alt = alt_limit;

	/* first, calculate the start and end of the current, or next, 
	   night-time observing period */
	night_limits(jd, -TWI_DEGREES, &night_start, &night_end);
	night_mid = (night_start + night_end)/2.0;
	jd_base = (floor(night_mid + 0.5) - 0.5);

	/* now calculate the rise and set times of the object */
	if ((alt_code = obj_limits(jd, ra, dec, alt, &utr, &uts, 1)) < 0) {
#ifdef VERBOSE
		fprintf(stderr, "%s.obj_jdlimits: can't calculate alt limits \n",
				progname);
#endif
		return(-1);
	}
	else if (alt_code == ALT_ALWAYS) {
		*flags = 1;
		*jds1 = night_start;
		*jde1 = night_end;
		return(0);
	}
	else if (alt_code == ALT_NEVER) {
		*flags = 0;
		return(0);
	}

	/* this object does rise and set. to compare its rise/set times to those
	   of night start/end, convert them to JD form, and put them into 
	   a 12-hour range around midnight by adding/subtracting 
	   23h 56m = 0.9972 days. */
	jdr = jd_base + utr/24.0;
	jds = jd_base + uts/24.0;
	if (jdr < night_mid - 0.5)
		jdr += 0.9972;
	else if (jdr > night_mid + 0.5)
		jdr -= 0.9972;
	if (jds < night_mid - 0.5)
		jds += 0.9972;
	else if (jds > night_mid + 0.5)
		jds -= 0.9972;
#ifdef DEBUG
	printf("jdr %12.2lf  jds %12.2lf  nst %12.2lf  nend %12.2lf \n",
			jdr, jds, night_start, night_end);
#endif

	/* now there are three cases to consider. In the first, the object rises
	   before the start of the night */
	if (jdr < night_start) {
		if (jds < night_start) {		/* sets before night falls */
			*flags = 0;
		}
		else if (jds < night_end) {		/* sets during the night */
			*flags = 1;
			*jds1 = night_start;
			*jde1 = jds;
		}
		else {							/* sets after end of night */
			*flags = 1;
			*jds1 = night_start;
			*jde1 = night_end;
		}
#ifdef DEBUG
		printf("  case 1: jdr < night_start \n");
#endif
		return(0);
	}	

	/* in the second case, the object rises after the night ends */
	if (jdr > night_end) {
		if (jds < night_start) {		/* sets before night falls */
			*flags = 0;
		}
		else if ((jds > night_start) && (jds < night_end)) {	/* dur night */
			*flags = 1;
			*jds1 = night_start;
			*jde1 = jds;
		}
		else if ((jds > night_end) && (jds < jdr)) {	/* up all night */
			*flags = 1;
			*jds1 = night_start;
			*jde1 = night_end;
		}
		else {							/* sets before night falls */
			*flags = 0;
		}
#ifdef DEBUG
		printf("  case 2: jdr > night_end \n");
#endif
		return(0);
	}

	/* in the last case, the object both rises and sets during the night.
	   it is possible for polar objects to visible TWICE during the night,
	   so the checks are more elaborate here */

#ifdef DEBUG
		printf("  case 3: ");
#endif
	*flags = 1;
	*jds1 = jdr;
	if ((jds > jdr) && (jds < night_end)) {		/* rises, sets during night */
		*jde1 = jds;
#ifdef DEBUG
		printf(" A\n");
#endif
	}
	else if ((jds > jdr) && (jds > night_end)) {	/* sets after dawn */
		*jde1 = night_end;
#ifdef DEBUG
		printf(" B\n");
#endif
	}
	else if ((jds < jdr) && (jds < night_start)) {	/* sets after dawn */
		*jde1 = night_end;
#ifdef DEBUG
		printf(" C\n");
#endif
	}
	else if ((jds < jdr) && (jds > night_start)) {	/* sets dur NEXT night */	
#ifdef DEBUG
		printf(" D\n");
#endif
		*flags = 2;
		*jds1 = night_start;
		*jde1 = jds;
		*jds2 = jdr;
		*jde2 = night_end;
	}
	else {			/* shouldn't get here */
#ifdef VERBOSE
		fprintf(stderr, "%s.obj_jdlimits: invalid (?) jdr %.2lf jds %.2lf\n",
				progname, jdr, jds);
#endif
		return(-1);
	}
	
	return(0);
}


	/* calculate the fractional Julian Date which corresponds to this (or the
	   coming) night at the given UT */

void
calc_jd(jdnow, ut, jdthen)
double jdnow, ut, *jdthen;
{
	double jds, jde, date_jd;

	night_limits(jdnow, -TWI_DEGREES, &jds, &jde);
	date_jd = floor(jdnow + 0.5) - 0.5;
	if ((*jdthen = date_jd + ut/24.0) < jds)
		*jdthen = date_jd + 1.0 + ut/24.0;
}