diff --git a/qpms/__init__.py b/qpms/__init__.py
index 3e52d07..7a423af 100644
--- a/qpms/__init__.py
+++ b/qpms/__init__.py
@@ -6,3 +6,4 @@ from .qpms_p import *
 from .lattices2d import *
 from .hexpoints import *
 from .tmatrices import *
+
diff --git a/qpms/qpms_p.py b/qpms/qpms_p.py
index 3c78cc7..22e7d68 100644
--- a/qpms/qpms_p.py
+++ b/qpms/qpms_p.py
@@ -783,11 +783,16 @@ def loadWfile(fileName, lMax = None, fatForm = True):
             'freqs_weirdunits' : freqs_weirdunits,
             }
 
-def processWfiles_sameKs(freqfilenames, destfilename, lMax = None):
-    '''Processes translation operator data in different files; each file is supposed to contain one frequency.
+def processWfiles_sameKs(freqfilenames, destfilename, lMax = None, f='npz'):
+    '''
+    Processes translation operator data in different files; each file is supposed to contain one frequency.
     The Ks in the different files are expected to be exactly the same and in the same order; 
     the result is sorted by frequencies and saved to npz file;
-    The representation is always "thin", i.e. the elements which are zero due to z-symmetry are skipped'''
+    The representation is always "thin", i.e. the elements which are zero due to z-symmetry are skipped
+    F is either 'npz' or 'd' 
+    'npz' creates npz archive, 'd' creates a directory with individual npy files, where the W-matrix data
+    are writen using memmap.
+    '''
     nparticles = 2 #TODO generalize
     
     um = 1e-6 # micrometre in SI units
@@ -808,7 +813,17 @@ def processWfiles_sameKs(freqfilenames, destfilename, lMax = None):
         nelem = lMax2nelem(lMax)
     ks_muster = np.array(data[:,3:5], copy=True)
     
-    allWs = np.empty((nfiles, nk_muster, nparticles, nparticles, 2, nelem, nelem,), dtype=complex)
+    if f == 'npz':
+        allWs = np.empty((nfiles, nk_muster, nparticles, nparticles, 2, nelem, nelem,), dtype=complex)
+    elif f == 'd':
+        j = os.path.join
+        d = destfilename
+        try:
+            os.mkdir(destfilename)
+        except FileExistsError as exc:
+            pass
+        allWs = np.lib.format.open_memmap(filename=j(d,"Wdata.npy"), dtype=complex, shape=(nfiles, nk_muster, nparticles, nparticles,2,nelem,nelem,),mode='w+')
+    else: raise ValueError("f has to be either 'npz' or 'd'")
     k0s = np.empty((nfiles,))
     freqs_weirdunits = np.empty((nfiles,))
     EeVs = np.empty((nfiles,))
@@ -832,12 +847,11 @@ def processWfiles_sameKs(freqfilenames, destfilename, lMax = None):
         k0s[succread] = data[0,2] # TODO check this as well?
         EeVs[succread] = data[0,1]
         succread += 1
+    allWs.close()
     freqs = freqs_weirdunits * c / um
 
-
-    #TODO TODO TODO sort everything by omega
-    
-    np.savez(destfilename,
+    if f == 'npz':
+        np.savez(destfilename,
              lMax = lMax,
              nelem = nelem,
              npart = nparticles,
@@ -850,16 +864,54 @@ def processWfiles_sameKs(freqfilenames, destfilename, lMax = None):
              ks = ks_muster,
              Wdata = allWs[:succread]
             )
+    elif f == 'd':
+        np.save(j(d,'lMax.npy'), lMax)
+        np.save(j(d,'nelem.npy'), nelem)
+        np.save(j(d,'npart.npy'), nparticles)
+        np.save(j(d,'nfreqs.npy'), succread)
+        np.save(j(d,'nk.npy'), nk_muster)
+        np.save(j(d,'freqs.npy'), freqs[:succread])
+        np.save(j(d,'freqs_weirdunits.npy'), freqs_weirdunits[:succread])
+        np.save(j(d,'EeVs_orig.npy'), EeVs[:succread])
+        np.save(j(d,'k0s.npy'), k0s[:succread])
+        np.save(j(d,'ks.npy'), ks_muster)
+        # Wdata already saved
     return
 
 
 
-def loadWfile_processed(fileName, lMax = None, fatForm = True):
-    data = np.load(fileName);
-    nk = data['nk'][()]
-    nfreqs = data['nfreqs'][()]
-    nparticles = data['npart'][()]
-    Ws = data['Wdata']
+def loadWfile_processed(fileName, lMax = None, fatForm = True, midk_halfwidth = None):
+    if os.path.isdir(fileName): # .npy files in a directory
+        p = filename
+        j = os.path.join
+        nk = np.load(j(p,'nk.npy'))[()]
+        nfreqs = np.load(j(p,'nfreqs.npy'))[()]
+        nparticles = np.load(j(p,'npart.npy'))[()]
+        Ws = np.load(j(p,'Wdata.npy'), mode='r')
+        ks = np.load(j(p,'ks.npy'))
+        freqs = np.load(j(p,'freqs.npy'))
+        k0s = np.load(j(p,'k0s.npy'))
+        EeVs_orig = np.load(j(p,'EeVs_orig.npy'))
+        freqs_weirdunits = np.load(j(p,'freqs_weirdunits.npy'))
+    else: # npz file
+        data = np.load(fileName, mode='r')
+        nk = data['nk'][()]
+        nfreqs = data['nfreqs'][()]
+        nparticles = data['npart'][()]
+        Ws = data['Wdata']
+        ks = data['ks']
+        freqs = data['freqs']
+        k0s = data['k0s']
+        EeVs_orig = data['EeVs_orig']
+        freqs_weirdunits = data['freqs_weirdunits']
+    if midk_halfwidth is not None:
+        k_mid_i = nk // 2
+        if midk_halfwidth < k_mid_i:
+            maxk_i = k_mid_i + midk_halfwidth
+            mink_i = k_mid_i - midk_halfwidth
+            nk = 2*midk_halfwidth+1
+            Ws = Ws[:,mink_i:(maxk_i+1)]
+            ks = ks[mink_i:(maxk_i+1)]
     if lMax is not None:
         nelem = lMax2nelem(lMax)
         Ws = Ws[...,:nelem,:nelem]
@@ -880,11 +932,11 @@ def loadWfile_processed(fileName, lMax = None, fatForm = True):
         'npart': nparticles,
         'nfreqs': nfreqs,
         'nk' : nk,
-        'freqs': data['freqs'],
-        'freqs_weirdunits': data['freqs_weirdunits'],
-        'EeVs_orig': data['EeVs_orig'],
-        'k0s': data['k0s'],
-        'ks': data['ks'],
+        'freqs': freqs,
+        'freqs_weirdunits': freqs_weirdunits,
+        'EeVs_orig': EeVs_orig,
+        'k0s': k0s,
+        'ks': ks,
         'Ws': Ws,
     }
     
