From 91016be038c69f23c52ab063a1fba37d92018c47 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Wed, 1 Feb 2023 17:29:50 +0100 Subject: [PATCH] ZH: view_beacon: plot noise and beacon traces --- .../view_beaconed_antenna.py | 44 ++++++++++++++----- 1 file changed, 34 insertions(+), 10 deletions(-) diff --git a/simulations/airshower_beacon_simulation/view_beaconed_antenna.py b/simulations/airshower_beacon_simulation/view_beaconed_antenna.py index 06990f8..5af931e 100755 --- a/simulations/airshower_beacon_simulation/view_beaconed_antenna.py +++ b/simulations/airshower_beacon_simulation/view_beaconed_antenna.py @@ -21,7 +21,7 @@ if __name__ == "__main__": from scriptlib import MyArgumentParser parser = MyArgumentParser() parser.add_argument('ant_idx', default=[72], nargs='*', type=int, help='Antenna Indices') - parser.add_argument('-p', '--polarisations', choices=['x', 'y', 'z', 'b', 'AxB'], action='append', help='Default: x,y,z') + parser.add_argument('-p', '--polarisations', choices=['x', 'y', 'z', 'b', 'AxB', 'n', 'b+n'], action='append', help='Default: x,y,z') parser.add_argument('--geom', action='store_true', help='Make a figure containg the geometry from tx to antenna(s)') parser.add_argument('--ft', action='store_true', help='Add FT strenghts of antenna traces') @@ -42,8 +42,11 @@ if __name__ == "__main__": #### fname_dir = path.dirname(fname) antennas_fname = path.join(fname_dir, beacon.antennas_fname) + tx_fname = path.join(fname_dir, beacon.tx_fname) f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) + _, __, txdata = beacon.read_tx_file(tx_fname) + beacon_amp = np.max(txdata['amplitudes'])# mu V/m idx = args.ant_idx @@ -66,7 +69,7 @@ if __name__ == "__main__": name_dist='.raw' - fig1, axs = plt.subplots(1+plot_ft_amplitude*2, figsize=figsize) + fig1, axs = plt.subplots(1+plot_ft_amplitude*1 +0*1, figsize=figsize) if not plot_ft_amplitude: axs = [axs] axs[0].set_xlabel('t [ns]') @@ -84,9 +87,10 @@ if __name__ == "__main__": axs[1].set_xlabel('f [GHz]') axs[1].set_ylabel('Power') - axs[2].set_ylabel("Phase") - axs[2].set_xlabel('f [GHz]') - axs[2].set_ylim(-np.pi,+np.pi) + if len(axs) > 2: + axs[2].set_ylabel("Phase") + axs[2].set_xlabel('f [GHz]') + axs[2].set_ylim(-np.pi,+np.pi) colorlist = [] for i in idx: @@ -102,11 +106,19 @@ if __name__ == "__main__": pattr = 'E'+str(p) if p == 'b': pattr = 'beacon' + elif p == 'n': + pattr = 'noise' elif p == 'AxB': pattr = 'E_AxB' + elif p =='b+n': + mydict[p] = getattr(ant,'noise') + beacon_amp*getattr(ant, 'beacon') + continue mydict[p] = getattr(ant, pattr) + if 'b' in mydict: + mydict['b'] *= beacon_amp + for j, (direction, trace) in enumerate(mydict.items()): l = axs[0].plot(ant.t, trace, label=f"$E_{{{direction}}}$ {ant.name}", alpha=0.7) @@ -118,26 +130,38 @@ if __name__ == "__main__": if not plot_ft_amplitude: continue - freqs = ft.fftfreq(n_samples, 1/samplerate)[:n_samples//2] - fft = 2*ft.fft(trace)[:n_samples//2]/n_samples + fft, freqs = lib.get_freq_spec(trace, 1/samplerate) - #axs[1].plot(freqs, np.abs(fft)**2, color=l[0].get_color()) + axs[1].plot(freqs, np.abs(fft)**2, color=l[0].get_color()) if True: cft = lib.direct_fourier_transform(f_beacon, ant.t, trace) - amp = 2*len(ant.t) * (cft[0]**2 + cft[1]**2) + amp = (cft[0]**2 + cft[1]**2) #axs[0].axhline(amp, color=l[0].get_color()) print(amp) phase = np.arctan2(cft[0],cft[1]) axs[1].plot(f_beacon, amp, color=l[0].get_color(), marker='3', alpha=0.8, ms=30) - axs[2].plot(f_beacon, phase, color=l[0].get_color(), marker='3', alpha=0.8, ms=30) + if len(axs) > 2: + axs[2].plot(f_beacon, phase, color=l[0].get_color(), marker='3', alpha=0.8, ms=30) if plot_ft_amplitude: fig1.legend(loc='center right', ncol=min(2, len(idx))) else: axs[0].legend(loc='upper right', ncol=min(3, len(idx))) + + # Keep trace plot symmetric around 0 + max_lim = max(np.abs(axs[0].get_ylim())) + axs[0].set_ylim(-max_lim, max_lim) + + # Keep spectrum between 0 and 100 MHz + if len(axs) > 1: + xlims = axs[1].get_xlim() + axs[1].set_xlim(max(0, xlims[0]), min(0.1, xlims[1])) + if False: # extra zoom + axs[1].set_xlim(f_beacon - 0.01, f_beacon + 0.01) + if fig_dir: fig1.savefig(path.join(fig_dir, path.basename(__file__) + f".trace{name_dist}.pdf"))