sidereal_time.c   sidereal_time.c 
skipping to change at line 17 skipping to change at line 17
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, US A. Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, US A.
Copyright (C) 2000 Liam Girdwood <liam@nova-ioe.org> Copyright (C) 2000 Liam Girdwood <liam@nova-ioe.org>
Copyright (C) 2015 Georg Zotti (deactivated old IAU-1980 Nutation functions , update to IAU-2000B)
*/ */
#include <math.h> #include <math.h>
#include <assert.h>
#include "precession.h"
#ifndef M_PI #ifndef M_PI
#define M_PI 3.14159265358979323846 #define M_PI 3.14159265358979323846
#endif #endif
/* puts a large angle in the correct range 0 - 360 degrees */ /* puts a large angle in the correct range 0 - 360 degrees */
double range_degrees(double d) double range_degrees(double d)
{ {
d = fmod( d, 360.); d = fmod( d, 360.);
if(d<0.) d += 360.; if(d<0.) d += 360.;
skipping to change at line 46 skipping to change at line 49
{ {
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: */ /* 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 /* Nutation is a periodic oscillation of the Earth's rotational axis around
's mean position.*/ its mean position.*/
/* GZ: These data from IAU1980 nutation are no longer in use: V0.14+ uses N
utation IAU-2000B, implemented in precession.c. */
/* /*
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.
*/ */
/* Structs no longer used
*
struct ln_nutation struct ln_nutation
{ {
double longitude; /*!< Nutation in longitude */ double deltaPsi; //!< DeltaPsi, Nutation in longitude
double obliquity; /*!< Nutation in obliquity */ double deltaEps; //!< DeltaEps, Nutation in obliquity
double ecliptic; /*!< Obliquity of the ecliptic */ double ecliptic; //!< epsilon, Mean Obliquity of the ecliptic
}; };
struct nutation_arguments struct nutation_arguments
{ {
double D; double D;
double M; double M;
double MM; double MM;
double F; double F;
double O; double O;
}; };
struct nutation_coefficients struct nutation_coefficients
{ {
double longitude1; double longitude1;
double longitude2; double longitude2;
double obliquity1; double obliquity1;
double obliquity2; double obliquity2;
}; };
*/
/* arguments and coefficients taken from table 21A on page 133 */ /* arguments and coefficients taken from table 21A on page 133 */
/* Table no longer used.
static const struct nutation_arguments arguments[TERMS] = { static const struct nutation_arguments arguments[TERMS] = {
{0.0, 0.0, 0.0, 0.0, 1.0}, {0.0, 0.0, 0.0, 0.0, 1.0},
{-2.0, 0.0, 0.0, 2.0, 2.0}, {-2.0, 0.0, 0.0, 2.0, 2.0},
{0.0, 0.0, 0.0, 2.0, 2.0}, {0.0, 0.0, 0.0, 2.0, 2.0},
{0.0, 0.0, 0.0, 0.0, 2.0}, {0.0, 0.0, 0.0, 0.0, 2.0},
{0.0, 1.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0, 0.0},
{0.0, 0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0, 0.0},
{-2.0, 1.0, 0.0, 2.0, 2.0}, {-2.0, 1.0, 0.0, 2.0, 2.0},
{0.0, 0.0, 0.0, 2.0, 1.0}, {0.0, 0.0, 0.0, 2.0, 1.0},
{0.0, 0.0, 1.0, 2.0, 2.0}, {0.0, 0.0, 1.0, 2.0, 2.0},
skipping to change at line 207 skipping to change at line 214
{-4.0, 0.0, 0.0, 0.0}, {-4.0, 0.0, 0.0, 0.0},
{-4.0, 0.0, 0.0, 0.0}, {-4.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}, {-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}}; {-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_JDE = 0.0, c_longitude = 0.0, c_obliquity = 0.0, c_eclip tic = 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 : JDE Julian Day (ET/TT), nutation Pointer to store nutation.
* Meeus, Astr. Alg. (1st ed., 1994), 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. /* GZ: Changed: ecliptic obliquity used to be constant J2000.0.
* If you don't compute this, you may as well forget about nutation! * If you don't compute this, you may as well forget about nutation!
* 2015: Laskar's 1986 formula replaced by Vondrak 2011. c_ecliptic is now
epsilon_A, Vondrak's obliquity of date.
* TODO: replace this whole function with Nutation IAU2000A or compatible v
ersion.
* FUNCTION IS UNUSED, MAYBE DELETE?
*/ */
void get_nutation (double JD, struct ln_nutation * nutation) /*
void get_nutation (double JDE, struct ln_nutation * nutation)
{ {
double D,M,MM,F,O,T; 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(JDE - c_JDE) > LN_NUTATION_EPOCH_THRESHOLD)
{ {
/* set the new epoch */ // set the new epoch
c_JD = JD; c_JDE = JDE;
/* set ecliptic. GZ: This is constant only, J2000.0. WRONG! // 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 = (JDE - 2451545.0)/36525;
/* GZotti: we don't need those. * / // 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 // D = 297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 18
9474.0; 9474.0;
M = 357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300 // M = 357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300
000.0; 000.0;
MM = 134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 5 // MM = 134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 5
6250.0; 6250.0;
F = 93.2719100 + 483202.017538 * T - 0.0036825 * T2 + T3 / 3 // F = 93.2719100 + 483202.017538 * T - 0.0036825 * T2 + T3 / 3
27270.0; 27270.0;
O = 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 4500 // O = 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 4500
00.0; 00.0;
// GZ: I prefer Horner's Scheme and its better numerical acc
uracy.
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;
/ * GZotti: Laskar's formula, only valid for J2000+/-10000 y ears! */ // GZotti: Laskar's formula (1986), only valid for J2000+/-1 0000 years!
if (fabs(T)<100) if (fabs(T)<100)
{ {
double U=T/100.0; 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 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; -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. */ else // Use IAU-1980 formula. This might be more stable in f araway 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=((0.001813*T-0.00059)*T-46.8150)*T+21.448;
} }
c_ecliptic/=60.0; c_ecliptic/=60.0;
c_ecliptic+=26.0; c_ecliptic+=26.0;
c_ecliptic/=60.0; c_ecliptic/=60.0;
c_ecliptic+=23.0; c_ecliptic+=23.0;
/* GZotti: I prefer Horner's Scheme and its better numerical // convert to radians
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 */
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++)
{ {
/* calc coefficients of sine and cosine */ // calc coefficients of sine and cosine
coeff_sine = coefficients[i].longitude1 + (coefficie nts[i].longitude2 * T); coeff_sine = coefficients[i].longitude1 + (coefficie nts[i].longitude2 * T);
coeff_cos = coefficients[i].obliquity1 + (coefficien ts[i].obliquity2 * T); coeff_cos = coefficients[i].obliquity1 + (coefficien ts[i].obliquity2 * T);
/* sum the arguments */ // sum the arguments
if (arguments[i].D != 0) if (arguments[i].D != 0)
{ {
c_longitude += coeff_sine * (sin (arguments[ i].D * D)); c_longitude += coeff_sine * (sin (arguments[ i].D * D));
c_obliquity += coeff_cos * (cos (arguments[i ].D * D)); c_obliquity += coeff_cos * (cos (arguments[i ].D * D));
} }
if (arguments[i].M != 0) if (arguments[i].M != 0)
{ {
c_longitude += coeff_sine * (sin (arguments[ i].M * M)); c_longitude += coeff_sine * (sin (arguments[ i].M * M));
c_obliquity += coeff_cos * (cos (arguments[i ].M * M)); c_obliquity += coeff_cos * (cos (arguments[i ].M * M));
} }
skipping to change at line 312 skipping to change at line 323
c_longitude += coeff_sine * (sin (arguments[ i].F * F)); c_longitude += coeff_sine * (sin (arguments[ i].F * F));
c_obliquity += coeff_cos * (cos (arguments[i ].F * F)); c_obliquity += coeff_cos * (cos (arguments[i ].F * F));
} }
if (arguments[i].O != 0) if (arguments[i].O != 0)
{ {
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.;
/* GZ: mean ecliptic should be still available! Addition mus // GZ: mean ecliptic should be still available! Addition mus
t be done where needed. */ t be done where needed.
/* c_ecliptic += c_obliquity; */ // c_ecliptic += c_obliquity;
} }
/* return results */ // return results
nutation->longitude = c_longitude; nutation->deltaPsi = c_longitude;
nutation->obliquity = c_obliquity; nutation->deltaEps = c_obliquity;
nutation->ecliptic = c_ecliptic; nutation->ecliptic = c_ecliptic;
} }
*/
/* Calculate the mean sidereal time at the meridian of Greenwich of a given /* Calculate the mean sidereal time at the meridian of Greenwich (GMST) of
date. a given date.
* returns apparent sidereal time (degree). * returns mean sidereal time (degrees).
* Formula 11.1, 11.4 pg 83 */ * Formula 11.1, 11.4 pg 83
double get_mean_sidereal_time (double JD) * MAKE SURE argument JD is UT, not TT!
* GZfix for V0.14: Replace by expression (43) given in:
* N. Capitaine, P.T.Wallace, J. Chapront: Expressions for IAU 2000 precess
ion quantities.
* A&A412, 567-586 (2003)
* DOI: 10.1051/0004-6361:20031539
*/
double get_mean_sidereal_time (double JD, double JDE)
{ {
double sidereal; // double sidereal;
double T; // double T, T1;
// T1 = (JD - 2451545.0);
T = (JD - 2451545.0) / 36525.0; // T= T1 * (1.0/ 36525.0);
//
// // calc mean angle. This is Meeus 11.4 (IAU 1982)
// sidereal = 280.46061837 + (360.98564736629 * (JD - 2451545.0)) + (0
.000387933 * T * T) - (T * T * T / 38710000.0);
// TODO GZ: verify that GMST1980 is still OK
//
double sidereal, UT1, t, tu;
UT1=(JD-floor(JD) +0.5) * 86400.; // time in seconds
t=(JDE-2451545.)/36525.;
tu=(JD-2451545.)/36525.;
sidereal= (((-0.000000002454*t-0.00000199708)*t-0.0000002926)*t+0.09
2772110)*t*t;
sidereal += (t-tu)*307.4771013;
sidereal += 8640184.79447825*tu + 24110.5493771;
sidereal += UT1;
// this is expressed in seconds. We need degrees.
// 1deg=4tempMin=240tempSec
sidereal *= 1./240.;
/* calc mean angle */ /* add again a convenient multiple of 360 degrees */
sidereal = 280.46061837 + (360.98564736629 * (JD - 2451545.0)) + (0.
000387933 * T * T) - (T * T * T / 38710000.0);
/* add a convenient multiple of 360 degrees */
sidereal = range_degrees (sidereal); sidereal = range_degrees (sidereal);
return sidereal; return sidereal;
} }
/* Calculate the apparent sidereal time at the meridian of Greenwich of a g iven date. /* Calculate the apparent sidereal time at the meridian of Greenwich of a g iven date.
* returns apparent sidereal time (degree). * returns apparent sidereal time (degrees).
* Formula 11.1, 11.4 pg 83 */ * Formula 11.1, 11.4 pg 83
double get_apparent_sidereal_time (double JD) * GZ modified for V0.14 to use nutation IAU-2000B
*/
double get_apparent_sidereal_time (double JD, double JDE)
{ {
double correction, sidereal; double meanSidereal = get_mean_sidereal_time (JD, JDE);
struct ln_nutation nutation;
/* get the mean sidereal time */
sidereal = get_mean_sidereal_time (JD);
/* add corrections for nutation in longitude and for the true obliqu // add corrections for nutation in longitude and for the true obliqu
ity of ity of the ecliptic
the ecliptic */ double deltaPsi, deltaEps;
get_nutation (JD, &nutation); getNutationAngles(JDE, &deltaPsi, &deltaEps);
/* GZ: This was the only place where this was used. I added the summ return meanSidereal+ (deltaPsi*cos(getPrecessionAngleVondrakEpsilon(
ation here. */ JDE) + deltaEps))*180./M_PI;
correction = (nutation.longitude * cos ((nutation.ecliptic+nutation.
obliquity)*M_PI/180.));
sidereal += correction;
return (sidereal);
} }
double get_mean_ecliptical_obliquity(double JDE) //// return value in degrees
{ //double get_mean_ecliptical_obliquity(double JDE)
struct ln_nutation nutation; //{
get_nutation(JDE, &nutation); //// struct ln_nutation nutation;
return nutation.ecliptic; //// get_nutation(JDE, &nutation);
} //// return nutation.ecliptic;
// return getPrecessionAngleVondrakEpsilon(JDE)*180.0/M_PI;
double get_nutation_longitude(double JDE) //}
{
struct ln_nutation nutation; //double get_nutation_longitude(double JDE)
get_nutation(JDE, &nutation); //{
return nutation.longitude; // struct ln_nutation nutation;
} // get_nutation(JDE, &nutation);
// return nutation.deltaPsi;
//}
 End of changes. 35 change blocks. 
90 lines changed or deleted 121 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/