/* Follow ORBITS, based on Spillar's newplanet.c, but            */
/* with better catalog handling.                                 */
/* Robert R. Howell   Jan. 15, 1995                              */
/* I may be able to eliminate a lot of the planet constants, but */
/* I will try the least changed code first.                      */

/* This modification of Steve Moshier's code by Earl Spillar
* I hacked the code so that it could work inside the wiro
* tracking environment.  The basic changes are were to gather info
* from the wiro shared memory area.
*/

/* 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 orbital perturbations given in
* Jean Meeus, "Astronomical Formulae for Calculators", 3rd ed.,
* Willmann-Bell, Inc., 1985.
*
* 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
*/
/* Here come the WIRO additions */

#include "../track/wirotypes.h"
#include "../track/wiro.h"
#include "../track/track.h"

#include <stdio.h>
#include <stdlib.h>

void log_entry( char *comment );

struct  wiro_memory *tinfo;

double showcor();
/* 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;

/* 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 */

/* approximate motion of right ascension and declination
* of object, in radians per day
*/
double dradt = 0.0;
double ddecdt = 0.0;
/* Data structures containing orbital elements of
* objects that orbit the sun.  See kep.h for the definition.
*/
#include "kep.h"

/* Space for star description read from a disc file.
*/
struct star fstar = {
0};

/* 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 = {
0};

/* Orbits for each planet.  The indicated orbital elements are
* not actually used, since they are now calculated
* from the Meeus expansions.  Magnitude and semidiameter
* are still used.
*/
/* Programs to compute perturbations. */
#if !BandS
int omercury();
#endif
int cmercury();
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,
#if BandS
0,
#else
omercury,
#endif
cmercury,
0.0,
0.0,
0.0
};

#if !BandS
int ovenus();
#endif
int cvenus();
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,
#if BandS
0,
#else
ovenus,
#endif
cvenus,
0.0,
0.0,
0.0
};

/* Meeus purturbations are invoked for earth if BandS = 0.
* If BandS = 1, the more accurate B&S expansion is computed.
* Fixed numerical values will be used if read in from a file
* named earth.orb.
* See kfiles.c, kep.h, oearth.c, and pearth.c.
*/
#if !BandS
int oearth();
#endif
int cearth();
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,
#if BandS
0,
#else
oearth,
#endif
cearth,
0.0,
0.0,
0.0
};
extern struct orbit earth;


int omars();
int cmars();
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,
omars,
cmars,
0.0,
0.0,
0.0
};

#if !BandS
int ojupiter();
#endif
int cjupiter();
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,
#if BandS
0,
cjupiter,
#else
ojupiter,
0,
#endif
0.0,
0.0,
0.0
};

int osaturn(), csaturn();
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,
osaturn,
csaturn,
0.0,
0.0,
0.0
};

int ouranus(), curanus();
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,
ouranus,
curanus,
0.0,
0.0,
0.0
};

int oneptune(), cneptune();
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,
oneptune,
cneptune,
0.0,
0.0,
0.0
};

/* Note there are no perturbation formulas for Pluto.
* The program automatically uses the given numerical
* values since the pointers to perturbation subroutines
* are null.
* For some reason the J2000 orbit given for Pluto in the AA does
* not give the same results as the "of date" orbit.  Yet, results
* are the same for the other planets.
*/
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,
0,
0,
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,
/*
otest,
ctest,
*/
0,0,
0.0,
0.0,
0.0
};


/* coordinates of object
*/
int objnum = 0;	/* I.D. number of object */
static double robject[3] = {
0.0}; /* position */
/* ecliptic polar coordinates:
* longitude, latitude in radians
* radius in au
*/
double obpolar[3] = {
0.0};

/* coordinates of Earth
*/
/* Heliocentric rectangular equatorial position
* of the earth at time TDT re equinox J2000
*/
double rearth[3] = {
0.0};
/* Corresponding polar coordinates of earth:
* longitude and latitude in radians, radius in au
*/
double eapolar[3] = {
0.0};

/* Julian date of ephemeris
*/
double JD = 0.0;
double TDT = 0.0;
double UT = 0.0;

/* 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] = {
0.0};

/* 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;
*/

/**********************************************************/
/* Main program starts here*/
main(argc, argv)
int argc;
char *argv[];
{
	struct orbit *elobject;	/* pointer to orbital elements of object */
	int i;
	double getdate(), gethms();
	char display_name[80];	/* Name to display for object: = number then name */		
	char cmd_line_for_log[80];

	if (argc < 2) {
		printf("Usage:   orbit  object_name\n");
		printf("   or    orbit  object_short_name\n");
		exit(1);
		}

	/* Open up wiro memory */
	tinfo = get_tinfo();
	
	kinit();

	prtflg = 1;
	JD = getdate(); /* date */
	JD += gethms();	 /* time of day , leave out for our stuff */
	printf("Circumstances - ");
	update(); /* find UT and ET */
	printf( "Julian day %.7lf\n", JD );

	elobject = &forbit;
	i = read_orbit_elements( argv[1], elobject, display_name );
	
	strcpy(cmd_line_for_log, "orbit ");
	strcat(cmd_line_for_log, argv[1]);
	log_entry(cmd_line_for_log);

	/* Always calculate heliocentric position of the earth */
	kepler( TDT, &earth, rearth, eapolar );

	/* calculate heliocentric position of the object */
	kepler( TDT, elobject, robject, obpolar );
	/* apply correction factors and print apparent place */
	reduce( elobject, robject, rearth );

	strncpy( tinfo->object_name, display_name, 40 );	/* Can be up to 80 */
	

}