/* Orbital elements and perturbations for the planet Mars
 * using formulas given by Meeus
 */
#include "planet.h"
#include "kep.h"


omars(e,J)
struct orbit *e;
double J;
{
double p;

e->epoch = J;
manoms(J);

e->M = M4;
e->a = 1.5236883;
e->ecc = ( -0.000000077*T + 0.000092064)*T + 0.09331290;
#if OFDATE
e->equinox = J;
e->i = ( 0.0000126*T - 0.0006750)*T + 1.850333;
f = ((0.00000414*T + 0.0001313)*T + 1.0697667)*T + 285.431761;
e->w = mod360(f);
f = (( -0.00000533*T - 0.0000014)*T + 0.7709917)*T + 48.786442;
e->W = mod360(f);
f = ( 0.0003107*T + 19141.69551)*T + 293.737334;
e->L = mod360(f);
#else
e->equinox = J2000;
e->i = ((-4.4e-8*T - 0.00002304)*T - 0.0081565)*T + 1.857866;
f = ((6.939e-6*T + 0.00046556)*T + 0.7387251)*T + 285.762379;
e->w = mod360(f);
f = (( -8.159e-6*T - 0.00064344)*T - 0.2941821)*T + 49.852347;
e->W = mod360(f);
f = e->M + e->w + e->W;
e->L = mod360(f);
#endif

/* long period perturbation */
f = (3.0*M5 - 8.0*M4 + 4.0*M)*DTR;
p = -0.01133*sin(f)
    -0.00933*cos(f);
e->L += p;
e->M += p;
return(0);
}




cmars(e)
struct orbit *e;
{
double f, p;
        
f = (M5 - M4 - 48.958)*DTR;
p =  0.00705*cos(f);
f = (2.0*M5 - M4 - 188.350)*DTR;
p += 0.00607*cos(f);
f = (2.0*M5 - 2.0*M4 - 191.897)*DTR;
p += 0.00445*cos(f);
f = (M - 2.0*M4 + 20.495)*DTR;
p += 0.00388*cos(f);
f = (M - M4 + 35.097)*DTR;
p += 0.00238*cos(f);
f = (2.0*M - 3.0*M4 + 158.638)*DTR;
p += 0.00204*cos(f);
f = (3.0*M4 - M2 - 57.602)*DTR;
p += 0.00177*cos(f);
f = (2.0*M - 4.0*M4 + 154.093)*DTR;
p += 0.00136*cos(f);
f = (M5 + 17.618)*DTR;
p += 0.00104*cos(f);
e->L += p*DTR;
        
f = (M5 - M4 + 41.1306)*DTR;
p = 0.000053227*cos(f);
f = (2.0*M5 - 2.0*M4 - 101.9847)*DTR;
p += 0.000050989*cos(f);
f = (2.0*M5 - M4 - 98.3292)*DTR;
p += 0.000038278*cos(f);
f = (M - M4 - 55.555)*DTR;
p += 0.000015996*cos(f);
f = (2.0*M - 3.0*M4 + 68.622)*DTR;
p += 0.000014764*cos(f);
f = (M5 - 2.0*M4 + 43.615)*DTR;
p += 0.000008966*cos(f);
f = (3.0*M5 - 2.0*M4 - 139.737)*DTR;
p += 0.000007914*cos(f);
f = (2.0*M5 - 3.0*M4 - 102.888)*DTR;
p += 0.000007004*cos(f);
f = (M - 2.0*M4 + 113.202)*DTR;
p += 0.000006620*cos(f);
f = (3.0*M5 - 3.0*M4 - 76.243)*DTR;
p += 0.000004930*cos(f);
f = (3.0*M - 5.0*M4 + 190.603)*DTR;
p += 0.000004693*cos(f);
f = (2.0*M - 4.0*M4 + 244.702)*DTR;
p += 0.000004571*cos(f);
f = (3.0*M5 - M4 - 115.828)*DTR;
p += 0.000004409*cos(f);
e->r += p;
return(0);
}