ericteunis/m-thesis-introduction
1
0
Fork 0
mirror of https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git synced 2025-01-02 00:13:30 +01:00
Code Issues Projects Releases Packages Wiki Activity

ZH: employ phase_comparison_figure function

Browse source
This commit is contained in:
Eric Teunis de Boone 2023-02-07 17:58:45 +01:00
parent e821ce00ca
commit 8fc7d4bc0c
4 changed files with 71 additions and 79 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

47
simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py
View file

@ -8,6 +8,7 @@ import numpy as np
import aa_generate_beacon as beacon import aa_generate_beacon as beacon
import lib import lib
from lib import figlib
if __name__ == "__main__": if __name__ == "__main__":
@ -108,45 +109,41 @@ if __name__ == "__main__":
for i in range(2): for i in range(2):
plot_residuals = i == 1 plot_residuals = i == 1
colors = ['blue', 'orange']
fig, axs = plt.subplots(2, 1, sharex=True, figsize=figsize) true_phases = my_phase_diffs
measured_phases = phase_diffs[:,1]
if True: hist_kwargs = {}
phase2time = lambda x: x/(2*np.pi*f_beacon) if plot_residuals:
time2phase = lambda x: 2*np.pi*x*f_beacon measured_phases = lib.phase_mod(measured_phases - true_phases)
secax = axs[0].secondary_xaxis('top', functions=(phase2time, time2phase)) hist_kwargs['histtype'] = 'stepfilled'
secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]')
fig = figlib.phase_comparison_figure(
measured_phases,
true_phases,
plot_residuals=plot_residuals,
f_beacon=f_beacon,
figsize=figsize,
hist_kwargs=hist_kwargs,
)
axs = fig.get_axes()
if plot_residuals: if plot_residuals:
phase_residuals = lib.phase_mod(phase_diffs[:,1] - my_phase_diffs)
fig.suptitle("Difference between Measured and Actual phase difference\n for Baselines (i,j" + (')' if not ref_ant_id else '='+str([ int(a.name) for a in ref_ants])+')')) fig.suptitle("Difference between Measured and Actual phase difference\n for Baselines (i,j" + (')' if not ref_ant_id else '='+str([ int(a.name) for a in ref_ants])+')'))
axs[-1].set_xlabel("Baseline Phase Residual $\\Delta\\varphi_{ij_{meas}} - \\Delta\\varphi_{ij_{true}}$ [rad]") axs[-1].set_xlabel("Baseline Phase Residual $\\Delta\\varphi_{ij_{meas}} - \\Delta\\varphi_{ij_{true}}$ [rad]")
else: else:
fig.suptitle("Comparison Measured and Actual phase difference\n for Baselines (i,j" + (')' if not ref_ant_id else '='+str([ int(a.name) for a in ref_ants])+')')) fig.suptitle("Comparison Measured and Actual phase difference\n for Baselines (i,j" + (')' if not ref_ant_id else '='+str([ int(a.name) for a in ref_ants])+')'))
axs[-1].set_xlabel("Baseline Phase $\\Delta\\varphi_{ij}$ [rad]") axs[-1].set_xlabel("Baseline Phase $\\Delta\\varphi_{ij}$ [rad]")
#
i=0 i=0
axs[i].set_ylabel("#") secax = axs[i].child_axes[0]
if plot_residuals: secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]')
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')
#
i=1 i=1
axs[i].set_ylabel("Baseline no.") axs[i].set_ylabel("Baseline no.")
if plot_residuals:
axs[i].plot(phase_residuals, np.arange(N_base), alpha=0.6, ls='none', marker='x', color=colors[0])
else:
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')
axs[i].legend()
fig.tight_layout()
if fig_dir: if fig_dir:
extra_name = "measured" extra_name = "measured"

42
simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py
View file

@ -10,6 +10,7 @@ import numpy as np
import aa_generate_beacon as beacon import aa_generate_beacon as beacon
import lib import lib
from lib import figlib
if __name__ == "__main__": if __name__ == "__main__":
@ -176,43 +177,36 @@ if __name__ == "__main__":
for i in range(2): for i in range(2):
plot_residuals = i == 1 plot_residuals = i == 1
colors = ['blue', 'orange'] true_phases = actual_antenna_phase_shifts
measured_phases = mean_clock_phase
fig, axs = plt.subplots(2, 1, sharex=True, figsize=figsize) hist_kwargs = {}
if plot_residuals:
measured_phases = lib.phase_mod(measured_phases - actual_antenna_phase_shifts)
hist_kwargs['histtype'] = 'stepfilled'
if True: fig = figlib.phase_comparison_figure(
phase2time = lambda x: x/(2*np.pi*f_beacon) measured_phases,
time2phase = lambda x: 2*np.pi*x*f_beacon true_phases,
secax = axs[0].secondary_xaxis('top', functions=(phase2time, time2phase)) plot_residuals=plot_residuals,
secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') f_beacon=f_beacon,
figsize=figsize,
hist_kwargs=hist_kwargs,
)
axs = fig.get_axes()
if plot_residuals: if plot_residuals:
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 Mean Phase Residual $\\Delta_\\varphi$") axs[-1].set_xlabel("Antenna Mean Phase Residual $\\Delta_\\varphi$")
else: else:
fig.suptitle("Comparison Measured and Actual phases (minus global phase)\n for Antenna $i$") fig.suptitle("Comparison Measured and Actual phases (minus global phase)\n for Antenna $i$")
axs[-1].set_xlabel("Antenna Mean Phase $\\varphi$") axs[-1].set_xlabel("Antenna Mean Phase $\\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_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')
i=1 i=1
axs[i].set_ylabel("Antenna no.") axs[i].set_ylabel("Antenna no.")
if plot_residuals: #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(phase_residuals, np.arange(N_ant), alpha=0.6, ls='none', marker='x', color=colors[0])
else:
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].legend()
fig.tight_layout() fig.tight_layout()
if fig_dir: if fig_dir:

