/* Orbital elements and perturbations for the earth
 * Meeus, chapter 18
 */

#include "kep.h"

extern double M, T;


oearth(e,J)
struct orbit *e;
double J;
{
double A, B, C, D, E, H, f;

e->epoch = J;
T = (J - J1900)/36525.0; /* centuries from 1900.0 */

/* mean anomaly of the earth (and sun)
 */
f = (( -3.3e-6*T - 0.000150)*T + 35999.04975)*T + 358.47583;
f = mod360(f);
e->M = f;
M = f;	/* save in global place */

/* eccentricity of the earth's orbit */
e->ecc = (-0.000000126*T - 0.0000418)*T + 0.01675104;

/* semimajor axis */
e->a = 1.0000002;

#if OFDATE
/* epoch of equinox */
e->equinox = J;
/* inclination */
e->i = 0.0;
/* ascending node */
e->W = 0.0;
/* geometric mean longitude of the earth */
f = (0.0003025*T + 36000.76892)*T + 99.69668; /*of sun: 279.69668*/
e->L = mod360(f);
/* Compute the argument of the perihelion from
 * Mean anomaly = Mean longitude - Longitude of perihelion
 */
f = e->L - e->M;	/* arg perihelion */
e->w = mod360(f);
#else
e->equinox = J2000;
e->i = ((1.4e-8*T - 9.27e-6)*T + 0.0130855)*T - 0.0130762;
f = ((3.333e-6*T + 0.00013610)*T + 0.5647920)*T + 287.511505;
e->w = mod360(f);
f = (( -2.8e-8*T + 7.94e-6)*T - 0.2416582)*T + 175.105679;
e->W = mod360(f);
f = e->w + e->M + e->W;
e->L = mod360(f);
#endif
return(0);
}





/* perturbations of the earth's orbit
 * added in after solving Kepler's equation
 */

cearth(e)
struct orbit *e;
{
double A, B, C, D, E, H, f;
double sin(), cos();

/* perturbations due to Venus: */
A = (22518.7541*T + 153.23)*DTR;
B = (45037.5082*T + 216.57)*DTR;
/* due to Jupiter: */
C = (32964.3577*T + 312.69)*DTR;
H = (65928.7155*T + 353.40)*DTR;
/* due to the moon: */
D = ((-0.00144*T + 445267.1142)*T + 350.74)*DTR;
/* long period perturbation: */
E = (20.20*T + 231.19)*DTR;

/* corrections to earth's longitude */
f =   0.00134*cos(A)
    + 0.00154*cos(B)
    + 0.00200*cos(C)
    + 0.00179*sin(D)
    + 0.00178*sin(E);
e->L += f*DTR;

/* corrections to earth's radius vector */
f =   0.00000543*sin(A)
    + 0.00001575*sin(B)
    + 0.00001627*sin(C)
    + 0.00003076*cos(D)
    + 0.00000927*sin(H);
e->r += f;
return(0);
}