From f69b74c8784f6bfda988b71ec27be161e74ade85 Mon Sep 17 00:00:00 2001 From: vsnever Date: Wed, 23 Nov 2022 22:01:12 +0300 Subject: [PATCH 1/9] Add cylindrical and periodic mappers. --- cherab/core/math/__init__.py | 13 +- cherab/core/math/mappers.pxd | 81 +++++-- cherab/core/math/mappers.pyx | 447 ++++++++++++++++++++++++++++++++++- 3 files changed, 514 insertions(+), 27 deletions(-) diff --git a/cherab/core/math/__init__.py b/cherab/core/math/__init__.py index b3e9a031..4f821332 100644 --- a/cherab/core/math/__init__.py +++ b/cherab/core/math/__init__.py @@ -1,6 +1,6 @@ -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -31,6 +31,11 @@ from .caching import Caching1D, Caching2D, Caching3D from .clamp import ClampOutput1D, ClampOutput2D, ClampOutput3D from .clamp import ClampInput1D, ClampInput2D, ClampInput3D -from .mappers import IsoMapper2D, IsoMapper3D, Swizzle2D, Swizzle3D, AxisymmetricMapper, VectorAxisymmetricMapper +from .mappers import IsoMapper2D, IsoMapper3D +from .mappers import Swizzle2D, Swizzle3D +from .mappers import AxisymmetricMapper, VectorAxisymmetricMapper +from .mappers import CylindricalMapper, VectorCylindricalMapper +from .mappers import PeriodicMapper1D, PeriodicMapper2D, PeriodicMapper3D +from .mappers import VectorPeriodicMapper1D, VectorPeriodicMapper2D, VectorPeriodicMapper3D from .mask import PolygonMask2D from .slice import Slice2D, Slice3D diff --git a/cherab/core/math/mappers.pxd b/cherab/core/math/mappers.pxd index 5b455d9d..4fee0cbf 100644 --- a/cherab/core/math/mappers.pxd +++ b/cherab/core/math/mappers.pxd @@ -1,6 +1,6 @@ -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -16,7 +16,11 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. -from cherab.core.math.function cimport Function1D, Function2D, Function3D, VectorFunction2D, VectorFunction3D +from libc.math cimport fmod +from raysect.core.math.function.float cimport Function1D, Function2D, Function3D +from raysect.core.math.function.vector3d cimport Function1D as VectorFunction1D +from raysect.core.math.function.vector3d cimport Function2D as VectorFunction2D +from raysect.core.math.function.vector3d cimport Function3D as VectorFunction3D from raysect.core cimport Vector3D @@ -26,8 +30,6 @@ cdef class IsoMapper2D(Function2D): readonly Function1D function1d readonly Function2D function2d - cdef double evaluate(self, double x, double y) except? -1e999 - cdef class IsoMapper3D(Function3D): @@ -35,15 +37,11 @@ cdef class IsoMapper3D(Function3D): readonly Function3D function3d readonly Function1D function1d - cdef double evaluate(self, double x, double y, double z) except? -1e999 - cdef class Swizzle2D(Function2D): cdef readonly Function2D function2d - cdef double evaluate(self, double x, double y) except? -1e999 - cdef class Swizzle3D(Function3D): @@ -51,18 +49,71 @@ cdef class Swizzle3D(Function3D): readonly Function3D function3d int shape[3] - cdef double evaluate(self, double x, double y, double z) except? -1e999 - cdef class AxisymmetricMapper(Function3D): cdef readonly Function2D function2d - cdef double evaluate(self, double x, double y, double z) except? -1e999 - cdef class VectorAxisymmetricMapper(VectorFunction3D): cdef readonly VectorFunction2D function2d - cdef Vector3D evaluate(self, double x, double y, double z) \ No newline at end of file + +cdef class CylindricalMapper(Function3D): + + cdef readonly Function3D function3d + + +cdef class VectorCylindricalMapper(VectorFunction3D): + + cdef readonly VectorFunction3D function3d + + +cdef inline double remainder(double x1, double x2) nogil: + if x2 == 0: + return x1 + x1 = fmod(x1, x2) + return x1 + x2 if (x1 < 0) else x1 + + +cdef class PeriodicMapper1D(Function1D): + + cdef: + readonly Function1D function1d + readonly double period + + +cdef class PeriodicMapper2D(Function2D): + + cdef: + readonly Function2D function2d + double period_x, period_y + + +cdef class PeriodicMapper3D(Function3D): + + cdef: + readonly Function3D function3d + readonly double period_x, period_y, period_z + + +cdef class VectorPeriodicMapper1D(VectorFunction1D): + + cdef: + readonly VectorFunction1D function1d + readonly double period + + +cdef class VectorPeriodicMapper2D(VectorFunction2D): + + cdef: + readonly VectorFunction2D function2d + readonly double period_x, period_y + + +cdef class VectorPeriodicMapper3D(VectorFunction3D): + + cdef: + readonly VectorFunction3D function3d + readonly double period_x, period_y, period_z diff --git a/cherab/core/math/mappers.pyx b/cherab/core/math/mappers.pyx index 5ff224a3..63f868a9 100644 --- a/cherab/core/math/mappers.pyx +++ b/cherab/core/math/mappers.pyx @@ -20,7 +20,10 @@ from libc.math cimport sqrt, atan2, M_PI -from cherab.core.math.function cimport autowrap_function1d, autowrap_function2d, autowrap_function3d, autowrap_vectorfunction2d +from raysect.core.math.function.float cimport autowrap_function1d, autowrap_function2d, autowrap_function3d +from raysect.core.math.function.vector3d cimport autowrap_function1d as autowrap_vectorfunction1d +from raysect.core.math.function.vector3d cimport autowrap_function2d as autowrap_vectorfunction2d +from raysect.core.math.function.vector3d cimport autowrap_function3d as autowrap_vectorfunction3d from raysect.core cimport rotate_z cimport cython @@ -251,13 +254,10 @@ cdef class AxisymmetricMapper(Function3D): def __init__(self, object function2d): if not callable(function2d): - raise TypeError("Function3D is not callable.") + raise TypeError("Function2D is not callable.") self.function2d = autowrap_function2d(function2d) - def __call__(self, double x, double y, double z): - return self.evaluate(x, y, z) - cdef double evaluate(self, double x, double y, double z) except? -1e999: """Return the value of function2d when it is y-axis symmetrically extended to the 3D space.""" @@ -299,13 +299,11 @@ cdef class VectorAxisymmetricMapper(VectorFunction3D): self.function2d = autowrap_vectorfunction2d(vectorfunction2d) - def __call__(self, double x, double y, double z): - return self.evaluate(x, y, z) - @cython.cdivision(True) cdef Vector3D evaluate(self, double x, double y, double z): """Return the value of function2d when it is y-axis symmetrically extended to the 3D space.""" + cdef double r, phi # convert to cylindrical coordinates phi = atan2(y, x) / M_PI * 180 @@ -313,3 +311,436 @@ cdef class VectorAxisymmetricMapper(VectorFunction3D): # perform axisymmetric rotation return self.function2d.evaluate(r, z).transform(rotate_z(phi)) + + +cdef class CylindricalMapper(Function3D): + """ + Converts Cartesian coordinates to cylindrical coordinates and calls a 3D function + defined in cylindrical coordinates, f(r, :math:`\\phi`, z). + + The angular coordinate is given in radians. + + Positive angular coordinate is measured counterclockwise from the xz plane. + + :param Function3D function3d: The function to be mapped. Must be defined + in the interval (:math:`-\\pi`, :math:`\\pi`] + on the angular axis. + + .. code-block:: pycon + + >>> from math import sqrt, cos + >>> from cherab.core.math import CylindricalMapper + >>> + >>> def my_func(r, phi, z): + >>> return r * cos(phi) + >>> + >>> f = CylindricalMapper(my_func) + >>> + >>> f(1, 0, 0) + 1.0 + >>> f(0.5 * sqrt(3), 0.5, 0) + 0.8660254037844385 + """ + + def __init__(self, object function3d): + + if not callable(function3d): + raise TypeError("Function3D is not callable.") + + self.function3d = autowrap_function3d(function3d) + + cdef double evaluate(self, double x, double y, double z) except? -1e999: + """ + Converts to cylindrical coordinates and evaluates the function + defined in cylindrical coordinates. + """ + cdef double r, phi + + r = sqrt(x * x + y * y) + phi = atan2(y, x) + + return self.function3d.evaluate(r, phi, z) + + +cdef class VectorCylindricalMapper(VectorFunction3D): + """ + Converts Cartesian coordinates to cylindrical coordinates, calls + a 3D vector function defined in cylindrical coordinates, f(r, :math:`\\phi`, z), + then converts the returned 3D vector to Cartesian coordinates. + + The angular coordinate is given in radians. + + Positive angular coordinate is measured counterclockwise from the xz plane. + + :param VectorFunction3D function3d: The function to be mapped. Must be defined + in the interval (:math:`-\\pi`, :math:`\\pi`] + on the angular axis. + + .. code-block:: pycon + + >>> from math import sqrt, cos + >>> from raysect.core.math import Vector3D + >>> from cherab.core.math import VectorCylindricalMapper + >>> + >>> def my_vec_func(r, phi, z): + >>> v = Vector3D(0, 1, 0) + >>> v.length = r * abs(cos(phi)) + >>> return v + >>> + >>> f = VectorCylindricalMapper(my_vec_func) + >>> + >>> f(1, 0, 0) + Vector3D(0.0, 1.0, 0.0) + >>> f(1/sqrt(2), 1/sqrt(2), 0) + Vector3D(-0.5, 0.5, 0.0) + """ + + def __init__(self, object function3d): + + if not callable(function3d): + raise TypeError("Function3D is not callable.") + + self.function3d = autowrap_vectorfunction3d(function3d) + + @cython.cdivision(True) + cdef Vector3D evaluate(self, double x, double y, double z): + """ + Converts to cylindrical coordinates, evaluates the vector function + defined in cylindrical coordinates and rotates the resulting vector + around z-axis. + """ + cdef double r, phi + + r = sqrt(x * x + y * y) + phi = atan2(y, x) + + return self.function3d.evaluate(r, phi, z).transform(rotate_z(phi / M_PI * 180)) + + +cdef class PeriodicMapper1D(Function1D): + """ + Maps a periodic 1D function into 1D space. + + :param Function1D function1d: The periodic 1D function to map defined + in the [0, period) interval. + :param double period: The period of the function. + + .. code-block:: pycon + + >>> from cherab.core.math import PeriodicMapper1D + >>> + >>> def f1(x): + >>> return x + >>> + >>> f2 = PeriodicMapper1D(f1, 1.) + >>> + >>> f2(1.5) + 0.5 + >>> f2(-0.3) + 0.7 + """ + + def __init__(self, object function1d, double period): + + if not callable(function1d): + raise TypeError("function1d is not callable.") + + self.function1d = autowrap_function1d(function1d) + + if period <= 0: + raise ValueError("Argument period must be positive.") + + self.period = period + + cdef double evaluate(self, double x) except? -1e999: + """Return the value of periodic function.""" + + return self.function1d.evaluate(remainder(x, self.period)) + + +cdef class PeriodicMapper2D(Function2D): + """ + Maps a periodic 2D function into 2D space. + + Set period_x/period_y to 0 if the function is not periodic along x/y axis. + + :param Function2D function2d: The periodic 2D function to map defined + in the ([0, period_x), [0, period_y)) intervals. + :param double period_x: The period of the function along x-axis. + 0 if not periodic. + :param double period_y: The period of the function along y-axis. + 0 if not periodic. + + .. code-block:: pycon + + >>> from cherab.core.math import PeriodicMapper2D + >>> + >>> def f1(x, y): + >>> return x * y + >>> + >>> f2 = PeriodicMapper2D(f1, 1., 1.) + >>> + >>> f2(1.5, 1.5) + 0.25 + >>> f2(-0.3, -1.3) + 0.49 + >>> + >>> f3 = PeriodicMapper2D(f1, 1., 0) + >>> + >>> f3(1.5, 1.5) + 0.75 + >>> f3(-0.3, -1.3) + -0.91 + """ + + def __init__(self, object function2d, double period_x, double period_y): + + if not callable(function2d): + raise TypeError("function2d is not callable.") + + self.function2d = autowrap_function2d(function2d) + + if period_x < 0: + raise ValueError("Argument period_x must be >= 0.") + if period_y < 0: + raise ValueError("Argument period_y must be >= 0.") + + self.period_x = period_x + self.period_y = period_y + + cdef double evaluate(self, double x, double y) except? -1e999: + """Return the value of periodic function.""" + + x = remainder(x, self.period_x) + y = remainder(y, self.period_y) + + return self.function2d.evaluate(x, y) + + +cdef class PeriodicMapper3D(Function3D): + """ + Maps a periodic 3D function into 3D space. + + Set period_x/period_y/period_z to 0 if the function is not periodic along x/y/z axis. + + :param Function3D function3d: The periodic 3D function to map defined in the + ([0, period_x), [0, period_y), [0, period_z)) intervals. + :param double period_x: The period of the function along x-axis. + 0 if not periodic. + :param double period_y: The period of the function along y-axis. + 0 if not periodic. + :param double period_z: The period of the function along z-axis. + 0 if not periodic. + + .. code-block:: pycon + + >>> from cherab.core.math import PeriodicMapper3D + >>> + >>> def f1(x, y, z): + >>> return x * y * z + >>> + >>> f2 = PeriodicMapper3D(f1, 1., 1., 1.) + >>> + >>> f2(1.5, 1.5, 1.5) + 0.125 + >>> f2(-0.3, -1.3, -2.3) + 0.343 + >>> + >>> f3 = PeriodicMapper3D(f1, 0, 1., 0) + >>> + >>> f3(1.5, 1.5, 1.5) + 1.125 + >>> f3(-0.3, -1.3, -0.3) + 0.063 + """ + + def __init__(self, object function3d, double period_x, double period_y, double period_z): + + if not callable(function3d): + raise TypeError("function2d is not callable.") + + self.function3d = autowrap_function3d(function3d) + + if period_x < 0: + raise ValueError("Argument period_x must be >= 0.") + if period_y < 0: + raise ValueError("Argument period_y must be >= 0.") + if period_z < 0: + raise ValueError("Argument period_z must be >= 0.") + + self.period_x = period_x + self.period_y = period_y + self.period_z = period_z + + cdef double evaluate(self, double x, double y, double z) except? -1e999: + """Return the value of periodic function.""" + + x = remainder(x, self.period_x) + y = remainder(y, self.period_y) + z = remainder(z, self.period_z) + + return self.function3d.evaluate(x, y, z) + + +cdef class VectorPeriodicMapper1D(VectorFunction1D): + """ + Maps a periodic 1D vector function into 1D space. + + :param VectorFunction1D function1d: The periodic 1D vector function to map + defined in the [0, period) interval. + :param double period: The period of the function. + + .. code-block:: pycon + + >>> from raysect.core.math import Vector3D + >>> from cherab.core.math import VectorPeriodicMapper1D + >>> + >>> def f1(x): + >>> return Vector3D(x, 0, 0) + >>> + >>> f2 = VectorPeriodicMapper1D(f1, 1.) + >>> + >>> f2(1.5) + Vector3D(0.5, 0, 0) + >>> f2(-0.3) + Vector3D(0.7, 0, 0) + """ + + def __init__(self, object function1d, double period): + + if not callable(function1d): + raise TypeError("function1d is not callable.") + + self.function1d = autowrap_vectorfunction1d(function1d) + + if period <= 0: + raise ValueError("Argument period must be positive.") + + self.period = period + + cdef Vector3D evaluate(self, double x): + """Return the value of periodic function.""" + + return self.function1d.evaluate(remainder(x, self.period)) + + +cdef class VectorPeriodicMapper2D(VectorFunction2D): + """ + Maps a periodic 2D vector function into 2D space. + + Set period_x/period_y to 0 if the function is not periodic along x/y axis. + + :param VectorFunction2D function2d: The periodic 2D vector function to map defined in + the ([0, period_x), [0, period_y)) intervals. + :param double period_x: The period of the function along x-axis. + 0 if not periodic. + :param double period_y: The period of the function along y-axis. + 0 if not periodic. + + .. code-block:: pycon + + >>> from cherab.core.math import VectorPeriodicMapper2D + >>> + >>> def f1(x, y): + >>> return Vector3D(x, y, 0) + >>> + >>> f2 = VectorPeriodicMapper2D(f1, 1., 1.) + >>> + >>> f2(1.5, 1.5) + Vector3D(0.5, 0.5, 0) + >>> f2(-0.3, -1.3) + Vector3D(0.7, 0.7, 0) + >>> + >>> f3 = VectorPeriodicMapper2D(f1, 1., 0) + >>> + >>> f3(1.5, 1.5) + Vector3D(0.5, 1.5, 0) + >>> f3(-0.3, -1.3) + Vector3D(0.7, -1.3, 0) + """ + + def __init__(self, object function2d, double period_x, double period_y): + + if not callable(function2d): + raise TypeError("function2d is not callable.") + + self.function2d = autowrap_vectorfunction2d(function2d) + + if period_x < 0: + raise ValueError("Argument period_x must be >= 0.") + if period_y < 0: + raise ValueError("Argument period_y must be >= 0.") + + self.period_x = period_x + self.period_y = period_y + + cdef Vector3D evaluate(self, double x, double y): + """Return the value of periodic function.""" + + x = remainder(x, self.period_x) + y = remainder(y, self.period_y) + + return self.function2d.evaluate(x, y) + + +cdef class VectorPeriodicMapper3D(VectorFunction3D): + """ + Maps a periodic 3D vector function into 3D space. + + Set period_x/period_y/period_z to 0 if the function is not periodic along x/y/z axis. + + :param VectorFunction3D function3d: The periodic 3D vector function to map defined in the + ([0, period_x), [0, period_y), [0, period_z)) intervals. + :param double period_x: The period of the function along x-axis. + 0 if not periodic. + :param double period_y: The period of the function along y-axis. + 0 if not periodic. + :param double period_z: The period of the function along z-axis. + 0 if not periodic. + + .. code-block:: pycon + + >>> from cherab.core.math import PeriodicMapper3D + >>> + >>> def f1(x, y, z): + >>> return Vector3D(x, y, z) + >>> + >>> f2 = VectorPeriodicMapper3D(f1, 1., 1., 1.) + >>> + >>> f2(1.5, 1.5, 1.5) + Vector3D(0.5, 0.5, 0.5) + >>> f2(-0.3, -1.3, -2.3) + Vector3D(0.7, 0.7, 0.7) + >>> + >>> f3 = VectorPeriodicMapper3D(f1, 0, 1., 0) + >>> + >>> f3(1.5, 0.5, 1.5) + Vector3D(1.5, 0.5, 1.5) + """ + + def __init__(self, object function3d, double period_x, double period_y, double period_z): + + if not callable(function3d): + raise TypeError("function2d is not callable.") + + self.function3d = autowrap_vectorfunction3d(function3d) + + if period_x < 0: + raise ValueError("Argument period_x must be >= 0.") + if period_y < 0: + raise ValueError("Argument period_y must be >= 0.") + if period_z < 0: + raise ValueError("Argument period_z must be >= 0.") + + self.period_x = period_x + self.period_y = period_y + self.period_z = period_z + + cdef Vector3D evaluate(self, double x, double y, double z): + """Return the value of periodic function.""" + + x = remainder(x, self.period_x) + y = remainder(y, self.period_y) + z = remainder(z, self.period_z) + + return self.function3d.evaluate(x, y, z) From 251b11b457cb1daa9cd5c82017a9bfd6c41aa94d Mon Sep 17 00:00:00 2001 From: vsnever Date: Fri, 25 Nov 2022 22:50:35 +0300 Subject: [PATCH 2/9] Add tests for periodic and cylindrical mappers. --- cherab/core/math/tests/test_mappers.py | 114 ++++++++++++++++++++++++- 1 file changed, 111 insertions(+), 3 deletions(-) diff --git a/cherab/core/math/tests/test_mappers.py b/cherab/core/math/tests/test_mappers.py index 777a6a85..69a5f771 100644 --- a/cherab/core/math/tests/test_mappers.py +++ b/cherab/core/math/tests/test_mappers.py @@ -1,6 +1,6 @@ -# Copyright 2016-2018 Euratom -# Copyright 2016-2018 United Kingdom Atomic Energy Authority -# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas # # Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the # European Commission - subsequent versions of the EUPL (the "Licence"); @@ -16,6 +16,7 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. +from raysect.core.math import Vector3D from cherab.core.math import mappers import numpy as np import unittest @@ -142,5 +143,112 @@ def test_axisymmetric_mapper_invalid_arg(self): """An error must be raised if the given argument is not callable.""" self.assertRaises(TypeError, mappers.AxisymmetricMapper, "blah") + def test_cylindrical_mapper(self): + """Cylindrical mapper.""" + def f3d(r, phi, z): return r * np.cos(phi) + z + cyl_func = mappers.CylindricalMapper(f3d) + self.assertAlmostEqual(cyl_func(1., 1., 0.5), + f3d(np.sqrt(2.), 0.25 * np.pi, 0.5), + places=10) + + def test_cylindrical_mapper_invalid_arg(self): + """An error must be raised if the given argument is not callable.""" + self.assertRaises(TypeError, mappers.CylindricalMapper, "blah") + + def test_vector_cylindrical_mapper(self): + """Cylindrical mapper.""" + def f3d(r, phi, z): return Vector3D(np.sin(phi), r * z, np.cos(phi)) + cyl_func = mappers.VectorCylindricalMapper(f3d) + vec1 = cyl_func(1., 1., 1.) + vec2 = Vector3D(-0.5, 1.5, 1 / np.sqrt(2)) + np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) + + def test_vector_cylindrical_mapper_invalid_arg(self): + """An error must be raised if the given argument is not callable.""" + self.assertRaises(TypeError, mappers.VectorCylindricalMapper, "blah") + + def test_periodic_mapper_1d(self): + """1D periodic mapper""" + period_func = mappers.PeriodicMapper1D(self.function1d, np.pi) + self.assertAlmostEqual(period_func(1.4 * np.pi), + self.function1d(0.4 * np.pi), + places=10) + self.assertAlmostEqual(period_func(-0.4 * np.pi), + self.function1d(0.6 * np.pi), + places=10) + + def test_periodic_mapper_1d_invalid_arg(self): + """1D periodic mapper. Invalid arguments.""" + # 1st argument is not callable + self.assertRaises(TypeError, mappers.PeriodicMapper1D, "blah", np.pi) + # period is not a number + self.assertRaises(TypeError, mappers.PeriodicMapper1D, self.function1d, "blah") + # period is negative + self.assertRaises(ValueError, mappers.PeriodicMapper1D, self.function1d, -1) + + def test_periodic_mapper_2d(self): + """2D periodic mapper""" + period_func = mappers.PeriodicMapper2D(self.function2d, 1, np.pi) + self.assertAlmostEqual(period_func(-0.4, 1.4 * np.pi), + self.function2d(0.6, 0.4 * np.pi), + places=10) + # Periodic only along x + period_func = mappers.PeriodicMapper2D(self.function2d, 1., 0) + self.assertAlmostEqual(period_func(-0.4, 1.4 * np.pi), + self.function2d(0.6, 1.4 * np.pi), + places=10) + # Periodic only along y + period_func = mappers.PeriodicMapper2D(self.function2d, 0, np.pi) + self.assertAlmostEqual(period_func(-0.4, 1.4 * np.pi), + self.function2d(-0.4, 0.4 * np.pi), + places=10) + + def test_periodic_mapper_2d_invalid_arg(self): + """2D periodic mapper. Invalid arguments.""" + # 1st argument is not callable + self.assertRaises(TypeError, mappers.PeriodicMapper2D, "blah", np.pi, np.pi) + # period is not a number + self.assertRaises(TypeError, mappers.PeriodicMapper2D, self.function2d, "blah", np.pi) + self.assertRaises(TypeError, mappers.PeriodicMapper2D, self.function2d, np.pi, "blah") + # period is negative + self.assertRaises(ValueError, mappers.PeriodicMapper2D, self.function2d, -1, np.pi) + self.assertRaises(ValueError, mappers.PeriodicMapper2D, self.function2d, np.pi, -1) + + def test_periodic_mapper_3d(self): + """3D periodic mapper""" + period_func = mappers.PeriodicMapper3D(self.function3d, 1, 1, 1) + self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), + self.function3d(0.6, 0.4, 0.1), + places=10) + # Periodic only along y and z + period_func = mappers.PeriodicMapper3D(self.function3d, 0, 1, 1) + self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), + self.function3d(-0.4, 0.4, 0.1), + places=10) + # Periodic only along x and z + period_func = mappers.PeriodicMapper3D(self.function3d, 1, 0, 1) + self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), + self.function3d(0.6, 1.4, 0.1), + places=10) + # Periodic only along x and y + period_func = mappers.PeriodicMapper3D(self.function3d, 1, 1, 0) + self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), + self.function3d(0.6, 0.4, 2.1), + places=10) + + def test_periodic_mapper_3d_invalid_arg(self): + """2D periodic mapper. Invalid arguments.""" + # 1st argument is not callable + self.assertRaises(TypeError, mappers.PeriodicMapper3D, "blah", np.pi, np.pi, np.pi) + # period is not a number + self.assertRaises(TypeError, mappers.PeriodicMapper3D, self.function3d, "blah", np.pi, np.pi) + self.assertRaises(TypeError, mappers.PeriodicMapper3D, self.function3d, np.pi, "blah", np.pi) + self.assertRaises(TypeError, mappers.PeriodicMapper3D, self.function3d, np.pi, np.pi, "blah") + # period is negative + self.assertRaises(ValueError, mappers.PeriodicMapper3D, self.function3d, -1, np.pi, np.pi) + self.assertRaises(ValueError, mappers.PeriodicMapper3D, self.function3d, np.pi, -1, np.pi) + self.assertRaises(ValueError, mappers.PeriodicMapper3D, self.function3d, np.pi, np.pi, -1) + + if __name__ == '__main__': unittest.main() \ No newline at end of file From c0a60eeda1d7c64921767aae9cc4c4ca449a9fcf Mon Sep 17 00:00:00 2001 From: vsnever Date: Fri, 25 Nov 2022 23:14:16 +0300 Subject: [PATCH 3/9] Update CHANGELOG. --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index dbe69a21..b8a3bc5d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,8 @@ New: * Add verbose parameter to SartOpencl solver (default is False). (#358) * Add Generomak core plasma profiles. (#360) * Add toroidal_mesh_from_polygon for making mesh for not fully-360 degrees axisymmetric elements. (#365) +* Add PeriodicMapperXD and VectorPeriodicMapperXD to support the data simulated with periodic boundary conditions. (#387) +* Add CylindricalMapper and VectorCylindricalMapper to map functions from cylindrical to Cartesian coordinates. (#387) Bug Fixes: ---------- From 73b87df644b09e8b536e14f60b95ab5976ba9f2f Mon Sep 17 00:00:00 2001 From: vsnever Date: Mon, 26 Dec 2022 14:06:40 +0300 Subject: [PATCH 4/9] Rename periodic/cylindrical mappers to periodic/cylindrical transforms and move them to cherab.core.math.transform module. --- cherab/core/math/__init__.pxd | 1 + cherab/core/math/__init__.py | 6 +- cherab/core/math/mappers.pxd | 62 --- cherab/core/math/mappers.pyx | 436 +-------------------- cherab/core/math/tests/test_mappers.py | 114 +----- cherab/core/math/tests/test_transform.py | 157 ++++++++ cherab/core/math/transform/__init__.pxd | 20 + cherab/core/math/transform/__init__.py | 21 + cherab/core/math/transform/cylindrical.pxd | 30 ++ cherab/core/math/transform/cylindrical.pyx | 128 ++++++ cherab/core/math/transform/periodic.pxd | 72 ++++ cherab/core/math/transform/periodic.pyx | 352 +++++++++++++++++ 12 files changed, 796 insertions(+), 603 deletions(-) create mode 100644 cherab/core/math/tests/test_transform.py create mode 100644 cherab/core/math/transform/__init__.pxd create mode 100644 cherab/core/math/transform/__init__.py create mode 100644 cherab/core/math/transform/cylindrical.pxd create mode 100644 cherab/core/math/transform/cylindrical.pyx create mode 100644 cherab/core/math/transform/periodic.pxd create mode 100644 cherab/core/math/transform/periodic.pyx diff --git a/cherab/core/math/__init__.pxd b/cherab/core/math/__init__.pxd index 99e92fbe..710ee381 100644 --- a/cherab/core/math/__init__.pxd +++ b/cherab/core/math/__init__.pxd @@ -25,3 +25,4 @@ from cherab.core.math.clamp cimport * from cherab.core.math.mappers cimport * from cherab.core.math.mask cimport * from cherab.core.math.slice cimport * +from cherab.core.math.transform cimport * diff --git a/cherab/core/math/__init__.py b/cherab/core/math/__init__.py index 4f821332..85336afa 100644 --- a/cherab/core/math/__init__.py +++ b/cherab/core/math/__init__.py @@ -34,8 +34,8 @@ from .mappers import IsoMapper2D, IsoMapper3D from .mappers import Swizzle2D, Swizzle3D from .mappers import AxisymmetricMapper, VectorAxisymmetricMapper -from .mappers import CylindricalMapper, VectorCylindricalMapper -from .mappers import PeriodicMapper1D, PeriodicMapper2D, PeriodicMapper3D -from .mappers import VectorPeriodicMapper1D, VectorPeriodicMapper2D, VectorPeriodicMapper3D from .mask import PolygonMask2D from .slice import Slice2D, Slice3D +from .transform import CylindricalTransform, VectorCylindricalTransform +from .transform import PeriodicTransform1D, PeriodicTransform2D, PeriodicTransform3D +from .transform import VectorPeriodicTransform1D, VectorPeriodicTransform2D, VectorPeriodicTransform3D diff --git a/cherab/core/math/mappers.pxd b/cherab/core/math/mappers.pxd index 4fee0cbf..3835e949 100644 --- a/cherab/core/math/mappers.pxd +++ b/cherab/core/math/mappers.pxd @@ -16,12 +16,9 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. -from libc.math cimport fmod from raysect.core.math.function.float cimport Function1D, Function2D, Function3D -from raysect.core.math.function.vector3d cimport Function1D as VectorFunction1D from raysect.core.math.function.vector3d cimport Function2D as VectorFunction2D from raysect.core.math.function.vector3d cimport Function3D as VectorFunction3D -from raysect.core cimport Vector3D cdef class IsoMapper2D(Function2D): @@ -58,62 +55,3 @@ cdef class AxisymmetricMapper(Function3D): cdef class VectorAxisymmetricMapper(VectorFunction3D): cdef readonly VectorFunction2D function2d - - -cdef class CylindricalMapper(Function3D): - - cdef readonly Function3D function3d - - -cdef class VectorCylindricalMapper(VectorFunction3D): - - cdef readonly VectorFunction3D function3d - - -cdef inline double remainder(double x1, double x2) nogil: - if x2 == 0: - return x1 - x1 = fmod(x1, x2) - return x1 + x2 if (x1 < 0) else x1 - - -cdef class PeriodicMapper1D(Function1D): - - cdef: - readonly Function1D function1d - readonly double period - - -cdef class PeriodicMapper2D(Function2D): - - cdef: - readonly Function2D function2d - double period_x, period_y - - -cdef class PeriodicMapper3D(Function3D): - - cdef: - readonly Function3D function3d - readonly double period_x, period_y, period_z - - -cdef class VectorPeriodicMapper1D(VectorFunction1D): - - cdef: - readonly VectorFunction1D function1d - readonly double period - - -cdef class VectorPeriodicMapper2D(VectorFunction2D): - - cdef: - readonly VectorFunction2D function2d - readonly double period_x, period_y - - -cdef class VectorPeriodicMapper3D(VectorFunction3D): - - cdef: - readonly VectorFunction3D function3d - readonly double period_x, period_y, period_z diff --git a/cherab/core/math/mappers.pyx b/cherab/core/math/mappers.pyx index 63f868a9..f1326f12 100644 --- a/cherab/core/math/mappers.pyx +++ b/cherab/core/math/mappers.pyx @@ -20,10 +20,9 @@ from libc.math cimport sqrt, atan2, M_PI +from raysect.core.math cimport Vector3D from raysect.core.math.function.float cimport autowrap_function1d, autowrap_function2d, autowrap_function3d -from raysect.core.math.function.vector3d cimport autowrap_function1d as autowrap_vectorfunction1d from raysect.core.math.function.vector3d cimport autowrap_function2d as autowrap_vectorfunction2d -from raysect.core.math.function.vector3d cimport autowrap_function3d as autowrap_vectorfunction3d from raysect.core cimport rotate_z cimport cython @@ -311,436 +310,3 @@ cdef class VectorAxisymmetricMapper(VectorFunction3D): # perform axisymmetric rotation return self.function2d.evaluate(r, z).transform(rotate_z(phi)) - - -cdef class CylindricalMapper(Function3D): - """ - Converts Cartesian coordinates to cylindrical coordinates and calls a 3D function - defined in cylindrical coordinates, f(r, :math:`\\phi`, z). - - The angular coordinate is given in radians. - - Positive angular coordinate is measured counterclockwise from the xz plane. - - :param Function3D function3d: The function to be mapped. Must be defined - in the interval (:math:`-\\pi`, :math:`\\pi`] - on the angular axis. - - .. code-block:: pycon - - >>> from math import sqrt, cos - >>> from cherab.core.math import CylindricalMapper - >>> - >>> def my_func(r, phi, z): - >>> return r * cos(phi) - >>> - >>> f = CylindricalMapper(my_func) - >>> - >>> f(1, 0, 0) - 1.0 - >>> f(0.5 * sqrt(3), 0.5, 0) - 0.8660254037844385 - """ - - def __init__(self, object function3d): - - if not callable(function3d): - raise TypeError("Function3D is not callable.") - - self.function3d = autowrap_function3d(function3d) - - cdef double evaluate(self, double x, double y, double z) except? -1e999: - """ - Converts to cylindrical coordinates and evaluates the function - defined in cylindrical coordinates. - """ - cdef double r, phi - - r = sqrt(x * x + y * y) - phi = atan2(y, x) - - return self.function3d.evaluate(r, phi, z) - - -cdef class VectorCylindricalMapper(VectorFunction3D): - """ - Converts Cartesian coordinates to cylindrical coordinates, calls - a 3D vector function defined in cylindrical coordinates, f(r, :math:`\\phi`, z), - then converts the returned 3D vector to Cartesian coordinates. - - The angular coordinate is given in radians. - - Positive angular coordinate is measured counterclockwise from the xz plane. - - :param VectorFunction3D function3d: The function to be mapped. Must be defined - in the interval (:math:`-\\pi`, :math:`\\pi`] - on the angular axis. - - .. code-block:: pycon - - >>> from math import sqrt, cos - >>> from raysect.core.math import Vector3D - >>> from cherab.core.math import VectorCylindricalMapper - >>> - >>> def my_vec_func(r, phi, z): - >>> v = Vector3D(0, 1, 0) - >>> v.length = r * abs(cos(phi)) - >>> return v - >>> - >>> f = VectorCylindricalMapper(my_vec_func) - >>> - >>> f(1, 0, 0) - Vector3D(0.0, 1.0, 0.0) - >>> f(1/sqrt(2), 1/sqrt(2), 0) - Vector3D(-0.5, 0.5, 0.0) - """ - - def __init__(self, object function3d): - - if not callable(function3d): - raise TypeError("Function3D is not callable.") - - self.function3d = autowrap_vectorfunction3d(function3d) - - @cython.cdivision(True) - cdef Vector3D evaluate(self, double x, double y, double z): - """ - Converts to cylindrical coordinates, evaluates the vector function - defined in cylindrical coordinates and rotates the resulting vector - around z-axis. - """ - cdef double r, phi - - r = sqrt(x * x + y * y) - phi = atan2(y, x) - - return self.function3d.evaluate(r, phi, z).transform(rotate_z(phi / M_PI * 180)) - - -cdef class PeriodicMapper1D(Function1D): - """ - Maps a periodic 1D function into 1D space. - - :param Function1D function1d: The periodic 1D function to map defined - in the [0, period) interval. - :param double period: The period of the function. - - .. code-block:: pycon - - >>> from cherab.core.math import PeriodicMapper1D - >>> - >>> def f1(x): - >>> return x - >>> - >>> f2 = PeriodicMapper1D(f1, 1.) - >>> - >>> f2(1.5) - 0.5 - >>> f2(-0.3) - 0.7 - """ - - def __init__(self, object function1d, double period): - - if not callable(function1d): - raise TypeError("function1d is not callable.") - - self.function1d = autowrap_function1d(function1d) - - if period <= 0: - raise ValueError("Argument period must be positive.") - - self.period = period - - cdef double evaluate(self, double x) except? -1e999: - """Return the value of periodic function.""" - - return self.function1d.evaluate(remainder(x, self.period)) - - -cdef class PeriodicMapper2D(Function2D): - """ - Maps a periodic 2D function into 2D space. - - Set period_x/period_y to 0 if the function is not periodic along x/y axis. - - :param Function2D function2d: The periodic 2D function to map defined - in the ([0, period_x), [0, period_y)) intervals. - :param double period_x: The period of the function along x-axis. - 0 if not periodic. - :param double period_y: The period of the function along y-axis. - 0 if not periodic. - - .. code-block:: pycon - - >>> from cherab.core.math import PeriodicMapper2D - >>> - >>> def f1(x, y): - >>> return x * y - >>> - >>> f2 = PeriodicMapper2D(f1, 1., 1.) - >>> - >>> f2(1.5, 1.5) - 0.25 - >>> f2(-0.3, -1.3) - 0.49 - >>> - >>> f3 = PeriodicMapper2D(f1, 1., 0) - >>> - >>> f3(1.5, 1.5) - 0.75 - >>> f3(-0.3, -1.3) - -0.91 - """ - - def __init__(self, object function2d, double period_x, double period_y): - - if not callable(function2d): - raise TypeError("function2d is not callable.") - - self.function2d = autowrap_function2d(function2d) - - if period_x < 0: - raise ValueError("Argument period_x must be >= 0.") - if period_y < 0: - raise ValueError("Argument period_y must be >= 0.") - - self.period_x = period_x - self.period_y = period_y - - cdef double evaluate(self, double x, double y) except? -1e999: - """Return the value of periodic function.""" - - x = remainder(x, self.period_x) - y = remainder(y, self.period_y) - - return self.function2d.evaluate(x, y) - - -cdef class PeriodicMapper3D(Function3D): - """ - Maps a periodic 3D function into 3D space. - - Set period_x/period_y/period_z to 0 if the function is not periodic along x/y/z axis. - - :param Function3D function3d: The periodic 3D function to map defined in the - ([0, period_x), [0, period_y), [0, period_z)) intervals. - :param double period_x: The period of the function along x-axis. - 0 if not periodic. - :param double period_y: The period of the function along y-axis. - 0 if not periodic. - :param double period_z: The period of the function along z-axis. - 0 if not periodic. - - .. code-block:: pycon - - >>> from cherab.core.math import PeriodicMapper3D - >>> - >>> def f1(x, y, z): - >>> return x * y * z - >>> - >>> f2 = PeriodicMapper3D(f1, 1., 1., 1.) - >>> - >>> f2(1.5, 1.5, 1.5) - 0.125 - >>> f2(-0.3, -1.3, -2.3) - 0.343 - >>> - >>> f3 = PeriodicMapper3D(f1, 0, 1., 0) - >>> - >>> f3(1.5, 1.5, 1.5) - 1.125 - >>> f3(-0.3, -1.3, -0.3) - 0.063 - """ - - def __init__(self, object function3d, double period_x, double period_y, double period_z): - - if not callable(function3d): - raise TypeError("function2d is not callable.") - - self.function3d = autowrap_function3d(function3d) - - if period_x < 0: - raise ValueError("Argument period_x must be >= 0.") - if period_y < 0: - raise ValueError("Argument period_y must be >= 0.") - if period_z < 0: - raise ValueError("Argument period_z must be >= 0.") - - self.period_x = period_x - self.period_y = period_y - self.period_z = period_z - - cdef double evaluate(self, double x, double y, double z) except? -1e999: - """Return the value of periodic function.""" - - x = remainder(x, self.period_x) - y = remainder(y, self.period_y) - z = remainder(z, self.period_z) - - return self.function3d.evaluate(x, y, z) - - -cdef class VectorPeriodicMapper1D(VectorFunction1D): - """ - Maps a periodic 1D vector function into 1D space. - - :param VectorFunction1D function1d: The periodic 1D vector function to map - defined in the [0, period) interval. - :param double period: The period of the function. - - .. code-block:: pycon - - >>> from raysect.core.math import Vector3D - >>> from cherab.core.math import VectorPeriodicMapper1D - >>> - >>> def f1(x): - >>> return Vector3D(x, 0, 0) - >>> - >>> f2 = VectorPeriodicMapper1D(f1, 1.) - >>> - >>> f2(1.5) - Vector3D(0.5, 0, 0) - >>> f2(-0.3) - Vector3D(0.7, 0, 0) - """ - - def __init__(self, object function1d, double period): - - if not callable(function1d): - raise TypeError("function1d is not callable.") - - self.function1d = autowrap_vectorfunction1d(function1d) - - if period <= 0: - raise ValueError("Argument period must be positive.") - - self.period = period - - cdef Vector3D evaluate(self, double x): - """Return the value of periodic function.""" - - return self.function1d.evaluate(remainder(x, self.period)) - - -cdef class VectorPeriodicMapper2D(VectorFunction2D): - """ - Maps a periodic 2D vector function into 2D space. - - Set period_x/period_y to 0 if the function is not periodic along x/y axis. - - :param VectorFunction2D function2d: The periodic 2D vector function to map defined in - the ([0, period_x), [0, period_y)) intervals. - :param double period_x: The period of the function along x-axis. - 0 if not periodic. - :param double period_y: The period of the function along y-axis. - 0 if not periodic. - - .. code-block:: pycon - - >>> from cherab.core.math import VectorPeriodicMapper2D - >>> - >>> def f1(x, y): - >>> return Vector3D(x, y, 0) - >>> - >>> f2 = VectorPeriodicMapper2D(f1, 1., 1.) - >>> - >>> f2(1.5, 1.5) - Vector3D(0.5, 0.5, 0) - >>> f2(-0.3, -1.3) - Vector3D(0.7, 0.7, 0) - >>> - >>> f3 = VectorPeriodicMapper2D(f1, 1., 0) - >>> - >>> f3(1.5, 1.5) - Vector3D(0.5, 1.5, 0) - >>> f3(-0.3, -1.3) - Vector3D(0.7, -1.3, 0) - """ - - def __init__(self, object function2d, double period_x, double period_y): - - if not callable(function2d): - raise TypeError("function2d is not callable.") - - self.function2d = autowrap_vectorfunction2d(function2d) - - if period_x < 0: - raise ValueError("Argument period_x must be >= 0.") - if period_y < 0: - raise ValueError("Argument period_y must be >= 0.") - - self.period_x = period_x - self.period_y = period_y - - cdef Vector3D evaluate(self, double x, double y): - """Return the value of periodic function.""" - - x = remainder(x, self.period_x) - y = remainder(y, self.period_y) - - return self.function2d.evaluate(x, y) - - -cdef class VectorPeriodicMapper3D(VectorFunction3D): - """ - Maps a periodic 3D vector function into 3D space. - - Set period_x/period_y/period_z to 0 if the function is not periodic along x/y/z axis. - - :param VectorFunction3D function3d: The periodic 3D vector function to map defined in the - ([0, period_x), [0, period_y), [0, period_z)) intervals. - :param double period_x: The period of the function along x-axis. - 0 if not periodic. - :param double period_y: The period of the function along y-axis. - 0 if not periodic. - :param double period_z: The period of the function along z-axis. - 0 if not periodic. - - .. code-block:: pycon - - >>> from cherab.core.math import PeriodicMapper3D - >>> - >>> def f1(x, y, z): - >>> return Vector3D(x, y, z) - >>> - >>> f2 = VectorPeriodicMapper3D(f1, 1., 1., 1.) - >>> - >>> f2(1.5, 1.5, 1.5) - Vector3D(0.5, 0.5, 0.5) - >>> f2(-0.3, -1.3, -2.3) - Vector3D(0.7, 0.7, 0.7) - >>> - >>> f3 = VectorPeriodicMapper3D(f1, 0, 1., 0) - >>> - >>> f3(1.5, 0.5, 1.5) - Vector3D(1.5, 0.5, 1.5) - """ - - def __init__(self, object function3d, double period_x, double period_y, double period_z): - - if not callable(function3d): - raise TypeError("function2d is not callable.") - - self.function3d = autowrap_vectorfunction3d(function3d) - - if period_x < 0: - raise ValueError("Argument period_x must be >= 0.") - if period_y < 0: - raise ValueError("Argument period_y must be >= 0.") - if period_z < 0: - raise ValueError("Argument period_z must be >= 0.") - - self.period_x = period_x - self.period_y = period_y - self.period_z = period_z - - cdef Vector3D evaluate(self, double x, double y, double z): - """Return the value of periodic function.""" - - x = remainder(x, self.period_x) - y = remainder(y, self.period_y) - z = remainder(z, self.period_z) - - return self.function3d.evaluate(x, y, z) diff --git a/cherab/core/math/tests/test_mappers.py b/cherab/core/math/tests/test_mappers.py index 69a5f771..a1d0eb73 100644 --- a/cherab/core/math/tests/test_mappers.py +++ b/cherab/core/math/tests/test_mappers.py @@ -36,6 +36,9 @@ def f2d(x, y): return x*np.sin(y) def f3d(x, y, z): return x*x*np.exp(y)-2*z*y self.function3d = f3d + def vecf2d(r, z): return Vector3D(0, r, z) + self.vectorfunction2d = vecf2d + def test_iso_mapper_2d(self): """Composition of a 1D and a 2D function.""" @@ -143,112 +146,17 @@ def test_axisymmetric_mapper_invalid_arg(self): """An error must be raised if the given argument is not callable.""" self.assertRaises(TypeError, mappers.AxisymmetricMapper, "blah") - def test_cylindrical_mapper(self): - """Cylindrical mapper.""" - def f3d(r, phi, z): return r * np.cos(phi) + z - cyl_func = mappers.CylindricalMapper(f3d) - self.assertAlmostEqual(cyl_func(1., 1., 0.5), - f3d(np.sqrt(2.), 0.25 * np.pi, 0.5), - places=10) - - def test_cylindrical_mapper_invalid_arg(self): - """An error must be raised if the given argument is not callable.""" - self.assertRaises(TypeError, mappers.CylindricalMapper, "blah") - - def test_vector_cylindrical_mapper(self): - """Cylindrical mapper.""" - def f3d(r, phi, z): return Vector3D(np.sin(phi), r * z, np.cos(phi)) - cyl_func = mappers.VectorCylindricalMapper(f3d) - vec1 = cyl_func(1., 1., 1.) - vec2 = Vector3D(-0.5, 1.5, 1 / np.sqrt(2)) + def test_vector_axisymmetric_mapper(self): + """Vector axisymmetric mapper.""" + symm_func = mappers.VectorAxisymmetricMapper(self.vectorfunction2d) + vec1 = symm_func(1., 1., 1.) + vec2 = Vector3D(-1., 1., 1.) np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) - def test_vector_cylindrical_mapper_invalid_arg(self): + def test_vector_axisymmetric_mapper_invalid_arg(self): """An error must be raised if the given argument is not callable.""" - self.assertRaises(TypeError, mappers.VectorCylindricalMapper, "blah") - - def test_periodic_mapper_1d(self): - """1D periodic mapper""" - period_func = mappers.PeriodicMapper1D(self.function1d, np.pi) - self.assertAlmostEqual(period_func(1.4 * np.pi), - self.function1d(0.4 * np.pi), - places=10) - self.assertAlmostEqual(period_func(-0.4 * np.pi), - self.function1d(0.6 * np.pi), - places=10) - - def test_periodic_mapper_1d_invalid_arg(self): - """1D periodic mapper. Invalid arguments.""" - # 1st argument is not callable - self.assertRaises(TypeError, mappers.PeriodicMapper1D, "blah", np.pi) - # period is not a number - self.assertRaises(TypeError, mappers.PeriodicMapper1D, self.function1d, "blah") - # period is negative - self.assertRaises(ValueError, mappers.PeriodicMapper1D, self.function1d, -1) - - def test_periodic_mapper_2d(self): - """2D periodic mapper""" - period_func = mappers.PeriodicMapper2D(self.function2d, 1, np.pi) - self.assertAlmostEqual(period_func(-0.4, 1.4 * np.pi), - self.function2d(0.6, 0.4 * np.pi), - places=10) - # Periodic only along x - period_func = mappers.PeriodicMapper2D(self.function2d, 1., 0) - self.assertAlmostEqual(period_func(-0.4, 1.4 * np.pi), - self.function2d(0.6, 1.4 * np.pi), - places=10) - # Periodic only along y - period_func = mappers.PeriodicMapper2D(self.function2d, 0, np.pi) - self.assertAlmostEqual(period_func(-0.4, 1.4 * np.pi), - self.function2d(-0.4, 0.4 * np.pi), - places=10) - - def test_periodic_mapper_2d_invalid_arg(self): - """2D periodic mapper. Invalid arguments.""" - # 1st argument is not callable - self.assertRaises(TypeError, mappers.PeriodicMapper2D, "blah", np.pi, np.pi) - # period is not a number - self.assertRaises(TypeError, mappers.PeriodicMapper2D, self.function2d, "blah", np.pi) - self.assertRaises(TypeError, mappers.PeriodicMapper2D, self.function2d, np.pi, "blah") - # period is negative - self.assertRaises(ValueError, mappers.PeriodicMapper2D, self.function2d, -1, np.pi) - self.assertRaises(ValueError, mappers.PeriodicMapper2D, self.function2d, np.pi, -1) - - def test_periodic_mapper_3d(self): - """3D periodic mapper""" - period_func = mappers.PeriodicMapper3D(self.function3d, 1, 1, 1) - self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), - self.function3d(0.6, 0.4, 0.1), - places=10) - # Periodic only along y and z - period_func = mappers.PeriodicMapper3D(self.function3d, 0, 1, 1) - self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), - self.function3d(-0.4, 0.4, 0.1), - places=10) - # Periodic only along x and z - period_func = mappers.PeriodicMapper3D(self.function3d, 1, 0, 1) - self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), - self.function3d(0.6, 1.4, 0.1), - places=10) - # Periodic only along x and y - period_func = mappers.PeriodicMapper3D(self.function3d, 1, 1, 0) - self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), - self.function3d(0.6, 0.4, 2.1), - places=10) - - def test_periodic_mapper_3d_invalid_arg(self): - """2D periodic mapper. Invalid arguments.""" - # 1st argument is not callable - self.assertRaises(TypeError, mappers.PeriodicMapper3D, "blah", np.pi, np.pi, np.pi) - # period is not a number - self.assertRaises(TypeError, mappers.PeriodicMapper3D, self.function3d, "blah", np.pi, np.pi) - self.assertRaises(TypeError, mappers.PeriodicMapper3D, self.function3d, np.pi, "blah", np.pi) - self.assertRaises(TypeError, mappers.PeriodicMapper3D, self.function3d, np.pi, np.pi, "blah") - # period is negative - self.assertRaises(ValueError, mappers.PeriodicMapper3D, self.function3d, -1, np.pi, np.pi) - self.assertRaises(ValueError, mappers.PeriodicMapper3D, self.function3d, np.pi, -1, np.pi) - self.assertRaises(ValueError, mappers.PeriodicMapper3D, self.function3d, np.pi, np.pi, -1) + self.assertRaises(TypeError, mappers.VectorAxisymmetricMapper, "blah") if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main() diff --git a/cherab/core/math/tests/test_transform.py b/cherab/core/math/tests/test_transform.py new file mode 100644 index 00000000..2ee4eefc --- /dev/null +++ b/cherab/core/math/tests/test_transform.py @@ -0,0 +1,157 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from raysect.core.math import Vector3D +from cherab.core.math import transform +import numpy as np +import unittest + + +class TestCylindricalTransform(unittest.TestCase): + """Cylindrical transform tests.""" + + def setUp(self): + """Initialisation with functions to map.""" + + def f3d(r, phi, z): return r * np.cos(phi) + z + self.function3d = f3d + + def vecf3d(r, phi, z): return Vector3D(np.sin(phi), r * z, np.cos(phi)) + self.vectorfunction3d = vecf3d + + def test_cylindrical_transform(self): + """Cylindrical transform.""" + cyl_func = transform.CylindricalTransform(self.function3d) + self.assertAlmostEqual(cyl_func(1., 1., 0.5), + self.function3d(np.sqrt(2.), 0.25 * np.pi, 0.5), + places=10) + + def test_cylindrical_transform_invalid_arg(self): + """An error must be raised if the given argument is not callable.""" + self.assertRaises(TypeError, transform.CylindricalTransform, "blah") + + def test_vector_cylindrical_transform(self): + """Cylindrical transform.""" + cyl_func = transform.VectorCylindricalTransform(self.vectorfunction3d) + vec1 = cyl_func(1., 1., 1.) + vec2 = Vector3D(-0.5, 1.5, 1 / np.sqrt(2)) + np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) + + def test_vector_cylindrical_transform_invalid_arg(self): + """An error must be raised if the given argument is not callable.""" + self.assertRaises(TypeError, transform.VectorCylindricalTransform, "blah") + + +class TestPeriodicTransform(unittest.TestCase): + """Periodic transform tests.""" + + def setUp(self): + """Initialisation with functions to map.""" + + def f1d(x): return x*np.cos(x-3) + self.function1d = f1d + def f2d(x, y): return x*np.sin(y) + self.function2d = f2d + def f3d(x, y, z): return x*x*np.exp(y)-2*z*y + self.function3d = f3d + + def test_periodic_transform_1d(self): + """1D periodic transform""" + period_func = transform.PeriodicTransform1D(self.function1d, np.pi) + self.assertAlmostEqual(period_func(1.4 * np.pi), + self.function1d(0.4 * np.pi), + places=10) + self.assertAlmostEqual(period_func(-0.4 * np.pi), + self.function1d(0.6 * np.pi), + places=10) + + def test_periodic_transform_1d_invalid_arg(self): + """1D periodic transform. Invalid arguments.""" + # 1st argument is not callable + self.assertRaises(TypeError, transform.PeriodicTransform1D, "blah", np.pi) + # period is not a number + self.assertRaises(TypeError, transform.PeriodicTransform1D, self.function1d, "blah") + # period is negative + self.assertRaises(ValueError, transform.PeriodicTransform1D, self.function1d, -1) + + def test_periodic_transform_2d(self): + """2D periodic transform""" + period_func = transform.PeriodicTransform2D(self.function2d, 1, np.pi) + self.assertAlmostEqual(period_func(-0.4, 1.4 * np.pi), + self.function2d(0.6, 0.4 * np.pi), + places=10) + # Periodic only along x + period_func = transform.PeriodicTransform2D(self.function2d, 1., 0) + self.assertAlmostEqual(period_func(-0.4, 1.4 * np.pi), + self.function2d(0.6, 1.4 * np.pi), + places=10) + # Periodic only along y + period_func = transform.PeriodicTransform2D(self.function2d, 0, np.pi) + self.assertAlmostEqual(period_func(-0.4, 1.4 * np.pi), + self.function2d(-0.4, 0.4 * np.pi), + places=10) + + def test_periodic_transform_2d_invalid_arg(self): + """2D periodic transform. Invalid arguments.""" + # 1st argument is not callable + self.assertRaises(TypeError, transform.PeriodicTransform2D, "blah", np.pi, np.pi) + # period is not a number + self.assertRaises(TypeError, transform.PeriodicTransform2D, self.function2d, "blah", np.pi) + self.assertRaises(TypeError, transform.PeriodicTransform2D, self.function2d, np.pi, "blah") + # period is negative + self.assertRaises(ValueError, transform.PeriodicTransform2D, self.function2d, -1, np.pi) + self.assertRaises(ValueError, transform.PeriodicTransform2D, self.function2d, np.pi, -1) + + def test_periodic_transform_3d(self): + """3D periodic transform""" + period_func = transform.PeriodicTransform3D(self.function3d, 1, 1, 1) + self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), + self.function3d(0.6, 0.4, 0.1), + places=10) + # Periodic only along y and z + period_func = transform.PeriodicTransform3D(self.function3d, 0, 1, 1) + self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), + self.function3d(-0.4, 0.4, 0.1), + places=10) + # Periodic only along x and z + period_func = transform.PeriodicTransform3D(self.function3d, 1, 0, 1) + self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), + self.function3d(0.6, 1.4, 0.1), + places=10) + # Periodic only along x and y + period_func = transform.PeriodicTransform3D(self.function3d, 1, 1, 0) + self.assertAlmostEqual(period_func(-0.4, 1.4, 2.1), + self.function3d(0.6, 0.4, 2.1), + places=10) + + def test_periodic_transform_3d_invalid_arg(self): + """2D periodic transform. Invalid arguments.""" + # 1st argument is not callable + self.assertRaises(TypeError, transform.PeriodicTransform3D, "blah", np.pi, np.pi, np.pi) + # period is not a number + self.assertRaises(TypeError, transform.PeriodicTransform3D, self.function3d, "blah", np.pi, np.pi) + self.assertRaises(TypeError, transform.PeriodicTransform3D, self.function3d, np.pi, "blah", np.pi) + self.assertRaises(TypeError, transform.PeriodicTransform3D, self.function3d, np.pi, np.pi, "blah") + # period is negative + self.assertRaises(ValueError, transform.PeriodicTransform3D, self.function3d, -1, np.pi, np.pi) + self.assertRaises(ValueError, transform.PeriodicTransform3D, self.function3d, np.pi, -1, np.pi) + self.assertRaises(ValueError, transform.PeriodicTransform3D, self.function3d, np.pi, np.pi, -1) + + +if __name__ == '__main__': + unittest.main() diff --git a/cherab/core/math/transform/__init__.pxd b/cherab/core/math/transform/__init__.pxd new file mode 100644 index 00000000..3ce39b6b --- /dev/null +++ b/cherab/core/math/transform/__init__.pxd @@ -0,0 +1,20 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from cherab.core.math.transform.periodic cimport * +from cherab.core.math.transform.cylindrical cimport * \ No newline at end of file diff --git a/cherab/core/math/transform/__init__.py b/cherab/core/math/transform/__init__.py new file mode 100644 index 00000000..168bb305 --- /dev/null +++ b/cherab/core/math/transform/__init__.py @@ -0,0 +1,21 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from .cylindrical import CylindricalTransform, VectorCylindricalTransform +from .periodic import PeriodicTransform1D, PeriodicTransform2D, PeriodicTransform3D +from .periodic import VectorPeriodicTransform1D, VectorPeriodicTransform2D, VectorPeriodicTransform3D diff --git a/cherab/core/math/transform/cylindrical.pxd b/cherab/core/math/transform/cylindrical.pxd new file mode 100644 index 00000000..504b2fd5 --- /dev/null +++ b/cherab/core/math/transform/cylindrical.pxd @@ -0,0 +1,30 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from raysect.core.math.function.float cimport Function3D +from raysect.core.math.function.vector3d cimport Function3D as VectorFunction3D + + +cdef class CylindricalTransform(Function3D): + + cdef readonly Function3D function3d + + +cdef class VectorCylindricalTransform(VectorFunction3D): + + cdef readonly VectorFunction3D function3d diff --git a/cherab/core/math/transform/cylindrical.pyx b/cherab/core/math/transform/cylindrical.pyx new file mode 100644 index 00000000..cd716cad --- /dev/null +++ b/cherab/core/math/transform/cylindrical.pyx @@ -0,0 +1,128 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from libc.math cimport sqrt, atan2, M_PI + +from raysect.core.math cimport Vector3D +from raysect.core.math.function.float cimport autowrap_function3d +from raysect.core.math.function.vector3d cimport autowrap_function3d as autowrap_vectorfunction3d +from raysect.core cimport rotate_z +cimport cython + +cdef class CylindricalTransform(Function3D): + """ + Converts Cartesian coordinates to cylindrical coordinates and calls a 3D function + defined in cylindrical coordinates, f(r, :math:`\\phi`, z). + + The angular coordinate is given in radians. + + Positive angular coordinate is measured counterclockwise from the xz plane. + + :param Function3D function3d: The function to be mapped. Must be defined + in the interval (:math:`-\\pi`, :math:`\\pi`] + on the angular axis. + + .. code-block:: pycon + + >>> from math import sqrt, cos + >>> from cherab.core.math import CylindricalTransform + >>> + >>> def my_func(r, phi, z): + >>> return r * cos(phi) + >>> + >>> f = CylindricalTransform(my_func) + >>> + >>> f(1, 0, 0) + 1.0 + >>> f(0.5 * sqrt(3), 0.5, 0) + 0.8660254037844385 + """ + + def __init__(self, object function3d): + + if not callable(function3d): + raise TypeError("Function3D is not callable.") + + self.function3d = autowrap_function3d(function3d) + + cdef double evaluate(self, double x, double y, double z) except? -1e999: + """ + Converts to cylindrical coordinates and evaluates the function + defined in cylindrical coordinates. + """ + cdef double r, phi + + r = sqrt(x * x + y * y) + phi = atan2(y, x) + + return self.function3d.evaluate(r, phi, z) + + +cdef class VectorCylindricalTransform(VectorFunction3D): + """ + Converts Cartesian coordinates to cylindrical coordinates, calls + a 3D vector function defined in cylindrical coordinates, f(r, :math:`\\phi`, z), + then converts the returned 3D vector to Cartesian coordinates. + + The angular coordinate is given in radians. + + Positive angular coordinate is measured counterclockwise from the xz plane. + + :param VectorFunction3D function3d: The function to be mapped. Must be defined + in the interval (:math:`-\\pi`, :math:`\\pi`] + on the angular axis. + + .. code-block:: pycon + + >>> from math import sqrt, cos + >>> from raysect.core.math import Vector3D + >>> from cherab.core.math import VectorCylindricalTransform + >>> + >>> def my_vec_func(r, phi, z): + >>> v = Vector3D(0, 1, 0) + >>> v.length = r * abs(cos(phi)) + >>> return v + >>> + >>> f = VectorCylindricalTransform(my_vec_func) + >>> + >>> f(1, 0, 0) + Vector3D(0.0, 1.0, 0.0) + >>> f(1/sqrt(2), 1/sqrt(2), 0) + Vector3D(-0.5, 0.5, 0.0) + """ + + def __init__(self, object function3d): + + if not callable(function3d): + raise TypeError("Function3D is not callable.") + + self.function3d = autowrap_vectorfunction3d(function3d) + + @cython.cdivision(True) + cdef Vector3D evaluate(self, double x, double y, double z): + """ + Converts to cylindrical coordinates, evaluates the vector function + defined in cylindrical coordinates and rotates the resulting vector + around z-axis. + """ + cdef double r, phi + + r = sqrt(x * x + y * y) + phi = atan2(y, x) + + return self.function3d.evaluate(r, phi, z).transform(rotate_z(phi / M_PI * 180)) diff --git a/cherab/core/math/transform/periodic.pxd b/cherab/core/math/transform/periodic.pxd new file mode 100644 index 00000000..e97e4602 --- /dev/null +++ b/cherab/core/math/transform/periodic.pxd @@ -0,0 +1,72 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from libc.math cimport fmod +from raysect.core.math.function.float cimport Function1D, Function2D, Function3D +from raysect.core.math.function.vector3d cimport Function1D as VectorFunction1D +from raysect.core.math.function.vector3d cimport Function2D as VectorFunction2D +from raysect.core.math.function.vector3d cimport Function3D as VectorFunction3D + + +cdef inline double remainder(double x1, double x2) nogil: + if x2 == 0: + return x1 + x1 = fmod(x1, x2) + return x1 + x2 if (x1 < 0) else x1 + + +cdef class PeriodicTransform1D(Function1D): + + cdef: + readonly Function1D function1d + readonly double period + + +cdef class PeriodicTransform2D(Function2D): + + cdef: + readonly Function2D function2d + double period_x, period_y + + +cdef class PeriodicTransform3D(Function3D): + + cdef: + readonly Function3D function3d + readonly double period_x, period_y, period_z + + +cdef class VectorPeriodicTransform1D(VectorFunction1D): + + cdef: + readonly VectorFunction1D function1d + readonly double period + + +cdef class VectorPeriodicTransform2D(VectorFunction2D): + + cdef: + readonly VectorFunction2D function2d + readonly double period_x, period_y + + +cdef class VectorPeriodicTransform3D(VectorFunction3D): + + cdef: + readonly VectorFunction3D function3d + readonly double period_x, period_y, period_z diff --git a/cherab/core/math/transform/periodic.pyx b/cherab/core/math/transform/periodic.pyx new file mode 100644 index 00000000..13869458 --- /dev/null +++ b/cherab/core/math/transform/periodic.pyx @@ -0,0 +1,352 @@ +# Copyright 2016-2022 Euratom +# Copyright 2016-2022 United Kingdom Atomic Energy Authority +# Copyright 2016-2022 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from raysect.core.math cimport Vector3D +from raysect.core.math.function.float cimport autowrap_function1d, autowrap_function2d, autowrap_function3d +from raysect.core.math.function.vector3d cimport autowrap_function1d as autowrap_vectorfunction1d +from raysect.core.math.function.vector3d cimport autowrap_function2d as autowrap_vectorfunction2d +from raysect.core.math.function.vector3d cimport autowrap_function3d as autowrap_vectorfunction3d + + +cdef class PeriodicTransform1D(Function1D): + """ + Extends a periodic 1D function to an infinite 1D space. + + :param Function1D function1d: The periodic 1D function defined + in the [0, period) interval. + :param double period: The period of the function. + + .. code-block:: pycon + + >>> from cherab.core.math import PeriodicTransform1D + >>> + >>> def f1(x): + >>> return x + >>> + >>> f2 = PeriodicTransform1D(f1, 1.) + >>> + >>> f2(1.5) + 0.5 + >>> f2(-0.3) + 0.7 + """ + + def __init__(self, object function1d, double period): + + if not callable(function1d): + raise TypeError("function1d is not callable.") + + self.function1d = autowrap_function1d(function1d) + + if period <= 0: + raise ValueError("Argument period must be positive.") + + self.period = period + + cdef double evaluate(self, double x) except? -1e999: + """Return the value of periodic function.""" + + return self.function1d.evaluate(remainder(x, self.period)) + + +cdef class PeriodicTransform2D(Function2D): + """ + Extends a periodic 2D function to an infinite 2D space. + + Set period_x/period_y to 0 if the function is not periodic along x/y axis. + + :param Function2D function2d: The periodic 2D function defined + in the ([0, period_x), [0, period_y)) intervals. + :param double period_x: The period of the function along x-axis. + 0 if not periodic. + :param double period_y: The period of the function along y-axis. + 0 if not periodic. + + .. code-block:: pycon + + >>> from cherab.core.math import PeriodicTransform2D + >>> + >>> def f1(x, y): + >>> return x * y + >>> + >>> f2 = PeriodicTransform2D(f1, 1., 1.) + >>> + >>> f2(1.5, 1.5) + 0.25 + >>> f2(-0.3, -1.3) + 0.49 + >>> + >>> f3 = PeriodicTransform2D(f1, 1., 0) + >>> + >>> f3(1.5, 1.5) + 0.75 + >>> f3(-0.3, -1.3) + -0.91 + """ + + def __init__(self, object function2d, double period_x, double period_y): + + if not callable(function2d): + raise TypeError("function2d is not callable.") + + self.function2d = autowrap_function2d(function2d) + + if period_x < 0: + raise ValueError("Argument period_x must be >= 0.") + if period_y < 0: + raise ValueError("Argument period_y must be >= 0.") + + self.period_x = period_x + self.period_y = period_y + + cdef double evaluate(self, double x, double y) except? -1e999: + """Return the value of periodic function.""" + + x = remainder(x, self.period_x) + y = remainder(y, self.period_y) + + return self.function2d.evaluate(x, y) + + +cdef class PeriodicTransform3D(Function3D): + """ + Extends a periodic 3D function to an infinite 3D space. + + Set period_x/period_y/period_z to 0 if the function is not periodic along x/y/z axis. + + :param Function3D function3d: The periodic 3D function defined in the + ([0, period_x), [0, period_y), [0, period_z)) intervals. + :param double period_x: The period of the function along x-axis. + 0 if not periodic. + :param double period_y: The period of the function along y-axis. + 0 if not periodic. + :param double period_z: The period of the function along z-axis. + 0 if not periodic. + + .. code-block:: pycon + + >>> from cherab.core.math import PeriodicTransform3D + >>> + >>> def f1(x, y, z): + >>> return x * y * z + >>> + >>> f2 = PeriodicTransform3D(f1, 1., 1., 1.) + >>> + >>> f2(1.5, 1.5, 1.5) + 0.125 + >>> f2(-0.3, -1.3, -2.3) + 0.343 + >>> + >>> f3 = PeriodicTransform3D(f1, 0, 1., 0) + >>> + >>> f3(1.5, 1.5, 1.5) + 1.125 + >>> f3(-0.3, -1.3, -0.3) + 0.063 + """ + + def __init__(self, object function3d, double period_x, double period_y, double period_z): + + if not callable(function3d): + raise TypeError("function2d is not callable.") + + self.function3d = autowrap_function3d(function3d) + + if period_x < 0: + raise ValueError("Argument period_x must be >= 0.") + if period_y < 0: + raise ValueError("Argument period_y must be >= 0.") + if period_z < 0: + raise ValueError("Argument period_z must be >= 0.") + + self.period_x = period_x + self.period_y = period_y + self.period_z = period_z + + cdef double evaluate(self, double x, double y, double z) except? -1e999: + """Return the value of periodic function.""" + + x = remainder(x, self.period_x) + y = remainder(y, self.period_y) + z = remainder(z, self.period_z) + + return self.function3d.evaluate(x, y, z) + + +cdef class VectorPeriodicTransform1D(VectorFunction1D): + """ + Extends a periodic 1D vector function to an infinite 1D space. + + :param VectorFunction1D function1d: The periodic 1D vector function + defined in the [0, period) interval. + :param double period: The period of the function. + + .. code-block:: pycon + + >>> from raysect.core.math import Vector3D + >>> from cherab.core.math import VectorPeriodicTransform1D + >>> + >>> def f1(x): + >>> return Vector3D(x, 0, 0) + >>> + >>> f2 = VectorPeriodicTransform1D(f1, 1.) + >>> + >>> f2(1.5) + Vector3D(0.5, 0, 0) + >>> f2(-0.3) + Vector3D(0.7, 0, 0) + """ + + def __init__(self, object function1d, double period): + + if not callable(function1d): + raise TypeError("function1d is not callable.") + + self.function1d = autowrap_vectorfunction1d(function1d) + + if period <= 0: + raise ValueError("Argument period must be positive.") + + self.period = period + + cdef Vector3D evaluate(self, double x): + """Return the value of periodic function.""" + + return self.function1d.evaluate(remainder(x, self.period)) + + +cdef class VectorPeriodicTransform2D(VectorFunction2D): + """ + Extends a periodic 2D vector function to an infinite 2D space. + + Set period_x/period_y to 0 if the function is not periodic along x/y axis. + + :param VectorFunction2D function2d: The periodic 2D vector function defined in + the ([0, period_x), [0, period_y)) intervals. + :param double period_x: The period of the function along x-axis. + 0 if not periodic. + :param double period_y: The period of the function along y-axis. + 0 if not periodic. + + .. code-block:: pycon + + >>> from cherab.core.math import VectorPeriodicTransform2D + >>> + >>> def f1(x, y): + >>> return Vector3D(x, y, 0) + >>> + >>> f2 = VectorPeriodicTransform2D(f1, 1., 1.) + >>> + >>> f2(1.5, 1.5) + Vector3D(0.5, 0.5, 0) + >>> f2(-0.3, -1.3) + Vector3D(0.7, 0.7, 0) + >>> + >>> f3 = VectorPeriodicTransform2D(f1, 1., 0) + >>> + >>> f3(1.5, 1.5) + Vector3D(0.5, 1.5, 0) + >>> f3(-0.3, -1.3) + Vector3D(0.7, -1.3, 0) + """ + + def __init__(self, object function2d, double period_x, double period_y): + + if not callable(function2d): + raise TypeError("function2d is not callable.") + + self.function2d = autowrap_vectorfunction2d(function2d) + + if period_x < 0: + raise ValueError("Argument period_x must be >= 0.") + if period_y < 0: + raise ValueError("Argument period_y must be >= 0.") + + self.period_x = period_x + self.period_y = period_y + + cdef Vector3D evaluate(self, double x, double y): + """Return the value of periodic function.""" + + x = remainder(x, self.period_x) + y = remainder(y, self.period_y) + + return self.function2d.evaluate(x, y) + + +cdef class VectorPeriodicTransform3D(VectorFunction3D): + """ + Extends a periodic 3D vector function to an infinite 3D space. + + Set period_x/period_y/period_z to 0 if the function is not periodic along x/y/z axis. + + :param VectorFunction3D function3d: The periodic 3D vector function defined in the + ([0, period_x), [0, period_y), [0, period_z)) intervals. + :param double period_x: The period of the function along x-axis. + 0 if not periodic. + :param double period_y: The period of the function along y-axis. + 0 if not periodic. + :param double period_z: The period of the function along z-axis. + 0 if not periodic. + + .. code-block:: pycon + + >>> from cherab.core.math import PeriodicTransform3D + >>> + >>> def f1(x, y, z): + >>> return Vector3D(x, y, z) + >>> + >>> f2 = VectorPeriodicTransform3D(f1, 1., 1., 1.) + >>> + >>> f2(1.5, 1.5, 1.5) + Vector3D(0.5, 0.5, 0.5) + >>> f2(-0.3, -1.3, -2.3) + Vector3D(0.7, 0.7, 0.7) + >>> + >>> f3 = VectorPeriodicTransform3D(f1, 0, 1., 0) + >>> + >>> f3(1.5, 0.5, 1.5) + Vector3D(1.5, 0.5, 1.5) + """ + + def __init__(self, object function3d, double period_x, double period_y, double period_z): + + if not callable(function3d): + raise TypeError("function2d is not callable.") + + self.function3d = autowrap_vectorfunction3d(function3d) + + if period_x < 0: + raise ValueError("Argument period_x must be >= 0.") + if period_y < 0: + raise ValueError("Argument period_y must be >= 0.") + if period_z < 0: + raise ValueError("Argument period_z must be >= 0.") + + self.period_x = period_x + self.period_y = period_y + self.period_z = period_z + + cdef Vector3D evaluate(self, double x, double y, double z): + """Return the value of periodic function.""" + + x = remainder(x, self.period_x) + y = remainder(y, self.period_y) + z = remainder(z, self.period_z) + + return self.function3d.evaluate(x, y, z) From 04b7813c82f577e5438cd985ebe6462cb86685ed Mon Sep 17 00:00:00 2001 From: vsnever Date: Wed, 28 Dec 2022 11:59:45 +0300 Subject: [PATCH 5/9] Add docs for math.transfrom sub-module. --- docs/source/math/math.rst | 1 + docs/source/math/transform.rst | 11 +++++++++++ 2 files changed, 12 insertions(+) create mode 100644 docs/source/math/transform.rst diff --git a/docs/source/math/math.rst b/docs/source/math/math.rst index eab9a1fa..836cbe84 100644 --- a/docs/source/math/math.rst +++ b/docs/source/math/math.rst @@ -18,6 +18,7 @@ utilities that Cherab provides for slicing, dicing and projecting these function function interpolators mappers + transform mask samplers slice diff --git a/docs/source/math/transform.rst b/docs/source/math/transform.rst new file mode 100644 index 00000000..f3d6cfa8 --- /dev/null +++ b/docs/source/math/transform.rst @@ -0,0 +1,11 @@ + +Transformations +--------------- + +These functions perform coordinate transformations as wrappers for other functions. Unlike mappers, these wrappers do not change the dimensionality of the wrapped functions. + +.. automodule:: cherab.core.math.transform.cylindrical + :members: + +.. automodule:: cherab.core.math.transform.periodic + :members: From 751f324960afdb7be2402a40f702c1721427e83a Mon Sep 17 00:00:00 2001 From: Vladislav Neverov Date: Tue, 14 Mar 2023 10:57:47 +0300 Subject: [PATCH 6/9] Add tests for VectorPeriodicTransform#D functions. --- cherab/core/math/tests/test_transform.py | 115 +++++++++++++++++++++-- 1 file changed, 109 insertions(+), 6 deletions(-) diff --git a/cherab/core/math/tests/test_transform.py b/cherab/core/math/tests/test_transform.py index 2ee4eefc..6fad98b0 100644 --- a/cherab/core/math/tests/test_transform.py +++ b/cherab/core/math/tests/test_transform.py @@ -28,10 +28,12 @@ class TestCylindricalTransform(unittest.TestCase): def setUp(self): """Initialisation with functions to map.""" - def f3d(r, phi, z): return r * np.cos(phi) + z + def f3d(r, phi, z): + return r * np.cos(phi) + z self.function3d = f3d - def vecf3d(r, phi, z): return Vector3D(np.sin(phi), r * z, np.cos(phi)) + def vecf3d(r, phi, z): + return Vector3D(np.sin(phi), r * z, np.cos(phi)) self.vectorfunction3d = vecf3d def test_cylindrical_transform(self): @@ -63,13 +65,30 @@ class TestPeriodicTransform(unittest.TestCase): def setUp(self): """Initialisation with functions to map.""" - def f1d(x): return x*np.cos(x-3) + def f1d(x): + return x * np.cos(x - 3) self.function1d = f1d - def f2d(x, y): return x*np.sin(y) + + def f2d(x, y): + return x * np.sin(y) self.function2d = f2d - def f3d(x, y, z): return x*x*np.exp(y)-2*z*y + + def f3d(x, y, z): + return x * x * np.exp(y) - 2 * z * y self.function3d = f3d + def vecf1d(x): + return Vector3D(x, x**2, x**3) + self.vectorfunction1d = vecf1d + + def vecf2d(x, y): + return Vector3D(x, y, x * y) + self.vectorfunction2d = vecf2d + + def vecf3d(x, y, z): + return Vector3D(x + y + z, (x + y) * z, x * y * z) + self.vectorfunction3d = vecf3d + def test_periodic_transform_1d(self): """1D periodic transform""" period_func = transform.PeriodicTransform1D(self.function1d, np.pi) @@ -140,7 +159,7 @@ def test_periodic_transform_3d(self): places=10) def test_periodic_transform_3d_invalid_arg(self): - """2D periodic transform. Invalid arguments.""" + """3D periodic transform. Invalid arguments.""" # 1st argument is not callable self.assertRaises(TypeError, transform.PeriodicTransform3D, "blah", np.pi, np.pi, np.pi) # period is not a number @@ -152,6 +171,90 @@ def test_periodic_transform_3d_invalid_arg(self): self.assertRaises(ValueError, transform.PeriodicTransform3D, self.function3d, np.pi, -1, np.pi) self.assertRaises(ValueError, transform.PeriodicTransform3D, self.function3d, np.pi, np.pi, -1) + def test_vector_periodic_transform_1d(self): + """1D vector periodic transform""" + period_func = transform.VectorPeriodicTransform1D(self.vectorfunction1d, 1) + vec1 = period_func(1.4) + vec2 = Vector3D(0.4, 0.16, 0.064) + np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) + + def test_vector_periodic_transform_1d_invalid_arg(self): + """1D vector periodic transform. Invalid arguments.""" + # 1st argument is not callable + self.assertRaises(TypeError, transform.VectorPeriodicTransform1D, "blah", 1.) + # period is not a number + self.assertRaises(TypeError, transform.VectorPeriodicTransform1D, self.vectorfunction1d, "blah") + # period is negative + self.assertRaises(ValueError, transform.VectorPeriodicTransform1D, self.vectorfunction1d, -1) + + def test_vector_periodic_transform_2d(self): + """2D vector periodic transform""" + period_func = transform.VectorPeriodicTransform2D(self.vectorfunction2d, 1, 1) + vec1 = period_func(-0.4, 1.6) + vec2 = Vector3D(0.6, 0.6, 0.36) + np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) + + # Periodic only along x + period_func = transform.VectorPeriodicTransform2D(self.vectorfunction2d, 1, 0) + vec1 = period_func(-0.4, 1.6) + vec2 = Vector3D(0.6, 1.6, 0.96) + np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) + + # Periodic only along y + period_func = transform.VectorPeriodicTransform2D(self.vectorfunction2d, 0, 1) + vec1 = period_func(-0.4, 1.6) + vec2 = Vector3D(-0.4, 0.6, -0.24) + np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) + + def test_vector_periodic_transform_2d_invalid_arg(self): + """2D vector periodic transform. Invalid arguments.""" + # 1st argument is not callable + self.assertRaises(TypeError, transform.VectorPeriodicTransform2D, "blah", 1, 1) + # period is not a number + self.assertRaises(TypeError, transform.VectorPeriodicTransform2D, self.function2d, "blah", 1) + self.assertRaises(TypeError, transform.VectorPeriodicTransform2D, self.function2d, 1, "blah") + # period is negative + self.assertRaises(ValueError, transform.VectorPeriodicTransform2D, self.function2d, -1, 1) + self.assertRaises(ValueError, transform.VectorPeriodicTransform2D, self.function2d, 1, -1) + + def test_vector_periodic_transform_3d(self): + """3D vector periodic transform""" + period_func = transform.VectorPeriodicTransform3D(self.vectorfunction3d, 1, 1, 1) + vec1 = period_func(-0.4, 1.6, 1.2) + vec2 = Vector3D(1.4, 0.24, 0.072) + np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) + + # Periodic along y and z + period_func = transform.VectorPeriodicTransform3D(self.vectorfunction3d, 0, 1, 1) + vec1 = period_func(-0.4, 1.6, 1.2) + vec2 = Vector3D(0.4, 0.04, -0.048) + np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) + + # Periodic along x and z + period_func = transform.VectorPeriodicTransform3D(self.vectorfunction3d, 1, 0, 1) + vec1 = period_func(-0.4, 1.6, 1.2) + vec2 = Vector3D(2.4, 0.44, 0.192) + np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) + + # Periodic along x and y + period_func = transform.VectorPeriodicTransform3D(self.vectorfunction3d, 1, 1, 0) + vec1 = period_func(-0.4, 1.6, 1.2) + vec2 = Vector3D(2.4, 1.44, 0.432) + np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) + + def test_vector_periodic_transform_3d_invalid_arg(self): + """2D vector periodic transform. Invalid arguments.""" + # 1st argument is not callable + self.assertRaises(TypeError, transform.VectorPeriodicTransform3D, "blah", 1, 1, 1) + # period is not a number + self.assertRaises(TypeError, transform.VectorPeriodicTransform3D, self.vectorfunction3d, "blah", 1, 1) + self.assertRaises(TypeError, transform.VectorPeriodicTransform3D, self.vectorfunction3d, 1, "blah", 1) + self.assertRaises(TypeError, transform.VectorPeriodicTransform3D, self.vectorfunction3d, 1, 1, "blah") + # period is negative + self.assertRaises(ValueError, transform.VectorPeriodicTransform3D, self.vectorfunction3d, -1, 1, 1) + self.assertRaises(ValueError, transform.VectorPeriodicTransform3D, self.vectorfunction3d, 1, -1, 1) + self.assertRaises(ValueError, transform.VectorPeriodicTransform3D, self.vectorfunction3d, 1, 1, -1) + if __name__ == '__main__': unittest.main() From ab14fb8b74a61e0745553b37f6de14071bdebb5e Mon Sep 17 00:00:00 2001 From: vsnever Date: Tue, 14 Mar 2023 12:57:22 +0300 Subject: [PATCH 7/9] Fixed copy-paste errors in TestPeriodicTransform. --- cherab/core/math/tests/test_transform.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cherab/core/math/tests/test_transform.py b/cherab/core/math/tests/test_transform.py index 6fad98b0..fc3a2445 100644 --- a/cherab/core/math/tests/test_transform.py +++ b/cherab/core/math/tests/test_transform.py @@ -211,11 +211,11 @@ def test_vector_periodic_transform_2d_invalid_arg(self): # 1st argument is not callable self.assertRaises(TypeError, transform.VectorPeriodicTransform2D, "blah", 1, 1) # period is not a number - self.assertRaises(TypeError, transform.VectorPeriodicTransform2D, self.function2d, "blah", 1) - self.assertRaises(TypeError, transform.VectorPeriodicTransform2D, self.function2d, 1, "blah") + self.assertRaises(TypeError, transform.VectorPeriodicTransform2D, self.vectorfunction2d, "blah", 1) + self.assertRaises(TypeError, transform.VectorPeriodicTransform2D, self.vectorfunction2d, 1, "blah") # period is negative - self.assertRaises(ValueError, transform.VectorPeriodicTransform2D, self.function2d, -1, 1) - self.assertRaises(ValueError, transform.VectorPeriodicTransform2D, self.function2d, 1, -1) + self.assertRaises(ValueError, transform.VectorPeriodicTransform2D, self.vectorfunction2d, -1, 1) + self.assertRaises(ValueError, transform.VectorPeriodicTransform2D, self.vectorfunction2d, 1, -1) def test_vector_periodic_transform_3d(self): """3D vector periodic transform""" @@ -243,7 +243,7 @@ def test_vector_periodic_transform_3d(self): np.testing.assert_almost_equal([vec1.x, vec1.y, vec1.z], [vec2.x, vec2.y, vec2.z], decimal=10) def test_vector_periodic_transform_3d_invalid_arg(self): - """2D vector periodic transform. Invalid arguments.""" + """3D vector periodic transform. Invalid arguments.""" # 1st argument is not callable self.assertRaises(TypeError, transform.VectorPeriodicTransform3D, "blah", 1, 1, 1) # period is not a number From 101456789d8ac1dcf47a8bfb17c0a6e8479358d5 Mon Sep 17 00:00:00 2001 From: "Lovell, Jack J" Date: Fri, 17 Jan 2020 13:37:52 +0000 Subject: [PATCH 8/9] Add analytic plasma demo using raysect's function framework This is roughly 25x faster than the pure python, class-based demo. --- .../analytic_plasma_function_framework.py | 165 ++++++++++++++++++ docs/source/conf.py | 2 + .../plasmas/analytic_function_plasma.rst | 22 ++- 3 files changed, 185 insertions(+), 4 deletions(-) create mode 100644 demos/plasmas/analytic_plasma_function_framework.py diff --git a/demos/plasmas/analytic_plasma_function_framework.py b/demos/plasmas/analytic_plasma_function_framework.py new file mode 100644 index 00000000..1459f28e --- /dev/null +++ b/demos/plasmas/analytic_plasma_function_framework.py @@ -0,0 +1,165 @@ + +# Copyright 2016-2018 Euratom +# Copyright 2016-2018 United Kingdom Atomic Energy Authority +# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + + +import numpy as np +import matplotlib.pyplot as plt +from scipy.constants import electron_mass, atomic_mass + +from raysect.core.math.function.float import Arg2D, Exp2D, Sqrt2D +from raysect.primitive import Cylinder +from raysect.optical import World, translate, Point3D, Vector3D, rotate_basis, Spectrum +from raysect.optical.observer import PinholeCamera, PowerPipeline2D + +from cherab.core import Species, Maxwellian, Plasma, Line +from cherab.core.math import sample3d, AxisymmetricMapper +from cherab.core.atomic import deuterium +from cherab.core.model import ExcitationLine +from cherab.openadas import OpenADAS + + +def NeutralFunction(peak_value, sigma, magnetic_axis, lcfs_radius=1): + """A neutral profile that is constant outside the plasma, + then exponentially decays inside the LCFS.""" + raxis = magnetic_axis[0] + zaxis = magnetic_axis[1] + radius_from_axis = Sqrt2D((Arg2D('x') - raxis)**2 + (Arg2D('y') - zaxis)**2) + scale = Exp2D(-((radius_from_axis - lcfs_radius)**2) / (2 * sigma**2)) + inside_lcfs = (radius_from_axis <= lcfs_radius) + # density = peak * scale * inside_lcfs + peak * (inside_lcfs - 1). + # Rearrange so inside_lcfs and scale are only called once each. + density = peak_value * (inside_lcfs * (scale - 1) + 1) + return AxisymmetricMapper(density) + + +def IonFunction(v_core, v_lcfs, magnetic_axis, p=4, q=3, lcfs_radius=1): + """An approximate toroidal plasma profile that follows a double + quadratic between the LCFS and magnetic axis.""" + r_axis = magnetic_axis[0] + z_axis = magnetic_axis[1] + radius_from_axis = Sqrt2D((Arg2D('x') - r_axis)**2 + (Arg2D('y') - z_axis)**2) + density = (v_core - v_lcfs) * (1 - (radius_from_axis / lcfs_radius)**p)**q + v_lcfs + inside_lcfs = (radius_from_axis <= lcfs_radius) + return AxisymmetricMapper(density * inside_lcfs) + + +# tunables +peak_density = 1e19 +peak_temperature = 2500 +magnetic_axis = (2.5, 0) + + +# setup scenegraph +world = World() + + +################### +# plasma creation # + +plasma = Plasma(parent=world) +plasma.atomic_data = OpenADAS(permit_extrapolation=True) +plasma.geometry = Cylinder(3.5, 2.2, transform=translate(0, 0, -1.1)) +plasma.geometry_transform = translate(0, 0, -1.1) + +# No net velocity for any species +zero_velocity = Vector3D(0, 0, 0) + +# define neutral species distribution +d0_density = NeutralFunction(peak_density, 0.1, magnetic_axis) +d0_temperature = 0.5 # constant 0.5eV temperature for all neutrals +d0_distribution = Maxwellian(d0_density, d0_temperature, zero_velocity, + deuterium.atomic_weight * atomic_mass) +d0_species = Species(deuterium, 0, d0_distribution) + +# define deuterium ion species distribution +d1_density = IonFunction(peak_density, 0, magnetic_axis) +d1_temperature = IonFunction(peak_temperature, 0, magnetic_axis) +d1_distribution = Maxwellian(d1_density, d1_temperature, zero_velocity, + deuterium.atomic_weight * atomic_mass) +d1_species = Species(deuterium, 1, d1_distribution) + +# define the electron distribution +e_density = IonFunction(peak_density, 0, magnetic_axis) +e_temperature = IonFunction(peak_temperature, 0, magnetic_axis) +e_distribution = Maxwellian(e_density, e_temperature, zero_velocity, electron_mass) + +# define species +plasma.b_field = Vector3D(1.0, 1.0, 1.0) +plasma.electron_distribution = e_distribution +plasma.composition = [d0_species, d1_species] + +# Add a balmer alpha line for visualisation purposes +d_alpha_excit = ExcitationLine(Line(deuterium, 0, (3, 2))) +plasma.models = [d_alpha_excit] + + +#################### +# Visualise Plasma # + +# Run some plots to check the distribution functions and emission profile are as expected +r, _, z, t_samples = sample3d(d1_temperature, (0, 4, 200), (0, 0, 1), (-2, 2, 200)) +plt.imshow(np.transpose(np.squeeze(t_samples)), extent=[0, 4, -2, 2]) +plt.colorbar() +plt.axis('equal') +plt.xlabel('r axis') +plt.ylabel('z axis') +plt.title("Ion temperature profile in r-z plane") + +plt.figure() +r, _, z, t_samples = sample3d(d1_temperature, (-4, 4, 400), (-4, 4, 400), (0, 0, 1)) +plt.imshow(np.transpose(np.squeeze(t_samples)), extent=[-4, 4, -4, 4]) +plt.colorbar() +plt.axis('equal') +plt.xlabel('x axis') +plt.ylabel('y axis') +plt.title("Ion temperature profile in x-y plane") + +plt.figure() +r, _, z, t_samples = sample3d(d0_density, (0, 4, 200), (0, 0, 1), (-2, 2, 200)) +plt.imshow(np.transpose(np.squeeze(t_samples)), extent=[0, 4, -2, 2]) +plt.colorbar() +plt.axis('equal') +plt.xlabel('r axis') +plt.ylabel('z axis') +plt.title("Neutral Density profile in r-z plane") + +plt.figure() +xrange = np.linspace(0, 4, 200) +yrange = np.linspace(-2, 2, 200) +d_alpha_rz_intensity = np.zeros((200, 200)) +direction = Vector3D(0, 1, 0) +for i, x in enumerate(xrange): + for j, y in enumerate(yrange): + emission = d_alpha_excit.emission(Point3D(x, 0.0, y), direction, Spectrum(650, 660, 1)) + d_alpha_rz_intensity[j, i] = emission.total() +plt.imshow(d_alpha_rz_intensity, extent=[0, 4, -2, 2], origin='lower') +plt.colorbar() +plt.xlabel('r axis') +plt.ylabel('z axis') +plt.title("D-alpha emission in R-Z") + + +camera = PinholeCamera((256, 256), pipelines=[PowerPipeline2D()], parent=world) +camera.transform = translate(2.5, -4.5, 0)*rotate_basis(Vector3D(0, 1, 0), Vector3D(0, 0, 1)) +camera.pixel_samples = 1 + +plt.ion() +camera.observe() +plt.ioff() +plt.show() diff --git a/docs/source/conf.py b/docs/source/conf.py index 60589037..d17ec793 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -39,6 +39,7 @@ extensions = [ 'sphinx.ext.autodoc', 'sphinx.ext.mathjax', + 'sphinx_tabs.tabs', ] # Add any paths that contain templates here, relative to this directory. @@ -142,6 +143,7 @@ html_context = {'css_files': [ '_static/theme_overrides.css', # override wide tables in RTD theme + '_static/tabs.css', ] } diff --git a/docs/source/demonstrations/plasmas/analytic_function_plasma.rst b/docs/source/demonstrations/plasmas/analytic_function_plasma.rst index e74ba874..627f0c66 100644 --- a/docs/source/demonstrations/plasmas/analytic_function_plasma.rst +++ b/docs/source/demonstrations/plasmas/analytic_function_plasma.rst @@ -10,11 +10,25 @@ shows how to use these functions in a plasma and visualises the results. Note that while it is possible to use pure python functions for development, they are typically ~100 times slower than their cython counterparts. Therefore, for use cases where speed is important -we recommend moving these functions to cython classes. For an example of a cython function, -see the `Gaussian `_ -cython class in the demos folder. +we recommend moving these functions to cython classes. An alternative solution which may not require +writing and compiling any additional cython code is to use Raysect's +`function framework `_ to build up +expressions which will be evaluated like Python functions. These will typically run slightly slower +than a hand-coded cython implementation but still significantly faster than a pure python +implementation. -.. literalinclude:: ../../../../demos/plasmas/analytic_plasma.py +Two examples are provided, one using a pure python implementation of analytic forms for neutral and +ion plasma species distributions, and one using objects from Raysect's function framework. + +.. tabs:: + + .. tab:: Pure python + + .. literalinclude:: ../../../../demos/plasmas/analytic_plasma.py + + .. tab:: Function framework + + .. literalinclude:: ../../../../demos/plasmas/analytic_plasma_function_framework.py .. figure:: analytic_plasma_slices.png :align: center From ac6863a5e548f6b06e4b60c5a66939e57e068472 Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Thu, 6 Apr 2023 15:48:16 +0100 Subject: [PATCH 9/9] Add a docs extras requirement for building documentation When building the docs, run `pip install cherab[docs]` to get the dependencies required for successfully running `./dev/build_docs.sh`. --- setup.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/setup.py b/setup.py index f0610a6c..5d9d7a8e 100644 --- a/setup.py +++ b/setup.py @@ -108,6 +108,10 @@ "matplotlib", "raysect==0.8.1", ], + extras_require={ + # Running ./dev/build_docs.sh runs setup.py, which requires cython. + "docs": ["cython", "sphinx", "sphinx-rtd-theme", "sphinx-tabs"], + }, packages=find_packages(), include_package_data=True, zip_safe=False,