From 2b97530afd8480953b64dca56ff19b9f26fcb148 Mon Sep 17 00:00:00 2001
From: Brandon Rhodes <brandon@rhodesmill.org>
Date: Mon, 3 Feb 2020 08:27:39 -0500
Subject: [PATCH] Fix #80 by adding lunar libration calculation

---
 skyfield/documentation/planetary.rst | 30 ++++++++++++++++++++++++++++
 skyfield/positionlib.py              | 17 ++++++++++------
 2 files changed, 41 insertions(+), 6 deletions(-)

diff --git a/skyfield/documentation/planetary.rst b/skyfield/documentation/planetary.rst
index def0a9ca2..866975efe 100644
--- a/skyfield/documentation/planetary.rst
+++ b/skyfield/documentation/planetary.rst
@@ -115,6 +115,36 @@ or in altitude and azimuth measured against the astronaut’s horizon.
     32deg 27' 09.7" degrees above the horizon
     118deg 12' 55.9" degrees around the horizon from north
 
+Computing lunar libration
+=========================
+
+The Moon’s libration is expressed
+as the latitude and longitude of the Moon location
+that is currently nearest the Earth.
+The convention seems to be that the simple geometric difference
+between the Earth’s and Moon’s positions are used,
+rather than the light-delayed position.
+Thus:
+
+.. testcode::
+
+    p = (earth - moon).at(t)
+    lat, lon, distance = p.frame_latlon(frame)
+    lon_degrees = (lon.degrees - 180.0) % 360.0 - 180.0
+    print('Libration in latitude: {:.3f}'.format(lat.degrees))
+    print('Libration in longitude: {:.3f}'.format(lon_degrees))
+
+.. testoutput::
+
+    Libration in latitude: -6.749
+    Libration in longitude: 1.520
+
+The only subtlety is that the libration longitude
+is not expressed as a number between 0° and 360°,
+as would be more usual for longitude,
+but instead as an offset positive or negative from zero,
+which the above code accomplishes with some quick subtraction and modulo.
+
 Computing a raw rotation matrix
 ===============================
 
diff --git a/skyfield/positionlib.py b/skyfield/positionlib.py
index 929112378..d03172ede 100644
--- a/skyfield/positionlib.py
+++ b/skyfield/positionlib.py
@@ -5,7 +5,7 @@
 from .data.spice import inertial_frames
 from .earthlib import compute_limb_angle, refract, reverse_terra
 from .functions import (
-    dots, from_polar, length_of, rot_x, rot_z, to_polar, _to_array,
+    _mxv, dots, from_polar, length_of, rot_x, rot_z, to_polar, _to_array,
 )
 from .relativity import add_aberration, add_deflection
 from .timelib import Time
@@ -310,11 +310,16 @@ def galactic_latlon(self):
     def frame_xyz(self, frame):
         """Express this position as an (x,y,z) vector in a particular frame."""
         R = frame.rotation_at(self.t)
-        # TODO: before documenting this routine, switch dot() to a real
-        # einsum multiply; and when doing so, make a central routine in
-        # functions.py for it instead of scattering yet more einsums
-        # everywhere.
-        return Distance(au=R.dot(self.position.au))
+        return Distance(au=_mxv(R, self.position.au))
+
+    def frame_latlon(self, frame):
+        """Return as longitude, latitude, and distance in the given frame."""
+        R = frame.rotation_at(self.t)
+        vector = _mxv(R, self.position.au)
+        d, lat, lon = to_polar(vector)
+        return (Angle(radians=lat, signed=True),
+                Angle(radians=lon),
+                Distance(au=d))
 
     # Aliases; maybe someday turn into deprecations with warnings?
     ecliptic_position = ecliptic_xyz