mirror of
https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git
synced 2024-12-22 11:33:32 +01:00
ZH: remove period script assuming event comes from transmitter
This commit is contained in:
parent
86f493bd1f
commit
7851a13065
1 changed files with 0 additions and 176 deletions
|
|
@ -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()
|
||||
Loading…
Reference in a new issue