diff --git a/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py b/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py index a83f717..4295b28 100755 --- a/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py +++ b/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py @@ -83,64 +83,119 @@ if __name__ == "__main__": # Plot True Phases at their locations if show_plots or fig_dir: - fig, ax = plt.subplots(figsize=figsize) - spatial_unit='m' - fig.suptitle('Clock phases\nf_beacon= {:2.0f}MHz'.format(f_beacon*1e3)) + actual_clock_phases = lib.phase_mod(np.array([ -2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas ])) + for i in range(2): + plot_residuals = i == 1 + spatial_unit='m' - antenna_locs = list(zip(*[(ant.x, ant.y) for ant in antennas])) - ax.set_xlabel('x' if spatial_unit is None else 'x [{}]'.format(spatial_unit)) - ax.set_ylabel('y' if spatial_unit is None else 'y [{}]'.format(spatial_unit)) - scatter_kwargs = {} - scatter_kwargs['cmap'] = 'inferno' - #scatter_kwargs['vmin'] = -np.pi - #scatter_kwargs['vmax'] = +np.pi - color_label='$\\varphi(\\sigma_t)$ [rad]' + antenna_locs = list(zip(*[(ant.x, ant.y) for ant in antennas])) - sc = ax.scatter(*antenna_locs, c=clock_phases, **scatter_kwargs) - fig.colorbar(sc, ax=ax, label=color_label) + scatter_kwargs = {} + scatter_kwargs['cmap'] = 'inferno' - if False: - for i, ant in enumerate(antennas): - ax.text(ant.x, ant.y, ant.name) + # Measurements + if not plot_residuals: + title='Clock phases' + color_label='$\\varphi(\\sigma_t)$ [rad]' - if not True: - ax.plot(tx.x, tx.y, 'X', color='k', markeredgecolor='white') + fname_extra='measured.' - if fig_dir: - fig.savefig(path.join(fig_dir, path.basename(__file__) + f".F{freq_name}.pdf")) + #scatter_kwargs['vmin'] = -np.pi + #scatter_kwargs['vmax'] = +np.pi - # Plot True Phases - Actual True Phases at their location - if show_plots or fig_dir: - fig, ax = plt.subplots(figsize=figsize) - fig.suptitle('Clock phase Residuals\nf_beacon={:2.0f}MHz'.format(f_beacon*1e3)) + # Plot Clock Phases - True Clock Phases at their location + else: + title='Clock phase Residuals' + color_label='$\\Delta\\varphi(\\sigma_t) = \\varphi_{meas} - \\varphi_{true}$ [rad]' - actual_clock_phases = np.array([ -2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas ]) + fname_extra='residuals.' - # Modify actual_clock_phases, the same way as clock_phases - # was modified - if remove_absolute_phase_offset_first_antenna or remove_absolute_phase_offset_minimum: - if remove_absolute_phase_offset_first_antenna: # just take the first phase - minimum_phase = actual_clock_phases[0] - else: # take the minimum - minimum_phase = np.min(actual_clock_phases, axis=-1) + # Modify actual_clock_phases, the same way as clock_phases + # was modified + if remove_absolute_phase_offset_first_antenna or remove_absolute_phase_offset_minimum: + if remove_absolute_phase_offset_first_antenna: # just take the first phase + minimum_phase = actual_clock_phases[0] + else: # take the minimum + minimum_phase = np.min(actual_clock_phases, axis=-1) - actual_clock_phases -= minimum_phase - actual_clock_phases = lib.phase_mod(actual_clock_phases) + actual_clock_phases -= minimum_phase + actual_clock_phases = lib.phase_mod(actual_clock_phases) - clock_phase_residuals = lib.phase_mod(clock_phases - actual_clock_phases) + clock_phase_residuals = lib.phase_mod(clock_phases - actual_clock_phases) - antenna_locs = list(zip(*[(ant.x, ant.y) for ant in antennas])) - ax.set_xlabel('x' if spatial_unit is None else 'x [{}]'.format(spatial_unit)) - ax.set_ylabel('y' if spatial_unit is None else 'y [{}]'.format(spatial_unit)) - scatter_kwargs = {} - scatter_kwargs['cmap'] = 'inferno' - color_label='$\\Delta\\varphi(\\sigma_t) = \\varphi_{meas} - \\varphi_{true}$ [rad]' + if not plot_residuals: + loc_c = clock_phases + else: + loc_c = clock_phase_residuals - sc = ax.scatter(*antenna_locs, c=clock_phase_residuals, **scatter_kwargs) - fig.colorbar(sc, ax=ax, label=color_label) + ## + ## Geometrical Plot + ## + fig, ax = plt.subplots(figsize=figsize) - if fig_dir: - fig.savefig(path.join(fig_dir, path.basename(__file__) + f".residual.F{freq_name}.pdf")) + ax.set_xlabel('x' if spatial_unit is None else 'x [{}]'.format(spatial_unit)) + ax.set_ylabel('y' if spatial_unit is None else 'y [{}]'.format(spatial_unit)) + + fig.suptitle(title+'\nf_beacon= {:2.0f}MHz'.format(f_beacon*1e3)) + + sc = ax.scatter(*antenna_locs, c=loc_c, **scatter_kwargs) + fig.colorbar(sc, ax=ax, label=color_label) + + if False: + for i, ant in enumerate(antennas): + ax.text(ant.x, ant.y, ant.name) + + if not True: + ax.plot(tx.x, tx.y, 'X', color='k', markeredgecolor='white') + + if fig_dir: + fig.tight_layout() + fig.savefig(path.join(fig_dir, path.basename(__file__) + f".geom.{fname_extra}F{freq_name}.pdf")) + + ## + ## Histogram + ## + fig, axs = plt.subplots(2, 1, sharex=True, figsize=figsize) + colors = ['blue', 'orange'] + + if True: + phase2time = lambda x: x/(2*np.pi*f_beacon) + time2phase = lambda x: 2*np.pi*x*f_beacon + secax = axs[0].secondary_xaxis('top', functions=(phase2time, time2phase)) + secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') + + if plot_residuals: + fig.suptitle("Difference between Measured and True Clock phases") + else: + fig.suptitle("Comparison Measured and True Clock Phases") + + axs[-1].set_xlabel(f'Antenna {title} {color_label}') + + # + hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step') + + i=0 + axs[i].set_ylabel("#") + axs[i].hist(loc_c, color=colors[0], label='Measured', ls='solid', **hist_kwargs) + if not plot_residuals: # also plot the true clock phases + axs[i].hist(actual_clock_phases, color=colors[1], label='Actual', ls='dashed', **hist_kwargs) + #axs[i].legend() + + # + plot_kwargs = dict(alpha=0.6, ls='none') + + i=1 + axs[i].set_ylabel("Antenna no.") + axs[i].plot(loc_c, np.arange(len(loc_c)), marker='x' if plot_residuals else '3', color=colors[0], label='Measured', **plot_kwargs) + + if not plot_residuals: # also plot the true clock phases + axs[i].plot(actual_clock_phases, np.arange(len(loc_c)), marker='4', color=colors[1], label='Actual', **plot_kwargs) + axs[i].legend() + + # Save figure + if fig_dir: + fig.tight_layout() + fig.savefig(path.join(fig_dir, path.basename(__file__) + f".{fname_extra}F{freq_name}.pdf")) print(f"True phases written to", antennas_fname)