Generate block-diagonal rotation matrix about ax.

    [copy handedness from somewhere]

    .. todo:: Complete mtx_rot docstring

def mtx_rot(ax, theta, reps=1):
    """ Generate block-diagonal rotation matrix about ax.

    [copy handedness from somewhere]

    .. todo:: Complete mtx_rot docstring

    """

    # Imports
    import numpy as np
    from scipy import linalg as spla
    from ..const import PRM

    # Ensure |ax| is large enough for confident directionality
    if spla.norm(ax) < PRM.ZERO_VEC_TOL:
        raise ValueError("Norm of 'ax' is too small.")
    ## end if

    # Ensure ax is a normalized np.float64 3-vector
    ax = make_nd_vec(ax, nd=3, t=np.float64, norm=True)

    # Ensure reps is a positive scalar integer
    if not np.isscalar(reps):
        raise ValueError("'reps' must be scalar.")
    ## end if
    if not np.issubdtype(type(reps), int):
        raise ValueError("'reps' must be an integer.")
    ## end if
    if not reps > 0:
        raise ValueError("'reps' must be a positive integer.")
    ## end if

    # Ensure theta is scalar
    if not np.isscalar(theta):
        raise ValueError("'theta' must be scalar.")
    ## end if

    # Assemble the modified Levi-Civita matrix
    mod_lc = np.array([ [0, -ax[2], ax[1]],
                    [ax[2], 0, -ax[0]],
                    [-ax[1], ax[0], 0] ], dtype=np.float64)

    # Compute the outer product of the axis vector
    ax_oprod = np.dot(ax.reshape((3,1)), ax.reshape((1,3)))

    # Construct the base matrix
    #  Will need to refer to external math to explain this.
    base_mtx = np.add(
                        np.add( (1.0 - np.cos(theta)) * ax_oprod,
                                            np.cos(theta) * np.eye(3)
                              ),
                        np.sin(theta) * mod_lc
                     )

    # Construct the block-diagonal replicated reflection matrix
    rot_mtx= spla.block_diag(*[base_mtx for i in range(reps)])

    # Return the result
    return rot_mtx