diff --git a/misc/dispersion2D-SVD.py b/misc/dispersion2D-SVD.py
index b76c4cc..8bc93f9 100755
--- a/misc/dispersion2D-SVD.py
+++ b/misc/dispersion2D-SVD.py
@@ -139,21 +139,6 @@ nelem = len(my)
 if pargs.lMax: #force commandline specified lMax
     TMatrices_orig = TMatrices_orig[...,0:nelem,:,0:nelem]
 
-ž = np.arange(2*nelem)
-tž = ž // nelem
-mž = my[ž%nelem]
-nž = ny[ž%nelem]
-TEž = ž[(mž+nž+tž) % 2 == 0]
-TMž = ž[(mž+nž+tž) % 2 == 1]
-
-č = np.arange(2*2*nelem)
-žč = č % (2* nelem)
-tč = tž[žč]
-mč = mž[žč]
-nč = nž[žč]
-TEč = č[(mč+nč+tč) % 2 == 0]
-TMč = č[(mč+nč+tč) % 2 == 1]
-
 TMatrices = np.array(np.broadcast_to(TMatrices_orig[:,nx,:,:,:,:],(len(freqs_orig),2,2,nelem,2,nelem)) )
 
 #TMatrices[:,:,:,:,:,ny==3] *= factor13inc
diff --git a/qpms/hexpoints.py b/qpms/hexpoints.py
index d53f7b7..03aba3a 100644
--- a/qpms/hexpoints.py
+++ b/qpms/hexpoints.py
@@ -436,8 +436,11 @@ def hexlattice_get_AB(lMax, k_hexside, maxlayer, circular=True, return_points =
     return d
 
 from scipy.constants import c
+from .timetrack import _time_b, _time_e
+from .qpms_p import symz_indexarrays
 
-def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, omega, klist, gaussianSigma=False, onlyNmin=0):
+def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, omega, klist, gaussianSigma=False, onlyNmin=0, verbose=False):
+    btime = _time_b(verbose)
     nelem = lMax * (lMax + 2)
     n2id = np.identity(2*nelem)
     n2id.shape = (2,nelem,2,nelem)
@@ -455,9 +458,11 @@ def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, ome
     d2u_translations = tdic['d2u_tr']*hexside*_s3
     
     if gaussianSigma:
+        sbtime = _time_b(verbose, step='Calculating gaussian envelope')
         unitcell_envelope = np.exp(-np.sum(tdic['self_tr']**2,axis=-1)/(2*gaussianSigma**2))
         u2d_envelope = np.exp(-np.sum(tdic['u2d_tr']**2,axis=-1)/(2*gaussianSigma**2))
         d2u_envelope = np.exp(-np.sum(tdic['d2u_tr']**2,axis=-1)/(2*gaussianSigma**2))
+        _time_e(sbtime, verbose, step='Calculating gaussian envelope')
     
     #TMatrices_om = TMatrices_interp(omega)
     if(not onlyNmin):
@@ -473,7 +478,8 @@ def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, ome
 
     leftmatrixlist = np.full((klist.shape[0],2,2,nelem,2,2,nelem),np.nan,dtype=complex)
     isNaNlist = np.zeros((klist.shape[0]), dtype=bool)
-    
+
+    sbtime = _time_b(verbose, step='Initializing matrices for SVD for a given list of k\'s.')
     # sem nějaká rozumná smyčka
     for ki in range(klist.shape[0]):
         k = klist[ki]
@@ -516,21 +522,21 @@ def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, ome
 
         leftmatrixlist[ki] = leftmatrix
 
-
     nnlist = np.logical_not(isNaNlist)
     leftmatrixlist_s = np.reshape(leftmatrixlist,(klist.shape[0], 2*2*nelem,2*2*nelem))[nnlist]
+    TEč, TMč = symz_indexarrays(lMax, 2)
     leftmatrixlist_TE = leftmatrixlist_s[np.ix_(np.arange(leftmatrixlist_s.shape[0]),TEč,TEč)]
     leftmatrixlist_TM = leftmatrixlist_s[np.ix_(np.arange(leftmatrixlist_s.shape[0]),TMč,TMč)]
 
+    _time_e(sbtime, verbose, step='Initializing matrices for SVD for a given list of k\'s.')
+    sbtime = _time_b(verbose, step='Calculating SVDs for a given list of k\'s.')
     if(not onlyNmin):
         svUfullTElist[nnlist], svSfullTElist[nnlist], svVfullTElist[nnlist] = np.linalg.svd(leftmatrixlist_TE, compute_uv=True)
         svUfullTMlist[nnlist], svSfullTMlist[nnlist], svVfullTMlist[nnlist] = np.linalg.svd(leftmatrixlist_TM, compute_uv=True)
+        _time_e(sbtime, verbose, step='Calculating SVDs for a given list of k\'s.')
         return ((svUfullTElist, svSfullTElist, svVfullTElist), (svUfullTMlist, svSfullTMlist, svVfullTMlist))
     else:
         minsvTElist[nnlist] = np.linalg.svd(leftmatrixlist_TE, compute_uv=False)[...,-onlyNmin:]
         minsvTMlist[nnlist] = np.linalg.svd(leftmatrixlist_TM, compute_uv=False)[...,-onlyNmin:]
-
-
-
-
-    
+        _time_e(sbtime, verbose, step='Calculating SVDs for a given list of k\'s.')
+        return (minsvTElist, minsvTMlist)
diff --git a/qpms/qpms_p.py b/qpms/qpms_p.py
index 714c8de..d852a6c 100644
--- a/qpms/qpms_p.py
+++ b/qpms/qpms_p.py
@@ -930,3 +930,41 @@ def apply_ndmatrix_left(matrix,tensor,axes):
     matrix = np.tensordot(matrix, tensor, axes=([-N+axn for axn in range(N)],axes))
     matrix = np.moveaxis(matrix, range(N), axes)
     return matrix
+
+
+def symz_indexarrays(lMax, npart = 1):
+    """
+    Returns indices that are used for separating the in-plane E ('TE' in the photonic crystal
+    jargon) and perpendicular E ('TM' in the photonic crystal jargon) modes
+    in the z-mirror symmetric systems.
+
+    Parameters
+    ----------
+    lMax : int
+        The maximum degree cutoff for the T-matrix to which these indices will be applied.
+
+    npart : int
+        Number of particles (TODO better description)
+
+    Returns
+    -------
+    TEč, TMč : (npart * 2 * nelem)-shaped bool ndarray
+        Mask arrays corresponding to the 'TE' and 'TM' modes, respectively.
+    """
+    nelem = lMax * (lMax + 2)
+    ž = np.arange(2*nelem) # single particle spherical wave indices
+    tž = ž // nelem # tž == 0: electric waves, tž == 1: magnetic waves
+    mž = my[ž%nelem]
+    nž = ny[ž%nelem]
+    TEž = ž[(mž+nž+tž) % 2 == 0]
+    TMž = ž[(mž+nž+tž) % 2 == 1]
+
+    č = np.arange(npart*2*nelem) # spherical wave indices for multiple particles (e.g. in a unit cell)
+    žč = č % (2* nelem)
+    tč = tž[žč]
+    mč = mž[žč]
+    nč = nž[žč]
+    TEč = č[(mč+nč+tč) % 2 == 0]
+    TMč = č[(mč+nč+tč) % 2 == 1]
+    return (TEč, TMč)
+