In [1]:
from isochrones.dartmouth import DartmouthModelGrid
from isochrones.mist import MISTModelGrid
In [2]:
g_dar = DartmouthModelGrid(['g','r','i','J','H','K'])
print(len(g_dar.df))
g_dar.df.head()
Out[2]:
In [3]:
DartmouthModelGrid.phot_systems
Out[3]:
In [4]:
DartmouthModelGrid.phot_bands
Out[4]:
In [5]:
DartmouthModelGrid.get_band('g')
Out[5]:
In [9]:
g_mist = MISTModelGrid(['G','B','V','J','H','K','W1','W2','W3'])
print(len(g_mist.df))
g_mist.df.head()
Out[9]:
In [12]:
MISTModelGrid.phot_systems
Out[12]:
In [6]:
MISTModelGrid.get_band('G')
Out[6]:
Isochrone objectIn order to usefully interface with these model grids, we have to interpolate. This is done using the (perhaps not greatly named) Isochrone object. There are two flavors of this object. One (the original) is based on scipy.interpolate.LinearNDInterpolater, which does interpolation based on a Delaunay triangulation of the points (essentially storing nearest neighbor information). For the size of the Dartmouth grids, this is marginally feasible, as the triangulation takes about 10 minutes to compute and about ~100Mb to store.
However, the MIST grids are ~10x larger than Dartmouth, and the scaling is bad, so the triangulation method becomes no longer feasible. So instead, I implemented a FastIsochrone version of this interpolation that takes advantage of the fact that two of the three dimensions are regular grids in order to quickly find the "box" surrounding a given point, and then I do inverse-distance nearest neighbor interpolation using this box of surrounding points. NB: it's important to do some version of "surrounding" interpolation like this (rather than straight-up k-nearest-neighbors with a KD tree, for example) because the distance between successive gridded values of feh and age are much larger than the distances between adjacent values of mass on the same track, so a k-nearest-neighbors approach always grabs points on the same metallicity & age track. Also, my implementation of the "find-the-surrounding-box" interpolation for a single point is about ~10x faster than the KDTree.
I do note, however, that I trust the Delaunay interpolation more near-ish the edges, where finding a proper surrounding box may be difficult---in other words, this might lead to less desirable results during, e.g., the fast stages of stellar end-of-life. (In fact, I'm not even sure how much I trust the Delaunay interpolation in locations where evolution is happening quickly.)
Below is a demo of the Isochrone objects for both the Dartmouth and MIST models. The general idea is that physical properties can get called as a function of (mass, age, feh), where age is really $\log_{10} (\rm{age}\,[{\rm Gyr}])$, and then magnitudes take those same three parameters plus distance and AV (V-band extinction).
In [1]:
from isochrones.dartmouth import Dartmouth_Isochrone
from isochrones.mist import MIST_Isochrone
In [5]:
dar = Dartmouth_Isochrone()
dar.bands #default bands
Out[5]:
In [2]:
mist = MIST_Isochrone()
mist.bands
Out[2]:
In [6]:
mass, age, feh = (1.01, 9.71, 0.01)
print(dar.logg(mass, age, feh), mist.logg(mass, age, feh))
print(dar.Teff(mass, age, feh), mist.Teff(mass, age, feh))
In [7]:
mist.mag['g'](mass, age, feh, 200, 0.2), mist.mag['V'](mass, age, feh, 200, 0.2)
Out[7]:
In [14]:
# Compare evaluation time
%timeit dar.logg(mass, age, feh)
%timeit mist.logg(mass, age, feh)
In [15]:
distance, AV = (200, 0.2)
pars = (mass, age, feh, distance, AV)
for ic in (dar, mist):
print('{} bands:'.format(ic.name))
for b in ic.bands:
print(b, ic.mag[b](*pars))
We can also compare two different interpolation schemes using the same grid:
In [16]:
from isochrones.dartmouth import Dartmouth_FastIsochrone
dar2 = Dartmouth_FastIsochrone(dar.bands)
for b in dar.bands:
print(b, dar.mag[b](*pars), dar2.mag[b](*pars))
In [17]:
from isochrones import StarModel
In [19]:
mod = StarModel(mist, Teff=(5700,100), logg=(4.5,0.1), feh=(0.0,0.1))
mod.fit(basename='spec_demo')
In [20]:
%matplotlib inline
mod.corner_physical();
OK, that was easy; now what if only broadband photometry is available?
In [21]:
mags = {b:(mist.mag[b](*pars), 0.02) for b in ['G', 'B','V','J','W1']}
mod2 = StarModel(mist, **mags)
In [22]:
mod2.fit(basename='phot_demo')
In [23]:
mod2.corner_physical();
In [24]:
# And, just to take a look at how the samples compare with the observations
mod2.corner_observed();
We can also include parallax, so we can see how much parallax might help in this particular situation.
In [25]:
mod3 = StarModel(mist, parallax=(5.0, 0.4), **mags)
mod3.fit(basename='plax_demo')
In [26]:
mod3.corner_physical();
In [30]:
# How much does the measurement of, e.g., radius improve when we have the parallax?
print(mod2.samples.mass_0_0.quantile([0.15, 0.85]))
print(mod3.samples.mass_0_0.quantile([0.15, 0.85]))
In [31]:
# What about measuring the AV?
print(mod2.samples.AV_0.quantile([0.15, 0.85]))
print(mod3.samples.AV_0.quantile([0.15, 0.85]))
Bottom line seems to be that parallax helps, but not enormously much over the multiband photometry on its own---though of course this is assuming that the stellar models are the truth!
In [32]:
!cat ../isochrones/tests/star2/star.ini
In [33]:
mod4 = StarModel.from_ini(mist, '../isochrones/tests/star2')
In [34]:
mod4.obs.print_ascii()
In [43]:
mod4.fit(basename='star2_demo')
mod4.corner_physical();
In [45]:
# OK, full disclosure, I totally made up that model. Here you can see it makes no sense...
mod4.corner_observed();
In [40]:
!cat ../isochrones/tests/star3/star.ini
In [41]:
mod5 = StarModel.from_ini(mist, '../isochrones/tests/star3')
mod5.obs.print_ascii()
In [42]:
mod5.fit(basename='star3_demo')
mod5.corner_physical();
In [ ]:
# nu_max, delta_nu, also density
# BASTA grids (and PARSEC)
In [2]:
from isochrones import get_ichrone
mist = get_ichrone('mist')
In [3]:
mass, age, feh, distance, AV = (0.95, 9.61, -0.2, 200, 0.2)
mist.mag['g'](mass, age, feh, distance, AV)
Out[3]:
In [4]:
mist.bands
Out[4]:
In [ ]: