Examples

iris Examples

---

Table of Contents

1. 1 Component H2O
2. 2 Components H2O
3. 2 Components H2O and CO2 with Thermal Broadening
4. 2 Components H2O with Keplerian Broadening
5. H2O with Power-law Temperature Gradient

---

1 Component H2O

Here we model a simple 1-component H2O slab, with excitation temperature of 800 K, 0.5 au $^2$ emitting area, and column density of $10^{18}$ cm $^{-2}$.

path_to_moldata = './'  # path where we want to save the HITRAN data

# Let's get the data for H2O 
setup.setup_linelists('H2O', 'H2O', 1, path_to_moldata)

# Define a fine wavelength grid (in micron) to evaluate opacities
fine_wgrid = np.arange(11.5,18.5,1e-5)
# Define a wavelength grid (in micron) to downsample the model
obs_wgrid = np.arange(12.0,18.0,0.002)
# The instrument resolving power
R = 3200

# Initialize slab
slab = iris.slab(molecules=['H2O'], wlow=11.0, whigh=19.0, path_to_moldata=path_to_moldata)

# Distance to source
distance = 120 # pc

# NOTE that we define an ARRAY for the temperature, column density, and area FOR EACH SPECIES.
# So if we wanted to add multiple components for a species just add more points to each array.

# Excitation temperatures for each molecule in K
T_ex =  np.array([np.array([800.0])])
# column densities in cm^-2
N_mol = np.array([np.array([1e18])])
# emitting areas in au^2
A_au =  np.array([np.array([0.5])])

# intrinsic line widths in km/s (line FWHM)
dV =  np.array([np.array([2.0])])

def compiled_slab(distance, T_ex, N_mol, A_au, dV, fine_wgrid, wavelength, R):

    slab.setup_disk(distance, T_ex, N_mol, A_au, dV) # initialize object
    slab.setup_grid(fine_wgrid, wavelength, R) # set up wavelength and model parameters
    slab.simulate() # make model

    return slab.downsampled_flux, slab.flux_model

'''Just-in-time compilation'''
compiled_slab_jit = jax.jit(compiled_slab)
downsampled_model, full_model = compiled_slab_jit(distance, T_ex, N_mol, A_au, dV, fine_wgrid, obs_wgrid, R)

2 Components H2O

For this we can re-use most of the code from Example 1, we just need to change the parameter matrices to model 2 components. In fact, we do not need to re-initialize the slab object or re-compile the code at all, since we are still using the same slab object.

We will model 2 components with temperatures 1000 and 500 K, $N_{mol}$ of $10^{18}$ and $10^{17}$ cm $^{-2}$, and emitting areas of 0.5 and 5.0 au $^{2}$

Note that in principle we could define different line widths for each component.

# Excitation temperatures for each molecule in K
T_ex =  np.array([np.array([1000.0, 500.0])])
# column densities in cm^-2
N_mol = np.array([np.array([1e18, 1e17])])
# emitting areas in au^2
A_au =  np.array([np.array([0.5, 5.0])])

# intrinsic line widths in km/s (line FWHM)
dV =  np.array([np.array([2.0, 2.0])])

downsampled_model, full_model = compiled_slab_jit(distance, T_ex, N_mol, A_au, dV, fine_wgrid, obs_wgrid, R)

2 Components H2O and CO2 with Thermal Broadening

We will model the lines more carefully now, specifying a thermal and turbulent component for the line widths. To start, we'll need some physical constants from astropy:

import astropy.units as u
import astropy.constants as const

# Boltzmann constant
kB_cgs = const.k_B.cgs.value
# Proton mass
mp_cgs = (1.0*u.M_p).cgs.value

path_to_moldata = './'  # path where we want to save the HITRAN data

# Let's get the data for H2O and CO2
setup.setup_linelists('H2O', 'H2O', 1, path_to_moldata)
setup.setup_linelists('CO2', 'CO2', 1, path_to_moldata)

Now we define the turbulent line standard deviation as 2 km/s, and calculate the thermal widths from the excitation temperatures:

$\sigma_{therm} = \sqrt{\frac{k_B T_{ex}}{\mu_X m_p}}$

Here $\mu_X$ is the molecular weight of each species, in this case 18 and 44 for H2O and CO2, respectively.

The total line standard deviation is the sum in quadrature:

$\sigma = \sqrt{\sigma_{therm}^2 + \sigma_{turb}^2}$

Remember we need to provide iris with a line FWHM, so

$\Delta V = 2 \sqrt{2 \log 2} \sigma$

# Define a fine wavelength grid (in micron) to evaluate opacities
fine_wgrid = np.arange(11.5,18.5,1e-5)
# Define a wavelength grid (in micron) to downsample the model
obs_wgrid = np.arange(12.0,18.0,0.002)
# The instrument resolving power
R = 3200

# Initialize slab
slab = iris.slab(molecules=['H2O', 'CO2'], wlow=11.0, whigh=19.0, path_to_moldata=path_to_moldata)

# Distance to source
distance = 120 # pc

# Molecular weights for each species
molecular_weights = np.array([18.0, 44.0])

# Excitation temperatures for each molecule in K
T_ex =  np.array([np.array([1000.0, 500.0]),
                  np.array([800.0, 400.0])])
# column densities in cm^-2
N_mol = np.array([np.array([1e18, 5e17]),
                  np.array([1e16, 5e15])])
# emitting areas in au^2
A_au =  np.array([np.array([0.5, 5.0]),
                  np.array([1.0, 10.0])])

# Reshape so that molecular_weights has the same dimensions as T_ex
molecular_weights = molecular_weights.reshape(2,1) * np.ones_like(T_ex)

# TURBULENT line widths in km/s (line standard deviation)
dV =  np.array([np.array([2.0, 2.0])])

# THERMAL line widths in km/s (line standard deviation)
thermal_dV = ((kB_cgs*(T_ex)/(molecular_weights*mp_cgs))**0.5) / 100000

# Add TURBULENT + THERMAL in quadrature, calculate FWHM
dV = 2.355 * ((thermal_dV)**2 + (dV)**2)**0.5

def compiled_slab(distance, T_ex, N_mol, A_au, dV, fine_wgrid, wavelength, R):

    slab.setup_disk(distance, T_ex, N_mol, A_au, dV) # initialize object
    slab.setup_grid(fine_wgrid, wavelength, R) # set up wavelength and model parameters
    slab.simulate() # make model

    return slab.downsampled_flux, slab.flux_model

'''Just-in-time compilation'''
compiled_slab_jit = jax.jit(compiled_slab)
downsampled_model, full_model = compiled_slab_jit(distance, T_ex, N_mol, A_au, dV, fine_wgrid, obs_wgrid, R)

2 Components H2O with Keplerian Broadening

H2O with Power-law Temperature Gradient

H2O with Varying Resolving Power

Bayesian Retrieval of H2O