diff --git a/simulations/airshower_beacon_simulation/aa_generate_beacon.py b/simulations/airshower_beacon_simulation/aa_generate_beacon.py index 64ea898..aff1098 100755 --- a/simulations/airshower_beacon_simulation/aa_generate_beacon.py +++ b/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 - remake_tx = True - - 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 + # Noise properties + noise_sigma = args.noise_sigma # mu V/m set to None to ignore # Beacon properties - if False: # slowest beacon to be found: - f_beacon = 10e-3 # GHz - low_bp = 5e-3 # GHz - high_bp = 80e-3 # GHz - else: # original wanted beacon - f_beacon = 51.53e-3 # GHz + beacon_amplitudes = np.array(args.beacon_amplitudes) # mu V/m + beacon_radiate_rsq = args.beacon_rsq_decay # beacon_amplitude is repaired for distance to 0,0,0 - # Bandpass for E field blockfilter - low_bp = 30e-3 # 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 + # Beacon properties + f_beacon = args.beacon_frequency # GHz # 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 - noise_sigma = 1e-4 # set to None to ignore + # Transmitter + remake_tx = True - # Disable block_filter - if False: + tx = Antenna(x=0,y=0,z=0,name='tx') # m + 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 = 10000 # ns = 10us - pre_N = 2000 # ns = 2us + new_N = int(args.new_trace_length) + pre_N = int(args.pre_trace_length) 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])