diff --git a/pvlib/ivtools/sdm/_fit_desoto_pvsyst_sandia.py b/pvlib/ivtools/sdm/_fit_desoto_pvsyst_sandia.py
index 90240c9144..e22f3be1ca 100644
--- a/pvlib/ivtools/sdm/_fit_desoto_pvsyst_sandia.py
+++ b/pvlib/ivtools/sdm/_fit_desoto_pvsyst_sandia.py
@@ -9,6 +9,7 @@
 
 from pvlib.pvsystem import singlediode, v_from_i
 from pvlib.ivtools.utils import rectify_iv_curve, _numdiff
+from pvlib.pvsystem import _pvsyst_Rsh
 
 
 def _initial_iv_params(ivcurves, ee, voc, isc, rsh, nnsvth):
@@ -167,7 +168,10 @@ def _extract_sdm_params(ee, tc, iph, io, rs, rsh, n, u, specs, const,
         # Find parameters for Rsh equation
 
         def fun_rsh(x, rshexp, ee, e0, rsh):
-            tf = np.log10(_rsh_pvsyst(x, R_sh_exp, ee, e0)) - np.log10(rsh)
+            tf = (
+                np.log10(_pvsyst_Rsh(ee, x[1], x[0], R_sh_exp, e0))
+                - np.log10(rsh)
+            )
             return tf
 
         x0 = np.array([grsh0, grshref])
@@ -293,20 +297,6 @@ def _update_io(voc, iph, io, rs, rsh, nnsvth):
     return new_io
 
 
-def _rsh_pvsyst(x, rshexp, g, go):
-    # computes rsh for PVsyst model where the parameters are in vector xL
-    # x[0] = Rsh0
-    # x[1] = Rshref
-
-    rsho = x[0]
-    rshref = x[1]
-
-    rshb = np.maximum(
-        (rshref - rsho * np.exp(-rshexp)) / (1. - np.exp(-rshexp)), 0.)
-    rsh = rshb + (rsho - rshb) * np.exp(-rshexp * g / go)
-    return rsh
-
-
 def _filter_params(ee, isc, io, rs, rsh):
     # Function _filter_params identifies bad parameter sets. A bad set contains
     # Nan, non-positive or imaginary values for parameters; Rs > Rsh; or data
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index b762372ae7..896791d86f 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -1948,31 +1948,17 @@ def calcparams_pvsyst(effective_irradiance, temp_cell,
 
     '''
 
-    # Boltzmann constant in J/K
-    k = constants.k
+    gamma = _pvsyst_gamma(temp_cell, gamma_ref, mu_gamma, temp_ref)
 
-    # elementary charge in coulomb
-    q = constants.e
+    nNsVth = _pvsyst_nNsVth(temp_cell, gamma, cells_in_series)
 
-    # reference temperature
-    Tref_K = temp_ref + 273.15
-    Tcell_K = temp_cell + 273.15
-
-    gamma = gamma_ref + mu_gamma * (Tcell_K - Tref_K)
-    nNsVth = gamma * k / q * cells_in_series * Tcell_K
-
-    IL = effective_irradiance / irrad_ref * \
-        (I_L_ref + alpha_sc * (Tcell_K - Tref_K))
+    IL = _pvsyst_IL(effective_irradiance, temp_cell, I_L_ref, alpha_sc,
+                    irrad_ref, temp_ref)
 
-    I0 = I_o_ref * ((Tcell_K / Tref_K) ** 3) * \
-        (np.exp((q * EgRef) / (k * gamma) * (1 / Tref_K - 1 / Tcell_K)))
+    I0 = _pvsyst_Io(temp_cell, gamma, I_o_ref, EgRef, temp_ref)
 
-    Rsh_tmp = \
-        (R_sh_ref - R_sh_0 * np.exp(-R_sh_exp)) / (1.0 - np.exp(-R_sh_exp))
-    Rsh_base = np.maximum(0.0, Rsh_tmp)
-
-    Rsh = Rsh_base + (R_sh_0 - Rsh_base) * \
-        np.exp(-R_sh_exp * effective_irradiance / irrad_ref)
+    Rsh = _pvsyst_Rsh(effective_irradiance, R_sh_ref, R_sh_0, R_sh_exp,
+                      irrad_ref)
 
     Rs = R_s
 
@@ -1992,6 +1978,54 @@ def calcparams_pvsyst(effective_irradiance, temp_cell,
     return tuple(pd.Series(a, index=index).rename(None) for a in out)
 
 
+def _pvsyst_Rsh(effective_irradiance, R_sh_ref, R_sh_0, R_sh_exp=5.5,
+                irrad_ref=1000):
+    Rsh_tmp = \
+        (R_sh_ref - R_sh_0 * np.exp(-R_sh_exp)) / (1.0 - np.exp(-R_sh_exp))
+    Rsh_base = np.maximum(0.0, Rsh_tmp)
+
+    Rsh = Rsh_base + (R_sh_0 - Rsh_base) * \
+        np.exp(-R_sh_exp * effective_irradiance / irrad_ref)
+
+    return Rsh
+
+
+def _pvsyst_IL(effective_irradiance, temp_cell, I_L_ref, alpha_sc,
+               irrad_ref=1000, temp_ref=25):
+    Tref_K = temp_ref + 273.15
+    Tcell_K = temp_cell + 273.15
+    IL = effective_irradiance / irrad_ref * \
+        (I_L_ref + alpha_sc * (Tcell_K - Tref_K))
+    return IL
+
+
+def _pvsyst_Io(temp_cell, gamma, I_o_ref, EgRef, temp_ref=25):
+    k = constants.k  # Boltzmann constant in J/K
+    q = constants.e  # elementary charge in coulomb
+
+    Tref_K = temp_ref + 273.15
+    Tcell_K = temp_cell + 273.15
+
+    Io = I_o_ref * ((Tcell_K / Tref_K) ** 3) * \
+        (np.exp((q * EgRef) / (k * gamma) * (1 / Tref_K - 1 / Tcell_K)))
+
+    return Io
+
+
+def _pvsyst_gamma(temp_cell, gamma_ref, mu_gamma, temp_ref=25):
+    gamma = gamma_ref + mu_gamma * (temp_cell - temp_ref)
+    return gamma
+
+
+def _pvsyst_nNsVth(temp_cell, gamma, cells_in_series):
+    k = constants.k  # Boltzmann constant in J/K
+    q = constants.e  # elementary charge in coulomb
+    Tcell_K = temp_cell + 273.15
+
+    nNsVth = gamma * k / q * cells_in_series * Tcell_K
+    return nNsVth
+
+
 def retrieve_sam(name=None, path=None):
     """
     Retrieve latest module and inverter info from a file bundled with pvlib,