Tutorial-ADM_Initial_Data-TOV.ipynb

{
 "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",
    "# [Polytropic TOV](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation) Initial Data\n",
    "\n",
    "## Authors: Phil Chang, Zach Etienne, & Leo Werneck\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module sets up initial data for a [TOV](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation) star in *spherical, isotropic coordinates*\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This module has been validated to exhibit convergence to zero of the Hamiltonian constraint violation at the expected order to the exact solution (see [start-to-finish TOV module](Tutorial-Start_to_Finish-BSSNCurvilinear-Setting_up_TOV_initial_data.ipynb) for full test). Note that convergence at the surface of the star is lower order due to the sharp drop to zero in $T^{\\mu\\nu}$.\n",
    "\n",
    "### NRPy+ Source Code for this module: [TOV/TOV_Solver.py](../edit/TOV/TOV_Solver.py)\n",
    "\n",
    "[comment]: <> (Introduction: TODO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows:\n",
    "\n",
    "1. [Step 1](#initializenrpy): Initialize core Python/NRPy+ modules\n",
    "1. [Step 2](#tov): The TOV Equations\n",
    "1. [Step 3](#code_validation): Code Validation against `TOV.TOV_Solver` NRPy+ module\n",
    "1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Initialize core Python/NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Import needed Python/NRPy+ modules\n",
    "import numpy as np                  # NumPy: A numerical methods module for Python\n",
    "import scipy.integrate as si        # SciPy: Python module for mathematics, science, and engineering applications\n",
    "import math, sys                    # Standard Python modules for math; multiplatform OS-level functions\n",
    "import TOV.Polytropic_EOSs as ppeos # NRPy+: Piecewise polytrope equation of state support"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='tov'></a>\n",
    "\n",
    "# Step 2: The TOV equations \\[Back to [top](#toc)\\]\n",
    "$$\\label{tov}$$\n",
    "\n",
    "The [TOV line element](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation) in terms of the *Schwarzschild coordinate* $r$ is written (in the $-+++$ form):\n",
    "$$\n",
    "ds^2 = - c^2 e^\\nu dt^2 + \\left(1 - \\frac{2Gm}{rc^2}\\right)^{-1} dr^2 + r^2 d\\Omega^2,\n",
    "$$\n",
    "where $m(r)$ is the mass-energy enclosed at a given $r$, and is equal to the total star's mass outside the stellar radius $r=R$.\n",
    "\n",
    "In terms of the *isotropic coordinate* $\\bar{r}$ with $G=c=1$ (i.e., the coordinate system and units we'd prefer to use), the ($-+++$ form) line element is written:\n",
    "$$\n",
    "ds^2 = - e^{\\nu} dt^2 + e^{4\\phi} \\left(d\\bar{r}^2 + \\bar{r}^2 d\\Omega^2\\right),\n",
    "$$\n",
    "where $\\phi$ here is the *conformal factor*.\n",
    "\n",
    "Setting components of the above line element equal to one another, we get (in $G=c=1$ units):\n",
    "\n",
    "\\begin{align}\n",
    "r^2 &= e^{4\\phi} \\bar{r}^2 \\implies e^{4\\phi} = \\frac{r^2}{\\bar{r}^2} \\\\\n",
    "\\left(1 - \\frac{2m}{r}\\right)^{-1} dr^2 &= e^{4\\phi} d\\bar{r}^2 \\\\\n",
    "\\implies \\frac{d\\bar{r}(r)}{dr} &= \\left(1 - \\frac{2m}{r} \\right)^{-1/2} \\frac{\\bar{r}(r)}{r}.\n",
    "\\end{align}\n",
    "\n",
    "The TOV equations provide radial ODEs for the pressure and $\\nu$ (from [the Wikipedia article on the TOV solution](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation)):\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dP}{dr} &= - \\frac{1}{r} \\left( \\frac{\\rho + P}{2} \\right) \\left(\\frac{2 m}{r} + 8 \\pi r^2 P\\right) \\left(1 - \\frac{2 m}{r}\\right)^{-1} \\\\\n",
    "\\frac{d \\nu}{d r} &= \\frac{1}{r}\\left(1 - \\frac{2 m}{r}\\right)^{-1} \\left(\\frac{2 m}{r} + 8 \\pi r^2 P\\right) \\\\\n",
    "\\end{align}\n",
    "\n",
    "Assuming a polytropic equation of state, which relates the pressure $P$ to the baryonic rest-mass density $\\rho_B$,\n",
    "\n",
    "$$\n",
    "P(\\rho_B) = K \\rho_B^\\Gamma,\n",
    "$$\n",
    "the specific internal energy will be given by\n",
    "$$\n",
    "\\epsilon = \\frac{P}{\\rho_B (\\Gamma - 1)},\n",
    "$$\n",
    "\n",
    "so the total mass-energy density $\\rho$ is given by\n",
    "$$\n",
    "\\rho = \\rho_B (1 + \\epsilon).\n",
    "$$\n",
    "\n",
    "Given this, the mass-energy $m(r)$ density is the solution to the ODE:\n",
    "$$\n",
    "\\frac{dm(r)}{dr} = 4\\pi r^2 \\rho(r)\n",
    "$$\n",
    "\n",
    "Thus the full set of ODEs that need to be solved is given by\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\begin{matrix}\n",
    "\\frac{dP}{dr} &=& - \\frac{1}{r} \\left( \\frac{\\rho + P}{2} \\right) \\left(\\frac{2 m}{r} + 8 \\pi r^2 P\\right) \\left(1 - \\frac{2 m}{r}\\right)^{-1} \\\\\n",
    "\\frac{d \\nu}{d r} &=& \\frac{1}{r}\\left(1 - \\frac{2 m}{r}\\right)^{-1} \\left(\\frac{2 m}{r} + 8 \\pi r^2 P\\right) \\\\\n",
    "\\frac{m(r)}{dr} &=& 4\\pi r^2 \\rho(r) \\\\\n",
    "\\frac{d\\bar{r}(r)}{dr} &=& \\left(1 - \\frac{2m}{r} \\right)^{-1/2} \\frac{\\bar{r}(r)}{r}\n",
    "\\end{matrix}\n",
    "}\\ .\n",
    "$$\n",
    "\n",
    "The following code solves these equations, and was largely written by Phil Chang."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1256 1256 1256 1256 1256 1256\n",
      "Just generated a TOV star with R_Schw = 9.566044579232513e-01 , M = 1.405030336771405e-01 , M/R_Schw = 1.468768334847266e-01 .\n"
     ]
    }
   ],
   "source": [
    "# Step 2: The TOV equations\n",
    "\n",
    "## TOV SOLVER FOR SINGLE AND PIECEWISE POLYTROPES\n",
    "## Authors: Phil Chang, Zachariah B. Etienne, Leo Werneck\n",
    "\n",
    "# Full documentation for this module may be found in the NRPy+ tutorial Jupyter notebook:\n",
    "#  Tutorial-Start_to_Finish-BSSNCurvilinear-Setting_up_TOV_initial_data.ipynb\n",
    "\n",
    "\n",
    "# Inputs:\n",
    "# * Output data file name\n",
    "# * rho_baryon_central, the central density of the TOV star.\n",
    "# * n, the polytropic equation of state index. n=1 models cold, degenerate neutron star matter.\n",
    "# * K_Polytrope, the polytropic constant.\n",
    "# * Verbose output toggle (default = True)\n",
    "\n",
    "# Output: An initial data file (default file name = \"outputTOVpolytrope.txt\") that well\n",
    "#         samples the (spherically symmetric) solution both inside and outside the star.\n",
    "#         It is up to the initial data module to perform the 1D interpolation to generate\n",
    "#         the solution at arbitrary radius. The file has the following columns:\n",
    "# Column 1: Schwarzschild radius\n",
    "# Column 2: rho(r), *total* mass-energy density (as opposed to baryonic rest-mass density)\n",
    "# Column 3: P(r), Pressure\n",
    "# Column 4: m(r), mass enclosed\n",
    "# Column 5: e^{nu(r)}, g_{tt}(r)\n",
    "# Column 6: e^{4 phi(r)}, conformal factor g_{rr}(r)\n",
    "# Column 7: rbar(r), Isotropic radius\n",
    "\n",
    "# rbar refers to the isotropic radius, and\n",
    "# R_Schw refers to the Schwarzschild radius\n",
    "\n",
    "def TOV_Solver(eos,\n",
    "               outfile = \"outputTOVpolytrope.txt\",\n",
    "               rho_baryon_central = 0.129285,\n",
    "               verbose = True,\n",
    "               return_M_and_RSchw = False,\n",
    "               accuracy = \"medium\",\n",
    "               integrator_type = \"default\",\n",
    "               no_output_File = False,\n",
    "               pressure_renormalization=1): # reset the pressure to stellar oscillations studies\n",
    "\n",
    "    def TOV_rhs(r_Schw, y) :\n",
    "    # In \\tilde units\n",
    "    #\n",
    "        P    = y[0]\n",
    "        m    = y[1]\n",
    "        # nu   = y[2] # nu is not needed as input into TOV_rhs\n",
    "        rbar = y[3]\n",
    "\n",
    "        # Compute rho_b and eps_cold, to be used below\n",
    "        # to compute rho_(total)\n",
    "        rho_baryon, eps_cold = ppeos.Polytrope_EOS__compute_rhob_and_eps_cold_from_P_cold(eos,P)\n",
    "\n",
    "#         with open(\"rhob_P_cold_and_eps_cold.dat\",\"a+\") as file:\n",
    "#             file.write(str(r_Schw).format(\"%.15e\")+\"  \"+str(rho_baryon).format(\"%.15e\")+\"  \"+str(P).format(\"%.15e\")+\"  \"+str(eps_cold).format(\"%.15e\")+\"\\n\")\n",
    "\n",
    "        # Compute rho, the *total* mass-energy density:\n",
    "        # .------------------------------.\n",
    "        # | rho = (1 + eps)*rho_(baryon) |\n",
    "        # .------------------------------.\n",
    "        # with eps = eps_cold, for the initial data.\n",
    "        rho = (1.0 + eps_cold)*rho_baryon\n",
    "\n",
    "#         m = 4*math.pi/3. * rho*r_Schw**3\n",
    "        if( r_Schw < 1e-4 or m <= 0.):\n",
    "            # From https://github.com/natj/tov/blob/master/tov.py#L33:\n",
    "            # dPdr = -cgs.G*(eden + P/cgs.c**2)*(m + 4.0*pi*r**3*P/cgs.c**2)\n",
    "            # dPdr = dPdr/(r*(r - 2.0*cgs.G*m/cgs.c**2))\n",
    "            dPdrSchw = -(rho + P)*(4.*math.pi/3.*r_Schw*rho + 4.*math.pi*r_Schw*P)/(1.-8.*math.pi*rho*r_Schw*r_Schw)\n",
    "            drbardrSchw = 1./(1. - 8.*math.pi*rho*r_Schw*r_Schw)**0.5\n",
    "        else:\n",
    "            dPdrSchw = -(rho + P)*(m + 4.*math.pi*r_Schw**3*P)/(r_Schw*r_Schw*(1.-2.*m/r_Schw))\n",
    "            drbardrSchw = 1./(1. - 2.*m/r_Schw)**0.5*rbar/r_Schw\n",
    "\n",
    "        dmdrSchw  =  4.*math.pi*r_Schw*r_Schw*rho\n",
    "        dnudrSchw = -2./(P + rho)*dPdrSchw\n",
    "        return [dPdrSchw, dmdrSchw, dnudrSchw, drbardrSchw]\n",
    "\n",
    "    def integrateStar( eos, P, dumpData = False ):\n",
    "\n",
    "        if accuracy == \"medium\":\n",
    "            min_step_size = 1e-5\n",
    "            max_step_size = 1e-2\n",
    "            integrator    = 'dop853'\n",
    "        elif accuracy == \"low\":\n",
    "            min_step_size = 1e-3\n",
    "            max_step_size = 1e-1\n",
    "            integrator    = 'dopri5'\n",
    "        elif accuracy == \"verylow\":\n",
    "            min_step_size = 1e-1\n",
    "            max_step_size = 5e-1\n",
    "            integrator    = 'dopri5'\n",
    "        elif accuracy == \"high\":\n",
    "            min_step_size = 1e-5\n",
    "            max_step_size = 1e-5\n",
    "            integrator    = 'dop853'\n",
    "        elif accuracy == \"veryhigh\":\n",
    "            min_step_size = 1e-7\n",
    "            max_step_size = 1e-6\n",
    "            integrator    = 'dop853'\n",
    "        else:\n",
    "            print(\"Unknown accuracy option: \"+str(accuracy))\n",
    "\n",
    "        if integrator_type == \"default\":\n",
    "            pass\n",
    "        else:\n",
    "            integrator = integrator_type\n",
    "\n",
    "        integrator = si.ode(TOV_rhs).set_integrator(integrator)#,rtol=1e-4,atol=1e-4)\n",
    "#         integrator = si.ode(TOV_rhs).set_integrator('dopri5',rtol=1e-4)\n",
    "        y0 = [P, 0., 0., 0.]\n",
    "        r_Schw = 0.\n",
    "        integrator.set_initial_value(y0,r_Schw)\n",
    "        dr_Schw = min_step_size\n",
    "        P = y0[0]\n",
    "\n",
    "        PArr      = []\n",
    "        r_SchwArr = []\n",
    "        mArr      = []\n",
    "        nuArr     = []\n",
    "        rbarArr   = []\n",
    "\n",
    "        while integrator.successful() and P > 1e-19*y0[0] :\n",
    "            P, m, nu, rbar = integrator.integrate(r_Schw + dr_Schw)\n",
    "            # Update the value of r_Schw to the latest integrated value\n",
    "            r_Schw += dr_Schw\n",
    "\n",
    "            dPdrSchw, dmdrSchw, dnudrSchw, drbardrSchw = TOV_rhs( r_Schw, [P,m,nu,rbar])\n",
    "            dr_Schw = 0.1*min(abs(P/dPdrSchw), abs(m/dmdrSchw))\n",
    "            dr_Schw = min(dr_Schw, max_step_size)\n",
    "            PArr.append(P)\n",
    "            r_SchwArr.append(r_Schw)\n",
    "            mArr.append(m)\n",
    "            nuArr.append(nu)\n",
    "            rbarArr.append(rbar)\n",
    "\n",
    "        M = mArr[-1]\n",
    "        R_Schw = r_SchwArr[-1]\n",
    "\n",
    "        if no_output_File == True:\n",
    "            return R_Schw, M\n",
    "\n",
    "        # Apply integration constant to ensure rbar is continuous across TOV surface\n",
    "        for ii in range(len(rbarArr)):\n",
    "            rbarArr[ii] *= 0.5*(np.sqrt(R_Schw*(R_Schw - 2.0*M)) + R_Schw - M) / rbarArr[-1]\n",
    "\n",
    "        nuArr_np = np.array(nuArr)\n",
    "        # Rescale solution to nu so that it satisfies BC: exp(nu(R))=exp(nutilde-nu(r=R)) * (1 - 2m(R)/R)\n",
    "        #   Thus, nu(R) = (nutilde - nu(r=R)) + log(1 - 2*m(R)/R)\n",
    "        nuArr_np = nuArr_np - nuArr_np[-1] + math.log(1.-2.*mArr[-1]/r_SchwArr[-1])\n",
    "\n",
    "        r_SchwArrExtend_np = 10.**(np.arange(0.01,5.0,0.01))*r_SchwArr[-1]\n",
    "\n",
    "        r_SchwArr.extend(r_SchwArrExtend_np)\n",
    "        mArr.extend(r_SchwArrExtend_np*0. + M)\n",
    "        PArr.extend(r_SchwArrExtend_np*0.)\n",
    "        exp2phiArr_np = np.append( np.exp(nuArr_np), 1. - 2.*M/r_SchwArrExtend_np)\n",
    "        nuArr.extend(np.log(1. - 2.*M/r_SchwArrExtend_np))\n",
    "        rbarArr.extend( 0.5*(np.sqrt(r_SchwArrExtend_np**2 - 2.*M*r_SchwArrExtend_np) + r_SchwArrExtend_np - M) )\n",
    "\n",
    "        #phiArr_np = np.append( np.exp(nuArr_np), 1. - 2.*M/r_SchwArrExtend_np)\n",
    "\n",
    "        # Appending to a Python array does what one would reasonably expect.\n",
    "        #   Appending to a numpy array allocates space for a new array with size+1,\n",
    "        #   then copies the data over... over and over... super inefficient.\n",
    "        r_SchwArr_np     = np.array(r_SchwArr)\n",
    "        PArr_np          = np.array(PArr)\n",
    "        rho_baryonArr_np = np.array(PArr) # This is just to initialize the array\n",
    "        for j in range(len(PArr_np)):\n",
    "            # Compute rho_b from P\n",
    "            rho_baryonArr_np[j] = ppeos.Polytrope_EOS__compute_rhob_from_P_cold(eos,PArr_np[j])\n",
    "\n",
    "        mArr_np               = np.array(mArr)\n",
    "        rbarArr_np            = np.array(rbarArr)\n",
    "        confFactor_exp4phi_np = (r_SchwArr_np/rbarArr_np)**2\n",
    "\n",
    "        # Compute the *total* mass-energy density (as opposed to the *baryonic* mass density)\n",
    "        rhoArr_np = []\n",
    "        for i in range(len(PArr)):\n",
    "            rho_baryon, eps_cold = ppeos.Polytrope_EOS__compute_rhob_and_eps_cold_from_P_cold(eos,PArr[i])\n",
    "            rho = (1.0 + eps_cold ) * rho_baryon\n",
    "            rhoArr_np.append(rho)\n",
    "\n",
    "        if verbose:\n",
    "            print(len(r_SchwArr_np),len(rhoArr_np),len(rho_baryonArr_np),len(PArr_np),len(mArr_np),len(exp2phiArr_np))\n",
    "\n",
    "        PArr_np *= pressure_renormalization # set for pressure renormalization studies\n",
    "\n",
    "        # Special thanks to Leonardo Werneck for pointing out this issue with zip()\n",
    "        if sys.version_info[0] < 3:\n",
    "            np.savetxt(outfile, zip(r_SchwArr_np,rhoArr_np,rho_baryonArr_np,PArr_np,mArr_np,exp2phiArr_np,confFactor_exp4phi_np,rbarArr_np),\n",
    "                       fmt=\"%.15e\")\n",
    "        else:\n",
    "            np.savetxt(outfile, list(zip(r_SchwArr_np,rhoArr_np,rho_baryonArr_np,PArr_np,mArr_np,exp2phiArr_np,confFactor_exp4phi_np,rbarArr_np)),\n",
    "                       fmt=\"%.15e\")\n",
    "\n",
    "        return R_Schw, M\n",
    "\n",
    "    # Set initial condition from rho_baryon_central\n",
    "    P_initial_condition = ppeos.Polytrope_EOS__compute_P_cold_from_rhob(eos, rho_baryon_central)\n",
    "\n",
    "    # Integrate the initial condition\n",
    "    R_Schw_TOV, M_TOV = integrateStar(eos, P_initial_condition, True)\n",
    "    if verbose:\n",
    "        print(\"Just generated a TOV star with R_Schw = %.15e , M = %.15e , M/R_Schw = %.15e .\" %(R_Schw_TOV,M_TOV,(M_TOV / R_Schw_TOV)))\n",
    "\n",
    "    if return_M_and_RSchw:\n",
    "        return M_TOV, R_Schw_TOV\n",
    "\n",
    "############################\n",
    "# Single polytrope example #\n",
    "############################\n",
    "# Set neos = 1 (single polytrope)\n",
    "neos = 1\n",
    "\n",
    "# Set rho_poly_tab (not needed for a single polytrope)\n",
    "rho_poly_tab = []\n",
    "\n",
    "# Set Gamma_poly_tab\n",
    "Gamma_poly_tab = [2.0]\n",
    "\n",
    "# Set K_poly_tab0\n",
    "K_poly_tab0 = 1. # ZACH NOTES: CHANGED FROM 100.\n",
    "\n",
    "# Set the eos quantities\n",
    "eos = ppeos.set_up_EOS_parameters__complete_set_of_input_variables(neos,rho_poly_tab,Gamma_poly_tab,K_poly_tab0)\n",
    "\n",
    "# Set initial condition (Pressure computed from central density)\n",
    "rho_baryon_central  = 0.129285\n",
    "\n",
    "M_TOV, R_Schw_TOV = TOV_Solver(eos,outfile=\"outputTOVpolytrope.txt\",rho_baryon_central=0.129285,verbose = True,\n",
    "                               return_M_and_RSchw=True,accuracy=\"medium\",integrator_type=\"default\",\n",
    "                               pressure_renormalization=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 3: Code Validation \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for these TOV initial data between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [TOV.TOV_Solver](../edit/TOV/TOV_Solver.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1256 1256 1256 1256 1256 1256\n",
      "Just generated a TOV star with\n",
      "* M        = 1.405030336771405e-01 ,\n",
      "* R_Schw   = 9.566044579232513e-01 ,\n",
      "* R_iso    = 8.100085557410308e-01 ,\n",
      "* M/R_Schw = 1.468768334847266e-01 \n",
      "\n",
      "TOV initial data test PASSED.\n"
     ]
    }
   ],
   "source": [
    "# Step 3: Code Validation against TOV.TOV_Solver module\n",
    "\n",
    "import filecmp\n",
    "\n",
    "import TOV.TOV_Solver as TOV\n",
    "TOV.TOV_Solver(eos,\n",
    "               outfile=\"outputTOVpolytrope-validation.txt\",\n",
    "               rho_baryon_central=0.129285,\n",
    "               verbose = True,\n",
    "               accuracy=\"medium\",\n",
    "               integrator_type=\"default\",\n",
    "               no_output_File = False)\n",
    "\n",
    "if filecmp.cmp('outputTOVpolytrope.txt',\n",
    "               'outputTOVpolytrope-validation.txt') == False:\n",
    "    print(\"ERROR: TOV initial data test FAILED!\")\n",
    "    sys.exit(1)\n",
    "else:\n",
    "    print(\"TOV initial data test PASSED.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 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",
    "[Tutorial-ADM_Initial_Data-TOV](Tutorial-ADM_Initial_Data-TOV.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": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ADM_Initial_Data-TOV.tex, and compiled LaTeX file to PDF\n",
      "    file Tutorial-ADM_Initial_Data-TOV.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ADM_Initial_Data-TOV\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}