#!/usr/bin/env python
# coding: utf-8

# # Overview of `isochrones`

# ## The `ModelGrid` object

# A `ModelGrid` gives an easy interface to retrieve a grid of stellar models, along with synthetic photometry in any available bands, as desired.  I have implemented the [dartmouth](http://stellar.dartmouth.edu/models/) and [MIST](http://waps.cfa.harvard.edu/MIST/) models, the implementation of which should provide guidance on doing other models.

# 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()


# In[3]:


DartmouthModelGrid.phot_systems


# In[4]:


DartmouthModelGrid.phot_bands


# In[5]:


DartmouthModelGrid.get_band('g')


# In[9]:


g_mist = MISTModelGrid(['G','B','V','J','H','K','W1','W2','W3'])
print(len(g_mist.df))
g_mist.df.head()


# In[12]:


MISTModelGrid.phot_systems


# In[6]:


MISTModelGrid.get_band('G')


# ## The `Isochrone` object
# 
# In 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


# In[2]:


mist = MIST_Isochrone()
mist.bands


# 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)


# In[14]:


# Compare evaluation time
get_ipython().run_line_magic('timeit', 'dar.logg(mass, age, feh)')
get_ipython().run_line_magic('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))


# ## The `StarModel` object
# 
# These model grids and interpolation schemes are great, but their purpose for existing is mostly to serve the `StarModel` object, which allows for fitting stellar properties.

# 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]:


get_ipython().run_line_magic('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! 

# ## More about `StarModel`
# 
# We can do much more with the `StarModel` object.  In particular, we can model binary stars.  Let's take a look at some examples.

# In[32]:


get_ipython().system('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]:


get_ipython().system('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)


# In[4]:


mist.bands


# In[ ]: