From d02be9fc0601ed15443652938a17c6c34c5bd8d0 Mon Sep 17 00:00:00 2001 From: rahuldhurkunde Date: Tue, 3 Feb 2026 18:40:44 +0000 Subject: [PATCH 1/3] Added eccentric_mchirp --- pycbc/conversions.py | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/pycbc/conversions.py b/pycbc/conversions.py index b86c402e843..c47a880881f 100644 --- a/pycbc/conversions.py +++ b/pycbc/conversions.py @@ -478,6 +478,47 @@ def mass2_from_tau0_tau3(tau0, tau3, f_lower): return mass2_from_mtotal_eta(mtotal, eta) +def ecc_mchirp_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 + + ecc_mchirp = mchirp * (1 + alpha * e2 + beta * e4 + gamma * e6) + + return ecc_mchirp + def lambda_tilde(mass1, mass2, lambda1, lambda2): """ The effective lambda parameter From 66622daf95092bcebf12c4f08bac510b9622d600 Mon Sep 17 00:00:00 2001 From: rahuldhurkunde Date: Thu, 5 Feb 2026 11:51:41 +0000 Subject: [PATCH 2/3] added the reverse mchirp_from_ecc_mchirp_ecc --- pycbc/conversions.py | 88 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) diff --git a/pycbc/conversions.py b/pycbc/conversions.py index c47a880881f..e7996ad232a 100644 --- a/pycbc/conversions.py +++ b/pycbc/conversions.py @@ -162,6 +162,34 @@ def hypertriangle(*params, bounds=(0, 1)): out_params = [out_params[i][0] for i in range(K)] return out_params +def solve_root(f, a, b, tol=1e-6, max_iter=100): + """Simple bisection root solver.""" + fa = f(a) + fb = f(b) + + if fa * fb > 0: + raise ValueError(f"Root not bracketed: f({a})={fa}, f({b})={fb}") + + if abs(fa) < tol: return a + if abs(fb) < tol: return b + + for _ in range(max_iter): + c = (a + b) / 2 + fc = f(c) + + if abs(fc) < tol or (b - a) / 2 < tol: + return c + + if fa * fc < 0: + b = c + fb = fc + else: + a = c + fa = fc + + return (a + b) / 2 + + # # ============================================================================= # @@ -519,6 +547,66 @@ def ecc_mchirp_from_mchirp_ecc(mchirp, ecc): return ecc_mchirp +def _mchirp_from_ecc_mchirp_ecc(ecc_mchirp, ecc): + """ + Inverse function to get mchirp from ecc_mchirp and ecc. + Solves ecc_mchirp_from_mchirp_ecc(mchirp, ecc) - ecc_mchirp = 0 + """ + def target(mchirp): + return ecc_mchirp_from_mchirp_ecc(mchirp, ecc) - ecc_mchirp + + # Search strategy: + # We know mchirp ~ ecc_mchirp. + # We try to bracket the root starting from this guess. + + guess = ecc_mchirp + # Initial bracket width factor + factor = 0.2 + + # Try finding a bracket + lower = guess * (1 - factor) + upper = guess * (1 + factor) + + # Limit max iterations for bracket expansion + max_iter = 10 + + try: + f_lower = target(lower) + f_upper = target(upper) + + # Expand if signs are the same + for _ in range(max_iter): + if f_lower * f_upper < 0: + break + + if abs(f_lower) < abs(f_upper): + # Root is likely lower + lower /= 1.5 + f_lower = target(lower) + else: + # Root is likely higher + upper *= 1.5 + f_upper = target(upper) + + # Safety checks for physical range (roughly 0 to 100 for this approx) + if upper > 100: + upper = 100 + f_upper = target(upper) + if f_lower * f_upper > 0: + # If still same sign at 100, we might be out of validity range + # or function turns over. + break + + if lower < 0.1: + lower = 0.1 + f_lower = target(lower) + + return solve_root(target, lower, upper) + except Exception: + return numpy.nan + +mchirp_from_ecc_mchirp_ecc = numpy.vectorize(_mchirp_from_ecc_mchirp_ecc) + def lambda_tilde(mass1, mass2, lambda1, lambda2): """ The effective lambda parameter From d2e1e750e7218e74f8c95a52ed24f34e71d9351c Mon Sep 17 00:00:00 2001 From: rahuldhurkunde Date: Wed, 25 Feb 2026 18:39:31 +0000 Subject: [PATCH 3/3] bring back the cubic term, update the inverse function --- pycbc/conversions.py | 136 ++++++++++++++++--------------------------- 1 file changed, 51 insertions(+), 85 deletions(-) diff --git a/pycbc/conversions.py b/pycbc/conversions.py index e7996ad232a..f3c1f1d055f 100644 --- a/pycbc/conversions.py +++ b/pycbc/conversions.py @@ -162,34 +162,6 @@ def hypertriangle(*params, bounds=(0, 1)): out_params = [out_params[i][0] for i in range(K)] return out_params -def solve_root(f, a, b, tol=1e-6, max_iter=100): - """Simple bisection root solver.""" - fa = f(a) - fb = f(b) - - if fa * fb > 0: - raise ValueError(f"Root not bracketed: f({a})={fa}, f({b})={fb}") - - if abs(fa) < tol: return a - if abs(fb) < tol: return b - - for _ in range(max_iter): - c = (a + b) / 2 - fc = f(c) - - if abs(fc) < tol or (b - a) / 2 < tol: - return c - - if fa * fc < 0: - b = c - fb = fc - else: - a = c - fa = fc - - return (a + b) / 2 - - # # ============================================================================= # @@ -506,7 +478,7 @@ def mass2_from_tau0_tau3(tau0, tau3, f_lower): return mass2_from_mtotal_eta(mtotal, eta) -def ecc_mchirp_from_mchirp_ecc(mchirp, ecc): +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 @@ -543,70 +515,63 @@ def ecc_mchirp_from_mchirp_ecc(mchirp, ecc): e4 = e2**2 e6 = e4 * e2 - ecc_mchirp = mchirp * (1 + alpha * e2 + beta * e4 + gamma * e6) + Emchirp = mchirp * (1 + alpha * e2 + beta * e4 + gamma * e6) - return ecc_mchirp + return Emchirp -def _mchirp_from_ecc_mchirp_ecc(ecc_mchirp, ecc): +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. """ - Inverse function to get mchirp from ecc_mchirp and ecc. - Solves ecc_mchirp_from_mchirp_ecc(mchirp, ecc) - ecc_mchirp = 0 - """ - def target(mchirp): - return ecc_mchirp_from_mchirp_ecc(mchirp, ecc) - ecc_mchirp + Emchirp, ecc, input_is_array = ensurearray(Emchirp, ecc) - # Search strategy: - # We know mchirp ~ ecc_mchirp. - # We try to bracket the root starting from this guess. + # 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 - guess = ecc_mchirp - # Initial bracket width factor - factor = 0.2 + e2 = ecc**2 + e4 = e2**2 + e6 = e4 * e2 - # Try finding a bracket - lower = guess * (1 - factor) - upper = guess * (1 + factor) + # 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) - # Limit max iterations for bracket expansion - max_iter = 10 + 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 - try: - f_lower = target(lower) - f_upper = target(upper) - - # Expand if signs are the same - for _ in range(max_iter): - if f_lower * f_upper < 0: - break - - if abs(f_lower) < abs(f_upper): - # Root is likely lower - lower /= 1.5 - f_lower = target(lower) - else: - # Root is likely higher - upper *= 1.5 - f_upper = target(upper) - - # Safety checks for physical range (roughly 0 to 100 for this approx) - if upper > 100: - upper = 100 - f_upper = target(upper) - if f_lower * f_upper > 0: - # If still same sign at 100, we might be out of validity range - # or function turns over. - break - - if lower < 0.1: - lower = 0.1 - f_lower = target(lower) - - return solve_root(target, lower, upper) - except Exception: - return numpy.nan - -mchirp_from_ecc_mchirp_ecc = numpy.vectorize(_mchirp_from_ecc_mchirp_ecc) + 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 @@ -1998,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',