Python numpy.arccos() Examples

def vector_angle(u, v, direction=None):
    '''
    vector_angle(u, v) yields the angle between the two vectors u and v. The optional argument 
    direction is by default None, which specifies that the smallest possible angle between the
    vectors be reported; if the vectors u and v are 2D vectors and direction parameters True and
    False specify the clockwise or counter-clockwise directions, respectively; if the vectors are
    3D vectors, then direction may be a 3D point that is not in the plane containing u, v, and the
    origin, and it specifies around which direction (u x v or v x u) the the counter-clockwise angle
    from u to v should be reported (the cross product vector that has a positive dot product with
    the direction argument is used as the rotation axis).
    '''
    if direction is None:
        return np.arccos(vector_angle_cos(u, v))
    elif direction is True:
        return np.arctan2(v[1], v[0]) - np.arctan2(u[1], u[0])
    elif direction is False:
        return np.arctan2(u[1], u[0]) - np.arctan2(v[1], v[0])
    else:
        axis1 = normalize(u)
        axis2 = normalize(np.cross(u, v))
        if np.dot(axis2, direction) < 0:
            axis2 = -axis2
        return np.arctan2(np.dot(axis2, v), np.dot(axis1, v)) 

def arccosine(x, null=(-np.inf, np.inf)):
    '''
    arccosine(x) is equivalent to acos(x) except that it also works on sparse arrays.

    The optional argument null (default, (-numpy.inf, numpy.inf)) may be specified to indicate what
    value(s) should be assigned when x < -1 or x > 1. If only one number is given, then it is used
    for both values; otherwise the first value corresponds to <-1 and the second to >1.  If null is
    None, then an error is raised when invalid values are encountered.
    '''
    if sps.issparse(x): x = x.toarray()
    else:               x = np.asarray(x)
    try:    (nln,nlp) = null
    except Exception: (nln,nlp) = (null,null)
    ii = None if nln is None else np.where(x < -1)
    jj = None if nlp is None else np.where(x > 1)
    if ii: x[ii] = 0
    if jj: x[jj] = 0
    x = np.arccos(x)
    if ii: x[ii] = nln
    if jj: x[jj] = nlp
    return x 

def rotate_camera_to_point_at(up_from, lookat_from, up_to, lookat_to):
  inputs = [up_from, lookat_from, up_to, lookat_to]
  for i in range(4):
    inputs[i] = normalize(np.array(inputs[i]).reshape((-1,)))
  up_from, lookat_from, up_to, lookat_to = inputs
  r1 = r_between(lookat_from, lookat_to)

  new_x = np.dot(r1, np.array([1, 0, 0]).reshape((-1, 1))).reshape((-1))
  to_x = normalize(np.cross(lookat_to, up_to))
  angle = np.arccos(np.dot(new_x, to_x))
  if angle > ANGLE_EPS:
    if angle < np.pi - ANGLE_EPS:
      ax = normalize(np.cross(new_x, to_x))
      flip = np.dot(lookat_to, ax)
      if flip > 0:
        r2 = get_r_matrix(lookat_to, angle)
      elif flip < 0:
        r2 = get_r_matrix(lookat_to, -1. * angle)
    else:
      # Angle of rotation is too close to 180 degrees, direction of rotation
      # does not matter.
      r2 = get_r_matrix(lookat_to, angle)
  else:
    r2 = np.eye(3)
  return np.dot(r2, r1) 

def subspace_disagreement_measure(self, Ps, Pt, Pst):
        """
        Get the best value for the number of subspaces
        For more details, read section 3.4 of the paper.
        **Parameters**
          Ps: Source subspace
          Pt: Target subspace
          Pst: Source + Target subspace
        """

        def compute_angles(A, B):
            _, S, _ = np.linalg.svd(np.dot(A.T, B))
            S[np.where(np.isclose(S, 1, atol=self.eps) == True)[0]] = 1
            return np.arccos(S)

        max_d = min(Ps.shape[1], Pt.shape[1], Pst.shape[1])
        alpha_d = compute_angles(Ps, Pst)
        beta_d = compute_angles(Pt, Pst)
        d = 0.5 * (np.sin(alpha_d) + np.sin(beta_d))
        return np.argmax(d) 

def inverse_method(self,N,d):
        
        t = np.linspace(1e-3,0.999,N)
        f = np.log( t / (1 - t) )
        f = f/f[0]
        
        psi= np.pi*f
        cosPsi = np.cos(psi)
        sinTheta = ( np.abs(cosPsi) + (1-np.abs(cosPsi))*np.random.rand(len(cosPsi)))
        
        theta = np.arcsin(sinTheta)
        theta = np.pi-theta + (2*theta - np.pi)*np.round(np.random.rand(len(t)))
        cosPhi = cosPsi/sinTheta
        phi = np.arccos(cosPhi)*(-1)**np.round(np.random.rand(len(t)))
        
        coords = SkyCoord(phi*u.rad,(np.pi/2-theta)*u.rad,d*np.ones(len(phi))*u.pc)

        return coords 

def make12(b):
    theta1 = numpy.arccos(1/numpy.sqrt(5))
    theta2 = (numpy.pi - theta1) * .5
    r = b/2./numpy.sin(theta1/2)
    rot72 = rotmatz(numpy.pi*2/5)
    s1 = numpy.sin(theta1)
    c1 = numpy.cos(theta1)
    p1 = numpy.array(( s1*r,  0,  c1*r))
    p2 = numpy.array((-s1*r,  0, -c1*r))
    coord = [(  0,  0,    r)]
    for i in range(5):
        coord.append(p1)
        p1 = numpy.dot(p1, rot72)
    for i in range(5):
        coord.append(p2)
        p2 = numpy.dot(p2, rot72)
    coord.append((  0,  0,  -r))
    return numpy.array(coord) 

def get_uv(self, xyz_vec):
        # Extract lens parameters of interest.
        fov_rad = self.lens.fov_deg * pi / 180
        fov_scale = np.float32(2 * self.lens.radius_px / fov_rad)
        # Normalize the input vector and rotate to match lens reference axes.
        xyz_rot = get_rotation_matrix(self.lens.center_qq) * matrix_norm(xyz_vec)
        # Convert to polar coordinates relative to lens boresight.
        # (In lens coordinates, unit vector's X axis gives boresight angle;
        #  normalize Y/Z to get a planar unit vector for the bearing.)
        # Note: Image +Y maps to 3D +Y, and image +X maps to 3D +Z.
        theta_rad = np.arccos(xyz_rot[0,:])
        proj_vec = matrix_norm(np.concatenate((xyz_rot[2,:], xyz_rot[1,:])))
        # Fisheye lens maps 3D angle to focal-plane radius.
        # TODO: Do we need a better model for lens distortion?
        rad_px = theta_rad * fov_scale
        # Convert back to focal-plane rectangular coordinates.
        uv = np.multiply(rad_px, proj_vec) + self.lens.center_px
        return np.asarray(uv + 0.5, dtype=int)

    # Given an 2xN array of UV pixel coordinates, check if each pixel is
    # within the fisheye field of view. Returns N-element boolean mask. 

def test_both_limit_angle_inactive(self, optimization_variables_avg):
        l, Q, A = optimization_variables_avg
        Q_in = 5e1 * np.eye(len(l))
        max_angle = 45
        m = 2e-3
        m1 = 4e-3
        f = .02
        x = optimization_methods.optimize_focality(l, Q, f,
                                                   max_el_current=m,
                                                   max_total_current=m1,
                                                   Qin=Q_in,
                                                   max_angle=max_angle)
        x_sp = optimize_focality(
            l, Q, f, max_el_current=m, max_total_current=m1,
            Qin=Q_in, max_angle=max_angle)

        assert np.linalg.norm(x, 1) <= 2 * m1 + 1e-4
        assert np.all(np.abs(x) <= m + 1e-4)
        assert np.isclose(np.sum(x), 0)
        assert np.isclose(l.dot(x), f)
        assert np.arccos(l.dot(x) / np.sqrt(x.dot(Q_in).dot(x))) <= np.deg2rad(max_angle)
        assert np.allclose(x.dot(Q.dot(x)), x_sp.dot(Q.dot(x_sp)), rtol=1e-4, atol=1e-5) 

