New effective eccentric mchirp parameter

pycbc/conversions.py

@@ -478,6 +478,100 @@ def mass2_from_tau0_tau3(tau0, tau3, f_lower):
     return mass2_from_mtotal_eta(mtotal, eta)
 
 
+def Emchirp_from_mchirp_ecc(mchirp, ecc):
+    """ The effective eccentric mchirp parameter as defined in
+    the equation 8 in https://arxiv.org/pdf/2107.14736.
+    Eccentricity must be defined at the dominant (2,2) mode GW
+    frequency of 10 Hz.
+    """
+    # Constants from Table 1
+    xi = 0.06110974175360381
+    delta = -0.4193723077257345
+
+    Xi_beta = 0.00801015132110059
+    Delta_beta = -2.14807199936756e-5
+    kappa_beta = 1.12702400406416e-8
+    zeta_beta = -1.9753003183066e-12
+
+    Xi_gamma = 0.024204222771565382
+    Delta_gamma = -6.261945897154536e-6
+    kappa_gamma = 1.1175104924576945e-8
+    zeta_gamma = -3.681726165703978e-12
+
+    # Calculate coefficients (Eq 9)
+    alpha = xi * mchirp + delta
+
+    beta = (Xi_beta * mchirp**2 +
+            Delta_beta * mchirp**4 +
+            kappa_beta * mchirp**6 +
+            zeta_beta * mchirp**8)
+
+    gamma = (Xi_gamma * mchirp**2 +
+             Delta_gamma * mchirp**4 +
+             kappa_gamma * mchirp**6 +
+             zeta_gamma * mchirp**8)
+
+    e2 = ecc**2
+    e4 = e2**2
+    e6 = e4 * e2
+
+    Emchirp = mchirp * (1 + alpha * e2 + beta * e4 + gamma * e6)
+
+    return Emchirp
+
+def mchirp_from_Emchirp_ecc(Emchirp, ecc):
+    """ The inverse of Emchirp_from_mchirp_ecc.
+    Returns the chirp mass given the eccentric chirp mass and eccentricity.
+    The eccentricity must be defined at 10 Hz.
+    """
+    Emchirp, ecc, input_is_array = ensurearray(Emchirp, ecc)
+
+    # Constants from Table 1 of https://arxiv.org/pdf/2107.14736
+    xi = 0.06110974175360381
+    delta = -0.4193723077257345
+    Xi_beta = 0.00801015132110059
+    Delta_beta = -2.14807199936756e-5
+    kappa_beta = 1.12702400406416e-8
+    zeta_beta = -1.9753003183066e-12
+    Xi_gamma = 0.024204222771565382
+    Delta_gamma = -6.261945897154536e-6
+    kappa_gamma = 1.1175104924576945e-8
+    zeta_gamma = -3.681726165703978e-12
+
+    e2 = ecc**2
+    e4 = e2**2
+    e6 = e4 * e2
+
+    # Use Newton's method to solve the equation for mchirp
+    # Initial guess using the quadratic approximation
+    A = xi * e2
+    B = 1 + delta * e2
+    C = -Emchirp
+    m = numpy.where(A > 0, (-B + numpy.sqrt(B**2 - 4*A*C)) / (2*A), Emchirp)
+
+    for _ in range(5):
+        m2 = m**2
+        m3 = m2 * m
+        m4 = m2**2
+        m5 = m4 * m
+        m6 = m3**2
+        m7 = m6 * m
+        m8 = m4**2
+
+        alpha = xi * m + delta
+        beta = Xi_beta * m2 + Delta_beta * m4 + kappa_beta * m6 + zeta_beta * m8
+        gamma = Xi_gamma * m2 + Delta_gamma * m4 + kappa_gamma * m6 + zeta_gamma * m8
+
+        f = m * (1 + alpha * e2 + beta * e4 + gamma * e6) - Emchirp
+
+        d_alpha = xi
+        d_beta = 2 * Xi_beta * m + 4 * Delta_beta * m3 + 6 * kappa_beta * m5 + 8 * zeta_beta * m7
+        d_gamma = 2 * Xi_gamma * m + 4 * Delta_gamma * m3 + 6 * kappa_gamma * m5 + 8 * zeta_gamma * m7
+
+        df = (1 + alpha * e2 + beta * e4 + gamma * e6) + m * (d_alpha * e2 + d_beta * e4 + d_gamma * e6)
+        m = m - f / df
+
+    return formatreturn(m, input_is_array)
 def lambda_tilde(mass1, mass2, lambda1, lambda2):
     """ The effective lambda parameter
 
@@ -1869,7 +1963,8 @@ def nltides_gw_phase_diff_isco(f_low, f0, amplitude, n, m1, m2):
     return formatreturn(phi_i - phi_l, input_is_array)
 
 
-__all__ = ['dquadmon_from_lambda', 'lambda_tilde',
+__all__ = ['Emchirp_from_mchirp_ecc', 'mchirp_from_Emchirp_ecc',
+           'dquadmon_from_lambda', 'lambda_tilde',
            'lambda_from_mass_tov_file', 'primary_mass',
            'secondary_mass', 'mtotal_from_mass1_mass2',
            'q_from_mass1_mass2', 'invq_from_mass1_mass2',