diff --git a/simulations/airshower_beacon_simulation/lib/lib.py b/simulations/airshower_beacon_simulation/lib/lib.py index 295916b..aa54b79 100644 --- a/simulations/airshower_beacon_simulation/lib/lib.py +++ b/simulations/airshower_beacon_simulation/lib/lib.py @@ -9,10 +9,16 @@ from earsim import Antenna c_light = 3e8*1e-9 # m/ns """ Beacon utils """ -def sine_beacon(f, t, t0=0, amplitude=1, baseline=0, phase=0): +def sine_beacon(f, t, t0=0, amplitude=1, baseline=0, phase=0, peak_at_t0=0): """ Return a sine appropriate as a beacon """ + baseline = baseline*np.ones_like(t) + + if peak_at_t0: # add peak near t0 + idx = (np.abs(t-t0)).argmin() + baseline[ max(0, idx-1):min(len(t), idx+1) ] += peak_at_t0 + amplitude + return amplitude * np.cos(2*np.pi*f*(t+t0) + phase) + baseline def phase_mod(phase, low=np.pi): diff --git a/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py b/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py index 4e86c34..b59fb71 100755 --- a/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py +++ b/simulations/airshower_beacon_simulation/lib/tests/test_beacon_fourier.py @@ -14,18 +14,20 @@ import lib seed = 12345 dt = 1 # ns -frequency = 45e-3 # GHz +frequency = 52.12345e-3 # GHz N = 5e2 +t_clock = 0 +beacon_amplitude = 1 t = np.arange(0, 10*int(1e3), dt, dtype=float) rng = np.random.default_rng(seed) -phase_res = np.zeros(int(N)) # Vary both the base time and the phase t_extra = 0 +phase_res = np.zeros(int(N)) +amp_res = np.zeros(int(N)) for i in range(int(N)): - # Change timebase t -= t_extra t_extra = (2*rng.uniform(size=1) - 1) *1e3 @@ -33,19 +35,74 @@ for i in range(int(N)): # Randomly phased beacon phase = lib.phase_mod(np.pi*(2*rng.uniform(size=1) -1)) # rad - beacon = lib.sine_beacon(frequency, t, t0=0, phase=phase) + beacon = beacon_amplitude * lib.sine_beacon(frequency, t, t0=0, phase=phase, peak_at_t0=False) - if True: # blank part of the beacon - blank_low, blank_high = 2*int(1e3), 4*int(1e3) - beacon[blank_low:blank_high] = 0 + if False: # blank part of the beacon + blank_low, blank_high = int(0.2*len(t)), int(0.4*len(t)) + t_mask = np.ones(len(t), dtype=bool) + t_mask[blank_low:blank_high] = False - measured = lib.find_beacon_in_traces([beacon], t, frequency, frequency_fit=False) + t = t[t_mask] + beacon = beacon[t_mask] - phase_res[i] = lib.phase_mod(measured[1][0] - phase) + # Introduce clock errors + t += t_clock + _, measured_phases, measured_amplitudes = lib.find_beacon_in_traces([beacon], t, frequency, frequency_fit=False) + t -= t_clock + + # Save residuals + phase_res[i] = lib.phase_mod(measured_phases[0] - phase) + amp_res[i] = measured_amplitudes[0] - beacon_amplitude + +### +## Present figures +### +phase2time = lambda x: x/(2*np.pi*frequency) +time2phase = lambda x: 2*np.pi*x*frequency fig, ax = plt.subplots() ax.set_title("Sine beacon phase determination\nwith random time shifts") ax.set_xlabel("$\\varphi_{meas} - \\varphi_{true}$ [rad]") ax.set_ylabel("#") ax.hist(phase_res, bins='sqrt') +if True: + ax.axvline( -1*lib.phase_mod(t_clock*frequency*2*np.pi), color='red', lw=5, alpha=0.8, label='true t_clock') +if True: + secax = ax.secondary_xaxis(1.0, functions=(phase2time, time2phase)) + secax.set_xlabel('Time [ns]') +ax.legend(title='N={:.1e}'.format(len(t))) +phase_xlims = ax.get_xlim() +fig.tight_layout() + +amp_xlims = None +if True: + fig, ax = plt.subplots() + ax.set_title("Amplitude") + ax.set_xlabel("$A_{meas} - 1$") + ax.set_ylabel("#") + ax.legend(title='N={:.1e}'.format(len(t))) + ax.hist(amp_res, bins='sqrt') + fig.tight_layout() + + amp_xlims = ax.get_xlim() + +if True: + fig, ax = plt.subplots() + ax.grid() + ax.set_xlabel("Phase Res [rad]") + ax.set_ylabel("Amplitude Res") + sc = ax.scatter( phase_res, amp_res ) + #fig.colorbar(sc, ax=ax) + + ax.set_xlim(phase_xlims) + if amp_xlims is not None: + ax.set_ylim(amp_xlims) + if True: + ax.axvline( -1*lib.phase_mod(t_clock*frequency*2*np.pi), color='red', lw=5, alpha=0.8, label='true t_clock') + if True: + secax = ax.secondary_xaxis(1.0, functions=(phase2time, time2phase)) + secax.set_xlabel('Time [ns]') + ax.legend(title='N={:.1e}'.format(len(t))) + fig.tight_layout() + plt.show() diff --git a/simulations/airshower_beacon_simulation/lib/tests/test_calculated_phase_measurement.py b/simulations/airshower_beacon_simulation/lib/tests/test_calculated_phase_measurement.py index 83cbb38..c1ec636 100755 --- a/simulations/airshower_beacon_simulation/lib/tests/test_calculated_phase_measurement.py +++ b/simulations/airshower_beacon_simulation/lib/tests/test_calculated_phase_measurement.py @@ -15,21 +15,22 @@ import lib seed = 12345 dt = 1 # ns -frequency = 45e-3 # GHz +frequency = 52.12345e-3 # GHz N = 5e2 +t_clock = 0 +beacon_amplitude = 1 c_light = lib.c_light -t_clock = 3/4 * 1/frequency t = np.arange(0, 10*int(1e3), dt, dtype=float) rng = np.random.default_rng(seed) -phase_in = lib.phase_mod(np.pi*(2*rng.uniform(size=1) -1)) # rad tx = Antenna(x=0,y=0,z=0,name='tx') -rx = Antenna(x=30,y=40,z=120,name='rx') +rx = Antenna(x=tx.x,y=tx.y,z=tx.z,name='rx') # Vary both the base time and the phase t_extra = 0 phase_res = np.zeros(int(N)) +amp_res = np.zeros(int(N)) for i in range(int(N)): # Change timebase t -= t_extra @@ -43,27 +44,76 @@ for i in range(int(N)): # Randomly phased beacon # at Antenna phase = lib.phase_mod(np.pi*(2*rng.uniform(size=1) -1)) # rad - beacon = lib.beacon_from(tx, rx, frequency, t, radiate_rsq=False, phase=phase, c_light=c_light) + beacon = beacon_amplitude * lib.beacon_from(tx, rx, frequency, t, t0=0, phase=phase, peak_at_t0=False, c_light=c_light, radiate_rsq=False) - if True: # blank part of the beacon - blank_low, blank_high = 2*int(1e3), 4*int(1e3) - beacon[blank_low:blank_high] = 0 + if False: # blank part of the beacon + blank_low, blank_high = int(0.2*len(t)), int(0.4*len(t)) + t_mask = np.ones(len(t), dtype=bool) + t_mask[blank_low:blank_high] = False + + t = t[t_mask] + beacon = beacon[t_mask] # Introduce clock errors t += t_clock - measured = lib.find_beacon_in_traces([beacon], t, frequency, frequency_fit=False) + _, measured_phases, measured_amplitudes = lib.find_beacon_in_traces([beacon], t, frequency, frequency_fit=False) t -= t_clock - calculated_phase = lib.remove_antenna_geometry_phase(tx, rx, frequency, measured[1][0], c_light=c_light) + calculated_phases = lib.remove_antenna_geometry_phase(tx, rx, frequency, measured_phases, c_light=c_light) - phase_res[i] = lib.phase_mod(calculated_phase - phase) + # Save residuals + phase_res[i] = lib.phase_mod(calculated_phases[0] - phase) + amp_res[i] = measured_amplitudes[0] - beacon_amplitude + +### +## Present figures +### +phase2time = lambda x: x/(2*np.pi*frequency) +time2phase = lambda x: 2*np.pi*x*frequency -# Make the figure fig, ax = plt.subplots() ax.set_title("Measured phase at Antenna - geometrical phase") ax.set_xlabel("$\\varphi_{meas} - \\varphi_{true}$ [rad]") ax.set_ylabel("#") -ax.hist(phase_res, bins='sqrt', density=False) -ax.axvline( -1*lib.phase_mod(t_clock*frequency*2*np.pi), color='red', lw=5, alpha=0.8, label='true t_clock') -ax.legend() +ax.hist(phase_res, bins='sqrt') +if True: + ax.axvline( -1*lib.phase_mod(t_clock*frequency*2*np.pi), color='red', lw=5, alpha=0.8, label='true t_clock') +if True: + secax = ax.secondary_xaxis(1.0, functions=(phase2time, time2phase)) + secax.set_xlabel('Time [ns]') +ax.legend(title='N={:.1e}'.format(len(t))) +phase_xlims = ax.get_xlim() +fig.tight_layout() + +amp_xlims = None +if True: + fig, ax = plt.subplots() + ax.set_title("Amplitude") + ax.set_xlabel("$A_{meas} - 1$") + ax.set_ylabel("#") + ax.legend(title='N={:.1e}'.format(len(t))) + ax.hist(amp_res, bins='sqrt') + fig.tight_layout() + + amp_xlims = ax.get_xlim() + +if True: + fig, ax = plt.subplots() + ax.grid() + ax.set_xlabel("Phase Res [rad]") + ax.set_ylabel("Amplitude Res") + sc = ax.scatter( phase_res, amp_res ) + #fig.colorbar(sc, ax=ax) + + ax.set_xlim(phase_xlims) + if amp_xlims is not None: + ax.set_ylim(amp_xlims) + if True: + ax.axvline( -1*lib.phase_mod(t_clock*frequency*2*np.pi), color='red', lw=5, alpha=0.8, label='true t_clock') + if True: + secax = ax.secondary_xaxis(1.0, functions=(phase2time, time2phase)) + secax.set_xlabel('Time [ns]') + ax.legend(title='N={:.1e}'.format(len(t))) + fig.tight_layout() + plt.show()