def test_both_limit_angle_limit(self, optimization_variables_avg):
        l, Q, A = optimization_variables_avg
        Q_in = 5e1 * np.eye(len(l))
        max_angle = 25
        m = 2e-3
        m1 = 4e-3
        f = .02
        x = optimization_methods.optimize_focality(l, Q, f,
                                                   max_el_current=m,
                                                   max_total_current=m1,
                                                   Qin=Q_in,
                                                   max_angle=max_angle)
        x_sp = optimize_focality(
            l, Q, f, max_el_current=m, max_total_current=m1,
            Qin=Q_in, max_angle=max_angle)
        assert np.linalg.norm(x, 1) <= 2 * m1 + 1e-4
        assert np.all(np.abs(x) <= m + 1e-4)
        assert np.isclose(np.sum(x), 0)
        assert np.isclose(l.dot(x), f)
        assert np.arccos(l.dot(x) / np.sqrt(x.dot(Q_in).dot(x))) <= np.deg2rad(max_angle)
        assert np.allclose(x.dot(Q.dot(x)), x_sp.dot(Q.dot(x_sp)), rtol=1e-4, atol=1e-5) 

def test_both_limit_angle_Q_iteration(self, optimization_variables_avg_QCQP):
        l, Q, A, Q_in = optimization_variables_avg_QCQP
        # l, Q, A = optimization_variables_avg
        # Q_in = 6 * np.eye(len(l)) + np.outer(l, l)
        max_angle = 20
        m = 2e-3
        m1 = 4e-3
        f = .01
        x = optimization_methods.optimize_focality(l, Q, f,
                                                   max_el_current=m,
                                                   max_total_current=m1,
                                                   Qin=Q_in,
                                                   max_angle=max_angle)
        x_sp = optimize_focality(
            l, Q, f, max_el_current=m, max_total_current=m1,
            Qin=Q_in, max_angle=max_angle)
        assert np.linalg.norm(x, 1) <= 2 * m1 + 1e-4
        assert np.all(np.abs(x) <= m + 1e-4)
        assert np.isclose(np.sum(x), 0)
        assert np.isclose(l.dot(x), f)
        assert np.arccos(l.dot(x) / np.sqrt(x.dot(Q_in).dot(x))) <= np.deg2rad(max_angle)
        assert np.allclose(x.dot(Q.dot(x)), x_sp.dot(Q.dot(x_sp)), rtol=1e-4, atol=1e-5) 

def test_both_limit_angle_infeasible_field(self, optimization_variables_avg_QCQP):
        l, Q, A, Q_in = optimization_variables_avg_QCQP
        max_angle = 15
        m = 2e-3
        m1 = 4e-3
        f = 2
        x = optimization_methods.optimize_focality(l, Q, f,
                                                   max_el_current=m,
                                                   max_total_current=m1,
                                                   Qin=Q_in,
                                                   max_angle=max_angle)

        assert np.linalg.norm(x, 1) <= 2 * m1 + 1e-4
        assert np.all(np.abs(x) <= m + 1e-4)
        assert np.isclose(np.sum(x), 0)
        assert np.arccos(l.dot(x) / np.sqrt(x.dot(Q_in).dot(x))) <= np.deg2rad(max_angle) 

def test_limit_nr_angle(seld, optimization_variables_avg_QCQP):
        l, Q, A, Q_in = optimization_variables_avg_QCQP
        max_angle = 15
        m = 2e-3
        m1 = 4e-3
        f = 2
        n = 4
        x = optimization_methods.optimize_focality(l, Q, f,
                                                   max_el_current=m,
                                                   max_total_current=m1,
                                                   Qin=Q_in,
                                                   max_angle=max_angle,
                                                   max_active_electrodes=n)

        assert np.linalg.norm(x, 1) <= 2 * m1 + 1e-4
        assert np.all(np.abs(x) <= m + 1e-4)
        assert np.isclose(np.sum(x), 0)
        assert np.linalg.norm(x, 0) <= n
        assert np.arccos(l.dot(x) / np.sqrt(x.dot(Q_in).dot(x))) <= np.deg2rad(max_angle) 

