m-thesis-introduction/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py

#!/usr/bin/env python3
# vim: fdm=indent ts=4

import h5py
from itertools import combinations, zip_longest
import matplotlib.pyplot as plt
import numpy as np

import aa_generate_beacon as beacon
import lib


if __name__ == "__main__":
    from os import path
    import sys

    import os
    import matplotlib
    if os.name == 'posix' and "DISPLAY" not in os.environ:
        matplotlib.use('Agg')

    fname = "ZH_airshower/mysim.sry"

    show_plots = True
    ref_ant_id = None # leave None for all baselines

    ####
    fname_dir = path.dirname(fname)
    antennas_fname = path.join(fname_dir, beacon.antennas_fname)
    time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname
    fig_dir = "./figures" # set None to disable saving

    basenames, time_diffs, f_beacons, true_phase_diffs, k_periods = beacon.read_baseline_time_diffs_hdf5(time_diffs_fname)

    f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)

    # TODO: allow multiple frequencies
    if (f_beacon != f_beacons).any():
        raise NotImplementedError

    N_base = len(basenames)
    N_ant = len(antennas)

    # reshape time_diffs into N_ant x N_ant array
    sigma_phase_matrix = np.full( (N_ant, N_ant), np.nan, dtype=float)
    #sigma_phase_matrix = np.arange(0, N_ant**2, dtype=float).reshape( (N_ant,N_ant))

    name2idx = lambda name: int(name)-1
    for i, b in enumerate(basenames):
        idx = (name2idx(b[0]), name2idx(b[1]))

        if idx[0] == idx[1]:
            pass

        sigma_phase_matrix[(idx[0], idx[1])] = true_phase_diffs[i]
        sigma_phase_matrix[(idx[1], idx[0])] = true_phase_diffs[i]

    # for each row j subtract the 0,j element from the whole row
    # and apply phase_mod
    first_row = sigma_phase_matrix[0]

    sigma_phase_matrix = sigma_phase_matrix - first_row[:,np.newaxis]
    sigma_phase_matrix = lib.phase_mod(sigma_phase_matrix)


    # Except for the first row, these are all separate measurements
    # Condense into phase offset per antenna
    if True: # do not use the first row
        my_mask = np.arange(1, len(sigma_phase_matrix), dtype=int)
    else:
        my_mask = np.arange(0, len(sigma_phase_matrix), dtype=int)

    mean_sigma_phase = np.nanmean(sigma_phase_matrix[my_mask], axis=0)
    std_sigma_phase  = np.nanstd( sigma_phase_matrix[my_mask], axis=0)

    # write into antenna hdf5
    with h5py.File(antennas_fname, 'a') as fp:
        group = fp['antennas']

        freq_name = None
        for i, ant in enumerate(antennas):
            h5ant = group[ant.name]

            h5beacon_info = h5ant['beacon_info']

            # find out freq_name
            if freq_name is None:
                freq_name = [ k for k in h5beacon_info.keys() if np.isclose(h5beacon_info[k].attrs['freq'], f_beacon)][0]

            h5attrs = h5beacon_info[freq_name].attrs

            idx = name2idx(ant.name)
            h5attrs['sigma_phase_mean'] = mean_sigma_phase[idx]
            h5attrs['sigma_phase_std']  =  std_sigma_phase[idx]


    ##############################
    # Compare actual time shifts #
    ##############################
    antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in antennas }

    if True:
        # show means and std
        fig, axs = plt.subplots(1,2, sharey=True)
        axs[0].set_title("")
        axs[0].set_xlabel("Antenna no.")
        axs[0].set_ylabel("Antenna Phase $\\Delta_\\varphi$")
        axs[1].set_xlabel("#")
        if True:
            forward = lambda x: x/(2*np.pi*f_beacon)
            inverse = lambda x: 2*np.pi*x*f_beacon
            secax = axs[1].secondary_yaxis('right', functions=(forward, inverse))
            secax.set_ylabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]')

        l = axs[0].errorbar(np.arange(N_ant), mean_sigma_phase, yerr=std_sigma_phase, marker='.', alpha=0.7)

        axs[1].hist(mean_sigma_phase, bins='sqrt', density=False, orientation='horizontal', color=l[0].get_color(), histtype='step')

        # Actual time shifts
        if True:
            actual_phase_shifts = [ -1*lib.phase_mod(2*np.pi*f_beacon*v) for k,v in antenna_time_shifts.items() ]
            antenna_names = [int(k) for k,v in antenna_time_shifts.items() ]

            # make sure to keep the same offset
            if True:
                phase_offset = mean_sigma_phase[0] - actual_phase_shifts[0]

                actual_phase_shifts += phase_offset

            l = axs[0].plot(antenna_names, actual_phase_shifts, ls='none', marker='3', alpha=0.8)

            if True:
                axs[1].hist(actual_phase_shifts, bins='sqrt', density=False, orientation='horizontal', ls='dashed', color=l[0].get_color(), histtype='step')

        fig.tight_layout()
        if fig_dir:
            fig.savefig(path.join(fig_dir, __file__ + f".residuals.pdf"))

    ##########################
    ##########################
    ##########################

    actual_time_shifts = []
    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_time_shifts.append(actual_time_shift)

    # unpack mean_sigma_phase back into a list of time diffs
    measured_time_diffs = []
    for i,b in enumerate(basenames):
        time0, time1 = mean_sigma_phase[name2idx(b[0])], mean_sigma_phase[name2idx(b[1])]
        measured_time_diffs.append(time1 - time0)

    # Make a plot
    if True:
        fig, ax = plt.subplots()
        ax.set_xlabel("Baseline no.")
        ax.set_ylabel("$\\Delta t$[ns]")
        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(np.arange(N_base), actual_time_shifts, marker='+', label='actual time shifts')
        ax.plot(np.arange(N_base), measured_time_diffs, marker='x', label='calculated')

        ax.legend()
        if fig_dir:
            fig.savefig(path.join(fig_dir, __file__ + f".calculated_shifts.pdf"))

    if show_plots:
        plt.show()