figure: 4antennasetup + iter clock fix

Also showcases how to remove time offsets iteratively

figures/beacon/src/four_antenna_setup.py

@@ -0,0 +1,228 @@
+#!/usr/bin/env python3
+
+__doc__ = \
+"""
+Create a figure showing the timing and geometry of 3+1 antennas,
+and the triangles between them.
+The additional antenna is to show that baselines are shared between
+triangles.
+
+In code (and on stdout), the antennas have time offsets which can be determined iteratively up to an overall offset.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+from itertools import combinations
+
+c_light = 3e8 # m/s
+
+def distance(x1, x2):
+    """
+    Calculate the Euclidean distance between two locations x1 and x2
+    """
+
+    assert type(x1) in [Antenna]
+    x1 = np.array([x1.x, x1.y, x1.z])
+
+    assert type(x2) in [Antenna]
+    x2 = np.array([x2.x, x2.y, x2.z])
+
+    return np.sqrt( np.sum( (x1-x2)**2 ) )
+
+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
+
+
+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.t0 = t0
+        self.name = name
+
+        self.offsets = []
+
+    def __repr__(self):
+        cls = self.__class__.__name__
+
+        return f'{cls}(x={self.x!r},y={self.y!r},z={self.z!r},t0={self.t0!r},name={self.name!r})'
+
+def antenna_triangles(antennas):
+    return combinations(antennas, 3)
+
+def antenna_baselines(antennas):
+    return combinations(antennas, 2)
+
+def plot_four_antenna_geometry(tx, ants, extra_ant=None, ax=None, line_kwargs={}, scatter_kwargs={}, scatter_zorder=5):
+
+    default_line_kwargs = dict( color='grey', lw=3, alpha=0.7)
+    default_scatter_kwargs = dict( color='grey', s=200)
+
+    if ax is None:
+        ax = plt.gca()
+        ax.set_aspect('equal')
+        ax.set_xlabel("x")
+        ax.set_ylabel("y")
+
+    line_kwargs = {**default_line_kwargs, **line_kwargs}
+    scatter_kwargs = {**default_scatter_kwargs, **scatter_kwargs}
+
+    # Plot Antennas + Tx
+    for i, ant in enumerate([tx] + ants + [extra_ant]):
+        ax.scatter(ant.x, ant.y, zorder=scatter_zorder, **scatter_kwargs)
+        ax.annotate(ant.name, (ant.x, ant.y), ha='center', va='center',zorder=scatter_zorder)
+
+    # Lines connecting Tx and ants
+    tmp_line_kwargs = line_kwargs
+    for ant in ants:
+        ax.plot([tx.x, ant.x], [tx.y, ant.y], **tmp_line_kwargs)
+
+    # Lines due to all Antennas (including extra_ant)
+    line_offset = 0.08*np.array([1,1])
+    for i, ant_triangle in enumerate(antenna_triangles(ants + [extra_ant])):
+        tmp_line_kwargs['color'] = None
+        tmp_line_kwargs['linestyle'] = '--'
+        tmp_line_kwargs['alpha'] = 0.4
+        for j, ant in enumerate(antenna_baselines(ant_triangle)):
+            a, b = ant[0], ant[1]
+            if j == 1: # fix ordering
+                a, b == b, a
+
+            dx, dy = (i-1)*line_offset
+    
+            l = ax.plot([ a.x+dx, b.x+dx], [a.y+dy, b.y+dy], **tmp_line_kwargs)
+            line_kwargs['color'] = l[0].get_color()
+
+    # Lines internal to ants triangle
+    tmp_line_kwargs = line_kwargs
+    tmp_line_kwargs['color'] = 'green'
+    tmp_line_kwargs['alpha'] = 0.7
+    for j, ant_pair in enumerate(combinations(ants,2)):
+        a, b = ant_pair[0], ant_pair[1]
+        if j == 1: # fix ordering
+            a, b = b, a
+
+        ax.plot([ a.x, b.x], [a.y, b.y], **tmp_line_kwargs)
+
+    return ax
+
+if __name__ == "__main__":
+    use_phase = False
+    correct_time = True
+
+    tx = Antenna(name="T", x=-8, y=2, t0=0)
+
+    ants = [
+        Antenna(name='1', x=0, y= 0, t0=1 ),
+        Antenna(name='2', x=2, y=-3, t0=4 ),
+        Antenna(name='3', x=1, y= 3, t0=10 ),
+        ]
+
+    extra_ant = Antenna(name='4', x=4, y=-1, t0=-6)
+
+
+    if False:#correct_time: # taken from the output of this script
+        """
+        Antenna Triangle(1,2,3)
+        j=0: 1,2
+        j=1: 3,1
+        j=2: 2,3
+        sigmas: [-2.99999999  9.         -6.00000001]
+        sigmas sum: -8.881784197001252e-16
+        (Antenna(x=0,y=0,z=0,t0=1,name='1'), Antenna(x=2,y=-3,z=0,t0=4,name='2'), Antenna(x=1,y=3,z=0,t0=10,name='3')) : [0.         2.99999999 9.        ]
+        Antenna Triangle(1,2,4)
+        sigmas: [-2.99999999 -7.00000001 10.        ]
+        sigmas sum: 0.0
+        (Antenna(x=0,y=0,z=0,t0=1,name='1'), Antenna(x=2,y=-3,z=0,t0=4,name='2'), Antenna(x=4,y=-1,z=0,t0=-6,name='4')) : [ 0.          2.99999999 -7.00000001]
+        Antenna Triangle(1,3,4)
+        sigmas: [-9.         -7.00000001 16.00000001]
+        sigmas sum: 0.0
+        (Antenna(x=0,y=0,z=0,t0=1,name='1'), Antenna(x=1,y=3,z=0,t0=10,name='3'), Antenna(x=4,y=-1,z=0,t0=-6,name='4')) : [ 0.          9.         -7.00000001]
+        Antenna Triangle(2,3,4)
+        sigmas: [ -6.00000001 -10.          16.00000001]
+        sigmas sum: 0.0
+        (Antenna(x=2,y=-3,z=0,t0=4,name='2'), Antenna(x=1,y=3,z=0,t0=10,name='3'), Antenna(x=4,y=-1,z=0,t0=-6,name='4')) : [  0.           6.00000001 -10.
+        """
+        print("Running with pre-corrected times")
+        ants[1].t0 -= 2.99999999
+        ants[2].t0 -= 9
+        extra_ant.t0 -= -7
+
+    ax = plot_four_antenna_geometry(tx, ants, extra_ant=extra_ant)
+
+    fig = ax.get_figure()
+    
+    ### Lol, show my calculations are right
+    for i, triangle in enumerate(antenna_triangles(ants + [extra_ant])):
+        print("Antenna Triangle({},{},{})".format(*[ant.name for ant in triangle]))
+
+        sigma = np.zeros((3))
+        for j, ant_pair in enumerate(antenna_baselines(triangle)):
+            a, b = ant_pair[0], ant_pair[1]
+            if j == 1: # fix ordering
+                a, b = b, a
+            
+            if i == 0: # print sigma pairing for first triangle
+                print('j={}: {},{}'.format(j, a.name, b.name))
+
+            phys_Dt = (distance(tx, a) - distance(tx, b))/c_light
+            meas_Dt = a.t0 - b.t0
+
+            if use_phase:
+                f_beacon = 50e6 # Hz
+                to_phase = lambda t: phase_mod(2*np.pi*f_beacon*t)
+
+                phys_Dt = to_phase(phys_Dt)
+                meas_Dt = to_phase(meas_Dt)
+                
+            sigma[j] = meas_Dt - phys_Dt
+            
+            if False:
+                print(
+                    "Dt'_{},{} = ".format(a.name, b.name)
+                    + "{}".format(meas_Dt)
+                    + " = {} - {}".format(a.t0, b.t0)
+                    )
+                print(
+                    "Dt_{},{} = ".format(a.name, b.name)
+                    + "{}".format(phys_Dt)
+                    + " = {} - {}".format(distance(tx, a), distance(tx, b))
+                    )
+                
+                print(
+                    "sigma_{},{} = ".format(a.name, b.name)
+                    + "{}".format(sigma[j])
+                    )
+        print("sigmas:", sigma)
+        if use_phase:
+            print("sigmas sum:", phase_mod(np.sum(sigma)))
+        else:
+            print("sigmas sum:", np.sum(sigma))
+
+        # Take the first antenna as reference
+        ref_idx = 0
+        ref_sigma = sigma[ref_idx]
+        sigma = sigma - ref_sigma
+
+        ant_sigma = 1/3*np.array([0, sigma[1] + sigma[2], 2*sigma[1] - sigma[2]])
+
+        for i in [1,2]:
+            triangle[i].offsets.append(-ant_sigma[i])
+            if correct_time:
+                triangle[i].backup_t0 = triangle[i].t0
+                triangle[i].t0 += -ant_sigma[i]
+
+    if correct_time:
+        print(ants + [extra_ant])
+        for i, ant in enumerate(ants + [extra_ant]):
+            print(i, ant.name, ant.offsets)
+
+    plt.show()