/*******************************************************************************
** SIDEREAL_TIME.C: Compute the mean sidereal time.                           **
*******************************************************************************/

#include <math.h>

/*******************************************************************************
** For easy formulation we invent some abbreviations.                         **
** They include RAD = PI/180, i.e., radians per circle, and some equinoctiae, **
** like JD1900 for the year 1900 and JD2000 for the year 2000, both in        **
** Julian date for the 1st of January, 0h.                                    **
*******************************************************************************/
#define RAD 0.01745329251994329651

#define JD1900 2415020.
#define JD2000 2451545.

#define JJ100 36525.

double sidereal_time(jd)
double jd;
   {
   double st,t,x;
   double range(),poly();

/*******************************************************************************
** Find time at UTC 0 hours.                                                  **
*******************************************************************************/
   x = 0.5 + floor(jd - 0.5);

/*******************************************************************************
** They are a few fractions of a second off, I don't know which one is right. **
** But I do: The NAVY uses the equinoctium to the year 2000, Meeus that to    **
** the year 1900. Hence the difference, due to precision.                     **
*******************************************************************************/
#ifdef YOU_BELIEVE_THE_NAVY

/*******************************************************************************
** Compute mean sidereal time according to NAVY Almanac.                      **
*******************************************************************************/
   t  = (x - JD2000) / JJ100;
   st = poly(0.2790572733,100.002139,0.0000010776,0.,t);
   st = 24. * (st - floor(st));

/*******************************************************************************
** Add in UTC with fudge factor.                                              **
*******************************************************************************/
   st += (jd - x) * 24. * 1.002737909;

/*******************************************************************************
** Compute correction for equation of equinoxes.                              **
*******************************************************************************/
   st -= 0.00029*sin(range(poly(125.04452,-1934.13626,-0.002071,0.,t))*RAD)/RAD;

#else
/*******************************************************************************
** Compute GAST according to Meeus Book.                                      **
*******************************************************************************/
   t  = (x - JD1900) / JJ100;
   st = poly(0.276919398,100.0021359,0.000001075,0.,t);
   st = 24. * (st - floor(st));
   
/*******************************************************************************
** Add in UTC with fudge factor.                                              **
*******************************************************************************/
   st += (jd - x) * 24. * 1.002737908;

#endif

   if (st > 24.)
      st -= 24.;

   return(st);
   }