diff --git a/simulations/airshower_beacon_simulation/ba_beacon_phases.py b/simulations/airshower_beacon_simulation/ba_beacon_phases.py
index b4fffac..74d66ba 100755
--- a/simulations/airshower_beacon_simulation/ba_beacon_phases.py
+++ b/simulations/airshower_beacon_simulation/ba_beacon_phases.py
@@ -16,7 +16,8 @@ import lib
 if __name__ == "__main__":
     from os import path
     import sys
-
+    import matplotlib
+    matplotlib.use('Agg')
 
     f_beacon_band = (49e-3,55e-3) #GHz
     allow_frequency_fitting = False
@@ -35,6 +36,8 @@ if __name__ == "__main__":
     fname_dir = path.dirname(fname)
     antennas_fname = path.join(fname_dir, beacon.antennas_fname)
 
+    fig_dir = "./figures" # set None to disable saving
+
     if not path.isfile(antennas_fname):
         print("Antenna file cannot be found, did you try generating a beacon?")
         sys.exit(1)
@@ -140,7 +143,7 @@ if __name__ == "__main__":
             # for reporting using plots
             found_data[i] = frequency, phase, amplitude
 
-            if show_plots and (i == 60 or i == 72):
+            if (show_plots or fig_dir) and (i == 60 or i == 72):
                 fig, ax = plt.subplots()
                 ax.set_title(f"Beacon at antenna {h5ant.attrs['name']}\nF:{frequency:.2e}, P:{phase:.4f}, A:{amplitude:.1e}")
                 ax.set_xlabel("t [ns]")
@@ -154,6 +157,9 @@ if __name__ == "__main__":
 
                 ax.legend()
 
+                if fig_dir:
+                    fig.savefig(path.join(fig_dir, __file__ + f".A{h5ant.attrs['name']}.pdf"))
+
             # save to file
             h5beacon_info = h5ant.require_group('beacon_info')
 
@@ -178,18 +184,23 @@ if __name__ == "__main__":
     print("Beacon Phases, Amplitudes and Frequencies written to", antennas_fname)
 
     # show histogram of found frequencies
-    if show_plots:
+    if show_plots or fig_dir:
         if True or allow_frequency_fitting:
             fig, ax = plt.subplots()
             ax.set_xlabel("Frequency")
             ax.set_ylabel("Counts")
             ax.axvline(f_beacon, ls='dashed', color='g')
             ax.hist(found_data[:,0], bins='sqrt', density=False)
+            if fig_dir:
+                fig.savefig(path.join(fig_dir, __file__ + f".hist_freq.pdf"))
 
         if True:
             fig, ax = plt.subplots()
             ax.set_xlabel("Amplitudes")
             ax.set_ylabel("Counts")
             ax.hist(found_data[:,2], bins='sqrt', density=False)
+            if fig_dir:
+                fig.savefig(path.join(fig_dir, __file__ + f".hist_amp.pdf"))
 
-    plt.show()
+    if show_plots:
+        plt.show()
diff --git a/simulations/airshower_beacon_simulation/bb_beacon_true_phases.py b/simulations/airshower_beacon_simulation/bb_beacon_true_phases.py
index 1eadd84..704ccb9 100755
--- a/simulations/airshower_beacon_simulation/bb_beacon_true_phases.py
+++ b/simulations/airshower_beacon_simulation/bb_beacon_true_phases.py
@@ -37,6 +37,9 @@ if __name__ == "__main__":
     from os import path
     import sys
 
+    import matplotlib
+    matplotlib.use('Agg')
+
     fname = "ZH_airshower/mysim.sry"
     c_light = 3e8*1e-9 # m/ns
     show_plots = True
@@ -48,6 +51,8 @@ if __name__ == "__main__":
     fname_dir = path.dirname(fname)
     antennas_fname = path.join(fname_dir, beacon.antennas_fname)
 
+    fig_dir = "./figures" # set None to disable saving
+
     if not path.isfile(antennas_fname):
         print("Antenna file cannot be found, did you try generating a beacon?")
         sys.exit(1)
@@ -89,7 +94,7 @@ if __name__ == "__main__":
                 h5beacon_freq.attrs['true_phase'] = true_phases[i]
 
         # Plot True Phases at their locations
-        if show_plots:
+        if show_plots or fig_dir:
             fig, ax = plt.subplots()
             spatial_unit=None
             fig.suptitle('True phases\nf_beacon= {:2.0f}MHz'.format(f_beacon*1e3))
@@ -106,4 +111,10 @@ if __name__ == "__main__":
             sc = ax.scatter(*antenna_locs, c=true_phases, **scatter_kwargs)
             fig.colorbar(sc, ax=ax, label=color_label)
 
+            if fig_dir:
+                fig.savefig(path.join(fig_dir, __file__ + f".F{freq_name}.pdf"))
+
+    print(f"True phases written to", antennas_fname)
+
+    if show_plots:
         plt.show()
diff --git a/simulations/airshower_beacon_simulation/bc_phase_sigmas.py b/simulations/airshower_beacon_simulation/bc_phase_sigmas.py
index b69bc9c..48b803e 100755
--- a/simulations/airshower_beacon_simulation/bc_phase_sigmas.py
+++ b/simulations/airshower_beacon_simulation/bc_phase_sigmas.py
@@ -14,6 +14,9 @@ if __name__ == "__main__":
     from os import path
     import sys
 
+    import matplotlib
+    matplotlib.use('Agg')
+
     fname = "ZH_airshower/mysim.sry"
     c_light = 3e8*1e-9
     show_plots = not True
@@ -24,6 +27,8 @@ if __name__ == "__main__":
     antennas_fname = path.join(fname_dir, beacon.antennas_fname)
     time_diffs_fname = 'time_diffs.hdf5' if not True else antennas_fname
 
+    fig_dir = "./figures" # set None to disable saving
+
     if not path.isfile(antennas_fname):
         print("Antenna file cannot be found, did you try generating a beacon?")
         sys.exit(1)
@@ -51,7 +56,7 @@ if __name__ == "__main__":
     # Get phase difference per baseline
     phase_diffs = np.empty( (len(baselines), 2) )
     for i, base in enumerate(baselines):
-        if i%50==0:
+        if i%100==0:
             print(i, "out of", len(baselines))
 
         # read f_beacon from the first antenna
@@ -117,4 +122,8 @@ if __name__ == "__main__":
         else:
             axs[i].plot(phase_residuals, np.arange(N_base), ls='none', marker='x')
 
+        if fig_dir:
+            fig.savefig(path.join(fig_dir, __file__ + f".F{freq_name}.pdf"))
+
+    if show_plots:
         plt.show()
diff --git a/simulations/airshower_beacon_simulation/ca_period_from_shower.py b/simulations/airshower_beacon_simulation/ca_period_from_shower.py
index 9a3997e..f167eef 100755
--- a/simulations/airshower_beacon_simulation/ca_period_from_shower.py
+++ b/simulations/airshower_beacon_simulation/ca_period_from_shower.py
@@ -10,6 +10,7 @@ import h5py
 from itertools import combinations, zip_longest, product
 import matplotlib.pyplot as plt
 import numpy as np
+from os import path
 from scipy.interpolate import interp1d
 
 from earsim import REvent
@@ -18,7 +19,7 @@ import aa_generate_beacon as beacon
 import lib
 from lib import rit
 
