This is nutate.c in view mode; [Download] [Up]
/************************************************************************/
/* nutate.c */
/* */
/* Nutation in longitude and obliquity */
/* computed at Julian date J. */
/* */
/* References: */
/* "Summary of 1980 IAU Theory of Nutation (Final Report of the */
/* IAU Working Group on Nutation)", P. K. Seidelmann et al., in */
/* Transactions of the IAU Vol. XVIII A, Reports on Astronomy, */
/* P. A. Wayman, ed.; D. Reidel Pub. Co., 1982. */
/* */
/* "Nutation and the Earth's Rotation", */
/* I.A.U. Symposium No. 78, May, 1977, page 256. */
/* I.A.U., 1980. */
/* */
/* Woolard, E.W., "A redevelopment of the theory of nutation", */
/* The Astronomical Journal, 58, 1-3 (1953). */
/* */
/* This program implements all of the 1980 IAU nutation series. */
/* Results checked at 100 points against the 1986 AA; all agreed. */
/* */
/* */
/* - S. L. Moshier, November 1987 */
/* October, 1992 - typo fixed in nutation matrix */
/************************************************************************/
/***** description
*
* $Id: nutate.c,v 1.3 1993/04/21 21:26:48 craig Exp $
*
*/
/***** modification history
*
* $Log: nutate.c,v $
* Revision 1.3 1993/04/21 21:26:48 craig
* Changed the path of the satellite.h include.
* Changed ecnsts to pcnsts.
*
* Revision 1.2 1993/04/21 15:22:07 craig
* First working version. Ran through indent and converted to ansi.
* Added hooks for working with the satellite programs.
*
*
*/
/***** include files *****/
#include <math.h>
#include "aaproto.h"
#include "satellite.h"
/***** global variables *****/
extern struct MCONSTANTS mcnsts;
extern struct PCONSTANTS pcnsts;
/* from epsiln.c */
extern double eps, coseps, sineps;
/* The answers are posted here by nutlo(): */
double jdnut = -1.0; /* time to which the nutation
applies */
double nutl = 0.0; /* nutation in longitude (radians) */
double nuto = 0.0; /* nutation in obliquity (radians) */
/* arrays to hold sines and cosines of multiple angles */
double ss[5][8];
double cc[5][8];
/***** local global variables *****/
/* Each term in the expansion has a trigonometric
* argument given by
* W = i*MM + j*MS + k*FF + l*DD + m*OM
* where the variables are defined below.
* The nutation in longitude is a sum of terms of the
* form (a + bT) * sin(W). The terms for nutation in obliquity
* are of the form (c + dT) * cos(W). The coefficients
* are arranged in the tabulation as follows:
*
* Coefficient:
* i j k l m a b c d
* 0, 0, 0, 0, 1, -171996, -1742, 92025, 89,
* The first line of the table, above, is done separately
* since two of the values do not fit into 16 bit integers.
* The values a and c are arc seconds times 10000. b and d
* are arc seconds per Julian century times 100000. i through m
* are integers. See the program for interpretation of MM, MS,
* etc., which are mean orbital elements of the Sun and Moon.
*
* If terms with coefficient less than X are omitted, the peak
* errors will be:
*
* omit error, omit error,
* a < longitude c < obliquity
* .0005" .0100" .0008" .0094"
* .0046 .0492 .0095 .0481
* .0123 .0880 .0224 .0905
* .0386 .1808 .0895 .1129
*/
static short nt[105 * 9] = {
0, 0, 0, 0, 2, 2062, 2, -895, 5,
-2, 0, 2, 0, 1, 46, 0, -24, 0,
2, 0, -2, 0, 0, 11, 0, 0, 0,
-2, 0, 2, 0, 2, -3, 0, 1, 0,
1, -1, 0, -1, 0, -3, 0, 0, 0,
0, -2, 2, -2, 1, -2, 0, 1, 0,
2, 0, -2, 0, 1, 1, 0, 0, 0,
0, 0, 2, -2, 2, -13187, -16, 5736, -31,
0, 1, 0, 0, 0, 1426, -34, 54, -1,
0, 1, 2, -2, 2, -517, 12, 224, -6,
0, -1, 2, -2, 2, 217, -5, -95, 3,
0, 0, 2, -2, 1, 129, 1, -70, 0,
2, 0, 0, -2, 0, 48, 0, 1, 0,
0, 0, 2, -2, 0, -22, 0, 0, 0,
0, 2, 0, 0, 0, 17, -1, 0, 0,
0, 1, 0, 0, 1, -15, 0, 9, 0,
0, 2, 2, -2, 2, -16, 1, 7, 0,
0, -1, 0, 0, 1, -12, 0, 6, 0,
-2, 0, 0, 2, 1, -6, 0, 3, 0,
0, -1, 2, -2, 1, -5, 0, 3, 0,
2, 0, 0, -2, 1, 4, 0, -2, 0,
0, 1, 2, -2, 1, 4, 0, -2, 0,
1, 0, 0, -1, 0, -4, 0, 0, 0,
2, 1, 0, -2, 0, 1, 0, 0, 0,
0, 0, -2, 2, 1, 1, 0, 0, 0,
0, 1, -2, 2, 0, -1, 0, 0, 0,
0, 1, 0, 0, 2, 1, 0, 0, 0,
-1, 0, 0, 1, 1, 1, 0, 0, 0,
0, 1, 2, -2, 0, -1, 0, 0, 0,
0, 0, 2, 0, 2, -2274, -2, 977, -5,
1, 0, 0, 0, 0, 712, 1, -7, 0,
0, 0, 2, 0, 1, -386, -4, 200, 0,
1, 0, 2, 0, 2, -301, 0, 129, -1,
1, 0, 0, -2, 0, -158, 0, -1, 0,
-1, 0, 2, 0, 2, 123, 0, -53, 0,
0, 0, 0, 2, 0, 63, 0, -2, 0,
1, 0, 0, 0, 1, 63, 1, -33, 0,
-1, 0, 0, 0, 1, -58, -1, 32, 0,
-1, 0, 2, 2, 2, -59, 0, 26, 0,
1, 0, 2, 0, 1, -51, 0, 27, 0,
0, 0, 2, 2, 2, -38, 0, 16, 0,
2, 0, 0, 0, 0, 29, 0, -1, 0,
1, 0, 2, -2, 2, 29, 0, -12, 0,
2, 0, 2, 0, 2, -31, 0, 13, 0,
0, 0, 2, 0, 0, 26, 0, -1, 0,
-1, 0, 2, 0, 1, 21, 0, -10, 0,
-1, 0, 0, 2, 1, 16, 0, -8, 0,
1, 0, 0, -2, 1, -13, 0, 7, 0,
-1, 0, 2, 2, 1, -10, 0, 5, 0,
1, 1, 0, -2, 0, -7, 0, 0, 0,
0, 1, 2, 0, 2, 7, 0, -3, 0,
0, -1, 2, 0, 2, -7, 0, 3, 0,
1, 0, 2, 2, 2, -8, 0, 3, 0,
1, 0, 0, 2, 0, 6, 0, 0, 0,
2, 0, 2, -2, 2, 6, 0, -3, 0,
0, 0, 0, 2, 1, -6, 0, 3, 0,
0, 0, 2, 2, 1, -7, 0, 3, 0,
1, 0, 2, -2, 1, 6, 0, -3, 0,
0, 0, 0, -2, 1, -5, 0, 3, 0,
1, -1, 0, 0, 0, 5, 0, 0, 0,
2, 0, 2, 0, 1, -5, 0, 3, 0,
0, 1, 0, -2, 0, -4, 0, 0, 0,
1, 0, -2, 0, 0, 4, 0, 0, 0,
0, 0, 0, 1, 0, -4, 0, 0, 0,
1, 1, 0, 0, 0, -3, 0, 0, 0,
1, 0, 2, 0, 0, 3, 0, 0, 0,
1, -1, 2, 0, 2, -3, 0, 1, 0,
-1, -1, 2, 2, 2, -3, 0, 1, 0,
-2, 0, 0, 0, 1, -2, 0, 1, 0,
3, 0, 2, 0, 2, -3, 0, 1, 0,
0, -1, 2, 2, 2, -3, 0, 1, 0,
1, 1, 2, 0, 2, 2, 0, -1, 0,
-1, 0, 2, -2, 1, -2, 0, 1, 0,
2, 0, 0, 0, 1, 2, 0, -1, 0,
1, 0, 0, 0, 2, -2, 0, 1, 0,
3, 0, 0, 0, 0, 2, 0, 0, 0,
0, 0, 2, 1, 2, 2, 0, -1, 0,
-1, 0, 0, 0, 2, 1, 0, -1, 0,
1, 0, 0, -4, 0, -1, 0, 0, 0,
-2, 0, 2, 2, 2, 1, 0, -1, 0,
-1, 0, 2, 4, 2, -2, 0, 1, 0,
2, 0, 0, -4, 0, -1, 0, 0, 0,
1, 1, 2, -2, 2, 1, 0, -1, 0,
1, 0, 2, 2, 1, -1, 0, 1, 0,
-2, 0, 2, 4, 2, -1, 0, 1, 0,
-1, 0, 4, 0, 2, 1, 0, 0, 0,
1, -1, 0, -2, 0, 1, 0, 0, 0,
2, 0, 2, -2, 1, 1, 0, -1, 0,
2, 0, 2, 2, 2, -1, 0, 0, 0,
1, 0, 0, 2, 1, -1, 0, 0, 0,
0, 0, 4, -2, 2, 1, 0, 0, 0,
3, 0, 2, -2, 2, 1, 0, 0, 0,
1, 0, 2, -2, 0, -1, 0, 0, 0,
0, 1, 2, 0, 1, 1, 0, 0, 0,
-1, -1, 0, 2, 1, 1, 0, 0, 0,
0, 0, -2, 0, 1, -1, 0, 0, 0,
0, 0, 2, -1, 2, -1, 0, 0, 0,
0, 1, 0, 2, 0, -1, 0, 0, 0,
1, 0, -2, -2, 0, -1, 0, 0, 0,
0, -1, 2, 0, 1, -1, 0, 0, 0,
1, 1, 0, -2, 1, -1, 0, 0, 0,
1, 0, -2, 2, 0, -1, 0, 0, 0,
2, 0, 0, 2, 0, 1, 0, 0, 0,
0, 0, 2, 4, 2, -1, 0, 0, 0,
0, 1, 0, 1, 0, 1, 0, 0, 0,
};
/*********/
/* nutlo */
/*********/
int nutlo (double J)
{
int i, j, k, k1, m;
short *p;
double f, g, T;
double MM, MS, FF, DD, OM;
double cu, su, cv, sv, sw;
double C, D;
if (jdnut == J)
{
return (0);
}
jdnut = J;
/* Julian centuries from 2000 January 1.5, barycentric dynamical time */
T = (J - pcnsts.J2000) / pcnsts.dapcen;
/* Fundamental arguments in the FK5 reference system. The
coefficients, originally given to 0.001", are converted here to
degrees. */
/* longitude of the mean ascending node of the lunar orbit on the
ecliptic, measured from the mean equinox of date */
OM = ((2.22e-6 * T + 0.00207833) * T - 1934.1362608) * T + 125.0445222;
OM = mcnsts.de2ra * mod360 (OM);
/* mean longitude of the Sun minus the mean longitude of the Sun's
perigee */
MS = ((-3.33e-6 * T - 0.0001603) * T + 35999.0503400) * T + 357.5277233;
MS = mcnsts.de2ra * mod360 (MS);
/* mean longitude of the Moon minus the mean longitude of the Moon's
perigee */
MM = ((1.778e-5 * T + 0.0086972) * T + 477198.8673981) * T + 134.9629814;
MM = mcnsts.de2ra * mod360 (MM);
/* mean longitude of the Moon minus the mean longitude of the Moon's
node */
FF = ((3.056e-6 * T - 0.00368250) * T + 483202.0175381) * T + 93.2719103;
FF = mcnsts.de2ra * mod360 (FF);
/* mean elongation of the Moon from the Sun. */
DD = ((5.278e-6 * T - 0.0019142) * T + 445267.1114800) * T + 297.8503631;
DD = mcnsts.de2ra * mod360 (DD);
/* Calculate sin( i*MM ), etc. for needed multiple angles */
sscc (0, MM, 3);
sscc (1, MS, 2);
sscc (2, FF, 4);
sscc (3, DD, 4);
sscc (4, OM, 2);
/* first terms, not in table: */
C = (-0.01742 * T - 17.1996) * ss[4][0]; /* sin(OM) */
D = (0.00089 * T + 9.2025) * cc[4][0]; /* cos(OM) */
p = &nt[0]; /* point to start of table */
for (i = 0; i < 105; i++)
{
/* argument of sine and cosine */
k1 = 0;
cv = 0.0;
sv = 0.0;
for (m = 0; m < 5; m++)
{
j = *p++;
if (j)
{
k = j;
if (j < 0)
{
k = -k;
}
su = ss[m][k - 1]; /* sin(k*angle) */
if (j < 0)
{
su = -su;
}
cu = cc[m][k - 1];
if (k1 == 0)
{ /* set first angle */
sv = su;
cv = cu;
k1 = 1;
}
else
{ /* combine angles */
sw = su * cv + cu * sv;
cv = cu * cv - su * sv;
sv = sw;
}
}
}
/* longitude coefficient */
f = *p++ * 0.0001;
if ((k = *p++) != 0)
{
f += 0.00001 * T * k;
}
/* obliquity coefficient */
g = *p++ * 0.0001;
if ((k = *p++) != 0)
{
g += 0.00001 * T * k;
}
/* accumulate the terms */
C += f * sv;
D += g * cv;
}
/*
printf( "nutation: in longitude %.3f\", in obliquity %.3f\"\n", C, D );
*/
/* Save answers, expressed in radians */
nutl = mcnsts.de2ra * C / 3600.0;
nuto = mcnsts.de2ra * D / 3600.0;
return (0);
}
/****************************************************************/
/* nutate */
/* */
/* Nutation -- AA page B20 */
/* using nutation in longitude and obliquity from nutlo() */
/* and obliquity of the ecliptic from epsiln() */
/* both calculated for Julian date J. */
/* */
/* p[] = equatorial rectangular position vector of object for */
/* mean ecliptic and equinox of date. */
/****************************************************************/
int nutate (double J, double *p)
{
int i;
double ce, se, cl, sl, sino, f;
double dp[3], p1[3];
nutlo (J); /* be sure we calculated nutl and
nuto */
epsiln (J); /* and also the obliquity of date */
f = eps + nuto;
ce = cos (f);
se = sin (f);
sino = sin (nuto);
cl = cos (nutl);
sl = sin (nutl);
/* Apply adjustment to equatorial rectangular coordinates of object.
This is a composite of three rotations: rotate about x axis to
ecliptic of date; rotate about new z axis by the nutation in
longitude; rotate about new x axis back to equator of date plus
nutation in obliquity. */
p1[0] = cl * p[0]
- sl * coseps * p[1]
- sl * sineps * p[2];
p1[1] = sl * ce * p[0]
+ (cl * coseps * ce + sineps * se) * p[1]
- (sino + (1.0 - cl) * sineps * ce) * p[2];
p1[2] = sl * se * p[0]
+ (sino + (cl - 1.0) * se * coseps) * p[1]
+ (cl * sineps * se + coseps * ce) * p[2];
for (i = 0; i < 3; i++)
{
dp[i] = p1[i] - p[i];
}
showcor ("nutation", p, dp);
for (i = 0; i < 3; i++)
{
p[i] = p1[i];
}
return (0);
}
/************************************************/
/* sscc */
/* */
/* Prepare lookup table of sin and cos ( i*Lj ) */
/* for required multiple angles */
/************************************************/
int sscc (int k, double arg, int n)
{
int i;
double cu, su, cv, sv, s;
su = sin (arg);
cu = cos (arg);
ss[k][0] = su; /* sin(L) */
cc[k][0] = cu; /* cos(L) */
sv = 2.0 * su * cu;
cv = cu * cu - su * su;
ss[k][1] = sv; /* sin(2L) */
cc[k][1] = cv;
for (i = 2; i < n; i++)
{
s = su * cv + cu * sv;
cv = cu * cv - su * sv;
sv = s;
ss[k][i] = sv; /* sin( i+1 L ) */
cc[k][i] = cv;
}
return (0);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.