Notebook

The Spinning Effective One-Body Initial Condition Solver¶

Author: Tyler Knowles¶

This module documents the reduced spinning effective one-body initial condition solver as numerically implemented in LALSuite's SEOBNRv3 gravitational waveform approximant. That is, we follow Section IV A of Buonanno, Chen, and Damour (2006).¶

Notebook Status: In progress

Validation Notes: This module is under active development -- do *not* use the resulting code for scientific applications. In the future, this module will be validated against the LALSuite SEOBNRv3/SEOBNRv3_opt code that was reviewed and approved for LIGO parameter estimation by the LIGO Scientific Collaboration.

Introduction¶

The Physical System of Interest¶

Consider two compact objects (e.g. black holes or neutron stars) with masses $m_{1}$, $m_{2}$ (in solar masses) and spin angular momenta ${\bf S}_{1}$, ${\bf S}_{2}$ in a binary system. The spinning effective one-body ("SEOB") Hamiltonian $H_{\rm real}$ (see BB2010 Equation (5.69)) describes the dynamics of this system. We seek initial conditions for nonadiabatic evolutions of such a system, and follow BCD2006 Section IV A.

To compute the initial conditions, we begin with the following system parameters:

the mass of each compact object, denoted $m_{1}$, $m_{2}$,
the spin vector of each compact object, denoted ${\bf S}_{1}$, ${\bf S}_{2}$, and
initial orbital frequency $f$.

We choose a right-handed spatial coordinate basis $\left\{ {\bf e}_{0}, {\bf e}_{1}, {\bf e}_{2} \right\}$ so that the initial separation vector ${\bf r}$ between the compact objects lies along the ${\bf e}_{0}$-axis and the orbital plane coincides with the ${\bf e}_{0}$, ${\bf e}_{1}$-plane. Assume that ${\bf S}_{1}$, ${\bf S}_{2}$ are written in this basis. Our goal is to produce initial dynamical variables

${\bf x} = \left( x, y, z \right)$, and
${\bf p} = \left( p_{x}, p_{y}, p_{z} \right)$.

We include below the physical parameters necessary to compute the initial conditions. Besides the physical parameters, we also need the Euler–Mascheroni constant $\gamma$ and the geomtrized solar mass $\mathcal{M}_{\odot}$, both hard-coded in LALSuite with the significant digits shown below. (The following links point directly to the appropriate LALSuite documentation: $\gamma$ and $\mathcal{M}_{\odot}$.)

Please note that throughout this notebook we adpot the following conventions:

$c = G = 1$ where $c$ is the speed of light in a vacuum and $G$ is Newton's gravitational constant,
$m_{1} \ge m_{2}$,
hatted vectors (e.g. $\hat{\bf L}_{N}$) usually denote scaled or unit vectors, and
the initial inclination angle $\iota$ of the system relative to some observer is chosen to be zero.

Please note that in BCD2006 the initial conditions are solved for given an initial separation; here we use a given initial frequency instead. The difference is in our approach to solving Equation (4.8). Our approach also differs from that found in LALSuite's SEOBNRv3 code XLALSimIMRSpinEOBInitialConditionsPrec() function (file: LALSimIMRSpinEOBInitialConditionsPrec.c) because we choose our intial coordinate system so that the inclination angle $\iota$ is zero and $m_{1} \ge m_{2}$.

Citations¶

Throughout this module, we refer to

Buonanno, Chen, and Damour (2006) as BCD2006,
Barausse and Buonanno (2010) as BB2010,
Taracchini, et. al. (2012) as T2012,
Damour, et. al. (2009) as DIN2009, and
Pan, et. al. (2014) as P2014.

LALSuite line numbers are taken from Git commit bba40f2 (see LALSuite's GitLab page).

In [1]:

# Initial condition solver for the spinning effective one-body formulation
# See https://arxiv.org/abs/gr-qc/0508067 Section IV A, which we refer to as BCD2006
# Import necessary NumPy, SymPy, and SEOBNR modules
import numpy as np
import os.path
from scipy.optimize import root
from scipy.interpolate import interp1d, interp2d
from numpy.linalg import norm
import SEOBNR.NQC_corrections as nqc
import SEOBNR.nqc_interp as nqi

# For testing, remove numpy and sympy expression files
# For now, do NOT regenerate CSE expressions
import shutil, os
import sys#TylerK: Add sys to get cmdline_helper from NRPy top directory; remove this line and next when debugged
sys.path.append('../')
!rm -r SEOBNR_Playground_Pycodes
outdir = os.path.join("SEOBNR_Playground_Pycodes/")
import cmdline_helper as cmd
cmd.mkdir(outdir)
with open(outdir+"__init__.py", "w") as file:
    file.write("")

# Input variables: will eventually structure this module as a function with the following input parameters
# m1, m2 given in solar masses, f in Hz, and spin in
m1 = 23.
m2 = 10.
f = 20.
S1 = np.array([0.01, 0.02, -0.03])
S2 = np.array([0.04, -0.05, 0.06])

# Initial conditions are computed with tortoise = 0; we later convert momentum if necessary
# See LALSuite's LALSimIMRSpinEOBInitialConditionsPrec.c Line 775 and the discussion
# preceeding Equation (14) of Taracchini, et. al. 2012 (https://arxiv.org/abs/1202.0790)
tortoise = 0.0

# The values of the following constants are from LALSuite (see LALSuite documentation at
# https://lscsoft.docs.ligo.org/lalsuite/lal/group___l_a_l_constants__h.html).
# Euler–Mascheroni constant $\gamma$
EMgamma = 0.577215664901532860606512090082402431
# Geomtrized solar mass $\mathcal{M}_{\odot}$
Msol = 4.925491025543575903411922162094833998e-6

#Convert the spins to dimensionless quantities
S1 *= m1*m1
S2 *= m2*m2

Table of Contents¶

$$\label{toc}$$

This notebook is organized as follows, matching the "steps" listed in BCD2006:

Step 1: Initial Coordinate Choice

Step 1.a: Mass terms
Step 1.b: Spin terms
Step 1.c: Normalized Orbital Angular Momenutm $\hat{\bf L}_{N}$
Step 1.d: Normalized Position $\hat{\bf r}$
Step 1.e: Normalized Velocity $\hat{\bf v}$
Note

Step 2: Compute ${\bf r}$, ${\bf p}_{r}$, ${\bf p}_{\theta}$, and ${\bf p}_{\phi}$

Step 2.a: $\omega$
Step 2.b: Initial Velocity $v$
Step 2.c: ${\bf S}_{\rm Kerr}$
Step 2.d: Root finding

Step 3: Rotate $\hat{\bf L} \to {\bf e}_{z}$

Note
Step 3.a: Normalize ${\bf q}$ and ${\bf p}$
Step 3.b: $\hat{\bf L}$
Step 3.c: Rotation matrix
Step 3.d: Rotate $\hat{\bf r}$
Step 3.e: Rotate $\hat{\bf v}$
Step 3.f: Rotate $\hat{\bf L}_{N}$
Step 3.g: Rotate ${\bf S}_{1}$
Step 3.h: Rotate ${\bf S}_{2}$
Step 3.i: Rotate $\hat{\bf S}_{1}$
Step 3.j: Rotate $\hat{\bf S}_{2}$
Step 3.k: Rotate ${\bf q}$
Step 3.l: Rotate ${\bf p}$

Step 4: Compute $\dot{\bf r}$

Step 4.a: Convert from Cartesian to Spherical Coordinates
Step 4.b: Second partial derivatives of $H_{\rm real}$
Stop 4.c: $\frac{ \partial E }{ \partial r }$
Step 4.e: $\boldsymbol{\sigma}^{*}$
Step 4.f: $H_{\rm real}$

Step 5: Invert the rotation of Step 3
Output: Output this notebook to $\LaTeX$-formatted PDF file
Validation: Perform validation checks against LALSuite's SEOBNRv3 code (commit bba40f2)

Step 1: Initial Coordinate Choice [Back to top]¶

$$\label{step1}$$

Step 1.a: Mass terms [Back to top]¶

$$\label{massterms}$$

Following the notation preceeding BCD2006 Equation (2.2), we define the total mass of the system $M$ and the symmetric mass ratio $\eta$:

\begin{align*} M &= m_{1} + m_{2} \\ \eta &= \frac{ m_{1} m_{2} }{ M^{2} } \end{align*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Lines 762--763.

In [2]:

# Binary system total mass $M$
M = m1 + m2

# Inverse mass terms used repeatedly when computing initial conditions
Minv = np.divide(1,M)
Msqinv = Minv*Minv

# Symmetric mass ratio $\eta$
eta = m1*m2*Msqinv

Step 1.b: Spin terms [Back to top]¶

$$\label{spinterms}$$

Since we assumed $G = c = 1$, we normalize and make the spin angular momenta dimensionless via:

\begin{align*} \hat{\bf S}_{1} &= \frac{ 1 }{ M^{2} } {\bf S}_{1} \\ \hat{\bf S}_{2} &= \frac{ 1 }{ M^{2} } {\bf S}_{2} \end{align*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Lines 768--771.

In [3]:

# Normalized, dimensionless spin vectors
S1hat = Msqinv*S1
S2hat = Msqinv*S2

Step 1.c: Normalized Orbital Angular Momenutm $\hat{\bf L}_{N}$ [Back to top]¶

$$\label{ln}$$

Since we assume that the initial separation vector ${\bf r}$ between $m_{1}$ and $m_{2}$ lies along the ${\bf e}_{0}$-axis and the initial orbital plane coincides with the ${\bf e}_{0},{\bf e}_{1}$-plane, the normalized inital orbital angular momentum vector $\hat{\bf L}_{N}$ is given by

\begin{equation*} \hat{\bf L}_{N} = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Lines 787--789.

In [4]:

# Normalized orbital angular momentum
LNhat = np.array([0., 0., 1.])

Step 1.d: Normalized Position $\hat{\bf r}$ [Back to top]¶

$$\label{rhat}$$

We assumed that the initial separation vector ${\bf r}$ lies along the ${\bf e}_{0}$-axis, so the normalized initial separation vector $\hat{\bf r}$ is given by

\begin{equation*} \hat{\bf r} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}. \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Lines 801--803.

In [5]:

# Normalized position vector
rhat = np.array([1., 0., 0.])

Step 1.e: Normalized Velocity $\hat{\bf v}$ [Back to top]¶

$$\label{vhat}$$

Given normalized orbital angular momentum ($\hat{\bf L}_{N}$) and normalized position ($\hat{\bf r}$), the normalized velocity vector ($\hat{\bf v}$) is given by

\begin{equation*} \hat{\bf v} = \frac{ \hat{\bf L}_{N} \times \hat{\bf r} }{ \left\lvert \hat{\bf L}_{N} \times \hat{\bf r} \right\rvert }. \end{equation*}

Given $\hat{\bf L}_{N} = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}$ and $\hat{\bf r} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$ it is clear that $\hat{\bf v} = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}$.

See LALSimIMRSpinEOBInitialConditionsPrec.c Lines 807--811.

In [6]:

# Normalized velocity vector
vhat = np.array([0., 1., 0.])

Note [Back to top]¶

$$\label{step1note}$$

Since we began assuming $\iota = 0$, we do not need to rotate $\hat{\bf r}$, $\hat{\bf v}$, $\hat{\bf L}_{N}$, ${\bf S}_{1}$, ${\bf S}_{2}$, $\hat{\bf S}_{1}$, or $\hat{\bf S}_{2}$ as is done at LALSimIMRSpinEOBInitialConditionsPrec.c Lines 840-847 (Step 1 of BCD2006 Section IV A). In particular, the rotation matrix in this case is the $3\times3$ identity matrix.

Step 2: Compute ${\bf r}$ and ${\bf p}$ in spherical coordinates [Back to top]¶

$$\label{step2}$$

We seek postion vector ${\bf r}$ and ${\bf p}$ assuming a spherical orbit without radiation reaction.

Step 2.a: Initial orbital frequency $\omega$ [Back to top]¶

$$\label{omega}$$

Noting that the plane of the polarization of the gravitational wave "rotates at twice the orbital rate" (see the "Effects of passing" section of this Wikipedia article), the initial orbital frequency is

\begin{equation*} \omega = M \mathcal{M}_{\odot} \pi f. \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 893.

In [7]:

# Omega: initial orbital angular frequency
omega = M*Msol*np.pi*f

Step 2.b: Initial Velocity $v$ [Back to top]¶

$$\label{velocity}$$

Is there a paper reference for this formula? Zach suggested Kepler's Laws, but a cursory look didn't reveal a convincing link.

\begin{equation*} v = \sqrt[3]{ \omega }. \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 894.

In [8]:

# v: initial velocity and velocity squared, since we use that quantity often
v = np.cbrt(omega)
vsq = v*v

Step 2.c: ${\bf S}_{\rm Kerr}$ [Back to top]¶

$$\label{skerr}$$

This cell may be unecessary because we compute a in the derivatives (and spins depned on time so $a$ is time-dependent!).

From BB2010 Equations (5.2), (5.64), and (5.67) we have

\begin{equation*} {\bf S}_{\rm Kerr} = {\bf S}_{1} + {\bf S}_{2}. \end{equation*}

Taking the square of BB2010 Equation (4.9),

\begin{equation*} a^{2} = \frac{ {\bf S}_{\rm Kerr} \cdot {\bf S}_{\rm Kerr} }{ M^{2} } \end{equation*}

so that

\begin{equation*} a = \sqrt{ a^{2} }. \end{equation*}

In [9]:

# Compute S_Kerr, the spin of the deformed Kerr background
# See https://arxiv.org/abs/0912.3517 Equations (5.2), (5.64), and (5.67)
SKerr = np.add(S1, S2)

# Normalize S_Kerr by total mass
SKerr *= Msqinv

# Compute a, which is a parameter in metric potentials of a Kerr spacetime
# See https://arxiv.org/abs/0912.3517 Equation (4.9)
asq = np.dot(SKerr,SKerr)
a = np.sqrt(asq)

Step 2.d: Root-finding [Back to top]¶

$$\label{rootfinding}$$

We will write components of the momentum vector ${\bf p}$ in spherical coordinates with components ${\bf p}_{r}$, ${\bf p}_{\theta}$, and ${\bf p}_{\phi}$. In the special case in which we find ourselves, we have (see BCD2006 Equations (4.7) and (4.9)):

\begin{align*} {\bf r}^{\theta} &= \frac{ \pi }{ 2 } \\ {\bf r}^{\phi} &= 0 \\ {\bf p}_{r} &= 0. \end{align*}

From BCD2006 Equations (4.8)--(4.9), we seek to solve

\begin{equation*} \begin{bmatrix} \frac{ \partial H }{ \partial {\bf r}^{r} } \\ \frac{ \partial H }{ \partial {\bf p}^{\theta} } \\ \frac{ \partial H }{ \partial {\bf p}^{\phi} } - \omega \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}. \end{equation*}

As the Hamiltonian is given in Cartesian coordinates, this requires computing $\frac{ \partial H }{ \partial {\bf r}^{0} }$, $\frac{ \partial H }{ \partial {\bf p}^{1} }$, and $\frac{ \partial H }{ \partial {\bf p}^{2} }$ and then converting to spherical coordinates. That is, using the chain rule and recalling $\phi = 0$ and $\theta = \frac{ \pi }{ 2 }$, we find

\begin{align*} \frac{\partial H}{\partial {\bf r}^{r}} &= \frac{\partial H}{\partial {\bf r}^{0}} - \frac{\frac{\partial H}{\partial {\bf p}^{1}}{\bf p}^{\phi}}{\left({\bf r}^{r}\right)^{2}} + \frac{\frac{\partial H}{\partial {\bf p}^{2}}{\bf p}^{\theta}}{\left({\bf r}^{r}\right)^{2}} \\ \frac{\partial H}{\partial {\bf p}^{\theta}} &= -\frac{\frac{\partial H}{\partial {\bf p}^{2}}}{{\bf r}^{r}} \\ \frac{\partial H}{\partial {\bf p}^{\phi}} &= \frac{\frac{\partial H}{\partial {\bf p}^{1}}}{{\bf r}^{r}}. \end{align*}

Validation note: Using input parameters \begin{equation*} {\rm -M\ 23\ -m\ 10\ -f\ 20\ -X\ 0.01\ -Y\ 0.02\ -Z\ -0.03\ -x\ 0.04\ -y\ -0.05\ -z\ 0.06} \end{equation*} in LALSuite (see LALSimIMRSpinEOBInitialConditionsPrec.c lines 1075--1077), we find that the output of the root finder is \begin{align*} {\bf q}[0] &= 2.129681018601393{\rm e}+01,\\ {\bf p}[1] &= 2.335391115414913{\rm e}-01,\\ {\bf p}[2] &= 2.780558832447243{\rm e}-06. \end{align*} SEOBNRv3_pert gives \begin{align*} {\bf q}[0] &= 2.129680598646680{\rm e}+01,\\ {\bf p}[1] &= 2.335390056385869{\rm e}-01,\\ {\bf p}[2] &= 2.646198198972773{\rm e}-06. \end{align*} Note that LALSuite therefore preserves 6, 6, and 1 signficant digits, respectively. This Jupyter notebook, using those same inputs, gives roots \begin{align*} {\bf q}[0] &= 21.29680702671223,\\ {\bf p}[1] &= 0.2335390343179406,\\ {\bf p}[2] &= 2.646124517885955{\rm e}-06. \end{align*} The perturbed inputs give \begin{align*} {\bf q}[0] &= 21.2968070267122,\\ {\bf p}[1] &= 0.23353903431794046,\\ {\bf p}[2] &= 2.6461245179212286{\rm e}-06. \end{align*} That is, we perserve 15, 15, and 10 significant digits, respectively. We therefore conclude that the root-finders agree to as many significant digits as possible.

In [10]:

# First run the notebook Tutorial-SEOBNR_Documentation, if that hasn't already been done
if not os.path.isfile("SEOBNR/Hamiltonian-Hreal_one_line_expressions.txt"):
    %run Tutorial-SEOBNR_Documentation.ipynb

# We need to reverse the expressions from the Hamiltonian notebook and perform CSE
with open('SEOBNR/SymPy_Hreal_on_bottom.txt', 'w') as output:
    for line in reversed(list(open("SEOBNR/Hamiltonian-Hreal_one_line_expressions.txt"))):
        output.write("%s\n" % line.rstrip())

# Check if a file of partial derivative expressions has already been generated.
# If not, generate them!
# TYLERK: revert Hamiltonian_and_derivs_playground to Hamiltonian_and_derivs after validation is complete
#if not os.path.isfile("SEOBNR_Playground_Pycodes/numpy_expressions.py"):
if not os.path.isfile("SEOBNR_Playground_Pycodes/newdHdx.py"):
    import SEOBNR.Hamiltonian_and_derivs_playground as Had
    Had.output_H_and_derivs()
    # TYLERK: For now, skip CSE (it takes a long time and we just want to validate derivatives)
    #import SEOBNR_Playground_Pycodes.sympy_expression as se
    #se.sympy_cse()
from SEOBNR_Playground_Pycodes.new_dHdx import new_compute_dHdx # For testing
from SEOBNR_Playground_Pycodes.new_dHdp2 import new_compute_dHdp2 # For testing
from SEOBNR_Playground_Pycodes.new_dHdp3 import new_compute_dHdp3 # For testing
from SEOBNR.constant_coeffs import compute_const_coeffs

KK, k0, k1, k2, k3, k4, k5, k5l, dSO, dSS = compute_const_coeffs(eta,EMgamma,a)
#The coefficients do agree with LALSuite!

# Inital root guess
root_guess = [np.divide(1,v*v), v*2, 0.001*200]
# This is the same initial guess given to GSL in LALSuite, but you won't know it unless you're
# careful about their scale factors (which are done and undone and done and undone...)

# Define the function of which we want to find the roots
def root_func(F):
    #Recompute Hamiltonian derivatives using latest minimization guess
    dHdx = new_compute_dHdx(m1, m2, eta, F[0], 0.0, 0.0, 0.0, F[1], F[2], S1hat[0], S1hat[1], S1hat[2], S2hat[0],
                          S2hat[1], S2hat[2], KK, k0, k1, dSO, dSS, 2, EMgamma)
    dHdp2 = new_compute_dHdp2(m1, m2, eta, F[0], 0.0, 0.0, 0.0, F[1], F[2], S1hat[0], S1hat[1], S1hat[2], S2hat[0],
                          S2hat[1], S2hat[2], KK, k0, k1, dSO, dSS, 1, EMgamma)
    dHdp3 = new_compute_dHdp3(m1, m2, eta, F[0], 0.0, 0.0, 0.0, F[1], F[2], S1hat[0], S1hat[1], S1hat[2], S2hat[0],
                          S2hat[1], S2hat[2], KK, k0, k1, dSO, dSS, 1, EMgamma)

    return [dHdx[0]/eta+(-dHdp3[0]*F[2]/eta-dHdp2[0]*F[1]/eta)/F[0], -dHdp3[0]/F[0]/eta, dHdp2[0]/F[0]/eta-omega]

# Find the roots of root_func
soln = root(root_func, root_guess, args=(), method='hybr', jac=None, tol=None, callback=None)
if not(soln.success):
    print("The root finder failed with error message: %s" % soln.message)
    sys.exit(1)

# Populate separation (q) and momentum (p) vectors with the results of root()
q = np.array([soln.x[0], 0., 0.])
p = np.array([0., soln.x[1], soln.x[2]])

Step 3: Rotate $\hat{\bf L} \to {\bf e}_{z}$ [Back to top]¶

$$\label{step3}$$

Note [Back to top]¶

$$\label{step3note}$$

At this point, LALSimIMRSpinEOBInitialConditionsPrec.c normalizes the Cartesian separation and momentum vectors constructed in Step 2. We already have a normalized separation vector $\hat{\bf r}$, so we skip that step.

Step 3.a: Normalize ${\bf q}$ and ${\bf p}$ [Back to top]¶

$$\label{phat}$$

Next we normalize the separation vector ${\bf q}$ and the position vector ${\bf p}$ we found in Step 2:

\begin{align*} \hat{\bf q} &= \frac{ {\bf q} }{ \left\lvert {\bf q} \right\rvert} \\ \hat{\bf p} &= \frac{ {\bf p} }{ \left\lvert {\bf p} \right\rvert}. \end{align*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1101.

In [11]:

# Normalize the separation and momentum vectors
qhat = q/norm(q)
phat = p/norm(p)

Step 3.b: $\hat{\bf L}$ [Back to top]¶

$$\label{lhat}$$

We compute the normalized relativistic angular momentum vector $\hat{\bf L}$:

\begin{equation*} \hat{\bf L} = \frac{ \hat{\bf r} \times \hat{\bf p} }{ \left\lvert \hat{\bf r} \times \hat{\bf p} \right\rvert }. \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Lines 1098--1100.

In [12]:

# Normalize the relativistic angular momentum vector
Lhat = np.cross(rhat,phat)
Lhat /= norm(Lhat)

Step 3.c: Rotation matrix [Back to top]¶

$$\label{rotate}$$

The rotation matrix from the $\left\{ \hat{\bf r}, {\bf v}, \hat{\bf L}_{N} \right\}$ frame to the $\left\{ \hat{\bf r}, {\bf p}, \hat{\bf L} \right\}$ frame is given by

\begin{equation*} \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf r}^{1} & \hat{\bf r}^{2} \\ \hat{\bf p}^{0} & \hat{\bf p}^{1} & \hat{\bf p}^{2} \\ \hat{\bf L}^{0} & \hat{\bf L}^{1} & \hat{\bf L}^{2}\end{bmatrix}. \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1107.

In [13]:

# Rotation matrix
rotate = np.array([rhat, phat, Lhat])

Step 3.d: Rotate $\hat{\bf r}$ [Back to top]¶

$$\label{rotatesrhat}$$

We now rotate $\hat{\bf r}$. We'll use primes to denote the rotated vector.

\begin{equation*} \hat{\bf r}^{\prime} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf r}^{1} & \hat{\bf r}^{2} \\ \hat{\bf p}^{0} & \hat{\bf p}^{1} & \hat{\bf p}^{2} \\ \hat{\bf L}^{0} & \hat{\bf L}^{1} & \hat{\bf L}^{2}\end{bmatrix} \begin{bmatrix} \hat{\bf r}^{0} \\ \hat{\bf r}^{1} \\ \hat{\bf r}^{2} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1112.

In [14]:

# Rotate the normalized separation vector
rhatprm = np.dot(rotate,rhat)

Step 3.e: Rotate $\hat{\bf v}$ [Back to top]¶

$$\label{rotatevhat}$$

We rotate $\hat{\bf v}$. We'll use primes to denote the rotated vector.

\begin{equation*} \hat{\bf v}^{\prime} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf r}^{1} & \hat{\bf r}^{2} \\ \hat{\bf p}^{0} & \hat{\bf p}^{1} & \hat{\bf p}^{2} \\ \hat{\bf L}^{0} & \hat{\bf L}^{1} & \hat{\bf L}^{2}\end{bmatrix} \begin{bmatrix} \hat{\bf v}^{0} \\ \hat{\bf v}^{1} \\ \hat{\bf v}^{2} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1113.

In [15]:

# Rotate the normalized velocity vector
vhatprm = np.dot(rotate, vhat)

Step 3.f: Rotate $\hat{\bf L}_{N}$ [Back to top]¶

$$\label{rotatelnhat}$$

We rotate $\hat{\bf L}_{N}$. We'll use primes to denote the rotated vector.

\begin{equation*} \hat{\bf L}_{N}^{\prime} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf r}^{1} & \hat{\bf r}^{2} \\ \hat{\bf p}^{0} & \hat{\bf p}^{1} & \hat{\bf p}^{2} \\ \hat{\bf L}^{0} & \hat{\bf L}^{1} & \hat{\bf L}^{2}\end{bmatrix} \begin{bmatrix} \hat{\bf L}_{N}^{0} \\ \hat{\bf L}_{N}^{1} \\ \hat{\bf L}_{N}^{2} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1114.

In [16]:

# Rotate the normalized angular momentum vector
LNhatprm = np.dot(rotate, LNhat)

Step 3.g: Rotate ${\bf S}_{1}$ [Back to top]¶

$$\label{rotates1}$$

We rotate ${\bf S}_{1}$. We'll use primes to denote the rotated vector.

\begin{equation*} {\bf S}_{1}^{\prime} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf r}^{1} & \hat{\bf r}^{2} \\ \hat{\bf p}^{0} & \hat{\bf p}^{1} & \hat{\bf p}^{2} \\ \hat{\bf L}^{0} & \hat{\bf L}^{1} & \hat{\bf L}^{2}\end{bmatrix} \begin{bmatrix} {\bf S}_{1}^{0} \\ {\bf S}_{1}^{1} \\ {\bf S}_{1}^{2} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1115.

In [17]:

# Rotate the S1 vector
S1prm = np.dot(rotate, S1)

Step 3.h: Rotate ${\bf S}_{2}$ [Back to top]¶

$$\label{rotates2}$$

We rotate ${\bf S}_{2}$. We'll use primes to denote the rotated vector.

\begin{equation*} {\bf S}_{2}^{\prime} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf r}^{1} & \hat{\bf r}^{2} \\ \hat{\bf p}^{0} & \hat{\bf p}^{1} & \hat{\bf p}^{2} \\ \hat{\bf L}^{0} & \hat{\bf L}^{1} & \hat{\bf L}^{2}\end{bmatrix} \begin{bmatrix} {\bf S}_{2}^{0} \\ {\bf S}_{2}^{1} \\ {\bf S}_{2}^{z} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1116.

In [18]:

# Rotate the S2 vector
S2prm = np.dot(rotate, S2)

Step 3.i: Rotate $\hat{\bf S}_{1}$ [Back to top]¶

$$\label{rotates1hat}$$

We rotate $\hat{\bf S}_{1}$. We'll use primes to denote the rotated vector.

\begin{equation*} \hat{\bf S}_{1}^{\prime} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf r}^{1} & \hat{\bf r}^{2} \\ \hat{\bf p}^{0} & \hat{\bf p}^{1} & \hat{\bf p}^{2} \\ \hat{\bf L}^{0} & \hat{\bf L}^{1} & \hat{\bf L}^{2}\end{bmatrix} \begin{bmatrix} \hat{\bf S}_{1}^{0} \\ \hat{\bf S}_{1}^{1} \\ \hat{\bf S}_{1}^{1} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1117.

In [19]:

# Rotate the normalized S1 vector
S1hatprm = np.dot(rotate, S1hat)

Step 3.j: Rotate $\hat{\bf S}_{2}$ [Back to top]¶

$$\label{rotates2hat\hat}$$

We rotate $\hat{\bf S}_{2}$. We'll use primes to denote the rotated vector.

\begin{equation*} \hat{\bf S}_{2}^{\prime} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf r}^{1} & \hat{\bf r}^{2} \\ \hat{\bf p}^{0} & \hat{\bf p}^{1} & \hat{\bf p}^{2} \\ \hat{\bf L}^{0} & \hat{\bf L}^{1} & \hat{\bf L}^{2}\end{bmatrix} \begin{bmatrix} \hat{\bf S}_{2}^{0} \\ \hat{\bf S}_{2}^{1} \\ \hat{\bf S}_{2}^{2} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1118.

In [20]:

# Rotate the normalized S2 vector
S2hatprm = np.dot(rotate, S2hat)

Step 3.k: Rotate ${\bf q}$ [Back to top]¶

$$\label{rotateq}$$

We rotate ${\bf q}$. We'll use primes to denote the rotated vector.

