mirror of
https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git
synced 2024-12-22 19:43:30 +01:00
94 lines
3.3 KiB
Python
94 lines
3.3 KiB
Python
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
|
|
from scipy.stats import norm
|
|
|
|
def phase_comparison_figure(
|
|
measured_phases,
|
|
true_phases,
|
|
plot_residuals=True,
|
|
f_beacon=None,
|
|
hist_kwargs={},
|
|
sc_kwargs={},
|
|
colors=['blue', 'orange'],
|
|
legend_on_scatter=True,
|
|
secondary_axis='time',
|
|
fit_gaussian=False,
|
|
**fig_kwargs
|
|
):
|
|
"""
|
|
Create a figure comparing measured_phase against true_phase
|
|
by both plotting the values, and the residuals.
|
|
"""
|
|
default_fig_kwargs = dict(sharex=True)
|
|
fig_kwargs = {**default_fig_kwargs, **fig_kwargs}
|
|
|
|
do_hist_plot = hist_kwargs is not False
|
|
do_scatter_plot = sc_kwargs is not False
|
|
|
|
fig, axs = plt.subplots(0+do_hist_plot+do_scatter_plot, 1, **fig_kwargs)
|
|
|
|
if not hasattr(axs, '__len__'):
|
|
axs = [axs]
|
|
|
|
if f_beacon and secondary_axis in ['phase', 'time']:
|
|
phase2time = lambda x: x/(2*np.pi*f_beacon)
|
|
time2phase = lambda x: 2*np.pi*x*f_beacon
|
|
|
|
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
|
|
if do_hist_plot:
|
|
i=0
|
|
default_hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step')
|
|
hist_kwargs = {**default_hist_kwargs, **hist_kwargs}
|
|
|
|
axs[i].set_ylabel("#")
|
|
|
|
_counts, _bins, _patches = axs[i].hist(measured_phases, color=colors[0], label='Measured', ls='solid', **hist_kwargs)
|
|
if not plot_residuals: # also plot the true clock phases
|
|
axs[i].hist(true_phases, color=colors[1], label='Actual', ls='dashed', **{**hist_kwargs, **dict(bins=_bins)})
|
|
|
|
if fit_gaussian:# Fit a gaussian to the histogram
|
|
gauss_kwargs = dict(color=colors[0], ls='dotted')
|
|
xs = np.linspace(min(_bins), max(_bins))
|
|
|
|
mean, std = norm.fit(measured_phases)
|
|
dx = _bins[1] - _bins[0]
|
|
scale = len(measured_phases)*dx
|
|
fit_ys = norm.pdf(xs, mean, std) * scale
|
|
|
|
axs[i].axvline(mean, ls='dotted', color='k', lw=2)
|
|
|
|
axs[i].plot( xs, fit_ys, label='fit', **gauss_kwargs)
|
|
|
|
# put mean, sigma and chisq on plot
|
|
text_str = f"$\\mu$ = {mean: .2e}\n$\\sigma$ = {std: .2e}"
|
|
text_kwargs = dict( transform=axs[i].transAxes, fontsize=14, verticalalignment='top')
|
|
axs[i].text(0.02, 0.95, text_str, **text_kwargs)
|
|
|
|
# Scatter plot
|
|
if do_scatter_plot:
|
|
i=1
|
|
default_sc_kwargs = dict(alpha=0.6, ls='none')
|
|
sc_kwargs = {**default_sc_kwargs, **sc_kwargs}
|
|
|
|
axs[i].set_ylabel("Antenna no.")
|
|
axs[i].plot(measured_phases, np.arange(len(measured_phases)), marker='x' if plot_residuals else '3', color=colors[0], label='Measured', **sc_kwargs)
|
|
|
|
if not plot_residuals: # also plot the true clock phases
|
|
axs[i].plot(true_phases, np.arange(len(true_phases)), marker='4', color=colors[1], label='Actual', **sc_kwargs)
|
|
|
|
if not plot_residuals and legend_on_scatter:
|
|
axs[i].legend()
|
|
|
|
fig.tight_layout()
|
|
|
|
return fig
|