/* Elements of the orbit of the planet Saturn * using formulas given by Meeus. * * This program seems to have much larger errors than * the ones for the other planets. Error of 90 arc seconds * in R.A. and 15 arc seconds in Dec observed re 1986 tables. */ #include "planet.h" #include "kep.h" static double sinpsi = 0.0; static double cospsi = 0.0; static double sin2psi = 0.0; static double cos2psi = 0.0; static double sin3psi = 0.0; static double cos3psi = 0.0; osaturn(e,J) struct orbit *e; double J; { double p, q; e->epoch = J; manoms(J); pjup(); sinpsi = sin(psi); cospsi = cos(psi); sin2psi = 2.0*sinpsi*cospsi; cos2psi = cospsi*cospsi - sinpsi*sinpsi; sin3psi = sinpsi*cos2psi + cospsi*sin2psi; cos3psi = cospsi*cos2psi - sinpsi*sin2psi; e->M = M6; #if OFDATE /* expansions for equinox of date */ e->equinox = J; /* inclination */ e->i = (( 0.00000004*T - 0.00001549)*T - 0.0039189)*T + 2.492519; /* arg perihelion */ f = (( 0.00000992*T + 0.00097854)*T + 1.0852207)*T + 338.307800; e->w = mod360(f); /* ascending node */ f = (( -0.00000531*T - 0.00015218)*T + 0.8731951)*T + 112.790414; e->W = mod360(f); /* mean longitude */ f = (( -0.0000058*T + 0.0003245)*T + 1223.509884)*T + 266.564377; e->L = mod360(f); #else /* J2000 coordinates */ e->equinox = J2000; /* inclination */ e->i = (( 2e-9*T - 5.017e-5)*T + 0.0024449)*T + 2.486204; /* arg perihelion */ f = (( 6.177e-6*T + 0.00070747)*T + 0.8220515)*T + 338.571353; e->w = mod360(f); /* ascending node */ f = (( -1.589e-6*T - 0.00018997)*T - 0.2599254)*T + 113.923406; e->W = mod360(f); /* mean longitude */ f = e->M + e->w + e->W; e->L = mod360(f); #endif /* eccentricity */ e->ecc = (( 0.00000000074*T - 0.000000728)*T - 0.00034550)*T + 0.05589232; /* mean distance */ e->a = 9.554747; /* perturbations in the mean longitude */ p = ((0.016714*nu + 0.018150)*nu - 0.814181)*sinV + (( -0.004100*nu + 0.160906)*nu - 0.010497)*cosV + 0.007581*sin2V - 0.007986*sinW; p = p - 0.148811*sinze - 0.040786*sin2ze - 0.015208*sin3ze - 0.006339*sin4ze; f = -0.006244 + (0.008931 + 0.002728*nu)*sinze - 0.016500*sin2ze - 0.005775*sin3ze + (0.081344 + 0.003206*nu)*cosze + 0.015019*cos2ze; p += f*sinQ; f = (0.085581 + 0.002494*nu)*sinze + (0.025328 - 0.003117*nu)*cosze + 0.014394*cos2ze + 0.006319*cos3ze; p += f*cosQ; f = 0.006369*sinze + 0.009156*sin2ze + 0.007525*sin3psi; p += f*sin2Q; f = -0.005236*cosze - 0.007736*cos2ze - 0.007528*cos3psi; p += f*cos2Q; e->L += p; /* perturbations in the perihelion */ q = ((-0.001533*nu + 0.007186)*nu + 0.077108)*sinV + ((-0.000536*nu - 0.014766)*nu + 0.045803)*cosV - 0.007075*sinze; f = -0.075825*sinze - 0.024839*sin2ze - 0.008631*sin3ze; q += f*sinQ; f = -0.072586 - 0.150383*cosze + 0.026897*cos2ze + 0.010053*cos3ze; q += f*cosQ; f = -(0.013597 +0.001719*nu)*sinze + (-0.007742 + 0.001517*nu)*cosze + (0.013586 - 0.001375*nu)*cos2ze; q += f*sin2Q; f = (-0.013667 + 0.001239*nu)*sinze + 0.011981*sin2ze + (0.014861 + 0.001136*nu)*cosze - (0.013064 + 0.001628*nu)*cos2ze; q += f*cos2Q; e->w += q; f = p - q/(e->ecc); e->M += f; esat(e); asat(e); return(0); } esat(e) struct orbit *e; { double sin4Q, cos4Q; double q; /* perturbations in the eccentricity */ q = ((0.0000091*nu + 0.0002548)*nu - 0.0007927)*sinV; q += (( -0.0000253*nu + 0.0001226)*nu + 0.0013381)*cosV; q += ( -0.0000121*nu + 0.0000248)*sin2V; q -= (0.0000091*nu + 0.0000305)*cos2V; q += 0.0000412*sin2ze; f = 0.0012415 + ( 0.0000390 - 0.0000617*nu)*sinze + ( 0.0000165 - 0.0000204*nu)*sin2ze; f = f + 0.0026599*cosze - 0.0004687*cos2ze - 0.0001870*cos3ze - 0.0000821*cos4ze - 0.0000377*cos5ze; f = f + 0.0000497*cos2psi; q += f*sinQ; f = 0.0000163 - 0.0000611*nu - 0.0012696*sinze - 0.0004200*sin2ze - 0.0001503*sin3ze - 0.0000619*sin4ze - 0.0000268*sin5ze; f = f - ( 0.0000282 + 0.0001306*nu)*cosze + (-0.0000086 + 0.0000230*nu)*cos2ze; f = f + 0.0000461*sin2psi; q += f*cosQ; f = -0.0000350 + (0.0002211 - 0.0000286*nu)*sinze - 0.0002208*sin2ze - 0.0000568*sin3ze - 0.0000346*sin4ze; f = f - (0.0002780 + 0.0000222*nu)*cosze + (0.0002022 + 0.0000263*nu)*cos2ze + 0.0000248*cos3ze + 0.0000242*sin3psi + 0.0000467*cos3psi; q += f*sin2Q; f = -0.0000490 - ( 0.0002842 + 0.0000279*nu)*sinze + ( 0.0000128 + 0.0000226*nu)*sin2ze; f = f + 0.0000224*sin3ze + (-0.0001594 + 0.0000282*nu)*cosze + ( 0.0002162 - 0.0000207*nu)*cos2ze; f = f + 0.0000561*cos3ze + 0.0000343*cos4ze + 0.0000469*sin3psi - 0.0000242*cos3psi; q += f*cos2Q; f = -0.0000205*sinze + 0.0000262*sin3ze; q += f*sin3Q; f = 0.0000208*cosze - 0.0000271*cos3ze; q += f*cos3Q; sin4Q = 2.0*sin2Q*cos2Q; cos4Q = cos2Q*cos2Q - sin2Q*sin2Q; q = q - 0.0000382*cos3ze*sin4Q - 0.0000376*sin3ze*cos4Q; e->ecc += q; } asat(e) struct orbit *e; { double q; /* perturbations in the semimajor axis */ q = 0.000572*nu*sinV + 0.002933*cosV + 0.033629*cosze - 0.003081*cos2ze - 0.001423*cos3ze - 0.000671*cos4ze - 0.000320*cos5ze; f = 0.001098 - 0.002812*sinze + 0.000688*sin2ze - 0.000393*sin3ze - 0.000228*sin4ze + 0.002138*cosze - 0.000999*cos2ze - 0.000642*cos3ze - 0.000325*cos4ze; q += f*sinQ; f = -0.000890 + 0.002206*sinze - 0.001590*sin2ze - 0.000647*sin3ze - 0.000344*sin4ze + 0.002885*cosze + ( 0.002172 + 0.000102*nu)*cos2ze + 0.000296*cos3ze; q += f*cosQ; f = -0.000267*sin2ze - 0.000778*cosze + 0.000495*cos2ze + 0.000250*cos3ze; q += f*sin2Q; f = -0.000856*sinze + 0.000441*sin2ze + 0.000296*cos2ze + 0.000211*cos3ze; q += f*cos2Q; f = -0.000427*sinze + 0.000398*sin3ze; q += f*sin3Q; f = 0.000344*cosze - 0.000427*cos3ze; q += f*cos3Q; e->a += q; } /* Corrections to the position to be applied * after solving Kepler's equation */ csaturn(e) struct orbit *e; { double p; /* perturbations to the heliocentric latitude */ f = 0.000747*sinQ +0.001069*cosQ; p = f*cosze; f = 0.002108*sin2ze +0.001261*cos2ze; p = p + f*sin2Q; f = 0.001236*sin2ze -0.002075*cos2ze; p = p + f*cos2Q; e->plat = p*DTR; return(0); }