#include <math.h>
/* This code was written by Stephen L. Moshier (www.moshier.net)
modified by onkeltorty in 2010*/
#ifndef PLUTO_H_INCLUDED
#define PLUTO_H_INCLUDED
double J2000 = 2451545.0; /* 2000 January 1.5 */
/* Conversion factors between degrees and radians */
double DTR = 1.7453292519943295769e-2;
double RTD = 5.7295779513082320877e1;
double RTS = 2.0626480624709635516e5; /* arc seconds per radian */
double STR = 4.8481368110953599359e-6; /* radians per arc second */
#define PI 3.14159265358979323846
#define TPI 6.2831853071795864769
#define modtp(x) (x - TPI * floor(x/TPI))
#define mod360(x) (x - 360.0 * floor(x/360.0))
static double psstab[153] = {0.0};
static double pcctab[153] = {0.0};
struct orbit
{
char obname[16]; /* name of the object */
double epoch; /* epoch of orbital elements */
double i; /* inclination */
double W; /* longitude of the ascending node */
double w; /* argument of the perihelion */
double a; /* mean distance (semimajor axis) */
double dm; /* daily motion */
double ecc; /* eccentricity */
double M; /* mean anomaly */
double equinox; /* epoch of equinox and ecliptic */
double mag; /* visual magnitude at 1AU from earth and sun */
double sdiam; /* equatorial semidiameter at 1au, arc seconds */
/* The following used by perterbation formulas: */
int (*oelmnt )(); /* address of program to compute elements */
int (*celmnt )(); /* program to correct longitude and radius */
double L; /* computed mean longitude */
double r; /* computed radius vector */
double plat; /* perturbation in ecliptic latitude */
};
#pragma warning(disable : 4244) //Compilerwarnung ausschalten
struct orbit orb = {0.0};
/* Coefficients for the mean longitudes of the planets, from:
* Bretagnon, P., "Theorie du mouvement de l'ensemble des
* planetes. Solution VSOP82," Astronomy and Astrophysics 114,
* 278-288 (1982). Pluto is *not* treated in that article.
*
* The highest degree term is listed first, the constant term
* last. Time is in hundreds of years from J2000. Angles
* are in radians.
*/
static double mls[] = {
3.100e-11, -9.34290e-8, 2608.79031415742, 4.40260884240, /* Mercury */
-3.038e-11, 2.87555e-8, 1021.32855462110, 3.17614669689, /* Venus */
7.3e-13, -9.91890e-8, 628.30758491800, 1.75347031435, /* etc. */
-5.057e-11, 4.54761e-8, 334.06124314923, 6.20348091341,
7.482e-11, -1.4837133e-6, 52.96909650946, 0.59954649739,
-3.3330e-10, 3.6659741e-6, 21.32990954380, 0.87401675650,
1.0450e-10, -8.48828e-8, 7.47815985673, 5.48129387159,
-4.340e-11, 1.66346e-8, 3.81330356378, 5.31188628676,
5.94686e-7, -1.975706e-4, 2.533920996, 4.16934609, /* Pluto */
};
/* Numerical tables for the orbit of Pluto, by Stephen L. Moshier.
*
* Accuracy over the interval 1400 B.C. to 3000 A.D., compared
* to the DE102:
*
* Peak Error R.M.S. Error
* Longitude 18" 3.5"
* Latitude 13" 2.8"
* Radius .00087au .00022au
*/
/* Polynomial coefficients for the heliocentric mean orbital
* elements of Pluto referenced to the equinox and ecliptic
* of J2000.
*
* The polynomials are least squares fits to a series of 105
* osculating orbital element sets derived from the DE102
* numerical integration.. The reference element sets were
* spaced 15,400 days apart starting at JD 1206200.5,
* extending from about 1400 B.C. to 3000 A.D.
*
* Differential corrections and a rotation to J2000 were applied
* to the DE102 orbital elements, as recommended by Newhall et al
* (Astronomy and Astrophysics 125, 150-167 (1983)).
*
* The highest degree term is listed first, the constant term last.
* Time is measured in thousands of years from J2000. Angles are
* in degrees, distance in au. If these mean elements are used without
* adding in the periodic terms given below, the maximum longitude
* error is about a tenth of a degree.
*
*/
static double elpluto[] = {
/*mean distance*/
0.000289, -0.001204, -0.011164, -0.006889, 0.043852, 39.544625,
/*inclination*/
0.000430, 0.002862, 0.005871, 0.002760, -0.002756, 17.139804,
/*node*/
-0.001108, -0.007646, -0.017331, -0.015068, -0.079643, 110.308843,
/*arg perihelion*/
-0.000018, 0.001940, 0.004109, -0.021246, -0.044277, 113.794573,
/*eccentricity*/
0.000010, 0.000011, -0.000115, -0.000111, 0.000651, 0.249084,
/*mean anomaly*/
-0.007562, -0.061952, -0.106116, -1.149059, 1452.063327, 374.813282,
};
/*perihelion*/
/*
-0.001126, -0.005707, -0.013222, -0.036313, -0.123921, 224.103416,
*/
/*daily motion*/
/*
-0.00000004, 0.00000018, 0.00000166, 0.00000101, -0.00000659, 0.00396357,
*/
/* Mean longitude of Pluto used for periodic perturbations: */
/*
-0.008689, -0.067659, -0.119338, -1.185373, 1451.939406, 238.916698,
*/
/* Data structure of an item in the perturbation table below.
* harms[] gives the harmonic of each planet's orbital frequency.
* Data arrays of polynomial coefficients corresponding to longitude,
* latitude, and radius follow this. Each perturbation term
* is to be calculated by
*
* (c T^2 + c T + c ) cos(w) + (s T^2 + s T + s ) sin(w)
* 2 1 0 2 1 0
*
* where w is the sum
*
* NTRIG-1
* -
* > harms[i] L[i]
* -
* i=0
*
* and the L[i] are the mean longitudes of the planets computed
* by the polynomial expansions given in the following reference:
* Bretagnon, P., "Theorie du mouvement de l'ensemble des
* planetes. Solution VSOP82," Astronomy and Astrophysics 114,
* 278-288 (1982). The mean longitude of Pluto is computed
* by the least squares fit given above.
*
* The time variable T is in thousands of years from J2000.
*/
/* Total number of perturbation terms in the table */
#define NAPP 26
/* Highest harmonic present */
#define NV 7
/* Number of planets represented in the array harms[] */
#define NTRIG 5
/* First planet (Mercury = 1) */
#define FTRIG 5
/* Number of polynomial terms (1 + degree of polynomial)*/
#define NPOL 3
/* Number of elements per row in the table */
#define PERTSIZ (NTRIG + 2*NPOL)
struct term {
char harms[NTRIG];
short longs[6*NPOL];
};
/* For each component (longitude, latitude, radius), the
* coefficients were derived by a simultaneous least squares fit
* of all the listed perturbations. The reference data were
* 8,059 positions from the DE102, spaced 200 days apart.
* Corrections given by Newhall et al were applied to the DE102.
*
* The combinations of planetary harmonics were determined by
* examining the Fourier spectrum of the error and selecting
* those combinations that matched the frequencies of the highest
* spectral peaks.
*
* The least squares procedure implicitly orthogonalizes all
* the basis vectors to take account of their cross correlations.
* This means that if the frequencies of two terms are spaced
* too closely, there is a numerical interaction between their
* coefficients. For example, the Uranus term {0, 0, 1, 0, 0}
* has two closely spaced lines on either side of it. The 4,400
* year baseline is not long enough to separate these frequencies
* as uncorrelated sine waves, so the physical intuition of distinct
* Fourier components becomes blurred.
*
* About 1 arc second of accuracy is lost by rounding the
* coefficients to the nearest 0.1".
*/
static struct term perts[] = {
/* J S U N P */
{{ 0, 0, 0, 0, 0},
/* (T^2 + T + 1)cos + (T^2 + T + 1)sin */
{ 325, 715, 406, 0, 0, 0, /* longitude, units of .1" */
-30, -70, -16, 0, 0, 0, /* latitude, units of .1" */
-1533, -128,-11459, 0, 0, 0}}, /* radius, units of 1" in au */
{{ 0, 0,-3, 4, 3}, /* frequency = 42.0 radians/1000yr */
{ -726, -2327, -1242, 715, 423, -1048,
153, 301, -16, -55, -1, 111,
-1917, -3775, 197, 1383, 2928, 435}},
{{ 0, 0,-4, 5, 4}, /* 71.0*/
{ -202, 590, 1054, -444, -1268, 79,
60, 40, -120, 38, 201, 80,
-369, 205, 1186, -289, -1391, -143}},
{{ 0, 0, 0, 1,-1}, /* 127.9*/
{ 44, -150, 140, -16, -130, -398,
-8, -20, 91, -5, 13, 59,
135, 507, 902, 65, -69, 110}},
{{ 0, 0,-2, 1, 5}, /* 152.7*/
{ -34, 54, 145, 82, 241, 62,
12, 21, -6, -10, -44, -28,
-74, -129, 78, 32, 233, 165}},
{{ 0, 0, 2, 0,-5}, /* 228.7*/
{ -377, -1228, -588, 160, -250, -313,
78, 25, -216, 107, 523, 367,
-215, 835, 618, -800, -2316, -787}},
{{ 0, 0, 0, 0, 1}, /* 253.4*/
{ -371, -785, -322, -919, -1993, -221,
-82, -217, -66, 467, 1204, 370,
1740, 4366, -1442, -744, -1063, -2304}},
{{ 0, 0, 2, 0,-7}, /* 278.1*/
{ 142, 32, -604, 199, 1076, 941,
-134, -276, 45, -63, -463, -356,
317, 2053, 2020, -282, -193, 1220}},
{{ 0, 0, 0, 1, 0}, /* 381.3*/
{ -60, -577, -582, 150, 296, -575,
46, 234, 131, -13, 32, 251,
-261, -1992, -1283, 262, 16, -2270}},
{{ 0, 0,-1, 1, 3}, /* 393.7*/
{ -83, -482, -63, -68, 162, 798,
19, 150, 114, 25, -1, -196,
-292, -1759, -621, -239, 400, 2532}},
{{ 0, 0,-2, 2, 1}, /* 479.6*/
{ 300, 2957, 258, 428, -690, -5147,
-286, 554, 3786, 246, 2190, 319,
-268, -5428, -5800, -988, -1795, 6806}},
{{ 0, 0, 1, 0,-1}, /* 494.4*/
{ 1562, 1625, -4258, -1021, -4956, -1376,
-894, -3893, -1040, -1045, -970, 3097,
-402, 4298, 7810, 2832, 8538, -1630}},
{{ 0, 0, 1,-4, 1}, /* 524.1*/
{ -260, -785, -185, 51, -237, -408,
28, 372, 358, -217, -513, 7,
310, 189, -725, 337, 1544, 802}},
{{ 0, 0, 0, 1, 1}, /* 634.7*/
{ 10, 42, 46, -9, -41, 14,
-7, -33, -6, -6, -24, -35,
8, 40, -51, 14, 56, 52}},
{{ 0, 0,-2, 0, 3}, /* 735.5*/
{ 1210, -259,-16372, -851, -9691, -8229,
580, 8002, 8184, 1054, 579,-12993,
-884, -5748, 238, -421, 3354, 12827}},
{{ 0, 0, 1, 0, 0}, /* 747.8*/
{ -4229, -9215, 12810, 6633, 18767, -4094,
5941, 16055, -4871, 3242, 6640,-10726,
-3187, -9671, -654, -4833,-12977, 8115}},
{{ 0, 0, 0,-4, 3}, /* 765.1*/
{ 1222, 6410, -3646, -910, 3674, 9531,
753, -3270, -8034, 1029, 5336, -3396,
-949, 1291, 8179, -810, -5909, 107}},
{{ 0, 0, 0,-4, 2}, /* 1018.5*/
{ 6, 6, -13, -4, -15, 11,
3, 15, -3, 5, 8, -8,
-2, -6, 19, -10, -16, 18}},
{{ 0, 1, 0, 0,-3}, /* 1372.8*/
{ 2, 2, 16, -2, -3, 3,
-2, -4, -4, -2, -2, 11,
1, 2, 0, 3, 5, -39}},
{{ 0, 1, 0, 0,-2}, /* 1626.2*/
{ -1, -2, -34, 1, 1, -34,
0, 1, 15, 1, 2, -6,
1, 0, -93, 0, -1, 90}},
{{ 0, 1, 0, 0,-1}, /* 1879.6*/
{ 1, 2, -1, 0, 0, 136,
0, 1, 6, -1, -2, 4,
-1, 0, 519, 1, 2, 1}},
{{ 0, 1, 0, 0, 0}, /* 2133.0*/
{ 0, 0, -13, -1, -1, 13,
-1, -2, 14, 0, 0, 5,
2, 3, 96, 0, 1, 96}},
{{ 1, 0, 0, 0, 0}, /* 5296.9*/
{ 0, 0, -23, -2, 0, 24,
0, 0, 30, 0, 0, 13,
3, 2, 175, -1, 0, 174}},
{{-1, 0, 0, 0, 3}, /* 4536.7*/
{ 0, 0, 34, 0, 1, -3,
1, 1, -8, 0, 0, -18,
-1, -1, 2, 0, 0, 64}},
{{ 1, 0, 0, 0,-2}, /* 4790.1*/
{ 0, -1, -64, 2, 0, -63,
0, 0, 28, 0, -1, -11,
4, 0, -165, 0, 1, 161}},
{{ 1, 0, 0, 0,-1}, /* 5043.5*/
{ 3, 2, 0, -1, 0, 249,
0, 0, 9, 0, 0, 7,
-5, -1, 941, -2, -2, -1}}
};
static double L[9] = {0.0};
/* Program to calculate the mean longitudes using
* the above coefficients.
*/
void longits( double J, int first, int nobjs, double L[] )
{
double f, T;
double *p;
int i;
T = (J - J2000)/36525.0;
p = &mls[ 4*(first-1) ];
for(i=0; i<nobjs; i++ )
{
f = *p++;
f = T * f + *p++;
f = T * f + *p++;
f = T * f + *p++;
L[i] = f;
}
}
/* Radians to degrees, minutes, seconds
*/
void dms(double x )
{
double s;
int d, m, is;
s = x * RTD;
if( s < 0.0 ) s = -s;
d = s;
s -= d;
s *= 60;
m = s;
s -= m;
s *= 60;
is = s + 0.5;
}
/*
Cephes Math Library Release 2.0: April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
double zatan2( double x, double y )
{
double z, w;
short code;
//double atan();
code = 0;
if( x < 0.0 )
code = 2;
if( y < 0.0 )
code |= 1;
if( x == 0.0 )
{
if( code & 1 )
return( 1.5*PI );
if( y == 0.0 )
return( 0.0 );
return( 0.5*PI );
}
if( y == 0.0 )
{
if( code & 2 )
return( PI );
return( 0.0 );
}
switch( code )
{
case 0: w = 0.0; break;
case 1: w = 2.0 * PI; break;
case 2:
case 3: w = PI; break;
}
z = atan( y/x );
return( w + z );
}
/* Routine to prepare lookup tables of sin and cos of i*L
* for required multiple angles
*/
void psscc( double *pss, double *pcc, double w, int n )
{
double cu, su, cv, sv, s;
int i;
su = sin( w );
cu = cos( w );
*pss++ = su; /* sin(L) */
*pcc++ = cu; /* cos(L) */
sv = 2.0*su*cu;
cv = cu*cu - su*su;
*pss++ = sv; /* sin(2L) */
*pcc++ = cv; /* cos(2L) */
for( i=2; i<n; i++ )
{
s = su*cv + cu*sv;
cv = cu*cv - su*sv;
sv = s;
*pss++ = sv; /* sin( (i+1)*L ) */
*pcc++ = cv; /* cos( (i+1)*L ) */
}
}
/* Subroutine to evaluate a polynomial.
*/
double epoly( double x, double *C, int n )
//double x; /* independent variable */
//double *C; /* coefficients, highest degree first */
//int n; /* degree of the polynomial */
{
double s;
int i;
s = *C++;
for( i=0; i<n; i++ )
{
s = s * x + *C++;
}
return( s );
}
/* Subroutine called by kepler() to set up
* the mean orbital elements for the requested epoch JD
*/
void opluto(struct orbit *o, double JD )
{
double x, T;
double *p;
T = (JD - J2000)/365250.0;
p = elpluto;
o->epoch = JD;
o->a = epoly( T, p, 5 );
p += 6;
o->i = epoly( T, p, 5 );
p += 6;
o->W = epoly( T, p, 5 );
p += 6;
o->w = epoly( T, p, 5 );
p += 6;
o->dm = 0.0;
o->ecc = epoly( T, p, 5 );
p += 6;
x = epoly( T, p, 5 );
o->M = mod360(x);
o->equinox = J2000;
o->r = 0.0;
o->plat = 0.0;
o->L = 0.0;
}
/* Routine to step through the table of perturbations
* and accumulate the terms.
*/
void chewtab( double T, struct term *table, int nobjs, int dpol, int nterms, double cc[], double ss[],
int nsc, int nans, double dvec[] )
//double T; /* time, in units appropriate to the table */
//struct term *table; /* table of perturbations */
//int nobjs; /* number of objects in table */
//int dpol; /* degree of polynomials in table */
//int nterms; /* number of perturbation terms in table */
//double cc[]; /* array of cos( j*L[i] ) */
//double ss[]; /* array of sin( j*L[i] ) */
//int nsc; /* highest harmonic listed in L[] */
//int nans; /* number of outputs */
//double dvec[]; /* output vector (e.g., longitude, latitude, radius) */
{
int i, ians, it, j, k, k1, m, nscm;
char *p;
short *pt;
double su, cu, sv, cv, c, s, y;
for( ians=0; ians<nans; ians++ )
dvec[ians] = 0.0;
cv = 1.0;
sv = 0.0;
for( it=0; it<nterms; it++ )
{
k1 = 0;
p = table->harms;
for( m=0; m<nobjs; m++ )
{
j = (int )*p++;
if( j )
{
k = j;
if( j < 0 )
k = -k;
nscm = nsc * m + k - 1;
su = ss[ nscm ]; /* sin(k*angle) */
if( j < 0 )
su = -su;
cu = cc[ nscm ];
if( k1 == 0 )
{ /* set first angle */
sv = su;
cv = cu;
k1 = 1;
}
else
{ /* combine angles */
y = su*cv + cu*sv;
cv = cu*cv - su*sv;
sv = y;
}
}
}
pt = table->longs;
for( ians=0; ians<nans; ians++ )
{
c = (double )(*pt++);
for( i=0; i<dpol; i++ )
{
c = c * T + (double )(*pt++);
}
s = (double )(*pt++);
for( i=0; i<dpol; i++ )
{
s = s * T + (double )(*pt++);
}
dvec[ians] += c * cv + s * sv;
}
++table;
}
}
/* Subroutine called by kepler(), after solving Kepler's equation,
* to apply perturbations to the spherical ecliptic coordinates.
*/
void cpluto( struct orbit *o )
{
double T;
double ppol[3];
double *pcc, *pss;
int i;
/* Set up table of sines and cosines of harmonics
* of each planet's frequency.
*/
longits( o->epoch, FTRIG, NTRIG, L );
pss = psstab;
pcc = pcctab;
for( i=0; i<NTRIG; i++ )
{
psscc( pss, pcc, L[i], NV );
pcc += NV;
pss += NV;
}
/* Calculate the perturbations. */
T = (o->epoch - J2000) / 365250.0;
chewtab( T, perts, NTRIG, NPOL-1, NAPP, pcctab, psstab, NV, 3, ppol );
o->L += 0.1 * STR * ppol[0];
o->plat = 0.1 * STR * ppol[1];
o->r += STR * ppol[2];
}
void kepler(double J, struct orbit *e, double polar[])
{
double alat, E, M, W, v, temp;
double inclination, ascnode, argperih;
double meandistance, eccent;
double r, coso, sino, sinW, cosv;
int k;
/* Compute orbital elements
*/
//(*(e->oelmnt) )(e,J);
opluto(e,J);
/* Decant the parameters from the data structure
*/
inclination = DTR * e->i;
ascnode = DTR * e->W;
argperih = DTR * e->w;
meandistance = e->a; /* semimajor axis */
eccent = e->ecc;
/* M is proportional to the area swept out by the radius
* vector of a circular orbit during the time between
* perihelion passage and Julian date J.
* It is the mean anomaly at time J.
*/
M = DTR * e->M;
M = modtp(M);
/* By Kepler's second law, M must be equal to
* the area swept out in the same time by an
* elliptical orbit of same total area.
* Integrate the ellipse expressed in polar coordinates
* r = a(1-e^2)/(1 + e cosW)
* with respect to the angle W to get an expression for the
* area swept out by the radius vector. The area is given
* by the mean anomaly; the angle is solved numerically.
*
* The answer is obtained in two steps. We first solve
* Kepler's equation
* M = E - eccent*sin(E)
* for the eccentric anomaly E. Then there is a
* closed form solution for W in terms of E.
*/
E = M; /* Initial guess is same as circular orbit. */
temp = 1.0;
k = 0;
do
{
if( ++k > 100 )
{
printf( "kepler() error\n" );
break;
}
/* The approximate area swept out in the ellipse */
temp = E - eccent * sin(E)
/* ...minus the area swept out in the circle */
- M;
/* ...should be zero. Use the derivative of the error
* to converge to solution by Newton's method.
*/
E -= temp/(1.0 - eccent*cos(E));
}
while( fabs(temp) > 1.0e-11 );
/* The exact formula for the area in the ellipse is
* 2.0*atan(c2*tan(0.5*W)) - c1*eccent*sin(W)/(1+e*cos(W))
* where
* c1 = sqrt( 1.0 - eccent*eccent )
* c2 = sqrt( (1.0-eccent)/(1.0+eccent) ).
* Substituting the following value of W
* yields the exact solution.
*/
temp = sqrt( (1.0+eccent)/(1.0-eccent) );
v = 2.0 * atan( temp * tan(0.5*E) );
/* The true anomaly.
*/
v = modtp(v);
/* Orbital longitude measured from node
* (argument of latitude)
*/
alat = v + argperih;
/* From the equation of the ellipse, get the
* radius from central focus to the object.
*/
cosv = cos( v );
r = meandistance*(1.0-eccent*eccent)/(1.0+eccent*cosv);
/* The heliocentric ecliptic longitude of the object
* is given by
* tan( longitude - ascnode ) = cos( inclination ) * tan( alat ).
*/
coso = cos( alat );
sino = sin( alat );
W = sino * cos( inclination );
E = zatan2( coso, W ) + ascnode;
/* The ecliptic latitude of the object
*/
sinW = sino * sin( inclination );
W = asin(sinW);
/* Apply perturbations
*/
e->L = E;
e->r = r;
//(*(e->celmnt) )(e);
cpluto(e);
E = e->L;
r = e->r;
W += e->plat;
/* Output the polar cooordinates
*/
polar[0] = E; /* longitude */
polar[1] = W; /* latitude */
polar[2] = r; /* radius */
}
#endif // PLUTO_H_INCLUDED
