# Demo entry 6351555

ISMLAB1

Submitted by anonymous on Mar 21, 2017 at 17:47
Language: Python. Code size: 6.5 kB.

```import numpy as np
import os,sys
from scipy.interpolate import interp1d
from scipy.interpolate import UnivariateSpline
from scipy.optimize import curve_fit
from scipy.optimize import fsolve
from scipy.special import wofz
import matplotlib.pyplot as pl
import matplotlib.cm as cm
from matplotlib import rc, rcParams
rcParams.update({'font.size': 20})
from itertools import cycle
from astropy.table import Table
import astropy.io.fits as pf
from astropy.cosmology import WMAP9 as cosmo

# 1. (a)

# Pure CO
ks = data[:,0] # cm^-1
absorbs = data[:,1] # 1

# mask the data to remove the second bump
masks = [ks > 2132, ks < 2141]

# a Guassian model
def gaus(x, a, x0, sigma):
return a*np.exp(-(x-x0)**2/2.0/sigma**2)

# a Lorentzian model
def lorentz(x, x0, gamma):
return gamma / np.pi / ((x-x0)**2 + gamma**2)

# a Voigt model
def voigt(x, a, x0, alpha, gamma):

sigma = alpha / np.sqrt(2 * np.log(2))

return a*np.real(wofz(((x-x0) + 1j*gamma)/sigma/np.sqrt(2))) / sigma /np.sqrt(2*np.pi)

# we fit the main absorption feature to a Voigt profile
popt, pcov = curve_fit(voigt, ks_fit, absorbs_fit, p0=[0.04,2139, 2, 2])

upperlim = 2147
lowerlim = 2132
"""
# plot the data and the fit
newks = np.linspace(lowerlim, upperlim, num = 500)
pl.plot(ks, absorbs, 'ro')
pl.plot(newks, voigt(newks, *popt), 'b-')
pl.xlabel('wavenumber (cm\$^{-1}\$)')
pl.ylabel('absorbance')
pl.xlim(lowerlim, upperlim)
pl.tight_layout()
pl.show()
"""

# integrate the curve_fit results and dose
int_results = quad(lambda x: voigt(x, *popt), lowerlim, upperlim)
line_area = int_results[0] # cm^-1
A_i = 1.4e-17 # cm/molec
N_s = 1e15 #  cm^-2
Dose = line_area*2.3/A_i/N_s # molec
print Dose # 15.81

# CO + H2O
ks1 = data1[:,0]
absorbs1 = data1[:,1]

# mask the data to remove the second bump
masks = [ks1 > 2120, ks1 < 2145]

popt1, pcov1 = curve_fit(voigt, ks1_fit, absorbs1_fit, p0=[0.01,2139, 2, 2])

upperlim = 2170
lowerlim = 2120
"""
newks = np.linspace(lowerlim, upperlim, num = 500)
pl.plot(ks1, absorbs1, 'ro')
pl.plot(newks, voigt(newks, *popt1), 'b-')
pl.xlabel('wavenumber (cm\$^{-1}\$)')
pl.ylabel('absorbance')
pl.xlim(lowerlim, upperlim)
pl.tight_layout()
pl.show()
"""

# integrate the curve_fit results and dose
int_results = quad(lambda x: voigt(x, *popt1), lowerlim, upperlim)
line_area = int_results[0] # cm^-1
A_i = 1.4e-17 # cm/molec
N_s = 1e15 #  cm^-2
Dose = line_area*2.3/A_i/N_s # molec
print Dose # 14.11
CO_dens = line_area*2.3/A_i # molec/cm^2
print CO_dens # 1.41e16

"""
Thus, our CO dose for the pure CO experiment is approximately 16.
Our CO dose for the CO+H2O experiment is approximately 14.
"""

# 1. (b)

# mask the part of the spectrum that is due to water absorption
masks = [ks1 > 2800, ks1 < 3560]

# do the fit
poptW, pcovW = curve_fit(voigt, ks1_fit, absorbs1_fit, p0=[0.03,3260, 200, 200])

# integrate the curve_fit results and density
upperlim = 3800
lowerlim = 2800
int_results = quad(lambda x: voigt(x, *poptW), lowerlim, upperlim)
line_area = int_results[0] # cm^-1
A_i = 2.2e-16 # cm/molec
H2O_dens = 2.3*line_area/A_i # molec/cm^2
print H2O_dens # 7.93e+16
# then the mixing ratio is
mix_ratio = CO_dens/H2O_dens
print mix_ratio # 0.178

"""
newks = np.linspace(lowerlim, upperlim, num = 500)
pl.plot(ks1, absorbs1, 'ro')
pl.plot(newks, voigt(newks, *poptW), 'b-')
pl.xlabel('wavenumber (cm\$^{-1}\$)')
pl.ylabel('absorbance')
pl.xlim(lowerlim, upperlim)
pl.tight_layout()
pl.show()
"""

"""
Thus, the CO/H2O mixing ratio is 0.178.
"""

# 1. (c) Calculate CO peak FWHM for both experiments

# pure CO
x0 = popt[1]
max_val = voigt(x0, *popt)

# define an equation to solve for the xs of fwhm
def myequ(x):
return voigt(x, *popt) - max_val/2

# use fsolve to solve the equation
x1 = fsolve(myequ, 2137)[0]
x2 = fsolve(myequ, 2139)[0]
print x1, x2 # 2137.29908082, 2139.26025961
fwhm = x2 - x1
print fwhm # 1.96117878759

# CO + H2O
x0 = popt1[1]
max_val = voigt(x0, *popt1)

# define an equation to solve for the xs of fwhm
def myequ(x):
return voigt(x, *popt1) - max_val/2

# use fsolve to solve the equation
x1 = fsolve(myequ, 2137)[0]
x2 = fsolve(myequ, 2139)[0]
print x1, x2 # 2132.21780321 2141.8357287
fwhm1 = x2 - x1
print fwhm1 # 9.61792548731

"""
Thus, the FWHMs of the CO absorption feature are 1.96 cm^-1
and 9.62 cm^-1 for the two experiments respectively.
The different in FWHM can be at least in part be attributed
to environmental effects due to the ice environment of the
second experiment, which change the energy level of the
absorption through diverse chemical bonds and many other
environmental effects.
"""

# 1. (d)

# pure CO

upperlim = 2150
lowerlim = 2130
# plot the data and the fit
pl.figure(1)
lines = ["-","--","-.",":"]
linecycler = cycle(lines)
colors = ["b", "y", "g", "k", "m"]
colorcycler = cycle(colors)

# plot the original
pl.plot(ks, absorbs, color='r', ls = next(linecycler), label='12K')

for T in np.linspace(15, 33, num = 10):
T = int(T)
ks = data[:,0] # cm^-1
absorbs = data[:,1] # 1

# calculate background
mask = [ks > 2130, ks < 2133]

pl.plot(ks, absorbs - absorbs_bg, color=next(colorcycler), ls = next(linecycler), label = str(T)+"K")

pl.xlabel('wavenumber (cm\$^{-1}\$)')
pl.ylabel('absorbance')
pl.xlim(lowerlim, upperlim)
pl.legend(loc='best',prop={'size':20})
pl.tight_layout()
pl.show()

"""
We see that as the temperature increases, the magnitude of the absorption also decreases.
This is due to the sublimation of CO, which leads to the reduction of CO surface density
along the IR LOS.
The width of the absorption feature also decreases, and I think this is due to the more
CO layers we have, the more inter-molecular energy structure we have, which leads to line
broadening. As we lose the layers, the energy levels get removed, thinning out the absorption
features.
"""
```

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.