jobovy / apogee

Tools for dealing with APOGEE data
BSD 3-Clause "New" or "Revised" License
43 stars 25 forks source link

Turbospectrum cm behaviour. #61

Closed drd13 closed 4 years ago

drd13 commented 4 years ago

Hi,

I'm using turbospectrum and Atlas9 to generate synthetic spectra and I've noticed what might be a bug. The parameter cm does not seem to modify the carbon in the spectra. This can, for example, be seen by using FERRE to retrieve synthetically generated spectra.

Here is a figure showing retrieved cm against cm input parameter. I've also attached the code required to reproduce it below.

image

Here is the same thing for am which shows the behavior I would have expected.

image

I am using the atlas9 synthetic spectra as obtained from http://research.iac.es/proyecto/ATLAS-APOGEE//. I imagine that this could maybe be due to a difference between whether the apogee colaborator atlas grid is used.

Code for reproducing cm figure

import apogee.modelspec.turbospec
from apogee.modelatm import atlas9
from apogee.modelspec import ferre
from apogee.tools import paramIndx
import matplotlib.pyplot as plt
import numpy as np

def get_parameters(metals,am,cm):
    atm= atlas9.Atlas9Atmosphere(teff=4750.,logg=2.5,metals=metals,am=am,cm=cm)
    ferre_synspec= apogee.modelspec.turbospec.synth(modelatm=atm,linelist='201404080919',lsf='all',cont='aspcap',vmacro=6.,isotopes='solar')
    errs = np.ones(np.shape(ferre_synspec[0]))*0.01
    ferre_fit = ferre.fit(ferre_synspec[0],errs,lib='GK',pca=True,sixd=True)
    abunds = ferre.elemfitall(ferre_synspec[0],errs,lib='GK',pca=True,sixd=True)
    return ferre_fit,abunds

fparams_list = []
abundances_list = []
cm_range = np.linspace(-1,1,3)
for cm in cm_range:
    fparams,abundances = get_parameters(metals=0,am=0,cm=cm)
    fparams_list.append(fparams)
    abundances_list.append(abundances)

fparams_array = np.concatenate(fparams_list)

plt.plot([-1,1],[-1,1])
plt.plot(cm_range,fparams_array[:,paramIndx("C")],label="C")
plt.xlabel("cm")
plt.ylabel("C retrieved")
plt.legend()
plt.show()

Code for reproducing am figure


fparams_list = []
abundances_list = []
am_range = np.linspace(-1,1,3)
for am in am_range:
    fparams,abundances = get_parameters(metals=0,am=am,cm=0)
    fparams_list.append(fparams)
    abundances_list.append(abundances)

fparams_array = np.concatenate(fparams_list)

plt.plot([-1,1],[-1,1])
plt.plot(am_range,fparams_array[:,paramIndx("alpha")],label="alpha")
plt.xlabel("am")
plt.ylabel("alpha retrieved")
plt.legend()
plt.show()
jobovy commented 4 years ago

Hi,

I think I understand why this happens: unlike MOOG, Turbospectrum does not use the atmosphere's abundances as the base of its model, so every abundance needs to be specified. While [a/Fe] needs to be separately passed to Turbospectrum, carbon is not, and if carbon was not specified as an abundance to be varied, its value was not passed to the line that takes into account the atmosphere's carbon content. I'm surprised that I made this mistake, because I was aware of this Turbospectrum feature when I wrote it, butI guess I never noticed.

I committed a fix that I think should address the issue, but I haven't tested it. Could you test it and let me know?

drd13 commented 4 years ago

Hi,

There were some issues related to tuples being immutable. After modifying

    if not 6 in [elem[0] for elem in args]:
        args= args+([6,0.])

to

if not 6 in [elem[0] for elem in args]:
        args= args+([6,0.],)

the output is what I would expect. image

I'll let you do the commit. Thank you for maintaining the package and the quick response!

jobovy commented 4 years ago

Oops, messed that up, thanks for checking! (I had even tested the line, but forgot the final comma in editing the actual code...). And thanks for a very well described issue, makes it much easier to fix the code!