def test_limit_nr_angle_change_Q(seld, optimization_variables_avg_QCQP):
        l, Q, A, Q_in = optimization_variables_avg_QCQP
        max_angle = 15
        m = 2e-3
        m1 = 4e-3
        f = .01
        n = 4
        x = optimization_methods.optimize_focality(l, Q, f,
                                                   max_el_current=m,
                                                   max_total_current=m1,
                                                   Qin=Q_in,
                                                   max_angle=max_angle,
                                                   max_active_electrodes=n)

        assert np.linalg.norm(x, 1) <= 2 * m1 + 1e-4
        assert np.all(np.abs(x) <= m + 1e-4)
        assert np.isclose(np.sum(x), 0)
        assert np.linalg.norm(x, 0) <= n
        assert np.arccos(l.dot(x) / np.sqrt(x.dot(Q_in).dot(x))) <= np.deg2rad(max_angle) 

def model_model_angles_btwn_axes(A, B, mxBasis):  # Note: default 'gm' basis
    """ Angle between the rotation axes of A and B (1-qubit gates)"""
    decomp = _tools.decompose_gate_matrix(A)
    decomp2 = _tools.decompose_gate_matrix(B)
    axisOfRotn = decomp.get('axis of rotation', None)
    rotnAngle = decomp.get('pi rotations', 'X')
    axisOfRotn2 = decomp2.get('axis of rotation', None)
    rotnAngle2 = decomp2.get('pi rotations', 'X')

    if rotnAngle == 'X' or abs(rotnAngle) < 1e-4 or \
       rotnAngle2 == 'X' or abs(rotnAngle2) < 1e-4:
        return _np.nan

    if axisOfRotn is None or axisOfRotn2 is None:
        return _np.nan

    real_dot = _np.clip(_np.real(_np.dot(axisOfRotn, axisOfRotn2)), -1.0, 1.0)
    return _np.arccos(abs(real_dot)) / _np.pi
    #Note: abs() allows axis to be off by 180 degrees -- if showing *angle* as
    #      well, must flip sign of angle of rotation if you allow axis to
    #      "reverse" by 180 degrees. 

def axisangle_from_rotm(R):
  # logarithm of rotation matrix
  # R = R.reshape(-1,3,3)
  # tr = np.trace(R, axis1=1, axis2=2)
  # phi = np.arccos(np.clip((tr - 1) / 2, -1, 1))
  # scale = np.zeros_like(phi)
  # div = 2 * np.sin(phi)
  # np.divide(phi, div, out=scale, where=np.abs(div) > 1e-6)
  # A = (R - R.transpose(0,2,1)) * scale.reshape(-1,1,1)
  # aa = np.stack((A[:,2,1], A[:,0,2], A[:,1,0]), axis=1)
  # return aa.squeeze()
  R = R.reshape(-1,3,3)
  omega = np.empty((R.shape[0], 3), dtype=R.dtype)
  omega[:,0] = R[:,2,1] - R[:,1,2]
  omega[:,1] = R[:,0,2] - R[:,2,0]
  omega[:,2] = R[:,1,0] - R[:,0,1]
  r = np.linalg.norm(omega, axis=1).reshape(-1,1)
  t = np.trace(R, axis1=1, axis2=2).reshape(-1,1)
  omega = np.arctan2(r, t-1) * omega
  aa = np.zeros_like(omega)
  np.divide(omega, r, out=aa, where=r != 0)
  return aa.squeeze() 

def quat_slerp_space(q0, q1, num=100, endpoint=True):
  q0 = q0.ravel()
  q1 = q1.ravel()
  dot = q0.dot(q1)
  if dot < 0:
    q1 *= -1
    dot *= -1
  t = np.linspace(0, 1, num=num, endpoint=endpoint, dtype=q0.dtype)
  t = t.reshape((-1,1))
  if dot > 0.9995:
    ret = q0 + t * (q1 - q0)
    return ret
  dot = np.clip(dot, -1, 1)
  theta0 = np.arccos(dot)
  theta = theta0 * t
  s0 = np.cos(theta) - dot * np.sin(theta) / np.sin(theta0)
  s1 = np.sin(theta) / np.sin(theta0)
  return (s0 * q0) + (s1 * q1) 

def test_branch_cuts(self):
        # check branch cuts and continuity on them
        _check_branch_cut(np.log,   -0.5, 1j, 1, -1, True)
        _check_branch_cut(np.log2,  -0.5, 1j, 1, -1, True)
        _check_branch_cut(np.log10, -0.5, 1j, 1, -1, True)
        _check_branch_cut(np.log1p, -1.5, 1j, 1, -1, True)
        _check_branch_cut(np.sqrt,  -0.5, 1j, 1, -1, True)

        _check_branch_cut(np.arcsin, [ -2, 2],   [1j, 1j], 1, -1, True)
        _check_branch_cut(np.arccos, [ -2, 2],   [1j, 1j], 1, -1, True)
        _check_branch_cut(np.arctan, [0-2j, 2j],  [1,  1], -1, 1, True)

        _check_branch_cut(np.arcsinh, [0-2j,  2j], [1,   1], -1, 1, True)
        _check_branch_cut(np.arccosh, [ -1, 0.5], [1j,  1j], 1, -1, True)
        _check_branch_cut(np.arctanh, [ -2,   2], [1j, 1j], 1, -1, True)

        # check against bogus branch cuts: assert continuity between quadrants
        _check_branch_cut(np.arcsin, [0-2j, 2j], [ 1,  1], 1, 1)
        _check_branch_cut(np.arccos, [0-2j, 2j], [ 1,  1], 1, 1)
        _check_branch_cut(np.arctan, [ -2,  2], [1j, 1j], 1, 1)

        _check_branch_cut(np.arcsinh, [ -2,  2, 0], [1j, 1j, 1], 1, 1)
        _check_branch_cut(np.arccosh, [0-2j, 2j, 2], [1,  1,  1j], 1, 1)
        _check_branch_cut(np.arctanh, [0-2j, 2j, 0], [1,  1,  1j], 1, 1) 

def test_branch_cuts_complex64(self):
        # check branch cuts and continuity on them
        _check_branch_cut(np.log,   -0.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.log2,  -0.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.log10, -0.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.log1p, -1.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.sqrt,  -0.5, 1j, 1, -1, True, np.complex64)

        _check_branch_cut(np.arcsin, [ -2, 2],   [1j, 1j], 1, -1, True, np.complex64)
        _check_branch_cut(np.arccos, [ -2, 2],   [1j, 1j], 1, -1, True, np.complex64)
        _check_branch_cut(np.arctan, [0-2j, 2j],  [1,  1], -1, 1, True, np.complex64)

        _check_branch_cut(np.arcsinh, [0-2j,  2j], [1,   1], -1, 1, True, np.complex64)
        _check_branch_cut(np.arccosh, [ -1, 0.5], [1j,  1j], 1, -1, True, np.complex64)
        _check_branch_cut(np.arctanh, [ -2,   2], [1j, 1j], 1, -1, True, np.complex64)

        # check against bogus branch cuts: assert continuity between quadrants
        _check_branch_cut(np.arcsin, [0-2j, 2j], [ 1,  1], 1, 1, False, np.complex64)
        _check_branch_cut(np.arccos, [0-2j, 2j], [ 1,  1], 1, 1, False, np.complex64)
        _check_branch_cut(np.arctan, [ -2,  2], [1j, 1j], 1, 1, False, np.complex64)

        _check_branch_cut(np.arcsinh, [ -2,  2, 0], [1j, 1j, 1], 1, 1, False, np.complex64)
        _check_branch_cut(np.arccosh, [0-2j, 2j, 2], [1,  1,  1j], 1, 1, False, np.complex64)
        _check_branch_cut(np.arctanh, [0-2j, 2j, 0], [1,  1,  1j], 1, 1, False, np.complex64) 

def test_against_cmath(self):
        import cmath

        points = [-1-1j, -1+1j, +1-1j, +1+1j]
        name_map = {'arcsin': 'asin', 'arccos': 'acos', 'arctan': 'atan',
                    'arcsinh': 'asinh', 'arccosh': 'acosh', 'arctanh': 'atanh'}
        atol = 4*np.finfo(complex).eps
        for func in self.funcs:
            fname = func.__name__.split('.')[-1]
            cname = name_map.get(fname, fname)
            try:
                cfunc = getattr(cmath, cname)
            except AttributeError:
                continue
            for p in points:
                a = complex(func(np.complex_(p)))
                b = cfunc(p)
                assert_(abs(a - b) < atol, "%s %s: %s; cmath: %s" % (fname, p, a, b)) 

def interp_z(z0, z1, ratio, interp='linear'):
    if interp == 'linear':
        z_t = (1 - ratio) * z0 + ratio * z1

    if interp == 'slerp':
        N = len(z0)
        z_t = []
        for i in range(N):
            z0_i = z0[i]
            z1_i = z1[i]
            z0_n = z0_i / np.linalg.norm(z0_i)
            z1_n = z1_i / np.linalg.norm(z1_i)
            omega = np.arccos(np.dot(z0_n, z1_n))
            sin_omega = np.sin(omega)
            if sin_omega == 0:
                z_i = interp_z(z0_i, z1_i, ratio, 'linear')
            else:
                z_i = np.sin((1 - ratio) * omega) / sin_omega * z0_i + np.sin(ratio * omega) / sin_omega * z1_i
            z_t.append(z_i[np.newaxis, ...])
        z_t = np.concatenate(z_t, axis=0)
    return z_t 

def polar_distance(x1, x2):
    """
    Given two arrays of numbers x1 and x2, pairs the cells that are the
    closest and provides the pairing matrix index: x1(index(1,:)) should be as
    close as possible to x2(index(2,:)). The function outputs the average of 
    the absolute value of the differences abs(x1(index(1,:))-x2(index(2,:))).
    :param x1: vector 1
    :param x2: vector 2
    :return: d: minimum distance between d
             index: the permutation matrix
    """
    x1 = np.reshape(x1, (1, -1), order='F')
    x2 = np.reshape(x2, (1, -1), order='F')
    N1 = x1.size
    N2 = x2.size
    diffmat = np.arccos(np.cos(x1 - np.reshape(x2, (-1, 1), order='F')))
    min_N1_N2 = np.min([N1, N2])
    index = np.zeros((min_N1_N2, 2), dtype=int)
    if min_N1_N2 > 1:
        for k in range(min_N1_N2):
            d2 = np.min(diffmat, axis=0)
            index2 = np.argmin(diffmat, axis=0)
            index1 = np.argmin(d2)
            index2 = index2[index1]
            index[k, :] = [index1, index2]
            diffmat[index2, :] = float('inf')
            diffmat[:, index1] = float('inf')
        d = np.mean(np.arccos(np.cos(x1[:, index[:, 0]] - x2[:, index[:, 1]])))
    else:
        d = np.min(diffmat)
        index = np.argmin(diffmat)
        if N1 == 1:
            index = np.array([1, index])
        else:
            index = np.array([index, 1])
    return d, index 

def angle(node_i,node_j,x):
        v=np.array([node_j.X-node_i.X,node_j.Y-node_i.Y,node_j.Z-node_i.Z])
        L1=np.sqrt(v.dot(v))
        L2=np.sqrt(x.dot(x))
        return np.arccos(v.dot(x)/L1/L2)

        #derivation 

