/* calculate the airmass for SN 1994I only, given JD on 
	     command line */

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

#define LATITUDE   37.9183  /* latitude of observatory (this for Leuschner) */
#define RASTR     "13:27:47.6"  /* declination of object in question */
#define DECSTR    "47:26:59"    /* declination of object in question */
#ifndef PI
#define PI 3.14159265359
#endif  PI
#define DEGTORAD(x)   ((x)*(PI/180.0))
#define RADTODEG(x)   ((x)*(180.0/PI))
#define MIN_ALT        5.0  /* minimum alt for which airmass calculated */
#define TINY        1.0e-6

#undef MAIN

#ifdef MAIN


main(argc, argv)
int argc;
char *argv[];
{
	double ra, dec, epoch, jra, jdec, airmass, alt, ha;
	double jd, lst;

	if (argc != 2) {
		fprintf(stderr, "usage: air94 JDdate\n");
		exit(1);
	}

	
	get_ra(&ra);
	get_dec(&dec);
	get_epoch(&epoch);
	to_j2000(epoch, ra, dec, &jra, &jdec);

	if (sscanf(argv[1], "%lf", &jd) != 1) {
		fprintf(stderr, "bad arg for JD: ..%s..\n", argv[1]);
		exit(1);
	}
	jd += 2440000.0;

	get_lst(jd, &lst);
	calc_airmass(jra, jdec, lst, &ha, &airmass, &alt);
	printf("%5.3lf\n", airmass);

}
#else
double
air94i(mjd)
double mjd;
{
	double ra, dec, epoch, jra, jdec, airmass, alt, ha;
	double jd, lst;
	
	get_ra(&ra);
	get_dec(&dec);
	get_epoch(&epoch);
	to_j2000(epoch, ra, dec, &jra, &jdec);

	jd = mjd + 2440000.0;

	get_lst(jd, &lst);
	calc_airmass(jra, jdec, lst, &ha, &airmass, &alt);
	return(airmass);

}
#endif









	/* get the RA, using the previous value if the user just
	   types <return>. */

get_ra(ra)
double *ra;
{
	char line[100];
	int nhh, nhm, num;
	double nhs, hhh;
	static int hh, hm;
	static double hs, oldra;

	sprintf(line, RASTR);
	if ((num = sscanf(line, "%d:%d:%lf", &nhh, &nhm, &nhs)) < 1) {
		*ra = oldra;
		return(0);
	}
	else if (num == 1) {
		hh = nhh;
		hm = 0;
		hs = 0.0;
	}
	else if (num == 2) {
		hh = nhh;
		hm = nhm;
		hs = 0.0;
	}
	else {
		hh = nhh;
		hm = nhm;
		hs = nhs;
	}
	hmstohhh(hh, hm, hs, &hhh);
	hhhtodeg(hhh, ra);	
	oldra = *ra;
}
	
	/* get the Dec, using the previous value if the user just
	   types <return>. */

get_dec(dec)
double *dec;
{
	char line[100];
	int nhh, nhm, num;
	double nhs, hhh;
	static int hh, hm;
	static double hs, olddec;

	sprintf(line, DECSTR);
	if ((num = sscanf(line, "%d:%d:%lf", &nhh, &nhm, &nhs)) < 1) {
		*dec = olddec;
		return(0);
	}
	else if (num == 1) {
		hh = nhh;
		hm = 0;
		hs = 0.0;
	}
	else if (num == 2) {
		hh = nhh;
		hm = nhm;
		hs = 0.0;
	}
	else {
		hh = nhh;
		hm = nhm;
		hs = nhs;
	}
	hmstohhh(hh, hm, hs, dec);
	olddec = *dec;
}

	/* get the Epoch, using the previous value if the user just
	   types <return>. */

get_epoch(epoch)
double *epoch;
{
	char line[100];
	int num;
	double newepoch;
	static double oldepoch;

	sprintf(line, "2000.0");
	if ((num = sscanf(line, "%lf", &newepoch)) < 1) {
		*epoch = oldepoch;
	}
	else if (num == 1) {
		*epoch = newepoch;
	}
	oldepoch = *epoch;
}

	/* get the JD - again using the previous value if desired */

this_get_jd(jd)
double *jd;
{
	char line[100];
	int nm, nd, ny, nh, num;
	static int day, month, year, hour, minute;
	double ut;
	static double oldjd;

	while (1) {
		printf("UT Date [MM DD YY] [%02d %02d %02d]: ", month, day, year);
		gets(line);
		if ((num = sscanf(line, "%d %d %d", &nm, &nd, &ny)) < 1) {
			break;
		}
		else if (num != 3) {
			fprintf(stderr, "you must enter all three values \n");
			continue;
		}
		else {
			year = ny + 1900;
			month = nm;
			day = nd;
			break;
		}
	}

	while (1) {
		printf("UT [%02d:%02d]: ", hour, minute);
		gets(line);
		if ((num = sscanf(line, "%d:%d", &nh, &nm)) < 1) {
			break;
		}
		else if (num != 2) {
			fprintf(stderr, "you must enter both values \n");
			continue;
		}
		else {
			hour = nh;
			minute = nm;
			break;
		}
	}

	ut = hour + ((double) minute)/60.0;
	get_jd(day, month, year, ut, jd);
	oldjd = *jd;
}

	/* calculate the airmass, altitude and HA of an object at the
	   given coordinates at the given LST.  The observatory latitude
	   is hardwired in here - see the #define at the top of the
	   program. */

calc_airmass(jra, jdec, lst, ha, airmass, alt)
double jra, jdec, lst, *ha, *airmass, *alt;
{
	double x, az;

	hhhtodeg(lst, &x);
	*ha = x - jra;
	degtohhh(*ha, &x);
	if (x > 12)
		x -= 24;
	else if (x < -12)
		x += 24;
	this_to_altaz(x, jdec, alt, &az, airmass);
}

	/* 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
this_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, z, xx;
	
	lat = LATITUDE;

	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 = 99.0;
	else {
		z = 1.0/a;
		xx = z - 0.0018167*(z - 1.0) - 0.002875*(z - 1.0)*(z - 1.0) 
				- 0.008083*(z - 1.0)*(z - 1.0)*(z - 1.0);
		*airmass = xx;
	}

	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);
}