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

/* Orbital elements
 */
ovenus(e,J)
struct orbit *e;
double J;
{
double p;

e->epoch = J;
manoms(J);
e->M = M2;
e->a = 0.7233316;
e->ecc = (0.000000091*T - 0.00004774)*T + 0.00682069;
#if OFDATE
e->equinox = J;
e->i = (-0.0000010*T + 0.0010058)*T + 3.393631;
f =  ( -0.0013864*T + 0.5081861)*T + 54.384186;
e->w = mod360(f);
f = ( 0.0004100*T + 0.8998500)*T + 75.779647;
e->W = mod360(f);
f = ( 0.0003097*T + 58519.21191)*T + 342.767053;
e->L = mod360(f);
#else
e->equinox = J2000;
e->i = ((1.8e-8*T - 3.250e-5)*T - 0.0007913)*T + 3.395459;
f =  ((-7.94e-7*T - 0.00114464)*T + 0.2892764)*T + 54.602827;
e->w = mod360(f);
f = ((7.69e-7*T - 0.00014010)*T - 0.2776656)*T + 76.957740;
e->W = mod360(f);
f = e->M + e->w + e->W;
e->L = mod360(f);
#endif
/* The following perturbation is applied
 * before solving Kepler's equation:
 */
f = (237.24 + 150.27*T)*DTR;
p = 0.00077 * sin(f);
e->L += p;
e->M += p;
return(0);
}


/* Perturbations added after solving Kepler's equation
 */        
cvenus(e)
struct orbit *e;
{
double p;

/* perturbation in longitude */
f = (2.0*M - 2.0*M2 - 148.225)*DTR;
p = 0.00313*cos(f);
f = (3.0*M - 3.0*M2 + 2.565)*DTR;
p += 0.00198*cos(f);
f = (M - M2 - 119.107)*DTR;
p += 0.00136*cos(f);
f = (3.0*M - 2.0*M2 - 135.912)*DTR;
p += 0.00096*cos(f);
f = (M5 - M2 - 208.087)*DTR;
p += 0.00082*cos(f);
e->L += p*DTR;
        
/* perturbation in radius vector */
f = (2.0*M - 2.0*M2 - 58.208)*DTR;
p = 0.000022501 * cos(f);
f = (3.0*M - 3.0*M2 + 92.577)*DTR;
p += 0.000019045 * cos(f);
f = (M5 - M2 - 118.090)*DTR;
p += 0.000006887 * cos(f);
f = (M - M2 - 29.110)*DTR;
p += 0.000005172 * cos(f);
f = (5.0*M - 4*M2 - 104.208)*DTR;
p += 0.000003620 * cos(f);
f = (4.0*M - 4.0*M2 + 63.513)*DTR;
p += 0.000003283 * cos(f);
f = (2.0*M5 - 2.0*M2 - 55.167)*DTR;
p += 0.000003074 * cos(f);
e->r += p;
return(0);
}