diff --git a/CHANGELOG.md b/CHANGELOG.md index dbe69a21..f68f62e3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,8 @@ New: * Add verbose parameter to SartOpencl solver (default is False). (#358) * Add Generomak core plasma profiles. (#360) * Add toroidal_mesh_from_polygon for making mesh for not fully-360 degrees axisymmetric elements. (#365) +* Add common spectroscopic instruments: Polychromator, SurveySpectrometer, CzernyTurnerSpectrometer. (#299) + Bug Fixes: ---------- diff --git a/cherab/tools/spectroscopy/__init__.py b/cherab/tools/spectroscopy/__init__.py new file mode 100644 index 00000000..23bbbe1f --- /dev/null +++ b/cherab/tools/spectroscopy/__init__.py @@ -0,0 +1,22 @@ + +# 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 .instrument import SpectroscopicInstrument +from .polychromator import PolychromatorFilter, TrapezoidalFilter, Polychromator +from .spectrometer import Spectrometer, CzernyTurnerSpectrometer diff --git a/cherab/tools/spectroscopy/instrument.py b/cherab/tools/spectroscopy/instrument.py new file mode 100644 index 00000000..0e0d666c --- /dev/null +++ b/cherab/tools/spectroscopy/instrument.py @@ -0,0 +1,118 @@ + +# 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. + +class SpectroscopicInstrument: + """ + Base class for spectroscopic instruments (spectrometers, polychromators, etc.). + This is an abstract class. + + :param str name: Instrument name. + + :ivar list pipeline_classes: The list of pipeline classes used with this instrument. + :ivar list pipeline_kwargs: The list of dicts with keywords passed to init methods of + pipeline classes used with this instrument. + :ivar float min_wavelength: Lower wavelength bound for spectral range. + :ivar float max_wavelength: Upper wavelength bound for spectral range. + :ivar int spectral_bins: The number of spectral samples over the wavelength range. + """ + + def __init__(self, name=''): + self._pipeline_classes = None + self.name = name + self._clear_spectral_settings() + + @property + def name(self): + # Instrument name. + return self._name + + @name.setter + def name(self, value): + self._name = str(value) + self._pipeline_kwargs = None + + @property + def pipeline_classes(self): + # The list of pipeline classes used with this instrument. + if self._pipeline_classes is None: + self._update_pipeline_classes() + + return self._pipeline_classes + + @property + def pipeline_kwargs(self): + # The list of dicts with keywords passed to init methods of + # pipeline classes used with this instrument. + if self._pipeline_kwargs is None: + self._update_pipeline_kwargs() + + return self._pipeline_kwargs + + def create_pipelines(self): + """ Returns a list of new pipelines created according to `pipeline_classes` + and keyword arguments.""" + if self._pipeline_classes is None: + self._update_pipeline_classes() + if self._pipeline_kwargs is None: + self._update_pipeline_kwargs() + + pipelines = [] + for PipelineClass, kwargs in zip(self._pipeline_classes, self._pipeline_kwargs): + pipeline = PipelineClass(**kwargs) + pipelines.append(pipeline) + + return pipelines + + @property + def min_wavelength(self): + # Lower wavelength bound for spectral range. + if self._min_wavelength is None: + self._update_spectral_settings() + + return self._min_wavelength + + @property + def max_wavelength(self): + # Upper wavelength bound for spectral range. + if self._max_wavelength is None: + self._update_spectral_settings() + + return self._max_wavelength + + @property + def spectral_bins(self): + # The number of spectral samples over the wavelength range. + if self._spectral_bins is None: + self._update_spectral_settings() + + return self._spectral_bins + + def _clear_spectral_settings(self): + self._min_wavelength = None + self._max_wavelength = None + self._spectral_bins = None + + def _update_spectral_settings(self): + raise NotImplementedError("To be defined in subclass.") + + def _update_pipeline_classes(self): + raise NotImplementedError("To be defined in subclass.") + + def _update_pipeline_kwargs(self): + raise NotImplementedError("To be defined in subclass.") diff --git a/cherab/tools/spectroscopy/polychromator.py b/cherab/tools/spectroscopy/polychromator.py new file mode 100644 index 00000000..56fb6e18 --- /dev/null +++ b/cherab/tools/spectroscopy/polychromator.py @@ -0,0 +1,221 @@ + +# 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. + +import numpy as np +from raysect.optical import InterpolatedSF +from raysect.optical.observer import RadiancePipeline0D + +from .instrument import SpectroscopicInstrument + + +class PolychromatorFilter(InterpolatedSF): + """ + Defines a polychromator filter as a Raysect's InterpolatedSF. + + :param object wavelengths: 1D array of wavelengths in nanometers. + :param object samples: 1D array of spectral samples. + :param bool normalise: True/false toggle for whether to normalise the + spectral function so its integral equals 1. + :param str name: Filter name (e.g. "H-alpha filter"). Default is ''. + + :ivar float min_wavelength: Lower wavelength bound of the filter's spectral range in nm. + :ivar float max_wavelength: Upper wavelength bound of the filter's spectral range in nm. + """ + + def __init__(self, wavelengths, samples, normalise=False, name=''): + + wavelengths = np.array(wavelengths, dtype=np.float64) + samples = np.array(samples, dtype=np.float64) + + if wavelengths.ndim != 1: + raise ValueError("Wavelength array must be 1D.") + + if samples.shape[0] != wavelengths.shape[0]: + raise ValueError("Wavelength and sample arrays must be the same length.") + + indices = np.argsort(wavelengths) + wavelengths = wavelengths[indices] + samples = samples[indices] + + self._min_wavelength = wavelengths[0] + self._max_wavelength = wavelengths[-1] + self._window = self._max_wavelength - self._min_wavelength + self._central_wavelength = 0.5 * (self._max_wavelength + self._min_wavelength) + + # setting the ends of the filter to zero, if they are not + if samples[0] != 0: + wavelengths = np.insert(wavelengths, 0, wavelengths[0] * (1. - 1.e-15)) + samples = np.insert(samples, 0, 0) + if samples[-1] != 0: + wavelengths = np.append(wavelengths, wavelengths[-1] * (1. + 1.e-15)) + samples = np.append(samples, 0) + + super().__init__(wavelengths, samples, normalise) + self._name = str(name) + + @property + def name(self): + # Filter name. + return self._name + + @property + def min_wavelength(self): + # Lower wavelength bound of the filter's spectral range in nm. + return self._min_wavelength + + @property + def max_wavelength(self): + # Upper wavelength bound of the filter's spectral range in nm. + return self._max_wavelength + + @property + def window(self): + # Size of the filtering window in nm. + return self._window + + @property + def central_wavelength(self): + # Central wavelength of the filter in nm. + return self._central_wavelength + + +class TrapezoidalFilter(PolychromatorFilter): + """ + Symmetrical trapezoidal polychromator filter. + + :param float wavelength: Central wavelength of the filter in nm. + :param float window: Size of the filtering window in nm. Default is 3. + :param float flat_top: Size of the flat top part of the filter in nm. + Default is None (equal to window). + :param str name: Filter name (e.g. "H-alpha filter"). Default is ''. + """ + + def __init__(self, central_wavelength, window=3., flat_top=None, name=''): + + if central_wavelength <= 0: + raise ValueError("Argument 'central_wavelength' must be positive.") + + if window <= 0: + raise ValueError("Argument 'window' must be positive.") + + flat_top = flat_top or window + + if flat_top <= 0: + raise ValueError("Argument 'flat_top' must be positive.") + if flat_top > window: + raise ValueError("Argument 'flat_top' must be less or equal than 'window'.") + + self._flat_top = flat_top + + if flat_top == window: + flat_top -= flat_top * 1.e-15 + + wavelengths = [central_wavelength - 0.5 * window, + central_wavelength - 0.5 * flat_top, + central_wavelength + 0.5 * flat_top, + central_wavelength + 0.5 * window] + samples = [0, 1, 1, 0] + super().__init__(wavelengths, samples, normalise=False, name=name) + + @property + def flat_top(self): + # Size of the flat top part of the filter in nm. + return self._flat_top + + +class Polychromator(SpectroscopicInstrument): + """ + A polychromator assembly with a set of different filters. + + :param list filters: List of the `PolychromatorFilter` instances. + :param int min_bins_per_window: Minimal number of spectral bins + per filtering window. Default is 10. + :param str name: Polychromator name. + + .. code-block:: pycon + + >>> from raysect.optical import World + >>> from raysect.optical.observer import FibreOptic + >>> from cherab.tools.spectroscopy import Polychromator, TrapezoidalFilter + >>> + >>> world = World() + >>> h_alpha_filter = TrapezoidalFilter(656.1, name='H-alpha filter') + >>> ciii_465nm_filter = TrapezoidalFilter(464.8, name='CIII 465 nm filter') + >>> polychromator = Polychromator([h_alpha_filter, ciii_465nm_filter], name='MyPolychromator') + >>> fibreoptic = FibreOptic(name="MyFibreOptic", parent=world) + >>> fibreoptic.min_wavelength = polychromator.min_wavelength + >>> fibreoptic.max_wavelength = polychromator.max_wavelength + >>> fibreoptic.spectral_bins = polychromator.spectral_bins + >>> fibreoptic.pipelines = polychromator.create_pipelines() + """ + + def __init__(self, filters, min_bins_per_window=10, name=''): + super().__init__(name) + self.min_bins_per_window = min_bins_per_window + self.filters = filters + + @property + def min_bins_per_window(self): + # Minimal number of spectral bins per filtering window. + return self._min_bins_per_window + + @min_bins_per_window.setter + def min_bins_per_window(self, value): + value = int(value) + if value <= 0: + raise ValueError("Attribute 'min_bins_per_window' must be positive.") + + self._min_bins_per_window = value + self._clear_spectral_settings() + + @property + def filters(self): + # List of the PolychromatorFilter instances. + return self._filters + + @filters.setter + def filters(self, value): + for poly_filter in value: + if not isinstance(poly_filter, PolychromatorFilter): + raise TypeError('Property filters must contain only PolychromatorFilter instances.') + + self._filters = value + self._clear_spectral_settings() + self._pipeline_classes = None + self._pipeline_kwargs = None + + def _update_pipeline_classes(self): + self._pipeline_classes = [RadiancePipeline0D for poly_filter in self._filters] + + def _update_pipeline_kwargs(self): + self._pipeline_kwargs = [{'name': self._name + ': ' + poly_filter.name, 'filter': poly_filter} for poly_filter in self._filters] + + def _update_spectral_settings(self): + + min_wavelength = np.inf + max_wavelength = 0 + step = np.inf + for poly_filter in self._filters: + step = min(step, poly_filter.window / self._min_bins_per_window) + min_wavelength = min(min_wavelength, poly_filter.min_wavelength) + max_wavelength = max(max_wavelength, poly_filter.max_wavelength) + + self._min_wavelength = min_wavelength + self._max_wavelength = max_wavelength + self._spectral_bins = int(np.ceil((max_wavelength - min_wavelength) / step)) diff --git a/cherab/tools/spectroscopy/spectrometer.py b/cherab/tools/spectroscopy/spectrometer.py new file mode 100644 index 00000000..585ac3f4 --- /dev/null +++ b/cherab/tools/spectroscopy/spectrometer.py @@ -0,0 +1,351 @@ + +# 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. + +import numpy as np +from raysect.optical import Spectrum +from raysect.optical.observer import SpectralRadiancePipeline0D + +from .instrument import SpectroscopicInstrument + + +class Spectrometer(SpectroscopicInstrument): + """ + Spectrometer that can accommodate multiple spectra. + + Spectrometer is initialized with a sequence of calibration arrays (one array per accommodated + spectrum) containing the wavelengths of the pixel borders. Namely, the values + :math:`w_{k}^{i}` and :math:`w_{k}^{i+1}` define the spectral range of the pixel :math:`p_i` + of the `k`-th spectrum. After the spectrum is ray-traced, it can be recalibrated with + `spectrometer.calibrate(spectrum)`. + + Note that Raysect cannot raytrace the spectra with non-constant spectral resolution. + Thus, the actual number of spectral bins of raytraced spectrum is defined with + `min_bins_per_pixel` attribute. + + :param tuple wavelength_to_pixel: Wavelength-to-pixel calibration arrays. + :param int min_bins_per_pixel: Minimal number of spectral bins + per pixel. Default is 1. + :param str name: Spectrometer name. + + :ivar tuple wavelengths: Central wavelengths of the pixels. + + .. code-block:: pycon + + >>> from raysect.optical import World, Spectrum + >>> from raysect.optical.observer import FibreOptic + >>> from cherab.tools.spectroscopy import Spectrometer + >>> from matplotlib import pyplot as plt + >>> + >>> wavelength_to_pixel = ([400., 400.5, 401.5, 402., 404.], + >>> [600., 600.5, 601.5, 602., 604., 607.]) + >>> spectrometer = Spectrometer(wavelength_to_pixel, min_bins_per_pixel=5, + >>> name='MySpectrometer') + >>> + >>> world = World() + >>> fibreoptic = FibreOptic(name="MyFibreOptic", parent=world) + >>> fibreoptic.min_wavelength = spectrometer.min_wavelength + >>> fibreoptic.max_wavelength = spectrometer.max_wavelength + >>> fibreoptic.spectral_bins = spectrometer.spectral_bins + >>> fibreoptic.pipelines = spectrometer.create_pipelines() + >>> ... + >>> fibreoptic.observe() + >>> spectrum = Spectrum(fibreoptic.min_wavelength, fibreoptic.max_wavelength, fibreoptic.spectral_bins) + >>> spectrum.samples[:] = fibreoptic.pipelines[0].mean + >>> calibrated_spectra = spectrometer.calibrate(spectrum) + >>> wavelengths = spectrometer.wavelengths + >>> + >>> plt.plot(wavelengths[0], calibrated_spectra[0]) + >>> plt.show() + """ + + def __init__(self, wavelength_to_pixel, min_bins_per_pixel=1, name=''): + + self.min_bins_per_pixel = min_bins_per_pixel + self.wavelength_to_pixel = wavelength_to_pixel + super().__init__(name) + + @property + def wavelength_to_pixel(self): + # Wavelength-to-pixel calibration arrays. + return self._wavelength_to_pixel + + @wavelength_to_pixel.setter + def wavelength_to_pixel(self, value): + _wavelength_to_pixel = [] + _wavelengths = [] + for wl2pix in value: + wl2pix = np.array(wl2pix, dtype=float) + if wl2pix.ndim != 1: + raise ValueError('Attribute wavelength_to_pixel must only contain one-dimensional arrays.') + if wl2pix.size < 2: + raise ValueError('Attribute wavelength_to_pixel must only contain arrays of at least 2 elements.') + if np.any(np.diff(wl2pix) <= 0): + raise ValueError('Attribute wavelength_to_pixel must only contain monotonically increasing arrays.') + wl2pix.flags.writeable = False + _wavelength_to_pixel.append(wl2pix) + wl_center = 0.5 * (wl2pix[1:] + wl2pix[:-1]) + wl_center.flags.writeable = False + _wavelengths.append(wl_center) + self._wavelength_to_pixel = tuple(_wavelength_to_pixel) + self._wavelengths = tuple(_wavelengths) + self._clear_spectral_settings() + + @property + def wavelengths(self): + # Central wavelengths of the pixels. + return self._wavelengths + + @property + def min_bins_per_pixel(self): + # Minimal number of spectral bins per pixel. + return self._min_bins_per_pixel + + @min_bins_per_pixel.setter + def min_bins_per_pixel(self, value): + value = int(value) + if value <= 0: + raise ValueError("Attribute 'min_bins_per_pixel' must be positive.") + + self._min_bins_per_pixel = value + self._clear_spectral_settings() + + def _update_pipeline_classes(self): + self._pipeline_classes = [SpectralRadiancePipeline0D] + + def _update_pipeline_kwargs(self): + self._pipeline_kwargs = [{'name': self._name}] + + def _update_spectral_settings(self): + self._min_wavelength = min(wl2pix[0] for wl2pix in self._wavelength_to_pixel) + self._max_wavelength = max(wl2pix[-1] for wl2pix in self._wavelength_to_pixel) + step = min(np.diff(wl2pix).min() for wl2pix in self._wavelength_to_pixel) / self._min_bins_per_pixel + self._spectral_bins = int(np.ceil((self._max_wavelength - self._min_wavelength) / step)) + + def calibrate(self, spectrum): + """ + Calibrates the spectrum according to the `wavelength_to_pixel` arrays + by averaging it over the pixel widths. + + :param Spectrum spectrum: Spectrum to calibrate. + + :returns: A tuple of calibrated spectra as ndarrays. + """ + if not isinstance(spectrum, Spectrum): + raise TypeError('Argument spectrum must be a Spectrum instance.') + if spectrum.min_wavelength > self.min_wavelength or spectrum.max_wavelength < self.max_wavelength: + raise ValueError('Unable to calibrate the spectrum. ' + 'The spectrum has narrower range ({}, {}) than the spectrometer ({}, {}).'.format(spectrum.min_wavelength, + spectrum.max_wavelength, + self.min_wavelength, + self.max_wavelength)) + calibrated_spectra = [] + for wl2pix in self.wavelength_to_pixel: + calibrated_spectrum = np.zeros(wl2pix.size - 1) + for i in range(wl2pix.size - 1): + calibrated_spectrum[i] = spectrum.integrate(wl2pix[i], wl2pix[i + 1]) / (wl2pix[i + 1] - wl2pix[i]) + calibrated_spectra.append(calibrated_spectrum) + + return calibrated_spectra + + +class CzernyTurnerSpectrometer(Spectrometer): + """ + Czerny-Turner spectrometer. + + The Czerny-Turner spectrometer is initialized with the parameters of the diffraction scheme + and a sequence of accommodated spectra, each of which is determined by the lower wavelength + bound and the number of pixels. + + This spectrometer automatically fills the wavelength-to-pixel calibration arrays + according to the parameters of the diffraction scheme. + + :param int diffraction_order: Diffraction order. + :param float grating: Diffraction grating in nm-1. + :param float focal_length: Focal length in nm. + :param float pixel_spacing: Pixel to pixel spacing on CCD in nm. + :param float diffraction_angle: Angle between incident and diffracted light in degrees. + :param tuple accommodated_spectra: A sequence of (`min_wavelength`, `pixels`) pairs, specifying + the lower wavelength bound and the number of pixels + of accommodated spectra. + :param int min_bins_per_pixel: Minimal number of spectral bins + per pixel. Default is 1. + :param str name: Spectrometer name. + + :ivar tuple wavelength_to_pixel: Wavelength-to-pixel calibration arrays. + + .. code-block:: pycon + + >>> from raysect.optical import World + >>> from raysect.optical.observer import FibreOptic + >>> from cherab.tools.spectroscopy import CzernyTurnerSpectrometer + >>> + >>> world = World() + >>> hires_spectrometer = CzernyTurnerSpectrometer(1, 2.e-3, 1.e9, 2.e4, 10., + >>> ((600., 512), (700., 128)), + >>> name='MySpectrometer') + >>> fibreoptic = FibreOptic(name="MyFibreOptic", parent=world) + >>> fibreoptic.min_wavelength = hires_spectrometer.min_wavelength + >>> fibreoptic.max_wavelength = hires_spectrometer.max_wavelength + >>> fibreoptic.spectral_bins = hires_spectrometer.spectral_bins + >>> fibreoptic.pipelines = hires_spectrometer.create_pipelines() + """ + + def __init__(self, diffraction_order, grating, focal_length, pixel_spacing, diffraction_angle, + accommodated_spectra, min_bins_per_pixel=1, name=''): + self._accommodated_spectra = None + self.diffraction_order = diffraction_order + self.grating = grating + self.focal_length = focal_length + self.pixel_spacing = pixel_spacing + self.diffraction_angle = diffraction_angle + self.accommodated_spectra = accommodated_spectra + self.min_bins_per_pixel = min_bins_per_pixel + self.name = name + + @property + def diffraction_order(self): + # Diffraction order. + return self._diffraction_order + + @diffraction_order.setter + def diffraction_order(self, value): + value = int(value) + if value <= 0: + raise ValueError("Attribute 'diffraction_order' must be positive.") + + self._diffraction_order = value + # resolution has changed, recalculating wavelength_to_pixel + self._update_wavelength_to_pixel() + + @property + def grating(self): + # Diffraction grating in nm-1. + return self._grating + + @grating.setter + def grating(self, value): + if value <= 0: + raise ValueError("Attribute 'grating' must be positive.") + + self._grating = value + # resolution has changed, recalculating wavelength_to_pixel + self._update_wavelength_to_pixel() + + @property + def focal_length(self): + # Focal length in nm. + return self._focal_length + + @focal_length.setter + def focal_length(self, value): + if value <= 0: + raise ValueError("Attribute 'focal_length' must be positive.") + + self._focal_length = value + # resolution has changed, recalculating wavelength_to_pixel + self._update_wavelength_to_pixel() + + @property + def pixel_spacing(self): + # Pixel to pixel spacing on CCD in nm. + return self._pixel_spacing + + @pixel_spacing.setter + def pixel_spacing(self, value): + if value <= 0: + raise ValueError("Attribute 'pixel_spacing' must be positive.") + + self._pixel_spacing = value + # resolution has changed, recalculating wavelength_to_pixel + self._update_wavelength_to_pixel() + + @property + def diffraction_angle(self): + # Angle between incident and diffracted light in degrees. + return np.rad2deg(self._diffraction_angle) + + @diffraction_angle.setter + def diffraction_angle(self, value): + if value <= 0: + raise ValueError("Attribute 'diffraction_angle' must be positive.") + + self._diffraction_angle = np.deg2rad(value) + # resolution has changed, recalculating wavelength_to_pixel + self._update_wavelength_to_pixel() + + @property + def accommodated_spectra(self): + return self._accommodated_spectra + + @accommodated_spectra.setter + def accommodated_spectra(self, value): + for min_wavelength, pixels in value: + if min_wavelength <= 0: + raise ValueError('The value of min_wavelength in accommodated_spectra must be positive.') + if pixels <= 0: + raise ValueError('The value of pixels in accommodated_spectra must be positive.') + self._accommodated_spectra = value + self._update_wavelength_to_pixel() + + def _update_wavelength_to_pixel(self): + + if self._accommodated_spectra is None: + return + + _wavelength_to_pixel = [] + _wavelengths = [] + for min_wavelength, pixels in self._accommodated_spectra: + pixels = int(pixels) + wl2pix = np.zeros(pixels + 1) + wl2pix[0] = min_wavelength + for i in range(1, pixels + 1): + wl2pix[i] = wl2pix[i - 1] + self.resolution(wl2pix[i - 1]) + wl2pix.flags.writeable = False + _wavelength_to_pixel.append(wl2pix) + wl_center = 0.5 * (wl2pix[1:] + wl2pix[:-1]) + wl_center.flags.writeable = False + _wavelengths.append(wl_center) + self._wavelength_to_pixel = tuple(_wavelength_to_pixel) + self._wavelengths = tuple(_wavelengths) + + self._clear_spectral_settings() + + @property + def wavelength_to_pixel(self): + # Wavelength-to-pixel calibration arrays. + return self._wavelength_to_pixel + + def resolution(self, wavelength): + """ + Calculates spectral resolution in nm for a given wavelength. + + :param wavelength: Wavelength in nm. + + :returns: Resolution in nm. + """ + grating = self._grating + m = self._diffraction_order + dxdp = self._pixel_spacing + angle = self._diffraction_angle + fl = self._focal_length + + p = 0.5 * m * grating * wavelength + _resolution = dxdp * (np.sqrt(np.cos(angle)**2 - p * p) - p * np.tan(angle)) / (m * fl * grating) + + return _resolution diff --git a/cherab/tools/tests/test_spectroscopic_instruments.py b/cherab/tools/tests/test_spectroscopic_instruments.py new file mode 100644 index 00000000..211f7830 --- /dev/null +++ b/cherab/tools/tests/test_spectroscopic_instruments.py @@ -0,0 +1,170 @@ +# 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. + +import unittest +import numpy as np + +from raysect.optical import Spectrum +from raysect.optical.observer.pipeline import RadiancePipeline0D, SpectralRadiancePipeline0D +from cherab.tools.spectroscopy import TrapezoidalFilter, PolychromatorFilter, Polychromator, CzernyTurnerSpectrometer, Spectrometer + + +class TestPolychromatorFilter(unittest.TestCase): + """ + Test for PolychromatorFilter class. + """ + + def test_spectrum(self): + wavelengths = [658, 654, 656] # unsorted + samples = [0.5, 0.5, 1] # non-zero at the ends + poly_filter = PolychromatorFilter(wavelengths, samples, name='test_filter') + wavelengths = np.linspace(653., 659., 7) + spectrum_true = np.array([0, 0.5, 0.75, 1., 0.75, 0.5, 0]) + spectrum_test = np.array([poly_filter(wvl) for wvl in wavelengths]) + self.assertTrue(np.all(spectrum_true == spectrum_test)) + + +class TestTrapezoidalFilter(unittest.TestCase): + """ + Test for TrapezoidalFilter class. + """ + + def test_spectrum(self): + wavelength = 500. + window = 6. + flat_top = 2. + poly_filter = TrapezoidalFilter(wavelength, window, flat_top, 'test_filter') + wavelengths = np.linspace(496., 504., 9) + spectrum_true = np.array([0, 0, 0.5, 1., 1., 1., 0.5, 0, 0]) + spectrum_test = np.array([poly_filter(wvl) for wvl in wavelengths]) + self.assertTrue(np.all(spectrum_true == spectrum_test)) + + +class TestPolychromator(unittest.TestCase): + """ + Test cases for Polychromator class. + """ + + poly_filters_default = (TrapezoidalFilter(400., 6., 2., 'filter 1'), + TrapezoidalFilter(700., 8., 4., 'filter 2')) + min_bins_per_window_default = 10 + + def test_pipeline_classes(self): + polychromator = Polychromator(self.poly_filters_default, self.min_bins_per_window_default, 'test polychromator') + pipeline_classes_true = [RadiancePipeline0D, RadiancePipeline0D] + self.assertSequenceEqual(pipeline_classes_true, polychromator.pipeline_classes) + + def test_pipeline_kwargs(self): + polychromator = Polychromator(self.poly_filters_default, self.min_bins_per_window_default, 'test polychromator') + pipeline_kwargs_true = [{'name': 'test polychromator: filter 1', 'filter': self.poly_filters_default[0]}, + {'name': 'test polychromator: filter 2', 'filter': self.poly_filters_default[1]}] + self.assertSequenceEqual(pipeline_kwargs_true, polychromator.pipeline_kwargs) + + def test_spectral_properties(self): + polychromator = Polychromator(self.poly_filters_default, self.min_bins_per_window_default) + min_wavelength_true = 397. + max_wavelength_true = 704. + spectral_bins_true = 512 + self.assertTrue(polychromator.min_wavelength == min_wavelength_true and + polychromator.max_wavelength == max_wavelength_true and + polychromator.spectral_bins == spectral_bins_true) + + def test_filter_change(self): + """ Checks if the spectral properties are updated correctly when the filters are replaced.""" + polychromator = Polychromator(self.poly_filters_default, self.min_bins_per_window_default) + polychromator.min_bins_per_window = 20 + polychromator.filters = [TrapezoidalFilter(500., 5., 2., 'filter 1'), + TrapezoidalFilter(600., 7., 4., 'filter 2')] + min_wavelength_true = 497.5 + max_wavelength_true = 603.5 + spectral_bins_true = 424 + self.assertTrue(polychromator.min_wavelength == min_wavelength_true and + polychromator.max_wavelength == max_wavelength_true and + polychromator.spectral_bins == spectral_bins_true) + + +class TestSpectrometer(unittest.TestCase): + """ + Test cases for Spectrometer class. + """ + + def test_pipeline_classes(self): + wavelength_to_pixel = ([400., 400.5],) + spectrometer = Spectrometer(wavelength_to_pixel, name='test spectrometer') + self.assertSequenceEqual([SpectralRadiancePipeline0D], spectrometer.pipeline_classes) + + def test_pipeline_kwargs(self): + wavelength_to_pixel = ([400., 400.5],) + spectrometer = Spectrometer(wavelength_to_pixel, name='test spectrometer') + self.assertSequenceEqual([{'name': 'test spectrometer'}], spectrometer.pipeline_kwargs) + + def test_spectral_properties(self): + wavelength_to_pixel = ([400., 400.5, 401.5, 402., 404.], [600., 600.5, 601.5, 602., 604., 607.]) + spectrometer = Spectrometer(wavelength_to_pixel, min_bins_per_pixel=2, name='test spectrometer') + min_wavelength_true = 400. + max_wavelength_true = 607. + spectra_bins_true = 828 + self.assertTrue(spectrometer.min_wavelength == min_wavelength_true and + spectrometer.max_wavelength == max_wavelength_true and + spectrometer.spectral_bins == spectra_bins_true) + + def test_calibration(self): + wavelength_to_pixel = ([400., 400.5, 401.5, 402., 404.],) + spectrometer = Spectrometer(wavelength_to_pixel, name='test spectrometer') + spectrum = Spectrum(399, 405, 12) + s, ds = np.linspace(0, 6., 13, retstep=True) + spectrum.samples[:] = s[:-1] + 0.5 * ds + calibrated_spectra = spectrometer.calibrate(spectrum) + self.assertTrue(np.all(calibrated_spectra[0] == np.array([1.25, 2., 2.75, 4.]))) + + +class TestCzernyTurnerSpectrometer(unittest.TestCase): + """ + Test cases for CzernyTurnerSpectrometer class. + """ + + diffraction_order = 1 + grating = 2.e-3 + focal_length = 1.e9 + pixel_spacing = 2.e4 + diffraction_angle = 10. + accommodated_spectra = ((400., 64), (500., 32)) + min_bins_per_pixel = 2 + + def test_resolution(self): + wavelengths = np.array([350., 550., 750.]) + resolutions_true = np.array([8.587997e-3, 7.199328e-3, 5.0599164e-3]) + spectrometer = CzernyTurnerSpectrometer(self.diffraction_order, self.grating, self.focal_length, self.pixel_spacing, + self.diffraction_angle, self.accommodated_spectra, name='test spectrometer') + resolutions = spectrometer.resolution(wavelengths) + self.assertTrue(np.all(np.abs(resolutions / resolutions_true - 1.) < 1.e-7)) + + def test_spectral_properties(self): + min_wavelength_true = 400 + max_wavelength_true = 500.24326 + spectra_bins_true = 26377 + spectrometer = CzernyTurnerSpectrometer(self.diffraction_order, self.grating, self.focal_length, self.pixel_spacing, + self.diffraction_angle, self.accommodated_spectra, + min_bins_per_pixel=self.min_bins_per_pixel, name='test spectrometer') + self.assertTrue(spectrometer.min_wavelength == min_wavelength_true and + spectrometer.spectral_bins == spectra_bins_true and + abs(spectrometer.max_wavelength - max_wavelength_true) < 1.e-5) + + +if __name__ == '__main__': + unittest.main() diff --git a/docs/source/tools/observers.rst b/docs/source/tools/observers.rst index e9a3a7e8..e93e2057 100644 --- a/docs/source/tools/observers.rst +++ b/docs/source/tools/observers.rst @@ -146,3 +146,4 @@ Spectroscopic prefix in class name. .. autoclass:: cherab.tools.observers.SpectroscopicFibreOptic :members + diff --git a/docs/source/tools/spectroscopy.rst b/docs/source/tools/spectroscopy.rst new file mode 100644 index 00000000..a116dd77 --- /dev/null +++ b/docs/source/tools/spectroscopy.rst @@ -0,0 +1,42 @@ + +Spectroscopy +============ + +The tools for plasma spectroscopy. + +.. _spectroscopy_instruments: + +Spectroscopic instruments +------------------------- + +Spectroscopic instruments such as polychromators and spectrometers +simplify the setup of properties of the observers and rendering pipelines. The instruments +are not connected to the scenegraph, so they cannot observe the world. However, the instruments +have properties, such as `min_wavelength`, `max_wavelength`, `spectral_bins`, +`pipeline_properties`, with which the observer can be configured. +The Cherab core package provides base classes for spectroscopic instruments, +so machine-specific packages can build more advance instruments from them, such as instruments +with spectral properties based on the actual experimental setup for a given shot/pulse. + +.. autoclass:: cherab.tools.spectroscopy.SpectroscopicInstrument + :members: + +.. autoclass:: cherab.tools.spectroscopy.PolychromatorFilter + :members: + +.. autoclass:: cherab.tools.spectroscopy.TrapezoidalFilter + :show-inheritance: + :members: + +.. autoclass:: cherab.tools.spectroscopy.Polychromator + :show-inheritance: + :members: + +.. autoclass:: cherab.tools.spectroscopy.Spectrometer + :show-inheritance: + :members: + +.. autoclass:: cherab.tools.spectroscopy.CzernyTurnerSpectrometer + :show-inheritance: + :members: + diff --git a/docs/source/tools/tools.rst b/docs/source/tools/tools.rst index 3b356121..19f58e95 100644 --- a/docs/source/tools/tools.rst +++ b/docs/source/tools/tools.rst @@ -8,6 +8,7 @@ Tools materials primitives observers + spectroscopy tomography utility