ericteunis/m-thesis-introduction
1
0
Fork 0
mirror of https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git synced 2024-12-22 03:23:34 +01:00
Code Issues Projects Releases Packages Wiki Activity

ZH: renaming phase variables III: sigma_phase_*->clock_phase_*

Browse source
This commit is contained in:
Eric Teunis de Boone 2023-01-19 17:13:29 +01:00
parent 6be1bb129f
commit 75999e6eb3
3 changed files with 26 additions and 26 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
simulations/airshower_beacon_simulation/aa_generate_beacon.py
View file

@ -33,7 +33,7 @@ def read_antenna_clock_repair_offsets(antennas, mode='all', freq_name=None):
# phase # phase
if mode in ['all', 'phases']: if mode in ['all', 'phases']:
clock_phase = ant.beacon_info[freq_name]['sigma_phase_mean'] clock_phase = ant.beacon_info[freq_name]['clock_phase_mean']
f_beacon = ant.beacon_info[freq_name]['freq'] f_beacon = ant.beacon_info[freq_name]['freq']
clock_phase_time = clock_phase/(2*np.pi*f_beacon) clock_phase_time = clock_phase/(2*np.pi*f_beacon)

48
simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py
View file

@ -48,19 +48,19 @@ if __name__ == "__main__":
N_ant = len(antennas) N_ant = len(antennas)
# reshape time_diffs into N_ant x N_ant matrix # reshape time_diffs into N_ant x N_ant matrix
sigma_phase_matrix = np.full( (N_ant, N_ant), np.nan, dtype=float) clock_phase_matrix = np.full( (N_ant, N_ant), np.nan, dtype=float)
## set i=i terms to 0 ## set i=i terms to 0
for i in range(N_ant): for i in range(N_ant):
sigma_phase_matrix[i,i] = 0 clock_phase_matrix[i,i] = 0
## fill matrix ## fill matrix
name2idx = lambda name: int(name)-1 name2idx = lambda name: int(name)-1
for i, b in enumerate(basenames): for i, b in enumerate(basenames):
idx = (name2idx(b[0]), name2idx(b[1])) idx = (name2idx(b[0]), name2idx(b[1]))
sigma_phase_matrix[(idx[0], idx[1])] = lib.phase_mod(clock_phase_diffs[i]) clock_phase_matrix[(idx[0], idx[1])] = lib.phase_mod(clock_phase_diffs[i])
sigma_phase_matrix[(idx[1], idx[0])] = lib.phase_mod(-1*clock_phase_diffs[i]) clock_phase_matrix[(idx[1], idx[0])] = lib.phase_mod(-1*clock_phase_diffs[i])
mat_kwargs = dict( mat_kwargs = dict(
norm = Normalize(vmin=-np.pi, vmax=+np.pi), norm = Normalize(vmin=-np.pi, vmax=+np.pi),
@ -74,7 +74,7 @@ if __name__ == "__main__":
ax.set_ylabel("Antenna no. i") ax.set_ylabel("Antenna no. i")
ax.set_xlabel("Antenna no. j") ax.set_xlabel("Antenna no. j")
im = ax.imshow(sigma_phase_matrix, interpolation='none', **mat_kwargs) im = ax.imshow(clock_phase_matrix, interpolation='none', **mat_kwargs)
fig.colorbar(im, ax=ax) fig.colorbar(im, ax=ax)
if fig_dir: if fig_dir:
@ -88,7 +88,7 @@ if __name__ == "__main__":
if True: if True:
# for each row j subtract the 0,j element from the whole row # for each row j subtract the 0,j element from the whole row
# and apply phase_mod # and apply phase_mod
first_row = -1*(sigma_phase_matrix[0,:] * np.ones_like(sigma_phase_matrix)).T first_row = -1*(clock_phase_matrix[0,:] * np.ones_like(clock_phase_matrix)).T
# Show subtraction Matrix as figure # Show subtraction Matrix as figure
if True: if True:
@ -105,18 +105,18 @@ if __name__ == "__main__":
plt.close(fig) plt.close(fig)
sigma_phase_matrix = sigma_phase_matrix - first_row clock_phase_matrix = clock_phase_matrix - first_row
sigma_phase_matrix = lib.phase_mod(sigma_phase_matrix) clock_phase_matrix = lib.phase_mod(clock_phase_matrix)
# Except for the first row, these are all separate measurements # Except for the first row, these are all separate measurements
# Condense into phase offset per antenna # Condense into phase offset per antenna
if True: # do not use the first row if True: # do not use the first row
my_mask = np.arange(1, len(sigma_phase_matrix), dtype=int) my_mask = np.arange(1, len(clock_phase_matrix), dtype=int)
else: else:
my_mask = np.arange(0, len(sigma_phase_matrix), dtype=int) my_mask = np.arange(0, len(clock_phase_matrix), dtype=int)
mean_sigma_phase = np.nanmean(sigma_phase_matrix[my_mask], axis=0) mean_clock_phase = np.nanmean(clock_phase_matrix[my_mask], axis=0)
std_sigma_phase = np.nanstd( sigma_phase_matrix[my_mask], axis=0) std_clock_phase = np.nanstd( clock_phase_matrix[my_mask], axis=0)
# Show resulting matrix as figure # Show resulting matrix as figure
@ -126,14 +126,14 @@ if __name__ == "__main__":
axs[0].set_ylabel("Antenna no. 0") axs[0].set_ylabel("Antenna no. 0")
axs[-1].set_xlabel("Antenna no. j") axs[-1].set_xlabel("Antenna no. j")
im = axs[0].imshow(sigma_phase_matrix, interpolation='none', **mat_kwargs) im = axs[0].imshow(clock_phase_matrix, interpolation='none', **mat_kwargs)
fig.colorbar(im, ax=axs) fig.colorbar(im, ax=axs)
axs[0].set_aspect('auto') axs[0].set_aspect('auto')
colours = [mat_kwargs['cmap'](mat_kwargs['norm'](x)) for x in mean_sigma_phase] colours = [mat_kwargs['cmap'](mat_kwargs['norm'](x)) for x in mean_clock_phase]
axs[1].set_ylabel("Mean Baseline Phase (0, j)[rad]") axs[1].set_ylabel("Mean Baseline Phase (0, j)[rad]")
axs[1].errorbar(np.arange(N_ant), mean_sigma_phase, yerr=std_sigma_phase, ls='none') axs[1].errorbar(np.arange(N_ant), mean_clock_phase, yerr=std_clock_phase, ls='none')
axs[1].scatter(np.arange(N_ant), mean_sigma_phase, c=colours,s=4) axs[1].scatter(np.arange(N_ant), mean_clock_phase, c=colours,s=4)
if fig_dir: if fig_dir:
fig.savefig(path.join(fig_dir, path.basename(__file__) + f".matrix.modified.pdf")) fig.savefig(path.join(fig_dir, path.basename(__file__) + f".matrix.modified.pdf"))
@ -157,8 +157,8 @@ if __name__ == "__main__":
h5attrs = h5beacon_info[freq_name].attrs h5attrs = h5beacon_info[freq_name].attrs
idx = name2idx(ant.name) idx = name2idx(ant.name)
h5attrs['sigma_phase_mean'] = mean_sigma_phase[idx] h5attrs['clock_phase_mean'] = mean_clock_phase[idx]
h5attrs['sigma_phase_std'] = std_sigma_phase[idx] h5attrs['clock_phase_std'] = std_clock_phase[idx]
############################## ##############################
@ -171,7 +171,7 @@ if __name__ == "__main__":
antenna_names = [int(k)-1 for k,v in actual_antenna_time_shifts.items() ] antenna_names = [int(k)-1 for k,v in actual_antenna_time_shifts.items() ]
# Make sure to shift all antennas by a global phase # Make sure to shift all antennas by a global phase
global_phase_shift = actual_antenna_phase_shifts[0] - mean_sigma_phase[0] global_phase_shift = actual_antenna_phase_shifts[0] - mean_clock_phase[0]
actual_antenna_phase_shifts = lib.phase_mod(actual_antenna_phase_shifts - global_phase_shift ) actual_antenna_phase_shifts = lib.phase_mod(actual_antenna_phase_shifts - global_phase_shift )
for i in range(2): for i in range(2):
@ -187,7 +187,7 @@ if __name__ == "__main__":
secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]')
if plot_residuals: if plot_residuals:
phase_residuals = lib.phase_mod(mean_sigma_phase - actual_antenna_phase_shifts) phase_residuals = lib.phase_mod(mean_clock_phase - actual_antenna_phase_shifts)
fig.suptitle("Difference between Measured and Actual phases (minus global phase)\n for Antenna $i$") fig.suptitle("Difference between Measured and Actual phases (minus global phase)\n for Antenna $i$")
axs[-1].set_xlabel("Antenna Phase Residual $\\Delta_\\varphi$") axs[-1].set_xlabel("Antenna Phase Residual $\\Delta_\\varphi$")
else: else:
@ -200,7 +200,7 @@ if __name__ == "__main__":
if plot_residuals: if plot_residuals:
axs[i].hist(phase_residuals, bins='sqrt', alpha=0.8, color=colors[0]) axs[i].hist(phase_residuals, bins='sqrt', alpha=0.8, color=colors[0])
else: 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(mean_clock_phase, bins='sqrt', density=False, alpha=0.8, color=colors[0], ls='solid' , histtype='step', label='Measured')
axs[i].hist(actual_antenna_phase_shifts, bins='sqrt', density=False, alpha=0.8, color=colors[1], ls='dashed', histtype='step', label='Actual') axs[i].hist(actual_antenna_phase_shifts, bins='sqrt', density=False, alpha=0.8, color=colors[1], ls='dashed', histtype='step', label='Actual')
@ -209,7 +209,7 @@ if __name__ == "__main__":
if plot_residuals: if plot_residuals:
axs[i].plot(phase_residuals, np.arange(N_ant), alpha=0.6, ls='none', marker='x', color=colors[0]) axs[i].plot(phase_residuals, np.arange(N_ant), alpha=0.6, ls='none', marker='x', color=colors[0])
else: 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].errorbar(mean_clock_phase, np.arange(N_ant), yerr=std_clock_phase, marker='4', alpha=0.7, ls='none', color=colors[0], label='Measured')
axs[i].plot(actual_antenna_phase_shifts, antenna_names, ls='none', marker='3', alpha=0.8, color=colors[1], label='Actual') axs[i].plot(actual_antenna_phase_shifts, antenna_names, ls='none', marker='3', alpha=0.8, color=colors[1], label='Actual')
axs[i].legend() axs[i].legend()
@ -231,10 +231,10 @@ if __name__ == "__main__":
actual_baseline_time_shifts.append(actual_baseline_time_shift) actual_baseline_time_shifts.append(actual_baseline_time_shift)
# unpack mean_sigma_phase back into a list of time diffs # unpack mean_clock_phase back into a list of time diffs
measured_baseline_time_diffs = [] measured_baseline_time_diffs = []
for i,b in enumerate(basenames): for i,b in enumerate(basenames):
phase0, phase1 = mean_sigma_phase[name2idx(b[0])], mean_sigma_phase[name2idx(b[1])] phase0, phase1 = mean_clock_phase[name2idx(b[0])], mean_clock_phase[name2idx(b[1])]
measured_baseline_time_diffs.append(lib.phase_mod(phase1 - phase0)/(2*np.pi*f_beacon)) measured_baseline_time_diffs.append(lib.phase_mod(phase1 - phase0)/(2*np.pi*f_beacon))
# Make a plot # Make a plot

2
simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py
View file

@ -51,7 +51,7 @@ if __name__ == "__main__":
# TODO: redo matrix sweeping for new timing?? # TODO: redo matrix sweeping for new timing??
measured_antenna_time_shifts = {} measured_antenna_time_shifts = {}
for i, ant in enumerate(antennas): for i, ant in enumerate(antennas):
clock_phase_time = ant.beacon_info[freq_name]['sigma_phase_mean']/(2*np.pi*f_beacon) clock_phase_time = ant.beacon_info[freq_name]['clock_phase_mean']/(2*np.pi*f_beacon)
best_k_time = ant.beacon_info[freq_name]['best_k_time'] best_k_time = ant.beacon_info[freq_name]['best_k_time']
total_clock_time = best_k_time + clock_phase_time total_clock_time = best_k_time + clock_phase_time