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/__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/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..e9950577 --- /dev/null +++ b/cherab/core/model/beam/thin_beam.pyx @@ -0,0 +1,329 @@ +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 + # there must be plasma present to configure + if not self._plasma: + return + # there must be atomic_data present to configure + if not self._atomic_data: + return + + # 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 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")