ZH: cli arguments for beacon generation

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])