ericteunis/m-thesis-introduction
1
0
Fork 0
mirror of https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git synced 2024-12-21 19:13:32 +01:00
Code Issues Projects Releases Packages Wiki Activity

ZH: dc_grid_power_time_fixes as used for Thesis Figures

Browse source
This commit is contained in:
Eric Teunis de Boone 2023-11-03 17:20:24 +01:00
parent 8e5a8b608d
commit 869ef6dc29
1 changed files with 42 additions and 19 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

61
airshower_beacon_simulation/dc_grid_power_time_fixes.py
View file

@ -18,14 +18,21 @@ import aa_generate_beacon as beacon
import lib import lib
from lib import rit from lib import rit
try:
from tqdm import tqdm
except:
tqdm = lambda x: x
def save_overlapping_traces_figure(test_location, ev, N_plot = 30, wx=200, title_extra=None, fname_distinguish='', fig_dir=None, location_in_plot_text_loc=(0.02, 0.95), location_text=None, **fig_kwargs): def save_overlapping_traces_figure(test_location, ev, N_plot = 30, wx=200, title_extra=None, fname_distinguish='', fig_dir=None, location_in_plot_text_loc=(0.02, 0.95), location_text=None, **fig_kwargs):
P, t_, a_, a_sum, t_sum = rit.pow_and_time(test_location, ev, dt=1) P, t_, a_, a_sum, t_sum = rit.pow_and_time(test_location, ev, dt=1)
fig, axs = plt.subplots(**fig_kwargs) fig, axs = plt.subplots(**fig_kwargs)
axs.set_title("Antenna traces" + (("\n" + title_extra) if title_extra is not None else '') ) #axs.set_title(
#"Antenna traces" + ("\n" if title_extra is not None else '') +
# title_extra )
axs.set_xlabel("Time [ns]") axs.set_xlabel("Time [ns]")
axs.set_ylabel("Amplitude [$\\mu V/m$]") axs.set_ylabel("Amplitude [$\\mu V/m$]")
if location_in_plot_text_loc: if False and location_in_plot_text_loc:
text_loc = location_in_plot_text_loc text_loc = location_in_plot_text_loc
if not location_text: if not location_text:
@ -37,6 +44,9 @@ def save_overlapping_traces_figure(test_location, ev, N_plot = 30, wx=200, title
a_max = [ np.amax(ant.E_AxB) for ant in ev.antennas ] a_max = [ np.amax(ant.E_AxB) for ant in ev.antennas ]
power_sort_idx = np.argsort(a_max) power_sort_idx = np.argsort(a_max)
if N_plot is None:
N_plot = len(a_max)
for i, idx in enumerate(reversed(power_sort_idx)): for i, idx in enumerate(reversed(power_sort_idx)):
if i >= N_plot: if i >= N_plot:
break break
@ -46,7 +56,7 @@ def save_overlapping_traces_figure(test_location, ev, N_plot = 30, wx=200, title
axs.plot(t_[idx], a_[idx], color='r', alpha=alpha, lw=2) axs.plot(t_[idx], a_[idx], color='r', alpha=alpha, lw=2)
axs.set_ylim([None, max(a_max)*1.2]) axs.set_ylim([None, max(a_max)*1.2])
if False: if True:
axs.set_xlim([-28000, -27300]) axs.set_xlim([-28000, -27300])
if fig_dir: if fig_dir:
@ -56,7 +66,7 @@ def save_overlapping_traces_figure(test_location, ev, N_plot = 30, wx=200, title
fig.tight_layout() fig.tight_layout()
fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.{case}.pdf')) fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.{case}.pdf'))
fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.{case}.png'), transparent=True) #fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.{case}.png'), transparent=True)
# Take center between t_low and t_high # Take center between t_low and t_high
if True: if True:
@ -82,7 +92,7 @@ def save_overlapping_traces_figure(test_location, ev, N_plot = 30, wx=200, title
axs.set_xlim(low_xlim, high_xlim) axs.set_xlim(low_xlim, high_xlim)
fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.zoomed.{case}.pdf')) fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.zoomed.{case}.pdf'))
fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.zoomed.{case}.png'), transparent=True) #fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.zoomed.{case}.png'), transparent=True)
return fig return fig
@ -95,9 +105,6 @@ if __name__ == "__main__":
if os.name == 'posix' and "DISPLAY" not in os.environ: if os.name == 'posix' and "DISPLAY" not in os.environ:
matplotlib.use('Agg') matplotlib.use('Agg')
if False: # change overall fontsize
plt.rc('font', size=25)
atm = AtmoCal() atm = AtmoCal()
from scriptlib import MyArgumentParser from scriptlib import MyArgumentParser
@ -122,7 +129,16 @@ if __name__ == "__main__":
if not wanted_cases or 'all' in wanted_cases: if not wanted_cases or 'all' in wanted_cases:
wanted_cases = valid_cases wanted_cases = valid_cases
figsize = (12,8) figsize = (6, 4.5)
figsize_map = (6, 4.8)
if True:
from matplotlib import rcParams
#rcParams["text.usetex"] = True
rcParams["font.family"] = "serif"
rcParams["font.size"] = "15"
rcParams["grid.linestyle"] = 'dotted'
rcParams["figure.figsize"] = figsize
figsize = rcParams['figure.figsize']
fig_dir = args.fig_dir fig_dir = args.fig_dir
show_plots = args.show_plots show_plots = args.show_plots
@ -162,14 +178,15 @@ if __name__ == "__main__":
scale4d = dXref*np.tan(np.deg2rad(4.)) scale4d = dXref*np.tan(np.deg2rad(4.))
scale02d = dXref*np.tan(np.deg2rad(0.2)) scale02d = dXref*np.tan(np.deg2rad(0.2))
Nx, Ny = 21, 21 Nx, Ny = 31, 31
#Nx, Ny = 15, 15
scales = { scales = {
'scale2d': scale2d, 'scale2d': scale2d,
'scale4d': scale4d, 'scale4d': scale4d,
'scale02d': scale02d, 'scale02d': scale02d,
} }
N_plot = 30 N_plot = None
trace_zoom_wx = 100 trace_zoom_wx = 100
plot_titling = { plot_titling = {
@ -180,7 +197,7 @@ if __name__ == "__main__":
} }
power_on_grid_sc_kwargs = dict( power_on_grid_sc_kwargs = dict(
s=250, s=90, # 250 for (10,8)
cmap='inferno' cmap='inferno'
) )
@ -293,6 +310,9 @@ if __name__ == "__main__":
xx, yy, p, maxp_loc = rit.shower_plane_slice(ev, X=X, Nx=1, Ny=1, wx=0, wy=0, zgr=zgr) xx, yy, p, maxp_loc = rit.shower_plane_slice(ev, X=X, Nx=1, Ny=1, wx=0, wy=0, zgr=zgr)
no_offset_maximum_power = p no_offset_maximum_power = p
if True and np.isclose(X, 400):
no_offset_maximum_power = 200
# we forced no_offset for determing the reference vmax # we forced no_offset for determing the reference vmax
# simply continue with the rest # simply continue with the rest
@ -313,28 +333,31 @@ if __name__ == "__main__":
loc = (0)*ev.uAxB + (0)*ev.uAxAxB + dX * ev.uA loc = (0)*ev.uAxB + (0)*ev.uAxAxB + dX * ev.uA
loc[0:1] = 0 loc[0:1] = 0
fig = save_overlapping_traces_figure(loc, ev, N_plot=N_plot, wx=trace_zoom_wx, title_extra = plot_titling[case], fname_distinguish=f'{case}.axis', location_text=f"on simulation axis, $X={X}$", fig_dir=fig_dir, figsize=figsize) fig = save_overlapping_traces_figure(loc, ev, N_plot=N_plot, wx=trace_zoom_wx, title_extra = plot_titling[case], fname_distinguish=f'{case}.axis.X{X}', location_text=f"on simulation axis, $X={X}$", fig_dir=fig_dir, figsize=figsize)
plt.close(fig) plt.close(fig)
if not True:
continue;
# #
# Measure power on grid # Measure power on grid
# and plot overlapping traces at position with highest power # and plot overlapping traces at position with highest power
# #
for scalename, scale in scales.items(): for scalename, scale in tqdm(scales.items()):
wx, wy = scale, scale wx, wy = scale, scale
print(f"Starting grid measurement for figure {case} with {scalename}") print(f"Starting grid measurement for figure {case} with {scalename}")
xx, yy, p, maxp_loc = rit.shower_plane_slice(ev, X=X, Nx=Nx, Ny=Nx, wx=wx, wy=wy, zgr=zgr) xx, yy, p, maxp_loc = rit.shower_plane_slice(ev, X=X, Nx=Nx, Ny=Nx, wx=wx, wy=wy, zgr=zgr)
fig, axs = rit.slice_figure(ev, X, xx, yy, p, mode='sp', scatter_kwargs={**dict( fig, axs = rit.slice_figure(ev, X, xx, yy, p, mode='sp', scatter_kwargs={**dict(
vmax=no_offset_maximum_power, vmax=no_offset_maximum_power,
vmin=0, vmin=0,
), **power_on_grid_sc_kwargs} ) ), **power_on_grid_sc_kwargs},
figsize=figsize_map)
suptitle = fig._suptitle.get_text() suptitle = fig._suptitle.get_text()
fig.suptitle("") fig.suptitle("")
axs.set_title("Shower plane slice\n" + plot_titling[case] + "\n" + suptitle) #axs.set_title("Shower plane slice\n" + plot_titling[case] + "\n" + suptitle)
axs.set_aspect('equal', 'datalim') axs.set_aspect('equal', 'datalim')
axs.legend(title=snr_str) #axs.legend(title=snr_str,loc='upper right')
axs.set_xlim(1.1*min(xx)/1e3, 1.1*max(xx)/1e3) axs.set_xlim(1.1*min(xx)/1e3, 1.1*max(xx)/1e3)
axs.set_ylim(1.1*min(yy)/1e3, 1.1*max(yy)/1e3) axs.set_ylim(1.1*min(yy)/1e3, 1.1*max(yy)/1e3)
@ -347,7 +370,7 @@ if __name__ == "__main__":
# #
# Plot overlapping traces at highest power of each scale # Plot overlapping traces at highest power of each scale
# #
fig = save_overlapping_traces_figure(maxp_loc, ev, N_plot=N_plot, wx=trace_zoom_wx, title_extra = plot_titling[case] + ', ' + scalename + ' best', fname_distinguish=scalename+'.best', fig_dir=fig_dir, figsize=figsize) fig = save_overlapping_traces_figure(maxp_loc, ev, N_plot=N_plot, wx=trace_zoom_wx, title_extra = plot_titling[case] + ', ' + scalename + ' best', fname_distinguish=scalename+f'.best.X{X}', fig_dir=fig_dir, figsize=figsize)
# #
# and plot overlapping traces at two other locations # and plot overlapping traces at two other locations
@ -357,7 +380,7 @@ if __name__ == "__main__":
# only add distance horizontally # only add distance horizontally
location = maxp_loc + np.sqrt(dist*1e3)*np.array([1,1,0]) location = maxp_loc + np.sqrt(dist*1e3)*np.array([1,1,0])
fig = save_overlapping_traces_figure(location, ev, N_plot=N_plot, wx=wx, title_extra = plot_titling[case] + ', ' + scalename + f', x + {dist}km', fname_distinguish=f'{scalename}.x{dist}', fig_dir=fig_dir, figsize=figsize) fig = save_overlapping_traces_figure(location, ev, N_plot=N_plot, wx=wx, title_extra = plot_titling[case] + ', ' + scalename + f', x + {dist}km', fname_distinguish=f'{scalename}.x{dist}.X{X}', fig_dir=fig_dir, figsize=figsize)
plt.close(fig) plt.close(fig)