diff --git a/docs/sphinx/source/reference/pv_modeling/iam.rst b/docs/sphinx/source/reference/pv_modeling/iam.rst
index 1871f9b4a2..5df44af2a2 100644
--- a/docs/sphinx/source/reference/pv_modeling/iam.rst
+++ b/docs/sphinx/source/reference/pv_modeling/iam.rst
@@ -17,3 +17,4 @@ Incident angle modifiers
    iam.marion_integrate
    iam.schlick
    iam.schlick_diffuse
+   iam.fresnel_ar
diff --git a/pvlib/iam.py b/pvlib/iam.py
index dfad91b6ff..53fdc98f40 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -19,6 +19,7 @@
     'ashrae': {'b'},
     'physical': {'n', 'K', 'L'},
     'martin_ruiz': {'a_r'},
+    'fresnel_ar': {'n_ar', 'n_air', 'n_glass'},
     'sapm': {'B0', 'B1', 'B2', 'B3', 'B4', 'B5'},
     'interp': set()
 }
@@ -750,6 +751,61 @@ def marion_integrate(function, surface_tilt, region, num=None):
     return Fd
 
 
+def fresnel_ar(aoi, n_ar=1.2, n_air=1.0, n_glass=1.56):
+    """Determine the incidence angle modifier using the Fresnel equations
+    for a surface with an anti-reflective coating.
+
+    Parameters
+    ----------
+    aoi : numeric
+        Angle of incidence between the module normal vector and the sun-beam
+        vector. [degrees]
+    n_ar : numeric, default 1.2
+        Index of refraction of the anti-reflective coating.
+    n_air : numeric, default 1.0
+        Index of refraction of the incident medium.
+    n_glass : numeric, default 1.56
+        Index of refraction of the glass cover.
+
+    Returns
+    -------
+    iam : numeric
+        The incident angle modifier.
+    """
+
+    aoi_input = aoi
+    aoi = np.asanyarray(aoi, dtype=float)
+
+    theta_ar = asind(n_air / n_ar * sind(aoi))
+    theta_glass = asind(n_ar / n_glass * sind(theta_ar))
+
+    rs1 = ((n_air * cosd(aoi) - n_ar * cosd(theta_ar)) /
+            (n_air * cosd(aoi) + n_ar * cosd(theta_ar))) ** 2
+    rp1 = ((n_air * cosd(theta_ar) - n_ar * cosd(aoi)) /
+            (n_air * cosd(theta_ar) + n_ar * cosd(aoi))) ** 2
+
+    rs2 = ((n_ar * cosd(theta_ar) - n_glass * cosd(theta_glass)) /
+            (n_ar * cosd(theta_ar) + n_glass * cosd(theta_glass))) ** 2
+    rp2 = ((n_ar * cosd(theta_glass) - n_glass * cosd(theta_ar)) /
+            (n_ar * cosd(theta_glass) + n_glass * cosd(theta_ar))) ** 2
+
+    rs = 1 - (1 - rs1) * (1 - rs2)
+    rp = 1 - (1 - rp1) * (1 - rp2)
+    refl = (rs + rp) / 2
+
+    r0_1 = ((n_air - n_ar) / (n_air + n_ar)) ** 2
+    r0_2 = ((n_ar - n_glass) / (n_ar + n_glass)) ** 2
+    r0 = 1 - (1 - r0_1) * (1 - r0_2)
+
+    iam = (1 - refl) / (1 - r0)
+    iam = np.where(np.abs(aoi) >= 90, 0.0, iam)
+
+    if isinstance(aoi_input, pd.Series):
+        iam = pd.Series(iam, index=aoi_input.index)
+
+    return iam
+
+
 def schlick(aoi):
     """
     Determine incidence angle modifier (IAM) for direct irradiance using the
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 48371ca961..851f59ec61 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -23,6 +23,37 @@
 from pvlib.tools import _build_kwargs, _build_args
 
 
+def _pvsyst_Rsh(effective_irradiance, R_sh_ref, R_sh_0, R_sh_exp, e0=1000):
+    """Simplified PVsyst shunt resistance model."""
+    effective_irradiance = np.asarray(effective_irradiance)
+    return R_sh_ref / (1 + (R_sh_ref / R_sh_0) *
+                       (effective_irradiance / e0) ** R_sh_exp)
+
+
+def _pvsyst_IL(effective_irradiance, alpha_sc, temp_cell, e0=1000, t0=25):
+    """Simplified PVsyst photocurrent model."""
+    return effective_irradiance / e0 * (1 + alpha_sc * (temp_cell - t0))
+
+
+def _pvsyst_Io(effective_irradiance, beta_voc, temp_cell, cells_in_series,
+               e0=1000, t0=25):
+    """Simplified PVsyst saturation current model."""
+    del effective_irradiance  # not used in simplified expression
+    return np.exp(beta_voc / cells_in_series * (temp_cell - t0))
+
+
+def _pvsyst_nNsVth(temp_cell, cells_in_series):
+    """Simplified PVsyst thermal voltage model."""
+    return cells_in_series * constants.k * (temp_cell + 273.15) / constants.e
+
+
+def _pvsyst_gamma(effective_irradiance, gamma_ref, mu_gamma, temp_cell,
+                  e0=1000, t0=25):
+    """Simplified PVsyst gamma model."""
+    del effective_irradiance, e0  # not used
+    return gamma_ref + mu_gamma * (temp_cell - t0)
+
+
 # a dict of required parameter names for each DC power model
 _DC_MODEL_PARAMS = {
     'sapm': {
@@ -419,7 +450,7 @@ def get_iam(self, aoi, iam_model='physical'):
 
         aoi_model : string, default 'physical'
             The IAM model to be used. Valid strings are 'physical', 'ashrae',
-            'martin_ruiz' and 'sapm'.
+            'martin_ruiz', 'fresnel_ar' and 'sapm'.
         Returns
         -------
         iam : numeric or tuple of numeric
@@ -1531,7 +1562,7 @@ def get_iam(self, aoi, iam_model='physical'):
             if `iam_model` is not a valid model name.
         """
         model = iam_model.lower()
