mamajek table

tomasstolker · Mar 13, 2019 · 79aa5bc · 79aa5bc
1 parent e3bedb0
commit 79aa5bc
Show file tree

Hide file tree

Showing 8 changed files with 205 additions and 34 deletions.
diff --git a/species/analysis/fit_spectrum.py b/species/analysis/fit_spectrum.py
@@ -23,9 +23,9 @@ def lnprob(param,
            objphot,
            specphot):
     '''
-    :param param: Values of the offset and scaling parameter.
+    :param param: Value of the scaling parameter.
     :type param: numpy.ndarray
-    :param bounds: Boundaries of the offset and scaling parameter
+    :param bounds: Boundary of the scaling parameter
     :type bounds: dict
     :param objphot:
     :type objphot:
@@ -39,8 +39,7 @@ def lnprob(param,
     global MIN_CHISQ
     global MIN_PARAM
 
-    if bounds['offset'][0] <= param[0] <= bounds['offset'][1] and \
-       bounds['scaling'][0] <= param[1] <= bounds['scaling'][1]:
+    if bounds['scaling'][0] <= param <= bounds['scaling'][1]:
         ln_prior = 0.
 
     else:
@@ -52,11 +51,11 @@ def lnprob(param,
     else:
         chisq = 0.
         for i, _ in enumerate(objphot):
-            chisq += (objphot[i][0]-(param[0]+param[1]*specphot[i]))**2 / objphot[i][1]**2
+            chisq += (objphot[i][0]-(param*specphot[i]))**2 / objphot[i][1]**2
 
         if chisq < MIN_CHISQ:
             MIN_CHISQ = chisq
-            MIN_PARAM = {'offset':param[0], 'scaling':param[1]}
+            MIN_PARAM = {'scaling':param}
 
         ln_prob = ln_prior - 0.5*chisq
 
@@ -112,7 +111,7 @@ def __init__(self,
             obj_phot = self.object.get_photometry(item)
             self.objphot.append((obj_phot[2], obj_phot[3]))
 
-        self.modelpar = ['offset', 'scaling']
+        self.modelpar = ['scaling']
 
     def run_mcmc(self,
                  nwalkers,
@@ -124,9 +123,8 @@ def run_mcmc(self,
         :type nwalkers: int
         :param nsteps: Number of steps for each walker.
         :type nsteps: int
-        :param guess: Guess of the offset and scaling parameter, as
-                      {'offset':(guess), 'scaling':(guess)}. Non-zero values work best.
-        :type guess: dict
+        :param guess: Guess of the scaling parameter.
+        :type guess: float
         :param tag: Database tag for the results.
         :type tag: int
 
@@ -139,12 +137,10 @@ def run_mcmc(self,
         sys.stdout.write('Running MCMC...')
         sys.stdout.flush()
 
-        ndim = 2
+        ndim = 1
 
         initial = np.zeros((nwalkers, ndim))
-
-        initial[:, 0] = guess['offset'] + np.random.normal(0, 1e-1*guess['offset'], nwalkers)
-        initial[:, 1] = guess['scaling'] + np.random.normal(0, 1e-1*guess['scaling'], nwalkers)
+        initial[:, 0] = guess['scaling'] + np.random.normal(0, 1e-1*guess['scaling'], nwalkers)
 
         sampler = emcee.EnsembleSampler(nwalkers=nwalkers,
                                         dim=ndim,

diff --git a/species/data/companions.py b/species/data/companions.py
@@ -134,9 +134,9 @@ def get_data():
                                    'Paranal/NACO.J':(12.47, 0.20), # Biller et al. 2010
                                    'Paranal/NACO.H':(11.93, 0.14), # Biller et al. 2010
                                    'Paranal/NACO.Ks':(11.53, 0.07), # Biller et al. 2010
-                                   # 'Paranal/NACO.Lp':(11.05, 0.18), # Beust et al. 2015
-                                   'Paranal/NACO.NB405':(9.81, 0.01), # Stolker et al. in prep.
-                                   'Paranal/NACO.Mp':(9.42, 0.02), # Stolker et al. in prep.
+                                   'Paranal/NACO.Lp':(11.05, 0.18), # Beust et al. 2015
+                                   'Paranal/NACO.NB405':(6.48+4.56, 0.01), # Stolker et al. in prep.
+                                   'Paranal/NACO.Mp':(6.55+4.66, 0.02), # Stolker et al. in prep.
                                    'Gemini/NICI.ED286':(11.68, 0.14), # Biller et al. 2010
                                    'Gemini/NIRI.H2S1v2-1-G0220':(11.39, 0.14)}}, # Biller et al. 2010
 

diff --git a/species/data/database.py b/species/data/database.py
@@ -17,7 +17,7 @@
 from species.core import box, constants
 from species.read import read_model, read_calibration
 from species.data import drift_phoenix, btnextgen, vega, irtf, spex, vlm_plx, leggett, \
-                         companions, filters, util
+                         companions, filters, util, mamajek
 
 
 warnings.simplefilter('ignore', UserWarning)
@@ -268,6 +268,9 @@ def add_photometry(self,
         elif library[0:7] == 'leggett':
             leggett.add_leggett(self.input_path, h5_file)
 
+        elif library[0:7] == 'mamajek':
+            mamajek.add_mamajek(self.input_path, h5_file)
+
         h5_file.close()
 
     def add_calibration(self,

diff --git a/species/data/mamajek.py b/species/data/mamajek.py
@@ -0,0 +1,106 @@
+'''
+Text
+'''
+
+import os
+import sys
+
+from urllib.request import urlretrieve
+
+import h5py
+import numpy as np
+import pandas as pd
+
+from species.data import util
+
+
+def add_mamajek(input_path,
+                database):
+    '''
+    Function for adding "A Modern Mean Dwarf Stellar Color and Effective Temperature Sequence".
+    http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt
+
+    :param input_path:
+    :type input_path: str
+    :param database:
+    :type database: h5py._hl.files.File
+
+    :return: None
+    '''
+
+    data_file = os.path.join(input_path, 'EEM_dwarf_UBVIJHK_colors_Teff.txt')
+
+    url = 'http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt'
+
+    if not os.path.isfile(data_file):
+        sys.stdout.write('Downloading Stellar Colors and Effective Temperatures (53 kB)...')
+        sys.stdout.flush()
+
+        urlretrieve(url, data_file)
+
+        sys.stdout.write(' [DONE]\n')
+        sys.stdout.flush()
+
+    sys.stdout.write('Adding Stellar Colors and Effective Temperatures...')
+    sys.stdout.flush()
+
+    group = 'photometry/mamajek'
+
+    database.create_group(group)
+
+    dataframe = pd.read_csv(data_file, delimiter=r"\s+", nrows=124, header=22, \
+                            na_values=['...', '....', '.....', ',..', '19.52:'],
+                            dtype={'SpT':str, 'Teff':float, 'logT':float, 'BCv':float, \
+                            'Mv':float, 'logL':float, 'B-V':float, 'Bt-Vt':float, 'U-B':float, \
+                            'V-Rc':float, 'V-Ic':float, 'V-Ks':float, 'J-H':float, 'H-Ks':float, \
+                            'Ks-W1':float, 'W1-W2':float, 'W1-W3':float, 'W1-W4':float, \
+                            'Msun':float, 'logAge':str, 'b-y':float, 'M_J':float, 'M_Ks':float, \
+                            'Mbol':float, 'i-z':float, 'z-Y':float, 'R_Rsun':float, 'SpT':str, \
+                            'G-V':float, 'Bp-Rp':float, 'M_G':float, 'G-Rp':float})
+
+    dataframe.columns = dataframe.columns.str.replace('#', '')
+
+    sptype = np.asarray(dataframe['SpT'])
+    sptype = util.update_sptype(sptype)
+
+    flag = np.repeat('star', sptype.size)
+    distance = np.repeat(10., sptype.size) # [pc]
+
+    v_mag = dataframe['Mv'] # [mag]
+    b_mag = dataframe['B-V']+dataframe['Mv'] # [mag]
+    u_mag = dataframe['U-B']+b_mag # [mag]
+
+    j_mag = dataframe['M_J'] # [mag]
+    ks_mag = dataframe['M_Ks'] # [mag]
+    h_mag = dataframe['H-Ks']+ks_mag # [mag]
+
+    w1_mag = -dataframe['Ks-W1']+ks_mag # [mag]
+    w2_mag = -dataframe['W1-W2']+w1_mag # [mag]
+    w3_mag = -dataframe['W1-W3']+w1_mag # [mag]
+    w4_mag = -dataframe['W1-W4']+w1_mag # [mag]
+
+    dtype = h5py.special_dtype(vlen=bytes)
+
+    dset = database.create_dataset(group+'/sptype', (np.size(sptype), ), dtype=dtype)
+    dset[...] = sptype
+
+    dset = database.create_dataset(group+'/flag', (np.size(flag), ), dtype=dtype)
+    dset[...] = flag
+
+    database.create_dataset(group+'/distance', data=distance, dtype='f')
+    database.create_dataset(group+'/Teff', data=dataframe['Teff'], dtype='f')
+    database.create_dataset(group+'/U', data=u_mag, dtype='f')
+    database.create_dataset(group+'/B', data=b_mag, dtype='f')
+    database.create_dataset(group+'/V', data=v_mag, dtype='f')
+    database.create_dataset(group+'/J', data=j_mag, dtype='f')
+    database.create_dataset(group+'/H', data=h_mag, dtype='f')
+    database.create_dataset(group+'/Ks', data=ks_mag, dtype='f')
+    database.create_dataset(group+'/W1', data=w1_mag, dtype='f')
+    database.create_dataset(group+'/W2', data=w2_mag, dtype='f')
+    database.create_dataset(group+'/W3', data=w3_mag, dtype='f')
+    database.create_dataset(group+'/W4', data=w4_mag, dtype='f')
+
+    sys.stdout.write(' [DONE]\n')
+    sys.stdout.flush()
+
+    database.close()
diff --git a/species/data/util.py b/species/data/util.py
@@ -16,7 +16,7 @@ def update_sptype(sptypes):
     :rtype: numpy.ndarray
     '''
 
