[fix/feat] Remove duplicate functions in grid_utils

posydon/unit_tests/utils/test_common_functions.py

@@ -89,7 +89,7 @@ def test_dir(self):
                     'THRESHOLD_HE_NAKED_ABUNDANCE', '__authors__',\
                     '__builtins__', '__cached__', '__doc__', '__file__',\
                     '__loader__', '__name__', '__package__', '__spec__',\
-                    'beaming', 'bondi_hoyle',\
+                    'beaming', 'beta_gw', 'bondi_hoyle',\
                     'calculate_H2recombination_energy',\
                     'calculate_Mejected_for_integrated_binding_energy',\
                     'calculate_Patton20_values_at_He_depl',\
@@ -233,6 +233,9 @@ def test_instance_histogram_sampler(self):
     def test_instance_read_histogram_from_file(self):
         assert isroutine(totest.read_histogram_from_file)
 
+    def test_instance_beta_gw(self):
+        assert isroutine(totest.beta_gw)
+
     def test_instance_inspiral_timescale_from_separation(self):
         assert isroutine(totest.inspiral_timescale_from_separation)
 
@@ -922,6 +925,19 @@ def test_read_histogram_from_file(self, csv_path_failing_3_data_lines,\
         assert np.allclose(arrays[0], np.array([0.2, 1.2, 2.2]))
         assert np.allclose(arrays[1], np.array([2.0, 2.0]))
 
+    def test_beta_gw(self):
+        # missing argument
+        with raises(TypeError, match="missing 2 required positional "\
+                                     +"arguments: 'star1_mass' and "\
+                                     +"'star2_mass'"):
+            totest.beta_gw()
+        # examples
+        tests = [(15.0, 30.0, approx(3.18232660295e-69, abs=6e-81)),\
+                 (30.0, 30.0, approx(8.48620427454e-69, abs=6e-81)),\
+                 (30.0, 60.0, approx(2.54586128236e-68, abs=6e-80))]
+        for (m1, m2, r) in tests:
+            assert totest.beta_gw(m1, m2) == r
+
     def test_inspiral_timescale_from_separation(self):
         # missing argument
         with raises(TypeError, match="missing 4 required positional "\

posydon/unit_tests/utils/test_gridutils.py

@@ -31,16 +31,17 @@ def test_dir(self):
         ## does not clear the warning registy correctly.
         if hasattr(totest, '__warningregistry__'):
             del totest.__warningregistry__
-        elements = {'LG_MTRANSFER_RATE_THRESHOLD', 'Msun', 'Pwarn', 'Rsun',\
-                    'T_merger_P', 'T_merger_a', '__authors__', '__builtins__',\
+        elements = {'LG_MTRANSFER_RATE_THRESHOLD', 'Pwarn',\
+                    '__authors__', '__builtins__',\
                     '__cached__', '__doc__', '__file__', '__loader__',\
                     '__name__', '__package__', '__spec__', 'add_field',\
-                    'beta_gw', 'cgrav', 'clean_inlist_file', 'clight',\
+                    'clean_inlist_file',\
                     'convert_output_to_table', 'find_index_nearest_neighbour',\
                     'find_nearest', 'fix_He_core', 'get_cell_edges',\
                     'get_final_proposed_points', 'get_new_grid_name', 'gzip',\
-                    'join_lists', 'kepler3_a', 'np', 'os', 'pd',\
-                    'read_EEP_data_file', 'read_MESA_data_file', 'secyear'}
+                    'inspiral_timescale_from_orbital_period',\
+                    'join_lists', 'np', 'os', 'pd',\
+                    'read_EEP_data_file', 'read_MESA_data_file'}
         totest_elements = set(dir(totest))
         missing_in_test = elements - totest_elements
         assert len(missing_in_test) == 0, "There are missing objects in "\
@@ -83,18 +84,6 @@ def test_instance_find_index_nearest_neighbour(self):
     def test_instance_get_final_proposed_points(self):
         assert isroutine(totest.get_final_proposed_points)
 
-    def test_instance_T_merger_P(self):
-        assert isroutine(totest.T_merger_P)
-
-    def test_instance_beta_gw(self):
-        assert isroutine(totest.beta_gw)
-
-    def test_instance_kepler3_a(self):
-        assert isroutine(totest.kepler3_a)
-
-    def test_instance_T_merger_a(self):
-        assert isroutine(totest.T_merger_a)
-
     def test_instance_convert_output_to_table(self):
         assert isroutine(totest.convert_output_to_table)
 
@@ -555,57 +544,6 @@ def test_get_final_proposed_points(self, capsys):
         assert np.allclose(mx, np.array([0.1, 0.3, 0.3, 0.5]))
         assert np.allclose(my, np.array([0.3, 0.3, 0.5, 0.5]))
 
-    def test_T_merger_P(self):
-        # missing argument
-        with raises(TypeError, match="missing 3 required positional "\
-                                     +"arguments: 'P', 'm1', and 'm2'"):
-            totest.T_merger_P()
-        # examples
-        tests = [(1.0, 15.0, 30.0, approx(0.37210532488, abs=6e-12)),\
-                 (2.0, 15.0, 30.0, approx(2.36272153666, abs=6e-12)),\
-                 (1.0, 30.0, 30.0, approx(0.20477745195, abs=6e-12)),\
-                 (1.0, 15.0, 60.0, approx(0.22058982311, abs=6e-12))]
-        for (P, m1, m2, r) in tests:
-            assert totest.T_merger_P(P, m1, m2) == r
-
-    def test_beta_gw(self):
-        # missing argument
-        with raises(TypeError, match="missing 2 required positional "\
-                                     +"arguments: 'm1' and 'm2'"):
-            totest.beta_gw()
-        # examples
-        tests = [(15.0, 30.0, approx(3.18232660295e-69, abs=6e-81)),\
-                 (30.0, 30.0, approx(8.48620427454e-69, abs=6e-81)),\
-                 (30.0, 60.0, approx(2.54586128236e-68, abs=6e-80))]
-        for (m1, m2, r) in tests:
-            assert totest.beta_gw(m1, m2) == r
-
-    def test_kepler3_a(self):
-        # missing argument
-        with raises(TypeError, match="missing 3 required positional "\
-                                     +"arguments: 'P', 'm1', and 'm2'"):
-            totest.kepler3_a()
-        # examples
-        tests = [(1.0, 15.0, 30.0, approx(14.9643417735, abs=6e-11)),\
-                 (2.0, 15.0, 30.0, approx(23.7544118733, abs=6e-11)),\
-                 (1.0, 30.0, 30.0, approx(16.4703892879, abs=6e-11)),\
-                 (1.0, 15.0, 60.0, approx(17.7421890201, abs=6e-11))]
-        for (P, m1, m2, r) in tests:
-            assert totest.kepler3_a(P, m1, m2) == r
-
-    def test_T_merger_a(self):
-        # missing argument
-        with raises(TypeError, match="missing 3 required positional "\
-                                     +"arguments: 'a', 'm1', and 'm2'"):
-            totest.T_merger_a()
-        # examples
-        tests = [(1.0, 15.0, 30.0, approx(7.42053829341e-06, abs=6e-18)),\
-                 (2.0, 15.0, 30.0, approx(1.18728612695e-04, abs=6e-16)),\
-                 (1.0, 30.0, 30.0, approx(2.78270186003e-06, abs=6e-18)),\
-                 (1.0, 15.0, 60.0, approx(2.22616148802e-06, abs=6e-18))]
-        for (a, m1, m2, r) in tests:
-            assert totest.T_merger_a(a, m1, m2) == r
-
     def test_convert_output_to_table(self, no_path, out_path,\
                                      MESA_BH_data_tight_orbit,\
                                      MESA_BH_data_wide_orbit,\

posydon/utils/common_functions.py

@@ -872,6 +872,31 @@ def read_histogram_from_file(path):
     return arrays
 
 
+def beta_gw(star1_mass, star2_mass):
+    """Evaluate Peters' beta coefficient (equation 5.9 from Peters 1964).
+
+    Parameters
+    ----------
+    star1_mass : float
+        Mass of the first star in solar masses.
+    star2_mass : float
+        Mass of the second star in solar masses.
+
+    Returns
+    -------
+    float
+        Peters' beta coefficient with masses given in solar units.
+        To obtain the full CGS beta (cm^4 s^-1), multiply the result
+        by ``const.Msun**3``.
+
+    References
+    ----------
+    .. [1] Peters 1964 Phys. Rev. 136, B1224
+
+    """
+    return (64.0 / 5.0) * const.standard_cgrav**3 / const.clight**5 * star1_mass * star2_mass * (star1_mass + star2_mass)
+
+
 def inspiral_timescale_from_separation(star1_mass, star2_mass,
                                        separation, eccentricity):
     """Compute the timescale of GW inspiral using the orbital separation.

posydon/utils/gridutils.py

@@ -6,9 +6,7 @@
 import numpy as np
 import pandas as pd
 
-from posydon.utils.constants import Msun, Rsun, clight
-from posydon.utils.constants import secyer as secyear
-from posydon.utils.constants import standard_cgrav as cgrav
+from posydon.utils.common_functions import inspiral_timescale_from_orbital_period
 from posydon.utils.limits_thresholds import LG_MTRANSFER_RATE_THRESHOLD
 from posydon.utils.posydonwarning import Pwarn
 
@@ -243,89 +241,6 @@ def get_final_proposed_points(proposed_x, grid_x, proposed_y, grid_y):
     return mapped_x, mapped_y
 
 
-def T_merger_P(P, m1, m2):
-    """Merger time given initial orbital period and masses of binary.
-
-    Parameters
-    ----------
-    P : float
-        Orbital period (days).
-    m1 : float
-        Mass of first star (Msun).
-    m2 : float
-        Mass of second star (Msun).
-
-    Returns
-    -------
-    float
-        Merger time (Gyr)
-
-    """
-    return T_merger_a(kepler3_a(P, m1, m2), m1, m2)
-
-
-def beta_gw(m1, m2):
-    """Evaluate the "beta" (equation 5.9) from Peters (1964).
-
-    Parameters
-    ----------
-    m1 : float
-        Mass of the first star.
-    m2 : type
-        Mass of the second star.
-
-    Returns
-    -------
-    float
-        Peters' beta in cgs units.
-
-    """
-    return 64. / 5. * cgrav**3 / clight**5 * m1 * m2 * (m1 + m2)
-
-
-def kepler3_a(P, m1, m2):
-    """Calculate the semimajor axis of a binary from its period and masses.
-
-    Parameters
-    ----------
-    P : float
-        Orbital period (days).
-    m1 : float
-        Mass of first star.
-    m2 : type
-        Mass of second star.
-
-    Returns
-    -------
-    float
-        Semi-major axis (Rsun) using Kepler's third law.
-
-    """
-    return ((P * 24.0 * 3600.0)**2.0
-            * cgrav * (m1 + m2) * Msun / (4.0 * np.pi**2))**(1.0 / 3.0) / Rsun
-
-
-def T_merger_a(a, m1, m2):
-    """Merger time given initial orbital separation and masses of binary.
-
-    Parameters
-    ----------
-    a : float
-        Orbital separation (Rsun).
-    m1 : float
-        Mass of first star (Msun).
-    m2 : float
-        Mass of second star (Msun).
-
-    Returns
-    -------
-    float
-        Merger time (Gyr) following Peters (1964), eq. (5.10).
-
-    """
-    return (a * Rsun)**4 / (4. * beta_gw(m1, m2) * Msun**3) / (secyear * 1.0e9)
-
-
 def convert_output_to_table(
     output_file, binary_history_file=None,
     star1_history_file=None, star2_history_file=None, column_names=[
@@ -443,10 +358,11 @@ def convert_output_to_table(
                 values["C_core_2(Msun)"] = hist["c_core_mass"].iloc[-1]
 
             if binary_history_file is not None:
-                tmerge = T_merger_P(
-                    binary_history["period_days"].iloc[-1],
+                tmerge = inspiral_timescale_from_orbital_period(
                     binary_history["star_1_mass"].iloc[-1],
-                    binary_history["star_2_mass"].iloc[-1])
+                    binary_history["star_2_mass"].iloc[-1],
+                    binary_history["period_days"].iloc[-1],
+                    0.0) / 1000.0
                 max_lg_mtransfer_rate = binary_history[
                     "lg_mtransfer_rate"].max()