diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml
index 8e054ecb08..09c022f51c 100644
--- a/.github/workflows/flake8.yml
+++ b/.github/workflows/flake8.yml
@@ -12,7 +12,7 @@ jobs:
         with:
           python-version: '3.11'
       - name: Install Flake8 5.0.4 linter
-        run: pip install flake8==5.0.4  # use this version for --diff option
+        run: pip install flake8==5.0.4  # use this version for --diff option, should match version in pyproject.toml
       - name: Setup Flake8 output matcher for PR annotations
         run: echo '::add-matcher::.github/workflows/flake8-linter-matcher.json'
       - name: Fetch pull request target branch
diff --git a/docs/examples/spectrum/spectral_factor.py b/docs/examples/spectrum/spectral_factor.py
index 83ee488fb4..a1638340c4 100644
--- a/docs/examples/spectrum/spectral_factor.py
+++ b/docs/examples/spectrum/spectral_factor.py
@@ -46,7 +46,8 @@
 # spectral factor functions:
 #
 # - :py:func:`~pvlib.spectrum.spectral_factor_sapm`, which requires only
-#   the absolute airmass, :math:`AM_a`
+#   the absolute airmass, :math:`AM_a` and five SAPM coefficients :math:`A0` -
+#   :math:`A4`
 # - :py:func:`~pvlib.spectrum.spectral_factor_pvspec`, which requires
 #   :math:`AM_a` and the clearsky index, :math:`k_c`
 # - :py:func:`~pvlib.spectrum.spectral_factor_firstsolar`, which requires
@@ -104,7 +105,14 @@
 module = pvlib.pvsystem.retrieve_sam('SandiaMod')['LG_LG290N1C_G3__2013_']
 #
 # Calculate M using the three models for an mc-Si PV module.
-m_sapm = pvlib.spectrum.spectral_factor_sapm(airmass_absolute, module)
+m_sapm = pvlib.spectrum.spectral_factor_sapm(
+    airmass_absolute,
+    module["A0"],
+    module["A1"],
+    module["A2"],
+    module["A3"],
+    module["A4"],
+)
 m_pvspec = pvlib.spectrum.spectral_factor_pvspec(airmass_absolute, kc,
                                                  'multisi')
 m_fs = pvlib.spectrum.spectral_factor_firstsolar(w, airmass_absolute,
diff --git a/pvlib/iam.py b/pvlib/iam.py
index 161de84589..eb9326245b 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -7,22 +7,14 @@
 IAM is typically a function of the angle of incidence (AOI) of the direct
 irradiance to the module's surface.
 """
+import functools
 
 import numpy as np
 import pandas as pd
-import functools
+import scipy.interpolate
 from scipy.optimize import minimize
-from pvlib.tools import cosd, sind, acosd
 
-# a dict of required parameter names for each IAM model
-# keys are the function names for the IAM models
-_IAM_MODEL_PARAMS = {
-    'ashrae': {'b'},
-    'physical': {'n', 'K', 'L'},
-    'martin_ruiz': {'a_r'},
-    'sapm': {'B0', 'B1', 'B2', 'B3', 'B4', 'B5'},
-    'interp': {'theta_ref', 'iam_ref'}
-}
+from pvlib.tools import cosd, sind, acosd
 
 
 def ashrae(aoi, b=0.05):
@@ -75,9 +67,11 @@ def ashrae(aoi, b=0.05):
 
     See Also
     --------
-    pvlib.iam.physical
-    pvlib.iam.martin_ruiz
     pvlib.iam.interp
+    pvlib.iam.martin_ruiz
+    pvlib.iam.physical
+    pvlib.iam.sapm
+    pvlib.iam.schlick
     """
 
     iam = 1 - b * (1 / np.cos(np.radians(aoi)) - 1)
@@ -92,7 +86,7 @@ def ashrae(aoi, b=0.05):
     return iam
 
 
-def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None):
+def physical(aoi, n=1.526, K=4.0, L=0.002, n_ar=None):
     r"""
     Determine the incidence angle modifier using refractive index ``n``,
     extinction coefficient ``K``, glazing thickness ``L`` and refractive
@@ -151,10 +145,11 @@ def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None):
 
     See Also
     --------
-    pvlib.iam.martin_ruiz
     pvlib.iam.ashrae
     pvlib.iam.interp
+    pvlib.iam.martin_ruiz
     pvlib.iam.sapm
+    pvlib.iam.schlick
     """
     n1, n3 = 1, n
     if n_ar is None or np.allclose(n_ar, n1):
@@ -286,11 +281,12 @@ def martin_ruiz(aoi, a_r=0.16):
 
     See Also
     --------
-    pvlib.iam.martin_ruiz_diffuse
-    pvlib.iam.physical
     pvlib.iam.ashrae
     pvlib.iam.interp
+    pvlib.iam.physical
     pvlib.iam.sapm
+    pvlib.iam.schlick
+    pvlib.iam.martin_ruiz_diffuse
     '''
     # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019
 
@@ -314,7 +310,7 @@ def martin_ruiz(aoi, a_r=0.16):
 
 def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
     '''
-    Determine the incidence angle modifiers (iam) for diffuse sky and
+    Determine the incidence angle modifiers (IAM) for diffuse sky and
     ground-reflected irradiance using the Martin and Ruiz incident angle model.
 
     Parameters
@@ -373,10 +369,8 @@ def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
     See Also
     --------
     pvlib.iam.martin_ruiz
-    pvlib.iam.physical
-    pvlib.iam.ashrae
-    pvlib.iam.interp
-    pvlib.iam.sapm
+    pvlib.iam.marion_diffuse
+    pvlib.iam.schlick_diffuse
     '''
     # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Oct. 2019
 
@@ -463,15 +457,14 @@ def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
 
     See Also
     --------
-    pvlib.iam.physical
     pvlib.iam.ashrae
     pvlib.iam.martin_ruiz
+    pvlib.iam.physical
     pvlib.iam.sapm
+    pvlib.iam.schlick
     '''
     # 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}
 
@@ -483,8 +476,9 @@ def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
         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')
+    interpolator = scipy.interpolate.interp1d(
+        theta_ref, iam_ref, kind=method, fill_value='extrapolate'
+    )
     aoi_input = aoi
 
     aoi = np.asanyarray(aoi)
@@ -501,9 +495,19 @@ def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
     return iam
 
 
-def sapm(aoi, module, upper=None):
+def sapm(aoi, B0, B1, B2, B3, B4, B5, upper=None):
     r"""
-    Determine the incidence angle modifier (IAM) using the SAPM model.
+    Caclulate the incidence angle modifier (IAM), :math:`f_2`, using the
+    Sandia Array Performance Model (SAPM).
+
+    The SAPM incidence angle modifier is part of the broader Sandia Array
+    Performance Model, which defines five points on an IV curve using empirical
+    module-specific coefficients. Module coefficients for the SAPM are
+    available in the SAPM database and can be retrieved for use through
+    :py:func:`pvlib.pvsystem.retrieve_sam()`. More details on the SAPM can be
+    found in [1]_, while a full description of the procedure to determine the
+    empirical model coefficients, including those for the SAPM incidence angle
+    modifier, can be found in [2]_.
 
     Parameters
     ----------
@@ -511,20 +515,46 @@ def sapm(aoi, module, upper=None):
         Angle of incidence in degrees. Negative input angles will return
         zeros.
 
-    module : dict-like
-        A dict or Series with the SAPM IAM model parameters.
-        See the :py:func:`sapm` notes section for more details.
+    B0 : float
+        The coefficient of the degree-0 polynomial term.
+
+    B1 : float
+        The coefficient of the degree-1 polynomial term.
+
+    B2 : float
+        The coefficient of the degree-2 polynomial term.
+
+    B3 : float
+        The coefficient of the degree-3 polynomial term.
+
+    B4 : float
+        The coefficient of the degree-4 polynomial term.
+
+    B5 : float
+        The coefficient of the degree-5 polynomial term.
 
     upper : float, optional
-        Upper limit on the results.
+        Upper limit on the results. None means no upper limiting.
 
     Returns
     -------
     iam : numeric
-        The SAPM angle of incidence loss coefficient, termed F2 in [1]_.
+        The SAPM angle of incidence loss coefficient, :math:`f_2` in [1]_.
 
     Notes
     -----
+    The SAPM spectral correction functions parameterises :math:`f_2` as a
+    fifth-order polynomial function of angle of incidence:
+
+    .. math::
+
+        f_2 = b_0 + b_1 AOI + b_2 AOI^2 + b_3 AOI^3 + b_4 AOI^4 + b_5 AOI^5.
+
+    where :math:`f_2` is the spectral mismatch factor, :math:`b_{0-5}` are
+    the module-specific coefficients, and :math:`AOI` is the angle of
+    incidence. More detail on how this incidence angle modifier function was
+    developed can be found in [3]_. Its measurement is described in [4]_.
+
     The SAPM [1]_ traditionally does not define an upper limit on the AOI
     loss function and values slightly exceeding 1 may exist for moderate
     angles of incidence (15-40 degrees). However, users may consider
@@ -532,30 +562,34 @@ def sapm(aoi, module, upper=None):
 
     References
     ----------
-    .. [1] King, D. et al, 2004, "Sandia Photovoltaic Array Performance
-       Model", SAND Report 3535, Sandia National Laboratories, Albuquerque,
-       NM.
-
-    .. [2] B.H. King et al, "Procedure to Determine Coefficients for the
-       Sandia Array Performance Model (SAPM)," SAND2016-5284, Sandia
-       National Laboratories (2016).
-
-    .. [3] B.H. King et al, "Recent Advancements in Outdoor Measurement
-       Techniques for Angle of Incidence Effects," 42nd IEEE PVSC (2015).
-       :doi:`10.1109/PVSC.2015.7355849`
+    .. [1] King, D., Kratochvil, J., and Boyson W. (2004), "Sandia
+           Photovoltaic Array Performance Model", (No. SAND2004-3535), Sandia
+           National Laboratories, Albuquerque, NM (United States).
+           :doi:`10.2172/919131`
+    .. [2] King, B., Hansen, C., Riley, D., Robinson, C., and Pratt, L.
+           (2016). Procedure to determine coefficients for the Sandia Array
+           Performance Model (SAPM) (No. SAND2016-5284). Sandia National
+           Laboratories, Albuquerque, NM (United States).
+           :doi:`10.2172/1256510`
+    .. [3] King, D., Kratochvil, J., and Boyson, W. "Measuring solar spectral
+           and angle-of-incidence effects on photovoltaic modules and solar
+           irradiance sensors." Conference Record of the 26th IEEE Potovoltaic
+           Specialists Conference (PVSC). IEEE, 1997.
+           :doi:`10.1109/PVSC.1997.654283`
+    .. [4] B.H. King et al, "Recent Advancements in Outdoor Measurement
+           Techniques for Angle of Incidence Effects," 42nd IEEE PVSC (2015).
+           :doi:`10.1109/PVSC.2015.7355849`
 
     See Also
     --------
-    pvlib.iam.physical
     pvlib.iam.ashrae
-    pvlib.iam.martin_ruiz
     pvlib.iam.interp
+    pvlib.iam.martin_ruiz
+    pvlib.iam.physical
+    pvlib.iam.schlick
     """
 
-    aoi_coeff = [module['B5'], module['B4'], module['B3'], module['B2'],
-                 module['B1'], module['B0']]
-
-    iam = np.polyval(aoi_coeff, aoi)
+    iam = np.polyval([B5, B4, B3, B2, B1, B0], aoi)
     iam = np.clip(iam, 0, upper)
     # nan tolerant masking
     aoi_lt_0 = np.full_like(aoi, False, dtype='bool')
@@ -575,9 +609,11 @@ def marion_diffuse(model, surface_tilt, **kwargs):
 
     Parameters
     ----------
-    model : str
+    model : str or callable
         The IAM function to evaluate across solid angle. Must be one of
-        `'ashrae', 'physical', 'martin_ruiz', 'sapm', 'schlick'`.
+        `'ashrae', `interp`, 'martin_ruiz', 'physical', 'sapm', 'schlick'`
+        or be callable. If callable, then must take numeric AOI in degrees
+        as input and return the fractional IAM.
 
     surface_tilt : numeric
         Surface tilt angles in decimal degrees.
@@ -602,6 +638,14 @@ def marion_diffuse(model, surface_tilt, **kwargs):
     See Also
     --------
     pvlib.iam.marion_integrate
+    pvlib.iam.ashrae
+    pvlib.iam.interp
+    pvlib.iam.martin_ruiz
+    pvlib.iam.physical
+    pvlib.iam.sapm
+    pvlib.iam.schlick
+    pvlib.iam.martin_ruiz_diffuse
+    pvlib.iam.schlick_diffuse
 
     References
     ----------
@@ -613,7 +657,7 @@ def marion_diffuse(model, surface_tilt, **kwargs):
     Examples
     --------
     >>> marion_diffuse('physical', surface_tilt=20)
-    {'sky': 0.9539178294437575,
+    {'sky': 0.953917829443757
      'horizon': 0.7652650139134007,
      'ground': 0.6387140117795903}
 
@@ -622,26 +666,28 @@ def marion_diffuse(model, surface_tilt, **kwargs):
      'horizon': array([0.86478428, 0.91825792]),
      'ground': array([0.77004435, 0.8522436 ])}
     """
+    if callable(model):
+        # A callable IAM model function was specified.
+        func = model
+    else:
+        # Check that a builtin IAM function was specified.
+        builtin_models = _get_builtin_models()
 
-    models = {
-        'physical': physical,
-        'ashrae': ashrae,
-        'sapm': sapm,
-        'martin_ruiz': martin_ruiz,
-        'schlick': schlick,
-    }
+        try:
+            model = builtin_models[model]
+        except KeyError as exc:
+            raise ValueError(
+                f'model must be one of: {builtin_models.keys()}'
+            ) from exc
 
-    try:
-        iam_model = models[model]
-    except KeyError:
-        raise ValueError('model must be one of: ' + str(list(models.keys())))
+        func = model["func"]
 
-    iam_function = functools.partial(iam_model, **kwargs)
-    iam = {}
-    for region in ['sky', 'horizon', 'ground']:
-        iam[region] = marion_integrate(iam_function, surface_tilt, region)
+    iam_function = functools.partial(func, **kwargs)
 
-    return iam
+    return {
+        region: marion_integrate(iam_function, surface_tilt, region)
+        for region in ['sky', 'horizon', 'ground']
+    }
 
 
 def marion_integrate(function, surface_tilt, region, num=None):
@@ -814,10 +860,6 @@ def schlick(aoi):
     iam : numeric
         The incident angle modifier.
 
-    See Also
-    --------
-    pvlib.iam.schlick_diffuse
-
     References
     ----------
     .. [1] Schlick, C. An inexpensive BRDF model for physically-based
@@ -827,6 +869,15 @@ def schlick(aoi):
        for Diffuse radiation on Inclined photovoltaic Surfaces (FEDIS)",
        Renewable and Sustainable Energy Reviews, vol. 161, 112362. June 2022.
        :doi:`10.1016/j.rser.2022.112362`
+
+    See Also
+    --------
+    pvlib.iam.ashrae
+    pvlib.iam.interp
+    pvlib.iam.martin_ruiz
+    pvlib.iam.physical
+    pvlib.iam.sapm
+    pvlib.iam.schlick_diffuse
     """
     iam = 1 - (1 - cosd(aoi)) ** 5
     iam = np.where(np.abs(aoi) >= 90.0, 0.0, iam)
