From 394fe4cb205c23a5958b58d61164489c2f85f76a Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Tue, 24 Jan 2023 22:13:12 +0300 Subject: [PATCH 1/8] Add numerical integration to Bremsstrahlung plasma model. --- cherab/core/model/plasma/bremsstrahlung.pxd | 21 ++- cherab/core/model/plasma/bremsstrahlung.pyx | 149 +++++++++++++------- 2 files changed, 109 insertions(+), 61 deletions(-) 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..e076478c 100644 --- a/cherab/core/model/plasma/bremsstrahlung.pyx +++ b/cherab/core/model/plasma/bremsstrahlung.pyx @@ -21,7 +21,7 @@ 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 @@ -33,6 +33,61 @@ cdef double PH_TO_J_FACTOR = PLANCK_CONSTANT * SPEED_OF_LIGHT * 1e9 cdef double EXP_FACTOR = PH_TO_J_FACTOR / ELEMENTARY_CHARGE +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 = 0.95e-19 * RECIP_4_PI * ni_gff_z2 * self.ne / (sqrt(self.te) * wvl) + radiance = pre_factor * exp(- EXP_FACTOR / (self.te * wvl)) * PH_TO_J_FACTOR + + # convert to W/m^3/str/nm + return radiance / wvl + + # todo: doppler shift? cdef class Bremsstrahlung(PlasmaModel): """ @@ -51,13 +106,17 @@ cdef class Bremsstrahlung(PlasmaModel): :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=GaussianQuadrature()): super().__init__(plasma, atomic_data) + self._brems_func = BremsFunction.__new__(BremsFunction) self.gaunt_factor = gaunt_factor + self.integrator = integrator # ensure that cache is initialised self._change() @@ -65,14 +124,24 @@ 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 + def __repr__(self): return '' @@ -84,13 +153,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 +168,30 @@ 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 + self._integrator.function = self._brems_func + + # 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 +200,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 +213,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 From 0578c09dfd90b02a88a97b47baac58e088a6362d Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Tue, 24 Jan 2023 23:08:04 +0300 Subject: [PATCH 2/8] Update changelog. --- CHANGELOG.md | 7 +++++++ cherab/core/model/plasma/bremsstrahlung.pyx | 4 ++-- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d7bf7c6d..276040eb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,13 @@ Project Changelog ================= +Release 1.5.0 (TBD) +------------------- + +New: +* Add numerical integration of Bremsstrahlung spectrum over a spectral bin. (#395) + + Release 1.4.0 (TBD) ------------------- diff --git a/cherab/core/model/plasma/bremsstrahlung.pyx b/cherab/core/model/plasma/bremsstrahlung.pyx index e076478c..01e896a9 100644 --- a/cherab/core/model/plasma/bremsstrahlung.pyx +++ b/cherab/core/model/plasma/bremsstrahlung.pyx @@ -110,13 +110,13 @@ cdef class Bremsstrahlung(PlasmaModel): over the spectral bin. Default is `GaussianQuadrature`. """ - def __init__(self, Plasma plasma=None, AtomicData atomic_data=None, FreeFreeGauntFactor gaunt_factor=None, Integrator1D integrator=GaussianQuadrature()): + 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 + self.integrator = integrator or GaussianQuadrature() # ensure that cache is initialised self._change() From bbcc03092350d3a3a316dfc7e30e70fe9ca61c09 Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Thu, 26 Jan 2023 12:54:53 +0300 Subject: [PATCH 3/8] Set the integrand in Bremsstrahlung in integrator setter and not in emission(). --- cherab/core/model/plasma/bremsstrahlung.pyx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/cherab/core/model/plasma/bremsstrahlung.pyx b/cherab/core/model/plasma/bremsstrahlung.pyx index 01e896a9..fa761eae 100644 --- a/cherab/core/model/plasma/bremsstrahlung.pyx +++ b/cherab/core/model/plasma/bremsstrahlung.pyx @@ -141,6 +141,7 @@ cdef class Bremsstrahlung(PlasmaModel): def integrator(self, Integrator1D value not None): self._integrator = value + self._integrator.function = self._brems_func def __repr__(self): return '' @@ -178,8 +179,6 @@ cdef class Bremsstrahlung(PlasmaModel): self._brems_func.species_density_mv[i] = species.distribution.density(point.x, point.y, point.z) i += 1 - self._integrator.function = self._brems_func - # add bremsstrahlung to spectrum lower_wavelength = spectrum.min_wavelength for i in range(spectrum.bins): From e4305474e8310d199f5e41f9ac592d2a7a56fd1f Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Mon, 1 May 2023 16:41:38 +0300 Subject: [PATCH 4/8] The numerical constant in the bremsstrahlung model is expressed in terms of physical constants as suggested by @skuba31. --- cherab/core/model/plasma/bremsstrahlung.pyx | 35 ++++++++++++++------- cherab/core/utility/constants.pxd | 1 + cherab/core/utility/constants.pyx | 1 + 3 files changed, 25 insertions(+), 12 deletions(-) diff --git a/cherab/core/model/plasma/bremsstrahlung.pyx b/cherab/core/model/plasma/bremsstrahlung.pyx index fa761eae..16b10aef 100644 --- a/cherab/core/model/plasma/bremsstrahlung.pyx +++ b/cherab/core/model/plasma/bremsstrahlung.pyx @@ -23,14 +23,17 @@ from raysect.optical cimport Spectrum, Point3D, Vector3D from cherab.core cimport Plasma, AtomicData 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): @@ -80,9 +83,9 @@ cdef class BremsFunction(Function1D): 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 = 0.95e-19 * RECIP_4_PI * ni_gff_z2 * self.ne / (sqrt(self.te) * wvl) - radiance = pre_factor * exp(- EXP_FACTOR / (self.te * wvl)) * PH_TO_J_FACTOR + # bremsstrahlung equation W/m^3/str + pre_factor = BREMS_CONST * ni_gff_z2 * self.ne / (sqrt(self.te) * wvl) + radiance = pre_factor * exp(- EXP_FACTOR / (self.te * wvl)) # convert to W/m^3/str/nm return radiance / wvl @@ -93,13 +96,21 @@ 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 .. 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}}, - - where the emission :math:`\\epsilon (\\lambda)` is in units of radiance (ph/s/sr/m^3/nm). + \\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}}}\\,, + + where the emission :math:`\\epsilon_{\\mathrm{ff}}(\\lambda)` is in units of radiance (W/sr/m^3/nm), + :math:`T_\\mathrm{e}` is in eV and :math:`\\lambda` is in nm. :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. 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 From 5f61aaec157a5f88b4b41c78636e7242527b45dc Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Mon, 1 May 2023 18:49:20 +0300 Subject: [PATCH 5/8] Updated bremsstrahlung model docstring. --- cherab/core/model/plasma/bremsstrahlung.pyx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cherab/core/model/plasma/bremsstrahlung.pyx b/cherab/core/model/plasma/bremsstrahlung.pyx index 16b10aef..cb8a88a6 100644 --- a/cherab/core/model/plasma/bremsstrahlung.pyx +++ b/cherab/core/model/plasma/bremsstrahlung.pyx @@ -112,6 +112,8 @@ cdef class Bremsstrahlung(PlasmaModel): where the emission :math:`\\epsilon_{\\mathrm{ff}}(\\lambda)` is in units of radiance (W/sr/m^3/nm), :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 From 00c70f30303eff0a3d42f05dcd6a36b3d45c4509 Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Tue, 2 May 2023 01:52:47 +0300 Subject: [PATCH 6/8] Added a test for Bremsstrahlung model. --- cherab/core/tests/test_bremsstrahlung.py | 94 ++++++++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100644 cherab/core/tests/test_bremsstrahlung.py 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() From 8e719ba31aed225e19b963b191a722d56452391a Mon Sep 17 00:00:00 2001 From: vsnever Date: Wed, 3 May 2023 13:28:36 +0300 Subject: [PATCH 7/8] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7f982618..19f1ae56 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ New: * 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) Release 1.4.0 (3 Feb 2023) From b3b0d5713a73eaa8e53a6be4d0a6fb605f135d2e Mon Sep 17 00:00:00 2001 From: vsnever Date: Wed, 3 May 2023 23:51:47 +0300 Subject: [PATCH 8/8] Made expression for bremsstrahlung look clearer and updated docstring in response to @skuba31 comments. --- cherab/core/model/plasma/bremsstrahlung.pyx | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/cherab/core/model/plasma/bremsstrahlung.pyx b/cherab/core/model/plasma/bremsstrahlung.pyx index cb8a88a6..404750a2 100644 --- a/cherab/core/model/plasma/bremsstrahlung.pyx +++ b/cherab/core/model/plasma/bremsstrahlung.pyx @@ -83,12 +83,11 @@ cdef class BremsFunction(Function1D): if ni > 0: ni_gff_z2 += ni * self.gaunt_factor.evaluate(z, self.te, wvl) * z * z - # bremsstrahlung equation W/m^3/str - pre_factor = BREMS_CONST * ni_gff_z2 * self.ne / (sqrt(self.te) * wvl) + # 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)) - # convert to W/m^3/str/nm - return radiance / wvl + return radiance # todo: doppler shift? @@ -101,16 +100,20 @@ cdef class Bremsstrahlung(PlasmaModel): 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_{\\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}}}\\,, + \\mathrm{e}^{-\\frac{10^9 hc}{e T_\\mathrm{e} \\lambda}}\\,, - where the emission :math:`\\epsilon_{\\mathrm{ff}}(\\lambda)` is in units of radiance (W/sr/m^3/nm), - :math:`T_\\mathrm{e}` is in eV and :math:`\\lambda` is in 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.