Source code for moldesign.mathutils

# Copyright 2016 Autodesk Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import numpy as np

from moldesign import units as u


[docs]def perpendicular(vec):
    assert vec.shape == (3,)
    direction = normalized(vec)
    if abs(direction[2]) < 0.9:
        cross_axis = np.array([0.0, 0.0, 1.0])
    else:
        cross_axis = np.array([0.0, 1.0, 0.0])
    perp = normalized(np.cross(direction, cross_axis))
    return perp


[docs]def normalized(vec):
    """ Return a vector normalized in L2.
    If vector is 0, return 0

    Args:
        vec (u.Vector): vector to be normalized

    Returns:
        u.Vector: normalized vector
    """
    if len(vec.shape) == 1:  # it's just a single column vector
        mag = vec.dot(vec)
        if mag == 0.0:
            return vec*0.0
        else:
            return vec/np.sqrt(mag)
    else:  # treat as list of vectors
        mag = (vec*vec).sum(axis=1)
        mag[mag == 0.0] = 1.0  # prevent div by 0
        return vec / np.sqrt(mag)[:, None]

[docs]def safe_arccos(costheta):
    """ Version of arccos that can handle numerical noise greater than 1.0
    """
    if hasattr(costheta, 'shape') and costheta.shape:  # vector version
        assert (np.abs(costheta)-1.0 < 1.0e-13).all()
        costheta[costheta > 1.0] = 1.0
        costheta[costheta < -1.0] = -1.0
        return np.arccos(costheta)

    else:
        if abs(costheta) > 1.0:
            assert abs(costheta) - 1.0 < 1.0e-14
            return u.pi
        else:
            return np.arccos(costheta)


[docs]def sub_angles(a, b):
    """ Subtract two angles, keeping the result within [-180,180)
    """
    c = a - b
    return (c + 180.0 * u.degrees) % (360.0 * u.degrees) - (180.0 * u.degrees)


[docs]def apply_4x4_transform(trans, vecs):
    """
    Applies a 4x4 transformation vector so one or more 3-D position vector
    :param trans:
    :param vecs:
    :return: transformed position vector
    """
    has_units = False
    if hasattr(vecs, 'get_units'):
        has_units = True
        units = vecs.get_units()
        vecs = vecs.magnitude
    if len(vecs.shape) == 1:
        v = np.ones(4)
        v[:3] = vecs
        vt = trans.dot(v)
        result = vt[:3] / vt[3]
    else:
        v = np.ones((4, len(vecs)))
        v[:3, :] = vecs.T
        vt = trans.dot(v)
        result = (vt[:3] / vt[3]).T
    if has_units:
        result = result * units
    return result