def from_spherical_coords(theta_phi, phi=None):
    """Return the quaternion corresponding to these spherical coordinates

    Assumes the spherical coordinates correspond to the quaternion R via

        R = exp(phi*z/2) * exp(theta*y/2)

    The angles naturally must be in radians for this to make any sense.

    Note that this quaternion rotates `z` onto the point with the given
    spherical coordinates, but also rotates `x` and `y` onto the usual basis
    vectors (theta and phi, respectively) at that point.

    Parameters
    ----------
    theta_phi: float or array of floats
        This argument may either contain an array with last dimension of
        size 2, where those two elements describe the (theta, phi) values in
        radians for each point; or it may contain just the theta values in
        radians, in which case the next argument must also be given.
    phi: None, float, or array of floats
        If this array is given, it must be able to broadcast against the
        first argument.

    Returns
    -------
    R: quaternion array
        If the second argument is not given to this function, the shape
        will be the same as the input shape except for the last dimension,
        which will be removed.  If the second argument is given, this
        output array will have the shape resulting from broadcasting the
        two input arrays against each other.

    """
    # Figure out the input angles from either type of input
    if phi is None:
        theta_phi = np.asarray(theta_phi, dtype=np.double)
        theta = theta_phi[..., 0]
        phi  = theta_phi[..., 1]
    else:
        theta = np.asarray(theta_phi, dtype=np.double)
        phi = np.asarray(phi, dtype=np.double)

    # Set up the output array
    R = np.empty(np.broadcast(theta, phi).shape + (4,), dtype=np.double)

    # Compute the actual values of the quaternion components
    R[..., 0] =  np.cos(phi/2)*np.cos(theta/2)  # scalar quaternion components
    R[..., 1] = -np.sin(phi/2)*np.sin(theta/2)  # x quaternion components
    R[..., 2] =  np.cos(phi/2)*np.sin(theta/2)  # y quaternion components
    R[..., 3] =  np.sin(phi/2)*np.cos(theta/2)  # z quaternion components

    return as_quat_array(R)
Return the quaternion corresponding to these spherical coordinates

    Assumes the spherical coordinates correspond to the quaternion R via

        R = exp(phi*z/2) * exp(theta*y/2)

    The angles naturally must be in radians for this to make any sense.

    Note that this quaternion rotates `z` onto the point with the given
    spherical coordinates, but also rotates `x` and `y` onto the usual basis
    vectors (theta and phi, respectively) at that point.

    Parameters
    ----------
    theta_phi: float or array of floats
        This argument may either contain an array with last dimension of
        size 2, where those two elements describe the (theta, phi) values in
        radians for each point; or it may contain just the theta values in
        radians, in which case the next argument must also be given.
    phi: None, float, or array of floats
        If this array is given, it must be able to broadcast against the
        first argument.

    Returns
    -------
    R: quaternion array
        If the second argument is not given to this function, the shape
        will be the same as the input shape except for the last dimension,
        which will be removed.  If the second argument is given, this
        output array will have the shape resulting from broadcasting the
        two input arrays against each other.