@@ -875,10 +926,6 @@ def schlick_diffuse(surface_tilt):
     iam_ground : numeric
         The incident angle modifier for ground-reflected diffuse.
 
-    See Also
-    --------
-    pvlib.iam.schlick
-
     Notes
     -----
     The analytical integration of the Schlick approximation was derived
@@ -913,6 +960,12 @@ def schlick_diffuse(surface_tilt):
        for Diffuse radiation on Inclined photovoltaic Surfaces (FEDIS)",
        Renewable and Sustainable Energy Reviews, vol. 161, 112362. June 2022.
        :doi:`10.1016/j.rser.2022.112362`
+
+    See Also
+    --------
+    pvlib.iam.schlick
+    pvlib.iam.marion_diffuse
+    pvlib.iam.martin_ruiz_diffuse
     """
     # these calculations are as in [2]_, but with the refractive index
     # weighting coefficient w set to 1.0 (so it is omitted)
@@ -943,27 +996,96 @@ def schlick_diffuse(surface_tilt):
     return cuk, cug
 
 
-def _get_model(model_name):
-    # check that model is implemented
-    model_dict = {'ashrae': ashrae, 'martin_ruiz': martin_ruiz,
-                  'physical': physical}
+def _get_builtin_models():
+    """
+    Get builtin IAM models' usage information.
+
+    Returns
+    -------
+    info : dict
+        A dictionary of dictionaries keyed by builtin IAM model name, with
+        each model dictionary containing:
+
+        * 'func': callable
+            The callable model function
+        * 'params_required': set of str
+            The model function's required parameters
+        * 'params_optional': set of str
+            The model function's optional parameters
+
+    See Also
+    --------
+    pvlib.iam.ashrae
+    pvlib.iam.interp
+    pvlib.iam.martin_ruiz
+    pvlib.iam.physical
+    pvlib.iam.sapm
+    pvlib.iam.schlick
+    """
+    return {
+        'ashrae': {
+            'func': ashrae,
+            'params_required': set(),
+            'params_optional': {'b'},
+        },
+        'interp': {
+            'func': interp,
+            'params_required': {'theta_ref', 'iam_ref'},
+            'params_optional': {'method', 'normalize'},
+        },
+        'martin_ruiz': {
+            'func': martin_ruiz,
+            'params_required': set(),
+            'params_optional': {'a_r'},
+        },
+        'physical': {
+            'func': physical,
+            'params_required': set(),
+            'params_optional': {'n', 'K', 'L', 'n_ar'},
+        },
+        'sapm': {
+            'func': sapm,
+            'params_required': {'B0', 'B1', 'B2', 'B3', 'B4', 'B5'},
+            'params_optional': {'upper'},
+        },
+        'schlick': {
+            'func': schlick,
+            'params_required': set(),
+            'params_optional': set(),
+        },
+    }
+
+
+def _get_fittable_or_convertable_model(builtin_model_name):
+    # check that model is implemented and fittable or convertable
+    implemented_builtin_models = {
+        'ashrae': ashrae, 'martin_ruiz': martin_ruiz, 'physical': physical,
+    }
+
     try:
-        model = model_dict[model_name]
-    except KeyError:
-        raise NotImplementedError(f"The {model_name} model has not been "
-                                  "implemented")
+        return implemented_builtin_models[builtin_model_name]
+    except KeyError as exc:
+        raise NotImplementedError(
+            f"No implementation for model {builtin_model_name}"
+        ) from exc
+
 
-    return model
+def _check_params(builtin_model_name, params):
+    # check that parameters passed in with IAM model belong to the model
+    handed_params = set(params.keys())
+    builtin_model = _get_builtin_models()[builtin_model_name]
+    expected_params = builtin_model["params_required"].union(
+        builtin_model["params_optional"]
+    )
 
+    if builtin_model_name == 'physical':
+        expected_params.remove('n_ar')  # Not required.
 
-def _check_params(model_name, params):
-    # check that the parameters passed in with the model
-    # belong to the model
-    exp_params = _IAM_MODEL_PARAMS[model_name]
-    if set(params.keys()) != exp_params:
-        raise ValueError(f"The {model_name} model was expecting to be passed "
-                         "{', '.join(list(exp_params))}, but "
-                         "was handed {', '.join(list(params.keys()))}")
+    if handed_params != expected_params:
+        raise ValueError(
+            f"The {builtin_model_name} model was expecting to be passed"
+            f"{expected_params}, but was handed {handed_params}"
+        )
 
 
 def _sin_weight(aoi):
@@ -1179,8 +1301,8 @@ def convert(source_name, source_params, target_name, weight=_sin_weight,
     pvlib.iam.martin_ruiz
     pvlib.iam.physical
     """
-    source = _get_model(source_name)
-    target = _get_model(target_name)
+    source = _get_fittable_or_convertable_model(source_name)
+    target = _get_fittable_or_convertable_model(target_name)
 
     aoi = np.linspace(0, 90, 91)
     _check_params(source_name, source_params)
@@ -1277,7 +1399,7 @@ def fit(measured_aoi, measured_iam, model_name, weight=_sin_weight, xtol=None):
     pvlib.iam.martin_ruiz
     pvlib.iam.physical
     """
-    target = _get_model(model_name)
+    target = _get_fittable_or_convertable_model(model_name)
 
     if model_name == "physical":
         bounds = [(0, 0.08), (1, 2)]
diff --git a/pvlib/modelchain.py b/pvlib/modelchain.py
index 8456aac114..365afd7b4f 100644
--- a/pvlib/modelchain.py
+++ b/pvlib/modelchain.py
@@ -6,14 +6,15 @@
 the time to read the source code for the module.
 """
 
+from dataclasses import dataclass, field
 from functools import partial
 import itertools
+from typing import Optional, Tuple, TypeVar, Union
 import warnings
+
 import pandas as pd
-from dataclasses import dataclass, field
-from typing import Union, Tuple, Optional, TypeVar
 
-from pvlib import pvsystem, iam
+from pvlib import pvsystem
 import pvlib.irradiance  # avoid name conflict with full import
 from pvlib.pvsystem import _DC_MODEL_PARAMS
 from pvlib.tools import _build_kwargs
@@ -289,12 +290,12 @@ class ModelChain:
     Parameters
     ----------
     system : PVSystem
-        A :py:class:`~pvlib.pvsystem.PVSystem` object that represents
-        the connected set of modules, inverters, etc.
+        A :py:class:`~pvlib.pvsystem.PVSystem` object that represents the
+        connected set of modules, inverters, etc.
 
     location : Location
-        A :py:class:`~pvlib.location.Location` object that represents
-        the physical location at which to evaluate the model.
+        A :py:class:`~pvlib.location.Location` object that represents the
+        physical location at which to evaluate the model.
 
     clearsky_model : str, default 'ineichen'
         Passed to location.get_clearsky. Only used when DNI is not found in
