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 oscillation of the rotational axis around /* Nutation is a periodic oscillation of the Earth's rotational axis around
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 Nutation in longitude double deltaPsi; //!< DeltaPsi, Nutation in longitude
double 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 */
= 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, 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 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 - > LN_NUTATION_EPOCH_THRESHOLD) if (fabs(JDE - c_JDE) > LN_NUTATION_EPOCH_THRESHOLD)
{ {
set the new epoch // set the new epoch
= 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 = - 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 (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 / 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;
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
for (i=0; i< TERMS; i++) for (i=0; i< TERMS; i++)
{ {
/ // 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
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
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
= c_longitude; nutation->deltaPsi = c_longitude;
= 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 sidereal time * 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 * 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 // double T, T1;
// T1 = (JD - 2451545.0);
= (JD - / // 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.;
/* /* add again a convenient multiple of 360 degrees */
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 * 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 * GZ modified for V0.14 to use nutation IAU-2000B
*/
double get_apparent_sidereal_time (double JD, double JDE)
{ {
double double meanSidereal = get_mean_sidereal_time (JD, JDE);
= get_mean_sidereal_time
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;
getNutationAngles(JDE, &deltaPsi, &deltaEps);
return meanSidereal+ (deltaPsi*cos(getPrecessionAngleVondrakEpsilon(
JDE) + deltaEps))*180./M_PI;
return
} }
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;
get_nutation_longitude(double JDE) //}
struct ln_nutation nutation; //double get_nutation_longitude(double JDE)
get_nutation(JDE, &nutation); //{
return // 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/