

In [1]:

    
try:
    reload(apogee.modelspec.turbospec)
    reload(apogee.modelspec.moog)
except NameError:
    import apogee.modelspec.turbospec
    import apogee.modelspec.moog
import apogee.modelspec.ferre
from apogee.modelatm import atlas9
from apogee.spec.plot import apStarWavegrid
import apogee.spec.plot as splot
import tempfile
import os
from galpy.util import bovy_plot

A comparison between Kurucz/Turbospectrum, Kurucz/MOOG, and Kurucz/ASSeT

Comparison between MOOG and Turbospectrum (true-)continuum-normalized spectra at high-resolution



In [2]:

    
teffs= [3500,4500]
loggs= [0,2,4]



In [3]:

    
moogspecs_hires= []
for teff in teffs:
    for logg in loggs:
        atm= atlas9.Atlas9Atmosphere(teff=teff,logg=logg)
        atmFile= tempfile.mktemp(dir=os.getcwd())+'.mod'
        atm.writeto(atmFile)
        try:
            moogspecs_hires.append(apogee.modelspec.moog.moogsynth(modelatm=atmFile,linelist='moog.201312161124.vac',
                                                                   wmin=15150.,wmax=16930.,dw=0.02,vmicro=2.,
                                                                   isotopes='solar')[1][0])
        finally:
            os.remove(atmFile)



In [4]:

    
turbospecs_hires= []
for teff in teffs:
    for logg in loggs:
        atm= atlas9.Atlas9Atmosphere(teff=teff,logg=logg)
        turbospecs_hires.append(apogee.modelspec.turbospec.turbosynth(modelatm=atm,linelist='turbospec.201312161124.new.vac',
                                                                      Hlinelist='Hlinedata.vac',
                                                                      wmin=15150.,wmax=16930.,dw=0.02,vmicro=2.,
                                                                      isotopes='solar')[1])



In [5]:

    
#Exclude region around Hydrogen lines
wav= numpy.linspace(15150.,16930.,89001)
hlines= [15004.970,15043.157,15086.906,15137.367,15196.005,15264.717,
             15345.992,15443.148,15560.708,15704.960,15884.888,16113.721,
             16411.681,16811.117]
cmpIndx= numpy.ones(len(wav),dtype='bool')
for hline in hlines:
    cmpIndx[numpy.fabs(wav-hline) < 10.]= False



In [6]:

    
cnt= 0
for teff in teffs:
    for logg in loggs:
        print teff, logg, numpy.nanmean((moogspecs_hires[cnt]-turbospecs_hires[cnt])[cmpIndx]), \
            numpy.nanstd((moogspecs_hires[cnt]-turbospecs_hires[cnt])[cmpIndx]), \
            numpy.nanmax((numpy.fabs(moogspecs_hires[cnt]-turbospecs_hires[cnt]))[cmpIndx])
        cnt+= 1









    



3500 0 -0.00590230753652 0.00698771920423 0.17033
3500 2 -0.00543194830511 0.00769627472406 0.10416
3500 4 -0.00481757765626 0.0115832688641 0.07014
4500 0 -0.00306272524746 0.0066908328459 0.29269
4500 2 -0.00248928508139 0.00401248852468 0.09589
4500 4 -0.00220162464241 0.00276247664386 0.05866

Comparison between Kurucz/Turbospectrum, Kurucz/MOOG, and Kurucz/ASSeT after LSF convolution and pseudo-continuum normalization



In [7]:

    
moogspecs= []
for teff in teffs:
    for logg in loggs:
        atm= atlas9.Atlas9Atmosphere(teff=teff,logg=logg)
        moogspecs.append(apogee.modelspec.moog.synth(modelatm=atm,linelist='moog.201312161124.vac',
                                                     vmicro=2.478-0.325*logg,
                                                     lsf='combo',cont='aspcap',vmacro=6.,isotopes='solar'))









    



/Library/Python/2.7/site-packages/numpy/core/fromnumeric.py:2507: VisibleDeprecationWarning: `rank` is deprecated; use the `ndim` attribute or function instead. To find the rank of a matrix see `numpy.linalg.matrix_rank`.
  VisibleDeprecationWarning)

/Users/bovy/Repos/apogee/apogee/modelspec/moog.py:179: RuntimeWarning: divide by zero encountered in divide
  out[cflux > 0.]/= cflux[cflux > 0.]



In [8]:

    
turbospecs= []
for teff in teffs:
    for logg in loggs:
        atm= atlas9.Atlas9Atmosphere(teff=teff,logg=logg)
        turbospecs.append(apogee.modelspec.turbospec.synth(modelatm=atm,linelist='turbospec.201312161124.new.vac',
                                                           Hlinelist='Hlinedata.vac',
                                                           vmicro=2.478-0.325*logg,
                                                           lsf='combo',cont='aspcap',vmacro=6.,isotopes='solar'))









    



/Users/bovy/Repos/apogee/apogee/modelspec/turbospec.py:154: RuntimeWarning: divide by zero encountered in divide
  out[cflux > 0.]/= cflux[cflux > 0.]



In [9]:

    
assetspecs= []
for teff in teffs:
    for logg in loggs:
        assetspecs.append(apogee.modelspec.ferre.interpolate(teff,logg,0.,0.,0.,0.))



In [10]:

    
#Exclude region around Hydrogen lines
wav= apStarWavegrid()
hlines= [15004.970,15043.157,15086.906,15137.367,15196.005,15264.717,
             15345.992,15443.148,15560.708,15704.960,15884.888,16113.721,
             16411.681,16811.117]
cmpIndx= numpy.ones(len(wav),dtype='bool')
for hline in hlines:
    cmpIndx[numpy.fabs(wav-hline) < 10.]= False



In [11]:

    
print "||= '''T,,eff,,''' =||= '''logg g''' =|||||||| '''ASSeT - Turbospec''' |||||||| '''ASSeT - MOOG''' |||||||| '''Turbospec - MOOG''' ||"
print '|| || || mean || std. dev. || robust std. dev. || max. diff || mean || std. dev. || robust std. dev. || max. diff || mean || std. dev. || robust std. dev. || max. diff ||' 
cnt= 0
for teff in teffs:
    for logg in loggs:
        print '|| {teff} || {logg} || {mat:.3f} || {sat:.3f} || {rat:.3f} || {dat:.3f} || {mam:.3f} || {sam:.3f} || {ram:.3f} || {dam:.3f} || {mtm:.3f} || {stm:.3f} || {rtm:.3f} || {dtm:.3f} ||'\
        .format(teff=teff,logg=logg,
                mat=numpy.nanmean((assetspecs[cnt]-turbospecs[cnt])[0,cmpIndx]),
                sat=numpy.nanstd((assetspecs[cnt]-turbospecs[cnt])[0,cmpIndx]),
                rat=1.4836*numpy.nanmedian((numpy.fabs(assetspecs[cnt]-turbospecs[cnt]-numpy.nanmedian((assetspecs[cnt]-turbospecs[cnt])[0,cmpIndx]))[0,cmpIndx])),
                dat=numpy.nanmax((numpy.fabs(assetspecs[cnt]-turbospecs[cnt]))[0,cmpIndx]),
                mam=numpy.nanmean((assetspecs[cnt]-moogspecs[cnt])[0,cmpIndx]),
                sam=numpy.nanstd((assetspecs[cnt]-moogspecs[cnt])[0,cmpIndx]),
                ram=1.4836*numpy.nanmedian((numpy.fabs(assetspecs[cnt]-moogspecs[cnt]-numpy.nanmedian((assetspecs[cnt]-moogspecs[cnt])[0,cmpIndx])))[0,cmpIndx]),
                dam=numpy.nanmax((numpy.fabs(assetspecs[cnt]-moogspecs[cnt]))[0,cmpIndx]),
                mtm=numpy.nanmean((turbospecs[cnt]-moogspecs[cnt])[0,cmpIndx]),
                stm=numpy.nanstd((turbospecs[cnt]-moogspecs[cnt])[0,cmpIndx]),
                rtm=1.4836*numpy.nanmedian((numpy.fabs(turbospecs[cnt]-moogspecs[cnt]-numpy.nanmedian((turbospecs[cnt]-moogspecs[cnt])[0,cmpIndx])))[0,cmpIndx]),
                dtm=numpy.nanmax((numpy.fabs(turbospecs[cnt]-moogspecs[cnt]))[0,cmpIndx]))
        cnt+= 1









    



