diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 8e9086432b..2a232bc3fc 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -55,6 +55,7 @@ Additional functions for quantities closely related to solar position.
    solarposition.calc_time
    solarposition.pyephem_earthsun_distance
    solarposition.nrel_earthsun_distance
+   solarposition.rise_set_transit_ephem
    spa.calculate_deltat
 
 The spa module contains the implementation of the built-in NREL SPA
diff --git a/docs/sphinx/source/whatsnew/v0.6.1.rst b/docs/sphinx/source/whatsnew/v0.6.1.rst
index 946490e571..e9cc42e51c 100644
--- a/docs/sphinx/source/whatsnew/v0.6.1.rst
+++ b/docs/sphinx/source/whatsnew/v0.6.1.rst
@@ -16,12 +16,15 @@ API Changes
 * Deprecated ``tmy``, ``tmy.readtmy2`` and ``tmy.readtmy3``;
   they will be removed in v0.7. Use the new :py:func:`pvlib.iotools.read_tmy2`
   and :py:func:`pvlib.iotools.read_tmy3` instead. (:issue:`261`)
+* Added kwarg `horizon` to :func:`~pvlib.solarposition.pyephem` and :func:`~pvlib.solarposition.calc_time` with default value `'+0:00'`
 
 
 Enhancements
 ~~~~~~~~~~~~
+* :func:`~pvlib.solarposition.rise_set_transit_ephem` returns sunrise, sunset and transit times using pyephem
 * Created :py:func:`pvlib.iotools.read_srml` and :py:func:`pvlib.iotools.read_srml_month_from_solardat`
   to read University of Oregon Solar Radiation Monitoring Laboratory data. (:issue:`589`)
+  
 
 Bug fixes
 ~~~~~~~~~
@@ -35,4 +38,5 @@ Testing
 Contributors
 ~~~~~~~~~~~~
 * Will Holmgren (:ghuser:`wholmgren`)
+* Cliff Hansen (:ghuser:`cwhanse`)
 * Leland Boeman (:ghuser:`lboeman`)
diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py
index 7f3ed3f8fb..b8b2e20274 100644
--- a/pvlib/solarposition.py
+++ b/pvlib/solarposition.py
@@ -6,6 +6,7 @@
 # Rob Andrews (@Calama-Consulting), Calama Consulting, 2014
 # Will Holmgren (@wholmgren), University of Arizona, 2014
 # Tony Lorenzo (@alorenzo175), University of Arizona, 2015
+# Cliff hansen (@cwhanse), Sandia National Laboratories, 2018
 
 from __future__ import division
 import os
