diff --git a/CHANGELOG.md b/CHANGELOG.md index dbe69a21..d6d704ff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ New: * Add group observer class for each of Raysect's 0D observers. (#332) * Add a demo for observer group handling and plotting. * Add verbose parameter to SartOpencl solver (default is False). (#358) +* Add Thomson Scattering model. (#97) * Add Generomak core plasma profiles. (#360) * Add toroidal_mesh_from_polygon for making mesh for not fully-360 degrees axisymmetric elements. (#365) diff --git a/cherab/core/laser/__init__.pxd b/cherab/core/laser/__init__.pxd new file mode 100644 index 00000000..2760a2ea --- /dev/null +++ b/cherab/core/laser/__init__.pxd @@ -0,0 +1,4 @@ +from cherab.core.laser.node cimport Laser +from cherab.core.laser.model cimport LaserModel +from cherab.core.laser.laserspectrum cimport LaserSpectrum +from cherab.core.laser.profile cimport LaserProfile \ No newline at end of file diff --git a/cherab/core/laser/__init__.py b/cherab/core/laser/__init__.py new file mode 100644 index 00000000..911a64b1 --- /dev/null +++ b/cherab/core/laser/__init__.py @@ -0,0 +1,4 @@ +from .node import Laser +from .model import LaserModel +from .laserspectrum import LaserSpectrum +from .profile import LaserProfile \ No newline at end of file diff --git a/cherab/core/laser/laserspectrum.pxd b/cherab/core/laser/laserspectrum.pxd new file mode 100644 index 00000000..a00945eb --- /dev/null +++ b/cherab/core/laser/laserspectrum.pxd @@ -0,0 +1,47 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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 cimport Function1D + +from cherab.core.utility.constants cimport SPEED_OF_LIGHT, PLANCK_CONSTANT + +from numpy cimport ndarray + + +cdef class LaserSpectrum(Function1D): + + cdef: + double _min_wavelength, _max_wavelength, _delta_wavelength + int _bins + ndarray _power, _power_spectral_density, _wavelengths # power_spectral_density [w/nm] + double[::1] power_mv, power_spectral_density_mv, wavelengths_mv + + cpdef double evaluate_integral(self, double lower_limit, double upper_limit) + + cpdef void _update_cache(self) + + cpdef double get_min_wavelenth(self) + + cpdef double get_max_wavelenth(self) + + cpdef int get_spectral_bins(self) + + cpdef double get_delta_wavelength(self) + + cpdef double _get_bin_power_spectral_density(self, double wavelength_lower, double wavelength_upper) \ No newline at end of file diff --git a/cherab/core/laser/laserspectrum.pyx b/cherab/core/laser/laserspectrum.pyx new file mode 100644 index 00000000..e196fba1 --- /dev/null +++ b/cherab/core/laser/laserspectrum.pyx @@ -0,0 +1,192 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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 cimport Function1D +from raysect.optical cimport Point3D, Vector3D + +from cherab.core.utility import Notifier +from cherab.core.utility.constants cimport SPEED_OF_LIGHT, PLANCK_CONSTANT + +import numpy as np +cimport numpy as np + + +cdef class LaserSpectrum(Function1D): + """ + Laser spectrum base class. + + This is an abstract class and cannot be used for observing. + + A 1D function holding information about the spectral properties + of a laser. The scattered spectrum is calculated as an iteration + over the laser spectrum. + + + .. warning:: + When adding a LaserSpectrum, a special care should be given + to the integral power of the laser spectrum. During the + scattering calculation, the spectral power can be multiplied + by the power spatial distribution [W * m ** -3] of the laser + power from the LaserProfile. If the integral power + of the LaserSpectrum is not 1, unexpected values + might be obtained. + + .. note:: + It is expected that majority of the fusion applications can + neglect the influence of the spectral shape of the + laser and can use laser spectrum with a single + bin, which approximates an infinitely narrow laser spectrum. + + :param float min_wavelength: The minimum wavelength of the laser + spectrum in nm. + :param float max_wavelength: The maximum wavelength of the laser + spectrum in nm. + :param int bins: The number of spectral bins of the laser spectrum. + :ivar float min_wavelength: The minimum wavelength of the laser + spectrum in nm. + :ivar float max_wavelength: The maximum wavelength of the laser + spectrum in nm. + :ivar int bins: The number of specral bins of the laser spectrum + :ivar ndarray wavelengths: The wavelengt coordinate vector in nm. + :ivar ndarray power_spectral_density: The values of the power + spectral density in W / nm. + :ivar float delta_wavelength: Spectral width of the bins in nm. + """ + + def __init__(self, double min_wavelength, double max_wavelength, int bins): + + super().__init__() + + self._check_wavelength_validity(min_wavelength, max_wavelength) + + self._min_wavelength = min_wavelength + self._max_wavelength = max_wavelength + + self.bins = bins + + @property + def min_wavelength(self): + return self._min_wavelength + + @min_wavelength.setter + def min_wavelength(self, double value): + + self._check_wavelength_validity(value, self.max_wavelength) + self._min_wavelength = value + self._update_cache() + + @property + def max_wavelength(self): + return self._max_wavelength + + @max_wavelength.setter + def max_wavelength(self, double value): + + self._check_wavelength_validity(self.min_wavelength, value) + self._max_wavelength = value + self._update_cache() + + @property + def bins(self): + return self._bins + + @bins.setter + def bins(self, int value): + if value <= 0: + raise ValueError("Value has to be larger than 0") + + self._bins = value + self._update_cache() + + @property + def wavelengths(self): + return self._wavelengths + + @property + def power_spectral_density(self): + return self._power_spectral_density + + @property + def delta_wavelength(self): + return self._delta_wavelength + + def _check_wavelength_validity(self, min_wavelength, max_wavelength): + + if min_wavelength <= 0: + raise ValueError("min_wavelength has to be larger than 0, but {} passed.".format(min_wavelength)) + if max_wavelength <= 0: + raise ValueError("min_wavelength has to be larger than 0, but {} passed.".format(max_wavelength)) + + if min_wavelength >= max_wavelength: + raise ValueError("min_wavelength has to be smaller than max_wavelength: min_wavelength={} > max_wavelength={}".format(min_wavelength, max_wavelength)) + + cpdef double get_min_wavelenth(self): + return self._min_wavelength + + cpdef double get_max_wavelenth(self): + return self._min_wavelength + + cpdef int get_spectral_bins(self): + return self._bins + + cpdef double get_delta_wavelength(self): + return self._delta_wavelength + + cpdef void _update_cache(self): + + cdef: + Py_ssize_t index + double delta_wvl_half, wvl_lower, wvl_upper, wvl + + self._delta_wavelength = (self._max_wavelength - self._min_wavelength) / self._bins + self._wavelengths = np.zeros(self.bins, dtype=np.double) + self.wavelengths_mv = self._wavelengths + + for index in range(self._bins): + self._wavelengths[index] = self._min_wavelength + (0.5 + index) * self._delta_wavelength + + self._power_spectral_density = np.zeros(self._bins, dtype=np.double) # power spectral density (PSD) + self.power_spectral_density_mv = self._power_spectral_density + + self._power = np.zeros(self._bins, dtype=np.double) # power in a spectral bin (PSD * delta wavelength) + self.power_mv = self._power + + delta_wvl_half = self._delta_wavelength * 0.5 + wvl_lower = self.wavelengths_mv[0] - delta_wvl_half + + for index in range(self._bins): + wvl = wvl_lower + delta_wvl_half + wvl_upper = wvl_lower + self._delta_wavelength + + self.power_spectral_density_mv[index] = self._get_bin_power_spectral_density(wvl_lower, wvl_upper) + self.power_mv[index] = self.power_spectral_density_mv[index] * self._delta_wavelength # power in the spectral bin for scattering calculations + + wvl_lower = wvl_upper + + cpdef double evaluate_integral(self, double lower_limit, double upper_limit): + raise NotImplementedError('Virtual method must be implemented in a sub-class.') + + cpdef double _get_bin_power_spectral_density(self, double wavelength_lower, double wavelength_upper): + """ + Returns the power spectral density in a bin. + + This method can be overidden if a better precision is needed. + For example for distributions with known cumulative distribution function. + """ + return 0.5 * (self.evaluate(wavelength_lower) + self.evaluate(wavelength_upper)) diff --git a/cherab/core/laser/material.pxd b/cherab/core/laser/material.pxd new file mode 100644 index 00000000..c91fa416 --- /dev/null +++ b/cherab/core/laser/material.pxd @@ -0,0 +1,31 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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.scenegraph._nodebase cimport _NodeBase +from raysect.core.math cimport AffineMatrix3D +from raysect.optical.material.emitter cimport InhomogeneousVolumeEmitter + +from cherab.core.laser.node cimport Laser + + +cdef class LaserMaterial(InhomogeneousVolumeEmitter): + + cdef: + AffineMatrix3D _laser_to_plasma, _laser_segment_to_laser_node + list _models diff --git a/cherab/core/laser/material.pyx b/cherab/core/laser/material.pyx new file mode 100644 index 00000000..4732e01a --- /dev/null +++ b/cherab/core/laser/material.pyx @@ -0,0 +1,66 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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.scenegraph._nodebase cimport _NodeBase +from raysect.optical cimport World, Primitive, Ray, Spectrum, Point3D, Vector3D, AffineMatrix3D +from raysect.optical.material.emitter cimport InhomogeneousVolumeEmitter +from raysect.optical.material.emitter.inhomogeneous cimport VolumeIntegrator + +from cherab.core.laser.node cimport Laser +from cherab.core.laser.model cimport LaserModel + + +cdef class LaserMaterial(InhomogeneousVolumeEmitter): + + def __init__(self, Laser laser not None, _NodeBase laser_segment not None, list models, VolumeIntegrator integrator not None): + + super().__init__(integrator) + + self._laser_segment_to_laser_node = laser_segment.to(laser) + self._laser_to_plasma = laser_segment.to(laser.plasma) + self.importance = laser.importance + + #validate and set models + for model in models: + if not isinstance(model, LaserModel): + raise TypeError("Model supplied to laser are not LaserMaterial is not LaserModel") + model.plasma = laser.plasma + model.laser_profile = laser.laser_profile + model.laser_spectrum = laser.laser_spectrum + + self._models = models + + cpdef Spectrum emission_function(self, Point3D point, Vector3D direction, Spectrum spectrum, + World world, Ray ray, Primitive primitive, + AffineMatrix3D to_local, AffineMatrix3D to_world): + + cdef: + Point3D point_plasma, point_laser + Vector3D direction_plasma, direction_laser + LaserModel model + + point_laser = point.transform(self._laser_segment_to_laser_node) + direction_laser = direction.transform(self._laser_segment_to_laser_node) # observation vector in the laser frame + point_plasma = point.transform(self._laser_to_plasma) + direction_plasma = direction.transform(self._laser_to_plasma) + + for model in self._models: + spectrum = model.emission(point_plasma, direction_plasma, point_laser, direction_laser, spectrum) + + return spectrum diff --git a/cherab/core/laser/model.pxd b/cherab/core/laser/model.pxd new file mode 100644 index 00000000..71cdf9a8 --- /dev/null +++ b/cherab/core/laser/model.pxd @@ -0,0 +1,37 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + + +from raysect.optical cimport Vector3D, Point3D +from raysect.optical.spectrum cimport Spectrum + +from cherab.core cimport Plasma +from cherab.core.laser.profile cimport LaserProfile +from cherab.core.laser.laserspectrum cimport LaserSpectrum + + +cdef class LaserModel: + cdef: + Plasma _plasma + LaserSpectrum _laser_spectrum + LaserProfile _laser_profile + + cpdef Spectrum emission(self, Point3D point_plasma, Vector3D observation_plasma, Point3D point_laser, + Vector3D observation_laser, Spectrum spectrum) + + cdef object __weakref__ diff --git a/cherab/core/laser/model.pyx b/cherab/core/laser/model.pyx new file mode 100644 index 00000000..c7d76cae --- /dev/null +++ b/cherab/core/laser/model.pyx @@ -0,0 +1,78 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + + +from raysect.optical cimport Vector3D, Point3D +from raysect.optical.spectrum cimport Spectrum + +from cherab.core cimport Plasma +from cherab.core.laser.profile cimport LaserProfile +from cherab.core.laser.laserspectrum cimport LaserSpectrum + + +cdef class LaserModel: + """ + Laser spectrum base class. + + This is an abstract class and cannot be used for observing. + + Calculates the contribution to a spectrum caused by a laser. + + :param laser_profile: LaserProfile object + :param plasma: Plasma object + :param laser_spectrum: LaserSpectrum object + + :ivar laser_profile: LaserProfile object + :ivar plasma: Plasma object + :ivar laser_spectrum: LaserSpectrum object + """ + def __init__(self, LaserProfile laser_profile=None, LaserSpectrum laser_spectrum=None, Plasma plasma=None): + + self._laser_profile = laser_profile + self._laser_spectrum = laser_spectrum + self._plasma = plasma + + cpdef Spectrum emission(self, Point3D point_plasma, Vector3D observation_plasma, Point3D point_laser, Vector3D observation_laser, + Spectrum spectrum): + + raise NotImplementedError('Virtual method must be implemented in a sub-class.') + + @property + def laser_profile(self): + return self._laser_profile + + @laser_profile.setter + def laser_profile(self, LaserProfile value): + self._laser_profile = value + + @property + def plasma(self): + return self._plasma + + @plasma.setter + def plasma(self, Plasma value): + self._plasma = value + + @property + def laser_spectrum(self): + return self._laser_spectrum + + @laser_spectrum.setter + def laser_spectrum(self, LaserSpectrum value): + + self._laser_spectrum = value diff --git a/cherab/core/laser/node.pxd b/cherab/core/laser/node.pxd new file mode 100644 index 00000000..6045fb4a --- /dev/null +++ b/cherab/core/laser/node.pxd @@ -0,0 +1,55 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + + +from raysect.optical cimport Point3D, Vector3D, Node, Spectrum, Primitive +from raysect.optical.material.emitter.inhomogeneous cimport VolumeIntegrator +from raysect.primitive cimport Cylinder + +from cherab.core.plasma cimport Plasma +from cherab.core.laser.profile cimport LaserProfile +from cherab.core.laser.laserspectrum cimport LaserSpectrum +from cherab.core.laser.model cimport LaserModel + + +cdef class ModelManager: + + cdef: + list _models + readonly object notifier + + cpdef object set(self, object models) + + cpdef object add(self, LaserModel model) + + cpdef object clear(self) + + +cdef class Laser(Node): + + cdef: + readonly object notifier + double _importance + Plasma _plasma + ModelManager _models + LaserProfile _laser_profile + LaserSpectrum _laser_spectrum + list _geometry + VolumeIntegrator _integrator + + cdef object __weakref__ \ No newline at end of file diff --git a/cherab/core/laser/node.pyx b/cherab/core/laser/node.pyx new file mode 100644 index 00000000..b0b821b5 --- /dev/null +++ b/cherab/core/laser/node.pyx @@ -0,0 +1,261 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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.primitive cimport Cylinder +from raysect.optical cimport World, AffineMatrix3D, Primitive, Ray +from raysect.optical.material.emitter.inhomogeneous cimport NumericalIntegrator +from raysect.core cimport translate, Material + +from cherab.core.laser.material cimport LaserMaterial +from cherab.core.laser.model cimport LaserModel +from cherab.core.laser.profile import LaserProfile +from cherab.core.laser.laserspectrum import LaserSpectrum +from cherab.core.utility import Notifier +from libc.math cimport M_PI + +from math import ceil + +cdef double DEGREES_TO_RADIANS = (M_PI / 180) + + +cdef class ModelManager: + + def __init__(self): + self._models = [] + self.notifier = Notifier() + + def __iter__(self): + return iter(self._models) + + cpdef object set(self, object models): + + # copy models and test it is an iterable + models = list(models) + + # check contents of list are laser models + for model in models: + if not isinstance(model, LaserModel): + raise TypeError('The model list must consist of only LaserModel objects.') + + self._models = models + self.notifier.notify() + + cpdef object add(self, LaserModel model): + + if not model: + raise ValueError('Model must not be None type.') + + self._models.append(model) + self.notifier.notify() + + cpdef object clear(self): + self._models = [] + self.notifier.notify() + + +cdef class Laser(Node): + """ + A scene-graph object representing a laser of laser light. + + The Cherab laser object holds basic information about the laser and connects + the components which are needed for the laser description. With specified + emission models it can contribute to observed radiation. + + The Laser object is a Raysect scene-graph node and lives in it's own + coordinate space. This coordinate space is defined relative to it's parent + scene-graph object by an AffineTransform. The beam parameters are defined + in the Laser object coordinate space. Models using the beam object must + convert any spatial coordinates into beam space before requesting values + from the Laser object. + + The main physical properties of the laser are defined by the three + attributes laser_spectrum, laser_profile and models. The laser_spectrum + has to be an instance of LaserSpectrum and defines the spectral properties + of the laser light. The laser_profile has to be an instance of LaserProfile + and it holds all the space related definitions as volumetric distribution + of laser light energy polarisation direction. In the models a list of LaserModels + can be stored, which calculate the contribution of the laser ligth to the observed + radiation. The models can cover various applications as for example + Thomson scattering. Please see the documentation of individual classes + for more detail. + + The shape of the laser (e.g. cylinder) and its parameters (e.g. radius) + is controled by the LaserProfile. + + The plasma reference has to be specified to attach the any models. + + :param Node parent: The parent node in the Raysect scene-graph. + See the Raysect documentation for more guidance. + :param AffineMatrix3D transform: The transform defining the spatial position + and orientation of this laser. See the Raysect documentation if you need + guidance on how to use AffineMatrix3D transforms. + :param str name: The name for this laser object. + :ivar Plasma plasma: The plasma instance with which this laser interacts. + :ivar float importance: The importance sampling factor. + :ivar LaserSpectrum laser_spectrum: The LaserSpectrum instance with which this laser interacts. + :ivar LaserProfile laser_profile: The LaserProfile instance with which this laser interacts. + :ivar ModelManager models: The manager class that sets and provides access to the + emission models for this laser. + :ivar VolumeIntegrator integrator: The configurable method for doing + volumetric integration through the laser along a Ray's path. Defaults to + a numerical integrator with 1mm step size, NumericalIntegrator(step=0.001). + """ + + def __init__(self, object parent=None, AffineMatrix3D transform=None, str name=None): + + super().__init__(parent, transform, name) + + self._set_init_values() + + self.notifier = Notifier() + + self._models = ModelManager() + + self._integrator = NumericalIntegrator(step=1e-3) + + self._importance = 1. + + def _set_init_values(self): + """ + Sets initial values of the laser shape to avoid errors. + """ + self._importance = 0. + self._geometry = [] + + @property + def plasma(self): + return self._plasma + + @plasma.setter + def plasma(self, Plasma value not None): + + #unregister from old plasma notifier + if self._plasma is not None: + self._plasma.notifier.remove(self._plasma_changed) + + self._plasma = value + self._plasma.notifier.add(self._plasma_changed) + + self._configure_materials() + + @property + def importance(self): + return self._importance + + @importance.setter + def importance(self, double value): + + self._importance = value + self._configure_materials() + + @property + def laser_spectrum(self): + return self._laser_spectrum + + @laser_spectrum.setter + def laser_spectrum(self, LaserSpectrum value): + self._laser_spectrum = value + self._configure_materials() + + @property + def laser_profile(self): + return self._laser_profile + + @laser_profile.setter + def laser_profile(self, LaserProfile value): + + if self._laser_profile is not None: + self._laser_profile.notifier.remove(self.configure_geometry) + + self._laser_profile = value + self._laser_profile.notifier.add(self.configure_geometry) + + self.configure_geometry() + + @property + def models(self): + return list(self._models) + + @models.setter + def models(self, value): + + # check necessary data is available + if not all([self._plasma, self._laser_profile, self._laser_spectrum]): + raise ValueError("The plasma, laser_profile and laser_spectrum must be set before before specifying any models.") + + self._models.set(value) + self._configure_materials() + + @property + def integrator(self): + return self._integrator + + @integrator.setter + def integrator(self, VolumeIntegrator value): + self._integrator = value + + for i in self._geometry: + i.material.integrator = value + + def configure_geometry(self): + """ + Reconfigure the laser primitives and materials. + """ + + self._build_geometry() + self._configure_materials() + + def _build_geometry(self): + """ + Delete and build new laser segments + """ + # remove old laser segments in any case + for i in self._geometry: + i.parent = None + self._geometry = [] + + # no point in adding segments if there is no model and profile + if self._laser_profile is None: + return + + # rebuild geometry + self._geometry = self._laser_profile.generate_geometry() + + for i in self._geometry: + i.parent = self + + def _configure_materials(self): + """ + Configure laser segment materials + """ + if not list(self._models) or self._plasma is None or self._laser_spectrum is None: + return + + for i in self._geometry: + i.material = LaserMaterial(self, i, list(self._models), self._integrator) + + def get_geometry(self): + return self._geometry + + def _plasma_changed(self): + """React to change of plasma and propagate the information.""" + self._configure_materials() + + def _modified(self): + self._configure_materials() diff --git a/cherab/core/laser/profile.pxd b/cherab/core/laser/profile.pxd new file mode 100644 index 00000000..c634a2f8 --- /dev/null +++ b/cherab/core/laser/profile.pxd @@ -0,0 +1,41 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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 cimport Function3D +from raysect.core.math.function.vector3d cimport Function3D as VectorFunction3D + +from raysect.optical cimport Spectrum, Point3D, Vector3D + +from cherab.core.laser.node cimport Laser + + +cdef class LaserProfile: + + cdef: + VectorFunction3D _polarization3d, _pointing3d + Function3D _energy_density3d + readonly object notifier + + cpdef Vector3D get_pointing(self, double x, double y, double z) + + cpdef Vector3D get_polarization(self, double x, double y, double z) + + cpdef double get_energy_density(self, double x, double y, double z) + + cpdef list generate_geometry(self) \ No newline at end of file diff --git a/cherab/core/laser/profile.pyx b/cherab/core/laser/profile.pyx new file mode 100644 index 00000000..6356d1cb --- /dev/null +++ b/cherab/core/laser/profile.pyx @@ -0,0 +1,163 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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 cimport Function3D +from raysect.core.math.function.vector3d cimport Function3D as VectorFunction3D + +from raysect.optical cimport SpectralFunction, Spectrum, InterpolatedSF, Point3D, Vector3D + +from cherab.core.laser.node cimport Laser +from cherab.core.utility import Notifier + + +cdef class LaserProfile: + """ + LaserProfile base class. + + This is an abstract class and cannot be used for observing. + + Provides information about spatial properties of the laser beam: + direction of the laser propagation (direction + of the Poynting vector), polarisation of the ligth as the direction + of the electric component vector and volumetric energy density of + the laser light. + + All the laser properties are evaluated in the frame of reference of + the laser. + + .. warning:: + When combining a LaserProfile with a LaserSpectrum for a laser, + a special care has to be given to obtain the correct power + of the scattered spectrum. Scattering models can multiply + both the spectral power density given by the LaserProfile and + the volumetric energy density given by the LaserProfile. + Combination of incompatible cases may yield incorrect + values of scattered power. + + :ivar Laser laser: The Laser scenegraph node the LaserProfile + is connected to. + """ + + def __init__(self): + + self.notifier = Notifier() + + def set_polarization_function(self, VectorFunction3D function): + """ + Assigns the 3D vector function describing the polarisation vector. + + The polarisation is given as the direction of the electric + component of the electromagnetic wave. + + The function is specified in the laser space. + + :param VectorFunction3D function: A 3D vector function describing + the polarisation vector. + """ + self._polarization3d = function + + def set_pointing_function(self, VectorFunction3D function): + """ + Assings the 3D vector function describing the direction of the laser propagation. + + The direction of the laser light propagation is the direction + of the Poynting vector. + + :param VectorFunction3D function: A 3D vector function describing + the laser light propagation direction + """ + self._pointing3d = function + + def set_energy_density_function(self, Function3D function): + """ + Assigns the 3D scalar function describing the laser energy distribution. + + The laser power distribution is the value of the volumetric + energy density of the laser light. + """ + self._energy_density3d = function + + cpdef Vector3D get_pointing(self, double x, double y, double z): + """ + Returns the laser light propagation direction. + + At the point (x, y, z) in the laser space. + + :param x: x coordinate in meters. + :param y: y coordinate in meters. + :param z: z coordinate in meters. + :return: Intensity in m^-3. + """ + + return self._pointing3d.evaluate(x, y, z) + + cpdef Vector3D get_polarization(self, double x, double y, double z): + """ + Returns a vector denoting the laser polarisation. + + The polarisation direction is the direction of the electric + component of the electromagnetic wave for the point (x, y, z) + in the laser space. + + :param x: x coordinate in meters. + :param y: y coordinate in meters. + :param z: z coordinate in meters. + :return: power density in Wm^-3. + """ + + return self._polarization3d(x, y, z) + + cpdef double get_energy_density(self, double x, double y, double z): + """ + Returns the volumetric energy density of the laser light in W*m^-3. + + At the point (x, y, z) in the laser space. + + :param x: x coordinate in meters in the laser frame. + :param y: y coordinate in meters in the laser frame. + :param z: z coordinate in meters in the laser frame. + :return: power density in W*m^-3. + """ + + return self._energy_density3d.evaluate(x, y, z) + + cpdef list generate_geometry(self): + """ + returns list of raysect primitives composing the laser geometry + + This method is called from the Laser instance to which the instance + of Profile is attached to. The Laser instance will be assigned as + the parent to the returned primitives in the Laser._configure method. + The Laser._configure method does not change any transforms. This is + why the returned primitives have to have their transforms already + initialised in the frame of the laser, when returned. + """ + + raise NotImplementedError("Virtual function density not defined.") + + def _change(self): + """ + Called if the laser properties change. + + If the model caches calculation data that would be invalidated if its + source data changes then this method may be overridden to clear the + cache. + """ + + pass diff --git a/cherab/core/laser/tests/__init__.py b/cherab/core/laser/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/cherab/core/laser/tests/test_laser.py b/cherab/core/laser/tests/test_laser.py new file mode 100644 index 00000000..f9f85732 --- /dev/null +++ b/cherab/core/laser/tests/test_laser.py @@ -0,0 +1,110 @@ +import unittest + +from raysect.optical import World +from raysect.optical.material.emitter.inhomogeneous import NumericalIntegrator + +from cherab.core import Plasma +from cherab.core.laser.node import Laser +from cherab.core.model.laser.laserspectrum import ConstantSpectrum +from cherab.core.model.laser.model import SeldenMatobaThomsonSpectrum +from cherab.core.model.laser.profile import UniformEnergyDensity + + +class TestLaser(unittest.TestCase): + + def test_laser_init(self): + """ + Test correct initialisation of a laser instance. + """ + + world = World() + laser = Laser(parent=world) + + with self.assertRaises(ValueError, msg="Model was attached before Plasma, Profile and LaserSpectrum were specified."): + laser.models = [SeldenMatobaThomsonSpectrum()] + + laser.laser_profile = UniformEnergyDensity() + with self.assertRaises(ValueError, msg="Model was attached before Plasma, Profile and LaserSpectrum were specified."): + laser.models = [SeldenMatobaThomsonSpectrum()] + + laser.laser_spectrum = ConstantSpectrum(min_wavelength=1059, max_wavelength=1061, bins=10) + with self.assertRaises(ValueError, msg="Model was attached before Plasma, Profile and LaserSpectrum were specified."): + laser.models = [SeldenMatobaThomsonSpectrum()] + + laser.plasma = Plasma(parent=world) + laser.models = [SeldenMatobaThomsonSpectrum()] + + def test_reference_change(self): + + world = World() + + laser_profile = UniformEnergyDensity(laser_length=1, laser_radius=0.1) + laser_spectrum = ConstantSpectrum(min_wavelength=1059, max_wavelength=1061, bins=10) + plasma = Plasma(parent=world) + models = [SeldenMatobaThomsonSpectrum()] + + laser_profile2 = UniformEnergyDensity() + laser_spectrum2 = ConstantSpectrum(min_wavelength=1059, max_wavelength=1061, bins=10) + plasma2 = Plasma(parent=world) + models2 = [SeldenMatobaThomsonSpectrum()] + + laser = Laser(parent=world) + + laser.laser_spectrum = laser_spectrum + laser.plasma = plasma + laser.laser_profile = laser_profile + laser.models = models + + for mod in list(laser.models): + self.assertIs(mod.laser_profile, laser_profile, msg="laser_profile reference in emission model" + "is not set correctly.") + self.assertIs(mod.plasma, plasma, msg="plasma reference in emission model" + "is not set correctly.") + self.assertIs(mod.laser_spectrum, laser_spectrum, msg="laser_spectrum reference in emission model" + "is not set correctly.") + + laser.laser_spectrum = laser_spectrum2 + laser.plasma = plasma2 + laser.laser_profile = laser_profile2 + + for mod in list(laser.models): + self.assertIs(mod.laser_profile, laser_profile2, msg="laser_profile reference in emission model" + "is not set correctly.") + self.assertIs(mod.plasma, plasma2, msg="plasma reference in emission model" + "is not set correctly.") + self.assertIs(mod.laser_spectrum, laser_spectrum2, msg="laser_spectrum reference in emission model" + "is not set correctly.") + + laser.models = models + models2 + + for mod in list(laser.models): + self.assertIs(mod.laser_profile, laser_profile2, msg="laser_profile reference in emission model" + "is not set correctly.") + self.assertIs(mod.plasma, plasma2, msg="plasma reference in emission model" + "is not set correctly.") + self.assertIs(mod.laser_spectrum, laser_spectrum2, msg="laser_spectrum reference in emission model" + "is not set correctly.") + + def test_integrator_change(self): + + world = World() + + laser_profile = UniformEnergyDensity(laser_length=1, laser_radius=0.1) + laser_spectrum = ConstantSpectrum(min_wavelength=1059, max_wavelength=1061, bins=10) + plasma = Plasma(parent=world) + models = [SeldenMatobaThomsonSpectrum()] + + laser = Laser(parent=world) + + laser.laser_spectrum = laser_spectrum + laser.plasma = plasma + laser.laser_profile = laser_profile + laser.models = models + + integrator = NumericalIntegrator(1e-4) + + laser.integrator = integrator + + for i in laser.get_geometry(): + self.assertIs(i.material.integrator, integrator, msg="Integrator not updated properly") + diff --git a/cherab/core/laser/tests/test_laserspectrum.py b/cherab/core/laser/tests/test_laserspectrum.py new file mode 100644 index 00000000..9473e43d --- /dev/null +++ b/cherab/core/laser/tests/test_laserspectrum.py @@ -0,0 +1,62 @@ +import unittest +import numpy as np + +from cherab.core.laser.laserspectrum import LaserSpectrum +from raysect.optical.spectrum import Spectrum + +class TestLaserSpectrum(unittest.TestCase): + + def test_laserspectrum_init(self): + # test min_wavelength boundaries + with self.assertRaises(ValueError, + msg="LaserSpectrum did not raise a ValueError with min_wavelength being zero."): + LaserSpectrum(0., 100, 200) + LaserSpectrum(-1, 100, 200) + + # test max_wavelength boundaries + with self.assertRaises(ValueError, + msg="LaserSpectrum did not raise a ValueError with max_wavelength being zero."): + LaserSpectrum(10, 0, 200) + LaserSpectrum(10, -1, 200) + + # test min_wavelength >= max_wavelength + with self.assertRaises(ValueError, + msg="LaserSpectrum did not raise a ValueError with max_wavelength < min_wavelength."): + LaserSpectrum(40, 30, 200) + LaserSpectrum(30, 30, 200) + + # test bins > 0 + with self.assertRaises(ValueError, + msg="LaserSpectrum did not raise a ValueError with max_wavelength < min_wavelength."): + LaserSpectrum(30, 30, 0) + LaserSpectrum(30, 30, -1) + + def test_laserspectrum_changes(self): + laser_spectrum = LaserSpectrum(100, 200, 100) + + # change min_wavelength to be larger or equal to max_wavelength + with self.assertRaises(ValueError, + msg="LaserSpectrum did not raise ValueError for min_wavelength change " + "with min_wavelength >= max_wavelength."): + laser_spectrum.min_wavelength = 300 + laser_spectrum.min_wavelength = 200 + + # change max_wavelength to be smaller than max_wavelength + with self.assertRaises(ValueError, + msg="LaserSpectrum did not raise ValueError for max_wavelength change " + "with min_wavelength > max_wavelength."): + laser_spectrum.max_wavelength = 50 + laser_spectrum.max_wavelength = 100 + + with self.assertRaises(ValueError, + msg="LaserSpectrum did not raise a ValueError with max_wavelength < min_wavelength."): + laser_spectrum.bins = -1 + laser_spectrum.bins = 0 + + # laser spectrum should have same behaviour as Spectrum from raysect.optical + spectrum = Spectrum(laser_spectrum.min_wavelength, laser_spectrum.max_wavelength, laser_spectrum.bins) + + # test caching of spectrum data, behaviour should be consistent with raysect.optical.spectrum.Spectrum + self.assertTrue(np.array_equal(laser_spectrum.wavelengths, spectrum.wavelengths), + "LaserSpectrum.wavelengths values are not equal to Spectrum.wavelengths " + "with same boundaries and number of bins") \ No newline at end of file diff --git a/cherab/core/model/laser/__init__.pxd b/cherab/core/model/laser/__init__.pxd new file mode 100644 index 00000000..3b0e3818 --- /dev/null +++ b/cherab/core/model/laser/__init__.pxd @@ -0,0 +1,4 @@ +from cherab.core.model.laser.laserspectrum import ConstantSpectrum, GaussianSpectrum +from cherab.core.model.laser.profile import UniformEnergyDensity, ConstantAxisymmetricGaussian +from cherab.core.model.laser.profile import ConstantBivariateGaussian, TrivariateGaussian, GaussianBeamAxisymmetric +from cherab.core.model.laser.model import SeldenMatobaThomsonSpectrum \ No newline at end of file diff --git a/cherab/core/model/laser/__init__.py b/cherab/core/model/laser/__init__.py new file mode 100644 index 00000000..930b55a7 --- /dev/null +++ b/cherab/core/model/laser/__init__.py @@ -0,0 +1,4 @@ +from .laserspectrum import ConstantSpectrum, GaussianSpectrum +from .profile import UniformEnergyDensity, ConstantAxisymmetricGaussian +from .profile import ConstantBivariateGaussian, TrivariateGaussian, GaussianBeamAxisymmetric +from .model import SeldenMatobaThomsonSpectrum diff --git a/cherab/core/model/laser/laserspectrum.pxd b/cherab/core/model/laser/laserspectrum.pxd new file mode 100644 index 00000000..e21d36fb --- /dev/null +++ b/cherab/core/model/laser/laserspectrum.pxd @@ -0,0 +1,34 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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.laser cimport LaserSpectrum + + +cdef class ConstantSpectrum(LaserSpectrum): + + cdef double evaluate(self, double x) except? -1e999 + + +cdef class GaussianSpectrum(LaserSpectrum): + + cdef: + double _stddev, _recip_stddev, _normalisation, _mean + double _norm_cdf + + cdef double evaluate(self, double x) except? -1e999 diff --git a/cherab/core/model/laser/laserspectrum.pyx b/cherab/core/model/laser/laserspectrum.pyx new file mode 100644 index 00000000..58c1cada --- /dev/null +++ b/cherab/core/model/laser/laserspectrum.pyx @@ -0,0 +1,133 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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.laser cimport LaserSpectrum +from libc.math cimport sqrt, exp, M_PI, erf, M_SQRT2 + + +cdef class ConstantSpectrum(LaserSpectrum): + """ + A laser spectrum with constant power. + + Has a constant, non-zero distribution of power spectral density + between the min_wavelength and max_wavelength. The integral value + of the power is 1 W. + + .. note:: + The ConstantSpectrum class is suitable for approximation + of an infinitely thin laser spectrum, e.g.: + ConstantSpectrum(1063.9, 1064.1, 1) + """ + + def __init__(self, double min_wavelength, double max_wavelength, int bins): + + super().__init__(min_wavelength, max_wavelength, bins) + + cdef double evaluate(self, double x) except? -1e999: + """ + Returns the spectral power density for the given wavelength. + + :param float x: Wavelength in nm. + + :return: Power spectral density in W/nm. + """ + + cdef: + double spectrum_width + int index + + if self._min_wavelength <= x <= self._max_wavelength: + return 1.0 / (self._max_wavelength - self._min_wavelength) + else: + return 0 + + +cdef class GaussianSpectrum(LaserSpectrum): + """ + A laser spectrum with a normally distributed power spectral density. + + Has a Gaussian-like spectral shape. The inegral value of power is 1 W. + + :param float mean: The mean value of the Gaussian distribution + of the laser spectrum in nm, can be thought of as the central + wavelength of the laser. + :param float stddev: Standard deviation of the Gaussian + distribution of the laser spectrum. + + :ivar float stddev: Standard deviation of the Gaussian + distribution of the laser spectrum. + :ivar float mean: The mean value of the Gaussian distribution + of the laser spectrum in nm, can be thought of as the central + wavelength of the laser. + """ + + def __init__(self, double min_wavelength, double max_wavelength, int bins, double mean, double stddev): + + self.stddev = stddev + self.mean = mean + super().__init__(min_wavelength, max_wavelength, bins) + + @property + def stddev(self): + return self._stddev + + @stddev.setter + def stddev(self, value): + if value <= 0: + raise ValueError("Value has to be larger than 0") + + self._stddev = value + self._recip_stddev = 1 / value + self._normalisation = 1 / (value * sqrt(2 * M_PI)) + self._norm_cdf = 1 / (value * M_SQRT2) + + @property + def mean(self): + return self._mean + + @mean.setter + def mean(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0") + + self._mean = value + + cdef double evaluate(self, double x) except? -1e999: + """ + Returns the spectral power density for the given wavelength. + + :param float x: Wavelength in nm. + + :return: Power spectral density in W/nm. + """ + return self._normalisation * exp(-0.5 * ((x - self._mean) * self._recip_stddev) ** 2) + + cpdef double _get_bin_power_spectral_density(self, double wavelength_lower, double wavelength_upper): + """ + Returns the power spectral density in a bin. + + Overrides the parent method to deliver better precision. + """ + + cdef: + double val_lower, val_upper + + val_lower = erf((wavelength_lower - self._mean) * self._norm_cdf) + val_upper = erf((wavelength_upper - self._mean) * self._norm_cdf) + return 0.5 * (val_upper - val_lower) / self._delta_wavelength \ No newline at end of file diff --git a/cherab/core/model/laser/math_functions.pxd b/cherab/core/model/laser/math_functions.pxd new file mode 100644 index 00000000..08827094 --- /dev/null +++ b/cherab/core/model/laser/math_functions.pxd @@ -0,0 +1,48 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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 cimport Function3D + +cdef class ConstantAxisymmetricGaussian3D(Function3D): + + cdef: + double _stddev, _normalisation, _kr + + cdef double evaluate(self, double x, double y, double z) except? -1e999 + + +cdef class ConstantBivariateGaussian3D(Function3D): + + cdef: + double _stddev_x, _stddev_y, _kx, _ky, _normalisation + + +cdef class TrivariateGaussian3D(Function3D): + + cdef: + double _mean_z, _stddev_x, _stddev_y, _stddev_z, _kx, _ky + double _kz, _normalisation + + +cdef class GaussianBeamModel(Function3D): + + cdef: + double _waist_z, _stddev_waist, _stddev_waist2, _wavelength, _rayleigh_range + + cdef double evaluate(self, double x, double y, double z) except? -1e999 diff --git a/cherab/core/model/laser/math_functions.pyx b/cherab/core/model/laser/math_functions.pyx new file mode 100644 index 00000000..1ed73a61 --- /dev/null +++ b/cherab/core/model/laser/math_functions.pyx @@ -0,0 +1,304 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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. + + +cimport cython + +from raysect.core.math.function.float cimport Function3D +from raysect.core.math.cython.utility cimport find_index + +from libc.math cimport sqrt, exp, pi + + +cdef class ConstantAxisymmetricGaussian3D(Function3D): + """ + A function with a 2D Gaussian in the x-y plane and equal standard deviations in x and y directions. + + .. math:: + F(x, y, z) = \\frac{1}{2 * \\pi \\sigma^2} exp\\left(-\\frac{x^2 + y^2}{2 * \\sigma^2}\\right) + + The function value has a Gaussian shape in the x-y plane with the standard deviations in + x and y direction being equal. The integral over an x-y plane is equal to 1 + and the mean values in x and y directions are equal to 0. + + :param float stddev: The standard deviation in both the x and y directions. + """ + + def __init__(self, stddev): + + super().__init__() + + self.stddev = stddev + + @property + def stddev(self): + + return self._stddev + + @stddev.setter + def stddev(self, value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._stddev = value + self._kr = -1 / (2 * value ** 2) + self._normalisation = 1 / (2 * pi * value ** 2) + + cdef double evaluate(self, double x, double y, double z) except? -1e999: + cdef: + double r2 + r2 = x ** 2 + y ** 2 + return self._normalisation * exp(r2 * self._kr) + + +cdef class ConstantBivariateGaussian3D(Function3D): + """ + A function with a 2D Gaussian in the x-y plane. + + .. math:: + F(x, y, z) = \\frac{1}{2 * \\pi \\sigma_x \\sigma_y} exp\\left(-\\frac{x^2 + y^2}{2 * \\sigma_x \\sigma_y}\\right) + + The function value has a Gaussian shape in the x-y plane. The integral over an x-y plane is equal to 1 + and the mean values in x and y directions are equal to 0. + The correlation between the standard deviations in x and y directions is equal to 0. + + :param float stddev_x: The standard deviation in the x directions. + :param float stddev_y: The standard deviation in the y directions. + """ + + def __init__(self, stddev_x, stddev_y): + + super().__init__() + self._init_params() + + self.stddev_x = stddev_x + self.stddev_y = stddev_y + + def _init_params(self): + self._stddev_x = 1 + self._stddev_y = 1 + + @property + def stddev_x(self): + return self._stddev_x + + @stddev_x.setter + def stddev_x(self, value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._stddev_x = value + + self._cache_constants() + + @property + def stddev_y(self): + return self._stddev_y + + @stddev_y.setter + def stddev_y(self, value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._stddev_y = value + + self._cache_constants() + + def _cache_constants(self): + self._kx = -1 / (2 * self._stddev_x ** 2) + self._ky = -1 / (2 * self._stddev_y ** 2) + self._normalisation = 1 / (2 * pi * self._stddev_x * self._stddev_y) + + cdef double evaluate(self, double x, double y, double z) except? -1e999: + return self._normalisation * exp(x ** 2 * self._kx + + y ** 2 * self._ky) + + +cdef class TrivariateGaussian3D(Function3D): + """ + A function with a 3D Gaussian shape. + + .. math:: + F(x, y, z) = \\frac{1}{\\sqrt{2 \\pi^3} \\sigma_x \\sigma_y \\sigma_z} exp\\left(-\\frac{x^2}{2 \\sigma_x^2} -\\frac{y^2}{2 \\sigma_y^2} - \\frac{(z - \\mu_z)^2}{2 \\sigma_z^2}\\right) + + The integral over the whole 3D space is equal to 1.The correlation between the standard deviations in x and y directions is equal to 0. The mean value in the + x and y directions are equal to 0. + + :param float mean_z: Mean value in the z direction. + :param float stddev_x: The standard deviation in the x directions. + :param float stddev_y: The standard deviation in the y directions. + :param float stddev_z: The standard deviation in the z directions. + """ + + def __init__(self, mean_z, stddev_x, stddev_y, stddev_z): + + super().__init__() + self._init_params() + + self.stddev_x = stddev_x + self.stddev_y = stddev_y + self.stddev_z = stddev_z + self.mean_z = mean_z + + def _init_params(self): + self._mean_z = 1 + self._stddev_x = 1 + self._stddev_y = 1 + self._stddev_z = 1 + + @property + def stddev_x(self): + return self._stddev_x + + @stddev_x.setter + def stddev_x(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._stddev_x = value + + self._cache_constants() + + @property + def stddev_y(self): + return self._stddev_y + + @stddev_y.setter + def stddev_y(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._stddev_y = value + + self._cache_constants() + + @property + def stddev_z(self): + return self._stddev_z + + @stddev_z.setter + def stddev_z(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._stddev_z = value + + self._cache_constants() + + @property + def mean_z(self): + return self._mean_z + + @mean_z.setter + def mean_z(self, double value): + self._mean_z = value + + def _cache_constants(self): + self._kx = -1 / (2 * self._stddev_x ** 2) + self._ky = -1 / (2 * self._stddev_y ** 2) + self._kz = -1 / (2 * self._stddev_z ** 2) + self._normalisation = 1 / (sqrt((2 * pi) ** 3) * self._stddev_x * self._stddev_y * self._stddev_z) + + cdef double evaluate(self, double x, double y, double z) except? -1e999: + return self._normalisation * exp(x ** 2 * self._kx + + y ** 2 * self._ky + + (z - self._mean_z) ** 2 * self._kz) + + +cdef class GaussianBeamModel(Function3D): + """ + A Gaussian beam function (https://en.wikipedia.org/wiki/Gaussian_beam) + + .. math:: + F(x, y, z) = \\frac{1}{2 \\pi \\sigma^2_z} exp\\left( -\\frac{x^2 + y^2}{2 \\sigma_z(z)^2 }\\right) + + where the standard deviation in the z direction + + .. math:: + \\sigma_z(z) = \\sigma_0 \\sqrt{1 + \\left(\\frac{z - z_0}{z_R}\\right)^2} + + is a function of position and the + + .. math:: + z_R = \\frac{\\pi \\omega_0^2 n}{\\lambda_l} + + is the Rayleigh range. + """ + + def __init__(self, double wavelength, double waist_z, double stddev_waist): + + # preset default values + self._wavelength = 1e3 + self._waist_z = 0 + self._stddev_waist = 1e-3 + + self.wavelength = wavelength + self.waist_z = waist_z + self.stddev_waist = stddev_waist + + @property + def wavelength(self): + return self._wavelength + + @wavelength.setter + def wavelength(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0, but {0} passed.".format(value)) + + self._wavelength = value + self._cache_constants() + + @property + def waist_z(self): + return self._waist_z + + @waist_z.setter + def waist_z(self, double value): + self._waist_z = value + + @property + def stddev_waist(self): + return self._stddev_waist + + @stddev_waist.setter + def stddev_waist(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0, but {0} passed.".format(value)) + + self._stddev_waist = value + self._stddev_waist2 = self._stddev_waist ** 2 + self._cache_constants() + + def _cache_constants(self): + + n = 1 # refractive index of vacuum + self._rayleigh_range = 2 * pi * n * self._stddev_waist2 / self._wavelength / 1e-9 + + @cython.cdivision(True) + cdef double evaluate(self, double x, double y, double z) except? -1e999: + + cdef: + double r2, stddev_z2, z_prime + + # shift to correct gaussiam beam model coords, it works with waist at z=0 + z_prime = z - self._waist_z + + r2 = x ** 2 + y ** 2 + stddev_z2 = self._stddev_waist2 * (1 + ((z_prime) / self._rayleigh_range) ** 2) + + return 1 / (2 * pi * stddev_z2) * exp(r2 / (-2 * stddev_z2)) diff --git a/cherab/core/model/laser/model.pxd b/cherab/core/model/laser/model.pxd new file mode 100644 index 00000000..27e0230f --- /dev/null +++ b/cherab/core/model/laser/model.pxd @@ -0,0 +1,40 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + + +from raysect.optical cimport Vector3D, Point3D +from raysect.optical.spectrum cimport Spectrum + +from cherab.core.laser cimport Laser, LaserModel, LaserSpectrum, LaserProfile + + +cdef class SeldenMatobaThomsonSpectrum(LaserModel): + + cdef: + double _CONST_ALPHA, _RATE_TS, _RECIP_M_PI + + cpdef Spectrum emission(self, Point3D point_plasma, Vector3D observation_plasma, Point3D point_laser, + Vector3D observation_laser, Spectrum spectrum) + + cdef double seldenmatoba_spectral_shape(self, double epsilon, double cos_theta, double alpha) + + cdef Spectrum _add_spectral_contribution(self, double ne, double te, double laser_energy, double angle_pointing, + double angle_polarization, double laser_wavelength, Spectrum spectrum) + + cpdef Spectrum calculate_spectrum(self, double ne, double te, double laser_energy_density, double laser_wavelength, + double observation_angle, double angle_polarization, Spectrum spectrum) \ No newline at end of file diff --git a/cherab/core/model/laser/model.pyx b/cherab/core/model/laser/model.pyx new file mode 100644 index 00000000..dbab8803 --- /dev/null +++ b/cherab/core/model/laser/model.pyx @@ -0,0 +1,220 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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 libc.math cimport exp, sqrt, cos, M_PI, sin +cimport cython + +from raysect.optical cimport Vector3D, Point3D +from raysect.optical.spectrum cimport Spectrum + +from cherab.core cimport Plasma +from cherab.core.laser cimport LaserModel, LaserProfile, LaserSpectrum +from cherab.core.utility.constants cimport DEGREES_TO_RADIANS +from cherab.core.utility.constants cimport SPEED_OF_LIGHT, ELECTRON_CLASSICAL_RADIUS, ELECTRON_REST_MASS, ELEMENTARY_CHARGE + + +cdef class SeldenMatobaThomsonSpectrum(LaserModel): + """ + Thomson Scattering based on Selden-Matoba. + + The class calculates Thomson scattering of the laser to the spectrum. The model of the scattered spectrum used is based on + the semi-empirical model by Selden and the Thomson scattering cross-section is taken from Matoba articles. The spectral contribution + of the scattered laser light c is calculated as a sum of contributions of all laser wavelengths + + .. math:: + c(\lambda) = c r_e^2 n_e cos^2\\theta \\sum_{\\lambda_L} \\frac{E_L(\\lambda_l) S(\\frac{\\lambda}{\\lambda_L} - 1, \\varphi, T_e)}{\\lambda_L}, + + + where :math:`\\lambda` is the spectrum's wavelength, :math:`r_e` is the classical electron radius, :math:`n_e` is the electron delsity, + :math:`\\theta` is the angle between the laser polarisation and scattering vectors, :math:`c` is the vacuum speed of light + :math:`\\lambda_L` is the laser wavelength, :math:`E_L` is the laser energy density, :math:`\\varphi` is the scattering angle and :math:`T_e` is the electron + temperature. The scattering function S is taken from the Matoba article. The multiplication by the speed of light is added to transfer the Thomson scattering + cross section into a reaction rate. + + .. seealso:: + The Prunty article provides a thorough introduction into the phyiscs of Thomson scattering. The articles by Selden and Matoba were used to build + this model. + + :Selden: `Selden, A.C., 1980. Simple analytic form of the relativistic Thomson scattering spectrum. Physics Letters A, 79(5-6), pp.405-406.` + :Matoba: `Matoba, T., et al., 1979. Analytical approximations in the theory of relativistic Thomson scattering for high temperature fusion plasma. + Japanese Journal of Applied Physics, 18(6), p.1127.` + :Prunty: `Prunty, S.L., 2014. A primer on the theory of Thomson scattering for high-temperature fusion plasmas. Physica Scripta, 89(12), p.128001.` + + """ + + def __init__(self, LaserProfile laser_profile=None, LaserSpectrum laser_spectrum=None, Plasma plasma=None): + + super().__init__(laser_profile, laser_spectrum, plasma) + + # Selden, A.C., 1980. Simple analytic form of the relativistic Thomson scattering spectrum. Physics Letters A, 79(5-6), pp.405-406. + self._CONST_ALPHA = ELECTRON_REST_MASS * SPEED_OF_LIGHT ** 2 / (2 * ELEMENTARY_CHARGE) #constant alpha, rewritten for Te in eV + + # from: Prunty, S. L. "A primer on the theory of Thomson scattering for high-temperature fusion plasmas." + # TS cross section equiation ~ 3.28 or + # Matoba, T., et al., 1979. Analytical approximations in the theory of relativistic Thomson scattering for high temperature fusion plasma. + # Japanese Journal of Applied Physics, 18(6), p.1127., TS cross section equiation 18 + # speed of light for correct normalisation of the scattered intensity calculation (from x-section to rate constant) + self._RATE_TS = ELECTRON_CLASSICAL_RADIUS ** 2 * SPEED_OF_LIGHT + + self._RECIP_M_PI = 1 / M_PI + + @cython.cdivision(True) + cdef double seldenmatoba_spectral_shape(self, double epsilon, double const_theta, double alpha): + + cdef: + double c, a, b + + # const_theta is 2 * (1 - cos(theta)) + + c = sqrt(alpha * self._RECIP_M_PI) * (1 - 15. / (16. * alpha) + 345. / (512. * alpha ** 2)) + a = (1 + epsilon) ** 3 * sqrt(const_theta * (1 + epsilon) + epsilon ** 2) + b = sqrt(1 + epsilon ** 2 / (const_theta * (1 + epsilon))) - 1 + + return c / a * exp(-2 * alpha * b) + + @cython.boundscheck(False) + @cython.wraparound(False) + cpdef Spectrum emission(self, Point3D point_plasma, Vector3D observation_plasma, Point3D point_laser, + Vector3D observation_laser, Spectrum spectrum): + cdef: + double angle_scattering, angle_pointing, angle_polarization + double te, ne, laser_energy_density, laser_energy + double plasma_x, plasma_y, plasma_z, laser_x, laser_y, laser_z + double[::1] laser_wavelength_mv, laser_spectrum_power_mv + int bins + Vector3D pointing_vector, polarisation_vector + Py_ssize_t index + + plasma_x = point_plasma.x + plasma_y = point_plasma.y + plasma_z = point_plasma.z + + # get electron parameters for the plasma point + te = self._plasma.get_electron_distribution().effective_temperature(plasma_x, plasma_y, plasma_z) + + #terminate early if electron temperature is 0 + if te <= 0: + return spectrum + + ne = self._plasma.get_electron_distribution().density(plasma_x, plasma_y, plasma_z) + + #terminate early if electron density is 0 + if ne <= 0: + return spectrum + + laser_x = point_laser.x + laser_y = point_laser.y + laser_z = point_laser.z + + #get laser volumetric power + laser_energy_density = self._laser_profile.get_energy_density(laser_x, laser_y, laser_z) + + #terminate early if laser power is 0 + if laser_energy_density == 0: + return spectrum + + pointing_vector = self._laser_profile.get_pointing(laser_x, laser_y, laser_z) + + #angle between observation and pointing vector + angle_pointing = observation_laser.angle(pointing_vector) # angle between observation and pointing vector of laser + + angle_scattering = (180. - angle_pointing) # scattering direction is the opposite to obervation direction + + # angle between polarisation and observation + polarisation_vector = self._laser_profile.get_polarization(laser_x, laser_y, laser_z) + angle_polarization = observation_laser.angle(polarisation_vector) # scattering direction is the opposite to obervation direction + + laser_wavelength_mv = self._laser_spectrum.wavelengths_mv + laser_spectrum_power_mv = self._laser_spectrum.power_mv # power in spectral bins (PSD * delta wavelength) + bins = self._laser_spectrum.get_spectral_bins() + + for index in range(bins): + laser_energy = laser_spectrum_power_mv[index] * laser_energy_density + if laser_energy > 0: + spectrum = self._add_spectral_contribution(ne, te, laser_energy, angle_scattering, + angle_polarization, laser_wavelength_mv[index], spectrum) + + return spectrum + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + cdef Spectrum _add_spectral_contribution(self, double ne, double te, double laser_energy, double angle_scattering, + double angle_polarization, double laser_wavelength, Spectrum spectrum): + + cdef: + int index, nbins + double alpha, epsilon, cos_anglescat, wavelength, min_wavelength, delta_wavelength + double const_theta, recip_laser_wavelength, scattered_power, spectrum_norm + double sin2_angle_pol + + alpha = self._CONST_ALPHA / te + # scattering angle of the photon = pi - observation_angle + cos_anglescat = cos(angle_scattering * DEGREES_TO_RADIANS) + + # pre-calculate constants for Selden-Matoba shape + const_theta = 2 * (1 - cos_anglescat) + + nbins = spectrum.bins + min_wavelength = spectrum.min_wavelength + delta_wavelength = spectrum.delta_wavelength + recip_laser_wavelength = 1 / laser_wavelength + + # dipole radiation has a cos ** 2 characteristic, here angle shifted by 90 deg + sin2_angle_pol = sin(angle_polarization * DEGREES_TO_RADIANS) ** 2 + + #from d_lambda to d_epsilon:d_epsilon = d_lambda / laser_wavelength + scattered_power = ne * self._RATE_TS * laser_energy * recip_laser_wavelength * sin2_angle_pol + for index in range(nbins): + wavelength = min_wavelength + (0.5 + index) * delta_wavelength + epsilon = (wavelength * recip_laser_wavelength) - 1 + spectrum_norm = self.seldenmatoba_spectral_shape(epsilon, const_theta, alpha) + spectrum.samples_mv[index] += spectrum_norm * scattered_power + + return spectrum + + cpdef Spectrum calculate_spectrum(self, double ne, double te, double laser_energy_density, double laser_wavelength, + double observation_angle, double angle_polarization, Spectrum spectrum): + """ + Calculates scattered spectrum for the given parameters. + + The method returns the Thomson scattered spectrum given the plasma parameters, without the need of specifying + plasma or laser. + + :param float ne: Plasma electron density in m**-3 + :param float te: Plasma electron temperature in eV + :param float laser_energy_density: Energy density of the laser light in J * m**-3 + :param float laser_wavelength: The laser light wavelength in nm + :param float observation_angle: The angle of observation is the angle between the observation direction and the direction + of the Poynting vector. + :param float angle_polarization: The angle between the observation direction and the polarisation direction of the laser light. + + :return: Spectrum + """ + # check for nonzero laser power, ne, te, wavelength + if ne <= 0 or te <= 0 or not laser_energy_density > 0: + return spectrum + if laser_wavelength <= 0: + raise ValueError("laser wavelength has to be larger than 0") + + angle_scattering = (180. - observation_angle) # scattering direction is the opposite to obervation direction + + return self._add_spectral_contribution(ne, te, laser_energy_density, angle_scattering, angle_polarization, laser_wavelength, spectrum) + + \ No newline at end of file diff --git a/cherab/core/model/laser/profile.pxd b/cherab/core/model/laser/profile.pxd new file mode 100644 index 00000000..12379335 --- /dev/null +++ b/cherab/core/model/laser/profile.pxd @@ -0,0 +1,57 @@ +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 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 cimport Function3D +from raysect.optical cimport Spectrum, Point3D, Vector3D + +from cherab.core.laser cimport LaserProfile + + +cdef class UniformEnergyDensity(LaserProfile): + + cdef: + double _energy_density, _laser_length, _laser_radius + + +cdef class ConstantAxisymmetricGaussian(LaserProfile): + + cdef: + double _stddev, _pulse_energy, _pulse_length, _laser_length, _laser_radius + Function3D _distribution + + +cdef class ConstantBivariateGaussian(LaserProfile): + + cdef: + double _stddev_x, _stddev_y, _pulse_energy, _pulse_length, _laser_length, _laser_radius + Function3D _distribution + + +cdef class TrivariateGaussian(LaserProfile): + + cdef: + double _stddev_x, _stddev_y, _stddev_z, _mean_z, _pulse_energy, _pulse_length, _laser_length, _laser_radius + Function3D _distribution + + +cdef class GaussianBeamAxisymmetric(LaserProfile): + + cdef: + double _pulse_energy, _pulse_length, _stddev_waist, _waist_z, _laser_wavelength, _laser_length, _laser_radius + Function3D _distribution diff --git a/cherab/core/model/laser/profile.pyx b/cherab/core/model/laser/profile.pyx new file mode 100644 index 00000000..82376980 --- /dev/null +++ b/cherab/core/model/laser/profile.pyx @@ -0,0 +1,765 @@ +from raysect.core.math.function.float import Constant3D +from raysect.core.math.function.vector3d cimport Constant3D as ConstantVector3D +from raysect.primitive import Cylinder +from raysect.optical cimport Spectrum, Vector3D, translate + +from cherab.core.laser cimport Laser, LaserProfile +from cherab.core.model.laser.math_functions cimport ConstantAxisymmetricGaussian3D, ConstantBivariateGaussian3D, TrivariateGaussian3D, GaussianBeamModel + +from cherab.core.utility.constants cimport SPEED_OF_LIGHT + +from libc.math cimport M_PI, sqrt, exp + + +cdef class UniformEnergyDensity(LaserProfile): + """ + LaserProfile with a constant volumetric energy density. + + Returns a laser with a cylindrical shape within which the laser volumentric energy density is constant. + The laser starts at z=0 and extends in the positive z direction. + + .. note: + The methods get_pointing, get_polarization and get_energy_density are not limited to the inside + of the laser cylinder. If called alone for position (x, y, z) outisde the laser cylinder, + they will still return non-zero values. + + In the following example, a laser of length of 2 m (extending from z=0 to z=2 m) with a radius of 3 cm + and volumetric energy density of 5 J*m^-3 and polarisation in the y direction is created: + + .. code-block:: pycon + + >>> from raysect.core import Vector3D + >>> from cherab.core.model.laser import UniformEnergyDensity + + >>> energy = 5 # energy density in J + >>> radius = 3e-2 # laser radius in m + >>> length = 2 # laser length in m + >>> polarisation = Vector3D(0, 1, 0) # polarisation direction + + # create the laser profile + >>> laser_profile = UniformEnergyDensity(energy, radius, length, polarisation) + + :param float energy_density: The volumetric energy density of the laser light. + :param float laser_length: The length of the laser cylinder. + :param float laser_radius: The radius of the laser cylinder. + :param Vector3D polarization: The direction of the laser polarization: + + :ivar float energy_density: The volumetric energy density of the laser light. + :ivar float laser_radius: The radius of the laser cylinder. + :ivar float laser_length: The length of the laser cylinder. + """ + + def __init__(self, double energy_density=1., double laser_length=1., double laser_radius=0.05, Vector3D polarization=Vector3D(0, 1, 0)): + super().__init__() + + self.set_polarization(polarization) + self.set_pointing_function(ConstantVector3D(Vector3D(0, 0, 1))) + self.energy_density = energy_density + + self._laser_radius = 0.05 + self._laser_length = 1. + + self.laser_radius = laser_radius + self.laser_length = laser_length + + def set_polarization(self, Vector3D value): + value = value.normalise() + self.set_polarization_function(ConstantVector3D(value)) + + @property + def laser_length(self): + return self._laser_length + + @laser_length.setter + def laser_length(self, value): + + if value <= 0: + raise ValueError("Laser length has to be larger than 0.") + + self._laser_length = value + self.notifier.notify() + + @property + def laser_radius(self): + return self._laser_radius + + @laser_radius.setter + def laser_radius(self, value): + + if value <= 0: + raise ValueError("Laser radius has to be larger than 0.") + + self._laser_radius = value + self.notifier.notify() + + @property + def energy_density(self): + return self._energy_density + + @energy_density.setter + def energy_density(self, value): + if value <= 0: + raise ValueError("Laser power density has to be larger than 0.") + + self._energy_density = value + funct = Constant3D(value) + self.set_energy_density_function(funct) + + cpdef list generate_geometry(self): + + return generate_segmented_cylinder(self.laser_radius, self.laser_length) + + +cdef class ConstantBivariateGaussian(LaserProfile): + """ + LaserProfile with a Gaussian-shaped volumetric energy density distribution in the xy plane + and constant pulse intensity. + + Returns a laser with a cylindrical shape and the propagation of the laser light in the positive z direction. + + The model imitates a laser beam with a uniform power output within a single pulse. This results + in the distribution of the energy density along the propagation direction of the laser (z-axis) to be also + uniform. The integral value of laser energy Exy in an x-y plane is given by + + .. math:: + E_{xy} = \\frac{E_p}{(c * \\tau)}, + + where Ep is the energy of the laser pulse, tau is the temporal pulse length and c is the speed of light in vacuum. + In an x-y plane, the volumetric energy density follows a bivariate Gaussian with a zero correlation: + + .. math:: + E(x, y) = \\frac{E_{xy}}{2 \\pi \\sigma_x \\sigma_y} exp\\left(-\\frac{x^2 + y^2}{2 \\sigma_x \\sigma_y}\\right). + + The sigma_x and sigma_y are standard deviations in x and y directions, respectively. + + .. note:: + The height of the cylinder, forming the laser beam, is given by the laser_length and is independent from the + temporal length of the laser pulse given by pulse_length. This gives the possibility to independently control + the size of the laser primitive and the value of the volumetric energy density. + + The methods get_pointing, get_polarization and get_energy_density are not limited to the inside + of the laser cylinder. If called for position (x, y, z) outisde the laser cylinder, they can still + return non-zero values. + + + The following example shows how to create a laser with sigma_x= 1 cm and sigma_y=2 cm, which makes the laser + profile in x-y plane to be elliptical. The pulse energy is 5 J and the laser temporal pulse length is 10 ns: + + .. code-block:: pycon + + >>> from raysect.core import Vector3D + >>> from cherab.core.model.laser import ConstantBivariateGaussian + + >>> radius = 3e-2 # laser radius in m + >>> length = 2 # laser length in m + >>> polarisation = Vector3D(0, 1, 0) # polarisation direction + >>> pulse_energy = 5 # energy in a laser pulse in J + >>> pulse_length = 1e-8 # pulse length in s + >>> width_x = 1e-2 # standard deviation in x direction in m + >>> width_y = 2e-2 # standard deviation in y direction in m + + # create the laser profile + >>> laser_profile = ConstantBivariateGaussian(pulse_energy, pulse_length, radius, length, width_x, width_y, polarisation) + + :param float pulse_energy: The energy of the laser in Joules delivered in a single laser pulse. + :param float pulse_length: The temporal length of the laser pulse in seconds. + :param float laser_length: The length of the laser cylinder. + :param float laser_radius: The radius of the laser cylinder. + :param float stddev_x: The standard deviation of the bivariate Gaussian distribution of the volumetric energy + density distribution of the laser light in the x axis in meters. + :param float stddev_y: The standard deviation of the bivariate Gaussian distribution of the volumetric energy + density distribution of the laser light in the y axis in meters. + :param Vector3D polarization: The direction of the laser polarization: + + :ivar float pulse_energy: The energy of the laser in Joules delivered in a single laser pulse. + :ivar float pulse_length: The temporal length of the laser pulse in seconds. + :ivar float laser_radius: The radius of the laser cylinder. + :ivar float laser_length: The length of the laser cylinder. + :ivar float stddev_x: The standard deviation of the bivariate Gaussian distribution of the volumetric energy + density distribution of the laser light in the x axis in meters. + :ivar float stddev_y: The standard deviation of the bivariate Gaussian distribution of the volumetric energy + density distribution of the laser light in the y axis in meters. + """ + def __init__(self, double pulse_energy=1, double pulse_length=1, double laser_radius=0.05, double laser_length=1., + double stddev_x=0.01, double stddev_y=0.01, Vector3D polarization=Vector3D(0, 1, 0)): + + super().__init__() + # set initial values + self._pulse_energy = 1 + self._pulse_length = 1 + self._stddev_x = 0.1 + self._stddev_y = 0.1 + + self._laser_radius = 0.05 + self._laser_length = 1 + + self.laser_radius = laser_radius + self.laser_length = laser_length + + self.set_polarization(polarization) + self.set_pointing_function(ConstantVector3D(Vector3D(0, 0, 1))) + + self.stddev_x = stddev_x + self.stddev_y = stddev_y + + self.pulse_energy = pulse_energy + self.pulse_length = pulse_length + + # set laser constants + self.set_polarization(polarization) + + def set_polarization(self, Vector3D value): + value = value.normalise() + self.set_polarization_function(ConstantVector3D(value)) + + @property + def laser_length(self): + return self._laser_length + + @laser_length.setter + def laser_length(self, value): + + if value <= 0: + raise ValueError("Laser length has to be larger than 0.") + + self._laser_length = value + self.notifier.notify() + + @property + def laser_radius(self): + return self._laser_radius + + @laser_radius.setter + def laser_radius(self, value): + + if value <= 0: + raise ValueError("Laser radius has to be larger than 0.") + + self._laser_radius = value + self.notifier.notify() + + @property + def pulse_energy(self): + return self._pulse_energy + + @pulse_energy.setter + def pulse_energy(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._pulse_energy = value + self.notifier.notify() + + @property + def pulse_length(self): + return self._pulse_length + + @pulse_length.setter + def pulse_length(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._pulse_length = value + self._function_changed() + + @property + def stddev_x(self): + return self._stddev_x + + @stddev_x.setter + def stddev_x(self, value): + if value <= 0: + raise ValueError("Standard deviation of the laser power has to be larger than 0.") + + self._stddev_x = value + self._function_changed() + + @property + def stddev_y(self): + return self._stddev_y + + @stddev_y.setter + def stddev_y(self, value): + if value <= 0: + raise ValueError("Standard deviation of the laser power has to be larger than 0.") + + self._stddev_y = value + self._function_changed() + + def _function_changed(self): + """ + Energy density should be returned in units [J/m ** 3]. Energy shape in xy + plane is defined by normal distribution (integral over xy plane for + constant z is 1). The units of such distribution are [m ** -2]. + In the z axis direction (direction of laser propagation), + the laser_energy is spread along the z axis using the velocity + of light SPEED_OF_LIGHT and the temporal duration of the pulse: + length = SPEED_OF_LIGTH * pulse_length. Combining the normal distribution with the normalisation + pulse_energy / length gives the units [J / m ** 3]. + """ + self._distribution = ConstantBivariateGaussian3D(self._stddev_x, self._stddev_y) + + length = SPEED_OF_LIGHT * self._pulse_length # convert from temporal to spatial length of pulse + normalisation = self._pulse_energy / length # normalisation to have correct spatial energy density [J / m**3] + + function = normalisation * self._distribution + self.set_energy_density_function(function) + + cpdef list generate_geometry(self): + + return generate_segmented_cylinder(self.laser_radius, self.laser_length) + + +cdef class TrivariateGaussian(LaserProfile): + """ + LaserProfile with a trivariate Gaussian-shaped volumetric energy density. + + Returns a laser with a cylindrical shape and the propagation of the laser light in the positive z direction. + This model imitates a laser beam with a Gaussian distribution of power output within a single pulse frozen in time: + + .. math:: + E(x, y, z) = \\frac{E_p}{\\sqrt{2 \\pi^3} \\sigma_x \\sigma_y \\sigma_z} exp\\left(-\\frac{x^2}{2 \\sigma_x^2} -\\frac{y^2}{2 \\sigma_y^2} -\\frac{(z - \\mu_z)^2}{2 \\sigma_z^2}\\right). + + + The sigma_x and sigma_y are standard deviations in x and y directions, respectively, and E_p is the energy deliverd by laser in a + single laser pulse. The mu_z is the mean of the distribution in the z direction and controls th position of the laser pulse along the z direction. + The standard deviation in z direction sigma_z is calculated from the pulse length tau_p, which is the + standard deviation of the Gaussian distributed ouput power of the laser within a single pulse: + + .. math:: + \\sigma_z = \\tau_p c. + + The c stands for the speed of light in vacuum. + + .. note:: + The height of the cylinder, forming the laser beam, is given by the laser_length and is independent from the + temporal length of the laser pulse given by pulse_length. This gives the possibility to independently control + the size of the laser primitive and the value of the volumetric energy density. + + The methods get_pointing, get_polarization and get_energy_density are not limited to the inside + of the laser cylinder. If called alone for position (x, y, z) outisde the laser cylinder, they can still + return non-zero values. + + + The following example shows how to create a laser with sigma_x = 1 cm and sigma_y = 2 cm, which makes the laser + profile in an x-y plane to be elliptical. The pulse energy is 5 J and the laser temporal pulse length is 10 ns. + The position of the laser pulse maximum mean_z is set to 0.5: + + .. code-block:: pycon + + >>> from raysect.core import Vector3D + >>> from cherab.core.model.laser import ConstantBivariateGaussian + + >>> radius = 3e-2 # laser radius in m + >>> length = 2 # laser length in m + >>> polarisation = Vector3D(0, 1, 0) # polarisation direction + >>> pulse_energy = 5 # energy in a laser pulse in J + >>> pulse_length = 1e-8 # pulse length in s + >>> pulse_z = 0.5 # position of the pulse mean + >>> width_x = 1e-2 # standard deviation in x direction in m + >>> width_y = 2e-2 # standard deviation in y direction in m + + # create the laser profile + >>> laser_profile = ConstantBivariateGaussian(pulse_energy, pulse_length, pulse_z, radius, length, width_x, width_y, polarisation) + + + :param float pulse_energy: The energy of the laser in Joules delivered in a single laser pulse. + :param float pulse_length: The standard deviation of the laser pulse length in the temporal domain. + :param float mean_z: Position of the mean value of the laser pulse in the z direction. Can be used to control the + position of the laser pulse along the laser propagation. + :param float laser_length: The length of the laser cylinder. + :param float laser_radius: The radius of the laser cylinder. + :param float stddev_x: The standard deviation of the bivariate Gaussian distribution of the volumetric energy + density distribution of the laser light in the x axis in meters. + :param float stddev_y: The standard deviation of the bivariate Gaussian distribution of the volumetric energy + density distribution of the laser light in the y axis in meters. + :param Vector3D polarization: The direction of the laser polarization. + + :ivar float pulse_energy: The energy of the laser in Joules delivered in a single laser pulse. + :ivar float pulse_length: The standard deviation of the laser pulse length in the temporal domain. + :ivar float mean_z: Position of the mean value of the laser pulse in the z direction. + Can be used to control the position of the laser pulse along the laser propagation. + :ivar float laser_radius: The radius of the laser cylinder. + :ivar float laser_length: The length of the laser cylinder. + :ivar float stddev_x: The standard deviation of the bivariate Gaussian distribution of the volumetric energy + density distribution of the laser light in the x axis in meters. + :ivar float stddev_y: The standard deviation of the bivariate Gaussian distribution of the volumetric energy + density distribution of the laser light in the y axis in meters. + """ + def __init__(self, double pulse_energy=1, double pulse_length=1, double mean_z=0, + double laser_length=1., double laser_radius=0.05, + double stddev_x=0.01, double stddev_y=0.01, + Vector3D polarization=Vector3D(0, 1, 0)): + + super().__init__() + # set initial values + self._pulse_energy = 1 + self._pulse_length = 1 + self._stddev_x = 0.1 + self._stddev_y = 0.1 + self._stddev_z = 1 + self._laser_radius = 0.05 + self._laser_length = 1 + self._mean_z = mean_z + + + self.laser_radius = laser_radius + self.laser_length = laser_length + self.stddev_x = stddev_x + self.stddev_y = stddev_y + self.mean_z = mean_z + + self.pulse_energy = pulse_energy + self.pulse_length = pulse_length + + self.set_polarization(polarization) + self.set_pointing_function(ConstantVector3D(Vector3D(0, 0, 1))) + + def set_polarization(self, Vector3D value): + value = value.normalise() + self.set_polarization_function(ConstantVector3D(value)) + + @property + def laser_length(self): + return self._laser_length + + @laser_length.setter + def laser_length(self, value): + + if value <= 0: + raise ValueError("Laser length has to be larger than 0.") + + self._laser_length = value + self.notifier.notify() + + @property + def laser_radius(self): + return self._laser_radius + + @laser_radius.setter + def laser_radius(self, value): + + if value <= 0: + raise ValueError("Laser radius has to be larger than 0.") + + self._laser_radius = value + self.notifier.notify() + + @property + def pulse_energy(self): + return self._pulse_energy + + @pulse_energy.setter + def pulse_energy(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._pulse_energy = value + self._function_changed() + + @property + def pulse_length(self): + return self._pulse_length + + @pulse_length.setter + def pulse_length(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._pulse_length = value + self._stddev_z = self._pulse_length * SPEED_OF_LIGHT + self._function_changed() + + @property + def stddev_x(self): + return self._stddev_x + + @stddev_x.setter + def stddev_x(self, value): + if value <= 0: + raise ValueError("Standard deviation of the laser power has to be larger than 0.") + + self._stddev_x = value + self._function_changed() + + @property + def stddev_y(self): + return self._stddev_y + + @stddev_y.setter + def stddev_y(self, value): + if value <= 0: + raise ValueError("Standard deviation of the laser power has to be larger than 0.") + + self._stddev_y = value + self._function_changed() + + @property + def mean_z(self): + return self._mean_z + + @mean_z.setter + def mean_z(self, double value): + self._mean_z = value + self._function_changed() + + def _function_changed(self): + """ + Energy density should be returned in units [J/m ** 3]. The integral value of the _distribution + is 1, thus multiplying _distribution by _pulse_energy gives correct values. + """ + + self._distribution = TrivariateGaussian3D(self._mean_z, self._stddev_x, self._stddev_y, + self._stddev_z) + + normalisation = self._pulse_energy + + function = normalisation * self._distribution + self.set_energy_density_function(function) + + cpdef list generate_geometry(self): + + return generate_segmented_cylinder(self.laser_radius, self.laser_length) + +cdef class GaussianBeamAxisymmetric(LaserProfile): + """ + LaserProfile with volumetric energy density following the Gaussian beam model. + + Returns a laser with a cylindrical shape and the propagation of the laser light in the positive z direction. This model implements + the axisymmetrical Gaussian beam model. It immitates a focused axis symmetrical laser beam with a uniform power ouput in a laser pulse. + The volumetric energy density is given by + + .. math:: + E(x, y, z) = \\frac{E_{xy}}{2 \\pi \\sigma^2(z)} exp\\left( -\\frac{x^2 + y^2}{2 \\sigma^2(z) }\\right) \\\\ + + where the sigma is the standard deviation of the Gaussian shape in the xy plane and is given by + + .. math:: + sigma(z) = \\sigma_0 \\sqrt{1 + \\left(\\frac{z - z_0}{z_R}\\right)^2}. + + The z_0 is the position of the beam focus and z_R is the Rayleigh length + + .. math:: + z_R = \\frac{\\pi \\omega_0^2 n}{\\lambda_l} + + where the omega_0 is the standard deviation in the xy plane in the focal point (beam waist) and lambda_l is the central wavelength of + the laser. The E_xy stand for the laser energy in an xy plane and is calculated as: + + .. math:: + E_{xy} = \\frac{E_p}{(c * \\tau)}, + + where the E_p is the energy in a single laser pulse and tau is the temporal pulse length. + + .. note:: + For more information about the Gaussian beam model see https://en.wikipedia.org/wiki/Gaussian_beam + + The methods get_pointing, get_polarization and get_energy_density are not limited to the inside + of the laser cylinder. If called alone for position (x, y, z) outisde the laser cylinder, they can still + return non-zero values. + + The following example shows how to create a laser with pulse energy 5J, pulse length 10 ns and with the laser cylinder primitive + being 2m long with 5 cm in diameter. The the standard deviation of the beam in the focal point (waist) is 5mm and the position of the + waist is z=50 cm. The laser wavelength is 1060 nm. + + .. code-block:: pycon + + >>> from raysect.core import Vector3D + >>> from cherab.core.model.laser import GaussianBeamAxisymmetric + + >>> radius = 5e-2 # laser radius in m + >>> length = 2 # laser length in m + >>> polarisation = Vector3D(0, 1, 0) # polarisation direction + >>> pulse_energy = 5 # energy in a laser pulse in J + >>> pulse_length = 1e-8 # pulse length in s + >>> waist_width = 5e-3 # standard deviation in the waist + >>> waist_z = 0.5 # position of the pulse mean + >>> width_x = 1e-2 # standard deviation in x direction in m + >>> width_y = 2e-2 # standard deviation in y direction in m + >>> laser_wlen = 1060 # laser wavelength in nm + + # create the laser profile + >>> laser_profile = GaussianBeamAxisymmetric(pulse_energy, pulse_length, length, radius, waist_z, waist_width, laser_wlen) + + :param float pulse_energy: The energy of the laser in Joules delivered in a single laser pulse. + :param float pulse_length: The temporal length of the laser pulse in seconds. + :param float laser_length: The length of the laser cylinder in meters. + :param float laser_radius: The radius of the laser cylinder in meters. + :param float waist_z: Position of the laser waist along the z axis in m. + :param float stddev_waist: The standard deviation of the laser width in the focal point (waist) in m. + :param float laser_wavelength: The central wavelength of the laser light in nanometers. + :param Vector3D polarization: The direction of the laser polarization. + + :ivar float pulse_energy: The energy of the laser in Joules delivered in a single laser pulse. + :ivar float pulse_length: The temporal length of the laser pulse in seconds. + :ivar float laser_length: The length of the laser cylinder in meters. + :ivar float laser_radius: The radius of the laser cylinder in meters. + :ivar float waist_z: Position of the laser waist along the z axis in m. + :ivar float stddev_waist: The standard deviation of the laser width in the focal point (waist) in m. + :ivar float laser_wavelength: The central wavelength of the laser light in nanometers. + :ivar Vector3D polarization: The direction of the laser polarization. + """ + + def __init__(self, double pulse_energy=1, double pulse_length=1, + double laser_length=1., double laser_radius=0.05, + double waist_z=0, double stddev_waist=0.01, + double laser_wavelength=1e3, Vector3D polarization=Vector3D(0, 1, 0)): + + super().__init__() + # set initial values + self._pulse_energy = 1 + self._pulse_length = 1 + self._stddev_waist = 0.1 + self._waist_z = waist_z + self._laser_wavelength = 1e3 + self._laser_radius = 0.05 + self._laser_length = 1 + + self.laser_length = laser_length + self.laser_radius = laser_radius + + self.set_polarization(polarization) + self.set_pointing_function(ConstantVector3D(Vector3D(0, 0, 1))) + + self.stddev_waist = stddev_waist + self.waist_z = waist_z + + self.pulse_energy = pulse_energy + self.pulse_length = pulse_length + self.laser_wavelength = laser_wavelength + + def set_polarization(self, Vector3D value): + value = value.normalise() + self.set_polarization_function(ConstantVector3D(value)) + + @property + def laser_length(self): + return self._laser_length + + @laser_length.setter + def laser_length(self, value): + + if value <= 0: + raise ValueError("Laser length has to be larger than 0.") + + self._laser_length = value + self.notifier.notify() + + @property + def laser_radius(self): + return self._laser_radius + + @laser_radius.setter + def laser_radius(self, value): + + if value <= 0: + raise ValueError("Laser radius has to be larger than 0.") + + self._laser_radius = value + self.notifier.notify() + + @property + def pulse_energy(self): + return self._pulse_energy + + @pulse_energy.setter + def pulse_energy(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._pulse_energy = value + self._function_changed() + + @property + def pulse_length(self): + return self._pulse_length + + @pulse_length.setter + def pulse_length(self, double value): + if value <= 0: + raise ValueError("Value has to be larger than 0.") + + self._pulse_length = value + self._function_changed() + + @property + def waist_z(self): + return self._waist_z + + @waist_z.setter + def waist_z(self, double value): + self._waist_z = value + self._function_changed() + + @property + def stddev_waist(self): + return self._stddev_waist + + @stddev_waist.setter + def stddev_waist(self, double value): + self._stddev_waist = value + self._function_changed() + + @property + def laser_wavelength(self): + return self._laser_wavelength + + @laser_wavelength.setter + def laser_wavelength(self, double value): + self._laser_wavelength = value + self._function_changed() + + def _function_changed(self): + """ + Energy density should be returned in units [J/m ** 3]. Energy shape in xy + plane is defined by normal distribution (integral over xy plane for + constant z is 1). The units of such distribution are [m ** -2]. + In the z axis direction (direction of laser propagation), + the laser_energy is spread along the z axis using the velocity + of light SPEED_OF_LIGHT and the temporal duration of the pulse: + length = SPEED_OF_LIGTH * pulse_length. Combining the normal distribution with the normalisation + pulse_energy / length gives the units [J / m ** 3]. + """ + + self._distribution = GaussianBeamModel(self._laser_wavelength, self._waist_z, self._stddev_waist) + # calculate volumetric power dentiy + length = SPEED_OF_LIGHT * self._pulse_length # convert from temporal to spatial length of pulse + normalisation = self._pulse_energy / length # normalisation to have correct spatial energy density [J / m**3] + + function = normalisation * self._distribution + self.set_energy_density_function(function) + + cpdef list generate_geometry(self): + + return generate_segmented_cylinder(self.laser_radius, self.laser_length) + + +def generate_segmented_cylinder(radius, length): + """ + Generates a segmented cylindrical laser geometry + + Approximates a long cylinder with a cylindrical segments to optimize + targetted and importance sampling. The height of a cylinder segments is roughly + 2 * cylinder radius. + + :return: List of cylinders + """ + + n_segments = int(length // (2 * radius)) # number of segments + geometry = [] + + #length of segment is either length / n_segments if length > radius or length if length < radius + if n_segments > 1: + segment_length = length / n_segments + for i in range(n_segments): + segment = Cylinder(name="Laser segment {0:d}".format(i), radius=radius, height=segment_length, + transform=translate(0, 0, i * segment_length)) + + geometry.append(segment) + elif 0 <= n_segments < 2: + segment = Cylinder(name="Laser segment {0:d}".format(0), radius=radius, height=length) + + geometry.append(segment) + else: + raise ValueError("Incorrect number of segments calculated.") + + return geometry \ No newline at end of file diff --git a/cherab/core/model/laser/tests/__init__.py b/cherab/core/model/laser/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/cherab/core/model/laser/tests/test_laserspectrum.py b/cherab/core/model/laser/tests/test_laserspectrum.py new file mode 100644 index 00000000..1f72e666 --- /dev/null +++ b/cherab/core/model/laser/tests/test_laserspectrum.py @@ -0,0 +1,46 @@ +import unittest +from scipy.integrate import nquad + +from cherab.core.model.laser.laserspectrum import GaussianSpectrum, ConstantSpectrum + +class TestLaserSpectrum(unittest.TestCase): + + def test_constantspectrum(self): + """ + Laser spectrum should be normalized, i.e. integral from minuns inf. to inf. should be one. + :return: + """ + min_wavelength = 1039.9 + max_wavelength = 1040.1 + bins = 10 + + spectrum = ConstantSpectrum(min_wavelength, max_wavelength, bins) + + # check if the power_spectral density is normalized + + integral = spectrum.power_spectral_density.sum() * spectrum.delta_wavelength + self.assertTrue(integral == 1, msg="Power spectral density is not normalised.") + + def test_gaussian_spectrum(self): + """ + Laser spectrum should be normalized, i.e. integral from minuns inf. to inf. should be one. + :return: + """ + min_wavelength = 1035 + max_wavelength = 1045 + bins = 100 + mean = 1040 + stddev = 0.5 + + spectrum = GaussianSpectrum(min_wavelength, max_wavelength, bins, mean, stddev) + integral = nquad(spectrum, [(min_wavelength, max_wavelength)])[0] + + # check if the power_spectral density is normalized + self.assertAlmostEqual(integral, 1., 8, msg="Power spectral density function is not normalised.") + + psd = spectrum.power_spectral_density + integral = 0 + for index in range(0, spectrum.bins - 1): + integral += (psd[index] + psd[index + 1]) / 2 * spectrum.delta_wavelength + + self.assertAlmostEqual(integral, 1, 8, msg="Power spectral density is not normalised.") diff --git a/cherab/core/model/laser/tests/test_model.py b/cherab/core/model/laser/tests/test_model.py new file mode 100644 index 00000000..992d5d52 --- /dev/null +++ b/cherab/core/model/laser/tests/test_model.py @@ -0,0 +1,104 @@ +import unittest +import numpy as np +from math import cos, sin, sqrt, radians, exp +from scipy.constants import pi, c, e, m_e, epsilon_0 + +from raysect.optical import World, Point3D, Vector3D, translate, Ray + +from cherab.core.laser import Laser +from cherab.core.model.laser import ConstantSpectrum, SeldenMatobaThomsonSpectrum, UniformEnergyDensity + +from cherab.tools.plasmas.slab import build_constant_slab_plasma + +class TestScatteringModel(unittest.TestCase): + laser_wavelength = 1040 + wavelength = np.linspace(600, 1200, 800) + scatangle = 45 + + def test_selden_matoba_scattered_spectrum(self): + + # calculate TS cross section and constants + r_e = e ** 2 / (4 * pi * epsilon_0 * m_e * c ** 2) # classical electron radius + + # re ** 2 is the cross section, c transforms xsection into rate constant + scat_const = r_e ** 2 * c + + ray_origin = Point3D(0, 0, 0) + + # angle of scattering + observation_angle = [45, 90, 135] + for obsangle in observation_angle: + + # pointing vector is in +z direction, angle of observation is 180 - obsangle + z = cos((obsangle) / 180 * pi) + x = sin((obsangle) / 180 * pi) + ray_direction = Vector3D(x, 0, z).normalise() + + # ray spectrum properties + min_wavelength = 600 + max_wavelength = 1200 + bins = 800 + + # plasma properties + e_density = 8e19 + t_e = [100, 1e3, 1e4] + + for vte in t_e: + # setup scene + world = World() + plasma = build_constant_slab_plasma(length=1, width=1, height=1, electron_density=e_density, + electron_temperature=vte, plasma_species=[], parent=world) + + # setup laser + laser_spectrum = ConstantSpectrum(1059, 1061, 1) + laser_profile = UniformEnergyDensity(laser_length=1, laser_radius=0.015) + scattering_model = SeldenMatobaThomsonSpectrum() + laser = Laser() + laser.parent = world + laser.transform = translate(0.05, 0, -0.5) + laser.laser_profile = laser_profile + laser.laser_spectrum = laser_spectrum + laser.plasma = plasma + laser.models = [scattering_model] + + # trace a single ray through the laser + ray = Ray(origin=ray_origin, direction=ray_direction, min_wavelength=min_wavelength, + max_wavelength=max_wavelength, bins=bins) + traced_spectrum = ray.trace(world) + + + # calculate spectrum ray-tracing should deliver + dl = 2 * laser_profile.laser_radius / sin((obsangle) / 180 * pi) # ray-laser cross section length + intensity_test = np.zeros_like(traced_spectrum.wavelengths) + + for vl in laser.laser_spectrum.wavelengths: + intensity_const = scat_const * e_density / vl * dl + for iwvl, vwvl in enumerate(traced_spectrum.wavelengths): + intensity_test[iwvl] += _selden_matoba_shape(vwvl, vte, obsangle, vl) * intensity_const + + for index, (vtest, vray) in enumerate(zip(intensity_test, traced_spectrum.samples)): + # skip the test for too low values of power spectral density, max is approx. 3e-5 + if vray > 1e-30: + rel_error = np.abs((vtest - vray) / vray) + self.assertLess(rel_error, 1e-7, + msg="Traced and test spectrum value do not match: " + "scattering angle = {0} deg, Te = {1} eV, wavelength = {2} nm." + .format(180 - obsangle, vte, traced_spectrum.wavelengths[index])) + + +def _selden_matoba_shape(wavelength, te, obsangle, laser_wavelength): + """ + Returns Selden-Matoba Spectral shape + """ + epsilon = wavelength / laser_wavelength - 1 + + alpha = m_e * c ** 2 / (2 * e * te) + + scat_angle = 180 - obsangle + cos_scat = cos(radians(scat_angle)) + + const_c = sqrt(alpha / pi) * (1 - 15 / 16 / alpha + 345 / 512 / alpha ** 2) + const_a = (1 + epsilon) ** 3 * sqrt(2 * (1 - cos_scat) * (1 + epsilon) + epsilon ** 2) + const_b = sqrt(1 + epsilon ** 2 / (2 * (1 - cos_scat) * (1 + epsilon))) - 1 + + return const_c / const_a * exp(-2 * alpha * const_b) diff --git a/cherab/core/model/laser/tests/test_profiles.py b/cherab/core/model/laser/tests/test_profiles.py new file mode 100644 index 00000000..6f2c9fc3 --- /dev/null +++ b/cherab/core/model/laser/tests/test_profiles.py @@ -0,0 +1,198 @@ +import unittest +import numpy as np +from math import exp, sqrt + +from raysect.core import Vector3D + +from cherab.core.model.laser.profile import UniformEnergyDensity, ConstantBivariateGaussian +from cherab.core.model.laser.profile import TrivariateGaussian, GaussianBeamAxisymmetric, generate_segmented_cylinder + +from scipy.integrate import nquad +from scipy.constants import c, pi + +class TestSegmentedCylinder(unittest.TestCase): + + def test_number_of_primitives(self): + + # for r > l there should be 1 cylinder segment only + primitives = generate_segmented_cylinder(radius=1, length=0.5) + self.assertEqual(len(primitives), 1, msg="Wrong nuber of laser segments, expected 1.") + + # for r < l tehre should be length // (2 * radius) segments + primitives = generate_segmented_cylinder(radius=0.5, length=10) + self.assertEqual(len(primitives), 10, msg="Wrong nuber of laser segments, expected 20.") + + +class TestLaserProfile(unittest.TestCase): + + def test_uniform_energy_density(self): + polarisation = Vector3D(1, 3, 8).normalise() + energy_density = 2 + model = UniformEnergyDensity(energy_density=energy_density, polarization=polarisation) + + # test polarisation + pol_model = model.get_polarization(1, 1, 1) + self.assertEqual(pol_model.x, polarisation.x, + msg="Model polarization x vector component does not agreee with input value.") + self.assertEqual(pol_model.y, polarisation.y, + msg="Model polarization y vector component does not agreee with input value.") + self.assertEqual(pol_model.z, polarisation.z, + msg="Model polarization z vector component does not agreee with input value.") + + # test power + self.assertEqual(model.get_energy_density(3, 4, 1), energy_density, + msg="Model power density distribution does not agree with input.") + + def test_bivariate_gaussian(self): + + pulse_energy = 2 + pulse_length = 1e-8 + stddev_x = 0.03 + stddev_y = 0.06 + polarisation = Vector3D(2, 3, 4).normalise() + model = ConstantBivariateGaussian(pulse_energy=pulse_energy, pulse_length=pulse_length, stddev_x=stddev_x, + stddev_y=stddev_y, polarization=polarisation) + + # test polarisation + pol_model = model.get_polarization(1, 1, 1) + self.assertEqual(pol_model.x, polarisation.x, + msg="Model polarization x vector component does not agreee with input value.") + self.assertEqual(pol_model.y, polarisation.y, + msg="Model polarization y vector component does not agreee with input value.") + self.assertEqual(pol_model.z, polarisation.z, + msg="Model polarization z vector component does not agreee with input value.") + + # Integrate over laser volume to check energy + xlim = [-5 * stddev_x, 5 * stddev_x] + ylim = [-5 * stddev_y, 5 * stddev_y] + zlim = [0, pulse_length * c] + energy_integrated = nquad(model.get_energy_density, [xlim, ylim, zlim])[0] + self.assertTrue(np.isclose(energy_integrated / pulse_energy, 1, 1e-3), + msg="Integrated laser energy of the model does not give results close to input energy") + + # Check laser power density profile + x = np.linspace(-3 * stddev_x, 3 * stddev_x, 30) + y = np.linspace(-3 * stddev_y, 3 * stddev_y, 30) + for vx in x: + for vy in y: + tmp = _constant_bivariate_gaussian2d(vx, vy, pulse_energy, pulse_length, stddev_x, stddev_y) + tmp2 = model.get_energy_density(vx, vy, 0) + self.assertTrue(np.isclose(tmp, tmp2, 1e-9), + msg="Model power density distribution for ({},{},{}) does not agree with input.".format( + vx, vy, 0)) + + def test_trivariate_gaussian(self): + + pulse_energy = 2 + pulse_length = 1e-8 + mean_z = 0 + stddev_x = 0.03 + stddev_y = 0.06 + polarisation = Vector3D(2, 3, 4).normalise() + model = TrivariateGaussian(pulse_energy=pulse_energy, pulse_length=pulse_length, mean_z=mean_z, + stddev_x=stddev_x, stddev_y=stddev_y, polarization=polarisation) + + # test polarisation + pol_model = model.get_polarization(1, 1, 1) + self.assertEqual(pol_model.x, polarisation.x, + msg="Model polarization x vector component does not agreee with input value.") + self.assertEqual(pol_model.y, polarisation.y, + msg="Model polarization y vector component does not agreee with input value.") + self.assertEqual(pol_model.z, polarisation.z, + msg="Model polarization z vector component does not agreee with input value.") + + # Integrate over laser volume to check energy + xlim = [-5 * stddev_x, 5 * stddev_x] + ylim = [-5 * stddev_y, 5 * stddev_y] + zlim = [mean_z - 5 * pulse_length * c, mean_z + 5 * pulse_length * c] + energy_integrated = nquad(model.get_energy_density, [xlim, ylim, zlim])[0] + + self.assertTrue(np.isclose(energy_integrated / pulse_energy, 1, 1e-3), + msg="Integrated laser energy of the model does not give results close to input energy") + + # Check laser power density profile + x = np.linspace(-3 * stddev_x, 3 * stddev_x, 30) + y = np.linspace(-3 * stddev_y, 3 * stddev_y, 30) + z = np.linspace(-3 * pulse_length * c, 3 * pulse_length * c, 30) + + for vx in x: + for vy in y: + for vz in z: + tmp = _constant_trivariate_gaussian3d(vx, vy, vz, pulse_energy, pulse_length, mean_z, stddev_x, + stddev_y) + tmp2 = model.get_energy_density(vx, vy, vz) + self.assertTrue(np.isclose(tmp, tmp2, 1e-9), + msg="Model power density distribution for ({},{},{})" + " does not agree with input.".format(vx, vy, vz)) + + def test_gaussianbeam(self): + + pulse_energy = 2 + pulse_length = 1e-9 + waist_z = 0 + stddev_waist = 0.003 + laser_wavelength = 1040 + polarisation = Vector3D(2, 3, 4).normalise() + + model = GaussianBeamAxisymmetric(pulse_energy=pulse_energy, pulse_length=pulse_length, waist_z=waist_z, + stddev_waist=stddev_waist, laser_wavelength=laser_wavelength, + polarization=polarisation) + + # test polarisation + pol_model = model.get_polarization(1, 1, 1) + self.assertEqual(pol_model.x, polarisation.x, + msg="Model polarization x vector component does not agreee with input value.") + self.assertEqual(pol_model.y, polarisation.y, + msg="Model polarization y vector component does not agreee with input value.") + self.assertEqual(pol_model.z, polarisation.z, + msg="Model polarization z vector component does not agreee with input value.") + + # Integrate over laser volume to check energy + xlim = [-20 * stddev_waist, 20 * stddev_waist] + zlim = [-1 * pulse_length / 2 * c, pulse_length / 2 * c] + energy_integrated = nquad(model.get_energy_density, [xlim, xlim, zlim])[0] + + self.assertTrue(np.isclose(energy_integrated / pulse_energy, 1, 1e-3), + msg="Integrated laser energy of the model does not give results close to input energy") + + # Check laser power density profile + x = np.linspace(-3 * stddev_waist, 3 * stddev_waist, 30) + y = np.linspace(-3 * stddev_waist, 3 * stddev_waist, 30) + z = np.linspace(-3 * pulse_length * c, 3 * pulse_length * c, 30) + + for vx in x: + for vy in y: + for vz in z: + tmp = _gaussian_beam_model(vx, vy, vz, pulse_energy, pulse_length, waist_z, stddev_waist, + laser_wavelength * 1e-9) + tmp2 = model.get_energy_density(vx, vy, vz) + self.assertTrue(np.isclose(tmp, tmp2, 1e-9), + msg="Model power density distribution for ({},{},{}) " + "does not agree with input.".format(vx, vy, vz)) + + +def _gaussian_beam_model(x, y, z, pulse_energy, pulse_length, waist_z, stddev_waist, wavelength): + laser_power_axis = pulse_energy / (pulse_length * c) + + n = 1 # refractive index + rayleigh_distance = 2 * pi * stddev_waist ** 2 * n / wavelength + + z_prime = z - waist_z + + stddev_z2 = stddev_waist ** 2 * (1 + (z_prime / rayleigh_distance) ** 2) + + r2 = x ** 2 + y ** 2 + + return laser_power_axis / (2 * pi * stddev_z2) * exp(r2 / (-2 * stddev_z2)) + + +def _constant_trivariate_gaussian3d(x, y, z, pulse_energy, pulse_length, mean_z, stddev_x, stddev_y): + stddev_z = pulse_length * c + return (pulse_energy / (sqrt((2 * pi) ** 3) * stddev_x * stddev_y * stddev_z) * + exp(-1 / 2 * ((x / stddev_x) ** 2 + (y / stddev_y) ** 2 + ((z - mean_z) / stddev_z) ** 2))) + + +def _constant_bivariate_gaussian2d(x, y, pulse_energy, pulse_length, stddev_x, stddev_y): + length = pulse_length * c + normalisation = pulse_energy / length + return normalisation / (stddev_x * stddev_y * 2 * pi) * exp(-1 / 2 * ((x / stddev_x) ** 2 + (y / stddev_y) ** 2)) diff --git a/demos/laser/laser_profile.py b/demos/laser/laser_profile.py new file mode 100644 index 00000000..fd5a08ad --- /dev/null +++ b/demos/laser/laser_profile.py @@ -0,0 +1,86 @@ +from random import gauss +from raysect.core import Vector3D + +from cherab.core.math.samplers import sample3d_grid +from cherab.core.model.laser.profile import (UniformEnergyDensity, ConstantBivariateGaussian, + TrivariateGaussian, GaussianBeamAxisymmetric) + +import numpy as np +import matplotlib.pyplot as plt + +def plot_profiles(profile, title=""): + """Plot the laser profile + + Produces 2d energy density plot in the x-z plane at y=0. + Produces plots of energy density profiles along x, y and z + directions. + """ + + # set coordinate grid and zero indices + x = np.linspace(-0.01, 0.01, 31) + z = np.linspace(-1, 1, 101) + x_zero = np.abs(x).argmin() + z_zero = np.abs(z).argmin() + + e_xz_profile = sample3d_grid(profile.get_energy_density, x, [0], z)[:, 0, :] + e_y_profile = sample3d_grid(profile.get_energy_density, [0], x, [0])[0, :, 0] + + # plot + fig = plt.figure(constrained_layout=True) + fig.suptitle(title) + spec = fig.add_gridspec(4, 3) + + # profile along y at z=0 + axz = fig.add_subplot(spec[0, 0:2]) + axz.plot(x, e_y_profile, color="C2", ls="dashed") + axz.set_xlabel("y [m]") + axz.set_ylabel("Energy [J/m^3]") + + # x-y plane + ax2d = fig.add_subplot(spec[1:3, 0:2]) + im = ax2d.pcolormesh(x, z, e_xz_profile.T) + fig.colorbar(im, ax=ax2d, label="Energy [J/m^3]") + ax2d.axhline(z[z_zero], color="C0", ls="dashed") + ax2d.axvline(x[x_zero], color="C1", ls="dashed") + ax2d.plot([x[x_zero]], z[z_zero], marker="x", color="C2") + ax2d.set_xlabel("x [m]") + ax2d.set_ylabel("z [m]") + + # profile along z at x=0 + axz = fig.add_subplot(spec[1:3, -1]) + axz.plot(e_xz_profile[x_zero, :], z, color="C0", ls="dashed") + axz.set_ylabel("z [m]") + axz.set_xlabel("Energy [J/m^3]") + + # profile along x at z=0 + axz = fig.add_subplot(spec[-1, 0:2]) + axz.plot(x, e_xz_profile[:, z_zero].T, color="C1", ls="dashed") + axz.set_xlabel("x [m]") + axz.set_ylabel("Energy [J/m^3]") + +# UniformEnergyDensity profile has constant parameters in the whole x, y, z space +uniform = UniformEnergyDensity(energy_density=1e3, laser_length=2., laser_radius=0.005, + polarization=Vector3D(0, 1, 0)) +plot_profiles(uniform, "Uniform Energy Profile") + +# Example of ConstantBivariteGaussian. With energy 2J, pulse temporal length 5 ns +# and different standard deviations in the x and y plane. The mean value in both +# dimensions is equal to 0. There is no correlation between x and y. + +gauss2d = ConstantBivariateGaussian(pulse_energy=2, pulse_length=5e-9, + stddev_x=2e-3, stddev_y=4e-3) +plot_profiles(gauss2d, "ConstantBivariateGaussian Energy Profile") + +# Example of TrivariageGaussian profile. With the pulse mean at z=0.2 +gauss3d = TrivariateGaussian(pulse_energy=2, pulse_length=2e-9, + mean_z=0.2, stddev_x=2e-3, stddev_y=4e-3 ) +plot_profiles(gauss3d, "TrivariateGaussian Energy Profile") + +# Example of GaussianBeamModel with wavelength 1e4 nm and 1 mm standard deviation +# in the waist. The wavelength is too high +# but was selected to produce nice plots in the defined x, y, z dimensions. + +gaussbeam = GaussianBeamAxisymmetric(stddev_waist=1e-3, waist_z=0.2, laser_wavelength=1e4) +plot_profiles(gaussbeam, "GaussianBeam Energy Profile") + +plt.show() \ No newline at end of file diff --git a/demos/laser/laser_spectrum.py b/demos/laser/laser_spectrum.py new file mode 100644 index 00000000..7904bbdc --- /dev/null +++ b/demos/laser/laser_spectrum.py @@ -0,0 +1,34 @@ +from cherab.core.model.laser import ConstantSpectrum, GaussianSpectrum + +import matplotlib.pyplot as plt + + +# construct a ConstantSpectrum with 10 spectral bins +constant_wide = ConstantSpectrum(min_wavelength=1059.9, max_wavelength=1060.1, bins=10) + +# plot the power_spectral_density attribute of the laser +_, ax = plt.subplots() +ax.plot(constant_wide.wavelengths, constant_wide.power_spectral_density) + +ax.set_xlabel("Wavelength [nm]") +ax.set_ylabel("W / nm") + +ax.set_title("Energy Spectral Density") +plt.show() + +# construct a narrow laser spectrum +constant_narrow = ConstantSpectrum(min_wavelength=1059.999, max_wavelength=1060.001, bins=1) +print("narow spectrum wavelengths: {}, power spectral density: {}".format(constant_narrow.wavelengths, + constant_narrow.power_spectral_density)) + +# construct a GaussianSpectrum with 20 bins +gaussian = GaussianSpectrum(min_wavelength=1059, max_wavelength=1061, bins=30, + mean=1060, stddev=0.3) + +_, ax = plt.subplots() +ax.plot(gaussian.wavelengths, gaussian.power_spectral_density) +ax.set_xlabel("Wavelength [nm]") +ax.set_ylabel("W / nm") + +ax.set_title("Energy Spectral Density") +plt.show() diff --git a/demos/laser/model_seldenmatoba.py b/demos/laser/model_seldenmatoba.py new file mode 100644 index 00000000..0491cd66 --- /dev/null +++ b/demos/laser/model_seldenmatoba.py @@ -0,0 +1,72 @@ +import matplotlib.pyplot as plt + +from raysect.optical.spectrum import Spectrum + +from cherab.core.model.laser import SeldenMatobaThomsonSpectrum +from cherab.core.utility import RecursiveDict + + +density = 3e19 # electron density +temperatures = [100, 2.e3] # electron temperatures in eV +laser_wavelength = 1060 # wavelength of the laser light in nm +laser_energy = 1 # energy density of the laser light in J/m^3 + +# angle between observation direction and electric field +angles_pol = [90, 45] + +# angles between propagation direction and observation direction +angles_obs = [45, 90, 135, 120] + +# define spectrum of observed scattering +min_wavelength = 650 +max_wavelength = 1300 +bins = 1000 + + +# calculate Thomson scattered spectra for the specified combinations of el. properties and observation +model = SeldenMatobaThomsonSpectrum() +scattered = RecursiveDict() + +for te in temperatures: + for ap in angles_pol: + for ao in angles_obs: + spectrum = Spectrum(min_wavelength, max_wavelength, bins) + scattered[te][ap][ao] = model.calculate_spectrum(density, te, laser_energy, laser_wavelength, + ao, ap, spectrum).samples +scattered = scattered.freeze() + +wvls = spectrum.wavelengths + +# plot temperature influence on scattered spectra +ap, ao = 90, 90 +_, ax = plt.subplots() +ax.set_title("Scattered spectrum, pol. angl.={:d}, obs. angl = {:d}".format(ap, ao)) +for te in temperatures: + rad = scattered[te][ap][ao] + ax.plot(wvls, rad, label="Te = {:3.1f} eV".format(te)) +ax.set_xlabel("Wavelength [nm]") +ax.set_ylabel("Spectral Radiance [W/m^3/sr]") +ax.legend() + +# plot influence of angle between polarisation and observation on scattered spectra +te, ao = 100, 90 +_, ax = plt.subplots() +ax.set_title("Scattered spectrum, te = {:3.1f}, obs. angl = {:d}".format(te, ao)) +for ap in angles_pol: + rad = scattered[te][ap][ao] + ax.plot(wvls, rad, label="pol. angl. = {:d} deg".format(ap)) +ax.set_xlabel("Wavelength [nm]") +ax.set_ylabel("Spectral Radiance [W/m^3/sr]") +ax.legend() + +# plot influence of observation angle on scattered spectra +te, ap = 2e3, 90 +_, ax = plt.subplots() +ax.set_title("Scattered spectrum, te = {:3.1f}, obs. angl = {:d}".format(te, ao)) +for ao in angles_obs: + rad = scattered[te][ap][ao] + ax.plot(wvls, rad, label="obs. angl. = {:d} deg".format(ao)) +ax.set_xlabel("Wavelength [nm]") +ax.set_ylabel("Spectral Radiance [W/m^3/sr]") +ax.legend() +plt.show() diff --git a/demos/laser/thomson_scattering.py b/demos/laser/thomson_scattering.py new file mode 100644 index 00000000..f9198d8e --- /dev/null +++ b/demos/laser/thomson_scattering.py @@ -0,0 +1,61 @@ +import numpy as np + +from raysect.optical import World, translate, Point3D, rotate_basis, Vector3D +from raysect.optical.observer import FibreOptic + +from cherab.core.model.laser import ConstantBivariateGaussian, ConstantSpectrum, SeldenMatobaThomsonSpectrum +from cherab.core.laser import Laser + +from cherab.generomak.plasma import get_core_plasma +from cherab.generomak.equilibrium import load_equilibrium + +import matplotlib.pyplot as plt + +world = World() + +plasma = get_core_plasma(parent=world) +equilibrium = load_equilibrium() + +# set up the laser +laser = Laser(name="Thomson Scattering laser", parent=world) +laser.transform = translate(equilibrium.magnetic_axis[0], 0, -2) +laser.plasma = plasma +laser.laser_profile = ConstantBivariateGaussian(pulse_energy=2, pulse_length=1e-8, + laser_length=4, laser_radius=1e-2, + stddev_x=3e-3, stddev_y=2e-3) +laser.laser_spectrum = ConstantSpectrum(min_wavelength=1059.9, max_wavelength=1060.1, bins=1) +laser.models = [SeldenMatobaThomsonSpectrum()] + +# generate points on laser in world space to measure on +laser_points = [Point3D(0, 0, round(z, 2)).transform(laser.to_root()) for z in np.linspace(2, 2.8, 5)] +te = [plasma.electron_distribution.effective_temperature(*point) for point in laser_points] +ne = [plasma.electron_distribution.density(*point) for point in laser_points] + +# place fibres 0.1m outside sep on lfs +fibre_position = Point3D(equilibrium.psin_to_r(1) + 0.1, 0, 0) + +# generate fibres list and observe +fibres = {} +for point in laser_points: + + direction = fibre_position.vector_to(point) + transform = translate(*fibre_position) * rotate_basis(direction, Vector3D(0, 0, 1)) + + fibre = FibreOptic(radius=1e-3, acceptance_angle=0.25, + parent=world, transform=transform, + min_wavelength=800, max_wavelength=1200, + spectral_bins=1000, + pixel_samples=1000) + + fibre.pipelines[0].display_progress = False + fibre.observe() + fibres[point.z] = fibre + +_, ax = plt.subplots() +for z, fibre in fibres.items(): + pipeline = fibre.pipelines[0] + ax.plot(pipeline.wavelengths, pipeline.samples.mean, label="z={:1.2f}m".format(z)) +ax.set_xlabel("wavelength [nm]") +ax.set_ylabel("spectral power W/nm") +ax.legend() +plt.show() diff --git a/docs/source/models/emission_models.rst b/docs/source/models/emission_models.rst index 73f76cc2..3a2c3b9c 100644 --- a/docs/source/models/emission_models.rst +++ b/docs/source/models/emission_models.rst @@ -13,3 +13,4 @@ Cherab contains a number of independent physics models for spectroscopic emissio brem/bremsstrahlung basic_line/basic_line_emission line_shapes/spectral_line_shapes + laser/laser diff --git a/docs/source/models/laser/laser.rst b/docs/source/models/laser/laser.rst new file mode 100644 index 00000000..da9f7f79 --- /dev/null +++ b/docs/source/models/laser/laser.rst @@ -0,0 +1,28 @@ +Laser +====== + +Laser Spectrum +-------------- + +.. autoclass:: cherab.core.model.laser.laserspectrum.ConstantSpectrum + :members: +.. autoclass:: cherab.core.model.laser.laserspectrum.GaussianSpectrum + :members: + +Laser Profile +------------- + +.. autoclass:: cherab.core.model.laser.profile.UniformEnergyDensity + :members: +.. autoclass:: cherab.core.model.laser.profile.ConstantBivariateGaussian + :members: +.. autoclass:: cherab.core.model.laser.profile.TrivariateGaussian + :members: +.. autoclass:: cherab.core.model.laser.profile.GaussianBeamAxisymmetric + :members: + +Laser Model +------------- + +.. autoclass:: cherab.core.model.laser.model.SeldenMatobaThomsonSpectrum + :members: diff --git a/docs/source/plasmas/laser.rst b/docs/source/plasmas/laser.rst new file mode 100644 index 00000000..6c94c6d8 --- /dev/null +++ b/docs/source/plasmas/laser.rst @@ -0,0 +1,20 @@ +Laser +====== + +Laser Spectrum +-------------- + +.. autoclass:: cherab.core.laser.laserspectrum.LaserSpectrum + :members: + +Laser Profile +------------- + +.. autoclass:: cherab.core.laser.profile.LaserProfile + :members: + +Laser Model +------------- + +.. autoclass:: cherab.core.laser.model.LaserModel + :members: diff --git a/docs/source/plasmas/plasmas.rst b/docs/source/plasmas/plasmas.rst index 07198e71..98e3d9c8 100644 --- a/docs/source/plasmas/plasmas.rst +++ b/docs/source/plasmas/plasmas.rst @@ -7,3 +7,4 @@ Plasmas core_plasma_classes equilibrium particle_beams + laser