diff --git a/CHANGELOG.md b/CHANGELOG.md index 4e5a2955..18dfc409 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,9 @@ New: * Add toroidal_mesh_from_polygon for making mesh for not fully-360 degrees axisymmetric elements. (#365) * Add common spectroscopic instruments: Polychromator, SurveySpectrometer, CzernyTurnerSpectrometer. (#299) * Add new classes for free-free Gaunt factors and improve accuracy of the Gaunt factor used in Bremsstrahlung emission model. (#352) +* Add GaussianQuadrature integration method for Function1D. (#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: ---------- 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..c451285f --- /dev/null +++ b/cherab/core/math/integrators/integrators1d.pxd @@ -0,0 +1,39 @@ +# 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, Function2D + + +cdef class Integrator1D: + + cdef: + Function1D function + + cdef double evaluate(self, double a, double b) except? -1e999 + + +cdef class GaussianQuadrature(Integrator1D): + + cdef: + int _min_order, _max_order + double _rtol + ndarray _roots, _weights + double[:] _roots_mv, _weights_mv + + cdef _build_cache(self) diff --git a/cherab/core/math/integrators/integrators1d.pyx b/cherab/core/math/integrators/integrators1d.pyx new file mode 100644 index 00000000..54ec4059 --- /dev/null +++ b/cherab/core/math/integrators/integrators1d.pyx @@ -0,0 +1,223 @@ +# 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 raysect.core.math.function.float cimport autowrap_function1d, Constant1D + +from libc.math cimport INFINITY +cimport cython + + +cdef class Integrator1D: + """ + Compute a definite integral of a one-dimensional function. + + :ivar Function1D integrand: A 1D function to integrate. + """ + + @property + def integrand(self): + """ + A 1D function to integrate. + + :rtype: int + """ + return self.function + + @integrand.setter + def integrand(self, object func not None): + + 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): + """ + Compute an integral of a one-dimensional function over a finite interval + using fixed-tolerance Gaussian quadrature. + (see Scipy `quadrature `). + + :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. + :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, 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.") + + 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 + self._build_cache() + + self.integrand = integrand + + 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_cache() + + @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_cache() + + @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_cache(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) + cdef double evaluate(self, double a, double b) except? -1e999: + """ + Integrates a one-dimensional function over a finite interval. + + :param double a: Lower limit of integration. + :param double b: Upper limit of integration. + + :returns: Gaussian quadrature approximation to integral. + """ + + cdef: + int order, i, ibegin + double newval, oldval, error, x, c, d + + 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 = c + d * self._roots_mv[i] + newval += self._weights_mv[i] * self.function.evaluate(x) + newval *= d + + error = abs(newval - oldval) + oldval = newval + + ibegin += order + + if error < self._rtol * abs(newval): + break + + return newval diff --git a/cherab/core/math/tests/test_integrators.py b/cherab/core/math/tests/test_integrators.py new file mode 100644 index 00000000..fa967165 --- /dev/null +++ b/cherab/core/math/tests/test_integrators.py @@ -0,0 +1,80 @@ +# 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(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 + 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 + 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.""" + quadrature = GaussianQuadrature(relative_tolerance=1.e-8) + a = -0.5 + b = 3. + quadrature.integrand = (2 / sqrt(pi)) * Exp1D(- Arg1D() * Arg1D()) + exact_integral = erf(b) - erf(a) + + self.assertAlmostEqual(quadrature(a, b), exact_integral, places=8) + + +if __name__ == '__main__': + unittest.main() diff --git a/cherab/core/model/lineshape.pxd b/cherab/core/model/lineshape.pxd index 711475b5..2591efe3 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 @@ -41,6 +42,7 @@ 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) diff --git a/cherab/core/model/lineshape.pyx b/cherab/core/model/lineshape.pyx index 5e346175..69278560 100644 --- a/cherab/core/model/lineshape.pyx +++ b/cherab/core/model/lineshape.pyx @@ -19,13 +19,17 @@ # 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 +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 @@ -78,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) @@ -146,14 +148,17 @@ 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.') @@ -289,6 +294,38 @@ 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, _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 self._norm / ((fabs(x - self._x0))**2.5 + self._a) + + cdef class StarkBroadenedLine(LineShapeModel): """ Parametrised Stark broadened line shape based on the Model Microfield Method (MMM). @@ -308,6 +345,9 @@ cdef class StarkBroadenedLine(LineShapeModel): :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 = { @@ -352,7 +392,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, + dict stark_model_coefficients=None, integrator=GaussianQuadrature()): stark_model_coefficients = stark_model_coefficients or self.STARK_MODEL_COEFFICIENTS_DEFAULT @@ -371,7 +412,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 +425,13 @@ 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 ne = self.plasma.get_electron_distribution().density(point.x, point.y, point.z) if ne <= 0.0: @@ -400,6 +443,8 @@ cdef class StarkBroadenedLine(LineShapeModel): lambda_1_2 = self._cij * ne**self._aij / (te**self._bij) + 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 if spectrum.max_wavelength < cutoff_lower_wavelength: @@ -413,28 +458,16 @@ 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 = spectrum.min_wavelength + start * spectrum.delta_wavelength - 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): + upper_wavelength = spectrum.min_wavelength + spectrum.delta_wavelength * (i + 1) - 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.evaluate(lower_wavelength, upper_wavelength) + spectrum.samples_mv[i] += radiance * bin_integral / spectrum.delta_wavelength lower_wavelength = upper_wavelength - lower_value = upper_value - - # 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__':