diff --git a/deerlab/fitmodel.py b/deerlab/fitmodel.py
index 9ebff6ac1..1425ceaf2 100644
--- a/deerlab/fitmodel.py
+++ b/deerlab/fitmodel.py
@@ -568,7 +568,7 @@ def separable_nonlinear_lsq_analysis(Vexp):
         prescales = [1 for V in Vexp]
         Vexp_ = [Vexp[i]/prescales[i] for i in range(nSignals)]
 
-        def scale_constraint(nonlinpar):
+        def scale_constraint(nonlinpar,linpar):
         # --------------------------------------------------------
             penalty = np.zeros(nSignals)
             for i in range(nSignals):
diff --git a/deerlab/snlls.py b/deerlab/snlls.py
index 3628a1b2d..ab75dcca3 100644
--- a/deerlab/snlls.py
+++ b/deerlab/snlls.py
@@ -87,9 +87,10 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx
         The regularization parameter can be manually specified by passing a scalar value
         instead of a string. The default ``'aic'``.
 
-    custom_penalty: callable 
-        Custom penalty function to impose upon the solution. Must return a vector to be
-        added to the residual vector. 
+    extrapenalty: callable 
+        Custom penalty function to impose upon the solution. Must take two inputs, a vector of non-linear parameters
+        and a vector of linear parameters, and return a vector to be added to the residual vector (``pen = fcn(pnonlin,plin)``).  
+        The square of the penalty is computed internally.
 
     alphareopt : float scalar, optional
         Relative parameter change threshold for reoptimizing the regularization parameter
@@ -383,7 +384,7 @@ def ResidualsFcn(p):
 
         # Compute residual from custom penalty
         if callable(extrapenalty):
-            penres = extrapenalty(p)
+            penres = extrapenalty(p,linfit)
             penres = np.atleast_1d(penres)
             res = np.concatenate((res,penres))
 
@@ -434,10 +435,11 @@ def ResidualsFcn(p):
         Jnonlin = Jacobian(ResidualsFcn,nonlinfit,lb,ub)
 
         # Jacobian (linear part)
-        Jlin = np.zeros((len(res),len(linfit)))
-        Jlin[:len(y),:] = scales_vec[:,np.newaxis]*Amodel(nonlinfit)
+        Jlin = scales_vec[:,np.newaxis]*Amodel(nonlinfit)
+        if callable(extrapenalty):
+            Jlin = np.concatenate((Jlin, Jacobian(lambda plin: extrapenalty(nonlinfit,plin),linfit,lbl,ubl)))
         if includeRegularization:
-            Jlin[len(res)-Nlin:,:] = reg_penalty(regtype, alpha, L, linfit, huberparam, Nnonlin)[1]
+            Jlin = np.concatenate((Jlin, reg_penalty(regtype, alpha, L, linfit, huberparam, Nnonlin)[1]))
 
         # Full Jacobian
         J = np.concatenate((Jnonlin,Jlin),axis=1)