Implement pvlib.irradiance.louche decomposition model

pvlib/irradiance.py

@@ -3117,3 +3117,70 @@
                                      'dhi': dhi,
                                      'dni': dni})
     return component_sum_df
+
+
+def pvl_louche(ghi, solar_zenith, datetime_or_doy):
+    """
+    Determine DNI and GHI from GHI using louche model.
+
+    Parameters
+    ----------
+    ghi : Series
+        Global horizontal irradiance. [W/m^2]
+
+    zenith : Series
+        True (not refraction-corrected) zenith angles in decimal
+        degrees. Angles must be >=0 and <=180.
+
+    datetime_or_doy : numeric, pandas.DatetimeIndex
+        Day of year or array of days of year e.g.
+        pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
+
+    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] Louche A, Notton G, Poggi P, Simmonnot G. Correlations for direct
+    normal and global horizontal irradiation on French Mediterranean site.
+    Solar Energy
+    1991;46:261-6
+
+    """
+    bool = np.logical_or(solar_zenith > 180, solar_zenith < 0)
+    solar_zenith = np.where(bool, np.NaN, solar_zenith)
+
+    if np.isscalar(datetime_or_doy):
+        bool = (np.any(datetime_or_doy > 366 or datetime_or_doy < 1, axis=0))
+        print(bool)
+        datetime_or_doy = np.where(bool, np.NaN, datetime_or_doy)
+
+    # this is the I0 calculation from the reference
+    # SSC uses solar constant = 1366.1
+    I0 = get_extra_radiation(datetime_or_doy)
+
+    I0h = I0*tools.cosd(solar_zenith)
+    Kt = ghi/I0h
+    pd.Series.clip(Kt, 0, inplace=True)
+    kb = -10.627*Kt**5 + 15.307*Kt**4 - 5.205 * \
+        Kt**3 + 0.994*Kt**2 - 0.059*Kt + 0.002
+    dni = kb*I0
+    dhi = ghi-dni*tools.cosd(solar_zenith)
+
+    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

pvlib/tests/test_irradiance.py

@@ -246,6 +246,7 @@ def test_haydavies_components(irrad_data, ephem_data, dni_et):
     assert_allclose(result['horizon'], expected['horizon'][-1], atol=1e-4)
     assert isinstance(result, dict)
 
+
 def test_reindl(irrad_data, ephem_data, dni_et):
     result = irradiance.reindl(
         40, 180, irrad_data['dhi'], irrad_data['dni'], irrad_data['ghi'],
@@ -903,8 +904,12 @@ def test_dirindex(times):
     assert np.allclose(out, expected_out, rtol=tolerance, atol=0,
                        equal_nan=True)
     tol_dirint = 0.2
-    assert np.allclose(out.values, dirint_close_values, rtol=tol_dirint, atol=0,
-                       equal_nan=True)
+    assert np.allclose(
+        out.values,
+        dirint_close_values,
+        rtol=tol_dirint,
+        atol=0,
+        equal_nan=True)
 
 
 def test_dirindex_min_cos_zenith_max_zenith():
@@ -1203,3 +1208,22 @@ def test_complete_irradiance():
                                        dhi=None,
                                        dni=i.dni,
                                        dni_clear=clearsky.dni)
+
+
+def test_pvl_louche():
+
+    index = pd.DatetimeIndex(['20190101']*3 + ['20190620']*3)
+    ghi = pd.Series([0, 50, 1000, 1000, 1000, 1000], index=index)
+    zenith = pd.Series([120, 85, 10, 10, 182, -2], index=index)
+    expected = pd.DataFrame(np.array(
+        [[2.827964,     1.413982,        -0.000000],
+         [130.089669, 38.661938, 0.405724],
+         [828.498650, 184.088106, 0.718133],
+         [887.407348, 126.074364, 0.768214],
+         [np.NaN, np.NaN, np.NaN],
+         [np.NaN, np.NaN, np.NaN]]),
+        columns=['dni', 'dhi', 'kt'], index=index)
+
+    out = irradiance.pvl_louche(ghi, zenith, index)
+
+    assert assert_frame_equal(out, expected) is None