Starshapes clarity kit

radiotools/atmosphere/cherenkov_radius.py

@@ -0,0 +1,182 @@
+# functions to estimate calculate the cherenkov radius of the radio emission of an air shower
+# code written by Felix Schlueter for the RadioAnalysis framework
+# added to radiotools by Lukas Guelzow
+
+import numpy as np
+import warnings
+
+from radiotools.atmosphere import models as atm
+
+
+def get_cherenkov_radius_model_from_depth(zenith, depth, obs_level, n0, model=None, at=None):
+ """ Calculates the radius of the (Cherenkov) cone with an apex at a given depth along a 
+ shower axis with a given zenith angle. The open angle of the cone equals the 
+ Cherenkov angle for a value of the refractive index at this position.
+
+ Paramter:
+
+ zenith : double
+ Zenith angle (in radian) under which a shower is observed 
+
+ depth : double
+ Slant depth (in g/cm^2), i.e., shower maximum of the observed shower
+
+ obs_level : double
+ Altitude (in meter) of the plane at which the shower is observed 
+
+ n0 : double
+ Refractive index at sea level (!= obs_level)
+
+ model : int
+ Model index for the atmospheric (density) profile model. Needed when no "at" is given
+
+ at : radiotools.atmosphere.models.Atmosphere
+ Atmospheric (density) profile model. Provides the density profile of the atmosphere in the typical 5-layer param.
+
+ Return : cherenkov Radius
+
+ """
+ if at is None:
+ at = atm.Atmosphere(model=model)
+
+ d = at.get_distance_xmax_geometric(zenith, depth, obs_level)
+ return get_cherenkov_radius_model_from_distance(zenith, d, obs_level, n0, at.model)
+
+
+def get_cherenkov_radius_model_from_height(zenith, height, obs_level, n0, model):
+ """ Calculates the radius of the (Cherenkov) cone with an apex at a given height above sea level on a
+ shower axis with a given zenith angle. The open angle of the cone equals the 
+ Cherenkov angle for a value of the refractive index at this position.
+
+ Paramter:
+
+ zenith : double
+ Zenith angle (in radian) under which a shower is observed 
+
+ height : double
+ Height above sea level (in m) of the apex, i.e., shower maximum of the observed shower
+
+ obs_level : double
+ Altitude (in meter) of the plane at which the shower is observed 
+
+ n0 : double
+ Refractive index at sea level (!= obs_level)
+
+ model : int
+ Model index for the atmospheric (density) profile model.
+
+ Return : cherenkov Radius
+
+ """
+
+ angle = get_cherenkov_angle_model(height, n0, model)
+ dmax = atm.get_distance_for_height_above_ground(
+ height - obs_level, zenith, observation_level=obs_level)
+ return cherenkov_radius(angle, dmax)
+
+
+def get_cherenkov_radius_model_from_distance(zenith, d, obs_level, n0, model):
+ """ Calculates the radius of the (Cherenkov) cone with an apex at a given distance from ground 
+ along the shower axis with a given zenith angle. The open angle of the cone equals the 
+ Cherenkov angle for a value of the refractive index at this position.
+
+ Paramter:
+
+ zenith : double
+ Zenith angle (in radian) under which a shower is observed 
+
+ d : double
+ Distance from ground to the apex, i.e., shower maximum of the observed shower along the shower axis (in m)
+
+ obs_level : double
+ Altitude (in meter) of the plane at which the shower is observed 
+
+ n0 : double
+ Refractive index at sea level (!= obs_level)
+
+ model : int
+ Model index for the atmospheric (density) profile model.
+
+ Return : cherenkov Radius
+
+ """
+ height = atm.get_height_above_ground(
+ d, zenith, observation_level=obs_level) + obs_level
+ angle = get_cherenkov_angle_model(height, n0, model)
+ return cherenkov_radius(angle, d)
+
+
+def get_cherenkov_angle_model(height, n0, model):
+ """ Return cherenkov angle for given height above sea level, 
+ refractive index at sea level and atmospheric model. 
+
+ Paramter:
+
+ height : double
+ Height above sea level (in m)
+
+ n0 : double
+ Refractive index at sea level (!= obs_level)
+
+ model : int
+ Model index for the atmospheric (density) profile model.
+
+ Return : cherenkov angle
+
+ """
+ n = atm.get_n(height, n0=n0, model=model)
+ return cherenkov_angle_model(n)
+
+
+def cherenkov_angle_model(n):
+ """ Return cherenkov angle for given refractive index.
+
+ Paramter:
+
+ n : double
+ Refractive index
+
+ Return : cherenkov angle
+
+ """
+ return np.arccos(1 / n)
+
+
+def cherenkov_radius(angle, d):
+ """ Return (cherenkov) radius
+
+ Paramters
+ ---------
+
+ angle : double
+ (Opening) angle of the cone (in rad)
+
+ d : double
+ Heigth of the cone (typically called distance, in meter)
+
+ Returns
+ -------
+
+ radius : double
+ (Cherenkov) radius
+
+ """
+ return np.tan(angle) * d
+
+
+# old: cherenkov_angle_from_density_refractivity(rho, dxmax, n_asl, rho_0, ...)
+# param of the cherenkov angle from star-shape simulations
+# here cherenkov radius refers to radius of strongest emission
+# is used to determine rmax in the RdIdealGrid simulations!
+def cherenkov_angle_param(
+ height, dist, n0, model,
+ a=9.48990456e-01, b=4.48698860e+03,
+ c=1.43097665e+00, d=2.46630811e+06):
+
+ n = atm.get_n(height, n0=n0, model=model)
+ A = a - (b / dist) ** (c) - dist / d
+
+ return A * cherenkov_angle_model(n)
+# # old param
+# def cherenkov_angle_from_density(x, A=0.24905864, B=0.92165234):
+# return np.deg2rad(A * np.log(x) + B)

radiotools/atmosphere/models.py

@@ -115,6 +115,20 @@
  'b': 1e4 * np.array([1111.70, 1128.64, 1413.98, 587.688, 1]),
  'c': 1e-2 * np.array([766099., 641716., 588082., 693300., 5430320300]),
  'h': 1e3 * np.array([7.6, 22.0, 40.4, 100.])
+ },
+ # Lenghu
+ # n @ sea level: 1.000276484920489
+ 40: {'a': 1e4 * np.array([-165.436862, -136.096417, -0.822638028, 0.000573230645, 0.01128292]),
+ 'b': 1e4 * np.array([1197.31182, 1170.37829, 1225.83319, 1331.36795, 1]),
+ 'c': 1e-2 * np.array([1050681.74, 1014055.01, 712892.522, 692298.059, 1.e9]),
+ 'h': 1e3 * np.array([3.703, 9.570, 26.816, 100.])
+ },
+ # Dunhuang
+ # n @ sea level: 1.000273455776266
+ 41: {'a': 1e4 * np.array([-213.042077, -116.247782, 0.00113274359, 0.000571786955, 0.01128292]),
+ 'b': 1e4 * np.array([1258.61938, 1170.17578, 1228.95805, 1228.88998, 1]),
+ 'c': 1e-2 * np.array([1065331.16, 949496.792, 696242.875, 969256.762, 1.e9]),
+ 'h': 1e3 * np.array([3.689, 9.378, 26.299, 100.])
  }
  }
 
@@ -817,7 +831,7 @@ def ftmp(d, zenith, xmax):
  x0 = get_distance_for_height_above_ground(self._get_vertical_height_flat(zenith[i], X[i]), zenith[i])
 
  # finding root e.g., distance for given xmax (when difference is 0)
- dxmax_geo = optimize.brentq(ftmp, -1e3, x0 + 1e4, xtol=1e-6, args=(zenith[i], X[i]))
+ dxmax_geo = optimize.brentq(ftmp, -1e3, x0 + 2e4, xtol=1e-6, args=(zenith[i], X[i]))
 
  height[i] = get_height_above_ground(
  dxmax_geo, zenith[i], observation_level=observation_level) + observation_level
@@ -845,7 +859,7 @@ def ftmp(d, zenith, xmax):
  x0 = get_distance_for_height_above_ground(self._get_vertical_height_flat(zenith[i], X[i]), zenith[i])
 
  # finding root e.g., distance for given xmax (when difference is 0)
- dxmax_geo = optimize.brentq(ftmp, -1e3, x0 + 1e4, xtol=1e-6, args=(zenith[i], X[i]))
+ dxmax_geo = optimize.brentq(ftmp, -1e3, x0 + 2e4, xtol=1e-6, args=(zenith[i], X[i]))
 
  height[i] = get_height_above_ground(
  dxmax_geo, zenith[i], observation_level=observation_level) + observation_level

radiotools/coordinatesystems.py

@@ -3,7 +3,7 @@
 
 import numpy as np
 import math
-from numpy.linalg import linalg
+import numpy.linalg as linalg
 import copy
 import sys
 
@@ -158,7 +158,7 @@ def transform_from_early_late(self, positions, core=None):
 
  _, nY = positions.shape
  if(nY != 3):
- sys.exit("Illeal position given")
+ sys.exit("Illegal position given")
  else:
  result = []
  for pos in positions: