lineshape extraction



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')









    



data opened from C:\Users\Blaise\Box Sync\Wright Shared\MX2\outside spectra\thin film abs.p
axis wm converted






    Out[2]:





<matplotlib.text.Text at 0x86aa1d0>

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]:





<matplotlib.legend.Legend at 0x882afb0>

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')









    



(474,)
(474,)






    Out[4]:





<matplotlib.legend.Legend at 0x8946970>

Now lets try to fit the data. This time using gaussians.



In [5]:

    
#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):
    g =  gauss(x, p)
    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')









    



[  1.81055001e-02   1.44717633e+04   3.03608091e+02]
[  1.89677339e-02   1.57011612e+04   4.38343173e+02]






    Out[5]:





<matplotlib.legend.Legend at 0x8ba01b0>



In [6]:

    
#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[6]:





<matplotlib.legend.Legend at 0x8c98150>



In [7]:

    
zi_A = gauss(xi, p[:3])
zi_B = gauss(xi, p[3:])
residual = zi - (zi_A + zi_B)

plt.plot(xi, zi_A, color = 'r', lw = 2)
plt.plot(xi, zi_B, color = 'b', lw = 2)
plt.plot(xi, zi, color = 'k', lw = 2)
plt.plot(xi, residual, 'k:', lw = 2)
plt.axvline(p[1], color = 'r', lw = 1)
plt.axvline(p[4], color = 'b', lw = 1)
plt.grid()
plt.xlabel('wn')
plt.savefig('gauss.png')

Now lorentzians.



In [8]:

    
#define a 2nd derivitive lorentzian function
def lorentz(x, p):
    # p = amp, center, FWHM
    a, center, FWHM = p
    x_var = (center - x) / (0.5*FWHM)
    return a/(1+x_var**2)
    
def d2lorentz(x, p):
    g =  lorentz(x, p)
    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 - (d2lorentz(x,p[:3]) + d2lorentz(x,p[3:]))

# initial guesses
A_guess = [0.005, 14500, 300]
B_guess = [0.005, 15500, 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 = d2lorentz(xi,p[:3]) + d2lorentz(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')









    



[  2.96701282e-02   1.44609978e+04   1.00182453e+03]
[  3.75533913e-02   1.56997728e+04   1.49161907e+03]






    Out[8]:





<matplotlib.legend.Legend at 0x8ed5cb0>



In [9]:

    
#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[9]:





<matplotlib.legend.Legend at 0x903f330>



In [10]:

    
zi_A = lorentz(xi, p[:3])
zi_B = lorentz(xi, p[3:])
residual = zi - (zi_A + zi_B)

plt.plot(xi, zi_A, color = 'r', lw = 2)
plt.plot(xi, zi_B, color = 'b', lw = 2)
plt.plot(xi, zi, color = 'k', lw = 2)
plt.plot(xi, residual, 'k:', lw = 2)
plt.axvline(p[1], color = 'r', lw = 1)
plt.axvline(p[4], color = 'b', lw = 1)
plt.grid()
plt.xlabel('wn')
plt.savefig('lorentz.png')