From 332965964c669c3ad447d3370de2bc0a6d989db3 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Wed, 26 Oct 2022 16:52:45 +0200 Subject: [PATCH 1/2] First working example BeamDistribution Added BeamDistribution class, with DistributionFunction as baseclass. BeamDistribution handles spatial properties of beam, e.g. geometry, density distribution etc. (these were removed from Beam node). BeamDistribution will allow BeamModels to adress arbitrary level of complexity, from the most complicated 6D integration during observation. Simpler realisation of BeamModels (the currently available ones) can use mehods of the class providing the moments of the distribution: bulk_velocity, temperature, density. SingleRayAttenuator Was made more stand-alone. It now has _transform and _transform_inv attributes which are transforms from and to beam frame. This is to allow simple re-using of SingleRayAttenuator for also "MultiRayAttenuator" BeamDistribuions ThinBeam is a replacement of the current beam model. BeamModels and BeamMaterial changed to work with the ThinBeam model and changes brought by the introduction of the BeamDistribution --- cherab/core/beam/__init__.pxd | 1 + cherab/core/beam/__init__.py | 1 + cherab/core/beam/distribution.pxd | 34 +++ cherab/core/beam/distribution.pyx | 89 ++++++ cherab/core/beam/material.pxd | 5 + cherab/core/beam/material.pyx | 23 +- cherab/core/beam/model.pxd | 8 +- cherab/core/beam/model.pyx | 65 ++-- cherab/core/beam/node.pxd | 32 +- cherab/core/beam/node.pyx | 321 +++----------------- cherab/core/distribution.pxd | 4 +- cherab/core/model/attenuator/singleray.pyx | 53 ++-- cherab/core/model/beam/beam_emission.pyx | 24 +- cherab/core/model/beam/charge_exchange.pyx | 31 +- cherab/core/model/beam/thin_beam.pxd | 37 +++ cherab/core/model/beam/thin_beam.pyx | 331 +++++++++++++++++++++ cherab/core/model/lineshape.pxd | 5 +- cherab/core/model/lineshape.pyx | 22 +- 18 files changed, 691 insertions(+), 395 deletions(-) create mode 100644 cherab/core/beam/distribution.pxd create mode 100644 cherab/core/beam/distribution.pyx create mode 100644 cherab/core/model/beam/thin_beam.pxd create mode 100644 cherab/core/model/beam/thin_beam.pyx diff --git a/cherab/core/beam/__init__.pxd b/cherab/core/beam/__init__.pxd index 6d6d6b0e..2007fadd 100644 --- a/cherab/core/beam/__init__.pxd +++ b/cherab/core/beam/__init__.pxd @@ -18,3 +18,4 @@ from cherab.core.beam.node cimport Beam from cherab.core.beam.model cimport BeamModel, BeamAttenuator +from cherab.core.beam.distribution cimport BeamDistribution \ No newline at end of file diff --git a/cherab/core/beam/__init__.py b/cherab/core/beam/__init__.py index 81d98ce3..1a66bfec 100644 --- a/cherab/core/beam/__init__.py +++ b/cherab/core/beam/__init__.py @@ -18,3 +18,4 @@ from .node import Beam from .model import BeamModel, BeamAttenuator +from cherab.core.beam.distribution import BeamDistribution diff --git a/cherab/core/beam/distribution.pxd b/cherab/core/beam/distribution.pxd new file mode 100644 index 00000000..6d384ec2 --- /dev/null +++ b/cherab/core/beam/distribution.pxd @@ -0,0 +1,34 @@ +# 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. + +from cherab.core.atomic cimport Element +from cherab.core.distribution cimport DistributionFunction +from cherab.core.beam.node cimport Beam + +cdef class BeamDistribution(DistributionFunction): + + + cdef: + Beam _beam + Element _element + + cdef Element get_element(self) + + cpdef Beam get_beam(self) + + cpdef list get_geometry(self) diff --git a/cherab/core/beam/distribution.pyx b/cherab/core/beam/distribution.pyx new file mode 100644 index 00000000..357545af --- /dev/null +++ b/cherab/core/beam/distribution.pyx @@ -0,0 +1,89 @@ +# 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. + + +from cherab.core.beam.node cimport Beam +from cherab.core.atomic cimport Element +from cherab.core.distribution cimport DistributionFunction, ZeroDistribution +from cherab.core.utility import Notifier + + +cdef class BeamDistribution(DistributionFunction): + + def __init__(self): + + super().__init__() + self._beam = None + + @property + def element(self): + return self._element + + @element.setter + def element(self, Element value not None): + self._element = value + self._element_changed() + self.notifier.notify() + + cdef Element get_element(self): + return self._element + + @property + def beam(self): + return self._beam + + @beam.setter + def beam(self, Beam value): + + self._beam = value + self._beam_changed() + + self.notifier.notify() + + cpdef Beam get_beam(self): + return self._beam + + cpdef list get_geometry(self): + """ + Get list of Primitives forming the beam geometry + """ + raise NotImplementedError('Virtual method must be implemented in a sub-class.') + + def _beam_changed(self): + """ + Reaction to _beam changes + + Virtual method call. + """ + pass + + def _element_changed(self): + """ + Reaction to _element change + + Virtual method call. + """ + pass + + def _modified(self): + """ + Called when distribution chages + + Virtual method call. + """ + self.notifier.notify() diff --git a/cherab/core/beam/material.pxd b/cherab/core/beam/material.pxd index 858dd421..65195476 100644 --- a/cherab/core/beam/material.pxd +++ b/cherab/core/beam/material.pxd @@ -16,11 +16,14 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. + +from raysect.optical cimport Primitive from raysect.optical.material.emitter cimport InhomogeneousVolumeEmitter from cherab.core.plasma cimport Plasma from cherab.core.beam cimport Beam from cherab.core.atomic cimport AtomicData +from cherab.core.beam.distribution cimport BeamDistribution cdef class BeamMaterial(InhomogeneousVolumeEmitter): @@ -29,5 +32,7 @@ cdef class BeamMaterial(InhomogeneousVolumeEmitter): Beam _beam Plasma _plasma AtomicData _atomic_data + Primitive _beam_primitive + BeamDistribution _distribution list _models diff --git a/cherab/core/beam/material.pyx b/cherab/core/beam/material.pyx index ce878285..cfff2187 100644 --- a/cherab/core/beam/material.pyx +++ b/cherab/core/beam/material.pyx @@ -20,18 +20,27 @@ from raysect.optical cimport World, Primitive, Ray, Spectrum, SpectralFunction, from raysect.optical.material.emitter cimport InhomogeneousVolumeEmitter from raysect.optical.material.emitter.inhomogeneous cimport VolumeIntegrator from cherab.core.beam.model cimport BeamModel +from cherab.core.beam.distribution cimport BeamDistribution cdef class BeamMaterial(InhomogeneousVolumeEmitter): - def __init__(self, Beam beam not None, Plasma plasma not None, AtomicData atomic_data not None, - list models not None, VolumeIntegrator integrator not None): + def __init__(self, + Beam beam not None, + BeamDistribution distribution not None, + Primitive beam_primitive not None, + Plasma plasma not None, + AtomicData atomic_data not None, + list models not None, + VolumeIntegrator integrator not None): super().__init__(integrator) self._beam = beam + self._distribution = distribution self._plasma = plasma self._atomic_data = atomic_data + self._beam_primitive = beam_primitive # validate for model in models: @@ -42,6 +51,7 @@ cdef class BeamMaterial(InhomogeneousVolumeEmitter): for model in models: model.beam = beam model.plasma = plasma + model.distribution = distribution model.atomic_data = atomic_data self._models = models @@ -53,20 +63,17 @@ cdef class BeamMaterial(InhomogeneousVolumeEmitter): cdef: BeamModel model Point3D plasma_point - Vector3D beam_direction, observation_direction - - beam_direction = self._beam.direction(point.x, point.y, point.z) + Vector3D observation_direction # transform points and directions # todo: cache this transform and rebuild if beam or plasma notifies - beam_to_plasma = self._beam.to(self._plasma) + beam_to_plasma = self._beam_primitive.to(self._plasma) plasma_point = point.transform(beam_to_plasma) - beam_direction = beam_direction.transform(beam_to_plasma) observation_direction = direction.transform(beam_to_plasma) # call each model and accumulate spectrum for model in self._models: - spectrum = model.emission(point, plasma_point, beam_direction, observation_direction, spectrum) + spectrum = model.emission(point, plasma_point, observation_direction, spectrum) return spectrum diff --git a/cherab/core/beam/model.pxd b/cherab/core/beam/model.pxd index 765b41c1..0f9be569 100644 --- a/cherab/core/beam/model.pxd +++ b/cherab/core/beam/model.pxd @@ -16,10 +16,12 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. +from raysect.core cimport AffineMatrix3D from raysect.optical cimport Spectrum, Point3D, Vector3D from cherab.core.plasma.node cimport Plasma from cherab.core.beam.node cimport Beam +from cherab.core.beam.distribution cimport BeamDistribution from cherab.core.atomic cimport AtomicData @@ -29,10 +31,11 @@ cdef class BeamModel: Plasma _plasma Beam _beam AtomicData _atomic_data + BeamDistribution _distribution cdef object __weakref__ - cpdef Spectrum emission(self, Point3D beam_point, Point3D plasma_point, Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum) + cpdef Spectrum emission(self, Point3D beam_point, Point3D plasma_point, Vector3D observation_direction, Spectrum spectrum) cdef class BeamAttenuator: @@ -40,8 +43,9 @@ cdef class BeamAttenuator: cdef: readonly object notifier Plasma _plasma - Beam _beam + BeamDistribution _distribution AtomicData _atomic_data + AffineMatrix3D _transform, _transform_inv cdef object __weakref__ diff --git a/cherab/core/beam/model.pyx b/cherab/core/beam/model.pyx index 4c84508d..115ad121 100644 --- a/cherab/core/beam/model.pyx +++ b/cherab/core/beam/model.pyx @@ -17,7 +17,7 @@ # under the Licence. from cherab.core.utility import Notifier - +from raysect.core cimport AffineMatrix3D cdef class BeamModel: @@ -71,6 +71,24 @@ cdef class BeamModel: # inform model source data has changed self._change() + @property + def distribution(self): + return self._distribution + + @distribution.setter + def distribution(self, BeamDistribution value not None): + + # disconnect from previous beam's notifications + if self._distribution: + self._distribution.notifier.remove(self._change) + + # attach to beam to inform model of changes to beam properties + self._distribution = value + self._distribution.notifier.add(self._change) + + # inform model source data has changed + self._change() + @property def atomic_data(self): return self._atomic_data @@ -83,7 +101,7 @@ cdef class BeamModel: # inform model source data has changed self._change() - cpdef Spectrum emission(self, Point3D beam_point, Point3D plasma_point, Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum): + cpdef Spectrum emission(self, Point3D beam_point, Point3D plasma_point, Vector3D observation_direction, Spectrum spectrum): """ Calculate the emission for a point on the beam in a specified direction. @@ -115,23 +133,36 @@ cdef class BeamModel: cdef class BeamAttenuator: - def __init__(self, Beam beam=None, Plasma plasma=None, AtomicData atomic_data=None): + def __init__(self, BeamDistribution distribution=None, Plasma plasma=None, AtomicData atomic_data=None): - # must notify beam if the attenuator properties change, affecting the density values + # must notify distribution if the attenuator properties change, affecting the density values self.notifier = Notifier() # configure - self._beam = beam + self._distribution = distribution self._plasma = plasma self._atomic_data = atomic_data + self._transform = AffineMatrix3D() + self._transform_inv = AffineMatrix3D() + # setup property change notifications for plasma if self._plasma: self._plasma.notifier.add(self._change) - # setup property change notifications for beam - if self._beam: - self._beam.notifier.add(self._change) + # setup property change notifications for distribution + if self._distribution: + self._distribution.notifier.add(self._change) + + @property + def transform(self): + return self._transform + + @transform.setter + def transform(self, AffineMatrix3D value): + self._transform = value + self._transform_inv = value.inverse() + self._change() @property def plasma(self): @@ -152,19 +183,19 @@ cdef class BeamAttenuator: self._change() @property - def beam(self): - return self._beam + def distribution(self): + return self._distribution - @beam.setter - def beam(self, Beam value not None): + @distribution.setter + def distribution(self, BeamDistribution value not None): # disconnect from previous beam's notifications - if self._beam: - self._beam.notifier.remove(self._change) + if self._distribution: + self._distribution.notifier.remove(self._change) - # attach to beam to inform model of changes to beam properties - self._beam = value - self._beam.notifier.add(self._change) + # attach to distribution to inform model of changes to distribution properties + self._distribution = value + self._distribution.notifier.add(self._change) # inform model source data has changed self._change() diff --git a/cherab/core/beam/node.pxd b/cherab/core/beam/node.pxd index d8159bbb..4c79f975 100644 --- a/cherab/core/beam/node.pxd +++ b/cherab/core/beam/node.pxd @@ -25,6 +25,7 @@ from cherab.core.atomic cimport AtomicData, Element from cherab.core.plasma cimport Plasma from cherab.core.beam.model cimport BeamAttenuator from cherab.core.beam.model cimport BeamModel +from cherab.core.beam.distribution cimport BeamDistribution cdef class ModelManager: @@ -44,40 +45,13 @@ cdef class Beam(Node): cdef: object notifier - Vector3D BEAM_AXIS - double _energy, _power, _temperature - Element _element - double _divergence_x, _divergence_y - double _length, _sigma Plasma _plasma + BeamDistribution _distribution AtomicData _atomic_data ModelManager _models - BeamAttenuator _attenuator - Primitive _geometry + list _geometry VolumeIntegrator _integrator cdef object __weakref__ - cpdef double density(self, double x, double y, double z) except? -1e999 - - cpdef Vector3D direction(self, double x, double y, double z) - - cdef double get_energy(self) - - cdef double get_power(self) - - cdef double get_temperature(self) - - cdef Element get_element(self) - - cdef double get_divergence_x(self) - - cdef double get_divergence_y(self) - - cdef double get_length(self) - - cdef double get_sigma(self) - cdef Plasma get_plasma(self) - - cdef int _modified(self) except -1 diff --git a/cherab/core/beam/node.pyx b/cherab/core/beam/node.pyx index 36877421..eeffb61e 100644 --- a/cherab/core/beam/node.pyx +++ b/cherab/core/beam/node.pyx @@ -28,6 +28,7 @@ from raysect.optical.material.emitter.inhomogeneous cimport NumericalIntegrator from cherab.core.beam.material cimport BeamMaterial from cherab.core.beam.model cimport BeamModel +from cherab.core.beam.distribution cimport BeamDistribution from cherab.core.atomic cimport AtomicData, Element from cherab.core.utility import Notifier from libc.math cimport tan, M_PI @@ -180,188 +181,18 @@ cdef class Beam(Node): # change reporting and tracking self.notifier = Notifier() - # beam properties - self.BEAM_AXIS = Vector3D(0.0, 0.0, 1.0) - self._energy = 0.0 # eV/amu - self._power = 0.0 # total beam power, W - self._temperature = 0.0 # Broadening of the beam (eV) - self._element = element = None # beam species, an Element object - self._divergence_x = 0.0 # beam divergence x (degrees) - self._divergence_y = 0.0 # beam divergence y (degrees) - self._length = 1.0 # m - self._sigma = 0.1 # m (gaussian beam width at origin) - # external data dependencies self._plasma = None self._atomic_data = None + self._distribution = None # setup emission model handler and trigger geometry rebuilding if the models change self._models = ModelManager() self._models.notifier.add(self._configure_geometry) - # beam attenuation model - self._attenuator = None - - # beam geometry - self._geometry = None - # emission model integrator self._integrator = NumericalIntegrator(step=0.001) - cpdef double density(self, double x, double y, double z) except? -1e999: - """ - Returns the bean density at the specified position in beam coordinates. - - Note: this function is only defined over the domain 0 < z < beam_length. - Outside of this range the density is clamped to zero. - - :param x: x coordinate in meters. - :param y: y coordinate in meters. - :param z: z coordinate in meters. - :return: Beam density in m^-3 - """ - - if self._attenuator is None: - raise ValueError('The beam must have an attenuator model to provide density values.') - - # todo: make z > length return 0 as a non-default toggle. It should throw an error by default to warn users they are requesting data outside the domain. - if z < 0 or z > self._length: - return 0 - - return self._attenuator.density(x, y, z) - - cpdef Vector3D direction(self, double x, double y, double z): - """ - Calculates the beam direction vector at a point in space. - - Note the values of the beam outside of the beam envelope should be - treated with caution. - - :param x: x coordinate in meters. - :param y: y coordinate in meters. - :param z: z coordinate in meters. - :return: Direction vector. - """ - - # if behind the beam just return the beam axis (for want of a better value) - if z <= 0: - return self.BEAM_AXIS - - # calculate direction from divergence - cdef double dx = tan(DEGREES_TO_RADIANS * self._divergence_x) - cdef double dy = tan(DEGREES_TO_RADIANS * self._divergence_y) - return new_vector3d(dx, dy, 1.0).normalise() - - @property - def energy(self): - return self._energy - - @energy.setter - def energy(self, double value): - if value < 0: - raise ValueError('Beam energy cannot be less than zero.') - self._energy = value - self.notifier.notify() - - cdef double get_energy(self): - return self._energy - - @property - def power(self): - return self._power - - @power.setter - def power(self, double value): - if value < 0: - raise ValueError('Beam power cannot be less than zero.') - self._power = value - self.notifier.notify() - - cdef double get_power(self): - return self._power - - @property - def temperature(self): - return self._temperature - - @temperature.setter - def temperature(self, double value): - if value < 0: - raise ValueError('Beam temperature cannot be less than zero.') - self._temperature = value - self.notifier.notify() - - cdef double get_temperature(self): - return self._temperature - - @property - def element(self): - return self._element - - @element.setter - def element(self, Element value not None): - self._element = value - self.notifier.notify() - - cdef Element get_element(self): - return self._element - - @property - def divergence_x(self): - return self._divergence_x - - @divergence_x.setter - def divergence_x(self, double value): - if value < 0: - raise ValueError('Beam x divergence cannot be less than zero.') - self._divergence_x = value - self.notifier.notify() - - cdef double get_divergence_x(self): - return self._divergence_x - - @property - def divergence_y(self): - return self._divergence_y - - @divergence_y.setter - def divergence_y(self, double value): - if value < 0: - raise ValueError('Beam y divergence cannot be less than zero.') - self._divergence_y = value - self.notifier.notify() - - cdef double get_divergence_y(self): - return self._divergence_y - - @property - def length(self): - return self._length - - @length.setter - def length(self, double value): - if value <= 0: - raise ValueError('Beam length must be greater than zero.') - self._length = value - self.notifier.notify() - - cdef double get_length(self): - return self._length - - @property - def sigma(self): - return self._sigma - - @sigma.setter - def sigma(self, double value): - if value <= 0: - raise ValueError('Beam sigma (width) must be greater than zero.') - self._sigma = value - self.notifier.notify() - - cdef double get_sigma(self): - return self._sigma - @property def atomic_data(self): return self._atomic_data @@ -369,8 +200,7 @@ cdef class Beam(Node): @atomic_data.setter def atomic_data(self, AtomicData value not None): self._atomic_data = value - self._configure_geometry() - self._configure_attenuator() + self.notifier.notify() @property def plasma(self): @@ -379,39 +209,11 @@ cdef class Beam(Node): @plasma.setter def plasma(self, Plasma value not None): self._plasma = value - self._configure_geometry() - self._configure_attenuator() + self.notifier.notify() cdef Plasma get_plasma(self): return self._plasma - @property - def attenuator(self): - return self._attenuator - - @attenuator.setter - def attenuator(self, BeamAttenuator value not None): - - # check necessary data is available - if not self._plasma: - raise ValueError('The beam must have a reference to a plasma object to be used with an attenuator.') - - if not self._atomic_data: - raise ValueError('The beam must have an atomic data source to be used with an emission model.') - - # disconnect from previous attenuator's notifications - if self._attenuator: - self._attenuator.notifier.remove(self._modified) - - self._attenuator = value - self._configure_attenuator() - - # connect to new attenuator's notifications - self._attenuator.notifier.add(self._modified) - - # attenuator supplies beam density, notify dependents there is a data change - self.notifier.notify() - @property def models(self): return self._models @@ -423,8 +225,8 @@ cdef class Beam(Node): if not self._plasma: raise ValueError('The beam must have a reference to a plasma object to be used with an emission model.') - if not self._attenuator: - raise ValueError('The beam must have an attenuator model to be used with an emission model.') + if not self._distribution: + raise ValueError('The beam must have a distribution to be used with an emission model.') if not self._atomic_data: raise ValueError('The beam must have an atomic data source to be used with an emission model.') @@ -441,7 +243,24 @@ cdef class Beam(Node): def integrator(self, VolumeIntegrator value): self._integrator = value self._configure_geometry() + + @property + def distribution(self): + return self._distribution + + @distribution.setter + def distribution(self, BeamDistribution value): + # stop notifications of the old distribution + if self._distribution: + self.notifier.remove(self._distribution._beam_changed) + + self._distribution = value + self._distribution.beam = self + self.notifier.add(self._distribution._beam_changed) + self._configure_geometry() + self.notifier.notify() + def _configure_geometry(self): # detach existing geometry @@ -457,93 +276,25 @@ cdef class Beam(Node): # check necessary data is available if not self._plasma: raise ValueError('The beam must have a reference to a plasma object to be used with an emission model.') - - if not self._attenuator: - raise ValueError('The beam must have an attenuator model to be used with an emission model.') + + if not self._distribution: + raise ValueError('The beam must have a distribution to be used with an emission model.') if not self._atomic_data: raise ValueError('The beam must have an atomic data source to be used with an emission model.') # build geometry to fit beam - self._geometry = self._generate_geometry() + self._geometry = self._distribution.get_geometry() # attach geometry to the beam - self._geometry.parent = self - self._geometry.name = 'Beam Geometry' - - # add plasma material - self._geometry.material = BeamMaterial(self, self._plasma, self._atomic_data, list(self._models), self.integrator) - - def _generate_geometry(self): - """ - Generate the bounding geometry for the beam model. - - Where possible the beam is bound by a cone as this offers the tightest - fitting bounding volume. To avoid numerical issues caused by creating - extremely long cones in low divergence cases, the geometry is switched - to a cylinder where the difference in volume between the cone and a - cylinder is less than 10%. - - :return: Beam geometry Primitive. - """ - - # number of beam sigma the bounding volume lies from the beam axis - num_sigma = self._attenuator.clamp_sigma - - # return Cylinder(NUM_SIGMA * self.sigma, height=self.length) - - # no divergence, use a cylinder - if self._divergence_x == 0 and self._divergence_y == 0: - return Cylinder(num_sigma * self.sigma, height=self.length) - - # rate of change of beam radius with z (using largest divergence) - drdz = tan(DEGREES_TO_RADIANS * max(self._divergence_x, self._divergence_y)) - - # radii of bounds at the beam origin (z=0) and the beam end (z=length) - radius_start = num_sigma * self.sigma - radius_end = radius_start + self.length * num_sigma * drdz - - # distance of the cone apex to the beam origin - distance_apex = radius_start / (num_sigma * drdz) - cone_height = self.length + distance_apex - - # calculate volumes - cylinder_volume = self.length * M_PI * radius_end**2 - cone_volume = M_PI * (cone_height * radius_end**2 - distance_apex * radius_start**2) / 3 - volume_ratio = cone_volume / cylinder_volume - - # if the volume difference is <10%, generate a cylinder - if volume_ratio > 0.9: - return Cylinder(num_sigma * self.sigma, height=self.length) - - # cone has to be rotated by 180 deg and shifted by beam length in the +z direction - cone_transform = translate(0, 0, self.length) * rotate_x(180) - - # create cone and cut off -z protrusion - return Intersect( - Cone(radius_end, cone_height, transform=cone_transform), - Cylinder(radius_end * 1.01, self.length * 1.01) - ) - - def _configure_attenuator(self): - - # there must be an attenuator present to configure - if not self._attenuator: - return - - # check necessary data is available - if not self._plasma: - raise ValueError('The beam must have a reference to a plasma object to be used with an attenuator.') - - if not self._atomic_data: - raise ValueError('The beam must have an atomic data source to be used with an emission model.') - - # setup attenuator - self._attenuator.beam = self - self._attenuator.plasma = self._plasma - self._attenuator.atomic_data = self._atomic_data - - cdef int _modified(self) except -1: + for segment in self._geometry: + segment.parent = self + # add plasma material + segment.material = BeamMaterial(self, self._distribution, segment, + self._plasma, self._atomic_data, + list(self._models), self.integrator) + + def _modified(self): """ Called when a scene-graph change occurs that modifies this Node's root transforms. This will occur if the Node's transform is modified, a @@ -552,4 +303,4 @@ cdef class Beam(Node): """ # beams section of the scene-graph has been modified, alert dependents - self.notifier.notify() + self.notifier.notify() \ No newline at end of file diff --git a/cherab/core/distribution.pxd b/cherab/core/distribution.pxd index 062d7a45..324eb87e 100644 --- a/cherab/core/distribution.pxd +++ b/cherab/core/distribution.pxd @@ -23,7 +23,9 @@ from cherab.core.math cimport Function3D, VectorFunction3D cdef class DistributionFunction: - cdef object notifier + cdef: + object __weakref__ + object notifier cdef double evaluate(self, double x, double y, double z, double vx, double vy, double vz) except? -1e999 diff --git a/cherab/core/model/attenuator/singleray.pyx b/cherab/core/model/attenuator/singleray.pyx index a2cfbde1..a0f9329e 100644 --- a/cherab/core/model/attenuator/singleray.pyx +++ b/cherab/core/model/attenuator/singleray.pyx @@ -27,7 +27,8 @@ 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.plasma cimport Plasma -from cherab.core.beam cimport Beam +from cherab.core.beam.node cimport Beam +from cherab.core.beam.distribution cimport BeamDistribution from cherab.core.species cimport Species from cherab.core.utility.constants cimport DEGREES_TO_RADIANS @@ -38,9 +39,10 @@ cimport cython # todo: attenuation calculation could be optimised further using memory views etc... cdef class SingleRayAttenuator(BeamAttenuator): - def __init__(self, double step=0.01, bint clamp_to_zero=False, double clamp_sigma=5.0, Beam beam=None, Plasma plasma=None, AtomicData atomic_data=None): + def __init__(self, double step=0.01, bint clamp_to_zero=False, double clamp_sigma=5.0, + BeamDistribution distribution=None, Plasma plasma=None, AtomicData atomic_data=None): - super().__init__(beam, plasma, atomic_data) + super().__init__(distribution, plasma, atomic_data) self._source_density = 0. self._density = None @@ -91,13 +93,18 @@ cdef class SingleRayAttenuator(BeamAttenuator): The point is specified in beam space. - :param x: x coordinate in meters. - :param y: y coordinate in meters. - :param z: z coordinate in meters. + :param x: x coordinate in the beam space in meters. + :param y: y coordinate in the beam space in meters. + :param z: z coordinate in the beam space in meters. :return: Density in m^-3. """ - cdef double sigma_x, sigma_y, norm_radius_sqr, gaussian_sample + cdef: + double sigma_x, sigma_y, norm_radius_sqr, gaussian_sample + Point3D att_pnt + + # transform x, y, z from beam into attenuator's reference frame + att_pnt = new_point3d(x, y, z).transform(self._transform_inv) # use cached data if available if self._stopping_data is None: @@ -107,11 +114,11 @@ cdef class SingleRayAttenuator(BeamAttenuator): self._calc_attenuation() # calculate beam width - sigma_x = self._beam.get_sigma() + z * self._tanxdiv - sigma_y = self._beam.get_sigma() + z * self._tanydiv + sigma_x = self._distribution.get_sigma() + att_pnt.z * self._tanxdiv + sigma_y = self._distribution.get_sigma() + att_pnt.z * self._tanydiv # normalised radius squared - norm_radius_sqr = ((x / sigma_x)**2 + (y / sigma_y)**2) + norm_radius_sqr = ((att_pnt.x / sigma_x)**2 + (att_pnt.y / sigma_y)**2) # clamp low densities to zero (beam models can skip their calculation if density is zero) # comparison is done using the squared radius to avoid a costly square root @@ -149,15 +156,17 @@ cdef class SingleRayAttenuator(BeamAttenuator): Vector3D direction int nbeam, i np.ndarray beam_z, xaxis, yaxis, zaxis, beam_density + Beam beam double bzv # calculate transform to plasma space - beam_to_plasma = self._beam.to(self._plasma) - direction = self._beam.BEAM_AXIS.transform(beam_to_plasma) + beam = self._distribution.get_beam() + att_to_plasma = beam.to(self._plasma).mul(self._transform) + direction = self._distribution.BEAM_AXIS.transform(att_to_plasma) # sample points along the beam - nbeam = max(1 + int(np.ceil(self._beam.length / self._step)), 4) - beam_z = np.linspace(0.0, self._beam.length, nbeam) + nbeam = max(1 + int(np.ceil(self._distribution.length / self._step)), 4) + beam_z = np.linspace(0.0, self._distribution.length, nbeam) xaxis = np.zeros(nbeam) yaxis = np.zeros(nbeam) @@ -165,17 +174,17 @@ cdef class SingleRayAttenuator(BeamAttenuator): for i, bzv in enumerate(beam_z): - paxis = new_point3d(0.0, 0.0, bzv).transform(beam_to_plasma) + paxis = new_point3d(0.0, 0.0, bzv).transform(att_to_plasma) xaxis[i] = paxis.x yaxis[i] = paxis.y zaxis[i] = paxis.z beam_density = self._beam_attenuation(beam_z, xaxis, yaxis, zaxis, - self._beam.energy, self._beam.power, - self._beam.element.atomic_weight, direction) + self._distribution.energy, self._distribution.power, + self._distribution.element.atomic_weight, direction) - self._tanxdiv = tan(DEGREES_TO_RADIANS * self._beam.divergence_x) - self._tanydiv = tan(DEGREES_TO_RADIANS * self._beam.divergence_y) + self._tanxdiv = tan(DEGREES_TO_RADIANS * self._distribution.divergence_x) + self._tanydiv = tan(DEGREES_TO_RADIANS * self._distribution.divergence_y) # a tiny degree of extrapolation is permitted to handle numerical accuracy issues with the end of the array self._density = Interpolator1DArray(beam_z, beam_density, 'linear', 'nearest', extrapolation_range=1e-9) @@ -207,7 +216,7 @@ cdef class SingleRayAttenuator(BeamAttenuator): naxis = axis.size - speed = EvAmuToMS.to(energy) + speed = self._distribution.get_speed() beam_velocity = direction.normalise() * speed beam_particle_rate = power / EvToJ.to(energy * mass) @@ -282,7 +291,7 @@ cdef class SingleRayAttenuator(BeamAttenuator): BeamStoppingRate stopping_coeff # sanity checks - if not self._beam: + if not self._distribution: raise ValueError("The beam attenuator is not connected to a beam object.") if not self._plasma: @@ -293,7 +302,7 @@ cdef class SingleRayAttenuator(BeamAttenuator): self._stopping_data = [] for species in self._plasma.composition: - stopping_coeff = self._atomic_data.beam_stopping_rate(self._beam.element, species.element, species.charge) + stopping_coeff = self._atomic_data.beam_stopping_rate(self._distribution.element, species.element, species.charge) self._stopping_data.append((species, stopping_coeff)) def _change(self): diff --git a/cherab/core/model/beam/beam_emission.pyx b/cherab/core/model/beam/beam_emission.pyx index 8ca55331..4e6ec5aa 100644 --- a/cherab/core/model/beam/beam_emission.pyx +++ b/cherab/core/model/beam/beam_emission.pyx @@ -96,26 +96,32 @@ cdef class BeamEmissionLine(BeamModel): @cython.wraparound(False) @cython.cdivision(True) @cython.initializedcheck(False) - cpdef Spectrum emission(self, Point3D beam_point, Point3D plasma_point, Vector3D beam_direction, + cpdef Spectrum emission(self, Point3D beam_point, Point3D plasma_point, Vector3D observation_direction, Spectrum spectrum): cdef: + double b_x, b_y, b_z double beam_density, rate, radiance - Vector3D beam_velocity + Vector3D beam_velocity, beam_direction # cache data on first run if self._rates_list is None: self._populate_cache() + # unpack beam point coordinates + b_x = beam_point.x + b_y = beam_point.y + b_z = beam_point.z + # obtain beam density from beam - beam_density = self._beam.density(beam_point.x, beam_point.y, beam_point.z) + beam_density = self._distribution.density(b_x, b_y, b_z) # abort calculation if beam density is zero if beam_density == 0.0: return spectrum # obtain beam velocity - beam_velocity = beam_direction.normalise().mul(evamu_to_ms(self._beam.get_energy())) + beam_velocity = self._distribution.bulk_velocity(b_x, b_y, b_z) # beam emission rate in W rate = self._beam_emission_rate(plasma_point.x, plasma_point.y, plasma_point.z, beam_velocity) @@ -124,7 +130,7 @@ cdef class BeamEmissionLine(BeamModel): radiance = RECIP_4_PI * beam_density * rate return self._lineshape.add_line(radiance, beam_point, plasma_point, - beam_direction, observation_direction, spectrum) + beam_velocity, observation_direction, spectrum) @cython.cdivision(True) cdef double _beam_emission_rate(self, double x, double y, double z, Vector3D beam_velocity) except? -1e999: @@ -165,7 +171,7 @@ cdef class BeamEmissionLine(BeamModel): # calculate mean beam interaction energy interaction_velocity = beam_velocity - target_velocity - interaction_speed = interaction_velocity.length + interaction_speed = interaction_velocity.get_length() interaction_energy = ms_to_evamu(interaction_speed) # species equivalent electron density @@ -185,6 +191,8 @@ cdef class BeamEmissionLine(BeamModel): # sanity checks if self._beam is None: raise RuntimeError("The emission model is not connected to a beam object.") + if self._distribution is None: + raise RuntimeError("The emission model is not connected to a distribution object.") if self._plasma is None: raise RuntimeError("The emission model is not connected to a plasma object.") if self._atomic_data is None: @@ -193,7 +201,7 @@ cdef class BeamEmissionLine(BeamModel): raise RuntimeError("The emission line has not been set.") # check specified emission line is consistent with attached beam object - beam_element = self._beam.element + beam_element = self._distribution.element transition = self._line.transition if beam_element != self._line.element: raise TypeError("The specified line element '{}' is incompatible with the attached neutral " @@ -210,8 +218,8 @@ cdef class BeamEmissionLine(BeamModel): rate = self._atomic_data.beam_emission_pec(beam_element, species.element, species.charge, transition) self._rates_list.append((species, rate)) + self._lineshape = BeamEmissionMultiplet(self._line, self._wavelength, self._beam, self._distribution, self._sigma_to_pi, # instance line shape renderer - self._lineshape = BeamEmissionMultiplet(self._line, self._wavelength, self._beam, self._sigma_to_pi, self._sigma1_to_sigma0, self._pi2_to_pi3, self._pi4_to_pi3) def _change(self): diff --git a/cherab/core/model/beam/charge_exchange.pyx b/cherab/core/model/beam/charge_exchange.pyx index 5cefa3d6..326fd8e9 100644 --- a/cherab/core/model/beam/charge_exchange.pyx +++ b/cherab/core/model/beam/charge_exchange.pyx @@ -79,11 +79,11 @@ cdef class BeamCXLine(BeamModel): @cython.wraparound(False) @cython.cdivision(True) @cython.initializedcheck(False) - cpdef Spectrum emission(self, Point3D beam_point, Point3D plasma_point, Vector3D beam_direction, + cpdef Spectrum emission(self, Point3D beam_point, Point3D plasma_point, Vector3D observation_direction, Spectrum spectrum): cdef: - double x, y, z + double p_x, p_y, p_z double donor_density double receiver_temperature, receiver_density, receiver_ion_mass, interaction_speed, interaction_energy, emission_rate Vector3D receiver_velocity, donor_velocity, interaction_velocity @@ -92,40 +92,45 @@ cdef class BeamCXLine(BeamModel): # cache data on first run if self._target_species is None: self._populate_cache() + + # extract for compact code + b_x = beam_point.x + b_y = beam_point.y + b_z = beam_point.z # obtain donor density from beam - donor_density = self._beam.density(beam_point.x, beam_point.y, beam_point.z) + donor_density = self._distribution.density(b_x, b_y, b_z) # abort calculation if donor density is zero if donor_density == 0.0: return spectrum # extract for more compact code - x = plasma_point.x - y = plasma_point.y - z = plasma_point.z + p_x = plasma_point.x + p_y = plasma_point.y + p_z = plasma_point.z # abort calculation if receiver density is zero - receiver_density = self._target_species.distribution.density(x, y, z) + receiver_density = self._target_species.distribution.density(p_x, p_y, p_z) if receiver_density == 0: return spectrum # abort calculation if receiver temperature is zero - receiver_temperature = self._target_species.distribution.effective_temperature(x, y, z) + receiver_temperature = self._target_species.distribution.effective_temperature(p_x, p_y, p_z) if receiver_temperature == 0: return spectrum - receiver_velocity = self._target_species.distribution.bulk_velocity(x, y, z) + receiver_velocity = self._target_species.distribution.bulk_velocity(p_x, p_y, p_z) receiver_ion_mass = self._target_species.element.atomic_weight - donor_velocity = beam_direction.normalise().mul(evamu_to_ms(self._beam.get_energy())) + donor_velocity = self._distribution.bulk_velocity(b_x, b_y, b_z) interaction_velocity = donor_velocity.sub(receiver_velocity) interaction_speed = interaction_velocity.get_length() interaction_energy = ms_to_evamu(interaction_speed) # calculate the composite charge-exchange emission coefficient - emission_rate = self._composite_cx_rate(x, y, z, interaction_energy, donor_velocity, receiver_temperature, receiver_density) + emission_rate = self._composite_cx_rate(p_x, p_y, p_z, interaction_energy, donor_velocity, receiver_temperature, receiver_density) # calculate emission line central wavelength, doppler shifted along observation direction natural_wavelength = self._wavelength @@ -272,6 +277,8 @@ cdef class BeamCXLine(BeamModel): # sanity checks if self._beam is None: raise RuntimeError("The emission model is not connected to a beam object.") + if self._distribution is None: + raise RuntimeError("The emission model is not connected to a distribution object.") if self._plasma is None: raise RuntimeError("The emission model is not connected to a plasma object.") if self._atomic_data is None: @@ -281,7 +288,7 @@ cdef class BeamCXLine(BeamModel): receiver_element = self._line.element receiver_charge = self._line.charge + 1 - donor_element = self._beam.element + donor_element = self._distribution.element transition = self._line.transition # locate target species diff --git a/cherab/core/model/beam/thin_beam.pxd b/cherab/core/model/beam/thin_beam.pxd new file mode 100644 index 00000000..5669e0c1 --- /dev/null +++ b/cherab/core/model/beam/thin_beam.pxd @@ -0,0 +1,37 @@ +from raysect.core cimport Vector3D +from raysect.optical cimport Primitive +from raysect.optical.material.emitter.inhomogeneous cimport NumericalIntegrator + +from cherab.core.beam.distribution cimport BeamDistribution +from cherab.core.atomic cimport AtomicData, Element +from cherab.core.plasma cimport Plasma +from cherab.core.beam.model cimport BeamAttenuator + +cdef class ThinBeam(BeamDistribution): + + cdef: + readonly Vector3D BEAM_AXIS + double _energy, _power, _temperature, _speed + double _divergence_x, _divergence_y + bint _z_outofbounds + double _length, _sigma + Plasma _plasma + AtomicData _atomic_data + BeamAttenuator _attenuator + Primitive _geometry + + cpdef double get_energy(self) + + cpdef double get_speed(self) + + cpdef double get_power(self) + + cpdef double get_divergence_x(self) + + cpdef double get_divergence_y(self) + + cpdef double get_length(self) + + cpdef double get_sigma(self) + + cpdef Plasma get_plasma(self) \ No newline at end of file diff --git a/cherab/core/model/beam/thin_beam.pyx b/cherab/core/model/beam/thin_beam.pyx new file mode 100644 index 00000000..21430edc --- /dev/null +++ b/cherab/core/model/beam/thin_beam.pyx @@ -0,0 +1,331 @@ +from raysect.core cimport Vector3D, new_vector3d +from raysect.primitive import Cylinder, Cone, Intersect + +from raysect.core cimport translate, rotate_x + +from cherab.core.beam.material cimport BeamMaterial +from cherab.core.beam.distribution cimport BeamDistribution +from cherab.core.atomic cimport AtomicData +from cherab.core.utility import Notifier, EvAmuToMS + +from libc.math cimport tan, M_PI + +cdef double DEGREES_TO_RADIANS = M_PI / 180 + + +cdef class ThinBeam(BeamDistribution): + + def __init__(self): + + super().__init__() + + # beam properties + self.BEAM_AXIS = Vector3D(0.0, 0.0, 1.0) + self._energy = 0.0 # eV/amu + self._power = 0.0 # total beam power, W + self._speed = 0.0 # speed of the particles (m/s) + self._temperature = 0.0 # Broadening of the beam (eV) + self._element = element = None # beam species, an Element object + self._divergence_x = 0.0 # beam divergence x (degrees) + self._divergence_y = 0.0 # beam divergence y (degrees) + self._length = 1.0 # m + self._sigma = 0.1 # m (gaussian beam width at origin) + self._z_outofbounds = False + + cpdef Vector3D bulk_velocity(self, double x, double y, double z): + """ + Evaluates the species' bulk velocity at the specified 3D coordinate. + + :param float x: position in meters + :param float y: position in meters + :param float z: position in meters + :return: velocity vector in m/s + :rtype: Vector3D + """ + # if behind the beam just return the beam axis (for want of a better value) + if z <= 0: + return self.BEAM_AXIS + + # calculate direction from divergence + cdef double dx = tan(DEGREES_TO_RADIANS * self._divergence_x) + cdef double dy = tan(DEGREES_TO_RADIANS * self._divergence_y) + + return new_vector3d(dx, dy, 1.0).normalise().mul(self._speed) + + cpdef double effective_temperature(self, double x, double y, double z) except? -1e999: + """ + Evaluates the species' effective temperature at the specified 3D coordinate. + + :param float x: position in meters + :param float y: position in meters + :param float z: position in meters + :return: temperature in eV + :rtype: float + """ + + return self._temperature + + cpdef double density(self, double x, double y, double z) except? -1e999: + """ + Evaluates the species' density at the specified 3D coordinate. + + :param float x: position in meters + :param float y: position in meters + :param float z: position in meters + :return: density in m^-3 + :rtype: float + """ + + if self._attenuator is None: + raise ValueError('The beam must have an attenuator model to provide density values.') + + # todo: make z > length return 0 as a non-default toggle. It should throw an error by default to warn users they are requesting data outside the domain. + if z < 0 or z > self._length: + if self._z_outofbounds: + return 0 + else: + raise ValueError("z coordinate in beam space out of bounds: z e (0, length)") + + return self._attenuator.density(x, y, z) + + @property + def z_outofbounds(self): + return self._z_outofbounds + + @z_outofbounds.setter + def z_outofbounds(self, bint value): + self._z_outofbounds = value + + @property + def temperature(self): + return self._temperature + + @temperature.setter + def temperature(self, double value): + + if value < 0: + raise ValueError('Temperature cannot be less than zero.') + + self._temperature = value + + @property + def energy(self): + return self._energy + + @energy.setter + def energy(self, double value): + if value < 0: + raise ValueError('Beam energy cannot be less than zero.') + self._energy = value + self.notifier.notify() + self._speed = EvAmuToMS.to(value) + + cpdef double get_energy(self): + return self._energy + + cpdef double get_speed(self): + return self._speed + + @property + def power(self): + return self._power + + @power.setter + def power(self, double value): + if value < 0: + raise ValueError('Beam power cannot be less than zero.') + self._power = value + self.notifier.notify() + + cpdef double get_power(self): + return self._power + + @property + def divergence_x(self): + return self._divergence_x + + @divergence_x.setter + def divergence_x(self, double value): + if value < 0: + raise ValueError('Beam x divergence cannot be less than zero.') + self._divergence_x = value + self.notifier.notify() + + cpdef double get_divergence_x(self): + return self._divergence_x + + @property + def divergence_y(self): + return self._divergence_y + + @divergence_y.setter + def divergence_y(self, double value): + if value < 0: + raise ValueError('Beam y divergence cannot be less than zero.') + self._divergence_y = value + self.notifier.notify() + + cpdef double get_divergence_y(self): + return self._divergence_y + + @property + def length(self): + return self._length + + @length.setter + def length(self, double value): + if value <= 0: + raise ValueError('Beam length must be greater than zero.') + self._length = value + self.notifier.notify() + + cpdef double get_length(self): + return self._length + + @property + def sigma(self): + return self._sigma + + @sigma.setter + def sigma(self, double value): + if value <= 0: + raise ValueError('Beam sigma (width) must be greater than zero.') + self._sigma = value + self.notifier.notify() + + cpdef double get_sigma(self): + return self._sigma + + @property + def atomic_data(self): + return self._atomic_data + + @atomic_data.setter + def atomic_data(self, AtomicData value not None): + self._atomic_data = value + self._configure_attenuator() + + @property + def plasma(self): + return self._plasma + + @plasma.setter + def plasma(self, Plasma value not None): + + self._plasma = value + self._configure_attenuator() + + cpdef Plasma get_plasma(self): + return self._plasma + + @property + def attenuator(self): + return self._attenuator + + @attenuator.setter + def attenuator(self, BeamAttenuator value not None): + + # disconnect from previous attenuator's notifications + if self._attenuator: + self._attenuator.notifier.remove(self._modified) + + self._attenuator = value + self._configure_attenuator() + + # connect to new attenuator's notifications + self._attenuator.notifier.add(self._modified) + + # attenuator supplies beam density, notify dependents there is a data change + self.notifier.notify() + + cpdef list get_geometry(self): + + # check necessary data is available + if not self._plasma: + raise ValueError('The beam must have a reference to a plasma object to be used with an emission model.') + + if not self._attenuator: + raise ValueError('The beam must have an attenuator model to be used with an emission model.') + + if not self._atomic_data: + raise ValueError('The beam must have an atomic data source to be used with an emission model.') + + return [self._generate_geometry()] + + def _generate_geometry(self): + """ + Generate the bounding geometry for the beam model. + + Where possible the beam is bound by a cone as this offers the tightest + fitting bounding volume. To avoid numerical issues caused by creating + extremely long cones in low divergence cases, the geometry is switched + to a cylinder where the difference in volume between the cone and a + cylinder is less than 10%. + + :return: Beam geometry Primitive. + """ + + # number of beam sigma the bounding volume lies from the beam axis + num_sigma = self._attenuator.clamp_sigma + + # return Cylinder(NUM_SIGMA * self.sigma, height=self.length) + + # no divergence, use a cylinder + if self._divergence_x == 0 and self._divergence_y == 0: + return Cylinder(num_sigma * self.sigma, height=self.length) + + # rate of change of beam radius with z (using largest divergence) + drdz = tan(DEGREES_TO_RADIANS * max(self._divergence_x, self._divergence_y)) + + # radii of bounds at the beam origin (z=0) and the beam end (z=length) + radius_start = num_sigma * self.sigma + radius_end = radius_start + self.length * num_sigma * drdz + + # distance of the cone apex to the beam origin + distance_apex = radius_start / (num_sigma * drdz) + cone_height = self.length + distance_apex + + # calculate volumes + cylinder_volume = self.length * M_PI * radius_end**2 + cone_volume = M_PI * (cone_height * radius_end**2 - distance_apex * radius_start**2) / 3 + volume_ratio = cone_volume / cylinder_volume + + # if the volume difference is <10%, generate a cylinder + if volume_ratio > 0.9: + return Cylinder(num_sigma * self.sigma, height=self.length) + + # cone has to be rotated by 180 deg and shifted by beam length in the +z direction + cone_transform = translate(0, 0, self.length) * rotate_x(180) + + # create cone and cut off -z protrusion + return Intersect( + Cone(radius_end, cone_height, transform=cone_transform), + Cylinder(radius_end * 1.01, self.length * 1.01) + ) + + def _configure_attenuator(self): + + # there must be an attenuator present to configure + if not self._attenuator: + return + + # check necessary data is available + if not self._plasma: + raise ValueError('The distribution must have a reference to a plasma object to be used with an attenuator.') + + if not self._atomic_data: + raise ValueError('The distribution must have an atomic data source to be used with an emission model.') + + + # setup attenuator + self._attenuator.distribution = self + self._attenuator.plasma = self._plasma + self._attenuator.atomic_data = self._atomic_data + + def _beam_changed(self): + """ + Reaction to _beam changes + """ + self._plasma = self._beam.plasma + self._atomic_data = self._beam.atomic_data + self._configure_attenuator() \ No newline at end of file diff --git a/cherab/core/model/lineshape.pxd b/cherab/core/model/lineshape.pxd index 711475b5..ce389ebb 100644 --- a/cherab/core/model/lineshape.pxd +++ b/cherab/core/model/lineshape.pxd @@ -25,6 +25,7 @@ from raysect.optical cimport Spectrum, Point3D, Vector3D from cherab.core cimport Line, Species, Plasma, Beam from cherab.core.math cimport Function1D, Function2D from cherab.core.atomic.zeeman cimport ZeemanStructure +from cherab.core.beam.distribution cimport BeamDistribution cpdef double doppler_shift(double wavelength, Vector3D observation_direction, Vector3D velocity) @@ -97,14 +98,16 @@ cdef class BeamLineShapeModel: Line line double wavelength Beam beam + BeamDistribution distribution cpdef Spectrum add_line(self, double radiance, Point3D beam_point, Point3D plasma_point, - Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum) + Vector3D beam_velocity, Vector3D observation_direction, Spectrum spectrum) cdef class BeamEmissionMultiplet(BeamLineShapeModel): cdef: + double _mstoevamu Function2D _sigma_to_pi Function1D _sigma1_to_sigma0, _pi2_to_pi3, _pi4_to_pi3 diff --git a/cherab/core/model/lineshape.pyx b/cherab/core/model/lineshape.pyx index 5e346175..f9147a4e 100644 --- a/cherab/core/model/lineshape.pyx +++ b/cherab/core/model/lineshape.pyx @@ -19,6 +19,7 @@ # under the Licence. import numpy as np +from scipy.constants import atomic_mass, elementary_charge cimport numpy as np from libc.math cimport sqrt, erf, M_SQRT2, floor, ceil, fabs from raysect.optical.spectrum cimport new_spectrum @@ -829,11 +830,12 @@ cdef class BeamLineShapeModel: :param Beam beam: The beam class that is emitting. """ - def __init__(self, Line line, double wavelength, Beam beam): + def __init__(self, Line line, double wavelength, Beam beam, BeamDistribution distribution): self.line = line self.wavelength = wavelength self.beam = beam + self.distribution = distribution cpdef Spectrum add_line(self, double radiance, Point3D beam_point, Point3D plasma_point, Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum): @@ -848,19 +850,20 @@ cdef class BeamEmissionMultiplet(BeamLineShapeModel): Produces Beam Emission Multiplet line shape, also known as the Motional Stark Effect spectrum. """ - def __init__(self, Line line, double wavelength, Beam beam, object sigma_to_pi, - object sigma1_to_sigma0, object pi2_to_pi3, object pi4_to_pi3): + def __init__(self, Line line, double wavelength, Beam beam, BeamDistribution distribution, + object sigma_to_pi, object sigma1_to_sigma0, object pi2_to_pi3, object pi4_to_pi3): - super().__init__(line, wavelength, beam) + super().__init__(line, wavelength, beam, distribution) self._sigma_to_pi = autowrap_function2d(sigma_to_pi) self._sigma1_to_sigma0 = autowrap_function1d(sigma1_to_sigma0) self._pi2_to_pi3 = autowrap_function1d(pi2_to_pi3) self._pi4_to_pi3 = autowrap_function1d(pi4_to_pi3) + self._mstoevamu = atomic_mass / (2 * elementary_charge) @cython.cdivision(True) cpdef Spectrum add_line(self, double radiance, Point3D beam_point, Point3D plasma_point, - Vector3D beam_direction, Vector3D observation_direction, Spectrum spectrum): + Vector3D beam_velocity, Vector3D observation_direction, Spectrum spectrum): cdef double x, y, z cdef Plasma plasma @@ -869,7 +872,7 @@ cdef class BeamEmissionMultiplet(BeamLineShapeModel): cdef double sigma_to_pi, d, intensity_sig, intensity_pi, e_field cdef double s1_to_s0, intensity_s0, intensity_s1 cdef double pi2_to_pi3, pi4_to_pi3, intensity_pi2, intensity_pi3, intensity_pi4 - cdef Vector3D b_field, beam_velocity + cdef Vector3D b_field # extract for more compact code x = plasma_point.x @@ -886,11 +889,10 @@ cdef class BeamEmissionMultiplet(BeamLineShapeModel): if ne <= 0.0: return spectrum - beam_energy = self.beam.get_energy() + beam_energy = beam_velocity.get_length() ** 2 * self._mstoevamu # calculate Stark splitting b_field = plasma.get_b_field().evaluate(x, y, z) - beam_velocity = beam_direction.normalise().mul(evamu_to_ms(beam_energy)) e_field = beam_velocity.cross(b_field).get_length() stark_split = fabs(STARK_SPLITTING_FACTOR * e_field) # TODO - calculate splitting factor? Reject other lines? @@ -899,8 +901,8 @@ cdef class BeamEmissionMultiplet(BeamLineShapeModel): central_wavelength = doppler_shift(natural_wavelength, observation_direction, beam_velocity) # calculate doppler broadening - beam_ion_mass = self.beam.get_element().atomic_weight - beam_temperature = self.beam.get_temperature() + beam_ion_mass = self.distribution.get_element().atomic_weight + beam_temperature = self.distribution.effective_temperature(beam_point.x, beam_point.y, beam_point.z) sigma = thermal_broadening(self.wavelength, beam_temperature, beam_ion_mass) # calculate relative intensities of sigma and pi lines From ac4053f27077d21dc97a82d26737fb7f92d35cd9 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Thu, 27 Oct 2022 09:52:55 +0200 Subject: [PATCH 2/2] Remove error raising from _configure attenuator removing errors should make the ThinBeam initialisation more user-friendly Adjust beam demos to make them work with the new architecture --- cherab/core/model/beam/__init__.py | 1 + cherab/core/model/beam/thin_beam.pyx | 10 ++-- demos/beam.py | 59 +++++++++++++---------- demos/plasmas/beam_into_slab.py | 70 ++++++++++++++++------------ 4 files changed, 80 insertions(+), 60 deletions(-) diff --git a/cherab/core/model/beam/__init__.py b/cherab/core/model/beam/__init__.py index aa1b1cff..fb5c3b45 100644 --- a/cherab/core/model/beam/__init__.py +++ b/cherab/core/model/beam/__init__.py @@ -21,3 +21,4 @@ from .charge_exchange import BeamCXLine from .beam_emission import BeamEmissionLine +from .thin_beam import ThinBeam diff --git a/cherab/core/model/beam/thin_beam.pyx b/cherab/core/model/beam/thin_beam.pyx index 21430edc..e9950577 100644 --- a/cherab/core/model/beam/thin_beam.pyx +++ b/cherab/core/model/beam/thin_beam.pyx @@ -308,14 +308,12 @@ cdef class ThinBeam(BeamDistribution): # there must be an attenuator present to configure if not self._attenuator: return - - # check necessary data is available + # there must be plasma present to configure if not self._plasma: - raise ValueError('The distribution must have a reference to a plasma object to be used with an attenuator.') - + return + # there must be atomic_data present to configure if not self._atomic_data: - raise ValueError('The distribution must have an atomic data source to be used with an emission model.') - + return # setup attenuator self._attenuator.distribution = self diff --git a/demos/beam.py b/demos/beam.py index 49bdc648..848fbe49 100644 --- a/demos/beam.py +++ b/demos/beam.py @@ -31,7 +31,7 @@ from cherab.core import Plasma, Beam, Species, Maxwellian from cherab.core.atomic import elements, Line from cherab.openadas import OpenADAS -from cherab.core.model import SingleRayAttenuator, BeamCXLine +from cherab.core.model import SingleRayAttenuator, BeamCXLine, ThinBeam from cherab.tools.plasmas import GaussianVolume integration_step = 0.02 @@ -74,17 +74,20 @@ plasma.composition = [d_species, he2_species, c6_species, ne10_species] # BEAM ------------------------------------------------------------------------ +distribution = ThinBeam() +distribution.energy = 60000 +distribution.power = 3e6 +distribution.element = elements.deuterium +distribution.sigma = 0.025 +distribution.divergence_x = 0.5 +distribution.divergence_y = 0.5 +distribution.length = 3.0 +distribution.attenuator = SingleRayAttenuator(clamp_to_zero=True) + beam = Beam(parent=world, transform=translate(1.0, 0.0, 0) * rotate(90, 0, 0)) beam.plasma = plasma beam.atomic_data = adas -beam.energy = 60000 -beam.power = 3e6 -beam.element = elements.deuterium -beam.sigma = 0.025 -beam.divergence_x = 0.5 -beam.divergence_y = 0.5 -beam.length = 3.0 -beam.attenuator = SingleRayAttenuator(clamp_to_zero=True) +beam.distribution = distribution beam.models = [ BeamCXLine(Line(elements.helium, 1, (4, 3))), BeamCXLine(Line(elements.helium, 1, (6, 4))), @@ -97,17 +100,20 @@ beam.integrator.step = integration_step beam.integrator.min_samples = 10 +distribution = ThinBeam() +distribution.energy = 60000 / 2 +distribution.power = 3e6 +distribution.element = elements.deuterium +distribution.sigma = 0.025 +distribution.divergence_x = 0.5 +distribution.divergence_y = 0.5 +distribution.length = 3.0 +distribution.attenuator = SingleRayAttenuator(clamp_to_zero=True) + beam = Beam(parent=world, transform=translate(1.0, 0.0, 0) * rotate(90, 0, 0)) beam.plasma = plasma beam.atomic_data = adas -beam.energy = 60000 / 2 -beam.power = 3e6 -beam.element = elements.deuterium -beam.sigma = 0.025 -beam.divergence_x = 0.5 -beam.divergence_y = 0.5 -beam.length = 3.0 -beam.attenuator = SingleRayAttenuator(clamp_to_zero=True) +beam.distribution = distribution beam.models = [ BeamCXLine(Line(elements.helium, 1, (4, 3))), BeamCXLine(Line(elements.helium, 1, (6, 4))), @@ -120,17 +126,20 @@ beam.integrator.step = integration_step beam.integrator.min_samples = 10 +distribution = ThinBeam() +distribution.energy = 60000 / 2 +distribution.power = 3e6 +distribution.element = elements.deuterium +distribution.sigma = 0.025 +distribution.divergence_x = 0.5 +distribution.divergence_y = 0.5 +distribution.length = 3.0 +distribution.attenuator = SingleRayAttenuator(clamp_to_zero=True) + beam = Beam(parent=world, transform=translate(1.0, 0.0, 0) * rotate(90, 0, 0)) beam.plasma = plasma beam.atomic_data = adas -beam.energy = 60000 / 3 -beam.power = 3e6 -beam.element = elements.deuterium -beam.sigma = 0.025 -beam.divergence_x = 0.5 -beam.divergence_y = 0.5 -beam.length = 3.0 -beam.attenuator = SingleRayAttenuator(clamp_to_zero=True) +beam.distribution = distribution beam.models = [ BeamCXLine(Line(elements.helium, 1, (4, 3))), BeamCXLine(Line(elements.helium, 1, (6, 4))), diff --git a/demos/plasmas/beam_into_slab.py b/demos/plasmas/beam_into_slab.py index 373ae7f9..b505b87c 100644 --- a/demos/plasmas/beam_into_slab.py +++ b/demos/plasmas/beam_into_slab.py @@ -9,7 +9,7 @@ from cherab.core import Beam from cherab.core.math import sample3d from cherab.core.atomic import hydrogen, deuterium, carbon, Line -from cherab.core.model import SingleRayAttenuator, BeamCXLine +from cherab.core.model import SingleRayAttenuator, BeamCXLine, ThinBeam from cherab.tools.plasmas.slab import build_slab_plasma from cherab.openadas import OpenADAS @@ -67,47 +67,59 @@ beam_energy = 50000 # keV +# set up full energy component +dist_full = ThinBeam() +dist_full.energy = beam_energy +dist_full.power = 3e6 +dist_full.element = deuterium +dist_full.sigma = 0.05 +dist_full.divergence_x = 0.5 +dist_full.divergence_y = 0.5 +dist_full.length = 3.0 +dist_full.attenuator = SingleRayAttenuator(clamp_to_zero=True) + beam_full = Beam(parent=world, transform=beam_transform) beam_full.plasma = plasma beam_full.atomic_data = adas -beam_full.energy = beam_energy -beam_full.power = 3e6 -beam_full.element = deuterium -beam_full.sigma = 0.05 -beam_full.divergence_x = 0.5 -beam_full.divergence_y = 0.5 -beam_full.length = 3.0 -beam_full.attenuator = SingleRayAttenuator(clamp_to_zero=True) +beam_full.distribution = dist_full beam_full.models = [BeamCXLine(Line(carbon, 5, (8, 7)))] beam_full.integrator.step = integration_step beam_full.integrator.min_samples = 10 +# set up 1/2 energy component +dist_half = ThinBeam() +dist_half.energy = beam_energy / 2 +dist_half.power = 3e6 +dist_half.element = deuterium +dist_half.sigma = 0.05 +dist_half.divergence_x = 0.5 +dist_half.divergence_y = 0.5 +dist_half.length = 3.0 +dist_half.attenuator = SingleRayAttenuator(clamp_to_zero=True) + beam_half = Beam(parent=world, transform=beam_transform) beam_half.plasma = plasma beam_half.atomic_data = adas -beam_half.energy = beam_energy / 2 -beam_half.power = 3e6 -beam_half.element = deuterium -beam_half.sigma = 0.05 -beam_half.divergence_x = 0.5 -beam_half.divergence_y = 0.5 -beam_half.length = 3.0 -beam_half.attenuator = SingleRayAttenuator(clamp_to_zero=True) +beam_half.distribution = dist_half beam_half.models = [BeamCXLine(Line(carbon, 5, (8, 7)))] beam_half.integrator.step = integration_step beam_half.integrator.min_samples = 10 +# set up 1/3 energy component +dist_third = ThinBeam() +dist_third.energy = beam_energy / 3 +dist_third.power = 3e6 +dist_third.element = deuterium +dist_third.sigma = 0.05 +dist_third.divergence_x = 0.5 +dist_third.divergence_y = 0.5 +dist_third.length = 3.0 +dist_third.attenuator = SingleRayAttenuator(clamp_to_zero=True) + beam_third = Beam(parent=world, transform=beam_transform) beam_third.plasma = plasma beam_third.atomic_data = adas -beam_third.energy = beam_energy / 3 -beam_third.power = 3e6 -beam_third.element = deuterium -beam_third.sigma = 0.05 -beam_third.divergence_x = 0.5 -beam_third.divergence_y = 0.5 -beam_third.length = 3.0 -beam_third.attenuator = SingleRayAttenuator(clamp_to_zero=True) +beam_third.distribution = dist_third beam_third.models = [BeamCXLine(Line(carbon, 5, (8, 7)))] beam_third.integrator.step = integration_step beam_third.integrator.min_samples = 10 @@ -118,7 +130,7 @@ plt.figure() -x, _, z, beam_density = sample3d(beam_full.density, (-0.5, 0.5, 200), (0, 0, 1), (0, 3, 200)) +x, _, z, beam_density = sample3d(beam_full.distribution.density, (-0.5, 0.5, 200), (0, 0, 1), (0, 3, 200)) plt.imshow(np.transpose(np.squeeze(beam_density)), extent=[-0.5, 0.5, 0, 3], origin='lower') plt.colorbar() plt.axis('equal') @@ -128,9 +140,9 @@ z = np.linspace(0, 3, 200) -beam_full_densities = [beam_full.density(0, 0, zz) for zz in z] -beam_half_densities = [beam_half.density(0, 0, zz) for zz in z] -beam_third_densities = [beam_third.density(0, 0, zz) for zz in z] +beam_full_densities = [beam_full.distribution.density(0, 0, zz) for zz in z] +beam_half_densities = [beam_half.distribution.density(0, 0, zz) for zz in z] +beam_third_densities = [beam_third.distribution.density(0, 0, zz) for zz in z] plt.figure() plt.plot(z, beam_full_densities, label="full energy") plt.plot(z, beam_half_densities, label="half energy")