ericteunis/m-thesis-introduction
1
0
Fork 0
mirror of https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git synced 2024-11-13 10:03:32 +01:00
Code Issues Projects Releases Packages Wiki Activity

ZH: renaming phase variables II: true_phase->clock_phase

Browse source
This commit is contained in:
Eric Teunis de Boone 2023-01-19 16:24:05 +01:00
parent 3ba45f4f52
commit 97ebdea54f
6 changed files with 41 additions and 41 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

14
simulations/airshower_beacon_simulation/aa_generate_beacon.py
View file

@ -208,13 +208,13 @@ def read_baseline_time_diffs_hdf5(fname):
dset = group[dset_name] dset = group[dset_name]
time_diffs = dset[:,0] time_diffs = dset[:,0]
f_beacon = dset[:,1] f_beacon = dset[:,1]
true_phase_diffs = dset[:,2] clock_phase_diffs = dset[:,2]
k_periods = dset[:,3] k_periods = dset[:,3]
return names, time_diffs, f_beacon, true_phase_diffs, k_periods return names, time_diffs, f_beacon, clock_phase_diffs, k_periods
def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods, f_beacon, time_diffs=None, overwrite=True): def write_baseline_time_diffs_hdf5(fname, baselines, clock_phase_diffs, k_periods, f_beacon, time_diffs=None, overwrite=True):
""" """
Write a combination of baselines, phase_diff, k_period and f_beacon to file. Write a combination of baselines, phase_diff, k_period and f_beacon to file.
@ -225,7 +225,7 @@ def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods
# this is a single baseline # this is a single baseline
N_baselines = 1 N_baselines = 1
baselines = [baselines] baselines = [baselines]
true_phase_diffs = [true_phase_diffs] clock_phase_diffs = [clock_phase_diffs]
k_periods = [k_periods] k_periods = [k_periods]
f_beacon = np.array([f_beacon]) f_beacon = np.array([f_beacon])
@ -237,9 +237,9 @@ def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods
f_beacon = np.array([f_beacon]*N_baselines) f_beacon = np.array([f_beacon]*N_baselines)
if time_diffs is None: if time_diffs is None:
time_diffs = k_periods/f_beacon + true_phase_diffs/(2*np.pi*f_beacon) time_diffs = k_periods/f_beacon + clock_phase_diffs/(2*np.pi*f_beacon)
assert len(baselines) == len(true_phase_diffs) == len(k_periods) == len(f_beacon) assert len(baselines) == len(clock_phase_diffs) == len(k_periods) == len(f_beacon)
with h5py.File(fname, 'a') as fp: with h5py.File(fname, 'a') as fp:
group_name = 'baseline_time_diffs' group_name = 'baseline_time_diffs'
@ -263,7 +263,7 @@ def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods
base_dset = group.create_dataset(base_dset_name, data=basenames) base_dset = group.create_dataset(base_dset_name, data=basenames)
data = np.vstack( (time_diffs, f_beacon, true_phase_diffs, k_periods) ).T data = np.vstack( (time_diffs, f_beacon, clock_phase_diffs, k_periods) ).T
dset = group.create_dataset(dset_name, data=data) dset = group.create_dataset(dset_name, data=data)
# }}} vim marker # }}} vim marker

30
simulations/airshower_beacon_simulation/bb_measure_true_phase.py
View file

@ -57,19 +57,19 @@ if __name__ == "__main__":
f_beacon = antennas[0].beacon_info[freq_name]['freq'] f_beacon = antennas[0].beacon_info[freq_name]['freq']
true_phases = lib.remove_antenna_geometry_phase(tx, antennas, f_beacon, beacon_phases, c_light=c_light) clock_phases = lib.remove_antenna_geometry_phase(tx, antennas, f_beacon, beacon_phases, c_light=c_light)
# Remove the phase from one antenna # Remove the phase from one antenna
# this is a free parameter # this is a free parameter
# (only required for absolute timing) # (only required for absolute timing)
if remove_absolute_phase_offset_first_antenna or remove_absolute_phase_offset_minimum: if remove_absolute_phase_offset_first_antenna or remove_absolute_phase_offset_minimum:
if remove_absolute_phase_offset_first_antenna: # just take the first phase if remove_absolute_phase_offset_first_antenna: # just take the first phase
minimum_phase = true_phases[0] minimum_phase = clock_phases[0]
else: # take the minimum else: # take the minimum
minimum_phase = np.min(true_phases, axis=-1) minimum_phase = np.min(clock_phases, axis=-1)
true_phases -= minimum_phase clock_phases -= minimum_phase
true_phases = lib.phase_mod(true_phases) clock_phases = lib.phase_mod(clock_phases)
# Save to antennas in file # Save to antennas in file
with h5py.File(antennas_fname, 'a') as fp: with h5py.File(antennas_fname, 'a') as fp:
@ -79,7 +79,7 @@ if __name__ == "__main__":
h5ant = fp['antennas'][ant.name] h5ant = fp['antennas'][ant.name]
h5beacon_freq = h5ant['beacon_info'][freq_name] h5beacon_freq = h5ant['beacon_info'][freq_name]
h5beacon_freq.attrs['true_phase'] = true_phases[i] h5beacon_freq.attrs['clock_phase'] = clock_phases[i]
# Plot True Phases at their locations # Plot True Phases at their locations
if show_plots or fig_dir: if show_plots or fig_dir:
@ -96,7 +96,7 @@ if __name__ == "__main__":
#scatter_kwargs['vmax'] = +np.pi #scatter_kwargs['vmax'] = +np.pi
color_label='$\\varphi(\\sigma_t)$ [rad]' color_label='$\\varphi(\\sigma_t)$ [rad]'
sc = ax.scatter(*antenna_locs, c=true_phases, **scatter_kwargs) sc = ax.scatter(*antenna_locs, c=clock_phases, **scatter_kwargs)
fig.colorbar(sc, ax=ax, label=color_label) fig.colorbar(sc, ax=ax, label=color_label)
if False: if False:
@ -114,20 +114,20 @@ if __name__ == "__main__":
fig, ax = plt.subplots() fig, ax = plt.subplots()
fig.suptitle('Clock phase Residuals\nf_beacon={:2.0f}MHz'.format(f_beacon*1e3)) fig.suptitle('Clock phase Residuals\nf_beacon={:2.0f}MHz'.format(f_beacon*1e3))
actual_true_phases = np.array([ -2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas ]) actual_clock_phases = np.array([ -2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas ])
# Modify actual_true_phases, the same way as true_phases # Modify actual_clock_phases, the same way as clock_phases
# was modified # was modified
if remove_absolute_phase_offset_first_antenna or remove_absolute_phase_offset_minimum: if remove_absolute_phase_offset_first_antenna or remove_absolute_phase_offset_minimum:
if remove_absolute_phase_offset_first_antenna: # just take the first phase if remove_absolute_phase_offset_first_antenna: # just take the first phase
minimum_phase = actual_true_phases[0] minimum_phase = actual_clock_phases[0]
else: # take the minimum else: # take the minimum
minimum_phase = np.min(actual_true_phases, axis=-1) minimum_phase = np.min(actual_clock_phases, axis=-1)
actual_true_phases -= minimum_phase actual_clock_phases -= minimum_phase
actual_true_phases = lib.phase_mod(actual_true_phases) actual_clock_phases = lib.phase_mod(actual_clock_phases)
true_phase_residuals = lib.phase_mod(true_phases - actual_true_phases) clock_phase_residuals = lib.phase_mod(clock_phases - actual_clock_phases)
antenna_locs = list(zip(*[(ant.x, ant.y) for ant in antennas])) antenna_locs = list(zip(*[(ant.x, ant.y) for ant in antennas]))
ax.set_xlabel('x' if spatial_unit is None else 'x [{}]'.format(spatial_unit)) ax.set_xlabel('x' if spatial_unit is None else 'x [{}]'.format(spatial_unit))
@ -136,7 +136,7 @@ if __name__ == "__main__":
scatter_kwargs['cmap'] = 'inferno' scatter_kwargs['cmap'] = 'inferno'
color_label='$\\Delta\\varphi(\\sigma_t) = \\varphi_{meas} - \\varphi_{true}$ [rad]' color_label='$\\Delta\\varphi(\\sigma_t) = \\varphi_{meas} - \\varphi_{true}$ [rad]'
sc = ax.scatter(*antenna_locs, c=true_phase_residuals, **scatter_kwargs) sc = ax.scatter(*antenna_locs, c=clock_phase_residuals, **scatter_kwargs)
fig.colorbar(sc, ax=ax, label=color_label) fig.colorbar(sc, ax=ax, label=color_label)
if fig_dir: if fig_dir:

