Source code for NuRadioReco.utilities.geometryUtilities

import numpy as np
from scipy import constants
from NuRadioReco.utilities import units
from numpy.lib import scimath as SM
import logging
logger = logging.getLogger('NuRadioReco.geometryUtilities')


[docs]def get_time_delay_from_direction(zenith, azimuth, positions, n=1.000293): """ Calculate the time delay between given positions for an arrival direction Parameters ---------- zenith: float [rad] Zenith angle in convention up = 0 azimuth: float [rad] Azimuth angle in convention East = 0, counter-clock-wise positions: array[N x 3] Positions on ground n: float (default: 1.000293) Index of reflection of propagation medium. By default, air is assumed """ shower_axis = np.array([np.sin(zenith) * np.cos(azimuth), np.sin(zenith) * np.sin(azimuth), np.cos(zenith)]) if positions.ndim == 1: times = -(1 / (constants.c / n)) * np.dot(shower_axis, positions) * units.s else: times = np.zeros(positions.shape[0]) for i in range(positions.shape[0]): times[i] = -(1 / (constants.c / n)) * np.dot(shower_axis, positions[i, :]) * units.s return times
[docs]def rot_z(angle): """Angle helper function""" m = np.array( [ [np.cos(angle), -1 * np.sin(angle), 0], [np.sin(angle), np.cos(angle), 0], [0, 0, 1] ] ) return m
[docs]def rot_x(angle): """Angle helper function""" m = np.array( [ [1, 0, 0], [0, np.cos(angle), -1 * np.sin(angle)], [0, np.sin(angle), np.cos(angle)] ] ) return m
[docs]def rot_y(angle): """Angle helper function""" m = np.array( [ [np.cos(angle), 0, np.sin(angle)], [0, 1, 0], [-1 * np.sin(angle), 0, np.cos(angle)] ] ) return m
[docs]def get_efield_in_spherical_coords(efield, theta, phi): """ Get 3D electric field from cartesian coordinates in spherical coordinates, using the arrival directions theta and phi """ ct = np.cos(theta) st = np.sin(theta) cp = np.cos(phi) sp = np.sin(phi) e1 = np.array([st * cp, st * sp, ct]) e2 = np.array([ct * cp, ct * sp, -st]) e3 = np.array([-sp, cp, 0]) # e1 /= linalg.norm(e1) # e2 /= linalg.norm(e2) # e3 /= linalg.norm(e3) transformation_matrix = np.matrix([e1, e2, e3]) # inverse_transformation_matrix = np.linalg.inv(transformation_matrix) efield_2 = np.squeeze(np.asarray(np.dot(transformation_matrix, efield))) return efield_2
[docs]def get_fresnel_angle(zenith_incoming, n_2=1.3, n_1=1.): """ Apply Snell's law for given zenith angle, when a signal travels from n1 to n2 """ t = n_1 / n_2 * np.sin(zenith_incoming) if t > 1: logger.debug('Fresnel refraction results in unphysical values, refraction from {n1} to {n2} with incoming angle {zenith:.1f}, returning None'.format(n1=n_1, n2=n_2, zenith=np.rad2deg(zenith_incoming))) return None else: if(zenith_incoming > 0.5 * np.pi): return np.pi - np.arcsin(t) return np.arcsin(t)
[docs]def get_fresnel_t_p(zenith_incoming, n_2=1.3, n_1=1.): """ returns the coefficient t which is the ratio of the transmitted wave's electric field amplitude to that of the incident wave for parallel polarization (p-wave) this polarization corresponds to the eTheta polarization parallel and perpendicular refers to the signal's polarization with respect to the 'plane of incident' which is defindes as: "the plane of incidence is the plane which contains the surface normal and the propagation vector of the incoming radiation." """ zenith_outgoing = get_fresnel_angle(zenith_incoming, n_2, n_1) if(zenith_outgoing is None): #check for total internal reflection return 0 t = 2 * n_1 * np.cos(zenith_incoming) / (n_1 * np.cos(zenith_outgoing) + n_2 * np.cos(zenith_incoming)) return t
[docs]def get_fresnel_t_s(zenith_incoming, n_2=1.3, n_1=1.): """ returns the coefficient t which is the ratio of the transmitted wave's electric field amplitude to that of the incident wave for perpendicular polarization (s-wave) this polarization corresponds to the ePhi polarization parallel and perpendicular refers to the signal's polarization with respect to the 'plane of incident' which is defindes as: "the plane of incidence is the plane which contains the surface normal and the propagation vector of the incoming radiation." """ zenith_outgoing = get_fresnel_angle(zenith_incoming, n_2, n_1) if(zenith_outgoing is None): #check for total internal reflection return 0 t = 2 * n_1 * np.cos(zenith_incoming) / (n_1 * np.cos(zenith_incoming) + n_2 * np.cos(zenith_outgoing)) return t
[docs]def get_fresnel_r_p(zenith_incoming, n_2=1.3, n_1=1.): """ returns the coefficient r which is the ratio of the reflected wave's electric field amplitude to that of the incident wave for parallel polarization (p-wave) this polarization corresponds to the eTheta polarization parallel and perpendicular refers to the signal's polarization with respect to the 'plane of incident' which is defindes as: "the plane of incidence is the plane which contains the surface normal and the propagation vector of the incoming radiation." """ n = n_2 / n_1 return (-n**2 * np.cos(zenith_incoming) + SM.sqrt(n**2 - np.sin(zenith_incoming)**2)) / \ (n**2 * np.cos(zenith_incoming) + SM.sqrt(n**2 - np.sin(zenith_incoming)**2))
[docs]def get_fresnel_r_s(zenith_incoming, n_2=1.3, n_1=1.): """ returns the coefficient r which is the ratio of the reflected wave's electric field amplitude to that of the incident wave for perpendicular polarization (s-wave) this polarization corresponds to the ePhi polarization parallel and perpendicular refers to the signal's polarization with respect to the 'plane of incident' which is defindes as: "the plane of incidence is the plane which contains the surface normal and the propagation vector of the incoming radiation." """ n = n_2 / n_1 return (np.cos(zenith_incoming) - SM.sqrt(n**2 - np.sin(zenith_incoming)**2)) / \ (np.cos(zenith_incoming) + SM.sqrt(n**2 - np.sin(zenith_incoming)**2))