47
simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py
View file

@ -10,6 +10,7 @@ import numpy as np
from os import path from os import path
import aa_generate_beacon as beacon import aa_generate_beacon as beacon
from lib import figlib
if __name__ == "__main__": if __name__ == "__main__":
import sys import sys
@ -76,15 +77,27 @@ if __name__ == "__main__":
for i in range(2): for i in range(2):
plot_residuals = i == 1 plot_residuals = i == 1
colors = ['blue', 'orange']
fig, axs = plt.subplots(2, 1, sharex=True, figsize=figsize)
if True: true_phases = actual_time_shifts
phase2time = lambda x: x/(2*np.pi*f_beacon) measured_phases = measured_time_shifts
time2phase = lambda x: 2*np.pi*x*f_beacon
secax = axs[0].secondary_xaxis('top', functions=(time2phase, phase2time)) hist_kwargs = {}
secax.set_xlabel('Phase $2\\pi t f_{beac}$ [rad]') if plot_residuals:
measured_phases = measured_phases - true_phases
hist_kwargs['histtype'] = 'stepfilled'
fig = figlib.phase_comparison_figure(
measured_phases,
true_phases,
plot_residuals=plot_residuals,
f_beacon=f_beacon,
figsize=figsize,
hist_kwargs=hist_kwargs,
secondary_axis='phase',
)
axs = fig.get_axes()
if plot_residuals: if plot_residuals:
time_shift_residuals = measured_time_shifts - actual_time_shifts time_shift_residuals = measured_time_shifts - actual_time_shifts
@ -94,26 +107,6 @@ if __name__ == "__main__":
fig.suptitle("Comparison Measured and Actual clock offset") fig.suptitle("Comparison Measured and Actual clock offset")
axs[-1].set_xlabel("Antenna Time Offset $t_c = \\left(\\frac{\\Delta\\varphi}{2\\pi} + k\\right) / f_{beac}$ [ns]") axs[-1].set_xlabel("Antenna Time Offset $t_c = \\left(\\frac{\\Delta\\varphi}{2\\pi} + k\\right) / f_{beac}$ [ns]")
i=0
axs[i].set_ylabel("#")
if plot_residuals:
axs[i].hist(time_shift_residuals, bins='sqrt', alpha=0.8, color=colors[0])
else:
axs[i].hist(measured_time_shifts, bins='sqrt', density=False, alpha=0.8, color=colors[0], ls='solid' , histtype='step', label='Measured')
axs[i].hist(actual_time_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(time_shift_residuals, np.arange(N_ant), alpha=0.6, ls='none', marker='x', color=colors[0])
else:
axs[i].errorbar(measured_time_shifts, np.arange(N_ant), yerr=None, marker='4', alpha=0.7, ls='none', color=colors[0], label='Measured')
axs[i].plot(actual_time_shifts, antenna_names, ls='none', marker='3', alpha=0.8, color=colors[1], label='Actual')
axs[i].legend()
fig.tight_layout()
if fig_dir: if fig_dir:
extra_name = "comparison" extra_name = "comparison"
if plot_residuals: if plot_residuals:

14
simulations/airshower_beacon_simulation/lib/figlib.py
View file

@ -10,6 +10,7 @@ def phase_comparison_figure(
sc_kwargs={}, sc_kwargs={},
colors=['blue', 'orange'], colors=['blue', 'orange'],
legend_on_scatter=True, legend_on_scatter=True,
secondary_axis='time',
**fig_kwargs **fig_kwargs
): ):
""" """
@ -27,11 +28,18 @@ def phase_comparison_figure(
if not hasattr(axs, '__len__'): if not hasattr(axs, '__len__'):
axs = [axs] axs = [axs]
if f_beacon: if f_beacon and secondary_axis in ['phase', 'time']:
phase2time = lambda x: x/(2*np.pi*f_beacon) phase2time = lambda x: x/(2*np.pi*f_beacon)
time2phase = lambda x: 2*np.pi*x*f_beacon time2phase = lambda x: 2*np.pi*x*f_beacon
secax = axs[0].secondary_xaxis('top', functions=(phase2time, time2phase))
secax.set_xlabel('Time $\\varphi/(2\\pi f_{beac})$ [ns]') if secondary_axis == 'time':
functions = (phase2time, time2phase)
label = 'Time $\\varphi/(2\\pi f_{beac})$ [ns]'
else:
functions = (time2phase, phase2time)
label = 'Phase $2\\pi t f_{beac}$ [rad]'
secax = axs[0].secondary_xaxis('top', functions=functions)
# Histogram # Histogram
if do_hist_plot: if do_hist_plot: