Generate block-diagonal reflection matrix about nv.

    reps must be >=1 and indicates the number of times the reflection
    matrix should be repeated along the block diagonal.  Typically this
    will be the number of atoms in a geometry.

    .. todo:: Complete mtx_refl docstring

def mtx_refl(nv, reps=1):
    """ Generate block-diagonal reflection matrix about nv.

    reps must be >=1 and indicates the number of times the reflection
    matrix should be repeated along the block diagonal.  Typically this
    will be the number of atoms in a geometry.

    .. todo:: Complete mtx_refl docstring

    """

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

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

    # Ensure nv is a normalized np.float64 3-vector
    nv = make_nd_vec(nv, 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

    # Initialize the single-point reflection transform matrix
    base_mtx = np.zeros(shape=(3,3), dtype=np.float64)

    # Construct the single-point transform matrix
    for i in range(3):
        for j in range(i,3):
            if i==j:
                base_mtx[i,j] = 1 - 2*nv[i]**2
            else:
                base_mtx[i,j] = base_mtx[j,i] = -2*nv[i]*nv[j]
            ## end if
        ## next j
    ## next i

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

    # Return the result
    return refl_mtx