# Source code for moldesign.mathutils

# Copyright 2016 Autodesk Inc.
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
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