Jupiter moons calculations are now consistent with Stellarium.

I translated the L1.2 FORTRAN code into C, and verified that the calculations match the Stellarium code I modified to produce EQJ coordinates. I still need to compare against JPL Horizons data.
2026-04-05 07:15:38 -04:00 · 2021-04-10 16:29:10 -04:00
parent 70a9fa5e27
commit dcbd6fe243
7 changed files with 939 additions and 516 deletions
--- a/generate/codegen.c
+++ b/generate/codegen.c
@@ -1387,12 +1387,26 @@ fail:

 typedef struct
 {
-    double mu;      /* mu = G(M+m), where M = Jupiter mass, m = moon mass. */
-    double al[2];   /* mean longitude coefficients */
+    double          mu;          /* mu = G(M+m), where M = Jupiter mass, m = moon mass. */
+    double          al[2];       /* mean longitude coefficients */
    vsop_series_t   a;
    vsop_series_t   l;
    vsop_series_t   z;
    vsop_series_t   zeta;
+}
+jupiter_moon_t;
+
+
+typedef struct
+{
+    /* angles that rotate Jupiter equatorial system to EQJ */
+    double          incl;
+    double          psi;
+
+    /* a rotation matrix calculated from incl and psi, for the convenience of the code generators. */
+    double          rot[3][3];
+
+    jupiter_moon_t  moon[NUM_JUPITER_MOONS];
    vsop_term_t     buffer[MAX_JUPITER_TERMS];
 }
 jupiter_moon_model_t;
@@ -1419,13 +1433,13 @@ static int ReadFixLine(char *line, size_t size, FILE *infile)
    do{if ((nscanned) != (required)) FAIL("LoadJupiterMoonModel(%s line %d): expected %d tokens, found %d\n", filename, lnum, required, nscanned);}while(0)


