{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
"<script>\n",
" window.dataLayer = window.dataLayer || [];\n",
" function gtag(){dataLayer.push(arguments);}\n",
" gtag('js', new Date());\n",
"\n",
" gtag('config', 'UA-59152712-8');\n",
"</script>\n",
"\n",
"# NRPyPN Shortcuts\n",
"\n",
"## Author: Zach Etienne\n",
"\n",
"## This notebook contains a number of useful shortcuts for constructing post-Newtonian expressions within NRPyPN/SymPy.\n",
"\n",
"**Notebook Status:** <font color='orange'><b> Validated through extensive use </b></font>\n",
"\n",
"**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). **Additional validation tests may have been performed, but are as yet, undocumented. (TODO)**\n",
"\n",
"### Corresponding Python module:\n",
"1. [NRPyPN_shortcuts.py](../edit/NRPyPN_shortcuts.py)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='toc'></a>\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"1. Part 1: [Declare global variables used as input to PN expressions](#globals)\n",
"1. Part 2: [Vector & arithmetic operations](#vectorops), including dot products (`dot(a,b)`) and cross products (`cross(a,b)`) (including the [Levi-Civita symbol](https://en.wikipedia.org/wiki/Levi-Civita_symbol)); also the shortcut `div(a,b)` to `sp.Rational(a,b)` for rational number `a/b\n",
"1. Part 3: [Numerical evaluation of SymPy expressions](#numeval) (`numeval(expr)`)\n",
"1. Part 4: [`eval__P_t__and__P_r()`: Compute $p_t$ (tangential orbital momentum) and $p_r$ (radial orbital momentum](#ptpr)\n",
"1. Part 5: [Validate against corresponding Python module](#validate)\n",
"1. Part 6: [LaTeX PDF output](#latex_pdf_output): $\\LaTeX$ PDF Output"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='globals'></a>\n",
"\n",
"# Part 1: Declare global variables used as input to PN expressions \\[Back to [top](#toc)\\]\n",
"$$\\label{globals}$$ \n",
"\n",
"Here, after initializing the module, we define \n",
"* `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.\n",
"* `S1U`=$\\mathbf{S}_1$ and `S2U`=$\\mathbf{S}_2$, the dimensionful spin vectors of the point masses\n",
"* `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\n",
"* `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\n",
"* `drdt`=$\\frac{dr}{dt}$, the radial infall rate of the binary system\n",
"* `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.\n",
"* `r`=`q`=$|\\mathbf{X}_1-\\mathbf{X}_2|$\n",
"* `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)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Overwriting NRPyPN_shortcuts-validation.py\n"
]
}
],
"source": [
"%%writefile NRPyPN_shortcuts-validation.py\n",
"# As documented in the NRPyPN notebook\n",
"# NRPyPN_shortcuts.ipynb, this Python script\n",
"# provides useful shortcuts for inputting\n",
"# post-Newtonian expressions into SymPy/NRPy+\n",
"\n",
"# Basic functions:\n",
"# dot(a,b): 3-vector dot product\n",
"# cross(a,b): 3-vector cross product\n",
"# div(a,b): a shortcut for SymPy's sp.Rational(a,b), to declare rational numbers\n",
"# num_eval(expr): Numerically evaluates NRPyPN expressions\n",
"\n",
"# Author: Zach Etienne\n",
"# zachetie **at** gmail **dot* com\n",
"\n",
"# Step 0: Add NRPy's directory to the path\n",
"# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n",
"import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
"import indexedexpNRPyPN as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
"\n",
"# Step 1: Declare several global variables used\n",
"# throughout NRPyPN\n",
"m1,m2 = sp.symbols('m1 m2',real=True)\n",
"S1U = ixp.declarerank1(\"S1U\")\n",
"S2U = ixp.declarerank1(\"S2U\")\n",
"pU = ixp.declarerank1(\"pU\")\n",
"nU = ixp.declarerank1(\"nU\")\n",
"\n",
"drdt = sp.symbols('drdt', real=True)\n",
"Pt, Pr = sp.symbols('Pt Pr', real=True)\n",
"# Some references use r, others use q to represent the\n",
"# distance between the two point masses. This is rather\n",
"# confusing since q is also used to represent the\n",
"# mass ratio m2/m1. However, q is the canonical position\n",
"# variable name in Hamiltonian mechanics, so both are\n",
"# well justified. It should be obvious which is which\n",
"# throughout NRPyPN.\n",
"r, q = sp.symbols('r q', real=True)\n",
"chi1U = ixp.declarerank1('chi1U')\n",
"chi2U = ixp.declarerank1('chi2U')\n",
"\n",
"# Euler-Mascheroni gamma constant:\n",
"gamma_EulerMascheroni = sp.EulerGamma\n",
"\n",
"# Derived quantities used in Damour et al papers:\n",
"n12U = ixp.zerorank1()\n",
"n21U = ixp.zerorank1()\n",
"p1U = ixp.zerorank1()\n",
"p2U = ixp.zerorank1()\n",
"for i in range(3):\n",
" n12U[i] = +nU[i]\n",
" n21U[i] = -nU[i]\n",
" p1U[i] = +pU[i]\n",
" p2U[i] = -pU[i]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='vectorops'></a>\n",
"\n",
"# Part 2: Vector operations, including dot products (`dot(a,b)`) and cross products (`cross(a,b)`) (including the [Levi-Civita symbol](https://en.wikipedia.org/wiki/Levi-Civita_symbol)); also the shortcut `div(a,b)` to `sp.Rational(a,b)` for rational number `a/b` \\[Back to [top](#toc)\\]\n",
"$$\\label{vectorops}$$ \n",
"\n",
"Here we define \n",
"* `dot(a,b)`=$\\mathbf{a}\\cdot\\mathbf{b}$, for 3-vectors $\\mathbf{a}$ and $\\mathbf{b}$\n",
"* `cross(a,b)`=$\\mathbf{a}\\times\\mathbf{b}$, for 3-vectors $\\mathbf{a}$ and $\\mathbf{b}$\n",
" * ... calling the corresponding `ixp.LeviCivitaSymbol_dim3_rank3()` function from [indexedexp.py](../../edit/indexedexp.py), which constructs the [Levi-Civita symbol](https://en.wikipedia.org/wiki/Levi-Civita_symbol) used for cross products\n",
"* `div(a,b)`=$a/b$, for integers $a$ and $b$. This is a convenient shortcut for SymPy's `sp.Rational(a,b)` function"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to NRPyPN_shortcuts-validation.py\n"
]
}
],
"source": [
"%%writefile -a NRPyPN_shortcuts-validation.py\n",
"\n",
"# Step 2.a: Define dot and cross product of vectors\n",
"def dot(vec1,vec2):\n",
" vec1_dot_vec2 = sp.sympify(0)\n",
" for i in range(3):\n",
" vec1_dot_vec2 += vec1[i]*vec2[i]\n",
" return vec1_dot_vec2\n",
"\n",
"def cross(vec1,vec2):\n",
" vec1_cross_vec2 = ixp.zerorank1()\n",
" LeviCivitaSymbol = ixp.LeviCivitaSymbol_dim3_rank3()\n",
" for i in range(3):\n",
" for j in range(3):\n",
" for k in range(3):\n",
" vec1_cross_vec2[i] += LeviCivitaSymbol[i][j][k]*vec1[j]*vec2[k]\n",
" return vec1_cross_vec2\n",
"\n",
"# Step 2.b: Construct rational numbers a/b via div(a,b)\n",
"def div(a,b):\n",
" return sp.Rational(a,b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='numeval'></a>\n",
"\n",
"# Part 3: `num_eval(expr)`: numerical evaluation of SymPy expressions \\[Back to [top](#toc)\\]\n",
"$$\\label{numeval}$$ \n",
"\n",
"* `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$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to NRPyPN_shortcuts-validation.py\n"
]
}
],
"source": [
"%%writefile -a NRPyPN_shortcuts-validation.py\n",
"\n",
"# Step 3: num_eval(expr), a means to numerically evaluate SymPy/NRPyPN\n",
"# expressions\n",
"def num_eval(expr,\n",
" qmassratio = 1.0, # must be >= 1\n",
" nr = 12.0, # Orbital separation\n",
" nchi1x = +0.,\n",
" nchi1y = +0.,\n",
" nchi1z = +0.,\n",
" nchi2x = +0.,\n",
" nchi2y = +0.,\n",
" nchi2z = +0.,\n",
" nPt=None, ndrdt=None):\n",
"\n",
" # DERIVED QUANTITIES BELOW\n",
" # We want m1+m2 = 1, so that\n",
" # m2/m1 = qmassratio\n",
" # We find below:\n",
" nm1 = 1/(1+qmassratio)\n",
" nm2 = qmassratio/(1+qmassratio)\n",
" # This way nm1+nm2 = (qmassratio+1)/(1+qmassratio) = 1 CHECK\n",
" # and nm2/nm1 = qmassratio CHECK\n",
"\n",
" nS1U0 = nchi1x*nm1**2\n",
" nS1U1 = nchi1y*nm1**2\n",
" nS1U2 = nchi1z*nm1**2\n",
" nS2U0 = nchi2x*nm2**2\n",
" nS2U1 = nchi2y*nm2**2\n",
" nS2U2 = nchi2z*nm2**2\n",
"\n",
" if nPt != None:\n",
" expr2 = expr.subs(Pt,nPt)\n",
" expr = expr2\n",
" if ndrdt != None:\n",
" expr2 = expr.subs(drdt,ndrdt)\n",
" expr = expr2\n",
" return expr\\\n",
".subs(m1,nm1).subs(m2,nm2)\\\n",
".subs(S1U[0],nS1U0).subs(S1U[1],nS1U1).subs(S1U[2],nS1U2)\\\n",
".subs(S2U[0],nS2U0).subs(S2U[1],nS2U1).subs(S2U[2],nS2U2)\\\n",
".subs(chi1U[0],nchi1x).subs(chi1U[1],nchi1y).subs(chi1U[2],nchi1z)\\\n",
".subs(chi2U[0],nchi2x).subs(chi2U[1],nchi2y).subs(chi2U[2],nchi2z)\\\n",
".subs(r,nr).subs(q,nr).subs(sp.pi,sp.N(sp.pi))\\\n",
".subs(gamma_EulerMascheroni,sp.N(sp.EulerGamma))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='ptpr'></a>\n",
"\n",
"# Part 4: `eval__P_t__and__P_r()`: Compute $p_t$ (tangential orbital momentum) and $p_r$ (radial orbital momentum \\[Back to [top](#toc)\\]\n",
"$$\\label{ptpr}$$ "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to NRPyPN_shortcuts-validation.py\n"
]
}
],
"source": [
"%%writefile -a NRPyPN_shortcuts-validation.py\n",
"\n",
"# Shortcut function to just evaluate P_t and P_r\n",
"# -= Inputs =-\n",
"# * qmassratio: mass ratio q>=1\n",
"# * nr: Orbital separation (total coordinate distance between black holes)\n",
"# * nchi1x,nchi1y,nchi1z: dimensionless spin vector for BH 1\n",
"# * nchi2x,nchi2y,nchi2z: dimensionless spin vector for BH 2\n",
"# -= Outputs =-\n",
"# The numerical values for\n",
"# * P_t: the tangential momentum\n",
"# * P_r: the radial momentum\n",
"def eval__P_t__and__P_r(qmassratio, nr,\n",
" nchi1x, nchi1y, nchi1z,\n",
" nchi2x, nchi2y, nchi2z):\n",
"\n",
" # Compute p_t, the tangential component of momentum\n",
" import PN_p_t as pt\n",
" pt.f_p_t(m1,m2, chi1U,chi2U, r)\n",
"\n",
" # Compute p_r, the radial component of momentum\n",
" import PN_p_r as pr\n",
" pr.f_p_r(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"\n",
" nPt = num_eval(pt.p_t,\n",
" qmassratio=qmassratio, nr=nr,\n",
" nchi1x=nchi1x,nchi1y=nchi1y,nchi1z=nchi1z,\n",
" nchi2x=nchi2x,nchi2y=nchi2y,nchi2z=nchi2z)\n",
"\n",
" nPr = num_eval(pr.p_r,\n",
" qmassratio = qmassratio, nr=nr,\n",
" nchi1x=nchi1x, nchi1y=nchi1y, nchi1z=nchi1z,\n",
" nchi2x=nchi2x, nchi2y=nchi2y, nchi2z=nchi2z, nPt=nPt)\n",
" return nPt, nPr\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='validate'></a>\n",
"\n",
"# Part 5: Validate against corresponding Python module \\[Back to [top](#toc)\\]\n",
"$$\\label{validate}$$ "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Printing difference between original NRPyPN_shortcuts.py and this code, NRPyPN_shortcuts-validation.py.\n",
"No difference. TEST PASSED!\n"
]
}
],
"source": [
"import difflib,sys\n",
"\n",
"print(\"Printing difference between original NRPyPN_shortcuts.py and this code, NRPyPN_shortcuts-validation.py.\")\n",
"\n",
"# Open the files to compare\n",
"with open(\"NRPyPN_shortcuts.py\") as file1, open(\"NRPyPN_shortcuts-validation.py\") as file2:\n",
" # Read the lines of each file,\n",
" file1_lines = file1.readlines()\n",
" file2_lines = file2.readlines()\n",
" num_diffs = 0\n",
" for line in difflib.unified_diff(file1_lines, file2_lines, fromfile=\"NRPyPN_shortcuts.py\", tofile=\"NRPyPN_shortcuts-validation.py\"):\n",
" sys.stdout.writelines(line)\n",
" num_diffs = num_diffs + 1\n",
" if num_diffs == 0:\n",
" print(\"No difference. TEST PASSED!\")\n",
" else:\n",
" print(\"ERROR: Disagreement found with .py file. See differences above.\")\n",
" sys.exit(1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='latex_pdf_output'></a>\n",
"\n",
"# Part 6: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
"$$\\label{latex_pdf_output}$$\n",
"\n",
"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\n",
"[NRPyPN_shortcuts.pdf](NRPyPN_shortcuts.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created NRPyPN_shortcuts.tex, and compiled LaTeX file to PDF file\n",
" NRPyPN_shortcuts.pdf\n"
]
}
],
"source": [
"import os,sys # Standard Python modules for multiplatform OS-level functions\n",
"import cmdline_helperNRPyPN as cmd # NRPy+: Multi-platform Python command-line interface\n",
"cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"NRPyPN_shortcuts\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.11"
}
},
"nbformat": 4,
"nbformat_minor": 4
}