diff --git a/images/re_vs_pdc_eve.png b/images/re_vs_pdc_eve.png
index e466e64..8ac4b10 100644
Binary files a/images/re_vs_pdc_eve.png and b/images/re_vs_pdc_eve.png differ
diff --git a/model/eve.py b/model/eve.py
index 21bc49b..88c81cb 100644
--- a/model/eve.py
+++ b/model/eve.py
@@ -28,3 +28,32 @@ class EveBS:
         else:
             self.obtained_data.append(None)
         return passed
+
+
+@dataclass
+class EvePNS:
+    k: int = 1
+    my_basises: typing.List[bool] = field(default_factory=list)
+    obtained_data: typing.List[typing.Optional[bool]] = field(default_factory=list)
+
+    def process_photons(self, photons: typing.List[typing.Tuple[bool, bool]]) -> typing.List[typing.Tuple[bool, bool]]:
+        basis = choice([False, True])
+        self.my_basises.append(basis)
+        clicks = [0, 0]
+        passed = []
+        captured = 0
+        for photon in photons:
+            if captured < self.k:
+                captured += 1
+                photon_bit, photon_basis = photon
+                if photon_basis == basis:
+                    clicks[photon_bit] += 1
+                else:
+                    clicks[choice([0, 1])] += 1
+            else:
+                passed.append(photon)
+        if any(clicks) and not all(clicks):
+            self.obtained_data.append(bool(clicks[1]))
+        else:
+            self.obtained_data.append(None)
+        return passed
diff --git a/model/main.py b/model/main.py
index b2861d5..1999e0f 100644
--- a/model/main.py
+++ b/model/main.py
@@ -1,8 +1,12 @@
+import json
+import os.path
+
 from alice import Alice
 from bob import Bob
 from channel import ChannelSym
 from threading import Thread
 import typing
+import sympy as sp
 
 from eve import EveBS
 
@@ -21,11 +25,17 @@ def get_e_and_r(alice, bob, n=1000) -> typing.Tuple[float, float]:
 
 
 channel_parameters = {'p_opt': 0.05, 'p_dc': 0.05, 'mu': 1, 'detector_sensitivity': 0.8, 'transmittance': 0.8}
+PAR_NAMES = {'p_opt': 'p_{opt}', 'p_dc': 'p_{dc}', 'mu': r'\mu', 'detector_sensitivity': r'\eta', 'transmittance': 't'}
+mu, t, p_dc, p_opt, eta = sp.symbols(r'\mu t p_{dc} p_{opt} \eta')
+e_theory = 0.5 * (p_dc + p_opt * mu * t * eta) / (2 * p_dc + mu * t * eta)
+q_theory = (mu * t * eta + 2 * p_dc) / 2 * (1 - e_theory * 2)
 
 
-def plot(parameter, values, add_eve: bool = True):
+def plot(parameter, values, add_eve: bool = False):
     data_r, data_e = [], []
-    n = 1000
+    data_r_theoretical, data_e_theoretical = [], []
+    parameters_sp = {PAR_NAMES[pname]: channel_parameters[pname] for pname in PAR_NAMES.keys()}
+    n = 10000
     for val in values:
         print(val)
         channel = ChannelSym(**{parameter: val})
@@ -37,12 +47,22 @@ def plot(parameter, values, add_eve: bool = True):
         e, r = get_e_and_r(alice, bob, n)
         data_r.append(r)
         data_e.append(e)
+        parameters_sp[PAR_NAMES.get(parameter, parameter)] = val
+        data_r_theoretical.append(float(q_theory.evalf(subs=parameters_sp)))
+        data_e_theoretical.append(float(e_theory.evalf(subs=parameters_sp)))
     from matplotlib import pyplot as plt
-    plt.plot(values, data_r)
-    plt.plot(values, data_e)
-    plt.legend(['R', 'E'])
-    plt.xlabel({'p_dc': '$p_{dc}$'}.get(parameter, parameter))
-    plt.title('$R$ and $E$ vs. ' + {'p_dc': '$p_{dc}$'}.get(parameter, parameter) + ' with Eve' * add_eve)
+    plt.plot(values, data_r, label='$R$')
+    plt.plot(values, data_e, label='$E$')
+    plt.plot(values, data_r_theoretical, label='$R_{theory}$')
+    plt.plot(values, data_e_theoretical, label='$E_{theory}$')
+    plt.legend()
+    if not os.path.exists('output'):
+        os.mkdir('output')
+    with open(f'output/dep_{parameter}.json', 'w') as f:
+        f.write(json.dumps([list(values), data_r, data_e, data_r_theoretical, data_e_theoretical]))
+    plt.xlabel('$' + PAR_NAMES.get(parameter, parameter) + '$')
+    plt.title('$R$ and $E$ vs. $' + PAR_NAMES.get(parameter, parameter) + '$' + ' with Eve' * add_eve)
+    plt.savefig(f'output/dep_{parameter}.png')
     plt.show()
 
 
@@ -75,3 +95,5 @@ def run_one():
 if __name__ == '__main__':
     import numpy as np
     plot('p_dc', np.arange(0, 0.1, 0.01))
+    plot('detector_sensitivity', np.arange(0.5, 1, 0.05))
+    plot('mu', np.arange(0.5, 1.5, 0.1))