#include <DipolePhaseFunction.hpp>
Public Member Functions  
void  initialize (Random *random, bool includePolarization=false) 
void  peeloffScattering (double &I, double &Q, double &U, double &V, Direction bfk, Direction bfkobs, Direction bfky, const StokesVector *sv) const 
Direction  performScattering (Direction bfk, StokesVector *sv) const 
Private Member Functions  
void  applyMueller (double theta, StokesVector *sv) const 
std::pair< double, double >  generateAnglesFromPhaseFunction (const StokesVector *sv) const 
double  generateCosineFromPhaseFunction () const 
double  phaseFunctionValue (double theta, double phi, const StokesVector *sv) const 
double  phaseFunctionValueForCosine (double costheta) const 
Private Attributes  
bool  _includePolarization 
Array  _phi1v 
Array  _phicv 
Array  _phisv 
Array  _phiv 
Random *  _random 
The DipolePhaseFunction helper class represents the dipole scattering phase function, with optional support for polarization by scattering. For example, in a wide wavelength range, scattering of photons by electrons can be described by elastic and wavelengthindependent Thomson scattering, which is in turn described using the dipole phase function. The corresponding Mueller matrix elements (and hence the phase function) can be expressed analytically as a function of just the scattering angle; see Bohren & Huffman (1998) or Wolf 2003 (Computer Physics Communications, 150, 99–115).
The implementation of this class is based on the analysis presented by Peest at al. 2017 (A&A, 601, A92).

private 
This function applies the Mueller matrix transformation for the specified scattering angle \(\theta\) to the given polarization state (which serves as both input and output for the function).
For scattering by spherical particles, the Mueller matrix has only four independent coefficients. For dipole scattering, these can be analytically expressed as
\[S_{11}=\cos^2\theta+1,\quad S_{12}=\cos^2\theta1,\quad S_{33}=2\cos\theta,\quad \text{ and } S_{34}=0.\]

private 
This function generates random scattering angles \(\theta\) and \(\phi\) sampled from the phase function \(\Phi(\theta,\phi)\) for the specified incoming polarization state. The results are returned as a pair of numbers in the order \(\theta\) and \(\phi\).
For dipole scattering, 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 by integrating over \(\phi\), yielding
\[ \Phi(\theta) \propto \cos^2\theta+1. \]
With the substitution \(t = \cos\theta\), the normalized probability distribution to be sampled can be written as
\[ p(t) = \frac{3}{8}\,(t^2+1),\text{ with } 1<t<1. \]
The corresponding cumulative distribution is
\[ P(t) = \int_{1}^t \frac{3}{8}\,(t'^2+1)\,\text{d}t' = \frac{1}{8}\,(t^3+3t+4). \]
The equation \({\cal{X}} = P(t)\), with \({\cal{X}}\) a uniform deviate, has a single realvalued solution,
\[ t = p  \frac{1}{p} \quad\text{with}\quad p = \left( 4{\cal{X}}2 +\sqrt{16{\cal{X}}^2 16{\cal{X}} +5} \right)^{1/3}\]
Once we have selected a random scattering angle \(\theta=\arccos(t)\), we sample a random azimuthal angle \(\phi\) from the normalized conditional distribution,
\[ \Phi_\theta(\phi) =\frac{\Phi(\theta,\phi)}{\int_0^{2\pi} \Phi(\theta,\phi')\,\text{d}\phi'} =\frac{1}{2\pi}\left(1+ P_{\text{L}}\,\frac{\cos^2\theta1}{\cos^2\theta+1}\cos 2(\phi  \gamma)\right). \]
This can be done by numerically solving the equation
\[ {\cal{X}} =\int_{0}^{\phi}\Phi_{\theta}(\phi')\,\text{d}\phi' =\frac{1}{2\pi} \left( \phi + P_{\text{L}}\,\frac{\cos^2\theta1}{\cos^2\theta+1} \sin\phi \cos(\phi  2\gamma)\right) \]
for \(\phi\), with \({\cal{X}}\) being a new uniform deviate.

private 
This function generates a random scattering angle cosine sampled from the phase function \(\Phi(\cos\theta)\).
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, which yields
\[ \Phi(\theta) \propto \cos^2\theta+1. \]
We can thus use the same procedure for sampling \(\theta\) as described for the phaseFunctionValue() function.
void DipolePhaseFunction::initialize  (  Random *  random, 
bool  includePolarization = false 

) 
This function caches a pointer to the simulation's random number generator and, if polarization is included, precalculates some data used for sampling from the corresponding phase function. The function must be called during setup (i.e. in singlethread mode). If includePolarization is omitted or set to false, calling the functions implementing the polarized phase function will cause undefined behavior.
void DipolePhaseFunction::peeloffScattering  (  double &  I, 
double &  Q,  
double &  U,  
double &  V,  
Direction  bfk,  
Direction  bfkobs,  
Direction  bfky,  
const StokesVector *  sv  
)  const 
This function calculates the contribution of a dipole scattering event to the peeloff photon luminosity and polarization state for the given geometry and polarization state. 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.
Direction DipolePhaseFunction::performScattering  (  Direction  bfk, 
StokesVector *  sv  
)  const 
Given the incoming photon packet direction and polarization state, this function calculates and returns a randomly sampled new propagation direction for a dipole scattering event, and if applicable (depending in the polarization flag passed to the constructor), updates the polarization state of the photon packet along the way.

private 
This function returns the value of the scattering phase function \(\Phi(\theta,\phi)\) for the specified scattering angles \(\theta\) and \(\phi\), and for the specified incoming polarization state. The phase function is normalized as
\[\int\Phi(\theta,\phi) \,\mathrm{d}\Omega =4\pi.\]
For dipole scattering, we can substitute the analytical expressions for the Mueller coefficients into the phase function, which leads to
\[ \Phi(\theta,\phi) = \frac{3}{4}\,\left[(\cos^2\theta+1) + P_{\text{L}}\,(\cos^2\theta1)\cos2(\phi  \gamma) \right], \]
where \(P_\text{L}\) is the linear polarization degree and \(\gamma\) the polarization angle of the incoming photon.

private 
This function returns the value of the scattering phase function \(\Phi(\cos\theta)\) for the specified scattering angle cosine \(\cos\theta\), where the phase function is normalized as
\[\int_{1}^1 \Phi(\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, which yields simply
\[ \Phi(\theta) = \frac{3}{4}\,(\cos^2\theta+1). \]