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: update baseline phase script

Browse source
This commit is contained in:
Eric Teunis de Boone 2022-12-19 21:34:49 +01:00
parent ecb671bee8
commit 929c6c7748
1 changed files with 46 additions and 37 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

83
simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py
View file

@ -21,7 +21,7 @@ if __name__ == "__main__":
fname = "ZH_airshower/mysim.sry"
c_light = 3e8*1e-9
show_plots = not True
show_plots = True
ref_ant_id = None # leave None for all baselines
####
@ -66,7 +66,6 @@ if __name__ == "__main__":
# Get true phase diffs
try:
true_phases = np.array([ant.beacon_info[freq_name]['true_phase'] for ant in base])
true_phases_diff = lib.phase_mod(lib.phase_mod(true_phases[1]) - lib.phase_mod(true_phases[0]))
except IndexError:
@ -79,59 +78,69 @@ if __name__ == "__main__":
beacon.write_baseline_time_diffs_hdf5(time_diffs_fname, baselines, phase_diffs[:,1], [0]*len(phase_diffs), phase_diffs[:,0])
# Read actual phases from antenna hdf5
actual_antenna_measured_phases = { a.name: 2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas }
##############################
# Compare actual time shifts #
##############################
actual_antenna_true_phases = { a.name: -2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas }
# Compare actual time shifts
my_phase_diffs = []
for i,b in enumerate(baselines):
actual_phase_measured_diff = lib.phase_mod( lib.phase_mod(actual_antenna_measured_phases[b[1].name]) - lib.phase_mod(actual_antenna_measured_phases[b[0].name]))
actual_true_phase_diff = lib.phase_mod( lib.phase_mod(actual_antenna_true_phases[b[1].name]) - lib.phase_mod(actual_antenna_true_phases[b[0].name]))
# remove phase due to time delay from transmitter difference
tds = np.array([ lib.geometry_time(tx, ant, c_light=c_light) for ant in base])
delta_td = tds[1] - tds[0]
delta_td_phase = lib.phase_mod(delta_td*2*np.pi*f_beacon)
actual_phase_diff = lib.phase_mod( actual_phase_measured_diff - delta_td_phase)
my_phase_diffs.append(actual_phase_diff)
this_actual_true_phase_diff = lib.phase_mod( actual_true_phase_diff )
my_phase_diffs.append(this_actual_true_phase_diff)
# Make a plot
if True:
N_base = len(baselines)
N_ant = len(antennas)
phase_residuals = lib.phase_mod(phase_diffs[:,1] - my_phase_diffs)
for i in range(2):
plot_residuals = i == 1
colors = ['blue', 'orange']
fig, axs = plt.subplots(2, 1, sharex=True)
axs[0].set_title("Measured phase difference - Actual phase difference")
axs[0].set_xlabel("Phase $\\Delta\\varphi = \\varphi_{meas} - \\varphi_{true}$")
axs[0].tick_params(bottom=True, labelbottom=True)
#axs[1].tick_params(top=True, labeltop=True, bottom=False, labelbottom=False)
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 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]')
i=0
axs[i].set_ylabel("#")
axs[i].hist(phase_residuals, bins='sqrt', density=False)
if plot_residuals:
phase_residuals = lib.phase_mod(phase_diffs[:,1] - my_phase_diffs)
i=1
axs[i].set_ylabel("Baseline\n combination #")
if not True:
axs[i].plot(my_phase_diffs, np.arange(N_base), ls='none', marker='+', label='actual time shifts')
l = axs[i].plot(phase_diffs[:,1], np.arange(N_base), ls='none', marker='x', label='calculated')
axs[0].set_title("Difference between Measured and Actual phase difference for Baseline (i,j)")
axs[1].set_xlabel("Baseline Phase Residual $\\varphi_{ij_{meas}} - \\varphi_{ij_{true}}$ [rad]")
else:
axs[0].set_title("Comparison Measured and Actual phase difference for Baseline (i,j)")
axs[1].set_xlabel("Baseline Phase $\\varphi_{ij}$ [rad]")
axs[i].legend()
i=0
axs[i].set_ylabel("#")
if plot_residuals:
axs[i].hist(phase_residuals, bins='sqrt', density=False, alpha=0.8, color=colors[0])
else:
axs[i].hist(phase_diffs[:,1], bins='sqrt', density=False, alpha=0.8, color=colors[0], ls='solid' , histtype='step', label='Measured')
axs[i].hist(my_phase_diffs, bins='sqrt', density=False, alpha=0.8, color=colors[1], ls='dashed', histtype='step', label='Actual')
else:
axs[i].plot(phase_residuals, np.arange(N_base), ls='none', marker='x')
fig.tight_layout()
i=1
axs[i].set_ylabel("Baseline no.")
if not plot_residuals:
axs[i].plot(phase_diffs[:,1], np.arange(N_base), alpha=0.8, color=colors[0], ls='none', marker='x', label='calculated')
axs[i].plot(my_phase_diffs, np.arange(N_base), alpha=0.8, color=colors[1], ls='none', marker='+', label='actual time shifts')
if fig_dir:
fig.savefig(path.join(fig_dir, __file__ + f".F{freq_name}.pdf"))
axs[i].legend()
else:
axs[i].plot(phase_residuals, np.arange(N_base), alpha=0.6, ls='none', marker='x', color=colors[0])
fig.tight_layout()
if fig_dir:
extra_name = "measured"
if plot_residuals:
extra_name = "differences"
fig.savefig(path.join(fig_dir, __file__ + f".{extra_name}.F{freq_name}.pdf"))
if show_plots:
plt.show()