Merge pull request #145 from michealroberts/feature/cosmology/get_angular_diameter_distance

feat: add get_angular_diameter_distance to cosmology in celerity module

Author: michealroberts

src/celerity/cosmology.py

@@ -111,8 +111,8 @@ def get_luminosity_distance(z: float, H: Measurement = H0_PLANCK_2018) -> Measur
     Luminosity distance D_L(z) = (1 + z) · D_C(z), in meters.
 
     Args:
-        z:      Redshift (z), must be >= 0.
-        H:      Measurement of the Hubble constant (H0) in km/s/Mpc with its associated uncertainty.
+        z: Redshift (z), must be >= 0.
+        H: Measurement of the Hubble constant (H0) in km/s/Mpc with its associated uncertainty.
 
     Returns:
         Measurement: Luminosity distance in meters with its associated uncertainty.
@@ -136,3 +136,37 @@ def get_luminosity_distance(z: float, H: Measurement = H0_PLANCK_2018) -> Measur
 
 
 # **************************************************************************************
+
+
+def get_angular_diameter_distance(
+    z: float, H: Measurement = H0_PLANCK_2018
+) -> Measurement:
+    """
+    Angular diameter distance D_A(z) = D_C(z) / (1 + z), in meters.
+
+    Args:
+        z: Redshift (z), must be >= 0.
+        H: Measurement of the Hubble constant (H0) in km/s/Mpc with its associated uncertainty.
+
+    Returns:
+        Measurement: Angular diameter distance in meters with its associated uncertainty.
+
+    Raises:
+        ValueError: If z is negative.
+        ValueError: If the Hubble constant is not positive.
+    """
+    if z < 0:
+        raise ValueError("Redshift must be non-negative.")
+
+    # Compute comoving distance (with its propagated uncertainty)
+    d = get_comoving_distance(z, H)
+
+    return Measurement(
+        {
+            "value": d["value"] / (1.0 + z),
+            "uncertainty": d.get("uncertainty", 0.0) / (1.0 + z),
+        }
+    )
+
+
+# **************************************************************************************

tests/test_cosmology.py

@@ -10,6 +10,7 @@
 
 from celerity.constants import H0_PLANCK_2018, H0_SH0ES_2022
 from celerity.cosmology import (
+    get_angular_diameter_distance,
     get_comoving_distance,
     get_hubble_parameter,
     get_luminosity_distance,
@@ -259,6 +260,111 @@ def test_shoes_z1000(self):
         )
 
 
+# **************************************************************************************
+
+
+class TestAngularDiameterDistance(unittest.TestCase):
+    def test_planck_z0(self):
+        d = get_angular_diameter_distance(0.0, H0_PLANCK_2018)
+        self.assertEqual(d["value"], 0.0)
+        self.assertEqual(d["uncertainty"], 0.0)
+
+    def test_shoes_z0(self):
+        d = get_angular_diameter_distance(0.0, H0_SH0ES_2022)
+        self.assertEqual(d["value"], 0.0)
+        self.assertEqual(d["uncertainty"], 0.0)
+
+    def test_planck_z0_1(self):
+        d = get_angular_diameter_distance(0.1, H0_PLANCK_2018)
+        expected_val = 1.21174e25
+        expected_unc = 8.6e22
+        self.assertAlmostEqual(d["value"], expected_val, delta=0.05 * expected_val)
+        self.assertAlmostEqual(
+            d["uncertainty"], expected_unc, delta=0.20 * expected_unc
+        )
+
+    def test_planck_z1(self):
+        d = get_angular_diameter_distance(1.0, H0_PLANCK_2018)
+        expected_val = 5.221447e25
+        expected_unc = 3.5475e23
+        self.assertAlmostEqual(d["value"], expected_val, delta=0.05 * expected_val)
+        self.assertAlmostEqual(
+            d["uncertainty"], expected_unc, delta=0.20 * expected_unc
+        )
+
+    def test_planck_z10(self):
+        d = get_angular_diameter_distance(10.0, H0_PLANCK_2018)
+        expected_val = 2.688876e25
+        expected_unc = 1.8232e23
+        self.assertAlmostEqual(d["value"], expected_val, delta=0.05 * expected_val)
+        self.assertAlmostEqual(
+            d["uncertainty"], expected_unc, delta=0.20 * expected_unc
+        )
+
+    def test_planck_z100(self):
+        d = get_angular_diameter_distance(100.0, H0_PLANCK_2018)
+        expected_val = 3.901548e24
+        expected_unc = 2.65067e22
+        self.assertAlmostEqual(d["value"], expected_val, delta=0.05 * expected_val)
+        self.assertAlmostEqual(
+            d["uncertainty"], expected_unc, delta=0.20 * expected_unc
+        )
+
+    def test_planck_z1000(self):
+        d = get_angular_diameter_distance(1000.0, H0_PLANCK_2018)
+        expected_val = 4.515514e23
+        expected_unc = 3.06893e21
+        self.assertAlmostEqual(d["value"], expected_val, delta=0.05 * expected_val)
+        self.assertAlmostEqual(
+            d["uncertainty"], expected_unc, delta=0.20 * expected_unc
+        )
+
+    def test_shoes_z0_1(self):
+        d = get_angular_diameter_distance(0.1, H0_SH0ES_2022)
+        expected_val = 1.12402e25
+        expected_unc = 1.634e23
+        self.assertAlmostEqual(d["value"], expected_val, delta=0.10 * expected_val)
+        self.assertAlmostEqual(
+            d["uncertainty"], expected_unc, delta=0.25 * expected_unc
+        )
+
+    def test_shoes_z1(self):
+        d = get_angular_diameter_distance(1.0, H0_SH0ES_2022)
+        expected_val = 4.842574e25
+        expected_unc = 6.9058e23
+        self.assertAlmostEqual(d["value"], expected_val, delta=0.10 * expected_val)
+        self.assertAlmostEqual(
+            d["uncertainty"], expected_unc, delta=0.25 * expected_unc
+        )
+
+    def test_shoes_z10(self):
+        d = get_angular_diameter_distance(10.0, H0_SH0ES_2022)
+        expected_val = 2.493742e25
+        expected_unc = 3.5518e23
+        self.assertAlmostEqual(d["value"], expected_val, delta=0.10 * expected_val)
+        self.assertAlmostEqual(
+            d["uncertainty"], expected_unc, delta=0.25 * expected_unc
+        )
+
+    def test_shoes_z100(self):
+        d = get_angular_diameter_distance(100.0, H0_SH0ES_2022)
+        expected_val = 3.618497e24
+        expected_unc = 5.15149e22
+        self.assertAlmostEqual(d["value"], expected_val, delta=0.10 * expected_val)
+        self.assertAlmostEqual(
+            d["uncertainty"], expected_unc, delta=0.25 * expected_unc
+        )
+
+    def test_shoes_z1000(self):
+        d = get_angular_diameter_distance(1000.0, H0_SH0ES_2022)
+        expected_val = 3.957134e23
+        expected_unc = 5.63253e21
+        self.assertAlmostEqual(d["value"], expected_val, delta=0.10 * expected_val)
+        self.assertAlmostEqual(
+            d["uncertainty"], expected_unc, delta=0.25 * expected_unc
+        )
+
+
 # **************************************************************************************
 
 if __name__ == "__main__":