/*
 *  photom: a package to perform photometric calibration
 *  Copyright (C) 2001  Michael William Richmond
 *
 *  Contact: Michael William Richmond
 *           Physics Department
 *           Rochester Institute of Technology
 *           85 Lomb Memorial Drive
 *           Rochester, NY  14623-5603
 *           E-mail: mwrsps@rit.edu
 *
 *  
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *  
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *  
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */



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

#include <stdio.h>
#include <math.h>
#include "photom_time.h"
#include "airmass.h"

#ifndef M_PI
#define M_PI 3.14159265359
#endif  
#define DEGTORAD(x)   ((x)*(M_PI/180.0))
#define RADTODEG(x)   ((x)*(180.0/M_PI))
#define MIN_ALT        5.0		/* minimum alt for which airmass calculated */
#define TINY        1.0e-6


static void calc_airmass(double jra, double jdec, double lst, double latitude,
                              double *ha, double *airmass, double *alt);
static int this_to_altaz(double ha, double dec, double latitude,
                             double *alt, double *az, double *airmass);



	/***********************************************************************
	 * PROCEDURE: tass_airmass
	 *
	 * DESCRIPTION: given an object's RA and Dec (J2000), 
	 *              and Julian Date of observation, calculate the object's
	 *              airmass.   
	 *                  - calc LST for location
	 *                  - then calc hour angle for object ...
	 *                  -  ... and convert to azimuth and elevation, then 
	 *                  - use airmass = secant(zenith angle) approximation
	 * 
 	 * RETURNS:
	 *     value of airmass
	 *                  
	 */

double
tass_airmass
	(
	double ra,           /* I: J2000 RA of object */
	double dec,          /* I: J2000 Dec of object */
	double jd,           /* I: Julian Date of observation */
	double latitude,     /* I: latitude of observatory (degrees) */
	double longitude     /* I: longitude of observatory  */
	                     /*        (degrees West of Greenwich) */
	)
{
	double lst;
	double airmass, alt, ha;

	/* this calculates LST in decimal hours */
	get_lst(jd, longitude, &lst);

	calc_airmass(ra, dec, lst, latitude, &ha, &airmass, &alt);
#ifdef DEBUG
	printf("HA %5.2lf   altitude %6.2lf   airmass %5.3lf\n", 
				ha, alt, airmass);
#endif

	return(airmass);
}


	/*****************************************************************
	 * PROCEDURE: calc_airmass 
	 *
	 * DESCRIPTION: Calculate the airmass, altitude and HA of an 
	 * object at the given coordinates, observed from the given 
	 * latitude at the given LST.  
	 *
	 * RETURNS:
	 *      nothing
	 */

static void
calc_airmass
	(
	double jra,           /* I: RA of object (decimal degrees) */
	double jdec,          /* I: Dec of object (decimal degrees) */
	double lst,           /* I: LST of observation (decimal hours) */
	double latitude,      /* I: latitude of observatory (decimal degrees) */
	double *ha,           /* O: hour angle of object (decimal degrees) */
	double *airmass,      /* O: airmass of object */
	double *alt           /* O: altitude of object (decimal degrees) */
	)
{
	double x, az;

	/* 
	 * the LST is in decimal hours, so we must convert into
	 * decimal degrees to perform some calculations 
	 */
	x = 15.0*lst;
	if (x == 360.0) {
		x = 0.0;
	}
	*ha = x - jra;

	/* must convert Hour Angle to decimal hours for the next step */
	x = *ha/15.0;
	if (x > 12)
		x -= 24;
	else if (x < -12)
		x += 24;
	this_to_altaz(x, jdec, latitude, alt, &az, airmass);
}

	/*********************************************************************
	 * PROCEDURE: this_to_altaz
	 *
	 * 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).
	 *
	 * Thanks to Dick Treffers for the original version of this code.
         *
	 * RETURNS:
	 *      0               if all OK, 
	 *      -1              if some error occurs
	 */

static int
this_to_altaz
	(
	double ha,           /* I: hour angle (decimal hours)  */
	double dec,          /* I: declination (decimal degrees) */
	double latitude,     /* I: latitude (decimal degrees) */
	double *alt,         /* O: altitude (decimal degrees) */
	double *az,          /* O: azimuth (decimal degrees) */
	double *airmass      /* O: airmass */
	)
{
	double lat, a;
	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);
}