Source code for dust_attenuation.radiative_transfer

# -*- coding: utf-8 -*-

import numpy as np
import astropy.units as u
import pkg_resources

from astropy.io import ascii
from astropy.modeling.tabular import tabular_model

from .baseclasses import BaseAtttauVModel
from .helpers import _test_valid_x_range


__all__ = ['WG00']

x_range_WG00 = [0.1, 3.0001]


[docs]class WG00(BaseAtttauVModel): """ Attenuation curve of Witt & Gordon (2000) Parameters ---------- tau_v: float optical depth in V band Raises ------ InputParameterError Input Av values outside of defined range Notes ----- From Witt & Gordon (2000, ApJ, Volume 528, pp. 799-816) Example: .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt import astropy.units as u from dust_attenuation.radiative_transfer import WG00 fig, ax = plt.subplots(1,2, figsize=(10,6)) # generate the curves and plot them # Use 1/microns for a better sampling x = np.arange(0.35,10.0,0.1)/u.micron x_Vband = 0.55 # microns tau_Vs = [0.25,0.4,1.1,17.0,46.0] for tau_V in tau_Vs[::-1]: att_model = WG00(tau_V = tau_V, geometry = 'cloudy', dust_type = 'mw', dust_distribution = 'clumpy') ax[0].plot(x,att_model(1/x),label=r'$\tau_V$ = %.2f mag' % (tau_V)) ax[1].plot(x,att_model(1/x)/att_model(x_Vband), label=r'$\tau_V$ = %.2f mag' % (tau_V)) ax[0].set_xlabel('$x$ [$\mu m^{-1}$]') ax[0].set_ylabel(r'$Att(x)$ [mag]') ax[1].set_xlabel('$x$ [$\mu m^{-1}$]') ax[1].set_ylabel(r'$Att(x)/Att_V$') ax[0].legend(loc='best') ax[1].legend(loc='best') fig.suptitle(r'CLOUDY / MW / clumpy model',size=15) plt.tight_layout() fig.subplots_adjust(top=0.88) plt.show() """ tau_V_range = [0.25, 50.0] x_range = x_range_WG00 def __init__(self, tau_V, geometry='dusty', dust_type='mw', dust_distribution='clumpy'): """ Load the attenuation curves for a given geometry, dust type and dust distribution. Parameters ---------- tau_V: float optical depth in V band geometry: string 'shell', 'cloudy' or 'dusty' dust_type: string 'mw' or 'smc' dust_distribution: string 'homogeneous' or 'clumpy' Returns ------- Attx: np array (float) Att(x) attenuation curve [mag] """ # Ensure strings are lower cases self.geometry = geometry.lower() self.dust_type = dust_type.lower() self.dust_distribution = dust_distribution.lower() data_path = pkg_resources.resource_filename('dust_attenuation', 'data/WG00/') data = ascii.read(data_path + self.geometry + '.txt', header_start=0) if self.dust_type == 'mw': start = 0 elif self.dust_type == 'smc': start = 25 # Column names tau_colname = 'tau' tau_att_colname = 'tau_att' fsca_colname = 'f(sca)' fdir_colname = 'f(dir)' fesc_colname = 'f(esc)' if self.dust_distribution == 'clumpy': tau_att_colname += '_c' fsca_colname += '_c' fdir_colname += '_c' fesc_colname += '_c' elif self.dust_distribution == 'homogeneous': tau_att_colname += '_h' fsca_colname += '_h' fdir_colname += '_h' fesc_colname += '_h' tau_att_list = [] tau_list = [] fsca_list = [] fdir_list = [] fesc_list = [] len_data = len(data['lambda']) # number of lines between 2 models steps = 25 counter = start while counter < len_data: tau_att_list.append( np.array(data[tau_att_colname][counter:counter+steps])) tau_list.append( np.array(data[tau_colname][counter:counter+steps])) fsca_list.append( np.array(data[fsca_colname][counter:counter+steps])) fdir_list.append( np.array(data[fdir_colname][counter:counter+steps])) fesc_list.append( np.array(data[fesc_colname][counter:counter+steps])) counter += int(2*steps) # Convert to np.array and take transpose to have (wvl, tau_V) tau_att_table = np.array(tau_att_list).T tau_table = np.array(tau_list).T fsca_table = np.array(fsca_list).T fdir_table = np.array(fdir_list).T fesc_table = np.array(fesc_list).T # wavelength grid. It is the same for all the models wvl = np.array(data['lambda'][0:25]) self.wvl_grid = wvl # Grid for the optical depth tau_V_grid = np.array([0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0]) # Create a 2D tabular model for tau_att and all flux fraction tab = tabular_model(2, name='2D_table') # Values corresponding to the x and y grid points gridpoints = (wvl, tau_V_grid) self.model = tab(gridpoints, lookup_table=tau_att_table, name='tau_att_WG00', bounds_error=False, fill_value=None, method='linear') self.tau = tab(gridpoints, lookup_table=tau_table, name='tau_WG00', bounds_error=False, fill_value=None, method='linear') self.fsca = tab(gridpoints, lookup_table=fsca_table, name='fsca_WG00', bounds_error=False, fill_value=None, method='linear') self.fdir = tab(gridpoints, lookup_table=fdir_table, name='fdir_WG00', bounds_error=False, fill_value=None, method='linear') self.fesc = tab(gridpoints, lookup_table=fesc_table, name='fesc_WG00', bounds_error=False, fill_value=None, method='linear') # In Python 2: super(WG00, self)  # In Python 3: super() but super(WG00, self) still works super(WG00, self).__init__(tau_V=tau_V)
[docs] def evaluate(self, x, tau_V): """ WG00 function Parameters ---------- x: float expects either x in units of wavelengths or frequency or assumes wavelengths in [micron] internally microns are used tau_V: float optical depth in V band Returns ------- Attx: np array (float) Att(x) attenuation curve [mag] Raises ------ ValueError Input x values outside of defined range """ # convert to wavenumbers (1/micron) if x input in units # otherwise, assume x in appropriate wavenumber units with u.add_enabled_equivalencies(u.spectral()): x_quant = u.Quantity(x, u.micron, dtype=np.float64) # strip the quantity to avoid needing to add units to all the # polynomical coefficients x = x_quant.value # check that the wavenumbers are within the defined range _test_valid_x_range(x, x_range_WG00, 'WG00') # setup the ax vectors n_x = len(x) xinterp = 1e4 * x yinterp = tau_V * np.ones(n_x) taux = self.model(xinterp, yinterp) # Convert optical depth to attenuation Attx = 1.086 * taux return Attx
[docs] def get_extinction(self, x, tau_V): """ Return the extinction at a given wavelength and V-band optical depth. Parameters ---------- x: float expects either x in units of wavelengths or frequency or assumes wavelengths in [micron] internally microns are used tau_V: float optical depth in V band Returns ------- ext: np array (float) ext(x) extinction curve [mag] Raises ------ ValueError Input x values outside of defined range """ # convert to wavenumbers (1/micron) if x input in units # otherwise, assume x in appropriate wavenumber units with u.add_enabled_equivalencies(u.spectral()): x_quant = u.Quantity(x, u.micron, dtype=np.float64) # strip the quantity to avoid needing to add units to all the # polynomical coefficients x = x_quant.value # check that the wavenumbers are within the defined range _test_valid_x_range(x, x_range_WG00, 'WG00') # setup the ax vectors x = np.atleast_1d(x) n_x = len(x) xinterp = 1e4 * x yinterp = tau_V * np.ones(n_x) return self.tau(xinterp, yinterp) * 1.086
[docs] def get_fsca(self, x, tau_V): """ Return the scattered flux fraction at a given wavelength and V-band optical depth. Parameters ---------- x: float expects either x in units of wavelengths or frequency or assumes wavelengths in [micron] internally microns are used tau_V: float optical depth in V band Returns ------- fsca: np array (float) fsca(x) scattered flux fraction Raises ------ ValueError Input x values outside of defined range """ # convert to wavenumbers (1/micron) if x input in units # otherwise, assume x in appropriate wavenumber units with u.add_enabled_equivalencies(u.spectral()): x_quant = u.Quantity(x, u.micron, dtype=np.float64) # strip the quantity to avoid needing to add units to all the # polynomical coefficients x = x_quant.value # check that the wavenumbers are within the defined range _test_valid_x_range(x, x_range_WG00, 'WG00') # setup the ax vectors x = np.atleast_1d(x) n_x = len(x) xinterp = 1e4 * x yinterp = tau_V * np.ones(n_x) return self.fsca(xinterp, yinterp)
[docs] def get_fdir(self, x, tau_V): """ Return the direct attenuated stellar flux fraction at a given wavelength and V-band optical depth. Parameters ---------- x: float expects either x in units of wavelengths or frequency or assumes wavelengths in [micron] internally microns are used tau_V: float optical depth in V band Returns ------- fsca: np array (float) fsca(x) scattered flux fraction Raises ------ ValueError Input x values outside of defined range """ # convert to wavenumbers (1/micron) if x input in units # otherwise, assume x in appropriate wavenumber units with u.add_enabled_equivalencies(u.spectral()): x_quant = u.Quantity(x, u.micron, dtype=np.float64) # strip the quantity to avoid needing to add units to all the # polynomical coefficients x = x_quant.value # check that the wavenumbers are within the defined range _test_valid_x_range(x, x_range_WG00, 'WG00') # setup the ax vectors x = np.atleast_1d(x) n_x = len(x) xinterp = 1e4 * x yinterp = tau_V * np.ones(n_x) return self.fdir(xinterp, yinterp)
[docs] def get_fesc(self, x, tau_V): """ Return the total escaping flux fraction at a given wavelength and V-band optical depth. Parameters ---------- x: float expects either x in units of wavelengths or frequency or assumes wavelengths in [micron] internally microns are used tau_V: float optical depth in V band Returns ------- fsca: np array (float) fsca(x) scattered flux fraction Raises ------ ValueError Input x values outside of defined range """ # convert to wavenumbers (1/micron) if x input in units # otherwise, assume x in appropriate wavenumber units with u.add_enabled_equivalencies(u.spectral()): x_quant = u.Quantity(x, u.micron, dtype=np.float64) # strip the quantity to avoid needing to add units to all the # polynomical coefficients x = x_quant.value # check that the wavenumbers are within the defined range _test_valid_x_range(x, x_range_WG00, 'WG00') # setup the ax vectors x = np.atleast_1d(x) n_x = len(x) xinterp = 1e4 * x yinterp = tau_V * np.ones(n_x) return self.fesc(xinterp, yinterp)
[docs] def get_albedo(self, x): """ Return the albedo in function of wavelength for the corresponding dust type (SMC or MW). The albedo gives the probability a photon is scattered from a dust grain. Parameters ---------- x: float expects either x in units of wavelengths or frequency or assumes wavelengths in [micron] internally microns are used Returns ------- albedo: np array (float) alb(x) albedo Raises ------ ValueError Input x values outside of defined range """ # convert to wavenumbers (1/micron) if x input in units # otherwise, assume x in appropriate wavenumber units with u.add_enabled_equivalencies(u.spectral()): x_quant = u.Quantity(x, u.micron, dtype=np.float64) # strip the quantity to avoid needing to add units to all the # polynomical coefficients x = x_quant.value # check that the wavenumbers are within the defined range _test_valid_x_range(x, x_range_WG00, 'WG00') # setup the ax vectors x = np.atleast_1d(x) alb_MW = np.array([0.320, 0.409, 0.481, 0.526, 0.542, 0.536, 0.503, 0.432, 0.371, 0.389, 0.437, 0.470, 0.486, 0.499, 0.506, 0.498, 0.502, 0.491, 0.481, 0.500, 0.473, 0.457, 0.448, 0.424, 0.400]) alb_SMC = np.array([0.400, 0.449, 0.473, 0.494, 0.508, 0.524, 0.529, 0.528, 0.523, 0.520, 0.516, 0.511, 0.505, 0.513, 0.515, 0.498, 0.494, 0.489, 0.484, 0.493, 0.475, 0.465, 0.439, 0.417, 0.400]) if self.dust_type == 'smc': albedo = alb_SMC elif self.dust_type == 'mw': albedo = alb_MW tab = tabular_model(1, name='Tabular1D') alb_fit = tab(self.wvl_grid, lookup_table=albedo, name='albedo', bounds_error=False, fill_value=None, method='linear') xinterp = 1e4 * x return alb_fit(xinterp)
[docs] def get_scattering_phase_function(self, x): """ Return the scattering phase function in function of wavelength for the corresponding dust type (SMC or MW). The scattering phase function gives the angle at which the photon scatters. Parameters ---------- x: float expects either x in units of wavelengths or frequency or assumes wavelengths in [micron] internally microns are used Returns ------- g: np array (float) g(x) scattering phase function Raises ------ ValueError Input x values outside of defined range """ # convert to wavenumbers (1/micron) if x input in units # otherwise, assume x in appropriate wavenumber units with u.add_enabled_equivalencies(u.spectral()): x_quant = u.Quantity(x, u.micron, dtype=np.float64) # strip the quantity to avoid needing to add units to all the # polynomical coefficients x = x_quant.value # check that the wavenumbers are within the defined range _test_valid_x_range(x, x_range_WG00, 'WG00') # setup the ax vectors x = np.atleast_1d(x) g_MW = np.array([0.800, 0.783, 0.767, 0.756, 0.745, 0.736, 0.727, 0.720, 0.712, 0.707, 0.702, 0.697, 0.691, 0.685, 0.678, 0.646, 0.624, 0.597, 0.563, 0.545, 0.533, 0.511, 0.480, 0.445, 0.420]) g_SMC = np.array([0.800, 0.783, 0.767, 0.756, 0.745, 0.736, 0.727, 0.720, 0.712, 0.707, 0.702, 0.697, 0.691, 0.685, 0.678, 0.646, 0.624, 0.597, 0.563, 0.545, 0.533, 0.511, 0.480, 0.445, 0.420]) if self.dust_type == 'smc': g = g_SMC elif self.dust_type == 'mw': g = g_MW tab = tabular_model(1, name='Tabular1D') g_fit = tab(self.wvl_grid, lookup_table=g, name='albedo', bounds_error=False, fill_value=None, method='linear') xinterp = 1e4 * x return g_fit(xinterp)