-        if model in ['ashrae', 'physical', 'martin_ruiz']:
+        if model in ['ashrae', 'physical', 'martin_ruiz', 'fresnel_ar']:
             param_names = iam._IAM_MODEL_PARAMS[model]
             kwargs = _build_kwargs(param_names, self.module_parameters)
             func = getattr(iam, model)
@@ -2900,7 +2931,7 @@ def singlediode(photocurrent, saturation_current, resistance_series,
 
 
 def max_power_point(photocurrent, saturation_current, resistance_series,
-                    resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.Inf,
+                    resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.inf,
                     method='brentq'):
     """
     Given the single diode equation coefficients, calculates the maximum power
@@ -2947,7 +2978,7 @@ def max_power_point(photocurrent, saturation_current, resistance_series,
     """
     i_mp, v_mp, p_mp = _singlediode.bishop88_mpp(
         photocurrent, saturation_current, resistance_series,
-        resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.Inf,
+        resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.inf,
         method=method.lower()
     )
     if isinstance(photocurrent, pd.Series):
diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index 9f5fd336ef..cd65675adc 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -56,7 +56,7 @@ def estimate_voc(photocurrent, saturation_current, nNsVth):
 
 def bishop88(diode_voltage, photocurrent, saturation_current,
              resistance_series, resistance_shunt, nNsVth, d2mutau=0,
-             NsVbi=np.Inf, breakdown_factor=0., breakdown_voltage=-5.5,
+             NsVbi=np.inf, breakdown_factor=0., breakdown_voltage=-5.5,
              breakdown_exp=3.28, gradients=False):
     r"""
     Explicit calculation of points on the IV curve described by the single
@@ -204,7 +204,7 @@ def bishop88(diode_voltage, photocurrent, saturation_current,
 
 def bishop88_i_from_v(voltage, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth,
-                      d2mutau=0, NsVbi=np.Inf, breakdown_factor=0.,
+                      d2mutau=0, NsVbi=np.inf, breakdown_factor=0.,
                       breakdown_voltage=-5.5, breakdown_exp=3.28,
                       method='newton'):
     """
@@ -292,7 +292,7 @@ def vd_from_brent(voc, v, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
 
 def bishop88_v_from_i(current, photocurrent, saturation_current,
                       resistance_series, resistance_shunt, nNsVth,
-                      d2mutau=0, NsVbi=np.Inf, breakdown_factor=0.,
+                      d2mutau=0, NsVbi=np.inf, breakdown_factor=0.,
                       breakdown_voltage=-5.5, breakdown_exp=3.28,
                       method='newton'):
     """
@@ -378,7 +378,7 @@ def vd_from_brent(voc, i, iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
 
 
 def bishop88_mpp(photocurrent, saturation_current, resistance_series,
-                 resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.Inf,
+                 resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.inf,
                  breakdown_factor=0., breakdown_voltage=-5.5,
                  breakdown_exp=3.28, method='newton'):
     """
diff --git a/pvlib/tests/test_singlediode.py b/pvlib/tests/test_singlediode.py
index ee20e4885e..a5186f4ac9 100644
--- a/pvlib/tests/test_singlediode.py
+++ b/pvlib/tests/test_singlediode.py
@@ -203,7 +203,7 @@ def test_pvsyst_recombination_loss(method, poa, temp_cell, expected, tol):
         # other conditions with breakdown model on and recombination model off
         (
             (1.e-4, -5.5, 3.28),
-            (0., np.Inf),
+            (0., np.inf),
             POA,
             TCELL,
             {
diff --git a/pvlib/tools.py b/pvlib/tools.py
index 229c5dd444..6f8e14995a 100644
--- a/pvlib/tools.py
+++ b/pvlib/tools.py
@@ -458,3 +458,8 @@ def _degrees_to_index(degrees, coordinate):
         index = int(np.around(index))
 
     return index
+
+
+def _first_order_centered_difference(func, x0, args=(), dx=1e-6):
+    """Return the first order centered finite difference of ``func`` at ``x0``."""
+    return (func(x0 + dx, *args) - func(x0 - dx, *args)) / (2 * dx)
diff --git a/tests/test_iam.py b/tests/test_iam.py
index f5ca231bd4..685b32ddf7 100644
--- a/tests/test_iam.py
+++ b/tests/test_iam.py
@@ -90,6 +90,18 @@ def test_physical_scalar():
     assert_allclose(iam, expected, equal_nan=True)
 
 
+def test_fresnel_ar():
+    aoi = np.array([0, 22.5, 45, 67.5, 90, 100, np.nan])
+    expected = np.array([1.0, 0.99972665, 0.99355573, 0.92808897, 0.0, 0.0,
+                         np.nan])
+    iam = _iam.fresnel_ar(aoi)
+    assert_allclose(iam, expected, atol=1e-7, equal_nan=True)
+
+    aoi = pd.Series(aoi)
+    iam_series = _iam.fresnel_ar(aoi)
+    assert_series_equal(iam_series, pd.Series(expected))
+
+
 def test_martin_ruiz():
 
     aoi = 45.
diff --git a/tests/test_pvsystem.py b/tests/test_pvsystem.py
index b58f9fd9e4..18530a22f2 100644
--- a/tests/test_pvsystem.py
+++ b/tests/test_pvsystem.py
@@ -27,6 +27,7 @@
     ('ashrae', {'b': 0.05}),
     ('physical', {'K': 4, 'L': 0.002, 'n': 1.526}),
     ('martin_ruiz', {'a_r': 0.16}),
+    ('fresnel_ar', {'n_ar': 1.2, 'n_air': 1.0, 'n_glass': 1.56}),
 ])
 def test_PVSystem_get_iam(mocker, iam_model, model_params):
     m = mocker.spy(_iam, iam_model)