MIST_Isochrone.isochrone() inconsistent with MISTModelGrid.df isochrone

[rickyegeland]

The code below produces two very different isochrones, where I expect them to be the same. This is true for the whole range of allowed (age, feh) which does not cause isochrone() to fail (that's another issue).

from isochrones.mist import MISTModelGrid, MIST_Isochrone
import numpy as np
from matplotlib import pyplot as plt

def test(age=9.0, feh=0.0):
    age_Gyr = 10.**age/1e9
    
    grid = MISTModelGrid(['B', 'V'])
    sel = np.isclose(grid.df.feh, feh) & np.isclose(grid.df.log10_isochrone_age_yr, age)
    iso1 = grid.df[sel]
    
    mist = MIST_Isochrone()
    iso2 = mist.isochrone(age, feh=feh)
    
    plt.figure()
    plt.plot(iso1.log_Teff, iso1.log_L, 'r-', label='grid.df')
    plt.plot(np.log10(iso2.Teff), iso2.logL, 'b-', label='mist.isochrone')
    plt.gca().invert_xaxis()
    plt.xlabel('log(Teff)')
    plt.ylabel('log(L)')
    plt.title("log(age)=%0.2f (%0.1f Gyr) [Fe/H]=%0.2f" % (age, age_Gyr, feh))
    plt.legend()
test()

[timothydmorton]

Thanks for raising these issues. I now have some time to actually go back and look into some of this a bit; I hope you haven't given up on it yet.

[msotov]

Hi,

I recently found out I have the same issue when using MIST_Isochrone().isochrone(age, feh). Doing something similar than what was posted previously by @rickyegeland:

from isochrones.mist import MIST_Isochrone, MISTModelGrid
import numpy as np
import matplotlib.pyplot as plt

mist = MIST_Isochrone()

iso1 = mist.isochrone(9.0, feh=0.0)

grid = MISTModelGrid(['B', 'V'])
sel = np.isclose(grid.df['[Fe/H]_init'], 0.0) & np.isclose(grid.df.log10_isochrone_age_yr, 9.0)
iso2 = grid.df[sel]

fig,(ax1, ax2) = plt.subplots(1,2, figsize=(15,6))

ax1.plot(iso2.log_Teff, iso2.log_g, label='grid.df,\nmass_range=%.3f-%.3f' % (min(iso2.initial_mass), max(iso2.initial_mass)))
ax1.plot(np.log10(iso1['Teff']), iso1['logg'], label='mist.isochrone,\nmass_range=%.3f-%.3f' % (min(iso1['M']), max(iso1['M'])))
ax2.plot(iso2.log_Teff, iso2.log_L, label='grid.df,\nmass_range=%.3f-%.3f' % (min(iso2.initial_mass), max(iso2.initial_mass)))
ax2.plot(np.log10(iso1['Teff']), iso1['logL'], label='mist.isochrone,\nmass_range=%.3f-%.3f' % (min(iso1['M']), max(iso1['M'])))
    
ax1.legend()
ax2.legend()
ax1.invert_yaxis()
ax1.invert_xaxis()
ax2.invert_xaxis()

ax1.set_xlabel('logTeff')
ax1.set_ylabel('logg')
ax2.set_xlabel('logTeff')
ax2.set_ylabel('logL')

plt.show()

Both models were initialized with the same values, but the results are very different. It looks like mist.isochrone() halts at the end of the main sequence.

I've been trying to use the code to fit giant stars, but the age always seems to be off, even when I add bounds to the age prior in StarModel to that it discards all the pre-main sequence stage. I wonder if this bug with mist.isochrone() has something to do with it? Maybe there's an issue with the interpolation within the FastIsochrone() class?

I also compared the evolutionary tracks with the ones you can download from the web interpolator in the MIST homepage. For Mstar=1.17 and [Fe/H]=0.17, I get

mist.evtrack gives me a very incorrect result, even the ages range for the evolutionary track from isochrones is weird, because the evolution of the star stops at the end of the main sequence, but moves no further. I think this problem with the interpolation (?) is what is causing the incorrect age determination of masses and ages for giant stars. Unfortunately I haven't been able to find the bug in the source code. I'd really appreciate it if you could look at this.

Thanks you!

[timothydmorton]

Hi! Thanks for writing. As this issue points out, the currently released version of isochrones (interpolating in mass/age/feh) does not do a good job interpolating beyond the MS.

If you have a bit of patience to wait it out, I have been working on major updates/improvements to this package that I hope to release sometime in the next few months that will completely fix this issue.

And if you'd rather not wait for the release and have an adventurous spirit, please feel free to try out the 'eep' branch of the code. This does interpolation in EEP/age/feh space, and thus fits for EEP instead of mass, which allows the interpolation to work well across the whole HR diagram.

This having been said, the 'eep' branch isn't actually where the next version will come from; rather, that's on the 'bolo' branch (which also corrects the implementation of extinction, and uses both isochrone grids and evolution track grids). However, that branch is really still in flux, but if you're interested in that let me know and I can give you some further guidance.

[vedantchandra]

@timothydmorton What's the fastest/most convenient way to interpolate from mass-age-metallicity directly to photometric magnitudes? If we are only considering MS stars, can we do away with the EEP step in return for a speed gain (and some small error)?

[timothydmorton]

Currently no; .generate() is right now the only way. But it appears that there are enough users interested in this that I should put some time into improving it!