Repository: astropy/astroquery

Issue:
Topographic position disallowed in vectors query but would be useful for tracking
With these example definitions
```
DSN_goldstone = '257@399'
SSN_millstone = {'lon':-71.490967, 'lat':42.617442, 'elevation':0.127}
epochs = {'start': '1998-01-23 06:00:00', 'stop': '1998-01-23 12:00:00', 'step': '1m'}
spacecraft='NEAR'
```
the call `Horizons(id=spacecraft, location=DSN_goldstone, epochs=epochs).vectors()` returns topocentric range and range rate representative of Doppler tracking from DSN Goldstone in this case. The equivalent call `Horizons(id=spacecraft, location=SSN_millstone, epochs=epochs).vectors()`, to compare radar range from MIT's Haystack is refused with `ValueError: cannot use topographic position in statevectors query`, solely because of its location argument format.

This seems inconsistent, and perhaps the first query should be also blocked for violating the spirit of the vectors query. However, I feel the Horizons topocentric range and range rate would be a reliable reference for cross-checking radar or telemetry tracking of satellites.

Discussion:
m-stclair: hi @bsipocz -- I've reimplemented these features, added basic test coverage for them, updated other tests as necessary, and updated docstrings. I'm assuming the fast CI failures where the pipeline can't find `tox` are because this is a draft; if not, I have no idea how to interrogate those. What else do you think needs to be done before I mark this as ready for review?
bsipocz: It's a weird and new failure about tox, the draft status of the PR shouldn't affect the CI this way. I'll look into it and fix it separately.

