def inertia_tensor(geom, masses):
    """Generate the 3x3 moment-of-inertia tensor.

    Compute the 3x3 moment-of-inertia tensor for the
    provided geometry and atomic masses.  Always recenters the
    geometry to the center of mass as the first step.

    Reference for inertia tensor: [Kro92]_, Eq. (2.26)

    .. todo:: Replace cite eventually with link to exposition in user guide.

    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.

    Returns
    -------
    tensor
        3 x 3 |npfloat_| --
        Moment of inertia tensor for the system

    Raises
    ------
    ~exceptions.ValueError
        If shapes of `geom` & `masses` are inconsistent

    """

    # Imports
    import numpy as np

    # 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'
    geom = ctr_geom(geom, masses)

    # Expand the masses if required. Shape should only ever be (N,) or (3N,),
    #  else would raise an exception within the above 'ctr_geom' call
    if geom.shape[0] == 3*masses.shape[0]:
        masses = masses.repeat(3)
    ## end if

    # Initialize the tensor matrix
    tensor = np.zeros((3,3))

    # Fill the matrix
    for i in range(3):
        for j in range(i,3):
            if i == j:
                # On-diagonal element; calculate indices to include
                ind = np.concatenate([np.array(list(map(lambda v: v % 3,
                                        range(i+1, i+3)))) + o for o in
                                        range(0,geom.shape[0],3)])

                # Calculate the tensor element
                tensor[i,i] = np.multiply(np.square(geom[ind]),
                                                        masses[ind]).sum()
            else:
                # Off-diagonal element; calculate the indices
                ind_i = np.array(range(i,geom.shape[0]+i,3))
                ind_j = np.array(range(j,geom.shape[0]+j,3))

                # Calculate the tensor element and its symmetric partner
                tensor[i,j] = np.multiply(
                        np.sqrt(np.multiply(masses[ind_i], masses[ind_j])) ,
                        np.multiply(geom[ind_i], geom[ind_j]) ).sum() * -1
                tensor[j,i] = tensor[i,j]
            ## end if
        ## next j
    ## next i

    # Return the tensor
    return tensor