diff --git a/qpms/symmetries.py b/qpms/symmetries.py
new file mode 100644
index 0000000..ca0633f
--- /dev/null
+++ b/qpms/symmetries.py
@@ -0,0 +1,200 @@
+from sympy.combinatorics import Permutation, PermutationGroup
+Permutation.print_cyclic = True
+import cmath
+from cmath import exp, pi
+from math import sqrt
+import numpy as np
+np.set_printoptions(linewidth=200)
+import qpms
+import numbers
+ň = None
+
+def grouprep_try(tdict, src, im, srcgens, imgens, immultop = None, imcmp = None):
+    tdict[src] = im
+    for i in range(len(srcgens)):
+        new_src = src * srcgens[i]
+        new_im = (im * imgens[i]) if (immultop is None) else immultop(im, imgens[i])
+        if new_src not in tdict.keys():
+            grouprep_try(tdict, new_src, new_im, srcgens, imgens, immultop, imcmp)
+        elif ((new_im != tdict[new_src]) if (imcmp is None) else (not imcmp(new_im, tdict[new_src]))): # check consistency
+            print(src, ' * ', srcgens[i], ' --> ', new_src)
+            print(im)
+            print(' * ')
+            print(imgens[i])
+            print(' --> ')
+            print(new_im)
+            print(' != ')
+            print(tdict[new_src])
+            raise ValueError("Homomorphism inconsistency detected")                
+    return
+
+# srcgroup is expected to be PermutationGroup and srcgens of the TODO
+# imcmp returns True if two elements of the image group are 'equal', otherwise False
+def generate_grouprep(srcgroup, im_identity, srcgens, imgens, immultop = None, imcmp = None):
+    sz = srcgens[0].size
+    for g in srcgens:
+        if g.size != sz:
+            raise ValueError('All the generators must have the same "size"')
+    tdict = dict()
+    grouprep_try(tdict, Permutation(sz-1), im_identity, srcgens, imgens, immultop = immultop, imcmp = imcmp)
+    if(srcgroup.order() != len(tdict.keys())): # basic check
+        raise ValueError('The supplied "generators" failed to generate the preimage group: ', 
+                         srcgroup.order(), " != ", len(tdict.keys()))
+    return tdict
+    
+epsilon = np.eye(2)
+alif = np.array(((-1/2,-sqrt(3)/2),(sqrt(3)/2,-1/2)))
+bih = np.array(((-1/2,sqrt(3)/2),(-sqrt(3)/2,-1/2)))
+lam = np.array(((1,0),(0,-1)))
+mim =  np.array(((-1/2,-sqrt(3)/2),(-sqrt(3)/2,1/2)))
+nun =  np.array(((-1/2,sqrt(3)/2),(sqrt(3)/2,1/2)))
+
+
+# Group D3h
+# Note that the size argument of permutations is necessary, otherwise e.g. c*c and  b*b would not be evaluated equal
+# N.B. the weird elements as Permutation(N) – it means identity permutation of size N+1.
+rot3_perm = Permutation(0,1,2, size=5) # C3 rotation
+xflip_perm = Permutation(0,2, size=5) # vertical mirror
+zflip_perm = Permutation(3,4, size=5) # horizontal mirror
+D3h_srcgens = [rot3_perm,xflip_perm,zflip_perm]
+D3h_permgroup = PermutationGroup(rot3_perm,xflip_perm,zflip_perm) # D3h
+
+#srcgens = [a,b,c]
+D3h_irreps = {
+    # Bradley, Cracknell p. 61
+    'E1' : generate_grouprep(D3h_permgroup, epsilon, D3h_srcgens, [alif, lam, epsilon], immultop = np.dot, imcmp = np.allclose),
+    'E2' : generate_grouprep(D3h_permgroup, epsilon, D3h_srcgens, [alif, lam, -epsilon], immultop = np.dot, imcmp = np.allclose),
+    # Bradley, Cracknell p. 59,
+    'A1p' : generate_grouprep(D3h_permgroup, 1, D3h_srcgens, [1,1,1]),
+    'A2p' : generate_grouprep(D3h_permgroup, 1, D3h_srcgens, [1,-1,1]),
+    'A1pp' : generate_grouprep(D3h_permgroup, 1, D3h_srcgens, [1,1,-1]),
+    'A2pp' : generate_grouprep(D3h_permgroup, 1, D3h_srcgens, [1,-1,-1]),
+}
+
+
+def mmult_tyty(a, b):
+        return(qpms.apply_ndmatrix_left(a, b, (-4,-3)))
+def mmult_ptypty(a, b):
+    return(qpms.apply_ndmatrix_left(a, b, (-6,-5,-4)))
+    
+#TODO lepší název fce
+def gen_point_D3h_svwf_rep(lMax):
+    my, ny = qpms.get_mn_y(lMax)
+    nelem = len(my)
+    C3_yy = qpms.WignerD_yy_fromvector(lMax, np.array([0,0,2*pi/3]))
+    C3_tyty = np.moveaxis(np.eye(2)[:,:,ň,ň] * C3_yy, 2,1)
+    zfl_tyty = qpms.zflip_tyty(lMax)
+    yfl_tyty = qpms.yflip_tyty(lMax)
+    xfl_tyty = qpms.xflip_tyty(lMax)
+    I_tyty = np.moveaxis(np.eye(2)[:,:,ň,ň] * np.eye(nelem), 2,1)
+    order = D3h_permgroup.order()
+    sphrep_full = generate_grouprep(D3h_permgroup, I_tyty, D3h_srcgens, [C3_tyty, xfl_tyty, zfl_tyty], 
+                           immultop = mmult_tyty, imcmp = np.allclose)
+    sphreps = dict()
+    for repkey, matrixrep in D3h_irreps.items():
+        arepmatrix = matrixrep[rot3_perm] # just one of the matrices to get the shape etc
+        if isinstance(arepmatrix, numbers.Number):
+            dim = 1 # repre dimension
+            preprocess = lambda x: np.array([[x]])
+        elif isinstance(arepmatrix, np.ndarray):
+            if(len(arepmatrix.shape)) != 2 or arepmatrix.shape[0] != arepmatrix.shape[1]:
+                raise ValueError("Arrays representing irrep matrices must be of square shape")
+            dim = arepmatrix.shape[0]
+            preprocess = lambda x: x
+        else: 
+            raise ValueError("Irrep is not a square array or number")
+        sphrep = np.zeros((dim,dim,2,nelem,2,nelem), dtype=complex)
+        for i in D3h_permgroup.elements:
+            sphrep += preprocess(matrixrep[i]).conj().transpose()[:,:,ň,ň,ň,ň] * sphrep_full[i]
+        sphrep *= dim / order
+        # clean the nonexact values here 
+        for x in [0, 0.5, -0.5, 0.5j, -0.5j]:
+            sphrep[np.isclose(sphrep,x)]=x
+        sphreps[repkey] = sphrep
+    return sphreps
+        
+def gen_hexlattice_Kpoint_svwf_rep(lMax, psi):
+    my, ny = qpms.get_mn_y(lMax)
+    nelem = len(my)
+    C3_yy = qpms.WignerD_yy_fromvector(lMax, np.array([0,0,2*pi/3]))
+    C3_tyty = np.moveaxis(np.eye(2)[:,:,ň,ň] * C3_yy, 2,1)
+    zfl_tyty = qpms.zflip_tyty(lMax)
+    yfl_tyty = qpms.yflip_tyty(lMax)
+    xfl_tyty = qpms.xflip_tyty(lMax)
+    I_tyty = np.moveaxis(np.eye(2)[:,:,ň,ň] * np.eye(nelem), 2,1)
+    hex_C3_K_ptypty = np.diag([exp(-psi*1j*2*pi/3),exp(+psi*1j*2*pi/3)])[:,ň,ň,:,ň,ň] * C3_tyty[ň,:,:,ň,:,:]
+    hex_zfl_ptypty = np.eye(2)[:,ň,ň,:,ň,ň] * zfl_tyty[ň,:,:,ň,:,:]
+    hex_xfl_ptypty = np.array([[0,1],[1,0]])[:,ň,ň,:,ň,ň] * xfl_tyty[ň,:,:,ň,:,:]
+    hex_I_ptypty = np.eye((2*2*nelem)).reshape((2,2,nelem,2,2,nelem))
+    order = D3h_permgroup.order()
+    hex_K_sphrep_full = generate_grouprep(D3h_permgroup, hex_I_ptypty, D3h_srcgens, [hex_C3_K_ptypty, hex_xfl_ptypty, hex_zfl_ptypty], 
+                           immultop = mmult_ptypty, imcmp = np.allclose)
+    hex_K_sphreps = dict()
+    for repkey, matrixrep in D3h_irreps.items():
+        arepmatrix = matrixrep[rot3_perm] # just one of the matrices to get the shape etc
+        if isinstance(arepmatrix, numbers.Number):
+            dim = 1 # repre dimension
+            preprocess = lambda x: np.array([[x]])
+        elif isinstance(arepmatrix, np.ndarray):
+            if(len(arepmatrix.shape)) != 2 or arepmatrix.shape[0] != arepmatrix.shape[1]:
+                raise ValueError("Arrays representing irrep matrices must be of square shape")
+            dim = arepmatrix.shape[0]
+            preprocess = lambda x: x
+        else: 
+            raise ValueError("Irrep is not a square array or number")
+        sphrep = np.zeros((dim,dim,2,2,nelem,2,2,nelem), dtype=complex)
+        for i in D3h_permgroup.elements:
+            sphrep += preprocess(matrixrep[i]).conj().transpose()[:,:,ň,ň,ň,ň,ň,ň] * hex_K_sphrep_full[i]
+        sphrep *= dim / order
+        # clean the nonexact values here 
+        for x in [0, 0.5, -0.5, 0.5j, -0.5j]:
+            sphrep[np.isclose(sphrep,x)]=x
+        hex_K_sphreps[repkey] = sphrep
+    return hex_K_sphreps        
+
+def normalize(v):
+    norm = np.linalg.norm(v.reshape((np.prod(v.shape),)), ord=2)
+    if norm == 0: 
+       return v*np.nan
+    return v / norm
+
+def gen_hexlattice_Kpoint_svwf_rep_projectors(lMax,psi):
+    nelem = lMax * (lMax+2)
+    projectors = dict()
+    for repi, W in gen_hexlattice_Kpoint_svwf_rep(lMax,psi).items():
+        totalvecs = 0
+        tmplist = list()
+        for p in (0,1):
+         for t in (0,1):
+          for y in range(nelem):
+           for ai in range(W.shape[0]):
+            for bi in range(W.shape[1]):
+                v = np.zeros((2,2,nelem))
+                v[p,t,y] = 1
+                #v = np.ones((2,2,nelem))
+                v1 = np.tensordot(W[ai,bi],v, axes = ([-3,-2,-1],[0,1,2])) 
+
+
+                if not np.allclose(v1,0):
+                    v1 = normalize(v1)
+                    for v2 in tmplist:
+                        dot = np.tensordot(v1.conjugate(),v2,axes = ([-3,-2,-1],[0,1,2]))
+                        if not np.allclose(dot,0):
+                            if not np.allclose(np.abs(dot),1):
+                                raise ValueError('You have to fix this piece of code.')# TODO maybe I should make sure that the absolute value is around 1
+                            break
+                    else:
+                        totalvecs += 1
+                        tmplist.append(v1)
+                        #for index, x in np.ndenumerate(v1):
+                        #    if x!=0:
+                        #        print(index, x)
+                        #print('----------')
+        theprojector = np.zeros((2,2,nelem,2,2,nelem), dtype = float) 
+        for v in tmplist:
+            theprojector += (v[:,:,:,ň,ň,ň] * v.conjugate()[ň,ň,ň,:,:,:]).real # TODO check is it possible to have imaginary elements?
+        for x in [0, 1, -1,sqrt(0.5),-sqrt(0.5),0.5,-0.5]:
+            theprojector[np.isclose(theprojector,x)]=x
+        projectors[repi] = theprojector
+    return projectors
+