"""
This module contains utility functions for geometry calculations.
Inparticular to calculate time delays between positions for a given arrival
direction and to simulate reflection and refraction at boundaries between different media.
"""
from NuRadioReco.utilities import units, ice
from numpy.lib import scimath as SM
from scipy import constants
import numpy as np
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
Returns
-------
times: np.array of floats
Time delays
"""
shower_axis = np.array([np.sin(zenith) * np.cos(azimuth), np.sin(zenith) * np.sin(azimuth), np.cos(zenith)])
times = -(1 / (constants.c / n)) * np.dot(positions, shower_axis) * 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
Parameters
----------
efield: np.array
Electric field in cartesian coordinates
theta: float
Zenith angle of the arriving signal
phi: float
Azimuth angle of the arriving signal
Returns
-------
np.array
Electric field in spherical coordinates
"""
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])
transformation_matrix = np.matrix([e1, e2, e3])
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
Parameters
----------
zenith_incoming: float
Zenith angle of the incoming signal
n_2: float
Refractive index of the medium the signal is transmitted into
n_1: float
Refractive index of the medium the signal is coming from
Returns
-------
float
Fresnel angle
"""
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 defined as: "the plane of incidence
is the plane which contains the surface normal and the propagation vector
of the incoming radiation."
Parameters
----------
zenith_incoming: float
Zenith angle of the incoming signal
n_2: float
Refractive index of the medium the signal is transmitted into
n_1: float
Refractive index of the medium the signal is coming from
Returns
-------
float
Fresnel coefficient t for theta (parallel) polarization
"""
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 defined as: "the plane of incidence
is the plane which contains the surface normal and the propagation vector
of the incoming radiation."
Parameters
----------
zenith_incoming: float
Zenith angle of the incoming signal
n_2: float
Refractive index of the medium the signal is transmitted into
n_1: float
Refractive index of the medium the signal is coming from
Returns
-------
float
Fresnel coefficient t for phi (perpendicular) polarization
"""
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 defined as: "the plane of incidence
is the plane which contains the surface normal and the propagation vector
of the incoming radiation."
Parameters
----------
zenith_incoming: float
Zenith angle of the incoming signal
n_2: float
Refractive index of the medium the signal is reflected from
n_1: float
Refractive index of the medium the signal is coming from
Returns
-------
float
Fresnel coefficient r for theta (parallel) polarization
"""
n = n_2 / n_1
return np.conjugate((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 defined as: "the plane of incidence
is the plane which contains the surface normal and the propagation vector
of the incoming radiation."
Parameters
----------
zenith_incoming: float
Zenith angle of the incoming signal
n_2: float
Refractive index of the medium the signal is reflected from
n_1: float
Refractive index of the medium the signal is coming from
Returns
-------
float
Fresnel coefficient r for phi (perpendicular) polarization
"""
n = n_2 / n_1
return np.conjugate((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)))
[docs]
def fresnel_factors_and_signal_zenith(detector, station, channel_id, zenith):
"""
Returns the zenith angle at the antenna and the fresnel coefficients t for theta (parallel)
and phi (perpendicular) polarization. Handles potential refraction into the firn if that
applies to the antenna position.
WARNING: for deeper channels this function might be inacccurate. Consider using raytracing.
parallel and perpendicular refers to the signal's polarization with respect
to the 'plane of incident' which is defined as: "the plane of incidence
is the plane which contains the surface normal and the propagation vector
of the incoming radiation.
Parameters
----------
detector: NuRadioReco.detector.Detector
Detector object
station: NuRadioReco.detector.Station
Station object
channel_id: int
Channel ID of the desired channel
zenith: float
Zenith angle of the incoming signal
Returns
-------
zenith_antenna: float
Zenith angle at the antenna (potentially including refraction)
t_theta: float
Fresnel transmission coefficient for theta polarization
t_phi: float
Fresnel transmission coefficient for phi polarization
"""
n_ice = ice.get_refractive_index(-0.01, detector.get_site(station.get_id()))
# no reflection/refraction at the ice-air boundary
zenith_antenna = zenith
t_theta = 1.
t_phi = 1.
position = detector.get_relative_position(station.get_id(), channel_id)
if position[2] < -3 * units.m:
logger.warning("This function might return inaccurate results for deep in-ice antennas. Consider using raytracing instead.")
# first check case if signal comes from above
if (zenith <= 0.5 * np.pi) and station.is_cosmic_ray() and (position[2] <= 0):
# is antenna below surface?
zenith_antenna = get_fresnel_angle(zenith, n_ice, 1)
t_theta = get_fresnel_t_p(zenith, n_ice, 1)
t_phi = get_fresnel_t_s(zenith, n_ice, 1)
logger.debug(("Channel {:d}: electric field is refracted into the firn. "
"theta {:.0f} -> {:.0f}. Transmission coefficient p (eTheta) "
"{:.2f} s (ePhi) {:.2f}".format(
channel_id, zenith / units.deg, zenith_antenna / units.deg, t_theta, t_phi)))
elif position[2] > 0:
# now the signal is coming from below, do we have an antenna above the surface?
zenith_antenna = get_fresnel_angle(zenith, 1., n_ice)
if zenith_antenna is None:
logger.warning("Fresnel reflection at air-firn boundary leads to unphysical results, no reconstruction possible")
return None, None, None
return zenith_antenna, t_theta, t_phi