||= '''T,,eff,,''' =||= '''logg g''' =|||||||| '''ASSeT - Turbospec''' |||||||| '''ASSeT - MOOG''' |||||||| '''Turbospec - MOOG''' ||
|| || || mean || std. dev. || robust std. dev. || max. diff || mean || std. dev. || robust std. dev. || max. diff || mean || std. dev. || robust std. dev. || max. diff ||
-c:15: RuntimeWarning: invalid value encountered in subtract

-c:16: RuntimeWarning: invalid value encountered in subtract

-c:17: RuntimeWarning: invalid value encountered in subtract

-c:18: RuntimeWarning: invalid value encountered in subtract

|| 3500 || 0 || 0.002 || 0.006 || 0.006 || 0.035 || 0.006 || 0.009 || 0.009 || 0.045 || 0.004 || 0.005 || 0.005 || 0.044 ||
|| 3500 || 2 || 0.000 || 0.004 || 0.003 || 0.016 || 0.004 || 0.006 || 0.005 || 0.036 || 0.003 || 0.005 || 0.004 || 0.029 ||
|| 3500 || 4 || -0.003 || 0.006 || 0.002 || 0.036 || 0.002 || 0.005 || 0.003 || 0.032 || 0.005 || 0.009 || 0.003 || 0.043 ||
|| 4500 || 0 || 0.002 || 0.004 || 0.003 || 0.077 || 0.007 || 0.007 || 0.006 || 0.042 || 0.005 || 0.004 || 0.004 || 0.091 ||
|| 4500 || 2 || 0.002 || 0.003 || 0.002 || 0.023 || 0.003 || 0.004 || 0.004 || 0.027 || 0.001 || 0.002 || 0.002 || 0.025 ||
|| 4500 || 4 || 0.001 || 0.003 || 0.002 || 0.016 || 0.002 || 0.004 || 0.003 || 0.023 || 0.001 || 0.002 || 0.002 || 0.025 ||

Analyze the spectra with the ASSeT library



In [12]:

    
fparam_asset= []
fparam_moog= []
cnt= 0
errs= numpy.ones(8575)*0.01
errs[True-cmpIndx]= 10. #remove areas around H lines
for teff in teffs:
    for logg in loggs:
        fparam_asset.append(apogee.modelspec.ferre.fit(assetspecs[cnt],errs))
        fparam_moog.append(apogee.modelspec.ferre.fit(moogspecs[cnt][0],errs))
        cnt+= 1



In [13]:

    
fparam_turbo= []
cnt= 0
errs= numpy.ones(8575)*0.01
errs[True-cmpIndx]= 10. #remove areas around H lines
for teff in teffs:
    for logg in loggs:
        fparam_turbo.append(apogee.modelspec.ferre.fit(turbospecs[cnt][0],errs))
        cnt+= 1



In [14]:

    
print "||= '''T,,eff,,''' =||= '''logg g''' =||||||||||||= '''ASSeT - Turbospec''' =||||||||||||= '''ASSeT - MOOG''' =||||||||||||= '''Turbospec - MOOG''' =||"
print '|| || || dT,,eff,, || dlog(g) || dmetals || d[C/Fe] || d[N/Fe] || d[a/Fe] || dT,,eff,, || dlog(g) || dmetals || d[C/Fe] || d[N/Fe] || d[a/Fe] || dT,,eff,, || dlog(g) || dmetals || d[C/Fe] || d[N/Fe] || d[a/Fe] ||' 
cnt= 0
for teff in teffs:
    for logg in loggs:
        thisline= '|| {teff} || {logg} || '.format(teff=teff,logg=logg)
        for fp in [fparam_asset[cnt]-fparam_turbo[cnt],fparam_asset[cnt]-fparam_moog[cnt],fparam_turbo[cnt]-fparam_moog[cnt]]:
            for indx in [0,1,3,4,5,6]:
                if indx == 0:
                    thisline+= '{fp:.0f} ||'.format(fp=fp[0,indx])
                else:
                    thisline+= '{fp:.2f} ||'.format(fp=fp[0,indx])
        print thisline
        cnt+= 1









    



||= '''T,,eff,,''' =||= '''logg g''' =||||||||||||= '''ASSeT - Turbospec''' =||||||||||||= '''ASSeT - MOOG''' =||||||||||||= '''Turbospec - MOOG''' =||
|| || || dT,,eff,, || dlog(g) || dmetals || d[C/Fe] || d[N/Fe] || d[a/Fe] || dT,,eff,, || dlog(g) || dmetals || d[C/Fe] || d[N/Fe] || d[a/Fe] || dT,,eff,, || dlog(g) || dmetals || d[C/Fe] || d[N/Fe] || d[a/Fe] ||
|| 3500 || 0 || 7 ||-0.00 ||-0.06 ||0.01 ||0.06 ||0.02 ||12 ||0.00 ||-0.10 ||-0.01 ||-0.02 ||-0.03 ||5 ||0.00 ||-0.05 ||-0.02 ||-0.08 ||-0.06 ||
|| 3500 || 2 || -10 ||-0.01 ||-0.04 ||0.01 ||0.04 ||0.01 ||5 ||0.13 ||-0.04 ||0.01 ||-0.04 ||-0.03 ||14 ||0.13 ||0.00 ||0.00 ||-0.08 ||-0.03 ||
|| 3500 || 4 || -25 ||-0.08 ||-0.00 ||-0.04 ||-0.68 ||-0.02 ||3 ||0.03 ||-0.02 ||-0.01 ||-0.05 ||-0.02 ||29 ||0.11 ||-0.02 ||0.03 ||0.62 ||0.00 ||
|| 4500 || 0 || -9 ||-0.02 ||-0.06 ||0.00 ||0.04 ||0.02 ||-8 ||0.01 ||-0.11 ||0.01 ||0.00 ||-0.01 ||1 ||0.03 ||-0.05 ||0.01 ||-0.04 ||-0.03 ||
|| 4500 || 2 || 9 ||-0.01 ||-0.04 ||0.02 ||0.03 ||0.02 ||5 ||0.04 ||-0.06 ||0.01 ||-0.01 ||-0.01 ||-4 ||0.06 ||-0.02 ||-0.01 ||-0.04 ||-0.03 ||
|| 4500 || 4 || 5 ||-0.03 ||-0.04 ||0.00 ||0.04 ||0.01 ||9 ||-0.02 ||-0.07 ||-0.01 ||0.04 ||0.01 ||4 ||0.00 ||-0.03 ||-0.01 ||0.00 ||-0.01 ||

Comparing individual elements derived from different syntheses



In [15]:

    
felem_asset= []
felem_moog= []
cnt= 0
errs= numpy.ones(8575)*0.01
errs[True-cmpIndx]= 10. #remove areas around H lines
for teff in teffs:
    for logg in loggs:
        felem_asset.append(apogee.modelspec.ferre.elemfitall(assetspecs[cnt],errs,fparam=fparam_asset[cnt]))
        felem_moog.append(apogee.modelspec.ferre.elemfitall(moogspecs[cnt][0],errs,fparam=fparam_moog[cnt]))
        cnt+= 1



In [16]:

    
felem_turbo= []
cnt= 0
errs= numpy.ones(8575)*0.01
errs[True-cmpIndx]= 10. #remove areas around H lines
for teff in teffs:
    for logg in loggs:
        felem_turbo.append(apogee.modelspec.ferre.elemfitall(turbospecs[cnt][0],errs,fparam=fparam_turbo[cnt]))
        cnt+= 1



