use horner's method in disc

docs/sphinx/source/whatsnew/v0.9.5.rst

@@ -26,13 +26,14 @@ Enhancements
   differentiable and bounded between zero and one. (:pull:`1179`)
 * Add ``model='gueymard2003'``, the airmass model used for REST and REST2,
   to :py:func:`~pvlib.atmosphere.get_relative_airmass`. (:pull:`1655`)
-
-
 * Added an optional ``model`` parameter to
   :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance` and
   :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa`
   to enable use of the hay-davies sky diffuse irradiance model
   instead of the default isotropic model. (:pull:`1668`)
+* Use `Horner's Method <https://en.wikipedia.org/wiki/Horner%27s_method>`_
+  to evaluate polynomials in :py:func:`~pvlib.irradiance.disc`, may
+  decrease runtime by 20%. (:issue:`1180`, :pull:`1183`)
 
 Bug fixes
 ~~~~~~~~~

pvlib/irradiance.py

@@ -1481,24 +1481,24 @@ def _disc_kn(clearness_index, airmass, max_airmass=12):
 
     am = np.minimum(am, max_airmass)  # GH 450
 
-    # powers of kt will be used repeatedly, so compute only once
-    kt2 = kt * kt  # about the same as kt ** 2
-    kt3 = kt2 * kt  # 5-10x faster than kt ** 3
-
-    bools = (kt <= 0.6)
-    a = np.where(bools,
-                 0.512 - 1.56*kt + 2.286*kt2 - 2.222*kt3,
-                 -5.743 + 21.77*kt - 27.49*kt2 + 11.56*kt3)
-    b = np.where(bools,
-                 0.37 + 0.962*kt,
-                 41.4 - 118.5*kt + 66.05*kt2 + 31.9*kt3)
-    c = np.where(bools,
-                 -0.28 + 0.932*kt - 2.048*kt2,
-                 -47.01 + 184.2*kt - 222.0*kt2 + 73.81*kt3)
+    is_cloudy = (kt <= 0.6)
+    # Use Horner's method to compute polynomials efficiently
+    a = np.where(
+        is_cloudy,
+        0.512 + kt*(-1.56 + kt*(2.286 - 2.222*kt)),
+        -5.743 + kt*(21.77 + kt*(-27.49 + 11.56*kt)))
+    b = np.where(
+        is_cloudy,
+        0.37 + 0.962*kt,
+        41.4 + kt*(-118.5 + kt*(66.05 + 31.9*kt)))
+    c = np.where(
+        is_cloudy,
+        -0.28 + kt*(0.932 - 2.048*kt),
+        -47.01 + kt*(184.2 + kt*(-222.0 + 73.81*kt)))
 
     delta_kn = a + b * np.exp(c*am)
 
-    Knc = 0.866 - 0.122*am + 0.0121*am**2 - 0.000653*am**3 + 1.4e-05*am**4
+    Knc = 0.866 + am*(-0.122 + am*(0.0121 + am*(-0.000653 + 1.4e-05*am)))
     Kn = Knc - delta_kn
     return Kn, am