# Source code for sfs.util

"""Various utility functions."""

import numpy as np
from . import defs

[docs]def rotation_matrix(n1, n2):
"""Compute rotation matrix for rotation from n1 to n2.

Parameters
----------
n1, n2 : (3,) array_like
Two vectors.  They don't have to be normalized.

Returns
-------
(3, 3) numpy.ndarray
Rotation matrix.

"""
n1 = normal_vector(n1)
n2 = normal_vector(n2)
I = np.identity(3)
if np.all(n1 == n2):
return I  # no rotation
elif np.all(n1 == -n2):
return -I  # flip
# TODO: check for *very close to* parallel vectors

# Algorithm from http://math.stackexchange.com/a/476311
v = v0, v1, v2 = np.cross(n1, n2)
s = np.linalg.norm(v)  # sine
c = np.inner(n1, n2)  # cosine
vx = np.matrix([[0, -v2, v1],
[v2, 0, -v0],
[-v1, v0, 0]])  # skew-symmetric cross-product matrix
return I + vx + vx**2 * (1 - c) / s**2

[docs]def wavenumber(omega, c=None):
"""Compute the wavenumber for a given radial frequency."""
if c is None:
c = defs.c
return omega / c

[docs]def direction_vector(alpha, beta=np.pi/2):
"""Compute normal vector from azimuth, colatitude."""
return sph2cart(alpha, beta, 1)

[docs]def sph2cart(alpha, beta, r):
"""Spherical to cartesian coordinates."""
x = r * np.cos(alpha) * np.sin(beta)
y = r * np.sin(alpha) * np.sin(beta)
z = r * np.cos(beta)

return x, y, z

[docs]def cart2sph(x, y, z):
"""Cartesian to spherical coordinates."""
alpha = np.arctan2(y, x)
beta = np.arccos(z / np.sqrt(x**2 + y**2))
r = np.sqrt(x**2 + y**2 + z**2)

return alpha, beta, r

[docs]def asarray_1d(a, **kwargs):
"""Squeeze the input and check if the result is one-dimensional.

Returns a converted to a :class:numpy.ndarray and stripped of
all singleton dimensions.  Scalars are "upgraded" to 1D arrays.
The result must have exactly one dimension.
If not, an error is raised.

"""
result = np.squeeze(np.asarray(a, **kwargs))
if result.ndim == 0:
result = result.reshape((1,))
elif result.ndim > 1:
raise ValueError("array must be one-dimensional")
return result

[docs]def asarray_of_rows(a, **kwargs):
"""Convert to 2D array, turn column vector into row vector.

Returns a converted to a :class:numpy.ndarray and stripped of
all singleton dimensions.  If the result has exactly one dimension,
it is re-shaped into a 2D row vector.

"""
result = np.squeeze(np.asarray(a, **kwargs))
if result.ndim == 1:
result = result.reshape(1, -1)
return result

[docs]def as_xyz_components(components, **kwargs):
"""Convert components to :class:XyzComponents of NumPy arrays.

The components are first converted to NumPy arrays (using
:func:numpy.asarray) which are then assembled into an
:class:XyzComponents object.

Parameters
----------
components : triple or pair of array_like
The values to be used as X, Y and Z arrays.  Z is optional.
**kwargs
All further arguments are forwarded to
:func:numpy.asarray, which is applied to the elements of
components.

"""
return XyzComponents([np.asarray(c, **kwargs) for c in components])

[docs]def strict_arange(start, stop, step=1, endpoint=False, dtype=None, **kwargs):
"""Like :func:numpy.arange, but compensating numeric errors.

Unlike :func:numpy.arange, but similar to :func:numpy.linspace,
providing endpoint=True includes both endpoints.

Parameters
----------
start, stop, step, dtype
See :func:numpy.arange.
endpoint
See :func:numpy.linspace.

.. note:: With endpoint=True, the difference between start
and end value must be an integer multiple of the
corresponding spacing value!
**kwargs
All further arguments are forwarded to :func:numpy.isclose.

Returns
-------
numpy.ndarray
Array of evenly spaced values.  See :func:numpy.arange.

"""
remainder = (stop - start) % step
if np.any(np.isclose(remainder, (0.0, step), **kwargs)):
if endpoint:
stop += step * 0.5
else:
stop -= step * 0.5
elif endpoint:
raise ValueError("Invalid stop value for endpoint=True")
return np.arange(start, stop, step, dtype)

[docs]def xyz_grid(x, y, z, spacing, endpoint=True, **kwargs):
"""Create a grid with given range and spacing.

Parameters
----------
x, y, z : float or pair of float
Inclusive range of the respective coordinate or a single value
if only a slice along this dimension is needed.
spacing : float or triple of float
Grid spacing.  If a single value is specified, it is used for
all dimensions, if multiple values are given, one value is used
per dimension.  If a dimension (x, y or z) has only a
single value, the corresponding spacing is ignored.
endpoint : bool, optional
If True (the default), the endpoint of each range is
included in the grid.  Use False to get a result similar to
:func:numpy.arange.  See :func:sfs.util.strict_arange.
**kwargs
All further arguments are forwarded to
:func:sfs.util.strict_arange.

Returns
-------
XyzComponents
A grid that can be used for sound field calculations.
See :class:sfs.util.XyzComponents.

--------
strict_arange, numpy.meshgrid

"""
if np.isscalar(spacing):
spacing = [spacing] * 3
ranges = []
scalars = []
for i, coord in enumerate([x, y, z]):
if np.isscalar(coord):
scalars.append((i, coord))
else:
start, stop = coord
ranges.append(strict_arange(start, stop, spacing[i],
endpoint=endpoint, **kwargs))
grid = np.meshgrid(*ranges, sparse=True, copy=False)
for i, s in scalars:
grid.insert(i, s)
return XyzComponents(grid)

[docs]def normalize(p, grid, xnorm):
"""Normalize sound field wrt position xnorm."""
return p / np.abs(probe(p, grid, xnorm))

[docs]def probe(p, grid, x):
"""Determine the value at position x in the sound field p."""
grid = as_xyz_components(grid)
x = asarray_1d(x)
r = np.linalg.norm(grid - x)
idx = np.unravel_index(r.argmin(), r.shape)
return p[idx]

"""Broadcast arguments to the same shape and then use :func:zip."""

[docs]def normal_vector(x):
"""Normalize a 1D vector."""
x = asarray_1d(x)
return x / np.linalg.norm(x)

[docs]def displacement(v, omega):
"""Particle displacement

.. math::

d(x, t) = \int_0^t v(x, t) dt

"""
return as_xyz_components(v) / (1j * omega)

[docs]def db(x, power=False):
"""Convert x to decibel.

Parameters
----------
x : array_like
Input data.  Values of 0 lead to negative infinity.
power : bool, optional
If power=False (the default), x is squared before
conversion.

"""
with np.errstate(divide='ignore'):
return 10 if power else 20 * np.log10(np.abs(x))

[docs]class XyzComponents(np.ndarray):
"""See __init__()."""

def __init__(self, components):
"""Triple (or pair) of components: x, y, and optionally z.

Instances of this class can be used to store coordinate grids
(either regular grids like in :func:sfs.util.xyz_grid or
arbitrary point clouds) or vector fields (e.g. particle
velocity).

This class is a subclass of :class:numpy.ndarray.  It is
one-dimensional (like a plain :class:list) and has a length of
3 (or 2, if no z-component is available).  It uses
dtype=object in order to be able to store other
numpy.ndarray\s of arbitrary shapes but also scalars, if
needed.  Because it is a NumPy array subclass, it can be used in
operations with scalars and "normal" NumPy arrays, as long as
they have a compatible shape.  Like any NumPy array, instances
of this class are iterable and can be used, e.g., in for-loops
incompatible shape, a plain numpy.ndarray with dtype=object
is returned.

To make sure the components are NumPy arrays themselves, use
:func:sfs.util.as_xyz_components.

Parameters
----------
components : (3,) or (2,) array_like
The values to be used as X, Y and Z data.  Z is optional.

"""
# This method does nothing, it's only here for the documentation!

def __new__(cls, components):
# object arrays cannot be created and populated in a single step:
obj = np.ndarray.__new__(cls, len(components), dtype=object)
for i, component in enumerate(components):
obj[i] = component
return obj

def __array_finalize__(self, obj):
if self.ndim == 0:
pass  # this is allowed, e.g. for np.inner()
elif self.ndim > 1 or len(self) not in (2, 3):
raise ValueError("XyzComponents can only have 2 or 3 components")

def __array_prepare__(self, obj, context=None):
if obj.ndim == 1 and len(obj) in (2, 3):
return obj.view(XyzComponents)
return obj

def __array_wrap__(self, obj, context=None):
if obj.ndim != 1 or len(obj) not in (2, 3):
return obj.view(np.ndarray)
return obj

def __getitem__(self, index):
if isinstance(index, slice):
start, stop, step = index.indices(len(self))
if start == 0 and stop in (2, 3) and step == 1:
return np.ndarray.__getitem__(self, index)
# Slices other than xy and xyz are "downgraded" to ndarray
return np.ndarray.__getitem__(self.view(np.ndarray), index)

def __repr__(self):
return 'XyzComponents(\n' + ',\n'.join(
'    {0}={1}'.format(name, repr(data).replace('\n', '\n      '))
for name, data in zip('xyz', self)) + ')'

def make_property(index, doc):

def getter(self):
return self[index]

def setter(self, value):
self[index] = value

return property(getter, setter, doc=doc)

x = make_property(0, doc='x-component.')
y = make_property(1, doc='y-component.')
z = make_property(2, doc='z-component (optional).')

del make_property

[docs]    def apply(self, func, *args, **kwargs):
"""Apply a function to each component.

The function func will be called once for each component,
passing the current component as first argument.  All further
arguments are passed after that.
The results are returned as a new :class:XyzComponents object.

"""
return XyzComponents([func(i, *args, **kwargs) for i in self])