diff --git a/species/analysis/fit_spectrum.py b/species/analysis/fit_spectrum.py
index 92731b1d..986efde4 100644
--- a/species/analysis/fit_spectrum.py
+++ b/species/analysis/fit_spectrum.py
@@ -177,7 +177,7 @@ def run_mcmc(self,
                 self.modelpar.append('scaling'+str(i))
                 self.bounds['scaling'+str(i)] = (0., 1e2)
 
-        with Pool(processes=cpu_count()) as pool:
+        with Pool(processes=cpu_count()):
             sampler = emcee.EnsembleSampler(nwalkers,
                                             ndim,
                                             lnprob,
diff --git a/species/data/ames_cond.py b/species/data/ames_cond.py
index 44c3355c..b94b1f20 100644
--- a/species/data/ames_cond.py
+++ b/species/data/ames_cond.py
@@ -95,7 +95,7 @@ def add_ames_cond(input_path,
                         continue
 
                 print_message = f'Adding AMES-Cond model spectra... {filename}'
-                print(f'\r{print_message:<75}', end='')
+                print(f'\r{print_message:<71}', end='')
 
                 data_wavel = []
                 data_flux = []
@@ -182,4 +182,4 @@ def add_ames_cond(input_path,
     data_util.write_data('ames-cond', ['teff', 'logg'], database, data_sorted)
 
     print_message = 'Adding AMES-Cond model spectra... [DONE]'
-    print(f'\r{print_message:<75}')
+    print(f'\r{print_message:<71}')
diff --git a/species/data/database.py b/species/data/database.py
index 680a4124..bbefb887 100644
--- a/species/data/database.py
+++ b/species/data/database.py
@@ -368,15 +368,13 @@ def add_object(self,
         app_mag : dict
             Apparent magnitudes. Not written if set to None.
         spectrum : dict
-            Dictionary with spectra and correlation matrices. Multiple spectra can be included and
+            Dictionary with spectra and covariance matrices. Multiple spectra can be included and
             the files have to be in the FITS or ASCII format. The spectra should have 3 columns
             with wavelength (micron), flux density (W m-2 micron-1), and error (W m-2 micron-1).
-            The correlation matrix should be 2D with the same number of wavelength points as the
-            spectrum. For example, {'sphere_ifs': ('ifs_spectrum.dat', 'ifs_correlation.fits')}.
-            No data is stored for a correlation matrix if set to None, for example, {'sphere_ifs':
+            The covariance matrix should be 2D with the same number of wavelength points as the
+            spectrum. For example, {'sphere_ifs': ('ifs_spectrum.dat', 'ifs_covariance.fits')}.
+            No covariance data is stored if set to None, for example, {'sphere_ifs':
             ('ifs_spectrum.dat', None)}. The ``spectrum`` parameter is ignored if set to None.
-            The correlation matrix and uncertainties are used to compute the covariance matrix.
-            The covariance matrix is stored in the database.
 
         Returns
         -------
@@ -480,11 +478,20 @@ def add_object(self,
                             if data.ndim == 2 and data.shape[0] == data.shape[1]:
                                 if key not in read_cov:
                                     if data.shape[0] == read_spec[key].shape[0]:
-                                        read_cov[key] = data_util.correlation_to_covariance(
-                                            data, read_spec[key][:, 2])
+                                        if np.all(np.diag(data) == 1.):
+                                            warnings.warn(f'The covariance matrix from {value[1]} '
+                                                          f'contains ones on the diagonal. '
+                                                          f'Converting this correlation matrix '
+                                                          f'into a covariance matrix.')
+
+                                            read_cov[key] = data_util.correlation_to_covariance(
+                                                data, read_spec[key][:, 2])
+
+                                        else:
+                                            read_cov[key] = data
 
                         if key not in read_cov:
-                            raise ValueError(f'The correlation matrix from {value[1]} can not be '
+                            raise ValueError(f'The covariance matrix from {value[1]} can not be '
                                              f'read. The data format should be 2D with the same '
                                              f'number of wavelength points as the spectrum.')
 
@@ -492,23 +499,32 @@ def add_object(self,
                     try:
                         data = np.loadtxt(value[1])
                     except UnicodeDecodeError:
-                        raise ValueError(f'The correlation matrix from {value[1]} can not be read. '
+                        raise ValueError(f'The covariance matrix from {value[1]} can not be read. '
                                          f'Please provide a FITS or ASCII file.')
 
                     if data.ndim != 2 or 3 not in data.shape:
-                        raise ValueError(f'The correlation matrix from {value[1]} can not be read. '
+                        raise ValueError(f'The covariance matrix from {value[1]} can not be read. '
                                          f'The data format should be 2D with the same number of '
                                          f'wavelength points as the spectrum.')
 
-                    read_cov[key] = data_util.correlation_to_covariance(
-                        data, read_spec[key][:, 2])
+                    if np.all(np.diag(data) == 1.):
+                        warnings.warn(f'The covariance matrix from {value[1]} contains ones on '
+                                      f'the diagonal. Converting this correlation matrix into a '
+                                      f'covariance matrix.')
+
+                        read_cov[key] = data_util.correlation_to_covariance(
+                            data, read_spec[key][:, 2])
+
+                    else:
+                        read_cov[key] = data
 
             for key, value in read_spec.items():
-                dset = h5_file.create_dataset(f'objects/{object_name}/spectrum/{key}/spectrum',
-                                              data=read_spec[key])
+                h5_file.create_dataset(f'objects/{object_name}/spectrum/{key}/spectrum',
+                                       data=read_spec[key])
 
-                dset = h5_file.create_dataset(f'objects/{object_name}/spectrum/{key}/covariance',
-                                              data=read_cov[key])
+                if read_cov[key] is not None:
+                    h5_file.create_dataset(f'objects/{object_name}/spectrum/{key}/covariance',
+                                           data=read_cov[key])
 
         print(' [DONE]')
 
@@ -1067,7 +1083,7 @@ def get_object(self,
             for item in h5_file[f'objects/{object_name}/spectrum']:
                 data_group = f'objects/{object_name}/spectrum/{item}'
 
-                if h5_file[f'{data_group}/covariance'].shape is None:
+                if f'{data_group}/covariance' not in h5_file:
                     spectrum[item] = (np.asarray(h5_file[f'{data_group}/spectrum']), None)
                 else:
                     spectrum[item] = (np.asarray(h5_file[f'{data_group}/spectrum']),
diff --git a/species/read/read_calibration.py b/species/read/read_calibration.py
index aa777ed6..3f08a2ca 100644
--- a/species/read/read_calibration.py
+++ b/species/read/read_calibration.py
@@ -10,7 +10,6 @@
 import numpy as np
 
 from scipy.optimize import curve_fit
-from scipy.interpolate import interp1d
 
 from species.analysis import photometry
 from species.core import box
diff --git a/species/read/read_model.py b/species/read/read_model.py
index d12bb1c9..e6e5f0f1 100644
--- a/species/read/read_model.py
+++ b/species/read/read_model.py
@@ -11,7 +11,7 @@
 import spectres
 import numpy as np
 
-from scipy.interpolate import RegularGridInterpolator, interp1d
+from scipy.interpolate import RegularGridInterpolator
 
 from species.analysis import photometry
 from species.core import box, constants
@@ -342,7 +342,7 @@ def get_model(self,
                                                  old_spec_wavs=self.wl_points,
                                                  spec_fluxes=flux)
 
-                    except:
+                    except ValueError:
                         index_error = True
 
                     if not index_error:
@@ -408,7 +408,7 @@ def get_data(self,
             Box with the model spectrum.
         """
 
-        for key, value in model_param.items():
+        for key in model_param:
             if key not in ['radius', 'distance', 'mass', 'luminosity']:
                 if key not in self.get_parameters():
                     raise ValueError(f'The \'{key}\' parameter is not required by \'{self.model}\''
diff --git a/species/read/read_object.py b/species/read/read_object.py
index e2d3d91e..b4731023 100644
--- a/species/read/read_object.py
+++ b/species/read/read_object.py
@@ -84,7 +84,7 @@ def get_spectrum(self):
                 for item in h5_file[f'objects/{self.object_name}/spectrum']:
                     data_group = f'objects/{self.object_name}/spectrum/{item}'
 
-                    if h5_file[f'{data_group}/covariance'].shape is None:
+                    if f'{data_group}/covariance' not in h5_file:
                         spectrum[item] = (np.asarray(h5_file[f'{data_group}/spectrum']), None)
                     else:
                         spectrum[item] = (np.asarray(h5_file[f'{data_group}/spectrum']),
@@ -95,26 +95,6 @@ def get_spectrum(self):
 
         return spectrum
 
-    # def get_instrument(self):
-    #     """
-    #     Function for extracting the instrument name of the spectrum.
-    #
-    #     Returns
-    #     -------
-    #     str
-    #         Instrument that was used for the spectrum.
-    #     """
-    #
-    #     with h5py.File(self.database, 'r') as h5_file:
-    #         if 'objects/'+self.object_name+'/spectrum' in h5_file:
-    #             dset = h5_file['objects/'+self.object_name+'/spectrum']
-    #             instrument = dset.attrs['instrument']
-    #
-    #         else:
-    #             instrument = None
-    #
-    #     return instrument
-
     def get_distance(self):
         """
         Function for extracting the distance to the object.
diff --git a/species/read/read_planck.py b/species/read/read_planck.py
index a6713002..b649d79a 100644
--- a/species/read/read_planck.py
+++ b/species/read/read_planck.py
@@ -87,6 +87,35 @@ def planck(wavel_points,
 
         return 1e-6 * 4.*math.pi * scaling * planck_1/planck_2  # [W m-2 micron-1]
 
+    @staticmethod
+    def update_parameters(model_param):
+        """
+        Internal function for updating the dictionary with model parameters.
+
+        Parameters
+        ----------
+        model_param : dict
+            Dictionary with the 'teff' (K), 'radius' (Rjup), and 'distance' (pc). The values
+            of 'teff' and 'radius' can be a single float, or a list with floats for a combination
+            of multiple Planck functions, e.g.
+            {'teff': [1500., 1000.], 'radius': [1., 2.], 'distance': 10.}.
+
+        Returns
+        -------
+        dict
+            Updated dictionary with model parameters.
+        """
+
+        updated_param = {}
+
+        for i, _ in enumerate(model_param['teff']):
+            updated_param[f'teff_{i}'] = model_param['teff'][i]
+            updated_param[f'radius_{i}'] = model_param['radius'][i]
+
+        updated_param['distance'] = model_param['distance']
+
+        return updated_param
+
     def get_spectrum(self,
                      model_param,
                      spec_res):
@@ -109,6 +138,9 @@ def get_spectrum(self,
             Box with the Planck spectrum.
         """
 
+        if 'teff' in model_param and isinstance(model_param['teff'], list):
+            model_param = self.update_parameters(model_param)
+
         wavel_points = [self.wavel_range[0]]
 
         while wavel_points[-1] <= self.wavel_range[1]:
@@ -163,6 +195,9 @@ def get_flux(self,
             Average flux density (W m-2 micron-1).
         """
 
+        if 'teff' in model_param and isinstance(model_param['teff'], list):
+            model_param = self.update_parameters(model_param)
+
         spectrum = self.get_spectrum(model_param, 100.)
 
         if synphot is None:
diff --git a/test/test_read/test_calibration.py b/test/test_read/test_calibration.py
index 7a5aad72..37be193b 100644
--- a/test/test_read/test_calibration.py
+++ b/test/test_read/test_calibration.py
@@ -1,7 +1,7 @@
 import os
 import shutil
-import pytest
 
+import pytest
 import numpy as np
 
 import species
diff --git a/test/test_read/test_color.py b/test/test_read/test_color.py
index a30f9a10..604271d2 100644
--- a/test/test_read/test_color.py
+++ b/test/test_read/test_color.py
@@ -1,7 +1,7 @@
 import os
 import shutil
-import pytest
 
+import pytest
 import numpy as np
 
 import species
diff --git a/test/test_read/test_filter.py b/test/test_read/test_filter.py
index 788e080a..9c406bee 100644
--- a/test/test_read/test_filter.py
+++ b/test/test_read/test_filter.py
@@ -1,7 +1,7 @@
 import os
 import shutil
-import pytest
 
+import pytest
 import numpy as np
 
 import species
diff --git a/test/test_read/test_isochrone.py b/test/test_read/test_isochrone.py
index 619e1e1e..7d5c70bc 100644
--- a/test/test_read/test_isochrone.py
+++ b/test/test_read/test_isochrone.py
@@ -1,8 +1,8 @@
 import os
 import shutil
-import pytest
 import urllib.request
 
+import pytest
 import numpy as np
 
 import species
diff --git a/test/test_read/test_model.py b/test/test_read/test_model.py
index 7f4b3f41..6129b171 100644
--- a/test/test_read/test_model.py
+++ b/test/test_read/test_model.py
@@ -1,7 +1,7 @@
 import os
-import pytest
 import shutil
 
+import pytest
 import numpy as np
 
 import species
diff --git a/test/test_read/test_object.py b/test/test_read/test_object.py
index 0d94a8dd..b10615a3 100644
--- a/test/test_read/test_object.py
+++ b/test/test_read/test_object.py
@@ -1,7 +1,7 @@
 import os
 import shutil
-import pytest
 
+import pytest
 import numpy as np
 
 import species
diff --git a/test/test_read/test_planck.py b/test/test_read/test_planck.py
index d9d0e706..e870c46f 100644
--- a/test/test_read/test_planck.py
+++ b/test/test_read/test_planck.py
@@ -1,7 +1,7 @@
 import os
 import shutil
-import pytest
 
+import pytest
 import numpy as np
 
 import species
diff --git a/test/test_read/test_spectrum.py b/test/test_read/test_spectrum.py
index dcd02d16..cc8b9d9a 100644
--- a/test/test_read/test_spectrum.py
+++ b/test/test_read/test_spectrum.py
@@ -1,7 +1,7 @@
 import os
 import shutil
-import pytest
 
+import pytest
 import numpy as np
 
 import species