From 2240d67c1c745760cfef9fcfbeeeaaa54bb4b0ce Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Fri, 18 Nov 2022 16:19:44 +0100 Subject: [PATCH] ZH: testable script to show phases at antennas --- .../aa_generate_beacon.py | 5 +- .../airshower_beacon_simulation/lib.py | 26 +++++++++ .../show_beacon_amplitude_antennas.py | 57 ++++++++++++++----- 3 files changed, 73 insertions(+), 15 deletions(-) diff --git a/simulations/airshower_beacon_simulation/aa_generate_beacon.py b/simulations/airshower_beacon_simulation/aa_generate_beacon.py index 5e7fca9..eef6861 100755 --- a/simulations/airshower_beacon_simulation/aa_generate_beacon.py +++ b/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)) diff --git a/simulations/airshower_beacon_simulation/lib.py b/simulations/airshower_beacon_simulation/lib.py index b167219..0fbaa61 100644 --- a/simulations/airshower_beacon_simulation/lib.py +++ b/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 diff --git a/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py b/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py index a682c97..70e2f62 100755 --- a/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py +++ b/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) - 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]) + 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([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: - # underlie a calculate phase field - xs = np.linspace( np.min(x), np.max(x), 50) - ys = np.linspace( np.min(y), np.max(y), 50) - phases, (xs, ys) = lib.phase_field_from_tx(xs, ys, tx, f_beacon, return_meshgrid=False) - sc2 = axs.scatter(xs, ys, c=phases, alpha=0.5, zorder=-5) - fig.colorbar(sc2, ax=axs) + if True: + # underlie a calculate phase field + if True: # only fill for antennas + xs = np.linspace( np.min(x), np.max(x), 4*np.ceil(len(antennas)**0.5) -1 ) + ys = np.linspace( np.min(y), np.max(y), 4*np.ceil(len(antennas)**0.5) -1) + 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) - sc = axs.scatter(x, y, c=vals, s=sizes) + 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, 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) - axs.plot(tx.x, tx.y, marker='X', color='k') + phase_diffs = vals - lib.phase_mod(phases) + phase_diffs = lib.phase_mod(phase_diffs) - fig.colorbar(sc, ax=axs, label=colorlabel) + 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()