C Libration: Include Moon's position in the return value.

Because I have to perform the expensive calculation to find
the Moon's ecliptic coordinates, I might as well return them
to the caller. This could help reduce calculation overhead
for some uses, and doesn't add any significant cost.

generate/ctest.c

@@ -4626,8 +4626,8 @@ static int Libration(const char *filename, int *ndata, double *var_lon, double *
     double phase, age, diam, dist, ra, dec, slon, slat, elon, elat, axisa;
     astro_time_t time;
     astro_libration_t lib;
-    double diff_elon, diff_elat;
-    double max_diff_elon = 0.0, max_diff_elat = 0.0;
+    double diff_elon, diff_elat, diff_distance;
+    double max_diff_elon = 0.0, max_diff_elat = 0.0, max_diff_distance = 0.0;
 
     infile = fopen(filename, "rt");
     if (infile == NULL)
@@ -4659,7 +4659,7 @@ static int Libration(const char *filename, int *ndata, double *var_lon, double *
                 FAIL("Libration(%s line %d): invalid month symbol '%s'\n", filename, lnum, mtext);
 
             time = Astronomy_MakeTime(year, month, day, hour, minute, 0.0);
-            lib = Astronomy_Libration(&time);
+            lib = Astronomy_Libration(time);
 
             diff_elon = 60.0 * ABS(lib.elon - elon);
             if (diff_elon > max_diff_elon)
@@ -4669,12 +4669,19 @@ static int Libration(const char *filename, int *ndata, double *var_lon, double *
             if (diff_elat > max_diff_elat)
                 max_diff_elat = diff_elat;
 
+            diff_distance = ABS(lib.dist_km - dist);
+            if (diff_distance > max_diff_distance)
+                max_diff_distance = diff_distance;
+
             if (diff_elon > 0.130)
                 FAIL("C Libration(%s line %d): EXCESSIVE diff_elon = %0.4lf arcmin\n", filename, lnum, diff_elon);
 
             if (diff_elat > 1.666)
                 FAIL("C Libration(%s line %d): EXCESSIVE diff_elat = %0.4lf arcmin\n", filename, lnum, diff_elat);
 
+            if (diff_distance > 53.9)
+                FAIL("C Libration(%s line %d): EXCESSIVE diff_distance = %0.3lf km\n", filename, lnum, diff_distance);
+
             /* Update sum-of-squared-errors. */
             *var_lon += diff_elon * diff_elon;
             *var_lat += diff_elat * diff_elat;
@@ -4683,8 +4690,8 @@ static int Libration(const char *filename, int *ndata, double *var_lon, double *
         }
     }
 
-    printf("C Libration(%s): PASS (%d test cases, max_diff_elon = %0.4lf arcmin, max_diff_elat = %0.4lf arcmin)\n", 
-        filename, count, max_diff_elon, max_diff_elat);
+    printf("C Libration(%s): PASS (%d test cases, max_diff_elon = %0.4lf arcmin, max_diff_elat = %0.4lf arcmin, max_diff_distance = %0.3lf km)\n", 
+        filename, count, max_diff_elon, max_diff_elat, max_diff_distance);
 
     *ndata += count;
     error = 0;

generate/template/astronomy.c

@@ -1877,25 +1877,33 @@ astro_vector_t Astronomy_GeoMoon(astro_time_t time)
  * of the Moon's fixed rotation rate, compared to its variable angular speed
  * of orbit around the Earth.
  *
+ * This function calculates a pair of perpendicular libration angles,
+ * one representing rotation of the Moon in eclitpic longitude `elon`, the other
+ * in ecliptic latitude `elat`, both relative to the Moon's mean Earth-facing position.
+ *
+ * This function also returns the geocentric position of the Moon
+ * expressed in ecliptic longitude `mlon`, ecliptic latitude `mlat`, and
+ * distance `dist_km` between the centers of the Earth and Moon expressed in kilometers.
+ *
  * @param time  The date and time for which to calculate libration angles.
- * @return The Moon's libration in latitude and longitude as seen from the Earth.
+ * @return The Moon's ecliptic position and libration angles as seen from the Earth.
  */
-astro_libration_t Astronomy_Libration(astro_time_t *time)
+astro_libration_t Astronomy_Libration(astro_time_t time)
 {
     astro_libration_t lib;
-    double geo_eclip_lon, geo_eclip_lat, distance_au;
     double t, t2, t3, t4;
     double f, omega, w, a, ldash, ldash2, bdash, bdash2;
     double k1, k2, m, mdash, d, e, rho, sigma, tau;
     double I = DEG2RAD * 1.54242;
 
-    CalcMoon(time->tt / 36525.0, &geo_eclip_lon, &geo_eclip_lat, &distance_au);
-
-    t = time->tt / 36525.0;
+    t = time.tt / 36525.0;
     t2 = t * t;
     t3 = t2 * t;
     t4 = t2 * t2;
 
+    CalcMoon(t, &lib.mlon, &lib.mlat, &lib.dist_km);
+    lib.dist_km *= KM_PER_AU;
+
     /* Moon's argument of latitude in radians. */
     f = DEG2RAD * NormalizeLongitude(93.2720950 + 483202.0175233*t - 0.0036539*t2 - t3/3526000 + t4/863310000);
 
@@ -1915,10 +1923,10 @@ astro_libration_t Astronomy_Libration(astro_time_t *time)
     e = 1.0 - 0.002516*t - 0.0000074*t2;
 
     /* Optical librations */
-    w = geo_eclip_lon - omega;
-    a = atan2(sin(w)*cos(geo_eclip_lat)*cos(I) - sin(geo_eclip_lat)*sin(I), cos(w)*cos(geo_eclip_lat));
+    w = lib.mlon - omega;
+    a = atan2(sin(w)*cos(lib.mlat)*cos(I) - sin(lib.mlat)*sin(I), cos(w)*cos(lib.mlat));
     ldash = LongitudeOffset(RAD2DEG * (a - f));
-    bdash = asin(-sin(w)*cos(geo_eclip_lat)*sin(I) - sin(geo_eclip_lat)*cos(I));
+    bdash = asin(-sin(w)*cos(lib.mlat)*sin(I) - sin(lib.mlat)*cos(I));
 
     /* Physical librations */
     k1 = DEG2RAD*(119.75 + 131.849*t);