diff --git a/m-metastatis/matrixcalc.py b/m-metastatis/matrixcalc.py
index 28d8d9cc149229cb2c07925837b2b80607d2c101..058edab871a073de49b0571b01fd898be3037410 100644
--- a/m-metastatis/matrixcalc.py
+++ b/m-metastatis/matrixcalc.py
@@ -1,13 +1,17 @@
 N=4
 M=2**N
 A= [[0] * M for i in range(M) ]  # waw! creatio of a two-dminesional array 
+DIST=[-1] * M
+
+preimage= [set() for i in range(M)]
+
 
 # local transition rule in thirds
 PTRANS=[0, 2, 2, 0]
 
 
 # from decimal to binarray
-def arrayconfig(iconfig):
+def config2array(iconfig):
     array=[]
     for i in range(N):
         bit=iconfig%2
@@ -15,6 +19,9 @@ def arrayconfig(iconfig):
         iconfig//=2
     return array
 
+arrayconfig=[config2array(i) for i in range(M)]
+
+
 # one active transiton in place pos
 def transAc(x, pos):
     y= x.copy()
@@ -26,26 +33,66 @@ def prettyprint(matrix):
     for i in range(len(matrix)):
         print(matrix[i])
 
-################################
-for iconfig in range(M):
-    x= arrayconfig(iconfig)
-    for pos in range(N):
-        T1= transAc(x, pos) if x[pos]==0 else x
-        T0= transAc(x, pos) if x[pos]==1 else x
-        index1= iconfig + 2 **pos if x[pos]==0 else iconfig
-        index0= iconfig - 2 **pos if x[pos]==1 else iconfig
-        s=x[pos-1]+x[pos]+x[(pos+1)%N]
-        part1=" ic:%d(%s) pos:%d S:%d"%(iconfig, x, pos, s)
-        parta= "i1:%d(%s) %d"%(index1, str(T1), PTRANS[s])
-        partb= "i0:%d(%s) %d"%(index0, str(T0),  3- PTRANS[s])
-        print(part1, parta, partb)
 
-        A[iconfig][index1]+= PTRANS[s]
-        A[iconfig][index0]+= 3 - PTRANS[s]
 
-    print("--")
+#after the transition matrix has been calculated
+def calculatePreimages():
+    for iconfig in range(M):
+        for jconfig in range(M):
+            if (iconfig!=jconfig) and A[iconfig][jconfig]>0:
+                preimage[jconfig].add(iconfig)
+
+
+def printPreimages():
+    for i in range(M):
+        out=""
+        preI=preimage[i]
+        for indexC in preI:
+            cfg=arrayconfig[indexC]
+            out+="%d:%s  "%(indexC,cfg)
+        print("%d <= %s"%(i,out))
+
+
+def calculateDist():
+    C=set()
+    C.add(0)
+    d=0
+    while C :
+        N=set()
+        for x in C:
+            DIST[x]= d
+            print("setting %d::%s at dist:%d"%(x, arrayconfig[x], d))
+            for y in preimage[x]:
+                if DIST[y]<0:
+                    N.add(y)
+                    
+        d+=1
+        C=N
+
+################################
+def calculateTransitionMatrix():
+    for iconfig in range(M):
+        x= arrayconfig[iconfig]
+        for pos in range(N):
+            T1= transAc(x, pos) if x[pos]==0 else x
+            T0= transAc(x, pos) if x[pos]==1 else x
+            index1= iconfig + 2 **pos if x[pos]==0 else iconfig
+            index0= iconfig - 2 **pos if x[pos]==1 else iconfig
+            s=x[pos-1]+x[pos]+x[(pos+1)%N]
+            part1=" ic:%d(%s) pos:%d S:%d"%(iconfig, x, pos, s)
+            parta= "i1:%d(%s) %d"%(index1, str(T1), PTRANS[s])
+            partb= "i0:%d(%s) %d"%(index0, str(T0),  3- PTRANS[s])
+            print(part1, parta, partb)
 
+            A[iconfig][index1]+= PTRANS[s]
+            A[iconfig][index0]+= 3 - PTRANS[s]
 
 
+print(arrayconfig)
 
-prettyprint(A)
\ No newline at end of file
+calculateTransitionMatrix()
+prettyprint(A)
+calculatePreimages()
+printPreimages()
+print("%%%%%%%%%%%%%%%")
+calculateDist()
\ No newline at end of file