55from diffpy .utils .parsers .loaddata import loadData
66
77
8- def _top_hat (x , slit_width ):
8+ def _top_hat (x , half_slit_width ):
99 """
1010 create a top-hat function, return 1.0 for values within the specified slit width and 0 otherwise
1111 """
12- return np .where ((x >= - slit_width ) & (x <= slit_width ), 1.0 , 0 )
12+ return np .where ((x >= - half_slit_width ) & (x <= half_slit_width ), 1.0 , 0 )
1313
1414
1515def _model_function (x , diameter , x0 , I0 , mud , slope ):
@@ -32,7 +32,7 @@ def _model_function(x, diameter, x0, I0, mud, slope):
3232 return (I0 - slope * x ) * np .exp (- mud / diameter * length )
3333
3434
35- def _extend_x_and_convolve (x , diameter , slit_width , x0 , I0 , mud , slope ):
35+ def _extend_x_and_convolve (x , diameter , half_slit_width , x0 , I0 , mud , slope ):
3636 """
3737 extend x values and I values for padding (so that we don't have tails in convolution), then perform convolution
3838 (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):
4242 x_right_pad = np .linspace (x .max (), x .max () + n_points * (x [1 ] - x [0 ]), n_points )
4343 x_extended = np .concatenate ([x_left_pad , x , x_right_pad ])
4444 I_extended = _model_function (x_extended , diameter , x0 , I0 , mud , slope )
45- kernel = _top_hat (x_extended - x_extended .mean (), slit_width )
45+ kernel = _top_hat (x_extended - x_extended .mean (), half_slit_width )
4646 I_convolved = I_extended # this takes care of the case where slit width is close to 0
4747 if kernel .sum () != 0 :
4848 kernel /= kernel .sum ()
@@ -56,8 +56,8 @@ def _objective_function(params, x, observed_data):
5656 compute the objective function for fitting a model to the observed/experimental data
5757 by minimizing the sum of squared residuals between the observed data and the convolved model data
5858 """
59- diameter , slit_width , x0 , I0 , mud , slope = params
60- convolved_model_data = _extend_x_and_convolve (x , diameter , slit_width , x0 , I0 , mud , slope )
59+ diameter , half_slit_width , x0 , I0 , mud , slope = params
60+ convolved_model_data = _extend_x_and_convolve (x , diameter , half_slit_width , x0 , I0 , mud , slope )
6161 residuals = observed_data - convolved_model_data
6262 return np .sum (residuals ** 2 )
6363
@@ -68,33 +68,37 @@ def _compute_single_mud(x_data, I_data):
6868 """
6969 bounds = [
7070 (1e-5 , x_data .max () - x_data .min ()), # diameter: [small positive value, upper bound]
71- (0 , (x_data .max () - x_data .min ()) / 2 ), # slit width: [0, upper bound]
71+ (0 , (x_data .max () - x_data .min ()) / 2 ), # half slit width: [0, upper bound]
7272 (x_data .min (), x_data .max ()), # x0: [min x, max x]
7373 (1e-5 , I_data .max ()), # I0: [small positive value, max observed intensity]
7474 (1e-5 , 20 ), # muD: [small positive value, upper bound]
75- (- 10000 , 10000 ), # slope: [lower bound, upper bound]
75+ (- 100000 , 100000 ), # slope: [lower bound, upper bound]
7676 ]
7777 result = dual_annealing (_objective_function , bounds , args = (x_data , I_data ))
78- diameter , slit_width , x0 , I0 , mud , slope = result .x
79- convolved_fitted_signal = _extend_x_and_convolve (x_data , diameter , slit_width , x0 , I0 , mud , slope )
78+ diameter , half_slit_width , x0 , I0 , mud , slope = result .x
79+ convolved_fitted_signal = _extend_x_and_convolve (x_data , diameter , half_slit_width , x0 , I0 , mud , slope )
8080 residuals = I_data - convolved_fitted_signal
8181 rmse = np .sqrt (np .mean (residuals ** 2 ))
8282 return mud , rmse
8383
8484
8585def compute_mud (filepath ):
8686 """
87- compute the best-fit mu*D value from a z-scan file
87+ compute the best-fit mu*D value from a z-scan file, removing the sample holder effect
88+
89+ This function loads z-scan data and fits it to a model
90+ that convolves a top-hat function with I = I0 * e^{-mu * l}.
91+ The fitting procedure is run multiple times, and we return the best-fit parameters based on the lowest rmse.
8892
8993 Parameters
9094 ----------
91- filepath str
92- the path to the z-scan file
95+ filepath: str
96+ The path to the z-scan file
9397
9498 Returns
9599 -------
96100 a float contains the best-fit mu*D value
97101 """
98102 x_data , I_data = loadData (filepath , unpack = True )
99- best_mud , _ = min ((_compute_single_mud (x_data , I_data ) for _ in range (10 )), key = lambda pair : pair [1 ])
103+ best_mud , _ = min ((_compute_single_mud (x_data , I_data ) for _ in range (20 )), key = lambda pair : pair [1 ])
100104 return best_mud
0 commit comments