{
 "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: Validated Post-Newtonian Expressions \n",
    "\n",
    "## This notebook validates Post-Newtonian expressions for input into Wolfram Mathematica, SymPy, or Highly-Optimized C Codes. It also provides low-eccentricity binary black hole initial data parameters for Bowen-York initial data (e.g., `TwoPunctures`).\n",
    "\n",
    "### Lead author: [Zachariah B. Etienne](https://etienneresearch.com) $\\leftarrow$ Please feel free to email contributions, comments, revisions, or errata!\n",
    "\n",
    "### Introduction\n",
    "\n",
    "[Post-Newtonian theory](https://en.wikipedia.org/wiki/Post-Newtonian_expansion) results in some of the longest and most complex mathematical expressions ever derived by humanity.\n",
    "\n",
    "These expressions form the core of most gravitational wave data analysis codes, but generally, the expressions are written in a format that is either inaccessible by others (i.e., closed-source) or written directly in e.g., C code and thus incompatible with symbolic algebra packages like [Wolfram Mathematica](https://www.wolfram.com/mathematica/) or the Python-based [SymPy](https://www.sympy.org/). \n",
    "\n",
    "Once in a symbolic algebra package, these expressions could be manipulated, extended, and output as *more* optimized C codes (e.g., using SymPy/NRPy+, thus speeding up gravitational wave data analysis).\n",
    "\n",
    "This repository aims to provide a trusted source for validated Post-Newtonian expressions useful for gravitational wave astronomy, using the open-source [SymPy](https://www.sympy.org/) computer algebra system so that expressions can be output into [Wolfram Mathematica](https://www.wolfram.com/mathematica/) or highly optimized C codes using the SymPy-based [NRPy+](http://nrpyplus.net/).\n",
    "\n",
    "***If you are unfamiliar with using Jupyter Notebooks, first review the official [Jupyter Notebook Basics Guide](https://nbviewer.jupyter.org/github/jupyter/notebook/blob/master/docs/source/examples/Notebook/Notebook%20Basics.ipynb).***\n",
    "\n",
    "### PART 0: Basic Functions for Expediting Equation Imports\n",
    "\n",
    "+ [NRPyPN_shortcuts](NRPyPN_shortcuts.ipynb)\n",
    "\n",
    "### PART 1: Post-Newtonian Hamiltonian for Binary Black Holes\n",
    "\n",
    "The Post-Newtonian Hamiltonian $H$ can be written as a sum of contributions:\n",
    "\\begin{array}\n",
    "\\ H = H_{\\rm Orb, NS} + H_{\\rm SO} + H_{\\rm SS} + H_{\\rm SSS},\n",
    "\\end{array}\n",
    "where \n",
    "* $H_{\\rm Orb, NS}$ is the non-spinning, purely orbital contribution\n",
    "* $H_{\\rm SO}$ accounts for spin-orbit coupling\n",
    "* $H_{\\rm SS}$ accounts for spin-spin contributions\n",
    "* $H_{\\rm SSS}$ accounts for spin-spin-spin contributions\n",
    "\n",
    "Click on any term below to open a notebook implementing that Hamiltonian component.\n",
    "#### $H_{\\rm Orb, NS}$, up to and including third post-Newtonian order (3PN)\n",
    "+ [$H_{\\rm Orb, NS}$](PN-Hamiltonian-Nonspinning.ipynb), as summarized in [Buonanno, Chen, and Damour (2006)](https://arxiv.org/abs/gr-qc/0508067) (see references therein for sources)\n",
    "\n",
    "#### $H_{\\rm SO}$, up to and including 3.5 post-Newtonian order (3.5PN)\n",
    "+ [The $H_{\\rm SO}$ notebook](PN-Hamiltonian-Spin-Orbit.ipynb) contains all spin-orbit coupling terms up to and including 3.5PN order:\n",
    "    * 1.5PN order (i.e., $H_{\\rm SO, 1.5PN}$), as summarized in [Buonanno, Chen, and Damour (2006)](https://arxiv.org/abs/gr-qc/0508067) (see references therein for sources)\n",
    "    * 2.5PN order (i.e., $H_{\\rm SO, 2.5PN}$), as derived in [Damour, Jaranowski, and Sch&#228;fer (2008)](https://arxiv.org/abs/0711.1048)\n",
    "    * 3.5PN order (i.e., $H_{\\rm SO, 3.5PN}$), as derived in [Hartung and Steinhoff (2011)](https://arxiv.org/abs/1104.3079)\n",
    "\n",
    "#### $H_{\\rm SS}$, up to and including third post-Newtonian order (3PN)\n",
    "+ [The $H_{\\rm SS}$ notebook](PN-Hamiltonian-Spin-Spin.ipynb) contains all spin-orbit coupling terms up to and including 3PN order\n",
    "    + 2PN order (i.e., $H_{S_1,S_2,{\\rm 2PN}}+H_{S_1^2,{\\rm 2PN}}+H_{S_2^2,{\\rm 2PN}}$), as summarized in [Buonanno, Chen, and Damour (2006)](https://arxiv.org/abs/gr-qc/0508067) (see references therein for sources)\n",
    "    + $S_1,S_2$ term at 3PN order (i.e., $H_{S_1,S_2,{\\rm 3PN}}$), as derived in [Steinhoff, Hergt, and Sch&#228;fer (2008a)](https://arxiv.org/abs/0712.1716)\n",
    "    + $S_1^2$ term at 3PN order (i.e., $H_{S_1^2,{\\rm 3PN}}$), as derived in [Steinhoff, Hergt, and Sch&#228;fer (2008b)](https://arxiv.org/abs/0809.2200)\n",
    "\n",
    "#### $H_{\\rm SSS}$, up to and including third post-Newtonian order (3PN)\n",
    "+ [$H_{SSS,{\\rm 3PN}}$](PN-Hamiltonian-SSS.ipynb), as derived in [Levi and Steinhoff (2015)](https://arxiv.org/abs/1410.2601)\n",
    "\n",
    "### Part 2: $\\frac{dE_{\\rm GW}}{dt}$, the gravitational wave flux, and $\\frac{dM}{dt}$, the tidal energy injected into the black holes\n",
    "\n",
    "+ [Gravitational-wave and mass fluxes $\\frac{dE_{\\rm GW}}{dt}$ and $\\frac{dM}{dt}$](PN-dE_GW_dt_and_dM_dt.ipynb), as reviewed by \n",
    "    * [Blanchet (2014)](https://link.springer.com/content/pdf/10.12942/lrr-2014-2.pdf), for the nonspinning terms\n",
    "    * [Ossokine *et al* (2015)](https://arxiv.org/abs/1502.01747), including precessing spin terms\n",
    "    * [Ajith *et al* (2007)](https://arxiv.org/abs/0709.0093), including tidal-heating injected energy into the black holes [Alvi (2001)](https://arxiv.org/abs/gr-qc/0107080)\n",
    "    \n",
    "### Part 3: Quasicircular Orbital Parameters $p_t$ and $p_r$\n",
    "\n",
    "+ [Tangential component of momentum $p_t$](PN-p_t.ipynb), as derived by \n",
    "    * [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036) up to and including 3.5PN order\n",
    "        * ... and validated against [Healy, Lousto, Nakano, and Zlochower (2017)](https://arxiv.org/abs/1702.00872), who derive the expression up to and including 3PN order\n",
    "+ [Orbital angular frequency $M\\Omega$](PN-MOmega.ipynb), as derived by \n",
    "    * [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036) up to and including 3.5PN order\n",
    "        * ... and validated against [Healy, Lousto, Nakano, and Zlochower (2017)](https://arxiv.org/abs/1702.00872), who derive the expression up to and including 3PN order\n",
    "+ [Radial component of momentum $p_r$](PN-p_r.ipynb), as derived by \n",
    "    * [Healy, Lousto, Nakano, and Zlochower (2017)](https://arxiv.org/abs/1702.00872) up to and including 3PN order\n",
    "    * [Ramos-Buades, Husa, and Pratten (2018)](https://arxiv.org/abs/1810.00036) up to and including 3.5PN order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fast $p_t$ and $p_r$ for binary black hole initial data\n",
    "\n",
    "**Assumptions:**\n",
    "\n",
    "The following code cell assumes the black holes are \n",
    "\n",
    "* situated initially on the $x$-axis, with black hole 1 at $x>0$, and black hole 2 at $x<0$\n",
    "* orbiting instantaneously in the $x$-$y$ plane, with the center-of-mass at $x=y=z=0$. \n",
    "\n",
    "Note that other initial coordinate configurations are possible; just permute the coordinates consistently with a right-handed coordinate system. For example,\n",
    "\n",
    "* perform the conversion $(x,y,z)\\to(z,x,y)$ if one desires the black holes to be on the $z-axis$ instantaneously orbiting on the $z$-$x$ plane; or \n",
    "* perform the conversion $(x,y,z)\\to(y,z,x)$ if one desires the black holes to be on the $y-axis$ instantaneously orbiting on the $y$-$z$ plane.\n",
    "\n",
    "**Instructions:** Set the 8 parameters for the binary black hole system in the code cell below, which make the above assumptions. Then after setting the parameters press *Shift+Enter*. The tangential and radial components of the orbital momentum vector will be computed and output.\n",
    "\n",
    "**If using `TwoPunctures`:** If you plan to use these parameters with TwoPunctures, you will need to set\n",
    "\n",
    "* `TwoPunctures::give_bare_mass = no`\n",
    "* `TwoPunctures::target_m_minus = ` $q/(1+q)$\n",
    "* `TwoPunctures::target_m_plus = ` $1/(1+q)$\n",
    "* `TwoPunctures::par_P_plus[0] = -|P_r|` (from code cell output below)\n",
    "* `TwoPunctures::par_P_plus[1] = +|P_t|` (from code cell output below)\n",
    "* `TwoPunctures::par_P_minus[0] = +|P_r|` (from code cell output below)\n",
    "* `TwoPunctures::par_P_minus[1] = -|P_t|` (from code cell output below)\n",
    "* `TwoPunctures::par_S_plus[0] = nchi1x*TwoPunctures::target_m_plus*TwoPunctures::target_m_plus` (from input into code cell below)\n",
    "* `TwoPunctures::par_S_plus[1] = nchi1y*TwoPunctures::target_m_plus*TwoPunctures::target_m_plus` (from input into code cell below)\n",
    "* `TwoPunctures::par_S_plus[2] = nchi1z*TwoPunctures::target_m_plus*TwoPunctures::target_m_plus` (from input into code cell below)\n",
    "* *similarly* for `TwoPunctures::par_S_minus[]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-16T04:40:33.745429Z",
     "iopub.status.busy": "2021-09-16T04:40:33.744975Z",
     "iopub.status.idle": "2021-09-16T04:40:37.681224Z",
     "shell.execute_reply": "2021-09-16T04:40:37.681632Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "If the black holes are situated on the x-axis, and orbiting instantaneously in the\n",
      "   xy-plane, then one could have P_y=+|P_t| and P_x=-|P_r| for the x>0 black hole,\n",
      "   and P_y=-|P_t| and P_x=+|P_r| for the x<0 black hole.  Other configurations are\n",
      "   possible provided one permutes the coordinates consistently with a right-handed\n",
      "   coordinate system.\n",
      "|P_t| =  0.0813467633310249\n",
      "|P_r| =  0.000553891860334773\n"
     ]
    }
   ],
   "source": [
    "#### INPUT PARAMETERS:\n",
    "qmassratio =  1.75  # m2/m1; by convention must be >= 1\n",
    "nr         =  10.8  # Orbital separation\n",
    "# Dimensionless spin parameters of each black hole, assuming\n",
    "#   the black holes are initially orbiting on the xy-plane\n",
    "#   and situated on the x-axis, with the center-of-mass at\n",
    "#   x=y=z=0. Black hole 1 is at x>0; and black hole 2 at x<0.\n",
    "nchi1x     = +0.   # chi_{1x}, x-component of dimensionless spin for black hole 1 (instantaneously at x>0)\n",
    "nchi1y     = +0.   # chi_{1y}, y-component of dimensionless spin for black hole 1 (instantaneously at x>0)\n",
    "nchi1z     = +0.6  # chi_{1z}, z-component of dimensionless spin for black hole 1 (instantaneously at x>0)\n",
    "nchi2x     = +0.   # chi_{2x}, x-component of dimensionless spin for black hole 2 (instantaneously at x<0)\n",
    "nchi2y     = +0.   # chi_{2y}, y-component of dimensionless spin for black hole 2 (instantaneously at x<0)\n",
    "nchi2z     = +0.6  # chi_{2z}, z-component of dimensionless spin for black hole 2 (instantaneously at x<0)\n",
    "\n",
    "##########################################################\n",
    "#### DON'T TOUCH; see output after running this cell. ####\n",
    "from NRPyPN_shortcuts import eval__P_t__and__P_r\n",
    "\n",
    "# Compute nPt, the tangential component of orbital momentum vector (output = numerical), as well as\n",
    "#         nPr, the radial component of orbital momentum vector (output = numerical)\n",
    "nPt, nPr = eval__P_t__and__P_r(qmassratio, nr,\n",
    "                               nchi1x, nchi1y, nchi1z,\n",
    "                               nchi2x, nchi2y, nchi2z)\n",
    "\n",
    "print(\"If the black holes are situated on the x-axis, and orbiting instantaneously in the\")\n",
    "print(\"   xy-plane, then one could have P_y=+|P_t| and P_x=-|P_r| for the x>0 black hole,\")\n",
    "print(\"   and P_y=-|P_t| and P_x=+|P_r| for the x<0 black hole.  Other configurations are\")\n",
    "print(\"   possible provided one permutes the coordinates consistently with a right-handed\")\n",
    "print(\"   coordinate system.\")\n",
    "print(\"|P_t| = \",nPt)\n",
    "print(\"|P_r| = \",nPr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "## Part 4: 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.pdf](NRPyPN.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": null,
   "metadata": {},
   "outputs": [],
   "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(\"NRPyNP\")"
   ]
  }
 ],
 "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
}