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

# Standard library
import datetime
import warnings

# Third-party
from astropy.coordinates import (EarthLocation, SkyCoord, AltAz, get_sun,
                                 get_moon, Angle, Latitude, Longitude,
                                 UnitSphericalRepresentation)
from astropy.extern.six import string_types
from astropy.utils import isiterable
from astropy.utils.compat.numpy import broadcast_to
import astropy.units as u
from astropy.time import Time
from astropy.utils import isiterable
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


__all__ = ["Observer", "MAGIC_TIME"]

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


def _generate_24hr_grid(t0, start, end, N, 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 : int
        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-1)],
                                    np.linspace(start, end, N)[1:-1],
                                    [end + 1/(N-1)]])*u.day
    else:
        time_grid = np.linspace(start, end, N)*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')

    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))

    @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=True):
        """
        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: bool
            If True, and the time and target objects cannot be broadcast,
            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. 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:
            # now we broadcast the targets array so that the first index
            # iterates over targets, any other indices over times
            if not target.isscalar:
                if time.isscalar:
                    target = target[:, np.newaxis]
                while target.ndim <= time.ndim:
                    target = target[:, np.newaxis]
        if 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=True):
        """
        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; defaults to `None` (optional)
            Celestial object(s) of interest. If ``target`` is `None`, return the
            `~astropy.coordinates.AltAz` frame without coordinates.

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

        grid: bool
            If True, and the time and target objects cannot be broadcast,
            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. 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:
            time, target = self._preprocess_inputs(time, target, grid)

        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):
        """
        Calculate the parallactic angle.

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

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

        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)

        # Eqn (14.1) of Meeus' Astronomical Algorithms
        LST = time.sidereal_time('mean', longitude=self.location.longitude)
        H = (LST - coordinate.ra).radian
        q = np.arctan(np.sin(H) /
                      (np.tan(self.location.latitude.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 = 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
        array_to_add = np.zeros(before_indices.shape[0])[:, np.newaxis].astype(int)
        array_to_add[1] = 1
        after_indices = before_indices + array_to_add

        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

        """
        slope = (alt_after-alt_before)/((jd_after - jd_before)*u.d)
        crossing_jd = (jd_after*u.d - ((alt_after - horizon)/slope))
        crossing_jd[np.isnan(crossing_jd)] = u.d*MAGIC_TIME.jd
        return np.squeeze(Time(crossing_jd, format='jd'))

    def _altitude_trig(self, LST, target):
        """
        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.

        Returns
        -------
        alt : `~astropy.unit.Quantity`
            Array of altitudes
        """
        LST, target = self._preprocess_inputs(LST, target)
        alt = np.arcsin(np.sin(self.location.latitude.radian) *
                        np.sin(target.dec) +
                        np.cos(self.location.latitude.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=150):
        """
        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 : int
            Number of altitudes to compute when searching for
            rise or set.

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

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

        altaz = self.altaz(times, target, grid=True)
        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=150):
        """
        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 : int
            Number of altitudes to compute when searching for
            rise or set.

        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, for_deriv=True)
        else:
            times = _generate_24hr_grid(time, -1, 0, N, 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=True)
        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)

        # Assemble arguments for function, depending on the function.
        if function == self._calc_riseset:
            args = lambda w: (time, target, w, rise_set, horizon)
        elif function == self._calc_transit:
            args = lambda w: (time, target, w, antitransit)
        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 = function(*args('next'))
            if which == 'next':
                return next_event

        if which == 'previous' or which == 'nearest':
            previous_event = function(*args('previous'))
            if which == 'previous':
                return previous_event

        if which == 'nearest':
            mask = abs(time - previous_event) < abs(time - next_event)
            return Time(np.where(mask, previous_event.utc.jd,
                        next_event.utc.jd), format='jd')


        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):
        """
        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 : coordinate object (i.e. `~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.)

        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")
        >>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_rise_time)) # doctest: +FLOAT_CMP
        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))

    @u.quantity_input(horizon=u.deg)
    def target_set_time(self, time, target, which='nearest', horizon=0*u.degree):
        """
        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 : coordinate object (i.e. `~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.)

        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")
        >>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_set_time)) # doctest: +FLOAT_CMP
        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))

    def target_meridian_transit_time(self, time, target, which='nearest'):
        """
        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 : coordinate object (i.e. `~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

        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")
        >>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_transit_time)) # doctest: +FLOAT_CMP
        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,
                                                rise_set='setting'))

    def target_meridian_antitransit_time(self, time, target, which='nearest'):
        """
        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 : coordinate object (i.e. `~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

        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")
        >>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_antitransit_time)) # doctest: +FLOAT_CMP
        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'))

    @u.quantity_input(horizon=u.deg)
    def sun_rise_time(self, time, which='nearest', horizon=0*u.degree):
        """
        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.)

        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")
        >>> 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)

    @u.quantity_input(horizon=u.deg)
    def sun_set_time(self, time, which='nearest', horizon=0*u.degree):
        """
        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.)

        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")
        >>> 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)

    def noon(self, time, which='nearest'):
        """
        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

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

    def midnight(self, time, which='nearest'):
        """
        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

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

    # Twilight convenience functions

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

        Parameters
        ----------
        time : `~astropy.time.Time` or other (see below)
            Time of observations. 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.

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

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

        Parameters
        ----------
        time : `~astropy.time.Time` or other (see below)
            Time of observations. 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.

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

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

        Parameters
        ----------
        time : `~astropy.time.Time` or other (see below)
            Time of observations. 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.

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

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

        Parameters
        ----------
        time : `~astropy.time.Time` or other (see below)
            Time of observations. 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

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

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

        Parameters
        ----------
        time : `~astropy.time.Time` or other (see below)
            Time of observations. 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.

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

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

        Parameters
        ----------
        time : `~astropy.time.Time` or other (see below)
            Time of observations. 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.

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

    # Moon-related methods.

    def moon_rise_time(self, time, **kwargs):
        """
        Returns the local moonrise time.

        The default moonrise returned is the next one to occur.

        Parameters
        ----------
        time : `~astropy.time.Time` or other (see below)

        Keywords: str, optional
            previous
            next
        """
        raise NotImplementedError()

    def moon_set_time(self, time, **kwargs):
        """
        Returns the local moonset time.

        The default moonset returned is the next one to occur.

        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).

        Keywords: str, optional
            previous
            next
        """
        raise NotImplementedError()

    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=2*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 2*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_moon(time, location=self.location, ephemeris=ephemeris)
        return self.altaz(time, moon, grid=False)

    @u.quantity_input(horizon=u.deg)
    def target_is_up(self, time, target, horizon=0*u.degree, return_altaz=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 : coordinate object (i.e. `~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.

        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)
        True
        >>> apo.target_is_up(time, [aldebaran, vega])
        [True, False]
        """
        if not isinstance(time, Time):
            time = Time(time)

        altaz = self.altaz(time, target)
        observable = altaz.alt > horizon
        if altaz.isscalar:
            observable = bool(observable)
        else:
            # TODO: simply return observable if we move to
            # a fully broadcasted API
            observable = [value for value in observable.flat]

        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:
            # TODO: simply return solar_altitude < horizon if we move to
            # a fully broadcasted API
            return [val for val in (solar_altitude < horizon).flat]

    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.longitude,
                                  model=model)

    def target_hour_angle(self, time, target):
        """
        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 : coordinate object (i.e. `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`) or list
            Target celestial object(s)

        Returns
        -------
        hour_angle : `~astropy.coordinates.Angle`
            The hour angle(s) of the target(s) at ``time``
        """
        time, target = self._preprocess_inputs(time, target)
        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)
        # workaround for NPY <= 1.8, otherwise np.where works even in scalar case
        if current_time.isscalar:
            start_time = current_time if night_mask else sun_set_time
        else:
            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