@@ -438,7 +439,26 @@ def get_sun_rise_set_transit(time, latitude, longitude, how='numpy',
     return result
 
 
-def _ephem_setup(latitude, longitude, altitude, pressure, temperature):
+def _ephem_convert_to_seconds_and_microseconds(date):
+    # utility from unreleased PyEphem 3.6.7.1
+    """Converts a PyEphem date into seconds"""
+    microseconds = int(round(24 * 60 * 60 * 1000000 * date))
+    seconds, microseconds = divmod(microseconds, 1000000)
+    seconds -= 2209032000  # difference between epoch 1900 and epoch 1970
+    return seconds, microseconds
+
+
+def _ephem_to_timezone(date, tzinfo):
+    # utility from unreleased PyEphem 3.6.7.1
+    """"Convert a PyEphem Date into a timezone aware python datetime"""
+    seconds, microseconds = _ephem_convert_to_seconds_and_microseconds(date)
+    date = dt.datetime.fromtimestamp(seconds, tzinfo)
+    date = date.replace(microsecond=microseconds)
+    return date
+
+
+def _ephem_setup(latitude, longitude, altitude, pressure, temperature,
+                 horizon):
     import ephem
     # initialize a PyEphem observer
     obs = ephem.Observer()
@@ -447,14 +467,97 @@ def _ephem_setup(latitude, longitude, altitude, pressure, temperature):
     obs.elevation = altitude
     obs.pressure = pressure / 100.  # convert to mBar
     obs.temp = temperature
+    obs.horizon = horizon
 
     # the PyEphem sun
     sun = ephem.Sun()
     return obs, sun
 
 
+def rise_set_transit_ephem(time, latitude, longitude,
+                           next_or_previous='next',
+                           altitude=0,
+                           pressure=101325, temperature=12, horizon='0:00'):
+    """
+    Calculate the next sunrise and sunset times using the PyEphem package.
+
+    Parameters
+    ----------
+    time : pandas.DatetimeIndex
+        Must be localized
+    latitude : float
+        positive is north of 0
+    longitude : float
+        positive is east of 0
+    next_or_previous : str
+        'next' or 'previous' sunrise and sunset relative to time
+    altitude : float, default 0
+        distance above sea level in meters.
+    pressure : int or float, optional, default 101325
+        air pressure in Pascals.
+    temperature : int or float, optional, default 12
+        air temperature in degrees C.
+    horizon : string, format +/-X:YY
+        arc degrees:arc minutes from geometrical horizon for sunrise and
+        sunset, e.g., horizon='+0:00' to use sun center crossing the
+        geometrical horizon to define sunrise and sunset,
+        horizon='-0:34' for when the sun's upper edge crosses the
+        geometrical horizon
+
+    Returns
+    -------
+    pandas.DataFrame
+        index is the same as input `time` argument
+        columns are 'sunrise', 'sunset', and 'transit'
+
+    See also
+    --------
+    pyephem
+    """
+
+    try:
+        import ephem
+    except ImportError:
+        raise ImportError('PyEphem must be installed')
+
+    # times must be localized
+    if not time.tz:
+        raise ValueError('rise_set_ephem: times must be localized')
+
+    obs, sun = _ephem_setup(latitude, longitude, altitude,
+                            pressure, temperature, horizon)
+    # create lists of sunrise and sunset time localized to time.tz
+    if next_or_previous.lower() == 'next':
+        rising = obs.next_rising
+        setting = obs.next_setting
+        transit = obs.next_transit
+    elif next_or_previous.lower() == 'previous':
+        rising = obs.previous_rising
+        setting = obs.previous_setting
+        transit = obs.previous_transit
+    else:
+        raise ValueError("next_or_previous must be either 'next' or" +
+                         " 'previous'")
+
+    sunrise = []
+    sunset = []
+    trans = []
+    for thetime in time:
+        thetime = thetime.to_pydatetime()
+        # pyephem drops timezone when converting to its internal datetime
+        # format, so handle timezone explicitly here
+        obs.date = ephem.Date(thetime - thetime.utcoffset())
+        sunrise.append(_ephem_to_timezone(rising(sun), time.tz))
+        sunset.append(_ephem_to_timezone(setting(sun), time.tz))
+        trans.append(_ephem_to_timezone(transit(sun), time.tz))
+
+    return pd.DataFrame(index=time, data={'sunrise': sunrise,
+                                          'sunset': sunset,
+                                          'transit': trans})
+
+
 def pyephem(time, latitude, longitude, altitude=0, pressure=101325,
-            temperature=12):
+            temperature=12, horizon='+0:00'):
     """
     Calculate the solar position using the PyEphem package.
 
@@ -463,17 +566,26 @@ def pyephem(time, latitude, longitude, altitude=0, pressure=101325,
     time : pandas.DatetimeIndex
         Localized or UTC.
     latitude : float
+        positive is north of 0
     longitude : float
+        positive is east of 0
     altitude : float, default 0
-        distance above sea level.
+        distance above sea level in meters.
     pressure : int or float, optional, default 101325
         air pressure in Pascals.
     temperature : int or float, optional, default 12
         air temperature in degrees C.
+    horizon : string, optional, default '+0:00'
+        arc degrees:arc minutes from geometrical horizon for sunrise and
+        sunset, e.g., horizon='+0:00' to use sun center crossing the
+        geometrical horizon to define sunrise and sunset,
+        horizon='-0:34' for when the sun's upper edge crosses the
+        geometrical horizon
 
     Returns
     -------
-    DataFrame
+    pandas.DataFrame
+        index is the same as input `time` argument
         The DataFrame will have the following columns:
         apparent_elevation, elevation,
         apparent_azimuth, azimuth,
@@ -499,7 +611,7 @@ def pyephem(time, latitude, longitude, altitude=0, pressure=101325,
     sun_coords = pd.DataFrame(index=time)
 
     obs, sun = _ephem_setup(latitude, longitude, altitude,
-                            pressure, temperature)
+                            pressure, temperature, horizon)
 
     # make and fill lists of the sun's altitude and azimuth
     # this is the pressure and temperature corrected apparent alt/az.
@@ -711,7 +823,8 @@ def ephemeris(time, latitude, longitude, pressure=101325, temperature=12):
 
 
 def calc_time(lower_bound, upper_bound, latitude, longitude, attribute, value,
-              altitude=0, pressure=101325, temperature=12, xtol=1.0e-12):
+              altitude=0, pressure=101325, temperature=12, horizon='+0:00',
+              xtol=1.0e-12):
     """
     Calculate the time between lower_bound and upper_bound
     where the attribute is equal to value. Uses PyEphem for
