add max_airmass=12 kwarg to disc (#626)

* add max_airmass=12 kwarg to disc

* more tests, clean up

* fix bad whatsnew merge

* more merge fixing

* move limit so it helps in gti_dirint. reword whatsnew

Author: wholmgren

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

@@ -21,6 +21,14 @@ API Changes
 * Added keyword argument ``horizon`` to :func:`~pvlib.solarposition.pyephem`
   and :func:`~pvlib.solarposition.calc_time` with default value ``'+0:00'``.
   (:issue:`588`)
+* Add `max_airmass` keyword argument to :py:func:`pvlib.irradiance.disc`.
+  Default value (`max_airmass=12`) is consistent with polynomial fit in
+  original paper describing the model. This change may result in different
+  output of functions that use the `disc` *Kn* calculation for times when
+  input zenith angles approach 90 degrees. This includes
+  :py:func:`pvlib.irradiance.dirint` and :py:func:`pvlib.irradiance.dirindex`
+  when `min_cos_zenith` and `max_zenith` kwargs are used, as well as
+  :py:func:`pvlib.irradiance.gti_dirint`. (:issue:`450`)
 * Changed key names for `components` returned from
   :py:func:`pvlib.clearsky.detect_clearsky`. (:issue:`596`)
 * Changed function name from `pvlib.solarposition.get_rise_set_transit`

pvlib/irradiance.py

@@ -1295,7 +1295,7 @@ def _kt_kt_prime_factor(airmass):
 
 
 def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
-         min_cos_zenith=0.065, max_zenith=87):
+         min_cos_zenith=0.065, max_zenith=87, max_airmass=12):
     """
     Estimate Direct Normal Irradiance from Global Horizontal Irradiance
     using the DISC model.
@@ -1338,6 +1338,11 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
         Maximum value of zenith to allow in DNI calculation. DNI will be
         set to 0 for times with zenith values greater than `max_zenith`.
 
+    max_airmass : numeric, default 12
+        Maximum value of the airmass to allow in Kn calculation.
+        Default value (12) comes from range over which Kn was fit
+        to airmass in the original paper.
+
     Returns
     -------
     output : OrderedDict or DataFrame
@@ -1376,7 +1381,7 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
     if pressure is not None:
         am = atmosphere.get_absolute_airmass(am, pressure)
 
-    Kn = _disc_kn(kt, am)
+    Kn, am = _disc_kn(kt, am, max_airmass=max_airmass)
     dni = Kn * I0
 
     bad_values = (solar_zenith > max_zenith) | (ghi < 0) | (dni < 0)
@@ -1393,16 +1398,29 @@ def disc(ghi, solar_zenith, datetime_or_doy, pressure=101325,
     return output
 
 
-def _disc_kn(clearness_index, airmass):
+def _disc_kn(clearness_index, airmass, max_airmass=12):
     """
     Calculate Kn for `disc`
+
+    Parameters
+    ----------
+    clearness_index : numeric
+    airmass : numeric
+    max_airmass : float
+        airmass > max_airmass is set to max_airmass before being used
+        in calculating Kn.
+
+    Returns
+    -------
+    Kn : numeric
+    am : numeric
+        airmass used in the calculation of Kn. am <= max_airmass.
     """
     # short names for equations
     kt = clearness_index
     am = airmass
 
-    # consider adding
-    # am = np.maximum(am, 12)  # GH 450
+    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
@@ -1423,7 +1441,7 @@ def _disc_kn(clearness_index, airmass):
 
     Knc = 0.866 - 0.122*am + 0.0121*am**2 - 0.000653*am**3 + 1.4e-05*am**4
     Kn = Knc - delta_kn
-    return Kn
+    return Kn, am
 
 
 def dirint(ghi, solar_zenith, times, pressure=101325., use_delta_kt_prime=True,
@@ -1934,7 +1952,7 @@ def _gti_dirint_lt_90(poa_global, aoi, aoi_lt_90, solar_zenith, solar_azimuth,
 
         # calculate kt and DNI from GTI
         kt = clearness_index(poa_global_i, aoi, I0)  # kt from Marion eqn 2
-        disc_dni = np.maximum(_disc_kn(kt, airmass) * I0, 0)
+        disc_dni = np.maximum(_disc_kn(kt, airmass)[0] * I0, 0)
         kt_prime = clearness_index_zenith_independent(kt, airmass)
         # dirint DNI in Marion eqn 3
         dni = _dirint_from_dni_ktprime(disc_dni, kt_prime, solar_zenith,
@@ -2017,7 +2035,7 @@ def _gti_dirint_gte_90(poa_global, aoi, solar_zenith, solar_azimuth,
     airmass = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966')
     airmass = atmosphere.get_absolute_airmass(airmass, pressure)
     kt = kt_prime_gte_90 * _kt_kt_prime_factor(airmass)
-    disc_dni = np.maximum(_disc_kn(kt, airmass) * I0, 0)
+    disc_dni = np.maximum(_disc_kn(kt, airmass)[0] * I0, 0)
 
     dni_gte_90 = _dirint_from_dni_ktprime(disc_dni, kt_prime, solar_zenith,
                                           False, temp_dew)

pvlib/test/test_irradiance.py

@@ -404,28 +404,57 @@ def test_disc_min_cos_zenith_max_zenith():
     times = pd.DatetimeIndex(['2016-07-19 06:11:00'], tz='America/Phoenix')
     out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times)
     expected = pd.DataFrame(np.array(
-        [[0.00000000e+00, 1.16046346e-02, 3.63954476e+01]]),
+        [[0.00000000e+00, 1.16046346e-02, 12.0]]),
         columns=columns, index=times)
     assert_frame_equal(out, expected)
 
+    # max_zenith and/or max_airmass keep these results reasonable
     out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
                           min_cos_zenith=0)
     expected = pd.DataFrame(np.array(
-        [[0.00000000e+00, 1.0, 3.63954476e+01]]),
+        [[0.00000000e+00, 1.0, 12.0]]),
         columns=columns, index=times)
     assert_frame_equal(out, expected)
 
+    # still get reasonable values because of max_airmass=12 limit
     out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
                           max_zenith=100)
     expected = pd.DataFrame(np.array(
-        [[6.68577449e+03, 1.16046346e-02, 3.63954476e+01]]),
+        [[0., 1.16046346e-02, 12.0]]),
         columns=columns, index=times)
     assert_frame_equal(out, expected)
 
+    # still get reasonable values because of max_airmass=12 limit
     out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
                           min_cos_zenith=0, max_zenith=100)
     expected = pd.DataFrame(np.array(
-        [[7.21238390e+03, 1.00000000e+00, 3.63954476e+01]]),
+        [[277.50185968, 1.0, 12.0]]),
+        columns=columns, index=times)
+    assert_frame_equal(out, expected)
+
+    # max_zenith keeps this result reasonable
+    out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
+                          min_cos_zenith=0, max_airmass=100)
+    expected = pd.DataFrame(np.array(
+        [[0.00000000e+00, 1.0, 36.39544757]]),
+        columns=columns, index=times)
+    assert_frame_equal(out, expected)
+
+    # allow zenith to be close to 90 and airmass to be infinite
+    # and we get crazy values
+    out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
+                          max_zenith=100, max_airmass=100)
+    expected = pd.DataFrame(np.array(
+        [[6.68577449e+03, 1.16046346e-02, 3.63954476e+01]]),
+        columns=columns, index=times)
+    assert_frame_equal(out, expected)
+
+    # allow min cos zenith to be 0, zenith to be close to 90,
+    # and airmass to be very big and we get even higher DNI values
+    out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
+                          min_cos_zenith=0, max_zenith=100, max_airmass=100)
+    expected = pd.DataFrame(np.array(
+        [[7.21238390e+03, 1., 3.63954476e+01]]),
         columns=columns, index=times)
     assert_frame_equal(out, expected)
 
@@ -497,13 +526,18 @@ def test_dirint_min_cos_zenith_max_zenith():
     expected = pd.Series([0.0, 0.0], index=times, name='dni')
     assert_series_equal(out, expected)
 
-    out = irradiance.dirint(ghi, solar_zenith, times, max_zenith=100)
-    expected = pd.Series([862.198, 848.387], index=times, name='dni')
+    out = irradiance.dirint(ghi, solar_zenith, times, max_zenith=90)
+    expected = pd.Series([0.0, 0.0], index=times, name='dni')
+    assert_series_equal(out, expected, check_less_precise=True)
+
+    out = irradiance.dirint(ghi, solar_zenith, times, min_cos_zenith=0,
+                            max_zenith=90)
+    expected = pd.Series([0.0, 144.264507], index=times, name='dni')
     assert_series_equal(out, expected, check_less_precise=True)
 
     out = irradiance.dirint(ghi, solar_zenith, times, min_cos_zenith=0,
                             max_zenith=100)
-    expected = pd.Series([147655.5994, 3749.8542], index=times, name='dni')
+    expected = pd.Series([0.0, 144.264507], index=times, name='dni')
     assert_series_equal(out, expected, check_less_precise=True)
 
 
@@ -670,13 +704,13 @@ def test_dirindex_min_cos_zenith_max_zenith():
     assert_series_equal(out, expected)
 
     out = irradiance.dirindex(ghi, ghi_clearsky, dni_clearsky, solar_zenith,
-                              times, max_zenith=100)
-    expected = pd.Series([0., 5.], index=times)
+                              times, max_zenith=90)
+    expected = pd.Series([nan, nan], index=times)
     assert_series_equal(out, expected)
 
     out = irradiance.dirindex(ghi, ghi_clearsky, dni_clearsky, solar_zenith,
                               times, min_cos_zenith=0, max_zenith=100)
-    expected = pd.Series([0., 5.], index=times)
+    expected = pd.Series([nan, 5.], index=times)
     assert_series_equal(out, expected)