From c6116eb52bb2fa2b23560947f45f340a1f91cb35 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <mark.mikofski@dnvgl.com>
Date: Thu, 25 Feb 2021 20:34:28 -0800
Subject: [PATCH 01/30] add Boland DF estimation

---
 pvlib/irradiance.py            | 95 ++++++++++++++++++++++++++++++++++
 pvlib/tests/test_irradiance.py | 16 ++++++
 2 files changed, 111 insertions(+)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 3ec6b213f9..0bec078c2d 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -2186,6 +2186,101 @@ def erbs(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
     return data
 
 
+def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
+    r"""
+    Estimate DNI and DHI from GHI using the Boland clearness index model.
+
+    The Boland model [1]_, [2]_ estimates the diffuse fraction DF from global
+    horizontal irradiance through an empirical relationship between DF
+    and the ratio of GHI to extraterrestrial irradiance, Kt. The
+    function uses the diffuse fraction to compute DHI as
+
+    .. math::
+
+        DHI = DF \times GHI
+
+    DNI is then estimated as
+
+    .. math::
+
+        DNI = (GHI - DHI)/\cos(Z)
+
+    where Z is the zenith angle.
+
+    Parameters
+    ----------
+    ghi: numeric
+        Global horizontal irradiance in W/m^2.
+    zenith: numeric
+        True (not refraction-corrected) zenith angles in decimal degrees.
+    datetime_or_doy : int, float, array, pd.DatetimeIndex
+        Day of year or array of days of year e.g.
+        pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
+    min_cos_zenith : numeric, default 0.065
+        Minimum value of cos(zenith) to allow when calculating global
+        clearness index `kt`. Equivalent to zenith = 86.273 degrees.
+    max_zenith : numeric, default 87
+        Maximum value of zenith to allow in DNI calculation. DNI will be
+        set to 0 for times with zenith values greater than `max_zenith`.
+
+    Returns
+    -------
+    data : OrderedDict or DataFrame
+        Contains the following keys/columns:
+
+            * ``dni``: the modeled direct normal irradiance in W/m^2.
+            * ``dhi``: the modeled diffuse horizontal irradiance in
+              W/m^2.
+            * ``kt``: Ratio of global to extraterrestrial irradiance
+              on a horizontal plane.
+
+    References
+    ----------
+    .. [1] John Boland, Lynne Scott, and Mark Luther, Modelling the diffuse
+       fraction of global solar radiation on a horizontal surface,
+       Environmetrics 12(2), pp 103-116, 2001,
+       :doi:`10.1002/1099-095X(200103)12:2%3C103::AID-ENV447%3E3.0.CO;2-2`
+    .. [2] J. Boland, B. Ridley (2008) Models of Diffuse Solar Fraction. In:
+       Badescu V. (eds) Modeling Solar Radiation at the Earth’s Surface.
+       Springer, Berlin, Heidelberg. :doi:`10.1007/978-3-540-77455-6_8`
+
+    See also
+    --------
+    dirint
+    disc
+    erbs
+    """
+
+    dni_extra = get_extra_radiation(datetime_or_doy)
+
+    kt = clearness_index(ghi, zenith, dni_extra, min_cos_zenith=min_cos_zenith,
+                         max_clearness_index=1)
+
+    # Boland equation
+    df = 1.0 / (1.0 + np.exp(-5.0 + 8.6 * kt))
+    # NOTE: [1] has different coefficients, for different time intervals
+    # 15-min: df = 1 / (1 + exp(8.645 * (kt - 0.613)))
+    # 1-hour: df = 1 / (1 + exp(7.997 * (kt - 0.586)))
+
+    dhi = df * ghi
+
+    dni = (ghi - dhi) / tools.cosd(zenith)
+    bad_values = (zenith > max_zenith) | (ghi < 0) | (dni < 0)
+    dni = np.where(bad_values, 0, dni)
+    # ensure that closure relationship remains valid
+    dhi = np.where(bad_values, ghi, dhi)
+
+    data = OrderedDict()
+    data['dni'] = dni
+    data['dhi'] = dhi
+    data['kt'] = kt
+
+    if isinstance(datetime_or_doy, pd.DatetimeIndex):
+        data = pd.DataFrame(data, index=datetime_or_doy)
+
+    return data
+
+
 def campbell_norman(zenith, transmittance, pressure=101325.0,
                     dni_extra=1367.0):
     '''
diff --git a/pvlib/tests/test_irradiance.py b/pvlib/tests/test_irradiance.py
index ba5821f750..e005de83b8 100644
--- a/pvlib/tests/test_irradiance.py
+++ b/pvlib/tests/test_irradiance.py
@@ -636,6 +636,22 @@ def test_erbs():
     assert_frame_equal(np.round(out, 0), np.round(expected, 0))
 
 
+def test_boland():
+    index = pd.DatetimeIndex(['20190101']*3 + ['20190620'])
+    ghi = pd.Series([0, 50, 1000, 1000], index=index)
+    zenith = pd.Series([120, 85, 10, 10], index=index)
+    expected = pd.DataFrame(np.array(
+        [[0.0,        0.0,        0.0],
+         [103.735879, 40.958822,  0.405724],
+         [776.006568, 235.782716, 0.718133],
+         [845.794317, 167.055199, 0.768214]]),
+        columns=['dni', 'dhi', 'kt'], index=index)
+
+    out = irradiance.boland(ghi, zenith, index)
+
+    assert_frame_equal(np.round(out, 0), np.round(expected, 0))
+
+
 def test_erbs_min_cos_zenith_max_zenith():
     # map out behavior under difficult conditions with various
     # limiting kwargs settings

From fc9fad21c48fc1c8064337d85267180e2d2b856c Mon Sep 17 00:00:00 2001
From: Mark Mikofski <mark.mikofski@dnvgl.com>
Date: Thu, 25 Feb 2021 20:42:43 -0800
Subject: [PATCH 02/30] update what's new, add to docs

---
 docs/sphinx/source/api.rst             | 1 +
 docs/sphinx/source/whatsnew/v0.9.0.rst | 4 ++++
 2 files changed, 5 insertions(+)

diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 8805d199a4..9f026f4e05 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -182,6 +182,7 @@ DNI estimation models
    irradiance.dirint
    irradiance.dirindex
    irradiance.erbs
+   irradiance.boland
    irradiance.campbell_norman
    irradiance.gti_dirint
 
diff --git a/docs/sphinx/source/whatsnew/v0.9.0.rst b/docs/sphinx/source/whatsnew/v0.9.0.rst
index 81e7a0c60b..6e037d2e22 100644
--- a/docs/sphinx/source/whatsnew/v0.9.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.0.rst
@@ -96,6 +96,10 @@ Enhancements
 * :py:meth:`~pvlib.pvsystem.PVSystem.get_ac` is added to calculate AC power
   from DC power. Use parameter ``model`` to specify which inverter model to use.
   (:pull:`1147`, :issue:`998`, :pull:`1150`)
+* :py:meth:`~pvlib.irradiance.boland` is another diffuse fraction, DF,
+  estimation method similar to Erbs but uses a single logistic exponential
+  correlation between DF and clearness index, kt, that is continuously
+  differentiable and bounded between zero and one. (:pull:`1179`)
 
 Bug fixes
 ~~~~~~~~~

From fc358b68a27cbf1e310d8deb060a0b6721eaebe6 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <mark.mikofski@dnvgl.com>
Date: Thu, 25 Feb 2021 21:58:29 -0800
Subject: [PATCH 03/30] add example

---
 docs/examples/plot_diffuse_fraction.py | 74 ++++++++++++++++++++++++++
 pvlib/irradiance.py                    | 17 ++----
 2 files changed, 79 insertions(+), 12 deletions(-)
 create mode 100644 docs/examples/plot_diffuse_fraction.py

diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/plot_diffuse_fraction.py
new file mode 100644
index 0000000000..80a6612177
--- /dev/null
+++ b/docs/examples/plot_diffuse_fraction.py
@@ -0,0 +1,74 @@
+"""
+Diffuse Fraction Estimation
+===========================
+
+Comparison of diffuse fraction estimation methods used to derive direct and
+diffuse components from measured global horizontal irradiance.
+"""
+# %%
+# PV systems are often tilted to optimize performance. Determining the total
+# irradiance incident on the plane of the array requires transposing the
+# diffuse component, because the entire isotropic sky dome is not visible to
+# the surface of the PV. However, irradiance sensors typically only measure
+# global horizontal irradiance, GHI, therefore correlations to estimate the
+# diffuse fraction the GHI can be used to resolve the diffuse and beam
+# components.
+
+# This example demonstrates how to use :py:meth:`pvlib.irradiance.erbs`,
+# :py:meth:`pvlib.irradiance.boland`, and several other methods of varying
+# complexity
+
+
+from datetime import datetime
+import pathlib
+from matplotlib import pyplot as plt
+import numpy as np
+import pandas as pd
+from pvlib.iotools import read_tmy3
+from pvlib.solarposition import get_solarposition
+from pvlib import irradiance, solarposition
+import pvlib
+
+# get full path to the data directory
+DATA_DIR = pathlib.Path(pvlib.__file__).parent / 'data'
+
+# get TMY3 data with rain
+greensboro, metadata = read_tmy3(DATA_DIR / '723170TYA.CSV', coerce_year=1990)
+solpos = get_solarposition(
+    greensboro.index, latitude=metadata['latitude'],
+    longitude=metadata['longitude'], altitude=metadata['altitude'],
+    pressure=greensboro.Pressure*100, temperature=greensboro.DryBulb)
+
+# %%
+# DISC is an NREL function from 1987 that uses an empirical relation between
+# GHI and clearness index
+out_disc = irradiance.disc(
+    greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100)
+out_disc = out_disc.rename(columns={'dni': 'dni_disc'})
+out_disc['dhi_disc'] = (
+    greensboro.GHI
+    - out_disc.dni_disc*np.cos(np.radians(solpos.apparent_zenith)))
+
+# %%
+# next is erbs
+out_erbs = irradiance.erbs(greensboro.GHI, solpos.zenith, greensboro.index)
+out_erbs = out_erbs.rename(columns={'dni': 'dni_erbs', 'dhi': 'dhi_erbs'})
+# %%
+JAN6AM,JAN6PM = '1990-01-04 00:00:00-05:00', '1990-01-07 23:59:59-05:00'
+f, ax = plt.subplots(2, 1, figsize=(8, 6), sharex=True)
+dni_renames = {'DNI': 'TMY', 'dni_disc': 'DISC', 'dni_erbs': 'Erbs'}
+dni = pd.concat(
+    [greensboro.DNI, out_disc.dni_disc, out_erbs.dni_erbs], axis=1)
+dni = dni.rename(columns=dni_renames)
+dni[JAN6AM:JAN6PM].plot(ax=ax[0])
+ax[0].grid(which="both")
+ax[0].set_ylabel('DNI $[W/m^2]$')
+dhi_renames = {'DHI': 'TMY', 'dhi_disc': 'DISC', 'dhi_erbs': 'Erbs'}
+dhi = pd.concat(
+    [greensboro.DHI, out_disc.dhi_disc, out_erbs.dhi_erbs], axis=1)
+dhi = dhi.rename(columns=dhi_renames)
+dhi[JAN6AM:JAN6PM].plot(ax=ax[1])
+ax[1].grid(which="both")
+ax[1].set_ylabel('DHI $[W/m^2]$')
+f.tight_layout()
+# %%
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 0bec078c2d..8da486c35e 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -2190,22 +2190,15 @@ def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
     r"""
     Estimate DNI and DHI from GHI using the Boland clearness index model.
 
-    The Boland model [1]_, [2]_ estimates the diffuse fraction DF from global
-    horizontal irradiance through an empirical relationship between DF
-    and the ratio of GHI to extraterrestrial irradiance, Kt. The
-    function uses the diffuse fraction to compute DHI as
+    The Boland model [1]_, [2]_ estimates the diffuse fraction, DF, from global
+    horizontal irradiance, GHI, through an empirical relationship between DF
+    and the ratio of GHI to extraterrestrial irradiance or clearness index, kt.
 
     .. math::
 
-        DHI = DF \times GHI
-
-    DNI is then estimated as
+        \mathit{DF} = \frac{1}{1 + \exp\left(-5 + 8.6 k_t\right)}
 
-    .. math::
-
-        DNI = (GHI - DHI)/\cos(Z)
-
-    where Z is the zenith angle.
+    where :math:`k_t` is the clearness index.
 
     Parameters
     ----------

From a69333bc04559ca33cca535f2af4b9c9cead6b81 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <mark.mikofski@dnvgl.com>
Date: Fri, 26 Feb 2021 01:17:28 -0800
Subject: [PATCH 04/30] finish example

---
 docs/examples/plot_diffuse_fraction.py | 172 ++++++++++++++++++++-----
 1 file changed, 142 insertions(+), 30 deletions(-)

diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/plot_diffuse_fraction.py
index 80a6612177..b1d8ca065c 100644
--- a/docs/examples/plot_diffuse_fraction.py
+++ b/docs/examples/plot_diffuse_fraction.py
@@ -5,43 +5,50 @@
 Comparison of diffuse fraction estimation methods used to derive direct and
 diffuse components from measured global horizontal irradiance.
 """
-# %%
-# PV systems are often tilted to optimize performance. Determining the total
-# irradiance incident on the plane of the array requires transposing the
-# diffuse component, because the entire isotropic sky dome is not visible to
-# the surface of the PV. However, irradiance sensors typically only measure
-# global horizontal irradiance, GHI, therefore correlations to estimate the
-# diffuse fraction the GHI can be used to resolve the diffuse and beam
-# components.
-
-# This example demonstrates how to use :py:meth:`pvlib.irradiance.erbs`,
-# :py:meth:`pvlib.irradiance.boland`, and several other methods of varying
-# complexity
 
+# %%
+# This example demonstrates how to use diffuse fraction estimation methods to
+# obtain direct and diffuse components from measured global horizontal
+# irradiance (GHI). PV systems are often tilted to optimize performance, so the
+# entire diffuse sky dome may not be visible to the PV surface. Determining the
+# total irradiance incident on the plane of the array requires transposing the
+# diffuse component, but irradiance sensors such as pyranometers typically only
+# measure GHI. Therefore pvlib provides several correlations to estimate the
+# diffuse fraction of the GHI, that can be used to resolve the diffuse and
+# direct components.
 
-from datetime import datetime
 import pathlib
 from matplotlib import pyplot as plt
 import numpy as np
 import pandas as pd
 from pvlib.iotools import read_tmy3
 from pvlib.solarposition import get_solarposition
-from pvlib import irradiance, solarposition
+from pvlib import irradiance
 import pvlib
 
-# get full path to the data directory
+# For this example we use the Greensboro, North Carolina, TMY3 file which is
+# in the pvlib data directory. TMY3 are made from the median months from years
+# of data measured from 1990 to 2010. Therefore we change the timestamps to a
+# common year, 1990.
 DATA_DIR = pathlib.Path(pvlib.__file__).parent / 'data'
-
-# get TMY3 data with rain
 greensboro, metadata = read_tmy3(DATA_DIR / '723170TYA.CSV', coerce_year=1990)
+
+# Many of the diffuse fraction estimation methods require the "true" zenith, so
+# we calculate the solar positions for the 1990 at Greensboro, NC.
 solpos = get_solarposition(
     greensboro.index, latitude=metadata['latitude'],
     longitude=metadata['longitude'], altitude=metadata['altitude'],
     pressure=greensboro.Pressure*100, temperature=greensboro.DryBulb)
 
 # %%
-# DISC is an NREL function from 1987 that uses an empirical relation between
-# GHI and clearness index
+# DISC
+# ----
+#
+# DISC :py:meth:`~pvlib.irradiance.disc` is an empirical correlation developed
+# at SERI (now NREL) in 1987. The direct normal irradiance (DNI) is related to
+# clearness index (kt) by two polynomials split at kt = 0.6, then combined with
+# an exponential relation with airmass.
+
 out_disc = irradiance.disc(
     greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100)
 out_disc = out_disc.rename(columns={'dni': 'dni_disc'})
@@ -50,25 +57,130 @@
     - out_disc.dni_disc*np.cos(np.radians(solpos.apparent_zenith)))
 
 # %%
-# next is erbs
+# DIRINT
+# ------
+#
+# DIRINT :py:meth:`~pvlib.irradiance.dirint` is a modification of DISC
+# developed by Richard Perez and Pierre Ineichen in 1992.
+
+dni_dirint = irradiance.dirint(
+    greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100,
+    temp_dew=greensboro.DewPoint)
+dhi_dirint = (
+    greensboro.GHI
+    - dni_dirint*np.cos(np.radians(solpos.apparent_zenith)))
+out_dirint = pd.DataFrame(
+    {'dni_dirint': dni_dirint, 'dhi_dirint': dhi_dirint},
+    index=greensboro.index)
+
+# %%
+# Erbs
+# ----
+#
+# The Erbs method, :py:meth:`~pvlib.irradiance.erbs` developed by Daryl Gregory
+# Erbs at the University of Wisconsin in 1982 is a piecewise correlation that
+# splits kt into 3 regions: linear for kt <= 0.22, a 4th order polynomial
+# between 0.22 < kt <= 0.8, and a horizontal line for kt > 0.8.
+
 out_erbs = irradiance.erbs(greensboro.GHI, solpos.zenith, greensboro.index)
 out_erbs = out_erbs.rename(columns={'dni': 'dni_erbs', 'dhi': 'dhi_erbs'})
+
+# %%
+# Boland
+# ----
+#
+# The Boland method, :py:meth:`~pvlib.irradiance.boland` is a single logistic
+# exponential correlation that is continuously differentiable and bounded
+# between zero and one.
+
+out_boland = irradiance.boland(greensboro.GHI, solpos.zenith, greensboro.index)
+out_boland = out_boland.rename(
+    columns={'dni': 'dni_boland', 'dhi': 'dhi_boland'})
+
+# %%
+# Combine everything together.
+
+dni_renames = {
+    'DNI': 'TMY', 'dni_disc': 'DISC', 'dni_dirint': 'DIRINT',
+    'dni_erbs': 'Erbs', 'dni_boland': 'Boland'}
+dni = [
+    greensboro.DNI, out_disc.dni_disc, out_dirint.dni_dirint,
+    out_erbs.dni_erbs, out_boland.dni_boland]
+dni = pd.concat(dni, axis=1).rename(columns=dni_renames)
+dhi_renames = {
+    'DHI': 'TMY', 'dhi_disc': 'DISC', 'dhi_dirint': 'DIRINT',
+    'dhi_erbs': 'Erbs', 'dhi_boland': 'Boland'}
+dhi = [
+    greensboro.DHI, out_disc.dhi_disc, out_dirint.dhi_dirint,
+    out_erbs.dhi_erbs, out_boland.dhi_boland]
+dhi = pd.concat(dhi, axis=1).rename(columns=dhi_renames)
+ghi_kt = pd.concat([greensboro.GHI/1000.0, out_erbs.kt], axis=1)
+
 # %%
-JAN6AM,JAN6PM = '1990-01-04 00:00:00-05:00', '1990-01-07 23:59:59-05:00'
-f, ax = plt.subplots(2, 1, figsize=(8, 6), sharex=True)
-dni_renames = {'DNI': 'TMY', 'dni_disc': 'DISC', 'dni_erbs': 'Erbs'}
-dni = pd.concat(
-    [greensboro.DNI, out_disc.dni_disc, out_erbs.dni_erbs], axis=1)
-dni = dni.rename(columns=dni_renames)
+# Finally, let's plot them for a few winter days and compare
+
+JAN6AM, JAN6PM = '1990-01-04 00:00:00-05:00', '1990-01-07 23:59:59-05:00'
+f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
 dni[JAN6AM:JAN6PM].plot(ax=ax[0])
 ax[0].grid(which="both")
 ax[0].set_ylabel('DNI $[W/m^2]$')
-dhi_renames = {'DHI': 'TMY', 'dhi_disc': 'DISC', 'dhi_erbs': 'Erbs'}
-dhi = pd.concat(
-    [greensboro.DHI, out_disc.dhi_disc, out_erbs.dhi_erbs], axis=1)
-dhi = dhi.rename(columns=dhi_renames)
+ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
 dhi[JAN6AM:JAN6PM].plot(ax=ax[1])
 ax[1].grid(which="both")
 ax[1].set_ylabel('DHI $[W/m^2]$')
+ghi_kt[JAN6AM:JAN6PM].plot(ax=ax[2])
+ax[2].grid(which='both')
+ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
+f.tight_layout()
+
+# %%
+# And a few spring days ...
+
+APR6AM, APR6PM = '1990-04-04 00:00:00-05:00', '1990-04-07 23:59:59-05:00'
+f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
+dni[APR6AM:APR6PM].plot(ax=ax[0])
+ax[0].grid(which="both")
+ax[0].set_ylabel('DNI $[W/m^2]$')
+ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
+dhi[APR6AM:APR6PM].plot(ax=ax[1])
+ax[1].grid(which="both")
+ax[1].set_ylabel('DHI $[W/m^2]$')
+ghi_kt[APR6AM:APR6PM].plot(ax=ax[2])
+ax[2].grid(which='both')
+ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
 f.tight_layout()
+
+# %%
+# And few summer days to finish off the seasons.
+
+JUL6AM, JUL6PM = '1990-07-04 00:00:00-05:00', '1990-07-07 23:59:59-05:00'
+f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
+dni[JUL6AM:JUL6PM].plot(ax=ax[0])
+ax[0].grid(which="both")
+ax[0].set_ylabel('DNI $[W/m^2]$')
+ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
+dhi[JUL6AM:JUL6PM].plot(ax=ax[1])
+ax[1].grid(which="both")
+ax[1].set_ylabel('DHI $[W/m^2]$')
+ghi_kt[JUL6AM:JUL6PM].plot(ax=ax[2])
+ax[2].grid(which='both')
+ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
+f.tight_layout()
+
 # %%
+# Conclusion
+# ----------
+# The Erbs and Boland are correlations with only kt, which is derived from the
+# horizontal component of the extra-terrestrial irradiance. Therefore at low
+# sun elevation (zenith ~ 90-deg), especially near sunset, this causes kt to
+# explode as the denominator approaches zero. This is controlled in pvlib by
+# setting ``min_cos_zenith`` and ``max_clearness_index`` which each have
+# reasonable defaults, but there are still concerning spikes at sunset for Jan.
+# 5th & 7th, April 4th, 5th, & 7th, and July 6th & 7th. The DISC & DIRINT
+# methods differ from Erbs and Boland be including airmass, which seems to
+# reduce DNI spikes over 1000[W/m^2], but still have errors at other times.
+#
+# Another difference is that DISC & DIRINT return DNI whereas Erbs & Boland
+# calculate the diffuse fraction which is then used to derive DNI from GHI and
+# the solar zenith, which exacerbates errors at low sun elevation due to the
+# relation: DNI = GHI*(1 - DF)/cos(zenith). 

From 2faf9ab8555171808c4781a5ad80ada7ea159aa3 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <mark.mikofski@dnvgl.com>
Date: Fri, 26 Feb 2021 01:19:33 -0800
Subject: [PATCH 05/30] trailing whitespace

---
 docs/examples/plot_diffuse_fraction.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/plot_diffuse_fraction.py
index b1d8ca065c..4d98600b64 100644
--- a/docs/examples/plot_diffuse_fraction.py
+++ b/docs/examples/plot_diffuse_fraction.py
@@ -183,4 +183,4 @@
 # Another difference is that DISC & DIRINT return DNI whereas Erbs & Boland
 # calculate the diffuse fraction which is then used to derive DNI from GHI and
 # the solar zenith, which exacerbates errors at low sun elevation due to the
-# relation: DNI = GHI*(1 - DF)/cos(zenith). 
+# relation: DNI = GHI*(1 - DF)/cos(zenith).

From 6c8cb94a6265d8511cdd0509979c7c37e9c602dd Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 16:28:56 -0800
Subject: [PATCH 06/30] respond to comments

- put Badescu textbook before paper
- revise wording defining kt
---
 pvlib/irradiance.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index b8a69d0b8c..fc198eec6e 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -2271,7 +2271,8 @@ def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
 
     The Boland model [1]_, [2]_ estimates the diffuse fraction, DF, from global
     horizontal irradiance, GHI, through an empirical relationship between DF
-    and the ratio of GHI to extraterrestrial irradiance or clearness index, kt.
+    and the clearness index, :math:`k_t`, the ratio of GHI to extraterrestrial
+    irradiance.
 
     .. math::
 
@@ -2308,13 +2309,13 @@ def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
 
     References
     ----------
-    .. [1] John Boland, Lynne Scott, and Mark Luther, Modelling the diffuse
+    .. [1] J. Boland, B. Ridley (2008) Models of Diffuse Solar Fraction. In:
+       Badescu V. (eds) Modeling Solar Radiation at the Earth’s Surface.
+       Springer, Berlin, Heidelberg. :doi:`10.1007/978-3-540-77455-6_8`
+    .. [2] John Boland, Lynne Scott, and Mark Luther, Modelling the diffuse
        fraction of global solar radiation on a horizontal surface,
        Environmetrics 12(2), pp 103-116, 2001,
        :doi:`10.1002/1099-095X(200103)12:2%3C103::AID-ENV447%3E3.0.CO;2-2`
-    .. [2] J. Boland, B. Ridley (2008) Models of Diffuse Solar Fraction. In:
-       Badescu V. (eds) Modeling Solar Radiation at the Earth’s Surface.
-       Springer, Berlin, Heidelberg. :doi:`10.1007/978-3-540-77455-6_8`
 
     See also
     --------

From 21da1898d33b2a660869503c41ee2684c77e8572 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 16:58:15 -0800
Subject: [PATCH 07/30] respond to comments

- remove redundant kt definition
- use intersphinx to link to pandas and numpy types
---
 pvlib/irradiance.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index fc198eec6e..f732164ef4 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -2278,15 +2278,14 @@ def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
 
         \mathit{DF} = \frac{1}{1 + \exp\left(-5 + 8.6 k_t\right)}
 
-    where :math:`k_t` is the clearness index.
 
     Parameters
     ----------
     ghi: numeric
-        Global horizontal irradiance in W/m^2.
+        Global horizontal irradiance. [W/m^2]
     zenith: numeric
         True (not refraction-corrected) zenith angles in decimal degrees.
-    datetime_or_doy : int, float, array, pd.DatetimeIndex
+    datetime_or_doy : int, float, numpy.ndarray, pandas.DatetimeIndex
         Day of year or array of days of year e.g.
         pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
     min_cos_zenith : numeric, default 0.065

From cfee27fa97531249f39586cb03b4bcf2494995b6 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 17:27:43 -0800
Subject: [PATCH 08/30] Update docs/examples/plot_diffuse_fraction.py

responding to comments, simplify wording L12-15 to omit "systems often tilted to optimize performance..."

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
---
 docs/examples/plot_diffuse_fraction.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/plot_diffuse_fraction.py
index 4d98600b64..84a17e78a2 100644
--- a/docs/examples/plot_diffuse_fraction.py
+++ b/docs/examples/plot_diffuse_fraction.py
@@ -9,10 +9,7 @@
 # %%
 # This example demonstrates how to use diffuse fraction estimation methods to
 # obtain direct and diffuse components from measured global horizontal
-# irradiance (GHI). PV systems are often tilted to optimize performance, so the
-# entire diffuse sky dome may not be visible to the PV surface. Determining the
-# total irradiance incident on the plane of the array requires transposing the
-# diffuse component, but irradiance sensors such as pyranometers typically only
+# irradiance (GHI). Irradiance sensors such as pyranometers typically only
 # measure GHI. Therefore pvlib provides several correlations to estimate the
 # diffuse fraction of the GHI, that can be used to resolve the diffuse and
 # direct components.

From 9cc772a0c21a457f0f773887c16b5cc98ea6e57a Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 17:39:27 -0800
Subject: [PATCH 09/30] Update docs/examples/plot_diffuse_fraction.py

respond to comments, reword decomposition example intro

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
---
 docs/examples/plot_diffuse_fraction.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/plot_diffuse_fraction.py
index 84a17e78a2..cd0f449188 100644
--- a/docs/examples/plot_diffuse_fraction.py
+++ b/docs/examples/plot_diffuse_fraction.py
@@ -10,9 +10,10 @@
 # This example demonstrates how to use diffuse fraction estimation methods to
 # obtain direct and diffuse components from measured global horizontal
 # irradiance (GHI). Irradiance sensors such as pyranometers typically only
-# measure GHI. Therefore pvlib provides several correlations to estimate the
-# diffuse fraction of the GHI, that can be used to resolve the diffuse and
-# direct components.
+# measure GHI. pvlib provides several functions that can be used
+# to separate GHI into the diffuse and direct components. The separate
+# components are needed to estimate the total irradiance on a tilted
+# surface.
 
 import pathlib
 from matplotlib import pyplot as plt

From 9bdd1bb8f55f81ba4a109c47c43101ba6a8b6ad8 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 17:51:55 -0800
Subject: [PATCH 10/30] respond to comments

- add intro before functions in example
---
 docs/examples/plot_diffuse_fraction.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/plot_diffuse_fraction.py
index cd0f449188..04973065f2 100644
--- a/docs/examples/plot_diffuse_fraction.py
+++ b/docs/examples/plot_diffuse_fraction.py
@@ -10,10 +10,9 @@
 # This example demonstrates how to use diffuse fraction estimation methods to
 # obtain direct and diffuse components from measured global horizontal
 # irradiance (GHI). Irradiance sensors such as pyranometers typically only
-# measure GHI. pvlib provides several functions that can be used
-# to separate GHI into the diffuse and direct components. The separate
-# components are needed to estimate the total irradiance on a tilted
-# surface.
+# measure GHI. pvlib provides several functions that can be used to separate
+# GHI into the diffuse and direct components. The separate components are
+# needed to estimate the total irradiance on a tilted surface.
 
 import pathlib
 from matplotlib import pyplot as plt
@@ -38,6 +37,10 @@
     longitude=metadata['longitude'], altitude=metadata['altitude'],
     pressure=greensboro.Pressure*100, temperature=greensboro.DryBulb)
 
+# %%
+# Methods for separating DHI into diffuse and direct components include:
+# **DISC**, **DIRINT**, **Erbs** and **Boland**.
+
 # %%
 # DISC
 # ----
@@ -112,7 +115,7 @@
     greensboro.DHI, out_disc.dhi_disc, out_dirint.dhi_dirint,
     out_erbs.dhi_erbs, out_boland.dhi_boland]
 dhi = pd.concat(dhi, axis=1).rename(columns=dhi_renames)
-ghi_kt = pd.concat([greensboro.GHI/1000.0, out_erbs.kt], axis=1)
+ghi_kt = pd.concat([greensboro.GHI/1367.0, out_erbs.kt], axis=1)
 
 # %%
 # Finally, let's plot them for a few winter days and compare

From cbb61bd6b57280109f862abb29ecf8860326b1ad Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 18:05:21 -0800
Subject: [PATCH 11/30] respond to comments in diffuse fraction example

- fix use :py:func: to xref functions in docs, instead of :meth:
- add comment to explain conversion of mbars to Pa
- use Gsc not E0 in kt comparison plot
---
 docs/examples/plot_diffuse_fraction.py | 15 ++++++++-------
 docs/sphinx/source/whatsnew/v0.9.5.rst |  2 +-
 2 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/plot_diffuse_fraction.py
index 04973065f2..f78e1cf567 100644
--- a/docs/examples/plot_diffuse_fraction.py
+++ b/docs/examples/plot_diffuse_fraction.py
@@ -35,7 +35,8 @@
 solpos = get_solarposition(
     greensboro.index, latitude=metadata['latitude'],
     longitude=metadata['longitude'], altitude=metadata['altitude'],
-    pressure=greensboro.Pressure*100, temperature=greensboro.DryBulb)
+    pressure=greensboro.Pressure*100,  # convert from millibar to Pa
+    temperature=greensboro.DryBulb)
 
 # %%
 # Methods for separating DHI into diffuse and direct components include:
@@ -45,7 +46,7 @@
 # DISC
 # ----
 #
-# DISC :py:meth:`~pvlib.irradiance.disc` is an empirical correlation developed
+# DISC :py:func:`~pvlib.irradiance.disc` is an empirical correlation developed
 # at SERI (now NREL) in 1987. The direct normal irradiance (DNI) is related to
 # clearness index (kt) by two polynomials split at kt = 0.6, then combined with
 # an exponential relation with airmass.
@@ -61,7 +62,7 @@
 # DIRINT
 # ------
 #
-# DIRINT :py:meth:`~pvlib.irradiance.dirint` is a modification of DISC
+# DIRINT :py:func:`~pvlib.irradiance.dirint` is a modification of DISC
 # developed by Richard Perez and Pierre Ineichen in 1992.
 
 dni_dirint = irradiance.dirint(
@@ -78,7 +79,7 @@
 # Erbs
 # ----
 #
-# The Erbs method, :py:meth:`~pvlib.irradiance.erbs` developed by Daryl Gregory
+# The Erbs method, :py:func:`~pvlib.irradiance.erbs` developed by Daryl Gregory
 # Erbs at the University of Wisconsin in 1982 is a piecewise correlation that
 # splits kt into 3 regions: linear for kt <= 0.22, a 4th order polynomial
 # between 0.22 < kt <= 0.8, and a horizontal line for kt > 0.8.
@@ -90,7 +91,7 @@
 # Boland
 # ----
 #
-# The Boland method, :py:meth:`~pvlib.irradiance.boland` is a single logistic
+# The Boland method, :py:func:`~pvlib.irradiance.boland` is a single logistic
 # exponential correlation that is continuously differentiable and bounded
 # between zero and one.
 
@@ -115,7 +116,7 @@
     greensboro.DHI, out_disc.dhi_disc, out_dirint.dhi_dirint,
     out_erbs.dhi_erbs, out_boland.dhi_boland]
 dhi = pd.concat(dhi, axis=1).rename(columns=dhi_renames)
-ghi_kt = pd.concat([greensboro.GHI/1367.0, out_erbs.kt], axis=1)
+ghi_kt = pd.concat([greensboro.GHI/1366.1, out_erbs.kt], axis=1)
 
 # %%
 # Finally, let's plot them for a few winter days and compare
@@ -131,7 +132,7 @@
 ax[1].set_ylabel('DHI $[W/m^2]$')
 ghi_kt[JAN6AM:JAN6PM].plot(ax=ax[2])
 ax[2].grid(which='both')
-ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
+ax[2].set_ylabel(r'$\frac{GHI}{G_{SC}}, k_t$')
 f.tight_layout()
 
 # %%
diff --git a/docs/sphinx/source/whatsnew/v0.9.5.rst b/docs/sphinx/source/whatsnew/v0.9.5.rst
index 5093554d8b..89f8deb861 100644
--- a/docs/sphinx/source/whatsnew/v0.9.5.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.5.rst
@@ -20,7 +20,7 @@ Enhancements
   :py:func:`pvlib.snow.loss_townsend` (:issue:`1636`, :pull:`1653`)
 * Added optional ``n_ar`` parameter to :py:func:`pvlib.iam.physical` to
   support an anti-reflective coating. (:issue:`1501`, :pull:`1616`)
-* :py:meth:`~pvlib.irradiance.boland` is another diffuse fraction, DF,
+* :py:func:`~pvlib.irradiance.boland` is another diffuse fraction, DF,
   estimation method similar to Erbs but uses a single logistic exponential
   correlation between DF and clearness index, kt, that is continuously
   differentiable and bounded between zero and one. (:pull:`1179`)

From c9ce7b87a6235c38b0e111c958583e4a518331dd Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 20:03:27 -0800
Subject: [PATCH 12/30] respond to comments

- add section to conclusions to warn users that TMY3 and NSRDB are also models,
 therefore not to assume differences are errors,
 and link to TMY3 & NSRDB documentation
- use implicit targets to link to DISC, DIRINT, Erbs, & Boland sections
- reverse change of scaled GHI by E0 to Gsc (aka: solar constant)
- use math directive for DNI formula
---
 docs/examples/plot_diffuse_fraction.py | 30 +++++++++++++++++++-------
 1 file changed, 22 insertions(+), 8 deletions(-)

diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/plot_diffuse_fraction.py
index f78e1cf567..9ab3d06300 100644
--- a/docs/examples/plot_diffuse_fraction.py
+++ b/docs/examples/plot_diffuse_fraction.py
@@ -39,8 +39,10 @@
     temperature=greensboro.DryBulb)
 
 # %%
+# pvlib Decomposition Functions
+# -----------------------------
 # Methods for separating DHI into diffuse and direct components include:
-# **DISC**, **DIRINT**, **Erbs** and **Boland**.
+# `DISC`_, `DIRINT`_, `Erbs`_, and `Boland`_.
 
 # %%
 # DISC
@@ -89,7 +91,7 @@
 
 # %%
 # Boland
-# ----
+# ------
 #
 # The Boland method, :py:func:`~pvlib.irradiance.boland` is a single logistic
 # exponential correlation that is continuously differentiable and bounded
@@ -103,20 +105,20 @@
 # Combine everything together.
 
 dni_renames = {
-    'DNI': 'TMY', 'dni_disc': 'DISC', 'dni_dirint': 'DIRINT',
+    'DNI': 'TMY3', 'dni_disc': 'DISC', 'dni_dirint': 'DIRINT',
     'dni_erbs': 'Erbs', 'dni_boland': 'Boland'}
 dni = [
     greensboro.DNI, out_disc.dni_disc, out_dirint.dni_dirint,
     out_erbs.dni_erbs, out_boland.dni_boland]
 dni = pd.concat(dni, axis=1).rename(columns=dni_renames)
 dhi_renames = {
-    'DHI': 'TMY', 'dhi_disc': 'DISC', 'dhi_dirint': 'DIRINT',
+    'DHI': 'TMY3', 'dhi_disc': 'DISC', 'dhi_dirint': 'DIRINT',
     'dhi_erbs': 'Erbs', 'dhi_boland': 'Boland'}
 dhi = [
     greensboro.DHI, out_disc.dhi_disc, out_dirint.dhi_dirint,
     out_erbs.dhi_erbs, out_boland.dhi_boland]
 dhi = pd.concat(dhi, axis=1).rename(columns=dhi_renames)
-ghi_kt = pd.concat([greensboro.GHI/1366.1, out_erbs.kt], axis=1)
+ghi_kt = pd.concat([greensboro.GHI/1000.0, out_erbs.kt], axis=1)
 
 # %%
 # Finally, let's plot them for a few winter days and compare
@@ -132,7 +134,7 @@
 ax[1].set_ylabel('DHI $[W/m^2]$')
 ghi_kt[JAN6AM:JAN6PM].plot(ax=ax[2])
 ax[2].grid(which='both')
-ax[2].set_ylabel(r'$\frac{GHI}{G_{SC}}, k_t$')
+ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
 f.tight_layout()
 
 # %%
@@ -172,6 +174,14 @@
 # %%
 # Conclusion
 # ----------
+# This example has compared several decomposition models to a TMY3 file for
+# Greensboro, North Carolina. However, DNI and DHI in the TMY3 file are
+# themselves the output of a model (either METSTAT or SUNY), and so differences
+# between *e.g.* DISC output and TMY3 shouldn't be regarded as an error in the
+# DISC model. Therefore, it's not a reasonable expectation to assume that the
+# four models should reproduce the TMY3 values. We refer those interested to
+# the `TMY3`_ and `NSRDB`_ user manuals.
+#
 # The Erbs and Boland are correlations with only kt, which is derived from the
 # horizontal component of the extra-terrestrial irradiance. Therefore at low
 # sun elevation (zenith ~ 90-deg), especially near sunset, this causes kt to
@@ -179,10 +189,14 @@
 # setting ``min_cos_zenith`` and ``max_clearness_index`` which each have
 # reasonable defaults, but there are still concerning spikes at sunset for Jan.
 # 5th & 7th, April 4th, 5th, & 7th, and July 6th & 7th. The DISC & DIRINT
-# methods differ from Erbs and Boland be including airmass, which seems to
+# methods differ from Erbs and Boland by including airmass, which seems to
 # reduce DNI spikes over 1000[W/m^2], but still have errors at other times.
 #
 # Another difference is that DISC & DIRINT return DNI whereas Erbs & Boland
 # calculate the diffuse fraction which is then used to derive DNI from GHI and
 # the solar zenith, which exacerbates errors at low sun elevation due to the
-# relation: DNI = GHI*(1 - DF)/cos(zenith).
+# relation:
+# :math:`DNI = GHI \frac{1 - \mathit{DF}}{\cos \left(\mathit{zenith} \right)}`.
+#
+# .. _TMY3: https://www.nrel.gov/docs/fy08osti/43156.pdf
+# .. _NSRDB: https://www.nrel.gov/docs/fy12osti/54824.pdf
\ No newline at end of file

From e23662dedc28e26dce60b473eeec73b0d2f10627 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 20:53:51 -0800
Subject: [PATCH 13/30] respond to comments

- revise & condense conclusions, don't refer to differences as errors
 without operational data
- add header before plots section, explain why comparing normalized GHI
- add subheaders for seasonal plots
---
 docs/examples/plot_diffuse_fraction.py | 55 +++++++++++++++-----------
 1 file changed, 33 insertions(+), 22 deletions(-)

diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/plot_diffuse_fraction.py
index 9ab3d06300..73acb4ce49 100644
--- a/docs/examples/plot_diffuse_fraction.py
+++ b/docs/examples/plot_diffuse_fraction.py
@@ -102,7 +102,15 @@
     columns={'dni': 'dni_boland', 'dhi': 'dhi_boland'})
 
 # %%
-# Combine everything together.
+# Comparison Plots
+# ----------------
+# In the plots below we compare the four decomposition models to the TMY3 file
+# for Greensboro, North Carolina. We also compare the clearness index, kt, with
+# GHI normalized by a reference irradiance, E0 = 1000 [W/m^2], to highlight
+# spikes caused when cosine of zenith approaches zero, particularly at sunset.
+#
+# First we combine the dataframes for the decomposition models and the TMY3
+# file together to make plotting easier.
 
 dni_renames = {
     'DNI': 'TMY3', 'dni_disc': 'DISC', 'dni_dirint': 'DIRINT',
@@ -121,6 +129,8 @@
 ghi_kt = pd.concat([greensboro.GHI/1000.0, out_erbs.kt], axis=1)
 
 # %%
+# Winter
+# ++++++
 # Finally, let's plot them for a few winter days and compare
 
 JAN6AM, JAN6PM = '1990-01-04 00:00:00-05:00', '1990-01-07 23:59:59-05:00'
@@ -138,6 +148,8 @@
 f.tight_layout()
 
 # %%
+# Spring
+# ++++++
 # And a few spring days ...
 
 APR6AM, APR6PM = '1990-04-04 00:00:00-05:00', '1990-04-07 23:59:59-05:00'
@@ -155,6 +167,8 @@
 f.tight_layout()
 
 # %%
+# Summer
+# ++++++
 # And few summer days to finish off the seasons.
 
 JUL6AM, JUL6PM = '1990-07-04 00:00:00-05:00', '1990-07-07 23:59:59-05:00'
@@ -174,29 +188,26 @@
 # %%
 # Conclusion
 # ----------
-# This example has compared several decomposition models to a TMY3 file for
-# Greensboro, North Carolina. However, DNI and DHI in the TMY3 file are
-# themselves the output of a model (either METSTAT or SUNY), and so differences
-# between *e.g.* DISC output and TMY3 shouldn't be regarded as an error in the
-# DISC model. Therefore, it's not a reasonable expectation to assume that the
-# four models should reproduce the TMY3 values. We refer those interested to
-# the `TMY3`_ and `NSRDB`_ user manuals.
+# This example compares several decomposition models to a TMY3 file for
+# Greensboro, North Carolina. However, DNI and DHI in TMY3 files are themselves
+# the output of models (either METSTAT or SUNY), and so differences between
+# *e.g.* DISC output and the TMY3 file shouldn't be regarded as errors, and
+# it's not a reasonable expectation to assume that the four models should
+# reproduce the TMY3 values. We refer those interested to the `TMY3`_ and
+# `NSRDB`_ user manuals.
 #
-# The Erbs and Boland are correlations with only kt, which is derived from the
-# horizontal component of the extra-terrestrial irradiance. Therefore at low
-# sun elevation (zenith ~ 90-deg), especially near sunset, this causes kt to
-# explode as the denominator approaches zero. This is controlled in pvlib by
-# setting ``min_cos_zenith`` and ``max_clearness_index`` which each have
-# reasonable defaults, but there are still concerning spikes at sunset for Jan.
-# 5th & 7th, April 4th, 5th, & 7th, and July 6th & 7th. The DISC & DIRINT
-# methods differ from Erbs and Boland by including airmass, which seems to
-# reduce DNI spikes over 1000[W/m^2], but still have errors at other times.
+# The Erbs and Boland models are correlations with only kt, which is derived
+# from the horizontal component of the extra-terrestrial irradiance. At low sun
+# elevation (zenith near 90 degrees), especially near sunset, kt can explode
+# because the denominator (extra-terrestrial irradiance) approaches zero. In
+# pvlib this behavior is moderated by ``min_cos_zenith`` and
+# ``max_clearness_index`` which each have reasonable defaults. Even so, near
+# sunset there are still spikes in kt and DNI from Erbs and Boland for Jan. 5th
+# & 7th, April 4th, 5th, & 7th, and July 6th & 7th.
 #
-# Another difference is that DISC & DIRINT return DNI whereas Erbs & Boland
-# calculate the diffuse fraction which is then used to derive DNI from GHI and
-# the solar zenith, which exacerbates errors at low sun elevation due to the
-# relation:
-# :math:`DNI = GHI \frac{1 - \mathit{DF}}{\cos \left(\mathit{zenith} \right)}`.
+# By contrast, the DISC and DIRINT methods estimate DNI first by means of
+# correlations, which include additional variables such as airmass. These methods
+# seem to reduce DNI spikes over 1000 [W/m^2].
 #
 # .. _TMY3: https://www.nrel.gov/docs/fy08osti/43156.pdf
 # .. _NSRDB: https://www.nrel.gov/docs/fy12osti/54824.pdf
\ No newline at end of file

From 90f4434e7988faed3813064587f3e55d3579e1ac Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 21:07:09 -0800
Subject: [PATCH 14/30] fix stickler

---
 docs/examples/plot_diffuse_fraction.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/plot_diffuse_fraction.py
index 73acb4ce49..94786c999a 100644
--- a/docs/examples/plot_diffuse_fraction.py
+++ b/docs/examples/plot_diffuse_fraction.py
@@ -206,8 +206,8 @@
 # & 7th, April 4th, 5th, & 7th, and July 6th & 7th.
 #
 # By contrast, the DISC and DIRINT methods estimate DNI first by means of
-# correlations, which include additional variables such as airmass. These methods
-# seem to reduce DNI spikes over 1000 [W/m^2].
+# correlations, which include additional variables such as airmass. These
+# methods seem to reduce DNI spikes over 1000 [W/m^2].
 #
 # .. _TMY3: https://www.nrel.gov/docs/fy08osti/43156.pdf
-# .. _NSRDB: https://www.nrel.gov/docs/fy12osti/54824.pdf
\ No newline at end of file
+# .. _NSRDB: https://www.nrel.gov/docs/fy12osti/54824.pdf

From ccf40a8eb641b540c0744dcb205bad7e4ebbc660 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 23:03:20 -0800
Subject: [PATCH 15/30] update readthedocs.yml to v2

---
 readthedocs.yml | 44 +++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 39 insertions(+), 5 deletions(-)

diff --git a/readthedocs.yml b/readthedocs.yml
index fb2d1374bb..0d2c5c6b7b 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -1,7 +1,41 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Set the version of Python and other tools you might need
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.9"
+
+# Build documentation in the docs/ directory with Sphinx
+sphinx:
+  configuration: docs/sphinx/source/conf.py
+
+# If using Sphinx, optionally build your docs in additional formats such as PDF
+# formats:
+#    - pdf
+
+# Optionally declare the Python requirements required to build your docs
 python:
-    version: 3
-    # only use the packages specified in setup.py
-    use_system_site_packages: false
-    pip_install: true
+  install:
+  # - requirements: docs/requirements.txt
+  - method: pip
+    path: .
     extra_requirements:
-        - doc
\ No newline at end of file
+    - docs
+
+    # python.version deprecated in v2 replaced by build.tools.python
+    # version: 3  
+
+    # setup_py_install replaced by python.install in v2
+    # only use the packages specified in setup.py
+
+    # renamed to python.system_packages in v2
+    # use_system_site_packages: false
+
+    # replaced by python.install in v2
+    # pip_install: true

From b1027bc30c0d3f545958b172d3d63ef108fdc64d Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 23:08:18 -0800
Subject: [PATCH 16/30] oops, extra requires is doc, no s

---
 readthedocs.yml | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/readthedocs.yml b/readthedocs.yml
index 0d2c5c6b7b..92e35f6def 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -22,11 +22,11 @@ sphinx:
 # Optionally declare the Python requirements required to build your docs
 python:
   install:
-  # - requirements: docs/requirements.txt
-  - method: pip
-    path: .
-    extra_requirements:
-    - docs
+    # - requirements: docs/requirements.txt
+    - method: pip
+      path: .
+      extra_requirements:
+        - doc
 
     # python.version deprecated in v2 replaced by build.tools.python
     # version: 3  

From 28fab4f4f88834a70762b12a7b42b0433b0d3937 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 23:18:46 -0800
Subject: [PATCH 17/30] use requirements to pin requirements

---
 docs/requirements.txt | 17 +++++++++++++++++
 readthedocs.yml       |  4 +++-
 2 files changed, 20 insertions(+), 1 deletion(-)
 create mode 100644 docs/requirements.txt

diff --git a/docs/requirements.txt b/docs/requirements.txt
new file mode 100644
index 0000000000..5cc3090539
--- /dev/null
+++ b/docs/requirements.txt
@@ -0,0 +1,17 @@
+numpy>=1.16.0
+pandas>=0.25.0
+pytz
+requests
+scipy>=1.4.0
+h5py
+ipython
+matplotlib
+sphinx==4.5.0
+pydata-sphinx-theme==0.8.1
+sphinx-gallery
+docutils==0.15.2
+pillow
+netcdf4
+siphon
+sphinx-toggleprompt>=0.0.5
+pvfactors==1.5.2
diff --git a/readthedocs.yml b/readthedocs.yml
index 92e35f6def..534b6bde70 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -22,7 +22,6 @@ sphinx:
 # Optionally declare the Python requirements required to build your docs
 python:
   install:
-    # - requirements: docs/requirements.txt
     - method: pip
       path: .
       extra_requirements:
@@ -39,3 +38,6 @@ python:
 
     # replaced by python.install in v2
     # pip_install: true
+
+    # pin requirements (recommended)
+    - requirements: docs/requirements.txt

From 9715084d00888953ec9b1b2af86a914910aada9d Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 26 Feb 2023 23:44:21 -0800
Subject: [PATCH 18/30] install pvlib again?

---
 readthedocs.yml | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/readthedocs.yml b/readthedocs.yml
index 534b6bde70..aa6ba529a9 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -41,3 +41,6 @@ python:
 
     # pin requirements (recommended)
     - requirements: docs/requirements.txt
+
+    - method: pip
+      path: .

From 1ae35cdc27ffe082c8f67868caff8a5ea28a554a Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Mon, 27 Feb 2023 00:54:50 -0800
Subject: [PATCH 19/30] Revert "update readthedocs.yml to v2"

This reverts commit ccf40a8eb641b540c0744dcb205bad7e4ebbc660.
---
 docs/requirements.txt | 17 ---------------
 readthedocs.yml       | 49 +++++--------------------------------------
 2 files changed, 5 insertions(+), 61 deletions(-)
 delete mode 100644 docs/requirements.txt

diff --git a/docs/requirements.txt b/docs/requirements.txt
deleted file mode 100644
index 5cc3090539..0000000000
--- a/docs/requirements.txt
+++ /dev/null
@@ -1,17 +0,0 @@
-numpy>=1.16.0
-pandas>=0.25.0
-pytz
-requests
-scipy>=1.4.0
-h5py
-ipython
-matplotlib
-sphinx==4.5.0
-pydata-sphinx-theme==0.8.1
-sphinx-gallery
-docutils==0.15.2
-pillow
-netcdf4
-siphon
-sphinx-toggleprompt>=0.0.5
-pvfactors==1.5.2
diff --git a/readthedocs.yml b/readthedocs.yml
index aa6ba529a9..fb2d1374bb 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -1,46 +1,7 @@
-# .readthedocs.yaml
-# Read the Docs configuration file
-# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
-
-# Required
-version: 2
-
-# Set the version of Python and other tools you might need
-build:
-  os: ubuntu-22.04
-  tools:
-    python: "3.9"
-
-# Build documentation in the docs/ directory with Sphinx
-sphinx:
-  configuration: docs/sphinx/source/conf.py
-
-# If using Sphinx, optionally build your docs in additional formats such as PDF
-# formats:
-#    - pdf
-
-# Optionally declare the Python requirements required to build your docs
 python:
-  install:
-    - method: pip
-      path: .
-      extra_requirements:
-        - doc
-
-    # python.version deprecated in v2 replaced by build.tools.python
-    # version: 3  
-
-    # setup_py_install replaced by python.install in v2
+    version: 3
     # only use the packages specified in setup.py
-
-    # renamed to python.system_packages in v2
-    # use_system_site_packages: false
-
-    # replaced by python.install in v2
-    # pip_install: true
-
-    # pin requirements (recommended)
-    - requirements: docs/requirements.txt
-
-    - method: pip
-      path: .
+    use_system_site_packages: false
+    pip_install: true
+    extra_requirements:
+        - doc
\ No newline at end of file

From ecdb50827e725da04e04bd94c1f0fce2a599da29 Mon Sep 17 00:00:00 2001
From: Kevin Anderson <kevin.anderson@nrel.gov>
Date: Mon, 27 Feb 2023 09:22:09 -0500
Subject: [PATCH 20/30] update RTD config to fix shallow clone issue

---
 readthedocs.yml | 31 +++++++++++++++++++++++++------
 1 file changed, 25 insertions(+), 6 deletions(-)

diff --git a/readthedocs.yml b/readthedocs.yml
index fb2d1374bb..dde255335c 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -1,7 +1,26 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+build:
+   os: ubuntu-20.04
+   tools:
+     python: "3.7"
+   jobs:
+     # fetch the full history so that setuptools_scm can determine the
+     # correct version string for long PRs with many commits
+     post_checkout:
+       - git fetch --unshallow
+
 python:
-    version: 3
-    # only use the packages specified in setup.py
-    use_system_site_packages: false
-    pip_install: true
-    extra_requirements:
-        - doc
\ No newline at end of file
+   # only use the packages specified in setup.py
+   system_packages: false
+
+   install:
+      - method: pip
+        path: .
+        extra_requirements:
+          - doc

From 27ff16ea00b4f03b0d6f6f3bb489f211e24f6f57 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Mon, 27 Feb 2023 08:17:28 -0800
Subject: [PATCH 21/30] Update pvlib/irradiance.py

Use Sphinx math role to render k_t in docstring for min cosine zenith

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
---
 pvlib/irradiance.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index f732164ef4..6fc5d79cfb 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -2290,7 +2290,7 @@ def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
         pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
     min_cos_zenith : numeric, default 0.065
         Minimum value of cos(zenith) to allow when calculating global
-        clearness index `kt`. Equivalent to zenith = 86.273 degrees.
+        clearness index :math:`k_t`. Equivalent to zenith = 86.273 degrees.
     max_zenith : numeric, default 87
         Maximum value of zenith to allow in DNI calculation. DNI will be
         set to 0 for times with zenith values greater than `max_zenith`.

From 8e923819316e528a6a11ddb83357a1ee4eeb93ef Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 5 Mar 2023 15:19:26 -0800
Subject: [PATCH 22/30] use solar_zenith in irradiance.boland docstring since
 #1403

---
 pvlib/irradiance.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 6fc5d79cfb..0a9d61fb51 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -2283,7 +2283,7 @@ def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
     ----------
     ghi: numeric
         Global horizontal irradiance. [W/m^2]
-    zenith: numeric
+    solar_zenith: numeric
         True (not refraction-corrected) zenith angles in decimal degrees.
     datetime_or_doy : int, float, numpy.ndarray, pandas.DatetimeIndex
         Day of year or array of days of year e.g.

From 9b7f6d8b977781a3f878bdee3468932132110b58 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 5 Mar 2023 15:33:29 -0800
Subject: [PATCH 23/30] update reference to modeling diffuse fraction by J.
 Boland, available online
 https://www.researchgate.net/publication/229873508_Modelling_the_diffuse_fraction_of_global_solar_radiation_on_a_horizontal_surface

---
 pvlib/irradiance.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 0a9d61fb51..eb6286e6f8 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -2330,7 +2330,7 @@ def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
 
     # Boland equation
     df = 1.0 / (1.0 + np.exp(-5.0 + 8.6 * kt))
-    # NOTE: [1] has different coefficients, for different time intervals
+    # NOTE: [2] has different coefficients, for different time intervals
     # 15-min: df = 1 / (1 + exp(8.645 * (kt - 0.613)))
     # 1-hour: df = 1 / (1 + exp(7.997 * (kt - 0.586)))
 

From dc75be31fcae7fdd9ebc6cd88d3b87061f0846c7 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 5 Mar 2023 20:32:37 -0800
Subject: [PATCH 24/30] let user pass boland coeffs as kwargs

- add a_coeff, b_coeff to docstring params
- update equation to match Boland paper, use A, B coeffs
- update coeffs to match Boland paper too (8.6-> 8.645 and 5 -> 5.3
- add note to explain different coeffs for different time intervals
- give coeffs for 15-min & 1-hr
- update zenith to solar_zenith everywhere
- update test expected to match output with new coeffs
---
 pvlib/irradiance.py            | 30 +++++++++++++++++++++++-------
 pvlib/tests/test_irradiance.py |  8 ++++----
 2 files changed, 27 insertions(+), 11 deletions(-)

diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index eb6286e6f8..d57610beed 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -2265,7 +2265,8 @@ def erbs(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
     return data
 
 
-def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
+def boland(ghi, solar_zenith, datetime_or_doy, a_coeff=8.645, b_coeff=0.613,
+           min_cos_zenith=0.065, max_zenith=87):
     r"""
     Estimate DNI and DHI from GHI using the Boland clearness index model.
 
@@ -2276,7 +2277,7 @@ def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
 
     .. math::
 
-        \mathit{DF} = \frac{1}{1 + \exp\left(-5 + 8.6 k_t\right)}
+        \mathit{DF} = \frac{1}{1 + \exp\left(A \left(k_t - B\right)\right)}
 
 
     Parameters
@@ -2288,6 +2289,10 @@ def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
     datetime_or_doy : int, float, numpy.ndarray, pandas.DatetimeIndex
         Day of year or array of days of year e.g.
         pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
+    a_coeff : float, default 8.645
+        logistic curve fit coefficient
+    b_coeff : float, default 0.613
+        logistic curve fit coefficient
     min_cos_zenith : numeric, default 0.065
         Minimum value of cos(zenith) to allow when calculating global
         clearness index :math:`k_t`. Equivalent to zenith = 86.273 degrees.
@@ -2321,23 +2326,34 @@ def boland(ghi, zenith, datetime_or_doy, min_cos_zenith=0.065, max_zenith=87):
     dirint
     disc
     erbs
+
+    Notes
+    -----
+    Boland diffuse fraction differs from other decomposition algorithms by use
+    of a logistic function to fit the entire range of clearness index,
+    :math:`k_t`. Parameters ``a_coeff`` and ``b_coeff`` are reported in [2]_
+    for different time intervals:
+
+    * 15-minute: ``a = 8.645`` and ``b = 0.613``
+    * 1-hour:  ``a = 7.997`` and ``b = 0.586``
     """
 
     dni_extra = get_extra_radiation(datetime_or_doy)
 
-    kt = clearness_index(ghi, zenith, dni_extra, min_cos_zenith=min_cos_zenith,
-                         max_clearness_index=1)
+    kt = clearness_index(
+        ghi, solar_zenith, dni_extra, min_cos_zenith=min_cos_zenith,
+        max_clearness_index=1)
 
     # Boland equation
-    df = 1.0 / (1.0 + np.exp(-5.0 + 8.6 * kt))
+    df = 1.0 / (1.0 + np.exp(a_coeff * (kt - b_coeff)))
     # NOTE: [2] has different coefficients, for different time intervals
     # 15-min: df = 1 / (1 + exp(8.645 * (kt - 0.613)))
     # 1-hour: df = 1 / (1 + exp(7.997 * (kt - 0.586)))
 
     dhi = df * ghi
 
-    dni = (ghi - dhi) / tools.cosd(zenith)
-    bad_values = (zenith > max_zenith) | (ghi < 0) | (dni < 0)
+    dni = (ghi - dhi) / tools.cosd(solar_zenith)
+    bad_values = (solar_zenith > max_zenith) | (ghi < 0) | (dni < 0)
     dni = np.where(bad_values, 0, dni)
     # ensure that closure relationship remains valid
     dhi = np.where(bad_values, ghi, dhi)
diff --git a/pvlib/tests/test_irradiance.py b/pvlib/tests/test_irradiance.py
index 732d574767..cffdd23e40 100644
--- a/pvlib/tests/test_irradiance.py
+++ b/pvlib/tests/test_irradiance.py
@@ -810,14 +810,14 @@ def test_boland():
     zenith = pd.Series([120, 85, 10, 10], index=index)
     expected = pd.DataFrame(np.array(
         [[0.0,        0.0,        0.0],
-         [103.735879, 40.958822,  0.405724],
-         [776.006568, 235.782716, 0.718133],
-         [845.794317, 167.055199, 0.768214]]),
+         [81.9448546, 42.8580353, 0.405723511],
+         [723.764990, 287.230626, 0.718132729],
+         [805.020419, 207.209650, 0.768214312]]),
         columns=['dni', 'dhi', 'kt'], index=index)
 
     out = irradiance.boland(ghi, zenith, index)
 
-    assert_frame_equal(np.round(out, 0), np.round(expected, 0))
+    assert np.allclose(out, expected)
 
 
 def test_erbs_min_cos_zenith_max_zenith():

From bee5bfd41b91781736dd11f7b43c1ec946cdf2de Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 5 Mar 2023 21:22:49 -0800
Subject: [PATCH 25/30] move plot diffuse fraction example to subdir

- create irradiance-decomposition subdir in gallery, add readme.rst
- use complete_irradiance() to get DHI using closure eqn's
---
 .../irradiance-decomposition/README.rst        |  3 +++
 .../plot_diffuse_fraction.py                   | 18 ++++++++++--------
 2 files changed, 13 insertions(+), 8 deletions(-)
 create mode 100644 docs/examples/irradiance-decomposition/README.rst
 rename docs/examples/{ => irradiance-decomposition}/plot_diffuse_fraction.py (93%)

diff --git a/docs/examples/irradiance-decomposition/README.rst b/docs/examples/irradiance-decomposition/README.rst
new file mode 100644
index 0000000000..95f18c9330
--- /dev/null
+++ b/docs/examples/irradiance-decomposition/README.rst
@@ -0,0 +1,3 @@
+Irradiance Decomposition
+------------------------
+
diff --git a/docs/examples/plot_diffuse_fraction.py b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
similarity index 93%
rename from docs/examples/plot_diffuse_fraction.py
rename to docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
index 94786c999a..6a090f1306 100644
--- a/docs/examples/plot_diffuse_fraction.py
+++ b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
@@ -55,10 +55,12 @@
 
 out_disc = irradiance.disc(
     greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100)
+# use "complete sum" AKA "closure" equations: DHI = GHI - DNI * cos(zenith)
+df = irradiance.complete_irradiance(
+    solar_zenith=solpos.apparent_zenith, ghi=greensboro.GHI, dni=out_disc.dni,
+    dhi=None)
 out_disc = out_disc.rename(columns={'dni': 'dni_disc'})
-out_disc['dhi_disc'] = (
-    greensboro.GHI
-    - out_disc.dni_disc*np.cos(np.radians(solpos.apparent_zenith)))
+out_disc['dhi_disc'] = df.dhi
 
 # %%
 # DIRINT
@@ -70,12 +72,12 @@
 dni_dirint = irradiance.dirint(
     greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100,
     temp_dew=greensboro.DewPoint)
-dhi_dirint = (
-    greensboro.GHI
-    - dni_dirint*np.cos(np.radians(solpos.apparent_zenith)))
+# use "complete sum" AKA "closure" equations: DHI = GHI - DNI * cos(zenith)
+df = irradiance.complete_irradiance(
+    solar_zenith=solpos.apparent_zenith, ghi=greensboro.GHI, dni=dni_dirint,
+    dhi=None)
 out_dirint = pd.DataFrame(
-    {'dni_dirint': dni_dirint, 'dhi_dirint': dhi_dirint},
-    index=greensboro.index)
+    {'dni_dirint': dni_dirint, 'dhi_dirint': df.dhi}, index=greensboro.index)
 
 # %%
 # Erbs

From 52ca13c1a9f06617adf1c7e7e17cc0a993777e20 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 5 Mar 2023 21:55:55 -0800
Subject: [PATCH 26/30] shift solar position calc to hour center

- in plot diffuse fraction examples
- add note to explain b/c TMY timestampe at hour end
- reset index to align with TMY3
- fix Jan, July, Apr timestamp names didn't match actual times
---
 .../plot_diffuse_fraction.py                  | 29 ++++++++++---------
 1 file changed, 16 insertions(+), 13 deletions(-)

diff --git a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
index 6a090f1306..1a143d8bf0 100644
--- a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
+++ b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
@@ -32,11 +32,14 @@
 
 # Many of the diffuse fraction estimation methods require the "true" zenith, so
 # we calculate the solar positions for the 1990 at Greensboro, NC.
+# NOTE: TMY3 files timestamps indicate the end of the hour, so shift indices
+# back 30-minutes to calculate solar position at center of the interval
 solpos = get_solarposition(
-    greensboro.index, latitude=metadata['latitude'],
+    greensboro.index.shift(freq="-30T"), latitude=metadata['latitude'],
     longitude=metadata['longitude'], altitude=metadata['altitude'],
     pressure=greensboro.Pressure*100,  # convert from millibar to Pa
     temperature=greensboro.DryBulb)
+solpos.index = greensboro.index  # reset index to end of the hour
 
 # %%
 # pvlib Decomposition Functions
@@ -135,16 +138,16 @@
 # ++++++
 # Finally, let's plot them for a few winter days and compare
 
-JAN6AM, JAN6PM = '1990-01-04 00:00:00-05:00', '1990-01-07 23:59:59-05:00'
+JAN_AM, JAN_PM = '1990-01-04 00:00:00-05:00', '1990-01-07 23:59:59-05:00'
 f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
-dni[JAN6AM:JAN6PM].plot(ax=ax[0])
+dni[JAN_AM:JAN_PM].plot(ax=ax[0])
 ax[0].grid(which="both")
 ax[0].set_ylabel('DNI $[W/m^2]$')
 ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
-dhi[JAN6AM:JAN6PM].plot(ax=ax[1])
+dhi[JAN_AM:JAN_PM].plot(ax=ax[1])
 ax[1].grid(which="both")
 ax[1].set_ylabel('DHI $[W/m^2]$')
-ghi_kt[JAN6AM:JAN6PM].plot(ax=ax[2])
+ghi_kt[JAN_AM:JAN_PM].plot(ax=ax[2])
 ax[2].grid(which='both')
 ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
 f.tight_layout()
@@ -154,16 +157,16 @@
 # ++++++
 # And a few spring days ...
 
-APR6AM, APR6PM = '1990-04-04 00:00:00-05:00', '1990-04-07 23:59:59-05:00'
+APR_AM, APR_PM = '1990-04-04 00:00:00-05:00', '1990-04-07 23:59:59-05:00'
 f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
-dni[APR6AM:APR6PM].plot(ax=ax[0])
+dni[APR_AM:APR_PM].plot(ax=ax[0])
 ax[0].grid(which="both")
 ax[0].set_ylabel('DNI $[W/m^2]$')
 ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
-dhi[APR6AM:APR6PM].plot(ax=ax[1])
+dhi[APR_AM:APR_PM].plot(ax=ax[1])
 ax[1].grid(which="both")
 ax[1].set_ylabel('DHI $[W/m^2]$')
-ghi_kt[APR6AM:APR6PM].plot(ax=ax[2])
+ghi_kt[APR_AM:APR_PM].plot(ax=ax[2])
 ax[2].grid(which='both')
 ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
 f.tight_layout()
@@ -173,16 +176,16 @@
 # ++++++
 # And few summer days to finish off the seasons.
 
-JUL6AM, JUL6PM = '1990-07-04 00:00:00-05:00', '1990-07-07 23:59:59-05:00'
+JUL_AM, JUL_PM = '1990-07-04 00:00:00-05:00', '1990-07-07 23:59:59-05:00'
 f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
-dni[JUL6AM:JUL6PM].plot(ax=ax[0])
+dni[JUL_AM:JUL_PM].plot(ax=ax[0])
 ax[0].grid(which="both")
 ax[0].set_ylabel('DNI $[W/m^2]$')
 ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
-dhi[JUL6AM:JUL6PM].plot(ax=ax[1])
+dhi[JUL_AM:JUL_PM].plot(ax=ax[1])
 ax[1].grid(which="both")
 ax[1].set_ylabel('DHI $[W/m^2]$')
-ghi_kt[JUL6AM:JUL6PM].plot(ax=ax[2])
+ghi_kt[JUL_AM:JUL_PM].plot(ax=ax[2])
 ax[2].grid(which='both')
 ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
 f.tight_layout()

From 1a3e280f2bbed053a1f1bc421778c902f95d8ae6 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 5 Mar 2023 22:09:48 -0800
Subject: [PATCH 27/30] don't need numpy in plot diffuse fraction examples

---
 docs/examples/irradiance-decomposition/plot_diffuse_fraction.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
index 1a143d8bf0..3b381cf046 100644
--- a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
+++ b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
@@ -16,7 +16,6 @@
 
 import pathlib
 from matplotlib import pyplot as plt
-import numpy as np
 import pandas as pd
 from pvlib.iotools import read_tmy3
 from pvlib.solarposition import get_solarposition

From 20e9a53b77287518b7a30ec7892f0689a014af52 Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 5 Mar 2023 22:23:02 -0800
Subject: [PATCH 28/30] minor edits in plot diffuse fraction example

- chnage df to be specific to disc or dirint
- change date names to Jan-04, Jan-07, etc. instead of _AM, _PM etc.
---
 .../plot_diffuse_fraction.py                  | 33 ++++++++++---------
 1 file changed, 17 insertions(+), 16 deletions(-)

diff --git a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
index 3b381cf046..dbbfd5a15d 100644
--- a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
+++ b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
@@ -58,11 +58,11 @@
 out_disc = irradiance.disc(
     greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100)
 # use "complete sum" AKA "closure" equations: DHI = GHI - DNI * cos(zenith)
-df = irradiance.complete_irradiance(
+df_disc = irradiance.complete_irradiance(
     solar_zenith=solpos.apparent_zenith, ghi=greensboro.GHI, dni=out_disc.dni,
     dhi=None)
 out_disc = out_disc.rename(columns={'dni': 'dni_disc'})
-out_disc['dhi_disc'] = df.dhi
+out_disc['dhi_disc'] = df_disc.dhi
 
 # %%
 # DIRINT
@@ -75,11 +75,12 @@
     greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100,
     temp_dew=greensboro.DewPoint)
 # use "complete sum" AKA "closure" equations: DHI = GHI - DNI * cos(zenith)
-df = irradiance.complete_irradiance(
+df_dirint = irradiance.complete_irradiance(
     solar_zenith=solpos.apparent_zenith, ghi=greensboro.GHI, dni=dni_dirint,
     dhi=None)
 out_dirint = pd.DataFrame(
-    {'dni_dirint': dni_dirint, 'dhi_dirint': df.dhi}, index=greensboro.index)
+    {'dni_dirint': dni_dirint, 'dhi_dirint': df_dirint.dhi},
+    index=greensboro.index)
 
 # %%
 # Erbs
@@ -137,16 +138,16 @@
 # ++++++
 # Finally, let's plot them for a few winter days and compare
 
-JAN_AM, JAN_PM = '1990-01-04 00:00:00-05:00', '1990-01-07 23:59:59-05:00'
+JAN04, JAN07 = '1990-01-04 00:00:00-05:00', '1990-01-07 23:59:59-05:00'
 f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
-dni[JAN_AM:JAN_PM].plot(ax=ax[0])
+dni[JAN04:JAN07].plot(ax=ax[0])
 ax[0].grid(which="both")
 ax[0].set_ylabel('DNI $[W/m^2]$')
 ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
-dhi[JAN_AM:JAN_PM].plot(ax=ax[1])
+dhi[JAN04:JAN07].plot(ax=ax[1])
 ax[1].grid(which="both")
 ax[1].set_ylabel('DHI $[W/m^2]$')
-ghi_kt[JAN_AM:JAN_PM].plot(ax=ax[2])
+ghi_kt[JAN04:JAN07].plot(ax=ax[2])
 ax[2].grid(which='both')
 ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
 f.tight_layout()
@@ -156,16 +157,16 @@
 # ++++++
 # And a few spring days ...
 
-APR_AM, APR_PM = '1990-04-04 00:00:00-05:00', '1990-04-07 23:59:59-05:00'
+APR04, APR07 = '1990-04-04 00:00:00-05:00', '1990-04-07 23:59:59-05:00'
 f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
-dni[APR_AM:APR_PM].plot(ax=ax[0])
+dni[APR04:APR07].plot(ax=ax[0])
 ax[0].grid(which="both")
 ax[0].set_ylabel('DNI $[W/m^2]$')
 ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
-dhi[APR_AM:APR_PM].plot(ax=ax[1])
+dhi[APR04:APR07].plot(ax=ax[1])
 ax[1].grid(which="both")
 ax[1].set_ylabel('DHI $[W/m^2]$')
-ghi_kt[APR_AM:APR_PM].plot(ax=ax[2])
+ghi_kt[APR04:APR07].plot(ax=ax[2])
 ax[2].grid(which='both')
 ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
 f.tight_layout()
@@ -175,16 +176,16 @@
 # ++++++
 # And few summer days to finish off the seasons.
 
-JUL_AM, JUL_PM = '1990-07-04 00:00:00-05:00', '1990-07-07 23:59:59-05:00'
+JUL04, JUL07 = '1990-07-04 00:00:00-05:00', '1990-07-07 23:59:59-05:00'
 f, ax = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
-dni[JUL_AM:JUL_PM].plot(ax=ax[0])
+dni[JUL04:JUL07].plot(ax=ax[0])
 ax[0].grid(which="both")
 ax[0].set_ylabel('DNI $[W/m^2]$')
 ax[0].set_title('Comparison of Diffuse Fraction Estimation Methods')
-dhi[JUL_AM:JUL_PM].plot(ax=ax[1])
+dhi[JUL04:JUL07].plot(ax=ax[1])
 ax[1].grid(which="both")
 ax[1].set_ylabel('DHI $[W/m^2]$')
-ghi_kt[JUL_AM:JUL_PM].plot(ax=ax[2])
+ghi_kt[JUL04:JUL07].plot(ax=ax[2])
 ax[2].grid(which='both')
 ax[2].set_ylabel(r'$\frac{GHI}{E0}, k_t$')
 f.tight_layout()

From 258ffd2b068c1b5c72508cfc5c7371b707b1aa5e Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 12 Mar 2023 16:27:18 -0700
Subject: [PATCH 29/30] Apply suggestions from code review by Adam

* use lower case coeffs in equation
* periods in docstrings
* some wordsmithing
* replace parameter types for datetime_or_doy with numeric
* consistent def for kt is ratio of GHI to horizontal extraterrestrial-irradiance

Co-authored-by: Adam R. Jensen <39184289+AdamRJensen@users.noreply.github.com>
---
 .../plot_diffuse_fraction.py                         | 12 ++++++------
 pvlib/irradiance.py                                  | 12 ++++++------
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
index dbbfd5a15d..df1a841c11 100644
--- a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
+++ b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
@@ -74,7 +74,7 @@
 dni_dirint = irradiance.dirint(
     greensboro.GHI, solpos.zenith, greensboro.index, greensboro.Pressure*100,
     temp_dew=greensboro.DewPoint)
-# use "complete sum" AKA "closure" equations: DHI = GHI - DNI * cos(zenith)
+# use "complete sum" AKA "closure" equation: DHI = GHI - DNI * cos(zenith)
 df_dirint = irradiance.complete_irradiance(
     solar_zenith=solpos.apparent_zenith, ghi=greensboro.GHI, dni=dni_dirint,
     dhi=None)
@@ -201,11 +201,11 @@
 # reproduce the TMY3 values. We refer those interested to the `TMY3`_ and
 # `NSRDB`_ user manuals.
 #
-# The Erbs and Boland models are correlations with only kt, which is derived
-# from the horizontal component of the extra-terrestrial irradiance. At low sun
-# elevation (zenith near 90 degrees), especially near sunset, kt can explode
-# because the denominator (extra-terrestrial irradiance) approaches zero. In
-# pvlib this behavior is moderated by ``min_cos_zenith`` and
+# The Erbs and Boland models are correlations only based on the clearness index kt,
+# which is the ratio of GHI to the the horizontal component of the extra-terrestrial
+# irradiance. At low sun elevation (zenith near 90 degrees), especially near sunset,
+# kt can explode because the denominator (extra-terrestrial irradiance) approaches
+# zero. In pvlib this behavior is moderated by ``min_cos_zenith`` and
 # ``max_clearness_index`` which each have reasonable defaults. Even so, near
 # sunset there are still spikes in kt and DNI from Erbs and Boland for Jan. 5th
 # & 7th, April 4th, 5th, & 7th, and July 6th & 7th.
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index d57610beed..79ce3c3904 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -2272,12 +2272,12 @@ def boland(ghi, solar_zenith, datetime_or_doy, a_coeff=8.645, b_coeff=0.613,
 
     The Boland model [1]_, [2]_ estimates the diffuse fraction, DF, from global
     horizontal irradiance, GHI, through an empirical relationship between DF
-    and the clearness index, :math:`k_t`, the ratio of GHI to extraterrestrial
-    irradiance.
+    and the clearness index, :math:`k_t`, the ratio of GHI to horizontal
+    extraterrestrial irradiance.
 
     .. math::
 
-        \mathit{DF} = \frac{1}{1 + \exp\left(A \left(k_t - B\right)\right)}
+        \mathit{DF} = \frac{1}{1 + \exp\left(a \left(k_t - b\right)\right)}
 
 
     Parameters
@@ -2286,13 +2286,13 @@ def boland(ghi, solar_zenith, datetime_or_doy, a_coeff=8.645, b_coeff=0.613,
         Global horizontal irradiance. [W/m^2]
     solar_zenith: numeric
         True (not refraction-corrected) zenith angles in decimal degrees.
-    datetime_or_doy : int, float, numpy.ndarray, pandas.DatetimeIndex
+    datetime_or_doy : numeric, pandas.DatetimeIndex
         Day of year or array of days of year e.g.
         pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
     a_coeff : float, default 8.645
-        logistic curve fit coefficient
+        Logistic curve fit coefficient.
     b_coeff : float, default 0.613
-        logistic curve fit coefficient
+        Logistic curve fit coefficient.
     min_cos_zenith : numeric, default 0.065
         Minimum value of cos(zenith) to allow when calculating global
         clearness index :math:`k_t`. Equivalent to zenith = 86.273 degrees.

From 1f6b9397b162e1960c07a46d770c3800457bcbbc Mon Sep 17 00:00:00 2001
From: Mark Mikofski <bwana.marko@yahoo.com>
Date: Sun, 12 Mar 2023 16:31:21 -0700
Subject: [PATCH 30/30] fix stickler - long lines in plot SF example

---
 .../plot_diffuse_fraction.py                    | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
index df1a841c11..dbd0406cf6 100644
--- a/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
+++ b/docs/examples/irradiance-decomposition/plot_diffuse_fraction.py
@@ -201,14 +201,15 @@
 # reproduce the TMY3 values. We refer those interested to the `TMY3`_ and
 # `NSRDB`_ user manuals.
 #
-# The Erbs and Boland models are correlations only based on the clearness index kt,
-# which is the ratio of GHI to the the horizontal component of the extra-terrestrial
-# irradiance. At low sun elevation (zenith near 90 degrees), especially near sunset,
-# kt can explode because the denominator (extra-terrestrial irradiance) approaches
-# zero. In pvlib this behavior is moderated by ``min_cos_zenith`` and
-# ``max_clearness_index`` which each have reasonable defaults. Even so, near
-# sunset there are still spikes in kt and DNI from Erbs and Boland for Jan. 5th
-# & 7th, April 4th, 5th, & 7th, and July 6th & 7th.
+# The Erbs and Boland models are correlations only based on the clearness index
+# kt, which is the ratio of GHI to the the horizontal component of the
+# extra-terrestrial irradiance. At low sun elevation (zenith near 90 degrees),
+# especially near sunset, kt can explode because the denominator
+# (extra-terrestrial irradiance) approaches zero. In pvlib this behavior is
+# moderated by ``min_cos_zenith`` and ``max_clearness_index`` which each have
+# reasonable defaults. Even so, near sunset there are still spikes in kt and
+# DNI from Erbs and Boland for Jan. 5th & 7th, April 4th, 5th, & 7th, and July
+# 6th & 7th.
 #
 # By contrast, the DISC and DIRINT methods estimate DNI first by means of
 # correlations, which include additional variables such as airmass. These