From 2cb299cab62d50430ffdfeb68f4d794bfccdf6a0 Mon Sep 17 00:00:00 2001 From: gschwefer 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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))