yt.utilities.math_utils module¶

yt.utilities.math_utils.compute_cylindrical_radius(CoM, L, P, V)[source]¶

Compute the radius for some data around an axis in cylindrical coordinates.

This is primarily for halo computations. Given some data, this computes the cylindrical radius for each point. This is accomplished by converting the reference frame of the center of mass of the halo.

Parameters:

CoM (array) – The center of mass in 3D.
L (array) – The angular momentum vector.
P (array) – The positions of the data to be modified (i.e. particle or grid cell positions). The array should be Nx3.
V (array) – The velocities of the data to be modified (i.e. particle or grid cell velocities). The array should be Nx3.

Returns:

cyl_r – An array N elements long that gives the radial velocity for each datum (particle).

Return type:

array

Examples

>>> CoM = np.array([0, 0, 0])
>>> L = np.array([0, 0, 1])
>>> P = np.array([[1, 0, 0], [1, 1, 1], [0, 0, 1], [1, 1, 0]])
>>> V = np.array([[0, 1, 10], [-1, -1, -1], [1, 1, 1], [1, -1, -1]])
>>> cyl_r = compute_cylindrical_radius(CoM, L, P, V)
>>> cyl_r
array([ 1.        ,  1.41421356,  0.        ,  1.41421356])

yt.utilities.math_utils.compute_parallel_velocity(CoM, L, P, V)[source]¶

Computes the parallel velocity for some data around an axis.

This is primarily for halo computations. Given some data, this computes the velocity component along the angular momentum vector. This is accomplished by converting the reference frame of the center of mass of the halo.

Parameters:

CoM (array) – The center of mass in 3D.
L (array) – The angular momentum vector.
P (array) – The positions of the data to be modified (i.e. particle or grid cell positions). The array should be Nx3.
V (array) – The velocities of the data to be modified (i.e. particle or grid cell velocities). The array should be Nx3.

Returns:

v – An array N elements long that gives the parallel velocity for each datum (particle).

Return type:

array

Examples

>>> CoM = np.array([0, 0, 0])
>>> L = np.array([0, 0, 1])
>>> P = np.array([[1, 0, 0], [1, 1, 1], [0, 0, 1], [1, 1, 0]])
>>> V = np.array([[0, 1, 10], [-1, -1, -1], [1, 1, 1], [1, -1, -1]])
>>> paraV = compute_parallel_velocity(CoM, L, P, V)
>>> paraV
array([10, -1,  1, -1])

yt.utilities.math_utils.compute_radial_velocity(CoM, L, P, V)[source]¶

Computes the radial velocity for some data around an axis.

This is primarily for halo computations. Given some data, this computes the radial velocity component for the data. This is accomplished by converting the reference frame of the center of mass of the halo.

Parameters:

CoM (array) – The center of mass in 3D.
L (array) – The angular momentum vector.
P (array) – The positions of the data to be modified (i.e. particle or grid cell positions). The array should be Nx3.
V (array) – The velocities of the data to be modified (i.e. particle or grid cell velocities). The array should be Nx3.

Returns:

v – An array N elements long that gives the radial velocity for each datum (particle).

Return type:

array

Examples

>>> CoM = np.array([0, 0, 0])
>>> L = np.array([0, 0, 1])
>>> P = np.array([[1, 0, 0], [1, 1, 1], [0, 0, 1], [1, 1, 0]])
>>> V = np.array([[0, 1, 10], [-1, -1, -1], [1, 1, 1], [1, -1, -1]])
>>> radV = compute_radial_velocity(CoM, L, P, V)
>>> radV
array([ 1.        ,  1.41421356 ,  0.        ,  0.])

yt.utilities.math_utils.compute_rotational_velocity(CoM, L, P, V)[source]¶

Computes the rotational velocity for some data around an axis.

This is primarily for halo computations. Given some data, this computes the circular rotational velocity of each point (particle) in reference to the axis defined by the angular momentum vector. This is accomplished by converting the reference frame of the center of mass of the halo.

Parameters:

CoM (array) – The center of mass in 3D.
L (array) – The angular momentum vector.
P (array) – The positions of the data to be modified (i.e. particle or grid cell positions). The array should be Nx3.
V (array) – The velocities of the data to be modified (i.e. particle or grid cell velocities). The array should be Nx3.

Returns:

v – An array N elements long that gives the circular rotational velocity for each datum (particle).

Return type:

array

Examples

>>> CoM = np.array([0, 0, 0])
>>> L = np.array([0, 0, 1])
>>> P = np.array([[1, 0, 0], [1, 1, 1], [0, 0, 1], [1, 1, 0]])
>>> V = np.array([[0, 1, 10], [-1, -1, -1], [1, 1, 1], [1, -1, -1]])
>>> circV = compute_rotational_velocity(CoM, L, P, V)
>>> circV
array([ 1.        ,  0.        ,  0.        ,  1.41421356])

yt.utilities.math_utils.compute_stddev_image(buff2, buff)[source]¶

This function computes the standard deviation of a weighted projection. It defines the standard deviation as sigma^2 = <v^2>-<v>^2, where the brackets indicate averages (with the weight) and v is the field in question.

There is the possibility that if the projection at a particular location is of a constant or a single cell/particle, then <v^2> = <v>^2 and instead of getting zero one gets roundoff error that results in sigma^2 < 0, which is unphysical.

We handle this here by performing the subtraction and checking that any and all negative values of sigma^2 can be attributed to roundoff and thus be safely set to zero. We error out if this is not the case. There are ways of computing the standard deviation that are designed to avoid this catastrophic cancellation, but this would require rather substantial and invasive changes to the projection machinery so for the time being it is avoided.

Parameters:

buff2 (ImageArray) – The image of the weighted averge of the field squared
buff (ImageArray) – The image of the weighted averge of the field

yt.utilities.math_utils.euclidean_dist(a, b)[source]¶

Find the Euclidean distance between two points.

Parameters:

a (array or list) – Either an ndim long list of coordinates corresponding to a single point or an (ndim, npoints) list of coordinates for many points in space.
b (array or list) – Either an ndim long list of coordinates corresponding to a single point or an (ndim, npoints) list of coordinates for many points in space.

Examples

>>> a = [0.1, 0.1, 0.1]
>>> b = [0.9, 0, 9, 0.9]
>>> period = 1.0
>>> dist = euclidean_dist(a, b)
>>> dist
1.38564064606

yt.utilities.math_utils.get_cyl_r(coords, normal)[source]¶

yt.utilities.math_utils.get_cyl_r_component(vectors, theta, normal)[source]¶

yt.utilities.math_utils.get_cyl_theta(coords, normal)[source]¶

yt.utilities.math_utils.get_cyl_theta_component(vectors, theta, normal)[source]¶

yt.utilities.math_utils.get_cyl_z(coords, normal)[source]¶

yt.utilities.math_utils.get_cyl_z_component(vectors, normal)[source]¶

yt.utilities.math_utils.get_lookat_matrix(eye, center, up)[source]¶

Given the position of a camera, the point it is looking at, and an up-direction. Computes the lookat matrix that moves all vectors such that the camera is at the origin of the coordinate system, looking down the z-axis.

Parameters:

eye (array_like) – The position of the camera. Must be 3D.
center (array_like) – The location that the camera is looking at. Must be 3D.
up (array_like) – The direction that is considered up for the camera. Must be 3D.

Returns:

lookat_matrix – A new 4x4 2D array in homogeneous coordinates. This matrix moves all vectors in the same way required to move the camera to the origin of the coordinate system, with it pointing down the negative z-axis.

Return type:

ndarray

yt.utilities.math_utils.get_orthographic_matrix(maxr, aspect, z_near, z_far)[source]¶

Given a field of view in radians, an aspect ratio, and a near and far plane distance, this routine computes the transformation matrix corresponding to perspective projection using homogeneous coordinates.

Parameters:

maxr (scalar) – should be max(|x|, |y|)
aspect (scalar) – The aspect ratio of width / height for the projection.
z_near (scalar) – The distance of the near plane from the camera.
z_far (scalar) – The distance of the far plane from the camera.

Returns:

persp_matrix – A new 4x4 2D array. Represents a perspective transformation in homogeneous coordinates. Note that this matrix does not actually perform the projection. After multiplying a 4D vector of the form (x_0, y_0, z_0, 1.0), the point will be transformed to some (x_1, y_1, z_1, w). The final projection is applied by performing a divide by w, that is (x_1/w, y_1/w, z_1/w, w/w). The matrix uses a row-major ordering, rather than the column major ordering typically used by OpenGL.

Return type:

ndarray

Notes

The usage of 4D homogeneous coordinates is for OpenGL and GPU hardware that automatically performs the divide by w operation. See the following for more details about the OpenGL perspective matrices.

http://www.scratchapixel.com/lessons/3d-basic-rendering/perspective-and-orthographic-projection-matrix/orthographic-projection-matrix https://www.tomdalling.com/blog/modern-opengl/explaining-homogenous-coordinates-and-projective-geometry/ http://www.songho.ca/opengl/gl_projectionmatrix.html

yt.utilities.math_utils.get_perspective_matrix(fovy, aspect, z_near, z_far)[source]¶

Parameters:

fovy (scalar) – The angle in degrees of the field of view.
aspect (scalar) – The aspect ratio of width / height for the projection.
z_near (scalar) – The distance of the near plane from the camera.
z_far (scalar) – The distance of the far plane from the camera.

Returns:

Return type:

ndarray

Notes

The usage of 4D homogeneous coordinates is for OpenGL and GPU hardware that automatically performs the divide by w operation. See the following for more details about the OpenGL perspective matrices.

https://www.tomdalling.com/blog/modern-opengl/explaining-homogenous-coordinates-and-projective-geometry/ http://www.songho.ca/opengl/gl_projectionmatrix.html

yt.utilities.math_utils.get_rotation_matrix(theta, rot_vector)[source]¶

Given an angle theta and a 3D vector rot_vector, this routine computes the rotation matrix corresponding to rotating theta radians about rot_vector.

Parameters:

theta (scalar) – The angle in radians.
rot_vector (array_like) – The axis of rotation. Must be 3D.

Returns:

rot_matrix – A new 3x3 2D array. This is the representation of a rotation of theta radians about rot_vector in the simulation box coordinate frame

Return type:

ndarray

See also

ortho_find

Examples

>>> a = [0, 1, 0]
>>> theta = 0.785398163  # pi/4
>>> rot = mu.get_rotation_matrix(theta, a)
>>> rot
array([[ 0.70710678,  0.        ,  0.70710678],
       [ 0.        ,  1.        ,  0.        ],
       [-0.70710678,  0.        ,  0.70710678]])
>>> np.dot(rot, a)
array([ 0.,  1.,  0.])
# since a is an eigenvector by construction
>>> np.dot(rot, [1, 0, 0])
array([ 0.70710678,  0.        , -0.70710678])

yt.utilities.math_utils.get_scale_matrix(dx, dy, dz)[source]¶

Given a scaling factor for each coordinate, returns a matrix that corresponds to the given scaling amounts.

Parameters:

dx (scalar) – A scaling factor for the x-coordinate.
dy (scalar) – A scaling factor for the y-coordinate.
dz (scalar) – A scaling factor for the z-coordinate.

Returns:

scale_matrix – A new 4x4 2D array. Represents a scaling by dx, dy, and dz in each coordinate respectively.

Return type:

ndarray

yt.utilities.math_utils.get_sph_phi(coords, normal)[source]¶

yt.utilities.math_utils.get_sph_phi_component(vectors, phi, normal)[source]¶

yt.utilities.math_utils.get_sph_r(coords)[source]¶

yt.utilities.math_utils.get_sph_r_component(vectors, theta, phi, normal)[source]¶

yt.utilities.math_utils.get_sph_theta(coords, normal)[source]¶

yt.utilities.math_utils.get_sph_theta_component(vectors, theta, phi, normal)[source]¶

yt.utilities.math_utils.get_translate_matrix(dx, dy, dz)[source]¶

Given a movement amount for each coordinate, creates a translation matrix that moves the vector by each amount.

Parameters:

dx (scalar) – A translation amount for the x-coordinate
dy (scalar) – A translation amount for the y-coordinate
dz (scalar) – A translation amount for the z-coordinate

Returns:

trans_matrix – A new 4x4 2D array. Represents a translation by dx, dy and dz in each coordinate respectively.

Return type:

