Irrlicht - from Noob to Pro: Pluto Moshier

#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