4626 lines
160 KiB
C
4626 lines
160 KiB
C
/* SWISSEPH
|
|
|
|
SWISSEPH modules that may be useful for other applications
|
|
e.g. chopt.c, venus.c, swetest.c
|
|
|
|
Authors: Dieter Koch and Alois Treindl, Astrodienst Zurich
|
|
|
|
coordinate transformations
|
|
obliquity of ecliptic
|
|
nutation
|
|
precession
|
|
delta t
|
|
sidereal time
|
|
setting or getting of tidal acceleration of moon
|
|
chebyshew interpolation
|
|
ephemeris file name generation
|
|
cyclic redundancy checksum CRC
|
|
modulo and normalization functions
|
|
|
|
**************************************************************/
|
|
/* Copyright (C) 1997 - 2008 Astrodienst AG, Switzerland. All rights reserved.
|
|
|
|
License conditions
|
|
------------------
|
|
|
|
This file is part of Swiss Ephemeris.
|
|
|
|
Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
|
|
or distributor accepts any responsibility for the consequences of using it,
|
|
or for whether it serves any particular purpose or works at all, unless he
|
|
or she says so in writing.
|
|
|
|
Swiss Ephemeris is made available by its authors under a dual licensing
|
|
system. The software developer, who uses any part of Swiss Ephemeris
|
|
in his or her software, must choose between one of the two license models,
|
|
which are
|
|
a) GNU public license version 2 or later
|
|
b) Swiss Ephemeris Professional License
|
|
|
|
The choice must be made before the software developer distributes software
|
|
containing parts of Swiss Ephemeris to others, and before any public
|
|
service using the developed software is activated.
|
|
|
|
If the developer choses the GNU GPL software license, he or she must fulfill
|
|
the conditions of that license, which includes the obligation to place his
|
|
or her whole software project under the GNU GPL or a compatible license.
|
|
See http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
|
|
|
|
If the developer choses the Swiss Ephemeris Professional license,
|
|
he must follow the instructions as found in http://www.astro.com/swisseph/
|
|
and purchase the Swiss Ephemeris Professional Edition from Astrodienst
|
|
and sign the corresponding license contract.
|
|
|
|
The License grants you the right to use, copy, modify and redistribute
|
|
Swiss Ephemeris, but only under certain conditions described in the License.
|
|
Among other things, the License requires that the copyright notices and
|
|
this notice be preserved on all copies.
|
|
|
|
Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
|
|
|
|
The authors of Swiss Ephemeris have no control or influence over any of
|
|
the derived works, i.e. over software or services created by other
|
|
programmers which use Swiss Ephemeris functions.
|
|
|
|
The names of the authors or of the copyright holder (Astrodienst) must not
|
|
be used for promoting any software, product or service which uses or contains
|
|
the Swiss Ephemeris. This copyright notice is the ONLY place where the
|
|
names of the authors can legally appear, except in cases where they have
|
|
given special permission in writing.
|
|
|
|
The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
|
|
for promoting such software, products or services.
|
|
*/
|
|
|
|
#include <string.h>
|
|
#include <ctype.h>
|
|
#include "swephexp.h"
|
|
#include "sweph.h"
|
|
#include "swephlib.h"
|
|
#if MSDOS
|
|
# include <process.h>
|
|
# define strdup _strdup
|
|
#endif
|
|
|
|
#ifdef TRACE
|
|
void swi_open_trace(char *serr);
|
|
TLS FILE *swi_fp_trace_c = NULL;
|
|
TLS FILE *swi_fp_trace_out = NULL;
|
|
TLS int32 swi_trace_count = 0;
|
|
#endif
|
|
|
|
static void init_crc32(void);
|
|
static int init_dt(void);
|
|
static double adjust_for_tidacc(double ans, double Y, double tid_acc, double tid_acc0, AS_BOOL adjust_after_1955);
|
|
static double deltat_espenak_meeus_1620(double tjd, double tid_acc);
|
|
static double deltat_stephenson_etc_2016(double tjd, double tid_acc);
|
|
static double deltat_longterm_morrison_stephenson(double tjd);
|
|
static double deltat_stephenson_morrison_2004_1600(double tjd, double tid_acc);
|
|
static double deltat_stephenson_morrison_1997_1600(double tjd, double tid_acc);
|
|
static double deltat_aa(double tjd, double tid_acc);
|
|
|
|
#define SEFLG_EPHMASK (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH)
|
|
|
|
/* Reduce x modulo 360 degrees
|
|
*/
|
|
double CALL_CONV swe_degnorm(double x)
|
|
{
|
|
double y;
|
|
y = fmod(x, 360.0);
|
|
if (fabs(y) < 1e-13) y = 0; /* Alois fix 11-dec-1999 */
|
|
if( y < 0.0 ) y += 360.0;
|
|
return(y);
|
|
}
|
|
|
|
/* Reduce x modulo TWOPI degrees
|
|
*/
|
|
double CALL_CONV swe_radnorm(double x)
|
|
{
|
|
double y;
|
|
y = fmod(x, TWOPI);
|
|
if (fabs(y) < 1e-13) y = 0; /* Alois fix 11-dec-1999 */
|
|
if( y < 0.0 ) y += TWOPI;
|
|
return(y);
|
|
}
|
|
|
|
double CALL_CONV swe_deg_midp(double x1, double x0)
|
|
{
|
|
double d, y;
|
|
d = swe_difdeg2n(x1, x0); /* arc from x0 to x1 */
|
|
y = swe_degnorm(x0 + d / 2);
|
|
return(y);
|
|
}
|
|
|
|
double CALL_CONV swe_rad_midp(double x1, double x0)
|
|
{
|
|
return DEGTORAD * swe_deg_midp(x1 * RADTODEG, x0 * RADTODEG);
|
|
}
|
|
|
|
/* Reduce x modulo 2*PI
|
|
*/
|
|
double swi_mod2PI(double x)
|
|
{
|
|
double y;
|
|
y = fmod(x, TWOPI);
|
|
if( y < 0.0 ) y += TWOPI;
|
|
return(y);
|
|
}
|
|
|
|
|
|
double swi_angnorm(double x)
|
|
{
|
|
if (x < 0.0 )
|
|
return x + TWOPI;
|
|
else if (x >= TWOPI)
|
|
return x - TWOPI;
|
|
else
|
|
return x;
|
|
}
|
|
|
|
void swi_cross_prod(double *a, double *b, double *x)
|
|
{
|
|
x[0] = a[1]*b[2] - a[2]*b[1];
|
|
x[1] = a[2]*b[0] - a[0]*b[2];
|
|
x[2] = a[0]*b[1] - a[1]*b[0];
|
|
}
|
|
|
|
/* Evaluates a given chebyshev series coef[0..ncf-1]
|
|
* with ncf terms at x in [-1,1]. Communications of the ACM, algorithm 446,
|
|
* April 1973 (vol. 16 no.4) by Dr. Roger Broucke.
|
|
*/
|
|
double swi_echeb(double x, double *coef, int ncf)
|
|
{
|
|
int j;
|
|
double x2, br, brp2, brpp;
|
|
x2 = x * 2.;
|
|
br = 0.;
|
|
brp2 = 0.; /* dummy assign to silence gcc warning */
|
|
brpp = 0.;
|
|
for (j = ncf - 1; j >= 0; j--) {
|
|
brp2 = brpp;
|
|
brpp = br;
|
|
br = x2 * brpp - brp2 + coef[j];
|
|
}
|
|
return (br - brp2) * .5;
|
|
}
|
|
|
|
/*
|
|
* evaluates derivative of chebyshev series, see echeb
|
|
*/
|
|
double swi_edcheb(double x, double *coef, int ncf)
|
|
{
|
|
double bjpl, xjpl;
|
|
int j;
|
|
double x2, bf, bj, dj, xj, bjp2, xjp2;
|
|
x2 = x * 2.;
|
|
bf = 0.; /* dummy assign to silence gcc warning */
|
|
bj = 0.; /* dummy assign to silence gcc warning */
|
|
xjp2 = 0.;
|
|
xjpl = 0.;
|
|
bjp2 = 0.;
|
|
bjpl = 0.;
|
|
for (j = ncf - 1; j >= 1; j--) {
|
|
dj = (double) (j + j);
|
|
xj = coef[j] * dj + xjp2;
|
|
bj = x2 * bjpl - bjp2 + xj;
|
|
bf = bjp2;
|
|
bjp2 = bjpl;
|
|
bjpl = bj;
|
|
xjp2 = xjpl;
|
|
xjpl = xj;
|
|
}
|
|
return (bj - bf) * .5;
|
|
}
|
|
|
|
/*
|
|
* conversion between ecliptical and equatorial polar coordinates.
|
|
* for users of SWISSEPH, not used by our routines.
|
|
* for ecl. to equ. eps must be negative.
|
|
* for equ. to ecl. eps must be positive.
|
|
* xpo, xpn are arrays of 3 doubles containing position.
|
|
* attention: input must be in degrees!
|
|
*/
|
|
void CALL_CONV swe_cotrans(double *xpo, double *xpn, double eps)
|
|
{
|
|
int i;
|
|
double x[6], e = eps * DEGTORAD;
|
|
for(i = 0; i <= 1; i++)
|
|
x[i] = xpo[i];
|
|
x[0] *= DEGTORAD;
|
|
x[1] *= DEGTORAD;
|
|
x[2] = 1;
|
|
for(i = 3; i <= 5; i++)
|
|
x[i] = 0;
|
|
swi_polcart(x, x);
|
|
swi_coortrf(x, x, e);
|
|
swi_cartpol(x, x);
|
|
xpn[0] = x[0] * RADTODEG;
|
|
xpn[1] = x[1] * RADTODEG;
|
|
xpn[2] = xpo[2];
|
|
}
|
|
|
|
/*
|
|
* conversion between ecliptical and equatorial polar coordinates
|
|
* with speed.
|
|
* for users of SWISSEPH, not used by our routines.
|
|
* for ecl. to equ. eps must be negative.
|
|
* for equ. to ecl. eps must be positive.
|
|
* xpo, xpn are arrays of 6 doubles containing position and speed.
|
|
* attention: input must be in degrees!
|
|
*/
|
|
void CALL_CONV swe_cotrans_sp(double *xpo, double *xpn, double eps)
|
|
{
|
|
int i;
|
|
double x[6], e = eps * DEGTORAD;
|
|
for (i = 0; i <= 5; i++)
|
|
x[i] = xpo[i];
|
|
x[0] *= DEGTORAD;
|
|
x[1] *= DEGTORAD;
|
|
x[2] = 1; /* avoids problems with polcart(), if x[2] = 0 */
|
|
x[3] *= DEGTORAD;
|
|
x[4] *= DEGTORAD;
|
|
swi_polcart_sp(x, x);
|
|
swi_coortrf(x, x, e);
|
|
swi_coortrf(x+3, x+3, e);
|
|
swi_cartpol_sp(x, xpn);
|
|
xpn[0] *= RADTODEG;
|
|
xpn[1] *= RADTODEG;
|
|
xpn[2] = xpo[2];
|
|
xpn[3] *= RADTODEG;
|
|
xpn[4] *= RADTODEG;
|
|
xpn[5] = xpo[5];
|
|
}
|
|
|
|
/*
|
|
* conversion between ecliptical and equatorial cartesian coordinates
|
|
* for ecl. to equ. eps must be negative
|
|
* for equ. to ecl. eps must be positive
|
|
*/
|
|
void swi_coortrf(double *xpo, double *xpn, double eps)
|
|
{
|
|
double sineps, coseps;
|
|
double x[3];
|
|
sineps = sin(eps);
|
|
coseps = cos(eps);
|
|
x[0] = xpo[0];
|
|
x[1] = xpo[1] * coseps + xpo[2] * sineps;
|
|
x[2] = -xpo[1] * sineps + xpo[2] * coseps;
|
|
xpn[0] = x[0];
|
|
xpn[1] = x[1];
|
|
xpn[2] = x[2];
|
|
}
|
|
|
|
/*
|
|
* conversion between ecliptical and equatorial cartesian coordinates
|
|
* sineps sin(eps)
|
|
* coseps cos(eps)
|
|
* for ecl. to equ. sineps must be -sin(eps)
|
|
*/
|
|
void swi_coortrf2(double *xpo, double *xpn, double sineps, double coseps)
|
|
{
|
|
double x[3];
|
|
x[0] = xpo[0];
|
|
x[1] = xpo[1] * coseps + xpo[2] * sineps;
|
|
x[2] = -xpo[1] * sineps + xpo[2] * coseps;
|
|
xpn[0] = x[0];
|
|
xpn[1] = x[1];
|
|
xpn[2] = x[2];
|
|
}
|
|
|
|
/* conversion of cartesian (x[3]) to polar coordinates (l[3]).
|
|
* x = l is allowed.
|
|
* if |x| = 0, then lon, lat and rad := 0.
|
|
*/
|
|
void swi_cartpol(double *x, double *l)
|
|
{
|
|
double rxy;
|
|
double ll[3];
|
|
if (x[0] == 0 && x[1] == 0 && x[2] == 0) {
|
|
l[0] = l[1] = l[2] = 0;
|
|
return;
|
|
}
|
|
rxy = x[0]*x[0] + x[1]*x[1];
|
|
ll[2] = sqrt(rxy + x[2]*x[2]);
|
|
rxy = sqrt(rxy);
|
|
ll[0] = atan2(x[1], x[0]);
|
|
if (ll[0] < 0.0) ll[0] += TWOPI;
|
|
if (rxy == 0) {
|
|
if (x[2] >= 0)
|
|
ll[1] = PI / 2;
|
|
else
|
|
ll[1] = -(PI / 2);
|
|
} else {
|
|
ll[1] = atan(x[2] / rxy);
|
|
}
|
|
l[0] = ll[0];
|
|
l[1] = ll[1];
|
|
l[2] = ll[2];
|
|
}
|
|
|
|
/* conversion from polar (l[3]) to cartesian coordinates (x[3]).
|
|
* x = l is allowed.
|
|
*/
|
|
void swi_polcart(double *l, double *x)
|
|
{
|
|
double xx[3];
|
|
double cosl1;
|
|
cosl1 = cos(l[1]);
|
|
xx[0] = l[2] * cosl1 * cos(l[0]);
|
|
xx[1] = l[2] * cosl1 * sin(l[0]);
|
|
xx[2] = l[2] * sin(l[1]);
|
|
x[0] = xx[0];
|
|
x[1] = xx[1];
|
|
x[2] = xx[2];
|
|
}
|
|
|
|
/* conversion of position and speed.
|
|
* from cartesian (x[6]) to polar coordinates (l[6]).
|
|
* x = l is allowed.
|
|
* if position is 0, function returns direction of
|
|
* motion.
|
|
*/
|
|
void swi_cartpol_sp(double *x, double *l)
|
|
{
|
|
int i;
|
|
double xx[6], ll[6];
|
|
double rxy, coslon, sinlon, coslat, sinlat;
|
|
/* zero position */
|
|
if (x[0] == 0 && x[1] == 0 && x[2] == 0) {
|
|
ll[0] = ll[1] = ll[3] = ll[4] = 0;
|
|
ll[5] = sqrt(square_sum((x+3)));
|
|
swi_cartpol(x+3, ll);
|
|
ll[2] = 0;
|
|
for (i = 0; i <= 5; i++)
|
|
l[i] = ll[i];
|
|
return;
|
|
}
|
|
/* zero speed */
|
|
if (x[3] == 0 && x[4] == 0 && x[5] == 0) {
|
|
l[3] = l[4] = l[5] = 0;
|
|
swi_cartpol(x, l);
|
|
return;
|
|
}
|
|
/* position */
|
|
rxy = x[0]*x[0] + x[1]*x[1];
|
|
ll[2] = sqrt(rxy + x[2]*x[2]);
|
|
rxy = sqrt(rxy);
|
|
ll[0] = atan2(x[1], x[0]);
|
|
if (ll[0] < 0.0) ll[0] += TWOPI;
|
|
ll[1] = atan(x[2] / rxy);
|
|
/* speed:
|
|
* 1. rotate coordinate system by longitude of position about z-axis,
|
|
* so that new x-axis = position radius projected onto x-y-plane.
|
|
* in the new coordinate system
|
|
* vy'/r = dlong/dt, where r = sqrt(x^2 +y^2).
|
|
* 2. rotate coordinate system by latitude about new y-axis.
|
|
* vz"/r = dlat/dt, where r = position radius.
|
|
* vx" = dr/dt
|
|
*/
|
|
coslon = x[0] / rxy; /* cos(l[0]); */
|
|
sinlon = x[1] / rxy; /* sin(l[0]); */
|
|
coslat = rxy / ll[2]; /* cos(l[1]); */
|
|
sinlat = x[2] / ll[2]; /* sin(ll[1]); */
|
|
xx[3] = x[3] * coslon + x[4] * sinlon;
|
|
xx[4] = -x[3] * sinlon + x[4] * coslon;
|
|
l[3] = xx[4] / rxy; /* speed in longitude */
|
|
xx[4] = -sinlat * xx[3] + coslat * x[5];
|
|
xx[5] = coslat * xx[3] + sinlat * x[5];
|
|
l[4] = xx[4] / ll[2]; /* speed in latitude */
|
|
l[5] = xx[5]; /* speed in radius */
|
|
l[0] = ll[0]; /* return position */
|
|
l[1] = ll[1];
|
|
l[2] = ll[2];
|
|
}
|
|
|
|
/* conversion of position and speed
|
|
* from polar (l[6]) to cartesian coordinates (x[6])
|
|
* x = l is allowed
|
|
* explanation s. swi_cartpol_sp()
|
|
*/
|
|
void swi_polcart_sp(double *l, double *x)
|
|
{
|
|
double sinlon, coslon, sinlat, coslat;
|
|
double xx[6], rxy, rxyz;
|
|
/* zero speed */
|
|
if (l[3] == 0 && l[4] == 0 && l[5] == 0) {
|
|
x[3] = x[4] = x[5] = 0;
|
|
swi_polcart(l, x);
|
|
return;
|
|
}
|
|
/* position */
|
|
coslon = cos(l[0]);
|
|
sinlon = sin(l[0]);
|
|
coslat = cos(l[1]);
|
|
sinlat = sin(l[1]);
|
|
xx[0] = l[2] * coslat * coslon;
|
|
xx[1] = l[2] * coslat * sinlon;
|
|
xx[2] = l[2] * sinlat;
|
|
/* speed; explanation s. swi_cartpol_sp(), same method the other way round*/
|
|
rxyz = l[2];
|
|
rxy = sqrt(xx[0] * xx[0] + xx[1] * xx[1]);
|
|
xx[5] = l[5];
|
|
xx[4] = l[4] * rxyz;
|
|
x[5] = sinlat * xx[5] + coslat * xx[4]; /* speed z */
|
|
xx[3] = coslat * xx[5] - sinlat * xx[4];
|
|
xx[4] = l[3] * rxy;
|
|
x[3] = coslon * xx[3] - sinlon * xx[4]; /* speed x */
|
|
x[4] = sinlon * xx[3] + coslon * xx[4]; /* speed y */
|
|
x[0] = xx[0]; /* return position */
|
|
x[1] = xx[1];
|
|
x[2] = xx[2];
|
|
}
|
|
|
|
double swi_dot_prod_unit(double *x, double *y)
|
|
{
|
|
double dop = x[0]*y[0]+x[1]*y[1]+x[2]*y[2];
|
|
double e1 = sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
|
|
double e2 = sqrt(y[0]*y[0]+y[1]*y[1]+y[2]*y[2]);
|
|
dop /= e1;
|
|
dop /= e2;
|
|
if (dop > 1)
|
|
dop = 1;
|
|
if (dop < -1)
|
|
dop = -1;
|
|
return dop;
|
|
}
|
|
|
|
/* functions for precession and ecliptic obliquity according to Vondrák et alii, 2011 */
|
|
#define AS2R (DEGTORAD / 3600.0)
|
|
#define D2PI TWOPI
|
|
#define EPS0 (84381.406 * AS2R)
|
|
#define NPOL_PEPS 4
|
|
#define NPER_PEPS 10
|
|
#define NPOL_PECL 4
|
|
#define NPER_PECL 8
|
|
#define NPOL_PEQU 4
|
|
#define NPER_PEQU 14
|
|
|
|
/* for pre_peps(): */
|
|
/* polynomials */
|
|
static const double pepol[NPOL_PEPS][2] = {
|
|
{+8134.017132, +84028.206305},
|
|
{+5043.0520035, +0.3624445},
|
|
{-0.00710733, -0.00004039},
|
|
{+0.000000271, -0.000000110}
|
|
};
|
|
|
|
/* periodics */
|
|
static const double peper[5][NPER_PEPS] = {
|
|
{+409.90, +396.15, +537.22, +402.90, +417.15, +288.92, +4043.00, +306.00, +277.00, +203.00},
|
|
{-6908.287473, -3198.706291, +1453.674527, -857.748557, +1173.231614, -156.981465, +371.836550, -216.619040, +193.691479, +11.891524},
|
|
{+753.872780, -247.805823, +379.471484, -53.880558, -90.109153, -353.600190, -63.115353, -28.248187, +17.703387, +38.911307},
|
|
{-2845.175469, +449.844989, -1255.915323, +886.736783, +418.887514, +997.912441, -240.979710, +76.541307, -36.788069, -170.964086},
|
|
{-1704.720302, -862.308358, +447.832178, -889.571909, +190.402846, -56.564991, -296.222622, -75.859952, +67.473503, +3.014055}
|
|
};
|
|
|
|
/* for pre_pecl(): */
|
|
/* polynomials */
|
|
static const double pqpol[NPOL_PECL][2] = {
|
|
{+5851.607687, -1600.886300},
|
|
{-0.1189000, +1.1689818},
|
|
{-0.00028913, -0.00000020},
|
|
{+0.000000101, -0.000000437}
|
|
};
|
|
|
|
/* periodics */
|
|
static const double pqper[5][NPER_PECL] = {
|
|
{708.15, 2309, 1620, 492.2, 1183, 622, 882, 547},
|
|
{-5486.751211, -17.127623, -617.517403, 413.44294, 78.614193, -180.732815, -87.676083, 46.140315},
|
|
// original publication A&A 534, A22 (2011):
|
|
//{-684.66156, 2446.28388, 399.671049, -356.652376, -186.387003, -316.80007, 198.296071, 101.135679},
|
|
// typo fixed according to A&A 541, C1 (2012)
|
|
{-684.66156, 2446.28388, 399.671049, -356.652376, -186.387003, -316.80007, 198.296701, 101.135679},
|
|
{667.66673, -2354.886252, -428.152441, 376.202861, 184.778874, 335.321713, -185.138669, -120.97283},
|
|
{-5523.863691, -549.74745, -310.998056, 421.535876, -36.776172, -145.278396, -34.74445, 22.885731}
|
|
};
|
|
|
|
/* for pre_pequ(): */
|
|
/* polynomials */
|
|
static const double xypol[NPOL_PEQU][2] = {
|
|
{+5453.282155, -73750.930350},
|
|
{+0.4252841, -0.7675452},
|
|
{-0.00037173, -0.00018725},
|
|
{-0.000000152, +0.000000231}
|
|
};
|
|
|
|
/* periodics */
|
|
static const double xyper[5][NPER_PEQU] = {
|
|
{256.75, 708.15, 274.2, 241.45, 2309, 492.2, 396.1, 288.9, 231.1, 1610, 620, 157.87, 220.3, 1200},
|
|
{-819.940624, -8444.676815, 2600.009459, 2755.17563, -167.659835, 871.855056, 44.769698, -512.313065, -819.415595, -538.071099, -189.793622, -402.922932, 179.516345, -9.814756},
|
|
{75004.344875, 624.033993, 1251.136893, -1102.212834, -2660.66498, 699.291817, 153.16722, -950.865637, 499.754645, -145.18821, 558.116553, -23.923029, -165.405086, 9.344131},
|
|
{81491.287984, 787.163481, 1251.296102, -1257.950837, -2966.79973, 639.744522, 131.600209, -445.040117, 584.522874, -89.756563, 524.42963, -13.549067, -210.157124, -44.919798},
|
|
{1558.515853, 7774.939698, -2219.534038, -2523.969396, 247.850422, -846.485643, -1393.124055, 368.526116, 749.045012, 444.704518, 235.934465, 374.049623, -171.33018, -22.899655}
|
|
};
|
|
|
|
void swi_ldp_peps(double tjd, double *dpre, double *deps)
|
|
{
|
|
int i;
|
|
int npol = NPOL_PEPS;
|
|
int nper = NPER_PEPS;
|
|
double t, p, q, w, a, s, c;
|
|
t = (tjd - J2000) / 36525.0;
|
|
p = 0;
|
|
q = 0;
|
|
/* periodic terms */
|
|
for (i = 0; i < nper; i++) {
|
|
w = D2PI * t;
|
|
a = w / peper[0][i];
|
|
s = sin(a);
|
|
c = cos(a);
|
|
p += c * peper[1][i] + s * peper[3][i];
|
|
q += c * peper[2][i] + s * peper[4][i];
|
|
}
|
|
/* polynomial terms */
|
|
w = 1;
|
|
for (i = 0; i < npol; i++) {
|
|
p += pepol[i][0] * w;
|
|
q += pepol[i][1] * w;
|
|
w *= t;
|
|
}
|
|
/* both to radians */
|
|
p *= AS2R;
|
|
q *= AS2R;
|
|
/* return */
|
|
if (dpre != NULL)
|
|
*dpre = p;
|
|
if (deps != NULL)
|
|
*deps = q;
|
|
//fprintf(stderr, "%.17f\n", *deps / DEGTORAD);
|
|
}
|
|
|
|
/*
|
|
* Long term high precision precession,
|
|
* according to Vondrak/Capitaine/Wallace, "New precession expressions, valid
|
|
* for long time intervals", in A&A 534, A22(2011).
|
|
*/
|
|
/* precession of the ecliptic */
|
|
static void pre_pecl(double tjd, double *vec)
|
|
{
|
|
int i;
|
|
int npol = NPOL_PECL;
|
|
int nper = NPER_PECL;
|
|
double t, p, q, w, a, s, c, z;
|
|
t = (tjd - J2000) / 36525.0;
|
|
p = 0;
|
|
q = 0;
|
|
/* periodic terms */
|
|
for (i = 0; i < nper; i++) {
|
|
w = D2PI * t;
|
|
a = w / pqper[0][i];
|
|
s = sin(a);
|
|
c = cos(a);
|
|
p += c * pqper[1][i] + s * pqper[3][i];
|
|
q += c * pqper[2][i] + s * pqper[4][i];
|
|
}
|
|
/* polynomial terms */
|
|
w = 1;
|
|
for (i = 0; i < npol; i++) {
|
|
p += pqpol[i][0] * w;
|
|
q += pqpol[i][1] * w;
|
|
w *= t;
|
|
}
|
|
/* both to radians */
|
|
p *= AS2R;
|
|
q *= AS2R;
|
|
/* ecliptic pole vector */
|
|
z = 1 - p * p - q * q;
|
|
if (z < 0)
|
|
z = 0;
|
|
else
|
|
z = sqrt(z);
|
|
s = sin(EPS0);
|
|
c = cos(EPS0);
|
|
vec[0] = p;
|
|
vec[1] = - q * c - z * s;
|
|
vec[2] = - q * s + z * c;
|
|
}
|
|
|
|
/* precession of the equator */
|
|
static void pre_pequ(double tjd, double *veq)
|
|
{
|
|
int i;
|
|
int npol = NPOL_PEQU;
|
|
int nper = NPER_PEQU;
|
|
double t, x, y, w, a, s, c;
|
|
t = (tjd - J2000) / 36525.0;
|
|
x = 0;
|
|
y = 0;
|
|
for (i = 0; i < nper; i++) {
|
|
w = D2PI * t;
|
|
a = w / xyper[0][i];
|
|
s = sin(a);
|
|
c = cos(a);
|
|
x += c * xyper[1][i] + s * xyper[3][i];
|
|
y += c * xyper[2][i] + s * xyper[4][i];
|
|
}
|
|
/* polynomial terms */
|
|
w = 1;
|
|
for (i = 0; i < npol; i++) {
|
|
x += xypol[i][0] * w;
|
|
y += xypol[i][1] * w;
|
|
w *= t;
|
|
}
|
|
x *= AS2R;
|
|
y *= AS2R;
|
|
/* equator pole vector */
|
|
veq[0] = x;
|
|
veq[1] = y;
|
|
w = x * x + y * y;
|
|
if (w < 1)
|
|
veq[2] = sqrt(1 - w);
|
|
else
|
|
veq[2] = 0;
|
|
}
|
|
|
|
#if 0
|
|
static void swi_cross_prod(double *a, double *b, double *x)
|
|
{
|
|
x[0] = a[1] * b[2] - a[2] * b[1];
|
|
x[1] = a[2] * b[0] - a[0] * b[2];
|
|
x[2] = a[0] * b[1] - a[1] * b[0];
|
|
}
|
|
#endif
|
|
|
|
/* precession matrix */
|
|
static void pre_pmat(double tjd, double *rp)
|
|
{
|
|
double peqr[3], pecl[3], v[3], w, eqx[3];
|
|
//tjd = 1219339.078000;
|
|
/*equator pole */
|
|
pre_pequ(tjd, peqr);
|
|
/* ecliptic pole */
|
|
pre_pecl(tjd, pecl);
|
|
// fprintf(stderr, "%.17f %.17f %.17f\n", peqr[0], peqr[1], peqr[2]);
|
|
// fprintf(stderr, "%.17f %.17f %.17f\n", pecl[0], pecl[1], pecl[2]);
|
|
/* equinox */
|
|
swi_cross_prod(peqr, pecl, v);
|
|
w = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
|
|
eqx[0] = v[0] / w;
|
|
eqx[1] = v[1] / w;
|
|
eqx[2] = v[2] / w;
|
|
swi_cross_prod(peqr, eqx, v);
|
|
rp[0] = eqx[0];
|
|
rp[1] = eqx[1];
|
|
rp[2] = eqx[2];
|
|
rp[3] = v[0];
|
|
rp[4] = v[1];
|
|
rp[5] = v[2];
|
|
rp[6] = peqr[0];
|
|
rp[7] = peqr[1];
|
|
rp[8] = peqr[2];
|
|
// int i;
|
|
// for (i = 0; i < 3; i++) {
|
|
// fprintf(stderr, "(%.17f %.17f %.17f)\n", rp[i*3], rp[i*3+1],rp[i*3+2]);
|
|
// } /**/
|
|
}
|
|
|
|
/* precession according to Owen 1990:
|
|
* Owen, William M., Jr., (JPL) "A Theory of the Earth's Precession
|
|
* Relative to the Invariable Plane of the Solar System", Ph.D.
|
|
* Dissertation, University of Florida, 1990.
|
|
* Implemented for time range -18000 to 14000.
|
|
*/
|
|
/*
|
|
* p. 177: central time Tc = -160, covering time span -200 <= T <= -120
|
|
* i.e. -14000 +- 40 centuries
|
|
* p. 178: central time Tc = -80, covering time span -120 <= T <= -40
|
|
* i.e. -6000 +- 40 centuries
|
|
* p. 179: central time Tc = 0, covering time span -40 <= T <= +40
|
|
* i.e. 2000 +- 40 centuries
|
|
* p. 180: central time Tc = 80, covering time span 40 <= T <= 120
|
|
* i.e. 10000 +- 40 centuries
|
|
* p. 181: central time Tc = 160, covering time span 120 <= T <= 200
|
|
* i.e. 10000 +- 40 centuries
|
|
*/
|
|
static const double owen_eps0_coef[5][10] = {
|
|
{23.699391439256386, 5.2330816033981775e-1, -5.6259493384864815e-2, -8.2033318431602032e-3, 6.6774163554156385e-4, 2.4931584012812606e-5, -3.1313623302407878e-6, 2.0343814827951515e-7, 2.9182026615852936e-8, -4.1118760893281951e-9,},
|
|
{24.124759551704588, -1.2094875596566286e-1, -8.3914869653015218e-2, 3.5357075322387405e-3, 6.4557467824807032e-4, -2.5092064378707704e-5, -1.7631607274450848e-6, 1.3363622791424094e-7, 1.5577817511054047e-8, -2.4613907093017122e-9,},
|
|
{23.439103144206208, -4.9386077073143590e-1, -2.3965445283267805e-4, 8.6637485629656489e-3, -5.2828151901367600e-5, -4.3951004595359217e-5, -1.1058785949914705e-6, 6.2431490022621172e-8, 3.4725376218710764e-8, 1.3658853127005757e-9,},
|
|
{22.724671295125046, -1.6041813558650337e-1, 7.0646783888132504e-2, 1.4967806745062837e-3, -6.6857270989190734e-4, 5.7578378071604775e-6, 3.3738508454638728e-6, -2.2917813537654764e-7, -2.1019907929218137e-8, 4.3139832091694682e-9,},
|
|
{22.914636050333696, 3.2123508304962416e-1, 3.6633220173792710e-2, -5.9228324767696043e-3, -1.882379107379328e-4, 3.2274552870236244e-5, 4.9052463646336507e-7, -5.9064298731578425e-8, -2.0485712675098837e-8, -6.2163304813908160e-10,},
|
|
};
|
|
|
|
static const double owen_psia_coef[5][10] = {
|
|
{-218.57864954903122, 51.752257487741612, 1.3304715765661958e-1, 9.2048123521890745e-2, -6.0877528127241278e-3, -7.0013893644531700e-5, -4.9217728385458495e-5, -1.8578234189053723e-6, 7.4396426162029877e-7, -5.9157528981843864e-9,},
|
|
{-111.94350527506128, 55.175558131675861, 4.7366115762797613e-1, -4.7701750975398538e-2, -9.2445765329325809e-3, 7.0962838707454917e-4, 1.5140455277814658e-4, -7.7813159018954928e-7, -2.4729402281953378e-6, -1.0898887008726418e-7,},
|
|
{-2.041452011529441e-1, 55.969995858494106, -1.9295093699770936e-1, -5.6819574830421158e-3, 1.1073687302518981e-2, -9.0868489896815619e-5, -1.1999773777895820e-4, 9.9748697306154409e-6, 5.7911493603430550e-7, -2.3647526839778175e-7,},
|
|
{111.61366860604471, 56.404525305162447, 4.4403302410703782e-1, 7.1490030578883907e-2, -4.9184559079790816e-3, -1.3912698949042046e-3, -6.8490613661884005e-5, 1.2394328562905297e-6, 1.7719847841480384e-6, 2.4889095220628068e-7,},
|
|
{228.40683531269390, 60.056143904919826, 2.9583200718478960e-2, -1.5710838319490748e-1, -7.0017356811600801e-3, 3.3009615142224537e-3, 2.0318123852537664e-4, -6.5840216067828310e-5, -5.9077673352976155e-6, 1.3983942185303064e-6,},
|
|
};
|
|
|
|
static const double owen_oma_coef[5][10] = {
|
|
{25.541291140949806, 2.377889511272162e-1, -3.7337334723142133e-1, 2.4579295485161534e-2, 4.3840999514263623e-3, -3.1126873333599556e-4, -9.8443045771748915e-6, -7.9403103080496923e-7, 1.0840116743893556e-9, 9.2865105216887919e-9,},
|
|
{24.429357654237926, -9.5205745947740161e-1, 8.6738296270534816e-2, 3.0061543426062955e-2, -4.1532480523019988e-3, -3.7920928393860939e-4, 3.5117012399609737e-5, 4.6811877283079217e-6, -8.1836046585546861e-8, -6.1803706664211173e-8,},
|
|
{23.450465062489337, -9.7259278279739817e-2, 1.1082286925130981e-2, -3.1469883339372219e-2, -1.0041906996819648e-4, 5.6455168475133958e-4, -8.4403910211030209e-6, -3.8269157371098435e-6, 3.1422585261198437e-7, 9.3481729116773404e-9,},
|
|
{22.581778052947806, -8.7069701538602037e-1, -9.8140710050197307e-2, 2.6025931340678079e-2, 4.8165322168786755e-3, -1.906558772193363e-4, -4.6838759635421777e-5, -1.6608525315998471e-6, -3.2347811293516124e-8, 2.8104728109642000e-9,},
|
|
{21.518861835737142, 2.0494789509441385e-1, 3.5193604846503161e-1, 1.5305977982348925e-2, -7.5015367726336455e-3, -4.0322553186065610e-4, 1.0655320434844041e-4, 7.1792339586935752e-6, -1.603874697543020e-6, -1.613563462813512e-7,},
|
|
};
|
|
|
|
static const double owen_chia_coef[5][10] = {
|
|
{8.2378850337329404e-1, -3.7443109739678667, 4.0143936898854026e-1, 8.1822830214590811e-2, -8.5978790792656293e-3, -2.8350488448426132e-5, -4.2474671728156727e-5, -1.6214840884656678e-6, 7.8560442001953050e-7, -1.032016641696707e-8,},
|
|
{-2.1726062070318606, 7.8470515033132925e-1, 4.4044931004195718e-1, -8.0671247169971653e-2, -8.9672662444325007e-3, 9.2248978383109719e-4, 1.5143472266372874e-4, -1.6387009056475679e-6, -2.4405558979328144e-6, -1.0148113464009015e-7,},
|
|
{-4.8518673570735556e-1, 1.0016737299946743e-1, -4.7074888613099918e-1, -5.8604054305076092e-3, 1.4300208240553435e-2, -6.7127991650300028e-5, -1.3703764889645475e-4, 9.0505213684444634e-6, 6.0368690647808607e-7, -2.2135404747652171e-7,},
|
|
{-2.0950740076326087, -9.4447359463206877e-1, 4.0940512860493755e-1, 1.0261699700263508e-1, -5.3133241571955160e-3, -1.6634631550720911e-3, -5.9477519536647907e-5, 2.9651387319208926e-6, 1.6434499452070584e-6, 2.3720647656961084e-7,},
|
|
{6.3315163285678715e-1, 3.5241082918420464, 2.1223076605364606e-1, -1.5648122502767368e-1, -9.1964075390801980e-3, 3.3896161239812411e-3, 2.1485178626085787e-4, -6.6261759864793735e-5, -5.9257969712852667e-6, 1.3918759086160525e-6,},
|
|
};
|
|
|
|
static void get_owen_t0_icof(double tjd, double *t0, int *icof)
|
|
{
|
|
int i, j = 0;
|
|
double t0s[5] = {-3392455.5, -470455.5, 2451544.5, 5373544.5, 8295544.5, };
|
|
*t0 = t0s[0];
|
|
for (i = 1; i < 5; i++) {
|
|
if (tjd < (t0s[i-1] + t0s[i]) / 2) {
|
|
;
|
|
} else {
|
|
*t0 = t0s[i];
|
|
j++;
|
|
}
|
|
}
|
|
*icof = j;
|
|
}
|
|
|
|
/* precession matrix Owen 1990 */
|
|
static void owen_pre_matrix(double tjd, double *rp, int iflag)
|
|
{
|
|
int i, icof = 0;
|
|
double eps0 = 0, chia = 0, psia = 0, oma = 0;
|
|
double coseps0, sineps0, coschia, sinchia, cospsia, sinpsia, cosoma, sinoma;
|
|
double k[10], tau[10];
|
|
double t0;
|
|
get_owen_t0_icof(tjd, &t0, &icof);
|
|
// fprintf(stderr, "%d, %.17f\n", icof, t0);
|
|
tau[0] = 0;
|
|
tau[1] = (tjd - t0) / 36525.0 / 40.0;
|
|
for (i = 2; i <= 9; i++) {
|
|
tau[i] = tau[1] * tau[i-1];
|
|
}
|
|
k[0] = 1;
|
|
k[1] = tau[1];
|
|
k[2] = 2 * tau[2] - 1;
|
|
k[3] = 4 * tau[3] - 3 * tau[1];
|
|
k[4] = 8 * tau[4] - 8 * tau[2] + 1;
|
|
k[5] = 16 * tau[5] - 20 * tau[3] + 5 * tau[1];
|
|
k[6] = 32 * tau[6] - 48 * tau[4] + 18 * tau[2] - 1;
|
|
k[7] = 64 * tau[7] - 112 * tau[5] + 56 * tau[3] - 7 * tau[1];
|
|
k[8] = 128 * tau[8] - 256 * tau[6] + 160 * tau[4] - 32 * tau[2] + 1;
|
|
k[9] = 256 * tau[9] - 576 * tau[7] + 432 * tau[5] - 120 * tau[3] + 9 * tau[1];
|
|
for (i = 0; i < 10; i++) {
|
|
//eps += (k[i] * owen_eps0_coef[icof][i]);
|
|
psia += (k[i] * owen_psia_coef[icof][i]);
|
|
oma += (k[i] * owen_oma_coef[icof][i]);
|
|
chia += (k[i] * owen_chia_coef[icof][i]);
|
|
}
|
|
if (iflag & (SEFLG_JPLHOR | SEFLG_JPLHOR_APPROX)) {
|
|
/*
|
|
* In comparison with JPL Horizons we have an almost constant offset
|
|
* almost constant offset in ecl. longitude of about -0.000019 deg.
|
|
* We fix this as follows: */
|
|
psia += -0.000018560;
|
|
}
|
|
eps0 = 84381.448 / 3600.0;
|
|
//fprintf(stderr, "e=%.17f, ps=%.17f, om=%.17f, ch=%.17f\n", eps0, psia, oma, chia);
|
|
eps0 *= DEGTORAD;
|
|
psia *= DEGTORAD;
|
|
chia *= DEGTORAD;
|
|
oma *= DEGTORAD;
|
|
coseps0 = cos(eps0);
|
|
sineps0 = sin(eps0);
|
|
coschia = cos(chia);
|
|
sinchia = sin(chia);
|
|
cospsia = cos(psia);
|
|
sinpsia = sin(psia);
|
|
cosoma = cos(oma);
|
|
sinoma = sin(oma);
|
|
rp[0] = coschia * cospsia + sinchia * cosoma * sinpsia;
|
|
rp[1] = (-coschia * sinpsia + sinchia * cosoma * cospsia) * coseps0 + sinchia * sinoma * sineps0;
|
|
rp[2] = (-coschia * sinpsia + sinchia * cosoma * cospsia) * sineps0 - sinchia * sinoma * coseps0;
|
|
rp[3] = -sinchia * cospsia + coschia * cosoma * sinpsia;
|
|
rp[4] = (sinchia * sinpsia + coschia * cosoma * cospsia) * coseps0 + coschia * sinoma * sineps0;
|
|
rp[5] = (sinchia * sinpsia + coschia * cosoma * cospsia) * sineps0 - coschia * sinoma * coseps0;
|
|
rp[6] = sinoma * sinpsia;
|
|
rp[7] = sinoma * cospsia * coseps0 - cosoma * sineps0;
|
|
rp[8] = sinoma * cospsia * sineps0 + cosoma * coseps0;
|
|
/*for (i = 0; i < 3; i++) {
|
|
fprintf(stderr, "(%.17f %.17f %.17f)\n", rp[i*3], rp[i*3+1],rp[i*3+2]);
|
|
} */
|
|
}
|
|
static void epsiln_owen_1986(double tjd, double *eps)
|
|
{
|
|
int i, icof = 0;
|
|
double k[10], tau[10];
|
|
double t0;
|
|
get_owen_t0_icof(tjd, &t0, &icof);
|
|
*eps = 0;
|
|
tau[0] = 0;
|
|
tau[1] = (tjd - t0) / 36525.0 / 40.0;
|
|
for (i = 2; i <= 9; i++) {
|
|
tau[i] = tau[1] * tau[i-1];
|
|
}
|
|
k[0] = 1;
|
|
k[1] = tau[1];
|
|
k[2] = 2 * tau[2] - 1;
|
|
k[3] = 4 * tau[3] - 3 * tau[1];
|
|
k[4] = 8 * tau[4] - 8 * tau[2] + 1;
|
|
k[5] = 16 * tau[5] - 20 * tau[3] + 5 * tau[1];
|
|
k[6] = 32 * tau[6] - 48 * tau[4] + 18 * tau[2] - 1;
|
|
k[7] = 64 * tau[7] - 112 * tau[5] + 56 * tau[3] - 7 * tau[1];
|
|
k[8] = 128 * tau[8] - 256 * tau[6] + 160 * tau[4] - 32 * tau[2] + 1;
|
|
k[9] = 256 * tau[9] - 576 * tau[7] + 432 * tau[5] - 120 * tau[3] + 9 * tau[1];
|
|
for (i = 0; i < 10; i++) {
|
|
*eps += (k[i] * owen_eps0_coef[icof][i]);
|
|
}
|
|
//fprintf(stderr, "eps=%.17f\n", *eps);
|
|
}
|
|
|
|
/* Obliquity of the ecliptic at Julian date J
|
|
*
|
|
* IAU Coefficients are from:
|
|
* J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
|
|
* "Expressions for the Precession Quantities Based upon the IAU
|
|
* (1976) System of Astronomical Constants," Astronomy and Astrophysics
|
|
* 58, 1-16 (1977).
|
|
*
|
|
* Before or after 200 years from J2000, the formula used is from:
|
|
* J. Laskar, "Secular terms of classical planetary theories
|
|
* using the results of general theory," Astronomy and Astrophysics
|
|
* 157, 59070 (1986).
|
|
*
|
|
* Bretagnon, P. et al.: 2003, "Expressions for Precession Consistent with
|
|
* the IAU 2000A Model". A&A 400,785
|
|
*B03 84381.4088 -46.836051*t -1667*10-7*t2 +199911*10-8*t3 -523*10-9*t4 -248*10-10*t5 -3*10-11*t6
|
|
*C03 84381.406 -46.836769*t -1831*10-7*t2 +20034*10-7*t3 -576*10-9*t4 -434*10-10*t5
|
|
*
|
|
* See precess and page B18 of the Astronomical Almanac.
|
|
*/
|
|
#define OFFSET_EPS_JPLHORIZONS (35.95)
|
|
#define DCOR_EPS_JPL_TJD0 2437846.5
|
|
#define NDCOR_EPS_JPL 51
|
|
double dcor_eps_jpl[] = {
|
|
36.726, 36.627, 36.595, 36.578, 36.640, 36.659, 36.731, 36.765,
|
|
36.662, 36.555, 36.335, 36.321, 36.354, 36.227, 36.289, 36.348, 36.257, 36.163,
|
|
35.979, 35.896, 35.842, 35.825, 35.912, 35.950, 36.093, 36.191, 36.009, 35.943,
|
|
35.875, 35.771, 35.788, 35.753, 35.822, 35.866, 35.771, 35.732, 35.543, 35.498,
|
|
35.449, 35.409, 35.497, 35.556, 35.672, 35.760, 35.596, 35.565, 35.510, 35.394,
|
|
35.385, 35.375, 35.415,
|
|
};
|
|
double swi_epsiln(double J, int32 iflag)
|
|
{
|
|
double T, eps;
|
|
double tofs, dofs, t0, t1;
|
|
int prec_model = swed.astro_models[SE_MODEL_PREC_LONGTERM];
|
|
int prec_model_short = swed.astro_models[SE_MODEL_PREC_SHORTTERM];
|
|
int jplhora_model = swed.astro_models[SE_MODEL_JPLHORA_MODE];
|
|
AS_BOOL is_jplhor = FALSE;
|
|
if (prec_model == 0) prec_model = SEMOD_PREC_DEFAULT;
|
|
if (prec_model_short == 0) prec_model_short = SEMOD_PREC_DEFAULT_SHORT;
|
|
if (jplhora_model == 0) jplhora_model = SEMOD_JPLHORA_DEFAULT;
|
|
if (iflag & SEFLG_JPLHOR)
|
|
is_jplhor = TRUE;
|
|
if ((iflag & SEFLG_JPLHOR_APPROX)
|
|
&& jplhora_model == SEMOD_JPLHORA_3
|
|
&& J <= HORIZONS_TJD0_DPSI_DEPS_IAU1980)
|
|
is_jplhor = TRUE;
|
|
T = (J - 2451545.0)/36525.0;
|
|
if (is_jplhor) {
|
|
if (J > 2378131.5 && J < 2525323.5) { // between 1.1.1799 and 1.1.2202
|
|
eps = (((1.813e-3*T-5.9e-4)*T-46.8150)*T+84381.448)*DEGTORAD/3600;
|
|
} else {
|
|
epsiln_owen_1986(J, &eps);
|
|
eps *= DEGTORAD;
|
|
}
|
|
} else if ((iflag & SEFLG_JPLHOR_APPROX) && jplhora_model == SEMOD_JPLHORA_2) {
|
|
eps = (((1.813e-3*T-5.9e-4)*T-46.8150)*T+84381.448)*DEGTORAD/3600;
|
|
} else if (prec_model_short == SEMOD_PREC_IAU_1976 && fabs(T) <= PREC_IAU_1976_CTIES ) {
|
|
eps = (((1.813e-3*T-5.9e-4)*T-46.8150)*T+84381.448)*DEGTORAD/3600;
|
|
} else if (prec_model == SEMOD_PREC_IAU_1976) {
|
|
eps = (((1.813e-3*T-5.9e-4)*T-46.8150)*T+84381.448)*DEGTORAD/3600;
|
|
} else if (prec_model_short == SEMOD_PREC_IAU_2000 && fabs(T) <= PREC_IAU_2000_CTIES ) {
|
|
eps = (((1.813e-3*T-5.9e-4)*T-46.84024)*T+84381.406)*DEGTORAD/3600;
|
|
} else if (prec_model == SEMOD_PREC_IAU_2000) {
|
|
eps = (((1.813e-3*T-5.9e-4)*T-46.84024)*T+84381.406)*DEGTORAD/3600;
|
|
} else if (prec_model_short == SEMOD_PREC_IAU_2006 && fabs(T) <= PREC_IAU_2006_CTIES) {
|
|
eps = (((((-4.34e-8 * T -5.76e-7) * T +2.0034e-3) * T -1.831e-4) * T -46.836769) * T + 84381.406) * DEGTORAD / 3600.0;
|
|
} else if (prec_model == SEMOD_PREC_NEWCOMB) {
|
|
double Tn = (J - 2396758.0)/36525.0;
|
|
eps = (0.0017 * Tn * Tn * Tn - 0.0085 * Tn * Tn - 46.837 * Tn + 84451.68) * DEGTORAD / 3600.0;
|
|
} else if (prec_model == SEMOD_PREC_IAU_2006) {
|
|
eps = (((((-4.34e-8 * T -5.76e-7) * T +2.0034e-3) * T -1.831e-4) * T -46.836769) * T + 84381.406) * DEGTORAD / 3600.0;
|
|
} else if (prec_model == SEMOD_PREC_BRETAGNON_2003) {
|
|
eps = ((((((-3e-11 * T - 2.48e-8) * T -5.23e-7) * T +1.99911e-3) * T -1.667e-4) * T -46.836051) * T + 84381.40880) * DEGTORAD / 3600.0;/* */
|
|
} else if (prec_model == SEMOD_PREC_SIMON_1994) {
|
|
eps = (((((2.5e-8 * T -5.1e-7) * T +1.9989e-3) * T -1.52e-4) * T -46.80927) * T + 84381.412) * DEGTORAD / 3600.0;/* */
|
|
} else if (prec_model == SEMOD_PREC_WILLIAMS_1994) {
|
|
eps = ((((-1.0e-6 * T +2.0e-3) * T -1.74e-4) * T -46.833960) * T + 84381.409) * DEGTORAD / 3600.0;/* */
|
|
} else if (prec_model == SEMOD_PREC_LASKAR_1986 || prec_model == SEMOD_PREC_WILL_EPS_LASK) {
|
|
T /= 10.0;
|
|
eps = ((((((((( 2.45e-10*T + 5.79e-9)*T + 2.787e-7)*T
|
|
+ 7.12e-7)*T - 3.905e-5)*T - 2.4967e-3)*T
|
|
- 5.138e-3)*T + 1.99925)*T - 0.0155)*T - 468.093)*T
|
|
+ 84381.448;
|
|
eps *= DEGTORAD/3600.0;
|
|
} else if (prec_model == SEMOD_PREC_OWEN_1990) {
|
|
epsiln_owen_1986(J, &eps);
|
|
eps *= DEGTORAD;
|
|
//fprintf(stderr, "epso=%.17f\n", eps);
|
|
} else { /* SEMOD_PREC_VONDRAK_2011 */
|
|
swi_ldp_peps(J, NULL, &eps);
|
|
if ((iflag & SEFLG_JPLHOR_APPROX) && jplhora_model != SEMOD_JPLHORA_2) {
|
|
tofs = (J - DCOR_EPS_JPL_TJD0) / 365.25;
|
|
dofs = OFFSET_EPS_JPLHORIZONS;
|
|
if (tofs < 0) {
|
|
tofs = 0;
|
|
dofs = dcor_eps_jpl[0];
|
|
} else if (tofs >= NDCOR_EPS_JPL - 1) {
|
|
tofs = NDCOR_EPS_JPL;
|
|
dofs = dcor_eps_jpl[NDCOR_EPS_JPL - 1];
|
|
} else {
|
|
t0 = (int) tofs;
|
|
t1 = t0 + 1;
|
|
dofs = dcor_eps_jpl[(int)t0];
|
|
dofs = (tofs - t0) * (dcor_eps_jpl[(int)t0] - dcor_eps_jpl[(int)t1]) + dcor_eps_jpl[(int)t0];
|
|
}
|
|
dofs /= (1000.0 * 3600.0);
|
|
eps += dofs * DEGTORAD;
|
|
}
|
|
//fprintf(stderr, "epsv=%.17f\n", eps);
|
|
}
|
|
return(eps);
|
|
}
|
|
|
|
/* Precession of the equinox and ecliptic
|
|
* from epoch Julian date J to or from J2000.0
|
|
*
|
|
* Original program by Steve Moshier.
|
|
* Changes in program structure and implementation of IAU 2003 (P03) and
|
|
* Vondrak 2011 by Dieter Koch.
|
|
*
|
|
* SEMOD_PREC_VONDRAK_2011
|
|
* J. Vondrák, N. Capitaine, and P. Wallace, "New precession expressions,
|
|
* valid for long time intervals", A&A 534, A22 (2011)
|
|
*
|
|
* SEMOD_PREC_IAU_2006
|
|
* N. Capitaine, P.T. Wallace, and J. Chapront, "Expressions for IAU 2000
|
|
* precession quantities", 2003, A&A 412, 567-586 (2003).
|
|
* This is a "short" term model, that can be combined with other models
|
|
*
|
|
* SEMOD_PREC_WILLIAMS_1994
|
|
* James G. Williams, "Contributions to the Earth's obliquity rate,
|
|
* precession, and nutation," Astron. J. 108, 711-724 (1994).
|
|
*
|
|
* SEMOD_PREC_SIMON_1994
|
|
* J. L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touze', G. Francou,
|
|
* and J. Laskar, "Numerical Expressions for precession formulae and
|
|
* mean elements for the Moon and the planets," Astronomy and Astrophysics
|
|
* 282, 663-683 (1994).
|
|
*
|
|
* SEMOD_PREC_IAU_1976
|
|
* IAU Coefficients are from:
|
|
* J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
|
|
* "Expressions for the Precession Quantities Based upon the IAU
|
|
* (1976) System of Astronomical Constants," Astronomy and
|
|
* Astrophysics 58, 1-16 (1977).
|
|
* This is a "short" term model, that can be combined with other models
|
|
*
|
|
* SEMOD_PREC_LASKAR_1986
|
|
* Newer formulas that cover a much longer time span are from:
|
|
* J. Laskar, "Secular terms of classical planetary theories
|
|
* using the results of general theory," Astronomy and Astrophysics
|
|
* 157, 59070 (1986).
|
|
*
|
|
* See also:
|
|
* P. Bretagnon and G. Francou, "Planetary theories in rectangular
|
|
* and spherical variables. VSOP87 solutions," Astronomy and
|
|
* Astrophysics 202, 309-315 (1988).
|
|
*
|
|
* Bretagnon and Francou's expansions for the node and inclination
|
|
* of the ecliptic were derived from Laskar's data but were truncated
|
|
* after the term in T**6. I have recomputed these expansions from
|
|
* Laskar's data, retaining powers up to T**10 in the result.
|
|
*
|
|
*/
|
|
|
|
static int precess_1(double *R, double J, int direction, int prec_method)
|
|
{
|
|
double T, Z = 0, z = 0, TH = 0;
|
|
int i;
|
|
double x[3];
|
|
double sinth, costh, sinZ, cosZ, sinz, cosz, A, B;
|
|
if( J == J2000 )
|
|
return(0);
|
|
T = (J - J2000)/36525.0;
|
|
if (prec_method == SEMOD_PREC_IAU_1976) {
|
|
Z = (( 0.017998*T + 0.30188)*T + 2306.2181)*T*DEGTORAD/3600;
|
|
z = (( 0.018203*T + 1.09468)*T + 2306.2181)*T*DEGTORAD/3600;
|
|
TH = ((-0.041833*T - 0.42665)*T + 2004.3109)*T*DEGTORAD/3600;
|
|
/*
|
|
* precession relative to ecliptic of start epoch is:
|
|
* pn = (5029.0966 + 2.22226*T-0.000042*T*T) * t + (1.11161-0.000127*T) * t * t - 0.000113*t*t*t;
|
|
* with: t = (tstart - tdate) / 36525.0
|
|
* T = (tstart - J2000) / 36525.0
|
|
*/
|
|
} else if (prec_method == SEMOD_PREC_IAU_2000) {
|
|
/* AA 2006 B28:*/
|
|
Z = (((((- 0.0000002*T - 0.0000327)*T + 0.0179663)*T + 0.3019015)*T + 2306.0809506)*T + 2.5976176)*DEGTORAD/3600;
|
|
z = (((((- 0.0000003*T - 0.000047)*T + 0.0182237)*T + 1.0947790)*T + 2306.0803226)*T - 2.5976176)*DEGTORAD/3600;
|
|
TH = ((((-0.0000001*T - 0.0000601)*T - 0.0418251)*T - 0.4269353)*T + 2004.1917476)*T*DEGTORAD/3600;
|
|
} else if (prec_method == SEMOD_PREC_IAU_2006) {
|
|
T = (J - J2000)/36525.0;
|
|
Z = (((((- 0.0000003173*T - 0.000005971)*T + 0.01801828)*T + 0.2988499)*T + 2306.083227)*T + 2.650545)*DEGTORAD/3600;
|
|
z = (((((- 0.0000002904*T - 0.000028596)*T + 0.01826837)*T + 1.0927348)*T + 2306.077181)*T - 2.650545)*DEGTORAD/3600;
|
|
TH = ((((-0.00000011274*T - 0.000007089)*T - 0.04182264)*T - 0.4294934)*T + 2004.191903)*T*DEGTORAD/3600;
|
|
} else if (prec_method == SEMOD_PREC_BRETAGNON_2003) {
|
|
Z = ((((((-0.00000000013*T - 0.0000003040)*T - 0.000005708)*T + 0.01801752)*T + 0.3023262)*T + 2306.080472)*T + 2.72767)*DEGTORAD/3600;
|
|
z = ((((((-0.00000000005*T - 0.0000002486)*T - 0.000028276)*T + 0.01826676)*T + 1.0956768)*T + 2306.076070)*T - 2.72767)*DEGTORAD/3600;
|
|
TH = ((((((0.000000000009*T + 0.00000000036)*T -0.0000001127)*T - 0.000007291)*T - 0.04182364)*T - 0.4266980)*T + 2004.190936)*T*DEGTORAD/3600;
|
|
#if 0
|
|
} else if (prec_method == SEMOD_PREC_NEWCOMB) {
|
|
double t1 = (J2000 - 2415020.3135) / 36524.2199;
|
|
double T = (J - J2000) / 36524.2199;
|
|
double T2 = T * T; double T3 = T2 * T;
|
|
Z = (2304.250 + 1.396 * t1) * T + 0.302 * T2 + 0.0179 * T3;
|
|
z = (2304.250 + 1.396 * t1) * T + 1.093 * T2 + 0.0192 * T3;
|
|
TH =(2004.682 - 0.853 * t1) * T - 0.426 * T2 - 0.0416 * T3;
|
|
Z *= (DEGTORAD/3600.0);
|
|
z *= (DEGTORAD/3600.0);
|
|
TH *= (DEGTORAD/3600.0);
|
|
#endif
|
|
#if 0
|
|
// from Newcomb, "Compendium" (1906), pp. 245f., relative to 1850
|
|
/* } else if (prec_method == SEMOD_PREC_NEWCOMB) {
|
|
double cties = 36524.2198782; // trop. centuries
|
|
double T = (J - B1850) / cties;
|
|
double T2 = T * T; double T3 = T2 * T;
|
|
double Z1 = 2303.56;
|
|
Z = 2303.56 * T + 0.3023 * T2 + 0.018 * T3;
|
|
z = 2303.55 * T + 1.094 * T2 + 0.018 * T3;
|
|
TH = 2005.11 * T - 0.43 * T2 - 0.041 * T3;
|
|
Z *= (DEGTORAD/3600.0);
|
|
z *= (DEGTORAD/3600.0);
|
|
TH *= (DEGTORAD/3600.0);
|
|
*/
|
|
#endif
|
|
#if 0
|
|
// Newcomb from Expl. supp. 61 pg. 38
|
|
// "Andoyar (Woolard and Clemence) expressions":
|
|
} else if (prec_method == SEMOD_PREC_NEWCOMB) {
|
|
double mills = 365242.198782; // trop. millennia
|
|
double t1 = (J2000 - B1850) / mills;
|
|
double t2 = (J - B1850) / mills;
|
|
double T = t2 - t1;
|
|
double T2 = T * T; double T3 = T2 * T;
|
|
double Z1 = 23035.545 + 139.720 * t1 + 0.060 * t1 * t1;
|
|
Z = Z1 * T + (30.240 - 0.270 * t1) * T2 + 17.995 * T3;
|
|
z = Z1 * T + (109.480 - 0.390 * t1) * T2 + 18.325 * T3;
|
|
TH = (20051.12 - 85.29 * t1 - 0.37 * t1 * t1) * T + (-42.65 - 0.37 * t1) * T2 - 41.80 * T3;
|
|
Z *= (DEGTORAD/3600.0);
|
|
z *= (DEGTORAD/3600.0);
|
|
TH *= (DEGTORAD/3600.0);
|
|
#endif
|
|
#if 1
|
|
// Newcomb according to Kinoshita 1975, very close to ExplSuppl/Andoyer;
|
|
// one additional digit.
|
|
} else if (prec_method == SEMOD_PREC_NEWCOMB) {
|
|
double mills = 365242.198782; // trop. millennia
|
|
double t1 = (J2000 - B1850) / mills;
|
|
double t2 = (J - B1850) / mills;
|
|
double T = t2 - t1;
|
|
double T2 = T * T; double T3 = T2 * T;
|
|
double Z1 = 23035.5548 + 139.720 * t1 + 0.069 * t1 * t1;
|
|
Z = Z1 * T + (30.242 - 0.269 * t1) * T2 + 17.996 * T3;
|
|
z = Z1 * T + (109.478 - 0.387 * t1) * T2 + 18.324 * T3;
|
|
TH = (20051.125 - 85.294 * t1 - 0.365 * t1 * t1) * T + (-42.647 - 0.365 * t1) * T2 - 41.802 * T3;
|
|
Z *= (DEGTORAD/3600.0);
|
|
z *= (DEGTORAD/3600.0);
|
|
TH *= (DEGTORAD/3600.0);
|
|
#endif
|
|
#if 0
|
|
// from Lieske, "Expressions for the Precession Quantities..." (1967), p. 20
|
|
} else if (prec_method == SEMOD_PREC_NEWCOMB) {
|
|
double cties = 36524.2198782; // trop. centuries
|
|
double t1 = (J2000 - J1900) / cties;
|
|
double t2 = (J - J1900) / cties;
|
|
double T = t2 - t1;
|
|
double T2 = T * T; double T3 = T2 * T;
|
|
double Z1 = 2304.253 + 1.3972 * t1 + 0.000125 * t1 * t1;
|
|
Z = Z1 * T + (0.3023 - 0.000211 * t1) * T2 + 0.0180 * T3;
|
|
z = Z1 * T + (1.0949 - 0.00046 * t1) * T2 + 0.0183 * T3;
|
|
TH = (2004.684 - 0.8532 * t1 - 0.000317 * t1 * t1) * T + (-0.4266 - 0.00032 * t1) * T2 - 0.0418 * T3;
|
|
Z *= (DEGTORAD/3600.0);
|
|
z *= (DEGTORAD/3600.0);
|
|
TH *= (DEGTORAD/3600.0);
|
|
#endif
|
|
} else {
|
|
return 0;
|
|
}
|
|
sinth = sin(TH);
|
|
costh = cos(TH);
|
|
sinZ = sin(Z);
|
|
cosZ = cos(Z);
|
|
sinz = sin(z);
|
|
cosz = cos(z);
|
|
A = cosZ*costh;
|
|
B = sinZ*costh;
|
|
if( direction < 0 ) { /* From J2000.0 to J */
|
|
x[0] = (A*cosz - sinZ*sinz)*R[0]
|
|
- (B*cosz + cosZ*sinz)*R[1]
|
|
- sinth*cosz*R[2];
|
|
x[1] = (A*sinz + sinZ*cosz)*R[0]
|
|
- (B*sinz - cosZ*cosz)*R[1]
|
|
- sinth*sinz*R[2];
|
|
x[2] = cosZ*sinth*R[0]
|
|
- sinZ*sinth*R[1]
|
|
+ costh*R[2];
|
|
} else { /* From J to J2000.0 */
|
|
x[0] = (A*cosz - sinZ*sinz)*R[0]
|
|
+ (A*sinz + sinZ*cosz)*R[1]
|
|
+ cosZ*sinth*R[2];
|
|
x[1] = - (B*cosz + cosZ*sinz)*R[0]
|
|
- (B*sinz - cosZ*cosz)*R[1]
|
|
- sinZ*sinth*R[2];
|
|
x[2] = - sinth*cosz*R[0]
|
|
- sinth*sinz*R[1]
|
|
+ costh*R[2];
|
|
}
|
|
for( i=0; i<3; i++ )
|
|
R[i] = x[i];
|
|
return(0);
|
|
}
|
|
|
|
/* In WILLIAMS and SIMON, Laskar's terms of order higher than t^4
|
|
have been retained, because Simon et al mention that the solution
|
|
is the same except for the lower order terms. */
|
|
|
|
/* SEMOD_PREC_WILLIAMS_1994 */
|
|
static const double pAcof_williams[] = {
|
|
-8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
|
|
-0.235316, 0.076, 110.5407, 50287.70000 };
|
|
static const double nodecof_williams[] = {
|
|
6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 1.9e-10,
|
|
-3.54e-9, -1.8103e-7, 1.26e-7, 7.436169e-5,
|
|
-0.04207794833, 3.052115282424};
|
|
static const double inclcof_williams[] = {
|
|
1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11,
|
|
-5.4000441e-11, 1.32115526e-9, -6.012e-7, -1.62442e-5,
|
|
0.00227850649, 0.0 };
|
|
|
|
/* SEMOD_PREC_SIMON_1994 */
|
|
/* Precession coefficients from Simon et al: */
|
|
static const double pAcof_simon[] = {
|
|
-8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
|
|
-0.235316, 0.07732, 111.2022, 50288.200 };
|
|
static const double nodecof_simon[] = {
|
|
6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 1.9e-10,
|
|
-3.54e-9, -1.8103e-7, 2.579e-8, 7.4379679e-5,
|
|
-0.0420782900, 3.0521126906};
|
|
static const double inclcof_simon[] = {
|
|
1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11,
|
|
-5.4000441e-11, 1.32115526e-9, -5.99908e-7, -1.624383e-5,
|
|
0.002278492868, 0.0 };
|
|
|
|
/* SEMOD_PREC_LASKAR_1986 */
|
|
/* Precession coefficients taken from Laskar's paper: */
|
|
static const double pAcof_laskar[] = {
|
|
-8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
|
|
-0.235316, 0.07732, 111.1971, 50290.966 };
|
|
/* Node and inclination of the earth's orbit computed from
|
|
* Laskar's data as done in Bretagnon and Francou's paper.
|
|
* Units are radians.
|
|
*/
|
|
static const double nodecof_laskar[] = {
|
|
6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 6.3190131e-10,
|
|
-3.48388152e-9, -1.813065896e-7, 2.75036225e-8, 7.4394531426e-5,
|
|
-0.042078604317, 3.052112654975 };
|
|
static const double inclcof_laskar[] = {
|
|
1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11,
|
|
-5.4000441e-11, 1.32115526e-9, -5.998737027e-7, -1.6242797091e-5,
|
|
0.002278495537, 0.0 };
|
|
|
|
static int precess_2(double *R, double J, int32 iflag, int direction, int prec_method)
|
|
{
|
|
int i;
|
|
double T, z;
|
|
double eps, sineps, coseps;
|
|
double x[3];
|
|
const double *p;
|
|
double A, B, pA, W;
|
|
const double *pAcof, *inclcof, *nodecof;
|
|
if( J == J2000 )
|
|
return(0);
|
|
if (prec_method == SEMOD_PREC_LASKAR_1986) {
|
|
pAcof = pAcof_laskar;
|
|
nodecof = nodecof_laskar;
|
|
inclcof = inclcof_laskar;
|
|
} else if (prec_method == SEMOD_PREC_SIMON_1994) {
|
|
pAcof = pAcof_simon;
|
|
nodecof = nodecof_simon;
|
|
inclcof = inclcof_simon;
|
|
} else if (prec_method == SEMOD_PREC_WILLIAMS_1994) {
|
|
pAcof = pAcof_williams;
|
|
nodecof = nodecof_williams;
|
|
inclcof = inclcof_williams;
|
|
} else { /* default, to satisfy compiler */
|
|
pAcof = pAcof_laskar;
|
|
nodecof = nodecof_laskar;
|
|
inclcof = inclcof_laskar;
|
|
}
|
|
T = (J - J2000)/36525.0;
|
|
/* Implementation by elementary rotations using Laskar's expansions.
|
|
* First rotate about the x axis from the initial equator
|
|
* to the ecliptic. (The input is equatorial.)
|
|
*/
|
|
if( direction == 1 )
|
|
eps = swi_epsiln(J, iflag); /* To J2000 */
|
|
else
|
|
eps = swi_epsiln(J2000, iflag); /* From J2000 */
|
|
sineps = sin(eps);
|
|
coseps = cos(eps);
|
|
x[0] = R[0];
|
|
z = coseps*R[1] + sineps*R[2];
|
|
x[2] = -sineps*R[1] + coseps*R[2];
|
|
x[1] = z;
|
|
/* Precession in longitude */
|
|
T /= 10.0; /* thousands of years */
|
|
p = pAcof;
|
|
pA = *p++;
|
|
for( i=0; i<9; i++ ) {
|
|
pA = pA * T + *p++;
|
|
}
|
|
pA *= DEGTORAD/3600 * T;
|
|
/* Node of the moving ecliptic on the J2000 ecliptic.
|
|
*/
|
|
p = nodecof;
|
|
W = *p++;
|
|
for( i=0; i<10; i++ )
|
|
W = W * T + *p++;
|
|
/* Rotate about z axis to the node.
|
|
*/
|
|
if( direction == 1 )
|
|
z = W + pA;
|
|
else
|
|
z = W;
|
|
B = cos(z);
|
|
A = sin(z);
|
|
z = B * x[0] + A * x[1];
|
|
x[1] = -A * x[0] + B * x[1];
|
|
x[0] = z;
|
|
/* Rotate about new x axis by the inclination of the moving
|
|
* ecliptic on the J2000 ecliptic.
|
|
*/
|
|
p = inclcof;
|
|
z = *p++;
|
|
for( i=0; i<10; i++ )
|
|
z = z * T + *p++;
|
|
if( direction == 1 )
|
|
z = -z;
|
|
B = cos(z);
|
|
A = sin(z);
|
|
z = B * x[1] + A * x[2];
|
|
x[2] = -A * x[1] + B * x[2];
|
|
x[1] = z;
|
|
/* Rotate about new z axis back from the node.
|
|
*/
|
|
if( direction == 1 )
|
|
z = -W;
|
|
else
|
|
z = -W - pA;
|
|
B = cos(z);
|
|
A = sin(z);
|
|
z = B * x[0] + A * x[1];
|
|
x[1] = -A * x[0] + B * x[1];
|
|
x[0] = z;
|
|
/* Rotate about x axis to final equator.
|
|
*/
|
|
if( direction == 1 )
|
|
eps = swi_epsiln(J2000, iflag);
|
|
else
|
|
eps = swi_epsiln(J, iflag);
|
|
sineps = sin(eps);
|
|
coseps = cos(eps);
|
|
z = coseps * x[1] - sineps * x[2];
|
|
x[2] = sineps * x[1] + coseps * x[2];
|
|
x[1] = z;
|
|
for( i=0; i<3; i++ )
|
|
R[i] = x[i];
|
|
return(0);
|
|
}
|
|
|
|
static int precess_3(double *R, double J, int direction, int iflag, int prec_meth)
|
|
{
|
|
//double T;
|
|
double x[3], pmat[9];
|
|
int i, j;
|
|
if( J == J2000 )
|
|
return(0);
|
|
/* Each precession angle is specified by a polynomial in
|
|
* T = Julian centuries from J2000.0. See AA page B18.
|
|
*/
|
|
//T = (J - J2000)/36525.0;
|
|
if (prec_meth == SEMOD_PREC_OWEN_1990)
|
|
owen_pre_matrix(J, pmat, iflag);
|
|
else
|
|
pre_pmat(J, pmat);
|
|
if (direction == -1) {
|
|
for (i = 0, j = 0; i <= 2; i++, j = i * 3) {
|
|
x[i] = R[0] * pmat[j + 0] +
|
|
R[1] * pmat[j + 1] +
|
|
R[2] * pmat[j + 2];
|
|
}
|
|
} else {
|
|
for (i = 0, j = 0; i <= 2; i++, j = i * 3) {
|
|
x[i] = R[0] * pmat[i + 0] +
|
|
R[1] * pmat[i + 3] +
|
|
R[2] * pmat[i + 6];
|
|
}
|
|
}
|
|
for (i = 0; i < 3; i++)
|
|
R[i] = x[i];
|
|
return(0);
|
|
}
|
|
|
|
/* Subroutine arguments:
|
|
*
|
|
* R = rectangular equatorial coordinate vector to be precessed.
|
|
* The result is written back into the input vector.
|
|
* J = Julian date
|
|
* direction =
|
|
* Precess from J to J2000: direction = 1
|
|
* Precess from J2000 to J: direction = -1
|
|
* Note that if you want to precess from J1 to J2, you would
|
|
* first go from J1 to J2000, then call the program again
|
|
* to go from J2000 to J2.
|
|
*/
|
|
int swi_precess(double *R, double J, int32 iflag, int direction )
|
|
{
|
|
double T = (J - J2000)/36525.0;
|
|
int prec_model = swed.astro_models[SE_MODEL_PREC_LONGTERM];
|
|
int prec_model_short = swed.astro_models[SE_MODEL_PREC_SHORTTERM];
|
|
int jplhora_model = swed.astro_models[SE_MODEL_JPLHORA_MODE];
|
|
AS_BOOL is_jplhor = FALSE;
|
|
if (prec_model == 0) prec_model = SEMOD_PREC_DEFAULT;
|
|
if (prec_model_short == 0) prec_model_short = SEMOD_PREC_DEFAULT_SHORT;
|
|
if (jplhora_model == 0) jplhora_model = SEMOD_JPLHORA_DEFAULT;
|
|
if (iflag & SEFLG_JPLHOR)
|
|
is_jplhor = TRUE;
|
|
if ((iflag & SEFLG_JPLHOR_APPROX)
|
|
&& jplhora_model == SEMOD_JPLHORA_3
|
|
&& J <= HORIZONS_TJD0_DPSI_DEPS_IAU1980)
|
|
is_jplhor = TRUE;
|
|
/* JPL Horizons uses precession IAU 1976 and nutation IAU 1980 plus
|
|
* some correction to nutation, arriving at extremely high precision */
|
|
if (is_jplhor) {
|
|
if (J > 2378131.5 && J < 2525323.5) { // between 1.1.1799 and 1.1.2202
|
|
return precess_1(R, J, direction, SEMOD_PREC_IAU_1976);
|
|
} else {
|
|
return precess_3(R, J, direction, iflag, SEMOD_PREC_OWEN_1990);
|
|
}
|
|
/* Use IAU 1976 formula for a few centuries. */
|
|
} else if (prec_model_short == SEMOD_PREC_IAU_1976 && fabs(T) <= PREC_IAU_1976_CTIES) {
|
|
return precess_1(R, J, direction, SEMOD_PREC_IAU_1976);
|
|
} else if (prec_model == SEMOD_PREC_IAU_1976) {
|
|
return precess_1(R, J, direction, SEMOD_PREC_IAU_1976);
|
|
/* Use IAU 2000 formula for a few centuries. */
|
|
} else if (prec_model_short == SEMOD_PREC_IAU_2000 && fabs(T) <= PREC_IAU_2000_CTIES) {
|
|
return precess_1(R, J, direction, SEMOD_PREC_IAU_2000);
|
|
} else if (prec_model == SEMOD_PREC_IAU_2000) {
|
|
return precess_1(R, J, direction, SEMOD_PREC_IAU_2000);
|
|
/* Use IAU 2006 formula for a few centuries. */
|
|
} else if (prec_model_short == SEMOD_PREC_IAU_2006 && fabs(T) <= PREC_IAU_2006_CTIES) {
|
|
return precess_1(R, J, direction, SEMOD_PREC_IAU_2006);
|
|
} else if (prec_model == SEMOD_PREC_IAU_2006) {
|
|
return precess_1(R, J, direction, SEMOD_PREC_IAU_2006);
|
|
} else if (prec_model == SEMOD_PREC_BRETAGNON_2003) {
|
|
return precess_1(R, J, direction, SEMOD_PREC_BRETAGNON_2003);
|
|
} else if (prec_model == SEMOD_PREC_NEWCOMB) {
|
|
return precess_1(R, J, direction, SEMOD_PREC_NEWCOMB);
|
|
} else if (prec_model == SEMOD_PREC_LASKAR_1986) {
|
|
return precess_2(R, J, iflag, direction, SEMOD_PREC_LASKAR_1986);
|
|
} else if (prec_model == SEMOD_PREC_SIMON_1994) {
|
|
return precess_2(R, J, iflag, direction, SEMOD_PREC_SIMON_1994);
|
|
} else if (prec_model == SEMOD_PREC_WILLIAMS_1994 || prec_model == SEMOD_PREC_WILL_EPS_LASK) {
|
|
return precess_2(R, J, iflag, direction, SEMOD_PREC_WILLIAMS_1994);
|
|
} else if (prec_model == SEMOD_PREC_OWEN_1990) {
|
|
return precess_3(R, J, direction, iflag, SEMOD_PREC_OWEN_1990);
|
|
} else { /* SEMOD_PREC_VONDRAK_2011 */
|
|
//int i;
|
|
//int retval = precess_3(R, J, direction, iflag, SEMOD_PREC_VONDRAK_2011);
|
|
// for (i = 0; i < 3; i++) printf("%.9f, ", R[i]);
|
|
return precess_3(R, J, direction, iflag, SEMOD_PREC_VONDRAK_2011);
|
|
}
|
|
}
|
|
|
|
/* 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
|
|
*
|
|
* - D. Koch, November 1995: small changes in structure,
|
|
* Corrections to IAU 1980 Series added from Expl. Suppl. p. 116
|
|
*
|
|
* 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 const short nt[] = {
|
|
/* LS and OC are units of 0.0001"
|
|
*LS2 and OC2 are units of 0.00001"
|
|
*MM,MS,FF,DD,OM, LS, LS2,OC, OC2 */
|
|
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,
|
|
#if 1
|
|
/*#if NUT_CORR_1987 switch is handled in function calc_nutation_iau1980() */
|
|
/* corrections to IAU 1980 nutation series by Herring 1987
|
|
* in 0.00001" !!!
|
|
* LS OC */
|
|
101, 0, 0, 0, 1,-725, 0, 213, 0,
|
|
101, 1, 0, 0, 0, 523, 0, 208, 0,
|
|
101, 0, 2,-2, 2, 102, 0, -41, 0,
|
|
101, 0, 2, 0, 2, -81, 0, 32, 0,
|
|
/* LC OS !!! */
|
|
102, 0, 0, 0, 1, 417, 0, 224, 0,
|
|
102, 1, 0, 0, 0, 61, 0, -24, 0,
|
|
102, 0, 2,-2, 2,-118, 0, -47, 0,
|
|
/*#endif*/
|
|
#endif
|
|
ENDMARK,
|
|
};
|
|
|
|
static int calc_nutation_iau1980(double J, double *nutlo)
|
|
{
|
|
/* arrays to hold sines and cosines of multiple angles */
|
|
double ss[5][8];
|
|
double cc[5][8];
|
|
double arg;
|
|
double args[5];
|
|
double f, g, T, T2;
|
|
double MM, MS, FF, DD, OM;
|
|
double cu, su, cv, sv, sw, s;
|
|
double C, D;
|
|
int i, j, k, k1, m, n;
|
|
int ns[5];
|
|
const short *p;
|
|
int nut_model = swed.astro_models[SE_MODEL_NUT];
|
|
if (nut_model == 0) nut_model = SEMOD_NUT_DEFAULT;
|
|
/* Julian centuries from 2000 January 1.5,
|
|
* barycentric dynamical time
|
|
*/
|
|
T = (J - 2451545.0) / 36525.0;
|
|
T2 = T * T;
|
|
/* 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 = -6962890.539 * T + 450160.280 + (0.008 * T + 7.455) * T2;
|
|
OM = swe_degnorm(OM/3600) * DEGTORAD;
|
|
/* mean longitude of the Sun minus the
|
|
* mean longitude of the Sun's perigee
|
|
*/
|
|
MS = 129596581.224 * T + 1287099.804 - (0.012 * T + 0.577) * T2;
|
|
MS = swe_degnorm(MS/3600) * DEGTORAD;
|
|
/* mean longitude of the Moon minus the
|
|
* mean longitude of the Moon's perigee
|
|
*/
|
|
MM = 1717915922.633 * T + 485866.733 + (0.064 * T + 31.310) * T2;
|
|
MM = swe_degnorm(MM/3600) * DEGTORAD;
|
|
/* mean longitude of the Moon minus the
|
|
* mean longitude of the Moon's node
|
|
*/
|
|
FF = 1739527263.137 * T + 335778.877 + (0.011 * T - 13.257) * T2;
|
|
FF = swe_degnorm(FF/3600) * DEGTORAD;
|
|
/* mean elongation of the Moon from the Sun.
|
|
*/
|
|
DD = 1602961601.328 * T + 1072261.307 + (0.019 * T - 6.891) * T2;
|
|
DD = swe_degnorm(DD/3600) * DEGTORAD;
|
|
args[0] = MM;
|
|
ns[0] = 3;
|
|
args[1] = MS;
|
|
ns[1] = 2;
|
|
args[2] = FF;
|
|
ns[2] = 4;
|
|
args[3] = DD;
|
|
ns[3] = 4;
|
|
args[4] = OM;
|
|
ns[4] = 2;
|
|
/* Calculate sin( i*MM ), etc. for needed multiple angles
|
|
*/
|
|
for (k = 0; k <= 4; k++) {
|
|
arg = args[k];
|
|
n = ns[k];
|
|
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;
|
|
}
|
|
}
|
|
/* 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) */
|
|
for(p = &nt[0]; *p != ENDMARK; p += 9) {
|
|
if (nut_model != SEMOD_NUT_IAU_CORR_1987 && (p[0] == 101 || p[0] == 102))
|
|
continue;
|
|
/* argument of sine and cosine */
|
|
k1 = 0;
|
|
cv = 0.0;
|
|
sv = 0.0;
|
|
for( m=0; m<5; m++ ) {
|
|
j = p[m];
|
|
if (j > 100)
|
|
j = 0; /* p[0] is a flag */
|
|
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, in 0.0001" */
|
|
f = p[5] * 0.0001;
|
|
if( p[6] != 0 )
|
|
f += 0.00001 * T * p[6];
|
|
/* obliquity coefficient, in 0.0001" */
|
|
g = p[7] * 0.0001;
|
|
if( p[8] != 0 )
|
|
g += 0.00001 * T * p[8];
|
|
if (*p >= 100) { /* coefficients in 0.00001" */
|
|
f *= 0.1;
|
|
g *= 0.1;
|
|
}
|
|
/* accumulate the terms */
|
|
if (*p != 102) {
|
|
C += f * sv;
|
|
D += g * cv;
|
|
}
|
|
else { /* cos for nutl and sin for nuto */
|
|
C += f * cv;
|
|
D += g * sv;
|
|
}
|
|
/*
|
|
if (i >= 105) {
|
|
printf("%4.10f, %4.10f\n",f*sv,g*cv);
|
|
}
|
|
*/
|
|
}
|
|
/*
|
|
printf("%4.10f, %4.10f, %4.10f, %4.10f\n",MS*RADTODEG,FF*RADTODEG,DD*RADTODEG,OM*RADTODEG);
|
|
printf( "nutation: in longitude %.9f\", in obliquity %.9f\"\n", C, D );
|
|
*/
|
|
/* Save answers, expressed in radians */
|
|
nutlo[0] = DEGTORAD * C / 3600.0;
|
|
nutlo[1] = DEGTORAD * D / 3600.0;
|
|
return(0);
|
|
}
|
|
|
|
/* Nutation IAU 2000A model
|
|
* (MHB2000 luni-solar and planetary nutation, without free core nutation)
|
|
*
|
|
* Function returns nutation in longitude and obliquity in radians with
|
|
* respect to the equinox of date. For the obliquity of the ecliptic
|
|
* the calculation of Lieske & al. (1977) must be used.
|
|
*
|
|
* The precision in recent years is about 0.001 arc seconds.
|
|
*
|
|
* The calculation includes luni-solar and planetary nutation.
|
|
* Free core nutation, which cannot be predicted, is omitted,
|
|
* the error being of the order of a few 0.0001 arc seconds.
|
|
*
|
|
* References:
|
|
*
|
|
* Capitaine, N., Wallace, P.T., Chapront, J., A & A 432, 366 (2005).
|
|
*
|
|
* Chapront, J., Chapront-Touze, M. & Francou, G., A & A 387, 700 (2002).
|
|
*
|
|
* Lieske, J.H., Lederle, T., Fricke, W. & Morando, B., "Expressions
|
|
* for the precession quantities based upon the IAU (1976) System of
|
|
* Astronomical Constants", A & A 58, 1-16 (1977).
|
|
*
|
|
* Mathews, P.M., Herring, T.A., Buffet, B.A., "Modeling of nutation
|
|
* and precession New nutation series for nonrigid Earth and
|
|
* insights into the Earth's interior", J.Geophys.Res., 107, B4,
|
|
* 2002.
|
|
*
|
|
* Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
|
|
* Francou, G., Laskar, J., A & A 282, 663-683 (1994).
|
|
*
|
|
* Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M., A & A Supp.
|
|
* Ser. 135, 111 (1999).
|
|
*
|
|
* Wallace, P.T., "Software for Implementing the IAU 2000
|
|
* Resolutions", in IERS Workshop 5.1 (2002).
|
|
*
|
|
* Nutation IAU 2000A series in:
|
|
* Kaplan, G.H., United States Naval Observatory Circular No. 179 (Oct. 2005)
|
|
* aa.usno.navy.mil/publications/docs/Circular_179.html
|
|
*
|
|
* MHB2000 code at
|
|
* - ftp://maia.usno.navy.mil/conv2000/chapter5/IAU2000A.
|
|
* - http://www.iau-sofa.rl.ac.uk/2005_0901/Downloads.html
|
|
*/
|
|
|
|
#include "swenut2000a.h"
|
|
static int calc_nutation_iau2000ab(double J, double *nutlo)
|
|
{
|
|
int i, j, k, inls;
|
|
double M, SM, F, D, OM;
|
|
double AL, ALSU, AF, AD, AOM, APA;
|
|
double ALME, ALVE, ALEA, ALMA, ALJU, ALSA, ALUR, ALNE;
|
|
double darg, sinarg, cosarg;
|
|
double dpsi = 0, deps = 0;
|
|
double T = (J - J2000 ) / 36525.0;
|
|
int nut_model = swed.astro_models[SE_MODEL_NUT];
|
|
if (nut_model == 0) nut_model = SEMOD_NUT_DEFAULT;
|
|
/* luni-solar nutation */
|
|
/* Fundamental arguments, Simon & al. (1994) */
|
|
/* Mean anomaly of the Moon. */
|
|
M = swe_degnorm(( 485868.249036 +
|
|
T*( 1717915923.2178 +
|
|
T*( 31.8792 +
|
|
T*( 0.051635 +
|
|
T*( - 0.00024470 ))))) / 3600.0) * DEGTORAD;
|
|
/* Mean anomaly of the Sun */
|
|
SM = swe_degnorm((1287104.79305 +
|
|
T*( 129596581.0481 +
|
|
T*( - 0.5532 +
|
|
T*( 0.000136 +
|
|
T*( - 0.00001149 ))))) / 3600.0) * DEGTORAD;
|
|
/* Mean argument of the latitude of the Moon. */
|
|
F = swe_degnorm(( 335779.526232 +
|
|
T*( 1739527262.8478 +
|
|
T*( - 12.7512 +
|
|
T*( - 0.001037 +
|
|
T*( 0.00000417 ))))) / 3600.0) * DEGTORAD;
|
|
/* Mean elongation of the Moon from the Sun. */
|
|
D = swe_degnorm((1072260.70369 +
|
|
T*( 1602961601.2090 +
|
|
T*( - 6.3706 +
|
|
T*( 0.006593 +
|
|
T*( - 0.00003169 ))))) / 3600.0) * DEGTORAD;
|
|
/* Mean longitude of the ascending node of the Moon. */
|
|
OM = swe_degnorm(( 450160.398036 +
|
|
T*( - 6962890.5431 +
|
|
T*( 7.4722 +
|
|
T*( 0.007702 +
|
|
T*( - 0.00005939 ))))) / 3600.0) * DEGTORAD;
|
|
/* luni-solar nutation series, in reverse order, starting with small terms */
|
|
if (nut_model == SEMOD_NUT_IAU_2000B)
|
|
inls = NLS_2000B;
|
|
else
|
|
inls = NLS;
|
|
for (i = inls - 1; i >= 0; i--) {
|
|
j = i * 5;
|
|
darg = swe_radnorm((double) nls[j + 0] * M +
|
|
(double) nls[j + 1] * SM +
|
|
(double) nls[j + 2] * F +
|
|
(double) nls[j + 3] * D +
|
|
(double) nls[j + 4] * OM);
|
|
sinarg = sin(darg);
|
|
cosarg = cos(darg);
|
|
k = i * 6;
|
|
dpsi += (cls[k+0] + cls[k+1] * T) * sinarg + cls[k+2] * cosarg;
|
|
deps += (cls[k+3] + cls[k+4] * T) * cosarg + cls[k+5] * sinarg;
|
|
}
|
|
nutlo[0] = dpsi * O1MAS2DEG;
|
|
nutlo[1] = deps * O1MAS2DEG;
|
|
if (nut_model == SEMOD_NUT_IAU_2000A) {
|
|
/* planetary nutation
|
|
* note: The MHB2000 code computes the luni-solar and planetary nutation
|
|
* in different routines, using slightly different Delaunay
|
|
* arguments in the two cases. This behaviour is faithfully
|
|
* reproduced here. Use of the Simon et al. expressions for both
|
|
* cases leads to negligible changes, well below 0.1 microarcsecond.*/
|
|
/* Mean anomaly of the Moon.*/
|
|
AL = swe_radnorm(2.35555598 + 8328.6914269554 * T);
|
|
/* Mean anomaly of the Sun.*/
|
|
ALSU = swe_radnorm(6.24006013 + 628.301955 * T);
|
|
/* Mean argument of the latitude of the Moon. */
|
|
AF = swe_radnorm(1.627905234 + 8433.466158131 * T);
|
|
/* Mean elongation of the Moon from the Sun. */
|
|
AD = swe_radnorm(5.198466741 + 7771.3771468121 * T);
|
|
/* Mean longitude of the ascending node of the Moon. */
|
|
AOM = swe_radnorm(2.18243920 - 33.757045 * T);
|
|
/* Planetary longitudes, Mercury through Neptune (Souchay et al. 1999). */
|
|
ALME = swe_radnorm(4.402608842 + 2608.7903141574 * T);
|
|
ALVE = swe_radnorm(3.176146697 + 1021.3285546211 * T);
|
|
ALEA = swe_radnorm(1.753470314 + 628.3075849991 * T);
|
|
ALMA = swe_radnorm(6.203480913 + 334.0612426700 * T);
|
|
ALJU = swe_radnorm(0.599546497 + 52.9690962641 * T);
|
|
ALSA = swe_radnorm(0.874016757 + 21.3299104960 * T);
|
|
ALUR = swe_radnorm(5.481293871 + 7.4781598567 * T);
|
|
ALNE = swe_radnorm(5.321159000 + 3.8127774000 * T);
|
|
/* General accumulated precession in longitude. */
|
|
APA = (0.02438175 + 0.00000538691 * T) * T;
|
|
/* planetary nutation series (in reverse order).*/
|
|
dpsi = 0;
|
|
deps = 0;
|
|
for (i = NPL - 1; i >= 0; i--) {
|
|
j = i * 14;
|
|
darg = swe_radnorm((double) npl[j + 0] * AL +
|
|
(double) npl[j + 1] * ALSU +
|
|
(double) npl[j + 2] * AF +
|
|
(double) npl[j + 3] * AD +
|
|
(double) npl[j + 4] * AOM +
|
|
(double) npl[j + 5] * ALME +
|
|
(double) npl[j + 6] * ALVE +
|
|
(double) npl[j + 7] * ALEA +
|
|
(double) npl[j + 8] * ALMA +
|
|
(double) npl[j + 9] * ALJU +
|
|
(double) npl[j +10] * ALSA +
|
|
(double) npl[j +11] * ALUR +
|
|
(double) npl[j +12] * ALNE +
|
|
(double) npl[j +13] * APA);
|
|
k = i * 4;
|
|
sinarg = sin(darg);
|
|
cosarg = cos(darg);
|
|
dpsi += (double) icpl[k+0] * sinarg + (double) icpl[k+1] * cosarg;
|
|
deps += (double) icpl[k+2] * sinarg + (double) icpl[k+3] * cosarg;
|
|
}
|
|
nutlo[0] += dpsi * O1MAS2DEG;
|
|
nutlo[1] += deps * O1MAS2DEG;
|
|
#if 1
|
|
/* changes required by adoption of P03 precession
|
|
* according to Capitaine et al. A & A 412, 366 (2005) = IAU 2006 */
|
|
dpsi = -8.1 * sin(OM) - 0.6 * sin(2 * F - 2 * D + 2 * OM);
|
|
dpsi += T * (47.8 * sin(OM) + 3.7 * sin(2 * F - 2 * D + 2 * OM) + 0.6 * sin(2 * F + 2 * OM) - 0.6 * sin(2 * OM));
|
|
deps = T * (-25.6 * cos(OM) - 1.6 * cos(2 * F - 2 * D + 2 * OM));
|
|
nutlo[0] += dpsi / (3600.0 * 1000000.0);
|
|
nutlo[1] += deps / (3600.0 * 1000000.0);
|
|
#endif
|
|
}
|
|
nutlo[0] *= DEGTORAD;
|
|
nutlo[1] *= DEGTORAD;
|
|
return 0;
|
|
}
|
|
|
|
/* an incomplete implementation of nutation Woolard 1953 */
|
|
static int calc_nutation_woolard(double J, double *nutlo)
|
|
{
|
|
double deps, dpsi;
|
|
double ls, ld; /* sun's mean longitude, moon's mean longitude */
|
|
double ms, md; /* sun's mean anomaly, moon's mean anomaly */
|
|
double nm; /* longitude of moon's ascending node */
|
|
double t, t2; /* number of Julian centuries of 36525 days since
|
|
* Jan 0.5 1900.
|
|
*/
|
|
double tls, tnm, tld; /* twice above */
|
|
double a, b; /* temps */
|
|
double mjd = J - J1900;
|
|
t = mjd/36525.;
|
|
t2 = t*t;
|
|
a = 100.0021358*t;
|
|
b = 360.*(a-(long)a);
|
|
ls = 279.697+.000303*t2+b;
|
|
a = 1336.855231*t;
|
|
b = 360.*(a-(long)a);
|
|
ld = 270.434-.001133*t2+b;
|
|
a = 99.99736056000026*t;
|
|
b = 360.*(a-(long)a);
|
|
ms = 358.476-.00015*t2+b;
|
|
a = 13255523.59*t;
|
|
b = 360.*(a-(long)a);
|
|
md = 296.105+.009192*t2+b;
|
|
a = 5.372616667*t;
|
|
b = 360.*(a-(long)a);
|
|
nm = 259.183+.002078*t2-b;
|
|
/* convert to radian forms for use with trig functions.
|
|
*/
|
|
tls = 2*ls * DEGTORAD;
|
|
nm = nm * DEGTORAD;
|
|
tnm = 2*nm;
|
|
ms = ms * DEGTORAD;
|
|
tld = 2*ld * DEGTORAD;
|
|
md = md * DEGTORAD;
|
|
/* find delta psi and eps, in arcseconds.
|
|
*/
|
|
dpsi = (-17.2327-.01737*t)*sin(nm)+(-1.2729-.00013*t)*sin(tls)
|
|
+.2088*sin(tnm)-.2037*sin(tld)+(.1261-.00031*t)*sin(ms)
|
|
+.0675*sin(md)-(.0497-.00012*t)*sin(tls+ms)
|
|
-.0342*sin(tld-nm)-.0261*sin(tld+md)+.0214*sin(tls-ms)
|
|
-.0149*sin(tls-tld+md)+.0124*sin(tls-nm)+.0114*sin(tld-md);
|
|
deps = (9.21+.00091*t)*cos(nm)+(.5522-.00029*t)*cos(tls)
|
|
-.0904*cos(tnm)+.0884*cos(tld)+.0216*cos(tls+ms)
|
|
+.0183*cos(tld-nm)+.0113*cos(tld+md)-.0093*cos(tls-ms)
|
|
-.0066*cos(tls-nm);
|
|
/* convert to radians.
|
|
*/
|
|
dpsi = dpsi/3600.0 * DEGTORAD;
|
|
deps = deps/3600.0 * DEGTORAD;
|
|
nutlo[1] = deps;
|
|
nutlo[0] = dpsi;
|
|
return OK;
|
|
}
|
|
|
|
static double bessel(double *v, int n, double t)
|
|
{
|
|
int i, iy, k;
|
|
double ans, p, B, d[6];
|
|
if (t <= 0) {
|
|
ans = v[0];
|
|
goto done;
|
|
}
|
|
if (t >= n - 1) {
|
|
ans = v[n - 1];
|
|
goto done;
|
|
}
|
|
p = floor(t);
|
|
iy = (int) t;
|
|
/* Zeroth order estimate is value at start of year */
|
|
ans = v[iy];
|
|
k = iy + 1;
|
|
if (k >= n)
|
|
goto done;
|
|
/* The fraction of tabulation interval */
|
|
p = t - p;
|
|
ans += p * (v[k] - v[iy]);
|
|
if( (iy - 1 < 0) || (iy + 2 >= n) )
|
|
goto done; /* can't do second differences */
|
|
/* Make table of first differences */
|
|
k = iy - 2;
|
|
for (i = 0; i < 5; i++) {
|
|
if((k < 0) || (k + 1 >= n))
|
|
d[i] = 0;
|
|
else
|
|
d[i] = v[k+1] - v[k];
|
|
k += 1;
|
|
}
|
|
/* Compute second differences */
|
|
for (i = 0; i < 4; i++ )
|
|
d[i] = d[i+1] - d[i];
|
|
B = 0.25 * p * (p - 1.0);
|
|
ans += B * (d[1] + d[2]);
|
|
#if DEMO
|
|
printf("B %.4lf, ans %.4lf\n", B, ans);
|
|
#endif
|
|
if (iy + 2 >= n)
|
|
goto done;
|
|
/* Compute third differences */
|
|
for (i = 0; i < 3; i++ )
|
|
d[i] = d[i + 1] - d[i];
|
|
B = 2.0 * B / 3.0;
|
|
ans += (p - 0.5) * B * d[1];
|
|
#if DEMO
|
|
printf("B %.4lf, ans %.4lf\n", B * (p - 0.5), ans);
|
|
#endif
|
|
if ((iy - 2 < 0) || (iy + 3 > n))
|
|
goto done;
|
|
/* Compute fourth differences */
|
|
for (i = 0; i < 2; i++)
|
|
d[i] = d[i + 1] - d[i];
|
|
B = 0.125 * B * (p + 1.0) * (p - 2.0);
|
|
ans += B * (d[0] + d[1]);
|
|
#if DEMO
|
|
printf("B %.4lf, ans %.4lf\n", B, ans);
|
|
#endif
|
|
done:
|
|
return ans;
|
|
}
|
|
|
|
static int calc_nutation(double J, int32 iflag, double *nutlo)
|
|
{
|
|
int n;
|
|
double dpsi, deps, J2;
|
|
int nut_model = swed.astro_models[SE_MODEL_NUT];
|
|
int jplhora_model = swed.astro_models[SE_MODEL_JPLHORA_MODE];
|
|
AS_BOOL is_jplhor = FALSE;
|
|
if (nut_model == 0) nut_model = SEMOD_NUT_DEFAULT;
|
|
if (jplhora_model == 0) jplhora_model = SEMOD_JPLHORA_DEFAULT;
|
|
if (iflag & SEFLG_JPLHOR)
|
|
is_jplhor = TRUE;
|
|
if ((iflag & SEFLG_JPLHOR_APPROX) &&
|
|
jplhora_model == SEMOD_JPLHORA_3
|
|
&& J <= HORIZONS_TJD0_DPSI_DEPS_IAU1980)
|
|
is_jplhor = TRUE;
|
|
if (is_jplhor) {
|
|
calc_nutation_iau1980(J, nutlo);
|
|
if (iflag & SEFLG_JPLHOR) {
|
|
n = (int) (swed.eop_tjd_end - swed.eop_tjd_beg + 0.000001);
|
|
J2 = J;
|
|
if (J < swed.eop_tjd_beg_horizons)
|
|
J2 = swed.eop_tjd_beg_horizons;
|
|
dpsi = bessel(swed.dpsi, n + 1, J2 - swed.eop_tjd_beg);
|
|
deps = bessel(swed.deps, n + 1, J2 - swed.eop_tjd_beg);
|
|
nutlo[0] += dpsi / 3600.0 * DEGTORAD;
|
|
nutlo[1] += deps / 3600.0 * DEGTORAD;
|
|
#if 0
|
|
printf("tjd=%f, dpsi=%f, deps=%f\n", J, dpsi * 1000, deps * 1000);
|
|
#endif
|
|
} else {
|
|
nutlo[0] += DPSI_IAU1980_TJD0 / 3600.0 * DEGTORAD;
|
|
nutlo[1] += DEPS_IAU1980_TJD0 / 3600.0 * DEGTORAD;
|
|
}
|
|
} else if (nut_model == SEMOD_NUT_IAU_1980 || nut_model == SEMOD_NUT_IAU_CORR_1987) {
|
|
calc_nutation_iau1980(J, nutlo);
|
|
} else if (nut_model == SEMOD_NUT_IAU_2000A || nut_model == SEMOD_NUT_IAU_2000B) {
|
|
calc_nutation_iau2000ab(J, nutlo);
|
|
if ((iflag & SEFLG_JPLHOR_APPROX) && jplhora_model == SEMOD_JPLHORA_2) {
|
|
nutlo[0] += -41.7750 / 3600.0 / 1000.0 * DEGTORAD;
|
|
nutlo[1] += -6.8192 / 3600.0 / 1000.0 * DEGTORAD;
|
|
}
|
|
} else if (nut_model == SEMOD_NUT_WOOLARD) {
|
|
calc_nutation_woolard(J, nutlo);
|
|
}
|
|
return OK;
|
|
}
|
|
|
|
static double quadratic_intp(double ym, double y0, double yp, double x)
|
|
{
|
|
double a, b, c, y;
|
|
c = y0;
|
|
b = (yp - ym) / 2.0;
|
|
a = (yp + ym) / 2.0 - c;
|
|
y = a * x * x + b * x + c;
|
|
return y;
|
|
}
|
|
|
|
int swi_nutation(double tjd, int32 iflag, double *nutlo)
|
|
{
|
|
int retc = OK;
|
|
double dnut[2], dx;
|
|
if (!swed.do_interpolate_nut) {
|
|
retc = calc_nutation(tjd, iflag, nutlo);
|
|
// from interpolation, with three data points in 1-day steps;
|
|
// maximum error is about 3 mas
|
|
} else {
|
|
// precalculated data points available
|
|
if (tjd < swed.interpol.tjd_nut2 && tjd > swed.interpol.tjd_nut0) {
|
|
dx = (tjd - swed.interpol.tjd_nut0) - 1.0;
|
|
nutlo[0] = quadratic_intp(swed.interpol.nut_dpsi0, swed.interpol.nut_dpsi1, swed.interpol.nut_dpsi2, dx);
|
|
nutlo[1] = quadratic_intp(swed.interpol.nut_deps0, swed.interpol.nut_deps1, swed.interpol.nut_deps2, dx);
|
|
} else {
|
|
swed.interpol.tjd_nut0 = tjd - 1.0; // one day earlier
|
|
swed.interpol.tjd_nut2 = tjd + 1.0; // one day later
|
|
retc = calc_nutation(swed.interpol.tjd_nut0, iflag, dnut);
|
|
if (retc == ERR) return ERR;
|
|
swed.interpol.nut_dpsi0 = dnut[0];
|
|
swed.interpol.nut_deps0 = dnut[1];
|
|
retc = calc_nutation(swed.interpol.tjd_nut2, iflag, dnut);
|
|
if (retc == ERR) return ERR;
|
|
swed.interpol.nut_dpsi2 = dnut[0];
|
|
swed.interpol.nut_deps2 = dnut[1];
|
|
retc = calc_nutation(tjd, iflag, nutlo);
|
|
if (retc == ERR) return ERR;
|
|
swed.interpol.nut_dpsi1 = nutlo[0];
|
|
swed.interpol.nut_deps1 = nutlo[1];
|
|
}
|
|
}
|
|
return retc;
|
|
}
|
|
|
|
#define OFFSET_JPLHORIZONS (-52.3)
|
|
#define DCOR_RA_JPL_TJD0 2437846.5
|
|
#define NDCOR_RA_JPL 51
|
|
double dcor_ra_jpl[] = {
|
|
-51.257, -51.103, -51.065, -51.503, -51.224, -50.796, -51.161, -51.181,
|
|
-50.932, -51.064, -51.182, -51.386, -51.416, -51.428, -51.586, -51.766, -52.038, -52.370,
|
|
-52.553, -52.397, -52.340, -52.676, -52.348, -51.964, -52.444, -52.364, -51.988, -52.212,
|
|
-52.370, -52.523, -52.541, -52.496, -52.590, -52.629, -52.788, -53.014, -53.053, -52.902,
|
|
-52.850, -53.087, -52.635, -52.185, -52.588, -52.292, -51.796, -51.961, -52.055, -52.134,
|
|
-52.165, -52.141, -52.255,
|
|
};
|
|
|
|
static void swi_approx_jplhor(double *x, double tjd, int32 iflag, AS_BOOL backward)
|
|
{
|
|
double t0, t1;
|
|
double t = (tjd - DCOR_RA_JPL_TJD0) / 365.25;
|
|
double dofs = OFFSET_JPLHORIZONS;
|
|
int jplhora_model = swed.astro_models[SE_MODEL_JPLHORA_MODE];
|
|
if (jplhora_model == 0) jplhora_model = SEMOD_JPLHORA_DEFAULT;
|
|
if (!(iflag & SEFLG_JPLHOR_APPROX))
|
|
return;
|
|
if (jplhora_model == SEMOD_JPLHORA_2)
|
|
return;
|
|
if (t < 0) {
|
|
t = 0;
|
|
dofs = dcor_ra_jpl[0];
|
|
} else if (t >= NDCOR_RA_JPL - 1) {
|
|
t = NDCOR_RA_JPL;
|
|
dofs = dcor_ra_jpl[NDCOR_RA_JPL - 1];
|
|
} else {
|
|
t0 = (int) t;
|
|
t1 = t0 + 1;
|
|
dofs = dcor_ra_jpl[(int)t0];
|
|
dofs = (t - t0) * (dcor_ra_jpl[(int)t0] - dcor_ra_jpl[(int)t1]) + dcor_ra_jpl[(int)t0];
|
|
}
|
|
dofs /= (1000.0 * 3600.0);
|
|
swi_cartpol(x, x);
|
|
if (backward)
|
|
x[0] -= dofs * DEGTORAD;
|
|
else
|
|
x[0] += dofs * DEGTORAD;
|
|
swi_polcart(x, x);
|
|
}
|
|
|
|
/* GCRS to J2000 */
|
|
void swi_bias(double *x, double tjd, int32 iflag, AS_BOOL backward)
|
|
{
|
|
#if 0
|
|
double DAS2R = 1.0 / 3600.0 * DEGTORAD;
|
|
double dpsi_bias = -0.041775 * DAS2R;
|
|
double deps_bias = -0.0068192 * DAS2R;
|
|
double dra0 = -0.0146 * DAS2R;
|
|
double deps2000 = 84381.448 * DAS2R;
|
|
#endif
|
|
double xx[6], rb[3][3];
|
|
int i;
|
|
int bias_model = swed.astro_models[SE_MODEL_BIAS];
|
|
int jplhora_model = swed.astro_models[SE_MODEL_JPLHORA_MODE];
|
|
if (bias_model == 0) bias_model = SEMOD_BIAS_DEFAULT;
|
|
if (jplhora_model == 0) jplhora_model = SEMOD_JPLHORA_DEFAULT;
|
|
if (bias_model == SEMOD_BIAS_NONE)
|
|
return;
|
|
if (iflag & SEFLG_JPLHOR_APPROX) {
|
|
if (jplhora_model == SEMOD_JPLHORA_2)
|
|
return;
|
|
if (jplhora_model == SEMOD_JPLHORA_3 && tjd < DPSI_DEPS_IAU1980_TJD0_HORIZONS)
|
|
return;
|
|
}
|
|
/* #if FRAME_BIAS_IAU2006 * frame bias 2006 */
|
|
if (bias_model == SEMOD_BIAS_IAU2006) {
|
|
rb[0][0] = +0.99999999999999412;
|
|
rb[1][0] = -0.00000007078368961;
|
|
rb[2][0] = +0.00000008056213978;
|
|
rb[0][1] = +0.00000007078368695;
|
|
rb[1][1] = +0.99999999999999700;
|
|
rb[2][1] = +0.00000003306428553;
|
|
rb[0][2] = -0.00000008056214212;
|
|
rb[1][2] = -0.00000003306427981;
|
|
rb[2][2] = +0.99999999999999634;
|
|
/* #else * frame bias 2000, makes no difference in result */
|
|
} else {
|
|
rb[0][0] = +0.9999999999999942;
|
|
rb[1][0] = -0.0000000707827974;
|
|
rb[2][0] = +0.0000000805621715;
|
|
rb[0][1] = +0.0000000707827948;
|
|
rb[1][1] = +0.9999999999999969;
|
|
rb[2][1] = +0.0000000330604145;
|
|
rb[0][2] = -0.0000000805621738;
|
|
rb[1][2] = -0.0000000330604088;
|
|
rb[2][2] = +0.9999999999999962;
|
|
}
|
|
/*#endif*/
|
|
#if 0
|
|
rb[0][0] = +0.9999999999999968;
|
|
rb[1][0] = +0.0000000000000000;
|
|
rb[2][0] = +0.0000000805621715;
|
|
rb[0][1] = -0.0000000000000027;
|
|
rb[1][1] = +0.9999999999999994;
|
|
rb[2][1] = +0.0000000330604145;
|
|
rb[0][2] = -0.0000000805621715;
|
|
rb[1][2] = -0.0000000330604145;
|
|
rb[2][2] = +0.9999999999999962;
|
|
#endif
|
|
if (backward) {
|
|
swi_approx_jplhor(x, tjd, iflag, TRUE);
|
|
for (i = 0; i <= 2; i++) {
|
|
xx[i] = x[0] * rb[i][0] +
|
|
x[1] * rb[i][1] +
|
|
x[2] * rb[i][2];
|
|
if (iflag & SEFLG_SPEED)
|
|
xx[i+3] = x[3] * rb[i][0] +
|
|
x[4] * rb[i][1] +
|
|
x[5] * rb[i][2];
|
|
}
|
|
} else {
|
|
for (i = 0; i <= 2; i++) {
|
|
xx[i] = x[0] * rb[0][i] +
|
|
x[1] * rb[1][i] +
|
|
x[2] * rb[2][i];
|
|
if (iflag & SEFLG_SPEED)
|
|
xx[i+3] = x[3] * rb[0][i] +
|
|
x[4] * rb[1][i] +
|
|
x[5] * rb[2][i];
|
|
}
|
|
swi_approx_jplhor(xx, tjd, iflag, FALSE);
|
|
}
|
|
for (i = 0; i <= 2; i++) x[i] = xx[i];
|
|
if (iflag & SEFLG_SPEED)
|
|
for (i = 3; i <= 5; i++) x[i] = xx[i];
|
|
}
|
|
|
|
/* GCRS to FK5 */
|
|
void swi_icrs2fk5(double *x, int32 iflag, AS_BOOL backward)
|
|
{
|
|
#if 0
|
|
double DAS2R = 1.0 / 3600.0 * DEGTORAD;
|
|
double dra0 = -0.0229 * DAS2R;
|
|
double dxi0 = 0.0091 * DAS2R;
|
|
double det0 = -0.0199 * DAS2R;
|
|
#endif
|
|
double xx[6], rb[3][3];
|
|
int i;
|
|
rb[0][0] = +0.9999999999999928;
|
|
rb[0][1] = +0.0000001110223287;
|
|
rb[0][2] = +0.0000000441180557;
|
|
rb[1][0] = -0.0000001110223330;
|
|
rb[1][1] = +0.9999999999999891;
|
|
rb[1][2] = +0.0000000964779176;
|
|
rb[2][0] = -0.0000000441180450;
|
|
rb[2][1] = -0.0000000964779225;
|
|
rb[2][2] = +0.9999999999999943;
|
|
if (backward) {
|
|
for (i = 0; i <= 2; i++) {
|
|
xx[i] = x[0] * rb[i][0] +
|
|
x[1] * rb[i][1] +
|
|
x[2] * rb[i][2];
|
|
if (iflag & SEFLG_SPEED)
|
|
xx[i+3] = x[3] * rb[i][0] +
|
|
x[4] * rb[i][1] +
|
|
x[5] * rb[i][2];
|
|
}
|
|
} else {
|
|
for (i = 0; i <= 2; i++) {
|
|
xx[i] = x[0] * rb[0][i] +
|
|
x[1] * rb[1][i] +
|
|
x[2] * rb[2][i];
|
|
if (iflag & SEFLG_SPEED)
|
|
xx[i+3] = x[3] * rb[0][i] +
|
|
x[4] * rb[1][i] +
|
|
x[5] * rb[2][i];
|
|
}
|
|
}
|
|
for (i = 0; i <= 5; i++) x[i] = xx[i];
|
|
}
|
|
|
|
/* DeltaT = Ephemeris Time - Universal Time, in days.
|
|
*
|
|
* Before 1955 we use the data developed by
|
|
* Stephenson, Morrison, and Hohenkerk (2016),
|
|
*
|
|
* 1955 - today + a couple of years:
|
|
* ---------------------------------
|
|
* The tabulated values of deltaT from the Astronomical
|
|
* Alamanc (AA 1997 etc. pp. K8-K9) are used. Some
|
|
* more recent values have been taken from IERS
|
|
* (http://maia.usno.navy.mil/ser7/deltat.data).
|
|
* Bessel's interpolation formula is implemented to obtain fourth
|
|
* order interpolated values at intermediate times.
|
|
* The values are adjusted depending on the ephemeris used
|
|
* and its inherent value of secular tidal acceleration ndot.
|
|
*
|
|
* future:
|
|
* ---------------------------------
|
|
* For the time after the last tabulated value, we use the formula
|
|
* of Stephenson (1997; p. 507), with a modification that avoids a jump
|
|
* at the end of the tabulated period. A linear term is added that
|
|
* makes a slow transition from the table to the formula over a period
|
|
* of 100 years. (Need not be updated, when table will be enlarged.)
|
|
*
|
|
* References:
|
|
*
|
|
* Stephenson, F. R., and L. V. Morrison, "Long-term changes
|
|
* in the rotation of the Earth: 700 B.C. to A.D. 1980,"
|
|
* Philosophical Transactions of the Royal Society of London
|
|
* Series A 313, 47-70 (1984)
|
|
*
|
|
* Borkowski, K. M., "ELP2000-85 and the Dynamical Time
|
|
* - Universal Time relation," Astronomy and Astrophysics
|
|
* 205, L8-L10 (1988)
|
|
* Borkowski's formula is derived from partly doubtful eclipses
|
|
* going back to 2137 BC and uses lunar position based on tidal
|
|
* coefficient of -23.9 arcsec/cy^2.
|
|
*
|
|
* Chapront-Touze, Michelle, and Jean Chapront, _Lunar Tables
|
|
* and Programs from 4000 B.C. to A.D. 8000_, Willmann-Bell 1991
|
|
* Their table agrees with the one here, but the entries are
|
|
* rounded to the nearest whole second.
|
|
*
|
|
* Stephenson, F. R., and M. A. Houlden, _Atlas of Historical
|
|
* Eclipse Maps_, Cambridge U. Press (1986)
|
|
*
|
|
* Stephenson, F.R. & Morrison, L.V., "Long-Term Fluctuations in
|
|
* the Earth's Rotation: 700 BC to AD 1990", Philosophical
|
|
* Transactions of the Royal Society of London,
|
|
* Ser. A, 351 (1995), 165-202.
|
|
*
|
|
* Stephenson, F. Richard, _Historical Eclipses and Earth's
|
|
* Rotation_, Cambridge U. Press (1997)
|
|
*
|
|
* Morrison, L. V., and F.R. Stephenson, "Historical Values of the Earth's
|
|
* Clock Error DT and the Calculation of Eclipses", JHA xxxv (2004),
|
|
* pp.327-336
|
|
*
|
|
* Stephenson, F.R., Morrison, L.V., and Hohenkerk, C.Y., "Measurement of the
|
|
* Earth's Rotation: 720 BC to AD 2015", Royal Society Proceedings A
|
|
* 7 Dec 2016,
|
|
* http://rspa.royalsocietypublishing.org/lookup/doi/10.1098/rspa.2016.0404
|
|
*
|
|
* Table from AA for 1620 through today
|
|
* Note, Stephenson and Morrison's table starts at the year 1630.
|
|
* The Chapronts' table does not agree with the Almanac prior to 1630.
|
|
* The actual accuracy decreases rapidly prior to 1780.
|
|
*
|
|
* Jean Meeus, Astronomical Algorithms, 2nd edition, 1998.
|
|
*
|
|
* For a comprehensive collection of publications and formulae, see:
|
|
* http://www.phys.uu.nl/~vgent/deltat/deltat_modern.htm
|
|
* http://www.phys.uu.nl/~vgent/deltat/deltat_old.htm
|
|
*
|
|
* For future values of delta t, the following data from the
|
|
* Earth Orientation Department of the US Naval Observatory can be used:
|
|
* (TAI-UTC) from: ftp://maia.usno.navy.mil/ser7/tai-utc.dat
|
|
* (UT1-UTC) from: ftp://maia.usno.navy.mil/ser7/finals.all (cols. 59-68)
|
|
* or: ftp://ftp.iers.org/products/eop/rapid/standard/finals.data
|
|
* file description in: ftp://maia.usno.navy.mil/ser7/readme.finals
|
|
* Delta T = TAI-UT1 + 32.184 sec = (TAI-UTC) - (UT1-UTC) + 32.184 sec
|
|
*
|
|
* Also, there is the following file:
|
|
* http://maia.usno.navy.mil/ser7/deltat.data, but it is about 3 months
|
|
* behind (on 3 feb 2009); and predictions:
|
|
* http://maia.usno.navy.mil/ser7/deltat.preds
|
|
*
|
|
* Last update of table dt[]: Dieter Koch, 18 dec 2013.
|
|
* ATTENTION: Whenever updating this table, do not forget to adjust
|
|
* the macros TABEND and TABSIZ !
|
|
*/
|
|
#define TABSTART 1620
|
|
#define TABEND 2028
|
|
#define TABSIZ (TABEND-TABSTART+1)
|
|
/* we make the table greater for additional values read from external file */
|
|
#define TABSIZ_SPACE (TABSIZ+100)
|
|
static TLS double dt[TABSIZ_SPACE] = {
|
|
/* 1620.0 - 1659.0 */
|
|
124.00, 119.00, 115.00, 110.00, 106.00, 102.00, 98.00, 95.00, 91.00, 88.00,
|
|
85.00, 82.00, 79.00, 77.00, 74.00, 72.00, 70.00, 67.00, 65.00, 63.00,
|
|
62.00, 60.00, 58.00, 57.00, 55.00, 54.00, 53.00, 51.00, 50.00, 49.00,
|
|
48.00, 47.00, 46.00, 45.00, 44.00, 43.00, 42.00, 41.00, 40.00, 38.00,
|
|
/* 1660.0 - 1699.0 */
|
|
37.00, 36.00, 35.00, 34.00, 33.00, 32.00, 31.00, 30.00, 28.00, 27.00,
|
|
26.00, 25.00, 24.00, 23.00, 22.00, 21.00, 20.00, 19.00, 18.00, 17.00,
|
|
16.00, 15.00, 14.00, 14.00, 13.00, 12.00, 12.00, 11.00, 11.00, 10.00,
|
|
10.00, 10.00, 9.00, 9.00, 9.00, 9.00, 9.00, 9.00, 9.00, 9.00,
|
|
/* 1700.0 - 1739.0 */
|
|
9.00, 9.00, 9.00, 9.00, 9.00, 9.00, 9.00, 9.00, 10.00, 10.00,
|
|
10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 11.00, 11.00, 11.00,
|
|
11.00, 11.00, 11.00, 11.00, 11.00, 11.00, 11.00, 11.00, 11.00, 11.00,
|
|
11.00, 11.00, 11.00, 11.00, 12.00, 12.00, 12.00, 12.00, 12.00, 12.00,
|
|
/* 1740.0 - 1779.0 */
|
|
12.00, 12.00, 12.00, 12.00, 13.00, 13.00, 13.00, 13.00, 13.00, 13.00,
|
|
13.00, 14.00, 14.00, 14.00, 14.00, 14.00, 14.00, 14.00, 15.00, 15.00,
|
|
15.00, 15.00, 15.00, 15.00, 15.00, 16.00, 16.00, 16.00, 16.00, 16.00,
|
|
16.00, 16.00, 16.00, 16.00, 16.00, 17.00, 17.00, 17.00, 17.00, 17.00,
|
|
/* 1780.0 - 1799.0 */
|
|
17.00, 17.00, 17.00, 17.00, 17.00, 17.00, 17.00, 17.00, 17.00, 17.00,
|
|
17.00, 17.00, 16.00, 16.00, 16.00, 16.00, 15.00, 15.00, 14.00, 14.00,
|
|
/* 1800.0 - 1819.0 */
|
|
13.70, 13.40, 13.10, 12.90, 12.70, 12.60, 12.50, 12.50, 12.50, 12.50,
|
|
12.50, 12.50, 12.50, 12.50, 12.50, 12.50, 12.50, 12.40, 12.30, 12.20,
|
|
/* 1820.0 - 1859.0 */
|
|
12.00, 11.70, 11.40, 11.10, 10.60, 10.20, 9.60, 9.10, 8.60, 8.00,
|
|
7.50, 7.00, 6.60, 6.30, 6.00, 5.80, 5.70, 5.60, 5.60, 5.60,
|
|
5.70, 5.80, 5.90, 6.10, 6.20, 6.30, 6.50, 6.60, 6.80, 6.90,
|
|
7.10, 7.20, 7.30, 7.40, 7.50, 7.60, 7.70, 7.70, 7.80, 7.80,
|
|
/* 1860.0 - 1899.0 */
|
|
7.88, 7.82, 7.54, 6.97, 6.40, 6.02, 5.41, 4.10, 2.92, 1.82,
|
|
1.61, .10, -1.02, -1.28, -2.69, -3.24, -3.64, -4.54, -4.71, -5.11,
|
|
-5.40, -5.42, -5.20, -5.46, -5.46, -5.79, -5.63, -5.64, -5.80, -5.66,
|
|
-5.87, -6.01, -6.19, -6.64, -6.44, -6.47, -6.09, -5.76, -4.66, -3.74,
|
|
/* 1900.0 - 1939.0 */
|
|
-2.72, -1.54, -.02, 1.24, 2.64, 3.86, 5.37, 6.14, 7.75, 9.13,
|
|
10.46, 11.53, 13.36, 14.65, 16.01, 17.20, 18.24, 19.06, 20.25, 20.95,
|
|
21.16, 22.25, 22.41, 23.03, 23.49, 23.62, 23.86, 24.49, 24.34, 24.08,
|
|
24.02, 24.00, 23.87, 23.95, 23.86, 23.93, 23.73, 23.92, 23.96, 24.02,
|
|
/* 1940.0 - 1949.0 */
|
|
24.33, 24.83, 25.30, 25.70, 26.24, 26.77, 27.28, 27.78, 28.25, 28.71,
|
|
/* 1950.0 - 1959.0 */
|
|
29.15, 29.57, 29.97, 30.36, 30.72, 31.07, 31.35, 31.68, 32.18, 32.68,
|
|
/* 1960.0 - 1969.0 */
|
|
33.15, 33.59, 34.00, 34.47, 35.03, 35.73, 36.54, 37.43, 38.29, 39.20,
|
|
/* 1970.0 - 1979.0 */
|
|
/* from 1974 on values (with 4-digit precision) were calculated from IERS data */
|
|
40.18, 41.17, 42.23, 43.37, 44.4841, 45.4761, 46.4567, 47.5214, 48.5344, 49.5862,
|
|
/* 1980.0 - 1989.0 */
|
|
50.5387, 51.3808, 52.1668, 52.9565, 53.7882, 54.3427, 54.8713, 55.3222, 55.8197, 56.3000,
|
|
/* 1990.0 - 1999.0 */
|
|
56.8553, 57.5653, 58.3092, 59.1218, 59.9845, 60.7854, 61.6287, 62.2951, 62.9659, 63.4673,
|
|
/* 2000.0 - 2009.0 */
|
|
63.8285, 64.0908, 64.2998, 64.4734, 64.5736, 64.6876, 64.8452, 65.1464, 65.4574, 65.7768,
|
|
/* 2010.0 - 2018.0 */
|
|
66.0699, 66.3246, 66.6030, 66.9069, 67.2810, 67.6439, 68.1024, 68.5927, 68.9676, 69.2202,
|
|
/* 2020.0 - */
|
|
69.3612,
|
|
/* Extrapolated values:
|
|
* 2021 - 2028 */
|
|
69.4271, 70.00, 70.50, 71.00, 71.50, 72.00, 72.50, 73.00,
|
|
};
|
|
|
|
#define TAB2_SIZ 27
|
|
#define TAB2_START (-1000)
|
|
#define TAB2_END 1600
|
|
#define TAB2_STEP 100
|
|
#define LTERM_EQUATION_YSTART 1820
|
|
#define LTERM_EQUATION_COEFF 32
|
|
/* Table for -1000 through 1600, from Morrison & Stephenson (2004). */
|
|
static const short dt2[TAB2_SIZ] = {
|
|
/*-1000 -900 -800 -700 -600 -500 -400 -300 -200 -100*/
|
|
25400,23700,22000,21000,19040,17190,15530,14080,12790,11640,
|
|
/* 0 100 200 300 400 500 600 700 800 900*/
|
|
10580, 9600, 8640, 7680, 6700, 5710, 4740, 3810, 2960, 2200,
|
|
/* 1000 1100 1200 1300 1400 1500 1600, */
|
|
1570, 1090, 740, 490, 320, 200, 120,
|
|
};
|
|
|
|
/* Table for -500 through 1600, from Stephenson & Morrison (1995).
|
|
*
|
|
* The first value for -550 has been added from Borkowski
|
|
* in order to make this table fit with the Borkowski formula
|
|
* for times before -550.
|
|
*/
|
|
#define TAB97_SIZ 43
|
|
#define TAB97_START (-500)
|
|
#define TAB97_END (1600)
|
|
#define TAB97_STEP (50)
|
|
static const short dt97[TAB97_SIZ] = {
|
|
/* -500 -450 -400 -350 -300 -250 -200 -150 -100 -50*/
|
|
16800,16000,15300,14600,14000,13400,12800,12200,11600,11100,
|
|
/* 0 50 100 150 200 250 300 350 400 450*/
|
|
10600,10100, 9600, 9100, 8600, 8200, 7700, 7200, 6700, 6200,
|
|
/* 500 550 600 650 700 750 800 850 900 950*/
|
|
5700, 5200, 4700, 4300, 3800, 3400, 3000, 2600, 2200, 1900,
|
|
/* 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450*/
|
|
1600, 1350, 1100, 900, 750, 600, 470, 380, 300, 230,
|
|
/* 1500 1550 1600 */
|
|
180, 140, 110,
|
|
};
|
|
|
|
/* returns DeltaT (ET - UT) in days
|
|
* double tjd = julian day in UT
|
|
* delta t is adjusted to the tidal acceleration that is compatible
|
|
* with the ephemeris flag contained in iflag and with the ephemeris
|
|
* files made accessible through swe_set_ephe_path() or swe_set_jplfile().
|
|
* If iflag = -1, then the default tidal acceleration is ussed (i.e.
|
|
* that of DE431).
|
|
*/
|
|
#define DEMO 0
|
|
static int32 calc_deltat(double tjd, int32 iflag, double *deltat, char *serr)
|
|
{
|
|
double ans = 0;
|
|
double B, Y, Ygreg, dd;
|
|
int iy;
|
|
int32 retc;
|
|
int deltat_model = swed.astro_models[SE_MODEL_DELTAT];
|
|
double tid_acc;
|
|
int32 denum, denumret;
|
|
int32 epheflag, otherflag;
|
|
//fprintf(stderr, "dmod=%f, %.f\n", (double) deltat_model, (double) SEMOD_DELTAT_DEFAULT);
|
|
if (deltat_model == 0) deltat_model = SEMOD_DELTAT_DEFAULT;
|
|
epheflag = iflag & SEFLG_EPHMASK;
|
|
otherflag = iflag & ~SEFLG_EPHMASK;
|
|
/* with iflag == -1, we use default tid_acc */
|
|
if (iflag == -1) {
|
|
retc = swi_get_tid_acc(tjd, 0, 9999, &denumret, &tid_acc, serr); /* for default tid_acc */
|
|
/* otherwise we use tid_acc consistent with epheflag */
|
|
} else {
|
|
denum = swed.jpldenum;
|
|
if (epheflag & SEFLG_SWIEPH) denum = swed.fidat[SEI_FILE_MOON].sweph_denum;
|
|
if (swi_init_swed_if_start() == 1 && !(epheflag & SEFLG_MOSEPH)) {
|
|
if (serr != NULL)
|
|
strcpy(serr, "Please call swe_set_ephe_path() or swe_set_jplfile() before calling swe_deltat_ex()");
|
|
retc = swi_set_tid_acc(tjd, epheflag, denum, NULL); /* _set_ saves tid_acc in swed */
|
|
} else {
|
|
retc = swi_set_tid_acc(tjd, epheflag, denum, serr); /* _set_ saves tid_acc in swed */
|
|
}
|
|
tid_acc = swed.tid_acc;
|
|
}
|
|
iflag = otherflag | retc;
|
|
Y = 2000.0 + (tjd - J2000)/365.25;
|
|
Ygreg = 2000.0 + (tjd - J2000)/365.2425;
|
|
/* Model for epochs before 1955, currently default in Swiss Ephemeris:
|
|
* Stephenson/Morrison/Hohenkerk 2016
|
|
* (we switch over to Astronomical Almanac K8-K9 and IERS at 1 Jan. 1955.
|
|
* To make the curve continuous we apply some linear term over
|
|
* 1000 days before this date.)
|
|
* Delta T according to Stephenson/Morrison/Hohenkerk 2016 is based on
|
|
* ancient, medieval, and modern observations of eclipses and occultations.
|
|
* Values of Deltat T before 1955 depend on this kind of observations.
|
|
* For more recent data we want to use the data provided by IERS
|
|
* (or Astronomical Almanac K8-K9).
|
|
*/
|
|
if (deltat_model == SEMOD_DELTAT_STEPHENSON_ETC_2016 && tjd < 2435108.5) { // tjd < 2432521.453645833) {
|
|
*deltat = deltat_stephenson_etc_2016(tjd, tid_acc);
|
|
if (tjd >= 2434108.5) {
|
|
*deltat += (1.0 - (2435108.5 - tjd) / 1000.0) * 0.6610218 / 86400.0;
|
|
}
|
|
return iflag;
|
|
}
|
|
/* Model used SE 1.77 - 2.05.01, for epochs before 1633:
|
|
* Polynomials by Espenak & Meeus 2006,
|
|
* derived from Stephenson & Morrison 2004.
|
|
* deltat_model == SEMOD_DELTAT_ESPENAK_MEEUS_2006:
|
|
* This method is used only for epochs before 1633. (For more recent
|
|
* epochs, we use the data provided by Astronomical Almanac K8-K9.)
|
|
*/
|
|
if (deltat_model == SEMOD_DELTAT_ESPENAK_MEEUS_2006 && tjd < 2317746.13090277789) {
|
|
*deltat = deltat_espenak_meeus_1620(tjd, tid_acc);
|
|
return iflag;
|
|
}
|
|
/* delta t model used in SE 1.72 - 1.76:
|
|
* Stephenson & Morrison 2004;
|
|
* before 1620 */
|
|
if (deltat_model == SEMOD_DELTAT_STEPHENSON_MORRISON_2004 && Y < TABSTART) {
|
|
// before 1600:
|
|
if (Y < TAB2_END) {
|
|
*deltat = deltat_stephenson_morrison_2004_1600(tjd, tid_acc);
|
|
return iflag;
|
|
} else {
|
|
/* between 1600 and 1620:
|
|
* linear interpolation between
|
|
* end of table dt2 and start of table dt */
|
|
if (Y >= TAB2_END) {
|
|
B = TABSTART - TAB2_END;
|
|
iy = (TAB2_END - TAB2_START) / TAB2_STEP;
|
|
dd = (Y - TAB2_END) / B;
|
|
ans = dt2[iy] + dd * (dt[0] - dt2[iy]);
|
|
ans = adjust_for_tidacc(ans, Ygreg, tid_acc, SE_TIDAL_26, FALSE);
|
|
*deltat = ans / 86400.0;
|
|
return iflag;
|
|
}
|
|
}
|
|
}
|
|
/* delta t model used in SE 1.64 - 1.71:
|
|
* Stephenson 1997;
|
|
* before 1620 */
|
|
if (deltat_model == SEMOD_DELTAT_STEPHENSON_1997 && Y < TABSTART) {
|
|
// before 1600:
|
|
if (Y < TAB97_END) {
|
|
*deltat = deltat_stephenson_morrison_1997_1600(tjd, tid_acc);
|
|
return iflag;
|
|
} else {
|
|
/* between 1600 and 1620:
|
|
* linear interpolation between
|
|
* end of table dt2 and start of table dt */
|
|
if (Y >= TAB97_END) {
|
|
B = TABSTART - TAB97_END;
|
|
iy = (TAB97_END - TAB97_START) / TAB97_STEP;
|
|
dd = (Y - TAB97_END) / B;
|
|
ans = dt97[iy] + dd * (dt[0] - dt97[iy]);
|
|
ans = adjust_for_tidacc(ans, Ygreg, tid_acc, SE_TIDAL_26, FALSE);
|
|
*deltat = ans / 86400.0;
|
|
return iflag;
|
|
}
|
|
}
|
|
}
|
|
/* delta t model used before SE 1.64:
|
|
* Stephenson/Morrison 1984 with Borkowski 1988;
|
|
* before 1620 */
|
|
if (deltat_model == SEMOD_DELTAT_STEPHENSON_MORRISON_1984 && Y < TABSTART) {
|
|
if( Y >= 948.0 ) {
|
|
/* Stephenson and Morrison, stated domain is 948 to 1600:
|
|
* 25.5(centuries from 1800)^2 - 1.9159(centuries from 1955)^2 */
|
|
B = 0.01 * (Y - 2000.0);
|
|
ans = (23.58 * B + 100.3) * B + 101.6;
|
|
} else {
|
|
/* Borkowski, before 948 and between 1600 and 1620 */
|
|
B = 0.01 * (Y - 2000.0) + 3.75;
|
|
ans = 35.0 * B * B + 40.;
|
|
}
|
|
*deltat = ans / 86400.0;
|
|
return iflag;
|
|
}
|
|
/* 1620 - today + a few years (tabend):
|
|
* Tabulated values of deltaT from Astronomical Almanac
|
|
* (AA 1997etc., pp. K8-K9) and from IERS
|
|
* (http://maia.usno.navy.mil/ser7/deltat.data).
|
|
*/
|
|
if (Y >= TABSTART) {
|
|
*deltat = deltat_aa(tjd, tid_acc);
|
|
return iflag;
|
|
}
|
|
#ifdef TRACE
|
|
swi_open_trace(NULL);
|
|
if (swi_trace_count < TRACE_COUNT_MAX) {
|
|
if (swi_fp_trace_c != NULL) {
|
|
fputs("\n/*SWE_DELTAT*/\n", swi_fp_trace_c);
|
|
fprintf(swi_fp_trace_c, " tjd = %.9f;", tjd);
|
|
fprintf(swi_fp_trace_c, " iflag = %d;", iflag);
|
|
fprintf(swi_fp_trace_c, " t = swe_deltat_ex(tjd, iflag, NULL);\n");
|
|
fputs(" printf(\"swe_deltat: %f\\t%f\\t\\n\", ", swi_fp_trace_c);
|
|
fputs("tjd, t);\n", swi_fp_trace_c);
|
|
fflush(swi_fp_trace_c);
|
|
}
|
|
if (swi_fp_trace_out != NULL) {
|
|
fprintf(swi_fp_trace_out, "swe_deltat: %f\t%f\t\n", tjd, ans);
|
|
fflush(swi_fp_trace_out);
|
|
}
|
|
}
|
|
#endif
|
|
*deltat = ans / 86400.0;
|
|
return iflag;
|
|
}
|
|
|
|
double CALL_CONV swe_deltat_ex(double tjd, int32 iflag, char *serr)
|
|
{
|
|
double deltat;
|
|
if (swed.delta_t_userdef_is_set)
|
|
return swed.delta_t_userdef;
|
|
if (serr != NULL)
|
|
*serr = '\0';
|
|
calc_deltat(tjd, iflag, &deltat, serr);
|
|
return deltat;
|
|
}
|
|
|
|
double CALL_CONV swe_deltat(double tjd)
|
|
{
|
|
int32 iflag = swi_guess_ephe_flag();
|
|
return swe_deltat_ex(tjd, iflag, NULL); /* with default tidal acceleration/default ephemeris */
|
|
}
|
|
|
|
/* The tabulated values of deltaT, in hundredths of a second,
|
|
* were taken from The Astronomical Almanac 1997etc., pp. K8-K9.
|
|
* Some more recent values are taken from IERS
|
|
* http://maia.usno.navy.mil/ser7/deltat.data .
|
|
* Bessel's interpolation formula is implemented to obtain fourth
|
|
* order interpolated values at intermediate times.
|
|
* The values are adjusted depending on the ephemeris used
|
|
* and its inherent value of secular tidal acceleration ndot.
|
|
* Note by Dieter Jan. 2017:
|
|
* Bessel interpolation assumes equidistant sampling points. However the
|
|
* sampling points are not equidistant, because they are for first January of
|
|
* every year and years can have either 365 or 366 days. The interpolation uses
|
|
* a step width of 365.25 days. As a consequence, in three out of four years
|
|
* the interpolation does not reproduce the exact values of the sampling points
|
|
* on the days they refer to. */
|
|
static double deltat_aa(double tjd, double tid_acc)
|
|
{
|
|
double ans = 0, ans2 = 0, ans3;
|
|
double p, B, B2, Y, dd;
|
|
double d[6];
|
|
int i, iy, k;
|
|
/* read additional values from swedelta.txt */
|
|
int tabsiz = init_dt();
|
|
int tabend = TABSTART + tabsiz - 1;
|
|
int deltat_model = swed.astro_models[SE_MODEL_DELTAT];
|
|
if (deltat_model == 0) deltat_model = SEMOD_DELTAT_DEFAULT;
|
|
Y = 2000.0 + (tjd - 2451544.5)/365.25;
|
|
if (Y <= tabend) {
|
|
/* Index into the table.
|
|
*/
|
|
p = floor(Y);
|
|
iy = (int) (p - TABSTART);
|
|
/* Zeroth order estimate is value at start of year */
|
|
ans = dt[iy];
|
|
k = iy + 1;
|
|
if( k >= tabsiz )
|
|
goto done; /* No data, can't go on. */
|
|
/* The fraction of tabulation interval */
|
|
p = Y - p;
|
|
/* First order interpolated value */
|
|
ans += p*(dt[k] - dt[iy]);
|
|
if( (iy-1 < 0) || (iy+2 >= tabsiz) )
|
|
goto done; /* can't do second differences */
|
|
/* Make table of first differences */
|
|
k = iy - 2;
|
|
for( i=0; i<5; i++ ) {
|
|
if( (k < 0) || (k+1 >= tabsiz) )
|
|
d[i] = 0;
|
|
else
|
|
d[i] = dt[k+1] - dt[k];
|
|
k += 1;
|
|
}
|
|
/* Compute second differences */
|
|
for( i=0; i<4; i++ )
|
|
d[i] = d[i+1] - d[i];
|
|
B = 0.25*p*(p-1.0);
|
|
ans += B*(d[1] + d[2]);
|
|
#if DEMO
|
|
printf( "B %.4lf, ans %.4lf\n", B, ans );
|
|
#endif
|
|
if( iy+2 >= tabsiz )
|
|
goto done;
|
|
/* Compute third differences */
|
|
for( i=0; i<3; i++ )
|
|
d[i] = d[i+1] - d[i];
|
|
B = 2.0*B/3.0;
|
|
ans += (p-0.5)*B*d[1];
|
|
#if DEMO
|
|
printf( "B %.4lf, ans %.4lf\n", B*(p-0.5), ans );
|
|
#endif
|
|
if( (iy-2 < 0) || (iy+3 > tabsiz) )
|
|
goto done;
|
|
/* Compute fourth differences */
|
|
for( i=0; i<2; i++ )
|
|
d[i] = d[i+1] - d[i];
|
|
B = 0.125*B*(p+1.0)*(p-2.0);
|
|
ans += B*(d[0] + d[1]);
|
|
#if DEMO
|
|
printf( "B %.4lf, ans %.4lf\n", B, ans );
|
|
#endif
|
|
done:
|
|
ans = adjust_for_tidacc(ans, Y, tid_acc, SE_TIDAL_26, FALSE);
|
|
return ans / 86400.0;
|
|
}
|
|
/* today - future:
|
|
* 3rd degree polynomial based on data given by
|
|
* Stephenson/Morrison/Hohenkerk 2016 here:
|
|
* http://astro.ukho.gov.uk/nao/lvm/
|
|
*/
|
|
if (deltat_model == SEMOD_DELTAT_STEPHENSON_ETC_2016) {
|
|
B = (Y - 2000);
|
|
if (Y < 2500) {
|
|
ans = B * B * B * 121.0 / 30000000.0 + B * B / 1250.0 + B * 521.0 / 3000.0 + 64.0;
|
|
/* for slow transition from tablulated data */
|
|
B2 = (tabend - 2000);
|
|
ans2 = B2 * B2 * B2 * 121.0 / 30000000.0 + B2 * B2 / 1250.0 + B2 * 521.0 / 3000.0 + 64.0;
|
|
/* we use a parable after 2500 */
|
|
} else {
|
|
B = 0.01 * (Y - 2000);
|
|
ans = B * B * 32.5 + 42.5;
|
|
}
|
|
/*
|
|
* Formula Stephenson (1997; p. 507),
|
|
* with modification to avoid jump at end of AA table,
|
|
* similar to what Meeus 1998 had suggested.
|
|
* Slow transition within 100 years.
|
|
*/
|
|
} else {
|
|
B = 0.01 * (Y - 1820);
|
|
ans = -20 + 31 * B * B;
|
|
/* for slow transition from tablulated data */
|
|
B2 = 0.01 * (tabend - 1820);
|
|
ans2 = -20 + 31 * B2 * B2;
|
|
}
|
|
/* slow transition from tabulated values to Stephenson formula: */
|
|
if (Y <= tabend+100) {
|
|
ans3 = dt[tabsiz-1];
|
|
dd = (ans2 - ans3);
|
|
ans += dd * (Y - (tabend + 100)) * 0.01;
|
|
}
|
|
return ans / 86400.0;
|
|
}
|
|
|
|
static double deltat_longterm_morrison_stephenson(double tjd)
|
|
{
|
|
double Ygreg = 2000.0 + (tjd - J2000)/365.2425;
|
|
double u = (Ygreg - 1820) / 100.0;
|
|
return (-20 + 32 * u * u);
|
|
}
|
|
|
|
static double deltat_stephenson_morrison_1997_1600(double tjd, double tid_acc)
|
|
{
|
|
double ans = 0, ans2, ans3;
|
|
double p, B, Y, dd;
|
|
int iy;
|
|
Y = 2000.0 + (tjd - J2000)/365.25;
|
|
/* before -500:
|
|
* formula by Stephenson (1997; p. 508) but adjusted to fit the starting
|
|
* point of table dt97 (Stephenson 1997). */
|
|
if( Y < TAB97_START ) {
|
|
B = (Y - 1735) * 0.01;
|
|
ans = -20 + 35 * B * B;
|
|
ans = adjust_for_tidacc(ans, Y, tid_acc, SE_TIDAL_26, FALSE);
|
|
/* transition from formula to table over 100 years */
|
|
if (Y >= TAB97_START - 100) {
|
|
/* starting value of table dt97: */
|
|
ans2 = adjust_for_tidacc(dt97[0], TAB97_START, tid_acc, SE_TIDAL_26, FALSE);
|
|
/* value of formula at epoch TAB97_START */
|
|
B = (TAB97_START - 1735) * 0.01;
|
|
ans3 = -20 + 35 * B * B;
|
|
ans3 = adjust_for_tidacc(ans3, Y, tid_acc, SE_TIDAL_26, FALSE);
|
|
dd = ans3 - ans2;
|
|
B = (Y - (TAB97_START - 100)) * 0.01;
|
|
/* fit to starting point of table dt97. */
|
|
ans = ans - dd * B;
|
|
}
|
|
}
|
|
/* between -500 and 1600:
|
|
* linear interpolation between values of table dt97 (Stephenson 1997) */
|
|
if (Y >= TAB97_START && Y < TAB2_END) {
|
|
p = floor(Y);
|
|
iy = (int) ((p - TAB97_START) / 50.0);
|
|
dd = (Y - (TAB97_START + 50 * iy)) / 50.0;
|
|
ans = dt97[iy] + (dt97[iy+1] - dt97[iy]) * dd;
|
|
/* correction for tidal acceleration used by our ephemeris */
|
|
ans = adjust_for_tidacc(ans, Y, tid_acc, SE_TIDAL_26, FALSE);
|
|
}
|
|
ans /= 86400.0;
|
|
return ans;
|
|
}
|
|
|
|
/* Stephenson & Morrison (2004) */
|
|
static double deltat_stephenson_morrison_2004_1600(double tjd, double tid_acc)
|
|
{
|
|
double ans = 0, ans2, ans3;
|
|
double p, B, dd;
|
|
double tjd0;
|
|
int iy;
|
|
double Y = 2000.0 + (tjd - J2000)/365.2425;
|
|
/* double Y = 2000.0 + (tjd - J2000)/365.25;*/
|
|
/* before -1000:
|
|
* formula by Stephenson & Morrison (2004; p. 335) but adjusted to fit the
|
|
* starting point of table dt2. */
|
|
if( Y < TAB2_START ) { // before -1000
|
|
ans = deltat_longterm_morrison_stephenson(tjd);
|
|
ans = adjust_for_tidacc(ans, Y, tid_acc, SE_TIDAL_26, FALSE);
|
|
/* transition from formula to table over 100 years */
|
|
if (Y >= TAB2_START - 100) {
|
|
/* starting value of table dt2: */
|
|
ans2 = adjust_for_tidacc(dt2[0], TAB2_START, tid_acc, SE_TIDAL_26, FALSE);
|
|
/* value of formula at epoch TAB2_START */
|
|
/* B = (TAB2_START - LTERM_EQUATION_YSTART) * 0.01;
|
|
ans3 = -20 + LTERM_EQUATION_COEFF * B * B;*/
|
|
tjd0 = (TAB2_START - 2000) * 365.2425 + J2000;
|
|
ans3 = deltat_longterm_morrison_stephenson(tjd0);
|
|
ans3 = adjust_for_tidacc(ans3, Y, tid_acc, SE_TIDAL_26, FALSE);
|
|
dd = ans3 - ans2;
|
|
B = (Y - (TAB2_START - 100)) * 0.01;
|
|
/* fit to starting point of table dt2. */
|
|
ans = ans - dd * B;
|
|
}
|
|
}
|
|
/* between -1000 and 1600:
|
|
* linear interpolation between values of table dt2 (Stephenson & Morrison 2004) */
|
|
if (Y >= TAB2_START && Y < TAB2_END) {
|
|
double Yjul = 2000 + (tjd - 2451557.5) / 365.25;
|
|
p = floor(Yjul);
|
|
iy = (int) ((p - TAB2_START) / TAB2_STEP);
|
|
dd = (Yjul - (TAB2_START + TAB2_STEP * iy)) / TAB2_STEP;
|
|
ans = dt2[iy] + (dt2[iy+1] - dt2[iy]) * dd;
|
|
/* correction for tidal acceleration used by our ephemeris */
|
|
ans = adjust_for_tidacc(ans, Y, tid_acc, SE_TIDAL_26, FALSE);
|
|
}
|
|
ans /= 86400.0;
|
|
return ans;
|
|
}
|
|
|
|
/*
|
|
* These coefficients represent the spline approximation discussed in the
|
|
* paper "Measurement of the Earth's Rotation: 720 BC to AD 2015",
|
|
* Stephenson, F.R., Morrison, L.V., and Hohenkerk, C.Y., published by
|
|
* Royal Society Proceedings A and available from their website at
|
|
* http://rspa.royalsocietypublishing.org/lookup/doi/10.1098/rspa.2016.0404.
|
|
* Year numbers have been replaced by Julian day numbers by D. Koch.
|
|
*/
|
|
#define NDTCF16 54
|
|
double dtcf16[NDTCF16][6] =
|
|
{
|
|
/*00*/ {1458085.5, 1867156.5, 20550.593,-21268.478, 11863.418, -4541.129}, /* ybeg=-720, yend= 400 */
|
|
/*01*/ {1867156.5, 2086302.5, 6604.404, -5981.266, -505.093, 1349.609}, /* ybeg= 400, yend=1000 */
|
|
/*02*/ {2086302.5, 2268923.5, 1467.654, -2452.187, 2460.927, -1183.759}, /* ybeg=1000, yend=1500 */
|
|
/*03*/ {2268923.5, 2305447.5, 292.635, -216.322, -43.614, 56.681}, /* ybeg=1500, yend=1600 */
|
|
/*04*/ {2305447.5, 2323710.5, 89.380, -66.754, 31.607, -10.497}, /* ybeg=1600, yend=1650 */
|
|
/*05*/ {2323710.5, 2349276.5, 43.736, -49.043, 0.227, 15.811}, /* ybeg=1650, yend=1720 */
|
|
/*06*/ {2349276.5, 2378496.5, 10.730, -1.321, 62.250, -52.946}, /* ybeg=1720, yend=1800 */
|
|
/*07*/ {2378496.5, 2382148.5, 18.714, -4.457, -1.509, 2.507}, /* ybeg=1800, yend=1810 */
|
|
/*08*/ {2382148.5, 2385800.5, 15.255, 0.046, 6.012, -4.634}, /* ybeg=1810, yend=1820 */
|
|
/*09*/ {2385800.5, 2389453.5, 16.679, -1.831, -7.889, 3.799}, /* ybeg=1820, yend=1830 */
|
|
/*10*/ {2389453.5, 2393105.5, 10.758, -6.211, 3.509, -0.388}, /* ybeg=1830, yend=1840 */
|
|
/*11*/ {2393105.5, 2396758.5, 7.668, -0.357, 2.345, -0.338}, /* ybeg=1840, yend=1850 */
|
|
/*12*/ {2396758.5, 2398584.5, 9.317, 1.659, 0.332, -0.932}, /* ybeg=1850, yend=1855 */
|
|
/*13*/ {2398584.5, 2400410.5, 10.376, -0.472, -2.463, 1.596}, /* ybeg=1855, yend=1860 */
|
|
/*14*/ {2400410.5, 2402237.5, 9.038, -0.610, 2.325, -2.497}, /* ybeg=1860, yend=1865 */
|
|
/*15*/ {2402237.5, 2404063.5, 8.256, -3.450, -5.166, 2.729}, /* ybeg=1865, yend=1870 */
|
|
/*16*/ {2404063.5, 2405889.5, 2.369, -5.596, 3.020, -0.919}, /* ybeg=1870, yend=1875 */
|
|
/*17*/ {2405889.5, 2407715.5, -1.126, -2.312, 0.264, -0.037}, /* ybeg=1875, yend=1880 */
|
|
/*18*/ {2407715.5, 2409542.5, -3.211, -1.894, 0.154, 0.562}, /* ybeg=1880, yend=1885 */
|
|
/*19*/ {2409542.5, 2411368.5, -4.388, 0.101, 1.841, -1.438}, /* ybeg=1885, yend=1890 */
|
|
/*20*/ {2411368.5, 2413194.5, -3.884, -0.531, -2.473, 1.870}, /* ybeg=1890, yend=1895 */
|
|
/*21*/ {2413194.5, 2415020.5, -5.017, 0.134, 3.138, -0.232}, /* ybeg=1895, yend=1900 */
|
|
/*22*/ {2415020.5, 2416846.5, -1.977, 5.715, 2.443, -1.257}, /* ybeg=1900, yend=1905 */
|
|
/*23*/ {2416846.5, 2418672.5, 4.923, 6.828, -1.329, 0.720}, /* ybeg=1905, yend=1910 */
|
|
/*24*/ {2418672.5, 2420498.5, 11.142, 6.330, 0.831, -0.825}, /* ybeg=1910, yend=1915 */
|
|
/*25*/ {2420498.5, 2422324.5, 17.479, 5.518, -1.643, 0.262}, /* ybeg=1915, yend=1920 */
|
|
/*26*/ {2422324.5, 2424151.5, 21.617, 3.020, -0.856, 0.008}, /* ybeg=1920, yend=1925 */
|
|
/*27*/ {2424151.5, 2425977.5, 23.789, 1.333, -0.831, 0.127}, /* ybeg=1925, yend=1930 */
|
|
/*28*/ {2425977.5, 2427803.5, 24.418, 0.052, -0.449, 0.142}, /* ybeg=1930, yend=1935 */
|
|
/*29*/ {2427803.5, 2429629.5, 24.164, -0.419, -0.022, 0.702}, /* ybeg=1935, yend=1940 */
|
|
/*30*/ {2429629.5, 2431456.5, 24.426, 1.645, 2.086, -1.106}, /* ybeg=1940, yend=1945 */
|
|
/*31*/ {2431456.5, 2433282.5, 27.050, 2.499, -1.232, 0.614}, /* ybeg=1945, yend=1950 */
|
|
/*32*/ {2433282.5, 2434378.5, 28.932, 1.127, 0.220, -0.277}, /* ybeg=1950, yend=1953 */
|
|
/*33*/ {2434378.5, 2435473.5, 30.002, 0.737, -0.610, 0.631}, /* ybeg=1953, yend=1956 */
|
|
/*34*/ {2435473.5, 2436569.5, 30.760, 1.409, 1.282, -0.799}, /* ybeg=1956, yend=1959 */
|
|
/*35*/ {2436569.5, 2437665.5, 32.652, 1.577, -1.115, 0.507}, /* ybeg=1959, yend=1962 */
|
|
/*36*/ {2437665.5, 2438761.5, 33.621, 0.868, 0.406, 0.199}, /* ybeg=1962, yend=1965 */
|
|
/*37*/ {2438761.5, 2439856.5, 35.093, 2.275, 1.002, -0.414}, /* ybeg=1965, yend=1968 */
|
|
/*38*/ {2439856.5, 2440952.5, 37.956, 3.035, -0.242, 0.202}, /* ybeg=1968, yend=1971 */
|
|
/*39*/ {2440952.5, 2442048.5, 40.951, 3.157, 0.364, -0.229}, /* ybeg=1971, yend=1974 */
|
|
/*40*/ {2442048.5, 2443144.5, 44.244, 3.198, -0.323, 0.172}, /* ybeg=1974, yend=1977 */
|
|
/*41*/ {2443144.5, 2444239.5, 47.291, 3.069, 0.193, -0.192}, /* ybeg=1977, yend=1980 */
|
|
/*42*/ {2444239.5, 2445335.5, 50.361, 2.878, -0.384, 0.081}, /* ybeg=1980, yend=1983 */
|
|
/*43*/ {2445335.5, 2446431.5, 52.936, 2.354, -0.140, -0.166}, /* ybeg=1983, yend=1986 */
|
|
/*44*/ {2446431.5, 2447527.5, 54.984, 1.577, -0.637, 0.448}, /* ybeg=1986, yend=1989 */
|
|
/*45*/ {2447527.5, 2448622.5, 56.373, 1.649, 0.709, -0.277}, /* ybeg=1989, yend=1992 */
|
|
/*46*/ {2448622.5, 2449718.5, 58.453, 2.235, -0.122, 0.111}, /* ybeg=1992, yend=1995 */
|
|
/*47*/ {2449718.5, 2450814.5, 60.677, 2.324, 0.212, -0.315}, /* ybeg=1995, yend=1998 */
|
|
/*48*/ {2450814.5, 2451910.5, 62.899, 1.804, -0.732, 0.112}, /* ybeg=1998, yend=2001 */
|
|
/*49*/ {2451910.5, 2453005.5, 64.082, 0.675, -0.396, 0.193}, /* ybeg=2001, yend=2004 */
|
|
/*50*/ {2453005.5, 2454101.5, 64.555, 0.463, 0.184, -0.008}, /* ybeg=2004, yend=2007 */
|
|
/*51*/ {2454101.5, 2455197.5, 65.194, 0.809, 0.161, -0.101}, /* ybeg=2007, yend=2010 */
|
|
/*52*/ {2455197.5, 2456293.5, 66.063, 0.828, -0.142, 0.168}, /* ybeg=2010, yend=2013 */
|
|
/*53*/ {2456293.5, 2457388.5, 66.917, 1.046, 0.360, -0.282}, /* ybeg=2013, yend=2016 */
|
|
};
|
|
static double deltat_stephenson_etc_2016(double tjd, double tid_acc)
|
|
{
|
|
double t, dt, Ygreg;
|
|
int i, irec = -1;
|
|
Ygreg = 2000.0 + (tjd - J2000)/365.2425;
|
|
// after the year -720 get value from spline curve
|
|
for (i = 0; i < NDTCF16; i++) {
|
|
if (tjd < dtcf16[i][0]) break;
|
|
if (tjd < dtcf16[i][1]) {
|
|
irec = i;
|
|
break;
|
|
}
|
|
}
|
|
if (irec >= 0) {
|
|
t = (tjd - dtcf16[irec][0]) / (dtcf16[irec][1] - dtcf16[irec][0]);
|
|
dt = dtcf16[irec][2] + dtcf16[irec][3] * t + dtcf16[irec][4] * t * t + dtcf16[irec][5] * t * t * t;
|
|
// for earlier epochs, use long term parabola
|
|
} else if (Ygreg < -720) {
|
|
t = (Ygreg - 1825) / 100.0;
|
|
dt = -320 + 32.5 * t * t;
|
|
dt -= 179.7337208; // to make curve continous on 1 Jan -720 (D. Koch)
|
|
// future
|
|
} else {
|
|
t = (Ygreg - 1825) / 100.0;
|
|
dt = -320 + 32.5 * t * t;
|
|
dt += 269.4790417; // to make curve continous on 1 Jan 2016 (D. Koch)
|
|
}
|
|
/* The parameter adjust_after_1955 must be TRUE here, because the
|
|
* Stephenson 2016 curve is based on occultation data alone,
|
|
* not on IERS data.
|
|
* Note, however, the current function deltat_stephenson_etc_2016()
|
|
* is called only for dates before 1 Jan 1955. */
|
|
dt = adjust_for_tidacc(dt, Ygreg, tid_acc, SE_TIDAL_STEPHENSON_2016, TRUE);
|
|
dt /= 86400.0;
|
|
return dt;
|
|
}
|
|
|
|
static double deltat_espenak_meeus_1620(double tjd, double tid_acc)
|
|
{
|
|
double ans = 0;
|
|
double Ygreg;
|
|
double u;
|
|
/* double Y = 2000.0 + (tjd - J2000)/365.25;*/
|
|
Ygreg = 2000.0 + (tjd - J2000)/365.2425;
|
|
if (Ygreg < -500) {
|
|
ans = deltat_longterm_morrison_stephenson(tjd);
|
|
} else if (Ygreg < 500) {
|
|
u = Ygreg / 100.0;
|
|
ans = (((((0.0090316521 * u + 0.022174192) * u - 0.1798452) * u - 5.952053) * u+ 33.78311) * u - 1014.41) * u + 10583.6;
|
|
} else if (Ygreg < 1600) {
|
|
u = (Ygreg - 1000) / 100.0;
|
|
ans = (((((0.0083572073 * u - 0.005050998) * u - 0.8503463) * u + 0.319781) * u + 71.23472) * u - 556.01) * u + 1574.2;
|
|
} else if (Ygreg < 1700) {
|
|
u = Ygreg - 1600;
|
|
ans = 120 - 0.9808 * u - 0.01532 * u * u + u * u * u / 7129.0;
|
|
} else if (Ygreg < 1800) {
|
|
u = Ygreg - 1700;
|
|
ans = (((-u / 1174000.0 + 0.00013336) * u - 0.0059285) * u + 0.1603) * u + 8.83;
|
|
} else if (Ygreg < 1860) {
|
|
u = Ygreg - 1800;
|
|
ans = ((((((0.000000000875 * u - 0.0000001699) * u + 0.0000121272) * u - 0.00037436) * u + 0.0041116) * u + 0.0068612) * u - 0.332447) * u + 13.72;
|
|
} else if (Ygreg < 1900) {
|
|
u = Ygreg - 1860;
|
|
ans = ((((u / 233174.0 - 0.0004473624) * u + 0.01680668) * u - 0.251754) * u + 0.5737) * u + 7.62;
|
|
} else if (Ygreg < 1920) {
|
|
u = Ygreg - 1900;
|
|
ans = (((-0.000197 * u + 0.0061966) * u - 0.0598939) * u + 1.494119) * u -2.79;
|
|
} else if (Ygreg < 1941) {
|
|
u = Ygreg - 1920;
|
|
ans = 21.20 + 0.84493 * u - 0.076100 * u * u + 0.0020936 * u * u * u;
|
|
} else if (Ygreg < 1961) {
|
|
u = Ygreg - 1950;
|
|
ans = 29.07 + 0.407 * u - u * u / 233.0 + u * u * u / 2547.0;
|
|
} else if (Ygreg < 1986) {
|
|
u = Ygreg - 1975;
|
|
ans = 45.45 + 1.067 * u - u * u / 260.0 - u * u * u / 718.0;
|
|
} else if (Ygreg < 2005) {
|
|
u = Ygreg - 2000;
|
|
ans = ((((0.00002373599 * u + 0.000651814) * u + 0.0017275) * u - 0.060374) * u + 0.3345) * u + 63.86;
|
|
}
|
|
ans = adjust_for_tidacc(ans, Ygreg, tid_acc, SE_TIDAL_26, FALSE);
|
|
ans /= 86400.0;
|
|
return ans;
|
|
}
|
|
|
|
/* Read delta t values from external file.
|
|
* record structure: year(whitespace)delta_t in 0.01 sec.
|
|
*/
|
|
static int init_dt(void)
|
|
{
|
|
FILE *fp;
|
|
int year;
|
|
int tab_index;
|
|
int tabsiz;
|
|
int i;
|
|
char s[AS_MAXCH];
|
|
char *sp;
|
|
if (!swed.init_dt_done) {
|
|
swed.init_dt_done = TRUE;
|
|
/* no error message if file is missing */
|
|
if ((fp = swi_fopen(-1, "swe_deltat.txt", swed.ephepath, NULL)) == NULL
|
|
&& (fp = swi_fopen(-1, "sedeltat.txt", swed.ephepath, NULL)) == NULL)
|
|
return TABSIZ;
|
|
while(fgets(s, AS_MAXCH, fp) != NULL) {
|
|
sp = s;
|
|
while (strchr(" \t", *sp) != NULL && *sp != '\0')
|
|
sp++; /* was *sp++ fixed by Alois 2-jul-2003 */
|
|
if (*sp == '#' || *sp == '\n')
|
|
continue;
|
|
year = atoi(s);
|
|
tab_index = year - TABSTART;
|
|
/* table space is limited. no error msg, if exceeded */
|
|
if (tab_index >= TABSIZ_SPACE)
|
|
continue;
|
|
sp += 4;
|
|
while (strchr(" \t", *sp) != NULL && *sp != '\0')
|
|
sp++; /* was *sp++ fixed by Alois 2-jul-2003 */
|
|
/*dt[tab_index] = (short) (atof(sp) * 100 + 0.5);*/
|
|
dt[tab_index] = atof(sp);
|
|
}
|
|
fclose(fp);
|
|
}
|
|
/* find table size */
|
|
tabsiz = 2001 - TABSTART + 1;
|
|
for (i = tabsiz - 1; i < TABSIZ_SPACE; i++) {
|
|
if (dt[i] == 0)
|
|
break;
|
|
else
|
|
tabsiz++;
|
|
}
|
|
tabsiz--;
|
|
return tabsiz;
|
|
}
|
|
|
|
/* Astronomical Almanac table is corrected by adding the expression
|
|
* -0.000091 (ndot + 26)(year-1955)^2 seconds
|
|
* to entries prior to 1955 (AA page K8), where ndot is the secular
|
|
* tidal term in the mean motion of the Moon.
|
|
*
|
|
* Entries after 1955 are referred to atomic time standards and
|
|
* are not affected by errors in Lunar or planetary theory.
|
|
*/
|
|
static double adjust_for_tidacc(double ans, double Y, double tid_acc, double tid_acc0, AS_BOOL adjust_after_1955)
|
|
{
|
|
double B;
|
|
if( Y < 1955.0 || adjust_after_1955) {
|
|
B = (Y - 1955.0);
|
|
ans += -0.000091 * (tid_acc - tid_acc0) * B * B;
|
|
}
|
|
return ans;
|
|
}
|
|
|
|
/* returns tidal acceleration used in swe_deltat() and swe_deltat_ex() */
|
|
double CALL_CONV swe_get_tid_acc()
|
|
{
|
|
return swed.tid_acc;
|
|
}
|
|
|
|
/* function sets tidal acceleration of the Moon.
|
|
* t_acc can be either
|
|
* - the value of the tidal acceleration in arcsec/cty^2
|
|
* of the Moon will be set consistent with that ephemeris.
|
|
* - SE_TIDAL_AUTOMATIC,
|
|
*/
|
|
void CALL_CONV swe_set_tid_acc(double t_acc)
|
|
{
|
|
if (t_acc == SE_TIDAL_AUTOMATIC) {
|
|
swed.tid_acc = SE_TIDAL_DEFAULT;
|
|
swed.is_tid_acc_manual = FALSE;
|
|
return;
|
|
}
|
|
swed.tid_acc = t_acc;
|
|
swed.is_tid_acc_manual = TRUE;
|
|
}
|
|
|
|
void CALL_CONV swe_set_delta_t_userdef(double dt)
|
|
{
|
|
if (dt == SE_DELTAT_AUTOMATIC) {
|
|
swed.delta_t_userdef_is_set = FALSE;
|
|
} else {
|
|
swed.delta_t_userdef_is_set = TRUE;
|
|
swed.delta_t_userdef = dt;
|
|
}
|
|
}
|
|
|
|
int32 swi_guess_ephe_flag()
|
|
{
|
|
int32 iflag = SEFLG_SWIEPH;
|
|
/* if jpl file is open, assume SEFLG_JPLEPH */
|
|
if (swed.jpl_file_is_open) {
|
|
iflag = SEFLG_JPLEPH;
|
|
} else {
|
|
iflag = SEFLG_SWIEPH;
|
|
}
|
|
return iflag;
|
|
}
|
|
|
|
int32 swi_get_tid_acc(double tjd_ut, int32 iflag, int32 denum, int32 *denumret, double *tid_acc, char *serr)
|
|
{
|
|
iflag &= SEFLG_EPHMASK;
|
|
if (swed.is_tid_acc_manual) {
|
|
*tid_acc = swed.tid_acc;
|
|
return iflag;
|
|
}
|
|
if (denum == 0) {
|
|
if (iflag & SEFLG_MOSEPH) {
|
|
*tid_acc = SE_TIDAL_DE404;
|
|
*denumret = 404;
|
|
return iflag;
|
|
}
|
|
if (iflag & SEFLG_JPLEPH) {
|
|
if (swed.jpl_file_is_open) {
|
|
denum = swed.jpldenum;
|
|
}
|
|
}
|
|
/* SEFLG_SWIEPH wanted or SEFLG_JPLEPH failed: */
|
|
if (iflag & SEFLG_SWIEPH) {
|
|
if (swed.fidat[SEI_FILE_MOON].fptr != NULL) {
|
|
denum = swed.fidat[SEI_FILE_MOON].sweph_denum;
|
|
}
|
|
}
|
|
}
|
|
switch(denum) {
|
|
case 200: *tid_acc = SE_TIDAL_DE200; break;
|
|
case 403: *tid_acc = SE_TIDAL_DE403; break;
|
|
case 404: *tid_acc = SE_TIDAL_DE404; break;
|
|
case 405: *tid_acc = SE_TIDAL_DE405; break;
|
|
case 406: *tid_acc = SE_TIDAL_DE406; break;
|
|
case 421: *tid_acc = SE_TIDAL_DE421; break;
|
|
case 422: *tid_acc = SE_TIDAL_DE422; break;
|
|
case 430: *tid_acc = SE_TIDAL_DE430; break;
|
|
case 431: *tid_acc = SE_TIDAL_DE431; break;
|
|
default: denum = SE_DE_NUMBER; *tid_acc = SE_TIDAL_DEFAULT; break;
|
|
}
|
|
*denumret = denum;
|
|
iflag &= SEFLG_EPHMASK;
|
|
return iflag;
|
|
}
|
|
|
|
int32 swi_set_tid_acc(double tjd_ut, int32 iflag, int32 denum, char *serr)
|
|
{
|
|
int32 retc = iflag;
|
|
int32 denumret;
|
|
/* manual tid_acc overrides automatic tid_acc */
|
|
if (swed.is_tid_acc_manual)
|
|
return retc;
|
|
retc = swi_get_tid_acc(tjd_ut, iflag, denum, &denumret, &(swed.tid_acc), serr);
|
|
#if TRACE
|
|
swi_open_trace(NULL);
|
|
if (swi_trace_count < TRACE_COUNT_MAX) {
|
|
if (swi_fp_trace_c != NULL) {
|
|
fputs("\n/*SWE_SET_TID_ACC*/\n", swi_fp_trace_c);
|
|
fprintf(swi_fp_trace_c, " t = %.9f;\n", swed.tid_acc);
|
|
fprintf(swi_fp_trace_c, " swe_set_tid_acc(t);\n");
|
|
fputs(" printf(\"swe_set_tid_acc: %f\\t\\n\", ", swi_fp_trace_c);
|
|
fputs("t);\n", swi_fp_trace_c);
|
|
fflush(swi_fp_trace_c);
|
|
}
|
|
if (swi_fp_trace_out != NULL) {
|
|
fprintf(swi_fp_trace_out, "swe_set_tid_acc: %f\t\n", swed.tid_acc);
|
|
fflush(swi_fp_trace_out);
|
|
}
|
|
}
|
|
#endif
|
|
return retc;
|
|
}
|
|
|
|
/*
|
|
* 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 calculated and then precessed to the equinox of date, using the
|
|
* default precession model of the Swiss Ephmeris. Afte that,
|
|
* sidereal time is derived.
|
|
* The algoritm provides exact agreement for epoch 1 Jan. 2003 with the
|
|
* definition of sidereal time as given in the IERS Convention 2010.
|
|
*/
|
|
/*#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
|
|
* ERA-based expression for for Greenwich Sidereal Time (GST) based
|
|
* on the IAU 2006 precession and IAU 2000A_R06 nutation
|
|
* ftp://maia.usno.navy.mil/conv2010/chapter5/tab5.2e.txt
|
|
*
|
|
* returns sidereal time in hours.
|
|
*
|
|
* program returns sidereal hours since sidereal midnight
|
|
* tjd julian day UT
|
|
* eps obliquity of ecliptic, degrees
|
|
* nut nutation, degrees
|
|
*/
|
|
/* C'_{s,j})_i C'_{c,j})_i */
|
|
#define SIDTNTERM 33
|
|
static const double stcf[SIDTNTERM * 2] = {
|
|
2640.96,-0.39,
|
|
63.52,-0.02,
|
|
11.75,0.01,
|
|
11.21,0.01,
|
|
-4.55,0.00,
|
|
2.02,0.00,
|
|
1.98,0.00,
|
|
-1.72,0.00,
|
|
-1.41,-0.01,
|
|
-1.26,-0.01,
|
|
-0.63,0.00,
|
|
-0.63,0.00,
|
|
0.46,0.00,
|
|
0.45,0.00,
|
|
0.36,0.00,
|
|
-0.24,-0.12,
|
|
0.32,0.00,
|
|
0.28,0.00,
|
|
0.27,0.00,
|
|
0.26,0.00,
|
|
-0.21,0.00,
|
|
0.19,0.00,
|
|
0.18,0.00,
|
|
-0.10,0.05,
|
|
0.15,0.00,
|
|
-0.14,0.00,
|
|
0.14,0.00,
|
|
-0.14,0.00,
|
|
0.14,0.00,
|
|
0.13,0.00,
|
|
-0.11,0.00,
|
|
0.11,0.00,
|
|
0.11,0.00,
|
|
};
|
|
#define SIDTNARG 14
|
|
/* l l' F D Om L_Me L_Ve L_E L_Ma L_J L_Sa L_U L_Ne p_A*/
|
|
static const int stfarg[SIDTNTERM * SIDTNARG] = {
|
|
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 2, -2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 2, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 2, -2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 2, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 1, 2, -2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 1, 2, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 4, -4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 1, -1, 1, 0, -8, 12, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
1, 0, 2, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
1, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 1, -2, 2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 1, -2, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 0, 0, 0, 0, 8, -13, 0, 0, 0, 0, 0, -1,
|
|
0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
2, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
1, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 1, 2, -2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
1, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 4, -2, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 2, -2, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
1, 0, -2, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
1, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
};
|
|
static double sidtime_non_polynomial_part(double tt)
|
|
{
|
|
int i, j;
|
|
double delm[SIDTNARG];
|
|
double dadd, darg;
|
|
/* L Mean anomaly of the Moon.*/
|
|
delm[0] = swe_radnorm(2.35555598 + 8328.6914269554 * tt);
|
|
/* LSU Mean anomaly of the Sun.*/
|
|
delm[1] = swe_radnorm(6.24006013 + 628.301955 * tt);
|
|
/* F Mean argument of the latitude of the Moon. */
|
|
delm[2] = swe_radnorm(1.627905234 + 8433.466158131 * tt);
|
|
/* D Mean elongation of the Moon from the Sun. */
|
|
delm[3] = swe_radnorm(5.198466741 + 7771.3771468121 * tt);
|
|
/* OM Mean longitude of the ascending node of the Moon. */
|
|
delm[4] = swe_radnorm(2.18243920 - 33.757045 * tt);
|
|
/* Planetary longitudes, Mercury through Neptune (Souchay et al. 1999).
|
|
* LME, LVE, LEA, LMA, LJU, LSA, LUR, LNE */
|
|
delm[5] = swe_radnorm(4.402608842 + 2608.7903141574 * tt);
|
|
delm[6] = swe_radnorm(3.176146697 + 1021.3285546211 * tt);
|
|
delm[7] = swe_radnorm(1.753470314 + 628.3075849991 * tt);
|
|
delm[8] = swe_radnorm(6.203480913 + 334.0612426700 * tt);
|
|
delm[9] = swe_radnorm(0.599546497 + 52.9690962641 * tt);
|
|
delm[10] = swe_radnorm(0.874016757 + 21.3299104960 * tt);
|
|
delm[11] = swe_radnorm(5.481293871 + 7.4781598567 * tt);
|
|
delm[12] = swe_radnorm(5.321159000 + 3.8127774000 * tt);
|
|
/* PA General accumulated precession in longitude. */
|
|
delm[13] = (0.02438175 + 0.00000538691 * tt) * tt;
|
|
dadd = -0.87 * sin(delm[4]) * tt;
|
|
for (i = 0; i < SIDTNTERM; i++) {
|
|
darg = 0;
|
|
for (j = 0; j < SIDTNARG; j++) {
|
|
darg += stfarg[i * SIDTNARG + j] * delm[j];
|
|
}
|
|
dadd += stcf[i * 2] * sin(darg) + stcf[i * 2 + 1] * cos(darg);
|
|
}
|
|
dadd /= (3600.0 * 1000000.0);
|
|
return dadd;
|
|
}
|
|
|
|
/*
|
|
* SEMOD_SIDT_IAU_2006
|
|
* N. Capitaine, P.T. Wallace, and J. Chapront, "Expressions for IAU 2000
|
|
* precession quantities", 2003, A&A 412, 567-586 (2003), p. 582.
|
|
* This is a "short" term model, that can be combined with other models
|
|
*/
|
|
/*#define SIDT_IERS_CONV_2010 TRUE*/
|
|
/* sidtime_long_term() is not used between the following two dates */
|
|
#define SIDT_LTERM_T0 2396758.5 /* 1 Jan 1850 */
|
|
#define SIDT_LTERM_T1 2469807.5 /* 1 Jan 2050 */
|
|
#define SIDT_LTERM_OFS0 (0.000378172 / 15.0)
|
|
#define SIDT_LTERM_OFS1 (0.001385646 / 15.0)
|
|
double CALL_CONV swe_sidtime0(double tjd, double eps, double nut)
|
|
{
|
|
double jd0; /* Julian day at midnight Universal Time */
|
|
double secs; /* Time of day, UT seconds since UT midnight */
|
|
double eqeq, jd, tu, tt, msday, jdrel;
|
|
double gmst, dadd;
|
|
int prec_model_short = swed.astro_models[SE_MODEL_PREC_SHORTTERM];
|
|
int sidt_model = swed.astro_models[SE_MODEL_SIDT];
|
|
if (prec_model_short == 0) prec_model_short = SEMOD_PREC_DEFAULT_SHORT;
|
|
if (sidt_model == 0) sidt_model = SEMOD_SIDT_DEFAULT;
|
|
swi_init_swed_if_start();
|
|
if (sidt_model == SEMOD_SIDT_LONGTERM) {
|
|
if (tjd <= SIDT_LTERM_T0 || tjd >= SIDT_LTERM_T1) {
|
|
gmst = sidtime_long_term(tjd, eps, nut);
|
|
if (tjd <= SIDT_LTERM_T0)
|
|
gmst -= SIDT_LTERM_OFS0;
|
|
else if (tjd >= SIDT_LTERM_T1)
|
|
gmst -= SIDT_LTERM_OFS1;
|
|
if (gmst >= 24) gmst -= 24;
|
|
if (gmst < 0) gmst += 24;
|
|
goto sidtime_done;
|
|
}
|
|
}
|
|
/* Julian day at given UT */
|
|
jd = tjd;
|
|
jd0 = floor(jd);
|
|
secs = tjd - jd0;
|
|
if( secs < 0.5 ) {
|
|
jd0 -= 0.5;
|
|
secs += 0.5;
|
|
} else {
|
|
jd0 += 0.5;
|
|
secs -= 0.5;
|
|
}
|
|
secs *= 86400.0;
|
|
tu = (jd0 - J2000)/36525.0; /* UT1 in centuries after J2000 */
|
|
if (sidt_model == SEMOD_SIDT_IERS_CONV_2010 || sidt_model == SEMOD_SIDT_LONGTERM) {
|
|
/* ERA-based expression for Greenwich Sidereal Time (GST) based
|
|
* on the IAU 2006 precession */
|
|
jdrel = tjd - J2000;
|
|
tt = (tjd + swe_deltat_ex(tjd, -1, NULL) - J2000) / 36525.0;
|
|
gmst = swe_degnorm((0.7790572732640 + 1.00273781191135448 * jdrel) * 360);
|
|
gmst += (0.014506 + tt * (4612.156534 + tt * (1.3915817 + tt * (-0.00000044 + tt * (-0.000029956 + tt * -0.0000000368))))) / 3600.0;
|
|
dadd = sidtime_non_polynomial_part(tt);
|
|
gmst = swe_degnorm(gmst + dadd);
|
|
/*printf("gmst iers=%f \n", gmst);*/
|
|
gmst = gmst / 15.0 * 3600.0;
|
|
/* sidt_model == SEMOD_SIDT_IAU_2006, older standards according to precession model */
|
|
} else if (sidt_model == SEMOD_SIDT_IAU_2006) {
|
|
tt = (jd0 + swe_deltat_ex(jd0, -1, NULL) - J2000)/36525.0; /* TT in centuries after J2000 */
|
|
gmst = (((-0.000000002454*tt - 0.00000199708)*tt - 0.0000002926)*tt + 0.092772110)*tt*tt + 307.4771013*(tt-tu) + 8640184.79447825*tu + 24110.5493771;
|
|
/* mean solar days per sidereal day at date tu;
|
|
* for the derivative of gmst, we can assume UT1 =~ TT */
|
|
msday = 1 + ((((-0.000000012270*tt - 0.00000798832)*tt - 0.0000008778)*tt + 0.185544220)*tt + 8640184.79447825)/(86400.*36525.);
|
|
gmst += msday * secs;
|
|
/* SEMOD_SIDT_IAU_1976 */
|
|
} else { /* IAU 1976 formula */
|
|
/* Greenwich Mean Sidereal Time at 0h UT of date */
|
|
gmst = (( -6.2e-6*tu + 9.3104e-2)*tu + 8640184.812866)*tu + 24110.54841;
|
|
/* mean solar days per sidereal day at date tu, = 1.00273790934 in 1986 */
|
|
msday = 1.0 + ((-1.86e-5*tu + 0.186208)*tu + 8640184.812866)/(86400.*36525.);
|
|
gmst += msday * secs;
|
|
}
|
|
/* Local apparent sidereal time at given UT at Greenwich */
|
|
eqeq = 240.0 * nut * cos(eps * DEGTORAD);
|
|
gmst = gmst + eqeq /* + 240.0*tlong */;
|
|
/* Sidereal seconds modulo 1 sidereal day */
|
|
gmst = gmst - 86400.0 * floor( gmst/86400.0 );
|
|
/* return in hours */
|
|
gmst /= 3600;
|
|
goto sidtime_done;
|
|
sidtime_done:
|
|
#ifdef TRACE
|
|
swi_open_trace(NULL);
|
|
if (swi_trace_count < TRACE_COUNT_MAX) {
|
|
if (swi_fp_trace_c != NULL) {
|
|
fputs("\n/*SWE_SIDTIME0*/\n", swi_fp_trace_c);
|
|
fprintf(swi_fp_trace_c, " tjd = %.9f;", tjd);
|
|
fprintf(swi_fp_trace_c, " eps = %.9f;", eps);
|
|
fprintf(swi_fp_trace_c, " nut = %.9f;\n", nut);
|
|
fprintf(swi_fp_trace_c, " t = swe_sidtime0(tjd, eps, nut);\n");
|
|
fputs(" printf(\"swe_sidtime0: %f\\tsidt = %f\\teps = %f\\tnut = %f\\t\\n\", ", swi_fp_trace_c);
|
|
fputs("tjd, t, eps, nut);\n", swi_fp_trace_c);
|
|
fflush(swi_fp_trace_c);
|
|
}
|
|
if (swi_fp_trace_out != NULL) {
|
|
fprintf(swi_fp_trace_out, "swe_sidtime0: %f\tsidt = %f\teps = %f\tnut = %f\t\n", tjd, gmst, eps, nut);
|
|
fflush(swi_fp_trace_out);
|
|
}
|
|
}
|
|
#endif
|
|
return gmst;
|
|
}
|
|
|
|
void CALL_CONV swe_set_interpolate_nut(AS_BOOL do_interpolate)
|
|
{
|
|
if (swed.do_interpolate_nut == do_interpolate)
|
|
return;
|
|
if (do_interpolate)
|
|
swed.do_interpolate_nut = TRUE;
|
|
else
|
|
swed.do_interpolate_nut = FALSE;
|
|
swed.interpol.tjd_nut0 = 0;
|
|
swed.interpol.tjd_nut2 = 0;
|
|
swed.interpol.nut_dpsi0 = 0;
|
|
swed.interpol.nut_dpsi1 = 0;
|
|
swed.interpol.nut_dpsi2 = 0;
|
|
swed.interpol.nut_deps0 = 0;
|
|
swed.interpol.nut_deps1 = 0;
|
|
swed.interpol.nut_deps2 = 0;
|
|
}
|
|
|
|
/* sidereal time, without eps and nut as parameters.
|
|
* tjd must be UT !!!
|
|
* for more informsation, see comment with swe_sidtime0()
|
|
*/
|
|
double CALL_CONV swe_sidtime(double tjd_ut)
|
|
{
|
|
int i;
|
|
double eps, nutlo[2], tsid;
|
|
double tjde;
|
|
/* delta t adjusted to default tidal acceleration of the moon */
|
|
tjde = tjd_ut + swe_deltat_ex(tjd_ut, -1, NULL);
|
|
swi_init_swed_if_start();
|
|
eps = swi_epsiln(tjde, 0) * RADTODEG;
|
|
swi_nutation(tjde, 0, nutlo);
|
|
for (i = 0; i < 2; i++)
|
|
nutlo[i] *= RADTODEG;
|
|
tsid = swe_sidtime0(tjd_ut, eps + nutlo[1], nutlo[0]);
|
|
return tsid;
|
|
}
|
|
|
|
/* SWISSEPH
|
|
* generates name of ephemeris file
|
|
* file name looks as follows:
|
|
* swephpl.m30, where
|
|
*
|
|
* "sweph" "swiss ephemeris"
|
|
* "pl","mo","as" planet, moon, or asteroid
|
|
* "m" or "_" BC or AD
|
|
*
|
|
* "30" start century
|
|
* tjd = ephemeris file for which julian day
|
|
* ipli = number of planet
|
|
* fname = ephemeris file name
|
|
*/
|
|
void swi_gen_filename(double tjd, int ipli, char *fname)
|
|
{
|
|
int icty;
|
|
int ncties = (int) NCTIES;
|
|
short gregflag;
|
|
int jmon, jday, jyear, sgn;
|
|
double jut;
|
|
char *sform;
|
|
switch(ipli) {
|
|
case SEI_MOON:
|
|
strcpy(fname, "semo");
|
|
break;
|
|
case SEI_EMB:
|
|
case SEI_MERCURY:
|
|
case SEI_VENUS:
|
|
case SEI_MARS:
|
|
case SEI_JUPITER:
|
|
case SEI_SATURN:
|
|
case SEI_URANUS:
|
|
case SEI_NEPTUNE:
|
|
case SEI_PLUTO:
|
|
case SEI_SUNBARY:
|
|
strcpy(fname, "sepl");
|
|
break;
|
|
case SEI_CERES:
|
|
case SEI_PALLAS:
|
|
case SEI_JUNO:
|
|
case SEI_VESTA:
|
|
case SEI_CHIRON:
|
|
case SEI_PHOLUS:
|
|
strcpy(fname, "seas");
|
|
break;
|
|
default: /* asteroid or planetary moon */
|
|
if (ipli > SE_PLMOON_OFFSET && ipli < SE_AST_OFFSET) {
|
|
sprintf(fname, "sat%ssepm%d.%s", DIR_GLUE, ipli, SE_FILE_SUFFIX);
|
|
} else {
|
|
sform = "ast%d%sse%05d.%s";
|
|
if (ipli - SE_AST_OFFSET > 99999)
|
|
sform = "ast%d%ss%06d.%s";
|
|
sprintf(fname, sform, (ipli - SE_AST_OFFSET) / 1000, DIR_GLUE, ipli - SE_AST_OFFSET, SE_FILE_SUFFIX);
|
|
}
|
|
return; /* asteroids or planetary moons: only one file 3000 bc - 3000 ad */
|
|
/* break; */
|
|
}
|
|
/* century of tjd */
|
|
/* if tjd > 1600 then gregorian calendar */
|
|
if (tjd >= 2305447.5) {
|
|
gregflag = TRUE;
|
|
swe_revjul(tjd, gregflag, &jyear, &jmon, &jday, &jut);
|
|
/* else julian calendar */
|
|
} else {
|
|
gregflag = FALSE;
|
|
swe_revjul(tjd, gregflag, &jyear, &jmon, &jday, &jut);
|
|
}
|
|
/* start century of file containing tjd */
|
|
if (jyear < 0)
|
|
sgn = -1;
|
|
else
|
|
sgn = 1;
|
|
icty = jyear / 100;
|
|
if (sgn < 0 && jyear % 100 != 0)
|
|
icty -=1;
|
|
while(icty % ncties != 0)
|
|
icty--;
|
|
#if 0
|
|
if (icty < BEG_YEAR / 100)
|
|
icty = BEG_YEAR / 100;
|
|
if (icty >= END_YEAR / 100)
|
|
icty = END_YEAR / 100 - ncties;
|
|
#endif
|
|
/* B.C. or A.D. */
|
|
if (icty < 0)
|
|
strcat(fname, "m");
|
|
else
|
|
strcat(fname, "_");
|
|
icty = abs(icty);
|
|
sprintf(fname + strlen(fname), "%02d.%s", icty, SE_FILE_SUFFIX);
|
|
#if 0
|
|
printf("fname %s\n", fname);
|
|
fflush(stdout);
|
|
#endif
|
|
}
|
|
|
|
/**************************************************************
|
|
cut the string s at any char in cutlist; put pointers to partial strings
|
|
into cpos[0..n-1], return number of partial strings;
|
|
if less than nmax fields are found, the first empty pointer is
|
|
set to NULL.
|
|
More than one character of cutlist in direct sequence count as one
|
|
separator only! cut_str_any("word,,,word2",","..) cuts only two parts,
|
|
cpos[0] = "word" and cpos[1] = "word2".
|
|
If more than nmax fields are found, nmax is returned and the
|
|
last field nmax-1 rmains un-cut.
|
|
**************************************************************/
|
|
int swi_cutstr(char *s, char *cutlist, char *cpos[], int nmax)
|
|
{
|
|
int n = 1;
|
|
cpos [0] = s;
|
|
while (*s != '\0') {
|
|
if ((strchr(cutlist, (int) *s) != NULL) && n < nmax) {
|
|
*s = '\0';
|
|
while (*(s + 1) != '\0' && strchr (cutlist, (int) *(s + 1)) != NULL) s++;
|
|
cpos[n++] = s + 1;
|
|
}
|
|
if (*s == '\n' || *s == '\r') { /* treat nl or cr like end of string */
|
|
*s = '\0';
|
|
break;
|
|
}
|
|
s++;
|
|
}
|
|
if (n < nmax) cpos[n] = NULL;
|
|
return (n);
|
|
} /* cutstr */
|
|
|
|
char *swi_right_trim(char *s)
|
|
{
|
|
char *sp = s + strlen(s) - 1;
|
|
// while (isspace((int)(unsigned char) *sp) && sp >= s)
|
|
while (sp >= s && isspace((int)(unsigned char) *sp))
|
|
*sp-- = '\0';
|
|
return s;
|
|
}
|
|
|
|
size_t swi_strnlen(const char *str, size_t n) {
|
|
const char * stop = (char *)memchr(str, '\0', n);
|
|
return stop ? stop - str : n;
|
|
}
|
|
|
|
/*
|
|
* The following C code (by Rob Warnock rpw3@sgi.com) does CRC-32 in
|
|
* BigEndian/BigEndian byte/bit order. That is, the data is sent most
|
|
* significant byte first, and each of the bits within a byte is sent most
|
|
* significant bit first, as in FDDI. You will need to twiddle with it to do
|
|
* Ethernet CRC, i.e., BigEndian/LittleEndian byte/bit order.
|
|
*
|
|
* The CRCs this code generates agree with the vendor-supplied Verilog models
|
|
* of several of the popular FDDI "MAC" chips.
|
|
*/
|
|
static TLS uint32 crc32_table[256];
|
|
/* Initialized first time "crc32()" is called. If you prefer, you can
|
|
* statically initialize it at compile time. [Another exercise.]
|
|
*/
|
|
|
|
uint32 swi_crc32(unsigned char *buf, int len)
|
|
{
|
|
unsigned char *p;
|
|
uint32 crc;
|
|
if (!crc32_table[1]) /* if not already done, */
|
|
init_crc32(); /* build table */
|
|
crc = 0xffffffff; /* preload shift register, per CRC-32 spec */
|
|
for (p = buf; len > 0; ++p, --len)
|
|
crc = (crc << 8) ^ crc32_table[(crc >> 24) ^ *p];
|
|
return ~crc; /* transmit complement, per CRC-32 spec */
|
|
}
|
|
|
|
/*
|
|
* Build auxiliary table for parallel byte-at-a-time CRC-32.
|
|
*/
|
|
#define CRC32_POLY 0x04c11db7 /* AUTODIN II, Ethernet, & FDDI */
|
|
|
|
static void init_crc32(void)
|
|
{
|
|
int32 i, j;
|
|
uint32 c;
|
|
for (i = 0; i < 256; ++i) {
|
|
for (c = i << 24, j = 8; j > 0; --j)
|
|
c = c & 0x80000000 ? (c << 1) ^ CRC32_POLY : (c << 1);
|
|
crc32_table[i] = c;
|
|
}
|
|
}
|
|
|
|
/*******************************************************
|
|
* other functions from swephlib.c;
|
|
* they are not needed for Swiss Ephemeris,
|
|
* but may be useful to former Placalc users.
|
|
********************************************************/
|
|
|
|
/************************************
|
|
normalize argument into interval [0..DEG360]
|
|
*************************************/
|
|
centisec CALL_CONV swe_csnorm(centisec p)
|
|
{
|
|
if (p < 0)
|
|
do { p += DEG360; } while (p < 0);
|
|
else if (p >= DEG360)
|
|
do { p -= DEG360; } while (p >= DEG360);
|
|
return (p);
|
|
}
|
|
|
|
/************************************
|
|
distance in centisecs p1 - p2
|
|
normalized to [0..360[
|
|
**************************************/
|
|
centisec CALL_CONV swe_difcsn (centisec p1, centisec p2)
|
|
{
|
|
return (swe_csnorm(p1 - p2));
|
|
}
|
|
|
|
double CALL_CONV swe_difdegn (double p1, double p2)
|
|
{
|
|
return (swe_degnorm(p1 - p2));
|
|
}
|
|
|
|
/************************************
|
|
distance in centisecs p1 - p2
|
|
normalized to [-180..180[
|
|
**************************************/
|
|
centisec CALL_CONV swe_difcs2n(centisec p1, centisec p2)
|
|
{ centisec dif;
|
|
dif = swe_csnorm(p1 - p2);
|
|
if (dif >= DEG180) return (dif - DEG360);
|
|
return (dif);
|
|
}
|
|
|
|
double CALL_CONV swe_difdeg2n(double p1, double p2)
|
|
{ double dif;
|
|
dif = swe_degnorm(p1 - p2);
|
|
if (dif >= 180.0) return (dif - 360.0);
|
|
return (dif);
|
|
}
|
|
|
|
double CALL_CONV swe_difrad2n(double p1, double p2)
|
|
{ double dif;
|
|
dif = swe_radnorm(p1 - p2);
|
|
if (dif >= TWOPI / 2) return (dif - TWOPI);
|
|
return (dif);
|
|
}
|
|
|
|
/*************************************
|
|
round second, but at 29.5959 always down
|
|
*************************************/
|
|
centisec CALL_CONV swe_csroundsec(centisec x)
|
|
{
|
|
centisec t;
|
|
t = (x + 50) / 100 *100L; /* round to seconds */
|
|
if (t > x && t % DEG30 == 0) /* was rounded up to next sign */
|
|
t = x / 100 * 100L; /* round last second of sign downwards */
|
|
return (t);
|
|
}
|
|
|
|
/*************************************
|
|
double to int32 with rounding, no overflow check
|
|
*************************************/
|
|
int32 CALL_CONV swe_d2l(double x)
|
|
{
|
|
if (x >=0)
|
|
return ((int32) (x + 0.5));
|
|
else
|
|
return (- (int32) (0.5 - x));
|
|
}
|
|
|
|
/*
|
|
* monday = 0, ... sunday = 6
|
|
*/
|
|
int CALL_CONV swe_day_of_week(double jd)
|
|
{
|
|
return (((int) floor (jd - 2433282 - 1.5) %7) + 7) % 7;
|
|
}
|
|
|
|
char *CALL_CONV swe_cs2timestr(CSEC t, int sep, AS_BOOL suppressZero, char *a)
|
|
/* does not suppress zeros in hours or minutes */
|
|
{
|
|
/* static char a[9];*/
|
|
centisec h,m,s;
|
|
strcpy (a, " ");
|
|
a[2] = a [5] = sep;
|
|
t = ((t + 50) / 100) % (24L *3600L); /* round to seconds */
|
|
s = t % 60L;
|
|
m = (t / 60) % 60L;
|
|
h = t / 3600 % 100L;
|
|
if (s == 0 && suppressZero)
|
|
a[5] = '\0';
|
|
else {
|
|
a [6] = (char) (s / 10 + '0');
|
|
a [7] = (char) (s % 10 + '0');
|
|
};
|
|
a [0] = (char) (h / 10 + '0');
|
|
a [1] = (char) (h % 10 + '0');
|
|
a [3] = (char) (m / 10 + '0');
|
|
a [4] = (char) (m % 10 + '0');
|
|
return (a);
|
|
} /* swe_cs2timestr() */
|
|
|
|
char *CALL_CONV swe_cs2lonlatstr(CSEC t, char pchar, char mchar, char *sp)
|
|
{
|
|
char a[10]; /* must be initialized at each call */
|
|
char *aa;
|
|
centisec h,m,s;
|
|
strcpy (a, " ' ");
|
|
/* mask dddEmm'ss" */
|
|
if (t < 0 ) pchar = mchar;
|
|
t = (ABS4 (t) + 50) / 100; /* round to seconds */
|
|
s = t % 60L;
|
|
m = t / 60 % 60L;
|
|
h = t / 3600 % 1000L;
|
|
if (s == 0)
|
|
a[6] = '\0'; /* cut off seconds */
|
|
else {
|
|
a [7] = (char) (s / 10 + '0');
|
|
a [8] = (char) (s % 10 + '0');
|
|
}
|
|
a [3] = pchar;
|
|
if (h > 99) a [0] = (char) (h / 100 + '0');
|
|
if (h > 9) a [1] = (char) (h % 100 / 10 + '0');
|
|
a [2] = (char) (h % 10 + '0');
|
|
a [4] = (char) (m / 10 + '0');
|
|
a [5] = (char) (m % 10 + '0');
|
|
aa = a;
|
|
while (*aa == ' ') aa++;
|
|
strcpy(sp, aa);
|
|
return (sp);
|
|
} /* swe_cs2lonlatstr() */
|
|
|
|
char *CALL_CONV swe_cs2degstr(CSEC t, char *a)
|
|
/* does suppress leading zeros in degrees */
|
|
{
|
|
/* char a[9]; must be initialized at each call */
|
|
centisec h,m,s;
|
|
t = t / 100 % (30L*3600L); /* truncate to seconds */
|
|
s = t % 60L;
|
|
m = t / 60 % 60L;
|
|
h = t / 3600 % 100L; /* only 0..99 degrees */
|
|
sprintf(a, "%2d%s%02d'%02d", h, ODEGREE_STRING, m, s);
|
|
return (a);
|
|
} /* swe_cs2degstr() */
|
|
|
|
/******************************************************************
|
|
* decimal degrees in zodiac to nakshatra position, deg, min, sec *
|
|
* for definition of input see function swe_split_deg().
|
|
* output:
|
|
* ideg degrees,
|
|
* imin minutes,
|
|
* isec seconds,
|
|
* dsecfr fraction of seconds
|
|
* inak nakshatra number;
|
|
******************************************************************/
|
|
static void split_deg_nakshatra(double ddeg, int32 roundflag, int32 *ideg, int32 *imin, int32 *isec, double *dsecfr, int32 *inak)
|
|
{
|
|
double dadd = 0;
|
|
double dnakshsize = 13.33333333333333;
|
|
double ddeghelp = fmod(ddeg, dnakshsize);
|
|
*inak = 1;
|
|
if (ddeg < 0) {
|
|
*inak = -1;
|
|
ddeg = 0;
|
|
}
|
|
// Sheoran "Vedic" ayanamsha: 0 Aries = 3°20 Ashvini
|
|
if ((swed.sidd.sid_mode & SE_SIDM_TRUE_SHEORAN) == SE_SIDM_TRUE_SHEORAN)
|
|
ddeg = swe_degnorm(ddeg + 3.33333333333333);
|
|
if (roundflag & SE_SPLIT_DEG_ROUND_DEG) {
|
|
dadd = 0.5;
|
|
} else if (roundflag & SE_SPLIT_DEG_ROUND_MIN) {
|
|
dadd = 0.5 / 60;
|
|
} else if (roundflag & SE_SPLIT_DEG_ROUND_SEC) {
|
|
dadd = 0.5 / 3600;
|
|
}
|
|
if (roundflag & SE_SPLIT_DEG_KEEP_DEG) {
|
|
if ((int32) (ddeghelp + dadd) - (int32) ddeghelp > 0)
|
|
dadd = 0;
|
|
} else if (roundflag & SE_SPLIT_DEG_KEEP_SIGN) {
|
|
if (ddeghelp + dadd >= dnakshsize)
|
|
dadd = 0;
|
|
}
|
|
ddeg += dadd;
|
|
*inak = (int32) (ddeg / dnakshsize);
|
|
if (*inak == 27) *inak = 0; // with rounding up from 359.9999
|
|
ddeg = fmod(ddeg, dnakshsize);
|
|
*ideg = (int32) ddeg;
|
|
ddeg -= *ideg;
|
|
*imin = (int32) (ddeg * 60);
|
|
ddeg -= *imin / 60.0;
|
|
*isec = (int32) (ddeg * 3600);
|
|
if (!(roundflag & (SE_SPLIT_DEG_ROUND_DEG | SE_SPLIT_DEG_ROUND_MIN | SE_SPLIT_DEG_ROUND_SEC))) {
|
|
*dsecfr = ddeg * 3600 - *isec;
|
|
}
|
|
} /* end split_deg_nakshtra */
|
|
|
|
/************************************************************
|
|
* splitting decimal degrees into (zod.sign,) deg, min, sec. *
|
|
* input:
|
|
* ddeg decimal degrees, ecliptic longitude
|
|
* roundflag by default there is no rounding. if rounding is
|
|
* required, the following bits can be set:
|
|
# define SE_SPLIT_DEG_ROUND_SEC 1
|
|
# define SE_SPLIT_DEG_ROUND_MIN 2
|
|
# define SE_SPLIT_DEG_ROUND_DEG 4
|
|
# define SE_SPLIT_DEG_ZODIACAL 8 * split into zodiac signs
|
|
# define SE_SPLIT_DEG_NAKSHATRA 1024 * split into nakshatras *
|
|
# define SE_SPLIT_DEG_KEEP_SIGN 16 * don't round to next zodiac sign,
|
|
* e.g. 29.9999998 will be rounded
|
|
* to 29°59'59" (or 29°59' or 29°) *
|
|
* or next nakshatra:
|
|
* e.g. 13.3333332 will be rounded
|
|
* to 13°19'59" (or 13°19' or 13°) *
|
|
# define SE_SPLIT_DEG_KEEP_DEG 32 * don't round to next degree
|
|
* e.g. 10.9999999 will be rounded
|
|
* to 10d59'59" (or 10d59' or 10d) *
|
|
* output:
|
|
* ideg degrees,
|
|
* imin minutes,
|
|
* isec seconds,
|
|
* dsecfr fraction of seconds
|
|
* isgn zodiac sign number;
|
|
* or +/- sign
|
|
*
|
|
*********************************************************/
|
|
void CALL_CONV swe_split_deg(double ddeg, int32 roundflag, int32 *ideg, int32 *imin, int32 *isec, double *dsecfr, int32 *isgn)
|
|
{
|
|
double dadd = 0;
|
|
*isgn = 1;
|
|
if (ddeg < 0) {
|
|
*isgn = -1;
|
|
ddeg = -ddeg;
|
|
} else if (roundflag & SE_SPLIT_DEG_NAKSHATRA) {
|
|
split_deg_nakshatra(ddeg, roundflag, ideg, imin, isec, dsecfr, isgn);
|
|
return;
|
|
}
|
|
if (roundflag & SE_SPLIT_DEG_ROUND_DEG) {
|
|
dadd = 0.5;
|
|
} else if (roundflag & SE_SPLIT_DEG_ROUND_MIN) {
|
|
dadd = 0.5 / 60.0;
|
|
} else if (roundflag & SE_SPLIT_DEG_ROUND_SEC) {
|
|
dadd = 0.5 / 3600.0;
|
|
}
|
|
if (roundflag & SE_SPLIT_DEG_KEEP_DEG) {
|
|
if ((int32) (ddeg + dadd) - (int32) ddeg > 0)
|
|
dadd = 0;
|
|
} else if (roundflag & SE_SPLIT_DEG_KEEP_SIGN) {
|
|
if (fmod(ddeg, 30) + dadd >= 30)
|
|
dadd = 0;
|
|
}
|
|
ddeg += dadd;
|
|
if (roundflag & SE_SPLIT_DEG_ZODIACAL) {
|
|
*isgn = (int32) (ddeg / 30);
|
|
if (*isgn == 12) // 360° = 0°
|
|
*isgn = 0;
|
|
ddeg = fmod(ddeg, 30);
|
|
}
|
|
*ideg = (int32) ddeg;
|
|
ddeg -= *ideg;
|
|
*imin = (int32) (ddeg * 60);
|
|
ddeg -= *imin / 60.0;
|
|
*isec = (int32) (ddeg * 3600);
|
|
if (!(roundflag & (SE_SPLIT_DEG_ROUND_DEG | SE_SPLIT_DEG_ROUND_MIN | SE_SPLIT_DEG_ROUND_SEC))) {
|
|
*dsecfr = ddeg * 3600 - *isec;
|
|
}
|
|
} /* end split_deg */
|
|
|
|
double swi_kepler(double E, double M, double ecce)
|
|
{
|
|
double dE = 1, E0;
|
|
double x;
|
|
/* simple formula for small eccentricities */
|
|
if (ecce < 0.4) {
|
|
while(dE > 1e-12) {
|
|
E0 = E;
|
|
E = M + ecce * sin(E0);
|
|
dE = fabs(E - E0);
|
|
}
|
|
/* complicated formula for high eccentricities */
|
|
} else {
|
|
while(dE > 1e-12) {
|
|
E0 = E;
|
|
/*
|
|
* Alois 21-jul-2000: workaround an optimizer problem in gcc
|
|
* swi_mod2PI sees very small negative argument e-322 and returns +2PI;
|
|
* we avoid swi_mod2PI for small x.
|
|
*/
|
|
x = (M + ecce * sin(E0) - E0) / (1 - ecce * cos(E0));
|
|
dE = fabs(x);
|
|
if (dE < 1e-2) {
|
|
E = E0 + x;
|
|
} else {
|
|
E = swi_mod2PI(E0 + x);
|
|
dE = fabs(E - E0);
|
|
}
|
|
}
|
|
}
|
|
return E;
|
|
}
|
|
|
|
void swi_FK4_FK5(double *xp, double tjd)
|
|
{
|
|
AS_BOOL correct_speed = TRUE;
|
|
if (xp[0] == 0 && xp[1] == 0 && xp[2] == 0)
|
|
return;
|
|
/* with zero speed, we assume that it should be really zero */
|
|
if (xp[3] == 0)
|
|
correct_speed = FALSE;
|
|
swi_cartpol_sp(xp, xp);
|
|
/* according to Expl.Suppl., p. 167f. */
|
|
xp[0] += (0.035 + 0.085 * (tjd - B1950) / 36524.2198782) / 3600 * 15 * DEGTORAD;
|
|
if (correct_speed)
|
|
xp[3] += (0.085 / 36524.2198782) / 3600 * 15 * DEGTORAD;
|
|
swi_polcart_sp(xp, xp);
|
|
}
|
|
|
|
void swi_FK5_FK4(double *xp, double tjd)
|
|
{
|
|
if (xp[0] == 0 && xp[1] == 0 && xp[2] == 0)
|
|
return;
|
|
swi_cartpol_sp(xp, xp);
|
|
/* according to Expl.Suppl., p. 167f. */
|
|
xp[0] -= (0.035 + 0.085 * (tjd - B1950) / 36524.2198782) / 3600 * 15 * DEGTORAD;
|
|
xp[3] -= (0.085 / 36524.2198782) / 3600 * 15 * DEGTORAD;
|
|
swi_polcart_sp(xp, xp);
|
|
}
|
|
|
|
/* function for inhouse testing only */
|
|
void set_astro_models(char *samod)
|
|
{
|
|
int *pmodel = &(swed.astro_models[0]);
|
|
char *sp, *sp2;
|
|
int i = 0;
|
|
swi_init_swed_if_start();
|
|
sp = samod;
|
|
pmodel[0] = atoi(sp);
|
|
i++;
|
|
while((sp2 = strchr(sp, ',')) != NULL && i < NSE_MODELS) {
|
|
sp = sp2 + 1;
|
|
pmodel[i] = atoi(sp);
|
|
i++;
|
|
}
|
|
}
|
|
|
|
|
|
/*
|
|
* Function for inhouse testing of old SE versions.
|
|
*
|
|
* Values of the following macros are defined in swephexp.h
|
|
* and must be identical to the ones used in the defines below
|
|
D1 SEMOD_DELTAT_STEPHENSON_MORRISON_1984
|
|
D2 SEMOD_DELTAT_STEPHENSON_1997
|
|
D3 SEMOD_DELTAT_STEPHENSON_MORRISON_2004
|
|
D4 SEMOD_DELTAT_ESPENAK_MEEUS_2006
|
|
D5 SEMOD_DELTAT_STEPHENSON_ETC_2016
|
|
|
|
P1 SEMOD_PREC_IAU_1976
|
|
P2 SEMOD_PREC_LASKAR_1986
|
|
P3 SEMOD_PREC_WILL_EPS_LASK
|
|
P4 SEMOD_PREC_WILLIAMS_1994
|
|
P5 SEMOD_PREC_SIMON_1994
|
|
P6 SEMOD_PREC_IAU_2000
|
|
P7 SEMOD_PREC_BRETAGNON_2003
|
|
P8 SEMOD_PREC_IAU_2006
|
|
P9 SEMOD_PREC_VONDRAK_2011
|
|
P10 SEMOD_PREC_OWEN_1990
|
|
P11 SEMOD_PREC_NEWCOMB
|
|
|
|
N1 SEMOD_NUT_IAU_1980
|
|
N2 SEMOD_NUT_IAU_CORR_1987
|
|
N3 SEMOD_NUT_IAU_2000A
|
|
N4 SEMOD_NUT_IAU_2000B
|
|
N5 SEMOD_NUT_WOOLARD
|
|
|
|
B1 SEMOD_BIAS_NONE
|
|
B2 SEMOD_BIAS_IAU2000
|
|
B3 SEMOD_BIAS_IAU2006
|
|
|
|
S1 SEMOD_SIDT_IAU_1976
|
|
S2 SEMOD_SIDT_IAU_2006
|
|
S3 SEMOD_SIDT_IERS_CONV_2010
|
|
S4 SEMOD_SIDT_LONGTERM
|
|
* D P P N B J J S
|
|
*/
|
|
# define AMODELS_SE_1_00 "1,3,1,1,1,0,0,1"
|
|
# define AMODELS_SE_1_64 "2,3,1,1,1,0,0,1"
|
|
# define AMODELS_SE_1_70 "2,8,8,4,2,0,0,2"
|
|
# define AMODELS_SE_1_72 "3,8,8,4,2,0,0,2"
|
|
# define AMODELS_SE_1_77 "4,8,8,4,2,0,0,2"
|
|
# define AMODELS_SE_1_78 "4,9,9,4,2,0,0,2"
|
|
# define AMODELS_SE_1_80 "4,9,9,4,3,0,0,1" /* note sid. time (S)! */
|
|
# define AMODELS_SE_2_00 "4,9,9,4,3,0,0,4"
|
|
# define AMODELS_SE_2_06 "5,9,9,4,3,0,0,4"
|
|
void CALL_CONV swe_set_astro_models(char *samod, int32 iflag)
|
|
{
|
|
double dversion;
|
|
char s[30], *sp;
|
|
swi_init_swed_if_start();
|
|
if (*samod != '\0' && isdigit(*samod)) {
|
|
set_astro_models(samod);
|
|
} else if (*samod == '\0' || strncmp(samod, "SE", 2) == 0) {
|
|
strncpy(s, samod, 20);
|
|
s[20] = '\0';
|
|
if ((sp = strchr(s + 5, '.')) != NULL) // remove second '.' in "SE2.05.01"
|
|
swi_strcpy(sp, sp+1);
|
|
if ((sp = strchr(s + 5, 'b')) != NULL) // remove 'b' in "SE2.05.02b04"
|
|
swi_strcpy(sp, sp+1);
|
|
dversion = atof(s + 2);
|
|
if (dversion == 0)
|
|
dversion = atof(SE_VERSION);
|
|
if (dversion >= 2.06) {
|
|
set_astro_models(AMODELS_SE_2_06);
|
|
} else if (dversion >= 2.01) {
|
|
set_astro_models(AMODELS_SE_2_00);
|
|
} else if (dversion >= 2.00) {
|
|
set_astro_models(AMODELS_SE_2_00);
|
|
if (swi_get_denum(SEI_SUN, iflag) == 431)
|
|
swe_set_tid_acc(SE_TIDAL_DE406);
|
|
} else if (dversion >= 1.80) {
|
|
set_astro_models(AMODELS_SE_1_80);
|
|
swe_set_tid_acc(SE_TIDAL_DE406);
|
|
} else if (dversion >= 1.78) {
|
|
set_astro_models(AMODELS_SE_1_78);
|
|
swe_set_tid_acc(SE_TIDAL_DE406);
|
|
} else if (dversion >= 1.77) {
|
|
set_astro_models(AMODELS_SE_1_77);
|
|
swe_set_tid_acc(SE_TIDAL_DE406);
|
|
} else if (dversion >= 1.72) {
|
|
set_astro_models(AMODELS_SE_1_72);
|
|
swe_set_tid_acc(-25.7376);
|
|
} else if (dversion >= 1.70) {
|
|
set_astro_models(AMODELS_SE_1_70);
|
|
swe_set_tid_acc(-25.7376);
|
|
} else if (dversion >= 1.64) {
|
|
set_astro_models(AMODELS_SE_1_64);
|
|
swe_set_tid_acc(-25.7376);
|
|
} else {
|
|
set_astro_models(AMODELS_SE_1_00);
|
|
swe_set_tid_acc(-25.7376);
|
|
}
|
|
}
|
|
}
|
|
|
|
/* function for inhouse testing only */
|
|
static void get_precession_model(int precmod, int32 iflag, char *s)
|
|
{
|
|
if (precmod == 0)
|
|
precmod = SEMOD_PREC_DEFAULT;
|
|
if (iflag & SEFLG_JPLEPH) {
|
|
if (iflag & SEFLG_JPLHOR) {
|
|
strcpy(s, "IAU 1976 (Lieske) / Owen 1990 before 1799");
|
|
return;
|
|
}
|
|
if (iflag & SEFLG_JPLHOR_APPROX) {
|
|
strcpy(s, "Vondrak 2011 / IAU 1976 (Lieske) before 1962 / Owen 1990 before 1799");
|
|
return;
|
|
}
|
|
}
|
|
switch(precmod) {
|
|
case SEMOD_PREC_IAU_1976:
|
|
strcpy(s, "IAU 1976 (Lieske)");
|
|
break;
|
|
case SEMOD_PREC_IAU_2000:
|
|
strcpy(s, "IAU 2000 (Lieske 1976, Mathews 2002)");
|
|
break;
|
|
case SEMOD_PREC_IAU_2006:
|
|
strcpy(s, "IAU 2006 (Capitaine & alii)");
|
|
break;
|
|
case SEMOD_PREC_BRETAGNON_2003:
|
|
strcpy(s, "Bretagnon 2003");
|
|
break;
|
|
case SEMOD_PREC_LASKAR_1986:
|
|
strcpy(s, "Laskar 1986");
|
|
break;
|
|
case SEMOD_PREC_SIMON_1994:
|
|
strcpy(s, "Simon 1994");
|
|
break;
|
|
case SEMOD_PREC_WILLIAMS_1994:
|
|
strcpy(s, "Williams 1994");
|
|
break;
|
|
case SEMOD_PREC_WILL_EPS_LASK:
|
|
strcpy(s, "Williams 1994 / Epsilon Laskar 1986");
|
|
break;
|
|
case SEMOD_PREC_OWEN_1990:
|
|
strcpy(s, "Owen 1990");
|
|
break;
|
|
case SEMOD_PREC_NEWCOMB:
|
|
strcpy(s, "Newcomb 1895");
|
|
break;
|
|
case SEMOD_PREC_VONDRAK_2011:
|
|
strcpy(s, "Vondrák 2011");
|
|
break;
|
|
default:
|
|
break;
|
|
}
|
|
}
|
|
|
|
/* function for inhouse testing only */
|
|
static void get_deltat_model(int dtmod, char *s)
|
|
{
|
|
if (dtmod == 0)
|
|
dtmod = SEMOD_DELTAT_DEFAULT;
|
|
switch(dtmod) {
|
|
case SEMOD_DELTAT_ESPENAK_MEEUS_2006:
|
|
strcpy(s, "Espenak/Meeus 2006 (before 1633)");
|
|
break;
|
|
case SEMOD_DELTAT_STEPHENSON_MORRISON_2004:
|
|
strcpy(s, "Stephenson/Morrison 2004 (before 1600)");
|
|
break;
|
|
case SEMOD_DELTAT_STEPHENSON_1997:
|
|
strcpy(s, "Stephenson 1997 (before 1600)");
|
|
break;
|
|
case SEMOD_DELTAT_STEPHENSON_MORRISON_1984:
|
|
strcpy(s, "Stephenson/Morrison 1984 (before 1600)");
|
|
break;
|
|
case SEMOD_DELTAT_STEPHENSON_ETC_2016:
|
|
strcpy(s, "Stephenson/Morrison/Hohenkerk 2016 (before 1955)");
|
|
break;
|
|
}
|
|
}
|
|
|
|
/* function for inhouse testing only */
|
|
static void get_nutation_model(int nutmod, int32 iflag, char *s)
|
|
{
|
|
int jplhormod = swed.astro_models[SE_MODEL_JPLHOR_MODE];
|
|
int jplhoramod = swed.astro_models[SE_MODEL_JPLHORA_MODE];
|
|
if (jplhormod == 0)
|
|
jplhormod = SEMOD_JPLHOR_DEFAULT;
|
|
if (jplhoramod == 0)
|
|
jplhoramod = SEMOD_JPLHORA_DEFAULT;
|
|
if (nutmod == 0)
|
|
nutmod = SEMOD_NUT_DEFAULT;
|
|
switch(nutmod) {
|
|
case SEMOD_NUT_WOOLARD:
|
|
strcpy(s, "Woolard 1953");
|
|
break;
|
|
case SEMOD_NUT_IAU_1980:
|
|
strcpy(s, "IAU 1980 (Wahr)");
|
|
break;
|
|
case SEMOD_NUT_IAU_CORR_1987:
|
|
strcpy(s, "Herring 1986");
|
|
break;
|
|
case SEMOD_NUT_IAU_2000A:
|
|
strcpy(s, "IAU 2000A (Mathews)");
|
|
break;
|
|
case SEMOD_NUT_IAU_2000B:
|
|
strcpy(s, "IAU 2000B (Mathews)");
|
|
break;
|
|
}
|
|
if (iflag & SEFLG_JPLEPH) {
|
|
if (iflag & SEFLG_JPLHOR)
|
|
strcpy(s, "IAU 1980 (Wahr)");
|
|
if (iflag & SEFLG_JPLHOR) {
|
|
strcat(s, "\n+ daily corrections to dpsi/deps 1962-today");
|
|
if (jplhormod == SEMOD_JPLHOR_LONG_AGREEMENT)
|
|
strcat(s, "\n good agreement with JPL Horizons between 1800 and today");
|
|
else
|
|
strcat(s, "\n defaults to SEFLG_JPLEPH_APPROX before 1962");
|
|
} else if (iflag & SEFLG_JPLHOR_APPROX){
|
|
strcat(s, "\n+ some corrections, approximating JPL Horizons");
|
|
if (jplhoramod == SEMOD_JPLHORA_1)
|
|
strcat(s, " (SEMOD_JPLHORA_1)");
|
|
else if (jplhoramod == SEMOD_JPLHORA_2)
|
|
strcat(s, " (SEMOD_JPLHORA_2)");
|
|
else
|
|
strcat(s, " (SEMOD_JPLHORA_3)");
|
|
}
|
|
}
|
|
}
|
|
|
|
/* function for inhouse testing only */
|
|
static void get_frame_bias_model(int biasmod, char *s)
|
|
{
|
|
if (biasmod == 0)
|
|
biasmod = SEMOD_BIAS_DEFAULT;
|
|
switch(biasmod) {
|
|
case SEMOD_BIAS_IAU2000:
|
|
strcpy(s, "IAU 2000");
|
|
break;
|
|
case SEMOD_BIAS_IAU2006:
|
|
strcpy(s, "IAU 2006");
|
|
break;
|
|
case SEMOD_BIAS_NONE:
|
|
strcpy(s, "none");
|
|
break;
|
|
}
|
|
}
|
|
|
|
/* function for inhouse testing only */
|
|
static void get_sidt_model(int sidtmod, char *s)
|
|
{
|
|
if (sidtmod == 0)
|
|
sidtmod = SEMOD_SIDT_DEFAULT;
|
|
switch(sidtmod) {
|
|
case SEMOD_SIDT_IAU_1976:
|
|
strcpy(s, "IAU 1976");
|
|
break;
|
|
case SEMOD_SIDT_IAU_2006:
|
|
strcpy(s, "IAU 2006 (Capitaine 2003)");
|
|
break;
|
|
case SEMOD_SIDT_IERS_CONV_2010:
|
|
strcpy(s, "IERS Convention 2010");
|
|
break;
|
|
case SEMOD_SIDT_LONGTERM:
|
|
strcpy(s, "IERS Convention 2010 + long-term extension by Astrodienst");
|
|
break;
|
|
}
|
|
}
|
|
|
|
/* function for inhouse testing only */
|
|
void CALL_CONV swe_get_astro_models(char *samod, char *sdet, int32 iflag)
|
|
{
|
|
int i, imod;
|
|
int *pmodel = &(swed.astro_models[0]);
|
|
char s[AS_MAXCH], samod0[AS_MAXCH];
|
|
AS_BOOL list_all_models = FALSE;
|
|
if (samod != NULL) {
|
|
if (strchr(samod, '+') != NULL)
|
|
list_all_models = TRUE;
|
|
swe_set_astro_models(samod, iflag);
|
|
}
|
|
*samod0 = '\0';
|
|
for (i = 0; i < NSE_MODELS; i++) {
|
|
imod = pmodel[i];
|
|
switch(i) {
|
|
case SE_MODEL_PREC_LONGTERM:
|
|
if (imod == SEMOD_PREC_DEFAULT) imod = 0; break;
|
|
case SE_MODEL_PREC_SHORTTERM:
|
|
if (imod == SEMOD_PREC_DEFAULT_SHORT) imod = 0; break;
|
|
case SE_MODEL_NUT:
|
|
if (imod == SEMOD_NUT_DEFAULT) imod = 0; break;
|
|
case SE_MODEL_SIDT:
|
|
if (imod == SEMOD_SIDT_DEFAULT) imod = 0; break;
|
|
case SE_MODEL_BIAS:
|
|
if (imod == SEMOD_BIAS_DEFAULT) imod = 0; break;
|
|
case SE_MODEL_JPLHOR_MODE:
|
|
if (imod == SEMOD_JPLHOR_DEFAULT) imod = 0; break;
|
|
case SE_MODEL_JPLHORA_MODE:
|
|
if (imod == SEMOD_JPLHORA_DEFAULT) imod = 0; break;
|
|
case SE_MODEL_DELTAT:
|
|
if (imod == SEMOD_DELTAT_DEFAULT) imod = 0; break;
|
|
}
|
|
sprintf(samod0 + strlen(samod0), "%d,", imod);
|
|
}
|
|
/*if (samod != NULL)
|
|
strcpy(samod, samod0);*/
|
|
*sdet = '\0';
|
|
if (sdet != NULL) {
|
|
/* JPL ephemeris number and tidal acceleration used with it */
|
|
sprintf(sdet + strlen(sdet), "JPL eph. %d; tidal acc. Moon used by SE: %.4f\n",
|
|
swi_get_denum(SEI_SUN, iflag), swe_get_tid_acc());
|
|
if (iflag & SEFLG_JPLEPH) {
|
|
if (iflag & SEFLG_JPLHOR)
|
|
strcat(sdet, "JPL Horizons method:\n");
|
|
if (iflag & SEFLG_JPLHOR_APPROX)
|
|
strcat(sdet, "JPL Horizons method (approximation):\n");
|
|
} else if (iflag & SEFLG_SWIEPH) {
|
|
strcat(sdet, "Swiss Ephemeris compressed files sepl*/semo*\n");
|
|
} else {
|
|
strcat(sdet, "Moshier semi-analytical approximation\n");
|
|
}
|
|
/* long-term Delta T calculation */
|
|
get_deltat_model(pmodel[SE_MODEL_DELTAT], s);
|
|
sprintf(sdet + strlen(sdet), "Delta T (long-term): %s\n", s);
|
|
/* precession model */
|
|
get_precession_model(pmodel[SE_MODEL_PREC_LONGTERM], iflag, s);
|
|
sprintf(sdet + strlen(sdet), "Precession: %s\n", s);
|
|
if (pmodel[SE_MODEL_PREC_LONGTERM] != pmodel[SE_MODEL_PREC_SHORTTERM] && !(iflag & (SEFLG_JPLHOR | SEFLG_JPLHOR_APPROX))) {
|
|
get_precession_model(pmodel[SE_MODEL_PREC_SHORTTERM], iflag, s);
|
|
sprintf(sdet + strlen(sdet), "+ short-term model: %s\n", s);
|
|
}
|
|
/* nutation */
|
|
get_nutation_model(pmodel[SE_MODEL_NUT], iflag, s);
|
|
sprintf(sdet + strlen(sdet), "Nutation: %s\n", s);
|
|
/* frame bias */
|
|
get_frame_bias_model(pmodel[SE_MODEL_BIAS], s);
|
|
sprintf(sdet + strlen(sdet), "Frame bias: %s\n", s);
|
|
/* sidereal time */
|
|
get_sidt_model(pmodel[SE_MODEL_SIDT], s);
|
|
sprintf(sdet + strlen(sdet), "Sid. time: %s\n", s);
|
|
/* swetest parameters */
|
|
sprintf(sdet + strlen(sdet), "swetest parameters: D P P N B J J S\n");
|
|
sprintf(sdet + strlen(sdet), " -amod%s", samod0);
|
|
sprintf(sdet + strlen(sdet), " -tidacc%f", swe_get_tid_acc());
|
|
strcat(sdet, "\n");
|
|
/* list all available astronomical models */
|
|
if (!list_all_models) {
|
|
sprintf(sdet + strlen(sdet), "For list of all available astronomical models, add a '+' to the version string\n(swetest parameter -amod%s+ or -amod%s+)\n", samod, samod0);
|
|
} else {
|
|
strcat(sdet, "DELTA T MODELS (D)\n");
|
|
for (i = 0; i <= SEMOD_NDELTAT; i++) {
|
|
if (i == SEMOD_DELTAT_DEFAULT) continue;
|
|
sprintf(sdet + strlen(sdet), " (%d)", i);
|
|
if (i == 0) sprintf(sdet + strlen(sdet), " (=%d)", SEMOD_DELTAT_DEFAULT);
|
|
get_deltat_model(i, s);
|
|
sprintf(sdet + strlen(sdet), ": %s\n", s);
|
|
}
|
|
strcat(sdet, "PRECESSION MODELS (P P) (long-term/short-term)\n");
|
|
for (i = 0; i <= SEMOD_NPREC; i++) {
|
|
if (i == SEMOD_PREC_DEFAULT) continue;
|
|
sprintf(sdet + strlen(sdet), " (%d)", i);
|
|
if (i == 0) sprintf(sdet + strlen(sdet), " (=%d)", SEMOD_PREC_DEFAULT);
|
|
get_precession_model(i, iflag, s);
|
|
sprintf(sdet + strlen(sdet), ": %s\n", s);
|
|
}
|
|
strcat(sdet, "NUTATION MODELS (N)\n");
|
|
for (i = 0; i <= SEMOD_NNUT; i++) {
|
|
if (i == SEMOD_NUT_DEFAULT) continue;
|
|
sprintf(sdet + strlen(sdet), " (%d)", i);
|
|
if (i == 0) sprintf(sdet + strlen(sdet), " (=%d)", SEMOD_NUT_DEFAULT);
|
|
get_nutation_model(i, iflag, s);
|
|
sprintf(sdet + strlen(sdet), ": %s\n", s);
|
|
}
|
|
strcat(sdet, "FRAME BIAS MODELS (B)\n");
|
|
for (i = 0; i <= SEMOD_NBIAS; i++) {
|
|
if (i == SEMOD_BIAS_DEFAULT) continue;
|
|
sprintf(sdet + strlen(sdet), " (%d)", i);
|
|
if (i == 0) sprintf(sdet + strlen(sdet), " (=%d)", SEMOD_BIAS_DEFAULT);
|
|
get_frame_bias_model(i, s);
|
|
sprintf(sdet + strlen(sdet), ": %s\n", s);
|
|
}
|
|
strcat(sdet, "JPL HORIZONS MODELS (J) (with SEFLG_JPLEPH|SEFLG_JPLHOR).\n");
|
|
strcat(sdet, " IAU 1980 (Wahr) + daily corrections to dpsi/deps 1962-today.\n");
|
|
strcat(sdet, " (0 (=1): between 1799 and 1962, dpsi/deps of 20-jan-1962 are used.\n");
|
|
strcat(sdet, " For times beyond the dpsi/deps table, the last tabulated values are used.\n");
|
|
strcat(sdet, " Beyond 1799 and 2201, precession Owen 1990 is used..\n");
|
|
strcat(sdet, " Documentation in swephexp.h under 'methods of JPL Horizons'\n");
|
|
strcat(sdet, "JPL HORIZONS APPROXIMATION (J) (with SEFLG_JPLEPH|SEFLG_JPLHORA)\n");
|
|
strcat(sdet, " Documentation in swephexp.h under 'methods of JPL Horizons'\n");
|
|
strcat(sdet, "SIDEREAL TIME MODELS (S)\n");
|
|
for (i = 0; i <= SEMOD_NSIDT; i++) {
|
|
if (i == SEMOD_SIDT_DEFAULT) continue;
|
|
sprintf(sdet + strlen(sdet), " (%d)", i);
|
|
if (i == 0) sprintf(sdet + strlen(sdet), " (=%d)", SEMOD_SIDT_DEFAULT);
|
|
get_sidt_model(i, s);
|
|
sprintf(sdet + strlen(sdet), ": %s\n", s);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
char *swi_strcpy(char *to, char *from)
|
|
{
|
|
char *sp, s[AS_MAXCH];
|
|
if (*from == '\0') {
|
|
*to = '\0';
|
|
return to;
|
|
}
|
|
if (strlen(from) < AS_MAXCH) {
|
|
strcpy(s, from);
|
|
strcpy(to, s);
|
|
} else {
|
|
sp = strdup(from);
|
|
if (sp == NULL) {
|
|
strcpy(to, from);
|
|
} else {
|
|
strcpy(to, sp);
|
|
free(sp);
|
|
}
|
|
}
|
|
return to;
|
|
}
|
|
|
|
#ifdef TRACE
|
|
void swi_open_trace(char *serr)
|
|
{
|
|
swi_trace_count++;
|
|
if (swi_trace_count >= TRACE_COUNT_MAX) {
|
|
if (swi_trace_count == TRACE_COUNT_MAX) {
|
|
if (serr != NULL)
|
|
sprintf(serr, "trace stopped, %d calls exceeded.", TRACE_COUNT_MAX);
|
|
if (swi_fp_trace_out != NULL)
|
|
fprintf(swi_fp_trace_out, "trace stopped, %d calls exceeded.\n", TRACE_COUNT_MAX);
|
|
if (swi_fp_trace_c != NULL)
|
|
fprintf(swi_fp_trace_c, "/* trace stopped, %d calls exceeded. */\n", TRACE_COUNT_MAX);
|
|
}
|
|
return;
|
|
}
|
|
if (swi_fp_trace_c == NULL) {
|
|
char fname[AS_MAXCH];
|
|
#if TRACE == 2
|
|
char *sp, *sp1;
|
|
int ipid;
|
|
#endif
|
|
/* remove(fname_trace_c); */
|
|
strcpy(fname, fname_trace_c);
|
|
#if TRACE == 2
|
|
sp = strchr(fname_trace_c, '.');
|
|
sp1 = strchr(fname, '.');
|
|
# if MSDOS
|
|
ipid = _getpid();
|
|
# else
|
|
ipid = getpid();
|
|
# endif
|
|
sprintf(sp1, "_%d%s", ipid, sp);
|
|
#endif
|
|
if ((swi_fp_trace_c = fopen(fname, FILE_A_ACCESS)) == NULL) {
|
|
if (serr != NULL) {
|
|
sprintf(serr, "could not open trace output file '%s'", fname);
|
|
}
|
|
} else {
|
|
fputs("#include \"sweodef.h\"\n", swi_fp_trace_c);
|
|
fputs("#include \"swephexp.h\"\n\n", swi_fp_trace_c);
|
|
fputs("void main()\n{\n", swi_fp_trace_c);
|
|
fputs(" double tjd, t, nut, eps; int i, ipl, retc; int32 iflag;\n", swi_fp_trace_c);
|
|
fputs(" double armc, geolat, cusp[12], ascmc[10]; int hsys;\n", swi_fp_trace_c);
|
|
fputs(" double xx[6]; int32 iflgret;\n", swi_fp_trace_c);
|
|
fputs(" char s[AS_MAXCH], star[AS_MAXCH], serr[AS_MAXCH];\n", swi_fp_trace_c);
|
|
fflush(swi_fp_trace_c);
|
|
}
|
|
}
|
|
if (swi_fp_trace_out == NULL) {
|
|
char fname[AS_MAXCH];
|
|
#if TRACE == 2
|
|
char *sp, *sp1;
|
|
int ipid;
|
|
#endif
|
|
/* remove(fname_trace_out); */
|
|
strcpy(fname, fname_trace_out);
|
|
#if TRACE == 2
|
|
sp = strchr(fname_trace_out, '.');
|
|
sp1 = strchr(fname, '.');
|
|
# if MSDOS
|
|
ipid = _getpid();
|
|
# else
|
|
ipid = getpid();
|
|
# endif
|
|
sprintf(sp1, "_%d%s", ipid, sp);
|
|
#endif
|
|
if ((swi_fp_trace_out = fopen(fname, FILE_A_ACCESS)) == NULL) {
|
|
if (serr != NULL) {
|
|
sprintf(serr, "could not open trace output file '%s'", fname);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
#endif
|