https://github.com/gwastro/pycbc
Tip revision: 678515e6ed9ced69a21d3ad49806d6029dfd595e authored by Alexander Harvey Nitz on 13 June 2015, 23:49:46 UTC
update version to 1.0.0
update version to 1.0.0
Tip revision: 678515e
pnutils.py
# Copyright (C) 2012 Alex Nitz
#
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# =============================================================================
#
# Preamble
#
# =============================================================================
#
"""This module contains convenience pN functions. This includes calculating conversions
between quantities.
"""
from __future__ import division
import lal, lalsimulation
from numpy import log
import numpy
from scipy.optimize import bisect
def chirp_distance(dist, mchirp, ref_mass=1.4):
return dist * (2.**(-1./5) * ref_mass / mchirp)**(5./6)
def mass1_mass2_to_mtotal_eta(mass1, mass2):
m_total = mass1 + mass2
eta = (mass1 * mass2) / (m_total * m_total)
return m_total,eta
def mtotal_eta_to_mass1_mass2(m_total, eta):
mass1 = 0.5 * m_total * (1.0 + (1.0 - 4.0 * eta)**0.5)
mass2 = 0.5 * m_total * (1.0 - (1.0 - 4.0 * eta)**0.5)
return mass1,mass2
def mass1_mass2_to_mchirp_eta(mass1, mass2):
m_total, eta = mass1_mass2_to_mtotal_eta(mass1, mass2)
m_chirp = m_total * eta**(3./5.)
return m_chirp,eta
def mchirp_eta_to_mass1_mass2(m_chirp, eta):
M = m_chirp / (eta**(3./5.))
return mtotal_eta_to_mass1_mass2(M, eta)
def mchirp_mass1_to_mass2(mchirp, mass1):
"""
This function takes a value of mchirp and one component mass and returns
the second component mass. As this is a cubic equation this requires
finding the roots and returning the one that is real.
Basically it can be shown that:
m2^3 - a(m2 + m1) = 0
where
a = Mc^5 / m1^3
this has 3 solutions but only one will be real.
"""
a = mchirp**5 / mass1**3
roots = numpy.roots([1,0,-a,-a*mass1])
# Find the real one
real_root = roots[(abs(roots - roots.real)).argmin()]
return real_root.real
def eta_mass1_to_mass2(eta, mass1, return_mass_heavier=False):
"""
This function takes values for eta and one component mass and returns the
second component mass. Similar to mchirp_mass1_to_mass2 this requires
finding the roots of a quadratic equation. Basically:
eta m2^2 + (2 eta - 1)m1 m2 + \eta m1^2 = 0
This has two solutions which correspond to mass1 being the heavier mass
or it being the lighter mass. By default the value corresponding to
mass1 > mass2 is returned. Use the return_mass_heavier kwarg to invert this
behaviour.
"""
roots = numpy.roots([eta, (2*eta - 1)*mass1, mass1*mass1*eta])
if return_mass_heavier==False:
return roots[roots.argmin()]
else:
return roots[roots.argmax()]
def A0(f_lower):
"""used in calculating chirp times: see Cokelaer, arxiv.org:0706.4437
appendix 1, also lalinspiral/python/sbank/tau0tau3.py
"""
return 5. / (256. * (lal.PI * f_lower)**(8./3.))
def A3(f_lower):
"""another parameter used for chirp times"""
return lal.PI / (8. * (lal.PI * f_lower)**(5./3.))
def mass1_mass2_to_tau0_tau3(mass1, mass2, f_lower):
m_total,eta = mass1_mass2_to_mtotal_eta(mass1, mass2)
# convert to seconds
m_total = m_total * lal.MTSUN_SI
# formulae from arxiv.org:0706.4437
tau0 = A0(f_lower) / (m_total**(5./3.) * eta)
tau3 = A3(f_lower) / (m_total**(2./3.) * eta)
return tau0,tau3
def tau0_tau3_to_mtotal_eta(tau0, tau3, f_lower):
m_total = (tau3 / A3(f_lower)) / (tau0 / A0(f_lower))
eta = m_total**(-2./3.) * (A3(f_lower) / tau3)
# convert back to solar mass units
return (m_total/lal.MTSUN_SI),eta
def tau0_tau3_to_mass1_mass2(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):
M, eta = mass1_mass2_to_mtotal_eta(mass1, mass2)
# get_beta_sigma_from_aligned_spins() takes
# the spin of the heaviest body first
heavy_spin = numpy.where(mass2 <= mass1, spin1z, spin2z)
light_spin = numpy.where(mass2 > mass1, spin1z, spin2z)
beta, sigma, gamma, xs = get_beta_sigma_from_aligned_spins(
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.
See <http://arxiv.org/pdf/0810.5336v3.pdf>.
Parameters
-----------
eta : float or numpy.array
Symmetric mass ratio of the input system(s)
spin1z : float or numpy.array
Spin(s) parallel to the orbit of the heaviest body(ies)
spin2z : float or numpy.array
Spin(s) parallel to the orbit of the smallest body(ies)
Returns
--------
beta : float or numpy.array
The 1.5PN spin combination
sigma : float or numpy.array
The 2PN spin combination
gamma : float or numpy.array
The 2.5PN spin combination
chis : float or numpy.array
(spin1z + spin2z) / 2.
"""
chiS = 0.5 * (spin1z + spin2z)
chiA = 0.5 * (spin1z - spin2z)
delta = (1 - 4 * eta) ** 0.5
spinspin = spin1z * spin2z
beta = (113. / 12. - 19. / 3. * eta) * chiS
beta += 113. / 12. * delta * chiA
sigma = eta / 48. * (474 * spinspin)
sigma += (1 - 2 * eta) * (81. / 16. * (chiS * chiS + chiA * chiA))
sigma += delta * (81. / 8. * (chiS * chiA))
gamma = (732985. / 2268. - 24260. / 81. * eta - \
340. / 9. * eta * eta) * chiS
gamma += (732985. / 2268. + 140. / 9. * eta) * delta * chiA
return beta, sigma, gamma, chiS
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 v**(3.0) / (M * lal.MTSUN_SI * lal.PI)
def frequency_to_velocity(f, M):
return (lal.PI * M * lal.MTSUN_SI * f)**(1.0/3.0)
def f_SchwarzISCO(M):
"""
Innermost stable circular orbit (ISCO) for a test particle
orbiting a Schwarzschild black hole
Parameters
----------
M : float or numpy.array
Total mass in solar mass units
Returns
-------
f : float or numpy.array
Frequency in Hz
"""
return velocity_to_frequency((1.0/6.0)**(0.5), M)
def f_BKLISCO(m1, m2):
"""
Mass ratio dependent ISCO derived from estimates of the final spin
of a merged black hole in a paper by Buonanno, Kidder, Lehner
(arXiv:0709.3839). See also arxiv:0801.4297v2 eq.(5)
Parameters
----------
m1 : float or numpy.array
First component mass in solar mass units
m2 : float or numpy.array
Second component mass in solar mass units
Returns
-------
f : float or numpy.array
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 )
def f_LightRing(M):
"""
Gravitational wave frequency corresponding to the light-ring orbit,
equal to 1/(3**(3/2) pi M) : see InspiralBankGeneration.c
Parameters
----------
M : float or numpy.array
Total mass in solar mass units
Returns
-------
f : float or numpy.array
Frequency in Hz
"""
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)
found to give good fit between stationary-phase templates and
numerical relativity waveforms [NB equal-mass & nonspinning!]
Equal to 1.07*omega_220/2*pi
Parameters
----------
M : float or numpy.array
Total mass in solar mass units
Returns
-------
f : float or numpy.array
Frequency in Hz
"""
return 1.07 * 0.5326 / (2*lal.PI * 0.955 * M * lal.MTSUN_SI)
def f_FRD(m1, m2):
"""
Fundamental RingDown frequency calculated from the Berti, Cardoso and
Will (gr-qc/0512160) value for the omega_220 QNM frequency using
mass-ratio dependent fits to the final BH mass and spin from Buonanno
et al. (arXiv:0706.3732) : see also InspiralBankGeneration.c
Parameters
----------
m1 : float or numpy.array
First component mass in solar mass units
m2 : float or numpy.array
Second component mass in solar mass units
Returns
-------
f : float or numpy.array
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)
def f_LRD(m1, m2):
"""
Lorentzian RingDown frequency = 1.2*FRD which captures part of
the Lorentzian tail from the decay of the QNMs
Parameters
----------
m1 : float or numpy.array
First component mass in solar mass units
m2 : float or numpy.array
Second component mass in solar mass units
Returns
-------
f : float or numpy.array
Frequency in Hz
"""
return 1.2 * f_FRD(m1, m2)
def _get_freq(freqfunc, m1, m2, s1z, s2z):
"""
Wrapper of the LALSimulation function returning the frequency
for a given frequency function and template parameters.
Parameters
----------
freqfunc : lalsimulation FrequencyFunction wrapped object e.g.
lalsimulation.fEOBNRv2RD
m1 : float-ish, i.e. castable to float
First component mass in solar masses
m2 : float-ish
Second component mass in solar masses
s1z : float-ish
First component dimensionless spin S_1/m_1^2 projected onto L
s2z : float-ish
Second component dimensionless spin S_2/m_2^2 projected onto L
Returns
-------
f : float
Frequency in Hz
"""
# 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))
# 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
for the given frequency function and template parameters.
Parameters
----------
freqfunc : string
Name of the frequency function to use, e.g., 'fEOBNRv2RD'
m1 : float or numpy.array
First component mass in solar masses
m2 : float or numpy.array
Second component mass in solar masses
s1z : float or numpy.array
First component dimensionless spin S_1/m_1^2 projected onto L
s2z : float or numpy.array
Second component dimensionless spin S_2/m_2^2 projected onto L
Returns
-------
f : float or numpy.array
Frequency in Hz
"""
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)
frequency for a given approximant an template parameters
Parameters
----------
approx : lalsimulation approximant wrapped object e.g.
lalsimulation.EOBNRv2
m1 : float-ish, i.e. castable to float
First component mass in solar masses
m2 : float-ish
Second component mass in solar masses
s1z : float-ish
First component dimensionless spin S_1/m_1^2 projected onto L
s2z : float-ish
Second component dimensionless spin S_2/m_2^2 projected onto L
Returns
-------
f : float
Frequency in Hz
"""
# 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))
# 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
(highest) frequency for a given approximant using given template
parameters.
NOTE: TaylorTx and TaylorFx are currently all given an ISCO cutoff !!
Parameters
----------
approx : string
Name of the approximant e.g. 'EOBNRv2'
m1 : float or numpy.array
First component mass in solar masses
m2 : float or numpy.array
Second component mass in solar masses
s1z : float or numpy.array
First component dimensionless spin S_1/m_1^2 projected onto L
s2z : float or numpy.array
Second component dimensionless spin S_2/m_2^2 projected onto L
Returns
-------
f : float or numpy.array
Frequency in Hz
"""
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 m1, m2, s1z, s2z
named_frequency_cutoffs = {
# functions depending on the total mass alone
"SchwarzISCO": lambda p: f_SchwarzISCO(p["m1"]+p["m2"]),
"LightRing" : lambda p: f_LightRing(p["m1"]+p["m2"]),
"ERD" : lambda p: f_ERD(p["m1"]+p["m2"]),
# functions depending on the 2 component masses
"BKLISCO" : lambda p: f_BKLISCO(p["m1"], p["m2"]),
"FRD" : lambda p: f_FRD(p["m1"], p["m2"]),
"LRD" : lambda p: f_LRD(p["m1"], p["m2"]),
# functions depending on 2 component masses and aligned spins
"MECO" : lambda p: meco_frequency(p["m1"], p["m2"],
p["s1z"], p["s2z"]),
"IMRPhenomBFinal": lambda p: get_freq("fIMRPhenomBFinal",
p["m1"], p["m2"],
p["s1z"], p["s2z"]),
"IMRPhenomCFinal": lambda p: get_freq("fIMRPhenomCFinal",
p["m1"], p["m2"],
p["s1z"], p["s2z"]),
"EOBNRv2RD" : lambda p: get_freq("fEOBNRv2RD", p["m1"], p["m2"],
p["s1z"], p["s2z"]),
"EOBNRv2HMRD" : lambda p: get_freq("fEOBNRv2HMRD", p["m1"], p["m2"],
p["s1z"], p["s2z"]),
"SEOBNRv1RD" : lambda p: get_freq("fSEOBNRv1RD", p["m1"], p["m2"],
p["s1z"], p["s2z"]),
"SEOBNRv1Peak": lambda p: get_freq("fSEOBNRv1Peak", p["m1"], p["m2"],
p["s1z"], p["s2z"]),
"SEOBNRv2RD" : lambda p: get_freq("fSEOBNRv2RD", p["m1"], p["m2"],
p["s1z"], p["s2z"]),
"SEOBNRv2Peak": lambda p: get_freq("fSEOBNRv2Peak", p["m1"], p["m2"],
p["s1z"], p["s2z"])
}
def frequency_cutoff_from_name(name, m1, m2, s1z, s2z):
"""
Returns the result of evaluating the frequency cutoff function
specified by 'name' on a template with given parameters.
Parameters
----------
name : string
Name of the cutoff function
m1 : float or numpy.array
First component mass in solar masses
m2 : float or numpy.array
Second component mass in solar masses
s1z : float or numpy.array
First component dimensionless spin S_1/m_1^2 projected onto L
s2z : float or numpy.array
Second component dimensionless spin S_2/m_2^2 projected onto L
Returns
-------
f : float or numpy.array
Frequency in Hz
"""
params = {"m1":m1, "m2":m2, "s1z":s1z, "s2z":s2z}
return named_frequency_cutoffs[name](params)
##############################This code was taken from Andy ###########
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
chisym = (chi1 + chi2) / 2.
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.
energy3 = 0.
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.
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)
Parameters
----------
m1 : float
First component mass in solar masses
m2 : float
Second component mass in solar masses
chi1 : float
First component dimensionless spin S_1/m_1^2 projected onto L
chi2 : float
Second component dimensionless spin S_2/m_2^2 projected onto L
Returns
-------
v : float
Velocity (dimensionless)
"""
energy0, 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)
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
chisym = (chi1 + chi2) / 2.
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
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)
return (dtdv0, dtdv2, dtdv3, dtdv4, dtdv5, dtdv6, dtdv6log, dtdv7)
def _dtdv_cutoff_velocity(m1, m2, chi1, chi2):
dtdv0, dtdv2, dtdv3, dtdv4, dtdv5, dtdv6, dtdv6log, dtdv7 = _dtdv_coeffs(m1, m2, chi1, chi2)
def dtdv_func(v):
return 1. + v * v * (dtdv2 + v * (dtdv3 \
+ v * (dtdv4
+ v * (dtdv5
+ v * ((dtdv6 + dtdv6log*3.*log(v))
+ v * dtdv7)))))
if dtdv_func(1.0) < 0.:
return bisect(dtdv_func, 0.05, 1.0)
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")
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")
elif spin_order == -1:
spin_order = implemented_spin_order
qmdef1 = 1.0
qmdef2 = 1.0
M = m1 + m2
dm = (m1-m2)/M
m1M = m1 / M
m2M = m2 / M
s1z = s1z * m1M * m1M
s2z = s2z * m2M * m2M
mchirp, eta = mass1_mass2_to_mchirp_eta(m1, m2)
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)
if phase_order >= 3:
ecof[3] = 0
if phase_order >= 4:
ecof[4] = (-81.0 + 57.0*eta - eta*eta) / 24.0
if phase_order >= 5:
ecof[5] = 0
if phase_order >= 6:
ecof[6] = - 675.0/64.0 + ( 34445.0/576.0 \
- 205.0/96.0 * lal.PI * lal.PI ) * eta \
- (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
ESS2 = 1.0 / eta
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
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)
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
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)
mchirp, eta = mass1_mass2_to_mchirp_eta(mass1, mass2)
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]
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)
return de
return bisect(test, 0.001, 1.0)
def t2_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)
t4_cutoff_velocity = meco_velocity
t4_cutoff_frequency = meco_frequency
