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

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

"""
Show how the Power changes when incorporating the
various clock offsets by plotting on a grid.
"""

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # required for projection='3d' on old matplotliblib versions
import numpy as np
from os import path
import joblib

from earsim import REvent
from atmocal import AtmoCal
import aa_generate_beacon as beacon
import lib
from lib import rit


if __name__ == "__main__":
    valid_cases = ['no_offset', 'repair_none', 'repair_phases', 'repair_all']

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

    atm = AtmoCal()

    from scriptlib import MyArgumentParser
    parser = MyArgumentParser()

    group = parser.add_argument_group('figures')
    for case in valid_cases:
        group.add_argument('--'+case.replace('_','-'), dest='figures', action='append_const', const=case)

    args = parser.parse_args()

    wanted_cases = args.figures

    if not wanted_cases or 'all' in wanted_cases:
        wanted_cases = valid_cases

    fname = "ZH_airshower/mysim.sry"

    fig_dir = args.fig_dir
    show_plots = args.show_plots

    remove_beacon_from_traces = True

    ####
    fname_dir = path.dirname(fname)
    antennas_fname = path.join(fname_dir, beacon.antennas_fname)
    pickle_fname = path.join(fname_dir, 'res.pkl')

    # create fig_dir
    if fig_dir:
        os.makedirs(fig_dir, exist_ok=True)

    # Read in antennas from file
    _, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
    # Read original REvent
    ev = REvent(fname)
    # .. patch in our antennas
    ev.antennas = antennas

    rit.set_pol_and_bp(ev)

    ##
    ## Setup grid
    ##
    X = 400
    zgr = 0 #not exact
    dXref = atm.distance_to_slant_depth(np.deg2rad(ev.zenith),750,zgr+ev.core[2])
    scale2d = dXref*np.tan(np.deg2rad(2.))
    scale4d = dXref*np.tan(np.deg2rad(4.))
    scale02d = dXref*np.tan(np.deg2rad(0.2))

    Nx, Ny = 21, 21
    scales = {
            'scale2d': scale2d,
            'scale4d': scale4d,
            'scale02d': scale02d,
            }

    # backup antenna times
    backup_antenna_t = [ ant.t for ant in ev.antennas ]
    backup_antenna_t_AxB = [ ant.t_AxB for ant in ev.antennas ]

    plot_titling = {
            'no_offset': "no clock offset",
            'repair_none': "unrepaired clock offset",
            'repair_phases': "phase resolved clock offsets repaired",
            'repair_all': "final measured clock offsets repaired"
            }

    # For now only implement using one freq_name
    freq_names = ev.antennas[0].beacon_info.keys()
    if len(freq_names) > 1:
        raise NotImplementedError

    freq_name = next(iter(freq_names))

    # Pre remove the beacon from the traces
    #
    # We need to remove it here, so we do not shoot ourselves in
    # the foot when changing to the various clock offsets.
    if remove_beacon_from_traces:
        for i, ant in enumerate(ev.antennas):
            beacon_phase = ant.beacon_info[freq_name]['beacon_phase']
            f = ant.beacon_info[freq_name]['freq']
            ampl_AxB = ant.beacon_info[freq_name]['amplitude']
            calc_beacon = lib.sine_beacon(f, ev.antennas[i].t, amplitude=ampl_AxB, phase=beacon_phase)

            # Only need to manipulate E_AxB
            ev.antennas[i].E_AxB -= calc_beacon

    with joblib.parallel_backend("loky"):
        for case in wanted_cases:
            print(f"Starting {case} figure")
            # Repair clock offsets with the measured offsets
            transl_modes = {'no_offset':'orig', 'repair_phases':'phases', 'repair_all':'all'}

            if case in transl_modes:
                transl_mode = transl_modes[case]
                measured_offsets = beacon.read_antenna_clock_repair_offsets(antennas, mode=transl_mode, freq_name=freq_name)
            else:
                measured_offsets = [0]*len(ev.antennas)

            for i, ant in enumerate(ev.antennas):
                total_clock_offset = measured_offsets[i]

                #ev.antennas[i].t =     backup_antenna_t[i] + total_clock_offset
                ev.antennas[i].t_AxB = backup_antenna_t_AxB[i] + total_clock_offset

            # Measure power on grid
            for scalename, scale in scales.items():
                wx, wy = scale, scale
                print(f"Starting grid measurement for figure {case} with {scalename}")
                xx, yy, p, maxp = rit.shower_plane_slice(ev, X=X, Nx=Nx, Ny=Nx, wx=wx, wy=wy)

                fig, axs = rit.slice_figure(ev, X, xx, yy, p, mode='sp')

                suptitle = fig._suptitle.get_text()
                fig.suptitle("")
                axs.set_title(suptitle +"\n" +plot_titling[case])
                #axs.set_aspect('equal', 'datalim')

                if fig_dir:
                    fig.tight_layout()
                    fig.savefig(path.join(fig_dir, path.basename(__file__) + f'.X{X}.{case}.{scalename}.pdf'))