#!/usr/bin/env python3 # vim: fdm=indent ts=4 import h5py from itertools import combinations, zip_longest import matplotlib.pyplot as plt import numpy as np import aa_generate_beacon as beacon import lib if __name__ == "__main__": from os import path import sys import os import matplotlib if os.name == 'posix' and "DISPLAY" not in os.environ: matplotlib.use('Agg') from scriptlib import MyArgumentParser parser = MyArgumentParser() args = parser.parse_args() figsize = (12,8) c_light = lib.c_light show_plots = args.show_plots remove_absolute_phase_offset_first_antenna = True # takes precedence remove_absolute_phase_offset_minimum = True #### fname_dir = args.data_dir antennas_fname = path.join(fname_dir, beacon.antennas_fname) fig_dir = args.fig_dir # set None to disable saving if not path.isfile(antennas_fname): print("Antenna file cannot be found, did you try generating a beacon?") sys.exit(1) # Read in antennas from file f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) # Make sure at least one beacon has been identified if not hasattr(antennas[0], 'beacon_info') or len(antennas[0].beacon_info) == 0: print(f"No analysed beacon found for {antennas[0].name}, try running the beacon phase analysis script first.") sys.exit(1) # N_beacon_freqs = len(antennas[0].beacon_info) for freq_name in antennas[0].beacon_info.keys(): beacon_phases = np.empty( (len(antennas)) ) for i, ant in enumerate(antennas): beacon_phases[i] = ant.beacon_info[freq_name]['beacon_phase'] f_beacon = antennas[0].beacon_info[freq_name]['freq'] clock_phases = lib.remove_antenna_geometry_phase(tx, antennas, f_beacon, beacon_phases, c_light=c_light) # Remove the phase from one antenna # this is a free parameter # (only required for absolute timing) 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 = clock_phases[0] else: # take the minimum minimum_phase = np.min(clock_phases, axis=-1) clock_phases -= minimum_phase clock_phases = lib.phase_mod(clock_phases) # Save to antennas in file with h5py.File(antennas_fname, 'a') as fp: h5group = fp['antennas'] for i, ant in enumerate(antennas): h5ant = fp['antennas'][ant.name] h5beacon_freq = h5ant['beacon_info'][freq_name] h5beacon_freq.attrs['clock_phase'] = clock_phases[i] # Plot True Phases at their locations if show_plots or fig_dir: 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])) scatter_kwargs = {} scatter_kwargs['cmap'] = 'inferno' # Measurements if not plot_residuals: title='Clock phases' color_label='$\\varphi(\\sigma_t)$ [rad]' fname_extra='measured.' #scatter_kwargs['vmin'] = -np.pi #scatter_kwargs['vmax'] = +np.pi # 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]' 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) 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) if not plot_residuals: loc_c = clock_phases else: loc_c = clock_phase_residuals ## ## Geometrical Plot ## fig, ax = plt.subplots(figsize=figsize) 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) if show_plots: plt.show()