From 58a3a9ff8bfbd9ede25c2d560c0a8153fac39c87 Mon Sep 17 00:00:00 2001 From: vsnever Date: Tue, 30 Jul 2024 12:49:24 +0200 Subject: [PATCH 1/5] Add the power radiated in spectral lines due to charge exchange with thermal neutral hydrogen to the TotalRadiatedPower model. --- CHANGELOG.md | 1 + cherab/core/atomic/rates.pxd | 2 +- cherab/core/atomic/rates.pyx | 2 +- .../model/plasma/total_radiated_power.pxd | 4 +- .../model/plasma/total_radiated_power.pyx | 66 ++++++++++++++++++- cherab/openadas/rates/radiated_power.pyx | 12 ++-- .../models/radiated_power/radiated_power.rst | 2 +- 7 files changed, 77 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c70837fb..91e0ba21 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ New: * **Beam dispersion calculation has changed from sigma(z) = sigma + z * tan(alpha) to sigma(z) = sqrt(sigma^2 + (z * tan(alpha))^2) for consistancy with the Gaussian beam model. Attention!!! The results of BES and CX spectroscopy are affected by this change. (#414)** * Improved beam direction calculation to allow for natural broadening of the BES line shape due to beam divergence. (#414) * Add kwargs to invert_regularised_nnls to pass them to scipy.optimize.nnls. (#438) +* Add the power radiated in spectral lines due to charge exchange with thermal neutral hydrogen to the TotalRadiatedPower model. (#370) Bug fixes: * Fix deprecated transforms being cached in LaserMaterial after laser.transform update (#420) diff --git a/cherab/core/atomic/rates.pxd b/cherab/core/atomic/rates.pxd index 1f844122..27cec455 100644 --- a/cherab/core/atomic/rates.pxd +++ b/cherab/core/atomic/rates.pxd @@ -84,7 +84,7 @@ cdef class _RadiatedPower: readonly Element element readonly int charge - cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999 + cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999 cdef class LineRadiationPower(_RadiatedPower): diff --git a/cherab/core/atomic/rates.pyx b/cherab/core/atomic/rates.pyx index f3a653db..852725a7 100644 --- a/cherab/core/atomic/rates.pyx +++ b/cherab/core/atomic/rates.pyx @@ -263,7 +263,7 @@ cdef class _RadiatedPower: """ return self.evaluate(electron_density, electron_temperature) - cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: + cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: """ Evaluate the radiated power at the given plasma conditions. diff --git a/cherab/core/model/plasma/total_radiated_power.pxd b/cherab/core/model/plasma/total_radiated_power.pxd index 1900d965..59499ba2 100644 --- a/cherab/core/model/plasma/total_radiated_power.pxd +++ b/cherab/core/model/plasma/total_radiated_power.pxd @@ -18,7 +18,7 @@ from cherab.core.atomic.elements cimport Element -from cherab.core.atomic.rates cimport LineRadiationPower, ContinuumPower +from cherab.core.atomic.rates cimport LineRadiationPower, ContinuumPower, CXRadiationPower from cherab.core.plasma cimport PlasmaModel from cherab.core.species cimport Species @@ -30,7 +30,9 @@ cdef class TotalRadiatedPower(PlasmaModel): Element _element int _charge Species _line_rad_species, _recom_species + list _hydrogen_species LineRadiationPower _plt_rate ContinuumPower _prb_rate + CXRadiationPower _prc_rate cdef int _populate_cache(self) except -1 diff --git a/cherab/core/model/plasma/total_radiated_power.pyx b/cherab/core/model/plasma/total_radiated_power.pyx index 28350d4e..12deda94 100644 --- a/cherab/core/model/plasma/total_radiated_power.pyx +++ b/cherab/core/model/plasma/total_radiated_power.pyx @@ -20,9 +20,40 @@ from raysect.optical cimport Spectrum, Point3D, Vector3D from cherab.core cimport Plasma, AtomicData from cherab.core.utility.constants cimport RECIP_4_PI +from cherab.core.atomic.elements import hydrogen, deuterium, tritium cdef class TotalRadiatedPower(PlasmaModel): + r""" + Emitter that calculates total power radiated by a given ion, which includes: + + - line power due to electron impact excitation, + - continuum and line power due to recombination and Bremsstrahlung, + - line power due to charge exchange with thermal neutral hydrogen and its isotopes. + + The emission calculated by this model is spectrally unresolved, + which means that the total radiated power will be spread of the entire + observable spectral range. + + .. math:: + \epsilon_{\mathrm{total}} = \frac{1}{4 \pi \Delta\lambda} \left( + n_{Z_\mathrm{i}} n_\mathrm{e} C_{\mathrm{excit}}(n_\mathrm{e}, T_\mathrm{e}) + + n_{Z_\mathrm{i} + 1} n_\mathrm{e} C_{\mathrm{recomb}}(n_\mathrm{e}, T_\mathrm{e}) + + n_{Z_\mathrm{i} + 1} n_\mathrm{hyd} C_{\mathrm{cx}}(n_\mathrm{e}, T_\mathrm{e}) \right) + + where :math:`n_{Z_\mathrm{i}}` is the target species density; + :math:`n_{Z_\mathrm{i} + 1}` is the recombining species density; + :math:`n_{\mathrm{hyd}}` is the total density of all hydrogen isotopes; + :math:`C_{\mathrm{excit}}, C_{\mathrm{recomb}}, C_{\mathrm{cx}}` are the radiated power + coefficients in :math:`W m^3` due to electron impact excitation, recombination + + Bremsstrahlung and charge exchange with thermal neutral hydrogen, respectively; + :math:`\Delta\lambda` is the observable spectral range. + + :param Element element: The atomic element/isotope. + :param int charge: The charge state of the element/isotope. + :param Plasma plasma: The plasma to which this emission model is attached. Default is None. + :param AtomicData atomic_data: The atomic data provider for this model. Default is None. + """ def __init__(self, Element element, int charge, Plasma plasma=None, AtomicData atomic_data=None): @@ -42,7 +73,9 @@ cdef class TotalRadiatedPower(PlasmaModel): cdef: int i - double ne, ni, ni_upper, te, plt_radiance, prb_radiance + double ne, ni, ni_upper, nhyd, te + double plt_radiance, prb_radiance, prc_radiance + Species hyd_species # cache data on first run if not self._cache_loaded: @@ -60,22 +93,37 @@ cdef class TotalRadiatedPower(PlasmaModel): ni_upper = self._recom_species.distribution.density(point.x, point.y, point.z) + nhyd = 0 + for hyd_species in self._hydrogen_species: + nhyd += hyd_species.distribution.density(point.x, point.y, point.z) + # add emission to spectrum if self._plt_rate and ni > 0: plt_radiance = RECIP_4_PI * self._plt_rate.evaluate(ne, te) * ne * ni / (spectrum.max_wavelength - spectrum.min_wavelength) else: plt_radiance = 0 + if self._prb_rate and ni_upper > 0: prb_radiance = RECIP_4_PI * self._prb_rate.evaluate(ne, te) * ne * ni_upper / (spectrum.max_wavelength - spectrum.min_wavelength) else: prb_radiance = 0 + + if self._prc_rate and ni_upper > 0 and nhyd > 0: + prc_radiance = RECIP_4_PI * self._prc_rate.evaluate(ne, te) * nhyd * ni_upper / (spectrum.max_wavelength - spectrum.min_wavelength) + else: + prc_radiance = 0 + for i in range(spectrum.bins): - spectrum.samples_mv[i] += plt_radiance + prb_radiance + spectrum.samples_mv[i] += plt_radiance + prb_radiance + prc_radiance return spectrum cdef int _populate_cache(self) except -1: + cdef: + Species hyd_species + Element hyd_isotope + # sanity checks if self._plasma is None: raise RuntimeError("The emission model is not connected to a plasma object.") @@ -100,6 +148,18 @@ cdef class TotalRadiatedPower(PlasmaModel): "recombination/continuum emission, (element={}, ionisation={})." "".format(self._element.symbol, self._charge+1)) + # cache hydrogen species and CX radiation rate + self._prc_rate = self._atomic_data.cx_radiated_power_rate(self._element, self._charge+1) + + self._hydrogen_species = [] + for hyd_isotope in (hydrogen, deuterium, tritium): + try: + hyd_species = self._plasma.composition.get(hyd_isotope, 0) + except ValueError: + pass + else: + self._hydrogen_species.append(hyd_species) + self._cache_loaded = True def _change(self): @@ -108,5 +168,7 @@ cdef class TotalRadiatedPower(PlasmaModel): self._cache_loaded = False self._line_rad_species = None self._recom_species = None + self._hydrogen_species = None self._plt_rate = None self._prb_rate = None + self._prc_rate = None diff --git a/cherab/openadas/rates/radiated_power.pyx b/cherab/openadas/rates/radiated_power.pyx index 570ced92..8e232a4e 100644 --- a/cherab/openadas/rates/radiated_power.pyx +++ b/cherab/openadas/rates/radiated_power.pyx @@ -46,7 +46,7 @@ cdef class LineRadiationPower(CoreLineRadiationPower): 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: + cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: # need to handle zeros, also density and temperature can become negative due to cubic interpolation if electron_density < 1.e-300: @@ -65,7 +65,7 @@ cdef class NullLineRadiationPower(CoreLineRadiationPower): Needed for use cases where the required atomic data is missing. """ - cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: + cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: return 0.0 @@ -92,7 +92,7 @@ cdef class ContinuumPower(CoreContinuumPower): 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: + cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: # need to handle zeros, also density and temperature can become negative due to cubic interpolation if electron_density < 1.e-300: @@ -111,7 +111,7 @@ cdef class NullContinuumPower(CoreContinuumPower): Needed for use cases where the required atomic data is missing. """ - cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: + cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: return 0.0 @@ -137,7 +137,7 @@ cdef class CXRadiationPower(CoreCXRadiationPower): 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: + cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: # need to handle zeros, also density and temperature can become negative due to cubic interpolation if electron_density < 1.e-300: @@ -156,5 +156,5 @@ cdef class NullCXRadiationPower(CoreCXRadiationPower): Needed for use cases where the required atomic data is missing. """ - cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: + cpdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: return 0.0 diff --git a/docs/source/models/radiated_power/radiated_power.rst b/docs/source/models/radiated_power/radiated_power.rst index 0b8594e1..d00fcec2 100644 --- a/docs/source/models/radiated_power/radiated_power.rst +++ b/docs/source/models/radiated_power/radiated_power.rst @@ -2,4 +2,4 @@ Total Radiated Power ==================== -Documentation for this model will go here soon... +.. autoclass:: cherab.core.model.plasma.total_radiated_power.TotalRadiatedPower From 1236237c659e95d0e60601924877fc5051b46573 Mon Sep 17 00:00:00 2001 From: vsnever Date: Tue, 30 Jul 2024 12:50:33 +0200 Subject: [PATCH 2/5] Add the TestCase for the TotalRadiatedPower model. --- .../core/tests/test_total_radiated_power.py | 142 ++++++++++++++++++ 1 file changed, 142 insertions(+) create mode 100644 cherab/core/tests/test_total_radiated_power.py diff --git a/cherab/core/tests/test_total_radiated_power.py b/cherab/core/tests/test_total_radiated_power.py new file mode 100644 index 00000000..2744e803 --- /dev/null +++ b/cherab/core/tests/test_total_radiated_power.py @@ -0,0 +1,142 @@ +# Copyright 2016-2023 Euratom +# Copyright 2016-2023 United Kingdom Atomic Energy Authority +# Copyright 2016-2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +import unittest + +import numpy as np + +from raysect.core import Point3D, Vector3D +from raysect.optical import Ray, World + +from cherab.core.atomic import AtomicData, LineRadiationPower, ContinuumPower, CXRadiationPower +from cherab.core.atomic import deuterium, hydrogen, nitrogen +from cherab.tools.plasmas.slab import build_constant_slab_plasma +from cherab.core.model import TotalRadiatedPower + +from cherab.core.utility import EvAmuToMS, EvToJ + + +class ConstantLineRadiationPower(LineRadiationPower): + """ + Constant line radiation power coefficient. + """ + + def __init__(self, value): + self.value = value + + def evaluate(self, density, temperature): + + return self.value + + +class ConstantContinuumPower(ContinuumPower): + """ + Constant continuum power coefficient. + """ + + def __init__(self, value): + self.value = value + + def evaluate(self, density, temperature): + + return self.value + + +class ConstantCXRadiationPower(CXRadiationPower): + """ + Constant charge exchange radiation power coefficient. + """ + + def __init__(self, value): + self.value = value + + def evaluate(self, density, temperature): + + return self.value + + +class MockAtomicData(AtomicData): + """Fake atomic data for test purpose.""" + + def line_radiated_power_rate(self, element, charge): + + return ConstantLineRadiationPower(1.e-32) + + def continuum_radiated_power_rate(self, element, charge): + + return ConstantContinuumPower(1.e-33) + + def cx_radiated_power_rate(self, element, charge): + + return ConstantCXRadiationPower(1.e-31) + + +class TestTotalRadiatedPower(unittest.TestCase): + + def setUp(self): + + self.world = World() + + self.atomic_data = MockAtomicData() + + plasma_species = [(deuterium, 0, 1.e18, 500., Vector3D(0, 0, 0)), + (hydrogen, 0, 1.e18, 500., Vector3D(0, 0, 0)), + (nitrogen, 6, 5.e18, 1100., Vector3D(0, 0, 0)), + (nitrogen, 7, 1.e19, 1100., Vector3D(0, 0, 0))] + self.plasma = build_constant_slab_plasma(length=1.2, width=1.2, height=1.2, + electron_density=1e19, electron_temperature=1000., + plasma_species=plasma_species) + self.plasma.parent = self.world + self.plasma.atomic_data = self.atomic_data + + def test_beam_density(self): + + self.plasma.models = [TotalRadiatedPower(nitrogen, 6)] + + # observing + origin = Point3D(1.5, 0, 0) + direction = Vector3D(-1, 0, 0) + ray = Ray(origin=origin, direction=direction, + min_wavelength=500., max_wavelength=550., bins=2) + radiated_power = ray.trace(self.world).total() + + # validating + ne = self.plasma.electron_distribution.density(0.5, 0.5, 0.5) + n_n6 = self.plasma.composition[(nitrogen, 6)].distribution.density(0.5, 0.5, 0.5) + n_n7 = self.plasma.composition[(nitrogen, 7)].distribution.density(0.5, 0.5, 0.5) + n_h0 = self.plasma.composition[(hydrogen, 0)].distribution.density(0.5, 0.5, 0.5) + n_d0 = self.plasma.composition[(deuterium, 0)].distribution.density(0.5, 0.5, 0.5) + + integration_length = 1.2 + + plt_rate = self.atomic_data.line_radiated_power_rate(nitrogen, 6).value + plt_radiance = 0.25 / np.pi * plt_rate * ne * n_n6 * integration_length + + prb_rate = self.atomic_data.continuum_radiated_power_rate(nitrogen, 7).value + prb_radiance = 0.25 / np.pi * prb_rate * ne * n_n7 * integration_length + + prc_rate = self.atomic_data.cx_radiated_power_rate(nitrogen, 7).value + prc_radiance = 0.25 / np.pi * prc_rate * (n_h0 + n_d0) * n_n7 * integration_length + + test_radiated_power = plt_radiance + prb_radiance + prc_radiance + + self.assertAlmostEqual(radiated_power / test_radiated_power, 1., delta=1e-8) + + +if __name__ == '__main__': + unittest.main() From 554a02f2b2f404c8c836fbe9efd5e433ce0a7df3 Mon Sep 17 00:00:00 2001 From: vsnever Date: Tue, 30 Jul 2024 23:35:22 +0200 Subject: [PATCH 3/5] Move power density -> spectral radiance conversion to a separate line. --- .../core/model/plasma/total_radiated_power.pyx | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/cherab/core/model/plasma/total_radiated_power.pyx b/cherab/core/model/plasma/total_radiated_power.pyx index 12deda94..fbf1a491 100644 --- a/cherab/core/model/plasma/total_radiated_power.pyx +++ b/cherab/core/model/plasma/total_radiated_power.pyx @@ -74,7 +74,7 @@ cdef class TotalRadiatedPower(PlasmaModel): cdef: int i double ne, ni, ni_upper, nhyd, te - double plt_radiance, prb_radiance, prc_radiance + double plt_power_density, prb_power_density, prc_power_density, radiance Species hyd_species # cache data on first run @@ -99,22 +99,24 @@ cdef class TotalRadiatedPower(PlasmaModel): # add emission to spectrum if self._plt_rate and ni > 0: - plt_radiance = RECIP_4_PI * self._plt_rate.evaluate(ne, te) * ne * ni / (spectrum.max_wavelength - spectrum.min_wavelength) + plt_power_density = self._plt_rate.evaluate(ne, te) * ne * ni else: - plt_radiance = 0 + plt_power_density = 0 if self._prb_rate and ni_upper > 0: - prb_radiance = RECIP_4_PI * self._prb_rate.evaluate(ne, te) * ne * ni_upper / (spectrum.max_wavelength - spectrum.min_wavelength) + prb_power_density = self._prb_rate.evaluate(ne, te) * ne * ni_upper else: - prb_radiance = 0 + prb_power_density = 0 if self._prc_rate and ni_upper > 0 and nhyd > 0: - prc_radiance = RECIP_4_PI * self._prc_rate.evaluate(ne, te) * nhyd * ni_upper / (spectrum.max_wavelength - spectrum.min_wavelength) + prc_power_density = self._prc_rate.evaluate(ne, te) * nhyd * ni_upper else: - prc_radiance = 0 + prc_power_density = 0 + + radiance = RECIP_4_PI * (plt_power_density + prb_power_density + prc_power_density) / (spectrum.max_wavelength - spectrum.min_wavelength) for i in range(spectrum.bins): - spectrum.samples_mv[i] += plt_radiance + prb_radiance + prc_radiance + spectrum.samples_mv[i] += radiance return spectrum From 3400bef492e6b56b772256b34c94ba249e6fea00 Mon Sep 17 00:00:00 2001 From: vsnever Date: Tue, 30 Jul 2024 23:46:48 +0200 Subject: [PATCH 4/5] Tidy up the code for the TotalRadiatedPower. --- .../model/plasma/total_radiated_power.pyx | 30 ++++++++----------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/cherab/core/model/plasma/total_radiated_power.pyx b/cherab/core/model/plasma/total_radiated_power.pyx index fbf1a491..30167d26 100644 --- a/cherab/core/model/plasma/total_radiated_power.pyx +++ b/cherab/core/model/plasma/total_radiated_power.pyx @@ -74,7 +74,7 @@ cdef class TotalRadiatedPower(PlasmaModel): cdef: int i double ne, ni, ni_upper, nhyd, te - double plt_power_density, prb_power_density, prc_power_density, radiance + double power_density, radiance Species hyd_species # cache data on first run @@ -98,22 +98,18 @@ cdef class TotalRadiatedPower(PlasmaModel): nhyd += hyd_species.distribution.density(point.x, point.y, point.z) # add emission to spectrum - if self._plt_rate and ni > 0: - plt_power_density = self._plt_rate.evaluate(ne, te) * ne * ni - else: - plt_power_density = 0 - - if self._prb_rate and ni_upper > 0: - prb_power_density = self._prb_rate.evaluate(ne, te) * ne * ni_upper - else: - prb_power_density = 0 - - if self._prc_rate and ni_upper > 0 and nhyd > 0: - prc_power_density = self._prc_rate.evaluate(ne, te) * nhyd * ni_upper - else: - prc_power_density = 0 - - radiance = RECIP_4_PI * (plt_power_density + prb_power_density + prc_power_density) / (spectrum.max_wavelength - spectrum.min_wavelength) + power_density = 0 + + if self._plt_rate and ni > 0: # excitation + power_density += self._plt_rate.evaluate(ne, te) * ne * ni + + if self._prb_rate and ni_upper > 0: # recombination + bremsstrahlung + power_density += self._prb_rate.evaluate(ne, te) * ne * ni_upper + + if self._prc_rate and ni_upper > 0 and nhyd > 0: # charge exchange + power_density += self._prc_rate.evaluate(ne, te) * nhyd * ni_upper + + radiance = RECIP_4_PI * power_density / (spectrum.max_wavelength - spectrum.min_wavelength) for i in range(spectrum.bins): spectrum.samples_mv[i] += radiance From fd4449bc2e8c4138012e315195dcd8250827f45d Mon Sep 17 00:00:00 2001 From: vsnever Date: Thu, 1 Aug 2024 13:51:15 +0200 Subject: [PATCH 5/5] Replace 'composition' Python property with get_composition() Cython getter in TotalRadiatedPower.populate_cache(). --- cherab/core/model/plasma/total_radiated_power.pyx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cherab/core/model/plasma/total_radiated_power.pyx b/cherab/core/model/plasma/total_radiated_power.pyx index 30167d26..84b42424 100644 --- a/cherab/core/model/plasma/total_radiated_power.pyx +++ b/cherab/core/model/plasma/total_radiated_power.pyx @@ -131,7 +131,7 @@ cdef class TotalRadiatedPower(PlasmaModel): # cache line radiation species and rate self._plt_rate = self._atomic_data.line_radiated_power_rate(self._element, self._charge) try: - self._line_rad_species = self._plasma.composition.get(self._element, self._charge) + self._line_rad_species = self._plasma.get_composition().get(self._element, self._charge) except ValueError: raise RuntimeError("The plasma object does not contain the required ion species for calculating" "total line radiaton, (element={}, ionisation={})." @@ -140,7 +140,7 @@ cdef class TotalRadiatedPower(PlasmaModel): # cache recombination species and radiation rate self._prb_rate = self._atomic_data.continuum_radiated_power_rate(self._element, self._charge+1) try: - self._recom_species = self._plasma.composition.get(self._element, self._charge+1) + self._recom_species = self._plasma.get_composition().get(self._element, self._charge+1) except ValueError: raise RuntimeError("The plasma object does not contain the required ion species for calculating" "recombination/continuum emission, (element={}, ionisation={})." @@ -152,7 +152,7 @@ cdef class TotalRadiatedPower(PlasmaModel): self._hydrogen_species = [] for hyd_isotope in (hydrogen, deuterium, tritium): try: - hyd_species = self._plasma.composition.get(hyd_isotope, 0) + hyd_species = self._plasma.get_composition().get(hyd_isotope, 0) except ValueError: pass else: