diff --git a/docs/sphinx/source/reference/pv_modeling.rst b/docs/sphinx/source/reference/pv_modeling.rst
index 4e57130711..797681ce23 100644
--- a/docs/sphinx/source/reference/pv_modeling.rst
+++ b/docs/sphinx/source/reference/pv_modeling.rst
@@ -186,6 +186,7 @@ Utilities for working with IV curve data
    :toctree: generated/
 
     ivtools.utils.rectify_iv_curve
+    ivtools.utils.astm_e1036
 
 Other
 -----
diff --git a/docs/sphinx/source/whatsnew/v0.9.4.rst b/docs/sphinx/source/whatsnew/v0.9.4.rst
index ecb659ccc9..cebb864758 100644
--- a/docs/sphinx/source/whatsnew/v0.9.4.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.4.rst
@@ -45,6 +45,8 @@ Enhancements
   in a simplified way, using the Faiman model as an example.
   :py:func:`~pvlib.temperature.faiman_rad`
   (:issue:`1594`, :pull:`1595`)
+* Add a function :py:func:`pvlib.ivtools.utils.astm_e1036` to perform ASTM E1036 extraction of IV
+  curve parameters (:pull:`1585`)
 
 Bug fixes
 ~~~~~~~~~
@@ -78,6 +80,7 @@ Contributors
 * Christian Orner (:ghuser:`chrisorner`)
 * Saurabh Aneja (:ghuser:`spaneja`)
 * Marcus Boumans (:ghuser:`bowie2211`)
+* Michael Deceglie (:ghuser:`mdeceglie`)
 * Yu Xie (:ghuser:`xieyupku`)
 * Anton Driesse (:ghuser:`adriesse`)
 * Cliff Hansen (:ghuser:`cwhanse`)
diff --git a/pvlib/ivtools/utils.py b/pvlib/ivtools/utils.py
index 17eefa31a0..554af24ecc 100644
--- a/pvlib/ivtools/utils.py
+++ b/pvlib/ivtools/utils.py
@@ -6,6 +6,7 @@
 
 import numpy as np
 import pandas as pd
+from numpy.polynomial.polynomial import Polynomial as Poly
 
 
 # A small number used to decide when a slope is equivalent to zero
@@ -423,3 +424,123 @@ def _schumaker_qspline(x, y):
     yhat = tmp2[:, 4]
     kflag = tmp2[:, 5]
     return t, c, yhat, kflag
+
+
+def astm_e1036(v, i, imax_limits=(0.75, 1.15), vmax_limits=(0.75, 1.15),
+               voc_points=3, isc_points=3, mp_fit_order=4):
+    '''
+    Extract photovoltaic IV parameters according to ASTM E1036. Assumes that
+    the power producing portion of the curve is in the first quadrant.
+
+    Parameters
+    ----------
+    v : array-like
+        Voltage points
+    i : array-like
+        Current points
+    imax_limits : tuple, default (0.75, 1.15)
+        Two-element tuple (low, high) specifying the fraction of estimated
+        Imp within which to fit a polynomial for max power calculation
+    vmax_limits : tuple, default (0.75, 1.15)
+        Two-element tuple (low, high) specifying the fraction of estimated
+        Vmp within which to fit a polynomial for max power calculation
+    voc_points : int, default 3
+        The number of points near open circuit to use for linear fit
+        and Voc calculation
+    isc_points : int, default 3
+        The number of points near short circuit to use for linear fit and
+        Isc calculation
+    mp_fit_order : int, default 4
+        The order of the polynomial fit of power vs. voltage near maximum
+        power
+
+
+    Returns
+    -------
+    dict
+        Results. The IV parameters are given by the keys 'voc', 'isc',
+        'vmp', 'imp', 'pmp', and 'ff'. The key 'mp_fit' gives the numpy
+        Polynomial object for the fit of power vs voltage near maximum
+        power.
+
+    References
+    ----------
+    .. [1] Standard Test Methods for Electrical Performance of Nonconcentrator
+       Terrestrial Photovoltaic Modules and Arrays Using Reference Cells,
+       ASTM E1036-15(2019), :doi:`10.1520/E1036-15R19`
+    '''
+
+    # Adapted from https://github.com/NREL/iv_params
+    # Copyright (c) 2022, Alliance for Sustainable Energy, LLC
+    # All rights reserved.
+
+    df = pd.DataFrame()
+    df['v'] = v
+    df['i'] = i
+    df['p'] = df['v'] * df['i']
+
+    # determine if we can use voc and isc estimates
+    i_min_ind = df['i'].abs().idxmin()
+    v_min_ind = df['v'].abs().idxmin()
+    voc_est = df['v'][i_min_ind]
+    isc_est = df['i'][v_min_ind]
+
+    # accept the estimates if they are close enough
+    # if not, perform a linear fit
+    if abs(df['i'][i_min_ind]) <= isc_est * 0.001:
+        voc = voc_est
+    else:
+        df['i_abs'] = df['i'].abs()
+        voc_df = df.nsmallest(voc_points, 'i_abs')
+        voc_fit = Poly.fit(voc_df['i'], voc_df['v'], 1)
+        voc = voc_fit(0)
+
+    if abs(df['v'][v_min_ind]) <= voc_est * 0.005:
+        isc = isc_est
+    else:
+        df['v_abs'] = df['v'].abs()
+        isc_df = df.nsmallest(isc_points, 'v_abs')
+        isc_fit = Poly.fit(isc_df['v'], isc_df['i'], 1)
+        isc = isc_fit(0)
+
+    # estimate max power point
+    max_index = df['p'].idxmax()
+    mp_est = df.loc[max_index]
+
+    # filter around max power
+    mask = (
+        (df['i'] >= imax_limits[0] * mp_est['i']) &
+        (df['i'] <= imax_limits[1] * mp_est['i']) &
+        (df['v'] >= vmax_limits[0] * mp_est['v']) &
+        (df['v'] <= vmax_limits[1] * mp_est['v'])
+    )
+    filtered = df[mask]
+
+    # fit polynomial and find max
+    mp_fit = Poly.fit(filtered['v'], filtered['p'], mp_fit_order)
+    # Note that this root finding procedure differs from
+    # the suggestion in the standard
+    roots = mp_fit.deriv().roots()
+    # only consider real roots
+    roots = roots.real[abs(roots.imag) < 1e-5]
+    # only consider roots in the relevant part of the domain
+    roots = roots[(roots < filtered['v'].max()) &
+                  (roots > filtered['v'].min())]
+    vmp = roots[np.argmax(mp_fit(roots))]
+    pmp = mp_fit(vmp)
+    # Imp isn't mentioned for update in the
+    # standard, but this seems to be in the intended spirit
+    imp = pmp / vmp
+
+    ff = pmp / (voc * isc)
+
+    result = {}
+    result['voc'] = voc
+    result['isc'] = isc
+    result['vmp'] = vmp
+    result['imp'] = imp
+    result['pmp'] = pmp
+    result['ff'] = ff
+    result['mp_fit'] = mp_fit
+
+    return result
diff --git a/pvlib/tests/ivtools/test_utils.py b/pvlib/tests/ivtools/test_utils.py
index 8f7826bdc2..d8a35e554d 100644
--- a/pvlib/tests/ivtools/test_utils.py
+++ b/pvlib/tests/ivtools/test_utils.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pandas as pd
 import pytest
-from pvlib.ivtools.utils import _numdiff, rectify_iv_curve
+from pvlib.ivtools.utils import _numdiff, rectify_iv_curve, astm_e1036
 from pvlib.ivtools.utils import _schumaker_qspline
 
 from ..conftest import DATA_DIR
@@ -76,3 +76,98 @@ def test__schmumaker_qspline(x, y, expected):
     np.testing.assert_allclose(t, expected[1], atol=0.0001)
     np.testing.assert_allclose(yhat, expected[2], atol=0.0001)
     np.testing.assert_allclose(kflag, expected[3], atol=0.0001)
+
+
+@pytest.fixture
+def i_array():
+    i = np.array([8.09403993, 8.09382549, 8.09361103, 8.09339656, 8.09318205,
+                  8.09296748, 8.09275275, 8.09253771, 8.09232204, 8.09210506,
+                  8.09188538, 8.09166014, 8.09142342, 8.09116305, 8.09085392,
+                  8.09044425, 8.08982734, 8.08878333, 8.08685945, 8.08312463,
+                  8.07566926, 8.06059856, 8.03005836, 7.96856869, 7.8469714,
+                  7.61489584, 7.19789314, 6.51138396, 5.49373476, 4.13267172,
+                  2.46021487, 0.52838624, -1.61055289])
+    return i
+
+
+@pytest.fixture
+def v_array():
+    v = np.array([-0.005, 0.015, 0.035, 0.055, 0.075, 0.095, 0.115, 0.135,
+                  0.155, 0.175, 0.195, 0.215, 0.235, 0.255, 0.275, 0.295,
+                  0.315, 0.335, 0.355, 0.375, 0.395, 0.415, 0.435, 0.455,
+                  0.475, 0.495, 0.515, 0.535, 0.555, 0.575, 0.595, 0.615,
+                  0.635])
+    return v
+
+
+# astm_e1036 tests
+def test_astm_e1036(v_array, i_array):
+    result = astm_e1036(v_array, i_array)
+    expected = {'voc': 0.6195097477985162,
+                'isc': 8.093986320386227,
+                'vmp': 0.494283417170082,
+                'imp': 7.626088301548568,
+                'pmp': 3.7694489853302127,
+                'ff': 0.7517393078504361}
+    fit = result.pop('mp_fit')
+    expected_fit = np.array(
+        [3.6260726, 0.49124176, -0.24644747, -0.26442383, -0.1223237])
+    assert fit.coef == pytest.approx(expected_fit)
+    assert result == pytest.approx(expected)
+
+
+def test_astm_e1036_fit_order(v_array, i_array):
+    result = astm_e1036(v_array, i_array, mp_fit_order=3)
+    fit = result.pop('mp_fit')
+    expected_fit = np.array(
+        [3.64081697, 0.49124176, -0.3720477, -0.26442383])
+    assert fit.coef == pytest.approx(expected_fit)
+
+
+def test_astm_e1036_est_isc_voc(v_array, i_array):
+    '''
+    Test the case in which Isc and Voc estimates are
+    valid without a linear fit
+    '''
+    v = v_array
+    i = i_array
+    v = np.append(v, [0.001, 0.6201])
+    i = np.append(i, [8.09397560e+00, 7.10653445e-04])
+    result = astm_e1036(v, i)
+    expected = {'voc': 0.6201,
+                'isc': 8.093975598317805,
+                'vmp': 0.494283417170082,
+                'imp': 7.626088301548568,
+                'pmp': 3.7694489853302127,
+                'ff': 0.751024747526615}
+    result.pop('mp_fit')
+    assert result == pytest.approx(expected)
+
+
+def test_astm_e1036_mpfit_limits(v_array, i_array):
+    result = astm_e1036(v_array,
+                        i_array,
+                        imax_limits=(0.85, 1.1),
+                        vmax_limits=(0.85, 1.1))
+    expected = {'voc': 0.6195097477985162,
+                'isc': 8.093986320386227,
+                'vmp': 0.49464214190725303,
+                'imp': 7.620032530519718,
+                'pmp': 3.769189212299219,
+                'ff': 0.7516875014460312}
+    result.pop('mp_fit')
+    assert result == pytest.approx(expected)
+
+
+def test_astm_e1036_fit_points(v_array, i_array):
+    i = i_array
+    i[3] = 8.1  # ensure an interesting change happens
+    result = astm_e1036(v_array, i, voc_points=4, isc_points=4)
+    expected = {'voc': 0.619337073271274,
+                'isc': 8.093160893325297,
+                'vmp': 0.494283417170082,
+                'imp': 7.626088301548568,
+                'pmp': 3.7694489853302127,
+                'ff': 0.7520255886236707}
+    result.pop('mp_fit')
+    assert result == pytest.approx(expected)