Extend figure to show phase gradient also

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]])
+            sc = ax.scatter(*antenna_locs, c=sigmas[ref_idx[i]], **scatter_kwargs)
-            fig.colorbar(sc, ax=ax, label='t' if time_unit is None else 't ['+time_unit+']')
+            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: