feat(time): Solar/lunar ephemerides

This adds six methods to `caput.time.Observer` to calculate
(solar/lunar) (transits/risings/settings).

These functions are originally from `ch_util.ephemeris` where
they originally worked only on the `chime` `Observer`.  They've
since been enhanced to use any `Observer`, so there's no reason
now not to have them in the `Observer` object directly, now.

CHANGE FROM CH_UTIL VERSION

There is one major change to these functions from what was in
`ch_util.ephemeris`:  instead of a `diameter` of 0.6 degrees
= 36 arcminutes, these methods now use a `diameter` of 100 arcminutes.

This means these functions now compute astronomical sun/moon rise/set,
which is conventionally taken to occur when the centre of the sun/moon
is 50 arcminutes below the horizon.

The 50 arcminute number comes from a solar/lunar radius of 16 arcminutes
plus 34 arcminutes of atmospheric refraction at the horizon.

This change fixes the 3-minute discrepancy Mateus noticed when these
functions were first ported from PyEphem (which _does_ account for
atmospheric refraction) to skyfield (which does _not_ account for the
refraction), here: chime-experiment/ch_util_private@b505cf0

caput/tests/test_time.py

@@ -443,3 +443,55 @@ def test_rise_set_times(chime, eph):
 
     assert risings == approx(precalc_times[precalc_risings], abs=2)
     assert settings == approx(precalc_times[~precalc_risings], abs=2)
+
+def test_solar_ephemerides(chime):
+    """Test solar ephemerides"""
+
+    dts = datetime(2020, 11, 5)
+    dte = datetime(2020, 11, 7)
+
+    times = chime.solar_transit(dts, dte)
+
+    # Calculated via the old version of `ch_util.ephemeris.solar_transit(dts, dte)``
+    precalc_times = [1604605326.0967, 1604691728.9071]
+    assert times == approx(precalc_times, abs=2)
+
+    # From old version of `ch_util.ephemeris`
+    #
+    # NB: these times are different than the ones used in test_rise_set_times
+    #     above, because there both the sun's diameter and atmospheric refraction
+    #     are taken to be zero.
+    precalc_risings = np.array([1604588047, 1604674544], dtype=np.float64)
+
+    risings = chime.solar_rising(dts, dte)
+    assert risings == approx(precalc_risings, abs=2)
+
+    # From old version of `ch_util.ephemeris`
+    precalc_settings = np.array([1604536261, 1604622568], dtype=np.float64)
+
+    settings = chime.solar_setting(dts, dte)
+    assert settings == approx(precalc_settings, abs=2)
+
+def test_lunar_ephemerides(chime):
+    """Test lunar ephemerides"""
+
+    dts = datetime(2020, 11, 5)
+    dte = datetime(2020, 11, 7)
+
+    times = chime.lunar_transit(dts, dte)
+
+    # Calculated via the old version of `ch_util.ephemeris.lunar_transit(dts, dte)``
+    precalc_times = [1604575728, 1604665306]
+    assert times == approx(precalc_times, abs=2)
+
+    # From old version of `ch_util.ephemeris`, with adjusted diameter
+    precalc_risings = np.array([1604545414, 1604634887], dtype=np.float64)
+
+    risings = chime.lunar_rising(dts, dte)
+    assert risings == approx(precalc_risings, abs=2)
+
+    # From old version of `ch_util.ephemeris`, with adjusted diameter
+    precalc_settings = np.array([1604606169, 1604695470], dtype=np.float64)
+
+    settings = chime.lunar_setting(dts, dte)
+    assert settings == approx(precalc_settings, abs=2)

caput/time.py

@@ -67,6 +67,7 @@
 Other skyfield helper functions:
 
 - :py:meth:`skyfield_star_from_ra_dec`
+- :py:meth:`skyfield_time_to_unix`
 
 Constants
 =========
@@ -503,6 +504,9 @@ def f(t):
     def rise_set_times(self, source, t0, t1=None, step=None, diameter=0.0):
         """Find all times a sources rises or sets in an interval.
 
+        Typical atmospheric refraction at the horizon is 34 arcminutes, but
+        this method does _not_ take that into account.
+
         Parameters
         ----------
         source : skyfield source
@@ -522,7 +526,9 @@ def rise_set_times(self, source, t0, t1=None, step=None, diameter=0.0):
         diameter : float
             The size of the source in degrees. Use this to ensure the whole source is
             below the horizon. Also, if the local horizon is higher (i.e. mountains),
-            this can be set to a negative value to account for this.
+            this can be set to a negative value to account for this.  You may also
+            use this parameter to account for atmospheric refraction, if desired,
+            by adding 68 arcminutes to the nominal diameter.
 
         Returns
         -------
@@ -539,6 +545,9 @@ def rise_set_times(self, source, t0, t1=None, step=None, diameter=0.0):
     def rise_times(self, source, t0, t1=None, step=None, diameter=0.0):
         """Find all times a sources rises in an interval.
 
