diff --git a/docs/sphinx/source/whatsnew/v0.7.0.rst b/docs/sphinx/source/whatsnew/v0.7.0.rst
index a5251d9047..505bd00825 100644
--- a/docs/sphinx/source/whatsnew/v0.7.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.7.0.rst
@@ -10,10 +10,16 @@ compatibility notes.
 **Python 2.7 support ended on June 1, 2019.** (:issue:`501`)
 **Minimum numpy version is now 1.10.4. Minimum pandas version is now 0.18.1.**
 
+Enhancements
+~~~~~~~~~~~~
+* Created two new incidence angle modifier functions: :py:func:`pvlib.pvsystem.iam_martin_ruiz`
+  and :py:func:`pvlib.pvsystem.iam_interp`. (:issue:`751`)
+
 Bug fixes
 ~~~~~~~~~
 * Fix handling of keyword arguments in `forecasts.get_processed_data`.
   (:issue:`745`)
+* Fix output as Series feature in :py:func:`pvlib.pvsystem.ashraeiam`.
 
 Testing
 ~~~~~~~
@@ -26,4 +32,5 @@ Contributors
 * Mark Campanellli (:ghuser:`markcampanelli`)
 * Will Holmgren (:ghuser:`wholmgren`)
 * Oscar Dowson (:ghuser:`odow`)
+* Anton Driesse (:ghuser:`adriesse`)
 * Alexander Morgan (:ghuser:`alexandermorgan`)
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 35701bd3a2..9000b9f611 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -17,7 +17,7 @@
 import pandas as pd
 
 from pvlib import atmosphere, irradiance, tools, singlediode as _singlediode
-from pvlib.tools import _build_kwargs
+from pvlib.tools import _build_kwargs, cosd
 from pvlib.location import Location
 
 
@@ -946,7 +946,7 @@ def ashraeiam(aoi, b=0.05):
     iam = np.where(aoi_gte_90, 0, iam)
     iam = np.maximum(0, iam)
 
-    if isinstance(iam, pd.Series):
+    if isinstance(aoi, pd.Series):
         iam = pd.Series(iam, index=aoi.index)
 
     return iam
@@ -1064,6 +1064,161 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):
     return iam
 
 
+def iam_martin_ruiz(aoi, a_r=0.16):
+    '''
+    Determine the incidence angle modifier (iam) using the Martin
+    and Ruiz incident angle model.
+
+    Parameters
+    ----------
+    aoi : numeric, degrees
+        The angle of incidence between the module normal vector and the
+        sun-beam vector in degrees. Theta must be a numeric scalar or vector.
+        iam is 0 where |aoi| > 90.
+
+    a_r : numeric
+        The angular losses coefficient described in equation 3 of [1].
+        This is an empirical dimensionless parameter. Values of a_r are
+        generally on the order of 0.08 to 0.25 for flat-plate PV modules.
+        a_r must be a positive numeric scalar or vector (same length as aoi).
+
+    Returns
+    -------
+    iam : numeric
+        The incident angle modifier(s)
+
+    Notes
+    -----
+    iam_martin_ruiz calculates the incidence angle modifier (iam)
+    as described by Martin and Ruiz in [1]. The information
+    required is the incident angle (aoi) and the angular losses
+    coefficient (a_r). Please note that [1] has a corrigendum which makes
+    the document much simpler to understand.
+
+    The incident angle modifier is defined as
+    [1-exp(-cos(aoi/ar))] / [1-exp(-1/ar)], which is
+    presented as AL(alpha) = 1 - IAM in equation 4 of [1]. Thus IAM is
+    equal to 1 at aoi = 0, and equal to 0 at aoi = 90.  This equation is only
+    valid for -90 <= aoi <= 90, therefore iam must be constrained to 0.0
+    beyond this range.
+
+    References
+    ----------
+    [1] N. Martin and J. M. Ruiz, "Calculation of the PV modules angular
+    losses under field conditions by means of an analytical model", Solar
+    Energy Materials & Solar Cells, vol. 70, pp. 25-38, 2001.
+
+    [2] N. Martin and J. M. Ruiz, "Corrigendum to 'Calculation of the PV
+    modules angular losses under field conditions by means of an
+    analytical model'", Solar Energy Materials & Solar Cells, vol. 110,
+    pp. 154, 2013.
+
+    See Also
+    --------
+    physicaliam
+    ashraeiam
+    iam_interp
+    '''
+    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019
+
+    aoi_input = aoi
+
+    aoi = np.asanyarray(aoi)
+    a_r = np.asanyarray(a_r)
+
+    if np.any(np.less_equal(a_r, 0)):
+        raise RuntimeError("The parameter 'a_r' cannot be zero or negative.")
+
+    with np.errstate(invalid='ignore'):
+        iam = (1 - np.exp(-cosd(aoi) / a_r)) / (1 - np.exp(-1 / a_r))
+        iam = np.where(np.abs(aoi) >= 90.0, 0.0, iam)
+
+    if isinstance(aoi_input, pd.Series):
+        iam = pd.Series(iam, index=aoi_input.index)
+
+    return iam
+
+
+def iam_interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
+    '''
+    Determine the incidence angle modifier (iam) by interpolating a set of
+    reference values, which are usually measured values.
+
+    Parameters
+    ----------
+    aoi : numeric, degrees
+        The angle of incidence between the module normal vector and the
+        sun-beam vector in degrees.
+
+    theta_ref : numeric, degrees
+        Vector of angles at which the iam is known.
+
+    iam_ref : numeric, unitless
+        iam values for each angle in theta_ref.
+
+    method : str, default 'linear'
+        Specifies the interpolation method.
+        Useful options are: 'linear', 'quadratic','cubic'.
+        See scipy.interpolate.interp1d for more options.
+
+    normalize : boolean
+        When true, the interpolated values are divided by the interpolated
+        value at zero degrees.  This ensures that the iam at normal
+        incidence is equal to 1.0.
+
+    Returns
+    -------
+    iam : numeric
+        The incident angle modifier(s)
+
+    Notes:
+    ------
+    theta_ref must have two or more points and may span any range of angles.
+    Typically there will be a dozen or more points in the range 0-90 degrees.
+    iam beyond the range of theta_ref are extrapolated, but constrained to be
+    non-negative.
+
+    The sign of aoi is ignored; only the magnitude is used.
+
+    See Also
+    --------
+    physicaliam
+    ashraeiam
+    iam_martin_ruiz
+    '''
+    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019
+
+    from scipy.interpolate import interp1d
+
+    # Scipy doesn't give the clearest feedback, so check number of points here.
+    MIN_REF_VALS = {'linear': 2, 'quadratic': 3, 'cubic': 4, 1: 2, 2: 3, 3: 4}
+
+    if len(theta_ref) < MIN_REF_VALS.get(method, 2):
+        raise ValueError("Too few reference points defined "
+                         "for interpolation method '%s'." % method)
+
+    if np.any(np.less(iam_ref, 0)):
+        raise ValueError("Negative value(s) found in 'iam_ref'. "
+                         "This is not physically possible.")
+
+    interpolator = interp1d(theta_ref, iam_ref, kind=method,
+                            fill_value='extrapolate')
+    aoi_input = aoi
+
+    aoi = np.asanyarray(aoi)
+    aoi = np.abs(aoi)
+    iam = interpolator(aoi)
+    iam = np.clip(iam, 0, None)
+
+    if normalize:
+        iam /= interpolator(0)
+
+    if isinstance(aoi_input, pd.Series):
+        iam = pd.Series(iam, index=aoi_input.index)
+
+    return iam
+
+
 def calcparams_desoto(effective_irradiance, temp_cell,
                       alpha_sc, a_ref, I_L_ref, I_o_ref, R_sh_ref, R_s,
                       EgRef=1.121, dEgdT=-0.0002677,
diff --git a/pvlib/test/test_pvsystem.py b/pvlib/test/test_pvsystem.py
index 17f70e9d12..7de8930f73 100644
--- a/pvlib/test/test_pvsystem.py
+++ b/pvlib/test/test_pvsystem.py
@@ -148,6 +148,88 @@ def test_PVSystem_physicaliam(mocker):
     assert iam < 1.
 
 
+def test_iam_martin_ruiz():
+
+    aoi = 45.
+    a_r = 0.16
+    expected = 0.98986965
+
+    # will fail of default values change
+    iam = pvsystem.iam_martin_ruiz(aoi)
+    assert_allclose(iam, expected)
+    # will fail of parameter names change
+    iam = pvsystem.iam_martin_ruiz(aoi=aoi, a_r=a_r)
+    assert_allclose(iam, expected)
+
+    a_r = 0.18
+    aoi = [-100, -60, 0, 60, 100, np.nan, np.inf]
+    expected = [0.0, 0.9414631, 1.0, 0.9414631, 0.0, np.nan, 0.0]
+
+    # check out of range of inputs as list
+    iam = pvsystem.iam_martin_ruiz(aoi, a_r)
+    assert_allclose(iam, expected, equal_nan=True)
+
+    # check out of range of inputs as array
+    iam = pvsystem.iam_martin_ruiz(np.array(aoi), a_r)
+    assert_allclose(iam, expected, equal_nan=True)
+
+    # check out of range of inputs as Series
+    aoi = pd.Series(aoi)
+    expected = pd.Series(expected)
+    iam = pvsystem.iam_martin_ruiz(aoi, a_r)
+    assert_series_equal(iam, expected)
+
+    # check exception clause
+    with pytest.raises(RuntimeError):
+        pvsystem.iam_martin_ruiz(0.0, a_r=0.0)
+
+
+@requires_scipy
+def test_iam_interp():
+
+    aoi_meas = [0.0, 45.0, 65.0, 75.0]
+    iam_meas = [1.0,  0.9,  0.8,  0.6]
+
+    # simple default linear method
+    aoi = 55.0
+    expected = 0.85
+    iam = pvsystem.iam_interp(aoi, aoi_meas, iam_meas)
+    assert_allclose(iam, expected)
+
+    # simple non-default method
+    aoi = 55.0
+    expected = 0.8878062
+    iam = pvsystem.iam_interp(aoi, aoi_meas, iam_meas, method='cubic')
+    assert_allclose(iam, expected)
+
+    # check with all reference values
+    aoi = aoi_meas
+    expected = iam_meas
+    iam = pvsystem.iam_interp(aoi, aoi_meas, iam_meas)
+    assert_allclose(iam, expected)
+
+    # check normalization and Series
+    aoi = pd.Series(aoi)
+    expected = pd.Series(expected)
+    iam_mult = np.multiply(0.9, iam_meas)
+    iam = pvsystem.iam_interp(aoi, aoi_meas, iam_mult, normalize=True)
+    assert_series_equal(iam, expected)
+
+    # check beyond reference values
+    aoi = [-45, 0, 45, 85, 90, 95, 100, 105, 110]
+    expected = [0.9, 1.0, 0.9, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0]
+    iam = pvsystem.iam_interp(aoi, aoi_meas, iam_meas)
+    assert_allclose(iam, expected)
+
+    # check exception clause
+    with pytest.raises(ValueError):
+        pvsystem.iam_interp(0.0, [0], [1])
+
+    # check exception clause
+    with pytest.raises(ValueError):
+        pvsystem.iam_interp(0.0, [0, 90], [1, -1])
+
+
 # if this completes successfully we'll be able to do more tests below.
 @pytest.fixture(scope="session")
 def sam_data():