This is kepler.c in view mode; [Download] [Up]
/************************************************************************/
/* kepler.c */
/* */
/* Program to solve Keplerian orbit given orbital parameters and the */
/* time. Returns Heliocentric equatorial rectangular coordinates of */
/* the object. */
/************************************************************************/
/***** description
*
* $Id: kepler.c,v 1.4 1993/04/22 19:43:55 craig Exp $
*
* This program detects several cases of given orbital elements.
* If a program for perturbations is pointed to, it is called
* to calculate all the elements.
* If there is no program, then the mean longitude is calculated
* from the mean anomaly and daily motion.
* If the daily motion is not given, it is calculated
* by Kepler's law.
* If the eccentricity is given to be 1.0, it means that
* meandistance is really the perihelion distance, as in a comet
* specification, and the orbit is parabolic.
*
* Reference: Taff, L.G., "Celestial Mechanics, A Computational
* Guide for the Practitioner." Wiley, 1985.
*/
/***** modification history
*
* $Log: kepler.c,v $
* Revision 1.4 1993/04/22 19:43:55 craig
* minor changes.
*
* Revision 1.3 1993/04/21 21:20:47 craig
* Changed the path of the satellite.h include.
* Changed ecnsts to pcnsts.
*
* Revision 1.2 1993/04/21 15:16:53 craig
* First working version. Ran through indent and converted to ansi.
* Added hooks for working with the satellite programs.
*
*
*/
/***** include files *****/
#include <math.h>
#include "aaproto.h"
#include "satellite.h"
/***** global variables *****/
/* from cnstinit.c */
extern struct orbit earth; /* orbital elements of the earth */
extern struct MCONSTANTS mcnsts;
extern struct PCONSTANTS pcnsts;
/* from epsiln.c */
extern double eps, coseps, sineps; /* obliquity of ecliptic */
/***** local function prototypes *****/
int embofs (double J, double *E, double *W, double *r);
/**********/
/* kepler */
/**********/
int kepler (double J, struct orbit *e, double *rect, double *polar)
{
int k;
double alat, E, M, W, temp;
double epoch, inclination, ascnode, argperih;
double meandistance, dailymotion, eccent, meananomaly;
double r, coso, sino, cosa;
/* Compute orbital elements if a program for doing so is supplied */
if (e->oelmnt)
{
k = (*(e->oelmnt)) (e, J);
if (k == -1)
{
goto dobs;
}
}
else if (e->celmnt)
{ /* call B & S algorithm */
dobs:
(*(e->celmnt)) (J, polar);
polar[0] = modtp (polar[0]);
E = polar[0]; /* longitude */
e->L = E;
W = polar[1]; /* latitude */
r = polar[2]; /* radius */
e->r = r;
e->epoch = J;
e->equinox = J;
goto kepdon;
}
/* Decant the parameters from the data structure */
epoch = e->epoch;
inclination = e->i;
ascnode = e->W * mcnsts.de2ra;
argperih = e->w;
meandistance = e->a; /* semimajor axis */
dailymotion = e->dm;
eccent = e->ecc;
meananomaly = e->M;
/* Check for parabolic orbit. */
if (eccent == 1.0)
{
/* meandistance = perihelion distance, q epoch = perihelion
passage date */
temp = meandistance * sqrt (meandistance);
W = (J - epoch) * pcnsts.k * 3. / sqrt (2.) / temp;
E = 0.0;
M = 1.0;
while (fabs (M) > 1.0e-11)
{
temp = E * E;
temp = (2.0 * E * temp + W) / (3.0 * (1.0 + temp));
M = temp - E;
if (temp != 0.0)
{
M /= temp;
}
E = temp;
}
r = meandistance * (1.0 + E * E);
M = atan (E);
M = 2.0 * M;
alat = M + mcnsts.de2ra * argperih;
goto parabcon;
}
if (eccent > 1.0)
{
/* The equation of the hyperbola in polar coordinates r, theta is
r = a(e^2 - 1)/(1 + e cos(theta)) so the perihelion distance q
= a(e-1), the "mean distance" a = q/(e-1). */
meandistance = meandistance / (eccent - 1.0);
temp = meandistance * sqrt (meandistance);
W = (J - epoch) * 0.01720209895 / temp;
/* solve M = -E + e sinh E */
E = W / (eccent - 1.0);
M = 1.0;
while (fabs (M) > 1.0e-11)
{
M = -E + eccent * sinh (E) - W;
E += M / (1.0 - eccent * cosh (E));
}
r = meandistance * (-1.0 + eccent * cosh (E));
temp = (eccent + 1.0) / (eccent - 1.0);
M = sqrt (temp) * tanh (0.5 * E);
M = 2.0 * atan (M);
alat = M + mcnsts.de2ra * argperih;
goto parabcon;
}
/* Calculate the daily motion, if it is not given. */
if (dailymotion == 0.0)
{
/* Assumes object in heliocentric orbit is massless. */
dailymotion = pcnsts.k * 180. / mcnsts.pi / (e->a * sqrt (e->a));
}
dailymotion *= (J - epoch);
/* 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 = mcnsts.de2ra * (meananomaly + dailymotion);
M = modtp (M);
/* If mean longitude was calculated, adjust it also for motion since
epoch of elements. */
if (e->L)
{
e->L += dailymotion;
e->L = mod360 (e->L);
}
/* 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;
do
{
/* The approximate area swept out in the ellipse */
/* ...minus the area swept out in the circle (M) */
temp = E - eccent * sin (E) - 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));
W = 2.0 * atan (temp * tan (0.5 * E));
/* The true anomaly. */
W = modtp (W);
meananomaly *= mcnsts.de2ra;
/* Orbital longitude measured from node (argument of latitude) */
if (e->L)
{
alat = (e->L) * mcnsts.de2ra + W - meananomaly - ascnode;
}
else
{
alat = W + mcnsts.de2ra * argperih; /* mean longitude not given */
}
/* From the equation of the ellipse, get the radius from central focus
to the object. */
r = meandistance * (1.0 - eccent * eccent) / (1.0 + eccent * cos (W));
parabcon:
/* The heliocentric ecliptic longitude of the object is given by tan(
longitude - ascnode ) = cos( inclination ) * tan( alat ). */
coso = cos (alat);
sino = sin (alat);
inclination *= mcnsts.de2ra;
W = sino * cos (inclination);
E = matan2 (W, coso) + ascnode;
/* The ecliptic latitude of the object */
W = sino * sin (inclination);
W = asin (W);
/* Apply perturbations, if a program is supplied. */
if (e->celmnt)
{
e->L = E;
e->r = r;
(*(e->celmnt)) (e);
E = e->L;
r = e->r;
W += e->plat;
}
/* If earth, Adjust from earth-moon barycenter to earth by AA page E2
unless orbital elements are calculated by formula. (The Meeus
perturbation formulas include this term for the moon.) */
if ((e == &earth) && (e->oelmnt == 0))
{
embofs (J, &E, &W, &r); /* see below */
}
/* Output the polar cooordinates */
polar[0] = E; /* longitude */
polar[1] = W; /* latitude */
polar[2] = r; /* radius */
kepdon:
/* Convert to rectangular coordinates, using the perturbed latitude. */
rect[2] = r * sin (W);
cosa = cos (W);
rect[1] = r * cosa * sin (E);
rect[0] = r * cosa * cos (E);
/* Convert from heliocentric ecliptic rectangular to heliocentric
equatorial rectangular coordinates by rotating eps radians about
the x axis. */
epsiln (e->equinox);
W = coseps * rect[1] - sineps * rect[2];
M = sineps * rect[1] + coseps * rect[2];
rect[1] = W;
rect[2] = M;
/* Precess the position to ecliptic and equinox of J2000.0 if not
already there. */
precess (rect, e->equinox, 1);
return (0);
}
/********************************************************/
/* embofs */
/* */
/* Adjust position from Earth-Moon barycenter to Earth */
/********************************************************/
/*
* J = Julian day number
* E = Earth's ecliptic longitude (radians)
* W = Earth's ecliptic latitude (radians)
* r = Earth's distance to the Sun (au)
*/
int embofs (double J, double *E, double *W, double *r)
{
double T, M, a, L, B, p;
double smp, cmp, s2mp, c2mp, s2d, c2d, sf, cf;
double s2f, sx, cx, xyz[3];
/* Short series for position of the Moon */
T = (J - pcnsts.J1900) / pcnsts.dapcen;
/* Mean anomaly of moon (MP) */
a = ((1.44e-5 * T + 0.009192) * T + 477198.8491) * T + 296.104608;
a *= mcnsts.de2ra;
smp = sin (a);
cmp = cos (a);
s2mp = 2.0 * smp * cmp; /* sin(2MP) */
c2mp = cmp * cmp - smp * smp; /* cos(2MP) */
/* Mean elongation of moon (D) */
a = ((1.9e-6 * T - 0.001436) * T + 445267.1142) * T + 350.737486;
a = 2.0 * mcnsts.de2ra * a;
s2d = sin (a);
c2d = cos (a);
/* Mean distance of moon from its ascending node (F) */
a = ((-3.e-7 * T - 0.003211) * T + 483202.0251) * T + 11.250889;
a *= mcnsts.de2ra;
sf = sin (a);
cf = cos (a);
s2f = 2.0 * sf * cf; /* sin(2F) */
sx = s2d * cmp - c2d * smp; /* sin(2D - MP) */
cx = c2d * cmp + s2d * smp; /* cos(2D - MP) */
/* Mean longitude of moon (LP) */
L = ((1.9e-6 * T - 0.001133) * T + 481267.8831) * T + 270.434164;
/* Mean anomaly of sun (M) */
M = ((-3.3e-6 * T - 1.50e-4) * T + 35999.0498) * T + 358.475833;
/* Ecliptic longitude of the moon */
L = L
+ 6.288750 * smp
+ 1.274018 * sx
+ 0.658309 * s2d
+ 0.213616 * s2mp
- 0.185596 * sin (mcnsts.de2ra * M)
- 0.114336 * s2f;
/* Ecliptic latitude of the moon */
a = smp * cf;
sx = cmp * sf;
B = 5.128189 * sf
+ 0.280606 * (a + sx) /* sin(MP+F) */
+ 0.277693 * (a - sx) /* sin(MP-F) */
+ 0.173238 * (s2d * cf - c2d * sf); /* sin(2D-F) */
B *= mcnsts.de2ra;
/* Parallax of the moon */
p = 0.950724
+ 0.051818 * cmp
+ 0.009531 * cx
+ 0.007843 * c2d
+ 0.002824 * c2mp;
p *= mcnsts.de2ra;
/* Elongation of Moon from Sun */
L = mcnsts.de2ra * mod360 (L) - (*E - mcnsts.pi);
/* The Moon's position has to be precessed to J2000 */
/* Distance in au */
a = 4.263523e-5 / sin (p);
/* Convert to rectangular ecliptic coordinates */
xyz[2] = a * sin (B);
cx = cos (B);
xyz[1] = a * cx * sin (L);
xyz[0] = a * cx * cos (L);
/* Convert to equatorial */
epsiln (J);
sf = coseps * xyz[1] - sineps * xyz[2];
cf = sineps * xyz[1] + coseps * xyz[2];
xyz[1] = sf;
xyz[2] = cf;
/* Precess the position
* to ecliptic and equinox of J2000.0
*/
precess (xyz, J, 1);
/* Rotate back to ecliptic coordinates */
sf = coseps * xyz[1] + sineps * xyz[2];
cf = -sineps * xyz[1] + coseps * xyz[2];
xyz[1] = sf;
xyz[2] = cf;
/* Convert to polar coordinates */
L = matan2 (xyz[1], xyz[0]);
B = asin (xyz[2] / a);
/* Distance from Earth to Earth-Moon barycenter */
a = 5.18041e-7 / p;
/* Adjust the coordinates of the Earth */
sx = a / *r; /* chord = R * dTheta */
*E += sx * sin (L);
*W -= sx * B;
*r += a * cos (L);
return (0);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.