Merge branch 'argparse' into main

simulations/airshower_beacon_simulation/Makefile

@@ -1,22 +1,32 @@
 .PHONY: all
 
-all: beacon clocks
+
+FIG_DIR := ./figures
+
+all: beacon clocks phases findks reconstruct
 
 
 beacon:
-	./aa_generate_beacon.py
+	./aa_generate_beacon.py > figures/aa.log
+	./ab_modify_clocks.py 0 > figures/ab.log
+	./ac_show_signal_to_noise.py --no-show-plots --fig-dir=${FIG_DIR}
 
 clocks:
-	./ab_modify_clocks.py
+	./ab_modify_clocks.py 15 --gaussian
 
 phases:
-	./ba_measure_beacon_phase.py
+	./ba_measure_beacon_phase.py --no-show-plots --fig-dir=${FIG_DIR}
-	./bb_measure_true_phase.py
+	./bb_measure_true_phase.py --no-show-plots --fig-dir=${FIG_DIR}
-	./bc_baseline_phase_deltas.py
+	./bc_baseline_phase_deltas.py --no-show-plots --fig-dir=${FIG_DIR}
-	./bd_antenna_phase_deltas.py
+	./bd_antenna_phase_deltas.py --no-show-plots --fig-dir=${FIG_DIR}
 
-period_multiples:
+findks:
-	./ca_periods_from_shower.py
+	./ca_period_from_shower.py --no-show-plots --fig-dir=${FIG_DIR}
+	./cb_report_measured_antenna_time_offsets.py --no-show-plots --fig-dir=${FIG_DIR}
+
+reconstruct:
+	./da_reconstruction.py --no-show-plots --fig-dir=${FIG_DIR}
+	./db_longitudinal_figure.py --no-show-plots --fig-dir=${FIG_DIR}
 
 dist-clean:
 	rm -f ZH_airshower/antennas.hdf5

simulations/airshower_beacon_simulation/README.md

@@ -30,4 +30,4 @@ Steps:
 
  3. Find k\*2\\pi phase offsets (periods) ([./ca_periods_from_showers.py])
 
- 4. Rewrite clocks
+ 4. Rewrite clocks and do interferometric reconstruction ([./da_reconstruction.py])

simulations/airshower_beacon_simulation/aa_generate_beacon.py

@@ -241,36 +241,44 @@ def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods
 
 if __name__ == "__main__":
     from os import path
+    from argparse import ArgumentParser
+
+    parser = ArgumentParser()
+    parser.add_argument('-n', '--noise-sigma', type=float, default=1e3, help='in [muV/m]')
+    parser.add_argument('-f', '--beacon-frequency', type=float, default=51.53e-3, help='The beacon\'s frequency [GHz]')
+
+    # Beacon Properties
+    parser.add_argument('-a', '--beacon-amplitudes', type=float, nargs=3, default=[1e3, 0, 0], help='in [muV/m]')
+    parser.add_argument('-d', '--beacon-rsq-decay', type=bool, default=True, help='Use Beacon amplitudes at 0,0,0')
+
+    # Bandpass
+    parser.add_argument('-p', '--use-passband', type=bool, default=True)
+    parser.add_argument('-l', '--passband-low',  type=float, default=30e-3, help='Lower frequency [GHz] of the passband filter. (set -1 for np.inf)')
+    parser.add_argument('-u', '--passband-high', type=float, default=80e-3, help='Upper frequency [GHz] of the passband filter. (set -1 for np.inf)')
+
+    # Trace length modification
+    parser.add_argument('-N', '--new-trace-length', type=float, help='resize airshower trace', default=1e4)
+    parser.add_argument('-P', '--pre-trace-length', type=float, help='amount of trace to prepend the airshower when resizing', default=2e3)
+
+    args = parser.parse_args()
+
+    ##
+    ## End of ArgumentParsing
+    ##
 
     rng = np.random.default_rng()
 
     fname = "ZH_airshower/mysim.sry"
 
-    # Transmitter
+    # Noise properties
-    remake_tx = True
+    noise_sigma = args.noise_sigma # mu V/m set to None to ignore
-
-    tx = Antenna(x=-2e3,y=0,z=0,name='tx') # m
-    if False:
-        # Move tx out a long way
-        tx.x, tx.y = -75e3, 75e3 # m
-    elif False:
-        # Move it to 0,0,0 (among the antennas)
-        tx.x, tx.y = 0, 0 #m
 
     # Beacon properties
-    if False: # slowest beacon to be found:
+    beacon_amplitudes = np.array(args.beacon_amplitudes) # mu V/m
-        f_beacon =  10e-3 # GHz
+    beacon_radiate_rsq = args.beacon_rsq_decay # beacon_amplitude is repaired for distance to 0,0,0
-        low_bp = 5e-3 # GHz
-        high_bp = 80e-3 # GHz
-    else: # original wanted beacon
-        f_beacon = 51.53e-3 # GHz
 
-        # Bandpass for E field blockfilter
+    # Beacon properties
-        low_bp = 30e-3 # GHz
+    f_beacon = args.beacon_frequency # GHz
-        high_bp = 80e-3 # GHz
-
-    beacon_amplitudes = 1e-6*np.array([1e3, 0, 0]) # mu V/m
-    beacon_radiate_rsq = True # beacon_amplitude is repaired for distance to 0,0,0
 
     # modify beacon power to be beacon_amplitude at 0,0,0
     if beacon_radiate_rsq:
@@ -280,11 +288,23 @@ if __name__ == "__main__":
 
         beacon_amplitudes *= ampl
 
-    # Noise properties
+    # Transmitter
-    noise_sigma = 1e-4 # set to None to ignore
+    remake_tx = True
 
-    # Disable block_filter
+    tx = Antenna(x=0,y=0,z=0,name='tx') # m
-    if False:
+    if True:
+        # Move tx out a long way
+        tx.x, tx.y = -75e3, 75e3 # m
+    elif False:
+        # Move it to 0,0,0 (among the antennas)
+        tx.x, tx.y = 0, 0 #m
+
+    # Bandpass for E field blockfilter
+    low_bp  = args.passband_low  if args.passband_low  >= 0 else np.inf # GHz
+    high_bp = args.passband_high if args.passband_high >= 0 else np.inf # GHz
+
+    # Enable/Disable block_filter
+    if not args.use_passband:
         block_filter = lambda x, dt, low, high: x
 
     ####
@@ -300,8 +320,9 @@ if __name__ == "__main__":
 
 
     print("Beacon amplitude at tx [muV/m]:", beacon_amplitudes)
+    print("Beacon amplitude at 0,0,0 [muV/m]:", beacon_amplitudes/ampl)
     print("Tx location:", [tx.x, tx.y, tx.z])
-    print("Noise sigma:", noise_sigma)
+    print("Noise sigma [muV/m]:", noise_sigma)
 
     # read in antennas
     ev = REvent(fname)
@@ -315,11 +336,10 @@ if __name__ == "__main__":
         if i%10 == 0:
             print(f"Beaconed antenna {i} out of", len(ev.antennas))
 
-        if not False: # modify trace lengths
+        if args.new_trace_length: # modify trace lengths
             N_samples = len(antenna.t)
-            #new_N = 2*N_samples
+            new_N = int(args.new_trace_length)
-            new_N = 10000 # ns = 10us
+            pre_N = int(args.pre_trace_length)
-            pre_N = 2000 # ns = 2us
             after_N = new_N - pre_N
 
             dt = antenna.t[1] - antenna.t[0]
@@ -335,12 +355,12 @@ if __name__ == "__main__":
             if i%10 == 0:
                 print(f"Modified trace lengths by {pre_N},{after_N-N_samples}")
 
-        beacon = lib.beacon_from(tx, antenna, f_beacon, antenna.t, c_light=c_light, radiate_rsq=beacon_radiate_rsq)
+        beacon = 1e-6 * lib.beacon_from(tx, antenna, f_beacon, antenna.t, c_light=c_light, radiate_rsq=beacon_radiate_rsq) # mu V/m
 
         # noise realisation
         noise_realisation = 0
         if noise_sigma is not None:
-            noise_realisation = rng.normal(0, noise_sigma, size=len(beacon))
+            noise_realisation = 1e-6 * rng.normal(0, noise_sigma, size=len(beacon)) # mu V/m
 
         # Collect all data to be saved (with the first 3 values the E fields)
         traces = np.array([antenna.Ex, antenna.Ey, antenna.Ez, beacon, noise_realisation])

simulations/airshower_beacon_simulation/ab_modify_clocks.py

@@ -19,7 +19,17 @@ if __name__ == "__main__":
     from os import path
     import sys
 
-    max_clock_offset = 100 if len(sys.argv) < 2 else float(sys.argv[1]) # ns
+    from argparse import ArgumentParser
+
+    parser = ArgumentParser()
+    parser.add_argument('max_clock_offset', nargs='?', type=float, default=25, help='(Default: %(default)d)')
+    parser.add_argument('--uniform', action='store_const', const='uniform', dest='dist_type')
+    parser.add_argument('--gaussian',   action='store_const', const='gauss', dest='dist_type')
+    parser.set_defaults(dist_type='gauss')
+
+    args = parser.parse_args()
+
+    max_clock_offset = args.max_clock_offset # ns
     remake_clock_offsets = True
 
     seed = 12345
@@ -54,7 +64,7 @@ if __name__ == "__main__":
         if True:
             print(f"Modifying clocks upto {max_clock_offset}ns.")
             clock_offsets = np.zeros( N_antennas )
-            if not True: # uniform
+            if args.dist_type == 'uniform': # uniform
                 print("Uniform distribution")
                 clock_offsets = max_clock_offset * (2*rng.uniform(size=N_antennas) - 1)
             else: # normal

simulations/airshower_beacon_simulation/ac_show_signal_to_noise.py

@@ -15,7 +15,6 @@ from earsim import REvent, block_filter
 import aa_generate_beacon as beacon
 import lib
 
-
 if __name__ == "__main__":
     from os import path
     import sys
@@ -24,10 +23,20 @@ if __name__ == "__main__":
     if os.name == 'posix' and "DISPLAY" not in os.environ:
         matplotlib.use('Agg')
 
+    from scriptlib import MyArgumentParser
+    parser = MyArgumentParser()
+
+    # Bandpass
+    parser.add_argument('-p', '--use-passband', type=bool, default=True)
+    parser.add_argument('-l', '--passband-low',  type=float, default=30e-3, help='Lower frequency [GHz] of the passband filter. (set -1 for np.inf)')
+    parser.add_argument('-u', '--passband-high', type=float, default=80e-3, help='Upper frequency [GHz] of the passband filter. (set -1 for np.inf)')
+
+    args = parser.parse_args()
+
     fname = "ZH_airshower/mysim.sry"
 
-    fig_dir = "./figures/"
+    fig_dir = args.fig_dir
-    show_plots = not False
+    show_plots = args.show_plots
 
     ####
     fname_dir = path.dirname(fname)
@@ -44,12 +53,15 @@ if __name__ == "__main__":
 
     # general properties
     dt = antennas[0].t[1] - antennas[0].t[0] # ns
-    pb = lib.passband(30e-3, 80e-3) # GHz
     beacon_pb = lib.passband(f_beacon-1e-3, f_beacon+1e-3) # GHz
-
     beacon_amp = np.max(txdata['amplitudes'])# mu V/m
 
