Source code for xrayutilities.math.vector

# This file is part of xrayutilities.
#
# xrayutilities is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
# Copyright (C) 2009 Eugen Wintersberger <eugen.wintersberger@desy.de>
# Copyright (C) 2010-2020 Dominik Kriegner <dominik.kriegner@gmail.com>

"""
module with vector operations for vectors of size 3,
since for so short vectors numpy does not give the best performance explicit
implementation of the equations is performed together with error checking to
ensure vectors of length 3.
"""

import math
import re

import numpy

from ..exception import InputError

circleSyntax = re.compile("[xyz][+-]")


def _checkvec(v):
    if isinstance(v, (list, tuple, numpy.ndarray)):
        vtmp = numpy.asarray(v, dtype=numpy.double)
    else:
        raise TypeError("Vector must be a list, tuple or numpy array")
    return vtmp


[docs]def VecNorm(v): """ Calculate the norm of a vector. Parameters ---------- v : list or array-like input vector(s), either one vector or an array of vectors with shape (n, 3) Returns ------- float or ndarray vector norm, either a single float or shape (n, ) """ if isinstance(v, numpy.ndarray): if len(v.shape) >= 2: return numpy.linalg.norm(v, axis=-1) if len(v) != 3: raise ValueError("Vector must be of length 3, but has length %d!" % len(v)) return math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
[docs]def VecUnit(v): """ Calculate the unit vector of v. Parameters ---------- v : list or array-like input vector(s), either one vector or an array of vectors with shape (n, 3) Returns ------- ndarray unit vector of `v`, either shape (3, ) or (n, 3) """ vtmp = _checkvec(v) if len(vtmp.shape) == 1: return vtmp / VecNorm(vtmp) else: return vtmp / VecNorm(vtmp)[..., numpy.newaxis]
[docs]def VecDot(v1, v2): """ Calculate the vector dot product. Parameters ---------- v1, v2 : list or array-like input vector(s), either one vector or an array of vectors with shape (n, 3) Returns ------- float or ndarray innter product of the vectors, either a single float or (n, ) """ if isinstance(v1, numpy.ndarray): if len(v1.shape) >= 2: return numpy.einsum('...i, ...i', v1, v2) if len(v1) != 3 or len(v2) != 3: raise ValueError("Vectors must be of size 3! (len(v1)=%d len(v2)=%d)" % (len(v1), len(v2))) return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]
[docs]def VecCross(v1, v2, out=None): """ Calculate the vector cross product. Parameters ---------- v1, v2 : list or array-like input vector(s), either one vector or an array of vectors with shape (n, 3) out : list or array-like, optional output vector Returns ------- ndarray cross product either of shape (3, ) or (n, 3) """ if isinstance(v1, numpy.ndarray): if len(v1.shape) >= 2 or len(v2.shape) >= 2: return numpy.cross(v1, v2) if len(v1) != 3 or len(v2) != 3: raise ValueError("Vectors must be of size 3! (len(v1)=%d len(v2)=%d)" % (len(v1), len(v2))) if out is None: out = numpy.empty(3) out[0] = v1[1] * v2[2] - v1[2] * v2[1] out[1] = v1[2] * v2[0] - v1[0] * v2[2] out[2] = v1[0] * v2[1] - v1[1] * v2[0] return out
[docs]def VecAngle(v1, v2, deg=False): """ calculate the angle between two vectors. The following formula is used v1.v2 = norm(v1)*norm(v2)*cos(alpha) alpha = acos((v1.v2)/(norm(v1)*norm(v2))) Parameters ---------- v1, v2 : list or array-like input vector(s), either one vector or an array of vectors with shape (n, 3) deg: bool, optional True: return result in degree, False: in radiants (default: False) Returns ------- float or ndarray the angle included by the two vectors `v1` and `v2`, either a single float or an array with shape (n, ) """ u1 = VecNorm(v1) u2 = VecNorm(v2) if isinstance(u1, numpy.ndarray) or isinstance(u2, numpy.ndarray): s = VecDot(v1, v2) / u1 / u2 s[numpy.abs(s) > 1.0] = numpy.sign(s[numpy.abs(s) > 1.0]) * 1.0 alpha = numpy.arccos(s) if deg: alpha = numpy.degrees(alpha) else: alpha = math.acos(max(min(1., VecDot(v1, v2) / u1 / u2), -1.)) if deg: alpha = math.degrees(alpha) return alpha
[docs]def distance(x, y, z, point, vec): """ calculate the distance between the point (x, y, z) and the line defined by the point and vector vec Parameters ---------- x : float or ndarray x coordinate(s) of the point(s) y : float or ndarray y coordinate(s) of the point(s) z : float or ndarray z coordinate(s) of the point(s) point : tuple, list or ndarray 3D point on the line to which the distance should be calculated vec : tuple, list or ndarray 3D vector defining the propergation direction of the line """ coords = numpy.vstack((numpy.ravel(x) - point[0], numpy.ravel(y) - point[1], numpy.ravel(z) - point[2])).T ret = VecNorm(VecCross(coords, numpy.asarray(vec)))/VecNorm(vec) if isinstance(x, numpy.ndarray): ret = ret.reshape(x.shape) else: ret = ret[0] return ret
[docs]def getVector(string): """ returns unit vector along a rotation axis given in the syntax 'x+' 'z-' or equivalents Parameters ---------- string: str vector string following the synthax [xyz][+-] Returns ------- ndarray vector along the given direction """ if len(string) != 2: raise InputError("wrong length of string for conversion given") if not circleSyntax.search(string): raise InputError("getVector: incorrect string syntax (%s)" % string) if string[0] == 'x': v = [1., 0, 0] elif string[0] == 'y': v = [0, 1., 0] elif string[0] == 'z': v = [0, 0, 1.] else: raise InputError("wrong first character of string given " "(needs to be one of x, y, z)") if string[1] == '+': v = numpy.asarray(v) * (+1) elif string[1] == '-': v = numpy.asarray(v) * (-1) else: raise InputError("wrong second character of string given " "(needs to be + or -)") return v
[docs]def getSyntax(vec): """ returns vector direction in the syntax 'x+' 'z-' or equivalents therefore works only for principle vectors of the coordinate system like e.g. [1, 0, 0] or [0, 2, 0] Parameters ---------- vec : list or array-like vector of length 3 Returns ------- str vector string following the synthax [xyz][+-] """ v = _checkvec(vec) if len(v) != 3: raise InputError("no valid 3D vector given") x = [1, 0, 0] y = [0, 1, 0] z = [0, 0, 1] if numpy.isclose(VecNorm(numpy.cross(numpy.cross(x, y), v)), 0): if v[2] >= 0: string = 'z+' else: string = 'z-' elif numpy.isclose(VecNorm(numpy.cross(numpy.cross(x, z), v)), 0): if v[1] >= 0: string = 'y+' else: string = 'y-' elif numpy.isclose(VecNorm(numpy.cross(numpy.cross(y, z), v)), 0): if v[0] >= 0: string = 'x+' else: string = 'x-' else: raise InputError("no valid 3D vector given") return string