add XRayIonicGasMix

SKIRT/CMakeLists.txt

@@ -8,6 +8,7 @@
 # ------------------------------------------------------------------
 
 # add all relevant subdirectories; each subdirectory defines a single target
+add_subdirectory(faddeeva)
 add_subdirectory(fitsio)
 add_subdirectory(voro)
 add_subdirectory(tetgen)

SKIRT/core/CMakeLists.txt

@@ -22,8 +22,8 @@ target_link_libraries(${TARGET} schema fundamentals)
 include_directories(../../SMILE/schema ../../SMILE/fundamentals)
 
 # add SKIRT library dependencies
-target_link_libraries(${TARGET} fitsio voro mpi utils tetgen)
-include_directories(../fitsio ../voro ../mpi ../utils)
+target_link_libraries(${TARGET} fitsio voro mpi utils tetgen faddeeva)
+include_directories(../fitsio ../voro ../mpi ../utils ../faddeeva)
 include_directories(SYSTEM ../tetgen)  # suppress warnings in tetgen header
 
 # adjust C++ compiler flags to our needs

SKIRT/core/LyaUtils.cpp → SKIRT/core/LyUtils.cpp

@@ -3,7 +3,7 @@
 ////       © Astronomical Observatory, Ghent University         ////
 ///////////////////////////////////////////////////////////////// */
 
-#include "LyaUtils.hpp"
+#include "LyUtils.hpp"
 #include "Configuration.hpp"
 #include "Constants.hpp"
 #include "Random.hpp"
@@ -13,32 +13,24 @@
 
 namespace
 {
-    constexpr double c = Constants::c();              // speed of light in vacuum
-    constexpr double kB = Constants::k();             // Boltzmann constant
-    constexpr double mp = Constants::Mproton();       // proton mass
-    constexpr double la = Constants::lambdaLya();     // central Lyman-alpha wavelength
-    constexpr double Aa = Constants::EinsteinALya();  // Einstein A coefficient for Lyman-alpha transition
+    constexpr double c = Constants::c();  // speed of light in vacuum
 }
 
 ////////////////////////////////////////////////////////////////////
 
-double LyaUtils::section(double lambda, double T)
+double LyUtils::section(double lambda, double center, double vth, double a, double g)
 {
-    double vth = sqrt(2. * kB / mp * T);                 // thermal velocity for T
-    double a = Aa * la / 4. / M_PI / vth;                // Voigt parameter
-    double x = (la - lambda) / lambda * c / vth;         // dimensionless frequency
-    double sigma0 = 3. * la * la * M_2_SQRTPI / 4. * a;  // cross section at line center
-    return sigma0 * VoigtProfile::value(a, x);           // cross section at given x
+    double x = (center - lambda) / lambda * c / vth;            // dimensionless frequency
+    double sigma0 = g * center * center * M_2_SQRTPI / 4. * a;  // cross section at line center
+    return sigma0 * VoigtProfile::value(a, x);                  // cross section at given x
 }
 
 ////////////////////////////////////////////////////////////////////
 
-std::pair<Vec, bool> LyaUtils::sampleAtomVelocity(double lambda, double T, double nH, Direction kin,
-                                                  Configuration* config, Random* random)
+Vec LyUtils::sampleAtomVelocity(double lambda, double center, double vth, double a, double T, double nH, Direction kin,
+                                Configuration* config, Random* random)
 {
-    double vth = sqrt(2. * kB / mp * T);          // thermal velocity for T
-    double a = Aa * la / 4. / M_PI / vth;         // Voigt parameter
-    double x = (la - lambda) / lambda * c / vth;  // dimensionless frequency
+    double x = (center - lambda) / lambda * c / vth;  // dimensionless frequency
 
     // generate two directions that are orthogonal to each other and to the incoming photon packet direction
     Direction k1 = (kin.kx() || kin.ky()) ? Direction(kin.ky(), -kin.kx(), 0., true) : Direction(1., 0., 0., false);
@@ -79,21 +71,16 @@ std::pair<Vec, bool> LyaUtils::sampleAtomVelocity(double lambda, double T, doubl
     // transform the dimensionless frequency into the rest frame of the atom
     x -= Vec::dot(u, kin);
 
-    // select the isotropic or the dipole phase function:
-    // all wing events and 1/3 of core events are dipole, and the remaining 2/3 core events are isotropic,
-    // where x=0.2 (in the atom frame) defines the transition between core and wings
-    bool dipole = abs(x) > 0.2 || random->uniform() < 1. / 3.;
-
     // scale the atom velocity from dimensionless to regular units
     u *= vth;
 
     // return the atom velocity and the phase function choice
-    return std::make_pair(u, dipole);
+    return u;
 }
 
 ////////////////////////////////////////////////////////////////////
 
-double LyaUtils::shiftWavelength(double lambda, const Vec& vatom, const Direction& kin, const Direction& kout)
+double LyUtils::shiftWavelength(double lambda, const Vec& vatom, const Direction& kin, const Direction& kout)
 {
     return lambda / (1 - Vec::dot(kin, vatom) / c) * (1 - Vec::dot(kout, vatom) / c);
 }

SKIRT/core/LyaUtils.hpp → SKIRT/core/LyUtils.hpp

@@ -3,11 +3,11 @@
 ////       © Astronomical Observatory, Ghent University         ////
 ///////////////////////////////////////////////////////////////// */
 
-#ifndef LYAUTILS_HPP
-#define LYAUTILS_HPP
+#ifndef LYUTILS_HPP
+#define LYUTILS_HPP
 
 #include "Direction.hpp"
-#include "Range.hpp"
+
 class Configuration;
 class Random;
 
@@ -102,12 +102,12 @@ class Random;
     purpose of this recipe, the scattering event is considered to occur in the core if the incoming
     dimensionless photon frequency in the rest frame of the interacting atom is smaller than a
     critical value, \f$|x|<0.2\f$. */
-namespace LyaUtils
+namespace LyUtils
 {
     /** This function returns the Lyman-alpha scattering cross section per hydrogen atom
         \f$\sigma_\alpha(\lambda, T)\f$ at the given photon wavelength and gas temperature, using
         the definition given in the class header. */
-    double section(double lambda, double T);
+    double section(double lambda, double center, double vth, double a, double g);
 
     /** This function draws a random hydrogen atom velocity as seen by an incoming photon from the
         appropriate probability distributions, reflecting the preference for photons to be
@@ -141,8 +141,8 @@ namespace LyaUtils
         - Return the atom velocity and a flag indicating the selected phase function.
 
         */
-    std::pair<Vec, bool> sampleAtomVelocity(double lambda, double T, double nH, Direction kin, Configuration* config,
-                                            Random* random);
+    Vec sampleAtomVelocity(double lambda, double center, double vth, double a, double T, double nH, Direction kin,
+                           Configuration* config, Random* random);
 
     /** This function returns the Doppler-shifted wavelength in the gas bulk rest frame after a
         Lyman-alpha scattering event, given the incoming wavelength in the gas bulk rest frame, the

SKIRT/core/LyaNeutralHydrogenGasMix.cpp

@@ -4,15 +4,26 @@
 ///////////////////////////////////////////////////////////////// */
 
 #include "LyaNeutralHydrogenGasMix.hpp"
-#include "Configuration.hpp"
 #include "Constants.hpp"
-#include "LyaUtils.hpp"
+#include "LyUtils.hpp"
 #include "MaterialState.hpp"
 #include "PhotonPacket.hpp"
 #include "Random.hpp"
 
 ////////////////////////////////////////////////////////////////////
 
+namespace
+{
+    // the combined Lya1 and Lya2 for hydrogen
+    constexpr double lyaA = Constants::EinsteinALya();
+    constexpr double lya = Constants::lambdaLya();
+    constexpr double g = 3.;
+    constexpr double kB = Constants::k();
+    constexpr double mp = Constants::Mproton();
+}
+
+////////////////////////////////////////////////////////////////////
+
 void LyaNeutralHydrogenGasMix::setupSelfBefore()
 {
     MaterialMix::setupSelfBefore();
@@ -87,6 +98,15 @@ double LyaNeutralHydrogenGasMix::mass() const
 
 ////////////////////////////////////////////////////////////////////
 
+double LyaNeutralHydrogenGasMix::section(double lambda, double T) const
+{
+    double vth = sqrt(2. * kB / mp * T);
+    double a = lyaA * lya / 4. / M_PI / vth;
+    return LyUtils::section(lambda, lya, vth, a, g);
+}
+
+////////////////////////////////////////////////////////////////////
+
 double LyaNeutralHydrogenGasMix::sectionAbs(double /*lambda*/) const
 {
     return 0.;
@@ -96,14 +116,14 @@ double LyaNeutralHydrogenGasMix::sectionAbs(double /*lambda*/) const
 
 double LyaNeutralHydrogenGasMix::sectionSca(double lambda) const
 {
-    return LyaUtils::section(lambda, defaultTemperature());
+    return section(lambda, defaultTemperature());
 }
 
 ////////////////////////////////////////////////////////////////////
 
 double LyaNeutralHydrogenGasMix::sectionExt(double lambda) const
 {
-    return LyaUtils::section(lambda, defaultTemperature());
+    return section(lambda, defaultTemperature());
 }
 
 ////////////////////////////////////////////////////////////////////
@@ -119,31 +139,52 @@ double LyaNeutralHydrogenGasMix::opacityAbs(double /*lambda*/, const MaterialSta
 double LyaNeutralHydrogenGasMix::opacitySca(double lambda, const MaterialState* state, const PhotonPacket* /*pp*/) const
 {
     double n = state->numberDensity();
-    return n > 0. ? n * LyaUtils::section(lambda, state->temperature()) : 0.;
+    return n > 0. ? n * section(lambda, state->temperature()) : 0.;
 }
 
 ////////////////////////////////////////////////////////////////////
 
 double LyaNeutralHydrogenGasMix::opacityExt(double lambda, const MaterialState* state, const PhotonPacket* /*pp*/) const
 {
     double n = state->numberDensity();
-    return n > 0. ? n * LyaUtils::section(lambda, state->temperature()) : 0.;
+    return n > 0. ? n * section(lambda, state->temperature()) : 0.;
 }
 
 ////////////////////////////////////////////////////////////////////
 
-bool LyaNeutralHydrogenGasMix::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda,
-                                                 Direction bfkobs, Direction bfky, const MaterialState* state,
-                                                 const PhotonPacket* pp) const
+void LyaNeutralHydrogenGasMix::setScatteringInfoIfNeeded(PhotonPacket* pp, const MaterialState* state,
+                                                         const double lambda) const
 {
-    // draw a random atom velocity & phase function, unless a previous peel-off stored this already
-    auto scatinfo = const_cast<PhotonPacket*>(pp)->getScatteringInfo();
+    auto scatinfo = pp->getScatteringInfo();
     if (!scatinfo->valid)
     {
         scatinfo->valid = true;
-        std::tie(scatinfo->velocity, scatinfo->dipole) = LyaUtils::sampleAtomVelocity(
-            lambda, state->temperature(), state->numberDensity(), pp->direction(), config(), random());
+
+        double T = state->temperature();
+        double nH = state->numberDensity();
+
+        double vth = sqrt(2. * kB / mp * T);
+        double a = lyaA * lya / 4. / M_PI / vth;
+        double x = (lya - lambda) / lambda * Constants::c() / vth;
+
+        // select the isotropic or the dipole phase function:
+        // all wing events and 1/3 of core events are dipole, and the remaining 2/3 core events are isotropic,
+        // where x=0.2 (in the atom frame) defines the transition between core and wings
+        scatinfo->dipole = abs(x) > 0.2 || random()->uniform() < 1. / 3.;
+
+        scatinfo->velocity =
+            LyUtils::sampleAtomVelocity(lambda, lya, vth, a, T, nH, pp->direction(), config(), random());
     }
+}
+
+////////////////////////////////////////////////////////////////////
+
+bool LyaNeutralHydrogenGasMix::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda,
+                                                 Direction bfkobs, Direction bfky, const MaterialState* state,
+                                                 const PhotonPacket* pp) const
+{
+    setScatteringInfoIfNeeded(const_cast<PhotonPacket*>(pp), state, lambda);
+    auto scatinfo = const_cast<PhotonPacket*>(pp)->getScatteringInfo();
 
     // add the contribution to the Stokes vector components depending on scattering type
     if (scatinfo->dipole)
@@ -158,7 +199,7 @@ bool LyaNeutralHydrogenGasMix::peeloffScattering(double& I, double& Q, double& U
     }
 
     // Doppler-shift the photon packet wavelength into and out of the atom frame
-    lambda = LyaUtils::shiftWavelength(lambda, scatinfo->velocity, pp->direction(), bfkobs);
+    lambda = LyUtils::shiftWavelength(lambda, scatinfo->velocity, pp->direction(), bfkobs);
 
     return false;
 }
@@ -167,14 +208,8 @@ bool LyaNeutralHydrogenGasMix::peeloffScattering(double& I, double& Q, double& U
 
 void LyaNeutralHydrogenGasMix::performScattering(double lambda, const MaterialState* state, PhotonPacket* pp) const
 {
-    // draw a random atom velocity & phase function, unless a peel-off stored this already
-    auto scatinfo = pp->getScatteringInfo();
-    if (!scatinfo->valid)
-    {
-        scatinfo->valid = true;
-        std::tie(scatinfo->velocity, scatinfo->dipole) = LyaUtils::sampleAtomVelocity(
-            lambda, state->temperature(), state->numberDensity(), pp->direction(), config(), random());
-    }
+    setScatteringInfoIfNeeded(const_cast<PhotonPacket*>(pp), state, lambda);
+    auto scatinfo = const_cast<PhotonPacket*>(pp)->getScatteringInfo();
 
     // draw the outgoing direction from the dipole or the isotropic phase function
     // and, if required, update the polarization state of the photon packet
@@ -190,7 +225,7 @@ void LyaNeutralHydrogenGasMix::performScattering(double lambda, const MaterialSt
     }
 
     // Doppler-shift the photon packet wavelength into and out of the atom frame
-    lambda = LyaUtils::shiftWavelength(lambda, scatinfo->velocity, pp->direction(), bfknew);
+    lambda = LyUtils::shiftWavelength(lambda, scatinfo->velocity, pp->direction(), bfknew);
 
     // execute the scattering event in the photon packet
     pp->scatter(bfknew, state->bulkVelocity(), lambda);

SKIRT/core/LyaNeutralHydrogenGasMix.hpp

@@ -11,9 +11,9 @@
 
 ////////////////////////////////////////////////////////////////////
 
-/** The LyaNeutralHydrogenGasMix class describes the material properties related to
-    Lyman-alpha line transfer for a population of neutral hydrogen atoms, including support for
-    polarization by scattering.
+/** The LyaNeutralHydrogenGasMix class describes the material properties related to Lyman-alpha
+    line transfer for a population of neutral hydrogen atoms, including support for polarization by
+    scattering.
 
     The spatial distributions for both the mass density and the temperature of the neutral hydrogen
     gas must be defined by the input model and are considered to be constant during the simulation.
@@ -89,8 +89,8 @@ class LyaNeutralHydrogenGasMix : public MaterialMix
         density, so this function returns a list containing these two items. */
     vector<StateVariable> specificStateVariableInfo() const override;
 
-    /** This function initializes the specific state variables requested by this fragmented dust
-        mix through the specificStateVariableInfo() function except for the number density. For the
+    /** This function initializes the specific state variables requested by this material mix
+        through the specificStateVariableInfo() function except for the number density. For the
         Lyman-alpha material mix, the function initializes the temperature to the specified
         imported temperature, or if this is not available, to the user-configured default
         temperature for this material mix. The metallicity and custom parameter arguments are
@@ -100,6 +100,11 @@ class LyaNeutralHydrogenGasMix : public MaterialMix
 
     //======== Low-level material properties =======
 
+private:
+    /** This private function calculates the cross section per hydrogen atom for the given
+        wavelength and temperature. */
+    double section(double lambda, double T) const;
+
 public:
     /** This function returns the mass of neutral hydrogen atom. */
     double mass() const override;
@@ -136,6 +141,15 @@ class LyaNeutralHydrogenGasMix : public MaterialMix
         Lyman-alpha material mix, the extinction opacity equals the scattering opacity. */
     double opacityExt(double lambda, const MaterialState* state, const PhotonPacket* pp) const override;
 
+private:
+    /** This private function draws an atom velocity and phase function and stores this information
+        in the photon packet's scattering information record, unless a previous peel-off stored
+        this already. The atom velocity sampling is delegated to the LyUtils::sampleAtomVelocity.
+        The phase function is chosen to be a dipole for all wing events (x > 0.2) and 1/3 of the
+        core events (x < 0.2). The remaining core events are isotropic. */
+    void setScatteringInfoIfNeeded(PhotonPacket* pp, const MaterialState* state, const double lambda) const;
+
+public:
     /** 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

SKIRT/core/MediumSystem.cpp

@@ -10,7 +10,6 @@
 #include "FatalError.hpp"
 #include "LockFree.hpp"
 #include "Log.hpp"
-#include "LyaUtils.hpp"
 #include "MaterialMix.hpp"
 #include "MaterialState.hpp"
 #include "NR.hpp"

SKIRT/core/SelectDustMixFamily.cpp

@@ -4,7 +4,6 @@
 ///////////////////////////////////////////////////////////////// */
 
 #include "SelectDustMixFamily.hpp"
-#include "Constants.hpp"
 
 ////////////////////////////////////////////////////////////////////
 

SKIRT/core/SimulationItemRegistry.cpp

@@ -293,6 +293,7 @@
 #include "VoronoiMeshSpatialGrid.hpp"
 #include "WeingartnerDraineDustMix.hpp"
 #include "XRayAtomicGasMix.hpp"
+#include "XRayIonicGasMix.hpp"
 #include "ZubkoDustMix.hpp"
 #include "ZubkoGraphiteGrainSizeDistribution.hpp"
 #include "ZubkoPAHGrainSizeDistribution.hpp"
@@ -597,6 +598,7 @@ SimulationItemRegistry::SimulationItemRegistry(string version, string format)
     ItemRegistry::add<SpinFlipAbsorptionMix>();
     ItemRegistry::add<SpinFlipHydrogenGasMix>();
     ItemRegistry::add<XRayAtomicGasMix>();
+    ItemRegistry::add<XRayIonicGasMix>();
     ItemRegistry::add<EmittingGasMix>();
     ItemRegistry::add<NonLTELineGasMix>();
     ItemRegistry::add<LyaNeutralHydrogenGasMix>();