diff --git a/CHANGELOG.md b/CHANGELOG.md index 1f61bd5a..86a64e3a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,13 @@ Project Changelog Release 1.5.0 (TBD) ------------------- +API changes: +* The line shape models are moved to a dedicated submodule. The user code should not be affected though. (#396) +* The line shape models now have AtomicData as a required parameter. +* The method show_supported_transitions() of StarkBroadenedLine and ParametrisedZeemanTriplet is removed. +* The argument stark_model_coefficients of StarkBroadenedLine is now a tuple instead of a dict. +* The argument line_parameters of ParametrisedZeemanTriplet is now a tuple instead of a dict. + New: * Support Raysect 0.8 * Add custom line shape support to BeamCXLine model. (#394) @@ -16,6 +23,7 @@ New: * **Beam dispersion calculation has changed from sigma(z) = sigma + z * tan(alpha) to sigma(z) = sqrt(sigma^2 + (z * tan(alpha))^2) for consistancy with the Gaussian beam model. Attention!!! The results of BES and CX spectroscopy are affected by this change. (#414)** * Improved beam direction calculation to allow for natural broadening of the BES line shape due to beam divergence. (#414) * Add kwargs to invert_regularised_nnls to pass them to scipy.optimize.nnls. (#438) +* StarkBroadenedLine now supports Doppler broadening and Zeeman splitting. (#393) * Add the power radiated in spectral lines due to charge exchange with thermal neutral hydrogen to the TotalRadiatedPower model. (#370) * Add thermal charge-exchange emission model. (#57) * PECs for C VI spectral lines for n <= 5 are now included in populate(). Rerun populate() after upgrading to 1.5 to update the atomic data repository. diff --git a/cherab/core/atomic/data/lineshape/stark/d.json b/cherab/core/atomic/data/lineshape/stark/d.json new file mode 100644 index 00000000..c0092760 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/stark/d.json @@ -0,0 +1,70 @@ +{ + "0": { + "3 -> 2": [ + 3.71e-18, + 0.7665, + 0.064 + ], + "4 -> 2": [ + 8.425e-18, + 0.7803, + 0.050 + ], + "5 -> 2": [ + 1.31e-15, + 0.6796, + 0.030 + ], + "6 -> 2": [ + 3.954e-16, + 0.7149, + 0.028 + ], + "7 -> 2": [ + 6.258e-16, + 0.712, + 0.029 + ], + "8 -> 2": [ + 7.378e-16, + 0.7159, + 0.032 + ], + "9 -> 2": [ + 8.947e-16, + 0.7177, + 0.033 + ], + "4 -> 3": [ + 1.330e-16, + 0.7449, + 0.045 + ], + "5 -> 3": [ + 6.64e-16, + 0.7356, + 0.044 + ], + "6 -> 3": [ + 2.481e-15, + 0.7118, + 0.016 + ], + "7 -> 3": [ + 3.270e-15, + 0.7137, + 0.029 + ], + "8 -> 3": [ + 4.343e-15, + 0.7133, + 0.032 + ], + "9 -> 3": [ + 5.588e-15, + 0.7165, + 0.033 + ] + }, + "reference": "B. Lomanowski, et al. Inferring divertor plasma properties from hydrogen Balmer and Paschen series spectroscopy in JET-ILW. Nuclear Fusion 55.12 (2015) 123028" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/stark/h.json b/cherab/core/atomic/data/lineshape/stark/h.json new file mode 100644 index 00000000..c0092760 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/stark/h.json @@ -0,0 +1,70 @@ +{ + "0": { + "3 -> 2": [ + 3.71e-18, + 0.7665, + 0.064 + ], + "4 -> 2": [ + 8.425e-18, + 0.7803, + 0.050 + ], + "5 -> 2": [ + 1.31e-15, + 0.6796, + 0.030 + ], + "6 -> 2": [ + 3.954e-16, + 0.7149, + 0.028 + ], + "7 -> 2": [ + 6.258e-16, + 0.712, + 0.029 + ], + "8 -> 2": [ + 7.378e-16, + 0.7159, + 0.032 + ], + "9 -> 2": [ + 8.947e-16, + 0.7177, + 0.033 + ], + "4 -> 3": [ + 1.330e-16, + 0.7449, + 0.045 + ], + "5 -> 3": [ + 6.64e-16, + 0.7356, + 0.044 + ], + "6 -> 3": [ + 2.481e-15, + 0.7118, + 0.016 + ], + "7 -> 3": [ + 3.270e-15, + 0.7137, + 0.029 + ], + "8 -> 3": [ + 4.343e-15, + 0.7133, + 0.032 + ], + "9 -> 3": [ + 5.588e-15, + 0.7165, + 0.033 + ] + }, + "reference": "B. Lomanowski, et al. Inferring divertor plasma properties from hydrogen Balmer and Paschen series spectroscopy in JET-ILW. Nuclear Fusion 55.12 (2015) 123028" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/stark/t.json b/cherab/core/atomic/data/lineshape/stark/t.json new file mode 100644 index 00000000..c0092760 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/stark/t.json @@ -0,0 +1,70 @@ +{ + "0": { + "3 -> 2": [ + 3.71e-18, + 0.7665, + 0.064 + ], + "4 -> 2": [ + 8.425e-18, + 0.7803, + 0.050 + ], + "5 -> 2": [ + 1.31e-15, + 0.6796, + 0.030 + ], + "6 -> 2": [ + 3.954e-16, + 0.7149, + 0.028 + ], + "7 -> 2": [ + 6.258e-16, + 0.712, + 0.029 + ], + "8 -> 2": [ + 7.378e-16, + 0.7159, + 0.032 + ], + "9 -> 2": [ + 8.947e-16, + 0.7177, + 0.033 + ], + "4 -> 3": [ + 1.330e-16, + 0.7449, + 0.045 + ], + "5 -> 3": [ + 6.64e-16, + 0.7356, + 0.044 + ], + "6 -> 3": [ + 2.481e-15, + 0.7118, + 0.016 + ], + "7 -> 3": [ + 3.270e-15, + 0.7137, + 0.029 + ], + "8 -> 3": [ + 4.343e-15, + 0.7133, + 0.032 + ], + "9 -> 3": [ + 5.588e-15, + 0.7165, + 0.033 + ] + }, + "reference": "B. Lomanowski, et al. Inferring divertor plasma properties from hydrogen Balmer and Paschen series spectroscopy in JET-ILW. Nuclear Fusion 55.12 (2015) 123028" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/zeeman/parametrised/b.json b/cherab/core/atomic/data/lineshape/zeeman/parametrised/b.json new file mode 100644 index 00000000..b63c7c85 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/zeeman/parametrised/b.json @@ -0,0 +1,35 @@ +{ + "4": { + "6 -> 5": [ + 0.0083423, + 2.0519, + -0.2960 + ], + "7 -> 6": [ + 0.0228379, + 1.6546, + -0.2941 + ], + "8 -> 6": [ + 0.0084065, + 1.8041, + -0.3177 + ], + "8 -> 7": [ + 0.0541883, + 1.4128, + -0.2966 + ], + "9 -> 7": [ + 0.0190781, + 1.5440, + -0.3211 + ], + "10 -> 8": [ + 0.0391914, + 1.3569, + -0.3252 + ] + }, + "reference": "A. Blom and C. Jupén. Parametrisation of the Zeeman effect for hydrogen-like spectra in high-temperature plasmas. Plasma Phys. Control. Fusion 44 (2002) 1229-1241" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/zeeman/parametrised/be.json b/cherab/core/atomic/data/lineshape/zeeman/parametrised/be.json new file mode 100644 index 00000000..a5e37caa --- /dev/null +++ b/cherab/core/atomic/data/lineshape/zeeman/parametrised/be.json @@ -0,0 +1,25 @@ +{ + "3": { + "5 -> 4": [ + 0.0060354, + 2.1245, + -0.3190 + ], + "6 -> 5": [ + 0.0202754, + 1.6538, + -0.3192 + ], + "7 -> 5": [ + 0.0078966, + 1.7017, + -0.3348 + ], + "8 -> 6": [ + 0.0205025, + 1.4581, + -0.3450 + ] + }, + "reference": "A. Blom and C. Jupén. Parametrisation of the Zeeman effect for hydrogen-like spectra in high-temperature plasmas. Plasma Phys. Control. Fusion 44 (2002) 1229-1241" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/zeeman/parametrised/c.json b/cherab/core/atomic/data/lineshape/zeeman/parametrised/c.json new file mode 100644 index 00000000..8e6ad067 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/zeeman/parametrised/c.json @@ -0,0 +1,45 @@ +{ + "5": { + "6 -> 5": [ + 0.0040900, + 2.4271, + -0.2818 + ], + "7 -> 6": [ + 0.0110398, + 1.9785, + -0.2816 + ], + "8 -> 6": [ + 0.0040747, + 2.1776, + -0.3035 + ], + "8 -> 7": [ + 0.0261405, + 1.6689, + -0.2815 + ], + "9 -> 7": [ + 0.0092096, + 1.8495, + -0.3049 + ], + "10 -> 8": [ + 0.0189020, + 1.6191, + -0.3078 + ], + "11 -> 8": [ + 0.0110428, + 1.6600, + -0.3162 + ], + "10 -> 9": [ + 0.0359009, + 1.4464, + -0.3104 + ] + }, + "reference": "A. Blom and C. Jupén. Parametrisation of the Zeeman effect for hydrogen-like spectra in high-temperature plasmas. Plasma Phys. Control. Fusion 44 (2002) 1229-1241" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/zeeman/parametrised/d.json b/cherab/core/atomic/data/lineshape/zeeman/parametrised/d.json new file mode 100644 index 00000000..3b7e8150 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/zeeman/parametrised/d.json @@ -0,0 +1,15 @@ +{ + "0": { + "3 -> 2": [ + 0.0402068, + 0.4384, + -0.5015 + ], + "4 -> 2": [ + 0.0220610, + 0.3702, + -0.5132 + ] + }, + "reference": "A. Blom and C. Jupén. Parametrisation of the Zeeman effect for hydrogen-like spectra in high-temperature plasmas. Plasma Phys. Control. Fusion 44 (2002) 1229-1241" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/zeeman/parametrised/h.json b/cherab/core/atomic/data/lineshape/zeeman/parametrised/h.json new file mode 100644 index 00000000..13745495 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/zeeman/parametrised/h.json @@ -0,0 +1,15 @@ +{ + "0": { + "3 -> 2": [ + 0.0402267, + 0.3415, + -0.5247 + ], + "4 -> 2": [ + 0.0220724, + 0.2837, + -0.5346 + ] + }, + "reference": "A. Blom and C. Jupén. Parametrisation of the Zeeman effect for hydrogen-like spectra in high-temperature plasmas. Plasma Phys. Control. Fusion 44 (2002) 1229-1241" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/zeeman/parametrised/he.json b/cherab/core/atomic/data/lineshape/zeeman/parametrised/he.json new file mode 100644 index 00000000..2ebd0625 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/zeeman/parametrised/he.json @@ -0,0 +1,25 @@ +{ + "1": { + "4 -> 3": [ + 0.0205206, + 1.6118, + -0.4838 + ], + "5 -> 3": [ + 0.0095879, + 1.4294, + -0.4975 + ], + "6 -> 4": [ + 0.0401955, + 1.0058, + -0.4918 + ], + "7 -> 4": [ + 0.0273521, + 0.9563, + -0.4981 + ] + }, + "reference": "A. Blom and C. Jupén. Parametrisation of the Zeeman effect for hydrogen-like spectra in high-temperature plasmas. Plasma Phys. Control. Fusion 44 (2002) 1229-1241" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/zeeman/parametrised/he3.json b/cherab/core/atomic/data/lineshape/zeeman/parametrised/he3.json new file mode 100644 index 00000000..04729f74 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/zeeman/parametrised/he3.json @@ -0,0 +1,25 @@ +{ + "1": { + "4 -> 3": [ + 0.0205200, + 1.4418, + -0.4892 + ], + "5 -> 3": [ + 0.0095879, + 1.2576, + -0.5001 + ], + "6 -> 4": [ + 0.0401980, + 0.8976, + -0.4971 + ], + "7 -> 4": [ + 0.0273538, + 0.8529, + -0.5039 + ] + }, + "reference": "A. Blom and C. Jupén. Parametrisation of the Zeeman effect for hydrogen-like spectra in high-temperature plasmas. Plasma Phys. Control. Fusion 44 (2002) 1229-1241" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/zeeman/parametrised/n.json b/cherab/core/atomic/data/lineshape/zeeman/parametrised/n.json new file mode 100644 index 00000000..d086f4a0 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/zeeman/parametrised/n.json @@ -0,0 +1,30 @@ +{ + "6": { + "7 -> 6": [ + 0.0060010, + 2.4789, + -0.2817 + ], + "8 -> 7": [ + 0.0141271, + 2.0249, + -0.2762 + ], + "9 -> 8": [ + 0.0300127, + 1.7415, + -0.2753 + ], + "10 -> 8": [ + 0.0102089, + 1.9464, + -0.2975 + ], + "11 -> 9": [ + 0.0193799, + 1.7133, + -0.2973 + ] + }, + "reference": "A. Blom and C. Jupén. Parametrisation of the Zeeman effect for hydrogen-like spectra in high-temperature plasmas. Plasma Phys. Control. Fusion 44 (2002) 1229-1241" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/zeeman/parametrised/ne.json b/cherab/core/atomic/data/lineshape/zeeman/parametrised/ne.json new file mode 100644 index 00000000..07e2c41b --- /dev/null +++ b/cherab/core/atomic/data/lineshape/zeeman/parametrised/ne.json @@ -0,0 +1,25 @@ +{ + "9": { + "9 -> 8": [ + 0.0072488, + 2.8838, + -0.2758 + ], + "10 -> 9": [ + 0.0141002, + 2.4755, + -0.2718 + ], + "11 -> 9": [ + 0.0046673, + 2.8410, + -0.2917 + ], + "11 -> 10": [ + 0.0257292, + 2.1890, + -0.2715 + ] + }, + "reference": "A. Blom and C. Jupén. Parametrisation of the Zeeman effect for hydrogen-like spectra in high-temperature plasmas. Plasma Phys. Control. Fusion 44 (2002) 1229-1241" +} \ No newline at end of file diff --git a/cherab/core/atomic/data/lineshape/zeeman/parametrised/o.json b/cherab/core/atomic/data/lineshape/zeeman/parametrised/o.json new file mode 100644 index 00000000..3f3663d7 --- /dev/null +++ b/cherab/core/atomic/data/lineshape/zeeman/parametrised/o.json @@ -0,0 +1,30 @@ +{ + "7": { + "8 -> 7": [ + 0.0083081, + 2.4263, + -0.2747 + ], + "9 -> 8": [ + 0.0176049, + 2.0652, + -0.2721 + ], + "10 -> 8": [ + 0.0059933, + 2.3445, + -0.2944 + ], + "10 -> 9": [ + 0.0343805, + 1.8122, + -0.2718 + ], + "11 -> 9": [ + 0.0113640, + 2.0268, + -0.2911 + ] + }, + "reference": "A. Blom and C. Jupén. Parametrisation of the Zeeman effect for hydrogen-like spectra in high-temperature plasmas. Plasma Phys. Control. Fusion 44 (2002) 1229-1241" +} \ No newline at end of file diff --git a/cherab/core/atomic/interface.pxd b/cherab/core/atomic/interface.pxd index f6b9b914..1dacfb57 100644 --- a/cherab/core/atomic/interface.pxd +++ b/cherab/core/atomic/interface.pxd @@ -59,5 +59,8 @@ cdef class AtomicData: cpdef ZeemanStructure zeeman_structure(self, Line line, object b_field=*) - cpdef FreeFreeGauntFactor free_free_gaunt_factor(self) + cpdef tuple zeeman_triplet_parameters(self, Line line) + + cpdef tuple stark_model_coefficients(self, Line line) + cpdef FreeFreeGauntFactor free_free_gaunt_factor(self) diff --git a/cherab/core/atomic/interface.pyx b/cherab/core/atomic/interface.pyx index 91417ef4..52c396d6 100644 --- a/cherab/core/atomic/interface.pyx +++ b/cherab/core/atomic/interface.pyx @@ -17,6 +17,8 @@ # under the Licence. from .gaunt import MaxwellianFreeFreeGauntFactor +import json +from os import path cdef class AtomicData: @@ -99,6 +101,10 @@ cdef class AtomicData: raise NotImplementedError("The recombination() virtual method is not implemented for this atomic data source.") cpdef ThermalCXPEC thermal_cx_pec(self, Element donor_ion, int donor_charge, Element receiver_ion, int receiver_charge, tuple transition): + """ + Thermal charge exchange photon emission coefficient for given donor and receiver species in W.m^3. + """ + raise NotImplementedError("The thermal_cx_pec() virtual method is not implemented for this atomic data source.") cpdef TotalRadiatedPower total_radiated_power(self, Element element): @@ -145,6 +151,50 @@ cdef class AtomicData: raise NotImplementedError("The zeeman_structure() virtual method is not implemented for this atomic data source.") + cpdef tuple zeeman_triplet_parameters(self, Line line): + """ + Returns Zeeman truplet parameters. See Table 1 in A. Blom and C. Jupén. + "Parametrisation of the Zeeman effect for hydrogen-like spectra in + high-temperature plasmas", Plasma Phys. Control. Fusion 44 (2002) `1229-1241 + `_. + """ + + symbol = line.element.symbol.lower() + upper, lower = line.transition + encoded_transition = '{} -> {}'.format(str(upper).lower(), str(lower).lower()) + + try: + with open(path.join(path.dirname(__file__), "data/lineshape/zeeman/parametrised/{}.json".format(symbol))) as f: + data = json.load(f) + coefficients = data[str(line.charge)][encoded_transition] + except (FileNotFoundError, KeyError): + raise RuntimeError('Requested Zeeman triplet parameters (element={}, charge={}, transition={})' + ' are not available.'.format(line.element.symbol, line.charge, line.transition)) + + return tuple(coefficients) + + cpdef tuple stark_model_coefficients(self, Line line): + """ + Returns Stark model coefficients. See Table 1 in B. Lomanowski, et al. + "Inferring divertor plasma properties from hydrogen Balmer + and Paschen series spectroscopy in JET-ILW." Nuclear Fusion 55.12 (2015) + `123028 `_. + """ + + symbol = line.element.symbol.lower() + upper, lower = line.transition + encoded_transition = '{} -> {}'.format(str(upper).lower(), str(lower).lower()) + + try: + with open(path.join(path.dirname(__file__), "data/lineshape/stark/{}.json".format(symbol))) as f: + data = json.load(f) + coefficients = data[str(line.charge)][encoded_transition] + except (FileNotFoundError, KeyError): + raise RuntimeError('Requested Stark model coefficients (element={}, charge={}, transition={})' + ' are not available.'.format(line.element.symbol, line.charge, line.transition)) + + return tuple(coefficients) + cpdef FreeFreeGauntFactor free_free_gaunt_factor(self): """ Returns the Maxwellian-averaged free-free Gaunt factor interpolated over the data diff --git a/cherab/core/model/beam/__init__.pxd b/cherab/core/model/beam/__init__.pxd index 8a3c6e22..c5a0ebbe 100644 --- a/cherab/core/model/beam/__init__.pxd +++ b/cherab/core/model/beam/__init__.pxd @@ -1 +1,2 @@ from cherab.core.model.beam.charge_exchange cimport BeamCXLine +from cherab.core.model.beam.beam_emission cimport BeamEmissionLine diff --git a/cherab/core/model/beam/beam_emission.pyx b/cherab/core/model/beam/beam_emission.pyx index 8ca55331..e00aaf5d 100644 --- a/cherab/core/model/beam/beam_emission.pyx +++ b/cherab/core/model/beam/beam_emission.pyx @@ -1,8 +1,8 @@ # cython: language_level=3 -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -22,7 +22,8 @@ cimport cython from libc.math cimport sqrt from raysect.core cimport Point3D, Vector3D -from cherab.core cimport Species, Plasma, Beam, Element, BeamEmissionPEC, Spectrum, AtomicData +from raysect.optical cimport Spectrum +from cherab.core cimport Species, Plasma, Beam, Element, BeamEmissionPEC, AtomicData from cherab.core.math.function cimport autowrap_function1d, autowrap_function2d from cherab.core.atomic.elements import Isotope, hydrogen from cherab.core.model.lineshape cimport BeamEmissionMultiplet @@ -211,8 +212,8 @@ cdef class BeamEmissionLine(BeamModel): self._rates_list.append((species, rate)) # instance line shape renderer - self._lineshape = BeamEmissionMultiplet(self._line, self._wavelength, self._beam, self._sigma_to_pi, - self._sigma1_to_sigma0, self._pi2_to_pi3, self._pi4_to_pi3) + self._lineshape = BeamEmissionMultiplet(self._line, self._wavelength, self._beam, self._atomic_data, + self._sigma_to_pi, self._sigma1_to_sigma0, self._pi2_to_pi3, self._pi4_to_pi3) def _change(self): diff --git a/cherab/core/model/beam/charge_exchange.pyx b/cherab/core/model/beam/charge_exchange.pyx index 9eb562a0..4b5c292e 100644 --- a/cherab/core/model/beam/charge_exchange.pyx +++ b/cherab/core/model/beam/charge_exchange.pyx @@ -362,7 +362,7 @@ cdef class BeamCXLine(BeamModel): self._excited_beam_data.append((rate, population_data)) # instance line shape renderer - self._lineshape = self._lineshape_class(self._line, self._wavelength, self._target_species, self._plasma, + self._lineshape = self._lineshape_class(self._line, self._wavelength, self._target_species, self._plasma, self._atomic_data, *self._lineshape_args, **self._lineshape_kwargs) def _change(self): diff --git a/cherab/core/model/lineshape.pxd b/cherab/core/model/lineshape.pxd deleted file mode 100644 index 2591efe3..00000000 --- a/cherab/core/model/lineshape.pxd +++ /dev/null @@ -1,112 +0,0 @@ -# cython: language_level=3 - -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 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 -cimport numpy as np - -from raysect.optical cimport Spectrum, Point3D, Vector3D -from cherab.core cimport Line, Species, Plasma, Beam -from cherab.core.math cimport Function1D, Function2D -from cherab.core.math.integrators cimport Integrator1D -from cherab.core.atomic.zeeman cimport ZeemanStructure - - -cpdef double doppler_shift(double wavelength, Vector3D observation_direction, Vector3D velocity) - -cpdef double thermal_broadening(double wavelength, double temperature, double atomic_weight) - -cpdef Spectrum add_gaussian_line(double radiance, double wavelength, double sigma, Spectrum spectrum) - - -cdef class LineShapeModel: - - cdef: - Line line - double wavelength - Species target_species - Plasma plasma - Integrator1D integrator - - cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum) - - -cdef class GaussianLine(LineShapeModel): - pass - - -cdef class MultipletLineShape(LineShapeModel): - - cdef: - int _number_of_lines - np.ndarray _multiplet - double[:,::1] _multiplet_mv - - -cdef class StarkBroadenedLine(LineShapeModel): - - cdef double _aij, _bij, _cij - - pass - - -cdef class ZeemanLineShapeModel(LineShapeModel): - - cdef double _polarisation - - pass - - -cdef class ZeemanTriplet(ZeemanLineShapeModel): - - pass - - -cdef class ParametrisedZeemanTriplet(ZeemanLineShapeModel): - - cdef double _alpha, _beta, _gamma - - pass - - -cdef class ZeemanMultiplet(ZeemanLineShapeModel): - - cdef ZeemanStructure _zeeman_structure - - pass - - -cdef class BeamLineShapeModel: - - cdef: - - Line line - double wavelength - Beam beam - - cpdef Spectrum add_line(self, double radiance, Point3D beam_point, Point3D plasma_point, - Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum) - - -cdef class BeamEmissionMultiplet(BeamLineShapeModel): - - cdef: - - Function2D _sigma_to_pi - Function1D _sigma1_to_sigma0, _pi2_to_pi3, _pi4_to_pi3 diff --git a/cherab/core/model/lineshape.pyx b/cherab/core/model/lineshape.pyx deleted file mode 100644 index 69278560..00000000 --- a/cherab/core/model/lineshape.pyx +++ /dev/null @@ -1,968 +0,0 @@ -# cython: language_level=3 - -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas -# -# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the -# European Commission - subsequent versions of the EUPL (the "Licence"); -# You may not use this work except in compliance with the Licence. -# You may obtain a copy of the Licence at: -# -# https://joinup.ec.europa.eu/software/page/eupl5 -# -# Unless required by applicable law or agreed to in writing, software distributed -# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR -# CONDITIONS OF ANY KIND, either express or implied. -# -# See the Licence for the specific language governing permissions and limitations -# under the Licence. - -import numpy as np -from scipy.special import hyp2f1 - -cimport numpy as np -from libc.math cimport sqrt, erf, M_SQRT2, floor, ceil, fabs -from raysect.optical.spectrum cimport new_spectrum -from raysect.core.math.function.float cimport Function1D - -from cherab.core cimport Plasma -from cherab.core.atomic.elements import hydrogen, deuterium, tritium, helium, helium3, beryllium, boron, carbon, nitrogen, oxygen, neon -from cherab.core.math.function cimport autowrap_function1d, autowrap_function2d -from cherab.core.math.integrators cimport GaussianQuadrature -from cherab.core.utility.constants cimport ATOMIC_MASS, ELEMENTARY_CHARGE, SPEED_OF_LIGHT - -cimport cython - -# required by numpy c-api -np.import_array() - - -cdef double RECIP_ATOMIC_MASS = 1 / ATOMIC_MASS - - -cdef double evamu_to_ms(double x): - return sqrt(2 * x * ELEMENTARY_CHARGE * RECIP_ATOMIC_MASS) - - -@cython.cdivision(True) -cpdef double doppler_shift(double wavelength, Vector3D observation_direction, Vector3D velocity): - """ - Calculates the Doppler shifted wavelength for a given velocity and observation direction. - - :param wavelength: The wavelength to Doppler shift in nanometers. - :param observation_direction: A Vector defining the direction of observation. - :param velocity: A Vector defining the relative velocity of the emitting source in m/s. - :return: The Doppler shifted wavelength in nanometers. - """ - cdef double projected_velocity - - # flow velocity projected on the direction of observation - observation_direction = observation_direction.normalise() - projected_velocity = velocity.dot(observation_direction) - - return wavelength * (1 + projected_velocity / SPEED_OF_LIGHT) - - -@cython.cdivision(True) -cpdef double thermal_broadening(double wavelength, double temperature, double atomic_weight): - """ - Returns the line width for a gaussian line as a standard deviation. - - :param wavelength: Central wavelength. - :param temperature: Temperature in eV. - :param atomic_weight: Atomic weight in AMU. - :return: Standard deviation of gaussian line. - """ - - # todo: add input sanity checks - return sqrt(temperature * ELEMENTARY_CHARGE / (atomic_weight * ATOMIC_MASS)) * wavelength / SPEED_OF_LIGHT - - -# the number of standard deviations outside the rest wavelength the line is considered to add negligible value (including a margin for safety) -DEF GAUSSIAN_CUTOFF_SIGMA = 10.0 - - -@cython.cdivision(True) -@cython.initializedcheck(False) -@cython.boundscheck(False) -@cython.wraparound(False) -cpdef Spectrum add_gaussian_line(double radiance, double wavelength, double sigma, Spectrum spectrum): - r""" - Adds a Gaussian line to the given spectrum and returns the new spectrum. - - The formula used is based on the following definite integral: - :math:`\frac{1}{\sigma \sqrt{2 \pi}} \int_{\lambda_0}^{\lambda_1} \exp(-\frac{(x-\mu)^2}{2\sigma^2}) dx = \frac{1}{2} \left[ -Erf(\frac{a-\mu}{\sqrt{2}\sigma}) +Erf(\frac{b-\mu}{\sqrt{2}\sigma}) \right]` - - :param float radiance: Intensity of the line in radiance. - :param float wavelength: central wavelength of the line in nm. - :param float sigma: width of the line in nm. - :param Spectrum spectrum: the current spectrum to which the gaussian line is added. - :return: - """ - - cdef double temp - cdef double cutoff_lower_wavelength, cutoff_upper_wavelength - cdef double lower_wavelength, upper_wavelength - cdef double lower_integral, upper_integral - cdef int start, end, i - - if sigma <= 0: - return spectrum - - # calculate and check end of limits - cutoff_lower_wavelength = wavelength - GAUSSIAN_CUTOFF_SIGMA * sigma - if spectrum.max_wavelength < cutoff_lower_wavelength: - return spectrum - - cutoff_upper_wavelength = wavelength + GAUSSIAN_CUTOFF_SIGMA * sigma - if spectrum.min_wavelength > cutoff_upper_wavelength: - return spectrum - - # locate range of bins where there is significant contribution from the gaussian (plus a health margin) - start = max(0, floor((cutoff_lower_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) - end = min(spectrum.bins, ceil((cutoff_upper_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) - - # add line to spectrum - temp = 1 / (M_SQRT2 * sigma) - lower_wavelength = spectrum.min_wavelength + start * spectrum.delta_wavelength - lower_integral = erf((lower_wavelength - wavelength) * temp) - for i in range(start, end): - - upper_wavelength = spectrum.min_wavelength + spectrum.delta_wavelength * (i + 1) - upper_integral = erf((upper_wavelength - wavelength) * temp) - - spectrum.samples_mv[i] += radiance * 0.5 * (upper_integral - lower_integral) / spectrum.delta_wavelength - - lower_wavelength = upper_wavelength - lower_integral = upper_integral - - return spectrum - - -cdef class LineShapeModel: - """ - A base class for building line shapes. - - :param Line line: The emission line object for this line shape. - :param float wavelength: The rest wavelength for this emission line. - :param Species target_species: The target plasma species that is emitting. - :param Plasma plasma: The emitting plasma object. - :param Integrator1D integrator: Integrator1D instance to integrate the line shape - over the spectral bin. Default is None. - """ - - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, Integrator1D integrator=None): - - self.line = line - self.wavelength = wavelength - self.target_species = target_species - self.plasma = plasma - self.integrator = integrator - - cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): - raise NotImplementedError('Child lineshape class must implement this method.') - - -cdef class GaussianLine(LineShapeModel): - """ - Produces Gaussian line shape. - - :param Line line: The emission line object for this line shape. - :param float wavelength: The rest wavelength for this emission line. - :param Species target_species: The target plasma species that is emitting. - :param Plasma plasma: The emitting plasma object. - - .. code-block:: pycon - - >>> from cherab.core.atomic import Line, deuterium - >>> from cherab.core.model import ExcitationLine, GaussianLine - >>> - >>> # Adding Gaussian line to the plasma model. - >>> d_alpha = Line(deuterium, 0, (3, 2)) - >>> excit = ExcitationLine(d_alpha, lineshape=GaussianLine) - >>> plasma.models.add(excit) - """ - - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma): - - super().__init__(line, wavelength, target_species, plasma) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - @cython.cdivision(True) - cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): - - cdef double ts, sigma, shifted_wavelength - cdef Vector3D ion_velocity - - ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) - if ts <= 0.0: - return spectrum - - ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) - - # calculate emission line central wavelength, doppler shifted along observation direction - shifted_wavelength = doppler_shift(self.wavelength, direction, ion_velocity) - - # calculate the line width - sigma = thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) - - return add_gaussian_line(radiance, shifted_wavelength, sigma, spectrum) - - -DEF MULTIPLET_WAVELENGTH = 0 -DEF MULTIPLET_RATIO = 1 - - -cdef class MultipletLineShape(LineShapeModel): - """ - Produces Multiplet line shapes. - - The lineshape radiance is calculated from a base PEC rate that is unresolved. This - radiance is then divided over a number of components as specified in the multiplet - argument. The multiplet components are specified with an Nx2 array where N is the - number of components in the multiplet. The first axis of the array contains the - wavelengths of each component, the second contains the line ratio for each component. - The component line ratios must sum to one. For example: - - :param Line line: The emission line object for the base rate radiance calculation. - :param float wavelength: The rest wavelength of the base emission line. - :param Species target_species: The target plasma species that is emitting. - :param Plasma plasma: The emitting plasma object. - :param multiplet: An Nx2 array that specifies the multiplet wavelengths and line ratios. - - .. code-block:: pycon - - >>> from cherab.core.atomic import Line, nitrogen - >>> from cherab.core.model import ExcitationLine, MultipletLineShape - >>> - >>> # multiplet specification in Nx2 array - >>> multiplet = [[403.509, 404.132, 404.354, 404.479, 405.692], [0.205, 0.562, 0.175, 0.029, 0.029]] - >>> - >>> # Adding the multiplet to the plasma model. - >>> nitrogen_II_404 = Line(nitrogen, 1, ("2s2 2p1 4f1 3G13.0", "2s2 2p1 3d1 3F10.0")) - >>> excit = ExcitationLine(nitrogen_II_404, lineshape=MultipletLineShape, lineshape_args=[multiplet]) - >>> plasma.models.add(excit) - """ - - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, - object multiplet): - - super().__init__(line, wavelength, target_species, plasma) - - multiplet = np.array(multiplet, dtype=np.float64) - - if not (len(multiplet.shape) == 2 and multiplet.shape[0] == 2): - raise ValueError("The multiplet specification must be an array of shape (Nx2).") - - if not multiplet[1,:].sum() == 1.0: - raise ValueError("The multiplet line ratios should sum to one.") - - self._number_of_lines = multiplet.shape[1] - self._multiplet = multiplet - self._multiplet_mv = self._multiplet - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): - - cdef double ts, sigma, shifted_wavelength, component_wavelength, component_radiance - cdef Vector3D ion_velocity - - ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) - if ts <= 0.0: - return spectrum - - ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) - - # calculate the line width - sigma = thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) - - for i in range(self._number_of_lines): - - component_wavelength = self._multiplet_mv[MULTIPLET_WAVELENGTH, i] - component_radiance = radiance * self._multiplet_mv[MULTIPLET_RATIO, i] - - # calculate emission line central wavelength, doppler shifted along observation direction - shifted_wavelength = doppler_shift(component_wavelength, direction, ion_velocity) - - spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) - - return spectrum - - -DEF LORENZIAN_CUTOFF_GAMMA = 50.0 - - -cdef class StarkFunction(Function1D): - """ - Normalised Stark function for the StarkBroadenedLine line shape. - """ - - cdef double _a, _x0, _norm - - STARK_NORM_COEFFICIENT = 4 * LORENZIAN_CUTOFF_GAMMA * hyp2f1(0.4, 1, 1.4, -(2 * LORENZIAN_CUTOFF_GAMMA)**2.5) - - def __init__(self, double wavelength, double lambda_1_2): - - if wavelength <= 0: - raise ValueError("Argument 'wavelength' must be positive.") - - if lambda_1_2 <= 0: - raise ValueError("Argument 'lambda_1_2' must be positive.") - - self._x0 = wavelength - self._a = (0.5 * lambda_1_2)**2.5 - # normalise, so the integral over x is equal to 1 in the limits - # (_x0 - LORENZIAN_CUTOFF_GAMMA * lambda_1_2, _x0 + LORENZIAN_CUTOFF_GAMMA * lambda_1_2) - self._norm = (0.5 * lambda_1_2)**1.5 / self.STARK_NORM_COEFFICIENT - - @cython.cdivision(True) - cdef double evaluate(self, double x) except? -1e999: - - return self._norm / ((fabs(x - self._x0))**2.5 + self._a) - - -cdef class StarkBroadenedLine(LineShapeModel): - """ - Parametrised Stark broadened line shape based on the Model Microfield Method (MMM). - Contains embedded atomic data in the form of fits to MMM. - Only Balmer and Paschen series are supported by default. - See B. Lomanowski, et al. "Inferring divertor plasma properties from hydrogen Balmer - and Paschen series spectroscopy in JET-ILW." Nuclear Fusion 55.12 (2015) - `123028 `_. - - Call `show_supported_transitions()` to see the list of supported transitions and - default model coefficients. - - :param Line line: The emission line object for this line shape. - :param float wavelength: The rest wavelength for this emission line. - :param Species target_species: The target plasma species that is emitting. - :param Plasma plasma: The emitting plasma object. - :param dict stark_model_coefficients: Alternative model coefficients in the form - {line_ij: (c_ij, a_ij, b_ij), ...}. - If None, the default model parameters will be used. - :param Integrator1D integrator: Integrator1D instance to integrate the line shape - over the spectral bin. Default is `GaussianQuadrature()`. - - """ - - STARK_MODEL_COEFFICIENTS_DEFAULT = { - Line(hydrogen, 0, (3, 2)): (3.71e-18, 0.7665, 0.064), - Line(hydrogen, 0, (4, 2)): (8.425e-18, 0.7803, 0.050), - Line(hydrogen, 0, (5, 2)): (1.31e-15, 0.6796, 0.030), - Line(hydrogen, 0, (6, 2)): (3.954e-16, 0.7149, 0.028), - Line(hydrogen, 0, (7, 2)): (6.258e-16, 0.712, 0.029), - Line(hydrogen, 0, (8, 2)): (7.378e-16, 0.7159, 0.032), - Line(hydrogen, 0, (9, 2)): (8.947e-16, 0.7177, 0.033), - Line(hydrogen, 0, (4, 3)): (1.330e-16, 0.7449, 0.045), - Line(hydrogen, 0, (5, 3)): (6.64e-16, 0.7356, 0.044), - Line(hydrogen, 0, (6, 3)): (2.481e-15, 0.7118, 0.016), - Line(hydrogen, 0, (7, 3)): (3.270e-15, 0.7137, 0.029), - Line(hydrogen, 0, (8, 3)): (4.343e-15, 0.7133, 0.032), - Line(hydrogen, 0, (9, 3)): (5.588e-15, 0.7165, 0.033), - Line(deuterium, 0, (3, 2)): (3.71e-18, 0.7665, 0.064), - Line(deuterium, 0, (4, 2)): (8.425e-18, 0.7803, 0.050), - Line(deuterium, 0, (5, 2)): (1.31e-15, 0.6796, 0.030), - Line(deuterium, 0, (6, 2)): (3.954e-16, 0.7149, 0.028), - Line(deuterium, 0, (7, 2)): (6.258e-16, 0.712, 0.029), - Line(deuterium, 0, (8, 2)): (7.378e-16, 0.7159, 0.032), - Line(deuterium, 0, (9, 2)): (8.947e-16, 0.7177, 0.033), - Line(deuterium, 0, (4, 3)): (1.330e-16, 0.7449, 0.045), - Line(deuterium, 0, (5, 3)): (6.64e-16, 0.7356, 0.044), - Line(deuterium, 0, (6, 3)): (2.481e-15, 0.7118, 0.016), - Line(deuterium, 0, (7, 3)): (3.270e-15, 0.7137, 0.029), - Line(deuterium, 0, (8, 3)): (4.343e-15, 0.7133, 0.032), - Line(deuterium, 0, (9, 3)): (5.588e-15, 0.7165, 0.033), - Line(tritium, 0, (3, 2)): (3.71e-18, 0.7665, 0.064), - Line(tritium, 0, (4, 2)): (8.425e-18, 0.7803, 0.050), - Line(tritium, 0, (5, 2)): (1.31e-15, 0.6796, 0.030), - Line(tritium, 0, (6, 2)): (3.954e-16, 0.7149, 0.028), - Line(tritium, 0, (7, 2)): (6.258e-16, 0.712, 0.029), - Line(tritium, 0, (8, 2)): (7.378e-16, 0.7159, 0.032), - Line(tritium, 0, (9, 2)): (8.947e-16, 0.7177, 0.033), - Line(tritium, 0, (4, 3)): (1.330e-16, 0.7449, 0.045), - Line(tritium, 0, (5, 3)): (6.64e-16, 0.7356, 0.044), - Line(tritium, 0, (6, 3)): (2.481e-15, 0.7118, 0.016), - Line(tritium, 0, (7, 3)): (3.270e-15, 0.7137, 0.029), - Line(tritium, 0, (8, 3)): (4.343e-15, 0.7133, 0.032), - Line(tritium, 0, (9, 3)): (5.588e-15, 0.7165, 0.033) - } - - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, - dict stark_model_coefficients=None, integrator=GaussianQuadrature()): - - stark_model_coefficients = stark_model_coefficients or self.STARK_MODEL_COEFFICIENTS_DEFAULT - - try: - # Fitted Stark Constants - cij, aij, bij = stark_model_coefficients[line] - if cij <= 0: - raise ValueError('Coefficient c_ij must be positive.') - if aij <= 0: - raise ValueError('Coefficient a_ij must be positive.') - if bij <= 0: - raise ValueError('Coefficient b_ij must be positive.') - self._aij = aij - self._bij = bij - self._cij = cij - except IndexError: - raise ValueError('Stark broadening coefficients for {} is not currently available.'.format(line)) - - super().__init__(line, wavelength, target_species, plasma, integrator) - - def show_supported_transitions(self): - """ Prints all supported transitions.""" - for line, coeff in self.STARK_MODEL_COEFFICIENTS_DEFAULT.items(): - print('{}: c_ij={}, a_ij={}, b_ij={}'.format(line, coeff[0], coeff[1], coeff[2])) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - @cython.cdivision(True) - cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): - - cdef: - double ne, te, lambda_1_2, lambda_5_2, wvl - double cutoff_lower_wavelength, cutoff_upper_wavelength - double lower_wavelength, upper_wavelength - double bin_integral - int start, end, i - Spectrum raw_lineshape - - ne = self.plasma.get_electron_distribution().density(point.x, point.y, point.z) - if ne <= 0.0: - return spectrum - - te = self.plasma.get_electron_distribution().effective_temperature(point.x, point.y, point.z) - if te <= 0.0: - return spectrum - - lambda_1_2 = self._cij * ne**self._aij / (te**self._bij) - - self.integrator.function = StarkFunction(self.wavelength, lambda_1_2) - - # calculate and check end of limits - cutoff_lower_wavelength = self.wavelength - LORENZIAN_CUTOFF_GAMMA * lambda_1_2 - if spectrum.max_wavelength < cutoff_lower_wavelength: - return spectrum - - cutoff_upper_wavelength = self.wavelength + LORENZIAN_CUTOFF_GAMMA * lambda_1_2 - if spectrum.min_wavelength > cutoff_upper_wavelength: - return spectrum - - # locate range of bins where there is significant contribution from the gaussian (plus a health margin) - start = max(0, floor((cutoff_lower_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) - end = min(spectrum.bins, ceil((cutoff_upper_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) - - # add line to spectrum - lower_wavelength = spectrum.min_wavelength + start * spectrum.delta_wavelength - - for i in range(start, end): - upper_wavelength = spectrum.min_wavelength + spectrum.delta_wavelength * (i + 1) - - bin_integral = self.integrator.evaluate(lower_wavelength, upper_wavelength) - spectrum.samples_mv[i] += radiance * bin_integral / spectrum.delta_wavelength - - lower_wavelength = upper_wavelength - - return spectrum - - -DEF BOHR_MAGNETON = 5.78838180123e-5 # in eV/T -DEF HC_EV_NM = 1239.8419738620933 # (Planck constant in eV s) x (speed of light in nm/s) - -DEF PI_POLARISATION = 0 -DEF SIGMA_POLARISATION = 1 -DEF SIGMA_PLUS_POLARISATION = 1 -DEF SIGMA_MINUS_POLARISATION = -1 -DEF NO_POLARISATION = 2 - - -cdef class ZeemanLineShapeModel(LineShapeModel): - r""" - A base class for building Zeeman line shapes. - - :param Line line: The emission line object for this line shape. - :param float wavelength: The rest wavelength for this emission line. - :param Species target_species: The target plasma species that is emitting. - :param Plasma plasma: The emitting plasma object. - :param str polarisation: Leaves only :math:`\pi`-/:math:`\sigma`-polarised components: - "pi" - leave only :math:`\pi`-polarised components, - "sigma" - leave only :math:`\sigma`-polarised components, - "no" - leave all components (default). - """ - - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, polarisation): - super().__init__(line, wavelength, target_species, plasma) - - self.polarisation = polarisation - - @property - def polarisation(self): - if self._polarisation == PI_POLARISATION: - return 'pi' - if self._polarisation == SIGMA_POLARISATION: - return 'sigma' - if self._polarisation == NO_POLARISATION: - return 'no' - - @polarisation.setter - def polarisation(self, value): - if value.lower() == 'pi': - self._polarisation = PI_POLARISATION - elif value.lower() == 'sigma': - self._polarisation = SIGMA_POLARISATION - elif value.lower() == 'no': - self._polarisation = NO_POLARISATION - else: - raise ValueError('Select between "pi", "sigma" or "no", {} is unsupported.'.format(value)) - - -cdef class ZeemanTriplet(ZeemanLineShapeModel): - r""" - Simple Doppler-Zeeman triplet (Paschen-Back effect). - - :param Line line: The emission line object for this line shape. - :param float wavelength: The rest wavelength for this emission line. - :param Species target_species: The target plasma species that is emitting. - :param Plasma plasma: The emitting plasma object. - :param str polarisation: Leaves only :math:`\pi`-/:math:`\sigma`-polarised components: - "pi" - leave central component, - "sigma" - leave side components, - "no" - all components (default). - """ - - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, polarisation='no'): - - super().__init__(line, wavelength, target_species, plasma, polarisation) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - @cython.cdivision(True) - cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): - - cdef double ts, sigma, shifted_wavelength, photon_energy, b_magn, component_radiance, cos_sqr, sin_sqr - cdef Vector3D ion_velocity, b_field - - ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) - if ts <= 0.0: - return spectrum - - ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) - - # calculate emission line central wavelength, doppler shifted along observation direction - shifted_wavelength = doppler_shift(self.wavelength, direction, ion_velocity) - - # calculate the line width - sigma = thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) - - # obtain magnetic field - b_field = self.plasma.get_b_field().evaluate(point.x, point.y, point.z) - b_magn = b_field.get_length() - - if b_magn == 0: - # no splitting if magnetic field strength is zero - if self._polarisation == NO_POLARISATION: - return add_gaussian_line(radiance, shifted_wavelength, sigma, spectrum) - - return add_gaussian_line(0.5 * radiance, shifted_wavelength, sigma, spectrum) - - # coefficients for intensities parallel and perpendicular to magnetic field - cos_sqr = (b_field.dot(direction.normalise()) / b_magn)**2 - sin_sqr = 1. - cos_sqr - - # adding pi component of the Zeeman triplet in case of NO_POLARISATION or PI_POLARISATION - if self._polarisation != SIGMA_POLARISATION: - component_radiance = 0.5 * sin_sqr * radiance - spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) - - # adding sigma +/- components of the Zeeman triplet in case of NO_POLARISATION or SIGMA_POLARISATION - if self._polarisation != PI_POLARISATION: - component_radiance = (0.25 * sin_sqr + 0.5 * cos_sqr) * radiance - - photon_energy = HC_EV_NM / self.wavelength - - shifted_wavelength = doppler_shift(HC_EV_NM / (photon_energy - BOHR_MAGNETON * b_magn), direction, ion_velocity) - spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) - - shifted_wavelength = doppler_shift(HC_EV_NM / (photon_energy + BOHR_MAGNETON * b_magn), direction, ion_velocity) - spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) - - return spectrum - - -cdef class ParametrisedZeemanTriplet(ZeemanLineShapeModel): - r""" - Parametrised Doppler-Zeeman triplet. It takes into account additional broadening due to - the line's fine structure without resolving the individual components of the fine - structure. The model is described with three parameters: :math:`\alpha`, - :math:`\beta` and :math:`\gamma`. - - The distance between :math:`\sigma^+` and :math:`\sigma^-` peaks: - :math:`\Delta \lambda_{\sigma} = \alpha B`, - where `B` is the magnetic field strength. - The ratio between Zeeman and thermal broadening line widths: - :math:`\frac{W_{Zeeman}}{W_{Doppler}} = \beta T^{\gamma}`, - where `T` is the species temperature in eV. - - Call `show_supported_transitions()` to see the list of supported transitions and - default parameters of the model. - - For details see A. Blom and C. Jupén, Parametrisation of the Zeeman effect - for hydrogen-like spectra in high-temperature plasmas, - Plasma Phys. Control. Fusion 44 (2002) `1229-1241 - `_. - - :param Line line: The emission line object for this line shape. - :param float wavelength: The rest wavelength for this emission line. - :param Species target_species: The target plasma species that is emitting. - :param Plasma plasma: The emitting plasma object. - :param dict line_parameters: Alternative parameters of the model in the form - {line_i: (alpha_i, beta_i, gamma_i), ...}. - If None, the default model parameters will be used. - :param str polarisation: Leaves only :math:`\pi`-/:math:`\sigma`-polarised components: - "pi" - leave central component, - "sigma" - leave side components, - "no" - all components (default). - """ - - LINE_PARAMETERS_DEFAULT = { # alpha, beta, gamma parameters for selected lines - Line(hydrogen, 0, (3, 2)): (0.0402267, 0.3415, -0.5247), - Line(hydrogen, 0, (4, 2)): (0.0220724, 0.2837, -0.5346), - Line(deuterium, 0, (3, 2)): (0.0402068, 0.4384, -0.5015), - Line(deuterium, 0, (4, 2)): (0.0220610, 0.3702, -0.5132), - Line(helium3, 1, (4, 3)): (0.0205200, 1.4418, -0.4892), - Line(helium3, 1, (5, 3)): (0.0095879, 1.2576, -0.5001), - Line(helium3, 1, (6, 4)): (0.0401980, 0.8976, -0.4971), - Line(helium3, 1, (7, 4)): (0.0273538, 0.8529, -0.5039), - Line(helium, 1, (4, 3)): (0.0205206, 1.6118, -0.4838), - Line(helium, 1, (5, 3)): (0.0095879, 1.4294, -0.4975), - Line(helium, 1, (6, 4)): (0.0401955, 1.0058, -0.4918), - Line(helium, 1, (7, 4)): (0.0273521, 0.9563, -0.4981), - Line(beryllium, 3, (5, 4)): (0.0060354, 2.1245, -0.3190), - Line(beryllium, 3, (6, 5)): (0.0202754, 1.6538, -0.3192), - Line(beryllium, 3, (7, 5)): (0.0078966, 1.7017, -0.3348), - Line(beryllium, 3, (8, 6)): (0.0205025, 1.4581, -0.3450), - Line(boron, 4, (6, 5)): (0.0083423, 2.0519, -0.2960), - Line(boron, 4, (7, 6)): (0.0228379, 1.6546, -0.2941), - Line(boron, 4, (8, 6)): (0.0084065, 1.8041, -0.3177), - Line(boron, 4, (8, 7)): (0.0541883, 1.4128, -0.2966), - Line(boron, 4, (9, 7)): (0.0190781, 1.5440, -0.3211), - Line(boron, 4, (10, 8)): (0.0391914, 1.3569, -0.3252), - Line(carbon, 5, (6, 5)): (0.0040900, 2.4271, -0.2818), - Line(carbon, 5, (7, 6)): (0.0110398, 1.9785, -0.2816), - Line(carbon, 5, (8, 6)): (0.0040747, 2.1776, -0.3035), - Line(carbon, 5, (8, 7)): (0.0261405, 1.6689, -0.2815), - Line(carbon, 5, (9, 7)): (0.0092096, 1.8495, -0.3049), - Line(carbon, 5, (10, 8)): (0.0189020, 1.6191, -0.3078), - Line(carbon, 5, (11, 8)): (0.0110428, 1.6600, -0.3162), - Line(carbon, 5, (10, 9)): (0.0359009, 1.4464, -0.3104), - Line(nitrogen, 6, (7, 6)): (0.0060010, 2.4789, -0.2817), - Line(nitrogen, 6, (8, 7)): (0.0141271, 2.0249, -0.2762), - Line(nitrogen, 6, (9, 8)): (0.0300127, 1.7415, -0.2753), - Line(nitrogen, 6, (10, 8)): (0.0102089, 1.9464, -0.2975), - Line(nitrogen, 6, (11, 9)): (0.0193799, 1.7133, -0.2973), - Line(oxygen, 7, (8, 7)): (0.0083081, 2.4263, -0.2747), - Line(oxygen, 7, (9, 8)): (0.0176049, 2.0652, -0.2721), - Line(oxygen, 7, (10, 8)): (0.0059933, 2.3445, -0.2944), - Line(oxygen, 7, (10, 9)): (0.0343805, 1.8122, -0.2718), - Line(oxygen, 7, (11, 9)): (0.0113640, 2.0268, -0.2911), - Line(neon, 9, (9, 8)): (0.0072488, 2.8838, -0.2758), - Line(neon, 9, (10, 9)): (0.0141002, 2.4755, -0.2718), - Line(neon, 9, (11, 9)): (0.0046673, 2.8410, -0.2917), - Line(neon, 9, (11, 10)): (0.0257292, 2.1890, -0.2715) - } - - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, dict line_parameters=None, polarisation='no'): - - super().__init__(line, wavelength, target_species, plasma, polarisation) - - line_parameters = line_parameters or self.LINE_PARAMETERS_DEFAULT - - try: - alpha, beta, gamma = line_parameters[self.line] - if alpha <= 0: - raise ValueError('Parameter alpha must be positive.') - if beta < 0: - raise ValueError('Parameter beta must be non-negative.') - self._alpha = alpha - self._beta = beta - self._gamma = gamma - - except KeyError: - raise ValueError('Data for {} is not available.'.format(self.line)) - - def show_supported_transitions(self): - """ Prints all supported transitions.""" - for line, param in self.LINE_PARAMETERS_DEFAULT.items(): - print('{}: alpha={}, beta={}, gamma={}'.format(line, param[0], param[1], param[2])) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - @cython.cdivision(True) - cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): - - cdef double ts, sigma, shifted_wavelength, b_magn, component_radiance, cos_sqr, sin_sqr - cdef Vector3D ion_velocity, b_field - - ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) - if ts <= 0.0: - return spectrum - - ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) - - # calculate emission line central wavelength, doppler shifted along observation direction - shifted_wavelength = doppler_shift(self.wavelength, direction, ion_velocity) - - # calculate the line width - sigma = thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) - - # fine structure broadening correction - sigma *= sqrt(1. + self._beta * self._beta * ts**(2. * self._gamma)) - - # obtain magnetic field - b_field = self.plasma.get_b_field().evaluate(point.x, point.y, point.z) - b_magn = b_field.get_length() - - if b_magn == 0: - # no splitting if magnetic filed strength is zero - if self._polarisation == NO_POLARISATION: - return add_gaussian_line(radiance, shifted_wavelength, sigma, spectrum) - - return add_gaussian_line(0.5 * radiance, shifted_wavelength, sigma, spectrum) - - # coefficients for intensities parallel and perpendicular to magnetic field - cos_sqr = (b_field.dot(direction.normalise()) / b_magn)**2 - sin_sqr = 1. - cos_sqr - - # adding pi component of the Zeeman triplet in case of NO_POLARISATION or PI_POLARISATION - if self._polarisation != SIGMA_POLARISATION: - component_radiance = 0.5 * sin_sqr * radiance - spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) - - # adding sigma +/- components of the Zeeman triplet in case of NO_POLARISATION or SIGMA_POLARISATION - if self._polarisation != PI_POLARISATION: - component_radiance = (0.25 * sin_sqr + 0.5 * cos_sqr) * radiance - - shifted_wavelength = doppler_shift(self.wavelength + 0.5 * self._alpha * b_magn, direction, ion_velocity) - spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) - shifted_wavelength = doppler_shift(self.wavelength - 0.5 * self._alpha * b_magn, direction, ion_velocity) - spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) - - return spectrum - - -cdef class ZeemanMultiplet(ZeemanLineShapeModel): - r""" - Doppler-Zeeman Multiplet. - - The lineshape radiance is calculated from a base PEC rate that is unresolved. This - radiance is then divided over a number of components as specified in the ``zeeman_structure`` - argument. The ``zeeman_structure`` specifies wavelengths and ratios of - :math:`\pi`-/:math:`\sigma`-polarised components as functions of the magnetic field strength. - These functions can be obtained using the output of the ADAS603 routines. - - :param Line line: The emission line object for the base rate radiance calculation. - :param float wavelength: The rest wavelength of the base emission line. - :param Species target_species: The target plasma species that is emitting. - :param Plasma plasma: The emitting plasma object. - :param zeeman_structure: A ``ZeemanStructure`` object that provides wavelengths and ratios - of :math:`\pi`-/:math:`\sigma^{+}`-/:math:`\sigma^{-}`-polarised - components for any given magnetic field strength. - :param str polarisation: Leaves only :math:`\pi`-/:math:`\sigma`-polarised components: - "pi" - leave only :math:`\pi`-polarised components, - "sigma" - leave only :math:`\sigma`-polarised components, - "no" - leave all components (default). - - """ - - def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, - ZeemanStructure zeeman_structure, polarisation='no'): - - super().__init__(line, wavelength, target_species, plasma, polarisation) - - self._zeeman_structure = zeeman_structure - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - @cython.cdivision(True) - cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): - - cdef int i - cdef double ts, sigma, shifted_wavelength, component_radiance - cdef Vector3D ion_velocity - cdef double[:, :] multiplet_pi_mv, multiplet_sigma_mv - - ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) - if ts <= 0.0: - return spectrum - - ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) - - # calculate the line width - sigma = thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) - - # obtain magnetic field - b_field = self.plasma.get_b_field().evaluate(point.x, point.y, point.z) - b_magn = b_field.get_length() - - if b_magn == 0: - # no splitting if magnetic filed strength is zero - shifted_wavelength = doppler_shift(self.wavelength, direction, ion_velocity) - if self._polarisation == NO_POLARISATION: - return add_gaussian_line(radiance, shifted_wavelength, sigma, spectrum) - - return add_gaussian_line(0.5 * radiance, shifted_wavelength, sigma, spectrum) - - # coefficients for intensities parallel and perpendicular to magnetic field - cos_sqr = (b_field.dot(direction.normalise()) / b_magn)**2 - sin_sqr = 1. - cos_sqr - - # adding pi components of the Zeeman multiplet in case of NO_POLARISATION or PI_POLARISATION - if self._polarisation != SIGMA_POLARISATION: - component_radiance = 0.5 * sin_sqr * radiance - multiplet_mv = self._zeeman_structure.evaluate(b_magn, PI_POLARISATION) - - for i in range(multiplet_mv.shape[1]): - shifted_wavelength = doppler_shift(multiplet_mv[MULTIPLET_WAVELENGTH, i], direction, ion_velocity) - spectrum = add_gaussian_line(component_radiance * multiplet_mv[MULTIPLET_RATIO, i], shifted_wavelength, sigma, spectrum) - - # adding sigma components of the Zeeman multiplet in case of NO_POLARISATION or SIGMA_POLARISATION - if self._polarisation != PI_POLARISATION: - component_radiance = (0.25 * sin_sqr + 0.5 * cos_sqr) * radiance - - multiplet_mv = self._zeeman_structure.evaluate(b_magn, SIGMA_PLUS_POLARISATION) - - for i in range(multiplet_mv.shape[1]): - shifted_wavelength = doppler_shift(multiplet_mv[MULTIPLET_WAVELENGTH, i], direction, ion_velocity) - spectrum = add_gaussian_line(component_radiance * multiplet_mv[MULTIPLET_RATIO, i], shifted_wavelength, sigma, spectrum) - - multiplet_mv = self._zeeman_structure.evaluate(b_magn, SIGMA_MINUS_POLARISATION) - - for i in range(multiplet_mv.shape[1]): - shifted_wavelength = doppler_shift(multiplet_mv[MULTIPLET_WAVELENGTH, i], direction, ion_velocity) - spectrum = add_gaussian_line(component_radiance * multiplet_mv[MULTIPLET_RATIO, i], shifted_wavelength, sigma, spectrum) - - return spectrum - - -cdef class BeamLineShapeModel: - """ - A base class for building beam emission line shapes. - - :param Line line: The emission line object for this line shape. - :param float wavelength: The rest wavelength for this emission line. - :param Beam beam: The beam class that is emitting. - """ - - def __init__(self, Line line, double wavelength, Beam beam): - - self.line = line - self.wavelength = wavelength - self.beam = beam - - cpdef Spectrum add_line(self, double radiance, Point3D beam_point, Point3D plasma_point, - Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum): - raise NotImplementedError('Child lineshape class must implement this method.') - - -DEF STARK_SPLITTING_FACTOR = 2.77e-8 - - -cdef class BeamEmissionMultiplet(BeamLineShapeModel): - """ - Produces Beam Emission Multiplet line shape, also known as the Motional Stark Effect spectrum. - """ - - def __init__(self, Line line, double wavelength, Beam beam, object sigma_to_pi, - object sigma1_to_sigma0, object pi2_to_pi3, object pi4_to_pi3): - - super().__init__(line, wavelength, beam) - - self._sigma_to_pi = autowrap_function2d(sigma_to_pi) - self._sigma1_to_sigma0 = autowrap_function1d(sigma1_to_sigma0) - self._pi2_to_pi3 = autowrap_function1d(pi2_to_pi3) - self._pi4_to_pi3 = autowrap_function1d(pi4_to_pi3) - - @cython.cdivision(True) - cpdef Spectrum add_line(self, double radiance, Point3D beam_point, Point3D plasma_point, - Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum): - - cdef double x, y, z - cdef Plasma plasma - cdef double te, ne, beam_energy, sigma, stark_split, beam_ion_mass, beam_temperature - cdef double natural_wavelength, central_wavelength - cdef double sigma_to_pi, d, intensity_sig, intensity_pi, e_field - cdef double s1_to_s0, intensity_s0, intensity_s1 - cdef double pi2_to_pi3, pi4_to_pi3, intensity_pi2, intensity_pi3, intensity_pi4 - cdef Vector3D b_field, beam_velocity - - # extract for more compact code - x = plasma_point.x - y = plasma_point.y - z = plasma_point.z - - plasma = self.beam.get_plasma() - - te = plasma.get_electron_distribution().effective_temperature(x, y, z) - if te <= 0.0: - return spectrum - - ne = plasma.get_electron_distribution().density(x, y, z) - if ne <= 0.0: - return spectrum - - beam_energy = self.beam.get_energy() - - # calculate Stark splitting - b_field = plasma.get_b_field().evaluate(x, y, z) - beam_velocity = beam_direction.normalise().mul(evamu_to_ms(beam_energy)) - e_field = beam_velocity.cross(b_field).get_length() - stark_split = fabs(STARK_SPLITTING_FACTOR * e_field) # TODO - calculate splitting factor? Reject other lines? - - # calculate emission line central wavelength, doppler shifted along observation direction - natural_wavelength = self.wavelength - central_wavelength = doppler_shift(natural_wavelength, observation_direction, beam_velocity) - - # calculate doppler broadening - beam_ion_mass = self.beam.get_element().atomic_weight - beam_temperature = self.beam.get_temperature() - sigma = thermal_broadening(self.wavelength, beam_temperature, beam_ion_mass) - - # calculate relative intensities of sigma and pi lines - sigma_to_pi = self._sigma_to_pi.evaluate(ne, beam_energy) - d = 1 / (1 + sigma_to_pi) - intensity_sig = sigma_to_pi * d * radiance - intensity_pi = 0.5 * d * radiance - - # add Sigma lines to output - s1_to_s0 = self._sigma1_to_sigma0.evaluate(ne) - intensity_s0 = 1 / (s1_to_s0 + 1) - intensity_s1 = 0.5 * s1_to_s0 * intensity_s0 - - spectrum = add_gaussian_line(intensity_sig * intensity_s0, central_wavelength, sigma, spectrum) - spectrum = add_gaussian_line(intensity_sig * intensity_s1, central_wavelength + stark_split, sigma, spectrum) - spectrum = add_gaussian_line(intensity_sig * intensity_s1, central_wavelength - stark_split, sigma, spectrum) - - # add Pi lines to output - pi2_to_pi3 = self._pi2_to_pi3.evaluate(ne) - pi4_to_pi3 = self._pi4_to_pi3.evaluate(ne) - intensity_pi3 = 1 / (1 + pi2_to_pi3 + pi4_to_pi3) - intensity_pi2 = pi2_to_pi3 * intensity_pi3 - intensity_pi4 = pi4_to_pi3 * intensity_pi3 - - spectrum = add_gaussian_line(intensity_pi * intensity_pi2, central_wavelength + 2 * stark_split, sigma, spectrum) - spectrum = add_gaussian_line(intensity_pi * intensity_pi2, central_wavelength - 2 * stark_split, sigma, spectrum) - spectrum = add_gaussian_line(intensity_pi * intensity_pi3, central_wavelength + 3 * stark_split, sigma, spectrum) - spectrum = add_gaussian_line(intensity_pi * intensity_pi3, central_wavelength - 3 * stark_split, sigma, spectrum) - spectrum = add_gaussian_line(intensity_pi * intensity_pi4, central_wavelength + 4 * stark_split, sigma, spectrum) - spectrum = add_gaussian_line(intensity_pi * intensity_pi4, central_wavelength - 4 * stark_split, sigma, spectrum) - - return spectrum diff --git a/cherab/core/model/lineshape/__init__.pxd b/cherab/core/model/lineshape/__init__.pxd new file mode 100644 index 00000000..1eeb4f7c --- /dev/null +++ b/cherab/core/model/lineshape/__init__.pxd @@ -0,0 +1,25 @@ +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from cherab.core.model.lineshape.beam cimport BeamLineShapeModel, BeamEmissionMultiplet +from cherab.core.model.lineshape.base cimport LineShapeModel +from cherab.core.model.lineshape.doppler cimport doppler_shift, thermal_broadening +from cherab.core.model.lineshape.gaussian cimport add_gaussian_line, GaussianLine +from cherab.core.model.lineshape.multiplet cimport MultipletLineShape +from cherab.core.model.lineshape.stark cimport StarkBroadenedLine +from cherab.core.model.lineshape.zeeman cimport ZeemanLineShapeModel, ZeemanTriplet, ParametrisedZeemanTriplet, ZeemanMultiplet diff --git a/cherab/core/model/lineshape/__init__.py b/cherab/core/model/lineshape/__init__.py new file mode 100644 index 00000000..30100403 --- /dev/null +++ b/cherab/core/model/lineshape/__init__.py @@ -0,0 +1,25 @@ +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from .beam import BeamLineShapeModel, BeamEmissionMultiplet +from .base import LineShapeModel +from .doppler import doppler_shift, thermal_broadening +from .gaussian import add_gaussian_line, GaussianLine +from .multiplet import MultipletLineShape +from .stark import add_lorentzian_line, StarkBroadenedLine +from .zeeman import ZeemanLineShapeModel, ZeemanTriplet, ParametrisedZeemanTriplet, ZeemanMultiplet diff --git a/cherab/core/model/lineshape/base.pxd b/cherab/core/model/lineshape/base.pxd new file mode 100644 index 00000000..b6267d45 --- /dev/null +++ b/cherab/core/model/lineshape/base.pxd @@ -0,0 +1,39 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from raysect.optical cimport Spectrum, Point3D, Vector3D +from cherab.core.atomic cimport Line +from cherab.core.species cimport Species +from cherab.core.plasma cimport Plasma +from cherab.core.atomic cimport AtomicData +from cherab.core.math.integrators cimport Integrator1D + + +cdef class LineShapeModel: + + cdef: + Line line + double wavelength + Species target_species + Plasma plasma + AtomicData atomic_data + Integrator1D integrator + + cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum) diff --git a/cherab/core/model/lineshape/base.pyx b/cherab/core/model/lineshape/base.pyx new file mode 100644 index 00000000..0b0cf7a2 --- /dev/null +++ b/cherab/core/model/lineshape/base.pyx @@ -0,0 +1,45 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + + +cdef class LineShapeModel: + """ + A base class for building line shapes. + + :param Line line: The emission line object for this line shape. + :param float wavelength: The rest wavelength for this emission line. + :param Species target_species: The target plasma species that is emitting. + :param Plasma plasma: The emitting plasma object. + :param AtomicData atomic_data: The atomic data provider. + :param Integrator1D integrator: Integrator1D instance to integrate the line shape + over the spectral bin. Default is None. + """ + + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, AtomicData atomic_data, Integrator1D integrator=None): + + self.line = line + self.wavelength = wavelength + self.target_species = target_species + self.plasma = plasma + self.atomic_data = atomic_data + self.integrator = integrator + + cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): + raise NotImplementedError('Child lineshape class must implement this method.') diff --git a/cherab/core/model/lineshape/beam/__init__.pxd b/cherab/core/model/lineshape/beam/__init__.pxd new file mode 100644 index 00000000..6172d6c3 --- /dev/null +++ b/cherab/core/model/lineshape/beam/__init__.pxd @@ -0,0 +1,20 @@ +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from cherab.core.model.lineshape.beam.base cimport * +from cherab.core.model.lineshape.beam.mse cimport * diff --git a/cherab/core/model/lineshape/beam/__init__.py b/cherab/core/model/lineshape/beam/__init__.py new file mode 100644 index 00000000..48d34582 --- /dev/null +++ b/cherab/core/model/lineshape/beam/__init__.py @@ -0,0 +1,20 @@ +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from .base import BeamLineShapeModel +from .mse import BeamEmissionMultiplet diff --git a/cherab/core/model/lineshape/beam/base.pxd b/cherab/core/model/lineshape/beam/base.pxd new file mode 100644 index 00000000..e5633446 --- /dev/null +++ b/cherab/core/model/lineshape/beam/base.pxd @@ -0,0 +1,37 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from raysect.optical cimport Spectrum, Point3D, Vector3D +from cherab.core.atomic cimport Line +from cherab.core.beam cimport Beam +from cherab.core.atomic cimport AtomicData + + +cdef class BeamLineShapeModel: + + cdef: + + Line line + double wavelength + Beam beam + AtomicData atomic_data + + cpdef Spectrum add_line(self, double radiance, Point3D beam_point, Point3D plasma_point, + Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum) diff --git a/cherab/core/model/lineshape/beam/base.pyx b/cherab/core/model/lineshape/beam/base.pyx new file mode 100644 index 00000000..d333dd6e --- /dev/null +++ b/cherab/core/model/lineshape/beam/base.pyx @@ -0,0 +1,41 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + + +cdef class BeamLineShapeModel: + """ + A base class for building beam emission line shapes. + + :param Line line: The emission line object for this line shape. + :param float wavelength: The rest wavelength for this emission line. + :param Beam beam: The beam class that is emitting. + :param AtomicData atomic_data: The atomic data provider. + """ + + def __init__(self, Line line, double wavelength, Beam beam, AtomicData atomic_data): + + self.line = line + self.wavelength = wavelength + self.beam = beam + self.atomic_data = atomic_data + + cpdef Spectrum add_line(self, double radiance, Point3D beam_point, Point3D plasma_point, + Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum): + raise NotImplementedError('Child lineshape class must implement this method.') diff --git a/cherab/core/model/lineshape/beam/mse.pxd b/cherab/core/model/lineshape/beam/mse.pxd new file mode 100644 index 00000000..a7b6131a --- /dev/null +++ b/cherab/core/model/lineshape/beam/mse.pxd @@ -0,0 +1,30 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from cherab.core.math cimport Function1D, Function2D +from cherab.core.model.lineshape.beam.base cimport BeamLineShapeModel + + +cdef class BeamEmissionMultiplet(BeamLineShapeModel): + + cdef: + + Function2D _sigma_to_pi + Function1D _sigma1_to_sigma0, _pi2_to_pi3, _pi4_to_pi3 diff --git a/cherab/core/model/lineshape/beam/mse.pyx b/cherab/core/model/lineshape/beam/mse.pyx new file mode 100644 index 00000000..9c64d435 --- /dev/null +++ b/cherab/core/model/lineshape/beam/mse.pyx @@ -0,0 +1,135 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from libc.math cimport fabs, sqrt +from raysect.optical cimport Spectrum, Point3D, Vector3D + +from cherab.core.plasma cimport Plasma +from cherab.core.beam cimport Beam +from cherab.core.atomic cimport AtomicData +from cherab.core.atomic cimport Line +from cherab.core.math.function cimport autowrap_function1d, autowrap_function2d +from cherab.core.utility.constants cimport ATOMIC_MASS, ELEMENTARY_CHARGE +from cherab.core.model.lineshape.gaussian cimport add_gaussian_line +from cherab.core.model.lineshape.doppler cimport thermal_broadening, doppler_shift + +cimport cython + + +cdef double RECIP_ATOMIC_MASS = 1 / ATOMIC_MASS + + +cdef double evamu_to_ms(double x): + return sqrt(2 * x * ELEMENTARY_CHARGE * RECIP_ATOMIC_MASS) + + +DEF STARK_SPLITTING_FACTOR = 2.77e-8 + + +cdef class BeamEmissionMultiplet(BeamLineShapeModel): + """ + Produces Beam Emission Multiplet line shape, also known as the Motional Stark Effect spectrum. + """ + + def __init__(self, Line line, double wavelength, Beam beam, AtomicData atomic_data, + object sigma_to_pi, object sigma1_to_sigma0, object pi2_to_pi3, object pi4_to_pi3): + + super().__init__(line, wavelength, beam, atomic_data) + + self._sigma_to_pi = autowrap_function2d(sigma_to_pi) + self._sigma1_to_sigma0 = autowrap_function1d(sigma1_to_sigma0) + self._pi2_to_pi3 = autowrap_function1d(pi2_to_pi3) + self._pi4_to_pi3 = autowrap_function1d(pi4_to_pi3) + + @cython.cdivision(True) + cpdef Spectrum add_line(self, double radiance, Point3D beam_point, Point3D plasma_point, + Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum): + + cdef double x, y, z + cdef Plasma plasma + cdef double te, ne, beam_energy, sigma, stark_split, beam_ion_mass, beam_temperature + cdef double natural_wavelength, central_wavelength + cdef double sigma_to_pi, d, intensity_sig, intensity_pi, e_field + cdef double s1_to_s0, intensity_s0, intensity_s1 + cdef double pi2_to_pi3, pi4_to_pi3, intensity_pi2, intensity_pi3, intensity_pi4 + cdef Vector3D b_field, beam_velocity + + # extract for more compact code + x = plasma_point.x + y = plasma_point.y + z = plasma_point.z + + plasma = self.beam.get_plasma() + + te = plasma.get_electron_distribution().effective_temperature(x, y, z) + if te <= 0.0: + return spectrum + + ne = plasma.get_electron_distribution().density(x, y, z) + if ne <= 0.0: + return spectrum + + beam_energy = self.beam.get_energy() + + # calculate Stark splitting + b_field = plasma.get_b_field().evaluate(x, y, z) + beam_velocity = beam_direction.normalise().mul(evamu_to_ms(beam_energy)) + e_field = beam_velocity.cross(b_field).get_length() + stark_split = fabs(STARK_SPLITTING_FACTOR * e_field) # TODO - calculate splitting factor? Reject other lines? + + # calculate emission line central wavelength, doppler shifted along observation direction + natural_wavelength = self.wavelength + central_wavelength = doppler_shift(natural_wavelength, observation_direction, beam_velocity) + + # calculate doppler broadening + beam_ion_mass = self.beam.get_element().atomic_weight + beam_temperature = self.beam.get_temperature() + sigma = thermal_broadening(self.wavelength, beam_temperature, beam_ion_mass) + + # calculate relative intensities of sigma and pi lines + sigma_to_pi = self._sigma_to_pi.evaluate(ne, beam_energy) + d = 1 / (1 + sigma_to_pi) + intensity_sig = sigma_to_pi * d * radiance + intensity_pi = 0.5 * d * radiance + + # add Sigma lines to output + s1_to_s0 = self._sigma1_to_sigma0.evaluate(ne) + intensity_s0 = 1 / (s1_to_s0 + 1) + intensity_s1 = 0.5 * s1_to_s0 * intensity_s0 + + spectrum = add_gaussian_line(intensity_sig * intensity_s0, central_wavelength, sigma, spectrum) + spectrum = add_gaussian_line(intensity_sig * intensity_s1, central_wavelength + stark_split, sigma, spectrum) + spectrum = add_gaussian_line(intensity_sig * intensity_s1, central_wavelength - stark_split, sigma, spectrum) + + # add Pi lines to output + pi2_to_pi3 = self._pi2_to_pi3.evaluate(ne) + pi4_to_pi3 = self._pi4_to_pi3.evaluate(ne) + intensity_pi3 = 1 / (1 + pi2_to_pi3 + pi4_to_pi3) + intensity_pi2 = pi2_to_pi3 * intensity_pi3 + intensity_pi4 = pi4_to_pi3 * intensity_pi3 + + spectrum = add_gaussian_line(intensity_pi * intensity_pi2, central_wavelength + 2 * stark_split, sigma, spectrum) + spectrum = add_gaussian_line(intensity_pi * intensity_pi2, central_wavelength - 2 * stark_split, sigma, spectrum) + spectrum = add_gaussian_line(intensity_pi * intensity_pi3, central_wavelength + 3 * stark_split, sigma, spectrum) + spectrum = add_gaussian_line(intensity_pi * intensity_pi3, central_wavelength - 3 * stark_split, sigma, spectrum) + spectrum = add_gaussian_line(intensity_pi * intensity_pi4, central_wavelength + 4 * stark_split, sigma, spectrum) + spectrum = add_gaussian_line(intensity_pi * intensity_pi4, central_wavelength - 4 * stark_split, sigma, spectrum) + + return spectrum diff --git a/cherab/core/model/lineshape/doppler.pxd b/cherab/core/model/lineshape/doppler.pxd new file mode 100644 index 00000000..3d4438fc --- /dev/null +++ b/cherab/core/model/lineshape/doppler.pxd @@ -0,0 +1,26 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from raysect.optical cimport Vector3D + + +cpdef double doppler_shift(double wavelength, Vector3D observation_direction, Vector3D velocity) + +cpdef double thermal_broadening(double wavelength, double temperature, double atomic_weight) diff --git a/cherab/core/model/lineshape/doppler.pyx b/cherab/core/model/lineshape/doppler.pyx new file mode 100644 index 00000000..93755733 --- /dev/null +++ b/cherab/core/model/lineshape/doppler.pyx @@ -0,0 +1,59 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from libc.math cimport sqrt + +from cherab.core.utility.constants cimport ATOMIC_MASS, ELEMENTARY_CHARGE, SPEED_OF_LIGHT + +cimport cython + + +@cython.cdivision(True) +cpdef double doppler_shift(double wavelength, Vector3D observation_direction, Vector3D velocity): + """ + Calculates the Doppler shifted wavelength for a given velocity and observation direction. + + :param wavelength: The wavelength to Doppler shift in nanometers. + :param observation_direction: A Vector defining the direction of observation. + :param velocity: A Vector defining the relative velocity of the emitting source in m/s. + :return: The Doppler shifted wavelength in nanometers. + """ + cdef double projected_velocity + + # flow velocity projected on the direction of observation + observation_direction = observation_direction.normalise() + projected_velocity = velocity.dot(observation_direction) + + return wavelength * (1 + projected_velocity / SPEED_OF_LIGHT) + + +@cython.cdivision(True) +cpdef double thermal_broadening(double wavelength, double temperature, double atomic_weight): + """ + Returns the line width for a gaussian line as a standard deviation. + + :param wavelength: Central wavelength. + :param temperature: Temperature in eV. + :param atomic_weight: Atomic weight in AMU. + :return: Standard deviation of gaussian line. + """ + + # todo: add input sanity checks + return sqrt(temperature * ELEMENTARY_CHARGE / (atomic_weight * ATOMIC_MASS)) * wavelength / SPEED_OF_LIGHT diff --git a/cherab/core/model/lineshape/gaussian.pxd b/cherab/core/model/lineshape/gaussian.pxd new file mode 100644 index 00000000..cd8d0091 --- /dev/null +++ b/cherab/core/model/lineshape/gaussian.pxd @@ -0,0 +1,29 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from raysect.optical cimport Spectrum +from cherab.core.model.lineshape.base cimport LineShapeModel + + +cpdef Spectrum add_gaussian_line(double radiance, double wavelength, double sigma, Spectrum spectrum) + + +cdef class GaussianLine(LineShapeModel): + pass diff --git a/cherab/core/model/lineshape/gaussian.pyx b/cherab/core/model/lineshape/gaussian.pyx new file mode 100644 index 00000000..ef854ea9 --- /dev/null +++ b/cherab/core/model/lineshape/gaussian.pyx @@ -0,0 +1,139 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from libc.math cimport erf, M_SQRT2, floor, ceil +from raysect.optical cimport Point3D, Vector3D + +from cherab.core.atomic cimport Line, AtomicData +from cherab.core.species cimport Species +from cherab.core.plasma cimport Plasma +from cherab.core.model.lineshape.doppler cimport doppler_shift, thermal_broadening + +cimport cython + + +# the number of standard deviations outside the rest wavelength the line is considered to add negligible value (including a margin for safety) +DEF GAUSSIAN_CUTOFF_SIGMA = 10.0 + + +@cython.cdivision(True) +@cython.initializedcheck(False) +@cython.boundscheck(False) +@cython.wraparound(False) +cpdef Spectrum add_gaussian_line(double radiance, double wavelength, double sigma, Spectrum spectrum): + r""" + Adds a Gaussian line to the given spectrum and returns the new spectrum. + + The formula used is based on the following definite integral: + :math:`\frac{1}{\sigma \sqrt{2 \pi}} \int_{\lambda_0}^{\lambda_1} \exp(-\frac{(x-\mu)^2}{2\sigma^2}) dx = \frac{1}{2} \left[ -Erf(\frac{a-\mu}{\sqrt{2}\sigma}) +Erf(\frac{b-\mu}{\sqrt{2}\sigma}) \right]` + + :param float radiance: Intensity of the line in radiance. + :param float wavelength: central wavelength of the line in nm. + :param float sigma: width of the line in nm. + :param Spectrum spectrum: the current spectrum to which the gaussian line is added. + :return: + """ + + cdef double temp + cdef double cutoff_lower_wavelength, cutoff_upper_wavelength + cdef double lower_wavelength, upper_wavelength + cdef double lower_integral, upper_integral + cdef int start, end, i + + if sigma <= 0: + return spectrum + + # calculate and check end of limits + cutoff_lower_wavelength = wavelength - GAUSSIAN_CUTOFF_SIGMA * sigma + if spectrum.max_wavelength < cutoff_lower_wavelength: + return spectrum + + cutoff_upper_wavelength = wavelength + GAUSSIAN_CUTOFF_SIGMA * sigma + if spectrum.min_wavelength > cutoff_upper_wavelength: + return spectrum + + # locate range of bins where there is significant contribution from the gaussian (plus a health margin) + start = max(0, floor((cutoff_lower_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) + end = min(spectrum.bins, ceil((cutoff_upper_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) + + # add line to spectrum + temp = 1 / (M_SQRT2 * sigma) + lower_wavelength = spectrum.min_wavelength + start * spectrum.delta_wavelength + lower_integral = erf((lower_wavelength - wavelength) * temp) + for i in range(start, end): + + upper_wavelength = spectrum.min_wavelength + spectrum.delta_wavelength * (i + 1) + upper_integral = erf((upper_wavelength - wavelength) * temp) + + spectrum.samples_mv[i] += radiance * 0.5 * (upper_integral - lower_integral) / spectrum.delta_wavelength + + lower_wavelength = upper_wavelength + lower_integral = upper_integral + + return spectrum + + +cdef class GaussianLine(LineShapeModel): + """ + Produces Gaussian line shape. + + :param Line line: The emission line object for this line shape. + :param float wavelength: The rest wavelength for this emission line. + :param Species target_species: The target plasma species that is emitting. + :param Plasma plasma: The emitting plasma object. + :param AtomicData atomic_data: The atomic data provider. + + .. code-block:: pycon + + >>> from cherab.core.atomic import Line, deuterium + >>> from cherab.core.model import ExcitationLine, GaussianLine + >>> + >>> # Adding Gaussian line to the plasma model. + >>> d_alpha = Line(deuterium, 0, (3, 2)) + >>> excit = ExcitationLine(d_alpha, lineshape=GaussianLine) + >>> plasma.models.add(excit) + """ + + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, AtomicData atomic_data): + + super().__init__(line, wavelength, target_species, plasma, atomic_data) + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + @cython.cdivision(True) + cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): + + cdef double ts, sigma, shifted_wavelength + cdef Vector3D ion_velocity + + ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) + if ts <= 0.0: + return spectrum + + ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) + + # calculate emission line central wavelength, doppler shifted along observation direction + shifted_wavelength = doppler_shift(self.wavelength, direction, ion_velocity) + + # calculate the line width + sigma = thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) + + return add_gaussian_line(radiance, shifted_wavelength, sigma, spectrum) diff --git a/cherab/core/model/lineshape/multiplet.pxd b/cherab/core/model/lineshape/multiplet.pxd new file mode 100644 index 00000000..577afadc --- /dev/null +++ b/cherab/core/model/lineshape/multiplet.pxd @@ -0,0 +1,31 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +cimport numpy as np + +from cherab.core.model.lineshape.base cimport LineShapeModel + + +cdef class MultipletLineShape(LineShapeModel): + + cdef: + int _number_of_lines + np.ndarray _multiplet + double[:, ::1] _multiplet_mv diff --git a/cherab/core/model/lineshape/multiplet.pyx b/cherab/core/model/lineshape/multiplet.pyx new file mode 100644 index 00000000..6f9e39ab --- /dev/null +++ b/cherab/core/model/lineshape/multiplet.pyx @@ -0,0 +1,117 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +import numpy as np + +from raysect.optical cimport Spectrum, Point3D, Vector3D + +from cherab.core.atomic cimport Line, AtomicData +from cherab.core.species cimport Species +from cherab.core.plasma cimport Plasma +from cherab.core.model.lineshape.doppler cimport doppler_shift, thermal_broadening +from cherab.core.model.lineshape.gaussian cimport add_gaussian_line + +cimport cython + +# required by numpy c-api +np.import_array() + + +DEF MULTIPLET_WAVELENGTH = 0 +DEF MULTIPLET_RATIO = 1 + + +cdef class MultipletLineShape(LineShapeModel): + """ + Produces Multiplet line shapes. + + The lineshape radiance is calculated from a base PEC rate that is unresolved. This + radiance is then divided over a number of components as specified in the multiplet + argument. The multiplet components are specified with an Nx2 array where N is the + number of components in the multiplet. The first axis of the array contains the + wavelengths of each component, the second contains the line ratio for each component. + The component line ratios must sum to one. For example: + + :param Line line: The emission line object for the base rate radiance calculation. + :param float wavelength: The rest wavelength of the base emission line. + :param Species target_species: The target plasma species that is emitting. + :param Plasma plasma: The emitting plasma object. + :param AtomicData atomic_data: The atomic data provider. + :param multiplet: An Nx2 array that specifies the multiplet wavelengths and line ratios. + + .. code-block:: pycon + + >>> from cherab.core.atomic import Line, nitrogen + >>> from cherab.core.model import ExcitationLine, MultipletLineShape + >>> + >>> # multiplet specification in Nx2 array + >>> multiplet = [[403.509, 404.132, 404.354, 404.479, 405.692], [0.205, 0.562, 0.175, 0.029, 0.029]] + >>> + >>> # Adding the multiplet to the plasma model. + >>> nitrogen_II_404 = Line(nitrogen, 1, ("2s2 2p1 4f1 3G13.0", "2s2 2p1 3d1 3F10.0")) + >>> excit = ExcitationLine(nitrogen_II_404, lineshape=MultipletLineShape, lineshape_args=[multiplet]) + >>> plasma.models.add(excit) + """ + + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, AtomicData atomic_data, + object multiplet): + + super().__init__(line, wavelength, target_species, plasma, atomic_data) + + multiplet = np.array(multiplet, dtype=np.float64) + + if not (len(multiplet.shape) == 2 and multiplet.shape[0] == 2): + raise ValueError("The multiplet specification must be an array of shape (Nx2).") + + if not multiplet[1, :].sum() == 1.0: + raise ValueError("The multiplet line ratios should sum to one.") + + self._number_of_lines = multiplet.shape[1] + self._multiplet = multiplet + self._multiplet_mv = self._multiplet + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): + + cdef double ts, sigma, shifted_wavelength, component_wavelength, component_radiance + cdef Vector3D ion_velocity + + ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) + if ts <= 0.0: + return spectrum + + ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) + + # calculate the line width + sigma = thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) + + for i in range(self._number_of_lines): + + component_wavelength = self._multiplet_mv[MULTIPLET_WAVELENGTH, i] + component_radiance = radiance * self._multiplet_mv[MULTIPLET_RATIO, i] + + # calculate emission line central wavelength, doppler shifted along observation direction + shifted_wavelength = doppler_shift(component_wavelength, direction, ion_velocity) + + spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) + + return spectrum diff --git a/cherab/core/model/lineshape/stark.pxd b/cherab/core/model/lineshape/stark.pxd new file mode 100644 index 00000000..8b4ac4b0 --- /dev/null +++ b/cherab/core/model/lineshape/stark.pxd @@ -0,0 +1,36 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from raysect.optical cimport Spectrum +from cherab.core.math.integrators cimport Integrator1D +from cherab.core.model.lineshape.zeeman cimport ZeemanLineShapeModel + + +cpdef Spectrum add_lorentzian_line(double radiance, double wavelength, double lambda_1_2, Spectrum spectrum, + Integrator1D integrator) + + +cdef class StarkBroadenedLine(ZeemanLineShapeModel): + + cdef: + double _aij, _bij, _cij + double _fwhm_poly_coeff_gauss[7] + double _fwhm_poly_coeff_lorentz[7] + double _weight_poly_coeff[6] diff --git a/cherab/core/model/lineshape/stark.pyx b/cherab/core/model/lineshape/stark.pyx new file mode 100644 index 00000000..9e333a5f --- /dev/null +++ b/cherab/core/model/lineshape/stark.pyx @@ -0,0 +1,348 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from scipy.special import hyp2f1 + +from libc.math cimport sqrt, floor, ceil, fabs, log, exp +from raysect.core.math.function.float cimport Function1D +from raysect.optical cimport Point3D, Vector3D + +from cherab.core.atomic cimport Line, AtomicData +from cherab.core.species cimport Species +from cherab.core.plasma cimport Plasma +from cherab.core.atomic.elements import hydrogen, deuterium, tritium +from cherab.core.math.function cimport autowrap_function1d, autowrap_function2d +from cherab.core.math.integrators cimport GaussianQuadrature +from cherab.core.utility.constants cimport BOHR_MAGNETON, HC_EV_NM +from cherab.core.model.lineshape.doppler cimport doppler_shift, thermal_broadening +from cherab.core.model.lineshape.gaussian cimport add_gaussian_line + + +cimport cython + + +DEF PI_POLARISATION = 0 +DEF SIGMA_POLARISATION = 1 +DEF SIGMA_PLUS_POLARISATION = 1 +DEF SIGMA_MINUS_POLARISATION = -1 +DEF NO_POLARISATION = 2 + +DEF LORENTZIAN_CUTOFF_GAMMA = 50.0 + +cdef double _SIGMA2FWHM = 2 * sqrt(2 * log(2)) + + +cdef class StarkFunction(Function1D): + """ + Normalised Stark function for the StarkBroadenedLine line shape. + + :param float wavelength: central wavelength of the line in nm. + :param float lambda_1_2: FWHM of the function in nm. + """ + + cdef double _a, _x0, _norm + + STARK_NORM_COEFFICIENT = 4 * LORENTZIAN_CUTOFF_GAMMA * hyp2f1(0.4, 1, 1.4, -(2 * LORENTZIAN_CUTOFF_GAMMA)**2.5) + + def __init__(self, double wavelength, double lambda_1_2): + + if wavelength <= 0: + raise ValueError("Argument 'wavelength' must be positive.") + + if lambda_1_2 <= 0: + raise ValueError("Argument 'lambda_1_2' must be positive.") + + self._x0 = wavelength + self._a = (0.5 * lambda_1_2)**2.5 + # normalise, so the integral over x is equal to 1 in the limits + # (_x0 - LORENTZIAN_CUTOFF_GAMMA * lambda_1_2, _x0 + LORENTZIAN_CUTOFF_GAMMA * lambda_1_2) + self._norm = (0.5 * lambda_1_2)**1.5 / self.STARK_NORM_COEFFICIENT + + @cython.cdivision(True) + cdef double evaluate(self, double x) except? -1e999: + + return self._norm / ((fabs(x - self._x0))**2.5 + self._a) + + +@cython.cdivision(True) +@cython.initializedcheck(False) +@cython.boundscheck(False) +@cython.wraparound(False) +cpdef Spectrum add_lorentzian_line(double radiance, double wavelength, double lambda_1_2, Spectrum spectrum, Integrator1D integrator): + r""" + Adds a modified Lorentzian line to the given spectrum and returns the new spectrum. + + The modified Lorentzian: + :math:`L(\lambda-\lambda_0, \Delta\lambda_{1/2}^{(L)})=\frac{C_0(\Delta\lambda_{1/2}^{(L)})^{3/2}}{(\lambda-\lambda_0)^{5/2}+(\frac{\Delta\lambda_{1/2}^{(L)}}{2})^{5/2}},` + + :math:`C_0=\frac{(1/2)^{3/2}}{4R{_2}F_1(\frac{2}{5},1,\frac{7}{5},-R^{5/2})},` + + where :math:`\Delta\lambda_{1/2}^{(L)}` is the line FWHM, :math:`{_2}F_1` is the hypergeometric function and :math:`R=50`. + The line shape is truncated at :math:`\lambdaR\Delta\lambda_{1/2}^{(L)}`. + + See B. Lomanowski, et al. "Inferring divertor plasma properties from hydrogen Balmer + and Paschen series spectroscopy in JET-ILW." Nuclear Fusion 55.12 (2015) + `123028 `_ for details. + + :param float radiance: Intensity of the line in radiance. + :param float wavelength: central wavelength of the line in nm. + :param float lambda_1_2: FWHM of the line shape in nm. + :param Spectrum spectrum: the current spectrum to which the Lorentzian line is added. + :param Integrator1D integrator: Integrator1D instance to integrate the line shape + over the spectral bin. + :return: + """ + + cdef double cutoff_lower_wavelength, cutoff_upper_wavelength + cdef double lower_wavelength, upper_wavelength + cdef double bin_integral + cdef int start, end, i + + if lambda_1_2 <= 0: + return spectrum + + integrator.function = StarkFunction(wavelength, lambda_1_2) + + # calculate and check end of limits + cutoff_lower_wavelength = wavelength - LORENTZIAN_CUTOFF_GAMMA * lambda_1_2 + if spectrum.max_wavelength < cutoff_lower_wavelength: + return spectrum + + cutoff_upper_wavelength = wavelength + LORENTZIAN_CUTOFF_GAMMA * lambda_1_2 + if spectrum.min_wavelength > cutoff_upper_wavelength: + return spectrum + + # locate range of bins where there is significant contribution from the gaussian (plus a health margin) + start = max(0, floor((cutoff_lower_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) + end = min(spectrum.bins, ceil((cutoff_upper_wavelength - spectrum.min_wavelength) / spectrum.delta_wavelength)) + + # add line to spectrum + lower_wavelength = spectrum.min_wavelength + start * spectrum.delta_wavelength + + for i in range(start, end): + upper_wavelength = spectrum.min_wavelength + spectrum.delta_wavelength * (i + 1) + + bin_integral = integrator.evaluate(lower_wavelength, upper_wavelength) + spectrum.samples_mv[i] += radiance * bin_integral / spectrum.delta_wavelength + + lower_wavelength = upper_wavelength + + return spectrum + + +cdef class StarkBroadenedLine(ZeemanLineShapeModel): + r""" + Parametrised Stark-Doppler-Zeeman line shape for Balmer and Paschen series based on + B. Lomanowski, et al. "Inferring divertor plasma properties from hydrogen Balmer + and Paschen series spectroscopy in JET-ILW." Nuclear Fusion 55.12 (2015) + `123028 `_. + + The following approximations are used: + + * The Zeeman and Stark effects are considered independently. + * Zeeman splitting is taken in the form of a simple triplet with a :math:`\pi`-component + centred at :math:`\lambda`, :math:`\sigma^{+}`-component at :math:`\frac{hc}{hc/\lambda -\mu B}` + and :math:`\sigma^{-}`-component at :math:`\frac{hc}{hc/\lambda +\mu B}`. + * The model of Stark broadening is obtained by fitting the Model Microfield Method (MMM). + * The convolution of Stark-Zeeman and Doppler profiles is replaced with the weighted sum + to speed-up calculations (so-called pseudo-Voigt profile). + + The Stark-broadened line shape is modelled as modified Lorentzian: + :math:`L(\lambda-\lambda_0, \Delta\lambda_{1/2}^{(L)})=\frac{C_0(\Delta\lambda_{1/2}^{(L)})^{3/2}}{(\lambda-\lambda_0)^{5/2}+(\frac{\Delta\lambda_{1/2}^{(L)}}{2})^{5/2}},` + + :math:`C_0=\frac{(1/2)^{3/2}}{4R{_2}F_1(\frac{2}{5},1,\frac{7}{5},-R^{5/2})},` + + where :math:`\Delta\lambda_{1/2}^{(L)}=c_{ij}\frac{n_e^{a_{ij}}}{T_e^{b_{ij}}}` is the line FWHM, + :math:`{_2}F_1` is the hypergeometric function and :math:`R=50`. + The line shape is truncated at :math:`\lambdaR\Delta\lambda_{1/2}^{(L)}`. + The :math:`a_{ij}`, :math:`b_{ij}` and :math:`c_{ij}` are the fitting coefficients. + + Each Zeeman component is modelled as a weighted sum of Stark (Lorentzian), :math:`L(\lambda-\lambda_0, \Delta\lambda_{1/2}^{(L)})`, + and Doppler (Gauss), :math:`G(\lambda'-\lambda_0, \Delta\lambda_{1/2}^{(G)})` profiles: + + :math:`\eta L(\lambda-\lambda_0, \Delta\lambda_{1/2}^{(V)}) + (1-\eta)G(\lambda-\lambda_0, \Delta\lambda_{1/2}^{(V)})`, + + with :math:`\Delta\lambda_{1/2}^{(V)}\equiv \Delta\lambda_{1/2}^{(V)}(\Delta\lambda_{1/2}^{(G)}, \Delta\lambda_{1/2}^{(L)})` + and :math:`\eta\equiv \eta(\Delta\lambda_{1/2}^{(G)}, \Delta\lambda_{1/2}^{(L)})`. + + Both :math:`\Delta\lambda_{1/2}^{(V)}` and :math:`\eta` are obtained by fitting the convolution: + + :math:`\int_{-\infty}^{\infty}G(\lambda'-\lambda_0, \Delta\lambda_{1/2}^{(G)})L(\lambda-\lambda'-\lambda_0, \Delta\lambda_{1/2}^{(L)})d\lambda'`. + + The :math:`\Delta\lambda_{1/2}^{(V)}` function is fitted as: + :math:`\Delta\lambda_{1/2}^{(V)}=\sum_{n=0}^{6}a_n(\frac{\Delta\lambda_{1/2}^{(L)}}{\Delta\lambda_{1/2}^{(G)}})^n` + for :math:`\frac{\Delta\lambda_{1/2}^{(L)}}{\Delta\lambda_{1/2}^{(G)}} \le 1` and + :math:`\Delta\lambda_{1/2}^{(V)}=\sum_{n=0}^{6}b_n(\frac{\Delta\lambda_{1/2}^{(G)}}{\Delta\lambda_{1/2}^{(L)}})^n` + for :math:`\frac{\Delta\lambda_{1/2}^{(L)}}{\Delta\lambda_{1/2}^{(G)}} > 1` with + + `a` = [1., 0.15882, 1.04388, -1.38281, 0.46251, 0.82325, -0.58026] and + + `b` = [1., 0, 0.57575, 0.37902, -0.42519, -0.31525, 0.31718]. + + While the :math:`\eta` function is fitted as: + :math:`\eta=exp(\sum_{n=0}^{5}c_n(ln(\frac{\Delta\lambda_{1/2}^{(L)}}{\Delta\lambda_{1/2}^{(V)}}))^n` + for :math:`0.01<\frac{\Delta\lambda_{1/2}^{(L)}}{\Delta\lambda_{1/2}^{(V)}}<0.999` with + + `c` = [5.14820e-04, 1.38821e+00, -9.60424e-02, -3.83995e-02, -7.40042e-03, -5.47626e-04]. + + :param Line line: The emission line object for this line shape. + :param float wavelength: The rest wavelength for this emission line. + :param Species target_species: The target plasma species that is emitting. + :param Plasma plasma: The emitting plasma object. + :param AtomicData atomic_data: The atomic data provider. + :param tuple stark_model_coefficients: Stark model coefficients in the form (c_ij, a_ij, b_ij). + Default is None (will use + `atomic_data.stark_model_coefficients`). + :param Integrator1D integrator: Integrator1D instance to integrate the line shape + over the spectral bin. Default is `GaussianQuadrature()`. + :param str polarisation: Leaves only :math:`\pi`-/:math:`\sigma`-polarised components: + "pi" - leave only :math:`\pi`-polarised components, + "sigma" - leave only :math:`\sigma`-polarised components, + "no" - leave all components (default). + """ + + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, AtomicData atomic_data, + tuple stark_model_coefficients=None, Integrator1D integrator=GaussianQuadrature(), polarisation='no'): + + super().__init__(line, wavelength, target_species, plasma, atomic_data, polarisation, integrator) + + try: + # Fitted Stark Constants + cij, aij, bij = stark_model_coefficients or self.atomic_data.stark_model_coefficients(line) + if cij <= 0: + raise ValueError('Coefficient c_ij must be positive.') + if aij <= 0: + raise ValueError('Coefficient a_ij must be positive.') + if bij <= 0: + raise ValueError('Coefficient b_ij must be positive.') + self._aij = aij + self._bij = bij + self._cij = cij + except IndexError: + raise ValueError('Stark broadening coefficients for {} is not currently available.'.format(line)) + + # polynomial coefficients in increasing powers + self._fwhm_poly_coeff_gauss = [1., 0, 0.57575, 0.37902, -0.42519, -0.31525, 0.31718] + self._fwhm_poly_coeff_lorentz = [1., 0.15882, 1.04388, -1.38281, 0.46251, 0.82325, -0.58026] + + self._weight_poly_coeff = [5.14820e-04, 1.38821e+00, -9.60424e-02, -3.83995e-02, -7.40042e-03, -5.47626e-04] + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + @cython.cdivision(True) + cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): + + cdef: + double ne, te, ts, shifted_wavelength, photon_energy, b_magn, comp_radiance, cos_sqr, sin_sqr + double sigma, fwhm_lorentz, fwhm_gauss, fwhm_full, fwhm_ratio, fwhm_lorentz_to_total, lorentz_weight, gauss_weight + cdef Vector3D ion_velocity, b_field + int i + + ne = self.plasma.get_electron_distribution().density(point.x, point.y, point.z) + + te = self.plasma.get_electron_distribution().effective_temperature(point.x, point.y, point.z) + + fwhm_lorentz = self._cij * ne**self._aij / (te**self._bij) if ne > 0 and te > 0 else 0 + + ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) + + fwhm_gauss = _SIGMA2FWHM * thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) if ts > 0 else 0 + + if fwhm_lorentz == 0 and fwhm_gauss == 0: + return spectrum + + # calculating full FWHM + if fwhm_gauss <= fwhm_lorentz: + fwhm_ratio = fwhm_gauss / fwhm_lorentz + fwhm_full = self._fwhm_poly_coeff_gauss[0] + for i in range(1, 7): + fwhm_full += self._fwhm_poly_coeff_gauss[i] * fwhm_ratio**i + fwhm_full *= fwhm_lorentz + else: + fwhm_ratio = fwhm_lorentz / fwhm_gauss + fwhm_full = self._fwhm_poly_coeff_lorentz[0] + for i in range(1, 7): + fwhm_full += self._fwhm_poly_coeff_lorentz[i] * fwhm_ratio**i + fwhm_full *= fwhm_gauss + + sigma = fwhm_full / _SIGMA2FWHM + + fwhm_lorentz_to_total = fwhm_lorentz / fwhm_full + + # calculating Lorentzian weight + if fwhm_lorentz_to_total < 0.01: + lorentz_weight = 0 + fwhm_full = 0 # force add_lorentzian_line() to immediately return + elif fwhm_lorentz_to_total > 0.999: + lorentz_weight = 1 + sigma = 0 # force add_gaussian_line() to immediately return + else: + lorentz_weight = self._weight_poly_coeff[0] + for i in range(1, 6): + lorentz_weight += self._weight_poly_coeff[i] * log(fwhm_lorentz_to_total)**i + lorentz_weight = exp(lorentz_weight) + + gauss_weight = 1 - lorentz_weight + + ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) + + # calculate emission line central wavelength, doppler shifted along observation direction + shifted_wavelength = doppler_shift(self.wavelength, direction, ion_velocity) + + # obtain magnetic field + b_field = self.plasma.get_b_field().evaluate(point.x, point.y, point.z) + b_magn = b_field.get_length() + + if b_magn == 0: + # no splitting if magnetic field strength is zero + if self._polarisation != NO_POLARISATION: + radiance *= 0.5 # pi or sigma polarisation, collecting only half of intensity + + spectrum = add_gaussian_line(gauss_weight * radiance, shifted_wavelength, sigma, spectrum) + spectrum = add_lorentzian_line(lorentz_weight * radiance, shifted_wavelength, fwhm_full, spectrum, self.integrator) + + return spectrum + + # coefficients for intensities parallel and perpendicular to magnetic field + cos_sqr = (b_field.dot(direction.normalise()) / b_magn)**2 + sin_sqr = 1. - cos_sqr + + # adding pi component of the Zeeman triplet in case of NO_POLARISATION or PI_POLARISATION + if self._polarisation != SIGMA_POLARISATION: + comp_radiance = 0.5 * sin_sqr * radiance + spectrum = add_gaussian_line(gauss_weight * comp_radiance, shifted_wavelength, sigma, spectrum) + spectrum = add_lorentzian_line(lorentz_weight * comp_radiance, shifted_wavelength, fwhm_full, spectrum, self.integrator) + + # adding sigma +/- components of the Zeeman triplet in case of NO_POLARISATION or SIGMA_POLARISATION + if self._polarisation != PI_POLARISATION: + comp_radiance = (0.25 * sin_sqr + 0.5 * cos_sqr) * radiance + + photon_energy = HC_EV_NM / self.wavelength + + shifted_wavelength = doppler_shift(HC_EV_NM / (photon_energy - BOHR_MAGNETON * b_magn), direction, ion_velocity) + spectrum = add_gaussian_line(gauss_weight * comp_radiance, shifted_wavelength, sigma, spectrum) + spectrum = add_lorentzian_line(lorentz_weight * comp_radiance, shifted_wavelength, fwhm_full, spectrum, self.integrator) + + shifted_wavelength = doppler_shift(HC_EV_NM / (photon_energy + BOHR_MAGNETON * b_magn), direction, ion_velocity) + spectrum = add_gaussian_line(gauss_weight * comp_radiance, shifted_wavelength, sigma, spectrum) + spectrum = add_lorentzian_line(lorentz_weight * comp_radiance, shifted_wavelength, fwhm_full, spectrum, self.integrator) + + return spectrum diff --git a/cherab/core/model/lineshape/zeeman.pxd b/cherab/core/model/lineshape/zeeman.pxd new file mode 100644 index 00000000..b67d2042 --- /dev/null +++ b/cherab/core/model/lineshape/zeeman.pxd @@ -0,0 +1,48 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from cherab.core.atomic.zeeman cimport ZeemanStructure +from cherab.core.model.lineshape.base cimport LineShapeModel + + +cdef class ZeemanLineShapeModel(LineShapeModel): + + cdef double _polarisation + + pass + + +cdef class ZeemanTriplet(ZeemanLineShapeModel): + + pass + + +cdef class ParametrisedZeemanTriplet(ZeemanLineShapeModel): + + cdef double _alpha, _beta, _gamma + + pass + + +cdef class ZeemanMultiplet(ZeemanLineShapeModel): + + cdef ZeemanStructure _zeeman_structure + + pass diff --git a/cherab/core/model/lineshape/zeeman.pyx b/cherab/core/model/lineshape/zeeman.pyx new file mode 100644 index 00000000..2be33020 --- /dev/null +++ b/cherab/core/model/lineshape/zeeman.pyx @@ -0,0 +1,365 @@ +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from libc.math cimport sqrt +from raysect.optical cimport Spectrum, Point3D, Vector3D + +from cherab.core.atomic cimport Line, AtomicData +from cherab.core.species cimport Species +from cherab.core.plasma cimport Plasma +from cherab.core.atomic.elements import hydrogen, deuterium, tritium, helium, helium3, beryllium, boron, carbon, nitrogen, oxygen, neon +from cherab.core.math.integrators cimport Integrator1D +from cherab.core.utility.constants cimport BOHR_MAGNETON, HC_EV_NM +from cherab.core.model.lineshape.doppler cimport doppler_shift, thermal_broadening +from cherab.core.model.lineshape.gaussian cimport add_gaussian_line + +cimport cython + + +DEF MULTIPLET_WAVELENGTH = 0 +DEF MULTIPLET_RATIO = 1 + +DEF PI_POLARISATION = 0 +DEF SIGMA_POLARISATION = 1 +DEF SIGMA_PLUS_POLARISATION = 1 +DEF SIGMA_MINUS_POLARISATION = -1 +DEF NO_POLARISATION = 2 + + +cdef class ZeemanLineShapeModel(LineShapeModel): + r""" + A base class for building Zeeman line shapes. + + :param Line line: The emission line object for this line shape. + :param float wavelength: The rest wavelength for this emission line. + :param Species target_species: The target plasma species that is emitting. + :param Plasma plasma: The emitting plasma object. + :param AtomicData atomic_data: The atomic data provider. + :param str polarisation: Leaves only :math:`\pi`-/:math:`\sigma`-polarised components: + "pi" - leave only :math:`\pi`-polarised components, + "sigma" - leave only :math:`\sigma`-polarised components, + "no" - leave all components (default). + :param Integrator1D integrator: Integrator1D instance to integrate the line shape + over the spectral bin. Default is None. + """ + + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, AtomicData atomic_data, + polarisation='no', Integrator1D integrator=None): + super().__init__(line, wavelength, target_species, plasma, atomic_data, integrator) + + self.polarisation = polarisation + + @property + def polarisation(self): + if self._polarisation == PI_POLARISATION: + return 'pi' + if self._polarisation == SIGMA_POLARISATION: + return 'sigma' + if self._polarisation == NO_POLARISATION: + return 'no' + + @polarisation.setter + def polarisation(self, value): + if value.lower() == 'pi': + self._polarisation = PI_POLARISATION + elif value.lower() == 'sigma': + self._polarisation = SIGMA_POLARISATION + elif value.lower() == 'no': + self._polarisation = NO_POLARISATION + else: + raise ValueError('Select between "pi", "sigma" or "no", {} is unsupported.'.format(value)) + + +cdef class ZeemanTriplet(ZeemanLineShapeModel): + r""" + Simple Doppler-Zeeman triplet (Paschen-Back effect). + + :param Line line: The emission line object for this line shape. + :param float wavelength: The rest wavelength for this emission line. + :param Species target_species: The target plasma species that is emitting. + :param Plasma plasma: The emitting plasma object. + :param AtomicData atomic_data: The atomic data provider. + :param str polarisation: Leaves only :math:`\pi`-/:math:`\sigma`-polarised components: + "pi" - leave central component, + "sigma" - leave side components, + "no" - all components (default). + """ + + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, AtomicData atomic_data, polarisation='no'): + + super().__init__(line, wavelength, target_species, plasma, atomic_data, polarisation) + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + @cython.cdivision(True) + cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): + + cdef double ts, sigma, shifted_wavelength, photon_energy, b_magn, component_radiance, cos_sqr, sin_sqr + cdef Vector3D ion_velocity, b_field + + ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) + if ts <= 0.0: + return spectrum + + ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) + + # calculate emission line central wavelength, doppler shifted along observation direction + shifted_wavelength = doppler_shift(self.wavelength, direction, ion_velocity) + + # calculate the line width + sigma = thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) + + # obtain magnetic field + b_field = self.plasma.get_b_field().evaluate(point.x, point.y, point.z) + b_magn = b_field.get_length() + + if b_magn == 0: + # no splitting if magnetic field strength is zero + if self._polarisation == NO_POLARISATION: + return add_gaussian_line(radiance, shifted_wavelength, sigma, spectrum) + + return add_gaussian_line(0.5 * radiance, shifted_wavelength, sigma, spectrum) + + # coefficients for intensities parallel and perpendicular to magnetic field + cos_sqr = (b_field.dot(direction.normalise()) / b_magn)**2 + sin_sqr = 1. - cos_sqr + + # adding pi component of the Zeeman triplet in case of NO_POLARISATION or PI_POLARISATION + if self._polarisation != SIGMA_POLARISATION: + component_radiance = 0.5 * sin_sqr * radiance + spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) + + # adding sigma +/- components of the Zeeman triplet in case of NO_POLARISATION or SIGMA_POLARISATION + if self._polarisation != PI_POLARISATION: + component_radiance = (0.25 * sin_sqr + 0.5 * cos_sqr) * radiance + + photon_energy = HC_EV_NM / self.wavelength + + shifted_wavelength = doppler_shift(HC_EV_NM / (photon_energy - BOHR_MAGNETON * b_magn), direction, ion_velocity) + spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) + + shifted_wavelength = doppler_shift(HC_EV_NM / (photon_energy + BOHR_MAGNETON * b_magn), direction, ion_velocity) + spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) + + return spectrum + + +cdef class ParametrisedZeemanTriplet(ZeemanLineShapeModel): + r""" + Parametrised Doppler-Zeeman triplet. It takes into account additional broadening due to + the line's fine structure without resolving the individual components of the fine + structure. The model is described with three parameters: :math:`\alpha`, + :math:`\beta` and :math:`\gamma`. + + The distance between :math:`\sigma^+` and :math:`\sigma^-` peaks: + :math:`\Delta \lambda_{\sigma} = \alpha B`, + where `B` is the magnetic field strength. + The ratio between Zeeman and thermal broadening line widths: + :math:`\frac{W_{Zeeman}}{W_{Doppler}} = \beta T^{\gamma}`, + where `T` is the species temperature in eV. + + For details see A. Blom and C. Jupén, Parametrisation of the Zeeman effect + for hydrogen-like spectra in high-temperature plasmas, + Plasma Phys. Control. Fusion 44 (2002) `1229-1241 + `_. + + :param Line line: The emission line object for this line shape. + :param float wavelength: The rest wavelength for this emission line. + :param Species target_species: The target plasma species that is emitting. + :param Plasma plasma: The emitting plasma object. + :param AtomicData atomic_data: The atomic data provider. + :param tuple line_parameters: Parameters of the model in the form (alpha, beta, gamma). + Default is None (will use `atomic_data.zeeman_triplet_parameters`). + :param str polarisation: Leaves only :math:`\pi`-/:math:`\sigma`-polarised components: + "pi" - leave central component, + "sigma" - leave side components, + "no" - all components (default). + """ + + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, AtomicData atomic_data, + tuple line_parameters=None, polarisation='no'): + + super().__init__(line, wavelength, target_species, plasma, atomic_data, polarisation) + + try: + alpha, beta, gamma = line_parameters or self.atomic_data.zeeman_triplet_parameters(line) + if alpha <= 0: + raise ValueError('Parameter alpha must be positive.') + if beta < 0: + raise ValueError('Parameter beta must be non-negative.') + self._alpha = alpha + self._beta = beta + self._gamma = gamma + + except KeyError: + raise ValueError('Data for {} is not available.'.format(self.line)) + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + @cython.cdivision(True) + cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): + + cdef double ts, sigma, shifted_wavelength, b_magn, component_radiance, cos_sqr, sin_sqr + cdef Vector3D ion_velocity, b_field + + ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) + if ts <= 0.0: + return spectrum + + ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) + + # calculate emission line central wavelength, doppler shifted along observation direction + shifted_wavelength = doppler_shift(self.wavelength, direction, ion_velocity) + + # calculate the line width + sigma = thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) + + # fine structure broadening correction + sigma *= sqrt(1. + self._beta * self._beta * ts**(2. * self._gamma)) + + # obtain magnetic field + b_field = self.plasma.get_b_field().evaluate(point.x, point.y, point.z) + b_magn = b_field.get_length() + + if b_magn == 0: + # no splitting if magnetic filed strength is zero + if self._polarisation == NO_POLARISATION: + return add_gaussian_line(radiance, shifted_wavelength, sigma, spectrum) + + return add_gaussian_line(0.5 * radiance, shifted_wavelength, sigma, spectrum) + + # coefficients for intensities parallel and perpendicular to magnetic field + cos_sqr = (b_field.dot(direction.normalise()) / b_magn)**2 + sin_sqr = 1. - cos_sqr + + # adding pi component of the Zeeman triplet in case of NO_POLARISATION or PI_POLARISATION + if self._polarisation != SIGMA_POLARISATION: + component_radiance = 0.5 * sin_sqr * radiance + spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) + + # adding sigma +/- components of the Zeeman triplet in case of NO_POLARISATION or SIGMA_POLARISATION + if self._polarisation != PI_POLARISATION: + component_radiance = (0.25 * sin_sqr + 0.5 * cos_sqr) * radiance + + shifted_wavelength = doppler_shift(self.wavelength + 0.5 * self._alpha * b_magn, direction, ion_velocity) + spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) + shifted_wavelength = doppler_shift(self.wavelength - 0.5 * self._alpha * b_magn, direction, ion_velocity) + spectrum = add_gaussian_line(component_radiance, shifted_wavelength, sigma, spectrum) + + return spectrum + + +cdef class ZeemanMultiplet(ZeemanLineShapeModel): + r""" + Doppler-Zeeman Multiplet. + + The lineshape radiance is calculated from a base PEC rate that is unresolved. This + radiance is then divided over a number of components as specified in the ``zeeman_structure`` + argument. The ``zeeman_structure`` specifies wavelengths and ratios of + :math:`\pi`-/:math:`\sigma`-polarised components as functions of the magnetic field strength. + These functions can be obtained using the output of the ADAS603 routines. + + :param Line line: The emission line object for the base rate radiance calculation. + :param float wavelength: The rest wavelength of the base emission line. + :param Species target_species: The target plasma species that is emitting. + :param Plasma plasma: The emitting plasma object. + :param AtomicData atomic_data: The atomic data provider. + :param zeeman_structure: A ``ZeemanStructure`` object that provides wavelengths and ratios + of :math:`\pi`-/:math:`\sigma^{+}`-/:math:`\sigma^{-}`-polarised + components for any given magnetic field strength. + Default is None (will use atomic_data.zeeman_structure). + :param str polarisation: Leaves only :math:`\pi`-/:math:`\sigma`-polarised components: + "pi" - leave only :math:`\pi`-polarised components, + "sigma" - leave only :math:`\sigma`-polarised components, + "no" - leave all components (default). + + """ + + def __init__(self, Line line, double wavelength, Species target_species, Plasma plasma, AtomicData atomic_data, + ZeemanStructure zeeman_structure=None, polarisation='no'): + + super().__init__(line, wavelength, target_species, plasma, atomic_data, polarisation) + + self._zeeman_structure = zeeman_structure or self.atomic_data.zeeman_structure(line) + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + @cython.cdivision(True) + cpdef Spectrum add_line(self, double radiance, Point3D point, Vector3D direction, Spectrum spectrum): + + cdef int i + cdef double ts, sigma, shifted_wavelength, component_radiance, b_magn, cos_sqr, sin_sqr + cdef Vector3D ion_velocity, b_field + cdef double[:, :] multiplet_mv + + ts = self.target_species.distribution.effective_temperature(point.x, point.y, point.z) + if ts <= 0.0: + return spectrum + + ion_velocity = self.target_species.distribution.bulk_velocity(point.x, point.y, point.z) + + # calculate the line width + sigma = thermal_broadening(self.wavelength, ts, self.line.element.atomic_weight) + + # obtain magnetic field + b_field = self.plasma.get_b_field().evaluate(point.x, point.y, point.z) + b_magn = b_field.get_length() + + if b_magn == 0: + # no splitting if magnetic filed strength is zero + shifted_wavelength = doppler_shift(self.wavelength, direction, ion_velocity) + if self._polarisation == NO_POLARISATION: + return add_gaussian_line(radiance, shifted_wavelength, sigma, spectrum) + + return add_gaussian_line(0.5 * radiance, shifted_wavelength, sigma, spectrum) + + # coefficients for intensities parallel and perpendicular to magnetic field + cos_sqr = (b_field.dot(direction.normalise()) / b_magn)**2 + sin_sqr = 1. - cos_sqr + + # adding pi components of the Zeeman multiplet in case of NO_POLARISATION or PI_POLARISATION + if self._polarisation != SIGMA_POLARISATION: + component_radiance = 0.5 * sin_sqr * radiance + multiplet_mv = self._zeeman_structure.evaluate(b_magn, PI_POLARISATION) + + for i in range(multiplet_mv.shape[1]): + shifted_wavelength = doppler_shift(multiplet_mv[MULTIPLET_WAVELENGTH, i], direction, ion_velocity) + spectrum = add_gaussian_line(component_radiance * multiplet_mv[MULTIPLET_RATIO, i], shifted_wavelength, sigma, spectrum) + + # adding sigma components of the Zeeman multiplet in case of NO_POLARISATION or SIGMA_POLARISATION + if self._polarisation != PI_POLARISATION: + component_radiance = (0.25 * sin_sqr + 0.5 * cos_sqr) * radiance + + multiplet_mv = self._zeeman_structure.evaluate(b_magn, SIGMA_PLUS_POLARISATION) + + for i in range(multiplet_mv.shape[1]): + shifted_wavelength = doppler_shift(multiplet_mv[MULTIPLET_WAVELENGTH, i], direction, ion_velocity) + spectrum = add_gaussian_line(component_radiance * multiplet_mv[MULTIPLET_RATIO, i], shifted_wavelength, sigma, spectrum) + + multiplet_mv = self._zeeman_structure.evaluate(b_magn, SIGMA_MINUS_POLARISATION) + + for i in range(multiplet_mv.shape[1]): + shifted_wavelength = doppler_shift(multiplet_mv[MULTIPLET_WAVELENGTH, i], direction, ion_velocity) + spectrum = add_gaussian_line(component_radiance * multiplet_mv[MULTIPLET_RATIO, i], shifted_wavelength, sigma, spectrum) + + return spectrum diff --git a/cherab/core/model/plasma/impact_excitation.pyx b/cherab/core/model/plasma/impact_excitation.pyx index 8a5d993c..b26336f1 100644 --- a/cherab/core/model/plasma/impact_excitation.pyx +++ b/cherab/core/model/plasma/impact_excitation.pyx @@ -1,6 +1,8 @@ -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -18,7 +20,7 @@ from raysect.optical cimport Spectrum, Point3D, Vector3D from cherab.core cimport Plasma, AtomicData -from cherab.core.model.lineshape cimport GaussianLine, LineShapeModel +from cherab.core.model.lineshape cimport GaussianLine from cherab.core.utility.constants cimport RECIP_4_PI @@ -123,7 +125,7 @@ cdef class ExcitationLine(PlasmaModel): # instance line shape renderer self._lineshape = self._lineshape_class(self._line, self._wavelength, self._target_species, self._plasma, - *self._lineshape_args, **self._lineshape_kwargs) + self._atomic_data, *self._lineshape_args, **self._lineshape_kwargs) def _change(self): diff --git a/cherab/core/model/plasma/recombination.pyx b/cherab/core/model/plasma/recombination.pyx index 004205a0..a33009d0 100644 --- a/cherab/core/model/plasma/recombination.pyx +++ b/cherab/core/model/plasma/recombination.pyx @@ -1,6 +1,8 @@ -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# cython: language_level=3 + +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -126,7 +128,7 @@ cdef class RecombinationLine(PlasmaModel): # instance line shape renderer self._lineshape = self._lineshape_class(self._line, self._wavelength, self._target_species, self._plasma, - *self._lineshape_args, **self._lineshape_kwargs) + self._atomic_data, *self._lineshape_args, **self._lineshape_kwargs) def _change(self): diff --git a/cherab/core/model/plasma/thermal_cx.pyx b/cherab/core/model/plasma/thermal_cx.pyx index 1325919c..88af9ae8 100644 --- a/cherab/core/model/plasma/thermal_cx.pyx +++ b/cherab/core/model/plasma/thermal_cx.pyx @@ -152,7 +152,7 @@ cdef class ThermalCXLine(PlasmaModel): # instance line shape renderer self._lineshape = self._lineshape_class(self._line, self._wavelength, self._target_species, self._plasma, - *self._lineshape_args, **self._lineshape_kwargs) + self._atomic_data, *self._lineshape_args, **self._lineshape_kwargs) def _change(self): diff --git a/cherab/core/tests/test_beamcxline.py b/cherab/core/tests/test_beamcxline.py index d9a0a4f4..a93cf946 100644 --- a/cherab/core/tests/test_beamcxline.py +++ b/cherab/core/tests/test_beamcxline.py @@ -128,7 +128,7 @@ def test_default_lineshape(self): target_species = self.plasma.composition.get(line.element, line.charge + 1) wavelength = self.atomic_data.wavelength(line.element, line.charge, line.transition) - gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma) + gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma, self.atomic_data) spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) spectrum = gaussian_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) @@ -160,7 +160,7 @@ def test_custom_lineshape(self): target_species = self.plasma.composition.get(line.element, line.charge + 1) wavelength = self.atomic_data.wavelength(line.element, line.charge, line.transition) - zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma) + zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma, self.atomic_data) spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) spectrum = zeeman_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) diff --git a/cherab/core/tests/test_line_emission.py b/cherab/core/tests/test_line_emission.py index b55af871..bf40aac0 100644 --- a/cherab/core/tests/test_line_emission.py +++ b/cherab/core/tests/test_line_emission.py @@ -125,7 +125,7 @@ def test_default_lineshape(self): ni = target_species.distribution.density(0.5, 0, 0) radiance = 0.25 / np.pi * rate * ni * ne * self.slab_length - gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma) + gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma, self.atomic_data) spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) spectrum = gaussian_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) @@ -154,7 +154,7 @@ def test_custom_lineshape(self): ni = target_species.distribution.density(0.5, 0, 0) radiance = 0.25 / np.pi * rate * ni * ne * self.slab_length - zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma) + zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma, self.atomic_data) spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) spectrum = zeeman_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) @@ -200,7 +200,7 @@ def test_default_lineshape(self): ni = target_species.distribution.density(0.5, 0, 0) radiance = 0.25 / np.pi * rate * ni * ne * self.slab_length - gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma) + gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma, self.atomic_data) spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) spectrum = gaussian_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) @@ -229,7 +229,7 @@ def test_custom_lineshape(self): ni = target_species.distribution.density(0.5, 0, 0) radiance = 0.25 / np.pi * rate * ni * ne * self.slab_length - zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma) + zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma, self.atomic_data) spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) spectrum = zeeman_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) @@ -279,7 +279,7 @@ def test_default_lineshape(self): receiver_density = target_species.distribution.density(0.5, 0, 0) radiance = 0.25 / np.pi * rate * receiver_density * donor_density * self.slab_length - gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma) + gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma, self.atomic_data) spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) spectrum = gaussian_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) @@ -311,7 +311,7 @@ def test_custom_lineshape(self): receiver_density = target_species.distribution.density(0.5, 0, 0) radiance = 0.25 / np.pi * rate * receiver_density * donor_density * self.slab_length - zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma) + zeeman_line = ZeemanTriplet(line, wavelength, target_species, self.plasma, self.atomic_data) spectrum = Spectrum(ray.min_wavelength, ray.max_wavelength, ray.bins) spectrum = zeeman_line.add_line(radiance, Point3D(0.5, 0, 0), direction, spectrum) diff --git a/cherab/core/tests/test_lineshapes.py b/cherab/core/tests/test_lineshapes.py index 97a32f1f..b6827256 100644 --- a/cherab/core/tests/test_lineshapes.py +++ b/cherab/core/tests/test_lineshapes.py @@ -1,6 +1,6 @@ -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -26,11 +26,12 @@ from raysect.core.math.function.float import Arg1D, Constant1D from raysect.optical import Spectrum -from cherab.core import Line +from cherab.core import Beam, Line, AtomicData from cherab.core.math.integrators import GaussianQuadrature from cherab.core.atomic import deuterium, nitrogen, ZeemanStructure from cherab.tools.plasmas.slab import build_constant_slab_plasma from cherab.core.model import GaussianLine, MultipletLineShape, StarkBroadenedLine, ZeemanTriplet, ParametrisedZeemanTriplet, ZeemanMultiplet +from cherab.core.model import BeamEmissionMultiplet ATOMIC_MASS = 1.66053906660e-27 @@ -51,13 +52,19 @@ def setUp(self): electron_temperature=20., plasma_species=plasma_species, b_field=Vector3D(0, 5., 0)) + self.atomic_data = AtomicData() + self.beam = Beam() + self.beam.plasma = self.plasma + self.beam.energy = 60000 + self.beam.temperature = 10 + self.beam.element = deuterium def test_gaussian_line(self): # setting up a line shape model line = Line(deuterium, 0, (3, 2)) # D-alpha line target_species = self.plasma.composition.get(line.element, line.charge) wavelength = 656.104 - gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma) + gaussian_line = GaussianLine(line, wavelength, target_species, self.plasma, self.atomic_data) # spectrum parameters min_wavelength = wavelength - 0.5 @@ -92,7 +99,7 @@ def test_multiplet_line_shape(self): target_species = self.plasma.composition.get(line.element, line.charge) multiplet = [[403.509, 404.132, 404.354, 404.479, 405.692], [0.205, 0.562, 0.175, 0.029, 0.029]] wavelength = 404.21 - multiplet_line = MultipletLineShape(line, wavelength, target_species, self.plasma, multiplet) + multiplet_line = MultipletLineShape(line, wavelength, target_species, self.plasma, self.atomic_data, multiplet) # spectrum parameters min_wavelength = min(multiplet[0]) - 0.5 @@ -128,7 +135,7 @@ def test_zeeman_triplet(self): line = Line(deuterium, 0, (3, 2)) # D-alpha line target_species = self.plasma.composition.get(line.element, line.charge) wavelength = 656.104 - triplet = ZeemanTriplet(line, wavelength, target_species, self.plasma) + triplet = ZeemanTriplet(line, wavelength, target_species, self.plasma, self.atomic_data) # spectrum parameters min_wavelength = wavelength - 0.5 @@ -179,7 +186,7 @@ def test_parametrised_zeeman_triplet(self): line = Line(deuterium, 0, (3, 2)) # D-alpha line target_species = self.plasma.composition.get(line.element, line.charge) wavelength = 656.104 - triplet = ParametrisedZeemanTriplet(line, wavelength, target_species, self.plasma) + triplet = ParametrisedZeemanTriplet(line, wavelength, target_species, self.plasma, self.atomic_data) # spectrum parameters min_wavelength = wavelength - 0.5 @@ -197,7 +204,7 @@ def test_parametrised_zeeman_triplet(self): spectrum[pol] = triplet.add_line(radiance, point, direction, spectrum[pol]) # validating - alpha, beta, gamma = triplet.LINE_PARAMETERS_DEFAULT[line] + alpha, beta, gamma = self.atomic_data.zeeman_triplet_parameters(line) temperature = target_species.distribution.effective_temperature(point.x, point.y, point.z) velocity = target_species.distribution.bulk_velocity(point.x, point.y, point.z) sigma = np.sqrt(temperature * ELEMENTARY_CHARGE / (line.element.atomic_weight * ATOMIC_MASS)) * wavelength / SPEED_OF_LIGHT @@ -238,7 +245,7 @@ def test_zeeman_multiplet(self): sigma_minus_components = [(HC_EV_NM / (photon_energy + BOHR_MAGNETON * Arg1D()), Constant1D(0.5))] zeeman_structure = ZeemanStructure(pi_components, sigma_plus_components, sigma_minus_components) - multiplet = ZeemanMultiplet(line, wavelength, target_species, self.plasma, zeeman_structure) + multiplet = ZeemanMultiplet(line, wavelength, target_species, self.plasma, self.atomic_data, zeeman_structure) # spectrum parameters min_wavelength = wavelength - 0.5 @@ -290,7 +297,7 @@ def test_stark_broadened_line(self): target_species = self.plasma.composition.get(line.element, line.charge) wavelength = 656.104 integrator = GaussianQuadrature(relative_tolerance=1.e-5) - stark_line = StarkBroadenedLine(line, wavelength, target_species, self.plasma, integrator=integrator) + stark_line = StarkBroadenedLine(line, wavelength, target_species, self.plasma, self.atomic_data, integrator=integrator) # spectrum parameters min_wavelength = wavelength - 0.2 @@ -305,25 +312,160 @@ def test_stark_broadened_line(self): spectrum = stark_line.add_line(radiance, point, direction, spectrum) # validating - cij, aij, bij = stark_line.STARK_MODEL_COEFFICIENTS_DEFAULT[line] + velocity = target_species.distribution.bulk_velocity(point.x, point.y, point.z) + doppler_factor = (1 + velocity.dot(direction.normalise()) / SPEED_OF_LIGHT) + + b_field = self.plasma.b_field(point.x, point.y, point.z) + b_magn = b_field.length + photon_energy = HC_EV_NM / wavelength + wl_sigma_plus = HC_EV_NM / (photon_energy - BOHR_MAGNETON * b_magn) + wl_sigma_minus = HC_EV_NM / (photon_energy + BOHR_MAGNETON * b_magn) + cos_sqr = (b_field.dot(direction.normalise()) / b_magn)**2 + sin_sqr = 1. - cos_sqr + + # Gaussian parameters + temperature = target_species.distribution.effective_temperature(point.x, point.y, point.z) + sigma = np.sqrt(temperature * ELEMENTARY_CHARGE / (line.element.atomic_weight * ATOMIC_MASS)) * wavelength / SPEED_OF_LIGHT + fwhm_gauss = 2 * np.sqrt(2 * np.log(2)) * sigma + + # Lorentzian parameters + cij, aij, bij = self.atomic_data.stark_model_coefficients(line) ne = self.plasma.electron_distribution.density(point.x, point.y, point.z) te = self.plasma.electron_distribution.effective_temperature(point.x, point.y, point.z) - lambda_1_2 = cij * ne**aij / (te**bij) + fwhm_lorentz = cij * ne**aij / (te**bij) + + # Total FWHM + if fwhm_gauss <= fwhm_lorentz: + fwhm_poly_coeff = [1., 0, 0.57575, 0.37902, -0.42519, -0.31525, 0.31718] + fwhm_ratio = fwhm_gauss / fwhm_lorentz + fwhm_full = fwhm_lorentz * np.poly1d(fwhm_poly_coeff[::-1])(fwhm_ratio) + else: + fwhm_poly_coeff = [1., 0.15882, 1.04388, -1.38281, 0.46251, 0.82325, -0.58026] + fwhm_ratio = fwhm_lorentz / fwhm_gauss + fwhm_full = fwhm_gauss * np.poly1d(fwhm_poly_coeff[::-1])(fwhm_ratio) + + wavelengths, delta = np.linspace(min_wavelength, max_wavelength, bins + 1, retstep=True) + + # Gaussian part + temp = 2 * np.sqrt(np.log(2)) / fwhm_full + erfs = erf((wavelengths - wavelength * doppler_factor) * temp) + gaussian = 0.25 * sin_sqr * (erfs[1:] - erfs[:-1]) / delta + erfs = erf((wavelengths - wl_sigma_plus * doppler_factor) * temp) + gaussian += 0.5 * (0.25 * sin_sqr + 0.5 * cos_sqr) * (erfs[1:] - erfs[:-1]) / delta + erfs = erf((wavelengths - wl_sigma_minus * doppler_factor) * temp) + gaussian += 0.5 * (0.25 * sin_sqr + 0.5 * cos_sqr) * (erfs[1:] - erfs[:-1]) / delta + # Lorentzian part lorenzian_cutoff_gamma = 50 stark_norm_coeff = 4 * lorenzian_cutoff_gamma * hyp2f1(0.4, 1, 1.4, -(2 * lorenzian_cutoff_gamma)**2.5) - norm = (0.5 * lambda_1_2)**1.5 / stark_norm_coeff + norm = (0.5 * fwhm_full)**1.5 / stark_norm_coeff - wavelengths, delta = np.linspace(min_wavelength, max_wavelength, bins + 1, retstep=True) + def stark_lineshape_pi(x): + return norm / ((np.abs(x - wavelength * doppler_factor))**2.5 + (0.5 * fwhm_full)**2.5) - def stark_lineshape(x): - return norm / ((np.abs(x - wavelength))**2.5 + (0.5 * lambda_1_2)**2.5) + def stark_lineshape_sigma_plus(x): + return norm / ((np.abs(x - wl_sigma_plus * doppler_factor))**2.5 + (0.5 * fwhm_full)**2.5) + + def stark_lineshape_sigma_minus(x): + return norm / ((np.abs(x - wl_sigma_minus * doppler_factor))**2.5 + (0.5 * fwhm_full)**2.5) + + weight_poly_coeff = [5.14820e-04, 1.38821e+00, -9.60424e-02, -3.83995e-02, -7.40042e-03, -5.47626e-04] + lorentz_weight = np.exp(np.poly1d(weight_poly_coeff[::-1])(np.log(fwhm_lorentz / fwhm_full))) for i in range(bins): - stark_bin = quadrature(stark_lineshape, wavelengths[i], wavelengths[i + 1], rtol=integrator.relative_tolerance)[0] / delta - self.assertAlmostEqual(stark_bin, spectrum.samples[i], delta=1e-9, + lorentz_bin = 0.5 * sin_sqr * quadrature(stark_lineshape_pi, wavelengths[i], wavelengths[i + 1], + rtol=integrator.relative_tolerance)[0] / delta + lorentz_bin += (0.25 * sin_sqr + 0.5 * cos_sqr) * quadrature(stark_lineshape_sigma_plus, wavelengths[i], wavelengths[i + 1], + rtol=integrator.relative_tolerance)[0] / delta + lorentz_bin += (0.25 * sin_sqr + 0.5 * cos_sqr) * quadrature(stark_lineshape_sigma_minus, wavelengths[i], wavelengths[i + 1], + rtol=integrator.relative_tolerance)[0] / delta + ref_value = lorentz_bin * lorentz_weight + gaussian[i] * (1. - lorentz_weight) + self.assertAlmostEqual(ref_value, spectrum.samples[i], delta=1e-9, msg='StarkBroadenedLine.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) + def test_beam_emission_multiplet(self): + # Test MSE line shape + # setting up a line shape model + line = Line(deuterium, 0, (3, 2)) # D-alpha line + wavelength = 656.104 + sigma_to_pi = 0.56 + sigma1_to_sigma0 = 0.7060001671878492 + pi2_to_pi3 = 0.3140003593919741 + pi4_to_pi3 = 0.7279994935840365 + mse_line = BeamEmissionMultiplet(line, wavelength, self.beam, self.atomic_data, + sigma_to_pi, sigma1_to_sigma0, pi2_to_pi3, pi4_to_pi3) + + # spectrum parameters + min_wavelength = wavelength - 3 + max_wavelength = wavelength + 3 + bins = 512 + point = Point3D(0.5, 0.5, 0.5) + direction = Vector3D(-1, 1, 0) / np.sqrt(2) + beam_direction = self.beam.direction(point.x, point.y, point.z) + + # obtaining spectrum + radiance = 1.0 + spectrum = Spectrum(min_wavelength, max_wavelength, bins) + spectrum = mse_line.add_line(radiance, point, point, beam_direction, direction, spectrum) + + # validating + + # calculate Stark splitting + b_field = self.plasma.b_field(point.x, point.y, point.z) + beam_velocity = beam_direction.normalise() * np.sqrt(2 * self.beam.energy * ELEMENTARY_CHARGE / ATOMIC_MASS) + e_field = beam_velocity.cross(b_field).length + STARK_SPLITTING_FACTOR = 2.77e-8 + stark_split = np.abs(STARK_SPLITTING_FACTOR * e_field) + + # calculate emission line central wavelength, doppler shifted along observation direction + central_wavelength = wavelength * (1 + beam_velocity.dot(direction.normalise()) / SPEED_OF_LIGHT) + + # calculate doppler broadening + beam_ion_mass = self.beam.element.atomic_weight + beam_temperature = self.beam.temperature + sigma = np.sqrt(beam_temperature * ELEMENTARY_CHARGE / (beam_ion_mass * ATOMIC_MASS)) * wavelength / SPEED_OF_LIGHT + temp = 1. / (np.sqrt(2.) * sigma) + + # calculate relative intensities of sigma and pi lines + d = 1 / (1 + sigma_to_pi) + intensity_sig = sigma_to_pi * d * radiance + intensity_pi = 0.5 * d * radiance + + wavelengths, delta = np.linspace(min_wavelength, max_wavelength, bins + 1, retstep=True) + + # add Sigma lines to output + intensity_s0 = 1 / (sigma1_to_sigma0 + 1) + intensity_s1 = 0.5 * sigma1_to_sigma0 * intensity_s0 + + erfs = erf((wavelengths - central_wavelength) * temp) + test_spectrum = 0.5 * intensity_sig * intensity_s0 * (erfs[1:] - erfs[:-1]) / delta + erfs = erf((wavelengths - central_wavelength - stark_split) * temp) + test_spectrum += 0.5 * intensity_sig * intensity_s1 * (erfs[1:] - erfs[:-1]) / delta + erfs = erf((wavelengths - central_wavelength + stark_split) * temp) + test_spectrum += 0.5 * intensity_sig * intensity_s1 * (erfs[1:] - erfs[:-1]) / delta + + # add Pi lines to output + intensity_pi3 = 1 / (1 + pi2_to_pi3 + pi4_to_pi3) + intensity_pi2 = pi2_to_pi3 * intensity_pi3 + intensity_pi4 = pi4_to_pi3 * intensity_pi3 + + erfs = erf((wavelengths - central_wavelength - 2 * stark_split) * temp) + test_spectrum += 0.5 * intensity_pi * intensity_pi2 * (erfs[1:] - erfs[:-1]) / delta + erfs = erf((wavelengths - central_wavelength + 2 * stark_split) * temp) + test_spectrum += 0.5 * intensity_pi * intensity_pi2 * (erfs[1:] - erfs[:-1]) / delta + erfs = erf((wavelengths - central_wavelength - 3 * stark_split) * temp) + test_spectrum += 0.5 * intensity_pi * intensity_pi3 * (erfs[1:] - erfs[:-1]) / delta + erfs = erf((wavelengths - central_wavelength + 3 * stark_split) * temp) + test_spectrum += 0.5 * intensity_pi * intensity_pi3 * (erfs[1:] - erfs[:-1]) / delta + erfs = erf((wavelengths - central_wavelength - 4 * stark_split) * temp) + test_spectrum += 0.5 * intensity_pi * intensity_pi4 * (erfs[1:] - erfs[:-1]) / delta + erfs = erf((wavelengths - central_wavelength + 4 * stark_split) * temp) + test_spectrum += 0.5 * intensity_pi * intensity_pi4 * (erfs[1:] - erfs[:-1]) / delta + + for i in range(bins): + self.assertAlmostEqual(test_spectrum[i], spectrum.samples[i], delta=1e-10, + msg='BeamEmissionMultiplet.add_line() method gives a wrong value at {} nm.'.format(wavelengths[i])) + if __name__ == '__main__': unittest.main() diff --git a/cherab/core/utility/constants.pxd b/cherab/core/utility/constants.pxd index 4c59b4a8..5ec0ca4a 100644 --- a/cherab/core/utility/constants.pxd +++ b/cherab/core/utility/constants.pxd @@ -27,7 +27,9 @@ cdef: double ELEMENTARY_CHARGE double SPEED_OF_LIGHT double PLANCK_CONSTANT + double HC_EV_NM double ELECTRON_CLASSICAL_RADIUS double ELECTRON_REST_MASS double RYDBERG_CONSTANT_EV double VACUUM_PERMITTIVITY + double BOHR_MAGNETON diff --git a/cherab/core/utility/constants.pyx b/cherab/core/utility/constants.pyx index d36b99f9..423d15e8 100644 --- a/cherab/core/utility/constants.pyx +++ b/cherab/core/utility/constants.pyx @@ -29,7 +29,9 @@ cdef: double ELEMENTARY_CHARGE = 1.602176634e-19 double SPEED_OF_LIGHT = 299792458.0 double PLANCK_CONSTANT = 6.62607015e-34 + double HC_EV_NM = 1239.8419738620933 # (Planck constant in eV s) x (speed of light in nm/s) double ELECTRON_CLASSICAL_RADIUS = 2.8179403262e-15 double ELECTRON_REST_MASS = 9.1093837015e-31 double RYDBERG_CONSTANT_EV = 13.605693122994 double VACUUM_PERMITTIVITY = 8.8541878128e-12 + double BOHR_MAGNETON = 5.78838180123e-5 # in eV/T diff --git a/demos/emission_models/stark_broadening.py b/demos/emission_models/stark_broadening.py index 090e929a..b3f2cf95 100755 --- a/demos/emission_models/stark_broadening.py +++ b/demos/emission_models/stark_broadening.py @@ -1,6 +1,6 @@ -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -34,7 +34,7 @@ # tunables -ion_density = 1e20 +ion_density = 2e20 sigma = 0.25 # setup scenegraph @@ -99,4 +99,5 @@ plt.xlabel('Wavelength (nm)') plt.ylabel('Radiance (W/m^2/str/nm)') plt.title('Observed Spectrum') +plt.yscale('log') plt.show() diff --git a/demos/emission_models/stark_zeeman.py b/demos/emission_models/stark_zeeman.py new file mode 100755 index 00000000..6dac025d --- /dev/null +++ b/demos/emission_models/stark_zeeman.py @@ -0,0 +1,95 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +""" +Calculates Balmer-alpha and Paschen-beta Stark-Zeeman spectral lines. + +Compare the figures with Figure 2 in B.A. Lomanowski at al, Nucl. Fusion 55 (2015) 123028, +https://iopscience.iop.org/article/10.1088/0029-5515/55/12/123028 +""" + + +# External imports +import matplotlib.pyplot as plt +from raysect.optical import World, Vector3D, Point3D, Ray + +# Cherab imports +from cherab.core import Line +from cherab.core.atomic.elements import deuterium +from cherab.core.model import ExcitationLine, StarkBroadenedLine +from cherab.openadas import OpenADAS +from cherab.tools.plasmas.slab import build_constant_slab_plasma + + +# tunables +ion_density = 2e20 +sigma = 0.25 + +# setup scenegraph +world = World() + +# create atomic data source +adas = OpenADAS(permit_extrapolation=True) + +# setup the Balmer and Paschen lines +balmer_alpha = Line(deuterium, 0, (3, 2)) +paschen_beta = Line(deuterium, 0, (5, 3)) + +# setup plasma for Balmer-alpha line +plasma_species = [(deuterium, 0, 1.e20, 0.3, Vector3D(0, 0, 0))] +plasma = build_constant_slab_plasma(length=1, width=1, height=1, electron_density=5e20, electron_temperature=1., + plasma_species=plasma_species, b_field=Vector3D(0, 3., 0), parent=world) +plasma.atomic_data = adas + +# add Balmer-alpha line to the plasma +plasma.models = [ExcitationLine(balmer_alpha, lineshape=StarkBroadenedLine)] + +# Ray-trace perpendicular to magnetic field and plot the results +r = Ray(origin=Point3D(-5, 0, 0), direction=Vector3D(1, 0, 0), + min_wavelength=655.9, max_wavelength=656.3, bins=200) +s = r.trace(world) + +plt.figure() +plt.plot(s.wavelengths, s.samples, ls='--', color='C3') +plt.xlabel('Wavelength (nm)') +plt.ylabel('Radiance (W/m^2/str/nm)') +plt.title('Balmer-alpha spectrum, Ne = 5e20 m-3, Te = 1 eV, B = 3 T') + +plasma.parent = None + +# setup plasma for Paschen-beta line +plasma_species = [(deuterium, 0, 1.e20, 0.3, Vector3D(0, 0, 0))] +plasma = build_constant_slab_plasma(length=1, width=1, height=1, electron_density=1e20, electron_temperature=1., + plasma_species=plasma_species, b_field=Vector3D(0, 3., 0), parent=world) +plasma.atomic_data = adas + +# add Paschen-beta line to the plasma +plasma.models = [ExcitationLine(paschen_beta, lineshape=StarkBroadenedLine)] + +# Ray-trace perpendicular to magnetic field and plot the results +r = Ray(origin=Point3D(-5, 0, 0), direction=Vector3D(1, 0, 0), + min_wavelength=1280.3, max_wavelength=1282.6, bins=200) +s = r.trace(world) + +plt.figure() +plt.plot(s.wavelengths, s.samples, ls='--', color='C3') +plt.xlabel('Wavelength (nm)') +plt.ylabel('Radiance (W/m^2/str/nm)') +plt.title('Paschen-beta spectrum, Ne = 1e20 m-3, Te = 1 eV, B = 3 T') + +plt.show() diff --git a/docs/source/demonstrations/demonstrations.rst b/docs/source/demonstrations/demonstrations.rst index 477f6620..0ee82ab3 100644 --- a/docs/source/demonstrations/demonstrations.rst +++ b/docs/source/demonstrations/demonstrations.rst @@ -153,16 +153,21 @@ Passive Spectroscopy - .. image:: ./passive_spectroscopy/multiplet_spectrum.png :height: 150px :width: 150px - * - :ref:`Stark Broadened Lines ` - - Specifying a Stark broadened lineshape instead of Doppler broadening. - - .. image:: ./passive_spectroscopy/stark_line_zoomed.png - :height: 150px - :width: 150px * - :ref:`Zeeman Spectroscopy ` - Specifying a Zeeman triplet or multiplet line shapes. - .. image:: ./passive_spectroscopy/zeeman_spectrum_45deg.png :height: 150px - :width: 150px + :width: 150px + * - :ref:`Stark Broadened Lines ` + - Specifying a Stark broadened lineshape. + - .. image:: ./passive_spectroscopy/stark_spectrum.png + :height: 150px + :width: 150px + * - :ref:`Stark-Zeeman Lines ` + - Modelling Stark-Zeeman lineshapes. + - .. image:: ./passive_spectroscopy/stark_zeeman_balmer_alpha.png + :height: 150px + :width: 150px Bolometry diff --git a/docs/source/demonstrations/passive_spectroscopy/stark_broadening.rst b/docs/source/demonstrations/passive_spectroscopy/stark_broadening.rst index 9a2f20e8..0830b7e1 100644 --- a/docs/source/demonstrations/passive_spectroscopy/stark_broadening.rst +++ b/docs/source/demonstrations/passive_spectroscopy/stark_broadening.rst @@ -14,9 +14,15 @@ electric field due to the presence of neighbouring ions. In normal tokamak operations this effect is negligible, except in the divertor region. It is possible to override the default doppler broadened line shape by -specifying the StarkBroadenedLine() lineshape class. In this example -we can see two stark broadened balmer series lines surrounding a -Nitrogen II multiplet feature. +specifying the StarkBroadenedLine() lineshape class. +This class suppors Balmer and Paschen series and is based on +the Stark-Doppler-Zeeman line shape model from B. Lomanowski, et al. +"Inferring divertor plasma properties from hydrogen Balmer +and Paschen series spectroscopy in JET-ILW." Nuclear Fusion 55.12 (2015) +`123028 `_. +In this example we can see two stark broadened balmer series lines surrounding a +Nitrogen II multiplet feature. The logarithmic scale is chosen to illustrate +the power-law decay of the Stark-broadened line wings. .. literalinclude:: ../../../../demos/emission_models/stark_broadening.py @@ -25,11 +31,4 @@ Nitrogen II multiplet feature. :width: 450px **Caption:** The observed spectrum with two stark broadened balmer lines - (397nm and 410nm) surrounding a NII multiplet feature. - -.. figure:: stark_line_zoomed.png - :align: center - :width: 450px - - **Caption:** A zoomed in view of the 410nm Balmer line revealing the - characteristic stark broadened lineshape. + (397nm and 410nm) surrounding a NII multiplet feature in the logarithmic scale. diff --git a/docs/source/demonstrations/passive_spectroscopy/stark_line_zoomed.png b/docs/source/demonstrations/passive_spectroscopy/stark_line_zoomed.png deleted file mode 100644 index 2a6ee04b..00000000 Binary files a/docs/source/demonstrations/passive_spectroscopy/stark_line_zoomed.png and /dev/null differ diff --git a/docs/source/demonstrations/passive_spectroscopy/stark_spectrum.png b/docs/source/demonstrations/passive_spectroscopy/stark_spectrum.png index 53230968..04d90589 100644 Binary files a/docs/source/demonstrations/passive_spectroscopy/stark_spectrum.png and b/docs/source/demonstrations/passive_spectroscopy/stark_spectrum.png differ diff --git a/docs/source/demonstrations/passive_spectroscopy/stark_zeeman.rst b/docs/source/demonstrations/passive_spectroscopy/stark_zeeman.rst new file mode 100644 index 00000000..eb07750e --- /dev/null +++ b/docs/source/demonstrations/passive_spectroscopy/stark_zeeman.rst @@ -0,0 +1,39 @@ + +.. _stark_zeeman: + +Stark-Zeeman line shape +======================= + +In B. Lomanowski, et al. "Inferring divertor plasma properties from hydrogen Balmer +and Paschen series spectroscopy in JET-ILW." Nuclear Fusion 55.12 (2015) +`123028 `_ Stark and Zeeman features +are treated independently. This is a simplification and for better accuracy these effects +must be taken into account jointly as in J. Rosato, Y. Marandet, R. Stamm, +Journal of Quantitative Spectroscopy & Radiative Transfer 187 (2017) +`333 `_. + +The StarkBroadenedLine() follows the Lomanowski's paper, but introduces a couple of additional +approximations: + +* Zeeman splitting is taken in the form of a simple triplet with a :math:`\pi`-component + centred at :math:`\lambda`, :math:`\sigma^{+}`-component at :math:`\frac{hc}{hc/\lambda -\mu B}` + and :math:`\sigma^{-}`-component at :math:`\frac{hc}{hc/\lambda +\mu B}`. +* The convolution of Stark-Zeeman and Doppler profiles is replaced with the weighted sum + to speed-up calculations (`pseudo-Voigt `_ approximation). + +This example calculates Balmer-alpha and Paschen-beta Stark-Zeeman spectral lines +for the same parameters of plasma as in Figure 2 in Lomanowski's paper. + +.. literalinclude:: ../../../../demos/emission_models/stark_zeeman.py + +.. figure:: stark_zeeman_balmer_alpha.png + :align: center + :width: 450px + + **Caption:** The Stark-Zeeman structure of the Balmer-alpha line. + +.. figure:: stark_zeeman_paschen_beta.png + :align: center + :width: 450px + + **Caption:** The Stark-Zeeman structure of the Paschen-beta line. diff --git a/docs/source/demonstrations/passive_spectroscopy/stark_zeeman_balmer_alpha.png b/docs/source/demonstrations/passive_spectroscopy/stark_zeeman_balmer_alpha.png new file mode 100644 index 00000000..532a8eab Binary files /dev/null and b/docs/source/demonstrations/passive_spectroscopy/stark_zeeman_balmer_alpha.png differ diff --git a/docs/source/demonstrations/passive_spectroscopy/stark_zeeman_paschen_beta.png b/docs/source/demonstrations/passive_spectroscopy/stark_zeeman_paschen_beta.png new file mode 100644 index 00000000..afa5044a Binary files /dev/null and b/docs/source/demonstrations/passive_spectroscopy/stark_zeeman_paschen_beta.png differ diff --git a/docs/source/models/line_shapes/spectral_line_shapes.rst b/docs/source/models/line_shapes/spectral_line_shapes.rst index 5dc6325e..0d7fac55 100644 --- a/docs/source/models/line_shapes/spectral_line_shapes.rst +++ b/docs/source/models/line_shapes/spectral_line_shapes.rst @@ -2,36 +2,46 @@ Spectral Line Shapes ==================== -Cherab contains Doppler-broadened, Stark-broadened and Doppler-Zeeman line shapes of -atomic spectra. Stark-Doppler and Stark-Doppler-Zeeman line shapes of hydrogen spectra -will be added in the future. +Cherab contains Doppler-broadened, Doppler-Zeeman and Stark-Doppler-Zeeman line shapes of +atomic spectra. **Assumption: Maxwellian distribution of emitting species is assumed.** **A general model of Doppler broadening will be implemented in the future.** -.. autoclass:: cherab.core.model.lineshape.add_gaussian_line +.. autoclass:: cherab.core.model.lineshape.doppler.doppler_shift -.. autoclass:: cherab.core.model.lineshape.LineShapeModel +.. autoclass:: cherab.core.model.lineshape.doppler.thermal_broadening + +.. autoclass:: cherab.core.model.lineshape.gaussian.add_gaussian_line + +.. autoclass:: cherab.core.model.lineshape.stark.add_lorentzian_line + +.. autoclass:: cherab.core.model.lineshape.base.LineShapeModel :members: -.. autoclass:: cherab.core.model.lineshape.GaussianLine +.. autoclass:: cherab.core.model.lineshape.gaussian.GaussianLine :members: -.. autoclass:: cherab.core.model.lineshape.MultipletLineShape +.. autoclass:: cherab.core.model.lineshape.multiplet.MultipletLineShape :members: -.. autoclass:: cherab.core.model.lineshape.StarkBroadenedLine +.. autoclass:: cherab.core.model.lineshape.zeeman.ZeemanLineShapeModel :members: -.. autoclass:: cherab.core.model.lineshape.ZeemanLineShapeModel +.. autoclass:: cherab.core.model.lineshape.zeeman.ZeemanTriplet :members: -.. autoclass:: cherab.core.model.lineshape.ZeemanTriplet +.. autoclass:: cherab.core.model.lineshape.zeeman.ParametrisedZeemanTriplet :members: -.. autoclass:: cherab.core.model.lineshape.ParametrisedZeemanTriplet +.. autoclass:: cherab.core.model.lineshape.zeeman.ZeemanMultiplet :members: -.. autoclass:: cherab.core.model.lineshape.ZeemanMultiplet +.. autoclass:: cherab.core.model.lineshape.stark.StarkBroadenedLine :members: +.. autoclass:: cherab.core.model.lineshape.beam.base.BeamLineShapeModel + :members: + +.. autoclass:: cherab.core.model.lineshape.beam.mse.BeamEmissionMultiplet + :members: