From 7fed0a3b7c1d44d501b64a7de19484ca0594ee1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Fri, 1 Sep 2023 11:18:15 +0200 Subject: [PATCH 1/3] Enable passing arrays of vectors to get_angle function --- radiotools/helper.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/radiotools/helper.py b/radiotools/helper.py index 59fa3cc..9a741a6 100644 --- a/radiotools/helper.py +++ b/radiotools/helper.py @@ -267,8 +267,13 @@ def get_angle(v1, v2): Returns: float or list of floats angle(s) between vector(s) """ - - arccos = np.dot(v1, v2) / (np.linalg.norm(v1.T, axis=0) * np.linalg.norm(v2.T, axis=0)) + + if v1.ndim == 2 and v2.ndim == 2: + arccos = np.array([ + np.dot(v1_, v2_) / (np.linalg.norm(v1_.T, axis=0) * np.linalg.norm(v2_.T, axis=0)) + for v1_, v2_ in zip(v1, v2)]) + else: + arccos = np.dot(v1, v2) / (np.linalg.norm(v1.T, axis=0) * np.linalg.norm(v2.T, axis=0)) # catch numerical overlow mask1 = arccos > 1 mask2 = arccos < -1 From 7a07b6ef9b085d2b4ec7cfffcda319e3615dc3ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Mon, 4 Mar 2024 15:39:23 +0100 Subject: [PATCH 2/3] use type=object for inhomogenous array --- radiotools/atmosphere/refractivity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/radiotools/atmosphere/refractivity.py b/radiotools/atmosphere/refractivity.py index 692d39f..2cb3cea 100644 --- a/radiotools/atmosphere/refractivity.py +++ b/radiotools/atmosphere/refractivity.py @@ -231,8 +231,8 @@ def get_cached_table_curved_atmosphere(self): self._distances.append(distances) self._refractivity_integrated_table.append(refractivity_integrated_table) - self._distances = np.array(self._distances) - self._refractivity_integrated_table = np.array(self._refractivity_integrated_table) + self._distances = np.array(self._distances, dtype=object) + self._refractivity_integrated_table = np.array(self._refractivity_integrated_table, dtype=object) data = {"refractivity_integrated_table": self._refractivity_integrated_table, "distance_increment": self._distance_increment, From b7c9f284fd7fa2c485b6c28d7213876d2ba49346 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Mon, 4 Mar 2024 15:48:23 +0100 Subject: [PATCH 3/3] Only cosmetics: renameing, whitespace, indentation, doc-strings, .... No functional change --- radiotools/atmosphere/refractivity.py | 78 +++++++++++++++------------ 1 file changed, 43 insertions(+), 35 deletions(-) diff --git a/radiotools/atmosphere/refractivity.py b/radiotools/atmosphere/refractivity.py index 2cb3cea..71eab6b 100644 --- a/radiotools/atmosphere/refractivity.py +++ b/radiotools/atmosphere/refractivity.py @@ -3,7 +3,6 @@ import numpy as np import sys -import copy import os @@ -20,14 +19,15 @@ def n_param_ZHAireS(h): Parametrisation for the refractivity N = n - 1, as funtion of height above sea level N(h) as used in ZHAireS. Parameters - ------------ + ---------- h : float Height over sea level in meter Returns - -------- + ------- + refractivity : float - N(h) + Refractive index for given height """ a = 325e-6 @@ -42,27 +42,29 @@ def get_refractivity_between_two_points_numerical(p1, p2, atm_model=None, refrac p1 and p2 need to be in the same coordiate system with the origin at sea level. Parameters - ------------ + ---------- + p1 : array (3,) coordinates in meter. p2 : array (3,) coordinates in meter. atm_model : int Number of the atmospheric model (from radiotools.atmosphere.models) which provides the desity profile rho(h) - to calculate the refractivity via Gladstone–Dale relation: N(h) = N(0) * rho(h) / rho(0). + to calculate the refractivity via Gladstone-Dale relation: N(h) = N(0) * rho(h) / rho(0). Is only used if "table" is not None (default: None) refractivity_at_sea_level : float - Refractivity at earth surface, i.e., N(0) (default: None). Necessary if refractivity is calculated - via Gladstone–Dale relation. Not necessary/used if a RefractivityTable is given. - table : RefractivityTable - If given, used to determine N(h). Instead of using Gladstone–Dale relation. (default: None) + Refractivity at earth surface, i.e., N(0) (default: None). Necessary if refractivity is calculated + via Gladstone-Dale relation. Not necessary/used if a RefractivityTable is given. + table : RefractivityTable + If given, used to determine N(h). Instead of using Gladstone-Dale relation. (default: None) debug : bool If True, prints out debug output (default: False). - Returns - -------- - integrated refractivity : float + ------- + + integrated_refractivity : float + Refractive index integrated along a straight line between two given points """ if table is None and (atm_model is None and refractivity_at_sea_level is None): @@ -105,7 +107,8 @@ def __init__(self, atm_model=atm.default_model, param=False, refractivity_at_sea """ Parameters - ------------ + ---------- + atm_model : int Number of the atmospheric model (from radiotools.atmosphere.models) which provides the desity profile rho(h) to calculate the refractivity via N(h) = N(0) * rho(h) / rho(0) (Only if param is False). @@ -141,7 +144,7 @@ def __init__(self, atm_model=atm.default_model, param=False, refractivity_at_sea self._height_increment = 10 self._max_heigth = 4e4 self._heights = np.arange(0, self._max_heigth, self._height_increment) - + rho0 = atm.get_density(0, allow_negative_heights=False, model=self._atm_model) if self._use_param: self._refractivity_table = np.array([n_param_ZHAireS(h) - 1 for h in self._heights]) @@ -165,21 +168,23 @@ def __init__(self, atm_model=atm.default_model, param=False, refractivity_at_sea def parse_gdas_file(self, gdas_file): with open(gdas_file, "r") as f: lines = f.readlines() - - atm_para = [l.strip("\n").split() for l in lines[1:5]] + + # atm_para = [line.strip("\n").split() for line in lines[1:5]] self._heights = np.zeros(len(lines) - 6) self._refractivity_table = np.zeros(len(lines) - 6) - for idx, l in enumerate(lines[6:]): - h, n = l.strip("\n").split() + for idx, line in enumerate(lines[6:]): + h, n = line.strip("\n").split() self._heights[idx] = float(h) self._refractivity_table[idx] = float(n) - 1 + self._height_increment = self._heights[1] - self._heights[0] self._max_heigth = np.amax(self._heights) - # just take the "layer" closet value to 0 (sea level) if its within 1m. This is stupid but should accurate enought. + # just take the "layer" closet value to 0 (sea level) if its within 1m. This is stupid but should accurate enought. null = np.argmin(np.abs(self._heights)) if np.abs(self._heights[null]) > 1: sys.exit("Could not find refractive index at sea level in gdas profile. stop...") + self._refractivity_at_sea_level = self._refractivity_table[null] @@ -191,6 +196,7 @@ def read_table_from_file(self, fname): self._distances = data["distances"] self._zeniths = data["zeniths"] + def get_cached_table_curved_atmosphere(self): """ Get table of integrated refractivity as function of zenith angle and distance. For curved atmosphere. @@ -214,7 +220,7 @@ def get_cached_table_curved_atmosphere(self): print("Write {} ...".format(fname)) # anchors for table binned in tan. - self._zeniths = np.arctan(np.linspace(np.tan(self._min_zenith), + self._zeniths = np.arctan(np.linspace(np.tan(self._min_zenith), np.tan(self._max_zenith), self._number_of_zenith_bins)) self._distances = [] self._refractivity_integrated_table = [] @@ -243,8 +249,8 @@ def get_cached_table_curved_atmosphere(self): def get_refractivity_for_height_tabulated(self, h): """ - Get refractivity as function of height above ground from pre-calculated table. - Linear interpolation between table anchors. + Get refractivity as function of height above ground from pre-calculated table. + Linear interpolation between table anchors. """ if h == 0: @@ -264,8 +270,8 @@ def get_refractivity_for_height_tabulated(self, h): def get_integrated_refractivity_for_height_tabulated(self, h): """ - Get integrated refractivity as function of height above ground from pre-calculated table. - Linear interpolation between table anchors. + Get integrated refractivity as function of height above ground from pre-calculated table. + Linear interpolation between table anchors. """ if h == 0: return self._refractivity_at_sea_level @@ -286,9 +292,9 @@ def get_integrated_refractivity_for_height_tabulated(self, h): def _get_integrated_refractivity_for_distance(self, d, zenith): """ - Get integrated refractivity as function of the zenith angle and distance from ground - (along a axis with the specified zenith angle) from pre-calculated table. - Takes clostest zenith angle from table to calculated the refractivity. + Get integrated refractivity as function of the zenith angle and distance from ground + (along a axis with the specified zenith angle) from pre-calculated table. + Takes clostest zenith angle from table to calculated the refractivity. """ if not self._curved: sys.exit("Table not available: please specifiy \"curved=True\"") @@ -327,10 +333,11 @@ def get_zenith_bin(self, zenith): def get_integrated_refractivity_for_distance(self, d, zenith): - """ Get integrated refractivity as function of the zenith angle and distance from ground - (along a axis with the specified zenith angle) from pre-calculated table. - If "_interpolate_zenith" is True, a linear interpolation in zenith angle is performed, - otherwise the clostest zenith angle bin is used for the calculation. + """ + Get integrated refractivity as function of the zenith angle and distance from ground + (along a axis with the specified zenith angle) from pre-calculated table. + If "_interpolate_zenith" is True, a linear interpolation in zenith angle is performed, + otherwise the clostest zenith angle bin is used for the calculation. """ if not self._curved: sys.exit("Table not available: please specifiy \"curved=True\"") @@ -420,9 +427,10 @@ def get_refractivity_between_two_points_numerical(self, p1, p2, debug=False): from scipy import constants from collections import defaultdict - tab_flat = RefractivityTable(atm_model=27, curved=False, number_of_zenith_bins=1000) - tab = RefractivityTable(atm_model=27, interpolate_zenith=True, curved=True, number_of_zenith_bins=1000) - at = atm.Atmosphere(27) + atm_model = 27 + tab_flat = RefractivityTable(atm_model=atm_model, curved=False, number_of_zenith_bins=1000) + tab = RefractivityTable(atm_model=atm_model, interpolate_zenith=True, curved=True, number_of_zenith_bins=1000) + at = atm.Atmosphere(atm_model) data = defaultdict(list)