diff --git a/docs/examples/iv-modeling/plot_villalva.py b/docs/examples/iv-modeling/plot_villalva.py
new file mode 100644
index 0000000000..65c1529e69
--- /dev/null
+++ b/docs/examples/iv-modeling/plot_villalva.py
@@ -0,0 +1,444 @@
+"""
+Villalva single-diode model
+===========================
+
+This example extracts single-diode model (SDM) parameters from manufacturer
+specifications using Villalva's iterative fitting method, validates the fitted
+STC I-V curve, recovers a known synthetic SDM parameter set, and calculates
+I-V curves at several irradiance and cell-temperature conditions.
+
+The implementation follows Villalva's thesis, Chapter 3 and Appendix A, and
+Villalva, Gazoli, and Ruppert Filho (2009).
+"""
+
+# %%
+# Imports
+# -------
+# The fitting function is exposed through ``pvlib.ivtools.sdm`` and the
+# operating-condition auxiliary equations are exposed through
+# ``pvlib.pvsystem``.
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+from scipy import constants
+
+from pvlib import pvsystem
+from pvlib.ivtools.sdm import fit_villalva
+
+
+# %%
+# Example module
+# --------------
+# The JA Solar JAM78D40-625/MB datasheet values are used at STC. The module
+# has 156 half-cells, represented here as 78 effective series junctions.
+# Villalva's basic fitting method also requires a user-specified dimensionless
+# diode ideality factor.
+
+module = {
+    "Name": "JA Solar JAM78D40-625/MB",
+    "P_mp_ref": 625.0,
+    "V_oc_ref": 55.49,
+    "I_sc_ref": 14.36,
+    "V_mp_ref": 46.37,
+    "I_mp_ref": 13.48,
+    "cells_in_series": 78,
+    "alpha_sc_rel": 0.00046,
+    "beta_voc_rel": -0.00260,
+    "irrad_ref": 1000.0,
+    "temp_ref": 25.0,
+}
+
+module["alpha_sc"] = module["I_sc_ref"] * module["alpha_sc_rel"]
+module["beta_voc"] = module["V_oc_ref"] * module["beta_voc_rel"]
+module["P_mp_fit"] = module["V_mp_ref"] * module["I_mp_ref"]
+
+print(pd.Series(module))
+
+
+# %%
+# Fit the STC parameters
+# ----------------------
+# The algorithm increments ``R_s`` and selects the valid iteration that
+# minimizes the absolute difference between modeled and reference maximum
+# power, where the reference target is ``V_mp * I_mp``.
+
+diode_factor = 1.0
+
+params, history = fit_villalva(
+    v_mp=module["V_mp_ref"],
+    i_mp=module["I_mp_ref"],
+    v_oc=module["V_oc_ref"],
+    i_sc=module["I_sc_ref"],
+    alpha_sc=module["alpha_sc"],
+    beta_voc=module["beta_voc"],
+    cells_in_series=module["cells_in_series"],
+    diode_factor=diode_factor,
+    temp_ref=module["temp_ref"],
+    irrad_ref=module["irrad_ref"],
+    rs_step=1e-4,
+    rs_max=0.5,
+)
+
+parameter_table = pd.DataFrame(
+    {
+        "Value": [
+            diode_factor,
+            params["a_ref"],
+            params["I_L_ref"],
+            params["I_o_ref"],
+            params["R_s"],
+            params["R_sh_ref"],
+        ],
+        "Unit": ["-", "V", "A", "A", "ohm", "ohm"],
+    },
+    index=[
+        "diode_factor",
+        "a_ref",
+        "I_L_ref",
+        "I_o_ref",
+        "R_s",
+        "R_sh_ref",
+    ],
+)
+
+print(parameter_table)
+
+
+# %%
+# Evolution of maximum power with series resistance
+# -------------------------------------------------
+# This plot shows the complete valid fitting history and the selected value of
+# ``R_s``.
+
+best_idx = history["abs_power_error"].idxmin()
+best = history.loc[best_idx]
+
+print(f"Best R_s = {best['R_s']:.6f} ohm")
+print(f"Modeled P_mp = {best['p_mp_model']:.6f} W")
+print(f"Reference P_mp = {best['p_mp_ref']:.6f} W")
+print(f"Absolute error = {best['abs_power_error']:.6e} W")
+
+plt.figure(figsize=(9, 5))
+plt.plot(history["R_s"], history["p_mp_model"], label="Modeled $P_{mp}$")
+plt.axhline(
+    module["P_mp_fit"],
+    linestyle="--",
+    label="Reference $V_{mp}I_{mp}$",
+)
+plt.scatter(
+    [best["R_s"]],
+    [best["p_mp_model"]],
+    zorder=3,
+    label=f"Best $R_s$ = {best['R_s']:.4f} $\\Omega$",
+)
+plt.axvline(best["R_s"], linestyle=":")
+plt.xlim(0, 0.26)
+plt.xlabel("$R_s$ [$\\Omega$]")
+plt.ylabel("$P_{mp}$ [W]")
+plt.title("Villalva fitting evolution: $P_{mp}$ vs. $R_s$")
+plt.grid(True)
+plt.legend()
+plt.tight_layout()
+plt.show()
+
+
+# %%
+# STC I-V curve and datasheet validation
+# --------------------------------------
+# The fitted parameters are passed directly to pvlib's existing single-diode
+# equation solvers. The key modeled points are compared with the datasheet
+# values used by the fitting procedure.
+
+stc_sde = {
+    "photocurrent": params["I_L_ref"],
+    "saturation_current": params["I_o_ref"],
+    "resistance_series": params["R_s"],
+    "resistance_shunt": params["R_sh_ref"],
+    "nNsVth": params["a_ref"],
+}
+
+stc = pvsystem.singlediode(method="lambertw", **stc_sde)
+
+datasheet = pd.Series(
+    {
+        "i_sc": module["I_sc_ref"],
+        "v_oc": module["V_oc_ref"],
+        "i_mp": module["I_mp_ref"],
+        "v_mp": module["V_mp_ref"],
+        "p_mp": module["P_mp_fit"],
+    }
+)
+
+model = pd.Series(
+    {
+        key: float(np.asarray(stc[key]))
+        for key in ["i_sc", "v_oc", "i_mp", "v_mp", "p_mp"]
+    }
+)
+
+validation = pd.DataFrame({"Datasheet": datasheet, "Villalva": model})
+validation["Error"] = validation["Villalva"] - validation["Datasheet"]
+validation["Error [%]"] = (
+    100 * validation["Error"] / validation["Datasheet"]
+)
+
+print(validation)
+
+v_stc = np.linspace(0.0, float(stc["v_oc"]), 200)
+i_stc = pvsystem.i_from_v(
+    voltage=v_stc,
+    method="lambertw",
+    **stc_sde,
+)
+
+plt.figure(figsize=(9, 6))
+plt.plot(v_stc, i_stc, label="Villalva model")
+plt.scatter(
+    [0.0, module["V_mp_ref"], module["V_oc_ref"]],
+    [module["I_sc_ref"], module["I_mp_ref"], 0.0],
+    label="Datasheet key points",
+    zorder=3,
+)
+plt.xlabel("Module voltage [V]")
+plt.ylabel("Module current [A]")
+plt.title(module["Name"] + "\nSTC I-V curve")
+plt.xlim(left=0)
+plt.ylim(bottom=0)
+plt.grid(True)
+plt.legend()
+plt.tight_layout()
+plt.show()
+
+
+# %%
+# Recover a known SDM parameter set
+# ---------------------------------
+# A second validation starts from a complete known SDM parameter set. Synthetic
+# STC key points are generated with ``pvsystem.singlediode`` and then supplied
+# to ``fit_villalva``. This checks how closely the fitting method recovers the
+# original parameters when its inputs are internally consistent with the SDE.
+
+known = {
+    "Name": "Canadian Solar CS5P-220M",
+    "N_s": 96,
+    "alpha_sc": 0.004539,
+    "beta_voc": -0.22216,
+    "a_ref": 2.6373,
+    "I_L_ref": 5.114,
+    "I_o_ref": 8.196e-10,
+    "R_s": 1.065,
+    "R_sh_ref": 381.68,
+    "temp_ref": 25.0,
+    "irrad_ref": 1000.0,
+}
+
+T_ref_K = known["temp_ref"] + 273.15
+
+diode_factor_known = known["a_ref"] / (
+    known["N_s"] * constants.k * T_ref_K / constants.e
+)
+
+synthetic = pvsystem.singlediode(
+    photocurrent=known["I_L_ref"],
+    saturation_current=known["I_o_ref"],
+    resistance_series=known["R_s"],
+    resistance_shunt=known["R_sh_ref"],
+    nNsVth=known["a_ref"],
+    method="lambertw",
+)
+
+synthetic_points = {
+    key: float(np.asarray(synthetic[key]))
+    for key in ["i_sc", "v_oc", "i_mp", "v_mp", "p_mp"]
+}
+
+print("Synthetic STC key points")
+print(pd.Series(synthetic_points))
+
+recovered, recovery_history = fit_villalva(
+    v_mp=synthetic_points["v_mp"],
+    i_mp=synthetic_points["i_mp"],
+    v_oc=synthetic_points["v_oc"],
+    i_sc=synthetic_points["i_sc"],
+    alpha_sc=known["alpha_sc"],
+    beta_voc=known["beta_voc"],
+    cells_in_series=known["N_s"],
+    diode_factor=diode_factor_known,
+    temp_ref=known["temp_ref"],
+    irrad_ref=known["irrad_ref"],
+    rs_step=1e-4,
+    rs_max=1.5,
+)
+
+parameter_names = ["I_L_ref", "I_o_ref", "R_s", "R_sh_ref", "a_ref"]
+recovery_table = pd.DataFrame(
+    {
+        "Starting value": [known[name] for name in parameter_names],
+        "Recovered value": [recovered[name] for name in parameter_names],
+    },
+    index=parameter_names,
+)
+recovery_table["Absolute error"] = (
+    recovery_table["Recovered value"] - recovery_table["Starting value"]
+)
+recovery_table["Relative error [%]"] = (
+    100
+    * recovery_table["Absolute error"]
+    / recovery_table["Starting value"]
+)
+
+print(f"Known Villalva diode ideality factor: {diode_factor_known:.9f}")
+print(recovery_table)
+
+v_validation = np.linspace(0.0, synthetic_points["v_oc"], 300)
+
+i_original = pvsystem.i_from_v(
+    voltage=v_validation,
+    photocurrent=known["I_L_ref"],
+    saturation_current=known["I_o_ref"],
+    resistance_series=known["R_s"],
+    resistance_shunt=known["R_sh_ref"],
+    nNsVth=known["a_ref"],
+    method="lambertw",
+)
+
+i_recovered = pvsystem.i_from_v(
+    voltage=v_validation,
+    photocurrent=recovered["I_L_ref"],
+    saturation_current=recovered["I_o_ref"],
+    resistance_series=recovered["R_s"],
+    resistance_shunt=recovered["R_sh_ref"],
+    nNsVth=recovered["a_ref"],
+    method="lambertw",
+)
+
+plt.figure(figsize=(9, 6))
+plt.plot(v_validation, i_original, label="Starting SDM")
+plt.plot(
+    v_validation,
+    i_recovered,
+    linestyle="--",
+    label="Recovered Villalva SDM",
+)
+plt.scatter(
+    [0.0, synthetic_points["v_mp"], synthetic_points["v_oc"]],
+    [synthetic_points["i_sc"], synthetic_points["i_mp"], 0.0],
+    label="Synthetic key points",
+    zorder=3,
+)
+plt.xlabel("Module voltage [V]")
+plt.ylabel("Module current [A]")
+plt.title("Villalva parameter-recovery validation")
+plt.xlim(left=0)
+plt.ylim(bottom=0)
+plt.grid(True)
+plt.legend()
+plt.tight_layout()
+plt.show()
+
+
+# %%
+# I-V curves under operating conditions
+# -------------------------------------
+# The fitted reference parameters are converted to operating-condition SDM
+# parameters with ``pvsystem.calcparams_villalva``. pvlib's existing
+# ``singlediode`` and ``i_from_v`` functions then solve the corresponding
+# I-V curves.
+
+cases = [
+    (1000, 55),
+    (800, 55),
+    (600, 55),
+    (400, 25),
+    (400, 40),
+    (400, 55),
+]
+
+conditions = pd.DataFrame(cases, columns=["Geff", "Tcell"])
+
+IL, I0, Rs, Rsh, nNsVth = pvsystem.calcparams_villalva(
+    effective_irradiance=conditions["Geff"],
+    temp_cell=conditions["Tcell"],
+    alpha_sc=params["alpha_sc"],
+    beta_voc=params["beta_voc"],
+    a_ref=params["a_ref"],
+    I_L_ref=params["I_L_ref"],
+    R_sh_ref=params["R_sh_ref"],
+    R_s=params["R_s"],
+    i_sc_ref=params["i_sc_ref"],
+    v_oc_ref=params["v_oc_ref"],
+    irrad_ref=params["irrad_ref"],
+    temp_ref=params["temp_ref"],
+)
+
+SDE_params = {
+    "photocurrent": IL,
+    "saturation_current": I0,
+    "resistance_series": Rs,
+    "resistance_shunt": Rsh,
+    "nNsVth": nNsVth,
+}
+
+curve_info = pvsystem.singlediode(method="lambertw", **SDE_params)
+
+v = pd.DataFrame(np.linspace(0.0, curve_info["v_oc"], 100))
+i = pd.DataFrame(
+    pvsystem.i_from_v(
+        voltage=v,
+        method="lambertw",
+        **SDE_params,
+    )
+)
+
+plt.figure(figsize=(9, 6))
+
+for idx, case in conditions.iterrows():
+    label = (
+        "$G_{eff}$ "
+        + f"{case['Geff']:.0f} W/m$^2$\n"
+        "$T_{cell}$ "
+        + f"{case['Tcell']:.0f} °C"
+    )
+    plt.plot(v[idx], i[idx], label=label)
+    plt.plot(
+        [curve_info["v_mp"][idx]],
+        [curve_info["i_mp"][idx]],
+        linestyle="",
+        marker="o",
+    )
+
+plt.xlim(left=0)
+plt.ylim(bottom=0)
+plt.xlabel("Module voltage [V]")
+plt.ylabel("Module current [A]")
+plt.title(module["Name"])
+plt.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))
+plt.grid(True)
+plt.tight_layout()
+plt.show()
+
+operating_results = pd.DataFrame(
+    {
+        "Geff [W/m2]": conditions["Geff"],
+        "Tcell [degC]": conditions["Tcell"],
+        "i_sc [A]": curve_info["i_sc"],
+        "v_oc [V]": curve_info["v_oc"],
+        "i_mp [A]": curve_info["i_mp"],
+        "v_mp [V]": curve_info["v_mp"],
+        "p_mp [W]": curve_info["p_mp"],
+    }
+)
+
+print(operating_results)
+
+
+# %%
+# API placement
+# -------------
+# The proposed public functions are ``pvlib.ivtools.sdm.fit_villalva`` for
+# reference-condition parameter extraction and
+# ``pvlib.pvsystem.calcparams_villalva`` for the operating-condition
+# auxiliary equations. The actual I-V solution remains in pvlib's existing
+# ``singlediode`` and ``i_from_v`` functions.
diff --git a/pvlib/ivtools/sdm/__init__.py b/pvlib/ivtools/sdm/__init__.py
index 2bb8b9876b..d5bcf3f58a 100644
--- a/pvlib/ivtools/sdm/__init__.py
+++ b/pvlib/ivtools/sdm/__init__.py
@@ -19,3 +19,8 @@
     fit_pvsyst_iec61853_sandia_2025,
     pvsyst_temperature_coeff,
 )
+
+from pvlib.ivtools.sdm.villalva import ( # noqa: F401
+    _villalva_params_at_rs,
+    fit_villalva,
+)
diff --git a/pvlib/ivtools/sdm/villalva.py b/pvlib/ivtools/sdm/villalva.py
new file mode 100644
index 0000000000..f2539b5804
--- /dev/null
+++ b/pvlib/ivtools/sdm/villalva.py
@@ -0,0 +1,252 @@
+"""Functions for fitting the Villalva single-diode model."""
+
+import numpy as np
+import pandas as pd
+
+from scipy import constants
+
+from pvlib import pvsystem
+
+
+def _villalva_params_at_rs(
+    resistance_series,
+    v_mp,
+    i_mp,
+    v_oc,
+    i_sc,
+    a_ref,
+):
+    """Calculate ``I_L``, ``I_o`` and ``R_sh`` for one candidate ``R_s``."""
+    exp_voc = np.expm1(v_oc / a_ref)
+    exp_vmp = np.expm1((v_mp + i_mp * resistance_series) / a_ref)
+
+    p_mp = v_mp * i_mp
+
+    # Algebraic combination of Villalva's equations:
+    #
+    # I_L = (R_sh + R_s) / R_sh * I_sc
+    #
+    # I_o = (I_L - V_oc / R_sh) / (exp(V_oc / a_ref) - 1)
+    #
+    # and the R_sh(R_s) relation obtained at the MPP.
+    a_term = (
+        v_mp * i_sc
+        - v_mp * i_sc * exp_vmp / exp_voc
+        - p_mp
+    )
+
+    b_term = (
+        v_mp * i_sc * resistance_series
+        - v_mp
+        * (i_sc * resistance_series - v_oc)
+        * exp_vmp
+        / exp_voc
+    )
+
+    resistance_shunt = (
+        v_mp * (v_mp + i_mp * resistance_series) - b_term
+    ) / a_term
+
+    photocurrent = (
+        (resistance_shunt + resistance_series)
+        / resistance_shunt
+        * i_sc
+    )
+
+    saturation_current = (
+        photocurrent - v_oc / resistance_shunt
+    ) / exp_voc
+
+    return photocurrent, saturation_current, resistance_shunt
+
+
+def fit_villalva(
+    v_mp,
+    i_mp,
+    v_oc,
+    i_sc,
+    alpha_sc,
+    beta_voc,
+    cells_in_series,
+    diode_factor,
+    temp_ref=25.0,
+    irrad_ref=1000.0,
+    rs_step=1e-4,
+    rs_max=None,
+):
+    r"""Fit Villalva single-diode model parameters from datasheet values.
+
+    Villalva's basic parameter-extraction method increments the series
+    resistance from zero. For each candidate :math:`R_s`, the corresponding
+    :math:`R_{sh}`, :math:`I_L`, and :math:`I_0` are calculated and the
+    maximum power of the resulting single-diode model is evaluated. The
+    selected solution minimizes the absolute difference between modeled and
+    reference maximum power.
+
+    The dimensionless diode ideality factor :math:`n` is supplied by the user. The
+    modified ideality factor at reference conditions is calculated as
+
+    .. math::
+
+       a_{ref} = n N_s k T_{ref} / q.
+
+    Parameters
+    ----------
+    v_mp : float
+        Maximum-power voltage at reference conditions. [V]
+    i_mp : float
+        Maximum-power current at reference conditions. [A]
+    v_oc : float
+        Open-circuit voltage at reference conditions. [V]
+    i_sc : float
+        Short-circuit current at reference conditions. [A]
+    alpha_sc : float
+        Short-circuit current temperature coefficient. [A/K]
+    beta_voc : float
+        Open-circuit voltage temperature coefficient. [V/K]
+    cells_in_series : int
+        Effective number of cells or junctions connected in series.
+    diode_factor : float
+        Dimensionless diode ideality factor used by the Villalva model.
+    temp_ref : float, default 25
+        Reference cell temperature. [°C]
+    irrad_ref : float, default 1000
+        Reference irradiance. [W/m²]
+    rs_step : float, default 1e-4
+        Increment used for the series-resistance sweep. [ohm]
+    rs_max : float, optional
+        Maximum series resistance considered. If ``None``, ``v_oc / i_sc``
+        is used. [ohm]
+
+    Returns
+    -------
+    params : dict
+        Fitted reference-condition Villalva parameters. The dictionary
+        contains ``I_L_ref``, ``I_o_ref``, ``R_s``, ``R_sh_ref``, ``a_ref``,
+        ``alpha_sc``, ``beta_voc``, ``i_sc_ref``, ``v_oc_ref``,
+        ``irrad_ref``, and ``temp_ref``.
+    history : pandas.DataFrame
+        Full fitting history for all valid values of ``R_s`` evaluated by the
+        algorithm.
+
+    Raises
+    ------
+    ValueError
+        If the minimum Villalva shunt resistance is not positive.
+    RuntimeError
+        If no valid Villalva solution is found.
+
+    Notes
+    -----
+    The returned reference parameters are intended for use with
+    :py:func:`pvlib.pvsystem.calcparams_villalva`. I-V curve points can then
+    be calculated with :py:func:`pvlib.pvsystem.singlediode` or
+    :py:func:`pvlib.pvsystem.i_from_v`.
+
+    References
+    ----------
+    .. [1] M. G. Villalva, "Three-phase electronic power converter for a grid-connected photovoltaic system," 
+       PhD Thesis, Unicamp, 2010. DOI: 10.47749/T/UNICAMP.2010.781324
+    .. [2] M. G. Villalva, J. R. Gazoli, and E. Ruppert Filho, 
+       "Comprehensive Approach to Modeling and Simulation of Photovoltaic Arrays," 
+       IEEE Transactions on Power Electronics, 2009. DOI: 10.1109/TPEL.2009.2013862
+    """
+    temp_ref_k = temp_ref + 273.15
+
+    # a_ref = a * Ns * kT/q
+    a_ref = (
+        diode_factor
+        * cells_in_series
+        * constants.k
+        * temp_ref_k
+        / constants.e
+    )
+
+    p_mp_ref = v_mp * i_mp
+
+    # Villalva minimum shunt resistance.
+    r_sh_min = v_mp / (i_sc - i_mp) - (v_oc - v_mp) / i_mp
+
+    if r_sh_min <= 0:
+        raise ValueError("Villalva R_sh_min must be positive.")
+
+    if rs_max is None:
+        rs_max = v_oc / i_sc
+
+    rows = []
+    rs_values = np.arange(0.0, rs_max + 0.5 * rs_step, rs_step)
+
+    for resistance_series in rs_values:
+        try:
+            photocurrent, saturation_current, resistance_shunt = (
+                _villalva_params_at_rs(
+                    resistance_series,
+                    v_mp,
+                    i_mp,
+                    v_oc,
+                    i_sc,
+                    a_ref,
+                )
+            )
+        except (FloatingPointError, ZeroDivisionError):
+            continue
+
+        if (
+            not np.isfinite(resistance_shunt)
+            or resistance_shunt < r_sh_min
+            or photocurrent <= 0
+            or saturation_current <= 0
+        ):
+            continue
+
+        try:
+            mpp = pvsystem.max_power_point(
+                photocurrent=photocurrent,
+                saturation_current=saturation_current,
+                resistance_series=resistance_series,
+                resistance_shunt=resistance_shunt,
+                nNsVth=a_ref,
+                method="brentq",
+            )
+        except (ValueError, RuntimeError):
+            continue
+
+        p_mp_model = float(np.asarray(mpp["p_mp"]))
+        power_error = p_mp_model - p_mp_ref
+
+        rows.append(
+            {
+                "R_s": resistance_series,
+                "R_sh_ref": resistance_shunt,
+                "I_L_ref": photocurrent,
+                "I_o_ref": saturation_current,
+                "v_mp_model": float(np.asarray(mpp["v_mp"])),
+                "i_mp_model": float(np.asarray(mpp["i_mp"])),
+                "p_mp_model": p_mp_model,
+                "p_mp_ref": p_mp_ref,
+                "power_error": power_error,
+                "abs_power_error": abs(power_error),
+            }
+        )
+
+    if not rows:
+        raise RuntimeError("No valid Villalva solution was found.")
+
+    history = pd.DataFrame(rows)
+    best = history.loc[history["abs_power_error"].idxmin()]
+
+    params = {
+        "I_L_ref": float(best["I_L_ref"]),
+        "I_o_ref": float(best["I_o_ref"]),
+        "R_s": float(best["R_s"]),
+        "R_sh_ref": float(best["R_sh_ref"]),
+        "a_ref": float(a_ref),
+        "alpha_sc": alpha_sc,
+        "beta_voc": beta_voc,
+        "i_sc_ref": i_sc,
+        "v_oc_ref": v_oc,
+        "irrad_ref": irrad_ref,
+        "temp_ref": temp_ref,
+    }
+
+    return params, history
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index 79d8a58d11..521d9f43f9 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1678,6 +1678,124 @@ def get_orientation(self, solar_zenith, solar_azimuth):
         )
         return tracking_data
 
+def calcparams_villalva(
+    effective_irradiance,
+    temp_cell,
+    alpha_sc,
+    beta_voc,
+    a_ref,
+    I_L_ref,
+    R_sh_ref,
+    R_s,
+    i_sc_ref,
+    v_oc_ref,
+    irrad_ref=1000.0,
+    temp_ref=25.0,
+):
+    r"""Calculate Villalva SDM parameters at operating conditions.
+
+    The Villalva auxiliary equations convert reference-condition parameters
+    into the five single-diode equation parameters required by
+    :py:func:`pvlib.pvsystem.singlediode`.
+
+    The photocurrent is modeled as
+
+    .. math::
+
+       I_L = (I_{L,ref} + \alpha_{sc}\Delta T)
+             \frac{G}{G_{ref}}.
+
+    The finite-:math:`R_{sh}` open-circuit relation is used for the
+    temperature-dependent saturation current so that the operating-condition
+    equations remain consistent with the fitted reference parameters. The
+    fitted :math:`R_s` and :math:`R_{sh}` are held constant with operating
+    condition, following the Villalva modeling assumption used here.
+
+    Parameters
+    ----------
+    effective_irradiance : numeric
+        Effective irradiance converted to photocurrent. [W/m²]
+    temp_cell : numeric
+        Cell temperature. [°C]
+    alpha_sc : float
+        Short-circuit current temperature coefficient. [A/K]
+    beta_voc : float
+        Open-circuit voltage temperature coefficient. [V/K]
+    a_ref : float
+        Product of diode ideality factor, cells in series, and cell thermal
+        voltage at the reference temperature. [V]
+    I_L_ref : float
+        Light-generated current at reference conditions. [A]
+    R_sh_ref : float
+        Shunt resistance at reference conditions. [ohm]
+    R_s : float
+        Series resistance. [ohm]
+    i_sc_ref : float
+        Short-circuit current at reference conditions. [A]
+    v_oc_ref : float
+        Open-circuit voltage at reference conditions. [V]
+    irrad_ref : float, default 1000
+        Reference irradiance. [W/m²]
+    temp_ref : float, default 25
+        Reference cell temperature. [°C]
+
+    Returns
+    -------
+    photocurrent : numeric
+        Light-generated current at the operating condition. [A]
+    saturation_current : numeric
+        Diode reverse saturation current at the operating condition. [A]
+    resistance_series : numeric
+        Series resistance at the operating condition. [ohm]
+    resistance_shunt : numeric
+        Shunt resistance at the operating condition. [ohm]
+    nNsVth : numeric
+        Product of diode ideality factor, cells in series, and cell thermal
+        voltage at the operating condition. [V]
+
+    References
+    ----------
+    .. [1] M. G. Villalva, PhD thesis, Chapter 3 and Appendix A.
+    .. [2] M. G. Villalva, J. R. Gazoli, and E. Ruppert Filho,
+       "Comprehensive Approach to Modeling and Simulation of Photovoltaic
+       Arrays," IEEE Transactions on Power Electronics, 2009.
+    """
+    temp_ref_k = temp_ref + 273.15
+    temp_cell_k = temp_cell + 273.15
+    delta_t = temp_cell - temp_ref
+
+    # n * Ns * Vth at operating temperature.
+    nNsVth = a_ref * temp_cell_k / temp_ref_k
+
+    # Villalva photocurrent equation.
+    photocurrent = (
+        I_L_ref + alpha_sc * delta_t
+    ) * effective_irradiance / irrad_ref
+
+    # Datasheet temperature coefficients.
+    i_sc = i_sc_ref + alpha_sc * delta_t
+    v_oc = v_oc_ref + beta_voc * delta_t
+
+    # Finite-Rsh form corresponding to Villalva's temperature-dependent
+    # saturation-current equation.
+    photocurrent_for_i0 = (
+        (R_sh_ref + R_s) / R_sh_ref * i_sc
+    )
+
+    saturation_current = (
+        photocurrent_for_i0 - v_oc / R_sh_ref
+    ) / np.expm1(v_oc / nNsVth)
+
+    resistance_series = R_s
+    resistance_shunt = R_sh_ref
+
+    return (
+        photocurrent,
+        saturation_current,
+        resistance_series,
+        resistance_shunt,
+        nNsVth,
+    )
 
 def calcparams_desoto(effective_irradiance, temp_cell,
                       alpha_sc, a_ref, I_L_ref, I_o_ref, R_sh_ref, R_s,