Merge pull request #455 from vsnever/enhancement/add_cx_to_total_radi…

…ated_power

Add emission due to CX with neutral hydrogen to the TotalRadiatedPower model

CHANGELOG.md

@@ -16,6 +16,7 @@ New:
 * **Beam dispersion calculation has changed from sigma(z) = sigma + z * tan(alpha) to sigma(z) = sqrt(sigma^2 + (z * tan(alpha))^2) for consistancy with the Gaussian beam model. Attention!!! The results of BES and CX spectroscopy are affected by this change. (#414)**
 * Improved beam direction calculation to allow for natural broadening of the BES line shape due to beam divergence. (#414)
 * Add kwargs to invert_regularised_nnls to pass them to scipy.optimize.nnls. (#438)
+* Add the power radiated in spectral lines due to charge exchange with thermal neutral hydrogen to the TotalRadiatedPower model. (#370)
 * Add thermal charge-exchange emission model. (#57)
 * PECs for C VI spectral lines for n <= 5 are now included in populate(). Rerun populate() after upgrading to 1.5 to update the atomic data repository.
 * All interpolated atomic rates now return 0 if plasma parameters <= 0, which matches the behaviour of emission models. (#450)

cherab/core/atomic/rates.pxd

@@ -86,7 +86,7 @@ cdef class _RadiatedPower:
         readonly Element element
         readonly int charge
 
-    cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999
+    cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999
 
 
 cdef class LineRadiationPower(_RadiatedPower):

cherab/core/atomic/rates.pyx

@@ -284,7 +284,7 @@ cdef class _RadiatedPower:
         """
         return self.evaluate(electron_density, electron_temperature)
 
-    cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
+    cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
         """
         Evaluate the radiated power at the given plasma conditions.
 

cherab/core/model/plasma/total_radiated_power.pxd

@@ -18,7 +18,7 @@
 
 
 from cherab.core.atomic.elements cimport Element
-from cherab.core.atomic.rates cimport LineRadiationPower, ContinuumPower
+from cherab.core.atomic.rates cimport LineRadiationPower, ContinuumPower, CXRadiationPower
 from cherab.core.plasma cimport PlasmaModel
 from cherab.core.species cimport Species
 
@@ -30,7 +30,9 @@ cdef class TotalRadiatedPower(PlasmaModel):
         Element _element
         int _charge
         Species _line_rad_species, _recom_species
+        list _hydrogen_species
         LineRadiationPower _plt_rate
         ContinuumPower _prb_rate
+        CXRadiationPower _prc_rate
 
     cdef int _populate_cache(self) except -1

cherab/core/model/plasma/total_radiated_power.pyx

@@ -20,9 +20,40 @@
 from raysect.optical cimport Spectrum, Point3D, Vector3D
 from cherab.core cimport Plasma, AtomicData
 from cherab.core.utility.constants cimport RECIP_4_PI
+from cherab.core.atomic.elements import hydrogen, deuterium, tritium
 
 
 cdef class TotalRadiatedPower(PlasmaModel):
+    r"""
+    Emitter that calculates total power radiated by a given ion, which includes:
+
+    - line power due to electron impact excitation,
+    - continuum and line power due to recombination and Bremsstrahlung,
+    - line power due to charge exchange with thermal neutral hydrogen and its isotopes.
+
+    The emission calculated by this model is spectrally unresolved,
+    which means that the total radiated power will be spread of the entire
+    observable spectral range.
+
+    .. math::
+        \epsilon_{\mathrm{total}} = \frac{1}{4 \pi \Delta\lambda} \left(
+        n_{Z_\mathrm{i}} n_\mathrm{e} C_{\mathrm{excit}}(n_\mathrm{e}, T_\mathrm{e}) + 
+        n_{Z_\mathrm{i} + 1} n_\mathrm{e} C_{\mathrm{recomb}}(n_\mathrm{e}, T_\mathrm{e}) +
+        n_{Z_\mathrm{i} + 1} n_\mathrm{hyd} C_{\mathrm{cx}}(n_\mathrm{e}, T_\mathrm{e}) \right)
+
+    where :math:`n_{Z_\mathrm{i}}` is the target species density;
+    :math:`n_{Z_\mathrm{i} + 1}` is the recombining species density;
+    :math:`n_{\mathrm{hyd}}` is the total density of all hydrogen isotopes;
+    :math:`C_{\mathrm{excit}}, C_{\mathrm{recomb}}, C_{\mathrm{cx}}` are the radiated power
+    coefficients in :math:`W m^3` due to electron impact excitation, recombination
+    + Bremsstrahlung and charge exchange with thermal neutral hydrogen, respectively;
+    :math:`\Delta\lambda` is the observable spectral range.
+
+    :param Element element: The atomic element/isotope.
+    :param int charge: The charge state of the element/isotope.
+    :param Plasma plasma: The plasma to which this emission model is attached. Default is None.
+    :param AtomicData atomic_data: The atomic data provider for this model. Default is None.
+    """
 
     def __init__(self, Element element, int charge, Plasma plasma=None, AtomicData atomic_data=None):
 
@@ -42,7 +73,9 @@ cdef class TotalRadiatedPower(PlasmaModel):
 
         cdef:
             int i
-            double ne, ni, ni_upper, te, plt_radiance, prb_radiance
+            double ne, ni, ni_upper, nhyd, te
+            double power_density, radiance
+            Species hyd_species
 
         # cache data on first run
         if not self._cache_loaded:
@@ -60,22 +93,35 @@ cdef class TotalRadiatedPower(PlasmaModel):
 
         ni_upper = self._recom_species.distribution.density(point.x, point.y, point.z)
 
+        nhyd = 0
+        for hyd_species in self._hydrogen_species:
+            nhyd += hyd_species.distribution.density(point.x, point.y, point.z)
+
         # add emission to spectrum
-        if self._plt_rate and ni > 0:
-            plt_radiance = RECIP_4_PI * self._plt_rate.evaluate(ne, te) * ne * ni / (spectrum.max_wavelength - spectrum.min_wavelength)
-        else:
-            plt_radiance = 0
-        if self._prb_rate and ni_upper > 0:
-            prb_radiance = RECIP_4_PI * self._prb_rate.evaluate(ne, te) * ne * ni_upper / (spectrum.max_wavelength - spectrum.min_wavelength)
-        else:
-            prb_radiance = 0
+        power_density = 0
+
+        if self._plt_rate and ni > 0:  # excitation
+            power_density += self._plt_rate.evaluate(ne, te) * ne * ni
+
+        if self._prb_rate and ni_upper > 0:  # recombination + bremsstrahlung
+            power_density += self._prb_rate.evaluate(ne, te) * ne * ni_upper
+
+        if self._prc_rate and ni_upper > 0 and nhyd > 0:  # charge exchange
+            power_density += self._prc_rate.evaluate(ne, te) * nhyd * ni_upper
+
+        radiance = RECIP_4_PI * power_density / (spectrum.max_wavelength - spectrum.min_wavelength)
+
         for i in range(spectrum.bins):
-            spectrum.samples_mv[i] += plt_radiance + prb_radiance
+            spectrum.samples_mv[i] += radiance
 
         return spectrum
 
     cdef int _populate_cache(self) except -1:
 
+        cdef:
+            Species hyd_species
+            Element hyd_isotope
+
         # sanity checks
         if self._plasma is None:
             raise RuntimeError("The emission model is not connected to a plasma object.")
@@ -85,7 +131,7 @@ cdef class TotalRadiatedPower(PlasmaModel):
         # cache line radiation species and rate
         self._plt_rate = self._atomic_data.line_radiated_power_rate(self._element, self._charge)
         try:
-            self._line_rad_species = self._plasma.composition.get(self._element, self._charge)
+            self._line_rad_species = self._plasma.get_composition().get(self._element, self._charge)
         except ValueError:
             raise RuntimeError("The plasma object does not contain the required ion species for calculating"
                                "total line radiaton, (element={}, ionisation={})."
@@ -94,12 +140,24 @@ cdef class TotalRadiatedPower(PlasmaModel):
         # cache recombination species and radiation rate
         self._prb_rate = self._atomic_data.continuum_radiated_power_rate(self._element, self._charge+1)
         try:
-            self._recom_species = self._plasma.composition.get(self._element, self._charge+1)
+            self._recom_species = self._plasma.get_composition().get(self._element, self._charge+1)
         except ValueError:
             raise RuntimeError("The plasma object does not contain the required ion species for calculating"
                                "recombination/continuum emission, (element={}, ionisation={})."
                                "".format(self._element.symbol, self._charge+1))
 
+        # cache hydrogen species and CX radiation rate
+        self._prc_rate = self._atomic_data.cx_radiated_power_rate(self._element, self._charge+1)
+
+        self._hydrogen_species = []
+        for hyd_isotope in (hydrogen, deuterium, tritium):
+            try:
+                hyd_species = self._plasma.get_composition().get(hyd_isotope, 0)
+            except ValueError:
+                pass
+            else:
+                self._hydrogen_species.append(hyd_species)
+
         self._cache_loaded = True
 
     def _change(self):
@@ -108,5 +166,7 @@ cdef class TotalRadiatedPower(PlasmaModel):
         self._cache_loaded = False
         self._line_rad_species = None
         self._recom_species = None
+        self._hydrogen_species = None
         self._plt_rate = None
         self._prb_rate = None
+        self._prc_rate = None

cherab/core/tests/test_total_radiated_power.py

@@ -0,0 +1,142 @@
+# Copyright 2016-2023 Euratom
+# Copyright 2016-2023 United Kingdom Atomic Energy Authority
+# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas
+#
+# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the
+# European Commission - subsequent versions of the EUPL (the "Licence");
+# You may not use this work except in compliance with the Licence.
+# You may obtain a copy of the Licence at:
+#
+# https://joinup.ec.europa.eu/software/page/eupl5
+#
+# Unless required by applicable law or agreed to in writing, software distributed
+# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR
+# CONDITIONS OF ANY KIND, either express or implied.
+#
+# See the Licence for the specific language governing permissions and limitations
+# under the Licence.
+
+import unittest
+
+import numpy as np
+
+from raysect.core import Point3D, Vector3D
+from raysect.optical import Ray, World
+
+from cherab.core.atomic import AtomicData, LineRadiationPower, ContinuumPower, CXRadiationPower
+from cherab.core.atomic import deuterium, hydrogen, nitrogen
+from cherab.tools.plasmas.slab import build_constant_slab_plasma
+from cherab.core.model import TotalRadiatedPower
+
+from cherab.core.utility import EvAmuToMS, EvToJ
+
+
+class ConstantLineRadiationPower(LineRadiationPower):
+    """
+    Constant line radiation power coefficient.
+    """
+
+    def __init__(self, value):
+        self.value = value
+
+    def evaluate(self, density, temperature):
+
+        return self.value
+
+
+class ConstantContinuumPower(ContinuumPower):
+    """
+    Constant continuum power coefficient.
+    """
+
+    def __init__(self, value):
+        self.value = value
+
+    def evaluate(self, density, temperature):
+
+        return self.value
+
+
+class ConstantCXRadiationPower(CXRadiationPower):
+    """
+    Constant charge exchange radiation power coefficient.
+    """
+
+    def __init__(self, value):
+        self.value = value
+
+    def evaluate(self, density, temperature):
+
+        return self.value
+
+
+class MockAtomicData(AtomicData):
+    """Fake atomic data for test purpose."""
+
+    def line_radiated_power_rate(self, element, charge):
+
+        return ConstantLineRadiationPower(1.e-32)
+
+    def continuum_radiated_power_rate(self, element, charge):
+
+        return ConstantContinuumPower(1.e-33)
+
+    def cx_radiated_power_rate(self, element, charge):
+
+        return ConstantCXRadiationPower(1.e-31)
+
+
+class TestTotalRadiatedPower(unittest.TestCase):
+
+    def setUp(self):
+
+        self.world = World()
+
+        self.atomic_data = MockAtomicData()
+
+        plasma_species = [(deuterium, 0, 1.e18, 500., Vector3D(0, 0, 0)),
+                          (hydrogen, 0, 1.e18, 500., Vector3D(0, 0, 0)),
+                          (nitrogen, 6, 5.e18, 1100., Vector3D(0, 0, 0)),
+                          (nitrogen, 7, 1.e19, 1100., Vector3D(0, 0, 0))]
+        self.plasma = build_constant_slab_plasma(length=1.2, width=1.2, height=1.2,
+                                                 electron_density=1e19, electron_temperature=1000.,
+                                                 plasma_species=plasma_species)
+        self.plasma.parent = self.world
+        self.plasma.atomic_data = self.atomic_data
+
+    def test_beam_density(self):
+
+        self.plasma.models = [TotalRadiatedPower(nitrogen, 6)]
+
+        # observing
+        origin = Point3D(1.5, 0, 0)
+        direction = Vector3D(-1, 0, 0)
+        ray = Ray(origin=origin, direction=direction,
+                  min_wavelength=500., max_wavelength=550., bins=2)
+        radiated_power = ray.trace(self.world).total()
+
+        # validating
+        ne = self.plasma.electron_distribution.density(0.5, 0.5, 0.5)
+        n_n6 = self.plasma.composition[(nitrogen, 6)].distribution.density(0.5, 0.5, 0.5)
+        n_n7 = self.plasma.composition[(nitrogen, 7)].distribution.density(0.5, 0.5, 0.5)
+        n_h0 = self.plasma.composition[(hydrogen, 0)].distribution.density(0.5, 0.5, 0.5)
+        n_d0 = self.plasma.composition[(deuterium, 0)].distribution.density(0.5, 0.5, 0.5)
+
+        integration_length = 1.2
+
+        plt_rate = self.atomic_data.line_radiated_power_rate(nitrogen, 6).value
+        plt_radiance = 0.25 / np.pi * plt_rate * ne * n_n6 * integration_length
+
+        prb_rate = self.atomic_data.continuum_radiated_power_rate(nitrogen, 7).value
+        prb_radiance = 0.25 / np.pi * prb_rate * ne * n_n7 * integration_length
+
+        prc_rate = self.atomic_data.cx_radiated_power_rate(nitrogen, 7).value
+        prc_radiance = 0.25 / np.pi * prc_rate * (n_h0 + n_d0) * n_n7 * integration_length
+
+        test_radiated_power = plt_radiance + prb_radiance + prc_radiance
+
+        self.assertAlmostEqual(radiated_power / test_radiated_power, 1., delta=1e-8)
+
+
+if __name__ == '__main__':
+    unittest.main()

cherab/openadas/rates/radiated_power.pyx

@@ -65,7 +65,7 @@ cdef class LineRadiationPower(CoreLineRadiationPower):
         extrapolation_type = 'nearest' if extrapolate else 'none'
         self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY)
 
-    cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
+    cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
         if electron_density <= 0 or electron_temperature <= 0:
@@ -81,7 +81,7 @@ cdef class NullLineRadiationPower(CoreLineRadiationPower):
     Needed for use cases where the required atomic data is missing.
     """
 
-    cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
+    cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
         return 0.0
 
 
@@ -127,7 +127,7 @@ cdef class ContinuumPower(CoreContinuumPower):
         extrapolation_type = 'nearest' if extrapolate else 'none'
         self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY)
 
-    cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
+    cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
         if electron_density <= 0 or electron_temperature <= 0:
@@ -143,7 +143,7 @@ cdef class NullContinuumPower(CoreContinuumPower):
     Needed for use cases where the required atomic data is missing.
     """
 
-    cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
+    cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
         return 0.0
 
 
@@ -188,7 +188,7 @@ cdef class CXRadiationPower(CoreCXRadiationPower):
         extrapolation_type = 'linear' if extrapolate else 'none'
         self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY)
 
-    cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
+    cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
 
         # need to handle zeros, also density and temperature can become negative due to cubic interpolation
         if electron_density <= 0 or electron_temperature <= 0:
@@ -204,5 +204,5 @@ cdef class NullCXRadiationPower(CoreCXRadiationPower):
     Needed for use cases where the required atomic data is missing.
     """
 
-    cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
+    cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999:
         return 0.0

docs/source/models/radiated_power/radiated_power.rst

@@ -2,4 +2,4 @@
 Total Radiated Power
 ====================
 
-Documentation for this model will go here soon...
+.. autoclass:: cherab.core.model.plasma.total_radiated_power.TotalRadiatedPower