-static int LoadJupiterMoonModel(jupiter_moon_model_t model[NUM_JUPITER_MOONS], const char *filename)
+static int LoadJupiterMoonModel(const char *filename, jupiter_moon_model_t *model)
 {
    int error;
    FILE *infile;
    int lnum;
    char line[200];
-    int nscanned, check, moon, nterms, tindex, var, next_term_index, read_mean_longitude;
+    int nscanned, check, mindex, nterms, tindex, var, next_term_index, read_mean_longitude;
    vsop_series_t *series = NULL;

    infile = fopen(filename, "rt");
@@ -1433,7 +1447,7 @@ static int LoadJupiterMoonModel(jupiter_moon_model_t model[NUM_JUPITER_MOONS], c
        FAIL("LoadJupiterMoonModel: Cannot open file for read: %s\n", filename);

    lnum = 0;
-    moon = -1;
+    mindex = -1;
    nterms = 0;
    tindex = 0;
    var = 3;
@@ -1445,9 +1459,33 @@ static int LoadJupiterMoonModel(jupiter_moon_model_t model[NUM_JUPITER_MOONS], c
        if (lnum == 19)
        {
            /* G(M+m), where M = Jupiter mass, m = moon mass. */
-            nscanned = sscanf(line, "%lf %lf %lf %lf", &model[0].mu, &model[1].mu, &model[2].mu, &model[3].mu);
+            nscanned = sscanf(line, "%lf %lf %lf %lf", &model->moon[0].mu, &model->moon[1].mu, &model->moon[2].mu, &model->moon[3].mu);
            REQUIRE_NSCAN(4);
        }
+        else if (lnum == 20)
+        {
+            double cp, sp, ci, si;
+
+            /* read angles that convert Jupiter equatorial coordinates into Earth J2000 equatorial coordinates. */
+            nscanned = sscanf(line, "%lf %lf", &model->psi, &model->incl);
+            REQUIRE_NSCAN(2);
+
+            /* calculate the corresponding rotation matrix. */
+            ci = cos(model->incl);
+            si = sin(model->incl);
+            cp = cos(model->psi);
+            sp = sin(model->psi);
+
+            model->rot[0][0] = +cp;
+            model->rot[1][0] = -sp*ci;
+            model->rot[2][0] = +sp*si;
+            model->rot[0][1] = +sp;
+            model->rot[1][1] = +cp*ci;
+            model->rot[2][1] = -cp*si;
+            model->rot[0][2] = 0.0;
+            model->rot[1][2] = +si;
+            model->rot[2][2] = +ci;
+        }
        else if (lnum >= 21)
        {
            if (tindex < nterms)
@@ -1455,15 +1493,11 @@ static int LoadJupiterMoonModel(jupiter_moon_model_t model[NUM_JUPITER_MOONS], c
                if (tindex == 0 && var == 1 && !read_mean_longitude)
                {
                    /* Special case: there is an extra row of mean longitude values here. */
-                    nscanned = sscanf(line, "%d %lf %lf",
-                        &check,
-                        &model[moon].al[0],
-                        &model[moon].al[1]);
-
-                    REQUIRE_NSCAN(3);
-                    if (check != 0)
-                        FAIL("LoadJupiterMoonModel(%s line %d): invalid check value %d (expected 0 for mean longitude line)\n", filename, lnum, check);
+                    nscanned = sscanf(line, "%lf %lf",
+                        &model->moon[mindex].al[0],
+                        &model->moon[mindex].al[1]);

+                    REQUIRE_NSCAN(2);
                    read_mean_longitude = 1;
                }
                else
@@ -1491,8 +1525,8 @@ static int LoadJupiterMoonModel(jupiter_moon_model_t model[NUM_JUPITER_MOONS], c
                {
                    /* Start the next moon. */
                    var = 0;
-                    ++moon;
-                    if (moon == NUM_JUPITER_MOONS)
+                    ++mindex;
+                    if (mindex == NUM_JUPITER_MOONS)
                        break;  /* We are done loading data! (We don't care about the Tchebycheff polynomials.) */
                }
                else
@@ -1504,10 +1538,10 @@ static int LoadJupiterMoonModel(jupiter_moon_model_t model[NUM_JUPITER_MOONS], c
                /* Select which variable corresponds to the value of 'var'. */
                switch (var)
                {
-                    case 0:  series = &model[moon].a;     break;
-                    case 1:  series = &model[moon].l;     break;
-                    case 2:  series = &model[moon].z;     break;
-                    case 3:  series = &model[moon].zeta;  break;
+                    case 0:  series = &model->moon[mindex].a;     break;
+                    case 1:  series = &model->moon[mindex].l;     break;
+                    case 2:  series = &model->moon[mindex].z;     break;
+                    case 3:  series = &model->moon[mindex].zeta;  break;
                    default:
                        FAIL("LoadJupiterMoonModel: Invalid var = %d\n", var);
                }
@@ -1528,8 +1562,8 @@ static int LoadJupiterMoonModel(jupiter_moon_model_t model[NUM_JUPITER_MOONS], c
        }
    }

-    if (moon != NUM_JUPITER_MOONS)
-        FAIL("LoadJupiterMoonModel(%s): only loaded %d moons.\n", filename, moon);
+    if (mindex != NUM_JUPITER_MOONS)
+        FAIL("LoadJupiterMoonModel(%s): only loaded %d moons.\n", filename, mindex);

    error = 0;
 fail:
@@ -1538,31 +1572,42 @@ fail:
 }


-static int JupiterMoons_C(cg_context_t *context, const jupiter_moon_model_t model[NUM_JUPITER_MOONS])
+static int JupiterMoons_C(cg_context_t *context, const jupiter_moon_model_t *model)
 {
-    int moon, var, i;
+    int mindex, var, i;
    const char *moon_name[] = { "Io", "Europa", "Ganymede", "Callisto" };
    const char *var_name[] = { "a", "l", "z", "zeta" };
    vsop_series_t series[4];

-    for (moon = 0; moon < NUM_JUPITER_MOONS; ++moon)
+    fprintf(context->outfile, "static const astro_rotation_t Rotation_JUP_EQJ =\n");
+    fprintf(context->outfile, "{\n");
+    fprintf(context->outfile, "    ASTRO_SUCCESS,\n");
+    fprintf(context->outfile, "    {\n");
+    fprintf(context->outfile, "        { %23.16le, %23.16le, %23.16le },\n", model->rot[0][0], model->rot[0][1], model->rot[0][2]);
+    fprintf(context->outfile, "        { %23.16le, %23.16le, %23.16le },\n", model->rot[1][0], model->rot[1][1], model->rot[1][2]);
+    fprintf(context->outfile, "        { %23.16le, %23.16le, %23.16le }\n",  model->rot[2][0], model->rot[2][1], model->rot[2][2]);
+    fprintf(context->outfile, "    }\n");
+    fprintf(context->outfile, "};\n\n");
+
+    for (mindex = 0; mindex < NUM_JUPITER_MOONS; ++mindex)
    {
-        series[0] = model[moon].a;
-        series[1] = model[moon].l;
-        series[2] = model[moon].z;
-        series[3] = model[moon].zeta;
+        series[0] = model->moon[mindex].a;
+        series[1] = model->moon[mindex].l;
+        series[2] = model->moon[mindex].z;
+        series[3] = model->moon[mindex].zeta;
        for (var = 0; var < NUM_JM_VARS; ++var)
        {
-            fprintf(context->outfile, "static const vsop_term_t jm_%s_%s[] =\n", moon_name[moon], var_name[var]);
+            int nterms = series[var].nterms_total;
+            fprintf(context->outfile, "static const vsop_term_t jm_%s_%s[] =\n", moon_name[mindex], var_name[var]);
            fprintf(context->outfile, "{\n");
-            for (i = 0; i < series->nterms_total; ++i)
+            for (i = 0; i < nterms; ++i)
            {
-                const vsop_term_t *term = &series->term[i];
+                const vsop_term_t *term = &series[var].term[i];
                fprintf(context->outfile, "    { %23.16le, %23.16le, %23.16le }%s\n",
                    term->amplitude,
                    term->phase,
                    term->frequency,
-                    (i + 1 < series->nterms_total) ? "," : "");
+                    (i+1 < nterms) ? "," : "");
            }
            fprintf(context->outfile, "};\n\n");
        }
@@ -1570,35 +1615,35 @@ static int JupiterMoons_C(cg_context_t *context, const jupiter_moon_model_t mode

    fprintf(context->outfile, "static jupiter_moon_t JupiterMoonModel[%d] =\n", NUM_JUPITER_MOONS);
    fprintf(context->outfile, "{\n");
-    for (moon = 0; moon < NUM_JUPITER_MOONS; ++moon)
+    for (mindex = 0; mindex < NUM_JUPITER_MOONS; ++mindex)
    {
-        series[0] = model[moon].a;
-        series[1] = model[moon].l;
-        series[2] = model[moon].z;
-        series[3] = model[moon].zeta;
-        fprintf(context->outfile, "    { %0.16le, {%0.16le, %0.16le}", model[moon].mu, model[moon].al[0], model[moon].al[1]);
+        series[0] = model->moon[mindex].a;
+        series[1] = model->moon[mindex].l;
+        series[2] = model->moon[mindex].z;
+        series[3] = model->moon[mindex].zeta;
+        fprintf(context->outfile, "    { %0.16le, {%0.16le, %0.16le}", model->moon[mindex].mu, model->moon[mindex].al[0], model->moon[mindex].al[1]);
        for (var = 0; var < NUM_JM_VARS; ++var)
-            fprintf(context->outfile, ", {%d, jm_%s_%s}", series[var].nterms_total, moon_name[moon], var_name[var]);
-        fprintf(context->outfile, " }%s\n", ((moon + 1 < NUM_JUPITER_MOONS) ? "," : ""));
+            fprintf(context->outfile, ", {%d, jm_%s_%s}", series[var].nterms_total, moon_name[mindex], var_name[var]);
+        fprintf(context->outfile, " }%s\n", ((mindex + 1 < NUM_JUPITER_MOONS) ? "," : ""));
    }
    fprintf(context->outfile, "}");
    return 0;
 }


-static int JupiterMoons_CSharp(cg_context_t *context, const jupiter_moon_model_t model[NUM_JUPITER_MOONS])
+static int JupiterMoons_CSharp(cg_context_t *context, const jupiter_moon_model_t *model)
 {
    return LogError(context, "JupiterMoons_CSharp: NOT YET IMPLEMENTED.\n");
 }


-static int JupiterMoons_JS(cg_context_t *context, const jupiter_moon_model_t model[NUM_JUPITER_MOONS])
+static int JupiterMoons_JS(cg_context_t *context, const jupiter_moon_model_t *model)
 {
    return LogError(context, "JupiterMoons_JS: NOT YET IMPLEMENTED.\n");
 }


-static int JupiterMoons_Python(cg_context_t *context, const jupiter_moon_model_t model[NUM_JUPITER_MOONS])
+static int JupiterMoons_Python(cg_context_t *context, const jupiter_moon_model_t *model)
 {
    return LogError(context, "JupiterMoons_Python: NOT YET IMPLEMENTED.\n");
 }
@@ -1609,11 +1654,11 @@ static int JupiterMoons(cg_context_t *context)
    int error;
    jupiter_moon_model_t *model;

-    model = calloc(NUM_JUPITER_MOONS, sizeof(jupiter_moon_model_t));
+    model = calloc(1, sizeof(jupiter_moon_model_t));
    if (model == NULL)
        FAIL("JupiterMoons: memory allocation failure!\n");

-    CHECK(LoadJupiterMoonModel(model, "jupiter_moons/fortran/BisL1.2.dat"));
+    CHECK(LoadJupiterMoonModel("jupiter_moons/fortran/BisL1.2.dat", model));

    switch (context->language)
    {
--- a/generate/ctest.c
+++ b/generate/ctest.c
@@ -3608,9 +3608,14 @@ static int JupiterMoons_CompareStellarium(void)
    char line[200];
    const char *inFileName = "jupiter_moons/horizons/stellarium.txt";
    double prev_tt, tt, pos[3], vel[3];
+    double dx, dy, dz, diff;
    int moon, nscanned;
    astro_time_t time;
    astro_jupiter_moons_t jm;
+    const double pos_tolerance = 1.0e-12;
+    const double vel_tolerance = 1.0e-12;
+
+    memset(&jm, 0, sizeof(jm));

    infile = fopen(inFileName, "rt");
    if (infile == NULL)
@@ -3624,15 +3629,39 @@ static int JupiterMoons_CompareStellarium(void)
        if (lnum == 1) continue;    /* ignore header line */
        nscanned = sscanf(line, "%lf %d %lf %lf %lf %lf %lf %lf", &tt, &moon, &pos[0], &pos[1], &pos[2], &vel[0], &vel[1], &vel[2]);
        if (nscanned != 8)
-            FAIL("JupiterMoons_CompareStellarium(%s line %d): error scanning data\n", inFileName, lnum);
+            FAIL("C JupiterMoons_CompareStellarium(%s line %d): error scanning data\n", inFileName, lnum);
+
+        if (moon < 0 || moon > 3)
+            FAIL("C JupiterMoons_CompareStellarium(%s line %d): invalid moon = %d\n", inFileName, lnum, moon);
+
        if (tt != prev_tt)
        {
            /* Calculate all 4 moons for this moment in time. */
            prev_tt = tt;
            time = AstroTerrestrialTime(tt);
            jm = Astronomy_JupiterMoons(time);
-            (void)jm;
+            CHECK_STATUS(jm.moon[0]);
+            CHECK_STATUS(jm.moon[1]);
+            CHECK_STATUS(jm.moon[2]);
+            CHECK_STATUS(jm.moon[3]);
        }
+
+        DEBUG("C JupiterMoons_CompareStellarium: tt=%6.1lf ref  pos=(%20.12le, %20.12le, %20.12le), vel=(%20.12le, %20.12le, %20.12le)\n", tt, pos[0], pos[1], pos[2], vel[0], vel[1], vel[2]);
+        DEBUG("              moon=%d                         calc pos=(%20.12le, %20.12le, %20.12le), vel=(%20.12le, %20.12le, %20.12le)\n\n", moon, jm.moon[moon].x, jm.moon[moon].y, jm.moon[moon].z, jm.moon[moon].vx, jm.moon[moon].vy, jm.moon[moon].vz);
+
+        dx = V(pos[0] - jm.moon[moon].x);
+        dy = V(pos[1] - jm.moon[moon].y);
+        dz = V(pos[2] - jm.moon[moon].z);
+        diff = sqrt(dx*dx + dy*dy + dz*dz);
+        if (diff > pos_tolerance)
+            FAIL("C JupiterMoons_CompareStellarium(%s line %d): EXCESSIVE POSITION ERROR: diff = %le\n", inFileName, lnum, diff);
+
+        dx = V(vel[0] - jm.moon[moon].vx);
+        dy = V(vel[1] - jm.moon[moon].vy);
+        dz = V(vel[2] - jm.moon[moon].vz);
+        diff = sqrt(dx*dx + dy*dy + dz*dz);
+        if (diff > vel_tolerance)
+            FAIL("C JupiterMoons_CompareStellarium(%s line %d): EXCESSIVE VELOCITY ERROR: diff = %le\n", inFileName, lnum, diff);
    }

    DEBUG("C JupiterMoons_Stellarium: PASS\n");
--- a/generate/jupiter_moons/fortran/L1.2.f
+++ b/generate/jupiter_moons/fortran/L1.2.f
@@ -441,10 +441,10 @@
        a= T1    !  -819.727638594856D0
        b= T2    !   812.721806990360D0
        x=(T/365.25D0-0.5D0*(b+a))/(0.5D0*(b-a))
-	if(abs(x).gt.1.d0) then
-	  write(*,*)'date',T,' out of interval of validity'
-	  return
-	end if
+        if(abs(x).gt.1.d0) then
+          write(*,*)'date',T,' out of interval of validity'
+	        return
+	      end if
        TN(0)=1.D0
        TN(1)=x
        do it=2,8
@@ -463,42 +463,43 @@
        kv=1
        s=0.d0
        do k=1,nbterm(ks,kv)
-         arg=phasc(k,ks,kv)+freqc(k,ks,kv)*T
-         s=s+ampl(k,ks,kv)*cos(arg)
-	end do
-	elem(1)=s
+          arg=phasc(k,ks,kv)+freqc(k,ks,kv)*T
+          s=s+ampl(k,ks,kv)*cos(arg)
+        end do
+        elem(1)=s

        kv=2
-         s=al(1,ks) +al(2,ks)*T
+        s = al(1,ks) + al(2,ks)*T
        do k=1,nbterm(ks,kv)
-         arg=phasc(k,ks,kv)+freqc(k,ks,kv)*T
-         s=s+ampl(k,ks,kv)*sin(arg)
-	end do
-	s=mod(s+val(1),pi2)
-	if(s.lt.0) s=s+pi2
-	elem(2)=s
+          arg=phasc(k,ks,kv)+freqc(k,ks,kv)*T
+          s=s+ampl(k,ks,kv)*sin(arg)
+      	end do
+        s=mod(s+val(1),pi2)
+        if(s.lt.0) s=s+pi2
+        elem(2)=s

        kv=3
-         s1=0.d0
-         s2=0.d0
+        s1=0.d0
+        s2=0.d0
        do k=1,nbterm(ks,kv)
-         arg=phasc(k,ks,kv)+freqc(k,ks,kv)*T
-         s1=s1+ampl(k,ks,kv)*cos(arg)
-         s2=s2+ampl(k,ks,kv)*sin(arg)
-	end do
-	elem(3)=s1+val(2)
-	elem(4)=s2+val(3)
+          arg=phasc(k,ks,kv)+freqc(k,ks,kv)*T
+          s1=s1+ampl(k,ks,kv)*cos(arg)
+          s2=s2+ampl(k,ks,kv)*sin(arg)
+	      end do
+	      elem(3)=s1+val(2)
+	      elem(4)=s2+val(3)

        kv=4
-         s1=0.d0
-         s2=0.d0
+        s1=0.d0
+        s2=0.d0
        do k=1,nbterm(ks,kv)
         arg=phasc(k,ks,kv)+freqc(k,ks,kv)*T
         s1=s1+ampl(k,ks,kv)*cos(arg)
         s2=s2+ampl(k,ks,kv)*sin(arg)
-	end do
-	elem(5)=s1+val(4)
-	elem(6)=s2+val(5)
+	      end do
+
+	      elem(5)=s1+val(4)
+	      elem(6)=s2+val(5)
 *
 *  computing cartesian coordinates from elements
 *
@@ -571,7 +572,7 @@
 *  OMEGA  is the longitude of ascending node on Fxy plane
 *  omega  is the argument of pericentre such as
 *         OMEGA + omega represents the longitude of pericentre
-*     and  L (= OMEGA+omega+M) repr<70>sents the mean longitude
+*     and  L (= OMEGA+omega+M) repr<70>sents the mean longitude
 *
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(PI2=3.141592653589793D0*2.d0)
--- a/generate/template/astronomy.c
+++ b/generate/template/astronomy.c
@@ -345,6 +345,16 @@ static astro_vector_t VecError(astro_status_t status, astro_time_t time)
    return vec;
 }

+static astro_state_vector_t StateVecError(astro_status_t status, astro_time_t time)
+{
+    astro_state_vector_t vec;
+    vec.x = vec.y = vec.z = NAN;
+    vec.vx = vec.vy = vec.vz = NAN;
+    vec.t = time;
+    vec.status = status;
+    return vec;
+}
+
 static astro_spherical_t SphereError(astro_status_t status)
 {
    astro_spherical_t sphere;
@@ -2329,13 +2339,130 @@ static astro_vector_t CalcPluto(astro_time_t time)

 $ASTRO_JUPITER_MOONS();

+static astro_state_vector_t JupiterMoon_elem2pv(astro_time_t time, double mu, double elem[6])
+{
+    /* Translation of FORTRAN subroutine ELEM2PV from: */
+    /* https://ftp.imcce.fr/pub/ephem/satel/galilean/L1/L1.2/ */
+    astro_state_vector_t state;
+    double EE, DE, CE, SE, DLE, RSAM1, ASR, PHI, PSI, X1, Y1, VX1, VY1, F2, P2, Q2, PQ;
+
+    const double A  = elem[0];
+    const double AL = elem[1];
+    const double K  = elem[2];
+    const double H  = elem[3];
+    const double Q  = elem[4];
+    const double P  = elem[5];
+
+    const double AN = sqrt(mu / (A*A*A));
+
+    EE = AL + K*sin(AL) - H*cos(AL);
+    do
+    {
+        CE = cos(EE);
+        SE = sin(EE);
+        DE = (AL - EE + K*SE - H*CE) / (1.0 - K*CE - H*SE);
+        EE += DE;
+    }
+    while (fabs(DE) >= 1.0e-12);
+
+    CE = cos(EE);
+    SE = sin(EE);
+    DLE = H*CE - K*SE;
+    RSAM1 = -K*CE - H*SE;
+    ASR = 1.0/(1.0 + RSAM1);
+    PHI = sqrt(1.0 - K*K - H*H);
+    PSI = 1.0/(1.0 + PHI);
+    X1 = A*(CE - K - PSI*H*DLE);
+    Y1 = A*(SE - H + PSI*K*DLE);
+    VX1 = AN*ASR*A*(-SE - PSI*H*RSAM1);
+    VY1 = AN*ASR*A*(+CE + PSI*K*RSAM1);
+    F2 = 2.0*sqrt(1.0 - Q*Q - P*P);
+    P2 = 1.0 - 2.0*P*P;
+    Q2 = 1.0 - 2.0*Q*Q;
+    PQ = 2.0*P*Q;
+
+    state.x = X1*P2 + Y1*PQ;
+    state.y = X1*PQ + Y1*Q2;
+    state.z = (Q*Y1 - X1*P)*F2;
+
+    state.vx = VX1*P2 + VY1*PQ;
+    state.vy = VX1*PQ + VY1*Q2;
+    state.vz = (Q*VY1 - VX1*P)*F2;
+
+    state.t = time;
+    state.status = ASTRO_SUCCESS;
+    return state;
+}
+
+static astro_state_vector_t CalcJupiterMoon(astro_time_t time, int mindex)
+{
+    /* This is a translation of FORTRAN code by Duriez, Lainey, and Vienne: */
+    /* https://ftp.imcce.fr/pub/ephem/satel/galilean/L1/L1.2/ */
+
+    astro_state_vector_t state;
+    int k;
+    double arg;
+    double elem[6];
+    const jupiter_moon_t *m = &JupiterMoonModel[mindex];
+    const double t = time.tt + 18262.5;     /* t = time since 1950-01-01T00:00:00Z */
+
+    /* Calculate 6 orbital elements at the given time t. */
+
+    elem[0] = 0.0;
+    for (k = 0; k < m->a.nterms; ++k)
+    {
+        arg = m->a.term[k].phase + (t * m->a.term[k].frequency);
+        elem[0] += m->a.term[k].amplitude * cos(arg);
+    }
+
+    elem[1] = m->al[0] + (t * m->al[1]);
+    for (k = 0; k < m->l.nterms; ++k)
+    {
+        arg = m->l.term[k].phase + (t * m->l.term[k].frequency);
+        elem[1] += m->l.term[k].amplitude * sin(arg);
+    }
+    elem[1] = fmod(elem[1], PI2);
+    if (elem[1] < 0.0)
+        elem[1] += PI2;
+
+    elem[2] = elem[3] = 0.0;
+    for (k = 0; k < m->z.nterms; ++k)
+    {
+        arg = m->z.term[k].phase + (t * m->z.term[k].frequency);
+        elem[2] += m->z.term[k].amplitude * cos(arg);
+        elem[3] += m->z.term[k].amplitude * sin(arg);
+    }
+
+    elem[4] = elem[5] = 0.0;
+    for (k = 0; k < m->zeta.nterms; ++k)
+    {
+        arg = m->zeta.term[k].phase + (t * m->zeta.term[k].frequency);
+        elem[4] += m->zeta.term[k].amplitude * cos(arg);
+        elem[5] += m->zeta.term[k].amplitude * sin(arg);
+    }
+
+    /* Convert the oribital elements into position vectors in the Jupiter equatorial system (JUP). */
+    state = JupiterMoon_elem2pv(time, m->mu, elem);
+
+    /* Re-orient position and velocity vectors from Jupiter-equatorial (JUP) to Earth-equatorial in J2000 (EQJ). */
+    return Astronomy_RotateState(Rotation_JUP_EQJ, state);
+}
+
+
 /**
- * @brief Calculates positions of Jupiter's largest 4 moons.
+ * @brief Calculates jovicentric positions of Jupiter's largest 4 moons.
 *
- * Calculates jovicentric position vectors for Jupiter's
- * moons Io, Europa, Ganymede, and Callisto, at the given date and time.
- * See #astro_jupiter_moons_t for more information about the representation
- * of the position vectors.
+ * Calculates position vectors for Jupiter's moons
+ * Io, Europa, Ganymede, and Callisto, at the given date and time.
+ * The position vectors are jovicentric, meaning their coordinate origin
+ * is the center of Jupiter. Their orientation is the Earth's equatorial
+ * system at the J2000 epoch, called `EQJ`. The vector components
+ * are expressed in astronomical units (AU).
+ *
+ * To convert to heliocentric vectors, call #Astronomy_HelioVector
+ * with `BODY_JUPITER` to get Jupiter's heliocentric position, then
+ * add the jovicentric vectors. Likewise, you can call #Astronomy_GeoVector
+ * with `BODY_JUPITER` to convert to geocentric vectors.
 *
 * @param time  The date and time for which to calculate the position vectors.
 * @return Position vectors of Jupiter's largest 4 moons, as described above.
@@ -2343,13 +2470,10 @@ $ASTRO_JUPITER_MOONS();
 astro_jupiter_moons_t Astronomy_JupiterMoons(astro_time_t time)
 {
    astro_jupiter_moons_t jm;
+    int mindex;

-    (void)JupiterMoonModel;
-
-    jm.moon[0] = VecError(ASTRO_NOT_INITIALIZED, time);
-    jm.moon[1] = VecError(ASTRO_NOT_INITIALIZED, time);
-    jm.moon[2] = VecError(ASTRO_NOT_INITIALIZED, time);
-    jm.moon[3] = VecError(ASTRO_NOT_INITIALIZED, time);
+    for (mindex = 0; mindex < NUM_JUPITER_MOONS; ++mindex)
+        jm.moon[mindex] = CalcJupiterMoon(time, mindex);

    return jm;
 }
@@ -5889,6 +6013,45 @@ astro_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t
 }


+/**
+ * @brief
+ *      Applies a rotation to a state vector, yielding a rotated vector.
+ *
+ * This function transforms a state vector in one orientation to a vector
+ * in another orientation.
+ *
+ * @param rotation
+ *      A rotation matrix that specifies how the orientation of the vector is to be changed.
+ *
+ * @param state
+ *      The state vector whose orientation is to be changed.
+ *      Both the position and velocity components are transformed.
+ *
+ * @return
+ *      A state vector in the orientation specified by `rotation`.
+ */
+astro_state_vector_t Astronomy_RotateState(astro_rotation_t rotation, astro_state_vector_t state)
+{
+    astro_state_vector_t target;
+
+    if (rotation.status != ASTRO_SUCCESS || state.status != ASTRO_SUCCESS)
+        return StateVecError(ASTRO_INVALID_PARAMETER, state.t);
+
+    target.status = ASTRO_SUCCESS;
+    target.t = state.t;
+
+    target.x = rotation.rot[0][0]*state.x + rotation.rot[1][0]*state.y + rotation.rot[2][0]*state.z;
+    target.y = rotation.rot[0][1]*state.x + rotation.rot[1][1]*state.y + rotation.rot[2][1]*state.z;
+    target.z = rotation.rot[0][2]*state.x + rotation.rot[1][2]*state.y + rotation.rot[2][2]*state.z;
+
+    target.vx = rotation.rot[0][0]*state.vx + rotation.rot[1][0]*state.vy + rotation.rot[2][0]*state.vz;
+    target.vy = rotation.rot[0][1]*state.vx + rotation.rot[1][1]*state.vy + rotation.rot[2][1]*state.vz;
+    target.vz = rotation.rot[0][2]*state.vx + rotation.rot[1][2]*state.vy + rotation.rot[2][2]*state.vz;
+
+    return target;
+}
+
+
 /**
 * @brief
 *      Calculates a rotation matrix from equatorial J2000 (EQJ) to ecliptic J2000 (ECL).
--- a/source/c/README.md
+++ b/source/c/README.md
@@ -822,11 +822,13 @@ Given a rotation matrix that performs some coordinate transform, this function r
 <a name="Astronomy_JupiterMoons"></a>
 ### Astronomy_JupiterMoons(time) &#8658; [`astro_jupiter_moons_t`](#astro_jupiter_moons_t)

-**Calculates positions of Jupiter's largest 4 moons.** 
+**Calculates jovicentric positions of Jupiter's largest 4 moons.** 



-Calculates jovicentric position vectors for Jupiter's moons Io, Europa, Ganymede, and Callisto, at the given date and time. See [`astro_jupiter_moons_t`](#astro_jupiter_moons_t) for more information about the representation of the position vectors.
+Calculates position vectors for Jupiter's moons Io, Europa, Ganymede, and Callisto, at the given date and time. The position vectors are jovicentric, meaning their coordinate origin is the center of Jupiter. Their orientation is the Earth's equatorial system at the J2000 epoch, called `EQJ`. The vector components are expressed in astronomical units (AU).
+
+To convert to heliocentric vectors, call [`Astronomy_HelioVector`](#Astronomy_HelioVector) with `BODY_JUPITER` to get Jupiter's heliocentric position, then add the jovicentric vectors. Likewise, you can call [`Astronomy_GeoVector`](#Astronomy_GeoVector) with `BODY_JUPITER` to convert to geocentric vectors.



@@ -1232,6 +1234,31 @@ Given an altitude angle and a refraction option, calculates the amount of "lift"

 Astronomy Engine uses dynamic memory allocation in only one place: it makes calculation of Pluto's orbit more efficient by caching 11 KB segments and recycling them. To force purging this cache and freeing all the dynamic memory, you can call this function at any time. It is always safe to call, although it will slow down the very next calculation of Pluto's position for a nearby time value. Calling this function before your program exits is optional, but it will be helpful for leak-checkers like valgrind. 

+---
+
+<a name="Astronomy_RotateState"></a>
+### Astronomy_RotateState(rotation, state) &#8658; [`astro_state_vector_t`](#astro_state_vector_t)
+
+**Applies a rotation to a state vector, yielding a rotated vector.** 
+
+
+
+This function transforms a state vector in one orientation to a vector in another orientation.
+
+
+
+**Returns:**  A state vector in the orientation specified by `rotation`. 
+
+
+
+| Type | Parameter | Description |
+| --- | --- | --- |
+| [`astro_rotation_t`](#astro_rotation_t) | `rotation` |  A rotation matrix that specifies how the orientation of the vector is to be changed. | 
+| [`astro_state_vector_t`](#astro_state_vector_t) | `state` |  The state vector whose orientation is to be changed. Both the position and velocity components are transformed. | 
+
+
+
+
 ---

 <a name="Astronomy_RotateVector"></a>
@@ -2897,13 +2924,13 @@ Returned by the functions [`Astronomy_Illumination`](#Astronomy_Illumination) an



-The [`Astronomy_JupiterMoons`](#Astronomy_JupiterMoons) function returns this struct to report position vectors for Jupiter's largest 4 moons Io, Europa, Ganymede, and Callisto. Each vector is relative to the center of Jupiter and is oriented in the EQJ system (that is, using Earth's equator at the J2000 epoch.) The vector components are expressed in astronomical units (AU).
+The [`Astronomy_JupiterMoons`](#Astronomy_JupiterMoons) function returns this struct to report position and velocity vectors for Jupiter's largest 4 moons Io, Europa, Ganymede, and Callisto. Each vector is relative to the center of Jupiter and is oriented in the EQJ system (that is, using Earth's equator at the J2000 epoch.) The positions are expressed in astronomical units (AU), and the velocities in AU/day.

 The following integer constants may be useful for indexing into the `moon` array: [`JM_IO`](#JM_IO), [`JM_EUROPA`](#JM_EUROPA), [`JM_GANYMEDE`](#JM_GANYMEDE), [`JM_CALLISTO`](#JM_CALLISTO). 

 | Type | Member | Description |
 | ---- | ------ | ----------- |
-| [`astro_vector_t`](#astro_vector_t) | `moon` |  Jovicentric coordinates of each moon, as described above.  |
+| [`astro_state_vector_t`](#astro_state_vector_t) | `moon` |  Jovicentric position and velocity of each moon, as described above.  |


 ---
@@ -3062,6 +3089,27 @@ You can create this structure directly, or you can call the convenience function
 | `double` | `dist` |  Distance in AU.  |


+---
+
+<a name="astro_state_vector_t"></a>
+### `astro_state_vector_t`
+
+**A state vector that contains a position (AU) and velocity (AU/day).** 
+
+
+
+| Type | Member | Description |
+| ---- | ------ | ----------- |
+| [`astro_status_t`](#astro_status_t) | `status` |  `ASTRO_SUCCESS` if this struct is valid; otherwise an error code.  |
+| `double` | `x` |  The Cartesian position x-coordinate of the vector in AU.  |
+| `double` | `y` |  The Cartesian position y-coordinate of the vector in AU.  |
+| `double` | `z` |  The Cartesian position z-coordinate of the vector in AU.  |
+| `double` | `vx` |  The Cartesian velocity x-coordinate of the vector in AU/day.  |
+| `double` | `vy` |  The Cartesian velocity y-coordinate of the vector in AU/day.  |
+| `double` | `vz` |  The Cartesian velocity z-coordinate of the vector in AU/day.  |
+| [`astro_time_t`](#astro_time_t) | `t` |  The date and time at which this state vector is valid.  |
+
+
 ---

 <a name="astro_time_t"></a>
--- a/source/c/astronomy.c
+++ b/source/c/astronomy.c
--- a/source/c/astronomy.h
+++ b/source/c/astronomy.h
@@ -217,6 +217,22 @@ typedef struct
 }
 astro_vector_t;

+/**
+ * @brief A state vector that contains a position (AU) and velocity (AU/day).
+ */
+typedef struct
+{
+    astro_status_t status;  /**< `ASTRO_SUCCESS` if this struct is valid; otherwise an error code. */
+    double x;               /**< The Cartesian position x-coordinate of the vector in AU. */
+    double y;               /**< The Cartesian position y-coordinate of the vector in AU. */
+    double z;               /**< The Cartesian position z-coordinate of the vector in AU. */
+    double vx;              /**< The Cartesian velocity x-coordinate of the vector in AU/day. */
+    double vy;              /**< The Cartesian velocity y-coordinate of the vector in AU/day. */
+    double vz;              /**< The Cartesian velocity z-coordinate of the vector in AU/day. */
+    astro_time_t t;         /**< The date and time at which this state vector is valid. */
+}
+astro_state_vector_t;
+
 /**
 * @brief Spherical coordinates: latitude, longitude, distance.
 */
@@ -865,18 +881,19 @@ astro_time_format_t;
 * @brief Holds the positions of Jupiter's major 4 moons.
 *
 * The #Astronomy_JupiterMoons function returns this struct
- * to report position vectors for Jupiter's largest 4 moons
+ * to report position and velocity vectors for Jupiter's largest 4 moons
 * Io, Europa, Ganymede, and Callisto. Each vector is relative
 * to the center of Jupiter and is oriented in the EQJ system
 * (that is, using Earth's equator at the J2000 epoch.)
- * The vector components are expressed in astronomical units (AU).
+ * The positions are expressed in astronomical units (AU),
+ * and the velocities in AU/day.
 *
 * The following integer constants may be useful for indexing
 * into the `moon` array: #JM_IO, #JM_EUROPA, #JM_GANYMEDE, #JM_CALLISTO.
 */
 typedef struct
 {
-    astro_vector_t moon[NUM_JUPITER_MOONS];     /**< Jovicentric coordinates of each moon, as described above. */
+    astro_state_vector_t moon[NUM_JUPITER_MOONS];     /**< Jovicentric position and velocity of each moon, as described above. */
 }
 astro_jupiter_moons_t;

@@ -987,6 +1004,7 @@ astro_equatorial_t Astronomy_EquatorFromVector(astro_vector_t vector);
 astro_vector_t Astronomy_VectorFromHorizon(astro_spherical_t sphere, astro_time_t time, astro_refraction_t refraction);
 astro_spherical_t Astronomy_HorizonFromVector(astro_vector_t vector, astro_refraction_t refraction);
 astro_vector_t Astronomy_RotateVector(astro_rotation_t rotation, astro_vector_t vector);
+astro_state_vector_t Astronomy_RotateState(astro_rotation_t rotation, astro_state_vector_t state);

 astro_rotation_t Astronomy_Rotation_EQD_EQJ(astro_time_t time);
 astro_rotation_t Astronomy_Rotation_EQD_ECL(astro_time_t time);