pneqns.py

Go to the documentation of this file.

"""
Code to use post-Newtonian equations [from Ajith, PRD 84, 084037 (2011), arXiv:1107.1267] to evolve spins in a binary black hole
 
P. Ajith, A. Gupta, and N. K. Johnson-McDaniel, 04.2016, based on earlier code
"""
 
from scipy.integrate import ode
import numpy as np
import lalsimulation as lalsim
from lal import PI, MTSUN_SI, MSUN_SI, GAMMA
from numpy import log
from numpy.linalg import norm
 
 
"""Preccesion frequency spins """
def Omega(v, m1, m2, S1, S2, Ln): #S1, S2 are full spin vector while Ln is a unit vector
 
    ''' Eqs.(3.8) of Ajith (2011) (http://arxiv.org/pdf/1107.1267v2.pdf)'''
 
    m = m1+m2
    deltam = m1-m2
    eta = (m1*m2)/(m1+m2)**2.
    return (v**5./m)*((3./4+eta/2.-3.*deltam/m/4.)*Ln+0.5*v/m**2.*(-3.0*np.dot((S2 +(m2/m1)*S1), Ln)*Ln+S2)+v**2.*(9./16.+5.*eta/4.-eta**2./24.-9.*deltam/m/16.+5.*deltam*eta/8./m)*Ln)
 
 
""" compute the re-expanded dEnergy/flux """
def denergy_by_flux(v, eta, delta, chiadL, chisdL, chiasqr, chissqr, chisdchia, order):
 
    ''' Eqs.(3.2) of Ajith (2011) http://arxiv.org/pdf/1107.1267v2.pdf '''
 
    # different powers of v
    v2 = v*v; v3 = v2*v; v4 = v3*v; v5 = v4*v; v6 = v5*v; v7 = v6*v; v9 = v7*v2
 
    # initialize the cofficients
    dEbF0 = dEbF2 = dEbF3 = dEbF4 = dEbF5 = dEbF6 = dEbF6L = dEbF7 = 0.
 
    if order >= 0:
            dEbF0 = -5./(32.*eta)
    if order >= 2:
            dEbF2 = 2.2113095238095237 + (11*eta)/4.
    if order >= 3:
            dEbF3 = (113*chiadL*delta)/12. + chisdL*(9.416666666666666 - (19.*eta)/3.) - 4.*PI
    if order >= 4:
            dEbF4 = 3.010315295099521 + (233*chisdchia*delta)/48. - (719.*chiadL*chisdL*delta)/48. + \
                    chiasqr*(2.4270833333333335 - 10.*eta) + chisdL*chisdL*(-7.489583333333333 - eta/24.) + \
                    chissqr*(2.4270833333333335 + (7.*eta)/24.) + (5429.*eta)/1008. + (617*eta*eta)/144. + \
                    chiadL*chiadL*(-7.489583333333333 + 30.*eta)
    if order >= 5:
                    dEbF5 = chiadL*delta*(72.71676587301587 + (7*eta)/2.) + chisdL*(72.71676587301587 - \
                                            (1213*eta)/18. - 17*eta*eta/2.) - (7729*PI)/672. + (13*eta*PI)/8.
    if order >= 6:
            dEbF6 = -115.2253249962622 - 15211*eta*eta/6912. + 25565*eta*eta*eta/5184. + \
                                    32*PI*PI/3. + eta*(258.1491854023631 - 451*PI*PI/48.) + (1712*GAMMA)/105.
            dEbF6L = 1712./105.
    if order >= 7:
            dEbF7 = (-15419335.*PI)/1.016064e6 - (75703.*eta*PI)/6048. + (14809.*eta*eta*PI)/3024.
 
    return (dEbF0/v9)*(1. + dEbF2*v2 + dEbF3*v3 + dEbF4*v4 + dEbF5*v5 + (dEbF6+dEbF6L*log(4.*v))*v6 + dEbF7*v7)
 
 
def precession_eqns(t, y_vec, m1, m2):
        """ The coupled set of ODEs containing the PN precession eqns as well as the evolution of dv/dt
            All these equations are listed in Ajith (2011) (http://arxiv.org/pdf/1107.1267v2.pdf).
        """
 
        # unpack the vector of variables
        v, S1x, S1y, S1z, S2x, S2y, S2z, Lx, Ly, Lz = y_vec
 
        m1_sqr = m1*m1
        m2_sqr = m2*m2
        m = m1+m2
        eta = m1*m2/(m1+m2)**2
        delta = (m1-m2)/(m1+m2)
 
        # spin and angular momentum vectors
        Ln = np.array([Lx, Ly, Lz])
        S1 = np.array([S1x, S1y, S1z])
        S2 = np.array([S2x, S2y, S2z])
        chi1 = S1/m1_sqr
        chi2 = S2/m2_sqr
        chi_s = (chi1+chi2)/2.
        chi_a = (chi1-chi2)/2.
 
        #print 'v = ', v , ' chi1 = ', chi1, ' chi2 = ', chi2
 
        # dot products
        chiadL = np.dot(chi_a, Ln)
        chisdL = np.dot(chi_s, Ln)
        chissqr = np.dot(chi_s, chi_s)
        chiasqr = np.dot(chi_a, chi_a)
        chisdchia = np.dot(chi_s, chi_a)
 
        # magnitude of the orb ang momentum. The prefactor of dLN/dt in Eq.(3.6) of Ajith (2011)
        Lmag = (eta*m**2./v)*(1.+(1.5+eta/6.)*v**2.)
 
        # precession freqs. Eqs.(3.8) of Ajith (2011)
        Omega1= Omega(v, m1, m2, S1, S2, Ln)
        Omega2= Omega(v, m2, m1, S2, S1, Ln)
 
        # spin precession eqns. Eqs.(3.7) of Ajith (2011)
        dS1_dt = np.cross(Omega1, S1)
        dS2_dt = np.cross(Omega2, S2)
 
        # ang momentum precession eqn. Eqs.(3.6) of Ajith (2011)
        dL_dt =-(dS1_dt+dS2_dt)/Lmag
 
        # orb frequency evolution. Eqs.(3.5) of Ajith (2011)
        dv_dt = -1./(m*denergy_by_flux(v, eta, delta, chiadL, chisdL, chiasqr, chissqr, chisdchia, 7))
 
        return  [dv_dt, dS1_dt[0], dS1_dt[1], dS1_dt[2], dS2_dt[0], dS2_dt[1], dS2_dt[2], dL_dt[0], dL_dt[1], dL_dt[2]]
 
def evolve_spins_dt(v0, m1, m2, chi1x, chi1y, chi1z, chi2x, chi2y, chi2z, Lnx, Lny, Lnz, v_final, dt):
        """ evolve the spins and orb ang momentum according to the PN precession eqns."""
 
        # Set maximum dv/dt for stopping condition
 
        dvdt_max = 1.e-4
 
        # Set maximum deviation from v_final allowed (unless dv/dt is greater than dvdt_max or negative, in which case just print a warning)
 
        delta_vf_max = 1.e-2
 
        # Set largest deviation from normalization of orbital angular momentum allowed
 
        delta_Lnorm_max = 1.e-8
 
        # Give indices for various quantities in the solution vector
 
        V_POS = 0 # Location of v (and thus dv/dt)
 
        LNX_POS = 7 # Location of the x component of Ln
        LNY_POS = 8 # Location of the y component of Ln
        LNZ_POS = 9 # Location of the z component of Ln
 
        # some sanity checks
        if m1<0:
            raise ValueError("ERROR: mass1 is negative")
        if m2<0:
            raise ValueError("ERROR: mass2 is negative")
        if (chi1x**2+chi1y**2+chi1z**2)>1.:
            raise ValueError("ERROR: magnitude of spin1 is greater than 1")
        if (chi2x**2+chi2y**2+chi2z**2)>1.:
            raise ValueError("ERROR: magnitude of spin2 is greater than 1")
        if dt<0:
            raise ValueError("ERROR: time step is negative")
 
        m1_sqr = m1*m1
        m2_sqr = m2*m2
 
        # pack the inputs
        y_vec0 = [v0, chi1x*m1_sqr, chi1y*m1_sqr, chi1z*m1_sqr, chi2x*m2_sqr, chi2y*m2_sqr, chi2z*m2_sqr, Lnx, Lny, Lnz]
        t0 = 0.
 
        # Defining chi_eff and Newtonian chirp time. T_MAX is set to twice of chirp time.
        chi_eff = lalsim.SimInspiralTaylorF2ReducedSpinComputeChi(m1, m2, chi1z, chi2z)
        T_MAX = 2.* lalsim.SimInspiralTaylorF2ReducedSpinChirpTime(v0**3./(PI*(m1+m2)*MTSUN_SI), m1*MSUN_SI, m2*MSUN_SI, chi_eff, 0)/MTSUN_SI
        R_TOL = 1e-6
        A_TOL = 1e-6
 
        #print("T_MAX = %f"%(T_MAX))
 
        # initialize the ODE solver
        backend = "dopri5"
        solver = ode(precession_eqns)
        solver.set_integrator(backend, atol=R_TOL, rtol=R_TOL)  # nsteps=1
        solver.set_initial_value(y_vec0, t0)
        solver.set_f_params(m1, m2)
 
        y_result = []
        t_output = []
        y_result.append(y_vec0)
        t_output.append(t0)
 
        # Set diagnostic quantities to some initial values that won't trigger any stopping conditions
 
        v_current = 0.
        dvdt_current = 1.e-5
        Lnorm_current = 1.
 
        # evolve the eqns
        while solver.successful() and solver.t < 2.*T_MAX and v_current <= v_final and dvdt_current > 0. and dvdt_current < dvdt_max and abs(Lnorm_current - 1.) < delta_Lnorm_max:
 
                solver.integrate(solver.t + dt, step=1)
                y_result.append(solver.y)
                t_output.append(solver.t)
 
                #print '... t = ', solver.t, ' y = ', solver.y
 
                v_current = solver.y[V_POS]
                dvdt_current = precession_eqns(solver.t,solver.y,m1,m2)[V_POS]
                Lnorm_current = np.sqrt((solver.y[LNX_POS])**2+(solver.y[LNY_POS])**2+(solver.y[LNZ_POS])**2)
 
        Y = np.array(y_result)
        t = np.array(t_output)
 
        # Check if the integration stopped more than delta_vf_max away from v_final. If so, check if this occurred because dv/dt became larger than dvdt_max or negative. In this case, print a warning and remove the offending data, else raise an error
 
        v_last = Y.T[V_POS][-1]
 
        if abs(v_last-v_final) > delta_vf_max:
            if dvdt_current <= 0 or dvdt_current > dvdt_max:
                Y = np.delete(Y, -1, axis = 0)
                print("Warning: Integration stopped at v_max = {0}, more than {1} different from v_final = {2}, because dv/dt became negative or exceeded {3}, with a value of {4}. The final value used is {5}.".format(v_last, delta_vf_max, v_final, dvdt_max, dvdt_current, precession_eqns(t,Y[-1],m1,m2)[V_POS]))
            else:
                raise ValueError("v_max = {0} is more than {1} different from v_final = {2} and dv/dt is positive and less than {3}".format(v_v[-1], delta_vf, v_final, dvdt_max))
 
        if abs(Lnorm_current - 1.) >= delta_Lnorm_max:
            raise ValueError("norm of Ln is more than {0:e} different from 1, with distance {1:e}".format(delta_Lnorm_max, abs(Lnorm_current - 1.)))
 
        v_v, S1x_v, S1y_v, S1z_v, S2x_v, S2y_v, S2z_v, Lx_v, Ly_v, Lz_v  = Y.T
 
        return v_v, S1x_v/m1_sqr, S1y_v/m1_sqr, S1z_v/m1_sqr, S2x_v/m2_sqr, S2y_v/m2_sqr, S2z_v/m2_sqr, Lx_v, Ly_v, Lz_v
 
 
def find_tilts_and_phi12_at_freq(v0, m1, m2, chi1x, chi1y, chi1z, chi2x, chi2y, chi2z, Lnx, Lny, Lnz, v_final, dt):
        """ given the spins and ang momentum at a given frequency, find the tilt and in-plane spin angles at a later frequency """
 
        print("v0 = %f, m1 = %f, m2 = %f, chi1x = %f, chi1y = %f, chi1z = %f, chi2x = %f, chi2y = %f, chi2z = %f"%(v0, m1, m2, chi1x, chi1y, chi1z, chi2x, chi2y, chi2z))
 
        # evolve the spins
        v_v, chi1x_v, chi1y_v, chi1z_v, chi2x_v, chi2y_v, chi2z_v, Lnx_v, Lny_v, Lnz_v = evolve_spins_dt(v0, m1, m2, chi1x, chi1y, chi1z, chi2x, chi2y, chi2z, Lnx, Lny, Lnz, v_final, dt)
 
        chi1_v = np.array([chi1x_v[-1], chi1y_v[-1], chi1z_v[-1]])
        chi2_v = np.array([chi2x_v[-1], chi2y_v[-1], chi2z_v[-1]])
 
        Ln_v = np.array([Lnx_v[-1], Lny_v[-1], Lnz_v[-1]])
 
        # norms and normalizing
        chi1_norm = norm(chi1_v)
        chi2_norm = norm(chi2_v)
 
        Ln_v /= norm(Ln_v)
 
        # dot products
        chi1dL_v = np.dot(chi1_v, Ln_v)
        chi2dL_v = np.dot(chi2_v, Ln_v)
 
        # in-plane spins
 
        chi1inplane = chi1_v - chi1dL_v*Ln_v
        chi2inplane = chi2_v - chi2dL_v*Ln_v
 
        # Defining cosine of tilts and phi12
        cos_tilt1 = chi1dL_v/chi1_norm
        cos_tilt2 = chi2dL_v/chi2_norm
        cos_phi12 = np.dot(chi1inplane,chi2inplane)/(norm(chi1inplane)*norm(chi2inplane))
 
        print("cos tilt1 = %f, cos tilt2 = %f, cos phi12 = %f"%(cos_tilt1, cos_tilt2, cos_phi12))
 
        return np.arccos(cos_tilt1), np.arccos(cos_tilt2), np.arccos(cos_phi12)