Taking data of binding energies of nucleii, a model fit of the Semi-Empirical Mass formula can be made.
Semi-Empirical Mass Formula:
Reading data from the file using numpy loadtxt, only the 1st, 2nd and 5th columns are of interest.
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display
import ipywidgets as widgets
from ipywidgets import interact, fixed
%matplotlib notebook
with open('data_SEMF.dat') as f:
lines = (line for line in f if not line.startswith('#'))
read_data = np.loadtxt(lines, delimiter = ' ', usecols=(0, 1, 5), unpack = True)
The curve_fit function has the following usage. If we have a model function, $f$, we assume the data, $y$, is simply $y = f(x) +$ errors then a curve fit can be performed in the following way:
Define the model function which takes $x$ and some free parameters and returns the model solution
Call the curve_fit function, passing it the model function, the x data, the y data and optionally a list of initial values for the fitting parameters
curve_fit will return 2 tuples: the first contains the least squared fit model parameters and the second contains the covariance matrix
Finally the model function can be called individually, passing it the fitted parameters, in order to obtain the fitted model
from scipy.optimize import curve_fit
def model_fit(data,term_mask,exclude):
A = data[0][exclude:]
Z = data[1][exclude:]
B = data[2][exclude:]
def _semi_empirical_mass(X, a_v, a_s, a_c, a_a, a_p):
model = np.zeros_like(X)
if term_mask[0]:
vol = a_v
model += vol
if term_mask[1]:
sa = -a_s*X**(-1./3.)
model += sa
if term_mask[2]:
es = -a_c*Z**2*X**(-4./3.)
model += es
if term_mask[3]:
asy = -a_a*((X-2*Z)**2)*X**(-2)
model += asy
if term_mask[4]:
N = X-Z
parity = np.zeros_like(X)
for i in range(len(parity)):
if N[i] % 2 == 0 and Z[i] % 2 == 0:
parity[i] = a_p*X[i]**(-3./4.)
elif N[i] % 2 == 1 and Z[i] % 2 == 1:
parity[i] = -a_p*X[i]**(-3./4.)
else:
parity[i] = 0.0
model += parity/X
return model
p0 = [15.8,18.0,0.72,23.5,33.5]
popt, pcov = curve_fit(_semi_empirical_mass, A, B, p0 = p0)
model = _semi_empirical_mass(A, popt[0], popt[1], popt[2], popt[3], popt[4])
ax.set_title((str(popt)))
line.set_data(A,model)
display(fig)
def widget_update(b_v = True,b_s = True,b_c = True,b_a = True,b_p = True, exclude=0):
term_mask = [b_v,b_s,b_c,b_a,b_p]
model_fit(read_data,term_mask,exclude)
Below is an interactive plot for the binding energy data and the fitted SEMF model. The fitted parameters are shown in the figure title and individual terms in the SEMF can be turned off to see their effect.
fig = plt.figure(figsize = (12,6))
ax = plt.subplot(111)
ax.set_xlabel('Atomic Number / A')
ax.set_ylabel('Binding Energy per Nucleon (MeV) / (B/A)')
ax.plot(read_data[0],read_data[2],'bx')
ax.set_xlim([0,max(read_data[0])])
line, = ax.plot([],[],'r-')
interact(widget_update, b_v = widgets.Checkbox(description='Volume Term',value=True),
b_s = widgets.Checkbox(description='Surface Area Term',value=True),
b_c = widgets.Checkbox(description='Coulomb Term',value=False),
b_a = widgets.Checkbox(description='Asymmetry Term',value=False),
b_p = widgets.Checkbox(description='Parity Term',value=False),
exclude = widgets.IntSlider(description='Exclude first n points',
min=0,max=50))
There is a data file called data_exercise.dat which contains a broad spectral line feature. It is proposed that the spectral line should either have a Lorentzian profile, due to lifetime broadening, or a Gaussian profile, due to thermal broadening. These have the following functional forms:
Lorentzian:
Write a program using curve_fit to fit a Gaussian and Lorentzian to the spectral feature and identify which model is a better fit. The frequency is given in THz and the emitting substance can assumed to be hydrogen.
For the chosen model, find the physical parameters e.g. decay time $\tau_s$ for lifetime broadening or temperature $T$ for thermal broadening
with open('data_exercise.dat') as f:
lines = (line for line in f if not line.startswith('#'))
read_data = np.loadtxt(lines, delimiter = ' ', unpack = True)
freq = read_data[0]
data = read_data[1]
fig3 = plt.figure(figsize = (12,6))
ax3 = plt.subplot(111)
ax3.set_xlim((np.min(freq),np.max(freq)))
ax3.plot(freq,data);
def gaussian(x,mu,sigma,amp):
p = amp*np.exp(-(x-mu)**2/(2*sigma**2))/np.sqrt(2*np.pi*sigma**2)
return p
def lorentzian(x,w_0,w_s,amp):
p = amp/((x-w_0)**2+w_s**2)
return p
p0 = [550.0,10.,1.]
popt_g, pcov = curve_fit(gaussian, freq, data, p0 = p0)
gauss_model = gaussian(freq,popt_g[0],popt_g[1],popt_g[2])
p0 = [550.0,10.,1.]
popt_l, pcov = curve_fit(lorentzian, freq, data, p0 = p0)
lorentz_model = lorentzian(freq,popt_l[0],popt_l[1],popt_l[2])
fig4 = plt.figure(figsize = (12,6))
ax4 = plt.subplot(111)
ax4.set_xlim((np.min(freq),np.max(freq)))
ax4.plot(freq,data,alpha = 0.5);
ax4.plot(freq,gauss_model,'r-');
ax4.plot(freq,lorentz_model,'r-.');
ax4.legend(('Data','Gaussian','Lorentzian'));
print('Central frequency and standard deviation (THz):')
print(popt_g[0],abs(popt_g[1]))
def temperature(mu,sigma):
c = 3*10**8
m_H = 1.67 * 10**(-27)
k_B = 1.38 * 10**(-23)
T = (m_H*(c**2)/k_B)*(sigma/mu)**2
return T
print('Temperature of gas (K):')
print(temperature(popt_g[0],popt_g[1]))
Central frequency and standard deviation (THz): (564.99947400407166, 0.0097453843618178701) Temperature of gas (K): 3240.27412104