Add pw and spectral correction code from pvlib matlab 1.3 (#115)

* add first solar spec correction. needs tests

* add calc_pw. needs tests

* add lower to coefs input

* clean up rebase. shorten long lines

* rename, clean docs, add pw test

* update fs coeffs, add fs tests

* doc fixes

* add 040 whatsnew

Author: wholmgren

docs/sphinx/source/whatsnew.rst

@@ -6,6 +6,7 @@ What's New
 
 These are new features and improvements of note in each release.
 
+.. include:: whatsnew/v0.4.0.txt
 .. include:: whatsnew/v0.3.3.txt
 .. include:: whatsnew/v0.3.2.txt
 .. include:: whatsnew/v0.3.1.txt

docs/sphinx/source/whatsnew/v0.4.0.txt

@@ -0,0 +1,41 @@
+.. _whatsnew_0400:
+
+v0.4.0 (July xx, 2016)
+-----------------------
+
+This is a major release from 0.3.3.
+We recommend that all users upgrade to this version after reviewing
+the API changes.
+
+
+API Changes
+~~~~~~~~~~~
+
+
+
+Enhancements
+~~~~~~~~~~~~
+
+* Adds the First Solar spectral correction model. (:issue:`115`)
+* Adds the Gueymard 1994 integrated precipitable water model. (:issue:`115`)
+
+
+Bug fixes
+~~~~~~~~~
+
+
+
+Documentation
+~~~~~~~~~~~~~
+
+
+
+Other
+~~~~~
+
+
+
+Code Contributors
+~~~~~~~~~~~~~~~~~
+
+* Will Holmgren

pvlib/atmosphere.py

@@ -242,3 +242,199 @@ def relativeairmass(zenith, model='kastenyoung1989'):
         am = np.nan if z > 90 else am
 
     return am
+
+
+def gueymard94_pw(temp_air, relative_humidity):
+    r"""
+    Calculates precipitable water (cm) from ambient air temperature (C)
+    and relatively humidity (%) using an empirical model. The
+    accuracy of this method is approximately 20% for moderate PW (1-3
+    cm) and less accurate otherwise.
+
+    The model was developed by expanding Eq. 1 in [2]_:
+
+    .. math::
+
+           w = 0.1 H_v \rho_v
+
+    using Eq. 2 in [2]_
+
+    .. math::
+
+           \rho_v = 216.7 R_H e_s /T
+
+    :math:`H_v` is the apparant water vapor scale height (km). The
+    expression for :math:`H_v` is Eq. 4 in [2]_:
+
+    .. math::
+
+           H_v = 0.4976 + 1.5265*T/273.15 + \exp(13.6897*T/273.15 - 14.9188*(T/273.15)^3)
+
+    :math:`\rho_v` is the surface water vapor density (g/m^3). In the
+    expression :math:`\rho_v`, :math:`e_s` is the saturation water vapor
+    pressure (millibar). The
+    expression for :math:`e_s` is Eq. 1 in [3]_
+
+    .. math::
+
+          e_s = \exp(22.330 - 49.140*(100/T) - 10.922*(100/T)^2 - 0.39015*T/100)
+
+    Parameters
+    ----------
+    temp_air : array-like
+        ambient air temperature at the surface (C)
+    relative_humidity : array-like
+        relative humidity at the surface (%)
+
+    Returns
+    -------
+    pw : array-like
+        precipitable water (cm)
+
+    References
+    ----------
+    .. [1] W. M. Keogh and A. W. Blakers, Accurate Measurement, Using Natural
+       Sunlight, of Silicon Solar Cells, Prog. in Photovoltaics: Res.
+       and Appl. 2004, vol 12, pp. 1-19 (DOI: 10.1002/pip.517)
+
+    .. [2] C. Gueymard, Analysis of Monthly Average Atmospheric Precipitable
+       Water and Turbidity in Canada and Northern United States,
+       Solar Energy vol 53(1), pp. 57-71, 1994.
+
+    .. [3] C. Gueymard, Assessment of the Accuracy and Computing Speed of
+       simplified saturation vapor equations using a new reference
+       dataset, J. of Applied Meteorology 1993, vol. 32(7), pp.
+       1294-1300.
+    """
+
+    T = temp_air + 273.15  # Convert to Kelvin
+    RH = relative_humidity
+
+    theta = T / 273.15
+
+    # Eq. 1 from Keogh and Blakers
+    pw = (
+        0.1 *
+        (0.4976 + 1.5265*theta + np.exp(13.6897*theta - 14.9188*(theta)**3)) *
+        (216.7*RH/(100*T)*np.exp(22.330 - 49.140*(100/T) -
+         10.922*(100/T)**2 - 0.39015*T/100)))
+
+    pw = np.maximum(pw, 0.1)
+
+    return pw
+
+
+def first_solar_spectral_correction(pw, airmass_absolute, module_type=None,
+                                    coefficients=None):
+    r"""
+    Spectral mismatch modifier based on precipitable water and absolute
+    (pressure corrected) airmass.
+
+    Estimates a spectral mismatch modifier M representing the effect on
+    module short circuit current of variation in the spectral
+    irradiance. M is estimated from absolute (pressure currected) air
+    mass, AMa, and precipitable water, Pwat, using the following
+    function:
+
+    .. math::
+
+        M = c_1 + c_2*AMa  + c_3*Pwat  + c_4*AMa^.5
+            + c_5*Pwat^.5 + c_6*AMa/Pwat
+
+    Default coefficients are determined for several cell types with
+    known quantum efficiency curves, by using the Simple Model of the
+    Atmospheric Radiative Transfer of Sunshine (SMARTS) [1]_. Using
+    SMARTS, spectrums are simulated with all combinations of AMa and
+    Pwat where:
+
+       * 0.5 cm <= Pwat <= 5 cm
+       * 0.8 <= AMa <= 4.75 (Pressure of 800 mbar and 1.01 <= AM <= 6)
+       * Spectral range is limited to that of CMP11 (280 nm to 2800 nm)
+       * spectrum simulated on a plane normal to the sun
+       * All other parameters fixed at G173 standard
+
+    From these simulated spectra, M is calculated using the known
+    quantum efficiency curves. Multiple linear regression is then
+    applied to fit Eq. 1 to determine the coefficients for each module.
+
+    Based on the PVLIB Matlab function ``pvl_FSspeccorr`` by Mitchell
+    Lee and Alex Panchula, at First Solar, 2015.
+
+    Parameters
+    ----------
+    pw : array-like
+        atmospheric precipitable water (cm).
+
+    airmass_absolute :
+        absolute (pressure corrected) airmass.
+
+    module_type : None or string
+        a string specifying a cell type. Can be lower or upper case
+        letters. Admits values of 'cdte', 'monosi', 'xsi', 'multisi',
+        'polysi'. If provided, this input selects coefficients for the
+        following default modules:
+
+            * 'cdte' - First Solar Series 4-2 CdTe modules.
+            * 'monosi', 'xsi' - First Solar TetraSun modules.
+            * 'multisi', 'polysi' - multi-crystalline silicon modules.
+
+        The module used to calculate the spectral correction
+        coefficients corresponds to the Mult-crystalline silicon
+        Manufacturer 2 Model C from [2]_.
+
+    coefficients : array-like
+        allows for entry of user defined spectral correction
+        coefficients. Coefficients must be of length 6. Derivation of
+        coefficients requires use of SMARTS and PV module quantum
+        efficiency curve. Useful for modeling PV module types which are
+        not included as defaults, or to fine tune the spectral
+        correction to a particular mono-Si, multi-Si, or CdTe PV module.
+        Note that the parameters for modules with very similar QE should
+        be similar, in most cases limiting the need for module specific
+        coefficients.
+
+    Returns
+    -------
+    modifier: array-like
+        spectral mismatch factor (unitless) which is can be multiplied
+        with broadband irradiance reaching a module's cells to estimate
+        effective irradiance, i.e., the irradiance that is converted to
+        electrical current.
+
+    References
+    ----------
+    .. [1] Gueymard, Christian. SMARTS2: a simple model of the atmospheric
+       radiative transfer of sunshine: algorithms and performance
+       assessment. Cocoa, FL: Florida Solar Energy Center, 1995.
+
+    .. [2] Marion, William F., et al. User's Manual for Data for Validating
+       Models for PV Module Performance. National Renewable Energy
+       Laboratory, 2014. http://www.nrel.gov/docs/fy14osti/61610.pdf
+    """
+
+    _coefficients = {}
+    _coefficients['cdte'] = (
+        0.87102, -0.040543, -0.00929202, 0.10052, 0.073062, -0.0034187)
+    _coefficients['monosi'] = (
+        0.86588, -0.021637, -0.0030218, 0.12081, 0.017514, -0.0012610)
+    _coefficients['xsi'] = _coefficients['monosi']
+    _coefficients['polysi'] = (
+        0.84674, -0.028568, -0.0051832, 0.13669, 0.029234, -0.0014207)
+    _coefficients['multisi'] = _coefficients['polysi']
+
+    if module_type is not None and coefficients is None:
+        coefficients = _coefficients[module_type.lower()]
+    elif module_type is None and coefficients is not None:
+        pass
+    else:
+        raise TypeError('ambiguous input, must supply only 1 of ' +
+                        'module_type and coefficients')
+
+    # Evaluate Spectral Shift
+    coeff = coefficients
+    AMa = airmass_absolute
+    modifier = (
+        coeff[0] + coeff[1]*AMa  + coeff[2]*pw  + coeff[3]*np.sqrt(AMa) +
+        + coeff[4]*np.sqrt(pw) + coeff[5]*AMa/pw)
+
+    return modifier

pvlib/test/test_atmosphere.py

@@ -1,13 +1,12 @@
-import logging
-pvl_logger = logging.getLogger('pvlib')
-
 import datetime
+import itertools
 
 import numpy as np
 import pandas as pd
 
 from nose.tools import raises
 from nose.tools import assert_almost_equals
+from numpy.testing import assert_allclose
 
 from pvlib.location import Location
 from pvlib import solarposition
@@ -65,3 +64,56 @@ def test_absoluteairmass_numeric():
 def test_absoluteairmass_nan():
     np.testing.assert_equal(np.nan, atmosphere.absoluteairmass(np.nan))
 
+
+def test_gueymard94_pw():
+    temp_air = np.array([0, 20, 40])
+    relative_humidity = np.array([0, 30, 100])
+    temps_humids = np.array(
+        list(itertools.product(temp_air, relative_humidity)))
+    pws = atmosphere.gueymard94_pw(temps_humids[:, 0], temps_humids[:, 1])
+
+    expected = np.array(
+        [  0.1       ,   0.33702061,   1.12340202,   0.1       ,
+         1.12040963,   3.73469877,   0.1       ,   3.44859767,  11.49532557])
+
+    assert_allclose(pws, expected, atol=0.01)
+
+
+def test_first_solar_spectral_correction():
+    ams = np.array([1, 3, 5])
+    pws = np.array([1, 3, 5])
+    ams, pws = np.meshgrid(ams, pws)
+
+    expect = {}
+    expect['cdte'] = np.array(
+        [[ 0.99134828,  0.97701063,  0.93975103],
+         [ 1.02852847,  1.01874908,  0.98604776],
+         [ 1.04722476,  1.03835703,  1.00656735]])
+    expect['monosi'] = np.array(
+        [[ 0.9782842 ,  1.02092726,  1.03602157],
+         [ 0.9859024 ,  1.0302268 ,  1.04700244],
+         [ 0.98885429,  1.03351495,  1.05062687]])
+    expect['polysi'] = np.array(
+        [[ 0.9774921 ,  1.01757872,  1.02649543],
+         [ 0.98947361,  1.0314545 ,  1.04226547],
+         [ 0.99403107,  1.03639082,  1.04758064]])
+
+    def run_fs_test(module_type):
+        out = atmosphere.first_solar_spectral_correction(pws, ams, module_type)
+        assert_allclose(out, expect[module_type], atol=0.001)
+
+    for module_type in expect.keys():
+        yield run_fs_test, module_type
+
+
+def test_first_solar_spectral_correction_supplied():
+    # use the cdte coeffs
+    coeffs = (0.87102, -0.040543, -0.00929202, 0.10052, 0.073062, -0.0034187)
+    out = atmosphere.first_solar_spectral_correction(1, 1, coefficients=coeffs)
+    expected = 0.99134828
+    assert_allclose(out, expected, atol=1e-3)
+
+
+@raises(TypeError)
+def test_first_solar_spectral_correction_ambiguous():
+    atmosphere.first_solar_spectral_correction(1, 1)