-    if True: # Apply filter to raw beacon/noise to compare with Filtered Traces
+    # General Bandpass
+    low_bp  = args.passband_low  if args.passband_low  >= 0 else np.inf # GHz
+    high_bp = args.passband_high if args.passband_high >= 0 else np.inf # GHz
+    pb = lib.passband(low_bp, high_bp) # GHz
+
+    if args.use_passband: # Apply filter to raw beacon/noise to compare with Filtered Traces
         myfilter = lambda x: block_filter(x, dt, pb[0], pb[1])
     else: # Compare raw beacon/noise with Filtered Traces
         myfilter = lambda x: x

simulations/airshower_beacon_simulation/ba_measure_beacon_phase.py

@@ -21,14 +21,19 @@ if __name__ == "__main__":
     if os.name == 'posix' and "DISPLAY" not in os.environ:
         matplotlib.use('Agg')
 
+    from scriptlib import MyArgumentParser
+    parser = MyArgumentParser()
+    parser.add_argument('--no-use-AxB', dest='use_AxB_trace', action='store_false')
+    args = parser.parse_args()
+
     f_beacon_band = (49e-3,55e-3) #GHz
     allow_frequency_fitting = False
     read_frequency_from_file = True
 
-    use_AxB_trace = True if len(sys.argv) < 2 else bool(int(sys.argv[1]))
+    use_AxB_trace = args.use_AxB_trace
     use_beacon_trace = True # only applicable if AxB = False
 
-    show_plots = True
+    show_plots = args.show_plots
 
     fname = "ZH_airshower/mysim.sry"
 
@@ -38,7 +43,7 @@ if __name__ == "__main__":
     fname_dir = path.dirname(fname)
     antennas_fname = path.join(fname_dir, beacon.antennas_fname)
 
-    fig_dir = "./figures" # set None to disable saving
+    fig_dir = args.fig_dir # set None to disable saving
 
     if not path.isfile(antennas_fname):
         print("Antenna file cannot be found, did you try generating a beacon?")

simulations/airshower_beacon_simulation/bb_measure_true_phase.py

@@ -18,9 +18,13 @@ if __name__ == "__main__":
     if os.name == 'posix' and "DISPLAY" not in os.environ:
         matplotlib.use('Agg')
 
+    from scriptlib import MyArgumentParser
+    parser = MyArgumentParser()
+    args = parser.parse_args()
+
     fname = "ZH_airshower/mysim.sry"
     c_light = lib.c_light
-    show_plots = True
+    show_plots = args.show_plots
 
     remove_absolute_phase_offset_first_antenna = True # takes precedence
     remove_absolute_phase_offset_minimum = True
@@ -29,7 +33,7 @@ if __name__ == "__main__":
     fname_dir = path.dirname(fname)
     antennas_fname = path.join(fname_dir, beacon.antennas_fname)
 
-    fig_dir = "./figures" # set None to disable saving
+    fig_dir = args.fig_dir # set None to disable saving
 
     if not path.isfile(antennas_fname):
         print("Antenna file cannot be found, did you try generating a beacon?")

simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py

@@ -19,17 +19,21 @@ if __name__ == "__main__":
     if os.name == 'posix' and "DISPLAY" not in os.environ:
         matplotlib.use('Agg')
 
+    from scriptlib import MyArgumentParser
+    parser = MyArgumentParser()
+    args = parser.parse_args()
+
     fname = "ZH_airshower/mysim.sry"
     c_light = 3e8*1e-9
-    show_plots = True
+    show_plots = args.show_plots
-    ref_ant_id = None if not True else [50] # leave None for all baselines
+    ref_ant_id = None if True else [50] # leave None for all baselines
 
     ####
     fname_dir = path.dirname(fname)
     antennas_fname = path.join(fname_dir, beacon.antennas_fname)
     time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname
 
-    fig_dir = "./figures" # set None to disable saving
+    fig_dir = args.fig_dir # set None to disable saving
 
     if not path.isfile(antennas_fname):
         print("Antenna file cannot be found, did you try generating a beacon?")

simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py

@@ -21,16 +21,20 @@ if __name__ == "__main__":
     if os.name == 'posix' and "DISPLAY" not in os.environ:
         matplotlib.use('Agg')
 
+    from scriptlib import MyArgumentParser
+    parser = MyArgumentParser()
+    args = parser.parse_args()
+
     fname = "ZH_airshower/mysim.sry"
 
-    show_plots = True
+    show_plots = args.show_plots
     ref_ant_id = None # leave None for all baselines
 
     ####
     fname_dir = path.dirname(fname)
     antennas_fname = path.join(fname_dir, beacon.antennas_fname)
     time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname
-    fig_dir = "./figures" # 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)
 

simulations/airshower_beacon_simulation/ca_period_from_shower.py

@@ -149,11 +149,15 @@ if __name__ == "__main__":
 
     atm = AtmoCal()
 
+    from scriptlib import MyArgumentParser
+    parser = MyArgumentParser(default_fig_dir="./figures/periods_from_shower_figures/")
+    args = parser.parse_args()
+
     fname = "ZH_airshower/mysim.sry"
 
-    fig_dir = "./figures/periods_from_shower_figures/"
+    fig_dir = args.fig_dir
     fig_subdir = path.join(fig_dir, 'shifts/')
-    show_plots = False
+    show_plots = args.show_plots
 
     allowed_ks = [ -2, -1, 0, 1, 2]
     Xref = 400

simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py

@@ -18,10 +18,14 @@ if __name__ == "__main__":
     if os.name == 'posix' and "DISPLAY" not in os.environ:
         matplotlib.use('Agg')
 
+    from scriptlib import MyArgumentParser
+    parser = MyArgumentParser()
+    args = parser.parse_args()
+
     fname = "ZH_airshower/mysim.sry"
 
-    fig_dir = "./figures" # set None to disable saving
+    fig_dir = args.fig_dir # set None to disable saving
-    show_plots = not False
+    show_plots = args.show_plots
 
     ####
     fname_dir = path.dirname(fname)

simulations/airshower_beacon_simulation/da_reconstruction.py

@@ -28,11 +28,15 @@ if __name__ == "__main__":
 
     atm = AtmoCal()
 
+    from scriptlib import MyArgumentParser
+    parser = MyArgumentParser()
+    args = parser.parse_args()
+
     fname = "ZH_airshower/mysim.sry"
 
-    fig_dir = "./figures/"
+    fig_dir = args.fig_dir
     fig_subdir = path.join(fig_dir, 'reconstruction')
-    show_plots = False
+    show_plots = args.show_plots
 
     ####
     fname_dir = path.dirname(fname)

simulations/airshower_beacon_simulation/db_longitudinal_figure.py

@@ -24,10 +24,14 @@ if __name__ == "__main__":
     if os.name == 'posix' and "DISPLAY" not in os.environ:
         matplotlib.use('Agg')
 
+    from scriptlib import MyArgumentParser
+    parser = MyArgumentParser()
+    args = parser.parse_args()
+
     fname = "ZH_airshower/mysim.sry"
 
-    fig_dir = "./figures/"
+    fig_dir = args.fig_dir
-    show_plots = False
+    show_plots = args.show_plots
 
     ####
     fname_dir = path.dirname(fname)

simulations/airshower_beacon_simulation/scriptlib.py

@@ -0,0 +1,24 @@
+"""
+Some preconfigured ArgumentParser
+"""
+
+from argparse import ArgumentParser
+
+def MyArgumentParser(default_fig_dir='./figures', default_show_plots=False, **kwargs):
+    """
+    A somewhat preconfigured ArgumentParser
+
+    Set show_plots=True to by default enable showing plots.
+    Likewise, set fig_dir=None to by default disable saving figures.
+    """
+    parser = ArgumentParser(**kwargs)
+
+    # Whether to show plots
+    parser.add_argument('--show-plots', action='store_true')
+    parser.add_argument('--no-show-plots', dest='show-plots', action='store_false')
+    parser.set_defaults(show_plots=default_show_plots)
+
+    # Figures directory
+    parser.add_argument('--fig-dir', type=str, default=default_fig_dir)
+
+    return parser