Function sidtime_long_term

swephlib.cpp:3290–3329

* The time range of DE431 requires a new calculation of sidereal time that * gives sensible results for the remote past and future. * The algorithm is based on the formula of the mean earth by Simon & alii, * "Precession formulae and mean elements for the Moon and the Planets", * A&A 282 (1994), p. 675/678. * The longitude of the mean earth relative to the mean equinox J2000 * is calc

/*#define SIDT_LTERM   TRUE
#if SIDT_LTERM*/
static double sidtime_long_term(double tjd_ut, double eps, double nut)
{
  double tsid = 0, tjd_et;
  double dlon, xs[6], xobl[6], dhour, nutlo[2];
  double dlt = AUNIT / CLIGHT / 86400.0;
  double t, t2, t3;
  tjd_et = tjd_ut + swe_deltat_ex(tjd_ut, -1, NULL);
  t = (tjd_et - J2000) / 365250.0;
  t2 = t * t; t3 = t * t2;
  /* mean longitude of earth J2000 */
  dlon = 100.46645683 + (1295977422.83429 * t - 2.04411 * t2 - 0.00523 * t3) / 3600.0;
  /* light time sun-earth */
  dlon = swe_degnorm(dlon - dlt * 360.0 / 365.2425);
  xs[0] = dlon * DEGTORAD; xs[1] = 0; xs[2] = 1;
  /* to mean equator J2000, cartesian */
  xobl[0] = 23.45; xobl[1] = 23.45;
  xobl[1] = swi_epsiln(J2000 + swe_deltat_ex(J2000, -1, NULL), 0) * RADTODEG;
  swi_polcart(xs, xs);
  swi_coortrf(xs, xs, -xobl[1] * DEGTORAD);
  /* precess to mean equinox of date */
  swi_precess(xs, tjd_et, 0, -1);
  /* to mean equinox of date */
  xobl[1] = swi_epsiln(tjd_et, 0) * RADTODEG;
  swi_nutation(tjd_et, 0, nutlo);
  xobl[0] = xobl[1] + nutlo[1] * RADTODEG;
  xobl[2] = nutlo[0] * RADTODEG;
  swi_coortrf(xs, xs, xobl[1] * DEGTORAD);
  swi_cartpol(xs, xs);
  xs[0] *= RADTODEG;
  dhour = fmod(tjd_ut - 0.5, 1) * 360;
  /* mean to true (if nut != 0) */ 
  if (eps == 0)
    xs[0] += xobl[2] * cos(xobl[0] * DEGTORAD);
  else
    xs[0] += nut * cos(eps * DEGTORAD);
  /* add hour */
  xs[0] = swe_degnorm(xs[0] + dhour);
  tsid = xs[0] / 15;
  return tsid;
}
/*#endif*/

/* Apparent Sidereal Time at Greenwich with equation of the equinoxes