Tutorial-GRHD_UnitConversion.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",
    "# Code and Physical Units\n",
    "\n",
    "## GRHD Units in terms of EOS \n",
    "\n",
    "$\\newcommand{\\rhoCode}{{\\tilde{\\rho}}}$\n",
    "\n",
    "$\\newcommand{\\MCode}{{\\tilde{M}}}$ $\\newcommand{\\rCode}{{\\tilde{r}}}$ $\\newcommand{\\PCode}{{\\tilde{P}}}$$\\newcommand{\\tCode}{{\\tilde{t}}}$$\\newcommand{\\Mfid}{{M_{\\rm fid}}}$$\\newcommand{\\MfidBar}{\\bar{M}_{\\rm fid}}$$\\newcommand{\\Mbar}{\\bar{M}}$\n",
    "$\\newcommand{\\rBar}{\\bar{r}}$$\\newcommand{\\tBar}{\\bar{t}}$\n",
    "In GRHD, we can set an equation of state of the form\n",
    "\\begin{equation}\n",
    "P = K\\rho^{1+1/n}\n",
    "\\end{equation}\n",
    "Taking $c_s^2 = \\partial P/\\partial \\rho = (1+1/n) K\\rho^{1/n}$.  This gives for some fidicial $\\rho_0$\n",
    "\\begin{equation}\n",
    "c_{s,0}^2 = \\left(1 + \\frac 1 n\\right)K\\rho_0^{1/n}.\n",
    "\\end{equation}\n",
    "Selecting $c_s^2 = c^2\\left(1 + 1/n\\right)$, we have\n",
    "\\begin{equation}\n",
    "\\rho_0 = \\left(\\frac {c^2}{K}\\right)^n\n",
    "\\end{equation}\n",
    "This is equivalent to setting the isothermal sound speed to $c$.  With this definition of $\\rho_0$, we can write\n",
    "\\begin{equation}\n",
    "P = \\rho_0c^2\\left(\\frac{\\rho}{\\rho_0}\\right)^{1+1/n}\n",
    "\\end{equation}\n",
    "which allows us to define the dimensionless density $\\rhoCode = \\rho/\\rho_0$ and dimensionless pressure $\\PCode = P/\\rho_0 c^2$\n",
    "\\begin{equation}\n",
    "\\PCode = \\rhoCode^{1+1/n},\n",
    "\\end{equation}\n",
    "where we adopt code units where $c=1$.  These dimensionless pressure and density are in $G=c=1$ units and can be used in GRHD code including inclusion in the spacetime solver via $T_{\\mu\\nu}$.  Note that this sets $K=1$ in these units.\n",
    "\n",
    "To find a dimensionless mass, $\\MCode$, dimensionless distance, $\\rCode$, and dimensionless time, $\\tCode$, we note\n",
    "$GM/rc^2$ is dimensionless\n",
    "\\begin{equation}\n",
    "\\frac{GM}{rc^2} = \\frac{G\\rho_0 r^2}{c^2} = \\frac{Gc^{2n-2}}{K^n}r^2 \\rightarrow \\rCode = \\frac{\\sqrt{G}c^{n-1}}{K^{n/2}} r = \\frac r {r_0},\n",
    "\\end{equation}\n",
    "where $r_0 = K^{n/2}/\\sqrt{G}c^{n-1}$.  Then\n",
    "\\begin{eqnarray}\n",
    "\\tCode &=& \\frac{t}{t_0} = \\frac{t}{r_0/c} = \\frac{\\sqrt{G}c^n}{K^{n/2}} t \\\\\n",
    "\\MCode &=& \\frac{M}{M_0} = \\frac{M}{\\rho_0 r_0^3} = M\\frac{K^n}{c^{2n}}\\frac{G^{3/2}c^{3(n-1)}}{K^{3n/2}} = \\frac{G^{3/2}c^{n-3}}{K^{n/2}} M,\n",
    "\\end{eqnarray}\n",
    "Hence, we have \n",
    "\\begin{eqnarray}\n",
    "\\rho_0 &=& \\left(\\frac{K}{c^2}\\right)^n\\\\\n",
    "r_0 &=& \\frac{c^{n+1}}{\\sqrt{G}K^{n/2}}\\\\\n",
    "t_0 &=& \\frac{c^{n}}{\\sqrt{G}K^{n/2}}\\\\\n",
    "M_0 &=& \\frac{c^{n+3}}{G^{3/2}K^{n/2}}\n",
    "\\end{eqnarray}\n",
    "\n",
    "## Mapping to SENR or any NR code\n",
    "\n",
    "So we will need a $\\Mfid$ which is define such that the (SENR) code units $\\MfidBar = 1$ or in other words in SENR codes units: \n",
    "\\begin{equation}\n",
    "\\Mbar = \\frac{M}{\\Mfid}\n",
    "\\end{equation}\n",
    "In these units:\n",
    "\\begin{eqnarray}\n",
    "\\rBar &=& \\frac{c^2}{G\\Mfid} r\\\\\n",
    "\\tBar &=& \\frac{c^3}{G\\Mfid} t\n",
    "\\end{eqnarray}\n",
    "At some level $\\Mfid$ is arbitrary, so we can select $M_0 = \\Mfid$.  In this case, this means that $\\rBar = \\rCode$, $\\tBar = \\tCode$, and $\\Mbar = \\MCode$, which fixes all the quantities. This comes at a cost the $\\bar{M}_{\\rm ADM}$ is not something nice like 1 or 2, but the choice is consistent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zach's comments:\n",
    "\n",
    "Sound speed $c_s$ is defined as\n",
    "\n",
    "$$\\frac{\\partial P}{\\partial \\rho} = c_s^2,$$\n",
    "\n",
    "so if we have a polytropic EOS, where\n",
    "\n",
    "$$P = K \\rho^{(1 + 1/n)},$$\n",
    "\n",
    "then\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial P}{\\partial \\rho} &= c_s^2 \\\\\n",
    "&= \\left(1 + \\frac{1}{n}\\right) K \\rho^{1/n}.\n",
    "\\end{align}\n",
    "\n",
    "Let's adopt the notation \n",
    "\n",
    "$$[\\rho] = \\text{\"the units of $\\rho$\"}$$\n",
    "\n",
    "Using this notation and the fact that $n$ is dimensionless, the above expression implies\n",
    "\n",
    "\\begin{align}\n",
    "\\left[\\rho^{1/n}\\right] &= \\left[\\frac{c_s^2}{K}\\right] \\\\\n",
    "\\implies \\left[\\rho\\right] &= \\left[\\frac{c_s^2}{K}\\right]^n\n",
    "\\end{align}\n",
    "\n",
    "I think you found the inverse to be true."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TOV Solver as illustration\n",
    "The TOV equations are \n",
    "\\begin{eqnarray}\n",
    "\\frac{dP}{dr} &=& -\\mu\\frac{GM}{r^2}\\left(1 + \\frac P {\\mu c^2}\\right)\\left(1 + \\frac {4\\pi r^3 P}{Mc^2}\\right)\\left(1 - \\frac {2GM}{rc^2}\\right)^{-1}\\\\\n",
    "\\frac{dM}{dr} &=& 4\\pi \\mu r^2,\n",
    "\\end{eqnarray}\n",
    "Here $M$ is the rest mass measure by a distant observer when we take $r\\rightarrow \\infty$.  Note this is different from the mass measured by integrating the density over the volume\n",
    "\\begin{equation}\n",
    "M' = \\int_0^{\\infty} \\frac{4\\pi r^2\\mu}{\\sqrt{1 - \\frac {2 GM}{rc^2}}} dr\n",
    "\\end{equation}\n",
    "Additionally and annoyingly, $\\mu = \\rho h$ is the mass-energy density.  Alot of the literature uses $\\rho$ for this, which is incredibly annoying. \n",
    "\n",
    "$\\newcommand{\\muCode}{{\\tilde{\\mu}}}$\n",
    "\n",
    "In dimensionless units they are \n",
    "\\begin{eqnarray}\n",
    "\\frac{d\\PCode}{d\\rCode} &=& -\\frac {\\left(\\muCode + \\PCode\\right)\\left(\\MCode + 4\\pi \\rCode^3 \\PCode\\right)}{\\rCode^2\\left(1 - \\frac {2\\MCode}{\\rCode}\\right)}\\\\\n",
    "\\frac{d\\MCode}{d\\rCode} &=& 4\\pi \\muCode\\rCode^2\n",
    "\\end{eqnarray}\n",
    "\n",
    "At this point, we need to discuss how to numerically integrate these models.  First we pick a central baryonic mass density $\\rhoCode_{0,c}$, then we compute a central pressure $\\PCode_c$ and central mass-energy density $\\muCode_c$.  At $\\rCode=0$, we assume that $\\muCode=\\muCode_c$ is a constant and so \n",
    "\\begin{eqnarray}\n",
    "\\frac{d\\PCode}{d\\rCode} &=& -\\frac {\\left(\\muCode_c + \\PCode_c\\right)\\left(\\MCode(\\rCode \\ll 1) + 4\\pi \\rCode^3 \\PCode_c\\right)}{\\rCode^2\\left(1 - \\frac {2\\MCode(\\rCode \\ll 1)}{\\rCode}\\right)}\\\\\n",
    "\\frac{d\\MCode}{d\\rCode} &=& 4\\pi \\muCode_c\\rCode^2 \\rightarrow \\MCode(\\rCode \\ll 1) = \\frac{4\\pi}{3} \\muCode_c \\rCode^3 \n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another dimensionless prescription\n",
    "\n",
    "Let consider an alternative formulation where rather than setting $K=1$, we set the characteristic mass $\\Mfid = M_0$.  In this case,\n",
    "\\begin{eqnarray}\n",
    "r_0 &=& \\frac{GM_0}{c^2} \\\\\n",
    "t_0 &=& \\frac{GM_0}{c^3} \\\\\n",
    "\\rho_0 &=& \\frac{M_0}{r_0^3} = \\frac{c^6}{G^3 M_0^2} = 6.17\\times 10^{17}\\left(\\frac {M_0} {1 M_{\\odot}}\\right)^{-2}\n",
    "\\end{eqnarray}\n",
    "In this case we can define $\\rhoCode = \\rho/\\rho_0$, $\\rCode = r/r_0$, $t_0 = t/t_0$. The only remaining thing to do is to define $\\PCode$.  Lets define $P_0'$ to be the pressure in dimensionful units at $\\rho_0$ (density in units of $1/M_0^2$):\n",
    "\\begin{equation}\n",
    "P = P_0'\\rhoCode^{1+1/n} \\rightarrow P_0' = K\\rho_0^{1+1/n},\n",
    "\\end{equation}\n",
    "So defining $P_0 = \\rho_0 c^2$, we have\n",
    "\\begin{equation}\n",
    "\\PCode = \\frac{P}{P_0} = \\frac{K\\rho_0^{1/n}}{c^2}\\rhoCode^{1+1/n} = \\PCode_0\\rhoCode^{1+1/n}\n",
    "\\end{equation}\n",
    "If we take $K=1$ and define $\\rho_0$ such that the $\\PCode_0 = 1$, we recover the results above.\n",
    "Finally for $\\muCode = \\rhoCode + \\PCode/n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## metric for TOV equation\n",
    "\n",
    "The metric for the TOV equation (taken) from wikipedia is\n",
    "\\begin{equation}\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",
    "\\end{equation}\n",
    "where $M$ is defined as above, the mass as measured by a distant observer.  The equation for $\\nu$ is\n",
    "\\begin{equation}\n",
    "\\frac{d\\nu}{dr} = -\\left(\\frac {2}{P +\\mu}\\right)\\frac{dP}{dr}\n",
    "\\end{equation}\n",
    "with the boundary condition\n",
    "\\begin{equation}\n",
    "\\exp(\\nu) = \\left(1-\\frac {2Gm(R)}{Rc^2}\\right)\n",
    "\\end{equation}\n",
    "\n",
    "Lets write this in dimensionless units:\n",
    "\\begin{equation}\n",
    "ds^2 = \\exp(\\nu) d\\tCode^2 - \\left(1 - \\frac{2\\MCode}{\\rCode}\\right)^{-1} d\\rCode^2 + \\rCode^2 d\\Omega^2\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "\\frac{d\\nu}{d\\rCode} = -\\left(\\frac {2}{\\PCode +\\muCode}\\right)\\frac{d\\PCode}{d\\rCode}\n",
    "\\end{equation}\n",
    "and BC:\n",
    "\\begin{equation}\n",
    "\\exp(\\nu) = \\left(1-\\frac {2\\MCode}{\\rCode}\\right)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD6CAYAAACxrrxPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3hVZb728e8vtIwIMkAQhmYZRBBFINKLBTkUBaVjVBRCKHrEYZDhHc4gyhwRHcSDIqGpBJQioiJSwpBIiwIBKcbIGFQUkEFQAigtyfP+kaiRAbOBvbP23rk/15Xr2mVlr3uRcPPwrGbOOUREJPRFeB1ARET8Q4UuIhImVOgiImFChS4iEiZU6CIiYUKFLiISJgosdDN72cwOmNnH53g/xsy2m9kOM0sxs/r+jykiIgWxgo5DN7PWwDEgwTlX7yzvNwfSnXPfm1kHYIxzrklBK65YsaK74oorLiy1iEgRtXnz5oPOuaizvVe8oG92zq0xsyt+4/2UfE8/BKr5EuqKK64gNTXVl0VFRCSPme0+13v+nkPvDyzz82eKiIgPChyh+8rMbiG30Fv+xjJxQBxAjRo1/LVqERHBTyN0M7sBmAF0cc4dOtdyzrlpzrlo51x0VNRZp4BEROQCXXShm1kNYBFwn3PuXxcfSURELkSBUy5mNhe4GahoZnuAx4ESAM65eGA0UAF4ycwAspxz0YEKLCIiZ+fLUS59Cng/Foj1WyIREbkgOlNURCRMqNBFJCStXr2anj17kpGR4XWUoKFCF5GQs3DhQtq1a8cbb7xBQkKC13GChgpdRELKlClT6NmzJ6dOnQJgz549HicKHip0EQkJzjkef/xxhgwZgnOODh06ALB3716PkwUPv50pKiISKNnZ2Tz00ENMnTqViIgIpk2bxk033cSyZcs0Qs9HhS4iQe3EiRPExMSwaNEiIiMjmT9/Pp07d+bQodyT0lXov1Chi0jQOnz4MF26dGHNmjVcdtllvPvuu7Rq1QqA8uXLExkZyZEjRzhy5Ahly5b1OK33NIcuIkFp3759tG7dmjVr1vCHP/yBdevW/VzmAGZGtWq5V+vWPHouFbqIBJ2dO3fSvHlzduzYQe3atUlJSaFevf+4v44K/QwqdBEJKhs3bqRFixbs3r2bJk2asG7dOmrWrHnWZX8qdM2j51Khi0jQWL58ObfccguHDh2iY8eOrFq1iooVK55zeRX6r6nQRSQozJ49mzvvvJMff/yRvn378vbbb1O6dOnf/J6qVasCKvSfqNBFxHP/+Mc/uP/++8nKyuIvf/kLr7zyCiVKlCjw+zRC/zUdtiginsnJyWHEiBFMmDABgIkTJ/Loo4/6/P0q9F9ToYuIJ06fPk2/fv2YM2cOJUqU4NVXX+Wee+45r8/QUS6/pkIXkUJ37NgxunfvzooVKyhdujSLFi2iXbt25/05lSpVonjx4hw8eJATJ04QGRkZgLShQ3PoIlKovv32W2699VZWrFhBVFQU77///gWVOUBERMTPO0Y1Slehi0gh+vLLL2nRogWbNm3iyiuvZP369URHX9wtiHWkyy9U6CJSKLZv307z5s357LPPqF+/PikpKdSqVeuiP1c7Rn+hQheRgFu9ejWtWrXim2++4eabb2b16tVUrlzZL5+tQv+FCl1EAmrRokX813/9F0eOHKF79+4sX76cyy67zG+fryNdfqFCF5GAiY+Pp3v37pw8eZKHHnqIefPmUapUKb+uQyP0X6jQRcTvnHOMGTOGwYMH45xj7NixvPDCCxQrVszv69JO0V/oOHQR8avs7GyGDBnCtGnTiIiIYOrUqcTGxgZsfRqh/6LAEbqZvWxmB8zs43O8b2Y2ycwyzGy7mTX0f0wRCQUnTpygR48eTJs2jcjISBYtWhTQMgeoUqUKZsb+/fs5ffp0QNcV7HyZcnkVaP8b73cAauV9xQFTLj6WiISaw4cP065dO9566y3KlSvHypUr6dKlS8DXW6JECSpXroxzjv379wd8fcGswEJ3zq0BvvuNRboACS7Xh0A5M6vir4AiEvx+ul3c2rVrqVq1KmvXrqVly5aFtn5Nu+Tyx07RqsDX+Z7vyXtNRIqA/LeLu/baa895u7hAUqHnKtSjXMwszsxSzSz122+/LcxVi0gA5L9dXNOmTVm3bh01atQo9Bw60iWXPwp9L1A93/Nqea/9B+fcNOdctHMuOioqyg+rFhGvvPfee7+6Xdw///lPKlSo4EkWjdBz+aPQFwP35x3t0hTIdM5944fPFZEgNXPmTLp06XJet4sLpFAq9C1btvDZZ58F5LN9OWxxLvABUNvM9phZfzMbZGaD8hZZCnwOZADTgSEBSSoinnPO8eSTTxIbG0t2djajRo3y+XZxgRQKp///8MMPDB8+nJtuuonY2FhycnL8vo4CTyxyzvUp4H0HPOS3RCISlLKyshgyZAjTp08nIiKCF198kcGDB3sdCwj+EfqKFSsYNGgQX375JRERETRs2JDTp0/7/TIIOlNURAr0448/0rt3b959910iIyOZO3cud911l9exfpb/Jhc5OTlERATHVU0OHDjAn/70J15//XUA6tevz/Tp07npppsCsr7g2GoRCVoHDx7k1ltv5d1336V8+fKsWrUqqMocIDIykgoVKpCVlcWBAwe8joNzjlmzZlGnTh1ef/11IiMjGT9+PJs2bQpYmYNG6CLyGz7//HPat2/PZ599Rs2aNVm+fDnXXnut17HOqlq1ahw6dIg9e/b47VrrFyIjI4OBAweSlJQEQNu2bYmPj+fqq68O+Lo1QheRs9qyZct/3GEoWMscvJ9HP336NE8//TTXX389SUlJVKhQgYSEBBITEwulzEEjdBE5i8TERLp168axY8e47bbbWLRoEWXLlvU61m/y8kiXjRs3MmDAALZv3w7Afffdx4QJEyjs8200QheRX0lISKBTp04cO3aMmJgYli5dGvRlDt6M0I8ePcrQoUNp2rQp27dv58orryQxMZGEhIRCL3NQoYtIHucc48aNo2/fvmRlZTFixAgSEhIoWbKk19F8Utin/y9ZsoTrrruOSZMmERERwYgRI/j444+5/fbbC2X9Z6MpFxEhOzubRx55hJdeegkz4/nnn+eRRx7xOtZ5KawR+v79+xk6dCgLFiwAoFGjRkyfPp0GDRoEdL2+UKGLFHHHjx8nJiaGt956i5IlSzJnzhx69OjhdazzFuhCz8nJYebMmYwYMYLDhw9TunRp/v73v/Pwww9TvHhwVGlwpBART3z33Xd07tyZ9evXU65cOd5++23atGnjdawLkn+nqHMOM/PbZ3/66acMHDiQNWvWANChQwemTJlCzZo1/bYOf9AcukgRtXv3blq2bMn69eupVq0a69atC9kyByhTpgxly5bl+PHjfP/99375zFOnTjF27Fjq16/PmjVrqFSpEnPnzuW9994LujIHFbpIkbRt2zaaNWtGeno69erV44MPPuC6667zOtZF8+e0S0pKCg0aNGD06NGcOnWKfv36kZ6eTu/evf06+vcnFbpIEZOYmEirVq345ptvaNOmDWvXrv25CEOdP450yczMZMiQIbRs2ZJPPvmEWrVqkZSUxMyZMylfvry/ogaECl2kCHnllVfo1KkTR48epWfPnixfvpxy5cp5HctvLnaE/tZbb1G3bl2mTJlCsWLFGDVqFNu3b+eWW27xZ8yAUaGLFAHOOcaMGUO/fv1+PsZ87ty5REZGeh3Nry600Pfu3UvXrl3p2rUr+/bto2nTpnz00Uf8/e9/D6k/Ix3lIhLmTp06RVxcHLNmzQq665j72/me/p+Tk0N8fDwjR47k6NGjlClThnHjxjFo0CCKFSsWyKgBoUIXCWOZmZl069aNVatWcckllzB//nzuuOMOr2MFzPmM0NPS0oiLiyMlJQWAzp07M3ny5JDen6BCFwlTX3/9NR07duTjjz/m8ssvZ8mSJURHR3sdK6B82Sl64sQJnnrqKZ5++mlOnz5NlSpVePHFF7n77ruD9ugVX6nQRcLQ1q1b6dSpE/v27ePaa69l2bJlXHHFFV7HCriCRuhr1qwhLi6OnTt3AjBw4ECefvrpsNkxrJ2iImFmxYoVtGrVin379tG6dWtSUlKKRJkDlC9fnsjISI4cOcKRI0d+fv37779nwIABtGnThp07d1KnTh3Wrl1LfHx82JQ5qNBFwsrLL7/886Vve/fuTWJiIr///e+9jlVozOw/LgGwYMEC6tSpw4wZMyhZsiRjxozho48+omXLlh6n9T8VukgYcM4xevRo+vfvT3Z2NiNHjuS1117z+13lQ8FPhf7BBx/QuXNnevXqxb///W9atmzJ1q1befzxx8P2z0Vz6CIh7tSpU8TGxjJ79mwiIiJ46aWXGDhwoNexPPNToffv3x+Ayy67jGeeeYbY2FgiIsJ7DKtCFwlhhw8fplu3biQlJVG6dGnmz59Pp06dvI7lqfyHHXbv3p1JkyZRpUoVDxMVHhW6SIj66quv6NixI2lpaVSuXJklS5bQqFEjr2N57oEHHmDXrl3ce++9dO7c2es4hUqFLhKCUlNTufPOO9m/fz916tRh2bJlQXk5Vy/Url3757sJFTU+TSiZWXsz22lmGWY28izv1zCzZDP7yMy2m1lH/0cVEYC3336b1q1bs3//fm655RbWr1+vMhfAh0I3s2LAZKADUBfoY2Z1z1jsf4AFzrkGQG/gJX8HFSnqnHNMmDCBrl27cvz4cR588EGWL19epA5LlN/mywi9MZDhnPvcOXcKmAd0OWMZB5TNe3wZsM9/EUXk9OnTDB48mOHDh+Oc46mnnmLmzJmULFnS62gSRHyZQ68KfJ3v+R6gyRnLjAESzey/gdJA27N9kJnFAXEANWrUON+sIkVSZmYmPXv2JDExkVKlSpGQkEDPnj29jiVByF8HZfYBXnXOVQM6ArPN7D8+2zk3zTkX7ZyLjoqK8tOqRcLXl19+SYsWLUhMTCQqKor3339fZS7n5Euh7wWq53teLe+1/PoDCwCccx8AkUBFfwQUKao2btxIkyZNSEtLo06dOmzYsIGmTZt6HUuCmC+FvgmoZWZXmllJcnd6Lj5jma+A2wDMrA65hf6tP4OKFCULFy6kTZs2HDhwgNtuu42UlBSuvPJKr2NJkCuw0J1zWcDDwAogndyjWdLM7Ekz++mo/T8DA8xsGzAXeMA55wIVWiRcOecYP348PXr04MSJE8TGxrJs2bKwuiKgBI5PJxY555YCS894bXS+x58ALfwbTaRo+elIlpkzZwLwzDPPMHz48JC/6YIUHp0pKhIE8l+T5Xe/+x1z5syha9euXseSEKNCF/HYrl27uPPOO0lPT+fyyy/n3Xff5aabbvI6loSg8L6WpEiQW716NY0bNyY9PZ169eqxYcMGlblcMBW6iEdmzJhB27Zt+e677+jUqZOuySIXTYUuUsiys7P505/+xIABA8jKyuLPf/4z77zzDmXLli34m0V+g+bQRQpRZmYmffr0YdmyZZQoUYL4+Hj69evndSwJEyp0kUKSf+dnhQoVWLRoEa1bt/Y6loQRTbmIFILVq1fTpEkT0tPTqVu3Lhs3blSZi9+p0EUCbObMmbRt25ZDhw7RsWNHPvjgA6666iqvY0kYUqGLBEh2djbDhg0jNjb2552fixcv1s5PCRjNoYsEgHZ+ihdU6CJ+pp2f4hVNuYj4kXZ+ipdU6CJ+4JxjypQp2vkpnlKhi1ykU6dOMXDgQIYMGUJWVhbDhw/Xzk/xhObQRS7C/v376datGykpKURGRjJjxgxiYmK8jiVFlApd5AKlpqZy1113sXfvXqpVq8bbb79No0aNvI4lRZimXEQuwJw5c2jVqhV79+6lRYsWpKamqszFcyp0kfPw0xz5fffdx4kTJ4iLiyMpKYnLL7/c62gimnIR8dX3339P7969SUxMpHjx4rzwwgsMGjTI61giP1Ohi/ggLS2NLl26sGvXLqKioli4cKGOL5egoykXkQK88847NG3alF27dtGgQQNSU1NV5hKUVOgi55CTk8OTTz7JXXfdxbFjx+jduzfr1q2jRo0aXkcTOStNuYicxdGjR3nggQdYtGgRZsa4ceMYMWIEZuZ1NJFzUqGLnOHTTz+la9eupKenc9lll/H666/TsWNHr2OJFMinKRcza29mO80sw8xGnmOZnmb2iZmlmdnr/o0pUjjeeustGjdu/KuLa6nMJVQUWOhmVgyYDHQA6gJ9zKzuGcvUAv4f0MI5dx3waACyigRMdnY2f/3rX+natStHjx6lR48ebNiwgWuuucbraCI+82XKpTGQ4Zz7HMDM5gFdgE/yLTMAmOyc+x7AOXfA30FFAuXgwYPcc889rFy5kmLFijF+/HiGDRum+XIJOb4UelXg63zP9wBNzljmGgAzWw8UA8Y455af+UFmFgfEATpSQILCli1b6Nq1K7t37yYqKor58+dzyy23eB1L5IL467DF4kAt4GagDzDdzMqduZBzbppzLto5Fx0VFeWnVYtcmFdeeYXmzZuze/duGjduzObNm1XmEtJ8KfS9QPV8z6vlvZbfHmCxc+60c+4L4F/kFrxI0Dl58iSDBw+mX79+nDx5kri4ONasWUP16tUL/maRIOZLoW8CapnZlWZWEugNLD5jmbfJHZ1jZhXJnYL53I85Rfxi79693HzzzcTHx1OqVClmzJjB1KlTKVWqlNfRRC5agXPozrksM3sYWEHu/PjLzrk0M3sSSHXOLc57r52ZfQJkA4855w4FMrjI+Vq9ejU9e/bkwIEDVK9enUWLFhEdHe11LBG/MeecJyuOjo52qampnqxbihbnHM8//zyPPfYY2dnZ3HbbbcydOxftx5FQZGabnXNnHYnoWi4S1o4ePUqfPn0YNmwY2dnZjBgxguXLl6vMJSzp1H8JWx9//DHdu3dn586dXHrppbz66qt069bN61giAaMRuoSlWbNm0bhxY3bu3Em9evVITU1VmUvYU6FLWDl+/DgDBgzggQce4Pjx49x///1s2LCB2rVrex1NJOA05SJhIyMjgx49erB161ZKlSrFiy++SP/+/XUKvxQZKnQJC4sWLeLBBx/kyJEjXH311SxcuJAbb7zR61gihUpTLhLSTp8+zbBhw+jWrRtHjhzh7rvvZvPmzSpzKZI0QpeQtWfPHnr16kVKSgrFixfnmWee4dFHH9UUixRZKnQJSYmJicTExHDw4EGqVq3KggULaN68udexRDylKRcJKdnZ2YwZM4b27dtz8OBB2rVrx0cffaQyF0EjdAkhBw4c4N5772XlypWYGU888QSjRo2iWLFiXkcTCQoqdAkJSUlJxMTEsH//fipWrMjrr7/O7bff7nUskaCiKRcJallZWYwePZq2bduyf/9+WrduzdatW1XmImehEboErT179nDPPfewdu1azIzRo0fzt7/9jeLF9Wsrcjb6myFB6b333qNv374cOnSIypUr89prr3Hrrbd6HUskqGnKRYLKqVOnGD58OHfccQeHDh2iXbt2bNu2TWUu4gON0CVofPHFF/Tu3ZuNGzdSrFgx/vd//5fHHnuMiAiNO0R8oUKXoLBw4UJiY2PJzMykRo0azJ07V8eWi5wnDX3EUydOnGDIkCH06NGDzMxM7rrrLp0oJHKBNEIXz3z66af06tWL7du3U7JkSf7xj3/w8MMP61osIhdIhS6FzjnH7NmzGTJkCD/88AN//OMfmT9/Pg0bNvQ6mkhI05SLFKrMzExiYmLo27cvP/zwA3369GHz5s0qcxE/0AhdCk1KSgr33HMPu3fvpnTp0kyaNIkHH3xQUywifqIRugRcVlYWTzzxBK1atWL37t00atSILVu20K9fP5W5iB9phC4BtXv3bmJiYli/fj1mxogRIxg7diwlS5b0OppI2PFphG5m7c1sp5llmNnI31ium5k5M4v2X0QJVfPmzaN+/fqsX7+eKlWqsHLlSsaPH68yFwmQAgvdzIoBk4EOQF2gj5nVPctyZYChwAZ/h5TQcvToUR588EH69OlDZmYmXbp0Yfv27dx2221eRxMJa76M0BsDGc65z51zp4B5QJezLDcWGA+c8GM+CTEbN26kQYMGvPrqq/zud79jypQpvPXWW1SsWNHraCJhz5dCrwp8ne/5nrzXfmZmDYHqzrn3/JhNQkh2djZPP/00LVq0YNeuXdxwww2kpqYyaNAg7fgUKSQXvVPUzCKA54AHfFg2DogDqFGjxsWuWoLEnj17uP/++0lOTgbg0UcfZdy4cURGRnqcTKRo8WWEvheonu95tbzXflIGqAe8b2ZfAk2BxWfbMeqcm+aci3bORUdFRV14agkab775JvXr1yc5OZlKlSqxdOlSJk6cqDIX8YAvhb4JqGVmV5pZSaA3sPinN51zmc65is65K5xzVwAfAp2dc6kBSSxBITMzk759+9K9e3e+++472rdvz/bt2+nQoYPX0USKrAIL3TmXBTwMrADSgQXOuTQze9LMOgc6oASfNWvWUL9+fRISEoiMjOTFF19k6dKlXH755V5HEynSfJpDd84tBZae8drocyx788XHkmB08uRJRo8ezbPPPotzjkaNGjFnzhyuvfZar6OJCDpTVHyUlpZGTEwM27ZtIyIiglGjRjF69GhKlCjhdTQRyaNCl9+Uk5PDpEmTGDlyJCdPnuTqq69m9uzZNGvWzOtoInIGFbqc0549e3jggQdYtWoVALGxsUycOJFLL73U42QicjYqdDmrefPmMXjwYA4fPkxUVBTTp0+nS5eznSAsIsFCl8+VXzl8+DAxMTH06dOHw4cPc8cdd7Bjxw6VuUgI0AhdfpacnEzfvn35+uuvueSSS5g4cSIDBgzQqfsiIUIjdOHHH39k6NCh3HrrrXz99dc0adKEbdu2ERcXpzIXCSEq9CLugw8+4MYbb2TSpEkUL16cJ598knXr1vHHP/7R62gicp405VJEnTx5kscff5xnn32WnJwc6tWrR0JCAg0aNPA6mohcII3Qi6AtW7YQHR3N+PHjARg5ciSpqakqc5EQpxF6EXL69GnGjRvH2LFjycrK4pprrmHWrFk0bdrU62gi4gcq9CIiLS2Nvn37snnzZgCGDh3KU089xSWXXOJxMhHxF025hLns7GyeffZZGjZsyObNm6lZsyZJSUk8//zzKnORMKMRehjLyMigb9++pKSkADBgwAAmTJhAmTJlPE4mIoGgEXoYysnJYfLkydSvX5+UlBSqVKnC0qVLmTZtmspcJIxphB5mvvrqK/r16/fzBbViYmKYNGkS5cuX9ziZiASaRuhhwjnHjBkzqFevHqtWraJixYosXLiQOXPmqMxFigiN0MPA7t27GTBgACtXrgTg7rvvJj4+nkqVKnmcTEQKk0boISwnJ4f4+Hjq1avHypUrqVChAnPnzuXNN99UmYsUQRqhh6gvvviC2NhYkpKSAOjWrRuTJ0/WjZpFijCN0EPMT0ewXH/99SQlJVGxYkUWLFjAwoULVeYiRZxG6CFk165d9O/fn9WrVwPQq1cvXnjhBaKiojxOJiLBQCP0EJCTk8P//d//cf3117N69WoqVarEm2++ybx581TmIvIzjdCD3L/+9S/69evH+vXrAbjnnnuYNGkSFSpU8DiZiAQbjdCDVHZ2NhMmTKB+/fqsX7+eypUr88477/Daa6+pzEXkrHwqdDNrb2Y7zSzDzEae5f1hZvaJmW03s1VmVtP/UYuOTz/9lJYtWzJ8+HBOnDjB/fffT1paGp07d/Y6mogEsQIL3cyKAZOBDkBdoI+Z1T1jsY+AaOfcDcBC4Bl/By0KsrKyGD9+PDfeeCMffvghf/jDH1iyZAmzZs3S2Z4iUiBfRuiNgQzn3OfOuVPAPKBL/gWcc8nOuR/znn4IVPNvzPC3fft2mjVrxsiRIzl58iT9+vUjLS2NTp06eR1NREKEL4VeFfg63/M9ea+dS39g2cWEKkpOnjzJ3/72Nxo1akRqairVq1dn2bJlzJw5k3LlynkdT0RCiF+PcjGze4FooM053o8D4gBq1Kjhz1WHpJSUFGJjY0lPTwfgoYceYty4cbrErYhcEF9G6HuB6vmeV8t77VfMrC0wCujsnDt5tg9yzk1zzkU756KL8vHTx44d45FHHqFly5akp6dTu3Zt1q5dy4svvqgyF5EL5kuhbwJqmdmVZlYS6A0szr+AmTUAppJb5gf8HzN8rFixgnr16vHCCy8QERHBX//6V7Zu3UrLli29jiYiIa7AKRfnXJaZPQysAIoBLzvn0szsSSDVObcYeBa4FHjDzAC+cs7pGLt8Dh06xLBhw0hISACgYcOGzJw5kxtvvNHjZCISLnyaQ3fOLQWWnvHa6HyP2/o5V9hwzrFw4UIefvhhDhw4QGRkJE888QTDhg2jeHGdqCsi/qNGCaB9+/YxZMgQ3nnnHQBat27NjBkzqFWrlsfJRCQc6dT/APjpdnB169blnXfeoUyZMsTHx5OcnKwyF5GA0QjdzzIyMoiLiyM5ORmAO+64gylTplCtms61EpHA0gjdT7KyspgwYQI33HADycnJVKxYkblz57J48WKVuYgUCo3Q/WDHjh3079+fTZs2AXDvvfcyceJEKlas6HEyESlKNEK/CMePH2fUqFE0bNiQTZs2Ub16dd577z1mz56tMheRQqcR+gVKTk4mLi6OjIwMzIwhQ4Ywbtw4ypYt63U0ESmiVOjn6dChQzz22GO88sorAFx33XVMmzaN5s2be5xMRIo6Tbn4yDnH3LlzqVOnDq+88golS5Zk7NixbNmyRWUuIkFBI3QffPnllwwePJjly5cD0KZNG6ZOnUrt2rU9TiYi8guN0H9DVlYWzz33HNdddx3Lly+nXLlyzJgxg+TkZJW5iAQdjdDPYcuWLQwYMIAtW7YA0KtXL55//nkqV67scTIRkbPTCP0MP/zwA4899hiNGzdmy5Yt1KhRgyVLljBv3jyVuYgENY3Q81mxYgWDBw/miy++ICIigkcffZSxY8dy6aWXeh1NRKRAKnTgwIEDDBs2jNdeew2A+vXrM336dG666SaPk4mI+K5IT7k453j11VepU6cOr732GpGRkYwfP55Nm5DVcpEAAAaZSURBVDapzEUk5BTZEXpGRgYDBw4kKSkJgLZt2xIfH8/VV1/tcTIRkQtT5Ebop0+f5umnn+b6668nKSmJChUqkJCQQGJiospcREJakRqhb9iwgQEDBrBjxw4A7rvvPp577jldSEtEwkKRGKEfPXqURx55hGbNmrFjxw6uuuoqEhMTSUhIUJmLSNgI+0J/9913qVu3Li+88AIRERH85S9/YceOHdx+++1eRxMR8auwnXL55ptvGDp0KG+88QYA0dHRTJ8+nRtvvNHjZCIigRF2I/ScnBymTZtGnTp1eOONNyhdujQTJ07kww8/VJmLSFgLqxH6p59+SlxcHGvXrgWgY8eOvPTSS9SsWdPjZCIigRcWI/STJ0/yxBNPUL9+fdauXUulSpWYN28eS5YsUZmLSJHhU6GbWXsz22lmGWY28izvlzKz+XnvbzCzK/wd9FzWrVtHgwYNGDNmDKdOnSI2Npb09HR69eqFmRVWDBERzxVY6GZWDJgMdADqAn3MrO4Zi/UHvnfO/RGYCIz3d9AzHT58mEGDBtGqVSvS09O55pprSE5OZvr06ZQvXz7QqxcRCTq+jNAbAxnOuc+dc6eAeUCXM5bpAszKe7wQuM0CNDx2zvHmm29St25dpk6dSvHixfmf//kftm3bxs033xyIVYqIhARfdopWBb7O93wP0ORcyzjnsswsE6gAHPRHyPwefPBBZs3K/bejWbNmTJs2jXr16vl7NSIiIadQd4qaWZyZpZpZ6rfffntBn9G8eXPKlCnDSy+9xLp161TmIiJ5fBmh7wWq53teLe+1sy2zx8yKA5cBh878IOfcNGAaQHR0tLuQwLGxsXTu3Fl3DxIROYMvI/RNQC0zu9LMSgK9gcVnLLMY6Jv3uDuQ5Jy7oMIuSEREhMpcROQsChyh582JPwysAIoBLzvn0szsSSDVObcYmAnMNrMM4DtyS19ERAqRT2eKOueWAkvPeG10vscngB7+jSYiIucjLM4UFRERFbqISNhQoYuIhAkVuohImFChi4iECQvQ4eIFr9jsW2D3BX57RQJwWYEgEa7bpu0KLdqu4FXTORd1tjc8K/SLYWapzrlor3MEQrhum7YrtGi7QpOmXEREwoQKXUQkTIRqoU/zOkAAheu2abtCi7YrBIXkHLqIiPynUB2hi4jIGYK60IP55tQXw4ftGmZmn5jZdjNbZWY1vch5vgrarnzLdTMzZ2Yhc7SBL9tmZj3zfm5pZvZ6YWe8ED78LtYws2Qz+yjv97GjFznPh5m9bGYHzOzjc7xvZjYpb5u3m1nDws4YMM65oPwi91K9u4CrgJLANqDuGcsMAeLzHvcG5nud20/bdQtwSd7jweGyXXnLlQHWAB8C0V7n9uPPrBbwEfD7vOeVvM7tp+2aBgzOe1wX+NLr3D5sV2ugIfDxOd7vCCwDDGgKbPA6s7++gnmEHlQ3p/ajArfLOZfsnPsx7+mH5N4lKtj58vMCGAuMB04UZriL5Mu2DQAmO+e+B3DOHSjkjBfCl+1yQNm8x5cB+wox3wVxzq0h974M59IFSHC5PgTKmVmVwkkXWMFc6Ge7OXXVcy3jnMsCfro5dTDzZbvy60/uaCLYFbhdef+1re6ce68wg/mBLz+za4BrzGy9mX1oZu0LLd2F82W7xgD3mtkecu+J8N+FEy2gzvfvYMjw6QYX4g0zuxeIBtp4neVimVkE8BzwgMdRAqU4udMuN5P7P6o1Zna9c+6wp6kuXh/gVefcBDNrRu6dyeo553K8Dib/KZhH6Odzc2p+6+bUQcaX7cLM2gKjgM7OuZOFlO1iFLRdZYB6wPtm9iW5c5eLQ2THqC8/sz3AYufcaefcF8C/yC34YObLdvUHFgA45z4AIsm9Hkoo8+nvYCgK5kIPqptT+1GB22VmDYCp5JZ5KMzFQgHb5ZzLdM5VdM5d4Zy7gtx9A52dc6nexD0vvvwuvk3u6Bwzq0juFMznhRnyAviyXV8BtwGYWR1yC/3bQk3pf4uB+/OOdmkKZDrnvvE6lF94vVf2t77I3Rv9L3L3xI/Ke+1JcosAcn+53gAygI3AVV5n9tN2/RP4N7A172ux15n9sV1nLPs+IXKUi48/MyN3SukTYAfQ2+vMftquusB6co+A2Qq08zqzD9s0F/gGOE3u/5z6A4OAQfl+VpPztnlHKP0eFvSlM0VFRMJEME+5iIjIeVChi4iECRW6iEiYUKGLiIQJFbqISJhQoYuIhAkVuohImFChi4iEif8PO2OG5zvi8dgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Just generated a TOV star with r= 0.9565681425227097 , m = 0.14050303285288188 , m/r = 0.1468824086931645 .\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import scipy.integrate as si\n",
    "import math\n",
    "import matplotlib.pyplot as pl\n",
    "n = 1.\n",
    "rho_central = 0.129285\n",
    "\n",
    "P0 = 1. # ZACH NOTES: CHANGED FROM 100.\n",
    "gamma = 1. + 1./n\n",
    "gam1 = gamma - 1.\n",
    "\n",
    "def pressure( rho) :\n",
    "    return P0*rho**gamma\n",
    "\n",
    "def rhs( r, y) :\n",
    "# In \\tilde units\n",
    "#\n",
    "    P = y[0]\n",
    "    m   = y[1]\n",
    "    nu  = y[2]\n",
    "    rbar = y[3]\n",
    "\n",
    "    rho = (P/P0)**(1./gamma)\n",
    "    mu = rho + P/gam1\n",
    "    dPdr = 0.\n",
    "    drbardr = 0.\n",
    "    if( r < 1e-4 or m <= 0.) :\n",
    "        m = 4*math.pi/3. * mu*r**3\n",
    "        dPdr = -(mu + P)*(4.*math.pi/3.*r*mu + 4.*math.pi*r*P)/(1.-8.*math.pi*mu*r*r)\n",
    "        drbardr = 1./(1. - 8.*math.pi*mu*r*r)**0.5\n",
    "    else :\n",
    "        dPdr = -(mu + P)*(m + 4.*math.pi*r**3*P)/(r*r*(1.-2.*m/r))\n",
    "        drbardr = 1./(1. - 2.*m/r)**0.5*rbar/r\n",
    "\n",
    "    dmdr = 4.*math.pi*r*r*mu\n",
    "    dnudr = -2./(P + mu)*dPdr\n",
    "    return [dPdr, dmdr, dnudr, drbardr]\n",
    "\n",
    "def integrateStar( P, showPlot = False, dumpData = False, compareFile=\"TOV/output_EinsteinToolkitTOVSolver.txt\") :\n",
    "    integrator = si.ode(rhs).set_integrator('dop853')\n",
    "    y0 = [P, 0., 0., 0.]\n",
    "    integrator.set_initial_value(y0,0.)\n",
    "    dr = 1e-5\n",
    "    P = y0[0]\n",
    "\n",
    "    PArr = []\n",
    "    rArr = []\n",
    "    mArr = []\n",
    "    nuArr = []\n",
    "    rbarArr = []\n",
    "    r = 0.\n",
    "\n",
    "    while integrator.successful() and P > 1e-9*y0[0] :\n",
    "        P, m, nu, rbar = integrator.integrate(r + dr)\n",
    "        r = integrator.t\n",
    "\n",
    "        dPdr, dmdr, dnudr, drbardr = rhs( r+dr, [P,m,nu,rbar])\n",
    "        dr = 0.1*min(abs(P/dPdr), abs(m/dmdr))\n",
    "        dr = min(dr, 1e-2)\n",
    "        PArr.append(P)\n",
    "        rArr.append(r)\n",
    "        mArr.append(m)\n",
    "        nuArr.append(nu)\n",
    "        rbarArr.append( rbar)\n",
    "\n",
    "    M = mArr[-1]\n",
    "    R = rArr[-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]/rArr[-1])\n",
    "\n",
    "    rArrExtend_np = 10.**(np.arange(0.01,5.0,0.01))*rArr[-1]\n",
    "\n",
    "    rArr.extend(rArrExtend_np)\n",
    "    mArr.extend(rArrExtend_np*0. + M)\n",
    "    PArr.extend(rArrExtend_np*0.)\n",
    "    phiArr_np = np.append( np.exp(nuArr_np), 1. - 2.*M/rArrExtend_np)\n",
    "    rbarArr.extend( 0.5*(np.sqrt(rArrExtend_np*rArrExtend_np-2.*M*rArrExtend_np) + rArrExtend_np - M))\n",
    "    # Appending a Python array does what one would reasonably expect.\n",
    "    #   Appending a numpy array allocates space for a new array with size+1,\n",
    "    #   then copies the data over... over and over... super inefficient.\n",
    "    mArr_np       = np.array(mArr)\n",
    "    rArr_np       = np.array(rArr)\n",
    "    PArr_np       = np.array(PArr)\n",
    "    rbarArr_np    = np.array(rbarArr)\n",
    "    rhoArr_np     = (PArr_np/P0)**(1./gamma)\n",
    "    confFactor_np = rArr_np/rbarArr_np\n",
    "    #confFactor_np = (1.0 / 12.0) * np.log(1.0/(1.0 - 2.0*mArr_np/rArr_np))\n",
    "    Grr_np        = 1.0/(1.0 - 2.0*mArr_np/rArr_np)\n",
    "    Gtt_np        = phiArr_np\n",
    "    if( showPlot) :\n",
    "        r,rbar,rprop,rho,m,phi = np.loadtxt( compareFile, usecols=[0,1,2,3,4,5],unpack=True)\n",
    "        pl.plot(rArr_np[rArr_np < r[-1]], rbarArr_np[rArr_np < r[-1]],lw=2,color=\"black\")\n",
    "        #pl.plot(r, rbar, lw=2,color=\"red\")\n",
    "\n",
    "        pl.show()\n",
    "\n",
    "\n",
    "        if( dumpData) :\n",
    "            np.savetxt( \"output.txt\", np.transpose([rArr_np,rhoArr_np,PArr_np,mArr_np,phiArr_np,confFactor_np,rbarArr_np]), fmt=\"%.15e\")\n",
    "            np.savetxt( \"metric.txt\", np.transpose([rArr_np, Grr_np, Gtt_np]),fmt=\"%.15e\")\n",
    "        #        np.savetxt( \"output.txt\", zip(rArr,rhoArr,mArr,phiArr), fmt=\"%12.7e\")\n",
    "\n",
    "#    return rArr[-1], mArr[-1], phiArr[-1]\n",
    "    return R, M\n",
    "\n",
    "mass = []\n",
    "radius = []\n",
    "\n",
    "R_TOV,M_TOV = integrateStar(pressure(rho_central), showPlot=True, dumpData=True)\n",
    "print(\"Just generated a TOV star with r= \"+str(R_TOV)+\" , m = \"+str(M_TOV)+\" , m/r = \"+str(M_TOV/R_TOV)+\" .\")\n",
    "\n",
    "#for rho0 in np.arange(0.01, 1., 0.01):\n",
    "#    r,m = integrateStar(pressure(rho0))\n",
    "#    mass.append(m)\n",
    "#    radius.append(r)\n",
    "\n",
    "#print(mass, radius)\n",
    "#pl.clf()\n",
    "#pl.plot(radius,mass)\n",
    "#pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAOxUlEQVR4nO3cf6jdd33H8efLZBHssLb2dmASlmzGjTscOk/jfiEdrjGF2YhWTBkudV2LGxnM/bOKg2AqqKCsKxRnNqtO0LYWt8UNF8o6kVEnOQm2Na1tr6HaxEGvTVU2/6jXvPfH/ZZd727u/d5fOTn5PB9w4J7P9/P93s+3n+Y+c87pbaoKSVJ7XjTqBUiSRsMASFKjDIAkNcoASFKjDIAkNcoASFKjNo56ActxxRVX1LZt20a9DEkaK8eOHfteVU3MHx+rAGzbto3hcDjqZUjSWEny7YXGfQtIkhplACSpUQZAkhplACSpUQZAkhplACSpUQZAkhplACSpUQZAkhplACSpUQZAkhplACSpUQZAkhplACSpUQZAkhplACSpUQZAkhplACSpUQZAkhplACSpUQZAkhplACSpUQZAkhplACSpUb0CkGR3kseTTCW5dYHjb0hyPMlMkuvnHduX5MnusW+Bcw8n+cbKb0GStBJLBiDJBuBO4FpgErghyeS8ad8BbgQ+O+/cy4EDwOuBncCBJJfNOf5W4L9XsX5J0gr1eQWwE5iqqpNV9TxwN7Bn7oSqeqqqHgbOzjv3TcD9VXWmqp4D7gd2AyT5WeDPgQ+s8h4kSSvQJwCbgafnPD/VjfWx2Lm3AR8FfrTYBZLckmSYZDg9Pd3z20qSljKSD4GTvAb4xar6h6XmVtWhqhpU1WBiYuI8rE6S2tAnAKeBrXOeb+nG+jjXub8BDJI8BfwH8KokX+55TUnSGugTgKPAjiTbk2wC9gKHe17/CLAryWXdh7+7gCNV9bGqekVVbQN+G3iiqq5e/vIlSSu1ZACqagbYz+wP88eAe6vqRJKDSa4DSHJVklPA24GPJznRnXuG2ff6j3aPg92YJGnEUlWjXkNvg8GghsPhqJchSWMlybGqGswf9zeBJalRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGmUAJKlRBkCSGtUrAEl2J3k8yVSSWxc4/oYkx5PMJLl+3rF9SZ7sHvu6sZck+Zck30xyIsmH1uZ2JEl9LRmAJBuAO4FrgUnghiST86Z9B7gR+Oy8cy8HDgCvB3YCB5Jc1h3+SFX9MvBa4LeSXLuK+5AkLVOfVwA7gamqOllVzwN3A3vmTqiqp6rqYeDsvHPfBNxfVWeq6jngfmB3Vf2oqv69O/d54DiwZZX3Iklahj4B2Aw8Pef5qW6sjyXPTfIy4M3Av/W8piRpDYz0Q+AkG4HPAXdU1clzzLklyTDJcHp6+vwuUJIuYn0CcBrYOuf5lm6sj6XOPQQ8WVW3n+sCVXWoqgZVNZiYmOj5bSVJS+kTgKPAjiTbk2wC9gKHe17/CLAryWXdh7+7ujGSfAC4FPiz5S9bkrRaSwagqmaA/cz+4H4MuLeqTiQ5mOQ6gCRXJTkFvB34eJIT3blngNuYjchR4GBVnUmyBXgfs/9V0fEkX0/yR+twf5Kkc0hVjXoNvQ0GgxoOh6NehiSNlSTHqmowf9zfBJakRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRvUKQJLdSR5PMpXk1gWOvyHJ8SQzSa6fd2xfkie7x745469L8kh3zTuSZPW3I0nqa8kAJNkA3AlcC0wCNySZnDftO8CNwGfnnXs5cAB4PbATOJDksu7wx4CbgR3dY/eK70KStGwbe8zZCUxV1UmAJHcDe4BHX5hQVU91x87OO/dNwP1VdaY7fj+wO8mXgZdW1X92438PvAX40mpu5lze/8UTPPrdH67HpSVp3U2+4qUcePOvrPl1+7wFtBl4es7zU91YH+c6d3P39ZLXTHJLkmGS4fT0dM9vK0laSp9XACNVVYeAQwCDwaBWco31KKckjbs+rwBOA1vnPN/SjfVxrnNPd1+v5JqSpDXQJwBHgR1JtifZBOwFDve8/hFgV5LLug9/dwFHquq/gB8m+fXuv/75A+CfVrB+SdIKLRmAqpoB9jP7w/wx4N6qOpHkYJLrAJJcleQU8Hbg40lOdOeeAW5jNiJHgYMvfCAM/Anwd8AU8C3W6QNgSdLCUrWit9VHYjAY1HA4HPUyJGmsJDlWVYP54/4msCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqMMgCQ1ygBIUqN6BSDJ7iSPJ5lKcusCx1+c5J7u+NeSbOvGNyX5ZJJHkjyU5Oo559zQjT+c5F+TXLFG9yRJ6mHJACTZANwJXAtMAjckmZw37Sbguap6JfBXwIe78ZsBqurVwDXAR5O8KMlG4K+B36mqXwUeBvavwf1Iknrq8wpgJzBVVSer6nngbmDPvDl7gE93X98HvDFJmA3GAwBV9QzwfWAApHtc0s17KfDdVd6LJGkZ+gRgM/D0nOenurEF51TVDPAD4OXAQ8B1STYm2Q68DthaVT8G/hh4hNkf/JPAJ1ZxH5KkZVrvD4HvYjYYQ+B24EHgJ0l+htkAvBZ4BbNvAb13oQskuSXJMMlwenp6nZcrSe3oE4DTwNY5z7d0YwvO6d7fvxR4tqpmquo9VfWaqtoDvAx4AngNQFV9q6oKuBf4zYW+eVUdqqpBVQ0mJiaWcWuSpMX0CcBRYEeS7Uk2AXuBw/PmHAb2dV9fDzxQVZXkJUkuAUhyDTBTVY8yG4zJJC/8RL8GeGyV9yJJWoaNS02oqpkk+4EjwAbgrqo6keQgMKyqw8y+f/+ZJFPAGWYjAXAlcCTJWWZ/6L+zu+Z3k7wf+EqSHwPfBm5c21uTJC0ms+/AjIfBYFDD4XDUy5CksZLkWFUN5o/7m8CS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1CgDIEmNMgCS1KheAUiyO8njSaaS3LrA8Rcnuac7/rUk27rxTUk+meSRJA8luXrOOZuSHEryRJJvJnnbGt2TJKmHjUtNSLIBuBO4BjgFHE1yuKoenTPtJuC5qnplkr3Ah4F3ADcDVNWrk1wJfCnJVVV1Fngf8ExVvSrJi4DL1/TOJEmL6vMKYCcwVVUnq+p54G5gz7w5e4BPd1/fB7wxSYBJ4AGAqnoG+D4w6Ob9IfDB7tjZqvream5EkrQ8fQKwGXh6zvNT3diCc6pqBvgB8HLgIeC6JBuTbAdeB2xN8rLuvNuSHE/y+SQ/t9A3T3JLkmGS4fT0dO8bkyQtbr0/BL6L2WAMgduBB4GfMPvW0xbgwar6NeCrwEcWukBVHaqqQVUNJiYm1nm5ktSOJT8DAE4DW+c839KNLTTnVJKNwKXAs1VVwHtemJTkQeAJ4FngR8AXukOfZ/ZzBEnSedLnFcBRYEeS7Uk2AXuBw/PmHAb2dV9fDzxQVZXkJUkuAUhyDTBTVY92YfgicHV3zhuBR5EknTdLvgKoqpkk+4EjwAbgrqo6keQgMKyqw8AngM8kmQLOMBsJgCuBI0nOMvsq4Z1zLv0X3Tm3A9PAu9bqpiRJS8vsX8bHw2AwqOFwOOplSNJYSXKsqgbzx/1NYElqlAGQpEYZAElqlAGQpEYZAElqlAGQpEYZAElqlAGQpEYZAElqlAGQpEYZAElqlAGQpEYZAElqlAGQpEYZAElqlAGQpEYZAElqlAGQpEYZAElqlAGQpEYZAElqlAGQpEYZAElqlAGQpEalqka9ht6STAPfnjN0KfCDBaYuNH4F8L11WtpynGvN5/t6yzmvz9zF5qzkmHu4tued7z081/wLYQ9b3L+fr6qJ/zdaVWP7AA71HQeGo17vYms+39dbznl95i42ZyXH3MPx3sNF9nXke9ji/p3rMe5vAX1xmeMXgrVe20qvt5zz+sxdbM5KjrmHa3ve+d5D929tz1vt/i1orN4CWo0kw6oajHodWjn3cPy5hxeWcX8FsByHRr0ArZp7OP7cwwtIM68AJEk/raVXAJKkOQyAJDXKAEhSowwAkOQtSf42yT1Jdo16PVqeJL+Q5BNJ7hv1WtRfkkuSfLr7s/f7o15Pi8Y+AEnuSvJMkm/MG9+d5PEkU0luXewaVfWPVXUz8G7gHeu5Xv20Ndq/k1V10/quVH0scz/fCtzX/dm77rwvVuMfAOBTwO65A0k2AHcC1wKTwA1JJpO8Osk/z3tcOefUv+zO0/nzKdZu/zR6n6LnfgJbgKe7aT85j2tUZ+OoF7BaVfWVJNvmDe8EpqrqJECSu4E9VfVB4PfmXyNJgA8BX6qq4+u7Ys21FvunC8dy9hM4xWwEvs7F8ZfRsXOx/kPfzP/9zQJm/0XbvMj8PwV+F7g+ybvXc2HqZVn7l+TlSf4GeG2S96734rRs59rPLwBvS/IxLuz/dcRFa+xfAayFqroDuGPU69DKVNWzzH5+ozFSVf8DvGvU62jZxfoK4DSwdc7zLd2YxoP7d3FxPy9QF2sAjgI7kmxPsgnYCxwe8ZrUn/t3cXE/L1BjH4AknwO+CvxSklNJbqqqGWA/cAR4DLi3qk6Mcp1amPt3cXE/x4v/MzhJatTYvwKQJK2MAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWqUAZCkRhkASWrU/wLKtRz4RyU5EQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate the Sedov Problem\n",
    "\n",
    "rArr_np = np.arange(0.01,5.,0.01)\n",
    "rbarArr_np = rArr_np\n",
    "rhoArr_np = np.ones(rArr_np.size)*0.1\n",
    "mArr_np = 4.*np.pi/3.*rArr_np**3*rhoArr_np\n",
    "PArr_np = rhoArr_np*1e-6\n",
    "PArr_np[rArr_np < 0.5] = 1e-2\n",
    "phiArr_np = np.ones(rArr_np.size)\n",
    "confFactor_np = rArr_np/rbarArr_np\n",
    "if sys.version_info[0] < 3:\n",
    "    np.savetxt( \"sedov.txt\", zip(rArr_np,rhoArr_np,PArr_np,mArr_np,phiArr_np,confFactor_np,rbarArr_np), fmt=\"%.15e\")\n",
    "else:\n",
    "    np.savetxt( \"sedov.txt\", list(zip(rArr_np,rhoArr_np,PArr_np,mArr_np,phiArr_np,confFactor_np,rbarArr_np)), fmt=\"%.15e\")\n",
    "\n",
    "pl.semilogx(rArr_np, rhoArr_np)\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Convert metric to be in terms of ADM quantities\n",
    "\n",
    "Above, the line element was written:\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 terms of $G=c=1$ units adopted by NRPy+, this becomes:\n",
    "$$\n",
    "ds^2 = - e^\\nu dt^2 + \\left(1 - \\frac{2M}{r}\\right)^{-1} dr^2 + r^2 d\\Omega^2.\n",
    "$$\n",
    "\n",
    "The ADM 3+1 line element for this diagonal metric in spherical coordinates is given by:\n",
    "$$\n",
    "ds^2 = (-\\alpha^2 + \\beta_k \\beta^k) dt^2 + \\gamma_{rr} dr^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/2} \\\\\n",
    "\\beta^k &= 0 \\\\\n",
    "\\gamma_{rr} &= \\left(1 - \\frac{2M}{r}\\right)^{-1}\\\\\n",
    "\\gamma_{\\theta\\theta} &= r^2 \\\\\n",
    "\\gamma_{\\phi\\phi} &= r^2 \\sin^2 \\theta \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert to Cartesian coordinates\n",
    "\n",
    "The above metric is given in spherical coordinates and we need everything in Cartesian coordinates.  Given this the \n",
    "transformation to Cartesian coordinates is \n",
    "\\begin{equation}\n",
    "g_{\\mu\\nu} = \\Lambda^{\\mu'}_{\\mu} \\Lambda^{\\nu'}_{\\nu} g_{\\mu'\\nu'},\n",
    "\\end{equation}\n",
    "where $\\Lambda^{\\mu'}_{\\mu}$ is the Jacobian defined as\n",
    "\\begin{equation}\n",
    "\\Lambda^{\\mu'}_{\\mu} = \\frac{\\partial x'^{\\mu'}}{\\partial x^{\\mu}}\n",
    "\\end{equation}\n",
    "In this particular case $x'$ is in spherical coordinates and $x$ is in Cartesian coordinates.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"NRPY+spherical_to_Cartesian_metric.h\"\n"
     ]
    }
   ],
   "source": [
    "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",
    "from outputC import outputC         # NRPy+: Basic C code output functionality\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\n",
    "\n",
    "# The ADM & BSSN formalisms only work in 3D; they are 3+1 decompositions of Einstein's equations.\n",
    "#    To implement axisymmetry or spherical symmetry, simply set all spatial derivatives in\n",
    "#    the relevant angular directions to zero; DO NOT SET DIM TO ANYTHING BUT 3.\n",
    "# Step 0: Set spatial dimension (must be 3 for BSSN)\n",
    "DIM = 3\n",
    "\n",
    "# Set the desired *output* coordinate system to Cylindrical:\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n",
    "rfm.reference_metric()\n",
    "\n",
    "CoordType_in = \"Spherical\"\n",
    "\n",
    "r_th_ph_or_Cart_xyz_of_xx = []\n",
    "if CoordType_in == \"Spherical\":\n",
    "    r_th_ph_or_Cart_xyz_of_xx = rfm.xxSph\n",
    "elif CoordType_in == \"Cartesian\":\n",
    "    r_th_ph_or_Cart_xyz_of_xx = rfm.xxCart\n",
    "\n",
    "Jac_dUSphorCart_dDrfmUD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        Jac_dUSphorCart_dDrfmUD[i][j] = sp.diff(r_th_ph_or_Cart_xyz_of_xx[i],rfm.xx[j])\n",
    "\n",
    "Jac_dUrfm_dDSphorCartUD, dummyDET = ixp.generic_matrix_inverter3x3(Jac_dUSphorCart_dDrfmUD)\n",
    "\n",
    "betaU   = ixp.zerorank1()\n",
    "gammaDD = ixp.zerorank2()\n",
    "gammaSphDD = ixp.zerorank2()\n",
    "grr, gthth, gphph = sp.symbols(\"grr gthth gphph\")\n",
    "gammaSphDD[0][0] = grr\n",
    "gammaSphDD[1][1] = gthth\n",
    "gammaSphDD[2][2] = gphph\n",
    "betaSphU = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        betaU[i] += Jac_dUrfm_dDSphorCartUD[i][j] * betaSphU[j]\n",
    "        for k in range(DIM):\n",
    "            for l in range(DIM):\n",
    "                gammaDD[i][j] += Jac_dUSphorCart_dDrfmUD[k][i]*Jac_dUSphorCart_dDrfmUD[l][j] * gammaSphDD[k][l]\n",
    "\n",
    "outputC([gammaDD[0][0], gammaDD[0][1], gammaDD[0][2], gammaDD[1][1], gammaDD[1][2], gammaDD[2][2]],\n",
    "              [\"mi.gamDDxx\", \"mi.gamDDxy\", \"mi.gamDDxz\", \"mi.gamDDyy\", \"mi.gamDDyz\",\"mi.gamDDzz\"], filename=\"NRPY+spherical_to_Cartesian_metric.h\")"
   ]
  }
 ],
 "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
}