\begin{equation*} {\bf r}^{\prime} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf r}^{1} & \hat{\bf r}^{2} \\ \hat{\bf p}^{0} & \hat{\bf p}^{1} & \hat{\bf p}^{2} \\ \hat{\bf L}^{0} & \hat{\bf L}^{1} & \hat{\bf L}^{2}\end{bmatrix} \begin{bmatrix} {\bf q}^{0} \\ {\bf q}^{1} \\ {\bf q}^{2} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1119.

In [21]:

# Rotate the separation vector
rprm = np.dot(rotate,q)

Step 3.l: Rotate ${\bf p}$ [Back to top]¶

$$\label{rotatep}$$

We rotate ${\bf p}$. We'll use primes to denote the rotated vector.

\begin{equation*} {\bf p}^{\prime} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf r}^{1} & \hat{\bf r}^{2} \\ \hat{\bf p}^{0} & \hat{\bf p}^{1} & \hat{\bf p}^{2} \\ \hat{\bf L}^{0} & \hat{\bf L}^{1} & \hat{\bf L}^{2}\end{bmatrix} \begin{bmatrix} {\bf p}^{0} \\ {\bf p}^{1} \\ {\bf p}^{2} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1120.

In [22]:

# Rotate the momentum vector
pprm = np.dot(rotate, p)

Step 4: Compute $\dot{\bf r}$ [Back to top]¶

$$\label{step4}$$

Step 4.a: Convert from Cartesian to Spherical Coordinates [Back to top]¶

$$\label{carttosph}$$

We convert position and momentum into spherical coordinates. In the special case where $\theta = \frac{ \pi }{ 2 }$ and $\phi = 0$, the spherical position vector ${\bf r} = \left( {\bf r}^{r}, {\bf r}^{\theta}, {\bf r}^{\phi} \right)$ is given by

\begin{align*} {\bf r}^{r} &= {\bf r}^{0} \\ {\bf r}^{\theta} &= \frac{ \pi }{ 2 } \\ {\bf r}^{\phi} &= 0 \end{align*}

and the spherical momentum vector ${\bf p} = \left( {\bf p}^{r}, {\bf p}^{\theta}, {\bf p}^{\phi} \right)$ is given by

\begin{align*} {\bf p}^{r} &= {\bf p}^{0} \\ {\bf p}^{\theta} &= - {\bf r}^{0}{\bf p}^{2} \\ {\bf p}^{\phi} &= {\bf r}^{0}{\bf p}^{1} \\ \end{align*}

LALSuite calls a Cartesian to spherical routine at LALSimIMRSpinEOBInitialConditionsPrec.c Line 1139, and the function itself is defined on Lines 243--285.

In [23]:

# Convert the separation vector from Cartesian to spherical coordinates
r = np.array([rprm[0], np.divide(np.pi,2.), 0.])
psph = np.array([pprm[0], -rprm[0]*pprm[2], rprm[0]*pprm[1]])

Step 4.b: Second partial derivatives of $H_{\rm real}$ [Back to top]¶

$$\label{seconderiv}$$

In order to include effects of radiation reaction, we need to compute BCD2006 Equation (4.14). This requires that we compute $\frac{ \partial H }{ \partial {\bf p}^{\phi} }$, $\frac{ \partial^{2} H_{\rm real} }{ \partial r^{2} }$, and $\frac{ \partial^{2} H_{\rm real} }{ \partial r \partial {\bf p}^{\phi} }$.

Recall that we are in the special case ${\bf r}^{\theta} = \frac{\pi}{2}$, ${\bf r}^{\phi} = 0$, so that conversion from Cartesian to spherical coordinates is given by the relations in Step 4.a. Then \begin{align*} \partial {\bf r}^{r} &= \partial {\bf r}^{0} \\ \partial {\bf p}^{r} &= \partial {\bf p}^{0} \\ \frac{ \partial {\bf p}^{z} }{ \partial {\bf r}^{r} } &= \frac{ {\bf p}^{\theta} }{ \left( {\bf r}^{r} \right)^{2} } \\ \frac{ \partial {\bf p}^{z} }{ \partial {\bf p}^{\theta} } &= -\frac{ 1 }{ {\bf r}^{r} } \\ \frac{ \partial {\bf p}^{y} }{ \partial {\bf r}^{r} } &= -\frac{ {\bf p}^{\phi} }{ \left( {\bf r}^{r} \right)^{2} } \\ \frac{ \partial {\bf p}^{y} }{ \partial {\bf p}^{\phi} } &= \frac{ 1 }{ {\bf r}^{r} }. \end{align*}

It follows that \begin{equation*} \frac{ \partial H }{ \partial {\bf r}^{r} } = \frac{ \partial H }{ \partial {\bf r}^{x} } \frac{ \partial {\bf r}^{0} }{ \partial {\bf r}^{r} } + \frac{ \partial H }{ \partial {\bf p}^{1} } \frac{ \partial {\bf p}^{1} }{ \partial {\bf r}^{r} } + \frac{ \partial H }{ \partial {\bf p}^{2} } \frac{ \partial {\bf p}^{2} }{ \partial {\bf r}^{r} } = \frac{ \partial H }{ \partial {\bf r}^{0} } - \frac{ \partial H }{ \partial {\bf p}^{1} } \frac{ {\bf p}^{\phi} }{ \left( {\bf r}^{r} \right)^{2} } + \frac{ \partial H }{ \partial {\bf p}^{2} } \frac{ {\bf p}^{\theta} }{ \left( {\bf r}^{r} \right)^{2} }, \end{equation*}

from which it follows that

\begin{align*} \frac{\partial^2 H}{\partial \left( {\bf r}^r \right)^2 } &= \frac{\partial^2 H}{\partial \left( {\bf r}^x \right)^2 } + \frac{\partial^2 H}{\partial \left( {\bf p}^z \right)^2 } \cdot \frac{ \left( {\bf p^\theta } \right)^2 }{ \left( {\bf r}^r \right)^4 } - \frac{\partial^2 H}{\partial {\bf p}^y \partial {\bf p}^z } \cdot \frac{ 2 {\bf p}^\theta * {\bf p}^\phi }{ {\bf r}^4 } + \frac{\partial^2 H}{\partial \left( {\bf p}^y \right)^2 } \cdot \frac{ \left( {\bf p^\phi } \right)^2 }{ \left( {\bf r}^r \right)^4 } + \frac{\partial^2 H}{\partial {\bf r}^x \partial {\bf p}^y } \cdot \frac{ 2 {\bf p}^\phi }{ {\bf r}^2 } \\ \frac{\partial^2 H}{\partial {\bf r}^r \partial {\bf p}^\phi } &= \frac{\partial^2 H}{\partial {\bf p}^y \partial {\bf p}^z } \cdot \frac{ {\bf p}^\theta }{ \left( {\bf r}^r \right)^3 } - \frac{\partial^2 H}{\partial \left( {\bf p}^y \right)^2 } \cdot \frac{ {\bf p}^\phi }{ \left( {\bf r}^r \right)^3 } + \frac{\partial^2 H}{\partial {\bf r}^x \partial {\bf p}^y } \cdot \frac{ 1 }{ {\bf r}^r } \end{align*}

Note: be sure that, following this, we use normalized spins.

In [24]:

## Import second derivatives of H from another function/routine
#if not os.path.isfile("SEOBNR_Playground_Pycodes/d2Hdx2.py"):
#    print("I'M BEING BAD")
##    import SEOBNR.Hamiltonian_second_derivs_playground as Hsd
##    Hsd.output_H_sec_derivs()
#    import SEOBNR.Hamiltonian_second_derivs_x_playground as Hsd
#    Hsd.output_H_sec_derivs()
#from SEOBNR_Playground_Pycodes.d2Hdx2 import compute_d2Hdx2 # For testing

##d2Hdx2 = compute_d2Hdx2(m1, m2, eta, 2.129681018601393e+01, 0.0, 0.0,
##                        0.0, 2.335391115580442e-01, -4.235164736271502e-22,
##                        S1[0], S1[1], S1[2],
##                        S2[0], S2[1], S2[2], KK, k0, k1, dSO, dSS, 2, EMgamma)
##d2Hdx2 = compute_d2Hdx2(m1, m2, eta, 2.129681018601393e+01, 0.0, 0.0,
##                        0.0, 2.335391115580442e-01, -4.235164736271502e-22,
##                        S1hat[0], S1hat[1], S1hat[2],
##                        S2hat[0], S2hat[1], S2hat[2], KK, k0, k1, dSO, dSS, 0, EMgamma)
#d2Hdx2 = compute_d2Hdx2(m1, m2, eta, 2.129680018601393e+01, 0.000000000000000e+00, 0.000000000000000e+00,
#                        0.000000000000000e+00, 2.335392212172933e-01, -4.235166724910499e-22,
#                        4.857667584940312e-03, 9.715161660389764e-03, -1.457311842632286e-02,
#                        3.673094582185491e-03, -4.591302628615413e-03, 5.509696538546906e-03, KK, k0, k1, dSO, dSS, 1, EMgamma)
#dHdx = new_compute_dHdx(m1, m2, eta, 2.129680018601393e+01, 0.000000000000000e+00, 0.000000000000000e+00,
#                        0.000000000000000e+00, 2.335392212172933e-01, -4.235166724910499e-22,
#                        4.857667584940312e-03, 9.715161660389764e-03, -1.457311842632286e-02,
#                        3.673094582185491e-03, -4.591302628615413e-03, 5.509696538546906e-03, KK, k0, k1, dSO, dSS, 1, EMgamma)
#dHdpy = new_compute_dHdp2(m1, m2, eta, 2.129680018601393e+01, 0.000000000000000e+00, 0.000000000000000e+00,
#                        0.000000000000000e+00, 2.335392212172933e-01, -4.235166724910499e-22,
#                        4.857667584940312e-03, 9.715161660389764e-03, -1.457311842632286e-02,
#                        3.673094582185491e-03, -4.591302628615413e-03, 5.509696538546906e-03, KK, k0, k1, dSO, dSS, 1, EMgamma)
#dHdpz = new_compute_dHdp3(m1, m2, eta, 2.129680018601393e+01, 0.000000000000000e+00, 0.000000000000000e+00,
#                        0.000000000000000e+00, 2.335392212172933e-01, -4.235166724910499e-22,
#                        4.857667584940312e-03, 9.715161660389764e-03, -1.457311842632286e-02,
#                        3.673094582185491e-03, -4.591302628615413e-03, 5.509696538546906e-03, KK, k0, k1, dSO, dSS, 1, EMgamma)