If the content is ready, then do mark it ready for review, and we can ask @mkelley to do a first round of review while I'm dealing with CI.
(Tests can be checked locally, including the remote tests, etc.)
bsipocz: @m-stclair - could you rebase? CI has been fixed, but I cannot step in as a maintainer to do the rebase as the PR was opened from an org fork rather than from yours.

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index e5351c15..82633349 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -100,6 +100,9 @@ jplhorizons
 
 - Optional keyword arguments are now keyword only. [#1802]
 
+- Topocentric coordinates can now be specified for both center and target in observer
+  and vector queries. [#2625]
+
 jplsbdb
 ^^^^^^^
 
diff --git a/astroquery/jplhorizons/__init__.py b/astroquery/jplhorizons/__init__.py
index 83cb0568..5ab76be8 100644
--- a/astroquery/jplhorizons/__init__.py
+++ b/astroquery/jplhorizons/__init__.py
@@ -50,7 +50,10 @@ class Conf(_config.ConfigNamespace):
                    'k2': ('k2', '---'),
                    'phasecoeff': ('phasecoeff', 'mag/deg'),
                    'solar_presence': ('solar_presence', '---'),
-                   'flags': ('flags', '---'),
+                   'lunar_presence': ('lunar_presence', '---'),
+                   'interfering_body': ('interfering_body', '---'),
+                   'illumination_flag': ('illumination_flag', '---'),
+                   'nearside_flag': ('nearside_flag', '---'),
                    'R.A._(ICRF)': ('RA', 'deg'),
                    'DEC_(ICRF)': ('DEC', 'deg'),
                    'R.A.___(ICRF)': ('RA', 'deg'),
diff --git a/astroquery/jplhorizons/core.py b/astroquery/jplhorizons/core.py
index 54ad23ea..7861425e 100644
--- a/astroquery/jplhorizons/core.py
+++ b/astroquery/jplhorizons/core.py
@@ -1,15 +1,15 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
-
 # 1. standard library imports
-from numpy import nan
-from numpy import isnan
-from numpy import ndarray
 from collections import OrderedDict
+from typing import Mapping
 import warnings
 
 # 2. third party imports
 from requests.exceptions import HTTPError
+from numpy import nan
+from numpy import isnan
+from numpy import ndarray
 from astropy.table import Table, Column
 from astropy.io import ascii
 from astropy.time import Time
@@ -51,8 +51,16 @@ class HorizonsClass(BaseQuery):
         Parameters
         ----------
 
-        id : str, required
-            Name, number, or designation of the object to be queried.
+        id : str or dict, required
+            Name, number, or designation of target object. Uses the same codes
+            as JPL Horizons. Arbitrary topocentric coordinates can be added
+            in a dict. The dict has to be of the form
+            {``'lon'``: longitude in deg (East positive, West
+            negative), ``'lat'``: latitude in deg (North positive, South
+            negative), ``'elevation'``: elevation in km above the reference
+            ellipsoid, [``'body'``: Horizons body ID of the central body;
+            optional; if this value is not provided it is assumed that this
+            location is on Earth]}.
 
         location : str or dict, optional
             Observer's location for ephemerides queries or center body name for
@@ -108,9 +116,16 @@ class HorizonsClass(BaseQuery):
         """
 
         super().__init__()
-        self.id = id
-        self.location = location
-
+        # check & format coordinate dictionaries for id and location; simply
+        # treat other values as given
+        if isinstance(id, Mapping):
+            self.id = self._prep_loc_dict(dict(id), "id")
+        else:
+            self.id = id
+        if isinstance(location, Mapping):
+            self.location = self._prep_loc_dict(dict(location), "location")
+        else:
+            self.location = location
         # check for epochs to be dict or list-like; else: make it a list
         if epochs is not None:
             if isinstance(epochs, (list, tuple, ndarray)):
@@ -535,16 +550,22 @@ class HorizonsClass(BaseQuery):
 
         URL = conf.horizons_server
 
-        # check for required information
+        # check for required information and assemble commanddline stub
         if self.id is None:
             raise ValueError("'id' parameter not set. Query aborted.")
+        elif isinstance(self.id, dict):
+            commandline = (
+                f"g:{self.id['lon']},{self.id['lat']},"
+                f"{self.id['elevation']}@{self.id['body']}"
+            )
+        else:
+            commandline = str(self.id)
         if self.location is None:
             self.location = '500@399'
         if self.epochs is None:
             self.epochs = Time.now().jd
+        # expand commandline based on self.id_type
 
-        # assemble commandline based on self.id_type
-        commandline = str(self.id)
         if self.id_type in ['designation', 'name',
                             'asteroid_name', 'comet_name']:
             commandline = ({'designation': 'DES=',
@@ -580,19 +601,7 @@ class HorizonsClass(BaseQuery):
             ('EXTRA_PREC', {True: 'YES', False: 'NO'}[extra_precision])])
 
         if isinstance(self.location, dict):
-            if ('lon' not in self.location or 'lat' not in self.location or
-                    'elevation' not in self.location):
-                raise ValueError(("'location' must contain lon, lat, "
-                                  "elevation"))
-
-            if 'body' not in self.location:
-                self.location['body'] = '399'
-            request_payload['CENTER'] = 'coord@{:s}'.format(
-                str(self.location['body']))
-            request_payload['COORD_TYPE'] = 'GEODETIC'
-            request_payload['SITE_COORD'] = "'{:f},{:f},{:f}'".format(
-                self.location['lon'], self.location['lat'],
-                self.location['elevation'])
+            request_payload = dict(**request_payload, **self._location_to_params(self.location))
         else:
             request_payload['CENTER'] = "'" + str(self.location) + "'"
 
@@ -1032,17 +1041,18 @@ class HorizonsClass(BaseQuery):
 
         URL = conf.horizons_server
 
-        # check for required information
+        # check for required information and assemble commandline stub
         if self.id is None:
             raise ValueError("'id' parameter not set. Query aborted.")
+        elif isinstance(self.id, dict):
+            commandline = "g:{lon},{lat},{elevation}@{body}".format(**self.id)
+        else:
+            commandline = str(self.id)
         if self.location is None:
             self.location = '500@10'
         if self.epochs is None:
             self.epochs = Time.now().jd
-
-        # assemble commandline based on self.id_type
-        commandline = str(self.id)
-
+        # expand commandline based on self.id_type
         if self.id_type in ['designation', 'name',
                             'asteroid_name', 'comet_name']:
             commandline = ({'designation': 'DES=',
@@ -1060,18 +1070,12 @@ class HorizonsClass(BaseQuery):
                 commandline += ' CAP{:s};'.format(closest_apparition)
             if no_fragments:
                 commandline += ' NOFRAG;'
-
-        if isinstance(self.location, dict):
-            raise ValueError(('cannot use topographic position in state'
-                              'vectors query'))
-
-        # configure request_payload for ephemerides query
+        # configure request_payload for vectors query
         request_payload = OrderedDict([
             ('format', 'text'),
             ('EPHEM_TYPE', 'VECTORS'),
             ('OUT_UNITS', 'AU-D'),
             ('COMMAND', '"' + commandline + '"'),
-            ('CENTER', ("'" + str(self.location) + "'")),
             ('CSV_FORMAT', ('"YES"')),
             ('REF_PLANE', {'ecliptic': 'ECLIPTIC',
                            'earth': 'FRAME',
@@ -1086,7 +1090,12 @@ class HorizonsClass(BaseQuery):
             ('VEC_DELTA_T', {True: 'YES', False: 'NO'}[delta_T]),
             ('OBJ_DATA', 'YES')]
         )
-
+        if isinstance(self.location, dict):
+            request_payload = dict(
+                **request_payload, **self._location_to_params(self.location)
+            )
+        else:
+            request_payload['CENTER'] = "'" + str(self.location) + "'"
         # parse self.epochs
         if isinstance(self.epochs, (list, tuple, ndarray)):
             request_payload['TLIST'] = "\n".join([str(epoch) for epoch in
@@ -1132,6 +1141,30 @@ class HorizonsClass(BaseQuery):
         return response
 
     # ---------------------------------- parser functions
+    @staticmethod
+    def _prep_loc_dict(loc_dict, attr_name):
+        """prepare coord specification dict for 'location' or 'id'"""
+        if {'lat', 'lon', 'elevation'} - loc_dict.keys():
+            raise ValueError(
+                f"dict values for '{attr_name}' must contain 'lat', 'lon', "
+                "'elevation' (and optionally 'body')"
+            )
+        if 'body' not in loc_dict:
+            loc_dict['body'] = 399
+        return loc_dict
+
+    @staticmethod
+    def _location_to_params(loc_dict):
+        """translate a 'location' dict to a dict of request parameters"""
+        loc_dict = {
+            "CENTER": f"coord@{loc_dict['body']}",
+            "COORD_TYPE": "GEODETIC",
+            "SITE_COORD": ",".join(
+                str(float(loc_dict[k])) for k in ['lon', 'lat', 'elevation']
+            )
+        }
+        loc_dict["SITE_COORD"] = f"'{loc_dict['SITE_COORD']}'"
+        return loc_dict
 
     def _parse_result(self, response, verbose=None):
         """
@@ -1181,14 +1214,18 @@ class HorizonsClass(BaseQuery):
         H, G = nan, nan
         M1, M2, k1, k2, phcof = nan, nan, nan, nan, nan
         headerline = []
+        centername = ''
         for idx, line in enumerate(src):
             # read in ephemerides header line; replace some field names
             if (self.query_type == 'ephemerides' and
                     "Date__(UT)__HR:MN" in line):
                 headerline = str(line).split(',')
                 headerline[2] = 'solar_presence'
-                headerline[3] = 'flags'
+                headerline[3] = "lunar_presence" if "Earth" in centername else "interfering_body"
                 headerline[-1] = '_dump'
+                if isinstance(self.id, dict) or str(self.id).startswith('g:'):
+                    headerline[4] = 'nearside_flag'
+                    headerline[5] = 'illumination_flag'
             # read in elements header line
             elif (self.query_type == 'elements' and
                   "JDTDB," in line):
@@ -1208,6 +1245,9 @@ class HorizonsClass(BaseQuery):
             # read in targetname
             if "Target body name" in line:
                 targetname = line[18:50].strip()
+            # read in center body name
+            if "Center body name" in line:
+                centername = line[18:50].strip()
             # read in H and G (if available)
             if "rotational period in hours)" in line:
                 HGline = src[idx + 2].split('=')
diff --git a/docs/jplhorizons/jplhorizons.rst b/docs/jplhorizons/jplhorizons.rst
index 4c8fc78b..9a6feba1 100644
--- a/docs/jplhorizons/jplhorizons.rst
+++ b/docs/jplhorizons/jplhorizons.rst
@@ -40,26 +40,59 @@ In order to query information for a specific Solar System body, a
    >>> print(obj)
    JPLHorizons instance "Ceres"; location=568, epochs=[2458133.33546], id_type=None
 
-``id`` refers to the target identifier and is mandatory; the exact
-string will be used in the query to the Horizons system.
-
-``location`` means either the observer's location (e.g., Horizons ephemerides
-query) or the body relative to which orbital elements are provided (e.g.,
-Horizons orbital elements or vectors query); the same codes as `used by Horizons
-<https://ssd.jpl.nasa.gov/horizons/manual.html#center>`_ are used here, which
-includes `MPC Observatory codes`_. The default is ``location=None``, which uses
-a geocentric location for ephemerides queries and the Sun as central body for
-orbital elements and state vector queries. User-defined topocentric locations
-for ephemerides queries can be provided, too, in the form of a dictionary. The
-dictionary has to be formatted as follows: {``'lon'``: longitude in degrees
-(East positive, West negative), ``'lat'``: latitude in degrees (North positive,
-South negative), ``'elevation'``: elevation in km above the reference
-ellipsoid}. In addition, ``'body'`` can be set to the Horizons body ID of the
-central body if different from Earth; by default, it is assumed that this
-location is on Earth if it has not been specifically set. The following example
-uses the coordinates of the `Statue of Liberty
+``id`` refers to the target object and is mandatory. ````str`` and ``int`` values
+are valid for all query types. ``Mapping`` (e.g. ``dict``) values are valid for
+observer (``ephemerides``) and vectors queries only. ``str`` or ``int`` values will be
+passed directly to Horizons. See the description of the ``id_type`` argument below
+for how Horizons interprets these values. See the paragraph below the description
+of the ``location`` argument for valid ``dict`` formatting.
+
+``location`` refers to the coordinate center for the ephemeris, which has
+slightly different physical interpretations depending on the query type:
+observer (``ephemerides``) queries: observer location
+vectors queries: coordinate origin for vectors
+elements queries: relative body for orbital elements
+
+``str`` and ``int`` values are valid for all query types. ``Mapping``
+(e.g. ``dict``) values are valid for observer (``ephemerides``) and vectors queries only. ``str`` or ``int``
+arguments will be passed directly to Horizons. See `this section of the Horizons
+manual`_ for how Horizons interprets coordinate center codes; also note that,
+unlike ``id``, these include (most) `MPC Observatory codes`_. See below for valid
+``dict`` formatting. The default is ``location=None``, which uses Earth body center
+for observer queries, and Sun body center for orbital elements and vectors queries.
+
+``dict``-like arguments to ``id`` or ``location`` define a topocentric location
+relative to a major body. Note that this is not possible for elements queries,
+and will only work for bodies with defined rotation models (Horizons does not
+have rotation models for many small or recently-discovered natural satellites).
+The dictionary has to be formatted as follows: {``'lon'``: longitude in degrees,
+``'lat'``: latitude in degrees (North positive, South negative), ``'elevation'``:
+elevation in km above the reference ellipsoid, ``body``: Horizons ID of center body,
+optional; default Earth}.
+
+Horizons always interprets longitude values as eastward. However, there are two
+major gotchas in this:
+1. For most prograde rotators, which is to say most major bodies, Horizons
+interprets west-longitude as positive and east-longitude as negative. However, values
+must still be entered in east-longitude, which means they must be negative; Horizons
+will raise an error if given any positive longitude value for such bodies. Instead enter
+the west-longitude - 360. For instance, a site on Mars (id code 499) at 30 degrees
+longitude, 30 degrees latitude, 0 km elevation should be specified as
+``{'body': 499, 'elevation': 0, 'lon': -330, 'lat': 30}``.
+2. This does not apply to the Earth, Moon, and Sun. Although they are prograde,
+Horizons interprets east-longitude as positive and west-longitude as negative for these
+bodies.
+
+Here is a complete list of retrograde major bodies in Horizons:  Venus (299), Arial (701),
+Umbriel (702), Titania (703), Oberon (704), Miranda (705), Cordelia (706), Ophelia (707),
+Bianca (708), Cressida (709), Desdemona (710), Juliet (711), Portia (712), Rosalind (713),
+Belinda (714), Puck (715), Uranus (799), Triton (801). All other major bodies are prograde.
+
+Two examples of usage for specified topocentric coordinates follow.
+
+1. This observer (``ephemerides``) query uses the coordinates of the `Statue of Liberty
 <https://www.google.com/maps/place/Statue+of+Liberty+National+Monument/@40.6892534,-74.0466891,17z/data=!3m1!4b1!4m5!3m4!1s0x89c25090129c363d:0x40c6a5770d25022b!8m2!3d40.6892494!4d-74.0445004>`_
-as the observer's location:
+as the observer's location, and Ceres as the target:
 
 .. doctest-remote-data::
 
@@ -72,6 +105,24 @@ as the observer's location:
     >>> print(obj)
     JPLHorizons instance "Ceres"; location={'lon': -74.0466891, 'lat': 40.6892534, 'elevation': 0.093}, epochs=[2458133.33546], id_type=None
 
+2. Specifying topocentric coordinates for both location and observer is often
+useful when performing geometric calculations for artificial satellites without
+completely-specified ephemeris data. For instance, published reduced data for the
+lunar satellite Chang'e-2 include orbital height and lat/lon. Although there is no published
+ephemeris for Chang'e-2, Horizons (combined with the fact that Chang'e-2 looked nadir),
+can be used to compute vectors from Chang'e-2 to specific points on the lunar surface.
+Here is an example of using ``jplhorizons`` to find the distance from Chang'e-2
+at a particular point in time to the center of the crater Double:
+
+.. doctest-remote-data::
+
+    >>> ce_2 = {'lon': 23.522, 'lat': 0.637, 'elevation': 181.2, 'body': 301}
+    >>> double = {'lon': 23.47, 'lat': 0.67, 'elevation': 0, 'body': 301}
+    >>> obj = Horizons(id=double, location=ce_2, epochs=2454483.84247)
+    >>> vecs = obj.vectors()
+    >>> distance_km = (vecs['x'] ** 2 + vecs['y'] ** 2 + vecs['z'] ** 2) ** 0.5 * 1.496e8
+    >>> print(f"{distance_km.value.data[0]:.3f}")
+    181.213
 
 ``epochs`` is either a scalar or list of Julian dates (floats or strings) in the
 case of discrete epochs, or, in the case of a range of epochs, a dictionary that
@@ -609,3 +660,4 @@ Reference/API
 .. _astropy units: http://docs.astropy.org/en/stable/units/index.html
 .. _Definition of Observer Table Quantities: https://ssd.jpl.nasa.gov/horizons/manual.html#observer-table
 .. _Horizons documentation: https://ssd.jpl.nasa.gov/horizons/manual.html#observer-table
+.. _this section of the Horizons manual: <https://ssd.jpl.nasa.gov/horizons/manual.html#center>
\ No newline at end of file


Test patch:
diff --git a/astroquery/jplhorizons/tests/test_jplhorizons.py b/astroquery/jplhorizons/tests/test_jplhorizons.py
index b8fbb4fe..129f13f0 100644
--- a/astroquery/jplhorizons/tests/test_jplhorizons.py
+++ b/astroquery/jplhorizons/tests/test_jplhorizons.py
@@ -83,7 +83,7 @@ def test_ephemerides_query(patch_request):
     assert res['targetname'] == "1 Ceres (A801 AA)"
     assert res['datetime_str'] == "2000-Jan-01 00:00:00.000"
     assert res['solar_presence'] == ""
-    assert res['flags'] == ""
+    assert res['lunar_presence'] == ""
     assert res['elongFlag'] == '/L'
     assert res['airmass'] == 999
 
@@ -256,22 +256,21 @@ def test_vectors_query_payload():
     res = jplhorizons.Horizons(id='Ceres', location='500@10',
                                epochs=2451544.5).vectors(
                                    get_query_payload=True)
-
     assert res == OrderedDict([
-        ('format', 'text'),
-        ('EPHEM_TYPE', 'VECTORS'),
-        ('OUT_UNITS', 'AU-D'),
-        ('COMMAND', '"Ceres"'),
-        ('CENTER', "'500@10'"),
-        ('CSV_FORMAT', '"YES"'),
-        ('REF_PLANE', 'ECLIPTIC'),
-        ('REF_SYSTEM', 'ICRF'),
-        ('TP_TYPE', 'ABSOLUTE'),
-        ('VEC_LABELS', 'YES'),
-        ('VEC_CORR', '"NONE"'),
-        ('VEC_DELTA_T', 'NO'),
-        ('OBJ_DATA', 'YES'),
-        ('TLIST', '2451544.5')])
+                ('format', 'text'),
+                ('EPHEM_TYPE', 'VECTORS'),
+                ('OUT_UNITS', 'AU-D'),
+                ('COMMAND', '"Ceres"'),
+                ('CSV_FORMAT', '"YES"'),
+                ('REF_PLANE', 'ECLIPTIC'),
+                ('REF_SYSTEM', 'ICRF'),
+                ('TP_TYPE', 'ABSOLUTE'),
+                ('VEC_LABELS', 'YES'),
+                ('VEC_CORR', '"NONE"'),
+                ('VEC_DELTA_T', 'NO'),
+                ('OBJ_DATA', 'YES'),
+                ('CENTER', "'500@10'"),
+                ('TLIST', '2451544.5')])
 
 
 def test_no_H(patch_request):
diff --git a/astroquery/jplhorizons/tests/test_jplhorizons_remote.py b/astroquery/jplhorizons/tests/test_jplhorizons_remote.py
index 4f69be68..07c8d378 100644
--- a/astroquery/jplhorizons/tests/test_jplhorizons_remote.py
+++ b/astroquery/jplhorizons/tests/test_jplhorizons_remote.py
@@ -1,10 +1,12 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
-
-import pytest
 from numpy.ma import is_masked
+import numpy as np
+import pytest
 
+from astropy.coordinates import spherical_to_cartesian
 from astropy.tests.helper import assert_quantity_allclose
+import astropy.units as u
 from astropy.utils.exceptions import AstropyDeprecationWarning
 
 from ... import jplhorizons
@@ -23,7 +25,7 @@ class TestHorizonsClass:
         assert res['targetname'] == "1 Ceres (A801 AA)"
         assert res['datetime_str'] == "2000-Jan-01 00:00:00.000"
         assert res['solar_presence'] == ""
-        assert res['flags'] == ""
+        assert res['lunar_presence'] == ""
         assert res['elongFlag'] == '/L'
         assert res['airmass'] == 999
 
@@ -75,7 +77,7 @@ class TestHorizonsClass:
         assert res['targetname'] == "1P/Halley"
         assert res['datetime_str'] == "2080-Jan-11 09:00"
         assert res['solar_presence'] == ""
-        assert res['flags'] == "m"
+        assert res['lunar_presence'] == "m"
         assert res['elongFlag'] == '/L'
 
         for value in ['H', 'G']:
@@ -98,7 +100,7 @@ class TestHorizonsClass:
         assert res['targetname'] == "73P/Schwassmann-Wachmann 3"
         assert res['datetime_str'] == "2080-Jan-01 00:00"
         assert res['solar_presence'] == "*"
-        assert res['flags'] == "m"
+        assert res['lunar_presence'] == "m"
         assert res['elongFlag'] == '/L'
 
         for value in ['H', 'G']:
@@ -123,7 +125,7 @@ class TestHorizonsClass:
         assert res['targetname'] == "167P/CINEOS"
         assert res['datetime_str'] == "2080-Jan-01 00:00"
         assert res['solar_presence'] == "*"
-        assert res['flags'] == "m"
+        assert res['lunar_presence'] == "m"
         assert res['elongFlag'] == '/T'
 
         for value in ['H', 'G', 'M1', 'k1']:
@@ -150,7 +152,7 @@ class TestHorizonsClass:
         assert res['targetname'] == "12P/Pons-Brooks"
         assert res['datetime_str'] == "2080-Jan-01 00:00"
         assert res['solar_presence'] == "*"
-        assert res['flags'] == "m"
+        assert res['lunar_presence'] == "m"
         assert res['elongFlag'] == '/L'
 
         for value in ['H', 'G', 'phasecoeff']:
@@ -405,3 +407,35 @@ class TestHorizonsClass:
         vec_highprec = obj.ephemerides(extra_precision=True)
 
         assert (vec_simple['RA'][0]-vec_highprec['RA'][0]) > 1e-7
+
+    def test_geodetic_queries(self):
+        """
+        black-box test for observer and vectors queries with geodetic
+        coordinates. checks spatial sensibility.
+        """
+        phobos = {'body': 401, 'lon': -30, 'lat': -20, 'elevation': 0}
+        deimos = {'body': 402, 'lon': -10, 'lat': -40, 'elevation': 0}
+        deimos_phobos = jplhorizons.Horizons(phobos, location=deimos, epochs=2.4e6)
+        phobos_deimos = jplhorizons.Horizons(deimos, location=phobos, epochs=2.4e6)
+        pd_eph, dp_eph = phobos_deimos.ephemerides(), deimos_phobos.ephemerides()
+        dp_xyz = spherical_to_cartesian(
+            dp_eph['delta'], dp_eph['DEC'], dp_eph['RA']
+        )
+        pd_xyz = spherical_to_cartesian(
+            pd_eph['delta'], pd_eph['DEC'], pd_eph['RA']
+        )
+        elementwise = [(dp_el + pd_el) for dp_el, pd_el in zip(dp_xyz, pd_xyz)]
+        eph_offset = (sum([off ** 2 for off in elementwise]) ** 0.5).to(u.km)
+        # horizons can do better than this, but we'd have to go to a little
+        # more trouble than is necessary for a software test...
+        assert np.isclose(eph_offset.value, 2.558895)
+        # ...and vectors queries are really what you're meant to use for
+        # this sort of thing.
+        pd_vec, dp_vec = phobos_deimos.vectors(), deimos_phobos.vectors()
+        vec_offset = np.sum(
+            (
+                pd_vec.as_array(names=('x', 'y', 'z')).view('f8')
+                + dp_vec.as_array(names=('x', 'y', 'z')).view('f8')
+            ) ** 2
+        )
+        assert np.isclose(vec_offset, 0)