m-thesis-introduction/simulations/airshower_beacon_simulation/aa_generate_beacon.py

#!/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 as copy

from earsim import REvent, Antenna
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, traces_key='traces'):
    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_key][0]
            antenna.Ex = ant[traces_key][1]
            antenna.Ey = ant[traces_key][2]
            antenna.Ez = ant[traces_key][3]
            if len(ant[traces_key]) > 4:
                antenna.beacon = ant[traces_key][4]

            if ant.attrs:
                antenna.attrs = {**ant.attrs}

            antennas.append(antenna)

    return f_beacon, tx, antennas

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

        for k,v in attrs_dict.items():
            ant_attrs[k] = v

        if name in ant_group:
            if not overwrite:
                raise NotImplementedError
            del ant_group[name]

        dset = ant_group.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

    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)

        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
        for j, _ in enumerate(beacon_amplitudes):
            E[j] += beacon_amplitudes[j]*beacon

        append_antenna_hdf5( antennas_fname, antenna, E, name='traces', prepend_time=True, attrs_dict=dict(t0=t0))

    print("Antenna HDF5 file written as " + str(antennas_fname))