This notebook demonstrates how to use the YSO SED models published in Robitaille (2017).
The published models include a tar file for each set of models. The name of each model set is composed of several characters that indicate which component is present. The characters, in order, are:
s (star)p (passive disk)p (power-law envelope) or u (Ulrich envelope)b (bipolar cavities)h (inner hole)m (ambient medium)i (interstellar dust).If a component is absent, a hyphen (-) is given instead.
Each tar file expands to give a directory with the same model set name. The format for each directory is described here.
The easiest way to access and fit these models in Python is to make use of the astropy and sedfitter packages.
Each model directory contains a parameters.fits file that includes the parameters for all the models. To read this, you can use for example the astropy.table package:
In [1]:
from astropy.table import Table
In [2]:
t = Table.read('sp--s-i/parameters.fits')
We can take a look at the first 15 rows of the table:
In [3]:
t[:15]
Out[3]:
The model name is a unique name that identifies each model and the viewing angle is indicated in the suffix (e.g. _01). The value of the inclination is also given in the inclination column. The remaining columns give the parameters for the models (which columns are present depends on the model set). The scattering column indicates whether scattered light is included in the SEDs (for some very optically thick models, scattering was disabled).
The easiest way to access the SEDs in Python is to use the SEDCube class from the sedfitter package to read in the flux.fits file for the model set you are interested in:
In [4]:
from sedfitter.sed import SEDCube
In [5]:
seds = SEDCube.read('sp--s-i/flux.fits')
This 'SED cube' is an efficient way to store the models fluxes in a single 3D array, where the three dimensions are the model, the aperture, and the wavelength.
The model names can be accessed with:
In [6]:
print(seds.names)
while the apertures, wavelengths, and frequencies can be accessed with:
In [7]:
print(seds.apertures)
In [8]:
print(seds.wav)
In [9]:
print(seds.nu)
A valid flag is used to indicate models that do not have complete/valid SEDs (for example because the model run did not complete):
In [10]:
print(seds.valid)
The fluxes and errors can be obtained using the val and
error attributes. We can check the shape of these arrays to check that they are indeed 3D arrays:
In [11]:
seds.val.shape
Out[11]:
In [12]:
seds.val.shape
Out[12]:
For this model set, there are 90000 models (10000 physical models times 9 inclinations), 20 apertures, and 200 wavelengths.
To access a specific SED, you can call seds.get_sed using a particular
model name:
In [13]:
sed = seds.get_sed('00p13Elr_03')
The wavelength, flux, and error can then be accessed with:
In [14]:
print(sed.wav)
In [15]:
print(sed.flux)
In [16]:
print(sed.error)
The SED is a 2D array with dimensions the number of apertures (20) and the number of wavelengths (200):
In [17]:
sed.flux.shape
Out[17]:
We can use this to visualize the SED:
In [18]:
%matplotlib inline
import matplotlib.pyplot as plt
In [19]:
_ = plt.loglog(sed.wav, sed.flux.transpose(), 'k-', alpha=0.5)
_ = plt.ylim(1e-2, 1e8)
To fit SEDs to observed data, you can also make use of the sedfitter package. What follows is a very short example - for more information on using the sedfitter package, be sure to read over the documentation.
To demonstrate this, we will fit the above models to the data for the NGC2264 source modelled in Robitaille (2017):
In [20]:
%cat data_ngc2264_20
We start off by setting up the list of filters/wavelengths and approximate aperture radii used:
In [21]:
from astropy import units as u
In [22]:
filters = ['BU', 'BB', 'BV', 'BR', 'BI', '2J', '2H', '2K', 'I1', 'I2',
5.580 * u.micron, 7.650 * u.micron, 9.95 * u.micron,
12.93 * u.micron, 17.72 * u.micron, 24.28 * u.micron,
29.95 * u.micron, 35.06 * u.micron,
'M2', 'M3', 'W1', 'W2']
In [23]:
apertures = [3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3.,
3., 3., 20., 30., 30., 30.] * u.arcsec
We also set up the extinction law used in Robitaille (2017):
In [24]:
from sedfitter.extinction import Extinction
In [25]:
extinction = Extinction.from_file('whitney.r550.par')
Finally, we run the fitting:
In [26]:
import sedfitter
In [27]:
sedfitter.fit('data_ngc2264_20', filters, apertures, 'sp--s-i',
'output_ngc2264_sp--s-i.fitinfo',
extinction_law=extinction,
distance_range=[0.869, 0.961] * u.kpc,
av_range=[0., 40.],
output_format=('F', 3.),
output_convolved=False, remove_resolved=True)
We now generate the SED plots with the data to examine the fit:
In [28]:
sedfitter.plot('output_ngc2264_sp--s-i.fitinfo',
output_dir='plots_sed', format='png',
plot_mode='A',
select_format=('F', 3.),
show_convolved=False, show_sed=True,
x_mode='M', x_range=(0.1, 2000),
y_mode='M', y_range=(1.e-14, 2e-8))
In [29]:
from IPython.display import Image
Image('plots_sed/20.png')
Out[29]: