/* Apply diurnal aberrations and calculate topocentric
 * altitude and azimuth, given the geocentric apparent
 * right ascension and declination.
 *
 * Ephemeris transit times can be obtained by modifying
 * aa.ini to declare TDT = UT.
 * 
 * This program assumes that deltaT, the difference
 * between Ephemeris Time (TDT) and Universal Time (UT),
 * has already been calculated.  The global variables
 * TDT and UT must already have been loaded with the
 * correct values of TDT and UT, respectively.  See deltat.c.
 *
 * Inputs are polar coordinates:
 * dist is the true geocentric distance in au.
 * ra and dec are in radians.
 *
 * J is the Julian date in UT measure.
 *
 * AA page B60 and D3.
 */

#include "kep.h"
extern double tlong, tlat, glat;


int altaz( pol, J )
double pol[3];
double J;
{
double dec, cosdec, sindec, lha, coslha, sinlha;
double ra, dist, last, alt, az, coslat, sinlat;
double N, D, x, y, z, TPI;

ra = pol[0];
dec = pol[1];
dist = pol[2];
TPI = 2.0*PI;

/* local apparent sidereal time, seconds converted to radians
 */
last = sidrlt( J, tlong ) * DTR/240.0;
lha = last - ra; /* local hour angle, radians */
if( prtflg )
	{
	printf( "Local apparent sidereal time " );
	hms( last );
	printf( "\n" );
	}
/* Display rate at which ra and dec are changing
 */
/*
 *if( prtflg )
 *	{
 *	x = RTS/24.0;
 *	N = x*dradt;
 *	D = x*ddecdt;
 *	if( N != 0.0 )
 *		printf( "dRA/dt %.2f\"/h, dDec/dt %.2f\"/h\n", N, D );
 *	}
 */

diurab( last, &ra, &dec );
/* Do rise, set, and transit times
   trnsit.c takes diurnal parallax into account,
   but not diurnal aberration.  */
lha = last - ra;
trnsit( J, lha, dec );

/* Diurnal parallax
 */
diurpx( last, &ra, &dec, dist );

/* Diurnal aberration
 */
/*diurab( last, &ra, &dec );*/

/* Convert ra and dec to altitude and azimuth
 */
cosdec = cos(dec);
sindec = sin(dec);
lha = last - ra;
coslha = cos(lha);
sinlha = sin(lha);

/* Use the geodetic latitude for altitude and azimuth */
x = DTR * glat;
coslat = cos(x);
sinlat = sin(x);

N = -cosdec * sinlha;
D =  sindec * coslat  -  cosdec * coslha * sinlat;
az = RTD * zatan2( D, N );
alt = sindec * sinlat  +  cosdec * coslha * coslat;
alt = RTD * asin(alt);

/* Correction for atmospheric refraction
 * unit = degrees
 */
D = refrac( alt );
alt += D;

/* Convert back to R.A. and Dec.
 */
y = sin(DTR*alt);
x = cos(DTR*alt);
z = cos(DTR*az);
sinlha = -x * sin(DTR*az);
coslha = y*coslat - x*z*sinlat;
sindec = y*sinlat + x*z*coslat;
lha = zatan2( coslha, sinlha );

y = ra; /* save previous values, before refrac() */
z = dec;
dec = asin( sindec );
ra = last - lha;
y = ra - y; /* change in ra */
while( y < -PI )
	y += TPI;
while( y > PI )
	y -= TPI;
y = RTS*y/15.0;
z = RTS*(dec - z);
if( prtflg )
	{
	printf( "atmospheric refraction %.3f deg  dRA %.3fs dDec %.2f\"\n",
		 D, y, z );
	printf( "Topocentric:  Altitude %.3f deg, ", alt );
	printf( "Azimuth %.3f deg\n", az );

	printf( "Topocentric: R.A." );
	hms( ra );
	printf( " Dec." );
	dms( dec );
	printf( "\n" );
	}
return(0);
}