/* calculate the airmass for some object(s) */

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

#define LATITUDE   37.3433		/* latitude of observatory (this for Lick) */
#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


char progname[] = "airmass";

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

	while (1) {
		get_ra(&ra);
		get_dec(&dec);
		get_epoch(&epoch);
		to_j2000(epoch, ra, dec, &jra, &jdec);
		this_get_jd(&jd);
		get_lst(jd, &lst);
		calc_airmass(jra, jdec, lst, &ha, &airmass, &alt);
		printf("HA %5.2lf   altitude %6.2lf   airmass %5.3lf\n", 
				ha, alt, airmass);
	}

}


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

	printf("object  RA [%02d:%02d:%04.1lf]: ", hh, hm, hs);
	gets(line);
	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;

	printf("object Dec [%02d:%02d:%04.1lf]: ", hh, hm, hs);
	gets(line);
	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;

	printf("object Epoch [%6.1lf]: ", oldepoch);
	gets(line);
	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);
}