Shadowband correction factor function and test

docs/source/documentation.rst

@@ -11,6 +11,7 @@ Code documentation
    plotting.irradiance_colormap_and_norm
    plotting.plot_google_maps
    plotting.plot_intraday_heatmap
+   instrument.shadowband_correction_factor
    plotting.plot_shading_heatmap
    horizon.get_horizon_mines
    quality.bsrn_limits

pyproject.toml

@@ -47,8 +47,8 @@ dynamic = ["version"]
 [project.optional-dependencies]
 test = ["pytest>=7", "pytest-cov", "packaging"]
 doc = [
-    "sphinx==9.1.0",
-    "myst-nb==1.4.0",
+    "sphinx==9.0.4",
+    "myst-nb==1.3.0",
     "sphinx-book-theme==1.1.4",
     "sphinx-copybutton==0.5.2",
     "pvlib==0.14.0",

src/solarpy/__init__.py

@@ -11,5 +11,6 @@
     quality,
     horizon,
     example,
+    instrument,
     iotools,
 )

src/solarpy/instrument/__init__.py

@@ -0,0 +1 @@
+from solarpy.instrument.shadowband import shadowband_correction_factor  # noqa: F401

src/solarpy/instrument/shadowband.py

@@ -0,0 +1,128 @@
+"""
+shadowband_correction.py
+------------------------
+Shadowband diffuse-irradiance correction factor for the Kipp & Zonen CM 121
+shadow-ring pyranometer accessory.
+
+Reference
+---------
+Kipp & Zonen *CM 121 B/C Shadow Ring Instruction Manual*, Appendix §6.1.
+https://www.kippzonen.com/Download/46/CM121-B-C-Shadow-Ring-Manual
+"""
+
+import numpy as np
+import pvlib as pv
+import pandas as pd
+import datetime
+
+
+def shadowband_correction_factor(
+    input_date,
+    latitude: float,
+    V: float = 0.185,
+):
+    """Calculate the daily shadowband correction factor for a Kipp & Zonen CM 121.
+
+    A shadow ring blocks a strip of the sky dome to prevent direct beam
+    radiation from reaching the pyranometer, so that only diffuse radiation
+    is measured.  Because the ring also blocks some of the diffuse sky, a
+    correction factor C > 1 must be applied to recover the true diffuse
+    irradiance::
+
+        G_d_true = C * G_d_measured
+
+    The factor depends on the ring geometry (angular half-width *V*), the
+    site latitude, and the solar declination on the measurement day.  It is
+    derived analytically by integrating the fraction of the isotropic sky
+    hemisphere obscured by the ring over the daily arc of the sun.
+
+    Parameters
+    ----------
+    input_date : scalar or array-like
+        Date(s) for which to compute the factor.  Accepts anything that
+        ``pandas.to_datetime`` understands: a ``pd.Timestamp``, a
+        ``datetime.date`` / ``datetime.datetime``, an ISO-8601 string, a
+        list/array of any of the above, or a ``pd.DatetimeIndex``.
+        Sub-daily time information is ignored; only the calendar date matters.
+    latitude : float
+        Observer latitude in decimal degrees.  Positive = North, negative = South.
+        Valid range: −90 to +90.
+    V : float, optional
+        Angular width of the shadow band in radians.  This is arc subtended by
+        the ring as seen from the sensor.  Default is 0.185 rad (~10.6°),
+        which matches the standard Kipp & Zonen CM 121 ring geometry.
+
+    Returns
+    -------
+    float or array-like
+        Dimensionless correction factor C ≥ 1.
+
+    Notes
+    -----
+    The correction factor is computed as:
+
+    .. math::
+
+        S = \\frac{2V \\cos\\delta}{\\pi}
+            \\bigl(U_0 \\sin\\phi\\sin\\delta
+                  + \\sin U_0 \\cos\\phi\\cos\\delta\\bigr)
+
+        C = \\frac{1}{1 - S}
+
+    where
+
+    * :math:`\\delta` – solar declination (Spencer 1971 approximation)
+    * :math:`\\phi`   – site latitude
+    * :math:`U_0`    – sunrise/sunset hour angle,
+      :math:`\\arccos(-\\tan\\phi\\,\\tan\\delta)`, clamped to [−1, 1] to
+      handle polar day/night conditions
+    """
+    # ------------------------------------------------------------------
+    # 1. Normalise input to pandas datetime so scalar and array paths
+    #    share a single code route below.
+    # ------------------------------------------------------------------
+    dates = pd.to_datetime(input_date)
+
+    # ------------------------------------------------------------------
+    # 2. Solar declination δ for each day of year.
+    #    Spencer (1971) approximation, returned in radians by pvlib.
+    # ------------------------------------------------------------------
+    day_of_year = dates.dayofyear
+    delta = pv.solarposition.declination_spencer71(day_of_year)  # radians
+
+    # ------------------------------------------------------------------
+    # 3. Site latitude φ in radians.
+    # ------------------------------------------------------------------
+    phi = np.deg2rad(latitude)
+
+    # ------------------------------------------------------------------
+    # 4. Sunrise/sunset hour angle U₀.
+    #    cos U₀ = −tan φ · tan δ.
+    #    Clamped to [−1, 1] so arccos is valid at the poles (midnight sun
+    #    / polar night), where the tangent product can exceed ±1.
+    # ------------------------------------------------------------------
+    cos_U0 = np.clip(-np.tan(phi) * np.tan(delta), -1.0, 1.0)
+    U0 = np.arccos(cos_U0)
+
+    # ------------------------------------------------------------------
+    # 5. Fraction S of the isotropic sky dome obscured by the ring.
+    #    Derived by integrating the ring's shadow across the daily solar arc
+    #    (Kipp & Zonen manual, Appendix §6.1).
+    # ------------------------------------------------------------------
+    S = (2 * V * np.cos(delta) / np.pi) * (
+        U0 * np.sin(phi) * np.sin(delta)
+        + np.sin(U0) * np.cos(phi) * np.cos(delta)
+    )
+
+    # ------------------------------------------------------------------
+    # 6. Correction factor C = 1 / (1 − S).
+    #    S is always < 1 for physically reasonable inputs, so C > 1.
+    # ------------------------------------------------------------------
+    C = 1.0 / (1.0 - S)
+
+    # ------------------------------------------------------------------
+    # 7. Return a plain float for scalar input, Series for array input.
+    # ------------------------------------------------------------------
+    if np.isscalar(input_date) or isinstance(input_date, (pd.Timestamp, datetime.date, datetime.datetime)):
+        return float(C)
+    return pd.Series(C, index=dates)

tests/test_instrument.py

@@ -0,0 +1,58 @@
+import solarpy
+import pandas as pd
+import datetime
+import pytest
+
+
+def test_shadowband_correction():
+    # Northern hemisphere, spring (near-zero declination)
+    assert solarpy.instrument.shadowband_correction_factor("2026-03-17", latitude=55.29007) == 1.0668908588370676
+    # Northern hemisphere, spring - pd.Timestamp
+    assert solarpy.instrument.shadowband_correction_factor(pd.Timestamp("2026-03-17"), latitude=55.29007) == 1.0668908588370676
+    # pd.Timestamp with time component - sub-daily info must be ignored
+    assert solarpy.instrument.shadowband_correction_factor(pd.Timestamp("2026-03-17 14:30:01"), latitude=55.29007) == pytest.approx(1.0668908588370676)
+    # datetime.date scalar
+    assert solarpy.instrument.shadowband_correction_factor(datetime.date(2026, 3, 17), latitude=55.29007) == pytest.approx(1.0668908588370676)
+    # datetime.datetime scalar
+    assert solarpy.instrument.shadowband_correction_factor(datetime.datetime(2026, 3, 17, 14, 30), latitude=55.29007) == pytest.approx(1.0668908588370676)
+
+    # Northern hemisphere, summer
+    assert solarpy.instrument.shadowband_correction_factor("2026-06-21", latitude=55.29007) == pytest.approx(1.1408347415382885)
+    # Northern hemisphere, autumn
+    assert solarpy.instrument.shadowband_correction_factor("2026-10-21", latitude=55.29007) == pytest.approx(1.0418060380205798)
+    # Northern hemisphere, winter
+    assert solarpy.instrument.shadowband_correction_factor("2026-12-21", latitude=55.29007) == pytest.approx(1.0126112472695903)
+
+    # Southern hemisphere
+    assert solarpy.instrument.shadowband_correction_factor("2026-10-21", latitude=-33.87) == pytest.approx(1.1282175055807726)
+
+def test_shadowband_correction_series():
+    # Single value series
+    index = pd.to_datetime(["2026-03-17"])
+    pd.testing.assert_series_equal(
+        solarpy.instrument.shadowband_correction_factor(index, latitude=55.29007),
+        pd.Series(1.0668908588370676, index=index))
+
+    # Multiple value series
+    index = pd.DatetimeIndex(["2026-03-17", "2026-10-21"])
+    pd.testing.assert_series_equal(
+    solarpy.instrument.shadowband_correction_factor(index, latitude=55.29007),
+    pd.Series([1.0668908588370676, 1.0418060380205798], index=index))
+
+def test_shadowband_correction_always_above_one():
+    # Test if shadowband correction factor C > 1 
+    dates = pd.DatetimeIndex(["2026-03-20", "2026-06-21", "2026-09-22", "2026-12-21"])
+    result = solarpy.instrument.shadowband_correction_factor(dates, latitude=55.29007)
+    assert (result > 1.0).all()
+
+
+def test_shadowband_correction_v_zero():
+    # With no ring there is nothing to correct for
+    assert solarpy.instrument.shadowband_correction_factor("2026-06-21", latitude=45.0, V=0.0) == pytest.approx(1.0)
+
+def test_shadowband_correction_polar():
+    # Test that the clip works and it follows the table provided by Kipp & Zonen
+    # At the pole in winter the sun never rises — ring blocks nothing, C = 1
+    assert solarpy.instrument.shadowband_correction_factor("2026-01-01", latitude=90.0) == pytest.approx(1.0)
+    # At the pole in summer the sun skims the horizon — C rises above 1
+    assert solarpy.instrument.shadowband_correction_factor("2026-06-21", latitude=90.0) > 1.0