text: """
  The ``clearsky`` module contains several methods
  to calculate clear sky GHI, DNI, and DHI.
  """
  
  import os
  from collections import OrderedDict
  import calendar
  
  import numpy as np
  import pandas as pd
  from scipy.optimize import minimize_scalar
  from scipy.linalg import hankel
  import h5py
  
  from pvlib import atmosphere, tools
  
  
  def ineichen(apparent_zenith, airmass_absolute, linke_turbidity,
               altitude=0, dni_extra=1364., perez_enhancement=False):
      '''
      Determine clear sky GHI, DNI, and DHI from Ineichen/Perez model.
  
      Implements the Ineichen and Perez clear sky model for global
      horizontal irradiance (GHI), direct normal irradiance (DNI), and
      calculates the clear-sky diffuse horizontal (DHI) component as the
      difference between GHI and DNI*cos(zenith) as presented in [1, 2]. A
      report on clear sky models found the Ineichen/Perez model to have
      excellent performance with a minimal input data set [3].
  
      Default values for monthly Linke turbidity provided by SoDa [4, 5].
  
      Parameters
      -----------
      apparent_zenith : numeric
          Refraction corrected solar zenith angle in degrees.
  
      airmass_absolute : numeric
          Pressure corrected airmass.
  
      linke_turbidity : numeric
          Linke Turbidity.
  
      altitude : numeric, default 0
          Altitude above sea level in meters.
  
      dni_extra : numeric, default 1364
          Extraterrestrial irradiance. The units of ``dni_extra``
          determine the units of the output.
  
      perez_enhancement : bool, default False
          Controls if the Perez enhancement factor should be applied.
          Setting to True may produce spurious results for times when
          the Sun is near the horizon and the airmass is high.
          See https://github.com/pvlib/pvlib-python/issues/435
  
      Returns
      -------
      clearsky : DataFrame (if Series input) or OrderedDict of arrays
          DataFrame/OrderedDict contains the columns/keys
          ``'dhi', 'dni', 'ghi'``.
  
      See also
      --------
      lookup_linke_turbidity
      pvlib.location.Location.get_clearsky
  
      References
      ----------
      .. [1] P. Ineichen and R. Perez, "A New airmass independent formulation for
         the Linke turbidity coefficient", Solar Energy, vol 73, pp. 151-157,
         2002.
  
      .. [2] R. Perez et. al., "A New Operational Model for Satellite-Derived
         Irradiances: Description and Validation", Solar Energy, vol 73, pp.
         307-317, 2002.
  
      .. [3] M. Reno, C. Hansen, and J. Stein, "Global Horizontal Irradiance
         Clear Sky Models: Implementation and Analysis", Sandia National
         Laboratories, SAND2012-2389, 2012.
  
      .. [4] http://www.soda-is.com/eng/services/climat_free_eng.php#c5 (obtained
         July 17, 2012).
  
      .. [5] J. Remund, et. al., "Worldwide Linke Turbidity Information", Proc.
         ISES Solar World Congress, June 2003. Goteborg, Sweden.
      '''
  
      # ghi is calculated using either the equations in [1] by setting
      # perez_enhancement=False (default behavior) or using the model
      # in [2] by setting perez_enhancement=True.
  
      # The NaN handling is a little subtle. The AM input is likely to
      # have NaNs that we'll want to map to 0s in the output. However, we
      # want NaNs in other inputs to propagate through to the output. This
      # is accomplished by judicious use and placement of np.maximum,
      # np.minimum, and np.fmax
  
      # use max so that nighttime values will result in 0s instead of
      # negatives. propagates nans.
      cos_zenith = np.maximum(tools.cosd(apparent_zenith), 0)
  
      tl = linke_turbidity
  
      fh1 = np.exp(-altitude/8000.)
      fh2 = np.exp(-altitude/1250.)
      cg1 = 5.09e-05 * altitude + 0.868
      cg2 = 3.92e-05 * altitude + 0.0387
  
      ghi = np.exp(-cg2*airmass_absolute*(fh1 + fh2*(tl - 1)))
  
      # https://github.com/pvlib/pvlib-python/issues/435
      if perez_enhancement:
          ghi *= np.exp(0.01*airmass_absolute**1.8)
  
      # use fmax to map airmass nans to 0s. multiply and divide by tl to
      # reinsert tl nans
      ghi = cg1 * dni_extra * cos_zenith * tl / tl * np.fmax(ghi, 0)
  
      # From [1] (Following [2] leads to 0.664 + 0.16268 / fh1)
      # See https://github.com/pvlib/pvlib-python/pull/808
      b = 0.664 + 0.163/fh1
      # BncI = "normal beam clear sky radiation"
      bnci = b * np.exp(-0.09 * airmass_absolute * (tl - 1))
      bnci = dni_extra * np.fmax(bnci, 0)
  
      # "empirical correction" SE 73, 157 & SE 73, 312.
      bnci_2 = ((1 - (0.1 - 0.2*np.exp(-tl))/(0.1 + 0.882/fh1)) /
                cos_zenith)
      bnci_2 = ghi * np.fmin(np.fmax(bnci_2, 0), 1e20)
  
      dni = np.minimum(bnci, bnci_2)
  
      dhi = ghi - dni*cos_zenith
  
      irrads = OrderedDict()
      irrads['ghi'] = ghi
      irrads['dni'] = dni
      irrads['dhi'] = dhi
  
      if isinstance(dni, pd.Series):
          irrads = pd.DataFrame.from_dict(irrads)
  
      return irrads
  
  
  def lookup_linke_turbidity(time, latitude, longitude, filepath=None,
                             interp_turbidity=True):
      """
      Look up the Linke Turibidity from the ``LinkeTurbidities.h5``
      data file supplied with pvlib.
  
      Parameters
      ----------
      time : pandas.DatetimeIndex
  
      latitude : float or int
  
      longitude : float or int
  
      filepath : None or string, default None
          The path to the ``.h5`` file.
  
      interp_turbidity : bool, default True
          If ``True``, interpolates the monthly Linke turbidity values
          found in ``LinkeTurbidities.h5`` to daily values.
  
      Returns
      -------
      turbidity : Series
      """
  
      # The .h5 file 'LinkeTurbidities.h5' contains a single 2160 x 4320 x 12
      # matrix of type uint8 called 'LinkeTurbidity'. The rows represent global
      # latitudes from 90 to -90 degrees; the columns represent global longitudes
      # from -180 to 180; and the depth (third dimension) represents months of
      # the year from January (1) to December (12). To determine the Linke
      # turbidity for a position on the Earth's surface for a given month do the
      # following: LT = LinkeTurbidity(LatitudeIndex, LongitudeIndex, month).
      # Note that the numbers within the matrix are 20 * Linke Turbidity,
      # so divide the number from the file by 20 to get the
      # turbidity.
  
      # The nodes of the grid are 5' (1/12=0.0833[arcdeg]) apart.
      # From Section 8 of Aerosol optical depth and Linke turbidity climatology
      # http://www.meteonorm.com/images/uploads/downloads/ieashc36_report_TL_AOD_climatologies.pdf
      # 1st row: 89.9583 S, 2nd row: 89.875 S
      # 1st column: 179.9583 W, 2nd column: 179.875 W
  
      if filepath is None:
          pvlib_path = os.path.dirname(os.path.abspath(__file__))
          filepath = os.path.join(pvlib_path, 'data', 'LinkeTurbidities.h5')
  
      latitude_index = _degrees_to_index(latitude, coordinate='latitude')
      longitude_index = _degrees_to_index(longitude, coordinate='longitude')
  
      with h5py.File(filepath, 'r') as lt_h5_file:
          lts = lt_h5_file['LinkeTurbidity'][latitude_index, longitude_index]
  
      if interp_turbidity:
          linke_turbidity = _interpolate_turbidity(lts, time)
      else:
          months = time.month - 1
          linke_turbidity = pd.Series(lts[months], index=time)
  
      linke_turbidity /= 20.
  
      return linke_turbidity
  
  
  def _is_leap_year(year):
      """Determine if a year is leap year.
  
      Parameters
      ----------
      year : numeric
  
      Returns
      -------
      isleap : array of bools
      """
      isleap = ((np.mod(year, 4) == 0) &
                ((np.mod(year, 100) != 0) | (np.mod(year, 400) == 0)))
      return isleap
  
  
  def _interpolate_turbidity(lts, time):
      """
      Interpolated monthly Linke turbidity onto daily values.
  
      Parameters
      ----------
      lts : np.array
          Monthly Linke turbidity values.
      time : pd.DatetimeIndex
          Times to be interpolated onto.
  
      Returns
      -------
      linke_turbidity : pd.Series
          The interpolated turbidity.
      """
      # Data covers 1 year. Assume that data corresponds to the value at the
      # middle of each month. This means that we need to add previous Dec and
      # next Jan to the array so that the interpolation will work for
      # Jan 1 - Jan 15 and Dec 16 - Dec 31.
      lts_concat = np.concatenate([[lts[-1]], lts, [lts[0]]])
  
      # handle leap years
      try:
          isleap = time.is_leap_year
      except AttributeError:
          year = time.year
          isleap = _is_leap_year(year)
  
      dayofyear = time.dayofyear
      days_leap = _calendar_month_middles(2016)
      days_no_leap = _calendar_month_middles(2015)
  
      # Then we map the month value to the day of year value.
      # Do it for both leap and non-leap years.
      lt_leap = np.interp(dayofyear, days_leap, lts_concat)
      lt_no_leap = np.interp(dayofyear, days_no_leap, lts_concat)
      linke_turbidity = np.where(isleap, lt_leap, lt_no_leap)
  
      linke_turbidity = pd.Series(linke_turbidity, index=time)
  
      return linke_turbidity
  
  
  def _calendar_month_middles(year):
      """List of middle day of each month, used by Linke turbidity lookup"""
      # remove mdays[0] since January starts at mdays[1]
      # make local copy of mdays since we need to change
      # February for leap years
      mdays = np.array(calendar.mdays[1:])
      ydays = 365
      # handle leap years
      if calendar.isleap(year):
          mdays[1] = mdays[1] + 1
          ydays = 366
      middles = np.concatenate(
          [[-calendar.mdays[-1] / 2.0],  # Dec last year
           np.cumsum(mdays) - np.array(mdays) / 2.,  # this year
           [ydays + calendar.mdays[1] / 2.0]])  # Jan next year
      return middles
  
  
  def _degrees_to_index(degrees, coordinate):
      """Transform input degrees to an output index integer. The Linke
      turbidity lookup tables have three dimensions, latitude, longitude, and
      month. Specify a degree value and either 'latitude' or 'longitude' to get
      the appropriate index number for the first two of these index numbers.
  
      Parameters
      ----------
      degrees : float or int
          Degrees of either latitude or longitude.
      coordinate : string
          Specify whether degrees arg is latitude or longitude. Must be set to
          either 'latitude' or 'longitude' or an error will be raised.
  
      Returns
      -------
      index : np.int16
          The latitude or longitude index number to use when looking up values
          in the Linke turbidity lookup table.
      """
      # Assign inputmin, inputmax, and outputmax based on degree type.
      if coordinate == 'latitude':
          inputmin = 90
          inputmax = -90
          outputmax = 2160
      elif coordinate == 'longitude':
          inputmin = -180
          inputmax = 180
          outputmax = 4320
      else:
          raise IndexError("coordinate must be 'latitude' or 'longitude'.")
  
      inputrange = inputmax - inputmin
      scale = outputmax/inputrange  # number of indices per degree
      center = inputmin + 1 / scale / 2  # shift to center of index
      outputmax -= 1  # shift index to zero indexing
      index = (degrees - center) * scale
      err = IndexError('Input, %g, is out of range (%g, %g).' %
                       (degrees, inputmin, inputmax))
  
      # If the index is still out of bounds after rounding, raise an error.
      # 0.500001 is used in comparisons instead of 0.5 to allow for a small
      # margin of error which can occur when dealing with floating point numbers.
      if index > outputmax:
          if index - outputmax <= 0.500001:
              index = outputmax
          else:
              raise err
      elif index < 0:
          if -index <= 0.500001:
              index = 0
          else:
              raise err
      # If the index wasn't set to outputmax or 0, round it and cast it as an
      # integer so it can be used in integer-based indexing.
      else:
          index = int(np.around(index))
  
      return index
  
  
  def haurwitz(apparent_zenith):
      '''
      Determine clear sky GHI using the Haurwitz model.
  
      Implements the Haurwitz clear sky model for global horizontal
      irradiance (GHI) as presented in [1, 2]. A report on clear
      sky models found the Haurwitz model to have the best performance
      in terms of average monthly error among models which require only
      zenith angle [3].
  
      Parameters
      ----------
      apparent_zenith : Series
          The apparent (refraction corrected) sun zenith angle
          in degrees.
  
      Returns
      -------
      ghi : DataFrame
          The modeled global horizonal irradiance in W/m^2 provided
          by the Haurwitz clear-sky model.
  
      References
      ----------
  
      .. [1] B. Haurwitz, "Insolation in Relation to Cloudiness and Cloud
         Density," Journal of Meteorology, vol. 2, pp. 154-166, 1945.
  
      .. [2] B. Haurwitz, "Insolation in Relation to Cloud Type," Journal of
         Meteorology, vol. 3, pp. 123-124, 1946.
  
      .. [3] M. Reno, C. Hansen, and J. Stein, "Global Horizontal Irradiance
         Clear Sky Models: Implementation and Analysis", Sandia National
         Laboratories, SAND2012-2389, 2012.
      '''
  
      cos_zenith = tools.cosd(apparent_zenith.values)
      clearsky_ghi = np.zeros_like(apparent_zenith.values)
      cos_zen_gte_0 = cos_zenith > 0
      clearsky_ghi[cos_zen_gte_0] = (1098.0 * cos_zenith[cos_zen_gte_0] *
                                     np.exp(-0.059/cos_zenith[cos_zen_gte_0]))
  
      df_out = pd.DataFrame(index=apparent_zenith.index,
                            data=clearsky_ghi,
                            columns=['ghi'])
  
      return df_out
  
  
  def simplified_solis(apparent_elevation, aod700=0.1, precipitable_water=1.,
                       pressure=101325., dni_extra=1364.):
      """
      Calculate the clear sky GHI, DNI, and DHI according to the
      simplified Solis model.
  
      Reference [1]_ describes the accuracy of the model as being 15, 20,
      and 18 W/m^2 for the beam, global, and diffuse components. Reference
      [2]_ provides comparisons with other clear sky models.
  
      Parameters
      ----------
      apparent_elevation : numeric
          The apparent elevation of the sun above the horizon (deg).
  
      aod700 : numeric, default 0.1
          The aerosol optical depth at 700 nm (unitless).
          Algorithm derived for values between 0 and 0.45.
  
      precipitable_water : numeric, default 1.0
          The precipitable water of the atmosphere (cm).
          Algorithm derived for values between 0.2 and 10 cm.
          Values less than 0.2 will be assumed to be equal to 0.2.
  
      pressure : numeric, default 101325.0
          The atmospheric pressure (Pascals).
          Algorithm derived for altitudes between sea level and 7000 m,
          or 101325 and 41000 Pascals.
  
      dni_extra : numeric, default 1364.0
          Extraterrestrial irradiance. The units of ``dni_extra``
          determine the units of the output.
  
      Returns
      -------
      clearsky : DataFrame (if Series input) or OrderedDict of arrays
          DataFrame/OrderedDict contains the columns/keys
          ``'dhi', 'dni', 'ghi'``.
  
      References
      ----------
      .. [1] P. Ineichen, "A broadband simplified version of the
         Solis clear sky model," Solar Energy, 82, 758-762 (2008).
  
      .. [2] P. Ineichen, "Validation of models that estimate the clear
         sky global and beam solar irradiance," Solar Energy, 132,
         332-344 (2016).
      """
  
      p = pressure
  
      w = precipitable_water
  
      # algorithm fails for pw < 0.2
      w = np.maximum(w, 0.2)
  
      # this algorithm is reasonably fast already, but it could be made
      # faster by precalculating the powers of aod700, the log(p/p0), and
      # the log(w) instead of repeating the calculations as needed in each
      # function
  
      i0p = _calc_i0p(dni_extra, w, aod700, p)
  
      taub = _calc_taub(w, aod700, p)
      b = _calc_b(w, aod700)
  
      taug = _calc_taug(w, aod700, p)
      g = _calc_g(w, aod700)
  
      taud = _calc_taud(w, aod700, p)
      d = _calc_d(aod700, p)
  
      # this prevents the creation of nans at night instead of 0s
      # it's also friendly to scalar and series inputs
      sin_elev = np.maximum(1.e-30, np.sin(np.radians(apparent_elevation)))
  
      dni = i0p * np.exp(-taub/sin_elev**b)
      ghi = i0p * np.exp(-taug/sin_elev**g) * sin_elev
      dhi = i0p * np.exp(-taud/sin_elev**d)
  
      irrads = OrderedDict()
      irrads['ghi'] = ghi
      irrads['dni'] = dni
      irrads['dhi'] = dhi
  
      if isinstance(dni, pd.Series):
          irrads = pd.DataFrame.from_dict(irrads)
  
      return irrads
  
  
  def _calc_i0p(i0, w, aod700, p):
      """Calculate the "enhanced extraterrestrial irradiance"."""
      p0 = 101325.
      io0 = 1.08 * w**0.0051
      i01 = 0.97 * w**0.032
      i02 = 0.12 * w**0.56
      i0p = i0 * (i02*aod700**2 + i01*aod700 + io0 + 0.071*np.log(p/p0))
  
      return i0p
  
  
  def _calc_taub(w, aod700, p):
      """Calculate the taub coefficient"""
      p0 = 101325.
      tb1 = 1.82 + 0.056*np.log(w) + 0.0071*np.log(w)**2
      tb0 = 0.33 + 0.045*np.log(w) + 0.0096*np.log(w)**2
      tbp = 0.0089*w + 0.13
  
      taub = tb1*aod700 + tb0 + tbp*np.log(p/p0)
  
      return taub
  
  
  def _calc_b(w, aod700):
      """Calculate the b coefficient."""
  
      b1 = 0.00925*aod700**2 + 0.0148*aod700 - 0.0172
      b0 = -0.7565*aod700**2 + 0.5057*aod700 + 0.4557
  
      b = b1 * np.log(w) + b0
  
      return b
  
  
  def _calc_taug(w, aod700, p):
      """Calculate the taug coefficient"""
      p0 = 101325.
      tg1 = 1.24 + 0.047*np.log(w) + 0.0061*np.log(w)**2
      tg0 = 0.27 + 0.043*np.log(w) + 0.0090*np.log(w)**2
      tgp = 0.0079*w + 0.1
      taug = tg1*aod700 + tg0 + tgp*np.log(p/p0)
  
      return taug
  
  
  def _calc_g(w, aod700):
      """Calculate the g coefficient."""
  
      g = -0.0147*np.log(w) - 0.3079*aod700**2 + 0.2846*aod700 + 0.3798
  
      return g
  
  
  def _calc_taud(w, aod700, p):
      """Calculate the taud coefficient."""
  
      # isscalar tests needed to ensure that the arrays will have the
      # right shape in the tds calculation.
      # there's probably a better way to do this.
  
      if np.isscalar(w) and np.isscalar(aod700):
          w = np.array([w])
          aod700 = np.array([aod700])
      elif np.isscalar(w):
          w = np.full_like(aod700, w)
      elif np.isscalar(aod700):
          aod700 = np.full_like(w, aod700)
  
      # set up nan-tolerant masks
      aod700_lt_0p05 = np.full_like(aod700, False, dtype='bool')
      np.less(aod700, 0.05, where=~np.isnan(aod700), out=aod700_lt_0p05)
      aod700_mask = np.array([aod700_lt_0p05, ~aod700_lt_0p05], dtype=np.int)
  
      # create tuples of coefficients for
      # aod700 < 0.05, aod700 >= 0.05
      td4 = 86*w - 13800, -0.21*w + 11.6
      td3 = -3.11*w + 79.4, 0.27*w - 20.7
      td2 = -0.23*w + 74.8, -0.134*w + 15.5
      td1 = 0.092*w - 8.86, 0.0554*w - 5.71
      td0 = 0.0042*w + 3.12, 0.0057*w + 2.94
      tdp = -0.83*(1+aod700)**(-17.2), -0.71*(1+aod700)**(-15.0)
  
      tds = (np.array([td0, td1, td2, td3, td4, tdp]) * aod700_mask).sum(axis=1)
  
      p0 = 101325.
      taud = (tds[4]*aod700**4 + tds[3]*aod700**3 + tds[2]*aod700**2 +
              tds[1]*aod700 + tds[0] + tds[5]*np.log(p/p0))
  
      # be polite about matching the output type to the input type(s)
      if len(taud) == 1:
          taud = taud[0]
  
      return taud
  
  
  def _calc_d(aod700, p):
      """Calculate the d coefficient."""
  
      p0 = 101325.
      dp = 1/(18 + 152*aod700)
      d = -0.337*aod700**2 + 0.63*aod700 + 0.116 + dp*np.log(p/p0)
  
      return d
  
  
  def _calc_stats(data, samples_per_window, sample_interval, H):
      """ Calculates statistics for each window, used by Reno-style clear
      sky detection functions. Does not return the line length statistic
      which is provided by _calc_windowed_stat and _line_length.
  
      Calculations are done on a sliding window defined by the Hankel matrix H.
      Columns in H define the indices for each window. Each window contains
      samples_per_window index values. The first window starts with index 0;
      the last window ends at the last index position in data.
  
      In the calculation of data_slope_nstd, a choice is made here where [1]_ is
      ambiguous. data_slope_nstd is the standard deviation of slopes divided by
      the mean GHI for each interval; see [1]_ Eq. 11. For intervals containing
      e.g. 10 values, there are 9 slope values in the standard deviation, and the
      mean is calculated using all 10 values. Eq. 11 in [1]_ is ambiguous if
      the mean should be calculated using 9 points (left ends of each slope)
      or all 10 points.
  
      Parameters
      ----------
      data : Series
      samples_per_window : int
          Number of data points in each window
      sample_interval : float
          Time in minutes in each sample interval
      H : 2D ndarray
          Hankel matrix defining the indices for each window.
  
      Returns
      -------
      data_mean : Series
          mean of data in each window
      data_max : Series
          maximum of data in each window
      data_slope_nstd : Series
          standard deviation of difference between data points in each window
      data_slope : Series
          difference between successive data points
  
      References
      ----------
      .. [1] Reno, M.J. and C.W. Hansen, "Identification of periods of clear
         sky irradiance in time series of GHI measurements" Renewable Energy,
         v90, p. 520-531, 2016.
      """
  
      data_mean = data.values[H].mean(axis=0)
      data_mean = _to_centered_series(data_mean, data.index, samples_per_window)
      data_max = data.values[H].max(axis=0)
      data_max = _to_centered_series(data_max, data.index, samples_per_window)
      # shift to get forward difference, .diff() is backward difference instead
      data_diff = data.diff().shift(-1)
      data_slope = data_diff / sample_interval
      data_slope_nstd = _slope_nstd_windowed(data_slope.values[:-1], data, H,
                                             samples_per_window, sample_interval)
      data_slope_nstd = data_slope_nstd
  
      return data_mean, data_max, data_slope_nstd, data_slope
  
  
  def _slope_nstd_windowed(slopes, data, H, samples_per_window, sample_interval):
      with np.errstate(divide='ignore', invalid='ignore'):
          nstd = slopes[H[:-1, ]].std(ddof=1, axis=0) \
              / data.values[H].mean(axis=0)
      return _to_centered_series(nstd, data.index, samples_per_window)
  
  
  def _max_diff_windowed(data, H, samples_per_window):
      raw = np.diff(data)
      raw = np.abs(raw[H[:-1, ]]).max(axis=0)
      return _to_centered_series(raw, data.index, samples_per_window)
  
  
  def _line_length_windowed(data, H, samples_per_window,
                            sample_interval):
      raw = np.sqrt(np.diff(data)**2. + sample_interval**2.)
      raw = np.sum(raw[H[:-1, ]], axis=0)
      return _to_centered_series(raw, data.index, samples_per_window)
  
  
  def _to_centered_series(vals, idx, samples_per_window):
      vals = np.pad(vals, ((0, len(idx) - len(vals)),), mode='constant',
                    constant_values=np.nan)
      shift = samples_per_window // 2  # align = 'center' only
      return pd.Series(index=idx, data=vals).shift(shift)
  
  
  def _get_sample_intervals(times, win_length):
      """ Calculates time interval and samples per window for Reno-style clear
      sky detection functions
      """
      deltas = np.diff(times.values) / np.timedelta64(1, '60s')
  
      # determine if we can proceed
      if times.inferred_freq and len(np.unique(deltas)) == 1:
          sample_interval = times[1] - times[0]
          sample_interval = sample_interval.seconds / 60  # in minutes
          samples_per_window = int(win_length / sample_interval)
          return sample_interval, samples_per_window
      else:
          raise NotImplementedError('algorithm does not yet support unequal '
                                    'times. consider resampling your data.')
  
  
  def _clear_sample_index(clear_windows, samples_per_window, align, H):
      """
      Returns indices of clear samples in clear windows
      """
      # H contains indices for each window, e.g. indices for the first window
      # are in first column of H.
      # clear_windows contains one boolean for each window and is aligned
      # by 'align', default to center
      # shift clear_windows.index to be aligned left (e.g. first value in the
      # left-most position) to line up with the first column of H.
  
      # commented if/else block for future align='left', 'right' capability
      # if align == 'right':
      #     shift = 1 - samples_per_window
      # elif align == 'center':
      #     shift = - (samples_per_window // 2)
      # else:
      #     shift = 0
      shift = -(samples_per_window // 2)
      idx = clear_windows.shift(shift)
      # drop rows at the end corresponding to windows past the end of data
      idx = idx.drop(clear_windows.index[1 - samples_per_window:])
      idx = idx.astype(bool)  # shift changed type to object
      clear_samples = np.unique(H[:, idx])
      return clear_samples
  
  
  def detect_clearsky(measured, clearsky, times=None, window_length=10,
                      mean_diff=75, max_diff=75,
                      lower_line_length=-5, upper_line_length=10,
                      var_diff=0.005, slope_dev=8, max_iterations=20,
                      return_components=False):
      """
      Detects clear sky times according to the algorithm developed by Reno
      and Hansen for GHI measurements. The algorithm [1]_ was designed and
      validated for analyzing GHI time series only. Users may attempt to
      apply it to other types of time series data using different filter
      settings, but should be skeptical of the results.
  
      The algorithm detects clear sky times by comparing statistics for a
      measured time series and an expected clearsky time series.
      Statistics are calculated using a sliding time window (e.g., 10
      minutes). An iterative algorithm identifies clear periods, uses the
      identified periods to estimate bias in the clearsky data, scales the
      clearsky data and repeats.
  
      Clear times are identified by meeting 5 criteria. Default values for
      these thresholds are appropriate for 10 minute windows of 1 minute
      GHI data.
  
      Parameters
      ----------
      measured : array or Series
          Time series of measured GHI. [W/m2]
      clearsky : array or Series
          Time series of the expected clearsky GHI. [W/m2]
      times : DatetimeIndex or None, default None.
          Times of measured and clearsky values. If None the index of measured
          will be used.
      window_length : int, default 10
          Length of sliding time window in minutes. Must be greater than 2
          periods.
      mean_diff : float, default 75
          Threshold value for agreement between mean values of measured
          and clearsky in each interval, see Eq. 6 in [1]. [W/m2]
      max_diff : float, default 75
          Threshold value for agreement between maxima of measured and
          clearsky values in each interval, see Eq. 7 in [1]. [W/m2]
      lower_line_length : float, default -5
          Lower limit of line length criterion from Eq. 8 in [1].
          Criterion satisfied when lower_line_length < line length difference
          < upper_line_length.
      upper_line_length : float, default 10
          Upper limit of line length criterion from Eq. 8 in [1].
      var_diff : float, default 0.005
          Threshold value in Hz for the agreement between normalized
          standard deviations of rate of change in irradiance, see Eqs. 9
          through 11 in [1].
      slope_dev : float, default 8
          Threshold value for agreement between the largest magnitude of
          change in successive values, see Eqs. 12 through 14 in [1].
      max_iterations : int, default 20
          Maximum number of times to apply a different scaling factor to
          the clearsky and redetermine clear_samples. Must be 1 or larger.
      return_components : bool, default False
          Controls if additional output should be returned. See below.
  
      Returns
      -------
      clear_samples : array or Series
          Boolean array or Series of whether or not the given time is
          clear. Return type is the same as the input type.
  
      components : OrderedDict, optional
          Dict of arrays of whether or not the given time window is clear
          for each condition. Only provided if return_components is True.
  
      alpha : scalar, optional
          Scaling factor applied to the clearsky_ghi to obtain the
          detected clear_samples. Only provided if return_components is
          True.
  
      Raises
      ------
      ValueError
          If measured is not a Series and times is not provided
      NotImplementedError
          If timestamps are not equally spaced
  
      References
      ----------
      .. [1] Reno, M.J. and C.W. Hansen, "Identification of periods of clear
         sky irradiance in time series of GHI measurements" Renewable Energy,
         v90, p. 520-531, 2016.
  
      Notes
      -----
      Initial implementation in MATLAB by Matthew Reno. Modifications for
      computational efficiency by Joshua Patrick and Curtis Martin. Ported
      to Python by Will Holmgren, Tony Lorenzo, and Cliff Hansen.
  
      Differences from MATLAB version:
  
          * no support for unequal times
          * automatically determines sample_interval
          * requires a reference clear sky series instead calculating one
            from a user supplied location and UTCoffset
          * parameters are controllable via keyword arguments
          * option to return individual test components and clearsky scaling
            parameter
          * uses centered windows (Matlab function uses left-aligned windows)
      """
  
      if times is None:
          try:
              times = measured.index
          except AttributeError:
              raise ValueError("times is required when measured is not a Series")
  
      # be polite about returning the same type as was input
      ispandas = isinstance(measured, pd.Series)
  
      # for internal use, need a Series
      if not ispandas:
          meas = pd.Series(measured, index=times)
      else:
          meas = measured
  
      if not isinstance(clearsky, pd.Series):
          clear = pd.Series(clearsky, index=times)
      else:
          clear = clearsky
  
      sample_interval, samples_per_window = _get_sample_intervals(times,
                                                                  window_length)
  
      # generate matrix of integers for creating windows with indexing
      H = hankel(np.arange(samples_per_window),
                 np.arange(samples_per_window-1, len(times)))
  
      # calculate measurement statistics
      meas_mean, meas_max, meas_slope_nstd, meas_slope = _calc_stats(
          meas, samples_per_window, sample_interval, H)
      meas_line_length = _line_length_windowed(
          meas, H, samples_per_window, sample_interval)
  
      # calculate clear sky statistics
      clear_mean, clear_max, _, clear_slope = _calc_stats(
          clear, samples_per_window, sample_interval, H)
  
      # find a scaling factor for the clear sky time series that minimizes the
      # RMSE between the clear times identified in the measured data and the
      # scaled clear sky time series. Optimization to determine the scaling
      # factor considers all identified clear times, which is different from [1]
      # where the scaling factor was determined from clear times on days with
      # at least 50% of the day being identified as clear.
      alpha = 1
      for iteration in range(max_iterations):
          scaled_clear = alpha * clear
          clear_line_length = _line_length_windowed(
              scaled_clear, H, samples_per_window, sample_interval)
  
          line_diff = meas_line_length - clear_line_length
          slope_max_diff = _max_diff_windowed(
              meas - scaled_clear, H, samples_per_window)
          # evaluate comparison criteria
          c1 = np.abs(meas_mean - alpha*clear_mean) < mean_diff
          c2 = np.abs(meas_max - alpha*clear_max) < max_diff
          c3 = (line_diff > lower_line_length) & (line_diff < upper_line_length)
          c4 = meas_slope_nstd < var_diff
          c5 = slope_max_diff < slope_dev
          c6 = (clear_mean != 0) & ~np.isnan(clear_mean)
          clear_windows = c1 & c2 & c3 & c4 & c5 & c6
  
          # create array to return
          clear_samples = np.full_like(meas, False, dtype='bool')
          # find the samples contained in any window classified as clear
          idx = _clear_sample_index(clear_windows, samples_per_window, 'center',
                                    H)
          clear_samples[idx] = True
  
          # find a new alpha
          previous_alpha = alpha
          clear_meas = meas[clear_samples]
          clear_clear = clear[clear_samples]
  
          def rmse(alpha):
              return np.sqrt(np.mean((clear_meas - alpha*clear_clear)**2))
  
          alpha = minimize_scalar(rmse).x
          if round(alpha*10000) == round(previous_alpha*10000):
              break
      else:
          import warnings
          warnings.warn('rescaling failed to converge after %s iterations'
                        % max_iterations, RuntimeWarning)
  
      # be polite about returning the same type as was input
      if ispandas:
          clear_samples = pd.Series(clear_samples, index=times)
  
      if return_components:
          components = OrderedDict()
          components['mean_diff_flag'] = c1
          components['max_diff_flag'] = c2
          components['line_length_flag'] = c3
          components['slope_nstd_flag'] = c4
          components['slope_max_flag'] = c5
          components['mean_nan_flag'] = c6
          components['windows'] = clear_windows
  
          components['mean_diff'] = np.abs(meas_mean - alpha * clear_mean)
          components['max_diff'] = np.abs(meas_max - alpha * clear_max)
          components['line_length'] = meas_line_length - clear_line_length
          components['slope_nstd'] = meas_slope_nstd
          components['slope_max'] = slope_max_diff
  
          return clear_samples, components, alpha
      else:
          return clear_samples
  
  
  def bird(zenith, airmass_relative, aod380, aod500, precipitable_water,
           ozone=0.3, pressure=101325., dni_extra=1364., asymmetry=0.85,
           albedo=0.2):
      """
      Bird Simple Clear Sky Broadband Solar Radiation Model
  
      Based on NREL Excel implementation by Daryl R. Myers [1, 2].
  
      Bird and Hulstrom define the zenith as the "angle between a line to
      the sun and the local zenith". There is no distinction in the paper
      between solar zenith and apparent (or refracted) zenith, but the
      relative airmass is defined using the Kasten 1966 expression, which
      requires apparent zenith. Although the formulation for calculated
      zenith is never explicitly defined in the report, since the purpose
      was to compare existing clear sky models with "rigorous radiative
      transfer models" (RTM) it is possible that apparent zenith was
      obtained as output from the RTM. However, the implentation presented
      in PVLIB is tested against the NREL Excel implementation by Daryl
      Myers which uses an analytical expression for solar zenith instead
      of apparent zenith.
  
      Parameters
      ----------
      zenith : numeric
          Solar or apparent zenith angle in degrees - see note above
      airmass_relative : numeric
          Relative airmass
      aod380 : numeric
          Aerosol optical depth [cm] measured at 380[nm]
      aod500 : numeric
          Aerosol optical depth [cm] measured at 500[nm]
      precipitable_water : numeric
          Precipitable water [cm]
      ozone : numeric
          Atmospheric ozone [cm], defaults to 0.3[cm]
      pressure : numeric
          Ambient pressure [Pa], defaults to 101325[Pa]
      dni_extra : numeric
          Extraterrestrial radiation [W/m^2], defaults to 1364[W/m^2]
      asymmetry : numeric
          Asymmetry factor, defaults to 0.85
      albedo : numeric
          Albedo, defaults to 0.2
  
      Returns
      -------
      clearsky : DataFrame (if Series input) or OrderedDict of arrays
          DataFrame/OrderedDict contains the columns/keys
          ``'dhi', 'dni', 'ghi', 'direct_horizontal'`` in  [W/m^2].
  
      See also
      --------
      pvlib.atmosphere.bird_hulstrom80_aod_bb
      pvlib.atmosphere.get_relative_airmass
  
      References
      ----------
      .. [1] R. E. Bird and R. L Hulstrom, "A Simplified Clear Sky model for
         Direct and Diffuse Insolation on Horizontal Surfaces" SERI Technical
         Report SERI/TR-642-761, Feb 1981. Solar Energy Research Institute,
         Golden, CO.
  
      .. [2] Daryl R. Myers, "Solar Radiation: Practical Modeling for Renewable
         Energy Applications", pp. 46-51 CRC Press (2013)
  
      .. [3] `NREL Bird Clear Sky Model <http://rredc.nrel.gov/solar/models/
         clearsky/>`_
  
      .. [4] `SERI/TR-642-761 <http://rredc.nrel.gov/solar/pubs/pdfs/
         tr-642-761.pdf>`_
  
      .. [5] `Error Reports <http://rredc.nrel.gov/solar/models/clearsky/
         error_reports.html>`_
      """
      etr = dni_extra  # extraradiation
      ze_rad = np.deg2rad(zenith)  # zenith in radians
      airmass = airmass_relative
      # Bird clear sky model
      am_press = atmosphere.get_absolute_airmass(airmass, pressure)
      t_rayleigh = (
          np.exp(-0.0903 * am_press ** 0.84 * (
              1.0 + am_press - am_press ** 1.01
          ))
      )
      am_o3 = ozone*airmass
      t_ozone = (
          1.0 - 0.1611 * am_o3 * (1.0 + 139.48 * am_o3) ** -0.3034 -
          0.002715 * am_o3 / (1.0 + 0.044 * am_o3 + 0.0003 * am_o3 ** 2.0)
      )
      t_gases = np.exp(-0.0127 * am_press ** 0.26)
      am_h2o = airmass * precipitable_water
      t_water = (
          1.0 - 2.4959 * am_h2o / (
              (1.0 + 79.034 * am_h2o) ** 0.6828 + 6.385 * am_h2o
          )
      )
      bird_huldstrom = atmosphere.bird_hulstrom80_aod_bb(aod380, aod500)
      t_aerosol = np.exp(
          -(bird_huldstrom ** 0.873) *
          (1.0 + bird_huldstrom - bird_huldstrom ** 0.7088) * airmass ** 0.9108
      )
      taa = 1.0 - 0.1 * (1.0 - airmass + airmass ** 1.06) * (1.0 - t_aerosol)
      rs = 0.0685 + (1.0 - asymmetry) * (1.0 - t_aerosol / taa)
      id_ = 0.9662 * etr * t_aerosol * t_water * t_gases * t_ozone * t_rayleigh
      ze_cos = np.where(zenith < 90, np.cos(ze_rad), 0.0)
      id_nh = id_ * ze_cos
      ias = (
          etr * ze_cos * 0.79 * t_ozone * t_gases * t_water * taa *
          (0.5 * (1.0 - t_rayleigh) + asymmetry * (1.0 - (t_aerosol / taa))) / (
              1.0 - airmass + airmass ** 1.02
          )
      )
      gh = (id_nh + ias) / (1.0 - albedo * rs)
      diffuse_horiz = gh - id_nh
      # TODO: be DRY, use decorator to wrap methods that need to return either
      # OrderedDict or DataFrame instead of repeating this boilerplate code
      irrads = OrderedDict()
      irrads['direct_horizontal'] = id_nh
      irrads['ghi'] = gh
      irrads['dni'] = id_
      irrads['dhi'] = diffuse_horiz
      if isinstance(irrads['dni'], pd.Series):
          irrads = pd.DataFrame.from_dict(irrads)
      return irrads
  