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

# # Padova Demo
# 
# Padova provides pythonic access to Padova group's stellar isochrone web app, [CMD](http://stev.oapd.inaf.it/cgi-bin/cmd). This notebook will briefly show you what Padova can do. First, some setup for this notebook:

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
get_ipython().run_line_magic('config', "InlineBackend.figure_format='retina'")
import numpy as np


# In[2]:


get_ipython().run_cell_magic('bash', '', 'pip install -U palettable\n')


# ****
# 
# ## Intro
# 
# Padova provides three classes tailored to the types of isochrone queries you want to make.
# 
# - `IsochroneRequest` will query for a single isochrone
# - `AgeGridRequest` will query for a set of isochrones with a single metallicity for a grid of ages
# - `MetallicityGridRequest` will query for a set of isochrones with a single age, for a grid of metallicities.

# In[3]:


from padova import IsochroneRequest, AgeGridRequest, MetallicityGridRequest


# ## Single isochrone requests
# 
# To obtain a single isochrone for a speific age and metallicity we'll create a `IsochroneRequest` instance and pass the stellar population settings: the age as $\log (A/\mathrm{yr})$, metallicity as a fraction the composition. We also specify the photometric system, (WFC3 filters). Since WFC3 isn't supported by the latest bolometric correction library, we use the older Girardi et al 2008 specification (`'odfnew'`).

# In[4]:


r = IsochroneRequest(log_age=9, z=0.02, photsys='wfc3_wide', photsys_version='odfnew')


# Now that the request has been computed and downloaded from the CMD API, we can inspect the isochrones. First, we can look at the raw result:

# In[5]:


r.data[:2000]


# The Padova package parses this result into Astropy-compatible tables. You can access the isochrone table with `isochone` attribute.

# In[6]:


isoc = r.isochrone
print("Isochrone of Z={0:.2f}, age={0:.3e}".format(isoc.z, isoc.age))
print("Available filters: {0}".format(" ".join(isoc.filter_names)))
print(isoc.meta)
print(isoc)


# We can plot up the isochrones as if they were in a nearby galaxy:

# In[7]:


from astropy.coordinates import Distance
import astropy.units as u
import matplotlib.pyplot as plt
import matplotlib as mpl

fig = plt.figure(figsize=(4, 4), dpi=300)
gs = mpl.gridspec.GridSpec(1, 1)
ax = fig.add_subplot(gs[0])

d = Distance(785 * u.kpc)
ax.plot(isoc['F475W'] - isoc['F814W'], isoc['F814W'] + d.distmod.value)
ax.set_xlim(-1, 4.)
ax.set_ylim(36., 16.)
ax.set_xlabel("F475W-F814W")
ax.set_ylabel("F814W")
fig.show()


# ## Requesting isochrones for a grid of ages
# 
# Rather than repeatedly calling the CMD server to generate isochrones over a grid, we can generate age or metallicity grids in one call. Here we'll request a set of isochrones spanning a grid of ages, for a single metallicity.

# In[8]:


r2 = AgeGridRequest(z=0.02,
                    min_log_age=7., max_log_age=10.1, delta_log_age=0.1,
                    photsys='wfc3_wide', photsys_version='odfnew')


# In[9]:


r2_isochrones = r2.isochrone_set


# In[10]:


from palettable.cubehelix import perceptual_rainbow_16
from astropy.coordinates import Distance
import astropy.units as u
import matplotlib.pyplot as plt
import matplotlib as mpl

fig = plt.figure(figsize=(4, 4), dpi=300)
gs = mpl.gridspec.GridSpec(1, 1)
ax = fig.add_subplot(gs[0])

scalar_map = mpl.cm.ScalarMappable(norm=mpl.colors.Normalize(vmin=7., vmax=10.1),
                                   cmap=perceptual_rainbow_16.mpl_colormap)
scalar_map.set_array(np.array([isoc.age for isoc in r2_isochrones]))

d = Distance(785 * u.kpc)
for isoc in r2_isochrones:
    ax.plot(isoc['F475W'] - isoc['F814W'], isoc['F814W'] + d.distmod.value,
            c=scalar_map.to_rgba(np.log10(isoc.age)))
cb = plt.colorbar(mappable=scalar_map, cax=None, ax=ax)
cb.set_label(r"$\log(A/\mathrm{yr})$")
ax.set_xlim(-1, 4.)
ax.set_ylim(36., 16.)
ax.set_xlabel("F475W-F814W")
ax.set_ylabel("F814W")
fig.show()


# ## Requesting isochrones for a grid of metallicities
# 
# This request behaves similarly to the age grid. Here's a complete demonstration:

# In[11]:


r3 = MetallicityGridRequest(log_age=9., min_z=0.0001, max_z=0.03, delta_z=0.001,
                            photsys='wfc3_wide', photsys_version='odfnew')
r3_isochrones = r3.isochrone_set


# In[12]:


from palettable.cubehelix import perceptual_rainbow_16
from astropy.coordinates import Distance
import astropy.units as u
import matplotlib.pyplot as plt
import matplotlib as mpl

fig = plt.figure(figsize=(4, 4), dpi=300)
gs = mpl.gridspec.GridSpec(1, 1)
ax = fig.add_subplot(gs[0])

scale_z = lambda z: np.log10(z / 0.019)
z_values = np.array([scale_z(isoc.z) for isoc in r3_isochrones])

scalar_map = mpl.cm.ScalarMappable(norm=mpl.colors.Normalize(vmin=z_values.min(), vmax=z_values.max()),
                                   cmap=perceptual_rainbow_16.mpl_colormap)
scalar_map.set_array(z_values)

d = Distance(785 * u.kpc)
for isoc in r3_isochrones:
    ax.plot(isoc['F475W'] - isoc['F814W'], isoc['F814W'] + d.distmod.value,
            c=scalar_map.to_rgba(scale_z(isoc.z)))
cb = plt.colorbar(mappable=scalar_map, cax=None, ax=ax)
cb.set_label(r"$\log(Z/Z_\odot)$")
ax.set_xlim(-1, 4.)
ax.set_ylim(36., 16.)
ax.set_xlabel("F475W-F814W")
ax.set_ylabel("F814W")
fig.show()


# ## Caching
# 
# Padova caches all results from the CMD web application. This means that you can repeatedly test and run your code without hammering the CMD servers, or experiencing any lag. Each request instance has a hash generated from the complete parameter set (including defaults):

# In[13]:


r.settings.__hash__()


# We can verify that the request we just made has been cached:

# In[14]:


from padova.resultcache import PadovaCache
cache = PadovaCache()
r.settings in cache


# Note that by default the cache files are located in `~/.padova_cache`.

# ## Exporting Isochrones to StarFISH
# 
# StarFISH is a Hess diagram simulation and fitting tool by Harris and Zaritsky. Padova can export isochrones into the format expected by StarFISH.

# In[15]:


isoc.export_for_starfish("here")


# In[ ]: