DustMix Class Reference

#include <DustMix.hpp>

Inheritance diagram for DustMix:

Public Types
enum class	ScatteringMode { HenyeyGreenstein , MaterialPhaseFunction , SphericalPolarization , SpheroidalPolarization }
Public Types inherited from MaterialMix
enum class	DynamicStateType { None , Primary , Secondary , PrimaryIfMergedIterations }
enum class	MaterialType { Dust , Electrons , Gas }

Public Member Functions
double	asymmpar (double lambda) const override
Array	emissionSpectrum (const MaterialState *state, const Array &Jv) const override
DisjointWavelengthGrid *	emissionWavelengthGrid () const override
Array	emissivity (const Array &Jv) const override
bool	hasContinuumEmission () const override
double	indicativeTemperature (const MaterialState *state, const Array &Jv) const override
double	mass () const override
MaterialType	materialType () const override
double	opacityAbs (double lambda, const MaterialState *state, const PhotonPacket *pp) const override
double	opacityExt (double lambda, const MaterialState *state, const PhotonPacket *pp) const override
double	opacitySca (double lambda, const MaterialState *state, const PhotonPacket *pp) const override
void	peeloffScattering (double &I, double &Q, double &U, double &V, double &lambda, Direction bfkobs, Direction bfky, const MaterialState *state, const PhotonPacket *pp) const override
void	performScattering (double lambda, const MaterialState *state, PhotonPacket *pp) const override
virtual ScatteringMode	scatteringMode () const
double	sectionAbs (double lambda) const override
double	sectionExt (double lambda) const override
const Array &	sectionsAbs (double lambda) const override
const Array &	sectionsAbspol (double lambda) const override
double	sectionSca (double lambda) const override
vector< StateVariable >	specificStateVariableInfo () const override
const Array &	thetaGrid () const override
Public Member Functions inherited from MaterialMix
virtual double	asymmpar (double lambda) const
virtual Array	emissionSpectrum (const MaterialState *state, const Array &Jv) const
virtual DisjointWavelengthGrid *	emissionWavelengthGrid () const
virtual Array	emissivity (const Array &Jv) const
virtual bool	hasContinuumEmission () const
virtual DynamicStateType	hasDynamicMediumState () const
virtual bool	hasExtraSpecificState () const
virtual bool	hasLineEmission () const
virtual bool	hasNegativeExtinction () const
virtual bool	hasPolarizedAbsorption () const
virtual bool	hasPolarizedEmission () const
virtual bool	hasPolarizedScattering () const
virtual bool	hasResonantScattering () const
virtual bool	hasScatteringDispersion () const
virtual bool	hasStochasticDustEmission () const
virtual double	indicativeTemperature (const MaterialState *state, const Array &Jv) const
virtual void	initializeSpecificState (MaterialState *state, double metallicity, double temperature, const Array &params) const
bool	isDust () const
bool	isElectrons () const
bool	isGas () const
virtual bool	isSpecificStateConverged (int numCells, int numUpdated, int numNotConverged, MaterialState *currentAggregate, MaterialState *previousAggregate) const
virtual Array	lineEmissionCenters () const
virtual Array	lineEmissionMasses () const
virtual Array	lineEmissionSpectrum (const MaterialState *state, const Array &Jv) const
virtual double	mass () const =0
virtual MaterialType	materialType () const =0
virtual double	opacityAbs (double lambda, const MaterialState *state, const PhotonPacket *pp) const =0
virtual double	opacityExt (double lambda, const MaterialState *state, const PhotonPacket *pp) const =0
virtual double	opacitySca (double lambda, const MaterialState *state, const PhotonPacket *pp) const =0
virtual vector< SnapshotParameter >	parameterInfo () const
virtual void	peeloffScattering (double &I, double &Q, double &U, double &V, double &lambda, Direction bfkobs, Direction bfky, const MaterialState *state, const PhotonPacket *pp) const =0
virtual void	performScattering (double lambda, const MaterialState *state, PhotonPacket *pp) const =0
virtual double	sectionAbs (double lambda) const =0
virtual double	sectionExt (double lambda) const =0
virtual const Array &	sectionsAbs (double lambda) const
virtual const Array &	sectionsAbspol (double lambda) const
virtual double	sectionSca (double lambda) const =0
virtual vector< StateVariable >	specificStateVariableInfo () const =0
virtual const Array &	thetaGrid () const
virtual UpdateStatus	updateSpecificState (MaterialState *state, const Array &Jv) const
Public Member Functions inherited from SimulationItem
template<class T >
T *	find (bool setup=true) const
template<class T >
T *	interface (int levels=-999999, bool setup=true) const
virtual string	itemName () const
void	setup ()
string	typeAndName () const
Public Member Functions inherited from Item
	Item (const Item &)=delete
virtual	~Item ()
void	addChild (Item *child)
const vector< Item * > &	children () const
virtual void	clearItemListProperty (const PropertyDef *property)
void	destroyChild (Item *child)
virtual bool	getBoolProperty (const PropertyDef *property) const
virtual vector< double >	getDoubleListProperty (const PropertyDef *property) const
virtual double	getDoubleProperty (const PropertyDef *property) const
virtual string	getEnumProperty (const PropertyDef *property) const
virtual int	getIntProperty (const PropertyDef *property) const
virtual vector< Item * >	getItemListProperty (const PropertyDef *property) const
virtual Item *	getItemProperty (const PropertyDef *property) const
virtual string	getStringProperty (const PropertyDef *property) const
int	getUtilityProperty (string name) const
virtual void	insertIntoItemListProperty (const PropertyDef *property, int index, Item *item)
Item &	operator= (const Item &)=delete
Item *	parent () const
virtual void	removeFromItemListProperty (const PropertyDef *property, int index)
virtual void	setBoolProperty (const PropertyDef *property, bool value)
virtual void	setDoubleListProperty (const PropertyDef *property, vector< double > value)
virtual void	setDoubleProperty (const PropertyDef *property, double value)
virtual void	setEnumProperty (const PropertyDef *property, string value)
virtual void	setIntProperty (const PropertyDef *property, int value)
virtual void	setItemProperty (const PropertyDef *property, Item *item)
virtual void	setStringProperty (const PropertyDef *property, string value)
void	setUtilityProperty (string name, int value)
virtual string	type () const

Protected Member Functions
	DustMix ()
virtual double	getOpticalProperties (const Array &lambdav, const Array &thetav, Array &sigmaabsv, Array &sigmascav, Array &asymmparv, Table< 2 > &S11vv, Table< 2 > &S12vv, Table< 2 > &S33vv, Table< 2 > &S34vv, ArrayTable< 2 > &sigmaabsvv, ArrayTable< 2 > &sigmaabspolvv)=0
void	informAvailableWavelengthRange (Range available)
virtual size_t	initializeExtraProperties (const Array &lambdav)
void	setupSelfAfter () override
Protected Member Functions inherited from MaterialMix
	MaterialMix ()
Configuration *	config () const
Random *	random () const
void	setupSelfBefore () override
Protected Member Functions inherited from SimulationItem
	SimulationItem ()
virtual bool	offersInterface (const std::type_info &interfaceTypeInfo) const
virtual void	setupSelfAfter ()
virtual void	setupSelfBefore ()
Protected Member Functions inherited from Item
	Item ()

Private Types
using	BaseType = MaterialMix
using	ItemType = DustMix

Private Member Functions
void	applyMueller (double lambda, double theta, StokesVector *sv) const
std::pair< double, double >	generateAnglesFromPhaseFunction (double lambda, const StokesVector *sv) const
double	generateCosineFromPhaseFunction (double lambda) const
int	indexForLambda (double lambda) const
int	indexForTheta (double theta) const
double	phaseFunctionValue (double lambda, double theta, double phi, const StokesVector *sv) const
double	phaseFunctionValueForCosine (double lambda, double costheta) const

Private Attributes
Array	_asymmparv
EquilibriumDustEmissionCalculator	_calc
Array	_lambdav
double	_mu
Array	_pfnormv
Array	_phi1v
Array	_phicv
Array	_phisv
Array	_phiv
Range	_required
Table< 2 >	_S11vv
Table< 2 >	_S12vv
Table< 2 >	_S33vv
Table< 2 >	_S34vv
ArrayTable< 2 >	_sigmaabspolvv
Array	_sigmaabsv
ArrayTable< 2 >	_sigmaabsvv
Array	_sigmaextv
Array	_sigmascav
Array	_thetav
ArrayTable< 2 >	_thetaXvv

Friends
class	ItemRegistry

Detailed Description

DustMix is the abstract base class for all classes representing the material properties of a dust medium. It implements the complete abstract interface defined by the MaterialMix base class for obtaining optical material properties. Depending on the scattering mode returned by the scatteringMode() function (which must be overridden by subclasses that support a scattering mode other than HenyeyGreenstein), the DustMix setup machinery requests the relevant optical dust properties from the subclass, and caches this information (possibly with some additional precalculated data) for later use. During the simulation, these properties can then be served up without further access to the subclass.

In all cases, the properties served through the MaterialMix base class interface correspond to the properties of a single grain population that is representative for the complete dust mix. For dust mixes described by multiple grain populations, the optical properties should be integrated over the grain size distribution and accumulated across all grain populations. In the context of tracking photon paths through a dusty medium, using these integrated absorption and scattering cross sections and Mueller matrix coefficients is mathematically exact. In other words, for scattering modes MaterialPhaseFunction and SphericalPolarization, the representative grain approach does not involve an approximation. However, the calculation of a representative scattering asymmetry parameter $g$ for use with the Henyey-Greenstein scattering mode does involve a non-exact averaging procedure. Because the Henyey-Greenstein scattering phase function is non-physical to begin with, using a single approximate $g$ value for the complete dust mix is usually considered to be acceptable. For more information on the calculation of representative grain properties in SKIRT, refer to the description of the MultiGrainDustMix class.

