ZH: testable script to show phases at antennas

simulations/airshower_beacon_simulation/aa_generate_beacon.py

@@ -167,8 +167,9 @@ if __name__ == "__main__":
 
     # make beacon per antenna
     for i, antenna in enumerate(ev.antennas):
-        beacon = lib.beacon_from(tx, antenna, f_beacon, antenna.t, radiate_rsq=beacon_radiate_rsq)
+        t0 = lib.distance(tx, antenna)/3e8 * 1e9 # ns
 
+        beacon = lib.beacon_from(tx, antenna, f_beacon, antenna.t, t0=t0, c_light=np.inf, radiate_rsq=beacon_radiate_rsq)
 
         E = np.array([antenna.Ex, antenna.Ey, antenna.Ez, beacon])
         append_antenna_hdf5( antennas_fname, antenna, E, name='orig_traces', prepend_time=True)
@@ -177,6 +178,6 @@ if __name__ == "__main__":
         for j, _ in enumerate(beacon_amplitudes):
             E[j] += beacon_amplitudes[j]*beacon
 
-        append_antenna_hdf5( antennas_fname, antenna, E, name='traces', prepend_time=True)
+        append_antenna_hdf5( antennas_fname, antenna, E, name='traces', prepend_time=True, attrs_dict=dict(t0=t0))
 
     print("Antenna HDF5 file written as " + str(antennas_fname))

simulations/airshower_beacon_simulation/lib.py

@@ -36,3 +36,29 @@ def beacon_from(tx, rx, f, t=0, t0=0, c_light=3e8, radiate_rsq=True, amplitude=1
         amplitude *= 1/(dist**2)
 
     return sine_beacon(f, t, t0=t0, amplitude=amplitude,**kwargs)
+
+
+def phase_field_from_tx(x, y, tx, f_beacon, c_light=3e8, t0=0, wrap_phase=True, return_meshgrid=True):
+    xs, ys = np.meshgrid(x, y, sparse=True)
+
+    x_distances = (tx.x - xs)**2
+    y_distances = (tx.y - ys)**2
+
+    dist = np.sqrt( x_distances + y_distances )
+
+    phase = (dist/c_light + t0) * f_beacon*2*np.pi
+    if wrap_phase:
+        phase = phase_mod(phase)
+
+    if return_meshgrid:
+        return phase, (xs, ys)
+    else:
+        return phase, (np.repeat(xs, len(ys), axis=0), np.repeat(ys, len(xs[0]), axis=1))
+
+
+def phase_mod(phase, low=np.pi):
+    """
+    Modulo phase such that it falls within the
+    interval $[-low, 2\pi - low)$.
+    """
+    return (phase + low) % (2*np.pi) - low

simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py

@@ -24,9 +24,16 @@ if __name__ == "__main__":
 
     f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
 
+    subtitle = ""
+    if True:
         beacon_frequencies = np.array([ant.attrs['beacon_freq'] for ant in antennas])
         beacon_amplitudes = np.array([ant.attrs['beacon_amplitude'] for ant in antennas])
-    beacon_phases = np.array([ant.attrs['beacon_phase'] for ant in antennas])
+        beacon_phases = np.array([lib.phase_mod(ant.attrs['beacon_phase']) for ant in antennas])
+    else:
+        subtitle = " Phases from t0"
+        beacon_frequencies = np.array([ f_beacon for ant in antennas ])
+        beacon_amplitudes = np.array([ 1 for ant in antennas ])
+        beacon_phases = np.array([ lib.phase_mod(ant.attrs['t0']*beacon_frequencies[i]*2*np.pi) for i, ant in enumerate(antennas)])
 
     #####
     sizes = 64
@@ -40,26 +47,50 @@ if __name__ == "__main__":
 
     x = [ a.x for a in antennas ]
     y = [ a.y for a in antennas ]
+
     #####
 
     fig, axs = plt.subplots()
-    axs.set_title("Amplitude at beacon frequency at each antenna")
+    axs.set_title(f"Amplitude at beacon frequency at each antenna\nf at A0: {beacon_frequencies[0]}GHz" + subtitle)
     axs.set_aspect('equal', 'datalim')
     axs.set_xlabel('[m]')
     axs.set_ylabel('[m]')
 
+    if True:
         if True:
             # underlie a calculate phase field
-        xs = np.linspace( np.min(x), np.max(x), 50)
+            if True: # only fill for antennas
-        ys = np.linspace( np.min(y), np.max(y), 50)
+                xs = np.linspace( np.min(x), np.max(x), 4*np.ceil(len(antennas)**0.5) -1 )
-        phases, (xs, ys) = lib.phase_field_from_tx(xs, ys, tx, f_beacon, return_meshgrid=False)
+                ys = np.linspace( np.min(y), np.max(y), 4*np.ceil(len(antennas)**0.5) -1)
-        sc2 = axs.scatter(xs, ys, c=phases, alpha=0.5, zorder=-5)
+            else: # make field from halfway the transmitter
+                xs = np.linspace( (tx.x - np.min(x))/2, np.max(x), 500)
+                ys = np.linspace( (tx.y - np.min(y))/2, np.max(y), 500)
+
+            phases, (xs, ys) = lib.phase_field_from_tx(xs, ys, tx, f_beacon*1e9,return_meshgrid=False)
+            sc2 = axs.scatter(xs, ys, c=phases, alpha=0.5, zorder=-5, cmap='Spectral_r')
             fig.colorbar(sc2, ax=axs)
     
-    sc = axs.scatter(x, y, c=vals, s=sizes)
+        sc = axs.scatter(x, y, c=vals, s=sizes, cmap='Spectral_r', edgecolors='k', marker='X')
     
         axs.plot(tx.x, tx.y, marker='X', color='k')
     
+        #for i, freq in enumerate(beacon_frequencies):
+        #    axs.text(f"{freq:.2e}", (x[i], y[i]))
+    
+        fig.colorbar(sc, ax=axs, label=colorlabel)
+    else:
+        phases, (xs, ys) = lib.phase_field_from_tx(x, y, tx, f_beacon*1e9, return_meshgrid=False)
+
+        phase_diffs = vals - lib.phase_mod(phases)
+        phase_diffs = lib.phase_mod(phase_diffs)
+
+        print(phases)
+
+        sc = axs.scatter(xs, ys, c=phase_diffs, s=sizes, cmap="Spectral_r")
+        axs.plot(tx.x, tx.y, marker='X', color='k')
+
         fig.colorbar(sc, ax=axs, label=colorlabel)
 
+
+    fig.savefig(__file__ + ".pdf")
     plt.show()