[Yt-svn] yt: 2 new changesets

[hg at spacepope.org]

hg Repository: yt
details:   yt/rev/daa638047115
changeset: 3597:daa638047115
user:      Britton Smith <brittonsmith at gmail.com>
date:
Sun Dec 12 15:52:18 2010 -0500
description:
Additional fix to mpi_catarray, thanks to Matt.

hg Repository: yt
details:   yt/rev/cf81ea02f843
changeset: 3598:cf81ea02f843
user:      Britton Smith <brittonsmith at gmail.com>
date:
Sun Dec 12 15:53:54 2010 -0500
description:
Merged.

diffstat:

 yt/data_objects/derived_quantities.py |    1 +
 yt/utilities/math_utils.py            |  106 ++++++++++++++++++++++++++++++++++
 2 files changed, 107 insertions(+), 0 deletions(-)

diffs (124 lines):

diff -r 68ef4be28392 -r cf81ea02f843 yt/data_objects/derived_quantities.py
--- a/yt/data_objects/derived_quantities.py	Thu Dec 09 18:15:18 2010 -0500
+++ b/yt/data_objects/derived_quantities.py	Sun Dec 12 15:53:54 2010 -0500
@@ -176,6 +176,7 @@
     :param field: The field to average
     :param weight: The field to weight by
     """
+    print 'quantity', field, data.pf
     num = (data[field] * data[weight]).sum()
     den = data[weight].sum()
     return num, den
diff -r 68ef4be28392 -r cf81ea02f843 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py	Thu Dec 09 18:15:18 2010 -0500
+++ b/yt/utilities/math_utils.py	Sun Dec 12 15:53:54 2010 -0500
@@ -409,3 +409,109 @@
     P[:,2] = 0
     return na.sqrt((P * P).sum(axis=1))
     
+def ortho_find(vec1):
+    r"""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`.
+
+    .. math:: 
+
+       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:
+
+    .. math:: 
+
+       vec1 \cdot vec2 = x_1 + (z_1 z_2) = 0 
+
+       \rightarrow z_2 = -(x_1 / z_1)
+
+    Of course, this will fail if `z1` = 0, in which case, let's say use 
+    `z2` = 1 and `x2` = 0:
+
+    .. math::
+    
+       \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:
+
+    .. math::
+
+       \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])
+    """
+    vec1 = na.array(vec1, dtype=na.float64)
+    # Normalize
+    norm = na.sqrt(na.vdot(vec1, vec1))
+    if norm == 0:
+        raise ValueError("Zero vector used as input.")
+    vec1 /= norm
+    x1 = vec1[0]
+    y1 = vec1[1]
+    z1 = vec1[2]
+    if z1 != 0:
+        x2 = 1.0
+        y2 = 0.0
+        z2 = -(x1 / z1)
+        norm2 = (1.0 + z2 ** 2.0) ** (0.5)
+    elif y1 != 0:
+        x2 = 0.0
+        z2 = 1.0
+        y2 = -(z1 / y1)
+        norm2 = (1.0 + y2 ** 2.0) ** (0.5)
+    else:
+        y2 = 1.0
+        z2 = 0.0
+        x2 = -(y1 / x1)
+        norm2 = (1.0 + z2 ** 2.0) ** (0.5)
+    vec2 = na.array([x2,y2,z2])
+    vec2 /= norm2
+    vec3 = na.cross(vec1, vec2)
+    return vec1, vec2, vec3