From 36d1e8136cb92b14ce52c1fc5b9ed71cc5eda88b Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Wed, 12 Oct 2022 15:46:50 +0200 Subject: [PATCH] Extend figure to show phase gradient also --- simulations/10_timegradient_per_ref_ant.py | 87 +++++++++++++++++++--- 1 file changed, 75 insertions(+), 12 deletions(-) diff --git a/simulations/10_timegradient_per_ref_ant.py b/simulations/10_timegradient_per_ref_ant.py index be2a1a5..a54471e 100755 --- a/simulations/10_timegradient_per_ref_ant.py +++ b/simulations/10_timegradient_per_ref_ant.py @@ -50,6 +50,13 @@ def geometry_time(dist, x2=None, c_light=c_light): return dist/c_light +def phase_mod(phase, low=np.pi): + """ + Modulo phase such that it falls within the + interval $[-low, 2\pi - low)$. + """ + return (phase + low) % (2*np.pi) - low + def antenna_triangles(antennas): return combinations(antennas, 3) @@ -75,7 +82,7 @@ def random_antenna(N_ant=1, antenna_ranges=[10e3,10e3,10e3], max_clock_skew=1): return antennas -def single_baseline_referenced_sigmas(tx, baseline, all_antennas): +def single_baseline_referenced_sigmas(tx, baseline, all_antennas, phase_func=None): N_ant = len(all_antennas) baseline_ts = np.array([b.t for b in baseline]) @@ -87,11 +94,14 @@ def single_baseline_referenced_sigmas(tx, baseline, all_antennas): for j, ant in enumerate(filter(not_baseline, all_antennas)): t_diff = ant.t - baseline_ts geo_diff = geometry_time(tx, ant) - baseline_geo - sigmas[j] = t_diff - geo_diff + if phase_func is not None: + sigmas[i] = phase_func(t_diff - geo_diff) + else: + sigmas[i] = t_diff - geo_diff return sigmas -def reference_antenna_sigmas(tx, ref_ant, all_antennas): +def reference_antenna_sigmas(tx, ref_ant, all_antennas, phase_func=None): N_ant = len(all_antennas) ref_geo = geometry_time(tx, ref_ant) @@ -103,28 +113,33 @@ def reference_antenna_sigmas(tx, ref_ant, all_antennas): t_diff = ant.t - ref_ant.t geo_diff = geometry_time(tx, ant) - ref_geo - sigmas[i] = t_diff - geo_diff + if phase_func is not None: + sigmas[i] = phase_func(t_diff - geo_diff) + else: + sigmas[i] = t_diff - geo_diff return sigmas -def all_sigmas_using_reference_antenna(tx, all_antennas): +def all_sigmas_using_reference_antenna(tx, all_antennas, phase_func=None): N_ant = len(all_antennas) sigmas = np.empty( (N_ant,N_ant) ) for i, ant in enumerate(all_antennas): - sigmas[i] = reference_antenna_sigmas(tx, ant, all_antennas) + sigmas[i] = reference_antenna_sigmas(tx, ant, all_antennas, phase_func=phase_func) return sigmas -def main(tx, antennas, spatial_unit=None, time_unit=None, ref_idx = [0, 1, -2, -1]): +def main(tx, antennas, spatial_unit=None, time_unit=None, ref_idx = [0, 1, -2, -1], plot_phase=False, remove_minimum=True, f_beacon=50e6, scatter_kwargs={}): # Use each baseline once as a reference # and loop over the remaining antennas N_ant = len(antennas) fig = None - baseline = [antennas[0], antennas[1]] + default_scatter_kwargs = {} + #for i, baseline in enumerate(antenna_baselines(antennas)): if False: + baseline = [antennas[0], antennas[1]] sigmas = single_baseline_referenced_sigmas(tx, baseline, antennas) print("Baseline {},{}".format(baseline[0].name, baseline[1].name)) print(sigmas) @@ -133,9 +148,47 @@ def main(tx, antennas, spatial_unit=None, time_unit=None, ref_idx = [0, 1, -2, - print() if True: - sigmas = all_sigmas_using_reference_antenna(tx, antennas) + if plot_phase: + phase_func = lambda t: phase_mod(2*np.pi* f_beacon * t) + color_label='$\\varphi$' + default_scatter_kwargs['cmap'] = 'Spectral_r' + default_scatter_kwargs['vmin'] = -np.pi + default_scatter_kwargs['vmax'] = +np.pi + else: + color_label='t' if time_unit is None else 't ['+time_unit+']' + phase_func = None + + scatter_kwargs = { **default_scatter_kwargs, **scatter_kwargs } + + sigmas = all_sigmas_using_reference_antenna(tx, antennas, phase_func=phase_func) + + if remove_minimum: + if True: + # actually use the time diffs with the first ref ant + # required for phase alignment + mins = sigmas[0] + else: + mins = -1*np.min(sigmas, axis=-1) + + sigmas = sigmas + mins[:, np.newaxis] + + if plot_phase: + # Redo the phase mod + sigmas = phase_mod(sigmas) fig, axs = plt.subplots(2,2, sharex=True, sharey=True) + title = "" + if remove_minimum: + title += '$\sigma_{0j}$ added' + if remove_minimum and plot_phase: + title += ', ' + if plot_phase: + t_scaler = 1 + if time_unit == 'ns': + t_scaler = 1e9 + title += 'f= {:2.0f}MHz'.format(f_beacon*t_scaler/1e6) + + fig.suptitle(title) antenna_locs = list(zip(*[(ant.x, ant.y) for ant in antennas])) for i, ax in enumerate(axs.flat): @@ -143,8 +196,8 @@ def main(tx, antennas, spatial_unit=None, time_unit=None, ref_idx = [0, 1, -2, - ax.set_xlabel('x' if spatial_unit is None else 'x [{}]'.format(spatial_unit)) ax.set_ylabel('y' if spatial_unit is None else 'y [{}]'.format(spatial_unit)) - sc = ax.scatter(*antenna_locs, c=sigmas[ref_idx[i]]) - fig.colorbar(sc, ax=ax, label='t' if time_unit is None else 't ['+time_unit+']') + sc = ax.scatter(*antenna_locs, c=sigmas[ref_idx[i]], **scatter_kwargs) + fig.colorbar(sc, ax=ax, label=color_label) ax.plot(antennas[ref_idx[i]].x, antennas[ref_idx[i]].y, 'rx') @@ -159,8 +212,17 @@ if __name__ == "__main__": parser = ArgumentParser(description=__doc__) parser.add_argument("fname", metavar="path/to/figure[/]", nargs="?", help="Location for generated figure, will append __file__ if a directory. If not supplied, figure is shown.") parser.add_argument('num_ant', help='Number of antennas to use', nargs='?', default=5, type=int) + parser.add_argument('--remove-min', action='store_true') + + command_group = parser.add_mutually_exclusive_group(required=False) + command_group.add_argument('--time', help='Calculate times (Default)', action='store_true') + command_group.add_argument('--phase', help='Calculate wrapped phases', action='store_true') args = parser.parse_args() + args.rm_minimum = True + + args.plot_phase = args.phase + del args.time, args.phase if args.fname == 'none': args.fname = None @@ -174,12 +236,13 @@ if __name__ == "__main__": ###### antenna_ranges = np.array([10*km,10*km,5*km]) antenna_max_clock_skew = 100*ns/ns # 0.1 us + f_beacon = 50e6*ns # 50 MHz tx = Antenna(name='tx', x=-300*km, y=200*km, z=0) antennas = random_antenna(args.num_ant, antenna_ranges, antenna_max_clock_skew) add_spatial_time_delay(tx, antennas) - fig, sigmas = main(tx, antennas, spatial_unit='m', time_unit='ns') + fig, sigmas = main(tx, antennas, spatial_unit='m', time_unit='ns', plot_phase=args.plot_phase, remove_minimum=args.rm_minimum, f_beacon=f_beacon) ###### Output if args.fname is not None: