Make instruments record fluorescence as secondary emission (SKIRT#256)

* detect emulation in configuration
* add no-op scattering to base class
* honor emulation in medium system
* honor emulation in instruments
* store rf according to phase

Author: petercamps

SKIRT/core/Configuration.cpp

@@ -11,7 +11,6 @@
 #include "MaterialMix.hpp"
 #include "MaterialWavelengthRangeInterface.hpp"
 #include "MonteCarloSimulation.hpp"
-#include "NR.hpp"
 #include "OligoWavelengthDistribution.hpp"
 #include "OligoWavelengthGrid.hpp"
 #include "ProbeSystem.hpp"
@@ -325,9 +324,13 @@ void Configuration::setupSelfBefore()
     _hasSingleConstantSectionMedium = numMedia == 1 && _hasConstantSectionMedium;
     _hasMultipleConstantSectionMedia = numMedia > 1 && _hasConstantSectionMedium;
 
-    // check for scattering dispersion
+    // check for scattering dispersion and for emulating secondary emission
     for (auto medium : ms->media())
+    {
         if (medium->mix()->hasScatteringDispersion()) _hasScatteringDispersion = true;
+        if (medium->mix()->scatteringEmulatesSecondaryEmission()) _scatteringEmulatesSecondaryEmission = true;
+    }
+    _needIndividualPeelOff = _hasScatteringDispersion | _scatteringEmulatesSecondaryEmission;
 
     // check for magnetic fields
     for (int h = 0; h != numMedia; ++h)

SKIRT/core/Configuration.hpp

@@ -209,9 +209,13 @@ class Configuration : public SimulationItem
         weight factors. */
     bool hasMultipleConstantSectionMedia() const { return _hasMultipleConstantSectionMedia; }
 
+    /** Returns true if a scattering interaction for one or more media may emulate secondary
+        emission, and false otherwise. */
+    bool scatteringEmulatesSecondaryEmission() const { return _scatteringEmulatesSecondaryEmission; }
+
     /** Returns true if a scattering interaction for one or more media may adjust the wavelength of
-        the interacting photon packet, and false otherwise. */
-    bool hasScatteringDispersion() const { return _hasScatteringDispersion; }
+        the interacting photon packet or may emulate secondary emission, and false otherwise. */
+    bool needIndividualPeelOff() const { return _needIndividualPeelOff; }
 
     /** Returns true if all media in the simulation support polarization, and false if none of the
         media do. A mixture of support and no support for polarization is not allowed and will
@@ -467,6 +471,8 @@ class Configuration : public SimulationItem
     bool _hasSingleConstantSectionMedium{false};
     bool _hasMultipleConstantSectionMedia{false};
     bool _hasScatteringDispersion{false};
+    bool _scatteringEmulatesSecondaryEmission{false};
+    bool _needIndividualPeelOff{false};
     bool _hasPolarization{false};
     bool _hasSpheroidalPolarization{false};
 

SKIRT/core/DustMix.cpp

@@ -427,7 +427,7 @@ namespace
 
 ////////////////////////////////////////////////////////////////////
 
-void DustMix::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs,
+bool DustMix::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs,
                                 Direction bfky, const MaterialState* /*state*/, const PhotonPacket* pp) const
 {
     switch (scatteringMode())
@@ -482,6 +482,7 @@ void DustMix::peeloffScattering(double& I, double& Q, double& U, double& V, doub
             break;
         }
     }
+    return false;
 }
 
 ////////////////////////////////////////////////////////////////////

SKIRT/core/DustMix.hpp

@@ -277,7 +277,7 @@ class DustMix : public MaterialMix
         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. */
-    void peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs, Direction bfky,
+    bool peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs, Direction bfky,
                            const MaterialState* state, const PhotonPacket* pp) const override;
 
     /** This function performs a scattering event on the specified photon packet in the spatial

SKIRT/core/ElectronMix.cpp

@@ -171,7 +171,7 @@ namespace
 
 ////////////////////////////////////////////////////////////////////
 
-void ElectronMix::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs,
+bool ElectronMix::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs,
                                     Direction bfky, const MaterialState* state, const PhotonPacket* pp) const
 {
     // if we have dispersion, adjust the incoming wavelength to the electron rest frame
@@ -197,6 +197,7 @@ void ElectronMix::peeloffScattering(double& I, double& Q, double& U, double& V,
     {
         lambda = PhotonPacket::shiftedEmissionWavelength(lambda, bfkobs, scatinfo->velocity);
     }
+    return false;
 }
 
 ////////////////////////////////////////////////////////////////////

SKIRT/core/ElectronMix.hpp

@@ -170,7 +170,7 @@ class ElectronMix : public MaterialMix
 
         If \em includePolarization has been set to true, the function supports polarization;
         otherwise it does not. */
-    void peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs, Direction bfky,
+    bool peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs, Direction bfky,
                            const MaterialState* state, const PhotonPacket* pp) const override;
 
     /** This function performs a scattering event on the specified photon packet in the spatial

SKIRT/core/FragmentDustMixDecorator.cpp

@@ -213,7 +213,7 @@ double FragmentDustMixDecorator::opacityExt(double lambda, const MaterialState*
 
 ////////////////////////////////////////////////////////////////////
 
-void FragmentDustMixDecorator::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda,
+bool FragmentDustMixDecorator::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda,
                                                  Direction bfkobs, Direction bfky, const MaterialState* state,
                                                  const PhotonPacket* pp) const
 {
@@ -240,6 +240,7 @@ void FragmentDustMixDecorator::peeloffScattering(double& I, double& Q, double& U
             V += Vf * wf;
         }
     }
+    return false;
 }
 
 ////////////////////////////////////////////////////////////////////

SKIRT/core/FragmentDustMixDecorator.hpp

@@ -184,7 +184,7 @@ class FragmentDustMixDecorator : public MaterialMix, public MultiGrainPopulation
         the fragment. The contributions to the Stokes vector components are stored in the \em I,
         \em Q, \em U, \em V arguments, which are guaranteed to be initialized to zero by the
         caller. For dust mixes, the wavelength remains unchanged. */
-    void peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs, Direction bfky,
+    bool peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs, Direction bfky,
                            const MaterialState* state, const PhotonPacket* pp) const override;
 
     /** This function performs a scattering event on the specified photon packet in the spatial

SKIRT/core/Instrument.cpp

@@ -5,7 +5,6 @@
 
 #include "Instrument.hpp"
 #include "Configuration.hpp"
-#include "FatalError.hpp"
 #include "FluxRecorder.hpp"
 
 ////////////////////////////////////////////////////////////////////
@@ -20,7 +19,7 @@ void Instrument::setupSelfBefore()
 
     // discover details about the simulation
     bool hasMedium = config->hasMedium();
-    bool hasMediumEmission = config->hasSecondaryEmission();
+    bool hasMediumEmission = config->hasSecondaryEmission() || config->scatteringEmulatesSecondaryEmission();
 
     // partially configure the flux recorder
     _recorder = new FluxRecorder(this);

SKIRT/core/LyaNeutralHydrogenGasMix.cpp

@@ -132,7 +132,7 @@ double LyaNeutralHydrogenGasMix::opacityExt(double lambda, const MaterialState*
 
 ////////////////////////////////////////////////////////////////////
 
-void LyaNeutralHydrogenGasMix::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda,
+bool LyaNeutralHydrogenGasMix::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda,
                                                  Direction bfkobs, Direction bfky, const MaterialState* state,
                                                  const PhotonPacket* pp) const
 {
@@ -159,6 +159,8 @@ void 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);
+
+    return false;
 }
 
 ////////////////////////////////////////////////////////////////////