diff --git a/docs/sphinx/source/reference/pv_modeling/parameters.rst b/docs/sphinx/source/reference/pv_modeling/parameters.rst
index 9b1817bd01..1f146eaccc 100644
--- a/docs/sphinx/source/reference/pv_modeling/parameters.rst
+++ b/docs/sphinx/source/reference/pv_modeling/parameters.rst
@@ -12,6 +12,7 @@ Functions for fitting single diode models
     ivtools.sdm.fit_cec_sam
     ivtools.sdm.fit_desoto
     ivtools.sdm.fit_pvsyst_sandia
+    ivtools.sdm.fit_pvsyst_iec61853_sandia_2025
     ivtools.sdm.fit_desoto_sandia
 
 Functions for fitting the single diode equation
diff --git a/docs/sphinx/source/whatsnew/v0.12.1.rst b/docs/sphinx/source/whatsnew/v0.12.1.rst
index f1c9b7e69d..b7fffd5996 100644
--- a/docs/sphinx/source/whatsnew/v0.12.1.rst
+++ b/docs/sphinx/source/whatsnew/v0.12.1.rst
@@ -14,8 +14,9 @@ Bug fixes
 
 Enhancements
 ~~~~~~~~~~~~
-* ``pvlib.ivtools.sdm`` is now a subpackage. (:issue:`2252`, :pull:`2256`)
-
+* :py:mod:`pvlib.ivtools.sdm` is now a subpackage. (:issue:`2252`, :pull:`2256`)
+* Add a function for estimating PVsyst SDM parameters from IEC 61853-1 matrix
+  data (:py:func:`~pvlib.ivtools.sdm.fit_pvsyst_iec61853_sandia_2025`). (:issue:`2185`, :pull:`2429`)
 
 Documentation
 ~~~~~~~~~~~~~
diff --git a/pvlib/ivtools/sdm/__init__.py b/pvlib/ivtools/sdm/__init__.py
index b45519fc16..8535f1f5e6 100644
--- a/pvlib/ivtools/sdm/__init__.py
+++ b/pvlib/ivtools/sdm/__init__.py
@@ -15,5 +15,6 @@
 
 from pvlib.ivtools.sdm.pvsyst import (  # noqa: F401
     fit_pvsyst_sandia,
+    fit_pvsyst_iec61853_sandia_2025,
     pvsyst_temperature_coeff,
 )
diff --git a/pvlib/ivtools/sdm/pvsyst.py b/pvlib/ivtools/sdm/pvsyst.py
index 6233c38c72..7e1f8c71ab 100644
--- a/pvlib/ivtools/sdm/pvsyst.py
+++ b/pvlib/ivtools/sdm/pvsyst.py
@@ -14,6 +14,9 @@
     _extract_sdm_params, _initial_iv_params, _update_iv_params
 )
 
+from pvlib.pvsystem import (
+    _pvsyst_Rsh, _pvsyst_IL, _pvsyst_Io, _pvsyst_nNsVth, _pvsyst_gamma
+)
 
 CONSTANTS = {'E0': 1000.0, 'T0': 25.0, 'k': constants.k, 'q': constants.e}
 
@@ -305,3 +308,323 @@ def maxp(temp_cell, irrad_ref, alpha_sc, gamma_ref, mu_gamma, I_L_ref,
     gamma_pdc = _first_order_centered_difference(maxp, x0=temp_ref, args=args)
 
     return gamma_pdc / pmp
+
+
+def fit_pvsyst_iec61853_sandia_2025(effective_irradiance, temp_cell,
+                                    i_sc, v_oc, i_mp, v_mp,
+                                    cells_in_series, EgRef=1.121,
+                                    alpha_sc=None, beta_mp=None,
+                                    R_s=None, r_sh_coeff=0.12,
+                                    min_Rsh_irradiance=None,
+                                    irradiance_tolerance=20,
+                                    temperature_tolerance=1):
+    """
+    Estimate parameters for the PVsyst module performance model using
+    IEC 61853-1 matrix measurements.
+
+    Parameters
+    ----------
+    effective_irradiance : array
+        Effective irradiance for each test condition [W/m²]
+    temp_cell : array
+        Cell temperature for each test condition [C]
+    i_sc : array
+        Short circuit current for each test condition [A]
+    v_oc : array
+        Open circuit voltage for each test condition [V]
+    i_mp : array
+        Current at maximum power point for each test condition [A]
+    v_mp : array
+        Voltage at maximum power point for each test condition [V]
+    cells_in_series : int
+        The number of cells connected in series.
+    EgRef : float, optional
+        The energy bandgap at reference temperature in units of eV.
+        1.121 eV for crystalline silicon. EgRef must be >0.
+    alpha_sc : float, optional
+        Temperature coefficient of short circuit current.  If not specified,
+        it will be estimated using the ``i_sc`` values at irradiance of
+        1000 W/m2. [A/K]
+    beta_mp : float, optional
+        Temperature coefficient of maximum power voltage.  If not specified,
+        it will be estimated using the ``v_mp`` values at irradiance of
+        1000 W/m2. [1/K]
+    R_s : float, optional
+        Series resistance value.  If not provided, a value will be estimated
+        from the input measurements. [ohm]
+    r_sh_coeff : float, default 0.12
+        Shunt resistance fitting coefficient.  The default value is taken
+        from [1]_.
+    min_Rsh_irradiance : float, optional
+        Irradiance threshold below which values are excluded when estimating
+        shunt resistance parameter values.  May be useful for modules
+        with problematic low-light measurements. [W/m²]
+    irradiance_tolerance : float, default 20
+        Tolerance for irradiance variation around the STC value.
+        The default value corresponds to a +/- 2% interval around the STC
+        value of 1000 W/m². [W/m²]
+    temperature_tolerance : float, default 1
+        Tolerance for temperature variation around the STC value.
+        The default value corresponds to a +/- 1 degree interval around the STC
+        value of 25 degrees. [C]
+
+    Returns
+    -------
+    dict
+        alpha_sc : float
+            short circuit current temperature coefficient [A/K]
+        gamma_ref : float
+            diode (ideality) factor at STC [unitless]
+        mu_gamma : float
+            temperature coefficient for diode (ideality) factor [1/K]
+        I_L_ref : float
+            light current at STC [A]
+        I_o_ref : float
+            dark current at STC [A]
+        R_sh_ref : float
+            shunt resistance at STC [ohm]
+        R_sh_0 : float
+            shunt resistance at zero irradiance [ohm]
+        R_sh_exp : float
+            exponential factor defining decrease in shunt resistance with
+            increasing effective irradiance
+        R_s : float
+            series resistance at STC [ohm]
+        cells_in_series : int
+            number of cells in series
+        EgRef : float
+            effective band gap at STC [eV]
+
+    See also
+    --------
+    pvlib.pvsystem.calcparams_pvsyst
+    pvlib.ivtools.sdm.fit_pvsyst_sandia
+
+    Notes
+    -----
+    Input arrays of operating conditions and electrical measurements must be
+    1-D with equal lengths.
+
+    Values supplied for ``alpha_sc``, ``beta_mp``, and ``R_s`` must be
+    consistent with the matrix data, as these values are used when estimating
+    other model parameters.
+
+    This method is non-iterative.  In some cases, it may be desirable to
+    refine the estimated parameter values using a numerical optimizer such as
+    the default method in ``scipy.optimize.minimize``.
+
+    References
+    ----------
+    .. [1] K. S. Anderson, C. W. Hansen, and M. Theristis, "A Noniterative
+       Method of Estimating Parameter Values for the PVsyst Version 6
+       Single-Diode Model From IEC 61853-1 Matrix Measurements," IEEE Journal
+       of Photovoltaics, vol. 15, 3, 2025. :doi:`10.1109/JPHOTOV.2025.3554338`
+    """
+
+    is_g_stc = np.isclose(effective_irradiance, 1000, rtol=0,
+                          atol=irradiance_tolerance)
+    is_t_stc = np.isclose(temp_cell, 25, rtol=0,
+                          atol=temperature_tolerance)
+
+    if alpha_sc is None:
+        mu_i_sc = _fit_tempco_pvsyst_iec61853_sandia_2025(i_sc[is_g_stc],
+                                                          temp_cell[is_g_stc])
+        i_sc_ref = float(i_sc[is_g_stc & is_t_stc].item())
+        alpha_sc = mu_i_sc * i_sc_ref
+
+    if beta_mp is None:
+        beta_mp = _fit_tempco_pvsyst_iec61853_sandia_2025(v_mp[is_g_stc],
+                                                          temp_cell[is_g_stc])
+
+    R_sh_ref, R_sh_0, R_sh_exp = \
+        _fit_shunt_resistances_pvsyst_iec61853_sandia_2025(
+            i_sc, i_mp, v_mp, effective_irradiance, temp_cell, beta_mp,
+            coeff=r_sh_coeff, min_irradiance=min_Rsh_irradiance)
+
+    if R_s is None:
+        R_s = _fit_series_resistance_pvsyst_iec61853_sandia_2025(v_oc, i_mp,
+                                                                 v_mp)
+
+    gamma_ref, mu_gamma = \
+        _fit_diode_ideality_factor_pvsyst_iec61853_sandia_2025(
+            i_sc[is_t_stc], v_oc[is_t_stc], i_mp[is_t_stc], v_mp[is_t_stc],
+            effective_irradiance[is_t_stc], temp_cell[is_t_stc],
+            R_sh_ref, R_sh_0, R_sh_exp, R_s, cells_in_series)
+
+    I_o_ref = _fit_saturation_current_pvsyst_iec61853_sandia_2025(
+        i_sc, v_oc, effective_irradiance, temp_cell, gamma_ref, mu_gamma,
+        R_sh_ref, R_sh_0, R_sh_exp, R_s, cells_in_series, EgRef
+    )
+
+    I_L_ref = _fit_photocurrent_pvsyst_iec61853_sandia_2025(
+        i_sc, effective_irradiance, temp_cell, alpha_sc,
+        gamma_ref, mu_gamma,
+        I_o_ref, R_sh_ref, R_sh_0, R_sh_exp, R_s, cells_in_series, EgRef
+    )
+
+    gamma_ref, mu_gamma = \
+        _fit_diode_ideality_factor_post_pvsyst_iec61853_sandia_2025(
+            i_mp, v_mp, effective_irradiance, temp_cell, alpha_sc, I_L_ref,
+            I_o_ref, R_sh_ref, R_sh_0, R_sh_exp, R_s, cells_in_series, EgRef)
+
+    fitted_params = dict(
+        alpha_sc=alpha_sc,
+        gamma_ref=gamma_ref,
+        mu_gamma=mu_gamma,
+        I_L_ref=I_L_ref,
+        I_o_ref=I_o_ref,
+        R_sh_ref=R_sh_ref,
+        R_sh_0=R_sh_0,
+        R_sh_exp=R_sh_exp,
+        R_s=R_s,
+        cells_in_series=cells_in_series,
+        EgRef=EgRef,
+    )
+    return fitted_params
+
+
+def _fit_tempco_pvsyst_iec61853_sandia_2025(values, temp_cell,
+                                            temp_cell_ref=25):
+    fit = np.polynomial.polynomial.Polynomial.fit(temp_cell, values, deg=1)
+    intercept, slope = fit.convert().coef
+    value_ref = intercept + slope*temp_cell_ref
+    return slope / value_ref
+
+
+def _fit_shunt_resistances_pvsyst_iec61853_sandia_2025(
+        i_sc, i_mp, v_mp, effective_irradiance, temp_cell,
+        beta_v_mp, coeff=0.2, min_irradiance=None):
+    if min_irradiance is None:
+        min_irradiance = 0
+
+    mask = effective_irradiance >= min_irradiance
+    i_sc = i_sc[mask]
+    i_mp = i_mp[mask]
+    v_mp = v_mp[mask]
+    effective_irradiance = effective_irradiance[mask]
+    temp_cell = temp_cell[mask]
+
+    # Equation 10
+    Rsh_est = (
+        (v_mp / (1 + beta_v_mp * (temp_cell - 25)))
+        / (coeff * (i_sc - i_mp))
+    )
+    Rshexp = 5.5
+
+    # Eq 11
+    y = Rsh_est
+    x = np.exp(-Rshexp * effective_irradiance / 1000)
+
+    fit = np.polynomial.polynomial.Polynomial.fit(x, y, deg=1)
+    intercept, slope = fit.convert().coef
+    Rshbase = intercept
+    Rsh0 = slope + Rshbase
+
+    # Eq 12
+    expRshexp = np.exp(-Rshexp)
+    Rshref = Rshbase * (1 - expRshexp) + Rsh0 * expRshexp
+
+    return Rshref, Rsh0, Rshexp
+
+
+def _fit_series_resistance_pvsyst_iec61853_sandia_2025(v_oc, i_mp, v_mp):
+    # Stein et al 2014, https://doi.org/10.1109/PVSC.2014.6925326
+
+    # Eq 13
+    x = np.array([np.ones(len(i_mp)), i_mp, np.log(i_mp), v_mp]).T
+    y = v_oc
+
+    coeff, _, _, _ = np.linalg.lstsq(x, y, rcond=None)
+    R_s = coeff[1]
+    return R_s
+
+
+def _fit_diode_ideality_factor_pvsyst_iec61853_sandia_2025(
+        i_sc, v_oc, i_mp, v_mp, effective_irradiance, temp_cell,
+        R_sh_ref, R_sh_0, R_sh_exp, R_s, cells_in_series):
+
+    NsVth = _pvsyst_nNsVth(temp_cell, gamma=1, cells_in_series=cells_in_series)
+    Rsh = _pvsyst_Rsh(effective_irradiance, R_sh_ref, R_sh_0, R_sh_exp)
+    term1 = (i_sc * (1 + R_s/Rsh) - v_oc / Rsh)  # Eq 15
+    term2 = (i_sc - i_mp) * (1 + R_s/Rsh) - v_mp / Rsh  # Eq 16
+
+    # Eq 14
+    x1 = NsVth * np.log(term2 / term1)
+
+    x = np.array([x1]).T
+    y = v_mp + i_mp*R_s - v_oc
+
+    coeff, _, _, _ = np.linalg.lstsq(x, y, rcond=None)
+    gamma_ref = coeff[0]
+    return gamma_ref, 0
+
+
+def _fit_saturation_current_pvsyst_iec61853_sandia_2025(
+        i_sc, v_oc, effective_irradiance, temp_cell, gamma_ref, mu_gamma,
+        R_sh_ref, R_sh_0, R_sh_exp, R_s, cells_in_series, EgRef):
+    R_sh = _pvsyst_Rsh(effective_irradiance, R_sh_ref, R_sh_0, R_sh_exp)
+    gamma = _pvsyst_gamma(temp_cell, gamma_ref, mu_gamma)
+    nNsVth = _pvsyst_nNsVth(temp_cell, gamma, cells_in_series)
+
+    # Eq 17
+    I_o_est = (i_sc * (1 + R_s/R_sh) - v_oc/R_sh) / (np.expm1(v_oc / nNsVth))
+    x = _pvsyst_Io(temp_cell, gamma, I_o_ref=1, EgRef=EgRef)
+
+    # Eq 18
+    log_I_o_ref = np.mean(np.log(I_o_est) - np.log(x))
+    I_o_ref = np.exp(log_I_o_ref)
+
+    return I_o_ref
+
+
+def _fit_photocurrent_pvsyst_iec61853_sandia_2025(
+        i_sc, effective_irradiance, temp_cell, alpha_sc, gamma_ref,
+        mu_gamma, I_o_ref, R_sh_ref, R_sh_0, R_sh_exp, R_s, cells_in_series,
+        EgRef):
+    R_sh = _pvsyst_Rsh(effective_irradiance, R_sh_ref, R_sh_0, R_sh_exp)
+    gamma = _pvsyst_gamma(temp_cell, gamma_ref, mu_gamma)
+    I_o = _pvsyst_Io(temp_cell, gamma, I_o_ref, EgRef)
+    nNsVth = _pvsyst_nNsVth(temp_cell, gamma, cells_in_series)
+
+    # Eq 19
+    I_L_est = i_sc + I_o * (np.expm1(i_sc * R_s / nNsVth)) + i_sc * R_s / R_sh
+
+    # Eq 20
+    x = np.array([effective_irradiance / 1000]).T
+    y = I_L_est - effective_irradiance / 1000 * alpha_sc * (temp_cell - 25)
+    coeff, _, _, _ = np.linalg.lstsq(x, y, rcond=None)
+    I_L_ref = coeff[0]
+    return I_L_ref
+
+
+def _fit_diode_ideality_factor_post_pvsyst_iec61853_sandia_2025(
+        i_mp, v_mp, effective_irradiance, temp_cell, alpha_sc, I_L_ref,
+        I_o_ref, R_sh_ref, R_sh_0, R_sh_exp, R_s, cells_in_series, EgRef):
+
+    Rsh = _pvsyst_Rsh(effective_irradiance, R_sh_ref, R_sh_0, R_sh_exp)
+    I_L = _pvsyst_IL(effective_irradiance, temp_cell, I_L_ref, alpha_sc)
+    NsVth = _pvsyst_nNsVth(temp_cell, gamma=1, cells_in_series=cells_in_series)
+
+    Tref_K = 25 + 273.15
+    Tcell_K = temp_cell + 273.15
+
+    # Eq 21
+    k = constants.k  # Boltzmann constant in J/K
+    q = constants.e  # elementary charge in coulomb
+    numerator = (
+        (q * EgRef / k) * (1/Tref_K - 1/Tcell_K)
+        + (v_mp + i_mp*R_s) / NsVth
+    )
+    denominator = (
+        np.log((I_L - i_mp - (v_mp+i_mp*R_s) / Rsh) / I_o_ref)
+        - 3 * np.log(Tcell_K / Tref_K)
+    )
+    gamma_est = numerator / denominator
+
+    # Eq 22
+    x = np.array([np.ones(len(i_mp)), temp_cell - 25]).T
+    y = gamma_est
+
+    coeff, _, _, _ = np.linalg.lstsq(x, y, rcond=None)
+    gamma_ref, mu_gamma = coeff
+    return gamma_ref, mu_gamma
diff --git a/tests/ivtools/sdm/test_pvsyst.py b/tests/ivtools/sdm/test_pvsyst.py
index 6927a8c1d3..2a1793633b 100644
--- a/tests/ivtools/sdm/test_pvsyst.py
+++ b/tests/ivtools/sdm/test_pvsyst.py
@@ -7,6 +7,8 @@
 
 from tests.conftest import requires_statsmodels, TESTS_DATA_DIR
 
+import pytest
+
 
 def _read_iv_curves_for_test(datafile, npts):
     """ read constants and npts IV curves from datafile """
@@ -198,3 +200,138 @@ def test_pvsyst_temperature_coeff():
         params['I_L_ref'], params['I_o_ref'], params['R_sh_ref'],
         params['R_sh_0'], params['R_s'], params['cells_in_series'])
     assert_allclose(gamma_pdc, expected, rtol=0.0005)
+
+
+@pytest.fixture
+def pvsyst_iec61853_table3():
+    # test cases in Table 3 of https://doi.org/10.1109/JPHOTOV.2025.3554338
+    parameters = ['alpha_sc', 'I_L_ref', 'I_o_ref', 'gamma_ref', 'mu_gamma',
+                  'R_sh_ref', 'R_sh_0', 'R_sh_exp', 'R_s', 'cells_in_series',
+                  'EgRef']
+    high_current = {
+        'true': [1e-3, 10, 1e-11, 1.0, -1e-4, 300, 1500, 5.5, 0.2, 60, 1.121],
+        'estimated': [9.993e-4, 9.997, 1.408e-10, 1.109, -1.666e-5, 543.7,
+                      6130, 5.5, 0.176, 60, 1.121]
+    }
+    high_voltage = {
+        'true': [1e-3, 1.8, 1e-9, 1.5, 1e-4, 3000, 10000, 5.5, 5.0, 108, 1.4],
+        'estimated': [9.983e-4, 1.799, 1.225e-8, 1.706, 2.662e-4, 4288, 47560,
+                      5.5, 4.292, 108, 1.4],
+    }
+    for d in [high_current, high_voltage]:
+        for key in d:
+            d[key] = dict(zip(parameters, np.array(d[key])))
+
+    return {
+        'high current': high_current, 'high voltage': high_voltage
+    }
+
+
+@pytest.fixture
+def iec61853_conditions():
+    ee = np.array([100, 100, 200, 200, 400, 400, 400, 600, 600, 600, 600,
+                   800, 800, 800, 800, 1000, 1000, 1000, 1000, 1100, 1100,
+                   1100], dtype=np.float64)
+    tc = np.array([15, 25, 15, 25, 15, 25, 50, 15, 25, 50, 75,
+                   15, 25, 50, 75, 15, 25, 50, 75, 25, 50, 75],
+                  dtype=np.float64)
+    return ee, tc
+
+
+def test_fit_pvsyst_iec61853_sandia_2025(pvsyst_iec61853_table3,
+                                         iec61853_conditions):
+    ee, tc = iec61853_conditions
+    for _, case in pvsyst_iec61853_table3.items():
+        true_params = case['true']
+        expected_params = case['estimated']
+
+        sde_params = pvsystem.calcparams_pvsyst(ee, tc, **true_params)
+        iv = pvsystem.singlediode(*sde_params)
+
+        fitted_params = sdm.fit_pvsyst_iec61853_sandia_2025(
+            ee, tc, iv['i_sc'], iv['v_oc'], iv['i_mp'], iv['v_mp'],
+            true_params['cells_in_series'], EgRef=true_params['EgRef'],
+        )
+        for key in expected_params.keys():
+            assert np.isclose(fitted_params[key], expected_params[key],
+                              atol=0, rtol=1e-3)
+
+
+def test_fit_pvsyst_iec61853_sandia_2025_optional(pvsyst_iec61853_table3,
+                                                  iec61853_conditions):
+    # verify that modifying each optional parameter results in different output
+    ee, tc = iec61853_conditions
+    case = pvsyst_iec61853_table3['high current']
+    true_params = case['true']
+    expected = case['estimated']
+    sde_params = pvsystem.calcparams_pvsyst(ee, tc, **true_params)
+    iv = pvsystem.singlediode(*sde_params)
+
+    main_inputs = dict(
+        effective_irradiance=ee, temp_cell=tc, i_sc=iv['i_sc'],
+        v_oc=iv['v_oc'], i_mp=iv['i_mp'], v_mp=iv['v_mp'],
+        cells_in_series=true_params['cells_in_series']
+    )
+
+    fitted_params = sdm.fit_pvsyst_iec61853_sandia_2025(**main_inputs,
+                                                        EgRef=1.0)
+    assert not np.isclose(fitted_params['I_o_ref'], expected['I_o_ref'],
+                          atol=0, rtol=1e-3)
+
+    fitted_params = sdm.fit_pvsyst_iec61853_sandia_2025(**main_inputs,
+                                                        alpha_sc=0.5e-3)
+    assert not np.isclose(fitted_params['alpha_sc'], expected['alpha_sc'],
+                          atol=0, rtol=1e-3)
+
+    fitted_params = sdm.fit_pvsyst_iec61853_sandia_2025(**main_inputs,
+                                                        beta_mp=-1e-4)
+    assert not np.isclose(fitted_params['R_sh_ref'], expected['R_sh_ref'],
+                          atol=0, rtol=1e-3)
+
+    fitted_params = sdm.fit_pvsyst_iec61853_sandia_2025(**main_inputs,
+                                                        r_sh_coeff=0.3)
+    assert not np.isclose(fitted_params['R_sh_ref'], expected['R_sh_ref'],
+                          atol=0, rtol=1e-3)
+
+    fitted_params = sdm.fit_pvsyst_iec61853_sandia_2025(**main_inputs,
+                                                        R_s=0.5)
+    assert not np.isclose(fitted_params['R_s'], expected['R_s'],
+                          atol=0, rtol=1e-3)
+
+    fitted_params = sdm.fit_pvsyst_iec61853_sandia_2025(**main_inputs,
+                                                        min_Rsh_irradiance=500)
+    assert not np.isclose(fitted_params['R_sh_ref'], expected['R_sh_ref'],
+                          atol=0, rtol=1e-3)
+
+
+def test_fit_pvsyst_iec61853_sandia_2025_tolerance(pvsyst_iec61853_table3,
+                                                   iec61853_conditions):
+    # verify that the *_tolerance parameters allow non-"perfect" irradiance
+    # and temperature values
+    ee, tc = iec61853_conditions
+    ee[ee == 1000] = 999
+    tc[tc == 25] = 25.1
+    case = pvsyst_iec61853_table3['high current']
+    true_params = case['true']
+    expected = case['estimated']
+    sde_params = pvsystem.calcparams_pvsyst(ee, tc, **true_params)
+    iv = pvsystem.singlediode(*sde_params)
+
+    inputs = dict(
+        effective_irradiance=ee, temp_cell=tc, i_sc=iv['i_sc'],
+        v_oc=iv['v_oc'], i_mp=iv['i_mp'], v_mp=iv['v_mp'],
+        cells_in_series=true_params['cells_in_series']
+    )
+    fitted_params = sdm.fit_pvsyst_iec61853_sandia_2025(**inputs)
+    # still get approximately the expected values
+    for key in expected.keys():
+        assert np.isclose(fitted_params[key], expected[key], atol=0, rtol=1e-2)
+
+    # but if the changes exceed the specified tolerance, then error:
+    with pytest.raises(ValueError, match='Coefficient array is empty'):
+        sdm.fit_pvsyst_iec61853_sandia_2025(**inputs, irradiance_tolerance=0.1)
+
+    with pytest.raises(ValueError,
+                       match='can only convert an array of size 1'):
+        sdm.fit_pvsyst_iec61853_sandia_2025(**inputs,
+                                            temperature_tolerance=0.1)