-    mlty = ('M', 'L', 'T', 'Y')
+    mlty = ('O', 'B', 'A', 'F', 'G', 'K', 'M', 'L', 'T', 'Y')
 
     for i, spt in enumerate(sptypes):
         if spt == 'None':

diff --git a/species/plot/plot_color.py b/species/plot/plot_color.py
@@ -94,8 +94,12 @@ def plot_color_magnitude(colorbox,
     if ylim:
         ax1.set_ylim(ylim[0], ylim[1])
 
+    if colorbox.object_type == 'star':
+        bounds = np.arange(0, 11, 1)
+    else:
+        bounds = np.arange(0, 8, 1)
+
     cmap = plt.cm.viridis
-    bounds = np.arange(0, 8, 1)
     norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
 
     sptype = colorbox.sptype
@@ -108,9 +112,13 @@ def plot_color_magnitude(colorbox,
     color = color[indices]
     magnitude = magnitude[indices]
 
-    spt_disc = util.sptype_discrete(sptype, color.shape)
+    if colorbox.object_type == 'star':
+        spt_disc = util.sptype_stellar(sptype, color.shape)
+        unique = np.arange(0, color.size, 1)
 
-    _, unique = np.unique(color, return_index=True)
+    else:
+        spt_disc = util.sptype_substellar(sptype, color.shape)
+        _, unique = np.unique(color, return_index=True)
 
     sptype = sptype[unique]
     color = color[unique]
@@ -124,8 +132,14 @@ def plot_color_magnitude(colorbox,
                   ticklocation='top', format='%.2f')
 
     cb.ax.tick_params(width=0.8, length=5, labelsize=10, direction='in', color='white')
-    cb.set_ticks(np.arange(0.5, 7., 1.))
-    cb.set_ticklabels(['M0-M4', 'M5-M9', 'L0-L4', 'L5-L9', 'T0-T4', 'T6-T8', 'Y1-Y2'])
+
+    if colorbox.object_type == 'star':
+        cb.set_ticks(np.arange(0.5, 10., 1.))
+        cb.set_ticklabels(['O', 'B', 'A', 'F', 'G', 'K', 'M', 'L', 'T', 'Y'])
+
+    else:
+        cb.set_ticks(np.arange(0.5, 7., 1.))
+        cb.set_ticklabels(['M0-M4', 'M5-M9', 'L0-L4', 'L5-L9', 'T0-T4', 'T6-T8', 'Y1-Y2'])
 
     if objects is not None:
         for item in objects:

diff --git a/species/plot/util.py b/species/plot/util.py
@@ -5,8 +5,8 @@
 import numpy as np
 
 
-def sptype_discrete(sptype,
-                    shape):
+def sptype_substellar(sptype,
+                      shape):
     '''
     :param sptype:
     :type sptype:
@@ -23,25 +23,77 @@ def sptype_discrete(sptype,
         sp = item[0:2]
 
         if sp in (np.string_('M0'), np.string_('M1'), np.string_('M2'), np.string_('M3'), np.string_('M4')):
-            spt_disc[i] = 0
+            spt_disc[i] = 0.5
 
         elif sp in (np.string_('M5'), np.string_('M6'), np.string_('M7'), np.string_('M8'), np.string_('M9')):
-            spt_disc[i] = 1
+            spt_disc[i] = 1.5
 
         elif sp in (np.string_('L0'), np.string_('L1'), np.string_('L2'), np.string_('L3'), np.string_('L4')):
-            spt_disc[i] = 2
+            spt_disc[i] = 2.5
 
         elif sp in (np.string_('L5'), np.string_('L6'), np.string_('L7'), np.string_('L8'), np.string_('L9')):
-            spt_disc[i] = 3
+            spt_disc[i] = 3.5
 
         elif sp in (np.string_('T0'), np.string_('T1'), np.string_('T2'), np.string_('T3'), np.string_('T4')):
-            spt_disc[i] = 4
+            spt_disc[i] = 4.5
 
         elif sp in (np.string_('T5'), np.string_('T6'), np.string_('T7'), np.string_('T8'), np.string_('T9')):
-            spt_disc[i] = 5
+            spt_disc[i] = 5.5
 
         elif np.string_('Y') in item:
-            spt_disc[i] = 6
+            spt_disc[i] = 6.5
+
+        else:
+            spt_disc[i] = np.nan
+            continue
+
+    return spt_disc
+
+
+def sptype_stellar(sptype,
+                   shape):
+    '''
+    :param sptype:
+    :type sptype:
+    :param shape:
+    :type shape:
+
+    :return::
+    :rtype: numpy.ndarray
+    '''
+
+    spt_disc = np.zeros(shape)
+
+    for i, item in enumerate(sptype):
+        if str(item)[2] == 'O':
+            spt_disc[i] = 0.5
+
+        elif str(item)[2] == 'B':
+            spt_disc[i] = 1.5
+
+        elif str(item)[2] == 'A':
+            spt_disc[i] = 2.5
+
+        elif str(item)[2] == 'F':
+            spt_disc[i] = 3.5
+
+        elif str(item)[2] == 'G':
+            spt_disc[i] = 4.5
+
+        elif str(item)[2] == 'K':
+            spt_disc[i] = 5.5
+
+        elif str(item)[2] == 'M':
+            spt_disc[i] = 6.5
+
+        elif str(item)[2] == 'L':
+            spt_disc[i] = 7.5
+
+        elif str(item)[2] == 'T':
+            spt_disc[i] = 8.5
+
+        elif str(item)[2] == 'Y':
+            spt_disc[i] = 9.5
 
         else:
             spt_disc[i] = np.nan

diff --git a/species/read/read_calibration.py b/species/read/read_calibration.py
@@ -104,7 +104,7 @@ def get_spectrum(self,
         h5_file.close()
 
         if parameters:
-            flux = parameters['offset'] + parameters['scaling']*flux
+            flux = parameters['scaling']*flux
 
         if extrapolate:
             def _power_law(wavelength, offset, scaling, power_index):