Fill delta and phase angle from rh and solar elongation

CHANGES.rst

@@ -7,6 +7,8 @@ New Features
 sbpy.data
 ^^^^^^^^^
 
+- Added ``Ephem.fill_delta_and_phase`` to calculate observer-target distance and
+  phase angle, given heliocentric distance and solar elongation. [#339]
 - Added ``Orbit.tisserand`` to calculate the Tisserand parameter of small
   body's orbits with respect to planets. [#325]
 - Added ``Orbit.D_criterion`` to evaluate the D-criterion between two sets

docs/sbpy/data/ephem.rst

@@ -2,14 +2,63 @@
 Using Ephem
 ===========
 
-As shown above (:ref:`How to use Data Containers`),
-`~sbpy.data.Ephem` objects can be created on the fly. However,
+As shown previously (:ref:`How to use Data Containers`), `~sbpy.data.Ephem`
+objects can be created on the fly:
+
+    >>> import astropy.units as u
+    >>> from sbpy.data import Ephem
+    >>> eph = Ephem.from_dict({'rh': [2] * u.au,
+    ...                        'solar_elongation': [90] * u.deg})
+    >>> print(eph)
+    <QTable length=1>
+       rh   solar_elongation
+       AU         deg       
+    float64     float64     
+    ------- ----------------
+        2.0             90.0
+
+However, this introduces the possibility of nonsensical values, as we can assign
+arbitrary phase angles and observer-target distances to this object:
+
+    >>> eph['delta'] = 1000 * u.au
+    >>> eph['phase'] = 0 * u.deg
+    >>> eph['comment'] = 'nonsense!'
+    >>> print(eph)
+    <QTable length=1>
+       rh   solar_elongation  delta   phase   comment 
+       AU         deg           AU     deg            
+    float64     float64      float64 float64    str9  
+    ------- ---------------- ------- ------- ---------
+        2.0             90.0  1000.0     0.0 nonsense!
+
+Instead, use the the `~sbpy.data.Ephem.fill_delta_and_phase` method, which will
+calculate the phase angle and observer-target distance from heliocentric
+distance and solar elongation.
+
+    >>> eph.fill_delta_and_phase(overwrite=True)  # allow overwriting delta and phase
+    >>> eph['comment'] = 'valid!'
+    >>> print(eph)
+    <QTable length=1>
+       rh   solar_elongation       delta              phase        comment
+       AU         deg                AU                deg                
+    float64     float64           float64            float64         str9 
+    ------- ---------------- ------------------ ------------------ -------
+        2.0             90.0 1.7320508075688772 30.000000000000004  valid!
+
+
+Remote services
+---------------
+
 `~sbpy.data.Ephem` can also be used to access ephemerides information
 from remote services with a largely uniform API.
 
-For instance, the following few lines will query
-ephemerides for asteroid Ceres on a given date and for the position of
-Mauna Kea Observatory (IAU observatory code ``568``) from the `JPL Horizons service <https://ssd.jpl.nasa.gov/horizons.cgi>`_:
+JPL Horizons (`sbpy.data.Ephem.from_horizons`)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+For instance, the following few lines will query ephemerides for asteroid Ceres
+on a given date and for the position of Mauna Kea Observatory (IAU observatory
+code ``568``) from the `JPL Horizons service
+<https://ssd.jpl.nasa.gov/horizons.cgi>`_:
 
     >>> from sbpy.data import Ephem
     >>> from astropy.time import Time
@@ -135,6 +184,9 @@ example:
     ---------- ------- ------- -------------- ... -------- ------- ---------
        1 Ceres    3.34    0.12              * ... 130.4303  9.2004 2458119.5
 
+Minor Planet Center (`sbpy.data.Ephem.from_mpc`)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
 Offering almost identical functionality, the
 `~sbpy.data.Ephem.from_mpc` method will retrieve ephemerides from the
 `Minor Planet Center <https://minorplanetcenter.net/>`_:
@@ -155,6 +207,9 @@ Offering almost identical functionality, the
             2P 2018-10-26 00:00:00.000 ...          81.0         -56.0
             2P 2018-10-23 00:00:00.000 ...          41.0         -41.0
 
+IMCCE Miriade (`sbpy.data.Ephem.from_miriade`)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
 Finally, `~sbpy.data.Ephem.from_miriade` will retrieve ephemerides
 from the `Miriade ephemeris generator
 <http://vo.imcce.fr/webservices/miriade/>`_ at `IMCCE
@@ -175,7 +230,10 @@ from the `Miriade ephemeris generator
          2P 2458415.5 329.83517041666664 ...    -0.055369    25.253586
          2P 2458416.5 329.76366666666667 ...    -0.051392   25.4700287
          2P 2458417.5  329.6967958333333 ...     -0.04743    25.677518
-
+
+Orbital elements (`sbpy.data.Ephem.from_oo`)
+--------------------------------------------
+
 Ephemerides can also be derived from `~sbpy.data.Orbit` objects using
 `sbpy`'s interface to `pyoorb
 <https://github.com/oorb/oorb/tree/master/python>`_ with the function
@@ -209,7 +267,7 @@ from the Discovery Channel Telescope:
        1 Ceres 239.48722955376766 ... 107.87239964547449  2458704.558989811
        1 Ceres 239.49204656314026 ... 107.88097065658197  2458704.600656478
 
-     
+
 The properties computed by pyoorb and listed in the resulting table
 are defined in the `pyoorb documentation
 <https://github.com/oorb/oorb/tree/master/python>`_. Note that this function requires pyoorb to be installed, which is not a requirement for `sbpy`.

sbpy/data/__init__.py

@@ -209,8 +209,9 @@ class Conf():
          'provenance': ['ephem', 'obs'],
          'dimension': dimensions.angle},
         {'description': 'Solar Elongation Angle',
-         'fieldnames': ['elong', 'solarelong',
-                        'solarelongation', 'elongation', 'Elongation'],
+         'fieldnames': ['elong', 'solarelong', 'solar_elong',
+                        'solar_elongation', 'solarelongation', 'elongation',
+                        'Elongation'],
          'provenance': ['ephem', 'obs'],
          'dimension': dimensions.angle},
         {'description': 'V-band Magnitude',

sbpy/data/core.py

@@ -647,6 +647,14 @@ def __setitem__(self, *args):
         """Refer cls.__setitem__ to self._table"""
         self.table.__setitem__(*args)
 
+    def __contains__(self, field_name):
+        """``True`` if the field name, or its equivalent, is present."""
+        try:
+            self._translate_columns(field_name)
+        except KeyError:
+            return False
+        return True
+
     def _translate_columns(self, target_colnames):
         """Translate target_colnames to the corresponding column names
         present in this object's table. Returns a list of actual column
@@ -658,22 +666,23 @@ def _translate_columns(self, target_colnames):
         if not isinstance(target_colnames, (list, ndarray, tuple)):
             target_colnames = [target_colnames]
 
-        translated_colnames = deepcopy(target_colnames)
-        for idx, colname in enumerate(target_colnames):
+        translated_colnames = []
+        for colname in target_colnames:
             # colname is already a column name in self.table
             if colname in self.field_names:
+                translated_colnames.append(colname)
                 continue
             # colname is an alternative column name
-            elif colname in sum(Conf.fieldnames, []):
-                for alt in Conf.fieldnames[Conf.fieldname_idx[colname]]:
-                    # translation available for colname
-                    if alt in self.field_names:
-                        translated_colnames[idx] = alt
-                        break
-            # colname is unknown, raise a KeyError
-            else:
-                raise KeyError('field {:s} not available.'.format(
-                    colname))
+            elif colname in Conf.fieldname_idx:
+                alternatives = Conf.fieldnames[Conf.fieldname_idx[colname]]
+                column_in_table = set(
+                    self.field_names).intersection(alternatives)
+                if len(column_in_table) > 0:
+                    # an alternative for colname is present in the table
+                    translated_colnames.append(column_in_table.pop())
+                    continue
+            # colname not in table, and cannot be translated, raise a KeyError
+            raise KeyError('field {:s} not available.'.format(colname))
 
         return translated_colnames
 

sbpy/data/ephem.py

@@ -853,3 +853,69 @@ def from_oo(cls, orbit, epochs=None, location='500', scope='full',
         ephem.table.remove_column('MJD')
 
         return ephem
+
+    def fill_delta_and_phase(self, observer_rh=1 * u.au, closest=True,
+                             overwrite=False):
+        """Compute observer-target distance and phase angle.
+
+        Requires fields for heliocentric distance and solar elongation.
+
+
+        Parameters
+        ----------
+        observer_rh : `~astropy.units.Quantity`, optional
+            Sun-observer distance, may be a single value or an array of values
+            of same length as the ephemeris object.
+
+        closest : bool, optional
+            Two solutions exist when the target heliocentric distance is closer
+            to the Sun than the observer.  If ``closest == True``, then the
+            solution with the smallest observer-target distance is returned.
+
+        overwrite : bool, optional
+            Set to ``True`` and, if present, the current observer-target
+            distance and phase angle fields will be overwritten.
+
+        """
+
+        # check for required fields
+        if not all(['rh' in self, 'solar_elongation' in self]):
+            raise ValueError("heliocentric distance and solar elongation are "
+                             "required but not present in this object.")
+
+        # delta or phase already exist?
+        if any(['delta' in self, 'phase' in self]) and not overwrite:
+            raise ValueError("delta and/or phase already present in this "
+                             "object, but overwrite is ``False``.")
+
+        # promote observer_rh to an array, as needed, and match units from rh
+        _observer_rh = np.broadcast_to(observer_rh,
+                                       self['rh'].shape,
+                                       subok=True)
+
+        # solve for phase using law of sines
+        phase = np.arcsin(_observer_rh / self['rh']
+                          * np.sin(self['solar_elongation']))
+
+        # carefully handle interior objects
+        sign = np.ones(len(self['rh']))
+        interior = self['rh'] < _observer_rh
+
+        if closest and any(interior):
+            phase[interior] = 180 * u.deg - phase[interior]
+            sign[interior] = -1
+
+        x = np.sqrt(self['rh']**2 - _observer_rh**2 *
+                    np.sin(self['solar_elongation'])**2)
+        delta = _observer_rh * np.cos(self['solar_elongation']) + sign * x
+
+        # carefully handle opposition
+        opposition = np.isclose(self['solar_elongation'], 180 * u.deg)
+        if any(opposition):
+            phase[opposition] = 0
+            delta[opposition] = (_observer_rh[opposition]
+                                 + self['rh'][opposition])
+
+        # save to data object
+        self['delta'] = delta.to(self['rh'].unit)
+        self['phase'] = phase.to(self['solar_elongation'].unit)

sbpy/data/tests/test_dataclass.py

@@ -1,6 +1,7 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import os
 from collections import OrderedDict
+from tempfile import NamedTemporaryFile
 from sbpy.data import dimensions
 import pytest
 from copy import deepcopy
@@ -417,9 +418,13 @@ def test_translate_columns(monkeypatch):
         {
             'fieldnames': ['zz', 'aa'],
             'dimension': dimensions.length
+        },
+        {
+            'fieldnames': ['heliocentric_distance', 'rh'],
+            'dimension': dimensions.length
         }
     ]
-    new_fieldnames = [['zz', 'aa']]
+    new_fieldnames = [info['fieldnames'] for info in new_fieldnames_info]
     new_fieldname_idx = {}
     for idx, field in enumerate(new_fieldnames):
         for alt in field:
@@ -436,9 +441,18 @@ def test_translate_columns(monkeypatch):
     assert tab._translate_columns(['aa', 'bb', 'cc']) == ['aa', 'bb', 'cc']
     assert tab._translate_columns(['zz', 'bb', 'cc']) == ['aa', 'bb', 'cc']
 
+    # x is not in the table
     with pytest.raises(KeyError):
         tab._translate_columns(['x'])
 
+    # heliocentric distance is not in the table
+    with pytest.raises(KeyError):
+        tab._translate_columns(['rh'])
+
+    # test translations via __contains__
+    assert all(col in tab for col in ['aa', 'bb', 'cc', 'zz'])
+    assert not any(col in tab for col in ['x', 'rh', 'heliocentric_distance'])
+
 
 def test_indexing():
     """make sure that indexing functionality is not compromised through
@@ -550,9 +564,10 @@ def test_io():
                                            ('cc', [7, 8, 9]*u.kg)]))
     tab.meta['test'] = 'stuff'
 
-    tab.to_file('dataclass_table.fits', format='fits', overwrite=True)
-
-    tab2 = DataClass.from_file('dataclass_table.fits', format='fits')
+    with NamedTemporaryFile(suffix='.fits', delete=False) as fits_file:
+        tab.to_file(fits_file, format='fits', overwrite=True)
+    tab2 = DataClass.from_file(fits_file.name, format='fits')
+    os.unlink(fits_file.name)
 
     assert all(tab.table == tab2.table)
     assert tab.meta == {key.lower(): val.lower()