@@ -311,31 +312,29 @@ class ModelChain:
 
     dc_model : str, or function, optional
         If not specified, the model will be inferred from the parameters that
-        are common to all of system.arrays[i].module_parameters.
-        Valid strings are 'sapm', 'desoto', 'cec', 'pvsyst', 'pvwatts'.
-        The ModelChain instance will be passed as the first argument
-        to a user-defined function.
+        are common to all of system.arrays[i].module_parameters. Valid strings
+        are 'sapm', 'desoto', 'cec', 'pvsyst', 'pvwatts'. The ModelChain
+        instance will be passed as the first argument to a user-defined
+        function.
 
     ac_model : str, or function, optional
         If not specified, the model will be inferred from the parameters that
-        are common to all of system.inverter_parameters.
-        Valid strings are 'sandia', 'adr', 'pvwatts'. The
-        ModelChain instance will be passed as the first argument to a
-        user-defined function.
+        are common to all of system.inverter_parameters. Valid strings are
+        'sandia', 'adr', 'pvwatts'. The ModelChain instance will be passed as
+        the first argument to a user-defined function.
 
     aoi_model : str, or function, optional
         If not specified, the model will be inferred from the parameters that
-        are common to all of system.arrays[i].module_parameters.
-        Valid strings are 'physical', 'ashrae', 'sapm', 'martin_ruiz',
-        'interp' and 'no_loss'. The ModelChain instance will be passed as the
-        first argument to a user-defined function.
+        are common to all of system.arrays[i].module_parameters. Valid strings
+        are 'ashrae', 'interp', 'martin_ruiz', 'physical', 'sapm', 'schlick',
+        'no_loss'. The ModelChain instance will be passed as the first
+        argument to a user-defined function.
 
     spectral_model : str, or function, optional
         If not specified, the model will be inferred from the parameters that
-        are common to all of system.arrays[i].module_parameters.
-        Valid strings are 'sapm', 'first_solar', 'no_loss'.
-        The ModelChain instance will be passed as the first argument to
-        a user-defined function.
+        are common to all of system.arrays[i].module_parameters. Valid strings
+        are 'sapm', 'first_solar', 'no_loss'. The ModelChain instance will be
+        passed as the first argument to a user-defined function.
 
     temperature_model : str or function, optional
         Valid strings are: 'sapm', 'pvsyst', 'faiman', 'fuentes', 'noct_sam'.
