sideral_time.c   sideral_time.c 
skipping to change at line 43 skipping to change at line 43
/* puts a large angle in the correct range 0 - 2PI radians */ /* puts a large angle in the correct range 0 - 2PI radians */
double range_radians (double r) double range_radians (double r)
{ {
r = fmod( r, 2.*M_PI ); r = fmod( r, 2.*M_PI );
if (r<0.) r += 2.*M_PI; if (r<0.) r += 2.*M_PI;
return r; return r;
} }
#define TERMS 63 #define TERMS 63
/* compute only every tenths of day: */
#define LN_NUTATION_EPOCH_THRESHOLD 0.1 #define LN_NUTATION_EPOCH_THRESHOLD 0.1
/* Nutation is a period oscillation of the Earths rotational axis around it 's mean position.*/ /* Nutation is a period oscillation of the Earths rotational axis around it 's mean position.*/
/* /*
Contains Nutation in longitude, obliquity and ecliptic obliquity. Contains Nutation in longitude, obliquity and ecliptic obliquity.
Angles are expressed in degrees. Angles are expressed in degrees.
*/ */
struct ln_nutation struct ln_nutation
{ {
skipping to change at line 212 skipping to change at line 213
{-3.0, 0.0, 0.0, 0.0}, {-3.0, 0.0, 0.0, 0.0},
{-3.0, 0.0, 0.0, 0.0}, {-3.0, 0.0, 0.0, 0.0},
{-3.0, 0.0, 0.0, 0.0}, {-3.0, 0.0, 0.0, 0.0},
{-3.0, 0.0, 0.0, 0.0}}; {-3.0, 0.0, 0.0, 0.0}};
/* cache values */ /* cache values */
static double c_JD = 0.0, c_longitude = 0.0, c_obliquity = 0.0, c_ecliptic = 0.0; static double c_JD = 0.0, c_longitude = 0.0, c_obliquity = 0.0, c_ecliptic = 0.0;
/* Calculate nutation of longitude and obliquity in degrees from Julian Eph emeris Day /* Calculate nutation of longitude and obliquity in degrees from Julian Eph emeris Day
* params : JD Julian Day, nutation Pointer to store nutation. * params : JD Julian Day, nutation Pointer to store nutation.
* Chapter 21 pg 131-134 Using Table 21A */ * Meeus, Astr. Alg. (1st ed., 1994), Chapter 21 pg 131-134 Using Table 21A
*/
/* GZ: Changed: ecliptic obliquity used to be constant J2000.0.
* If you don't compute this, you may as well forget about nutation!
*/
void get_nutation (double JD, struct ln_nutation * nutation) void get_nutation (double JD, struct ln_nutation * nutation)
{ {
double D,M,MM,F,O,T,T2,T3; double D,M,MM,F,O,T;
double coeff_sine, coeff_cos; double coeff_sine, coeff_cos;
int i; int i;
/* should we bother recalculating nutation */ /* should we bother recalculating nutation */
if (fabs(JD - c_JD) > LN_NUTATION_EPOCH_THRESHOLD) if (fabs(JD - c_JD) > LN_NUTATION_EPOCH_THRESHOLD)
{ {
/* set the new epoch */ /* set the new epoch */
c_JD = JD; c_JD = JD;
/* set ecliptic */ /* set ecliptic. GZ: This is constant only, J2000.0. WRONG!
c_ecliptic = 23.0 + 26.0 / 60.0 + 27.407 / 3600.0; */
/* c_ecliptic = 23.0 + 26.0 / 60.0 + 27.407 / 3600.0; */
/* calc T */ /* calc T */
T = (JD - 2451545.0)/36525; T = (JD - 2451545.0)/36525;
/* GZotti: we don't need those. * /
T2 = T * T; T2 = T * T;
T3 = T2 * T; T3 = T2 * T;
/* calculate D,M,M',F and Omega */ / * calculate D,M,M',F and Omega * /
D = 297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 18 9474.0; D = 297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 18 9474.0;
M = 357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300 000.0; M = 357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300 000.0;
MM = 134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 5 6250.0; MM = 134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 5 6250.0;
F = 93.2719100 + 483202.017538 * T - 0.0036825 * T2 + T3 / 3 27270.0; F = 93.2719100 + 483202.017538 * T - 0.0036825 * T2 + T3 / 3 27270.0;
O = 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 4500 00.0; O = 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 4500 00.0;
/ * GZotti: Laskar's formula, only valid for J2000+/-10000 y
ears! */
if (fabs(T)<100)
{
double U=T/100.0;
c_ecliptic=((((((((((2.45*U+ 5.79)*U +27.87)*U +7.12)*U
-39.05)*U -249.67)*U
-51.38)*U +1999.25)*U -1.55)*U)-4680.93)
*U +21.448;
}
else /* Use IAU formula. This might be more stable in farawa
y times, but is less accurate in any case for |U|<1. */
{
c_ecliptic=((0.001813*T-0.00059)*T-46.8150)*T+21.448;
}
c_ecliptic/=60.0;
c_ecliptic+=26.0;
c_ecliptic/=60.0;
c_ecliptic+=23.0;
/* GZotti: I prefer Horner's Scheme and its better numerical
accuracy. */
D = ((T / 189474.0 - 0.0019142)*T + 445267.111480)*T + 297.8
5036;
M = ((T / 300000.0 - 0.0001603)*T + 35999.050340)*T + 357.5
2772;
MM= ((T / 56250.0 + 0.0086972)*T + 477198.867398)*T + 134.9
6298;
F = ((T / 327270.0 - 0.0036825)*T + 483202.017538)*T + 93.2
7191;
O = ((T / 450000.0 + 0.0020708)*T - 1934.136261)*T + 125.0
4452;
/* convert to radians */ /* convert to radians */
D *= M_PI/180.; D *= M_PI/180.;
M *= M_PI/180.; M *= M_PI/180.;
MM *= M_PI/180.; MM *= M_PI/180.;
F *= M_PI/180.; F *= M_PI/180.;
O *= M_PI/180.; O *= M_PI/180.;
/* calc sum of terms in table 21A */ /* calc sum of terms in table 21A */
for (i=0; i< TERMS; i++) for (i=0; i< TERMS; i++)
{ {
skipping to change at line 287 skipping to change at line 316
{ {
c_longitude += coeff_sine * (sin (arguments[ i].O * O)); c_longitude += coeff_sine * (sin (arguments[ i].O * O));
c_obliquity += coeff_cos * (cos (arguments[i ].O * O)); c_obliquity += coeff_cos * (cos (arguments[i ].O * O));
} }
} }
/* change to degrees */ /* change to degrees */
c_longitude /= 36000000.; c_longitude /= 36000000.;
c_obliquity /= 36000000.; c_obliquity /= 36000000.;
c_ecliptic += c_obliquity; /* GZ: mean ecliptic should be still available! Addition mus
} t be done where needed. */
/* c_ecliptic += c_obliquity; */
}
/* return results */ /* return results */
nutation->longitude = c_longitude; nutation->longitude = c_longitude;
nutation->obliquity = c_obliquity; nutation->obliquity = c_obliquity;
nutation->ecliptic = c_ecliptic; nutation->ecliptic = c_ecliptic;
} }
/* Calculate the mean sidereal time at the meridian of Greenwich of a given date. /* Calculate the mean sidereal time at the meridian of Greenwich of a given date.
* returns apparent sidereal time (degree). * returns apparent sidereal time (degree).
* Formula 11.1, 11.4 pg 83 */ * Formula 11.1, 11.4 pg 83 */
skipping to change at line 330 skipping to change at line 360
double correction, sidereal; double correction, sidereal;
struct ln_nutation nutation; struct ln_nutation nutation;
/* get the mean sidereal time */ /* get the mean sidereal time */
sidereal = get_mean_sidereal_time (JD); sidereal = get_mean_sidereal_time (JD);
/* add corrections for nutation in longitude and for the true obliquity of /* add corrections for nutation in longitude and for the true obliquity of
the ecliptic */ the ecliptic */
get_nutation (JD, &nutation); get_nutation (JD, &nutation);
correction = (nutation.longitude * cos (nutation.obliquity*M_PI/180.)); /* GZ: This was the only place where this was used. I added the summatio
n here. */
correction = (nutation.longitude * cos ((nutation.ecliptic+nutation.obli
quity)*M_PI/180.));
sidereal += correction; sidereal += correction;
return (sidereal); return (sidereal);
} }
double get_mean_ecliptical_obliquity(double JDE)
{
struct ln_nutation nutation;
get_nutation(JDE, &nutation);
return nutation.ecliptic;
}
 End of changes. 11 change blocks. 
9 lines changed or deleted 55 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/