@@ -736,6 +849,12 @@ def calc_time(lower_bound, upper_bound, latitude, longitude, attribute, value,
         atmospheric correction.
     temperature : int or float, optional, default 12
         Air temperature in degrees C.
+    horizon : string, optional, default '+0:00'
+        arc degrees:arc minutes from geometrical horizon for sunrise and
+        sunset, e.g., horizon='+0:00' to use sun center crossing the
+        geometrical horizon to define sunrise and sunset,
+        horizon='-0:34' for when the sun's upper edge crosses the
+        geometrical horizon
     xtol : float, optional, default 1.0e-12
         The allowed error in the result from value
 
@@ -758,7 +877,7 @@ def calc_time(lower_bound, upper_bound, latitude, longitude, attribute, value,
         raise ImportError('The calc_time function requires scipy')
 
     obs, sun = _ephem_setup(latitude, longitude, altitude,
-                            pressure, temperature)
+                            pressure, temperature, horizon)
 
     def compute_attr(thetime, target, attr):
         obs.date = thetime
diff --git a/pvlib/test/test_location.py b/pvlib/test/test_location.py
index b37a6cecc0..ef2b44c180 100644
--- a/pvlib/test/test_location.py
+++ b/pvlib/test/test_location.py
@@ -18,7 +18,7 @@
 import pvlib
 from pvlib.location import Location
 
-from test_solarposition import expected_solpos
+from test_solarposition import expected_solpos, golden_mst
 
 from conftest import requires_scipy
 
@@ -252,8 +252,7 @@ def test_from_tmy_2():
     assert_frame_equal(loc.tmy_data, data)
 
 
-def test_get_solarposition(expected_solpos):
-    from test_solarposition import golden_mst
+def test_get_solarposition(expected_solpos, golden_mst):
     times = pd.date_range(datetime.datetime(2003,10,17,12,30,30),
                           periods=1, freq='D', tz=golden_mst.tz)
     ephem_data = golden_mst.get_solarposition(times, temperature=11)
diff --git a/pvlib/test/test_solarposition.py b/pvlib/test/test_solarposition.py
index b99e5e7ff4..6c407a7a9b 100644
--- a/pvlib/test/test_solarposition.py
+++ b/pvlib/test/test_solarposition.py
@@ -21,14 +21,21 @@
                       end=datetime.datetime(2014,6,26), freq='15Min')
 
 tus = Location(32.2, -111, 'US/Arizona', 700) # no DST issues possible
-# In 2003, DST in US was from April 6 to October 26
-golden_mst = Location(39.742476, -105.1786, 'MST', 1830.14) # no DST issues possible
-golden = Location(39.742476, -105.1786, 'America/Denver', 1830.14) # DST issues possible
-
 times_localized = times.tz_localize(tus.tz)
 
 tol = 5
 
+
+@pytest.fixture()
+def golden():
+    return Location(39.742476, -105.1786, 'America/Denver', 1830.14)
+
+
+@pytest.fixture()
+def golden_mst():
+    return Location(39.742476, -105.1786, 'MST', 1830.14)
+
+
 @pytest.fixture()
 def expected_solpos():
     return pd.DataFrame({'elevation': 39.872046,
@@ -37,6 +44,7 @@ def expected_solpos():
                          'apparent_elevation': 39.888378},
                         index=['2003-10-17T12:30:30Z'])
 
+
 @pytest.fixture()
 def expected_solpos_multi():
     return pd.DataFrame({'elevation': [39.872046, 39.505196],
@@ -45,12 +53,61 @@ def expected_solpos_multi():
                          'apparent_elevation': [39.888378, 39.521740]},
                         index=[['2003-10-17T12:30:30Z', '2003-10-18T12:30:30Z']])
 
+
+@pytest.fixture()
+def expected_rise_set_spa():
+    # for Golden, CO, from NREL SPA website
+    times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2),
+                              datetime.datetime(2015, 8, 2),
+                              ]).tz_localize('MST')
+    sunrise = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 7, 21, 55),
+                                datetime.datetime(2015, 8, 2, 5, 0, 27)
+                                ]).tz_localize('MST').tolist()
+    sunset = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 16, 47, 43),
+                               datetime.datetime(2015, 8, 2, 19, 13, 58)
+                               ]).tz_localize('MST').tolist()
+    transit = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 12, 4, 45),
+                                datetime.datetime(2015, 8, 2, 12, 6, 58)
+                                ]).tz_localize('MST').tolist()
+    return pd.DataFrame({'transit': transit,
+                         'sunrise': sunrise,
+                         'sunset': sunset},
+                        index=times)
+
+
+@pytest.fixture()
+def expected_rise_set_ephem():
+    # for Golden, CO, from USNO websites
+    times = pd.DatetimeIndex([datetime.datetime(2015, 1, 1),
+                              datetime.datetime(2015, 1, 2),
+                              datetime.datetime(2015, 1, 3),
+                              datetime.datetime(2015, 8, 2),
+                              ]).tz_localize('MST')
+    sunrise = pd.DatetimeIndex([datetime.datetime(2015, 1, 1, 7, 22, 0),
+                                datetime.datetime(2015, 1, 2, 7, 22, 0),
+                                datetime.datetime(2015, 1, 3, 7, 22, 0),
+                                datetime.datetime(2015, 8, 2, 5, 0, 0)
+                                ]).tz_localize('MST').tolist()
+    sunset = pd.DatetimeIndex([datetime.datetime(2015, 1, 1, 16, 47, 0),
+                               datetime.datetime(2015, 1, 2, 16, 48, 0),
+                               datetime.datetime(2015, 1, 3, 16, 49, 0),
+                               datetime.datetime(2015, 8, 2, 19, 13, 0)
+                               ]).tz_localize('MST').tolist()
+    transit = pd.DatetimeIndex([datetime.datetime(2015, 1, 1, 12, 4, 0),
+                                datetime.datetime(2015, 1, 2, 12, 5, 0),
+                                datetime.datetime(2015, 1, 3, 12, 5, 0),
+                                datetime.datetime(2015, 8, 2, 12, 7, 0)
+                                ]).tz_localize('MST').tolist()
+    return pd.DataFrame({'sunrise': sunrise,
+                         'sunset': sunset,
+                         'transit': transit},
+                        index=times)
 # the physical tests are run at the same time as the NREL SPA test.
 # pyephem reproduces the NREL result to 2 decimal places.
 # this doesn't mean that one code is better than the other.
 
 @requires_spa_c
-def test_spa_c_physical(expected_solpos):
+def test_spa_c_physical(expected_solpos, golden_mst):
     times = pd.date_range(datetime.datetime(2003,10,17,12,30,30),
                           periods=1, freq='D', tz=golden_mst.tz)
     ephem_data = solarposition.spa_c(times, golden_mst.latitude,
@@ -62,7 +119,7 @@ def test_spa_c_physical(expected_solpos):
 
 
 @requires_spa_c
-def test_spa_c_physical_dst(expected_solpos):
+def test_spa_c_physical_dst(expected_solpos, golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=1, freq='D', tz=golden.tz)
     ephem_data = solarposition.spa_c(times, golden.latitude,
@@ -72,7 +129,8 @@ def test_spa_c_physical_dst(expected_solpos):
     expected_solpos.index = times
     assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])
 
-def test_spa_python_numpy_physical(expected_solpos):
+
+def test_spa_python_numpy_physical(expected_solpos, golden_mst):
     times = pd.date_range(datetime.datetime(2003,10,17,12,30,30),
                           periods=1, freq='D', tz=golden_mst.tz)
     ephem_data = solarposition.spa_python(times, golden_mst.latitude,
@@ -84,7 +142,8 @@ def test_spa_python_numpy_physical(expected_solpos):
     expected_solpos.index = times
     assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])
 