The representative grain properties offered by this class are insufficient to accurately calculate dust emission spectra for the dust mixture. This is so because the emission spectrum is a nonlinear function of (among many other things) the grain size, and thus a single grain cannot accurately represent a population with a (potentialy large) range of grain sizes. Again, refer to the description of the MultiGrainDustMix class for more information.

The implementation of polarization by scattering in this class is based on the analysis presented by Peest at al. 2017 (A&A, 601, A92).

Extreme forward (or backward) Henyey-Greenstein scattering

In the X-ray wavelength range, dust grains exhibit extreme forward scattering with asymmetry parameter $g$ values very close to unity. Calculating the value of and sampling from the Henyey-Greenstein (HG) phase function becomes numerically unstable for $|g|>1-{10}^{-6}$. Therefore, the implementation clips any larger $|g|$ values to that limit.

More annoyingly, the bias weights for scattering peel-off photon packets can become very large, causing unacceptably high noise levels in the fluxes recorded by instruments. Consider the HG phase function,

$${\mathrm{\Phi }}_{\mathrm{H}\mathrm{G}}(\cos \theta )=\frac{1-g^2}{(1+g^2-2g\cos \theta {)}^{3/2}}.$$

Given the angle $\theta$ between the direction of the incoming photon packet and a given SKIRT instrument, the bias weight for a peel-off photon packet sent to the instrument is given by $w={\mathrm{\Phi }}_{\mathrm{H}\mathrm{G}}(\cos \theta )$. For $g>0$, the largest bias weights are associated with peel-off photon packets emitted in the original propagation direction ( $\theta =0$). Experiments show that the weigths should be kept under $\approx {10}^3$ to achieve acceptable noise levels with a practical number of photon packets. Incidentally, for $g=0.95$ we find $w\approx 780$, indicating the maximum asymmetry parameter value we could support given this criterion. For $g=1-{10}^{-6}$, we find $w\approx 2\times {10}^{12}$.

To address this situation, we average the peel-off bias weight for a photon packet sent to an instrument over a portion of the phase function, as opposed to taking the value at a given angle. This introduces a certain amount of blur in recorded fluxes, while simultaneously decreasing the level of noise. The effect is reminiscent of that of an actual instrument's point spread function.

We introduce the notation $\mathrm{\Psi }(\alpha ,\beta )$ for the definite integral of the HG phase function between two angles $0\le \alpha \le \beta \le \pi$, and assuming $0<|g|<1$,

$$\mathrm{\Psi }(\alpha ,\beta )\equiv {\int }_{\cos \beta }^{\cos \alpha }{\mathrm{\Phi }}_{\mathrm{H}\mathrm{G}}(\cos \theta )\mathrm{d}\cos \theta =\frac{(1-g^2)}{g}\frac{(t_{\beta }-t_{\alpha })}{t_{\beta }{t}_{\alpha }},$$

with

$$t_x=\sqrt{1+g^2-2g\cos x},x=\alpha ,\beta .$$

We similarly use $\ell (\alpha ,\beta )$ for the length of the interval between the cosines of two angles $0\le \alpha \le \beta \le \pi$,

$$\ell (\alpha ,\beta )\equiv \cos \alpha -\cos \beta .$$

We further denote the half-opening angle of the instrument's solid angle as $\delta$, with $0<\delta \ll \pi /2$. The averaged peel-off bias weight can then be written as,

$$⟨w⟩=\{\begin{array}{ll}{\displaystyle \frac{\mathrm{\Psi }(\theta -\delta ,\theta +\delta )}{\ell (\theta -\delta ,\theta +\delta )}}& \mathrm{f}\mathrm{o}\mathrm{r}\delta \le \theta \le \pi -\delta ,\\ {\displaystyle \frac{\mathrm{\Psi }(0,\delta -\theta )+\mathrm{\Psi }(0,\delta +\theta )}{\ell (0,\delta -\theta )+\ell (0,\delta +\theta )}}& \mathrm{f}\mathrm{o}\mathrm{r}\theta <\delta ,\\ {\displaystyle \frac{\mathrm{\Psi }(\theta -\delta ,\pi )+\mathrm{\Psi }(2\pi -\theta -\delta ,\pi )}{\ell (\theta -\delta ,\pi )+\ell (2\pi -\theta -\delta ,\pi )}}& \mathrm{f}\mathrm{o}\mathrm{r}\theta >\pi -\delta .\end{array}$$

The last two expressions are needed to handle the cases where the averaging interval straddles zero or $\pi$ because the definitions above are not valid for angles outside that range.

Investigation of these expressions shows that a half-opening angle of $\delta =4^{\circ }$ dampens the maximum averaged bias weight for all $g$ to $⟨w⟩\approx 820$. For $g=0.999$, the maximum weight reaches $⟨w⟩\approx 810$ and remains saturated from there on. In conclusion, the implementation uses the regular bias weight for $|g|\le 0.95$ and the averaged bias weight for larger values, with $\delta =4^{\circ }$. As a result, recorded fluxes will be blurred accordingly.

Member Enumeration Documentation

ScatteringMode

enum class DustMix::ScatteringMode

strong

This enumeration lists the scattering modes supported by the dust mix hierarchy.

Constructor & Destructor Documentation

DustMix()

DustMix::DustMix ( )

inlineprotected

Default constructor for abstract Item subclass DustMix : "a dust mix" .

Member Function Documentation

applyMueller()

void DustMix::applyMueller ( double  lambda, double  theta, StokesVector *  sv  ) const

private

This function is used with the SphericalPolarization scattering mode. It applies the Mueller matrix transformation for the specified wavelength $\lambda$ and scattering angle $\theta$ to the given polarization state (which serves as both input and output for the function).

For scattering by spherical grains, the Mueller matrix has only four independent coefficients, namely $S_{11}$, $S_{12}$, $S_{33}$, and $S_{34}$, which depend on both the photon wavelength $\lambda$ and the scattering angle $\theta$. These coefficients are obtained from the tables pre-computed during setup.

asymmpar()

double DustMix::asymmpar ( double  lambda ) const

overridevirtual

This function returns the scattering asymmetry parameter $g_{\lambda }=⟨\cos \theta ⟩$ at wavelength $\lambda$, which is used with the HenyeyGreenstein scattering mode. It retrieves the requested value from a table that was pre-computed during setup.

Reimplemented from MaterialMix.

emissionSpectrum()

Array DustMix::emissionSpectrum ( const MaterialState *  state, const Array &  Jv  ) const

overridevirtual

This function returns the emission spectrum (radiated power per unit of solid angle) in the spatial cell and medium component represented by the specified material state and the receiving material mix when it would be embedded in the specified radiation field. For a dust mix, it returns the result of the emissivity() function multiplied by the hydrogen number density retrieved from the material state.

Reimplemented from MaterialMix.

emissionWavelengthGrid()

DisjointWavelengthGrid * DustMix::emissionWavelengthGrid ( ) const

overridevirtual

This function returns the wavelength grid on which dust emission is discretized, i.e. the wavelength grid returned by the Configuration::dustEmissionWLG() function.

Reimplemented from MaterialMix.

emissivity()

Array DustMix::emissivity ( const Array &  Jv ) const

overridevirtual

This function returns the emissivity spectrum ${\epsilon }_{{\ell }^{\prime }}$ (radiated power per unit of solid angle and per hydrogen atom) of the dust mix (or rather of the representative grain population corresponding to the dust mix) when it would be embedded in a given radiation field, assuming that the dust grains are in local thermal equilibrium. The input and output arrays are discretized on the wavelength grids returned by the Configuration::radiationFieldWLG() and Configuration::dustEmissionWLG() functions, repectively.

The equilibrium emissivity of a representative grain population in an embedding radiation field $J_{\lambda }$ can be calculated as

$${\epsilon }_{\lambda }={\varsigma }_{\lambda ,b}^{\text{abs}}{B}_{\lambda }(T_{\text{eq}})$$

with $\mu$ the total dust mass of the dust mix, ${\varsigma }_{\lambda }^{\text{abs}}$ the absorption cross section of the representative grain, and $T_{\text{eq}}$ the equilibrium temperature of that grain, obtained from the energy balance equation as described for the DustMix::indicativeTemperature() function.

The behavior of this function is undefined if the simulation does not track the radiation field, because in that case setup does not calculate the information on which this function relies.

Reimplemented from MaterialMix.

Reimplemented in MultiGrainDustMix.

generateAnglesFromPhaseFunction()

std::pair< double, double > DustMix::generateAnglesFromPhaseFunction ( double  lambda, const StokesVector *  sv  ) const

private

This function is used with the SphericalPolarization scattering mode. It generates random scattering angles $\theta$ and $\varphi$ sampled from the phase function ${\mathrm{\Phi }}_{\lambda }(\theta ,\varphi )$ at wavelength $\lambda$, and for the specified incoming polarization state. The results are returned as a pair of numbers in the order $\theta$ and $\varphi$.

For scattering by spherical grains, we sample from the phase function listed for the phaseFunctionValue() function using the conditional probability technique. We reduce the phase function to the marginal distribution $\mathrm{\Phi }(\theta )$,

$$\mathrm{\Phi }(\theta )={\int }_0^{2\pi }\mathrm{\Phi }(\theta ,\varphi )\text{d}\varphi =2\pi N{S}_{11}=\frac{S_{11}}{{\int }_0^{\pi }{S}_{11}\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }}.$$

We sample a random $\theta$ value from this distribution through numerical inversion, that is to say, by solving the equation

$$\mathcal{X}=\frac{{\int }_0^{\theta }{S}_{11}\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }}{{\int }_0^{\pi }{S}_{11}\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }}$$

for $\theta$, where $\mathcal{X}$ is a uniform deviate.

Once we have selected a random scattering angle $\theta$, we sample a random azimuthal angle $\varphi$ from the normalized conditional distribution,

$${\mathrm{\Phi }}_{\theta }(\varphi )=\frac{\mathrm{\Phi }(\theta ,\varphi )}{{\int }_0^{2\pi }\mathrm{\Phi }(\theta ,{\varphi }^{\prime })\text{d}{\varphi }^{\prime }}=\frac{1}{2\pi }(1+P_{\text{L}}\frac{S_{12}}{S_{11}}\cos 2(\varphi -\gamma )),$$

where the ratio $S_{12}/S_{11}$ depends on the scattering angle $\theta$ in addition to wavelength.

This can again be done through numerical inversion, by solving the equation

$$\mathcal{X}={\int }_0^{\varphi }{\mathrm{\Phi }}_{\theta }({\varphi }^{\prime })\text{d}{\varphi }^{\prime }=\frac{1}{2\pi }(\varphi +P_{\text{L}}\frac{S_{12}}{S_{11}}\sin \varphi \cos (\varphi -2\gamma ))$$

for $\varphi$, with $\mathcal{X}$ being a new uniform deviate.

generateCosineFromPhaseFunction()

double DustMix::generateCosineFromPhaseFunction ( double  lambda ) const

private

This function is used with the MaterialPhaseFunction scattering mode, which assumes that the scattering phase function depends only on the cosine of the scattering angle. The function generates a random scattering angle cosine sampled from the phase function ${\mathrm{\Phi }}_{\lambda }(\cos \theta )$ at wavelength $\lambda$.

The phase function for unpolarized radiation is obtained from the general polarized case described for the phaseFunctionValue() function by setting the linear polarization degree $P_{\text{L}}$ to zero. The simplified function no longer depends on $\varphi$, and uses only the first Mueller matrix coefficient $S_{11}$. To sample the scattering angle $\theta$ from this phase function, we follow the same procedure as described for the generateAnglesFromPhaseFunction() function, with $P_{\text{L}}=0$.

getOpticalProperties()

virtual double DustMix::getOpticalProperties ( const Array &  lambdav, const Array &  thetav, Array &  sigmaabsv, Array &  sigmascav, Array &  asymmparv, Table< 2 > &  S11vv, Table< 2 > &  S12vv, Table< 2 > &  S33vv, Table< 2 > &  S34vv, ArrayTable< 2 > &  sigmaabsvv, ArrayTable< 2 > &  sigmaabspolvv  )

protectedpure virtual

This function must be implemented in each subclass to obtain the representative grain optical properties for the dust mix. Although the function arguments provide for all properties that may be required, the list of actually required properties for a particular subclass depends on the scattering mode advertised by the scatteringMode() function.

The first two arguments respectively specify the wavelength grid and (if applicable) the scattering angle grid on which the properties must be tabulated. The function must store the absorption and scattering cross sections and the asymmetry parameter into the three subsequent output arrays, and the Mueller coefficients into the four subsequent output tables (indexed on wavelength and angle in that order). These arrays and tables will already have the appropriate size (corresponding to the input wavelength grids) when the function gets called. Finally, the function must return the dust mass per hydrogen atom for the dust mix.

For the HenyeyGreenstein scattering mode, the function must output the cross sections and the asymmetry parameter, and it must leave the Mueller coeficient tables untouched. For scattering modes other than HenyeyGreenstein, the asymmetry parameter output array may be filled with "relevant" values (not used for calculations but possibly exposed to the user through a Probe), or it may be left untouched. For the SphericalPolarization scattering mode, the function must fill all four Mueller coeficient tables. For the MaterialPhaseFunction scattering mode, the function must fill only the first table and it must leave the other tables untouched.

For the SpheroidalPolarization mode, the function must also fill the sigmaabsvv and sigmaabspolvv tables.

Implemented in MultiGrainDustMix, SingleGrainDustMix, and TabulatedDustMix.

hasContinuumEmission()

bool DustMix::hasContinuumEmission ( ) const

overridevirtual

This function returns true for this class because all dust mixes support secondary continuum emission.

Reimplemented from MaterialMix.

indexForLambda()

int DustMix::indexForLambda ( double  lambda ) const

private

This function returns the index in the private wavelength grid corresponding to the specified wavelength. The parameters for converting a wavelength to the appropriate index are stored in data members during setup.

indexForTheta()

int DustMix::indexForTheta ( double  theta ) const

private

This function returns the index in the private scattering angle grid corresponding to the specified scattering angle. The parameters for converting a scattering angle to the appropriate index are built-in constants.

indicativeTemperature()

double DustMix::indicativeTemperature ( const MaterialState *  state, const Array &  Jv  ) const

overridevirtual

This function returns an indicative temperature of the material mix when it would be embedded in a given radiation field. For dust mixes, it returns the equilibrium temperature $T_{\text{eq}}$ of the dust mix (or rather of the representative grain population corresponding to the dust mix) when it would be embedded in the radiation field specified by the mean intensities $(J_{\lambda }{)}_{\ell }$, which must be discretized on the simulation's radiation field wavelength grid as returned by the Configuration::radiationFieldWLG() function. If the specified Jv array is empty (because the simulation does not track the radiation field), this function returns zero.

The equilibrium temperature is obtained from the energy balance equation,

$${\int }_0^{\mathrm{\infty }}{\varsigma }^{\text{abs}}(\lambda )J_{\lambda }(\lambda )\text{d}\lambda ={\int }_0^{\mathrm{\infty }}{\varsigma }^{\text{abs}}(\lambda )B_{\lambda }(T_{\text{eq}},\lambda )\text{d}\lambda ,$$

where the left-hand side is integrated over the radiation field wavelength grid, and the right-hand side is precalculated during setup for a range of temperatures through integration over a built-in wavelength grid.

Reimplemented from MaterialMix.

informAvailableWavelengthRange()

void DustMix::informAvailableWavelengthRange ( Range  available )

protected

This function logs a warning message if the given range is smaller than the simulation wavelength range. It can (but does not have to) be called from subclasses that support dust properties for a limited wavelength range.

initializeExtraProperties()

virtual size_t DustMix::initializeExtraProperties ( const Array &  lambdav )

protectedvirtual

This function can be implemented in a subclass to initialize dust properties that are required to offer additional functionality. The argument specifies the wavelength grid on which the properties may be tabulated (i.e. the same grid as passed to the getOpticalProperties() function.

The function is called by the DustMix class during setup after the getOpticalProperties() function has been called and its results have been processed. The function returns the number of bytes allocated by the subclass to support the extra features (this number is used for logging purposes). The default implementation of this function does nothing and returns zero.

Reimplemented in MultiGrainDustMix.

mass()

double DustMix::mass ( ) const

overridevirtual

This function returns the dust mass $\mu$ per hydrogen atom for this dust mix. It returns a value that was pre-computed during setup.

Implements MaterialMix.

materialType()

MaterialType DustMix::materialType ( ) const

overridevirtual

This function returns the fundamental material type represented by this material mix, in other words it returns MaterialType::Dust.

Implements MaterialMix.

opacityAbs()

double DustMix::opacityAbs ( double  lambda, const MaterialState *  state, const PhotonPacket *  pp  ) const

overridevirtual

This function returns the absorption opacity $k^{\text{abs}}=n{\varsigma }^{\text{abs}}$ for the given wavelength and material state. The photon properties are not used.

Implements MaterialMix.

opacityExt()

double DustMix::opacityExt ( double  lambda, const MaterialState *  state, const PhotonPacket *  pp  ) const

overridevirtual

This function returns the extinction opacity $k^{\text{ext}}=k^{\text{abs}}+k^{\text{sca}}$ for the given wavelength and material state. The photon properties are not used.

Implements MaterialMix.

opacitySca()

double DustMix::opacitySca ( double  lambda, const MaterialState *  state, const PhotonPacket *  pp  ) const

overridevirtual

This function returns the scattering opacity $k^{\text{sca}}=n{\varsigma }^{\text{sca}}$ for the given wavelength and material state. The photon properties are not used.

Implements MaterialMix.

peeloffScattering()

void DustMix::peeloffScattering ( double &  I, double &  Q, double &  U, double &  V, double &  lambda, Direction  bfkobs, Direction  bfky, const MaterialState *  state, const PhotonPacket *  pp  ) const

overridevirtual

This function calculates the contribution of the medium component associated with this material mix to the peel-off photon luminosity, polarization state, and wavelength shift for the given wavelength, geometry, material state, and photon properties. The contributions to the Stokes vector components are stored in the I, Q, U, V arguments, which are guaranteed to be initialized to zero by the caller. For dust mixes, the wavelength remains unchanged.

Evaluation of the phase function depends on the scattering mode supported by supported by the dust mix, as defined by each subclass. For the most basic mode, the material mix provides a value for the scattering asymmetry parameter $g=⟨\cos \theta ⟩$. A value of $g=0$ corresponds to isotropic scattering. Other values $-1\le g\le 1$ are substituted in the Henyey-Greenstein phase function,

$$\mathrm{\Phi }(\cos \theta )=\frac{1-g^2}{(1+g^2-2g\cos \theta {)}^{3/2}}.$$

For other scattering modes, the phase function provided by the material mix is invoked instead.

In case polarization is supported in the current simulation configuration, the polarization state of the peel off photon packet is adjusted as well. The adjusted Stokes vector for a particular medium component is obtained as follows. The function rotates the Stokes vector from the reference direction in the previous scattering plane into the peel-off scattering plane, applies the Mueller matrix on the Stokes vector, and further rotates the Stokes vector from the reference direction in the peel-off scattering plane to the x-axis of the instrument to which the peel-off photon packet is headed.

Implements MaterialMix.

performScattering()

void DustMix::performScattering ( double  lambda, const MaterialState *  state, PhotonPacket *  pp  ) const

overridevirtual

This function performs a scattering event on the specified photon packet in the spatial cell and medium component represented by the specified material state and the receiving material mix.

For dust mixes, the operation depends on the scattering mode supported by the dust mix, as defined by each subclass. For the most basic mode, the dust mix provides a value for the scattering asymmetry parameter $g=⟨\cos \theta ⟩$. For the value $g=0$, corresponding to isotropic scattering, a new direction is generated uniformly on the unit sphere. For other values $-1\le g\le 1$, a scattering angle $\theta$ is sampled from the Henyey-Greenstein phase function,

$$\mathrm{\Phi }(\cos \theta )=\frac{1-g^2}{(1+g^2-2g\cos \theta {)}^{3/2}}.$$

This can be accomplished as follows. Substituting $\mu =\cos \theta$, the probability distribution for $\mu$ (normalized to unity) becomes

$$p(\mu )\text{d}\mu =\frac{1}{2}\frac{1-g^2}{(1+g^2-2g\mu {)}^{3/2}}\text{d}\mu -1\le \mu \le 1$$

We can use the transformation method to sample from this distribution. Given a uniform deviate $\mathcal{X}$, we need to solve

$$\mathcal{X}={\int }_{-1}^{\mu }p({\mu }^{\prime })\text{d}{\mu }^{\prime }$$

Performing the integration and solving for $\mu$ yields

$$\cos \theta =\mu =\frac{1+g^2-f^2}{2g}\text{with}f=\frac{1-g^2}{1-g+2g\mathcal{X}}\text{for}g\ne 0$$

For other scattering modes, a function provided by the dust mix is invoked instead to obtain a random scattering direction for the photon packet.

In case polarization is supported in the current simulation configuration, the polarization state of the photon packet is adjusted as well. Note that all media must either support polarization or not support it, mixing these support levels is not allowed. Compliance with this requirement is verified during setup of the simulation. The adjusted Stokes vector is obtained as follows, again using the randomly selected medium component. After obtaining the sampled scattering angles $\theta$ and $\varphi$ from the material mix, the Stokes vector of the photon packet is rotated into the scattering plane and transformed by applying the Mueller matrix. Finally, the new direction is computed from the previously sampled $\theta$ and $\varphi$ angles.

Implements MaterialMix.

phaseFunctionValue()

double DustMix::phaseFunctionValue ( double  lambda, double  theta, double  phi, const StokesVector *  sv  ) const

private

This function is used with the SphericalPolarization scattering mode. It returns the value of the scattering phase function ${\mathrm{\Phi }}_{\lambda }(\theta ,\varphi )$ at wavelength $\lambda$ for the specified scattering angles $\theta$ and $\varphi$, and for the specified incoming polarization state. The phase function is normalized as

$$\int {\mathrm{\Phi }}_{\lambda }(\theta ,\varphi )\mathrm{d}\mathrm{\Omega }=4\pi .$$

The phase function for scattering by spherical grains can be written as

$$\mathrm{\Phi }(\theta ,\varphi )=N{S}_{11}(1+P_{\text{L}}\frac{S_{12}}{S_{11}}\cos 2(\varphi -\gamma ))$$

with

$$N=\frac{1}{2\pi {\int }_0^{\pi }{S}_{11}\sin \theta \text{d}\theta },$$

where $P_{\text{L}}$ is the linear polarization degree and $\gamma$ the polarization angle of the incoming photon, and where the Mueller matrix coefficients $S_{xx}$ depend on both the photon wavelength $\lambda$ and the scattering angle $\theta$.

phaseFunctionValueForCosine()

double DustMix::phaseFunctionValueForCosine ( double  lambda, double  costheta  ) const

private

This function is used with the MaterialPhaseFunction scattering mode, which assumes that the scattering phase function depends only on the cosine of the scattering angle. The function returns the value of the scattering phase function ${\mathrm{\Phi }}_{\lambda }(\cos \theta )$ at wavelength $\lambda$ for the specified scattering angle cosine $\cos \theta$, where the phase function is normalized as

$${\int }_{-1}^1{\mathrm{\Phi }}_{\lambda }(\cos \theta )\mathrm{d}\cos \theta =2.$$

The phase function for unpolarized radiation is obtained from the general polarized case described for the phaseFunctionValue() function by setting the linear polarization degree $P_{\text{L}}$ to zero. The simplified function uses only the first Mueller matrix coefficient $S_{11}$.

scatteringMode()

virtual ScatteringMode DustMix::scatteringMode ( ) const

virtual

This function returns the scattering mode supported by this dust mix. In the current implementation, this can be one of the following modes:

• HenyeyGreenstein: the value returned by the asymmpar() function serves as the assymmetry parameter $g$ for the Henyey-Greenstein phase function. For a value of $g=0$, the Henyey-Greenstein phase function describes isotropic scattering.
• MaterialPhaseFunction: this material type implements a custom phase function that depends only on the cosine of the scattering angle, for unpolarized radiation. Specifically, the phaseFunctionValueForCosine() and generateCosineFromPhaseFunction() functions are used to obtain the value of the phase function and to sample a scattering angle from it.
• SphericalPolarization: this material type supports polarization through scattering by spherical particles. In this mode, the phase function depends on the polarization state of the incoming radiation, and the polarization state of the outgoing radiation must be updated appropriately. The phaseFunctionValue() and generateAnglesFromPhaseFunction() functions are used to obtain the value of the phase function and to sample a scattering angle from it, and the applyMueller() function is used to updated the polarization state.
• SpheroidalPolarization: this material type supports polarization through scattering, absorption and emission by nonspherical, spheroidal particles. Currently, only emission is implemented and all other areas of the code treat spheroidal particles as if they were spherical.

The implementation of this function in this base class returns the HenyeyGreenstein scattering mode as a default value. Subclasses that support another scattering mode must override this function and return the appropriate value.

Reimplemented in ConfigurableDustMix, DustMixFragment, MeanPinteBenchmarkDustMix, and MeanTrustBenchmarkDustMix.

sectionAbs()

double DustMix::sectionAbs ( double  lambda ) const

overridevirtual

This function returns the absorption cross section per entity ${\varsigma }_{\lambda }^{\text{abs}}$ at wavelength $\lambda$. It retrieves the requested value from a table that was pre-computed during setup.

Implements MaterialMix.

sectionExt()

double DustMix::sectionExt ( double  lambda ) const

overridevirtual

This function returns the total extinction cross section per entity ${\varsigma }_{\lambda }^{\text{ext}}={\varsigma }_{\lambda }^{\text{abs}}+{\varsigma }_{\lambda }^{\text{sca}}$ at wavelength $\lambda$. It retrieves the requested value from a table that was pre-computed during setup.

Implements MaterialMix.

sectionsAbs()

const Array & DustMix::sectionsAbs ( double  lambda ) const

overridevirtual

This function is intended for use with the SpheroidalPolarization mode. It returns the absorption cross sections per entity ${\varsigma }_{\lambda }^{\text{abs}}(\theta )$ at wavelength $\lambda$ as a function of the emission angle $\theta$, discretized on the grid returned by the thetaGrid() function.

Reimplemented from MaterialMix.

sectionsAbspol()

const Array & DustMix::sectionsAbspol ( double  lambda ) const

overridevirtual

This function is intended for use with the SpheroidalPolarization mode. It returns the linear polarization absorption cross sections per entity ${\varsigma }_{\lambda }^{\text{abspol}}(\theta )$ at wavelength $\lambda$ as a function of the emission angle $\theta$, discretized on the grid returned by the thetaGrid() function.

Reimplemented from MaterialMix.

sectionSca()

double DustMix::sectionSca ( double  lambda ) const

overridevirtual

This function returns the scattering cross section per entity ${\varsigma }_{\lambda }^{\text{sca}}$ at wavelength $\lambda$. It retrieves the requested value from a table that was pre-computed during setup.

Implements MaterialMix.

setupSelfAfter()

void DustMix::setupSelfAfter ( )

overrideprotectedvirtual

This function obtains and/or precalculates information used to serve optical properties for the dust mix to the simulation.

First, it obtains the optical properties of the dust mix from the subclass on a wavelength grid with a range spanning all wavelengths in the simulation, and with sufficient resolution so that no further interpolation is needed. This includes the absorption and scattering cross sections, the scattering asymmetry parameter, and/or the Mueller matrix coefficients as required by the implemented scattering mode.

If the simulation wavelength range extends beyond 10 cm, it is cut off at that value and any dust cross sections beyond 10 cm are forced to zero. None of the currently provided dust mixes offers optical properties beyond 10 cm, and historically any values outside the supported range are clamped to the nearest available value. This leads to substantially overestimated dust extinction in the radio wavelength range. Hence this "hack".

Furthermore, if the simulation tracks the radiation field, this function precalculates the Planck-integrated absorption cross sections on an appropriate temperature grid. This information is used to obtain the equilibrium temperature of the material mix (or rather, of its representative grain population) in a given embedding radiation field.

Reimplemented from SimulationItem.

Reimplemented in DustMixFragment.

specificStateVariableInfo()

vector< StateVariable > DustMix::specificStateVariableInfo ( ) const

overridevirtual

This function returns a list of StateVariable objects describing the specific state variables used by the receiving material mix. See the description of the MaterialMix::specificStateVariableInfo() function for more information.

Dust mixes require just the standard specific state variable of type numberDensity , so this function returns a list containing a single item.

Implements MaterialMix.

thetaGrid()

const Array & DustMix::thetaGrid ( ) const

overridevirtual

This function is intended for use with the SpheroidalPolarization mode. It returns the grid used for discretizing quantities that are a function of the scattering/emission angle $\theta$. The same grid is returned by all material mixes that have SpheroidalPolarization mode.

Reimplemented from MaterialMix.

---

The documentation for this class was generated from the following file:

• DustMix.hpp