Merge pull request #433 from vsnever/enhancement/beam_direction_and_d…

…ensity

Improve density and direction calculation in the neutral beam model

CHANGELOG.md

@@ -13,6 +13,8 @@ New:
 * Replace the coarse numerical constant in the Bremsstrahlung model with an exact expression. (#409)
 * Add the kind attribute to RayTransferPipelineXD that determines whether the ray transfer matrix is multiplied by sensitivity ('power') or not ('radiance'). (#412)
 * Improved parsing of metadata from the ADAS ADF15 'bnd' files for H-like ions. Raises a runtime error if the metadata cannot be parsed. (#424)
+* **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)
 
 Bug fixes:
 * Fix deprecated transforms being cached in LaserMaterial after laser.transform update (#420)

cherab/core/beam/node.pxd

@@ -47,7 +47,7 @@ cdef class Beam(Node):
         Vector3D BEAM_AXIS
         double _energy, _power, _temperature
         Element _element
-        double _divergence_x, _divergence_y
+        double _divergence_x, _divergence_y, _tanxdiv, _tanydiv
         double _length, _sigma
         Plasma _plasma
         AtomicData _atomic_data

cherab/core/beam/node.pyx

@@ -30,7 +30,7 @@ from cherab.core.beam.material cimport BeamMaterial
 from cherab.core.beam.model cimport BeamModel
 from cherab.core.atomic cimport AtomicData, Element
 from cherab.core.utility import Notifier
-from libc.math cimport tan, M_PI
+from libc.math cimport tan, M_PI, sqrt
 
 
 cdef double DEGREES_TO_RADIANS = M_PI / 180
@@ -144,6 +144,7 @@ cdef class Beam(Node):
     :ivar float sigma: The Gaussian beam width at the origin in m.
     :ivar float temperature: The broadening of the beam (eV).
 
+
     .. code-block:: pycon
 
        >>> # This example shows how to initialise and populate a basic beam
@@ -188,6 +189,8 @@ cdef class Beam(Node):
         self._element = element = None             # beam species, an Element object
         self._divergence_x = 0.0                   # beam divergence x (degrees)
         self._divergence_y = 0.0                   # beam divergence y (degrees)
+        self._tanxdiv = 0.0                        # tan(DEGREES_TO_RADIANS * divergence_x)
+        self._tanydiv = 0.0                        # tan(DEGREES_TO_RADIANS * divergence_y)
         self._length = 1.0                         # m
         self._sigma = 0.1                          # m (gaussian beam width at origin)
 
@@ -231,8 +234,25 @@ cdef class Beam(Node):
         return self._attenuator.density(x, y, z)
 
     cpdef Vector3D direction(self, double x, double y, double z):
-        """
-        Calculates the beam direction vector at a point in space.
+        r"""
+        Calculates the beam direction vector at a point in beam coordinate space.
+
+        The beam direction (non-normalised) is calculated as follows (z > 0):
+
+        .. math::
+            e_x = x\frac{(ztg(\alpha_x))^2}{\sigma^2 + (ztg(\alpha_x))^2},
+
+            e_y = y\frac{(ztg(\alpha_y))^2}{\sigma^2 + (ztg(\alpha_y))^2},
+
+            e_z = z,
+
+        where :math:`\sigma` is the Gaussian beam deviation at origin,
+        :math:`\alpha_x` and :math:`\alpha_y` are the beam divergence angles
+        in the x and y dimensions respectively.
+
+        For z <= 0 the beam direction is (0, 0, 1).
+
+        The function returns normalised beam direction.
         
         Note the values of the beam outside of the beam envelope should be
         treated with caution.
@@ -248,9 +268,15 @@ cdef class Beam(Node):
             return self.BEAM_AXIS
 
         # calculate direction from divergence
-        cdef double dx = tan(DEGREES_TO_RADIANS * self._divergence_x)
-        cdef double dy = tan(DEGREES_TO_RADIANS * self._divergence_y)
-        return new_vector3d(dx, dy, 1.0).normalise()
+        cdef double z_tanx_sqr = z * z * self._tanxdiv * self._tanxdiv
+        cdef double z_tany_sqr = z * z * self._tanydiv * self._tanydiv
+        cdef double sigma_sqr = self._sigma * self._sigma
+        cdef double sigma_x_sqr = sigma_sqr + z_tanx_sqr
+        cdef double sigma_y_sqr = sigma_sqr + z_tany_sqr
+        cdef double ex = x * z_tanx_sqr / sigma_x_sqr
+        cdef double ey = y * z_tany_sqr / sigma_y_sqr
+
+        return new_vector3d(ex, ey, z).normalise()
 
     @property
     def energy(self):
@@ -315,6 +341,7 @@ cdef class Beam(Node):
         if value < 0:
             raise ValueError('Beam x divergence cannot be less than zero.')
         self._divergence_x = value
+        self._tanxdiv = tan(DEGREES_TO_RADIANS * value)
         self.notifier.notify()
 
     cdef double get_divergence_x(self):
@@ -329,6 +356,7 @@ cdef class Beam(Node):
         if value < 0:
             raise ValueError('Beam y divergence cannot be less than zero.')
         self._divergence_y = value
+        self._tanydiv = tan(DEGREES_TO_RADIANS * value)
         self.notifier.notify()
 
     cdef double get_divergence_y(self):
@@ -501,10 +529,10 @@ cdef class Beam(Node):
 
         # radii of bounds at the beam origin (z=0) and the beam end (z=length)
         radius_start = num_sigma * self.sigma
-        radius_end = radius_start + self.length * num_sigma * drdz
+        radius_end = num_sigma * sqrt(self.sigma**2 + self.length**2 * drdz**2)
 
         # distance of the cone apex to the beam origin
-        distance_apex = radius_start / (num_sigma * drdz)
+        distance_apex = radius_start * self.length / (radius_end - radius_start)
         cone_height = self.length + distance_apex
 
         # calculate volumes

cherab/core/model/attenuator/singleray.pyx

@@ -37,6 +37,21 @@ cimport cython
 
 # todo: attenuation calculation could be optimised further using memory views etc...
 cdef class SingleRayAttenuator(BeamAttenuator):
+    r"""
+    Calculates beam attenuation in the single-ray approximation.
+    Attenuation is calculated along the beam axis and extrapolated across the beam.
+
+    :param double step: Distance between sample points along the beam axis in meters
+        for beam stopping calculation. Defaults to 0.01.
+    :param bint clamp_to_zero: Omptimises beam density calculation.
+        If True, the beam density outside the clamping range is zero. Defaults to False.
+    :param double clamp_sigma: The clamping range as a factor of beam :math:`\sigma(z)`.
+        Defaults to 5.
+    :param Beam beam: The beam instance to which this attenuator is attached. Defaults to None.
+    :param Plasma plasma: The plasma instance with which this beam interacts. Defaults to None.
+    :param AtomicData atomic_data: The atomic data provider class for this attenuator.
+        Defaults to None.
+    """
 
     def __init__(self, double step=0.01, bint clamp_to_zero=False, double clamp_sigma=5.0, Beam beam=None, Plasma plasma=None, AtomicData atomic_data=None):
 
@@ -86,18 +101,41 @@ cdef class SingleRayAttenuator(BeamAttenuator):
 
     @cython.cdivision(True)
     cpdef double density(self, double x, double y, double z) except? -1e999:
-        """
-        Returns the beam density at the specified point.
-        
-        The point is specified in beam space.
+        r"""
+        Returns the beam density at the specified point in beam coordinate space.
+        The beam density is calculated as follows:
+
+        .. math::
+            n(x, y, z) = \frac{R}{2\pi v_0 \sigma_x\sigma_y} exp\left(-\frac{1}{2}\left(\frac{x^2}{\sigma_x^2}+\frac{y^2}{\sigma_y^2}\right)\right)exp\left(-\int_{0}^{z}\frac{S(z')}{v_0}dz'\right),
+
+            \sigma_x = \sqrt{\sigma^2 + (ztg(\alpha_x))^2}\hspace{0.5cm}\sigma_y = \sqrt{\sigma^2 + (ztg(\alpha_y))^2},
+
+        where :math:`R=\frac{P}{E}` is the particle rate of the beam defined as the power
+        of the beam divided by the kinetic energy of the single particle, :math:`v_0=\sqrt{2E/m}`
+        is the particle speed, :math:`\sigma` is the Gaussian beam deviation at origin,
+        :math:`\alpha_x` and :math:`\alpha_y` are the beam divergence angles in the x and y
+        dimensions respectively, :math:`S(z)` is the composite beam attenuation coefficient due to
+        collisional-radiative interaction with the plasma species:
+
+        .. math::
+            S(z) = \sum_{i=1}^{N}Z_i n_i S_i(E_{int}, n_{i,e}^{(eq)}, T_i),
+
+            n_{i,e}^{(eq)} = \frac{1}{Z_i}\sum_{j=1}^{N}Z_j^2 n_j.
+
+        Here :math:`Z_i` is the charge of the i-th type of plasma ions,
+        :math:`n_i` is density of the i-th type of plasma ions, :math:`N` is the number of type of plasma
+        ions, :math:`E_{int}` is the kinetic energy of the beam atoms in the frame of reference where
+        ions of the i-th type are at rest, :math:`T_{i}` is the temperature of ions of the i-th type.
+
+        The values of partial beam attenuation coefficients, :math:`S_i`, are provided by the atomic data source.
         
         :param x: x coordinate in meters.
         :param y: y coordinate in meters. 
         :param z: z coordinate in meters.
         :return: Density in m^-3. 
         """
 
-        cdef double sigma_x, sigma_y, norm_radius_sqr, gaussian_sample
+        cdef double sigma0_sqr, sigma_x, sigma_y, norm_radius_sqr, gaussian_sample
 
         # use cached data if available
         if self._stopping_data is None:
@@ -107,8 +145,9 @@ cdef class SingleRayAttenuator(BeamAttenuator):
             self._calc_attenuation()
 
         # calculate beam width
-        sigma_x = self._beam.get_sigma() + z * self._tanxdiv
-        sigma_y = self._beam.get_sigma() + z * self._tanydiv
+        sigma0_sqr = self._beam.get_sigma()**2
+        sigma_x = sqrt(sigma0_sqr + (z * self._tanxdiv)**2)
+        sigma_y = sqrt(sigma0_sqr + (z * self._tanydiv)**2)
 
         # normalised radius squared
         norm_radius_sqr = ((x / sigma_x)**2 + (y / sigma_y)**2)

cherab/core/tests/test_beam.py

@@ -0,0 +1,156 @@
+# 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 World, Vector3D, translate
+
+from cherab.core import Beam
+from cherab.core.atomic import AtomicData, BeamStoppingRate
+from cherab.core.atomic import deuterium
+from cherab.tools.plasmas.slab import build_constant_slab_plasma
+from cherab.core.model import SingleRayAttenuator
+
+from cherab.core.utility import EvAmuToMS, EvToJ
+
+
+class ConstantBeamStoppingRate(BeamStoppingRate):
+    """
+    Constant beam CX PEC for test purpose.
+    """
+
+    def __init__(self, donor_metastable, value):
+        self.donor_metastable = donor_metastable
+        self.value = value
+
+    def evaluate(self, energy, density, temperature):
+
+        return self.value
+
+
+class MockAtomicData(AtomicData):
+    """Fake atomic data for test purpose."""
+
+    def beam_stopping_rate(self, beam_ion, plasma_ion, charge):
+
+        return ConstantBeamStoppingRate(1, 1.e-13)
+
+
+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 test_beam_density(self):
+
+        z0 = 0.8
+        x0, y0 = 0.5, 0.5
+
+        density_on_axis = self.beam.density(0, 0, z0)
+        density_off_axis = self.beam.density(x0, y0, z0)
+        density_outside_beam = self.beam.density(0, 0, -1)
+
+        # validating
+
+        speed = EvAmuToMS.to(self.beam.energy)
+        # constant stopping rate
+        stopping_rate = self.atomic_data.beam_stopping_rate(deuterium, deuterium, 1)(0, 0, 0)
+        attenuation_factor = np.exp(-z0 * self.plasma_density * stopping_rate / speed)
+
+        beam_particle_rate = self.beam.power / EvToJ.to(self.beam.energy * deuterium.atomic_weight)
+
+        sigma0_sqr = self.beam.sigma**2
+        tanxdiv = np.tan(np.deg2rad(self.beam.divergence_x))
+        tanydiv = np.tan(np.deg2rad(self.beam.divergence_y))
+        sigma_x = np.sqrt(sigma0_sqr + (z0 * tanxdiv)**2)
+        sigma_y = np.sqrt(sigma0_sqr + (z0 * tanydiv)**2)
+
+        norm_radius_sqr = ((x0 / sigma_x)**2 + (y0 / sigma_y)**2)
+
+        gaussian_sample_on_axis = 1. / (2 * np.pi * sigma_x * sigma_y)
+        gaussian_sample_off_axis = np.exp(-0.5 * norm_radius_sqr) / (2 * np.pi * sigma_x * sigma_y)
+
+        test_density_on_axis = beam_particle_rate / speed * gaussian_sample_on_axis * attenuation_factor
+        test_density_off_axis = beam_particle_rate / speed * gaussian_sample_off_axis * attenuation_factor
+
+        self.assertAlmostEqual(density_on_axis / test_density_on_axis, 1., delta=1.e-12,
+                               msg='Beam.density() gives a wrong value on the beam axis.')
+        self.assertAlmostEqual(density_off_axis / test_density_off_axis, 1., delta=1.e-12,
+                               msg='Beam.density() gives a wrong value off the beam axis.')
+        self.assertEqual(density_outside_beam, 0,
+                         msg='Beam.density() gives a non-zero value outside beam.')
+
+    def test_beam_direction(self):
+        # setting up the model
+
+        z0 = 0.8
+        x0, y0 = 0.5, 0.5
+
+        direction_on_axis = self.beam.direction(0, 0, z0)
+        direction_off_axis = self.beam.direction(x0, y0, z0)
+        direction_outside_beam = self.beam.direction(0, 0, -1)
+
+        # validating
+
+        sigma0_sqr = self.beam.sigma**2
+        z_tanx_sqr = (z0 * np.tan(np.deg2rad(self.beam.divergence_x)))**2
+        z_tany_sqr = (z0 * np.tan(np.deg2rad(self.beam.divergence_y)))**2
+
+        ex = x0 * z_tanx_sqr / (sigma0_sqr + z_tanx_sqr)
+        ey = y0 * z_tany_sqr / (sigma0_sqr + z_tany_sqr)
+        ez = z0
+
+        test_direction_off_axis = Vector3D(ex, ey, ez).normalise()
+
+        self.assertEqual(direction_on_axis, Vector3D(0, 0, 1),
+                         msg='Beam.density() gives a wrong value on the beam axis.')
+        for v, test_v in zip(direction_off_axis, test_direction_off_axis):
+            self.assertAlmostEqual(v, test_v, delta=1.e-12,
+                                   msg='Beam.direction() gives a wrong value off the beam axis.')
+        self.assertEqual(direction_outside_beam, Vector3D(0, 0, 1),
+                         msg='Beam.density() gives a non-zero value outside beam.')
+
+
+if __name__ == '__main__':
+    unittest.main()

cherab/core/tests/test_beamcxline.py

@@ -58,7 +58,7 @@ def evaluate(self, energy, density, temperature):
         return self.value
 
 
-class TestAtomicData(AtomicData):
+class MockAtomicData(AtomicData):
     """Fake atomic data for test purpose."""
 
     def beam_cx_pec(self, donor_ion, receiver_ion, receiver_charge, transition):
@@ -78,7 +78,7 @@ class TestBeamCXLine(unittest.TestCase):
 
     world = World()
 
-    atomic_data = TestAtomicData()
+    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.,

demos/emission_models/beam_emission_spectrum.py

@@ -69,9 +69,9 @@
 beam_energy = 110000  # keV
 beam_current = 10  # A
 beam_sigma = 0.05
-beam_divergence = 0.5
+beam_divergence = 1.3  # degree
 beam_length = 3.0
-beam_temperature = 20
+beam_temperature = 1.0
 
 bes_full_model = BeamEmissionLine(Line(deuterium, 0, (3, 2)),
                                   sigma_to_pi=SIGMA_TO_PI, sigma1_to_sigma0=SIGMA1_TO_SIGMA0,

docs/source/models/beam/beam_attenuator.rst

@@ -0,0 +1,10 @@
+
+Beam Attenuation
+----------------
+
+
+Single Ray Attenuator
+^^^^^^^^^^^^^^^^^^^^^
+
+.. autoclass:: cherab.core.model.attenuator.singleray.SingleRayAttenuator
+   :members: