m-thesis-introduction/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py

149 lines
5.7 KiB
Python
Raw Permalink Normal View History

#!/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
2022-12-08 14:41:33 +01:00
import os
2022-12-05 17:48:58 +01:00
import matplotlib
2022-12-08 14:41:33 +01:00
if os.name == 'posix' and "DISPLAY" not in os.environ:
matplotlib.use('Agg')
2022-12-05 17:48:58 +01:00
from scriptlib import MyArgumentParser
parser = MyArgumentParser()
args = parser.parse_args()
2023-02-02 08:57:03 +01:00
figsize = (12,8)
2022-12-15 14:40:20 +01:00
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)
2023-01-12 14:49:54 +01:00
fig_dir = args.fig_dir # set None to disable saving
2022-12-05 17:48:58 +01:00
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)) )
2022-12-15 12:56:20 +01:00
for i, ant in enumerate(antennas):
beacon_phases[i] = ant.beacon_info[freq_name]['beacon_phase']
2022-12-15 12:56:20 +01:00
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
2022-12-05 17:48:58 +01:00
if show_plots or fig_dir:
2023-02-01 14:13:26 +01:00
fig, ax = plt.subplots(figsize=figsize)
spatial_unit='m'
fig.suptitle('Clock phases\nf_beacon= {:2.0f}MHz'.format(f_beacon*1e3))
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]'
sc = ax.scatter(*antenna_locs, c=clock_phases, **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')
2022-12-05 17:48:58 +01:00
if fig_dir:
fig.savefig(path.join(fig_dir, path.basename(__file__) + f".F{freq_name}.pdf"))
2022-12-05 17:48:58 +01:00
# Plot True Phases - Actual True Phases at their location
if show_plots or fig_dir:
2023-02-01 14:13:26 +01:00
fig, ax = plt.subplots(figsize=figsize)
fig.suptitle('Clock phase Residuals\nf_beacon={:2.0f}MHz'.format(f_beacon*1e3))
actual_clock_phases = np.array([ -2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas ])
# 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)
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]'
sc = ax.scatter(*antenna_locs, c=clock_phase_residuals, **scatter_kwargs)
fig.colorbar(sc, ax=ax, label=color_label)
if fig_dir:
fig.savefig(path.join(fig_dir, path.basename(__file__) + f".residual.F{freq_name}.pdf"))
2022-12-05 17:48:58 +01:00
print(f"True phases written to", antennas_fname)
if show_plots:
2022-11-25 12:12:04 +01:00
plt.show()