Fitting data with a grid of model spectra#

In this tutorial, we will fit spectra and photometric fluxes of beta Pic b with the synthetic spectra from a radiative-convective equilibrium atmosphere model. Here we will use DRIFT-PHOENIX but there are also several other models supported by species.

Getting started#

We start by setting the library path of MultiNest (see installation instructions), which is required for the Bayesian inference.

import os
os.environ['DYLD_LIBRARY_PATH'] = '/Users/tomasstolker/applications/MultiNest/lib'

We also need to import urllib and the species toolkit.

import species
import urllib.request

We initiate species by running the SpeciesInit class. By doing so, both the configuration file and the HDF5 database are created in the working folder.

Initiating species v0.5.1... [DONE]
Creating species_config.ini... [DONE]
Database: /Users/tomasstolker/applications/species/docs/tutorials/species_database.hdf5
Data folder: /Users/tomasstolker/applications/species/docs/tutorials/data
Working folder: /Users/tomasstolker/applications/species/docs/tutorials
Creating species_database.hdf5... [DONE]
Creating data folder... [DONE]
<species.core.init.SpeciesInit at 0x1867b1480>

Let’s now download the GRAVITY \(K\) band spectrum of beta Pic b that was published by Gravity Collaboration et al. 2020 and we will also make use of the GPI \(YJH\) band spectra from Chilcote et al. 2017.

('betapicb_gpi_h.dat', <http.client.HTTPMessage at 0x1867b19c0>)

Adding model spectra to the database#

Data will be added to the database by first creating an instance of Database.

database = species.Database()

The Database object has a number of methods to read and write data. For adding model spectra, we use the add_model method which will download the grid of DRIFT-PHOENIX spectra to the data_folder (see configuration file) except if the data was already present. After downloading, it will import the spectra into the database. We will only require a limited range of \(T_\mathrm{eff}\) since we know the approximate temperature of beta Pic b.

database.add_model(model='drift-phoenix', teff_range=(1500., 2000.))
Downloading DRIFT-PHOENIX model spectra (229 MB)... [DONE]
Unpacking DRIFT-PHOENIX model spectra (229 MB)... [DONE]
Please cite Helling et al. (2008) when using DRIFT-PHOENIX in a publication
Reference URL:
Wavelength range (um) = 0.1 - 50
Spectral resolution = 4000
Teff range (K) = 1500.0 - 2000.0
Adding DRIFT-PHOENIX model spectra... [DONE]
Grid points stored in the database:
   - Teff = [1500. 1600. 1700. 1800. 1900. 2000.]
   - log(g) = [3.  3.5 4.  4.5 5.  5.5]
   - [Fe/H] = [-0.6 -0.3 -0.   0.3]

Adding photometry and spectra#

We will also add photometric data of beta Pic b that have been adopted in species from various references. To get an overview of all photometric data that were stored in the database, we can use the list_companions method of the Database object (i.e. database.list_companions() in this case).

To add the photometric data and parallax of beta Pic b to the database, we use the add_companion method, which will also convert the magnitudes into fluxes with a flux-calibrated spectrum of Vega and the filter profiles (downloaded from the SVO Filter Profile Service). It is also possible to manually add photometric data with the app_mag and flux_density arguments of add_object.

database.add_companion(name='beta Pic b')
Getting GPI_YJHK spectrum of beta Pic b... [DONE]
IMPORTANT: Please cite Chilcote et al. 2017, AJ, 153, 182
           when making use of this spectrum in a publication
Getting GRAVITY spectrum of beta Pic b... [DONE]
IMPORTANT: Please cite Gravity Collaboration et al. 2020, A&A, 633, 110
           when making use of this spectrum in a publication
Downloading Vega spectrum (270 kB)... [DONE]
Adding Vega spectrum... [DONE]
Adding filter: Magellan/VisAO.Ys... [DONE]
Adding filter: Paranal/NACO.J... [DONE]
Adding filter: Gemini/NICI.ED286... [DONE]
Adding filter: Paranal/NACO.H... [DONE]
Adding filter: Paranal/NACO.Ks... [DONE]
Adding filter: Paranal/NACO.NB374... [DONE]
Adding filter: Paranal/NACO.Lp... [DONE]
Adding filter: Paranal/NACO.NB405... [DONE]
Adding filter: Paranal/NACO.Mp... [DONE]
Adding object: beta Pic b
   - Parallax (mas) = 51.44 +/- 0.12
   - Magellan/VisAO.Ys:
      - Apparent magnitude = 15.53 +/- 0.34
      - Flux (W m-2 um-1) = 4.25e-15 +/- 1.35e-15
   - Paranal/NACO.J:
      - Apparent magnitude = 14.11 +/- 0.21
      - Flux (W m-2 um-1) = 6.87e-15 +/- 1.34e-15
   - Gemini/NICI.ED286:
      - Apparent magnitude = 13.18 +/- 0.15
      - Flux (W m-2 um-1) = 6.99e-15 +/- 9.69e-16
   - Paranal/NACO.H:
      - Apparent magnitude = 13.32 +/- 0.14
      - Flux (W m-2 um-1) = 5.47e-15 +/- 7.08e-16
   - Paranal/NACO.Ks:
      - Apparent magnitude = 12.64 +/- 0.11
      - Flux (W m-2 um-1) = 4.04e-15 +/- 4.10e-16
   - Paranal/NACO.NB374:
      - Apparent magnitude = 11.25 +/- 0.23
      - Flux (W m-2 um-1) = 1.69e-15 +/- 3.61e-16
   - Paranal/NACO.Lp:
      - Apparent magnitude = 11.30 +/- 0.06
      - Flux (W m-2 um-1) = 1.59e-15 +/- 8.79e-17
   - Paranal/NACO.NB405:
      - Apparent magnitude = 10.98 +/- 0.05
      - Flux (W m-2 um-1) = 1.61e-15 +/- 7.42e-17
   - Paranal/NACO.Mp:
      - Apparent magnitude = 11.10 +/- 0.12
      - Flux (W m-2 um-1) = 7.86e-16 +/- 8.70e-17
   - Spectrum:
      - Database tag: GPI_YJHK
      - Filename: ./data/companion_data/betapicb_gpi_yjhk.dat
      - Data shape: (152, 3)
      - Wavelength range (um): 0.98 - 2.37
      - Mean flux (W m-2 um-1): 5.40e-15
      - Mean error (W m-2 um-1): 6.44e-16
   - GRAVITY spectrum:
      - Object: Unknown
      - Database tag: GRAVITY
      - Filename: ./data/companion_data/BetaPictorisb_2018-09-22.fits
      - Data shape: (237, 3)
      - Wavelength range (um): 1.97 - 2.49
      - Mean flux (W m-2 um-1): 4.65e-15
      - Mean error (W m-2 um-1): 1.00e-16
   - GRAVITY covariance matrix:
      - Object: Unknown
      - Database tag: GRAVITY
      - Filename: ./data/companion_data/BetaPictorisb_2018-09-22.fits
      - Data shape: (237, 237)
   - Spectral resolution:
      - GPI_YJHK: 40.0
      - GRAVITY: 500.0
/Users/tomasstolker/applications/species/species/data/ UserWarning: Found 33 fluxes with NaN in the data of GPI_YJHK. Removing the spectral fluxes that contain a NaN.

Next, we will add the spectra with the add_object method. Make sure to use the same database name as argument of object_name (i.e. beta Pic b in this case) such that the spectra are added to the object name for which previously photometric data had been stored with add_companion.

The argument of spectrum is a dictionary with the names of the spectra as key and a tuple with three elements as value. The first element is a text or FITS file with the spectrum, the second element is the optional covariance matrix, and the third value is the mandatory spectral resolution. The latter is used to smooth the model spectra to the instrument resolution.

For this example, we have adopted the VLTI/GRAVITY \(K\) band spectrum from Gravity Collaboration et al. 2020 and the Gemini/GPI \(YJH\) band spectra from Chilcote et al. 2017. Data that have been stored in the FITS format by the ExoGRAVITY collaboration is automatically recognized so the same FITS filename can be provided as first and second element in the dictionary.

database.add_object('beta Pic b',
                    spectrum={'GRAVITY': ('BetaPictorisb_2018-09-22.fits', 'BetaPictorisb_2018-09-22.fits', 500.),
                              'GPI_Y': ('betapicb_gpi_y.dat', None, 40.),
                              'GPI_J': ('betapicb_gpi_j.dat', None, 40.),
                              'GPI_H': ('betapicb_gpi_h.dat', None, 40.)},
Adding object: beta Pic b
   - GRAVITY spectrum:
      - Object: Unknown
      - Database tag: GRAVITY
      - Filename: BetaPictorisb_2018-09-22.fits
      - Data shape: (237, 3)
      - Wavelength range (um): 1.97 - 2.49
      - Mean flux (W m-2 um-1): 4.65e-15
      - Mean error (W m-2 um-1): 1.00e-16
   - Spectrum:
      - Database tag: GPI_Y
      - Filename: betapicb_gpi_y.dat
      - Data shape: (29, 3)
      - Wavelength range (um): 0.98 - 1.13
      - Mean flux (W m-2 um-1): 5.40e-15
      - Mean error (W m-2 um-1): 1.87e-15
   - Spectrum:
      - Database tag: GPI_J
      - Filename: betapicb_gpi_j.dat
      - Data shape: (32, 3)
      - Wavelength range (um): 1.13 - 1.34
      - Mean flux (W m-2 um-1): 6.84e-15
      - Mean error (W m-2 um-1): 4.79e-16
   - Spectrum:
      - Database tag: GPI_H
      - Filename: betapicb_gpi_h.dat
      - Data shape: (34, 3)
      - Wavelength range (um): 1.51 - 1.79
      - Mean flux (W m-2 um-1): 5.63e-15
      - Mean error (W m-2 um-1): 3.86e-16
   - GRAVITY covariance matrix:
      - Object: Unknown
      - Database tag: GRAVITY
      - Filename: BetaPictorisb_2018-09-22.fits
      - Data shape: (237, 237)
   - Spectral resolution:
      - GRAVITY: 500.0
      - GPI_Y: 40.0
      - GPI_J: 40.0
      - GPI_H: 40.0

The use of add_object is incremental so when rerunning the method it will add new data in case a filter or spectrum name hadn’t been used before, or it will overwrite existing data in case a filter or spectrum name does already exist. To remove existing data of an object, it is easiest to remove the group of data entirely from the database with the delete_data method of Database (e.g. delete_data('objects/beta Pic b') and rerun add_companion and/or add_object afterwards.

Bayesian inference with nested sampling#

To fit the data with the grid of model spectra (i.e. grid retrieval), we start by creating an instance of FitModel, which requires the database tag of the object (beta Pic b) and the atmosphere model (drift-phoenix).

The dictionary of the bounds parameter contains the prior boundaries for the various parameters. By default, the parameters of the atmosphere model are set to the full range of the available spectra (e.g. \(\log(g)\) and \(\mathrm{[Fe/H]}\) in the example below).

In addition to the parameters of the atmosphere model, several optional parameters can be fitted, for example to account for calibration systematics, extinction, and excess emission from a disk. More details on the various parameters can be found in the API documentation of FitModel. As example, we fit a scaling parameter for the GPI \(H\) band spectrum.

The arguments of inc_phot and inc_spec are either a boolean or a list of photometry and spectra that were previously stored in the database for the object beta Pic b. We will use all the spectra but only the photometric data in the \(L\) and \(M\) bands.

To account for correlated noise, it is possible to apply the approach from Wang et al. (2020) to estimate the covariances by using a Gaussian process. This will add two additional parameters per spectrum (correlation length and amplitude) but is not mandatory. For the example, we only select the GPI \(H\) band spectrum, to limit the computation time.

For simplicity we create a list with the 4 filter names that we want to use.

inc_phot = ['Paranal/NACO.Lp', 'Paranal/NACO.NB374', 'Paranal/NACO.NB405', 'Paranal/NACO.Mp']

We now create the instance of FitModel. By doing so, the grid of model spectra will be interpolated to the wavelengths of the data.

fit = species.FitModel(object_name='beta Pic b',
                       bounds={'teff': (1500., 2000.),
                               'radius': (0.5, 2.),
                               'GPI_H': ((0.5, 1.5), None)},
Getting object: beta Pic b... [DONE]
Interpolating Paranal/NACO.Lp... [DONE]
Interpolating Paranal/NACO.NB374... [DONE]
Interpolating Paranal/NACO.NB405... [DONE]
Interpolating Paranal/NACO.Mp... [DONE]
Interpolating GPI_H... [DONE]
Interpolating GPI_J... [DONE]
Interpolating GPI_Y... [DONE]
Interpolating GRAVITY... [DONE]
Fitting 6 parameters:
   - teff
   - logg
   - feh
   - radius
   - parallax
   - scaling_GPI_H
Prior boundaries:
   - teff = (1500.0, 2000.0)
   - radius = (0.5, 2.0)
   - logg = (3.0, 5.5)
   - feh = (-0.6, 0.3)
   - scaling_GPI_H = (0.5, 1.5)
Weights for the log-likelihood function:
   - GPI_H = 1.00e+00
   - GPI_J = 1.00e+00
   - GPI_Y = 1.00e+00
   - GRAVITY = 1.00e+00
   - Paranal/NACO.Lp = 1.00e+00
   - Paranal/NACO.NB374 = 1.00e+00
   - Paranal/NACO.NB405 = 1.00e+00
   - Paranal/NACO.Mp = 1.00e+00

We are now ready to sample the posterior distribution by either using MultiNest or UltraNest with the run_multinest) and run_ultranest) methods of FitModel. Both are nested sampling algorithms which are powerful in sampling multi-modal distributions and will estimate the marginalized likelihood (i.e. model evidence), which enables pair-wise model comparison through the Bayes factor.

In this example, we will use the run_multinest method since UltraNest is somewhat more straightforward to install than MultiNest. We specify the database tag were the posterior samples will be stored, the number of live points, the output folder where MultiNest will write its data, and we use the mass of beta Pic b from Gravity Collaboration et al. 2020 as Gaussian prior on \(\log(g)\).

To speed up the computation, it is possible to run the nested sampling in parallel (e.g. with mpirun) to benefit from the multiprocessing support by UltraNest and MultiNest. In that case it is important that any functions of species that will write to the Database will be commented out since simultaneous writing to the HDF5 database by different processes is not possible. It is therefore recommended to first add all the required data to the database and then only run SpeciesInit, FitModel, and the sampler (run_multinest or run_ultranest) in parallel with MPI.

                  prior={'mass': (9., 1.6)})
