Source code for

import numpy as np
from collections import defaultdict
from scipy.interpolate import RegularGridInterpolator
import pkg_resources

from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy import units as u

import beast.plotting.plot_compare_spec_type as plot_match

[docs] def compare_spec_type( phot_cat_file, beast_stats_file, spec_ra, spec_dec, spec_type, spec_subtype, lumin_class, match_radius=1.0, bright_filter="F475W", output_filebase=None, plot=False ): """ BEAST verification: compare spectrally-typed stars to BEAST fits Parameters ---------- phot_cat_file : string name of the file with the photometry beast_stats_file : string name of the file with the BEAST stats spec_ra, spec_dec : array of floats coordinates for the spectrally typed stars (in degrees) spec_type : list of strings spectral type for each star (O/B/A/F/G/K/M) spec_subtype : list of floats spectral subtype for each star (e.g., 7 for a B7 star, or 3.5 for a G3.5 star) lumin_class : list of strings luminosity class for each star (I/II/III/IV/V) match_radius : float (default=1) radius (arcsec) for searching for the BEAST star matching the spectrally typed star bright_filter : string (default='F475W') the star within match_radius that is brightest in this filter will be chosen as the match output_filebase : string (default=None) Prefix for saving the file of match info ("_spectype_match.fits" will be appended). If None, will return the table rather than saving it. plot : boolean (default=False) Plot stars found to be a match. Returns ------- spec_match : dict if output_filebase is None, a dictionary with match info is returned. """ # read in the photometry catalog beast_phot = ra_col = [x for x in beast_phot.colnames if "RA" in x.upper()][0] dec_col = [x for x in beast_phot.colnames if "DEC" in x.upper()][0] beast_phot_catalog = SkyCoord( ra=beast_phot[ra_col] *, dec=beast_phot[dec_col] * ) # read in BEAST results beast_stats =, hdu=1) beast_stats_catalog = SkyCoord( ra=beast_stats["RA"] *, dec=beast_stats["DEC"] * ) # setup the effective temperature table teff_function = setup_teff_table() # setup the surface gravity table logg_function = setup_logg_table() # dictionary to hold the matching results spec_match = defaultdict(list) for i in range(len(spec_ra)): # save star info spec_match["spec_ra"].append(spec_ra[i]) spec_match["spec_dec"].append(spec_dec[i]) spec_match["spec_type"].append( "{0} {1} {2}".format(spec_type[i], spec_subtype[i], lumin_class[i]) ) # Finding distances to every star in the catalog will take a while. # Therefore, first check to see if the spectrally typed star is even in # the field, by finding the distance for the closest star. If it's a # small enough number, then proceed with calculating all separations. c = SkyCoord(spec_ra[i], spec_dec[i], unit="deg") _, sep_test, _ = c.match_to_catalog_sky(beast_stats_catalog) # print(sep_test[0].arcmin) # if the star is close enough, continue matching if sep_test[0].arcsec <= match_radius: # calculate separations for full catalog sep = c.separation(beast_stats_catalog) # find sources within user-defined radius small_sep_ind = np.where(sep.arcsec < match_radius)[0] # get photometry for those sources phot_list = np.zeros(len(small_sep_ind)) phot_ind_list = np.zeros(len(small_sep_ind)).astype(int) for j in range(len(small_sep_ind)): phot_ind, _, _ = beast_stats_catalog[ small_sep_ind[j] ].match_to_catalog_sky(beast_phot_catalog) phot_ind_list[j] = phot_ind phot_ref_col = bright_filter + "_RATE" if phot_ref_col.lower() in beast_phot.colnames: phot_list[j] = beast_phot[phot_ref_col.lower()][phot_ind] elif phot_ref_col.upper() in beast_phot.colnames: phot_list[j] = beast_phot[phot_ref_col.upper()][phot_ind] else: raise ValueError("{} not in catalog file".format(bright_filter)) # find brightest match best_ind = small_sep_ind[phot_list == np.max(phot_list)][0] best_phot_ind = phot_ind_list[phot_list == np.max(phot_list)][0] # grab physical parameters for that source # - temperature teff_p16 = 10 ** beast_stats["logT_p16"][best_ind] teff_p50 = 10 ** beast_stats["logT_p50"][best_ind] teff_p84 = 10 ** beast_stats["logT_p84"][best_ind] # - log(g) logg_p16 = beast_stats["logg_p16"][best_ind] logg_p50 = beast_stats["logg_p50"][best_ind] logg_p84 = beast_stats["logg_p84"][best_ind] # grab physical parameters for the spectral type star_teff = lookup_param( teff_function, spec_type[i], spec_subtype[i], lumin_class[i] ) star_logg = lookup_param( logg_function, spec_type[i], spec_subtype[i], lumin_class[i] ) # calculate the number of sigmas away from the "true" value # - temperature if star_teff > teff_p50: teff_sigma = (star_teff - teff_p50) / (teff_p84 - teff_p50) else: teff_sigma = -(teff_p50 - star_teff) / (teff_p50 - teff_p16) # - log(g) if star_logg > logg_p50: logg_sigma = (star_logg - logg_p50) / (logg_p84 - logg_p50) else: logg_sigma = -(logg_p50 - star_logg) / (logg_p50 - logg_p16) # save the results spec_match["spec_teff"].append(star_teff) spec_match["spec_logg"].append(star_logg) spec_match["phot_cat_ind"].append(best_phot_ind) spec_match["stats_cat_ind"].append(best_ind) spec_match["beast_teff_p50"].append(teff_p50) spec_match["beast_teff_p16"].append(teff_p16) spec_match["beast_teff_p84"].append(teff_p84) spec_match["beast_logg_p50"].append(logg_p50) spec_match["beast_logg_p16"].append(logg_p16) spec_match["beast_logg_p84"].append(logg_p84) spec_match["teff_sigma"].append(teff_sigma) spec_match["logg_sigma"].append(logg_sigma) # if there's not a close enough match, put in placeholder values else: spec_match["spec_teff"].append(np.nan) spec_match["spec_logg"].append(np.nan) spec_match["phot_cat_ind"].append(np.nan) spec_match["stats_cat_ind"].append(np.nan) spec_match["beast_teff_p50"].append(np.nan) spec_match["beast_teff_p16"].append(np.nan) spec_match["beast_teff_p84"].append(np.nan) spec_match["beast_logg_p50"].append(np.nan) spec_match["beast_logg_p16"].append(np.nan) spec_match["beast_logg_p84"].append(np.nan) spec_match["teff_sigma"].append(np.nan) spec_match["logg_sigma"].append(np.nan) # write out the table if output_filebase is not None: spec_match_table = output_filebase + "_spectype_match.fits" Table(spec_match).write(spec_match_table, overwrite=True) if plot: plot_match.plot_compare_spec_type(spec_match_table, savefig='png') else: return dict(spec_match)
[docs] def lookup_param(input_function, spec_type, spec_subtype, lumin_class): """ For given spectral classification, evaluate a function to find the physical parameter (e.g., T_eff, log(g)) Note: the `.item()` at the end is because it would otherwise return, e.g., `array(99.9)` instead of `99.9`. Parameters ---------- input_function : RegularGridInterpolator object function to interpolate across a table spec_type : string spectral type for the star (O/B/A/F/G/K/M) spec_subtype : float spectral subtype for the star (e.g., 7 for a B7 star, or 3.5 for a G3.5 star) lumin_class : string luminosity class for the star (I/II/III/IV/V) Returns ------- function_value : float the physical parameter for the star """ return input_function( (_lumin_class_to_number(lumin_class), _spectype_to_id(spec_type, spec_subtype)) ).item()
[docs] def setup_teff_table(): """ Read in and interpolate the T_eff table Ref: Table B.3 and B.4 in Stellar Spectral Classification Returns ------- teff_interp_function : RegularGridInterpolator object function to interpolate across the T_eff table """ # read in the table data_path = pkg_resources.resource_filename("beast", "tools/data/") teff_table = + "effective_temperature.txt", format="ascii") # make each row a number rather than a letter+number teff_table["row_id"] = np.zeros(len(teff_table)) for i in range(len(teff_table)): teff_table["row_id"][i] = _spectype_to_id( teff_table["spec_type"][i], teff_table["spec_subtype"][i] ) # interpolate the columns with missing data nans = np.isnan(teff_table["I"]) teff_table["I"][nans] = np.interp( teff_table["row_id"][nans], teff_table["row_id"][~nans], teff_table["I"][~nans] ) nans = np.isnan(teff_table["III"]) teff_table["III"][nans] = np.interp( teff_table["row_id"][nans], teff_table["row_id"][~nans], teff_table["III"][~nans], ) # now make a general interpolator teff_interp_function = RegularGridInterpolator( (np.array([1, 3, 5]), teff_table["row_id"]), np.array([teff_table["I"], teff_table["III"], teff_table["V"]]), ) return teff_interp_function
[docs] def setup_logg_table(): """ Read in and interpolate the log(g) table Ref: Table 15.8 in Allen's Astrophysical Quantities Returns ------- logg_interp_function : RegularGridInterpolator object function to interpolate across the log(g) table """ # read in the table data_path = pkg_resources.resource_filename("beast", "tools/data/") logg_table = + "surface_gravity.txt", format="ascii") # make each row a number rather than a letter+number logg_table["row_id"] = np.zeros(len(logg_table)) for i in range(len(logg_table)): logg_table["row_id"][i] = _spectype_to_id( logg_table["spec_type"][i], logg_table["spec_subtype"][i] ) # interpolate the columns with missing data nans = np.isnan(logg_table["I"]) logg_table["I"][nans] = np.interp( logg_table["row_id"][nans], logg_table["row_id"][~nans], logg_table["I"][~nans] ) nans = np.isnan(logg_table["III"]) logg_table["III"][nans] = np.interp( logg_table["row_id"][nans], logg_table["row_id"][~nans], logg_table["III"][~nans], ) # convert from g/g_sun units to cm/s^2 logg_sun = 27444 # cm/s^2 for col in ["I", "III", "V"]: logg_table[col] = np.log10(10 ** logg_table[col] * logg_sun) # now make a general interpolator logg_interp_function = RegularGridInterpolator( (np.array([1, 3, 5]), logg_table["row_id"]), np.array([logg_table["I"], logg_table["III"], logg_table["V"]]), ) return logg_interp_function
def _spectype_to_id(spec_type, spec_subtype): """ Convert a spectral type to a numerical ID, so they can be interpolated """ if (spec_subtype < 0) or (spec_subtype >= 10): raise ValueError("spec_subtype must be >=0 and <10") return ["O", "B", "A", "F", "G", "K", "M"].index(spec_type) * 10 + spec_subtype def _lumin_class_to_number(lumin_class): """ Convert the luminosity class from roman numeral to integer """ return ["I", "II", "III", "IV", "V"].index(lumin_class) + 1