From 0b736d4121fc1d05e116dfd38cfa1a4a5a86ead1 Mon Sep 17 00:00:00 2001
From: Abhishek Parikh <abhishek.parikh2@dnv.com>
Date: Thu, 24 Jun 2021 21:19:50 -0700
Subject: [PATCH 01/23] first experimental commit

---
 pvlib/snow.py | 68 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 68 insertions(+)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index bf40c7b995..a3116f5a46 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -185,3 +185,71 @@ def dc_loss_nrel(snow_coverage, num_strings):
        Available at https://www.nrel.gov/docs/fy18osti/67399.pdf
     '''
     return np.ceil(snow_coverage * num_strings) / num_strings
+
+def townsend_Se(S, N):
+    '''
+    Calculates effective snow for a given month based upon the total snowfall
+    received in a month in inches and the number of events where snowfall is greater
+    than 1 inch
+
+    Parameters
+    ----------
+    S : numeric
+        Snowfall in inches received in a month
+
+    N: numeric
+        Number of snowfall events with snowfall > 1"
+
+    Returns
+    -------
+    effective_snowfall : numeric
+        Effective snowfall as defined in the townsend model 
+
+    References
+    ----------
+    .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An 
+            update from two winters of measurements in the SIERRA. Conference 
+            Record of the IEEE Photovoltaic Specialists Conference. 
+            003231-003236. 10.1109/PVSC.2011.6186627. 
+       Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
+
+    '''
+    return (0.5 * S * (1 + 1/N))
+
+def townsend_snow_loss_model(x):
+    '''
+    Loss, % = C1 x Se’ x cos2
+    (tilt) x GIT x RH / TAIR2
+    / POA0.67Eqn. 3
+
+    Parameters
+    ----------
+    snow_coverage : numeric
+        The fraction of row slant height covered by snow at each time step.
+
+    num_strings: int
+        The number of parallel-connected strings along a row slant height.
+
+    Returns
+    -------
+    loss : numeric
+        fraction of DC capacity loss due to snow coverage at each time step.
+
+    References
+    ----------
+    .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An 
+            update from two winters of measurements in the SIERRA. Conference 
+            Record of the IEEE Photovoltaic Specialists Conference. 
+            003231-003236. 10.1109/PVSC.2011.6186627. 
+       Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
+    '''
+    C1 = 5.7e04
+    C2 = 0.51
+
+
+    gamma = [R*Se’*cos(tilt)]/[(H2 – Se’2)/2*tan(P)] 
+    
+    GIT = 1 - C2 * np.exp(-gamma)
+    loss = C1 * Se_ * (np.cos(tilt))**2 * GIT * RH / T_air**2 / POA**0.67
+    
+    return x**2

From 3f2c40dfd98c0858e3f27169e36a4954ec7b39a6 Mon Sep 17 00:00:00 2001
From: Abhishek Parikh <abhishek.parikh2@dnv.com>
Date: Fri, 2 Jul 2021 09:16:35 -0700
Subject: [PATCH 02/23] Added numpy array support

---
 pvlib/snow.py            | 63 ++++++++++++++++++++++++++++++----------
 pvlib/tests/test_snow.py | 23 +++++++++++++++
 2 files changed, 71 insertions(+), 15 deletions(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index a3116f5a46..e6f0fa1adc 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -214,26 +214,52 @@ def townsend_Se(S, N):
        Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
 
     '''
-    return (0.5 * S * (1 + 1/N))
+    return(np.where(N>0, 0.5 * S * (1 + 1/N), 0))
 
-def townsend_snow_loss_model(x):
+def townsend_snow_loss_model(S, N, tilt, RH, T_air, POA, R, H, P=40):
     '''
-    Loss, % = C1 x Se’ x cos2
-    (tilt) x GIT x RH / TAIR2
-    / POA0.67Eqn. 3
+    Calculates monthly snow loss based on a generalized monthly snow loss model discussed in [1]_.
 
     Parameters
     ----------
-    snow_coverage : numeric
-        The fraction of row slant height covered by snow at each time step.
+    S : numeric
+        Inches of snow received in the current month
 
-    num_strings: int
-        The number of parallel-connected strings along a row slant height.
+    N : numeric
+        Number of snowfall events with snowfall > 1"
+
+    tilt : numeric
+        Array tilt in degrees
+
+    RH : numeric
+        Relative humidity in percentage
+
+    T_air : numeric
+        Ambient temperature in celcius
+
+    POA : numeric
+        Plane of array irradiance in kWh/m2/month
+
+    R : numeric
+        Row length in the slanted plane of array dimension in inches
+
+    H : numeric
+        Drop height from array edge to ground in inches
+
+    P : numeric
+        piled snow angle, assumed to stabilize at 40° , the midpoint of 
+        25°-55° avalanching slope angles
+
+    S_prev : numeric
+        Inches of snow received in the previous month
+
+    N_prev : numeric
+        Number of 1" or greater snow events in the previous month
 
     Returns
     -------
     loss : numeric
-        fraction of DC capacity loss due to snow coverage at each time step.
+        Average monthly DC capacity loss in percentage due to snow coverage
 
     References
     ----------
@@ -243,13 +269,20 @@ def townsend_snow_loss_model(x):
             003231-003236. 10.1109/PVSC.2011.6186627. 
        Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
     '''
+
     C1 = 5.7e04
     C2 = 0.51
-
-
-    gamma = [R*Se’*cos(tilt)]/[(H2 – Se’2)/2*tan(P)] 
+    
+    S_prev = np.roll(S,1)
+    N_prev = np.roll(N,1)
+    
+    Se = townsend_Se(S, N)
+    Se_prev = townsend_Se(S_prev, N_prev)
+    
+    Se_weighted = 1/3 * Se_prev + 2/3 * Se
+    gamma = (R * Se_weighted * np.cos(np.deg2rad(tilt)))/(np.clip((H**2 - Se_weighted**2),a_min=0.01,a_max=None)/2/np.tan(np.deg2rad(P)))
     
     GIT = 1 - C2 * np.exp(-gamma)
-    loss = C1 * Se_ * (np.cos(tilt))**2 * GIT * RH / T_air**2 / POA**0.67
+    loss = C1 * Se_weighted * (np.cos(np.deg2rad(tilt)))**2 * GIT * RH / (T_air+273.15)**2 / POA**0.67
     
-    return x**2
+    return (np.round(loss,2))
\ No newline at end of file
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index 3db56c8f61..77c1f72399 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -95,3 +95,26 @@ def test_dc_loss_nrel():
     expected = pd.Series([1, 1, .5, .625, .25, .5, 0])
     actual = snow.dc_loss_nrel(snow_coverage, num_strings)
     assert_series_equal(expected, actual)
+
+def test_townsend_Se():
+    S = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
+    N = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
+    expected = np.array([7.5, 7.5, 5, 0, 0, 0, 0, 0, 0, 0, 3.75, 6.66666667])
+    actual = snow.townsend_Se(S, N)
+    np.testing.assert_allclose(expected, actual, rtol=1e-07)
+
+def test_townsend_snow_loss_model():
+    S = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
+    N = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
+    tilt = 20
+    RH = np.array([80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80])
+    T_air = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+    POA = np.array([350, 350, 350, 350, 350, 350, 350, 350, 350, 350, 350, 350])
+    P = 40
+    R = 100
+    H = 10
+    expected = np.array([7.7, 7.99, 6.22, 1.72, 0, 0, 0, 0, 0, 0, 2.64, 6.07])
+    actual = snow.townsend_snow_loss_model(S,N,tilt,RH,T_air,POA,R,
+                                            H,P)
+    np.testing.assert_allclose(expected, actual, rtol=1e-07)
+

From 974d516d44dd508f2456e81334bc6a9c5b772087 Mon Sep 17 00:00:00 2001
From: Abhishek Parikh <abhishek.parikh2@dnv.com>
Date: Thu, 15 Jul 2021 20:16:35 -0700
Subject: [PATCH 03/23] changed var names

---
 pvlib/snow.py            | 42 ++++++++++++++++++----------------------
 pvlib/tests/test_snow.py | 18 ++++++++---------
 2 files changed, 28 insertions(+), 32 deletions(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index e6f0fa1adc..2d839c006c 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -216,46 +216,40 @@ def townsend_Se(S, N):
     '''
     return(np.where(N>0, 0.5 * S * (1 + 1/N), 0))
 
-def townsend_snow_loss_model(S, N, tilt, RH, T_air, POA, R, H, P=40):
+def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air, poa_global, row_len, H, P=40):
     '''
     Calculates monthly snow loss based on a generalized monthly snow loss model discussed in [1]_.
 
     Parameters
     ----------
-    S : numeric
-        Inches of snow received in the current month
+    snow_total : numeric
+        Inches of snow received in the current month. Referred as S in the paper
 
-    N : numeric
-        Number of snowfall events with snowfall > 1"
+    snow_events : numeric
+        Number of snowfall events with snowfall > 1". Referred as N in the paper
 
     tilt : numeric
         Array tilt in degrees
 
-    RH : numeric
+    relative_humidity : numeric
         Relative humidity in percentage
 
-    T_air : numeric
+    temp_air : numeric
         Ambient temperature in celcius
 
-    POA : numeric
+    poa_global : numeric
         Plane of array irradiance in kWh/m2/month
 
-    R : numeric
+    row_len : float
         Row length in the slanted plane of array dimension in inches
 
-    H : numeric
+    H : float
         Drop height from array edge to ground in inches
 
-    P : numeric
+    P : float
         piled snow angle, assumed to stabilize at 40° , the midpoint of 
         25°-55° avalanching slope angles
 
-    S_prev : numeric
-        Inches of snow received in the previous month
-
-    N_prev : numeric
-        Number of 1" or greater snow events in the previous month
-
     Returns
     -------
     loss : numeric
@@ -273,16 +267,18 @@ def townsend_snow_loss_model(S, N, tilt, RH, T_air, POA, R, H, P=40):
     C1 = 5.7e04
     C2 = 0.51
     
-    S_prev = np.roll(S,1)
-    N_prev = np.roll(N,1)
+    snow_total_prev = np.roll(snow_total,1)
+    snow_events_prev = np.roll(snow_events,1)
     
-    Se = townsend_Se(S, N)
-    Se_prev = townsend_Se(S_prev, N_prev)
+    Se = townsend_Se(snow_total, snow_events)
+    Se_prev = townsend_Se(snow_total_prev, snow_events_prev)
     
     Se_weighted = 1/3 * Se_prev + 2/3 * Se
-    gamma = (R * Se_weighted * np.cos(np.deg2rad(tilt)))/(np.clip((H**2 - Se_weighted**2),a_min=0.01,a_max=None)/2/np.tan(np.deg2rad(P)))
+    gamma = (row_len * Se_weighted * np.cos(np.deg2rad(tilt)))/ \
+        (np.clip((H**2 - Se_weighted**2),a_min=0.01,a_max=None)/2/np.tan(np.deg2rad(P)))
     
     GIT = 1 - C2 * np.exp(-gamma)
-    loss = C1 * Se_weighted * (np.cos(np.deg2rad(tilt)))**2 * GIT * RH / (T_air+273.15)**2 / POA**0.67
+    loss = C1 * Se_weighted * (np.cos(np.deg2rad(tilt)))**2 * GIT * relative_humidity / \
+        (temp_air+273.15)**2 / poa_global**0.67
     
     return (np.round(loss,2))
\ No newline at end of file
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index 77c1f72399..4aa30b62e8 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -103,18 +103,18 @@ def test_townsend_Se():
     actual = snow.townsend_Se(S, N)
     np.testing.assert_allclose(expected, actual, rtol=1e-07)
 
-def test_townsend_snow_loss_model():
-    S = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
-    N = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
+def test_loss_townsend():
+    snow_total = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
+    snow_events = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
     tilt = 20
-    RH = np.array([80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80])
-    T_air = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
-    POA = np.array([350, 350, 350, 350, 350, 350, 350, 350, 350, 350, 350, 350])
+    relative_humidity = np.array([80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80])
+    temp_air = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+    poa_global = np.array([350, 350, 350, 350, 350, 350, 350, 350, 350, 350, 350, 350])
     P = 40
-    R = 100
+    row_len = 100
     H = 10
     expected = np.array([7.7, 7.99, 6.22, 1.72, 0, 0, 0, 0, 0, 0, 2.64, 6.07])
-    actual = snow.townsend_snow_loss_model(S,N,tilt,RH,T_air,POA,R,
-                                            H,P)
+    actual = snow.loss_townsend(snow_total,snow_events,tilt,relative_humidity,temp_air,
+                                            poa_global,row_len,H,P)
     np.testing.assert_allclose(expected, actual, rtol=1e-07)
 

From 655fa79574962e3a5cdd8897d3936fa3d6b04c4b Mon Sep 17 00:00:00 2001
From: Abhishek Parikh <abhishek.parikh2@dnv.com>
Date: Mon, 19 Jul 2021 17:06:29 -0700
Subject: [PATCH 04/23] changed townsend_Se to private

---
 docs/sphinx/source/api.rst |  1 +
 pvlib/snow.py              | 25 ++++++++++++++-----------
 pvlib/tests/test_snow.py   |  4 ++--
 3 files changed, 17 insertions(+), 13 deletions(-)

diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 0e5438a7fc..6bf04af1af 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -403,6 +403,7 @@ Snow
    snow.coverage_nrel
    snow.fully_covered_nrel
    snow.dc_loss_nrel
+   snow.loss_townsend
 
 Soiling
 -------
diff --git a/pvlib/snow.py b/pvlib/snow.py
index 2d839c006c..0efb62e006 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -186,11 +186,12 @@ def dc_loss_nrel(snow_coverage, num_strings):
     '''
     return np.ceil(snow_coverage * num_strings) / num_strings
 
-def townsend_Se(S, N):
+
+def _townsend_Se(S, N):
     '''
     Calculates effective snow for a given month based upon the total snowfall
-    received in a month in inches and the number of events where snowfall is greater
-    than 1 inch
+    received in a month in inches and the number of events where snowfall is 
+    greater than 1 inch
 
     Parameters
     ----------
@@ -210,12 +211,13 @@ def townsend_Se(S, N):
     .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An 
             update from two winters of measurements in the SIERRA. Conference 
             Record of the IEEE Photovoltaic Specialists Conference. 
-            003231-003236. 10.1109/PVSC.2011.6186627. 
+            003231-003236. :doi:`10.1109/PVSC.2011.6186627`
        Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
 
     '''
     return(np.where(N>0, 0.5 * S * (1 + 1/N), 0))
 
+
 def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air, poa_global, row_len, H, P=40):
     '''
     Calculates monthly snow loss based on a generalized monthly snow loss model discussed in [1]_.
@@ -235,10 +237,10 @@ def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air, po
         Relative humidity in percentage
 
     temp_air : numeric
-        Ambient temperature in celcius
+        Ambient temperature [C]
 
     poa_global : numeric
-        Plane of array irradiance in kWh/m2/month
+        Plane of array insolation in kWh/m2/month
 
     row_len : float
         Row length in the slanted plane of array dimension in inches
@@ -270,15 +272,16 @@ def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air, po
     snow_total_prev = np.roll(snow_total,1)
     snow_events_prev = np.roll(snow_events,1)
     
-    Se = townsend_Se(snow_total, snow_events)
-    Se_prev = townsend_Se(snow_total_prev, snow_events_prev)
+    Se = _townsend_Se(snow_total, snow_events)
+    Se_prev = _townsend_Se(snow_total_prev, snow_events_prev)
     
     Se_weighted = 1/3 * Se_prev + 2/3 * Se
     gamma = (row_len * Se_weighted * np.cos(np.deg2rad(tilt)))/ \
-        (np.clip((H**2 - Se_weighted**2),a_min=0.01,a_max=None)/2/np.tan(np.deg2rad(P)))
+        (np.clip((H**2 - Se_weighted**2),a_min=0.01,a_max=None)/2/ \
+            np.tan(np.deg2rad(P)))
     
     GIT = 1 - C2 * np.exp(-gamma)
-    loss = C1 * Se_weighted * (np.cos(np.deg2rad(tilt)))**2 * GIT * relative_humidity / \
-        (temp_air+273.15)**2 / poa_global**0.67
+    loss = C1 * Se_weighted * (np.cos(np.deg2rad(tilt)))**2 * GIT * \
+        relative_humidity / (temp_air+273.15)**2 / poa_global**0.67
     
     return (np.round(loss,2))
\ No newline at end of file
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index 4aa30b62e8..cf7ad879b7 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -96,11 +96,11 @@ def test_dc_loss_nrel():
     actual = snow.dc_loss_nrel(snow_coverage, num_strings)
     assert_series_equal(expected, actual)
 
-def test_townsend_Se():
+def test__townsend_Se():
     S = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
     N = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
     expected = np.array([7.5, 7.5, 5, 0, 0, 0, 0, 0, 0, 0, 3.75, 6.66666667])
-    actual = snow.townsend_Se(S, N)
+    actual = snow._townsend_Se(S, N)
     np.testing.assert_allclose(expected, actual, rtol=1e-07)
 
 def test_loss_townsend():

From c0c8b4d866d830b2ef0c1342c509846a0502443d Mon Sep 17 00:00:00 2001
From: abhisheksparikh <abhi.parikh305@gmail.com>
Date: Wed, 16 Mar 2022 21:35:23 -0700
Subject: [PATCH 05/23] removed snow.loss_townsend in api.rst

---
 docs/sphinx/source/api.rst | 1 -
 1 file changed, 1 deletion(-)

diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 6bf04af1af..0e5438a7fc 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -403,7 +403,6 @@ Snow
    snow.coverage_nrel
    snow.fully_covered_nrel
    snow.dc_loss_nrel
-   snow.loss_townsend
 
 Soiling
 -------

From 78716f7bfe504197df305e4cc7b7891933f4b075 Mon Sep 17 00:00:00 2001
From: abhisheksparikh <abhi.parikh305@gmail.com>
Date: Wed, 16 Mar 2022 21:52:57 -0700
Subject: [PATCH 06/23] added loss_townsend description in
 effects_on_pv_system_output.rst

---
 docs/sphinx/source/reference/effects_on_pv_system_output.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/sphinx/source/reference/effects_on_pv_system_output.rst b/docs/sphinx/source/reference/effects_on_pv_system_output.rst
index 92efa946b4..62c1035860 100644
--- a/docs/sphinx/source/reference/effects_on_pv_system_output.rst
+++ b/docs/sphinx/source/reference/effects_on_pv_system_output.rst
@@ -21,6 +21,7 @@ Snow
    snow.coverage_nrel
    snow.fully_covered_nrel
    snow.dc_loss_nrel
+   snow.loss_townsend
 
 Soiling
 -------

From 5e06a40f289fc098583c705875b07dddd29286b3 Mon Sep 17 00:00:00 2001
From: abhisheksparikh <abhi.parikh305@gmail.com>
Date: Wed, 16 Mar 2022 22:13:15 -0700
Subject: [PATCH 07/23] fixed stickler checks

---
 pvlib/snow.py            | 58 +++++++++++++++++++++-------------------
 pvlib/tests/test_snow.py | 14 ++++++----
 2 files changed, 40 insertions(+), 32 deletions(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index 0efb62e006..db96bb1887 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -190,7 +190,7 @@ def dc_loss_nrel(snow_coverage, num_strings):
 def _townsend_Se(S, N):
     '''
     Calculates effective snow for a given month based upon the total snowfall
-    received in a month in inches and the number of events where snowfall is 
+    received in a month in inches and the number of events where snowfall is
     greater than 1 inch
 
     Parameters
@@ -204,31 +204,35 @@ def _townsend_Se(S, N):
     Returns
     -------
     effective_snowfall : numeric
-        Effective snowfall as defined in the townsend model 
+        Effective snowfall as defined in the townsend model
 
     References
     ----------
-    .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An 
-            update from two winters of measurements in the SIERRA. Conference 
-            Record of the IEEE Photovoltaic Specialists Conference. 
+    .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An
+            update from two winters of measurements in the SIERRA. Conference
+            Record of the IEEE Photovoltaic Specialists Conference.
             003231-003236. :doi:`10.1109/PVSC.2011.6186627`
        Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
 
-    '''
-    return(np.where(N>0, 0.5 * S * (1 + 1/N), 0))
+    '''# noqa
+    return(np.where(N > 0, 0.5 * S * (1 + 1/N), 0))
 
 
-def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air, poa_global, row_len, H, P=40):
+def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air,
+                  poa_global, row_len, H, P=40):
     '''
-    Calculates monthly snow loss based on a generalized monthly snow loss model discussed in [1]_.
+    Calculates monthly snow loss based on a generalized monthly snow loss model
+    discussed in [1]_.
 
     Parameters
     ----------
     snow_total : numeric
-        Inches of snow received in the current month. Referred as S in the paper
+        Inches of snow received in the current month. Referred as S in the
+        paper
 
     snow_events : numeric
-        Number of snowfall events with snowfall > 1". Referred as N in the paper
+        Number of snowfall events with snowfall > 1". Referred as N in the
+        paper
 
     tilt : numeric
         Array tilt in degrees
@@ -249,7 +253,7 @@ def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air, po
         Drop height from array edge to ground in inches
 
     P : float
-        piled snow angle, assumed to stabilize at 40° , the midpoint of 
+        piled snow angle, assumed to stabilize at 40° , the midpoint of
         25°-55° avalanching slope angles
 
     Returns
@@ -259,29 +263,29 @@ def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air, po
 
     References
     ----------
-    .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An 
-            update from two winters of measurements in the SIERRA. Conference 
-            Record of the IEEE Photovoltaic Specialists Conference. 
-            003231-003236. 10.1109/PVSC.2011.6186627. 
+    .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An
+            update from two winters of measurements in the SIERRA. Conference
+            Record of the IEEE Photovoltaic Specialists Conference.
+            003231-003236. 10.1109/PVSC.2011.6186627.
        Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
-    '''
+    '''# noqa
 
     C1 = 5.7e04
     C2 = 0.51
-    
-    snow_total_prev = np.roll(snow_total,1)
-    snow_events_prev = np.roll(snow_events,1)
-    
+
+    snow_total_prev = np.roll(snow_total, 1)
+    snow_events_prev = np.roll(snow_events, 1)
+
     Se = _townsend_Se(snow_total, snow_events)
     Se_prev = _townsend_Se(snow_total_prev, snow_events_prev)
-    
+
     Se_weighted = 1/3 * Se_prev + 2/3 * Se
-    gamma = (row_len * Se_weighted * np.cos(np.deg2rad(tilt)))/ \
-        (np.clip((H**2 - Se_weighted**2),a_min=0.01,a_max=None)/2/ \
+    gamma = (row_len * Se_weighted * np.cos(np.deg2rad(tilt))) / \
+        (np.clip((H**2 - Se_weighted**2), a_min=0.01, a_max=None) / 2 /
             np.tan(np.deg2rad(P)))
-    
+
     GIT = 1 - C2 * np.exp(-gamma)
     loss = C1 * Se_weighted * (np.cos(np.deg2rad(tilt)))**2 * GIT * \
         relative_humidity / (temp_air+273.15)**2 / poa_global**0.67
-    
-    return (np.round(loss,2))
\ No newline at end of file
+
+    return (np.round(loss, 2))
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index cf7ad879b7..c580ef9d07 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -96,6 +96,7 @@ def test_dc_loss_nrel():
     actual = snow.dc_loss_nrel(snow_coverage, num_strings)
     assert_series_equal(expected, actual)
 
+
 def test__townsend_Se():
     S = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
     N = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
@@ -103,18 +104,21 @@ def test__townsend_Se():
     actual = snow._townsend_Se(S, N)
     np.testing.assert_allclose(expected, actual, rtol=1e-07)
 
+
 def test_loss_townsend():
     snow_total = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
     snow_events = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
     tilt = 20
-    relative_humidity = np.array([80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80])
+    relative_humidity = np.array([80, 80, 80, 80, 80, 80, 80, 80, 80, 80,
+                                  80, 80])
     temp_air = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
-    poa_global = np.array([350, 350, 350, 350, 350, 350, 350, 350, 350, 350, 350, 350])
+    poa_global = np.array([350, 350, 350, 350, 350, 350, 350, 350, 350, 350,
+                           350, 350])
     P = 40
     row_len = 100
     H = 10
     expected = np.array([7.7, 7.99, 6.22, 1.72, 0, 0, 0, 0, 0, 0, 2.64, 6.07])
-    actual = snow.loss_townsend(snow_total,snow_events,tilt,relative_humidity,temp_air,
-                                            poa_global,row_len,H,P)
+    actual = snow.loss_townsend(snow_total, snow_events, tilt,
+                                relative_humidity, temp_air,
+                                poa_global, row_len, H, P)
     np.testing.assert_allclose(expected, actual, rtol=1e-07)
-

From a5c53151067e03327991732e226b22925901a89c Mon Sep 17 00:00:00 2001
From: abhisheksparikh <abhi.parikh305@gmail.com>
Date: Fri, 18 Mar 2022 00:34:03 -0700
Subject: [PATCH 08/23] removed rounding of loss and changed to 0-1 range

Several other small changes - variable names change, comment change - in response to Kevin's review notes
---
 docs/sphinx/source/whatsnew/v0.9.1.rst |  3 ++
 pvlib/snow.py                          | 65 +++++++++++++-------------
 pvlib/tests/test_snow.py               | 18 +++----
 3 files changed, 46 insertions(+), 40 deletions(-)

diff --git a/docs/sphinx/source/whatsnew/v0.9.1.rst b/docs/sphinx/source/whatsnew/v0.9.1.rst
index 8aa1abb71f..30ee76fb1a 100644
--- a/docs/sphinx/source/whatsnew/v0.9.1.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.1.rst
@@ -26,6 +26,8 @@ Enhancements
 * Added ``map_variables`` option to :func:`~pvlib.iotools.read_crn` (:pull:`1368`)
 * Added :py:func:`pvlib.temperature.prilliman` for modeling cell temperature
   at short time steps (:issue:`1081`, :pull:`1391`)
+* Added Townsend Powers Snow loss model in :py:func:`pvlib.snow`
+  (:issue:`1246`, :pull:`1251`)
 
 Bug fixes
 ~~~~~~~~~
@@ -78,3 +80,4 @@ Contributors
 * Saurabh Aneja (:ghuser:`spaneja`)
 * Jack Kelly (:ghuser:`JackKelly`)
 * Somasree Majumder(:ghuser:`soma2000-lang`)
+* Abhishek Parikh (:ghuser:`abhisheksparikh`)
diff --git a/pvlib/snow.py b/pvlib/snow.py
index db96bb1887..c003c5531a 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -5,7 +5,7 @@
 
 import numpy as np
 import pandas as pd
-from pvlib.tools import sind
+from pvlib.tools import sind, cosd, tand
 
 
 def _time_delta_in_hours(times):
@@ -196,10 +196,10 @@ def _townsend_Se(S, N):
     Parameters
     ----------
     S : numeric
-        Snowfall in inches received in a month
+        Snowfall in inches received in a month [in]
 
     N: numeric
-        Number of snowfall events with snowfall > 1"
+        Number of snowfall events with snowfall > 1" [-]
 
     Returns
     -------
@@ -209,17 +209,18 @@ def _townsend_Se(S, N):
     References
     ----------
     .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An
-            update from two winters of measurements in the SIERRA. Conference
-            Record of the IEEE Photovoltaic Specialists Conference.
-            003231-003236. :doi:`10.1109/PVSC.2011.6186627`
+       update from two winters of measurements in the SIERRA. Conference
+       Record of the IEEE Photovoltaic Specialists Conference.
+       003231-003236. :doi:`10.1109/PVSC.2011.6186627`
        Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
 
-    '''# noqa
+    '''  # noqa: E501
     return(np.where(N > 0, 0.5 * S * (1 + 1/N), 0))
 
 
-def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air,
-                  poa_global, row_len, H, P=40):
+def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
+                  temp_air, poa_global, slant_height, lower_edge_drop_height,
+                  angle_of_repose=40):
     '''
     Calculates monthly snow loss based on a generalized monthly snow loss model
     discussed in [1]_.
@@ -228,47 +229,47 @@ def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air,
     ----------
     snow_total : numeric
         Inches of snow received in the current month. Referred as S in the
-        paper
+        paper [in]
 
     snow_events : numeric
         Number of snowfall events with snowfall > 1". Referred as N in the
-        paper
+        paper [-]
 
-    tilt : numeric
-        Array tilt in degrees
+    surface_tilt : numeric
+        Array surface_tilt [deg]
 
     relative_humidity : numeric
-        Relative humidity in percentage
+        Relative humidity [%]
 
     temp_air : numeric
-        Ambient temperature [C]
+        Ambient temperature [°C]
 
     poa_global : numeric
-        Plane of array insolation in kWh/m2/month
+        Plane of array insolation [kWh/m2/month]
 
-    row_len : float
-        Row length in the slanted plane of array dimension in inches
+    slant_height : float
+        Row length in the slanted plane of array dimension [in]
 
-    H : float
-        Drop height from array edge to ground in inches
+    lower_edge_drop_height : float
+        Drop height from array edge to ground [in]
 
     P : float
         piled snow angle, assumed to stabilize at 40° , the midpoint of
-        25°-55° avalanching slope angles
+        25°-55° avalanching slope angles [deg]
 
     Returns
     -------
     loss : numeric
-        Average monthly DC capacity loss in percentage due to snow coverage
+        Average monthly DC capacity loss due to snow coverage [%]
 
     References
     ----------
     .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An
-            update from two winters of measurements in the SIERRA. Conference
-            Record of the IEEE Photovoltaic Specialists Conference.
-            003231-003236. 10.1109/PVSC.2011.6186627.
+       update from two winters of measurements in the SIERRA. Conference
+       Record of the IEEE Photovoltaic Specialists Conference.
+       003231-003236. 10.1109/PVSC.2011.6186627.
        Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
-    '''# noqa
+    '''  # noqa: E501
 
     C1 = 5.7e04
     C2 = 0.51
@@ -280,12 +281,12 @@ def loss_townsend(snow_total, snow_events, tilt, relative_humidity, temp_air,
     Se_prev = _townsend_Se(snow_total_prev, snow_events_prev)
 
     Se_weighted = 1/3 * Se_prev + 2/3 * Se
-    gamma = (row_len * Se_weighted * np.cos(np.deg2rad(tilt))) / \
-        (np.clip((H**2 - Se_weighted**2), a_min=0.01, a_max=None) / 2 /
-            np.tan(np.deg2rad(P)))
+    gamma = (slant_height * Se_weighted * cosd(surface_tilt)) / \
+        (np.clip((lower_edge_drop_height**2 - Se_weighted**2), a_min=0.01,
+         a_max=None) / 2 / tand(angle_of_repose))
 
     GIT = 1 - C2 * np.exp(-gamma)
-    loss = C1 * Se_weighted * (np.cos(np.deg2rad(tilt)))**2 * GIT * \
-        relative_humidity / (temp_air+273.15)**2 / poa_global**0.67
+    loss = (C1 * Se_weighted * (cosd(surface_tilt))**2 * GIT *
+            relative_humidity / (temp_air+273.15)**2 / poa_global**0.67) / 100
 
-    return (np.round(loss, 2))
+    return loss
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index c580ef9d07..bc71d5537a 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -108,17 +108,19 @@ def test__townsend_Se():
 def test_loss_townsend():
     snow_total = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
     snow_events = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
-    tilt = 20
+    surface_tilt = 20
     relative_humidity = np.array([80, 80, 80, 80, 80, 80, 80, 80, 80, 80,
                                   80, 80])
     temp_air = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
     poa_global = np.array([350, 350, 350, 350, 350, 350, 350, 350, 350, 350,
                            350, 350])
-    P = 40
-    row_len = 100
-    H = 10
-    expected = np.array([7.7, 7.99, 6.22, 1.72, 0, 0, 0, 0, 0, 0, 2.64, 6.07])
-    actual = snow.loss_townsend(snow_total, snow_events, tilt,
+    angle_of_repose = 40
+    slant_height = 100
+    lower_edge_drop_height = 10
+    expected = np.array([0.07696253, 0.07992262, 0.06216201, 0.01715392, 0, 0,
+                         0, 0, 0, 0, 0.02643821, 0.06068194])
+    actual = snow.loss_townsend(snow_total, snow_events, surface_tilt,
                                 relative_humidity, temp_air,
-                                poa_global, row_len, H, P)
-    np.testing.assert_allclose(expected, actual, rtol=1e-07)
+                                poa_global, slant_height,
+                                lower_edge_drop_height, angle_of_repose)
+    np.testing.assert_allclose(expected, actual, rtol=1e-05)

From c2d6d4f450e096f6554d6467d6aeb9a27db20c60 Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Thu, 26 May 2022 12:40:04 -0700
Subject: [PATCH 09/23] implementing changes suggested in PR #1251

---
 docs/sphinx/source/whatsnew/v0.9.1.rst |   4 +-
 docs/sphinx/source/whatsnew/v0.9.2.rst |   4 +
 pvlib/snow.py                          | 105 ++++++++++++++-----------
 pvlib/tests/test_snow.py               |   6 +-
 4 files changed, 69 insertions(+), 50 deletions(-)

diff --git a/docs/sphinx/source/whatsnew/v0.9.1.rst b/docs/sphinx/source/whatsnew/v0.9.1.rst
index a2f0875105..07f0153be4 100644
--- a/docs/sphinx/source/whatsnew/v0.9.1.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.1.rst
@@ -1,3 +1,4 @@
+
 .. _whatsnew_0910:
 
 v0.9.1 (March 29, 2022)
@@ -25,8 +26,6 @@ Enhancements
 * Added ``map_variables`` option to :func:`~pvlib.iotools.read_crn` (:pull:`1368`)
 * Added :py:func:`pvlib.temperature.prilliman` for modeling cell temperature
   at short time steps (:issue:`1081`, :pull:`1391`)
-* Added Townsend Powers Snow loss model in :py:func:`pvlib.snow`
-  (:issue:`1246`, :pull:`1251`)
 
 Bug fixes
 ~~~~~~~~~
@@ -85,7 +84,6 @@ Contributors
 * Jack Kelly (:ghuser:`JackKelly`)
 * Somasree Majumder(:ghuser:`soma2000-lang`)
 * Naman Priyadarshi (:ghuser:`Naman-Priyadarshi`)
-* Abhishek Parikh (:ghuser:`abhisheksparikh`)
 
 * Will Holmgren (:ghuser:`wholmgren`)
 * Mark Mikofksi (:ghuser:`mikofski`)
diff --git a/docs/sphinx/source/whatsnew/v0.9.2.rst b/docs/sphinx/source/whatsnew/v0.9.2.rst
index 05b7ed58ab..ac72453d14 100644
--- a/docs/sphinx/source/whatsnew/v0.9.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.2.rst
@@ -18,6 +18,8 @@ Bug fixes
   where passing localized timezones with large UTC offsets could return
   rise/set/transit times for the wrong day in recent versions of ``ephem``
   (:issue:`1449`, :pull:`1448`)
+* Added Townsend-Powers monthly snow loss model: :py:func:`pvlib.snow.loss_townsend`
+  (:issue:`1246`, :pull:`1251`)
 
 
 Testing
@@ -40,3 +42,5 @@ Contributors
 * Naman Priyadarshi (:ghuser:`Naman-Priyadarshi`)
 * Chencheng Luo (:ghuser:`roger-lcc`)
 * Prajwal Borkar (:ghuser:`PrajwalBorkar`) 
+* Abhishek Parikh (:ghuser:`abhisheksparikh`)
+* Taos Transue (:ghuser:`reepoi`)
diff --git a/pvlib/snow.py b/pvlib/snow.py
index c003c5531a..7f267d395a 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -187,24 +187,23 @@ def dc_loss_nrel(snow_coverage, num_strings):
     return np.ceil(snow_coverage * num_strings) / num_strings
 
 
-def _townsend_Se(S, N):
+def _townsend_effective_snow(snow_load, snow_events):
     '''
-    Calculates effective snow for a given month based upon the total snowfall
-    received in a month in inches and the number of events where snowfall is
-    greater than 1 inch
+    Calculates effective snow using the total snowfall in inches received
+    each month and the number of snowfall events each month.
 
     Parameters
     ----------
-    S : numeric
-        Snowfall in inches received in a month [in]
+    snow_load : array-like
+        Inches of snow received each month. Referred to as S in the paper [in]
 
-    N: numeric
-        Number of snowfall events with snowfall > 1" [-]
+    snow_events : array-like
+        Number of snowfall events each month. Referred to as N in the paper [-]
 
     Returns
     -------
-    effective_snowfall : numeric
-        Effective snowfall as defined in the townsend model
+    effective_snowfall : array-like
+        Effective snowfall as defined in the Townsend model
 
     References
     ----------
@@ -215,37 +214,36 @@ def _townsend_Se(S, N):
        Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
 
     '''  # noqa: E501
-    return(np.where(N > 0, 0.5 * S * (1 + 1/N), 0))
+    snow_events_no_zeros = np.maximum(snow_events, 1)
+    effective_snow = 0.5 * snow_load * (1 + 1 / snow_events_no_zeros)
+    return np.where(snow_events > 0, effective_snow, 0)
 
 
 def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
                   temp_air, poa_global, slant_height, lower_edge_drop_height,
                   angle_of_repose=40):
     '''
-    Calculates monthly snow loss based on a generalized monthly snow loss model
-    discussed in [1]_.
+    Calculates monthly snow loss based on the Townsend monthly snow loss model [1]_.
 
     Parameters
     ----------
-    snow_total : numeric
-        Inches of snow received in the current month. Referred as S in the
-        paper [in]
+    snow_total : array-like
+        Inches of snow received each month. Referred to as S in the paper [in]
 
-    snow_events : numeric
-        Number of snowfall events with snowfall > 1". Referred as N in the
-        paper [-]
+    snow_events : array-like
+        Number of snowfall events each month. Referred to as N in the paper [-]
 
-    surface_tilt : numeric
-        Array surface_tilt [deg]
+    surface_tilt : float
+        Tilt angle of the array [deg]
 
-    relative_humidity : numeric
-        Relative humidity [%]
+    relative_humidity : array-like
+        Monthly average relative humidity [%]
 
-    temp_air : numeric
-        Ambient temperature [°C]
+    temp_air : array-like
+        Monthly average ambient temperature [C]
 
-    poa_global : numeric
-        Plane of array insolation [kWh/m2/month]
+    poa_global : array-like
+        Monthly plane of array insolation [kWh/m2/month]
 
     slant_height : float
         Row length in the slanted plane of array dimension [in]
@@ -253,14 +251,19 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
     lower_edge_drop_height : float
         Drop height from array edge to ground [in]
 
-    P : float
-        piled snow angle, assumed to stabilize at 40° , the midpoint of
+    angle_of_repose : float, default 40
+        piled snow angle, assumed to stabilize at 40°, the midpoint of
         25°-55° avalanching slope angles [deg]
 
     Returns
     -------
-    loss : numeric
-        Average monthly DC capacity loss due to snow coverage [%]
+    loss : array-like
+        Monthly average DC capacity loss fraction due to snow coverage
+
+    Notes
+    -----
+    This model has not been validated for tracking arrays; however, for tracking
+    arrays [1]_ suggests using the maximum rotation angle in place of surface_tilt.
 
     References
     ----------
@@ -277,16 +280,30 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
     snow_total_prev = np.roll(snow_total, 1)
     snow_events_prev = np.roll(snow_events, 1)
 
-    Se = _townsend_Se(snow_total, snow_events)
-    Se_prev = _townsend_Se(snow_total_prev, snow_events_prev)
-
-    Se_weighted = 1/3 * Se_prev + 2/3 * Se
-    gamma = (slant_height * Se_weighted * cosd(surface_tilt)) / \
-        (np.clip((lower_edge_drop_height**2 - Se_weighted**2), a_min=0.01,
-         a_max=None) / 2 / tand(angle_of_repose))
-
-    GIT = 1 - C2 * np.exp(-gamma)
-    loss = (C1 * Se_weighted * (cosd(surface_tilt))**2 * GIT *
-            relative_humidity / (temp_air+273.15)**2 / poa_global**0.67) / 100
-
-    return loss
+    effective_snow = _townsend_effective_snow(snow_total, snow_events)
+    effective_snow_prev = _townsend_effective_snow(snow_total_prev, snow_events_prev)
+    effective_snow_weighted = 1 / 3 * effective_snow_prev + 2 / 3 * effective_snow
+
+    drop_height_clipped = np.maximum(lower_edge_drop_height, 0.01)
+    gamma = (
+        slant_height
+        * effective_snow_weighted
+        * cosd(surface_tilt)
+        / (drop_height_clipped**2 - effective_snow_weighted**2)
+        * 2
+        * tand(angle_of_repose)
+    )
+
+    ground_interference_term = 1 - C2 * np.exp(-gamma)
+    temp_air_kelvin = temp_air + 273.15
+    loss_percentage = (
+        C1
+        * effective_snow_weighted
+        * cosd(surface_tilt)**2
+        * ground_interference_term
+        * relative_humidity
+        / temp_air_kelvin**2
+        / poa_global**0.67
+    )
+
+    return loss_percentage / 100
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index bc71d5537a..cb72ae58b8 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -98,10 +98,10 @@ def test_dc_loss_nrel():
 
 
 def test__townsend_Se():
-    S = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
-    N = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
+    snow_load = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
+    snow_events = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
     expected = np.array([7.5, 7.5, 5, 0, 0, 0, 0, 0, 0, 0, 3.75, 6.66666667])
-    actual = snow._townsend_Se(S, N)
+    actual = snow._townsend_effective_snow(snow_load, snow_events)
     np.testing.assert_allclose(expected, actual, rtol=1e-07)
 
 

From 68c85f7dcd61a66ba39dc1adeb5558c5bf436dac Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Thu, 26 May 2022 12:46:14 -0700
Subject: [PATCH 10/23] removing new line

---
 docs/sphinx/source/whatsnew/v0.9.1.rst | 1 -
 1 file changed, 1 deletion(-)

diff --git a/docs/sphinx/source/whatsnew/v0.9.1.rst b/docs/sphinx/source/whatsnew/v0.9.1.rst
index 07f0153be4..cbde204237 100644
--- a/docs/sphinx/source/whatsnew/v0.9.1.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.1.rst
@@ -84,7 +84,6 @@ Contributors
 * Jack Kelly (:ghuser:`JackKelly`)
 * Somasree Majumder(:ghuser:`soma2000-lang`)
 * Naman Priyadarshi (:ghuser:`Naman-Priyadarshi`)
-
 * Will Holmgren (:ghuser:`wholmgren`)
 * Mark Mikofksi (:ghuser:`mikofski`)
 * Will Hobbs (:ghuser:`williamhobbs`)

From e5118ef2185c55b534feca20a98395e37f9c16d3 Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Thu, 26 May 2022 12:47:28 -0700
Subject: [PATCH 11/23] removing new line

---
 docs/sphinx/source/whatsnew/v0.9.1.rst | 1 -
 1 file changed, 1 deletion(-)

diff --git a/docs/sphinx/source/whatsnew/v0.9.1.rst b/docs/sphinx/source/whatsnew/v0.9.1.rst
index cbde204237..2440f215c7 100644
--- a/docs/sphinx/source/whatsnew/v0.9.1.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.1.rst
@@ -1,4 +1,3 @@
-
 .. _whatsnew_0910:
 
 v0.9.1 (March 29, 2022)

From c253f5943b46a04b9322b19e6006f2187cda2e1c Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Tue, 7 Jun 2022 09:41:18 -0700
Subject: [PATCH 12/23] remove new line

---
 pvlib/snow.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index 7f267d395a..d789cc1c38 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -212,7 +212,6 @@ def _townsend_effective_snow(snow_load, snow_events):
        Record of the IEEE Photovoltaic Specialists Conference.
        003231-003236. :doi:`10.1109/PVSC.2011.6186627`
        Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
-
     '''  # noqa: E501
     snow_events_no_zeros = np.maximum(snow_events, 1)
     effective_snow = 0.5 * snow_load * (1 + 1 / snow_events_no_zeros)

From e2c1f76aa5d9d2bac10e5bd4d7f733bab8bc5caa Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Tue, 7 Jun 2022 09:45:32 -0700
Subject: [PATCH 13/23] Se to effective_snow

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

diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index cb72ae58b8..14a2ed035a 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -97,7 +97,7 @@ def test_dc_loss_nrel():
     assert_series_equal(expected, actual)
 
 
-def test__townsend_Se():
+def test__townsend_effective_snow():
     snow_load = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
     snow_events = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
     expected = np.array([7.5, 7.5, 5, 0, 0, 0, 0, 0, 0, 0, 3.75, 6.66666667])

From 4783556e22d1f67e6b9251d42305699c7f32099b Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Tue, 7 Jun 2022 09:57:08 -0700
Subject: [PATCH 14/23] adding PR number to whatsnew

---
 docs/sphinx/source/whatsnew/v0.9.2.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/sphinx/source/whatsnew/v0.9.2.rst b/docs/sphinx/source/whatsnew/v0.9.2.rst
index ac72453d14..476a8d6edf 100644
--- a/docs/sphinx/source/whatsnew/v0.9.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.2.rst
@@ -19,7 +19,7 @@ Bug fixes
   rise/set/transit times for the wrong day in recent versions of ``ephem``
   (:issue:`1449`, :pull:`1448`)
 * Added Townsend-Powers monthly snow loss model: :py:func:`pvlib.snow.loss_townsend`
-  (:issue:`1246`, :pull:`1251`)
+  (:issue:`1246`, :pull:`1251`, :pull:`1468`)
 
 
 Testing

From 588115ed36656eb50629efd10cc5608847df2cd5 Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Tue, 7 Jun 2022 10:04:14 -0700
Subject: [PATCH 15/23] address stickler line too long

---
 pvlib/snow.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index d789cc1c38..304caf656f 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -280,8 +280,14 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
     snow_events_prev = np.roll(snow_events, 1)
 
     effective_snow = _townsend_effective_snow(snow_total, snow_events)
-    effective_snow_prev = _townsend_effective_snow(snow_total_prev, snow_events_prev)
-    effective_snow_weighted = 1 / 3 * effective_snow_prev + 2 / 3 * effective_snow
+    effective_snow_prev = _townsend_effective_snow(
+        snow_total_prev,
+        snow_events_prev
+    )
+    effective_snow_weighted = (
+        1 / 3 * effective_snow_prev
+        + 2 / 3 * effective_snow
+    )
 
     drop_height_clipped = np.maximum(lower_edge_drop_height, 0.01)
     gamma = (

From 36ae2640cf76399ca2cf19577e286d884cb45b41 Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Mon, 13 Jun 2022 09:53:55 -0700
Subject: [PATCH 16/23] remove links and noqa E501, and fix long lines

---
 docs/sphinx/source/whatsnew/v0.9.2.rst |  5 ++---
 pvlib/snow.py                          | 16 ++++++++--------
 2 files changed, 10 insertions(+), 11 deletions(-)

diff --git a/docs/sphinx/source/whatsnew/v0.9.2.rst b/docs/sphinx/source/whatsnew/v0.9.2.rst
index 476a8d6edf..ba74184bda 100644
--- a/docs/sphinx/source/whatsnew/v0.9.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.2.rst
@@ -8,6 +8,8 @@ Deprecations
 
 Enhancements
 ~~~~~~~~~~~~
+* Added Townsend-Powers monthly snow loss model: :py:func:`pvlib.snow.loss_townsend`
+  (:issue:`1246`, :pull:`1251`, :pull:`1468`)
 
 Bug fixes
 ~~~~~~~~~
@@ -18,9 +20,6 @@ Bug fixes
   where passing localized timezones with large UTC offsets could return
   rise/set/transit times for the wrong day in recent versions of ``ephem``
   (:issue:`1449`, :pull:`1448`)
-* Added Townsend-Powers monthly snow loss model: :py:func:`pvlib.snow.loss_townsend`
-  (:issue:`1246`, :pull:`1251`, :pull:`1468`)
-
 
 Testing
 ~~~~~~~
diff --git a/pvlib/snow.py b/pvlib/snow.py
index 304caf656f..366c0adea6 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -211,8 +211,7 @@ def _townsend_effective_snow(snow_load, snow_events):
        update from two winters of measurements in the SIERRA. Conference
        Record of the IEEE Photovoltaic Specialists Conference.
        003231-003236. :doi:`10.1109/PVSC.2011.6186627`
-       Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
-    '''  # noqa: E501
+    '''
     snow_events_no_zeros = np.maximum(snow_events, 1)
     effective_snow = 0.5 * snow_load * (1 + 1 / snow_events_no_zeros)
     return np.where(snow_events > 0, effective_snow, 0)
@@ -222,7 +221,8 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
                   temp_air, poa_global, slant_height, lower_edge_drop_height,
                   angle_of_repose=40):
     '''
-    Calculates monthly snow loss based on the Townsend monthly snow loss model [1]_.
+    Calculates monthly snow loss based on the Townsend monthly snow loss
+    model [1]_.
 
     Parameters
     ----------
@@ -242,7 +242,7 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
         Monthly average ambient temperature [C]
 
     poa_global : array-like
-        Monthly plane of array insolation [kWh/m2/month]
+        Monthly plane of array insolation [kWh/m2]
 
     slant_height : float
         Row length in the slanted plane of array dimension [in]
@@ -261,8 +261,9 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
 
     Notes
     -----
-    This model has not been validated for tracking arrays; however, for tracking
-    arrays [1]_ suggests using the maximum rotation angle in place of surface_tilt.
+    This model has not been validated for tracking arrays; however, for
+    tracking arrays [1]_ suggests using the maximum rotation angle in place
+    of surface_tilt.
 
     References
     ----------
@@ -270,8 +271,7 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
        update from two winters of measurements in the SIERRA. Conference
        Record of the IEEE Photovoltaic Specialists Conference.
        003231-003236. 10.1109/PVSC.2011.6186627.
-       Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
-    '''  # noqa: E501
+    '''
 
     C1 = 5.7e04
     C2 = 0.51

From d9dad348e7aa8583ab17dc23bb9f1efc7990ad6f Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Mon, 13 Jun 2022 12:39:35 -0700
Subject: [PATCH 17/23] converting to metric system

---
 pvlib/snow.py            | 29 +++++++++++++++--------------
 pvlib/tests/test_snow.py | 19 +++++++++++--------
 2 files changed, 26 insertions(+), 22 deletions(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index 366c0adea6..e351d85b22 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -187,23 +187,23 @@ def dc_loss_nrel(snow_coverage, num_strings):
     return np.ceil(snow_coverage * num_strings) / num_strings
 
 
-def _townsend_effective_snow(snow_load, snow_events):
+def _townsend_effective_snow(snow_total, snow_events):
     '''
-    Calculates effective snow using the total snowfall in inches received
-    each month and the number of snowfall events each month.
+    Calculates effective snow using the total snowfall received each month and
+    the number of snowfall events each month.
 
     Parameters
     ----------
-    snow_load : array-like
-        Inches of snow received each month. Referred to as S in the paper [in]
+    snow_total : array-like
+        Snow received each month. Referred to as S in [1]_ [cm]
 
     snow_events : array-like
-        Number of snowfall events each month. Referred to as N in the paper [-]
+        Number of snowfall events each month. Referred to as N in [1]_ [-]
 
     Returns
     -------
     effective_snowfall : array-like
-        Effective snowfall as defined in the Townsend model
+        Effective snowfall as defined in the Townsend model [cm]
 
     References
     ----------
@@ -213,7 +213,7 @@ def _townsend_effective_snow(snow_load, snow_events):
        003231-003236. :doi:`10.1109/PVSC.2011.6186627`
     '''
     snow_events_no_zeros = np.maximum(snow_events, 1)
-    effective_snow = 0.5 * snow_load * (1 + 1 / snow_events_no_zeros)
+    effective_snow = 0.5 * snow_total * (1 + 1 / snow_events_no_zeros)
     return np.where(snow_events > 0, effective_snow, 0)
 
 
@@ -227,10 +227,10 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
     Parameters
     ----------
     snow_total : array-like
-        Inches of snow received each month. Referred to as S in the paper [in]
+        Snow received each month. Referred to as S in [1]_ [cm]
 
     snow_events : array-like
-        Number of snowfall events each month. Referred to as N in the paper [-]
+        Number of snowfall events each month. Referred to as N in [1]_ [-]
 
     surface_tilt : float
         Tilt angle of the array [deg]
@@ -245,10 +245,10 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
         Monthly plane of array insolation [kWh/m2]
 
     slant_height : float
-        Row length in the slanted plane of array dimension [in]
+        Row length in the slanted plane of array dimension [m]
 
     lower_edge_drop_height : float
-        Drop height from array edge to ground [in]
+        Drop height from array edge to ground [m]
 
     angle_of_repose : float, default 40
         piled snow angle, assumed to stabilize at 40°, the midpoint of
@@ -288,13 +288,14 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
         1 / 3 * effective_snow_prev
         + 2 / 3 * effective_snow
     )
+    effective_snow_weighted_m = effective_snow_weighted / 100
 
     drop_height_clipped = np.maximum(lower_edge_drop_height, 0.01)
     gamma = (
         slant_height
-        * effective_snow_weighted
+        * effective_snow_weighted_m
         * cosd(surface_tilt)
-        / (drop_height_clipped**2 - effective_snow_weighted**2)
+        / (drop_height_clipped**2 - effective_snow_weighted_m**2)
         * 2
         * tand(angle_of_repose)
     )
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index 14a2ed035a..6e8f7d9d23 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -98,15 +98,18 @@ def test_dc_loss_nrel():
 
 
 def test__townsend_effective_snow():
-    snow_load = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
+    snow_total = np.array([25.4, 25.4, 12.7, 2.54, 0, 0, 0, 0, 0, 0, 12.7,
+                           25.4])
     snow_events = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
-    expected = np.array([7.5, 7.5, 5, 0, 0, 0, 0, 0, 0, 0, 3.75, 6.66666667])
-    actual = snow._townsend_effective_snow(snow_load, snow_events)
+    expected = np.array([19.05, 19.05, 12.7, 0, 0, 0, 0, 0, 0, 0, 9.525,
+                         254 / 15])
+    actual = snow._townsend_effective_snow(snow_total, snow_events)
     np.testing.assert_allclose(expected, actual, rtol=1e-07)
 
 
 def test_loss_townsend():
-    snow_total = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
+    snow_total = np.array([25.4, 25.4, 12.7, 2.54, 0, 0, 0, 0, 0, 0, 12.7,
+                           25.4])
     snow_events = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
     surface_tilt = 20
     relative_humidity = np.array([80, 80, 80, 80, 80, 80, 80, 80, 80, 80,
@@ -115,10 +118,10 @@ def test_loss_townsend():
     poa_global = np.array([350, 350, 350, 350, 350, 350, 350, 350, 350, 350,
                            350, 350])
     angle_of_repose = 40
-    slant_height = 100
-    lower_edge_drop_height = 10
-    expected = np.array([0.07696253, 0.07992262, 0.06216201, 0.01715392, 0, 0,
-                         0, 0, 0, 0, 0.02643821, 0.06068194])
+    slant_height = 2.54
+    lower_edge_drop_height = 0.254
+    expected = np.array([0.195485, 0.203003, 0.157892, 0.043571, 0, 0, 0, 0,
+                         0, 0, 0.067153, 0.154132])
     actual = snow.loss_townsend(snow_total, snow_events, surface_tilt,
                                 relative_humidity, temp_air,
                                 poa_global, slant_height,

From 1d5eb8835bef1533631e4ce2fc1a562ccdd7a934 Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Mon, 13 Jun 2022 12:53:30 -0700
Subject: [PATCH 18/23] poa_global from kWh/m2 to Wh/m2

---
 pvlib/snow.py            | 5 +++--
 pvlib/tests/test_snow.py | 4 ++--
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index e351d85b22..22ba01fda4 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -242,7 +242,7 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
         Monthly average ambient temperature [C]
 
     poa_global : array-like
-        Monthly plane of array insolation [kWh/m2]
+        Monthly plane of array insolation [Wh/m2]
 
     slant_height : float
         Row length in the slanted plane of array dimension [m]
@@ -302,6 +302,7 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
 
     ground_interference_term = 1 - C2 * np.exp(-gamma)
     temp_air_kelvin = temp_air + 273.15
+    poa_global_kWh = poa_global / 1000
     loss_percentage = (
         C1
         * effective_snow_weighted
@@ -309,7 +310,7 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
         * ground_interference_term
         * relative_humidity
         / temp_air_kelvin**2
-        / poa_global**0.67
+        / poa_global_kWh**0.67
     )
 
     return loss_percentage / 100
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index 6e8f7d9d23..5da02383c5 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -115,8 +115,8 @@ def test_loss_townsend():
     relative_humidity = np.array([80, 80, 80, 80, 80, 80, 80, 80, 80, 80,
                                   80, 80])
     temp_air = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
-    poa_global = np.array([350, 350, 350, 350, 350, 350, 350, 350, 350, 350,
-                           350, 350])
+    poa_global = np.array([350000, 350000, 350000, 350000, 350000, 350000,
+                           350000, 350000, 350000, 350000, 350000, 350000])
     angle_of_repose = 40
     slant_height = 2.54
     lower_edge_drop_height = 0.254

From bbe5bc33f8e6fff0d0e405484a7bd70863cc6c45 Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Mon, 13 Jun 2022 13:30:49 -0700
Subject: [PATCH 19/23] changing returned loss from kWh to Wh

---
 pvlib/snow.py            | 3 +--
 pvlib/tests/test_snow.py | 4 ++--
 2 files changed, 3 insertions(+), 4 deletions(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index 22ba01fda4..cd0afba2e3 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -302,7 +302,6 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
 
     ground_interference_term = 1 - C2 * np.exp(-gamma)
     temp_air_kelvin = temp_air + 273.15
-    poa_global_kWh = poa_global / 1000
     loss_percentage = (
         C1
         * effective_snow_weighted
@@ -310,7 +309,7 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
         * ground_interference_term
         * relative_humidity
         / temp_air_kelvin**2
-        / poa_global_kWh**0.67
+        / poa_global**0.67
     )
 
     return loss_percentage / 100
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index 5da02383c5..52677a854e 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -120,8 +120,8 @@ def test_loss_townsend():
     angle_of_repose = 40
     slant_height = 2.54
     lower_edge_drop_height = 0.254
-    expected = np.array([0.195485, 0.203003, 0.157892, 0.043571, 0, 0, 0, 0,
-                         0, 0, 0.067153, 0.154132])
+    expected = np.array([0.00191035, 0.00198382, 0.00154297, 0.00042579, 0, 0,
+                         0, 0, 0, 0, 0.00065624, 0.00150623])
     actual = snow.loss_townsend(snow_total, snow_events, surface_tilt,
                                 relative_humidity, temp_air,
                                 poa_global, slant_height,

From b3ab5aa977076834c4d06d95b9c59f714885b6f8 Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Mon, 13 Jun 2022 17:42:15 -0700
Subject: [PATCH 20/23] fixing capacity loss calculation units to keep correct
 C1 value

---
 pvlib/snow.py            | 6 ++++--
 pvlib/tests/test_snow.py | 4 ++--
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index cd0afba2e3..adfb9c2851 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -302,14 +302,16 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
 
     ground_interference_term = 1 - C2 * np.exp(-gamma)
     temp_air_kelvin = temp_air + 273.15
+    effective_snow_weighted_in = effective_snow_weighted / 2.54
+    poa_global_kWh = poa_global / 1000
     loss_percentage = (
         C1
-        * effective_snow_weighted
+        * effective_snow_weighted_in
         * cosd(surface_tilt)**2
         * ground_interference_term
         * relative_humidity
         / temp_air_kelvin**2
-        / poa_global**0.67
+        / poa_global_kWh**0.67
     )
 
     return loss_percentage / 100
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index 52677a854e..f770232005 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -120,8 +120,8 @@ def test_loss_townsend():
     angle_of_repose = 40
     slant_height = 2.54
     lower_edge_drop_height = 0.254
-    expected = np.array([0.00191035, 0.00198382, 0.00154297, 0.00042579, 0, 0,
-                         0, 0, 0, 0, 0.00065624, 0.00150623])
+    expected = np.array([0.07696253, 0.07992262, 0.06216201, 0.01715392, 0, 0,
+                         0, 0, 0, 0, 0.02643821, 0.06068194])
     actual = snow.loss_townsend(snow_total, snow_events, surface_tilt,
                                 relative_humidity, temp_air,
                                 poa_global, slant_height,

From b759bcfe8f000221d22ba7b1fb092b838677e94b Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Mon, 20 Jun 2022 08:50:49 -0700
Subject: [PATCH 21/23] convert relative humidity from percent to fraction

---
 pvlib/snow.py            | 3 ++-
 pvlib/tests/test_snow.py | 5 +++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index adfb9c2851..1b7a4a3249 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -301,6 +301,7 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
     )
 
     ground_interference_term = 1 - C2 * np.exp(-gamma)
+    relative_humidity_fraction = relative_humidity / 100
     temp_air_kelvin = temp_air + 273.15
     effective_snow_weighted_in = effective_snow_weighted / 2.54
     poa_global_kWh = poa_global / 1000
@@ -309,7 +310,7 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
         * effective_snow_weighted_in
         * cosd(surface_tilt)**2
         * ground_interference_term
-        * relative_humidity
+        * relative_humidity_fraction
         / temp_air_kelvin**2
         / poa_global_kWh**0.67
     )
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index f770232005..6e9782e349 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -120,8 +120,9 @@ def test_loss_townsend():
     angle_of_repose = 40
     slant_height = 2.54
     lower_edge_drop_height = 0.254
-    expected = np.array([0.07696253, 0.07992262, 0.06216201, 0.01715392, 0, 0,
-                         0, 0, 0, 0, 0.02643821, 0.06068194])
+    expected = np.array([0.0007696253, 0.0007992262, 0.0006216201,
+                         0.0001715392, 0, 0, 0, 0, 0, 0, 0.0002643821,
+                         0.0006068194])
     actual = snow.loss_townsend(snow_total, snow_events, surface_tilt,
                                 relative_humidity, temp_air,
                                 poa_global, slant_height,

From f46fc00cb67f737e0dc3fd78ef4e359f178cd30a Mon Sep 17 00:00:00 2001
From: reepoi <spacechips@rocketmail.com>
Date: Wed, 22 Jun 2022 15:39:55 -0700
Subject: [PATCH 22/23] neatening docstrings and adjusting variable names

---
 pvlib/snow.py            | 50 ++++++++++++++++++++--------------------
 pvlib/tests/test_snow.py |  4 ++--
 2 files changed, 27 insertions(+), 27 deletions(-)

diff --git a/pvlib/snow.py b/pvlib/snow.py
index 1b7a4a3249..579a812c38 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -195,22 +195,22 @@ def _townsend_effective_snow(snow_total, snow_events):
     Parameters
     ----------
     snow_total : array-like
-        Snow received each month. Referred to as S in [1]_ [cm]
+        Snow received each month. Referred to as S in [1]_. [cm]
 
     snow_events : array-like
-        Number of snowfall events each month. Referred to as N in [1]_ [-]
+        Number of snowfall events each month. Referred to as N in [1]_. [-]
 
     Returns
     -------
     effective_snowfall : array-like
-        Effective snowfall as defined in the Townsend model [cm]
+        Effective snowfall as defined in the Townsend model. [cm]
 
     References
     ----------
     .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An
-       update from two winters of measurements in the SIERRA. Conference
-       Record of the IEEE Photovoltaic Specialists Conference.
-       003231-003236. :doi:`10.1109/PVSC.2011.6186627`
+       update from two winters of measurements in the SIERRA. 37th IEEE
+       Photovoltaic Specialists Conference, Seattle, WA, USA.
+       :doi:`10.1109/PVSC.2011.6186627`
     '''
     snow_events_no_zeros = np.maximum(snow_events, 1)
     effective_snow = 0.5 * snow_total * (1 + 1 / snow_events_no_zeros)
@@ -218,7 +218,7 @@ def _townsend_effective_snow(snow_total, snow_events):
 
 
 def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
-                  temp_air, poa_global, slant_height, lower_edge_drop_height,
+                  temp_air, poa_global, slant_height, lower_edge_height,
                   angle_of_repose=40):
     '''
     Calculates monthly snow loss based on the Townsend monthly snow loss
@@ -227,50 +227,50 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
     Parameters
     ----------
     snow_total : array-like
-        Snow received each month. Referred to as S in [1]_ [cm]
+        Snow received each month. Referred to as S in [1]_. [cm]
 
     snow_events : array-like
-        Number of snowfall events each month. Referred to as N in [1]_ [-]
+        Number of snowfall events each month. Referred to as N in [1]_. [-]
 
     surface_tilt : float
-        Tilt angle of the array [deg]
+        Tilt angle of the array. [deg]
 
     relative_humidity : array-like
-        Monthly average relative humidity [%]
+        Monthly average relative humidity. [%]
 
     temp_air : array-like
-        Monthly average ambient temperature [C]
+        Monthly average ambient temperature. [C]
 
     poa_global : array-like
-        Monthly plane of array insolation [Wh/m2]
+        Monthly plane of array insolation. [Wh/m2]
 
     slant_height : float
-        Row length in the slanted plane of array dimension [m]
+        Row length in the slanted plane of array dimension. [m]
 
-    lower_edge_drop_height : float
-        Drop height from array edge to ground [m]
+    lower_edge_height : float
+        Distance from array lower edge to the ground. [m]
 
     angle_of_repose : float, default 40
-        piled snow angle, assumed to stabilize at 40°, the midpoint of
-        25°-55° avalanching slope angles [deg]
+        Piled snow angle, assumed to stabilize at 40°, the midpoint of
+        25°-55° avalanching slope angles. [deg]
 
     Returns
     -------
     loss : array-like
-        Monthly average DC capacity loss fraction due to snow coverage
+        Monthly average DC capacity loss fraction due to snow coverage.
 
     Notes
     -----
     This model has not been validated for tracking arrays; however, for
     tracking arrays [1]_ suggests using the maximum rotation angle in place
-    of surface_tilt.
+    of ``surface_tilt``.
 
     References
     ----------
     .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An
-       update from two winters of measurements in the SIERRA. Conference
-       Record of the IEEE Photovoltaic Specialists Conference.
-       003231-003236. 10.1109/PVSC.2011.6186627.
+       update from two winters of measurements in the SIERRA. 37th IEEE
+       Photovoltaic Specialists Conference, Seattle, WA, USA.
+       :doi:`10.1109/PVSC.2011.6186627`
     '''
 
     C1 = 5.7e04
@@ -290,12 +290,12 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
     )
     effective_snow_weighted_m = effective_snow_weighted / 100
 
-    drop_height_clipped = np.maximum(lower_edge_drop_height, 0.01)
+    lower_edge_height_clipped = np.maximum(lower_edge_height, 0.01)
     gamma = (
         slant_height
         * effective_snow_weighted_m
         * cosd(surface_tilt)
-        / (drop_height_clipped**2 - effective_snow_weighted_m**2)
+        / (lower_edge_height_clipped**2 - effective_snow_weighted_m**2)
         * 2
         * tand(angle_of_repose)
     )
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index 6e9782e349..fd728fe80a 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -119,12 +119,12 @@ def test_loss_townsend():
                            350000, 350000, 350000, 350000, 350000, 350000])
     angle_of_repose = 40
     slant_height = 2.54
-    lower_edge_drop_height = 0.254
+    lower_edge_height = 0.254
     expected = np.array([0.0007696253, 0.0007992262, 0.0006216201,
                          0.0001715392, 0, 0, 0, 0, 0, 0, 0.0002643821,
                          0.0006068194])
     actual = snow.loss_townsend(snow_total, snow_events, surface_tilt,
                                 relative_humidity, temp_air,
                                 poa_global, slant_height,
-                                lower_edge_drop_height, angle_of_repose)
+                                lower_edge_height, angle_of_repose)
     np.testing.assert_allclose(expected, actual, rtol=1e-05)

From 0d1542f4e09680c19af8d0989b39c3af1dba49c3 Mon Sep 17 00:00:00 2001
From: Taos Transue <spacechips@rocketmail.com>
Date: Sun, 11 Sep 2022 08:59:25 -0700
Subject: [PATCH 23/23] changing eqn 3 percentage loss to fractional loss and
 adding comment explanation

---
 docs/sphinx/source/whatsnew/v0.9.2.rst |  4 ----
 docs/sphinx/source/whatsnew/v0.9.3.rst |  5 +++++
 pvlib/snow.py                          | 10 ++++++++--
 pvlib/tests/test_snow.py               |  5 ++---
 4 files changed, 15 insertions(+), 9 deletions(-)

diff --git a/docs/sphinx/source/whatsnew/v0.9.2.rst b/docs/sphinx/source/whatsnew/v0.9.2.rst
index 6da4ee555e..d8221e48db 100644
--- a/docs/sphinx/source/whatsnew/v0.9.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.2.rst
@@ -5,8 +5,6 @@ v0.9.2 (August 19, 2022)
 
 Enhancements
 ~~~~~~~~~~~~
-* Added Townsend-Powers monthly snow loss model: :py:func:`pvlib.snow.loss_townsend`
-  (:issue:`1246`, :pull:`1251`, :pull:`1468`)
 * albedo can now be provided as a column in the `weather` DataFrame input to
   :py:meth:`pvlib.modelchain.ModelChain.run_model`. (:issue:`1387`, :pull:`1478`)
 * albedo is now available as an input to :py:meth:`pvlib.pvsystem.PVSystem.get_irradiance`
@@ -74,8 +72,6 @@ Contributors
 * Prajwal Borkar (:ghuser:`PrajwalBorkar`) 
 * Cliff Hansen (:ghuser:`cwhanse`)
 * Kevin Anderson (:ghuser:`kanderso-nrel`)
-* Abhishek Parikh (:ghuser:`abhisheksparikh`)
-* Taos Transue (:ghuser:`reepoi`)
 * Cliff Hansen (:ghuser:`cwhanse`)
 * Jules Chéron (:ghuser:`jules-ch`)
 * Kurt Rhee (:ghuser:`kurt-rhee`)
diff --git a/docs/sphinx/source/whatsnew/v0.9.3.rst b/docs/sphinx/source/whatsnew/v0.9.3.rst
index 2f3c643d05..3bd6c45e4f 100644
--- a/docs/sphinx/source/whatsnew/v0.9.3.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.3.rst
@@ -12,6 +12,9 @@ Enhancements
 * Low resolution altitude lookup map
   :py:func:`~pvlib.location.lookup_altitude`
   (:issue:`1516`, :pull:`1518`)
+* Added Townsend-Powers monthly snow loss model:
+  :py:func:`pvlib.snow.loss_townsend`
+  (:issue:`1246`, :pull:`1251`, :pull:`1468`)
 
 Bug fixes
 ~~~~~~~~~
@@ -37,3 +40,5 @@ Contributors
 ~~~~~~~~~~~~
 * João Guilherme (:ghuser:`joaoguilhermeS`)
 * Nicolas Martinez (:ghuser:`nicomt`)
+* Abhishek Parikh (:ghuser:`abhisheksparikh`)
+* Taos Transue (:ghuser:`reepoi`)
diff --git a/pvlib/snow.py b/pvlib/snow.py
index 579a812c38..c06b317a2f 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -305,7 +305,13 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
     temp_air_kelvin = temp_air + 273.15
     effective_snow_weighted_in = effective_snow_weighted / 2.54
     poa_global_kWh = poa_global / 1000
-    loss_percentage = (
+
+    # Calculate Eqn. 3 in the reference.
+    # Although the reference says Eqn. 3 calculates percentage loss, the y-axis
+    # of Figure 7 indicates Eqn. 3 calculates fractional loss. Since the slope
+    # of the line in Figure 7 is the same as C1 in Eqn. 3, it is assumed that
+    # Eqn. 3 calculates fractional loss.
+    loss_fraction = (
         C1
         * effective_snow_weighted_in
         * cosd(surface_tilt)**2
@@ -315,4 +321,4 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
         / poa_global_kWh**0.67
     )
 
-    return loss_percentage / 100
+    return np.clip(loss_fraction, 0, 1)
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index fd728fe80a..11f0b21e83 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -120,9 +120,8 @@ def test_loss_townsend():
     angle_of_repose = 40
     slant_height = 2.54
     lower_edge_height = 0.254
-    expected = np.array([0.0007696253, 0.0007992262, 0.0006216201,
-                         0.0001715392, 0, 0, 0, 0, 0, 0, 0.0002643821,
-                         0.0006068194])
+    expected = np.array([0.07696253, 0.07992262, 0.06216201, 0.01715392, 0, 0,
+                         0, 0, 0, 0, 0.02643821, 0.06068194])
     actual = snow.loss_townsend(snow_total, snow_events, surface_tilt,
                                 relative_humidity, temp_air,
                                 poa_global, slant_height,