def get_orientation_matrix(arrayAxis, alignToAxis):
    """
    Get the rotation matrix that aligns arrayAxis to alignToAxis

    :Parameters:
        #. arrayAxis (list, tuple, numpy.ndarray): xyzArray axis.
        #. alignToAxis (list, tuple, numpy.ndarray): The axis to align to.
    """
    # normalize alignToAxis
    alignToAxisNorm = np.linalg.norm(alignToAxis)
    assert alignToAxisNorm>0, LOGGER.error("alignToAxis returned 0 norm")
    alignToAxis = np.array(alignToAxis, dtype=FLOAT_TYPE)/alignToAxisNorm
    # normalize arrayAxis
    arrayAxisNorm = np.linalg.norm(arrayAxis)
    assert arrayAxisNorm>0, LOGGER.error("arrayAxis returned 0 norm")
    arrayAxis = np.array(arrayAxis, dtype=FLOAT_TYPE)/arrayAxisNorm
    # calculate rotationAngle
    dotProduct = np.dot(arrayAxis, alignToAxis)
    if np.abs(dotProduct-1) <= PRECISION :
        rotationAngle = 0
    elif np.abs(dotProduct+1) <= PRECISION :
        rotationAngle = PI
    else:
        rotationAngle = np.arccos( dotProduct )
    if np.isnan(rotationAngle) or np.abs(rotationAngle) <= PRECISION :
        return np.array([[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]]).astype(FLOAT_TYPE)
    # calculate rotation axis.
    if np.abs(rotationAngle-PI) <= PRECISION:
        rotationAxis = get_random_perpendicular_vector(arrayAxis)
    else:
        rotationAxis = np.cross(alignToAxis, arrayAxis)
    #rotationAxis /= np.linalg.norm(rotationAxis)
    # calculate rotation matrix
    return get_rotation_matrix(rotationAxis, rotationAngle) 

def orient(xyzArray, arrayAxis, alignToAxis):
    """
    Rotates xyzArray using the rotation matrix that rotates and aligns
    arrayAxis to alignToAXis.

    :Parameters:
        #. xyzArray (numpy.ndarray): The xyz (N,3) array to rotate.
        #. arrayAxis (list, tuple, numpy.ndarray): xyzArray axis.
        #. alignToAxis (list, tuple, numpy.ndarray): The axis to align to.
    """
    # normalize alignToAxis
    alignToAxisNorm = np.linalg.norm(alignToAxis)
    assert alignToAxisNorm>0, LOGGER.error("alignToAxis returned 0 norm")
    alignToAxis = np.array(alignToAxis, dtype=FLOAT_TYPE)/alignToAxisNorm
    # normalize arrayAxis
    arrayAxisNorm = np.linalg.norm(arrayAxis)
    assert arrayAxisNorm>0, LOGGER.error("arrayAxis returned 0 norm")
    arrayAxis = np.array(arrayAxis, dtype=FLOAT_TYPE)/arrayAxisNorm
    # calculate rotationAngle
    dotProduct = np.dot(arrayAxis, alignToAxis)
    if np.abs(dotProduct-1) <= PRECISION :
        rotationAngle = 0
    elif np.abs(dotProduct+1) <= PRECISION :
        rotationAngle = PI
    else:
        rotationAngle = np.arccos( dotProduct )
    if np.isnan(rotationAngle) or np.abs(rotationAngle) <= PRECISION :
        return xyzArray
    # calculate rotation axis.
    if np.abs(rotationAngle-PI) <= PRECISION:
        rotationAxis = get_random_perpendicular_vector(arrayAxis)
    else:
        rotationAxis = np.cross(alignToAxis, arrayAxis)
    #rotationAxis /= np.linalg.norm(rotationAxis)
    # calculate rotation matrix
    rotationMatrix = get_rotation_matrix(rotationAxis, rotationAngle)
    # rotate and return
    return rotate(xyzArray , rotationMatrix) 

def r_between(v_from_, v_to_):
  v_from = normalize(v_from_)
  v_to = normalize(v_to_)
  ax = normalize(np.cross(v_from, v_to))
  angle = np.arccos(np.dot(v_from, v_to))
  return get_r_matrix(ax, angle) 

