Examples
- 1 Component H2O
- 2 Components H2O
- 2 Components H2O and CO2 with Thermal Broadening
- 2 Components H2O and CO2 with Keplerian Broadening
Here we model a simple 1-component H2O slab, with excitation temperature of 800 K, 0.5 au
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)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,
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)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.valuepath_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:
Here
The total line standard deviation is the sum in quadrature:
Remember we need to provide iris with a line FWHM, so
# 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)In this Example, we'll do something similar as in Example 3 but we will also include the effects of Keplerian line broadening.
We start by setting up the slab with thermal and turbulent broadening only:
# 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])])
downsampled_model, full_model = compiled_slab_jit(distance, T_ex, N_mol, A_au, dV, fine_wgrid, obs_wgrid, R)
# 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.5Now, we need to specify the disk inclination, stellar mass, and innermost radius for each molecule. iris will automatically calculate the outer radii from the emitting areas provided, even in models with multiple components.
# Disk Inclination in degree
inc = 45.0
# Stellar Mass in M solar
M_star = 1.0
# Inner Radii in au
r_in = np.array([0.01, 0.1]) # 0.01 au for H2O, 0.1 au for CO2Now we can make the slab model with Keplerian line profiles using slab.simulate_keplerian():
def compiled_slab_keplerian(distance, T_ex, N_mol, A_au, dV, inc, M_star, r_in, fine_wgrid, wavelength, R):
slab.setup_disk(distance, T_ex=T_ex, N_mol=N_mol, A_au=A_au, dV=dV, inc=inc, M_star=M_star, r_in=r_in) # initialize object
slab.setup_grid(fine_wgrid, wavelength, R) # set up wavelength and model parameters
slab.simulate_keplerian() # make KEPLERIAN model
return slab.downsampled_flux, slab.flux_model
'''Just-in-time compilation'''
compiled_slab_keplerian_jit = jax.jit(compiled_slab_keplerian)
downsampled_model, full_model = compiled_slab_keplerian_jit(distance, T_ex, N_mol, A_au, dV, inc, M_star, r_in, fine_wgrid, obs_wgrid, R)Below we can see a comparison between the same (full resolution) slab model with and without Keplerian broadening:
