When one observes the spectrum of a source of light, the x-axis of the spectrum is in pixels. To convert this to wavelengths, it is used a spectrum of a well know source of light, e.g., the light of a lamp of ThAr. The spectrum of the lamp has known emission lines, thus it is possible to convert from pixel to wavelngth.
From theory, the emission lines of the lamp spectrum should be infinitely narrow. Unfortunately, this does not happens because of the optics of the instrument, which broadens the spectrum. When studying spectral lines of unknown sources, this broadeing effect has to be considered.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from instprofile import *
I have copied a spectrum of a real lamp and it can be seen in the plot below.
In [2]:
lamp_spec = np.loadtxt('lamp.dat')
fig = plt.figure(figsize=[15,9])
ax = fig.add_subplot(111)
ax.plot(lamp_spec[:,0], lamp_spec[:,1])
ax.set_xlabel('Wavelength ($\AA$)')
ax.set_ylabel('Flux')
_ = ax.set_ylim([0,140000])
Let's first start with a slice of this spectrum and try to see if it detects the peaks. The value of the delta threhsold in flux is very important. It defines which peaks are chosen. If the values increases, only the strong lines are picked.
In [3]:
cond1 = lamp_spec[:,0] > 3725
cond2 = lamp_spec[:,0] < 3735
lamp_slice = lamp_spec[cond1 & cond2]
delta_threshold = 5000 # Threhold in flux between peaks
peaks = peakdet(lamp_slice[:,1], delta_threshold, lamp_slice[:,0])[0]
#Fit gaussians
gauss = lambda p, x: p[0]*np.exp(-(p[1]-x)**2/(2*p[2]**2)) #1d Gaussian func
y_vector = np.zeros_like(lamp_slice[:,1])
for param in find_fit(lamp_slice[:,0], lamp_slice[:,1], delta_threshold):
y_vector += gauss([param[0], param[1], param[2]], lamp_slice[:,0])
# Plot
fig = plt.figure(figsize=[15,9])
ax = fig.add_subplot(111)
ax.plot(lamp_slice[:,0], lamp_slice[:,1], label='Subset of the lamp spectrum')
ax.plot(peaks[:,0], peaks[:,1], 'ro', label ='Peaks found')
ax.plot(lamp_slice[:,0], y_vector, 'g--',linewidth=3, label='Gaussian fit')
ax.set_xlabel('Wavelength ($\AA$)')
ax.set_ylabel('Flux')
_ = ax.legend()
In [4]:
gauss = lambda p, x: p[0]*np.exp(-(p[1]-x)**2/(2*p[2]**2)) #1d Gaussian func
y_vector = np.zeros_like(lamp_slice[:,1])
for param in find_fit(lamp_slice[:,0], lamp_slice[:,1], delta_threshold):
y_vector += gauss([param[0], param[1], param[2]], lamp_slice[:,0])
The detection is OK. We can now obtain the instrumental profile of the entire spectrum.
In [5]:
iprof = inst_profile(lamp_spec[:,0], lamp_spec[:,1], delta_threshold)
fig = plt.figure(figsize=[15,9])
ax = fig.add_subplot(111)
ax.plot(iprof[:,0], iprof[:,1], 'gs')
ax.set_xlabel('Wavelength ($\AA$)')
ax.set_ylabel('FWHM')
_ = ax.set_ylim([0,55])
But note that it is not perfect and can obtain very strange values of FWHM. We can manually set an upper limit.
In [6]:
iprof = inst_profile(lamp_spec[:,0], lamp_spec[:,1], delta_threshold, 0.25)
fig = plt.figure(figsize=[15,9])
ax = fig.add_subplot(111)
ax.plot(iprof[:,0], iprof[:,1], 'gs')
ax.hlines(0.25, 3710, 3790, 'r', '--')
ax.text(3775, 0.255, 'Upper limit', fontsize=16, color='red')
ax.set_xlabel('Wavelength ($\AA$)')
ax.set_ylabel('FWHM')
_ = ax.set_ylim([0,0.3])