In [79]:
try:
reload(apogee.modelatm.atlas9)
reload(apogee.modelspec.moog)
reload(apogee.modelspec.ferre)
except NameError:
import apogee.modelatm.atlas9
import apogee.modelspec.moog
import apogee.modelspec.ferre
import apogee.spec.cannon as apcannon
import apogee.spec.plot as splot
import apogee.spec.window as apwindow
from apogee.tools import paramIndx, elemIndx
from matplotlib import cm
from galpy.util import bovy_plot
Convenience function to generate a model spectrum
In [2]:
def simspec(ofe=0.,teff=4000.,logg=2.5,metals=0.,cfe=0.,nfe=0.,isotopes='solar',solarscaledatm=False):
# Set up a model atmosphere for this point
if solarscaledatm:
atm= apogee.modelatm.atlas9.Atlas9Atmosphere(teff=teff,logg=logg,metals=metals,am=0.,cm=0.)
abu= ([8,ofe],)
else:
atm= apogee.modelatm.atlas9.Atlas9Atmosphere(teff=teff,logg=logg,metals=metals,am=ofe,cm=cfe)
# Keep the other alpha elements the same
abu= ([8,0.],[12,-ofe],[14,-ofe],[16,-ofe],[20,-ofe],[22,-ofe])
synspec= apogee.modelspec.moog.synth(*abu,
modelatm=atm,linelist='moog.201405291900.vac',lsf='combo',
cont='aspcap',vmacro=6.,isotopes=isotopes,
vmicro=2.478-0.325*logg)
return synspec
Convenience function to analyze a simulated spectrum using the MOOG library
In [3]:
def fit(synspec,initcannon=True,blockOwindows=False):
# initcannon=True only runs a single FERRE fit initialized close to the best fit
# blockOwindows=True removes the O windows from the fparam fit
# Assume SNR=200
errs= 0.005*numpy.ones_like(synspec)
if blockOwindows:
errs[apwindow.tophat('O',apStarWavegrid=True)]= 10.
fparam= apogee.modelspec.ferre.fit(synspec,errs,initcannon=initcannon,lib='msGK',dr='current')
felem= apogee.modelspec.ferre.elemfitall(synspec,0.005*numpy.ones_like(synspec),fparam=fparam,lib='msGK',dr='current')
return (fparam,felem)
In [80]:
def cannon(synspec,blockOwindows=False):
# blockOwindows=True removes the O windows from the analysis
# Assume SNR=200
errs= 0.005*numpy.ones_like(synspec)
if blockOwindows:
errs[apwindow.tophat('O',apStarWavegrid=True)]= 10.
labels= apcannon.polylabels(synspec,errs)
return labels
A function to plot the results
In [100]:
def plot_results(var,fparam,felem,fparam_ref,labels,var_label=r'$\mathrm{O}$'):
figsize(18,20)
# First do the parameters
params= ['TEFF','LOGG','METALS','C','N','ALPHA']
param_labels= [r'$\Delta T_{\mathrm{eff}}\,(\mathrm{K})$',r'$\Delta \log g$',r'$\Delta [\mathrm{M/H}]$',
r'$\Delta [\mathrm{C/M}]$',r'$\Delta [\mathrm{N/M}$]',r'$\Delta [\alpha/\mathrm{M}$]']
cannonlabelIndx= [0,1,2,None,None,3]
ranges= [[-300.,300.],[-1.,1.],[-0.4,0.4],[-0.4,0.4],[-0.4,0.4],[-0.4,0.4]]
for ii, param in enumerate(params):
subplot(7,3,ii+1)
ty= fparam[:,paramIndx(param)]-fparam_ref[paramIndx(param)]
tdy= numpy.amax(ty)-numpy.amin(ty)
yrange= ranges[ii]
bovy_plot.bovy_plot(var,ty,'ko',xlabel=var_label,ylabel=param_labels[ii],
xrange=[var[0]-0.1,var[-1]+0.1],yrange=yrange,gcf=True)
if not cannonlabelIndx[ii] is None:
ty= labels[:,cannonlabelIndx[ii]]-fparam_ref[paramIndx(param)]
bovy_plot.bovy_plot(var,ty,'ro',overplot=True)
# Now do the elements
elems= ['C','N','O','Na','Mg','Al','Si','S','K','Ca','Ti','V','Mn','Fe','Ni']
for ii, elem in enumerate(elems):
subplot(7,3,ii+7)
ty= felem[elem]
if not elem.lower() == 'o':
ylabel= r'$\Delta \mathrm{%s}$' % elem
else:
ylabel= r'$\mathrm{Recovered\ %s}$' % elem
bovy_plot.bovy_plot(var,ty,'ko',xlabel=var_label,ylabel=ylabel,
xrange=[var[0]-0.1,var[-1]+0.1],yrange=[-0.4,0.4],gcf=True)
#tight_layout()
Let's generate a series of spectra at 4000 K
In [5]:
teff= 4000.
ofes= [-0.3,-0.2,-0.1,0, 0.1, 0.2, 0.3]
synspecs40= []
for ofe in ofes:
synspecs40.append(simspec(ofe=ofe,teff=teff)[0])
Let's take a look at these spectra
In [7]:
colormap= cm.coolwarm
def color(ofe):
return colormap((ofe-ofes[0])/(ofes[-1]-ofes[0]))
splot.windows(synspecs40[0],'O1',color=color(ofes[0]))
for ii in range(1,len(synspecs40)): splot.windows(synspecs40[ii],'O1',color=color(ofes[ii]),overplot=True)
Now we can analyze these simulated spectra
In [8]:
results_fparam= []
results_felem= []
for ii in range(len(ofes)):
tresult= fit(synspecs40[ii])
results_fparam.append(tresult[0])
results_felem.append(tresult[1])
In [94]:
results_cannon= []
for ii in range(len(ofes)):
results_cannon.append(cannon(synspecs40[ii]))
In [90]:
fparam_ref= [4000.,2.5,0.,0.,0.,0.,0.]
results_fparam= numpy.reshape(numpy.array(results_fparam),(len(ofes),7))
newresults_felem= {}
for elem in results_felem[0]: newresults_felem[elem]= numpy.array([results_felem[ii][elem] for ii in range(len(ofes))])
plot_results(ofes,results_fparam,newresults_felem,fparam_ref,numpy.reshape(numpy.array(results_cannon),(len(ofes),4)),
var_label=r'$\mathrm{O}$')
Same for a solar-scaled atmosphere (the model Jon ran)
In [10]:
teff= 4000.
ofes= [-0.3,-0.2,-0.1,0, 0.1, 0.2, 0.3]
synspecs40ss= []
for ofe in ofes:
synspecs40ss.append(simspec(ofe=ofe,teff=teff,solarscaledatm=True)[0])
In [11]:
results_fparam_ss= []
results_felem_ss= []
for ii in range(len(ofes)):
tresult= fit(synspecs40ss[ii])
results_fparam_ss.append(tresult[0])
results_felem_ss.append(tresult[1])
In [92]:
results_cannon_ss= []
for ii in range(len(ofes)):
results_cannon_ss.append(cannon(synspecs40ss[ii]))
In [93]:
fparam_ref= [4000.,2.5,0.,0.,0.,0.,0.]
results_fparam_ss= numpy.reshape(numpy.array(results_fparam_ss),(len(ofes),7))
newresults_felem_ss= {}
for elem in results_felem_ss[0]: newresults_felem_ss[elem]= numpy.array([results_felem_ss[ii][elem] for ii in range(len(ofes))])
plot_results(ofes,results_fparam_ss,newresults_felem_ss,fparam_ref,numpy.reshape(numpy.array(results_cannon_ss),(len(ofes),4)),
var_label=r'$\mathrm{O}$')
In [13]:
teff= 4500.
ofes= [-0.3,-0.2,-0.1,0, 0.1, 0.2, 0.3]
synspecs45= []
for ofe in ofes:
synspecs45.append(simspec(ofe=ofe,teff=teff)[0])
In [14]:
results_fparam_45= []
results_felem_45= []
for ii in range(len(ofes)):
tresult= fit(synspecs45[ii])
results_fparam_45.append(tresult[0])
results_felem_45.append(tresult[1])
In [97]:
results_cannon_45= []
for ii in range(len(ofes)):
results_cannon_45.append(cannon(synspecs45[ii]))
In [98]:
fparam_ref= [4500.,2.5,0.,0.,0.,0.,0.]
results_fparam_45= numpy.reshape(numpy.array(results_fparam_45),(len(ofes),7))
newresults_felem_45= {}
for elem in results_felem_45[0]: newresults_felem_45[elem]= numpy.array([results_felem_45[ii][elem] for ii in range(len(ofes))])
plot_results(ofes,results_fparam_45,newresults_felem_45,fparam_ref,numpy.reshape(numpy.array(results_cannon_45),(len(ofes),4)),
var_label=r'$\mathrm{O}$')
Block out the [O/H] windows in the fparam fit
In [28]:
results_fparam_45bO= []
results_felem_45bO= []
for ii in range(len(ofes)):
tresult= fit(synspecs45[ii],blockOwindows=True)
results_fparam_45bO.append(tresult[0])
results_felem_45bO.append(tresult[1])
In [95]:
results_cannon_45bO= []
for ii in range(len(ofes)):
results_cannon_45bO.append(cannon(synspecs45[ii],blockOwindows=True))
In [99]:
fparam_ref= [4500.,2.5,0.,0.,0.,0.,0.]
results_fparam_45bO= numpy.reshape(numpy.array(results_fparam_45bO),(len(ofes),7))
newresults_felem_45bO= {}
for elem in results_felem_45bO[0]: newresults_felem_45bO[elem]= numpy.array([results_felem_45bO[ii][elem] for ii in range(len(ofes))])
plot_results(ofes,results_fparam_45bO,newresults_felem_45bO,fparam_ref,numpy.reshape(numpy.array(results_cannon_45bO),(len(ofes),4)),
var_label=r'$\mathrm{O}$')
Let's do MCMC on the lowest [O/Fe] synthetic spectrum:
In [26]:
chain= apogee.modelspec.ferre.mcmc(synspecs45[0],0.005*numpy.ones_like(synspecs45[0]),
initsig=[100.,0.25,0.1,0.1,0.1,0.1,0.1],
fparam=numpy.reshape(results_fparam_45[0],(1,7)),nsamples=100)
In [24]:
figsize(8,8)
bovy_plot.scatterplot(chain[:,0],chain[:,-1],'k,')
In [27]:
figsize(8,8)
bovy_plot.scatterplot(chain[:,0],chain[:,-1],'k,')
So, the $T_{\mathrm{eff}} =4600$ K really seems to be the better fit. Let's compare the best-fit spectrum and the input fparam to the data
In [64]:
ispec= apogee.modelspec.ferre.interpolate([4500.,results_fparam_45[0][0]],[2.5,results_fparam_45[0][1]],
[0.,results_fparam_45[0][3]],[-0.2,results_fparam_45[0][6]],
[0.,results_fparam_45[0][5]],[0.,results_fparam_45[0][4]],
lib='msGK',dr='current',apStarWavegrid=True)
In [65]:
print numpy.nanstd(ispec[1]-synspecs45[0])
print numpy.nanstd(ispec[0]-synspecs45[0])
The first scatter is smaller, so $T_{\mathrm{eff}} = 4600$ K is really the better fit. Below the next two plots I compute the scatter for the well-recovered case of [O/H] = 0, which shows that this level of residual scatter appears typical. A comparison of the spectra and the data follows (blue is the best-fit, green is with the correct $T_{\mathrm{eff}} = 4600$):
In [66]:
splot.waveregions(ispec[1]-synspecs45[0],yrange=[-0.05,0.05],cleanZero=False,labelLines=False,ylabel=r'$\Delta f/f_c$')
splot.waveregions(ispec[0]-synspecs45[0],cleanZero=False,overplot=True)
In [67]:
splot.windows(ispec[1]-synspecs45[0],'O1',yrange=[-0.05,0.05],cleanZero=False,labelLines=False,ylabel=r'$\Delta f/f_c$')
splot.windows(ispec[0]-synspecs45[0],'O1',cleanZero=False,overplot=True)
For reference, here's what this looks like for a good fit (index 4 is the zero [O/H] offset):
In [68]:
ispecgf= apogee.modelspec.ferre.interpolate([4500.,results_fparam_45[4][0]],[2.5,results_fparam_45[4][1]],
[0.,results_fparam_45[4][3]],[0.,results_fparam_45[4][6]],
[0.,results_fparam_45[4][5]],[0.,results_fparam_45[4][4]],
lib='msGK',dr='current',apStarWavegrid=True)
In [69]:
print numpy.nanstd(ispecgf[1]-synspecs45[4])
print numpy.nanstd(ispecgf[0]-synspecs45[4])
In [70]:
splot.waveregions(ispecgf[1]-synspecs45[4],yrange=[-0.05,0.05],cleanZero=False,labelLines=False,ylabel=r'$\Delta f/f_c$')
splot.waveregions(ispecgf[0]-synspecs45[4],cleanZero=False,overplot=True)
In [71]:
splot.windows(ispecgf[1]-synspecs45[4],'O1',yrange=[-0.05,0.05],cleanZero=False,labelLines=False,ylabel=r'$\Delta f/f_c$')
splot.windows(ispecgf[0]-synspecs45[4],'O1',cleanZero=False,overplot=True)
In [ ]: