In [1]:
# NB: the recommended way to run isoclassify is through the command-line interface (CLI)
# using examples such as given in examples/examples.csv (see README)
# below is an example to run isoclassify from ipython, which requires some hacking by
# pre-defining functions but allows direct interaction with posteriors.
In [1]:
# required packages
%matplotlib inline
%reload_ext autoreload
%autoreload 2
import os
import copy
import glob
import h5py,pdb
import numpy as np
from matplotlib import pylab as plt
import pandas as pd
import astropy.units as units
from astropy.coordinates import SkyCoord
from dustmaps.bayestar import BayestarWebQuery
import mwdust
from isoclassify.direct import classify as classify_direct
from isoclassify.grid import classify as classify_grid
from isoclassify import DATADIR
In [2]:
# load models
fn = os.path.join(DATADIR,'mesa.h5')
file = h5py.File(fn,'r+', driver='core', backing_store=False)
model = {'age':np.array(file['age']),\
'mass':np.array(file['mass']),\
'feh':np.array(file['feh']),\
'teff':np.array(file['teff']),\
'logg':np.array(file['logg']),\
'rad':np.array(file['rad']),\
'lum':np.array(file['rad']),\
'rho':np.array(file['rho']),\
'dage':np.array(file['dage']),\
'dmass':np.array(file['dmass']),\
'dfeh':np.array(file['dfeh']),\
'eep':np.array(file['eep']),\
'bmag':np.array(file['bmag']),\
'vmag':np.array(file['vmag']),\
'btmag':np.array(file['btmag']),\
'vtmag':np.array(file['vtmag']),\
'gmag':np.array(file['gmag']),\
'rmag':np.array(file['rmag']),\
'imag':np.array(file['imag']),\
'zmag':np.array(file['zmag']),\
'jmag':np.array(file['jmag']),\
'hmag':np.array(file['hmag']),\
'kmag':np.array(file['kmag']),\
'd51mag':np.array(file['d51mag']),\
'gamag':np.array(file['gamag']),\
'fdnu':np.array(file['fdnu']),\
'avs':np.zeros(len(np.array(file['gamag']))),\
'dis':np.zeros(len(np.array(file['gamag'])))}
#ebf.read(os.path.join(DATADIR,'mesa.ebf'))
# prelims to manipulate some model variables (to be automated soon ...)
#pdb.set_trace()
model['rho'] = np.log10(model['rho'])
model['lum'] = model['rad']**2*(model['teff']/5777.)**4
# 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 [3]:
# define class that contains observables
x = classify_grid.obsdata()
In [4]:
# add [Teff, logg, FeH] and [sigma_Teff, sigma_logg, sigma_FeH]
x.addspec([5777.,4.44,0.0],[60.,0.07,0.04])
In [5]:
# perform classification based on those inputs
paras = classify_grid.classify(input=x, model=model, dustmodel=0,plot=1)
In [6]:
# print results for radius
paras.rad,paras.radep,paras.radem
Out[6]:
In [7]:
# plot the radius posterior
plt.plot(paras.radpx,paras.radpy)
Out[7]:
In [8]:
# now let's add some more observables, e.g. a parallax
x.addplx(0.07,0.007)
In [9]:
# and some JHK photometry
x.addjhk([4.38,4.04,4.00],[0.03,0.03,0.03])
In [15]:
# using photometry requires some treatment of reddening and extinction. the following functions
# required for this
def query_dustmodel_coords(ra,dec):
reddenMap = BayestarWebQuery(version='bayestar2017')
sightLines = SkyCoord(ra*units.deg,dec*units.deg,frame='icrs')
reddenContainer = reddenMap(sightLines,mode='best')
del reddenMap # To clear reddenMap from memory
distanceSamples = np.array([0.06309573,0.07943284,0.1,0.12589255,0.15848933,0.19952627,0.25118864,0.31622776,0.3981072,0.50118726,0.6309574,0.7943282 ,1.,1.2589258,1.5848933,1.9952621,2.511887,3.1622777,3.981073,5.011873,6.3095727,7.943284,10.,12.589258,15.848933,19.952621,25.11887,31.622776,39.81073,50.11873,63.095726])*1000. # In pc, from bayestar2017 map distance samples
dustModelDF = pd.DataFrame({'ra': [ra], 'dec': [dec]})
for index in range(len(reddenContainer)):
dustModelDF['av_'+str(round(distanceSamples[index],6))] = reddenContainer[index]
return dustModelDF
def query_dustmodel_coords_allsky(ra,dec):
reddenMap = mwdust.Combined15()
sightLines = SkyCoord(ra*units.deg,dec*units.deg,frame='galactic')
distanceSamples = np.array([0.06309573,0.07943284,0.1,0.12589255,0.15848933,0.19952627,0.25118864,0.31622776,0.3981072,0.50118726,0.6309574,0.7943282 ,1.,1.2589258,1.5848933,1.9952621,2.511887,3.1622777,3.981073,5.011873,6.3095727,7.943284,10.,12.589258,15.848933,19.952621,25.11887,31.622776,39.81073,50.11873,63.095726])*1000. # In pc, from bayestar2017 map distance samples
reddenContainer=reddenMap(sightLines.l.value,sightLines.b.value,distanceSamples/1000.)
del reddenMap # To clear reddenMap from memory
dustModelDF = pd.DataFrame({'ra': [ra], 'dec': [dec]})
for index in range(len(reddenContainer)):
dustModelDF['av_'+str(round(distanceSamples[index],6))] = reddenContainer[index]
return dustModelDF
def extinction(law):
if (law == 'cardelli'):
out = {
"ab":4.1708789,
"av":3.1071930,
"abt":4.3358221,
"avt":3.2867038,
"ag":3.8281101,
"ar":2.7386468,
"ai":2.1109662,
"az":1.4975613,
"aj":0.89326176,
"ah":0.56273418,
"ak":0.35666104,
"aga":2.4623915
}
if (law == 'schlafly11'):
out = {
"ab":3.626,
"av":2.742,
"abt":4.5309214,
"avt":3.1026801,
"ag":3.303,
"ar":2.285,
"ai":1.698,
"az":1.263,
"aj":0.77510388,
"ah":0.50818384,
"ak":0.33957048,
"aga":1.9139634
}
if (law == 'schlafly16'):
# see http://argonaut.skymaps.info/usage under "Gray Component". this is a lower limit.
grayoffset=0.063
out = {
"ab":3.6060565+grayoffset,
"av":2.9197679+grayoffset,
"abt":3.7204173+grayoffset,
"avt":3.0353634+grayoffset,
"ag":3.384+grayoffset,
"ar":2.483+grayoffset,
"ai":1.838+grayoffset,
"az":1.414+grayoffset,
"aj":0.650+grayoffset,
"ah":0.327+grayoffset,
"ak":0.161+grayoffset,
"aga":2.2203186+grayoffset
}
return out
In [16]:
# if we don't want to use a reddening map, isoclassify fits for Av. However, we need to
# define an extinction law
ext = extinction('cardelli')
In [17]:
# perform classification
paras = classify_grid.classify(input=x, model=model, dustmodel=0,plot=1,ext=ext)
In [18]:
# if we want to use a reddening map, we'll need to add coordinates
x.addcoords(243.9052932171665,-08.3694394)
In [19]:
# and define a dustmodel. this is the "allsky" model by Bovy
dustmodel = query_dustmodel_coords_allsky(x.ra,x.dec)
In [20]:
# redo fit, this time using the dustmodel
paras = classify_grid.classify(input=x, model=model, dustmodel=dustmodel,plot=1,ext=ext)
In [ ]: