{
"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",
"# $p_r$, the radial component of the momentum vector, up to and including 3.5 PN order\n",
"\n",
"## Author: Zach Etienne\n",
"\n",
"## This notebook constructs the radial component of the momentum vector, up to and including 3.5 post-Newtonian order, and performs numerous validation tests. \n",
"\n",
"**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
"\n",
"**Validation Notes:** All PN expressions in this notebook were transcribed twice by hand on separate occasions, and expressions were corrected as needed to ensure consistency with published work. Published work was cross-validated and typo(s) in published work were corrected. In addition, all expressions in this notebook were validated against those in the Mathematica notebook used by [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036) (thanks to Toni Ramos-Buades for sharing this!) Finally, 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.**ß\n",
"\n",
"### This notebook exists as the following Python module:\n",
"1. [PN_p_r.py](../../edit/NRPyPN/PN_p_r.py)\n",
"\n",
"### This notebook & corresponding Python module depend on the following NRPy+/NRPyPN Python modules:\n",
"1. [indexedexp.py](../../edit/indexedexp.py): [**documentation+tutorial**](../Tutorial-Indexed_Expressions.ipynb)\n",
"1. [NRPyPN_shortcuts.py](../../edit/NRPyPN/NRPyPN_shortcuts.py): [**documentation**](NRPyPN_shortcuts.ipynb)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='toc'></a>\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"1. [Part 0](#imports): Import needed Python modules\n",
"1. [Part 1](#strategy): Basic strategy for computing $p_r$\n",
" 1. [Part 1.a](#fullhamiltonian): Construct the full Hamiltonian consistent with a binary orbiting instantaneously in the $x$-$y$ plane\n",
" 1. [Part 1.b](#dr_dt) Computing $\\frac{dr}{dt}$\n",
" 1. [Part 1.c](#p_r) Computing $p_r\\left(\\frac{dr}{dt},...\\right)$ using two approaches\n",
" 1. [Part 1.c.i](#approach1_fullham) Approach 1: Use the full Hamiltonian\n",
" 1. [Part 1.c.ii](#approach2_fullham) Approach 2 (default approach): use the full Hamiltonian, with higher-order terms over the course of derivation removed\n",
"1. [Part 2](#code_validation): Validation Tests\n",
" 1. [Part 2.a](#code_validation_trans): Validation of transcribed versions for $p_r$\n",
" 1. [Part 2.b](#code_validation_2approaches): Comparison of two approaches for computing $p_r$\n",
" 1. [Part 2.c](#code_validationpub): Validation against trusted numerical values (i.e., in Table V of [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036))\n",
" 1. [Part 2.d](#code_validationpython): Validation against corresponding Python module\n",
"1. Part 2: [LaTeX PDF output](#latex_pdf_output): $\\LaTeX$ PDF Output"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='imports'></a>\n",
"\n",
"# Part 0: Import needed Python modules \\[Back to [top](#toc)\\]\n",
"$$\\label{imports}$$ "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:03.837101Z",
"iopub.status.busy": "2021-03-07T17:36:03.836272Z",
"iopub.status.idle": "2021-03-07T17:36:04.167721Z",
"shell.execute_reply": "2021-03-07T17:36:04.167066Z"
}
},
"outputs": [],
"source": [
"# 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",
"from NRPyPN_shortcuts import Pt,Pr,nU,div # NRPyPN: shortcuts for e.g., vector operations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='strategy'></a>\n",
"\n",
"# Part 1: Basic strategy for computing $p_r$ \\[Back to [top](#toc)\\]\n",
"$$\\label{strategy}$$ \n",
"\n",
"\n",
"<a id='fullhamiltonian'></a>\n",
"\n",
"## Part 1.a: Construct the full Hamiltonian consistent with a binary orbiting instantaneously in the $x$-$y$ plane \\[Back to [top](#toc)\\]\n",
"$$\\label{fullhamiltonian}$$ \n",
"\n",
"As detailed in [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036), the basic strategy for computing $p_r$ first involves constructing the full Hamiltonian expression assuming that, consistent with a binary system initially on the $x$-axis and orbiting instantaneously on the $x$-$y$ plane, the momenta and normal vectors are given by:\n",
"\n",
"\\begin{align}\n",
"P_1^x &= -p_r\\\\\n",
"P_2^x &= +p_r\\\\\n",
"P_1^y &= +p_t\\\\\n",
"P_2^y &= -p_t\\\\\n",
"P_1^z = P_2^z &= 0\\\\\n",
"\\mathbf{n} &= (1,0,0)\n",
"\\end{align}\n",
"\n",
"Now let's construct the full Hamiltonian, applying these assumptions as we go:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:04.180222Z",
"iopub.status.busy": "2021-03-07T17:36:04.179611Z",
"iopub.status.idle": "2021-03-07T17:36:04.181733Z",
"shell.execute_reply": "2021-03-07T17:36:04.182174Z"
}
},
"outputs": [],
"source": [
"# Step 1: Construct full Hamiltonian\n",
"# expression for a binary instantaneously\n",
"# orbiting on the xy plane, store\n",
"# result to Htot_xyplane_binary\n",
"def f_Htot_xyplane_binary(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r):\n",
" def make_replacements(expr):\n",
" zero = sp.sympify(0)\n",
" one = sp.sympify(1)\n",
" return expr.subs(p1U[1], Pt).subs(p2U[1], -Pt).subs(p1U[2], zero).subs(p2U[2], zero).subs(p1U[0], -Pr).subs(\n",
" p2U[0], Pr) \\\n",
" .subs(nU[0], one).subs(nU[1], zero).subs(nU[2], zero)\n",
"\n",
" import PN_Hamiltonian_NS as H_NS\n",
" H_NS.f_H_Newt__H_NS_1PN__H_NS_2PN(m1, m2, p1U, n12U, r)\n",
" H_NS.f_H_NS_3PN(m1, m2, p1U, n12U, r)\n",
"\n",
" global Htot_xyplane_binary\n",
" Htot_xyplane_binary = make_replacements(+H_NS.H_Newt\n",
" + H_NS.H_NS_1PN\n",
" + H_NS.H_NS_2PN\n",
" + H_NS.H_NS_3PN)\n",
"\n",
" import PN_Hamiltonian_SO as H_SO\n",
" H_SO.f_H_SO_1p5PN(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)\n",
" H_SO.f_H_SO_2p5PN(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)\n",
" H_SO.f_H_SO_3p5PN(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)\n",
" Htot_xyplane_binary += make_replacements(+H_SO.H_SO_1p5PN\n",
" + H_SO.H_SO_2p5PN\n",
" + H_SO.H_SO_3p5PN)\n",
"\n",
" import PN_Hamiltonian_SS as H_SS\n",
" H_SS.f_H_SS_2PN(m1, m2, S1U, S2U, nU, r)\n",
" H_SS.f_H_SS_S1S2_3PN(m1, m2, n12U, S1U, S2U, p1U, p2U, r)\n",
" H_SS.f_H_SS_S1sq_S2sq_3PN(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)\n",
" Htot_xyplane_binary += make_replacements(+H_SS.H_SS_2PN\n",
" + H_SS.H_SS_S1S2_3PN\n",
" + H_SS.H_SS_S1sq_S2sq_3PN)\n",
"\n",
" import PN_Hamiltonian_SSS as H_SSS\n",
" H_SSS.f_H_SSS_3PN(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)\n",
" Htot_xyplane_binary += make_replacements(+H_SSS.H_SSS_3PN)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='dr_dt'></a>\n",
"\n",
"## Part 1.b: Computing $\\frac{dr}{dt}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{dr_dt}$$ \n",
"\n",
"Hamilton's equations of motion imply that\n",
"\n",
"$$\n",
"\\frac{dr}{dt} = \\frac{\\partial H}{\\partial p_r}.\n",
"$$\n",
"\n",
"Next we Taylor expand $\\frac{\\partial H}{\\partial p_r}$ in powers of $p_r$, about $p_r=0$:\n",
"\n",
"\\begin{align}\n",
"\\frac{dr}{dt} = \\frac{\\partial H}{\\partial p_r} = \\left.\\frac{\\partial H}{\\partial p_r}\\right|_{p_r=0} \n",
"+ \\frac{1}{1!} p_r \\left.\\frac{\\partial^2 H}{\\partial p_r^2}\\right|_{p_r=0} \n",
"+ \\frac{1}{2!} p_r^2 \\left.\\frac{\\partial^3 H}{\\partial p_r^3}\\right|_{p_r=0} \n",
"+ \\frac{1}{3!} p_r^3 \\left.\\frac{\\partial^4 H}{\\partial p_r^4}\\right|_{p_r=0} \n",
"+ \\mathcal{O}(p_r^4)\n",
"\\end{align}\n",
"so to first order we get\n",
"\n",
"$$\n",
"p_r \\approx \\left(\\frac{dr}{dt} - \\left.\\frac{\\partial H}{\\partial p_r}\\right|_{p_r=0} \\right) \\left( \\left.\\frac{\\partial^2 H}{\\partial p_r^2}\\right|_{p_r=0} \\right)^{-1}.\n",
"$$\n",
"\n",
"Given the input masses and spins, we can compute $p_t$ using the formula given in [this NRPyPN notebook](PN-p_t.ipynb) (from [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036)), and the above will lead us with one equation and two unknowns: $p_r$ and $\\frac{dr}{dt}$. This is an equation we can derive directly from our Hamiltonian expression, and compare the output to the expression derived to 3.5PN order in [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036).\n",
"\n",
"Let's next construct an expression for $\\frac{dr}{dt}$ in terms of known quantities, via\n",
"\n",
"$$\n",
"\\frac{dr}{dt}=\\left(\\frac{dE_{\\rm GW}}{dt}+\\frac{dM}{dt}\\right)\\left[\\frac{dH_{\\rm circ}}{dr}\\right]^{-1},\n",
"$$\n",
"where\n",
"\n",
"$$\n",
"\\frac{dH_{\\rm circ}(r,p_t(r))}{dr} = \\frac{\\partial H(p_r=0)}{\\partial r} \n",
"+ \\frac{\\partial H(p_r=0)}{\\partial p_t} \\frac{\\partial p_t}{\\partial r},\n",
"$$\n",
"\n",
"and the total energy flux\n",
"$$\n",
"\\mathcal{E}(M\\Omega,...) = \\left(\\frac{dE_{\\rm GW}}{dt}+\\frac{dM}{dt}\\right) \n",
"$$ \n",
"\n",
"is given in terms of input parameters (e.g., masses and spins), plus $M\\Omega$, which we also constructed in terms of the input masses, spins, and orbital separation $r$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:04.189662Z",
"iopub.status.busy": "2021-03-07T17:36:04.188981Z",
"iopub.status.idle": "2021-03-07T17:36:04.191200Z",
"shell.execute_reply": "2021-03-07T17:36:04.191660Z"
}
},
"outputs": [],
"source": [
"# Function for computing dr/dt\n",
"def f_dr_dt(Htot_xyplane_binary, m1,m2, n12U, chi1U,chi2U, S1U,S2U, r):\n",
" # First compute p_t\n",
" import PN_p_t as pt\n",
" pt.f_p_t(m1,m2, chi1U,chi2U, r)\n",
"\n",
" # Then compute dH_{circ}/dr = partial_H(p_r=0)/partial_r\n",
" # + partial_H(p_r=0)/partial_{p_t} partial_{p_t}/partial_r\n",
" dHcirc_dr = (+sp.diff(Htot_xyplane_binary.subs(Pr,sp.sympify(0)),r)\n",
" +sp.diff(Htot_xyplane_binary.subs(Pr,sp.sympify(0)),Pt)*sp.diff(pt.p_t,r))\n",
"\n",
" # Then compute M\\Omega\n",
" import PN_MOmega as MOm\n",
" MOm.f_MOmega(m1,m2, chi1U,chi2U, r)\n",
"\n",
" # Next compute dE_GW_dt_plus_dM_dt\n",
" import PN_dE_GW_dt_and_dM_dt as dEdt\n",
" dEdt.f_dE_GW_dt_and_dM_dt(MOm.MOmega, m1,m2, n12U, S1U,S2U)\n",
"\n",
" # Finally, compute dr/dt\n",
" global dr_dt\n",
" dr_dt = dEdt.dE_GW_dt_plus_dM_dt / dHcirc_dr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='p_r'></a>\n",
"\n",
"## Part 1.c: Computing $p_r\\left(\\frac{dr}{dt},...\\right)$ using two approaches \\[Back to [top](#toc)\\]\n",
"$$\\label{p_r}$$ \n",
"\n",
"\n",
"<a id='approach1_fullham'></a>\n",
"\n",
"### Part 1.c.i: Approach 1: Use the full Hamiltonian \\[Back to [top](#toc)\\]\n",
"$$\\label{approach1_fullham}$$ \n",
"\n",
"First, we use the approximation based on the Hamiltonian computed above:\n",
"\n",
"$$\n",
"p_r \\approx \\left(\\frac{dr}{dt} - \\left.\\frac{\\partial H}{\\partial p_r}\\right|_{p_r=0} \\right) \\left[ \\left.\\frac{\\partial^2 H}{\\partial p_r^2}\\right|_{p_r=0} \\right]^{-1}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:04.198517Z",
"iopub.status.busy": "2021-03-07T17:36:04.197850Z",
"iopub.status.idle": "2021-03-07T17:36:04.200066Z",
"shell.execute_reply": "2021-03-07T17:36:04.200524Z"
}
},
"outputs": [],
"source": [
"# Next we compute p_r as a function of dr_dt (unknown) and known quantities using\n",
"# p_r \\approx [dr/dt - (partial_H/partial_{p_r})|_{p_r=0}] * [(partial^2_{H}/partial_{p_r^2})|_{p_r=0}]^{-1}\n",
"def f_p_r_fullHam(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r):\n",
" f_Htot_xyplane_binary(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)\n",
" f_dr_dt(Htot_xyplane_binary, m1,m2, n12U, chi1U,chi2U, S1U,S2U, r)\n",
"\n",
" dHdpr_przero = sp.diff(Htot_xyplane_binary,Pr).subs(Pr,sp.sympify(0))\n",
" d2Hdpr2_przero = sp.diff(sp.diff(Htot_xyplane_binary,Pr),Pr).subs(Pr,sp.sympify(0))\n",
" global p_r_fullHam\n",
" p_r_fullHam = (dr_dt - dHdpr_przero)/(d2Hdpr2_przero)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='approach2_default'></a>\n",
"\n",
"### Part 1.c.ii: Approach 2 (default approach): use the full Hamiltonian, with higher-order terms over the course of derivation removed \\[Back to [top](#toc)\\]\n",
"$$\\label{approach2_default}$$ \n",
"\n",
"This approach (Eq 2.16 of [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036)) uses the same approximation as Approach 1\n",
"$$\n",
"p_r \\approx \\left(\\frac{dr}{dt} - \\left.\\frac{\\partial H}{\\partial p_r}\\right|_{p_r=0} \\right) \\left[ \\left.\\frac{\\partial^2 H}{\\partial p_r^2}\\right|_{p_r=0} \\right]^{-1},\n",
"$$\n",
"except it throws away higher-order terms.\n",
"\n",
"To reduce the possibility of copying error, the equation for $p_r$ is taken directly from the arXiv LaTeX source code of Eq 2.16 in [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036), and only mildly formatted to (1) improve presentation in Jupyter notebooks and (2) to ensure some degree of consistency in notation across different terms in other NRPyPN notebooks:\n",
"\n",
"\\begin{equation}\n",
"\\begin{split}\n",
"p_r &= \\left[-\\frac{dr}{dt} \\right.\\\\\n",
"&\\quad\\quad \\left. +\\frac{1}{r^{7/2}} \\left( -\\frac{(6 q+13) q^2 \\chi _{1x} \\chi _{2 y}}{4 (q+1)^4}-\\frac{(6 q+1) q^2 \\chi _{2 x} \\chi _{2 y}}{4 (q+1)^4}+\\chi _{1y} \\left(-\\frac{q (q+6) \\chi _{1x}}{4 (q+1)^4}-\\frac{q (13 q+6) \\chi _{2 x}}{4 (q+1)^4}\\right)\\right) \\right.\\\\\n",
"&\\quad\\quad \\left. +\\frac{1}{r^4} \\left( \\chi _{1z} \\left(\\frac{3 q (5 q+2) \\chi _{1x} \\chi _{2 y}}{2 (q+1)^4} -\\frac{3 q^2 (2 q+5) \\chi _{2 x} \\chi _{2 y}}{2 (q+1)^4}\\right)+\\chi _{1y} \\chi _{2 z} \\left(\\frac{3 q^2 (2 q+5) \\chi _{2 x}}{2 (q+1)^4}-\\frac{3 q (5 q+2) \\chi _{1x}}{2 (q+1)^4}\\right)\\right)\\right] \\\\\n",
"& \\times \\left[ -\\frac{(q+1)^2}{q}-\\frac{1 \\left(-7 q^2-15 q-7\\right)}{2 q r} \\right.\n",
" \\\\\n",
"&\\quad\\quad \\left. -\\frac{47 q^4+229 q^3+363 q^2+229 q+47}{8 q (q+1)^2 r^2} -\\frac{1}{r^{5/2}}\\left( \\frac{\\left(4 q^2+11 q+12\\right) \\chi _{1z}}{4 q (q+1)}+\\frac{\\left(12 q^2+11 q+4\\right) \\chi _{2 z}}{4 (q+1)} \\right) \\right. \\\\\n",
"&\\quad\\quad \\left. \\left.- \\frac{1}{r^{7/2}} \\left( \\frac{\\left(-53 q^5-357 q^4-1097 q^3-1486 q^2-842 q-144\\right) \\chi _{1z}}{16 q (q+1)^4}+\\frac{\\left(-144 q^5-842 q^4-1486 q^3-1097 q^2-357 q-53\\right) \\chi _{2 z}}{16 (q+1)^4} \\right)\\right. \\right. \\\\\n",
"&\\quad\\quad \\left. -\\frac{1}{r^3} \\left(\\frac{\\left(q^2+9 q+9\\right) \\chi _{1x}^2}{2 q (q+1)^2}+\\frac{\\left(3 q^2+5 q+3\\right) \\chi _{2 x} \\chi _{1x}}{(q+1)^2}+\\frac{\\left(3 q^2+8 q+3\\right) \\chi _{1y} \\chi _{2 y}}{2 (q+1)^2}-\\frac{9 q^2 \\chi _{2 y}^2}{4 (q+1)}+\\frac{\\left(3 q^2+8 q+3\\right) \\chi _{1z} \\chi _{2 z}}{2 (q+1)^2}-\\frac{9 q^2 \\chi _{2 z}^2}{4 (q+1)} \\right. \\right.\n",
"\\\\\n",
"&\\quad\\quad\\quad \\left. \\left. +\\frac{\\left(9 q^3+9 q^2+q\\right) \\chi _{2 x}^2}{2 (q+1)^2}+\\frac{-363 q^6-2608 q^5-7324 q^4-10161 q^3-7324 q^2-2608 q-363}{48 q (q+1)^4}-\\frac{9 \\chi _{1y}^2}{4 q (q+1)}-\\frac{9 \\chi _{1z}^2}{4 q (q+1)}-\\frac{\\pi ^2}{16} \\right) \\right]^{-1}. \n",
"\\end{split}\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:04.221654Z",
"iopub.status.busy": "2021-03-07T17:36:04.220653Z",
"iopub.status.idle": "2021-03-07T17:36:04.222952Z",
"shell.execute_reply": "2021-03-07T17:36:04.223443Z"
}
},
"outputs": [],
"source": [
"# Here's the Ramos-Buades, Husa, and Pratten (2018)\n",
"# approach for computing p_r.\n",
"# Transcribed from Eq 2.18 of\n",
"# Ramos-Buades, Husa, and Pratten (2018),\n",
"# https://arxiv.org/abs/1810.00036\n",
"def f_p_r(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r):\n",
" q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.\n",
" f_Htot_xyplane_binary(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)\n",
" f_dr_dt(Htot_xyplane_binary, m1,m2, n12U, chi1U,chi2U, S1U,S2U, r)\n",
" chi1x = chi1U[0]\n",
" chi1y = chi1U[1]\n",
" chi1z = chi1U[2]\n",
" chi2x = chi2U[0]\n",
" chi2y = chi2U[1]\n",
" chi2z = chi2U[2]\n",
" p_r_num = (-dr_dt\n",
" +(-(6*q+13)*q**2*chi1x*chi2y/(4*(q+1)**4)\n",
" -(6*q+ 1)*q**2*chi2x*chi2y/(4*(q+1)**4)\n",
" +chi1y*(-q*( q+6)*chi1x/(4*(q+1)**4)\n",
" -q*(13*q+6)*chi2x/(4*(q+1)**4)))/r**div(7,2)\n",
" +(+chi1z*(+3*q *(5*q+2)*chi1x*chi2y/(2*(q+1)**4)\n",
" -3*q**2*(2*q+5)*chi2x*chi2y/(2*(q+1)**4))\n",
" +chi1y*chi2z*(+3*q**2*(2*q+5)*chi2x/(2*(q+1)**4)\n",
" -3*q *(5*q+2)*chi1x/(2*(q+1)**4)))/r**4)\n",
" p_r_den = (-(q+1)**2/q - (-7*q**2-15*q-7)/(2*q*r)\n",
" -(47*q**4 + 229*q**3 + 363*q**2 + 229*q + 47)/(8*q*(q+1)**2*r**2)\n",
" -(+( 4*q**2 + 11*q + 12)*chi1z/(4*q*(q+1))\n",
" +(12*q**2 + 11*q + 4)*chi2z/(4* (q+1)))/r**div(5,2)\n",
" -(+(- 53*q**5 - 357*q**4 - 1097*q**3 - 1486*q**2 - 842*q - 144)*chi1z/(16*q*(q+1)**4)\n",
" +(-144*q**5 - 842*q**4 - 1486*q**3 - 1097*q**2 - 357*q - 53)*chi2z/(16 *(q+1)**4))/r**div(7,2)\n",
" -(+( q**2 + 9*q + 9)*chi1x**2/(2*q*(q+1)**2)\n",
" +(3*q**2 + 5*q + 3)*chi2x*chi1x/((q+1)**2)\n",
" +(3*q**2 + 8*q + 3)*chi1y*chi2y/(2*(q+1)**2)\n",
" -9*q**2*chi2y**2/(4*(q+1))\n",
" +(3*q**2 + 8*q + 3)*chi1z*chi2z/(2*(q+1)**2)\n",
" -9*q**2*chi2z**2/(4*(q+1))\n",
" +(9*q**3 + 9*q**2 + q)*chi2x**2/(2*(q+1)**2)\n",
" +(-363*q**6 - 2608*q**5 - 7324*q**4 - 10161*q**3 - 7324*q**2 - 2608*q - 363)/(48*q*(q+1)**4)\n",
" -9*chi1y**2/(4*q*(q+1))\n",
" -9*chi1z**2/(4*q*(q+1)) - sp.pi**2/16)/r**3)\n",
" global p_r\n",
" p_r = p_r_num/p_r_den"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:04.274711Z",
"iopub.status.busy": "2021-03-07T17:36:04.273890Z",
"iopub.status.idle": "2021-03-07T17:36:04.276276Z",
"shell.execute_reply": "2021-03-07T17:36:04.276759Z"
}
},
"outputs": [],
"source": [
"# Second version, used for validation purposes only.\n",
"def f_p_r_RHP2018v2(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r):\n",
" q = m2/m1\n",
" f_Htot_xyplane_binary(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)\n",
" f_dr_dt(Htot_xyplane_binary, m1,m2, n12U, chi1U,chi2U, S1U,S2U, r)\n",
" chi1x = chi1U[0]\n",
" chi1y = chi1U[1]\n",
" chi1z = chi1U[2]\n",
" chi2x = chi2U[0]\n",
" chi2y = chi2U[1]\n",
" chi2z = chi2U[2]\n",
" p_r_num = (-dr_dt\n",
" +(-(6*q+13)*q**2*chi1x*chi2y/(4*(q+1)**4)\n",
" -(6*q+1)*q**2*chi2x*chi2y/(4*(q+1)**4)\n",
" +chi1y*(-q*(q+6)*chi1x/(4*(q+1)**4) - q*(13*q+6)*chi2x/(4*(q+1)**4)))/r**div(7,2)\n",
" +(+chi1z*(+div(3,2)*(q* (5*q+2)*chi1x*chi2y)/(q+1)**4\n",
" -div(3,2)*(q**2*(2*q+5)*chi2x*chi2y)/(q+1)**4)\n",
" +chi1y*chi2z*(+div(3,2)*(q**2*(2*q+5)*chi2x)/(q+1)**4\n",
" -div(3,2)*(q* (5*q+2)*chi1x)/(q+1)**4))/r**4)\n",
" p_r_den = (-(q+1)**2/q\n",
" -(-7*q**2-15*q-7)/(2*q*r)\n",
" -(47*q**4 + 229*q**3 + 363*q**2 + 229*q + 47)/(8*q*(q+1)**2*r**2)\n",
" -( (4*q**2+11*q+12)*chi1z/(4*q*(q+1)) + (12*q**2+11*q+4)*chi2z/(4*(q+1)) )/r**div(5,2)\n",
" -(+(-53 *q**5 - 357*q**4 - 1097*q**3 - 1486*q**2 - 842*q - 144)*chi1z/(16*q*(q+1)**4)\n",
" +(-144*q**5 - 842*q**4 - 1486*q**3 - 1097*q**2 - 357*q - 53 )*chi2z/(16* (q+1)**4) )/r**div(7,2)\n",
" -(+( q**2+9*q+9)*chi1x**2/(2*q*(q+1)**2)\n",
" +(3*q**2+5*q+3)*chi2x*chi1x/((q+1)**2)\n",
" +(3*q**2+8*q+3)*chi1y*chi2y/(2*(q+1)**2)\n",
" -(9*q**2*chi2y**2/(4*(q+1)))\n",
" +(3*q**2+8*q+3)*chi1z*chi2z/(2*(q+1)**2)\n",
" -(9*q**2*chi2z**2/(4*(q+1)))\n",
" +(9*q**3+9*q**2+q)*chi2x**2/(2*(q+1)**2)\n",
" +(-363*q**6 - 2608*q**5 - 7324*q**4 - 10161*q**3 - 7324*q**2 - 2608*q - 363)/(48*q*(q+1)**4)\n",
" -(9*chi1y**2)/(4*q*(q+1))\n",
" -(9*chi1z**2)/(4*q*(q+1))\n",
" -sp.pi**2/16) / r**3)\n",
" global p_r_RHP2018v2\n",
" p_r_RHP2018v2 = p_r_num/p_r_den"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:04.300423Z",
"iopub.status.busy": "2021-03-07T17:36:04.299440Z",
"iopub.status.idle": "2021-03-07T17:36:04.302440Z",
"shell.execute_reply": "2021-03-07T17:36:04.301860Z"
}
},
"outputs": [],
"source": [
"# Third version, directly from Toni Ramos-Buades' Mathematica notebook (thanks Toni!)\n",
"def f_p_r_RHP2018v3(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r):\n",
" q = m2/m1\n",
" f_Htot_xyplane_binary(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)\n",
" f_dr_dt(Htot_xyplane_binary, m1,m2, n12U, chi1U,chi2U, S1U,S2U, r)\n",
" chi1x = chi1U[0]\n",
" chi1y = chi1U[1]\n",
" chi1z = chi1U[2]\n",
" chi2x = chi2U[0]\n",
" chi2y = chi2U[1]\n",
" chi2z = chi2U[2]\n",
" Pi = sp.pi\n",
" p_r_num = (-dr_dt\n",
" +((chi1y*chi2z*((3*chi2x*q**2*(5 + 2*q))/(2*(1 + q)**4) - (3*chi1x*q*(2 + 5*q))/(2*(1 + q)**4)) + chi1z*((-3*chi2x*chi2y*q**2*(5 + 2*q))/(2*(1 + q)**4) + (3*chi1x*chi2y*q*(2 + 5*q))/(2*(1 + q)**4)))/\n",
" r**4 + (-(chi2x*chi2y*q**2*(1 + 6*q))/(4*(1 + q)**4) - (chi1x*chi2y*q**2*(13 + 6*q))/(4*(1 + q)**4) + chi1y*(-(chi1x*q*(6 + q))/(4*(1 + q)**4) - (chi2x*q*(6 + 13*q))/(4*(1 + q)**4)))/r**(sp.Rational(7,2))))\n",
" p_r_den = (-((1 + q)**2/q) - ((chi2z*(-53 - 357*q - 1097*q**2 - 1486*q**3 - 842*q**4 - 144*q**5))/(16*(1 + q)**4) + (chi1z*(-144 - 842*q - 1486*q**2 - 1097*q**3 - 357*q**4 - 53*q**5))/(16*q*(1 + q)**4))/r**(sp.Rational(7,2)) -\n",
" (-Pi**2/16 - (9*chi1y**2)/(4*q*(1 + q)) - (9*chi1z**2)/(4*q*(1 + q)) - (9*chi2y**2*q**2)/(4*(1 + q)) - (9*chi2z**2*q**2)/(4*(1 + q)) + (chi1x**2*(9 + 9*q + q**2))/(2*q*(1 + q)**2) +\n",
" (chi1x*chi2x*(3 + 5*q + 3*q**2))/(1 + q)**2 + (chi1y*chi2y*(3 + 8*q + 3*q**2))/(2*(1 + q)**2) + (chi1z*chi2z*(3 + 8*q + 3*q**2))/(2*(1 + q)**2) + (chi2x**2*(q + 9*q**2 + 9*q**3))/(2*(1 + q)**2) +\n",
" (-363 - 2608*q - 7324*q**2 - 10161*q**3 - 7324*q**4 - 2608*q**5 - 363*q**6)/(48*q*(1 + q)**4))/r**3 - ((chi1z*(12 + 11*q + 4*q**2))/(4*q*(1 + q)) + (chi2z*(4 + 11*q + 12*q**2))/(4*(1 + q)))/r**(sp.Rational(5,2)) -\n",
" (47 + 229*q + 363*q**2 + 229*q**3 + 47*q**4)/(8*q*(1 + q)**2*r**2) - (-7 - 15*q - 7*q**2)/(2*q*r))\n",
" global p_r_RHP2018v3\n",
" p_r_RHP2018v3 = p_r_num/p_r_den"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validation'></a>\n",
"\n",
"# Part 2: Validation Tests \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validation}$$ \n",
"\n",
"<a id='code_validation_trans'></a>\n",
"\n",
"## Part 2.a: Validation of transcribed versions for $p_r$ \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validation_trans}$$ \n",
"\n",
"As a code validation check, we verify agreement between \n",
"* the SymPy expressions transcribed from the cited published work on two separate occasions, and"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:04.311022Z",
"iopub.status.busy": "2021-03-07T17:36:04.310139Z",
"iopub.status.idle": "2021-03-07T17:36:23.374386Z",
"shell.execute_reply": "2021-03-07T17:36:23.374849Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"TRANSCRIPTION TEST PASSES\n"
]
}
],
"source": [
"from NRPyPN_shortcuts import m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r # NRPyPN: import needed input variables\n",
"\n",
"def error(varname):\n",
" print(\"ERROR: When comparing Python module & notebook, \"+varname+\" was found not to match.\")\n",
" sys.exit(1)\n",
"\n",
"# Validation against second transcription of the expressions:\n",
"f_p_r( m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"f_p_r_RHP2018v2(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"f_p_r_RHP2018v3(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"\n",
"if sp.simplify(p_r - p_r_RHP2018v2) != 0: error(\"p_r_RHP2018v2\")\n",
"if sp.simplify(p_r_RHP2018v2 - p_r_RHP2018v3) != 0: error(\"p_r_RHP2018v3\")\n",
"print(\"TRANSCRIPTION TEST PASSES\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validation_2approaches'></a>\n",
"\n",
"## Part 2.b: Comparison of two approaches for computing $p_r$ \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validation_2approaches}$$ \n",
"\n",
"As a validation check, we inject random inputs into the two approaches for computing $p_r$ outlined above and compute the errors."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:23.442368Z",
"iopub.status.busy": "2021-03-07T17:36:23.406517Z",
"iopub.status.idle": "2021-03-07T17:36:39.197683Z",
"shell.execute_reply": "2021-03-07T17:36:39.198186Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0001 | 3.9019e-03 | 1.0678e-04 1.0720e-04 | 6.91 12.27 -0.09 -0.27 0.01 -0.10 0.31 -0.22\n",
"Average relative error: 3.90192e-03\n"
]
}
],
"source": [
"import random # Standard Python module: provides random number generation functionality\n",
"from NRPyPN_shortcuts import num_eval, gamma_EulerMascheroni # NRPyPN: import numerical evaluation routine & gamma constant\n",
"\n",
"def eval_random(i,trusted, other, p_t):\n",
" random.seed(i)\n",
"\n",
" qmassratio = 1.0 + 7*random.random() # must be >= 1\n",
" nr = 10. + 3*random.random() # Orbital separation\n",
" # Choose spins so that the total spin magnitude does not exceed 1.\n",
" nchi1x = -0.55 + 1.1*random.random()\n",
" nchi1y = -0.55 + 1.1*random.random()\n",
" nchi1z = -0.55 + 1.1*random.random()\n",
" nchi2x = -0.55 + 1.1*random.random()\n",
" nchi2y = -0.55 + 1.1*random.random()\n",
" nchi2z = -0.55 + 1.1*random.random()\n",
"\n",
" nPt = num_eval(p_t,\n",
" qmassratio=qmassratio, nr=nr,\n",
" nchi1x=nchi1x,nchi1y=nchi1y,nchi1z=nchi1z,\n",
" nchi2x=nchi2x,nchi2y=nchi2y,nchi2z=nchi2z)\n",
" trusted_result = num_eval(trusted,qmassratio=qmassratio, nr=nr,\n",
" nchi1x=nchi1x,nchi1y=nchi1y,nchi1z=nchi1z,\n",
" nchi2x=nchi2x,nchi2y=nchi2y,nchi2z=nchi2z,\n",
" nPt=nPt)\n",
" other_result= num_eval(other,qmassratio=qmassratio, nr=nr,\n",
" nchi1x=nchi1x,nchi1y=nchi1y,nchi1z=nchi1z,\n",
" nchi2x=nchi2x,nchi2y=nchi2y,nchi2z=nchi2z,\n",
" nPt=nPt)\n",
"\n",
" relerror = abs((trusted_result-other_result)/trusted_result)\n",
" print(\"%04d\" % (i+1), \"|\", #\"%.2f\" % int(count+1)/int(i+1),\n",
" \"%.4e\" % (relerror),\"|\",\"%.4e\" % trusted_result,\"%.4e\" % other_result,\"|\",\n",
" \"%.2f\" % (qmassratio),\"%.2f\" % nr,\n",
" \"%.2f\" % (nchi1x),\"%.2f\" % (nchi1y),\"%.2f\" % (nchi1z),\n",
" \"%.2f\" % (nchi2x),\"%.2f\" % (nchi2y),\"%.2f\" % (nchi2z))\n",
" return relerror\n",
"\n",
"import PN_p_t as pt\n",
"pt.f_p_t(m1,m2, chi1U,chi2U, r)\n",
"\n",
"f_p_r_fullHam(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"f_p_r( m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"\n",
"num_trials = 1\n",
"relerror_array = []\n",
"if num_trials > 100:\n",
" def master_func(i):\n",
" return eval_random(i,p_r_RHP2018,p_r,pt.p_t)\n",
"\n",
" import multiprocessing\n",
" # ixp.zerorank1(DIM=1000) #0.0\n",
" pool = multiprocessing.Pool()\n",
" relerror_array = pool.map(master_func,range(1000))\n",
" pool.terminate()\n",
" pool.join()\n",
"else:\n",
" for i in range(num_trials):\n",
" relerror_array.append(eval_random(i,p_r,p_r_fullHam,pt.p_t))\n",
"\n",
"avg_relerror = 0\n",
"for i in range(len(relerror_array)):\n",
" avg_relerror += relerror_array[i]\n",
"avg_relerror /= len(relerror_array)\n",
"print(\"Average relative error: \",\"%.5e\" % (avg_relerror))\n",
"\n",
"# Without 3.5PN term in p_r denominator (comparing\n",
"# Eq 2.16 of https://arxiv.org/pdf/1810.00036.pdf with\n",
"# Eq 16 of https://arxiv.org/pdf/1702.00872.pdf):\n",
"# average relative error in p_r over 1000 random inputs = 8.24534e-03\n",
"\n",
"# With 3.5PN term in p_r:\n",
"# average relative error in p_r over 1000 random inputs = 8.19325e-03"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validationpub'></a>\n",
"\n",
"## Part 2.c: Validation against trusted numerical values (i.e., in Table V of [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036)) \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validationpub}$$ \n",
"\n",
"We will measure success in this validation to be within the error bar of the two values of $p_r$ given for each iteration, preferably being closer to the second iteration's value than the first."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:39.272467Z",
"iopub.status.busy": "2021-03-07T17:36:39.236641Z",
"iopub.status.idle": "2021-03-07T17:36:52.177613Z",
"shell.execute_reply": "2021-03-07T17:36:52.176984Z"
}
},
"outputs": [],
"source": [
"import PN_p_t as pt\n",
"pt.f_p_t(m1,m2, chi1U,chi2U, r)\n",
"f_p_r_fullHam(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"f_p_r( m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"import PN_p_r as pr\n",
"pr.f_p_r(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"p_r = pr.p_r"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:52.189467Z",
"iopub.status.busy": "2021-03-07T17:36:52.188506Z",
"iopub.status.idle": "2021-03-07T17:36:52.190879Z",
"shell.execute_reply": "2021-03-07T17:36:52.191396Z"
}
},
"outputs": [],
"source": [
"# Useful function for comparing published & NRPyPN results\n",
"def compare_pub_NPN(desc,pub0,pub1, qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z):\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",
" p_r_approach1 = num_eval(p_r_fullHam, qmassratio = qmassratio, nr=nr,\n",
" nchi1x=nchi1x, nchi1y=nchi1y, nchi1z=nchi1z,\n",
" nchi2x=nchi2x, nchi2y=nchi2y, nchi2z=nchi2z,\n",
" nPt=nPt)\n",
"\n",
" p_r_approach2 = num_eval(p_r, qmassratio = qmassratio, nr=nr,\n",
" nchi1x=nchi1x, nchi1y=nchi1y, nchi1z=nchi1z,\n",
" nchi2x=nchi2x, nchi2y=nchi2y, nchi2z=nchi2z,\n",
" nPt=nPt)\n",
"\n",
" print(\"##################################################\")\n",
" print(\" \"+desc)\n",
" print(\"##################################################\")\n",
" print(\"p_r = %.9f\" % pub0, \" <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)\")\n",
" print(\"p_r = %.9f\" % p_r_approach1,\" <- Result from NRPyPN, approach 1\")\n",
" print(\"p_r = %.9f\" % p_r_approach2,\" <- Result from NRPyPN, approach 2\")\n",
" print(\"p_r = %.9f\" % pub1, \" <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)\")\n",
" relerrorpct1 = abs((pub0-p_r_approach1)/pub0)*100\n",
" relerrorpct2 = abs((pub0-p_r_approach2)/pub0)*100\n",
" relerror01pct = abs((pub0-pub1)/pub1)*100\n",
" strrelerror1pct = \"%.4f\" % (relerrorpct1)\n",
" strrelerror2pct = \"%.4f\" % (relerrorpct2)\n",
" strrelerror01pct = \"%.4f\" % (relerror01pct)\n",
" print(\"Relative error between published iteration 0 vs it. 1: \"+strrelerror01pct+\"%\")\n",
" resultstring1 = \"Relative error between NRPyPN & published, approach 1: \"+strrelerror1pct+\"%\"\n",
" resultstring2 = \"Relative error between NRPyPN & published, approach 2: \"+strrelerror2pct+\"%\"\n",
" if relerrorpct1 > relerror01pct:\n",
" resultstring1 += \" < \"+strrelerror01pct+\"% <--- NOT GOOD!\"\n",
" resultstring2 += \" < \"+strrelerror01pct+\"% <--- NOT GOOD!\"\n",
" else:\n",
" resultstring1 += \" < \"+strrelerror01pct+\"% <--- GOOD AGREEMENT!\"\n",
" resultstring2 += \" < \"+strrelerror01pct+\"% <--- GOOD AGREEMENT!\"\n",
" print(resultstring1)\n",
" print(resultstring2+\"\\n\")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:52.220417Z",
"iopub.status.busy": "2021-03-07T17:36:52.205057Z",
"iopub.status.idle": "2021-03-07T17:36:57.529881Z",
"shell.execute_reply": "2021-03-07T17:36:57.530400Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"##################################################\n",
" Case: q=1, nonspinning, initial separation 12\n",
"##################################################\n",
"p_r = 0.000538330 <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)\n",
"p_r = 0.000539171 <- Result from NRPyPN, approach 1\n",
"p_r = 0.000539860 <- Result from NRPyPN, approach 2\n",
"p_r = 0.000468113 <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)\n",
"Relative error between published iteration 0 vs it. 1: 15.0000%\n",
"Relative error between NRPyPN & published, approach 1: 0.1562% < 15.0000% <--- GOOD AGREEMENT!\n",
"Relative error between NRPyPN & published, approach 2: 0.2843% < 15.0000% <--- GOOD AGREEMENT!\n",
"\n"
]
}
],
"source": [
"# 1. Let's consider the case:\n",
"# * Mass ratio q=1, chi1=chi2=(0,0,0), radial separation r=12\n",
"pub_result0 = 0.53833e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036\n",
"pub_result1 = 0.468113e-3 # Expected result 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036\n",
"qmassratio = 1.0 # must be >= 1\n",
"nr = 12. # Orbital separation\n",
"# Choose spins so that the total spin magnitude does not exceed 1.\n",
"nchi1x = +0.\n",
"nchi1y = +0.\n",
"nchi1z = +0.\n",
"nchi2x = +0.\n",
"nchi2y = +0.\n",
"nchi2z = +0.\n",
"\n",
"compare_pub_NPN(\"Case: q=1, nonspinning, initial separation 12\",pub_result0,pub_result1,\n",
" qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:36:57.558973Z",
"iopub.status.busy": "2021-03-07T17:36:57.543483Z",
"iopub.status.idle": "2021-03-07T17:37:03.242359Z",
"shell.execute_reply": "2021-03-07T17:37:03.242839Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"##################################################\n",
" Case: q=1.5, chi1z=-0.6, chi2z=0.6, initial separation 10.8\n",
"##################################################\n",
"p_r = 0.000699185 <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)\n",
"p_r = 0.000675944 <- Result from NRPyPN, approach 1\n",
"p_r = 0.000677020 <- Result from NRPyPN, approach 2\n",
"p_r = 0.000641051 <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)\n",
"Relative error between published iteration 0 vs it. 1: 9.0685%\n",
"Relative error between NRPyPN & published, approach 1: 3.3240% < 9.0685% <--- GOOD AGREEMENT!\n",
"Relative error between NRPyPN & published, approach 2: 3.1702% < 9.0685% <--- GOOD AGREEMENT!\n",
"\n"
]
}
],
"source": [
"# 2. Let's consider the case:\n",
"# * Mass ratio q=1.5, chi1= (0,0,-0.6); chi2=(0,0,0.6), radial separation r=10.8\n",
"pub_result0 = 0.699185e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036\n",
"pub_result1 = 0.641051e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036\n",
"qmassratio = 1.5 # must be >= 1\n",
"nr = 10.8 # Orbital separation\n",
"nchi1x = +0.\n",
"nchi1y = +0.\n",
"nchi1z = -0.6\n",
"nchi2x = +0.\n",
"nchi2y = +0.\n",
"nchi2z = +0.6\n",
"compare_pub_NPN(\"Case: q=1.5, chi1z=-0.6, chi2z=0.6, initial separation 10.8\",pub_result0,pub_result1,\n",
" qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:37:03.317819Z",
"iopub.status.busy": "2021-03-07T17:37:03.281876Z",
"iopub.status.idle": "2021-03-07T17:37:08.963778Z",
"shell.execute_reply": "2021-03-07T17:37:08.963221Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"##################################################\n",
" Case: q=4.0, chi1z=-0.8, chi2z=0.8, initial separation 11.0\n",
"##################################################\n",
"p_r = 0.000336564 <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)\n",
"p_r = 0.000251240 <- Result from NRPyPN, approach 1\n",
"p_r = 0.000251143 <- Result from NRPyPN, approach 2\n",
"p_r = 0.000247080 <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)\n",
"Relative error between published iteration 0 vs it. 1: 36.2166%\n",
"Relative error between NRPyPN & published, approach 1: 25.3515% < 36.2166% <--- GOOD AGREEMENT!\n",
"Relative error between NRPyPN & published, approach 2: 25.3804% < 36.2166% <--- GOOD AGREEMENT!\n",
"\n"
]
}
],
"source": [
"# 3. Let's consider the case:\n",
"# * Mass ratio q=4, chi1= (0,0,-0.8); chi2=(0,0,0.8), radial separation r=11\n",
"pub_result0 = 0.336564e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036\n",
"pub_result1 = 0.24708e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036\n",
"qmassratio = 4.0 # must be >= 1\n",
"nr = 11.0 # Orbital separation\n",
"nchi1x = +0.\n",
"nchi1y = +0.\n",
"nchi1z = -0.8\n",
"nchi2x = +0.\n",
"nchi2y = +0.\n",
"nchi2z = +0.8\n",
"compare_pub_NPN(\"Case: q=4.0, chi1z=-0.8, chi2z=0.8, initial separation 11.0\",pub_result0,pub_result1,\n",
" qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:37:08.987731Z",
"iopub.status.busy": "2021-03-07T17:37:08.977390Z",
"iopub.status.idle": "2021-03-07T17:37:14.572470Z",
"shell.execute_reply": "2021-03-07T17:37:14.571852Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"##################################################\n",
" Case: q=2.0, chi2x=-0.3535, chi2y=+0.3535, chi2z=+0.5, initial separation 10.8\n",
"##################################################\n",
"p_r = 0.000531374 <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)\n",
"p_r = 0.000545357 <- Result from NRPyPN, approach 1\n",
"p_r = 0.000542626 <- Result from NRPyPN, approach 2\n",
"p_r = 0.000448148 <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)\n",
"Relative error between published iteration 0 vs it. 1: 18.5711%\n",
"Relative error between NRPyPN & published, approach 1: 2.6314% < 18.5711% <--- GOOD AGREEMENT!\n",
"Relative error between NRPyPN & published, approach 2: 2.1175% < 18.5711% <--- GOOD AGREEMENT!\n",
"\n"
]
}
],
"source": [
"# 4. Let's consider the case:\n",
"# * Mass ratio q=2, chi1= (0,0,0); chi2=(−0.3535, 0.3535, 0.5), radial separation r=10.8\n",
"pub_result0 = 0.531374e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036\n",
"pub_result1 = 0.448148e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036\n",
"qmassratio = 2.0 # must be >= 1\n",
"nr = 10.8 # Orbital separation\n",
"nchi1x = +0.\n",
"nchi1y = +0.\n",
"nchi1z = +0.\n",
"nchi2x = -0.3535\n",
"nchi2y = +0.3535\n",
"nchi2z = +0.5\n",
"compare_pub_NPN(\"Case: q=2.0, chi2x=-0.3535, chi2y=+0.3535, chi2z=+0.5, initial separation 10.8\",pub_result0,pub_result1,\n",
" qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:37:14.647165Z",
"iopub.status.busy": "2021-03-07T17:37:14.611224Z",
"iopub.status.idle": "2021-03-07T17:37:20.332384Z",
"shell.execute_reply": "2021-03-07T17:37:20.332884Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"##################################################\n",
" \n",
"Case: q=8.0, chi1z=chi2z=+0.5, initial separation 11\n",
"\n",
"Note: This one is weird. Clearly the value in the table\n",
" has a typo, such that the p_r and p_t values\n",
" should be interchanged; p_t is about 20% the\n",
" next smallest value in the table, and the\n",
" parameters aren't that different. We therefore\n",
" assume that this is the case, and nonetheless\n",
" find agreement for p_t with the published result\n",
" to about 0.07%, and p_r to about 5%. Aiven that\n",
" the table values seem to be clearly wrong, this\n",
" level of agreement is an encouraging sign.\n",
"\n",
"##################################################\n",
"p_r = 0.000102969 <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)\n",
"p_r = 0.000097694 <- Result from NRPyPN, approach 1\n",
"p_r = 0.000097574 <- Result from NRPyPN, approach 2\n",
"p_r = 0.000139132 <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)\n",
"Relative error between published iteration 0 vs it. 1: 25.9919%\n",
"Relative error between NRPyPN & published, approach 1: 5.1225% < 25.9919% <--- GOOD AGREEMENT!\n",
"Relative error between NRPyPN & published, approach 2: 5.2399% < 25.9919% <--- GOOD AGREEMENT!\n",
"\n"
]
}
],
"source": [
"# 5. Let's consider the case:\n",
"# * Mass ratio q=8, chi1= (0, 0, 0.5); chi2=(0, 0, 0.5), radial separation r=11\n",
"pub_result0 = 0.102969e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036\n",
"pub_result1 = 0.139132e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036\n",
"qmassratio = 8.0 # must be >= 1\n",
"nr = 11.0 # Orbital separation\n",
"nchi1x = +0.\n",
"nchi1y = +0.\n",
"nchi1z = +0.5\n",
"nchi2x = +0.\n",
"nchi2y = +0.\n",
"nchi2z = +0.5\n",
"compare_pub_NPN(\"\"\"\n",
"Case: q=8.0, chi1z=chi2z=+0.5, initial separation 11\n",
"\n",
"Note: This one is weird. Clearly the value in the table\n",
" has a typo, such that the p_r and p_t values\n",
" should be interchanged; p_t is about 20% the\n",
" next smallest value in the table, and the\n",
" parameters aren't that different. We therefore\n",
" assume that this is the case, and nonetheless\n",
" find agreement for p_t with the published result\n",
" to about 0.07%, and p_r to about 5%. Aiven that\n",
" the table values seem to be clearly wrong, this\n",
" level of agreement is an encouraging sign.\n",
"\"\"\",pub_result0,pub_result1,\n",
" qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validationpython'></a>\n",
"\n",
"## Part 2.d: Validation against corresponding Python module \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validationpython}$$ \n",
"\n",
"Here we verify agreement between \n",
"\n",
"* the SymPy expressions generated in this notebook, and the corresponding Python module."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:37:20.407519Z",
"iopub.status.busy": "2021-03-07T17:37:20.371784Z",
"iopub.status.idle": "2021-03-07T17:37:29.145471Z",
"shell.execute_reply": "2021-03-07T17:37:29.144787Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PYTHON MODULE TEST PASSES\n"
]
}
],
"source": [
"def error(varname):\n",
" print(\"ERROR: When comparing Python module & notebook, \"+varname+\" was found not to match.\")\n",
" sys.exit(1)\n",
"\n",
"# Validation against Python module\n",
"import PN_p_r as pr\n",
"f_p_r( m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"pr.f_p_r(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)\n",
"\n",
"if sp.simplify(p_r - pr.p_r) != 0: error(\"pr.p_r\")\n",
"print(\"PYTHON MODULE TEST PASSES\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='latex_pdf_output'></a>\n",
"\n",
"# Part 3: 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",
"[PN-p_r.pdf](PN-p_r.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": 18,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:37:29.150203Z",
"iopub.status.busy": "2021-03-07T17:37:29.149537Z",
"iopub.status.idle": "2021-03-07T17:37:33.447885Z",
"shell.execute_reply": "2021-03-07T17:37:33.447186Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created PN-p_r.tex, and compiled LaTeX file to PDF file PN-p_r.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(\"PN-p_r\")"
]
}
],
"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
}