In [17]:

    
figsize(24,24)
cnt= 0
icnt= 0
for teff in teffs:
    for logg in loggs:
        for dxh in [dict(zip(felem_asset[cnt].keys(),[felem_asset[cnt][k]-felem_turbo[cnt][k] for k in felem_asset[cnt].keys()])),
                   dict(zip(felem_asset[cnt].keys(),[felem_asset[cnt][k]-felem_moog[cnt][k] for k in felem_asset[cnt].keys()])),
                   dict(zip(felem_asset[cnt].keys(),[felem_turbo[cnt][k]-felem_moog[cnt][k] for k in felem_asset[cnt].keys()]))]:
            subplot(6,3,icnt+1)
            splot.elements(dxh,'ko',gcf=True,yrange=[-0.2,0.2])
            if icnt == 0:
                bovy_plot.bovy_text(r'$\mathrm{ASSeT-Turbospec}$',title=True,size=18.)
            elif icnt == 1:
                bovy_plot.bovy_text(r'$\mathrm{ASSeT-MOOG}$',title=True,size=18.)
            elif icnt == 2:
                bovy_plot.bovy_text(r'$\mathrm{Turbospec-MOOG}$',title=True,size=18.)
            if icnt % 3 == 0:
                bovy_plot.bovy_text(r'$T_{{\mathrm{{eff}}}} = {teff}, \log g = {logg}$'.format(teff=teff,logg=logg),
                                    top_left=True,size=18.)
            icnt+= 1
        cnt+= 1
tight_layout()
bovy_plot.bovy_end_print('/Users/bovy/Desktop/indivElement_turboMoogAsset.png',dpi=100)

Plot the spectra in great detail



In [18]:

    
from matplotlib.backends.backend_pdf import PdfPages
# We'll chop each detector into four
startindxs= [188,988,1755,2600,3605,4150,4827,5500,6382,6830,7343,7900]
endindxs= [988,1755,2600,3250,4150,4827,5500,6100,6830,7343,7900,8325]
with PdfPages('/Users/bovy/Desktop/indivElement_turboMoogAsset_spec.pdf') as pdf:
    cnt= 0
    for teff in teffs:
        for logg in loggs:
            for ii in range(len(startindxs)):
                splot.waveregions(assetspecs[cnt],startindxs=[startindxs[ii]],endindxs=[endindxs[ii]],labelLines=True,
                                  labelTeff=teff,labellogg=logg,
                                  yrange=[-0.1,1.2],color='b')
                splot.waveregions(turbospecs[cnt][0],startindxs=[startindxs[ii]],endindxs=[endindxs[ii]],overplot=True,
                                  color='g')
                splot.waveregions(moogspecs[cnt][0],startindxs=[startindxs[ii]],endindxs=[endindxs[ii]],overplot=True,
                                  color='gold')
                splot.waveregions(assetspecs[cnt]-turbospecs[cnt][0]+0.2,startindxs=[startindxs[ii]],endindxs=[endindxs[ii]],
                                  overplot=True,color='b',cleanZero=False)              
                splot.waveregions(assetspecs[cnt]-moogspecs[cnt][0]+0.1,startindxs=[startindxs[ii]],endindxs=[endindxs[ii]],
                                  overplot=True,color='g',cleanZero=False)              
                splot.waveregions(turbospecs[cnt][0]-moogspecs[cnt][0],startindxs=[startindxs[ii]],endindxs=[endindxs[ii]],
                                  overplot=True,color='gold',cleanZero=False)   
                if ii == 0:
                    lsize= 10.
                    bovy_plot.bovy_text(137.,0.925,r'$\mathrm{ASSeT}$',color='b',size=lsize)
                    bovy_plot.bovy_text(137.,0.8,r'$\mathrm{Turbospec}$',color='g',size=lsize)
                    bovy_plot.bovy_text(137.,0.675,r'$\mathrm{MOOG}$',color='gold',size=lsize)
                    bovy_plot.bovy_text(137.,0.2,r'$\mathrm{ASSeT-Turbospec}$',color='b',size=lsize)
                    bovy_plot.bovy_text(137.,0.075,r'$\mathrm{ASSeT-MOOG}$',color='g',size=lsize)
                    bovy_plot.bovy_text(137.,-0.025,r'$\mathrm{Turbospec-MOOG}$',color='gold',size=lsize)
                pdf.savefig()
                close()
            cnt+= 1









    



/Users/bovy/Repos/apogee/apogee/spec/plot.py:176: RuntimeWarning: invalid value encountered in less_equal
  if cleanZero:

/Library/Python/2.7/site-packages/numpy/lib/nanfunctions.py:220: RuntimeWarning: All-NaN axis encountered
  warnings.warn("All-NaN axis encountered", RuntimeWarning)

-c:21: RuntimeWarning: invalid value encountered in subtract



In [18]: