Add ThermalCXLine model

CHANGELOG.md

@@ -16,6 +16,8 @@ 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 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)
 
 Bug fixes:

cherab/core/atomic/interface.pxd

@@ -45,6 +45,8 @@ cdef class AtomicData:
 
     cpdef RecombinationPEC recombination_pec(self, Element ion, int charge, tuple transition)
 
+    cpdef ThermalCXPEC thermal_cx_pec(self, Element donor_ion, int donor_charge, Element receiver_ion, int receiver_charge, tuple transition)
+
     cpdef TotalRadiatedPower total_radiated_power(self, Element element)
 
     cpdef LineRadiationPower line_radiated_power_rate(self, Element element, int charge)

cherab/core/atomic/interface.pyx

@@ -98,6 +98,9 @@ cdef class AtomicData:
 
         raise NotImplementedError("The recombination() virtual method is not implemented for this atomic data source.")
 
+    cpdef ThermalCXPEC thermal_cx_pec(self, Element donor_ion, int donor_charge, Element receiver_ion, int receiver_charge, tuple transition):
+        raise NotImplementedError("The thermal_cx_pec() virtual method is not implemented for this atomic data source.")
+
     cpdef TotalRadiatedPower total_radiated_power(self, Element element):
         """
         The total (summed over all charge states) radiated power

cherab/core/atomic/rates.pxd

@@ -46,8 +46,8 @@ cdef class RecombinationPEC(_PECRate):
     pass
 
 
-cdef class ThermalCXPEC(_PECRate):
-    pass
+cdef class ThermalCXPEC:
+    cpdef double evaluate(self, double electron_density, double electron_temperature, double donor_temperature) except? -1e999
 
 
 cdef class BeamCXPEC:

cherab/core/atomic/rates.pyx

@@ -130,11 +130,27 @@ cdef class RecombinationPEC(_PECRate):
     pass
 
 
-cdef class ThermalCXPEC(_PECRate):
+cdef class ThermalCXPEC:
     """
     Thermal charge exchange rate coefficient.
     """
-    pass
+
+    def __call__(self, double electron_density, double electron_temperature, donor_temperature):
+        """Returns a CX photon emissivity coefficient at the specified plasma conditions.
+
+        This function just wraps the cython evaluate() method.
+        """
+        return self.evaluate(electron_density, electron_temperature, donor_temperature)
+
+    cpdef double evaluate(self, double electron_density, double electron_temperature, double donor_temperature) except? -1e999:
+        """Returns a CX photon emissivity coefficient at given conditions.
+
+        :param electron_density: Electron density in m^-3.
+        :param electron_temperature: Electron temperature in eV.
+        :param donor_temperature: Donor temperature in eV.
+        :return: The effective CX PEC rate in W/m^3.
+        """
+        raise NotImplementedError("The evaluate() virtual method must be implemented.")
 
 
 cdef class BeamCXPEC:

cherab/core/model/plasma/__init__.pxd

@@ -20,4 +20,5 @@
 from cherab.core.model.plasma.bremsstrahlung cimport Bremsstrahlung
 from cherab.core.model.plasma.impact_excitation cimport ExcitationLine
 from cherab.core.model.plasma.recombination cimport RecombinationLine
+from cherab.core.model.plasma.thermal_cx cimport ThermalCXLine
 from cherab.core.model.plasma.total_radiated_power cimport TotalRadiatedPower

cherab/core/model/plasma/__init__.py

@@ -20,4 +20,5 @@
 from .bremsstrahlung import Bremsstrahlung
 from .impact_excitation import ExcitationLine
 from .recombination import RecombinationLine
+from .thermal_cx import ThermalCXLine
 from .total_radiated_power import TotalRadiatedPower

cherab/core/model/plasma/thermal_cx.pxd

@@ -0,0 +1,35 @@
+# Copyright 2016-2018 Euratom
+# Copyright 2016-2018 United Kingdom Atomic Energy Authority
+# Copyright 2016-2018 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.
+
+from cherab.core.atomic cimport Line
+from cherab.core.plasma cimport PlasmaModel
+from cherab.core.species cimport Species
+from cherab.core.model.lineshape cimport LineShapeModel
+
+
+cdef class ThermalCXLine(PlasmaModel):
+
+    cdef:
+        Line _line
+        double _wavelength
+        Species _target_species
+        list _rates
+        LineShapeModel _lineshape
+        object _lineshape_class, _lineshape_args, _lineshape_kwargs
+
+    cdef int _populate_cache(self) except -1

cherab/core/model/plasma/thermal_cx.pyx

@@ -0,0 +1,163 @@
+# Copyright 2016-2018 Euratom
+# Copyright 2016-2018 United Kingdom Atomic Energy Authority
+# Copyright 2016-2018 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.
+
+from raysect.optical cimport Spectrum, Point3D, Vector3D
+from cherab.core cimport Plasma, AtomicData
+from cherab.core.atomic cimport ThermalCXPEC
+from cherab.core.model.lineshape cimport GaussianLine, LineShapeModel
+from cherab.core.utility.constants cimport RECIP_4_PI
+
+
+cdef class ThermalCXLine(PlasmaModel):
+    r"""
+    Emitter that calculates spectral line emission from a plasma object
+    as a result of thermal charge exchange of the target species with the donor species.
+
+    .. math::
+        \epsilon_{\mathrm{CX}}(\lambda) = \frac{1}{4 \pi} n_{Z_\mathrm{i} + 1}
+        \sum_j{n_{Z_\mathrm{j}} \mathrm{PEC}_{\mathrm{cx}}(n_\mathrm{e}, T_\mathrm{e}, T_{Z_\mathrm{j}})}
+        f(\lambda),
+
+    where :math:`n_{Z_\mathrm{i} + 1}` is the receiver species density,
+    :math:`n_{Z_\mathrm{j}}` is the donor species density,
+    :math:`\mathrm{PEC}_{\mathrm{cx}}` is the thermal CX photon emission coefficient
+    for the specified spectral line of the :math:`Z_\mathrm{i}` ion,
+    :math:`T_{Z_\mathrm{j}}` is the donor species temperature,
+    :math:`f(\lambda)` is the normalised spectral line shape,
+
+    :param Line line: Spectroscopic emission line object.
+    :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.
+    :param object lineshape: Line shape model class. Default is None (GaussianLine).
+    :param object lineshape_args: A list of line shape model arguments. Default is None.
+    :param object lineshape_kwargs: A dictionary of line shape model keyword arguments. Default is None.
+
+    :ivar Plasma plasma: The plasma to which this emission model is attached.
+    :ivar AtomicData atomic_data: The atomic data provider for this model.
+    """
+
+    def __init__(self, Line line, Plasma plasma=None, AtomicData atomic_data=None, object lineshape=None,
+                 object lineshape_args=None, object lineshape_kwargs=None):
+
+        super().__init__(plasma, atomic_data)
+
+        self._line = line
+
+        self._lineshape_class = lineshape or GaussianLine
+        if not issubclass(self._lineshape_class, LineShapeModel):
+            raise TypeError("The attribute lineshape must be a subclass of LineShapeModel.")
+
+        if lineshape_args:
+            self._lineshape_args = lineshape_args
+        else:
+            self._lineshape_args = []
+        if lineshape_kwargs:
+            self._lineshape_kwargs = lineshape_kwargs
+        else:
+            self._lineshape_kwargs = {}
+
+        # ensure that cache is initialised
+        self._change()
+
+    def __repr__(self):
+        return '<ThermalCXLine: element={}, charge={}, transition={}>'.format(self._line.element.name, self._line.charge, self._line.transition)
+
+    cpdef Spectrum emission(self, Point3D point, Vector3D direction, Spectrum spectrum):
+
+        cdef:
+            double ne, te, receiver_density, donor_density, donor_temperature, weighted_rate, radiance
+            Species species
+            ThermalCXPEC rate
+
+        # cache data on first run
+        if self._target_species is None:
+            self._populate_cache()
+
+        ne = self._plasma.get_electron_distribution().density(point.x, point.y, point.z)
+        if ne <= 0.0:
+            return spectrum
+
+        te = self._plasma.get_electron_distribution().effective_temperature(point.x, point.y, point.z)
+        if te <= 0.0:
+            return spectrum
+
+        receiver_density = self._target_species.distribution.density(point.x, point.y, point.z)
+        if receiver_density <= 0.0:
+            return spectrum
+
+        # obtain composite CX PEC by iterating over all possible CX donors
+        weighted_rate = 0
+        for species, rate in self._rates:
+            donor_density = species.distribution.density(point.x, point.y, point.z)
+            donor_temperature = species.distribution.effective_temperature(point.x, point.y, point.z)
+            weighted_rate += donor_density * rate.evaluate(ne, te, donor_temperature)
+
+        # add emission line to spectrum
+        radiance = RECIP_4_PI * weighted_rate * receiver_density
+        return self._lineshape.add_line(radiance, point, direction, spectrum)
+
+    cdef int _populate_cache(self) except -1:
+
+        cdef:
+            int receiver_charge
+            Species species
+            ThermalCXPEC rate
+
+        # sanity checks
+        if self._plasma is None:
+            raise RuntimeError("The emission model is not connected to a plasma object.")
+        if self._atomic_data is None:
+            raise RuntimeError("The emission model is not connected to an atomic data source.")
+
+        if self._line is None:
+            raise RuntimeError("The emission line has not been set.")
+
+        # locate target species
+        receiver_charge = self._line.charge + 1
+        try:
+            self._target_species = self._plasma.composition.get(self._line.element, receiver_charge)
+        except ValueError:
+            raise RuntimeError("The plasma object does not contain the ion species for the specified CX line "
+                               "(element={}, ionisation={}).".format(self._line.element.symbol, receiver_charge))
+
+        # obtain rate functions
+        self._rates = []
+        # iterate over all posible electron donors in plasma composition
+        # and for each donor, cache the PEC rate function for the CX reaction with this receiver
+        for species in self._plasma.composition:
+            # exclude the receiver species from the list of donors and omit fully ionised species
+            if species != self._target_species and species.charge < species.element.atomic_number:
+                rate = self._atomic_data.thermal_cx_pec(species.element, species.charge,  # donor
+                                                        self._line.element, receiver_charge,  # receiver
+                                                        self._line.transition)
+                self._rates.append((species, rate))
+
+        # identify wavelength
+        self._wavelength = self._atomic_data.wavelength(self._line.element, self._line.charge, self._line.transition)
+
+        # instance line shape renderer
+        self._lineshape = self._lineshape_class(self._line, self._wavelength, self._target_species, self._plasma,
+                                                *self._lineshape_args, **self._lineshape_kwargs)
+
+    def _change(self):
+
+        # clear cache to force regeneration on first use
+        self._target_species = None
+        self._wavelength = 0.0
+        self._rates = None
+        self._lineshape = None