# Source code for halotools.utils.vector_utilities

r"""
A set of vector calculations to aid in rotation calculations
"""

from __future__ import (division, print_function, absolute_import,
unicode_literals)
import numpy as np

__all__=['elementwise_dot',
'elementwise_norm',
'normalized_vectors',
'angles_between_list_of_vectors',
'vectors_normal_to_planes',
'project_onto_plane']
__author__ = ['Duncan Campbell', 'Andrew Hearin']

[docs]def normalized_vectors(vectors):
r"""
Return a unit-vector for each n-dimensional vector in the input list of n-dimensional points.

Parameters
----------
x : ndarray
Numpy array of shape (npts, ndim) storing a collection of n-dimensional points

Returns
-------
normed_x : ndarray
Numpy array of shape (npts, ndim)

Examples
--------
Let's create a set of semi-random 3D unit vectors.

>>> npts = int(1e3)
>>> ndim = 3
>>> x = np.random.random((npts, ndim))
>>> normed_x = normalized_vectors(x)
"""

vectors = np.atleast_2d(vectors)
npts = vectors.shape

with np.errstate(divide='ignore', invalid='ignore'):
return vectors/elementwise_norm(vectors).reshape((npts, -1))

[docs]def elementwise_norm(x):
r"""
Calculate the normalization of each element in a list of n-dimensional points.

Parameters
----------
x : ndarray
Numpy array of shape (npts, ndim) storing a collection of n-dimensional points

Returns
-------
result : ndarray
Numpy array of shape (npts, ) storing the norm of each n-dimensional point in x.

Examples
--------
>>> npts = int(1e3)
>>> ndim = 3
>>> x = np.random.random((npts, ndim))
>>> norms = elementwise_norm(x)
"""

x = np.atleast_2d(x)
return np.sqrt(np.sum(x**2, axis=1))

[docs]def elementwise_dot(x, y):
r"""
Calculate the dot product between
each pair of elements in two input lists of n-dimensional points.

Parameters
----------
x : ndarray
Numpy array of shape (npts, ndim) storing a collection of n-dimensional vectors

y : ndarray
Numpy array of shape (npts, ndim) storing a collection of n-dimensional vectors

Returns
-------
result : ndarray
Numpy array of shape (npts, ) storing the dot product between each
pair of corresponding vectors in x and y.

Examples
--------
Let's create two sets of semi-random 3D vectors, x1 and x2.

>>> npts = int(1e3)
>>> ndim = 3
>>> x1 = np.random.random((npts, ndim))
>>> x2 = np.random.random((npts, ndim))

We then can find the dot product between each pair of vectors in x1 and x2.

>>> dots = elementwise_dot(x1, x2)
"""

x = np.atleast_2d(x)
y = np.atleast_2d(y)
return np.sum(x*y, axis=1)

[docs]def angles_between_list_of_vectors(v0, v1, tol=1e-3, vn=None):
r""" Calculate the angle between a collection of n-dimensional vectors

Parameters
----------
v0 : ndarray
Numpy array of shape (npts, ndim) storing a collection of ndim-D vectors
Note that the normalization of v0 will be ignored.

v1 : ndarray
Numpy array of shape (npts, ndim) storing a collection of ndim-D vectors
Note that the normalization of v1 will be ignored.

tol : float, optional
Acceptable numerical error for errors in angle.
This variable is only used to round off numerical noise that otherwise
causes exceptions to be raised by the inverse cosine function.
Default is 0.001.

n1 : ndarray
normal vector

Returns
-------
angles : ndarray
Numpy array of shape (npts, ) storing the angles between each pair of
corresponding points in v0 and v1.

Returned values are in units of radians spanning [0, pi].

Examples
--------
Let's create two sets of semi-random 3D unit vectors.

>>> npts = int(1e4)
>>> ndim = 3
>>> v1 = np.random.random((npts, ndim))
>>> v2 = np.random.random((npts, ndim))

We then can find the angle between each pair of vectors in v1 and v2.

>>> angles = angles_between_list_of_vectors(v1, v2)
"""

dot = elementwise_dot(normalized_vectors(v0), normalized_vectors(v1))

if vn is None:
#  Protect against tiny numerical excesses beyond the range [-1 ,1]
mask1 = (dot > 1) & (dot < 1 + tol)
mask2 = (dot < -1) & (dot > -1 - tol)
a = np.arccos(dot)
else:
cross = np.cross(v0,v1)
a = np.arctan2(elementwise_dot(cross, vn), dot)

return a

[docs]def vectors_normal_to_planes(x, y):
r"""
Given a collection of 3d vectors x and y, return a collection of
3d unit-vectors that are orthogonal to x and y.

Parameters
----------
x : ndarray
Numpy array of shape (npts, 3) storing a collection of 3d vectors

Note that the normalization of x will be ignored.

y : ndarray
Numpy array of shape (npts, 3) storing a collection of 3d vectors

Note that the normalization of y will be ignored.

Returns
-------
z : ndarray
Numpy array of shape (npts, 3). Each 3d vector in z will be orthogonal
to the corresponding vector in x and y.

Examples
--------
Define a set of random 3D vectors

>>> npts = int(1e4)
>>> x = np.random.random((npts, 3))
>>> y = np.random.random((npts, 3))

now calculate a thrid set of vectors to a corresponding pair in x and y.

>>> normed_z = angles_between_list_of_vectors(x, y)
"""
return normalized_vectors(np.cross(x, y))

[docs]def project_onto_plane(x1, x2):
r"""
Given a collection of vectors, x1 and x2, project each vector
in x1 onto the plane normal to the corresponding vector x2.

Parameters
----------
x1 : ndarray
Numpy array of shape (npts, 3) storing a collection of 3d points

x2 : ndarray
Numpy array of shape (npts, 3) storing a collection of 3d points

Returns
-------
result : ndarray
Numpy array of shape (npts, 3) storing a collection of 3d points

Examples
--------
Define a set of random 3D vectors.

>>> npts = int(1e4)
>>> x1 = np.random.random((npts, 3))
>>> x2 = np.random.random((npts, 3))

Find the projection of each vector in x1 onto a plane defined by
each vector in x2.

>>> x3 = project_onto_plane(x1, x2)
"""

n = normalized_vectors(x2)
d = elementwise_dot(x1,n)

return x1 - d[:,np.newaxis]*n