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",
    "# [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 notebook establishes initial data for a [TOV](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation) star in spherical, isotropic coordinates. The notebook presents a module for setting up ADM metric quantities and the stress-energy tensor $T^{\\mu\\nu}$ from interpolated TOV data and demonstrates reliability through validation against the Python `TOV.TOV_Ccodegen_library` module. Detailed steps for reading and interpolating TOV data files, while avoiding the Gibbs phenomenon at the stellar surface, are outlined.\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-TOV_initial_data.ipynb) for a 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](#c_code_generation): C-code routines (library) for reading TOV data files\n",
    "    1. [Step 4.a](#id_persist): Declare `ID_persist`, the C `struct` data type that will store data read from file\n",
    "    1. [Step 4.b](#read_data_file): `TOV_read_data_file_set_ID_persist()`: Read TOV data file and store data to the `ID_persist` struct\n",
    "    1. [Step 4.c](#interp_data_file): `TOV_interpolate_1D()`: Interpolate TOV data to any desired distance from the center of the star\n",
    "    1. [Step 4.d](#tov_id_function): `TOV_ID_function()`: Driver routine for setting up ADM metric quantities and $T^{\\mu\\nu}$ from interpolated TOV data\n",
    "        1. [Step 4.d.i](#adm_ito_tov): ADM metric quantities, written in terms of the TOV solution\n",
    "        1. [Step 4.d.ii](#tmunu_ito_tov): Stress-energy tensor $T^{\\mu\\nu}$, written in terms of the TOV solution\n",
    "        1. [Step 4.d.iii](#tov_id_function_itself): The `TOV_ID_function()` itself\n",
    "1. [Step 5](#ccodegen_library_validation): Validate C code generation against `TOV.TOV_Ccodegen_library` Python module \n",
    "1. [Step 6](#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": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:06:26.375685Z",
     "iopub.status.busy": "2021-03-07T17:06:26.374015Z",
     "iopub.status.idle": "2021-03-07T17:06:26.609669Z",
     "shell.execute_reply": "2021-03-07T17:06:26.610190Z"
    }
   },
   "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\n",
    "from outputC import outputC,add_to_Cfunction_dict  # NRPy+: Core C code output module\n",
    "import sympy as sp                  # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import NRPy_param_funcs as par      # NRPy+: Parameter interface\n",
    "import indexedexp as ixp            # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm      # NRPy+: Reference metric 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": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:06:26.774454Z",
     "iopub.status.busy": "2021-03-07T17:06:26.639256Z",
     "iopub.status.idle": "2021-03-07T17:06:26.791454Z",
     "shell.execute_reply": "2021-03-07T17:06:26.791961Z"
    }
   },
   "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-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 ):\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)\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": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:06:26.937141Z",
     "iopub.status.busy": "2021-03-07T17:06:26.797807Z",
     "iopub.status.idle": "2021-03-07T17:06:26.954030Z",
     "shell.execute_reply": "2021-03-07T17:06:26.954518Z"
    },
    "scrolled": true
   },
   "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='c_code_generation'></a>\n",
    "\n",
    "# Step 4: C-code routines (library) for reading TOV data files \\[Back to [top](#toc)\\]\n",
    "$$\\label{c_code_generation}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='id_persist'></a>\n",
    "\n",
    "## Step 4.a: Declare `ID_persist`, the C `struct` data type that will store data read from file \\[Back to [top](#toc)\\]\n",
    "$$\\label{id_persist}$$\n",
    "\n",
    "The TOV data file stores various physical quantities describing the star's pressure, density, and gravitational fields, *all as a function of distance from the center of the star*, `rbar`. `Rbar` is the radius of the star.\n",
    "\n",
    "The `ID_persist` struct takes the following form\n",
    "\n",
    "```c\n",
    "  REAL Rbar;      // rbar corresponding to the outermost radius at which density is nonzero\n",
    "  int Rbar_idx;   // Index (line of data file) corresponding to Rbar\n",
    "  int interp_stencil_size;  // Lagrange polynomial interpolation stencil size\n",
    "  int numlines_in_file;     // Number of radii stored in the file\n",
    "  // The following arrays store stellar information at all numlines_in_file radii:\n",
    "  REAL *restrict r_Schw_arr; // Stellar radial coordinate in units of Schwarzschild radius\n",
    "  REAL *restrict rho_arr;    // Mass-energy density\n",
    "  REAL *restrict rho_baryon_arr; // Baryonic mass density\n",
    "  REAL *restrict P_arr;  // Pressure\n",
    "  REAL *restrict M_arr;  // Integrated rest mass\n",
    "  REAL *restrict expnu_arr;    // Metric quantity\n",
    "  REAL *restrict exp4phi_arr;  // Metric quantity\n",
    "  REAL *restrict rbar_arr;  // Rbar coordinate\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ID_persist_str():\n",
    "    return r\"\"\"\n",
    "  REAL Rbar;      // rbar corresponding to the outermost radius at which density is nonzero\n",
    "  int Rbar_idx;   // Index (line of data file) corresponding to Rbar\n",
    "  int interp_stencil_size;  // Lagrange polynomial interpolation stencil size\n",
    "  int numlines_in_file;     // Number of radii stored in the file\n",
    "  // The following arrays store stellar information at all numlines_in_file radii:\n",
    "  REAL *restrict r_Schw_arr; // Stellar radial coordinate in units of Schwarzschild radius\n",
    "  REAL *restrict rho_arr;    // Mass-energy density\n",
    "  REAL *restrict rho_baryon_arr; // Baryonic mass density\n",
    "  REAL *restrict P_arr;  // Pressure\n",
    "  REAL *restrict M_arr;  // Integrated rest mass\n",
    "  REAL *restrict expnu_arr;    // Metric quantity\n",
    "  REAL *restrict exp4phi_arr;  // Metric quantity\n",
    "  REAL *restrict rbar_arr;  // Rbar coordinate\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='read_data_file'></a>\n",
    "\n",
    "## Step 4.b: `TOV_read_data_file_set_ID_persist()`: Read TOV data file and store data to the `ID_persist` struct \\[Back to [top](#toc)\\]\n",
    "$$\\label{read_data_file}$$\n",
    "\n",
    "The following function also sets the `interp_stencil_size` parameter, indicating the total size of the Lagrange polynomial interpolation stencil. The default of 12 reflects a quite high interpolation order, corresponding to an 11th-order polynomial being fit through the data at various radii to estimate stellar quantities at a single desired (arbitrary) distance from the center of the star.\n",
    "\n",
    "As hydrodynamic data (like density and pressure) at the stellar radius sharply drop to zero, the interpolation algorithm offsets the center of the stencil so that the interpolation never crosses the stellar surface. This avoids the [Gibbs phenomenon](https://en.wikipedia.org/wiki/Gibbs_phenomenon), ensuring super high fidelity of the output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_TOV_read_data_file_set_ID_persist(interp_stencil_size=12):\n",
    "    includes = [\"NRPy_basic_defines.h\"]\n",
    "    desc = \"Returns the number of lines in a TOV data file.\"\n",
    "    c_type = \"void\"\n",
    "    name = \"TOV_read_data_file_set_ID_persist\"\n",
    "    params = \"const char *input_filename, ID_persist_struct *ID_persist\"\n",
    "    body = r\"\"\"\n",
    "  char filename[100];\n",
    "  snprintf(filename, 100, input_filename);\n",
    "  FILE *TOV_solution_datafile = fopen(filename, \"r\");\n",
    "  if(TOV_solution_datafile == NULL) {\n",
    "    fprintf(stderr,\"ERROR: could not open TOV solution data file %s\\n\",filename);\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  // First set the interpolation stencil size. Note that the interpolation is set up to\n",
    "  //   avoid interpolations across the star's surface. Hence we can use a super high\n",
    "  //   order interpolant.\n",
    "  ID_persist->interp_stencil_size = \"\"\"+str(interp_stencil_size)+r\"\"\";\n",
    "\n",
    "  int numlines_in_file = 0;\n",
    "  {\n",
    "    char * line = NULL;\n",
    "\n",
    "    size_t len = 0;\n",
    "    ssize_t read;\n",
    "    while ((read = getline(&line, &len, TOV_solution_datafile)) != -1) {\n",
    "      numlines_in_file++;\n",
    "    }\n",
    "    rewind(TOV_solution_datafile);\n",
    "\n",
    "    free(line);\n",
    "  }\n",
    "  ID_persist->numlines_in_file = numlines_in_file;\n",
    "\n",
    "  // Now that numlines_in_file is set, we can now allocate memory for all arrays.\n",
    "  {\n",
    "    ID_persist->r_Schw_arr     = (REAL *restrict)malloc(sizeof(REAL)*numlines_in_file);\n",
    "    ID_persist->rho_arr        = (REAL *restrict)malloc(sizeof(REAL)*numlines_in_file);\n",
    "    ID_persist->rho_baryon_arr = (REAL *restrict)malloc(sizeof(REAL)*numlines_in_file);\n",
    "    ID_persist->P_arr          = (REAL *restrict)malloc(sizeof(REAL)*numlines_in_file);\n",
    "    ID_persist->M_arr          = (REAL *restrict)malloc(sizeof(REAL)*numlines_in_file);\n",
    "    ID_persist->expnu_arr      = (REAL *restrict)malloc(sizeof(REAL)*numlines_in_file);\n",
    "    ID_persist->exp4phi_arr    = (REAL *restrict)malloc(sizeof(REAL)*numlines_in_file);\n",
    "    ID_persist->rbar_arr       = (REAL *restrict)malloc(sizeof(REAL)*numlines_in_file);\n",
    "  }\n",
    "\n",
    "  {\n",
    "    char * line = NULL;\n",
    "\n",
    "    size_t len = 0;\n",
    "    ssize_t read;\n",
    "\n",
    "    int which_line = 0;\n",
    "    while ((read = getline(&line, &len, TOV_solution_datafile)) != -1) {\n",
    "      // Define the line delimiters (i.e., the stuff that goes between the data on a given\n",
    "      //     line of data.  Here, we define both spaces \" \" and tabs \"\\t\" as data delimiters.\n",
    "      const char delimiters[] = \" \\t\";\n",
    "\n",
    "      // Now we define \"token\", a pointer to the first column of data\n",
    "      char *token;\n",
    "\n",
    "      // Each successive time we call strtok(NULL,blah), we read in a new column of data from\n",
    "      //     the originally defined character array, as pointed to by token.\n",
    "\n",
    "      token=strtok(line, delimiters); if(token==NULL) { fprintf(stderr, \"Error reading %s\\n\", filename); exit(1); }\n",
    "      ID_persist->r_Schw_arr[which_line]     = strtod(token, NULL); token = strtok( NULL, delimiters );\n",
    "      ID_persist->rho_arr[which_line]        = strtod(token, NULL); token = strtok( NULL, delimiters );\n",
    "      ID_persist->rho_baryon_arr[which_line] = strtod(token, NULL); token = strtok( NULL, delimiters );\n",
    "      ID_persist->P_arr[which_line]          = strtod(token, NULL); token = strtok( NULL, delimiters );\n",
    "      ID_persist->M_arr[which_line]          = strtod(token, NULL); token = strtok( NULL, delimiters );\n",
    "      ID_persist->expnu_arr[which_line]      = strtod(token, NULL); token = strtok( NULL, delimiters );\n",
    "      ID_persist->exp4phi_arr[which_line]    = strtod(token, NULL); token = strtok( NULL, delimiters );\n",
    "      ID_persist->rbar_arr[which_line]       = strtod(token, NULL);\n",
    "\n",
    "      which_line++;\n",
    "    }\n",
    "    free(line);\n",
    "\n",
    "    fclose(TOV_solution_datafile);\n",
    "  }\n",
    "\n",
    "  {\n",
    "    // Finally set Rbar and Rbar_idx\n",
    "    ID_persist->Rbar     = -100.0;\n",
    "    ID_persist->Rbar_idx = -100;\n",
    "    for(int i=1;i<numlines_in_file;i++) {\n",
    "      if(ID_persist->rho_arr[i-1] > 0  &&  ID_persist->rho_arr[i] == 0) {\n",
    "        ID_persist->Rbar = ID_persist->rbar_arr[i-1];\n",
    "        ID_persist->Rbar_idx = i-1;\n",
    "      }\n",
    "    }\n",
    "    if(ID_persist->Rbar < 0) {\n",
    "      fprintf(stderr,\"Error: could not find rbar=Rbar (i.e., the surface of the star) from data file.\\n\");\n",
    "      exit(1);\n",
    "    }\n",
    "  }\n",
    "\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        enableCparameters=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='interp_data_file'></a>\n",
    "\n",
    "## Step 4.c: `TOV_interpolate_1D()`: Interpolate TOV data to any desired distance from the center of the star \\[Back to [top](#toc)\\]\n",
    "$$\\label{interp_data_file}$$\n",
    "\n",
    "The TOV solution stored in the data file samples only a very limited number of points inside and outside the star. However we'll need data at arbitrary radii between these points. `TOV_interpolate_1D()` provides this service.\n",
    "\n",
    "`TOV_interpolate_1D()` applies 1-dimensional Lagrange polynomial interpolation with stencil size `ID_persist->interp_stencil_size` (default 12) to obtain TOV data between sampled points in the ordered data file (each point corresponds to a specific radius).\n",
    "\n",
    "Once a desired output radius `rrbar` is chosen, `TOV_interpolate_1D()` calls the included bisection algorithm `bisection_idx_finder()` to find the closest sample point with radius `rbar`, to desired output point `rrbar`. If the stencil crosses the star's surface, it then adjusts the stencil to avoid the [Gibbs phenomenon](https://en.wikipedia.org/wiki/Gibbs_phenomenon) associated with Lagrange interpolating data with a kink."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_TOV_interpolate_1D():\n",
    "    includes=[\"NRPy_basic_defines.h\"]\n",
    "    prefunc = r\"\"\"\n",
    "// Find interpolation index using Bisection root-finding algorithm:\n",
    "static inline int bisection_idx_finder(const REAL rrbar, const int numlines_in_file, const REAL *restrict rbar_arr) {\n",
    "  int x1 = 0;\n",
    "  int x2 = numlines_in_file-1;\n",
    "  REAL y1 = rrbar-rbar_arr[x1];\n",
    "  REAL y2 = rrbar-rbar_arr[x2];\n",
    "  if(y1*y2 >= 0) {\n",
    "    fprintf(stderr,\"INTERPOLATION BRACKETING ERROR %e | %e %e\\n\",rrbar,y1,y2);\n",
    "    exit(1);\n",
    "  }\n",
    "  for(int i=0;i<numlines_in_file;i++) {\n",
    "    int x_midpoint = (x1+x2)/2;\n",
    "    REAL y_midpoint = rrbar-rbar_arr[x_midpoint];\n",
    "    if(y_midpoint*y1 < 0) {\n",
    "      x2 = x_midpoint;\n",
    "      y2 = y_midpoint;\n",
    "    } else {\n",
    "      x1 = x_midpoint;\n",
    "      y1 = y_midpoint;\n",
    "    }\n",
    "    if( abs(x2-x1) == 1 ) {\n",
    "      // If rbar_arr[x1] is closer to rrbar than rbar_arr[x2] then return x1:\n",
    "      if(fabs(rrbar-rbar_arr[x1]) < fabs(rrbar-rbar_arr[x2])) return x1;\n",
    "      // Otherwiser return x2:\n",
    "      return x2;\n",
    "    }\n",
    "  }\n",
    "  fprintf(stderr,\"INTERPOLATION BRACKETING ERROR: DID NOT CONVERGE.\\n\");\n",
    "  exit(1);\n",
    "}\n",
    "\"\"\"\n",
    "    desc = \"\"\"Read a TOV solution from data file and perform\n",
    "1D interpolation of the solution to a desired radius.\n",
    "\n",
    "Author: Zachariah B. Etienne\n",
    "        zachetie **at** gmail **dot* com\n",
    "\"\"\"\n",
    "    c_type=\"void\"\n",
    "    name=\"TOV_interpolate_1D\"\n",
    "    params=\"\"\"REAL rrbar,const ID_persist_struct *ID_persist,\n",
    "                        REAL *restrict rho,REAL *restrict rho_baryon,REAL *restrict P,REAL *restrict M,REAL *restrict expnu,REAL *restrict exp4phi\"\"\"\n",
    "    body = r\"\"\"\n",
    "  // First unpack ID_persist struct (this should be pretty quick relative to the rest of the routine):\n",
    "  const REAL Rbar               = ID_persist->Rbar;\n",
    "  const int Rbar_idx            = ID_persist->Rbar_idx;\n",
    "  const int interp_stencil_size = ID_persist->interp_stencil_size;\n",
    "  const int numlines_in_file    = ID_persist->numlines_in_file;\n",
    "  const REAL *restrict rbar_arr       = ID_persist->rbar_arr;\n",
    "  const REAL *restrict r_Schw_arr     = ID_persist->r_Schw_arr;\n",
    "  const REAL *restrict rho_arr        = ID_persist->rho_arr;\n",
    "  const REAL *restrict rho_baryon_arr = ID_persist->rho_baryon_arr;\n",
    "  const REAL *restrict P_arr          = ID_persist->P_arr;\n",
    "  const REAL *restrict M_arr          = ID_persist->M_arr;\n",
    "  const REAL *restrict expnu_arr      = ID_persist->expnu_arr;\n",
    "  const REAL *restrict exp4phi_arr    = ID_persist->exp4phi_arr;\n",
    "\n",
    "  // For this case, we know that for all functions, f(r) = f(-r)\n",
    "  if(rrbar < 0) rrbar = -rrbar;\n",
    "\n",
    "  // First find the central interpolation stencil index:\n",
    "  int idx = bisection_idx_finder(rrbar,numlines_in_file,rbar_arr);\n",
    "\n",
    "\n",
    "  int idxmin = MAX(0,idx-interp_stencil_size/2-1);\n",
    "\n",
    "  // -= Do not allow the interpolation stencil to cross the star's surface =-\n",
    "  // max index is when idxmin + (interp_stencil_size-1) = Rbar_idx\n",
    "  //  -> idxmin at most can be Rbar_idx - interp_stencil_size + 1\n",
    "  if(rrbar < Rbar) {\n",
    "    idxmin = MIN(idxmin,Rbar_idx - interp_stencil_size + 1);\n",
    "  } else {\n",
    "    idxmin = MAX(idxmin,Rbar_idx+1);\n",
    "    idxmin = MIN(idxmin,numlines_in_file - interp_stencil_size + 1);\n",
    "  }\n",
    "  // Now perform the Lagrange polynomial interpolation:\n",
    "\n",
    "  // First set the interpolation coefficients:\n",
    "  REAL rbar_sample[interp_stencil_size];\n",
    "  for(int i=idxmin;i<idxmin+interp_stencil_size;i++) {\n",
    "    rbar_sample[i-idxmin] = rbar_arr[i];\n",
    "  }\n",
    "  REAL l_i_of_r[interp_stencil_size];\n",
    "  for(int i=0;i<interp_stencil_size;i++) {\n",
    "    REAL numer = 1.0;\n",
    "    REAL denom = 1.0;\n",
    "    for(int j=0;j<i;j++) {\n",
    "      numer *= rrbar - rbar_sample[j];\n",
    "      denom *= rbar_sample[i] - rbar_sample[j];\n",
    "    }\n",
    "    for(int j=i+1;j<interp_stencil_size;j++) {\n",
    "      numer *= rrbar - rbar_sample[j];\n",
    "      denom *= rbar_sample[i] - rbar_sample[j];\n",
    "    }\n",
    "    l_i_of_r[i] = numer/denom;\n",
    "  }\n",
    "\n",
    "  // Then perform the interpolation:\n",
    "  *rho = 0.0;\n",
    "  *rho_baryon = 0.0;\n",
    "  *P = 0.0;\n",
    "  *M = 0.0;\n",
    "  *expnu = 0.0;\n",
    "  *exp4phi = 0.0;\n",
    "\n",
    "  REAL r_Schw = 0.0;\n",
    "  for(int i=idxmin;i<idxmin+interp_stencil_size;i++) {\n",
    "    r_Schw      += l_i_of_r[i-idxmin] * r_Schw_arr[i];\n",
    "    *rho        += l_i_of_r[i-idxmin] * rho_arr[i];\n",
    "    *rho_baryon += l_i_of_r[i-idxmin] * rho_baryon_arr[i];\n",
    "    *P          += l_i_of_r[i-idxmin] * P_arr[i];\n",
    "    *M          += l_i_of_r[i-idxmin] * M_arr[i];\n",
    "    *expnu      += l_i_of_r[i-idxmin] * expnu_arr[i];\n",
    "    *exp4phi    += l_i_of_r[i-idxmin] * exp4phi_arr[i];\n",
    "  }\n",
    "\n",
    "  if(rrbar > Rbar) {\n",
    "    *rho        = 0;\n",
    "    *rho_baryon = 0;\n",
    "    *P          = 0;\n",
    "    *M          = M_arr[Rbar_idx+1];\n",
    "    *expnu      = 1. - 2.*(*M) / r_Schw;\n",
    "    *exp4phi    = pow(r_Schw / rrbar,2.0);\n",
    "  }\n",
    "\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        prefunc=prefunc,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        enableCparameters=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='tov_id_function'></a>\n",
    "\n",
    "## Step 4.d: `TOV_ID_function()`: Driver routine for setting up ADM metric quantities and $T^{\\mu\\nu}$ from interpolated TOV data \\[Back to [top](#toc)\\]\n",
    "$$\\label{tov_id_function}$$\n",
    "\n",
    "`TOV_ID_function()` is a function with standard input and output parameters, providing the necessary linkage to the NRPy+ BSSN initial data infrastructure. This infrastructure takes as input ADM metric data and $T^{\\mu\\nu}$ in the spherical or Cartesian basis. TOV data, of course, are provided in the spherical basis.\n",
    "\n",
    "In other words, `TOV_ID_function()` takes as input the desired Cartesian point $(x,y,z)$ at which we desire initial data, and calls the `TOV_interpolate_1D()` routine to get raw data at that point. `TOV_ID_function()` converts the interpolated quantities to ADM metric and $T^{\\mu\\nu}$, stored in the standardized `initial_data` struct.\n",
    "\n",
    "Conversions between data within the TOV solution output file and ADM/$T^{\\mu\\nu}$ are as follows:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='adm_ito_tov'></a>\n",
    "\n",
    "### Step 4.d.i: ADM metric quantities, written in terms of the TOV solution \\[Back to [top](#toc)\\]\n",
    "$$\\label{adm_ito_tov}$$\n",
    "\n",
    "The [TOV line element](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation) in *Schwarzschild coordinates* 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",
    "\n",
    "In *isotropic coordinates* with $G=c=1$ (i.e., the coordinate system 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",
    "The ADM 3+1 line element for this diagonal metric in isotropic spherical coordinates is given by:\n",
    "$$\n",
    "ds^2 = (-\\alpha^2 + \\beta_k \\beta^k) dt^2 + \\gamma_{\\bar{r}\\bar{r}} d\\bar{r}^2 + \\gamma_{\\theta\\theta} d\\theta^2+ \\gamma_{\\phi\\phi} d\\phi^2,\n",
    "$$\n",
    "\n",
    "from which we can immediately read off the ADM quantities:\n",
    "\\begin{align}\n",
    "\\alpha &= e^{\\nu(\\bar{r})/2} \\\\\n",
    "\\beta^k &= 0 \\\\\n",
    "\\gamma_{\\bar{r}\\bar{r}} &= e^{4\\phi}\\\\\n",
    "\\gamma_{\\theta\\theta} &= e^{4\\phi} \\bar{r}^2 \\\\\n",
    "\\gamma_{\\phi\\phi} &= e^{4\\phi} \\bar{r}^2 \\sin^2 \\theta \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ADM_quantities_ito_TOV_soln(rbar, theta, expnu, exp4phi):\n",
    "    # in TOV ID, betaU=BU=KDD=0\n",
    "    alpha = sp.sqrt(expnu)\n",
    "\n",
    "    betaU = ixp.zerorank1()\n",
    "    BU = ixp.zerorank1()\n",
    "\n",
    "    gammaDD = ixp.zerorank2(DIM=3)\n",
    "    gammaDD[0][0] = exp4phi\n",
    "    gammaDD[1][1] = exp4phi * rbar ** 2\n",
    "    gammaDD[2][2] = exp4phi * rbar ** 2 * sp.sin(theta) ** 2\n",
    "\n",
    "    KDD = ixp.zerorank2()\n",
    "\n",
    "    return alpha, betaU, BU, gammaDD, KDD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='tmunu_ito_tov'></a>\n",
    "\n",
    "### Step 4.d.ii: Stress-energy tensor $T^{\\mu\\nu}$, written in terms of the TOV solution \\[Back to [top](#toc)\\]\n",
    "$$\\label{tmunu_ito_tov}$$\n",
    "\n",
    "\n",
    "We will also need the stress-energy tensor $T^{\\mu\\nu}$. [As discussed here](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation), the stress-energy tensor is diagonal:\n",
    "\n",
    "\\begin{align}\n",
    "T^t_t &= -\\rho \\\\\n",
    "T^i_j &= P \\delta^i_j \\\\\n",
    "\\text{All other components of }T^\\mu_\\nu &= 0.\n",
    "\\end{align}\n",
    "\n",
    "Since $\\beta^i=0$ the inverse metric expression simplifies to (Eq. 4.49 in [Gourgoulhon](https://arxiv.org/pdf/gr-qc/0703035.pdf)):\n",
    "$$\n",
    "g^{\\mu\\nu} = \\begin{pmatrix} \n",
    "-\\frac{1}{\\alpha^2} & \\frac{\\beta^i}{\\alpha^2} \\\\\n",
    "\\frac{\\beta^i}{\\alpha^2} & \\gamma^{ij} - \\frac{\\beta^i\\beta^j}{\\alpha^2}\n",
    "\\end{pmatrix} =\n",
    "\\begin{pmatrix} \n",
    "-\\frac{1}{\\alpha^2} & 0 \\\\\n",
    "0 & \\gamma^{ij}\n",
    "\\end{pmatrix},\n",
    "$$\n",
    "\n",
    "and since the 3-metric is diagonal we get\n",
    "\n",
    "\\begin{align}\n",
    "\\gamma^{\\bar{r}\\bar{r}} &= e^{-4\\phi}\\\\\n",
    "\\gamma^{\\theta\\theta} &= e^{-4\\phi}\\frac{1}{\\bar{r}^2} \\\\\n",
    "\\gamma^{\\phi\\phi} &= e^{-4\\phi}\\frac{1}{\\bar{r}^2 \\sin^2 \\theta}.\n",
    "\\end{align}\n",
    "\n",
    "Thus raising $T^\\mu_\\nu$ yields a diagonal $T^{\\mu\\nu}$\n",
    "\n",
    "\\begin{align}\n",
    "T^{tt} &= -g^{tt} \\rho = \\frac{1}{\\alpha^2} \\rho = e^{-\\nu(\\bar{r})} \\rho \\\\\n",
    "T^{\\bar{r}\\bar{r}} &= g^{\\bar{r}\\bar{r}} P = \\frac{1}{e^{4 \\phi}} P \\\\\n",
    "T^{\\theta\\theta} &= g^{\\theta\\theta} P = \\frac{1}{e^{4 \\phi}\\bar{r}^2} P\\\\\n",
    "T^{\\phi\\phi} &= g^{\\phi\\phi} P = \\frac{1}{e^{4\\phi}\\bar{r}^2 \\sin^2 \\theta} P \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def T4UU_ito_TOV_soln(rbar, theta, rho, P, expnu, exp4phi):\n",
    "\n",
    "    T4UU = ixp.zerorank2(DIM=4)\n",
    "    T4UU[0][0] = rho/expnu\n",
    "    T4UU[1][1] = P/exp4phi\n",
    "    T4UU[2][2] = P/(exp4phi*rbar**2)\n",
    "    T4UU[3][3] = P/(exp4phi*rbar**2*sp.sin(theta)**2)\n",
    "\n",
    "    return T4UU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='tov_id_function_itself'></a>\n",
    "\n",
    "### Step 4.d.iii: The `TOV_ID_function()` itself \\[Back to [top](#toc)\\]\n",
    "$$\\label{tov_id_function_itself}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_TOV_ID_function():\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    desc = \"\"\"This function sets TOV initial data at a single point xCart[3] = (x,y,z),\n",
    "outputting to initial_data_struct *restrict initial_data\"\"\"\n",
    "    c_type = \"void\"\n",
    "    name = \"TOV_ID_function\"\n",
    "    params = \"const paramstruct *params, const REAL xCart[3], const ID_persist_struct *restrict ID_persist, initial_data_struct *restrict initial_data\"\n",
    "    body  = \"\"\"  // First set r(=rbar), theta, phi in terms of xCart[3]:\n",
    "  const REAL Cartx=xCart[0], Carty=xCart[1], Cartz=xCart[2];\"\"\"\n",
    "    desired_rfm_coord = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
    "    # \"Spherical\" is the native basis for TOV initial data.\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\", \"Spherical\")\n",
    "    rfm.reference_metric()\n",
    "    body += r\"\"\"\n",
    "  REAL rbar, theta, phi;\n",
    "  {\n",
    "\"\"\" + outputC(rfm.Cart_to_xx[:3], [\"rbar\", \"theta\", \"phi\"], filename=\"returnstring\",\n",
    "              params=\"outCverbose=False,includebraces=False,preindent=2\")\n",
    "    body += r\"\"\"\n",
    "  }\n",
    "\n",
    "  // Next set gamma_{ij} in spherical basis\n",
    "  REAL rho,rho_baryon,P,M,expnu,exp4phi;\n",
    "  TOV_interpolate_1D(rbar,ID_persist, &rho,&rho_baryon,&P,&M,&expnu,&exp4phi);\n",
    "\"\"\"\n",
    "    rbar, theta, rho, P, expnu, exp4phi = sp.symbols('rbar theta rho P expnu exp4phi', real=True)\n",
    "    alpha, betaU, BU, gammaDD, KDD = ADM_quantities_ito_TOV_soln(rbar, theta, expnu, exp4phi)\n",
    "    T4UU = T4UU_ito_TOV_soln(rbar, theta, rho, P, expnu, exp4phi)\n",
    "\n",
    "    list_of_output_exprs = [alpha]\n",
    "    list_of_output_varnames = [\"initial_data->alpha\"]\n",
    "    for i in range(3):\n",
    "        list_of_output_exprs += [betaU[i]]\n",
    "        list_of_output_varnames += [\"initial_data->betaSphorCartU\" + str(i)]\n",
    "        list_of_output_exprs += [BU[i]]\n",
    "        list_of_output_varnames += [\"initial_data->BSphorCartU\" + str(i)]\n",
    "        for j in range(i, 3):\n",
    "            list_of_output_exprs += [gammaDD[i][j]]\n",
    "            list_of_output_varnames += [\"initial_data->gammaSphorCartDD\" + str(i) + str(j)]\n",
    "            list_of_output_exprs += [KDD[i][j]]\n",
    "            list_of_output_varnames += [\"initial_data->KSphorCartDD\" + str(i) + str(j)]\n",
    "    for mu in range(4):\n",
    "        for nu in range(mu, 4):\n",
    "            list_of_output_exprs += [T4UU[mu][nu]]\n",
    "            list_of_output_varnames += [\"initial_data->T4SphorCartUU\" + str(mu) + str(nu)]\n",
    "\n",
    "    # Sort the outputs before calling outputC()\n",
    "    # https://stackoverflow.com/questions/9764298/is-it-possible-to-sort-two-listswhich-reference-each-other-in-the-exact-same-w\n",
    "    list_of_output_varnames, list_of_output_exprs = (list(t) for t in zip(*sorted(zip(list_of_output_varnames, list_of_output_exprs))))\n",
    "\n",
    "    body += outputC(list_of_output_exprs, list_of_output_varnames,\n",
    "                    filename=\"returnstring\", params=\"outCverbose=False,includebraces=False,preindent=1\")\n",
    "\n",
    "    # Restore CoordSystem:\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\", desired_rfm_coord)\n",
    "    rfm.reference_metric()\n",
    "\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        enableCparameters=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ccodegen_library_validation'></a>\n",
    "\n",
    "# Step 5: Validate C code generation against `TOV.TOV_Ccodegen_library` Python module \\[Back to [top](#toc)\\]\n",
    "$$\\label{ccodegen_library_validation}$$\n",
    "\n",
    "To ensure this documentation remains maintained in lockstep with development of the separate `TOV.TOV_Ccodegen_library` Python module, we register all the C functions both here and from the Python module and compare them. If they are found to be *identical*, the test passes. Otherwise it fails and the notebook will exit with an error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PASS! ALL FUNCTIONS ARE IDENTICAL\n"
     ]
    }
   ],
   "source": [
    "import TOV.TOV_Ccodegen_library as TCL\n",
    "import sys\n",
    "\n",
    "funclist = [(\"ID_persist_str\", ID_persist_str, TCL.ID_persist_str),\n",
    "            (\"add_to_Cfunction_dict_TOV_read_data_file_set_ID_persist\", add_to_Cfunction_dict_TOV_read_data_file_set_ID_persist, TCL.add_to_Cfunction_dict_TOV_read_data_file_set_ID_persist),\n",
    "            (\"add_to_Cfunction_dict_TOV_interpolate_1D\", add_to_Cfunction_dict_TOV_interpolate_1D, TCL.add_to_Cfunction_dict_TOV_interpolate_1D),\n",
    "            (\"ADM_quantities_ito_TOV_soln\", ADM_quantities_ito_TOV_soln, TCL.ADM_quantities_ito_TOV_soln),\n",
    "            (\"T4UU_ito_TOV_soln\", T4UU_ito_TOV_soln, TCL.T4UU_ito_TOV_soln),\n",
    "            (\"ADM_quantities_ito_TOV_soln\", ADM_quantities_ito_TOV_soln, TCL.ADM_quantities_ito_TOV_soln),\n",
    "            (\"add_to_Cfunction_dict_TOV_ID_function\", add_to_Cfunction_dict_TOV_ID_function, TCL.add_to_Cfunction_dict_TOV_ID_function)\n",
    "           ]\n",
    "\n",
    "if sys.version_info.major >= 3:\n",
    "    import inspect\n",
    "    for func in funclist:\n",
    "        # https://stackoverflow.com/questions/20059011/check-if-two-python-functions-are-equal\n",
    "        if inspect.getsource(func[1]) != inspect.getsource(func[2]):\n",
    "            print(\"inspect.getsource(func[1]):\")\n",
    "            print(inspect.getsource(func[1]))\n",
    "            print(\"inspect.getsource(func[2]):\")\n",
    "            print(inspect.getsource(func[2]))\n",
    "            print(\"ERROR: function \" + func[0] + \" is not the same as the Ccodegen_library version!\")\n",
    "            sys.exit(1)\n",
    "\n",
    "    print(\"PASS! ALL FUNCTIONS ARE IDENTICAL\")\n",
    "else:\n",
    "    print(\"SORRY CANNOT CHECK FUNCTION IDENTITY WITH PYTHON 2. PLEASE update your Python installation.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 6: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[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": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:06:26.960083Z",
     "iopub.status.busy": "2021-03-07T17:06:26.958715Z",
     "iopub.status.idle": "2021-03-07T17:06:30.249332Z",
     "shell.execute_reply": "2021-03-07T17:06:30.250235Z"
    }
   },
   "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 (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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}