From fb89490be8c6e08e7756ecf9499d40bc0e732f18 Mon Sep 17 00:00:00 2001 From: afedynitch Date: Mon, 27 Jun 2022 15:06:47 +0800 Subject: [PATCH 1/2] Fix: For ICatm detector zenith for cache --- MCEq/geometry/density_profiles.py | 722 +++++++++++++++++------------- 1 file changed, 402 insertions(+), 320 deletions(-) diff --git a/MCEq/geometry/density_profiles.py b/MCEq/geometry/density_profiles.py index 7e6910f..5530ae8 100644 --- a/MCEq/geometry/density_profiles.py +++ b/MCEq/geometry/density_profiles.py @@ -1,4 +1,3 @@ - from abc import ABCMeta, abstractmethod from six import with_metaclass from os.path import join @@ -29,11 +28,12 @@ class EarthsAtmosphere(with_metaclass(ABCMeta)): def __init__(self, *args, **kwargs): from MCEq.geometry.geometry import EarthGeometry - self.geom = kwargs.pop('geometry', EarthGeometry()) + + self.geom = kwargs.pop("geometry", EarthGeometry()) self.thrad = None self.theta_deg = None self._max_den = config.max_density - self.max_theta = 90. + self.max_theta = 90.0 self.location = None self.season = None @@ -67,16 +67,20 @@ def calculate_density_spline(self, n_steps=2000): from scipy.interpolate import UnivariateSpline if self.theta_deg is None: - raise Exception('zenith angle not set') + raise Exception("zenith angle not set") else: info( - 5, 'Calculating spline of rho(X) for zenith {0:4.1f} degrees.'. - format(self.theta_deg)) + 5, + "Calculating spline of rho(X) for zenith {0:4.1f} degrees.".format( + self.theta_deg + ), + ) thrad = self.thrad path_length = self.geom.l(thrad) vec_rho_l = np.vectorize( - lambda delta_l: self.get_density(self.geom.h(delta_l, thrad))) + lambda delta_l: self.get_density(self.geom.h(delta_l, thrad)) + ) dl_vec = np.linspace(0, path_length, n_steps) now = time() @@ -85,7 +89,7 @@ def calculate_density_spline(self, n_steps=2000): X_int = cumtrapz(vec_rho_l(dl_vec), dl_vec) # dl_vec = dl_vec[1:] - info(5, '.. took {0:1.2f}s'.format(time() - now)) + info(5, ".. took {0:1.2f}s".format(time() - now)) # Save depth value at h_obs self._max_X = X_int[-1] @@ -97,45 +101,44 @@ def calculate_density_spline(self, n_steps=2000): self._s_h2X = UnivariateSpline(h_intp, np.log(X_intp), k=2, s=0.0) self._s_X2rho = UnivariateSpline(X_int, vec_rho_l(dl_vec), k=2, s=0.0) - self._s_lX2h = UnivariateSpline(np.log(X_intp)[::-1], - h_intp[::-1], - k=2, - s=0.0) + self._s_lX2h = UnivariateSpline(np.log(X_intp)[::-1], h_intp[::-1], k=2, s=0.0) @property def max_X(self): """Depth at altitude 0.""" - if not hasattr(self, '_max_X'): + if not hasattr(self, "_max_X"): self.set_theta(0) return self._max_X @property def max_den(self): """Density at altitude 0.""" - if not hasattr(self, '_max_den'): + if not hasattr(self, "_max_den"): self.set_theta(0) return self._max_den @property def s_h2X(self): """Spline for conversion from altitude to depth.""" - if not hasattr(self, '_s_h2X'): + if not hasattr(self, "_s_h2X"): self.set_theta(0) return self._s_h2X + @property def s_X2rho(self): """Spline for conversion from depth to density.""" - if not hasattr(self, '_s_X2rho'): + if not hasattr(self, "_s_X2rho"): self.set_theta(0) return self._s_X2rho + @property def s_lX2h(self): """Spline for conversion from depth to altitude.""" - if not hasattr(self, '_s_lX2h'): + if not hasattr(self, "_s_lX2h"): self.set_theta(0) return self._s_lX2h - def set_theta(self, theta_deg, force_spline_calc=False): + def set_theta(self, theta_deg): """Configures geometry and initiates spline calculation for :math:`\\rho(X)`. @@ -148,11 +151,9 @@ def set_theta(self, theta_deg, force_spline_calc=False): Args: theta_deg (float): zenith angle :math:`\\theta` at detector - force_spline_calc (bool): forces (re-)calculation of the - spline for each call """ - if theta_deg < 0. or theta_deg > self.max_theta: - raise Exception('Zenith angle not in allowed range.') + if theta_deg < 0.0 or theta_deg > self.max_theta: + raise Exception("Zenith angle not in allowed range.") self.thrad = theta_rad(theta_deg) self.theta_deg = theta_deg @@ -171,7 +172,7 @@ def r_X2rho(self, X): float: :math:`1/\\rho` in cm**3/g """ - return 1. / self.s_X2rho(X) + return 1.0 / self.s_X2rho(X) def h2X(self, h): """Returns the depth along path as function of height above @@ -221,9 +222,9 @@ def X2rho(self, X): return self.s_X2rho(X) def moliere_air(self, h_cm): - """Returns the Moliere unit of air for US standard atmosphere. """ + """Returns the Moliere unit of air for US standard atmosphere.""" - return 9.3 / (self.get_density(h_cm) * 100.) + return 9.3 / (self.get_density(h_cm) * 100.0) def nref_rel_air(self, h_cm): """Returns the refractive index - 1 in air (density parametrization @@ -233,17 +234,15 @@ def nref_rel_air(self, h_cm): return 0.000283 * self.get_density(h_cm) / self.get_density(0) def gamma_cherenkov_air(self, h_cm): - """Returns the Lorentz factor gamma of Cherenkov threshold in air (MeV). - """ + """Returns the Lorentz factor gamma of Cherenkov threshold in air (MeV).""" nrel = self.nref_rel_air(h_cm) - return (1. + nrel) / np.sqrt(2. * nrel + nrel**2) + return (1.0 + nrel) / np.sqrt(2.0 * nrel + nrel**2) def theta_cherenkov_air(self, h_cm): - """Returns the Cherenkov angle in air (degrees). - """ + """Returns the Cherenkov angle in air (degrees).""" - return np.arccos(1. / (1. + self.nref_rel_air(h_cm))) * 180. / np.pi + return np.arccos(1.0 / (1.0 + self.nref_rel_air(h_cm))) * 180.0 / np.pi class CorsikaAtmosphere(EarthsAtmosphere): @@ -264,6 +263,7 @@ class CorsikaAtmosphere(EarthsAtmosphere): location (str): see :func:`init_parameters` season (str,optional): see :func:`init_parameters` """ + _atm_param = None def __init__(self, location, season=None): @@ -271,21 +271,24 @@ def __init__(self, location, season=None): ("USStd", None), ("BK_USStd", None), ("Karlsruhe", None), - ("ANTARES/KM3NeT-ORCA", 'Summer'), - ("ANTARES/KM3NeT-ORCA", 'Winter'), - ("KM3NeT-ARCA", 'Summer'), - ("KM3NeT-ARCA", 'Winter'), - ("KM3NeT",None), - ('SouthPole','December'), - ('PL_SouthPole','January'), - ('PL_SouthPole','August'), + ("ANTARES/KM3NeT-ORCA", "Summer"), + ("ANTARES/KM3NeT-ORCA", "Winter"), + ("KM3NeT-ARCA", "Summer"), + ("KM3NeT-ARCA", "Winter"), + ("KM3NeT", None), + ("SouthPole", "December"), + ("PL_SouthPole", "January"), + ("PL_SouthPole", "August"), ] - assert (location, season) in cka_atmospheres, \ - '{0}/{1} not available for CorsikaAtmsophere'.format( - location, season - ) + assert ( + location, + season, + ) in cka_atmospheres, "{0}/{1} not available for CorsikaAtmsophere".format( + location, season + ) self.init_parameters(location, season) import MCEq.geometry.corsikaatm.corsikaatm as corsika_acc + self.corsika_acc = corsika_acc EarthsAtmosphere.__init__(self) @@ -324,108 +327,147 @@ def init_parameters(self, location, season): if location == "USStd": _aatm = np.array([-186.5562, -94.919, 0.61289, 0.0, 0.01128292]) _batm = np.array([1222.6562, 1144.9069, 1305.5948, 540.1778, 1.0]) - _catm = np.array([994186.38, 878153.55, 636143.04, 772170., 1.0e9]) + _catm = np.array([994186.38, 878153.55, 636143.04, 772170.0, 1.0e9]) _thickl = np.array( - [1036.102549, 631.100309, 271.700230, 3.039494, 0.001280]) + [1036.102549, 631.100309, 271.700230, 3.039494, 0.001280] + ) _hlay = np.array([0, 4.0e5, 1.0e6, 4.0e6, 1.0e7]) elif location == "BK_USStd": - _aatm = np.array([ - -149.801663, -57.932486, 0.63631894, 4.3545369e-4, 0.01128292 - ]) + _aatm = np.array( + [-149.801663, -57.932486, 0.63631894, 4.3545369e-4, 0.01128292] + ) _batm = np.array([1183.6071, 1143.0425, 1322.9748, 655.69307, 1.0]) - _catm = np.array( - [954248.34, 800005.34, 629568.93, 737521.77, 1.0e9]) + _catm = np.array([954248.34, 800005.34, 629568.93, 737521.77, 1.0e9]) _thickl = np.array( - [1033.804941, 418.557770, 216.981635, 4.344861, 0.001280]) + [1033.804941, 418.557770, 216.981635, 4.344861, 0.001280] + ) _hlay = np.array([0.0, 7.0e5, 1.14e6, 3.7e6, 1.0e7]) elif location == "Karlsruhe": - _aatm = np.array( - [-118.1277, -154.258, 0.4191499, 5.4094056e-4, 0.01128292]) + _aatm = np.array([-118.1277, -154.258, 0.4191499, 5.4094056e-4, 0.01128292]) _batm = np.array([1173.9861, 1205.7625, 1386.7807, 555.8935, 1.0]) - _catm = np.array([919546., 963267.92, 614315., 739059.6, 1.0e9]) + _catm = np.array([919546.0, 963267.92, 614315.0, 739059.6, 1.0e9]) _thickl = np.array( - [1055.858707, 641.755364, 272.720974, 2.480633, 0.001280]) + [1055.858707, 641.755364, 272.720974, 2.480633, 0.001280] + ) _hlay = np.array([0.0, 4.0e5, 1.0e6, 4.0e6, 1.0e7]) - elif location == "KM3NeT": # averaged over detector and season - _aatm = np.array([-141.31449999999998, -8.256029999999999, 0.6132505, -0.025998975, 0.4024275]) - _batm = np.array([1153.0349999999999, 1263.3325, 1257.0724999999998, 404.85974999999996, 1.0]) + elif location == "KM3NeT": # averaged over detector and season + _aatm = np.array( + [ + -141.31449999999998, + -8.256029999999999, + 0.6132505, + -0.025998975, + 0.4024275, + ] + ) + _batm = np.array( + [ + 1153.0349999999999, + 1263.3325, + 1257.0724999999998, + 404.85974999999996, + 1.0, + ] + ) _catm = np.array([967990.75, 668591.75, 636790.0, 814070.75, 21426175.0]) - _thickl = np.array([1011.8521512499999, 275.84507575000003, 51.0230705, 2.983134, 0.21927724999999998]) + _thickl = np.array( + [ + 1011.8521512499999, + 275.84507575000003, + 51.0230705, + 2.983134, + 0.21927724999999998, + ] + ) _hlay = np.array([0.0, 993750.0, 2081250.0, 4150000.0, 6877500.0]) elif location == "ANTARES/KM3NeT-ORCA": - if season == 'Summer': + if season == "Summer": _aatm = np.array([-158.85, -5.38682, 0.889893, -0.0286665, 0.50035]) _batm = np.array([1145.62, 1176.79, 1248.92, 415.543, 1.0]) _catm = np.array([998469.0, 677398.0, 636790.0, 823489.0, 16090500.0]) - _thickl = np.array([986.951713, 306.4668, 40.546793, 4.288721, 0.277182]) + _thickl = np.array( + [986.951713, 306.4668, 40.546793, 4.288721, 0.277182] + ) _hlay = np.array([0, 9.0e5, 22.0e5, 38.0e5, 68.2e5]) - elif season == 'Winter': + elif season == "Winter": _aatm = np.array([-132.16, -2.4787, 0.298031, -0.0220264, 0.348021]) _batm = np.array([1120.45, 1203.97, 1163.28, 360.027, 1.0]) _catm = np.array([933697.0, 643957.0, 636790.0, 804486.0, 23109000.0]) - _thickl = np.array([988.431172, 273.033464, 37.185105, 1.162987, 0.192998]) + _thickl = np.array( + [988.431172, 273.033464, 37.185105, 1.162987, 0.192998] + ) _hlay = np.array([0, 9.5e5, 22.0e5, 47.0e5, 68.2e5]) elif location == "KM3NeT-ARCA": - if season == 'Summer': + if season == "Summer": _aatm = np.array([-157.857, -28.7524, 0.790275, -0.0286999, 0.481114]) _batm = np.array([1190.44, 1171.0, 1344.78, 445.357, 1.0]) _catm = np.array([1006100.0, 758614.0, 636790.0, 817384.0, 16886800.0]) - _thickl = np.array([1032.679434, 328.978681, 80.601135, 4.420745, 0.264112]) + _thickl = np.array( + [1032.679434, 328.978681, 80.601135, 4.420745, 0.264112] + ) _hlay = np.array([0, 9.0e5, 18.0e5, 38.0e5, 68.2e5]) - elif season == 'Winter': + elif season == "Winter": _aatm = np.array([-116.391, 3.5938, 0.474803, -0.0246031, 0.280225]) _batm = np.array([1155.63, 1501.57, 1271.31, 398.512, 1.0]) _catm = np.array([933697.0, 594398.0, 636790.0, 810924.0, 29618400.0]) - _thickl = np.array([1039.346286, 194.901358, 45.759249, 2.060083, 0.142817]) + _thickl = np.array( + [1039.346286, 194.901358, 45.759249, 2.060083, 0.142817] + ) _hlay = np.array([0, 12.25e5, 21.25e5, 43.0e5, 70.5e5]) - elif location == 'SouthPole': - if season == 'December': - _aatm = np.array( - [-128.601, -39.5548, 1.13088, -0.00264960, 0.00192534]) + elif location == "SouthPole": + if season == "December": + _aatm = np.array([-128.601, -39.5548, 1.13088, -0.00264960, 0.00192534]) _batm = np.array([1139.99, 1073.82, 1052.96, 492.503, 1.0]) - _catm = np.array( - [861913., 744955., 675928., 829627., 5.8587010e9]) + _catm = np.array([861913.0, 744955.0, 675928.0, 829627.0, 5.8587010e9]) _thickl = np.array( - [1011.398804, 588.128367, 240.955360, 3.964546, 0.000218]) + [1011.398804, 588.128367, 240.955360, 3.964546, 0.000218] + ) _hlay = np.array([0.0, 4.0e5, 1.0e6, 4.0e6, 1.0e7]) elif season == "June": _aatm = np.array( - [-163.331, -65.3713, 0.402903, -0.000479198, 0.00188667]) + [-163.331, -65.3713, 0.402903, -0.000479198, 0.00188667] + ) _batm = np.array([1183.70, 1108.06, 1424.02, 207.595, 1.0]) - _catm = np.array( - [875221., 753213., 545846., 793043., 5.9787908e9]) + _catm = np.array([875221.0, 753213.0, 545846.0, 793043.0, 5.9787908e9]) _thickl = np.array( - [1020.370363, 586.143464, 228.374393, 1.338258, 0.000214]) + [1020.370363, 586.143464, 228.374393, 1.338258, 0.000214] + ) _hlay = np.array([0.0, 4.0e5, 1.0e6, 4.0e6, 1.0e7]) else: - raise Exception('CorsikaAtmosphere(): Season "' + season + - '" not parameterized for location SouthPole.') - elif location == 'PL_SouthPole': - if season == 'January': - _aatm = np.array( - [-113.139, -7930635, -54.3888, -0.0, 0.00421033]) + raise Exception( + 'CorsikaAtmosphere(): Season "' + + season + + '" not parameterized for location SouthPole.' + ) + elif location == "PL_SouthPole": + if season == "January": + _aatm = np.array([-113.139, -7930635, -54.3888, -0.0, 0.00421033]) _batm = np.array([1133.10, 1101.20, 1085.00, 1098.00, 1.0]) - _catm = np.array( - [861730., 826340., 790950., 682800., 2.6798156e9]) - _thickl = np.array([ - 1019.966898, 718.071682, 498.659703, 340.222344, 0.000478 - ]) + _catm = np.array([861730.0, 826340.0, 790950.0, 682800.0, 2.6798156e9]) + _thickl = np.array( + [1019.966898, 718.071682, 498.659703, 340.222344, 0.000478] + ) _hlay = np.array([0.0, 2.67e5, 5.33e5, 8.0e5, 1.0e7]) elif season == "August": - _aatm = np.array( - [-59.0293, -21.5794, -7.14839, 0.0, 0.000190175]) + _aatm = np.array([-59.0293, -21.5794, -7.14839, 0.0, 0.000190175]) _batm = np.array([1079.0, 1071.90, 1182.0, 1647.1, 1.0]) - _catm = np.array( - [764170., 699910., 635650., 551010., 59.329575e9]) + _catm = np.array([764170.0, 699910.0, 635650.0, 551010.0, 59.329575e9]) _thickl = np.array( - [1019.946057, 391.739652, 138.023515, 43.687992, 0.000022]) + [1019.946057, 391.739652, 138.023515, 43.687992, 0.000022] + ) _hlay = np.array([0.0, 6.67e5, 13.33e5, 2.0e6, 1.0e7]) else: - raise Exception('CorsikaAtmosphere(): Season "' + season + - '" not parameterized for location SouthPole.') + raise Exception( + 'CorsikaAtmosphere(): Season "' + + season + + '" not parameterized for location SouthPole.' + ) else: - raise Exception("CorsikaAtmosphere:init_parameters(): Location " + - str(location) + " not parameterized.") + raise Exception( + "CorsikaAtmosphere:init_parameters(): Location " + + str(location) + + " not parameterized." + ) self._atm_param = np.array([_aatm, _batm, _catm, _thickl, _hlay]) @@ -458,7 +500,7 @@ def depth2height(self, x_v): return height def get_density(self, h_cm): - """ Returns the density of air in g/cm**3. + """Returns the density of air in g/cm**3. Uses the optimized module function :func:`corsika_get_density_jit`. @@ -472,7 +514,7 @@ def get_density(self, h_cm): # return corsika_get_density_jit(h_cm, self._atm_param) def get_mass_overburden(self, h_cm): - """ Returns the mass overburden in atmosphere in g/cm**2. + """Returns the mass overburden in atmosphere in g/cm**2. Uses the optimized module function :func:`corsika_get_m_overburden_jit` @@ -517,15 +559,16 @@ def calc_thickl(self): """ from scipy.integrate import quad + thickl = [] for h in self._atm_param[4]: - thickl.append('{0:4.6f}'.format( - quad(self.get_density, h, 112.8e5, epsrel=1e-4)[0])) - info(5, '_thickl = np.array([' + ', '.join(thickl) + '])') + thickl.append( + "{0:4.6f}".format(quad(self.get_density, h, 112.8e5, epsrel=1e-4)[0]) + ) + info(5, "_thickl = np.array([" + ", ".join(thickl) + "])") return thickl - class IsothermalAtmosphere(EarthsAtmosphere): """Isothermal model of the atmosphere. @@ -545,7 +588,7 @@ class IsothermalAtmosphere(EarthsAtmosphere): X0 (float): Ground level overburden """ - def __init__(self, location, season, hiso_km=6.3, X0=1300.): + def __init__(self, location, season, hiso_km=6.3, X0=1300.0): self.hiso_cm = hiso_km * 1e5 self.X0 = X0 self.location = location @@ -554,7 +597,7 @@ def __init__(self, location, season, hiso_km=6.3, X0=1300.): EarthsAtmosphere.__init__(self) def get_density(self, h_cm): - """ Returns the density of air in g/cm**3. + """Returns the density of air in g/cm**3. Args: h_cm (float): height in cm @@ -565,7 +608,7 @@ def get_density(self, h_cm): return self.X0 / self.hiso_cm * np.exp(-h_cm / self.hiso_cm) def get_mass_overburden(self, h_cm): - """ Returns the mass overburden in atmosphere in g/cm**2. + """Returns the mass overburden in atmosphere in g/cm**2. Args: h_cm (float): height in cm @@ -575,6 +618,7 @@ def get_mass_overburden(self, h_cm): """ return self.X0 * np.exp(-h_cm / self.hiso_cm) + class MSIS00Atmosphere(EarthsAtmosphere): """Wrapper class for a python interface to the NRLMSISE-00 model. @@ -592,31 +636,26 @@ class MSIS00Atmosphere(EarthsAtmosphere): season (str,optional): see :func:`init_parameters` """ - def __init__(self, - location, - season=None, - doy=None, - use_loc_altitudes=False): + def __init__(self, location, season=None, doy=None, use_loc_altitudes=False): from MCEq.geometry.nrlmsise00_mceq import cNRLMSISE00 msis_atmospheres = [ - 'SouthPole', - 'Karlsruhe', - 'Geneva', - 'Tokyo', - 'SanGrasso', - 'TelAviv', - 'KSC', - 'SoudanMine', - 'Tsukuba', - 'LynnLake', - 'PeaceRiver', - 'FtSumner' + "SouthPole", + "Karlsruhe", + "Geneva", + "Tokyo", + "SanGrasso", + "TelAviv", + "KSC", + "SoudanMine", + "Tsukuba", + "LynnLake", + "PeaceRiver", + "FtSumner", ] - assert location in msis_atmospheres, \ - '{0} not available for MSIS00Atmosphere'.format( - location - ) + assert ( + location in msis_atmospheres + ), "{0} not available for MSIS00Atmosphere".format(location) self._msis = cNRLMSISE00() @@ -644,11 +683,11 @@ def init_parameters(self, location, season, doy, use_loc_altitudes): # Clear cached value to force spline recalculation self.theta_deg = None if use_loc_altitudes: - info(0, 'Using loc altitude', self._msis.alt_surface, 'cm') + info(0, "Using loc altitude", self._msis.alt_surface, "cm") self.geom.h_obs = self._msis.alt_surface def get_density(self, h_cm): - """ Returns the density of air in g/cm**3. + """Returns the density of air in g/cm**3. Wraps around ctypes calls to the NRLMSISE-00 C library. @@ -661,7 +700,7 @@ def get_density(self, h_cm): return self._msis.get_density(h_cm) def set_location(self, location): - """ Changes MSIS location by strings defined in _msis_wrapper. + """Changes MSIS location by strings defined in _msis_wrapper. Args: location (str): location as defined in :class:`NRLMSISE-00.` @@ -670,7 +709,7 @@ def set_location(self, location): self._msis.set_location(location) def set_season(self, month): - """ Changes MSIS location by month strings defined in _msis_wrapper. + """Changes MSIS location by month strings defined in _msis_wrapper. Args: location (str): month as defined in :class:`NRLMSISE-00.` @@ -679,7 +718,7 @@ def set_season(self, month): self._msis.set_season(month) def set_doy(self, day_of_year): - """ Changes MSIS season by day of year. + """Changes MSIS season by day of year. Args: day_of_year (int): 1. Jan.=0, 1.Feb=32 @@ -688,7 +727,7 @@ def set_doy(self, day_of_year): self._msis.set_doy(day_of_year) def get_temperature(self, h_cm): - """ Returns the temperature of air in K. + """Returns the temperature of air in K. Wraps around ctypes calls to the NRLMSISE-00 C library. @@ -712,26 +751,28 @@ class AIRSAtmosphere(EarthsAtmosphere): """ def __init__(self, location, season, extrapolate=True, *args, **kwargs): - if location != 'SouthPole': - raise Exception(self.__class__.__name__ + - "(): Only South Pole location supported. " + - location) + if location != "SouthPole": + raise Exception( + self.__class__.__name__ + + "(): Only South Pole location supported. " + + location + ) self.extrapolate = extrapolate self.month2doy = { - 'January': 1, - 'February': 32, - 'March': 60, - 'April': 91, - 'May': 121, - 'June': 152, - 'July': 182, - 'August': 213, - 'September': 244, - 'October': 274, - 'November': 305, - 'December': 335 + "January": 1, + "February": 32, + "March": 60, + "April": 91, + "May": 121, + "June": 152, + "July": 182, + "August": 213, + "September": 244, + "October": 274, + "November": 305, + "December": 335, } self.season = season @@ -750,18 +791,20 @@ def init_parameters(self, location, **kwargs): from os import path def bytespdate2num(b): - return datestr2num(b.decode('utf-8')) + return datestr2num(b.decode("utf-8")) - data_path = (join( - path.expanduser('~'), - 'OneDrive/Dokumente/projects/atmospheric_variations/')) + data_path = join( + path.expanduser("~"), "OneDrive/Dokumente/projects/atmospheric_variations/" + ) - if 'table_path' in kwargs: - data_path = kwargs['table_path'] + if "table_path" in kwargs: + data_path = kwargs["table_path"] - files = [('dens', 'airs_amsu_dens_180_daily.txt'), - ('temp', 'airs_amsu_temp_180_daily.txt'), - ('alti', 'airs_amsu_alti_180_daily.txt')] + files = [ + ("dens", "airs_amsu_dens_180_daily.txt"), + ("temp", "airs_amsu_temp_180_daily.txt"), + ("alti", "airs_amsu_alti_180_daily.txt"), + ] data_collection = {} @@ -772,12 +815,12 @@ def bytespdate2num(b): IC79_idx_2 = None for d_key, fname in files: - fname = data_path + 'tables/' + fname + fname = data_path + "tables/" + fname # tabf = open(fname).read() - tab = np.loadtxt(fname, - converters={0: bytespdate2num}, - usecols=[0] + list(range(2, 27))) + tab = np.loadtxt( + fname, converters={0: bytespdate2num}, usecols=[0] + list(range(2, 27)) + ) # with open(fname, 'r') as f: # comline = f.readline() # p_levels = [ @@ -791,41 +834,42 @@ def bytespdate2num(b): elif date.year == 2011: IC79_idx_2 = di surf_val = tab[:, 1] - cols = tab[:, min_press_idx + 2:] + cols = tab[:, min_press_idx + 2 :] data_collection[d_key] = (dates, surf_val, cols) self.interp_tab_d = {} self.interp_tab_t = {} self.dates = {} - dates = data_collection['alti'][0] + dates = data_collection["alti"][0] - msis = MSIS00Atmosphere(location, 'January') + msis = MSIS00Atmosphere(location, "January") for didx, date in enumerate(dates): - h_vec = np.array(data_collection['alti'][2][didx, :] * 1e2) - d_vec = np.array(data_collection['dens'][2][didx, :]) - t_vec = np.array(data_collection['temp'][2][didx, :]) + h_vec = np.array(data_collection["alti"][2][didx, :] * 1e2) + d_vec = np.array(data_collection["dens"][2][didx, :]) + t_vec = np.array(data_collection["temp"][2][didx, :]) if self.extrapolate: # Extrapolate using msis h_extra = np.linspace(h_vec[-1], self.geom.h_atm * 1e2, 250) msis._msis.set_doy(self._get_y_doy(date)[1] - 1) msis_extra_d = np.array([msis.get_density(h) for h in h_extra]) - msis_extra_t = np.array( - [msis.get_temperature(h) for h in h_extra]) + msis_extra_t = np.array([msis.get_temperature(h) for h in h_extra]) # Interpolate last few altitude bins ninterp = 5 for ni in range(ninterp): - cl = (1 - np.exp(-ninterp + ni + 1)) - ch = (1 - np.exp(-ni)) - norm = 1. / (cl + ch) - d_vec[-ni - - 1] = (d_vec[-ni - 1] * cl * norm + - msis.get_density(h_vec[-ni - 1]) * ch * norm) + cl = 1 - np.exp(-ninterp + ni + 1) + ch = 1 - np.exp(-ni) + norm = 1.0 / (cl + ch) + d_vec[-ni - 1] = ( + d_vec[-ni - 1] * cl * norm + + msis.get_density(h_vec[-ni - 1]) * ch * norm + ) t_vec[-ni - 1] = ( - t_vec[-ni - 1] * cl * norm + - msis.get_temperature(h_vec[-ni - 1]) * ch * norm) + t_vec[-ni - 1] * cl * norm + + msis.get_temperature(h_vec[-ni - 1]) * ch * norm + ) # Merge the two datasets h_vec = np.hstack([h_vec[:-1], h_extra]) @@ -867,12 +911,15 @@ def set_season(self, month): def set_IC79_day(self, IC79_day): import datetime + if IC79_day > self.IC79_days: - raise Exception(self.__class__.__name__ + - "::set_IC79_day(): IC79_day above range.") - target_day = self._get_y_doy(self.dates[self.IC79_start] + - datetime.timedelta(days=IC79_day)) - info(2, 'setting IC79_day', IC79_day) + raise Exception( + self.__class__.__name__ + "::set_IC79_day(): IC79_day above range." + ) + target_day = self._get_y_doy( + self.dates[self.IC79_start] + datetime.timedelta(days=IC79_day) + ) + info(2, "setting IC79_day", IC79_day) self.h, self.dens = self.interp_tab_d[target_day] _, self.temp = self.interp_tab_t[target_day] self.date = self.dates[target_day] @@ -883,7 +930,7 @@ def _get_y_doy(self, date): return date.timetuple().tm_year, date.timetuple().tm_yday def get_density(self, h_cm): - """ Returns the density of air in g/cm**3. + """Returns the density of air in g/cm**3. Interpolates table at requested value for previously set year and day of year (doy). @@ -903,7 +950,7 @@ def get_density(self, h_cm): return ret def get_temperature(self, h_cm): - """ Returns the temperature in K. + """Returns the temperature in K. Interpolates table at requested value for previously set year and day of year (doy). @@ -933,17 +980,17 @@ class MSIS00IceCubeCentered(MSIS00Atmosphere): """ def __init__(self, location, season): - if location != 'SouthPole': - info(2, 'location forced to the South Pole') - location = 'SouthPole' + if location != "SouthPole": + info(2, "location forced to the South Pole") + location = "SouthPole" MSIS00Atmosphere.__init__(self, location, season) # Allow for upgoing zenith angles - self.max_theta = 180. + self.max_theta = 180.0 def latitude(self, det_zenith_deg): - """ Returns the geographic latitude of the shower impact point. + """Returns the geographic latitude of the shower impact point. Assumes a spherical earth. The detector is 1948m under the surface. @@ -959,29 +1006,40 @@ def latitude(self, det_zenith_deg): r = self.geom.r_E d = 1948 # m - theta_rad = det_zenith_deg / 180. * np.pi + theta_rad = det_zenith_deg / 180.0 * np.pi - x = (np.sqrt(2. * r * d + ((r - d) * np.cos(theta_rad))**2 - d**2) - - (r - d) * np.cos(theta_rad)) + x = np.sqrt(2.0 * r * d + ((r - d) * np.cos(theta_rad)) ** 2 - d**2) - ( + r - d + ) * np.cos(theta_rad) - return -90. + np.arctan2(x * np.sin(theta_rad), - r - d + x * np.cos(theta_rad)) / np.pi * 180. + return ( + -90.0 + + np.arctan2(x * np.sin(theta_rad), r - d + x * np.cos(theta_rad)) + / np.pi + * 180.0 + ) - def set_theta(self, theta_deg, force_spline_calc=True): + def set_theta(self, theta_deg): - self._msis.set_location_coord(longitude=0., - latitude=self.latitude(theta_deg)) + self._msis.set_location_coord(longitude=0.0, latitude=self.latitude(theta_deg)) info( - 1, 'latitude = {0:5.2f} for zenith angle = {1:5.2f}'.format( - self.latitude(theta_deg), theta_deg)) - if theta_deg > 90.: + 1, + "latitude = {0:5.2f} for zenith angle = {1:5.2f}".format( + self.latitude(theta_deg), theta_deg + ), + ) + downgoing_theta_deg = theta_deg + if theta_deg > 90.0: + downgoing_theta_deg = 180.0 - theta_deg info( - 1, 'theta = {0:5.2f} below horizon. using theta = {1:5.2f}'. - format(theta_deg, 180. - theta_deg)) - theta_deg = 180. - theta_deg - MSIS00Atmosphere.set_theta(self, - theta_deg, - force_spline_calc=force_spline_calc) + 1, + "theta = {0:5.2f} below horizon. using theta = {1:5.2f}".format( + theta_deg, downgoing_theta_deg + ), + ) + MSIS00Atmosphere.set_theta(self, downgoing_theta_deg) + + self.theta_deg = theta_deg class GeneralizedTarget(object): @@ -1007,10 +1065,11 @@ class GeneralizedTarget(object): """ def __init__( - self, - len_target=config.len_target * 1e2, # cm - env_density=config.env_density, # g/cm3 - env_name=config.env_name): + self, + len_target=config.len_target * 1e2, # cm + env_density=config.env_density, # g/cm3 + env_name=config.env_name, + ): self.len_target = len_target self.env_density = env_density @@ -1023,16 +1082,15 @@ def max_den(self): def reset(self): """Resets material list to defaults.""" - self.mat_list = [[ - 0., self.len_target, self.env_density, self.env_name - ]] + self.mat_list = [[0.0, self.len_target, self.env_density, self.env_name]] self._update_variables() def _update_variables(self): """Updates internal variables. Not needed to call by user.""" - self.start_bounds, self.end_bounds, \ - self.densities = list(zip(*self.mat_list))[:-1] + self.start_bounds, self.end_bounds, self.densities = list(zip(*self.mat_list))[ + :-1 + ] self.densities = np.array(self.densities) self.start_bounds = np.array(self.start_bounds) self.end_bounds = np.array(self.end_bounds) @@ -1046,9 +1104,10 @@ def set_length(self, new_length_cm): """ if new_length_cm < self.mat_list[-1][0]: raise Exception( - "GeneralizedTarget::set_length(): " + - "can not set length below lower boundary of last " + - "material.") + "GeneralizedTarget::set_length(): " + + "can not set length below lower boundary of last " + + "material." + ) self.len_target = new_length_cm self.mat_list[-1][1] = new_length_cm self._update_variables() @@ -1066,30 +1125,39 @@ def add_material(self, start_position_cm, density, name): Exception: If requested start_position_cm is not properly defined. """ - if start_position_cm < 0. or start_position_cm > self.len_target: - raise Exception("GeneralizedTarget::add_material(): " + - "distance exceeds target dimensions.") - elif (start_position_cm == self.mat_list[-1][0] - and self.mat_list[-1][-1] == self.env_name): - self.mat_list[-1] = [ - start_position_cm, self.len_target, density, name - ] + if start_position_cm < 0.0 or start_position_cm > self.len_target: + raise Exception( + "GeneralizedTarget::add_material(): " + + "distance exceeds target dimensions." + ) + elif ( + start_position_cm == self.mat_list[-1][0] + and self.mat_list[-1][-1] == self.env_name + ): + self.mat_list[-1] = [start_position_cm, self.len_target, density, name] elif start_position_cm <= self.mat_list[-1][0]: - raise Exception("GeneralizedTarget::add_material(): " + - "start_position_cm is ahead of previous material.") + raise Exception( + "GeneralizedTarget::add_material(): " + + "start_position_cm is ahead of previous material." + ) else: self.mat_list[-1][1] = start_position_cm - self.mat_list.append( - [start_position_cm, self.len_target, density, name]) + self.mat_list.append([start_position_cm, self.len_target, density, name]) - info(2, - ("{0}::add_material(): Material '{1}' added. " + - "location on path {2} to {3} m").format(self.__class__.__name__, - name, - self.mat_list[-1][0], - self.mat_list[-1][1])) + info( + 2, + ( + "{0}::add_material(): Material '{1}' added. " + + "location on path {2} to {3} m" + ).format( + self.__class__.__name__, + name, + self.mat_list[-1][0], + self.mat_list[-1][1], + ), + ) self._update_variables() @@ -1101,8 +1169,10 @@ def set_theta(self, *args): NotImplementedError: always """ - raise NotImplementedError('GeneralizedTarget::set_theta(): Method' + - 'not defined for this target class.') + raise NotImplementedError( + "GeneralizedTarget::set_theta(): Method" + + "not defined for this target class." + ) def _integrate(self): """Walks through material list and computes the depth along the @@ -1114,36 +1184,37 @@ def _integrate(self): """ from scipy.interpolate import UnivariateSpline + self.density_depth = None - self.knots = [0.] - self.X_int = [0.] + self.knots = [0.0] + self.X_int = [0.0] for start, end, density, _ in self.mat_list: self.knots.append(end) self.X_int.append(density * (end - start) + self.X_int[-1]) - self._s_X2h = UnivariateSpline(self.X_int, self.knots, k=1, s=0.) - self._s_h2X = UnivariateSpline(self.knots, self.X_int, k=1, s=0.) + self._s_X2h = UnivariateSpline(self.X_int, self.knots, k=1, s=0.0) + self._s_h2X = UnivariateSpline(self.knots, self.X_int, k=1, s=0.0) self._max_X = self.X_int[-1] @property def s_X2h(self): """Spline for depth at distance.""" - if not hasattr(self, '_s_X2h'): + if not hasattr(self, "_s_X2h"): self._integrate() return self._s_X2h @property def s_h2X(self): """Spline for distance at depth.""" - if not hasattr(self, '_s_h2X'): + if not hasattr(self, "_s_h2X"): self._integrate() return self._s_h2X @property def max_X(self): """Maximal depth of target.""" - if not hasattr(self, '_max_X'): + if not hasattr(self, "_max_X"): self._integrate() return self._max_X @@ -1163,12 +1234,15 @@ def get_density_X(self, X): # allow for some small constant extrapolation for odepack solvers if X[-1] > self.max_X and X[-1] < self.max_X * 1.003: X[-1] = self.max_X - if np.min(X) < 0. or np.max(X) > self.max_X: + if np.min(X) < 0.0 or np.max(X) > self.max_X: # return self.get_density(self.s_X2h(self.max_X)) - info(0, 'Depth {0:4.3f} exceeds target dimensions {1:4.3f}'.format( - np.max(X), self.max_X - )) - raise Exception('Invalid input') + info( + 0, + "Depth {0:4.3f} exceeds target dimensions {1:4.3f}".format( + np.max(X), self.max_X + ), + ) + raise Exception("Invalid input") return self.get_density(self.s_X2h(X)) @@ -1182,7 +1256,7 @@ def r_X2rho(self, X): float: :math:`1/\\rho` in cm**3/g """ - return 1. / self.get_density_X(X) + return 1.0 / self.get_density_X(X) def get_density(self, l_cm): """Returns the density in g/cm**3 as a function of position l in cm. @@ -1200,12 +1274,13 @@ def get_density(self, l_cm): res = np.zeros_like(l_cm) if np.min(l_cm) < 0 or np.max(l_cm) > self.len_target: - raise Exception("GeneralizedTarget::get_density(): " + - "requested position exceeds target legth.") + raise Exception( + "GeneralizedTarget::get_density(): " + + "requested position exceeds target legth." + ) for i, li in enumerate(l_cm): bi = 0 - while not (li >= self.start_bounds[bi] - and li <= self.end_bounds[bi]): + while not (li >= self.start_bounds[bi] and li <= self.end_bounds[bi]): bi += 1 res[i] = self.densities[bi] return res @@ -1228,110 +1303,117 @@ def draw_materials(self, axes=None, logx=False): xend = mat[1] alpha = 0.188 * mat[2] / max(self.densities) + 0.248 if alpha > 1: - alpha = 1. - elif alpha < 0.: - alpha = 0. - axes.fill_between((xstart, xend), (ymax, ymax), (0., 0.), - label=mat[2], - facecolor='grey', - alpha=alpha) + alpha = 1.0 + elif alpha < 0.0: + alpha = 0.0 + axes.fill_between( + (xstart, xend), + (ymax, ymax), + (0.0, 0.0), + label=mat[2], + facecolor="grey", + alpha=alpha, + ) # axes.text(0.5e-2 * (xstart + xend), 0.5 * ymax, str(nm)) - axes.plot([xl for xl in self.knots], self.X_int, lw=1.7, color='r') + axes.plot([xl for xl in self.knots], self.X_int, lw=1.7, color="r") if logx: - axes.set_xscale('log', nonposx='clip') + axes.set_xscale("log", nonposx="clip") - axes.set_ylim(0., ymax) - axes.set_xlabel('distance in target (cm)') - axes.set_ylabel(r'depth X (g/cm$^2)$') + axes.set_ylim(0.0, ymax) + axes.set_xlabel("distance in target (cm)") + axes.set_ylabel(r"depth X (g/cm$^2)$") self.print_table(min_dbg_lev=2) def print_table(self, min_dbg_lev=0): - """Prints table of materials to standard output. - """ + """Prints table of materials to standard output.""" - templ = '{0:^3} | {1:15} | {2:^9.3g} | {3:^9.3g} | {4:^8.5g}' + templ = "{0:^3} | {1:15} | {2:^9.3g} | {3:^9.3g} | {4:^8.5g}" info( min_dbg_lev, - '********************* List of materials ***********************', - no_caller=True) - head = '{0:3} | {1:15} | {2:9} | {3:9} | {4:9}'.format( - 'no', 'name', 'start [cm]', 'end [cm]', 'density [g/cm**3]') - info(min_dbg_lev, '-' * len(head), no_caller=True) + "********************* List of materials ***********************", + no_caller=True, + ) + head = "{0:3} | {1:15} | {2:9} | {3:9} | {4:9}".format( + "no", "name", "start [cm]", "end [cm]", "density [g/cm**3]" + ) + info(min_dbg_lev, "-" * len(head), no_caller=True) info(min_dbg_lev, head, no_caller=True) - info(min_dbg_lev, '-' * len(head), no_caller=True) + info(min_dbg_lev, "-" * len(head), no_caller=True) for nm, mat in enumerate(self.mat_list): - info(min_dbg_lev, - templ.format(nm, mat[3], mat[0], mat[1], mat[2]), - no_caller=True) + info( + min_dbg_lev, + templ.format(nm, mat[3], mat[0], mat[1], mat[2]), + no_caller=True, + ) -if __name__ == '__main__': +if __name__ == "__main__": import matplotlib.pyplot as plt plt.figure(figsize=(5, 4)) - plt.title('CORSIKA atmospheres') + plt.title("CORSIKA atmospheres") cka_atmospheres = [ ("USStd", None), ("BK_USStd", None), ("Karlsruhe", None), - ("ANTARES/KM3NeT-ORCA", 'Summer'), - ("ANTARES/KM3NeT-ORCA", 'Winter'), - ("KM3NeT-ARCA", 'Summer'), - ("KM3NeT-ARCA", 'Winter'), + ("ANTARES/KM3NeT-ORCA", "Summer"), + ("ANTARES/KM3NeT-ORCA", "Winter"), + ("KM3NeT-ARCA", "Summer"), + ("KM3NeT-ARCA", "Winter"), ("KM3NeT", None), - ('SouthPole','December'), - ('PL_SouthPole','January'), - ('PL_SouthPole','August'), + ("SouthPole", "December"), + ("PL_SouthPole", "January"), + ("PL_SouthPole", "August"), ] cka_surf_100 = [] for loc, season in cka_atmospheres: cka_obj = CorsikaAtmosphere(loc, season) cka_obj.set_theta(0.0) x_vec = np.linspace(0, cka_obj.max_X, 5000) - plt.plot(x_vec, - 1 / cka_obj.r_X2rho(x_vec), - lw=1.5, - label='{0}/{1}'.format(loc, season) if season is not None - else '{0}'.format(loc)) - cka_surf_100.append((cka_obj.max_X, 1. / cka_obj.r_X2rho(100.))) + plt.plot( + x_vec, + 1 / cka_obj.r_X2rho(x_vec), + lw=1.5, + label="{0}/{1}".format(loc, season) + if season is not None + else "{0}".format(loc), + ) + cka_surf_100.append((cka_obj.max_X, 1.0 / cka_obj.r_X2rho(100.0))) print(cka_surf_100) - plt.ylabel(r'Density $\rho$ (g/cm$^3$)') - plt.xlabel(r'Depth (g/cm$^2$)') - plt.legend(loc='upper left') + plt.ylabel(r"Density $\rho$ (g/cm$^3$)") + plt.xlabel(r"Depth (g/cm$^2$)") + plt.legend(loc="upper left") plt.tight_layout() plt.figure(figsize=(5, 4)) - plt.title('NRLMSISE-00 atmospheres') + plt.title("NRLMSISE-00 atmospheres") msis_atmospheres = [ - ('SouthPole', "January"), - ('Karlsruhe', "January"), - ('Geneva', "January"), - ('Tokyo', "January"), - ('SanGrasso', "January"), - ('TelAviv', "January"), - ('KSC', "January"), - ('SoudanMine', "January"), - ('Tsukuba', "January"), - ('LynnLake', "January"), - ('PeaceRiver', "January"), - ('FtSumner', "January") + ("SouthPole", "January"), + ("Karlsruhe", "January"), + ("Geneva", "January"), + ("Tokyo", "January"), + ("SanGrasso", "January"), + ("TelAviv", "January"), + ("KSC", "January"), + ("SoudanMine", "January"), + ("Tsukuba", "January"), + ("LynnLake", "January"), + ("PeaceRiver", "January"), + ("FtSumner", "January"), ] msis_surf_100 = [] for loc, season in msis_atmospheres: msis_obj = MSIS00Atmosphere(loc, season) msis_obj.set_theta(0.0) x_vec = np.linspace(0, msis_obj.max_X, 5000) - plt.plot(x_vec, - 1 / msis_obj.r_X2rho(x_vec), - lw=1.5, - label='{0}'.format(loc)) - msis_surf_100.append((msis_obj.max_X, 1. / msis_obj.r_X2rho(100.))) + plt.plot(x_vec, 1 / msis_obj.r_X2rho(x_vec), lw=1.5, label="{0}".format(loc)) + msis_surf_100.append((msis_obj.max_X, 1.0 / msis_obj.r_X2rho(100.0))) print(msis_surf_100) - plt.ylabel(r'Density $\rho$ (g/cm$^3$)') - plt.xlabel(r'Depth (g/cm$^2$)') - plt.legend(loc='upper left') + plt.ylabel(r"Density $\rho$ (g/cm$^3$)") + plt.xlabel(r"Depth (g/cm$^2$)") + plt.legend(loc="upper left") plt.tight_layout() plt.show() From d024fdef67149e9f7784cd80398b8474e43b4333 Mon Sep 17 00:00:00 2001 From: afedynitch Date: Mon, 27 Jun 2022 17:29:19 +0800 Subject: [PATCH 2/2] Version and changelog updated --- CHANGELOG.md | 17 ++++++++++++++--- MCEq/version.py | 2 +- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index db72b22..a88f97a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,16 +1,24 @@ -Changes in MCEq since moving from MCEq_classic to the 1.X.X versions: +# Changes in MCEq since moving from MCEq_classic to the 1.X.X versions: + +Version 1.2.6: + +- Fixed a bug where alternating changes between up- and downgoing zenith angles are not detected correctly and the same results returned Version 1.2.5: + - Migrated to github actions for CI. Thx @jncots Version 1.2.3: + - Binary wheel for aarch64 on linux. Thanks to @obidev Version 1.2.2: + - Added wheels for Python 3.9 on 64 bit systems - Removed binary wheels for 32bit due to lack of h5py wheels and mainstream has transitioned to 64bit. 32bit users can build MCEq from source. Version 1.2.0: + - New data tables: physics will be affected mostly low energies < 30 GeV and minor corrections can be visible for particle ratios at higher energies. [See dedicated doc page](http://mceq.readthedocs.org/en/latest/v12v11_diff.html). @@ -26,7 +34,7 @@ corrections can be visible for particle ratios at higher energies. Version 1.1.3: -- Added atmospheres for KM3NeT by @Kakiczi (https://github.com/Kakiczi) +- Added atmospheres for KM3NeT by @Kakiczi () - new keyword for MCEqRun "build_matrices"=False (to prevent matrix building on init) - Equivalent projectile mappings separated for SIBYLL 2.1 and 2.3 @@ -56,20 +64,24 @@ may have generated some confusion. Other changes include: - tests for atmospheres Version 1.0.9: + - disable_decays flag in advanced options fixed - threshold energy not used in n_mu, n_e - new generic function 'n_particles' for arbitrary particle types - new config option dedx_material Version 1.0.8: + - Fixed a Python3 compatibility issue in density profiles - Cross checked and corrected the functionality of "disabled particles" in config file - Version tagged for paper submission Version 1.0.6 and 1.0.7: + - A few typos corrected Version 1.0.5: + - Check added to make sure depth grids are strictly increasing - Tutorial updated to reflect this fact - New advanced variable in config "stability_margin" @@ -93,4 +105,3 @@ General remark:: one based on DifferentialOperators that are not part of the "C Matrix". With the new structure it became simpler to exchange physical models or track particle decays of any kind. - diff --git a/MCEq/version.py b/MCEq/version.py index 09964d6..593fb6f 100644 --- a/MCEq/version.py +++ b/MCEq/version.py @@ -1 +1 @@ -__version__ = '1.2.5' +__version__ = '1.2.6'