diff --git a/CHANGELOG.md b/CHANGELOG.md index 39560584..035c5b8c 100644 --- a/CHANGELOG.md +++ b/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: diff --git a/cherab/core/atomic/interface.pxd b/cherab/core/atomic/interface.pxd index 66d8faf8..f6b9b914 100644 --- a/cherab/core/atomic/interface.pxd +++ b/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) diff --git a/cherab/core/atomic/interface.pyx b/cherab/core/atomic/interface.pyx index 3e44d52c..91417ef4 100644 --- a/cherab/core/atomic/interface.pyx +++ b/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 diff --git a/cherab/core/atomic/rates.pxd b/cherab/core/atomic/rates.pxd index ac8d69ce..41871191 100644 --- a/cherab/core/atomic/rates.pxd +++ b/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: diff --git a/cherab/core/atomic/rates.pyx b/cherab/core/atomic/rates.pyx index 52510e4b..37225e73 100644 --- a/cherab/core/atomic/rates.pyx +++ b/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: diff --git a/cherab/core/model/plasma/__init__.pxd b/cherab/core/model/plasma/__init__.pxd index 6971f828..c398fd30 100644 --- a/cherab/core/model/plasma/__init__.pxd +++ b/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 diff --git a/cherab/core/model/plasma/__init__.py b/cherab/core/model/plasma/__init__.py index 52766032..436060b5 100644 --- a/cherab/core/model/plasma/__init__.py +++ b/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 diff --git a/cherab/core/model/plasma/thermal_cx.pxd b/cherab/core/model/plasma/thermal_cx.pxd new file mode 100644 index 00000000..636fddaa --- /dev/null +++ b/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 diff --git a/cherab/core/model/plasma/thermal_cx.pyx b/cherab/core/model/plasma/thermal_cx.pyx new file mode 100644 index 00000000..1325919c --- /dev/null +++ b/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 ''.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 diff --git a/cherab/core/tests/test_line_emission.py b/cherab/core/tests/test_line_emission.py new file mode 100644 index 00000000..b55af871 --- /dev/null +++ b/cherab/core/tests/test_line_emission.py @@ -0,0 +1,324 @@ +# 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, translate +from raysect.optical import World, Spectrum, Ray + +from cherab.core.atomic import Line, AtomicData, ImpactExcitationPEC, RecombinationPEC, ThermalCXPEC +from cherab.core.atomic import deuterium, carbon +from cherab.tools.plasmas.slab import build_constant_slab_plasma +from cherab.core.model import ExcitationLine, RecombinationLine, ThermalCXLine, GaussianLine, ZeemanTriplet + + +class ConstantImpactExcitationPEC(ImpactExcitationPEC): + """ + Constant electron impact excitation PEC for test purpose. + """ + + def __init__(self, value): + self.value = value + + def evaluate(self, density, temperature): + + return self.value + + +class ConstantRecombinationPEC(RecombinationPEC): + """ + Constant recombination PEC for test purpose. + """ + + def __init__(self, value): + self.value = value + + def evaluate(self, density, temperature): + + return self.value + + +class ConstantThermalCXPEC(ThermalCXPEC): + """ + Constant recombination PEC for test purpose. + """ + + def __init__(self, value): + self.value = value + + def evaluate(self, electron_density, electron_temperature, donor_temperature): + + return self.value + + +class MockAtomicData(AtomicData): + """Fake atomic data for test purpose.""" + + def impact_excitation_pec(self, ion, charge, transition): + + return ConstantImpactExcitationPEC(1.4e-39) + + def recombination_pec(self, ion, charge, transition): + + return ConstantRecombinationPEC(8.e-41) + + def thermal_cx_pec(self, donor_ion, donor_charge, receiver_ion, receiver_charge, transition): + + return ConstantThermalCXPEC(1.2e-46) + + def wavelength(self, ion, charge, transition): + + return 529.27 + + +class TestExcitationLine(unittest.TestCase): + + def setUp(self): + + self.world = World() + + self.atomic_data = MockAtomicData() + + plasma_species = [(carbon, 5, 2.e18, 800., Vector3D(0, 0, 0))] + self.slab_length = 1.2 + self.plasma = build_constant_slab_plasma(length=self.slab_length, width=1, height=1, + electron_density=1e19, electron_temperature=1000., + plasma_species=plasma_species, b_field=Vector3D(0, 10., 0)) + self.plasma.atomic_data = self.atomic_data + self.plasma.parent = self.world + + def test_default_lineshape(self): + # setting up the model + line = Line(carbon, 5, (8, 7)) + self.plasma.models = [ExcitationLine(line)] + wavelength = self.atomic_data.wavelength(line.element, line.charge, line.transition) + + # observing + origin = Point3D(1.5, 0, 0) + direction = Vector3D(-1, 0, 0) + ray = Ray(origin=origin, direction=direction, + min_wavelength=wavelength - 1.5, max_wavelength=wavelength + 1.5, bins=512) + excit_spectrum = ray.trace(self.world) + + # validating + ne = self.plasma.electron_distribution.density(0.5, 0, 0) # constant slab + te = self.plasma.electron_distribution.effective_temperature(0.5, 0, 0) + rate = self.atomic_data.impact_excitation_pec(line.element, line.charge, line.transition)(ne, te) + target_species = self.plasma.composition.get(line.element, line.charge) + ni = target_species.distribution.density(0.5, 0, 0) + radiance = 0.25 / np.pi * rate * ni * ne * self.slab_length + + gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma) + spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) + spectrum = gaussian_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) + + for i in range(ray.bins): + self.assertAlmostEqual(excit_spectrum.samples[i], spectrum.samples[i], delta=1e-8, + msg='ExcitationLine model gives a wrong value at {} nm.'.format(spectrum.wavelengths[i])) + + def test_custom_lineshape(self): + # setting up the model + line = Line(carbon, 5, (8, 7)) + self.plasma.models = [ExcitationLine(line, lineshape=ZeemanTriplet)] + wavelength = self.atomic_data.wavelength(line.element, line.charge, line.transition) + + # observing + origin = Point3D(1.5, 0, 0) + direction = Vector3D(-1, 0, 0) + ray = Ray(origin=origin, direction=direction, + min_wavelength=wavelength - 1.5, max_wavelength=wavelength + 1.5, bins=512) + excit_spectrum = ray.trace(self.world) + + # validating + ne = self.plasma.electron_distribution.density(0.5, 0, 0) # constant slab + te = self.plasma.electron_distribution.effective_temperature(0.5, 0, 0) + rate = self.atomic_data.impact_excitation_pec(line.element, line.charge, line.transition)(ne, te) + target_species = self.plasma.composition.get(line.element, line.charge) + ni = target_species.distribution.density(0.5, 0, 0) + radiance = 0.25 / np.pi * rate * ni * ne * self.slab_length + + zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma) + spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) + spectrum = zeeman_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) + + for i in range(ray.bins): + self.assertAlmostEqual(excit_spectrum.samples[i], spectrum.samples[i], delta=1e-8, + msg='ExcitationLine model gives a wrong value at {} nm.'.format(spectrum.wavelengths[i])) + + +class TestRecombinationLine(unittest.TestCase): + + def setUp(self): + + self.world = World() + + self.atomic_data = MockAtomicData() + + plasma_species = [(carbon, 6, 1.67e18, 800., Vector3D(0, 0, 0))] + self.slab_length = 1.2 + self.plasma = build_constant_slab_plasma(length=self.slab_length, width=1, height=1, + electron_density=1e19, electron_temperature=1000., + plasma_species=plasma_species, b_field=Vector3D(0, 10., 0)) + self.plasma.atomic_data = self.atomic_data + self.plasma.parent = self.world + + def test_default_lineshape(self): + # setting up the model + line = Line(carbon, 5, (8, 7)) + self.plasma.models = [RecombinationLine(line)] + wavelength = self.atomic_data.wavelength(line.element, line.charge, line.transition) + + # observing + origin = Point3D(1.5, 0, 0) + direction = Vector3D(-1, 0, 0) + ray = Ray(origin=origin, direction=direction, + min_wavelength=wavelength - 1.5, max_wavelength=wavelength + 1.5, bins=512) + recomb_spectrum = ray.trace(self.world) + + # validating + ne = self.plasma.electron_distribution.density(0.5, 0, 0) # constant slab + te = self.plasma.electron_distribution.effective_temperature(0.5, 0, 0) + rate = self.atomic_data.recombination_pec(line.element, line.charge, line.transition)(ne, te) + target_species = self.plasma.composition.get(line.element, line.charge + 1) + ni = target_species.distribution.density(0.5, 0, 0) + radiance = 0.25 / np.pi * rate * ni * ne * self.slab_length + + gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma) + spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) + spectrum = gaussian_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) + + for i in range(ray.bins): + self.assertAlmostEqual(recomb_spectrum.samples[i], spectrum.samples[i], delta=1e-8, + msg='RecombinationLine model gives a wrong value at {} nm.'.format(spectrum.wavelengths[i])) + + def test_custom_lineshape(self): + # setting up the model + line = Line(carbon, 5, (8, 7)) + self.plasma.models = [RecombinationLine(line, lineshape=ZeemanTriplet)] + wavelength = self.atomic_data.wavelength(line.element, line.charge, line.transition) + + # observing + origin = Point3D(1.5, 0, 0) + direction = Vector3D(-1, 0, 0) + ray = Ray(origin=origin, direction=direction, + min_wavelength=wavelength - 1.5, max_wavelength=wavelength + 1.5, bins=512) + recomb_spectrum = ray.trace(self.world) + + # validating + ne = self.plasma.electron_distribution.density(0.5, 0, 0) # constant slab + te = self.plasma.electron_distribution.effective_temperature(0.5, 0, 0) + rate = self.atomic_data.recombination_pec(line.element, line.charge, line.transition)(ne, te) + target_species = self.plasma.composition.get(line.element, line.charge + 1) + ni = target_species.distribution.density(0.5, 0, 0) + radiance = 0.25 / np.pi * rate * ni * ne * self.slab_length + + zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma) + spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) + spectrum = zeeman_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) + + for i in range(ray.bins): + self.assertAlmostEqual(recomb_spectrum.samples[i], spectrum.samples[i], delta=1e-8, + msg='RecombinationLine model gives a wrong value at {} nm.'.format(spectrum.wavelengths[i])) + + +class TestThermalCXLine(unittest.TestCase): + + def setUp(self): + + self.world = World() + + self.atomic_data = MockAtomicData() + + plasma_species = [(carbon, 6, 1.67e18, 800., Vector3D(0, 0, 0)), + (deuterium, 0, 1.e19, 100., Vector3D(0, 0, 0))] + self.slab_length = 1.2 + self.plasma = build_constant_slab_plasma(length=self.slab_length, width=1, height=1, + electron_density=1e19, electron_temperature=1000., + plasma_species=plasma_species, b_field=Vector3D(0, 10., 0)) + self.plasma.atomic_data = self.atomic_data + self.plasma.parent = self.world + + def test_default_lineshape(self): + # setting up the model + line = Line(carbon, 5, (8, 7)) + self.plasma.models = [ThermalCXLine(line)] + wavelength = self.atomic_data.wavelength(line.element, line.charge, line.transition) + + # observing + origin = Point3D(1.5, 0, 0) + direction = Vector3D(-1, 0, 0) + ray = Ray(origin=origin, direction=direction, + min_wavelength=wavelength - 1.5, max_wavelength=wavelength + 1.5, bins=512) + thermalcx_spectrum = ray.trace(self.world) + + # validating + ne = self.plasma.electron_distribution.density(0.5, 0, 0) # constant slab + te = self.plasma.electron_distribution.effective_temperature(0.5, 0, 0) + donor_species = self.plasma.composition.get(deuterium, 0) + donor_density = donor_species.distribution.density(0.5, 0, 0) + donor_temperature = donor_species.distribution.effective_temperature(0.5, 0, 0) + rate = self.atomic_data.thermal_cx_pec(deuterium, 0, line.element, line.charge, line.transition)(ne, te, donor_temperature) + target_species = self.plasma.composition.get(line.element, line.charge + 1) + receiver_density = target_species.distribution.density(0.5, 0, 0) + radiance = 0.25 / np.pi * rate * receiver_density * donor_density * self.slab_length + + gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma) + spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) + spectrum = gaussian_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) + + for i in range(ray.bins): + self.assertAlmostEqual(thermalcx_spectrum.samples[i], spectrum.samples[i], delta=1e-8, + msg='ThermalCXLine model gives a wrong value at {} nm.'.format(spectrum.wavelengths[i])) + + def test_custom_lineshape(self): + # setting up the model + line = Line(carbon, 5, (8, 7)) + self.plasma.models = [ThermalCXLine(line, lineshape=ZeemanTriplet)] + wavelength = self.atomic_data.wavelength(line.element, line.charge, line.transition) + + # observing + origin = Point3D(1.5, 0, 0) + direction = Vector3D(-1, 0, 0) + ray = Ray(origin=origin, direction=direction, + min_wavelength=wavelength - 1.5, max_wavelength=wavelength + 1.5, bins=512) + thermalcx_spectrum = ray.trace(self.world) + + # validating + ne = self.plasma.electron_distribution.density(0.5, 0, 0) # constant slab + te = self.plasma.electron_distribution.effective_temperature(0.5, 0, 0) + donor_species = self.plasma.composition.get(deuterium, 0) + donor_density = donor_species.distribution.density(0.5, 0, 0) + donor_temperature = donor_species.distribution.effective_temperature(0.5, 0, 0) + rate = self.atomic_data.thermal_cx_pec(deuterium, 0, line.element, line.charge, line.transition)(ne, te, donor_temperature) + target_species = self.plasma.composition.get(line.element, line.charge + 1) + receiver_density = target_species.distribution.density(0.5, 0, 0) + radiance = 0.25 / np.pi * rate * receiver_density * donor_density * self.slab_length + + zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma) + spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) + spectrum = zeeman_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) + + for i in range(ray.bins): + self.assertAlmostEqual(thermalcx_spectrum.samples[i], spectrum.samples[i], delta=1e-8, + msg='ThermalCXLine model gives a wrong value at {} nm.'.format(spectrum.wavelengths[i])) + + +if __name__ == '__main__': + unittest.main() diff --git a/cherab/openadas/install.py b/cherab/openadas/install.py index d634cde7..53e68899 100644 --- a/cherab/openadas/install.py +++ b/cherab/openadas/install.py @@ -1,7 +1,7 @@ -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# 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"); @@ -21,9 +21,11 @@ import os import urllib.parse import urllib.request +import numpy as np from cherab.openadas import repository from cherab.openadas.parse import * +from cherab.core.atomic import hydrogen from cherab.core.utility import RecursiveDict, PerCm3ToPerM3, Cm3ToM3 @@ -264,6 +266,13 @@ def install_adf15(element, ionisation, file_path, download=False, repository_pat # decode file and write out rates rates, wavelengths = parse_adf15(element, ionisation, path, header_format=header_format) + + if 'thermalcx' in rates: + cx_rates = rates.pop('thermalcx') + # CX rates for Tdon = Trec (2D function of Ne, Te) + # converting to 3D function of Ne, Te, Tdon + repository.update_pec_thermal_cx_rates(_thermalcx_adf15_2dto3d_converter(cx_rates)) + repository.update_pec_rates(rates, repository_path) repository.update_wavelengths(wavelengths, repository_path) @@ -402,3 +411,23 @@ def _notation_adf11_adas2cherab(rate_adas, filetype): return rate_cherab +def _thermalcx_adf15_2dto3d_converter(rates): + """ + Converts thermal CX PEC rates parsed from a standard ADF 15 file + to the format supported by the repository. + + In the standard ADF 15 file, it is assumed that the donor is H0 and Tdon = Trec. + """ + new_rates = RecursiveDict() + for element, charge_states in rates.items(): + for charge, transitions in charge_states.items(): + for transition, rate in transitions.items(): + data = np.empty((len(rate['ne']), len(rate['te']), 2)) + data[:, :, :] = rate['rate'][:, :, None] + new_rate = {'ne': rate['ne'], + 'te': rate['te'], + 'td': np.array([0.01, 10000]), + 'rate': data} + new_rates[hydrogen][0][element][charge + 1][transition] = new_rate + + return new_rates diff --git a/cherab/openadas/openadas.py b/cherab/openadas/openadas.py index 7e850bc4..5d42935a 100644 --- a/cherab/openadas/openadas.py +++ b/cherab/openadas/openadas.py @@ -386,6 +386,48 @@ def recombination_pec(self, ion, charge, transition): return RecombinationPEC(wavelength, data, extrapolate=self._permit_extrapolation) + def thermal_cx_pec(self, donor_element, donor_charge, receiver_element, receiver_charge, transition): + """ + Thermal CX photon emission coefficient for a given species. + + Open-ADAS data is interpolated with cubic spline in log-log space. + Nearest neighbour extrapolation is used when permit_extrapolation is True. + + :param donor_element: Element object defining the donor ion type. + :param donor_charge: Charge state of the donor ion. + :param receiver_element: Element object defining the receiver ion type. + :param receiver_charge: Charge state of the receiver ion. + :param transition: Tuple containing (initial level, final level) of the receiver + in charge state receiver_charge - 1. + :return: Thermal charge exchange photon emission coefficient in W.m^3 + as a function of electron density, electron temperature and donor temperature. + """ + + # extract elements from isotopes because there are no isotope rates in ADAS + if isinstance(donor_element, Isotope): + donor_element = donor_element.element + + if isinstance(receiver_element, Isotope): + receiver_element = receiver_element.element + + try: + # read thermal CX rate from json file in the repository + data = repository.get_pec_thermal_cx_rate(donor_element, donor_charge, + receiver_element, receiver_charge, + transition, + repository_path=self._data_path) + + except RuntimeError: + if self._missing_rates_return_null: + return NullThermalCXPEC() + raise + + # obtain isotope's rest wavelength for a given transition + # the wavelength is used ot convert the PEC from photons/s/m3 to W/m3 + wavelength = self.wavelength(receiver_element, receiver_charge - 1, transition) + + return ThermalCXPEC(wavelength, data, extrapolate=self._permit_extrapolation) + def line_radiated_power_rate(self, ion, charge): """ Line radiated power coefficient for a given species. diff --git a/cherab/openadas/parse/adf15.py b/cherab/openadas/parse/adf15.py index 51bbcf6b..269edd0a 100644 --- a/cherab/openadas/parse/adf15.py +++ b/cherab/openadas/parse/adf15.py @@ -1,6 +1,6 @@ -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# 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"); @@ -126,7 +126,7 @@ def _scrape_metadata_hydrogen(file, element, charge): elif rate_type_adas == 'RECOM': rate_type = 'recombination' elif rate_type_adas == 'CHEXC': - rate_type = 'cx_thermal' + rate_type = 'thermalcx' else: raise ValueError("Unrecognised rate type - {}".format(rate_type_adas)) @@ -169,7 +169,7 @@ def _scrape_metadata_hydrogen_like(file, element, charge): elif rate_type_adas == 'RECOM': rate_type = 'recombination' elif rate_type_adas == 'CHEXC': - rate_type = 'cx_thermal' + rate_type = 'thermalcx' else: raise ValueError("Unrecognised rate type - {}".format(rate_type_adas)) @@ -237,7 +237,7 @@ def _scrape_metadata_full(file, element, charge): elif rate_type_adas == 'RECOM': rate_type = 'recombination' elif rate_type_adas == 'CHEXC': - rate_type = 'cx_thermal' + rate_type = 'thermalcx' else: raise ValueError("Unrecognised rate type - {}".format(rate_type_adas)) diff --git a/cherab/openadas/rates/pec.pxd b/cherab/openadas/rates/pec.pxd index 46a54cbe..57bc7560 100644 --- a/cherab/openadas/rates/pec.pxd +++ b/cherab/openadas/rates/pec.pxd @@ -19,7 +19,7 @@ from cherab.core cimport ImpactExcitationPEC as CoreImpactExcitationPEC from cherab.core cimport RecombinationPEC as CoreRecombinationPEC from cherab.core cimport ThermalCXPEC as CoreThermalCXPEC -from cherab.core.math cimport Function2D +from cherab.core.math cimport Function2D, Function3D cdef class ImpactExcitationPEC(CoreImpactExcitationPEC): @@ -48,8 +48,14 @@ cdef class NullRecombinationPEC(CoreRecombinationPEC): pass -# cdef class CXThermalRate(CoreCXThermalRate): -# pass -# -# cdef class ThermalCXRate(CoreThermalCXRate): -# pass +cdef class ThermalCXPEC(CoreThermalCXPEC): + + cdef: + readonly dict raw_data + readonly double wavelength + readonly tuple density_range, temperature_range, donor_temperature_range + Function3D _rate + + +cdef class NullThermalCXPEC(CoreThermalCXPEC): + pass diff --git a/cherab/openadas/rates/pec.pyx b/cherab/openadas/rates/pec.pyx index 60844452..b10b0547 100644 --- a/cherab/openadas/rates/pec.pyx +++ b/cherab/openadas/rates/pec.pyx @@ -20,7 +20,7 @@ import numpy as np from libc.math cimport INFINITY, log10 -from raysect.core.math.function.float cimport Interpolator2DArray +from raysect.core.math.function.float cimport Interpolator2DArray, Interpolator3DArray from cherab.core.utility.conversion import PhotonToJ @@ -148,3 +148,54 @@ cdef class NullRecombinationPEC(CoreRecombinationPEC): cpdef double evaluate(self, double density, double temperature) except? -1e999: return 0.0 + + +cdef class ThermalCXPEC(CoreThermalCXPEC): + + def __init__(self, double wavelength, dict data, extrapolate=False): + """ + :param wavelength: Resting wavelength of corresponding emission line in nm. + :param data: Dictionary containing rate data. + :param extrapolate: Enable nearest-neighbour extrapolation (default=False). + """ + + self.wavelength = wavelength + self.raw_data = data + + # unpack + ne = data['ne'] + te = data['te'] + td = data['td'] + rate = data['rate'] + + # pre-convert data to W m^3 from Photons s^-1 cm^3 prior to interpolation + rate = np.log10(PhotonToJ.to(rate, wavelength)) + + # store limits of data + self.density_range = ne.min(), ne.max() + self.temperature_range = te.min(), te.max() + self.donor_temperature_range = td.min(), td.max() + + # interpolate rate + # using nearest extrapolation to avoid infinite values at 0 for some rates + extrapolation_type = 'nearest' if extrapolate else 'none' + self._rate = Interpolator3DArray(np.log10(ne), np.log10(te), np.log10(td), rate, 'cubic', extrapolation_type, INFINITY, INFINITY, INFINITY) + + cpdef double evaluate(self, double electron_density, double electron_temperature, double donor_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 or donor_temperature <= 0: + return 0 + + # calculate rate and convert from log10 space to linear space + return 10 ** self._rate.evaluate(log10(electron_density), log10(electron_temperature), log10(donor_temperature)) + + +cdef class NullThermalCXPEC(CoreThermalCXPEC): + """ + A PEC rate that always returns zero. + Needed for use cases where the required atomic data is missing. + """ + + cpdef double evaluate(self, double electron_density, double electron_temperature, double donor_temperature) except? -1e999: + return 0.0 diff --git a/cherab/openadas/repository/create.py b/cherab/openadas/repository/create.py index 6b9a4f33..b42e5651 100644 --- a/cherab/openadas/repository/create.py +++ b/cherab/openadas/repository/create.py @@ -147,6 +147,7 @@ def populate(download=True, repository_path=None, adas_path=None): (carbon, 0, 'adf15/pec96#c/pec96#c_vsu#c0.dat'), (carbon, 1, 'adf15/pec96#c/pec96#c_vsu#c1.dat'), (carbon, 2, 'adf15/pec96#c/pec96#c_vsu#c2.dat'), + (carbon, 5, 'adf15/pec96#c/pec96#c_pju#c5.dat'), # (neon, 0, 'adf15/pec96#ne/pec96#ne_pju#ne0.dat'), #TODO: OPENADAS DATA CORRUPT # (neon, 1, 'adf15/pec96#ne/pec96#ne_pju#ne1.dat'), #TODO: OPENADAS DATA CORRUPT (nitrogen, 0, 'adf15/pec96#n/pec96#n_vsu#n0.dat'), diff --git a/cherab/openadas/repository/pec.py b/cherab/openadas/repository/pec.py index 13fc055e..ef933b23 100644 --- a/cherab/openadas/repository/pec.py +++ b/cherab/openadas/repository/pec.py @@ -90,40 +90,37 @@ def add_pec_recombination_rate(element, charge, transition, rate, repository_pat }, repository_path) -def add_pec_thermalcx_rate(element, charge, transition, rate, repository_path=None): +def add_pec_thermal_cx_rate(donor_element, donor_charge, receiver_element, receiver_charge, transition, rate, repository_path=None): """ Adds a single PEC thermal charge exchange rate to the repository. - If adding multiple rate, consider using the update_pec_rates() function + If adding multiple rate, consider using the update_pec_thermal_rates() function instead. The update function avoid repeatedly opening and closing the rate files. - :param element: Plasma species (Element/Isotope). - :param charge: Charge of the plasma species. + :param donor_element: Electron donor plasma species (Element/Isotope). + :param donor_charge: Electron donor charge. + :param receiver_element: Electron receiver plasma species (Element/Isotope). + :param receiver_charge: Electron receiver charge. :param transition: Tuple containing (initial level, final level). :param rate: Thermal CX PEC dictionary containing the following entries: | 'ne': array-like of size (N) with electron density in m^-3, | 'te': array-like of size (M) with electron temperature in eV, - | 'rate': array-like of size (N, M) with thermal CX PEC in photon.m^3.s^-1. + | 'td': array-like of size (K) with donor temperature in eV, + | 'rate': array-like of size (N, M, K) with thermal CX PEC in photon.m^3.s^-1. :param repository_path: Path to the atomic data repository. """ + rates2update = RecursiveDict() + rates2update[donor_element][donor_charge][receiver_element][receiver_charge][transition] = rate - update_pec_rates({ - 'thermalcx': { - element: { - charge: { - transition: rate - } - } - } - }, repository_path) + update_pec_thermal_cx_rates(rates2update.freeze(), repository_path) def update_pec_rates(rates, repository_path=None): """ - Updates the PEC files /pec///.json. + Updates excitation and recombination PEC files /pec///.json. in the atomic data repository. File contains multiple PECs, indexed by the transition. @@ -131,7 +128,7 @@ def update_pec_rates(rates, repository_path=None): :param rates: Dictionary in the form: | { : { : { : { : } } } }, where - | is the one of the following PEC types: 'excitation', 'recombination', 'thermalcx'. + | is the one of the following PEC types: 'excitation', 'recombination'. | is the plasma species (Element/Isotope). | is the charge of the plasma species. | is the tuple containing (initial level, final level). @@ -145,8 +142,7 @@ def update_pec_rates(rates, repository_path=None): valid_classes = [ 'excitation', - 'recombination', - 'thermalcx' + 'recombination' ] repository_path = repository_path or DEFAULT_REPOSITORY_PATH @@ -210,6 +206,99 @@ def update_pec_rates(rates, repository_path=None): json.dump(content, f, indent=2, sort_keys=True) +def update_pec_thermal_cx_rates(rates, repository_path=None): + """ + Updates thermal CX PEC files /pec/thermal_cx////.json + in the atomic data repository. + + File contains multiple PECs, indexed by the transition. + + :param rates: Dictionary in the form: + + | { : { : { : { : { : } } } } }, where + | is the electron donor species (Element/Isotope). + | is the electron donor charge. + | is the electron receiver species (Element/Isotope). + | is the electron receiver charge. + | is the tuple containing (initial level, final level). + | is the PEC dictionary containing the following entries: + | 'ne': array-like of size (N) with electron density in m^-3, + | 'te': array-like of size (M) with electron temperature in eV, + | 'td': array-like of size (K) with donor temperature in eV, + | 'rate': array-like of size (N, M, K) with PEC in photon.m^3.s^-1. + + :param repository_path: Path to the atomic data repository. + """ + repository_path = repository_path or DEFAULT_REPOSITORY_PATH + + for donor_element, donor_charge_states in rates.items(): + for donor_charge, receiver_elements in donor_charge_states.items(): + for receiver_element, receiver_charge_states in receiver_elements.items(): + for receiver_charge, transitions in receiver_charge_states.items(): + + # sanitise and validate + if not isinstance(donor_element, Element): + raise TypeError('The donor element must be an Element object.') + + if not valid_charge(donor_element, donor_charge + 1): + raise ValueError('Donor charge state is equal to or larger than the number of protons in the element.') + + if not isinstance(receiver_element, Element): + raise TypeError('The receiver element must be an Element object.') + + if not valid_charge(receiver_element, receiver_charge): + raise ValueError('Receiver charge state is larger than the number of protons in the element.') + + rate_path = 'pec/thermal_cx/{0}/{1}/{2}/{3}.json'.format(donor_element.symbol.lower(), donor_charge, + receiver_element.symbol.lower(), receiver_charge) + path = os.path.join(repository_path, rate_path) + + # read in any existing rates + try: + with open(path, 'r') as f: + content = RecursiveDict.from_dict(json.load(f)) + except FileNotFoundError: + content = RecursiveDict() + + # add/replace data for a transition + for transition, data in transitions.items(): + key = encode_transition(transition) + + # sanitise/validate data + ne = np.array(data['ne'], np.float64) + te = np.array(data['te'], np.float64) + td = np.array(data['td'], np.float64) + rate = np.array(data['rate'], np.float64) + + if ne.ndim != 1: + raise ValueError('Electron density array must be a 1D array.') + + if te.ndim != 1: + raise ValueError('Electron temperature array must be a 1D array.') + + if td.ndim != 1: + raise ValueError('Donor temperature array must be a 1D array.') + + if (ne.shape[0], te.shape[0], td.shape[0]) != rate.shape: + raise ValueError('Electron density, electron temperature, donor temperature' + ' and rate data arrays have inconsistent sizes.') + + content[key] = { + 'ne': ne.tolist(), + 'te': te.tolist(), + 'td': td.tolist(), + 'rate': rate.tolist() + } + + # create directory structure if missing + directory = os.path.dirname(path) + os.makedirs(directory, exist_ok=True) + + # write new data + with open(path, 'w') as f: + json.dump(content.freeze(), f, indent=2, sort_keys=True) + + def get_pec_excitation_rate(element, charge, transition, repository_path=None): """ Reads the excitation PEC from the repository for the given @@ -252,13 +341,35 @@ def get_pec_recombination_rate(element, charge, transition, repository_path=None return _get_pec_rate('recombination', element, charge, transition, repository_path) -def get_pec_thermalcx_rate(element, charge, transition, repository_path=None): +def _get_pec_rate(cls, element, charge, transition, repository_path=None): + + repository_path = repository_path or DEFAULT_REPOSITORY_PATH + path = os.path.join(repository_path, 'pec/{}/{}/{}.json'.format(cls, element.symbol.lower(), charge)) + try: + with open(path, 'r') as f: + content = json.load(f) + d = content[encode_transition(transition)] + except (FileNotFoundError, KeyError): + raise RuntimeError('Requested PEC rate (class={}, element={}, charge={}, transition={})' + ' is not available.'.format(cls, element.symbol, charge, transition)) + + # convert to numpy arrays + d['ne'] = np.array(d['ne'], np.float64) + d['te'] = np.array(d['te'], np.float64) + d['rate'] = np.array(d['rate'], np.float64) + + return d + + +def get_pec_thermal_cx_rate(donor_element, donor_charge, receiver_element, receiver_charge, transition, repository_path=None): """ Reads the thermal charge exchange PEC from the repository for the given - element, charge and transition. + donor element, donor charge, receiver element, receiver charge and transition. - :param element: Plasma species (Element/Isotope). - :param charge: Charge of the plasma species. + :param donor_element: Electron donor plasma species (Element/Isotope). + :param donor_charge: Electron donor charge. + :param receiver_element: Electron receiver plasma species (Element/Isotope). + :param receiver_charge: Electron receiver charge. :param transition: Tuple containing (initial level, final level). :param repository_path: Path to the atomic data repository. @@ -266,28 +377,29 @@ def get_pec_thermalcx_rate(element, charge, transition, repository_path=None): | 'ne': 1D array of size (N) with electron density in m^-3, | 'te': 1D array of size (M) with electron temperature in eV, - | 'rate': 2D array of size (N, M) with thermal CX PEC in photon.m^3.s^-1. + | 'td': 1D array of size (K) with donor temperature in eV, + | 'rate': 2D array of size (N, M, K) with thermal CX PEC in photon.m^3.s^-1. """ - return _get_pec_rate('thermalcx', element, charge, transition, repository_path) - - -def _get_pec_rate(cls, element, charge, transition, repository_path=None): repository_path = repository_path or DEFAULT_REPOSITORY_PATH - path = os.path.join(repository_path, 'pec/{}/{}/{}.json'.format(cls, element.symbol.lower(), charge)) + + rate_path = 'pec/thermal_cx/{0}/{1}/{2}/{3}.json'.format(donor_element.symbol.lower(), donor_charge, + receiver_element.symbol.lower(), receiver_charge) + path = os.path.join(repository_path, rate_path) try: with open(path, 'r') as f: content = json.load(f) d = content[encode_transition(transition)] except (FileNotFoundError, KeyError): - raise RuntimeError('Requested PEC rate (class={}, element={}, charge={}, transition={})' - ' is not available.'.format(cls, element.symbol, charge, transition)) + raise RuntimeError('Requested thermal charge-exchange PEC (donor={}, donor charge={}, receiver={}, receiver charge={})' + ' is not available.' + ''.format(donor_element.symbol, donor_charge, receiver_element.symbol, receiver_charge)) # convert to numpy arrays d['ne'] = np.array(d['ne'], np.float64) d['te'] = np.array(d['te'], np.float64) + d['td'] = np.array(d['td'], np.float64) d['rate'] = np.array(d['rate'], np.float64) return d - diff --git a/demos/emission_models/charge_exchange.py b/demos/emission_models/charge_exchange.py index 11d95386..02d61235 100644 --- a/demos/emission_models/charge_exchange.py +++ b/demos/emission_models/charge_exchange.py @@ -122,7 +122,7 @@ integration_step = 0.0025 beam_transform = translate(-0.5, 0.0, 0) * rotate_basis(Vector3D(1, 0, 0), Vector3D(0, 0, 1)) -beam_energy = 50000 # keV +beam_energy = 50000 # eV beam_full = Beam(parent=world, transform=beam_transform) beam_full.plasma = plasma diff --git a/demos/emission_models/thermal_charge_exchange.py b/demos/emission_models/thermal_charge_exchange.py new file mode 100644 index 00000000..42c7b88a --- /dev/null +++ b/demos/emission_models/thermal_charge_exchange.py @@ -0,0 +1,126 @@ + +import matplotlib.pyplot as plt +from scipy.constants import electron_mass, atomic_mass + +from raysect.core import translate, rotate_basis, Point3D, Vector3D +from raysect.primitive import Box +from raysect.optical import World, Ray + +from cherab.core import Plasma, Beam, Maxwellian, Species +from cherab.core.math import ScalarToVectorFunction3D +from cherab.core.atomic import hydrogen, deuterium, carbon, Line +from cherab.core.model import SingleRayAttenuator, BeamCXLine, ThermalCXLine +from cherab.tools.plasmas.slab import NeutralFunction, IonFunction +from cherab.openadas import OpenADAS +from cherab.openadas.install import install_adf15 + + +############### +# Make Plasma # + +width = 1.0 +length = 1.0 +height = 3.0 +peak_density = 5e19 +edge_density = 1e18 +pedestal_top = 1 +neutral_temperature = 0.5 +peak_temperature = 2500 +edge_temperature = 25 +impurities = [(carbon, 6, 0.005)] + +world = World() +adas = OpenADAS(permit_extrapolation=True, missing_rates_return_null=True) + +plasma = Plasma(parent=world) +plasma.atomic_data = adas + +# plasma slab along x direction +plasma.geometry = Box(Point3D(0, -width / 2, -height / 2), Point3D(length, width / 2, height / 2)) + +species = [] + +# make a non-zero velocity profile for the plasma +vy_profile = IonFunction(1E5, 0, pedestal_top=pedestal_top) +velocity_profile = ScalarToVectorFunction3D(0, vy_profile, 0) + +# define neutral species distribution +h0_density = NeutralFunction(peak_density, 0.1, pedestal_top=pedestal_top) +h0_temperature = neutral_temperature +h0_distribution = Maxwellian(h0_density, h0_temperature, velocity_profile, + hydrogen.atomic_weight * atomic_mass) +species.append(Species(hydrogen, 0, h0_distribution)) + +# define hydrogen ion species distribution +h1_density = IonFunction(peak_density, edge_density, pedestal_top=pedestal_top) +h1_temperature = IonFunction(peak_temperature, edge_temperature, pedestal_top=pedestal_top) +h1_distribution = Maxwellian(h1_density, h1_temperature, velocity_profile, + hydrogen.atomic_weight * atomic_mass) +species.append(Species(hydrogen, 1, h1_distribution)) + +# add impurities +if impurities: + for impurity, ionisation, concentration in impurities: + imp_density = IonFunction(peak_density * concentration, edge_density * concentration, pedestal_top=pedestal_top) + imp_temperature = IonFunction(peak_temperature, edge_temperature, pedestal_top=pedestal_top) + imp_distribution = Maxwellian(imp_density, imp_temperature, velocity_profile, + impurity.atomic_weight * atomic_mass) + species.append(Species(impurity, ionisation, imp_distribution)) + +# define the electron distribution +e_density = IonFunction(peak_density, edge_density, pedestal_top=pedestal_top) +e_temperature = IonFunction(peak_temperature, edge_temperature, pedestal_top=pedestal_top) +e_distribution = Maxwellian(e_density, e_temperature, velocity_profile, electron_mass) + +# define species +plasma.electron_distribution = e_distribution +plasma.composition = species + +# add thermal CX emission model +cVI_5_4 = Line(carbon, 5, (5, 4)) +plasma.models = [ThermalCXLine(cVI_5_4)] + +# trace thermal CX spectrum along y direction +ray = Ray(origin=Point3D(0.4, -3.5, 0), direction=Vector3D(0, 1, 0), + min_wavelength=112.3, max_wavelength=112.7, bins=512) +thermal_cx_spectrum = ray.trace(world) + +########################### +# Inject beam into plasma # + +integration_step = 0.0025 +# injected along x direction +beam_transform = translate(-0.5, 0.0, 0) * rotate_basis(Vector3D(1, 0, 0), Vector3D(0, 0, 1)) +beam_energy = 50000 # eV + +beam = Beam(parent=world, transform=beam_transform) +beam.plasma = plasma +beam.atomic_data = adas +beam.energy = beam_energy +beam.power = 3e6 +beam.element = deuterium +beam.sigma = 0.05 +beam.divergence_x = 0.5 +beam.divergence_y = 0.5 +beam.length = 3.0 +beam.attenuator = SingleRayAttenuator(clamp_to_zero=True) +beam.integrator.step = integration_step +beam.integrator.min_samples = 10 + +# remove thermal CX model and add beam CX model +plasma.models = [] +beam.models = [BeamCXLine(cVI_5_4)] + +# trace the spectrum again +beam_cx_spectrum = ray.trace(world) + +# plot the spectra +plt.figure() +plt.plot(thermal_cx_spectrum.wavelengths, thermal_cx_spectrum.samples, label='thermal CX') +plt.plot(beam_cx_spectrum.wavelengths, beam_cx_spectrum.samples, label='beam CX') +plt.legend() +plt.xlabel('Wavelength (nm)') +plt.ylabel('Radiance (W/m^2/str/nm)') +plt.title('Sampled CX spectra') + +plt.show() diff --git a/docs/source/models/basic_line/basic_line_emission.rst b/docs/source/models/basic_line/basic_line_emission.rst index f8c80c60..b71fc03c 100644 --- a/docs/source/models/basic_line/basic_line_emission.rst +++ b/docs/source/models/basic_line/basic_line_emission.rst @@ -5,3 +5,5 @@ Basic Line Emission .. autoclass:: cherab.core.model.plasma.impact_excitation.ExcitationLine .. autoclass:: cherab.core.model.plasma.recombination.RecombinationLine + +.. autoclass:: cherab.core.model.plasma.thermal_cx.ThermalCXLine