Running nested sampling with MultiNest...
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    6
/Users/tomasstolker/.pyenv/versions/3.10.0/envs/species3.10/lib/python3.10/site-packages/pymultinest/ DeprecationWarning: inspect.getargspec() is deprecated since Python 3.0, use inspect.signature() or inspect.getfullargspec()
  nargs = len(inspect.getargspec(LogLikelihood).args) - inspect.ismethod(LogLikelihood)
  analysing data from multinest/.txt ln(ev)=   11453.765410571104      +/-  0.18739327712282511
 Total Likelihood Evaluations:        45127
 Sampling finished. Exiting MultiNest

Nested sampling global log-evidence: 11453.77 +/- 0.19
Nested importance sampling global log-evidence: 11452.05 +/- 0.01
Sample with the highest likelihood:
   - Log-likelihood = 11473.79
   - teff = 1717.68
   - logg = 3.95
   - feh = 0.13
   - radius = 1.49
   - parallax = 51.34
   - scaling_GPI_H = 1.12
Integrated autocorrelation time:
   - teff: 3.74
   - logg: 4.48
   - feh: 9.94
   - radius: 3.89
   - parallax: 9.36
   - scaling_GPI_H: 1.64

Plotting the posterior samples#

The samples from the parameter estimation have been stored in the database. We can now run the plot_posterior function to plot the 1D and 2D projections of the posterior distributions by making use of Here we specify the database tag with the results from run_ultranest and we also include the derived posterior for the bolometric luminosity and mass.

                       offset=(-0.3 , -0.3),
                       title_fmt=['.0f', '.2f', '.2f', '.2f', '.2f', '.2f', '.3f', '.1f'],
Median sample:
   - teff = 1.72e+03
   - logg = 3.97e+00
   - feh = 1.29e-01
   - radius = 1.48e+00
   - parallax = 5.14e+01
   - scaling_GPI_H = 1.13e+00
Plotting the posterior... [DONE]

Extracting boxes with results#

We will now extract Box objects with the results. These will be combined in a plot at the end of the tutorial. We start by selecting 30 random samples from the posterior and use these parameters to interpolate the DRIFT-PHOENIX grid. To do so, we use the get_mcmc_spectra method of the Database. We will smooth the spectra to a resolution of \(R=500\) to match with the resolution of the GRAVITY spectrum.

samples = database.get_mcmc_spectra(tag='betapic',
Getting MCMC spectra: 100%|██████████████████████████████████████████| 30/30 [00:00<00:00, 153.62it/s]

The function returned a list with 30 ModelBox objects. Let’s have a look at the content of the first box by running the open_box method. This method can be used for all boxes to easily inspect the attributes of a Box object.

Opening ModelBox...
model = drift-phoenix
type = None
wavelength = [ 0.1         0.10002476  0.10004953 ... 49.97524871 49.98762282
 50.        ]
flux = [2.28135989e-69 2.10383420e-69 1.87528020e-69 ... 2.07169466e-19
 2.06913400e-19 2.06730648e-19]
parameters = {'teff': 1721.6267865756652, 'logg': 3.972139346124936, 'feh': 0.12856649608628445, 'radius': 1.4783947036643863, 'parallax': 51.499197452238754, 'luminosity': 0.00017469381970180441, 'mass': 7.908494362248091}
quantity = flux
contribution = None
bol_flux = None

Instead of using the plotting functionalities of species, the attributes from a Box can also be directly extracted for further processing. For example, to print the wavelengths of the first ModelBox that was returned by get_mcmc_spectra:

[ 0.1         0.10002476  0.10004953 ... 49.97524871 49.98762282
 50.        ]

Next, we get the median parameter values with the get_median_sample method of the Database and adopt these as best-fit parameters. Similarly, the get_probable_sample method can be used for extracting the sample with the maximum likelihood.

best = database.get_median_sample(tag='betapic')

The best-fit model spectrum is now extracted from the DRIFT-PHOENIX grid by using the functionalities of the ReadModel class.

read_model = species.ReadModel(model='drift-phoenix',

The grid is interpolated at the best-fit parameters with get_model. The argument of model_param contains the parameter dictionary that was returned by get_median_sample. The spectrum will be smoothed to a resolution of \(R=500\). To extract the spectrum without smoothing at the resolution as stored in the database is achieved by setting spec_res=None and smooth=False.

modelbox = read_model.get_model(model_param=best,
/Users/tomasstolker/applications/species/species/read/ UserWarning: The 'scaling_GPI_H' parameter is not required by 'drift-phoenix' so the parameter will be ignored. The mandatory parameters are ['teff', 'logg', 'feh'].

Some warnings are printed because the model_param dictionary contained calibration parameters (that were fitted), but these are not required when reading a model spectrum.

We will also extract all spectra and photometric fluxes of beta Pic b with the get_object method which will return the output in an ObjectBox.

objectbox = database.get_object(object_name='beta Pic b',
Getting object: beta Pic b... [DONE]

We then run the ObjectBox through the update_spectra function with the best-fit parameters. This function applies the flux scaling and could also inflate the errors in case the model_param dictionary contains error inflation parameters.

objectbox = species.update_spectra(objectbox=objectbox,
Scaling the flux of GPI_H: 1.13... [DONE]
/Users/tomasstolker/applications/species/species/util/ DeprecationWarning: The update_spectra function is deprecated and will be removed in a future release. Please use the update_objectbox function instead.

Next, we calculate the residuals of the best-fit model with get_residuals. This function will interpolate the grid at the specified parameters, calculate synthetic photometry for the filters that are provided as argument of inc_phot, and smooth and resample the spectra to the resolution and wavelengths of the spectra that are provided as argument of inc_spec (all spectra of the ObjectBox in this case). The returned residuals will be stored in a ResidualsBox.

residuals = species.get_residuals(datatype='model',
Calculating synthetic photometry... [DONE]
Calculating residuals... [DONE]
Residuals (sigma):
   - Paranal/NACO.Lp: -4.75
   - Paranal/NACO.NB374: -1.12
   - Paranal/NACO.NB405: -0.76
   - Paranal/NACO.Mp: -0.81
   - GPI_H: min: -1.17, max: 2.07
   - GPI_J: min: -1.91, max: 3.99
   - GPI_Y: min: -1.36, max: 1.37
   - GRAVITY: min: -5.69, max: 5.11
Reduced chi2 = 2.79
Number of degrees of freedom = 330

Finally, we compute synthetic photometry from the best-fit model spectrum with the multi_photometry function. Since we will plot all available photometric data of the ObjectBox, we will calculate synthetic photometry for all the names of the filters attribute of the ObjectBox.

synphot = species.multi_photometry(datatype='model',
Calculating synthetic photometry... [DONE]

Let’s have a look inside the SynphotBox.

Opening SynphotBox...
name = synphot
flux = {'Gemini/NICI.ED286': 6.313100131660149e-15, 'Magellan/VisAO.Ys': 3.261934978977361e-15, 'Paranal/NACO.H': 6.023674177528963e-15, 'Paranal/NACO.J': 6.382730078161801e-15, 'Paranal/NACO.Ks': 4.9832619510599846e-15, 'Paranal/NACO.Lp': 2.007715791664406e-15, 'Paranal/NACO.Mp': 8.564256078585903e-16, 'Paranal/NACO.NB374': 2.0982081194013644e-15, 'Paranal/NACO.NB405': 1.6663501715060354e-15}

Plotting the spectral energy distribution#

Now that we have gather all the Box objects with the results, we can pass them as list to the boxes parameter of the plot_spectrum function. The ResidualsBox is separately provided as argument of residuals. We also include a list of filter names for which the transmission profiles will be plotted. The arguments of residuals and filters can also be set to None to not include these data in the plot.

The somewhat complex part of plot_spectrum is the optional plot_kwargs parameter. The argument is a list with the same length as the list of boxes. Each item contains a dictionary with keyword arguments that can be used to fine-tune the styling of the plot. The item of the SynphotBox can be set to None because it is automatically adjusted to the styles that are used for the ObjectBox.

species.plot_spectrum(boxes=[samples, modelbox, objectbox, synphot],
                      plot_kwargs=[{'ls': '-', 'lw': 0.2, 'color': 'gray'},
                                   {'ls': '-', 'lw': 1., 'color': 'black'},
                                   {'GRAVITY': {'marker': 'o', 'ms': 5., 'mew': 0., 'color': 'seagreen', 'ls': 'none', 'alpha': 0.5, 'label': 'VLTI/GRAVITY'},
                                    'GPI_Y': {'marker': 'o', 'ms': 5., 'mew': 0., 'color': 'darkslateblue', 'ls': 'none', 'alpha': 0.5, 'label': 'Gemini/GPI'},
                                    'GPI_J': {'marker': 'o', 'ms': 5., 'mew': 0., 'color': 'darkslateblue', 'ls': 'none', 'alpha': 0.5},
                                    'GPI_H': {'marker': 'o', 'ms': 5., 'mew': 0., 'color': 'darkslateblue', 'ls': 'none', 'alpha': 0.5},
                                    'Magellan/VisAO.Ys': {'marker': 's', 'ms': 5., 'color': 'goldenrod', 'ls': 'none', 'label': 'Magellan/VisAO'},
                                    'Gemini/NICI.ED286': {'marker': 's', 'ms': 5., 'color': 'teal', 'ls': 'none', 'label': 'Gemini/NICI'},
                                    'Paranal/NACO.J': {'marker': 's', 'ms': 5., 'color': 'tomato', 'ls': 'none', 'label': 'VLT/NACO'},
                                    'Paranal/NACO.H': {'marker': 's', 'ms': 5., 'color': 'tomato', 'ls': 'none'},
                                    'Paranal/NACO.Ks': {'marker': 's', 'ms': 5., 'color': 'tomato', 'ls': 'none'},
                                    'Paranal/NACO.NB374': {'marker': 's', 'ms': 5., 'color': 'tomato', 'ls': 'none'},
                                    'Paranal/NACO.Lp': {'marker': 's', 'ms': 5., 'color': 'tomato', 'ls': 'none'},
                                    'Paranal/NACO.NB405': {'marker': 's', 'markersize': 5., 'color': 'tomato', 'ls': 'none'},
                                    'Paranal/NACO.Mp': {'marker': 's', 'markersize': 5., 'color': 'tomato', 'ls': 'none'}},
                      xlim=(0.8, 5.2),
                      ylim=(-1.15e-15, 1.15e-14),
                      ylim_res=(-7., 7.),
                      scale=('linear', 'linear'),
                      offset=(-0.4, -0.05),
                      legend=[{'loc': 'lower left', 'frameon': False, 'fontsize': 11.},
                              {'loc': 'upper right', 'frameon': False, 'fontsize': 12.}],
                      figsize=(8., 4.),
                      quantity='flux density',
Plotting spectrum... [DONE]