I’m relatively new to Python GUI and have been working on a simple project relating to my physics classes.
I’ve been trying to simulate the Rutherford Scattering experimental data using Monte Carlo method, as well as using ipywidgets for an interactive simulation where you can use the sliders to change the values in the equation.
The issue that is occurring is that when I try and do a curve fit of the simulated data, changing the parameters of N_i values/Kinetic energy/Thickness/Element of Foil, the curve fit sometimes works but when you change the parameters to something different, it flatlines after 10 degrees (x_values). Then if you change the parameters back to how it was, it will sometimes fail or work again.
I searched online and it says that the issue may arise from the initial parameters being off but when I attempt to use different values or update the code any way, it still causes the curve fit to break. This is the closest I can get the code to match that of the Rutherford Scattering experiment.
This is the code that I’m working with. Every time the Run button is pressed, it plots a simulated data of the Rutherford Equation, which has been consistent for varying N_i values. Then you press the Analytical Curve button, which would give the theoretical curve of the equation. Then the Fit Curve button should give a curve that matches closely with the simulated data. There are sliders for N_i values, the kinetic energy, foil thickness and the element the foil is made of.
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import ipywidgets as widgets
from IPython.display import display, Latex
import time
# Constants
Z1 = 2 # Atomic number of alpha particle
e = 1.602e-19 # Elementary charge (C)
epsilon_0 = 8.854e-12 # Vacuum permittivity (F/m)
M_alpha = 6.644e-27 # Mass of alpha particle (kg)
# Simulation Function: Rutherford's scattering formula
def rutherford_scattering(theta, N_i, a, t, r, Z2, EK):
"""Compute N(theta) based on Rutherford's law."""
first = (N_i * a * t) / (16 * r**2)
second = ((Z1 * Z2 * e**2) / (4 * np.pi * epsilon_0 * EK))**2
return first * second * (1 / np.sin(theta / 2)**4)
# Run the simulation (Monte Carlo simulation)
def monte_carlo_simulation(analytical, theta, EK_slider_value, Ni_slider_value, noise_level=0.2):
"""Add random noise to the analytical data to simulate experimental results."""
# Apply more noise at larger angles (factor based on angle)
angle_factor = 1 / np.sin(theta / 2)**4 # Using angle factor for noise
# Scale noise with respect to kinetic energy (lower noise at higher KE)
# Adjust the scaling factor based on a more realistic energy value in Joules (not multiplied by 1e-13 here)
scaled_noise_level = 1 / np.sqrt(EK_slider_value) # Inversely proportional to KE
# Noise scaling based on the number of particles (Ni)
noise_scaling = 1 / np.sqrt(Ni_slider_value) # More noise for smaller Ni
# Combine noise scaling
total_noise_level = scaled_noise_level + noise_scaling + angle_factor
# Generate noise (normally distributed)
noise = np.random.normal(0, total_noise_level, size=theta.shape) # Generate noise for each data point
# Add noise to the analytical values
simulated_data = analytical + noise
# Ensure that N(theta) does not go below zero
simulated_data = np.maximum(simulated_data, 0) # Clamp values to be >= 0
return simulated_data
# Widgets for input parameters
Ni_slider = widgets.IntSlider(value=13000, min=6000, max=20000, step=500, description="N<sub>i</sub>: No. alpha particles", style={'description_width': 'initial'}, layout={'width': '400px'})
t_input = widgets.FloatText(value=400, min=200, max=1400, step=100, description="t: thickness of foil (nm)", style={'description_width': 'initial'})
EK_slider = widgets.FloatSlider(value=5, min=1, max=10, step=1, description="E<sub>K</sub>: Kinetic energy (J)", style={'description_width': 'initial'}, layout={'width': '400px'})
Z2_dropdown = widgets.Dropdown(
options={
'Gold (79)': 79,
'Silver (47)': 47,
'Copper (29)': 29,
'Aluminum (13)': 13,
'Iron (26)': 26,
'Titanium (22)': 22,
'Carbon (6)': 6,
'Oxygen (8)': 8,
'Lead (82)': 82,
'Mercury (80)': 80
},
value=79,
description="Z<sub>2</sub>: Element of foil",
style={'description_width': 'initial'}
)
# Buttons for running simulation, curve fitting, and displaying results
simulate_button = widgets.Button(description="Run Monte Carlo Simulation")
add_analytical_button = widgets.Button(description="Add Analytical Curve")
fit_button = widgets.Button(description="Fit Curve")
# Output plots and error messages
output_plot = widgets.Output()
output_error = widgets.Output()
output_analytical = widgets.Output()
# Layout for widgets
ui = widgets.VBox([
widgets.HBox([Ni_slider]),
widgets.HBox([EK_slider]),
widgets.VBox([t_input]),
widgets.VBox([Z2_dropdown]),
widgets.HBox([simulate_button, add_analytical_button, fit_button]),
output_plot, # Main plot for simulation
output_analytical, # Dedicated plot for analytical curve
output_error
])
# Global variables for storing simulated data and analytical data
theta_values = np.linspace(np.radians(1), np.radians(179), 500)
simulated_data = None
analytical_data = None
# Debounce function to prevent rapid successive calls
def debounce(f, wait=1):
"""Debounce a function to delay execution until user stops adjusting sliders."""
last_call_time = [0] # Use mutable object to retain state between calls
def wrapped_function(change):
now = time.time()
if now - last_call_time[0] > wait:
last_call_time[0] = now
f(change)
return wrapped_function
@debounce # Updated functions with debouncing applied
# Run the simulation (Monte Carlo simulation)
def run_simulation(change):
global simulated_data
with output_plot:
output_plot.clear_output() # Clear previous content in the output area
# Read values from widgets
Ni = Ni_slider.value
t_nm = t_input.value # Thickness in nanometers
t = t_nm * 1e-9 # Convert thickness to meters
EK = EK_slider.value * 1e-13 # Ensure the kinetic energy is properly converted to joules
Z2 = Z2_dropdown.value
r = 0.1 # Fixed distance from foil (m)
a = 1e28 # Approx. atoms/m^3 for dense materials
# Compute the theoretical Rutherford scattering distribution (analytical curve)
analytical = rutherford_scattering(theta_values, Ni, a, t, r, Z2, EK)
global analytical_data
analytical_data = analytical # Store analytical data for future use
# Simulate the data by adding noise to the theoretical model
simulated_data = monte_carlo_simulation(analytical_data, theta_values, EK_slider.value, Ni_slider.value, noise_level=0.1)
# Scatter plot: Simulated data
plt.figure(figsize=(8, 6))
plt.plot(np.degrees(theta_values), simulated_data, '.', color='green', label="Simulated Data")
plt.xlabel("Scattering Angle (degrees)")
plt.ylabel(r"$N(theta)$")
plt.yscale('log')
plt.title("Simulated Data")
plt.grid()
plt.tight_layout()
plt.legend()
plt.show()
@debounce
# Add analytical Rutherford scattering curve to the simulated data
def add_analytical_curve(change):
global simulated_data, analytical_data
if analytical_data is None or simulated_data is None:
with output_error:
output_error.clear_output()
print("Run the simulation first to generate simulated data and analytical curve.")
return # Ensure both simulated and analytical data exist
# Use the main output area to display both curves together
with output_plot:
output_plot.clear_output() # Clear previous content in the main plot area
# Create the plot
plt.figure(figsize=(8, 6))
plt.plot(np.degrees(theta_values), simulated_data, '.', color='green', label="Simulated Data")
plt.plot(np.degrees(theta_values), analytical_data, '-', color='blue', label="Analytical Curve")
plt.xlabel("Scattering Angle (degrees)")
plt.ylabel(r"$N(theta)$")
plt.yscale('log')
plt.title("Simulated Data with Analytical Curve")
plt.grid()
plt.tight_layout()
plt.legend()
plt.show()
# Curve fitting function
def fit_curve_function(theta, A, B):
"""Model function for fitting a curve to the data."""
return A / np.sin(theta / 2)**4 + B
@debounce
# Fit the simulated data
def fit_simulation(change):
global simulated_data, analytical_data
if analytical_data is None or simulated_data is None:
with output_error:
output_error.clear_output()
print("Run the simulation first to generate simulated data and analytical curve.")
return
try:
# Better initial parameter estimates
A_init = np.max(simulated_data) * (np.sin(np.pi/4))**4
B_init = np.min(simulated_data)
# Try multiple initial guesses if the first fit fails
initial_guesses = [
[A_init, B_init],
[A_init * 10, B_init],
[A_init * 0.1, B_init],
[A_init, 0]
]
best_fit = None
best_residual = np.inf
for guess in initial_guesses:
try:
# Wider bounds
bounds = ([0, 0], [np.inf, np.max(simulated_data)])
# Add method='trf' and increase max iterations
popt_try, pcov_try = curve_fit(
fit_curve_function,
theta_values,
simulated_data,
p0=guess,
bounds=bounds,
method='trf',
maxfev=20000
)
# Calculate fit quality
fitted = fit_curve_function(theta_values, *popt_try)
residual = np.sum((simulated_data - fitted)**2)
if residual < best_residual:
best_residual = residual
best_fit = (popt_try, pcov_try)
except:
continue
if best_fit is None:
raise RuntimeError("Failed to find a good fit with any initial parameters")
popt, pcov = best_fit
# Rest of your plotting code remains the same
A, B = popt
fitted_curve = fit_curve_function(theta_values, *popt)
with output_plot:
output_plot.clear_output()
plt.figure(figsize=(8, 6))
plt.plot(np.degrees(theta_values), simulated_data, '.', color='green', label="Simulated Data")
plt.plot(np.degrees(theta_values), analytical_data, '-', color='blue', label="Analytical Curve")
plt.plot(np.degrees(theta_values), fitted_curve, '--', color='red', label="Fitted Curve")
plt.xlabel("Scattering Angle (degrees)")
plt.ylabel(r"$N(theta)$")
plt.yscale('log')
plt.title("Curve Fitting to Simulated Data")
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()
# Display fit parameters and covariance matrix
with output_error:
output_error.clear_output()
print(f"Fitted Parameters: A = {A:.3e}, B = {B:.3e}")
print(f"Covariance Matrix: n{pcov}")
except Exception as e:
with output_error:
print(f"Error during curve fitting: {e}")
# Connect buttons
simulate_button.on_click(run_simulation)
add_analytical_button.on_click(add_analytical_curve)
fit_button.on_click(fit_simulation)
# Display the widgets
display(ui)
This code results in the first two buttons working completely fine, the sliders showing relatively close plots to what it would be in real life. However, sometimes when Fit Curve is pressed, the curve fit flatlines past 10 degrees. It will work for some values but then not work when put on those values again.
Is there something wrong with the format or the layout that’s causing this issue? Or is a logarithmic issue since my y-values are log for the large values? I’ve tried as much as I can and I am unsure where to go further.
Thank you
edit: I’ve also noticed the issue with the noise levels for KE. Usually at high KE, there needs to be less noise as there is less scatter, and for low KE there needs to be more noise for more scatter. I’ve noticed it’s doing the exact opposite, yet I’ve done the inverse and it still doesn’t change it so I’m assuming that its accumulating the noise from somewhere else
Keanna is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
Or is a logarithmic issue since my y-values are log for the large values?
Sort of, yeah. When you fit a function like this that spans from 10^7 to 10-1, curve_fit() cares about a mistake of 1 on a value of y=10^7 just as much as it cares about a mistake of 1 on a value of y=1. In other words, if it can improve the fit for the values at y=10^7 by 0.001%, this is worth making the curve fit at y=10^-1 terrible.
One way of fixing this is to switch from an absolute error to a relative error, by fitting the logarithm of your function to the logarithm of your data.
Example:
log_offset = 0.001
popt_try, pcov_try = curve_fit(
lambda *args: np.log(fit_curve_function(*args) + log_offset),
theta_values,
np.log(simulated_data + log_offset),
p0=guess,
bounds=bounds,
method='trf',
maxfev=20000
)
Note: I am actually applying the function np.log(x + log_offset), because your data contains zeros. The log_offset is a parameter which you should probably tune for your application. It controls how large of a value an input of zero has when transformed into log-space. Alternatively, you can drop zeros from your data.
Here is what the fit looks like without this change:
Before change:
After change:
Once we log-transform, we get what looks visually better fit, but if we compare the residual before and after, we can see that the new fit is actually worse, in an absolute sense.
absolute error residual SSE 1.1637203281333856e+16
relative error residual SSE 1.922495566518763e+16
1