#!/usr/bin/env python3 # vim: fdm=indent ts=4 """ Add a beacon measurement on top of the simulated airshower. """ import numpy as np import json import h5py import os.path as path from copy import deepcopy from earsim import REvent, Antenna, block_filter import lib tx_fname = 'tx.json' antennas_fname = 'antennas.hdf5' def write_tx_file(fname, tx, f_beacon): with open(fname, 'w') as fp: return json.dump( dict( f_beacon=f_beacon, tx=dict( x=tx.x, y=tx.y, z=tx.z, name=tx.name ) ), fp ) def read_tx_file(fname): with open(fname, 'r') as fp: data = json.load(fp) f_beacon = data['f_beacon'] tx = Antenna(**data['tx']) return tx, f_beacon def read_beacon_hdf5(fname, **h5ant_kwargs): with h5py.File(fname, 'r') as h5: tx = Antenna_from_h5ant(h5['tx'], traces_key=None) f_beacon = tx.attrs['f_beacon'] antennas = [] for k, h5ant in h5['antennas'].items(): ant = Antenna_from_h5ant(h5ant, **h5ant_kwargs) antennas.append(ant) return f_beacon, tx, antennas def Antenna_from_h5ant(h5ant, traces_key='traces', raise_exception=True, read_AxB=True, read_beacon_info=True): mydict = { k:h5ant.attrs.get(k) for k in ['x', 'y', 'z', 'name'] } ant = Antenna(**mydict) if h5ant.attrs: ant.attrs = {**h5ant.attrs} # Traces if traces_key is None: pass elif traces_key not in h5ant: if raise_exception: raise ValueError("Traces_key not in file") else: ant.t = deepcopy(h5ant[traces_key][0]) ant.Ex = deepcopy(h5ant[traces_key][1]) ant.Ey = deepcopy(h5ant[traces_key][2]) ant.Ez = deepcopy(h5ant[traces_key][3]) if len(h5ant[traces_key]) > 4: ant.beacon = deepcopy(h5ant[traces_key][4]) # E_AxB if read_AxB and 'E_AxB' in h5ant: ant.t_AxB = deepcopy(h5ant['E_AxB'][0]) ant.E_AxB = deepcopy(h5ant['E_AxB'][1]) # Beacons if read_beacon_info and 'beacon_info' in h5ant: h5beacon = h5ant['beacon_info'] beacon_info = {} for name in h5beacon.keys(): beacon_info[name] = dict(h5beacon[name].attrs) ant.beacon_info = beacon_info return ant def init_antenna_hdf5(fname, tx = None, f_beacon = None): with h5py.File(fname, 'w') as fp: if tx is not None or f_beacon is not None: tx_group = fp.create_group('tx') tx_attrs = tx_group.attrs if f_beacon is not None: tx_attrs['f_beacon'] = f_beacon if tx is not None: tx_attrs['x'] = tx.x tx_attrs['y'] = tx.y tx_attrs['z'] = tx.z tx_attrs['name'] = tx.name return fname def append_antenna_hdf5(fname, antenna, columns = [], name='traces', prepend_time=True, overwrite=True, attrs_dict={}): if not overwrite: raise NotImplementedError with h5py.File(fname, 'a') as fp: if 'antennas' in fp.keys(): if not overwrite: raise NotImplementedError group = fp['antennas'] else: group = fp.create_group('antennas') if antenna.name in group: if not overwrite: raise NotImplementedError h5ant = group[antenna.name] else: h5ant = group.create_group(antenna.name) h5ant_attrs = h5ant.attrs h5ant_attrs['x'] = antenna.x h5ant_attrs['y'] = antenna.y h5ant_attrs['z'] = antenna.z h5ant_attrs['name'] = antenna.name for k,v in attrs_dict.items(): h5ant_attrs[k] = v if name in h5ant: if not overwrite: raise NotImplementedError del h5ant[name] dset = h5ant.create_dataset(name, (len(columns) + 1*prepend_time, len(columns[0])), dtype='f') if prepend_time: dset[0] = antenna.t for i, col in enumerate(columns, 1*prepend_time): dset[i] = col if __name__ == "__main__": from os import path remake_tx = True fname = "ZH_airshower/mysim.sry" tx = Antenna(x=-500,y=0,z=0,name='tx') # m f_beacon = 51.53e-3 # GHz # Bandpass for E field blockfilter low_bp = 0.03 # GHz high_bp = 0.08 # GHz beacon_amplitudes = 1e-6*np.array([1e2, 0, 0]) # mu V/m beacon_radiate_rsq = True if beacon_radiate_rsq: # Move tx out, and magnify beacon_amplitude (at tx) tx = Antenna(x=-20e3,y=0,z=0,name='tx') # m dist = lib.distance(tx, Antenna(x=0, y=0, z=0)) ampl = max(1, dist**2) beacon_amplitudes *= ampl #### fname_dir = path.dirname(fname) tx_fname = path.join(fname_dir, tx_fname) antennas_fname = path.join(fname_dir, antennas_fname) if not path.isfile(tx_fname) or remake_tx: write_tx_file(tx_fname, tx, f_beacon) else: tx, f_beacon = read_tx_file(tx_fname) # read in antennas ev = REvent(fname) N_antennas = len(ev.antennas) # initialize hdf5 file init_antenna_hdf5(antennas_fname, tx, f_beacon) # make beacon per antenna for i, antenna in enumerate(ev.antennas): t0 = lib.distance(tx, antenna)/3e8 * 1e9 # ns beacon = lib.beacon_from(tx, antenna, f_beacon, antenna.t, t0=t0, c_light=np.inf, radiate_rsq=beacon_radiate_rsq) traces = np.array([antenna.Ex, antenna.Ey, antenna.Ez, beacon]) # add to relevant polarisation # and apply block filter dt = antenna.t[1] - antenna.t[0] for j, _ in enumerate(beacon_amplitudes): traces[j] = block_filter(traces[j], dt, low_bp, high_bp) append_antenna_hdf5( antennas_fname, antenna, traces, name='traces', prepend_time=True, attrs_dict=dict(t0=t0)) # Save E field in E_AxB E = [np.dot(ev.uAxB,[ex,ey,ez]) for ex,ey,ez in zip(traces[0], traces[1], traces[2])] append_antenna_hdf5( antennas_fname, antenna, [E], name='E_AxB', prepend_time=True) print("Antenna HDF5 file written as " + str(antennas_fname))