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


/* Calculations needed for Jupiter through Neptune:
 */
pjup()
{

nu = T/5.0 + 0.1;
P = 237.47555 + 3034.9061*T;
Q = 265.91650 + 1222.1139*T;
S = 243.51721 + 428.4677*T;
f = 5.0*Q - 2.0*P;
V = mod360(f)*DTR;
f = 2.0*P - 6.0*Q + 3.0*S;
W = mod360(f)*DTR;
f = Q - P;
ze = mod360(f)*DTR;
f = S - Q;
psi = mod360(f)*DTR;
P = mod360(P)*DTR;
Q = mod360(Q)*DTR;
S = mod360(S)*DTR;

sinV = sin(V);        
cosV = cos(V);
/* Use trigonometry to speed up sin, cos of multiple angles.
 * This reduces the relative accuracy, but we only need
 * absolute accuracy to about 8 places.
 */
sin2V = 2.0*sinV*cosV;
cos2V = cosV*cosV - sinV*sinV;
sinW = sin(W);
cosze = cos(ze);
sinze = sin(ze);

sin2ze = 2.0*sinze*cosze;
cos2ze = cosze*cosze - sinze*sinze;

sin3ze = sinze*cos2ze + cosze*sin2ze;
cos3ze = cosze*cos2ze - sinze*sin2ze;

sin4ze = sinze*cos3ze + cosze*sin3ze;
cos4ze = cosze*cos3ze - sinze*sin3ze;

sin5ze = sinze*cos4ze + cosze*sin4ze;
cos5ze = cosze*cos4ze - sinze*sin4ze;

sinQ = sin(Q);
cosQ = cos(Q);

sin2Q = 2.0*sinQ*cosQ;
cos2Q = cosQ*cosQ - sinQ*sinQ;

sin3Q = sinQ*cos2Q + cosQ*sin2Q;
cos3Q = cosQ*cos2Q - sinQ*sin2Q;
}


/* Elements and perturbations for Jupiter
 */
ojupiter(e,J)
struct orbit *e;
double J;
{

e->epoch = J;
manoms(J);
pjup();
e->M = M5;
e->a = 5.202561;
e->ecc = (( -0.0000000017*T - 0.0000004676)*T + 0.000164180)*T + 0.04833475;
#if OFDATE
e->equinox = J;
e->i = ( 0.0000039*T - 0.0056961)*T + 1.308736;
f = (( 0.00000508*T + 0.00070405)*T + 0.5994317)*T + 273.277558;
e->w = mod360(f);
f = (( -0.00000851*T + 0.00035222)*T + 1.0105300)*T + 99.443414;
e->W = mod360(f);
f = (( -0.00000165*T + 0.0003347)*T + 3036.301986)*T + 238.049257;
e->L = mod360(f);
#else
e->equinox = J2000;
e->i = ((1.27e-7*T + 2.942e-5)*T - 0.0022374)*T + 1.305288;
f = ((8.999e-6*T - 0.00021857)*T + 0.0478404)*T + 273.829584;
e->w = mod360(f);
f = (( -1.2460e-5*T + 0.00096672)*T + 0.1659357)*T + 100.287838;
e->W = mod360(f);
f = e->M + e->w + e->W;
e->L = mod360(f);
#endif
djup(e);
return(0);
}



/* Program split arbitrarily to get it through
 * Microsoft C compiler:
 */
djup(e)
struct orbit *e;
{
double p, q;

/* perturbations in mean longitude */
p = ((-0.004692*nu - 0.010281)*nu + 0.331364)*sinV
   +((0.002075*nu - 0.064436)*nu + 0.003228)*cosV
   -(( -0.000489*nu + 0.000275)*nu + 0.003083)*sin2V
   +0.002472*sinW
   +0.013619*sinze
   +0.018472*sin2ze
   +0.006717*sin3ze
   +0.002775*sin4ze;
f = (0.007275 - 0.001253*nu)*sinze
    +0.006417*sin2ze
    +0.002439*sin3ze
    -(0.033839 + 0.001125*nu)*cosze
    -0.003767*cos2ze;
p += f*sinQ;
f = -(0.035681 + 0.001208*nu)*sinze
     -0.004261*sin2ze
     +0.002178
     +(-0.006333 + 0.001161*nu)*cosze
     -0.006675*cos2ze
     -0.002664*cos3ze;
p += f*cosQ;   
f =  -0.002572*sinze
     -0.003567*sin2ze;
p += f*sin2Q;
f =  0.002094*cosze
    +0.003342*cos2ze;
p += f*cos2Q;

/* perturbations in the perihelion */
q = (0.007192 - 0.003147*nu)*sinV
   +((0.000197*nu - 0.000675)*nu - 0.020428)*cosV;
f = (0.007269 + 0.000672*nu)*sinze
    -0.004344
    +0.034036*cosze
    +0.005614*cos2ze
    +0.002964*cos3ze;
q += f*sinQ;
f =  0.037761*sinze
    +0.006158*sin2ze
    -0.006603*cosze;
q += f*cosQ;
f = -0.005356*sinze
    +0.002722*sin2ze
    +0.004483*cosze
    -0.002642*cos2ze;
q += f*sin2Q;
f =  0.004403*sinze
    -0.002536*sin2ze
    +0.005547*cosze
    -0.002689*cos2ze;
q += f*cos2Q;        
e->w += q;
#if DEBUG
printf( "pperih %.5e\n", q );
#endif
f = p - (q/(e->ecc));
e->M += f;
#if DEBUG
printf( "pmanom %.5e\n", f );
#endif
e->L += p;
#if DEBUG
printf( "plong %.5e\n", p );
#endif

/* perturbations in the eccentricity */
q = ((-.0000043*nu + .0000130)*nu + .0003606)*sinV;
q += (.0001289 - .0000580*nu)*cosV;
f = -0.0006764*sinze
    -0.0001110*sin2ze
    -0.0000224*sin3ze
    -0.0000204
    +( 0.0001284 + 0.0000116*nu)*cosze
    +0.0000188*cos2ze;
q += f*sinQ;
f = (0.0001460 + 0.0000130*nu)*sinze
    +0.0000224*sin2ze
    -0.0000817
    +0.0006074*cosze
    +0.0000992*cos2ze
    +0.0000508*cos3ze
    +0.0000230*cos4ze
    +0.0000108*cos5ze;
q += f*cosQ;
f = -(0.0000956 + 0.0000073*nu)*sinze
    +0.0000448*sin2ze
    +0.0000137*sin3ze
    +( -0.0000997 + 0.0000108*nu)*cosze
    +0.0000480*cos2ze
    +0.0000148*cos3ze;
q += f*sin2Q;
f = ( -0.0000956 +0.0000099*nu)*sinze
    +0.0000490*sin2ze
    +0.0000158*sin3ze
    +0.0000179
    +( 0.0001024 + 0.0000075*nu)*cosze
    -0.0000437*cos2ze
    -0.0000132*cos3ze;
q += f*cos2Q;     
e->ecc += q;

/* perturbations in the semimajor axis */
q = -0.000263*cosV
    +0.000205*cosze
    +0.000693*cos2ze
    +0.000312*cos3ze
    +0.000147*cos4ze;
f =  0.000299*sinze
    +0.000181*cos2ze;
q += f*sinQ;
f =  0.000204*sin2ze
    +0.000111*sin3ze
    -0.000337*cosze
    -0.000111*cos2ze;
q += f*cosQ;        
e->a += q;
}

cjupiter(e)
struct orbit *e;
{

e->plat = 0.0;
return(0);
}