ndarray

yt.utilities.math_utils.modify_reference_frame(CoM, L, P=None, V=None)[source]¶

Rotates and translates data into a new reference frame to make calculations easier.

This is primarily useful for calculations of halo data. The data is translated into the center of mass frame. Next, it is rotated such that the angular momentum vector for the data is aligned with the z-axis. Put another way, if one calculates the angular momentum vector on the data that comes out of this function, it will always be along the positive z-axis. If the center of mass is re-calculated, it will be at the origin.

Parameters:

CoM (array) – The center of mass in 3D.
L (array) – The angular momentum vector.
Optional
--------
P (array) – The positions of the data to be modified (i.e. particle or grid cell positions). The array should be Nx3.
V (array) – The velocities of the data to be modified (i.e. particle or grid cell velocities). The array should be Nx3.

Returns:

L (array) – The angular momentum vector equal to [0, 0, 1] modulo machine error.
P (array) – The modified positional data. Only returned if P is not None
V (array) – The modified velocity data. Only returned if V is not None

Examples

>>> CoM = np.array([0.5, 0.5, 0.5])
>>> L = np.array([1, 0, 0])
>>> P = np.array([[1, 0.5, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0.5], [0, 0, 0]])
>>> V = p.copy()
>>> LL, PP, VV = modify_reference_frame(CoM, L, P, V)
>>> LL
array([  6.12323400e-17,   0.00000000e+00,   1.00000000e+00])
>>> PP
array([[  3.06161700e-17,   0.00000000e+00,   5.00000000e-01],
       [ -3.06161700e-17,   0.00000000e+00,  -5.00000000e-01],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  5.00000000e-01,  -5.00000000e-01,  -5.00000000e-01]])
>>> VV
array([[ -5.00000000e-01,   5.00000000e-01,   1.00000000e+00],
       [ -5.00000000e-01,   5.00000000e-01,   3.06161700e-17],
       [ -5.00000000e-01,   5.00000000e-01,   5.00000000e-01],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00]])

yt.utilities.math_utils.normalize_vector(vector)[source]¶

yt.utilities.math_utils.ortho_find(vec1)[source]¶

Find two complementary orthonormal vectors to a given vector.

For any given non-zero vector, there are infinite pairs of vectors orthonormal to it. This function gives you one arbitrary pair from that set along with the normalized version of the original vector.

Parameters:

vec1 (array_like) – An array or list to represent a 3-vector.

Returns:

vec1 (array) – The original 3-vector array after having been normalized.
vec2 (array) – A 3-vector array which is orthonormal to vec1.
vec3 (array) – A 3-vector array which is orthonormal to vec1 and vec2.

Raises:

ValueError – If input vector is the zero vector.

Notes

Our initial vector is vec1 which consists of 3 components: x1, y1, and z1. ortho_find determines a vector, vec2, which is orthonormal to vec1 by finding a vector which has a zero-value dot-product with vec1.

\[vec1 \cdot vec2 = x_1 x_2 + y_1 y_2 + z_1 z_2 = 0\]

As a starting point, we arbitrarily choose vec2 to have x2 = 1, y2 = 0:

\[ \begin{align}\begin{aligned}vec1 \cdot vec2 = x_1 + (z_1 z_2) = 0\\\rightarrow z_2 = -(x_1 / z_1)\end{aligned}\end{align} \]

Of course, this will fail if z1 = 0, in which case, let’s say use z2 = 1 and x2 = 0:

\[\rightarrow y_2 = -(z_1 / y_1)\]

Similarly, if y1 = 0, this case will fail, in which case we use y2 = 1 and z2 = 0:

\[\rightarrow x_2 = -(y_1 / x_1)\]

Since we don’t allow vec1 to be zero, all cases are accounted for.

Producing vec3, the complementary orthonormal vector to vec1 and vec2 is accomplished by simply taking the cross product of vec1 and vec2.

Examples

>>> a = [1.0, 2.0, 3.0]
>>> a, b, c = ortho_find(a)
>>> a
array([ 0.26726124,  0.53452248,  0.80178373])
>>> b
array([ 0.9486833 ,  0.        , -0.31622777])
>>> c
array([-0.16903085,  0.84515425, -0.50709255])

yt.utilities.math_utils.periodic_dist(a, b, period, periodicity=(True, True, True))[source]¶

Find the Euclidean periodic distance between two sets of points.

Parameters:

a (array or list) – Either an ndim long list of coordinates corresponding to a single point or an (ndim, npoints) list of coordinates for many points in space.
b (array of list) – Either an ndim long list of coordinates corresponding to a single point or an (ndim, npoints) list of coordinates for many points in space.
period (float or array or list) – If the volume is symmetrically periodic, this can be a single float, otherwise an array or list of floats giving the periodic size of the volume for each dimension.
periodicity (An ndim-element tuple of booleans) – If an entry is true, the domain is assumed to be periodic along that direction.

Examples

>>> a = [0.1, 0.1, 0.1]
>>> b = [0.9, 0, 9, 0.9]
>>> period = 1.0
>>> dist = periodic_dist(a, b, 1.0)
>>> dist
0.346410161514

yt.utilities.math_utils.periodic_position(pos, ds)[source]¶

Assuming periodicity, find the periodic position within the domain.

Parameters:

pos (array) – An array of floats.
ds (Dataset) – A simulation static output.

Examples

>>> a = np.array([1.1, 0.5, 0.5])
>>> data = {"Density": np.ones([32, 32, 32])}
>>> ds = load_uniform_grid(data, [32, 32, 32], 1.0)
>>> ppos = periodic_position(a, ds)
>>> ppos
array([ 0.1,  0.5,  0.5])

yt.utilities.math_utils.periodic_ray(start, end, left=None, right=None)[source]¶

Break up periodic ray into non-periodic segments. Accepts start and end points of periodic ray as YTArrays. Accepts optional left and right edges of periodic volume as YTArrays. Returns a list of lists of coordinates, where each element of the top-most list is a 2-list of start coords and end coords of the non-periodic ray:

[[[x0start,y0start,z0start], [x0end, y0end, z0end]],: [[x1start,y1start,z1start], [x1end, y1end, z1end]], …,]

Parameters:

start (array) – The starting coordinate of the ray.
end (array) – The ending coordinate of the ray.
left (optional, array) – The left corner of the periodic domain. If not given, an array of zeros with same size as the starting coordinate us used.
right (optional, array) – The right corner of the periodic domain. If not given, an array of ones with same size as the starting coordinate us used.

Examples

>>> import yt
>>> start = yt.YTArray([0.5, 0.5, 0.5])
>>> end = yt.YTArray([1.25, 1.25, 1.25])
>>> periodic_ray(start, end)
[
    [
        YTArray([0.5, 0.5, 0.5]) (dimensionless),
        YTArray([1., 1., 1.]) (dimensionless)
    ],
    [
        YTArray([0., 0., 0.]) (dimensionless),
        YTArray([0.25, 0.25, 0.25]) (dimensionless)
    ]
 ]

yt.utilities.math_utils.quartiles(a, axis=None, out=None, overwrite_input=False)[source]¶

Compute the quartile values (25% and 75%) along the specified axis in the same way that the numpy.median calculates the median (50%) value alone a specified axis. Check numpy.median for details, as it is virtually the same algorithm.

Returns an array of the quartiles of the array elements [lower quartile, upper quartile].

Parameters:

a (array_like) – Input array or object that can be converted to an array.
axis ({None, int}, optional) – Axis along which the quartiles are computed. The default (axis=None) is to compute the quartiles along a flattened version of the array.
out (ndarray, optional) – Alternative output array in which to place the result. It must have the same shape and buffer length as the expected output, but the type (of the output) will be cast if necessary.
overwrite_input ({False, True}, optional) – If True, then allow use of memory of input array (a) for calculations. The input array will be modified by the call to quartiles. This will save memory when you do not need to preserve the contents of the input array. Treat the input as undefined, but it will probably be fully or partially sorted. Default is False. Note that, if overwrite_input is True and the input is not already an ndarray, an error will be raised.

Returns:

quartiles – A new 2D array holding the result (unless out is specified, in which case that array is returned instead). If the input contains integers, or floats of smaller precision than 64, then the output data-type is float64. Otherwise, the output data-type is the same as that of the input.

Return type:

ndarray