ZH: improve some plots

This commit is contained in:
Eric Teunis de Boone 2022-11-22 15:37:55 +01:00
parent 9debde9d02
commit f6b474373d
3 changed files with 50 additions and 22 deletions

View file

@ -134,7 +134,7 @@ if __name__ == "__main__":
myt = np.linspace(min(traces[0]), max(traces[0]), 10*len(traces[0])) myt = np.linspace(min(traces[0]), max(traces[0]), 10*len(traces[0]))
ax.plot(t_trace, traces[-1], marker='.', label='trace') ax.plot(t_trace, traces[-1], marker='.', label='trace')
ax.plot(myt, lib.sine_beacon(frequency, myt, amplitude=amplitude, phase=phase), ls='dashed', label='simulated') ax.plot(myt, lib.sine_beacon(frequency, myt, amplitude=amplitude, phase=phase), ls='dashed', label='simulated')
ax.set_title(f"Beacon at antenna {h5ant}\nF:{frequency}, P:{phase}, A:{amplitude}") ax.set_title(f"Beacon at antenna {h5ant.attrs['name']}\nF:{frequency:.2e}, P:{phase:.4f}, A:{amplitude:.1e}")
ax.legend() ax.legend()
# save to file # save to file

View file

@ -18,17 +18,24 @@ if __name__ == "__main__":
fname = "ZH_airshower/mysim.sry" fname = "ZH_airshower/mysim.sry"
plot_phase_field = True
plot_tx = False
#### ####
fname_dir = path.dirname(fname) fname_dir = path.dirname(fname)
antennas_fname = path.join(fname_dir, beacon.antennas_fname) antennas_fname = path.join(fname_dir, beacon.antennas_fname)
f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
subtitle = "" pretitle = ""
if True: if True:
beacon_frequencies = np.array([ant.attrs['beacon_freq'] for ant in antennas]) freq_name = list(antennas[0].beacon_info.keys())[0]
beacon_amplitudes = np.array([ant.attrs['beacon_amplitude'] for ant in antennas]) beacon_frequencies = np.array([ant.beacon_info[freq_name]['freq'] for ant in antennas])
beacon_phases = np.array([lib.phase_mod(ant.attrs['beacon_phase_measured']) for ant in antennas]) beacon_amplitudes = np.array([ant.beacon_info[freq_name]['amplitude'] for ant in antennas])
beacon_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['phase'] for ant in antennas]))
if 'true_phase' in antennas[0].beacon_info[freq_name]:
beacon_true_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['true_phase'] for ant in antennas]))
else: else:
subtitle = " Phases from t0" subtitle = " Phases from t0"
beacon_frequencies = np.array([ f_beacon for ant in antennas ]) beacon_frequencies = np.array([ f_beacon for ant in antennas ])
@ -41,6 +48,11 @@ if __name__ == "__main__":
vals = beacon_phases vals = beacon_phases
colorlabel = '$\\varphi$' colorlabel = '$\\varphi$'
sizes = 64*(beacon_amplitudes/np.max(beacon_amplitudes))**2 sizes = 64*(beacon_amplitudes/np.max(beacon_amplitudes))**2
elif True: # True Phases
vals = beacon_true_phases
colorlabel = '$\\sigma_\\varphi$'
plot_phase_field = False
plot_tx = False
else: else:
vals = beacon_amplitudes vals = beacon_amplitudes
colorlabel = "[$\\mu$V/m]" colorlabel = "[$\\mu$V/m]"
@ -51,13 +63,14 @@ if __name__ == "__main__":
##### #####
fig, axs = plt.subplots() fig, axs = plt.subplots()
axs.set_title(f"Amplitude at beacon frequency at each antenna\nf at A0: {beacon_frequencies[0]}GHz" + subtitle) axs.set_title(pretitle + f"Beacon frequency at A0: {beacon_frequencies[0]:.3e}GHz")
axs.set_aspect('equal', 'datalim') axs.set_aspect('equal', 'datalim')
axs.set_xlabel('[m]') axs.set_xlabel('[m]')
axs.set_ylabel('[m]') axs.set_ylabel('[m]')
if True: if True:
if True: vmax, vmin = None, None
if plot_phase_field:
# underlie a calculate phase field # underlie a calculate phase field
if True: # only fill for antennas if True: # only fill for antennas
xs = np.linspace( np.min(x), np.max(x), 4*np.ceil(len(antennas)**0.5) -1 ) xs = np.linspace( np.min(x), np.max(x), 4*np.ceil(len(antennas)**0.5) -1 )
@ -67,11 +80,18 @@ if __name__ == "__main__":
ys = np.linspace( (tx.y - np.min(y))/2, np.max(y), 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) 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')
vmax, vmin = max(np.max(phases), np.max(vals)), min(np.min(phases), np.min(vals))
sc2 = axs.scatter(xs, ys, c=phases, alpha=0.5, zorder=-5, cmap='Spectral_r', vmin=vmin, vmax=vmax)
if False:
# do not share the same colours
fig.colorbar(sc2, ax=axs) fig.colorbar(sc2, ax=axs)
vmax, vmin = None, None
sc = axs.scatter(x, y, c=vals, s=sizes, cmap='Spectral_r', edgecolors='k', marker='X') sc = axs.scatter(x, y, c=vals, s=sizes, cmap='Spectral_r', edgecolors='k', marker='X', vmin=vmin, vmax=vmax)
if plot_tx:
axs.plot(tx.x, tx.y, marker='X', color='k') axs.plot(tx.x, tx.y, marker='X', color='k')
#for i, freq in enumerate(beacon_frequencies): #for i, freq in enumerate(beacon_frequencies):
@ -91,6 +111,6 @@ if __name__ == "__main__":
fig.colorbar(sc, ax=axs, label=colorlabel) fig.colorbar(sc, ax=axs, label=colorlabel)
if not True:
fig.savefig(__file__ + ".pdf") fig.savefig(__file__ + ".pdf")
plt.show() plt.show()

View file

@ -58,13 +58,15 @@ if __name__ == "__main__":
axs[0].axvline(ant.t[0], color='k', alpha=0.5) axs[0].axvline(ant.t[0], color='k', alpha=0.5)
if not True: if True: # total E field
mydict = dict(AxB=ant.E_AxB)
elif False: # polarisations
mydict = dict(x=ant.Ex, y=ant.Ex, z=ant.Ez) mydict = dict(x=ant.Ex, y=ant.Ex, z=ant.Ez)
else: else: # beacon
mydict = dict(E=ant.E_AxB, b=ant.beacon) mydict = dict(b=ant.beacon)
for j, (direction, trace) in enumerate(mydict.items()): for j, (direction, trace) in enumerate(mydict.items()):
l = axs[0].plot(ant.t, trace, label=f"E{direction} {ant.name}") l = axs[0].plot(ant.t, trace, label=f"$E_{{{direction}}}$ {ant.name}", alpha=0.8)
#if False and j == 0 and 't0' in ant.attrs: #if False and j == 0 and 't0' in ant.attrs:
# axs[0].axvline(ant.attrs['t0'], color=l[0].get_color(), alpha=0.5) # axs[0].axvline(ant.attrs['t0'], color=l[0].get_color(), alpha=0.5)
@ -87,19 +89,25 @@ if __name__ == "__main__":
print(amp) print(amp)
phase = np.arctan2(cft[0],cft[1]) phase = np.arctan2(cft[0],cft[1])
axs[1].plot(f_beacon, amp, color=l[0].get_color(), marker='3') 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') axs[2].plot(f_beacon, phase, color=l[0].get_color(), marker='3', alpha=0.8, ms=30)
if plot_ft_amplitude: if plot_ft_amplitude:
fig1.legend(loc='center right', ncol=min(3, len(idx))) fig1.legend(loc='center right', ncol=min(2, len(idx)))
else: else:
fig1.legend(loc='upper right', ncol=min(3, len(idx))) fig1.legend(loc='upper right', ncol=min(3, len(idx)))
if plot_geometry: if plot_geometry:
if len(mydict) == 1:
geom_colorlist = colorlist
else:
# only take the colour belonging to mydict[0]
geom_colorlist = [ colorlist[len(mydict)*(i)] for i in range(len(colorlist)//len(mydict)) ]
fig2, axs2 = plt.subplots(1) fig2, axs2 = plt.subplots(1)
plot_antenna_geometry(antennas, ax=axs2, plot_max_values=False, color='grey', plot_names=False) plot_antenna_geometry(antennas, ax=axs2, plot_max_values=False, color='grey', plot_names=False)
for j, _ in enumerate(mydict): plot_antenna_geometry([ antennas[i] for i in idx], ax=axs2, colors=geom_colorlist, plot_max_values=False)
plot_antenna_geometry([ antennas[i] for i in idx], ax=axs2, colors=colorlist[j + len(colorlist)//len(mydict)], plot_max_values=False)
axs2.plot(tx.x, tx.y, marker='X', color='k') axs2.plot(tx.x, tx.y, marker='X', color='k')
axs2.set_title("Geometry with selected antennas") axs2.set_title("Geometry with selected antennas")