Merge branch 'development' into feature/thermal_cx_line

CHANGELOG.md

@@ -18,10 +18,13 @@ New:
 * 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:
 * Fix deprecated transforms being cached in LaserMaterial after laser.transform update (#420)
 * Fix IRVB calculate sensitivity method.
+* Fix missing donor_metastable attribute in the core BeamCXPEC class (#411).
+* **Fix the receiver ion density being passed to the BeamCXPEC instead of the total ion density in the BeamCXLine. Also fix incorrect BeamCXPEC dosctrings. Attention!!! The results of CX spectroscopy are affected by this change. (#441)**
 
 Release 1.4.0 (3 Feb 2023)
 -------------------

cherab/core/atomic/interface.pyx

@@ -1,6 +1,6 @@
-# Copyright 2016-2022 Euratom
-# Copyright 2016-2022 United Kingdom Atomic Energy Authority
-# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas
+# Copyright 2016-2024 Euratom
+# Copyright 2016-2024 United Kingdom Atomic Energy Authority
+# Copyright 2016-2024 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");
@@ -29,73 +29,120 @@ cdef class AtomicData:
 
     cpdef double wavelength(self, Element ion, int charge, tuple transition):
         """
-        Returns the natural wavelength of the specified transition in nm.
+        The natural wavelength of the specified transition in nm.
         """
 
         raise NotImplementedError("The wavelength() virtual method is not implemented for this atomic data source.")
 
     cpdef IonisationRate ionisation_rate(self, Element ion, int charge):
+        """
+        Electron impact ionisation rate for a given species in m^3/s.
+        """
+
         raise NotImplementedError("The ionisation_rate() virtual method is not implemented for this atomic data source.")
 
     cpdef RecombinationRate recombination_rate(self, Element ion, int charge):
+        """
+        Recombination rate for a given species in m^3/s.
+        """
+
         raise NotImplementedError("The recombination_rate() virtual method is not implemented for this atomic data source.")
 
     cpdef ThermalCXRate thermal_cx_rate(self, Element donor_ion, int donor_charge, Element receiver_ion, int receiver_charge):
+        """
+        Thermal charge exchange effective rate coefficient for a given donor and receiver species in m^3/s.
+        """
+
         raise NotImplementedError("The thermal_cx_rate() virtual method is not implemented for this atomic data source.")
 
     cpdef list beam_cx_pec(self, Element donor_ion, Element receiver_ion, int receiver_charge, tuple transition):
         """
-        Returns a list of applicable charge exchange emission rates in W.m^3.
+        A list of Effective charge exchange photon emission coefficient for a given donor (beam) in W.m^3.
         """
 
         raise NotImplementedError("The cxs_rates() virtual method is not implemented for this atomic data source.")
 
     cpdef BeamStoppingRate beam_stopping_rate(self, Element beam_ion, Element plasma_ion, int charge):
         """
-        Returns a list of applicable beam stopping coefficients in m^3/s.
+        Beam stopping coefficient for a given beam and target species in m^3/s.
         """
 
         raise NotImplementedError("The beam_stopping() virtual method is not implemented for this atomic data source.")
 
     cpdef BeamPopulationRate beam_population_rate(self, Element beam_ion, int metastable, Element plasma_ion, int charge):
         """
-        Returns a list of applicable dimensionless beam population coefficients.
+        Dimensionless Beam population coefficient for a given beam and target species.
         """
 
         raise NotImplementedError("The beam_population() virtual method is not implemented for this atomic data source.")
 
     cpdef BeamEmissionPEC beam_emission_pec(self, Element beam_ion, Element plasma_ion, int charge, tuple transition):
         """
-        Returns a list of applicable beam emission coefficients in W.m^3.
+        The beam photon emission coefficient for a given beam and target species
+        and a given transition in W.m^3.
         """
 
         raise NotImplementedError("The beam_emission() virtual method is not implemented for this atomic data source.")
 
     cpdef ImpactExcitationPEC impact_excitation_pec(self, Element ion, int charge, tuple transition):
+        """
+        Electron impact excitation photon emission coefficient for a given species in W.m^3.
+        """
+
         raise NotImplementedError("The impact_excitation() virtual method is not implemented for this atomic data source.")
 
     cpdef RecombinationPEC recombination_pec(self, Element ion, int charge, tuple transition):
+        """
+        Recombination photon emission coefficient for a given species in W.m^3.
+        """
+
         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
+        in equilibrium conditions for a given species in W.m^3.
+        """
+
         raise NotImplementedError("The total_radiated_power() virtual method is not implemented for this atomic data source.")
 
     cpdef LineRadiationPower line_radiated_power_rate(self, Element element, int charge):
+        """
+        Line radiated power coefficient for a given species in W.m^3.
+        """
+
         raise NotImplementedError("The line_radiated_power_rate() virtual method is not implemented for this atomic data source.")
 
     cpdef ContinuumPower continuum_radiated_power_rate(self, Element element, int charge):
+        """
+        Continuum radiated power coefficient for a given species in W.m^3.
+        """
+
         raise NotImplementedError("The continuum_radiated_power_rate() virtual method is not implemented for this atomic data source.")
 
     cpdef CXRadiationPower cx_radiated_power_rate(self, Element element, int charge):
+        """
+        Charge exchange radiated power coefficient for a given species in W.m^3.
+        """
+
         raise NotImplementedError("The cx_radiated_power_rate() virtual method is not implemented for this atomic data source.")
 
     cpdef FractionalAbundance fractional_abundance(self, Element ion, int charge):
+        """
+        Fractional abundance of a given species in thermodynamic equilibrium.
+        """
+
         raise NotImplementedError("The fractional_abundance() virtual method is not implemented for this atomic data source.")
 
     cpdef ZeemanStructure zeeman_structure(self, Line line, object b_field=None):
+        r"""
+        Wavelengths and ratios of :math:`\pi`-/:math:`\sigma`-polarised Zeeman components
+        for any given value of magnetic field strength.
+        """
+
         raise NotImplementedError("The zeeman_structure() virtual method is not implemented for this atomic data source.")
 
     cpdef FreeFreeGauntFactor free_free_gaunt_factor(self):

cherab/core/atomic/rates.pxd

@@ -51,6 +51,8 @@ cdef class ThermalCXPEC:
 
 
 cdef class BeamCXPEC:
+    cdef readonly int donor_metastable
+
     cpdef double evaluate(self, double energy, double temperature, double density, double z_effective, double b_field) except? -1e999
 
 

cherab/core/atomic/rates.pyx

@@ -160,8 +160,13 @@ cdef class BeamCXPEC:
     transition :math:`n\rightarrow n'` of ion :math:`Z^{(\alpha+1)+}` with electron donor
     :math:`H^0` in metastable state :math:`m_{i}`. Equivalent to
     :math:`q^{eff}_{n\rightarrow n'}` in `adf12 <http://open.adas.ac.uk/adf12>_`.
+
+    :param donor_metastable: The metastable state of the donor species for which the rate data applies.
     """
 
+    def __init__(self, int donor_metastable):
+        self.donor_metastable = donor_metastable
+
     def __call__(self, double energy, double temperature, double density, double z_effective, double b_field):
         """Evaluates the Beam CX rate at the given plasma conditions.
 
@@ -174,7 +179,7 @@ cdef class BeamCXPEC:
 
         :param float energy: Interaction energy in eV/amu.
         :param float temperature: Receiver ion temperature in eV.
-        :param float density: Receiver ion density in m^-3
+        :param float density: Plasma total ion density in m^-3
         :param float z_effective: Plasma Z-effective.
         :param float b_field: Magnetic field magnitude in Tesla.
         :return: The effective rate

cherab/core/model/beam/charge_exchange.pxd

@@ -36,7 +36,7 @@ cdef class BeamCXLine(BeamModel):
         object _lineshape_class, _lineshape_args, _lineshape_kwargs
 
     cdef double _composite_cx_rate(self, double x, double y, double z, double interaction_energy,
-                                          Vector3D donor_velocity, double receiver_temperature, double receiver_density) except? -1e999
+                                          Vector3D donor_velocity, double receiver_temperature) except? -1e999
 
     cdef double _beam_population(self, double x, double y, double z, Vector3D beam_velocity, list population_data) except? -1e999
 

cherab/core/model/beam/charge_exchange.pyx

@@ -159,7 +159,7 @@ cdef class BeamCXLine(BeamModel):
         interaction_energy = ms_to_evamu(interaction_speed)
 
         # calculate the composite charge-exchange emission coefficient
-        emission_rate = self._composite_cx_rate(x, y, z, interaction_energy, donor_velocity, receiver_temperature, receiver_density)
+        emission_rate = self._composite_cx_rate(x, y, z, interaction_energy, donor_velocity, receiver_temperature)
 
         # spectral line emission in W/m^3/str
         radiance = RECIP_4_PI * donor_density * receiver_density * emission_rate
@@ -170,7 +170,7 @@ cdef class BeamCXLine(BeamModel):
     @cython.wraparound(False)
     @cython.cdivision(True)
     cdef double _composite_cx_rate(self, double x, double y, double z, double interaction_energy,
-                                   Vector3D donor_velocity, double receiver_temperature, double receiver_density) except? -1e999:
+                                   Vector3D donor_velocity, double receiver_temperature) except? -1e999:
         """
         Performs a beam population weighted average of the effective cx rates.
 
@@ -189,23 +189,26 @@ cdef class BeamCXLine(BeamModel):
         :param interaction_energy: The donor-receiver interaction energy in eV/amu.
         :param donor_velocity: A Vector defining the donor particle velocity in m/s.
         :param receiver_temperature: The receiver species temperature in eV.
-        :param receiver_density: The receiver species density in m^-3
         :return: The composite charge exchange rate in W.m^3.
         """
 
         cdef:
-            double z_effective, b_field, rate, total_population, population, effective_rate
+            double ion_density, z_effective
+            double b_field
+            double rate, effective_rate
+            double population, total_population
             BeamCXPEC cx_rate
             list population_data
 
-        # calculate z_effective and the B-field magnitude
+        # calculate ion density, z_effective and the B-field magnitude
+        ion_density = self._plasma.ion_density(x, y, z)
         z_effective = self._plasma.z_effective(x, y, z)
         b_field = self._plasma.get_b_field().evaluate(x, y, z).get_length()
 
         # rate for the ground state (metastable = 1)
         rate = self._ground_beam_rate.evaluate(interaction_energy,
                                                receiver_temperature,
-                                               receiver_density,
+                                               ion_density,
                                                z_effective,
                                                b_field)
 
@@ -220,7 +223,7 @@ cdef class BeamCXLine(BeamModel):
 
             effective_rate = cx_rate.evaluate(interaction_energy,
                                               receiver_temperature,
-                                              receiver_density,
+                                              ion_density,
                                               z_effective,
                                               b_field)
 

cherab/core/tests/test_beam.py

@@ -55,32 +55,38 @@ def beam_stopping_rate(self, beam_ion, plasma_ion, charge):
 
 class TestBeam(unittest.TestCase):
 
-    atomic_data = MockAtomicData()
-
-    world = World()
-
-    plasma_density = 1.e19
-    plasma_temperature = 1.e3
-    plasma_species = [(deuterium, 1, plasma_density, plasma_temperature, Vector3D(0, 0, 0))]
-    plasma = build_constant_slab_plasma(length=1, width=1, height=1, electron_density=plasma_density,
-                                        electron_temperature=plasma_temperature,
-                                        plasma_species=plasma_species)
-    plasma.atomic_data = atomic_data
-    plasma.parent = world
-
-    beam = Beam(transform=translate(0.5, 0, 0))
-    beam.atomic_data = atomic_data
-    beam.plasma = plasma
-    beam.attenuator = SingleRayAttenuator(clamp_to_zero=True)
-    beam.energy = 50000
-    beam.power = 1e6
-    beam.temperature = 10
-    beam.element = deuterium
-    beam.parent = world
-    beam.sigma = 0.2
-    beam.divergence_x = 1.
-    beam.divergence_y = 2.
-    beam.length = 10.
+    def setUp(self):
+
+        self.atomic_data = MockAtomicData()
+
+        self.world = World()
+
+        self.plasma_density = 1.e19
+        self.plasma_temperature = 1.e3
+        plasma_species = [(deuterium, 1, self.plasma_density, self.plasma_temperature, Vector3D(0, 0, 0))]
+        plasma = build_constant_slab_plasma(length=1, width=1, height=1,
+                                            electron_density=self.plasma_density,
+                                            electron_temperature=self.plasma_temperature,
+                                            plasma_species=plasma_species)
+        plasma.atomic_data = self.atomic_data
+        plasma.parent = self.world
+
+        beam = Beam(transform=translate(0.5, 0, 0))
+        beam.atomic_data = self.atomic_data
+        beam.plasma = plasma
+        beam.attenuator = SingleRayAttenuator(clamp_to_zero=True)
+        beam.energy = 50000
+        beam.power = 1e6
+        beam.temperature = 10
+        beam.element = deuterium
+        beam.parent = self.world
+        beam.sigma = 0.2
+        beam.divergence_x = 1.
+        beam.divergence_y = 2.
+        beam.length = 10.
+
+        self.plasma = plasma
+        self.beam = beam
 
     def test_beam_density(self):
 

cherab/core/tests/test_beamcxline.py

@@ -36,7 +36,7 @@ class ConstantBeamCXPEC(BeamCXPEC):
     """
 
     def __init__(self, donor_metastable, value):
-        self.donor_metastable = donor_metastable
+        super().__init__(donor_metastable)
         self.value = value
 
     def evaluate(self, energy, temperature, density, z_effective, b_field):
@@ -76,25 +76,33 @@ def wavelength(self, ion, charge, transition):
 
 class TestBeamCXLine(unittest.TestCase):
 
-    world = World()
-
-    atomic_data = MockAtomicData()
-
-    plasma_species = [(deuterium, 1, 1.e19, 200., Vector3D(0, 0, 0))]
-    plasma = build_constant_slab_plasma(length=1, width=1, height=1, electron_density=1e19, electron_temperature=200.,
-                                        plasma_species=plasma_species, b_field=Vector3D(0, 10., 0))
-    plasma.atomic_data = atomic_data
-    plasma.parent = world
-
-    beam = Beam(transform=translate(0.5, 0, 0))
-    beam.atomic_data = atomic_data
-    beam.plasma = plasma
-    beam.attenuator = SingleRayAttenuator(clamp_to_zero=True)
-    beam.energy = 50000
-    beam.power = 1e6
-    beam.temperature = 10
-    beam.element = deuterium
-    beam.parent = world
+    def setUp(self):
+
+        self.world = World()
+
+        self.atomic_data = MockAtomicData()
+
+        plasma_species = [(deuterium, 1, 1.e19, 200., Vector3D(0, 0, 0))]
+        plasma = build_constant_slab_plasma(length=1, width=1, height=1,
+                                            electron_density=1e19,
+                                            electron_temperature=200.,
+                                            plasma_species=plasma_species,
+                                            b_field=Vector3D(0, 10., 0))
+        plasma.atomic_data = self.atomic_data
+        plasma.parent = self.world
+
+        beam = Beam(transform=translate(0.5, 0, 0))
+        beam.atomic_data = self.atomic_data
+        beam.plasma = plasma
+        beam.attenuator = SingleRayAttenuator(clamp_to_zero=True)
+        beam.energy = 50000
+        beam.power = 1e6
+        beam.temperature = 10
+        beam.element = deuterium
+        beam.parent = self.world
+
+        self.plasma = plasma
+        self.beam = beam
 
     def test_default_lineshape(self):
         # setting up the model