-def test_spa_python_numpy_physical_dst(expected_solpos):
+
+def test_spa_python_numpy_physical_dst(expected_solpos, golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=1, freq='D', tz=golden.tz)
     ephem_data = solarposition.spa_python(times, golden.latitude,
@@ -98,7 +157,7 @@ def test_spa_python_numpy_physical_dst(expected_solpos):
 
 
 @requires_numba
-def test_spa_python_numba_physical(expected_solpos):
+def test_spa_python_numba_physical(expected_solpos, golden_mst):
     times = pd.date_range(datetime.datetime(2003,10,17,12,30,30),
                           periods=1, freq='D', tz=golden_mst.tz)
     ephem_data = solarposition.spa_python(times, golden_mst.latitude,
@@ -112,7 +171,7 @@ def test_spa_python_numba_physical(expected_solpos):
 
 
 @requires_numba
-def test_spa_python_numba_physical_dst(expected_solpos):
+def test_spa_python_numba_physical_dst(expected_solpos, golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=1, freq='D', tz=golden.tz)
     ephem_data = solarposition.spa_python(times, golden.latitude,
@@ -125,7 +184,7 @@ def test_spa_python_numba_physical_dst(expected_solpos):
 
 
 @needs_pandas_0_17
-def test_get_sun_rise_set_transit():
+def test_get_sun_rise_set_transit(expected_rise_set_spa, golden):
     south = Location(-35.0, 0.0, tz='UTC')
     times = pd.DatetimeIndex([datetime.datetime(1996, 7, 5, 0),
                               datetime.datetime(2004, 12, 4, 0)]
@@ -134,53 +193,186 @@ def test_get_sun_rise_set_transit():
                                 datetime.datetime(2004, 12, 4, 4, 38, 57)]
                                ).tz_localize('UTC').tolist()
     sunset = pd.DatetimeIndex([datetime.datetime(1996, 7, 5, 17, 1, 4),
-                               datetime.datetime(2004, 12, 4, 19, 2, 2)]
+                               datetime.datetime(2004, 12, 4, 19, 2, 3)]
                               ).tz_localize('UTC').tolist()
+    frame = pd.DataFrame({'sunrise': sunrise, 'sunset': sunset}, index=times)
+
     result = solarposition.get_sun_rise_set_transit(times, south.latitude,
                                                     south.longitude,
-                                                    delta_t=64.0)
-    frame = pd.DataFrame({'sunrise':sunrise, 'sunset':sunset}, index=times)
+                                                    delta_t=65.0)
     result_rounded = pd.DataFrame(index=result.index)
     # need to iterate because to_datetime does not accept 2D data
     # the rounding fails on pandas < 0.17
     for col, data in result.iteritems():
-        result_rounded[col] = pd.to_datetime(
-            np.floor(data.values.astype(np.int64) / 1e9)*1e9, utc=True)
+        result_rounded[col] = data.dt.round('1s')
 
     del result_rounded['transit']
     assert_frame_equal(frame, result_rounded)
 
+    # test for Golden, CO compare to NREL SPA
+    result = solarposition.get_sun_rise_set_transit(
+        expected_rise_set_spa.index, golden.latitude, golden.longitude,
+        delta_t=65.0)
 
-    # tests from USNO
-    # Golden
-    golden = Location(39.0, -105.0, tz='MST')
-    times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2),
-                              datetime.datetime(2015, 8, 2),]
-                             ).tz_localize('MST')
-    sunrise = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 7, 19, 2),
-                                datetime.datetime(2015, 8, 2, 5, 1, 26)
-                                ]).tz_localize('MST').tolist()
-    sunset = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 16, 49, 10),
-                               datetime.datetime(2015, 8, 2, 19, 11, 31)
-                               ]).tz_localize('MST').tolist()
-    result = solarposition.get_sun_rise_set_transit(times, golden.latitude,
-                                                    golden.longitude,
-                                                    delta_t=64.0)
-    frame = pd.DataFrame({'sunrise':sunrise, 'sunset':sunset}, index=times)
+    # round to nearest minute
     result_rounded = pd.DataFrame(index=result.index)
     # need to iterate because to_datetime does not accept 2D data
