home *** CD-ROM | disk | FTP | other *** search
- /*******************************************************************************
- ** 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);
- }
-