{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\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:** Validated \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": [ "\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": [ "\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": [ "\n", "\n", "# Part 1: Basic strategy for computing $p_r$ \\[Back to [top](#toc)\\]\n", "$$\\label{strategy}$$ \n", "\n", "\n", "\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": [ "\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": [ "\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", "\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": [ "\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": [ "\n", "\n", "# Part 2: Validation Tests \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation}$$ \n", "\n", "\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": [ "\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": [ "\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": [ "\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": [ "\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 }