-    # the rounding fails on pandas < 0.17
     for col, data in result.iteritems():
-        result_rounded[col] = (pd.to_datetime(
-            np.floor(data.values.astype(np.int64) / 1e9)*1e9, utc=True)
-            .tz_convert('MST'))
+        result_rounded[col] = data.dt.round('s').tz_convert('MST')
+
+    assert_frame_equal(expected_rise_set_spa, result_rounded)
+
+
+@requires_ephem
+def test_rise_set_transit_ephem(expected_rise_set_ephem, golden):
+    # test for Golden, CO compare to USNO, using local midnight
+    result = solarposition.rise_set_transit_ephem(
+        expected_rise_set_ephem.index, golden.latitude, golden.longitude,
+        next_or_previous='next', altitude=golden.altitude, pressure=0,
+        temperature=11, horizon='-0:34')
+    # round to nearest minute
+    result_rounded = pd.DataFrame(index=result.index)
+    for col, data in result.iteritems():
+        result_rounded[col] = data.dt.round('min').tz_convert('MST')
+    assert_frame_equal(expected_rise_set_ephem, result_rounded)
+
+    # test next sunrise/sunset with times
+    times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 3, 0, 0),
+                              datetime.datetime(2015, 1, 2, 10, 15, 0),
+                              datetime.datetime(2015, 1, 2, 15, 3, 0),
+                              datetime.datetime(2015, 1, 2, 21, 6, 7)
+                              ]).tz_localize('MST')
+    expected = pd.DataFrame(index=times,
+                            columns=['sunrise', 'sunset'],
+                            dtype='datetime64[ns]')
+    expected['sunrise'] = pd.Series(index=times, data=
+       [expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunrise'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'sunrise'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'sunrise'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'sunrise']])
+    expected['sunset'] = pd.Series(index=times, data=
+       [expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunset'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunset'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunset'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'sunset']])
+    expected['transit'] = pd.Series(index=times, data=
+       [expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'transit'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'transit'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'transit'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'transit']])
+
+    result = solarposition.rise_set_transit_ephem(times,
+                                                  golden.latitude,
+                                                  golden.longitude,
+                                                  next_or_previous='next',
+                                                  altitude=golden.altitude,
+                                                  pressure=0,
+                                                  temperature=11,
+                                                  horizon='-0:34')
+    # round to nearest minute
+    result_rounded = pd.DataFrame(index=result.index)
+    for col, data in result.iteritems():
+        result_rounded[col] = data.dt.round('min').tz_convert('MST')
+    assert_frame_equal(expected, result_rounded)
+
+    # test previous sunrise/sunset with times
+    times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 3, 0, 0),
+                              datetime.datetime(2015, 1, 2, 10, 15, 0),
+                              datetime.datetime(2015, 1, 3, 3, 0, 0),
+                              datetime.datetime(2015, 1, 3, 13, 6, 7)
+                              ]).tz_localize('MST')
+    expected = pd.DataFrame(index=times,
+                            columns=['sunrise', 'sunset'],
+                            dtype='datetime64[ns]')
+    expected['sunrise'] = pd.Series(index=times, data=
+       [expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 1), 'sunrise'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunrise'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunrise'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'sunrise']])
+    expected['sunset'] = pd.Series(index=times, data=
+       [expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 1), 'sunset'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 1), 'sunset'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunset'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'sunset']])
+    expected['transit'] = pd.Series(index=times, data=
+       [expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 1), 'transit'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 1), 'transit'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 2), 'transit'],
+        expected_rise_set_ephem.loc[datetime.datetime(2015, 1, 3), 'transit']])
+
+    result = solarposition.rise_set_transit_ephem(times,
+                                                  golden.latitude,
+                                                  golden.longitude,
+                                                  next_or_previous='previous',
+                                                  altitude=golden.altitude,
+                                                  pressure=0,
+                                                  temperature=11,
+                                                  horizon='-0:34')
+    # round to nearest minute
+    result_rounded = pd.DataFrame(index=result.index)
+    for col, data in result.iteritems():
+        result_rounded[col] = data.dt.round('min').tz_convert('MST')
+    assert_frame_equal(expected, result_rounded)
+
+    # test with different timezone
+    times = times.tz_convert('UTC')
+    expected = expected.tz_convert('UTC')  # resuse result from previous
+    for col, data in expected.iteritems():
+        expected[col] = data.dt.tz_convert('UTC')
+    result = solarposition.rise_set_transit_ephem(times,
+                                                  golden.latitude,
+                                                  golden.longitude,
+                                                  next_or_previous='previous',
+                                                  altitude=golden.altitude,
+                                                  pressure=0,
+                                                  temperature=11,
+                                                  horizon='-0:34')
+    # round to nearest minute
+    result_rounded = pd.DataFrame(index=result.index)
+    for col, data in result.iteritems():
+        result_rounded[col] = data.dt.round('min').tz_convert(times.tz)
+    assert_frame_equal(expected, result_rounded)
+
+
+@requires_ephem
+def test_rise_set_transit_ephem_error(expected_rise_set_ephem, golden):
+    with pytest.raises(ValueError):
+        solarposition.rise_set_transit_ephem(expected_rise_set_ephem.index,
+                                             golden.latitude,
+                                             golden.longitude,
+                                             next_or_previous='other')
+    tz_naive = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 3, 0, 0)])
+    with pytest.raises(ValueError):
+        solarposition.rise_set_transit_ephem(tz_naive,
+                                             golden.latitude,
+                                             golden.longitude,
+                                             next_or_previous='next')
 
-    del result_rounded['transit']
-    assert_frame_equal(frame, result_rounded)
+
+@requires_ephem
+def test_rise_set_transit_ephem_horizon(golden):
+    times = pd.DatetimeIndex([datetime.datetime(2016, 1, 3, 0, 0, 0)
+                              ]).tz_localize('MST')
+    # center of sun disk
+    center = solarposition.rise_set_transit_ephem(times,
+                                                  latitude=golden.latitude,
+                                                  longitude=golden.longitude)
+    edge = solarposition.rise_set_transit_ephem(times,
+                                                latitude=golden.latitude,
+                                                longitude=golden.longitude,
+                                                horizon='-0:34')
+    result_rounded = (edge['sunrise'] - center['sunrise']).dt.round('min')
+
+    sunrise_delta = datetime.datetime(2016, 1, 3, 7, 17, 11) - \
+        datetime.datetime(2016, 1, 3, 7, 21, 33)
+    expected = pd.Series(index=times,
+                         data=sunrise_delta,
+                         name='sunrise').dt.round('min')
+    assert_series_equal(expected, result_rounded)
 
 
 @requires_ephem
-def test_pyephem_physical(expected_solpos):
+def test_pyephem_physical(expected_solpos, golden_mst):
     times = pd.date_range(datetime.datetime(2003,10,17,12,30,30),
                           periods=1, freq='D', tz=golden_mst.tz)
     ephem_data = solarposition.pyephem(times, golden_mst.latitude,
@@ -190,8 +382,9 @@ def test_pyephem_physical(expected_solpos):
     assert_frame_equal(expected_solpos.round(2),
                        ephem_data[expected_solpos.columns].round(2))
 
+
 @requires_ephem
-def test_pyephem_physical_dst(expected_solpos):
+def test_pyephem_physical_dst(expected_solpos, golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30), periods=1,
                           freq='D', tz=golden.tz)
     ephem_data = solarposition.pyephem(times, golden.latitude,
@@ -201,6 +394,7 @@ def test_pyephem_physical_dst(expected_solpos):
     assert_frame_equal(expected_solpos.round(2),
                        ephem_data[expected_solpos.columns].round(2))
 
+
 @requires_ephem
 def test_calc_time():
     import pytz
@@ -227,6 +421,7 @@ def test_calc_time():
     assert_allclose((az.replace(second=0, microsecond=0) -
                           epoch_dt).total_seconds(), actual_timestamp)
 
+
 @requires_ephem
 def test_earthsun_distance():
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
@@ -235,7 +430,7 @@ def test_earthsun_distance():
     assert_allclose(1, distance, atol=0.1)
 
 
-def test_ephemeris_physical(expected_solpos):
+def test_ephemeris_physical(expected_solpos, golden_mst):
     times = pd.date_range(datetime.datetime(2003,10,17,12,30,30),
                           periods=1, freq='D', tz=golden_mst.tz)
     ephem_data = solarposition.ephemeris(times, golden_mst.latitude,
@@ -248,7 +443,7 @@ def test_ephemeris_physical(expected_solpos):
     assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])
 
 
-def test_ephemeris_physical_dst(expected_solpos):
+def test_ephemeris_physical_dst(expected_solpos, golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=1, freq='D', tz=golden.tz)
     ephem_data = solarposition.ephemeris(times, golden.latitude,
@@ -260,7 +455,7 @@ def test_ephemeris_physical_dst(expected_solpos):
     assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])
 
 
-def test_ephemeris_physical_no_tz(expected_solpos):
+def test_ephemeris_physical_no_tz(expected_solpos, golden_mst):
     times = pd.date_range(datetime.datetime(2003,10,17,19,30,30),
                           periods=1, freq='D')
     ephem_data = solarposition.ephemeris(times, golden_mst.latitude,
@@ -273,7 +468,7 @@ def test_ephemeris_physical_no_tz(expected_solpos):
     assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])
 
 
-def test_get_solarposition_error():
+def test_get_solarposition_error(golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=1, freq='D', tz=golden.tz)
     with pytest.raises(ValueError):
@@ -294,7 +489,7 @@ def test_get_solarposition_error():
                  'equation_of_time', 'zenith'],
         index=expected_solpos().index))
     ])
-def test_get_solarposition_pressure(pressure, expected):
+def test_get_solarposition_pressure(pressure, expected, golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=1, freq='D', tz=golden.tz)
     ephem_data = solarposition.get_solarposition(times, golden.latitude,
@@ -309,16 +504,16 @@ def test_get_solarposition_pressure(pressure, expected):
 
 
 @pytest.mark.parametrize(
-    "altitude, expected", [
-    (golden.altitude, expected_solpos()),
+    "altitude, expected",
+    [(golden().altitude, expected_solpos()),
     (2000, pd.DataFrame(
         np.array([[  39.88788,   50.11212,  194.34024,   39.87205,   14.64151,
                      50.12795]]),
-        columns=['apparent_elevation', 'apparent_zenith', 'azimuth', 'elevation',
-                 'equation_of_time', 'zenith'],
+        columns=['apparent_elevation', 'apparent_zenith', 'azimuth',
+                 'elevation', 'equation_of_time', 'zenith'],
         index=expected_solpos().index))
     ])
-def test_get_solarposition_altitude(altitude, expected):
+def test_get_solarposition_altitude(altitude, expected, golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=1, freq='D', tz=golden.tz)
     ephem_data = solarposition.get_solarposition(times, golden.latitude,
@@ -340,7 +535,7 @@ def test_get_solarposition_altitude(altitude, expected):
     ((None, 'nrel_numba', expected_solpos_multi())),
     (67.0, 'nrel_numba', expected_solpos_multi())
     ])
-def test_get_solarposition_deltat(delta_t, method, expected):
+def test_get_solarposition_deltat(delta_t, method, expected, golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=2, freq='D', tz=golden.tz)
     ephem_data = solarposition.get_solarposition(times, golden.latitude,
@@ -356,7 +551,7 @@ def test_get_solarposition_deltat(delta_t, method, expected):
     assert_frame_equal(this_expected, ephem_data[this_expected.columns])
 
 
-def test_get_solarposition_no_kwargs(expected_solpos):
+def test_get_solarposition_no_kwargs(expected_solpos, golden):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=1, freq='D', tz=golden.tz)
     ephem_data = solarposition.get_solarposition(times, golden.latitude,
@@ -368,7 +563,7 @@ def test_get_solarposition_no_kwargs(expected_solpos):
 
 
 @requires_ephem
-def test_get_solarposition_method_pyephem(expected_solpos):
+def test_get_solarposition_method_pyephem(expected_solpos, golden):
     times = pd.date_range(datetime.datetime(2003, 10, 17, 13, 30, 30),
                           periods=1, freq='D', tz=golden.tz)
     ephem_data = solarposition.get_solarposition(times, golden.latitude,