Notebook

NRPyPN Shortcuts¶

Author: Zach Etienne¶

This notebook contains a number of useful shortcuts for constructing post-Newtonian expressions within NRPyPN/SymPy.¶

Notebook Status: Validated through extensive use

Validation Notes: This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented below. Additional validation tests may have been performed, but are as yet, undocumented. (TODO)

Corresponding Python module:¶

NRPyPN_shortcuts.py

Table of Contents¶

$\label{toc}$

Part 1: Declare global variables used as input to PN expressions
Part 2: Vector & arithmetic operations, including dot products (dot(a,b)) and cross products (cross(a,b)) (including the Levi-Civita symbol); also the shortcut div(a,b) to sp.Rational(a,b) for rational number `a/b
Part 3: Numerical evaluation of SymPy expressions (numeval(expr))
Part 4: eval__P_t__and__P_r(): Compute $p_t$ (tangential orbital momentum) and $p_r$ (radial orbital momentum
Part 5: Validate against corresponding Python module
Part 6: LaTeX PDF output: $\LaTeX$ PDF Output

Part 1: Declare global variables used as input to PN expressions [Back to top]¶

$\label{globals}$

Here, after initializing the module, we define

m1 and m2, the masses of the point particles. We will require that m1+m2=1, using $G=c=1$ units, as all results can be rescaled appropriately with total mass.
S1U= $\mathbf{S}_1$ and S2U= $\mathbf{S}_2$ , the dimensionful spin vectors of the point masses
pU= $\mathbf{p} = \mathbf{P}_1/\mu = -\mathbf{P}_2/\mu$ , where $\mathbf{P}_i$ is the momentum vector of mass $i$ with respect to the center of mass, and $\mu=m_1 m_2 / (m_1+m_2)^2$ is the reduced mass
nU= $(\mathbf{X}_1-\mathbf{X}_2)/|\mathbf{X}_1-\mathbf{X}_2|$ , where $\mathbf{X}_i$ is the position vector of mass $i$ with respect to the center of mass
drdt= $\frac{dr}{dt}$ , the radial infall rate of the binary system
Pt= $p_t$ , Pr= $p_r$ , the tangential and radial components of the momentum vector for the binary, with respect to its center of mass.
r=q= $|\mathbf{X}_1-\mathbf{X}_2|$
chi1U= $\mathbf{\chi}_1=\mathbf{S}_1/m_2^2$ and chi2U= $\mathbf{\chi}_2=\mathbf{S}_2/m_2^2$ , the dimensionless spin vectors of the point masses (each with physically permissible values between -1 and +1)

In [1]:

%%writefile NRPyPN_shortcuts-validation.py
# As documented in the NRPyPN notebook
# NRPyPN_shortcuts.ipynb, this Python script
# provides useful shortcuts for inputting
# post-Newtonian expressions into SymPy/NRPy+

# Basic functions:
# dot(a,b): 3-vector dot product
# cross(a,b): 3-vector cross product
# div(a,b): a shortcut for SymPy's sp.Rational(a,b), to declare rational numbers
# num_eval(expr): Numerically evaluates NRPyPN expressions

# Author:  Zach Etienne
#          zachetie **at** gmail **dot* com

# Step 0: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends
import indexedexpNRPyPN as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support

# Step 1: Declare several global variables used
#         throughout NRPyPN
m1,m2 = sp.symbols('m1 m2',real=True)
S1U = ixp.declarerank1("S1U")
S2U = ixp.declarerank1("S2U")
pU = ixp.declarerank1("pU")
nU = ixp.declarerank1("nU")

drdt = sp.symbols('drdt', real=True)
Pt, Pr = sp.symbols('Pt Pr', real=True)
# Some references use r, others use q to represent the
#   distance between the two point masses. This is rather
#   confusing since q is also used to represent the
#   mass ratio m2/m1. However, q is the canonical position
#   variable name in Hamiltonian mechanics, so both are
#   well justified. It should be obvious which is which
#   throughout NRPyPN.
r, q = sp.symbols('r q', real=True)
chi1U = ixp.declarerank1('chi1U')
chi2U = ixp.declarerank1('chi2U')

# Euler-Mascheroni gamma constant:
gamma_EulerMascheroni = sp.EulerGamma

# Derived quantities used in Damour et al papers:
n12U = ixp.zerorank1()
n21U = ixp.zerorank1()
p1U = ixp.zerorank1()
p2U = ixp.zerorank1()
for i in range(3):
    n12U[i] = +nU[i]
    n21U[i] = -nU[i]
    p1U[i]     = +pU[i]
    p2U[i]     = -pU[i]

Overwriting NRPyPN_shortcuts-validation.py

Part 2: Vector operations, including dot products (`dot(a,b)`) and cross products (`cross(a,b)`) (including the Levi-Civita symbol); also the shortcut `div(a,b)` to `sp.Rational(a,b)` for rational number `a/b` [Back to top]¶

$\label{vectorops}$

Here we define

dot(a,b)= $\mathbf{a}\cdot\mathbf{b}$ , for 3-vectors $\mathbf{a}$ and $\mathbf{b}$
cross(a,b)=a×b, for 3-vectors a and b
- ... calling the corresponding ixp.LeviCivitaSymbol_dim3_rank3() function from indexedexp.py, which constructs the Levi-Civita symbol used for cross products
div(a,b)= $a/b$ , for integers $a$ and $b$ . This is a convenient shortcut for SymPy's sp.Rational(a,b) function

In [2]:

%%writefile -a NRPyPN_shortcuts-validation.py

# Step 2.a: Define dot and cross product of vectors
def dot(vec1,vec2):
    vec1_dot_vec2 = sp.sympify(0)
    for i in range(3):
        vec1_dot_vec2 += vec1[i]*vec2[i]
    return vec1_dot_vec2

def cross(vec1,vec2):
    vec1_cross_vec2 = ixp.zerorank1()
    LeviCivitaSymbol = ixp.LeviCivitaSymbol_dim3_rank3()
    for i in range(3):
        for j in range(3):
            for k in range(3):
                vec1_cross_vec2[i] += LeviCivitaSymbol[i][j][k]*vec1[j]*vec2[k]
    return vec1_cross_vec2

# Step 2.b: Construct rational numbers a/b via div(a,b)
def div(a,b):
    return sp.Rational(a,b)

Appending to NRPyPN_shortcuts-validation.py

Part 3: `num_eval(expr)`: numerical evaluation of SymPy expressions [Back to top]¶

$\label{numeval}$

num_eval(expr, qmassratio, nr, nchi1x, nchi1y, nchi1z, nchi2x, nchi2y, nchi2z, nPt=None, ndrdt=None) numerically evaluates expression expr with options to also specify numerical values for the tangential momentum $p_t$ and the radial infall rate $\frac{dr}{dt}$ . Default options evaluate a $q=1$ mass-ratio, nonspinning binary with separation $r=12$ .

In [3]:

%%writefile -a NRPyPN_shortcuts-validation.py

# Step 3: num_eval(expr), a means to numerically evaluate SymPy/NRPyPN
#         expressions
def num_eval(expr,
             qmassratio =  1.0,  # must be >= 1
             nr         = 12.0,  # Orbital separation
             nchi1x     = +0.,
             nchi1y     = +0.,
             nchi1z     = +0.,
             nchi2x     = +0.,
             nchi2y     = +0.,
             nchi2z     = +0.,
             nPt=None, ndrdt=None):

    # DERIVED QUANTITIES BELOW
    # We want m1+m2 = 1, so that
    #         m2/m1 = qmassratio
    # We find below:
    nm1   =          1/(1+qmassratio)
    nm2   = qmassratio/(1+qmassratio)
    # This way nm1+nm2 = (qmassratio+1)/(1+qmassratio) = 1 CHECK
    #      and nm2/nm1 = qmassratio                        CHECK

    nS1U0 = nchi1x*nm1**2
    nS1U1 = nchi1y*nm1**2
    nS1U2 = nchi1z*nm1**2
    nS2U0 = nchi2x*nm2**2
    nS2U1 = nchi2y*nm2**2
    nS2U2 = nchi2z*nm2**2

    if nPt != None:
        expr2 = expr.subs(Pt,nPt)
        expr  = expr2
    if ndrdt != None:
        expr2 = expr.subs(drdt,ndrdt)
        expr  = expr2
    return expr\
.subs(m1,nm1).subs(m2,nm2)\
.subs(S1U[0],nS1U0).subs(S1U[1],nS1U1).subs(S1U[2],nS1U2)\
.subs(S2U[0],nS2U0).subs(S2U[1],nS2U1).subs(S2U[2],nS2U2)\
.subs(chi1U[0],nchi1x).subs(chi1U[1],nchi1y).subs(chi1U[2],nchi1z)\
.subs(chi2U[0],nchi2x).subs(chi2U[1],nchi2y).subs(chi2U[2],nchi2z)\
.subs(r,nr).subs(q,nr).subs(sp.pi,sp.N(sp.pi))\
.subs(gamma_EulerMascheroni,sp.N(sp.EulerGamma))

Appending to NRPyPN_shortcuts-validation.py

Part 4: `eval__P_tandP_r()`: Compute $p_t$ (tangential orbital momentum) and $p_r$ (radial orbital momentum [Back to top]¶

$\label{ptpr}$

In [4]:

%%writefile -a NRPyPN_shortcuts-validation.py

# Shortcut function to just evaluate P_t and P_r
#  -= Inputs =-
#  * qmassratio: mass ratio q>=1
#  * nr: Orbital separation (total coordinate distance between black holes)
#  * nchi1x,nchi1y,nchi1z: dimensionless spin vector for BH 1
#  * nchi2x,nchi2y,nchi2z: dimensionless spin vector for BH 2
#  -= Outputs =-
# The numerical values for
#  * P_t: the tangential momentum
#  * P_r: the radial momentum
def eval__P_t__and__P_r(qmassratio, nr,
                        nchi1x, nchi1y, nchi1z,
                        nchi2x, nchi2y, nchi2z):

    # Compute p_t, the tangential component of momentum
    import PN_p_t as pt
    pt.f_p_t(m1,m2, chi1U,chi2U, r)

    # Compute p_r, the radial component of momentum
    import PN_p_r as pr
    pr.f_p_r(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)

    nPt = num_eval(pt.p_t,
                   qmassratio=qmassratio, nr=nr,
                   nchi1x=nchi1x,nchi1y=nchi1y,nchi1z=nchi1z,
                   nchi2x=nchi2x,nchi2y=nchi2y,nchi2z=nchi2z)

    nPr = num_eval(pr.p_r,
                    qmassratio = qmassratio, nr=nr,
                    nchi1x=nchi1x, nchi1y=nchi1y, nchi1z=nchi1z,
                    nchi2x=nchi2x, nchi2y=nchi2y, nchi2z=nchi2z,  nPt=nPt)
    return nPt, nPr

Appending to NRPyPN_shortcuts-validation.py

Part 5: Validate against corresponding Python module [Back to top]¶

$\label{validate}$

In [5]:

import difflib,sys

print("Printing difference between original NRPyPN_shortcuts.py and this code, NRPyPN_shortcuts-validation.py.")

# Open the files to compare
with open("NRPyPN_shortcuts.py") as file1, open("NRPyPN_shortcuts-validation.py") as file2:
    # Read the lines of each file,
    file1_lines = file1.readlines()
    file2_lines = file2.readlines()
    num_diffs = 0
    for line in difflib.unified_diff(file1_lines, file2_lines, fromfile="NRPyPN_shortcuts.py", tofile="NRPyPN_shortcuts-validation.py"):
        sys.stdout.writelines(line)
        num_diffs = num_diffs + 1
    if num_diffs == 0:
        print("No difference. TEST PASSED!")
    else:
        print("ERROR: Disagreement found with .py file. See differences above.")
        sys.exit(1)

Printing difference between original NRPyPN_shortcuts.py and this code, NRPyPN_shortcuts-validation.py.
No difference. TEST PASSED!

Part 6: 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 NRPyPN_shortcuts.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

In [6]:

import os,sys                    # Standard Python modules for multiplatform OS-level functions
import cmdline_helperNRPyPN as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("NRPyPN_shortcuts")

Created NRPyPN_shortcuts.tex, and compiled LaTeX file to PDF file
    NRPyPN_shortcuts.pdf