#include <cstdio>
#include <cstdlib>
using namespace std;
#include "body.h"
#include "xpUtil.h"
/*
The Galilean satellite ephemerides E5 are described in Lieske,
Astron. Astrophys. Suppl. Ser., 129, 205-217 (1998)
*/
static void
computeArguments(const double t,
double &l1, double &l2, double &l3, double &l4,
double &om1, double &om2, double &om3, double &om4,
double &psi, double &Gp, double &G)
{
// mean longitudes
l1 = (106.077187 + 203.48895579033 * t) * deg_to_rad;
l2 = (175.731615 + 101.37472473479 * t) * deg_to_rad;
l3 = (120.558829 + 50.31760920702 * t) * deg_to_rad;
l4 = ( 84.444587 + 21.57107117668 * t) * deg_to_rad;
// proper nodes
om1 = (312.334566 - 0.13279385940 * t) * deg_to_rad;
om2 = (100.441116 - 0.03263063731 * t) * deg_to_rad;
om3 = (119.194241 - 0.00717703155 * t) * deg_to_rad;
om4 = (322.618633 - 0.00175933880 * t) * deg_to_rad;
// longitude of origin of coordingate (Jupiter's pole)
psi = (316.518203 - 2.08362E-06 * t) * deg_to_rad;
// mean anomaly of Saturn
Gp = (31.978528 + 0.03345973390 * t) * deg_to_rad;
// mean anomaly of Jupiter
G = (30.237557 + 0.08309257010 * t) * deg_to_rad;
}
void
jupsat(const double jd, const body b,
double &X, double &Y, double &Z)
{
const double t = jd - 2443000.5;
// mean longitudes
double l1, l2, l3, l4;
// proper nodes
double om1, om2, om3, om4;
// longitude of origin of coordinates (Jupiter's pole)
double psi;
// mean anomaly of Saturn
double Gp;
// mean anomaly of Jupiter
double G;
computeArguments(t, l1, l2, l3, l4, om1, om2, om3, om4, psi, Gp, G);
// free libration
const double phi = (199.676608 + 0.17379190461 * t) * deg_to_rad;
// periapse longitudes
const double pi1 = ( 97.088086 + 0.16138586144 * t) * deg_to_rad;
const double pi2 = (154.866335 + 0.04726306609 * t) * deg_to_rad;
const double pi3 = (188.184037 + 0.00712733949 * t) * deg_to_rad;
const double pi4 = (335.286807 + 0.00183999637 * t) * deg_to_rad;
// longitude of perihelion of jupiter
const double PIj = 13.469942 * deg_to_rad;
// phase angles
const double phi1 = 188.374346 * deg_to_rad;
const double phi2 = 52.224824 * deg_to_rad;
const double phi3 = 257.184000 * deg_to_rad;
const double phi4 = 149.152605 * deg_to_rad;
// semimajor axes, in AU
const double axis1 = 2.819353E-3;
const double axis2 = 4.485883E-3;
const double axis3 = 7.155366E-3;
const double axis4 = 12.585464E-3;
// common factors
const double PIG2 = (PIj + G) * 2;
double xi = 0;
double upsilon = 0;
double zeta = 0;
double radius, lon;
switch (b)
{
case IO:
{
lon = l1;
const int NXI = 10;
double XI[NXI][2] = {
{ 170, l1 - l2 },
{ 106, l1 - l3 },
{ -2, l1 - pi1 },
{ -2, l1 - pi2 },
{ -387, l1 - pi3 },
{ -214, l1 - pi4 },
{ -66, l1 + pi3 - PIG2 },
{ -41339, 2*(l1 - l2) },
{ 3, 2*(l1 - l3) },
{ -131, 4*(l1-l2) }
};
for (int i = 0; i < NXI; i++)
xi += XI[i][0] * cos(XI[i][1]);
xi *= 1e-7;
radius = axis1 * (1 + xi);
const int NU = 41;
double U[NU][2] = {
{ -26, 2 * psi - PIG2 },
{ -553, 2*(psi - PIj) },
{ -240, om3 + psi - PIG2 },
{ 92, psi - om2 },
{ -72, psi - om3 },
{ -49, psi - om4 },
{ -325, G },
{ 65, 2*G },
{ -33, 5*Gp - 2*G + phi2 },
{ -27, om3 - om4 },
{ 145, om2 - om3 },
{ 30, om2 - om4 },
{ -38, pi4 - PIj },
{ -6071, pi3 - pi4 },
{ 282, pi2 - pi3 },
{ 156, pi2 - pi4 },
{ -38, pi1 - pi3 },
{ -25, pi1 - pi4 },
{ -27, pi1 + pi4 - PIG2 },
{ -1176, pi1 + pi3 - PIG2 },
{ 1288, phi },
{ 39, 3*l3 - 7*l4 + 4*pi4 },
{ -32, 3*l3 - 7*l4 + pi3 + 3*pi4 },
{ -1162, l1 - 2*l2 + pi4 },
{ -1887, l1 - 2*l2 + pi3 },
{ -1244, l1 - 2*l2 + pi2 },
{ 38, l1 - 2*l2 + pi1 },
{ -617, l1 - l2 },
{ -270, l1 - l3 },
{ -26, l1 - l4 },
{ 4, l1 - pi1 },
{ 5, l1 - pi2 },
{ 776, l1 - pi3 },
{ 462, l1 - pi4 },
{ 149, l1 + pi3 - PIG2 },
{ 21, 2*l1 - 4*l2 + om2 + om3 },
{ -200, 2*l1 - 4*l2 + 2*om2 },
{ 82483, 2*(l1 - l2) },
{ -35, 2*(l1 - l3) },
{ -3, 3*l1 - 4*l2 + pi3 },
{ 276, 4*(l1 - l2) }
};
for (int i = 0; i < NU; i++)
upsilon += U[i][0] * sin(U[i][1]);
upsilon *= 1e-7;
// now use the "time completed" series
const double n = 203.48895579033 * deg_to_rad;
computeArguments(t + upsilon/n, l1, l2, l3, l4,
om1, om2, om3, om4, psi, Gp, G);
const int NZ = 7;
double Z[NZ][2] = {
{ 46, l1 + psi - 2*PIj - 2*G},
{ 6393, l1 - om1 },
{ 1825, l1 - om2 },
{ 329, l1 - om3 },
{ 93, l1 - om4 },
{ -311, l1 - psi },
{ 75, 3*l1 - 4*l2 + om2 }
};
for (int i = 0; i < NZ; i++)
zeta += Z[i][0] * sin(Z[i][1]);
zeta *= 1e-7;
}
break;
case EUROPA:
{
lon = l2;
const int NXI = 24;
double XI[NXI][2] = {
{ -18, om2 - om3 },
{ -27, 2*l3 - PIG2 },
{ 553, l2 - l3 },
{ 45, l2 - l4 },
{ -102, l2 - pi1 },
{ -1442, l2 - pi2 },
{ -3116, l2 - pi3 },
{ -1744, l2 - pi4 },
{ -15, l2 - PIj - G },
{ -64, 2*(l2 - l4) },
{ 164, 2*(l2 - om2) },
{ 18, 2*l2 - om2 - om3 },
{ -54, 5*(l2 - l3) },
{ -30, l1 - 2*l2 + pi4 },
{ -67, l1 - 2*l2 + pi3 },
{ 93848, l1 - l2 },
{ 48, l1 - 2*l3 + pi4 },
{ 107, l1 - 2*l3 + pi3 },
{ -19, l1 - 2*l3 + pi2 },
{ 523, l1 - l3 },
{ 30, l1 - pi3 },
{ -290, 2*(l1 - l2) },
{ -91, 2*(l1 - l3) },
{ 22, 4*(l1 - l2) }
};
for (int i = 0; i < NXI; i++)
xi += XI[i][0] * cos(XI[i][1]);
xi *= 1e-7;
radius = axis2 * (1 + xi);
const int NU = 66;
double U[NU][2] = {
{ 98, 2*psi - PIG2 },
{ -1353, 2*(psi - PIj) },
{ 551, psi + om3 - PIG2 },
{ 26, psi + om2 - PIG2 },
{ 31, psi - om2 },
{ 255, psi - om3 },
{ 218, psi - om4 },
{ -1845, G},
{ -253, 2*G },
{ 18, 2*(Gp - G) + phi4 },
{ 19, 2*Gp - G + phi1 },
{ -15, 5*Gp - 3*G + phi1 },
{ -150, 5*G - 2*G + phi2 },
{ 102, om3 - om4 },
{ 56, om2 - om3 },
{ 72, pi4 - PIj },
{ 2259, pi3 - pi4 },
{ -24, pi3 - pi4 + om3 - om4 },
{ -23, pi2 - pi3 },
{ -36, pi2 - pi4 },
{ -31, pi1 - pi2 },
{ 4, pi1 - pi3 },
{ 111, pi1 - pi4 },
{ -354, pi1 + pi3 - PIG2 },
{ -3103, phi },
{ 55, 2*l3 - PIG2 },
{ -111, 3*l3 - 7*l4 + 4*pi4 },
{ 91, 3*l3 - 7*l4 + pi3 + 3*pi4 },
{ -25, 3*l3 - 7*l4 + 2*pi3 + 2*pi4 },
{ -1994, l2 - l3 },
{ -137, l2 - l4 },
{ 1, l2 - pi1 },
{ 2886, l2 - pi2 },
{ 6250, l2 - pi3 },
{ 3463, l2 - pi4 },
{ 30, l2 - PIj - G },
{ -18, 2*l2 - 3*l3 + pi4 },
{ -39, 2*l2 - 3*l3 + pi3 },
{ 98, 2*(l2 - l4) },
{ -164, 2*(l2 - om2) },
{ -18, 2*l2 - om2 - om3 },
{ 72, 5*(l2 - l3) },
{ 30, l1 - 2*l2 - pi3 + PIG2 },
{ 4180, l1 - 2*l2 + pi4, },
{ 7428, l1 - 2*l2 + pi3, },
{ -2329, l1 - 2*l2 + pi2, },
{ -19, l1 - 2*l2 + pi1, },
{ -185835, l1 - l2 },
{ -110, l1 - 2*l3 + pi4, },
{ -200, l1 - 2*l3 + pi3, },
{ 39, l1 - 2*l3 + pi2, },
{ -16, l1 - 2*l3 + pi1, },
{ -803, l1 - l3 },
{ -19, l1 - pi2 },
{ -75, l1 - pi3 },
{ -31, l1 - pi4 },
{ -9, 2*l1 - 4*l2 + om3 + psi },
{ 4, 2*l1 - 4*l2 + 2*om3 },
{ -14, 2*l1 - 4*l2 + om2 + om3 },
{ 150, 2*l1 - 4*l2 + 2*om2 },
{ -11, 2*l1 - 4*l2 + PIG2 },
{ -9, 2*l1 - 4*l2 + pi3 + pi4, },
{ -8, 2*l1 - 4*l2 + 2*pi3 },
{ 915, 2*(l1 - l2) },
{ 96, 2*(l1 - l3) },
{ -18, 4*(l1 - l2) }
};
for (int i = 0; i < NU; i++)
upsilon += U[i][0] * sin(U[i][1]);
upsilon *= 1e-7;
// now use the "time completed" series
const double n = 101.37472473479 * deg_to_rad;
computeArguments(t + upsilon/n, l1, l2, l3, l4,
om1, om2, om3, om4, psi, Gp, G);
const int NZ = 11;
double Z[NZ][2] = {
{ 17, l2 + psi - 2*(PIj - G) - G },
{ 143, l2 + psi - 2*(PIj - G) },
{ -144, l2 - om1 },
{ 81004, l2 - om2 },
{ 4512, l2 - om3 },
{ 1160, l2 - om4 },
{ -19, l2 - psi - G },
{ -3284, l2 - psi },
{ 35, l2 - psi + G },
{ -28, l1 - 2*l3 + om3 },
{ 272, l1 - 2*l3 + om2 }
};
for (int i = 0; i < NZ; i++)
zeta += Z[i][0] * sin(Z[i][1]);
zeta *= 1e-7;
}
break;
case GANYMEDE:
{
lon = l3;
const int NXI = 31;
double XI[NXI][2] = {
{ 24, psi - om3 },
{ -9, om3 - om4 },
{ 10, pi3 - pi4 },
{ 294, l3 - l4 },
{ 18, l3 - pi2 },
{ -14388, l3 - pi3 },
{ -7919, l3 - pi4 },
{ -23, l3 - PIj - G },
{ -20, l3 + pi4 - PIG2 },
{ -51, l3 + pi3 - PIG2 },
{ 39, 2*l3 - 3*l4 + pi4 },
{ -1761, 2*(l3 - l4) },
{ -11, 2*(l3 - pi3) },
{ -10, 2*(l3 - pi3 - pi4) },
{ -27, 2*l3 - PIG2 },
{ 24, 2*(l3 - om3) },
{ 9, 2 * l3 - om3 - om4 },
{ -24, 2 * l3 - om3 - psi },
{ -16, 3*l3 - 4*l4 + pi4 },
{ -156, 3*(l3 - l4) },
{ -42, 4*(l3 - l4) },
{ -11, 5*(l3 - l4) },
{ 6342, l2 - l3 },
{ 9, l2 - pi3 },
{ 39, 2*l2 - 3*l3 + pi4 },
{ 70, 2*l2 - 3*l3 + pi3 },
{ 10, l1 - 2*l2 + pi4 },
{ 20, l1 - 2*l2 + pi3 },
{ -153, l1 - l2 },
{ 156, l1 - l3 },
{ 11, 2*(l1 - l2) }
};
for (int i = 0; i < NXI; i++)
xi += XI[i][0] * cos(XI[i][1]);
xi *= 1e-7;
radius = axis3 * (1 + xi);
const int NU = 75;
double U[NU][2] = {
{ 10, psi - pi3 + pi4 - om3 },
{ 28, 2*psi - PIG2 },
{ -1770, 2*(psi - PIj) },
{ -48, psi + om3 - PIG2 },
{ 14, psi - om2 },
{ 411, psi - om3 },
{ 345, psi - om4 },
{ -2338, G },
{ -66, 2*G },
{ 10, Gp - G + phi3 },
{ 22, 2*(Gp - G) + phi4 },
{ 26, 2*Gp - G + phi1 },
{ 11, 3*Gp - 2*G + phi2 + phi3 },
{ 9, 3*Gp - G + phi1 - phi2 },
{ -19, 5*Gp - 3*G + phi1 },
{ -208, 5*Gp - 2*G + phi2 },
{ 159, om3 - om4 },
{ 21, om2 - om3 },
{ 121, pi4 - PIj },
{ 6604, pi3 - pi4 },
{ -65, pi3 - pi4 + om3 - om4 },
{ -88, pi2 - pi3 },
{ -72, pi2 - pi4 },
{ -26, pi1 - pi3 },
{ -9, pi1 - pi4 },
{ 16, pi1 + pi4 - PIG2 },
{ 125, pi1 + pi3 - PIG2 },
{ 307, phi },
{ -10, l4 - pi4 },
{ -100, l3 - 2*l4 + pi4 },
{ 83, l3 - 2*l4 + pi3 },
{ -944, l3 - l4 },
{ -37, l3 - pi2 },
{ 28780, l3 - pi3 },
{ 15849, l3 - pi4 },
{ 7, l3 - pi4 + om3 - om4 },
{ 46, l3 - PIj - G },
{ 51, l3 + pi4 - PIG2 },
{ 11, l3 + pi3 - PIG2 - G },
{ 97, l3 + pi3 - PIG2 },
{ 1, l3 + pi1 - PIG2 },
{ -101, 2*l3 - 3*l4 + pi4 },
{ 13, 2*l3 - 3*l4 + pi3 },
{ 3222, 2*(l3 - l4) },
{ 29, 2*(l3 - pi3) },
{ 25, 2*l3 - pi3 - pi4 },
{ 37, 2*l3 - PIG2 },
{ -24, 2*(l3 - om3) },
{ -9, 2*l3 - om3 - om4 },
{ 24, 2*l3 - om3 - psi },
{ -174, 3*l3 - 7*l4 + 4*pi4 },
{ 140, 3*l3 - 7*l4 + pi3 + 3*pi4 },
{ -55, 3*l3 - 7*l4 + 2*pi3 + 2*pi4 },
{ 27, 3*l3 - 4*l4 + pi4 },
{ 227, 3*(l3 - l4) },
{ 53, 4*(l3 - l4) },
{ 13, 5*(l3 - l4) },
{ 42, l2 - 3*l3 + 2*l4 },
{ -12055, l2 - l3 },
{ -24, l2 - pi3 },
{ -10, l2 - pi4 },
{ -79, 2*l2 - 3*l3 + pi4 },
{ -131, 2*l2 - 3*l3 + pi3 },
{ -665, l1 - 2*l2 + pi4 },
{ -1228, l1 - 2*l2 + pi3 },
{ 1082, l1 - 2*l2 + pi2 },
{ 90, l1 - 2*l2 + pi1 },
{ 190, l1 - l2 },
{ 218, l1 - l3 },
{ 2, 2*l1 - 4*l2 + om3 + psi },
{ -4, 2*l1 - 4*l2 + 2*om3 },
{ 3, 2*l1 - 4*l2 + 2*om2 },
{ 2, 2*l1 - 4*l2 + pi3 + pi4 },
{ 2, 2*l1 - 4*l2 + 2*pi3 },
{ -13, 2*(l1 - l2) }
};
for (int i = 0; i < NU; i++)
upsilon += U[i][0] * sin(U[i][1]);
upsilon *= 1e-7;
// now use the "time completed" series
const double n = 50.31760920702 * deg_to_rad;
computeArguments(t + upsilon/n, l1, l2, l3, l4,
om1, om2, om3, om4, psi, Gp, G);
const int NZ = 13;
double Z[NZ][2] = {
{ 37, l2 + psi - 2*(PIj - G) - G },
{ 321, l2 + psi - 2*(PIj - G) },
{ -15, l2 + psi - 2*PIj - G },
{ -45, l3 - 2*PIj + psi },
{ -2797, l3 - om2 },
{ 32402, l3 - om3 },
{ 6847, l3 - om4 },
{ -45, l3 - psi - G },
{ -16911, l3 - psi },
{ 51, l3 - psi + G },
{ 10, 2*l2 - 3*l3 + psi },
{ -21, 2*l2 - 3*l3 + om3 },
{ 30, 2*l2 - 3*l3 + om2 }
};
for (int i = 0; i < NZ; i++)
zeta += Z[i][0] * sin(Z[i][1]);
zeta *= 1e-7;
}
break;
case CALLISTO:
{
lon = l4;
const int NXI = 43;
double XI[NXI][2] = {
{ -19, psi - om3 },
{ 167, psi - om4 },
{ 11, G },
{ 12, om3 - om4 },
{ -13, pi3 - pi4 },
{ 1621, l4 - pi3 },
{ -24, l4 - pi4 + 2*(psi - PIj) },
{ -17, l4 - pi4 - G },
{ -73546, l4 - pi4 },
{ 15, l4 - pi4 + G },
{ 30, l4 - pi4 + 2*(PIj - psi) },
{ -5, l4 - PIj + 2*G },
{ -89, l4 - PIj - G },
{ 182, l4 - PIj },
{ -6, l4 + pi4 - 2*PIj - 4*G },
{ -62, l4 + pi4 - 2*PIj - 3*G },
{ -543, l4 + pi4 - 2*PIj - 2*G },
{ 27, l4 + pi4 - 2*PIj - G },
{ 6, l4 + pi4 - 2*PIj },
{ 6, l4 + pi4 - om4 - psi },
{ -9, l4 + pi3 - 2*pi4 },
{ 14, l4 + pi3 - PIG2 },
{ 13, 2*l4 - pi3 - pi4 },
{ -271, 2*(l4 - pi4) },
{ -25, 2*l4 - PIG2 - G },
{ -155, 2*l4 - PIG2 },
{ -12, 2*l4 - om3 - om4 },
{ 19, 2*l4 - om3 - psi },
{ 48, 2*(l4 - om4) },
{ -167, 2*l4 - om4 - psi },
{ 142, 2*(l4 - psi) },
{ -22, l3 - 2*l4 + pi4 },
{ 20, l3 - 2*l4 + pi3 },
{ 974, l3 - l4 },
{ 24, 2*l3 - 3*l4 + pi4 },
{ 177, 2*(l3 - l4) },
{ 4, 3*l3 - 4*l4 + pi4 },
{ 42, 3*(l3 - l4) },
{ 14, 4*(l3 - l4) },
{ 5, 5*(l3 - l4) },
{ -8, l2 - 3*l3 + 2*l4 },
{ 92, l2 - l4 },
{ 105, l1 - l4 }
};
for (int i = 0; i < NXI; i++)
xi += XI[i][0] * cos(XI[i][1]);
xi *= 1e-7;
radius = axis4 * (1 + xi);
const int NU = 86;
double U[NU][2] = {
{ 8, 2*psi - pi3 - pi4 },
{ -9, psi - pi3 - pi4 + om4 },
{ 27, psi - pi3 + pi4 - om4 },
{ -409, 2*(psi - pi4) },
{ 310, psi - 2*pi4 + om4 },
{ -19, psi - 2*pi4 + om3 },
{ 8, 2*psi - pi4 - PIj },
{ -5, psi - pi4 - PIj + om4 },
{ 63, psi - pi4 + PIj - om4 },
{ 8, 2*psi - PIG2 - G },
{ 73, 2*psi - PIG2 },
{ -5768, 2*(psi - PIj) },
{ 16, psi + om4 - PIG2 },
{ -97, psi - om3 },
{ 152, 2*(psi - om4) },
{ 2070, psi - om4 },
{ -5604, G },
{ -204, 2*G },
{ -10, 3*G },
{ 24, Gp - G + phi3 },
{ 11, Gp + phi1 - 2*phi2 },
{ 52, 2*(Gp - G) + phi4 },
{ 61, 2*Gp - G + phi1 },
{ 25, 3*Gp - 2*G + phi2 + phi3 },
{ 21, 3*Gp - G + phi1 - phi2 },
{ -45, 5*Gp - 3*G + phi1 },
{ -495, 5*Gp - 3*G + phi2 },
{ -44, om3 - om4 },
{ 5, pi4 - PIj - G },
{ 234, pi4 - PIj },
{ 11, 2*pi4 - PIG2 },
{ -10, 2*pi4 - om3 - om4 },
{ 68, 2*(pi4 - om4) },
{ -13, pi3 - pi4 - om4 + psi },
{ -5988, pi3 - pi4 },
{ -47, pi3 - pi4 + om3 - om4 },
{ -3249, l4 - pi3 },
{ 48, l4 - pi4 + 2*(psi - PIj) },
{ 10, l4 - pi4 - om4 + psi },
{ 33, l4 - pi4 - G },
{ 147108, l4 - pi4 },
{ -31, l4 - pi4 + G },
{ -6, l4 - pi4 + om4 - psi },
{ -61, l4 - pi4 + 2*(PIj - psi) },
{ 10, l4 - PIj - 2*G },
{ 178, l4 - PIj - G },
{ -363, l4 - PIj },
{ 5, l4 + pi4 - 2*PIj - 5*Gp + 2*G - phi1 },
{ 12, l4 + pi4 - 2*PIj - 4*G },
{ 124, l4 + pi4 - 2*PIj - 3*G },
{ 1088, l4 + pi4 - 2*PIj - 2*G },
{ -55, l4 + pi4 - 2*PIj - G },
{ -12, l4 + pi4 - 2*PIj },
{ -13, l4 + pi4 - om4 - psi },
{ 6, l4 + pi4 - 2*psi },
{ 17, l4 + pi3 - 2*pi4 },
{ -28, l4 + pi3 - PIG2 },
{ -33, 2*l4 - pi3 - pi4 },
{ 676, 2*(l4 - pi4) },
{ 36, 2*(l4 - PIj - G) - G },
{ 218, 2*(l4 - PIj - G) },
{ -5, 2*(l4 - PIj) - G },
{ 12, 2*l4 - om3 - om4 },
{ -19, 2*l4 - om3 - psi },
{ -48, 2*(l4 - om4) },
{ 167, 2*l4 - om4 - psi },
{ -142, 2*(l4 - psi) },
{ 148, l3 - 2*l4 + pi4 },
{ -94, l3 - 2*l4 + pi3 },
{ -390, l3 - l4 },
{ 9, 2*l3 - 4*l4 + 2*pi4 },
{ -37, 2*l3 - 3*l4 + pi4 },
{ 6, 2*l3 - 3*l4 + pi3 },
{ -195, 2*(l3 - l4) },
{ 6, 3*l3 - 7*l4 + 2*pi4 + om4 + psi },
{ 187, 3*l3 - 7*l4 + 4*pi4 },
{ -149, 3*l3 - 7*l4 + pi3 + 3*pi4 },
{ 51, 3*l3 - 7*l4 + 2*(pi3 + pi4) },
{ -10, 3*l3 - 7*l4 + 3*pi3 + pi4 },
{ 6, 3*(l3 - 2*l4 + pi4) },
{ -8, 3*l3 - 4*l4 + pi4 },
{ -41, 3*(l3 - l4) },
{ -13, 4*(l3 - l4) },
{ -44, l2 - 3*l3 + 2*l4 },
{ 89, l2 - l4 },
{ 106, l1 - l4 }
};
for (int i = 0; i < NU; i++)
upsilon += U[i][0] * sin(U[i][1]);
upsilon *= 1e-7;
// now use the "time completed" series
const double n = 21.57107117668 * deg_to_rad;
computeArguments(t + upsilon/n, l1, l2, l3, l4,
om1, om2, om3, om4, psi, Gp, G);
const int NZ = 18;
double Z[NZ][2] = {
{ 8, l4 - 2*PIj - om4 - 2*psi },
{ 8, l4 - 2*PIj + psi - 4*G },
{ 88, l4 - 2*PIj + psi - 3*G },
{ 773, l4 - 2*PIj + psi - 2*G },
{ -38, l4 - 2*PIj + psi - G },
{ 5, l4 - 2*PIj + psi },
{ 9, l4 - om1 },
{ -17, l4 - om2 },
{ -5112, l4 - om3 },
{ -7, l4 - om4 - G },
{ 44134, l4 - om4 },
{ 7, l4 - om4 + G },
{ -102, l4 - psi - G },
{ -76579, l4 - psi },
{ 104, l4 - psi + G },
{ -10, l4 - psi + 5*Gp - 2*G + phi2 },
{ -11, l3 - 2*l4 + psi },
{ 7, l3 - 2*l4 + om4 }
};
for (int i = 0; i < NZ; i++)
zeta += Z[i][0] * sin(Z[i][1]);
zeta *= 1e-7;
}
break;
default:
xpExit("Unknown Jupiter satellite\n", __FILE__, __LINE__);
}
// Jupiter equatorial coordinates
X = radius * cos(lon - psi + upsilon);
Y = radius * sin(lon - psi + upsilon);
Z = radius * zeta;
// rotate to Jupiter's orbital plane
const double I = 3.10401 * deg_to_rad;
rotateX(X, Y, Z, -I);
// rotate towards ascending node of Jupiter's equator on its
// orbital plane
const double OM = 99.95326 * deg_to_rad;
rotateZ(X, Y, Z, OM - psi);
// rotate to ecliptic
const double J = 1.30691 * deg_to_rad;
rotateX(X, Y, Z, -J);
// rotate towards ascending node of Jupiter's orbit on ecliptic
rotateZ(X, Y, Z, -OM);
// rotate to earth equator B1950
const double eps = 23.4457889 * deg_to_rad;
rotateX(X, Y, Z, -eps);
// precess to J2000
precessB1950J2000(X, Y, Z);
}
syntax highlighted by Code2HTML, v. 0.9.1