diff --git a/skyfield/earthlib.py b/skyfield/earthlib.py
index aa2fb6ca..30dd4cbe 100644
--- a/skyfield/earthlib.py
+++ b/skyfield/earthlib.py
@@ -1,7 +1,7 @@
 """Formulae for specific earth behaviors and effects."""
 
 from numpy import (abs, arcsin, arccos, arctan2, array, clip, cos,
-                   minimum, pi, sin, sqrt, tan, where, zeros_like)
+                   minimum, pi, sin, sqrt, tan, where, zeros_like, dot)
 
 from .constants import (AU_M, ANGVEL, DAY_S, DEG2RAD, ERAD,
                         IERS_2010_INVERSE_EARTH_FLATTENING, RAD2DEG, T0, tau)
@@ -91,7 +91,7 @@ def compute_limb_angle(position_au, observer_au):
 
     # Compute zenith distance of observed object.
 
-    coszd = dots(position_au, observer_au) / (disobj * disobs)
+    coszd = dot(observer_au, position_au) / (disobj * disobs)
     coszd = clip(coszd, -1.0, 1.0)
     zdobj = arccos(coszd)
 
diff --git a/skyfield/relativity.py b/skyfield/relativity.py
index d1320015..77aa30fc 100644
--- a/skyfield/relativity.py
+++ b/skyfield/relativity.py
@@ -1,4 +1,4 @@
-from numpy import abs, einsum, sqrt, where
+from numpy import abs, einsum, sqrt, where, diagonal
 
 from .constants import C, AU_M, C_AUDAY, GS
 from .functions import _AVOID_DIVIDE_BY_ZERO, dots, length_of
@@ -110,6 +110,9 @@ def light_time_difference(position, observer_position):
     # Light-time returned is the projection of vector 'pos_obs' onto the
     # unit vector 'u1' (formed from 'pos1'), divided by the speed of light.
 
+    if observer_position.shape == (3,3):
+        observer_position = diagonal(observer_position)
+
     diflt = einsum('a...,a...', u1, observer_position) / C_AUDAY
     return diflt
 
diff --git a/skyfield/vectorlib.py b/skyfield/vectorlib.py
index a571eb5e..2986f8e1 100644
--- a/skyfield/vectorlib.py
+++ b/skyfield/vectorlib.py
@@ -5,7 +5,7 @@
 from .constants import C_AUDAY
 from .descriptorlib import reify
 from .errors import DeprecationError
-from .functions import length_of
+from .functions import length_of, _reconcile
 from .positionlib import build_position
 from .timelib import Time
 
@@ -215,8 +215,10 @@ def _at(self, t):
             p2, v2, _, message = vf._at(t)
             if vf.center == 399:
                 gcrs_position = -p
-            p += p2
-            v += v2
+            p, p2 = _reconcile(p, p2)
+            p = p + p2
+            v, v2 = _reconcile(v, v2)
+            v = v + v2
         if vfs[0].center == 0 and vf.center == 399:
             gcrs_position = p2
         return p, v, gcrs_position, message
@@ -241,6 +243,7 @@ def _correct_for_light_travel_time(observer, target):
 
     tposition, tvelocity, gcrs_position, message = target._at(t)
 
+    tposition, cposition = _reconcile(tposition, cposition)
     distance = length_of(tposition - cposition)
     light_time0 = 0.0
     for i in range(10):
@@ -259,6 +262,7 @@ def _correct_for_light_travel_time(observer, target):
         light_time0 = light_time
     else:
         raise ValueError('light-travel time failed to converge')
+    tvelocity, cvelocity = _reconcile(tvelocity, cvelocity)
     return tposition - cposition, tvelocity - cvelocity, t, light_time
 
 def _jpl_name(target):