notebook.community

Edit and run



In [1]:

    
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

from isochrones.dartmouth import Dartmouth_Isochrone
from isochrones import StarModel



In [2]:

    
dar = Dartmouth_Isochrone()



In [3]:

    
mass,age,feh = (1,9.7,0.0)
dar(mass, age, feh, return_df=False) #return a dictionary instead of DataFrame









    Out[3]:





{'Teff': 5787.9238187659339,
 'age': 9.7,
 'logL': array(0.016997207862528108),
 'logg': array(4.424630134394198),
 'mag': {'B': array(5.446314300786797),
  'D51': array(4.930972308028574),
  'H': array(3.288792194844888),
  'I': array(4.052893551792754),
  'J': array(3.604890258678507),
  'K': array(3.2521177420306944),
  'Kepler': array(4.517391615626373),
  'R': array(4.384175974816311),
  'U': array(5.698555620262865),
  'V': array(4.754756162064822),
  'W1': array(3.237289453916152),
  'W2': array(3.265105094726214),
  'W3': array(3.2256690889791524),
  'g': array(5.063454478475286),
  'i': array(4.459568004606947),
  'r': array(4.567647939278613),
  'z': array(4.4508446461643665)},
 'mass': array(1.0),
 'radius': 1.0160123968403469}



In [4]:

    
ages = np.linspace(9,9.8,200)
radii1 = dar.radius(1,ages,0.0)
radii2 = dar.radius(1,ages,0.2)
radii3 = dar.radius(1,ages,-0.2)
plt.plot(ages,radii1,label='feh=0.0')
plt.plot(ages,radii2,label='feh=0.2')
plt.plot(ages,radii3,label='feh=-0.2')
plt.xlabel('log(age)')
plt.ylabel('radius')
plt.legend(loc='upper left');



In [5]:

    
ev = dar.evtrack(1.0,dage=0.01)
plt.plot((ev['Teff']),ev['logL'],'.')
plt.xlabel('Teff')
plt.ylabel('logL');



In [6]:

    
mod = StarModel(dar,Teff=(5770,80),logg=(4.44,0.10),feh=(0.0,0.1))

Find best point estimate for model:



In [7]:

    
mod.maxlike()









    Out[7]:





array([  9.90154610e-01,   9.70841639e+00,  -5.77333630e-03])

Fit for uncertainties using MCMC:



In [8]:

    
mod.fit_mcmc()









    Out[8]:





<emcee.ensemble.EnsembleSampler at 0x7fba7d685cd0>

Let's take a look at some parameters:



In [9]:

    
mod.plot_samples('radius')
mod.plot_samples('age')









    



---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-9-a677741e10ea> in <module>()
----> 1 mod.plot_samples('radius')
      2 mod.plot_samples('age')

/u/tdm/anaconda/lib/python2.7/site-packages/isochrones-0.8.2-py2.7.egg/isochrones/starmodel.pyc in plot_samples(self, prop, fig, label, histtype, bins, lw, **kwargs)
    728         """
    729         setfig(fig)
--> 730         samples,stats = self.prop_samples(prop)
    731         fig = plt.hist(samples,bins=bins,normed=True,
    732                  histtype=histtype,lw=lw,**kwargs)

/u/tdm/anaconda/lib/python2.7/site-packages/isochrones-0.8.2-py2.7.egg/isochrones/starmodel.pyc in prop_samples(self, prop, return_values, conf)
    703             lo = med - sorted[lo_ind]
    704             hi = sorted[hi_ind] - med
--> 705             return samples.values, (med,lo,hi)
    706         else:
    707             return samples.values

AttributeError: 'numpy.ndarray' object has no attribute 'values'





    





<matplotlib.figure.Figure at 0x7fba71bbaa90>

We can also look at triangle plots (if you have the triangle package installed):



In [10]:

    
mod.triangle_plots()









    Out[10]:





(<matplotlib.figure.Figure at 0x7fba78987590>,
 <matplotlib.figure.Figure at 0x7fba78987e50>)

Kepler-93 (KOI-69) has published values: $M = 0.91 \pm 0.06$, $R = 0.92 \pm 0.02$, based on an asteroseismic $\log g$:



In [3]:

    
mags = {'H': (8.4459999999999997,0.08),
 'J': (8.7710000000000008,0.08),
 'K': (8.3699999999999992,0.08),
 'g':(10.33,0.08),
 'r':(9.83,0.08)}
specprops = dict(Teff=(5626,64),logg=(4.47,0.07),feh=(-.16,0.1))



In [4]:

    
smod_phot = StarModel(dar,**mags)
smod_phot.fit_mcmc()









    Out[4]:





<emcee.ensemble.EnsembleSampler at 0x109b31e90>



In [7]:

    
smod_phot.prop_triangle();



In [5]:

    
smod_phot.triangle(params=['mass','radius'], truths=[0.91,0.92]);



In [14]:

    
smod_spec = StarModel(dar,**specprops)
smod_spec.fit_mcmc()









    Out[14]:





<emcee.ensemble.EnsembleSampler at 0x7fba71bbab50>



In [15]:

    
smod_spec.triangle(params=['mass','radius'], truths=[0.91,0.92]);



In [16]:

    
allprops = {'H': (8.4459999999999997,0.08),
 'J': (8.7710000000000008,0.08),
 'K': (8.3699999999999992,0.08),
 'g':(10.33,0.08),
 'r':(9.83,0.08), 
 'Teff':(5626,64),
  'logg':(4.47,0.07),
  'feh':(-0.16,0.1)}

smod_all = StarModel(dar,**allprops)
smod_all.fit_mcmc()









    Out[16]:





<emcee.ensemble.EnsembleSampler at 0x7fba6e14db50>



In [17]:

    
smod_all.triangle(params=['mass','radius'], truths=[0.91,0.92]);



In [ ]: