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

#define PF  (3.14159/180.0)				/* conversion degrees to radians */
#define SMALL       1.0e-6				/* a small number */

#define TELESCOPE_SLEW_RATE		1.0		/* in degrees per second */
#define TELESCOPE_SETTLE_TIME	1.0		/* in seconds 			*/
#define TELESCOPE_LATITUDE		38.0	/* in degrees 			*/
#define TELESCOPE_DOME_PENALTY	15.0	/* zing em if they move it */
#define TELESCOPE_DOME_RATE		2.0		/* degees/per second 	*/ 
#define TELESCOPE_DOME_SETTLE	4.0

#define EXTINCTION				.1		/* magnitudes per airmass */
#define READOUT_TIME			10.0	/* CCD readout time 	*/
#define BIG_SECZ				50.0
#define REF_AIR                 1.2		/* constant used in 'largest_ha()' */


	/* given an hour angle and declination, calculate and return
	   the altitude, azimuth and (approximate) airmass of an object.

	   the input hour angle and dec should be expressed in decimal degrees
	   (i.e. HA=0 at meridian, HA=-90 when rising, HA=+90 when setting),
	   and the output altitude and azimuth are in decimal degrees (alt=0
	   means on the horizon, alt=90 means at zenith, az=0 means due
	   north, az=90 due East, az=180 due South, az=270 due West). 

	   if pointing directly at zenith, report azimuth=180. */

altaz(ha, dec, alt, az, psecz)
double ha, dec, *alt, *az, *psecz;
{
	static double sl,cl,latitude;
	double sh,sd,ch,cd,y,temp;

	if(latitude==0.){				/* calculate for once and for all */
		latitude=TELESCOPE_LATITUDE;
		sl=sin(latitude*PF);
		cl=cos(latitude*PF);
	}

	sh=sin(ha*PF);	sd=sin(dec*PF);	
	ch=cos(ha*PF);	cd=cos(dec*PF);	

	temp=sl*sd+cl*cd*ch;			/* the cosine of the zenith angle */
	*alt = 90.0 - acos(temp)/PF;
	
	if(temp>(1.0/BIG_SECZ))
		*psecz=1.0/temp;
	else
		*psecz=BIG_SECZ;

	y=ch*sl-sd*cl/cd;				/*calc telescope azimuth */

	if((fabs(sh)+fabs(y))<SMALL)	/* check for zenith condition */
		*az=180.;
	else
		*az=atan2(sh,y)/PF+180.;
}

	/* given a Dec for an object, return the largest (absolute)
	   value of Hour Angle at which the object may be observed. This will
	   depend on the telescope hard limits in HA, and also on the airmass
	   limits on wishes to place on observations. 

	   The current scheme is: only observe an object when its airmass would
	   be less than

	           max airmass = sqrt(airmass*airmass + REF_AIR*REF_AIR) 
	*/

double
largest_ha(dec)
double dec;
{
	double ha, alt, az, airmass, dha, min_airmass;
	
	ha = 0.0;
	dha = 1.0;		/* make larger for faster, smaller for more precise */

	altaz(ha, dec, &alt, &az, &airmass);
	min_airmass = airmass;
	for (ha = dha; ha < MAX_HA_LIMIT; ha += dha) {
		altaz(ha, dec, &alt, &az, &airmass);
		if (airmass > sqrt(min_airmass*min_airmass + REF_AIR*REF_AIR))
			return(ha - dha);
	}
	return(MAX_HA_LIMIT);
}