#!/usr/bin/env python
# coding: utf-8

# # Model spectra and synthetic photometry

# In this tutorial, we will have a look at some spectra of the DRIFT-PHOENIX atmospheric model. The spectra are first downloaded and added to the database. Then we will use the functionalities of [ReadModel](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_model.ReadModel) to extract a spectrum and calculate a photometric flux.

# ## Getting started

# We start by importing the required Python modules.

# In[1]:


import species


# Then we initialize [species](https://species.readthedocs.io/en/latest/species.html) with [SpeciesInit](https://species.readthedocs.io/en/latest/species.core.html#species.core.init.SpeciesInit), which creates a default configuration file and the HDF5 database.

# In[2]:


species.SpeciesInit()


# ## Adding model spectra to the database

# To store the spectra, we first create an instance of [Database](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database).

# In[3]:


database = species.Database()


# Let's check which atmospheric models are available by running the [available_models](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database.available_models) method of the [Database](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database) object.

# In[4]:


_ = database.available_models()


# Next, we will import the model spectra with the [add_model](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database.add_model) method of [Database](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database). This step will automatically download the DRIFT-PHOENIX spectra (R = 2000) to the [data_folder](https://species.readthedocs.io/en/latest/configuration.html). The dowloaded data will then be unpacked and added to the database. We restrict the temperature range to 1300-1700 K, so not all spectra are added to the databse.

# In[5]:


database.add_model(model='drift-phoenix', teff_range=(1300., 1700.))


# Two of the grid points were missing in the original data and have been added with a linear, multidimensional interpolation.

# ## Interpolating the model grid

# We will read the spectra from the database by creating an instance of [ReadModel](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_model.ReadModel).

# In[6]:


read_model = species.ReadModel(model='drift-phoenix', wavel_range=(0.5, 10.))


# Let's see what the grid boundaries are from the spectra that are stored in the database.

# In[7]:


read_model.get_bounds()


# We will now interpolate the grid in the ($T_\mathrm{eff}$, $\log(g)$, $[\mathrm{Fe}/\mathrm{H}]$) space at some specific parameter values, which need to be provided in a dictionary. The radius and distance are optional, otherwise the emitted flux is given at the planet surface.

# In[8]:


model_param = {'teff':1510., 'logg':4.1, 'feh':0.1, 'radius': 1., 'distance': 100.}


# To interpolate a spectrum, we use the [get_model](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_model.ReadModel.get_model) method and provide the parameter dictionary, and also an optional spectral resolution. Together with `smooth=True`, the spectrum will be smoothed (but not resampeld) to the given spectral resolution.

# In[9]:


model_box = read_model.get_model(model_param=model_param, spec_res=200., smooth=True)


# The data is stored in a [ModelBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ModelBox). Let's have a look at its content.

# In[10]:


model_box.open_box()


# We will now use the same atmospheric parameters but we will add some [ISM extinction](https://species.readthedocs.io/en/latest/species.util.html?#species.util.dust_util.ism_extinction) with the relation from [Cardelli et al. (1989)](https://ui.adsabs.harvard.edu/abs/1989ApJ...345..245C/abstract). Therefore, we add the V band extinction and reddening parameters to the dictionary.

# In[11]:


model_param = {'teff':1510., 'logg':4.1, 'feh':0.1, 'radius': 1., 'distance': 100., 'ism_ext': 5., 'ism_red': 3.}


# We use again the [get_model](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_model.ReadModel.get_model) method and store the result in a different [ModelBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ModelBox).

# In[12]:


model_ext = read_model.get_model(model_param=model_param, spec_res=200., smooth=True)


# The two boxes with the model spectra are provided to the [plot_spectrum](https://species.readthedocs.io/en/latest/species.plot.html#species.plot.plot_spectrum.plot_spectrum). We also include some filter profiles to indicate where the telluric windows are.

# In[13]:


species.plot_spectrum(boxes=[model_box, model_ext],
                      filters=['MKO/NSFCam.J', 'MKO/NSFCam.H', 'MKO/NSFCam.K', 'MKO/NSFCam.Lp', 'MKO/NSFCam.Mp'],
                      offset=(-0.15, -0.04),
                      xlim=(0.8, 5.),
                      ylim=(0., 5.5e-17),
                      legend={'loc': 'lower right', 'frameon': False, 'fontsize': 10.},
                      figsize=(7., 3.),
                      output=None)


# ## Extracting a spectrum at a grid point

# It is also possible to extract a spectrum at one of the grid points, which doesn't require any interpolation. Let's check with the [get_points](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_model.ReadModel.get_points) method what parameter values are stored in the database.

# In[13]:


read_model.get_points()


# We create a dictionary with values at one of the grid points.

# In[14]:


model_param = {'teff':1500., 'logg':4., 'feh':0.}


# And now use the [get_data](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_model.ReadModel.get_data) method to extract a spectrum.

# In[15]:


model_full = read_model.get_data(model_param)


# Let's make another plot with [plot_spectrum](https://species.readthedocs.io/en/latest/species.plot.html#species.plot.plot_spectrum.plot_spectrum). The spectrum is now shown at the spectral resolution as stored in the database ($R = 2000$).

# In[17]:


species.plot_spectrum(boxes=[model_full],
                      filters=None,
                      offset=(-0.12, -0.05),
                      xlim=(0.8, 5.),
                      ylim=(0., 1e5),
                      legend={'loc': 'upper right', 'fontsize': 10.},
                      figsize=(7., 3.),
                      output=None)


# ## Calculating synthetic photometry

# The [ReadModel](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_model.ReadModel) class can also be used for calculating photometric fluxes and magnitudes. To do so, we create a new instance and set the `filter_name` argument to the [VLT/NACO M' filter](http://svo2.cab.inta-csic.es/svo/theory/fps/index.php?id=Paranal/NACO.Mp&&mode=browse&gname=Paranal&gname2=NACO#filter). This will automatically downloadd and addd the filter profile. 

# In[17]:


read_model = species.ReadModel(model='drift-phoenix', filter_name='Paranal/NACO.Mp')


# We create again a dictionary with the parameters but now run the [get_flux](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_model.ReadModel.get_flux) method, which returns the flux in W m-2 um-1.

# In[18]:


model_param = {'teff':1510., 'logg':4.1, 'feh':0.1, 'radius': 1., 'distance': 100.}
flux = read_model.get_flux(model_param)
print(f'Flux (W m-2 um-1) = {flux[0]:.2e}')


# Since we provided a radius and distance, the emitted flux at the planet surface has been scaled by $(\mathrm{radius}/\mathrm{distance})^2$.

# Similarly, we can use the [get_magnitude](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_model.ReadModel.get_magnitude) method to calculate the magnitude for the NACO M' filter. Note that the returned absolute magnitude is set to `None` if the parameter dictionary does not contain a radius and distance.

# In[19]:


app_mag, abs_mag = read_model.get_magnitude(model_param)
print(f'Apparent magnitude = {app_mag:.2f}')
print(f'Absolute magnitude = {abs_mag:.2f}')


# As expected, at a distance of 100 pc, the difference between the apparent and absolute magnitude is 5.