In [17]:
try:
reload(apogee.spec.lsf)
except NameError:
import apogee.spec.lsf
import apogee.modelspec.turbospec
from apogee.modelatm import atlas9
from apogee.spec.plot import apStarWavegrid
import apogee.spec.plot
import time
apogee.spec.lsf contains a LSF deconvolution technique that efficiently solves the $\chi^2$ minimization of a spectrum a 3x higher resolution than the native apStar resolution
$\chi^2 = \sum_\lambda (o_\lambda - \sum_\bar\lambda LSF_{\lambda\,\bar{\lambda}}\cdot o_\bar\lambda)^2/\sigma_\lambda + \epsilon\,\sum_\tilde{\lambda} (m_\tilde{\lambda}-m_{\tilde{\lambda}+\mathrm{d}\tilde{\lambda}})^2\,.$
That is, the observed spectrum $o_\lambda$ is modeled as the convolution of the model spectrum $m_\tilde{\lambda}$, which lives at 3x higher resolution, with the LSF (first term) and differences between adjacent pixels are penalized by the second term. $\epsilon$ is a free parameter that sets the smoothness of the high-resolution model spectrum. Ideally this should probably be a function of wavelength and perhaps even $T_{\mathrm{eff}}$ etc.
Following Nick Troup, we will generate a model spectrum with $T_{\mathrm{eff}} = 4000$ K, $\log g = 1$ and solar abundances using turbospectrum, convolve it with the LSF of fiber 80 (instead of 150, because 80 is more different from the combo LSF)[, and add S/N=100 errors, commented out]
In [101]:
x= numpy.linspace(-7.,7.,43)
lsf80= apogee.spec.lsf.eval(x,fiber=80)
In [107]:
atm= atlas9.Atlas9Atmosphere(teff=4000.,logg=1.,metals=0.0,am=0.0,cm=0.0)
synspec80= apogee.modelspec.turbospec.synth(modelatm=atm,
linelist='turbospec.201312161124',Hlinelist='Hlinedata',air=True,
xlsf=x,lsf=lsf80,cont='aspcap',vmacro=6.,isotopes='arcturus')
#synspec80[0]+= numpy.random.normal(size=len(synspec150[0]))*0.01
We also generate the 'goal' spectrum, which is convolved with the combo LSF
In [8]:
synspec= apogee.modelspec.turbospec.synth(modelatm=atm,
linelist='turbospec.201312161124',Hlinelist='Hlinedata',air=True,
lsf='combo',cont='aspcap',vmacro=6.,isotopes='arcturus')
We deconvolve the synthetic spectrum that was convolved with the LSF from fiber 80 (using very small uncertainties)
In [130]:
start= time.time()
synspec80_deconvolved= apogee.spec.lsf.deconvolve(synspec80[0],numpy.ones_like(synspec80[0])*0.01,lsf=lsf150,
eps=100.)
print "Deconvolution took {0:.2f} seconds".format(time.time()-start)
We then reconvolve with the combo LSF
In [131]:
# Need to add this somewhere
dx= x[1]-x[0]
hires= int(round(1./dx))
l10wav= numpy.log10(apStarWavegrid())
dowav= l10wav[1]-l10wav[0]
wav= 10.**numpy.arange(l10wav[0],l10wav[-1]+dowav/hires,dowav/hires)
xc, lsfcombo= apogee.spec.lsf._load_precomp(fiber='combo')
synspec80_reconvolved= apogee.spec.lsf.convolve(wav,synspec80_deconvolved,lsf=lsfcombo,xlsf=xc)
Now we make some plots
In [132]:
for panel in apogee.spec.plot.highres(synspec80_deconvolved[::3],synspec80[0],synspec80_reconvolved[0],synspec[0],
(synspec80[0]-synspec[0])*10.+0.25,
(synspec80_reconvolved[0]-synspec[0])*10.+0.25,
yrange=[0.1,1.2],
color=['c','r','g','b','k','y'],labelLines=True,cleanZero=False,fig_width=14.):
show()
The difference between the spectra
In [133]:
print "80-combo = {0:.4f}+/-{1:.4f}, re/deconvolved 80-combo = {2:.4f}+/-{3:.4f}".format(numpy.nanmedian(synspec80[0]-synspec[0]),
numpy.nanstd(synspec80[0]-synspec[0]),
numpy.nanmedian(synspec80_reconvolved[0]-synspec[0]),
numpy.nanstd(synspec80_reconvolved[0]-synspec[0]))
print "80-combo range = {0:.4f}, re/deconvolved 80-combo range = {1:.4f}".format(numpy.ptp((synspec80[0]-synspec[0])[True-numpy.isnan(synspec150[0]-synspec[0])]),
numpy.ptp((synspec80_reconvolved[0]-synspec[0])[True-numpy.isnan(synspec150_reconvolved[0]-synspec[0])]))
In [ ]: