{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\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": [
""
]
},
"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": [
""
]
},
"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
}