#print("dHdx = %.15e" % (dHdx/eta))
#print("dHdpy = %.15e" % (dHdpy/eta))
#print("dHdpz = %.15e" % (dHdpz/eta))
#print("r = %.15e" % r[0])
#print("ptheta = %.15e" % psph[1])
#print("pphi = %.15e" % psph[2])
#print("dHdr is %.15e" % (dHdx - dHdpy*psph[2]/r[0]/r[0] + dHdpz*psph[1]/r[0]/r[0]))
#print("dHdr term2 is %.15e" % (dHdpy*psph[2]/r[0]/r[0]))
#print("dHdr term3 is %.15e" % (dHdpz*psph[1]/r[0]/r[0]))
#print("d2Hdx2 = %.15e" % (d2Hdx2/eta))

#sys.exit(1) # Might as well exit here for now.

Step 4.c: $\frac{ \partial E }{ \partial r }$ [Back to top]¶

$$\label{dedr}$$

We seek to compute $\frac{ \partial H }{\partial r}$, and BCD2006 uses the convention $H \equiv E$. (see BCD2006 Equation (3.7)). From BCD2006 Equation Equation (4.14) (noting that this equation applies in spherical coordinates when ${\bf r}$ is directed along the ${\bf e}_{0}$ axis),

\begin{equation*} \frac{ \partial E }{ \partial r } = -\frac{ \frac{ \partial H }{ \partial {\bf p}^{\phi} } \frac{ \partial^{2} H }{ \left(\partial {\bf r}^{r} \right)^{2} } }{ \frac{ \partial^{2} H }{ \partial {\bf r}^{r} \partial {\bf p}^{\phi} } }. \end{equation*}

In [25]:

## Time derivative of Hamiltonain with respect to separation magnitude r
#dEdr = -dHdpphi*d2Hdr2/d2Hdrdpphi

Step 4.e: $\boldsymbol{\sigma}^{*}$ [Back to top]¶

$$\label{sigmastar}$$

From BB2010 Equation (5.3),

\begin{equation*} \boldsymbol{\sigma}^{*} = \frac{ m_{2} }{ m_{1} } {\bf S}_{1} + \frac{ m_{1} }{ m_{2} }{\bf S}_{2}. \end{equation*}

In [26]:

## Spin combination sigmastar
#sigmastar = np.add(np.divide(m2,m1)*S1, np.divide(m1,m2)*S2)

Step 4.f: $H_{\rm real}$ [Back to top]¶

$$\label{hreal}$$

We now compute $H_{\rm real}$ (LALSimIMRSpinEOBInitialConditionsPrec.c Line 1217). To do so, we need to restructure the output of Tutorial-SEOBNR_Documentation by first making sure each expression is on a single line and then reversing the lines.

In [27]:

#import SEOBNR_Playground_Pycodes.Hreal_on_bottom as Ham
##All inputs agree with LALSuite
##eta, KK, tortoise, and dSO all agree with LALSuite to 16 significant digits
##Hard-code other inputs so we know they agree exactly with LALSuite
##LALSuite command used: ./lalsimulation/src/lalsim-inspiral -a SEOBNRv3 -M 23 -m 10 -f 20 -X 0.01 -Y 0.02 -Z -0.03 -x 0.04 -y -0.05 -z 0.06

#Hreal = Ham.compute_Hreal(m1=m1, m2=m2, EMgamma=EMgamma, tortoise=1, dSO=dSO, dSS=dSS,
#                          x=rprm[0], y=rprm[1], z=rprm[2], p1=pprm[0], p2=pprm[1], p3=pprm[2],
#                          S1x=S1hatprm[0], S1y=S1hatprm[1], S1z=S1hatprm[2],
#                          S2x=S2hatprm[0], S2y=S2hatprm[1], S2z=S2hatprm[2])
#print(Hreal)#TylerK
##Hreal = Ham.compute_Hreal(m1, m2, EMgamma, 1, dSO, dSS,
##                          2.129681018601393e+01, 0.000000000000000e+00, 0.000000000000000e+00,
##                          0.000000000000000e+00, 2.335391115580442e-01, -4.235164736271502e-22,
##                          4.857667584940312e-03, 9.715161660389764e-03, -1.457311842632286e-02,
##                          3.673094582185491e-03, -4.591302628615413e-03, 5.509696538546906e-03)

In [28]:

##Temporary validation code block: all hard-coded values from LALSuite!
#Hreal_valid = Ham.compute_Hreal(m1=23., m2=10., EMgamma=EMgamma, tortoise=1,
#                                        dSO=-7.966696593617955e+01, dSS=1.261873764525631e+01,
#                          x=2.129681018601393e+01, y=0.000000000000000e+00, z=0.000000000000000e+00,
#                          p1=0.000000000000000e+00, p2=2.335391115580442e-01, p3=-4.235164736271502e-22,
#                          S1x=4.857667584940312e-03, S1y=9.715161660389764e-03, S1z=-1.457311842632286e-02,
#                          S2x=3.673094582185491e-03, S2y=-4.591302628615413e-03, S2z=5.509696538546906e-03)
#print(Hreal_valid)#TylerK
#Hreal_valid = Ham.compute_Hreal()
#print(Hreal_valid)#TylerK
#if(np.abs(Hreal_valid-9.952429072947245e-01)>1e-14):
#    print("ERROR. You have broken the Hamiltonian computation!")
#    sys.exit(1)

Polar data [Back to top]¶

$$\label{polardata}$$

At LALSimIMRSpinEOBInitialConditionsPrec.c Lines 1234--1238, we set the following polar data ${\bf P}$:

\begin{align*} {\bf P}^{0} &= {\bf r}^{r} \\ {\bf P}^{1} &= 0 \\ {\bf P}^{2} &= {\bf p}^{r} \\ {\bf P}^{3} &= {\bf p}^{\phi} \end{align*}

In [29]:

## Populate a vector of polar coordinate values
#polar = np.array([r[0], 0., psph[0], psph[1]])

vPhiKepler [Back to top]¶

$$\label{vphikepler}$$

From T2012 Equation (A2),

\begin{equation*} {\rm vPhiKepler} = \frac{ 1 }{ \omega^{2} \left( {\bf r}^{r} \right)^{3} }. \end{equation*}

See LALSimIMRSpinEOBFactorizedWaveformPrec_v3opt.c Lines 113 and 1271--1315. Note that SEOBNRv3_opt recalculates $\omega$, but I think the $\omega$ above is based on a circular orbit and therefore the recalcuation is unnecessary.

In [30]:

## Keplarian velocity
#vPhiKepler = 1./(omega*omega*r[0]*r[0]*r[0])

${\bf r} \times {\bf p}$ [Back to top]¶

$$\label{rcrossp}$$

We'll use the notation

\begin{equation*} {\rm rcrossp} = {\bf r}^{\prime} \times {\bf p}^{\prime}. \end{equation*}

See LALSimIMRSpinEOBFactorizedWaveformPrec_v3opt.c Lines 170--172.

In [31]:

## r cross p
#rcrossp = np.cross(rprm,pprm)

vPhi [Back to top]¶

$$\label{vphi}$$

We'll use the notation ( paper reference?)

\begin{equation*} {\rm vPhi} = \omega {\bf r}^{r} \sqrt[3]{\rm vPhiKepler}. \end{equation*}

See LALSimIMRSpinEOBFactorizedWaveformPrec_v3opt.c Lines 185 and 190.

In [32]:

## Keplarian velocity
#vPhi = omega*r[0]*np.cbrt(vPhiKepler)

${\bf S}_{i} \cdot \hat{\bf L}$ [Back to top]¶

$$\label{sidotl}$$

We compute ${\bf S}_{1} \cdot \hat{\bf L}$ and ${\bf S}_{2} \cdot \hat{\bf L}$.

See LALSimIMRSpinEOBFactorizedFluxPrec_v3opt.c lines 131--134.

In [33]:

## S dot L
#s1dotL = np.dot(S1,Lhat)
#s2dotL = np.dot(S2,Lhat)

$\boldsymbol{\chi}_{\rm S}$, $\boldsymbol{\chi}_{\rm A}$ [Back to top]¶

$$\label{chii}$$

From P2014 Equation 17, we have

\begin{align*} \chi_{\rm S} = \frac{1}{2} \left( {\bf S}_{1} + {\bf S}_{2} \right) \cdot \hat{\bf L} \\ \chi_{\rm A} = \frac{1}{2} \left( {\bf S}_{1} - {\bf S}_{2} \right) \cdot \hat{\bf L} \end{align*}

In [34]:

## Spin combinations chiS and chiA
#chiS = 0.5*(s1dotL + s2dotL)
#chiA = 0.5*(s1dotL - s2dotL)

$\hat{m}_{i}$ [Back to top]¶

$$\label{mihat}$$

We scale the masses $m_{1}$, $m_{2}$ by total mass. See LALSimIMREOBNewtonianMultipole.c Lines 540--541.

\begin{align*} \hat{m}_{1} = \frac{ m_{1} }{ M } \\ \hat{m}_{2} = \frac{ m_{2} }{ M } \\ \end{align*}

In [35]:

## Normalized mass
#mhat1 = m1*Minv
#mhat2 = m2*Minv

Newtonian multipole [Back to top]¶

$$\label{newtonianmultipole}$$

The Newtonian multipolar waveform is given in DIN2009 Equation (4). For a given $(\ell, m)$ we define

\begin{align*} \epsilon &= \left( \ell + m \right) {\rm mod } 2 \\ n &= \left( i m \right)^{\ell} \frac{ 8 \pi }{ \left( 2 \ell + 1 \right)!! } \sqrt{ \frac{ \left( \ell + 1 \right) \left( \ell + 2 \right) }{ \ell \left( \ell - 1 \right) } } \end{align*}

along with the associated Legendre function evaluated at zero. See LALSimIMRSpinEOBFactorizedWaveformPrec_v3opt.c Line 206 and LALSimIMREOBNewtonianMultipole.c Lines 205, 210, 290, and 309--506.

In [36]:

## Compute Newtonian multipole
## Compute the associated Legendre function of degree l and order m at x=0
#def AssociatedLegendre(l,m):
#    if l==1:
#        if m==1:
#            return -1.
#        else:
#            print("You used a bad (l,m)")
#    if l==2:
#        if m==2:
#            return 3.
#        elif m==1:
#            return 0.
#        else:
#            print("You used a bad (l,m)")
#    if l==3:
#        if m==3:
#            return 15.
#        elif m==2:
#            return 0.
#        elif m==1:
#            return 1.5
#        else:
#            print("You used a bad (l,m)")
#    if l==4:
#        if m==4:
#            return 105.
#        elif m==3:
#            return 0.
#        elif m==2:
#            return -7.5
#        elif m==1:
#            return 0.
#        else:
#            print("You used a bad (l,m)")
#    if l==5:
#        if m==5:
#            return -945.
#        elif m==4:
#            return 0.
#        elif m==3:
#            return 52.5
#        elif m==2:
#            return 0.
#        elif m==1:
#            return -1.875
#        else:
#            print("You used a bad (l,m)")
#    if l==6:
#        if m==6:
#            return 10395.
#        elif m==5:
#            return 0.
#        elif m==4:
#            return -472.5
#        elif m==3:
#            return 0.
#        elif m==2:
#            return 13.125
#        elif m==1:
#            return 0.
#        else:
#            print("You used a bad (l,m)")
#    if l==7:
#        if m==7:
#            return -135135.
#        elif m==6:
#            return 0.
#        elif m==5:
#            return 5197.5
#        elif m==4:
#            return 0.
#        elif m==3:
#            return -118.125
#        elif m==2:
#            return 0.
#        elif m==1:
#            return 2.1875
#        else:
#            print("You used a bad (l,m)")
#    if l==8:
#        if m==8:
#            return 2027025.
#        elif m==7:
#            return 0.
#        elif m==6:
#            return -67567.5
#        elif m==5:
#            return 0.
#        elif m==4:
#            return 1299.375
#        elif m==3:
#            return 0.
#        elif m==2:
#            return -19.6875
#        elif m==1:
#            return 0.
#        else:
#            print("You used a bad (l,m)")

## Compute the prefix for the Newtonian multipole
#def NewtonianPrefix(m1,m2,l,m,epsilon,eta):
#    Mtot = m1 + m2
#    m1hat = np.divide(m1,Mtot)
#    m2hat = np.divide(m2,Mtot)
#    if (m%2)==0:
#        sign = 1
#    else:
#        sign = -1
#    lpepm1 = l + epsilon - 1
#    if (m1!=m2) or sign==1:
#        c = np.power(m2hat,lpepm1) + sign*np.power(m1hat,lpepm1)
#    else:
#        if l==2 or l==3:
#            c = -1.
#        elif l==4 or l==5:
#            c = -0.5
#        else:
#            c = 0.
#    n = np.power(complex(0,m), l)
#    doubfact = doublefactorial(2*l+1)
#    if epsilon==0:
#        n *= 8.*np.divide(np.pi,doubfact)
#        n *= np.sqrt(np.divide((l+1)*(l+2),l*(l-1)))
#    elif epsilon==1:
#        n = -n
#        n *= 16.j*np.divide(np.pi,doubfact)
#        n *= np.sqrt( np.divide((2*l+1)* (l+2) * (l*l - m*m),(2*l - 1) * (l+1) * l * (l-1)) )
#    else:
#        print("Epsilon must be 0 or 1")
#        exit()
#    return n*eta*c

## Function to compute a double factorial; see https://en.wikipedia.org/wiki/Double_factorial
#def doublefactorial(n):
#     if n <= 0:
#         return 1
#     else:
#         return n * doublefactorial(n-2)

hLMTab [Back to top]¶

$$\label{hlmtab}$$

In order to compute flux, we need to build the matrix "hLMTab". See T2012 Equation (17) and the Appendix, along with this private LIGO doc.

In [37]:

## The following populates a matrix T_{lm} of resummed leading-order logarithms of tail effects
#deltam = np.divide(m1 - m2,m1 + m2)
#flux = 0.

#fa1 = interp1d(nqi.domain, nqi.a1Range, kind='cubic')
#fa2 = interp1d(nqi.domain, nqi.a2Range, kind='cubic')
#fa3 = interp1d(nqi.domain, nqi.a3Range, kind='cubic')
#fb1 = interp1d(nqi.domain, nqi.b1Range, kind='cubic')
#fb2 = interp1d(nqi.domain, nqi.b2Range, kind='cubic')

#a1 = fa1(eta)
#a2 = fa2(eta)
#a3 = fa3(eta)
#b1 = -fb1(eta)
#b2 = -fb2(eta)

#fa3sAmax = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Amaxa3sVal, kind='cubic')
#fa4Amax = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Amaxa4Val, kind='cubic')
#fa5Amax = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Amaxa5Val, kind='cubic')
#fb3Amax = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Amaxb3Val, kind='cubic')
#fb4Amax = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Amaxb4Val, kind='cubic')

#fa3sAmed = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Ameda3sVal, kind='cubic')
#fa4Amed = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Ameda4Val, kind='cubic')
#fa5Amed = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Ameda5Val, kind='cubic')
#fb3Amed = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Amedb3Val, kind='cubic')
#fb4Amed = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Amedb4Val, kind='cubic')

#fa3sAmin = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Amina3sVal, kind='cubic')
#fa4Amin = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Amina4Val, kind='cubic')
#fa5Amin = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Amina5Val, kind='cubic')
#fb3Amin = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Aminb3Val, kind='cubic')
#fb4Amin = interp2d(nqi.aDomain, nqi.etaDomain, nqi.Aminb4Val, kind='cubic')

#chiAmaxCoeffs = [fa3sAmax(a,eta), fa4Amax(a,eta), fa5Amax(a,eta), fb3Amax(a,eta), fb4Amax(a,eta)]
#chiAmedCoeffs = [fa3sAmed(a,eta), fa4Amed(a,eta), fa5Amed(a,eta), fb3Amed(a,eta), fb4Amed(a,eta)]
#chiAminCoeffs = [fa3sAmin(a,eta), fa4Amin(a,eta), fa5Amin(a,eta), fb3Amin(a,eta), fb4Amin(a,eta)]

#chi = a/(1. - 2.*eta)

#if eta < 1.0e-15:
#    chiAmax = np.divide(chi + 1.,2.)
#    chiAmin = np.divide(chi - 1.,2.)
#else:
#    if chi <= 0:
#        chiAmax = (1. + chi)*(1. - 2.*eta)/(1.+ deltam - 2.*eta)
#        if (1. + deltam - 2.*eta + 2.*chi*(1. - 2.*eta))/(1. - deltam - 2.*eta) < 1.:
#            chiAmin = -(1. + chi)*(1. - 2.*eta)/(1. - deltam - 2.*eta)
#        else:
#            chiAmin = -(1. - chi)*(1. - 2.*eta)/(1. + deltam - 2.*eta)
#    else:
#        chiAmin = -(1. - chi)*(1. - 2.*eta)/(1. + deltam - 2.*eta)
#        if -(1. + deltam - 2.*eta - 2.*chi*(1. - 2.*eta))/(1. - deltam - 2.*eta) > -1.:
#            chiAmax = (1. - chi)*(1. - 2.*eta)/(1. - deltam - 2.*eta)
#        else:
#            chiAmax = (1. + chi)*(1. - 2.*eta)/(1. + deltam - 2.*eta)

#chiAmed = np.divide(chiAmax + chiAmin,2.)
#if chiAmax < 1.0e-15:
#    cmax = 1.0
#    cmed = 0.0
#    cmin = 0.0
#else:
#    cmax =  (chiA - chiAmed)*(chiA - chiAmin)/(chiAmax - chiAmed)/(chiAmax - chiAmin)
#    cmed = -(chiA - chiAmax)*(chiA - chiAmin)/(chiAmax - chiAmed)/(chiAmed - chiAmin)
#    cmin =  (chiA - chiAmax)*(chiA - chiAmed)/(chiAmax - chiAmin)/(chiAmed - chiAmin)

#nqcmax = chiAmaxCoeffs[0]
#nqcmed = chiAmedCoeffs[0]
#nqcmin = chiAminCoeffs[0]
#a3S = cmax*nqcmax + cmed*nqcmed + cmin*nqcmin

#nqcmax = chiAmaxCoeffs[1]
#nqcmed = chiAmedCoeffs[1]
#nqcmin = chiAminCoeffs[1]
#a4  = cmax*nqcmax + cmed*nqcmed + cmin*nqcmin

#nqcmax = chiAmaxCoeffs[2]
#nqcmed = chiAmedCoeffs[2]
#nqcmin = chiAminCoeffs[2]
#a5  = cmax*nqcmax + cmed*nqcmed + cmin*nqcmin

#nqcmax = chiAmaxCoeffs[3]
#nqcmed = chiAmedCoeffs[3]
#nqcmin = chiAminCoeffs[3]
#b3  = cmax*nqcmax + cmed*nqcmed + cmin*nqcmin
#nqcmax = chiAmaxCoeffs[4]
#nqcmed = chiAmedCoeffs[4]
#nqcmin = chiAminCoeffs[4]
#b4  = cmax*nqcmax + cmed*nqcmed + cmin*nqcmin

#rsq = polar[0]*polar[0]
#sqrtr = np.sqrt(polar[0])
#prsq = polar[2]*polar[2]

#mag = 1. + (prsq/(rsq*omega*omega))*(a1 + a2/polar[0] + (a3 + a3S)/(polar[0]*sqrtr) + a4/rsq + a5/(rsq*sqrtr))
#phase = b1*polar[2]/(polar[0]*omega) + prsq*polar[2]/(polar[0]*omega)*(b2 + b3/sqrtr + b4/polar[0])
#nqc = complex(mag*np.cos(phase),0)
#nqc += complex(0,mag*np.sin(phase))

#import factorized_modes as fm
#for l in range(2, 9):
#    for m in range(1, l+1):
#        epsilon = (l + m) % 2
#        legendre = AssociatedLegendre(l-epsilon,m)*np.sqrt(2*l+1*np.divide(np.math.factorial(l-m),4)*np.pi*np.math.factorial(l+m))
#        #Note that LALSimIMREOBNewtonianMultipole.c Line 74 atrributes the
#        #Newtonian prefix calculations to https://arxiv.org/abs/1106.1021v2
#        prefix = NewtonianPrefix(m1,m2,l,m,epsilon,eta)
#        multipole = prefix*legendre*np.power(vPhi*vPhi,(l+epsilon)/2.)
#        if ((l+m)%2)==0:
#            Slm = (Hreal*Hreal - 1.)/(2.*eta) + 1.
#        else:
#            Slm = v*psph[2]
#        eulerlog = EMgamma + np.log(2.*m*v)
#        k = m*omega
#        Hrealk = Hreal * k
#        Hrealksq4 = 4. * Hrealk*Hrealk
#        Hrealk4pi = 4. * np.pi *Hrealk
#        Tlmprefac = np.sqrt(Hrealk4pi/(1.-np.exp(-Hrealk4pi)))/np.math.factorial(l)
#        Tlmprodfac = 1.
#        for i in range(1,l+1):
#            Tlmprodfac *= Hrealksq4 + (i*i)
#        Tlm = Tlmprefac*np.sqrt(Tlmprodfac)
#        auxflm = 0.
#        if l==2:
#            if m==2:
#                rholm = 1 + vsq * (fm.rho22v2 + v*(fm.rho22v3 + v*(fm.rho22v4 + v*(fm.rho22v5 + v*(fm.rho22v6
#                                + fm.rho22v6l*eulerlog + v*(fm.rho22v7 + v*(fm.rho22v8 + fm.rho22v8l*eulerlog
#                                + (fm.rho22v10 + fm.rho22v10l*eulerlog)*vsq)))))))
#            elif m==1:
#                rholm = 1. + v * (fm.rho21v1 + v*(fm.rho21v2 + v*(fm.rho21v3 + v*(fm.rho21v4 + v*(fm.rho21v5
#                                + v*(fm.rho21v6 + fm.rho21v6l*eulerlog + v*(fm.rho21v7 + fm.rho21v7l*eulerlog
#                                + v*(fm.rho21v8 + fm.rho21v8l*eulerlog + (fm.rho21v10 + fm.rho21v10l*eulerlog)*vsq))))))))
#                auxflm = v*fm.f21v1 + vsq*v*fm.f21v3
#            else:
#                print("You used a bad (l,m)")
#        elif l==3:
#            if m==3:
#                rholm = 1. + vsq*(fm.rho33v2 + v*(fm.rho33v3 + v*(fm.rho33v4 + v*(fm.rho33v5 + v*(fm.rho33v6
#                                + fm.rho33v6l*eulerlog + v*(fm.rho33v7 + (fm.rho33v8 + fm.rho33v8l*eulerlog)*v))))))
#                auxflm = v*vsq*fm.f33v3;
#            elif m==2:
#                rholm = 1. + v*(fm.rho32v + v*(fm.rho32v2 + v*(fm.rho32v3 + v*(fm.rho32v4 + v*(fm.rho32v5
#                                + v*(fm.rho32v6 + fm.rho32v6l*eulerlog + (fm.rho32v8 + fm.rho32v8l*eulerlog)*vsq))))))
#            elif m==1:
#                rholm = 1. + vsq*(fm.rho31v2 + v*(fm.rho31v3 + v*(fm.rho31v4 + v*(fm.rho31v5 + v*(fm.rho31v6
#                                + fm.rho31v6l*eulerlog + v*(fm.rho31v7 + (fm.rho31v8 + fm.rho31v8l*eulerlog)*v))))))
#                auxflm = v*vsq*fm.f31v3
#            else:
#                print("You used a bad (l,m)")
#        elif l==4:
#            if m==4:
#                rholm = 1. + vsq*(fm.rho44v2 + v*(fm.rho44v3 + v*(fm.rho44v4 + v*(fm.rho44v5 + (fm.rho44v6
#                                + fm.rho44v6l*eulerlog)*v))))
#            elif m==3:
#                rholm = 1. + v*(fm.rho43v + v*(fm.rho43v2 + vsq*(fm.rho43v4 + v*(fm.rho43v5 + (fm.rho43v6
#                                + fm.rho43v6l*eulerlog)*v))))
#                auxflm = v*fm.f43v
#            elif m==2:
#                rholm = 1. + vsq*(fm.rho42v2 + v*(fm.rho42v3 + v*(fm.rho42v4 + v*(fm.rho42v5 + (fm.rho42v6
#                                + fm.rho42v6l*eulerlog)*v))))
#            elif m==1:
#                rholm = 1. + v*(fm.rho41v + v*(fm.rho41v2 + vsq*(fm.rho41v4 + v*(fm.rho41v5 + (fm.rho41v6
#                                + fm.rho41v6l*eulerlog)*v))))
#                auxflm = v*fm.f41v
#            else:
#                print("You used a bad (l,m)")
#        elif l==5:
#            if m==5:
#                rholm = 1. + vsq*(fm.rho55v2 + v*(fm.rho55v3 + v*(fm.rho55v4 + v*(fm.rho55v5 + fm.rho55v6*v))))
#            elif m==4:
#                rholm = 1. + vsq*(fm.rho54v2 + v*(fm.rho54v3 + fm.rho54v4*v))
#            elif m==3:
#                rholm = 1. + vsq*(fm.rho53v2 + v*(fm.rho53v3 + v*(fm.rho53v4 + fm.rho53v5*v)))
#            elif m==2:
#                rholm = 1. + vsq*(fm.rho52v2 + v*(fm.rho52v3 + fm.rho52v4*v))
#            elif m==1:
#                rholm = 1. + vsq*(fm.rho51v2 + v*(fm.rho51v3 + v*(fm.rho51v4 + fm.rho51v5*v)))
#            else:
#                print("You used a bad (l,m)")
#        elif l==6:
#            if m==6:
#                rholm = 1. + vsq*(fm.rho66v2 + v*(fm.rho66v3 + fm.rho66v4*v))
#            elif m==5:
#                rholm = 1. + vsq*(fm.rho65v2 + fm.rho65v3*v)
#            elif m==4:
#                rholm = 1. + vsq*(fm.rho64v2 + v*(fm.rho64v3 + fm.rho64v4*v))
#            elif m==3:
#                rholm = 1. + vsq*(fm.rho63v2 + fm.rho63v3*v)
#            elif m==2:
#                rholm = 1. + vsq*(fm.rho62v2 + v*(fm.rho62v3 + fm.rho62v4*v))
#            elif m==1:
#                rholm = 1. + vsq*(fm.rho61v2 + fm.rho61v3*v)
#            else:
#                print("You used a bad (l,m)")
#        elif l==7:
#            if m==7:
#                rholm = 1. + vsq*(fm.rho77v2 + fm.rho77v3*v)
#            elif m==6:
#                rholm = 1. + fm.rho76v2*vsq
#            elif m==5:
#                rholm = 1. + vsq*(fm.rho75v2 + fm.rho75v3*v)
#            elif m==4:
#                rholm = 1. + fm.rho74v2*vsq
#            elif m==3:
#                rholm = 1. + vsq*(fm.rho73v2 + fm.rho73v3*v)
#            elif m==2:
#                rholm = 1. + fm.rho72v2*vsq
#            elif m==1:
#                rholm = 1. + vsq*(fm.rho71v2 + fm.rho71v3*v)
#            else:
#                print("You used a bad (l,m)")
#        elif l==8:
#            if m==8:
#                rholm = 1. + fm.rho88v2*vsq
#            elif m==7:
#                rholm = 1. + fm.rho87v2*vsq
#            elif m==6:
#                rholm = 1. + fm.rho86v2*vsq
#            elif m==5:
#                rholm = 1. + fm.rho85v2*vsq
#            elif m==4:
#                rholm = 1. + fm.rho84v2*vsq
#            elif m==3:
#                rholm = 1. + fm.rho83v2*vsq
#            elif m==2:
#                rholm = 1. + fm.rho82v2*vsq
#            elif m==1:
#                rholm = 1. + fm.rho81v2*vsq
#            else:
#                print("You used a bad (l,m)")
#        else:
#            print("You used a bad (l,m)")
#        rholmPowl = np.power(rholm,l)
#        if eta==0.25 and (m % 2):
#            rholmPowl = auxflm
#        else:
#            rholmPowl += auxflm
#        hlm = Tlm*Slm*rholmPowl*multipole
#        if (m*m*omega*omega*hlm*hlm) > 5.:
#            hlm *= nqc
#        flux += m*m*omega*omega*hlm*hlm
#        if omega*omega > 1 or flux > 5:
#            flux = 0.
#        flux *= np.divide(8.,np.pi)
#flux /= eta
#rdot = -flux/dEdr

#pr = rdot/(dHdpr/px)

Step 5: Invert the rotation of Step 3 [Back to top]¶

$$\label{step5}$$

Inverse Rotation Matrix [Back to top]¶

$$\label{invrotationmatrix}$$

The matrix to invert the rotation applied in Step 3 is:

\begin{equation*} \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf p}^{0} & \hat{\bf L}^{0} \\ \hat{\bf r}^{1} & \hat{\bf p}^{1} & \hat{\bf L}^{1} \\ \hat{\bf r}^{2} & \hat{\bf p}^{2} & \hat{\bf L}^{2}\end{bmatrix}. \end{equation*}

To see that this is indeed the correct matrix inverse, note that by construction $\hat{\bf q}$, $\hat{\bf p}$, and $\hat{\bf L}$ are all unit vectors orthogonal to one another. See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1107.

In [38]:

#invert00 = rhat0
#invert01 = phat0
#invert02 = Lhat0
#invert10 = rhat1
#invert11 = phat1
#invert12 = Lhat1
#invert20 = rhat2
#invert21 = phat2
#invert22 = Lhat2

Rotate $\hat{\bf r}^{\prime}$ [Back to top]¶

$$\label{invrotaterhat}$$

We rotate $\hat{\bf r}^{\prime}$ and call the new separation vector ${\bf r}$.

\begin{equation*} \hat{\bf r} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf p}^{0} & \hat{\bf L}^{0} \\ \hat{\bf r}^{1} & \hat{\bf p}^{1} & \hat{\bf L}^{1} \\ \hat{\bf r}^{2} & \hat{\bf p}^{2} & \hat{\bf L}^{2} \end{bmatrix} \begin{bmatrix} \hat{\bf r}^{\prime 0} \\ \hat{\bf r}^{\prime 1} \\ \hat{\bf r}^{\prime 2} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1315.

In [39]:

#rhat0 = rhat0*rhatprm0 + phat0*rhatprm1 + Lhat0*rhatprm2
#rhat1 = rhat1*rhatprm0 + phat1*rhatprm1 + Lhat1*rhatprm2
#rhat0 = rhat2*rhatprm0 + phat2*rhatprm1 + Lhat2*rhatprm2

Rotate $\hat{\bf v}^{\prime}$ [Back to top]¶

$$\label{invrotatevhat}$$

We rotate $\hat{\bf v}^{\prime}$ and call the new separation vector ${\bf v}$.

\begin{equation*} \hat{\bf v} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf p}^{0} & \hat{\bf L}^{0} \\ \hat{\bf r}^{1} & \hat{\bf p}^{1} & \hat{\bf L}^{1} \\ \hat{\bf r}^{2} & \hat{\bf p}^{2} & \hat{\bf L}^{2} \end{bmatrix} \begin{bmatrix} \hat{\bf v}^{\prime 0} \\ \hat{\bf v}^{\prime 1} \\ \hat{\bf v}^{\prime 2} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1316.

In [40]:

#vhat0 = rhat0*vhatprm0 + phat0*vhatprm1 + Lhat0*vhatprm2
#vhat1 = rhat1*vhatprm0 + phat1*vhatprm1 + Lhat1*vhatprm2
#vhat2 = rhat2*vhatprm0 + phat2*vhatprm1 + Lhat2*vhatprm2

Rotate $\hat{\bf L}_{N}^{\prime}$ [Back to top]¶

$$\label{invrotatelnhat}$$

We rotate $\hat{\bf L}_{N}^{\prime}$ and call the new separation vector ${\bf L}_{N}$.

\begin{equation*} \hat{\bf L}_{N} = \begin{bmatrix} \hat{\bf r}^{0} & \hat{\bf p}^{0} & \hat{\bf L}^{0} \\ \hat{\bf r}^{1} & \hat{\bf p}^{1} & \hat{\bf L}^{1} \\ \hat{\bf r}^{2} & \hat{\bf p}^{2} & \hat{\bf L}^{2} \end{bmatrix} \begin{bmatrix} \hat{\bf L}_{N}^{\prime 0} \\ \hat{\bf L}_{N}^{\prime 1} \\ \hat{\bf L}_{N}^{\prime 2} \end{bmatrix} \end{equation*}

See LALSimIMRSpinEOBInitialConditionsPrec.c Line 1317.

In [41]:

#LNhat0 = rhat0*LNhatprm0 + phat0*LNhatprm1 + Lhat0*LNhatprm2
#LNhat1 = rhat1*LNhatprm0 + phat1*LNhatprm1 + Lhat1*LNhatprm2
#LNhat2 = rhat2*LNhatprm0 + phat2*LNhatprm1 + Lhat2*LNhatprm2

Tortoise Conversion Matrix [Back to top]¶

$$\label{tortoise_matrix}$$

We're now back to LALSpinPrecHcapRvecDerivative_v3opt.c, Lines 92--96.

From Pan, Buonanno, Buchman, et. al. (2010) Equation (A3) the matrix for the coordinate conversion to tortoise coordinates is

\begin{align*} \begin{pmatrix} 1 + \frac{ x^{2} }{ r^{2} } \left( \xi - 1 \right) & \frac{ x y }{ r^{2} } \left( \xi - 1 \right) & \frac{ x z }{ r^{2} } \left( \xi - 1 \right) \\ \frac{ x y }{ r^{2} } \left( \xi - 1 \right) & 1 + \frac{ y^{2} }{ r^{2} } \left( \xi - 1 \right) & \frac{ y z }{ r^{2} } \left( \xi - 1 \right) \\ \frac{ x z }{ r^{2} } \left( \xi - 1 \right) & \frac{ y z }{ r^{2} } \left( \xi - 1 \right) & 1 + \frac{ z^{2} }{ r^{2} } \left( \xi - 1 \right) \end{pmatrix} \end{align*}

In [42]:

#ximinus1 = xi - 1
#toTort = sp.Array([[1 + x*x*ximinus1/(r*r), x*y*ximinus1/(r*r), x*z*ximinus1/(r*r)],
#                   [x*y*ximinus1/(r*r), 1 + y*y*ximinus1/(r*r), y*z*ximinus1/(r*r)],
#                   [x*z*ximinus1/(r*r), y*z*ximinus1/(r*r), 1 + z*z*ximinus1/(r*r)]])

Output: Output this notebook to $\LaTeX$-formatted PDF file [Back to top]¶

$$\label{latex_pdf_output}$$

The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename Tutorial-SEOBNR_Initial_Conditions.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

In [43]:

cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-SEOBNR_Initial_Conditions")

Created Tutorial-SEOBNR_Initial_Conditions.tex, and compiled LaTeX file to
    PDF file Tutorial-SEOBNR_Initial_Conditions.pdf

Validation: Perform validation checks against LALSuite's SEOBNRv3 code (commit bba40f2) [Back to top]¶

$$\label{validation}$$

In [44]:

## Validation Cell
## Here we perform a validation check by comparing the derivative values to hard-coded values produced by SEOBNRv3
## in LALSuite.  If this check fails, y'all done sump'tin wrong!

#derivative_list = [dHdx,dHdy,dHdz,dHdpx,dHdpy,dHdpz,dHds1x,dHds1y,dHds1z,dHds2x,dHds2y,dHds2z]

#for q in derivative_list:
#    from SEOBNR_Playground_Pycodes.new_q import new_compute_q
#from SEOBNR.constant_coeffs import compute_const_coeffs

#KK, k0, k1, k2, k3, k4, k5, k5l, dSO, dSS = compute_const_coeffs(eta,EMgamma,a)
## The coefficients do agree with LALSuite!

#tortoise = 1 #Only for testing
#Hreal = compute_Hreal(m1, m2, eta, 10.0, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, tortoise, EMgamma)
#Hreal_pert = compute_Hreal(m1, m2, eta, 10.0*(1.+1e-15), 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, tortoise, EMgamma)
#
#termbyterm_dHdx = new_compute_dHdx(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 2, EMgamma)
#termbyterm_dHdy = new_compute_dHdy(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 2, EMgamma)
#termbyterm_dHdz = new_compute_dHdz(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 2, EMgamma)
#termbyterm_dHdpx = new_compute_dHdpx(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 1, EMgamma)
#termbyterm_dHdpy = new_compute_dHdpy(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 1, EMgamma)
#termbyterm_dHdpz = new_compute_dHdpz(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 1, EMgamma)
#termbyterm_dHds1x = new_compute_dHds1x(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 1, EMgamma)
#termbyterm_dHds1y = new_compute_dHds1y(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 1, EMgamma)
#termbyterm_dHds1z = new_compute_dHds1z(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 1, EMgamma)
#termbyterm_dHds2x = new_compute_dHds2x(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 1, EMgamma)
#termbyterm_dHds2y = new_compute_dHds2y(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 1, EMgamma)
#termbyterm_dHds2z = new_compute_dHds2z(m1, m2, eta, 10, 11.0, 12.0,
#                      0.01, 0.02, 0.03,
#                      0.004, 0.005, -0.006,
#                      0.007, -0.008, 0.009,
#                      KK, k0, k1, dSO, dSS, 1, EMgamma)

#print("exact Hreal = %.15e" % Hreal)
#print("pertd Hreal = %.15e" % Hreal_pert)
#print("relative diff in Hreal = %.15e\n" % (np.abs(Hreal - Hreal_pert)/np.abs(Hreal)))

#print("new term-by-term computation of dHdx = %.15e\n" % (termbyterm_dHdx[0]))
#print("new term-by-term computation of dHdy = %.15e\n" % termbyterm_dHdy[0])
#print("new term-by-term computation of dHdz = %.15e\n" % termbyterm_dHdz[0])
#print("new term-by-term computation of dHdpx = %.15e\n" % termbyterm_dHdpx[0])
#print("new term-by-term computation of dHdpy = %.15e\n" % termbyterm_dHdpy[0])
#print("new term-by-term computation of dHdpz = %.15e\n" % termbyterm_dHdpz[0])
#print("new term-by-term computation of dHds1x = %.15e\n" % termbyterm_dHds1x[0])
#print("new term-by-term computation of dHds1y = %.15e\n" % termbyterm_dHds1y[0])
#print("new term-by-term computation of dHds1z = %.15e\n" % termbyterm_dHds1z[0])
#print("new term-by-term computation of dHds2x = %.15e\n" % termbyterm_dHds2x[0])
#print("new term-by-term computation of dHds2y = %.15e\n" % termbyterm_dHds2y[0])
#print("new term-by-term computation of dHds2z = %.15e\n" % termbyterm_dHds2z[0])

In [ ]: