fit tools stellar spectrum

species/__init__.py

@@ -3,14 +3,15 @@
 from species.analysis.fit_spectrum import FitSpectrum
 
 from species.analysis.photometry import SyntheticPhotometry, \
-                                        apparent_to_absolute
+                                        apparent_to_absolute, \
+                                        multi_photometry
+
 
 from species.read.read_calibration import ReadCalibration
 
 from species.read.read_filter import ReadFilter
 
 from species.read.read_model import ReadModel, \
-                                    multi_photometry, \
                                     get_mass, \
                                     add_luminosity
 

species/analysis/fit_spectrum.py

@@ -2,19 +2,16 @@
 Text
 '''
 
-import os
 import sys
 import math
-import configparser
 
-import h5py
 import emcee
 import progress.bar
 import numpy as np
 
 from species.analysis import photometry
 from species.data import database
-from species.read import read_model, read_object, read_calibration
+from species.read import read_object, read_calibration
 
 
 MIN_CHISQ = np.inf
@@ -26,10 +23,10 @@ def lnprob(param,
            objphot,
            specphot):
     '''
-    :param param: Parameter values.
+    :param param: Values of the offset and scaling parameter.
     :type param: numpy.ndarray
-    :param bounds: Parameter boundaries.
-    :type bounds: tuple(float, float)
+    :param bounds: Boundaries of the offset and scaling parameter
+    :type bounds: dict
     :param objphot:
     :type objphot:
     :param specphot:
@@ -42,7 +39,8 @@ def lnprob(param,
     global MIN_CHISQ
     global MIN_PARAM
 
-    if bounds[0] <= param <= bounds[1]:
+    if bounds['offset'][0] <= param[0] <= bounds['offset'][1] and \
+       bounds['scaling'][0] <= param[1] <= bounds['scaling'][1]:
         ln_prior = 0.
 
     else:
@@ -53,12 +51,12 @@ def lnprob(param,
 
     else:
         chisq = 0.
-        for i, item in enumerate(objphot):
-            chisq += (objphot[i][0]-param*specphot[0])**2 / objphot[i][1]**2
+        for i, _ in enumerate(objphot):
+            chisq += (objphot[i][0]-(param[0]+param[1]*specphot[i]))**2 / objphot[i][1]**2
 
         if chisq < MIN_CHISQ:
             MIN_CHISQ = chisq
-            MIN_PARAM = {'scaling':param}
+            MIN_PARAM = {'offset':param[0], 'scaling':param[1]}
 
         ln_prob = ln_prior - 0.5*chisq
 
@@ -83,8 +81,9 @@ def __init__(self,
         :type filters: tuple(str, )
         :param spectrum: Calibration spectrum.
         :type spectrum: str
-        :param bounds: Range of the scaling parameter (min, max).
-        :type bounds: tuple(float, float)
+        :param bounds: Boundaries of the offset and scaling parameter, as
+                       {'offset':(min, max), 'scaling':(min, max)}.
+        :type bounds: dict
 
         :return: None
         '''
@@ -113,7 +112,7 @@ def __init__(self,
             obj_phot = self.object.get_photometry(item)
             self.objphot.append((obj_phot[2], obj_phot[3]))
 
-        self.modelpar = ['scaling']
+        self.modelpar = ['offset', 'scaling']
 
     def run_mcmc(self,
                  nwalkers,
@@ -125,8 +124,9 @@ def run_mcmc(self,
         :type nwalkers: int
         :param nsteps: Number of steps for each walker.
         :type nsteps: int
-        :param guess: Guess of the scaling factor.
-        :type guess: float
+        :param guess: Guess of the offset and scaling parameter, as
+                      {'offset':(guess), 'scaling':(guess)}. Non-zero values work best.
+        :type guess: dict
         :param tag: Database tag for the results.
         :type tag: int
 
@@ -139,10 +139,12 @@ def run_mcmc(self,
         sys.stdout.write('Running MCMC...')
         sys.stdout.flush()
 
-        ndim = 1
+        ndim = 2
 
         initial = np.zeros((nwalkers, ndim))
-        initial[:, 0] = guess + np.random.normal(0, 1e-3*guess, nwalkers)
+
+        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)
 
         sampler = emcee.EnsembleSampler(nwalkers=nwalkers,
                                         dim=ndim,
@@ -156,7 +158,7 @@ def run_mcmc(self,
                                    max=nsteps,
                                    suffix='%(percent)d%%')
 
-        for i, _ in enumerate(sampler.sample(initial, iterations=nsteps)):
+        for _ in sampler.sample(initial, iterations=nsteps):
             progbar.next()
 
         progbar.finish()

species/analysis/photometry.py

@@ -11,7 +11,41 @@
 from scipy.integrate import simps
 
 from species.data import database
-from species.read import read_filter, read_calibration
+from species.core import box
+from species.read import read_filter, read_calibration, read_model
+
+
+def multi_photometry(datatype,
+                     spectrum,
+                     filters,
+                     parameters):
+    '''
+    :param datatype: Data type ('model' or 'calibration').
+    :type datatype: str
+    :param spectrum: Spectrum name (e.g., 'drift-phoenix').
+    :type spectrum: str
+    :param filters: Filter IDs.
+    :type filters: tuple(str, )
+    :param parameters: Model parameter values.
+    :type parameters: dict
+
+    :return: Box with synthetic photometry.
+    :rtype: species.core.box.SynphotBox
+    '''
+
+    flux = {}
+
+    if datatype == 'model':
+        for item in filters:
+            readmodel = read_model.ReadModel(spectrum, item)
+            flux[item] = readmodel.get_photometry(parameters, ('specres', 100.))
+
+    elif datatype == 'calibration':
+        for item in filters:
+            readcalib = read_calibration.ReadCalibration(spectrum, item)
+            flux[item] = readcalib.get_photometry(parameters)
+
+    return box.create_box('synphot', name='synphot', flux=flux)
 
 
 def apparent_to_absolute(app_mag,

species/core/box.py

@@ -9,8 +9,8 @@ def open_box(box):
     '''
 
     for item in box.__dict__.keys():
-        print(item, '=', box.__dict__[item], '\n')
-
+        sys.stdout.write(str(item)+' = '+str(box.__dict__[item])+'\n')
+        sys.stdout.flush()
 
 def create_box(boxtype, **kwargs):
     '''

species/data/database.py

@@ -485,7 +485,7 @@ def add_samples(self,
         try:
             int_auto = emcee.autocorr.integrated_time(sampler.flatchain)
 
-            sys.stdout.write('Integrated autocorrelation time ='+str(int_auto)+'\n')
+            sys.stdout.write('Integrated autocorrelation time = '+str(int_auto)+'\n')
             sys.stdout.flush()
 
         except emcee.autocorr.AutocorrError:
@@ -585,14 +585,13 @@ def get_mcmc_spectra(self,
                 model_par['distance'] = distance
 
             if spectrum_type == 'model':
-                box = readmodel.get_model(model_par, sampling)
+                specbox = readmodel.get_model(model_par, sampling)
             elif spectrum_type == 'calibration':
-                box = readcalib.get_spectrum()
-                box.flux *= model_par['scaling']
+                specbox = readcalib.get_spectrum(model_par)
 
             box.type = 'mcmc'
 
-            boxes.append(box)
+            boxes.append(specbox)
 
         h5_file.close()
 

species/plot/plot_spectrum.py

@@ -234,7 +234,6 @@ def plot_spectrum(boxes,
                         elif item in ('$\log\,g$', '$R$', '$M$', '[Fe/H]'):
                             value = '{:.2f}'.format(par_val[i])
                         elif item == '$L$':
-                            # print(item, par_val[i], par_key[i])
                             value = '{0:.1e}'.format(par_val[i])
                         else:
                             continue

species/read/read_calibration.py

@@ -10,6 +10,7 @@
 
 from scipy.interpolate import interp1d
 
+from species.analysis import photometry
 from species.core import box
 from species.read import read_filter
 
@@ -48,9 +49,26 @@ def __init__(self,
 
         self.database = config['species']['database']
 
+    def interpolate(self):
+        '''
+        :return: Interpolated spectrum.
+        :rtype: scipy.interpolate.interpolate.interp1d
+        '''
+
+        calibbox = self.get_spectrum()
+
+        return interp1d(calibbox.wavelength,
+                        calibbox.flux,
+                        kind='cubic',
+                        bounds_error=False,
+                        fill_value=float('nan'))
+
     def get_spectrum(self,
+                     model_par=None,
                      negative=False):
         '''
+        :param model_par: Model parameter values. Not used if set to None.
+        :type model_par: dict
         :param negative: Include negative values.
         :type negative: bool
 
@@ -69,6 +87,9 @@ def get_spectrum(self,
 
         h5_file.close()
 
+        if model_par:
+            data[1, ] = model_par['offset'] + model_par['scaling']*data[1, ]
+
         return box.create_box(boxtype='spectrum',
                               spectrum='calibration',
                               wavelength=data[0, ],
@@ -78,16 +99,18 @@ def get_spectrum(self,
                               sptype=None,
                               distance=None)
 
-    def interpolate(self):
+    def get_photometry(self,
+                       model_par,
+                       synphot=None):
         '''
-        :return: Interpolated spectrum.
-        :rtype: scipy.interpolate.interpolate.interp1d
+        :param model_par: Model parameter values.
+        :type model_par: dict
+
+        :return: Average flux density (W m-2 micron-1).
+        :rtype: float
         '''
 
-        calibbox = self.get_spectrum()
+        specbox = self.get_spectrum(model_par)
+        synphot = photometry.SyntheticPhotometry(self.filter_name)
 
-        return interp1d(calibbox.wavelength,
-                        calibbox.flux,
-                        kind='cubic',
-                        bounds_error=False,
-                        fill_value=float('nan'))
+        return synphot.spectrum_to_photometry(specbox.wavelength, specbox.flux)

species/read/read_model.py

@@ -19,29 +19,6 @@
 from species.read import read_filter
 
 
-def multi_photometry(model,
-                     filters,
-                     parameters):
-    '''
-    :param model: Model name.
-    :type model: str
-    :param filters: Filter IDs.
-    :type filters: tuple(str, )
-    :param parameters: Model parameter values.
-    :type parameters: dict
-
-    :return: Box with synthetic photometry.
-    :rtype: species.core.box.SynphotBox
-    '''
-
-    flux = {}
-    for item in filters:
-        readmodel = ReadModel(model, item)
-        flux[item] = readmodel.get_photometry(parameters, ('specres', 100.))
-
-    return box.create_box('synphot', name='synphot', flux=flux)
-
-
 def get_mass(model_par):
     '''
     :param model_par: Model parameter values. Should contain the surface gravity and radius.