From 2cb299cab62d50430ffdfeb68f4d794bfccdf6a0 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Wed, 25 Oct 2023 12:11:42 +0200
Subject: [PATCH 01/26] abstract function and exponential profile
---
src/ctapipe/atmosphere.py | 20 ++++++++++++++++++++
1 file changed, 20 insertions(+)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index 884e03965a5..40e03f1446d 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -58,6 +58,20 @@ def integral(self, height: u.Quantity) -> u.Quantity:
"""
+ @abc.abstractmethod
+ def height_from_overburden(self, overburden: u.Quantity) -> u.Quantity:
+ r"""Get the height a.s.l. from the mass overburden in the atmosphere.
+ Inverse of the integral function
+
+ ..
+
+ Returns
+ -------
+ u.Quantity["m"]:
+ Height a.s.l. for given overburden
+
+ """
+
def line_of_sight_integral(
self, distance: u.Quantity, zenith_angle=0 * u.deg, output_units=u.g / u.cm**2
):
@@ -192,6 +206,12 @@ def integral(
self.scale_density * self.scale_height * np.exp(-height / self.scale_height)
)
+ @u.quantity_input(overburden=u.g / (u.cm * u.cm))
+ def height_from_overburden(self, overburden) -> u.Quantity:
+ return -self.scale_height * np.log(
+ overburden / (self.scale_height * self.scale_density)
+ )
+
class TableAtmosphereDensityProfile(AtmosphereDensityProfile):
"""Tabular profile from a table that has both the density and it's integral
From 79a9bf34ec2f427f2b07cb88ac32fcff8d6bed32 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Wed, 25 Oct 2023 13:41:02 +0200
Subject: [PATCH 02/26] implemented in tabulated profile
---
src/ctapipe/atmosphere.py | 10 ++++++++++
1 file changed, 10 insertions(+)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index 40e03f1446d..0126823860c 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -278,6 +278,12 @@ def __init__(self, table: Table):
kind="cubic",
)
+ self._height_interp = interp1d(
+ np.log10(self.table["column_density"].to("g cm-2").value),
+ self.table["height"].to("km").value,
+ kind="cubic",
+ )
+
# ensure it can be read back
self.table.meta["TAB_TYPE"] = "ctapipe.atmosphere.TableAtmosphereDensityProfile"
self.table.meta["TAB_VER"] = 1
@@ -294,6 +300,10 @@ def integral(self, height) -> u.Quantity:
10 ** self._col_density_interp(height.to_value(u.km)), u.g / u.cm**2
)
+ @u.quantity_input(overburden=u.g / (u.cm * u.cm))
+ def height_from_overburden(self, overburden) -> u.Quantity:
+ return u.Quantity(self._height_interp(overburden), u.km)
+
def __repr__(self):
return (
f"{self.__class__.__name__}(meta={self.table.meta}, rows={len(self.table)})"
From 9a328b8774027b9b7ed71acd86f1d912f4aeea2d Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Thu, 26 Oct 2023 10:14:09 +0200
Subject: [PATCH 03/26] implemented for five layer atmosphere
---
src/ctapipe/atmosphere.py | 31 ++++++++++++++++++++++++++++++-
1 file changed, 30 insertions(+), 1 deletion(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index 0126823860c..f6ac6f2e907 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -302,7 +302,9 @@ def integral(self, height) -> u.Quantity:
@u.quantity_input(overburden=u.g / (u.cm * u.cm))
def height_from_overburden(self, overburden) -> u.Quantity:
- return u.Quantity(self._height_interp(overburden), u.km)
+ return u.Quantity(
+ self._height_interp(np.log10(overburden.to("g cm-2").value)), u.km
+ )
def __repr__(self):
return (
@@ -320,6 +322,11 @@ def _exponential(h, a, b, c):
return a + b * np.exp(-h / c)
+def _inv_exponential(T, a, b, c):
+ "inverse function for eponetial atmosphere"
+ return -c * np.log((T - a) / b)
+
+
def _d_exponential(h, a, b, c):
"""derivative of exponential atmosphere"""
return -b / c * np.exp(-h / c)
@@ -330,6 +337,11 @@ def _linear(h, a, b, c):
return a - b * h / c
+def _inv_linear(T, a, b, c):
+ "inverse function for linear atmosphere"
+ return -b / c * (T - a)
+
+
def _d_linear(h, a, b, c):
"""derivative of linear atmosphere"""
return -b / c
@@ -365,6 +377,10 @@ def __init__(self, table: Table):
partial(f, a=param_a[i], b=param_b[i], c=param_c[i])
for i, f in enumerate([_exponential] * 4 + [_linear])
]
+ self._inv_funcs = [
+ partial(f, a=param_a[i], b=param_b[i], c=param_c[i])
+ for i, f in enumerate([_inv_linear] + 4 * [_inv_exponential])
+ ]
self._d_funcs = [
partial(f, a=param_a[i], b=param_b[i], c=param_c[i])
for i, f in enumerate([_d_exponential] * 4 + [_d_linear])
@@ -420,6 +436,19 @@ def integral(self, height) -> u.Quantity:
)
).to(u.g / u.cm**2)
+ @u.quantity_input(overburden=u.g / (u.cm * u.cm))
+ def height_from_overburden(self, overburden) -> u.Quantity:
+ overburdens_at_heights = self.integral(self.table["height"].to(u.km))
+ which_func = np.digitize(overburden, overburdens_at_heights) - 1
+ condlist = [which_func == i for i in range(5)]
+ return u.Quantity(
+ np.piecewise(
+ x=overburden,
+ condlist=condlist,
+ funclist=self._inv_funcs,
+ )
+ ).to(u.km)
+
def __repr__(self):
return (
f"{self.__class__.__name__}(meta={self.table.meta}, rows={len(self.table)})"
From e9dd466b5fb16659150640fb4bb97be92e3e5a60 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Thu, 26 Oct 2023 10:14:50 +0200
Subject: [PATCH 04/26] circle test implemented and passing
---
src/ctapipe/tests/test_atmosphere.py | 43 ++++++++++++++++++++++++++++
1 file changed, 43 insertions(+)
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index 2fc0648488c..500809299db 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -133,3 +133,46 @@ def test_against_reference():
0,
atol=1e-5,
)
+
+
+def test_height_overburden_circle(table_profile):
+ """
+ Test that successive application of height to overburden
+ and overburden to height functions return original values
+ """
+
+ # Five-layer atmosphere
+ fit_reference = np.array(
+ [
+ [0.00 * 100000, -140.508, 1178.05, 994186, 0],
+ [9.75 * 100000, -18.4377, 1265.08, 708915, 0],
+ [19.00 * 100000, 0.217565, 1349.22, 636143, 0],
+ [46.00 * 100000, -0.000201796, 703.745, 721128, 0],
+ [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],
+ ]
+ )
+
+ profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
+
+ circle_height_5_layer = profile_5.height_from_overburden(
+ profile_5.integral(47 * u.km)
+ )
+
+ assert np.allclose(circle_height_5_layer, 47 * u.km, rtol=0.0001)
+
+ # Exponential atmosphere
+ density_model = atmo.ExponentialAtmosphereDensityProfile(
+ scale_height=10 * u.km, scale_density=0.00125 * u.g / u.cm**3
+ )
+
+ circle_height_exponential = density_model.height_from_overburden(
+ density_model.integral(47 * u.km)
+ )
+
+ assert np.allclose(circle_height_exponential, 47 * u.km, rtol=0.0001)
+
+ circle_height_table = table_profile.height_from_overburden(
+ table_profile.integral(47 * u.km)
+ )
+
+ assert np.allclose(circle_height_table, 47 * u.km, rtol=0.0001)
From e968f8764e813e8f677a964f5784c21615575b6e Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Thu, 26 Oct 2023 10:28:51 +0200
Subject: [PATCH 05/26] Add changelog
---
docs/changes/2422.feature.rst | 2 ++
1 file changed, 2 insertions(+)
create mode 100644 docs/changes/2422.feature.rst
diff --git a/docs/changes/2422.feature.rst b/docs/changes/2422.feature.rst
new file mode 100644
index 00000000000..f755f4caf4a
--- /dev/null
+++ b/docs/changes/2422.feature.rst
@@ -0,0 +1,2 @@
+Implement the overburden-to height a.s.l. transformation fucntion in the atmosphere module
+and test that round-trip returns original value.
\ No newline at end of file
From f4162e859cb8eb3dc900de13d10ce8e495175a58 Mon Sep 17 00:00:00 2001
From: Georg Schwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 27 Oct 2023 13:24:52 +0200
Subject: [PATCH 06/26] fixed five layer profile function
---
src/ctapipe/atmosphere.py | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index f6ac6f2e907..d8f91abd36f 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -323,7 +323,7 @@ def _exponential(h, a, b, c):
def _inv_exponential(T, a, b, c):
- "inverse function for eponetial atmosphere"
+ "inverse function for exponetial atmosphere"
return -c * np.log((T - a) / b)
@@ -339,7 +339,7 @@ def _linear(h, a, b, c):
def _inv_linear(T, a, b, c):
"inverse function for linear atmosphere"
- return -b / c * (T - a)
+ return -c / b * (T - a)
def _d_linear(h, a, b, c):
@@ -378,7 +378,7 @@ def __init__(self, table: Table):
for i, f in enumerate([_exponential] * 4 + [_linear])
]
self._inv_funcs = [
- partial(f, a=param_a[i], b=param_b[i], c=param_c[i])
+ partial(f, a=param_a[4-i], b=param_b[4-i], c=param_c[4-i])
for i, f in enumerate([_inv_linear] + 4 * [_inv_exponential])
]
self._d_funcs = [
@@ -438,8 +438,8 @@ def integral(self, height) -> u.Quantity:
@u.quantity_input(overburden=u.g / (u.cm * u.cm))
def height_from_overburden(self, overburden) -> u.Quantity:
- overburdens_at_heights = self.integral(self.table["height"].to(u.km))
- which_func = np.digitize(overburden, overburdens_at_heights) - 1
+ overburdens_at_heights = np.flip(self.integral(self.table["height"].to(u.km)))
+ which_func = np.digitize(overburden, overburdens_at_heights)
condlist = [which_func == i for i in range(5)]
return u.Quantity(
np.piecewise(
From 22766636f1cee90390edc8e67c28b6d5f11bc88d Mon Sep 17 00:00:00 2001
From: Georg Schwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 27 Oct 2023 13:25:10 +0200
Subject: [PATCH 07/26] unit test civers all five layers now
---
src/ctapipe/tests/test_atmosphere.py | 14 ++++++++++----
1 file changed, 10 insertions(+), 4 deletions(-)
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index 500809299db..ae2ae56692d 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -154,11 +154,17 @@ def test_height_overburden_circle(table_profile):
profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
- circle_height_5_layer = profile_5.height_from_overburden(
- profile_5.integral(47 * u.km)
- )
+ layer_5_heights=u.Quantity([5,15,30,70,110]*u.km)
+
+ for height in layer_5_heights:
+
+ circle_height_5_layer = profile_5.height_from_overburden(
+ profile_5.integral(height)
+ )
+
+ assert np.allclose(circle_height_5_layer, height, rtol=0.005)
+
- assert np.allclose(circle_height_5_layer, 47 * u.km, rtol=0.0001)
# Exponential atmosphere
density_model = atmo.ExponentialAtmosphereDensityProfile(
From 45fc98ce61673d209cad0eda3bb592993d29cd34 Mon Sep 17 00:00:00 2001
From: Georg Schwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 27 Oct 2023 13:28:42 +0200
Subject: [PATCH 08/26] ran black
---
src/ctapipe/atmosphere.py | 2 +-
src/ctapipe/tests/test_atmosphere.py | 5 +----
2 files changed, 2 insertions(+), 5 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index d8f91abd36f..51a8e70435c 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -378,7 +378,7 @@ def __init__(self, table: Table):
for i, f in enumerate([_exponential] * 4 + [_linear])
]
self._inv_funcs = [
- partial(f, a=param_a[4-i], b=param_b[4-i], c=param_c[4-i])
+ partial(f, a=param_a[4 - i], b=param_b[4 - i], c=param_c[4 - i])
for i, f in enumerate([_inv_linear] + 4 * [_inv_exponential])
]
self._d_funcs = [
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index ae2ae56692d..75afd59a9c4 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -154,18 +154,15 @@ def test_height_overburden_circle(table_profile):
profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
- layer_5_heights=u.Quantity([5,15,30,70,110]*u.km)
+ layer_5_heights = u.Quantity([5, 15, 30, 70, 110] * u.km)
for height in layer_5_heights:
-
circle_height_5_layer = profile_5.height_from_overburden(
profile_5.integral(height)
)
assert np.allclose(circle_height_5_layer, height, rtol=0.005)
-
-
# Exponential atmosphere
density_model = atmo.ExponentialAtmosphereDensityProfile(
scale_height=10 * u.km, scale_density=0.00125 * u.g / u.cm**3
From 4a384fb6b08fe5337f21c303b42b03efd4cd3b89 Mon Sep 17 00:00:00 2001
From: Maximilian Linhoff <maximilian.linhoff@tu-dortmund.de>
Date: Fri, 27 Oct 2023 19:42:20 +0200
Subject: [PATCH 09/26] Improve unit handling in atmosphere module, use fill
values for interpolation
---
src/ctapipe/atmosphere.py | 78 +++++++++++++++-------------
src/ctapipe/tests/test_atmosphere.py | 17 +++++-
2 files changed, 57 insertions(+), 38 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index 51a8e70435c..3b94edbef35 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -15,7 +15,7 @@
import numpy as np
from astropy import units as u
-from astropy.table import Table
+from astropy.table import QTable
from scipy.interpolate import interp1d
__all__ = [
@@ -29,6 +29,9 @@
1,
}
+GRAMMAGE_UNIT = u.g / u.cm**2
+DENSITY_UNIT = u.g / u.cm**3
+
class AtmosphereDensityProfile(abc.ABC):
"""
@@ -73,7 +76,10 @@ def height_from_overburden(self, overburden: u.Quantity) -> u.Quantity:
"""
def line_of_sight_integral(
- self, distance: u.Quantity, zenith_angle=0 * u.deg, output_units=u.g / u.cm**2
+ self,
+ distance: u.Quantity,
+ zenith_angle=0 * u.deg,
+ output_units=GRAMMAGE_UNIT,
):
r"""Line-of-sight integral from the shower distance to infinity, along
the direction specified by the zenith angle. This is sometimes called
@@ -144,7 +150,7 @@ def peek(self):
return fig, axis
@classmethod
- def from_table(cls, table: Table):
+ def from_table(cls, table: QTable):
"""return a subclass of AtmosphereDensityProfile from a serialized
table"""
@@ -191,7 +197,7 @@ class ExponentialAtmosphereDensityProfile(AtmosphereDensityProfile):
"""
scale_height: u.Quantity = 8 * u.km
- scale_density: u.Quantity = 0.00125 * u.g / (u.cm**3)
+ scale_density: u.Quantity = 0.00125 * u.g / u.cm**3
@u.quantity_input(height=u.m)
def __call__(self, height) -> u.Quantity:
@@ -206,7 +212,7 @@ def integral(
self.scale_density * self.scale_height * np.exp(-height / self.scale_height)
)
- @u.quantity_input(overburden=u.g / (u.cm * u.cm))
+ @u.quantity_input(overburden=u.g / u.cm**2)
def height_from_overburden(self, overburden) -> u.Quantity:
return -self.scale_height * np.log(
overburden / (self.scale_height * self.scale_density)
@@ -247,11 +253,11 @@ class TableAtmosphereDensityProfile(AtmosphereDensityProfile):
load a TableAtmosphereDensityProfile from a supported EventSource
"""
- def __init__(self, table: Table):
+ def __init__(self, table: QTable):
"""
Parameters
----------
- table: Table
+ table : QTable
Table with columns `height`, `density`, and `column_density`
"""
@@ -259,29 +265,32 @@ def __init__(self, table: Table):
if col not in table.colnames:
raise ValueError(f"Missing expected column: {col} in table")
- self.table = table[
- (table["height"] >= 0)
- & (table["density"] > 0)
- & (table["column_density"] > 0)
- ]
+ valid = (
+ (table["height"].value >= 0)
+ & (table["density"].value > 0)
+ & (table["column_density"].value > 0)
+ )
+ self.table = QTable(table[valid], copy=False)
# interpolation is done in log-y to minimize spline wobble
+ log_density = np.log10(self.table["density"].to_value(DENSITY_UNIT))
+ log_column_density = np.log10(
+ self.table["column_density"].to_value(GRAMMAGE_UNIT)
+ )
+ height_km = self.table["height"].to_value(u.km)
- self._density_interp = interp1d(
- self.table["height"].to("km").value,
- np.log10(self.table["density"].to("g cm-3").value),
+ interp_kwargs = dict(
kind="cubic",
+ bounds_error=False,
+ )
+ self._density_interp = interp1d(
+ height_km, log_density, fill_value=(np.nan, -np.inf), **interp_kwargs
)
self._col_density_interp = interp1d(
- self.table["height"].to("km").value,
- np.log10(self.table["column_density"].to("g cm-2").value),
- kind="cubic",
+ height_km, log_column_density, fill_value=(np.nan, -np.inf), **interp_kwargs
)
-
self._height_interp = interp1d(
- np.log10(self.table["column_density"].to("g cm-2").value),
- self.table["height"].to("km").value,
- kind="cubic",
+ log_column_density, height_km, fill_value=(np.inf, np.nan), **interp_kwargs
)
# ensure it can be read back
@@ -290,21 +299,18 @@ def __init__(self, table: Table):
@u.quantity_input(height=u.m)
def __call__(self, height) -> u.Quantity:
- return u.Quantity(
- 10 ** self._density_interp(height.to_value(u.km)), u.g / u.cm**3
- )
+ log_density = self._density_interp(height.to_value(u.km))
+ return u.Quantity(10**log_density, DENSITY_UNIT, copy=False)
@u.quantity_input(height=u.m)
def integral(self, height) -> u.Quantity:
- return u.Quantity(
- 10 ** self._col_density_interp(height.to_value(u.km)), u.g / u.cm**2
- )
+ log_col_density = self._col_density_interp(height.to_value(u.km))
+ return u.Quantity(10**log_col_density, GRAMMAGE_UNIT, copy=False)
@u.quantity_input(overburden=u.g / (u.cm * u.cm))
def height_from_overburden(self, overburden) -> u.Quantity:
- return u.Quantity(
- self._height_interp(np.log10(overburden.to("g cm-2").value)), u.km
- )
+ log_overburden = np.log10(overburden.to_value(GRAMMAGE_UNIT))
+ return u.Quantity(self._height_interp(log_overburden), u.km, copy=False)
def __repr__(self):
return (
@@ -365,12 +371,12 @@ class FiveLayerAtmosphereDensityProfile(AtmosphereDensityProfile):
A User’s Guide", 2021, Appendix F
"""
- def __init__(self, table: Table):
+ def __init__(self, table: QTable):
self.table = table
- param_a = self.table["a"].to("g/cm2")
- param_b = self.table["b"].to("g/cm2")
- param_c = self.table["c"].to("km")
+ param_a = self.table["a"].to(GRAMMAGE_UNIT)
+ param_b = self.table["b"].to(GRAMMAGE_UNIT)
+ param_c = self.table["c"].to(u.km)
# build list of column density functions and their derivatives:
self._funcs = [
@@ -396,7 +402,7 @@ def from_array(cls, array: np.ndarray, metadata: dict = None):
if array.shape != (5, 5):
raise ValueError("expected ndarray with shape (5,5)")
- table = Table(
+ table = QTable(
array,
names=["height", "a", "b", "c", "1/c"],
units=["cm", "g/cm2", "g/cm2", "cm", "cm-1"],
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index 75afd59a9c4..6009506251e 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -26,8 +26,8 @@ def get_simtel_profile_from_eventsource():
return source.atmosphere_density_profile
-@pytest.fixture(scope="session", name="table_profile")
-def fixture_table_profile():
+@pytest.fixture(scope="session")
+def table_profile():
"""a table profile for testing"""
return get_simtel_profile_from_eventsource()
@@ -179,3 +179,16 @@ def test_height_overburden_circle(table_profile):
)
assert np.allclose(circle_height_table, 47 * u.km, rtol=0.0001)
+
+
+def test_out_of_range(table_profile):
+ with pytest.warns(RuntimeWarning, match="divide by zero"):
+ assert np.isposinf(table_profile.height_from_overburden(0 * u.g / u.cm**2))
+
+ assert np.isnan(table_profile.height_from_overburden(2000 * u.g / u.cm**2))
+
+ assert table_profile(150 * u.km).value == 0.0
+ assert np.isnan(table_profile(-0.1 * u.km))
+
+ assert table_profile.integral(150 * u.km).value == 0.0
+ assert np.isnan(table_profile.integral(-0.1 * u.km))
From 804b74b68166f79fcc3bb5f9c0e0d4a3734f5f18 Mon Sep 17 00:00:00 2001
From: Georg Schwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Mon, 30 Oct 2023 11:42:32 +0100
Subject: [PATCH 10/26] Edge cases and tests for exp and 5 layer profiles
---
src/ctapipe/atmosphere.py | 46 +++++++++++++++++++---------
src/ctapipe/tests/test_atmosphere.py | 42 ++++++++++++++++++++++++-
2 files changed, 73 insertions(+), 15 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index 3b94edbef35..e3731e24bad 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -201,21 +201,32 @@ class ExponentialAtmosphereDensityProfile(AtmosphereDensityProfile):
@u.quantity_input(height=u.m)
def __call__(self, height) -> u.Quantity:
- return self.scale_density * np.exp(-height / self.scale_height)
+ return np.where(
+ height >= 0 * u.m,
+ self.scale_density * np.exp(-height / self.scale_height),
+ np.nan,
+ )
@u.quantity_input(height=u.m)
def integral(
self,
height,
) -> u.Quantity:
- return (
- self.scale_density * self.scale_height * np.exp(-height / self.scale_height)
+ return np.where(
+ height >= 0 * u.m,
+ self.scale_density
+ * self.scale_height
+ * np.exp(-height / self.scale_height),
+ np.nan,
)
@u.quantity_input(overburden=u.g / u.cm**2)
def height_from_overburden(self, overburden) -> u.Quantity:
- return -self.scale_height * np.log(
- overburden / (self.scale_height * self.scale_density)
+ return np.where(
+ overburden <= self.scale_height * self.scale_density,
+ -self.scale_height
+ * np.log(overburden / (self.scale_height * self.scale_density)),
+ np.nan,
)
@@ -340,12 +351,12 @@ def _d_exponential(h, a, b, c):
def _linear(h, a, b, c):
"""linear atmosphere"""
- return a - b * h / c
+ return np.where(h < a * c / b, a - b * h / c, 0)
def _inv_linear(T, a, b, c):
"inverse function for linear atmosphere"
- return -c / b * (T - a)
+ return np.where(T > 0, -c / b * (T - a), np.inf)
def _d_linear(h, a, b, c):
@@ -353,6 +364,10 @@ def _d_linear(h, a, b, c):
return -b / c
+def _nan_func(x, unit):
+ return np.nan * unit
+
+
class FiveLayerAtmosphereDensityProfile(AtmosphereDensityProfile):
r"""
CORSIKA 5-layer atmosphere model
@@ -383,14 +398,17 @@ def __init__(self, table: QTable):
partial(f, a=param_a[i], b=param_b[i], c=param_c[i])
for i, f in enumerate([_exponential] * 4 + [_linear])
]
+ self._funcs.insert(0, partial(_nan_func, unit=GRAMMAGE_UNIT))
self._inv_funcs = [
partial(f, a=param_a[4 - i], b=param_b[4 - i], c=param_c[4 - i])
for i, f in enumerate([_inv_linear] + 4 * [_inv_exponential])
]
+ self._inv_funcs.append(partial(_nan_func, unit=u.m))
self._d_funcs = [
partial(f, a=param_a[i], b=param_b[i], c=param_c[i])
for i, f in enumerate([_d_exponential] * 4 + [_d_linear])
]
+ self._d_funcs.insert(0, partial(_nan_func, unit=DENSITY_UNIT))
@classmethod
def from_array(cls, array: np.ndarray, metadata: dict = None):
@@ -419,8 +437,8 @@ def from_array(cls, array: np.ndarray, metadata: dict = None):
@u.quantity_input(height=u.m)
def __call__(self, height) -> u.Quantity:
- which_func = np.digitize(height, self.table["height"]) - 1
- condlist = [which_func == i for i in range(5)]
+ which_func = np.digitize(height, self.table["height"])
+ condlist = [which_func == i for i in range(6)]
return u.Quantity(
-1
* np.piecewise(
@@ -428,25 +446,25 @@ def __call__(self, height) -> u.Quantity:
condlist=condlist,
funclist=self._d_funcs,
)
- ).to(u.g / u.cm**3)
+ ).to(DENSITY_UNIT)
@u.quantity_input(height=u.m)
def integral(self, height) -> u.Quantity:
- which_func = np.digitize(height, self.table["height"]) - 1
- condlist = [which_func == i for i in range(5)]
+ which_func = np.digitize(height, self.table["height"])
+ condlist = [which_func == i for i in range(6)]
return u.Quantity(
np.piecewise(
x=height,
condlist=condlist,
funclist=self._funcs,
)
- ).to(u.g / u.cm**2)
+ ).to(GRAMMAGE_UNIT)
@u.quantity_input(overburden=u.g / (u.cm * u.cm))
def height_from_overburden(self, overburden) -> u.Quantity:
overburdens_at_heights = np.flip(self.integral(self.table["height"].to(u.km)))
which_func = np.digitize(overburden, overburdens_at_heights)
- condlist = [which_func == i for i in range(5)]
+ condlist = [which_func == i for i in range(6)]
return u.Quantity(
np.piecewise(
x=overburden,
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index 6009506251e..2c46c84d10e 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -181,7 +181,7 @@ def test_height_overburden_circle(table_profile):
assert np.allclose(circle_height_table, 47 * u.km, rtol=0.0001)
-def test_out_of_range(table_profile):
+def test_out_of_range_table(table_profile):
with pytest.warns(RuntimeWarning, match="divide by zero"):
assert np.isposinf(table_profile.height_from_overburden(0 * u.g / u.cm**2))
@@ -192,3 +192,43 @@ def test_out_of_range(table_profile):
assert table_profile.integral(150 * u.km).value == 0.0
assert np.isnan(table_profile.integral(-0.1 * u.km))
+
+
+def test_out_of_range_exponential():
+ density_model = atmo.ExponentialAtmosphereDensityProfile(
+ scale_height=10 * u.km, scale_density=0.00125 * u.g / u.cm**3
+ )
+
+ with pytest.warns(RuntimeWarning, match="divide by zero"):
+ assert np.isposinf(density_model.height_from_overburden(0 * u.g / u.cm**2))
+
+ assert np.isnan(density_model.height_from_overburden(2000 * u.g / u.cm**2))
+
+ assert np.isnan(density_model(-0.1 * u.km))
+
+ assert np.isnan(density_model.integral(-0.1 * u.km))
+
+
+def test_out_of_range_five_layer():
+ # Five-layer atmosphere
+ fit_reference = np.array(
+ [
+ [0.00 * 100000, -140.508, 1178.05, 994186, 0],
+ [9.75 * 100000, -18.4377, 1265.08, 708915, 0],
+ [19.00 * 100000, 0.217565, 1349.22, 636143, 0],
+ [46.00 * 100000, -0.000201796, 703.745, 721128, 0],
+ [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],
+ ]
+ )
+
+ profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
+
+ assert np.isposinf(profile_5.height_from_overburden(0 * u.g / u.cm**2))
+
+ assert np.isnan(profile_5.height_from_overburden(2000 * u.g / u.cm**2))
+
+ assert np.allclose(profile_5(150 * u.km).value, 0.0, atol=1e-9)
+ assert np.isnan(profile_5(-0.1 * u.km))
+
+ assert np.allclose(profile_5.integral(150 * u.km).value, 0.0, atol=1e-9)
+ assert np.isnan(profile_5.integral(-0.1 * u.km))
From 7f3da524eb6678c840fb1be7e9d4dd7f028d5598 Mon Sep 17 00:00:00 2001
From: Georg Schwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Thu, 2 Nov 2023 09:13:37 +0100
Subject: [PATCH 11/26] Implement function to get height from slant depth
---
src/ctapipe/atmosphere.py | 24 ++++++++++++++++++++++++
1 file changed, 24 insertions(+)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index e3731e24bad..5955657621e 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -102,7 +102,31 @@ def line_of_sight_integral(
return (
self.integral(distance * np.cos(zenith_angle)) / np.cos(zenith_angle)
).to(output_units)
+
+ def height_from_slant_depth(
+ self,
+ slant_depth: u.Quantity,
+ zenith_angle=0 * u.deg,
+ output_units=u.m,
+ ):
+ r"""Calculates height a.s.l. in the atmosphere from traversed slant depth
+ taking into account the shower zenith angle.
+
+ Parameters
+ ----------
+ slant_depth: u.Quantity["grammage"]
+ line-of-site distance from observer to point
+ zenith_angle: u.Quantity["angle"]
+ zenith angle of observation
+ output_units: u.Unit
+ unit to output (must be convertible to m)
+ """
+
+ return (
+ self.height_from_overburden(slant_depth * np.cos(zenith_angle))
+ ).to(output_units)
+
def peek(self):
"""
Draw quick plot of profile
From 6a19c116fd0facd669d44535179d7896464cc836 Mon Sep 17 00:00:00 2001
From: Georg Schwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Thu, 2 Nov 2023 09:17:33 +0100
Subject: [PATCH 12/26] ran black
---
src/ctapipe/atmosphere.py | 22 +++++++++++-----------
1 file changed, 11 insertions(+), 11 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index 5955657621e..f6817d0ec9f 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -29,8 +29,8 @@
1,
}
-GRAMMAGE_UNIT = u.g / u.cm**2
-DENSITY_UNIT = u.g / u.cm**3
+GRAMMAGE_UNIT = u.g / u.cm ** 2
+DENSITY_UNIT = u.g / u.cm ** 3
class AtmosphereDensityProfile(abc.ABC):
@@ -102,7 +102,7 @@ def line_of_sight_integral(
return (
self.integral(distance * np.cos(zenith_angle)) / np.cos(zenith_angle)
).to(output_units)
-
+
def height_from_slant_depth(
self,
slant_depth: u.Quantity,
@@ -123,10 +123,10 @@ def height_from_slant_depth(
"""
- return (
- self.height_from_overburden(slant_depth * np.cos(zenith_angle))
- ).to(output_units)
-
+ return (self.height_from_overburden(slant_depth * np.cos(zenith_angle))).to(
+ output_units
+ )
+
def peek(self):
"""
Draw quick plot of profile
@@ -221,7 +221,7 @@ class ExponentialAtmosphereDensityProfile(AtmosphereDensityProfile):
"""
scale_height: u.Quantity = 8 * u.km
- scale_density: u.Quantity = 0.00125 * u.g / u.cm**3
+ scale_density: u.Quantity = 0.00125 * u.g / u.cm ** 3
@u.quantity_input(height=u.m)
def __call__(self, height) -> u.Quantity:
@@ -244,7 +244,7 @@ def integral(
np.nan,
)
- @u.quantity_input(overburden=u.g / u.cm**2)
+ @u.quantity_input(overburden=u.g / u.cm ** 2)
def height_from_overburden(self, overburden) -> u.Quantity:
return np.where(
overburden <= self.scale_height * self.scale_density,
@@ -335,12 +335,12 @@ def __init__(self, table: QTable):
@u.quantity_input(height=u.m)
def __call__(self, height) -> u.Quantity:
log_density = self._density_interp(height.to_value(u.km))
- return u.Quantity(10**log_density, DENSITY_UNIT, copy=False)
+ return u.Quantity(10 ** log_density, DENSITY_UNIT, copy=False)
@u.quantity_input(height=u.m)
def integral(self, height) -> u.Quantity:
log_col_density = self._col_density_interp(height.to_value(u.km))
- return u.Quantity(10**log_col_density, GRAMMAGE_UNIT, copy=False)
+ return u.Quantity(10 ** log_col_density, GRAMMAGE_UNIT, copy=False)
@u.quantity_input(overburden=u.g / (u.cm * u.cm))
def height_from_overburden(self, overburden) -> u.Quantity:
From a4fd39b4522856271c990d677d4c503613953d76 Mon Sep 17 00:00:00 2001
From: Georg Schwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Thu, 2 Nov 2023 17:06:50 +0100
Subject: [PATCH 13/26] Formtted again
---
src/ctapipe/atmosphere.py | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index f6817d0ec9f..c769a0d1d00 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -29,8 +29,8 @@
1,
}
-GRAMMAGE_UNIT = u.g / u.cm ** 2
-DENSITY_UNIT = u.g / u.cm ** 3
+GRAMMAGE_UNIT = u.g / u.cm**2
+DENSITY_UNIT = u.g / u.cm**3
class AtmosphereDensityProfile(abc.ABC):
@@ -244,7 +244,7 @@ def integral(
np.nan,
)
- @u.quantity_input(overburden=u.g / u.cm ** 2)
+ @u.quantity_input(overburden=u.g / u.cm**2)
def height_from_overburden(self, overburden) -> u.Quantity:
return np.where(
overburden <= self.scale_height * self.scale_density,
@@ -335,12 +335,12 @@ def __init__(self, table: QTable):
@u.quantity_input(height=u.m)
def __call__(self, height) -> u.Quantity:
log_density = self._density_interp(height.to_value(u.km))
- return u.Quantity(10 ** log_density, DENSITY_UNIT, copy=False)
+ return u.Quantity(10**log_density, DENSITY_UNIT, copy=False)
@u.quantity_input(height=u.m)
def integral(self, height) -> u.Quantity:
log_col_density = self._col_density_interp(height.to_value(u.km))
- return u.Quantity(10 ** log_col_density, GRAMMAGE_UNIT, copy=False)
+ return u.Quantity(10**log_col_density, GRAMMAGE_UNIT, copy=False)
@u.quantity_input(overburden=u.g / (u.cm * u.cm))
def height_from_overburden(self, overburden) -> u.Quantity:
From 99bf1ab675594869cf8e3722393f154be762c7c6 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 3 Nov 2023 10:47:59 +0100
Subject: [PATCH 14/26] correct pre-commit hooks
---
src/ctapipe/atmosphere.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index c769a0d1d00..07293e046e8 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -221,7 +221,7 @@ class ExponentialAtmosphereDensityProfile(AtmosphereDensityProfile):
"""
scale_height: u.Quantity = 8 * u.km
- scale_density: u.Quantity = 0.00125 * u.g / u.cm ** 3
+ scale_density: u.Quantity = 0.00125 * u.g / u.cm**3
@u.quantity_input(height=u.m)
def __call__(self, height) -> u.Quantity:
From 9c74c348bf26bce9daaeccec929952d8b49e5771 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 3 Nov 2023 11:47:49 +0100
Subject: [PATCH 15/26] Include observer altitude in LoS integral
---
src/ctapipe/atmosphere.py | 8 ++++++--
1 file changed, 6 insertions(+), 2 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index 07293e046e8..f469127b00a 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -79,6 +79,7 @@ def line_of_sight_integral(
self,
distance: u.Quantity,
zenith_angle=0 * u.deg,
+ observer_altitude=0 * u.m,
output_units=GRAMMAGE_UNIT,
):
r"""Line-of-sight integral from the shower distance to infinity, along
@@ -86,7 +87,7 @@ def line_of_sight_integral(
the *slant depth*. The atmosphere here is assumed to be Cartesian, the
curvature of the Earth is not taken into account.
- .. math:: X(h, \Psi) = \int_{h}^{\infty} \rho(h' \cos{\Psi}) dh'
+ .. math:: X(d, \Psi) = \int_{h_{obs} + d\cos{\Psi}}^{\infty} \rho(h') dh' / \cos{\Psi}
Parameters
----------
@@ -94,13 +95,16 @@ def line_of_sight_integral(
line-of-site distance from observer to point
zenith_angle: u.Quantity["angle"]
zenith angle of observation
+ observer_altitude: u.Quantity["length]
+ Height a.s.l. of the observer
output_units: u.Unit
unit to output (must be convertible to g/cm2)
"""
return (
- self.integral(distance * np.cos(zenith_angle)) / np.cos(zenith_angle)
+ self.integral(distance * np.cos(zenith_angle) + observer_altitude)
+ / np.cos(zenith_angle)
).to(output_units)
def height_from_slant_depth(
From 941c8eb500191fd0a10810175488690bec7d1a9e Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 3 Nov 2023 11:48:36 +0100
Subject: [PATCH 16/26] Test LoS integral height and zenith dependence
---
src/ctapipe/tests/test_atmosphere.py | 14 ++++++++++++++
1 file changed, 14 insertions(+)
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index 2c46c84d10e..2ca7ef79cf8 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -135,6 +135,20 @@ def test_against_reference():
)
+def test_line_of_sight_integral(table_profile):
+ """
+ Test observer altitude and zenith angle dependence of LoS integral function
+ """
+
+ assert table_profile.line_of_sight_integral(
+ 10 * u.km, observer_altitude=2000 * u.m
+ ) < table_profile.line_of_sight_integral(10 * u.km, observer_altitude=1000 * u.m)
+
+ assert table_profile.line_of_sight_integral(
+ 10 * u.km, zenith_angle=30 * u.deg
+ ) > table_profile.line_of_sight_integral(10 * u.km, zenith_angle=0 * u.deg)
+
+
def test_height_overburden_circle(table_profile):
"""
Test that successive application of height to overburden
From 51f13be5ed9a8be4089467b1ef33540693477078 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Mon, 6 Nov 2023 18:59:24 +0100
Subject: [PATCH 17/26] Test now using global units from atmosphere module
---
src/ctapipe/tests/test_atmosphere.py | 22 +++++++++++-----------
1 file changed, 11 insertions(+), 11 deletions(-)
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index 2ca7ef79cf8..7a2a51e3961 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -65,7 +65,7 @@ def test_models(density_model):
# check we can also compute the integral
column_density = density_model.integral(10 * u.km)
- assert column_density.unit.is_equivalent(u.g / u.cm**2)
+ assert column_density.unit.is_equivalent(atmo.GRAMMAGE_UNIT)
assert np.isclose(
density_model.integral(1 * u.km), density_model.integral(1000 * u.m)
@@ -79,9 +79,9 @@ def test_exponential_model():
"""check exponential models"""
density_model = atmo.ExponentialAtmosphereDensityProfile(
- scale_height=10 * u.m, scale_density=0.00125 * u.g / u.cm**3
+ scale_height=10 * u.m, scale_density=0.00125 * atmo.DENSITY_UNIT
)
- assert np.isclose(density_model(1_000_000 * u.km), 0 * u.g / u.cm**3)
+ assert np.isclose(density_model(1_000_000 * u.km), 0 * atmo.DENSITY_UNIT)
assert np.isclose(density_model(0 * u.km), density_model.scale_density)
@@ -179,7 +179,7 @@ def test_height_overburden_circle(table_profile):
# Exponential atmosphere
density_model = atmo.ExponentialAtmosphereDensityProfile(
- scale_height=10 * u.km, scale_density=0.00125 * u.g / u.cm**3
+ scale_height=10 * u.km, scale_density=0.00125 * atmo.DENSITY_UNIT
)
circle_height_exponential = density_model.height_from_overburden(
@@ -197,9 +197,9 @@ def test_height_overburden_circle(table_profile):
def test_out_of_range_table(table_profile):
with pytest.warns(RuntimeWarning, match="divide by zero"):
- assert np.isposinf(table_profile.height_from_overburden(0 * u.g / u.cm**2))
+ assert np.isposinf(table_profile.height_from_overburden(0 * atmo.GRAMMAGE_UNIT))
- assert np.isnan(table_profile.height_from_overburden(2000 * u.g / u.cm**2))
+ assert np.isnan(table_profile.height_from_overburden(2000 * atmo.GRAMMAGE_UNIT))
assert table_profile(150 * u.km).value == 0.0
assert np.isnan(table_profile(-0.1 * u.km))
@@ -210,13 +210,13 @@ def test_out_of_range_table(table_profile):
def test_out_of_range_exponential():
density_model = atmo.ExponentialAtmosphereDensityProfile(
- scale_height=10 * u.km, scale_density=0.00125 * u.g / u.cm**3
+ scale_height=10 * u.km, scale_density=0.00125 * atmo.DENSITY_UNIT
)
with pytest.warns(RuntimeWarning, match="divide by zero"):
- assert np.isposinf(density_model.height_from_overburden(0 * u.g / u.cm**2))
+ assert np.isposinf(density_model.height_from_overburden(0 * atmo.GRAMMAGE_UNIT))
- assert np.isnan(density_model.height_from_overburden(2000 * u.g / u.cm**2))
+ assert np.isnan(density_model.height_from_overburden(2000 * atmo.GRAMMAGE_UNIT))
assert np.isnan(density_model(-0.1 * u.km))
@@ -237,9 +237,9 @@ def test_out_of_range_five_layer():
profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
- assert np.isposinf(profile_5.height_from_overburden(0 * u.g / u.cm**2))
+ assert np.isposinf(profile_5.height_from_overburden(0 * atmo.GRAMMAGE_UNIT))
- assert np.isnan(profile_5.height_from_overburden(2000 * u.g / u.cm**2))
+ assert np.isnan(profile_5.height_from_overburden(2000 * atmo.GRAMMAGE_UNIT))
assert np.allclose(profile_5(150 * u.km).value, 0.0, atol=1e-9)
assert np.isnan(profile_5(-0.1 * u.km))
From a950213997f9fe2398526e1d6007b6c199057168 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Tue, 7 Nov 2023 09:32:17 +0100
Subject: [PATCH 18/26] Test value of LoS integral too
---
src/ctapipe/tests/test_atmosphere.py | 23 +++++++++++++++++++++++
1 file changed, 23 insertions(+)
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index 7a2a51e3961..d8ce59a7fb4 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -148,6 +148,29 @@ def test_line_of_sight_integral(table_profile):
10 * u.km, zenith_angle=30 * u.deg
) > table_profile.line_of_sight_integral(10 * u.km, zenith_angle=0 * u.deg)
+ los_integral_h_10 = table_profile.line_of_sight_integral(
+ 0 * u.km, observer_altitude=table_profile.table["height"].to("km")[5:15]
+ )
+
+ assert np.allclose(
+ los_integral_h_10,
+ table_profile.table["column_density"].to("g cm-2")[5:15],
+ rtol=1e-3,
+ )
+
+ los_integral_z_20 = table_profile.line_of_sight_integral(
+ table_profile.table["height"].to("km")[5:15] / np.cos(np.deg2rad(20)),
+ observer_altitude=0 * u.m,
+ zenith_angle=20 * u.deg,
+ )
+
+ assert np.allclose(
+ los_integral_z_20,
+ table_profile.table["column_density"].to("g cm-2")[5:15]
+ / np.cos(np.deg2rad(20)),
+ rtol=1e-3,
+ )
+
def test_height_overburden_circle(table_profile):
"""
From 3c2304bd59777ddfb5519373a0d085d06d0c1867 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 23 Feb 2024 14:53:48 +0100
Subject: [PATCH 19/26] classes take Tables, work with QTables internally
---
src/ctapipe/atmosphere.py | 21 ++++++++++++---------
1 file changed, 12 insertions(+), 9 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index f469127b00a..6c0ddaf7632 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -15,7 +15,7 @@
import numpy as np
from astropy import units as u
-from astropy.table import QTable
+from astropy.table import QTable, Table
from scipy.interpolate import interp1d
__all__ = [
@@ -178,10 +178,12 @@ def peek(self):
return fig, axis
@classmethod
- def from_table(cls, table: QTable):
+ def from_table(cls, table: Table):
"""return a subclass of AtmosphereDensityProfile from a serialized
table"""
+ table = QTable(table, copy=False)
+
if "TAB_TYPE" not in table.meta:
raise ValueError("expected a TAB_TYPE metadata field")
@@ -292,11 +294,11 @@ class TableAtmosphereDensityProfile(AtmosphereDensityProfile):
load a TableAtmosphereDensityProfile from a supported EventSource
"""
- def __init__(self, table: QTable):
+ def __init__(self, table: Table):
"""
Parameters
----------
- table : QTable
+ table : Table
Table with columns `height`, `density`, and `column_density`
"""
@@ -305,9 +307,9 @@ def __init__(self, table: QTable):
raise ValueError(f"Missing expected column: {col} in table")
valid = (
- (table["height"].value >= 0)
- & (table["density"].value > 0)
- & (table["column_density"].value > 0)
+ (table["height"] >= 0)
+ & (table["density"] > 0)
+ & (table["column_density"] > 0)
)
self.table = QTable(table[valid], copy=False)
@@ -414,8 +416,9 @@ class FiveLayerAtmosphereDensityProfile(AtmosphereDensityProfile):
A User’s Guide", 2021, Appendix F
"""
- def __init__(self, table: QTable):
- self.table = table
+ def __init__(self, table: Table):
+
+ self.table = QTable(table, copy=False)
param_a = self.table["a"].to(GRAMMAGE_UNIT)
param_b = self.table["b"].to(GRAMMAGE_UNIT)
From 1067eb066c8d31b5a42170c562a900be9546a260 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 23 Feb 2024 15:43:30 +0100
Subject: [PATCH 20/26] Implement slant_depth_from_height function
---
src/ctapipe/atmosphere.py | 37 +++++------
src/ctapipe/tests/test_atmosphere.py | 93 ++++++++++++++++++++++++----
2 files changed, 96 insertions(+), 34 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index 6c0ddaf7632..ff265909f20 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -75,37 +75,32 @@ def height_from_overburden(self, overburden: u.Quantity) -> u.Quantity:
"""
- def line_of_sight_integral(
+ def slant_depth_from_height(
self,
- distance: u.Quantity,
+ height: u.Quantity,
zenith_angle=0 * u.deg,
- observer_altitude=0 * u.m,
output_units=GRAMMAGE_UNIT,
):
- r"""Line-of-sight integral from the shower distance to infinity, along
+ r"""Line-of-sight integral from height to infinity, along
the direction specified by the zenith angle. This is sometimes called
the *slant depth*. The atmosphere here is assumed to be Cartesian, the
- curvature of the Earth is not taken into account.
+ curvature of the Earth is not taken into account. Inverse of
+ height_from_slant_depth function.
- .. math:: X(d, \Psi) = \int_{h_{obs} + d\cos{\Psi}}^{\infty} \rho(h') dh' / \cos{\Psi}
+ .. math:: X(h, \Psi) = \int_{h}^{\infty} \rho(h') dh' / \cos{\Psi}
Parameters
----------
- distance: u.Quantity["length"]
- line-of-site distance from observer to point
+ height: u.Quantity["length"]
+ height a.s.l. at which to start integral
zenith_angle: u.Quantity["angle"]
zenith angle of observation
- observer_altitude: u.Quantity["length]
- Height a.s.l. of the observer
output_units: u.Unit
unit to output (must be convertible to g/cm2)
"""
- return (
- self.integral(distance * np.cos(zenith_angle) + observer_altitude)
- / np.cos(zenith_angle)
- ).to(output_units)
+ return (self.integral(height) / np.cos(zenith_angle)).to(output_units)
def height_from_slant_depth(
self,
@@ -150,21 +145,21 @@ def peek(self):
axis[0].set_ylabel(f"Density / {density.unit.to_string('latex')}")
axis[0].grid(True)
- distance = np.linspace(1, 100, 500) * u.km
+ heights = np.linspace(1, 100, 500) * u.km
for zenith_angle in [0, 40, 50, 70] * u.deg:
- column_density = self.line_of_sight_integral(distance, zenith_angle)
- axis[1].plot(distance, column_density, label=f"$\\Psi$={zenith_angle}")
+ column_density = self.slant_depth_from_height(heights, zenith_angle)
+ axis[1].plot(heights, column_density, label=f"$\\Psi$={zenith_angle}")
axis[1].legend(loc="best")
- axis[1].set_xlabel(f"Distance / {distance.unit.to_string('latex')}")
+ axis[1].set_xlabel(f"Distance / {heights.unit.to_string('latex')}")
axis[1].set_ylabel(f"Column Density / {column_density.unit.to_string('latex')}")
axis[1].set_yscale("log")
axis[1].grid(True)
zenith_angle = np.linspace(0, 80, 20) * u.deg
- for distance in [0, 5, 10, 20] * u.km:
- column_density = self.line_of_sight_integral(distance, zenith_angle)
- axis[2].plot(zenith_angle, column_density, label=f"Height={distance}")
+ for height in [0, 5, 10, 20] * u.km:
+ column_density = self.slant_depth_from_height(height, zenith_angle)
+ axis[2].plot(zenith_angle, column_density, label=f"Height={height}")
axis[2].legend(loc="best")
axis[2].set_xlabel(
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index d8ce59a7fb4..d067b683991 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -129,38 +129,37 @@ def test_against_reference():
1.0 - profile_5(height) / reference_table["rho_5"], 0, atol=1e-5
)
np.testing.assert_allclose(
- 1.0 - profile_5.line_of_sight_integral(height) / reference_table["thick_5"],
+ 1.0 - profile_5.slant_depth_from_height(height) / reference_table["thick_5"],
0,
atol=1e-5,
)
-def test_line_of_sight_integral(table_profile):
+def test_slant_depth_from_height(table_profile):
"""
Test observer altitude and zenith angle dependence of LoS integral function
"""
- assert table_profile.line_of_sight_integral(
- 10 * u.km, observer_altitude=2000 * u.m
- ) < table_profile.line_of_sight_integral(10 * u.km, observer_altitude=1000 * u.m)
+ assert table_profile.slant_depth_from_height(
+ 20 * u.km
+ ) < table_profile.slant_depth_from_height(10 * u.km)
- assert table_profile.line_of_sight_integral(
+ assert table_profile.slant_depth_from_height(
10 * u.km, zenith_angle=30 * u.deg
- ) > table_profile.line_of_sight_integral(10 * u.km, zenith_angle=0 * u.deg)
+ ) > table_profile.slant_depth_from_height(10 * u.km, zenith_angle=0 * u.deg)
- los_integral_h_10 = table_profile.line_of_sight_integral(
- 0 * u.km, observer_altitude=table_profile.table["height"].to("km")[5:15]
+ los_integral_z_0 = table_profile.slant_depth_from_height(
+ table_profile.table["height"].to("km")[5:15]
)
assert np.allclose(
- los_integral_h_10,
+ los_integral_z_0,
table_profile.table["column_density"].to("g cm-2")[5:15],
rtol=1e-3,
)
- los_integral_z_20 = table_profile.line_of_sight_integral(
- table_profile.table["height"].to("km")[5:15] / np.cos(np.deg2rad(20)),
- observer_altitude=0 * u.m,
+ los_integral_z_20 = table_profile.slant_depth_from_height(
+ table_profile.table["height"].to("km")[5:15],
zenith_angle=20 * u.deg,
)
@@ -218,6 +217,74 @@ def test_height_overburden_circle(table_profile):
assert np.allclose(circle_height_table, 47 * u.km, rtol=0.0001)
+def test_height_slant_depth_circle(table_profile):
+ """
+ Test that successive application of height to overburden
+ and overburden to height functions return original values
+ """
+
+ # Five-layer atmosphere
+ fit_reference = np.array(
+ [
+ [0.00 * 100000, -140.508, 1178.05, 994186, 0],
+ [9.75 * 100000, -18.4377, 1265.08, 708915, 0],
+ [19.00 * 100000, 0.217565, 1349.22, 636143, 0],
+ [46.00 * 100000, -0.000201796, 703.745, 721128, 0],
+ [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],
+ ]
+ )
+
+ profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
+
+ layer_5_heights = u.Quantity([5, 15, 30, 70, 110] * u.km)
+
+ for height in layer_5_heights:
+ circle_height_5_layer = profile_5.height_from_slant_depth(
+ profile_5.slant_depth_from_height(height)
+ )
+
+ assert np.allclose(circle_height_5_layer, height, rtol=0.005)
+
+ for height in layer_5_heights:
+ circle_height_5_layer_z_20 = profile_5.height_from_slant_depth(
+ profile_5.slant_depth_from_height(height, zenith_angle=20 * u.deg),
+ zenith_angle=20 * u.deg,
+ )
+
+ assert np.allclose(circle_height_5_layer_z_20, height, rtol=0.005)
+
+ # Exponential atmosphere
+ density_model = atmo.ExponentialAtmosphereDensityProfile(
+ scale_height=10 * u.km, scale_density=0.00125 * atmo.DENSITY_UNIT
+ )
+
+ circle_height_exponential = density_model.height_from_slant_depth(
+ density_model.slant_depth_from_height(47 * u.km)
+ )
+
+ assert np.allclose(circle_height_exponential, 47 * u.km, rtol=0.0001)
+
+ circle_height_exponential_z_20 = density_model.height_from_slant_depth(
+ density_model.slant_depth_from_height(47 * u.km, zenith_angle=20 * u.deg),
+ zenith_angle=20 * u.deg,
+ )
+
+ assert np.allclose(circle_height_exponential_z_20, 47 * u.km, rtol=0.0001)
+
+ circle_height_table = table_profile.height_from_slant_depth(
+ table_profile.slant_depth_from_height(47 * u.km)
+ )
+
+ assert np.allclose(circle_height_table, 47 * u.km, rtol=0.0001)
+
+ circle_height_table_z_20 = table_profile.height_from_slant_depth(
+ table_profile.slant_depth_from_height(47 * u.km, zenith_angle=20 * u.deg),
+ zenith_angle=20 * u.deg,
+ )
+
+ assert np.allclose(circle_height_table_z_20, 47 * u.km, rtol=0.0001)
+
+
def test_out_of_range_table(table_profile):
with pytest.warns(RuntimeWarning, match="divide by zero"):
assert np.isposinf(table_profile.height_from_overburden(0 * atmo.GRAMMAGE_UNIT))
From aa1fb9c7ffc55b40aaacf2dc38ff8e773dcebaca Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Thu, 29 Feb 2024 11:27:57 +0100
Subject: [PATCH 21/26] remove superfluous dots
---
src/ctapipe/atmosphere.py | 2 --
1 file changed, 2 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index ff265909f20..071c7d613b7 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -66,8 +66,6 @@ def height_from_overburden(self, overburden: u.Quantity) -> u.Quantity:
r"""Get the height a.s.l. from the mass overburden in the atmosphere.
Inverse of the integral function
- ..
-
Returns
-------
u.Quantity["m"]:
From 4ced98a65759b60495cce3bdb86c4e0c90e726d1 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Thu, 4 Apr 2024 18:03:46 +0200
Subject: [PATCH 22/26] Rename test functions
---
src/ctapipe/tests/test_atmosphere.py | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index d067b683991..021ba969007 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -171,7 +171,7 @@ def test_slant_depth_from_height(table_profile):
)
-def test_height_overburden_circle(table_profile):
+def test_height_overburden_roundtrip(table_profile):
"""
Test that successive application of height to overburden
and overburden to height functions return original values
@@ -217,7 +217,7 @@ def test_height_overburden_circle(table_profile):
assert np.allclose(circle_height_table, 47 * u.km, rtol=0.0001)
-def test_height_slant_depth_circle(table_profile):
+def test_height_slant_depth_roundtrip(table_profile):
"""
Test that successive application of height to overburden
and overburden to height functions return original values
From 5869a428bc656026dcb618c89400a4ba5cdb8465 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 5 Apr 2024 09:28:22 +0200
Subject: [PATCH 23/26] Implement warning for very large zenith angles
---
src/ctapipe/atmosphere.py | 19 +++++++++++++++++--
1 file changed, 17 insertions(+), 2 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index 071c7d613b7..a536e0b535e 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -10,6 +10,7 @@
"""
import abc
+import warnings
from dataclasses import dataclass
from functools import partial
@@ -33,6 +34,13 @@
DENSITY_UNIT = u.g / u.cm**3
+class SlantDepthZenithRangeWarning(UserWarning):
+ """
+ Issued when the zenith angle of an event is beyond
+ the range where the slant depth computation is correct (approx. 70 deg)
+ """
+
+
class AtmosphereDensityProfile(abc.ABC):
"""
Base class for models of atmosphere density.
@@ -82,8 +90,9 @@ def slant_depth_from_height(
r"""Line-of-sight integral from height to infinity, along
the direction specified by the zenith angle. This is sometimes called
the *slant depth*. The atmosphere here is assumed to be Cartesian, the
- curvature of the Earth is not taken into account. Inverse of
- height_from_slant_depth function.
+ curvature of the Earth is not taken into account. This approximation breaks down
+ for large zenith angles (>70 deg), in which case this function
+ does not give correct results. Inverse of height_from_slant_depth function.
.. math:: X(h, \Psi) = \int_{h}^{\infty} \rho(h') dh' / \cos{\Psi}
@@ -98,6 +107,12 @@ def slant_depth_from_height(
"""
+ if zenith_angle > 70 * u.deg:
+ warnings.warn(
+ "Zenith angle too high for accurate slant depth",
+ SlantDepthZenithRangeWarning,
+ )
+
return (self.integral(height) / np.cos(zenith_angle)).to(output_units)
def height_from_slant_depth(
From 8653377a3d437058a361d3b98f1c1c63398152a1 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 5 Apr 2024 11:33:44 +0200
Subject: [PATCH 24/26] Include linter changes
---
docs/changes/2422.feature.rst | 4 ++--
src/ctapipe/atmosphere.py | 3 +--
2 files changed, 3 insertions(+), 4 deletions(-)
diff --git a/docs/changes/2422.feature.rst b/docs/changes/2422.feature.rst
index f755f4caf4a..8eb28e7f3fd 100644
--- a/docs/changes/2422.feature.rst
+++ b/docs/changes/2422.feature.rst
@@ -1,2 +1,2 @@
-Implement the overburden-to height a.s.l. transformation fucntion in the atmosphere module
-and test that round-trip returns original value.
\ No newline at end of file
+Implement the overburden-to height a.s.l. transformation function in the atmosphere module
+and test that round-trip returns original value.
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index a536e0b535e..b9f29e4c674 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -378,7 +378,7 @@ def _exponential(h, a, b, c):
def _inv_exponential(T, a, b, c):
- "inverse function for exponetial atmosphere"
+ "inverse function for exponential atmosphere"
return -c * np.log((T - a) / b)
@@ -425,7 +425,6 @@ class FiveLayerAtmosphereDensityProfile(AtmosphereDensityProfile):
"""
def __init__(self, table: Table):
-
self.table = QTable(table, copy=False)
param_a = self.table["a"].to(GRAMMAGE_UNIT)
From 5cd88e552d5a37139b076b2b95a83c28f875929e Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 5 Apr 2024 15:22:42 +0200
Subject: [PATCH 25/26] Allow for array of zenith angles in 70 deg check
---
src/ctapipe/atmosphere.py | 10 ++++++++--
1 file changed, 8 insertions(+), 2 deletions(-)
diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py
index b9f29e4c674..0fcd865054c 100644
--- a/src/ctapipe/atmosphere.py
+++ b/src/ctapipe/atmosphere.py
@@ -107,7 +107,7 @@ def slant_depth_from_height(
"""
- if zenith_angle > 70 * u.deg:
+ if np.any(zenith_angle > 70 * u.deg):
warnings.warn(
"Zenith angle too high for accurate slant depth",
SlantDepthZenithRangeWarning,
@@ -135,6 +135,12 @@ def height_from_slant_depth(
"""
+ if np.any(zenith_angle > 70 * u.deg):
+ warnings.warn(
+ "Zenith angle too high for accurate slant depth",
+ SlantDepthZenithRangeWarning,
+ )
+
return (self.height_from_overburden(slant_depth * np.cos(zenith_angle))).to(
output_units
)
@@ -169,7 +175,7 @@ def peek(self):
axis[1].set_yscale("log")
axis[1].grid(True)
- zenith_angle = np.linspace(0, 80, 20) * u.deg
+ zenith_angle = np.linspace(0, 70, 20) * u.deg
for height in [0, 5, 10, 20] * u.km:
column_density = self.slant_depth_from_height(height, zenith_angle)
axis[2].plot(zenith_angle, column_density, label=f"Height={height}")
From 079968742301801312ad3c369f3913cbba8817b4 Mon Sep 17 00:00:00 2001
From: gschwefer <georg.schwefer@mpi-hd.mpg.de>
Date: Fri, 5 Apr 2024 16:40:39 +0200
Subject: [PATCH 26/26] Implement 5 layer atmosphere as fixture
---
src/ctapipe/tests/test_atmosphere.py | 103 ++++++++++-----------------
1 file changed, 36 insertions(+), 67 deletions(-)
diff --git a/src/ctapipe/tests/test_atmosphere.py b/src/ctapipe/tests/test_atmosphere.py
index 021ba969007..64013336451 100644
--- a/src/ctapipe/tests/test_atmosphere.py
+++ b/src/ctapipe/tests/test_atmosphere.py
@@ -32,6 +32,23 @@ def table_profile():
return get_simtel_profile_from_eventsource()
+@pytest.fixture(scope="session")
+def layer5_profile():
+ """a 5 layer profile for testing"""
+ fit_reference = np.array(
+ [
+ [0.00 * 100000, -140.508, 1178.05, 994186, 0],
+ [9.75 * 100000, -18.4377, 1265.08, 708915, 0],
+ [19.00 * 100000, 0.217565, 1349.22, 636143, 0],
+ [46.00 * 100000, -0.000201796, 703.745, 721128, 0],
+ [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],
+ ]
+ )
+
+ profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
+ return profile_5
+
+
def get_simtel_fivelayer_profile():
"""
get a sample 3-layer profile
@@ -98,7 +115,7 @@ def test_table_model_interpolation(table_profile):
assert np.isfinite(table_profile.integral(height_fine)).all()
-def test_against_reference():
+def test_against_reference(layer5_profile):
"""
Test five-layer and table methods against a reference analysis from
SimTelArray. Data from communication with K. Bernloehr.
@@ -111,25 +128,14 @@ def test_against_reference():
"atmosphere_profile_comparison_from_simtelarray"
)
- fit_reference = np.array(
- [
- [0.00 * 100000, -140.508, 1178.05, 994186, 0],
- [9.75 * 100000, -18.4377, 1265.08, 708915, 0],
- [19.00 * 100000, 0.217565, 1349.22, 636143, 0],
- [46.00 * 100000, -0.000201796, 703.745, 721128, 0],
- [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],
- ]
- )
-
- profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
-
height = reference_table["Altitude_km"].to("km")
np.testing.assert_allclose(
- 1.0 - profile_5(height) / reference_table["rho_5"], 0, atol=1e-5
+ 1.0 - layer5_profile(height) / reference_table["rho_5"], 0, atol=1e-5
)
np.testing.assert_allclose(
- 1.0 - profile_5.slant_depth_from_height(height) / reference_table["thick_5"],
+ 1.0
+ - layer5_profile.slant_depth_from_height(height) / reference_table["thick_5"],
0,
atol=1e-5,
)
@@ -171,30 +177,17 @@ def test_slant_depth_from_height(table_profile):
)
-def test_height_overburden_roundtrip(table_profile):
+def test_height_overburden_roundtrip(table_profile, layer5_profile):
"""
Test that successive application of height to overburden
and overburden to height functions return original values
"""
- # Five-layer atmosphere
- fit_reference = np.array(
- [
- [0.00 * 100000, -140.508, 1178.05, 994186, 0],
- [9.75 * 100000, -18.4377, 1265.08, 708915, 0],
- [19.00 * 100000, 0.217565, 1349.22, 636143, 0],
- [46.00 * 100000, -0.000201796, 703.745, 721128, 0],
- [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],
- ]
- )
-
- profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
-
layer_5_heights = u.Quantity([5, 15, 30, 70, 110] * u.km)
for height in layer_5_heights:
- circle_height_5_layer = profile_5.height_from_overburden(
- profile_5.integral(height)
+ circle_height_5_layer = layer5_profile.height_from_overburden(
+ layer5_profile.integral(height)
)
assert np.allclose(circle_height_5_layer, height, rtol=0.005)
@@ -217,37 +210,24 @@ def test_height_overburden_roundtrip(table_profile):
assert np.allclose(circle_height_table, 47 * u.km, rtol=0.0001)
-def test_height_slant_depth_roundtrip(table_profile):
+def test_height_slant_depth_roundtrip(table_profile, layer5_profile):
"""
Test that successive application of height to overburden
and overburden to height functions return original values
"""
- # Five-layer atmosphere
- fit_reference = np.array(
- [
- [0.00 * 100000, -140.508, 1178.05, 994186, 0],
- [9.75 * 100000, -18.4377, 1265.08, 708915, 0],
- [19.00 * 100000, 0.217565, 1349.22, 636143, 0],
- [46.00 * 100000, -0.000201796, 703.745, 721128, 0],
- [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],
- ]
- )
-
- profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
-
layer_5_heights = u.Quantity([5, 15, 30, 70, 110] * u.km)
for height in layer_5_heights:
- circle_height_5_layer = profile_5.height_from_slant_depth(
- profile_5.slant_depth_from_height(height)
+ circle_height_5_layer = layer5_profile.height_from_slant_depth(
+ layer5_profile.slant_depth_from_height(height)
)
assert np.allclose(circle_height_5_layer, height, rtol=0.005)
for height in layer_5_heights:
- circle_height_5_layer_z_20 = profile_5.height_from_slant_depth(
- profile_5.slant_depth_from_height(height, zenith_angle=20 * u.deg),
+ circle_height_5_layer_z_20 = layer5_profile.height_from_slant_depth(
+ layer5_profile.slant_depth_from_height(height, zenith_angle=20 * u.deg),
zenith_angle=20 * u.deg,
)
@@ -313,26 +293,15 @@ def test_out_of_range_exponential():
assert np.isnan(density_model.integral(-0.1 * u.km))
-def test_out_of_range_five_layer():
+def test_out_of_range_five_layer(layer5_profile):
# Five-layer atmosphere
- fit_reference = np.array(
- [
- [0.00 * 100000, -140.508, 1178.05, 994186, 0],
- [9.75 * 100000, -18.4377, 1265.08, 708915, 0],
- [19.00 * 100000, 0.217565, 1349.22, 636143, 0],
- [46.00 * 100000, -0.000201796, 703.745, 721128, 0],
- [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],
- ]
- )
-
- profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
- assert np.isposinf(profile_5.height_from_overburden(0 * atmo.GRAMMAGE_UNIT))
+ assert np.isposinf(layer5_profile.height_from_overburden(0 * atmo.GRAMMAGE_UNIT))
- assert np.isnan(profile_5.height_from_overburden(2000 * atmo.GRAMMAGE_UNIT))
+ assert np.isnan(layer5_profile.height_from_overburden(2000 * atmo.GRAMMAGE_UNIT))
- assert np.allclose(profile_5(150 * u.km).value, 0.0, atol=1e-9)
- assert np.isnan(profile_5(-0.1 * u.km))
+ assert np.allclose(layer5_profile(150 * u.km).value, 0.0, atol=1e-9)
+ assert np.isnan(layer5_profile(-0.1 * u.km))
- assert np.allclose(profile_5.integral(150 * u.km).value, 0.0, atol=1e-9)
- assert np.isnan(profile_5.integral(-0.1 * u.km))
+ assert np.allclose(layer5_profile.integral(150 * u.km).value, 0.0, atol=1e-9)
+ assert np.isnan(layer5_profile.integral(-0.1 * u.km))