-def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_sample_shifts, dt=None, sample_shift_first_trace=0, plot_iteration_with_shifted_trace=None):
+def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_sample_shifts, dt=None, sample_shift_first_trace=0, plot_iteration_with_shifted_trace=None, fig_dir=None, fig_distinguish=None):
     """
     Find the best sample_shift for each antenna by summing the antenna traces
     and seeing how to get the best alignment.
@@ -63,6 +64,7 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp
             best_sample_shifts[i] = sample_shift_first_trace
             continue
 
+        # init figure
         if i == plot_iteration_with_shifted_trace:
             fig, ax = plt.subplots()
             ax.set_title("Traces at ({:.1f},{:.1f},{:.1f})".format(*test_loc))
@@ -84,18 +86,26 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp
         best_sample_shifts[i] = allowed_sample_shifts[best_idx]
         a_sum += np.roll(a_int, best_sample_shifts[i])
 
+        # cleanup figure
+        if i == plot_iteration_with_shifted_trace:
+            if fig_dir: # note this is a global variable here
+                fig.savefig(path.join(fig_dir, __file__ + '.loc{:.1f}-{:.1f}-{:.1f}'.format(*test_loc) + f'.i{i}{fig_distinguish}.pdf'))
+            plt.close(fig)
+
     # Return ks
     return best_sample_shifts, np.max(a_sum)
 
 if __name__ == "__main__":
-    from os import path
     import sys
+    import matplotlib
+    matplotlib.use('Agg')
 
     atm = AtmoCal()
 
     fname = "ZH_airshower/mysim.sry"
 
-    savepath = "./periods_from_shower_figures/"
+    fig_dir = "./figures/periods_from_shower_figures/"
+    fig_subdir = path.join(fig_dir, 'shifts/')
 
     allowed_ks = np.arange(-2, 3, 1)
     Xref = 400
@@ -173,15 +183,19 @@ if __name__ == "__main__":
         xx = []
         yy = []
         N_loc = len(maxima_per_loc)
+
         for i, (x_, y_) in enumerate(product(x,y)):
+            tmp_fig_subdir = None
             if i % 10 ==0:
                 print(f"Testing location {i} out of {N_loc}")
+                tmp_fig_subdir = fig_subdir
+
             test_loc = (x_+xoff)* ev.uAxB + (y_+yoff)*ev.uAxAxB + dXref *ev.uA
             xx.append(x_+xoff)
             yy.append(y_+yoff)
 
             # Find best k for each antenna
-            shifts, maximum = find_best_sample_shifts_summing_at_location(test_loc, ev.antennas, allowed_sample_shifts, dt=dt)
+            shifts, maximum = find_best_sample_shifts_summing_at_location(test_loc, ev.antennas, allowed_sample_shifts, dt=dt, fig_dir=tmp_fig_subdir, plot_iteration_with_shifted_trace=len(ev.antennas)-1, fig_distinguish=f".run{r}")
 
             # Translate sample shifts back into period multiple k
             ks = np.rint(shifts*f_beacon*dt)
@@ -194,7 +208,7 @@ if __name__ == "__main__":
         locs = list(zip(xx, yy))
 
         ## Save maxima to file
-        np.savetxt(savepath + __file__+f'.maxima.run{r}.txt', np.column_stack((locs, maxima_per_loc, ks_per_loc)) )
+        np.savetxt(fig_dir + __file__+f'.maxima.run{r}.txt', np.column_stack((locs, maxima_per_loc, ks_per_loc)) )
 
         if True: #plot maximum at test locations
             fig, axs = plt.subplots()
@@ -204,11 +218,12 @@ if __name__ == "__main__":
             sc = axs.scatter(xx/1e3, yy/1e3, c=maxima_per_loc, cmap='Spectral_r', alpha=0.6)
             fig.colorbar(sc, ax=axs)
 
-            fig.savefig(savepath + __file__+f'.maxima.run{r}.pdf')
+            if fig_dir:
+                fig.savefig(path.join(fig_dir, __file__+f'.maxima.run{r}.pdf'))
 
         ## Save ks to file
         best_idx = np.argmax(maxima_per_loc)
-        np.savetxt(savepath + __file__+f'.bestk.run{r}.txt', ks_per_loc[best_idx] )
+        np.savetxt(fig_dir + __file__+f'.bestk.run{r}.txt', ks_per_loc[best_idx] )
         print('Best k:', ks_per_loc[best_idx])
 
         # Abort if no improvement
diff --git a/simulations/airshower_beacon_simulation/figures/.gitignore b/simulations/airshower_beacon_simulation/figures/.gitignore
new file mode 100644
index 0000000..d6b7ef3
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/figures/.gitignore
@@ -0,0 +1,2 @@
+*
+!.gitignore
diff --git a/simulations/airshower_beacon_simulation/figures/periods_from_shower_figures/shifts/.gitignore b/simulations/airshower_beacon_simulation/figures/periods_from_shower_figures/shifts/.gitignore
new file mode 100644
index 0000000..d6b7ef3
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/figures/periods_from_shower_figures/shifts/.gitignore
@@ -0,0 +1,2 @@
+*
+!.gitignore