From b202dbb2ac2762c0f9259044e98dab8b5b1b7754 Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Tue, 30 Aug 2022 20:35:12 +0300 Subject: [PATCH 1/9] Add Gaussian quadrature integrator and tests. --- cherab/core/math/integrators/__init__.pxd | 20 ++ cherab/core/math/integrators/__init__.py | 19 ++ .../core/math/integrators/integrators1d.pxd | 36 ++++ .../core/math/integrators/integrators1d.pyx | 199 ++++++++++++++++++ cherab/core/math/tests/test_integrators.py | 78 +++++++ 5 files changed, 352 insertions(+) create mode 100644 cherab/core/math/integrators/__init__.pxd create mode 100644 cherab/core/math/integrators/__init__.py create mode 100644 cherab/core/math/integrators/integrators1d.pxd create mode 100644 cherab/core/math/integrators/integrators1d.pyx create mode 100644 cherab/core/math/tests/test_integrators.py diff --git a/cherab/core/math/integrators/__init__.pxd b/cherab/core/math/integrators/__init__.pxd new file mode 100644 index 00000000..db06ef43 --- /dev/null +++ b/cherab/core/math/integrators/__init__.pxd @@ -0,0 +1,20 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 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.math.integrators.integrators1d cimport Integrator1D, GaussianQuadrature + diff --git a/cherab/core/math/integrators/__init__.py b/cherab/core/math/integrators/__init__.py new file mode 100644 index 00000000..86b7d58d --- /dev/null +++ b/cherab/core/math/integrators/__init__.py @@ -0,0 +1,19 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 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 .integrators1d import Integrator1D, GaussianQuadrature diff --git a/cherab/core/math/integrators/integrators1d.pxd b/cherab/core/math/integrators/integrators1d.pxd new file mode 100644 index 00000000..56ec342b --- /dev/null +++ b/cherab/core/math/integrators/integrators1d.pxd @@ -0,0 +1,36 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 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 numpy cimport ndarray +from raysect.core.math.function.float cimport Function1D + + +cdef class Integrator1D: + + cpdef (double, double) integrate(self, Function1D func, double a, double b) + + +cdef class GaussianQuadrature(Integrator1D): + + cdef: + int _min_order, _max_order + double _rtol + ndarray _roots, _weights + double[:] _roots_mv, _weights_mv + + cdef _build_cash(self) diff --git a/cherab/core/math/integrators/integrators1d.pyx b/cherab/core/math/integrators/integrators1d.pyx new file mode 100644 index 00000000..177be590 --- /dev/null +++ b/cherab/core/math/integrators/integrators1d.pyx @@ -0,0 +1,199 @@ +# cython: language_level=3 + +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 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 numpy as np +from scipy.special import roots_legendre + +from libc.math cimport INFINITY +cimport cython + + +cdef class Integrator1D: + """ + Compute a definite integral of a one-dimensional function. + """ + + cpdef (double, double) integrate(self, Function1D func, double a, double b): + """ + Integrates a one-dimensional function over the given interval. + + :param Function1D func: A function to integrate. + :param double a: Lower limit of integration. + :param double b: Upper limit of integration. + + :returns: Two-element tuple containing the integral of func from a to b + and an estimate of the absolute error in the result. + """ + + raise NotImplementedError("The integrate() virtual method must be implemented.") + + +cdef class GaussianQuadrature(Integrator1D): + """ + Compute an integral of a one-dimensional function over a finite interval + using fixed-tolerance Gaussian quadrature. + (see Scipy `quadrature `). + + :param double relative_tolerance: Iteration stops when relative error between + last two iterates is less than this value. Default is 1.e-7. + :param int max_order: Maximum order on Gaussian quadrature. Default is 50. + :param int min_order: Minimum order on Gaussian quadrature. Default is 1. + + :ivar double relative_tolerance: Iteration stops when relative error between + last two iterates is less than this value. + :ivar int max_order: Maximum order on Gaussian quadrature. + :ivar int min_order: Minimum order on Gaussian quadrature. + """ + + def __init__(self, double relative_tolerance=1.e-7, int max_order=50, int min_order=1): + + if min_order < 1 or max_order < 1: + raise ValueError("Order of Gaussian quadrature must be >= 1.") + + if min_order >= max_order: + raise ValueError("Minimum order of Gaussian quadrature must be less than the maximum order.") + + self._min_order = min_order + self._max_order = max_order + self._build_cash() + + self.relative_tolerance = relative_tolerance + + @property + def min_order(self): + """ + Minimum order on Gaussian quadrature. + + :rtype: int + """ + return self._min_order + + @min_order.setter + def min_order(self, int value): + + if value < 1: + raise ValueError("Order of Gaussian quadrature must be >= 1.") + + if value > self._max_order: + raise ValueError("Minimum order of Gaussian quadrature must be less than or equal to the maximum order.") + + self._min_order = value + + self._build_cash() + + @property + def max_order(self): + """ + Maximum order on Gaussian quadrature. + + :rtype: float + """ + return self._max_order + + @max_order.setter + def max_order(self, int value): + + if value < 1: + raise ValueError("Order of Gaussian quadrature must be >= 1.") + + if value < self._min_order: + raise ValueError("Maximum order of Gaussian quadrature must be greater than or equal to the minimum order.") + + self._max_order = value + + self._build_cash() + + @property + def relative_tolerance(self): + """ + Iteration stops when relative error between last two iterates is less than this value. + + :rtype: double + """ + return self._rtol + + @relative_tolerance.setter + def relative_tolerance(self, double value): + + if value <= 0: + raise ValueError("Relative tolerance must be positive.") + + self._rtol = value + + cdef _build_cash(self): + """ + Caches the roots and weights of the Gauss-Legendre quadrature. + """ + + cdef: + int order, n, i + + n = (self._max_order + self._min_order) * (self._max_order - self._min_order + 1) // 2 + + self._roots = np.zeros(n, dtype=np.float64) + self._weights = np.zeros(n, dtype=np.float64) + + i = 0 + for order in range(self._min_order, self._max_order + 1): + self._roots[i:i + order], self._weights[i:i + order] = roots_legendre(order) + i += order + + self._roots_mv = self._roots + self._weights_mv = self._weights + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + @cython.initializedcheck(False) + cpdef (double, double) integrate(self, Function1D func, double a, double b): + """ + Integrates a one-dimensional function over a finite interval. + + :param Function1D func: A function to integrate. + :param double a: Lower limit of integration. + :param double b: Upper limit of integration. + + :returns: Two-element tuple containing Gaussian quadrature approximation to integral + and difference between last two estimates of the integral. + """ + + cdef: + int order, i, ibegin + double newval, oldval, error, x + + oldval = INFINITY + ibegin = 0 + + for order in range(self._min_order, self._max_order): + newval = 0 + for i in range(ibegin, ibegin + order): + x = a + 0.5 * (b - a) * (self._roots_mv[i] + 1.) + newval += self._weights_mv[i] * func.evaluate(x) + newval *= 0.5 * (b - a) + + error = abs(newval - oldval) + oldval = newval + + ibegin += order + + if error < self._rtol * abs(newval): + break + + return newval, error diff --git a/cherab/core/math/tests/test_integrators.py b/cherab/core/math/tests/test_integrators.py new file mode 100644 index 00000000..4e563cb6 --- /dev/null +++ b/cherab/core/math/tests/test_integrators.py @@ -0,0 +1,78 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 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.core.math.function.float import Exp1D, Arg1D +from cherab.core.math.integrators import GaussianQuadrature +from math import sqrt, pi +from scipy.special import erf +import unittest + + +class TestGaussianQuadrature(unittest.TestCase): + """Gaussian quadrature integrator tests.""" + + def test_properties(self): + """Test property assignment.""" + min_order = 3 + max_order = 30 + reltol = 1.e-6 + quadrature = GaussianQuadrature(relative_tolerance=reltol, max_order=max_order, min_order=min_order) + + self.assertEqual(quadrature.relative_tolerance, reltol) + self.assertEqual(quadrature.max_order, max_order) + self.assertEqual(quadrature.min_order, min_order) + + min_order = 0 + max_order = 2 # < min_order + reltol = -1 + + with self.assertRaises(ValueError): + quadrature.max_order = max_order + + with self.assertRaises(ValueError): + quadrature.min_order = min_order + + with self.assertRaises(ValueError): + quadrature.relative_tolerance = reltol + + min_order = 1 + max_order = 20 + reltol = 1.e-5 + + quadrature.relative_tolerance = reltol + quadrature.min_order = min_order + quadrature.max_order = max_order + + self.assertEqual(quadrature.relative_tolerance, reltol) + self.assertEqual(quadrature.min_order, min_order) + self.assertEqual(quadrature.max_order, max_order) + + def test_integrate(self): + """Test integration.""" + func = (2 / sqrt(pi)) * Exp1D(- Arg1D() * Arg1D()) + quadrature = GaussianQuadrature() + a = -0.5 + b = 3. + result, error = quadrature.integrate(func, a, b) + exact_integral = erf(b) - erf(a) + + self.assertAlmostEqual(result, exact_integral, places=7) + + +if __name__ == '__main__': + unittest.main() From a32a6a716ae4a8ce8b46103163c0c8b793cf9fa2 Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Wed, 31 Aug 2022 19:11:34 +0300 Subject: [PATCH 2/9] Fix an error in GaussianQuadrature. Increase default relative tolerance value. --- cherab/core/math/integrators/integrators1d.pyx | 10 +++++----- cherab/core/math/tests/test_integrators.py | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/cherab/core/math/integrators/integrators1d.pyx b/cherab/core/math/integrators/integrators1d.pyx index 177be590..7f9ce421 100644 --- a/cherab/core/math/integrators/integrators1d.pyx +++ b/cherab/core/math/integrators/integrators1d.pyx @@ -52,7 +52,7 @@ cdef class GaussianQuadrature(Integrator1D): (see Scipy `quadrature `). :param double relative_tolerance: Iteration stops when relative error between - last two iterates is less than this value. Default is 1.e-7. + last two iterates is less than this value. Default is 1.e-5. :param int max_order: Maximum order on Gaussian quadrature. Default is 50. :param int min_order: Minimum order on Gaussian quadrature. Default is 1. @@ -62,13 +62,13 @@ cdef class GaussianQuadrature(Integrator1D): :ivar int min_order: Minimum order on Gaussian quadrature. """ - def __init__(self, double relative_tolerance=1.e-7, int max_order=50, int min_order=1): + def __init__(self, double relative_tolerance=1.e-5, int max_order=50, int min_order=1): if min_order < 1 or max_order < 1: raise ValueError("Order of Gaussian quadrature must be >= 1.") - if min_order >= max_order: - raise ValueError("Minimum order of Gaussian quadrature must be less than the maximum order.") + if min_order > max_order: + raise ValueError("Minimum order of Gaussian quadrature must be less than or equal to the maximum order.") self._min_order = min_order self._max_order = max_order @@ -181,7 +181,7 @@ cdef class GaussianQuadrature(Integrator1D): oldval = INFINITY ibegin = 0 - for order in range(self._min_order, self._max_order): + for order in range(self._min_order, self._max_order + 1): newval = 0 for i in range(ibegin, ibegin + order): x = a + 0.5 * (b - a) * (self._roots_mv[i] + 1.) diff --git a/cherab/core/math/tests/test_integrators.py b/cherab/core/math/tests/test_integrators.py index 4e563cb6..a60a4eea 100644 --- a/cherab/core/math/tests/test_integrators.py +++ b/cherab/core/math/tests/test_integrators.py @@ -65,13 +65,13 @@ def test_properties(self): def test_integrate(self): """Test integration.""" func = (2 / sqrt(pi)) * Exp1D(- Arg1D() * Arg1D()) - quadrature = GaussianQuadrature() + quadrature = GaussianQuadrature(relative_tolerance=1.e-8) a = -0.5 b = 3. result, error = quadrature.integrate(func, a, b) exact_integral = erf(b) - erf(a) - self.assertAlmostEqual(result, exact_integral, places=7) + self.assertAlmostEqual(result, exact_integral, places=8) if __name__ == '__main__': From c4a604d5133ea302a33bd61253aa3af781029c7a Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Wed, 31 Aug 2022 19:13:34 +0300 Subject: [PATCH 3/9] Add NumericallyIntegrableLineShapeModel class. Update StarkBroadenedLine to use numerical integration over the spectral bin. --- cherab/core/model/lineshape.pxd | 8 +++- cherab/core/model/lineshape.pyx | 65 ++++++++++++++++++++++++++------- 2 files changed, 58 insertions(+), 15 deletions(-) diff --git a/cherab/core/model/lineshape.pxd b/cherab/core/model/lineshape.pxd index 711475b5..050d589e 100644 --- a/cherab/core/model/lineshape.pxd +++ b/cherab/core/model/lineshape.pxd @@ -24,6 +24,7 @@ cimport numpy as np from raysect.optical cimport Spectrum, Point3D, Vector3D from cherab.core cimport Line, Species, Plasma, Beam from cherab.core.math cimport Function1D, Function2D +from cherab.core.math.integrators cimport Integrator1D from cherab.core.atomic.zeeman cimport ZeemanStructure @@ -45,6 +46,11 @@ cdef class LineShapeModel: cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum) +cdef class NumericallyIntegrableLineShapeModel(LineShapeModel): + + cdef Integrator1D integrator + + cdef class GaussianLine(LineShapeModel): pass @@ -57,7 +63,7 @@ cdef class MultipletLineShape(LineShapeModel): double[:,::1] _multiplet_mv -cdef class StarkBroadenedLine(LineShapeModel): +cdef class StarkBroadenedLine(NumericallyIntegrableLineShapeModel): cdef double _aij, _bij, _cij diff --git a/cherab/core/model/lineshape.pyx b/cherab/core/model/lineshape.pyx index 5e346175..c0bb9711 100644 --- a/cherab/core/model/lineshape.pyx +++ b/cherab/core/model/lineshape.pyx @@ -22,10 +22,12 @@ import numpy as np cimport numpy as np from libc.math cimport sqrt, erf, M_SQRT2, floor, ceil, fabs from raysect.optical.spectrum cimport new_spectrum +from raysect.core.math.function.float cimport Function1D from cherab.core cimport Plasma from cherab.core.atomic.elements import hydrogen, deuterium, tritium, helium, helium3, beryllium, boron, carbon, nitrogen, oxygen, neon from cherab.core.math.function cimport autowrap_function1d, autowrap_function2d +from cherab.core.math.integrators cimport GaussianQuadrature from cherab.core.utility.constants cimport ATOMIC_MASS, ELEMENTARY_CHARGE, SPEED_OF_LIGHT cimport cython @@ -159,6 +161,24 @@ cdef class LineShapeModel: raise NotImplementedError('Child lineshape class must implement this method.') +cdef class NumericallyIntegrableLineShapeModel(LineShapeModel): + """ + Line shape model to be numerically integrated over the spectral bin. + + :param Line line: The emission line object for this line shape. + :param float wavelength: The rest wavelength for this emission line. + :param Species target_species: The target plasma species that is emitting. + :param Plasma plasma: The emitting plasma object. + :param Integrator1D integrator: Integrator1D instance to integrate the line shape + over the spectral bin. Default is `GaussianQuadrature()`. + """ + + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, Integrator1D integrator=None): + + super().__init__(line, wavelength, target_species, plasma) + self.integrator = integrator or GaussianQuadrature() + + cdef class GaussianLine(LineShapeModel): """ Produces Gaussian line shape. @@ -289,7 +309,21 @@ cdef class MultipletLineShape(LineShapeModel): return spectrum -cdef class StarkBroadenedLine(LineShapeModel): +cdef class StarkFunction(Function1D): + + cdef double _a, _x0 + + def __init__(self, double wavelength, double lambda_1_2): + self._x0 = wavelength + self._a = (0.5 * lambda_1_2)**2.5 + + @cython.cdivision(True) + cdef double evaluate(self, double x) except? -1e999: + + return 1. / ((fabs(x - self._x0))**2.5 + self._a) + + +cdef class StarkBroadenedLine(NumericallyIntegrableLineShapeModel): """ Parametrised Stark broadened line shape based on the Model Microfield Method (MMM). Contains embedded atomic data in the form of fits to MMM. @@ -352,7 +386,8 @@ cdef class StarkBroadenedLine(LineShapeModel): Line(tritium, 0, (9, 3)): (5.588e-15, 0.7165, 0.033) } - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, dict stark_model_coefficients=None): + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, Integrator1D integrator=None, + dict stark_model_coefficients=None): stark_model_coefficients = stark_model_coefficients or self.STARK_MODEL_COEFFICIENTS_DEFAULT @@ -371,7 +406,7 @@ cdef class StarkBroadenedLine(LineShapeModel): except IndexError: raise ValueError('Stark broadening coefficients for {} is not currently available.'.format(line)) - super().__init__(line, wavelength, target_species, plasma) + super().__init__(line, wavelength, target_species, plasma, integrator) def show_supported_transitions(self): """ Prints all supported transitions.""" @@ -384,11 +419,14 @@ cdef class StarkBroadenedLine(LineShapeModel): @cython.cdivision(True) cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): - cdef double ne, te, lambda_1_2, lambda_5_2, wvl - cdef double cutoff_lower_wavelength, cutoff_upper_wavelength - cdef double lower_value, lower_wavelength, upper_value, upper_wavelength - cdef int start, end, i - cdef Spectrum raw_lineshape + cdef: + double ne, te, lambda_1_2, lambda_5_2, wvl + double cutoff_lower_wavelength, cutoff_upper_wavelength + double lower_wavelength, upper_wavelength + double bin_integral + int start, end, i + Spectrum raw_lineshape + StarkFunction stark_funcion ne = self.plasma.get_electron_distribution().density(point.x, point.y, point.z) if ne <= 0.0: @@ -400,6 +438,8 @@ cdef class StarkBroadenedLine(LineShapeModel): lambda_1_2 = self._cij * ne**self._aij / (te**self._bij) + stark_funcion = StarkFunction(self.wavelength, lambda_1_2) + # calculate and check end of limits cutoff_lower_wavelength = self.wavelength - LORENZIAN_CUTOFF_GAMMA * lambda_1_2 if spectrum.max_wavelength < cutoff_lower_wavelength: @@ -413,21 +453,18 @@ cdef class StarkBroadenedLine(LineShapeModel): start = max(0, floor((cutoff_lower_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) end = min(spectrum.bins, ceil((cutoff_upper_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) - # TODO - replace with cumulative integrals # add line to spectrum raw_lineshape = spectrum.new_spectrum() lower_wavelength = raw_lineshape.min_wavelength + start * raw_lineshape.delta_wavelength - lower_value = 1 / ((fabs(lower_wavelength - self.wavelength))**2.5 + (0.5 * lambda_1_2)**2.5) - for i in range(start, end): + for i in range(start, end): upper_wavelength = raw_lineshape.min_wavelength + raw_lineshape.delta_wavelength * (i + 1) - upper_value = 1 / ((fabs(upper_wavelength - self.wavelength))**2.5 + (0.5 * lambda_1_2)**2.5) - raw_lineshape.samples_mv[i] += 0.5 * (upper_value + lower_value) + bin_integral, _ = self.integrator.integrate(stark_funcion, lower_wavelength, upper_wavelength) + raw_lineshape.samples_mv[i] += bin_integral lower_wavelength = upper_wavelength - lower_value = upper_value # perform normalisation raw_lineshape.div_scalar(raw_lineshape.total()) From 6308e0fac89ec830697e25c000c6867ceba6edc1 Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Wed, 31 Aug 2022 19:55:23 +0300 Subject: [PATCH 4/9] Update changelog. --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 27b5ecac..1dfaac01 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,9 @@ New: * Add verbose parameter to SartOpencl solver (default is False). (#358) * Add Generomak core plasma profiles. (#360) * Add toroidal_mesh_from_polygon for making mesh for not fully-360 degrees axisymmetric elements. (#365) +* Add GaussianQuadrature integration method for Function1D. (#366) +* Add NumericallyIntegrableLineShapeModel for lineshapes that cannot be analytically integrated over a spectral bin. (#366) +* Add a numerical integration of StarkBroadenedLine over the spectral bin. (#366) Bug Fixes: ---------- From 526cc17342a6c7a3b910f138603c492e30adcc10 Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Sun, 11 Sep 2022 03:05:02 +0300 Subject: [PATCH 5/9] Make Integrator1D a Function2D instance. --- .../core/math/integrators/integrators1d.pxd | 7 +-- .../core/math/integrators/integrators1d.pyx | 54 ++++++++++++------- cherab/core/math/tests/test_integrators.py | 14 +++-- 3 files changed, 48 insertions(+), 27 deletions(-) diff --git a/cherab/core/math/integrators/integrators1d.pxd b/cherab/core/math/integrators/integrators1d.pxd index 56ec342b..8109e6d7 100644 --- a/cherab/core/math/integrators/integrators1d.pxd +++ b/cherab/core/math/integrators/integrators1d.pxd @@ -17,12 +17,13 @@ # under the Licence. from numpy cimport ndarray -from raysect.core.math.function.float cimport Function1D +from raysect.core.math.function.float cimport Function1D, Function2D -cdef class Integrator1D: +cdef class Integrator1D(Function2D): - cpdef (double, double) integrate(self, Function1D func, double a, double b) + cdef: + Function1D function cdef class GaussianQuadrature(Integrator1D): diff --git a/cherab/core/math/integrators/integrators1d.pyx b/cherab/core/math/integrators/integrators1d.pyx index 7f9ce421..1e69efd4 100644 --- a/cherab/core/math/integrators/integrators1d.pyx +++ b/cherab/core/math/integrators/integrators1d.pyx @@ -21,28 +21,35 @@ import numpy as np from scipy.special import roots_legendre +from raysect.core.math.function.float cimport autowrap_function1d + from libc.math cimport INFINITY cimport cython -cdef class Integrator1D: +cdef class Integrator1D(Function2D): """ Compute a definite integral of a one-dimensional function. + + :ivar Function1D integrand: A 1D function to integrate. """ - cpdef (double, double) integrate(self, Function1D func, double a, double b): + @property + def integrand(self): """ - Integrates a one-dimensional function over the given interval. - - :param Function1D func: A function to integrate. - :param double a: Lower limit of integration. - :param double b: Upper limit of integration. + A 1D function to integrate. - :returns: Two-element tuple containing the integral of func from a to b - and an estimate of the absolute error in the result. + :rtype: int """ + return self.function + + @integrand.setter + def integrand(self, object func): - raise NotImplementedError("The integrate() virtual method must be implemented.") + if func is None: + self.function = None + else: + self.function = autowrap_function1d(func) cdef class GaussianQuadrature(Integrator1D): @@ -51,18 +58,20 @@ cdef class GaussianQuadrature(Integrator1D): using fixed-tolerance Gaussian quadrature. (see Scipy `quadrature `). + :param object integrand: A 1D function to integrate. :param double relative_tolerance: Iteration stops when relative error between last two iterates is less than this value. Default is 1.e-5. :param int max_order: Maximum order on Gaussian quadrature. Default is 50. :param int min_order: Minimum order on Gaussian quadrature. Default is 1. + :ivar Function1D integrand: A 1D function to integrate. :ivar double relative_tolerance: Iteration stops when relative error between last two iterates is less than this value. :ivar int max_order: Maximum order on Gaussian quadrature. :ivar int min_order: Minimum order on Gaussian quadrature. """ - def __init__(self, double relative_tolerance=1.e-5, int max_order=50, int min_order=1): + def __init__(self, object integrand=None, double relative_tolerance=1.e-5, int max_order=50, int min_order=1): if min_order < 1 or max_order < 1: raise ValueError("Order of Gaussian quadrature must be >= 1.") @@ -74,6 +83,8 @@ cdef class GaussianQuadrature(Integrator1D): self._max_order = max_order self._build_cash() + self.integrand = integrand + self.relative_tolerance = relative_tolerance @property @@ -162,31 +173,34 @@ cdef class GaussianQuadrature(Integrator1D): @cython.wraparound(False) @cython.cdivision(True) @cython.initializedcheck(False) - cpdef (double, double) integrate(self, Function1D func, double a, double b): + cdef double evaluate(self, double a, double b) except? -1e999: """ Integrates a one-dimensional function over a finite interval. - :param Function1D func: A function to integrate. :param double a: Lower limit of integration. :param double b: Upper limit of integration. - :returns: Two-element tuple containing Gaussian quadrature approximation to integral - and difference between last two estimates of the integral. + :returns: Gaussian quadrature approximation to integral. """ cdef: int order, i, ibegin - double newval, oldval, error, x + double newval, oldval, error, x, c, d + + if self.function is None: + raise AttributeError("Integrand is not set.") oldval = INFINITY ibegin = 0 + c = 0.5 * (a + b) + d = 0.5 * (b - a) for order in range(self._min_order, self._max_order + 1): newval = 0 for i in range(ibegin, ibegin + order): - x = a + 0.5 * (b - a) * (self._roots_mv[i] + 1.) - newval += self._weights_mv[i] * func.evaluate(x) - newval *= 0.5 * (b - a) + x = c + d * self._roots_mv[i] + newval += self._weights_mv[i] * self.function.evaluate(x) + newval *= d error = abs(newval - oldval) oldval = newval @@ -196,4 +210,4 @@ cdef class GaussianQuadrature(Integrator1D): if error < self._rtol * abs(newval): break - return newval, error + return newval diff --git a/cherab/core/math/tests/test_integrators.py b/cherab/core/math/tests/test_integrators.py index a60a4eea..a62ffafd 100644 --- a/cherab/core/math/tests/test_integrators.py +++ b/cherab/core/math/tests/test_integrators.py @@ -31,11 +31,12 @@ def test_properties(self): min_order = 3 max_order = 30 reltol = 1.e-6 - quadrature = GaussianQuadrature(relative_tolerance=reltol, max_order=max_order, min_order=min_order) + quadrature = GaussianQuadrature(integrand=Arg1D, relative_tolerance=reltol, max_order=max_order, min_order=min_order) self.assertEqual(quadrature.relative_tolerance, reltol) self.assertEqual(quadrature.max_order, max_order) self.assertEqual(quadrature.min_order, min_order) + self.assertEqual(quadrature.integrand, Arg1D) min_order = 0 max_order = 2 # < min_order @@ -57,21 +58,26 @@ def test_properties(self): quadrature.relative_tolerance = reltol quadrature.min_order = min_order quadrature.max_order = max_order + quadrature.integrand = Exp1D self.assertEqual(quadrature.relative_tolerance, reltol) self.assertEqual(quadrature.min_order, min_order) self.assertEqual(quadrature.max_order, max_order) + self.assertEqual(quadrature.integrand, Exp1D) def test_integrate(self): """Test integration.""" - func = (2 / sqrt(pi)) * Exp1D(- Arg1D() * Arg1D()) quadrature = GaussianQuadrature(relative_tolerance=1.e-8) a = -0.5 b = 3. - result, error = quadrature.integrate(func, a, b) + + with self.assertRaises(AttributeError): # integrand is not set + quadrature(a, b) + + quadrature.integrand = (2 / sqrt(pi)) * Exp1D(- Arg1D() * Arg1D()) exact_integral = erf(b) - erf(a) - self.assertAlmostEqual(result, exact_integral, places=8) + self.assertAlmostEqual(quadrature(a, b), exact_integral, places=8) if __name__ == '__main__': From 86d908b4fe4cb7f523921992925c1e019914cedc Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Sun, 11 Sep 2022 03:08:35 +0300 Subject: [PATCH 6/9] Remove NumericallyIntegratedLineShapeModel. Add integrator to LineShapeModel. --- cherab/core/model/lineshape.pxd | 8 ++----- cherab/core/model/lineshape.pyx | 37 +++++++++++---------------------- 2 files changed, 14 insertions(+), 31 deletions(-) diff --git a/cherab/core/model/lineshape.pxd b/cherab/core/model/lineshape.pxd index 050d589e..2591efe3 100644 --- a/cherab/core/model/lineshape.pxd +++ b/cherab/core/model/lineshape.pxd @@ -42,15 +42,11 @@ cdef class LineShapeModel: double wavelength Species target_species Plasma plasma + Integrator1D integrator cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum) -cdef class NumericallyIntegrableLineShapeModel(LineShapeModel): - - cdef Integrator1D integrator - - cdef class GaussianLine(LineShapeModel): pass @@ -63,7 +59,7 @@ cdef class MultipletLineShape(LineShapeModel): double[:,::1] _multiplet_mv -cdef class StarkBroadenedLine(NumericallyIntegrableLineShapeModel): +cdef class StarkBroadenedLine(LineShapeModel): cdef double _aij, _bij, _cij diff --git a/cherab/core/model/lineshape.pyx b/cherab/core/model/lineshape.pyx index c0bb9711..96c36495 100644 --- a/cherab/core/model/lineshape.pyx +++ b/cherab/core/model/lineshape.pyx @@ -148,37 +148,22 @@ cdef class LineShapeModel: :param float wavelength: The rest wavelength for this emission line. :param Species target_species: The target plasma species that is emitting. :param Plasma plasma: The emitting plasma object. + :param Integrator1D integrator: Integrator1D instance to integrate the line shape + over the spectral bin. Default is None. """ - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma): + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, Integrator1D integrator=None): self.line = line self.wavelength = wavelength self.target_species = target_species self.plasma = plasma + self.integrator = integrator cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): raise NotImplementedError('Child lineshape class must implement this method.') -cdef class NumericallyIntegrableLineShapeModel(LineShapeModel): - """ - Line shape model to be numerically integrated over the spectral bin. - - :param Line line: The emission line object for this line shape. - :param float wavelength: The rest wavelength for this emission line. - :param Species target_species: The target plasma species that is emitting. - :param Plasma plasma: The emitting plasma object. - :param Integrator1D integrator: Integrator1D instance to integrate the line shape - over the spectral bin. Default is `GaussianQuadrature()`. - """ - - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, Integrator1D integrator=None): - - super().__init__(line, wavelength, target_species, plasma) - self.integrator = integrator or GaussianQuadrature() - - cdef class GaussianLine(LineShapeModel): """ Produces Gaussian line shape. @@ -323,7 +308,7 @@ cdef class StarkFunction(Function1D): return 1. / ((fabs(x - self._x0))**2.5 + self._a) -cdef class StarkBroadenedLine(NumericallyIntegrableLineShapeModel): +cdef class StarkBroadenedLine(LineShapeModel): """ Parametrised Stark broadened line shape based on the Model Microfield Method (MMM). Contains embedded atomic data in the form of fits to MMM. @@ -342,6 +327,9 @@ cdef class StarkBroadenedLine(NumericallyIntegrableLineShapeModel): :param dict stark_model_coefficients: Alternative model coefficients in the form {line_ij: (c_ij, a_ij, b_ij), ...}. If None, the default model parameters will be used. + :param Integrator1D integrator: Integrator1D instance to integrate the line shape + over the spectral bin. Default is `GaussianQuadrature()`. + """ STARK_MODEL_COEFFICIENTS_DEFAULT = { @@ -386,8 +374,8 @@ cdef class StarkBroadenedLine(NumericallyIntegrableLineShapeModel): Line(tritium, 0, (9, 3)): (5.588e-15, 0.7165, 0.033) } - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, Integrator1D integrator=None, - dict stark_model_coefficients=None): + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, + dict stark_model_coefficients=None, integrator=GaussianQuadrature()): stark_model_coefficients = stark_model_coefficients or self.STARK_MODEL_COEFFICIENTS_DEFAULT @@ -426,7 +414,6 @@ cdef class StarkBroadenedLine(NumericallyIntegrableLineShapeModel): double bin_integral int start, end, i Spectrum raw_lineshape - StarkFunction stark_funcion ne = self.plasma.get_electron_distribution().density(point.x, point.y, point.z) if ne <= 0.0: @@ -438,7 +425,7 @@ cdef class StarkBroadenedLine(NumericallyIntegrableLineShapeModel): lambda_1_2 = self._cij * ne**self._aij / (te**self._bij) - stark_funcion = StarkFunction(self.wavelength, lambda_1_2) + self.integrator.function = StarkFunction(self.wavelength, lambda_1_2) # calculate and check end of limits cutoff_lower_wavelength = self.wavelength - LORENZIAN_CUTOFF_GAMMA * lambda_1_2 @@ -461,7 +448,7 @@ cdef class StarkBroadenedLine(NumericallyIntegrableLineShapeModel): for i in range(start, end): upper_wavelength = raw_lineshape.min_wavelength + raw_lineshape.delta_wavelength * (i + 1) - bin_integral, _ = self.integrator.integrate(stark_funcion, lower_wavelength, upper_wavelength) + bin_integral = self.integrator.evaluate(lower_wavelength, upper_wavelength) raw_lineshape.samples_mv[i] += bin_integral lower_wavelength = upper_wavelength From 6a915f9edc1c584b3c89856efabbea0da400cec3 Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Mon, 17 Oct 2022 18:59:44 +0300 Subject: [PATCH 7/9] Remove Function2D inheritance for Integrator1D. Forbid setting the integrand to None. --- .../core/math/integrators/integrators1d.pxd | 6 ++- .../core/math/integrators/integrators1d.pyx | 42 ++++++++++++------- cherab/core/math/tests/test_integrators.py | 4 -- 3 files changed, 30 insertions(+), 22 deletions(-) diff --git a/cherab/core/math/integrators/integrators1d.pxd b/cherab/core/math/integrators/integrators1d.pxd index 8109e6d7..c451285f 100644 --- a/cherab/core/math/integrators/integrators1d.pxd +++ b/cherab/core/math/integrators/integrators1d.pxd @@ -20,11 +20,13 @@ from numpy cimport ndarray from raysect.core.math.function.float cimport Function1D, Function2D -cdef class Integrator1D(Function2D): +cdef class Integrator1D: cdef: Function1D function + cdef double evaluate(self, double a, double b) except? -1e999 + cdef class GaussianQuadrature(Integrator1D): @@ -34,4 +36,4 @@ cdef class GaussianQuadrature(Integrator1D): ndarray _roots, _weights double[:] _roots_mv, _weights_mv - cdef _build_cash(self) + cdef _build_cache(self) diff --git a/cherab/core/math/integrators/integrators1d.pyx b/cherab/core/math/integrators/integrators1d.pyx index 1e69efd4..54ec4059 100644 --- a/cherab/core/math/integrators/integrators1d.pyx +++ b/cherab/core/math/integrators/integrators1d.pyx @@ -21,13 +21,13 @@ import numpy as np from scipy.special import roots_legendre -from raysect.core.math.function.float cimport autowrap_function1d +from raysect.core.math.function.float cimport autowrap_function1d, Constant1D from libc.math cimport INFINITY cimport cython -cdef class Integrator1D(Function2D): +cdef class Integrator1D: """ Compute a definite integral of a one-dimensional function. @@ -44,12 +44,25 @@ cdef class Integrator1D(Function2D): return self.function @integrand.setter - def integrand(self, object func): + def integrand(self, object func not None): - if func is None: - self.function = None - else: - self.function = autowrap_function1d(func) + self.function = autowrap_function1d(func) + + cdef double evaluate(self, double a, double b) except? -1e999: + + raise NotImplementedError("The evaluate() method has not been implemented.") + + def __call__(self, double a, double b): + """ + Integrates a one-dimensional function over a finite interval. + + :param double a: Lower limit of integration. + :param double b: Upper limit of integration. + + :returns: Definite integral of a one-dimensional function. + """ + + return self.evaluate(a, b) cdef class GaussianQuadrature(Integrator1D): @@ -58,7 +71,7 @@ cdef class GaussianQuadrature(Integrator1D): using fixed-tolerance Gaussian quadrature. (see Scipy `quadrature `). - :param object integrand: A 1D function to integrate. + :param object integrand: A 1D function to integrate. Default is Constant1D(0). :param double relative_tolerance: Iteration stops when relative error between last two iterates is less than this value. Default is 1.e-5. :param int max_order: Maximum order on Gaussian quadrature. Default is 50. @@ -71,7 +84,7 @@ cdef class GaussianQuadrature(Integrator1D): :ivar int min_order: Minimum order on Gaussian quadrature. """ - def __init__(self, object integrand=None, double relative_tolerance=1.e-5, int max_order=50, int min_order=1): + def __init__(self, object integrand=Constant1D(0), double relative_tolerance=1.e-5, int max_order=50, int min_order=1): if min_order < 1 or max_order < 1: raise ValueError("Order of Gaussian quadrature must be >= 1.") @@ -81,7 +94,7 @@ cdef class GaussianQuadrature(Integrator1D): self._min_order = min_order self._max_order = max_order - self._build_cash() + self._build_cache() self.integrand = integrand @@ -107,7 +120,7 @@ cdef class GaussianQuadrature(Integrator1D): self._min_order = value - self._build_cash() + self._build_cache() @property def max_order(self): @@ -129,7 +142,7 @@ cdef class GaussianQuadrature(Integrator1D): self._max_order = value - self._build_cash() + self._build_cache() @property def relative_tolerance(self): @@ -148,7 +161,7 @@ cdef class GaussianQuadrature(Integrator1D): self._rtol = value - cdef _build_cash(self): + cdef _build_cache(self): """ Caches the roots and weights of the Gauss-Legendre quadrature. """ @@ -187,9 +200,6 @@ cdef class GaussianQuadrature(Integrator1D): int order, i, ibegin double newval, oldval, error, x, c, d - if self.function is None: - raise AttributeError("Integrand is not set.") - oldval = INFINITY ibegin = 0 c = 0.5 * (a + b) diff --git a/cherab/core/math/tests/test_integrators.py b/cherab/core/math/tests/test_integrators.py index a62ffafd..fa967165 100644 --- a/cherab/core/math/tests/test_integrators.py +++ b/cherab/core/math/tests/test_integrators.py @@ -70,10 +70,6 @@ def test_integrate(self): quadrature = GaussianQuadrature(relative_tolerance=1.e-8) a = -0.5 b = 3. - - with self.assertRaises(AttributeError): # integrand is not set - quadrature(a, b) - quadrature.integrand = (2 / sqrt(pi)) * Exp1D(- Arg1D() * Arg1D()) exact_integral = erf(b) - erf(a) From e62f0adde3eae62cdf949a2db8013f352c5167ea Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Mon, 17 Oct 2022 21:52:05 +0300 Subject: [PATCH 8/9] Add initialisation checks to StarkFunction and make it return a normalised line shape. --- cherab/core/model/lineshape.pyx | 41 +++++++++++++++++----------- cherab/core/tests/test_lineshapes.py | 32 +++++++++++++--------- 2 files changed, 44 insertions(+), 29 deletions(-) diff --git a/cherab/core/model/lineshape.pyx b/cherab/core/model/lineshape.pyx index 96c36495..69278560 100644 --- a/cherab/core/model/lineshape.pyx +++ b/cherab/core/model/lineshape.pyx @@ -19,6 +19,8 @@ # under the Licence. import numpy as np +from scipy.special import hyp2f1 + cimport numpy as np from libc.math cimport sqrt, erf, M_SQRT2, floor, ceil, fabs from raysect.optical.spectrum cimport new_spectrum @@ -80,8 +82,6 @@ cpdef double thermal_broadening(double wavelength, double temperature, double at # the number of standard deviations outside the rest wavelength the line is considered to add negligible value (including a margin for safety) DEF GAUSSIAN_CUTOFF_SIGMA = 10.0 -DEF LORENZIAN_CUTOFF_GAMMA = 50.0 - @cython.cdivision(True) @cython.initializedcheck(False) @@ -294,18 +294,36 @@ cdef class MultipletLineShape(LineShapeModel): return spectrum +DEF LORENZIAN_CUTOFF_GAMMA = 50.0 + + cdef class StarkFunction(Function1D): + """ + Normalised Stark function for the StarkBroadenedLine line shape. + """ - cdef double _a, _x0 + cdef double _a, _x0, _norm + + STARK_NORM_COEFFICIENT = 4 * LORENZIAN_CUTOFF_GAMMA * hyp2f1(0.4, 1, 1.4, -(2 * LORENZIAN_CUTOFF_GAMMA)**2.5) def __init__(self, double wavelength, double lambda_1_2): + + if wavelength <= 0: + raise ValueError("Argument 'wavelength' must be positive.") + + if lambda_1_2 <= 0: + raise ValueError("Argument 'lambda_1_2' must be positive.") + self._x0 = wavelength self._a = (0.5 * lambda_1_2)**2.5 + # normalise, so the integral over x is equal to 1 in the limits + # (_x0 - LORENZIAN_CUTOFF_GAMMA * lambda_1_2, _x0 + LORENZIAN_CUTOFF_GAMMA * lambda_1_2) + self._norm = (0.5 * lambda_1_2)**1.5 / self.STARK_NORM_COEFFICIENT @cython.cdivision(True) cdef double evaluate(self, double x) except? -1e999: - return 1. / ((fabs(x - self._x0))**2.5 + self._a) + return self._norm / ((fabs(x - self._x0))**2.5 + self._a) cdef class StarkBroadenedLine(LineShapeModel): @@ -441,25 +459,16 @@ cdef class StarkBroadenedLine(LineShapeModel): end = min(spectrum.bins, ceil((cutoff_upper_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) # add line to spectrum - raw_lineshape = spectrum.new_spectrum() - - lower_wavelength = raw_lineshape.min_wavelength + start * raw_lineshape.delta_wavelength + lower_wavelength = spectrum.min_wavelength + start * spectrum.delta_wavelength for i in range(start, end): - upper_wavelength = raw_lineshape.min_wavelength + raw_lineshape.delta_wavelength * (i + 1) + upper_wavelength = spectrum.min_wavelength + spectrum.delta_wavelength * (i + 1) bin_integral = self.integrator.evaluate(lower_wavelength, upper_wavelength) - raw_lineshape.samples_mv[i] += bin_integral + spectrum.samples_mv[i] += radiance * bin_integral / spectrum.delta_wavelength lower_wavelength = upper_wavelength - # perform normalisation - raw_lineshape.div_scalar(raw_lineshape.total()) - - for i in range(start, end): - # Radiance ??? - spectrum.samples_mv[i] += radiance * raw_lineshape.samples_mv[i] - return spectrum diff --git a/cherab/core/tests/test_lineshapes.py b/cherab/core/tests/test_lineshapes.py index c29f1f2e..da4dc41d 100644 --- a/cherab/core/tests/test_lineshapes.py +++ b/cherab/core/tests/test_lineshapes.py @@ -18,17 +18,17 @@ import unittest -import os import numpy as np -from scipy.special import erf +from scipy.special import erf, hyp2f1 +from scipy.integrate import quadrature from raysect.core import Point3D, Vector3D from raysect.core.math.function.float import Arg1D, Constant1D from raysect.optical import Spectrum from cherab.core import Line +from cherab.core.math.integrators import GaussianQuadrature from cherab.core.atomic import deuterium, nitrogen, ZeemanStructure -from cherab.openadas import OpenADAS from cherab.tools.plasmas.slab import build_constant_slab_plasma from cherab.core.model import GaussianLine, MultipletLineShape, StarkBroadenedLine, ZeemanTriplet, ParametrisedZeemanTriplet, ZeemanMultiplet @@ -116,7 +116,7 @@ def test_multiplet_line_shape(self): for i in range(bins): self.assertAlmostEqual(multi_gaussian[i], spectrum.samples[i], delta=1e-10, - msg='GaussianLine.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) + msg='MultipletLineShape.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) def test_zeeman_triplet(self): # setting up a line shape model @@ -167,7 +167,7 @@ def test_zeeman_triplet(self): for pol in ('no', 'pi', 'sigma'): for i in range(bins): self.assertAlmostEqual(tri_gaussian[pol][i], spectrum[pol].samples[i], delta=1e-10, - msg='GaussianLine.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) + msg='ZeemanTriplet.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) def test_parametrised_zeeman_triplet(self): # setting up a line shape model @@ -219,7 +219,7 @@ def test_parametrised_zeeman_triplet(self): for pol in ('no', 'pi', 'sigma'): for i in range(bins): self.assertAlmostEqual(tri_gaussian[pol][i], spectrum[pol].samples[i], delta=1e-10, - msg='GaussianLine.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) + msg='ParametrisedZeemanTriplet.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) def test_zeeman_multiplet(self): # setting up a line shape model @@ -277,14 +277,15 @@ def test_zeeman_multiplet(self): for pol in ('no', 'pi', 'sigma'): for i in range(bins): self.assertAlmostEqual(tri_gaussian[pol][i], spectrum[pol].samples[i], delta=1e-10, - msg='GaussianLine.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) + msg='ZeemanMultiplet.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) def test_stark_broadened_line(self): # setting up a line shape model line = Line(deuterium, 0, (6, 2)) # D-delta line target_species = self.plasma.composition.get(line.element, line.charge) wavelength = 656.104 - stark_line = StarkBroadenedLine(line, wavelength, target_species, self.plasma) + integrator = GaussianQuadrature(relative_tolerance=1.e-5) + stark_line = StarkBroadenedLine(line, wavelength, target_species, self.plasma, integrator=integrator) # spectrum parameters min_wavelength = wavelength - 0.2 @@ -304,14 +305,19 @@ def test_stark_broadened_line(self): te = self.plasma.electron_distribution.effective_temperature(point.x, point.y, point.z) lambda_1_2 = cij * ne**aij / (te**bij) + lorenzian_cutoff_gamma = 50 + stark_norm_coeff = 4 * lorenzian_cutoff_gamma * hyp2f1(0.4, 1, 1.4, -(2 * lorenzian_cutoff_gamma)**2.5) + norm = (0.5 * lambda_1_2)**1.5 / stark_norm_coeff + wavelengths, delta = np.linspace(min_wavelength, max_wavelength, bins + 1, retstep=True) - stark_lineshape = 1 / ((np.abs(wavelengths - wavelength))**2.5 + (0.5 * lambda_1_2)**2.5) - stark_lineshape = 0.5 * (stark_lineshape[1:] + stark_lineshape[:-1]) - stark_lineshape /= stark_lineshape.sum() * delta + + def stark_lineshape(x): + return norm / ((np.abs(x - wavelength))**2.5 + (0.5 * lambda_1_2)**2.5) for i in range(bins): - self.assertAlmostEqual(stark_lineshape[i], spectrum.samples[i], delta=1e-10, - msg='GaussianLine.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) + stark_bin = quadrature(stark_lineshape, wavelengths[i], wavelengths[i + 1], rtol=integrator.relative_tolerance)[0] / delta + self.assertAlmostEqual(stark_bin, spectrum.samples[i], delta=1e-9, + msg='StarkBroadenedLine.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) if __name__ == '__main__': From 308208ceccfa8af6cfbab4c8aa308234696a6137 Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Tue, 18 Oct 2022 23:14:03 +0300 Subject: [PATCH 9/9] Update CHANGELOG.md --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1dfaac01..edc760b8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,7 +15,7 @@ New: * Add Generomak core plasma profiles. (#360) * Add toroidal_mesh_from_polygon for making mesh for not fully-360 degrees axisymmetric elements. (#365) * Add GaussianQuadrature integration method for Function1D. (#366) -* Add NumericallyIntegrableLineShapeModel for lineshapes that cannot be analytically integrated over a spectral bin. (#366) +* Add integrator attribute to LineShapeModel to use with lineshapes that cannot be analytically integrated over a spectral bin. (#366) * Add a numerical integration of StarkBroadenedLine over the spectral bin. (#366) Bug Fixes: