Repository: pvlib/pvlib-python
Pull request #1074 (2020-10-05 16:43:18)

Refactor detect_clearsky
 - [x] I am familiar with the [contributing guidelines](https://pvlib-python.readthedocs.io/en/latest/contributing.html)
 - [x] Tests added
 - [x] Adds description and name entries in the appropriate "what's new" file in [`docs/sphinx/source/whatsnew`](https://github.com/pvlib/pvlib-python/tree/master/docs/sphinx/source/whatsnew) for all changes. Includes link to the GitHub Issue with `` :issue:`num` `` or this Pull Request with `` :pull:`num` ``. Includes contributor name and/or GitHub username (link with `` :ghuser:`user` ``).
 - [x] New code is fully documented. Includes [numpydoc](https://numpydoc.readthedocs.io/en/latest/format.html) compliant docstrings, examples, and comments where necessary.
 - [x] Pull request is nearly complete and ready for detailed review.
 - [x] Maintainer: Appropriate GitHub Labels and Milestone are assigned to the Pull Request and linked Issue.

Refactor of `clearsky.detect_clearsky` in preparation for adding BrightSun detect clearsky algorithm. Major changes:
- uses pd.Series.rolling instead of Hankel matrix to compute criteria over sliding windows. This change accomodates centered windows, i.e., aligning the criteria for a window to a midpoint rather than to the left endpoint (default behavior in this code).
- calculation of criteria delegated to helper functions.
- new kwarg `align` (default 'left') otherwise function signature is unchanged.


cwhanse: if `clearsky.detect_clearsky` is of interest please take a look, in particular, at `_calc_stats` and the other helper functions.

I need to add some tests for the non-trivial helper functions, and I'm sure there's room for performance improvement. The `align` feature is added anticipating the BrightSun method which uses centered windows. I think that the BrightSun method could be implemented using numpy arrays and the Hankel matrix, but using pandas.Series.rolling seems more transparent to me.

We could drop `align='right'` as an option (the current code's behavior is `align='left'` but `pandas.rolling` only offers `center=False (default, means right-aligned) and `center=True` (centered). Maintaining the current behavior with `align='left'` by default wasn't too hard. I'm not sure how much value the `align` options offer, but there they are. `align='center'` is more in line (sorry) with what Matt and I had in mind when we wrote the algorithm, but we didn't code it that way in Matlab.
wholmgren: Is there a compelling use case for `align='left'`? It sounds to me like we should consider it a bug in the initial implementation, and always use `center=True` (without exposing a kwarg). 

I will carefully review when I have the time. I want to understand why this implementation works but the previous attempts to use pandas did not work.
cwhanse: `align='left'` is here to preserve the old behavior, personally I agree with `align='center'` as the better choice for a default.

I also think we could drop support for `array` inputs and lose the `times` argument.
wholmgren: To put it differently: why should we preserve the old behavior? It sounds like a bug. People that want the old behavior can use an old, buggier version of pvlib.

> I also think we could drop support for array inputs and lose the times argument.

agree
cwhanse: I guess I'm unwilling to call the left alignment a bug - it's not wrong, rather, it's an undesirable behavior. I'm OK changing the default to 'center'. And open to not offering 'left', 'right' as options.
wholmgren: Other than consistency with previous versions, when might a user prefer 'left' or 'right'? Without a compelling reason, I'd hard-code center and then let users ask/contribute if they come up with a good reason for a kwarg.
cwhanse: > Other than consistency with previous versions, when might a user prefer 'left' or 'right'? 

Consistency is the only reason I came up with. The option to align 'left', 'center', 'right' falls out of the combination of numpy and pandas defaults.

> I'd hard-code center 

I'm OK with this.


cwhanse: > To put it differently: why should we preserve the old behavior? It sounds like a bug. People that want the old behavior can use an old, buggier version of pvlib.
> 
> > I also think we could drop support for array inputs and lose the times argument.
> 
> agree

@wholmgren prefer to deprecate removing the `times` argument and restricting input type to Series?
wholmgren: Ok with me but how do we gracefully deprecate a required argument from the middle of a signature?
cwhanse: > > @wholmgren prefer to deprecate removing the `times` argument and restricting input type to Series?
> 
> Ok with me but how do we gracefully deprecate a required argument from the middle of a signature? I think in the past we've only done this kind of change when moving a function to a new name or module.


"gracefully" idk, maybe accept array types and the `times` argument, and use the `times` argument but raise a warning if it is supplied.

Or release it as a breaking API change.

My guess is that the array-type users are in the minority and the Series-type users would appreciate losing the `times` parameter.
wholmgren: I'm not opposed to delaying this to 0.9 and making it a breaking change. We can do a non-breaking change if we change both `times` and `window_length` to kwargs. Maybe someone else has a better idea.
cwhanse: Let's leave the parameters unchanged. In v0.9 I'd want to rename this function to `detect_clearsky_reno` to make room for `detect_clearsky_brightsun`. We can adjust the parameters in the renamed function.
cwhanse: > * Consider dropping the newly unnecessary `times` argument on line 643 in `clearsky.rst`.

We decided not to drop the 'times' argument in this refactor, but rather to rename this function in v0.9 to `detect_clearsky_reno` and drop the 'times' argument then.

> * I notice that the reference performed the rescaling "using only clear periods on days when greater than 50% of the day was classified as clear".  Just pointing it out in case it's worth doing the same here or documenting the difference.

Adding a comment to this effect.
kanderso-nrel: > 
> 
> > * Consider dropping the newly unnecessary `times` argument on line 643 in `clearsky.rst`.
> 
> We decided not to drop the 'times' argument in this refactor, but rather to rename this function in v0.9 to `detect_clearsky_reno` and drop the 'times' argument then.

I'm on board with that plan for 0.9 -- to clarify, this suggestion was about this [documentation page](https://pvlib-python.readthedocs.io/en/stable/clearsky.html#detect-clearsky) where the inputs are Series and passing `cs.index` isn't needed anymore.  OK with me to just leave it until 0.9 when the parameter goes away for good. 
wholmgren: https://github.com/cwhanse/pvlib-python/pull/6 adds an asv test for this function. Results on my machine for the latest commit 5927d73 using `asv run -b 'detect_clearsky.DetectClear.time_detect_clearsky' HEAD`

```
[  0.14%] · For pvlib-python commit 5927d738 <pr/cwhanse/1074>:
[  0.14%] ·· Building for conda-py3.6-ephem3.7.6.0-numba0.36.1-numpy1.12.0-pandas0.22.0-pytables3.6.1-scipy1.2.0....
[  0.14%] ·· Benchmarking conda-py3.6-ephem3.7.6.0-numba0.36.1-numpy1.12.0-pandas0.22.0-pytables3.6.1-scipy1.2.0
[  0.17%] ··· Running (detect_clearsky.DetectClear.time_detect_clearsky--).
[  0.21%] ··· detect_clearsky.DetectClear.time_detect_clearsky                                                                                                                                                                      1.27±0.02s
[  0.21%] ·· Building for conda-py3.8-ephem-numba-numpy-pandas-pytables-scipy....
[  0.21%] ·· Benchmarking conda-py3.8-ephem-numba-numpy-pandas-pytables-scipy
[  0.24%] ··· Running (detect_clearsky.DetectClear.time_detect_clearsky--).
[  0.28%] ··· detect_clearsky.DetectClear.time_detect_clearsky                                                                                                                                                                      13.5±0.03s
```

disconcerting result for the python 3.8 environment.
cwhanse: thanks @wholmgren  I profiled the code and nearly all the time is being spent in those two '.apply' methods. That's the good news, I know where to focus. I don't think this should be merged with such a drastic loss of performance. 
kanderso-nrel: I would guess that the order of magnitude slowdown between the two asv environments is from pandas introducing the optional `raw` parameter in v0.23.0, making `Rolling.apply()` switch from passing arrays to passing Series by default. 
wholmgren: > I would guess that the order of magnitude slowdown between the two asv environments is from pandas introducing the optional raw parameter in v0.23.0, making Rolling.apply() switch from passing arrays to passing Series by default.

That's a good guess. If true, we should bump the min pandas requirement so that we can pass `raw=True` without concern for the version. 

@cwhanse I'm fine with holding up 0.8.1 for a little bit while you experiment. I just want it released by the time we'd want to merge the Array PR.

Patch:
diff --git a/benchmarks/benchmarks/detect_clearsky.py b/benchmarks/benchmarks/detect_clearsky.py
new file mode 100644
index 0000000..7f2039f
--- /dev/null
+++ b/benchmarks/benchmarks/detect_clearsky.py
@@ -0,0 +1,34 @@
+"""
+ASV benchmarks for detect clear sky function.
+"""
+
+import pandas as pd
+from pvlib import clearsky, solarposition
+import numpy as np
+
+
+class DetectClear:
+    params = [1, 10, 100]  # number of days
+    param_names = ['ndays']
+
+    def setup(self, ndays):
+        self.times = pd.date_range(start='20180601', freq='1min',
+                                   periods=1440*ndays)
+        self.lat = 35.1
+        self.lon = -106.6
+        self.solar_position = solarposition.get_solarposition(
+            self.times, self.lat, self.lon)
+        clearsky_df = clearsky.simplified_solis(
+            self.solar_position['apparent_elevation'])
+        self.clearsky = clearsky_df['ghi']
+        measured_dni = clearsky_df['dni'].where(
+            (self.times.hour % 2).astype(bool), 0)
+        cos_zen = np.cos(np.deg2rad(self.solar_position['apparent_zenith']))
+        self.measured = measured_dni * cos_zen + clearsky_df['dhi']
+        self.measured *= 0.98
+        self.window_length = 10
+
+    def time_detect_clearsky(self, ndays):
+        clearsky.detect_clearsky(
+            self.measured, self.clearsky, self.times, self.window_length
+        )
diff --git a/docs/examples/plot_spectrl2_fig51A.py b/docs/examples/plot_spectrl2_fig51A.py
new file mode 100644
index 0000000..6f3470a
--- /dev/null
+++ b/docs/examples/plot_spectrl2_fig51A.py
@@ -0,0 +1,97 @@
+"""
+Modeling Spectral Irradiance
+============================
+
+Recreating Figure 5-1A from the SPECTRL2 NREL Technical Report.
+"""
+
+# %%
+# This example shows how to model the spectral distribution of irradiance
+# based on atmospheric conditions. The spectral distribution of irradiance is
+# the power content at each wavelength band in the solar spectrum and is
+# affected by various scattering and absorption mechanisms in the atmosphere.
+# This example recreates an example figure from the SPECTRL2 NREL Technical
+# Report [1]_. The figure shows modeled spectra at hourly intervals across
+# a single morning.
+#
+# References
+# ----------
+# .. [1] Bird, R, and Riordan, C., 1984, "Simple solar spectral model for
+#    direct and diffuse irradiance on horizontal and tilted planes at the
+#    earth's surface for cloudless atmospheres", NREL Technical Report
+#    TR-215-2436 doi:10.2172/5986936.
+
+# %%
+# The SPECTRL2 model has several inputs; some can be calculated with pvlib,
+# but other must come from a weather dataset. In this case, these weather
+# parameters are example assumptions taken from the technical report.
+
+from pvlib import spectrum, solarposition, irradiance, atmosphere
+import pandas as pd
+import matplotlib.pyplot as plt
+
+# assumptions from the technical report:
+lat = 37
+lon = -100
+tilt = 37
+azimuth = 180
+pressure = 101300  # sea level, roughly
+water_vapor_content = 0.5  # cm
+tau500 = 0.1
+ozone = 0.31  # atm-cm
+albedo = 0.2
+
+times = pd.date_range('1984-03-20 06:17', freq='h', periods=6, tz='Etc/GMT+7')
+solpos = solarposition.get_solarposition(times, lat, lon)
+aoi = irradiance.aoi(tilt, azimuth, solpos.apparent_zenith, solpos.azimuth)
+
+# The technical report uses the 'kasten1966' airmass model, but later
+# versions of SPECTRL2 use 'kastenyoung1989'.  Here we use 'kasten1966'
+# for consistency with the technical report.
+relative_airmass = atmosphere.get_relative_airmass(solpos.apparent_zenith,
+                                                   model='kasten1966')
+
+# %%
+# With all the necessary inputs in hand we can model spectral irradiance using
+# :py:func:`pvlib.spectrum.spectrl2`.  Note that because we are calculating
+# the spectra for more than one set of conditions, we will get back 2-D
+# arrays (one dimension for wavelength, one for time).
+
+spectra = spectrum.spectrl2(
+    apparent_zenith=solpos.apparent_zenith,
+    aoi=aoi,
+    surface_tilt=tilt,
+    ground_albedo=albedo,
+    surface_pressure=pressure,
+    relative_airmass=relative_airmass,
+    precipitable_water=water_vapor_content,
+    ozone=ozone,
+    aerosol_turbidity_500nm=tau500,
+)
+
+# %%
+# The ``poa_global`` array represents the total spectral irradiance on our
+# hypothetical solar panel. Let's plot it against wavelength to recreate
+# Figure 5-1A:
+
+plt.figure()
+plt.plot(spectra['wavelength'], spectra['poa_global'])
+plt.xlim(200, 2700)
+plt.ylim(0, 1.8)
+plt.title(r"Day 80 1984, $\tau=0.1$, Wv=0.5 cm")
+plt.ylabel(r"Irradiance ($W m^{-2} nm^{-1}$)")
+plt.xlabel(r"Wavelength ($nm$)")
+time_labels = times.strftime("%H:%M %p")
+labels = [
+    "AM {:0.02f}, Z{:0.02f}, {}".format(*vals)
+    for vals in zip(relative_airmass, solpos.apparent_zenith, time_labels)
+]
+plt.legend(labels)
+plt.show()
+
+# %%
+# Note that the airmass and zenith values do not exactly match the values in
+# the technical report; this is because airmass is estimated from solar
+# position and the solar position calculation in the technical report does not
+# exactly match the one used here.  However, the differences are minor enough
+# to not materially change the spectra.
diff --git a/docs/sphinx/Makefile b/docs/sphinx/Makefile
index aedab24..0c96b4e 100644
--- a/docs/sphinx/Makefile
+++ b/docs/sphinx/Makefile
@@ -48,6 +48,9 @@ help:
 
 clean:
 	rm -rf $(BUILDDIR)/*
+	rm -rf source/generated
+	rm -rf source/auto_examples
+	rm -rf source/savefig
 
 html:
 	$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
diff --git a/docs/sphinx/make.bat b/docs/sphinx/make.bat
index 9534b01..d12b216 100644
--- a/docs/sphinx/make.bat
+++ b/docs/sphinx/make.bat
@@ -12,6 +12,17 @@ set BUILDDIR=build
 
 if "%1" == "" goto help
 
+
+if "%1" == "clean" (
+    REM override the default `make clean` behavior of sphinx-build;
+    REM this lets us clean out the various build files in sphinx/source/
+	for /d %%i in (%BUILDDIR%\*) do rmdir /q /s %%i
+	rmdir /q /s %SOURCEDIR%\generated >nul 2>&1
+	rmdir /q /s %SOURCEDIR%\auto_examples >nul 2>&1
+	rmdir /q /s %SOURCEDIR%\savefig >nul 2>&1
+	goto end
+)
+
 %SPHINXBUILD% >NUL 2>NUL
 if errorlevel 9009 (
 	echo.
diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 7db1b09..2ec4e3f 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -417,6 +417,13 @@ Shading
    shading.masking_angle_passias
    shading.sky_diffuse_passias
 
+Spectrum
+--------
+
+.. autosummary::
+   :toctree: generated/
+
+   spectrum.spectrl2
 
 Tracking
 ========
diff --git a/docs/sphinx/source/whatsnew.rst b/docs/sphinx/source/whatsnew.rst
index 0e71fd7..271c898 100644
--- a/docs/sphinx/source/whatsnew.rst
+++ b/docs/sphinx/source/whatsnew.rst
@@ -6,6 +6,7 @@ What's New
 
 These are new features and improvements of note in each release.
 
+.. include:: whatsnew/v0.8.1.rst
 .. include:: whatsnew/v0.8.0.rst
 .. include:: whatsnew/v0.7.2.rst
 .. include:: whatsnew/v0.7.1.rst
diff --git a/docs/sphinx/source/whatsnew/v0.8.1.rst b/docs/sphinx/source/whatsnew/v0.8.1.rst
index 0c1b588..4e80d00 100644
--- a/docs/sphinx/source/whatsnew/v0.8.1.rst
+++ b/docs/sphinx/source/whatsnew/v0.8.1.rst
@@ -6,17 +6,25 @@ v0.8.1 (MONTH DAY YEAR)
 Breaking changes
 ~~~~~~~~~~~~~~~~
 
-
 Deprecations
 ~~~~~~~~~~~~
-* ``pvlib.irradiance.liujordan`` is deprecated.
+* ``pvlib.irradiance.liujordan`` is deprecated. :py:func:`pvlib.irradiance.campbell_norman`
+  replaces ``pvlib.irradiance.liujordan``.
 
 Enhancements
 ~~~~~~~~~~~~
+* Add a numpy-based implementation of the SPECTRL2 spectral irradiance model
+  :py:func:`pvlib.spectrum.spectrl2` (:pull:`1062`)
 * Create :py:func:`~pvlib.pvsystem.PVSystem.fuentes_celltemp` and add ``temperature_model='fuentes'``
   option to :py:class:`~pvlib.modelchain.ModelChain`. (:pull:`1042`) (:issue:`1073`)
 * Added :py:func:`pvlib.temperature.ross` for cell temperature modeling using
   only NOCT. (:pull:`1045`)
+* :py:func:`pvlib.clearsky.detect_clearsky` now uses centered rolling windows
+  instead of left-aligned rolling windows. (:pull:`1074`)
+* The 'times' and 'window_length` parameters are now optional kwargs in
+  :py:func:`pvlib.clearsky.detect_clearsky`. If omitted, 'times' is set
+  equal to the index of parameter 'measured', and 'window_length' is set to
+  10 minutes. (:pull:`1074`)
 * Added :py:func:`pvlib.inverter.sandia_multi` and :py:func:`pvlib.inverter.pvwatts_multi`
   for modeling inverters with multiple MPPTs (:issue:`457`, :pull:`1085`, :pull:`1106`)
 * Added optional ``attributes`` parameter to :py:func:`pvlib.iotools.get_psm3`
@@ -42,6 +50,8 @@ Testing
 * Add airspeed velocity performance testing configuration and a few benchmarks.
   (:issue:`419`, :pull:`1049`, :pull:`1059`)
 * Add Python 3.9 CI configurations. (:issue:`1102`, :pull:`1112`)
+* Update ``test_pvgis.py`` to be more flexible about the PVGIS copyright notice
+  (:pull:`1121`)
 
 Documentation
 ~~~~~~~~~~~~~
@@ -49,6 +59,7 @@ Documentation
   (:issue:`1055`, :pull:`1075`)
 * Add gallery example about backtracking on sloped terrain. (:pull:`1077`)
 * Add toggle button for code prompts to make copying code easier (:pull:`1096`)
+* Add project urls to setup.py for pypi page (:pull:`1119`)
 
 Requirements
 ~~~~~~~~~~~~
@@ -62,3 +73,5 @@ Contributors
 * Cliff Hansen (:ghuser:`cwhanse`)
 * Will Vining (:ghuser:`wfvining`)
 * Michael Jurasovic (:ghuser:`jurasofish`)
+* Aziz Ben Othman (:ghuser:`AzizCode92`)
+
diff --git a/pvlib/__init__.py b/pvlib/__init__.py
index 55ffa51..ff6b375 100644
--- a/pvlib/__init__.py
+++ b/pvlib/__init__.py
@@ -20,6 +20,7 @@ from pvlib import (  # noqa: F401
     soiling,
     solarposition,
     spa,
+    spectrum,
     temperature,
     tools,
     tracking,
diff --git a/pvlib/clearsky.py b/pvlib/clearsky.py
index a6794e2..857c5dd 100644
--- a/pvlib/clearsky.py
+++ b/pvlib/clearsky.py
@@ -9,6 +9,8 @@ import calendar
 
 import numpy as np
 import pandas as pd
+from scipy.optimize import minimize_scalar
+from scipy.linalg import hankel
 
 from pvlib import atmosphere, tools
 
@@ -597,7 +599,118 @@ def _calc_d(aod700, p):
     return d
 
 
-def detect_clearsky(measured, clearsky, times, window_length,
+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
+
+    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
+    """
+
+    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, H, samples_per_window)
+    data_slope_nstd = data_slope_nstd
+
+    return data_mean, data_max, data_slope_nstd, data_slope
+
+
+def _slope_nstd_windowed(data, H, samples_per_window):
+    with np.errstate(divide='ignore', invalid='ignore'):
+        raw = np.diff(data)
+        raw = raw[H[:-1, ]].std(ddof=1, axis=0) / data.values[H].mean(axis=0)
+    return _to_centered_series(raw, 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,
@@ -623,24 +736,25 @@ def detect_clearsky(measured, clearsky, times, window_length,
     Parameters
     ----------
     measured : array or Series
-        Time series of measured values.
+        Time series of measured GHI. [W/m2]
     clearsky : array or Series
-        Time series of the expected clearsky values.
-    times : DatetimeIndex
-        Times of measured and clearsky values.
-    window_length : int
+        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].
+        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].
+        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
+        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
@@ -671,6 +785,13 @@ def detect_clearsky(measured, clearsky, times, window_length,
         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
@@ -692,78 +813,81 @@ def detect_clearsky(measured, clearsky, times, window_length,
         * 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)
     """
 
-    # calculate deltas in units of minutes (matches input window_length units)
-    deltas = np.diff(times.values) / np.timedelta64(1, '60s')
+    if times is None:
+        try:
+            times = measured.index
+        except AttributeError:
+            raise ValueError("times is required when measured is not a Series")
 
-    # determine the unique deltas and if we can proceed
-    unique_deltas = np.unique(deltas)
-    if len(unique_deltas) == 1:
-        sample_interval = unique_deltas[0]
+    # 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:
-        raise NotImplementedError('algorithm does not yet support unequal '
-                                  'times. consider resampling your data.')
+        meas = measured
 
-    intervals_per_window = int(window_length / sample_interval)
+    if not isinstance(clearsky, pd.Series):
+        clear = pd.Series(clearsky, index=times)
+    else:
+        clear = clearsky
 
-    # generate matrix of integers for creating windows with indexing
-    from scipy.linalg import hankel
-    H = hankel(np.arange(intervals_per_window),                   # noqa: N806
-               np.arange(intervals_per_window - 1, len(times)))
+    sample_interval, samples_per_window = _get_sample_intervals(times,
+                                                                window_length)
 
-    # convert pandas input to numpy array, but save knowledge of input state
-    # so we can return a series if that's what was originally provided
-    ispandas = isinstance(measured, pd.Series)
-    measured = np.asarray(measured)
-    clearsky = np.asarray(clearsky)
+    # 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 = np.mean(measured[H], axis=0)
-    meas_max = np.max(measured[H], axis=0)
-    meas_diff = np.diff(measured[H], n=1, axis=0)
-    meas_slope = np.diff(measured[H], n=1, axis=0) / sample_interval
-    # matlab std function normalizes by N-1, so set ddof=1 here
-    meas_slope_nstd = np.std(meas_slope, axis=0, ddof=1) / meas_mean
-    meas_line_length = np.sum(np.sqrt(
-        meas_diff * meas_diff +
-        sample_interval * sample_interval), axis=0)
+    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 = np.mean(clearsky[H], axis=0)
-    clear_max = np.max(clearsky[H], axis=0)
-    clear_diff = np.diff(clearsky[H], n=1, axis=0)
-    clear_slope = np.diff(clearsky[H], n=1, axis=0) / sample_interval
-
-    from scipy.optimize import minimize_scalar
-
+    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):
-        clear_line_length = np.sum(np.sqrt(
-            alpha * alpha * clear_diff * clear_diff +
-            sample_interval * sample_interval), axis=0)
+        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 = np.max(np.abs(meas_slope -
-                           alpha * clear_slope), axis=0) < slope_dev
+        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(measured, False, dtype='bool')
+        clear_samples = np.full_like(meas, False, dtype='bool')
         # find the samples contained in any window classified as clear
-        clear_samples[np.unique(H[:, clear_windows])] = True
+        idx = _clear_sample_index(clear_windows, samples_per_window, 'center',
+                                  H)
+        clear_samples[idx] = True
 
         # find a new alpha
         previous_alpha = alpha
-        clear_meas = measured[clear_samples]
-        clear_clear = clearsky[clear_samples]
+        clear_meas = meas[clear_samples]
+        clear_clear = clear[clear_samples]
 
         def rmse(alpha):
             return np.sqrt(np.mean((clear_meas - alpha*clear_clear)**2))
@@ -773,7 +897,7 @@ def detect_clearsky(measured, clearsky, times, window_length,
             break
     else:
         import warnings
-        warnings.warn('failed to converge after %s iterations'
+        warnings.warn('rescaling failed to converge after %s iterations'
                       % max_iterations, RuntimeWarning)
 
     # be polite about returning the same type as was input
@@ -794,8 +918,7 @@ def detect_clearsky(measured, clearsky, times, window_length,
         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'] = (np.max(
-            meas_slope - alpha * clear_slope, axis=0))
+        components['slope_max'] = slope_max_diff
 
         return clear_samples, components, alpha
     else:
diff --git a/pvlib/data/spectrl2_example_spectra.csv b/pvlib/data/spectrl2_example_spectra.csv
new file mode 100644
index 0000000..7c1d8b1
--- /dev/null
+++ b/pvlib/data/spectrl2_example_spectra.csv
@@ -0,0 +1,123 @@
+,wavelength,specdif,specdir,specetr,specglo
+0,0.30000001192092896,0.7665966153144836,0.40335652232170105,541.6846923828125,1.036954402923584
+1,0.3050000071525574,11.298118591308594,6.824652671813965,564.326416015625,15.872478485107422
+2,0.3100000023841858,36.487388610839844,25.05140495300293,628.7140502929688,53.278594970703125
+3,0.3149999976158142,80.22676849365234,62.04690170288086,700.1771850585938,121.81494903564453
+4,0.3199999928474426,108.57009887695312,93.82110595703125,722.8189697265625,171.45558166503906
+5,0.32499998807907104,155.0250701904297,148.59210205078125,841.8905639648438,254.62191772460938
+6,0.33000001311302185,198.05923461914062,209.17189025878906,972.2830200195312,338.2608947753906
+7,0.33500000834465027,198.03289794921875,229.05174255371094,941.959228515625,351.5594482421875
+8,0.3400000035762787,192.4607391357422,242.4566650390625,910.3212890625,354.97216796875
+9,0.3449999988079071,194.3870391845703,265.3816833496094,921.1368408203125,372.26446533203125
+10,0.3499999940395355,206.6869659423828,304.3901062011719,986.0298461914062,410.7105712890625
+11,0.36000001430511475,202.11509704589844,342.1354675292969,986.4341430664062,431.4382629394531
+12,0.3700000047683716,225.76776123046875,433.04852294921875,1131.988525390625,516.0272216796875
+13,0.3799999952316284,216.52438354492188,464.911865234375,1115.7147216796875,528.140869140625
+14,0.38999998569488525,197.29737854003906,469.2872619628906,1044.959228515625,511.84661865234375
+15,0.4000000059604645,274.6373291015625,717.1521606445312,1495.0657958984375,755.322998046875
+16,0.4099999964237213,307.39105224609375,874.312255859375,1719.664306640625,893.4163818359375
+17,0.41999998688697815,306.0761413574219,941.7630004882812,1759.1864013671875,937.311767578125
+18,0.4300000071525574,271.7977294921875,899.2123413085938,1604.3326416015625,874.5128784179688
+19,0.4399999976158142,306.446533203125,1084.28076171875,1856.8291015625,1033.20751953125
+20,0.44999998807907104,325.4994812011719,1225.8232421875,2026.642578125,1147.132080078125
+21,0.46000000834465027,317.2213134765625,1288.9774169921875,2065.052734375,1181.184326171875
+22,0.4699999988079071,295.31744384765625,1289.505615234375,2008.4482421875,1159.634521484375
+23,0.47999998927116394,288.3025817871094,1347.8497314453125,2048.880126953125,1191.7259521484375
+24,0.49000000953674316,258.0177917480469,1287.21337890625,1916.4659423828125,1120.7984619140625
+25,0.5,248.51402282714844,1318.9627685546875,1929.6063232421875,1132.5753173828125
+26,0.5099999904632568,240.04953002929688,1351.5888671875,1947.800537109375,1145.9791259765625
+27,0.5199999809265137,218.66915893554688,1302.779541015625,1850.7642822265625,1091.88330078125
+28,0.5299999713897705,215.91168212890625,1357.883544921875,1911.4119873046875,1126.0604248046875
+29,0.5400000214576721,207.6764678955078,1375.676513671875,1918.487548828125,1129.7513427734375
+30,0.550000011920929,198.75257873535156,1383.860595703125,1912.4227294921875,1126.31298828125
+31,0.5699999928474426,177.2172393798828,1355.662841796875,1859.8614501953125,1085.8775634765625
+32,0.5929999947547913,154.48672485351562,1309.30224609375,1787.0843505859375,1032.0728759765625
+33,0.6100000143051147,146.30194091796875,1325.1925048828125,1746.6524658203125,1034.538818359375
+34,0.6299999952316284,134.20530700683594,1312.0911865234375,1675.89697265625,1013.6607666015625
+35,0.656000018119812,115.96925354003906,1243.6280517578125,1540.450439453125,949.535888671875
+36,0.6675999760627747,113.64875793457031,1267.23583984375,1547.5260009765625,963.0390014648438
+37,0.6899999976158142,92.26598358154297,1114.307373046875,1435.327880859375,839.1528930664062
+38,0.7099999785423279,94.02474975585938,1197.5704345703125,1414.1011962890625,896.7203369140625
+39,0.7179999947547913,79.2312240600586,1045.2738037109375,1388.831298828125,779.8469848632812
+40,0.724399983882904,75.92426300048828,1022.4166870117188,1387.820556640625,761.2196044921875
+41,0.7400000095367432,80.37458038330078,1116.2152099609375,1312.010986328125,828.5402221679688
+42,0.7524999976158142,77.84077453613281,1116.2403564453125,1282.697998046875,826.0232543945312
+43,0.7574999928474426,75.56704711914062,1098.0296630859375,1258.4388427734375,811.54345703125
+44,0.762499988079071,45.746158599853516,696.9605102539062,1236.201416015625,512.89794921875
+45,0.7674999833106995,63.194129943847656,952.3276977539062,1218.007080078125,701.5109252929688
+46,0.7799999713897705,68.51187133789062,1054.7071533203125,1195.7696533203125,775.4505615234375
+47,0.800000011920929,62.87863540649414,1017.720458984375,1160.391845703125,745.0262451171875
+48,0.8159999847412109,51.39439010620117,871.8419799804688,1102.776611328125,635.7639770507812
+49,0.8237000107765198,47.48876953125,822.1177368164062,1073.4635009765625,598.5297241210938
+50,0.8314999938011169,50.44326400756836,882.2304077148438,1049.2044677734375,641.77587890625
+51,0.8399999737739563,50.50870895385742,897.87939453125,1033.03173828125,652.3304443359375
+52,0.8600000143051147,49.128204345703125,909.5646362304688,1009.4802856445312,658.7822265625
+53,0.8799999952316284,44.81117248535156,865.2081909179688,957.4243774414062,624.7343139648438
+54,0.9049999713897705,30.195653915405273,625.3776245117188,902.8414916992188,449.3675537109375
+55,0.9150000214576721,30.289173126220703,637.6793212890625,877.5715942382812,457.70654296875
+56,0.925000011920929,28.448226928710938,610.4505004882812,838.656005859375,437.61492919921875
+57,0.9300000071525574,19.592674255371094,432.5026550292969,839.262451171875,309.48626708984375
+58,0.9369999766349792,14.277915954589844,322.7142639160156,822.7865600585938,230.5836181640625
+59,0.9480000138282776,14.852182388305664,341.4164123535156,795.39404296875,243.69338989257812
+60,0.9649999737739563,25.331485748291016,583.9312744140625,776.59326171875,416.72314453125
+61,0.9800000190734863,27.016454696655273,635.7760009765625,775.2792358398438,453.1580810546875
+62,0.9934999942779541,28.801536560058594,689.1533203125,765.7777099609375,490.72039794921875
+63,1.0399999618530273,24.99117660522461,644.7619018554688,695.5275268554688,457.1557922363281
+64,1.0700000524520874,22.25499153137207,602.015625,647.6159057617188,425.76806640625
+65,1.100000023841858,16.722328186035156,479.10223388671875,612.7435302734375,337.8502197265625
+66,1.1200000047683716,5.06715202331543,155.4837188720703,592.224365234375,109.28324127197266
+67,1.1299999952316284,6.635468006134033,205.1408233642578,576.3549194335938,144.13522338867188
+68,1.1449999809265137,6.193835735321045,195.7523651123047,570.1890258789062,137.40078735351562
+69,1.1610000133514404,11.748560905456543,370.4906005859375,550.0742797851562,260.07733154296875
+70,1.1699999570846558,12.715644836425781,403.7008972167969,539.15771484375,283.30426025390625
+71,1.2000000476837158,12.658197402954102,416.2677001953125,507.014404296875,291.66998291015625
+72,1.2400000095367432,12.723058700561523,438.0545959472656,482.6542663574219,306.33795166015625
+73,1.2699999809265137,10.594919204711914,380.4599609375,447.4786376953125,265.6058349609375
+74,1.2899999618530273,10.970724105834961,402.2275085449219,444.7494812011719,280.5718078613281
+75,1.3200000524520874,8.522184371948242,326.9190979003906,421.2990417480469,227.64627075195312
+76,1.350000023841858,1.6226987838745117,67.14704132080078,395.6248779296875,46.62935256958008
+77,1.3949999809265137,0.12126490473747253,5.323389053344727,362.7740478515625,3.689373254776001
+78,1.4424999952316284,1.2762092351913452,58.43701171875,331.0351257324219,40.44478988647461
+79,1.462499976158142,2.304751396179199,106.74402618408203,320.92718505859375,73.8520736694336
+80,1.4769999980926514,2.212883234024048,104.03244018554688,310.6170654296875,71.94271850585938
+81,1.496999979019165,4.278558731079102,201.01097106933594,303.6426086425781,139.0102081298828
+82,1.5199999809265137,5.802607536315918,274.1861877441406,295.9605407714844,189.58140563964844
+83,1.5390000343322754,5.494048118591309,264.2200927734375,278.47381591796875,182.59288024902344
+84,1.5579999685287476,5.016831874847412,246.71206665039062,275.0371398925781,170.38055419921875
+85,1.5779999494552612,4.857605457305908,243.11094665527344,262.09893798828125,167.80760192871094
+86,1.5920000076293945,4.522202968597412,229.53021240234375,249.5651092529297,158.3694305419922
+87,1.6100000143051147,4.286829471588135,221.66986083984375,246.63380432128906,152.86550903320312
+88,1.6299999952316284,4.491518497467041,235.7653350830078,246.12840270996094,162.51797485351562
+89,1.6460000276565552,4.248502254486084,226.42987060546875,237.33450317382812,156.01766967773438
+90,1.6779999732971191,3.856564998626709,211.73875427246094,222.88014221191406,145.77871704101562
+91,1.7400000095367432,2.889831066131592,168.52455139160156,192.85955810546875,115.84679412841797
+92,1.7999999523162842,0.6722216010093689,42.90081787109375,172.94691467285156,29.427356719970703
+93,1.8600000143051147,0.03204738348722458,2.1758835315704346,146.0597686767578,1.4904770851135254
+94,1.9199999570846558,0.10780706256628036,7.664745330810547,137.16477966308594,5.2452569007873535
+95,1.9600000381469727,0.32299190759658813,23.525781631469727,124.32769775390625,16.09161949157715
+96,1.9850000143051147,1.270032286643982,91.32678985595703,125.13633728027344,62.483646392822266
+97,2.005000114440918,0.4182654023170471,31.37992286682129,114.21975708007812,21.451290130615234
+98,2.0350000858306885,1.2233593463897705,90.81049346923828,109.67118072509766,62.09091567993164
+99,2.065000057220459,0.9671005010604858,73.7654800415039,98.55244445800781,50.40989685058594
+100,2.0999999046325684,1.0698728561401367,83.00049591064453,93.39739227294922,56.70262145996094
+101,2.1480000019073486,0.966998815536499,77.45055389404297,83.2894515991211,52.87978744506836
+102,2.197999954223633,0.8404589295387268,69.74111938476562,75.4052505493164,47.58584213256836
+103,2.2699999809265137,0.7222429513931274,63.00203323364258,69.0372543334961,42.95062255859375
+104,2.359999895095825,0.5515680909156799,51.417625427246094,64.48867797851562,35.015262603759766
+105,2.450000047683716,0.17017853260040283,17.249914169311523,50.03431701660156,11.732279777526855
+106,2.5,0.05053719878196716,5.350593090057373,49.023521423339844,3.6368796825408936
+107,2.5999999046325684,3.1101667907762476e-09,3.513878539251891e-07,39.016658782958984,2.3863492515374674e-07
+108,2.700000047683716,2.7611552367440284e-12,3.305597739977628e-10,36.99506759643555,2.2432547486239685e-10
+109,2.799999952316284,2.0488089447212587e-08,2.5939837087207707e-06,32.34541702270508,1.7591577261555358e-06
+110,2.9000000953674316,0.007262714207172394,0.9694051146507263,28.4033203125,0.657025933265686
+111,3.0,0.02586463652551174,3.6192307472229004,25.067697525024414,2.4517266750335693
+112,3.0999999046325684,0.02316739596426487,3.4096922874450684,22.33855438232422,2.308582067489624
+113,3.200000047683716,0.032580118626356125,5.011773109436035,19.811569213867188,3.3918216228485107
+114,3.299999952316284,0.021878913044929504,3.5411267280578613,17.688899993896484,2.395390272140503
+115,3.4000000953674316,0.04572540521621704,7.63726282119751,15.86946964263916,5.1647539138793945
+116,3.5,0.0644257664680481,11.073113441467285,14.25220012664795,7.48640251159668
+117,3.5999999046325684,0.058482903987169266,10.483891487121582,12.837087631225586,7.085522174835205
+118,3.700000047683716,0.05231606587767601,9.779257774353027,11.624134063720703,6.607059955596924
+119,3.799999952316284,0.047291696071624756,9.203333854675293,10.512260437011719,6.216011047363281
+120,3.9000000953674316,0.04069847613573074,8.261123657226562,9.602545738220215,5.577882766723633
+121,4.0,0.039725467562675476,8.346587181091309,8.692831039428711,5.634192943572998
diff --git a/pvlib/spectrum/__init__.py b/pvlib/spectrum/__init__.py
new file mode 100644
index 0000000..36b3503
--- /dev/null
+++ b/pvlib/spectrum/__init__.py
@@ -0,0 +1,1 @@
+from pvlib.spectrum.spectrl2 import spectrl2  # noqa: F401
diff --git a/pvlib/spectrum/spectrl2.py b/pvlib/spectrum/spectrl2.py
new file mode 100644
index 0000000..15f006b
--- /dev/null
+++ b/pvlib/spectrum/spectrl2.py
@@ -0,0 +1,385 @@
+r"""
+The ``spectrl2`` module implements the Bird Simple Spectral Model.
+"""
+
+import pvlib
+from pvlib.tools import cosd
+import numpy as np
+import pandas as pd
+
+# SPECTRL2 extraterrestrial spectrum and atmospheric absorption coefficients
+_SPECTRL2_COEFFS = np.zeros(122, dtype=np.dtype([
+    ('wavelength', 'float64'),
+    ('spectral_irradiance_et', 'float64'),
+    ('water_vapor_absorption', 'float64'),
+    ('ozone_absorption', 'float64'),
+    ('mixed_absorption', 'float64'),
+]))
+_SPECTRL2_COEFFS['wavelength'] = [  # nm
+    300.0, 305.0, 310.0, 315.0, 320.0, 325.0, 330.0, 335.0, 340.0, 345.0,
+    350.0, 360.0, 370.0, 380.0, 390.0, 400.0, 410.0, 420.0, 430.0, 440.0,
+    450.0, 460.0, 470.0, 480.0, 490.0, 500.0, 510.0, 520.0, 530.0, 540.0,
+    550.0, 570.0, 593.0, 610.0, 630.0, 656.0, 667.6, 690.0, 710.0, 718.0,
+    724.4, 740.0, 752.5, 757.5, 762.5, 767.5, 780.0, 800.0, 816.0, 823.7,
+    831.5, 840.0, 860.0, 880.0, 905.0, 915.0, 925.0, 930.0, 937.0, 948.0,
+    965.0, 980.0, 993.5, 1040.0, 1070.0, 1100.0, 1120.0, 1130.0, 1145.0,
+    1161.0, 1170.0, 1200.0, 1240.0, 1270.0, 1290.0, 1320.0, 1350.0, 1395.0,
+    1442.5, 1462.5, 1477.0, 1497.0, 1520.0, 1539.0, 1558.0, 1578.0, 1592.0,
+    1610.0, 1630.0, 1646.0, 1678.0, 1740.0, 1800.0, 1860.0, 1920.0, 1960.0,
+    1985.0, 2005.0, 2035.0, 2065.0, 2100.0, 2148.0, 2198.0, 2270.0, 2360.0,
+    2450.0, 2500.0, 2600.0, 2700.0, 2800.0, 2900.0, 3000.0, 3100.0, 3200.0,
+    3300.0, 3400.0, 3500.0, 3600.0, 3700.0, 3800.0, 3900.0, 4000.0
+]
+_SPECTRL2_COEFFS['spectral_irradiance_et'] = [  # W/m^2/nm
+    0.5359, 0.5583, 0.622, 0.6927, 0.7151, 0.8329, 0.9619, 0.9319, 0.9006,
+    0.9113, 0.9755, 0.9759, 1.1199, 1.1038, 1.0338, 1.4791, 1.7013, 1.7404,
+    1.5872, 1.837, 2.005, 2.043, 1.987, 2.027, 1.896, 1.909, 1.927, 1.831,
+    1.891, 1.898, 1.892, 1.84, 1.768, 1.728, 1.658, 1.524, 1.531, 1.42,
+    1.399, 1.374, 1.373, 1.298, 1.269, 1.245, 1.223, 1.205, 1.183, 1.148,
+    1.091, 1.062, 1.038, 1.022, 0.9987, 0.9472, 0.8932, 0.8682, 0.8297,
+    0.8303, 0.814, 0.7869, 0.7683, 0.767, 0.7576, 0.6881, 0.6407, 0.6062,
+    0.5859, 0.5702, 0.5641, 0.5442, 0.5334, 0.5016, 0.4775, 0.4427, 0.44,
+    0.4168, 0.3914, 0.3589, 0.3275, 0.3175, 0.3073, 0.3004, 0.2928, 0.2755,
+    0.2721, 0.2593, 0.2469, 0.244, 0.2435, 0.2348, 0.2205, 0.1908, 0.1711,
+    0.1445, 0.1357, 0.123, 0.1238, 0.113, 0.1085, 0.0975, 0.0924, 0.0824,
+    0.0746, 0.0683, 0.0638, 0.0495, 0.0485, 0.0386, 0.0366, 0.032, 0.0281,
+    0.0248, 0.0221, 0.0196, 0.0175, 0.0157, 0.0141, 0.0127, 0.0115, 0.0104,
+    0.0095, 0.0086
+]
+_SPECTRL2_COEFFS['water_vapor_absorption'] = [
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.075, 0.0, 0.0, 0.0, 0.0, 0.016, 0.0125, 1.8, 2.5, 0.061,
+    0.0008, 0.0001, 1e-05, 1e-05, 0.0006, 0.036, 1.6, 2.5, 0.5, 0.155, 1e-05,
+    0.0026, 7.0, 5.0, 5.0, 27.0, 55.0, 45.0, 4.0, 1.48, 0.1, 1e-05, 0.001, 3.2,
+    115.0, 70.0, 75.0, 10.0, 5.0, 2.0, 0.002, 0.002, 0.1, 4.0, 200.0, 1000.0,
+    185.0, 80.0, 80.0, 12.0, 0.16, 0.002, 0.0005, 0.0001, 1e-05, 0.0001, 0.001,
+    0.01, 0.036, 1.1, 130.0, 1000.0, 500.0, 100.0, 4.0, 2.9, 1.0, 0.4, 0.22,
+    0.25, 0.33, 0.5, 4.0, 80.0, 310.0, 15000.0, 22000.0, 8000.0, 650.0, 240.0,
+    230.0, 100.0, 120.0, 19.5, 3.6, 3.1, 2.5, 1.4, 0.17, 0.0045
+]
+_SPECTRL2_COEFFS['ozone_absorption'] = [
+    10.0, 4.8, 2.7, 1.35, 0.8, 0.38, 0.16, 0.075, 0.04, 0.019, 0.007, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.003, 0.006, 0.009, 0.014, 0.021, 0.03,
+    0.04, 0.048, 0.063, 0.075, 0.085, 0.12, 0.119, 0.12, 0.09, 0.065, 0.051,
+    0.028, 0.018, 0.015, 0.012, 0.01, 0.008, 0.007, 0.006, 0.005, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+]
+_SPECTRL2_COEFFS['mixed_absorption'] = [
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.15, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0,
+    0.35, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.05, 0.3,
+    0.02, 0.0002, 0.00011, 1e-05, 0.05, 0.011, 0.005, 0.0006, 0.0, 0.005, 0.13,
+    0.04, 0.06, 0.13, 0.001, 0.0014, 0.0001, 1e-05, 1e-05, 0.0001, 0.001, 4.3,
+    0.2, 21.0, 0.13, 1.0, 0.08, 0.001, 0.00038, 0.001, 0.0005, 0.00015,
+    0.00014, 0.00066, 100.0, 150.0, 0.13, 0.0095, 0.001, 0.8, 1.9, 1.3, 0.075,
+    0.01, 0.00195, 0.004, 0.29, 0.025
+]
+
+
+def _spectrl2_transmittances(apparent_zenith, relative_airmass,
+                             surface_pressure, precipitable_water, ozone,
+                             optical_thickness, scattering_albedo, dayofyear):
+    """
+    Calculate transmittance factors from Section 2 of Bird and Riordan 1984.
+
+    Parameters
+    ----------
+    apparent_zenith, relative_airmass, surface_pressure, precipitable_water,
+    ozone, dayofyear: float or 1d np.array
+        One value per timestamp
+    optical_thickness, scattering_albedo: np.ndarray
+        Array with shape (122, N) where N is either 1 or len(apparent_zenith)
+
+    Returns
+    -------
+    earth_sun_distance_correction: float or 1d np.array
+        Same shape/type as apparent_zenith
+    rayleigh_transmittance, aerosol_transmittance, vapor_transmittance,
+    ozone_transmittance, mixed_transmittance, aerosol_scattering,
+    aerosol_absorption: np.ndarray
+        Array with shape (122, N) where N is len(apparent_zenith)
+    """
+    # add a dimension so that each ndarray is 2d with shape (122, 1)
+    wavelength = _SPECTRL2_COEFFS['wavelength'][:, np.newaxis]
+    vapor_coeff = _SPECTRL2_COEFFS['water_vapor_absorption'][:, np.newaxis]
+    ozone_coeff = _SPECTRL2_COEFFS['ozone_absorption'][:, np.newaxis]
+    mixed_coeff = _SPECTRL2_COEFFS['mixed_absorption'][:, np.newaxis]
+
+    # ET spectral irradiance correction for earth-sun distance seasonality.
+    # Note that we only want the distance correction coefficient, so set
+    # solar_constant=1:
+    earth_sun_distance_correction = \
+        pvlib.irradiance.get_extra_radiation(dayofyear, method='spencer',
+                                             solar_constant=1)  # Eq 2-2, 2-3
+    # Rayleigh scattering
+    # note: 101300 is used for consistentcy with reference; can't use
+    # atmosphere.get_absolute_airmass because it uses 101325
+    airmass = relative_airmass * surface_pressure / 101300
+    wavelength_um = wavelength / 1000
+    rayleigh_transmittance = np.exp(
+        # Note: the report uses 1.335 but spectrl2_2.c uses 1.3366
+        # -airmass / (wavelength_um**4 * (115.6406 - 1.335 / wavelength_um**2))
+        -airmass / (wavelength_um**4 * (115.6406 - 1.3366 / wavelength_um**2))
+    )  # Eq 2-4
+
+    # Aerosol scattering and absorption, Eq 2-6
+    aerosol_transmittance = np.exp(-optical_thickness * relative_airmass)
+
+    # Water vapor absorption, Eq 2-8
+    aWM = vapor_coeff * precipitable_water * relative_airmass
+    vapor_transmittance = np.exp(-0.2385 * aWM / (1 + 20.07 * aWM)**0.45)
+
+    # Ozone absorption
+    ozone_max_height = 22
+    h0_norm = ozone_max_height / 6370
+    ozone_mass_numerator = (1 + h0_norm)
+    ozone_mass_denominator = np.sqrt(cosd(apparent_zenith)**2 + 2 * h0_norm)
+    ozone_mass = ozone_mass_numerator / ozone_mass_denominator  # Eq 2-10
+    ozone_transmittance = np.exp(-ozone_coeff * ozone * ozone_mass)  # Eq 2-9
+
+    # Mixed gas absorption, Eq 2-11
+    aM = mixed_coeff * airmass
+    # Note: the report uses 118.93, but spectrl2_2.c uses 118.3
+    # mixed_transmittance = np.exp(-1.41 * aM / (1 + 118.93 * aM)**0.45)
+    mixed_transmittance = np.exp(-1.41 * aM / (1 + 118.3 * aM)**0.45)
+
+    # split out aerosol components for diffuse irradiance calcs
+    aerosol_scattering = np.exp(
+        -scattering_albedo * optical_thickness * relative_airmass
+    )  # Eq 3-9
+
+    aerosol_absorption = np.exp(
+        -(1 - scattering_albedo) * optical_thickness * relative_airmass
+    )  # Eq 3-10
+
+    return (
+        earth_sun_distance_correction,
+        rayleigh_transmittance,
+        aerosol_transmittance,
+        vapor_transmittance,
+        ozone_transmittance,
+        mixed_transmittance,
+        aerosol_scattering,
+        aerosol_absorption,
+    )
+
+
+def spectrl2(apparent_zenith, aoi, surface_tilt, ground_albedo,
+             surface_pressure, relative_airmass, precipitable_water, ozone,
+             aerosol_turbidity_500nm, dayofyear=None,
+             scattering_albedo_400nm=0.945, alpha=1.14,
+             wavelength_variation_factor=0.095, aerosol_asymmetry_factor=0.65):
+    """
+    Estimate spectral irradiance using the Bird Simple Spectral Model
+    (SPECTRL2).
+
+    The Bird Simple Spectral Model [1]_ produces terrestrial spectra between
+    300 and 4000 nm with a resolution of approximately 10 nm. Direct and
+    diffuse spectral irradiance are modeled for horizontal and tilted surfaces
+    under cloudless skies. SPECTRL2 models radiative transmission, absorption,
+    and scattering due to atmospheric aerosol, water vapor, and ozone content.
+
+    Parameters
+    ----------
+    apparent_zenith : numeric
+        Solar zenith angle [degrees]
+    aoi : numeric
+        Angle of incidence of the solar vector on the panel [degrees]
+    surface_tilt : numeric
+        Panel tilt from horizontal [degrees]
+    ground_albedo : numeric
+        Albedo [0-1] of the ground surface. Can be provided as a scalar value
+        if albedo is not spectrally-dependent, or as a 122xN matrix where
+        the first dimension spans the wavelength range and the second spans
+        the number of simulations. [unitless]
+    surface_pressure : numeric
+        Surface pressure [Pa]
+    relative_airmass : numeric
+        Relative airmass. The airmass model used in [1]_ is the `'kasten1966'`
+        model, while a later implementation by NREL uses the
+        `'kastenyoung1989'` model. [unitless]
+    precipitable_water : numeric
+        Atmospheric water vapor content [cm]
+    ozone : numeric
+        Atmospheric ozone content [atm-cm]
+    aerosol_turbidity_500nm : numeric
+        Aerosol turbidity at 500 nm [unitless]
+    dayofyear : numeric, optional
+        The day of year [1-365].  Must be provided if ``apparent_zenith`` is
+        not a pandas Series.
+    scattering_albedo_400nm : numeric, default 0.945
+        Aerosol single scattering albedo at 400nm. The default value of 0.945
+        is suggested in [1]_ for a rural aerosol model. [unitless]
+    alpha : numeric, default 1.14
+        Angstrom turbidity exponent. The default value of 1.14 is suggested
+        in [1]_ for a rural aerosol model. [unitless]
+    wavelength_variation_factor : numeric, default 0.095
+        Wavelength variation factor [unitless]
+    aerosol_asymmetry_factor : numeric, default 0.65
+        Aerosol asymmetry factor (mean cosine of scattering angle) [unitless]
+
+    Returns
+    -------
+    spectra : dict
+        A dict of arrays.  With the exception of `wavelength`, which has length
+        122, each array has shape (122, N) where N is the length of the
+        input ``apparent_zenith``.  All values are spectral irradiance
+        with units W/m^2/nm except for `wavelength`, which is in nanometers.
+
+            * wavelength
+            * dni_extra
+            * dhi
+            * dni
+            * poa_sky_diffuse
+            * poa_ground_diffuse
+            * poa_direct
+            * poa_global
+
+    Notes
+    -----
+    NREL's C implementation ``spectrl2_2.c`` [2]_ of the model differs in
+    several ways from the original report [1]_.  The report itself also has
+    a few differences between the in-text equations and the code appendix.
+    The list of known differences is shown below.  Note that this
+    implementation follows ``spectrl2_2.c``.
+
+    =================== ========== ========== ===============
+    Equation            Report     Appendix   spectrl2_2.c
+    =================== ========== ========== ===============
+    2-4                 1.335      1.335      1.3366
+    2-11                118.93     118.93     118.3
+    3-8                 To'        Tu'        Tu'
+    3-5, 3-6, 3-7, 3-1  double Cs  single Cs  single Cs
+    2-5                 kasten1966 kasten1966 kastenyoung1989
+    =================== ========== ========== ===============
+
+    References
+    ----------
+    .. [1] Bird, R, and Riordan, C., 1984, "Simple solar spectral model for
+       direct and diffuse irradiance on horizontal and tilted planes at the
+       earth's surface for cloudless atmospheres", NREL Technical Report
+       TR-215-2436 doi:10.2172/5986936.
+    .. [2] Bird Simple Spectral Model: spectrl2_2.c.
+       https://www.nrel.gov/grid/solar-resource/spectral.html
+    """
+    # values need to be np arrays for broadcasting, so unwrap Series if needed:
+    is_pandas = isinstance(apparent_zenith, pd.Series)
+    if is_pandas:
+        original_index = apparent_zenith.index
+        (apparent_zenith, aoi, surface_tilt, ground_albedo, surface_pressure,
+         relative_airmass, precipitable_water, ozone, aerosol_turbidity_500nm,
+         scattering_albedo_400nm, alpha, wavelength_variation_factor,
+         aerosol_asymmetry_factor) = \
+            tuple(map(np.asanyarray, [
+                apparent_zenith, aoi, surface_tilt, ground_albedo,
+                surface_pressure, relative_airmass, precipitable_water, ozone,
+                aerosol_turbidity_500nm, scattering_albedo_400nm, alpha,
+                wavelength_variation_factor, aerosol_asymmetry_factor]))
+
+        dayofyear = original_index.dayofyear.values
+
+    if not is_pandas and dayofyear is None:
+        raise ValueError('dayofyear must be specified if not using pandas '
+                         'Series inputs')
+
+    # add a dimension so that each ndarray is 2d with shape (122, 1)
+    wavelength = _SPECTRL2_COEFFS['wavelength'][:, np.newaxis]
+    spectrum_et = _SPECTRL2_COEFFS['spectral_irradiance_et'][:, np.newaxis]
+
+    optical_thickness = \
+        pvlib.atmosphere.angstrom_aod_at_lambda(aod0=aerosol_turbidity_500nm,
+                                                lambda0=500, alpha=alpha,
+                                                lambda1=wavelength)  # Eq 2-7
+
+    # Eq 3-16
+    scattering_albedo = scattering_albedo_400nm * \
+        np.exp(-wavelength_variation_factor * np.log(wavelength / 400)**2)
+
+    spectrl2 = _spectrl2_transmittances(apparent_zenith, relative_airmass,
+                                        surface_pressure, precipitable_water,
+                                        ozone, optical_thickness,
+                                        scattering_albedo, dayofyear)
+    D, Tr, Ta, Tw, To, Tu, Tas, Taa = spectrl2
+
+    spectrum_et_adj = spectrum_et * D
+    # spectrum of direct irradiance, Eq 2-1
+    Id = spectrum_et_adj * Tr * Ta * Tw * To * Tu
+
+    cosZ = cosd(apparent_zenith)
+    # Eq 3-17
+    Cs = np.where(wavelength <= 450, ((wavelength + 550)/1000)**1.8, 1.0)
+    ALG = np.log(1 - aerosol_asymmetry_factor)  # Eq 3-14
+    BFS = ALG * (0.0783 + ALG * (-0.3824 - ALG * 0.5874))  # Eq 3-13
+    AFS = ALG * (1.459 + ALG * (0.1595 + ALG * 0.4129))  # Eq 3-12
+    Fs = 1 - 0.5 * np.exp((AFS + BFS * cosZ) * cosZ)  # Eq 3-11
+    Fsp = 1 - 0.5 * np.exp((AFS + BFS / 1.8) / 1.8)  # Eq 3.15
+
+    # evaluate the "primed terms" -- transmittances evaluated at airmass=1.8
+    primes = _spectrl2_transmittances(apparent_zenith, 1.8,
+                                      surface_pressure, precipitable_water,
+                                      ozone, optical_thickness,
+                                      scattering_albedo, dayofyear)
+    _, Trp, Tap, Twp, Top, Tup, Tasp, Taap = primes
+
+    # Note: not sure what the correct form of this equation is.
+    # The first coefficient is To' in Eq 3-8 but Tu' in the code appendix.
+    # spectrl2_2.c uses Tu'.
+    sky_reflectivity = (
+        # Top * Twp * Taap * (0.5 * (1-Trp) + (1-Fsp) * Trp * (1-Tasp))
+        Tup * Twp * Taap * (0.5 * (1-Trp) + (1-Fsp) * Trp * (1-Tasp))
+    )  # Eq 3-8
+
+    # a common factor for 3-5 and 3-6
+    common_factor = spectrum_et_adj * cosZ * To * Tu * Tw * Taa
+    # Note: spectrl2_2.c differs from the report in how the Cs value is used.
+    # The two commented out lines match the report, while the following match
+    # spectrl2_2.c. With regard to Cs, the equations in the report and
+    # spectrl12_2.c are algebraically equivalent.
+    # Ir = common_factor * (1 - Tr**0.95) * 0.5 * Cs  # Eq 3-5
+    # Ia = common_factor * Tr**1.5 * (1 - Tas) * Fs * Cs  # Eq 3-6
+    Ir = common_factor * (1 - Tr**0.95) * 0.5  # Eq 3-5
+    Ia = common_factor * Tr**1.5 * (1 - Tas) * Fs  # Eq 3-6
+
+    rs = sky_reflectivity
+    rg = ground_albedo
+    Ig = (Id * cosZ + Ir + Ia) * rs * rg / (1 - rs * rg)  # Eq 3-7
+
+    # total scattered irradiance
+    # Note: see discussion about Cs above.
+    # Is = Ir + Ia + Ig  # Eq 3-1
+    Is = (Ir + Ia + Ig) * Cs  # Eq 3-1
+
+    # calculate spectral irradiance on a tilted surface, Eq 3-18
+    Ibeam = Id * cosd(aoi)
+
+    # don't need surface_azimuth if we provide projection_ratio
+    projection_ratio = cosd(aoi) / cosZ
+    Isky = pvlib.irradiance.haydavies(surface_tilt=surface_tilt,
+                                      surface_azimuth=None,
+                                      dhi=Is,
+                                      dni=Id,
+                                      dni_extra=spectrum_et_adj,
+                                      projection_ratio=projection_ratio)
+
+    ghi = Id * cosZ + Is
+    Iground = pvlib.irradiance.get_ground_diffuse(surface_tilt, ghi, albedo=rg)
+
+    Itilt = Ibeam + Isky + Iground
+    wavelength_1d = wavelength.reshape(-1)  # only needs 1 dimension
+    return {
+        'wavelength': wavelength_1d,
+        'dni_extra': spectrum_et_adj,
+        'dhi': Is,
+        'dni': Id,
+        'poa_sky_diffuse': Isky,
+        'poa_ground_diffuse': Iground,
+        'poa_direct': Ibeam,
+        'poa_global': Itilt,
+    }
diff --git a/setup.py b/setup.py
index b89a3c9..e2bacea 100755
--- a/setup.py
+++ b/setup.py
@@ -73,6 +73,12 @@ setuptools_kwargs = {
     'python_requires': '>=3.6'
 }
 
+PROJECT_URLS = {
+    "Bug Tracker": "https://github.com/pvlib/pvlib-python/issues",
+    "Documentation": "https://pvlib-python.readthedocs.io/",
+    "Source Code": "https://github.com/pvlib/pvlib-python",
+}
+
 # set up pvlib packages to be installed and extensions to be compiled
 PACKAGES = ['pvlib']
 
@@ -109,5 +115,6 @@ setup(name=DISTNAME,
       maintainer_email=MAINTAINER_EMAIL,
       license=LICENSE,
       url=URL,
+      project_urls=PROJECT_URLS,
       classifiers=CLASSIFIERS,
       **setuptools_kwargs)


Test patch:
diff --git a/pvlib/tests/iotools/test_pvgis.py b/pvlib/tests/iotools/test_pvgis.py
index 1957639..0dce093 100644
--- a/pvlib/tests/iotools/test_pvgis.py
+++ b/pvlib/tests/iotools/test_pvgis.py
@@ -159,7 +159,8 @@ def _compare_pvgis_tmy_csv(expected, month_year_expected, inputs_expected,
     for meta_value in meta:
         if not meta_value:
             continue
-        if meta_value == 'PVGIS (c) European Communities, 2001-2020':
+        # don't check end year because it changes every year
+        if meta_value[:-4] == 'PVGIS (c) European Communities, 2001-':
             continue
         assert meta_value in csv_meta
 
diff --git a/pvlib/tests/test_clearsky.py b/pvlib/tests/test_clearsky.py
index 02613fc..fd886f2 100644
--- a/pvlib/tests/test_clearsky.py
+++ b/pvlib/tests/test_clearsky.py
@@ -4,6 +4,7 @@ import numpy as np
 from numpy import nan
 import pandas as pd
 import pytz
+from scipy.linalg import hankel
 
 import pytest
 from numpy.testing import assert_allclose
@@ -549,7 +550,15 @@ def detect_clearsky_data():
 def test_detect_clearsky(detect_clearsky_data):
     expected, cs = detect_clearsky_data
     clear_samples = clearsky.detect_clearsky(
-        expected['GHI'], cs['ghi'], cs.index, 10)
+        expected['GHI'], cs['ghi'], times=cs.index, window_length=10)
+    assert_series_equal(expected['Clear or not'], clear_samples,
+                        check_dtype=False, check_names=False)
+
+
+def test_detect_clearsky_defaults(detect_clearsky_data):
+    expected, cs = detect_clearsky_data
+    clear_samples = clearsky.detect_clearsky(
+        expected['GHI'], cs['ghi'])
     assert_series_equal(expected['Clear or not'], clear_samples,
                         check_dtype=False, check_names=False)
 
@@ -557,7 +566,8 @@ def test_detect_clearsky(detect_clearsky_data):
 def test_detect_clearsky_components(detect_clearsky_data):
     expected, cs = detect_clearsky_data
     clear_samples, components, alpha = clearsky.detect_clearsky(
-        expected['GHI'], cs['ghi'], cs.index, 10, return_components=True)
+        expected['GHI'], cs['ghi'], times=cs.index, window_length=10,
+        return_components=True)
     assert_series_equal(expected['Clear or not'], clear_samples,
                         check_dtype=False, check_names=False)
     assert isinstance(components, OrderedDict)
@@ -569,11 +579,11 @@ def test_detect_clearsky_iterations(detect_clearsky_data):
     alpha = 1.0448
     with pytest.warns(RuntimeWarning):
         clear_samples = clearsky.detect_clearsky(
-            expected['GHI'], cs['ghi']*alpha, cs.index, 10, max_iterations=1)
-    assert (clear_samples[:'2012-04-01 10:41:00'] == True).all()
-    assert (clear_samples['2012-04-01 10:42:00':] == False).all()
+            expected['GHI'], cs['ghi']*alpha, max_iterations=1)
+    assert clear_samples[:'2012-04-01 10:41:00'].all()
+    assert not clear_samples['2012-04-01 10:42:00':].any()  # expected False
     clear_samples = clearsky.detect_clearsky(
-            expected['GHI'], cs['ghi']*alpha, cs.index, 10, max_iterations=20)
+        expected['GHI'], cs['ghi']*alpha, max_iterations=20)
     assert_series_equal(expected['Clear or not'], clear_samples,
                         check_dtype=False, check_names=False)
 
@@ -581,7 +591,7 @@ def test_detect_clearsky_iterations(detect_clearsky_data):
 def test_detect_clearsky_kwargs(detect_clearsky_data):
     expected, cs = detect_clearsky_data
     clear_samples = clearsky.detect_clearsky(
-        expected['GHI'], cs['ghi'], cs.index, 10,
+        expected['GHI'], cs['ghi'], times=cs.index, window_length=10,
         mean_diff=1000, max_diff=1000, lower_line_length=-1000,
         upper_line_length=1000, var_diff=10, slope_dev=1000)
     assert clear_samples.all()
@@ -590,7 +600,7 @@ def test_detect_clearsky_kwargs(detect_clearsky_data):
 def test_detect_clearsky_window(detect_clearsky_data):
     expected, cs = detect_clearsky_data
     clear_samples = clearsky.detect_clearsky(
-        expected['GHI'], cs['ghi'], cs.index, 3)
+        expected['GHI'], cs['ghi'], window_length=3)
     expected = expected['Clear or not'].copy()
     expected.iloc[-3:] = True
     assert_series_equal(expected, clear_samples,
@@ -600,7 +610,8 @@ def test_detect_clearsky_window(detect_clearsky_data):
 def test_detect_clearsky_arrays(detect_clearsky_data):
     expected, cs = detect_clearsky_data
     clear_samples = clearsky.detect_clearsky(
-        expected['GHI'].values, cs['ghi'].values, cs.index, 10)
+        expected['GHI'].values, cs['ghi'].values, times=cs.index,
+        window_length=10)
     assert isinstance(clear_samples, np.ndarray)
     assert (clear_samples == expected['Clear or not'].values).all()
 
@@ -611,8 +622,72 @@ def test_detect_clearsky_irregular_times(detect_clearsky_data):
     times[0] += 10**9
     times = pd.DatetimeIndex(times)
     with pytest.raises(NotImplementedError):
-        clear_samples = clearsky.detect_clearsky(
-            expected['GHI'].values, cs['ghi'].values, times, 10)
+        clearsky.detect_clearsky(expected['GHI'].values, cs['ghi'].values,
+                                 times, 10)
+
+
+def test_detect_clearsky_missing_index(detect_clearsky_data):
+    expected, cs = detect_clearsky_data
+    with pytest.raises(ValueError):
+        clearsky.detect_clearsky(expected['GHI'].values, cs['ghi'].values)
+
+
+@pytest.fixture
+def detect_clearsky_helper_data():
+    samples_per_window = 3
+    sample_interval = 1
+    x = pd.Series(np.arange(0, 7)**2.)
+    # line length between adjacent points
+    sqt = pd.Series(np.sqrt(np.array([np.nan, 2., 10., 26., 50., 82, 122.])))
+    H = hankel(np.arange(samples_per_window),
+               np.arange(samples_per_window-1, len(sqt)))
+    return x, samples_per_window, sample_interval, H
+
+
+def test__line_length_windowed(detect_clearsky_helper_data):
+    x, samples_per_window, sample_interval, H = detect_clearsky_helper_data
+    # sqt is hand-calculated assuming window=3
+    # line length between adjacent points
+    sqt = pd.Series(np.sqrt(np.array([np.nan, 2., 10., 26., 50., 82, 122.])))
+    expected = {}
+    expected['line_length'] = sqt + sqt.shift(-1)
+    result = clearsky._line_length_windowed(
+        x, H, samples_per_window, sample_interval)
+    assert_series_equal(result, expected['line_length'])
+
+
+def test__max_diff_windowed(detect_clearsky_helper_data):
+    x, samples_per_window, sample_interval, H = detect_clearsky_helper_data
+    expected = {}
+    expected['max_diff'] = pd.Series(
+        data=[np.nan, 3., 5., 7., 9., 11., np.nan], index=x.index)
+    result = clearsky._max_diff_windowed(x, H, samples_per_window)
+    assert_series_equal(result, expected['max_diff'])
+
+
+def test__calc_stats(detect_clearsky_helper_data):
+    x, samples_per_window, sample_interval, H = detect_clearsky_helper_data
+    # stats are hand-computed assuming window = 3, sample_interval = 1,
+    # and right-aligned labels
+    mean_x = pd.Series(np.array([np.nan, np.nan, 5, 14, 29, 50, 77]) / 3.)
+    max_x = pd.Series(np.array([np.nan, np.nan, 4, 9, 16, 25, 36]))
+    diff_std = np.array([np.nan, np.nan, np.sqrt(2), np.sqrt(2), np.sqrt(2),
+                         np.sqrt(2), np.sqrt(2)])
+    slope_nstd = diff_std / mean_x
+    slope = x.diff().shift(-1)
+    expected = {}
+    expected['mean'] = mean_x.shift(-1)  # shift to align to center
+    expected['max'] = max_x.shift(-1)
+    # slope between adjacent points
+    expected['slope'] = slope
+    expected['slope_nstd'] = slope_nstd.shift(-1)
+    result = clearsky._calc_stats(
+        x, samples_per_window, sample_interval, H)
+    res_mean, res_max, res_slope_nstd, res_slope = result
+    assert_series_equal(res_mean, expected['mean'])
+    assert_series_equal(res_max, expected['max'])
+    assert_series_equal(res_slope_nstd, expected['slope_nstd'])
+    assert_series_equal(res_slope, expected['slope'])
 
 
 def test_bird():
diff --git a/pvlib/tests/test_spectrum.py b/pvlib/tests/test_spectrum.py
new file mode 100644
index 0000000..a28ceca
--- /dev/null
+++ b/pvlib/tests/test_spectrum.py
@@ -0,0 +1,94 @@
+import pytest
+from numpy.testing import assert_allclose
+import pandas as pd
+import numpy as np
+from pvlib import spectrum
+from conftest import DATA_DIR
+
+SPECTRL2_TEST_DATA = DATA_DIR / 'spectrl2_example_spectra.csv'
+
+
+@pytest.fixture
+def spectrl2_data():
+    # reference spectra generated with solar_utils==0.3
+    """
+    expected = solar_utils.spectrl2(
+        units=1,
+        location=[40, -80, -5],
+        datetime=[2020, 3, 15, 10, 45, 59],
+        weather=[1013, 15],
+        orientation=[0, 180],
+        atmospheric_conditions=[1.14, 0.65, 0.344, 0.1, 1.42],
+        albedo=[0.3, 0.7, 0.8, 1.3, 2.5, 4.0] + [0.2]*6,
+    )
+    """
+    kwargs = {
+        'surface_tilt': 0,
+        'relative_airmass': 1.4899535986910446,
+        'apparent_zenith': 47.912086486816406,
+        'aoi': 47.91208648681641,
+        'ground_albedo': 0.2,
+        'surface_pressure': 101300,
+        'ozone': 0.344,
+        'precipitable_water': 1.42,
+        'aerosol_turbidity_500nm': 0.1,
+        'dayofyear': 75
+    }
+    df = pd.read_csv(SPECTRL2_TEST_DATA)
+    # convert um to nm
+    df['wavelength'] *= 1000
+    df[['specdif', 'specdir', 'specetr', 'specglo']] /= 1000
+    return kwargs, df
+
+
+def test_spectrl2(spectrl2_data):
+    # compare against output from solar_utils wrapper around NREL spectrl2_2.c
+    kwargs, expected = spectrl2_data
+    actual = spectrum.spectrl2(**kwargs)
+    assert_allclose(expected['wavelength'].values, actual['wavelength'])
+    assert_allclose(expected['specdif'].values, actual['dhi'].ravel(),
+                    atol=7e-5)
+    assert_allclose(expected['specdir'].values, actual['dni'].ravel(),
+                    atol=1.5e-4)
+    assert_allclose(expected['specetr'], actual['dni_extra'].ravel(),
+                    atol=2e-4)
+    assert_allclose(expected['specglo'], actual['poa_global'].ravel(),
+                    atol=1e-4)
+
+
+def test_spectrl2_array(spectrl2_data):
+    # test that supplying arrays instead of scalars works
+    kwargs, expected = spectrl2_data
+    kwargs = {k: np.array([v, v, v]) for k, v in kwargs.items()}
+    actual = spectrum.spectrl2(**kwargs)
+
+    assert actual['wavelength'].shape == (122,)
+
+    keys = ['dni_extra', 'dhi', 'dni', 'poa_sky_diffuse', 'poa_ground_diffuse',
+            'poa_direct', 'poa_global']
+    for key in keys:
+        assert actual[key].shape == (122, 3)
+
+
+def test_spectrl2_series(spectrl2_data):
+    # test that supplying Series instead of scalars works
+    kwargs, expected = spectrl2_data
+    kwargs.pop('dayofyear')
+    index = pd.to_datetime(['2020-03-15 10:45:59']*3)
+    kwargs = {k: pd.Series([v, v, v], index=index) for k, v in kwargs.items()}
+    actual = spectrum.spectrl2(**kwargs)
+
+    assert actual['wavelength'].shape == (122,)
+
+    keys = ['dni_extra', 'dhi', 'dni', 'poa_sky_diffuse', 'poa_ground_diffuse',
+            'poa_direct', 'poa_global']
+    for key in keys:
+        assert actual[key].shape == (122, 3)
+
+
+def test_dayofyear_missing(spectrl2_data):
+    # test that not specifying dayofyear with non-pandas inputs raises error
+    kwargs, expected = spectrl2_data
+    kwargs.pop('dayofyear')
+    with pytest.raises(ValueError, match='dayofyear must be specified'):
+        _ = spectrum.spectrl2(**kwargs)