diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..89fd415
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "simulations/airshower_beacon_simulation/earsim"]
+	path = simulations/airshower_beacon_simulation/earsim
+	url = https://gitlab.com/harmscho/earsim.git
diff --git a/simulations/airshower_beacon_simulation/README.md b/simulations/airshower_beacon_simulation/README.md
new file mode 100644
index 0000000..c670c8c
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/README.md
@@ -0,0 +1,7 @@
+# Airshower + Beacon simulation
+
+Simulate receiving an airshower and beacon in an array.
+Using the beacon, the clocks can be calibrated.
+
+The ZHaires simulated airshower is stored in [./ZH_airshower](./ZH_airshower).
+The produced files can be read using [./earsim](./earsim).
diff --git a/simulations/airshower_beacon_simulation/ZH_airshower/.gitignore b/simulations/airshower_beacon_simulation/ZH_airshower/.gitignore
new file mode 100644
index 0000000..d6b7ef3
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/ZH_airshower/.gitignore
@@ -0,0 +1,2 @@
+*
+!.gitignore
diff --git a/simulations/airshower_beacon_simulation/earsim b/simulations/airshower_beacon_simulation/earsim
new file mode 160000
index 0000000..47febcf
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/earsim
@@ -0,0 +1 @@
+Subproject commit 47febcf2803d2baadcc25ed61977ebb3aa49bc14
diff --git a/simulations/airshower_beacon_simulation/view_ant0.py b/simulations/airshower_beacon_simulation/view_ant0.py
new file mode 100755
index 0000000..2904c3a
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/view_ant0.py
@@ -0,0 +1,80 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from earsim import REvent
+
+
+def plot_antenna_Efields(antenna, ax=None, plot_Ex=True, plot_Ey=True, plot_Ez=True, **kwargs):
+    """Show waveforms from an antenna"""
+    if ax is None:
+        ax = plt.gca()
+
+    default_kwargs = dict(alpha=0.6)
+
+    kwargs = {**default_kwargs, **kwargs}
+
+    i = antenna.name
+
+    if plot_Ex:
+        ax.plot(antenna.t, antenna.Ex, label="E{x}".format(x='x', i=i), **kwargs)
+    if plot_Ey:
+        ax.plot(antenna.t, antenna.Ey, label="E{x}".format(x='y', i=i), **kwargs)
+    if plot_Ez:
+        ax.plot(antenna.t, antenna.Ez, label="E{x}".format(x='z', i=i), **kwargs)
+
+    ax.set_xlabel("[ns]")
+    ax.set_ylabel("[$\mu$V/m]")
+    ax.set_title("Antenna {i}".format(i=i))
+    ax.legend()
+
+    return ax
+
+def plot_antenna_max_values(antennas, ax=None, log_max=False, plot_names=True, **kwargs):
+    """Show the max values at each antenna."""
+    default_kwargs = dict(
+            cmap='Spectral_r'
+            )
+    
+    kwargs = {**default_kwargs, **kwargs}
+
+    x = [ a.x for a in antennas ]
+    y = [ a.y for a in antennas ]
+    
+    max_vals = np.asarray([ np.max([np.max(a.Ex), np.max(a.Ey), np.max(a.Ez)]) for a in antennas ])
+
+    if log_max:
+        max_vals = np.log10(max_vals)
+
+
+    if ax is None:
+        ax = plt.gca()
+
+    sc = ax.scatter(x, y, c=max_vals, **kwargs)
+    fig = ax.get_figure()
+    fig.colorbar(sc, ax=ax, label="[$\mu$V/m]")
+
+    if plot_names:
+        [ ax.annotate(a.name, (a.x, a.y), ha='center',va='center') for a in antennas ]
+    
+    ax.set_title("Maximum E field at each antenna")
+    ax.set_ylabel("[m]")
+    ax.set_xlabel("[m]")
+
+
+    return ax
+
+if __name__ == "__main__":
+    fname = "ZH_airshower/mysim.sry"
+    i = 0
+
+    ev = REvent(fname)
+
+    fig, ax1 = plt.subplots()
+    plot_antenna_Efields(ev.antennas[i], ax=ax1)
+
+    fig2, ax2 = plt.subplots()
+    plot_antenna_max_values(ev.antennas, ax=ax2)
+
+    plt.show()