Skip to main content
`Ctrl`+`K`

# Source code for astroplan.observer

    # Licensed under a 3-clause BSD style license - see LICENSE.rst
    from __future__ import (absolute_import, division, print_function,
                            unicode_literals)
    string_types = (str,)

    # Standard library
    import sys
    import datetime
    import warnings
    # Third-party
    from astropy.coordinates import (EarthLocation, SkyCoord, AltAz, get_sun,
                                     get_body, Angle, Longitude)
    import astropy.units as u
    from astropy.time import Time
    from astropy.utils.exceptions import AstropyDeprecationWarning
    import numpy as np
    import pytz

    # Package
    from .exceptions import TargetNeverUpWarning, TargetAlwaysUpWarning
    from .moon import moon_illumination, moon_phase_angle
    from .target import get_skycoord, SunFlag, MoonFlag

    __all__ = ["Observer"]

    MAGIC_TIME = Time(-999, format='jd')

    # Handle deprecated MAGIC_TIME variable
    def deprecation_wrap_module(mod, deprecated):
        """Return a wrapped object that warns about deprecated accesses"""
        deprecated = set(deprecated)

        class DeprecateWrapper(object):
            def __getattr__(self, attr):
                if attr in deprecated:
                    warnmsg = ("`MAGIC_TIME` will be deprecated in future versions "
                               "of astroplan. Use masked Time objects instead.")
                    warnings.warn(warnmsg, AstropyDeprecationWarning)
                return getattr(mod, attr)

        return DeprecateWrapper()

    sys.modules[__name__] = deprecation_wrap_module(sys.modules[__name__],
                                                    deprecated=['MAGIC_TIME'])

    def _process_nans_in_jds(jds):
        """
        Some functions calculate times for events that won't happen, yielding nans. This wrapper
        manages vectors of (potentially) invalid JDs that must be passed to the astropy.time.Time
        constructor. Returns a masked Time object.
        """
        if np.isscalar(jds):
            jds = np.array([jds])
        not_finite = ~np.isfinite(jds)
        masked_jds = np.ma.masked_array(jds.to(u.day).value if hasattr(jds, 'unit') else jds.copy())
        masked_jds[not_finite] = np.ma.masked
        return masked_jds

    def _generate_24hr_grid(t0, start, end, n_grid_points, for_deriv=False):
        """
        Generate a nearly linearly spaced grid of time durations.

        The midpoints of these grid points will span times from ``t0``+``start``
        to ``t0``+``end``, including the end points, which is useful when taking
        numerical derivatives.

        Parameters
        ----------
        t0 : `~astropy.time.Time`
            Time queried for, grid will be built from or up to this time.

        start : float
            Number of days before/after ``t0`` to start the grid.

        end : float
            Number of days before/after ``t0`` to end the grid.

        n_grid_points : int (optional)
            Number of grid points to generate

        for_deriv : bool
            Generate time series for taking numerical derivative (modify
            bounds)?

        Returns
        -------
        `~astropy.time.Time`
        """

        if for_deriv:
            time_grid = np.concatenate([[start - 1 / (n_grid_points - 1)],
                                        np.linspace(start, end,
                                                    n_grid_points)[1:-1],
                                        [end + 1 / (n_grid_points - 1)]]) * u.day
        else:
            time_grid = np.linspace(start, end, n_grid_points) * u.day

        # broadcast so grid is first index, and remaining shape of t0
        # falls in later indices. e.g. if t0 is shape (10), time_grid
        # will be shape (N, 10). If t0 is shape (5, 2), time_grid is (N, 5, 2)
        while time_grid.ndim <= t0.ndim:
            time_grid = time_grid[:, np.newaxis]
        # we want to avoid 1D grids since we always want to broadcast against targets
        if time_grid.ndim == 1:
            time_grid = time_grid[:, np.newaxis]
        return t0 + time_grid

    class Observer(object):

        """
        A container class for information about an observer's location and
        environment.

        Examples
        --------
        We can create an observer at Subaru Observatory in Hawaii two ways. First,
        locations for some observatories are stored in astroplan, and these can be
        accessed by name, like so:

        >>> from astroplan import Observer
        >>> subaru = Observer.at_site("Subaru", timezone="US/Hawaii")

        To find out which observatories can be accessed by name, check out
        `~astropy.coordinates.EarthLocation.get_site_names`.

        Next, you can initialize an observer by specifying the location with
        `~astropy.coordinates.EarthLocation`:

        >>> from astropy.coordinates import EarthLocation
        >>> import astropy.units as u
        >>> location = EarthLocation.from_geodetic(-155.4761*u.deg, 19.825*u.deg,
        ...                                        4139*u.m)
        >>> subaru = Observer(location=location, name="Subaru", timezone="US/Hawaii")

        You can also create an observer without an
        `~astropy.coordinates.EarthLocation`:

        >>> from astroplan import Observer
        >>> import astropy.units as u
        >>> subaru = Observer(longitude=-155.4761*u.deg, latitude=19.825*u.deg,
        ...                   elevation=0*u.m, name="Subaru", timezone="US/Hawaii")

        """
        @u.quantity_input(elevation=u.m)
        def __init__(self, location=None, timezone='UTC', name=None, latitude=None,
                     longitude=None, elevation=0*u.m, pressure=None,
                     relative_humidity=None, temperature=None, description=None):
            """
            Parameters
            ----------
            location : `~astropy.coordinates.EarthLocation`
                The location (latitude, longitude, elevation) of the observatory.

            timezone : str or `datetime.tzinfo` (optional)
                The local timezone to assume. If a string, it will be passed
                through ``pytz.timezone()`` to produce the timezone object.

            name : str
                A short name for the telescope, observatory or location.

            latitude : float, str, `~astropy.units.Quantity` (optional)
                The latitude of the observing location. Should be valid input for
                initializing a `~astropy.coordinates.Latitude` object.

            longitude : float, str, `~astropy.units.Quantity` (optional)
                The longitude of the observing location. Should be valid input for
                initializing a `~astropy.coordinates.Longitude` object.

            elevation : `~astropy.units.Quantity` (optional), default = 0 meters
                The elevation of the observing location, with respect to sea
                level. Defaults to sea level.

            pressure : `~astropy.units.Quantity` (optional)
                The ambient pressure. Defaults to zero (i.e. no atmosphere).

            relative_humidity : float (optional)
                The ambient relative humidity.

            temperature : `~astropy.units.Quantity` (optional)
                The ambient temperature.

            description : str (optional)
                A short description of the telescope, observatory or observing
                location.
            """

            self.name = name
            self.pressure = pressure
            self.temperature = temperature
            self.relative_humidity = relative_humidity

            # If lat/long given instead of EarthLocation, convert them
            # to EarthLocation
            if location is None and (latitude is not None and
                                     longitude is not None):
                self.location = EarthLocation.from_geodetic(longitude, latitude,
                                                            elevation)

            elif isinstance(location, EarthLocation):
                self.location = location

            else:
                raise TypeError('Observatory location must be specified with '
                                'either (1) an instance of '
                                'astropy.coordinates.EarthLocation or (2) '
                                'latitude and longitude in degrees as '
                                'accepted by astropy.coordinates.Latitude and '
                                'astropy.coordinates.Latitude.')

            # Accept various timezone inputs, default to UTC
            if isinstance(timezone, datetime.tzinfo):
                self.timezone = timezone
            elif isinstance(timezone, string_types):
                self.timezone = pytz.timezone(timezone)
            else:
                raise TypeError('timezone keyword should be a string, or an '
                                'instance of datetime.tzinfo')

        @property
        def longitude(self):
            """The longitude of the observing location, derived from the location."""
            return self.location.lon

        @property
        def latitude(self):
            """The latitude of the observing location, derived from the location."""
            return self.location.lat

        @property
        def elevation(self):
            """The elevation of the observing location with respect to sea level."""
            return self.location.height

        def __repr__(self):
            """
            String representation of the `~astroplan.Observer` object.

            Examples
            --------

            >>> from astroplan import Observer
            >>> keck = Observer.at_site("Keck", timezone="US/Hawaii")
            >>> print(keck)                                    # doctest: +FLOAT_CMP
            <Observer: name='Keck',
                location (lon, lat, el)=(-155.478333333 deg, 19.8283333333 deg, 4160.0 m),
                timezone=<DstTzInfo 'US/Hawaii' LMT-1 day, 13:29:00 STD>>
            """
            class_name = self.__class__.__name__
            attr_names = ['name', 'location', 'timezone', 'pressure', 'temperature',
                          'relative_humidity']
            attr_values = [getattr(self, attr) for attr in attr_names]
            attributes_strings = []
            for name, value in zip(attr_names, attr_values):
                if value is not None:
                    # Format location for easy readability
                    if name == 'location':
                        formatted_loc = ["{} {}".format(i.value, i.unit)
                                         for i in value.to_geodetic()]
                        attributes_strings.append(
                            "{} (lon, lat, el)=({})".format(
                                name, ", ".join(formatted_loc)))
                    else:
                        if name != 'name':
                            value = repr(value)
                        else:
                            value = "'{}'".format(value)
                        attributes_strings.append("{}={}".format(name, value))
            return "<{}: {}>".format(class_name, ",\n    ".join(attributes_strings))

        def _key(self):
            """
            Generate a tuple of the attributes that determine uniqueness of
            `~astroplan.Observer` objects.

            Returns
            -------
            key : tuple

            Examples
            --------

            >>> from astroplan import Observer
            >>> keck = Observer.at_site("Keck", timezone="US/Hawaii")
            >>> keck._key()
            ('Keck', None, None, None, <Longitude -155.47833333 deg>,
                <Latitude 19.82833333 deg>, <Quantity 4160. m>,
                <DstTzInfo 'US/Hawaii' LMT-1 day, 13:29:00 STD>)
            """

            return (self.name,
                    self.pressure,
                    self.temperature,
                    self.relative_humidity,
                    self.longitude,
                    self.latitude,
                    self.elevation,
                    self.timezone,)

        def __hash__(self):
            """
            Hash the `~astroplan.Observer` object.

            Examples
            --------

            >>> from astroplan import Observer
            >>> keck = Observer.at_site("Keck", timezone="US/Hawaii")
            >>> hash(keck)
            -3872382927731250571
            """

            return hash(self._key())

        def __eq__(self, other):
            """
            Equality check for `~astroplan.Observer` objects.

            Examples
            --------

            >>> from astroplan import Observer
            >>> keck = Observer.at_site("Keck", timezone="US/Hawaii")
            >>> keck2 = Observer.at_site("Keck", timezone="US/Hawaii")
            >>> keck == keck2
            True
            """

            if isinstance(other, Observer):
                return self._key() == other._key()
            else:
                return NotImplemented

        def __ne__(self, other):
            """
            Inequality check for `~astroplan.Observer` objects.

            Examples
            --------

            >>> from astroplan import Observer
            >>> keck = Observer.at_site("Keck", timezone="US/Hawaii")
            >>> kpno = Observer.at_site("KPNO", timezone="US/Arizona")
            >>> keck != kpno
            True
            """

            return not self.__eq__(other)

        @classmethod
        def at_site(cls, site_name, **kwargs):
            """
            Initialize an `~astroplan.observer.Observer` object with a site name.

            Extra keyword arguments are passed to the `~astroplan.Observer`
            constructor (see `~astroplan.Observer` for available keyword
            arguments).

            Parameters
            ----------
            site_name : str
                Observatory name, must be resolvable with
                `~astropy.coordinates.EarthLocation.get_site_names`.

            Returns
            -------
            `~astroplan.observer.Observer`
                Observer object.

            Examples
            --------
            Initialize an observer at Kitt Peak National Observatory:

            >>> from astroplan import Observer
            >>> import astropy.units as u
            >>> kpno_generic = Observer.at_site('kpno')
            >>> kpno_today = Observer.at_site('kpno', pressure=1*u.bar, temperature=0*u.deg_C)
            """
            name = kwargs.pop('name', site_name)
            if 'location' in kwargs:
                raise ValueError("Location kwarg should not be used if "
                                 "initializing an Observer with Observer.at_site()")
            return cls(location=EarthLocation.of_site(site_name), name=name, **kwargs)

        def astropy_time_to_datetime(self, astropy_time):
            """
            Convert the `~astropy.time.Time` object ``astropy_time`` to a
            localized `~datetime.datetime` object.

            Timezones localized with `pytz`_.

            .. _pytz: https://pypi.python.org/pypi/pytz/

            Parameters
            ----------
            astropy_time : `~astropy.time.Time`
                Scalar or list-like.

            Returns
            -------
            `~datetime.datetime`
                Localized datetime, where the timezone of the datetime is
                set by the ``timezone`` keyword argument of the
                `~astroplan.Observer` constructor.

            Examples
            --------
            Convert an astropy time to a localized `~datetime.datetime`:

            >>> from astroplan import Observer
            >>> from astropy.time import Time
            >>> subaru = Observer.at_site("Subaru", timezone="US/Hawaii")
            >>> astropy_time = Time('1999-12-31 06:00:00')
            >>> print(subaru.astropy_time_to_datetime(astropy_time))
            1999-12-30 20:00:00-10:00
            """

            if not astropy_time.isscalar:
                return [self.astropy_time_to_datetime(t) for t in astropy_time]

            # Convert astropy.time.Time to a UTC localized datetime (aware)
            utc_datetime = pytz.utc.localize(astropy_time.utc.datetime)

            # Convert UTC to local timezone
            return self.timezone.normalize(utc_datetime)

        def datetime_to_astropy_time(self, date_time):
            """
            Convert the `~datetime.datetime` object ``date_time`` to a
            `~astropy.time.Time` object.

            Timezones localized with `pytz`_. If the ``date_time`` is naive, the
            implied timezone is the ``timezone`` structure of ``self``.

            Parameters
            ----------
            date_time : `~datetime.datetime` or list-like

            Returns
            -------
            `~astropy.time.Time`
                Astropy time object (no timezone information preserved).

            Examples
            --------
            Convert a localized `~datetime.datetime` to a `~astropy.time.Time`
            object. Non-localized datetimes are assumed to be UTC.
            <Time object: scale='utc' format='datetime' value=1999-12-31 06:00:00>

            >>> from astroplan import Observer
            >>> import datetime
            >>> import pytz
            >>> subaru = Observer.at_site("Subaru", timezone="US/Hawaii")
            >>> hi_date_time = datetime.datetime(2005, 6, 21, 20, 0, 0, 0)
            >>> subaru.datetime_to_astropy_time(hi_date_time)
            <Time object: scale='utc' format='datetime' value=2005-06-22 06:00:00>
            >>> utc_date_time = datetime.datetime(2005, 6, 22, 6, 0, 0, 0,
            ...                                   tzinfo=pytz.timezone("UTC"))
            >>> subaru.datetime_to_astropy_time(utc_date_time)
            <Time object: scale='utc' format='datetime' value=2005-06-22 06:00:00>
            """

            if hasattr(date_time, '__iter__'):
                return Time([self.datetime_to_astropy_time(t) for t in date_time])

            # For timezone-naive datetimes, assign local timezone
            if date_time.tzinfo is None:
                date_time = self.timezone.localize(date_time)

            return Time(date_time, location=self.location)

        def _is_broadcastable(self, shp1, shp2):
            """Test if two shape tuples are broadcastable"""
            if shp1 == shp2:
                return True
            for a, b in zip(shp1[::-1], shp2[::-1]):
                if a == 1 or b == 1 or a == b:
                    pass
                else:
                    return False
            return True

        def _preprocess_inputs(self, time, target=None, grid_times_targets=False):
            """
            Preprocess time and target inputs

            This routine takes the inputs for time and target and attempts to
            return a single `~astropy.time.Time` and `~astropy.coordinates.SkyCoord`
            for each argument, which may be non-scalar if necessary.

            time : `~astropy.time.Time` or other (see below)
                The time(s) to use in the calculation. It can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            target : `~astroplan.FixedTarget`, `~astropy.coordinates.SkyCoord`, or list
                The target(s) to use in the calculation.

            grid_times_targets: bool
                If True, the target object will have extra dimensions packed onto
                the end, so that calculations with M targets and N times will
                return an (M, N) shaped result. Otherwise, we rely on broadcasting
                the shapes together using standard numpy rules. Useful for grid
                searches for rise/set times etc.
            """
            # make sure we have a non-scalar time
            if not isinstance(time, Time):
                time = Time(time)

            if target is None:
                return time, None

            # convert any kind of target argument to non-scalar SkyCoord
            target = get_skycoord(target)
            if grid_times_targets:
                if target.isscalar:
                    # ensure we have a (1, 1) shape coord
                    target = SkyCoord(np.tile(target, 1))[:, np.newaxis]
                else:
                    while target.ndim <= time.ndim:
                        target = target[:, np.newaxis]

            elif not self._is_broadcastable(target.shape, time.shape):
                raise ValueError('Time and Target arguments cannot be broadcast '
                                 'against each other with shapes {} and {}'
                                 .format(time.shape, target.shape))
            return time, target

        def altaz(self, time, target=None, obswl=None, grid_times_targets=False):
            """
            Get an `~astropy.coordinates.AltAz` frame or coordinate.

            If ``target`` is None, generates an altitude/azimuth frame. Otherwise,
            calculates the transformation to that frame for the requested ``target``.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                The time at which the observation is taking place. Will be used as
                the ``obstime`` attribute in the resulting frame or coordinate. This
                will be passed in as the first argument to the `~astropy.time.Time`
                initializer, so it can be anything that `~astropy.time.Time` will
                accept (including a `~astropy.time.Time` object)

            target : `~astroplan.FixedTarget`, `~astropy.coordinates.SkyCoord`, or list (optional)
                Celestial object(s) of interest. If ``target`` is `None`, returns
                the `~astropy.coordinates.AltAz` frame without coordinates.

            obswl : `~astropy.units.Quantity` (optional)
                Wavelength of the observation used in the calculation.

            grid_times_targets: bool (optional)
                If True, the target object will have extra dimensions packed
                onto the end, so that calculations with M targets and N times
                will return an (M, N) shaped result. Otherwise, we rely on
                broadcasting the shapes together using standard numpy
                rules. Useful for grid searches for rise/set times etc.

            Returns
            -------
            `~astropy.coordinates.AltAz`
                If ``target`` is `None`, returns `~astropy.coordinates.AltAz`
                frame. If ``target`` is not `None`, returns the ``target``
                transformed to the `~astropy.coordinates.AltAz` frame.

            Examples
            --------
            Create an instance of the `~astropy.coordinates.AltAz` frame for an
            observer at Apache Point Observatory at a particular time:

            >>> from astroplan import Observer
            >>> from astropy.time import Time
            >>> from astropy.coordinates import SkyCoord
            >>> apo = Observer.at_site("APO")
            >>> time = Time('2001-02-03 04:05:06')
            >>> target = SkyCoord(0*u.deg, 0*u.deg)
            >>> altaz_frame = apo.altaz(time)

            Now transform the target's coordinates to the alt/az frame:

            >>> target_altaz = target.transform_to(altaz_frame) # doctest: +SKIP

            Alternatively, construct an alt/az frame and transform the target to
            that frame all in one step:

            >>> target_altaz = apo.altaz(time, target) # doctest: +SKIP
            """
            if target is not None:
                if target is MoonFlag:
                    target = get_body("moon", time, location=self.location)
                elif target is SunFlag:
                    target = get_sun(time)

                time, target = self._preprocess_inputs(time, target, grid_times_targets)

            altaz_frame = AltAz(location=self.location, obstime=time,
                                pressure=self.pressure, obswl=obswl,
                                temperature=self.temperature,
                                relative_humidity=self.relative_humidity)
            if target is None:
                # Return just the frame
                return altaz_frame
            else:
                return target.transform_to(altaz_frame)

        def parallactic_angle(self, time, target, grid_times_targets=False,
                              kind='mean', model=None):
            """
            Calculate the parallactic angle.

            Parameters
            ----------
            time : `~astropy.time.Time`
                Observation time.

            target : `~astroplan.FixedTarget` or `~astropy.coordinates.SkyCoord` or list
                Target celestial object(s).

            grid_times_targets: bool
                If True, the target object will have extra dimensions packed onto
                the end, so that calculations with M targets and N times will
                return an (M, N) shaped result. Otherwise, we rely on broadcasting
                the shapes together using standard numpy rules.

            kind : str
                Argument to the `~astropy.time.Time.sidereal_time` function, which may
                be either ``'mean'`` or ``'apparent'``, i.e., accounting for precession only,
                or also for nutation.

            model : str or None; optional
                Precession (and nutation) model to use in the call to
                `~astropy.time.Time.sidereal_time`, see docs for that method for details.

            Returns
            -------
            `~astropy.coordinates.Angle`
                Parallactic angle.

            Notes
            -----
            The parallactic angle is the angle between the great circle that
            intersects a celestial object and the zenith, and the object's hour
            circle [1]_.

            .. [1] https://en.wikipedia.org/wiki/Parallactic_angle

            """
            time, coordinate = self._preprocess_inputs(time, target, grid_times_targets)

            # Eqn (14.1) of Meeus' Astronomical Algorithms
            LST = time.sidereal_time(kind, longitude=self.location.lon, model=model)
            H = (LST - coordinate.ra).radian
            q = np.arctan2(np.sin(H),
                           (np.tan(self.location.lat.radian) *
                            np.cos(coordinate.dec.radian) -
                            np.sin(coordinate.dec.radian)*np.cos(H)))*u.rad

            return Angle(q)

        # Sun-related methods.
        @u.quantity_input(horizon=u.deg)
        def _horiz_cross(self, t, alt, rise_set, horizon=0*u.degree):
            """
            Find time ``t`` when values in array ``a`` go from
            negative to positive or positive to negative (exclude endpoints)

            ``return_limits`` will return nearest times to zero-crossing.

            Parameters
            ----------
            t : `~astropy.time.Time`
                Grid of N times, any shape. Search grid along first axis, e.g (N, ...)
            alt : `~astropy.units.Quantity`
                Grid of altitudes
                Depending on broadcasting we either have ndim >=3 and
                M targets along first axis, e.g (M, N, ...), or
                ndim = 2 and targets/times in last axis
            rise_set : {"rising",  "setting"}
                Calculate either rising or setting across the horizon
            horizon : float
                Number of degrees above/below actual horizon to use
                for calculating rise/set times (i.e.,
                -6 deg horizon = civil twilight, etc.)

            Returns
            -------
            Returns the lower and upper limits on the time and altitudes
            of the horizon crossing. The altitude limits have shape (M, ...) and the
            time limits have shape (...). These arrays aresuitable for interpolation
            to find the horizon crossing time.
            """
            # handle different cases by enforcing standard shapes on
            # the altitude grid
            finesse_time_indexes = False
            if alt.ndim == 1:
                raise ValueError('Must supply more at least a 2D grid of altitudes')
            elif alt.ndim == 2:
                # TODO: this test for ndim=2 doesn't work. if times is e.g (2,5)
                # then alt will have ndim=3, but shape (100, 2, 5) so grid
                # is in first index...
                ntargets = alt.shape[1]
                ngrid = alt.shape[0]
                unit = alt.unit
                alt = np.broadcast_to(alt, (ntargets, ngrid, ntargets)).T
                alt = alt*unit
                extra_dimension_added = True
                if t.shape[1] == 1:
                    finesse_time_indexes = True
            else:
                extra_dimension_added = False
            output_shape = (alt.shape[0],) + alt.shape[2:]

            if rise_set == 'rising':
                # Find index where altitude goes from below to above horizon
                condition = ((alt[:, :-1, ...] < horizon) *
                             (alt[:, 1:, ...] > horizon))
            elif rise_set == 'setting':
                # Find index where altitude goes from above to below horizon
                condition = ((alt[:, :-1, ...] > horizon) *
                             (alt[:, 1:, ...] < horizon))

            noncrossing_indices = np.sum(condition, axis=1, dtype=np.intp) < 1
            alt_lims1 = u.Quantity(np.zeros(output_shape), unit=u.deg)
            alt_lims2 = u.Quantity(np.zeros(output_shape), unit=u.deg)
            jd_lims1 = np.zeros(output_shape)
            jd_lims2 = np.zeros(output_shape)
            if np.any(noncrossing_indices):
                for target_index in set(np.where(noncrossing_indices)[0]):
                    warnmsg = ('Target with index {} does not cross horizon={} '
                               'within 24 hours'.format(target_index, horizon))
                    if (alt[target_index, ...] > horizon).all():
                        warnings.warn(warnmsg, TargetAlwaysUpWarning)
                    else:
                        warnings.warn(warnmsg, TargetNeverUpWarning)

                alt_lims1[np.nonzero(noncrossing_indices)] = np.nan
                alt_lims2[np.nonzero(noncrossing_indices)] = np.nan
                jd_lims1[np.nonzero(noncrossing_indices)] = np.nan
                jd_lims2[np.nonzero(noncrossing_indices)] = np.nan

            before_indices = np.array(np.nonzero(condition))
            # we want to add an vector like (0, 1, ...) to get after indices
            after_indices = before_indices.copy()
            after_indices[1, :] += 1

            al1 = alt[tuple(before_indices)]
            al2 = alt[tuple(after_indices)]
            # slice the time in the same way, but delete the object index
            before_time_index_tuple = np.delete(before_indices, 0, 0)
            after_time_index_tuple = np.delete(after_indices, 0, 0)
            if finesse_time_indexes:
                before_time_index_tuple[1:] = 0
                after_time_index_tuple[1:] = 0
            tl1 = t[tuple(before_time_index_tuple)]
            tl2 = t[tuple(after_time_index_tuple)]

            alt_lims1[tuple(np.delete(before_indices, 1, 0))] = al1
            alt_lims2[tuple(np.delete(before_indices, 1, 0))] = al2
            jd_lims1[tuple(np.delete(before_indices, 1, 0))] = tl1.utc.jd
            jd_lims2[tuple(np.delete(before_indices, 1, 0))] = tl2.utc.jd

            if extra_dimension_added:
                return (alt_lims1.diagonal(), alt_lims2.diagonal(),
                        jd_lims1.diagonal(), jd_lims2.diagonal())
            else:
                return alt_lims1, alt_lims2, jd_lims1, jd_lims2

        @u.quantity_input(horizon=u.deg)
        def _two_point_interp(self, jd_before, jd_after,
                              alt_before, alt_after, horizon=0*u.deg):
            """
            Do linear interpolation between two ``altitudes`` at
            two ``times`` to determine the time where the altitude
            goes through zero.

            Parameters
            ----------
            jd_before : `float`
                JD(UTC) before crossing event

            jd_after : `float`
                JD(UTC) after crossing event

            alt_before : `~astropy.units.Quantity`
                altitude before crossing event

            alt_after : `~astropy.units.Quantity`
                altitude after crossing event

            horizon : `~astropy.units.Quantity`
                Solve for the time when the altitude is equal to
                reference_alt.

            Returns
            -------
            t : `~astropy.time.Time`
                Time when target crosses the horizon

            """
            # Approximate the horizon-crossing time:
            slope = (alt_after-alt_before) / ((jd_after - jd_before) * u.d)
            crossing_jd = (jd_after * u.d - ((alt_after - horizon) / slope))

            # Determine whether or not there are NaNs in the crossing_jd array which
            # represent computations where no horizon crossing was found:
            times = Time(_process_nans_in_jds(crossing_jd), format='jd')

            return np.squeeze(times)

        def _altitude_trig(self, LST, target, grid_times_targets=False):
            """
            Calculate the altitude of ``target`` at local sidereal times ``LST``.

            This method provides a factor of ~3 speed up over calling `altaz`, and
            inherently does *not* take the atmosphere into account.

            Parameters
            ----------
            LST : `~astropy.time.Time`
                Local sidereal times (array)

            target : {`~astropy.coordinates.SkyCoord`, `FixedTarget`} or similar
                Target celestial object's coordinates.

            grid_times_targets: bool
                If True, the target object will have extra dimensions packed onto
                the end, so that calculations with M targets and N times will
                return an (M, N) shaped result. Otherwise, we rely on broadcasting
                the shapes together using standard numpy rules.

            Returns
            -------
            alt : `~astropy.unit.Quantity`
                Array of altitudes
            """
            LST, target = self._preprocess_inputs(LST, target, grid_times_targets)
            alt = np.arcsin(np.sin(self.location.lat.radian) *
                            np.sin(target.dec) +
                            np.cos(self.location.lat.radian) *
                            np.cos(target.dec) *
                            np.cos(LST.radian - target.ra.radian))
            return alt

        def _calc_riseset(self, time, target, prev_next, rise_set, horizon,
                          n_grid_points=150, grid_times_targets=False):
            """
            Time at next rise/set of ``target``.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument to
                the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            target : `~astropy.coordinates.SkyCoord`
                Position of target or multiple positions of that target
                at multiple times (if target moves, like the Sun)

            prev_next : str - either 'previous' or 'next'
                Test next rise/set or previous rise/set

            rise_set : str - either 'rising' or 'setting'
                Compute prev/next rise or prev/next set

            horizon : `~astropy.units.Quantity`
                Degrees above/below actual horizon to use
                for calculating rise/set times (i.e.,
                -6 deg horizon = civil twilight, etc.)

            n_grid_points : int (optional)
                Number of altitudes to compute when searching for
                rise or set.

            grid_times_targets: bool
                If True, the target object will have extra dimensions packed onto
                the end, so that calculations with M targets and N times will
                return an (M, N) shaped result. Otherwise, we rely on broadcasting
                the shapes together using standard numpy rules.

            Returns
            -------
            ret1 : `~astropy.time.Time`
                Time of rise/set
            """
            if not isinstance(time, Time):
                time = Time(time)

            if prev_next == 'next':
                start = 0
                end = (1 + (target.approx_sidereal_drift.to(u.day).value
                            if hasattr(target, 'approx_sidereal_drift') else 0))
            else:
                start = (-1 - (target.approx_sidereal_drift.to(u.day).value
                               if hasattr(target, 'approx_sidereal_drift') else 0))
                end = 0

            times = _generate_24hr_grid(time, start, end, n_grid_points)

            altaz = self.altaz(times, target,
                               grid_times_targets=grid_times_targets)

            altitudes = altaz.alt

            al1, al2, jd1, jd2 = self._horiz_cross(times, altitudes, rise_set,
                                                   horizon)
            return self._two_point_interp(jd1, jd2, al1, al2,
                                          horizon=horizon)

        def _calc_transit(self, time, target, prev_next, antitransit=False,
                          n_grid_points=150, grid_times_targets=False):
            """
            Time at next transit of the meridian of `target`.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument to
                the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            target : `~astropy.coordinates.SkyCoord`
                Position of target or multiple positions of that target
                at multiple times (if target moves, like the Sun)

            prev_next : str - either 'previous' or 'next'
                Test next rise/set or previous rise/set

            antitransit : bool
                Toggle compute antitransit (below horizon, equivalent to midnight
                for the Sun)

            n_grid_points : int (optional)
                Number of altitudes to compute when searching for
                rise or set.

            grid_times_targets: bool
                If True, the target object will have extra dimensions packed onto
                the end, so that calculations with M targets and N times will
                return an (M, N) shaped result. Otherwise, we rely on broadcasting
                the shapes together using standard numpy rules.

            Returns
            -------
            ret1 : `~astropy.time.Time`
                Time of transit/antitransit
            """
            # TODO FIX BROADCASTING HERE
            if not isinstance(time, Time):
                time = Time(time)

            if prev_next == 'next':
                times = _generate_24hr_grid(time, 0, 1, n_grid_points,
                                            for_deriv=True)
            else:
                times = _generate_24hr_grid(time, -1, 0, n_grid_points,
                                            for_deriv=True)

            # The derivative of the altitude with respect to time is increasing
            # from negative to positive values at the anti-transit of the meridian
            if antitransit:
                rise_set = 'rising'
            else:
                rise_set = 'setting'

            altaz = self.altaz(times, target,
                               grid_times_targets=grid_times_targets)

            altitudes = altaz.alt
            if altitudes.ndim > 2:
                # shape is (M, N, ...) where M is targets and N is grid
                d_altitudes = altitudes.diff(axis=1)
            else:
                # shape is (N, M) where M is targets and N is grid
                d_altitudes = altitudes.diff(axis=0)

            dt = Time((times.jd[1:] + times.jd[:-1])/2, format='jd')

            horizon = 0*u.degree  # Find when derivative passes through zero
            al1, al2, jd1, jd2 = self._horiz_cross(dt, d_altitudes,
                                                   rise_set, horizon)
            return self._two_point_interp(jd1, jd2, al1, al2,
                                          horizon=horizon)

        def _determine_which_event(self, function, args_dict):
            """
            Run through the next/previous/nearest permutations of the solutions
            to `function(time, ...)`, and return the previous/next/nearest one
            specified by the args stored in args_dict.
            """
            time = args_dict.pop('time', None)
            target = args_dict.pop('target', None)
            which = args_dict.pop('which', None)
            horizon = args_dict.pop('horizon', None)
            rise_set = args_dict.pop('rise_set', None)
            antitransit = args_dict.pop('antitransit', None)
            grid_times_targets = args_dict.pop('grid_times_targets', False)
            n_grid_points = args_dict.pop('n_grid_points', 150)

            # Assemble arguments for function, depending on the function.
            if function == self._calc_riseset:
                def event_function(w):
                    return function(time, target, w, rise_set, horizon,
                                    grid_times_targets=grid_times_targets,
                                    n_grid_points=n_grid_points)
            elif function == self._calc_transit:
                def event_function(w):
                    return function(time, target, w, antitransit=antitransit,
                                    grid_times_targets=grid_times_targets,
                                    n_grid_points=n_grid_points)
            else:
                raise ValueError('Function {} not supported in '
                                 '_determine_which_event.'.format(function))

            if not isinstance(time, Time):
                time = Time(time)

            if which == 'next' or which == 'nearest':
                next_event = event_function('next')
                if which == 'next':
                    return next_event

            if which == 'previous' or which == 'nearest':
                previous_event = event_function('previous')
                if which == 'previous':
                    return previous_event

            if which == 'nearest':
                try:
                    mask = abs(time - previous_event) < abs(time - next_event)
                    ma = np.where(mask, previous_event.utc.jd, next_event.utc.jd)
                except TypeError:
                    # encountered if time is scalar & nan
                    ma = next_event.utc.jd if not next_event.isscalar else [next_event.utc.jd]
                return Time(_process_nans_in_jds(ma), format='jd', scale='utc')

            raise ValueError('"which" kwarg must be "next", "previous" or '
                             '"nearest".')

        @u.quantity_input(horizon=u.deg)
        def target_rise_time(self, time, target, which='nearest',
                             horizon=0*u.degree, grid_times_targets=False, n_grid_points=150):
            """
            Calculate rise time.

            Compute time of the next/previous/nearest rise of the ``target``
            object, where "rise" is defined as the time when the ``target``
            transitions from altitudes below the ``horizon`` to above the
            ``horizon``.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
                Target celestial object(s)

            which : {'next', 'previous', 'nearest'}
                Choose which sunrise relative to the present ``time`` would you
                like to calculate

            horizon : `~astropy.units.Quantity` (optional), default = zero degrees
                Degrees above/below actual horizon to use
                for calculating rise/set times (i.e.,
                -6 deg horizon = civil twilight, etc.)

            grid_times_targets: bool
                If True, the target object will have extra dimensions packed
                onto the end, so that calculations with M targets and N times
                will return an (M, N) shaped result. Otherwise, we rely on
                broadcasting the shapes together using standard numpy rules.

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields rise time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Rise time of target

            Examples
            --------
            Calculate the rise time of Rigel at Keck Observatory:

            >>> from astroplan import Observer, FixedTarget
            >>> from astropy.time import Time
            >>> time = Time("2001-02-03 04:05:06")
            >>> target = FixedTarget.from_name("Rigel")
            >>> keck = Observer.at_site("Keck")
            >>> rigel_rise_time = keck.target_rise_time(time, target, which="next") # doctest: +SKIP
            >>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_rise_time)) # doctest: +SKIP
            ISO: 2001-02-04 00:51:23.330, JD: 2451944.53569
            """
            return self._determine_which_event(self._calc_riseset,
                                               dict(time=time, target=target,
                                                    which=which, rise_set='rising',
                                                    horizon=horizon,
                                                    n_grid_points=n_grid_points,
                                                    grid_times_targets=grid_times_targets))

        @u.quantity_input(horizon=u.deg)
        def target_set_time(self, time, target, which='nearest', horizon=0*u.degree,
                            grid_times_targets=False, n_grid_points=150):
            """
            Calculate set time.

            Compute time of the next/previous/nearest set of ``target``, where
            "set" is defined as when the ``target`` transitions from altitudes
            above ``horizon`` to below ``horizon``.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
                Target celestial object(s)

            which : {'next', 'previous', 'nearest'}
                Choose which sunset relative to the present ``time`` would you
                like to calculate

            horizon : `~astropy.units.Quantity` (optional), default = zero degrees
                Degrees above/below actual horizon to use
                for calculating rise/set times (i.e.,
                -6 deg horizon = civil twilight, etc.)

            grid_times_targets: bool
                If True, the target object will have extra dimensions packed
                onto the end, so that calculations with M targets and N times
                will return an (M, N) shaped result. Otherwise, we rely on
                broadcasting the shapes together using standard numpy rules.

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields set time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Set time of target.

            Examples
            --------
            Calculate the set time of Rigel at Keck Observatory:

            >>> from astroplan import Observer, FixedTarget
            >>> from astropy.time import Time
            >>> time = Time("2001-02-03 04:05:06")
            >>> target = FixedTarget.from_name("Rigel")
            >>> keck = Observer.at_site("Keck")
            >>> rigel_set_time = keck.target_set_time(time, target, which="next") # doctest: +SKIP
            >>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_set_time)) # doctest: +SKIP
            ISO: 2001-02-03 12:29:34.768, JD: 2451944.02054
            """
            return self._determine_which_event(self._calc_riseset,
                                               dict(time=time, target=target,
                                                    which=which,
                                                    rise_set='setting',
                                                    horizon=horizon,
                                                    n_grid_points=n_grid_points,
                                                    grid_times_targets=grid_times_targets))

        def target_meridian_transit_time(self, time, target, which='nearest',
                                         grid_times_targets=False, n_grid_points=150):
            """
            Calculate time at the transit of the meridian.

            Compute time of the next/previous/nearest transit of the ``target``
            object.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
                Target celestial object(s)

            which : {'next', 'previous', 'nearest'}
                Choose which sunrise relative to the present ``time`` would you
                like to calculate

            grid_times_targets: bool
                If True, the target object will have extra dimensions packed
                onto the end, so that calculations with M targets and N times
                will return an (M, N) shaped result. Otherwise, we rely on
                broadcasting the shapes together using standard numpy rules.

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields rise time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Transit time of target

            Examples
            --------
            Calculate the meridian transit time of Rigel at Keck Observatory:

            >>> from astroplan import Observer, FixedTarget
            >>> from astropy.time import Time
            >>> time = Time("2001-02-03 04:05:06")
            >>> target = FixedTarget.from_name("Rigel")
            >>> keck = Observer.at_site("Keck")
            >>> rigel_transit_time = keck.target_meridian_transit_time(time, target,
            ...                                                        which="next") # doctest: +SKIP
            >>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_transit_time)) # doctest: +SKIP
            ISO: 2001-02-03 06:42:26.863, JD: 2451943.77948
            """
            return self._determine_which_event(self._calc_transit,
                                               dict(time=time, target=target,
                                                    which=which,
                                                    n_grid_points=n_grid_points,
                                                    rise_set='setting',
                                                    grid_times_targets=grid_times_targets))

        def target_meridian_antitransit_time(self, time, target, which='nearest',
                                             grid_times_targets=False, n_grid_points=150):
            """
            Calculate time at the antitransit of the meridian.

            Compute time of the next/previous/nearest antitransit of the ``target``
            object.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
                Target celestial object(s)

            which : {'next', 'previous', 'nearest'}
                Choose which sunrise relative to the present ``time`` would you
                like to calculate

            grid_times_targets : bool
                If True, the target object will have extra dimensions packed onto the end,
                so that calculations with M targets and N times will return an (M, N)
                shaped result. Otherwise, we rely on broadcasting the shapes together
                using standard numpy rules.

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields rise time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Antitransit time of target

            Examples
            --------
            Calculate the meridian anti-transit time of Rigel at Keck Observatory:

            >>> from astroplan import Observer, FixedTarget
            >>> from astropy.time import Time
            >>> time = Time("2001-02-03 04:05:06")
            >>> target = FixedTarget.from_name("Rigel")
            >>> keck = Observer.at_site("Keck")
            >>> rigel_antitransit_time = keck.target_meridian_antitransit_time(
            ...     time, target, which="next") # doctest: +SKIP
            >>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_antitransit_time)) # doctest: +SKIP
            ISO: 2001-02-03 18:40:29.761, JD: 2451944.27812

            """
            return self._determine_which_event(self._calc_transit,
                                               dict(time=time, target=target,
                                                    which=which, antitransit=True,
                                                    rise_set='setting',
                                                    n_grid_points=n_grid_points,
                                                    grid_times_targets=grid_times_targets))

        @u.quantity_input(horizon=u.deg)
        def sun_rise_time(self, time, which='nearest', horizon=0*u.degree, n_grid_points=150):
            """
            Time of sunrise.

            Compute time of the next/previous/nearest sunrise, where
            sunrise is defined as when the Sun transitions from altitudes
            below ``horizon`` to above ``horizon``.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which sunrise relative to the present ``time`` would you
                like to calculate.

            horizon : `~astropy.units.Quantity` (optional), default = zero degrees
                Degrees above/below actual horizon to use
                for calculating rise/set times (i.e.,
                -6 deg horizon = civil twilight, etc.)

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields rise time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Time of sunrise

            Examples
            --------
            Calculate the time of the previous sunrise at Apache Point Observatory:

            >>> from astroplan import Observer
            >>> from astropy.time import Time
            >>> apo = Observer.at_site("APO")
            >>> time = Time('2001-02-03 04:05:06')
            >>> sun_rise = apo.sun_rise_time(time, which="previous") # doctest: +SKIP
            >>> print("ISO: {0.iso}, JD: {0.jd}".format(sun_rise)) # doctest: +SKIP
            ISO: 2001-02-02 14:02:50.554, JD: 2451943.08531
            """
            return self.target_rise_time(time, get_sun(time), which, horizon,
                                         n_grid_points=n_grid_points)

        @u.quantity_input(horizon=u.deg)
        def sun_set_time(self, time, which='nearest', horizon=0*u.degree, n_grid_points=150):
            """
            Time of sunset.

            Compute time of the next/previous/nearest sunset, where
            sunset is defined as when the Sun transitions from altitudes
            below ``horizon`` to above ``horizon``.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which sunset relative to the present ``time`` would you
                like to calculate

            horizon : `~astropy.units.Quantity` (optional), default = zero degrees
                Degrees above/below actual horizon to use
                for calculating rise/set times (i.e.,
                -6 deg horizon = civil twilight, etc.)

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields set time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Time of sunset

            Examples
            --------
            Calculate the time of the next sunset at Apache Point Observatory:

            >>> from astroplan import Observer
            >>> from astropy.time import Time
            >>> apo = Observer.at_site("APO")
            >>> time = Time('2001-02-03 04:05:06')
            >>> sun_set = apo.sun_set_time(time, which="next") # doctest: +SKIP
            >>> print("ISO: {0.iso}, JD: {0.jd}".format(sun_set)) # doctest: +SKIP
            ISO: 2001-02-04 00:35:42.102, JD: 2451944.52479
            """
            return self.target_set_time(time, get_sun(time), which, horizon,
                                        n_grid_points=n_grid_points)

        def noon(self, time, which='nearest', n_grid_points=150):
            """
            Time at solar noon.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which noon relative to the present ``time`` would you
                like to calculate

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields noon time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Time at solar noon
            """
            return self.target_meridian_transit_time(time, get_sun(time), which,
                                                     n_grid_points=n_grid_points)

        def midnight(self, time, which='nearest', n_grid_points=150):
            """
            Time at solar midnight.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which noon relative to the present ``time`` would you
                like to calculate

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields midnight time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Time at solar midnight
            """
            return self.target_meridian_antitransit_time(time, get_sun(time), which,
                                                         n_grid_points=n_grid_points)

        # Twilight convenience functions

        def twilight_evening_astronomical(self, time, which='nearest', n_grid_points=150):
            """
            Time at evening astronomical (-18 degree) twilight.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which twilight relative to the present ``time`` would you
                like to calculate. Default is nearest.

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields twilight time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Time of twilight
            """
            return self.sun_set_time(time, which, horizon=-18*u.degree,
                                     n_grid_points=n_grid_points)

        def twilight_evening_nautical(self, time, which='nearest', n_grid_points=150):
            """
            Time at evening nautical (-12 degree) twilight.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which twilight relative to the present ``time`` would you
                like to calculate. Default is nearest.

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields twilight time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Time of twilight
            """
            return self.sun_set_time(time, which, horizon=-12*u.degree,
                                     n_grid_points=n_grid_points)

        def twilight_evening_civil(self, time, which='nearest', n_grid_points=150):
            """
            Time at evening civil (-6 degree) twilight.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which twilight relative to the present ``time`` would you
                like to calculate. Default is nearest.

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields twilight time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Time of twilight
            """
            return self.sun_set_time(time, which, horizon=-6*u.degree,
                                     n_grid_points=n_grid_points)

        def twilight_morning_astronomical(self, time, which='nearest', n_grid_points=150):
            """
            Time at morning astronomical (-18 degree) twilight.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which twilight relative to the present ``time`` would you
                like to calculate

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields twilight time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Time of twilight
            """
            return self.sun_rise_time(time, which, horizon=-18*u.degree,
                                      n_grid_points=n_grid_points)

        def twilight_morning_nautical(self, time, which='nearest', n_grid_points=150):
            """
            Time at morning nautical (-12 degree) twilight.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which twilight relative to the present ``time`` would you
                like to calculate. Default is nearest.

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields twilight time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Time of twilight
            """
            return self.sun_rise_time(time, which, horizon=-12*u.degree,
                                      n_grid_points=n_grid_points)

        def twilight_morning_civil(self, time, which='nearest', n_grid_points=150):
            """
            Time at morning civil (-6 degree) twilight.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which twilight relative to the present ``time`` would you
                like to calculate. Default is nearest.

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields twilight time precisions better than one minute.

            Returns
            -------
            `~astropy.time.Time`
                Time of sunset
            """
            return self.sun_rise_time(time, which, horizon=-6*u.degree,
                                      n_grid_points=n_grid_points)

        # Moon-related methods.

        def moon_rise_time(self, time, which='nearest', horizon=0*u.deg, n_grid_points=150):
            """
            Returns the local moon rise time.

            Compute time of the next/previous/nearest moon rise, where
            moon rise is defined as the time when the moon transitions from
            altitudes below ``horizon`` to above ``horizon``.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which moon rise relative to the present ``time`` would you
                like to calculate.

            horizon : `~astropy.units.Quantity` (optional), default = zero degrees
                Degrees above/below actual horizon to use
                for calculating rise/set times (i.e.,
                -6 deg horizon = civil twilight, etc.)

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields rise time precisions better than one minute.
            """
            return self.target_rise_time(time, MoonFlag, which, horizon,
                                         n_grid_points=n_grid_points)

        def moon_set_time(self, time, which='nearest', horizon=0*u.deg, n_grid_points=150):
            """
            Returns the local moon set time.

            Compute time of the next/previous/nearest moon set, where
            moon set is defined as the time when the moon transitions from
            altitudes below ``horizon`` to above ``horizon``.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            which : {'next', 'previous', 'nearest'}
                Choose which moon set relative to the present ``time`` would you
                like to calculate.

            horizon : `~astropy.units.Quantity` (optional), default = zero degrees
                Degrees above/below actual horizon to use
                for calculating set/set times (i.e.,
                -6 deg horizon = civil twilight, etc.)

            n_grid_points : int (optional)
                The number of grid points on which to search for the horizon
                crossings of the target over a 24 hour period, default is 150 which
                yields set time precisions better than one minute.
            """
            return self.target_set_time(time, MoonFlag, which, horizon,
                                        n_grid_points=n_grid_points)

        def moon_illumination(self, time):
            """
            Calculate the illuminated fraction of the moon.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                This will be passed in as the first argument to
                the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object).

            Returns
            -------
            float
                Fraction of lunar surface illuminated

            Examples
            --------
            How much of the lunar surface is illuminated at 2015-08-29 18:35 UTC,
            which we happen to know is the time of a full moon?

            >>> from astroplan import Observer
            >>> from astropy.time import Time
            >>> apo = Observer.at_site("APO")
            >>> time = Time("2015-08-29 18:35")
            >>> apo.moon_illumination(time) # doctest: +SKIP
            array([ 0.99972487])
            """
            if not isinstance(time, Time):
                time = Time(time)

            return moon_illumination(time)

        def moon_phase(self, time=None):
            """
            Calculate lunar orbital phase.

            For example, phase=pi is "new", phase=0 is "full".

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                This will be passed in as the first argument to
                the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object).

            Returns
            -------
            moon_phase_angle : float
                Orbital phase angle of the moon where pi corresponds to new moon,
                zero corresponds to full moon.

            Examples
            --------
            Calculate the phase of the moon at 2015-08-29 18:35 UTC. Near zero
            radians corresponds to a nearly full moon.

            >>> from astroplan import Observer
            >>> from astropy.time import Time
            >>> apo = Observer.at_site('APO')
            >>> time = Time('2015-08-29 18:35')
            >>> apo.moon_phase(time) # doctest: +SKIP
            <Quantity [ 0.03317537] rad>
            """
            if time is not None and not isinstance(time, Time):
                time = Time(time)

            return moon_phase_angle(time)

        def moon_altaz(self, time, ephemeris=None):
            """
            Returns the position of the moon in alt/az.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                This will be passed in as the first argument to
                the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object).

            ephemeris : str, optional
                Ephemeris to use.  If not given, use the one set with
                ``astropy.coordinates.solar_system_ephemeris.set`` (which is
                set to 'builtin' by default).

            Returns
            -------
            altaz : `~astropy.coordinates.SkyCoord`
                Position of the moon transformed to altitude and azimuth

            Examples
            --------
            Calculate the altitude and azimuth of the moon at Apache Point
            Observatory:

            >>> from astroplan import Observer
            >>> from astropy.time import Time
            >>> apo = Observer.at_site("APO")
            >>> time = Time("2015-08-29 18:35")
            >>> altaz_moon = apo.moon_altaz(time) # doctest: +SKIP
            >>> print("alt: {0.alt}, az: {0.az}".format(altaz_moon)) # doctest: +SKIP
            alt: -63.72706397691421 deg, az: 345.3640380598265 deg
            """
            if not isinstance(time, Time):
                time = Time(time)

            moon = get_body("moon", time, location=self.location, ephemeris=ephemeris)
            return self.altaz(time, moon)

        def sun_altaz(self, time):
            """
            Returns the position of the Sun in alt/az.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                This will be passed in as the first argument to
                the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object).

            ephemeris : str, optional
                Ephemeris to use.  If not given, use the one set with
                ``astropy.coordinates.solar_system_ephemeris.set`` (which is
                set to 'builtin' by default).

            Returns
            -------
            altaz : `~astropy.coordinates.SkyCoord`
                Position of the sun transformed to altitude and azimuth
            """
            if not isinstance(time, Time):
                time = Time(time)

            sun = get_sun(time)
            return self.altaz(time, sun)

        @u.quantity_input(horizon=u.deg)
        def target_is_up(self, time, target, horizon=0*u.degree,
                         return_altaz=False, grid_times_targets=False):
            """
            Is ``target`` above ``horizon`` at this ``time``?

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
                Target celestial object(s)

            horizon : `~astropy.units.Quantity` (optional), default = zero degrees
                Degrees above/below actual horizon to use
                for calculating rise/set times (i.e.,
                -6 deg horizon = civil twilight, etc.)

            return_altaz : bool (optional)
                Also return the '~astropy.coordinates.AltAz' coordinate.

            grid_times_targets: bool
                If True, the target object will have extra dimensions packed
                onto the end, so that calculations with M targets and N times
                will return an (M, N) shaped result. Otherwise, we rely on
                broadcasting the shapes together using standard numpy rules.

            Returns
            -------
            observable : boolean or np.ndarray(bool)
                True if ``target`` is above ``horizon`` at ``time``, else False.

            Examples
            --------
            Are Aldebaran and Vega above the horizon at Apache Point Observatory
            at 2015-08-29 18:35 UTC?

            >>> from astroplan import Observer, FixedTarget
            >>> from astropy.time import Time
            >>> apo = Observer.at_site("APO")
            >>> time = Time("2015-08-29 18:35")
            >>> aldebaran = FixedTarget.from_name("Aldebaran")
            >>> vega = FixedTarget.from_name("Vega")
            >>> apo.target_is_up(time, aldebaran) # doctest: +SKIP
            True
            >>> apo.target_is_up(time, [aldebaran, vega]) # doctest: +SKIP
            array([ True, False], dtype=bool)
            """
            if not isinstance(time, Time):
                time = Time(time)

            altaz = self.altaz(time, target, grid_times_targets=grid_times_targets)
            observable = altaz.alt > horizon

            if altaz.isscalar:
                observable = bool(observable)

            if not return_altaz:
                return observable
            else:
                return observable, altaz

        @u.quantity_input(horizon=u.deg)
        def is_night(self, time, horizon=0*u.deg, obswl=None):
            """
            Is the Sun below ``horizon`` at ``time``?

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            horizon : `~astropy.units.Quantity` (optional), default = zero degrees
                Degrees above/below actual horizon to use
                for calculating day/night (i.e.,
                -6 deg horizon = civil twilight, etc.)

            obswl : `~astropy.units.Quantity` (optional)
                Wavelength of the observation used in the calculation

            Returns
            -------
            sun_below_horizon : bool or np.ndarray(bool)
                `True` if sun is below ``horizon`` at ``time``, else `False`.

            Examples
            --------
            Is it "nighttime" (i.e. is the Sun below ``horizon``) at Apache Point
            Observatory at 2015-08-29 18:35 UTC?

            >>> from astroplan import Observer
            >>> from astropy.time import Time
            >>> apo = Observer.at_site("APO")
            >>> time = Time("2015-08-29 18:35")
            >>> apo.is_night(time) # doctest: +SKIP
            False
            """
            if not isinstance(time, Time):
                time = Time(time)

            solar_altitude = self.altaz(time, target=get_sun(time), obswl=obswl).alt

            if solar_altitude.isscalar:
                return bool(solar_altitude < horizon)
            else:
                return solar_altitude < horizon

        def local_sidereal_time(self, time, kind='apparent', model=None):
            """
            Convert ``time`` to local sidereal time for observer.

            This is a thin wrapper around the `~astropy.time.Time.sidereal_time`
            method.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            kind : {'mean', 'apparent'} (optional)
                Passed to the ``kind`` argument of
                `~astropy.time.Time.sidereal_time`

            model : str or `None`; optional
                The precession/nutation model to assume - see
                `~astropy.time.Time.sidereal_time` for more details.

            Returns
            -------
            `~astropy.coordinates.Longitude`
                Local sidereal time.
            """
            if not isinstance(time, Time):
                time = Time(time)

            return time.sidereal_time(kind, longitude=self.location.lon,
                                      model=model)

        def target_hour_angle(self, time, target, grid_times_targets=False):
            """
            Calculate the local hour angle of ``target`` at ``time``.

            Parameters
            ----------
            time : `~astropy.time.Time` or other (see below)
                Time of observation. This will be passed in as the first argument
                to the `~astropy.time.Time` initializer, so it can be anything that
                `~astropy.time.Time` will accept (including a `~astropy.time.Time`
                object)

            target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
                Target celestial object(s)

            grid_times_targets: bool
                If True, the target object will have extra dimensions packed
                onto the end, so that calculations with M targets and N times
                will return an (M, N) shaped result. Otherwise, we rely on
                broadcasting the shapes together using standard numpy rules.

            Returns
            -------
            hour_angle : `~astropy.coordinates.Angle`
                The hour angle(s) of the target(s) at ``time``
            """
            time, target = self._preprocess_inputs(time, target, grid_times_targets)
            return Longitude(self.local_sidereal_time(time) - target.ra)

        @u.quantity_input(horizon=u.degree)
        def tonight(self, time=None, horizon=0 * u.degree, obswl=None):
            """
            Return a time range corresponding to the nearest night

            This will return a range of `~astropy.time.Time` corresponding to the
            beginning and ending of the night. If in the middle of a given night,
            return times from `~astropy.time.Time.now` until the nearest
            `~astroplan.Observer.sun_rise_time`

            Parameters
            ----------
            time : `~astropy.time.Time` (optional), default = `~astropy.time.Time.now`
                The start time for tonight, which is allowed to be arbitrary.
                See description above for behavior
            horizon : `~astropy.units.Quantity` (optional), default = zero degrees
                Degrees above/below actual horizon to use for calculating rise/set
                times (e.g., -6 deg horizon = civil twilight, etc.)
            obswl : `~astropy.units.Quantity` (optional)
                Wavelength of the observation used in the calculation

            Returns
            -------
            times : `~astropy.time.Time`
                A tuple of times corresponding to the start and end of current night
            """
            current_time = Time.now() if time is None else time
            night_mask = self.is_night(current_time, horizon=horizon, obswl=obswl)
            sun_set_time = self.sun_set_time(current_time, which='next',
                                             horizon=horizon)

            start_time = np.where(night_mask, current_time, sun_set_time)
            # np.where gives us a list of start Times - convert to Time object
            if not isinstance(start_time, Time):
                start_time = Time(start_time)
            end_time = self.sun_rise_time(start_time, which='next', horizon=horizon)

            return start_time, end_time