From 68093e5677f1df5f2e7dd4385815c24f182c34de Mon Sep 17 00:00:00 2001 From: Bhooshan Gadre Date: Mon, 26 Jun 2023 04:55:30 -0700 Subject: [PATCH] removing assertion --- pycbc/pnutils.py | 426 ++++++++++++++++++++++++++++++----------------- 1 file changed, 271 insertions(+), 155 deletions(-) diff --git a/pycbc/pnutils.py b/pycbc/pnutils.py index 02cfa87c445..2a88e466204 100644 --- a/pycbc/pnutils.py +++ b/pycbc/pnutils.py @@ -15,7 +15,6 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - # # ============================================================================= # @@ -33,28 +32,34 @@ lalsimulation = libutils.import_optional('lalsimulation') + def nearest_larger_binary_number(input_len): """ Return the nearest binary number larger than input_len. """ return int(2**numpy.ceil(numpy.log2(input_len))) + def chirp_distance(dist, mchirp, ref_mass=1.4): return conversions.chirp_distance(dist, mchirp, ref_mass=ref_mass) + def mass1_mass2_to_mtotal_eta(mass1, mass2): m_total = conversions.mtotal_from_mass1_mass2(mass1, mass2) eta = conversions.eta_from_mass1_mass2(mass1, mass2) - return m_total,eta + return m_total, eta + def mtotal_eta_to_mass1_mass2(m_total, eta): mass1 = conversions.mass1_from_mtotal_eta(m_total, eta) mass2 = conversions.mass2_from_mtotal_eta(m_total, eta) - return mass1,mass2 + return mass1, mass2 + def mass1_mass2_to_mchirp_eta(mass1, mass2): m_chirp = conversions.mchirp_from_mass1_mass2(mass1, mass2) eta = conversions.eta_from_mass1_mass2(mass1, mass2) - return m_chirp,eta + return m_chirp, eta + def mchirp_eta_to_mass1_mass2(m_chirp, eta): mtotal = conversions.mtotal_from_mchirp_eta(m_chirp, eta) @@ -62,6 +67,7 @@ def mchirp_eta_to_mass1_mass2(m_chirp, eta): mass2 = conversions.mass2_from_mtotal_eta(mtotal, eta) return mass1, mass2 + def mchirp_mass1_to_mass2(mchirp, mass1): """ This function takes a value of mchirp and one component mass and returns @@ -79,6 +85,7 @@ def mchirp_mass1_to_mass2(mchirp, mass1): """ return conversions.mass2_from_mchirp_mass1(mchirp, mass1) + def eta_mass1_to_mass2(eta, mass1, return_mass_heavier=False, force_real=True): """ This function takes values for eta and one component mass and returns the @@ -92,8 +99,12 @@ def eta_mass1_to_mass2(eta, mass1, return_mass_heavier=False, force_real=True): mass1 > mass2 is returned. Use the return_mass_heavier kwarg to invert this behaviour. """ - return conversions.mass_from_knownmass_eta(mass1, eta, - known_is_secondary=return_mass_heavier, force_real=force_real) + return conversions.mass_from_knownmass_eta( + mass1, + eta, + known_is_secondary=return_mass_heavier, + force_real=force_real) + def mchirp_q_to_mass1_mass2(mchirp, q): """ This function takes a value of mchirp and the mass ratio @@ -110,32 +121,38 @@ def mchirp_q_to_mass1_mass2(mchirp, q): mass2 = conversions.mass2_from_mchirp_eta(mchirp, eta) return mass1, mass2 + def A0(f_lower): """used in calculating chirp times: see Cokelaer, arxiv.org:0706.4437 appendix 1, also lalinspiral/python/sbank/tau0tau3.py """ return conversions._a0(f_lower) + def A3(f_lower): """another parameter used for chirp times""" return conversions._a3(f_lower) + def mass1_mass2_to_tau0_tau3(mass1, mass2, f_lower): tau0 = conversions.tau0_from_mass1_mass2(mass1, mass2, f_lower) tau3 = conversions.tau3_from_mass1_mass2(mass1, mass2, f_lower) - return tau0,tau3 + return tau0, tau3 + def tau0_tau3_to_mtotal_eta(tau0, tau3, f_lower): mtotal = conversions.mtotal_from_tau0_tau3(tau0, tau3, f_lower) eta = conversions.eta_from_tau0_tau3(tau0, tau3, f_lower) return mtotal, eta + def tau0_tau3_to_mass1_mass2(tau0, tau3, f_lower): - m_total,eta = tau0_tau3_to_mtotal_eta(tau0, tau3, f_lower) + m_total, eta = tau0_tau3_to_mtotal_eta(tau0, tau3, f_lower) return mtotal_eta_to_mass1_mass2(m_total, eta) -def mass1_mass2_spin1z_spin2z_to_beta_sigma_gamma(mass1, mass2, - spin1z, spin2z): + +def mass1_mass2_spin1z_spin2z_to_beta_sigma_gamma(mass1, mass2, spin1z, + spin2z): _, eta = mass1_mass2_to_mtotal_eta(mass1, mass2) # get_beta_sigma_from_aligned_spins() takes # the spin of the heaviest body first @@ -145,6 +162,7 @@ def mass1_mass2_spin1z_spin2z_to_beta_sigma_gamma(mass1, mass2, eta, heavy_spin, light_spin) return beta, sigma, gamma + def get_beta_sigma_from_aligned_spins(eta, spin1z, spin2z): """ Calculate the various PN spin combinations from the masses and spins. @@ -172,7 +190,7 @@ def get_beta_sigma_from_aligned_spins(eta, spin1z, spin2z): """ chiS = 0.5 * (spin1z + spin2z) chiA = 0.5 * (spin1z - spin2z) - delta = (1 - 4 * eta) ** 0.5 + delta = (1 - 4 * eta)**0.5 spinspin = spin1z * spin2z beta = (113. / 12. - 19. / 3. * eta) * chiS beta += 113. / 12. * delta * chiA @@ -184,21 +202,27 @@ def get_beta_sigma_from_aligned_spins(eta, spin1z, spin2z): gamma += (732985. / 2268. + 140. / 9. * eta) * delta * chiA return beta, sigma, gamma + def solar_mass_to_kg(solar_masses): return solar_masses * lal.MSUN_SI + def parsecs_to_meters(distance): return distance * lal.PC_SI + def megaparsecs_to_meters(distance): return parsecs_to_meters(distance) * 1e6 + def velocity_to_frequency(v, M): return conversions.velocity_to_frequency(v, M) + def frequency_to_velocity(f, M): return conversions.frequency_to_velocity(f, M) + def f_SchwarzISCO(M): """ Innermost stable circular orbit (ISCO) for a test particle @@ -216,6 +240,7 @@ def f_SchwarzISCO(M): """ return conversions.f_schwarzchild_isco(M) + def f_BKLISCO(m1, m2): """ Mass ratio dependent ISCO derived from estimates of the final spin @@ -235,8 +260,10 @@ def f_BKLISCO(m1, m2): Frequency in Hz """ # q is defined to be in [0,1] for this formula - q = numpy.minimum(m1/m2, m2/m1) - return f_SchwarzISCO(m1+m2) * ( 1 + 2.8*q - 2.6*q*q + 0.8*q*q*q ) + q = numpy.minimum(m1 / m2, m2 / m1) + return f_SchwarzISCO(m1 + m2) * (1 + 2.8 * q - 2.6 * q * q + + 0.8 * q * q * q) + def f_LightRing(M): """ @@ -255,6 +282,7 @@ def f_LightRing(M): """ return 1.0 / (3.0**(1.5) * lal.PI * M * lal.MTSUN_SI) + def f_ERD(M): """ Effective RingDown frequency studied in Pan et al. (arXiv:0704.1964) @@ -272,7 +300,8 @@ def f_ERD(M): f : float or numpy.array Frequency in Hz """ - return 1.07 * 0.5326 / (2*lal.PI * 0.955 * M * lal.MTSUN_SI) + return 1.07 * 0.5326 / (2 * lal.PI * 0.955 * M * lal.MTSUN_SI) + def f_FRD(m1, m2): """ @@ -294,9 +323,10 @@ def f_FRD(m1, m2): Frequency in Hz """ m_total, eta = mass1_mass2_to_mtotal_eta(m1, m2) - tmp = ( (1. - 0.63*(1. - 3.4641016*eta + 2.9*eta**2)**(0.3)) / - (1. - 0.057191*eta - 0.498*eta**2) ) - return tmp / (2.*lal.PI * m_total*lal.MTSUN_SI) + tmp = ((1. - 0.63 * (1. - 3.4641016 * eta + 2.9 * eta**2)**(0.3)) / + (1. - 0.057191 * eta - 0.498 * eta**2)) + return tmp / (2. * lal.PI * m_total * lal.MTSUN_SI) + def f_LRD(m1, m2): """ @@ -317,6 +347,7 @@ def f_LRD(m1, m2): """ return 1.2 * f_FRD(m1, m2) + def _get_freq(freqfunc, m1, m2, s1z, s2z): """ Wrapper of the LALSimulation function returning the frequency @@ -343,12 +374,15 @@ def _get_freq(freqfunc, m1, m2, s1z, s2z): # Convert to SI units for lalsimulation m1kg = float(m1) * lal.MSUN_SI m2kg = float(m2) * lal.MSUN_SI - return lalsimulation.SimInspiralGetFrequency( - m1kg, m2kg, 0, 0, float(s1z), 0, 0, float(s2z), int(freqfunc)) + return lalsimulation.SimInspiralGetFrequency(m1kg, m2kg, 0, 0, + float(s1z), 0, 0, float(s2z), + int(freqfunc)) + # vectorize to enable calls with numpy arrays _vec_get_freq = numpy.vectorize(_get_freq) + def get_freq(freqfunc, m1, m2, s1z, s2z): """ Returns the LALSimulation function which evaluates the frequency @@ -375,6 +409,7 @@ def get_freq(freqfunc, m1, m2, s1z, s2z): lalsim_ffunc = getattr(lalsimulation, freqfunc) return _vec_get_freq(lalsim_ffunc, m1, m2, s1z, s2z) + def _get_final_freq(approx, m1, m2, s1z, s2z): """ Wrapper of the LALSimulation function returning the final (highest) @@ -401,12 +436,14 @@ def _get_final_freq(approx, m1, m2, s1z, s2z): # Convert to SI units for lalsimulation m1kg = float(m1) * lal.MSUN_SI m2kg = float(m2) * lal.MSUN_SI - return lalsimulation.SimInspiralGetFinalFreq( - m1kg, m2kg, 0, 0, float(s1z), 0, 0, float(s2z), int(approx)) + return lalsimulation.SimInspiralGetFinalFreq(m1kg, m2kg, 0, 0, float(s1z), + 0, 0, float(s2z), int(approx)) + # vectorize to enable calls with numpy arrays _vec_get_final_freq = numpy.vectorize(_get_final_freq) + def get_final_freq(approx, m1, m2, s1z, s2z): """ Returns the LALSimulation function which evaluates the final @@ -435,50 +472,68 @@ def get_final_freq(approx, m1, m2, s1z, s2z): lalsim_approx = lalsimulation.GetApproximantFromString(approx) return _vec_get_final_freq(lalsim_approx, m1, m2, s1z, s2z) + # Dictionary of functions with uniform API taking a # parameter dict indexed on mass1, mass2, spin1z, spin2z named_frequency_cutoffs = { # functions depending on the total mass alone - "SchwarzISCO": lambda p: f_SchwarzISCO(p["mass1"]+p["mass2"]), - "LightRing" : lambda p: f_LightRing(p["mass1"]+p["mass2"]), - "ERD" : lambda p: f_ERD(p["mass1"]+p["mass2"]), + "SchwarzISCO": + lambda p: f_SchwarzISCO(p["mass1"] + p["mass2"]), + "LightRing": + lambda p: f_LightRing(p["mass1"] + p["mass2"]), + "ERD": + lambda p: f_ERD(p["mass1"] + p["mass2"]), # functions depending on the 2 component masses - "BKLISCO" : lambda p: f_BKLISCO(p["mass1"], p["mass2"]), - "FRD" : lambda p: f_FRD(p["mass1"], p["mass2"]), - "LRD" : lambda p: f_LRD(p["mass1"], p["mass2"]), + "BKLISCO": + lambda p: f_BKLISCO(p["mass1"], p["mass2"]), + "FRD": + lambda p: f_FRD(p["mass1"], p["mass2"]), + "LRD": + lambda p: f_LRD(p["mass1"], p["mass2"]), # functions depending on 2 component masses and aligned spins - "MECO" : lambda p: meco_frequency(p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "HybridMECO" : lambda p: hybrid_meco_frequency( + "MECO": + lambda p: meco_frequency(p["mass1"], p["mass2"], p["spin1z"], p["spin2z"]), + "HybridMECO": + lambda p: hybrid_meco_frequency( p["mass1"], p["mass2"], p["spin1z"], p["spin2z"], qm1=None, qm2=None), - "IMRPhenomBFinal": lambda p: get_freq("fIMRPhenomBFinal", - p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "IMRPhenomCFinal": lambda p: get_freq("fIMRPhenomCFinal", - p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "IMRPhenomDPeak": lambda p: get_freq("fIMRPhenomDPeak", - p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "EOBNRv2RD" : lambda p: get_freq("fEOBNRv2RD", p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "EOBNRv2HMRD" : lambda p: get_freq("fEOBNRv2HMRD", p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "SEOBNRv1RD" : lambda p: get_freq("fSEOBNRv1RD", p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "SEOBNRv1Peak": lambda p: get_freq("fSEOBNRv1Peak", p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "SEOBNRv2RD": lambda p: get_freq("fSEOBNRv2RD", p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "SEOBNRv2Peak": lambda p: get_freq("fSEOBNRv2Peak", p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "SEOBNRv4RD": lambda p: get_freq("fSEOBNRv4RD", p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "SEOBNRv4Peak": lambda p: get_freq("fSEOBNRv4Peak", p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]), - "SEOBNRv5Peak": lambda p: get_freq("fSEOBNRv5Peak", p["mass1"], p["mass2"], - p["spin1z"], p["spin2z"]) - } + "IMRPhenomBFinal": + lambda p: get_freq("fIMRPhenomBFinal", p["mass1"], p["mass2"], p["spin1z"], + p["spin2z"]), + "IMRPhenomCFinal": + lambda p: get_freq("fIMRPhenomCFinal", p["mass1"], p["mass2"], p["spin1z"], + p["spin2z"]), + "IMRPhenomDPeak": + lambda p: get_freq("fIMRPhenomDPeak", p["mass1"], p["mass2"], p["spin1z"], + p["spin2z"]), + "EOBNRv2RD": + lambda p: get_freq("fEOBNRv2RD", p["mass1"], p["mass2"], p["spin1z"], p[ + "spin2z"]), + "EOBNRv2HMRD": + lambda p: get_freq("fEOBNRv2HMRD", p["mass1"], p["mass2"], p["spin1z"], p[ + "spin2z"]), + "SEOBNRv1RD": + lambda p: get_freq("fSEOBNRv1RD", p["mass1"], p["mass2"], p["spin1z"], p[ + "spin2z"]), + "SEOBNRv1Peak": + lambda p: get_freq("fSEOBNRv1Peak", p["mass1"], p["mass2"], p["spin1z"], p[ + "spin2z"]), + "SEOBNRv2RD": + lambda p: get_freq("fSEOBNRv2RD", p["mass1"], p["mass2"], p["spin1z"], p[ + "spin2z"]), + "SEOBNRv2Peak": + lambda p: get_freq("fSEOBNRv2Peak", p["mass1"], p["mass2"], p["spin1z"], p[ + "spin2z"]), + "SEOBNRv4RD": + lambda p: get_freq("fSEOBNRv4RD", p["mass1"], p["mass2"], p["spin1z"], p[ + "spin2z"]), + "SEOBNRv4Peak": + lambda p: get_freq("fSEOBNRv4Peak", p["mass1"], p["mass2"], p["spin1z"], p[ + "spin2z"]), + "SEOBNRv5Peak": + lambda p: get_freq("fSEOBNRv5Peak", p["mass1"], p["mass2"], p["spin1z"], p[ + "spin2z"]) +} + def frequency_cutoff_from_name(name, m1, m2, s1z, s2z): """ @@ -503,9 +558,10 @@ def frequency_cutoff_from_name(name, m1, m2, s1z, s2z): f : float or numpy.array Frequency in Hz """ - params = {"mass1":m1, "mass2":m2, "spin1z":s1z, "spin2z":s2z} + params = {"mass1": m1, "mass2": m2, "spin1z": s1z, "spin2z": s2z} return named_frequency_cutoffs[name](params) + def _get_imr_duration(m1, m2, s1z, s2z, f_low, approximant="SEOBNRv4"): """Wrapper of lalsimulation template duration approximate formula""" m1, m2, s1z, s2z, f_low = float(m1), float(m2), float(s1z), float(s2z),\ @@ -513,44 +569,52 @@ def _get_imr_duration(m1, m2, s1z, s2z, f_low, approximant="SEOBNRv4"): if approximant == "SEOBNRv2": chi = lalsimulation.SimIMRPhenomBComputeChi(m1, m2, s1z, s2z) time_length = lalsimulation.SimIMRSEOBNRv2ChirpTimeSingleSpin( - m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, chi, f_low) + m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, chi, f_low) elif approximant == 'IMRPhenomXAS': time_length = lalsimulation.SimIMRPhenomXASDuration( - m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, s1z, s2z, f_low) + m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, s1z, s2z, f_low) elif approximant == "IMRPhenomD": time_length = lalsimulation.SimIMRPhenomDChirpTime( - m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, s1z, s2z, f_low) + m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, s1z, s2z, f_low) elif approximant == "SEOBNRv4": # NB for no clear reason this function has f_low as first argument time_length = lalsimulation.SimIMRSEOBNRv4ROMTimeOfFrequency( - f_low, m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, s1z, s2z) + f_low, m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, s1z, s2z) elif approximant == 'SPAtmplt' or approximant == 'TaylorF2': chi = lalsimulation.SimInspiralTaylorF2ReducedSpinComputeChi( - m1, m2, s1z, s2z - ) + m1, m2, s1z, s2z) time_length = lalsimulation.SimInspiralTaylorF2ReducedSpinChirpTime( - f_low, m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, chi, -1 - ) + f_low, m1 * lal.MSUN_SI, m2 * lal.MSUN_SI, chi, -1) else: raise RuntimeError("I can't calculate a duration for %s" % approximant) - assert time_length > 0, f"Got a negative IMR duration of {time_length} for {approximant} with params (m1, m2, s1z, s2z, f_low) = ({m1}, {m2}, {s1z}, {s2z}, {f_low})" # FIXME Add an extra factor of 1.1 for 'safety' since the duration # functions are approximate return time_length * 1.1 + get_imr_duration = numpy.vectorize(_get_imr_duration) -def get_inspiral_tf(tc, mass1, mass2, spin1, spin2, f_low, n_points=100, - pn_2order=7, approximant='TaylorF2'): + +def get_inspiral_tf(tc, + mass1, + mass2, + spin1, + spin2, + f_low, + n_points=100, + pn_2order=7, + approximant='TaylorF2'): """Compute the time-frequency evolution of an inspiral signal. Return a tuple of time and frequency vectors tracking the evolution of an inspiral signal in the time-frequency plane. """ + # handle param-dependent approximant specification class Params: pass + params = Params() params.mass1 = mass1 params.mass2 = mass2 @@ -569,27 +633,32 @@ class Params: f_high = f_SchwarzISCO(mass1 + mass2) track_f = numpy.logspace(numpy.log10(f_low), numpy.log10(f_high), n_points) - track_t = numpy.array([findchirp_chirptime(float(mass1), float(mass2), - float(f), pn_2order) for f in track_f]) - elif approximant in ['SEOBNRv2', 'SEOBNRv2_ROM_DoubleSpin', - 'SEOBNRv2_ROM_DoubleSpin_HI']: + track_t = numpy.array([ + findchirp_chirptime(float(mass1), float(mass2), float(f), + pn_2order) for f in track_f + ]) + elif approximant in [ + 'SEOBNRv2', 'SEOBNRv2_ROM_DoubleSpin', 'SEOBNRv2_ROM_DoubleSpin_HI' + ]: f_high = get_final_freq('SEOBNRv2', mass1, mass2, spin1, spin2) track_f = numpy.logspace(numpy.log10(f_low), numpy.log10(f_high), n_points) # use HI function as it has wider freq range validity track_t = numpy.array([ - lalsimulation.SimIMRSEOBNRv2ROMDoubleSpinHITimeOfFrequency(f, - solar_mass_to_kg(mass1), solar_mass_to_kg(mass2), - float(spin1), float(spin2)) for f in track_f]) + lalsimulation.SimIMRSEOBNRv2ROMDoubleSpinHITimeOfFrequency( + f, solar_mass_to_kg(mass1), solar_mass_to_kg(mass2), + float(spin1), float(spin2)) for f in track_f + ]) elif approximant in ['SEOBNRv4', 'SEOBNRv4_ROM']: f_high = get_final_freq('SEOBNRv4', mass1, mass2, spin1, spin2) # use frequency below final freq in case of rounding error - track_f = numpy.logspace(numpy.log10(f_low), numpy.log10(0.999*f_high), - n_points) + track_f = numpy.logspace(numpy.log10(f_low), + numpy.log10(0.999 * f_high), n_points) track_t = numpy.array([ - lalsimulation.SimIMRSEOBNRv4ROMTimeOfFrequency( - f, solar_mass_to_kg(mass1), solar_mass_to_kg(mass2), - float(spin1), float(spin2)) for f in track_f]) + lalsimulation.SimIMRSEOBNRv4ROMTimeOfFrequency( + f, solar_mass_to_kg(mass1), solar_mass_to_kg(mass2), + float(spin1), float(spin2)) for f in track_f + ]) else: raise ValueError('Approximant ' + approximant + ' not supported') return (tc - track_t, track_f) @@ -602,30 +671,31 @@ def _energy_coeffs(m1, m2, chi1, chi2): """ Return the center-of-mass energy coefficients up to 3.0pN (2.5pN spin) """ mtot = m1 + m2 - eta = m1*m2 / (mtot*mtot) - chi = (m1*chi1 + m2*chi2) / mtot + eta = m1 * m2 / (mtot * mtot) + chi = (m1 * chi1 + m2 * chi2) / mtot chisym = (chi1 + chi2) / 2. - beta = (113.*chi - 76.*eta*chisym)/12. - sigma12 = 79.*eta*chi1*chi2/8. + beta = (113. * chi - 76. * eta * chisym) / 12. + sigma12 = 79. * eta * chi1 * chi2 / 8. sigmaqm = 81.*m1*m1*chi1*chi1/(16.*mtot*mtot) \ + 81.*m2*m2*chi2*chi2/(16.*mtot*mtot) - energy0 = -0.5*eta - energy2 = -0.75 - eta/12. + energy0 = -0.5 * eta + energy2 = -0.75 - eta / 12. energy3 = 0. - energy4 = -3.375 + (19*eta)/8. - pow(eta,2)/24. + energy4 = -3.375 + (19 * eta) / 8. - pow(eta, 2) / 24. energy5 = 0. energy6 = -10.546875 - (155*pow(eta,2))/96. - (35*pow(eta,3))/5184. \ + eta*(59.80034722222222 - (205*pow(lal.PI,2))/96.) - energy3 += (32*beta)/113. + (52*chisym*eta)/113. + energy3 += (32 * beta) / 113. + (52 * chisym * eta) / 113. - energy4 += (-16*sigma12)/79. - (16*sigmaqm)/81. + energy4 += (-16 * sigma12) / 79. - (16 * sigmaqm) / 81. energy5 += (96*beta)/113. + ((-124*beta)/339. - (522*chisym)/113.)*eta \ - (710*chisym*pow(eta,2))/339. return (energy0, energy2, energy3, energy4, energy5, energy6) + def meco_velocity(m1, m2, chi1, chi2): """ Returns the velocity of the minimum energy cutoff for 3.5pN (2.5pN spin) @@ -648,44 +718,55 @@ def meco_velocity(m1, m2, chi1, chi2): """ _, energy2, energy3, energy4, energy5, energy6 = \ _energy_coeffs(m1, m2, chi1, chi2) + def eprime(v): return 2. + v * v * (4.*energy2 + v * (5.*energy3 \ + v * (6.*energy4 + v * (7.*energy5 + 8.*energy6 * v)))) + return bisect(eprime, 0.05, 1.0) + def _meco_frequency(m1, m2, chi1, chi2): """Returns the frequency of the minimum energy cutoff for 3.5pN (2.5pN spin) """ - return velocity_to_frequency(meco_velocity(m1, m2, chi1, chi2), m1+m2) + return velocity_to_frequency(meco_velocity(m1, m2, chi1, chi2), m1 + m2) + meco_frequency = numpy.vectorize(_meco_frequency) + def _dtdv_coeffs(m1, m2, chi1, chi2): """ Returns the dt/dv coefficients up to 3.5pN (2.5pN spin) """ mtot = m1 + m2 - eta = m1*m2 / (mtot*mtot) - chi = (m1*chi1 + m2*chi2) / mtot + eta = m1 * m2 / (mtot * mtot) + chi = (m1 * chi1 + m2 * chi2) / mtot chisym = (chi1 + chi2) / 2. - beta = (113.*chi - 76.*eta*chisym)/12. - sigma12 = 79.*eta*chi1*chi2/8. + beta = (113. * chi - 76. * eta * chisym) / 12. + sigma12 = 79. * eta * chi1 * chi2 / 8. sigmaqm = 81.*m1*m1*chi1*chi1/(16.*mtot*mtot) \ + 81.*m2*m2*chi2*chi2/(16.*mtot*mtot) - dtdv0 = 1. # FIXME: Wrong but doesn't matter for now. - dtdv2 = (1./336.) * (743. + 924.*eta) + dtdv0 = 1. # FIXME: Wrong but doesn't matter for now. + dtdv2 = (1. / 336.) * (743. + 924. * eta) dtdv3 = -4. * lal.PI + beta - dtdv4 = (3058673. + 5472432.*eta + 4353552.*eta*eta)/1016064. - sigma12 - sigmaqm - dtdv5 = (1./672.) * lal.PI * (-7729. + 1092.*eta) + (146597.*beta/18984. + 42.*beta*eta/113. - 417307.*chisym*eta/18984. - 1389.*chisym*eta*eta/226.) - dtdv6 = 22.065 + 165.416*eta - 2.20067*eta*eta + 4.93152*eta*eta*eta - dtdv6log = 1712./315. - dtdv7 = (lal.PI/1016064.) * (-15419335. - 12718104.*eta + 4975824.*eta*eta) + dtdv4 = (3058673. + 5472432. * eta + + 4353552. * eta * eta) / 1016064. - sigma12 - sigmaqm + dtdv5 = (1. / 672.) * lal.PI * (-7729. + 1092. * eta) + ( + 146597. * beta / 18984. + 42. * beta * eta / 113. - + 417307. * chisym * eta / 18984. - 1389. * chisym * eta * eta / 226.) + dtdv6 = 22.065 + 165.416 * eta - 2.20067 * eta * eta + 4.93152 * eta * eta * eta + dtdv6log = 1712. / 315. + dtdv7 = (lal.PI / 1016064.) * (-15419335. - 12718104. * eta + + 4975824. * eta * eta) return (dtdv0, dtdv2, dtdv3, dtdv4, dtdv5, dtdv6, dtdv6log, dtdv7) + def _dtdv_cutoff_velocity(m1, m2, chi1, chi2): - _, dtdv2, dtdv3, dtdv4, dtdv5, dtdv6, dtdv6log, dtdv7 = _dtdv_coeffs(m1, m2, chi1, chi2) + _, dtdv2, dtdv3, dtdv4, dtdv5, dtdv6, dtdv6log, dtdv7 = _dtdv_coeffs( + m1, m2, chi1, chi2) def dtdv_func(v): x = dtdv7 @@ -701,18 +782,21 @@ def dtdv_func(v): else: return 1.0 + def energy_coefficients(m1, m2, s1z=0, s2z=0, phase_order=-1, spin_order=-1): """ Return the energy coefficients. This assumes that the system has aligned spins only. """ implemented_phase_order = 7 implemented_spin_order = 7 if phase_order > implemented_phase_order: - raise ValueError("pN coeffiecients of that order have not been implemented") + raise ValueError( + "pN coeffiecients of that order have not been implemented") elif phase_order == -1: phase_order = implemented_phase_order if spin_order > implemented_spin_order: - raise ValueError("pN coeffiecients of that order have not been implemented") + raise ValueError( + "pN coeffiecients of that order have not been implemented") elif spin_order == -1: spin_order = implemented_spin_order @@ -720,7 +804,7 @@ def energy_coefficients(m1, m2, s1z=0, s2z=0, phase_order=-1, spin_order=-1): qmdef2 = 1.0 M = m1 + m2 - dm = (m1-m2)/M + dm = (m1 - m2) / M m1M = m1 / M m2M = m2 / M @@ -729,18 +813,18 @@ def energy_coefficients(m1, m2, s1z=0, s2z=0, phase_order=-1, spin_order=-1): _, eta = mass1_mass2_to_mchirp_eta(m1, m2) - ecof = numpy.zeros(phase_order+1) + ecof = numpy.zeros(phase_order + 1) # Orbital terms if phase_order >= 0: ecof[0] = 1.0 if phase_order >= 1: ecof[1] = 0 if phase_order >= 2: - ecof[2] = -(1.0/12.0) * (9.0 + eta) + ecof[2] = -(1.0 / 12.0) * (9.0 + eta) if phase_order >= 3: ecof[3] = 0 if phase_order >= 4: - ecof[4] = (-81.0 + 57.0*eta - eta*eta) / 24.0 + ecof[4] = (-81.0 + 57.0 * eta - eta * eta) / 24.0 if phase_order >= 5: ecof[5] = 0 if phase_order >= 6: @@ -749,50 +833,55 @@ def energy_coefficients(m1, m2, s1z=0, s2z=0, phase_order=-1, spin_order=-1): - (155.0/96.0) *eta * eta - 35.0/5184.0 * eta * eta # Spin terms - ESO15s1 = 8.0/3.0 + 2.0*m2/m1 - ESO15s2 = 8.0/3.0 + 2.0*m1/m2 + ESO15s1 = 8.0 / 3.0 + 2.0 * m2 / m1 + ESO15s2 = 8.0 / 3.0 + 2.0 * m1 / m2 ESS2 = 1.0 / eta - EQM2s1 = qmdef1/2.0/m1M/m1M - EQM2s1L = -qmdef1*3.0/2.0/m1M/m1M + EQM2s1 = qmdef1 / 2.0 / m1M / m1M + EQM2s1L = -qmdef1 * 3.0 / 2.0 / m1M / m1M #EQM2s2 = qmdef2/2.0/m2M/m2M - EQM2s2L = -qmdef2*3.0/2.0/m2M/m2M + EQM2s2L = -qmdef2 * 3.0 / 2.0 / m2M / m2M - ESO25s1 = 11.0 - 61.0*eta/9.0 + (dm/m1M) * (-3.0 + 10.*eta/3.0) - ESO25s2 = 11.0 - 61.0*eta/9.0 + (dm/m2M) * (3.0 - 10.*eta/3.0) + ESO25s1 = 11.0 - 61.0 * eta / 9.0 + (dm / m1M) * (-3.0 + 10. * eta / 3.0) + ESO25s2 = 11.0 - 61.0 * eta / 9.0 + (dm / m2M) * (3.0 - 10. * eta / 3.0) - ESO35s1 = 135.0/4.0 - 367.0*eta/4.0 + 29.0*eta*eta/12.0 + (dm/m1M) * (-27.0/4.0 + 39.0*eta - 5.0*eta*eta/4.0) - ESO35s2 = 135.0/4.0 - 367.0*eta/4.0 + 29.0*eta*eta/12.0 - (dm/m2M) * (-27.0/4.0 + 39.0*eta - 5.0*eta*eta/4.0) + ESO35s1 = 135.0 / 4.0 - 367.0 * eta / 4.0 + 29.0 * eta * eta / 12.0 + ( + dm / m1M) * (-27.0 / 4.0 + 39.0 * eta - 5.0 * eta * eta / 4.0) + ESO35s2 = 135.0 / 4.0 - 367.0 * eta / 4.0 + 29.0 * eta * eta / 12.0 - ( + dm / m2M) * (-27.0 / 4.0 + 39.0 * eta - 5.0 * eta * eta / 4.0) - if spin_order >=3: + if spin_order >= 3: ecof[3] += ESO15s1 * s1z + ESO15s2 * s2z - if spin_order >=4: - ecof[4] += ESS2 * (s1z*s2z - 3.0*s1z*s2z) - ecof[4] += EQM2s1*s1z*s1z + EQM2s1*s2z*s2z + EQM2s1L*s1z*s1z + EQM2s2L*s2z*s2z - if spin_order >=5: - ecof[5] = ESO25s1*s1z + ESO25s2*s2z - if spin_order >=7: - ecof[7] += ESO35s1*s1z + ESO35s2*s2z + if spin_order >= 4: + ecof[4] += ESS2 * (s1z * s2z - 3.0 * s1z * s2z) + ecof[ + 4] += EQM2s1 * s1z * s1z + EQM2s1 * s2z * s2z + EQM2s1L * s1z * s1z + EQM2s2L * s2z * s2z + if spin_order >= 5: + ecof[5] = ESO25s1 * s1z + ESO25s2 * s2z + if spin_order >= 7: + ecof[7] += ESO35s1 * s1z + ESO35s2 * s2z return ecof + def energy(v, mass1, mass2, s1z=0, s2z=0, phase_order=-1, spin_order=-1): ecof = energy_coefficients(mass1, mass2, s1z, s2z, phase_order, spin_order) _, eta = mass1_mass2_to_mchirp_eta(mass1, mass2) - amp = - (1.0/2.0) * eta + amp = -(1.0 / 2.0) * eta e = 0.0 for i in numpy.arange(0, len(ecof), 1): - e += v**(i+2.0) * ecof[i] + e += v**(i + 2.0) * ecof[i] return e * amp + def meco2(m1, m2, s1z=0, s2z=0, phase_order=-1, spin_order=-1): ecof = energy_coefficients(m1, m2, s1z, s2z, phase_order, spin_order) def test(v): de = 0 for i in numpy.arange(0, len(ecof), 1): - de += v**(i+1.0)* ecof[i] * (i + 2) + de += v**(i + 1.0) * ecof[i] * (i + 2) return de @@ -800,10 +889,14 @@ def test(v): def t2_cutoff_velocity(m1, m2, chi1, chi2): - return min(meco_velocity(m1,m2,chi1,chi2), _dtdv_cutoff_velocity(m1,m2,chi1,chi2)) + return min(meco_velocity(m1, m2, chi1, chi2), + _dtdv_cutoff_velocity(m1, m2, chi1, chi2)) + def t2_cutoff_frequency(m1, m2, chi1, chi2): - return velocity_to_frequency(t2_cutoff_velocity(m1, m2, chi1, chi2), m1 + m2) + return velocity_to_frequency(t2_cutoff_velocity(m1, m2, chi1, chi2), + m1 + m2) + t4_cutoff_velocity = meco_velocity t4_cutoff_frequency = meco_frequency @@ -814,7 +907,7 @@ def t2_cutoff_frequency(m1, m2, chi1, chi2): def kerr_lightring(v, chi): """Return the function whose first root defines the Kerr light ring""" - return 1 + chi * v**3 - 3 * v**2 * (1 - chi * v**3)**(1./3) + return 1 + chi * v**3 - 3 * v**2 * (1 - chi * v**3)**(1. / 3) def kerr_lightring_velocity(chi): @@ -929,7 +1022,9 @@ def hybrid_meco_velocity(m1, m2, chi1, chi2, qm1=None, qm2=None): chi = (chi1 * m1 + chi2 * m2) / (m1 + m2) vmax = kerr_lightring_velocity(chi) - 0.01 - return minimize(hybridEnergy, 0.2, args=(m1, m2, chi1, chi2, qm1, qm2), + return minimize(hybridEnergy, + 0.2, + args=(m1, m2, chi1, chi2, qm1, qm2), bounds=[(0.1, vmax)]).x.item() @@ -963,12 +1058,20 @@ def hybrid_meco_frequency(m1, m2, chi1, chi2, qm1=None, qm2=None): if qm2 is None: qm2 = 1 - return velocity_to_frequency(hybrid_meco_velocity(m1, m2, chi1, chi2, qm1, qm2), m1 + m2) + return velocity_to_frequency( + hybrid_meco_velocity(m1, m2, chi1, chi2, qm1, qm2), m1 + m2) -def jframe_to_l0frame(mass1, mass2, f_ref, phiref=0., thetajn=0., phijl=0., - spin1_a=0., spin2_a=0., - spin1_polar=0., spin2_polar=0., +def jframe_to_l0frame(mass1, + mass2, + f_ref, + phiref=0., + thetajn=0., + phijl=0., + spin1_a=0., + spin2_a=0., + spin1_polar=0., + spin2_polar=0., spin12_deltaphi=0.): """Converts J-frame parameters into L0 frame. @@ -1031,18 +1134,29 @@ def jframe_to_l0frame(mass1, mass2, f_ref, phiref=0., thetajn=0., phijl=0., thetajn, phijl, spin1_polar, spin2_polar, spin12_deltaphi, spin1_a, spin2_a, mass1*lal.MSUN_SI, mass2*lal.MSUN_SI, f_ref, phiref) - out = {'inclination': inclination, - 'spin1x': spin1x, - 'spin1y': spin1y, - 'spin1z': spin1z, - 'spin2x': spin2x, - 'spin2y': spin2y, - 'spin2z': spin2z} + out = { + 'inclination': inclination, + 'spin1x': spin1x, + 'spin1y': spin1y, + 'spin1z': spin1z, + 'spin2x': spin2x, + 'spin2y': spin2y, + 'spin2z': spin2z + } return out -def l0frame_to_jframe(mass1, mass2, f_ref, phiref=0., inclination=0., - spin1x=0., spin1y=0., spin1z=0., - spin2x=0., spin2y=0., spin2z=0.): + +def l0frame_to_jframe(mass1, + mass2, + f_ref, + phiref=0., + inclination=0., + spin1x=0., + spin1y=0., + spin1z=0., + spin2x=0., + spin2y=0., + spin2z=0.): """Converts L0-frame parameters to J-frame. Parameters @@ -1107,11 +1221,13 @@ def l0frame_to_jframe(mass1, mass2, f_ref, phiref=0., inclination=0., lalsimulation.SimInspiralTransformPrecessingWvf2PE( inclination, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z, mass1, mass2, f_ref, phiref) - out = {'thetajn': thetajn, - 'phijl': phijl, - 'spin1_polar': s1pol, - 'spin2_polar': s2pol, - 'spin12_deltaphi': s12_deltaphi, - 'spin1_a': spin1_a, - 'spin2_a': spin2_a} + out = { + 'thetajn': thetajn, + 'phijl': phijl, + 'spin1_polar': s1pol, + 'spin2_polar': s2pol, + 'spin12_deltaphi': s12_deltaphi, + 'spin1_a': spin1_a, + 'spin2_a': spin2_a + } return out