ZH: script to setup the beacons

simulations/airshower_beacon_simulation/generate_beacon.py

@@ -0,0 +1,138 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import json
+import h5py
+
+from earsim import REvent, Antenna
+
+import lib
+
+tx_fname = 'tx.json'
+clocks_fname = 'clocks.csv'
+beacon_fname = 'beacons.hdf5'
+
+def write_beacon_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_beacon_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):
+    with h5py.File(fname, 'r') as h5:
+        tx_attrs = h5['tx'].attrs
+        f_beacon = tx_attrs.get('f_beacon')
+
+        mydict = { k:tx_attrs.get(k) for k in ['x', 'y', 'z', 'name'] }
+        tx = Antenna(**mydict)
+        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]
+
+            antennas.append(antenna)
+
+    return f_beacon, tx, antennas
+
+def init_beacon_hdf5(fname, tx, f_beacon):
+    with h5py.File(fname, 'w') as fp:
+        tx_group = fp.create_group('tx')
+        tx_group['f_beacon'] = f_beacon
+        tx_group['x'] = tx.x
+        tx_group['y'] = tx.y
+        tx_group['z'] = tx.z
+        tx_group['name'] = tx.name
+
+    return fname
+
+def append_beacon_hdf5(fname, antenna, *datasets, clock_offset=None, overwrite=False):
+
+    with h5py.File(fname, 'a') as fp:
+        if 'antennas' not in fp.keys():
+            group = fp.create_group('antennas')
+        else:
+            group = fp['antennas']
+
+        if overwrite:
+            raise NotImplementedError
+
+        ant_group = group.create_group(antenna.name)
+        ant_group['x'] = antenna.x
+        ant_group['y'] = antenna.y
+        ant_group['z'] = antenna.z
+        ant_group['name'] = antenna.name
+        if clock_offset is not None:
+            ant_group['clock_offset'] = clock_offset
+
+        dset = ant_group.create_dataset('traces', (len(datasets)+1, len(antenna.t)), dtype='f')
+
+        dset[0] = antenna.t
+        for i, mydset in enumerate(datasets,1):
+            dset[i] = mydset
+
+
+if __name__ == "__main__":
+    from os import path
+
+    max_clock_offset = 100# ns
+    remake_tx = False
+    remake_clock_offsets = False
+
+    seed = 12345
+    rng = np.random.default_rng(seed)
+
+    fname = "ZH_airshower/mysim.sry"
+    tx = Antenna(x=-500,y=0,z=0,name='tx')
+    f_beacon = 50e-3 # GHz
+
+    ####
+    fname_dir = path.dirname(fname)
+    tx_fname = path.join(fname_dir, tx_fname)
+    clocks_fname = path.join(fname_dir, clocks_fname)
+    beacon_fname = path.join(fname_dir, beacon_fname)
+
+    if not path.isfile(tx_fname) or remake_tx:
+        write_beacon_file(tx_fname, tx, f_beacon)
+    else:
+        tx, f_beacon = read_beacon_file(tx_fname)
+
+
+    # read in antennas
+    ev = REvent(fname)
+    N_antennas = len(ev.antennas)
+
+    if not path.isfile(clocks_fname) or remake_clock_offsets:
+        clock_offsets = max_clock_offset * (2*rng.uniform(size=N_antennas) - 1)
+        np.savetxt(clocks_fname, clock_offsets)
+    else:
+        clock_offsets = np.loadtxt(clocks_fname)
+
+    # make beacon per antenna
+    init_beacon_hdf5(beacon_fname, tx, f_beacon)
+    for i, antenna in enumerate(ev.antennas):
+        append_beacon_hdf5(
+                beacon_fname,
+                antenna,
+                lib.beacon_from(tx, antenna, f_beacon, antenna.t, t0=clock_offsets[i]),
+                clock_offset=clock_offsets[i],
+            )

simulations/airshower_beacon_simulation/lib.py

@@ -0,0 +1,29 @@
+"""
+Library for this simulation
+"""
+import numpy as np
+
+from earsim import Antenna
+
+def sine_beacon(f, t, t0=0, amplitude=1):
+    return amplitude * np.sin(2*np.pi*f*(t-t0))
+
+
+def distance(x1, x2):
+    """
+    Calculate the Euclidean distance between two locations x1 and x2
+    """
+
+    assert type(x1) in [Antenna]
+    x1 = np.array([x1.x, x1.y, x1.z])
+
+    assert type(x2) in [Antenna]
+    x2 = np.array([x2.x, x2.y, x2.z])
+
+    return np.sqrt( np.sum( (x1-x2)**2 ) )
+
+def beacon_from(tx, rx, f, t=0, t0=0, c_light=3e8, **kwargs):
+    dist = distance(tx,rx)/c_light
+    t0 = t0 + dist/c_light
+
+    return sine_beacon(f, t, t0=t0, **kwargs)