Principal axes and moments of inertia for the indicated geometry.

    Calculated by :func:`scipy.linalg.eigh`, since the moment of inertia tensor
    is symmetric (real-Hermitian) by construction.  More convenient to
    compute both the axes and moments at the same time since the eigenvectors
    must be processed to ensure repeatable results.

    The principal axes (inertia tensor eigenvectors) are processed in a
    fashion to ensure repeatable, **identical** generation, including
    orientation AND directionality.

    .. todo:: Add ref to exposition in webdocs once written up.

    Parameters
    ----------
    geom
        length-3N |npfloat_| --
        Coordinates of the atoms

    masses
        length-N OR length-3N |npfloat_| --
        Atomic masses of the atoms. Length-3N option is to allow calculation of
        a per-coordinate perturbed value.

    on_tol
        |npfloat_|,  optional --
        Tolerance for deviation from unity/zero for principal axis dot products
        within which axes are considered orthonormal. Default is
        :data:`opan.const.DEF.ORTHONORM_TOL`.

    Returns
    -------
    moments
        length-3 |npfloat_| --
        Principal inertial moments, sorted in increasing order
        :math:`\\left(0 \\leq I_A \\leq I_B \\leq I_C\\right)`

    axes
        3 x 3 |npfloat_| --
        Principal axes, as column vectors, sorted with the principal moments
        and processed for repeatability. The axis corresponding to
        ``moments[i]`` is retrieved as ``axes[:,i]``

    top
        :class:`~opan.const.EnumTopType` --
        Detected molecular top type

def principals(geom, masses, on_tol=_DEF.ORTHONORM_TOL):
    """Principal axes and moments of inertia for the indicated geometry.

    Calculated by :func:`scipy.linalg.eigh`, since the moment of inertia tensor
    is symmetric (real-Hermitian) by construction.  More convenient to
    compute both the axes and moments at the same time since the eigenvectors
    must be processed to ensure repeatable results.

    The principal axes (inertia tensor eigenvectors) are processed in a
    fashion to ensure repeatable, **identical** generation, including
    orientation AND directionality.

    .. todo:: Add ref to exposition in webdocs once written up.

    Parameters
    ----------
    geom
        length-3N |npfloat_| --
        Coordinates of the atoms

    masses
        length-N OR length-3N |npfloat_| --
        Atomic masses of the atoms. Length-3N option is to allow calculation of
        a per-coordinate perturbed value.

    on_tol
        |npfloat_|,  optional --
        Tolerance for deviation from unity/zero for principal axis dot products
        within which axes are considered orthonormal. Default is
        :data:`opan.const.DEF.ORTHONORM_TOL`.

    Returns
    -------
    moments
        length-3 |npfloat_| --
        Principal inertial moments, sorted in increasing order
        :math:`\\left(0 \\leq I_A \\leq I_B \\leq I_C\\right)`

    axes
        3 x 3 |npfloat_| --
        Principal axes, as column vectors, sorted with the principal moments
        and processed for repeatability. The axis corresponding to
        ``moments[i]`` is retrieved as ``axes[:,i]``

    top
        :class:`~opan.const.EnumTopType` --
        Detected molecular top type

    """

    # Imports
    import numpy as np
    from scipy import linalg as spla
    from ..const import PRM, EnumTopType as ETT
    from ..error import InertiaError, VectorError
    from .vector import rej, parallel_check as prlchk
    from .vector import orthonorm_check as orthchk

    # Center the geometry. Takes care of any improper shapes of geom or
    #  masses via the internal call to 'ctr_mass' within the call to
    #  'ctr_geom'.  Will need the centered geometry eventually anyways.
    geom = ctr_geom(geom, masses)

    # Get the inertia tensor
    tensor = inertia_tensor(geom, masses)

    # Orthogonalize and store eigenvalues/-vectors. eigh documentation says it
    #  will return ordered eigenvalues.... Store eigenvalues directly to
    #  the return variable; since eigenvectors probably need work, store them
    #  to a holding variable.
    moments, vecs = spla.eigh(tensor)

    # 'fail' init for 'top
    top = None

    # Detect top type; start with error check
    if moments[0] < -PRM.ZERO_MOMENT_TOL:  # pragma: no cover
        # Invalid moment; raise error. Should be impossible!
        raise InertiaError(InertiaError.NEG_MOMENT,
                    "Negative principal inertial moment", "")
    elif moments[0] < PRM.ZERO_MOMENT_TOL:
        # Zero first moment. Check whether others are too
        if all(moments < PRM.ZERO_MOMENT_TOL):
            top = ETT.ATOM
        else:
            top = ETT.LINEAR
        ## end if
    else:
        if abs((moments[1] / moments[0]) - 1.0) < PRM.EQUAL_MOMENT_TOL:
            # Spherical or oblate symmetrical
            if abs((moments[2] / moments[1]) - 1.0) < PRM.EQUAL_MOMENT_TOL:
                top = ETT.SPHERICAL
            else:
                top = ETT.SYMM_OBL
            ## end if
        else:
            # Prolate symmetrical or Asymmetric
            if abs((moments[2] / moments[1]) - 1.0) < PRM.EQUAL_MOMENT_TOL:
                top = ETT.SYMM_PROL
            else:
                top = ETT.ASYMM
            ## end if
        ## end if
    ## end if

    # Check for nothing assigned (this should never occur!)
    if top is None:  # pragma: no cover
        raise InertiaError(InertiaError.TOP_TYPE,
                    "Unrecognized molecular top type","")
    ## end if

    # Initialize the axes
    axes = np.zeros((3,3))

    # Define the axes depending on the top type
    if top == ETT.ATOM:
        # Just use the coordinate axes
        axes = np.identity(3, dtype=np.float_)

    elif top == ETT.LINEAR:
        # Zero-moment (molecular) axis always pointed toward the first atom,
        #  or the second if the first is at center-of-mass
        if spla.norm(geom[0:3]) >= PRM.ZERO_VEC_TOL:
            axes[:,0] = geom[0:3] / spla.norm(geom[0:3])
        else:  # pragma: no cover (assume alt case ok)
            axes[:,0] = geom[3:6] / spla.norm(geom[3:6])
        ## end if

        # Second axis is the normalized rejection of the x-axis on the first
        #  axis, unless the molecule lies along the x-axis in which case it
        #  is taken as the normalized rejection of the y-axis on the first vec.
        if prlchk(axes[:,0], np.array([1.,0.,0.])):
            # Too nearly (anti-)parallel
            axes[:,1] = rej(np.array([0.,1.,0.]), axes[:,0])
        else:  # pragma: no cover (assume alt case ok)
            # Sufficiently non-(anti-)parallel
            axes[:,1] = rej(np.array([1.,0.,0.]), axes[:,0])
        ## end if
        axes[:,1] /= spla.norm(axes[:,1])

        # Third axis is the first crossed with the second
        axes[:,2] = np.cross(axes[:,0], axes[:,1])

    elif top == ETT.ASYMM:
        # Vectors should already be orthonormal; following error should
        #  never occur
        if not orthchk(vecs, tol=on_tol):  # pragma: no cover
            raise VectorError(VectorError.ORTHONORM,
                         "'eigh' produced non-orthonormal axes", "")
        ## end if

        # Duplicate the vectors to the axes object
        axes = vecs.copy()

        # Orient first two axes to have positive dot products with their
        #  respective first non-zero, non-orthogonal atomic displacements.
        #  Possibly fragile to some sort of highly unusual geometry.
        axes[:,0] *= np.sign(np.dot(vecs[:,0], _fadn_orth(vecs[:,0], geom)))
        axes[:,1] *= np.sign(np.dot(vecs[:,1], _fadn_orth(vecs[:,1], geom)))

        # Orient the third axis such that a3 = a1 {cross} a2
        axes[:,2] *= np.sign(np.dot(axes[:,2], np.cross(axes[:,0], axes[:,1])))

    elif top == ETT.SYMM_OBL:
        # First axis is taken as the normalized rejection of the first
        #  non-(anti)parallel atomic displacement onto the third eigenvector.
        axes[:,0] = rej(_fadn_par(vecs[:,2], geom), vecs[:,2])
        axes[:,0] /= spla.norm(axes[:,0])

        # Try to take the third axis directionality as that giving a positive
        #  dot product with the first non-orthogonal atomic displacement.
        # A planar system will cause an error in the _fadn_orth call, that
        #  is trapped here.
        # If planar, take the third axis as the normalized cross product
        #  of the first and second nonzero atomic displacements.

        try:
            axes[:,2] = vecs[:,2] * np.sign(np.dot(vecs[:,2],
                                                _fadn_orth(vecs[:,2], geom)))
        except InertiaError as IE:
            # Check that typecode is as expected for error from planar system.
            if not IE.tc == InertiaError.BAD_GEOM:  # pragma: no cover
                raise
            ## end if

            # Select the appropriate displacements to define the third axis
            if spla.norm(geom[0:3]) < PRM.ZERO_VEC_TOL:
                # First displacement is zero
                axes[:,2] = np.cross(geom[3:6], geom[6:9]) # pragma: no cover
            elif spla.norm(geom[3:6]) < PRM.ZERO_VEC_TOL:
                # Second displacement is zero
                axes[:,2] = np.cross(geom[0:3], geom[6:9]) # pragma: no cover
            else:
                # First and second displacements are okay
                axes[:,2] = np.cross(geom[0:3], geom[3:6])
            ## end if

            # Regardless of which calculation, normalize the vector
        finally:
            axes[:,2] /= spla.norm(axes[:,2])
        ## end try

        # Second axis is the third axis crossed with the first
        axes[:,1] = np.cross(axes[:,2], axes[:,0])

    elif top == ETT.SYMM_PROL:
        # Special case of prolate symmetric is linear, which is separately
        #  detected and already addressed.
        # First (non-degenerate) axis is just taken to have a positive dot
        #  product with its first non-orthogonal nonzero displacement
        axes[:,0] = vecs[:,0] * np.sign(np.dot(vecs[:,0],
                                                    _fadn_orth(vecs[:,0], geom)))

        # Second (first degenerate) axis is the normalized rejection of the
        #  first non-parallel displacement onto the first axis
        axes[:,1] = rej(_fadn_par(axes[:,0], geom), axes[:,0])
        axes[:,1] /= spla.norm(axes[:,1])

        # Third (second degenerate) axis is just the first axis crossed with
        #  the second.
        axes[:,2] = np.cross(axes[:,0], axes[:,1])

    elif top == ETT.SPHERICAL:
        # No preferred orientation -- ALL vectors are degenerate axes
        # First axis is the first nonzero displacement, normalized
        axes[:,0] = geom[3:6] if spla.norm(geom[0:3]) < PRM.ZERO_VEC_TOL \
                                                                else geom[0:3]
        axes[:,0] /= spla.norm(axes[:,0])

        # Second axis is the normalized rejection onto the first axis of the
        #  first nonzero non-parallel displacement from that first axis
        axes[:,1] = rej(_fadn_par(axes[:,0], geom), axes[:,0])
        axes[:,1] /= spla.norm(axes[:,1])

        # Third axis is the first crossed with the second
        axes[:,2] = np.cross(axes[:,0], axes[:,1])

    ## end if

    # Reconfirm orthonormality. Again, the error should never occur.
    if not orthchk(axes, tol=on_tol): # pragma: no cover
        raise VectorError(VectorError.ORTHONORM,
                    "Axis conditioning broke orthonormality","")
    ## end if

    # Return the moments, axes, and top type
    return moments, axes, top