# Path: attitude/geom/util.py
# def dot(*matrices):
#     return reduce(N.dot, matrices)
# 
# def vector(*args):
#     """
#     A single vector in any coordinate basis,
#     as a numpy array.
#     """
#     return N.array(args)
# 
# def augment(vec):
#     """
#     Augment a vector in any orthonormal basis
#     with a trailing 1 to form a homogeneous
#     coordinate vector in that coordinate system.
#     """
#     return N.append(vec,[1])
# 
# class Plane(N.ndarray):
#     def hessian_normal(plane):
#         """
#         Return the Hessian Normal form of a plane
#         (ax + by + cz + d = 0) where [a,b,c] forms
#         the unit normal vector and d is the distance
#         to the origin."""
#         return plane/N.linalg.norm(plane[:-1])
# 
#     def offset(plane):
#         """
#         Returns the offset of the plane from the
#         origin or an arbitrary point.
#         """
#         v = plane.hessian_normal()
#         return -N.array(v[:-1]*v[-1])
# 
#     def normal(plane):
#         v = plane.hessian_normal()
#         return v[:-1]
# 
# def angle(v1,v2, cos=False):
#     """
#     Find the angle between two vectors.
# 
#     :param cos: If True, the cosine of the
#     angle will be returned. False by default.
#     """
#     n = (norm(v1)*norm(v2))
#     _ = dot(v1,v2)/n
#     return _ if cos else N.arccos(_)
# 
# def perpendicular_vector(n):
#     """
#     Get a random vector perpendicular
#     to the given vector
#     """
#     dim = len(n)
#     if dim == 2:
#         return n[::-1]
#     # More complex in 3d
#     for ix in range(dim):
#         _ = N.zeros(dim)
#         # Try to keep axes near the global projection
#         # by finding vectors perpendicular to higher-
#         # index axes first. This may or may not be worth
#         # doing.
#         _[dim-ix-1] = 1
#         v1 = N.cross(n,_)
#         if N.linalg.norm(v1) != 0:
#             return v1
#     raise ValueError("Cannot find perpendicular vector")


import numpy as N
from .util import (dot, vector, augment, Plane,
                   angle, perpendicular_vector)
        ax = list(axes)
        #ax[-1] *= -1  # Not sure what is going on here...
        arr = N.diag(ax + [-1])
        return arr.view(cls)

    @classmethod
    def from_semiaxes(cls,axes):
        """
        Get axis-aligned elliptical conic from axis lenths
        This can be converted into a hyperbola by getting the dual conic
        """
        ax = list(1/N.array(axes)**2)
        #ax[-1] *= -1  # Not sure what is going on here...
        arr = N.diag(ax + [-1])
        return arr.view(cls)

    def center(conic):
        # (https://en.wikipedia.org/wiki/Matrix_representation_of_conic_sections#Center)
        ec = N.linalg.inv(conic[:-1,:-1])
        eo = -conic[:-1,-1]
        return dot(ec,eo.T)

    def contains(ell, p, shell_only=False):
        """
        Check to see whether point is inside
        conic.

        :param exact: Only solutions exactly on conic
          are considered (default: False).
        """

# Path: attitude/geom/util.py
# def dot(*matrices):
#     return reduce(N.dot, matrices)
# 
# def vector(*args):
#     """
#     A single vector in any coordinate basis,
#     as a numpy array.
#     """
#     return N.array(args)
# 
# def augment(vec):
#     """
#     Augment a vector in any orthonormal basis
#     with a trailing 1 to form a homogeneous
#     coordinate vector in that coordinate system.
#     """
#     return N.append(vec,[1])
# 
# class Plane(N.ndarray):
#     def hessian_normal(plane):
#         """
#         Return the Hessian Normal form of a plane
#         (ax + by + cz + d = 0) where [a,b,c] forms
#         the unit normal vector and d is the distance
#         to the origin."""
#         return plane/N.linalg.norm(plane[:-1])
# 
#     def offset(plane):
#         """
#         Returns the offset of the plane from the
#         origin or an arbitrary point.
#         """
#         v = plane.hessian_normal()
#         return -N.array(v[:-1]*v[-1])
# 
#     def normal(plane):
#         v = plane.hessian_normal()
#         return v[:-1]
# 
# def angle(v1,v2, cos=False):
#     """
#     Find the angle between two vectors.
# 
#     :param cos: If True, the cosine of the
#     angle will be returned. False by default.
#     """
#     n = (norm(v1)*norm(v2))
#     _ = dot(v1,v2)/n
#     return _ if cos else N.arccos(_)
# 
# def perpendicular_vector(n):
#     """
#     Get a random vector perpendicular
#     to the given vector
#     """
#     dim = len(n)
#     if dim == 2:
#         return n[::-1]
#     # More complex in 3d
#     for ix in range(dim):
#         _ = N.zeros(dim)
#         # Try to keep axes near the global projection
#         # by finding vectors perpendicular to higher-
#         # index axes first. This may or may not be worth
#         # doing.
#         _[dim-ix-1] = 1
#         v1 = N.cross(n,_)
#         if N.linalg.norm(v1) != 0:
#             return v1
#     raise ValueError("Cannot find perpendicular vector")


# Path: attitude/geom/conics.py
import numpy as N
from .util import (dot, vector, augment, Plane,
                   angle, perpendicular_vector)

        ax = list(axes)
        #ax[-1] *= -1  # Not sure what is going on here...
        arr = N.diag(ax + [-1])
        return arr.view(cls)

    @classmethod
    def from_semiaxes(cls,axes):
        """
        Get axis-aligned elliptical conic from axis lenths
        This can be converted into a hyperbola by getting the dual conic
        """
        ax = list(1/N.array(axes)**2)
        #ax[-1] *= -1  # Not sure what is going on here...
        arr = N.diag(ax + [-1])
        return arr.view(cls)

    def center(conic):
        # (https://en.wikipedia.org/wiki/Matrix_representation_of_conic_sections#Center)
        ec = N.linalg.inv(conic[:-1,:-1])
        eo = -conic[:-1,-1]
        return dot(ec,eo.T)

    def contains(ell, p, shell_only=False):
        """
        Check to see whether point is inside
        conic.

        :param exact: Only solutions exactly on conic
          are considered (default: False).
        """

        v = augment(p)