+        Typical atmospheric refraction at the horizon is 34 arcminutes, but
+        this method does _not_ take that into account.
+
         Parameters
         ----------
         source : skyfield source
@@ -557,7 +566,9 @@ def rise_times(self, source, t0, t1=None, step=None, diameter=0.0):
         diameter : float
             The size of the source in degrees. Use this to ensure the whole source is
             below the horizon. Also, if the local horizon is higher (i.e. mountains),
-            this can be set to a negative value to account for this.
+            this can be set to a negative value to account for this.  You may also
+            use this parameter to account for atmospheric refraction, if desired,
+            by adding 68 arcminutes to the nominal diameter.
 
         Returns
         -------
@@ -571,6 +582,9 @@ def rise_times(self, source, t0, t1=None, step=None, diameter=0.0):
     def set_times(self, source, t0, t1=None, step=None, diameter=0.0):
         """Find all times a sources sets in an interval.
 
+        Typical atmospheric refraction at the horizon is 34 arcminutes, but
+        this method does _not_ take that into account.
+
         Parameters
         ----------
         source : skyfield source
@@ -589,7 +603,9 @@ def set_times(self, source, t0, t1=None, step=None, diameter=0.0):
         diameter : float
             The size of the source in degrees. Use this to ensure the whole source is
             below the horizon. Also, if the local horizon is higher (i.e. mountains),
-            this can be set to a negative value to account for this.
+            this can be set to a negative value to account for this.  You may also
+            use this parameter to account for atmospheric refraction, if desired,
+            by adding 68 arcminutes to the nominal diameter.
 
         Returns
         -------
@@ -600,6 +616,192 @@ def set_times(self, source, t0, t1=None, step=None, diameter=0.0):
             source, t0, t1, step, diameter, skip_rise=True, skip_set=False
         )[0]
 
+    def solar_transit(self, t0, t1=None, step=None, lower=False, return_dec=False):
+        """Find the Solar transits between two times.
+
+        Parameters
+        ----------
+        t0 : float unix time, or datetime
+            The start time to search for. Any type that be converted to a UNIX time
+            by caput.
+        t1 : float unix time, or datetime, optional
+            The end time of the search interval. If not set, this is 1 day after the
+            start time `t0`.
+        step : float or None, optional
+            The initial search step in days. This is used to find the approximate
+            times of transit, and should be set to something less than the spacing
+            between events.  If None is passed, an initial search step of 0.2 days,
+            or else one fifth of the specified interval, is used, whichever is smaller.
+        lower : bool, optional
+            By default this only returns the upper (regular) transit. This will cause
+            lower transits to be returned instead.
+        return_dec : bool, optional
+            If set, also return the declination of the source at transit.
+
+        Returns
+        -------
+        times : np.ndarray
+            Solar transit times as UNIX epoch times.
+        """
+        return self.transit_times(
+            skyfield_wrapper.ephemeris["sun"], t0, t1, step, lower, return_dec
+        )
+
+    def lunar_transit(self, t0, t1=None, step=None, lower=False, return_dec=False):
+        """Find the Lunar transits between two times.
+
+        Parameters
+        ----------
+        t0 : float unix time, or datetime
+            The start time to search for. Any type that be converted to a UNIX time
+            by caput.
+        t1 : float unix time, or datetime, optional
+            The end time of the search interval. If not set, this is 1 day after the
+            start time `t0`.
+        step : float or None, optional
+            The initial search step in days. This is used to find the approximate
+            times of transit, and should be set to something less than the spacing
+            between events.  If None is passed, an initial search step of 0.2 days,
+            or else one fifth of the specified interval, is used, whichever is smaller.
+        lower : bool, optional
+            By default this only returns the upper (regular) transit. This will cause
+            lower transits to be returned instead.
+        return_dec : bool, optional
+            If set, also return the declination of the source at transit.
+
+        Returns
+        -------
+        times : np.ndarray
+            Lunar transit times as UNIX epoch times.
+        """
+        return self.transit_times(
+            skyfield_wrapper.ephemeris["moon"], t0, t1, step, lower, return_dec
+        )
+
+    def solar_setting(self, t0, t1=None, step=None):
+        """Find the Solar settings between two times.
+
+        This method calculates the conventional astronomical sunset, which
+        occurs when the centre of the sun is 50 arcminutes below the horizon.
+        This accounts for a solar diameter of 32 arcminutes, plus 34 arcminutes
+        of atmospheric refraction at the horizon.
+
+        Parameters
+        ----------
+        t0 : float unix time, or datetime
+            The start time to search for. Any type that be converted to a UNIX time
+            by caput.
+        t1 : float unix time, or datetime, optional
+            The end time of the search interval. If not set, this is 1 day after the
+            start time `t0`.
+        step : float or None, optional
+            The initial search step in days. This is used to find the approximate
+            times of setting, and should be set to something less than the spacing
+            between events.  If None is passed, an initial search step of 0.2 days,
+            or else one fifth of the specified interval, is used, whichever is smaller.
+
+        Returns
+        -------
+        times : np.ndarray
+            Solar setting times as UNIX epoch times.
+        """
+        return self.set_times(
+            skyfield_wrapper.ephemeris["sun"], t0, t1, step, diameter=100. / 60,
+        )
+
+    def lunar_setting(self, t0, t1=None, step=None):
+        """Find the Lunar settings between two times.
+
+        This method calculates the conventional astronomical moonset, which
+        occurs when the centre of the moon is 50 arcminutes below the horizon.
+        This accounts for a lunar diameter of 32 arcminutes, plus 34 arcminutes
+        of atmospheric refraction at the horizon.
+
+        Parameters
+        ----------
+        t0 : float unix time, or datetime
+            The start time to search for. Any type that be converted to a UNIX time
+            by caput.
+        t1 : float unix time, or datetime, optional
+            The end time of the search interval. If not set, this is 1 day after the
+            start time `t0`.
+        step : float or None, optional
+            The initial search step in days. This is used to find the approximate
+            times of setting, and should be set to something less than the spacing
+            between events.  If None is passed, an initial search step of 0.2 days,
+            or else one fifth of the specified interval, is used, whichever is smaller.
+
+        Returns
+        -------
+        times : np.ndarray
+            Lunar setting times as UNIX epoch times.
+        """
+        return self.set_times(
+            skyfield_wrapper.ephemeris["moon"], t0, t1, step, diameter=100. / 60,
+        )
+
+    def solar_rising(self, t0, t1=None, step=None):
+        """Find the Solar risings between two times.
+
+        This method calculates the conventional astronomical sunrise, which
+        occurs when the centre of the sun is 50 arcminutes below the horizon.
+        This accounts for a solar diameter of 32 arcminutes, plus 34 arcminutes
+        of atmospheric refraction at the horizon.
+
+        Parameters
+        ----------
+        t0 : float unix time, or datetime
+            The start time to search for. Any type that be converted to a UNIX time
+            by caput.
+        t1 : float unix time, or datetime, optional
+            The end time of the search interval. If not set, this is 1 day after the
+            start time `t0`.
+        step : float or None, optional
+            The initial search step in days. This is used to find the approximate
+            times of rising, and should be set to something less than the spacing
+            between events.  If None is passed, an initial search step of 0.2 days,
+            or else one fifth of the specified interval, is used, whichever is smaller.
+
+        Returns
+        -------
+        times : np.ndarray
+            Solar rising times as UNIX epoch times.
+        """
+        return self.rise_times(
+            skyfield_wrapper.ephemeris["sun"], t0, t1, step, diameter=100. / 60,
+        )
+
+    def lunar_rising(self, t0, t1=None, step=None):
+        """Find the Lunar risings between two times.
+
+        This method calculates the conventional astronomical moonrise, which
+        occurs when the centre of the moon is 50 arcminutes below the horizon.
+        This accounts for a lunar diameter of 32 arcminutes, plus 34 arcminutes
+        of atmospheric refraction at the horizon.
+
+        Parameters
+        ----------
+        t0 : float unix time, or datetime
+            The start time to search for. Any type that be converted to a UNIX time
+            by caput.
+        t1 : float unix time, or datetime, optional
+            The end time of the search interval. If not set, this is 1 day after the
+            start time `t0`.
+        step : float or None, optional
+            The initial search step in days. This is used to find the approximate
+            times of rising, and should be set to something less than the spacing
+            between events.  If None is passed, an initial search step of 0.2 days,
+            or else one fifth of the specified interval, is used, whichever is smaller.
+
+        Returns
+        -------
+        times : np.ndarray
+            Lunar rising times as UNIX epoch times.
+        """
+        return self.rise_times(
+            skyfield_wrapper.ephemeris["moon"], t0, t1, step, diameter=100. / 60,
+        )
+
     def _sr_work(self, source, t0, t1, step, diameter, skip_rise=False, skip_set=False):
         # A work routine factoring out common functionality