diff --git a/CHANGELOG.md b/CHANGELOG.md index 9af7646d..de0692ab 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,8 @@ New: * Add custom line shape support to BeamCXLine model. (#394) * Add PeriodicTransformXD and VectorPeriodicTransformXD functions to support the data simulated with periodic boundary conditions. (#387) * Add CylindricalTransform and VectorCylindricalTransform to transform functions from cylindrical to Cartesian coordinates. (#387) +* Add numerical integration of Bremsstrahlung spectrum over a spectral bin. (#395) +* 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) diff --git a/cherab/core/model/plasma/bremsstrahlung.pxd b/cherab/core/model/plasma/bremsstrahlung.pxd index f8251b28..44025c8c 100644 --- a/cherab/core/model/plasma/bremsstrahlung.pxd +++ b/cherab/core/model/plasma/bremsstrahlung.pxd @@ -19,20 +19,27 @@ # cython: language_level=3 from numpy cimport ndarray +from cherab.core.math cimport Function1D +from cherab.core.math.integrators cimport Integrator1D from cherab.core.atomic cimport FreeFreeGauntFactor from cherab.core.plasma cimport PlasmaModel -cdef class Bremsstrahlung(PlasmaModel): +cdef class BremsFunction(Function1D): cdef: - FreeFreeGauntFactor _gaunt_factor - bint _user_provided_gaunt_factor - ndarray _species_charge, _species_density - double[::1] _species_density_mv, _species_charge_mv + double ne, te + FreeFreeGauntFactor gaunt_factor + ndarray species_density, species_charge + double[::1] species_density_mv + double[::1] species_charge_mv - cdef double _bremsstrahlung(self, double wvl, double te, double ne) - cdef int _populate_cache(self) except -1 +cdef class Bremsstrahlung(PlasmaModel): + cdef: + BremsFunction _brems_func + bint _user_provided_gaunt_factor + Integrator1D _integrator + cdef int _populate_cache(self) except -1 diff --git a/cherab/core/model/plasma/bremsstrahlung.pyx b/cherab/core/model/plasma/bremsstrahlung.pyx index 2ba7245b..404750a2 100644 --- a/cherab/core/model/plasma/bremsstrahlung.pyx +++ b/cherab/core/model/plasma/bremsstrahlung.pyx @@ -21,16 +21,73 @@ import numpy as np from raysect.optical cimport Spectrum, Point3D, Vector3D from cherab.core cimport Plasma, AtomicData -from cherab.core.atomic cimport FreeFreeGauntFactor +from cherab.core.math.integrators cimport GaussianQuadrature from cherab.core.species cimport Species -from cherab.core.utility.constants cimport RECIP_4_PI, ELEMENTARY_CHARGE, SPEED_OF_LIGHT, PLANCK_CONSTANT -from libc.math cimport sqrt, log, exp +from cherab.core.utility.constants cimport RECIP_4_PI, ELEMENTARY_CHARGE, SPEED_OF_LIGHT, PLANCK_CONSTANT, ELECTRON_REST_MASS, VACUUM_PERMITTIVITY +from libc.math cimport sqrt, log, exp, M_PI cimport cython -cdef double PH_TO_J_FACTOR = PLANCK_CONSTANT * SPEED_OF_LIGHT * 1e9 +cdef double EXP_FACTOR = PLANCK_CONSTANT * SPEED_OF_LIGHT * 1e9 / ELEMENTARY_CHARGE -cdef double EXP_FACTOR = PH_TO_J_FACTOR / ELEMENTARY_CHARGE +cdef double BREMS_CONST = (ELEMENTARY_CHARGE**2 * RECIP_4_PI / VACUUM_PERMITTIVITY)**3 +BREMS_CONST *= 32 * M_PI**2 / (3 * sqrt(3) * ELECTRON_REST_MASS**2 * SPEED_OF_LIGHT**3) +BREMS_CONST *= sqrt(2 * ELECTRON_REST_MASS / (M_PI * ELEMENTARY_CHARGE)) +BREMS_CONST *= SPEED_OF_LIGHT * 1e9 * RECIP_4_PI + + +cdef class BremsFunction(Function1D): + """ + Calculates bremsstrahlung spectrum. + + :param FreeFreeGauntFactor gaunt_factor: Free-free Gaunt factor as a function of Z, Te and + wavelength. + :param object species_density: Array-like object wiyh ions' density in m-3. + :param object species_charge: Array-like object wiyh ions' charge. + :param double ne: Electron density in m-3. + :param double te: Electron temperature in eV. + """ + + def __init__(self, FreeFreeGauntFactor gaunt_factor, object species_density, object species_charge, double ne, double te): + + if ne <= 0: + raise ValueError("Argument ne must be positive.") + self.ne = ne + + if te <= 0: + raise ValueError("Argument te must be positive.") + self.te = te + + self.gaunt_factor = gaunt_factor + self.species_density = np.asarray(species_density, dtype=np.float64) # copied if type does not match + self.species_density_mv = self.species_density + self.species_charge = np.asarray(species_charge, dtype=np.float64) # copied if type does not match + self.species_charge_mv = self.species_charge + + @cython.cdivision(True) + @cython.boundscheck(False) + @cython.wraparound(False) + cdef double evaluate(self, double wvl) except? -1e999: + """ + :param double wvl: Wavelength in nm. + :return: + """ + + cdef double ni_gff_z2, radiance, pre_factor, ni, z + cdef int i + + ni_gff_z2 = 0 + for i in range(self.species_charge_mv.shape[0]): + z = self.species_charge_mv[i] + ni = self.species_density_mv[i] + if ni > 0: + ni_gff_z2 += ni * self.gaunt_factor.evaluate(z, self.te, wvl) * z * z + + # bremsstrahlung equation W/m^3/str/nm + pre_factor = BREMS_CONST / (sqrt(self.te) * wvl * wvl) * self.ne * ni_gff_z2 + radiance = pre_factor * exp(- EXP_FACTOR / (self.te * wvl)) + + return radiance # todo: doppler shift? @@ -38,26 +95,44 @@ cdef class Bremsstrahlung(PlasmaModel): """ Emitter that calculates bremsstrahlung emission from a plasma object. - The bremmstrahlung formula implemented is equation 2 from M. Beurskens, - et. al., 'ITER LIDAR performance analysis', Rev. Sci. Instrum. 79, 10E727 (2008), + The bremmstrahlung formula implemented is equation 5.3.40 + from I. H. Hutchinson, 'Principles of Plasma Diagnostics', second edition, + Cambridge University Press, 2002, ISBN: 9780511613630, + https://doi.org/10.1017/CBO9780511613630 + + Note that in eq. 5.3.40, the emissivity :math:`j(\\nu)` is given in (W/m^3/sr/Hz) with respect + to frequency, :math:`\\nu`. Here, the emissivity :math:`\\epsilon_{\\mathrm{ff}}(\\lambda)` + is given in (W/m^3/nm/sr) with respect to wavelength, :math:`\\lambda = \\frac{10^{9} c}{\\nu}`, + and taking into account that :math:`d\\nu=-\\frac{10^{9} c}{\\lambda^2}d\\lambda`. .. math:: - \\epsilon (\\lambda) = \\frac{0.95 \\times 10^{-19}}{\\lambda 4 \\pi} \\sum_{i} \\left(g_{ff}(Z_i, T_e, \\lambda) n_i Z_i^2\\right) n_e T_e^{-1/2} \\times \\exp{\\frac{-hc}{\\lambda T_e}}, + \\epsilon_{\\mathrm{ff}}(\\lambda) = \\left( \\frac{e^2}{4 \\pi \\varepsilon_0} \\right)^3 + \\frac{32 \\pi^2}{3 \\sqrt{3} m_\\mathrm{e}^2 c^3} + \\sqrt{\\frac{2 m_\\mathrm{e}^3}{\\pi e T_\\mathrm{e}}} + \\frac{10^{9} c}{4 \\pi \\lambda^2} + n_\\mathrm{e} \\sum_i \\left( n_\\mathrm{i} g_\\mathrm{ff} (Z_\\mathrm{i}, T_\\mathrm{e}, \\lambda) Z_\\mathrm{i}^2 \\right) + \\mathrm{e}^{-\\frac{10^9 hc}{e T_\\mathrm{e} \\lambda}}\\,, - where the emission :math:`\\epsilon (\\lambda)` is in units of radiance (ph/s/sr/m^3/nm). + where :math:`T_\\mathrm{e}` is in eV and :math:`\\lambda` is in nm. + + :math:`g_\\mathrm{ff} (Z_\\mathrm{i}, T_\\mathrm{e}, \\lambda)` is the free-free Gaunt factor. :ivar Plasma plasma: The plasma to which this emission model is attached. Default is None. :ivar AtomicData atomic_data: The atomic data provider for this model. Default is None. :ivar FreeFreeGauntFactor gaunt_factor: Free-free Gaunt factor as a function of Z, Te and wavelength. If not provided, the `atomic_data` is used. + :ivar Integrator1D integrator: Integrator1D instance to integrate Bremsstrahlung radiation + over the spectral bin. Default is `GaussianQuadrature`. """ - def __init__(self, Plasma plasma=None, AtomicData atomic_data=None, FreeFreeGauntFactor gaunt_factor=None): + def __init__(self, Plasma plasma=None, AtomicData atomic_data=None, FreeFreeGauntFactor gaunt_factor=None, Integrator1D integrator=None): super().__init__(plasma, atomic_data) + self._brems_func = BremsFunction.__new__(BremsFunction) self.gaunt_factor = gaunt_factor + self.integrator = integrator or GaussianQuadrature() # ensure that cache is initialised self._change() @@ -65,14 +140,25 @@ cdef class Bremsstrahlung(PlasmaModel): @property def gaunt_factor(self): - return self._gaunt_factor + return self._brems_func.gaunt_factor @gaunt_factor.setter def gaunt_factor(self, value): - self._gaunt_factor = value + self._brems_func.gaunt_factor = value self._user_provided_gaunt_factor = True if value else False + @property + def integrator(self): + + return self._integrator + + @integrator.setter + def integrator(self, Integrator1D value not None): + + self._integrator = value + self._integrator.function = self._brems_func + def __repr__(self): return '' @@ -84,13 +170,12 @@ cdef class Bremsstrahlung(PlasmaModel): cdef: double ne, te - double lower_wavelength, upper_wavelength - double lower_sample, upper_sample + double lower_wavelength, upper_wavelength, bin_integral Species species int i # cache data on first run - if self._species_charge is None: + if self._brems_func.species_charge is None: self._populate_cache() ne = self._plasma.get_electron_distribution().density(point.x, point.y, point.z) @@ -100,57 +185,28 @@ cdef class Bremsstrahlung(PlasmaModel): if te <= 0: return spectrum + self._brems_func.ne = ne + self._brems_func.te = te + # collect densities of charged species i = 0 for species in self._plasma.get_composition(): if species.charge > 0: - self._species_density_mv[i] = species.distribution.density(point.x, point.y, point.z) + self._brems_func.species_density_mv[i] = species.distribution.density(point.x, point.y, point.z) i += 1 - # numerically integrate using trapezium rule - # todo: add sub-sampling to increase numerical accuracy + # add bremsstrahlung to spectrum lower_wavelength = spectrum.min_wavelength - lower_sample = self._bremsstrahlung(lower_wavelength, te, ne) for i in range(spectrum.bins): - upper_wavelength = spectrum.min_wavelength + spectrum.delta_wavelength * (i + 1) - upper_sample = self._bremsstrahlung(upper_wavelength, te, ne) - spectrum.samples_mv[i] += 0.5 * (lower_sample + upper_sample) + bin_integral = self._integrator.evaluate(lower_wavelength, upper_wavelength) + spectrum.samples_mv[i] += bin_integral / spectrum.delta_wavelength lower_wavelength = upper_wavelength - lower_sample = upper_sample return spectrum - @cython.cdivision(True) - @cython.boundscheck(False) - @cython.wraparound(False) - cdef double _bremsstrahlung(self, double wvl, double te, double ne): - """ - :param double wvl: Wavelength in nm. - :param double te: Electron temperature in eV - :param double ne: Electron dencity in m^-3 - :return: - """ - - cdef double ni_gff_z2, radiance, pre_factor, ni, z - cdef int i - - ni_gff_z2 = 0 - for i in range(self._species_charge_mv.shape[0]): - z = self._species_charge_mv[i] - ni = self._species_density_mv[i] - if ni > 0: - ni_gff_z2 += ni * self._gaunt_factor.evaluate(z, te, wvl) * z * z - - # bremsstrahlung equation W/m^3/str/nm - pre_factor = 0.95e-19 * RECIP_4_PI * ni_gff_z2 * ne / (sqrt(te) * wvl) - radiance = pre_factor * exp(- EXP_FACTOR / (te * wvl)) * PH_TO_J_FACTOR - - # convert to W/m^3/str/nm - return radiance / wvl - cdef int _populate_cache(self) except -1: cdef list species_charge @@ -159,12 +215,12 @@ cdef class Bremsstrahlung(PlasmaModel): if self._plasma is None: raise RuntimeError("The emission model is not connected to a plasma object.") - if self._gaunt_factor is None: + if self._brems_func.gaunt_factor is None: if self._atomic_data is None: raise RuntimeError("The emission model is not connected to an atomic data source.") # initialise Gaunt factor on first run using the atomic data - self._gaunt_factor = self._atomic_data.free_free_gaunt_factor() + self._brems_func.gaunt_factor = self._atomic_data.free_free_gaunt_factor() species_charge = [] for species in self._plasma.get_composition(): @@ -172,18 +228,18 @@ cdef class Bremsstrahlung(PlasmaModel): species_charge.append(species.charge) # Gaunt factor takes Z as double to support Zeff, so caching Z as float64 - self._species_charge = np.array(species_charge, dtype=np.float64) - self._species_charge_mv = self._species_charge + self._brems_func.species_charge = np.array(species_charge, dtype=np.float64) + self._brems_func.species_charge_mv = self._brems_func.species_charge - self._species_density = np.zeros_like(self._species_charge) - self._species_density_mv = self._species_density + self._brems_func.species_density = np.zeros_like(self._brems_func.species_charge) + self._brems_func.species_density_mv = self._brems_func.species_density def _change(self): # clear cache to force regeneration on first use if not self._user_provided_gaunt_factor: - self._gaunt_factor = None - self._species_charge = None - self._species_charge_mv = None - self._species_density = None - self._species_density_mv = None + self._brems_func.gaunt_factor = None + self._brems_func.species_charge = None + self._brems_func.species_charge_mv = None + self._brems_func.species_density = None + self._brems_func.species_density_mv = None diff --git a/cherab/core/tests/test_bremsstrahlung.py b/cherab/core/tests/test_bremsstrahlung.py new file mode 100644 index 00000000..5373776b --- /dev/null +++ b/cherab/core/tests/test_bremsstrahlung.py @@ -0,0 +1,94 @@ +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +import unittest + +import numpy as np + +from raysect.core import Point3D, Vector3D +from raysect.optical import World, Ray + +from cherab.core.atomic import AtomicData, MaxwellianFreeFreeGauntFactor +from cherab.core.math.integrators import GaussianQuadrature +from cherab.core.atomic import deuterium, nitrogen +from cherab.tools.plasmas.slab import build_constant_slab_plasma +from cherab.core.model import Bremsstrahlung + +import scipy.constants as const + + +class TestBremsstrahlung(unittest.TestCase): + + world = World() + + plasma_species = [(deuterium, 1, 1.e19, 2000., Vector3D(0, 0, 0)), (nitrogen, 7, 1.e18, 2000., Vector3D(0, 0, 0))] + plasma = build_constant_slab_plasma(length=1, width=1, height=1, electron_density=1e19, electron_temperature=2000., + plasma_species=plasma_species) + plasma.parent = world + plasma.atomic_data = AtomicData() + + def test_bremsstrahlung_model(self): + # setting up the model + gaunt_factor = MaxwellianFreeFreeGauntFactor() + bremsstrahlung = Bremsstrahlung(gaunt_factor=gaunt_factor) + self.plasma.models = [bremsstrahlung] + + # observing + origin = Point3D(1.5, 0, 0) + direction = Vector3D(-1, 0, 0) + ray = Ray(origin=origin, direction=direction, + min_wavelength=400., max_wavelength=800., bins=128) + brems_spectrum = ray.trace(self.world) + + # validating + brems_const = (const.e**2 * 0.25 / np.pi / const.epsilon_0)**3 + brems_const *= 32 * np.pi**2 / (3 * np.sqrt(3) * const.m_e**2 * const.c**3) + brems_const *= np.sqrt(2 * const.m_e / (np.pi * const.e)) + brems_const *= const.c * 1e9 * 0.25 / np.pi + exp_factor = const.h * const.c * 1.e9 / const. e + + ne = self.plasma.electron_distribution.density(0.5, 0, 0) + te = self.plasma.electron_distribution.effective_temperature(0.5, 0, 0) + + def brems_func(wvl): + ni_gff_z2 = 0 + for species in self.plasma.composition: + z = species.charge + ni = self.plasma.composition[(species.element, species.charge)].distribution.density(0.5, 0, 0) + ni_gff_z2 += ni * gaunt_factor(z, te, wvl) * z * z + + return brems_const * ni_gff_z2 * ne / (np.sqrt(te) * wvl * wvl) * np.exp(- exp_factor / (te * wvl)) + + integrator = GaussianQuadrature(brems_func) + + test_samples = np.zeros(brems_spectrum.bins) + delta_wavelength = (brems_spectrum.max_wavelength - brems_spectrum.min_wavelength) / brems_spectrum.bins + lower_wavelength = brems_spectrum.min_wavelength + for i in range(brems_spectrum.bins): + upper_wavelength = brems_spectrum.min_wavelength + delta_wavelength * (i + 1) + bin_integral = integrator(lower_wavelength, upper_wavelength) + test_samples[i] = bin_integral / delta_wavelength + lower_wavelength = upper_wavelength + + for i in range(brems_spectrum.bins): + self.assertAlmostEqual(brems_spectrum.samples[i], test_samples[i], delta=1e-10, + msg='BeamCXLine model gives a wrong value at {} nm.'.format(brems_spectrum.wavelengths[i])) + + +if __name__ == '__main__': + unittest.main() diff --git a/cherab/core/utility/constants.pxd b/cherab/core/utility/constants.pxd index 9cce304d..4c59b4a8 100644 --- a/cherab/core/utility/constants.pxd +++ b/cherab/core/utility/constants.pxd @@ -30,3 +30,4 @@ cdef: double ELECTRON_CLASSICAL_RADIUS double ELECTRON_REST_MASS double RYDBERG_CONSTANT_EV + double VACUUM_PERMITTIVITY diff --git a/cherab/core/utility/constants.pyx b/cherab/core/utility/constants.pyx index 03e70617..d36b99f9 100644 --- a/cherab/core/utility/constants.pyx +++ b/cherab/core/utility/constants.pyx @@ -32,3 +32,4 @@ cdef: double ELECTRON_CLASSICAL_RADIUS = 2.8179403262e-15 double ELECTRON_REST_MASS = 9.1093837015e-31 double RYDBERG_CONSTANT_EV = 13.605693122994 + double VACUUM_PERMITTIVITY = 8.8541878128e-12