20
simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py
View file

@ -70,30 +70,30 @@ if __name__ == "__main__":
# Get true phase diffs # Get true phase diffs
try: try:
true_phases = np.array([ant.beacon_info[freq_name]['true_phase'] for ant in base]) clock_phases = np.array([ant.beacon_info[freq_name]['clock_phase'] for ant in base])
true_phases_diff = lib.phase_mod(lib.phase_mod(true_phases[1]) - lib.phase_mod(true_phases[0])) clock_phases_diff = lib.phase_mod(lib.phase_mod(clock_phases[1]) - lib.phase_mod(clock_phases[0]))
except IndexError: except IndexError:
# true_phase not determined yet # clock_phase not determined yet
print(f"Missing true_phases for {freq_name} in baseline {base[0].name},{base[1].name}") print(f"Missing clock_phases for {freq_name} in baseline {base[0].name},{base[1].name}")
true_phases_diff = np.nan clock_phases_diff = np.nan
# save phase difference with antenna names # save phase difference with antenna names
phase_diffs[i] = [f_beacon, true_phases_diff] phase_diffs[i] = [f_beacon, clock_phases_diff]
beacon.write_baseline_time_diffs_hdf5(time_diffs_fname, baselines, phase_diffs[:,1], [0]*len(phase_diffs), phase_diffs[:,0]) beacon.write_baseline_time_diffs_hdf5(time_diffs_fname, baselines, phase_diffs[:,1], [0]*len(phase_diffs), phase_diffs[:,0])
############################## ##############################
# Compare actual time shifts # # Compare actual time shifts #
############################## ##############################
actual_antenna_true_phases = { a.name: -2*np.pi*a.attrs['clock_offset']*f_beacon for a in sorted(antennas, key=lambda a: int(a.name)) } actual_antenna_clock_phases = { a.name: -2*np.pi*a.attrs['clock_offset']*f_beacon for a in sorted(antennas, key=lambda a: int(a.name)) }
# Compare actual time shifts # Compare actual time shifts
my_phase_diffs = [] my_phase_diffs = []
for i,b in enumerate(baselines): for i,b in enumerate(baselines):
actual_true_phase_diff = lib.phase_mod( lib.phase_mod(actual_antenna_true_phases[b[1].name]) - lib.phase_mod(actual_antenna_true_phases[b[0].name])) actual_clock_phase_diff = lib.phase_mod( lib.phase_mod(actual_antenna_clock_phases[b[1].name]) - lib.phase_mod(actual_antenna_clock_phases[b[0].name]))
this_actual_true_phase_diff = lib.phase_mod( actual_true_phase_diff ) this_actual_clock_phase_diff = lib.phase_mod( actual_clock_phase_diff )
my_phase_diffs.append(this_actual_true_phase_diff) my_phase_diffs.append(this_actual_clock_phase_diff)
# Make a plot # Make a plot
if True: if True:

6
simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py
View file

@ -36,7 +36,7 @@ if __name__ == "__main__":
time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname
fig_dir = args.fig_dir # set None to disable saving fig_dir = args.fig_dir # set None to disable saving
basenames, time_diffs, f_beacons, true_phase_diffs, k_periods = beacon.read_baseline_time_diffs_hdf5(time_diffs_fname) basenames, time_diffs, f_beacons, clock_phase_diffs, k_periods = beacon.read_baseline_time_diffs_hdf5(time_diffs_fname)
f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
@ -59,8 +59,8 @@ if __name__ == "__main__":
for i, b in enumerate(basenames): for i, b in enumerate(basenames):
idx = (name2idx(b[0]), name2idx(b[1])) idx = (name2idx(b[0]), name2idx(b[1]))
sigma_phase_matrix[(idx[0], idx[1])] = lib.phase_mod(true_phase_diffs[i]) sigma_phase_matrix[(idx[0], idx[1])] = lib.phase_mod(clock_phase_diffs[i])
sigma_phase_matrix[(idx[1], idx[0])] = lib.phase_mod(-1*true_phase_diffs[i]) sigma_phase_matrix[(idx[1], idx[0])] = lib.phase_mod(-1*clock_phase_diffs[i])
mat_kwargs = dict( mat_kwargs = dict(
norm = Normalize(vmin=-np.pi, vmax=+np.pi), norm = Normalize(vmin=-np.pi, vmax=+np.pi),

6
simulations/airshower_beacon_simulation/lib/lib.py
View file

@ -70,16 +70,16 @@ def remove_antenna_geometry_phase(tx, antennas, f_beacon, measured_phases=None,
if not hasattr(measured_phases, '__len__'): if not hasattr(measured_phases, '__len__'):
measured_phases = [measured_phases] measured_phases = [measured_phases]
true_phases = np.empty( (len(antennas)) ) clock_phases = np.empty( (len(antennas)) )
for i, ant in enumerate(antennas): for i, ant in enumerate(antennas):
measured_phase = measured_phases[i] measured_phase = measured_phases[i]
geom_time = geometry_time(tx, ant, c_light=c_light) geom_time = geometry_time(tx, ant, c_light=c_light)
geom_phase = geom_time * 2*np.pi*f_beacon geom_phase = geom_time * 2*np.pi*f_beacon
true_phases[i] = phase_mod(measured_phase - geom_phase) clock_phases[i] = phase_mod(measured_phase - geom_phase)
return true_phases return clock_phases
""" Fourier """ """ Fourier """

6
simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py
View file

@ -34,8 +34,8 @@ if __name__ == "__main__":
beacon_amplitudes = np.array([ant.beacon_info[freq_name]['amplitude'] for ant in antennas]) beacon_amplitudes = np.array([ant.beacon_info[freq_name]['amplitude'] for ant in antennas])
beacon_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['beacon_phase'] for ant in antennas])) beacon_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['beacon_phase'] for ant in antennas]))
if False and 'true_phase' in antennas[0].beacon_info[freq_name]: if False and 'clock_phase' in antennas[0].beacon_info[freq_name]:
beacon_true_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['true_phase'] for ant in antennas])) beacon_clock_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['clock_phase'] for ant in antennas]))
else: else:
subtitle = " Phases from t0" subtitle = " Phases from t0"
beacon_frequencies = np.array([ f_beacon for ant in antennas ]) beacon_frequencies = np.array([ f_beacon for ant in antennas ])
@ -49,7 +49,7 @@ if __name__ == "__main__":
colorlabel = '$\\varphi$' colorlabel = '$\\varphi$'
sizes = 64*(beacon_amplitudes/np.max(beacon_amplitudes))**2 sizes = 64*(beacon_amplitudes/np.max(beacon_amplitudes))**2
elif True: # True Phases elif True: # True Phases
vals = beacon_true_phases vals = beacon_clock_phases
colorlabel = '$\\sigma_\\varphi$' colorlabel = '$\\sigma_\\varphi$'
plot_phase_field = False plot_phase_field = False
plot_tx = False plot_tx = False