{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"# Code and Physical Units\n",
"\n",
"## GRHD Units in terms of EOS \n",
"\n",
"## This notebook introduces an equation of state for general relativistic hydrodynamics (GRHD) and demonstrates the conversion of dimensional quantities to dimensionless units.\n",
"\n",
"**Notebook Status:** Not Validated \n",
"\n",
"\n",
"## Introduction\n",
"$\\newcommand{\\rhoCode}{{\\tilde{\\rho}}}$\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 fiducial $\\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 defined 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": [
"\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"The module is organized as follows:\n",
"\n",
"1. [Step 1](#comments): Zach's comments\n",
"1. [Step 2](#tov_solver): TOV solver as illustration\n",
" 1. [Step 2.a](#prescription): Another dimensionless prescription\n",
" 1. [Step 2.b](#metric): Metric for TOV equation\n",
"1. [Step 3](#adm_conversion): Convert metric to be in terms of ADM quantities\n",
"1. [Step 4](#cartesian_conversion): Convert to Cartesian coordinates\n",
"1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 1: Zach's comments \\[Back to [top](#toc)\\]\n",
"$$\\label{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": [
"\n",
"\n",
"# Step 2: TOV Solver as illustration \\[Back to [top](#toc)\\]\n",
"$$\\label{tov_solver}$$\n",
"\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 measured 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. A lot 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": [
"\n",
"\n",
"## Step 2.a: Another dimensionless prescription \\[Back to [top](#toc)\\]\n",
"$$\\label{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": [
"\n",
"\n",
"## Step 2.b: Metric for TOV equation \\[Back to [top](#toc)\\]\n",
"$$\\label{metric}$$\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",
"Let's 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": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:39.127986Z",
"iopub.status.busy": "2021-03-07T17:11:39.126708Z",
"iopub.status.idle": "2021-03-07T17:11:39.830488Z",
"shell.execute_reply": "2021-03-07T17:11:39.831151Z"
},
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD6CAYAAACxrrxPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3hVZb728e8vtIwIMkAQhmYZRBBFINKLBTkUBaVjVBRCKHrEYZDhHc4gyhwRHcSDIqGpBJQioiJSwpBIiwIBKcbIGFQUkEFQAigtyfP+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": [
""
]
},
"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": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:39.839347Z",
"iopub.status.busy": "2021-03-07T17:11:39.838663Z",
"iopub.status.idle": "2021-03-07T17:11:40.370639Z",
"shell.execute_reply": "2021-03-07T17:11:40.371338Z"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAOxUlEQVR4nO3cf6jdd33H8efLZBHssLb2dmASlmzGjTscOk/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": [
""
]
},
"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,
"jupyter": {
"outputs_hidden": true
}
},
"source": [
"\n",
"\n",
"# Step 3: Convert metric to be in terms of ADM quantities \\[Back to [top](#toc)\\]\n",
"$$\\label{adm_conversion}$$\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": [
"\n",
"\n",
"# Step 4: Convert to Cartesian coordinates \\[Back to [top](#toc)\\]\n",
"$$\\label{cartesian_conversion}$$\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": {
"execution": {
"iopub.execute_input": "2021-03-07T17:11:40.377788Z",
"iopub.status.busy": "2021-03-07T17:11:40.376714Z",
"iopub.status.idle": "2021-03-07T17:11:40.975660Z",
"shell.execute_reply": "2021-03-07T17:11:40.975000Z"
}
},
"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.xx_to_Cart\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\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 5: 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-GRHD_UnitConversion.pdf](Tutorial-GRHD_UnitConversion.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": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-GRHD_UnitConversion.tex, and compiled LaTeX file to PDF\n",
" file Tutorial-GRHD_UnitConversion.pdf\n"
]
}
],
"source": [
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-GRHD_UnitConversion\")"
]
}
],
"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
}