diff --git a/simulation/.gitignore b/simulation/.gitignore
index fbca225..4d80fa1 100644
--- a/simulation/.gitignore
+++ b/simulation/.gitignore
@@ -1 +1,4 @@
 results/
+env/
+*.csv
+
diff --git a/simulation/plot.py b/simulation/plot.py
new file mode 100644
index 0000000..8b5ba07
--- /dev/null
+++ b/simulation/plot.py
@@ -0,0 +1,37 @@
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+import numpy as numpy
+import pandas as pd
+
+def getcolor(n):
+    return ['#aa0000', '#00aa00', '#0000aa', '#aaaa00'][int((i - 12) / 2)]
+
+df = pd.read_csv('pretty.csv', delimiter='\s', engine='python')
+
+x = df['k'].unique()
+
+for i in df['i_min'].unique():
+    rows = df.loc[df['i_min'] == i]
+    plt.errorbar(
+        x,
+        rows['mean'],
+        yerr=rows['ci'],
+        marker='.',
+        capsize=4,
+        ms=3,
+        mec='r',
+        mfc='r',
+        linestyle='-',
+        linewidth=1,
+        color=getcolor(i),
+        ecolor='red'
+    )
+
+handles = []
+for i in df['i_min'].unique():
+    handles.append(mpatches.Patch(color=getcolor(i), label='Imin = ' + str(i)))
+
+plt.legend(handles=handles)
+
+plt.show()
+
diff --git a/simulation/pretty.py b/simulation/pretty.py
new file mode 100644
index 0000000..16d306a
--- /dev/null
+++ b/simulation/pretty.py
@@ -0,0 +1,138 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import csv
+import os
+import re
+import statistics
+
+### PARAMETERS ###
+RESULTSDIR='results/'
+
+### FUNCTIONS ###
+def compute_confidence_interval(samples):   # 95%
+    stdev = statistics.stdev(samples)
+    ci = 1.96 * stdev / np.sqrt(len(samples))
+    return ci
+
+### GLOBAL VARIABLES ###
+all_route_discovery = dict()
+
+### MAIN ###
+for filename in os.listdir(RESULTSDIR):
+    df = pd.read_csv(RESULTSDIR + '/' + filename,
+        delimiter='\s|\t',
+        index_col=0,
+        nrows=4,
+        engine='python',
+        header=None).T
+
+    kappa = int(df['KAPPA'].to_string().split()[1])
+    i_min = int(df['I_MIN'].to_string().split()[1])
+    i_max = int(df['I_MAX'].to_string().split()[1])
+    key = str(kappa) + ':' + str(i_min) + ':' + str(i_max)
+
+    if (not key in all_route_discovery):
+        all_route_discovery[key] = { 'samples': [], 'k': kappa, 'i_min': i_min }
+
+    with open(RESULTSDIR + '/' + filename) as file:
+        for line in file:
+            if (re.match('ALL-ROUTE-DISCOVERY', line)):
+                t = int(line.split()[1]) / (10**6)
+                all_route_discovery[key]['samples'].append(t)
+                break
+
+# now, for every k, i_min, i_max
+# all_route_discovery contains an array with all samples of the measured value
+
+print('k i_min mean ci')
+for simulation in all_route_discovery:
+    all_route_discovery[simulation]['mean'] = np.mean(all_route_discovery[simulation]['samples'])
+    all_route_discovery[simulation]['ci']   = compute_confidence_interval(all_route_discovery[simulation]['samples'])
+
+    print(
+        #simulation,
+        all_route_discovery[simulation]['k'],
+        all_route_discovery[simulation]['i_min'],
+        all_route_discovery[simulation]['mean'],
+        all_route_discovery[simulation]['ci']
+    )
+
+exit(0)
+
+'''
+d = []
+f = []
+
+# compute fd in our interval of interest
+d = np.arange(0, (M * 2**0.5 / 2))
+for i in d:
+    f.append(fd(i))
+
+# compute integral of fd
+integral = []
+last = 0
+for i in range(0, len(f)):
+    part = f[i]
+    last = last + part
+    integral.append(last)
+
+print("normalization fd: ", np.trapz(f))
+
+# make plots
+fig, ax1 = plt.subplots()
+
+color = 'tab:red'
+ax1.set_xlabel('d [m]')
+ax1.set_ylabel('fD(d)')
+ax1.plot(d, f, 'c,', label='fD')
+
+ax2 = ax1.twinx()
+
+color = 'tab:blue'
+ax2.set_ylabel('FD(d)')
+ax2.plot(d, integral, 'b,', label="FD")
+
+plt.legend()
+plt.show()
+
+def fs(s):
+    if s >= 0 and s <= (T * M**2 ) / 4:
+        return np.pi / ( T * M**2 )
+    elif s >= (T * M**2 ) / 4 and s <= (T * M**2) / 2:
+        return (np.pi / (T * M**2)) - (4 / (T * M**2)) * np.arccos(M/2 * ((T / s)**0.5))
+    else:
+        return 0
+
+### Distribution of Service Time
+s = T * d**2
+f = []
+for i in s:
+    f.append(fs(i))
+
+print("normalization fs: ", np.trapz(f, s))
+
+# compute integral of fs
+integral = [ 0 ]
+last = 0
+for i in range(1, len(f)):
+    part = f[i] * (s[i] - s[i - 1])
+    last = last + part
+    integral.append(last)
+
+# make plots
+fig, ax1 = plt.subplots()
+color = 'tab:red'
+ax1.set_xlabel('s [s]')
+ax1.set_ylabel('fS(s)')
+ax1.plot(s, f, 'c,', label='fS')
+
+ax2 = ax1.twinx()
+
+color = 'tab:blue'
+ax2.set_ylabel('FD(d)')
+ax2.plot(s, integral, 'b,', label='FS')
+
+plt.legend()
+plt.show()
+'''
\ No newline at end of file
diff --git a/simulation/requirements.txt b/simulation/requirements.txt
new file mode 100644
index 0000000..cfb5eff
--- /dev/null
+++ b/simulation/requirements.txt
@@ -0,0 +1,10 @@
+cycler==0.10.0
+kiwisolver==1.1.0
+matplotlib==3.1.2
+numpy==1.17.4
+pandas==0.25.3
+pkg-resources==0.0.0
+pyparsing==2.4.5
+python-dateutil==2.8.1
+pytz==2019.3
+six==1.13.0