diff --git a/simulations/airshower_beacon_simulation/ba_beacon_phases.py b/simulations/airshower_beacon_simulation/ba_beacon_phases.py index 7919a0e..047f888 100755 --- a/simulations/airshower_beacon_simulation/ba_beacon_phases.py +++ b/simulations/airshower_beacon_simulation/ba_beacon_phases.py @@ -12,11 +12,13 @@ if __name__ == "__main__": from os import path import sys - + f_beacon_band = (49e-3,55e-3) #GHz - allow_frequency_fitting = True + allow_frequency_fitting = False read_frequency_from_file = True + show_plots = True + fname = "ZH_airshower/mysim.sry" #### @@ -65,40 +67,72 @@ if __name__ == "__main__": sys.exit(1) traces = ant_group['traces'] + t_trace = traces[0] - 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 - ) + if not True: + # only take the Beacon trace + test_traces = traces[4] + orients = ['B'] + else: + test_traces = traces[1:] + orients = ['Ex', 'Ey', 'Ez', 'B'] - # only take Ex for now - frequency = freqs[-1] - phase = phases[-1] - amplitude = amps[-1] + if True: + freqs, phases, amps = lib.find_beacon_in_traces( + test_traces, t_trace, + f_beacon_estimate=f_beacon, + frequency_fit=allow_frequency_fitting, + f_beacon_estimate_band=f_beacon_estimate_band + ) + else: + # Testing + freqs = [f_beacon] + t0 = ant_group.attrs['t0'] - print(frequency, phase, amplitude) + phases = [ 2*np.pi*t0*f_beacon ] + amps = [ 3e-7 ] + + # choose highest amp + #idx = np.argmax(amps, axis=1) + idx = 0 + + frequency = freqs[idx] + phase = phases[idx] + amplitude = amps[idx] + orientation = orients[idx] + + if show_plots and (i == 60 or i == 72): + fig, ax = plt.subplots() + trace_amp = max(traces[-1]) - min(traces[-1]) + + myt = np.linspace(min(traces[0]), max(traces[0]), 10*len(traces[0])) + ax.plot(t_trace, traces[-1], marker='.', label='trace') + ax.plot(myt, lib.sine_beacon(frequency, myt, amplitude=amplitude, phase=phase), ls='dashed', label='simulated') + ax.set_title(f"Beacon at antenna {ant_group}\nF:{frequency}, P:{phase}, A:{amplitude}") + ax.legend() ant_group.attrs['beacon_freq'] = frequency ant_group.attrs['beacon_phase'] = phase ant_group.attrs['beacon_amplitude'] = amplitude - ant_group.attrs['beacon_orientation'] = 'Ex' + ant_group.attrs['beacon_orientation'] = orientation found_data[i] = frequency, phase, amplitude + print("Beacon Phases, Amplitudes and Frequencies written to", antennas_fname) + # show histogram of found frequencies - if True: + if show_plots: 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) + ax.axvline(f_beacon, ls='dashed', color='g') + ax.hist(found_data[:,0], bins='sqrt', 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) + ax.hist(found_data[:,2], bins='sqrt', density=False) plt.show() diff --git a/simulations/airshower_beacon_simulation/lib.py b/simulations/airshower_beacon_simulation/lib.py index 42ac13a..7b14641 100644 --- a/simulations/airshower_beacon_simulation/lib.py +++ b/simulations/airshower_beacon_simulation/lib.py @@ -11,7 +11,7 @@ def sine_beacon(f, t, t0=0, amplitude=1, baseline=0, phase=0): """ Return a sine appropriate as a beacon """ - return amplitude * np.cos(2*np.pi*f*(t+t0) + phase) + baseline + return amplitude * np.sin(2*np.pi*f*(t+t0) + phase) + baseline def distance(x1, x2): @@ -74,6 +74,9 @@ def direct_fourier_transform(freqs, time, samplesets_iterable): return np.dot(c_k, samplesets_iterable), np.dot(s_k, samplesets_iterable) def phase_field_from_tx(x, y, tx, f_beacon, c_light=3e8, t0=0, wrap_phase=True, return_meshgrid=True): + + assert type(tx) in [Antenna] + xs, ys = np.meshgrid(x, y, sparse=True) x_distances = (tx.x - xs)**2 @@ -128,8 +131,8 @@ def find_beacon_in_traces( real, imag = ft_amp amps = 1/n_samples * ( real**2 + imag**2)**0.5 - # find frequency peak and surrounding - # bins valid for parabola fitting + # find frequency peak and surrounding bins + # that are valid for the parabola fit max_amp_idx = np.argmax(amps) max_amp = amps[max_amp_idx] @@ -170,9 +173,12 @@ def find_beacon_in_traces( freq = deriv.roots()[0] frequencies[i] = freq - else: + + else: # no frequency finding frequencies[:] = f_beacon_estimate + + n_samples = len(t_trace) # evaluate fourier transform at freq for each trace for i, freq in enumerate(frequencies): if freq is np.nan: @@ -183,6 +189,6 @@ def find_beacon_in_traces( real, imag = direct_fourier_transform(freq, t_trace, traces[i]) phases[i] = np.arctan2(real, imag) - amplitudes[i] = 2/len(t_trace) * (real**2 + imag**2)**0.5 + amplitudes[i] = 2/n_samples * (real**2 + imag**2)**0.5 return frequencies, phases, amplitudes