diff --git a/pycbc/conversions.py b/pycbc/conversions.py index b86c402e843..f3c1f1d055f 100644 --- a/pycbc/conversions.py +++ b/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',