def angle_diff(true_q, pred_q):
  angles = 2 * (
      180.0 /
      np.pi) * np.arccos(np.abs(np.sum(np.multiply(pred_q, true_q), axis=1)))
  return angles 

def compute_angle(points1, points2, points3, *, degrees: bool = False) -> np.ndarray:
    """
    Computes the angle (p1, p2 [vertex], p3) between the provided points on a per-row basis.

    Parameters
    ----------
    points1 : np.ndarray
        The first list of points, can be 1D or 2D
    points2 : np.ndarray
        The second list of points, can be 1D or 2D
    points3 : np.ndarray
        The third list of points, can be 1D or 2D
    degrees : bool, options
        Returns the angle in degrees rather than radians if True

    Returns
    -------
    angles : np.ndarray
        The angle between the three points in radians

    Notes
    -----
    Units are not considered inside these expressions, please preconvert to the same units before using.
    """

    points1 = np.atleast_2d(points1)
    points2 = np.atleast_2d(points2)
    points3 = np.atleast_2d(points3)

    v12 = points1 - points2
    v23 = points2 - points3

    denom = _norm(v12) * _norm(v23)
    cosine_angle = np.einsum("ij,ij->i", v12, v23) / denom

    angle = np.pi - np.arccos(cosine_angle)

    if degrees:
        return np.degrees(angle)
    else:
        return angle 

def principal_angles(self, Ps, Pt):
        """
        Compute the principal angles between source (:math:`P_s`) and target (:math:`P_t`) subspaces in a Grassman which is defined as the following:

        :math:`d^{2}(P_s, P_t) = \sum_{i}( \theta_i^{2} )`,

        """
        # S = cos(theta_1, theta_2, ..., theta_n)
        _, S, _ = np.linalg.svd(np.dot(Ps.T, Pt))
        thetas_squared = np.arccos(S) ** 2

        return np.sum(thetas_squared) 

def test_set_planet_phase(self):
        """
        Test that set_planet_phase places planets at the correct phase angle
        
        Because some implementations depend on a specific planet population,
        there needs to be additional logic in the setup
        """
        whitelist = ['KnownRVPlanetsUniverse']
        
        for mod in self.allmods:
            if mod.__name__ in whitelist:
                continue
            with RedirectStreams(stdout=self.dev_null):
                spec = copy.deepcopy(self.spec)
                spec['modules']['PlanetPhysicalModel']='FortneyMarleyCahoyMix1'
                spec['modules']['StarCatalog']='EXOCAT1'
                if 'Kepler' in mod.__name__:
                    spec['modules']['PlanetPopulation']='KeplerLike1'
                elif 'SAG13' in mod.__name__:
                    spec['modules']['PlanetPopulation']='SAG13'
                    spec['Rprange'] = [1,10]
                elif 'DulzPlavchan' in mod.__name__:
                    spec['modules']['PlanetPopulation'] = 'DulzPlavchan'
                    
                obj = mod(**spec)
                
            # attempt to set planet phase to pi/4
            obj.set_planet_phase(np.pi/4.)
            betas = np.arccos(obj.r[:,2]/obj.d)
            val1 = np.abs(betas.to('rad').value - np.pi/4.)
            val2 = np.abs(betas.to('rad').value - np.pi/2.)
            inds1 = np.where(val1 < 1e-4)[0]
            inds2 = np.where(val2 < 1e-4)[0]
            num = len(inds1) + len(inds2)
                        
            self.assertTrue(num == obj.nPlans,"Phase angles not set correctly") 

def initialXYZpoints(num_pts=30):
    """ Quick and unprecise way of distributing points on a sphere
    """
    indices = arange(0, num_pts, dtype=float) + 0.5
    phi = arccos(1 - 2*indices/num_pts)
    theta = pi * (1 + 5**0.5) * indices
    x, y, z = cos(theta) * sin(phi), sin(theta) * sin(phi), cos(phi)
    v = np.asarray([[x[i], y[i], z[i]] for i in np.arange(len(x))]) # an array of each point on the sphere
    d = np.linalg.norm(v,axis=1) # used to ensure the length of each vector is 1
    return x, y, z, v