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

static double sinG = 0.0;
static double cosG = 0.0;

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

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

f = (-0.000070*T + 218.46134)*T + 37.73063;
M8 = mod360(f);
e->M = M8;
e->a = 30.10957;
e->ecc = ( -0.000000002*T + 0.000006330)*T + 0.00899704;
#if OFDATE
e->equinox = J;
e->i = ( -0.0000091*T - 0.0095436)*T + 1.779242;
f = (( 0.000004113*T + 0.00014095)*T + 0.3256394)*T + 276.045975;
e->w = mod360(f);
f = (( -0.000004718*T + 0.00024987)*T + 1.0989350)*T + 130.681389;
e->W = mod360(f);
f = (( -0.00000060*T + 0.0003205)*T + 219.885914)*T + 84.457994;
e->L = mod360(f);
#else
e->equinox = J2000;
e->i = ((1.8e-8*T - 2.27e-6)*T - 0.0000144)*T + 1.769715;
f = ((2.226e-6*T + 3.849e-5)*T + 0.0368127)*T + 276.335328;
e->w = mod360(f);
f = (( -2.858e-6*T + 4.428e-5)*T - 0.0084187)*T + 131.788486;
e->W = mod360(f);
f = e->M + e->w + e->W;
e->L = mod360(f);
#endif

G = 83.76922 + 218.4901*T;
G = mod360(G)*DTR;
H = modtp(2.0*G - S);
ze = modtp(G - P);
eta = modtp(G - Q);
th = modtp(G - S);
sinth = sin(th);
costh = cos(th);
sin2th = 2.0*sinth*costh;
cos2th = costh*costh - sinth*sinth;
sin1H = sin(H);
cos1H = cos(H);
sin2H = 2.0*sin1H*cos1H;
cos2H = cos1H*cos1H - sin1H*sin1H;
sinG = sin(G);
cosG = cos(G);

/* perturbations in mean longitude */        
p = (-0.589833 + 0.001089*nu)*sin1H;
p = p +(-0.056094 + 0.004658*nu)*cos1H;
p = p -0.024286*sin2H;
e->L += p;

/* perturbations in the perihelion */
q = 0.024039*sin1H
   +0.006206*sin2H
   -0.025303*cos1H
   -0.005992*cos2H;
e->w += q;
e->M = M8 + p - q/(e->ecc);       

/* perturbations in eccentricity */
q = .0004389*sin1H
   +.0001129*sin2H
   +.0004262*cos1H
   +.0001089*cos2H;
e->ecc += q;

/* correction to the semimajor axis */
q = -0.000817*sin1H
    +0.008189*cos1H
    +0.000781*cos2H;
e->a += q;        
return(0);
}



cneptune(e)
struct orbit *e;
{
double q;

/* correction to the true longitude */
q = -0.009556*sin(ze)
    -0.005178*sin(eta)
    +0.002572*sin2th;
q = q -0.002972*cos2th*sinG;
q = q -0.002833*sin2th*cosG;
e->L += q*DTR;

/* correction to the heliocentric latitude */
q = 0.000336*cos2th*sinG;
q = q + 0.000364*sin2th*cosG;
e->plat = q * DTR;

/* correction to the radius vector */
q = -0.040596
    +0.004992*cos(ze)
    +0.002744*cos(eta)
    +0.002044*costh
    +0.001051*cos2th;
e->r += q;
return(0);
}