C BaryState: first version, test in progress.

Implemented the C function Astronomy_BaryState().
Used JPL Horizons to generate some test data.
Started work on the C unit test for BaryState,
but it is not yet finished. This is just a good
checkpoint for this work in progress.

generate/barystate/README.md

@@ -0,0 +1,30 @@
+### Data for testing barycentric state calculation.
+
+This [JPL Horizons](https://ssd.jpl.nasa.gov/horizons.cgi) data is for testing calculation of the barycentric state
+(position and velocity vectors) using the Astronomy Engine `BaryState` function.
+There are test files for all the supported bodies.
+The files use varying start time, stop time, and time step values,
+but they are otherwise configured the same.
+
+As an example, here is the batch data for generating `Sun.txt`:
+
+```
+!$$SOF
+COMMAND= '10'
+CENTER= '500@0'
+MAKE_EPHEM= 'YES'
+TABLE_TYPE= 'VECTORS'
+START_TIME= '1900-01-01'
+STOP_TIME= '2100-01-01'
+STEP_SIZE= '20 d'
+OUT_UNITS= 'AU-D'
+REF_PLANE= 'FRAME'
+REF_SYSTEM= 'J2000'
+VECT_CORR= 'NONE'
+VEC_LABELS= 'YES'
+VEC_DELTA_T= 'NO'
+CSV_FORMAT= 'NO'
+OBJ_DATA= 'YES'
+VEC_TABLE= '2'
+!$$EOF
+```

generate/ctest.c

@@ -132,6 +132,7 @@ static int GeoidTest(void);
 static int JupiterMoonsTest(void);
 static int Issue103(void);
 static int AberrationTest(void);
+static int BaryStateTest(void);
 
 typedef int (* unit_test_func_t) (void);
 
@@ -145,6 +146,7 @@ unit_test_t;
 static unit_test_t UnitTests[] =
 {
     {"aberration",              AberrationTest},
+    {"barystate",               BaryStateTest},
     {"check",                   AstroCheck},
     {"constellation",           ConstellationTest},
     {"earth_apsis",             EarthApsis},
@@ -3974,7 +3976,7 @@ static int GeoidTestCase(astro_time_t time, astro_observer_t observer, astro_equ
     }
 
     h_diff = ABS(vobs.height - observer.height);
-    DEBUG("C GeoidTestCase: vobs=(lat=%lf, lon=%lf, height=%lf), lat_diff=%lg, lon_diff=%lg, h_diff=%lg\n", 
+    DEBUG("C GeoidTestCase: vobs=(lat=%lf, lon=%lf, height=%lf), lat_diff=%lg, lon_diff=%lg, h_diff=%lg\n",
         vobs.latitude, vobs.longitude, vobs.height, lat_diff, lon_diff, h_diff);
 
     if (lat_diff > 1.0e-6)
@@ -4319,3 +4321,126 @@ fail:
 }
 
 /*-----------------------------------------------------------------------------------------------------------*/
+
+static int VerifyBaryState(
+    astro_body_t body,
+    const char *filename,
+    int lnum,
+    astro_time_t time,
+    double pos[3],
+    double vel[3])
+{
+    return 0;
+}
+
+static int BaryStateBody(astro_body_t body, const char *filename)
+{
+    int error, lnum, nscanned;
+    int found_begin = 0;
+    int found_end = 0;
+    int count = 0;
+    int part = 0;
+    FILE *infile = NULL;
+    char line[100];
+    astro_time_t time;
+    double jd, pos[3], vel[3];
+
+    infile = fopen(filename, "rt");
+    if (infile == NULL)
+        FAIL("C BaryStateBody: Cannot open input file: %s\n", filename);
+
+    lnum = 0;
+    while (!found_end && ReadLine(line, sizeof(line), infile, filename, lnum))
+    {
+        ++lnum;
+        if (!found_begin)
+        {
+            if (strlen(line) >= 5 && !memcmp(line, "$$SOE", 5))
+                found_begin = 1;
+        }
+        else
+        {
+            /*
+                Input comes in triplets of lines:
+
+                2444249.500000000 = A.D. 1980-Jan-11 00:00:00.0000 TDB
+                 X =-3.314860345089456E-01 Y = 8.463418210972562E-01 Z = 3.667227830514760E-01
+                 VX=-1.642704711077836E-02 VY=-5.494770742558920E-03 VZ=-2.383170237527642E-03
+
+                Track which of these 3 cases we are in using the 'part' variable...
+            */
+
+            switch (part)
+            {
+            case 0:
+                if (strlen(line) >= 5 && !memcmp(line, "$$EOE", 5))
+                {
+                    found_end = 1;
+                }
+                else
+                {
+                    nscanned = sscanf(line, "%lf", &jd);
+                    if (nscanned != 1)
+                        FAIL("C BaryStateBody(%s line %d) ERROR reading Julian date.\n", filename, lnum);
+                    V(jd);
+
+                    /* Convert julian day value to astro_time_t. */
+                    time = Astronomy_TimeFromDays(jd - 2451545.0);
+                }
+                break;
+
+            case 1:
+                nscanned = sscanf(line, " X =%lf Y =%lf Z =%lf", &pos[0], &pos[1], &pos[2]);
+                if (nscanned != 3)
+                    FAIL("C BaryStateBody(%s line %d) ERROR reading position vector.\n", filename, lnum);
+                V(pos[0]);
+                V(pos[1]);
+                V(pos[2]);
+                break;
+
+            case 2:
+                nscanned = sscanf(line, " VX=%lf VY=%lf VZ=%lf", &vel[0], &vel[1], &vel[2]);
+                if (nscanned != 3)
+                    FAIL("C BaryStateBody(%s line %d) ERROR reading velocity vector.\n", filename, lnum);
+                V(vel[0]);
+                V(vel[1]);
+                V(vel[2]);
+                CHECK(VerifyBaryState(body, filename, lnum, time, pos, vel));
+                ++count;
+                break;
+
+            default:
+                FAIL("C BaryStateBody: INTERNAL ERROR : part=%d\n", part);
+            }
+
+            part = (part + 1) % 3;
+        }
+    }
+
+    DEBUG("C BaryStateBody(%s): PASS - Tested %d cases.\n", filename, count);
+    error = 0;
+fail:
+    if (infile != NULL) fclose(infile);
+    return error;
+}
+
+static int BaryStateTest(void)
+{
+    int error;  /* set as a side-effect of CHECK macro */
+
+    CHECK(BaryStateBody(BODY_SUN,     "barystate/Sun.txt"));
+    CHECK(BaryStateBody(BODY_MERCURY, "barystate/Mercury.txt"));
+    CHECK(BaryStateBody(BODY_VENUS,   "barystate/Venus.txt"));
+    CHECK(BaryStateBody(BODY_EARTH,   "barystate/Earth.txt"));
+    CHECK(BaryStateBody(BODY_MARS,    "barystate/Mars.txt"));
+    CHECK(BaryStateBody(BODY_JUPITER, "barystate/Jupiter.txt"));
+    CHECK(BaryStateBody(BODY_SATURN,  "barystate/Saturn.txt"));
+    CHECK(BaryStateBody(BODY_URANUS,  "barystate/Uranus.txt"));
+    CHECK(BaryStateBody(BODY_NEPTUNE, "barystate/Neptune.txt"));
+
+    printf("C BaryStateTest: PASS\n");
+fail:
+    return error;
+}
+
+/*-----------------------------------------------------------------------------------------------------------*/

generate/template/astronomy.c

@@ -2852,6 +2852,99 @@ finished:
     return vector;
 }
 
+static astro_state_vector_t ExportState(body_state_t terse, astro_time_t time)
+{
+    astro_state_vector_t state;
+
+    state.status = ASTRO_SUCCESS;
+    state.x = terse.r.x;
+    state.y = terse.r.y;
+    state.z = terse.r.z;
+    state.vx = terse.v.x;
+    state.vy = terse.v.y;
+    state.vz = terse.v.z;
+    state.t = time;
+
+    return state;
+}
+
+/**
+ * @brief  Calculates barycentric position and velocity vectors for the given body.
+ *
+ * Given a body and a time, calculates the barycentric position and velocity
+ * vectors for the center of that body at that time.
+ * The vectors are expressed in equatorial J2000 coordinates (EQJ).
+ *
+ * @param body
+ *      The celestial body whose barycentric state vector is to be calculated.
+ *      Supported values are `BODY_SUN`, `BODY_SSB`, and all planets except Pluto:
+ *      `BODY_MERCURY`, `BODY_VENUS`, `BODY_EARTH`, `BODY_MARS`, `BODY_JUPITER`,
+ *      `BODY_SATURN`, `BODY_URANUS`, `BODY_NEPTUNE`.
+ * @param time
+ *      The date and time for which to calculate position and velocity.
+ * @return
+ *      A structure that contains barycentric position and velocity vectors.
+ */
+astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time)
+{
+    astro_state_vector_t state;
+    body_state_t bary[5];
+    body_state_t planet;
+
+    if (body == BODY_SSB)
+    {
+        /* Trivial case: the solar system barycenter itself. */
+        state.status = ASTRO_SUCCESS;
+        state.x = state.y = state.z = 0.0;
+        state.vx = state.vy = state.vz = 0.0;
+        state.t = time;
+        return state;
+    }
+
+    /*
+        Find the barycentric position and velocity for the 5 major bodies:
+        bary[0] = Sun
+        bary[1] = Jupiter
+        bary[2] = Saturn
+        bary[3] = Uranus
+        bary[4] = Neptune
+    */
+    MajorBodyBary(bary, time.tt);
+
+    switch (body)
+    {
+    /* If the caller is asking for one of the major bodies, we can immediately return the answer. */
+    case BODY_SUN:      return ExportState(bary[0], time);
+    case BODY_JUPITER:  return ExportState(bary[1], time);
+    case BODY_SATURN:   return ExportState(bary[2], time);
+    case BODY_URANUS:   return ExportState(bary[3], time);
+    case BODY_NEPTUNE:  return ExportState(bary[4], time);
+
+    /* Otherwise, we need to calculate the heliocentric state of the given body */
+    /* and add the Sun's heliocentric position to obtain the body's barycentric state. */
+    /* BarySun + HelioBody = BaryBody */
+    /* Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. */
+    case BODY_MERCURY:
+    case BODY_VENUS:
+    case BODY_EARTH:
+    case BODY_MARS:
+        planet = CalcVsopPosVel(&vsop[body], time.tt);
+        state.x  = bary[0].r.x + planet.r.x;
+        state.y  = bary[0].r.y + planet.r.y;
+        state.z  = bary[0].r.z + planet.r.z;
+        state.vx = bary[0].v.x + planet.v.x;
+        state.vy = bary[0].v.y + planet.v.y;
+        state.vz = bary[0].v.z + planet.v.z;
+        state.t  = time;
+        state.status = ASTRO_SUCCESS;
+        return state;
+
+    /* FIXFIXFIX - add support for BODY_MOON, BODY_EMB, BODY_PLUTO. */
+    default:
+        return StateVecError(ASTRO_INVALID_BODY, time);
+    }
+}
+
 /**
  * @brief   Calculates equatorial coordinates of a celestial body as seen by an observer on the Earth's surface.
  *

source/c/README.md

@@ -259,6 +259,31 @@ This function calculates the angular separation between the given body and the S
 
 
 
+---
+
+<a name="Astronomy_BaryState"></a>
+### Astronomy_BaryState(body, time) &#8658; [`astro_state_vector_t`](#astro_state_vector_t)
+
+**Calculates barycentric position and velocity vectors for the given body.** 
+
+
+
+Given a body and a time, calculates the barycentric position and velocity vectors for the center of that body at that time. The vectors are expressed in equatorial J2000 coordinates (EQJ).
+
+
+
+**Returns:**  A structure that contains barycentric position and velocity vectors. 
+
+
+
+| Type | Parameter | Description |
+| --- | --- | --- |
+| [`astro_body_t`](#astro_body_t) | `body` |  The celestial body whose barycentric state vector is to be calculated. Supported values are `BODY_SUN`, `BODY_SSB`, and all planets except Pluto: `BODY_MERCURY`, `BODY_VENUS`, `BODY_EARTH`, `BODY_MARS`, `BODY_JUPITER`, `BODY_SATURN`, `BODY_URANUS`, `BODY_NEPTUNE`.  | 
+| [`astro_time_t`](#astro_time_t) | `time` |  The date and time for which to calculate position and velocity.  | 
+
+
+
+
 ---
 
 <a name="Astronomy_BodyCode"></a>

source/c/astronomy.c

@@ -4078,6 +4078,99 @@ finished:
     return vector;
 }
 
+static astro_state_vector_t ExportState(body_state_t terse, astro_time_t time)
+{
+    astro_state_vector_t state;
+
+    state.status = ASTRO_SUCCESS;
+    state.x = terse.r.x;
+    state.y = terse.r.y;
+    state.z = terse.r.z;
+    state.vx = terse.v.x;
+    state.vy = terse.v.y;
+    state.vz = terse.v.z;
+    state.t = time;
+
+    return state;
+}
+
+/**
+ * @brief  Calculates barycentric position and velocity vectors for the given body.
+ *
+ * Given a body and a time, calculates the barycentric position and velocity
+ * vectors for the center of that body at that time.
+ * The vectors are expressed in equatorial J2000 coordinates (EQJ).
+ *
+ * @param body
+ *      The celestial body whose barycentric state vector is to be calculated.
+ *      Supported values are `BODY_SUN`, `BODY_SSB`, and all planets except Pluto:
+ *      `BODY_MERCURY`, `BODY_VENUS`, `BODY_EARTH`, `BODY_MARS`, `BODY_JUPITER`,
+ *      `BODY_SATURN`, `BODY_URANUS`, `BODY_NEPTUNE`.
+ * @param time
+ *      The date and time for which to calculate position and velocity.
+ * @return
+ *      A structure that contains barycentric position and velocity vectors.
+ */
+astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time)
+{
+    astro_state_vector_t state;
+    body_state_t bary[5];
+    body_state_t planet;
+
+    if (body == BODY_SSB)
+    {
+        /* Trivial case: the solar system barycenter itself. */
+        state.status = ASTRO_SUCCESS;
+        state.x = state.y = state.z = 0.0;
+        state.vx = state.vy = state.vz = 0.0;
+        state.t = time;
+        return state;
+    }
+
+    /*
+        Find the barycentric position and velocity for the 5 major bodies:
+        bary[0] = Sun
+        bary[1] = Jupiter
+        bary[2] = Saturn
+        bary[3] = Uranus
+        bary[4] = Neptune
+    */
+    MajorBodyBary(bary, time.tt);
+
+    switch (body)
+    {
+    /* If the caller is asking for one of the major bodies, we can immediately return the answer. */
+    case BODY_SUN:      return ExportState(bary[0], time);
+    case BODY_JUPITER:  return ExportState(bary[1], time);
+    case BODY_SATURN:   return ExportState(bary[2], time);
+    case BODY_URANUS:   return ExportState(bary[3], time);
+    case BODY_NEPTUNE:  return ExportState(bary[4], time);
+
+    /* Otherwise, we need to calculate the heliocentric state of the given body */
+    /* and add the Sun's heliocentric position to obtain the body's barycentric state. */
+    /* BarySun + HelioBody = BaryBody */
+    /* Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. */
+    case BODY_MERCURY:
+    case BODY_VENUS:
+    case BODY_EARTH:
+    case BODY_MARS:
+        planet = CalcVsopPosVel(&vsop[body], time.tt);
+        state.x  = bary[0].r.x + planet.r.x;
+        state.y  = bary[0].r.y + planet.r.y;
+        state.z  = bary[0].r.z + planet.r.z;
+        state.vx = bary[0].v.x + planet.v.x;
+        state.vy = bary[0].v.y + planet.v.y;
+        state.vz = bary[0].v.z + planet.v.z;
+        state.t  = time;
+        state.status = ASTRO_SUCCESS;
+        return state;
+
+    /* FIXFIXFIX - add support for BODY_MOON, BODY_EMB, BODY_PLUTO. */
+    default:
+        return StateVecError(ASTRO_INVALID_BODY, time);
+    }
+}
+
 /**
  * @brief   Calculates equatorial coordinates of a celestial body as seen by an observer on the Earth's surface.
  *

source/c/astronomy.h

@@ -962,6 +962,7 @@ astro_func_result_t Astronomy_HelioDistance(astro_body_t body, astro_time_t time
 astro_vector_t Astronomy_HelioVector(astro_body_t body, astro_time_t time);
 astro_vector_t Astronomy_GeoVector(astro_body_t body, astro_time_t time, astro_aberration_t aberration);
 astro_vector_t Astronomy_GeoMoon(astro_time_t time);
+astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time);
 astro_jupiter_moons_t Astronomy_JupiterMoons(astro_time_t time);
 
 astro_equatorial_t Astronomy_Equator(