from astropy.table import Table

t = Table.read('sp--s-i/parameters.fits')

"""
Explanation: Using the Robitaille (2017) YSO SED models
This notebook demonstrates how to use the YSO SED models published in Robitaille (2017).
Format
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.
Accessing Parameters
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:
End of explanation
"""


t[:15]

"""
Explanation: We can take a look at the first 15 rows of the table:
End of explanation
"""


from sedfitter.sed import SEDCube

seds = SEDCube.read('sp--s-i/flux.fits')

"""
Explanation: 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).
Accessing SEDs
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:
End of explanation
"""


print(seds.names)

"""
Explanation: 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:
End of explanation
"""


print(seds.apertures)

print(seds.wav)

print(seds.nu)

"""
Explanation: while the apertures, wavelengths, and frequencies can be accessed with:
End of explanation
"""


print(seds.valid)

"""
Explanation: A valid flag is used to indicate models that do not have complete/valid SEDs (for example because the model run did not complete):
End of explanation
"""


seds.val.shape

seds.val.shape

"""
Explanation: 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:
End of explanation
"""


sed = seds.get_sed('00p13Elr_03')

"""
Explanation: 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:
End of explanation
"""


print(sed.wav)

print(sed.flux)

print(sed.error)

"""
Explanation: The wavelength, flux, and error can then be accessed with:
End of explanation
"""


sed.flux.shape

"""
Explanation: The SED is a 2D array with dimensions the number of apertures (20) and the number of wavelengths (200):
End of explanation
"""


%matplotlib inline
import matplotlib.pyplot as plt

_ = plt.loglog(sed.wav, sed.flux.transpose(), 'k-', alpha=0.5)
_ = plt.ylim(1e-2, 1e8)

"""
Explanation: We can use this to visualize the SED:
End of explanation
"""


%cat data_ngc2264_20

"""
Explanation: Fitting SEDs to data
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):
End of explanation
"""


from astropy import units as u

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']

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

"""
Explanation: We start off by setting up  the list of filters/wavelengths and approximate aperture radii used:
End of explanation
"""


from sedfitter.extinction import Extinction

extinction = Extinction.from_file('whitney.r550.par')

"""
Explanation: We also set up the extinction law used in Robitaille (2017):
End of explanation
"""


import sedfitter

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)

"""
Explanation: Finally, we run the fitting:
End of explanation
"""


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))

from IPython.display import Image
Image('plots_sed/20.png')

"""
Explanation: We now generate the SED plots with the data to examine the fit:
End of explanation
"""


