Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Calculate stellar mass and radius using the Torres calibration
import numpy as np
try: # numba will provide a ~2x speedup
from numba import jit
except ImportError: # but we can do without it
jit = lambda fn: fn
@jit
def massTorres(teff, erteff, logg, erlogg, feh, erfeh,
ntrials=10000, corrected=True, add_intrinsic=True):
"""
Calculate stellar mass using the Torres et al. (2010) callibration.
Parameters
----------
teff, erteff : floats
Effective temperature and associated uncertainty.
logg, erlogg : floats
Surface gravity and associated uncertainty.
feh, erfeh : floats
Metallicity [Fe/H] and associated uncertainty.
ntrials : int, optional
Number of Monte Carlo trials for the uncertainty calculation.
corrected : bool, optional
Correct the mass for the offset relative to isochrone-derived masses.
add_intrinsic : bool, optional
Add (quadratically) the intrinsic error of the calibration
(0.027 in log mass).
Returns
-------
meanMass, sigMass : floats
Estimate for the stellar mass and associated uncertainty.
"""
randomteff = teff + erteff * np.random.randn(ntrials)
randomlogg = logg + erlogg * np.random.randn(ntrials)
randomfeh = feh + erfeh * np.random.randn(ntrials)
# Parameters for the Torres calibration
a1 = 1.5689
a2 = 1.3787
a3 = 0.4243
a4 = 1.139
a5 = -0.1425
a6 = 0.01969
a7 = 0.1010
X = np.log10(randomteff) - 4.1
logMass = a1 + a2*X + a3*X**2 + a4*X**3 + a5*randomlogg**2 \
+ a6*randomlogg**3 + a7*randomfeh
meanlogMass = np.mean(logMass)
siglogMass = np.sum((logMass - meanlogMass)**2) / (ntrials - 1)
if add_intrinsic:
siglogMass = np.sqrt(0.027*0.027 + siglogMass)
meanMass = 10**meanlogMass
sigMass = 10**(meanlogMass + siglogMass) - meanMass
if corrected:
# correction comes from Santos+(2013), the SWEET-Cat paper
corrected_meanMass = 0.791 * meanMass**2 - 0.575 * meanMass + 0.701
return corrected_meanMass, sigMass
else:
return meanMass, sigMass
@jit
def radiusTorres(teff, erteff, logg, erlogg, feh, erfeh,
ntrials=10000, add_intrinsic=True):
"""
Calculate stellar radius using the Torres et al. (2010) callibration.
Parameters
----------
teff, erteff : floats
Effective temperature and associated uncertainty.
logg, erlogg : floats
Surface gravity and associated uncertainty.
feh, erfeh : floats
Metallicity [Fe/H] and associated uncertainty.
ntrials : int, optional
Number of Monte Carlo trials for the uncertainty calculation.
add_intrinsic : bool, optional
Add (quadratically) the intrinsic error of the calibration
(0.014 in log radius).
Returns
-------
meanRadius, sigRadius : floats
Estimate for the stellar radius and associated uncertainty.
"""
randomteff = teff + erteff * np.random.randn(ntrials)
randomlogg = logg + erlogg * np.random.randn(ntrials)
randomfeh = feh + erfeh * np.random.randn(ntrials)
# Parameters for the Torres calibration
b1 = 2.4427
b2 = 0.6679
b3 = 0.1771
b4 = 0.705
b5 = -0.21415
b6 = 0.02306
b7 = 0.04173
X = np.log10(randomteff) - 4.1
logRadius = b1 + b2*X + b3*X**2 + b4*X**3 + b5*randomlogg**2 \
+ b6*randomlogg**3 + b7*randomfeh
meanlogRadius = np.mean(logRadius)
siglogRadius = np.sum((logRadius - meanlogRadius)**2) / (ntrials - 1)
if add_intrinsic:
siglogRadius = np.sqrt(0.014*0.014 + siglogRadius)
meanRadius = 10**meanlogRadius
sigRadius = 10**(meanlogRadius + siglogRadius) - meanRadius
return meanRadius, sigRadius
@jason-neal

This comment has been minimized.

Copy link

@jason-neal jason-neal commented Jul 20, 2018

I like the jit lambda trick.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment