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 11:33:32 +01:00
Code Issues Projects Releases Packages Wiki Activity

ZH: figures showing baseline diff reconstruction in time domain

Browse source
This commit is contained in:
Eric Teunis de Boone 2022-12-22 18:00:34 +01:00
parent 00d73175f7
commit 61b5b8a8d6
1 changed files with 36 additions and 22 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

48
simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py
View file

@ -160,11 +160,11 @@ if __name__ == "__main__":
############################## ##############################
# Compare actual time shifts # # Compare actual time shifts #
############################## ##############################
antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in sorted(antennas, key=lambda a: int(a.name)) } actual_antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in sorted(antennas, key=lambda a: int(a.name)) }
if True: if True:
actual_phase_shifts = [ -1*lib.phase_mod(2*np.pi*f_beacon*v) for k,v in antenna_time_shifts.items() ] actual_antenna_phase_shifts = [ -1*lib.phase_mod(2*np.pi*f_beacon*v) for k,v in actual_antenna_time_shifts.items() ]
antenna_names = [int(k)-1 for k,v in antenna_time_shifts.items() ] antenna_names = [int(k)-1 for k,v in actual_antenna_time_shifts.items() ]
for i in range(2): for i in range(2):
plot_residuals = i == 1 plot_residuals = i == 1
@ -179,7 +179,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_phase_shifts) phase_residuals = lib.phase_mod(mean_sigma_phase - actual_antenna_phase_shifts)
fig.suptitle("Difference between Measured and Actual phases\n for Antenna $i$") fig.suptitle("Difference between Measured and Actual phases\n for Antenna $i$")
axs[-1].set_xlabel("Antenna Phase Residual $\\Delta_\\varphi$") axs[-1].set_xlabel("Antenna Phase Residual $\\Delta_\\varphi$")
else: else:
@ -193,7 +193,7 @@ if __name__ == "__main__":
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_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') axs[i].hist(actual_antenna_phase_shifts, bins='sqrt', density=False, alpha=0.8, color=colors[1], ls='dashed', histtype='step', label='Actual')
i=1 i=1
@ -202,7 +202,7 @@ if __name__ == "__main__":
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_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].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()
fig.tight_layout() fig.tight_layout()
@ -217,31 +217,45 @@ if __name__ == "__main__":
########################## ##########################
########################## ##########################
actual_time_shifts = [] actual_baseline_time_shifts = []
for i,b in enumerate(basenames): 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_baseline_time_shift = actual_antenna_time_shifts[b[0]] - actual_antenna_time_shifts[b[1]]
actual_time_shifts.append(actual_time_shift) actual_baseline_time_shifts.append(actual_baseline_time_shift)
# unpack mean_sigma_phase back into a list of time diffs # unpack mean_sigma_phase back into a list of time diffs
measured_time_diffs = [] measured_baseline_time_diffs = []
for i,b in enumerate(basenames): for i,b in enumerate(basenames):
time0, time1 = mean_sigma_phase[name2idx(b[0])], mean_sigma_phase[name2idx(b[1])] phase0, phase1 = mean_sigma_phase[name2idx(b[0])], mean_sigma_phase[name2idx(b[1])]
measured_time_diffs.append(time1 - time0) measured_baseline_time_diffs.append(lib.phase_mod(phase1 - phase0)/(2*np.pi*f_beacon))
# Make a plot # Make a plot
if True: if True:
for i in range(2):
fig, ax = plt.subplots() fig, ax = plt.subplots()
ax.set_title("Baseline Time difference reconstruction" + ( '' if i == 0 else ' (wrapped time)'))
ax.set_xlabel("Baseline no.") ax.set_xlabel("Baseline no.")
ax.set_ylabel("$\\Delta t$[ns]") ax.set_ylabel("Time $\\Delta t$ [ns]")
if True:
forward = lambda x: x/(2*np.pi*f_beacon)
inverse = lambda x: 2*np.pi*x*f_beacon
secax = ax.secondary_yaxis('right', functions=(inverse, forward))
secax.set_ylabel('Phase $\\Delta \\varphi$ [rad]')
if True: # indicate single beacon period span 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((-1, -1), (-1/(2*f_beacon), 1/(2*f_beacon)), marker='3', ms=10, label='1/f_beacon')
ax.plot(np.arange(N_base), actual_time_shifts, marker='+', label='actual time shifts') if i == 0:
ax.plot(np.arange(N_base), measured_time_diffs, marker='x', label='calculated') ax.plot(np.arange(N_base), actual_baseline_time_shifts, marker='+', label='actual time shifts')
else:
ax.plot(np.arange(N_base), (actual_baseline_time_shifts+1/(2*f_beacon))%(1/f_beacon) - 1/(2*f_beacon), marker='+', label='actual time shifts')
ax.plot(np.arange(N_base), measured_baseline_time_diffs, marker='x', label='calculated')
ax.legend() ax.legend()
if fig_dir: if fig_dir:
fig.savefig(path.join(fig_dir, __file__ + f".calculated_shifts.pdf")) extra_name = ''
if i == 1:
extra_name = '.wrapped'
fig.savefig(path.join(fig_dir, __file__ + f".time_comparison{extra_name}.pdf"))
if show_plots: if show_plots:
plt.show() plt.show()