diff --git a/figures/beacon/.gitignore b/figures/beacon/.gitignore new file mode 100644 index 0000000..6633554 --- /dev/null +++ b/figures/beacon/.gitignore @@ -0,0 +1,2 @@ +*.png +*.pdf diff --git a/figures/beacon/Makefile b/figures/beacon/Makefile index 46b9688..0521137 100644 --- a/figures/beacon/Makefile +++ b/figures/beacon/Makefile @@ -3,11 +3,19 @@ all: dist dist: \ sine_beacon.pdf sine_beacon.png \ - ttl_beacon.pdf ttl_beacon.png + ttl_beacon.pdf ttl_beacon.png \ + field_singleleft_time.pdf field_singleleft_time.png \ + field_singleleft_phase.pdf field_singleleft_phase.png \ + field_singlecenter_time.pdf field_singlecenter_time.png \ + field_singlecenter_phase.pdf field_singlecenter_phase.png dist-clean: rm -v sine_beacon.* rm -v ttl_beacon.* + rm -v field_singleleft_time.* + rm -v field_singleleft_phase.* + rm -v field_singlecenter_time.* + rm -v field_singlecenter_phase.* beacon_spatial_time_difference_setup.pdf: src/beacon_spatial_time_difference_setup.py $< $@ @@ -18,3 +26,12 @@ sine_beacon.%: src/single_beacon.py ttl_beacon.%: src/single_beacon.py $< --periods 2 --with-rates ttl $@ +field_singleleft_time.%: src/beacon_field.py + $< --time $@ single-left +field_singleleft_phase.%: src/beacon_field.py + $< --phase $@ single-left + +field_singlecenter_time.%: src/beacon_field.py + $< --time $@ single-center +field_singlecenter_phase.%: src/beacon_field.py + $< --phase $@ single-center diff --git a/figures/beacon/src/beacon_field.py b/figures/beacon/src/beacon_field.py new file mode 100755 index 0000000..1ea5633 --- /dev/null +++ b/figures/beacon/src/beacon_field.py @@ -0,0 +1,273 @@ +#!/usr/bin/env python3 + +__doc__ = \ +""" +Generate a figure showing a value combining +the delays between a transmitter and a set of antennas +""" + + +import matplotlib.pyplot as plt +import numpy as np +from itertools import chain, combinations, product + +c_light = 3e8 +f_beacon = 50e6 + +class Antenna: + """ + Simple Antenna class + """ + def __init__(self,x=0,y=0,z=0,t0=0,name=""): + self.x = x + self.y = y + self.z = z + self.t = t0 + self.name = name + + def __repr__(self): + cls = self.__class__.__name__ + + return f'{cls}(x={self.x!r},y={self.y!r},z={self.z!r},t0={self.t!r},name={self.name!r})' + +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 dist(a,b): + return ((a.x-b.x)**2+(a.y-b.y)**2+(a.z-b.z)**2)**0.5 + +def phase(a,b,f=f_beacon,wrap=False): + t = dist(a,b)/c_light + phase = t*f*2*np.pi + + if wrap: + phase = phase_mod(phase) + + return phase + +def grid_plot(grid, ax=None, **plot_kwargs): + if ax is None: + ax = plt.gca() + + x = [a.x for a in grid] + y = [a.y for a in grid] + l = [a.name for a in grid] + ax.plot(x,y,'kx', **plot_kwargs) + for x_,y_,l_ in zip(x,y,l): + ax.annotate(l_,(x_,y_)) + +def antenna_combinations(ants, ref_ant=None): + if ref_ant is not None: # use only one reference antenna for the baselines + ref = ref_ant + ant_combi = [ (ref, j) for j in ants if j is not ref ] + else: # use all baselines + ant_combi = combinations(ants, 2) + + return list(ant_combi) + +def calculate_field(field, tx, ant_combi, ref_ant=None, calculate_phase=True): + if plot_phase: + p_ij = [ (phase(tx,i) - phase(tx,j)) for i,j in ant_combi ] + else: + t_ij = [ (dist(tx,i) - dist(tx, j))/c_light for i,j in ant_combi ] + + val = [] + xx = [] + yy = [] + + xs = field[0] + ys = field[1] + + for x in xs: + for y in ys: + tx_t = Antenna(x=x,y=y) + + if plot_phase: # phase + dp_ij = np.array([ (phase(tx_t,i) - phase(tx_t,j)) - p_ij[k] for k,(i,j) in enumerate(ant_combi) ]) + + if True: # phase wrap + dp_ij = phase_mod(dp_ij) + + val.append(np.sum(dp_ij**2)**0.5) + else: # time + dt_ij = np.array([ (dist(tx_t,i) - dist(tx_t, j))/c_light - t_ij[k] for k,(i,j) in enumerate(ant_combi) ]) + + val.append(np.sum(dt_ij**2)**0.5) + + xx.append(x) + yy.append(y) + + return np.array(xx), np.array(yy), np.array(val) + +def plot_field(tx, ants, xx, yy, val, ax=None, ref_ant=None, color_label='$\\left( t - \\tau \\right)^2$', plot_phase=None, mask=None,**scatter_kwargs): + if ax is None: + ax = plt.gca() + + default_scatter_kwargs = {} + default_scatter_kwargs['cmap'] = 'Spectral_r' + + if mask is not None: + val[mask] = np.nan + + if plot_phase is None: + pass + elif plot_phase: + default_scatter_kwargs['vmax'] = len(ants)*np.pi + #default_scatter_kwargs['cmap'] = 'gray' + pass + + scatter_kwargs = {**default_scatter_kwargs, **scatter_kwargs} + + grid_plot([tx] + ants, ax=ax) + if ref_ant is not None: + ax.set_title("Single reference antenna: {}, f: {}MHz".format(ref_ant.name, f_beacon/1e6)) + else: + ax.set_title("All baselines f: {} MHz".format(f_beacon/1e6)) + + sc = ax.scatter(xx,yy,c=val, **scatter_kwargs) + fig = ax.get_figure() + fig.colorbar(sc, ax=ax, label=color_label) + + return ax + +def square_grid(dx=1, N_x=10, dy=None, N_y=None, x_start=0, y_start=0): + N_y = N_x if N_y is None else N_y + dy = dx if dy is None else dy + + + print([(N_x,N_y), (dx,dy), (x_start,y_start)]) + + return product( + (x_start + n*dx for n in range(N_x)), + (y_start + n*dy for n in range(N_y)), + ) + +def triangular_grid(dx=1, N_x=10, dy=None, N_y=None, x_start=0, y_start=0): + + return chain( + square_grid(dx=dx,dy=dy,N_x=N_x,N_y=N_y,x_start=x_start,y_start=y_start), + square_grid(dx=dx,dy=dy,N_x=N_x,N_y=N_y,x_start=x_start + dx/2,y_start=y_start + dy/2), + ) + + + +if __name__ == "__main__": + from argparse import ArgumentParser + import os.path as path + + 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("type", choices=['single-left', 'single-center', 'square', 'tri', 'preset']) + + + command_group = parser.add_mutually_exclusive_group(required=True) + command_group.add_argument('--time', help='Use the time difference for the field', action='store_false') + command_group.add_argument('--phase', help='Use wrapped phase for the field', action='store_true') + + parser.add_argument('--ref', dest='ref_ant', metavar='ref_antenna', help='Number of antenna to use as reference') + + args = parser.parse_args() + + if args.fname is not None and path.isdir(args.fname): + args.fname = path.join(args.fname, path.splitext(path.basename(__file__))[0] + ".pdf") + + args.plot_phase = args.phase or args.time + del args.time, args.phase + + + if 'single' in args.type: # single baseline + ### Field + x_low, x_high, N_x = -300, 300, 81 + y_low, y_high, N_y = -300, 300, 81 + + ### Geometry + ants = [ + Antenna(x=-50,y=0,z=0,name="a"), + Antenna(x=50,y=0,z=0,name="b"), + ] + + if args.type == 'single-center': + tx = Antenna(x=-000,y=200,z=0,name="tx") + else: + tx = Antenna(x=-200,y=200,z=0,name="tx") + + elif args.type == 'square' or args.type == 'tri': # from grid definition + ### Field + x_low, x_high, N_x = -1800, 1800, 161 + y_low, y_high, N_y = -x_low, -x_high, N_x + + ### Geometry + tx = Antenna(x=-800,y=300,z=0,name="tx") + x_start, dx, ant_N_x = 0, 50, 2 + y_start, dy, ant_N_y = 0, dx, ant_N_x + + if args.type == 'square': # square grid + grid_func = square_grid + elif args.type == 'tri': # triangular + grid_func = triangular_grid + + grid = grid_func(dx=dx, dy=dy, N_x=ant_N_x, N_y=ant_N_y, x_start=x_start, y_start=y_start) + + + ants = [ Antenna(x=x,y=y,z=0,name=i) for i, (x,y) in enumerate(grid) ] + else: + ### Field + x_low, x_high, N_x = -400, 400, 161 + y_low, y_high, N_y = -x_low, -x_high, N_x + + ### Geometry + tx = Antenna(x=-300,y=300,z=0,name="tx") + ants = [ + Antenna(x=100,y=0,z=0,name="a"), + Antenna(x=0,y=-50,z=0,name="b"), + Antenna(x=+50,y=-180,z=0,name="c"), + Antenna(x=125,y=180,z=0,name="d"), + ] + + ### + ### Options + ### + plot_phase = args.plot_phase + ref_ant = args.ref_ant + + + ant_combi = antenna_combinations(ants, ref_ant=ref_ant) + + xs = np.linspace(x_low, x_high, N_x) + ys = np.linspace(y_low, y_high, N_y) + + xx, yy, val = calculate_field((xs, ys), tx, ant_combi, calculate_phase=plot_phase) + + mask = None + if False and plot_phase: + mask = abs(val) > np.pi + + + kwargs = {} + if plot_phase: + color_label='$\\sqrt{ \\sum \\left(\\varphi(x) - \\Delta \\varphi\\right)^2}$' + else: + color_label='$\\sqrt{ \\sum \\left(t(x) - \\Delta t\\right)^2}$ [ns]' + val *= 1e9 + kwargs['vmax'] = 100 + + ax = plot_field(tx, ants, xx, yy, val, ax=None, ref_ant=ref_ant, mask=mask, color_label=color_label, **kwargs) + +# if plot_phase: +# N_lowest = np.min(len(ant_combi)-1, 10) +# lowest_idx = np.argpartition(val, N_lowest)[:N_lowest] +# +# print(lowest_idx) +# print(val[lowest_idx]) +# print( list(zip(np.array(xx)[lowest_idx], np.array(yy)[lowest_idx])) ) + + if args.fname is not None: + plt.savefig(args.fname) + else: + plt.show()