/* 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 <string.h>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/* Do we need to include /sys/time.h here in place of timeb.h ?? */
#include "../track/wirotypes.h"
#include "../track/track.h"
#include "../track/wtrack.h"
#include "../track/wiro.h"
extern struct wiro_memory *tinfo;
#include "kep.h"
extern double tlong, tlat, glat;
double sidrlt(), refrac();

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



/* Display rate at which ra and dec are changing
 */

	x = RTS/24.0; /*  arcseconds /radian divided by hours */
	N = x*dradt;
	D = x*ddecdt;
	if( N != 0.0 )
		printf( "dRA/dt %.2lf\"/h, dDec/dt %.2lf\"/h\n", N, D );

/* Load rates into  telesope variables */
	timecopy(CURRENT,&TONOBJ);       /* Time we started on object */
/*	V_RA = dradt * RTD / 15.0 ;
	V_DEC = ddecdt * RTD; */

	V_RA = N / (15.0 * 3600.0);
	V_DEC = D / 3600.0;


/* Do rise, set, and transit times
 */
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 );
D = 0;  /* Turn off refraction so we don't count it twice!  */
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 );

/* OK, before we finish our topocentric stuff, ship off these things 
   to the telescope.  The rest is for grins only.
*/

ORIG_RA = ra * RTD / 15.0;
ORIG_DEC = dec * RTD;

tinfo->motion_type = SOLAR;
HP = 0.0 ;   /* Turn off horizontal parralax - this code does it allready */

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 %.3lfs dDec %.2lf\"\n",
		 D, y, z );
	}
printf( "Topocentric:  Altitude %.3lf deg, ", alt );
printf( "Azimuth %.3lf deg\n", az );

printf( "Topocentric: R.A." );
hms( ra ); 

printf( " Dec." );
dms( dec );
printf( "\n" );
printf( "The constant SOLAR is defined as %d \n", SOLAR);
}