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

144 lines
4.6 KiB
Python
Raw Normal View History

#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
import numpy.fft as ft
import aa_generate_beacon as beacon
from view_orig_ant0 import plot_antenna_geometry
import lib
from earsim import Antenna
if __name__ == "__main__":
from os import path
import sys
import matplotlib
import os
if os.name == 'posix' and "DISPLAY" not in os.environ:
matplotlib.use('Agg')
from scriptlib import MyArgumentParser
parser = MyArgumentParser()
parser.add_argument('ant_idx', default=[72], nargs='*', help='Antenna Indices')
parser.add_argument('-p', '--polarisations', choices=['x', 'y', 'z', 'b', 'AxB'], action='append', help='Default: x,y,z')
parser.add_argument('--geom', action='store_true', help='Make a figure containg the geometry from tx to antenna(s)')
parser.add_argument('--ft', action='store_true', help='Add FT strenghts of antenna traces')
args = parser.parse_args()
fname = "ZH_airshower/mysim.sry"
plot_ft_amplitude = args.ft
plot_geometry = args.geom
fig_dir = args.fig_dir
show_plots = args.show_plots
if not args.polarisations:
args.polarisations = ['x','y', 'z']
2022-11-18 19:36:48 +01:00
####
fname_dir = path.dirname(fname)
antennas_fname = path.join(fname_dir, beacon.antennas_fname)
f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
idx = args.ant_idx
if not idx:
if not True:
idx = [0, 1, len(antennas)//2, len(antennas)//2+1, -2, -1]
elif not True:
idx = np.arange(1, 20, 2, dtype=int)
elif True:
# center 6 antennas
names = [55, 56, 57, 65, 66, 45, 46]
idx = [ i for i, ant in enumerate(antennas) if int(ant.name) in names ]
fig1, axs = plt.subplots(1+plot_ft_amplitude*2)
if not plot_ft_amplitude:
axs = [axs]
axs[0].set_xlabel('t [ns]')
axs[0].set_ylabel('[$\mu$V/m]')
if plot_ft_amplitude:
axs[1].set_xlabel('f [GHz]')
axs[1].set_ylabel('Power')
axs[2].set_ylabel("Phase")
axs[2].set_xlabel('f [GHz]')
axs[2].set_ylim(-np.pi,+np.pi)
colorlist = []
for i in idx:
ant = antennas[i]
n_samples = len(ant.t)
samplerate = (ant.t[-1] - ant.t[0])/n_samples
axs[0].axvline(ant.t[0], color='k', alpha=0.5)
mydict = {}
for p in args.polarisations:
pattr = 'E'+str(p)
if p == 'b':
pattr = 'beacon'
elif p == 'AxB':
pattr = 'E_AxB'
mydict[p] = getattr(ant, pattr)
for j, (direction, trace) in enumerate(mydict.items()):
l = axs[0].plot(ant.t, trace, label=f"$E_{{{direction}}}$ {ant.name}", alpha=0.7)
#if False and j == 0 and 't0' in ant.attrs:
# axs[0].axvline(ant.attrs['t0'], color=l[0].get_color(), alpha=0.5)
colorlist.append(l[0].get_color())
if not plot_ft_amplitude:
continue
freqs = ft.fftfreq(n_samples, 1/samplerate)[:n_samples//2]
fft = 2*ft.fft(trace)[:n_samples//2]/n_samples
2022-11-18 19:36:48 +01:00
#axs[1].plot(freqs, np.abs(fft)**2, color=l[0].get_color())
if True:
cft = lib.direct_fourier_transform(f_beacon, ant.t, trace)
amp = 2*len(ant.t) * (cft[0]**2 + cft[1]**2)
#axs[0].axhline(amp, color=l[0].get_color())
print(amp)
phase = np.arctan2(cft[0],cft[1])
2022-11-22 15:37:55 +01:00
axs[1].plot(f_beacon, amp, color=l[0].get_color(), marker='3', alpha=0.8, ms=30)
axs[2].plot(f_beacon, phase, color=l[0].get_color(), marker='3', alpha=0.8, ms=30)
if plot_ft_amplitude:
2022-11-22 15:37:55 +01:00
fig1.legend(loc='center right', ncol=min(2, len(idx)))
else:
axs[0].legend(loc='upper right', ncol=min(3, len(idx)))
if fig_dir:
fig1.savefig(path.join(fig_dir, path.basename(__file__) + f".trace.pdf"))
if plot_geometry:
2022-11-22 15:37:55 +01:00
if len(mydict) == 1:
geom_colorlist = colorlist
else:
# only take the colour belonging to mydict[0]
geom_colorlist = [ colorlist[len(mydict)*(i)] for i in range(len(colorlist)//len(mydict)) ]
fig2, axs2 = plt.subplots(1)
plot_antenna_geometry(antennas, ax=axs2, plot_max_values=False, color='grey', plot_names=False)
2022-11-22 15:37:55 +01:00
plot_antenna_geometry([ antennas[i] for i in idx], ax=axs2, colors=geom_colorlist, plot_max_values=False)
axs2.plot(tx.x, tx.y, marker='X', color='k')
2022-11-18 19:36:48 +01:00
axs2.set_title("Geometry with selected antennas")
if fig_dir:
fig2.savefig(path.join(fig_dir, path.basename(__file__) + f".geom.pdf"))
if show_plots:
plt.show()