diff --git a/config.txt b/config.txt
new file mode 100644
index 0000000..5c60c0d
--- /dev/null
+++ b/config.txt
@@ -0,0 +1,20 @@
+# Simulation part
+num_walker=20000
+# Motion model part
+delta=161e3
+eta=0.0
+# Distribution part
+tau=1e-1
+# Spectrum part
+dwell_time=1e-6
+num_acq=4096
+techo_start=0e-6
+techo_stop=40e-6
+techo_steps=5
+# STE part
+tevo_start=0e-6
+tevo_stop=60e-6
+tevo_steps=121
+tmix_start=1e-6
+tmix_stop=1e0
+tmix_steps=51
diff --git a/test.py b/test.py
new file mode 100644
index 0000000..37a062d
--- /dev/null
+++ b/test.py
@@ -0,0 +1,93 @@
+import pathlib
+import subprocess
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+
+def run_sims(taus):
+    for tau in taus:
+        with pathlib.Path('./config.txt').open('a') as f:
+            f.write(f'tau={tau}\n')
+        subprocess.run(['./cmake-build-debug/rwsim', '-ste', './config.txt'])
+
+
+def dampening(x: np.ndarray, apod: float) -> np.ndarray:
+    return np.exp(-apod * x)
+
+
+def pulse_attn(freq: np.ndarray, t_pulse: float):
+    # cf. Schmitt-Rohr/Spieß eq. 2.126; omega_1 * t_p = pi/2
+    pi_half_squared = np.pi**2 / 4
+    omega = 2 * np.pi * freq
+
+    numerator = np.sin(np.sqrt(pi_half_squared + omega**2 * t_pulse**2 / 2))
+    denominator = np.sqrt(pi_half_squared + omega**2 * t_pulse**2 / 4)
+
+    return np.pi * numerator / denominator / 2
+
+
+def post_process_spectrum(taus, apod, tpulse):
+    reduction_factor = np.zeros((taus.size, 5)) # hard-coded t_echo :(
+
+    for i, tau in enumerate(taus):
+        try:
+            raw_data = np.loadtxt(f'fid_tau={tau:.6e}.dat')
+        except OSError:
+            continue
+
+        t = raw_data[:, 0]
+        timesignal = raw_data[:, 1:]
+
+        timesignal *= dampening(t, apod)[:, None]
+        timesignal[0, :] /= 2
+
+        # FT to spectrum
+        freq = np.fft.fftshift(np.fft.fftfreq(t.size, d=1e-6))
+        spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real
+        spec *= pulse_attn(freq, t_pulse=tpulse)[:, None]
+
+        reduction_factor[i, :] = 2*timesignal[0, :]
+
+        plt.plot(freq, spec)
+        plt.show()
+
+    plt.semilogx(taus, reduction_factor, '.')
+    plt.show()
+
+
+def post_process_ste(taus):
+    for i, tau in enumerate(tau_values):
+
+        try:
+            raw_data_cc = np.loadtxt(f'coscos_tau={tau:.6e}.dat')
+            raw_data_ss = np.loadtxt(f'sinsin_tau={tau:.6e}.dat')
+        except OSError:
+            continue
+
+        t_mix = raw_data_cc[:, 0]
+
+        plt.semilogx(t_mix, raw_data_cc[:, 1:])
+        plt.show()
+
+        plt.semilogx(t_mix, raw_data_ss[:, 1:])
+        plt.show()
+
+        plt.plot(raw_data_cc[0, 1:])
+        plt.plot(raw_data_ss[0, 1:])
+        plt.show()
+
+        plt.plot(raw_data_cc[-1, 1:]/raw_data_cc[0, 1:])
+        plt.plot(raw_data_ss[-1, 1:]/raw_data_ss[0, 1:])
+        plt.show()
+
+
+if __name__ == '__main__':
+    tau_values = np.logspace(-3, -2, 1)
+    lb = 2e3
+    pulse_length = 2e-6
+
+    run_sims(tau_values)
+    post_process_ste(tau_values)
+    # post_process_spectrum(tau_values, lb, pulse_length)
\ No newline at end of file