diff --git a/simulations/airshower_beacon_simulation/ca_beacon_periods.py b/simulations/airshower_beacon_simulation/ca_beacon_periods.py deleted file mode 100755 index 307b1cc..0000000 --- a/simulations/airshower_beacon_simulation/ca_beacon_periods.py +++ /dev/null @@ -1,176 +0,0 @@ -#!/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 - - fname = "ZH_airshower/mysim.sry" - - show_plots = True - save_result = 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 = antennas_fname - - # Read in antennas from file - f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) - - # run over all baselines - if ref_ant_id is None: - print("Doing all baselines") - baselines = list(combinations(antennas,2)) - # use ref_ant - else: - ref_ant = antennas[ref_ant_id] - print(f"Doing all baselines with {ref_ant.name}") - baselines = list(zip_longest([], antennas, fillvalue=ref_ant)) - - freq_names = antennas[0].beacon_info.keys() - if len(freq_names) > 1: - raise NotImplementedError - - freq_name = next(iter(freq_names)) - - # Determine integer multiple of periods to shift - # and True phase differences - time_diffs = np.empty( (len(baselines), 3) ) - for i, base in enumerate(baselines): - if i%50==0: - print(i, "out of", len(baselines)) - - # which traces to keep track of - traces = [ base[0].E_AxB, base[1].E_AxB ] - - # read f_beacon from the first antenna - f_beacon = base[0].beacon_info[freq_name]['freq'] - - # how many samples do we need to shift - sample_shifts, maxima = lib.coherence_sum_maxima(traces[0], traces[1], periodic=False) - best_sample_shift = sample_shifts[np.argmax(maxima)] - - # turn sample_shift into time - sampling_dt = (base[1].t[1] - base[1].t[0]) # ns - delta_t_coherence = sampling_dt*best_sample_shift # ns - - # get the amount of periods to move - k_period, t_rest = np.divmod(delta_t_coherence, 1/f_beacon) - - # always keep the reference before traces[1] - if t_rest < 0: # np.divmod already does this - k_period -= 1 - t_rest = 1/f_beacon + t_rest - - # 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[0]) - lib.phase_mod(true_phases[1])) - except IndexError: - # true_phase not determined yet - print(f"Missing true_phases for {freq_name} in baseline {base[0].name},{base[1].name}") - true_phases_diff = np.nan - - # save k_period with antenna names - time_diffs[i] = [true_phases_diff, k_period, f_beacon] - - # Plotting for one or two iterations - if show_plots and (i in [ 1, 57 ] or k_period > 3): - # More than three periods is quite much so report it - - # Show correlation maxima plot - if not True: - fig, ax = plt.subplots() - ax.set_title(f"Correlation Maxima {i}") - ax.set_xlabel("k") - ax.set_ylabel("Maximum correlation") - ax.plot(ks, maxima) - ax.plot(best_k, maxima[max_idx], marker='X') - - # Delta t due to the beacon - # Note that we want to show some overlapping waveforms - # Therefore we use the phase from the original waveforms - # and not the true_phases (we lost t_d information there) - phases = np.array([ant.beacon_info[freq_name]['phase'] for ant in base]) - phases_diff = lib.phase_mod(lib.phase_mod(phases[0]) - lib.phase_mod(phases[1])) - delta_t_beacon = phases_diff/(2*np.pi*f_beacon) - - - # Do we make it a shared plot with both the - # signal and the beacon? - beacons = None - if True: - try: - beacons = [ base[0].beacon, base[1].beacon ] - except: - # No beacon waveforms available.. - pass - - # Start the figure - fig, axs = plt.subplots(1+(beacons is not None), 1, sharex=True) - if beacons is None: - axs = [axs] - - print('i',i,f"B({base[0].name},{base[1].name})", 'k[T]',k_period, 'rest[ns]',t_rest, 'T[ns]',1/f_beacon, 'dT_coher[ns]', delta_t_coherence, 'dT_beac[ns]', delta_t_beacon) - - fig.suptitle( - ", ".join([ - f"$\\Delta$t_beacon [ns]: {delta_t_beacon:.2f}", - f"$\\Delta\\varphi$: {phases_diff:.4f}", - f"$\\Delta\\sigma_\\varphi$: {true_phases_diff:.4f}", - f"", - ]) - ) - - - axs[-1].set_xlabel('t [ns]') - axs[0].set_ylabel('Amplitude [a.u.]') - - # plot vertical lines indicating f_beacon - min_t, max_t = min(base[0].t[0], base[1].t[0]), max(base[0].t[-1], base[1].t[-1]) - N_lines = int( (max_t - min_t)*f_beacon) +1 - for i, t in enumerate(np.arange(N_lines)/f_beacon): - for ax in axs: - ax.axvline( min_t + t, color='k', alpha=0.3) - - # Plot traces - l1 = axs[0].plot(base[0].t, traces[0], label=f'Ref: {base[0].name}', alpha=0.8) - l2 = axs[0].plot(base[1].t, traces[1], label=f'Orig: {base[1].name}', alpha=0.3, marker='+', ms=5) - axs[0].plot(base[0].t + k_period/f_beacon + t_rest, traces[1], label='Coherence', alpha=0.3, marker='x', ms=5) - l3 = axs[0].plot(base[1].t - delta_t_beacon +k_period/f_beacon, traces[1], label=f'$\\Delta t_{{\\varphi}}$ + ($k={k_period:.0f}$)T', alpha=0.8) - - axs[0].legend(fancybox=True, framealpha=0.5) - - # Plot beacon if available - if beacons is not None: - ax = axs[1] - ax.set_title("Original Beacons") - ax.plot(base[0].t, beacons[0], label=f'Ref: {base[0].name}', alpha=0.8, color=l1[0].get_color()) - ax.plot(base[1].t, beacons[1], label=f'Orig: {base[1].name}', alpha=0.3, marker='+', ms=5, color=l2[0].get_color()) - ax.plot(base[1].t -delta_t_beacon +k_period/f_beacon, beacons[1], label=f'$\\Delta t_{{\\varphi}}$ + ($k={k_period:.0f}$)T', alpha=0.8, color=l3[0].get_color()) - - if False: - if False: - fig.savefig(__file__ + f"_i{i}_k{k_period}_zoomed_out.pdf") - ax.set_xlim(base[0].t[0]-1/f_beacon, base[0].t[0] + 5/f_beacon) - - fig.savefig(__file__ + f"_i{i}_k{k_period}.pdf") - - # Report back to CLI - print("Period Multiples resolved from", antennas_fname) - if save_result: - # Save integer periods to antennas - beacon.write_baseline_time_diffs_hdf5(time_diffs_fname, baselines, time_diffs[:,0], time_diffs[:,1], time_diffs[:,2]) - print("Timediffs saved to", time_diffs_fname) - - plt.show()