diff --git a/simulations/airshower_beacon_simulation/ca_period_from_shower.py b/simulations/airshower_beacon_simulation/ca_period_from_shower.py
index 9d9305a..bb54ae6 100755
--- a/simulations/airshower_beacon_simulation/ca_period_from_shower.py
+++ b/simulations/airshower_beacon_simulation/ca_period_from_shower.py
@@ -339,19 +339,23 @@ if __name__ == "__main__":
 
         if True: #plot maximum at test locations
             fig, axs = plt.subplots()
-            axs.set_title(f"Optimizing signal strength, Grid Run {r}")
+            axs.set_title(f"Optimizing signal strength varying k per antenna,\n Grid Run {r}")
             axs.set_ylabel("vxvxB [km]")
             axs.set_xlabel("  vxB [km]")
             axs.set_aspect('equal', 'datalim')
             sc = axs.scatter(xx/1e3, yy/1e3, c=maxima_per_loc, cmap='Spectral_r', alpha=0.6)
-            fig.colorbar(sc, ax=axs)
+            fig.colorbar(sc, ax=axs, label='Max Amplitude [$\\mu V/m$]')
+
+            # indicate maximum value
+            idx = np.argmax(maxima_per_loc)
+            axs.plot(xx[idx]/1e3, yy[idx]/1e3, 'bx', ms=30)
 
             if fig_dir:
                 old_xlims = axs.get_xlim()
                 old_ylims = axs.get_ylim()
                 fig.tight_layout()
                 fig.savefig(path.join(fig_dir, __file__+f'.maxima.run{r}.pdf'))
-                if True:
+                if False:
                     axs.plot(tx.x/1e3, tx.y/1e3, marker='X', color='k')
                     fig.tight_layout()
                     fig.savefig(path.join(fig_dir, __file__+f'.maxima.run{r}.with_tx.pdf'))
@@ -372,19 +376,51 @@ if __name__ == "__main__":
         ## Do a small reconstruction of the shower for best ks
         if True:
             print("Reconstructing for best k")
-            _, __, p, ___ = rit.shower_plane_slice(ev, X=Xref, Nx=len(x), Ny=len(y), wx=scale2d, wy=scale2d, xoff=xoff, yoff=yoff, zgr=0)
 
-            fig, axs = plt.subplots()
-            axs.set_title(f"Shower slice for best k, Grid Run {r}")
-            axs.set_ylabel("vxvxB [km]")
-            axs.set_xlabel("  vxB [km]")
-            axs.set_aspect('equal', 'datalim')
-            sc = axs.scatter(xx/1e3, yy/1e3, c=p, cmap='Spectral_r', alpha=0.6)
-            fig.colorbar(sc, ax=axs)
+            for j in range(2):
+                power_reconstruction = j==1
 
-            if fig_dir:
-                fig.tight_layout()
-                fig.savefig(path.join(fig_dir, __file__+f'.reconstruction.run{r}.pdf'))
+                if power_reconstruction: # Do power reconstruction
+                    # backup antenna times
+                    backup_times = [ ant.t_AxB for ant in ev.antennas ]
+
+                    # incorporate ks into timing
+                    for i, ant in enumerate(ev.antennas):
+                        ev.antennas[i].t_AxB = ant.t_AxB - best_k[i] * 1/f_beacon
+
+                    xx, yy, p, ___ = rit.shower_plane_slice(ev, X=Xref, Nx=len(x), Ny=len(y), wx=x[-1]-x[0], wy=y[-1]-y[0], xoff=xoff, yoff=yoff, zgr=0)
+
+                    # repair antenna times
+                    for i, backup_t_AxB in enumerate(backup_times):
+                        ev.antennas[i].t_AxB = backup_t_AxB
+                else: # get maximum amplitude at each location
+                    maxima = np.empty( len(locs) )
+                    for i, loc in enumerate(locs):
+                        test_loc = loc[0]* ev.uAxB + loc[1]*ev.uAxAxB + dXref *ev.uA
+                        P, t_, a_, a_sum, t_sum = rit.pow_and_time(test_loc, ev, dt=dt)
+                        maxima[i] = np.max(a_sum)
+
+                fig, axs = plt.subplots()
+                axs.set_title(f"Shower slice for best k, Grid Run {r}")
+                axs.set_ylabel("vxvxB [km]")
+                axs.set_xlabel("  vxB [km]")
+                axs.set_aspect('equal', 'datalim')
+                if power_reconstruction:
+                    sc = axs.scatter(xx/1e3, yy/1e3, c=p, cmap='Spectral_r', alpha=0.6)
+                    fig.colorbar(sc, ax=axs, label='Power')
+                else:
+                    sc = axs.scatter(xx/1e3, yy/1e3, c=maxima, cmap='Spectral_r', alpha=0.6)
+                    fig.colorbar(sc, ax=axs, label='Max Amplitude [$\\mu V/m$]')
+
+
+                if fig_dir:
+                    if power_reconstruction:
+                        fname_extra = "power"
+                    else:
+                        fname_extra = "max_amp"
+
+                    fig.tight_layout()
+                    fig.savefig(path.join(fig_dir, __file__+f'.reconstruction.run{r}.{fname_extra}.pdf'))
 
         # Abort if no improvement
         if ( r!= 0 and (old_ks_per_loc == ks_per_loc[best_idx]).all() ):