/* This program calculates orbits of planetary bodies and reduces
* the coordinates of planets or stars to geocentric and topocentric
* place. An effort has been made to use rigorous methods throughout.
*
* References to AA page numbers are to The Astronomical Almanac, 1986
* published by the U.S. Government Printing Office.
*
* The program uses as a default the PLAN404 approximations to DE404
* for planetary positions.
*
* Warning! Your atan2() function may not work the same as the one
* assumed by this program.
* atan2(x,y) computes atan(y/x), result between 0 and 2pi.
*
* S. L. Moshier, November, 1987
*/
/* 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 */
double PI = 3.14159265358979323846;
extern double PI;
/* Standard epochs. Note Julian epochs (J) are measured in
* years of 365.25 days.
*/
double J2000 = 2451545.0; /* 2000 January 1.5 */
double B1950 = 2433282.423; /* 1950 January 0.923 Besselian epoch */
double J1900 = 2415020.0; /* 1900 January 0, 12h UT */
/* Data structures containing orbital elements of
* objects that orbit the sun. See kep.h for the definition.
*/
#include "kep.h"
#include <stdlib.h>
#ifdef _MSC_VER
#if _MSC_VER >= 1000
#include <stdlib.h>
#include <string.h>
#endif
#endif
/* approximate motion of right ascension and declination
* of object, in radians per day
*/
double FAR dradt;
double FAR ddecdt;
/* Space for star description read from a disc file.
*/
struct star fstar;
/* Space for orbit read from a disc file. Entering 99 for the
* planet number yields a prompt for a file name containg ASCII
* strings specifying the elements.
*/
struct orbit forbit;
/* Orbits for each planet. The indicated orbital elements are
* not actually used, since the positions are are now calculated
* from a formula. Magnitude and semidiameter are still used.
*/
/* Programs to compute perturbations. */
extern struct plantbl mer404, ven404, ear404, mar404;
extern struct plantbl jup404, sat404, ura404, nep404, plu404;
struct orbit mercury = {
"Mercury ",
2446800.5, /* January 5.0, 1987 */
7.0048,
48.177,
29.074,
0.387098,
4.09236,
0.205628,
198.7199,
2446800.5,
-0.42,
3.36,
&mer404,
0.0,
0.0,
0.0
};
struct orbit venus = {
"Venus ",
2446800.5,
3.3946,
76.561,
54.889,
0.723329,
1.60214,
0.006757,
9.0369,
2446800.5,
/* Note the calculated apparent visual magnitude for Venus
* is not very accurate.
*/
-4.40,
8.34,
&ven404,
0.0,
0.0,
0.0
};
/* Fixed numerical values will be used for earth if read in from a file
* named earth.orb. See kfiles.c, kep.h.
*/
struct orbit earth = {
"Earth ",
2446800.5,
0.0,
0.0,
102.884,
0.999999,
0.985611,
0.016713,
1.1791,
2446800.5,
-3.86,
0.0,
&ear404,
0.0,
0.0,
0.0
};
extern struct orbit earth;
struct orbit mars = {
"Mars ",
2446800.5,
1.8498,
49.457,
286.343,
1.523710,
0.524023,
0.093472,
53.1893,
2446800.5,
-1.52,
4.68,
&mar404,
0.0,
0.0,
0.0
};
struct orbit jupiter = {
"Jupiter ",
2446800.5,
1.3051,
100.358,
275.129,
5.20265,
0.0830948,
0.048100,
344.5086,
2446800.5,
-9.40,
98.44,
&jup404,
0.0,
0.0,
0.0
};
struct orbit saturn = {
"Saturn ",
2446800.5,
2.4858,
113.555,
337.969,
9.54050,
0.0334510,
0.052786,
159.6327,
2446800.5,
-8.88,
82.73,
&sat404,
0.0,
0.0,
0.0
};
struct orbit uranus = {
"Uranus ",
2446800.5,
0.7738,
73.994,
98.746,
19.2233,
0.0116943,
0.045682,
84.8516,
2446800.5,
-7.19,
35.02,
&ura404,
0.0,
0.0,
0.0
};
struct orbit neptune = {
"Neptune ",
2446800.5,
1.7697,
131.677,
250.623,
30.1631,
0.00594978,
0.009019,
254.2568,
2446800.5,
-6.87,
33.50,
&nep404,
0.0,
0.0,
0.0
};
struct orbit pluto = {
"Pluto ",
2446640.5,
17.1346,
110.204,
114.21,
39.4633,
0.00397570,
0.248662,
355.0554,
2446640.5,
-1.0,
2.07,
&plu404,
0.0,
0.0,
0.0
};
/*
int otest(), ctest();
*/
struct orbit test = {
"Test orbit ",
2446800.5,
1.8498,
49.457,
286.343,
1.523710,
0.524023,
0.093472,
53.1893,
2446800.5,
-1.52,
4.68,
0,
0.0,
0.0,
0.0
};
/* coordinates of object
*/
int objnum = 0; /* I.D. number of object */
double robject[3] = {0.0}; /* position */
/* ecliptic polar coordinates:
* longitude, latitude in radians
* radius in au
*/
double FAR obpolar[3];
/* coordinates of Earth
*/
/* Heliocentric rectangular equatorial position
* of the earth at time TDT re equinox J2000
*/
double FAR rearth[3];
/* Corresponding polar coordinates of earth:
* longitude and latitude in radians, radius in au
*/
double FAR eapolar[3];
/* Julian date of ephemeris
*/
double JD;
double TDT;
double UT;
/* flag = 0 if TDT assumed = UT,
* = 1 if input time is TDT,
* = 2 if input time is UT.
*/
int jdflag = 0;
/* correction vector, saved for display */
double dp[3];
/* display formats for printf()
*/
extern char *intfmt, *dblfmt;
/* display enable flag
*/
int prtflg = 1;
/* Tabulation parameters
*/
static double djd = 1.0;
static int ntab = 1;
struct orbit *elobject; /* pointer to orbital elements of object */
/* Main program starts here.
*/
int main()
{
int i;
double zgetdate(), gethms();
kinit();
loop:
prtflg = 1;
printf( "Enter starting date of tabulation\n" );
JD = zgetdate(); /* date */
JD += gethms(); /* time of day */
update(); /* find UT and ET */
printf( "Julian day %.7f\n", JD );
getnum( "Enter interval between tabulations in days", &djd, dblfmt );
getnum( "Number of tabulations to display", &ntab, intfmt );
if( ntab <= 0 )
ntab = 1;
loop1:
getnum( "Planet number 0-9 or 88 to read star, 99 to read orbit",
&objnum, intfmt );
switch(objnum)
{
case -1: exit(0);
case 0: elobject = 0;
printf( "\n The Sun\n" );
break;
case 1: elobject = &mercury; break;
case 2: elobject = &venus; break;
case 3: elobject = 0;
printf( "\n The Moon\n" );
break;
case 4: elobject = &mars; break;
case 5: elobject = &jupiter; break;
case 6: elobject = &saturn; break;
case 7: elobject = &uranus; break;
case 8: elobject = &neptune; break;
case 9: elobject = &pluto; break;
case 10: elobject = &test; break;
case 88:
elobject = (struct orbit *)&fstar;
i = getstar( (struct star *) elobject );
if( i == 1 )
goto loop1;
if( i == 0 )
break;
goto operr;
case 99:
elobject = &forbit;
i = getorbit( elobject );
if( i == 1 )
goto loop1;
if( i == 0 )
break;
default:
operr: printf( "Operator error.\n" );
goto loop;
}
if( elobject == (struct orbit *)&fstar )
showcname( &elobject->obname[0] );
else if( elobject )
printf( "\n %s\n", &elobject->obname[0] );
for( i=0; i<ntab; i++ )
{
/* print Julian date */
printf( "\nJD %.2f, ", JD );
update();
/* Always calculate heliocentric position of the earth */
kepler( TDT, &earth, rearth, eapolar );
switch( objnum )
{
case 0: dosun(); iter_trnsit( dosun ); break;
case 3: domoon(); iter_trnsit( domoon ); break;
case 88: dostar(); iter_trnsit( dostar ); break;
default: doplanet(); iter_trnsit( doplanet ); break;
}
printf( "\n" );
JD += djd;
}
goto loop;
#ifdef _MSC_VER
return 0;
#endif
}
syntax highlighted by Code2HTML, v. 0.9.1