From bb776d358c720542c3019004ae5b4b40c80eaf8e Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Mon, 14 Nov 2022 20:49:35 +0100 Subject: [PATCH] WUotD --- .../aa_generate_beacon.py | 46 ++-- .../ab_modify_clocks.py | 8 +- .../ba_beacon_phases.py | 197 ++++++++++++++++++ .../show_beacon_amplitude_antennas.py | 65 ++++++ 4 files changed, 302 insertions(+), 14 deletions(-) create mode 100755 simulations/airshower_beacon_simulation/ba_beacon_phases.py create mode 100755 simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py diff --git a/simulations/airshower_beacon_simulation/aa_generate_beacon.py b/simulations/airshower_beacon_simulation/aa_generate_beacon.py index f8939e2..433bad8 100755 --- a/simulations/airshower_beacon_simulation/aa_generate_beacon.py +++ b/simulations/airshower_beacon_simulation/aa_generate_beacon.py @@ -1,4 +1,10 @@ #!/usr/bin/env python3 +# vim: fdm=indent ts=4 + +""" +Add a beacon measurement on top of the +simulated airshower. +""" import numpy as np import json @@ -36,7 +42,7 @@ def read_tx_file(fname): return tx, f_beacon -def read_beacon_hdf5(fname): +def read_beacon_hdf5(fname, traces_key='traces'): with h5py.File(fname, 'r') as h5: tx_attrs = h5['tx'].attrs f_beacon = tx_attrs.get('f_beacon') @@ -48,12 +54,15 @@ def read_beacon_hdf5(fname): for k, ant in h5['antennas'].items(): mydict = { k:ant.attrs.get(k) for k in ['x', 'y', 'z', 'name'] } antenna = Antenna(**mydict) - antenna.t = ant['traces'][0] - antenna.Ex = ant['traces'][1] - antenna.Ey = ant['traces'][2] - antenna.Ez = ant['traces'][3] - if len(ant['traces']) > 4: - antenna.beacon = ant['traces'][4] + antenna.t = ant[traces_key][0] + antenna.Ex = ant[traces_key][1] + antenna.Ey = ant[traces_key][2] + antenna.Ez = ant[traces_key][3] + if len(ant[traces_key]) > 4: + antenna.beacon = ant[traces_key][4] + + if ant.attrs: + antenna.attrs = {**ant.attrs} antennas.append(antenna) @@ -118,13 +127,23 @@ def append_antenna_hdf5(fname, antenna, columns = [], name='traces', prepend_tim if __name__ == "__main__": from os import path - remake_tx = False + remake_tx = True fname = "ZH_airshower/mysim.sry" - tx = Antenna(x=-500,y=0,z=0,name='tx') - f_beacon = 50e-3 # GHz + tx = Antenna(x=-500,y=0,z=0,name='tx') # m + f_beacon = 51.53e-3 # GHz beacon_amplitudes = 1e-6*np.array([1e2, 0, 0]) # mu V/m + beacon_radiate_rsq = True + + if beacon_radiate_rsq: + # Move tx out, and magnify beacon_amplitude (at tx) + tx = Antenna(x=-20e3,y=0,z=0,name='tx') # m + dist = lib.distance(tx, Antenna(x=0, y=0, z=0)) + + ampl = dist**2 + + beacon_amplitudes *= ampl #### fname_dir = path.dirname(fname) @@ -145,14 +164,15 @@ if __name__ == "__main__": # make beacon per antenna for i, antenna in enumerate(ev.antennas): - beacon = lib.beacon_from(tx, antenna, f_beacon, antenna.t) + beacon = lib.beacon_from(tx, antenna, f_beacon, antenna.t, radiate_rsq=beacon_radiate_rsq) + E = np.array([antenna.Ex, antenna.Ey, antenna.Ez, beacon]) append_antenna_hdf5( antennas_fname, antenna, E, name='orig_traces', prepend_time=True) # add to relevant polarisation - for i, _ in enumerate(beacon_amplitudes): - E[i] += beacon_amplitudes[i]*beacon + for j, _ in enumerate(beacon_amplitudes): + E[j] += beacon_amplitudes[j]*beacon append_antenna_hdf5( antennas_fname, antenna, E, name='traces', prepend_time=True) diff --git a/simulations/airshower_beacon_simulation/ab_modify_clocks.py b/simulations/airshower_beacon_simulation/ab_modify_clocks.py index e5418b5..578758e 100755 --- a/simulations/airshower_beacon_simulation/ab_modify_clocks.py +++ b/simulations/airshower_beacon_simulation/ab_modify_clocks.py @@ -1,4 +1,10 @@ #!/usr/bin/env python3 +# vim: fdm=indent ts=4 + +""" +Add a uniformly sampled time offset +to the clock of each antenna. +""" import numpy as np import json @@ -21,7 +27,7 @@ if __name__ == "__main__": import sys max_clock_offset = 100# ns - remake_clock_offsets = False + remake_clock_offsets = True seed = 12345 rng = np.random.default_rng(seed) diff --git a/simulations/airshower_beacon_simulation/ba_beacon_phases.py b/simulations/airshower_beacon_simulation/ba_beacon_phases.py new file mode 100755 index 0000000..5780946 --- /dev/null +++ b/simulations/airshower_beacon_simulation/ba_beacon_phases.py @@ -0,0 +1,197 @@ +#!/usr/bin/env python3 +# vim: fdm=indent ts=4 + +import numpy as np +import h5py + +import lib + +import aa_generate_beacon as beacon + +from lib import direct_fourier_transform + +from numpy.polynomial import Polynomial + +def find_beacon_in_traces( + traces, + t_trace, + f_beacon_estimate = 50e6, + frequency_fit = False, + N_test_freqs = 5e2, + f_beacon_estimate_band = 0.01, + amp_cut = 0.8 + ): + """ + f_beacon_band is inclusive + + traces is [trace, trace, trace, .. ] + """ + + amplitudes = np.zeros(len(traces)) + phases = np.zeros(len(traces)) + frequencies = np.zeros(len(traces)) + + if frequency_fit: # fit frequency + test_freqs = f_beacon_estimate + f_beacon_estimate_band * np.linspace(-1, 1, int(N_test_freqs)+1) + ft_amp_gen = direct_fourier_transform(test_freqs, t_trace, (x for x in traces)) + + n_samples = len(t_trace) + + for i, ft_amp in enumerate(ft_amp_gen): + real, imag = ft_amp + amps = 1/n_samples * ( real**2 + imag**2)**0.5 + + # find frequency peak and surrounding + # bins valid for parabola fitting + max_amp_idx = np.argmax(amps) + max_amp = amps[max_amp_idx] + + if True: + frequencies[i] = test_freqs[max_amp_idx] + continue + + valid_mask = amps >= amp_cut*max_amp + + if True: # make sure not to use other peaks + lower_mask = valid_mask[0:max_amp_idx] + upper_mask = valid_mask[max_amp_idx:] + + if any(lower_mask): + lower_end = np.argmin(lower_mask[::-1]) + else: + lower_end = max_amp_idx + + if any(upper_mask): + upper_end = np.argmin(upper_mask) + else: + upper_end = 0 + + + valid_mask[0:(max_amp_idx - lower_end)] = False + valid_mask[(max_amp_idx + upper_end):] = False + + if all(~valid_mask): + frequencies[i] = np.nan + continue + + # fit Parabola + parafit = Polynomial.fit(test_freqs[valid_mask], amps[valid_mask], 2) + func = parafit.convert() + + # find frequency where derivative is 0 + deriv = func.deriv(1) + freq = deriv.roots()[0] + + frequencies[i] = freq + else: + frequencies[:] = f_beacon_estimate + + # evaluate fourier transform at freq for each trace + for i, freq in enumerate(frequencies): + if freq is np.nan: + phases[i] = np.nan + amplitudes[i] = np.nan + continue + + real, imag = direct_fourier_transform(freq, t_trace, traces[i]) + + phases[i] = np.arctan2(real, imag) + amplitudes[i] = 1/len(t_trace) * (real**2 + imag**2)**0.5 + + return frequencies, phases, amplitudes + +if __name__ == "__main__": + from os import path + import sys + + + f_beacon_band = (49e-3,55e-3) #GHz + allow_frequency_fitting = True + read_frequency_from_file = True + + fname = "ZH_airshower/mysim.sry" + + #### + fname_dir = path.dirname(fname) + antennas_fname = path.join(fname_dir, beacon.antennas_fname) + + if not path.isfile(antennas_fname): + print("Antenna file cannot be found, did you try generating a beacon?") + sys.exit(1) + + # read in antennas + with h5py.File(antennas_fname, 'a') as fp: + if 'antennas' not in fp.keys(): + print("Antenna file corrupted? no antennas") + sys.exit(1) + + group = fp['antennas'] + + f_beacon = None + if read_frequency_from_file and 'tx' in fp: + tx = fp['tx'] + if 'f_beacon' in tx.attrs: + f_beacon = tx.attrs['f_beacon'] + else: + print("No frequency found in file.") + sys.exit(2) + f_beacon_estimate_band = 0.01*f_beacon + + elif allow_frequency_fitting: + f_beacon_estimate_band = (f_beacon_band[1] - f_beacon_band[0])/2 + f_beacon = f_beacon_band[1] - f_beacon_estimate_band + else: + print("Not allowed to fit frequency and no tx group found in file.") + sys.exit(2) + + N_antennas = len(group.keys()) + # just for funzies + found_data = np.zeros((N_antennas, 3)) + + # Determine frequency and phase + for i, name in enumerate(group.keys()): + ant_group = group[name] + + if 'traces' not in ant_group.keys(): + print(f"Antenna file corrupted? no 'traces' in {name}") + sys.exit(1) + + traces = ant_group['traces'] + + freqs, phases, amps = find_beacon_in_traces( + traces[1:-1], traces[0], + f_beacon_estimate=f_beacon, + frequency_fit=allow_frequency_fitting, + f_beacon_estimate_band=f_beacon_estimate_band + ) + + # only take Ex for now + frequency = freqs[-1] + phase = phases[-1] + amplitude = amps[-1] + + print(frequency, phase, amplitude) + + ant_group.attrs['beacon_freq'] = frequency + ant_group.attrs['beacon_phase'] = phase + ant_group.attrs['beacon_amplitude'] = amplitude + ant_group.attrs['beacon_orientation'] = 'Ex' + + found_data[i] = frequency, phase, amplitude + + # show histogram of found frequencies + if True: + import matplotlib.pyplot as plt + if True or allow_frequency_fitting: + fig, ax = plt.subplots() + ax.set_xlabel("Frequency") + ax.set_ylabel("Counts") + ax.hist(found_data[:,0], bins='auto', density=False) + + if True: + fig, ax = plt.subplots() + ax.set_xlabel("Amplitudes") + ax.set_ylabel("Counts") + ax.hist(found_data[:,2], bins='auto', density=False) + + plt.show() diff --git a/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py b/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py new file mode 100755 index 0000000..a682c97 --- /dev/null +++ b/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python3 +# vim: fdm=indent ts=4 + +__doc__ = \ +""" +Show the beacon amplitude per antenna. +""" + +import numpy as np +import h5py +import matplotlib.pyplot as plt + +import aa_generate_beacon as beacon +import lib + +if __name__ == "__main__": + import os.path as path + + fname = "ZH_airshower/mysim.sry" + + #### + fname_dir = path.dirname(fname) + antennas_fname = path.join(fname_dir, beacon.antennas_fname) + + f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) + + beacon_frequencies = np.array([ant.attrs['beacon_freq'] for ant in antennas]) + beacon_amplitudes = np.array([ant.attrs['beacon_amplitude'] for ant in antennas]) + beacon_phases = np.array([ant.attrs['beacon_phase'] for ant in antennas]) + + ##### + sizes = 64 + if True: + vals = beacon_phases + colorlabel = '$\\varphi$' + sizes = 64*(beacon_amplitudes/np.max(beacon_amplitudes))**2 + else: + vals = beacon_amplitudes + colorlabel = "[$\\mu$V/m]" + + x = [ a.x for a in antennas ] + y = [ a.y for a in antennas ] + ##### + + fig, axs = plt.subplots() + axs.set_title("Amplitude at beacon frequency at each antenna") + axs.set_aspect('equal', 'datalim') + axs.set_xlabel('[m]') + axs.set_ylabel('[m]') + + if True: + # underlie a calculate phase field + xs = np.linspace( np.min(x), np.max(x), 50) + ys = np.linspace( np.min(y), np.max(y), 50) + phases, (xs, ys) = lib.phase_field_from_tx(xs, ys, tx, f_beacon, return_meshgrid=False) + sc2 = axs.scatter(xs, ys, c=phases, alpha=0.5, zorder=-5) + fig.colorbar(sc2, ax=axs) + + sc = axs.scatter(x, y, c=vals, s=sizes) + + axs.plot(tx.x, tx.y, marker='X', color='k') + + fig.colorbar(sc, ax=axs, label=colorlabel) + + plt.show()