diff --git a/SKIRT/CMakeLists.txt b/SKIRT/CMakeLists.txt index 59d52eca..7be8507b 100644 --- a/SKIRT/CMakeLists.txt +++ b/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) diff --git a/SKIRT/core/CMakeLists.txt b/SKIRT/core/CMakeLists.txt index dd630882..d110f8dd 100644 --- a/SKIRT/core/CMakeLists.txt +++ b/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 diff --git a/SKIRT/core/LyaUtils.cpp b/SKIRT/core/LyUtils.cpp similarity index 60% rename from SKIRT/core/LyaUtils.cpp rename to SKIRT/core/LyUtils.cpp index 4452e3e6..98592e58 100644 --- a/SKIRT/core/LyaUtils.cpp +++ b/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 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 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); } diff --git a/SKIRT/core/LyaUtils.hpp b/SKIRT/core/LyUtils.hpp similarity index 96% rename from SKIRT/core/LyaUtils.hpp rename to SKIRT/core/LyUtils.hpp index fb90c2bf..1b0eaf53 100644 --- a/SKIRT/core/LyaUtils.hpp +++ b/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 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 diff --git a/SKIRT/core/LyaNeutralHydrogenGasMix.cpp b/SKIRT/core/LyaNeutralHydrogenGasMix.cpp index 7963f4f2..c2000f20 100644 --- a/SKIRT/core/LyaNeutralHydrogenGasMix.cpp +++ b/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,7 +139,7 @@ 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.; } //////////////////////////////////////////////////////////////////// @@ -127,23 +147,44 @@ double LyaNeutralHydrogenGasMix::opacitySca(double lambda, const MaterialState* 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(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(pp), state, lambda); + auto scatinfo = const_cast(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(pp), state, lambda); + auto scatinfo = const_cast(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); diff --git a/SKIRT/core/LyaNeutralHydrogenGasMix.hpp b/SKIRT/core/LyaNeutralHydrogenGasMix.hpp index dec954db..58626efd 100644 --- a/SKIRT/core/LyaNeutralHydrogenGasMix.hpp +++ b/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 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 diff --git a/SKIRT/core/MediumSystem.cpp b/SKIRT/core/MediumSystem.cpp index 79e19f1d..5815a67c 100644 --- a/SKIRT/core/MediumSystem.cpp +++ b/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" diff --git a/SKIRT/core/SelectDustMixFamily.cpp b/SKIRT/core/SelectDustMixFamily.cpp index b7a3c9cb..2eba82f3 100644 --- a/SKIRT/core/SelectDustMixFamily.cpp +++ b/SKIRT/core/SelectDustMixFamily.cpp @@ -4,7 +4,6 @@ ///////////////////////////////////////////////////////////////// */ #include "SelectDustMixFamily.hpp" -#include "Constants.hpp" //////////////////////////////////////////////////////////////////// diff --git a/SKIRT/core/SimulationItemRegistry.cpp b/SKIRT/core/SimulationItemRegistry.cpp index cea5df46..5f05b8a5 100644 --- a/SKIRT/core/SimulationItemRegistry.cpp +++ b/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(); ItemRegistry::add(); ItemRegistry::add(); + ItemRegistry::add(); ItemRegistry::add(); ItemRegistry::add(); ItemRegistry::add(); diff --git a/SKIRT/core/VoigtProfile.cpp b/SKIRT/core/VoigtProfile.cpp index 4b38a9a8..38d47e1a 100644 --- a/SKIRT/core/VoigtProfile.cpp +++ b/SKIRT/core/VoigtProfile.cpp @@ -6,34 +6,14 @@ #include "VoigtProfile.hpp" #include "FatalError.hpp" #include "Random.hpp" +#include "Faddeeva.hh" //////////////////////////////////////////////////////////////////// double VoigtProfile::value(double a, double x) { - // coefficients for the approximation function (Table A1, Smith+15) - constexpr double A0 = 15.75328153963877; - constexpr double A1 = 286.9341762324778; - constexpr double A2 = 19.05706700907019; - constexpr double A3 = 28.22644017233441; - constexpr double A4 = 9.526399802414186; - constexpr double A5 = 35.29217026286130; - constexpr double A6 = 0.8681020834678775; - constexpr double B0 = 0.0003300469163682737; - constexpr double B1 = 0.5403095364583999; - constexpr double B2 = 2.676724102580895; - constexpr double B3 = 12.82026082606220; - constexpr double B4 = 3.21166435627278; - constexpr double B5 = 32.032981933420; - constexpr double B6 = 9.0328158696; - constexpr double B7 = 23.7489999060; - constexpr double B8 = 1.82106170570; - - // calculation of the approximation (Appendix A1, Smith+15) - double z = x * x; - if (z <= 3.0) return exp(-z) * (1.0 - a * (A0 + A1 / (z - A2 + A3 / (z - A4 + A5 / (z - A6))))); - if (z < 25.0) return exp(-z) + a * (B0 + B1 / (z - B2 + B3 / (z + B4 + B5 / (z - B6 + B7 / (z - B8))))); - return 0.5 * M_2_SQRTPI * a / (z - 1.5 - 1.5 / (z - 3.5 - 5.0 / (z - 5.5))); + std::complex z(x, a); + return std::real(Faddeeva::w(z, 0.001)); } //////////////////////////////////////////////////////////////////// @@ -70,7 +50,7 @@ double VoigtProfile::sample(double a, double x, Random* random) double p = (theta0 + M_PI_2) / ((1. - exp(-u0 * u0)) * theta0 + (1. + exp(-u0 * u0)) * M_PI_2); // perform the rejection method loop for a maximum number of attempts - int n = 10000; + int n = 100000; // not the most efficient for certain a,x combinations (eg. Zn+29 at T=3K) while (n--) { // determine which one of the two comparison functions to use diff --git a/SKIRT/core/VoigtProfile.hpp b/SKIRT/core/VoigtProfile.hpp index 73cb97bf..9e82a032 100644 --- a/SKIRT/core/VoigtProfile.hpp +++ b/SKIRT/core/VoigtProfile.hpp @@ -6,7 +6,6 @@ #ifndef VOIGTPROFILE_HPP #define VOIGTPROFILE_HPP -#include "Basics.hpp" class Random; //////////////////////////////////////////////////////////////////// diff --git a/SKIRT/core/XRayAtomicGasMix.hpp b/SKIRT/core/XRayAtomicGasMix.hpp index 129c99d7..80393838 100644 --- a/SKIRT/core/XRayAtomicGasMix.hpp +++ b/SKIRT/core/XRayAtomicGasMix.hpp @@ -88,9 +88,9 @@ implementation accurately reflects the physical processes, apart from approximations caused by discretization, tabulation, and interpolation. - Finally, \em H2fraction sets the fraction of hydrogen atoms locked up in molecular H2, for which - the Rayleigh scattering cross section is enhanced due to coherent scattering on the H2 molecule, - following the exact same implementation as the MONACO code (Odaka+11). + Finally, \em H2fraction sets the fraction of hydrogen atoms locked up in molecular H2, for + which the Rayleigh scattering cross section is enhanced due to coherent scattering on the H2 + molecule, following the exact same implementation as the MONACO code (Odaka+11). Configuring the simulation @@ -102,13 +102,15 @@ with a high-resolution wavelength grid, or separate instruments can be configured with wavelength grids to resolve specific features of interest. - To study fluorescence lines seperated from the background continuum, configure the instrument - to record flux components seperately and consider the scattered flux component, which includes - the fluorescence lines in addition to any flux scattered from other media components. To ensure - that low-intensity lines are properly included in this flux, set the advanced property \em - minWeightReduction in the \c PhotonPacketOptions section to a value of \c 1e10 or so. With the - default value of \c 1e4, low-intensity fluorescence photon packets are killed before having a - chance to register in the instruments. + Although there is no secondary emission phase, fluorescent emission, which is modelled as + scattering, does contribute to the secondary emission spectrum. To study this process + seperately from the background continuum, configure the instrument to record flux components + seperately and consider the secondary emission component, which includes the fluorescence lines + in addition to any flux scattered from these lines. To ensure that low-intensity lines are + properly included in this flux, set the advanced property \em minWeightReduction in the \c + PhotonPacketOptions section to a value of \c 1e10 or so. With the default value of \c 1e4, + low-intensity fluorescence photon packets are killed before having a chance to register in the + instruments. The input model must define the spatial distribution of the hydrogen number density \f$n_\mathrm{H} = n_\mathrm{HI} + n_\mathrm{HII} + 2\,n_\mathrm{H2}\f$, i.e. including atomic, @@ -131,9 +133,10 @@ of Anders & Grevesse (1989), the default abundance table in Xspec. In the default list, the abundance of hydrogen is set to unity. However, it is acceptable to specify a hydrogen abundance lower than one, for example to model an ionized hydrogen fraction. The first element - in the abundance list represents all hydrogen atoms, both atomic and those in H2 (i.e. 2 per H2). - The \em H2fraction (configured in the ski file) then sets the fraction of these hydrogen atoms - that are locked up in molecular H2 (for which the Rayleigh scattering cross section is enhanced). + in the abundance list represents all hydrogen atoms, both atomic and those in H2 (i.e. 2 per + H2). The \em H2fraction (configured in the ski file) then sets the fraction of these hydrogen + atoms that are locked up in molecular H2 (for which the Rayleigh scattering cross section is + enhanced). Photo-absorption cross section @@ -203,9 +206,9 @@ available as a table (for each H atom in H2, we have \f$2\times \sigma_{RSA, H}(E)\f$). The normalised scattering phase function for element Z is now given by \f[ \Phi_{RSA, Z}(\theta, E)= \frac{3}{4}\, \frac{\sigma_T}{\sigma_{RSA, Z}(E)}\Big[ 1 + \cos^2\theta \Big] \cdot \Big[ - \big(F_Z(q) + F'_Z(E)\big)^2 + {F''_Z}^2(E)\Big], \f] with the same atomic form factors \f$F_Z(q)\f$ - as before, and tabulated real and imaginary anomalous scattering functions \f$F'_Z(E)\f$ and - \f$F''_Z(E)\f$. + \big(F_Z(q) + F'_Z(E)\big)^2 + {F''_Z}^2(E)\Big], \f] with the same atomic form factors + \f$F_Z(q)\f$ as before, and tabulated real and imaginary anomalous scattering functions + \f$F'_Z(E)\f$ and \f$F''_Z(E)\f$. Electron scattering - photon energy shift @@ -466,7 +469,7 @@ class XRayAtomicGasMix : public MaterialMix class ScatteringHelper; private: - // all data members are precalculated in setupSelfAfter() + // all data members are precalculated in setupSelfBefore() // wavelength grid (shifted to the left of the actually sampled points to approximate rounding) Array _lambdav; // indexed on ell diff --git a/SKIRT/core/XRayIonicGasMix.cpp b/SKIRT/core/XRayIonicGasMix.cpp new file mode 100644 index 00000000..b6ff5c25 --- /dev/null +++ b/SKIRT/core/XRayIonicGasMix.cpp @@ -0,0 +1,1021 @@ +/*////////////////////////////////////////////////////////////////// +//// The SKIRT project -- advanced radiative transfer //// +//// © Astronomical Observatory, Ghent University //// +///////////////////////////////////////////////////////////////// */ +#include "XRayIonicGasMix.hpp" +#include "AtomUtils.hpp" +#include "ComptonPhaseFunction.hpp" +#include "Configuration.hpp" +#include "Constants.hpp" +#include "DipolePhaseFunction.hpp" +#include "FatalError.hpp" +#include "LyUtils.hpp" +#include "MaterialState.hpp" +#include "NR.hpp" +#include "PhotonPacket.hpp" +#include "Random.hpp" +#include "Range.hpp" +#include "StoredTable.hpp" +#include "StringUtils.hpp" +#include "TextInFile.hpp" +#include +#include + +//////////////////////////////////////////////////////////////////// + +namespace +{ + // ---- common helper functions ---- + + static constexpr int numAtoms = 30; // maximum atomic number used in this class + + // convert photon energy in eV to and from wavelength in m (same conversion in both directions) + constexpr double wavelengthToFromEnergy(double x) + { + constexpr double front = Constants::h() * Constants::c() / Constants::Qelectron(); + return front / x; + } + + // ---- hardcoded configuration constants ---- + + // wavelength range over which our cross sections may be nonzero + constexpr Range nonZeroRange(wavelengthToFromEnergy(500e3), wavelengthToFromEnergy(4.3)); + + // number of wavelengths per dex in high-resolution grid + constexpr size_t numWavelengthsPerDex = 2500; + + // load data from resource file with N columns into a vector of structs of type S that can be constructed + // from an array with N elements, and return that vector + template vector loadStruct(const SimulationItem* item, string filename, string description) + { + vector result; + TextInFile infile(item, filename, description, true); + for (int i = 0; i != N; ++i) infile.addColumn(string()); + Array row; + while (infile.readRow(row)) result.emplace_back(row); + return result; + } + + // resource data for photo-absorption + struct PhotoAbsorbResource + { + PhotoAbsorbResource(const Array& a) + : Z(a[0]), N(a[1]), n(a[2]), l(a[3]), Eth(a[4]), E0(a[5]), sigma0(a[6]), ya(a[7]), P(a[8]), yw(a[9]) + {} + + int ionIndex{-1}; // index of the ion + int Z; // atomic number + int N; // number of electrons + int n; // principal quantum number of the shell + int l; // orbital quantum number of the subshell + double Eth; // subshell ionization threshold energy (eV) + double E0; // fit parameter (eV) + double sigma0; // fit parameter (Mb = 10^-22 m^2) + double ya; // fit parameter (1) + double P; // fit parameter (1) + double yw; // fit parameter (1) + static constexpr double Emax = 5e5; // maximum energy for validity of the formula (eV) + + double Es; + double sigmamax; + + // return photo-absorption cross section in m2 for given energy in eV and cross section parameters, + // without taking into account thermal dispersion + double photoAbsorbSection(double E) const + { + if (E < Eth || E >= Emax) return 0.; + + double x = E / E0; + double y = x; + double xm1 = x - 1.; + double Q = 5.5 + l - 0.5 * P; + double F = (xm1 * xm1 + yw * yw) * std::pow(y, -Q) * std::pow(1. + std::sqrt(y / ya), -P); + return 1e-22 * sigma0 * F; // from Mb to m2 + } + + // return photo-absorption cross section in m2 for given energy in eV and cross section parameters, + // approximating thermal dispersion by replacing the steep threshold transition by a sigmoid + // error function with given parameters (dispersion and maximum value) + double photoAbsorbThermalSection(double E) const + { + if (E <= Eth - 2. * Es) return 0.; + if (E >= Eth + 2. * Es) return photoAbsorbSection(E); + return sigmamax * (0.5 + 0.5 * std::erf((E - Eth) / Es)); + } + }; + + // resource data for fluorescence + struct FluorescenceResource + { + FluorescenceResource(const Array& a) : Z(a[0]), N(a[1]), n(a[2]), l(a[3]), omega(a[4]), E(a[5]), W(a[6]) {} + + int paIndex{-1}; // index of the photo-absorption transition + int Z; // atomic number + int N; // number of electrons + int n; // principal quantum number of the shell with the hole + int l; // orbital quantum number of the subshell with the hole + double omega; // fluorescence yield (1) + double E; // (central) energy of the emitted photon (eV) + double W; // FWHM of the Lorentz shape for the emitted photon (eV), or zero + }; + + // resource data for Lyman-series + struct LymanResource + { + LymanResource(const Array& a) : Z(a[0]), index(a[1]), lamA(a[2]), lam(a[3]) {} + int ionIndex{-1}; // index of the ion + double sprob{-1}; // scatter probability after resonant scattering (<=1) + double vth{-1}; // sqrt(2) * thermal velocity (m/s) + int Z; // atomic number + int index; // Lyman index (alpha1/2, alpha3/2, beta1/2, ...) + double lamA; // wavelength * Einstein A (m/s) + double lam; // wavelength (m) + + double section(double lambda) const + { + double a = lamA / (4. * M_PI * vth); + + double g = (index % 2) + 1.; + return LyUtils::section(lambda, lam, vth, a, g); + } + }; + + // resource data for Lyman branching (incoherent scattering) + struct LymanBranchResource + { + LymanBranchResource(const Array& a) : Z(a[0]), upper(a[1]), lower(a[2]), prob(a[3]) {} + int Z; // atomic number + int upper; // upper Lyman index + int lower; // lower Lyman index + double prob; // branching probability + }; +} + +//////////////////////////////////////////////////////////////////// + +// ---- base class for scattering helpers ---- + +class XRayIonicGasMix::ScatteringHelper +{ +public: + virtual ~ScatteringHelper() {} + + // return scattering cross section for atom in m2 + virtual double sectionSca(double lambda, int Z) const = 0; + + // peel-off unpolarized scattering event: override this in helpers that don't support polarization + virtual void peeloffScattering(double& /*I*/, double& /*lambda*/, int /*Z*/, Direction /*bfk*/, + Direction /*bfkobs*/) const + { + // default implementation does nothing + } + + // perform unpolarized scattering event: override this in helpers that don't support polarization + virtual Direction performScattering(double& /*lambda*/, int /*Z*/, Direction /*bfk*/) const + { + // default implementation returns null vector + return Direction(); + } + + // peel-off polarized scattering event: override this in helpers that do support polarization + virtual void peeloffScattering(double& I, double& /*Q*/, double& /*U*/, double& /*V*/, double& lambda, int Z, + Direction bfk, Direction bfkobs, Direction /*bfky*/, + const StokesVector* /*sv*/) const + { + // default implementation calls unpolarized version + peeloffScattering(I, lambda, Z, bfk, bfkobs); + } + + // perform polarized scattering event: override this in helpers that do support polarization + virtual Direction performScattering(double& lambda, int Z, Direction bfk, StokesVector* /*sv*/) const + { + // default implementation calls unpolarized version + return performScattering(lambda, Z, bfk); + } +}; + +//////////////////////////////////////////////////////////////////// + +// ---- no scattering helper ---- + +namespace +{ + // this helper does nothing; it is used as a stub in case there is no scattering of a given type + class NoScatteringHelper : public XRayIonicGasMix::ScatteringHelper + { + public: + NoScatteringHelper(SimulationItem* /*item*/) {} + + double sectionSca(double /*lambda*/, int /*Z*/) const override { return 0.; } + }; +} + +//////////////////////////////////////////////////////////////////// + +// ---- free-electron Compton scattering helper ---- + +namespace +{ + // transition wavelength from Compton to Thomson scattering + constexpr double comptonWL = wavelengthToFromEnergy(100.); // 0.1 keV or 12.4 nm + + // this helper forwards all calls to an external helper class for regular Compton scattering + // (or Thomson scattering for lower energies, because Compton becomes numerically unstable) + class FreeComptonHelper : public XRayIonicGasMix::ScatteringHelper + { + private: + ComptonPhaseFunction _cpf; + DipolePhaseFunction _dpf; + + public: + FreeComptonHelper(SimulationItem* item) + { + auto random = item->find(); + _cpf.initialize(random); + _dpf.initialize(random); + } + + double sectionSca(double lambda, int Z) const override + { + double sigma = Z * Constants::sigmaThomson(); + if (lambda < comptonWL) sigma *= _cpf.sectionSca(lambda); + return sigma; + } + + void peeloffScattering(double& I, double& lambda, int /*Z*/, Direction bfk, Direction bfkobs) const override + { + if (lambda < comptonWL) + { + double Q, U, V; + _cpf.peeloffScattering(I, Q, U, V, lambda, bfk, bfkobs, Direction(), nullptr); + } + else + { + double Q, U, V; + _dpf.peeloffScattering(I, Q, U, V, bfk, bfkobs, Direction(), nullptr); + } + } + + Direction performScattering(double& lambda, int /*Z*/, Direction bfk) const override + { + return lambda < comptonWL ? _cpf.performScattering(lambda, bfk, nullptr) + : _dpf.performScattering(bfk, nullptr); + } + }; +} + +//////////////////////////////////////////////////////////////////// + +// ---- free-electron Compton with polarization scattering helper ---- + +namespace +{ + // this helper forwards all calls to an external helper class for Compton scattering + // (or Thomson scattering for lower energies) with support for polarization + class FreeComptonWithPolarizationHelper : public XRayIonicGasMix::ScatteringHelper + { + private: + ComptonPhaseFunction _cpf; + DipolePhaseFunction _dpf; + + public: + FreeComptonWithPolarizationHelper(SimulationItem* item) + { + auto random = item->find(); + _cpf.initialize(random, true); + _dpf.initialize(random, true); + } + + double sectionSca(double lambda, int Z) const override + { + double sigma = Z * Constants::sigmaThomson(); + if (lambda < comptonWL) sigma *= _cpf.sectionSca(lambda); + return sigma; + } + + void peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, int /*Z*/, Direction bfk, + Direction bfkobs, Direction bfky, const StokesVector* sv) const override + { + lambda < comptonWL ? _cpf.peeloffScattering(I, Q, U, V, lambda, bfk, bfkobs, bfky, sv) + : _dpf.peeloffScattering(I, Q, U, V, bfk, bfkobs, bfky, sv); + } + + Direction performScattering(double& lambda, int /*Z*/, Direction bfk, StokesVector* sv) const override + { + return lambda < comptonWL ? _cpf.performScattering(lambda, bfk, sv) : _dpf.performScattering(bfk, sv); + } + }; +} + +//////////////////////////////////////////////////////////////////// + +void XRayIonicGasMix::setupSelfBefore() +{ + MaterialMix::setupSelfBefore(); + + auto config = find(); + + // ------------ parse user properties ------------ + + // parse all required ions from the ions property + string ionString = StringUtils::squeeze(ions()); + if (ionString.empty()) throw FATALERROR("No ions specified"); + for (string ion : StringUtils::split(ionString, ",")) + { + int Z, N; + std::tie(Z, N) = AtomUtils::parseIon(ion); + _ionParamv.emplace_back(Z, N); + } + _numIons = _ionParamv.size(); + + if (_numIons != static_cast(abundances().size())) + throw FATALERROR("Number of ions and abundances do not match"); + + // create scattering helpers depending on the user-configured implementation type + switch (electronScattering()) + { + case ElectronScattering::None: _com = new NoScatteringHelper(this); break; + case ElectronScattering::Free: _com = new FreeComptonHelper(this); break; + case ElectronScattering::FreeWithPolarization: _com = new FreeComptonWithPolarizationHelper(this); break; + } + if (resonantScattering()) + { + _dpf = new DipolePhaseFunction(); + _dpf->initialize(random(), true); + } + + // resources that are maintained during the setup + vector usedPar; + vector usedFlr; + vector usedLyr; + + // Use nested scope to load and preprocess resources and discard unused resources + { + // ------------ load full resources ------------ + + // photo-absorption data + auto paResource = loadStruct(this, "Ionic_PA.txt", "photoabsorption data"); + // fluorescence data + auto flResource = loadStruct(this, "Ionic_FL.txt", "fluorescence data"); + // generic Lyman series data + auto lyResource = loadStruct(this, "Ionic_LY.txt", "lyman series data"); + // Lyman recombination temperature-dependent yields + StoredTable<3> lyyResource(this, "Ionic_LY_Y.stab", "Z(1),Ly(1),T(K)", "Y(1)"); + // Lyman branching probabilities + vector lybResource; + if (resonantScattering()) + lybResource = loadStruct(this, "Ionic_LY_B.txt", "branching probabilities"); + + // ------------ preprocess resources ------------ + + // Lyman recombination can be modelled as fluorescence following an inner shell PA (n,l)=(1,0) + // This ignores the cascade and only models the transition back to the inner shell. + // There is no PA data for (n,l)=(1,0) so we can simply add them without worrying about duplicates. + flResource.reserve(flResource.size() + lyResource.size()); + for (const auto& lyr : lyResource) + { + double Z = lyr.Z; + double Ly = lyr.index; + double E = wavelengthToFromEnergy(lyr.lam); + double omega = lyyResource(Z, Ly, temperature()); + Array params = {Z, 1, 1, 0, omega, E, 0.}; // ensure order is correct here! + flResource.emplace_back(params); + } + + // ------------ discard unused resources ------------ + + // for each (used) ion save all the photo-absorption, fluorescence and resonant Lyman transitions + + // for each ION + for (int i = 0; i < _numIons; i++) + { + auto& ion = _ionParamv[i]; + + // add all PA for this ION + for (auto& pa : paResource) + { + if (pa.Z == ion.Z && pa.N == ion.N) + { + pa.ionIndex = i; + usedPar.push_back(pa); + + // add all FL for this PA + for (auto& fl : flResource) + { + if (fl.Z == pa.Z && fl.N == pa.N && fl.n == pa.n && fl.l == pa.l) + { + int p = usedPar.size() - 1; + fl.paIndex = p; + usedFlr.push_back(fl); + } + } + } + } + + // add RS for this (H-like) ION + if (resonantScattering() && ion.N == 1) + { + for (auto& ly : lyResource) + { + if (ly.Z == ion.Z) + { + ly.ionIndex = i; + usedLyr.push_back(ly); + } + } + } + } + _numFluo = usedFlr.size(); + _numLym = usedLyr.size(); + + // ------------ postprocess used resources ------------ + + // calculate the persistent thermal velocities + // doesn't actually use any resources, but stores this for all 30 atomic numbers + _vthermv.resize(numAtoms, 0.); + for (int Z = 1; Z <= numAtoms; Z++) _vthermv[Z - 1] = sqrt(Constants::k() * temperature() / AtomUtils::mass(Z)); + + // Photo-absorption + // calculate the parameters for the sigmoid function approximating the convolution with a Gaussian + // at the threshold energy for each cross section record, and store the result into a temporary vector; + // the information includes the thermal energy dispersion at the threshold energy and + // the intrinsic cross section at the threshold energy plus twice this energy dispersion + for (auto& upa : usedPar) + { + auto& ion = _ionParamv[upa.ionIndex]; + upa.Es = upa.Eth * vtherm(ion.Z) / Constants::c(); + upa.sigmamax = upa.photoAbsorbSection(upa.Eth + 2. * upa.Es); + } + + // Resonant scattering + // Calculate the total branching probability for each resonant Lyman transition. + // This is the probability that a new photon will be emitted after a 'resonant absorption' event. + // This value can be lower than 1 because low energy photons are ignored and thus not re-emitted. + // There is some code duplication since we do this later in the calculating of the persistent data. + if (resonantScattering()) + { + for (auto& uly : usedLyr) + { + // strore thermal velocity for convenience + uly.vth = M_SQRT2 * vtherm(uly.Z); + + // total branching probability + uly.sprob = 0.; + // add up all the probabilities for each current->lower branches + for (auto& b : lybResource) + { + if (b.Z == uly.Z && b.upper == uly.index) uly.sprob += b.prob; + if (uly.sprob == 1.) break; // speed up since branching matrix has a lot of 1s and 0s + } + } + } + + // ------------ calculate/store persistent data ------------ + + // The persistent data is the data that is needed beyond the setup (scattering) + // No changes should be made to the usedFlr, usedLyr, or lybr arrays after this point. + // The usedFlr and usedLyr need to be in the same order as the persistent params! + + // Fluorescence + // Store the Z, wavelength, and width of each fluorescence transition. + // These are needed when scattering photons after a photo-absorption event. + _fluorescenceParamv.resize(_numFluo); + for (int f = 0; f != _numFluo; ++f) + { + const auto& ufl = usedFlr[f]; + auto& flp = _fluorescenceParamv[f]; + + flp.Z = ufl.Z; + flp.lambda = wavelengthToFromEnergy(ufl.E); + flp.width = ufl.W / 2.; // convert from FWHM to HWHM + } + + // Resonant scattering + // Store the Z, Lyman index, wavelength, Voigt parameter and the cumulative branching for each resonant transition. + // These are needed to sample atom velocities and to determine the branch to scatter to. + if (resonantScattering()) + { + _lymanParamv.resize(_numLym); + for (int l = 0; l != _numLym; ++l) + { + const auto& uly = usedLyr[l]; + auto& lyp = _lymanParamv[l]; + + double vth = M_SQRT2 * vtherm(uly.Z); + + lyp.Z = uly.Z; + lyp.index = uly.index; + lyp.lambda = uly.lam; + lyp.a = uly.lamA / (4. * M_PI * vth); + + int Z = lyp.Z; + int upper = lyp.index; + + // branching probability + Array pLyl(0., upper + 1); // can only decay to lower Lyman index + for (auto& b : lybResource) + { + // for each Lyman Z, upper level, store all the lower probabilities + if (Z == b.Z && upper == b.upper) pLyl[b.lower] = b.prob; + } + + // store the cumulative + NR::cdf(lyp.cumbranchingv, pLyl); + } + } + } + + // ------------ wavelength grid ------------ + + // construct a wavelength grid for sampling cross sections containing a merged set of grid points + // in the relevant wavelength range (intersection of simulation range and nonzero range): + // - a fine grid in log space that provides sufficient resolution for most applications + // - all specific wavelengths mentioned in the configuration of the simulation (grids, normalizations, ...) + // ensuring that the cross sections are calculated at exactly these wavelengths + // - 7 extra wavelength points around the threshold energies for all transitions, + // placed at -2, -4/3, -2/3, 0, 2/3, 4/3, 2 times the thermal energy dispersion + + // we first gather all the wavelength points, in arbitrary order, and then sort them + vector lambdav; + lambdav.reserve(5 * numWavelengthsPerDex); + + // get the relevant range (intersection of simulation range and nonzero range) + Range range = config->simulationWavelengthRange(); + range.intersect(nonZeroRange); + + // add a fine grid in log space; + // use integer multiples as logarithmic grid points so that the grid is stable for changing wavelength ranges + constexpr double numPerDex = numWavelengthsPerDex; // converted to double to avoid casting + int minLambdaSerial = std::floor(numPerDex * log10(range.min())); + int maxLambdaSerial = std::ceil(numPerDex * log10(range.max())); + for (int k = minLambdaSerial; k <= maxLambdaSerial; ++k) lambdav.push_back(pow(10., k / numPerDex)); + + // add the wavelengths mentioned in the configuration of the simulation + for (double lambda : config->simulationWavelengths()) + if (range.contains(lambda)) lambdav.push_back(lambda); + + // add wavelength points around the threshold energies for all transitions + for (const auto& upa : usedPar) + { + double Es = upa.Es; + for (double delta : {-2., -4. / 3., -2. / 3., 0., 2. / 3., 4. / 3., 2.}) + { + double lambda = wavelengthToFromEnergy(upa.Eth + delta * Es); + if (range.contains(lambda)) lambdav.push_back(lambda); + } + } + + // add the fluorescence emission wavelengths + for (const auto& ufl : usedFlr) + { + double lambda = wavelengthToFromEnergy(ufl.E); + if (range.contains(lambda)) lambdav.push_back(lambda); + } + + // add the outer wavelengths of our nonzero range, plus an extra just outside of that range, + // so that there are always at least three points and thus two bins in the grid + lambdav.push_back(nonZeroRange.min()); + lambdav.push_back(nonZeroRange.max()); + lambdav.push_back(nonZeroRange.max() * 1.000001); // this wavelength point is never actually used + + // sort the wavelengths and remove duplicates + NR::unique(lambdav); + int numLambda = lambdav.size(); + + // derive a wavelength grid that will be used for converting a wavelength to an index in the above array; + // the grid points are shifted to the left of the actual sample points to approximate rounding + _lambdav.resize(numLambda); + _lambdav[0] = lambdav[0]; + for (int ell = 1; ell != numLambda; ++ell) + { + _lambdav[ell] = sqrt(lambdav[ell] * lambdav[ell - 1]); + } + + // ------------ extinction ------------ + + // calculate the extinction cross section at every wavelength; to guarantee that the cross section is zero + // for wavelengths outside our range, leave the values for the three outer wavelength points at zero + _sigmaextv.resize(numLambda); + for (int ell = 1; ell < numLambda - 2; ++ell) + { + double lambda = lambdav[ell]; + double sigma = 0.; + + // electron scattering + for (int i = 0; i < _numIons; i++) + { + const auto& ion = _ionParamv[i]; + sigma += _com->sectionSca(lambda, ion.Z) * _abundances[i]; + } + + // photo-absorption and fluorescence + for (const auto& upa : usedPar) + { + double E = wavelengthToFromEnergy(lambda); + sigma += upa.photoAbsorbThermalSection(E) * _abundances[upa.ionIndex]; + } + + // resonant scattering + for (const auto& uly : usedLyr) + { + sigma += uly.section(lambda) * _abundances[uly.ionIndex]; + } + + _sigmaextv[ell] = sigma; + } + + // ------------ scattering ------------ + + // make room for the scattering cross section and the cumulative fluorescence/scattering probabilities + _sigmascav.resize(numLambda); + _cumsigmascavv.resize(numLambda, 0); + + // provide temporary array for the non-normalized fluorescence/scattering contributions (at the current wavelength) + int numInteractions = _numIons + _numFluo + _numLym; + Array sections(numInteractions); + + // calculate the above for every wavelength; as before, leave the values for the outer wavelength points at zero + for (int ell = 1; ell < numLambda - 2; ++ell) + { + double lambda = lambdav[ell]; + double E = wavelengthToFromEnergy(lambda); + + // electron scattering + for (int i = 0; i < _numIons; i++) + { + const auto& ion = _ionParamv[i]; + + sections[i] = _com->sectionSca(lambda, ion.Z) * _abundances[i]; + } + + // fluorescence: iterate over both cross section and fluorescence parameter sets in sync + for (int f = 0; f < _numFluo; f++) + { + const auto& ufl = usedFlr[f]; + const auto& upa = usedPar[ufl.paIndex]; + + double section = upa.photoAbsorbThermalSection(E) * _abundances[upa.ionIndex] * ufl.omega; + sections[_numIons + f] = section; + } + + // resonant scattering + for (int l = 0; l < _numLym; l++) + { + const auto& uly = usedLyr[l]; + + double section = uly.section(lambda) * _abundances[uly.ionIndex] * uly.sprob; + sections[_numIons + _numFluo + l] = section; + } + + // determine the normalized cumulative probability distribution and the cross section + _sigmascav[ell] = NR::cdf(_cumsigmascavv[ell], sections); + } +} + +//////////////////////////////////////////////////////////////////// + +XRayIonicGasMix::~XRayIonicGasMix() +{ + delete _com; + if (resonantScattering()) delete _dpf; +} + +//////////////////////////////////////////////////////////////////// + +int XRayIonicGasMix::indexForLambda(double lambda) const +{ + return NR::locateFail(_lambdav, lambda); +} + +//////////////////////////////////////////////////////////////////// + +double XRayIonicGasMix::vtherm(int Z) const +{ + return _vthermv[Z - 1]; +} + +//////////////////////////////////////////////////////////////////// + +MaterialMix::MaterialType XRayIonicGasMix::materialType() const +{ + return MaterialMix::MaterialType::Gas; +} + +//////////////////////////////////////////////////////////////////// + +bool XRayIonicGasMix::hasPolarizedScattering() const +{ + return electronScattering() == ElectronScattering::FreeWithPolarization || resonantScattering(); +} + +//////////////////////////////////////////////////////////////////// + +bool XRayIonicGasMix::hasResonantScattering() const +{ + return resonantScattering(); +} + +//////////////////////////////////////////////////////////////////// + +bool XRayIonicGasMix::hasScatteringDispersion() const +{ + return true; +} + +//////////////////////////////////////////////////////////////////// + +bool XRayIonicGasMix::scatteringEmulatesSecondaryEmission() const +{ + return true; +} + +//////////////////////////////////////////////////////////////////// + +vector XRayIonicGasMix::specificStateVariableInfo() const +{ + return vector{StateVariable::numberDensity()}; +} + +//////////////////////////////////////////////////////////////////// + +double XRayIonicGasMix::mass() const +{ + return Constants::Mproton(); +} + +//////////////////////////////////////////////////////////////////// + +double XRayIonicGasMix::sectionAbs(double lambda) const +{ + int index = indexForLambda(lambda); + if (index < 0) return 0.; + return _sigmaextv[index] - _sigmascav[index]; +} + +//////////////////////////////////////////////////////////////////// + +double XRayIonicGasMix::sectionSca(double lambda) const +{ + int index = indexForLambda(lambda); + if (index < 0) return 0.; + return _sigmascav[index]; +} + +//////////////////////////////////////////////////////////////////// + +double XRayIonicGasMix::sectionExt(double lambda) const +{ + int index = indexForLambda(lambda); + if (index < 0) return 0.; + return _sigmaextv[index]; +} + +//////////////////////////////////////////////////////////////////// + +double XRayIonicGasMix::opacityAbs(double lambda, const MaterialState* state, const PhotonPacket* /*pp*/) const +{ + double number = state->numberDensity(); + return number > 0. ? sectionAbs(lambda) * number : 0.; +} + +//////////////////////////////////////////////////////////////////// + +double XRayIonicGasMix::opacitySca(double lambda, const MaterialState* state, const PhotonPacket* /*pp*/) const +{ + double number = state->numberDensity(); + return number > 0. ? sectionSca(lambda) * number : 0.; +} + +//////////////////////////////////////////////////////////////////// + +double XRayIonicGasMix::opacityExt(double lambda, const MaterialState* state, const PhotonPacket* /*pp*/) const +{ + double number = state->numberDensity(); + return number > 0. ? sectionExt(lambda) * number : 0.; +} + +//////////////////////////////////////////////////////////////////// + +void XRayIonicGasMix::setScatteringInfoIfNeeded(PhotonPacket* pp, const MaterialState* state, const double lambda) const +{ + auto scatinfo = pp->getScatteringInfo(); + + if (!scatinfo->valid) + { + scatinfo->valid = true; + + // indexForLambda should never be able to fail here, check needed anyway? + scatinfo->species = NR::locateClip(_cumsigmascavv[indexForLambda(lambda)], random()->uniform()); + + // Compton scattering + if (scatinfo->species < _numIons) + { + int i = scatinfo->species; + const auto& ion = _ionParamv[i]; + scatinfo->velocity = vtherm(ion.Z) * random()->maxwell(); + } + // Fluorescenct emission (scattering) + else if (scatinfo->species < _numIons + _numFluo) + { + int f = scatinfo->species - _numIons; + auto flp = _fluorescenceParamv[f]; + double lambda = flp.lambda; + double width = flp.width; + + scatinfo->velocity = vtherm(flp.Z) * random()->maxwell(); + + if (width == 0.) + { + // for a zero-width line, simply copy the central wavelength + scatinfo->lambda = lambda; + } + else + { + // otherwise sample a wavelength from the Lorentz line shape in energy space; + // the tails of the Lorentz distribition are very long, occasionaly resulting in negative energies; + // therefore we loop until the sampled wavelength is meaningful + while (true) + { + double center = wavelengthToFromEnergy(lambda); + scatinfo->lambda = wavelengthToFromEnergy(center + width * random()->lorentz()); + if (nonZeroRange.contains(scatinfo->lambda)) break; + } + } + } + // Resonant Lyman scattering + else + { + int l = scatinfo->species - _numIons - _numFluo; + const auto& ulyp = _lymanParamv[l]; // upper Lyman + + int upper = ulyp.index; + int lower = NR::locateFail(ulyp.cumbranchingv, random()->uniform()); + + if (lower == -1) throw FATALERROR("Sampling from Lyman branching probability has failed"); + + double center; + double a; + double vth; + bool J32; // true if Ju=3/2 -> happens at odd Lyman index + + // if coherent (no branching) + if (lower == upper) + { + vth = M_SQRT2 * vtherm(ulyp.Z); + a = ulyp.a; + center = ulyp.lambda; + J32 = upper % 2 == 1; + + scatinfo->lambda = 0.; // explicitly don't use + } + // if incoherent (branching) + else + { + int index = l - (upper - lower); // index of the lower branching + if (index < 0 || index >= _numLym) throw FATALERROR("upper/lower index out of range"); + + const auto& llyp = _lymanParamv[index]; // lower Lyman + + // set parameters to those of the lower branching + vth = M_SQRT2 * vtherm(llyp.Z); + a = llyp.a; + center = llyp.lambda; + J32 = false; // branching is isotropic + + scatinfo->lambda = llyp.lambda; + } + + // if J32 -> Lya1, Lyb1, ... -> 50/50 dipole/isotropic + // if J12 -> Lya2, Lyb2, ... -> 100 isotropic + scatinfo->dipole = J32 ? random()->uniform() < 0.5 : false; + + // sample a atom velocity from Voigt profile + scatinfo->velocity = LyUtils::sampleAtomVelocity( + lambda, center, vth, a, temperature(), state->numberDensity(), pp->direction(), config(), random()); + } + } +} + +//////////////////////////////////////////////////////////////////// + +bool XRayIonicGasMix::peeloffScattering(double& I, double& Q, double& U, double& V, double& lambda, Direction bfkobs, + Direction bfky, const MaterialState* state, const PhotonPacket* pp) const +{ + // draw a random scattering channel and atom velocity, unless a previous peel-off stored this already + setScatteringInfoIfNeeded(const_cast(pp), state, lambda); + auto scatinfo = const_cast(pp)->getScatteringInfo(); + + // Compton scattering in electron rest frame; with support for polarization if enabled + if (scatinfo->species < _numIons) + { + int i = scatinfo->species - _numIons; + const auto& ion = _ionParamv[i]; + // transform the wavelength into the rest frame of the electron + lambda = PhotonPacket::shiftedReceptionWavelength(lambda, pp->direction(), scatinfo->velocity); + _com->peeloffScattering(I, Q, U, V, lambda, ion.Z, pp->direction(), bfkobs, bfky, pp); + lambda = PhotonPacket::shiftedEmissionWavelength(lambda, bfkobs, scatinfo->velocity); + return false; + } + + // fluorescence + else if (scatinfo->species < _numIons + _numFluo) + { + // unpolarized isotropic emission; the bias weight is trivially 1 and there is no contribution to Q, U, V + I = 1.; + + // update the photon packet wavelength to the (possibly sampled) wavelength of this fluorescence transition + lambda = scatinfo->lambda; + lambda = PhotonPacket::shiftedEmissionWavelength(lambda, bfkobs, scatinfo->velocity); + return true; + } + + // resonant scattering + else + { + if (scatinfo->dipole) + _dpf->peeloffScattering(I, Q, U, V, pp->direction(), bfkobs, bfky, pp); + else + // isotropic scattering removes polarization, so the contribution is trivially 1 + I = 1.; + + // if coherent (no branching) + if (scatinfo->lambda == 0.) + lambda = LyUtils::shiftWavelength(lambda, scatinfo->velocity, pp->direction(), bfkobs); + // if incoherent (branching) + else + lambda = PhotonPacket::shiftedEmissionWavelength(scatinfo->lambda, bfkobs, scatinfo->velocity); + + return false; + } +} + +//////////////////////////////////////////////////////////////////// + +void XRayIonicGasMix::performScattering(double lambda, const MaterialState* state, PhotonPacket* pp) const +{ + // draw a random fluorescence channel and atom velocity, unless a previous peel-off stored this already + setScatteringInfoIfNeeded(const_cast(pp), state, lambda); + auto scatinfo = const_cast(pp)->getScatteringInfo(); + + // room for the outgoing direction + Direction bfknew; + + // Compton scattering, with support for polarization if enabled: + // determine the new propagation direction and wavelength, and if polarized, update the stokes vector + if (scatinfo->species < _numIons) + { + int i = scatinfo->species - _numIons; + const auto& ion = _ionParamv[i]; + lambda = PhotonPacket::shiftedReceptionWavelength(lambda, pp->direction(), scatinfo->velocity); + bfknew = _com->performScattering(lambda, ion.Z, pp->direction(), pp); + lambda = PhotonPacket::shiftedEmissionWavelength(lambda, bfknew, scatinfo->velocity); + } + + // fluorescence, always unpolarized and isotropic + else if (scatinfo->species < _numIons + _numFluo) + { + // clear the stokes vector (only relevant if polarization support is enabled) + pp->setUnpolarized(); + + // update the photon packet wavelength to the (possibly sampled) wavelength of this fluorescence transition + lambda = scatinfo->lambda; + bfknew = random()->direction(); + lambda = PhotonPacket::shiftedEmissionWavelength(lambda, bfknew, scatinfo->velocity); + + // indicate that this packet emulates secondary emission; + pp->setEmulatedSecondaryOrigin(state->mediumIndex()); + } + + // resonant scattering + else + { + if (scatinfo->dipole) + { + bfknew = _dpf->performScattering(pp->direction(), pp); + } + else + { + bfknew = random()->direction(); + pp->setUnpolarized(); + } + + // if coherent (no branching) + if (scatinfo->lambda == 0.) + lambda = LyUtils::shiftWavelength(lambda, scatinfo->velocity, pp->direction(), bfknew); + // if incoherent (branching) + else + lambda = PhotonPacket::shiftedEmissionWavelength(scatinfo->lambda, bfknew, scatinfo->velocity); + } + + // execute the scattering event in the photon packet + pp->scatter(bfknew, state->bulkVelocity(), lambda); +} + +//////////////////////////////////////////////////////////////////// + +double XRayIonicGasMix::indicativeTemperature(const MaterialState* /*state*/, const Array& /*Jv*/) const +{ + return temperature(); +} + +//////////////////////////////////////////////////////////////////// diff --git a/SKIRT/core/XRayIonicGasMix.hpp b/SKIRT/core/XRayIonicGasMix.hpp new file mode 100644 index 00000000..723b25c2 --- /dev/null +++ b/SKIRT/core/XRayIonicGasMix.hpp @@ -0,0 +1,497 @@ +/*////////////////////////////////////////////////////////////////// +//// The SKIRT project -- advanced radiative transfer //// +//// © Astronomical Observatory, Ghent University //// +///////////////////////////////////////////////////////////////// */ + +#ifndef XRAYIONICGASMIX_HPP +#define XRAYIONICGASMIX_HPP + +#include "ArrayTable.hpp" +#include "DipolePhaseFunction.hpp" +#include "MaterialMix.hpp" +#include "PhotonPacket.hpp" + +//////////////////////////////////////////////////////////////////// + +/** The XRayIonicGasMix class describes the material properties of a partially ionised gas in the + X-ray wavelength range, taking into account the effects of photo-absorption, fluorescence, + Lyman recombination, electron scattering, and optionally resonant Lyman resonant scattering. It + is largely an extension of the XRayAtomicGasMix class, supporting all ions instead of just the + neutral atoms. To avoid the use of this material mix outside of the regime for which it has + been designed, all cross sections are forced to zero below 4.3 eV and above 500 keV + (corresponding approximately to a wavelength range from 2.5 pm to 290 nm). + + The class assumes a gas containing a mixture of (ionized) elements with atomic numbers ranging + from 1 (hydrogen) to 30 (zinc). Including fully ionized ions, this results in a total of 495 + possible ions. Not all ions have to be used, the user can specify what ions are included with + the \em ions property. The spatial density distribution of the gas is determined by the product + of the number density and the abundances of the ions in use. The temperature of the gas can + also be configured by the user as a constant property. In other words, the abundances and the + temperature are considered to be spatially constant (for a given medium component), while the + overall density can obviously vary across space as usual. + + Photo-absorption and fluorescence + + Photo-absorption by an ion is the process where the energy of a photon is used to liberate a + bound electron from one of the electron shells of the atom. This class supports + photo-absorption by any of the 495 ions in the gas for any of the electron shells, i.e. up to + the K, L, M, or N shell depending on the ion. + + Fluorescence (in this context) is the process where an electron from a higher energy level + "falls" into an empty space created by photo-absorption, emitting a new photon with a different + energy. For each electron shell and for each possible fluorescence transition towards that + shell, the \em yield defines the probability that such fluorescence event occurs after an + electron has been liberated in that shell. This class supports the following fluorescent lines, + all with energies above \f$4.3 \, \text{eV}\f$: K\f$_{\alpha2}\f$, K\f$_{\alpha1}\f$, + K\f$_{\beta3}\f$, K\f$_{\beta1}\f$, L\f$_{\beta4}\f$, L\f$_{\beta3}\f$, L\f$_{1,2}\f$, + L\f$_{1,3}\f$, L\f$_{\eta}\f$, L\f$_{l}\f$, L\f$_{\gamma5}\f$, L\f$_{\beta6}\f$, + L\f$_{\beta1}\f$, L\f$_{\alpha2}\f$, L\f$_{\alpha1}\f$, M\f$_{1,2}\f$, M\f$_{1,3}\f$, + M\f$_{2,4}\f$, M\f$_{3,4}\f$, M\f$_{3,5}\f$ M\f$_{2}\f$N\f$_{1}\f$, M\f$_{3}\f$N\f$_{1}\f$ + (these transitions are not all the same as the ones in the XRayAtomicGasMix class!). + + Unlike the XRayAtomicGasMix, the fluorescent lines never have an intrinsic shape. + + Because fluorescence only occurs as the result of a photo-absorption event, this class + implements fluorescence as a form of scattering (where the wavelength of the photon being + scattered changes) as opposed to emission. This allows both photo-absorption and fluorescence + to be treated during primary emission. A possible drawback is that the weaker fluorescence + lines will be represented by a fairly small number of photon packets. + + Lyman recombination + + For hydrogen-like ions, photo-absorption cannot be followed by fluorescence since the ion will + have no electrons left. Instead it can be followed by a radiative recombination event either to + the ground state or to an excited level, which will then cascade down to the ground state. This + class assumes instantaneous recombination and only implements the Lyman series photons produced + from cascades from the level \f$i\f$ to the ground state. This will thus only result in line + emission and no continuum emission. This allows for this process to be modelled the same way as + fluorescence. + + Scattering by bound electrons + + Electrons bound to atoms or free in the gas scatter incoming X-ray photons. Currently this + class can only model electrons as free using inelastic (Compton) scattering. We completely + ignore elastic (Rayleigh) scattering and inelastic (bound-Compton) scattering. The user can + select one of three implementations for electron scattering, these are: + + - \em None: ignore scattering by all electrons. + + - \em Free: use free-electron Compton scattering for all electrons. This assumes all electrons + (bound or free) scatter using the free-electron Compton scattering cross section. + + - \em FreeWithPolarization: use free-electron Compton scattering for all electrons with support + for polarization. + + This means that all electrons are treated the same way. A material mix with Fe+0 will thus have + the same scattering cross section as a material mix with Fe+26. + + Lyman resonant scattering + + Lyman resonant scattering is the process where a photon is absorbed and re-emitted by a bound + electron in an ion, promoting it from the ground state (\f$n=1\f$) to an excited level (\f$n + \leq 10\f$) and back again. This class implements these transitions for all hydrogen-like ions + up to atomic number 30. There are 9 Lyman-series lines included, each with fine-structure + splitting, resulting in a total of 18 transitions to the ground state (e.g., Ly\f${\alpha1}\f$, + Ly\f${\alpha2}\f$, Ly\f$_{\beta1}\f$, etc.). + + Lyman-series transitions may occur either as a direct transition from and to the ground state + or via a radiative cascade. Direct transitions correspond to coherent scattering events, while + cascades are incoherent and erase information about the initial photon packet. when a cascade + occurs, we ignore all intermediate (non-Lyman) radiative decays, keeping only the final (Lyman) + transition to the ground state. This means we can also model the cascade as a single scattering + event. The branching data, determining the probability of an (in)coherent transition, + is obtained from the SPEX database Kaastra et al. (2024). + + Configuring the simulation + + In addition to a medium component configured with the material mix represented by this class, + simulations will usually include primary sources and possibly a dust medium. There is, however, + no need to include secondary emission, so the simulation mode can be set to "ExtinctionOnly" + and there is no need to store the radiation field. Only if the user option \em + resonantScattering is enabled does the simulation mode need to be set to "LyaExtinctionOnly". + The resulting continuum spectrum and absorption and emission features can be recorded by a + single instrument configured with a high-resolution wavelength grid, or separate instruments + can be configured with wavelength grids to resolve specific features of interest. + + Although there is no secondary emission phase, fluorescent emission and lyman recombination, + which are modelled as scattering, do contribute to the secondary emission spectrum. To study + these processes seperately from the background continuum, configure the instrument to record + flux components seperately and consider the secondary emission component, which includes these + lines in addition to any flux scattered from them. To ensure that low-intensity lines are + properly included in this flux, set the advanced property \em minWeightReduction in the \c + PhotonPacketOptions section to a value of \c 1e10 or so. With the default value of \c 1e4, + low-intensity line photon packets are killed before having a chance to register in the + instruments. + + The input model must define the spatial distribution of the number density. The use can then + define the abundances of the ions in the gas in this material mix. The spatial density is + simply the product of the number density and the abundances of the ions in use. + + If this material mix is associated with a subclass of ImportedMedium, the spatial density + distribution is read from an input file. In that case, the ski file attributes \em + importMetallicity, \em importTemperature, and \em importVariableMixParams must be left at + 'false'. For example, if bulk velocities are also imported for this medium component (i.e. \em + importVelocity is 'true'), the column order would be \f[ ..., n, v_\mathrm{x}, v_\mathrm{y}, + v_\mathrm{z} \f] + + The relative abundances of the present ions in the gas and the temperature of the gas are + configured in the ski file as constant properties. In other words, the abundances and the + temperature are considered to be spatially constant (for a given medium component). + + Photo-absorption cross section + + The total photo-absorption cross section per hydrogen atom for this material mix is obtained by + accumulating the photo-absorption cross sections for all shells and for all individual + elements, weighted by element abundance, and convolved with a Gaussian profile reflecting each + element's thermal velocity. Because the abundances and the temperature are fixed, this + calculation can be performed during setup and the result stored, discretized on a + high-resolution wavelength grid for later retrieval. + + Verner and Yakovlev (1995, 1996) provide analytic fits to the photo-absorption cross sections + \f$\sigma_{ph}(E)\f$ as a function of photon energy \f$E\f$ for all the ions up to Z=30: + + \f[\begin{aligned} \sigma_{ph}(E) &= \begin{cases} 0 & E < E_\mathrm{th} \\ \sigma_0 \, F(y) & + E_\mathrm{th} \le E < E_\mathrm{max}\\ 0 & E_\mathrm{max} \le E \end{cases}, \\ y &= + \frac{E}{E_0}\\ F(y) &= \left[(y-1)^2+y_{\rm w}^2 \right]y^{-Q} \left(1+ \sqrt{(y/y_{\rm a})} + \right )^{-P}, \\ Q&=5.5+l-0.5P, \end{aligned} \f] with \f$E_\mathrm{th}\f$ the tabulated + ionization threshold energy, \f$E_\mathrm{max}\f$ the constant maximum energy (500 keV) for the + formula to be valid, \f$\sigma_0\f$, \f$E_0\f$, \f$y_{\rm w}\f$, \f$y_{\rm a}\f$ and \f$P\f$ + five tabulated fitting parameters, and \f$l\f$ the subshell orbital quantum number (\f$l=0, 1, + 2, 3\f$ for s, p, d, f orbitals respectively). + + Fluorescence cross section + + The total fluorescence cross section per hydrogen atom for this material mix is obtained + similarly, but now including only the K, L and M shell photo-absorption cross sections for each + ion, multiplied by the appropriate fluorescence yields in addition to the element abundance. + These yields are obtained from Kaastra & Mewe (1993) for all ions up to Z=30, from neutral down + to B-like. + + Lyman recombination cross section + + The total Lyman recombination cross section per hydrogen atom for this material mix is obtained + the exact same was as for fluorescence. The yields are now temperature-dependent however, but + since the temperature is assumed to be constant this is simply read during the setup. The + temperature-dependent yields are obtained from Mao & Kaastra (2016). + + Electron scattering - cross section and phase function + + As described above, this class provides a free-electron Compton scattering implementation. Note + that the cross section must be multiplied by the abundance of the corresponding ion. Since we + assume all electrons to be free, the cross section for an ion with atomic number \f$Z\f$ is + simply given by \f$Z\,\sigma_\mathrm{C}\f$, where \f$\sigma_\mathrm{C}\f$ is the + (wavelength-dependent) Compton cross section for a single free electron. The implementation of + the scattering events is delegated to the ComptonPhaseFunction class; see there for more + information on the cross section and phase function for free-electron Compton scattering. + + Electron scattering - photon energy shift + + Compton scattering is inelastic, meaning that the photon transfers a fraction of its energy to + the electron involved in the interaction. The energy shift \f$E'/E\f$ is given by the Compton + factor \f$C(\theta, E)\f$ defined above. The implementation is delegated to the + ComptonPhaseFunction class. + + Lyman resonant scattering - cross section and phase function + + The resonant Lyman transitions are broadened both thermally and intrinsically. The resulting + cross section is described by a Voigt profile, the implementation of this is delegated to the + \em LyUtils and \em VoigtProfile namespaces. The data for the line energies and Einstein A + coefficients are obtained from the SPEX database Kaastra et al. (2024). + + If the resonant scattering occurs for the fine-structure levels Ly\f$_{\alpha1}\f$, + Ly\f$_{\beta1}\f$, Ly\f$_{\gamma1}\f$, etc. the resulting photon has a 50% chance to be emitted + isotropically and 50% chance for a dipole phase function. This dipole phase function always + includes polarization, the \em ElectronScattering property does not affect it in any way. The + implementation of this phase function can be found in the \em DipolePhaseFunction class. For + the other fine-structure levels Ly\f$_{\alpha2}\f$, Ly\f$_{\beta2}\f$, etc. the resulting + photon is always emitted isotropically. + + If the scattering is coherent (no-cascade) the wavelength shift is determined by: \f[ \lambda' + = \lambda \frac{(1 - \boldsymbol{k}_\mathrm{out} \cdot \boldsymbol{v}_\mathrm{atom} / c)}{(1 - + \boldsymbol{k}_\mathrm{in} \cdot \boldsymbol{v}_\mathrm{atom} / c)} \f] otherwise (cascade) the + wavelength shift is only determined by the outgoing direction: \f[ \lambda' = \lambda (1 - + \boldsymbol{k}_\mathrm{out} \cdot \boldsymbol{v}_\mathrm{atom} / c) \f] + + Performing scattering + + The function performing an actual scattering event randomly selects one of the supported + scattering channels (i.e. scattering by an electron, fluorescent line emission following a + photo-absorption event or scattering by a resonant Lyman line). The relative probabilities for + these transitions as a function of incoming photon packet wavelength are also calculated during + setup. The selected transition determines the scattering mechanism. For electrons, Compton + scattering is used. For fluorescence, the emission direction is isotropic, and the outgoing + wavelength is the fluorescence wavelength. For the resonant Lyman lines, the direction can be + either isotropic and unpolarized or following a dipole phase function for both direction and + polarization. The outcome is determined by the Lyman index... lya1, lyb1, etc. + + Thermal dispersion + + The thermal dispersion appropriate for a given interaction depends on the (constant) + temperature configured for this material mix and on the mass of the interacting atom. For + electron scattering and fluorescence events, the implementation is straightforward. Once a + channel has been randomly selected, the magnitude of the interacting atom's thermal velocity is + taken from a precalculated table and a velocity vector is sampled from the Maxwell + distribution. + + For photo-absorption, the situation is much more involved. In principle, the full cross section + curve for each ionization transition must be convolved with a Gaussian kernel of appropriate + width. In practice, the cross section curves are very smooth except for the step at the + threshold energy. The effect of the convolution is therefore limited to energies near the + threshold energy, replacing the infinitely steep step by a sigmoid function. Considering that + the convolution of a step function with a Gaussion is given by the error function, this class + uses the following approximation. With \f$E_\mathrm{th}\f$ the threshold energy and + \f$E_\mathrm{s}\f$ the energy dispersion corresponding to the thermal velocity of the atom + being ionized, the cross section near \f$E_\mathrm{th}\f$ is replaced by + + \f[ \sigma_{ph}'(E) = \frac{\sigma_{ph}(E_\mathrm{th} +2 E_\mathrm{s})} {2} \, \left[ 1+ + \mathrm{erf} ( \frac{E - E_\mathrm{th}} {E_\mathrm{s}} ) \right] \qquad \mathrm{for} \; + E_\mathrm{th} -2 E_\mathrm{s} < E < E_\mathrm{th} +2 E_\mathrm{s}. \f] + + The sigmoid is scaled by the value of the actual cross section at the end of the interval, + achieving good continuity at that point in most cases, including all important transitions. + + The thermal broadening present in the resonant scattering is already built into the Voigt + profile, its implementation can be found in the \em LyUtils and \em VoigtProfile namespaces. + */ +class XRayIonicGasMix : public MaterialMix +{ + /** The enumeration type indicating the implementation used for scattering by electrons. */ + ENUM_DEF(ElectronScattering, None, Free, FreeWithPolarization) + ENUM_VAL(ElectronScattering, None, "ignore electron") + ENUM_VAL(ElectronScattering, Free, "use free-electron Compton scattering for all electrons") + ENUM_VAL(ElectronScattering, FreeWithPolarization, + "use free-electron Compton scattering with support for polarization") + ENUM_END() + + ITEM_CONCRETE(XRayIonicGasMix, MaterialMix, + "An ionised gas mix supporting photo-absorption, fluorescence and resonant Lyman scattering for " + "X-ray wavelengths") + ATTRIBUTE_TYPE_INSERT(XRayIonicGasMix, "GasMix") + + PROPERTY_STRING(ions, "the names of the ions for each element (e.g. H,He+,Fe+25,..)") + + PROPERTY_DOUBLE_LIST(abundances, "the abundances of the ions in the same order as the ions property") + + PROPERTY_DOUBLE(temperature, "the temperature of the gas in K") + ATTRIBUTE_QUANTITY(temperature, "temperature") + ATTRIBUTE_MIN_VALUE(temperature, "[3") + ATTRIBUTE_MAX_VALUE(temperature, "1e9]") + ATTRIBUTE_DEFAULT_VALUE(temperature, "1e4") + ATTRIBUTE_DISPLAYED_IF(temperature, "Level2") + + PROPERTY_ENUM(electronScattering, ElectronScattering, "implementation of scattering by electrons") + ATTRIBUTE_DEFAULT_VALUE(electronScattering, "Free") + ATTRIBUTE_DISPLAYED_IF(electronScattering, "Level3") + + PROPERTY_BOOL(resonantScattering, "enable Lyman resonant scattering for all hydrogen-like ions") + ATTRIBUTE_DEFAULT_VALUE(resonantScattering, "false") + ATTRIBUTE_DISPLAYED_IF(resonantScattering, "Level2") + ATTRIBUTE_RELEVANT_IF(resonantScattering, "simulationModeLyaExtinctionOnly") + + ITEM_END() + + //============= Construction - Setup - Destruction ============= + +public: + /** This function precalculates relevant cross sections and relative contributions over a + high-resolution wavelength grid. It also stores presistent data that is used during the + simulation to perform scattering. */ + void setupSelfBefore() override; + + /** The destructor destructs the phase function helpers that were created during setup. */ + ~XRayIonicGasMix(); + + //======== Private support functions ======= + +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. */ + int indexForLambda(double lambda) const; + + /** This function returns the precalculated thermal velocity of the atom with the specified + atomic number \f$Z\f$. */ + double vtherm(int Z) const; + + //============= Capabilities ============= + +public: + /** This function returns the fundamental material type represented by this material mix, which + is MaterialType::Gas. */ + MaterialType materialType() const override; + + /** This function returns true if this material mix supports polarization during scattering + events. For this class, the function returns true if the \em scatterBoundElectrons property + has been set to \c FreeWithPolarization or is the \em resonantScattering property is set to + true. */ + bool hasPolarizedScattering() const override; + + /** This function returns true if scattering for this material mix is resonant. For this class, + the function returns true if the \em resonantScattering property is set to true. */ + bool hasResonantScattering() const override; + + /** This function returns true, indicating that a scattering interaction for this material mix + may (and usually does) adjust the wavelength of the interacting photon packet. */ + bool hasScatteringDispersion() const override; + + /** This function returns true, indicating that a scattering interaction for this material mix + may emulate secondary emission. This is used to implement fluorescence as scattering. */ + bool scatteringEmulatesSecondaryEmission() const override; + + //============= Medium state setup ============= + + /** This function returns a list of StateVariable objects describing the specific state + variables used by the receiving material mix. For this class, the function returns just the + descriptor for the number density. */ + vector specificStateVariableInfo() const override; + + //============= Low-level material properties ============= + + /** This function returns the mass of a hydrogen atom. */ + double mass() const override; + + /** This function returns the absorption (i.e. \f$\sigma_\mathrm{ext} - \sigma_\mathrm{sca}\f$) + cross section per hydrogen atom at the given wavelength and using the abundances and + temperature configured for this material mix. */ + double sectionAbs(double lambda) const override; + + /** This function returns the scattering (i.e. electron scattering, fluorescence, resonant + scattering) cross section per hydrogen atom at the given wavelength and using the + abundances and temperature configured for this material mix. */ + double sectionSca(double lambda) const override; + + /** This function returns the extinction cross section per hydrogen atom at the given + wavelength and using the abundances and temperature configured for this material mix. */ + double sectionExt(double lambda) const override; + + //============= High-level photon life cycle ============= + + /** This function returns the absorption (i.e. \f$\sigma_\mathrm{ext} - \sigma_\mathrm{sca}\f$) + opacity at the given wavelength and material state, using the abundances and temperature + configured for this material mix. The photon packet properties are not used. */ + double opacityAbs(double lambda, const MaterialState* state, const PhotonPacket* pp) const override; + + /** This function returns the scattering (i.e. electron scattering, fluorescence, resonant + scattering) opacity at the given wavelength and material state, using the abundances and + temperature configured for this material mix. The photon packet properties are not used. */ + double opacitySca(double lambda, const MaterialState* state, const PhotonPacket* pp) const override; + + /** This function returns the extinction opacity at the given wavelength and material state, + using the abundances and temperature configured for this material mix. The photon packet + properties are not used. */ + double opacityExt(double lambda, const MaterialState* state, const PhotonPacket* pp) const override; + +private: + /** This private function draws a random scattering channel and atom velocity and stores this + information in the photon packet's scattering information record, unless a previous + peel-off stored this already. For fluorescence transitions that support a nonzero line + width, the function also draws a random wavelength from the line shape. For a resonant + scattering event it determines whether a cascade occrus and if the phase function is + isotropic or a dipole. */ + 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 + 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. If there is + wavelength shift, the new wavelength value replaces the incoming value of the \em lambda + argument. + + The function first calls the private setScatteringInfoIfNeeded() function to establish a + random scattering channel and atom velocity for this event. In case the selected channel is + electron scattering, the peel-off bias weight and wavelength shift are determined by the + Compton phase function and the selected atom velocity. For fluorescence, the peel-off bias + weight is trivially one because emission is isotropic and unpolarized, and the outgoing + wavelength is determined by Doppler-shifting the rest wavelength of the selected + fluorescence transition for the selected atom velocity. For resonant scattering the + peel-off bias weight and polarization are determined by its isotropic or dipole phase + function, and the wavelength is shifted accordingly depending on whether the scattering is + coherent or through a cascade. */ + 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 + cell and medium component represented by the specified material state and the receiving + material mix. It first calls the private setScatteringInfoIfNeeded() function to establish + a random scattering channel and atom velocity for this event. + + In case the selected channel is electron scattering, the outgoing direction and adjusted + wavelength are determined by the Compton phase function and the selected atom velocity. For + fluorescence, emission is isotropic, so the outgoing direction is randomly chosen from the + isotropic distribution. The outgoing wavelength is determined by Doppler-shifting the rest + wavelength of the selected fluorescence transition for the selected atom velocity. For + resonant scattering, the emission direction and polarization are determined by its + isotropic or dipole phase function. The wavelength is shifted depending on whether the + scattering is coherent or occurs through a cascade. Coherent scattering depends on both the + incoming and outgoing directions, whereas a cascade removes the incoming photon information + and depends only on the outgoing direction. */ + void performScattering(double lambda, const MaterialState* state, PhotonPacket* pp) const override; + + //======== Temperature ======= + +public: + /** This function returns an indicative temperature of the material mix when it would be + embedded in a given radiation field. The implementation in this class ignores the radiation + field and returns the (spatially constant) temperature configured for this material mix. */ + double indicativeTemperature(const MaterialState* state, const Array& Jv) const override; + + //======================== Data Members ======================== + +public: + // base class for electron scattering helpers (public because we derive from it in anonymous namespace) + class ScatteringHelper; + +private: + // all of the data below is calculated in the setupSelfBefore(), but is persistent for use after setup to perform scattering + + struct IonParam + { + IonParam(int Z, int N) : Z(Z), N(N) {} + + int Z; // atomic number + int N; // number of electrons + }; + // Compton scattering -> IonParam + // Fluorescence (+Lyman RC) -> FluorescenceParam + struct FluorescenceParam + { + int Z; // atomic number + double lambda; // wavelength (m) + double width; // width (eV) + }; + // Resonant scattering -> LymanParam + struct LymanParam + { + int Z; // atomic number + int index; // Lyman index (alpha1/2, alpha3/2, beta1/2, ...) + double lambda; // wavelength (m) + double a; // Voigt parameter + Array cumbranchingv; // normalized cumulative branching + }; + + int _numIons; // number of ions + int _numFluo; // number of fluorescence (+Lyman RC) transitions + int _numLym; // number of Lyman resonant scattering transitions + + // persistent data for scattering + vector _ionParamv; + vector _fluorescenceParamv; + vector _lymanParamv; + vector _vthermv; // indexed on Z + + // wavelength grid (shifted to the left of the actually sampled points to approximate rounding) + Array _lambdav; + // cross sections + Array _sigmaextv; // indexed on lambdav + Array _sigmascav; // indexed on lambdav + ArrayTable<2> _cumsigmascavv; // indexed on lambdav, interaction (electron + fluorescence + resonant) + + // compton-electron scattering helpers depending on the configured implementation + ScatteringHelper* _com{nullptr}; // Compton scattering helper + // dipole phase function for resonant scattering + DipolePhaseFunction* _dpf{nullptr}; +}; + +#endif diff --git a/SKIRT/faddeeva/CMakeLists.txt b/SKIRT/faddeeva/CMakeLists.txt new file mode 100644 index 00000000..f4bc3b6e --- /dev/null +++ b/SKIRT/faddeeva/CMakeLists.txt @@ -0,0 +1,26 @@ +# ////////////////////////////////////////////////////////////////// +# /// The SKIRT project -- advanced radiative transfer /// +# /// © Astronomical Observatory, Ghent University /// +# ////////////////////////////////////////////////////////////////// + +# ------------------------------------------------------------------ +# Calculates the Faddeeva function. +# The package is taken from http://ab-initio.mit.edu/faddeeva// +# This is specifically used to evaluate the Voigt profile by setting H(x, a) = Real(w(x+ia)). +#--------------------------------------------------------------------- + +# set the target name +set(TARGET faddeeva) + +# list the source files in this directory +file(GLOB SOURCES "*.cc") +file(GLOB HEADERS "*.hh") + +# create the library target +add_library(${TARGET} STATIC ${SOURCES} ${HEADERS}) + +# adjust C++ compiler flags to our needs +set(NO_EXTRA_WARNINGS true) # to avoid warnings in the faddeeva code (warnings still present) +include("../../SMILE/build/CompilerFlags.cmake") + + diff --git a/SKIRT/faddeeva/Faddeeva.cc b/SKIRT/faddeeva/Faddeeva.cc new file mode 100644 index 00000000..d553facc --- /dev/null +++ b/SKIRT/faddeeva/Faddeeva.cc @@ -0,0 +1,4585 @@ +// -*- mode:c++; tab-width:2; indent-tabs-mode:nil; -*- + +/* Copyright (c) 2012 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* (Note that this file can be compiled with either C++, in which + case it uses C++ std::complex, or C, in which case it + uses C99 double complex.) */ + +/* Available at: http://ab-initio.mit.edu/Faddeeva + + Computes various error functions (erf, erfc, erfi, erfcx), + including the Dawson integral, in the complex plane, based + on algorithms for the computation of the Faddeeva function + w(z) = exp(-z^2) * erfc(-i*z). + Given w(z), the error functions are mostly straightforward + to compute, except for certain regions where we have to + switch to Taylor expansions to avoid cancellation errors + [e.g. near the origin for erf(z)]. + + To compute the Faddeeva function, we use a combination of two + algorithms: + + For sufficiently large |z|, we use a continued-fraction expansion + for w(z) similar to those described in: + + Walter Gautschi, "Efficient computation of the complex error + function," SIAM J. Numer. Anal. 7(1), pp. 187-198 (1970) + + G. P. M. Poppe and C. M. J. Wijers, "More efficient computation + of the complex error function," ACM Trans. Math. Soft. 16(1), + pp. 38-46 (1990). + + Unlike those papers, however, we switch to a completely different + algorithm for smaller |z|: + + Mofreh R. Zaghloul and Ahmed N. Ali, "Algorithm 916: Computing the + Faddeyeva and Voigt Functions," ACM Trans. Math. Soft. 38(2), 15 + (2011). + + (I initially used this algorithm for all z, but it turned out to be + significantly slower than the continued-fraction expansion for + larger |z|. On the other hand, it is competitive for smaller |z|, + and is significantly more accurate than the Poppe & Wijers code + in some regions, e.g. in the vicinity of z=1+1i.) + + Note that this is an INDEPENDENT RE-IMPLEMENTATION of these algorithms, + based on the description in the papers ONLY. In particular, I did + not refer to the authors' Fortran or Matlab implementations, respectively, + (which are under restrictive ACM copyright terms and therefore unusable + in free/open-source software). + + Steven G. Johnson, Massachusetts Institute of Technology + http://math.mit.edu/~stevenj + October 2012. + + -- Note that Algorithm 916 assumes that the erfc(x) function, + or rather the scaled function erfcx(x) = exp(x*x)*erfc(x), + is supplied for REAL arguments x. I originally used an + erfcx routine derived from DERFC in SLATEC, but I have + since replaced it with a much faster routine written by + me which uses a combination of continued-fraction expansions + and a lookup table of Chebyshev polynomials. For speed, + I implemented a similar algorithm for Im[w(x)] of real x, + since this comes up frequently in the other error functions. + + A small test program is included the end, which checks + the w(z) etc. results against several known values. To compile + the test function, compile with -DTEST_FADDEEVA (that is, + #define TEST_FADDEEVA). + + If HAVE_CONFIG_H is #defined (e.g. by compiling with -DHAVE_CONFIG_H), + then we #include "config.h", which is assumed to be a GNU autoconf-style + header defining HAVE_* macros to indicate the presence of features. In + particular, if HAVE_ISNAN and HAVE_ISINF are #defined, we use those + functions in math.h instead of defining our own, and if HAVE_ERF and/or + HAVE_ERFC are defined we use those functions from for erf and + erfc of real arguments, respectively, instead of defining our own. + + REVISION HISTORY: + 4 October 2012: Initial public release (SGJ) + 5 October 2012: Revised (SGJ) to fix spelling error, + start summation for large x at round(x/a) (> 1) + rather than ceil(x/a) as in the original + paper, which should slightly improve performance + (and, apparently, slightly improves accuracy) + 19 October 2012: Revised (SGJ) to fix bugs for large x, large -y, + and 15 1e154. + Set relerr argument to min(relerr,0.1). + 27 October 2012: Enhance accuracy in Re[w(z)] taken by itself, + by switching to Alg. 916 in a region near + the real-z axis where continued fractions + have poor relative accuracy in Re[w(z)]. Thanks + to M. Zaghloul for the tip. + 29 October 2012: Replace SLATEC-derived erfcx routine with + completely rewritten code by me, using a very + different algorithm which is much faster. + 30 October 2012: Implemented special-case code for real z + (where real part is exp(-x^2) and imag part is + Dawson integral), using algorithm similar to erfx. + Export ImFaddeeva_w function to make Dawson's + integral directly accessible. + 3 November 2012: Provide implementations of erf, erfc, erfcx, + and Dawson functions in Faddeeva:: namespace, + in addition to Faddeeva::w. Provide header + file Faddeeva.hh. + 4 November 2012: Slightly faster erf for real arguments. + Updated MATLAB and Octave plugins. + 27 November 2012: Support compilation with either C++ or + plain C (using C99 complex numbers). + For real x, use standard-library erf(x) + and erfc(x) if available (for C99 or C++11). + #include "config.h" if HAVE_CONFIG_H is #defined. + 15 December 2012: Portability fixes (copysign, Inf/NaN creation), + use CMPLX/__builtin_complex if available in C, + slight accuracy improvements to erf and dawson + functions near the origin. Use gnulib functions + if GNULIB_NAMESPACE is defined. + 18 December 2012: Slight tweaks (remove recomputation of x*x in Dawson) + 12 May 2015: Bugfix for systems lacking copysign function. +*/ + +///////////////////////////////////////////////////////////////////////// +/* If this file is compiled as a part of a larger project, + support using an autoconf-style config.h header file + (with various "HAVE_*" #defines to indicate features) + if HAVE_CONFIG_H is #defined (in GNU autotools style). */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +///////////////////////////////////////////////////////////////////////// +// macros to allow us to use either C++ or C (with C99 features) + +#ifdef __cplusplus + +#include "Faddeeva.hh" + +#include +#include +#include +using namespace std; + +// use std::numeric_limits, since 1./0. and 0./0. fail with some compilers (MS) +#define Inf numeric_limits::infinity() +#define NaN numeric_limits::quiet_NaN() + +typedef complex cmplx; + +// Use C-like complex syntax, since the C syntax is more restrictive +#define cexp(z) exp(z) +#define creal(z) real(z) +#define cimag(z) imag(z) +#define cpolar(r, t) polar(r, t) + +#define C(a, b) cmplx(a, b) + +#define FADDEEVA(name) Faddeeva::name +#define FADDEEVA_RE(name) Faddeeva::name + +// isnan/isinf were introduced in C++11 +#if (__cplusplus < 201103L) && (!defined(HAVE_ISNAN) || !defined(HAVE_ISINF)) +static inline bool my_isnan(double x) { return x != x; } +#define isnan my_isnan +static inline bool my_isinf(double x) { return 1 / x == 0.; } +#define isinf my_isinf +#elif (__cplusplus >= 201103L) +// g++ gets confused between the C and C++ isnan/isinf functions +#define isnan std::isnan +#define isinf std::isinf +#endif + +// copysign was introduced in C++11 (and is also in POSIX and C99) +#if defined(_WIN32) || defined(__WIN32__) +#define copysign _copysign // of course MS had to be different +#elif defined(GNULIB_NAMESPACE) // we are using using gnulib +#define copysign GNULIB_NAMESPACE::copysign +#elif (__cplusplus < 201103L) && !defined(HAVE_COPYSIGN) && \ + !defined(__linux__) && !(defined(__APPLE__) && defined(__MACH__)) && \ + !defined(_AIX) +static inline double my_copysign(double x, double y) { + return x < 0 != y < 0 ? -x : x; +} +#define copysign my_copysign +#endif + +// If we are using the gnulib (e.g. in the GNU Octave sources), +// gnulib generates a link warning if we use ::floor instead of gnulib::floor. +// This warning is completely innocuous because the only difference between +// gnulib::floor and the system ::floor (and only on ancient OSF systems) +// has to do with floor(-0), which doesn't occur in the usage below, but +// the Octave developers prefer that we silence the warning. +#ifdef GNULIB_NAMESPACE +#define floor GNULIB_NAMESPACE::floor +#endif + +#else // !__cplusplus, i.e. pure C (requires C99 features) + +#include "Faddeeva.h" + +#define _GNU_SOURCE // enable GNU libc NAN extension if possible + +#include +#include + +typedef double complex cmplx; + +#define FADDEEVA(name) Faddeeva_##name +#define FADDEEVA_RE(name) Faddeeva_##name##_re + +/* Constructing complex numbers like 0+i*NaN is problematic in C99 + without the C11 CMPLX macro, because 0.+I*NAN may give NaN+i*NAN if + I is a complex (rather than imaginary) constant. For some reason, + however, it works fine in (pre-4.7) gcc if I define Inf and NaN as + 1/0 and 0/0 (and only if I compile with optimization -O1 or more), + but not if I use the INFINITY or NAN macros. */ + +/* __builtin_complex was introduced in gcc 4.7, but the C11 CMPLX macro + may not be defined unless we are using a recent (2012) version of + glibc and compile with -std=c11... note that icc lies about being + gcc and probably doesn't have this builtin(?), so exclude icc explicitly */ +#if !defined(CMPLX) && \ + (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) && \ + !(defined(__ICC) || defined(__INTEL_COMPILER)) +#define CMPLX(a, b) __builtin_complex((double)(a), (double)(b)) +#endif + +#ifdef CMPLX // C11 +#define C(a, b) CMPLX(a, b) +#define Inf INFINITY // C99 infinity +#ifdef NAN // GNU libc extension +#define NaN NAN +#else +#define NaN (0. / 0.) // NaN +#endif +#else +#define C(a, b) ((a) + I * (b)) +#define Inf (1. / 0.) +#define NaN (0. / 0.) +#endif + +static inline cmplx cpolar(double r, double t) { + if (r == 0.0 && !isnan(t)) + return 0.0; + else + return C(r * cos(t), r * sin(t)); +} + +#endif // !__cplusplus, i.e. pure C (requires C99 features) + +///////////////////////////////////////////////////////////////////////// +// Auxiliary routines to compute other special functions based on w(z) + +// compute erfcx(z) = exp(z^2) erfz(z) +cmplx FADDEEVA(erfcx)(cmplx z, double relerr) { + return FADDEEVA(w)(C(-cimag(z), creal(z)), relerr); +} + +// compute the error function erf(x) +double FADDEEVA_RE(erf)(double x) { +#if !defined(__cplusplus) + return erf(x); // C99 supplies erf in math.h +#elif (__cplusplus >= 201103L) || defined(HAVE_ERF) + return ::erf(x); // C++11 supplies std::erf in cmath +#else + double mx2 = -x * x; + if (mx2 < -750) // underflow + return (x >= 0 ? 1.0 : -1.0); + + if (x >= 0) { + if (x < 8e-2) + goto taylor; + return 1.0 - exp(mx2) * FADDEEVA_RE(erfcx)(x); + } else { // x < 0 + if (x > -8e-2) + goto taylor; + return exp(mx2) * FADDEEVA_RE(erfcx)(-x) - 1.0; + } + + // Use Taylor series for small |x|, to avoid cancellation inaccuracy + // erf(x) = 2/sqrt(pi) * x * (1 - x^2/3 + x^4/10 - x^6/42 + x^8/216 + ...) +taylor: + return x * (1.1283791670955125739 + + mx2 * (0.37612638903183752464 + + mx2 * (0.11283791670955125739 + + mx2 * (0.026866170645131251760 + + mx2 * 0.0052239776254421878422)))); +#endif +} + +// compute the error function erf(z) +cmplx FADDEEVA(erf)(cmplx z, double relerr) { + double x = creal(z), y = cimag(z); + + if (y == 0) + return C(FADDEEVA_RE(erf)(x), + y); // preserve sign of 0 + if (x == 0) // handle separately for speed & handling of y = Inf or NaN + return C(x, // preserve sign of 0 + /* handle y -> Inf limit manually, since + exp(y^2) -> Inf but Im[w(y)] -> 0, so + IEEE will give us a NaN when it should be Inf */ + y * y > 720 ? (y > 0 ? Inf : -Inf) + : exp(y * y) * FADDEEVA(w_im)(y)); + + double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow + double mIm_z2 = -2 * x * y; // Im(-z^2) + if (mRe_z2 < -750) // underflow + return (x >= 0 ? 1.0 : -1.0); + + /* Handle positive and negative x via different formulas, + using the mirror symmetries of w, to avoid overflow/underflow + problems from multiplying exponentially large and small quantities. */ + if (x >= 0) { + if (x < 8e-2) { + if (fabs(y) < 1e-2) + goto taylor; + else if (fabs(mIm_z2) < 5e-3 && x < 5e-3) + goto taylor_erfi; + } + /* don't use complex exp function, since that will produce spurious NaN + values when multiplying w in an overflow situation. */ + return 1.0 - exp(mRe_z2) * (C(cos(mIm_z2), sin(mIm_z2)) * + FADDEEVA(w)(C(-y, x), relerr)); + } else { // x < 0 + if (x > -8e-2) { // duplicate from above to avoid fabs(x) call + if (fabs(y) < 1e-2) + goto taylor; + else if (fabs(mIm_z2) < 5e-3 && x > -5e-3) + goto taylor_erfi; + } else if (isnan(x)) + return C(NaN, y == 0 ? 0 : NaN); + /* don't use complex exp function, since that will produce spurious NaN + values when multiplying w in an overflow situation. */ + return exp(mRe_z2) * + (C(cos(mIm_z2), sin(mIm_z2)) * FADDEEVA(w)(C(y, -x), relerr)) - + 1.0; + } + + // Use Taylor series for small |z|, to avoid cancellation inaccuracy + // erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...) +taylor: { + cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2 + return z * (1.1283791670955125739 + + mz2 * (0.37612638903183752464 + + mz2 * (0.11283791670955125739 + + mz2 * (0.026866170645131251760 + + mz2 * 0.0052239776254421878422)))); +} + + /* for small |x| and small |xy|, + use Taylor series to avoid cancellation inaccuracy: + erf(x+iy) = erf(iy) + + 2*exp(y^2)/sqrt(pi) * + [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ... + - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ] + where: + erf(iy) = exp(y^2) * Im[w(y)] + */ +taylor_erfi: { + double x2 = x * x, y2 = y * y; + double expy2 = exp(y2); + return C( + expy2 * x * + (1.1283791670955125739 - + x2 * (0.37612638903183752464 + 0.75225277806367504925 * y2) + + x2 * x2 * + (0.11283791670955125739 + + y2 * (0.45135166683820502956 + 0.15045055561273500986 * y2))), + expy2 * (FADDEEVA(w_im)(y) - x2 * y * + (1.1283791670955125739 - + x2 * (0.56418958354775628695 + + 0.37612638903183752464 * y2)))); +} +} + +// erfi(z) = -i erf(iz) +cmplx FADDEEVA(erfi)(cmplx z, double relerr) { + cmplx e = FADDEEVA(erf)(C(-cimag(z), creal(z)), relerr); + return C(cimag(e), -creal(e)); +} + +// erfi(x) = -i erf(ix) +double FADDEEVA_RE(erfi)(double x) { + return x * x > 720 ? (x > 0 ? Inf : -Inf) : exp(x * x) * FADDEEVA(w_im)(x); +} + +// erfc(x) = 1 - erf(x) +double FADDEEVA_RE(erfc)(double x) { +#if !defined(__cplusplus) + return erfc(x); // C99 supplies erfc in math.h +#elif (__cplusplus >= 201103L) || defined(HAVE_ERFC) + return ::erfc(x); // C++11 supplies std::erfc in cmath +#else + if (x * x > 750) // underflow + return (x >= 0 ? 0.0 : 2.0); + return x >= 0 ? exp(-x * x) * FADDEEVA_RE(erfcx)(x) + : 2. - exp(-x * x) * FADDEEVA_RE(erfcx)(-x); +#endif +} + +// erfc(z) = 1 - erf(z) +cmplx FADDEEVA(erfc)(cmplx z, double relerr) { + double x = creal(z), y = cimag(z); + + if (x == 0.) + return C(1, + /* handle y -> Inf limit manually, since + exp(y^2) -> Inf but Im[w(y)] -> 0, so + IEEE will give us a NaN when it should be Inf */ + y * y > 720 ? (y > 0 ? -Inf : Inf) + : -exp(y * y) * FADDEEVA(w_im)(y)); + if (y == 0.) { + if (x * x > 750) // underflow + return C(x >= 0 ? 0.0 : 2.0, + -y); // preserve sign of 0 + return C(x >= 0 ? exp(-x * x) * FADDEEVA_RE(erfcx)(x) + : 2. - exp(-x * x) * FADDEEVA_RE(erfcx)(-x), + -y); // preserve sign of zero + } + + double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow + double mIm_z2 = -2 * x * y; // Im(-z^2) + if (mRe_z2 < -750) // underflow + return (x >= 0 ? 0.0 : 2.0); + + if (x >= 0) + return cexp(C(mRe_z2, mIm_z2)) * FADDEEVA(w)(C(-y, x), relerr); + else + return 2.0 - cexp(C(mRe_z2, mIm_z2)) * FADDEEVA(w)(C(y, -x), relerr); +} + +// compute Dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x) +double FADDEEVA_RE(Dawson)(double x) { + const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2 + return spi2 * FADDEEVA(w_im)(x); +} + +// compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z) +cmplx FADDEEVA(Dawson)(cmplx z, double relerr) { + const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2 + double x = creal(z), y = cimag(z); + + // handle axes separately for speed & proper handling of x or y = Inf or NaN + if (y == 0) + return C(spi2 * FADDEEVA(w_im)(x), + -y); // preserve sign of 0 + if (x == 0) { + double y2 = y * y; + if (y2 < 2.5e-5) { // Taylor expansion + return C(x, // preserve sign of 0 + y * (1. + y2 * (0.6666666666666666666666666666666666666667 + + y2 * 0.26666666666666666666666666666666666667))); + } + return C(x, // preserve sign of 0 + spi2 * (y >= 0 ? exp(y2) - FADDEEVA_RE(erfcx)(y) + : FADDEEVA_RE(erfcx)(-y) - exp(y2))); + } + + double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow + double mIm_z2 = -2 * x * y; // Im(-z^2) + cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2 + + /* Handle positive and negative x via different formulas, + using the mirror symmetries of w, to avoid overflow/underflow + problems from multiplying exponentially large and small quantities. */ + if (y >= 0) { + if (y < 5e-3) { + if (fabs(x) < 5e-3) + goto taylor; + else if (fabs(mIm_z2) < 5e-3) + goto taylor_realaxis; + } + cmplx res = cexp(mz2) - FADDEEVA(w)(z, relerr); + return spi2 * C(-cimag(res), creal(res)); + } else { // y < 0 + if (y > -5e-3) { // duplicate from above to avoid fabs(x) call + if (fabs(x) < 5e-3) + goto taylor; + else if (fabs(mIm_z2) < 5e-3) + goto taylor_realaxis; + } else if (isnan(y)) + return C(x == 0 ? 0 : NaN, NaN); + cmplx res = FADDEEVA(w)(-z, relerr) - cexp(mz2); + return spi2 * C(-cimag(res), creal(res)); + } + + // Use Taylor series for small |z|, to avoid cancellation inaccuracy + // dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ... +taylor: + return z * (1. + mz2 * (0.6666666666666666666666666666666666666667 + + mz2 * 0.2666666666666666666666666666666666666667)); + + /* for small |y| and small |xy|, + use Taylor series to avoid cancellation inaccuracy: + dawson(x + iy) + = D + y^2 (D + x - 2Dx^2) + + y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3) + + iy [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3) + + y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ... + where D = dawson(x) + + However, for large |x|, 2Dx -> 1 which gives cancellation problems in + this series (many of the leading terms cancel). So, for large |x|, + we need to substitute a continued-fraction expansion for D. + + dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...)))))) + + The 6 terms shown here seems to be the minimum needed to be + accurate as soon as the simpler Taylor expansion above starts + breaking down. Using this 6-term expansion, factoring out the + denominator, and simplifying with Maple, we obtain: + + Re dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / x + = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4 + Im dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / y + = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4 + + Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction + expansion for the real part, and a 2-term expansion for the imaginary + part. (This avoids overflow problems for huge |x|.) This yields: + + Re dawson(x + iy) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x) + Im dawson(x + iy) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1) + + */ +taylor_realaxis: { + double x2 = x * x; + if (x2 > 1600) { // |x| > 40 + double y2 = y * y; + if (x2 > 25e14) { // |x| > 5e7 + double xy2 = (x * y) * (x * y); + return C( + (0.5 + y2 * (0.5 + 0.25 * y2 - 0.16666666666666666667 * xy2)) / x, + y * + (-1 + + y2 * (-0.66666666666666666667 + 0.13333333333333333333 * xy2 - + 0.26666666666666666667 * y2)) / + (2 * x2 - 1)); + } + return (1. / (-15 + x2 * (90 + x2 * (-60 + 8 * x2)))) * + C(x * (33 + x2 * (-28 + 4 * x2) + y2 * (18 - 4 * x2 + 4 * y2)), + y * (-15 + x2 * (24 - 4 * x2) + y2 * (4 * x2 - 10 - 4 * y2))); + } else { + double D = spi2 * FADDEEVA(w_im)(x); + double y2 = y * y; + return C( + D + y2 * (D + x - 2 * D * x2) + + y2 * y2 * + (D * (0.5 - x2 * (2 - 0.66666666666666666667 * x2)) + + x * (0.83333333333333333333 - 0.33333333333333333333 * x2)), + y * (1 - 2 * D * x + + y2 * 0.66666666666666666667 * (1 - x2 - D * x * (3 - 2 * x2)) + + y2 * y2 * + (0.26666666666666666667 - + x2 * (0.6 - 0.13333333333333333333 * x2) - + D * x * + (1 - x2 * (1.3333333333333333333 - + 0.26666666666666666667 * x2))))); + } +} +} + +///////////////////////////////////////////////////////////////////////// + +// return sinc(x) = sin(x)/x, given both x and sin(x) +// [since we only use this in cases where sin(x) has already been computed] +static inline double sinc(double x, double sinx) { + return fabs(x) < 1e-4 ? 1 - (0.1666666666666666666667) * x * x : sinx / x; +} + +// sinh(x) via Taylor series, accurate to machine precision for |x| < 1e-2 +static inline double sinh_taylor(double x) { + return x * (1 + (x * x) * (0.1666666666666666666667 + + 0.00833333333333333333333 * (x * x))); +} + +static inline double sqr(double x) { return x * x; } + +// precomputed table of expa2n2[n-1] = exp(-a2*n*n) +// for double-precision a2 = 0.26865... in FADDEEVA(w), below. +static const double expa2n2[] = { + 7.64405281671221563e-01, + 3.41424527166548425e-01, + 8.91072646929412548e-02, + 1.35887299055460086e-02, + 1.21085455253437481e-03, + 6.30452613933449404e-05, + 1.91805156577114683e-06, + 3.40969447714832381e-08, + 3.54175089099469393e-10, + 2.14965079583260682e-12, + 7.62368911833724354e-15, + 1.57982797110681093e-17, + 1.91294189103582677e-20, + 1.35344656764205340e-23, + 5.59535712428588720e-27, + 1.35164257972401769e-30, + 1.90784582843501167e-34, + 1.57351920291442930e-38, + 7.58312432328032845e-43, + 2.13536275438697082e-47, + 3.51352063787195769e-52, + 3.37800830266396920e-57, + 1.89769439468301000e-62, + 6.22929926072668851e-68, + 1.19481172006938722e-73, + 1.33908181133005953e-79, + 8.76924303483223939e-86, + 3.35555576166254986e-92, + 7.50264110688173024e-99, + 9.80192200745410268e-106, + 7.48265412822268959e-113, + 3.33770122566809425e-120, + 8.69934598159861140e-128, + 1.32486951484088852e-135, + 1.17898144201315253e-143, + 6.13039120236180012e-152, + 1.86258785950822098e-160, + 3.30668408201432783e-169, + 3.43017280887946235e-178, + 2.07915397775808219e-187, + 7.36384545323984966e-197, + 1.52394760394085741e-206, + 1.84281935046532100e-216, + 1.30209553802992923e-226, + 5.37588903521080531e-237, + 1.29689584599763145e-247, + 1.82813078022866562e-258, + 1.50576355348684241e-269, + 7.24692320799294194e-281, + 2.03797051314726829e-292, + 3.34880215927873807e-304, + 0.0 // underflow (also prevents reads past array end, below) +}; + +///////////////////////////////////////////////////////////////////////// + +cmplx FADDEEVA(w)(cmplx z, double relerr) { + if (creal(z) == 0.0) + return C(FADDEEVA_RE(erfcx)(cimag(z)), + creal(z)); // give correct sign of 0 in cimag(w) + else if (cimag(z) == 0) + return C(exp(-sqr(creal(z))), FADDEEVA(w_im)(creal(z))); + + double a, a2, c; + if (relerr <= DBL_EPSILON) { + relerr = DBL_EPSILON; + a = 0.518321480430085929872; // pi / sqrt(-log(eps*0.5)) + c = 0.329973702884629072537; // (2/pi) * a; + a2 = 0.268657157075235951582; // a^2 + } else { + const double pi = 3.14159265358979323846264338327950288419716939937510582; + if (relerr > 0.1) + relerr = 0.1; // not sensible to compute < 1 digit + a = pi / sqrt(-log(relerr * 0.5)); + c = (2 / pi) * a; + a2 = a * a; + } + const double x = fabs(creal(z)); + const double y = cimag(z), ya = fabs(y); + + cmplx ret = 0.; // return value + + double sum1 = 0, sum2 = 0, sum3 = 0, sum4 = 0, sum5 = 0; + +#define USE_CONTINUED_FRACTION 1 // 1 to use continued fraction for large |z| + +#if USE_CONTINUED_FRACTION + if (ya > 7 || (x > 6 // continued fraction is faster + /* As pointed out by M. Zaghloul, the continued + fraction seems to give a large relative error in + Re w(z) for |x| ~ 6 and small |y|, so use + algorithm 816 in this region: */ + && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28))) { + + /* Poppe & Wijers suggest using a number of terms + nu = 3 + 1442 / (26*rho + 77) + where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4. + (They only use this expansion for rho >= 1, but rho a little less + than 1 seems okay too.) + Instead, I did my own fit to a slightly different function + that avoids the hypotenuse calculation, using NLopt to minimize + the sum of the squares of the errors in nu with the constraint + that the estimated nu be >= minimum nu to attain machine precision. + I also separate the regions where nu == 2 and nu == 1. */ + const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) + double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0 + if (x + ya > 4000) { // nu <= 2 + if (x + ya > 1e7) { // nu == 1, w(z) = i/sqrt(pi) / z + // scale to avoid overflow + if (x > ya) { + double yax = ya / xs; + double denom = ispi / (xs + yax * ya); + ret = C(denom * yax, denom); + } else if (isinf(ya)) + return ((isnan(x) || y < 0) ? C(NaN, NaN) : C(0, 0)); + else { + double xya = xs / ya; + double denom = ispi / (xya * xs + ya); + ret = C(denom, denom * xya); + } + } else { // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5) + double dr = xs * xs - ya * ya - 0.5, di = 2 * xs * ya; + double denom = ispi / (dr * dr + di * di); + ret = C(denom * (xs * di - ya * dr), denom * (xs * dr + ya * di)); + } + } else { // compute nu(z) estimate and do general continued fraction + const double c0 = 3.9, c1 = 11.398, c2 = 0.08254, c3 = 0.1421, + c4 = 0.2023; // fit + double nu = floor(c0 + c1 / (c2 * x + c3 * ya + c4)); + double wr = xs, wi = ya; + for (nu = 0.5 * (nu - 1); nu > 0.4; nu -= 0.5) { + // w <- z - nu/w: + double denom = nu / (wr * wr + wi * wi); + wr = xs - wr * denom; + wi = ya + wi * denom; + } + { // w(z) = i/sqrt(pi) / w: + double denom = ispi / (wr * wr + wi * wi); + ret = C(denom * wi, denom * wr); + } + } + if (y < 0) { + // use w(z) = 2.0*exp(-z*z) - w(-z), + // but be careful of overflow in exp(-z*z) + // = exp(-(xs*xs-ya*ya) -2*i*xs*ya) + return 2.0 * cexp(C((ya - xs) * (xs + ya), 2 * xs * y)) - ret; + } else + return ret; + } +#else // !USE_CONTINUED_FRACTION + if (x + ya > 1e7) { // w(z) = i/sqrt(pi) / z, to machine precision + const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) + double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0 + // scale to avoid overflow + if (x > ya) { + double yax = ya / xs; + double denom = ispi / (xs + yax * ya); + ret = C(denom * yax, denom); + } else { + double xya = xs / ya; + double denom = ispi / (xya * xs + ya); + ret = C(denom, denom * xya); + } + if (y < 0) { + // use w(z) = 2.0*exp(-z*z) - w(-z), + // but be careful of overflow in exp(-z*z) + // = exp(-(xs*xs-ya*ya) -2*i*xs*ya) + return 2.0 * cexp(C((ya - xs) * (xs + ya), 2 * xs * y)) - ret; + } else + return ret; + } +#endif // !USE_CONTINUED_FRACTION + + /* Note: The test that seems to be suggested in the paper is x < + sqrt(-log(DBL_MIN)), about 26.6, since otherwise exp(-x^2) + underflows to zero and sum1,sum2,sum4 are zero. However, long + before this occurs, the sum1,sum2,sum4 contributions are + negligible in double precision; I find that this happens for x > + about 6, for all y. On the other hand, I find that the case + where we compute all of the sums is faster (at least with the + precomputed expa2n2 table) until about x=10. Furthermore, if we + try to compute all of the sums for x > 20, I find that we + sometimes run into numerical problems because underflow/overflow + problems start to appear in the various coefficients of the sums, + below. Therefore, we use x < 10 here. */ + else if (x < 10) { + double prod2ax = 1, prodm2ax = 1; + double expx2; + + if (isnan(y)) + return C(y, y); + + /* Somewhat ugly copy-and-paste duplication here, but I see significant + speedups from using the special-case code with the precomputed + exponential, and the x < 5e-4 special case is needed for accuracy. */ + + if (relerr == DBL_EPSILON) { // use precomputed exp(-a2*(n*n)) table + if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4 + const double x2 = x * x; + expx2 = 1 - x2 * (1 - 0.5 * x2); // exp(-x*x) via Taylor + // compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision + const double ax2 = 1.036642960860171859744 * x; // 2*a*x + const double exp2ax = + 1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667 * ax2)); + const double expm2ax = + 1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667 * ax2)); + for (int n = 1; 1; ++n) { + const double coef = expa2n2[n - 1] * expx2 / (a2 * (n * n) + y * y); + prod2ax *= exp2ax; + prodm2ax *= expm2ax; + sum1 += coef; + sum2 += coef * prodm2ax; + sum3 += coef * prod2ax; + + // really = sum5 - sum4 + sum5 += coef * (2 * a) * n * sinh_taylor((2 * a) * n * x); + + // test convergence via sum3 + if (coef * prod2ax < relerr * sum3) + break; + } + } else { // x > 5e-4, compute sum4 and sum5 separately + expx2 = exp(-x * x); + const double exp2ax = exp((2 * a) * x), expm2ax = 1 / exp2ax; + for (int n = 1; 1; ++n) { + const double coef = expa2n2[n - 1] * expx2 / (a2 * (n * n) + y * y); + prod2ax *= exp2ax; + prodm2ax *= expm2ax; + sum1 += coef; + sum2 += coef * prodm2ax; + sum4 += (coef * prodm2ax) * (a * n); + sum3 += coef * prod2ax; + sum5 += (coef * prod2ax) * (a * n); + // test convergence via sum5, since this sum has the slowest decay + if ((coef * prod2ax) * (a * n) < relerr * sum5) + break; + } + } + } else { // relerr != DBL_EPSILON, compute exp(-a2*(n*n)) on the fly + const double exp2ax = exp((2 * a) * x), expm2ax = 1 / exp2ax; + if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4 + const double x2 = x * x; + expx2 = 1 - x2 * (1 - 0.5 * x2); // exp(-x*x) via Taylor + for (int n = 1; 1; ++n) { + const double coef = + exp(-a2 * (n * n)) * expx2 / (a2 * (n * n) + y * y); + prod2ax *= exp2ax; + prodm2ax *= expm2ax; + sum1 += coef; + sum2 += coef * prodm2ax; + sum3 += coef * prod2ax; + + // really = sum5 - sum4 + sum5 += coef * (2 * a) * n * sinh_taylor((2 * a) * n * x); + + // test convergence via sum3 + if (coef * prod2ax < relerr * sum3) + break; + } + } else { // x > 5e-4, compute sum4 and sum5 separately + expx2 = exp(-x * x); + for (int n = 1; 1; ++n) { + const double coef = + exp(-a2 * (n * n)) * expx2 / (a2 * (n * n) + y * y); + prod2ax *= exp2ax; + prodm2ax *= expm2ax; + sum1 += coef; + sum2 += coef * prodm2ax; + sum4 += (coef * prodm2ax) * (a * n); + sum3 += coef * prod2ax; + sum5 += (coef * prod2ax) * (a * n); + // test convergence via sum5, since this sum has the slowest decay + if ((coef * prod2ax) * (a * n) < relerr * sum5) + break; + } + } + } + const double expx2erfcxy = // avoid spurious overflow for large negative y + y > -6 // for y < -6, erfcx(y) = 2*exp(y*y) to double precision + ? expx2 * FADDEEVA_RE(erfcx)(y) + : 2 * exp(y * y - x * x); + if (y > 5) { // imaginary terms cancel + const double sinxy = sin(x * y); + ret = (expx2erfcxy - c * y * sum1) * cos(2 * x * y) + + (c * x * expx2) * sinxy * sinc(x * y, sinxy); + } else { + double xs = creal(z); + const double sinxy = sin(xs * y); + const double sin2xy = sin(2 * xs * y), cos2xy = cos(2 * xs * y); + const double coef1 = expx2erfcxy - c * y * sum1; + const double coef2 = c * xs * expx2; + ret = C(coef1 * cos2xy + coef2 * sinxy * sinc(xs * y, sinxy), + coef2 * sinc(2 * xs * y, sin2xy) - coef1 * sin2xy); + } + } else { // x large: only sum3 & sum5 contribute (see above note) + if (isnan(x)) + return C(x, x); + if (isnan(y)) + return C(y, y); + +#if USE_CONTINUED_FRACTION + ret = exp(-x * x); // |y| < 1e-10, so we only need exp(-x*x) term +#else + if (y < 0) { + /* erfcx(y) ~ 2*exp(y*y) + (< 1) if y < 0, so + erfcx(y)*exp(-x*x) ~ 2*exp(y*y-x*x) term may not be negligible + if y*y - x*x > -36 or so. So, compute this term just in case. + We also need the -exp(-x*x) term to compute Re[w] accurately + in the case where y is very small. */ + ret = cpolar(2 * exp(y * y - x * x) - exp(-x * x), -2 * creal(z) * y); + } else + ret = exp(-x * x); // not negligible in real part if y very small +#endif + // (round instead of ceil as in original paper; note that x/a > 1 here) + double n0 = floor(x / a + 0.5); // sum in both directions, starting at n0 + double dx = a * n0 - x; + sum3 = exp(-dx * dx) / (a2 * (n0 * n0) + y * y); + sum5 = a * n0 * sum3; + double exp1 = exp(4 * a * dx), exp1dn = 1; + int dn; + for (dn = 1; n0 - dn > 0; ++dn) { // loop over n0-dn and n0+dn terms + double np = n0 + dn, nm = n0 - dn; + double tp = exp(-sqr(a * dn + dx)); + double tm = tp * (exp1dn *= exp1); // trick to get tm from tp + tp /= (a2 * (np * np) + y * y); + tm /= (a2 * (nm * nm) + y * y); + sum3 += tp + tm; + sum5 += a * (np * tp + nm * tm); + if (a * (np * tp + nm * tm) < relerr * sum5) + goto finish; + } + while (1) { // loop over n0+dn terms only (since n0-dn <= 0) + double np = n0 + dn++; + double tp = exp(-sqr(a * dn + dx)) / (a2 * (np * np) + y * y); + sum3 += tp; + sum5 += a * np * tp; + if (a * np * tp < relerr * sum5) + goto finish; + } + } +finish: + return ret + C((0.5 * c) * y * (sum2 + sum3), + (0.5 * c) * copysign(sum5 - sum4, creal(z))); +} + +///////////////////////////////////////////////////////////////////////// + +/* erfcx(x) = exp(x^2) erfc(x) function, for real x, written by + Steven G. Johnson, October 2012. + + This function combines a few different ideas. + + First, for x > 50, it uses a continued-fraction expansion (same as + for the Faddeeva function, but with algebraic simplifications for z=i*x). + + Second, for 0 <= x <= 50, it uses Chebyshev polynomial approximations, + but with two twists: + + a) It maps x to y = 4 / (4+x) in [0,1]. This simple transformation, + inspired by a similar transformation in the octave-forge/specfun + erfcx by Soren Hauberg, results in much faster Chebyshev convergence + than other simple transformations I have examined. + + b) Instead of using a single Chebyshev polynomial for the entire + [0,1] y interval, we break the interval up into 100 equal + subintervals, with a switch/lookup table, and use much lower + degree Chebyshev polynomials in each subinterval. This greatly + improves performance in my tests. + + For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x), + with the usual checks for overflow etcetera. + + Performance-wise, it seems to be substantially faster than either + the SLATEC DERFC function [or an erfcx function derived therefrom] + or Cody's CALERF function (from netlib.org/specfun), while + retaining near machine precision in accuracy. */ + +/* Given y100=100*y, where y = 4/(4+x) for x >= 0, compute erfc(x). + + Uses a look-up table of 100 different Chebyshev polynomials + for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated + with the help of Maple and a little shell script. This allows + the Chebyshev polynomials to be of significantly lower degree (about 1/4) + compared to fitting the whole [0,1] interval with a single polynomial. */ +static double erfcx_y100(double y100) { + switch ((int)y100) { + case 0: { + double t = 2 * y100 - 1; + return 0.70878032454106438663e-3 + + (0.71234091047026302958e-3 + + (0.35779077297597742384e-5 + + (0.17403143962587937815e-7 + + (0.81710660047307788845e-10 + + (0.36885022360434957634e-12 + 0.15917038551111111111e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 1: { + double t = 2 * y100 - 3; + return 0.21479143208285144230e-2 + + (0.72686402367379996033e-3 + + (0.36843175430938995552e-5 + + (0.18071841272149201685e-7 + + (0.85496449296040325555e-10 + + (0.38852037518534291510e-12 + 0.16868473576888888889e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 2: { + double t = 2 * y100 - 5; + return 0.36165255935630175090e-2 + + (0.74182092323555510862e-3 + + (0.37948319957528242260e-5 + + (0.18771627021793087350e-7 + + (0.89484715122415089123e-10 + + (0.40935858517772440862e-12 + 0.17872061464888888889e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 3: { + double t = 2 * y100 - 7; + return 0.51154983860031979264e-2 + + (0.75722840734791660540e-3 + + (0.39096425726735703941e-5 + + (0.19504168704300468210e-7 + + (0.93687503063178993915e-10 + + (0.43143925959079664747e-12 + 0.18939926435555555556e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 4: { + double t = 2 * y100 - 9; + return 0.66457513172673049824e-2 + + (0.77310406054447454920e-3 + + (0.40289510589399439385e-5 + + (0.20271233238288381092e-7 + + (0.98117631321709100264e-10 + + (0.45484207406017752971e-12 + 0.20076352213333333333e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 5: { + double t = 2 * y100 - 11; + return 0.82082389970241207883e-2 + + (0.78946629611881710721e-3 + + (0.41529701552622656574e-5 + + (0.21074693344544655714e-7 + + (0.10278874108587317989e-9 + + (0.47965201390613339638e-12 + 0.21285907413333333333e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 6: { + double t = 2 * y100 - 13; + return 0.98039537275352193165e-2 + + (0.80633440108342840956e-3 + + (0.42819241329736982942e-5 + + (0.21916534346907168612e-7 + + (0.10771535136565470914e-9 + + (0.50595972623692822410e-12 + 0.22573462684444444444e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 7: { + double t = 2 * y100 - 15; + return 0.11433927298290302370e-1 + + (0.82372858383196561209e-3 + + (0.44160495311765438816e-5 + + (0.22798861426211986056e-7 + + (0.11291291745879239736e-9 + + (0.53386189365816880454e-12 + 0.23944209546666666667e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 8: { + double t = 2 * y100 - 17; + return 0.13099232878814653979e-1 + + (0.84167002467906968214e-3 + + (0.45555958988457506002e-5 + + (0.23723907357214175198e-7 + + (0.11839789326602695603e-9 + + (0.56346163067550237877e-12 + 0.25403679644444444444e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 9: { + double t = 2 * y100 - 19; + return 0.14800987015587535621e-1 + + (0.86018092946345943214e-3 + + (0.47008265848816866105e-5 + + (0.24694040760197315333e-7 + + (0.12418779768752299093e-9 + + (0.59486890370320261949e-12 + 0.26957764568888888889e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 10: { + double t = 2 * y100 - 21; + return 0.16540351739394069380e-1 + + (0.87928458641241463952e-3 + + (0.48520195793001753903e-5 + + (0.25711774900881709176e-7 + + (0.13030128534230822419e-9 + + (0.62820097586874779402e-12 + 0.28612737351111111111e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 11: { + double t = 2 * y100 - 23; + return 0.18318536789842392647e-1 + + (0.89900542647891721692e-3 + + (0.50094684089553365810e-5 + + (0.26779777074218070482e-7 + + (0.13675822186304615566e-9 + + (0.66358287745352705725e-12 + 0.30375273884444444444e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 12: { + double t = 2 * y100 - 25; + return 0.20136801964214276775e-1 + + (0.91936908737673676012e-3 + + (0.51734830914104276820e-5 + + (0.27900878609710432673e-7 + + (0.14357976402809042257e-9 + + (0.70114790311043728387e-12 + 0.32252476000000000000e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 13: { + double t = 2 * y100 - 27; + return 0.21996459598282740954e-1 + + (0.94040248155366777784e-3 + + (0.53443911508041164739e-5 + + (0.29078085538049374673e-7 + + (0.15078844500329731137e-9 + + (0.74103813647499204269e-12 + 0.34251892320000000000e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 14: { + double t = 2 * y100 - 29; + return 0.23898877187226319502e-1 + + (0.96213386835900177540e-3 + + (0.55225386998049012752e-5 + + (0.30314589961047687059e-7 + + (0.15840826497296335264e-9 + + (0.78340500472414454395e-12 + 0.36381553564444444445e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 15: { + double t = 2 * y100 - 31; + return 0.25845480155298518485e-1 + + (0.98459293067820123389e-3 + + (0.57082915920051843672e-5 + + (0.31613782169164830118e-7 + + (0.16646478745529630813e-9 + + (0.82840985928785407942e-12 + 0.38649975768888888890e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 16: { + double t = 2 * y100 - 33; + return 0.27837754783474696598e-1 + + (0.10078108563256892757e-2 + + (0.59020366493792212221e-5 + + (0.32979263553246520417e-7 + + (0.17498524159268458073e-9 + + (0.87622459124842525110e-12 + 0.41066206488888888890e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 17: { + double t = 2 * y100 - 35; + return 0.29877251304899307550e-1 + + (0.10318204245057349310e-2 + + (0.61041829697162055093e-5 + + (0.34414860359542720579e-7 + + (0.18399863072934089607e-9 + + (0.92703227366365046533e-12 + 0.43639844053333333334e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 18: { + double t = 2 * y100 - 37; + return 0.31965587178596443475e-1 + + (0.10566560976716574401e-2 + + (0.63151633192414586770e-5 + + (0.35924638339521924242e-7 + + (0.19353584758781174038e-9 + + (0.98102783859889264382e-12 + 0.46381060817777777779e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 19: { + double t = 2 * y100 - 39; + return 0.34104450552588334840e-1 + + (0.10823541191350532574e-2 + + (0.65354356159553934436e-5 + + (0.37512918348533521149e-7 + + (0.20362979635817883229e-9 + + (0.10384187833037282363e-11 + 0.49300625262222222221e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 20: { + double t = 2 * y100 - 41; + return 0.36295603928292425716e-1 + + (0.11089526167995268200e-2 + + (0.67654845095518363577e-5 + + (0.39184292949913591646e-7 + + (0.21431552202133775150e-9 + + (0.10994259106646731797e-11 + 0.52409949102222222221e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 21: { + double t = 2 * y100 - 43; + return 0.38540888038840509795e-1 + + (0.11364917134175420009e-2 + + (0.70058230641246312003e-5 + + (0.40943644083718586939e-7 + + (0.22563034723692881631e-9 + + (0.11642841011361992885e-11 + 0.55721092871111111110e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 22: { + double t = 2 * y100 - 45; + return 0.40842225954785960651e-1 + + (0.11650136437945673891e-2 + + (0.72569945502343006619e-5 + + (0.42796161861855042273e-7 + + (0.23761401711005024162e-9 + + (0.12332431172381557035e-11 + 0.59246802364444444445e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 23: { + double t = 2 * y100 - 47; + return 0.43201627431540222422e-1 + + (0.11945628793917272199e-2 + + (0.75195743532849206263e-5 + + (0.44747364553960993492e-7 + + (0.25030885216472953674e-9 + + (0.13065684400300476484e-11 + 0.63000532853333333334e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 24: { + double t = 2 * y100 - 49; + return 0.45621193513810471438e-1 + + (0.12251862608067529503e-2 + + (0.77941720055551920319e-5 + + (0.46803119830954460212e-7 + + (0.26375990983978426273e-9 + + (0.13845421370977119765e-11 + 0.66996477404444444445e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 25: { + double t = 2 * y100 - 51; + return 0.48103121413299865517e-1 + + (0.12569331386432195113e-2 + + (0.80814333496367673980e-5 + + (0.48969667335682018324e-7 + + (0.27801515481905748484e-9 + + (0.14674637611609884208e-11 + 0.71249589351111111110e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 26: { + double t = 2 * y100 - 53; + return 0.50649709676983338501e-1 + + (0.12898555233099055810e-2 + + (0.83820428414568799654e-5 + + (0.51253642652551838659e-7 + + (0.29312563849675507232e-9 + + (0.15556512782814827846e-11 + 0.75775607822222222221e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 27: { + double t = 2 * y100 - 55; + return 0.53263363664388864181e-1 + + (0.13240082443256975769e-2 + + (0.86967260015007658418e-5 + + (0.53662102750396795566e-7 + + (0.30914568786634796807e-9 + + (0.16494420240828493176e-11 + 0.80591079644444444445e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 28: { + double t = 2 * y100 - 57; + return 0.55946601353500013794e-1 + + (0.13594491197408190706e-2 + + (0.90262520233016380987e-5 + + (0.56202552975056695376e-7 + + (0.32613310410503135996e-9 + + (0.17491936862246367398e-11 + 0.85713381688888888890e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 29: { + double t = 2 * y100 - 59; + return 0.58702059496154081813e-1 + + (0.13962391363223647892e-2 + + (0.93714365487312784270e-5 + + (0.58882975670265286526e-7 + + (0.34414937110591753387e-9 + + (0.18552853109751857859e-11 + 0.91160736711111111110e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 30: { + double t = 2 * y100 - 61; + return 0.61532500145144778048e-1 + + (0.14344426411912015247e-2 + + (0.97331446201016809696e-5 + + (0.61711860507347175097e-7 + + (0.36325987418295300221e-9 + + (0.19681183310134518232e-11 + 0.96952238400000000000e-14 * t) * + t) * + t) * + t) * + t) * + t; + } + case 31: { + double t = 2 * y100 - 63; + return 0.64440817576653297993e-1 + + (0.14741275456383131151e-2 + + (0.10112293819576437838e-4 + + (0.64698236605933246196e-7 + + (0.38353412915303665586e-9 + + (0.20881176114385120186e-11 + 0.10310784480000000000e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 32: { + double t = 2 * y100 - 65; + return 0.67430045633130393282e-1 + + (0.15153655418916540370e-2 + + (0.10509857606888328667e-4 + + (0.67851706529363332855e-7 + + (0.40504602194811140006e-9 + + (0.22157325110542534469e-11 + 0.10964842115555555556e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 33: { + double t = 2 * y100 - 67; + return 0.70503365513338850709e-1 + + (0.15582323336495709827e-2 + + (0.10926868866865231089e-4 + + (0.71182482239613507542e-7 + + (0.42787405890153386710e-9 + + (0.23514379522274416437e-11 + 0.11659571751111111111e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 34: { + double t = 2 * y100 - 69; + return 0.73664114037944596353e-1 + + (0.16028078812438820413e-2 + + (0.11364423678778207991e-4 + + (0.74701423097423182009e-7 + + (0.45210162777476488324e-9 + + (0.24957355004088569134e-11 + 0.12397238257777777778e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 35: { + double t = 2 * y100 - 71; + return 0.76915792420819562379e-1 + + (0.16491766623447889354e-2 + + (0.11823685320041302169e-4 + + (0.78420075993781544386e-7 + + (0.47781726956916478925e-9 + + (0.26491544403815724749e-11 + 0.13180196462222222222e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 36: { + double t = 2 * y100 - 73; + return 0.80262075578094612819e-1 + + (0.16974279491709504117e-2 + + (0.12305888517309891674e-4 + + (0.82350717698979042290e-7 + + (0.50511496109857113929e-9 + + (0.28122528497626897696e-11 + 0.14010889635555555556e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 37: { + double t = 2 * y100 - 75; + return 0.83706822008980357446e-1 + + (0.17476561032212656962e-2 + + (0.12812343958540763368e-4 + + (0.86506399515036435592e-7 + + (0.53409440823869467453e-9 + + (0.29856186620887555043e-11 + 0.14891851591111111111e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 38: { + double t = 2 * y100 - 77; + return 0.87254084284461718231e-1 + + (0.17999608886001962327e-2 + + (0.13344443080089492218e-4 + + (0.90900994316429008631e-7 + + (0.56486134972616465316e-9 + + (0.31698707080033956934e-11 + 0.15825697795555555556e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 39: { + double t = 2 * y100 - 79; + return 0.90908120182172748487e-1 + + (0.18544478050657699758e-2 + + (0.13903663143426120077e-4 + + (0.95549246062549906177e-7 + + (0.59752787125242054315e-9 + + (0.33656597366099099413e-11 + 0.16815130613333333333e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 40: { + double t = 2 * y100 - 81; + return 0.94673404508075481121e-1 + + (0.19112284419887303347e-2 + + (0.14491572616545004930e-4 + + (0.10046682186333613697e-6 + + (0.63221272959791000515e-9 + + (0.35736693975589130818e-11 + 0.17862931591111111111e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 41: { + double t = 2 * y100 - 83; + return 0.98554641648004456555e-1 + + (0.19704208544725622126e-2 + + (0.15109836875625443935e-4 + + (0.10567036667675984067e-6 + + (0.66904168640019354565e-9 + + (0.37946171850824333014e-11 + 0.18971959040000000000e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 42: { + double t = 2 * y100 - 85; + return 0.10255677889470089531e0 + + (0.20321499629472857418e-2 + + (0.15760224242962179564e-4 + + (0.11117756071353507391e-6 + + (0.70814785110097658502e-9 + + (0.40292553276632563925e-11 + 0.20145143075555555556e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 43: { + double t = 2 * y100 - 87; + return 0.10668502059865093318e0 + + (0.20965479776148731610e-2 + + (0.16444612377624983565e-4 + + (0.11700717962026152749e-6 + + (0.74967203250938418991e-9 + + (0.42783716186085922176e-11 + 0.21385479360000000000e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 44: { + double t = 2 * y100 - 89; + return 0.11094484319386444474e0 + + (0.21637548491908170841e-2 + + (0.17164995035719657111e-4 + + (0.12317915750735938089e-6 + + (0.79376309831499633734e-9 + + (0.45427901763106353914e-11 + 0.22696025653333333333e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 45: { + double t = 2 * y100 - 91; + return 0.11534201115268804714e0 + + (0.22339187474546420375e-2 + + (0.17923489217504226813e-4 + + (0.12971465288245997681e-6 + + (0.84057834180389073587e-9 + + (0.48233721206418027227e-11 + 0.24079890062222222222e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 46: { + double t = 2 * y100 - 93; + return 0.11988259392684094740e0 + + (0.23071965691918689601e-2 + + (0.18722342718958935446e-4 + + (0.13663611754337957520e-6 + + (0.89028385488493287005e-9 + + (0.51210161569225846701e-11 + 0.25540227111111111111e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 47: { + double t = 2 * y100 - 95; + return 0.12457298393509812907e0 + + (0.23837544771809575380e-2 + + (0.19563942105711612475e-4 + + (0.14396736847739470782e-6 + + (0.94305490646459247016e-9 + + (0.54366590583134218096e-11 + 0.27080225920000000000e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 48: { + double t = 2 * y100 - 97; + return 0.12941991566142438816e0 + + (0.24637684719508859484e-2 + + (0.20450821127475879816e-4 + + (0.15173366280523906622e-6 + + (0.99907632506389027739e-9 + + (0.57712760311351625221e-11 + 0.28703099555555555556e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 49: { + double t = 2 * y100 - 99; + return 0.13443048593088696613e0 + + (0.25474249981080823877e-2 + + (0.21385669591362915223e-4 + + (0.15996177579900443030e-6 + + (0.10585428844575134013e-8 + + (0.61258809536787882989e-11 + 0.30412080142222222222e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 50: { + double t = 2 * y100 - 101; + return 0.13961217543434561353e0 + + (0.26349215871051761416e-2 + + (0.22371342712572567744e-4 + + (0.16868008199296822247e-6 + + (0.11216596910444996246e-8 + + (0.65015264753090890662e-11 + 0.32210394506666666666e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 51: { + double t = 2 * y100 - 103; + return 0.14497287157673800690e0 + + (0.27264675383982439814e-2 + + (0.23410870961050950197e-4 + + (0.17791863939526376477e-6 + + (0.11886425714330958106e-8 + + (0.68993039665054288034e-11 + 0.34101266222222222221e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 52: { + double t = 2 * y100 - 105; + return 0.15052089272774618151e0 + + (0.28222846410136238008e-2 + + (0.24507470422713397006e-4 + + (0.18770927679626136909e-6 + + (0.12597184587583370712e-8 + + (0.73203433049229821618e-11 + 0.36087889048888888890e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 53: { + double t = 2 * y100 - 107; + return 0.15626501395774612325e0 + + (0.29226079376196624949e-2 + + (0.25664553693768450545e-4 + + (0.19808568415654461964e-6 + + (0.13351257759815557897e-8 + + (0.77658124891046760667e-11 + 0.38173420035555555555e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 54: { + double t = 2 * y100 - 109; + return 0.16221449434620737567e0 + + (0.30276865332726475672e-2 + + (0.26885741326534564336e-4 + + (0.20908350604346384143e-6 + + (0.14151148144240728728e-8 + + (0.82369170665974313027e-11 + 0.40360957457777777779e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 55: { + double t = 2 * y100 - 111; + return 0.16837910595412130659e0 + + (0.31377844510793082301e-2 + + (0.28174873844911175026e-4 + + (0.22074043807045782387e-6 + + (0.14999481055996090039e-8 + + (0.87348993661930809254e-11 + 0.42653528977777777779e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 56: { + double t = 2 * y100 - 113; + return 0.17476916455659369953e0 + + (0.32531815370903068316e-2 + + (0.29536024347344364074e-4 + + (0.23309632627767074202e-6 + + (0.15899007843582444846e-8 + + (0.92610375235427359475e-11 + 0.45054073102222222221e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 57: { + double t = 2 * y100 - 115; + return 0.18139556223643701364e0 + + (0.33741744168096996041e-2 + + (0.30973511714709500836e-4 + + (0.24619326937592290996e-6 + + (0.16852609412267750744e-8 + + (0.98166442942854895573e-11 + 0.47565418097777777779e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 58: { + double t = 2 * y100 - 117; + return 0.18826980194443664549e0 + + (0.35010775057740317997e-2 + + (0.32491914440014267480e-4 + + (0.26007572375886319028e-6 + + (0.17863299617388376116e-8 + + (0.10403065638343878679e-10 + 0.50190265831111111110e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 59: { + double t = 2 * y100 - 119; + return 0.19540403413693967350e0 + + (0.36342240767211326315e-2 + + (0.34096085096200907289e-4 + + (0.27479061117017637474e-6 + + (0.18934228504790032826e-8 + + (0.11021679075323598664e-10 + 0.52931171733333333334e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 60: { + double t = 2 * y100 - 121; + return 0.20281109560651886959e0 + + (0.37739673859323597060e-2 + + (0.35791165457592409054e-4 + + (0.29038742889416172404e-6 + + (0.20068685374849001770e-8 + + (0.11673891799578381999e-10 + 0.55790523093333333334e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 61: { + double t = 2 * y100 - 123; + return 0.21050455062669334978e0 + + (0.39206818613925652425e-2 + + (0.37582602289680101704e-4 + + (0.30691836231886877385e-6 + + (0.21270101645763677824e-8 + + (0.12361138551062899455e-10 + 0.58770520160000000000e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 62: { + double t = 2 * y100 - 125; + return 0.21849873453703332479e0 + + (0.40747643554689586041e-2 + + (0.39476163820986711501e-4 + + (0.32443839970139918836e-6 + + (0.22542053491518680200e-8 + + (0.13084879235290858490e-10 + 0.61873153262222222221e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 63: { + double t = 2 * y100 - 127; + return 0.22680879990043229327e0 + + (0.42366354648628516935e-2 + + (0.41477956909656896779e-4 + + (0.34300544894502810002e-6 + + (0.23888264229264067658e-8 + + (0.13846596292818514601e-10 + 0.65100183751111111110e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 64: { + double t = 2 * y100 - 129; + return 0.23545076536988703937e0 + + (0.44067409206365170888e-2 + + (0.43594444916224700881e-4 + + (0.36268045617760415178e-6 + + (0.25312606430853202748e-8 + + (0.14647791812837903061e-10 + 0.68453122631111111110e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 65: { + double t = 2 * y100 - 131; + return 0.24444156740777432838e0 + + (0.45855530511605787178e-2 + + (0.45832466292683085475e-4 + + (0.38352752590033030472e-6 + + (0.26819103733055603460e-8 + + (0.15489984390884756993e-10 + 0.71933206364444444445e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 66: { + double t = 2 * y100 - 133; + return 0.25379911500634264643e0 + + (0.47735723208650032167e-2 + + (0.48199253896534185372e-4 + + (0.40561404245564732314e-6 + + (0.28411932320871165585e-8 + + (0.16374705736458320149e-10 + 0.75541379822222222221e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 67: { + double t = 2 * y100 - 135; + return 0.26354234756393613032e0 + + (0.49713289477083781266e-2 + + (0.50702455036930367504e-4 + + (0.42901079254268185722e-6 + + (0.30095422058900481753e-8 + + (0.17303497025347342498e-10 + 0.79278273368888888890e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 68: { + double t = 2 * y100 - 137; + return 0.27369129607732343398e0 + + (0.51793846023052643767e-2 + + (0.53350152258326602629e-4 + + (0.45379208848865015485e-6 + + (0.31874057245814381257e-8 + + (0.18277905010245111046e-10 + 0.83144182364444444445e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 69: { + double t = 2 * y100 - 139; + return 0.28426714781640316172e0 + + (0.53983341916695141966e-2 + + (0.56150884865255810638e-4 + + (0.48003589196494734238e-6 + + (0.33752476967570796349e-8 + + (0.19299477888083469086e-10 + 0.87139049137777777779e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 70: { + double t = 2 * y100 - 141; + return 0.29529231465348519920e0 + + (0.56288077305420795663e-2 + + (0.59113671189913307427e-4 + + (0.50782393781744840482e-6 + + (0.35735475025851713168e-8 + + (0.20369760937017070382e-10 + 0.91262442613333333334e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 71: { + double t = 2 * y100 - 143; + return 0.30679050522528838613e0 + + (0.58714723032745403331e-2 + + (0.62248031602197686791e-4 + + (0.53724185766200945789e-6 + + (0.37827999418960232678e-8 + + (0.21490291930444538307e-10 + 0.95513539182222222221e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 72: { + double t = 2 * y100 - 145; + return 0.31878680111173319425e0 + + (0.61270341192339103514e-2 + + (0.65564012259707640976e-4 + + (0.56837930287837738996e-6 + + (0.40035151353392378882e-8 + + (0.22662596341239294792e-10 + 0.99891109760000000000e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 73: { + double t = 2 * y100 - 147; + return 0.33130773722152622027e0 + + (0.63962406646798080903e-2 + + (0.69072209592942396666e-4 + + (0.60133006661885941812e-6 + + (0.42362183765883466691e-8 + + (0.23888182347073698382e-10 + 0.10439349811555555556e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 74: { + double t = 2 * y100 - 149; + return 0.34438138658041336523e0 + + (0.66798829540414007258e-2 + + (0.72783795518603561144e-4 + + (0.63619220443228800680e-6 + + (0.44814499336514453364e-8 + + (0.25168535651285475274e-10 + 0.10901861383111111111e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 75: { + double t = 2 * y100 - 151; + return 0.35803744972380175583e0 + + (0.69787978834882685031e-2 + + (0.76710543371454822497e-4 + + (0.67306815308917386747e-6 + + (0.47397647975845228205e-8 + + (0.26505114141143050509e-10 + 0.11376390933333333333e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 76: { + double t = 2 * y100 - 153; + return 0.37230734890119724188e0 + + (0.72938706896461381003e-2 + + (0.80864854542670714092e-4 + + (0.71206484718062688779e-6 + + (0.50117323769745883805e-8 + + (0.27899342394100074165e-10 + 0.11862637614222222222e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 77: { + double t = 2 * y100 - 155; + return 0.38722432730555448223e0 + + (0.76260375162549802745e-2 + + (0.85259785810004603848e-4 + + (0.75329383305171327677e-6 + + (0.52979361368388119355e-8 + + (0.29352606054164086709e-10 + 0.12360253370666666667e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 78: { + double t = 2 * y100 - 157; + return 0.40282355354616940667e0 + + (0.79762880915029728079e-2 + + (0.89909077342438246452e-4 + + (0.79687137961956194579e-6 + + (0.55989731807360403195e-8 + + (0.30866246101464869050e-10 + 0.12868841946666666667e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 79: { + double t = 2 * y100 - 159; + return 0.41914223158913787649e0 + + (0.83456685186950463538e-2 + + (0.94827181359250161335e-4 + + (0.84291858561783141014e-6 + + (0.59154537751083485684e-8 + + (0.32441553034347469291e-10 + 0.13387957943111111111e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 80: { + double t = 2 * y100 - 161; + return 0.43621971639463786896e0 + + (0.87352841828289495773e-2 + + (0.10002929142066799966e-3 + + (0.89156148280219880024e-6 + + (0.62480008150788597147e-8 + + (0.34079760983458878910e-10 + 0.13917107176888888889e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 81: { + double t = 2 * y100 - 163; + return 0.45409763548534330981e0 + + (0.91463027755548240654e-2 + + (0.10553137232446167258e-3 + + (0.94293113464638623798e-6 + + (0.65972492312219959885e-8 + + (0.35782041795476563662e-10 + 0.14455745872000000000e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 82: { + double t = 2 * y100 - 165; + return 0.47282001668512331468e0 + + (0.95799574408860463394e-2 + + (0.11135019058000067469e-3 + + (0.99716373005509038080e-6 + + (0.69638453369956970347e-8 + + (0.37549499088161345850e-10 + 0.15003280712888888889e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 83: { + double t = 2 * y100 - 167; + return 0.49243342227179841649e0 + + (0.10037550043909497071e-1 + + (0.11750334542845234952e-3 + + (0.10544006716188967172e-5 + + (0.73484461168242224872e-8 + + (0.39383162326435752965e-10 + 0.15559069118222222222e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 84: { + double t = 2 * y100 - 169; + return 0.51298708979209258326e0 + + (0.10520454564612427224e-1 + + (0.12400930037494996655e-3 + + (0.11147886579371265246e-5 + + (0.77517184550568711454e-8 + + (0.41283980931872622611e-10 + 0.16122419680000000000e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 85: { + double t = 2 * y100 - 171; + return 0.53453307979101369843e0 + + (0.11030120618800726938e-1 + + (0.13088741519572269581e-3 + + (0.11784797595374515432e-5 + + (0.81743383063044825400e-8 + + (0.43252818449517081051e-10 + 0.16692592640000000000e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 86: { + double t = 2 * y100 - 173; + return 0.55712643071169299478e0 + + (0.11568077107929735233e-1 + + (0.13815797838036651289e-3 + + (0.12456314879260904558e-5 + + (0.86169898078969313597e-8 + + (0.45290446811539652525e-10 + 0.17268801084444444444e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 87: { + double t = 2 * y100 - 175; + return 0.58082532122519320968e0 + + (0.12135935999503877077e-1 + + (0.14584223996665838559e-3 + + (0.13164068573095710742e-5 + + (0.90803643355106020163e-8 + + (0.47397540713124619155e-10 + 0.17850211608888888889e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 88: { + double t = 2 * y100 - 177; + return 0.60569124025293375554e0 + + (0.12735396239525550361e-1 + + (0.15396244472258863344e-3 + + (0.13909744385382818253e-5 + + (0.95651595032306228245e-8 + + (0.49574672127669041550e-10 + 0.18435945564444444444e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 89: { + double t = 2 * y100 - 179; + return 0.63178916494715716894e0 + + (0.13368247798287030927e-1 + + (0.16254186562762076141e-3 + + (0.14695084048334056083e-5 + + (0.10072078109604152350e-7 + + (0.51822304995680707483e-10 + 0.19025081422222222222e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 90: { + double t = 2 * y100 - 181; + return 0.65918774689725319200e0 + + (0.14036375850601992063e-1 + + (0.17160483760259706354e-3 + + (0.15521885688723188371e-5 + + (0.10601827031535280590e-7 + + (0.54140790105837520499e-10 + 0.19616655146666666667e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 91: { + double t = 2 * y100 - 183; + return 0.68795950683174433822e0 + + (0.14741765091365869084e-1 + + (0.18117679143520433835e-3 + + (0.16392004108230585213e-5 + + (0.11155116068018043001e-7 + + (0.56530360194925690374e-10 + 0.20209663662222222222e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 92: { + double t = 2 * y100 - 185; + return 0.71818103808729967036e0 + + (0.15486504187117112279e-1 + + (0.19128428784550923217e-3 + + (0.17307350969359975848e-5 + + (0.11732656736113607751e-7 + + (0.58991125287563833603e-10 + 0.20803065333333333333e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 93: { + double t = 2 * y100 - 187; + return 0.74993321911726254661e0 + + (0.16272790364044783382e-1 + + (0.20195505163377912645e-3 + + (0.18269894883203346953e-5 + + (0.12335161021630225535e-7 + + (0.61523068312169087227e-10 + 0.21395783431111111111e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 94: { + double t = 2 * y100 - 189; + return 0.78330143531283492729e0 + + (0.17102934132652429240e-1 + + (0.21321800585063327041e-3 + + (0.19281661395543913713e-5 + + (0.12963340087354341574e-7 + + (0.64126040998066348872e-10 + 0.21986708942222222222e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 95: { + double t = 2 * y100 - 191; + return 0.81837581041023811832e0 + + (0.17979364149044223802e-1 + + (0.22510330592753129006e-3 + + (0.20344732868018175389e-5 + + (0.13617902941839949718e-7 + + (0.66799760083972474642e-10 + 0.22574701262222222222e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 96: { + double t = 2 * y100 - 193; + return 0.85525144775685126237e0 + + (0.18904632212547561026e-1 + + (0.23764237370371255638e-3 + + (0.21461248251306387979e-5 + + (0.14299555071870523786e-7 + + (0.69543803864694171934e-10 + 0.23158593688888888889e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 97: { + double t = 2 * y100 - 195; + return 0.89402868170849933734e0 + + (0.19881418399127202569e-1 + + (0.25086793128395995798e-3 + + (0.22633402747585233180e-5 + + (0.15008997042116532283e-7 + + (0.72357609075043941261e-10 + 0.23737194737777777778e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 98: { + double t = 2 * y100 - 197; + return 0.93481333942870796363e0 + + (0.20912536329780368893e-1 + + (0.26481403465998477969e-3 + + (0.23863447359754921676e-5 + + (0.15746923065472184451e-7 + + (0.75240468141720143653e-10 + 0.24309291271111111111e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 99: { + double t = 2 * y100 - 199; + return 0.97771701335885035464e0 + + (0.22000938572830479551e-1 + + (0.27951610702682383001e-3 + + (0.25153688325245314530e-5 + + (0.16514019547822821453e-7 + + (0.78191526829368231251e-10 + 0.24873652355555555556e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + } + // we only get here if y = 1, i.e. |x| < 4*eps, in which case + // erfcx is within 1e-15 of 1.. + return 1.0; +} + +double FADDEEVA_RE(erfcx)(double x) { + if (x >= 0) { + if (x > 50) { // continued-fraction expansion is faster + const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) + if (x > 5e7) // 1-term expansion, important to avoid overflow + return ispi / x; + /* 5-term expansion (rely on compiler for CSE), simplified from: + ispi / (x+0.5/(x+1/(x+1.5/(x+2/x)))) */ + return ispi * ((x * x) * (x * x + 4.5) + 2) / + (x * ((x * x) * (x * x + 5) + 3.75)); + } + return erfcx_y100(400 / (4 + x)); + } else + return x < -26.7 ? HUGE_VAL + : (x < -6.1 ? 2 * exp(x * x) + : 2 * exp(x * x) - erfcx_y100(400 / (4 - x))); +} + +///////////////////////////////////////////////////////////////////////// +/* Compute a scaled Dawson integral + FADDEEVA(w_im)(x) = 2*Dawson(x)/sqrt(pi) + equivalent to the imaginary part w(x) for real x. + + Uses methods similar to the erfcx calculation above: continued fractions + for large |x|, a lookup table of Chebyshev polynomials for smaller |x|, + and finally a Taylor expansion for |x|<0.01. + + Steven G. Johnson, October 2012. */ + +/* Given y100=100*y, where y = 1/(1+x) for x >= 0, compute w_im(x). + + Uses a look-up table of 100 different Chebyshev polynomials + for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated + with the help of Maple and a little shell script. This allows + the Chebyshev polynomials to be of significantly lower degree (about 1/30) + compared to fitting the whole [0,1] interval with a single polynomial. */ +static double w_im_y100(double y100, double x) { + switch ((int)y100) { + case 0: { + double t = 2 * y100 - 1; + return 0.28351593328822191546e-2 + + (0.28494783221378400759e-2 + + (0.14427470563276734183e-4 + + (0.10939723080231588129e-6 + + (0.92474307943275042045e-9 + + (0.89128907666450075245e-11 + 0.92974121935111111110e-13 * t) * + t) * + t) * + t) * + t) * + t; + } + case 1: { + double t = 2 * y100 - 3; + return 0.85927161243940350562e-2 + + (0.29085312941641339862e-2 + + (0.15106783707725582090e-4 + + (0.11716709978531327367e-6 + + (0.10197387816021040024e-8 + + (0.10122678863073360769e-10 + 0.10917479678400000000e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 2: { + double t = 2 * y100 - 5; + return 0.14471159831187703054e-1 + + (0.29703978970263836210e-2 + + (0.15835096760173030976e-4 + + (0.12574803383199211596e-6 + + (0.11278672159518415848e-8 + + (0.11547462300333495797e-10 + 0.12894535335111111111e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 3: { + double t = 2 * y100 - 7; + return 0.20476320420324610618e-1 + + (0.30352843012898665856e-2 + + (0.16617609387003727409e-4 + + (0.13525429711163116103e-6 + + (0.12515095552507169013e-8 + + (0.13235687543603382345e-10 + 0.15326595042666666667e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 4: { + double t = 2 * y100 - 9; + return 0.26614461952489004566e-1 + + (0.31034189276234947088e-2 + + (0.17460268109986214274e-4 + + (0.14582130824485709573e-6 + + (0.13935959083809746345e-8 + + (0.15249438072998932900e-10 + 0.18344741882133333333e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 5: { + double t = 2 * y100 - 11; + return 0.32892330248093586215e-1 + + (0.31750557067975068584e-2 + + (0.18369907582308672632e-4 + + (0.15761063702089457882e-6 + + (0.15577638230480894382e-8 + + (0.17663868462699097951e-10 + + (0.22126732680711111111e-12 + 0.30273474177737853668e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 6: { + double t = 2 * y100 - 13; + return 0.39317207681134336024e-1 + + (0.32504779701937539333e-2 + + (0.19354426046513400534e-4 + + (0.17081646971321290539e-6 + + (0.17485733959327106250e-8 + + (0.20593687304921961410e-10 + + (0.26917401949155555556e-12 + 0.38562123837725712270e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 7: { + double t = 2 * y100 - 15; + return 0.45896976511367738235e-1 + + (0.33300031273110976165e-2 + + (0.20423005398039037313e-4 + + (0.18567412470376467303e-6 + + (0.19718038363586588213e-8 + + (0.24175006536781219807e-10 + + (0.33059982791466666666e-12 + 0.49756574284439426165e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 8: { + double t = 2 * y100 - 17; + return 0.52640192524848962855e-1 + + (0.34139883358846720806e-2 + + (0.21586390240603337337e-4 + + (0.20247136501568904646e-6 + + (0.22348696948197102935e-8 + + (0.28597516301950162548e-10 + + (0.41045502119111111110e-12 + 0.65151614515238361946e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 9: { + double t = 2 * y100 - 19; + return 0.59556171228656770456e-1 + + (0.35028374386648914444e-2 + + (0.22857246150998562824e-4 + + (0.22156372146525190679e-6 + + (0.25474171590893813583e-8 + + (0.34122390890697400584e-10 + + (0.51593189879111111110e-12 + 0.86775076853908006938e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 10: { + double t = 2 * y100 - 21; + return 0.66655089485108212551e-1 + + (0.35970095381271285568e-2 + + (0.24250626164318672928e-4 + + (0.24339561521785040536e-6 + + (0.29221990406518411415e-8 + + (0.41117013527967776467e-10 + + (0.65786450716444444445e-12 + 0.11791885745450623331e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 11: { + double t = 2 * y100 - 23; + return 0.73948106345519174661e-1 + + (0.36970297216569341748e-2 + + (0.25784588137312868792e-4 + + (0.26853012002366752770e-6 + + (0.33763958861206729592e-8 + + (0.50111549981376976397e-10 + + (0.85313857496888888890e-12 + 0.16417079927706899860e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 12: { + double t = 2 * y100 - 25; + return 0.81447508065002963203e-1 + + (0.38035026606492705117e-2 + + (0.27481027572231851896e-4 + + (0.29769200731832331364e-6 + + (0.39336816287457655076e-8 + + (0.61895471132038157624e-10 + + (0.11292303213511111111e-11 + 0.23558532213703884304e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 13: { + double t = 2 * y100 - 27; + return 0.89166884027582716628e-1 + + (0.39171301322438946014e-2 + + (0.29366827260422311668e-4 + + (0.33183204390350724895e-6 + + (0.46276006281647330524e-8 + + (0.77692631378169813324e-10 + + (0.15335153258844444444e-11 + 0.35183103415916026911e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 14: { + double t = 2 * y100 - 29; + return 0.97121342888032322019e-1 + + (0.40387340353207909514e-2 + + (0.31475490395950776930e-4 + + (0.37222714227125135042e-6 + + (0.55074373178613809996e-8 + + (0.99509175283990337944e-10 + + (0.21552645758222222222e-11 + 0.55728651431872687605e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 15: { + double t = 2 * y100 - 31; + return 0.10532778218603311137e0 + + (0.41692873614065380607e-2 + + (0.33849549774889456984e-4 + + (0.42064596193692630143e-6 + + (0.66494579697622432987e-8 + + (0.13094103581931802337e-9 + + (0.31896187409777777778e-11 + 0.97271974184476560742e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 16: { + double t = 2 * y100 - 33; + return 0.11380523107427108222e0 + + (0.43099572287871821013e-2 + + (0.36544324341565929930e-4 + + (0.47965044028581857764e-6 + + (0.81819034238463698796e-8 + + (0.17934133239549647357e-9 + + (0.50956666166186293627e-11 + + (0.18850487318190638010e-12 + 0.79697813173519853340e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 17: { + double t = 2 * y100 - 35; + return 0.12257529703447467345e0 + + (0.44621675710026986366e-2 + + (0.39634304721292440285e-4 + + (0.55321553769873381819e-6 + + (0.10343619428848520870e-7 + + (0.26033830170470368088e-9 + + (0.87743837749108025357e-11 + + (0.34427092430230063401e-12 + 0.10205506615709843189e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 18: { + double t = 2 * y100 - 37; + return 0.13166276955656699478e0 + + (0.46276970481783001803e-2 + + (0.43225026380496399310e-4 + + (0.64799164020016902656e-6 + + (0.13580082794704641782e-7 + + (0.39839800853954313927e-9 + + (0.14431142411840000000e-10 + 0.42193457308830027541e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 19: { + double t = 2 * y100 - 39; + return 0.14109647869803356475e0 + + (0.48088424418545347758e-2 + + (0.47474504753352150205e-4 + + (0.77509866468724360352e-6 + + (0.18536851570794291724e-7 + + (0.60146623257887570439e-9 + + (0.18533978397305276318e-10 + + (0.41033845938901048380e-13 - 0.46160680279304825485e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 20: { + double t = 2 * y100 - 41; + return 0.15091057940548936603e0 + + (0.50086864672004685703e-2 + + (0.52622482832192230762e-4 + + (0.95034664722040355212e-6 + + (0.25614261331144718769e-7 + + (0.80183196716888606252e-9 + + (0.12282524750534352272e-10 + (-0.10531774117332273617e-11 - + 0.86157181395039646412e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 21: { + double t = 2 * y100 - 43; + return 0.16114648116017010770e0 + + (0.52314661581655369795e-2 + + (0.59005534545908331315e-4 + + (0.11885518333915387760e-5 + + (0.33975801443239949256e-7 + + (0.82111547144080388610e-9 + (-0.12357674017312854138e-10 + + (-0.24355112256914479176e-11 - + 0.75155506863572930844e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 22: { + double t = 2 * y100 - 45; + return 0.17185551279680451144e0 + + (0.54829002967599420860e-2 + + (0.67013226658738082118e-4 + + (0.14897400671425088807e-5 + + (0.40690283917126153701e-7 + + (0.44060872913473778318e-9 + + (-0.52641873433280000000e-10 - 0.30940587864543343124e-11 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 23: { + double t = 2 * y100 - 47; + return 0.18310194559815257381e0 + + (0.57701559375966953174e-2 + + (0.76948789401735193483e-4 + + (0.18227569842290822512e-5 + + (0.41092208344387212276e-7 + + (-0.44009499965694442143e-9 + (-0.92195414685628803451e-10 + + (-0.22657389705721753299e-11 + + 0.10004784908106839254e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 24: { + double t = 2 * y100 - 49; + return 0.19496527191546630345e0 + + (0.61010853144364724856e-2 + + (0.88812881056342004864e-4 + + (0.21180686746360261031e-5 + + (0.30652145555130049203e-7 + + (-0.16841328574105890409e-8 + + (-0.11008129460612823934e-9 + (-0.12180794204544515779e-12 + + 0.15703325634590334097e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 25: { + double t = 2 * y100 - 51; + return 0.20754006813966575720e0 + + (0.64825787724922073908e-2 + + (0.10209599627522311893e-3 + + (0.22785233392557600468e-5 + + (0.73495224449907568402e-8 + + (-0.29442705974150112783e-8 + + (-0.94082603434315016546e-10 + + (0.23609990400179321267e-11 + 0.14141908654269023788e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 26: { + double t = 2 * y100 - 53; + return 0.22093185554845172146e0 + + (0.69182878150187964499e-2 + + (0.11568723331156335712e-3 + + (0.22060577946323627739e-5 + + (-0.26929730679360840096e-7 + + (-0.38176506152362058013e-8 + + (-0.47399503861054459243e-10 + + (0.40953700187172127264e-11 + 0.69157730376118511127e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 27: { + double t = 2 * y100 - 55; + return 0.23524827304057813918e0 + + (0.74063350762008734520e-2 + + (0.12796333874615790348e-3 + + (0.18327267316171054273e-5 + + (-0.66742910737957100098e-7 + + (-0.40204740975496797870e-8 + + (0.14515984139495745330e-10 + + (0.44921608954536047975e-11 - 0.18583341338983776219e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 28: { + double t = 2 * y100 - 57; + return 0.25058626331812744775e0 + + (0.79377285151602061328e-2 + + (0.13704268650417478346e-3 + + (0.11427511739544695861e-5 + + (-0.10485442447768377485e-6 + + (-0.34850364756499369763e-8 + + (0.72656453829502179208e-10 + + (0.36195460197779299406e-11 - 0.84882136022200714710e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 29: { + double t = 2 * y100 - 59; + return 0.26701724900280689785e0 + + (0.84959936119625864274e-2 + + (0.14112359443938883232e-3 + + (0.17800427288596909634e-6 + + (-0.13443492107643109071e-6 + + (-0.23512456315677680293e-8 + + (0.11245846264695936769e-9 + + (0.19850501334649565404e-11 - 0.11284666134635050832e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 30: { + double t = 2 * y100 - 61; + return 0.28457293586253654144e0 + + (0.90581563892650431899e-2 + + (0.13880520331140646738e-3 + + (-0.97262302362522896157e-6 + + (-0.15077100040254187366e-6 + + (-0.88574317464577116689e-9 + + (0.12760311125637474581e-9 + + (0.20155151018282695055e-12 - 0.10514169375181734921e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 31: { + double t = 2 * y100 - 63; + return 0.30323425595617385705e0 + + (0.95968346790597422934e-2 + + (0.12931067776725883939e-3 + + (-0.21938741702795543986e-5 + + (-0.15202888584907373963e-6 + + (0.61788350541116331411e-9 + + (0.11957835742791248256e-9 + (-0.12598179834007710908e-11 - + 0.75151817129574614194e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 32: { + double t = 2 * y100 - 65; + return 0.32292521181517384379e0 + + (0.10082957727001199408e-1 + + (0.11257589426154962226e-3 + + (-0.33670890319327881129e-5 + + (-0.13910529040004008158e-6 + + (0.19170714373047512945e-8 + + (0.94840222377720494290e-10 + (-0.21650018351795353201e-11 - + 0.37875211678024922689e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 33: { + double t = 2 * y100 - 67; + return 0.34351233557911753862e0 + + (0.10488575435572745309e-1 + + (0.89209444197248726614e-4 + + (-0.43893459576483345364e-5 + + (-0.11488595830450424419e-6 + + (0.28599494117122464806e-8 + + (0.61537542799857777779e-10 - 0.24935749227658002212e-11 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 34: { + double t = 2 * y100 - 69; + return 0.36480946642143669093e0 + + (0.10789304203431861366e-1 + + (0.60357993745283076834e-4 + + (-0.51855862174130669389e-5 + + (-0.83291664087289801313e-7 + + (0.33898011178582671546e-8 + + (0.27082948188277716482e-10 + (-0.23603379397408694974e-11 + + 0.19328087692252869842e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 35: { + double t = 2 * y100 - 71; + return 0.38658679935694939199e0 + + (0.10966119158288804999e-1 + + (0.27521612041849561426e-4 + + (-0.57132774537670953638e-5 + + (-0.48404772799207914899e-7 + + (0.35268354132474570493e-8 + (-0.32383477652514618094e-11 + + (-0.19334202915190442501e-11 + + 0.32333189861286460270e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 36: { + double t = 2 * y100 - 73; + return 0.40858275583808707870e0 + + (0.11006378016848466550e-1 + + (-0.76396376685213286033e-5 + + (-0.59609835484245791439e-5 + + (-0.13834610033859313213e-7 + + (0.33406952974861448790e-8 + (-0.26474915974296612559e-10 + + (-0.13750229270354351983e-11 + + 0.36169366979417390637e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 37: { + double t = 2 * y100 - 75; + return 0.43051714914006682977e0 + + (0.10904106549500816155e-1 + + (-0.43477527256787216909e-4 + + (-0.59429739547798343948e-5 + + (0.17639200194091885949e-7 + + (0.29235991689639918688e-8 + (-0.41718791216277812879e-10 + + (-0.81023337739508049606e-12 + + 0.33618915934461994428e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 38: { + double t = 2 * y100 - 77; + return 0.45210428135559607406e0 + + (0.10659670756384400554e-1 + + (-0.78488639913256978087e-4 + + (-0.56919860886214735936e-5 + + (0.44181850467477733407e-7 + + (0.23694306174312688151e-8 + (-0.49492621596685443247e-10 + + (-0.31827275712126287222e-12 + + 0.27494438742721623654e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 39: { + double t = 2 * y100 - 79; + return 0.47306491195005224077e0 + + (0.10279006119745977570e-1 + + (-0.11140268171830478306e-3 + + (-0.52518035247451432069e-5 + + (0.64846898158889479518e-7 + + (0.17603624837787337662e-8 + + (-0.51129481592926104316e-10 + + (0.62674584974141049511e-13 + 0.20055478560829935356e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 40: { + double t = 2 * y100 - 81; + return 0.49313638965719857647e0 + + (0.97725799114772017662e-2 + + (-0.14122854267291533334e-3 + + (-0.46707252568834951907e-5 + + (0.79421347979319449524e-7 + + (0.11603027184324708643e-8 + + (-0.48269605844397175946e-10 + + (0.32477251431748571219e-12 + 0.12831052634143527985e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 41: { + double t = 2 * y100 - 83; + return 0.51208057433416004042e0 + + (0.91542422354009224951e-2 + + (-0.16726530230228647275e-3 + + (-0.39964621752527649409e-5 + + (0.88232252903213171454e-7 + + (0.61343113364949928501e-9 + + (-0.42516755603130443051e-10 + + (0.47910437172240209262e-12 + 0.66784341874437478953e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 42: { + double t = 2 * y100 - 85; + return 0.52968945458607484524e0 + + (0.84400880445116786088e-2 + + (-0.18908729783854258774e-3 + + (-0.32725905467782951931e-5 + + (0.91956190588652090659e-7 + + (0.14593989152420122909e-9 + + (-0.35239490687644444445e-10 + 0.54613829888448694898e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 43: { + double t = 2 * y100 - 87; + return 0.54578857454330070965e0 + + (0.76474155195880295311e-2 + + (-0.20651230590808213884e-3 + + (-0.25364339140543131706e-5 + + (0.91455367999510681979e-7 + + (-0.23061359005297528898e-9 + + (-0.27512928625244444444e-10 + 0.54895806008493285579e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 44: { + double t = 2 * y100 - 89; + return 0.56023851910298493910e0 + + (0.67938321739997196804e-2 + + (-0.21956066613331411760e-3 + + (-0.18181127670443266395e-5 + + (0.87650335075416845987e-7 + + (-0.51548062050366615977e-9 + + (-0.20068462174044444444e-10 + 0.50912654909758187264e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 45: { + double t = 2 * y100 - 91; + return 0.57293478057455721150e0 + + (0.58965321010394044087e-2 + + (-0.22841145229276575597e-3 + + (-0.11404605562013443659e-5 + + (0.81430290992322326296e-7 + + (-0.71512447242755357629e-9 + + (-0.13372664928000000000e-10 + 0.44461498336689298148e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 46: { + double t = 2 * y100 - 93; + return 0.58380635448407827360e0 + + (0.49717469530842831182e-2 + + (-0.23336001540009645365e-3 + + (-0.51952064448608850822e-6 + + (0.73596577815411080511e-7 + + (-0.84020916763091566035e-9 + + (-0.76700972702222222221e-11 + 0.36914462807972467044e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 47: { + double t = 2 * y100 - 95; + return 0.59281340237769489597e0 + + (0.40343592069379730568e-2 + + (-0.23477963738658326185e-3 + + (0.34615944987790224234e-7 + + (0.64832803248395814574e-7 + + (-0.90329163587627007971e-9 + + (-0.30421940400000000000e-11 + 0.29237386653743536669e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 48: { + double t = 2 * y100 - 97; + return 0.59994428743114271918e0 + + (0.30976579788271744329e-2 + + (-0.23308875765700082835e-3 + + (0.51681681023846925160e-6 + + (0.55694594264948268169e-7 + + (-0.91719117313243464652e-9 + + (0.53982743680000000000e-12 + 0.22050829296187771142e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 49: { + double t = 2 * y100 - 99; + return 0.60521224471819875444e0 + + (0.21732138012345456060e-2 + + (-0.22872428969625997456e-3 + + (0.92588959922653404233e-6 + + (0.46612665806531930684e-7 + + (-0.89393722514414153351e-9 + + (0.31718550353777777778e-11 + 0.15705458816080549117e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 50: { + double t = 2 * y100 - 101; + return 0.60865189969791123620e0 + + (0.12708480848877451719e-2 + + (-0.22212090111534847166e-3 + + (0.12636236031532793467e-5 + + (0.37904037100232937574e-7 + + (-0.84417089968101223519e-9 + + (0.49843180828444444445e-11 + 0.10355439441049048273e-12 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 51: { + double t = 2 * y100 - 103; + return 0.61031580103499200191e0 + + (0.39867436055861038223e-3 + + (-0.21369573439579869291e-3 + + (0.15339402129026183670e-5 + + (0.29787479206646594442e-7 + + (-0.77687792914228632974e-9 + + (0.61192452741333333334e-11 + 0.60216691829459295780e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 52: { + double t = 2 * y100 - 105; + return 0.61027109047879835868e0 + + (-0.43680904508059878254e-3 + + (-0.20383783788303894442e-3 + + (0.17421743090883439959e-5 + + (0.22400425572175715576e-7 + + (-0.69934719320045128997e-9 + + (0.67152759655111111110e-11 + 0.26419960042578359995e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 53: { + double t = 2 * y100 - 107; + return 0.60859639489217430521e0 + + (-0.12305921390962936873e-2 + + (-0.19290150253894682629e-3 + + (0.18944904654478310128e-5 + + (0.15815530398618149110e-7 + + (-0.61726850580964876070e-9 + 0.68987888999111111110e-11 * t) * + t) * + t) * + t) * + t) * + t; + } + case 54: { + double t = 2 * y100 - 109; + return 0.60537899426486075181e0 + + (-0.19790062241395705751e-2 + + (-0.18120271393047062253e-3 + + (0.19974264162313241405e-5 + + (0.10055795094298172492e-7 + + (-0.53491997919318263593e-9 + + (0.67794550295111111110e-11 - 0.17059208095741511603e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 55: { + double t = 2 * y100 - 111; + return 0.60071229457904110537e0 + + (-0.26795676776166354354e-2 + + (-0.16901799553627508781e-3 + + (0.20575498324332621581e-5 + + (0.51077165074461745053e-8 + + (-0.45536079828057221858e-9 + + (0.64488005516444444445e-11 - 0.29311677573152766338e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 56: { + double t = 2 * y100 - 113; + return 0.59469361520112714738e0 + + (-0.33308208190600993470e-2 + + (-0.15658501295912405679e-3 + + (0.20812116912895417272e-5 + + (0.93227468760614182021e-9 + + (-0.38066673740116080415e-9 + + (0.59806790359111111110e-11 - 0.36887077278950440597e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 57: { + double t = 2 * y100 - 115; + return 0.58742228631775388268e0 + + (-0.39321858196059227251e-2 + + (-0.14410441141450122535e-3 + + (0.20743790018404020716e-5 + + (-0.25261903811221913762e-8 + + (-0.31212416519526924318e-9 + + (0.54328422462222222221e-11 - 0.40864152484979815972e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 58: { + double t = 2 * y100 - 117; + return 0.57899804200033018447e0 + + (-0.44838157005618913447e-2 + + (-0.13174245966501437965e-3 + + (0.20425306888294362674e-5 + + (-0.53330296023875447782e-8 + + (-0.25041289435539821014e-9 + + (0.48490437205333333334e-11 - 0.42162206939169045177e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 59: { + double t = 2 * y100 - 119; + return 0.56951968796931245974e0 + + (-0.49864649488074868952e-2 + + (-0.11963416583477567125e-3 + + (0.19906021780991036425e-5 + + (-0.75580140299436494248e-8 + + (-0.19576060961919820491e-9 + + (0.42613011928888888890e-11 - 0.41539443304115604377e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 60: { + double t = 2 * y100 - 121; + return 0.55908401930063918964e0 + + (-0.54413711036826877753e-2 + + (-0.10788661102511914628e-3 + + (0.19229663322982839331e-5 + + (-0.92714731195118129616e-8 + + (-0.14807038677197394186e-9 + + (0.36920870298666666666e-11 - 0.39603726688419162617e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 61: { + double t = 2 * y100 - 123; + return 0.54778496152925675315e0 + + (-0.58501497933213396670e-2 + + (-0.96582314317855227421e-4 + + (0.18434405235069270228e-5 + + (-0.10541580254317078711e-7 + + (-0.10702303407788943498e-9 + + (0.31563175582222222222e-11 - 0.36829748079110481422e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 62: { + double t = 2 * y100 - 125; + return 0.53571290831682823999e0 + + (-0.62147030670760791791e-2 + + (-0.85782497917111760790e-4 + + (0.17553116363443470478e-5 + + (-0.11432547349815541084e-7 + + (-0.72157091369041330520e-10 + + (0.26630811607111111111e-11 - 0.33578660425893164084e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 63: { + double t = 2 * y100 - 127; + return 0.52295422962048434978e0 + + (-0.65371404367776320720e-2 + + (-0.75530164941473343780e-4 + + (0.16613725797181276790e-5 + + (-0.12003521296598910761e-7 + + (-0.42929753689181106171e-10 + + (0.22170894940444444444e-11 - 0.30117697501065110505e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 64: { + double t = 2 * y100 - 129; + return 0.50959092577577886140e0 + + (-0.68197117603118591766e-2 + + (-0.65852936198953623307e-4 + + (0.15639654113906716939e-5 + + (-0.12308007991056524902e-7 + + (-0.18761997536910939570e-10 + + (0.18198628922666666667e-11 - 0.26638355362285200932e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 65: { + double t = 2 * y100 - 131; + return 0.49570040481823167970e0 + + (-0.70647509397614398066e-2 + + (-0.56765617728962588218e-4 + + (0.14650274449141448497e-5 + + (-0.12393681471984051132e-7 + + (0.92904351801168955424e-12 + + (0.14706755960177777778e-11 - 0.23272455351266325318e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 66: { + double t = 2 * y100 - 133; + return 0.48135536250935238066e0 + + (-0.72746293327402359783e-2 + + (-0.48272489495730030780e-4 + + (0.13661377309113939689e-5 + + (-0.12302464447599382189e-7 + + (0.16707760028737074907e-10 + + (0.11672928324444444444e-11 - 0.20105801424709924499e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 67: { + double t = 2 * y100 - 135; + return 0.46662374675511439448e0 + + (-0.74517177649528487002e-2 + + (-0.40369318744279128718e-4 + + (0.12685621118898535407e-5 + + (-0.12070791463315156250e-7 + + (0.29105507892605823871e-10 + + (0.90653314645333333334e-12 - 0.17189503312102982646e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 68: { + double t = 2 * y100 - 137; + return 0.45156879030168268778e0 + + (-0.75983560650033817497e-2 + + (-0.33045110380705139759e-4 + + (0.11732956732035040896e-5 + + (-0.11729986947158201869e-7 + + (0.38611905704166441308e-10 + + (0.68468768305777777779e-12 - 0.14549134330396754575e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 69: { + double t = 2 * y100 - 139; + return 0.43624909769330896904e0 + + (-0.77168291040309554679e-2 + + (-0.26283612321339907756e-4 + + (0.10811018836893550820e-5 + + (-0.11306707563739851552e-7 + + (0.45670446788529607380e-10 + + (0.49782492549333333334e-12 - 0.12191983967561779442e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 70: { + double t = 2 * y100 - 141; + return 0.42071877443548481181e0 + + (-0.78093484015052730097e-2 + + (-0.20064596897224934705e-4 + + (0.99254806680671890766e-6 + + (-0.10823412088884741451e-7 + + (0.50677203326904716247e-10 + + (0.34200547594666666666e-12 - 0.10112698698356194618e-13 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 71: { + double t = 2 * y100 - 143; + return 0.40502758809710844280e0 + + (-0.78780384460872937555e-2 + + (-0.14364940764532853112e-4 + + (0.90803709228265217384e-6 + + (-0.10298832847014466907e-7 + + (0.53981671221969478551e-10 + + (0.21342751381333333333e-12 - 0.82975901848387729274e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 72: { + double t = 2 * y100 - 145; + return 0.38922115269731446690e0 + + (-0.79249269708242064120e-2 + + (-0.91595258799106970453e-5 + + (0.82783535102217576495e-6 + + (-0.97484311059617744437e-8 + + (0.55889029041660225629e-10 + + (0.10851981336888888889e-12 - 0.67278553237853459757e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 73: { + double t = 2 * y100 - 147; + return 0.37334112915460307335e0 + + (-0.79519385109223148791e-2 + + (-0.44219833548840469752e-5 + + (0.75209719038240314732e-6 + + (-0.91848251458553190451e-8 + + (0.56663266668051433844e-10 + + (0.23995894257777777778e-13 - 0.53819475285389344313e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 74: { + double t = 2 * y100 - 149; + return 0.35742543583374223085e0 + + (-0.79608906571527956177e-2 + + (-0.12530071050975781198e-6 + + (0.68088605744900552505e-6 + + (-0.86181844090844164075e-8 + + (0.56530784203816176153e-10 + + (-0.43120012248888888890e-13 - 0.42372603392496813810e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 75: { + double t = 2 * y100 - 151; + return 0.34150846431979618536e0 + + (-0.79534924968773806029e-2 + + (0.37576885610891515813e-5 + + (0.61419263633090524326e-6 + + (-0.80565865409945960125e-8 + + (0.55684175248749269411e-10 + + (-0.95486860764444444445e-13 - 0.32712946432984510595e-14 * t) * + t) * + t) * + t) * + t) * + t) * + t; + } + case 76: { + double t = 2 * y100 - 153; + return 0.32562129649136346824e0 + + (-0.79313448067948884309e-2 + + (0.72539159933545300034e-5 + + (0.55195028297415503083e-6 + + (-0.75063365335570475258e-8 + + (0.54281686749699595941e-10 - 0.13545424295111111111e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 77: { + double t = 2 * y100 - 155; + return 0.30979191977078391864e0 + + (-0.78959416264207333695e-2 + + (0.10389774377677210794e-4 + + (0.49404804463196316464e-6 + + (-0.69722488229411164685e-8 + + (0.52469254655951393842e-10 - 0.16507860650666666667e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 78: { + double t = 2 * y100 - 157; + return 0.29404543811214459904e0 + + (-0.78486728990364155356e-2 + + (0.13190885683106990459e-4 + + (0.44034158861387909694e-6 + + (-0.64578942561562616481e-8 + + (0.50354306498006928984e-10 - 0.18614473550222222222e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 79: { + double t = 2 * y100 - 159; + return 0.27840427686253660515e0 + + (-0.77908279176252742013e-2 + + (0.15681928798708548349e-4 + + (0.39066226205099807573e-6 + + (-0.59658144820660420814e-8 + + (0.48030086420373141763e-10 - 0.20018995173333333333e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 80: { + double t = 2 * y100 - 161; + return 0.26288838011163800908e0 + + (-0.77235993576119469018e-2 + + (0.17886516796198660969e-4 + + (0.34482457073472497720e-6 + + (-0.54977066551955420066e-8 + + (0.45572749379147269213e-10 - 0.20852924954666666667e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 81: { + double t = 2 * y100 - 163; + return 0.24751539954181029717e0 + + (-0.76480877165290370975e-2 + + (0.19827114835033977049e-4 + + (0.30263228619976332110e-6 + + (-0.50545814570120129947e-8 + + (0.43043879374212005966e-10 - 0.21228012028444444444e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 82: { + double t = 2 * y100 - 165; + return 0.23230087411688914593e0 + + (-0.75653060136384041587e-2 + + (0.21524991113020016415e-4 + + (0.26388338542539382413e-6 + + (-0.46368974069671446622e-8 + + (0.40492715758206515307e-10 - 0.21238627815111111111e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 83: { + double t = 2 * y100 - 167; + return 0.21725840021297341931e0 + + (-0.74761846305979730439e-2 + + (0.23000194404129495243e-4 + + (0.22837400135642906796e-6 + + (-0.42446743058417541277e-8 + + (0.37958104071765923728e-10 - 0.20963978568888888889e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 84: { + double t = 2 * y100 - 169; + return 0.20239979200788191491e0 + + (-0.73815761980493466516e-2 + + (0.24271552727631854013e-4 + + (0.19590154043390012843e-6 + + (-0.38775884642456551753e-8 + + (0.35470192372162901168e-10 - 0.20470131678222222222e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 85: { + double t = 2 * y100 - 171; + return 0.18773523211558098962e0 + + (-0.72822604530339834448e-2 + + (0.25356688567841293697e-4 + + (0.16626710297744290016e-6 + + (-0.35350521468015310830e-8 + + (0.33051896213898864306e-10 - 0.19811844544000000000e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 86: { + double t = 2 * y100 - 173; + return 0.17327341258479649442e0 + + (-0.71789490089142761950e-2 + + (0.26272046822383820476e-4 + + (0.13927732375657362345e-6 + + (-0.32162794266956859603e-8 + + (0.30720156036105652035e-10 - 0.19034196304000000000e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 87: { + double t = 2 * y100 - 175; + return 0.15902166648328672043e0 + + (-0.70722899934245504034e-2 + + (0.27032932310132226025e-4 + + (0.11474573347816568279e-6 + + (-0.29203404091754665063e-8 + + (0.28487010262547971859e-10 - 0.18174029063111111111e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 88: { + double t = 2 * y100 - 177; + return 0.14498609036610283865e0 + + (-0.69628725220045029273e-2 + + (0.27653554229160596221e-4 + + (0.92493727167393036470e-7 + + (-0.26462055548683583849e-8 + + (0.26360506250989943739e-10 - 0.17261211260444444444e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 89: { + double t = 2 * y100 - 179; + return 0.13117165798208050667e0 + + (-0.68512309830281084723e-2 + + (0.28147075431133863774e-4 + + (0.72351212437979583441e-7 + + (-0.23927816200314358570e-8 + + (0.24345469651209833155e-10 - 0.16319736960000000000e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 90: { + double t = 2 * y100 - 181; + return 0.11758232561160626306e0 + + (-0.67378491192463392927e-2 + + (0.28525664781722907847e-4 + + (0.54156999310046790024e-7 + + (-0.21589405340123827823e-8 + + (0.22444150951727334619e-10 - 0.15368675584000000000e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 91: { + double t = 2 * y100 - 183; + return 0.10422112945361673560e0 + + (-0.66231638959845581564e-2 + + (0.28800551216363918088e-4 + + (0.37758983397952149613e-7 + + (-0.19435423557038933431e-8 + + (0.20656766125421362458e-10 - 0.14422990012444444444e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 92: { + double t = 2 * y100 - 185; + return 0.91090275493541084785e-1 + + (-0.65075691516115160062e-2 + + (0.28982078385527224867e-4 + + (0.23014165807643012781e-7 + + (-0.17454532910249875958e-8 + + (0.18981946442680092373e-10 - 0.13494234691555555556e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 93: { + double t = 2 * y100 - 187; + return 0.78191222288771379358e-1 + + (-0.63914190297303976434e-2 + + (0.29079759021299682675e-4 + + (0.97885458059415717014e-8 + + (-0.15635596116134296819e-8 + + (0.17417110744051331974e-10 - 0.12591151763555555556e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 94: { + double t = 2 * y100 - 189; + return 0.65524757106147402224e-1 + + (-0.62750311956082444159e-2 + + (0.29102328354323449795e-4 + + (-0.20430838882727954582e-8 + + (-0.13967781903855367270e-8 + + (0.15958771833747057569e-10 - 0.11720175765333333333e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 95: { + double t = 2 * y100 - 191; + return 0.53091065838453612773e-1 + + (-0.61586898417077043662e-2 + + (0.29057796072960100710e-4 + + (-0.12597414620517987536e-7 + + (-0.12440642607426861943e-8 + + (0.14602787128447932137e-10 - 0.10885859114666666667e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 96: { + double t = 2 * y100 - 193; + return 0.40889797115352738582e-1 + + (-0.60426484889413678200e-2 + + (0.28953496450191694606e-4 + + (-0.21982952021823718400e-7 + + (-0.11044169117553026211e-8 + + (0.13344562332430552171e-10 - 0.10091231402844444444e-12 * t) * + t) * + t) * + t) * + t) * + t; + } + case 97: + case 98: + case 99: + case 100: { // use Taylor expansion for small x (|x| <= 0.0309...) + // (2/sqrt(pi)) * (x - 2/3 x^3 + 4/15 x^5 - 8/105 x^7 + 16/945 x^9) + double x2 = x * x; + return x * (1.1283791670955125739 - + x2 * (0.75225277806367504925 - + x2 * (0.30090111122547001970 - + x2 * (0.085971746064420005629 - + x2 * 0.016931216931216931217)))); + } + } + /* Since 0 <= y100 < 101, this is only reached if x is NaN, + in which case we should return NaN. */ + return NaN; +} + +double FADDEEVA(w_im)(double x) { + if (x >= 0) { + if (x > 45) { // continued-fraction expansion is faster + const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) + if (x > 5e7) // 1-term expansion, important to avoid overflow + return ispi / x; + /* 5-term expansion (rely on compiler for CSE), simplified from: + ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */ + return ispi * ((x * x) * (x * x - 4.5) + 2) / + (x * ((x * x) * (x * x - 5) + 3.75)); + } + return w_im_y100(100 / (1 + x), x); + } else { // = -FADDEEVA(w_im)(-x) + if (x < -45) { // continued-fraction expansion is faster + const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) + if (x < -5e7) // 1-term expansion, important to avoid overflow + return ispi / x; + /* 5-term expansion (rely on compiler for CSE), simplified from: + ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */ + return ispi * ((x * x) * (x * x - 4.5) + 2) / + (x * ((x * x) * (x * x - 5) + 3.75)); + } + return -w_im_y100(100 / (1 - x), -x); + } +} + +///////////////////////////////////////////////////////////////////////// + +// Compile with -DTEST_FADDEEVA to compile a little test program +#ifdef TEST_FADDEEVA + +#ifdef __cplusplus +#include +#else +#include +#endif + +// compute relative error |b-a|/|a|, handling case of NaN and Inf, +static double relerr(double a, double b) { + if (isnan(a) || isnan(b) || isinf(a) || isinf(b)) { + if ((isnan(a) && !isnan(b)) || (!isnan(a) && isnan(b)) || + (isinf(a) && !isinf(b)) || (!isinf(a) && isinf(b)) || + (isinf(a) && isinf(b) && a * b < 0)) + return Inf; // "infinite" error + return 0; // matching infinity/nan results counted as zero error + } + if (a == 0) + return b == 0 ? 0 : Inf; + else + return fabs((b - a) / a); +} + +int main(void) { + double errmax_all = 0; + { + printf("############# w(z) tests #############\n"); +#define NTST 57 // define instead of const for C compatibility + cmplx z[NTST] = {C(624.2, -0.26123), + C(-0.4, 3.), + C(0.6, 2.), + C(-1., 1.), + C(-1., -9.), + C(-1., 9.), + C(-0.0000000234545, 1.1234), + C(-3., 5.1), + C(-53, 30.1), + C(0.0, 0.12345), + C(11, 1), + C(-22, -2), + C(9, -28), + C(21, -33), + C(1e5, 1e5), + C(1e14, 1e14), + C(-3001, -1000), + C(1e160, -1e159), + C(-6.01, 0.01), + C(-0.7, -0.7), + C(2.611780000000000e+01, 4.540909610972489e+03), + C(0.8e7, 0.3e7), + C(-20, -19.8081), + C(1e-16, -1.1e-16), + C(2.3e-8, 1.3e-8), + C(6.3, -1e-13), + C(6.3, 1e-20), + C(1e-20, 6.3), + C(1e-20, 16.3), + C(9, 1e-300), + C(6.01, 0.11), + C(8.01, 1.01e-10), + C(28.01, 1e-300), + C(10.01, 1e-200), + C(10.01, -1e-200), + C(10.01, 0.99e-10), + C(10.01, -0.99e-10), + C(1e-20, 7.01), + C(-1, 7.01), + C(5.99, 7.01), + C(1, 0), + C(55, 0), + C(-0.1, 0), + C(1e-20, 0), + C(0, 5e-14), + C(0, 51), + C(Inf, 0), + C(-Inf, 0), + C(0, Inf), + C(0, -Inf), + C(Inf, Inf), + C(Inf, -Inf), + C(NaN, NaN), + C(NaN, 0), + C(0, NaN), + C(NaN, Inf), + C(Inf, NaN)}; + cmplx w[NTST] = { + /* w(z), computed with WolframAlpha + ... note that WolframAlpha is problematic + some of the above inputs, so I had to + use the continued-fraction expansion + in WolframAlpha in some cases, or switch + to Maple */ + C(-3.78270245518980507452677445620103199303131110e-7, + 0.000903861276433172057331093754199933411710053155), + C(0.1764906227004816847297495349730234591778719532788, + -0.02146550539468457616788719893991501311573031095617), + C(0.2410250715772692146133539023007113781272362309451, + 0.06087579663428089745895459735240964093522265589350), + C(0.30474420525691259245713884106959496013413834051768, + -0.20821893820283162728743734725471561394145872072738), + C(7.317131068972378096865595229600561710140617977e34, + 8.321873499714402777186848353320412813066170427e34), + C(0.0615698507236323685519612934241429530190806818395, + -0.00676005783716575013073036218018565206070072304635), + C(0.3960793007699874918961319170187598400134746631, + -5.593152259116644920546186222529802777409274656e-9), + C(0.08217199226739447943295069917990417630675021771804, + -0.04701291087643609891018366143118110965272615832184), + C(0.00457246000350281640952328010227885008541748668738, + -0.00804900791411691821818731763401840373998654987934), + C(0.8746342859608052666092782112565360755791467973338452, 0.), + C(0.00468190164965444174367477874864366058339647648741, + 0.0510735563901306197993676329845149741675029197050), + C(-0.0023193175200187620902125853834909543869428763219, + -0.025460054739731556004902057663500272721780776336), + C(9.11463368405637174660562096516414499772662584e304, + 3.97101807145263333769664875189354358563218932e305), + C(-4.4927207857715598976165541011143706155432296e281, + -2.8019591213423077494444700357168707775769028e281), + C(2.820947917809305132678577516325951485807107151e-6, + 2.820947917668257736791638444590253942253354058e-6), + C(2.82094791773878143474039725787438662716372268e-15, + 2.82094791773878143474039725773333923127678361e-15), + C(-0.0000563851289696244350147899376081488003110150498, + -0.000169211755126812174631861529808288295454992688), + C(-5.586035480670854326218608431294778077663867e-162, + 5.586035480670854326218608431294778077663867e-161), + C(0.00016318325137140451888255634399123461580248456, + -0.095232456573009287370728788146686162555021209999), + C(0.69504753678406939989115375989939096800793577783885, + -1.8916411171103639136680830887017670616339912024317), + C(0.0001242418269653279656612334210746733213167234822, + 7.145975826320186888508563111992099992116786763e-7), + C(2.318587329648353318615800865959225429377529825e-8, + 6.182899545728857485721417893323317843200933380e-8), + C(-0.0133426877243506022053521927604277115767311800303, + -0.0148087097143220769493341484176979826888871576145), + C(1.00000000000000012412170838050638522857747934, + 1.12837916709551279389615890312156495593616433e-16), + C(0.9999999853310704677583504063775310832036830015, + 2.595272024519678881897196435157270184030360773e-8), + C(-1.4731421795638279504242963027196663601154624e-15, + 0.090727659684127365236479098488823462473074709), + C(5.79246077884410284575834156425396800754409308e-18, + 0.0907276596841273652364790985059772809093822374), + C(0.0884658993528521953466533278764830881245144368, + 1.37088352495749125283269718778582613192166760e-22), + C(0.0345480845419190424370085249304184266813447878, + 2.11161102895179044968099038990446187626075258e-23), + C(6.63967719958073440070225527042829242391918213e-36, + 0.0630820900592582863713653132559743161572639353), + C(0.00179435233208702644891092397579091030658500743634, + 0.0951983814805270647939647438459699953990788064762), + C(9.09760377102097999924241322094863528771095448e-13, + 0.0709979210725138550986782242355007611074966717), + C(7.2049510279742166460047102593255688682910274423e-304, + 0.0201552956479526953866611812593266285000876784321), + C(3.04543604652250734193622967873276113872279682e-44, + 0.0566481651760675042930042117726713294607499165), + C(3.04543604652250734193622967873276113872279682e-44, + 0.0566481651760675042930042117726713294607499165), + C(0.5659928732065273429286988428080855057102069081e-12, + 0.056648165176067504292998527162143030538756683302), + C(-0.56599287320652734292869884280802459698927645e-12, + 0.0566481651760675042929985271621430305387566833029), + C(0.0796884251721652215687859778119964009569455462, + 1.11474461817561675017794941973556302717225126e-22), + C(0.07817195821247357458545539935996687005781943386550, + -0.01093913670103576690766705513142246633056714279654), + C(0.04670032980990449912809326141164730850466208439937, + 0.03944038961933534137558064191650437353429669886545), + C(0.36787944117144232159552377016146086744581113103176, + 0.60715770584139372911503823580074492116122092866515), + C(0, 0.010259688805536830986089913987516716056946786526145), + C(0.99004983374916805357390597718003655777207908125383, + -0.11208866436449538036721343053869621153527769495574), + C(0.99999999999999999999999999999999999999990000, + 1.12837916709551257389615890312154517168802603e-20), + C(0.999999999999943581041645226871305192054749891144158, 0), + C(0.0110604154853277201542582159216317923453996211744250, 0), + C(0, 0), + C(0, 0), + C(0, 0), + C(Inf, 0), + C(0, 0), + C(NaN, NaN), + C(NaN, NaN), + C(NaN, NaN), + C(NaN, 0), + C(NaN, NaN), + C(NaN, NaN)}; + double errmax = 0; + for (int i = 0; i < NTST; ++i) { + cmplx fw = FADDEEVA(w)(z[i], 0.); + double re_err = relerr(creal(w[i]), creal(fw)); + double im_err = relerr(cimag(w[i]), cimag(fw)); + printf( + "w(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n", + creal(z[i]), cimag(z[i]), creal(fw), cimag(fw), creal(w[i]), + cimag(w[i]), re_err, im_err); + if (re_err > errmax) + errmax = re_err; + if (im_err > errmax) + errmax = im_err; + } + if (errmax > 1e-13) { + printf("FAILURE -- relative error %g too large!\n", errmax); + return 1; + } + printf("SUCCESS (max relative error = %g)\n", errmax); + if (errmax > errmax_all) + errmax_all = errmax; + } + { +#undef NTST +#define NTST 41 // define instead of const for C compatibility + cmplx z[NTST] = {C(1, 2), + C(-1, 2), + C(1, -2), + C(-1, -2), + C(9, -28), + C(21, -33), + C(1e3, 1e3), + C(-3001, -1000), + C(1e160, -1e159), + C(5.1e-3, 1e-8), + C(-4.9e-3, 4.95e-3), + C(4.9e-3, 0.5), + C(4.9e-4, -0.5e1), + C(-4.9e-5, -0.5e2), + C(5.1e-3, 0.5), + C(5.1e-4, -0.5e1), + C(-5.1e-5, -0.5e2), + C(1e-6, 2e-6), + C(0, 2e-6), + C(0, 2), + C(0, 20), + C(0, 200), + C(Inf, 0), + C(-Inf, 0), + C(0, Inf), + C(0, -Inf), + C(Inf, Inf), + C(Inf, -Inf), + C(NaN, NaN), + C(NaN, 0), + C(0, NaN), + C(NaN, Inf), + C(Inf, NaN), + C(1e-3, NaN), + C(7e-2, 7e-2), + C(7e-2, -7e-4), + C(-9e-2, 7e-4), + C(-9e-2, 9e-2), + C(-7e-4, 9e-2), + C(7e-2, 0.9e-2), + C(7e-2, 1.1e-2)}; + cmplx w[NTST] = {// erf(z[i]), evaluated with Maple + C(-0.5366435657785650339917955593141927494421, + -5.049143703447034669543036958614140565553), + C(0.5366435657785650339917955593141927494421, + -5.049143703447034669543036958614140565553), + C(-0.5366435657785650339917955593141927494421, + 5.049143703447034669543036958614140565553), + C(0.5366435657785650339917955593141927494421, + 5.049143703447034669543036958614140565553), + C(0.3359473673830576996788000505817956637777e304, + -0.1999896139679880888755589794455069208455e304), + C(0.3584459971462946066523939204836760283645e278, + 0.3818954885257184373734213077678011282505e280), + C(0.9996020422657148639102150147542224526887, + 0.00002801044116908227889681753993542916894856), + C(-1, 0), + C(1, 0), + C(0.005754683859034800134412990541076554934877, + 0.1128349818335058741511924929801267822634e-7), + C(-0.005529149142341821193633460286828381876955, + 0.005585388387864706679609092447916333443570), + C(0.007099365669981359632319829148438283865814, + 0.6149347012854211635026981277569074001219), + C(0.3981176338702323417718189922039863062440e8, + -0.8298176341665249121085423917575122140650e10), + C(-Inf, -Inf), + C(0.007389128308257135427153919483147229573895, + 0.6149332524601658796226417164791221815139), + C(0.4143671923267934479245651547534414976991e8, + -0.8298168216818314211557046346850921446950e10), + C(-Inf, -Inf), + C(0.1128379167099649964175513742247082845155e-5, + 0.2256758334191777400570377193451519478895e-5), + C(0, 0.2256758334194034158904576117253481476197e-5), + C(0, 18.56480241457555259870429191324101719886), + C(0, 0.1474797539628786202447733153131835124599e173), + C(0, Inf), + C(1, 0), + C(-1, 0), + C(0, Inf), + C(0, -Inf), + C(NaN, NaN), + C(NaN, NaN), + C(NaN, NaN), + C(NaN, 0), + C(0, NaN), + C(NaN, NaN), + C(NaN, NaN), + C(NaN, NaN), + C(0.07924380404615782687930591956705225541145, + 0.07872776218046681145537914954027729115247), + C(0.07885775828512276968931773651224684454495, + -0.0007860046704118224342390725280161272277506), + C(-0.1012806432747198859687963080684978759881, + 0.0007834934747022035607566216654982820299469), + C(-0.1020998418798097910247132140051062512527, + 0.1010030778892310851309082083238896270340), + C(-0.0007962891763147907785684591823889484764272, + 0.1018289385936278171741809237435404896152), + C(0.07886408666470478681566329888615410479530, + 0.01010604288780868961492224347707949372245), + C(0.07886723099940260286824654364807981336591, + 0.01235199327873258197931147306290916629654)}; +#define TST(f, isc) \ + printf("############# " #f "(z) tests #############\n"); \ + double errmax = 0; \ + for (int i = 0; i < NTST; ++i) { \ + cmplx fw = FADDEEVA(f)(z[i], 0.); \ + double re_err = relerr(creal(w[i]), creal(fw)); \ + double im_err = relerr(cimag(w[i]), cimag(fw)); \ + printf(#f \ + "(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n", \ + creal(z[i]), cimag(z[i]), creal(fw), cimag(fw), creal(w[i]), \ + cimag(w[i]), re_err, im_err); \ + if (re_err > errmax) \ + errmax = re_err; \ + if (im_err > errmax) \ + errmax = im_err; \ + } \ + if (errmax > 1e-13) { \ + printf("FAILURE -- relative error %g too large!\n", errmax); \ + return 1; \ + } \ + printf("Checking " #f "(x) special case...\n"); \ + for (int i = 0; i < 10000; ++i) { \ + double x = pow(10., -300. + i * 600. / (10000 - 1)); \ + double re_err = \ + relerr(FADDEEVA_RE(f)(x), creal(FADDEEVA(f)(C(x, x * isc), 0.))); \ + if (re_err > errmax) \ + errmax = re_err; \ + re_err = \ + relerr(FADDEEVA_RE(f)(-x), creal(FADDEEVA(f)(C(-x, x * isc), 0.))); \ + if (re_err > errmax) \ + errmax = re_err; \ + } \ + { \ + double re_err = \ + relerr(FADDEEVA_RE(f)(Inf), creal(FADDEEVA(f)(C(Inf, 0.), 0.))); \ + if (re_err > errmax) \ + errmax = re_err; \ + re_err = \ + relerr(FADDEEVA_RE(f)(-Inf), creal(FADDEEVA(f)(C(-Inf, 0.), 0.))); \ + if (re_err > errmax) \ + errmax = re_err; \ + re_err = relerr(FADDEEVA_RE(f)(NaN), creal(FADDEEVA(f)(C(NaN, 0.), 0.))); \ + if (re_err > errmax) \ + errmax = re_err; \ + } \ + if (errmax > 1e-13) { \ + printf("FAILURE -- relative error %g too large!\n", errmax); \ + return 1; \ + } \ + printf("SUCCESS (max relative error = %g)\n", errmax); \ + if (errmax > errmax_all) \ + errmax_all = errmax + + TST(erf, 1e-20); + } + { + // since erfi just calls through to erf, just one test should + // be sufficient to make sure I didn't screw up the signs or something +#undef NTST +#define NTST 1 // define instead of const for C compatibility + cmplx z[NTST] = {C(1.234, 0.5678)}; + cmplx w[NTST] = {// erfi(z[i]), computed with Maple + C(1.081032284405373149432716643834106923212, + 1.926775520840916645838949402886591180834)}; + TST(erfi, 0); + } + { + // since erfcx just calls through to w, just one test should + // be sufficient to make sure I didn't screw up the signs or something +#undef NTST +#define NTST 1 // define instead of const for C compatibility + cmplx z[NTST] = {C(1.234, 0.5678)}; + cmplx w[NTST] = {// erfcx(z[i]), computed with Maple + C(0.3382187479799972294747793561190487832579, + -0.1116077470811648467464927471872945833154)}; + TST(erfcx, 0); + } + { +#undef NTST +#define NTST 30 // define instead of const for C compatibility + cmplx z[NTST] = { + C(1, 2), C(-1, 2), C(1, -2), C(-1, -2), + C(9, -28), C(21, -33), C(1e3, 1e3), C(-3001, -1000), + C(1e160, -1e159), C(5.1e-3, 1e-8), C(0, 2e-6), C(0, 2), + C(0, 20), C(0, 200), C(2e-6, 0), C(2, 0), + C(20, 0), C(200, 0), C(Inf, 0), C(-Inf, 0), + C(0, Inf), C(0, -Inf), C(Inf, Inf), C(Inf, -Inf), + C(NaN, NaN), C(NaN, 0), C(0, NaN), C(NaN, Inf), + C(Inf, NaN), C(88, 0)}; + cmplx w[NTST] = {// erfc(z[i]), evaluated with Maple + C(1.536643565778565033991795559314192749442, + 5.049143703447034669543036958614140565553), + C(0.4633564342214349660082044406858072505579, + 5.049143703447034669543036958614140565553), + C(1.536643565778565033991795559314192749442, + -5.049143703447034669543036958614140565553), + C(0.4633564342214349660082044406858072505579, + -5.049143703447034669543036958614140565553), + C(-0.3359473673830576996788000505817956637777e304, + 0.1999896139679880888755589794455069208455e304), + C(-0.3584459971462946066523939204836760283645e278, + -0.3818954885257184373734213077678011282505e280), + C(0.0003979577342851360897849852457775473112748, + -0.00002801044116908227889681753993542916894856), + C(2, 0), + C(0, 0), + C(0.9942453161409651998655870094589234450651, + -0.1128349818335058741511924929801267822634e-7), + C(1, -0.2256758334194034158904576117253481476197e-5), + C(1, -18.56480241457555259870429191324101719886), + C(1, -0.1474797539628786202447733153131835124599e173), + C(1, -Inf), + C(0.9999977432416658119838633199332831406314, 0), + C(0.004677734981047265837930743632747071389108, 0), + C(0.5395865611607900928934999167905345604088e-175, 0), + C(0, 0), + C(0, 0), + C(2, 0), + C(1, -Inf), + C(1, Inf), + C(NaN, NaN), + C(NaN, NaN), + C(NaN, NaN), + C(NaN, 0), + C(1, NaN), + C(NaN, NaN), + C(NaN, NaN), + C(0, 0)}; + TST(erfc, 1e-20); + } + { +#undef NTST +#define NTST 48 // define instead of const for C compatibility + cmplx z[NTST] = {C(2, 1), C(-2, 1), + C(2, -1), C(-2, -1), + C(-28, 9), C(33, -21), + C(1e3, 1e3), C(-1000, -3001), + C(1e-8, 5.1e-3), C(4.95e-3, -4.9e-3), + C(5.1e-3, 5.1e-3), C(0.5, 4.9e-3), + C(-0.5e1, 4.9e-4), C(-0.5e2, -4.9e-5), + C(0.5e3, 4.9e-6), C(0.5, 5.1e-3), + C(-0.5e1, 5.1e-4), C(-0.5e2, -5.1e-5), + C(1e-6, 2e-6), C(2e-6, 0), + C(2, 0), C(20, 0), + C(200, 0), C(0, 4.9e-3), + C(0, -5.1e-3), C(0, 2e-6), + C(0, -2), C(0, 20), + C(0, -200), C(Inf, 0), + C(-Inf, 0), C(0, Inf), + C(0, -Inf), C(Inf, Inf), + C(Inf, -Inf), C(NaN, NaN), + C(NaN, 0), C(0, NaN), + C(NaN, Inf), C(Inf, NaN), + C(39, 6.4e-5), C(41, 6.09e-5), + C(4.9e7, 5e-11), C(5.1e7, 4.8e-11), + C(1e9, 2.4e-12), C(1e11, 2.4e-14), + C(1e13, 2.4e-16), C(1e300, 2.4e-303)}; + cmplx w[NTST] = {// dawson(z[i]), evaluated with Maple + C(0.1635394094345355614904345232875688576839, + -0.1531245755371229803585918112683241066853), + C(-0.1635394094345355614904345232875688576839, + -0.1531245755371229803585918112683241066853), + C(0.1635394094345355614904345232875688576839, + 0.1531245755371229803585918112683241066853), + C(-0.1635394094345355614904345232875688576839, + 0.1531245755371229803585918112683241066853), + C(-0.01619082256681596362895875232699626384420, + -0.005210224203359059109181555401330902819419), + C(0.01078377080978103125464543240346760257008, + 0.006866888783433775382193630944275682670599), + C(-0.5808616819196736225612296471081337245459, + 0.6688593905505562263387760667171706325749), + C(Inf, -Inf), + C(0.1000052020902036118082966385855563526705e-7, + 0.005100088434920073153418834680320146441685), + C(0.004950156837581592745389973960217444687524, + -0.004899838305155226382584756154100963570500), + C(0.005100176864319675957314822982399286703798, + 0.005099823128319785355949825238269336481254), + C(0.4244534840871830045021143490355372016428, + 0.002820278933186814021399602648373095266538), + C(-0.1021340733271046543881236523269967674156, + -0.00001045696456072005761498961861088944159916), + C(-0.01000200120119206748855061636187197886859, + 0.9805885888237419500266621041508714123763e-8), + C(0.001000002000012000023960527532953151819595, + -0.9800058800588007290937355024646722133204e-11), + C(0.4244549085628511778373438768121222815752, + 0.002935393851311701428647152230552122898291), + C(-0.1021340732357117208743299813648493928105, + -0.00001088377943049851799938998805451564893540), + C(-0.01000200120119126652710792390331206563616, + 0.1020612612857282306892368985525393707486e-7), + C(0.1000000000007333333333344266666666664457e-5, + 0.2000000000001333333333323199999999978819e-5), + C(0.1999999999994666666666675199999999990248e-5, 0), + C(0.3013403889237919660346644392864226952119, 0), + C(0.02503136792640367194699495234782353186858, 0), + C(0.002500031251171948248596912483183760683918, 0), + C(0, 0.004900078433419939164774792850907128053308), + C(0, -0.005100088434920074173454208832365950009419), + C(0, 0.2000000000005333333333341866666666676419e-5), + C(0, -48.16001211429122974789822893525016528191), + C(0, 0.4627407029504443513654142715903005954668e174), + C(0, -Inf), + C(0, 0), + C(-0, 0), + C(0, Inf), + C(0, -Inf), + C(NaN, NaN), + C(NaN, NaN), + C(NaN, NaN), + C(NaN, 0), + C(0, NaN), + C(NaN, NaN), + C(NaN, NaN), + C(0.01282473148489433743567240624939698290584, + -0.2105957276516618621447832572909153498104e-7), + C(0.01219875253423634378984109995893708152885, + -0.1813040560401824664088425926165834355953e-7), + C(0.1020408163265306334945473399689037886997e-7, + -0.1041232819658476285651490827866174985330e-25), + C(0.9803921568627452865036825956835185367356e-8, + -0.9227220299884665067601095648451913375754e-26), + C(0.5000000000000000002500000000000000003750e-9, + -0.1200000000000000001800000188712838420241e-29), + C(5.00000000000000000000025000000000000000000003e-12, + -1.20000000000000000000018000000000000000000004e-36), + C(5.00000000000000000000000002500000000000000000e-14, + -1.20000000000000000000000001800000000000000000e-42), + C(5e-301, 0)}; + TST(Dawson, 1e-20); + } + printf("#####################################\n"); + printf("SUCCESS (max relative error = %g)\n", errmax_all); +} + +#endif diff --git a/SKIRT/faddeeva/Faddeeva.hh b/SKIRT/faddeeva/Faddeeva.hh new file mode 100644 index 00000000..c4a2e971 --- /dev/null +++ b/SKIRT/faddeeva/Faddeeva.hh @@ -0,0 +1,62 @@ +/* Copyright (c) 2012 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* Available at: http://ab-initio.mit.edu/Faddeeva + + Header file for Faddeeva.cc; see that file for more information. */ + +#ifndef FADDEEVA_HH +#define FADDEEVA_HH 1 + +#include + +namespace Faddeeva { + +// compute w(z) = exp(-z^2) erfc(-iz) [ Faddeeva / scaled complex error func ] +extern std::complex w(std::complex z,double relerr=0); +extern double w_im(double x); // special-case code for Im[w(x)] of real x + +// Various functions that we can compute with the help of w(z) + +// compute erfcx(z) = exp(z^2) erfc(z) +extern std::complex erfcx(std::complex z, double relerr=0); +extern double erfcx(double x); // special case for real x + +// compute erf(z), the error function of complex arguments +extern std::complex erf(std::complex z, double relerr=0); +extern double erf(double x); // special case for real x + +// compute erfi(z) = -i erf(iz), the imaginary error function +extern std::complex erfi(std::complex z, double relerr=0); +extern double erfi(double x); // special case for real x + +// compute erfc(z) = 1 - erf(z), the complementary error function +extern std::complex erfc(std::complex z, double relerr=0); +extern double erfc(double x); // special case for real x + +// compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z) +extern std::complex Dawson(std::complex z, double relerr=0); +extern double Dawson(double x); // special case for real x + +} // namespace Faddeeva + +#endif // FADDEEVA_HH diff --git a/SKIRT/utils/AtomUtils.cpp b/SKIRT/utils/AtomUtils.cpp new file mode 100644 index 00000000..472b65d0 --- /dev/null +++ b/SKIRT/utils/AtomUtils.cpp @@ -0,0 +1,85 @@ +/*////////////////////////////////////////////////////////////////// +//// The SKIRT project -- advanced radiative transfer //// +//// © Astronomical Observatory, Ghent University //// +///////////////////////////////////////////////////////////////// */ + +#include "AtomUtils.hpp" +#include "Constants.hpp" +#include "FatalError.hpp" +#include "StringUtils.hpp" +#include +#include + +//////////////////////////////////////////////////////////////////// + +namespace +{ + static const std::map atomMap = { + {"H", 1}, {"He", 2}, {"Li", 3}, {"Be", 4}, {"B", 5}, {"C", 6}, {"N", 7}, {"O", 8}, + {"F", 9}, {"Ne", 10}, {"Na", 11}, {"Mg", 12}, {"Al", 13}, {"Si", 14}, {"P", 15}, {"S", 16}, + {"Cl", 17}, {"Ar", 18}, {"K", 19}, {"Ca", 20}, {"Sc", 21}, {"Ti", 22}, {"V", 23}, {"Cr", 24}, + {"Mn", 25}, {"Fe", 26}, {"Co", 27}, {"Ni", 28}, {"Cu", 29}, {"Zn", 30}}; + + static const vector masses = {1.0079, 4.0026, 6.941, 9.01218, 10.81, 12.011, 14.0067, 15.9994, + 18.9984, 20.179, 22.9898, 24.305, 26.9815, 28.0855, 30.9738, 32.06, + 35.453, 39.948, 39.0983, 40.08, 44.9559, 47.9, 50.9415, 51.996, + 54.938, 55.847, 58.9332, 58.7, 63.546, 65.38}; +} + +//////////////////////////////////////////////////////////////////// + +short AtomUtils::atomToZ(string element) +{ + return atomMap.at(element); +} + +//////////////////////////////////////////////////////////////////// + +double AtomUtils::mass(short Z) +{ + return masses[Z - 1] * Constants::amu(); +} + +//////////////////////////////////////////////////////////////////// + +std::pair AtomUtils::parseIon(string ion) +{ + // read ions + ion = StringUtils::squeeze(ion); + // split ion string into (atom, plus, ionization) + std::regex pattern("([A-Za-z]+)(\\+?)([0-9]*)"); + std::smatch match; + + if (!std::regex_match(ion, match, pattern)) throw FATALERROR("Could not parse ion format: " + ion); + + int Z, N; + string element = match[1].str(); + string plus = match[2].str(); + string number = match[3].str(); + + auto it = atomMap.find(element); + if (it == atomMap.end()) throw FATALERROR("No element found with name " + element + " for ion: " + ion); + Z = it->second; + + // enforce format to avoid ambiguities + if (number.empty()) + { + if (plus.empty()) // eg. Fe + N = Z; + else // eg. Fe+ + N = Z - 1; + } + else + { + if (!plus.empty()) // eg. Fe+14 + N = Z - StringUtils::toInt(number); + else // eg. Fe14 + throw FATALERROR("Ion format should contain a '+' eg. Fe+14 for ion: " + ion); + } + + if (N < 0 || N > Z) throw FATALERROR("Invalid ionization degree for ion: " + ion); + + return std::make_pair(Z, N); +} + +//////////////////////////////////////////////////////////////////// diff --git a/SKIRT/utils/AtomUtils.hpp b/SKIRT/utils/AtomUtils.hpp new file mode 100644 index 00000000..3a543bd7 --- /dev/null +++ b/SKIRT/utils/AtomUtils.hpp @@ -0,0 +1,33 @@ +/*////////////////////////////////////////////////////////////////// +//// The SKIRT project -- advanced radiative transfer //// +//// © Astronomical Observatory, Ghent University //// +///////////////////////////////////////////////////////////////// */ + +#ifndef ATOMUTILS_HPP +#define ATOMUTILS_HPP + +#include "Basics.hpp" + +/** This static class provides utility functions related to the treatment of atomic and ionic + species. */ +class AtomUtils final +{ +public: + /** This function returns the atomic number of the specified elements. */ + static short atomToZ(string element); + + /** This function returns the mass of the specified atomic number in SI units. */ + static double mass(short Z); + + /** This function returns a pair of the atomic number and number of electrons (Z,N) of the + specified ion string. The ion string can have the following formats, as listed in the table + below. Any other format will throw an error. + + + +
FormatReturns
'Z'\f$(Z,N)=(Z,Z)\f$
'Z+'\f$(Z,N)=(Z,Z-1)\f$
'Z+I'\f$(Z,N)=(Z,Z-I)\f$
The number of + electrons, N, can range from 0 to Z, any other value will result in an error. */ + static std::pair parseIon(string ion); +}; + +#endif