#!/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') fname = "ZH_airshower/mysim.sry" show_plots = True ref_ant_id = None # leave None for all baselines #### fname_dir = path.dirname(fname) antennas_fname = path.join(fname_dir, beacon.antennas_fname) time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname fig_dir = "./figures" # set None to disable saving basenames, time_diffs, f_beacons, true_phase_diffs, k_periods = beacon.read_baseline_time_diffs_hdf5(time_diffs_fname) f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) # TODO: allow multiple frequencies if (f_beacon != f_beacons).any(): raise NotImplementedError N_base = len(basenames) N_ant = len(antennas) # reshape time_diffs into N_ant x N_ant array sigma_phase_matrix = np.full( (N_ant, N_ant), np.nan, dtype=float) #sigma_phase_matrix = np.arange(0, N_ant**2, dtype=float).reshape( (N_ant,N_ant)) name2idx = lambda name: int(name)-1 for i, b in enumerate(basenames): idx = (name2idx(b[0]), name2idx(b[1])) if idx[0] == idx[1]: pass sigma_phase_matrix[(idx[0], idx[1])] = true_phase_diffs[i] sigma_phase_matrix[(idx[1], idx[0])] = true_phase_diffs[i] # for each row j subtract the 0,j element from the whole row # and apply phase_mod first_row = sigma_phase_matrix[0] sigma_phase_matrix = sigma_phase_matrix - first_row[:,np.newaxis] sigma_phase_matrix = lib.phase_mod(sigma_phase_matrix) # Except for the first row, these are all separate measurements # Condense into phase offset per antenna if True: # do not use the first row my_mask = np.arange(1, len(sigma_phase_matrix), dtype=int) else: my_mask = np.arange(0, len(sigma_phase_matrix), dtype=int) mean_sigma_phase = np.nanmean(sigma_phase_matrix[my_mask], axis=0) std_sigma_phase = np.nanstd( sigma_phase_matrix[my_mask], axis=0) # write into antenna hdf5 with h5py.File(antennas_fname, 'a') as fp: group = fp['antennas'] freq_name = None for i, ant in enumerate(antennas): h5ant = group[ant.name] h5beacon_info = h5ant['beacon_info'] # find out freq_name if freq_name is None: freq_name = [ k for k in h5beacon_info.keys() if np.isclose(h5beacon_info[k].attrs['freq'], f_beacon)][0] h5attrs = h5beacon_info[freq_name].attrs idx = name2idx(ant.name) h5attrs['sigma_phase_mean'] = mean_sigma_phase[idx] h5attrs['sigma_phase_std'] = std_sigma_phase[idx] ############################## # Compare actual time shifts # ############################## antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in sorted(antennas, key=lambda a: int(a.name)) } if True: actual_phase_shifts = [ -1*lib.phase_mod(2*np.pi*f_beacon*v) for k,v in antenna_time_shifts.items() ] antenna_names = [int(k)-1 for k,v in antenna_time_shifts.items() ] for i in range(2): plot_residuals = i == 1 colors = ['blue', 'orange'] fig, axs = plt.subplots(2, 1, sharex=True) if True: forward = lambda x: x/(2*np.pi*f_beacon) inverse = lambda x: 2*np.pi*x*f_beacon secax = axs[0].secondary_xaxis('top', functions=(forward, inverse)) secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') if plot_residuals: phase_residuals = lib.phase_mod(mean_sigma_phase - actual_phase_shifts) fig.suptitle("Difference between Measured and Actual phases\n for Antenna $i$") axs[-1].set_xlabel("Antenna Phase Residual $\\Delta_\\varphi$") else: fig.suptitle("Comparison Measured and Actual phases\n for Antenna $i$") axs[-1].set_xlabel("Antenna Phase $\\Delta_\\varphi$") i=0 axs[i].set_ylabel("#") if plot_residuals: axs[i].hist(phase_residuals, bins='sqrt', alpha=0.8, color=colors[0]) else: axs[i].hist(mean_sigma_phase, bins='sqrt', density=False, alpha=0.8, color=colors[0], ls='solid' , histtype='step', label='Measured') axs[i].hist(actual_phase_shifts, bins='sqrt', density=False, alpha=0.8, color=colors[1], ls='dashed', histtype='step', label='Actual') i=1 axs[i].set_ylabel("Antenna no.") if plot_residuals: axs[i].plot(phase_residuals, np.arange(N_ant), alpha=0.6, ls='none', marker='x', color=colors[0]) else: axs[i].errorbar(mean_sigma_phase, np.arange(N_ant), yerr=std_sigma_phase, marker='4', alpha=0.7, ls='none', color=colors[0], label='Measured') axs[i].plot(actual_phase_shifts, antenna_names, ls='none', marker='3', alpha=0.8, color=colors[1], label='Actual') axs[i].legend() fig.tight_layout() if fig_dir: extra_name = "measured" if plot_residuals: extra_name = "residuals" fig.savefig(path.join(fig_dir, __file__ + f".phase.{extra_name}.pdf")) ########################## ########################## ########################## actual_time_shifts = [] for i,b in enumerate(basenames): actual_time_shift = lib.phase_mod(lib.phase_mod(antenna_time_shifts[b[0]]*2*np.pi*f_beacon) - lib.phase_mod(antenna_time_shifts[b[1]]*2*np.pi*f_beacon)) actual_time_shifts.append(actual_time_shift) # unpack mean_sigma_phase back into a list of time diffs measured_time_diffs = [] for i,b in enumerate(basenames): time0, time1 = mean_sigma_phase[name2idx(b[0])], mean_sigma_phase[name2idx(b[1])] measured_time_diffs.append(time1 - time0) # Make a plot if True: fig, ax = plt.subplots() ax.set_xlabel("Baseline no.") ax.set_ylabel("$\\Delta t$[ns]") if True: # indicate single beacon period span ax.plot((-1, -1), (0, 1/f_beacon), marker='3', ms=10, label='1/f_beacon') ax.plot(np.arange(N_base), actual_time_shifts, marker='+', label='actual time shifts') ax.plot(np.arange(N_base), measured_time_diffs, marker='x', label='calculated') ax.legend() if fig_dir: fig.savefig(path.join(fig_dir, __file__ + f".calculated_shifts.pdf")) if show_plots: plt.show()