diff --git a/cherab/core/math/interpolators/interpolators1d.pyx b/cherab/core/math/interpolators/interpolators1d.pyx index 1890bf12..773f28db 100644 --- a/cherab/core/math/interpolators/interpolators1d.pyx +++ b/cherab/core/math/interpolators/interpolators1d.pyx @@ -270,6 +270,9 @@ cdef class _Interpolate1DBase(Function1D): cdef class Interpolate1DLinear(_Interpolate1DBase): """ + .. deprecated:: 1.3.0 + Use `raysect.core.math.function.float.Interpolator1DArray` instead. + Interpolates 1D data using linear interpolation. Inherits from Function1D, implements `__call__(x)`. @@ -363,6 +366,9 @@ cdef class Interpolate1DLinear(_Interpolate1DBase): cdef class Interpolate1DCubic(_Interpolate1DBase): """ + .. deprecated:: 1.3.0 + Use `raysect.core.math.function.float.Interpolator1DArray` instead. + Interpolates 1D data using cubic interpolation. Inherits from Function1D, implements `__call__(x)`. diff --git a/cherab/core/math/interpolators/interpolators2d.pyx b/cherab/core/math/interpolators/interpolators2d.pyx index 08cab0ba..911cf25b 100644 --- a/cherab/core/math/interpolators/interpolators2d.pyx +++ b/cherab/core/math/interpolators/interpolators2d.pyx @@ -333,6 +333,9 @@ cdef class _Interpolate2DBase(Function2D): cdef class Interpolate2DLinear(_Interpolate2DBase): """ + .. deprecated:: 1.3.0 + Use `raysect.core.math.function.float.Interpolator2DArray` instead. + Interpolates 2D data using linear interpolation. Inherits from Function2D, implements `__call__(x, y)`. @@ -496,6 +499,9 @@ cdef class Interpolate2DLinear(_Interpolate2DBase): cdef class Interpolate2DCubic(_Interpolate2DBase): """ + .. deprecated:: 1.3.0 + Use `raysect.core.math.function.float.Interpolator2DArray` instead. + Interpolates 2D data using cubic interpolation. Inherits from Function2D, implements `__call__(x, y)`. diff --git a/cherab/core/math/interpolators/interpolators3d.pyx b/cherab/core/math/interpolators/interpolators3d.pyx index b38011d3..183c2335 100644 --- a/cherab/core/math/interpolators/interpolators3d.pyx +++ b/cherab/core/math/interpolators/interpolators3d.pyx @@ -330,6 +330,9 @@ cdef class _Interpolate3DBase(Function3D): cdef class Interpolate3DLinear(_Interpolate3DBase): """ + .. deprecated:: 1.3.0 + Use `raysect.core.math.function.float.Interpolator3DArray` instead. + Interpolates 3D data using linear interpolation. Inherits from Function3D, implements `__call__(x, y, z)`. @@ -489,6 +492,9 @@ cdef class Interpolate3DLinear(_Interpolate3DBase): cdef class Interpolate3DCubic(_Interpolate3DBase): """ + .. deprecated:: 1.3.0 + Use `raysect.core.math.function.float.Interpolator3DArray` instead. + Interpolates 3D data using cubic interpolation. Inherits from Function3D, implements `__call__(x, y, z)`. diff --git a/cherab/core/math/mappers.pyx b/cherab/core/math/mappers.pyx index 8022bbc2..5ff224a3 100644 --- a/cherab/core/math/mappers.pyx +++ b/cherab/core/math/mappers.pyx @@ -37,7 +37,8 @@ cdef class IsoMapper2D(Function2D): .. code-block:: pycon - >>> from cherab.core.math import IsoMapper2D, Interpolate1DCubic + >>> from raysect.core.math.function.float import Interpolator1DArray + >>> from cherab.core.math import IsoMapper2D >>> from cherab.tools.equilibrium import example_equilibrium >>> >>> equilibrium = example_equilibrium() @@ -45,7 +46,7 @@ cdef class IsoMapper2D(Function2D): >>> # extract the 2D psi function >>> psi_n = equilibrium.psi_normalised >>> # make a 1D psi profile - >>> profile = Interpolate1DCubic([0, 0.5, 0.9, 1.0], [2500, 2000, 1000, 0]) + >>> profile = Interpolator1DArray([0, 0.5, 0.9, 1.0], [2500, 2000, 1000, 0], 'cubic', 'none', 0) >>> # perform the flux function mapping >>> f = IsoMapper2D(psi_n, profile) >>> @@ -81,7 +82,8 @@ cdef class IsoMapper3D(Function3D): .. code-block:: pycon - >>> from cherab.core.math import IsoMapper2D, Interpolate1DCubic, AxisymmetricMapper + >>> from raysect.core.math.function.float import Interpolator1DArray + >>> from cherab.core.math import IsoMapper2D, AxisymmetricMapper >>> from cherab.tools.equilibrium import example_equilibrium >>> >>> equilibrium = example_equilibrium() @@ -90,7 +92,7 @@ cdef class IsoMapper3D(Function3D): >>> psi_n = equilibrium.psi_normalised >>> psi_n_3d = AxisymmetricMapper(psi_n) >>> # make a 1D psi profile - >>> profile = Interpolate1DCubic([0, 0.5, 0.9, 1.0], [2500, 2000, 1000, 0]) + >>> profile = Interpolator1DArray([0, 0.5, 0.9, 1.0], [2500, 2000, 1000, 0], 'cubic', 'none', 0) >>> # perform the flux function mapping >>> f = IsoMapper3D(psi_n_3d, profile) >>> diff --git a/cherab/core/model/attenuator/singleray.pxd b/cherab/core/model/attenuator/singleray.pxd index 6681f5fc..66ace689 100644 --- a/cherab/core/model/attenuator/singleray.pxd +++ b/cherab/core/model/attenuator/singleray.pxd @@ -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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); diff --git a/cherab/core/model/attenuator/singleray.pyx b/cherab/core/model/attenuator/singleray.pyx index e9e5cd93..a2cfbde1 100644 --- a/cherab/core/model/attenuator/singleray.pyx +++ b/cherab/core/model/attenuator/singleray.pyx @@ -23,9 +23,9 @@ import numpy as np cimport numpy as np from raysect.optical cimport AffineMatrix3D, Point3D, Vector3D, new_point3d +from raysect.core.math.function.float cimport Interpolator1DArray from cherab.core.utility import EvAmuToMS, EvToJ from cherab.core.atomic cimport BeamStoppingRate, AtomicData -from cherab.core.math cimport Interpolate1DLinear from cherab.core.plasma cimport Plasma from cherab.core.beam cimport Beam from cherab.core.species cimport Species @@ -178,7 +178,7 @@ cdef class SingleRayAttenuator(BeamAttenuator): self._tanydiv = tan(DEGREES_TO_RADIANS * self._beam.divergence_y) # a tiny degree of extrapolation is permitted to handle numerical accuracy issues with the end of the array - self._density = Interpolate1DLinear(beam_z, beam_density, extrapolate=True, extrapolation_range=1e-9) + self._density = Interpolator1DArray(beam_z, beam_density, 'linear', 'nearest', extrapolation_range=1e-9) @cython.cdivision(True) @cython.boundscheck(False) diff --git a/cherab/openadas/rates/atomic.pxd b/cherab/openadas/rates/atomic.pxd index cd595743..738aa5d4 100644 --- a/cherab/openadas/rates/atomic.pxd +++ b/cherab/openadas/rates/atomic.pxd @@ -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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -19,19 +19,18 @@ from cherab.core cimport IonisationRate as CoreIonisationRate from cherab.core cimport RecombinationRate as CoreRecombinationRate from cherab.core cimport ThermalCXRate as CoreThermalCXRate -from cherab.core.math cimport Interpolate2DCubic +from cherab.core.math cimport Function2D cdef class IonisationRate(CoreIonisationRate): cdef: readonly dict raw_data - readonly double wavelength readonly tuple density_range, temperature_range - Interpolate2DCubic _rate + Function2D _rate -cdef class NullImpactExcitationRate(CoreIonisationRate): +cdef class NullIonisationRate(CoreIonisationRate): pass @@ -39,9 +38,8 @@ cdef class RecombinationRate(CoreRecombinationRate): cdef: readonly dict raw_data - readonly double wavelength readonly tuple density_range, temperature_range - Interpolate2DCubic _rate + Function2D _rate cdef class NullRecombinationRate(CoreRecombinationRate): @@ -52,9 +50,8 @@ cdef class ThermalCXRate(CoreThermalCXRate): cdef: readonly dict raw_data - readonly double wavelength readonly tuple density_range, temperature_range - Interpolate2DCubic _rate + Function2D _rate cdef class NullThermalCXRate(CoreThermalCXRate): diff --git a/cherab/openadas/rates/atomic.pyx b/cherab/openadas/rates/atomic.pyx index 6a7a766e..04500124 100644 --- a/cherab/openadas/rates/atomic.pyx +++ b/cherab/openadas/rates/atomic.pyx @@ -1,7 +1,7 @@ -# 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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -18,6 +18,9 @@ # under the Licence. import numpy as np +from libc.math cimport INFINITY, log10 + +from raysect.core.math.function.float cimport Interpolator2DArray cdef class IonisationRate(CoreIonisationRate): @@ -33,21 +36,28 @@ cdef class IonisationRate(CoreIonisationRate): # unpack ne = data['ne'] te = data['te'] - rate = np.log10(data['rate']) + rate = np.log10(data['rate']) # store limits of data self.density_range = ne.min(), ne.max() self.temperature_range = te.min(), te.max() # interpolate rate - self._rate = Interpolate2DCubic( - ne, te, rate, extrapolate=extrapolate, extrapolation_type="quadratic" - ) + # using nearest extrapolation to avoid infinite values at 0 for some rates + extrapolation_type = 'nearest' if extrapolate else 'none' + self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY) cpdef double evaluate(self, double density, double temperature) except? -1e999: + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if density < 1.e-300: + density = 1.e-300 + + if temperature < 1.e-300: + temperature = 1.e-300 + # calculate rate and convert from log10 space to linear space - return 10 ** self._rate.evaluate(density, temperature) + return 10 ** self._rate.evaluate(log10(density), log10(temperature)) cdef class NullIonisationRate(CoreIonisationRate): @@ -73,22 +83,28 @@ cdef class RecombinationRate(CoreRecombinationRate): # unpack ne = data['ne'] te = data['te'] - rate = np.log10(data['rate']) + rate = np.log10(data['rate']) # store limits of data self.density_range = ne.min(), ne.max() self.temperature_range = te.min(), te.max() # interpolate rate - self._rate = Interpolate2DCubic( - ne, te, rate, extrapolate=extrapolate, extrapolation_type="quadratic" - ) - + # using nearest extrapolation to avoid infinite values at 0 for some rates + extrapolation_type = 'nearest' if extrapolate else 'none' + self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY) cpdef double evaluate(self, double density, double temperature) except? -1e999: + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if density < 1.e-300: + density = 1.e-300 + + if temperature < 1.e-300: + temperature = 1.e-300 + # calculate rate and convert from log10 space to linear space - return 10 ** self._rate.evaluate(density, temperature) + return 10 ** self._rate.evaluate(log10(density), log10(temperature)) cdef class NullRecombinationRate(CoreRecombinationRate): @@ -114,21 +130,27 @@ cdef class ThermalCXRate(CoreThermalCXRate): # unpack ne = data['ne'] te = data['te'] - rate = np.log10(data['rate']) + rate = np.log10(data['rate']) # store limits of data self.density_range = ne.min(), ne.max() self.temperature_range = te.min(), te.max() # interpolate rate - self._rate = Interpolate2DCubic( - ne, te, rate, extrapolate=extrapolate, extrapolation_type="quadratic" - ) + extrapolation_type = 'linear' if extrapolate else 'none' + self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY) cpdef double evaluate(self, double density, double temperature) except? -1e999: + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if density < 1.e-300: + density = 1.e-300 + + if temperature < 1.e-300: + temperature = 1.e-300 + # calculate rate and convert from log10 space to linear space - return 10 ** self._rate.evaluate(density, temperature) + return 10 ** self._rate.evaluate(log10(density), log10(temperature)) cdef class NullThermalCXRate(CoreThermalCXRate): @@ -138,4 +160,4 @@ cdef class NullThermalCXRate(CoreThermalCXRate): """ cpdef double evaluate(self, double density, double temperature) except? -1e999: - return 0.0 \ No newline at end of file + return 0.0 diff --git a/cherab/openadas/rates/beam.pxd b/cherab/openadas/rates/beam.pxd index f990d095..4f23f2bf 100644 --- a/cherab/openadas/rates/beam.pxd +++ b/cherab/openadas/rates/beam.pxd @@ -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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); diff --git a/cherab/openadas/rates/beam.pyx b/cherab/openadas/rates/beam.pyx index 5b4da2db..f40af8dd 100644 --- a/cherab/openadas/rates/beam.pyx +++ b/cherab/openadas/rates/beam.pyx @@ -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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -16,10 +16,14 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. +import numpy as np + from cherab.core.utility.conversion import PhotonToJ cimport cython -from cherab.core.math cimport Interpolate1DCubic, Interpolate2DCubic +from libc.math cimport INFINITY, log10 +from raysect.core.math.function.float cimport Interpolator1DArray, Interpolator2DArray, Constant2D, Arg2D +from cherab.core.math cimport IsoMapper2D # todo: clarify variables @@ -40,8 +44,8 @@ cdef class BeamStoppingRate(CoreBeamStoppingRate): e = data["e"] # eV/amu n = data["n"] # m^-3 t = data["t"] # eV - sen = data["sen"] # m^3/s - st = data["st"] / data["sref"] # dimensionless + sen = np.log10(data["sen"]) # m^3/s + st = np.log10(data["st"] / data["sref"]) # dimensionless # store limits of data self.beam_energy_range = e.min(), e.max() @@ -49,8 +53,17 @@ cdef class BeamStoppingRate(CoreBeamStoppingRate): self.temperature_range = t.min(), t.max() # interpolate - self._npl_eb = Interpolate2DCubic(e, n, sen, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") - self._tp = Interpolate1DCubic(t, st, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") + extrapolation_type_2d = 'linear' if extrapolate else 'none' + extrapolation_type_1d = 'quadratic' if extrapolate else 'none' + if len(e) == 1 and len(n) == 1: + self._npl_eb = Constant2D(sen[0, 0]) + elif len(e) == 1: + self._npl_eb = IsoMapper2D(Arg2D('y'), Interpolator1DArray(np.log10(n), sen[0], 'cubic', extrapolation_type_1d, INFINITY)) + elif len(n) == 1: + self._npl_eb = IsoMapper2D(Arg2D('x'), Interpolator1DArray(np.log10(e), sen[:, 0], 'cubic', extrapolation_type_1d, INFINITY)) + else: + self._npl_eb = Interpolator2DArray(np.log10(e), np.log10(n), sen, 'cubic', extrapolation_type_2d, INFINITY, INFINITY) + self._tp = Interpolator1DArray(np.log10(t), st, 'cubic', extrapolation_type_1d, INFINITY) cpdef double evaluate(self, double energy, double density, double temperature) except? -1e999: """ @@ -64,17 +77,18 @@ cdef class BeamStoppingRate(CoreBeamStoppingRate): :return: The beam stopping coefficient in m^3.s^-1 """ - cdef double val + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if energy < 1.e-300: + energy = 1.e-300 - val = self._npl_eb.evaluate(energy, density) - if val <= 0: - return 0.0 + if density < 1.e-300: + density = 1.e-300 - val *= self._tp.evaluate(temperature) - if val <= 0: - return 0.0 + if temperature < 1.e-300: + temperature = 1.e-300 - return val + # calculate rate and convert from log10 space to linear space + return 10 ** (self._npl_eb.evaluate(log10(energy), log10(density)) + self._tp.evaluate(log10(temperature))) cdef class NullBeamStoppingRate(CoreBeamStoppingRate): @@ -104,8 +118,8 @@ cdef class BeamPopulationRate(CoreBeamPopulationRate): e = data["e"] # eV/amu n = data["n"] # m^-3 t = data["t"] # eV - sen = data["sen"] # dimensionless - st = data["st"] / data["sref"] # dimensionless + sen = np.log10(data["sen"]) # dimensionless + st = np.log10(data["st"] / data["sref"]) # dimensionless # store limits of data self.beam_energy_range = e.min(), e.max() @@ -113,8 +127,17 @@ cdef class BeamPopulationRate(CoreBeamPopulationRate): self.temperature_range = t.min(), t.max() # interpolate - self._npl_eb = Interpolate2DCubic(e, n, sen, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") - self._tp = Interpolate1DCubic(t, st, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") + extrapolation_type_2d = 'linear' if extrapolate else 'none' + extrapolation_type_1d = 'quadratic' if extrapolate else 'none' + if len(e) == 1 and len(n) == 1: + self._npl_eb = Constant2D(sen[0, 0]) + elif len(e) == 1: + self._npl_eb = IsoMapper2D(Arg2D('y'), Interpolator1DArray(np.log10(n), sen[0], 'cubic', extrapolation_type_1d, INFINITY)) + elif len(n) == 1: + self._npl_eb = IsoMapper2D(Arg2D('x'), Interpolator1DArray(np.log10(e), sen[:, 0], 'cubic', extrapolation_type_1d, INFINITY)) + else: + self._npl_eb = Interpolator2DArray(np.log10(e), np.log10(n), sen, 'cubic', extrapolation_type_2d, INFINITY, INFINITY) + self._tp = Interpolator1DArray(np.log10(t), st, 'cubic', extrapolation_type_1d, INFINITY) cpdef double evaluate(self, double energy, double density, double temperature) except? -1e999: """ @@ -128,17 +151,18 @@ cdef class BeamPopulationRate(CoreBeamPopulationRate): :return: The beam population coefficient in dimensionless units. """ - cdef double val + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if energy < 1.e-300: + energy = 1.e-300 - val = self._npl_eb.evaluate(energy, density) - if val <= 0: - return 0.0 + if density < 1.e-300: + density = 1.e-300 - val *= self._tp.evaluate(temperature) - if val <= 0: - return 0.0 + if temperature < 1.e-300: + temperature = 1.e-300 - return val + # calculate rate and convert from log10 space to linear space + return 10 ** (self._npl_eb.evaluate(log10(energy), log10(density)) + self._tp.evaluate(log10(temperature))) cdef class NullBeamPopulationRate(CoreBeamPopulationRate): @@ -170,8 +194,8 @@ cdef class BeamEmissionPEC(CoreBeamEmissionPEC): e = data["e"] # eV/amu n = data["n"] # m^-3 t = data["t"] # eV - sen = PhotonToJ.to(data["sen"], wavelength) # W.m^3/s - st = data["st"] / data["sref"] # dimensionless + sen = np.log10(PhotonToJ.to(data["sen"], wavelength)) # W.m^3/s + st = np.log10(data["st"] / data["sref"]) # dimensionless # store limits of data self.beam_energy_range = e.min(), e.max() @@ -179,8 +203,17 @@ cdef class BeamEmissionPEC(CoreBeamEmissionPEC): self.temperature_range = t.min(), t.max() # interpolate - self._npl_eb = Interpolate2DCubic(e, n, sen, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") - self._tp = Interpolate1DCubic(t, st, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") + extrapolation_type_2d = 'linear' if extrapolate else 'none' + extrapolation_type_1d = 'quadratic' if extrapolate else 'none' + if len(e) == 1 and len(n) == 1: + self._npl_eb = Constant2D(sen[0, 0]) + elif len(e) == 1: + self._npl_eb = IsoMapper2D(Arg2D('y'), Interpolator1DArray(np.log10(n), sen[0], 'cubic', extrapolation_type_1d, INFINITY)) + elif len(n) == 1: + self._npl_eb = IsoMapper2D(Arg2D('x'), Interpolator1DArray(np.log10(e), sen[:, 0], 'cubic', extrapolation_type_1d, INFINITY)) + else: + self._npl_eb = Interpolator2DArray(np.log10(e), np.log10(n), sen, 'cubic', extrapolation_type_2d, INFINITY, INFINITY) + self._tp = Interpolator1DArray(np.log10(t), st, 'cubic', extrapolation_type_1d, INFINITY) cpdef double evaluate(self, double energy, double density, double temperature) except? -1e999: """ @@ -194,17 +227,18 @@ cdef class BeamEmissionPEC(CoreBeamEmissionPEC): :return: The beam emission coefficient in m^3.s^-1 """ - cdef double val + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if energy < 1.e-300: + energy = 1.e-300 - val = self._npl_eb.evaluate(energy, density) - if val <= 0: - return 0.0 + if density < 1.e-300: + density = 1.e-300 - val *= self._tp.evaluate(temperature) - if val <= 0: - return 0.0 + if temperature < 1.e-300: + temperature = 1.e-300 - return val + # calculate rate and convert from log10 space to linear space + return 10 ** (self._npl_eb.evaluate(log10(energy), log10(density)) + self._tp.evaluate(log10(temperature))) cdef class NullBeamEmissionPEC(CoreBeamEmissionPEC): diff --git a/cherab/openadas/rates/cx.pxd b/cherab/openadas/rates/cx.pxd index 47fdf75d..393027d0 100644 --- a/cherab/openadas/rates/cx.pxd +++ b/cherab/openadas/rates/cx.pxd @@ -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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -17,7 +17,7 @@ # under the Licence. from cherab.core cimport BeamCXPEC as CoreBeamCXPEC -from cherab.core.math cimport Function1D, Function2D +from cherab.core.math cimport Function1D cdef class BeamCXPEC(CoreBeamCXPEC): diff --git a/cherab/openadas/rates/cx.pyx b/cherab/openadas/rates/cx.pyx index ece84cf4..cb827a8b 100644 --- a/cherab/openadas/rates/cx.pyx +++ b/cherab/openadas/rates/cx.pyx @@ -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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -16,10 +16,12 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. +import numpy as np from cherab.core.utility.conversion import Cm3ToM3, PerCm3ToPerM3, PhotonToJ cimport cython -from cherab.core.math cimport Interpolate1DCubic +from libc.math cimport INFINITY, log10 +from raysect.core.math.function.float cimport Interpolator1DArray, Constant1D cdef class BeamCXPEC(CoreBeamCXPEC): @@ -47,7 +49,7 @@ cdef class BeamCXPEC(CoreBeamCXPEC): bmag = data["b"] # Tesla qref = data["qref"] # m^3/s - qeb = PhotonToJ.to(data["qeb"], wavelength) # W.m^3 + qeb = np.log10(PhotonToJ.to(data["qeb"], wavelength)) # W.m^3 qti = data["qti"] / qref # dimensionless qni = data["qni"] / qref # dimensionless qzeff = data["qz"] / qref # dimensionless @@ -61,11 +63,13 @@ cdef class BeamCXPEC(CoreBeamCXPEC): self.b_field_range = bmag.min(), bmag.max() # interpolate the rate data - self._eb = Interpolate1DCubic(eb, qeb, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") - self._ti = Interpolate1DCubic(ti, qti, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") - self._ni = Interpolate1DCubic(ni, qni, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") - self._zeff = Interpolate1DCubic(zeff, qzeff, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") - self._b = Interpolate1DCubic(bmag, qbmag, tolerate_single_value=True, extrapolate=extrapolate, extrapolation_type="quadratic") + extrapolation_type_log = 'quadratic' if extrapolate else 'none' + extrapolation_type = 'nearest' if extrapolate else 'none' + self._eb = Interpolator1DArray(np.log10(eb), qeb, 'cubic', extrapolation_type_log, INFINITY) if len(qeb) > 1 else Constant1D(qeb[0]) + self._ti = Interpolator1DArray(ti, qti, 'cubic', extrapolation_type, INFINITY) if len(qti) > 1 else Constant1D(qti[0]) + self._ni = Interpolator1DArray(ni, qni, 'cubic', extrapolation_type, INFINITY) if len(qni) > 1 else Constant1D(qni[0]) + self._zeff = Interpolator1DArray(zeff, qzeff, 'cubic', extrapolation_type, INFINITY) if len(qzeff) > 1 else Constant1D(qzeff[0]) + self._b = Interpolator1DArray(bmag, qbmag, 'cubic', extrapolation_type, INFINITY) if len(qbmag) > 1 else Constant1D(qbmag[0]) cpdef double evaluate(self, double energy, double temperature, double density, double z_effective, double b_field) except? -1e999: """ @@ -83,9 +87,11 @@ cdef class BeamCXPEC(CoreBeamCXPEC): cdef double rate - rate = self._eb.evaluate(energy) - if rate <= 0: - return 0.0 + # need to handle zeros for log-log interpolation + if energy < 1.e-300: + energy = 1.e-300 + + rate = 10 ** self._eb.evaluate(log10(energy)) rate *= self._ti.evaluate(temperature) if rate <= 0: diff --git a/cherab/openadas/rates/pec.pxd b/cherab/openadas/rates/pec.pxd index 435b1721..46a54cbe 100644 --- a/cherab/openadas/rates/pec.pxd +++ b/cherab/openadas/rates/pec.pxd @@ -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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -19,7 +19,7 @@ from cherab.core cimport ImpactExcitationPEC as CoreImpactExcitationPEC from cherab.core cimport RecombinationPEC as CoreRecombinationPEC from cherab.core cimport ThermalCXPEC as CoreThermalCXPEC -from cherab.core.math cimport Interpolate2DCubic +from cherab.core.math cimport Function2D cdef class ImpactExcitationPEC(CoreImpactExcitationPEC): @@ -28,7 +28,7 @@ cdef class ImpactExcitationPEC(CoreImpactExcitationPEC): readonly dict raw_data readonly double wavelength readonly tuple density_range, temperature_range - Interpolate2DCubic _rate + Function2D _rate cdef class NullImpactExcitationPEC(CoreImpactExcitationPEC): @@ -41,7 +41,7 @@ cdef class RecombinationPEC(CoreRecombinationPEC): readonly dict raw_data readonly double wavelength readonly tuple density_range, temperature_range - Interpolate2DCubic _rate + Function2D _rate cdef class NullRecombinationPEC(CoreRecombinationPEC): diff --git a/cherab/openadas/rates/pec.pyx b/cherab/openadas/rates/pec.pyx index c2f3bc4e..eafc6c64 100644 --- a/cherab/openadas/rates/pec.pyx +++ b/cherab/openadas/rates/pec.pyx @@ -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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -17,9 +17,10 @@ # under the Licence. import numpy as np -cimport numpy as np -import matplotlib.pyplot as plt +from libc.math cimport INFINITY, log10 + +from raysect.core.math.function.float cimport Interpolator2DArray from cherab.core.utility.conversion import PhotonToJ @@ -38,24 +39,31 @@ cdef class ImpactExcitationPEC(CoreImpactExcitationPEC): # unpack ne = data['ne'] te = data['te'] - rate = data['rate'] + rate = data['rate'] # pre-convert data to W m^3 from Photons s^-1 cm^3 prior to interpolation - rate = PhotonToJ.to(rate, wavelength) + rate = np.log10(PhotonToJ.to(rate, wavelength)) # store limits of data self.density_range = ne.min(), ne.max() self.temperature_range = te.min(), te.max() # interpolate rate - self._rate = Interpolate2DCubic( - ne, te, rate, extrapolate=extrapolate, extrapolation_type="quadratic" - ) + # using nearest extrapolation to avoid infinite values at 0 for some rates + extrapolation_type = 'nearest' if extrapolate else 'none' + self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY) cpdef double evaluate(self, double density, double temperature) except? -1e999: - # prevent -ve values (possible if extrapolation enabled) - return max(0, self._rate.evaluate(density, temperature)) + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if density < 1.e-300: + density = 1.e-300 + + if temperature < 1.e-300: + temperature = 1.e-300 + + # calculate rate and convert from log10 space to linear space + return 10 ** self._rate.evaluate(log10(density), log10(temperature)) cdef class NullImpactExcitationPEC(CoreImpactExcitationPEC): @@ -83,25 +91,31 @@ cdef class RecombinationPEC(CoreRecombinationPEC): # unpack ne = data['ne'] te = data['te'] - rate = data['rate'] + rate = data['rate'] # pre-convert data to W m^3 from Photons s^-1 cm^3 prior to interpolation - rate = PhotonToJ.to(rate, wavelength) + rate = np.log10(PhotonToJ.to(rate, wavelength)) # store limits of data self.density_range = ne.min(), ne.max() self.temperature_range = te.min(), te.max() # interpolate rate - self._rate = Interpolate2DCubic( - ne, te, rate, extrapolate=extrapolate, extrapolation_type="quadratic" - ) - + # using nearest extrapolation to avoid infinite values at 0 for some rates + extrapolation_type = 'nearest' if extrapolate else 'none' + self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY) cpdef double evaluate(self, double density, double temperature) except? -1e999: - # prevent -ve values (possible if extrapolation enabled) - return max(0, self._rate.evaluate(density, temperature)) + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if density < 1.e-300: + density = 1.e-300 + + if temperature < 1.e-300: + temperature = 1.e-300 + + # calculate rate and convert from log10 space to linear space + return 10 ** self._rate.evaluate(log10(density), log10(temperature)) cdef class NullRecombinationPEC(CoreRecombinationPEC): @@ -112,39 +126,3 @@ cdef class NullRecombinationPEC(CoreRecombinationPEC): cpdef double evaluate(self, double density, double temperature) except? -1e999: return 0.0 - - - -# cdef class ThermalCXRate(CoreThermalCXRate): -# -# def __init__(self, double wavelength ,dict rate_data, extrapolate=False): -# pass -# cpdef double evaluate(self, double density, double temperature): -# pass -# -# def plot(self, density_list=None, x_limit=None, y_limit=None): -# -# plt.figure() -# for dens in density_list: -# rates = [self.__call__(dens, temp) for temp in self._temperature] -# plt.loglog(self._temperature, rates, label='{:.4G} m$^{{-3}}$'.format(dens)) -# -# if x_limit: -# plt.xlim(x_limit) -# if y_limit: -# plt.ylim(y_limit) -# -# plt.legend(loc=4) -# plt.xlabel('Electron Temperature (eV)') -# plt.ylabel('$PEC$ (m$^3$ s$^{-1}$)') -# plt.title('Photon emissivity coefficient') -# -# -# cdef class NullThermalCXRate(CoreThermalCXRate): -# """ -# A PEC rate that always returns zero. -# Needed for use cases where the required atomic data is missing. -# """ -# -# cpdef double evaluate(self, double density, double temperature) except? -1e999: -# return 0.0 diff --git a/cherab/openadas/rates/radiated_power.pxd b/cherab/openadas/rates/radiated_power.pxd index 438aafb5..73fa09aa 100644 --- a/cherab/openadas/rates/radiated_power.pxd +++ b/cherab/openadas/rates/radiated_power.pxd @@ -1,7 +1,7 @@ -# 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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -18,7 +18,7 @@ # under the Licence. -from cherab.core.math.interpolators.interpolators2d cimport Interpolate2DCubic +from cherab.core.math cimport Function2D from cherab.core.atomic.rates cimport LineRadiationPower as CoreLineRadiationPower from cherab.core.atomic.rates cimport ContinuumPower as CoreContinuumPower from cherab.core.atomic.rates cimport CXRadiationPower as CoreCXRadiationPower @@ -29,7 +29,7 @@ cdef class LineRadiationPower(CoreLineRadiationPower): cdef: readonly dict raw_data readonly tuple density_range, temperature_range - Interpolate2DCubic _rate + Function2D _rate cdef class NullLineRadiationPower(CoreLineRadiationPower): @@ -41,7 +41,7 @@ cdef class ContinuumPower(CoreContinuumPower): cdef: readonly dict raw_data readonly tuple density_range, temperature_range - Interpolate2DCubic _rate + Function2D _rate cdef class NullContinuumPower(CoreContinuumPower): @@ -53,7 +53,7 @@ cdef class CXRadiationPower(CoreCXRadiationPower): cdef: readonly dict raw_data readonly tuple density_range, temperature_range - Interpolate2DCubic _rate + Function2D _rate cdef class NullCXRadiationPower(CoreCXRadiationPower): diff --git a/cherab/openadas/rates/radiated_power.pyx b/cherab/openadas/rates/radiated_power.pyx index 50111281..570ced92 100644 --- a/cherab/openadas/rates/radiated_power.pyx +++ b/cherab/openadas/rates/radiated_power.pyx @@ -1,7 +1,7 @@ -# 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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -17,6 +17,11 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. +import numpy as np +from libc.math cimport INFINITY, log10 + +from raysect.core.math.function.float cimport Interpolator2DArray + cdef class LineRadiationPower(CoreLineRadiationPower): """Base class for radiated powers.""" @@ -30,21 +35,28 @@ cdef class LineRadiationPower(CoreLineRadiationPower): # unpack ne = data['ne'] te = data['te'] - rate = data['rate'] + rate = np.log10(data['rate']) # store limits of data self.density_range = ne.min(), ne.max() self.temperature_range = te.min(), te.max() # interpolate rate - self._rate = Interpolate2DCubic( - ne, te, rate, extrapolate=extrapolate, extrapolation_type="quadratic" - ) + # using nearest extrapolation to avoid infinite values at 0 for some rates + extrapolation_type = 'nearest' if extrapolate else 'none' + self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY) cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: - # prevent -ve values (possible if extrapolation enabled) - return max(0, self._rate.evaluate(electron_density, electron_temperature)) + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if electron_density < 1.e-300: + electron_density = 1.e-300 + + if electron_temperature < 1.e-300: + electron_temperature = 1.e-300 + + # calculate rate and convert from log10 space to linear space + return 10 ** self._rate.evaluate(log10(electron_density), log10(electron_temperature)) cdef class NullLineRadiationPower(CoreLineRadiationPower): @@ -69,21 +81,28 @@ cdef class ContinuumPower(CoreContinuumPower): # unpack ne = data['ne'] te = data['te'] - rate = data['rate'] + rate = np.log10(data['rate']) # store limits of data self.density_range = ne.min(), ne.max() self.temperature_range = te.min(), te.max() # interpolate rate - self._rate = Interpolate2DCubic( - ne, te, rate, extrapolate=extrapolate, extrapolation_type="quadratic" - ) + # using nearest extrapolation to avoid infinite values at 0 for some rates + extrapolation_type = 'nearest' if extrapolate else 'none' + self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY) cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: - # prevent -ve values (possible if extrapolation enabled) - return max(0, self._rate.evaluate(electron_density, electron_temperature)) + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if electron_density < 1.e-300: + electron_density = 1.e-300 + + if electron_temperature < 1.e-300: + electron_temperature = 1.e-300 + + # calculate rate and convert from log10 space to linear space + return 10 ** self._rate.evaluate(log10(electron_density), log10(electron_temperature)) cdef class NullContinuumPower(CoreContinuumPower): @@ -108,21 +127,27 @@ cdef class CXRadiationPower(CoreCXRadiationPower): # unpack ne = data['ne'] te = data['te'] - rate = data['rate'] + rate = np.log10(data['rate']) # store limits of data self.density_range = ne.min(), ne.max() self.temperature_range = te.min(), te.max() # interpolate rate - self._rate = Interpolate2DCubic( - ne, te, rate, extrapolate=extrapolate, extrapolation_type="quadratic" - ) + extrapolation_type = 'linear' if extrapolate else 'none' + self._rate = Interpolator2DArray(np.log10(ne), np.log10(te), rate, 'cubic', extrapolation_type, INFINITY, INFINITY) cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: - # prevent -ve values (possible if extrapolation enabled) - return max(0, self._rate.evaluate(electron_density, electron_temperature)) + # need to handle zeros, also density and temperature can become negative due to cubic interpolation + if electron_density < 1.e-300: + electron_density = 1.e-300 + + if electron_temperature < 1.e-300: + electron_temperature = 1.e-300 + + # calculate rate and convert from log10 space to linear space + return 10 ** self._rate.evaluate(log10(electron_density), log10(electron_temperature)) cdef class NullCXRadiationPower(CoreCXRadiationPower): diff --git a/cherab/tools/equilibrium/__init__.py b/cherab/tools/equilibrium/__init__.py index f4045c75..22809771 100644 --- a/cherab/tools/equilibrium/__init__.py +++ b/cherab/tools/equilibrium/__init__.py @@ -21,4 +21,3 @@ from .plot import plot_equilibrium from .eqdsk import import_eqdsk from .example import example_equilibrium -# from .equ import import_equ_psi diff --git a/cherab/tools/equilibrium/efit.pyx b/cherab/tools/equilibrium/efit.pyx index 92044bf8..955f6d39 100644 --- a/cherab/tools/equilibrium/efit.pyx +++ b/cherab/tools/equilibrium/efit.pyx @@ -26,11 +26,11 @@ cimport cython from raysect.optical cimport Vector3D, Point2D, new_vector3d from raysect.core.math.function.float cimport Blend2D as ScalarBlend2D from raysect.core.math.function.vector3d cimport Blend2D as VectorBlend2D +from raysect.core.math.function.float cimport Interpolator1DArray, Interpolator2DArray from cherab.core.math cimport Function1D, autowrap_function1d from cherab.core.math cimport Function2D, autowrap_function2d from cherab.core.math cimport VectorFunction2D, autowrap_vectorfunction2d, ConstantVector2D -from cherab.core.math cimport Interpolate1DCubic, Interpolate2DCubic from cherab.core.math cimport PolygonMask2D from cherab.core.math cimport IsoMapper2D, AxisymmetricMapper, VectorAxisymmetricMapper from cherab.core.math cimport ClampOutput2D @@ -109,18 +109,18 @@ cdef class EFITEquilibrium: self.psi_data = psi # interpolate poloidal flux grid data - self.psi = Interpolate2DCubic(r, z, psi) + self.psi = Interpolator2DArray(r, z, psi, 'cubic', 'none', 0, 0) self.psi_axis = psi_axis self.psi_lcfs = psi_lcfs - self.psi_normalised = ClampOutput2D(Interpolate2DCubic(r, z, (psi - psi_axis) / (psi_lcfs - psi_axis)), min=0) + self.psi_normalised = ClampOutput2D(Interpolator2DArray(r, z, (psi - psi_axis) / (psi_lcfs - psi_axis), 'cubic', 'none', 0, 0), min=0) # store equilibrium attributes self.r_range = r.min(), r.max() self.z_range = z.min(), z.max() self._b_vacuum_magnitude = b_vacuum_magnitude self._b_vacuum_radius = b_vacuum_radius - self._f_profile = Interpolate1DCubic(f_profile[0, :], f_profile[1, :]) - self.q = Interpolate1DCubic(q_profile[0, :], q_profile[1, :]) + self._f_profile = Interpolator1DArray(f_profile[0, :], f_profile[1, :], 'cubic', 'none', 0) + self.q = Interpolator1DArray(q_profile[0, :], q_profile[1, :], 'cubic', 'none', 0) # populate points self._process_points(magnetic_axis, x_points, strike_points) @@ -187,8 +187,8 @@ cdef class EFITEquilibrium: dix_dr = 1.0 / np.gradient(r, edge_order=2) diy_dz = 1.0 / np.gradient(z, edge_order=2) - dpsi_dr = Interpolate2DCubic(r, z, dpsi_dix * dix_dr[:, np.newaxis]) - dpsi_dz = Interpolate2DCubic(r, z, dpsi_diy * diy_dz[np.newaxis, :]) + dpsi_dr = Interpolator2DArray(r, z, dpsi_dix * dix_dr[:, np.newaxis], 'cubic', 'none', 0, 0) + dpsi_dz = Interpolator2DArray(r, z, dpsi_diy * diy_dz[np.newaxis, :], 'cubic', 'none', 0, 0) return dpsi_dr, dpsi_dz @@ -213,7 +213,7 @@ cdef class EFITEquilibrium: return # interpolate sampled data, allowing a small bit of extrapolation to cope with numerical sampling accuracy - self.psin_to_r = Interpolate1DCubic(psin, r, extrapolate=True, extrapolation_range=SAMPLE_RESOLUTION, extrapolation_type='quadratic') + self.psin_to_r = Interpolator1DArray(psin, r, 'cubic', 'quadratic', SAMPLE_RESOLUTION) def map2d(self, object profile, double value_outside_lcfs=0.0): """ @@ -243,7 +243,7 @@ cdef class EFITEquilibrium: profile = autowrap_function1d(profile) else: profile = np.array(profile, np.float64) - profile = Interpolate1DCubic(profile[0, :], profile[1, :]) + profile = Interpolator1DArray(profile[0, :], profile[1, :], 'cubic', 'none', 0) # map around equilibrium f = IsoMapper2D(self.psi_normalised, profile) @@ -321,21 +321,21 @@ cdef class EFITEquilibrium: toroidal = autowrap_function1d(toroidal) else: toroidal = np.array(toroidal, np.float64) - toroidal = Interpolate1DCubic(toroidal[0, :], toroidal[1, :]) + toroidal = Interpolator1DArray(toroidal[0, :], toroidal[1, :], 'cubic', 'none', 0) # convert poloidal data to 1d function if not already a function object if isinstance(poloidal, Function1D) or callable(poloidal): poloidal = autowrap_function1d(poloidal) else: poloidal = np.array(poloidal, np.float64) - poloidal = Interpolate1DCubic(poloidal[0, :], poloidal[1, :]) + poloidal = Interpolator1DArray(poloidal[0, :], poloidal[1, :], 'cubic', 'none', 0) # convert normal data to 1d function if not already a function object if isinstance(normal, Function1D) or callable(normal): normal = autowrap_function1d(normal) else: normal = np.array(normal, np.float64) - normal = Interpolate1DCubic(normal[0, :], normal[1, :]) + normal = Interpolator1DArray(normal[0, :], normal[1, :], 'cubic', 'none', 0) v = FluxCoordToCartesian(self.b_field, self.psi_normalised, toroidal, poloidal, normal) diff --git a/cherab/tools/equilibrium/equ.py b/cherab/tools/equilibrium/equ.py deleted file mode 100644 index 2e0252a8..00000000 --- a/cherab/tools/equilibrium/equ.py +++ /dev/null @@ -1,121 +0,0 @@ -# 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. - -# todo: this code is now incomplete it should return an equilibrium object -# import re -# import numpy as np -# from cherab.core.math import Interpolate2DCubic -# -# -# def import_equ_psi(equ_file): -# """ -# Imports a 2D Psi grid from an equ mesh equilibrium. -# -# :param str equ_file: The .EQU mesh generator equilibrium file. -# :return: A 2D Psi(r,z) function -# :rtype: Interpolate2DCubic -# """ -# -# fh = open(equ_file, 'r') -# file_lines = fh.readlines() -# fh.close() -# -# # Load r array -# line_i = 0 -# while True: -# line = file_lines[line_i] -# match = re.match('^\s*r\(1:jm\);', line) -# -# if not match: -# line_i += 1 -# continue -# -# line_i += 1 -# r_values = [] -# line = file_lines[line_i] -# while line and not line.isspace(): -# print(line) -# values = line.split() -# print(values) -# for v in values: -# r_values.append(float(v)) -# -# line_i += 1 -# line = file_lines[line_i] -# -# jm = len(r_values) -# break -# r_values = np.array(r_values) -# -# # Load z array -# while True: -# line = file_lines[line_i] -# match = re.match('^\s*z\(1:km\);', line) -# -# if not match: -# line_i += 1 -# continue -# -# line_i += 1 -# z_values = [] -# line = file_lines[line_i] -# while not re.match('^\s*$', line): -# values = line.split() -# for v in values: -# z_values.append(float(v)) -# -# line_i += 1 -# line = file_lines[line_i] -# -# km = len(z_values) -# break -# z_values = np.array(z_values) -# -# # Load (r, z) array -# while True: -# line = file_lines[line_i] -# match = re.match('^\s*\(\(psi\(j,k\)-psib,j=1,jm\),k=1,km\)', line) -# -# if not match: -# line_i += 1 -# continue -# -# line_i += 1 -# psi_values = [] -# line = file_lines[line_i] -# while not re.match('^\s*$', line): -# values = line.split() -# for v in values: -# psi_values.append(float(v)) -# -# line_i += 1 -# try: -# line = file_lines[line_i] -# except IndexError: -# break -# -# if len(psi_values) != km * jm: -# raise ValueError("Number of values in r, z array not equal to (km, jm).") -# -# break -# -# psi_raw = np.array(psi_values) -# psi = psi_raw.reshape((km, jm)) -# psi = np.swapaxes(psi, 0, 1) -# -# return Interpolate2DCubic(r_values, z_values, psi) diff --git a/cherab/tools/plasmas/ionisation_balance.py b/cherab/tools/plasmas/ionisation_balance.py index bf14e37d..07a7ea90 100644 --- a/cherab/tools/plasmas/ionisation_balance.py +++ b/cherab/tools/plasmas/ionisation_balance.py @@ -1,7 +1,7 @@ -# 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-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -21,10 +21,11 @@ import numpy as np from scipy.optimize import lsq_linear from collections.abc import Iterable +from raysect.core.math.function.float import Function1D, Function2D, Interpolator1DArray, Interpolator2DArray from cherab.core import AtomicData from cherab.core.atomic import Element -from cherab.core.math import Interpolate1DLinear, Function1D, Function2D, Interpolate2DLinear, AxisymmetricMapper +from cherab.core.math import AxisymmetricMapper from cherab.tools.equilibrium import EFITEquilibrium @@ -626,7 +627,7 @@ def interpolators1d_fractional(atomic_data: AtomicData, element: Element, free_v # use profiles to create interpolators for profiles fractional_interpolators = {} for key, item in fractional_profiles.items(): - fractional_interpolators[key] = Interpolate1DLinear(free_variable, item) + fractional_interpolators[key] = Interpolator1DArray(free_variable, item, 'linear', 'none', 0) return fractional_interpolators @@ -657,7 +658,7 @@ def interpolators2d_fractional(atomic_data: AtomicData, element: Element, free_v # use profiles to create interpolators for profiles fractional_interpolators = {} for key, item in fractional_profiles.items(): - fractional_interpolators[key] = Interpolate2DLinear(*free_variable, item) + fractional_interpolators[key] = Interpolator2DArray(*free_variable, item, 'linear', 'none', 0, 0) return fractional_interpolators @@ -690,7 +691,7 @@ def interpolators1d_from_elementdensity(atomic_data: AtomicData, element: Elemen density_interpolators = {} for key, value in densities.items(): - density_interpolators[key] = Interpolate1DLinear(free_variable, value) + density_interpolators[key] = Interpolator1DArray(free_variable, value, 'linear', 'none', 0) return density_interpolators @@ -725,7 +726,7 @@ def interpolators1d_match_plasma_neutrality(atomic_data: AtomicData, element: El # use profiles to create interpolators for profiles density_interpolators = {} for key, item in density_profiles.items(): - density_interpolators[key] = Interpolate1DLinear(free_variable, item) + density_interpolators[key] = Interpolator1DArray(free_variable, item, 'linear', 'none', 0) return density_interpolators @@ -757,7 +758,7 @@ def interpolators2d_from_elementdensity(atomic_data: AtomicData, element: Elemen density_interpolators = {} for key, value in densities.items(): - density_interpolators[key] = Interpolate2DLinear(*free_variable, value) + density_interpolators[key] = Interpolator2DArray(*free_variable, value, 'linear', 'none', 0, 0) return density_interpolators @@ -792,7 +793,7 @@ def interpolators2d_match_plasma_neutrality(atomic_data: AtomicData, element: El # use profiles to create interpolators for profiles density_interpolators = {} for key, item in density_profiles.items(): - density_interpolators[key] = Interpolate2DLinear(*free_variable, item) + density_interpolators[key] = Interpolator2DArray(*free_variable, item, 'linear', 'none', 0, 0) return density_interpolators diff --git a/cherab/tools/tests/data/atomic_rates_mockup/ionisation/h.json b/cherab/tools/tests/data/atomic_rates_mockup/ionisation/h.json new file mode 100644 index 00000000..abf36b6a --- /dev/null +++ b/cherab/tools/tests/data/atomic_rates_mockup/ionisation/h.json @@ -0,0 +1,79 @@ +{ + "0": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-43, + 1.0e-20, + 1.0e-14, + 3.0e-14, + 2.0e-14, + 1.0e-14 + ], + [ + 1.2e-43, + 1.2e-20, + 1.2e-14, + 3.2e-14, + 2.2e-14, + 1.2e-14 + ], + [ + 1.4e-43, + 1.4e-20, + 1.4e-14, + 3.4e-14, + 2.4e-14, + 1.4e-14 + ], + [ + 1.8e-43, + 1.8e-20, + 1.8e-14, + 3.8e-14, + 2.8e-14, + 1.8e-14 + ], + [ + 2.0e-43, + 2.0e-20, + 2.0e-14, + 4.0e-14, + 3.0e-14, + 2.0e-14 + ], + [ + 2.4e-43, + 2.4e-20, + 2.4e-14, + 4.4e-14, + 3.4e-14, + 2.4e-14 + ], + [ + 2.4e-43, + 2.4e-20, + 2.4e-14, + 4.4e-14, + 3.4e-14, + 2.4e-14 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + } +} diff --git a/cherab/tools/tests/data/atomic_rates_mockup/ionisation/he.json b/cherab/tools/tests/data/atomic_rates_mockup/ionisation/he.json new file mode 100644 index 00000000..b8130e8c --- /dev/null +++ b/cherab/tools/tests/data/atomic_rates_mockup/ionisation/he.json @@ -0,0 +1,156 @@ +{ + "0": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.5e-50, + 1.5e-20, + 1.5e-14, + 2.0e-14, + 2.5e-14, + 2.6e-14 + ], + [ + 2.0e-50, + 2.0e-20, + 2.0e-14, + 2.5e-14, + 3.0e-14, + 3.1e-14 + ], + [ + 2.5e-50, + 2.5e-20, + 2.5e-14, + 3.0e-14, + 3.5e-14, + 3.6e-14 + ], + [ + 3.0e-50, + 3.0e-20, + 3.0e-14, + 3.5e-14, + 4.0e-14, + 4.1e-14 + ], + [ + 3.5e-50, + 3.5e-20, + 3.5e-14, + 4.0e-14, + 4.5e-14, + 4.6e-14 + ], + [ + 4.0e-50, + 4.0e-20, + 4.0e-14, + 4.5e-14, + 5.0e-14, + 5.1e-14 + ], + [ + 4.5e-50, + 4.5e-20, + 4.5e-14, + 5.0e-14, + 5.5e-14, + 5.6e-14 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "1": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.5e-51, + 1.5e-21, + 1.5e-15, + 2.0e-15, + 2.5e-15, + 2.6e-15 + ], + [ + 2.0e-51, + 2.0e-21, + 2.0e-15, + 2.5e-15, + 3.0e-15, + 3.1e-15 + ], + [ + 2.5e-51, + 2.5e-21, + 2.5e-15, + 3.0e-15, + 3.5e-15, + 3.6e-15 + ], + [ + 3.0e-51, + 3.0e-21, + 3.0e-15, + 3.5e-15, + 4.0e-15, + 4.1e-15 + ], + [ + 3.5e-51, + 3.5e-21, + 3.5e-15, + 4.0e-15, + 4.5e-15, + 4.6e-15 + ], + [ + 4.0e-51, + 4.0e-21, + 4.0e-15, + 4.5e-15, + 5.0e-15, + 5.1e-15 + ], + [ + 4.5e-51, + 4.5e-21, + 4.5e-15, + 5.0e-15, + 5.5e-15, + 5.6e-15 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + } +} diff --git a/cherab/tools/tests/data/atomic_rates_mockup/ionisation/ne.json b/cherab/tools/tests/data/atomic_rates_mockup/ionisation/ne.json new file mode 100644 index 00000000..3bdd23a4 --- /dev/null +++ b/cherab/tools/tests/data/atomic_rates_mockup/ionisation/ne.json @@ -0,0 +1,772 @@ +{ + "0": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-40, + 1.0e-18, + 1.0e-14, + 5.0e-14, + 1.0e-13, + 5.0e-13 + ], + [ + 1.2e-40, + 1.2e-18, + 1.2e-14, + 5.2e-14, + 1.2e-13, + 5.2e-13 + ], + [ + 1.4e-40, + 1.4e-18, + 1.4e-14, + 5.4e-14, + 1.4e-13, + 5.4e-13 + ], + [ + 1.6e-40, + 1.6e-18, + 1.6e-14, + 5.6e-14, + 1.6e-13, + 5.6e-13 + ], + [ + 1.7e-40, + 1.7e-18, + 1.7e-14, + 5.7e-14, + 1.7e-13, + 5.7e-13 + ], + [ + 1.8e-40, + 1.8e-18, + 1.8e-14, + 5.8e-14, + 1.8e-13, + 5.8e-13 + ], + [ + 1.9e-40, + 1.9e-18, + 1.9e-14, + 5.9e-14, + 1.9e-13, + 5.9e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "1": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-41, + 1.0e-19, + 1.0e-15, + 5.0e-15, + 1.0e-14, + 5.0e-14 + ], + [ + 1.2e-41, + 1.2e-19, + 1.2e-15, + 5.2e-15, + 1.2e-14, + 5.2e-14 + ], + [ + 1.4e-41, + 1.4e-19, + 1.4e-15, + 5.4e-15, + 1.4e-14, + 5.4e-14 + ], + [ + 1.6e-41, + 1.6e-19, + 1.6e-15, + 5.6e-15, + 1.6e-14, + 5.6e-14 + ], + [ + 1.7e-41, + 1.7e-19, + 1.7e-15, + 5.7e-15, + 1.7e-14, + 5.7e-14 + ], + [ + 1.8e-41, + 1.8e-19, + 1.8e-15, + 5.8e-15, + 1.8e-14, + 5.8e-14 + ], + [ + 1.9e-41, + 1.9e-19, + 1.9e-15, + 5.9e-15, + 1.9e-14, + 5.9e-14 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "2": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 5.0e-42, + 5.0e-20, + 5.0e-16, + 1.0e-15, + 5.0e-15, + 1.0e-14 + ], + [ + 5.2e-42, + 5.2e-20, + 5.2e-16, + 1.2e-15, + 5.2e-15, + 1.2e-14 + ], + [ + 5.4e-42, + 5.4e-20, + 5.4e-16, + 1.4e-15, + 5.4e-15, + 1.4e-14 + ], + [ + 5.6e-42, + 5.6e-20, + 5.6e-16, + 1.6e-15, + 5.6e-15, + 1.6e-14 + ], + [ + 5.7e-42, + 5.7e-20, + 5.7e-16, + 1.7e-15, + 5.7e-15, + 1.7e-14 + ], + [ + 5.8e-42, + 5.8e-20, + 5.8e-16, + 1.8e-15, + 5.8e-15, + 1.8e-14 + ], + [ + 5.9e-42, + 5.9e-20, + 5.9e-16, + 1.9e-15, + 5.9e-15, + 1.9e-14 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "3": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-42, + 1.0e-20, + 1.0e-16, + 5.0e-16, + 1.0e-15, + 5.0e-15 + ], + [ + 1.2e-42, + 1.2e-20, + 1.2e-16, + 5.2e-16, + 1.2e-15, + 5.2e-15 + ], + [ + 1.4e-42, + 1.4e-20, + 1.4e-16, + 5.4e-16, + 1.4e-15, + 5.4e-15 + ], + [ + 1.6e-42, + 1.6e-20, + 1.6e-16, + 5.6e-16, + 1.6e-15, + 5.6e-15 + ], + [ + 1.7e-42, + 1.7e-20, + 1.7e-16, + 5.7e-16, + 1.7e-15, + 5.7e-15 + ], + [ + 1.8e-42, + 1.8e-20, + 1.8e-16, + 5.8e-16, + 1.8e-15, + 5.8e-15 + ], + [ + 1.9e-42, + 1.9e-20, + 1.9e-16, + 5.9e-16, + 1.9e-15, + 5.9e-15 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "4": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 5.0e-43, + 5.0e-21, + 5.0e-17, + 1.0e-16, + 5.0e-16, + 1.0e-17 + ], + [ + 5.2e-43, + 5.2e-21, + 5.2e-17, + 1.2e-16, + 5.2e-16, + 1.2e-15 + ], + [ + 5.4e-43, + 5.4e-21, + 5.4e-17, + 1.4e-16, + 5.4e-16, + 1.4e-15 + ], + [ + 5.6e-43, + 5.6e-21, + 5.6e-17, + 1.6e-16, + 5.6e-16, + 1.6e-15 + ], + [ + 5.7e-43, + 5.7e-21, + 5.7e-17, + 1.7e-16, + 5.7e-16, + 1.7e-15 + ], + [ + 5.8e-43, + 5.8e-21, + 5.8e-17, + 1.8e-16, + 5.8e-16, + 1.8e-15 + ], + [ + 5.9e-43, + 5.9e-21, + 5.9e-17, + 1.9e-16, + 5.9e-16, + 1.9e-15 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "5": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-43, + 1.0e-21, + 1.0e-17, + 5.0e-17, + 1.0e-16, + 5.0e-16 + ], + [ + 1.2e-43, + 1.2e-21, + 1.2e-17, + 5.2e-17, + 1.2e-16, + 5.2e-16 + ], + [ + 1.4e-43, + 1.4e-21, + 1.4e-17, + 5.4e-17, + 1.4e-16, + 5.4e-16 + ], + [ + 1.6e-43, + 1.6e-21, + 1.6e-17, + 5.6e-17, + 1.6e-16, + 5.6e-16 + ], + [ + 1.7e-43, + 1.7e-21, + 1.7e-17, + 5.7e-17, + 1.7e-16, + 5.7e-16 + ], + [ + 1.8e-43, + 1.8e-21, + 1.8e-17, + 5.8e-17, + 1.8e-16, + 5.8e-16 + ], + [ + 1.9e-43, + 1.9e-21, + 1.9e-17, + 5.9e-17, + 1.9e-16, + 5.9e-16 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "6": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 5.0e-44, + 5.0e-22, + 5.0e-18, + 1.0e-17, + 5.0e-17, + 1.0e-16 + ], + [ + 5.2e-44, + 5.2e-22, + 5.2e-18, + 1.2e-17, + 5.2e-17, + 1.2e-16 + ], + [ + 5.4e-44, + 5.4e-22, + 5.4e-18, + 1.4e-17, + 5.4e-17, + 1.4e-16 + ], + [ + 5.6e-44, + 5.6e-22, + 5.6e-18, + 1.6e-17, + 5.6e-17, + 1.6e-16 + ], + [ + 5.7e-44, + 5.7e-22, + 5.7e-18, + 1.7e-17, + 5.7e-17, + 1.7e-16 + ], + [ + 5.8e-44, + 5.8e-22, + 5.8e-18, + 1.8e-17, + 5.8e-17, + 1.8e-16 + ], + [ + 5.9e-44, + 5.9e-22, + 5.9e-18, + 1.9e-17, + 5.9e-17, + 1.9e-16 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "7": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-44, + 1.0e-22, + 1.0e-18, + 5.0e-18, + 1.0e-17, + 5.0e-17 + ], + [ + 1.2e-44, + 1.2e-22, + 1.2e-18, + 5.2e-18, + 1.2e-17, + 5.2e-17 + ], + [ + 1.4e-44, + 1.4e-22, + 1.4e-18, + 5.4e-18, + 1.4e-17, + 5.4e-17 + ], + [ + 1.6e-44, + 1.6e-22, + 1.6e-18, + 5.6e-18, + 1.6e-17, + 5.6e-17 + ], + [ + 1.7e-44, + 1.7e-22, + 1.7e-18, + 5.7e-18, + 1.7e-17, + 5.7e-17 + ], + [ + 1.8e-44, + 1.8e-22, + 1.8e-18, + 5.8e-18, + 1.8e-17, + 5.8e-17 + ], + [ + 1.9e-44, + 1.9e-22, + 1.9e-18, + 5.9e-18, + 1.9e-17, + 5.9e-17 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "8": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 5.0e-45, + 5.0e-23, + 5.0e-19, + 1.0e-18, + 5.0e-18, + 1.0e-17 + ], + [ + 5.2e-45, + 5.2e-23, + 5.2e-19, + 1.2e-18, + 5.2e-18, + 1.2e-17 + ], + [ + 5.4e-45, + 5.4e-23, + 5.4e-19, + 1.4e-18, + 5.4e-18, + 1.4e-17 + ], + [ + 5.6e-45, + 5.6e-23, + 5.6e-19, + 1.6e-18, + 5.6e-18, + 1.6e-17 + ], + [ + 5.7e-45, + 5.7e-23, + 5.7e-19, + 1.7e-18, + 5.7e-18, + 1.7e-17 + ], + [ + 5.8e-45, + 5.8e-23, + 5.8e-19, + 1.8e-18, + 5.8e-18, + 1.8e-17 + ], + [ + 5.9e-45, + 5.9e-23, + 5.9e-19, + 1.9e-18, + 5.9e-18, + 1.9e-17 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "9": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-45, + 1.0e-23, + 1.0e-19, + 5.0e-19, + 1.0e-18, + 5.0e-18 + ], + [ + 1.2e-45, + 1.2e-23, + 1.2e-19, + 5.2e-19, + 1.2e-18, + 5.2e-18 + ], + [ + 1.4e-45, + 1.4e-23, + 1.4e-19, + 5.4e-19, + 1.4e-18, + 5.4e-18 + ], + [ + 1.6e-45, + 1.6e-23, + 1.6e-19, + 5.6e-19, + 1.6e-18, + 5.6e-18 + ], + [ + 1.7e-45, + 1.7e-23, + 1.7e-19, + 5.7e-19, + 1.7e-18, + 5.7e-18 + ], + [ + 1.8e-45, + 1.8e-23, + 1.8e-19, + 5.8e-19, + 1.8e-18, + 5.8e-18 + ], + [ + 1.9e-45, + 1.9e-23, + 1.9e-19, + 5.9e-19, + 1.9e-18, + 5.9e-18 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + } +} diff --git a/cherab/tools/tests/data/atomic_rates_mockup/recombination/h.json b/cherab/tools/tests/data/atomic_rates_mockup/recombination/h.json new file mode 100644 index 00000000..8d562332 --- /dev/null +++ b/cherab/tools/tests/data/atomic_rates_mockup/recombination/h.json @@ -0,0 +1,79 @@ +{ + "1": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-18, + 1.0e-19, + 1.0e-20, + 1.0e-21, + 1.0e-22, + 1.0e-23 + ], + [ + 1.5e-18, + 1.5e-19, + 1.5e-20, + 1.5e-21, + 1.5e-22, + 1.5e-23 + ], + [ + 2.0e-18, + 2.0e-19, + 2.0e-20, + 2.0e-21, + 2.0e-22, + 2.0e-23 + ], + [ + 2.5e-18, + 2.5e-19, + 2.5e-20, + 2.5e-21, + 2.5e-22, + 2.5e-23 + ], + [ + 3.0e-18, + 3.0e-19, + 3.0e-20, + 3.0e-21, + 3.0e-22, + 3.0e-23 + ], + [ + 3.5e-18, + 3.5e-19, + 3.5e-20, + 3.5e-21, + 3.5e-22, + 3.5e-23 + ], + [ + 4.5e-18, + 4.5e-19, + 4.5e-20, + 4.5e-21, + 4.5e-22, + 4.5e-23 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + } +} diff --git a/cherab/tools/tests/data/atomic_rates_mockup/recombination/he.json b/cherab/tools/tests/data/atomic_rates_mockup/recombination/he.json new file mode 100644 index 00000000..a48ce01d --- /dev/null +++ b/cherab/tools/tests/data/atomic_rates_mockup/recombination/he.json @@ -0,0 +1,156 @@ +{ + "1": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.2e-18, + 5.2e-19, + 1.2e-19, + 5.2e-20, + 1.2e-20, + 1.2e-22 + ], + [ + 2.2e-18, + 1.2e-18, + 2.2e-19, + 1.2e-19, + 2.2e-20, + 2.2e-22 + ], + [ + 3.2e-18, + 2.2e-18, + 3.2e-19, + 2.2e-19, + 3.2e-20, + 3.2e-22 + ], + [ + 4.2e-18, + 3.2e-18, + 4.2e-19, + 3.2e-19, + 4.2e-20, + 4.2e-22 + ], + [ + 5.2e-18, + 4.2e-18, + 5.2e-19, + 4.2e-19, + 5.2e-20, + 5.2e-22 + ], + [ + 6.2e-18, + 4.2e-18, + 6.2e-19, + 4.2e-19, + 6.2e-20, + 6.2e-22 + ], + [ + 7.2e-18, + 5.2e-18, + 7.2e-19, + 5.2e-19, + 7.2e-20, + 7.2e-22 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, +"2": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.2e-19, + 5.2e-20, + 1.2e-20, + 5.2e-21, + 1.2e-21, + 1.2e-23 + ], + [ + 2.2e-19, + 1.2e-19, + 2.2e-20, + 1.2e-20, + 2.2e-21, + 2.2e-23 + ], + [ + 3.2e-19, + 2.2e-19, + 3.2e-20, + 2.2e-20, + 3.2e-21, + 3.2e-23 + ], + [ + 4.2e-19, + 3.2e-19, + 4.2e-20, + 3.2e-20, + 4.2e-21, + 4.2e-23 + ], + [ + 5.2e-19, + 4.2e-19, + 5.2e-20, + 4.2e-20, + 5.2e-21, + 5.2e-23 + ], + [ + 6.2e-19, + 4.2e-19, + 6.2e-20, + 4.2e-20, + 6.2e-21, + 6.2e-23 + ], + [ + 7.2e-19, + 5.2e-19, + 7.2e-20, + 5.2e-20, + 7.2e-21, + 7.2e-23 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + } +} diff --git a/cherab/tools/tests/data/atomic_rates_mockup/recombination/ne.json b/cherab/tools/tests/data/atomic_rates_mockup/recombination/ne.json new file mode 100644 index 00000000..2452250b --- /dev/null +++ b/cherab/tools/tests/data/atomic_rates_mockup/recombination/ne.json @@ -0,0 +1,772 @@ +{ + "1": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-19, + 5.0e-20, + 1.0e-20, + 5.0e-21, + 1.0e-21, + 1.0e-22 + ], + [ + 1.2e-19, + 5.2e-20, + 1.2e-20, + 5.2e-21, + 1.2e-21, + 1.2e-22 + ], + [ + 1.4e-19, + 5.4e-20, + 1.4e-20, + 5.4e-21, + 1.4e-21, + 1.4e-22 + ], + [ + 1.6e-19, + 5.6e-20, + 1.6e-20, + 5.6e-21, + 1.6e-21, + 1.6e-22 + ], + [ + 1.7e-19, + 5.7e-20, + 1.7e-20, + 5.7e-21, + 1.7e-21, + 1.7e-22 + ], + [ + 1.8e-19, + 5.8e-20, + 1.8e-20, + 5.8e-21, + 1.8e-21, + 1.8e-22 + ], + [ + 1.9e-19, + 5.9e-20, + 1.9e-20, + 5.9e-21, + 1.9e-21, + 1.9e-22 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "2": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 2.0e-19, + 1.0e-19, + 2.0e-20, + 1.0e-20, + 2.0e-21, + 2.0e-22 + ], + [ + 2.2e-19, + 1.1e-19, + 2.2e-20, + 1.1e-20, + 2.2e-21, + 2.2e-22 + ], + [ + 2.5e-19, + 1.2e-19, + 2.5e-20, + 1.2e-20, + 2.5e-21, + 2.5e-22 + ], + [ + 3.0e-19, + 1.5e-19, + 3.0e-20, + 1.5e-20, + 3.0e-21, + 3.0e-22 + ], + [ + 3.5e-19, + 1.7e-19, + 3.5e-20, + 1.7e-20, + 3.5e-21, + 3.5e-22 + ], + [ + 4.0e-19, + 2.0e-19, + 4.0e-20, + 2.0e-20, + 3.0e-21, + 3.0e-22 + ], + [ + 4.5e-19, + 2.2e-19, + 4.5e-20, + 2.2e-20, + 4.5e-21, + 2.2e-22 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "3": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 5.0e-19, + 1.0e-19, + 5.0e-20, + 1.0e-20, + 5.0e-21, + 5.0e-22 + ], + [ + 5.5e-19, + 1.5e-19, + 5.5e-20, + 1.5e-20, + 5.5e-21, + 5.5e-22 + ], + [ + 6.0e-19, + 3.0e-19, + 6.0e-20, + 3.0e-20, + 6.0e-21, + 6.0e-22 + ], + [ + 6.5e-19, + 3.2e-19, + 6.5e-20, + 3.2e-20, + 6.5e-21, + 6.5e-22 + ], + [ + 7.0e-19, + 3.5e-19, + 7.0e-20, + 3.5e-20, + 7.0e-21, + 7.0e-22 + ], + [ + 8.0e-19, + 4.0e-19, + 8.0e-20, + 4.0e-20, + 8.0e-21, + 8.0e-22 + ], + [ + 9.0e-19, + 4.5e-19, + 9.0e-20, + 4.5e-20, + 9.0e-21, + 9.0e-22 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "4": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-18, + 5.0e-19, + 1.0e-19, + 5.0e-20, + 1.0e-20, + 1.0e-21 + ], + [ + 1.2e-18, + 5.2e-19, + 1.2e-19, + 5.2e-20, + 1.2e-20, + 1.2e-21 + ], + [ + 1.4e-19, + 5.4e-19, + 1.4e-19, + 5.4e-20, + 1.4e-20, + 1.4e-21 + ], + [ + 1.6e-18, + 5.6e-19, + 1.6e-19, + 5.6e-20, + 1.6e-20, + 1.6e-21 + ], + [ + 1.7e-18, + 5.7e-19, + 1.7e-19, + 5.7e-20, + 1.7e-20, + 1.7e-21 + ], + [ + 1.8e-18, + 5.8e-19, + 1.8e-19, + 5.8e-20, + 1.8e-20, + 1.8e-21 + ], + [ + 1.9e-18, + 5.9e-19, + 1.9e-19, + 5.9e-20, + 1.9e-20, + 1.9e-21 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "5": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 2.0e-18, + 1.0e-18, + 2.0e-19, + 1.0e-19, + 2.0e-20, + 2.0e-21 + ], + [ + 2.2e-18, + 1.1e-18, + 2.2e-19, + 1.1e-19, + 2.2e-20, + 2.2e-21 + ], + [ + 2.5e-18, + 1.2e-18, + 2.5e-19, + 1.2e-19, + 2.5e-20, + 2.5e-21 + ], + [ + 3.0e-18, + 1.5e-18, + 3.0e-19, + 1.5e-19, + 3.0e-20, + 3.0e-21 + ], + [ + 3.5e-18, + 1.7e-18, + 3.5e-19, + 1.7e-19, + 3.5e-20, + 3.5e-21 + ], + [ + 4.0e-18, + 2.0e-18, + 4.0e-19, + 2.0e-19, + 3.0e-20, + 3.0e-21 + ], + [ + 4.5e-18, + 2.2e-18, + 4.5e-19, + 2.2e-19, + 4.5e-20, + 2.2e-21 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "6": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 5.0e-18, + 1.0e-18, + 5.0e-19, + 1.0e-19, + 5.0e-20, + 5.0e-21 + ], + [ + 5.5e-18, + 1.5e-18, + 5.5e-19, + 1.5e-19, + 5.5e-20, + 5.5e-21 + ], + [ + 6.0e-18, + 3.0e-18, + 6.0e-19, + 3.0e-19, + 6.0e-20, + 6.0e-21 + ], + [ + 6.5e-18, + 3.2e-18, + 6.5e-19, + 3.2e-19, + 6.5e-20, + 6.5e-21 + ], + [ + 7.0e-18, + 3.5e-18, + 7.0e-19, + 3.5e-19, + 7.0e-20, + 7.0e-21 + ], + [ + 8.0e-18, + 4.0e-18, + 8.0e-19, + 4.0e-19, + 8.0e-20, + 8.0e-21 + ], + [ + 9.0e-18, + 4.5e-18, + 9.0e-19, + 4.5e-19, + 9.0e-20, + 9.0e-21 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "7": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-17, + 5.0e-18, + 1.0e-18, + 5.0e-19, + 1.0e-19, + 1.0e-20 + ], + [ + 1.2e-17, + 5.2e-18, + 1.2e-18, + 5.2e-19, + 1.2e-19, + 1.2e-20 + ], + [ + 1.4e-17, + 5.4e-18, + 1.4e-18, + 5.4e-19, + 1.4e-19, + 1.4e-20 + ], + [ + 1.6e-17, + 5.6e-18, + 1.6e-18, + 5.6e-19, + 1.6e-19, + 1.6e-20 + ], + [ + 1.7e-17, + 5.7e-18, + 1.7e-18, + 5.7e-19, + 1.7e-19, + 1.7e-20 + ], + [ + 1.8e-17, + 5.8e-18, + 1.8e-18, + 5.8e-19, + 1.8e-19, + 1.8e-20 + ], + [ + 1.9e-17, + 5.9e-18, + 1.9e-18, + 5.9e-19, + 1.9e-19, + 1.9e-20 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "8": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 2.0e-17, + 1.0e-17, + 2.0e-18, + 1.0e-18, + 2.0e-19, + 2.0e-20 + ], + [ + 2.2e-17, + 1.1e-17, + 2.2e-18, + 1.1e-18, + 2.2e-19, + 2.2e-20 + ], + [ + 2.5e-17, + 1.2e-17, + 2.5e-18, + 1.2e-18, + 2.5e-19, + 2.5e-20 + ], + [ + 3.0e-17, + 1.5e-17, + 3.0e-18, + 1.5e-18, + 3.0e-19, + 3.0e-20 + ], + [ + 3.5e-17, + 1.7e-17, + 3.5e-18, + 1.7e-18, + 3.5e-19, + 3.5e-20 + ], + [ + 4.0e-17, + 2.0e-17, + 4.0e-18, + 2.0e-18, + 3.0e-19, + 3.0e-20 + ], + [ + 4.5e-17, + 2.2e-17, + 4.5e-18, + 2.2e-18, + 4.5e-19, + 2.2e-20 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "9": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 5.0e-17, + 1.0e-17, + 5.0e-18, + 1.0e-18, + 5.0e-19, + 5.0e-20 + ], + [ + 5.5e-17, + 1.5e-17, + 5.5e-18, + 1.5e-18, + 5.5e-19, + 5.5e-20 + ], + [ + 6.0e-17, + 3.0e-17, + 6.0e-18, + 3.0e-18, + 6.0e-19, + 6.0e-20 + ], + [ + 6.5e-17, + 3.2e-17, + 6.5e-18, + 3.2e-18, + 6.5e-19, + 6.5e-20 + ], + [ + 7.0e-17, + 3.5e-17, + 7.0e-18, + 3.5e-18, + 7.0e-19, + 7.0e-20 + ], + [ + 8.0e-17, + 4.0e-17, + 8.0e-18, + 4.0e-18, + 8.0e-19, + 8.0e-20 + ], + [ + 9.0e-17, + 4.5e-17, + 9.0e-18, + 4.5e-18, + 9.0e-19, + 9.0e-20 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "10": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-16, + 5.0e-17, + 1.0e-17, + 5.0e-18, + 1.0e-18, + 1.0e-19 + ], + [ + 1.2e-16, + 5.2e-17, + 1.2e-17, + 5.2e-18, + 1.2e-18, + 1.2e-19 + ], + [ + 1.4e-16, + 5.4e-17, + 1.4e-17, + 5.4e-18, + 1.4e-18, + 1.4e-19 + ], + [ + 1.6e-16, + 5.6e-17, + 1.6e-17, + 5.6e-18, + 1.6e-18, + 1.6e-19 + ], + [ + 1.7e-16, + 5.7e-17, + 1.7e-17, + 5.7e-18, + 1.7e-18, + 1.7e-19 + ], + [ + 1.8e-16, + 5.8e-17, + 1.8e-17, + 5.8e-18, + 1.8e-18, + 1.8e-19 + ], + [ + 1.9e-16, + 5.9e-17, + 1.9e-17, + 5.9e-18, + 1.9e-18, + 1.9e-19 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + } +} diff --git a/cherab/tools/tests/data/atomic_rates_mockup/thermal_cx/h/0/h.json b/cherab/tools/tests/data/atomic_rates_mockup/thermal_cx/h/0/h.json new file mode 100644 index 00000000..52614635 --- /dev/null +++ b/cherab/tools/tests/data/atomic_rates_mockup/thermal_cx/h/0/h.json @@ -0,0 +1,79 @@ +{ + "1": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 5.0e-15, + 7.5e-15, + 1.0e-14, + 2.5e-14, + 5.0e-14, + 1.0e-13 + ], + [ + 5.0e-15, + 7.5e-15, + 1.0e-14, + 2.5e-14, + 5.0e-14, + 1.0e-13 + ], + [ + 5.0e-15, + 7.5e-15, + 1.0e-14, + 2.5e-14, + 5.0e-14, + 1.0e-13 + ], + [ + 5.0e-15, + 7.5e-15, + 1.0e-14, + 2.5e-14, + 5.0e-14, + 1.0e-13 + ], + [ + 5.0e-15, + 7.5e-15, + 1.0e-14, + 2.5e-14, + 5.0e-14, + 1.0e-13 + ], + [ + 5.0e-15, + 7.5e-15, + 1.0e-14, + 2.5e-14, + 5.0e-14, + 1.0e-13 + ], + [ + 5.0e-15, + 7.5e-15, + 1.0e-14, + 2.5e-14, + 5.0e-14, + 1.0e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + } +} diff --git a/cherab/tools/tests/data/atomic_rates_mockup/thermal_cx/h/0/he.json b/cherab/tools/tests/data/atomic_rates_mockup/thermal_cx/h/0/he.json new file mode 100644 index 00000000..7db581d2 --- /dev/null +++ b/cherab/tools/tests/data/atomic_rates_mockup/thermal_cx/h/0/he.json @@ -0,0 +1,156 @@ +{ + "1": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80 + ], + [ + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80 + ], + [ + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80 + ], + [ + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80 + ], + [ + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80 + ], + [ + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80 + ], + [ + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80, + 2.0e-80 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, +"2": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80 + ], + [ + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80 + ], + [ + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80 + ], + [ + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80 + ], + [ + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80 + ], + [ + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80 + ], + [ + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80, + 1.0e-80 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + } +} diff --git a/cherab/tools/tests/data/atomic_rates_mockup/thermal_cx/h/0/ne.json b/cherab/tools/tests/data/atomic_rates_mockup/thermal_cx/h/0/ne.json new file mode 100644 index 00000000..f0dcc81a --- /dev/null +++ b/cherab/tools/tests/data/atomic_rates_mockup/thermal_cx/h/0/ne.json @@ -0,0 +1,772 @@ +{ + "1": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-16, + 5.0e-16, + 1.0e-15, + 5.0e-15, + 1.0e-14, + 1.0e-13 + ], + [ + 1.0e-16, + 5.0e-16, + 1.0e-15, + 5.0e-15, + 1.0e-14, + 1.0e-13 + ], + [ + 1.0e-16, + 5.0e-16, + 1.0e-15, + 5.0e-15, + 1.0e-14, + 1.0e-13 + ], + [ + 1.0e-16, + 5.0e-16, + 1.0e-15, + 5.0e-15, + 1.0e-14, + 1.0e-13 + ], + [ + 1.0e-16, + 5.0e-16, + 1.0e-15, + 5.0e-15, + 1.0e-14, + 1.0e-13 + ], + [ + 1.0e-16, + 5.0e-16, + 1.0e-15, + 5.0e-15, + 1.0e-14, + 1.0e-13 + ], + [ + 1.0e-16, + 5.0e-16, + 1.0e-15, + 5.0e-15, + 1.0e-14, + 1.0e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "2": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 2.0e-16, + 1.0e-15, + 2.0e-15, + 1.0e-14, + 2.0e-14, + 2.0e-13 + ], + [ + 2.0e-16, + 1.0e-15, + 2.0e-15, + 1.0e-14, + 2.0e-14, + 2.0e-13 + ], + [ + 2.0e-16, + 1.0e-15, + 2.0e-15, + 1.0e-14, + 2.0e-14, + 2.0e-13 + ], + [ + 2.0e-16, + 1.0e-15, + 2.0e-15, + 1.0e-14, + 2.0e-14, + 2.0e-13 + ], + [ + 2.0e-16, + 1.0e-15, + 2.0e-15, + 1.0e-14, + 2.0e-14, + 2.0e-13 + ], + [ + 2.0e-16, + 1.0e-15, + 2.0e-15, + 1.0e-14, + 2.0e-14, + 2.0e-13 + ], + [ + 2.0e-16, + 1.0e-15, + 2.0e-15, + 1.0e-14, + 2.0e-14, + 2.0e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "3": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 3.0e-16, + 1.5e-16, + 3.0e-15, + 1.5e-14, + 3.0e-14, + 3.0e-13 + ], + [ + 3.0e-16, + 1.5e-16, + 3.0e-15, + 1.5e-14, + 3.0e-14, + 3.0e-13 + ], + [ + 3.0e-16, + 1.5e-16, + 3.0e-15, + 1.5e-14, + 3.0e-14, + 3.0e-13 + ], + [ + 3.0e-16, + 1.5e-16, + 3.0e-15, + 1.5e-14, + 3.0e-14, + 3.0e-13 + ], + [ + 3.0e-16, + 1.5e-16, + 3.0e-15, + 1.5e-14, + 3.0e-14, + 3.0e-13 + ], + [ + 3.0e-16, + 1.5e-16, + 3.0e-15, + 1.5e-14, + 3.0e-14, + 3.0e-13 + ], + [ + 3.0e-16, + 1.5e-16, + 3.0e-15, + 1.5e-14, + 3.0e-14, + 3.0e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "4": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 4.0e-16, + 2.0e-15, + 4.0e-15, + 2.0e-14, + 4.0e-14, + 4.0e-13 + ], + [ + 4.0e-16, + 2.0e-15, + 4.0e-15, + 2.0e-14, + 4.0e-14, + 4.0e-13 + ], + [ + 4.0e-16, + 2.0e-15, + 4.0e-15, + 2.0e-14, + 4.0e-14, + 4.0e-13 + ], + [ + 4.0e-16, + 2.0e-15, + 4.0e-15, + 2.0e-14, + 4.0e-14, + 4.0e-13 + ], + [ + 4.0e-16, + 2.0e-15, + 4.0e-15, + 2.0e-14, + 4.0e-14, + 4.0e-13 + ], + [ + 4.0e-16, + 2.0e-15, + 4.0e-15, + 2.0e-14, + 4.0e-14, + 4.0e-13 + ], + [ + 4.0e-16, + 2.0e-15, + 4.0e-15, + 2.0e-14, + 4.0e-14, + 4.0e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "5": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 5.0e-16, + 2.5e-15, + 5.0e-15, + 2.5e-14, + 5.0e-14, + 5.0e-13 + ], + [ + 5.0e-16, + 2.5e-15, + 5.0e-15, + 2.5e-14, + 5.0e-14, + 5.0e-13 + ], + [ + 5.0e-16, + 2.5e-15, + 5.0e-15, + 2.5e-14, + 5.0e-14, + 5.0e-13 + ], + [ + 5.0e-16, + 2.5e-15, + 5.0e-15, + 2.5e-14, + 5.0e-14, + 5.0e-13 + ], + [ + 5.0e-16, + 2.5e-15, + 5.0e-15, + 2.5e-14, + 5.0e-14, + 5.0e-13 + ], + [ + 5.0e-16, + 2.5e-15, + 5.0e-15, + 2.5e-14, + 5.0e-14, + 5.0e-13 + ], + [ + 5.0e-16, + 2.5e-15, + 5.0e-15, + 2.5e-14, + 5.0e-14, + 5.0e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "6": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 6.0e-16, + 3.0e-15, + 6.0e-15, + 3.0e-14, + 6.0e-14, + 6.0e-13 + ], + [ + 6.0e-16, + 3.0e-15, + 6.0e-15, + 3.0e-14, + 6.0e-14, + 6.0e-13 + ], + [ + 6.0e-16, + 3.0e-15, + 6.0e-15, + 3.0e-14, + 6.0e-14, + 6.0e-13 + ], + [ + 6.0e-16, + 3.0e-15, + 6.0e-15, + 3.0e-14, + 6.0e-14, + 6.0e-13 + ], + [ + 6.0e-16, + 3.0e-15, + 6.0e-15, + 3.0e-14, + 6.0e-14, + 6.0e-13 + ], + [ + 6.0e-16, + 3.0e-15, + 6.0e-15, + 3.0e-14, + 6.0e-14, + 6.0e-13 + ], + [ + 6.0e-16, + 3.0e-15, + 6.0e-15, + 3.0e-14, + 6.0e-14, + 6.0e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "7": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 7.0e-16, + 3.5e-15, + 7.0e-15, + 3.5e-14, + 7.0e-14, + 7.0e-13 + ], + [ + 7.0e-16, + 3.5e-15, + 7.0e-15, + 3.5e-14, + 7.0e-14, + 7.0e-13 + ], + [ + 7.0e-16, + 3.5e-15, + 7.0e-15, + 3.5e-14, + 7.0e-14, + 7.0e-13 + ], + [ + 7.0e-16, + 3.5e-15, + 7.0e-15, + 3.5e-14, + 7.0e-14, + 7.0e-13 + ], + [ + 7.0e-16, + 3.5e-15, + 7.0e-15, + 3.5e-14, + 7.0e-14, + 7.0e-13 + ], + [ + 7.0e-16, + 3.5e-15, + 7.0e-15, + 3.5e-14, + 7.0e-14, + 7.0e-13 + ], + [ + 7.0e-16, + 3.5e-15, + 7.0e-15, + 3.5e-14, + 7.0e-14, + 7.0e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "8": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 8.0e-16, + 4.0e-15, + 8.0e-15, + 4.0e-14, + 8.0e-14, + 8.0e-13 + ], + [ + 8.0e-16, + 4.0e-15, + 8.0e-15, + 4.0e-14, + 8.0e-14, + 8.0e-13 + ], + [ + 8.0e-16, + 4.0e-15, + 8.0e-15, + 4.0e-14, + 8.0e-14, + 8.0e-13 + ], + [ + 8.0e-16, + 4.0e-15, + 8.0e-15, + 4.0e-14, + 8.0e-14, + 8.0e-13 + ], + [ + 8.0e-16, + 4.0e-15, + 8.0e-15, + 4.0e-14, + 8.0e-14, + 8.0e-13 + ], + [ + 8.0e-16, + 4.0e-15, + 8.0e-15, + 4.0e-14, + 8.0e-14, + 8.0e-13 + ], + [ + 8.0e-16, + 4.0e-15, + 8.0e-15, + 4.0e-14, + 8.0e-14, + 8.0e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "9": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 9.0e-16, + 4.5e-15, + 9.0e-15, + 4.5e-14, + 9.0e-14, + 9.0e-13 + ], + [ + 9.0e-16, + 4.5e-15, + 9.0e-15, + 4.5e-14, + 9.0e-14, + 9.0e-13 + ], + [ + 9.0e-16, + 4.5e-15, + 9.0e-15, + 4.5e-14, + 9.0e-14, + 9.0e-13 + ], + [ + 9.0e-16, + 4.5e-15, + 9.0e-15, + 4.5e-14, + 9.0e-14, + 9.0e-13 + ], + [ + 9.0e-16, + 4.5e-15, + 9.0e-15, + 4.5e-14, + 9.0e-14, + 9.0e-13 + ], + [ + 9.0e-16, + 4.5e-15, + 9.0e-15, + 4.5e-14, + 9.0e-14, + 9.0e-13 + ], + [ + 9.0e-16, + 4.5e-15, + 9.0e-15, + 4.5e-14, + 9.0e-14, + 9.0e-13 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + }, + "10": { + "ne": [ + 1.0e+15, + 1.0e+16, + 1.0e+17, + 1.0e+18, + 1.0e+19, + 1.0e+20, + 1.0e+21 + ], + "rate": [ + [ + 1.0e-15, + 5.0e-15, + 1.0e-14, + 5.0e-14, + 1.0e-13, + 1.0e-12 + ], + [ + 1.0e-15, + 5.0e-15, + 1.0e-14, + 5.0e-14, + 1.0e-13, + 1.0e-12 + ], + [ + 1.0e-15, + 5.0e-15, + 1.0e-14, + 5.0e-14, + 1.0e-13, + 1.0e-12 + ], + [ + 1.0e-15, + 5.0e-15, + 1.0e-14, + 5.0e-14, + 1.0e-13, + 1.0e-12 + ], + [ + 1.0e-15, + 5.0e-15, + 1.0e-14, + 5.0e-14, + 1.0e-13, + 1.0e-12 + ], + [ + 1.0e-15, + 5.0e-15, + 1.0e-14, + 5.0e-14, + 1.0e-13, + 1.0e-12 + ], + [ + 1.0e-15, + 5.0e-15, + 1.0e-14, + 5.0e-14, + 1.0e-13, + 1.0e-12 + ] + ], + "te": [ + 0.1, + 1, + 10.0, + 100.0, + 1000.0, + 10000.0 + ] + } +} diff --git a/cherab/tools/tests/test_ionization_balance.py b/cherab/tools/tests/test_ionization_balance.py index 2fa65f22..77d4ba72 100644 --- a/cherab/tools/tests/test_ionization_balance.py +++ b/cherab/tools/tests/test_ionization_balance.py @@ -1,985 +1,989 @@ -# -# # 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 unittest -# from collections.abc import Iterable -# import numpy as np -# -# from cherab.core.atomic import neon, hydrogen, helium -# from cherab.core.math import Interpolate1DCubic, Interpolate2DCubic, Function1D, Function2D, AxisymmetricMapper -# from cherab.openadas import OpenADAS -# from cherab.tools.plasmas.ionisation_balance import (fractional_abundance, from_elementdensity, match_plasma_neutrality, -# interpolators1d_fractional, interpolators1d_from_elementdensity, -# interpolators1d_match_plasma_neutrality, -# interpolators2d_fractional, interpolators2d_from_elementdensity, -# interpolators2d_match_plasma_neutrality, -# abundance_axisymmetric_mapper, -# equilibrium_map3d_fractional, equilibrium_map3d_from_elementdensity, -# equilibrium_map3d_match_plasma_neutrality) -# -# from cherab.tools.equilibrium import example_equilibrium -# -# -# def double_parabola(r, centre, edge, p, q): -# return (centre - edge) * np.power((1 - np.power((r - r.min()) / (r.max() - r.min()), p)), q) + edge -# -# -# def normal(x, mu, sd, height=1, offset=0): -# return height * np.exp(-1 * np.power(x - mu, 2) / (2 * sd ** 2)) + offset -# -# -# def get_electron_density_profile(abundances): -# n_e = np.zeros((abundances[0].shape[1])) -# for abundance in abundances: -# for rownumber, row in enumerate(abundance.T): -# n_e[rownumber] += np.sum(row * np.arange(row.shape[0])) -# -# return n_e -# -# -# def get_electron_density_spot(densities): -# n_e = 0 -# for spec in densities: -# for index, value in enumerate(spec): -# n_e += index * value -# -# return n_e -# -# -# def exp_decay(r, lamb, max_val): -# return max_val * np.exp((r - r.max()) * lamb) -# -# -# class TestIonizationBalance1D(unittest.TestCase): -# -# # create plasma profiles and interpolators -# # 1d profiles -# psin_1d = np.linspace(0, 1.1, 15, endpoint=True) -# psin_1d_detailed = np.linspace(0, 1.1, 50, endpoint=True) -# -# t_e_profile_1d = double_parabola(psin_1d, 5000, 10, 2, 2) -# n_e_profile_1d = double_parabola(psin_1d, 6e19, 5e18, 2, 2) -# -# t_element_profile_1d = double_parabola(psin_1d, 1500, 40, 2, 2) -# n_element_profile_1d = double_parabola(psin_1d, 1e17, 1e17, 2, 2) + normal(psin_1d, 0.9, 0.1, 5e17) -# n_element2_profile_1d = double_parabola(psin_1d, 5e17, 1e17, 2, 2) -# -# n_tcx_donor_profile_1d = exp_decay(psin_1d, 10, 3e16) -# -# t_e_1d = Interpolate1DCubic(psin_1d, t_e_profile_1d) -# n_e_1d = Interpolate1DCubic(psin_1d, n_e_profile_1d) -# -# t_element_1d = Interpolate1DCubic(psin_1d, t_element_profile_1d) -# n_element_1d = Interpolate1DCubic(psin_1d, n_element_profile_1d) -# n_element2_1d = Interpolate1DCubic(psin_1d, n_element2_profile_1d) -# -# n_tcx_donor_1d = Interpolate1DCubic(psin_1d, n_tcx_donor_profile_1d) -# -# # denser psi array to test interpolators -# n_e_profile_detailed_1d = np.zeros_like(psin_1d_detailed) -# for index, value in enumerate(psin_1d_detailed): -# n_e_profile_detailed_1d[index] = n_e_1d(value) -# -# # 2d profiles -# -# r = np.linspace(-1, 1, 6) -# z = np.linspace(-1, 1, 8) -# -# psin_2d = np.zeros((*r.shape, *z.shape)) -# -# for index0, value0 in enumerate(r): -# for index1, value1 in enumerate(z): -# if np.sqrt(value0 ** 2 + value1 ** 2) < psin_1d.max(): -# psin_2d[index0, index1] = np.sqrt(value0 ** 2 + value1 ** 2) -# else: -# psin_2d[index0, index1] = psin_1d.max() -# -# t_e_profile_2d = np.zeros_like(psin_2d) -# for index in np.ndindex(*t_e_profile_2d.shape): -# t_e_profile_2d[index] = t_e_1d(psin_2d[index]) -# -# t_e_2d = Interpolate2DCubic(r, z, t_e_profile_2d) -# -# n_e_profile_2d = np.zeros_like(psin_2d) -# for index in np.ndindex(*n_e_profile_2d.shape): -# n_e_profile_2d[index] = n_e_1d(psin_2d[index]) -# -# n_e_2d = Interpolate2DCubic(r, z, n_e_profile_2d) -# -# t_element_profile_2d = np.zeros_like(psin_2d) -# for index in np.ndindex(*t_element_profile_2d.shape): -# t_element_profile_2d[index] = t_element_1d(psin_2d[index]) -# -# t_element_2d = Interpolate2DCubic(r, z, t_element_profile_2d) -# -# n_element_profile_2d = np.zeros_like(psin_2d) -# for index in np.ndindex(*n_element_profile_2d.shape): -# n_element_profile_2d[index] = n_element_1d(psin_2d[index]) -# -# n_element_2d = Interpolate2DCubic(r, z, n_element_profile_2d) -# -# n_element2_profile_2d = np.zeros_like(psin_2d) -# for index in np.ndindex(*n_element2_profile_2d.shape): -# n_element2_profile_2d[index] = n_element2_1d(psin_2d[index]) -# -# n_element2_2d = Interpolate2DCubic(r, z, n_element2_profile_2d) -# -# n_tcx_donor_profile_2d = np.zeros_like(psin_2d) -# for index in np.ndindex(*n_tcx_donor_profile_2d.shape): -# n_tcx_donor_profile_2d[index] = n_tcx_donor_1d(psin_2d[index]) -# -# n_tcx_donor_2d = Interpolate2DCubic(r, z, n_element_profile_2d) -# -# # define psi for single-point tests -# psi_value = 0.9 -# -# # load adas atomic database and define elements -# atomic_data = OpenADAS(permit_extrapolation=True) -# -# element = neon -# element2 = helium -# element_bulk = hydrogen -# -# tcx_donor = hydrogen -# -# TOLERANCE = 1e-3 -# -# def sumup_fractions(self, fractions): -# -# if isinstance(fractions, dict): -# if isinstance(fractions[0], Iterable): -# total = np.zeros_like(fractions[0]) -# else: -# total = 0 -# -# for index, values in fractions.items(): -# total += values -# -# elif isinstance(fractions, np.ndarray): -# total = np.zeros_like(fractions[0, ...]) -# -# for index in np.ndindex(fractions.shape): -# total[index] += fractions[index] -# -# return total -# -# def sumup_electrons(self, densities): -# if isinstance(densities, dict): -# if isinstance(densities[0], Iterable): -# total = np.zeros_like(densities[0]) -# else: -# total = 0 -# -# for index, values in densities.items(): -# total += values * index -# -# elif isinstance(densities, np.ndarray): -# total = np.zeros_like(densities[0, ...]) -# -# for index in np.ndindex(densities.shape): -# total[index] += densities[index] * index[0] -# -# return total -# -# def evaluate_interpolators(self, interpolators, free_variable): -# -# profiles = {} -# for key, item in interpolators.items(): -# if isinstance(item, Function1D): -# profiles[key] = np.zeros_like(free_variable) -# for index in np.ndindex(*free_variable.shape): -# profiles[key][index] = item(free_variable[index]) -# elif isinstance(item, Function2D): -# profiles[key] = np.zeros((*free_variable[0].shape, *free_variable[1].shape)) -# for index0, value0 in enumerate(free_variable[0]): -# for index1, value1 in enumerate(free_variable[1]): -# profiles[key][index0, index1] = item(value0, value1) -# elif isinstance(item, AxisymmetricMapper): -# profiles[key] = np.zeros_like(free_variable) -# for index, value in enumerate(free_variable): -# profiles[key][index] = item(value, 0, 0) -# -# return profiles -# -# def test_fractional_0d_from_0d(self): -# """ -# test fractional abundance calculation with float numbers as inputs -# """ -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d(self.psi_value), -# self.t_e_1d(self.psi_value)) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(1 - fraction_sum < self.TOLERANCE) -# -# def test_fractional_0d_from_0d_tcx(self): -# """ -# test fractional abundance calculation with thermal cx and float numbers as inputs -# """ -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d(self.psi_value), -# self.t_e_1d(self.psi_value), -# self.tcx_donor, self.n_tcx_donor_1d(self.psi_value), 0) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(1 - fraction_sum < self.TOLERANCE) -# -# def test_fractional_0d_from_interpolators(self): -# """ -# test interpolators and free_variable as inputs -# """ -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_1d, -# free_variable=self.psi_value) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(1 - fraction_sum < self.TOLERANCE) -# -# def test_fractional_0d_from_interpolators_tcx(self): -# """ -# test interpolators and free_variable as inputs -# """ -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_1d, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0, free_variable=self.psi_value) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(1 - fraction_sum < self.TOLERANCE) -# -# def test_fractional_0d_from_mixed(self): -# """ -# test mixed types of inputs -# """ -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d(self.psi_value), -# self.t_e_1d, -# free_variable=self.psi_value) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(1 - fraction_sum < self.TOLERANCE) -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, -# self.t_e_1d(self.psi_value), -# free_variable=self.psi_value) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(1 - fraction_sum < self.TOLERANCE) -# -# def test_fractional_1d_from_1d(self): -# """ -# test calculation of 1d fractional profiles with 1d iterables as inputs -# """ -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_1d, -# self.t_e_profile_1d) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# def test_fractional_1d_from_1d_tcx(self): -# """ -# test calculation of 1d fractional profiles with 1d iterables as inputs -# """ -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_1d, -# self.t_e_profile_1d, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_profile_1d, -# tcx_donor_charge=0) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# def test_fractional_1d_from_interpolators(self): -# """ -# test calculation of 1d fractional profiles with 1d interpolators as inputs -# """ -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_1d, -# free_variable=self.psin_1d) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# def test_fractional_1d_from_interpolators_tcx(self): -# """ -# test calculation of 1d fractional profiles with 1d iterables as inputs -# """ -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_1d, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0, free_variable=self.psin_1d) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# def test_fractional_from_mixed(self): -# """ -# test calculation of 1d fractional profiles with mixed types as inputs -# """ -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_1d, self.t_e_1d, -# free_variable=self.psin_1d) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_profile_1d, -# free_variable=self.psin_1d) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# def test_fractional_from_mixed_tcx(self): -# """ -# test calculation of 1d fractional profiles with mixed types as inputs -# """ -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_1d, self.t_e_1d, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0, free_variable=self.psin_1d) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_profile_1d, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0, free_variable=self.psin_1d) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_profile_1d, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_profile_1d, -# tcx_donor_charge=0, free_variable=self.psin_1d) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# def test_fractional_inetrpolators_1d(self): -# """ -# test calculation of 1d fractional interpolators -# """ -# -# interpolators_fractional = interpolators1d_fractional(self.atomic_data, self.element, self.psin_1d, -# self.n_e_profile_1d, self.t_e_1d) -# -# profiles = self.evaluate_interpolators(interpolators_fractional, self.psin_1d) -# fraction_sum = self.sumup_fractions(profiles) -# -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# def test_fractional_inetrpolators_1d_tcx(self): -# """ -# test calculation of 1d fractional interpolators with thermal cx -# """ -# -# interpolators_fractional = interpolators1d_fractional(self.atomic_data, self.element, self.psin_1d, -# self.n_e_profile_1d, self.t_e_1d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_profile_1d, -# tcx_donor_charge=0) -# -# profiles = self.evaluate_interpolators(interpolators_fractional, self.psin_1d) -# fraction_sum = self.sumup_fractions(profiles) -# -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# def test_balance_0d_elementdensity(self): -# """ -# test calculation of ionization balance -# """ -# # test with floats as input -# densities = from_elementdensity(self.atomic_data, self.element, self.n_element_1d(self.psi_value), -# self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value)) -# total = self.sumup_fractions(densities) -# self.assertTrue(np.isclose(total, self.n_element_1d(self.psi_value), rtol=self.TOLERANCE)) -# -# # test with interpolators -# densities = from_elementdensity(self.atomic_data, self.element, self.n_element_1d(self.psi_value), -# self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value)) -# total = self.sumup_fractions(densities) -# self.assertTrue(np.isclose(total, self.n_element_1d(self.psi_value), rtol=self.TOLERANCE)) -# -# # test with mixed parameters -# densities = from_elementdensity(self.atomic_data, self.element, self.n_element_1d, -# self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value), -# free_variable=self.psi_value) -# total = self.sumup_fractions(densities) -# self.assertTrue(np.isclose(total, self.n_element_1d(self.psi_value), rtol=self.TOLERANCE)) -# -# def test_balance_0d_plasma_neutrality(self): -# """test matching of plasma neutrality""" -# -# densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_1d, -# self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value), -# free_variable=self.psi_value) -# -# densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_1d, -# self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value), -# free_variable=self.psi_value) -# -# densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], -# self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value)) -# -# total = self.sumup_electrons(densities_1) -# total += self.sumup_electrons(densities_2) -# total += self.sumup_electrons(densities_3) -# -# self.assertTrue(np.isclose(total, self.n_e_1d(self.psi_value), rtol=self.TOLERANCE)) -# -# def test_balance_0d_plasma_neutrality_tcx(self): -# """test matching of plasma neutrality""" -# -# densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_1d, -# self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value), -# free_variable=self.psi_value, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, tcx_donor_charge=0, -# ) -# -# densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_1d, -# self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value), -# free_variable=self.psi_value, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, tcx_donor_charge=0, -# ) -# -# densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], -# self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value)) -# -# total = self.sumup_electrons(densities_1) -# total += self.sumup_electrons(densities_2) -# total += self.sumup_electrons(densities_3) -# -# self.assertTrue(np.isclose(total, self.n_e_1d(self.psi_value), rtol=self.TOLERANCE)) -# -# def test_balance_1d_plasma_neutrality(self): -# """test matching of plasma neutrality for 1d profiles""" -# -# densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_1d, -# self.n_e_1d, self.t_e_profile_1d, -# free_variable=self.psin_1d) -# -# densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_1d, -# self.n_e_profile_1d, self.t_e_1d, -# free_variable=self.psin_1d) -# -# densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], -# self.n_e_1d, self.t_e_profile_1d, -# free_variable=self.psin_1d) -# -# total = self.sumup_electrons(densities_1) -# total += self.sumup_electrons(densities_2) -# total += self.sumup_electrons(densities_3) -# -# self.assertTrue(np.allclose(total, self.n_e_profile_1d, rtol=self.TOLERANCE)) -# -# def test_balance_1d_plasma_neutrality_tcx(self): -# """test matching of plasma neutrality for 1d profiles with thermal cx""" -# -# densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_1d, -# self.n_e_1d, self.t_e_profile_1d, -# free_variable=self.psin_1d, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, tcx_donor_charge=0 -# ) -# -# densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_1d, -# self.n_e_profile_1d, self.t_e_1d, -# free_variable=self.psin_1d, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, tcx_donor_charge=0) -# -# densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], -# self.n_e_1d, self.t_e_profile_1d, -# free_variable=self.psin_1d, -# tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0) -# -# total = self.sumup_electrons(densities_1) -# total += self.sumup_electrons(densities_2) -# total += self.sumup_electrons(densities_3) -# -# self.assertTrue(np.allclose(total, self.n_e_profile_1d, rtol=self.TOLERANCE)) -# -# def test_balance_1d_interpolators_from_element_density(self): -# """ -# test calculation of 1d interpolators of charge stage densities -# """ -# -# interpolators_abundance = interpolators1d_from_elementdensity(self.atomic_data, self.element, self.psin_1d, -# self.n_element_1d, -# self.n_e_profile_1d, self.t_e_1d) -# -# profiles = self.evaluate_interpolators(interpolators_abundance, self.psin_1d) -# fraction_sum = self.sumup_fractions(profiles) -# -# self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_1d, rtol=self.TOLERANCE)) -# -# def test_balance_1d_interpolators_from_element_density_tcx(self): -# """ -# test calculation of 1d interpolators of ion charge state densities -# """ -# -# interpolators_abundance = interpolators1d_from_elementdensity(self.atomic_data, self.element, self.psin_1d, -# self.n_element_1d, -# self.n_e_profile_1d, self.t_e_1d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0) -# -# profiles = self.evaluate_interpolators(interpolators_abundance, self.psin_1d) -# fraction_sum = self.sumup_fractions(profiles) -# -# self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_1d, rtol=self.TOLERANCE)) -# -# def test_balance_1d_interpolators_plasma_neutrality(self): -# """ -# test calulation of 1d interpolators for ion charge state densities using plasma neutrality condition. -# """ -# -# interpolators_abundance_1 = interpolators1d_from_elementdensity(self.atomic_data, self.element, self.psin_1d, -# self.n_element_1d, -# self.n_e_profile_1d, self.t_e_1d) -# -# interpolators_abundance_2 = interpolators1d_from_elementdensity(self.atomic_data, self.element2, self.psin_1d, -# self.n_element2_1d, -# self.n_e_profile_1d, self.t_e_1d) -# -# interpolators_abundance_3 = interpolators1d_match_plasma_neutrality(self.atomic_data, self.element, -# self.psin_1d, -# [interpolators_abundance_1, -# interpolators_abundance_2], -# self.n_e_profile_1d, self.t_e_1d) -# -# profiles1 = self.evaluate_interpolators(interpolators_abundance_1, self.psin_1d) -# profiles2 = self.evaluate_interpolators(interpolators_abundance_2, self.psin_1d) -# profiles3 = self.evaluate_interpolators(interpolators_abundance_3, self.psin_1d) -# -# total = self.sumup_electrons(profiles1) -# total += self.sumup_electrons(profiles2) -# total += self.sumup_electrons(profiles3) -# -# self.assertTrue(np.allclose(total, self.n_e_profile_1d, rtol=self.TOLERANCE)) -# -# def test_balance_1d_interpolators_plasma_neutrality_tcx(self): -# """ -# test calulation of 1d interpolators for ion charge state densities using plasma neutrality condition. -# """ -# -# interpolators_abundance_1 = interpolators1d_from_elementdensity(self.atomic_data, self.element, self.psin_1d, -# self.n_element_1d, -# self.n_e_profile_1d, self.t_e_1d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0) -# -# interpolators_abundance_2 = interpolators1d_from_elementdensity(self.atomic_data, self.element2, self.psin_1d, -# self.n_element2_1d, -# self.n_e_profile_1d, self.t_e_1d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0) -# -# interpolators_abundance_3 = interpolators1d_match_plasma_neutrality(self.atomic_data, self.element, -# self.psin_1d, -# [interpolators_abundance_1, -# interpolators_abundance_2], -# self.n_e_profile_1d, self.t_e_1d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0) -# -# profiles1 = self.evaluate_interpolators(interpolators_abundance_1, self.psin_1d) -# profiles2 = self.evaluate_interpolators(interpolators_abundance_2, self.psin_1d) -# profiles3 = self.evaluate_interpolators(interpolators_abundance_3, self.psin_1d) -# -# total = self.sumup_electrons(profiles1) -# total += self.sumup_electrons(profiles2) -# total += self.sumup_electrons(profiles3) -# -# self.assertTrue(np.allclose(total, self.n_e_profile_1d, rtol=self.TOLERANCE)) -# -# def test_fractional_2d_from_2d(self): -# """ -# test fractional abundance 2d profile calculation with arrays as inputs -# """ -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_2d, -# self.t_e_profile_2d) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) -# -# def test_fractional_2d_from_2d_tcx(self): -# """ -# test fractional abundance 2d profile calculation with arrays as inputs with thermal cx -# """ -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_2d, -# self.t_e_profile_2d, tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_profile_2d, tcx_donor_charge=0) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) -# -# def test_fractional_2d_from_interpolators(self): -# """ -# test fractional abundance 2d profile calculation with interpolators as inputs -# """ -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_2d, -# self.t_e_2d, free_variable=(self.r, self.z)) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) -# -# def test_fractional_2d_from_interpolators_tcx(self): -# """ -# test fractional abundance 2d profile calculation with interpolators as inputs with thermal cx -# """ -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_2d, -# self.t_e_2d, tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, -# free_variable=(self.r, self.z)) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) -# -# def test_fractional_2d_from_mixed(self): -# """ -# test fractional abundance 2d profile calculation with mixed inputs -# """ -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_2d, -# self.t_e_2d, free_variable=(self.r, self.z)) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) -# -# def test_fractional_2d_from_mixed_tcx(self): -# """ -# test fractional abundance 2d profile calculation with mixed inputs with thermal cx -# """ -# abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_2d, -# self.t_e_2d, tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, -# free_variable=(self.r, self.z)) -# -# fraction_sum = self.sumup_fractions(abundance_fractional) -# self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) -# -# def test_balance_2d_elementdensity(self): -# """ -# test abundance 2d profile calculation from element density -# """ -# abundance = from_elementdensity(self.atomic_data, self.element, self.n_element_profile_2d, self.n_e_2d, -# self.t_e_profile_2d, free_variable=(self.r, self.z)) -# -# fraction_sum = self.sumup_fractions(abundance) -# self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_2d, rtol=self.TOLERANCE)) -# -# def test_balance_2d_elementdensity_tcx(self): -# """ -# test abundance 2d profile calculation from element density -# """ -# abundance = from_elementdensity(self.atomic_data, self.element, self.n_element_profile_2d, self.n_e_2d, -# self.t_e_profile_2d, tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, -# free_variable=(self.r, self.z)) -# -# fraction_sum = self.sumup_fractions(abundance) -# self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_2d, rtol=self.TOLERANCE)) -# -# def test_balance_2d_plasma_neutrality(self): -# """test matching of plasma neutrality for 2d profiles""" -# -# densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_2d, -# self.n_e_2d, self.t_e_profile_2d, -# free_variable=(self.r, self.z)) -# -# densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_2d, -# self.n_e_profile_2d, self.t_e_2d, -# free_variable=(self.r, self.z)) -# -# densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], -# self.n_e_2d, self.t_e_profile_2d, -# free_variable=(self.r, self.z)) -# -# total = self.sumup_electrons(densities_1) -# total += self.sumup_electrons(densities_2) -# total += self.sumup_electrons(densities_3) -# -# self.assertTrue(np.allclose(total, self.n_e_profile_2d, rtol=self.TOLERANCE)) -# -# def test_balance_2d_plasma_neutrality_tcx(self): -# """test matching of plasma neutrality for 2d profiles""" -# -# densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_2d, -# self.n_e_2d, self.t_e_profile_2d, tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, -# free_variable=(self.r, self.z)) -# -# densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_2d, -# self.n_e_profile_2d, self.t_e_2d, tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, -# free_variable=(self.r, self.z)) -# -# densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], -# self.n_e_2d, self.t_e_profile_2d, tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, -# free_variable=(self.r, self.z)) -# -# total = self.sumup_electrons(densities_1) -# total += self.sumup_electrons(densities_2) -# total += self.sumup_electrons(densities_3) -# -# self.assertTrue(np.allclose(total, self.n_e_profile_2d, rtol=self.TOLERANCE)) -# -# def test_fractional_inetrpolators_2d(self): -# """ -# test calculation of 1d fractional interpolators -# """ -# -# interpolators_fractional = interpolators2d_fractional(self.atomic_data, self.element, (self.r, self.z), -# self.n_e_profile_2d, self.t_e_2d) -# -# profiles = self.evaluate_interpolators(interpolators_fractional, (self.r, self.z)) -# fraction_sum = self.sumup_fractions(profiles) -# -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# def test_fractional_inetrpolators_2d_tcx(self): -# """ -# test calculation of 1d fractional interpolators with thermal cx -# """ -# -# interpolators_fractional = interpolators2d_fractional(self.atomic_data, self.element, (self.r, self.z), -# self.n_e_profile_2d, self.t_e_2d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_profile_2d, -# tcx_donor_charge=0) -# -# profiles = self.evaluate_interpolators(interpolators_fractional, (self.r, self.z)) -# fraction_sum = self.sumup_fractions(profiles) -# -# self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) -# -# def test_balance_2d_interpolators_from_element_density(self): -# """ -# test calculation of 2d interpolators of charge stage densities -# """ -# -# interpolators_abundance = interpolators2d_from_elementdensity(self.atomic_data, self.element, (self.r, self.z), -# self.n_element_2d, -# self.n_e_profile_2d, self.t_e_2d) -# -# profiles = self.evaluate_interpolators(interpolators_abundance, (self.r, self.z)) -# fraction_sum = self.sumup_fractions(profiles) -# -# self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_2d, rtol=self.TOLERANCE)) -# -# def test_balance_2d_interpolators_from_element_density_tcx(self): -# """ -# test calculation of 2d interpolators of ion charge state densities -# """ -# -# interpolators_abundance = interpolators2d_from_elementdensity(self.atomic_data, self.element, (self.r, self.z), -# self.n_element_2d, -# self.n_e_profile_2d, self.t_e_2d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_2d, -# tcx_donor_charge=0) -# -# profiles = self.evaluate_interpolators(interpolators_abundance, (self.r, self.z)) -# fraction_sum = self.sumup_fractions(profiles) -# -# self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_2d, rtol=self.TOLERANCE)) -# -# def test_balance_2d_interpolators_plasma_neutrality(self): -# """ -# test calulation of 2d interpolators for ion charge state densities using plasma neutrality condition. -# """ -# -# interpolators_abundance_1 = interpolators2d_from_elementdensity(self.atomic_data, self.element, -# (self.r, self.z), -# self.n_element_2d, -# self.n_e_profile_2d, self.t_e_2d) -# -# interpolators_abundance_2 = interpolators2d_from_elementdensity(self.atomic_data, self.element2, -# (self.r, self.z), -# self.n_element2_2d, -# self.n_e_profile_2d, self.t_e_2d) -# -# interpolators_abundance_3 = interpolators2d_match_plasma_neutrality(self.atomic_data, self.element, -# (self.r, self.z), -# [interpolators_abundance_1, -# interpolators_abundance_2], -# self.n_e_profile_2d, self.t_e_2d) -# -# profiles1 = self.evaluate_interpolators(interpolators_abundance_1, (self.r, self.z)) -# profiles2 = self.evaluate_interpolators(interpolators_abundance_2, (self.r, self.z)) -# profiles3 = self.evaluate_interpolators(interpolators_abundance_3, (self.r, self.z)) -# -# total = self.sumup_electrons(profiles1) -# total += self.sumup_electrons(profiles2) -# total += self.sumup_electrons(profiles3) -# -# self.assertTrue(np.allclose(total, self.n_e_profile_2d, rtol=self.TOLERANCE)) -# -# def test_balance_2d_interpolators_plasma_neutrality_tcx(self): -# """ -# test calulation of 2d interpolators for ion charge state densities using plasma neutrality condition. -# """ -# -# interpolators_abundance_1 = interpolators2d_from_elementdensity(self.atomic_data, self.element, -# (self.r, self.z), -# self.n_element_2d, -# self.n_e_profile_2d, self.t_e_2d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_2d, -# tcx_donor_charge=0) -# -# interpolators_abundance_2 = interpolators2d_from_elementdensity(self.atomic_data, self.element2, -# (self.r, self.z), -# self.n_element2_2d, -# self.n_e_profile_2d, self.t_e_2d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_2d, -# tcx_donor_charge=0) -# -# interpolators_abundance_3 = interpolators2d_match_plasma_neutrality(self.atomic_data, self.element, -# (self.r, self.z), -# [interpolators_abundance_1, -# interpolators_abundance_2], -# self.n_e_profile_2d, self.t_e_2d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_2d, -# tcx_donor_charge=0) -# -# profiles1 = self.evaluate_interpolators(interpolators_abundance_1, (self.r, self.z)) -# profiles2 = self.evaluate_interpolators(interpolators_abundance_2, (self.r, self.z)) -# profiles3 = self.evaluate_interpolators(interpolators_abundance_3, (self.r, self.z)) -# -# total = self.sumup_electrons(profiles1) -# total += self.sumup_electrons(profiles2) -# total += self.sumup_electrons(profiles3) -# -# self.assertTrue(np.allclose(total, self.n_e_profile_2d, rtol=self.TOLERANCE)) -# -# def test_axisymmetric_mapper(self): -# "test axisymmetric mapping" -# -# #fractional abundance -# interpolators_fractional = interpolators2d_fractional(self.atomic_data, self.element, (self.r, self.z), -# self.n_e_profile_2d, self.t_e_2d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_profile_2d, -# tcx_donor_charge=0) -# -# mappers = abundance_axisymmetric_mapper(interpolators_fractional) -# -# -# -# profile = self.evaluate_interpolators(mappers, self.r) -# -# total = self.sumup_fractions(profile) -# -# -# self.assertTrue(np.allclose(total, 1, rtol=self.TOLERANCE)) -# -# def test_equilibrium_map3d_fractional(self): -# """ -# test calculation of fractional abundance and application of map3d functionality of equilibrium -# """ -# -# equilibrium = example_equilibrium() -# -# mapper = equilibrium_map3d_fractional(self.atomic_data, self.element, equilibrium, self.psin_1d, -# self.n_e_profile_1d, self.t_e_profile_1d, self.tcx_donor, -# self.n_tcx_donor_profile_1d, tcx_donor_charge=0) -# psin_1d = np.linspace(0, 0.99, 10) -# r = np.zeros_like(psin_1d) -# for index, value in enumerate(psin_1d): -# r[index] = equilibrium.psin_to_r(value) -# -# profile = self.evaluate_interpolators(mapper, r) -# -# total = self.sumup_fractions(profile) -# -# self.assertTrue(np.allclose(total, 1, rtol=self.TOLERANCE)) -# -# def test_equilibrium_map3d_from_elementdensity(self): -# """ -# test calculation of abundance and application of map3d functionality of equilibrium -# """ -# -# equilibrium = example_equilibrium() -# -# mapper = equilibrium_map3d_from_elementdensity(self.atomic_data, self.element, equilibrium, self.psin_1d, -# self.n_element_1d, self.n_e_1d, -# self.t_e_1d, self.tcx_donor, self.n_tcx_donor_1d, -# tcx_donor_charge=0) -# -# inlcfs = np.where(self.psin_1d < 1)[0] -# psin_1d = self.psin_1d[inlcfs] -# n = np.zeros_like(psin_1d) -# r = np.zeros_like(psin_1d) -# for index, value in enumerate(psin_1d): -# n[index] = self.n_element_1d(value) -# r[index] = equilibrium.psin_to_r(value) -# -# profile = self.evaluate_interpolators(mapper, r) -# total = self.sumup_fractions(profile) -# -# self.assertTrue(np.allclose(total, n, rtol=self.TOLERANCE)) -# -# def test_equilibrium_map3d_plasma_neutrality(self): -# """ -# test calculation of abundance and application of map3d functionality of equilibrium -# """ -# -# equilibrium = example_equilibrium() -# -# interpolators_element1 = interpolators1d_from_elementdensity(self.atomic_data, self.element, self.psin_1d, -# self.n_element_1d, -# self.n_e_1d, self.t_e_1d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0) -# -# interpolators_element2 = interpolators1d_from_elementdensity(self.atomic_data, self.element2, self.psin_1d, -# self.n_element2_1d, -# self.n_e_1d, self.t_e_1d, -# tcx_donor=self.tcx_donor, -# tcx_donor_n=self.n_tcx_donor_1d, -# tcx_donor_charge=0) -# -# mapper = equilibrium_map3d_match_plasma_neutrality(self.atomic_data, self.element_bulk, equilibrium, self.psin_1d, -# [interpolators_element1, interpolators_element2], -# self.n_e_1d, self.t_e_1d, self.tcx_donor, -# self.n_tcx_donor_1d, 0) -# -# inlcfs = np.where(self.psin_1d < 1)[0] -# psin_1d = self.psin_1d[inlcfs] -# n = np.zeros_like(psin_1d) -# r = np.zeros_like(psin_1d) -# for index, value in enumerate(psin_1d): -# n[index] = self.n_e_1d(value) -# r[index] = equilibrium.psin_to_r(value) -# -# profile1 = self.evaluate_interpolators(interpolators_element1, psin_1d) -# profile2 = self.evaluate_interpolators(interpolators_element2, psin_1d) -# profile3 = self.evaluate_interpolators(mapper, r) -# -# total = self.sumup_electrons(profile1) -# total += self.sumup_electrons(profile2) -# total += self.sumup_electrons(profile3) -# -# self.assertTrue(np.allclose(total, n, rtol=self.TOLERANCE)) + +# Copyright 2016-2021 Euratom +# Copyright 2016-2021 United Kingdom Atomic Energy Authority +# Copyright 2016-2021 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + + +import unittest +import os +from collections.abc import Iterable +import numpy as np +from raysect.core.math.function.float import Function1D, Function2D, Interpolator1DArray, Interpolator2DArray + +from cherab.core.atomic import neon, hydrogen, helium +from cherab.core.math import AxisymmetricMapper +from cherab.openadas import OpenADAS +from cherab.tools.plasmas.ionisation_balance import (fractional_abundance, from_elementdensity, match_plasma_neutrality, + interpolators1d_fractional, interpolators1d_from_elementdensity, + interpolators1d_match_plasma_neutrality, + interpolators2d_fractional, interpolators2d_from_elementdensity, + interpolators2d_match_plasma_neutrality, + abundance_axisymmetric_mapper, + equilibrium_map3d_fractional, equilibrium_map3d_from_elementdensity, + equilibrium_map3d_match_plasma_neutrality) + +from cherab.tools.equilibrium import example_equilibrium + + +def double_parabola(r, centre, edge, p, q): + return (centre - edge) * np.power((1 - np.power((r - r.min()) / (r.max() - r.min()), p)), q) + edge + + +def normal(x, mu, sd, height=1, offset=0): + return height * np.exp(-1 * np.power(x - mu, 2) / (2 * sd ** 2)) + offset + + +def get_electron_density_profile(abundances): + n_e = np.zeros((abundances[0].shape[1])) + for abundance in abundances: + for rownumber, row in enumerate(abundance.T): + n_e[rownumber] += np.sum(row * np.arange(row.shape[0])) + + return n_e + + +def get_electron_density_spot(densities): + n_e = 0 + for spec in densities: + for index, value in enumerate(spec): + n_e += index * value + + return n_e + + +def exp_decay(r, lamb, max_val): + return max_val * np.exp((r - r.max()) * lamb) + + +class TestIonizationBalance1D(unittest.TestCase): + + # create plasma profiles and interpolators + # 1d profiles + psin_1d = np.linspace(0, 1.1, 15, endpoint=True) + psin_1d_detailed = np.linspace(0, 1.1, 50, endpoint=True) + + t_e_profile_1d = double_parabola(psin_1d, 5000, 10, 2, 2) + n_e_profile_1d = double_parabola(psin_1d, 6e19, 5e18, 2, 2) + + t_element_profile_1d = double_parabola(psin_1d, 1500, 40, 2, 2) + n_element_profile_1d = double_parabola(psin_1d, 1e17, 1e17, 2, 2) + normal(psin_1d, 0.9, 0.1, 5e17) + n_element2_profile_1d = double_parabola(psin_1d, 5e17, 1e17, 2, 2) + + n_tcx_donor_profile_1d = exp_decay(psin_1d, 10, 3e16) + + t_e_1d = Interpolator1DArray(psin_1d, t_e_profile_1d, 'cubic', 'none', 0) + n_e_1d = Interpolator1DArray(psin_1d, n_e_profile_1d, 'cubic', 'none', 0) + + t_element_1d = Interpolator1DArray(psin_1d, t_element_profile_1d, 'cubic', 'none', 0) + n_element_1d = Interpolator1DArray(psin_1d, n_element_profile_1d, 'cubic', 'none', 0) + n_element2_1d = Interpolator1DArray(psin_1d, n_element2_profile_1d, 'cubic', 'none', 0) + + n_tcx_donor_1d = Interpolator1DArray(psin_1d, n_tcx_donor_profile_1d, 'cubic', 'none', 0) + + # denser psi array to test interpolators + n_e_profile_detailed_1d = np.zeros_like(psin_1d_detailed) + for index, value in enumerate(psin_1d_detailed): + n_e_profile_detailed_1d[index] = n_e_1d(value) + + # 2d profiles + + r = np.linspace(-1, 1, 6) + z = np.linspace(-1, 1, 8) + + psin_2d = np.zeros((*r.shape, *z.shape)) + + for index0, value0 in enumerate(r): + for index1, value1 in enumerate(z): + if np.sqrt(value0 ** 2 + value1 ** 2) < psin_1d.max(): + psin_2d[index0, index1] = np.sqrt(value0 ** 2 + value1 ** 2) + else: + psin_2d[index0, index1] = psin_1d.max() + + t_e_profile_2d = np.zeros_like(psin_2d) + for index in np.ndindex(*t_e_profile_2d.shape): + t_e_profile_2d[index] = t_e_1d(psin_2d[index]) + + t_e_2d = Interpolator2DArray(r, z, t_e_profile_2d, 'cubic', 'none', 0, 0) + + n_e_profile_2d = np.zeros_like(psin_2d) + for index in np.ndindex(*n_e_profile_2d.shape): + n_e_profile_2d[index] = n_e_1d(psin_2d[index]) + + n_e_2d = Interpolator2DArray(r, z, n_e_profile_2d, 'cubic', 'none', 0, 0) + + t_element_profile_2d = np.zeros_like(psin_2d) + for index in np.ndindex(*t_element_profile_2d.shape): + t_element_profile_2d[index] = t_element_1d(psin_2d[index]) + + t_element_2d = Interpolator2DArray(r, z, t_element_profile_2d, 'cubic', 'none', 0, 0) + + n_element_profile_2d = np.zeros_like(psin_2d) + for index in np.ndindex(*n_element_profile_2d.shape): + n_element_profile_2d[index] = n_element_1d(psin_2d[index]) + + n_element_2d = Interpolator2DArray(r, z, n_element_profile_2d, 'cubic', 'none', 0, 0) + + n_element2_profile_2d = np.zeros_like(psin_2d) + for index in np.ndindex(*n_element2_profile_2d.shape): + n_element2_profile_2d[index] = n_element2_1d(psin_2d[index]) + + n_element2_2d = Interpolator2DArray(r, z, n_element2_profile_2d, 'cubic', 'none', 0, 0) + + n_tcx_donor_profile_2d = np.zeros_like(psin_2d) + for index in np.ndindex(*n_tcx_donor_profile_2d.shape): + n_tcx_donor_profile_2d[index] = n_tcx_donor_1d(psin_2d[index]) + + n_tcx_donor_2d = Interpolator2DArray(r, z, n_element_profile_2d, 'cubic', 'none', 0, 0) + + # define psi for single-point tests + psi_value = 0.9 + + # load adas atomic database and define elements + repository_path = os.path.join(os.path.dirname(__file__), 'data/atomic_rates_mockup') + atomic_data = OpenADAS(data_path=repository_path, permit_extrapolation=True) + + element = neon + element2 = helium + element_bulk = hydrogen + + tcx_donor = hydrogen + + TOLERANCE = 1e-3 + + def sumup_fractions(self, fractions): + + if isinstance(fractions, dict): + if isinstance(fractions[0], Iterable): + total = np.zeros_like(fractions[0]) + else: + total = 0 + + for index, values in fractions.items(): + total += values + + elif isinstance(fractions, np.ndarray): + total = np.zeros_like(fractions[0, ...]) + + for index in np.ndindex(fractions.shape): + total[index] += fractions[index] + + return total + + def sumup_electrons(self, densities): + if isinstance(densities, dict): + if isinstance(densities[0], Iterable): + total = np.zeros_like(densities[0]) + else: + total = 0 + + for index, values in densities.items(): + total += values * index + + elif isinstance(densities, np.ndarray): + total = np.zeros_like(densities[0, ...]) + + for index in np.ndindex(densities.shape): + total[index] += densities[index] * index[0] + + return total + + def evaluate_interpolators(self, interpolators, free_variable): + + profiles = {} + for key, item in interpolators.items(): + if isinstance(item, Function1D): + profiles[key] = np.zeros_like(free_variable) + for index in np.ndindex(*free_variable.shape): + profiles[key][index] = item(free_variable[index]) + elif isinstance(item, Function2D): + profiles[key] = np.zeros((*free_variable[0].shape, *free_variable[1].shape)) + for index0, value0 in enumerate(free_variable[0]): + for index1, value1 in enumerate(free_variable[1]): + profiles[key][index0, index1] = item(value0, value1) + elif isinstance(item, AxisymmetricMapper): + profiles[key] = np.zeros_like(free_variable) + for index, value in enumerate(free_variable): + profiles[key][index] = item(value, 0, 0) + + return profiles + + def test_fractional_0d_from_0d(self): + """ + test fractional abundance calculation with float numbers as inputs + """ + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d(self.psi_value), + self.t_e_1d(self.psi_value)) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(1 - fraction_sum < self.TOLERANCE) + + def test_fractional_0d_from_0d_tcx(self): + """ + test fractional abundance calculation with thermal cx and float numbers as inputs + """ + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d(self.psi_value), + self.t_e_1d(self.psi_value), + self.tcx_donor, self.n_tcx_donor_1d(self.psi_value), 0) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(1 - fraction_sum < self.TOLERANCE) + + def test_fractional_0d_from_interpolators(self): + """ + test interpolators and free_variable as inputs + """ + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_1d, + free_variable=self.psi_value) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(1 - fraction_sum < self.TOLERANCE) + + def test_fractional_0d_from_interpolators_tcx(self): + """ + test interpolators and free_variable as inputs + """ + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_1d, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0, free_variable=self.psi_value) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(1 - fraction_sum < self.TOLERANCE) + + def test_fractional_0d_from_mixed(self): + """ + test mixed types of inputs + """ + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d(self.psi_value), + self.t_e_1d, + free_variable=self.psi_value) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(1 - fraction_sum < self.TOLERANCE) + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, + self.t_e_1d(self.psi_value), + free_variable=self.psi_value) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(1 - fraction_sum < self.TOLERANCE) + + def test_fractional_1d_from_1d(self): + """ + test calculation of 1d fractional profiles with 1d iterables as inputs + """ + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_1d, + self.t_e_profile_1d) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + def test_fractional_1d_from_1d_tcx(self): + """ + test calculation of 1d fractional profiles with 1d iterables as inputs + """ + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_1d, + self.t_e_profile_1d, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_profile_1d, + tcx_donor_charge=0) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + def test_fractional_1d_from_interpolators(self): + """ + test calculation of 1d fractional profiles with 1d interpolators as inputs + """ + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_1d, + free_variable=self.psin_1d) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + def test_fractional_1d_from_interpolators_tcx(self): + """ + test calculation of 1d fractional profiles with 1d iterables as inputs + """ + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_1d, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0, free_variable=self.psin_1d) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + def test_fractional_from_mixed(self): + """ + test calculation of 1d fractional profiles with mixed types as inputs + """ + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_1d, self.t_e_1d, + free_variable=self.psin_1d) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_profile_1d, + free_variable=self.psin_1d) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + def test_fractional_from_mixed_tcx(self): + """ + test calculation of 1d fractional profiles with mixed types as inputs + """ + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_1d, self.t_e_1d, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0, free_variable=self.psin_1d) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_profile_1d, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0, free_variable=self.psin_1d) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_1d, self.t_e_profile_1d, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_profile_1d, + tcx_donor_charge=0, free_variable=self.psin_1d) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + def test_fractional_inetrpolators_1d(self): + """ + test calculation of 1d fractional interpolators + """ + + interpolators_fractional = interpolators1d_fractional(self.atomic_data, self.element, self.psin_1d, + self.n_e_profile_1d, self.t_e_1d) + + profiles = self.evaluate_interpolators(interpolators_fractional, self.psin_1d) + fraction_sum = self.sumup_fractions(profiles) + + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + def test_fractional_inetrpolators_1d_tcx(self): + """ + test calculation of 1d fractional interpolators with thermal cx + """ + + interpolators_fractional = interpolators1d_fractional(self.atomic_data, self.element, self.psin_1d, + self.n_e_profile_1d, self.t_e_1d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_profile_1d, + tcx_donor_charge=0) + + profiles = self.evaluate_interpolators(interpolators_fractional, self.psin_1d) + fraction_sum = self.sumup_fractions(profiles) + + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + def test_balance_0d_elementdensity(self): + """ + test calculation of ionization balance + """ + # test with floats as input + densities = from_elementdensity(self.atomic_data, self.element, self.n_element_1d(self.psi_value), + self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value)) + total = self.sumup_fractions(densities) + self.assertTrue(np.isclose(total, self.n_element_1d(self.psi_value), rtol=self.TOLERANCE)) + + # test with interpolators + densities = from_elementdensity(self.atomic_data, self.element, self.n_element_1d(self.psi_value), + self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value)) + total = self.sumup_fractions(densities) + self.assertTrue(np.isclose(total, self.n_element_1d(self.psi_value), rtol=self.TOLERANCE)) + + # test with mixed parameters + densities = from_elementdensity(self.atomic_data, self.element, self.n_element_1d, + self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value), + free_variable=self.psi_value) + total = self.sumup_fractions(densities) + self.assertTrue(np.isclose(total, self.n_element_1d(self.psi_value), rtol=self.TOLERANCE)) + + def test_balance_0d_plasma_neutrality(self): + """test matching of plasma neutrality""" + + densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_1d, + self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value), + free_variable=self.psi_value) + + densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_1d, + self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value), + free_variable=self.psi_value) + + densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], + self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value)) + + total = self.sumup_electrons(densities_1) + total += self.sumup_electrons(densities_2) + total += self.sumup_electrons(densities_3) + + self.assertTrue(np.isclose(total, self.n_e_1d(self.psi_value), rtol=self.TOLERANCE)) + + def test_balance_0d_plasma_neutrality_tcx(self): + """test matching of plasma neutrality""" + + densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_1d, + self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value), + free_variable=self.psi_value, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, tcx_donor_charge=0, + ) + + densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_1d, + self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value), + free_variable=self.psi_value, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, tcx_donor_charge=0, + ) + + densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], + self.n_e_1d(self.psi_value), self.t_e_1d(self.psi_value)) + + total = self.sumup_electrons(densities_1) + total += self.sumup_electrons(densities_2) + total += self.sumup_electrons(densities_3) + + self.assertTrue(np.isclose(total, self.n_e_1d(self.psi_value), rtol=self.TOLERANCE)) + + def test_balance_1d_plasma_neutrality(self): + """test matching of plasma neutrality for 1d profiles""" + + densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_1d, + self.n_e_1d, self.t_e_profile_1d, + free_variable=self.psin_1d) + + densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_1d, + self.n_e_profile_1d, self.t_e_1d, + free_variable=self.psin_1d) + + densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], + self.n_e_1d, self.t_e_profile_1d, + free_variable=self.psin_1d) + + total = self.sumup_electrons(densities_1) + total += self.sumup_electrons(densities_2) + total += self.sumup_electrons(densities_3) + + self.assertTrue(np.allclose(total, self.n_e_profile_1d, rtol=self.TOLERANCE)) + + def test_balance_1d_plasma_neutrality_tcx(self): + """test matching of plasma neutrality for 1d profiles with thermal cx""" + + densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_1d, + self.n_e_1d, self.t_e_profile_1d, + free_variable=self.psin_1d, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, tcx_donor_charge=0 + ) + + densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_1d, + self.n_e_profile_1d, self.t_e_1d, + free_variable=self.psin_1d, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, tcx_donor_charge=0) + + densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], + self.n_e_1d, self.t_e_profile_1d, + free_variable=self.psin_1d, + tcx_donor=self.tcx_donor, tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0) + + total = self.sumup_electrons(densities_1) + total += self.sumup_electrons(densities_2) + total += self.sumup_electrons(densities_3) + + self.assertTrue(np.allclose(total, self.n_e_profile_1d, rtol=self.TOLERANCE)) + + def test_balance_1d_interpolators_from_element_density(self): + """ + test calculation of 1d interpolators of charge stage densities + """ + + interpolators_abundance = interpolators1d_from_elementdensity(self.atomic_data, self.element, self.psin_1d, + self.n_element_1d, + self.n_e_profile_1d, self.t_e_1d) + + profiles = self.evaluate_interpolators(interpolators_abundance, self.psin_1d) + fraction_sum = self.sumup_fractions(profiles) + + self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_1d, rtol=self.TOLERANCE)) + + def test_balance_1d_interpolators_from_element_density_tcx(self): + """ + test calculation of 1d interpolators of ion charge state densities + """ + + interpolators_abundance = interpolators1d_from_elementdensity(self.atomic_data, self.element, self.psin_1d, + self.n_element_1d, + self.n_e_profile_1d, self.t_e_1d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0) + + profiles = self.evaluate_interpolators(interpolators_abundance, self.psin_1d) + fraction_sum = self.sumup_fractions(profiles) + + self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_1d, rtol=self.TOLERANCE)) + + def test_balance_1d_interpolators_plasma_neutrality(self): + """ + test calulation of 1d interpolators for ion charge state densities using plasma neutrality condition. + """ + + interpolators_abundance_1 = interpolators1d_from_elementdensity(self.atomic_data, self.element, self.psin_1d, + self.n_element_1d, + self.n_e_profile_1d, self.t_e_1d) + + interpolators_abundance_2 = interpolators1d_from_elementdensity(self.atomic_data, self.element2, self.psin_1d, + self.n_element2_1d, + self.n_e_profile_1d, self.t_e_1d) + + interpolators_abundance_3 = interpolators1d_match_plasma_neutrality(self.atomic_data, self.element, + self.psin_1d, + [interpolators_abundance_1, + interpolators_abundance_2], + self.n_e_profile_1d, self.t_e_1d) + + profiles1 = self.evaluate_interpolators(interpolators_abundance_1, self.psin_1d) + profiles2 = self.evaluate_interpolators(interpolators_abundance_2, self.psin_1d) + profiles3 = self.evaluate_interpolators(interpolators_abundance_3, self.psin_1d) + + total = self.sumup_electrons(profiles1) + total += self.sumup_electrons(profiles2) + total += self.sumup_electrons(profiles3) + + self.assertTrue(np.allclose(total, self.n_e_profile_1d, rtol=self.TOLERANCE)) + + def test_balance_1d_interpolators_plasma_neutrality_tcx(self): + """ + test calulation of 1d interpolators for ion charge state densities using plasma neutrality condition. + """ + + interpolators_abundance_1 = interpolators1d_from_elementdensity(self.atomic_data, self.element, self.psin_1d, + self.n_element_1d, + self.n_e_profile_1d, self.t_e_1d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0) + + interpolators_abundance_2 = interpolators1d_from_elementdensity(self.atomic_data, self.element2, self.psin_1d, + self.n_element2_1d, + self.n_e_profile_1d, self.t_e_1d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0) + + interpolators_abundance_3 = interpolators1d_match_plasma_neutrality(self.atomic_data, self.element, + self.psin_1d, + [interpolators_abundance_1, + interpolators_abundance_2], + self.n_e_profile_1d, self.t_e_1d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0) + + profiles1 = self.evaluate_interpolators(interpolators_abundance_1, self.psin_1d) + profiles2 = self.evaluate_interpolators(interpolators_abundance_2, self.psin_1d) + profiles3 = self.evaluate_interpolators(interpolators_abundance_3, self.psin_1d) + + total = self.sumup_electrons(profiles1) + total += self.sumup_electrons(profiles2) + total += self.sumup_electrons(profiles3) + + self.assertTrue(np.allclose(total, self.n_e_profile_1d, rtol=self.TOLERANCE)) + + def test_fractional_2d_from_2d(self): + """ + test fractional abundance 2d profile calculation with arrays as inputs + """ + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_2d, + self.t_e_profile_2d) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) + + def test_fractional_2d_from_2d_tcx(self): + """ + test fractional abundance 2d profile calculation with arrays as inputs with thermal cx + """ + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_2d, + self.t_e_profile_2d, tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_profile_2d, tcx_donor_charge=0) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) + + def test_fractional_2d_from_interpolators(self): + """ + test fractional abundance 2d profile calculation with interpolators as inputs + """ + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_2d, + self.t_e_2d, free_variable=(self.r, self.z)) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) + + def test_fractional_2d_from_interpolators_tcx(self): + """ + test fractional abundance 2d profile calculation with interpolators as inputs with thermal cx + """ + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_2d, + self.t_e_2d, tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, + free_variable=(self.r, self.z)) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) + + def test_fractional_2d_from_mixed(self): + """ + test fractional abundance 2d profile calculation with mixed inputs + """ + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_2d, + self.t_e_2d, free_variable=(self.r, self.z)) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) + + def test_fractional_2d_from_mixed_tcx(self): + """ + test fractional abundance 2d profile calculation with mixed inputs with thermal cx + """ + abundance_fractional = fractional_abundance(self.atomic_data, self.element, self.n_e_profile_2d, + self.t_e_2d, tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, + free_variable=(self.r, self.z)) + + fraction_sum = self.sumup_fractions(abundance_fractional) + self.assertTrue(np.allclose(fraction_sum, np.ones_like(self.n_e_profile_2d), rtol=self.TOLERANCE)) + + def test_balance_2d_elementdensity(self): + """ + test abundance 2d profile calculation from element density + """ + abundance = from_elementdensity(self.atomic_data, self.element, self.n_element_profile_2d, self.n_e_2d, + self.t_e_profile_2d, free_variable=(self.r, self.z)) + + fraction_sum = self.sumup_fractions(abundance) + self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_2d, rtol=self.TOLERANCE)) + + def test_balance_2d_elementdensity_tcx(self): + """ + test abundance 2d profile calculation from element density + """ + abundance = from_elementdensity(self.atomic_data, self.element, self.n_element_profile_2d, self.n_e_2d, + self.t_e_profile_2d, tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, + free_variable=(self.r, self.z)) + + fraction_sum = self.sumup_fractions(abundance) + self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_2d, rtol=self.TOLERANCE)) + + def test_balance_2d_plasma_neutrality(self): + """test matching of plasma neutrality for 2d profiles""" + + densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_2d, + self.n_e_2d, self.t_e_profile_2d, + free_variable=(self.r, self.z)) + + densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_2d, + self.n_e_profile_2d, self.t_e_2d, + free_variable=(self.r, self.z)) + + densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], + self.n_e_2d, self.t_e_profile_2d, + free_variable=(self.r, self.z)) + + total = self.sumup_electrons(densities_1) + total += self.sumup_electrons(densities_2) + total += self.sumup_electrons(densities_3) + + self.assertTrue(np.allclose(total, self.n_e_profile_2d, rtol=self.TOLERANCE)) + + def test_balance_2d_plasma_neutrality_tcx(self): + """test matching of plasma neutrality for 2d profiles""" + + densities_1 = from_elementdensity(self.atomic_data, self.element, self.n_element_2d, + self.n_e_2d, self.t_e_profile_2d, tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, + free_variable=(self.r, self.z)) + + densities_2 = from_elementdensity(self.atomic_data, self.element2, self.n_element2_2d, + self.n_e_profile_2d, self.t_e_2d, tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, + free_variable=(self.r, self.z)) + + densities_3 = match_plasma_neutrality(self.atomic_data, self.element_bulk, [densities_1, densities_2], + self.n_e_2d, self.t_e_profile_2d, tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_2d, tcx_donor_charge=0, + free_variable=(self.r, self.z)) + + total = self.sumup_electrons(densities_1) + total += self.sumup_electrons(densities_2) + total += self.sumup_electrons(densities_3) + + self.assertTrue(np.allclose(total, self.n_e_profile_2d, rtol=self.TOLERANCE)) + + def test_fractional_inetrpolators_2d(self): + """ + test calculation of 1d fractional interpolators + """ + + interpolators_fractional = interpolators2d_fractional(self.atomic_data, self.element, (self.r, self.z), + self.n_e_profile_2d, self.t_e_2d) + + profiles = self.evaluate_interpolators(interpolators_fractional, (self.r, self.z)) + fraction_sum = self.sumup_fractions(profiles) + + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + def test_fractional_inetrpolators_2d_tcx(self): + """ + test calculation of 1d fractional interpolators with thermal cx + """ + + interpolators_fractional = interpolators2d_fractional(self.atomic_data, self.element, (self.r, self.z), + self.n_e_profile_2d, self.t_e_2d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_profile_2d, + tcx_donor_charge=0) + + profiles = self.evaluate_interpolators(interpolators_fractional, (self.r, self.z)) + fraction_sum = self.sumup_fractions(profiles) + + self.assertTrue(np.allclose(fraction_sum, 1, atol=self.TOLERANCE)) + + def test_balance_2d_interpolators_from_element_density(self): + """ + test calculation of 2d interpolators of charge stage densities + """ + + interpolators_abundance = interpolators2d_from_elementdensity(self.atomic_data, self.element, (self.r, self.z), + self.n_element_2d, + self.n_e_profile_2d, self.t_e_2d) + + profiles = self.evaluate_interpolators(interpolators_abundance, (self.r, self.z)) + fraction_sum = self.sumup_fractions(profiles) + + self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_2d, rtol=self.TOLERANCE)) + + def test_balance_2d_interpolators_from_element_density_tcx(self): + """ + test calculation of 2d interpolators of ion charge state densities + """ + + interpolators_abundance = interpolators2d_from_elementdensity(self.atomic_data, self.element, (self.r, self.z), + self.n_element_2d, + self.n_e_profile_2d, self.t_e_2d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_2d, + tcx_donor_charge=0) + + profiles = self.evaluate_interpolators(interpolators_abundance, (self.r, self.z)) + fraction_sum = self.sumup_fractions(profiles) + + self.assertTrue(np.allclose(fraction_sum, self.n_element_profile_2d, rtol=self.TOLERANCE)) + + def test_balance_2d_interpolators_plasma_neutrality(self): + """ + test calulation of 2d interpolators for ion charge state densities using plasma neutrality condition. + """ + + interpolators_abundance_1 = interpolators2d_from_elementdensity(self.atomic_data, self.element, + (self.r, self.z), + self.n_element_2d, + self.n_e_profile_2d, self.t_e_2d) + + interpolators_abundance_2 = interpolators2d_from_elementdensity(self.atomic_data, self.element2, + (self.r, self.z), + self.n_element2_2d, + self.n_e_profile_2d, self.t_e_2d) + + interpolators_abundance_3 = interpolators2d_match_plasma_neutrality(self.atomic_data, self.element, + (self.r, self.z), + [interpolators_abundance_1, + interpolators_abundance_2], + self.n_e_profile_2d, self.t_e_2d) + + profiles1 = self.evaluate_interpolators(interpolators_abundance_1, (self.r, self.z)) + profiles2 = self.evaluate_interpolators(interpolators_abundance_2, (self.r, self.z)) + profiles3 = self.evaluate_interpolators(interpolators_abundance_3, (self.r, self.z)) + + total = self.sumup_electrons(profiles1) + total += self.sumup_electrons(profiles2) + total += self.sumup_electrons(profiles3) + + self.assertTrue(np.allclose(total, self.n_e_profile_2d, rtol=self.TOLERANCE)) + + def test_balance_2d_interpolators_plasma_neutrality_tcx(self): + """ + test calulation of 2d interpolators for ion charge state densities using plasma neutrality condition. + """ + + interpolators_abundance_1 = interpolators2d_from_elementdensity(self.atomic_data, self.element, + (self.r, self.z), + self.n_element_2d, + self.n_e_profile_2d, self.t_e_2d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_2d, + tcx_donor_charge=0) + + interpolators_abundance_2 = interpolators2d_from_elementdensity(self.atomic_data, self.element2, + (self.r, self.z), + self.n_element2_2d, + self.n_e_profile_2d, self.t_e_2d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_2d, + tcx_donor_charge=0) + + interpolators_abundance_3 = interpolators2d_match_plasma_neutrality(self.atomic_data, self.element, + (self.r, self.z), + [interpolators_abundance_1, + interpolators_abundance_2], + self.n_e_profile_2d, self.t_e_2d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_2d, + tcx_donor_charge=0) + + profiles1 = self.evaluate_interpolators(interpolators_abundance_1, (self.r, self.z)) + profiles2 = self.evaluate_interpolators(interpolators_abundance_2, (self.r, self.z)) + profiles3 = self.evaluate_interpolators(interpolators_abundance_3, (self.r, self.z)) + + total = self.sumup_electrons(profiles1) + total += self.sumup_electrons(profiles2) + total += self.sumup_electrons(profiles3) + + self.assertTrue(np.allclose(total, self.n_e_profile_2d, rtol=self.TOLERANCE)) + + def test_axisymmetric_mapper(self): + "test axisymmetric mapping" + + # fractional abundance + interpolators_fractional = interpolators2d_fractional(self.atomic_data, self.element, (self.r, self.z), + self.n_e_profile_2d, self.t_e_2d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_profile_2d, + tcx_donor_charge=0) + + mappers = abundance_axisymmetric_mapper(interpolators_fractional) + + profile = self.evaluate_interpolators(mappers, self.r) + + total = self.sumup_fractions(profile) + + self.assertTrue(np.allclose(total, 1, rtol=self.TOLERANCE)) + + def test_equilibrium_map3d_fractional(self): + """ + test calculation of fractional abundance and application of map3d functionality of equilibrium + """ + + equilibrium = example_equilibrium() + + mapper = equilibrium_map3d_fractional(self.atomic_data, self.element, equilibrium, self.psin_1d, + self.n_e_profile_1d, self.t_e_profile_1d, self.tcx_donor, + self.n_tcx_donor_profile_1d, tcx_donor_charge=0) + psin_1d = np.linspace(0, 0.99, 10) + r = np.zeros_like(psin_1d) + for index, value in enumerate(psin_1d): + r[index] = equilibrium.psin_to_r(value) + + profile = self.evaluate_interpolators(mapper, r) + + total = self.sumup_fractions(profile) + + self.assertTrue(np.allclose(total, 1, rtol=self.TOLERANCE)) + + def test_equilibrium_map3d_from_elementdensity(self): + """ + test calculation of abundance and application of map3d functionality of equilibrium + """ + + equilibrium = example_equilibrium() + + mapper = equilibrium_map3d_from_elementdensity(self.atomic_data, self.element, equilibrium, self.psin_1d, + self.n_element_1d, self.n_e_1d, + self.t_e_1d, self.tcx_donor, self.n_tcx_donor_1d, + tcx_donor_charge=0) + + inlcfs = np.where(self.psin_1d < 1)[0] + psin_1d = self.psin_1d[inlcfs] + n = np.zeros_like(psin_1d) + r = np.zeros_like(psin_1d) + for index, value in enumerate(psin_1d): + n[index] = self.n_element_1d(value) + r[index] = equilibrium.psin_to_r(value) + + profile = self.evaluate_interpolators(mapper, r) + total = self.sumup_fractions(profile) + + self.assertTrue(np.allclose(total, n, rtol=self.TOLERANCE)) + + def test_equilibrium_map3d_plasma_neutrality(self): + """ + test calculation of abundance and application of map3d functionality of equilibrium + """ + + equilibrium = example_equilibrium() + + interpolators_element1 = interpolators1d_from_elementdensity(self.atomic_data, self.element, self.psin_1d, + self.n_element_1d, + self.n_e_1d, self.t_e_1d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0) + + interpolators_element2 = interpolators1d_from_elementdensity(self.atomic_data, self.element2, self.psin_1d, + self.n_element2_1d, + self.n_e_1d, self.t_e_1d, + tcx_donor=self.tcx_donor, + tcx_donor_n=self.n_tcx_donor_1d, + tcx_donor_charge=0) + + mapper = equilibrium_map3d_match_plasma_neutrality(self.atomic_data, self.element_bulk, equilibrium, self.psin_1d, + [interpolators_element1, interpolators_element2], + self.n_e_1d, self.t_e_1d, self.tcx_donor, + self.n_tcx_donor_1d, 0) + + inlcfs = np.where(self.psin_1d < 1)[0] + psin_1d = self.psin_1d[inlcfs] + n = np.zeros_like(psin_1d) + r = np.zeros_like(psin_1d) + for index, value in enumerate(psin_1d): + n[index] = self.n_e_1d(value) + r[index] = equilibrium.psin_to_r(value) + + profile1 = self.evaluate_interpolators(interpolators_element1, psin_1d) + profile2 = self.evaluate_interpolators(interpolators_element2, psin_1d) + profile3 = self.evaluate_interpolators(mapper, r) + + total = self.sumup_electrons(profile1) + total += self.sumup_electrons(profile2) + total += self.sumup_electrons(profile3) + + self.assertTrue(np.allclose(total, n, rtol=self.TOLERANCE)) + + +if __name__ == "__main__": + unittest.main() diff --git a/demos/equilibrium/equilibrium.py b/demos/equilibrium/equilibrium.py index 55b64925..d962a35a 100644 --- a/demos/equilibrium/equilibrium.py +++ b/demos/equilibrium/equilibrium.py @@ -29,7 +29,8 @@ import matplotlib.pyplot as plt -from cherab.core.math import sample2d, Slice3D, Interpolate1DCubic +from raysect.core.math.function.float import Interpolator1DArray +from cherab.core.math import sample2d, Slice3D from cherab.tools.equilibrium import example_equilibrium, plot_equilibrium @@ -107,7 +108,7 @@ # In this example we interpolate the temperature data manually to produce a continuous function # and then map that function around the flux surfaces to give the same result -te_psin = Interpolate1DCubic([0, 0.5, 0.9, 1.0], [2500, 2000, 1000, 0]) +te_psin = Interpolator1DArray([0, 0.5, 0.9, 1.0], [2500, 2000, 1000, 0], 'cubic', 'none', 0) # map to produce 2D and 3D temperature profiles temperature_2d = equilibrium.map2d(te_psin) diff --git a/demos/plasmas/ionisation_balance_1d.py b/demos/plasmas/ionisation_balance_1d.py index 35250b05..de7b8274 100644 --- a/demos/plasmas/ionisation_balance_1d.py +++ b/demos/plasmas/ionisation_balance_1d.py @@ -3,9 +3,9 @@ import matplotlib._color_data as mcd import matplotlib.pyplot as plt import numpy as np +from raysect.core.math.function.float import Interpolator1DArray from cherab.core.atomic import neon, hydrogen, helium -from cherab.core.math import Interpolate1DCubic from cherab.openadas import OpenADAS from cherab.tools.plasmas.ionisation_balance import (fractional_abundance, interpolators1d_fractional, from_elementdensity, @@ -67,13 +67,13 @@ def exp_decay(r, lamb, max_val): n_element2_profile = double_parabola(psin_1d, 5e17, 1e17, 2, 2) n_tcx_donor_profile = exp_decay(psin_1d, 10, 3e16) -t_e = Interpolate1DCubic(psin_1d, t_e_profile) -n_e = Interpolate1DCubic(psin_1d, n_e_profile) +t_e = Interpolator1DArray(psin_1d, t_e_profile, 'cubic', 'none', 0) +n_e = Interpolator1DArray(psin_1d, n_e_profile, 'cubic', 'none', 0) -t_element = Interpolate1DCubic(psin_1d, t_element_profile) -n_element = Interpolate1DCubic(psin_1d, n_element_profile) -n_element2 = Interpolate1DCubic(psin_1d, n_element2_profile) -n_tcx_donor = Interpolate1DCubic(psin_1d, n_tcx_donor_profile) +t_element = Interpolator1DArray(psin_1d, t_element_profile, 'cubic', 'none', 0) +n_element = Interpolator1DArray(psin_1d, n_element_profile, 'cubic', 'none', 0) +n_element2 = Interpolator1DArray(psin_1d, n_element2_profile, 'cubic', 'none', 0) +n_tcx_donor = Interpolator1DArray(psin_1d, n_tcx_donor_profile, 'cubic', 'none', 0) # load adas atomic database and define elements adas = OpenADAS(permit_extrapolation=True) @@ -96,7 +96,7 @@ def exp_decay(r, lamb, max_val): ax.plot(psin_1d, abundance_fractional_profile_tcx[key], "--", label="{0} {1}+ (tcx)".format(element.symbol, key), color=colors[key]) -ax.legend(loc=6) +ax.legend(loc=6, ncol=2) ax.set_xlabel("$\Psi_n$") ax.set_ylabel("fractional abundance [a.u.]") plt.title('Fractional Abundance VS $\Psi_n$') @@ -116,7 +116,7 @@ def exp_decay(r, lamb, max_val): ax.plot(psin_1d, density_element_profiles_tcx[key], "--", label="{0} {1}+ (tcx)".format(element.symbol, key), color=colors[key]) -ax.legend(loc=6) +ax.legend(loc=2, ncol=2) ax.set_xlabel("$\Psi_n$") ax.set_ylabel("ion density [m$^{-3}]$") @@ -155,7 +155,7 @@ def exp_decay(r, lamb, max_val): ax.plot(psin_1d, n_e_profile, "kx", label="input n_e") ax.plot(psin_1d, n_e_recalculated, "k-", label="recalculated n_e") -ax.legend(loc=6) +ax.legend(loc=1, ncol=2) ax.set_xlabel("$\Psi_n$") ax.set_ylabel("ion density [m$^{-3}]$") @@ -177,3 +177,5 @@ def exp_decay(r, lamb, max_val): n_e, t_e, tcx_donor=donor_element, tcx_donor_n=n_tcx_donor, tcx_donor_charge=0) + +plt.show() diff --git a/demos/plasmas/ionisation_balance_2d.py b/demos/plasmas/ionisation_balance_2d.py index 849fa7e7..7a8a9ca9 100644 --- a/demos/plasmas/ionisation_balance_2d.py +++ b/demos/plasmas/ionisation_balance_2d.py @@ -2,8 +2,8 @@ import matplotlib.pyplot as plt import numpy as np +from raysect.core.math.function.float import Interpolator1DArray, Interpolator2DArray from cherab.core.atomic import neon, hydrogen, helium -from cherab.core.math import Interpolate1DCubic, Interpolate2DCubic from cherab.openadas import OpenADAS from cherab.tools.equilibrium import example_equilibrium from cherab.tools.plasmas.ionisation_balance import (fractional_abundance, equilibrium_map3d_fractional, @@ -106,13 +106,13 @@ def exp_decay(r, lamb, max_val, lim=None): n_element2_profile_1d = double_parabola(psin_1d, 5e17, 1e17, 2, 2, 1) n_tcx_donor_profile_1d = exp_decay(psin_1d, 10, 3e16, 1) -t_e_1d = Interpolate1DCubic(psin_1d, t_e_profile_1d) -n_e_1d = Interpolate1DCubic(psin_1d, n_e_profile_1d) +t_e_1d = Interpolator1DArray(psin_1d, t_e_profile_1d, 'cubic', 'none', 0) +n_e_1d = Interpolator1DArray(psin_1d, n_e_profile_1d, 'cubic', 'none', 0) -t_element_1d = Interpolate1DCubic(psin_1d, t_element_profile_1d) -n_element_1d = Interpolate1DCubic(psin_1d, n_element_profile_1d) -n_element2_1d = Interpolate1DCubic(psin_1d, n_element2_profile_1d) -n_tcx_donor_1d = Interpolate1DCubic(psin_1d, n_tcx_donor_profile_1d) +t_element_1d = Interpolator1DArray(psin_1d, t_element_profile_1d, 'cubic', 'none', 0) +n_element_1d = Interpolator1DArray(psin_1d, n_element_profile_1d, 'cubic', 'none', 0) +n_element2_1d = Interpolator1DArray(psin_1d, n_element2_profile_1d, 'cubic', 'none', 0) +n_tcx_donor_1d = Interpolator1DArray(psin_1d, n_tcx_donor_profile_1d, 'cubic', 'none', 0) psin_2d = np.zeros(equilibrium.psi_data.shape) @@ -124,37 +124,37 @@ def exp_decay(r, lamb, max_val, lim=None): for index in np.ndindex(*t_e_profile_2d.shape): t_e_profile_2d[index] = t_e_1d(psin_2d[index]) -t_e_2d = Interpolate2DCubic(equilibrium.r_data, equilibrium.z_data, t_e_profile_2d) +t_e_2d = Interpolator2DArray(equilibrium.r_data, equilibrium.z_data, t_e_profile_2d, 'cubic', 'none', 0, 0) n_e_profile_2d = np.zeros_like(psin_2d) for index in np.ndindex(*n_e_profile_2d.shape): n_e_profile_2d[index] = n_e_1d(psin_2d[index]) -n_e_2d = Interpolate2DCubic(equilibrium.r_data, equilibrium.z_data, n_e_profile_2d) +n_e_2d = Interpolator2DArray(equilibrium.r_data, equilibrium.z_data, n_e_profile_2d, 'cubic', 'none', 0, 0) t_element_profile_2d = np.zeros_like(psin_2d) for index in np.ndindex(*t_element_profile_2d.shape): t_element_profile_2d[index] = t_element_1d(psin_2d[index]) -t_element_2d = Interpolate2DCubic(equilibrium.r_data, equilibrium.z_data, t_element_profile_2d) +t_element_2d = Interpolator2DArray(equilibrium.r_data, equilibrium.z_data, t_element_profile_2d, 'cubic', 'none', 0, 0) n_element_profile_2d = np.zeros_like(psin_2d) for index in np.ndindex(*n_element_profile_2d.shape): n_element_profile_2d[index] = n_element_1d(psin_2d[index]) -n_element_2d = Interpolate2DCubic(equilibrium.r_data, equilibrium.z_data, n_element_profile_2d) +n_element_2d = Interpolator2DArray(equilibrium.r_data, equilibrium.z_data, n_element_profile_2d, 'cubic', 'none', 0, 0) n_element2_profile_2d = np.zeros_like(psin_2d) for index in np.ndindex(*n_element2_profile_2d.shape): n_element2_profile_2d[index] = n_element2_1d(psin_2d[index]) -n_element2_2d = Interpolate2DCubic(equilibrium.r_data, equilibrium.z_data, n_element2_profile_2d) +n_element2_2d = Interpolator2DArray(equilibrium.r_data, equilibrium.z_data, n_element2_profile_2d, 'cubic', 'none', 0, 0) n_tcx_donor_profile_2d = np.zeros_like(psin_2d) for index in np.ndindex(*n_tcx_donor_profile_2d.shape): n_tcx_donor_profile_2d[index] = n_tcx_donor_1d(psin_2d[index]) -n_tcx_donor_2d = Interpolate2DCubic(equilibrium.r_data, equilibrium.z_data, n_element_profile_2d) +n_tcx_donor_2d = Interpolator2DArray(equilibrium.r_data, equilibrium.z_data, n_element_profile_2d, 'cubic', 'none', 0, 0) ######################################################################################################################## @@ -252,6 +252,7 @@ def exp_decay(r, lamb, max_val, lim=None): ax.contour(equilibrium.r_data, equilibrium.z_data, equilibrium.psi_data.T, colors="white") ax.plot(equilibrium.limiter_polygon[:, 0], equilibrium.limiter_polygon[:, 1], "k-") ax.plot(equilibrium.lcfs_polygon[:, 0], equilibrium.lcfs_polygon[:, 1], "r-") - ax.set_title("{} {}+".format(element.symbol, key)) + ax.set_title("{} {}+".format(element_bulk.symbol, key)) ax.set_aspect(1) +plt.show()