ZH: backup traces in hdf5

simulations/airshower_beacon_simulation/generate_beacon.py

@@ -4,11 +4,9 @@ import numpy as np
 import json
 import h5py
 import os.path as path
-
 from copy import deepcopy as copy
 
 from earsim import REvent, Antenna
-
 import lib
 
 tx_fname = 'tx.json'
@@ -46,14 +44,15 @@ def read_beacon_hdf5(fname):
         mydict = { k:tx_attrs.get(k) for k in ['x', 'y', 'z', 'name'] }
         tx = Antenna(**mydict)
 
-        print(f_beacon, tx)
-
         antennas = []
         for k, ant in h5['antennas'].items():
             mydict = { k:ant.attrs.get(k) for k in ['x', 'y', 'z', 'name'] }
             antenna = Antenna(**mydict)
             antenna.t = ant['traces'][0]
-            antenna.beacon = ant['traces'][1]
+            antenna.Ex = ant['traces'][1]
+            antenna.Ey = ant['traces'][2]
+            antenna.Ez = ant['traces'][3]
+            antenna.beacon = ant['traces'][4]
 
             antennas.append(antenna)
 
@@ -76,69 +75,45 @@ def init_antenna_hdf5(fname, tx = None, f_beacon = None):
 
     return fname
 
-def append_antenna_hdf5(fname, antenna, datasets = [], overwrite=False):
+def append_antenna_hdf5(fname, antenna, columns = [], name='traces', prepend_time=True, overwrite=True):
+    if not overwrite:
+        raise NotImplementedError
+
     with h5py.File(fname, 'a') as fp:
-        if 'antennas' not in fp.keys():
+        if 'antennas' in fp.keys():
-            group = fp.create_group('antennas')
+            if not overwrite:
-        else:
+                raise NotImplementedError
             group = fp['antennas']
+        else:
+            group = fp.create_group('antennas')
 
-        if overwrite:
-            raise NotImplementedError
 
-        ant_group = group.create_group(antenna.name)
+        if antenna.name in group:
+            if not overwrite:
+                raise NotImplementedError
+            ant_group = group[antenna.name]
+        else:
+            ant_group = group.create_group(antenna.name)
+
         ant_attrs = ant_group.attrs
         ant_attrs['x'] = antenna.x
         ant_attrs['y'] = antenna.y
         ant_attrs['z'] = antenna.z
         ant_attrs['name'] = antenna.name
 
-        dset = ant_group.create_dataset('traces', (len(datasets)+1, len(antenna.t)), dtype='f')
+        if name in ant_group:
+            if not overwrite:
+                raise NotImplementedError
+            del ant_group[name]
 
-        dset[0] = antenna.t
+        dset = ant_group.create_dataset(name, (len(columns) + 1*prepend_time, len(columns[0])), dtype='f')
-        for i, mydset in enumerate(datasets,1):
-            dset[i] = mydset
 
-def BeaconedREvent(fname, fname_beacon=None, beacon_amplitude=(1,0,0)):
+        if prepend_time:
+            dset[0] = antenna.t
 
-    if fname_beacon is None:
+        for i, col in enumerate(columns, 1*prepend_time):
-        fname_beacon = path.join(path.dirname(fname), beacon_fname)
+            dset[i] = col
-
-    f_beacon, tx, antennas = read_beacon_hdf5(fname_beacon)
-
-    ev = REvent(fname)
-    ev.tx = tx
-    ev.tx.f_beacon = f_beacon
-
-    merge_beacon_into_REvent(antennas, ev, beacon_amplitude=beacon_amplitude)
     
-    return ev
-    
-def merge_beacon_into_REvent(beacons, ev, beacon_amplitude=(1,0,0)):
-    assert len(beacons) == len(ev.antennas)
-
-    if not hasattr(beacon_amplitude, '__len__'):
-        beacon_amplitude = np.repeat(beacon_amplitude,3)
-
-    assert len(beacon_amplitude) == 3
-
-    for i, beacon in enumerate(beacons):
-        ant = ev.antennas[i]
-        assert len(beacon.t) == len(ant.t)
-
-        # copy original Efields
-        ev.antennas[i].orig_Ex = copy(ev.antennas[i].Ex)
-        ev.antennas[i].orig_Ey = copy(ev.antennas[i].Ey)
-        ev.antennas[i].orig_Ez = copy(ev.antennas[i].Ez)
-
-        ev.antennas[i].beacon = beacon.beacon
-        if beacon_amplitude[0] != 0:
-            ev.antennas[i].Ex += beacon_amplitude[0]*beacon.beacon
-        if beacon_amplitude[1] != 0:
-            ev.antennas[i].Ey += beacon_amplitude[1]*beacon.beacon
-        if beacon_amplitude[2] != 0:
-            ev.antennas[i].Ez += beacon_amplitude[2]*beacon.beacon
-
 if __name__ == "__main__":
     from os import path
 
@@ -148,6 +123,8 @@ if __name__ == "__main__":
     tx = Antenna(x=-500,y=0,z=0,name='tx')
     f_beacon = 50e-3 # GHz
 
+    beacon_amplitudes = 1e-6*np.array([1e2, 0, 0]) # mu V/m
+
     ####
     fname_dir = path.dirname(fname)
     tx_fname = path.join(fname_dir, tx_fname)
@@ -169,11 +146,13 @@ if __name__ == "__main__":
     for i, antenna in enumerate(ev.antennas):
         beacon = lib.beacon_from(tx, antenna, f_beacon, antenna.t)
 
-        E = [antenna.Ex, antenna.Ey, antenna.Ez ]
+        E = np.array([antenna.Ex, antenna.Ey, antenna.Ez, beacon])
+        append_antenna_hdf5( antennas_fname, antenna, E, name='orig_traces', prepend_time=True)
 
         # add to relevant polarisation
-        E[0] += beacon
+        for i, _ in enumerate(beacon_amplitudes):
+            E[i] += beacon_amplitudes[i]*beacon
 
-        append_antenna_hdf5( antennas_fname, antenna, (E[0], E[1], E[2] ))
+        append_antenna_hdf5( antennas_fname, antenna, E, name='traces', prepend_time=True)
 
     print("Antenna HDF5 file written as " + str(antennas_fname))

simulations/airshower_beacon_simulation/modify_clocks.py

@@ -45,7 +45,7 @@ if __name__ == "__main__":
     # read in antennas
     with h5py.File(antennas_fname, 'a') as fp:
         if 'antennas' not in fp.keys():
-            print("Antenna file corrupted?")
+            print("Antenna file corrupted? no antennas")
             sys.exit(1)
 
         group = fp['antennas']
@@ -59,6 +59,10 @@ if __name__ == "__main__":
             ant_group = group[name]
             clk_offset = clock_offsets[i]
 
+            if 'traces' not in ant_group.keys():
+                print(f"Antenna file corrupted? no 'traces' in {name}")
+                sys.exit(1)
+
             ant_attrs = ant_group.attrs
             if 'clock_offset' in ant_attrs:
                 tmp_offset = ant_attrs['clock_offset']