/* 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; double g_az, g_el; 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 ); g_az = az; g_el = alt; // printf( "Topocentric: R.A." ); // hms( ra ); // printf( " Dec." ); // dms( dec ); // printf( "\n" ); } return(0); }