In [1]:
#import external modules
import os
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
matplotlib.rcParams['font.size'] = 8
matplotlib.rcParams['savefig.dpi'] = 1.5 * matplotlib.rcParams['savefig.dpi']
from scipy.optimize import leastsq
import WrightTools as wt
box_path = wt.kit.get_box_path()
In [2]:
#get absorbance data
abs_data_path = os.path.join(box_path, 'MX2', 'outside spectra', 'thin film abs.p')
abs_data = wt.data.from_pickle(abs_data_path)
abs_data.convert('wn')
cutoff = 325
xi = abs_data.axes[0].points[cutoff:]
zi = abs_data.channels[0].values[cutoff:]
plt.plot(xi, zi)
plt.grid()
plt.xlim(12500, 19500)
plt.title('raw absorbance')
plt.xlabel('wn')
Out[2]:
We need to lightly smooth the data in order to get a nice second derivative.
In [3]:
#smooth data
xi_raw = xi.copy()
zi_raw = zi.copy()
dat1 = np.array([xi, zi])
n = 10
for i in range(n, len(xi)-n):
#change the x value to the average
window = dat1[1][i-n:i+n].copy()
dat1[1][i] = window.mean()
xi, zi = dat1[:][:,n:-n]
plt.plot(xi_raw, zi_raw)
plt.plot(xi, zi)
plt.grid()
plt.xlim(12500, 19500)
plt.xlabel('wn')
plt.legend(['before smoothing', 'after smoothing'], loc = 'upper left')
Out[3]:
We can now take the second derivative of both smoothed and raw data to see what our smoothing work accomplished.
In [4]:
#differentiate data
xi_diff_raw, zi_diff_raw = wt.kit.diff(xi_raw, zi_raw, order = 2)
xi_diff, zi_diff = wt.kit.diff(xi, zi, order = 2)
print zi_diff.shape
print zi.shape
plt.plot(xi_raw, zi_diff_raw)
plt.plot(xi, zi_diff)
plt.grid()
plt.xlim(12500, 19500)
plt.xlabel('wn')
plt.legend(['before smoothing', 'after smoothing'], loc = 'upper left')
Out[4]:
Now lets try to fit the data.
In [19]:
#define a 2nd derivitive gaussian function
def gauss(x, p):
# p = amp, mu, sigma
a, mu, sigma = p
return a*np.exp(-(x-mu)**2 / (2*np.abs(sigma)**2))
def d2gauss(x, p):
a, mu, sigma = p
g = a*np.exp(-(x-mu)**2 / (2*np.abs(sigma)**2))
diff = wt.kit.diff(x, g, order = 2)[1]
return np.nan_to_num(diff)
#error function is sum of two such functions
f_err = lambda p, x, y: y - (d2gauss(x,p[:3]) + d2gauss(x,p[3:]))
# initial guesses
A_guess = [0.015, 14400, 250]
B_guess = [0.010, 15700, 300]
p0 = A_guess + B_guess
#fit
trunc = 10
out = leastsq(f_err, p0, args=(xi[:-trunc], zi_diff[:-trunc]))
p = out[0]
print p[:3]
print p[3:]
zfit = d2gauss(xi,p[:3]) + d2gauss(xi,p[3:])
plt.plot(xi, zi_diff)
plt.plot(xi, zfit, lw = 2)
plt.axvline(p[1], color = 'r')
plt.axvline(p[4], color = 'c')
plt.grid()
plt.xlim(12500, 19500)
plt.xlabel('wn')
plt.legend(['data', 'fit', 'A', 'B'], loc = 'lower right')
Out[19]:
In [20]:
#residual
residual = zi_diff - zfit
plt.plot(xi_diff, zi_diff)
plt.plot(xi_diff, residual, lw = 1)
plt.axvline(p[1], color = 'r')
plt.axvline(p[4], color = 'c')
plt.grid()
plt.xlim(12500, 19500)
plt.xlabel('wn')
plt.legend(['data', 'residual', 'A', 'B'], loc = 'lower right')
Out[20]:
In [ ]:
zi_A =