@@ -386,7 +385,6 @@ def __init__(self, system, location,
 
         self.results = ModelChainResult()
 
-
     @classmethod
     def with_pvwatts(cls, system, location,
                      clearsky_model='ineichen',
@@ -770,62 +768,98 @@ def aoi_model(self, model):
             model = model.lower()
             if model == 'ashrae':
                 self._aoi_model = self.ashrae_aoi_loss
-            elif model == 'physical':
-                self._aoi_model = self.physical_aoi_loss
-            elif model == 'sapm':
-                self._aoi_model = self.sapm_aoi_loss
-            elif model == 'martin_ruiz':
-                self._aoi_model = self.martin_ruiz_aoi_loss
             elif model == 'interp':
                 self._aoi_model = self.interp_aoi_loss
+            elif model == 'martin_ruiz':
+                self._aoi_model = self.martin_ruiz_aoi_loss
             elif model == 'no_loss':
                 self._aoi_model = self.no_aoi_loss
+            elif model == 'physical':
+                self._aoi_model = self.physical_aoi_loss
+            elif model == 'sapm':
+                self._aoi_model = self.sapm_aoi_loss
+            elif model == 'schlick':
+                self._aoi_model = self.schlick_aoi_loss
             else:
                 raise ValueError(model + ' is not a valid aoi loss model')
-        else:
+        elif callable(model):
             self._aoi_model = partial(model, self)
+        else:
+            raise ValueError(model + ' is not a valid aoi loss model')
 
     def infer_aoi_model(self):
+        """
+        Infer AOI model by checking for all required model paramters or at
+        least one optional model parameter in module_parameter collected
+        across all arrays.
+        """
         module_parameters = tuple(
-            array.module_parameters for array in self.system.arrays)
+            array.module_parameters for array in self.system.arrays
+        )
         params = _common_keys(module_parameters)
-        if iam._IAM_MODEL_PARAMS['physical'] <= params:
-            return self.physical_aoi_loss
-        elif iam._IAM_MODEL_PARAMS['sapm'] <= params:
-            return self.sapm_aoi_loss
-        elif iam._IAM_MODEL_PARAMS['ashrae'] <= params:
+        builtin_models = pvlib.iam._get_builtin_models()
+
+        if (builtin_models['ashrae']["params_required"] and (
+            builtin_models['ashrae']["params_required"] <= params)) or (
+            not builtin_models['ashrae']["params_required"] and
+            (builtin_models['ashrae']["params_optional"] & params)
+        ):
             return self.ashrae_aoi_loss
-        elif iam._IAM_MODEL_PARAMS['martin_ruiz'] <= params:
-            return self.martin_ruiz_aoi_loss
-        elif iam._IAM_MODEL_PARAMS['interp'] <= params:
+
+        if (builtin_models['interp']["params_required"] and (
+            builtin_models['interp']["params_required"] <= params)) or (
+            not builtin_models['interp']["params_required"] and
+            (builtin_models['interp']["params_optional"] & params)
+        ):
             return self.interp_aoi_loss
-        else:
-            raise ValueError('could not infer AOI model from '
-                             'system.arrays[i].module_parameters. Check that '
-                             'the module_parameters for all Arrays in '
-                             'system.arrays contain parameters for the '
-                             'physical, aoi, ashrae, martin_ruiz or interp '
-                             'model; explicitly set the model with the '
-                             'aoi_model kwarg; or set aoi_model="no_loss".')
+
+        if (builtin_models['martin_ruiz']["params_required"] and (
+            builtin_models['martin_ruiz']["params_required"] <= params)) or (
+            not builtin_models['martin_ruiz']["params_required"] and
+            (builtin_models['martin_ruiz']["params_optional"] & params)
+        ):
+            return self.martin_ruiz_aoi_loss
+
+        if (builtin_models['physical']["params_required"] and (
+            builtin_models['physical']["params_required"] <= params)) or (
+            not builtin_models['physical']["params_required"] and
+            (builtin_models['physical']["params_optional"] & params)
+        ):
+            return self.physical_aoi_loss
+
+        if (builtin_models['sapm']["params_required"] and (
+            builtin_models['sapm']["params_required"] <= params)) or (
+            not builtin_models['sapm']["params_required"] and
+            (builtin_models['sapm']["params_optional"] & params)
+        ):
+            return self.sapm_aoi_loss
+
+        # schlick model has no parameters to distinguish, esp. from no_loss.
+
+        raise ValueError(
+            'could not infer AOI model from '
+            'system.arrays[i].module_parameters. Check that the '
+            'module_parameters for all Arrays in system.arrays contain '
+            'parameters for the ashrae, interp, martin_ruiz, physical, or '
+            'sapm model; explicitly set the model (esp. the parameter-free '
+            'schlick) with the aoi_model kwarg or set aoi_model="no_loss".'
+        )
 
     def ashrae_aoi_loss(self):
         self.results.aoi_modifier = self.system.get_iam(
-            self.results.aoi,
-            iam_model='ashrae'
+            self.results.aoi, iam_model='ashrae'
         )
         return self
 
     def physical_aoi_loss(self):
         self.results.aoi_modifier = self.system.get_iam(
-            self.results.aoi,
-            iam_model='physical'
+            self.results.aoi, iam_model='physical'
         )
         return self
 
     def sapm_aoi_loss(self):
         self.results.aoi_modifier = self.system.get_iam(
-            self.results.aoi,
-            iam_model='sapm'
+            self.results.aoi, iam_model='sapm'
         )
         return self
 
@@ -837,8 +871,13 @@ def martin_ruiz_aoi_loss(self):
 
     def interp_aoi_loss(self):
         self.results.aoi_modifier = self.system.get_iam(
-            self.results.aoi,
-            iam_model='interp'
+            self.results.aoi, iam_model='interp'
+        )
+        return self
+
+    def schlick_aoi_loss(self):
+        self.results.aoi_modifier = self.system.get_iam(
+            self.results.aoi, iam_model='schlick'
         )
         return self
 
@@ -867,8 +906,10 @@ def spectral_model(self, model):
                 self._spectral_model = self.no_spectral_loss
             else:
                 raise ValueError(model + ' is not a valid spectral loss model')
-        else:
+        elif callable(model):
             self._spectral_model = partial(model, self)
+        else:
+            raise ValueError(model + ' is not a valid spectral loss model')
 
     def infer_spectral_model(self):
         """Infer spectral model from system attributes."""
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index f99f0275af..f4c76723fb 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -8,7 +8,6 @@
 import io
 import itertools
 from pathlib import Path
-import inspect
 from urllib.request import urlopen
 import numpy as np
 from scipy import constants
@@ -399,8 +398,8 @@ def get_iam(self, aoi, iam_model='physical'):
             The angle of incidence in degrees.
 
         aoi_model : string, default 'physical'
-            The IAM model to be used. Valid strings are 'physical', 'ashrae',
-            'martin_ruiz', 'sapm' and 'interp'.
+            The IAM model to be used. Valid strings are 'ashrae', 'interp',
+            'martin_ruiz', 'physical', 'sapm', and  'schlick'.
         Returns
         -------
         iam : numeric or tuple of numeric
@@ -634,9 +633,14 @@ def sapm_spectral_loss(self, airmass_absolute):
             The SAPM spectral loss coefficient.
         """
         return tuple(
-            spectrum.spectral_factor_sapm(airmass_absolute,
-                                          array.module_parameters)
-            for array in self.arrays
+            spectrum.spectral_factor_sapm(
+                airmass_absolute,
+                array.module_parameters["A0"],
+                array.module_parameters["A1"],
+                array.module_parameters["A2"],
+                array.module_parameters["A3"],
+                array.module_parameters["A4"],
+            ) for array in self.arrays
         )
 
     @_unwrap_single_value
@@ -1161,12 +1165,15 @@ def get_irradiance(self, solar_zenith, solar_azimuth, dni, ghi, dhi,
 
     def get_iam(self, aoi, iam_model='physical'):
         """
-        Determine the incidence angle modifier using the method specified by
+        Determine the incidence-angle modifier (IAM) for the given
+        angle of incidence (AOI) using the method specified by
         ``iam_model``.
 
         Parameters for the selected IAM model are expected to be in
-        ``Array.module_parameters``. Default parameters are available for
-        the 'physical', 'ashrae' and 'martin_ruiz' models.
+        ``Array.module_parameters``. A minimal set of default parameters
+        are available for the 'ashrae', 'martin_ruiz', and 'physical'
+        models, but not for 'interp', 'sapm'. The 'schlick' model does not
+        take any parameters.
 
         Parameters
         ----------
@@ -1174,13 +1181,13 @@ def get_iam(self, aoi, iam_model='physical'):
             The angle of incidence in degrees.
 
         aoi_model : string, default 'physical'
-            The IAM model to be used. Valid strings are 'physical', 'ashrae',
-            'martin_ruiz', 'sapm' and 'interp'.
+            The IAM model to be used. Valid strings are 'ashrae', 'interp',
+            'martin_ruiz', 'physical', 'sapm', and 'schlick'.
 
         Returns
         -------
         iam : numeric
-            The AOI modifier.
+            The IAM at the specified AOI.
 
         Raises
         ------
@@ -1188,18 +1195,27 @@ 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', 'interp']:
-            func = getattr(iam, model)  # get function at pvlib.iam
-            # get all parameters from function signature to retrieve them from
-            # module_parameters if present
-            params = set(inspect.signature(func).parameters.keys())
-            params.discard('aoi')  # exclude aoi so it can't be repeated
-            kwargs = _build_kwargs(params, self.module_parameters)
-            return func(aoi, **kwargs)
-        elif model == 'sapm':
-            return iam.sapm(aoi, self.module_parameters)
-        else:
-            raise ValueError(model + ' is not a valid IAM model')
+
+        try:
+            model_info = iam._get_builtin_models()[model]
+        except KeyError as exc:
+            raise ValueError(f'{iam_model} is not a valid iam_model') from exc
+
+        for param in model_info["params_required"]:
+            if param not in self.module_parameters:
+                raise KeyError(f"{param} is missing in module_parameters")
+
+        params_required = _build_kwargs(
+            model_info["params_required"], self.module_parameters
+        )
+
+        params_optional = _build_kwargs(
+            model_info["params_optional"], self.module_parameters
+        )
+
+        return model_info["func"](
+            aoi, **params_required, **params_optional
+        )
 
     def get_cell_temperature(self, poa_global, temp_air, wind_speed, model,
                              effective_irradiance=None):
@@ -2379,8 +2395,23 @@ def sapm_effective_irradiance(poa_direct, poa_diffuse, airmass_absolute, aoi,
     pvlib.pvsystem.sapm
     """
 
-    F1 = spectrum.spectral_factor_sapm(airmass_absolute, module)
-    F2 = iam.sapm(aoi, module)
+    F1 = spectrum.spectral_factor_sapm(
+        airmass_absolute,
+        module["A0"],
+        module["A1"],
+        module["A2"],
+        module["A3"],
+        module["A4"],
+    )
+    F2 = iam.sapm(
+        aoi,
+        module["B0"],
+        module["B1"],
+        module["B2"],
+        module["B3"],
+        module["B4"],
+        module["B5"],
+    )
 
     Ee = F1 * (poa_direct * F2 + module['FD'] * poa_diffuse)
 
diff --git a/pvlib/spectrum/mismatch.py b/pvlib/spectrum/mismatch.py
index 3afc210e73..dfde1a5a53 100644
--- a/pvlib/spectrum/mismatch.py
+++ b/pvlib/spectrum/mismatch.py
@@ -295,20 +295,19 @@ def spectral_factor_firstsolar(precipitable_water, airmass_absolute,
     return modifier
 
 
-def spectral_factor_sapm(airmass_absolute, module):
+def spectral_factor_sapm(airmass_absolute, A0, A1, A2, A3, A4):
     """
-    Calculates the spectral mismatch factor, :math:`f_1`,
-    using the Sandia Array Performance Model approach.
+    Calculates the spectral mismatch factor, :math:`f_1`, using the Sandia
+    Array Performance Model (SAPM) approach.
 
-    The SAPM spectral factor function is part of the broader Sandia Array
+    The SAPM incidence angle modifier is part of the broader Sandia Array
     Performance Model, which defines five points on an IV curve using empirical
     module-specific coefficients. Module coefficients for the SAPM are
-    available in the SAPM database and can be retrieved for use in the
-    ``module`` parameter through
-    :py:func:`pvlib.pvsystem.retrieve_sam()`. More details on the
-    SAPM can be found in [1]_, while a full description of the procedure to
-    determine the empirical model coefficients, including those for the SAPM
-    spectral correction, can be found in [2]_.
+    available in the SAPM database and can be retrieved for use through
+    :py:func:`pvlib.pvsystem.retrieve_sam()`. More details on the SAPM can be
+    found in [1]_, while a full description of the procedure to determine the
+    empirical model coefficients, including those for the SAPM spectral
+    correction, can be found in [2]_.
 
     Parameters
     ----------
@@ -317,10 +316,20 @@ def spectral_factor_sapm(airmass_absolute, module):
 
         Note: ``np.nan`` airmass values will result in 0 output.
 
-    module : dict-like
-        A dict, Series, or DataFrame defining the SAPM parameters.
-        Must contain keys `'A0'` through `'A4'`.
-        See the :py:func:`pvlib.pvsystem.sapm` notes section for more details.
+    A0 : float
+        The coefficient of the degree-0 polynomial term.
+
+    A1 : float
+        The coefficient of the degree-1 polynomial term.
+
+    A2 : float
+        The coefficient of the degree-2 polynomial term.
+
+    A3 : float
+        The coefficient of the degree-3 polynomial term.
+
+    A4 : float
+        The coefficient of the degree-4 polynomial term.
 
     Returns
     -------
@@ -330,7 +339,7 @@ def spectral_factor_sapm(airmass_absolute, module):
     Notes
     -----
     The SAPM spectral correction functions parameterises :math:`f_1` as a
-    fourth order polynomial function of absolute air mass:
+    fourth-order polynomial function of absolute air mass:
 
     .. math::
 
@@ -361,10 +370,7 @@ def spectral_factor_sapm(airmass_absolute, module):
 
     """
 
-    am_coeff = [module['A4'], module['A3'], module['A2'], module['A1'],
-                module['A0']]
-
-    spectral_loss = np.polyval(am_coeff, airmass_absolute)
+    spectral_loss = np.polyval([A4, A3, A2, A1, A0], airmass_absolute)
 
     spectral_loss = np.where(np.isnan(spectral_loss), 0, spectral_loss)
 
diff --git a/pvlib/tests/spectrum/test_mismatch.py b/pvlib/tests/spectrum/test_mismatch.py
index 5397a81f46..992462b5c6 100644
--- a/pvlib/tests/spectrum/test_mismatch.py
+++ b/pvlib/tests/spectrum/test_mismatch.py
@@ -148,7 +148,14 @@ def test_spectral_factor_firstsolar_range():
 ])
 def test_spectral_factor_sapm(sapm_module_params, airmass, expected):
 
-    out = spectrum.spectral_factor_sapm(airmass, sapm_module_params)
+    out = spectrum.spectral_factor_sapm(
+        airmass,
+        sapm_module_params["A0"],
+        sapm_module_params["A1"],
+        sapm_module_params["A2"],
+        sapm_module_params["A3"],
+        sapm_module_params["A4"],
+    )
 
     if isinstance(airmass, pd.Series):
         assert_series_equal(out, expected, check_less_precise=4)
diff --git a/pvlib/tests/test_iam.py b/pvlib/tests/test_iam.py
index f5ca231bd4..eb76cf5984 100644
--- a/pvlib/tests/test_iam.py
+++ b/pvlib/tests/test_iam.py
@@ -3,15 +3,46 @@
 
 @author: cwhanse
 """
+import inspect
 
 import numpy as np
+from numpy.testing import assert_allclose
 import pandas as pd
-
 import pytest
-from .conftest import assert_series_equal
-from numpy.testing import assert_allclose
+import scipy.interpolate
 
 from pvlib import iam as _iam
+from pvlib.tests.conftest import assert_series_equal
+
+
+def test_get_builtin_models():
+    builtin_models = _iam._get_builtin_models()
+
+    models = set(builtin_models.keys())
+    models_expected = {
+        "ashrae", "interp", "martin_ruiz", "physical", "sapm", "schlick"
+    }
+    assert set(models) == models_expected
+
+    for model in models:
+        builtin_model = builtin_models[model]
+
+        params_required_expected = set(
+            k for k, v in inspect.signature(
+                builtin_model["func"]
+            ).parameters.items() if v.default is inspect.Parameter.empty
+        )
+        assert builtin_model["params_required"].union(
+            {"aoi"}
+        ) == params_required_expected, model
+
+        params_optional_expected = set(
+            k for k, v in inspect.signature(
+                builtin_model["func"]
+            ).parameters.items() if v.default is not inspect.Parameter.empty
+        )
+        assert builtin_model["params_optional"] == \
+            params_optional_expected, model
 
 
 def test_ashrae():
@@ -222,7 +253,15 @@ def test_iam_interp():
 ])
 def test_sapm(sapm_module_params, aoi, expected):
 
-    out = _iam.sapm(aoi, sapm_module_params)
+    out = _iam.sapm(
+        aoi,
+        sapm_module_params["B0"],
+        sapm_module_params["B1"],
+        sapm_module_params["B2"],
+        sapm_module_params["B3"],
+        sapm_module_params["B4"],
+        sapm_module_params["B5"],
+    )
 
     if isinstance(aoi, pd.Series):
         assert_series_equal(out, expected, check_less_precise=4)
@@ -232,13 +271,38 @@ def test_sapm(sapm_module_params, aoi, expected):
 
 def test_sapm_limits():
     module_parameters = {'B0': 5, 'B1': 0, 'B2': 0, 'B3': 0, 'B4': 0, 'B5': 0}
-    assert _iam.sapm(1, module_parameters) == 5
+    assert _iam.sapm(
+        1,
+        module_parameters["B0"],
+        module_parameters["B1"],
+        module_parameters["B2"],
+        module_parameters["B3"],
+        module_parameters["B4"],
+        module_parameters["B5"],
+    ) == 5
 
     module_parameters = {'B0': 5, 'B1': 0, 'B2': 0, 'B3': 0, 'B4': 0, 'B5': 0}
-    assert _iam.sapm(1, module_parameters, upper=1) == 1
+    assert _iam.sapm(
+        1,
+        module_parameters["B0"],
+        module_parameters["B1"],
+        module_parameters["B2"],
+        module_parameters["B3"],
+        module_parameters["B4"],
+        module_parameters["B5"],
+        upper=1,
+    ) == 1
 
     module_parameters = {'B0': -5, 'B1': 0, 'B2': 0, 'B3': 0, 'B4': 0, 'B5': 0}
-    assert _iam.sapm(1, module_parameters) == 0
+    assert _iam.sapm(
+        1,
+        module_parameters["B0"],
+        module_parameters["B1"],
+        module_parameters["B2"],
+        module_parameters["B3"],
+        module_parameters["B4"],
+        module_parameters["B5"],
+    ) == 0
 
 
 def test_marion_diffuse_model(mocker):
@@ -265,10 +329,14 @@ def test_marion_diffuse_model(mocker):
     assert physical_spy.call_count == 3
 
     for k, v in ashrae_expected.items():
-        assert_allclose(ashrae_actual[k], v)
+        assert_allclose(
+            ashrae_actual[k], v, err_msg=f"ashrae component {k}"
+        )
 
     for k, v in physical_expected.items():
-        assert_allclose(physical_actual[k], v)
+        assert_allclose(
+            physical_actual[k], v, err_msg=f"physical component {k}"
+        )
 
 
 def test_marion_diffuse_kwargs():
@@ -281,11 +349,73 @@ def test_marion_diffuse_kwargs():
     actual = _iam.marion_diffuse('ashrae', 20, b=0.04)
 
     for k, v in expected.items():
-        assert_allclose(actual[k], v)
+        assert_allclose(actual[k], v, err_msg=f"component {k}")
+
+
+# IEC 61853-2 measurement data for tests.
+AOI_DATA = np.array([
+    0, 10, 20, 30, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90,
+])
+IAM_DATA = np.array([
+    1.0000, 1.0000, 1.0014, 1.0002, 0.9984, 0.9941, 0.9911, 0.9815, 0.9631,
+    0.9352, 0.8922, 0.8134, 0.6778, 0.4541, 0.0000,
+])
+
+
+def test_marion_diffuse_iam_function_without_kwargs():
+    """
+    Test PCHIP-interpolated custom IAM function from an IEC 61853-2 IAM
+    measurement, without any kwargs and with array input.
+    """
+    pchip = scipy.interpolate.PchipInterpolator(
+        AOI_DATA, IAM_DATA, extrapolate=False
+    )
+
+    # Custom IAM function without kwargs, using IEC 61853-2 measurement data.
+    def iam_pchip(aoi):
+        """
+        Compute unitless incident-angle modifier as a function of aoi, in
+        degrees.
+        """
+        iam_ = pchip(aoi)
+        iam_[90 < aoi] = 0.0
+
+        return iam_
+
+    expected = {
+        'sky': np.array([0.95671526, 0.96967113, 0.95672627, 0.88137573]),
+        'horizon': np.array([0.03718587, 0.94953826, 0.97722834, 0.94862772]),
+        'ground': np.array([0., 0.88137573, 0.95672627, 0.96967113]),
+    }
+    actual = _iam.marion_diffuse(iam_pchip, np.array([0.0, 45, 90, 135]))
+
+    for k, v in expected.items():
+        assert_allclose(actual[k], v, rtol=2e-7, err_msg=f"component {k}")
+
+
+def test_marion_diffuse_iam_with_kwargs():
+    """
+    Test custom IAM function (using iam.interp) with kwargs and scalar input.
+    """
+    expected = {
+        'sky': 0.9688222974371822,
+        'horizon': 0.94824614710175,
+        'ground': 0.878266931831978,
+    }
+    actual = _iam.marion_diffuse(
+        _iam.interp,
+        45.0,
+        theta_ref=AOI_DATA,
+        iam_ref=IAM_DATA,
+        normalize=False,
+    )
+
+    for k, v in expected.items():
+        assert_allclose(actual[k], v, err_msg=f"component {k}")
 
 
 def test_marion_diffuse_invalid():
-    with pytest.raises(ValueError):
+    with pytest.raises(ValueError, match='model must be one of: '):
         _iam.marion_diffuse('not_a_model', 20)
 
 
@@ -481,7 +611,9 @@ def scaled_weight(aoi):
 
 
 def test_convert_model_not_implemented():
-    with pytest.raises(NotImplementedError, match='model has not been'):
+    with pytest.raises(
+        NotImplementedError, match='No implementation for model foo'
+    ):
         _iam.convert('ashrae', {'b': 0.1}, 'foo')
 
 
@@ -532,7 +664,9 @@ def scaled_weight(aoi):
 
 
 def test_fit_model_not_implemented():
-    with pytest.raises(NotImplementedError, match='model has not been'):
+    with pytest.raises(
+        NotImplementedError, match='No implementation for model foo'
+    ):
         _iam.fit(np.array([0, 10]), np.array([1, 0.99]), 'foo')
 
 
diff --git a/pvlib/tests/test_modelchain.py b/pvlib/tests/test_modelchain.py
index dcbd820f16..e37aebfd99 100644
--- a/pvlib/tests/test_modelchain.py
+++ b/pvlib/tests/test_modelchain.py
@@ -2,17 +2,14 @@
 
 import numpy as np
 import pandas as pd
+import pytest
 
 from pvlib import iam, modelchain, pvsystem, temperature, inverter
 from pvlib.modelchain import ModelChain
 from pvlib.pvsystem import PVSystem
 from pvlib.location import Location
-from pvlib._deprecation import pvlibDeprecationWarning
 
 from .conftest import assert_series_equal, assert_frame_equal
-import pytest
-
-from .conftest import fail_on_pvlib_version
 
 
 @pytest.fixture(scope='function')
@@ -1428,7 +1425,7 @@ def constant_aoi_loss(mc):
 
 
 @pytest.mark.parametrize('aoi_model', [
-    'sapm', 'ashrae', 'physical', 'martin_ruiz'
+    'ashrae', 'martin_ruiz', 'physical', 'sapm', 'schlick'  # not interp
 ])
 def test_aoi_models(sapm_dc_snl_ac_system, location, aoi_model,
                     weather, mocker):
@@ -1444,7 +1441,7 @@ def test_aoi_models(sapm_dc_snl_ac_system, location, aoi_model,
 
 
 @pytest.mark.parametrize('aoi_model', [
-    'sapm', 'ashrae', 'physical', 'martin_ruiz'
+    'ashrae', 'martin_ruiz', 'physical', 'sapm', 'schlick'  # not interp
 ])
 def test_aoi_models_singleon_weather_single_array(
         sapm_dc_snl_ac_system, location, aoi_model, weather):
@@ -1502,27 +1499,66 @@ def test_aoi_model_user_func(sapm_dc_snl_ac_system, location, weather, mocker):
     assert mc.results.ac.iloc[1] < 1
 
 
-@pytest.mark.parametrize('aoi_model', [
-    'sapm', 'ashrae', 'physical', 'martin_ruiz', 'interp'
-])
+@pytest.mark.parametrize(
+    'aoi_model',
+    # "schlick" omitted, cannot be distinguished from "no_loss" AOI model.
+    [
+        {'params': {'b': 0.025}},
+        {'params': {'iam_ref': (1., 0.85), 'theta_ref': (0., 80.)}},
+        {'params': {'a_r': 0.1}},
+        {'params': {'n': 1.5}},
+        {
+            'params': {
+                'B0': 1,
+                'B1': -0.002438,
+                'B2': 0.0003103,
+                'B3': -0.00001246,
+                'B4': 2.11E-07,
+                'B5': -1.36E-09,
+            }
+        },
+    ],
+    ids=['ashrae', 'interp', 'martin_ruiz', 'physical', 'sapm'],
+)
 def test_infer_aoi_model(location, system_no_aoi, aoi_model):
-    for k in iam._IAM_MODEL_PARAMS[aoi_model]:
-        system_no_aoi.arrays[0].module_parameters.update({k: 1.0})
+    system_no_aoi.arrays[0].module_parameters.update(aoi_model["params"])
     mc = ModelChain(system_no_aoi, location, spectral_model='no_loss')
     assert isinstance(mc, ModelChain)
 
 
-@pytest.mark.parametrize('aoi_model,model_kwargs', [
-    # model_kwargs has both required and optional kwargs; test all
-    ('physical',
-     {'n': 1.526, 'K': 4.0, 'L': 0.002,  # required
-      'n_ar': 1.8}),  # extra
-    ('interp',
-     {'theta_ref': (0, 75, 85, 90), 'iam_ref': (1, 0.8, 0.42, 0),  # required
-      'method': 'cubic', 'normalize': False})])  # extra
-def test_infer_aoi_model_with_extra_params(location, system_no_aoi, aoi_model,
-                                           model_kwargs, weather, mocker):
-    # test extra parameters not defined at iam._IAM_MODEL_PARAMS are passed
+@pytest.mark.parametrize(
+    'aoi_model,model_kwargs',
+    [
+        (
+            'sapm',
+            {
+                # required
+                'B0': 1,
+                'B1': -0.002438,
+                'B2': 0.0003103,
+                'B3': -0.00001246,
+                'B4': 2.11E-07,
+                'B5': -1.36E-09,
+                # optional
+                'upper': 1.,
+            }
+        ),
+        (
+            'interp',
+            {
+                # required
+                'theta_ref': (0, 75, 85, 90),
+                'iam_ref': (1, 0.8, 0.42, 0),
+                # optional
+                'method': 'cubic',
+                'normalize': False,
+            }
+        )
+    ]
+)
+def test_infer_aoi_model_with_required_and_optional_params(
+    location, system_no_aoi, aoi_model, model_kwargs, weather, mocker
+):
     m = mocker.spy(iam, aoi_model)
     system_no_aoi.arrays[0].module_parameters.update(**model_kwargs)
     mc = ModelChain(system_no_aoi, location, spectral_model='no_loss')
@@ -1533,8 +1569,7 @@ def test_infer_aoi_model_with_extra_params(location, system_no_aoi, aoi_model,
 
 
 def test_infer_aoi_model_invalid(location, system_no_aoi):
-    exc_text = 'could not infer AOI model'
-    with pytest.raises(ValueError, match=exc_text):
+    with pytest.raises(ValueError, match='could not infer AOI model'):
         ModelChain(system_no_aoi, location, spectral_model='no_loss')
 
 
diff --git a/pvlib/tests/test_pvsystem.py b/pvlib/tests/test_pvsystem.py
index fd482c5127..9372234aff 100644
--- a/pvlib/tests/test_pvsystem.py
+++ b/pvlib/tests/test_pvsystem.py
@@ -25,16 +25,56 @@
 
 @pytest.mark.parametrize('iam_model,model_params', [
     ('ashrae', {'b': 0.05}),
-    ('physical', {'K': 4, 'L': 0.002, 'n': 1.526}),
+    ('interp', {'iam_ref': (1., 0.85), 'theta_ref': (0., 80.)}),
     ('martin_ruiz', {'a_r': 0.16}),
+    ('physical', {'K': 4, 'L': 0.002, 'n': 1.526}),
+    (
+        'sapm',
+        {
+            k: v for k, v in pvsystem.retrieve_sam(
+                'SandiaMod')['Canadian_Solar_CS5P_220M___2009_'].items()
+            if k in ('B0', 'B1', 'B2', 'B3', 'B4', 'B5')
+        }
+    ),
+    ('schlick', {}),
 ])
 def test_PVSystem_get_iam(mocker, iam_model, model_params):
     m = mocker.spy(_iam, iam_model)
     system = pvsystem.PVSystem(module_parameters=model_params)
-    thetas = 1
+    thetas = 45
     iam = system.get_iam(thetas, iam_model=iam_model)
     m.assert_called_with(thetas, **model_params)
-    assert iam < 1.
+    assert 0 < iam < 1
+
+
+def test_PVSystem_get_iam_unsupported_model():
+    iam_model = 'foobar'
+    system = pvsystem.PVSystem()
+
+    with pytest.raises(
+        ValueError, match=f'{iam_model} is not a valid iam_model'
+    ):
+        system.get_iam(45, iam_model=iam_model)
+
+
+@pytest.mark.parametrize('iam_model,model_params', [
+    ('interp', {'iam_ref': (1., 0.85)}),  # missing theta_ref
+    (
+        'sapm',
+        {
+            k: v for k, v in pvsystem.retrieve_sam(
+                'SandiaMod')['Canadian_Solar_CS5P_220M___2009_'].items()
+            if k in ('B0', 'B1', 'B2', 'B3', 'B4')  # missing B5
+        }
+    ),
+])
+def test_PVSystem_get_iam_missing_required_param(iam_model, model_params):
+    system = pvsystem.PVSystem(module_parameters=model_params)
+
+    with pytest.raises(
+        KeyError, match="is missing in module_parameters"
+    ):
+        system.get_iam(45, iam_model=iam_model)
 
 
 def test_PVSystem_multi_array_get_iam():
@@ -58,7 +98,15 @@ def test_PVSystem_get_iam_sapm(sapm_module_params, mocker):
     mocker.spy(_iam, 'sapm')
     aoi = 0
     out = system.get_iam(aoi, 'sapm')
-    _iam.sapm.assert_called_once_with(aoi, sapm_module_params)
+    _iam.sapm.assert_called_once_with(
+        aoi,
+        sapm_module_params["B0"],
+        sapm_module_params["B1"],
+        sapm_module_params["B2"],
+        sapm_module_params["B3"],
+        sapm_module_params["B4"],
+        sapm_module_params["B5"],
+    )
     assert_allclose(out, 1.0, atol=0.01)
 
 
@@ -236,7 +284,14 @@ def test_PVSystem_multi_array_sapm(sapm_module_params):
 def test_sapm_spectral_loss_deprecated(sapm_module_params):
     with pytest.warns(pvlibDeprecationWarning,
                       match='Use pvlib.spectrum.spectral_factor_sapm'):
-        pvsystem.sapm_spectral_loss(1, sapm_module_params)
+        pvsystem.sapm_spectral_loss(
+            1,
+            sapm_module_params["A0"],
+            sapm_module_params["A1"],
+            sapm_module_params["A2"],
+            sapm_module_params["A3"],
+            sapm_module_params["A4"],
+        )
 
 
 def test_PVSystem_sapm_spectral_loss(sapm_module_params, mocker):
@@ -244,8 +299,14 @@ def test_PVSystem_sapm_spectral_loss(sapm_module_params, mocker):
     system = pvsystem.PVSystem(module_parameters=sapm_module_params)
     airmass = 2
     out = system.sapm_spectral_loss(airmass)
-    spectrum.spectral_factor_sapm.assert_called_once_with(airmass,
-                                                          sapm_module_params)
+    spectrum.spectral_factor_sapm.assert_called_once_with(
+        airmass,
+        sapm_module_params["A0"],
+        sapm_module_params["A1"],
+        sapm_module_params["A2"],
+        sapm_module_params["A3"],
+        sapm_module_params["A4"],
+    )
     assert_allclose(out, 1, atol=0.5)
 
 
diff --git a/pyproject.toml b/pyproject.toml
index 8464b7ce23..d1583421e1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -9,7 +9,7 @@ description = "A set of functions and classes for simulating the performance of
 authors = [
     { name = "pvlib python Developers", email = "pvlib-admin@googlegroups.com" },
 ]
-requires-python = ">=3.9"
+requires-python = ">=3.9, <3.13"
 dependencies = [
     'numpy >= 1.19.3',
     'pandas >= 1.3.0',
@@ -49,9 +49,9 @@ dynamic = ["version"]
 optional = [
     'cython',
     'ephem',
-    'nrel-pysam',
-    'numba >= 0.17.0',
-    'solarfactors',
+    'nrel-pysam; python_version < "3.13"',
+    'numba >= 0.17.0; python_version < "3.13"',
+    'solarfactors; python_version < "3.12"',
     'statsmodels',
 ]
 doc = [
@@ -65,10 +65,11 @@ doc = [
     'pillow',
     'sphinx-toggleprompt == 0.5.2',
     'sphinx-favicon',
-    'solarfactors',
+    'solarfactors; python_version < "3.12"',
     'sphinx-hoverxref',
 ]
 test = [
+    'flake8==5.0.4',  # Should match version in github workflow flake8.yml.
     'pytest',
     'pytest-cov',
     'pytest-mock',