In [1]:
%matplotlib inline
%reload_ext autoreload
%autoreload 2
import numpy as np
import matplotlib.pyplot as plt
from classify_grid import *
import os, ebf
from astropy.io import ascii
import time
#import mwdust
In [2]:
# load MIST models
homedir=os.path.expanduser('~/')
model=ebf.read('/Users/petigura/code/isoclassify/mesa.ebf')
In [3]:
df = pd.DataFrame(model)
In [ ]:
In [9]:
# stefan boltzman law
from astropy import constants as c
from astropy import units as u
%pylab
rstar = np.array(df.rad) * c.R_sun
teff = np.array(df.teff)* u.K
Lbol = 4 * pi * rstar**2 * c.sigma_sb * teff**4
L0 = 3.0128e28 * u.W
df['Mbol'] = np.log10(Lbol/L0) / (-0.4)
df['BCK'] = df.kmag - df.Mbol
df['BCV'] = df.vmag - df.Mbol
df['vmk'] = df.vmag - df.kmag
#cut = df.query('0.99 < mass < 1.01 and 4.4 < age < 4.6 and feh==0.0')
cut = df.query('5680 < teff < 5720 and 3 < age < 6')
cols = 'age mass rad teff logg feh BCK BCV vmk'.split()
cut[cols]
cut.BCK.hist()
Out[9]:
In [79]:
teff0 = 5770
logg0 = 4.44
fe0 = 0.0
dteff = 60
dlogg = 0.10
dfeh = 0.04
#df = df[df.age==4.5]
df['dist'] = np.sqrt( ((df.teff - teff0)/dteff)**2 + ((df.logg - logg0)/dlogg)**2 + ((df.feh - fe0)/dfeh)**2)
ref = df.loc[df.dist.idxmin()]
df['dist'] = np.sqrt( ((df.teff - teff0+dteff)/dteff)**2 + ((df.logg - logg0)/dlogg)**2 + ((df.feh - fe0)/dfeh)**2)
ref_teff = df.loc[df.dist.idxmin()]
df['dist'] = np.sqrt( ((df.teff - teff0)/dteff)**2 + ((df.logg - logg0+ dlogg)/dlogg)**2 + ((df.feh - fe0)/dfeh)**2)
ref_logg = df.loc[df.dist.idxmin()]
df['dist'] = np.sqrt( ((df.teff - teff0)/dteff)**2 + ((df.logg - logg0+ dlogg)/dlogg)**2 + ((df.feh - fe0+dfeh)/dfeh)**2)
ref_feh = df.loc[df.dist.idxmin()]
print ref[cols]
print ref_teff[cols]
print ref_logg[cols]
print ref_feh[cols]
changing temperature by 60~K -> changed v-k by -0.04 mag, which corresponds to a change in the bolometric correction of 0.04 mag, which is consistent with what I see.
In [81]:
0.9907*0.04 - 0.0395*0.04
Out[81]:
In [3]:
# parelims to manipulate some model variables (to be automated soon ...)
model['rho']=np.log10(model['rho'])
# next line turns off Dnu scaling relation corrections
model['fdnu'][:]=1.
model['avs']=np.zeros(len(model['teff']))
model['dis']=np.zeros(len(model['teff']))
In [5]:
# next 2 lines allow to use a reddening model (needs galactic coordinates)
#x.addcoords(338.3683920,-9.0227690)
#dustmodel = mwdust.Combined15()
In [6]:
# initilize class with observables
In [67]:
x=obsdata()
# add any combiantion of observables
# Teff, logg, FeH + uncertainties
x.addspec([5777.,4.44,0.0],[60.,0.07,0.04])
# numax & Dnu + uncertainties
#x.addseismo([1240.,63.5],[70.,1.5])
# 2MASS photometry
x.addjhk([-99,-99,-99],[0,0,0.02])
x.addplx(1,0.001)
# Sloan photometry
#x.addgriz([11.776,11.354,11.238,11.178],[0.02,0.02,0.02,0.02])
paras=classify(input=x,model=model,dustmodel=0)
In [16]:
paras=classify(input=x,model=model,dustmodel=0.)
In [64]:
# add any combiantion of observables
# Teff, logg, FeH + uncertainties
#x.addspec([5109.23,3.49595,0.0380301],[60.,0.10,0.04])
x.addspec([5109.23,3.49595,0.0380301],[60.,0.1,0.04])
# numax & Dnu + uncertainties
#x.addseismo([1240.,63.5],[70.,1.5])
# 2MASS photometry
x.addjhk([-99,-99,10.0],[0,0,0.02])
# Sloan photometry
#x.addgriz([11.776,11.354,11.238,11.178],[0.02,0.02,0.02,0.02])
x.addplx(1./1372.,1./372.*0.03)
# run classification
%pylab inline
paras=classify(input=x,model=model,dustmodel=0.)
gcf().set_tight_layout(True)
In [62]:
paras.
Out[62]:
In [17]:
from scipy.interpolate import interp1d
In [18]:
interp1d?
In [9]:
# print mass median +/- 1 sigma, plot posterior
print paras.feh,paras.fehep,paras.fehem
plt.plot(paras.fehpx,paras.fehpy)
Out[9]:
In [ ]:
# plot teff posterior
plt.plot(paras.teffpx,paras.teffpy)
In [ ]:
In [ ]:
# print age median +/- 1 sigma, plot posterior
print paras.age,paras.ageep,paras.ageem
plt.plot(paras.agepx,paras.agepy)
In [ ]:
# print mass median +/- 1 sigma, plot posterior
print paras.mass,paras.massep,paras.massem
plt.plot(paras.masspx,paras.masspy)
In [ ]:
# print mass median +/- 1 sigma, plot posterior
print paras.feh,paras.fehep,paras.fehem
plt.plot(paras.fehpx,paras.fehpy)
In [ ]:
# delete numax & Dnu constraint
x.addseismo([-99.,-99.],[70.,1.5])
# add parallax with a 3% uncertainty
x.addplx(1./372.,1./372.*0.03)
In [ ]:
# re-run classification
paras=classify(input=x,model=model,dustmodel=0.,doplot=0)
In [ ]:
# print age median +/- 1 sigma, plot posterior
print paras.age,paras.ageep,paras.ageem
plt.plot(paras.agepx,paras.agepy)
In [ ]:
In [73]:
10 / np.log(10) * 60 / 5700
Out[73]:
In [ ]: