rename variables in mud_calculator.py and write more informative docstring

news/muD.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* no news added - minor edits in mud_calculator.py
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

src/diffpy/labpdfproc/mud_calculator.py

@@ -5,11 +5,11 @@
 from diffpy.utils.parsers.loaddata import loadData
 
 
-def _top_hat(x, slit_width):
+def _top_hat(x, half_slit_width):
     """
     create a top-hat function, return 1.0 for values within the specified slit width and 0 otherwise
     """
-    return np.where((x >= -slit_width) & (x <= slit_width), 1.0, 0)
+    return np.where((x >= -half_slit_width) & (x <= half_slit_width), 1.0, 0)
 
 
 def _model_function(x, diameter, x0, I0, mud, slope):
@@ -32,7 +32,7 @@ def _model_function(x, diameter, x0, I0, mud, slope):
     return (I0 - slope * x) * np.exp(-mud / diameter * length)
 
 
-def _extend_x_and_convolve(x, diameter, slit_width, x0, I0, mud, slope):
+def _extend_x_and_convolve(x, diameter, half_slit_width, x0, I0, mud, slope):
     """
     extend x values and I values for padding (so that we don't have tails in convolution), then perform convolution
     (note that the convolved I values are the same as modeled I values if slit width is close to 0)
@@ -42,7 +42,7 @@ def _extend_x_and_convolve(x, diameter, slit_width, x0, I0, mud, slope):
     x_right_pad = np.linspace(x.max(), x.max() + n_points * (x[1] - x[0]), n_points)
     x_extended = np.concatenate([x_left_pad, x, x_right_pad])
     I_extended = _model_function(x_extended, diameter, x0, I0, mud, slope)
-    kernel = _top_hat(x_extended - x_extended.mean(), slit_width)
+    kernel = _top_hat(x_extended - x_extended.mean(), half_slit_width)
     I_convolved = I_extended  # this takes care of the case where slit width is close to 0
     if kernel.sum() != 0:
         kernel /= kernel.sum()
@@ -56,8 +56,8 @@ def _objective_function(params, x, observed_data):
     compute the objective function for fitting a model to the observed/experimental data
     by minimizing the sum of squared residuals between the observed data and the convolved model data
     """
-    diameter, slit_width, x0, I0, mud, slope = params
-    convolved_model_data = _extend_x_and_convolve(x, diameter, slit_width, x0, I0, mud, slope)
+    diameter, half_slit_width, x0, I0, mud, slope = params
+    convolved_model_data = _extend_x_and_convolve(x, diameter, half_slit_width, x0, I0, mud, slope)
     residuals = observed_data - convolved_model_data
     return np.sum(residuals**2)
 
@@ -68,33 +68,37 @@ def _compute_single_mud(x_data, I_data):
     """
     bounds = [
         (1e-5, x_data.max() - x_data.min()),  # diameter: [small positive value, upper bound]
-        (0, (x_data.max() - x_data.min()) / 2),  # slit width: [0, upper bound]
+        (0, (x_data.max() - x_data.min()) / 2),  # half slit width: [0, upper bound]
         (x_data.min(), x_data.max()),  # x0: [min x, max x]
         (1e-5, I_data.max()),  # I0: [small positive value, max observed intensity]
         (1e-5, 20),  # muD: [small positive value, upper bound]
-        (-10000, 10000),  # slope: [lower bound, upper bound]
+        (-100000, 100000),  # slope: [lower bound, upper bound]
     ]
     result = dual_annealing(_objective_function, bounds, args=(x_data, I_data))
-    diameter, slit_width, x0, I0, mud, slope = result.x
-    convolved_fitted_signal = _extend_x_and_convolve(x_data, diameter, slit_width, x0, I0, mud, slope)
+    diameter, half_slit_width, x0, I0, mud, slope = result.x
+    convolved_fitted_signal = _extend_x_and_convolve(x_data, diameter, half_slit_width, x0, I0, mud, slope)
     residuals = I_data - convolved_fitted_signal
     rmse = np.sqrt(np.mean(residuals**2))
     return mud, rmse
 
 
 def compute_mud(filepath):
-    """
-    compute the best-fit mu*D value from a z-scan file
+    """Compute the best-fit mu*D value from a z-scan file, removing the sample holder effect.
+
+    This function loads z-scan data and fits it to a model
+    that convolves a top-hat function with I = I0 * e^{-mu * l}.
+    The fitting procedure is run multiple times, and we return the best-fit parameters based on the lowest rmse.
 
     Parameters
     ----------
-    filepath str
-        the path to the z-scan file
+    filepath : str
+        The path to the z-scan file.
 
     Returns
     -------
-    a float contains the best-fit mu*D value
+    mu*D : float
+        The best-fit mu*D value.
     """
     x_data, I_data = loadData(filepath, unpack=True)
-    best_mud, _ = min((_compute_single_mud(x_data, I_data) for _ in range(10)), key=lambda pair: pair[1])
+    best_mud, _ = min((_compute_single_mud(x_data, I_data) for _ in range(20)), key=lambda pair: pair[1])
     return best_mud