Merge branch 'master' into fit_calib

species/__init__.py

@@ -65,6 +65,7 @@
 
 from species.util.read_util import add_luminosity, \
                                    get_mass, \
+                                   get_radius, \
                                    powerlaw_spectrum, \
                                    gaussian_spectrum, \
                                    update_spectra

species/analysis/compare_spectra.py

@@ -310,7 +310,7 @@ def compare_model(self,
         """
         Method for finding the best fitting spectrum from a grid of atmospheric model spectra by
         evaluating the goodness-of-fit statistic from Cushing et al. (2008). Currently, this method
-        only supports model grids with only :math:`T_\mathrm{eff}` and :math:`\log(g)` as free
+        only supports model grids with only :math:`T_\\mathrm{eff}` and :math:`\\log(g)` as free
         parameters (e.g. BT-Settl). Please create an issue on Github if support for models with
         more than two parameters is required.
 

species/analysis/emission_line.py

@@ -4,13 +4,19 @@
 
 import configparser
 import os
+import warnings
 
 from typing import Dict, Optional, Tuple, Union
 
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
-import ultranest
+
+try:
+    import ultranest
+except:
+    warnings.warn('UltraNest could not be imported. '
+                  'Perhaps because cython was not correctly compiled?')
 
 from astropy import units as u
 from astropy.modeling.fitting import LinearLSQFitter

species/analysis/fit_model.py

@@ -12,7 +12,12 @@
 import emcee
 import numpy as np
 import spectres
-import ultranest
+
+try:
+    import ultranest
+except:
+    warnings.warn('UltraNest could not be imported. '
+                  'Perhaps because cython was not correctly compiled?')
 
 # Installation of MultiNest is not possible on readthedocs
 try:

species/core/box.py

@@ -27,8 +27,13 @@ def create_box(boxtype,
         box.filter_mag = kwargs['filter_mag']
         box.color = kwargs['color']
         box.magnitude = kwargs['magnitude']
-        box.sptype = kwargs['sptype']
         box.names = kwargs['names']
+        if 'sptype' in kwargs:
+            box.sptype = kwargs['sptype']
+        if 'mass' in kwargs:
+            box.mass = kwargs['mass']
+        if 'radius' in kwargs:
+            box.radius = kwargs['radius']
 
     if boxtype == 'colorcolor':
         box = ColorColorBox()
@@ -37,8 +42,13 @@ def create_box(boxtype,
         box.filters = kwargs['filters']
         box.color1 = kwargs['color1']
         box.color2 = kwargs['color2']
-        box.sptype = kwargs['sptype']
         box.names = kwargs['names']
+        if 'sptype' in kwargs:
+            box.sptype = kwargs['sptype']
+        if 'mass' in kwargs:
+            box.mass = kwargs['mass']
+        if 'radius' in kwargs:
+            box.radius = kwargs['radius']
 
     elif boxtype == 'isochrone':
         box = IsochroneBox()
@@ -63,6 +73,7 @@ def create_box(boxtype,
         box = ObjectBox()
         box.name = kwargs['name']
         box.filters = kwargs['filters']
+        box.mean_wavel = kwargs['mean_wavel']
         box.magnitude = kwargs['magnitude']
         box.flux = kwargs['flux']
         box.distance = kwargs['distance']
@@ -171,8 +182,10 @@ def __init__(self):
         self.filter_mag = None
         self.color = None
         self.magnitude = None
-        self.sptype = None
         self.names = None
+        self.sptype = None
+        self.mass = None
+        self.radius = None
 
 
 class ColorColorBox(Box):
@@ -193,8 +206,10 @@ def __init__(self):
         self.filters = None
         self.color1 = None
         self.color2 = None
-        self.sptype = None
         self.names = None
+        self.sptype = None
+        self.mass = None
+        self.radius = None
 
 
 class IsochroneBox(Box):
@@ -303,6 +318,7 @@ def __init__(self):
 
         self.name = None
         self.filters = None
+        self.mean_wavel = None
         self.magnitude = None
         self.flux = None
         self.distance = None

species/core/constants.py

@@ -12,3 +12,4 @@
 M_EARTH = 5.9722e24  # (kg)
 R_EARTH = 6.3781e6  # (m)
 SIGMA_SB = 5.670374419e-8  # (W m−2 K−4)
+ATOMIC_MASS = 1.66053906660e-27  # (kg)

species/data/database.py

@@ -4,6 +4,7 @@
 
 import configparser
 import os
+import json
 import warnings
 
 from typing import Dict, List, Optional, Tuple, Union
@@ -103,6 +104,8 @@ def list_companions() -> None:
             None
         """
 
+        spec_data = companions.get_spec_data()
+
         for planet_name, planet_dict in companions.get_data().items():
             distance = planet_dict['distance']
             app_mag = planet_dict['app_mag']
@@ -113,6 +116,10 @@ def list_companions() -> None:
             for mag_name, mag_dict in app_mag.items():
                 print(f'{mag_name} (mag) = {mag_dict[0]} +/- {mag_dict[1]}')
 
+            if planet_name in spec_data:
+                for key, value in spec_data[planet_name].items():
+                    print(f'{key} spectrum from {value[3]}')
+
             print()
 
     @typechecked
@@ -145,8 +152,9 @@ def delete_data(self,
     def add_companion(self,
                       name: Union[Optional[str], Optional[List[str]]] = None) -> None:
         """
-        Function for adding the magnitudes of directly imaged planets and brown dwarfs from
-        :class:`~species.data.companions.get_data` to the database.
+        Function for adding the magnitudes and spectra of directly imaged planets and brown dwarfs
+        from :class:`~species.data.companions.get_data` and
+        :class:`~species.data.companions.get_comp_spec`to the database.
 
         Parameters
         ----------
@@ -170,9 +178,12 @@ def add_companion(self,
             name = data.keys()
 
         for item in name:
+            spec_dict = companions.companion_spectra(self.input_path, item)
+
             self.add_object(object_name=item,
                             distance=data[item]['distance'],
-                            app_mag=data[item]['app_mag'])
+                            app_mag=data[item]['app_mag'],
+                            spectrum=spec_dict)
 
     @typechecked
     def add_dust(self) -> None:
@@ -886,8 +897,8 @@ def add_object(self,
                 print(f'      - Filename: {value[0]}')
                 print(f'      - Data shape: {read_spec[key].shape}')
                 print(f'      - Wavelength range (um): {wavelength[0]:.2f} - {wavelength[-1]:.2f}')
-                print(f'      - Mean flux (W m-2 um-1): {np.mean(flux):.2e}')
-                print(f'      - Mean error (W m-2 um-1): {np.mean(error):.2e}')
+                print(f'      - Mean flux (W m-2 um-1): {np.nanmean(flux):.2e}')
+                print(f'      - Mean error (W m-2 um-1): {np.nanmean(error):.2e}')
 
                 if isinstance(deredden, float):
                     print(f'      - Dereddening A_V: {deredden}')
@@ -920,7 +931,7 @@ def add_object(self,
                             for i, hdu_item in enumerate(hdulist):
                                 data = np.asarray(hdu_item.data)
 
-                                corr_warn = f'The covariance matrix from {value[1]} contains ' \
+                                corr_warn = f'The matrix from {value[1]} contains ' \
                                             f'ones along the diagonal. Converting this ' \
                                             f'correlation matrix into a covariance matrix.'
 
@@ -957,7 +968,7 @@ def add_object(self,
                     print('   - Covariance matrix:')
 
                     if np.all(np.diag(data) == 1.):
-                        warnings.warn(f'The covariance matrix from {value[1]} contains ones on '
+                        warnings.warn(f'The matrix from {value[1]} contains ones on '
                                       f'the diagonal. Converting this correlation matrix into a '
                                       f'covariance matrix.')
 
@@ -1448,7 +1459,8 @@ def get_compare_sample(self,
             model_param['distance'] = dset.attrs['distance']
 
             if n_spec_name == 1:
-                model_param['radius'] = dset.attrs[f'radius_{item}']
+                spec_name = dset.attrs['spec_name0']
+                model_param['radius'] = dset.attrs[f'radius_{spec_name}']
 
             else:
                 if spec_fix is None:
@@ -1787,6 +1799,7 @@ def get_object(self,
 
             magnitude = {}
             flux = {}
+            mean_wavel = {}
 
             for observatory in dset.keys():
                 if observatory not in ['distance', 'spectrum']:
@@ -1797,13 +1810,17 @@ def get_object(self,
                             magnitude[name] = dset[name][0:2]
                             flux[name] = dset[name][2:4]
 
+                            filter_trans = read_filter.ReadFilter(name)
+                            mean_wavel[name] = filter_trans.mean_wavelength()
+
             phot_filters = list(magnitude.keys())
 
         else:
 
             magnitude = None
             flux = None
             phot_filters = None
+            mean_wavel = None
 
         if inc_spec and f'objects/{object_name}/spectrum' in h5_file:
             spectrum = {}
@@ -1835,6 +1852,7 @@ def get_object(self,
         return box.create_box('object',
                               name=object_name,
                               filters=phot_filters,
+                              mean_wavel=mean_wavel,
                               magnitude=magnitude,
                               flux=flux,
                               distance=distance,
@@ -1845,7 +1863,7 @@ def get_samples(self,
                     tag: str,
                     burnin: Optional[int] = None,
                     random: Optional[int] = None,
-                    out_file: Optional[str] = None) -> box.SamplesBox:
+                    json_file: Optional[str] = None) -> box.SamplesBox:
         """
         Parameters
         ----------
@@ -1858,10 +1876,9 @@ def get_samples(self,
         random : int, None
             Number of random samples to select. All samples (with the burnin excluded) are
             selected if set to ``None``.
-        out_file : str, None
-            Output file to store the posterior samples. The data will be stored in a FITS file if
-            the argument of ``out_file`` ends with `.fits`. Otherwise, the data will be written to
-            a text file. The data will not be written to a file if the argument is set to ``None``.
+        json_file : str, None
+            JSON file to store the posterior samples. The data will not be written if the argument
+            is set to ``None``.
 
         Returns
         -------
@@ -1916,18 +1933,14 @@ def get_samples(self,
 
         median_sample = self.get_median_sample(tag, burnin)
 
-        if out_file is not None:
-            header = ''
-            for i, item in enumerate(param):
-                header += f'{item}'
-                if i != len(param) - 1:
-                    header += ' - '
+        if json_file is not None:
+            samples_dict = {}
 
-            if out_file.endswith('.fits'):
-                fits.writeto(out_file, samples, overwrite=True)
+            for i, item in enumerate(param):
+                samples_dict[item] = list(samples[:, i])
 
-            else:
-                np.savetxt(out_file, samples, header=header)
+            with open(json_file, 'w') as out_file:
+                json.dump(samples_dict, out_file, indent=4)
 
         return box.create_box('samples',
                               spectrum=spectrum,
@@ -2091,12 +2104,15 @@ def add_comparison(self,
 
             # Indices of the best-fit model
             best_index = np.unravel_index(goodness_sum.argmin(), goodness_sum.shape)
+            dset.attrs['best_fit'] = goodness_sum[best_index]
 
             print('Best-fit parameters:')
+            print(f'   - Goodness-of-fit = {goodness_sum[best_index]:.2e}')
 
             for i, item in enumerate(model_param):
                 best_param = coord_points[i][best_index[i]]
                 dset.attrs[f'best_param{i}'] = best_param
+                print(f'   - {item} = {best_param}')
 
             for i, item in enumerate(spec_name):
                 scaling = flux_scaling[best_index[0], best_index[1], best_index[2], i]
@@ -2105,7 +2121,7 @@ def add_comparison(self,
                 radius /= constants.R_JUP  # (Rjup)
 
                 dset.attrs[f'radius_{item}'] = radius
-                print(f' -   {item} radius (Rjup) = {radius:.2f}')
+                print(f'   - {item} radius (Rjup) = {radius:.2f}')
 
                 dset.attrs[f'scaling_{item}'] = scaling
-                print(f' -   {item} scaling = {scaling:.2e}')
+                print(f'   - {item} scaling = {scaling:.2e}')

species/data/filters.py

@@ -70,6 +70,20 @@ def download_filter(filter_id: str) -> Tuple[Optional[np.ndarray],
 
         os.remove('VisAO_zp_filter_curve.dat')
 
+    elif filter_id == 'Keck/NIRC2.NB_4.05':
+        # The filter profile of Br_alpha has been digitized from
+        # https://www2.keck.hawaii.edu/inst/nirc2/filters.html
+
+        url = 'https://home.strw.leidenuniv.nl/~stolker/species/Keck_NIRC2.NB_4.05.dat'
+        urllib.request.urlretrieve(url, 'Keck_NIRC2.NB_4.05.dat')
+
+        wavelength, transmission = np.loadtxt('Keck_NIRC2.NB_4.05.dat', unpack=True)
+
+        # Not sure if energy- or photon-counting detector
+        det_type = 'photon'
+
+        os.remove('Keck_NIRC2.NB_4.05.dat')
+
     elif filter_id in ['LCO/VisAO.Ys', 'Magellan/VisAO.Ys']:
         url = 'https://xwcl.science/magao/visao/VisAO_Ys_filter_curve.dat'
         urllib.request.urlretrieve(url, 'VisAO_Ys_filter_curve.dat')