Signal Propagation

Propagation module

The modules for the raytracing are stored in the folder SignalProp. All the propagation effects (attenuation, focussing) are also taken account for in the raytracer module itself. The configuration of the raytracing and propagation effects are specified in the config.yaml file under propagation with the following attributes:

module: [string] the ray tracing method to use
ice_model: [string] the description of the refractive index of the ice and all its special effects
attenuation_model: [string] the description of the attenuation of the ice
attenuation_ice: [boolean] whether the attenuation due to the propagation through the ice should be applied. Note: The 1/R amplitude scaling will be applied in either case.
n_freq: [int] the number of frequencies where the attenuation length is calculated for. The remaining frequencies will be determined from a linear interpolation between the reference frequencies. The reference frequencies are equally spaced over the complete frequency range.
focusing: [boolean] whether the focusing effect should be applied.
focusing_limit: [float] the maximum amplification factor of the focusing correction
n_reflections: [int] the maximum number of reflections off a reflective layer at the bottom of the ice layer

Below you find the default settings of the config file.

propagation:
  module: analytic
  ice_model: southpole_2015
  attenuation_model: SP1
  attenuate_ice: True
  n_freq: 25
  focusing: False
  focusing_limit: 2
  n_reflections: 0
  birefringence: False

How to implement new ice-models and information on all the available ice-models and attenuation/birefringence models can be found in the documentation.

Ray tracing

Ray tracing is the module to calculate the trajectory of the emitted radiation. Depending on the ice model one wants to use, the user can specify which ray tracer method should be used by NuRadioMC. This can be done in the config.yaml file by setting the propagation module to the desired module name.

propagation:
  module: [module name]

One has to keep in mind that it is possible some methods only work using specific classes of ice models.

Analytical ray tracer

The analytical ray tracer of NuRadioMC is a method that can only handle the so called simple ice models. Simple ice models are media with a planar geometry and a refractive index with an exponential profile depending on the depth in the ice. The module name to be used in the config file for this method is analytic. The analytic method is implemented in both python but also in c++ for more rapid solving which will automatically be used when all dependencies are available on the users machine. Below you’ll find an overview of the mathematics and an example of the possible rays. The details of this method can be found in the NuRadioMC paper.

Take the following ice model:

\[n(z) = n_{ice} - \Delta_n \exp(\frac{z}{z_0}) [Eq. (1)]\]

where z is the depth and $n_{ice}$, $\Delta_n$, $z_0$ are the parameters of the model. The ray trajectories in a planar medium with an exponential refractive index can be calculated analytically using the optical variational principle - the time it takes for the ray to go from emitter to observer must be a stationary point. Not necessarily a minimum, as it is commonly said. Take, for instance, a ray reflected on a spherical mirror.

According to the variational principle, the ray path in a medium given by Eq. (1) can be written as:

\[y(z) = \pm z_0 (n_{ice}^2 C_0^2 - 1)^{1/2} \cdot \ln\left(\frac{\gamma}{2 (d(\gamma^2 - b\gamma + d)^{1/2} - b\gamma + 2d}\right) + C_1,\]

where y is the horizontal coordinate, $\gamma = \Delta_n \exp(z/z_0)$, $b = 2n_{ice}$, and $d = n_{ice}^2 - C_0^{-2}$. $C_0$ is an integration constant related to the angle at launch position and $C_1$ is another integration constant that gives the starting point.

The ray path can be expressed in closed form, as well as the travel time and the path length. However, the calculation of the frequency-dependent attenuation length must be done numerically.

GSL speed-boost

A C++ version of the ray tracer exists, which can be called from Python, and is significantly faster. The C++ version needs to be compiled on the user’s system, which is attempted automatically if this is not already the case. However, this requires that the GSL library is installed, and the $GSLDIR variable is defined and points to the GNU Scientific Library directory.

export GSLDIR=/path/to/my/GNU_Scientific_Library

Once GSLDIR is configured, the user can also compile it by hand executing the following instruction in the SignalProp/CPPAnalyticRayTracing folder:

python setup.py build_ext --inplace

The C++ analytic ray tracer can also be used standalone, using the MakeFile included. For instructions, please refer to the README included in the NuRadioMC/SignalProp/CPPAnalyticRayTracing directory.

RadioPropa numerical ray tracer (in development)

For ice models other then the simple ones, one need a numerical ray tracer which is provided by the RadioPropa method. This method uses the RadioPropa package which is written in c++. Information on the installation of RadioPropa can found on https://github.com/nu-radio/RadioPropa. The module name for this method is radiopropa.

propagation:
  module: radiopropa

RadioPropa is a modular ray tracing code that solves the eikonal equation for a ray fired at a certain place in a certain direction using a Runge-Kutta method in arbitrary refractivity fields. The implemented NuRadio ray tracer uses this to scan a certain section of the ice in a iterative manner to see whether a channel will be hit or not as shown below

For now, this method can be used for a refractive index with any profile depending on the depth (only z, no x or y dependence) in the ice and some additional features like discontinuities or reflective/transmissive layers. In the future, more effect and the handling of more complex profiles will become available.

Example scripts

How to calculate an analytic ray path

The following code shows how to perform a analytic ray tracing and extract information on the solutions, such as trajectory, travel time, or attenuation.

from NuRadioMC.SignalProp import propagation
from NuRadioMC.SignalProp.analyticraytracing import solution_types, ray_tracing_2D
from NuRadioMC.utilities import medium
from NuRadioReco.utilities import units
import matplotlib.pyplot as plt
import numpy as np

prop = propagation.get_propagation_module('analytic')

ref_index_model = 'greenland_simple'
ice = medium.get_ice_model(ref_index_model)

# Let us work on the y = 0 plane
initial_point = np.array( [70, 0, -300] ) * units.m
final_point = np.array( [100, 0, -30] ) * units.m
attenuation_model = 'GL1'

# This function creates a ray tracing instance refracted index, attenuation model,
# number of frequencies # used for integrating the attenuation and interpolate afterwards,
# and the number of allowed reflections.
rays = prop(ice, attenuation_model,
            n_frequencies_integration=25,
            n_reflections=0)

rays.set_start_and_end_point(initial_point,final_point)
rays.find_solutions()

for i_solution in range(rays.get_number_of_solutions()):

    solution_int = rays.get_solution_type(i_solution)
    solution_type = solution_types[solution_int]

    # To plot the ray path, we can use the 2D ray tracing class, which works on
    # a plane. Since we have been working on the y = 0 plane, we can construct
    # the 2D vectors without translations or rotations. Just ignore the y component.
    rays_2D = ray_tracing_2D(ice, attenuation_model)
    initial_point_2D = np.array( [initial_point[0], initial_point[2]] )
    final_point_2D = np.array( [final_point[0], final_point[2]] )
    C_0 = rays.get_results()[i_solution]['C0']

    xx, zz = rays_2D.get_path(initial_point_2D, final_point_2D, C_0)
    plt.plot(xx, zz, label=solution_type)

    # We can also get the 3D receiving vector at the observer position, for instance
    receive_vector = rays.get_receive_vector(i_solution)
    # Or the path length
    path_length = rays.get_path_length(i_solution)
    # And the travel time
    travel_time = rays.get_travel_time(i_solution)

plt.xlabel('horizontal coordinate [m]')
plt.ylabel('vertical coordinate [m]')
plt.legend()
plt.show()

# We can also calculate the attenuation for a set of frequencies

sampling_rate_detector = 1 * units.GHz
nyquist_frequency = 0.5 * sampling_rate_detector
frequencies = np.linspace(50 * units.MHz, nyquist_frequency, 100)

for i_solution in range(rays.get_number_of_solutions()):

    solution_int = rays.get_solution_type(i_solution)
    solution_type = solution_types[solution_int]

    attenuation = rays.get_attenuation(i_solution, frequencies, nyquist_frequency)

    plt.plot(frequencies/units.MHz, attenuation, label=solution_type)

plt.xlabel('Frequency [MHz]')
plt.ylabel('Attenuation factor')
plt.ylim((0,1))
plt.legend()
plt.show()

How to calculate an radiopropa ray path

The following code shows how to perform a ray tracing and extract information on the solutions, such as trajectory, travel time, or attenuation.

from NuRadioMC.SignalProp import propagation
from NuRadioMC.SignalProp.simple_radiopropa_tracer import solution_types, ray_tracing
from NuRadioMC.utilities import medium
from NuRadioReco.utilities import units
import matplotlib.pyplot as plt
import numpy as np

prop = propagation.get_propagation_module('radiopropa')

ref_index_model = 'greenland_simple'
ice = medium.get_ice_model(ref_index_model)

# Let us work on the y = 0 plane
initial_point = np.array( [70, 0, -300] ) * units.m
final_point = np.array( [100, 0, -30] ) * units.m
attenuation_model = 'GL1'

# This function creates a ray tracing instance refracted index, attenuation model,
# number of frequencies # used for integrating the attenuation and interpolate afterwards,
# and the number of allowed reflections.
rays = prop(ice, attenuation_model,
            n_frequencies_integration=25,
            n_reflections=0)

rays.set_start_and_end_point(initial_point,final_point)
rays.find_solutions()

for i_solution in range(rays.get_number_of_solutions()):

    solution_int = rays.get_solution_type(i_solution)
    solution_type = solution_types[solution_int]

    path = rays.get_path(i_solution)
    # We can calculate the azimuthal angle phi to rotate the
    # 3D path into the 2D plane of the points. This is only
    # necessary if we are not working in the y=0 plane
    launch_vector = rays.get_launch_vector(i_solution))
    phi = np.arctan(launch_vector[1]/launch_vector[0])
    plt.plot(path[:,0]/np.cos(phi), path[:,2], label=solution_type)

    # We can also get the 3D receiving vector at the observer position, for instance
    receive_vector = rays.get_receive_vector(i_solution)
    # Or the path length
    path_length = rays.get_path_length(i_solution)
    # And the travel time
    travel_time = rays.get_travel_time(i_solution)

plt.xlabel('horizontal coordinate [m]')
plt.ylabel('vertical coordinate [m]')
plt.legend()
plt.show()

# We can also calculate the attenuation for a set of frequencies

sampling_rate_detector = 1 * units.GHz
nyquist_frequency = 0.5 * sampling_rate_detector
frequencies = np.linspace(50 * units.MHz, nyquist_frequency, 100)

for i_solution in range(rays.get_number_of_solutions()):

    solution_int = rays.get_solution_type(i_solution)
    solution_type = solution_types[solution_int]

    attenuation = rays.get_attenuation(i_solution, frequencies, nyquist_frequency)

    plt.plot(frequencies/units.MHz, attenuation, label=solution_type)

plt.xlabel('Frequency [MHz]')
plt.ylabel('Attenuation factor')
plt.ylim((0,1))
plt.legend()
plt.show()