examples/BHEvolutionWithSF.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ff3f7338",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<script>Jupyter.notebook.kernel.restart()</script>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Isotropic Schwarzschild BH example\n",
    "# see further details in https://github.com/GRChombo/engrenage/wiki/Running-the-black-hole-example\n",
    "\n",
    "# restart the kernel to clear past work\n",
    "# (can also do this manually from the Kernel options above)\n",
    "from IPython.core.display import HTML\n",
    "HTML(\"<script>Jupyter.notebook.kernel.restart()</script>\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5ae36db4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the required python modules\n",
    "import numpy as np\n",
    "from scipy.interpolate import interp1d\n",
    "from scipy.integrate import odeint\n",
    "from scipy.integrate import solve_ivp\n",
    "import time\n",
    "import random\n",
    "import sys\n",
    "from tqdm import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# homemade code\n",
    "sys.path.append('../')\n",
    "from source.rhsevolution import *              # go here to look at how the evolution works\n",
    "from source.bhinitialconditionswithSF import *       # go here to change the initial conditions\n",
    "from source.hamdiagnostic import *             # go here to change the Ham constraint diagnostic\n",
    "from source.fluxdiagnostic import *             # go here to change the Flux constraint diagnostic\n",
    "from source.Grid import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e0b8dc67",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiAAAAHACAYAAACbJdMDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqDUlEQVR4nO3df3RU9Z3/8ddkMpkwMQlCSEIgBGjxICA/I66gNSDBqljYH2xXtER7ds/RRQRzWgXrtwaLhNpTy3FVKtiq/UGhFrHUFWRQAX9UYAMoiougqMgPY7aYiYxMhuTz/cPNrCEhZJI7n5sfz4eHI/O5n8/nfvLO6LzOvXfu9RhjjAAAACxKcnsBAACg+yGAAAAA6wggAADAOgIIAACwjgACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOtcCyDbtm3Tddddp7y8PHk8Hj377LMJ3+eRI0d04403qnfv3goEAho9erQqKioSvl8AANCYawHk5MmTGjVqlB5++GEr+ztx4oQmTpwon8+nDRs2aN++ffr5z3+unj17Wtk/AAD4P56O8DA6j8ejdevWacaMGbG22tpa3XPPPfr973+vzz//XCNGjNBPf/pTFRUVtWkfCxYs0GuvvaZXXnnFmUUDAIA267DXgNx888167bXXtHr1ar311luaOXOmvv3tb+vAgQNtmm/9+vUqLCzUzJkzlZ2drTFjxmjlypUOrxoAALRGhzwC8v7772vIkCH65JNPlJeXF+s3ZcoUjR8/XkuWLIl7H6mpqZKk0tJSzZw5Uzt27ND8+fP12GOPafbs2Y78HAAAoHWS3V5Ac3bt2iVjjC644IJG7ZFIRL1795Ykffjhhxo0aFCL88yZMyd2jUl9fb0KCwtj4WXMmDF65513tHz5cgIIAACWdcgAUl9fL6/Xq4qKCnm93kbbzjvvPElSv3799O6777Y4z/nnnx/7e9++fTVs2LBG2y+88EKtXbvWoVUDAIDW6pABZMyYMaqrq1NlZaUuv/zyZvv4fD4NHTq01XNOnDhR+/fvb9T23nvvqaCgoF1rBQAA8XMtgHzxxRc6ePBg7PWhQ4e0Z88e9erVSxdccIFuuOEGzZ49Wz//+c81ZswYVVVV6aWXXtJFF12ka665Ju793XHHHZowYYKWLFmif/7nf9aOHTu0YsUKrVixwskfCwAAtIJrF6Fu2bJFkyZNatJeUlKiJ598UtFoVIsXL9ZvfvMbHTlyRL1799all16qRYsW6aKLLmrTPp977jktXLhQBw4c0KBBg1RaWqp/+7d/a++PAgAA4tQhvgUDAAC6lw57HxAAANB1EUAAAIB11i9Cra+v19GjR5Weni6Px2N79wAAoA2MMaqpqVFeXp6Sktp//MJ6ADl69Kjy8/Nt7xYAADjg8OHD6t+/f7vnsR5A0tPTJX31A2RkZDg2bzQa1aZNmzR16lT5fD7H5kXLqLs7qLs7qLs7qLs7zqx7KBRSfn5+7HO8vawHkIbTLhkZGY4HkEAgoIyMDN6gFlF3d1B3d1B3d1B3d5yt7k5dPsFFqAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOusP4wuIYyRak/KWxeRak9KhocVWRONUnc3UHd3UHd3UPf4+QKSQw+NS5SuEUCiYfl+VqBpkvSW24vpXnwSdXcBdXcHdXcHdW+Du49KKWlur6JFnIIBAADWdY0jIL6Aoj/8SC+8sElXXTVVPh+H6GyJRqPU3QXU3R3U3R3UvQ18AbdXcE5dI4B4PFJKmuq8/q8OOfEGtccTpe5uoO7uoO7uoO5dUtynYI4cOaIbb7xRvXv3ViAQ0OjRo1VRUZGItQEAgC4qriMgJ06c0MSJEzVp0iRt2LBB2dnZev/999WzZ88ELQ8AAHRFcQWQn/70p8rPz9cTTzwRaxs4cKDTawIAAF1cXKdg1q9fr8LCQs2cOVPZ2dkaM2aMVq5cmai1AQCALiquIyAffPCBli9frtLSUt19993asWOHbr/9dvn9fs2ePbvZMZFIRJFIJPY6FApJ+uqq5mg02o6lN9Ywl5Nz4tyouzuouzuouzuouzvOrLvT9fcYY0xrO6ekpKiwsFCvv/56rO3222/Xzp079de//rXZMWVlZVq0aFGT9lWrVikQ6PhfEwIAAFI4HNasWbNUXV2tjIyMds8X1xGQvn37atiwYY3aLrzwQq1du/asYxYuXKjS0tLY61AopPz8fE2dOtWRH6BBNBpVMBhUcXEx3xO3iLq7g7q7g7q7g7q748y6N5zBcEpcAWTixInav39/o7b33ntPBQUFZx3j9/vl9/ubtPt8voS8kRI1L1pG3d1B3d1B3d1B3d3RUHenax/XRah33HGH3njjDS1ZskQHDx7UqlWrtGLFCs2ZM8fRRQEAgK4trgBy8cUXa926dfrDH/6gESNG6Cc/+YmWLVumG264IVHrAwAAXVDct2KfNm2apk2bloi1AACAboKn4QIAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsiyuAlJWVyePxNPqTm5ubqLUBAIAuKjneAcOHD9fmzZtjr71er6MLAgAAXV/cASQ5OZmjHgAAoF3iDiAHDhxQXl6e/H6/LrnkEi1ZskSDBw8+a/9IJKJIJBJ7HQqFJEnRaFTRaLQNS25ew1xOzolzo+7uoO7uoO7uoO7uOLPuTtffY4wxre28YcMGhcNhXXDBBfr000+1ePFi/fd//7feeecd9e7du9kxZWVlWrRoUZP2VatWKRAItH3lAADAmnA4rFmzZqm6uloZGRntni+uAHKmkydP6hvf+IbuvPNOlZaWNtunuSMg+fn5qqqqcuQHaBCNRhUMBlVcXCyfz+fYvGgZdXcHdXcHdXcHdXfHmXUPhULKyspyLIDEfQrm69LS0nTRRRfpwIEDZ+3j9/vl9/ubtPt8voS8kRI1L1pG3d1B3d1B3d1B3d3RUHena9+u+4BEIhG9++676tu3r1PrAQAA3UBcAeQHP/iBtm7dqkOHDmn79u36p3/6J4VCIZWUlCRqfQAAoAuK6xTMJ598ouuvv15VVVXq06eP/u7v/k5vvPGGCgoKErU+AADQBcUVQFavXp2odQAAgG6EZ8EAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCAAAsI4AAgAArGtXACkvL5fH49H8+fMdWg4AAOgO2hxAdu7cqRUrVmjkyJFOrgcAAHQDbQogX3zxhW644QatXLlS559/vtNrAgAAXVxyWwbNmTNH1157raZMmaLFixe32DcSiSgSicReh0IhSVI0GlU0Gm3L7pvVMJeTc+LcqLs7qLs7qLs7qLs7zqy70/WPO4CsXr1au3bt0s6dO1vVv7y8XIsWLWrSvmnTJgUCgXh3f07BYNDxOXFu1N0d1N0d1N0d1N0dDXUPh8OOzusxxpjWdj58+LAKCwu1adMmjRo1SpJUVFSk0aNHa9myZc2Oae4ISH5+vqqqqpSRkdG+1X9NNBpVMBhUcXGxfD6fY/OiZdTdHdTdHdTdHdTdHWfWPRQKKSsrS9XV1Y58fsd1BKSiokKVlZUaN25crK2urk7btm3Tww8/rEgkIq/X22iM3++X3+9vMpfP50vIGylR86Jl1N0d1N0d1N0d1N0dDXV3uvZxBZArr7xSe/fubdR28803a+jQobrrrruahA8AAIDmxBVA0tPTNWLEiEZtaWlp6t27d5N2AACAs+FOqAAAwLo2fQ3367Zs2eLAMgAAQHfCERAAAGAdAQQAAFhHAAEAANYRQAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFAAACAdQQQAABgXbvvhAoAQGdVX1+v2tpat5fRIfh8PqsPlSWAAAC6pdraWh06dEj19fVuL6XD6Nmzp3Jzc+XxeBK+LwIIAKDbMcbo2LFj8nq9ys/PV1JS974iwRijcDisyspKSVLfvn0Tvk8CCACg2zl9+rTC4bDy8vIUCATcXk6H0KNHD0lSZWWlsrOzE76/7h35AADdUl1dnSQpJSXF5ZV0LA1hLBqNJnxfBBAAQLdl41qHzsRmPQggAADAOgIIAACwjgACAACsI4AAAADrCCAAAHQSAwcO1LJlyxq1jR49WmVlZa6spz24DwgAoNszxujLaJ0r++7h83bLb+MQQAAA3d6X0ToN+/ELrux7331XKZDS/T6OOQUDAACs636RCwCAM/TwebXvvqtc23drJSUlyRjTqM3GXUsTgQACAOj2PB5PpzgN0qdPHx07diz2OhQK6dChQy6uqO04BQMAQCcxefJk/fa3v9Urr7yit99+WyUlJfJ6W38EpSPp+HEPAABIkhYuXKgPPvhA06ZNU2Zmpn7yk5902iMgBBAAADqJjIwMrVmzplFbSUmJS6tpH07BAAAA6wggAADAOgIIAACwjgACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCAAAncRNN92kGTNmNGr705/+pNTUVD3wwAPuLKqNeBYMAACd1OOPP645c+bokUce0b/+67+6vZy4cAQEAIBO6IEHHtBtt92mVatWdbrwIXEEBAAAyRgpGnZn376A5PHENWTBggV65JFH9Nxzz2nKlCkJWlhiEUAAAIiGpSV57uz77qNSSlqru2/YsEF//vOf9eKLL2ry5MkJXFhixXUKZvny5Ro5cqQyMjKUkZGhSy+9VBs2bEjU2gAAwBlGjhypgQMH6sc//rFqamrcXk6bxXUEpH///lq6dKm++c1vSpKeeuopTZ8+Xbt379bw4cMTskAAABLOF/jqSIRb+45Dv379tHbtWk2aNEnf/va3tXHjRqWnpydocYkTVwC57rrrGr2+//77tXz5cr3xxhsEEABA5+XxxHUaxG0DBgzQ1q1bNWnSJE2dOlUvvPCCMjIy3F5WXNr8LZi6ujqtXr1aJ0+e1KWXXurkmgAAwDn0799fW7Zs0f/8z/9o6tSpqq6udntJcYn7ItS9e/fq0ksv1alTp3Teeedp3bp1GjZs2Fn7RyIRRSKR2OtQKCRJikajikajbVhy8xrmcnJOnBt1dwd1dwd1d0ci6h6NRmWMUX19verr6x2bN9GMMbF1S1Lfvn318ssv68orr1RxcbE2btyonj17tnn++vp6GWMUjUZj+0jU+95jjDHxDKitrdXHH3+szz//XGvXrtXjjz+urVu3njWElJWVadGiRU3aV61apUAgvvNeAAA4ITk5Wbm5ucrPz1dKSorby+kwamtrdfjwYR0/flynT59utC0cDmvWrFmqrq525HRP3AHkTFOmTNE3vvENPfbYY81ub+4ISH5+vqqqqhw9XxWNRhUMBlVcXCyfz+fYvGgZdXcHdXcHdXdHIup+6tQpHT58WAMHDlRqaqojc3YFp06d0ocffqj8/Hx5vd5GdQ+FQsrKynIsgLT7PiDGmEYB40x+v19+v79Ju8/nS8h/wImaFy2j7u6g7u6g7u5wsu51dXXyeDxKSkpSUhI3BW+QlJQkj8cjn88nr9cr6f/q7vR7Pq4Acvfdd+vqq69Wfn6+ampqtHr1am3ZskUbN250dFEAAKBriyuAfPrpp/re976nY8eOKTMzUyNHjtTGjRtVXFycqPUBAIAuKK4A8qtf/SpR6wAAAN0IJ74AAIB1BBAAAGAdAQQAAFhHAAEAANYRQAAA6CSKioo0f/58t5fhiHbfiAwAgO6qrr5Ouyp36bPwZ+oT6KOx2WPlTfK6vaxOgQACAEAbbP5os5buWKpPw5/G2nICOVowfoGmFExxcWX/JxqNNrmDaXNtbuAUDAAAcdr80WaVbiltFD4kqTJcqdItpdr80eaE7bu+vl533nmnevXqpdzcXJWVlcW2eTwe/fKXv9T06dOVlpamxYsXq6ysTKNHj9avf/1rDR48WH6/X+18DJwjCCAAgG7PGKNwNNyqPzWRGpXvKJdR0w9x87//LN2xVDWRmlbNF28YeOqpp5SWlqbt27frgQce0H333adgMBjbfu+992r69Onau3evvv/970uSDh48qD/+8Y9au3at9uzZ065aOYVTMACAbu/L01/qklWXODbfp+FPNWH1hFb13T5ruwK+QKvnHjlypO69915J0pAhQ/Twww/rxRdfjD0WZdasWbHg0aC2tla//e1v1adPn1bvJ9E4AgIAQCcycuTIRq/79u2rysrK2OvCwsImYwoKCjpU+JA4AgIAgHok99D2Wdtb1bfi0wr9+4v/fs5+j175qMbljGvVvuNx5gWkHo9H9fX1sddpaWlNxjTX5jYCCACg2/N4PK0+DTIhb4JyAjmqDFc2ex2IRx7lBHI0IW8CX8ltAadgAACIgzfJqwXjF0j6Kmx8XcPru8bfRfg4BwIIAABxmlIwRQ8WPajsQHaj9pxAjh4serDD3AekI+MUDAAAbTClYIom5U+yeifULVu2NGl79tlnY39v7iu9ZWVlje4V0lEQQAAAaCNvklcX517s9jI6JU7BAAAA6wggAADAOgIIAACwjgACAACsI4AAALqtjvBU2I7k63dUTTS+BQMA6HZ8Pp88Ho8+++wz9enTRx6P59yDujBjjGpra/XZZ58pKSlJKSkpqqurS+g+CSAAgG7H6/Wqf//++uSTT/Thhx+6vZwOIxAIaMCAAUpKSiKAAACQCOedd56GDBmiaDTq9lI6BK/Xq+TkZGtHgwggAIBuy+v1yuvlmS1u4CJUAABgHQEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFhHAAEAANYRQAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFgXVwApLy/XxRdfrPT0dGVnZ2vGjBnav39/otYGAAC6qLgCyNatWzVnzhy98cYbCgaDOn36tKZOnaqTJ08man0AAKALSo6n88aNGxu9fuKJJ5Sdna2Kigp961vfcnRhAACg64orgJypurpaktSrV6+z9olEIopEIrHXoVBIkhSNRhWNRtuz+0Ya5nJyTpwbdXcHdXcHdXcHdXfHmXV3uv4eY4xpy0BjjKZPn64TJ07olVdeOWu/srIyLVq0qEn7qlWrFAgE2rJrAABgWTgc1qxZs1RdXa2MjIx2z9fmADJnzhz953/+p1599VX179//rP2aOwKSn5+vqqoqR36ABtFoVMFgUMXFxfL5fI7Ni5ZRd3dQd3dQd3dQd3ecWfdQKKSsrCzHAkibTsHMnTtX69ev17Zt21oMH5Lk9/vl9/ubtPt8voS8kRI1L1pG3d1B3d1B3d1B3d3RUHenax9XADHGaO7cuVq3bp22bNmiQYMGOboYAADQPcQVQObMmaNVq1bpz3/+s9LT03X8+HFJUmZmpnr06JGQBQIAgK4nrvuALF++XNXV1SoqKlLfvn1jf9asWZOo9QEAgC4o7lMwAAAA7cWzYAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFhHAAEAANYRQAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFhHAAEAANYRQAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFhHAAEAANYRQAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWxR1Atm3bpuuuu055eXnyeDx69tlnE7AsAADQlcUdQE6ePKlRo0bp4YcfTsR6AABAN5Ac74Crr75aV199dSLW0mbGGIVrTytSJ4VrT8tnPG4vqduIRqm7G6i7O6i7O6h7/Hr4vPJ4Onat4g4g8YpEIopEIrHXoVBIkhSNRhWNRh3ZR7j2tEb95CVJybpzx0uOzIl4UHd3UHd3UHd3UPd4vPn/JiuQ0r6P+IbP6DP/7ZSEB5Dy8nItWrSoSfumTZsUCAQc2UekTrLwowAA0Cm88MIm+b3OzBUMBiVJ4XDYmQn/l8cYY9o82OPRunXrNGPGjLP2ae4ISH5+vqqqqpSRkdHWXTdijFEoHNFLL72kyZMny+cjjNgSjZ6m7i6g7u6g7u6g7vFz4hRMNBpVMBhUcXGxfD6fQqGQsrKyVF1d7cjnd8J/k36/X36/v0m7z+eTz+dzbD+ZHo/8XikzLdXRedGyaDRK3V1A3d1B3d1B3d3V8HntdO25DwgAALAu7iMgX3zxhQ4ePBh7fejQIe3Zs0e9evXSgAEDHF0cAADomuIOIP/1X/+lSZMmxV6XlpZKkkpKSvTkk086tjAAANB1xR1AioqK1I7rVgEAALgGBAAA2EcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFhHAAEAANYRQAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFhHAAEAANYRQAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFhHAAEAANYRQAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWEUAAAIB1yW0Z9Oijj+pnP/uZjh07puHDh2vZsmW6/PLLnV5bq31Z+6UWv7ZY6z9fr3tW3dNom1de+Tw+JSclK8mTJG+SV6neVKV6U+VP9suX5FO0PqoUb4pSvamSpEhdRKnJqeqV2qvZ/Xk8HuWm5aqnv6eqI9U6Hj4uY0ysraa2RsZj1DOlp85PPV8nIidUfaq6UVt1bbUyUzLb9O/z/eerT6CPjDGq+rJKJyInGrX97dTf1LtHb0f+3ifQR6OyRunNqjf1WfizJq8zfZl6P/q+Nn64UbnpuS327RPoo7HZYyVJuyp3Od5nbPZYeZO8qquva9KeqL7eJK8knXVbW8aca1vD9g+iH8TqHu/YlrbbmqO1fRLRr619j9cc1wfRD1RXXyeffO2e10Z/m2PcGNfesR1hfEeaw6a4A8iaNWs0f/58Pfroo5o4caIee+wxXX311dq3b58GDBiQiDW26PaXbtfLh18+6/Y61anO1El1FhfVxSR5klRv6s/6WpL0euv6ZqZkSh6pOlLteJ+cQI6uGXSNnj/0vD4Nf5rwvjmBHC0Yv0CStHTH0ibb2jLmXNumFEzR5o82q3xHuSrDlbG6xzO2pe2SztnHiTla2ycR/Zzo+9z657Rw/MJ2zWujv80xboxr79iOML4jzWGbxxhj4hlwySWXaOzYsVq+fHms7cILL9SMGTNUXl5+zvGhUEiZmZmqrq5WRkZG/Cv+mnOFDyCRPPLIKK7/fFocc65tknTT8Jv05DtPNunX3rGS9GDRg5Kk0i2lZ+3jxByt7dMQdpzsJ6lD9LXR3+aYRI+7Iu8KPf/887rmmmvk8/niGnuu4OPm+I40R3Oi0Wijujv5+S3FGUBqa2sVCAT09NNP6+///u9j7fPmzdOePXu0devWc87h1A/wZe2XGv+H8W0eD3RGzR59cmhsdo9sSVLll5UJnaM1fXICOfrTtD/pH5/7x6+O9DjQb9131kmSZqyf4WrfhlNziewvydoYG+OevuZpbQ5u1lVXXaXk5OR27bO9a3ZqvK05PPIoJ5Cjjf+4Me7TMR0qgBw9elT9+vXTa6+9pgkTJsTalyxZoqeeekr79+9vMiYSiSgSicReh0Ih5efnq6qqql0/wNKdS/XHA39s83gAALqLFVeuUGFOYVxjotGogsGgiouLYwEkKyvLsQDSpotQPR5Po9fGmCZtDcrLy7Vo0aIm7Zs2bVIgEGjL7iVJFTUVbR4LAEB3EvxrUJUpZz/S0uLYYFCSFA6HnVxSfAEkKytLXq9Xx48fb9ReWVmpnJycZscsXLhQpaWlsdcNR0CmTp3argT11s639P6B99s8HkDHdvuo2/XQmw851u8/iv5DkjR3y1xX+47NHqtdlbsS2l+StTE2xv3isl/oxN4Tmjx5cuwakLbus71rdmq87TmKLy125AiIk+IKICkpKRo3bpyCwWCja0CCwaCmT5/e7Bi/3y+/39+k3efzNbqYKF4/KPwBp2DQ7SR5kmSMifvi13ON9cjz1bUZHqkyXHnW+ds7R2v75ARyVDKiRGsOrHGs3+X5X90qICeQ42pfb5JXl+dfntD+kqyNsTHusv6X6YW3X1BGj4zY50Zb99neNTs13vYc4/PGt/kruQ2f1+35zG5O3DciKy0t1eOPP65f//rXevfdd3XHHXfo448/1i233OLows6lR0oPTcqfZHWfwNc1XGF+5t/bOuZc2zzyqGRYScLGLrhkQeyrvImao7V97hp/l1KSUxzt503yypvkdb2vpIT3tznGjXHtHdsRxnekOdwSdwD57ne/q2XLlum+++7T6NGjtW3bNj3//PMqKChIxPpa9NDkhwghFiR5klp8HU/fnik9lenPTEif3ECubh5+s3ICOVb65gRy9IuiX+gXRb9QdiC73WPOte3BogdVWliqB4seVJ9AnzaPPdv2KQVTNKVgSot9nJijtX0kOd7Pqb7Zgex2zWujv80xboxr79iOML4jzeGGuO8D0l5Of41H+tqdUD9e32Qbd0JN/J1Qt2/friGjhnAnVIt3Qj0VOaXlf1mub476JndCbUO/tvY9XnNcB988qFuvu1Wp/tR2z2ujv80xiRp35tdBndpnRxnfkeb4ug71NVwnJCKASOd+gyIxqLs7qLs7qLs7qLs7Eh1AeBgdAACwjgACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCAAAsC6up+E6oeHGq04/1jcajSocDisUCnGnPIuouzuouzuouzuouzvOrHvD57ZTN1C3HkBqamokSfn5+bZ3DQAA2qmmpkaZmZnn7ngO1p8FU19fr6NHjyo9PV0eT+seYd4aoVBI+fn5Onz4sKPPmEHLqLs7qLs7qLs7qLs7zqy7MUY1NTXKy8tTUlL7r+CwfgQkKSlJ/fv3T9j8GRkZvEFdQN3dQd3dQd3dQd3d8fW6O3HkowEXoQIAAOsIIAAAwLouE0D8fr/uvfde+f1+t5fSrVB3d1B3d1B3d1B3dyS67tYvQgUAAOgyR0AAAEDnQQABAADWEUAAAIB1XSaAPProoxo0aJBSU1M1btw4vfLKK24vqcsoLy/XxRdfrPT0dGVnZ2vGjBnav39/oz7GGJWVlSkvL089evRQUVGR3nnnHZdW3DWVl5fL4/Fo/vz5sTbqnhhHjhzRjTfeqN69eysQCGj06NGqqKiIbafuzjt9+rTuueceDRo0SD169NDgwYN13333qb6+PtaHurfftm3bdN111ykvL08ej0fPPvtso+2tqXEkEtHcuXOVlZWltLQ0fec739Enn3wS/2JMF7B69Wrj8/nMypUrzb59+8y8efNMWlqa+eijj9xeWpdw1VVXmSeeeMK8/fbbZs+ePebaa681AwYMMF988UWsz9KlS016erpZu3at2bt3r/nud79r+vbta0KhkIsr7zp27NhhBg4caEaOHGnmzZsXa6fuzvvb3/5mCgoKzE033WS2b99uDh06ZDZv3mwOHjwY60Pdnbd48WLTu3dv89xzz5lDhw6Zp59+2px33nlm2bJlsT7Uvf2ef/5586Mf/cisXbvWSDLr1q1rtL01Nb7llltMv379TDAYNLt27TKTJk0yo0aNMqdPn45rLV0igIwfP97ccsstjdqGDh1qFixY4NKKurbKykojyWzdutUYY0x9fb3Jzc01S5cujfU5deqUyczMNL/85S/dWmaXUVNTY4YMGWKCwaC54oorYgGEuifGXXfdZS677LKzbqfuiXHttdea73//+43a/uEf/sHceOONxhjqnghnBpDW1Pjzzz83Pp/PrF69OtbnyJEjJikpyWzcuDGu/Xf6UzC1tbWqqKjQ1KlTG7VPnTpVr7/+ukur6tqqq6slSb169ZIkHTp0SMePH2/0O/D7/briiiv4HThgzpw5uvbaazVlypRG7dQ9MdavX6/CwkLNnDlT2dnZGjNmjFauXBnbTt0T47LLLtOLL76o9957T5L05ptv6tVXX9U111wjibrb0JoaV1RUKBqNNuqTl5enESNGxP17sP4sGKdVVVWprq5OOTk5jdpzcnJ0/Phxl1bVdRljVFpaqssuu0wjRoyQpFidm/sdfPTRR9bX2JWsXr1au3bt0s6dO5tso+6J8cEHH2j58uUqLS3V3XffrR07duj222+X3+/X7NmzqXuC3HXXXaqurtbQoUPl9XpVV1en+++/X9dff70k3u82tKbGx48fV0pKis4///wmfeL9zO30AaTBmU/WNcY4+rRdfOW2227TW2+9pVdffbXJNn4Hzjp8+LDmzZunTZs2KTU19az9qLuz6uvrVVhYqCVLlkiSxowZo3feeUfLly/X7NmzY/2ou7PWrFmj3/3ud1q1apWGDx+uPXv2aP78+crLy1NJSUmsH3VPvLbUuC2/h05/CiYrK0ter7dJ8qqsrGyS4tA+c+fO1fr16/Xyyy83eqJxbm6uJPE7cFhFRYUqKys1btw4JScnKzk5WVu3btVDDz2k5OTkWG2pu7P69u2rYcOGNWq78MIL9fHHH0vi/Z4oP/zhD7VgwQL9y7/8iy666CJ973vf0x133KHy8nJJ1N2G1tQ4NzdXtbW1OnHixFn7tFanDyApKSkaN26cgsFgo/ZgMKgJEya4tKquxRij2267Tc8884xeeuklDRo0qNH2QYMGKTc3t9HvoLa2Vlu3buV30A5XXnml9u7dqz179sT+FBYW6oYbbtCePXs0ePBg6p4AEydObPI18/fee08FBQWSeL8nSjgcVlJS448kr9cb+xoudU+81tR43Lhx8vl8jfocO3ZMb7/9dvy/hzZdOtvBNHwN91e/+pXZt2+fmT9/vklLSzMffvih20vrEm699VaTmZlptmzZYo4dOxb7Ew6HY32WLl1qMjMzzTPPPGP27t1rrr/+er4elwBf/xaMMdQ9EXbs2GGSk5PN/fffbw4cOGB+//vfm0AgYH73u9/F+lB355WUlJh+/frFvob7zDPPmKysLHPnnXfG+lD39qupqTG7d+82u3fvNpLMgw8+aHbv3h27bUVranzLLbeY/v37m82bN5tdu3aZyZMnd9+v4RpjzCOPPGIKCgpMSkqKGTt2bOwromg/Sc3+eeKJJ2J96uvrzb333mtyc3ON3+833/rWt8zevXvdW3QXdWYAoe6J8Ze//MWMGDHC+P1+M3ToULNixYpG26m780KhkJk3b54ZMGCASU1NNYMHDzY/+tGPTCQSifWh7u338ssvN/v/85KSEmNM62r85Zdfmttuu8306tXL9OjRw0ybNs18/PHHca+Fp+ECAADrOv01IAAAoPMhgAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAAADrCCAAEqK2ttbtJQDowJLdXgCArqGoqEgjRoxQSkqKfvOb32j48OHaunWr28sC0EERQAA45qmnntKtt96q1157TTznEkBLeBouAEcUFRWpurpau3fvdnspADoBrgEB4JjCwkK3lwCgkyCAAHBMWlqa20sA0EkQQAAAgHUEEAAAYB0BBAAAWMe3YAAAgHUcAQEAANYRQAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFj3/wHd4Y6IDTdYBgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Input parameters for grid and evolution here\n",
    "max_r = 96.0 # outer edge of the grid (including ghosts)\n",
    "num_points_r = 128 # total (including ghosts)\n",
    "log_factor = 1.04 # increase in dr interval at each grid step\n",
    "my_grid = Grid(max_r, num_points_r, log_factor)\n",
    "r = my_grid.r_vector\n",
    "\n",
    "initial_state = get_initial_state(my_grid)\n",
    "\n",
    "#unpackage the vector for readability\n",
    "(initial_u, initial_v , initial_phi, initial_hrr, initial_htt, initial_hpp, \n",
    " initial_K, initial_arr, initial_att, initial_app, \n",
    " initial_lambdar, initial_shiftr, initial_br, initial_lapse) = np.array_split(initial_state, NUM_VARS)\n",
    "\n",
    "#plot initial conditions\n",
    "plt.xlabel('r')\n",
    "plt.plot(r, initial_u, label='u')\n",
    "#plt.plot(r, initial_v, label='v')\n",
    "#plt.plot(r, initial_arr, label='arr')\n",
    "#plt.plot(r, initial_att, label='att')\n",
    "plt.plot(r, initial_K, label='K')\n",
    "#plt.plot(r, initial_phi, label='phi')\n",
    "plt.plot(r, initial_hrr, '-o', label='hrr') # zero, but plot as dots to see the grid\n",
    "#plt.plot(r, initial_lapse, label='lapse')\n",
    "#plt.plot(r, initial_lambdar, label='lambdar')\n",
    "plt.legend(loc='best')\n",
    "#plt.xlim(-1.2,4.2)\n",
    "#plt.ylim(-0.001,0.001)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5515ca2e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAAHACAYAAACLRZ3fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABMwElEQVR4nO3de1xVVeL///cBDoeLggJxM1Q0U1MqglK00qaksnuT2lTUTGZfx2nKmPk0MdWMOp90aj5jdp+xcbKy0sf8zKkmMrHULPF+S23UTEMNRBQBReDA2b8/8Bw9cvEAGw4HXs/Hg4fsfdZeZ+3F7e3aa69tMQzDEAAAAFrEz9sNAAAA6AgIVQAAACYgVAEAAJiAUAUAAGACQhUAAIAJCFUAAAAmIFQBAACYgFAFAABgAkIVAACACQhVAAAAJiBUmejLL7/ULbfcovj4eFksFv373/9u9fc8ePCg7rvvPkVGRiokJESXXnqpNmzY0OrvCwAA3BGqTHTixAldcskleuWVV9rk/YqLizV8+HBZrVZ9+umn2rFjh/7617+qW7dubfL+AADgNAsPVG4dFotFixYt0u233+7aV1VVpaefflrvvvuujh07psGDB+u5557TyJEjm/UeTz75pL7++mutXLnSnEYDAIBmY6SqDf3iF7/Q119/rfnz52vr1q0aM2aMbrjhBu3evbtZ9X300UdKTU3VmDFjFB0dreTkZL3xxhsmtxoAAHiCkapWcvZI1Z49e9SvXz8dOHBA8fHxrnLXXXedrrjiCk2fPr3J7xEUFCRJyszM1JgxY7R27VpNnjxZf//733X//febch4AAMAzAd5uQGexceNGGYahCy+80G1/ZWWlIiMjJUn79u1TYmJio/X86le/cs3ZcjgcSk1NdQWy5ORkbd++Xa+//jqhCgCANkaoaiMOh0P+/v7asGGD/P393V7r0qWLJKlHjx769ttvG62ne/furs/j4uJ00UUXub0+cOBALVy40KRWAwAATxGq2khycrJqampUWFioq666qt4yVqtVAwYM8LjO4cOHa+fOnW77du3apV69erWorQAAoOkIVSY6fvy4vvvuO9f23r17tXnzZkVEROjCCy/Uvffeq/vvv19//etflZycrKKiIn3xxRdKSkrS6NGjm/x+jz/+uIYNG6bp06dr7NixWrt2rWbPnq3Zs2ebeVoAAMADTFQ30fLly3XNNdfU2f/AAw9o7ty5stvt+t///V+9/fbbOnjwoCIjI5WWlqapU6cqKSmpWe/5n//8R1lZWdq9e7cSExOVmZmpCRMmtPRUAABAE/lUqPryyy/1l7/8RRs2bFB+fn6ddaDqs2LFCmVmZmr79u2Kj4/XE088oYkTJ7qVWbhwoZ555hnt2bNHffv21bPPPqs77rijFc8EAAB0ND61TlVTVyzfu3evRo8erauuukqbNm3S73//ez366KNuE7lzc3M1btw4ZWRkaMuWLcrIyNDYsWO1Zs2a1joNAADQAfnUSNWZ6lux/Gy/+93v9NFHH7ndUTdx4kRt2bJFubm5kqRx48aptLRUn376qavMDTfcoO7du+v9999vtfYDAICOpUNPVM/NzVV6errbvuuvv15z5syR3W6X1WpVbm6uHn/88TplZs2a1WC9lZWVqqysdG07HA4dPXpUkZGRslgspp4DAABoHYZhqKysTPHx8fLza/nFuw4dqgoKChQTE+O2LyYmRtXV1SoqKlJcXFyDZQoKChqsd8aMGZo6dWqrtBkAALSt/fv36/zzz29xPR06VEmqM3LkvNp55v76yjQ24pSVlaXMzEzXdklJiXr27Km9e/eqa9euZjRbkjTx3U1as7fYtb3p6Z+YVjcaZ7fbtWzZMl1zzTWyWq3ebk6nQb+3PfrcO+h37zi738vKypSYmGja3+4OHapiY2PrjDgVFhYqICDA9WiYhsqcPXp1JpvNJpvNVmd/RESEwsLCTGj5qfcJ6So/2+nLjM42o/XZ7XaFhIQoMjKSX3htiH5ve/S5d9Dv3nF2vzv73qypOz51919TpaWlKScnx23fkiVLlJqa6urIhsoMGzaszdoJAAB8n0+NVDW2YnnPnj2VlZWlgwcP6u2335ZUe6ffK6+84loQMzc3V3PmzHG7q++xxx7T1Vdfreeee0633XabPvzwQy1dulRfffVVm58fAADwXT41UrV+/XolJycrOTlZkpSZmank5GT94Q9/kCTl5+crLy/PVT4xMVHZ2dlavny5Lr30Uv3pT3/SSy+9pJ/+9KeuMsOGDdP8+fP15ptv6uKLL9bcuXO1YMECDRkypG1PDgAA+DSfGqkaOXKkGltWa+7cuXX2jRgxQhs3bmy03rvuukt33XVXS5sHAAA6MZ8aqQIAAGivCFUAAAAmIFQBAACYgFAFAABgAkIVAACACQhVAAAAJiBUAQAAmIBQBQAAYAJCFQAAgAkIVQAAACYgVAEAAJiAUAUAAGACQhUAAIAJCFUAAAAmIFQBAACYgFAFAABgAkIVAACACQhVAAAAJiBUAQAAmIBQBQAAYAJCFQAAgAkIVQAAACYgVAEAAJiAUAUAAGACQhUAAIAJCFUAAAAmIFQBAACYgFAFAABgAkIVAACACQhVAAAAJiBUAQAAmIBQBQAAYAJCFQAAgAkIVQAAACYgVAEAAJjA50LVa6+9psTERAUFBSklJUUrV65ssOzPf/5zWSyWOh+DBg1ylZk7d269ZSoqKtridAAAQAfhU6FqwYIFmjx5sp566ilt2rRJV111lW688Ubl5eXVW/7FF19Ufn6+62P//v2KiIjQmDFj3MqFhYW5lcvPz1dQUFBbnBIAAOggfCpUzZw5U+PHj9dDDz2kgQMHatasWUpISNDrr79eb/nw8HDFxsa6PtavX6/i4mL94he/cCtnsVjcysXGxrbF6QAAgA4kwNsN8FRVVZU2bNigJ5980m1/enq6Vq1a5VEdc+bM0XXXXadevXq57T9+/Lh69eqlmpoaXXrppfrTn/6k5OTkBuuprKxUZWWla7u0tFSSZLfbZbfbPT2lc3IYDrdtM+tG45x9TZ+3Lfq97dHn3kG/e8fZ/W52//tMqCoqKlJNTY1iYmLc9sfExKigoOCcx+fn5+vTTz/Ve++957Z/wIABmjt3rpKSklRaWqoXX3xRw4cP15YtW9SvX79665oxY4amTp1aZ/+SJUsUEhLShLNqXFGRn84cTMzOzjatbngmJyfH203olOj3tkefewf97h3Ofi8vLze1Xp8JVU4Wi8Vt2zCMOvvqM3fuXHXr1k2333672/6hQ4dq6NChru3hw4frsssu08svv6yXXnqp3rqysrKUmZnp2i4tLVVCQoLS09MVFhbWhLNp3ILC9dpVctS1PXr0aNPqRuPsdrtycnI0atQoWa1Wbzen06Df2x597h30u3ec3e/OK01m8ZlQFRUVJX9//zqjUoWFhXVGr85mGIb++c9/KiMjQ4GBgY2W9fPz0+WXX67du3c3WMZms8lms9XZb7VaTf3h8LO4T3njB6/tmf01hWfo97ZHn3sH/e4dzn43u+99ZqJ6YGCgUlJS6gyV5uTkaNiwYY0eu2LFCn333XcaP378Od/HMAxt3rxZcXFxLWovAADoXHxmpEqSMjMzlZGRodTUVKWlpWn27NnKy8vTxIkTJdVeljt48KDefvttt+PmzJmjIUOGaPDgwXXqnDp1qoYOHap+/fqptLRUL730kjZv3qxXX321Tc4JAAB0DD4VqsaNG6cjR45o2rRpys/P1+DBg5Wdne26my8/P7/OmlUlJSVauHChXnzxxXrrPHbsmB5++GEVFBQoPDxcycnJ+vLLL3XFFVe0+vkAAICOw6dClSRNmjRJkyZNqve1uXPn1tkXHh7e6Oz+F154QS+88IJZzQMAAJ2Uz8ypAgAAaM8IVQAAACYgVAEAAJiAUAUAAGACQhUAAIAJCFUAAAAmIFQBAACYgFAFAABgAkIVAACACQhVAAAAJiBUAQAAmIBQBQAAYAJCFQAAgAkIVQAAACYgVAEAAJiAUAUAAGACQhUAAIAJCFUAAAAmIFQBAACYgFAFAABgAkIVAACACQhVAAAAJiBUAQAAmIBQBQAAYAJCFQAAgAkIVQAAACYgVAEAAJiAUAUAAGACQhUAAIAJCFUAAKBVfP1dkXb8WOrtZrQZQhUAADDd/qPlypizRhlz1qi6xuHt5rQJQhUAADDd9h9L5TCkIyeqtL2TjFYRqgAAgOl2Hypzfb76+yNebEnbIVQBAADT7S487vo8l1AFAADQPLvOGKlat/eo7J1gXpXPharXXntNiYmJCgoKUkpKilauXNlg2eXLl8tisdT5+O9//+tWbuHChbroootks9l00UUXadGiRa19GgAAdFjVNQ59X3RCkmT1t+hEVY2+OVji5Va1Pp8KVQsWLNDkyZP11FNPadOmTbrqqqt04403Ki8vr9Hjdu7cqfz8fNdHv379XK/l5uZq3LhxysjI0JYtW5SRkaGxY8dqzZo1rX06AAB0SPuLT6qq2qEgq59GXBgtqXPMq/KpUDVz5kyNHz9eDz30kAYOHKhZs2YpISFBr7/+eqPHRUdHKzY21vXh7+/vem3WrFkaNWqUsrKyNGDAAGVlZenaa6/VrFmzWvlsAADomJyX/vqe10XDL4iUJOXuIVS1G1VVVdqwYYPS09Pd9qenp2vVqlWNHpucnKy4uDhde+21WrZsmdtrubm5deq8/vrrz1knAACo33enJqlfGNNVaX1rQ9X6fcWqqu7Y86oCvN0ATxUVFammpkYxMTFu+2NiYlRQUFDvMXFxcZo9e7ZSUlJUWVmpd955R9dee62WL1+uq6++WpJUUFDQpDolqbKyUpWVla7t0tLa9Tfsdrvsdnuzzq8+DsP9m8/MutE4Z1/T522Lfm979Ll3dPR+/29+7fypPpHBSuwepO4hVhWX27VxX5FSenX3WrvO7nez+99nQpWTxWJx2zYMo84+p/79+6t///6u7bS0NO3fv1//93//5wpVTa1TkmbMmKGpU6fW2b9kyRKFhIR4dB6eKCry05mDidnZ2abVDc/k5OR4uwmdEv3e9uhz7+io/b7xO39JFpXs36nFi/+rnkF+Ki730zufrdah8w1vN8/V7+Xl5abW6zOhKioqSv7+/nVGkAoLC+uMNDVm6NChmjdvnms7Nja2yXVmZWUpMzPTtV1aWqqEhASlp6crLCzM47acy4LC9dpVctS1PXr0aNPqRuPsdrtycnI0atQoWa1Wbzen06Df2x597h0dud9rHIaeWPe5JIfuvnGEekWG6Ghknrb8578qtp6n0aNTvda2s/vdeaXJLD4TqgIDA5WSkqKcnBzdcccdrv05OTm67bbbPK5n06ZNiouLc22npaUpJydHjz/+uGvfkiVLNGzYsAbrsNlsstlsdfZbrVZTfzj8LO5T3jraD54vMPtrCs/Q722PPveOjtjvB4tOqLLaIVuAnxKjw+TvZ9GV/aIl/Vcb8o7JYfGTLcD/nPW0Jme/m933PhOqJCkzM1MZGRlKTU1VWlqaZs+erby8PE2cOFFS7QjSwYMH9fbbb0uqvbOvd+/eGjRokKqqqjRv3jwtXLhQCxcudNX52GOP6eqrr9Zzzz2n2267TR9++KGWLl2qr776yivnCACAL3OupN73vC7y96udSnNBdBdFdQlU0fEqbdlfoisSI7zZxFbjU6Fq3LhxOnLkiKZNm6b8/HwNHjxY2dnZ6tWrlyQpPz/fbc2qqqoq/fa3v9XBgwcVHBysQYMG6ZNPPnG7jDZs2DDNnz9fTz/9tJ555hn17dtXCxYs0JAhQ9r8/AAA8HXO5RT6xXRx7bNYLBrSJ1KfbM1X7p4jhKr2YtKkSZo0aVK9r82dO9dt+4knntATTzxxzjrvuusu3XXXXWY0DwCATs25nEK/6C5u+9Ocoer7Ij2mfvUd6vN8Zp0qAADQ/u0udI5UdXXb71yvamPeMVXYa9q8XW2BUAUAAEzhcBgNjlT1iQpVdFebqqod2pR3zAuta32EKgAAYIoDxSdVYXcoMMBPPSPc1220WCwa2ufUI2s66HMACVUAAMAUzknqfaJCFeBfN2I4LwGu7qDPASRUAQAAUziXUzh7PpVT2qmRqk37i3WyquPNqyJUAQAAUzgnqV941nwqp16RIYoNC5K9xtDGvOK2bFqbIFQBAABT7D7kHKmqP1RZLBbXJcDcDngJkFAFAABazO3OvwYu/0mnLwF2xMnqhCoAANBiB4+d1El7jaz+FvU6686/MzlHqrbsP6byquq2al6bIFQBAIAWc86n6hPVpd47/5zO7x6sHt2CVe0wtH5fx5pXRagCAAAtdq75VE4deb0qQhUAAGixXc5QFd3wfCqnjjpZnVAFAABa7DvncgrnGKmSpKF9IiRJ3xws0fHKjjOvilAFAABaxOEwzlj489yh6vzuIUqICFaNw9C6fUdbu3lthlAFAABa5MeSkyqvOnXnX2SoR8c4l1boSI+sIVQBAIAWcY5SJUaFytrInX9ncs2r6kCT1QlVAACgRXafepCyJ5PUnZx3AG47WKLSCnurtKutEaoAAECLeLqcwpniwoPVOzJEDkNat7djzKsiVAEAgBbZVej5cgpn6mhLKxCqAABAsxmGoe+cl/+aMFIlnb4EuHovoQoAAHRy+SUVOlFVowA/i3p7eOefk/MOwO0/lqqk3PfnVRGqAABAs+06NUrVOypUgQFNixXRYUHqc16oDENa0wFGqwhVAACg2b47NZ/Kk5XU65PWgZ4DSKgCAADN5hypuqCJk9SdnJPVV3/v+3cAEqoAAECzuR5PE928kaohibWh6tv8UhWfqDKtXd5AqAIAAM1Se+ef8/Jf80aqzutqcwUyX59XRagCAADNUlBaobLKavn7WdQ7KqTZ9XSU9aoIVQAAoFmcK6n3jgyRLcC/2fW4Hq7s4/OqCFUAAKBZdjXjmX/1GXIqVO08VKYjxytb3C5vIVQBAIBmcS6n0NSV1M8WERqoAbG1wcyXR6sIVQAAoFlcd/41c5L6mVyPrPHh9aoIVQAAoMkMwzjj8l/LRqqkMyarE6oAAEBnUlhWqbKKavlZpD7nNe2Zf/UZkhghi6X2kmJhWYUJLWx7hCoAANBkrmf+RYa26M4/p24hgRoYGybJd+dVEaoAAECTOZdTuMCES39Opx9Z45uXAH0uVL322mtKTExUUFCQUlJStHLlygbLfvDBBxo1apTOO+88hYWFKS0tTZ999plbmblz58pisdT5qKjwzaFHAADawu7Clq2kXh/XZHUfXQTUp0LVggULNHnyZD311FPatGmTrrrqKt14443Ky8urt/yXX36pUaNGKTs7Wxs2bNA111yjW265RZs2bXIrFxYWpvz8fLePoKCgtjglAAB80m7nJPUWLqdwpisSI+Rnkb4vOqFDpb43uOFToWrmzJkaP368HnroIQ0cOFCzZs1SQkKCXn/99XrLz5o1S0888YQuv/xy9evXT9OnT1e/fv308ccfu5WzWCyKjY11+wAAAPUzDMM1UmXm5b/wYKsGxYdL8s1H1vhMqKqqqtKGDRuUnp7utj89PV2rVq3yqA6Hw6GysjJFRES47T9+/Lh69eql888/XzfffHOdkSwAAHDa4bJKlZy0y88i9T3PvFAl+fa8qgBvN8BTRUVFqqmpUUxMjNv+mJgYFRQUeFTHX//6V504cUJjx4517RswYIDmzp2rpKQklZaW6sUXX9Tw4cO1ZcsW9evXr956KisrVVl5ehn90tJSSZLdbpfdbm/qqTXIYTjcts2sG41z9jV93rbo97ZHn3uHr/f7tz8ekyQldA+Rvxyy2x2NH9AEqT3DNVvSqj1FpvfP2f1udv0+E6qcLBaL27ZhGHX21ef999/XlClT9OGHHyo6Otq1f+jQoRo6dKhre/jw4brsssv08ssv66WXXqq3rhkzZmjq1Kl19i9ZskQhIc1/SvfZior8dOZgYnZ2tml1wzM5OTnebkKnRL+3PfrcO3y137/Mt0jyV5hx3PS/TRXVkp/8lXf0pN5dlK3uNlOrl3S638vLy02t12dCVVRUlPz9/euMShUWFtYZvTrbggULNH78eP3rX//Sdddd12hZPz8/XX755dq9e3eDZbKyspSZmenaLi0tVUJCgtLT0xUWFubB2XhmQeF67So5vVbH6NGjTasbjbPb7crJydGoUaNktVq93ZxOg35ve/S5d/h6v+d+tEPad0DDk/pq9Kj6r+q0xLv5q7X1QKlCEy/V6EvjTav37H53Xmkyi8+EqsDAQKWkpCgnJ0d33HGHa39OTo5uu+22Bo97//339eCDD+r999/XTTfddM73MQxDmzdvVlJSUoNlbDabbLa60dlqtZr6w+FncZ/y5os/eL7O7K8pPEO/tz363Dt8td+/P1w7wtM/LqxV2j+s73naeqBUa/cd05jLe5lev7PfzW67z4QqScrMzFRGRoZSU1OVlpam2bNnKy8vTxMnTpRUO4J08OBBvf3225JqA9X999+vF198UUOHDnWNcgUHBys8vPbugqlTp2ro0KHq16+fSktL9dJLL2nz5s169dVXvXOSAAC0Y4ZhaFeh85l/5q1RdaahfSL0txV7fO45gD4VqsaNG6cjR45o2rRpys/P1+DBg5Wdna1evWpTbH5+vtuaVX//+99VXV2tX/3qV/rVr37l2v/AAw9o7ty5kqRjx47p4YcfVkFBgcLDw5WcnKwvv/xSV1xxRZueGwAAvqDoeJWOldtlaYU7/5wu7x2hAD+LDhSf1P6j5UqIMG++cmvyqVAlSZMmTdKkSZPqfc0ZlJyWL19+zvpeeOEFvfDCCya0DACAjm/3qVGqnhEhCg5s+TP/6hNqC9DF54drY94xrf7+iM+EKp9ZpwoAAHif85l//Uxc9LM+zkfW+NIlQEIVAADwmHOk6oJWmk/l5FoEdM8RGYbRqu9lFkIVAADwmHOk6kITn/lXn9ReEbL6W/RjSYXyjpq7nlRraXaoqqqq0s6dO1VdXW1mewAAQDvmfOZfa9355xQc6K9LE7pJ8p1H1jQ5VJWXl2v8+PEKCQnRoEGDXHfbPfroo/rzn/9segMBAED7cOR4pY6eqJLFYu6DlBvimlflIw9XbnKoysrK0pYtW7R8+XIFBQW59l933XVasGCBqY0DAADtx65Tl/7O7x7canf+nSntjMnqvjCvqslLKvz73//WggULNHToULdn7l100UXas2ePqY0DAADtx3etvOjn2S7r1V2B/n46VFqpfUfKlRgV2ibv21xNHqk6fPiw2wOJnU6cOOHRg40BAIBvcs2nauVJ6k5BVn8l9+wmyTcuATY5VF1++eX65JNPXNvOIPXGG28oLS3NvJYBAIB2Zdehth2pkqTrB8Xq5ovj1NMHFgBt8uW/GTNm6IYbbtCOHTtUXV2tF198Udu3b1dubq5WrFjRGm0EAADtwHeFbbOcwpkevDJRDyqxzd6vJZo8UjVs2DB9/fXXKi8vV9++fbVkyRLFxMQoNzdXKSkprdFGAADgZUdPVKnoeJWk1nvmn69r1rP/kpKS9NZbb5ndFgAA0E7tPnXpr0e3YIXafO7RwW2iyb3iXJeqIT179mx2YwAAQPu02wuX/nxNk0NV7969G73Lr6ampkUNAgAA7Y9zpKpfTNtNUvc1TQ5VmzZtctu22+3atGmTZs6cqWeffda0hgEAgPbj9ONpGKlqSJND1SWXXFJnX2pqquLj4/WXv/xFd955pykNAwAA7YdzNXVGqhrW7Acqn+3CCy/UunXrzKoOAAC0E8UnqlR0vFJS2zzzz1c1eaSqtLTUbdswDOXn52vKlCnq16+faQ0DAADtw3eHa0epenQLVhfu/GtQk3umW7dudSaqG4ahhIQEzZ8/37SGAQCA9sG5kjqjVI1rcqhatmyZ27afn5/OO+88XXDBBQoIIL0CANDR7D7EcgqeaHIKGjFiRGu0AwAAtFO7C9v+mX++yKNQ9dFHH3lc4a233trsxgAAgPbHOVJ1ASNVjfIoVN1+++0eVWaxWFj8EwCADqSk3K7Csto7/1ijqnEehSqHw9Ha7QAAAO2Q89JfXHiQugZZvdya9s20daoAAEDH41pJnUU/z6lZt+udOHFCK1asUF5enqqqqtxee/TRR01pGAAA8D7ncgpc+ju3Zj37b/To0SovL9eJEycUERGhoqIihYSEKDo6mlAFAEAH8h3P/PNYky//Pf7447rlllt09OhRBQcHa/Xq1frhhx+UkpKi//u//2uNNgIAAC/ZzTP/PNbkULV582b95je/kb+/v/z9/VVZWamEhAQ9//zz+v3vf98abQQAAF5QctKugtIKSaym7okmhyqr1ep6TE1MTIzy8vIkSeHh4a7PAQCA73Ne+osNC1J4MHf+nUuT51QlJydr/fr1uvDCC3XNNdfoD3/4g4qKivTOO+8oKSmpNdoIAAC8YLdzkjqLfnqkySNV06dPV1xcnCTpT3/6kyIjI/XLX/5ShYWFmj17tukNBAAA3uFcToFLf55p8khVamqq6/PzzjtP2dnZpjYIAAC0D85QdSGT1D3S5JGqqVOnas+ePa3RFgAA0I7sZo2qJmlyqFq4cKEuvPBCDR06VK+88ooOHz7cGu0CAABeVFZhV35J7Z1//aIZqfJEk0PV1q1btXXrVv3kJz/RzJkz1aNHD40ePVrvvfeeysvLW6ONbl577TUlJiYqKChIKSkpWrlyZaPlV6xYoZSUFAUFBalPnz7629/+VqfMwoULddFFF8lms+miiy7SokWLWqv5AAD4BOelv+iuNoWHcOefJ5r17L9BgwZp+vTp+v7777Vs2TIlJiZq8uTJio2NNbt9bhYsWKDJkyfrqaee0qZNm3TVVVfpxhtvbHAph71792r06NG66qqrtGnTJv3+97/Xo48+qoULF7rK5Obmaty4ccrIyNCWLVuUkZGhsWPHas2aNa16LgAAtGffuRb95NKfp1r8QOXQ0FAFBwcrMDBQdrvdjDY1aObMmRo/frweeughDRw4ULNmzVJCQoJef/31esv/7W9/U8+ePTVr1iwNHDhQDz30kB588EG3ld9nzZqlUaNGKSsrSwMGDFBWVpauvfZazZo1q1XPBQCA9mx3oXM+FZf+PNWsByrv3btX7733nt59913t2rVLV199taZMmaIxY8aY3T6XqqoqbdiwQU8++aTb/vT0dK1atareY3Jzc5Wenu627/rrr9ecOXNkt9tltVqVm5urxx9/vE4ZQhUANJ1hGPXsq6dcU46vt1x9ddZfa31lPS3XlDrrb2f951Ntt+tkde28pYCaRtrpYZua0see+jafNaqaqsmhKi0tTWvXrlVSUpJ+8Ytf6J577lGPHj1ao21uioqKVFNTo5iYGLf9MTExKigoqPeYgoKCestXV1erqKhIcXFxDZZpqE5JqqysVGVlpWu7tLRUkmS3200drXMYDrft1h4JxGnOvm7vfW4Yhuw1huw1DlXVOGSvMVRV7ZC9xqHqGkPVDkMOo/bfGoehaodDDodObTtqX3dI1Q6Halxlzj7GkMNRd7vGYchhSA6jtrzh+rz2l7lxar/DcH4uV7kzt898vabGoUOFfvrw6EZJlrr1GrV/Tuqr12HUvrHj1F8SQ873cm7LtWGc+tT5x8kwTu9zHW2c/qNUX9m6dZ3+mhiuNzzdFw29r9u/p963vrIyGqrr9Pue/r6o53ulge+f2n/99fjqJY2Uq2cnTBCgJ9ct83YjzikxIrjd/y701Nm/280+ryaHqmuuuUb/+Mc/NGjQIFMb4innI3KcDMOos+9c5c/e39Q6Z8yYoalTp9bZv2TJEoWEhDTc+CYqKvLTmVdoWROs7eXk5DT5GMOQqhzSyWqpokaqdEj2GqnSYZHdIVXV1L5e5frccsbnp/dXGxZVO6RqQ6o5699qQ6p2SDVGw9+nvstPOlrk7UZ0MhaCE+oVF2zo4De5OrzD2y0xl/N3u9k32DU5VE2fPt3UBngqKipK/v7+dUaQCgsL64w0OcXGxtZbPiAgQJGRkY2WaahOScrKylJmZqZru7S0VAkJCUpPT1dYWFiTzqsxCwrXa1fJUdf26NGjTasbjbPb7crJydGVI3+isipDR09U6ciJKh09UaWjJ+w6cqJKxeVVKj1ZrdIKu8oqqlVaUa3jFdUqq6xWjcM7f6H8LFJggJ+s/n4K8LMowM8if9e/fvL3k/xPfR7gZ5GfnxTg5+cq42c5fcyZH2fW43dGWX8/i/wstf8x8bNIFrlv+1ksstSz7Wc5Xe7MbcPh0H//+60GXXSRAvz93Y6rr7xbvZJ0xnuc2nT9B8n5+qmSte069W/tS6ePO/362WXPrstZ1tKMuk6Xra8uV9kzjlUDbT7z2DPf82z1RfCa6mqt+PJLjbj6agVYrQ2W87S++gp6Xl89x3r4vvX/P9jc+lrSvrMLVtvt+vzzz3XtddfKGmBtqFgj7Wt+PzdVY4MMvsb5u33UqFGyWq2uK01madacKm8IDAxUSkqKcnJydMcdd7j25+Tk6Lbbbqv3mLS0NH388cdu+5YsWaLU1FRZT/3ySEtLU05Ojtu8qiVLlmjYsGENtsVms8lms9XZb7VaXfWawc/ifh+BmXVDqqp26EBxufJLKnTw2EnlH6tQfslJHTx2Uj8eO6kDR/xVmftls+sP8LOoS1CAQqz+Cg6s/QixBigo0F8hVn+FBPq7Pne+Huzcb/WXLcBfgQEWBfr7y+pvcYUl26l/A8/4N/DUv/5+vv3Lz263K7t4h0YP6cX3exux2+0KD5TiI7rQ523IHuCnAD8pNMhGv3uB8++12X3vM6FKkjIzM5WRkaHU1FSlpaVp9uzZysvL08SJEyXVjiAdPHhQb7/9tiRp4sSJeuWVV5SZmakJEyYoNzdXc+bM0fvvv++q87HHHtPVV1+t5557Trfddps+/PBDLV26VF999ZVXzhHmKzlp157Dx7Wn8Lj2HD6h7wqP6/vDx/XD0fJzjCjVBpRAfz9FhAYqsktg7b+hgYrsYlNEaKDCg60KC7aqa1CAwoKsCgsKcG0HW/071P/wAACN86lQNW7cOB05ckTTpk1Tfn6+Bg8erOzsbPXq1UuSlJ+f77ZmVWJiorKzs/X444/r1VdfVXx8vF566SX99Kc/dZUZNmyY5s+fr6efflrPPPOM+vbtqwULFmjIkCFtfn5oucrqGu34sVSb9x/T5v3HtCnvmPKONnzNPNjqr/huQYrvFqz48GDFdwtWXLcgRXexatfmNbrrplHq3iWYcAQAOCefClWSNGnSJE2aNKne1+bOnVtn34gRI7Rx48ZG67zrrrt01113mdE8tLEKe43W7TuqlbuLtHbvUe34sVRVNY465WLDgtQ3OlR9z+uiC6K7qO95tR8xYbZ6A5PdblfJTqlrkJVABQDwSLNCVUVFhbZu3arCwkI5HO5/wG699VZTGgY05MjxSi399pAWbyvQqj1HVFnt/j0YERqoSxO66dKEbkru2U0Xn99N4cHMWQAAtK4mh6rFixfr/vvvV1FR3VueLRaLampqTGkYcKYKe40+216g/2/DAX39XZHOnAoVGxakq/pFafgFUUru2U09I0IYXQIAtLkmh6pHHnlEY8aM0R/+8IdGlx0AzLD/aLne/Hqf/rVhv8oqql37B/cI0w2DYjXqolhdGNOFEAUA8Lomh6rCwkJlZmYSqNCq/ltQqpc+363F2wpco1I9ugXrpynn66eX9VCvyFDvNhAAgLM0OVTdddddWr58ufr27dsa7UEn98ORE/rrkl36eOuPrhWer+oXpfFXJurqfufJz8fXYQIAdFxNDlWvvPKKxowZo5UrVyopKanOwlmPPvqoaY1D51Fhr9HfVuzRa8v3qOrUxPObkuL062sv0IBY81apBwCgtTQ5VL333nv67LPPFBwcrOXLl9d5hh6hCk21Ma9YmQs2a9+R2vWkrrwgSk/eOECDe4R7uWUAAHiuyaHq6aef1rRp0/Tkk0/Kz8/v3AcADaiucejVZXv00he7VeMwFBNm0x9uHqTRSbFMPAcA+Jwmh6qqqiqNGzeOQIUWKTlp1yPvbdTK3bVLc9x2abz+dPtghQWxnhQAwDc1ORk98MADWrBgQWu0BZ3EvqITuvO1r7Vyd5GCrf56YdwlevHuZAIVAMCnNXmkqqamRs8//7w+++wzXXzxxXUmqs+cOdO0xqHj2VlQpnveWK0jJ6oUFx6kfzyQqkHxzJ0CAPi+Joeqb775RsnJyZKkbdu2ub3GPBg0ZsePpbpvzhodPVGli+LCNPcXlys6LMjbzQIAwBRNDlXLli1rjXagg9tbdEL3/mO1isvtuvj8cL3z4BCFh3C5DwDQcTTrgcpAUxw9UaVfvLlWxeV2JfUI1zvjh/CAYwBAh9OsULVu3Tr961//Ul5enqqqqtxe++CDD0xpGDqGyuoaPfz2eu07Uq7zuwfrnz+/nEAFAOiQmnz33/z58zV8+HDt2LFDixYtkt1u144dO/TFF18oPJwJx3D33Kc7tf6HYoUFBWjuLy7XeV1t3m4SAACtosmhavr06XrhhRf0n//8R4GBgXrxxRf17bffauzYserZs2drtBE+aumOQ/rn13slSTPHXqoLort6uUUAALSeJoeqPXv26KabbpIk2Ww2nThxQhaLRY8//rhmz55tegPhmwpLK/Q//98WSdKDwxN13UUxXm4RAACtq8mhKiIiQmVlZZKkHj16uJZVOHbsmMrLy81tHXzW1P/sUHG5XYPiw/S7G/t7uzkAALS6Jk9Uv+qqq5STk6OkpCSNHTtWjz32mL744gvl5OTo2muvbY02wscs+2+hPtmaL38/i56/62LZAvy93SQAAFpdk0PVK6+8ooqKCklSVlaWrFarvvrqK91555165plnTG8gfMvJqho9/e/a0csHh/dmtXQAQKfR5FAVERHh+tzPz09PPPGEnnjiCVMbBd81d9U+HTx2UvHhQXp81IXebg4AAG3G41BVWlrqUbmwsLBmNwa+raTcrteXfydJ+k16f4UEsrYsAKDz8PivXrdu3Rp9tp9hGLJYLKqpqTGlYfA9r6/Yo9KKavWP6arbk3t4uzkAALQpj0PVmc/8MwxDo0eP1j/+8Q/16MEfT0hHjldq7qraNameuKG//P14uDYAoHPxOFSNGDHCbdvf319Dhw5Vnz59TG8UfM/buT+owu7QJeeH6ycDor3dHAAA2lyT16kCznayqkbvrP5BkvTw1X0bvUwMAEBHRahCiy3ceEBHT1QpISJY1w9i5XQAQOfUolDFiAQMw9Cbp57vN354ogL8yekAgM7J4zlVd955p9t2RUWFJk6cqNDQULf9H3zwgTktg09Y/0Ox9hw+oZBAf92VmuDt5gAA4DUeh6rwcPeVse+77z7TGwPfs2DdfknSTUlx6mJjXSoAQOfl8V/BN998szXbAR9UVmHXJ1vzJUnjLmeUCgDQuTEBBs32n635OmmvUZ/zQpXSq7u3mwMAgFcRqtBsH2/5UZI0JiWBmxYAAJ0eoQrNcuR4pVZ/f0RS7XwqAAA6O0IVmuWz7YfkMKTBPcLUMzLE280BAMDrfCZUFRcXKyMjQ+Hh4QoPD1dGRoaOHTvWYHm73a7f/e53SkpKUmhoqOLj43X//ffrxx9/dCs3cuRIWSwWt4+77767lc/G9326rXaC+o2DGaUCAEDyoVB1zz33aPPmzVq8eLEWL16szZs3KyMjo8Hy5eXl2rhxo5555hlt3LhRH3zwgXbt2qVbb721TtkJEyYoPz/f9fH3v/+9NU/F5xWfqNKqPbWX/kZz6Q8AAElNWFLBm7799lstXrxYq1ev1pAhQyRJb7zxhtLS0rRz507179+/zjHh4eHKyclx2/fyyy/riiuuUF5ennr27OnaHxISotjY2NY9iQ7ky92HVeMwNCC2qxKjQs99AAAAnYBPhKrc3FyFh4e7ApUkDR06VOHh4Vq1alW9oao+JSUlslgs6tatm9v+d999V/PmzVNMTIxuvPFG/fGPf1TXrl0brKeyslKVlZWu7dLSUkm1lxztdnsTzqxxDsPhtm1m3S2x7NtDkqSrLohsN20ym/O8Our5tVf0e9ujz72DfveOs/vd7P73iVBVUFCg6OjoOvujo6NVUFDgUR0VFRV68skndc899ygsLMy1/95771ViYqJiY2O1bds2ZWVlacuWLXVGuc40Y8YMTZ06tc7+JUuWKCTEvEnbRUV+OvMKbXZ2tml1N5fDkJbu8Jdkke3oHmVnf+ftJrWqxr4P0Hro97ZHn3sH/e4dzn4vLy83tV6vhqopU6bUG07OtG7dOkn1P7zZMAyP1key2+26++675XA49Nprr7m9NmHCBNfngwcPVr9+/ZSamqqNGzfqsssuq7e+rKwsZWZmurZLS0uVkJCg9PR0t8DWUgsK12tXyVHX9ujRo02ru7m2/1iq46tXKyTQX78cc50CA3xmWl6T2O125eTkaNSoUbJard5uTqdBv7c9+tw76HfvOLvfnVeazOLVUPXII4+c80673r17a+vWrTp06FCd1w4fPqyYmJhGj7fb7Ro7dqz27t2rL7744pyh57LLLpPVatXu3bsbDFU2m002m63OfqvVauoPh5/FPbC0hx+8r78vliQN6xul0OC6fdDRmP01hWfo97ZHn3sH/e4dzn43u++9GqqioqIUFRV1znJpaWkqKSnR2rVrdcUVV0iS1qxZo5KSEg0bNqzB45yBavfu3Vq2bJkiIyPP+V7bt2+X3W5XXBx3tdVn5e7DkqQR/c/zcksAAGhffOLazcCBA3XDDTdowoQJWr16tVavXq0JEybo5ptvdpukPmDAAC1atEiSVF1drbvuukvr16/Xu+++q5qaGhUUFKigoEBVVVWSpD179mjatGlav3699u3bp+zsbI0ZM0bJyckaPny4V861PausrtHGvGOSpOF9zx1QAQDoTHwiVEm1d+glJSUpPT1d6enpuvjii/XOO++4ldm5c6dKSkokSQcOHNBHH32kAwcO6NJLL1VcXJzrY9WqVZKkwMBAff7557r++uvVv39/Pfroo0pPT9fSpUvl7+/f5ufY3m09UKKqaoeiugSylAIAAGfxibv/JCkiIkLz5s1rtIxhGK7Pe/fu7bZdn4SEBK1YscKU9nUGa/fWTpq/vHcED1AGAOAsPjNSBe9bt+90qAIAAO4IVfBIjcPQhn21d/4RqgAAqItQBY/8t6BUZZXVCg3018C4hlebBwCgsyJUwSObTt31d1mv7grw59sGAICz8dcRHvnmQO1dlZec3827DQEAoJ0iVMEjWw/Whqqk88O93BIAANonQhXOqcJeo12HyiRJFxOqAACoF6EK57Qjv1Q1DkNRXWyKDQvydnMAAGiXCFU4J+d8qovPD2fRTwAAGkCowjltPRWqknpw6Q8AgIYQqnBO2w4SqgAAOBdCFRpVVe3QnsPHJUkD48O83BoAANovQhUate/ICVU7DHWxBSg+nEnqAAA0hFCFRu0sqF1K4cKYLkxSBwCgEYQqNMq5PlX/WC79AQDQGEIVGuUcqeof08XLLQEAoH0jVKFRO0+NVF0Y29XLLQEAoH0jVKFB5VXVyjtaLknqH0OoAgCgMYQqNOi7wuMyDCkyNFCRXWzebg4AAO0aoQoN2lt0QpLUN5r5VAAAnAuhCg36/nBtqOoTFerllgAA0P4RqtAg50hVb0IVAADnRKhCg/YdqQ1ViYQqAADOiVCFehmGob2HCVUAAHiKUIV6HTlRpbLKalksUs+IEG83BwCAdo9QhXo551PFhwcryOrv5dYAAND+EapQL+elvz7ncekPAABPEKpQr72nJqn3jiRUAQDgCUIV6pV3pPbxNL0imU8FAIAnCFWo14FjJyVJ53cnVAEA4AlCFep1sLh2pOr87sFebgkAAL6BUIU6TlbVqOh4lSQpgZEqAAA8QqhCHQeP1Y5SdbUFKCw4wMutAQDANxCqUMf+4tr5VD26B8tisXi5NQAA+AZCFeo4UMwkdQAAmspnQlVxcbEyMjIUHh6u8PBwZWRk6NixY40e8/Of/1wWi8XtY+jQoW5lKisr9etf/1pRUVEKDQ3VrbfeqgMHDrTimbR/B5ikDgBAk/lMqLrnnnu0efNmLV68WIsXL9bmzZuVkZFxzuNuuOEG5efnuz6ys7PdXp88ebIWLVqk+fPn66uvvtLx48d18803q6amprVOpd076BqpIlQBAOApn5iF/O2332rx4sVavXq1hgwZIkl64403lJaWpp07d6p///4NHmuz2RQbG1vvayUlJZozZ47eeecdXXfddZKkefPmKSEhQUuXLtX1119v/sn4AC7/AQDQdD4RqnJzcxUeHu4KVJI0dOhQhYeHa9WqVY2GquXLlys6OlrdunXTiBEj9Oyzzyo6OlqStGHDBtntdqWnp7vKx8fHa/DgwVq1alWDoaqyslKVlZWu7dLSUkmS3W6X3W5v0bmeyWE43LbNrLsxzst/sV2tbfae7Y3zvDvr+XsL/d726HPvoN+94+x+N7v/fSJUFRQUuILQmaKjo1VQUNDgcTfeeKPGjBmjXr16ae/evXrmmWf0k5/8RBs2bJDNZlNBQYECAwPVvXt3t+NiYmIarXfGjBmaOnVqnf1LlixRSIh5oztFRX468wrt2ZcuW0O1Qyo6XvttsW3dV/rB2upv2a7l5OR4uwmdEv3e9uhz76DfvcPZ7+Xl5abW69VQNWXKlHrDyZnWrVsnSfXe2m8YRqO3/I8bN871+eDBg5WamqpevXrpk08+0Z133tngceeqNysrS5mZma7t0tJSJSQkKD09XWFhYY2eT1MsKFyvXSVHXdujR482re6GHDx2UlqzUlZ/i8bccqP8/Drnkgp2u105OTkaNWqUrNZOnizbEP3e9uhz76DfvePsfndeaTKLV0PVI488orvvvrvRMr1799bWrVt16NChOq8dPnxYMTExHr9fXFycevXqpd27d0uSYmNjVVVVpeLiYrfRqsLCQg0bNqzBemw2m2w2W539VqvV1B8OP4v7fQRt8YN3pLxMkhQTFiSbLbDV36+9M/trCs/Q722PPvcO+t07nP1udt97NVRFRUUpKirqnOXS0tJUUlKitWvX6oorrpAkrVmzRiUlJY2Gn7MdOXJE+/fvV1xcnCQpJSVFVqtVOTk5Gjt2rCQpPz9f27Zt0/PPP9+MM/J9BSW1c8Viw4K83BIAAHyLTyypMHDgQN1www2aMGGCVq9erdWrV2vChAm6+eab3SapDxgwQIsWLZIkHT9+XL/97W+Vm5urffv2afny5brlllsUFRWlO+64Q5IUHh6u8ePH6ze/+Y0+//xzbdq0Sffdd5+SkpJcdwN2NgWlFZKkmHBCFQAATeETE9Ul6d1339Wjjz7qulPv1ltv1SuvvOJWZufOnSopKZEk+fv765tvvtHbb7+tY8eOKS4uTtdcc40WLFigrl27uo554YUXFBAQoLFjx+rkyZO69tprNXfuXPn7+7fdybUjh06FKkaqAABoGp8JVREREZo3b16jZQzDcH0eHByszz777Jz1BgUF6eWXX9bLL7/c4jZ2BAUlhCoAAJrDJy7/oe04L//FcvkPAIAmIVTBzSFCFQAAzUKogothGFz+AwCgmQhVcCk5aVdlde2jcaLD6q7DBQAAGkaogotzPlVEaKBsAZ3z7kcAAJqLUAWXorIqSdJ5XRilAgCgqQhVcCk6XruaelRXHk8DAEBTEarg4gpVjFQBANBkhCq4HC4jVAEA0FyEKrgcZqQKAIBmI1TBpeh47UT1qC7MqQIAoKkIVXApcl7+68pIFQAATUWogotzojpLKgAA0HSEKkiSHA5DR044L/8RqgAAaCpCFSRJx07aVeMwJEmRzKkCAKDJCFWQdPrSX7cQq6z+fFsAANBU/PWEpDMmqXPpDwCAZiFUQdKZa1Rx6Q8AgOYgVEGSdPTUJPXIUEaqAABoDkIVJEnF5XZJUkQoI1UAADQHoQqSpOJTI1XdQ6xebgkAAL6JUAVJ0tHyU6GKkSoAAJqFUAVJ0jFnqAohVAEA0ByEKkiSjp6onVPFSBUAAM1DqIKkM0eqmFMFAEBzEKogwzBcSypw+Q8AgOYhVEEn7TWqrHZI4vIfAADNRaiCa42qQH8/hQb6e7k1AAD4JkIVTq9RFWqVxWLxcmsAAPBNhCqomOUUAABoMUIVmKQOAIAJCFXQsXLnGlUspwAAQHMRqsBIFQAAJiBUwbXwZzcW/gQAoNkIVVBpRbUkKTyYUAUAQHP5TKgqLi5WRkaGwsPDFR4eroyMDB07dqzRYywWS70ff/nLX1xlRo4cWef1u+++u5XPpn0pOVk7p4pQBQBA8wV4uwGeuueee3TgwAEtXrxYkvTwww8rIyNDH3/8cYPH5Ofnu21/+umnGj9+vH7605+67Z8wYYKmTZvm2g4ODjax5e0foQoAgJbziVD17bffavHixVq9erWGDBkiSXrjjTeUlpamnTt3qn///vUeFxsb67b94Ycf6pprrlGfPn3c9oeEhNQp25mUngpVYUGEKgAAmssnLv/l5uYqPDzcFagkaejQoQoPD9eqVas8quPQoUP65JNPNH78+Dqvvfvuu4qKitKgQYP029/+VmVlZaa13Rc4R6rCGKkCAKDZfGKkqqCgQNHR0XX2R0dHq6CgwKM63nrrLXXt2lV33nmn2/57771XiYmJio2N1bZt25SVlaUtW7YoJyenwboqKytVWVnp2i4tLZUk2e122e12j9rjCYfhcNs2s+4zlVbU1htqtbTae/gaZz/QH22Lfm979Ll30O/ecXa/m93/Xg1VU6ZM0dSpUxsts27dOkmq95l0hmF4/Ky6f/7zn7r33nsVFBTktn/ChAmuzwcPHqx+/fopNTVVGzdu1GWXXVZvXTNmzKi33UuWLFFISIhH7fFEUZGfzhxMzM7ONq1up2qHVGGv/TZYvXKZQnwiZredxsI1Wg/93vboc++g373D2e/l5eWm1uvVP6GPPPLIOe+06927t7Zu3apDhw7Vee3w4cOKiYk55/usXLlSO3fu1IIFC85Z9rLLLpPVatXu3bsbDFVZWVnKzMx0bZeWliohIUHp6ekKCws753t4akHheu0qOeraHj16tGl1Ox0uq5TWrJDFIt15843y8+OBylLt/15ycnI0atQoWa1cFm0r9Hvbo8+9g373jrP73XmlySxeDVVRUVGKioo6Z7m0tDSVlJRo7dq1uuKKKyRJa9asUUlJiYYNG3bO4+fMmaOUlBRdcskl5yy7fft22e12xcXFNVjGZrPJZrPV2W+1Wk394fCzuE95a40fvPLq2suYXW0BstlYUf1sZn9N4Rn6ve3R595Bv3uHs9/N7nufmKg+cOBA3XDDDZowYYJWr16t1atXa8KECbr55pvd7vwbMGCAFi1a5HZsaWmp/vWvf+mhhx6qU++ePXs0bdo0rV+/Xvv27VN2drbGjBmj5ORkDR8+vNXPqz1wLafAauoAALSIT4QqqfYOvaSkJKWnpys9PV0XX3yx3nnnHbcyO3fuVElJidu++fPnyzAM/exnP6tTZ2BgoD7//HNdf/316t+/vx599FGlp6dr6dKl8vf3b9XzaS+ck9RZowoAgJbxmWnJERERmjdvXqNlDMOos+/hhx/Www8/XG/5hIQErVixwpT2+SrWqAIAwBw+M1KF1sFq6gAAmINQ1cmVEqoAADAFoaqTYzV1AADMQajq5Lj8BwCAOQhVnVzpyWpJUliQz9yzAABAu0So6uSOV9aGqq7c/QcAQIsQqjq5slOhqouNkSoAAFqCUNXJHT+1+GcXLv8BANAihKpO7jgjVQAAmIJQ1cmVVTgnqjOnCgCAliBUdWI1DkPlVTWSuPwHAEBLEao6MeelP0kKtXWOB0gDANBaCFWdmDNUBQb4yRZAqAIAoCUIVZ3Y8VPzqboySR0AgBYjVHVixytZTgEAALMQqjqx0gqWUwAAwCyEqk7sOKEKAADTEKo6MZ77BwCAeQhVnZhrojpzqgAAaDFCVSfGw5QBADAPoaoTc82pYqQKAIAWI1R1YmUVp5ZUYKQKAIAWI1R1YqcnqhOqAABoKUJVJ3acOVUAAJiGUNWJOUNVKKEKAIAWI1R1YuWVNZIYqQIAwAyEqk7sRFXtSFVIoL+XWwIAgO8jVHVi5VW1I1Vc/gMAoOUIVZ3YiUpGqgAAMAuhqpOqrnGostohSQoNZKQKAICWIlR1UuX2GtfnITZGqgAAaClCVSflvPMvwM+iQH++DQAAaCn+mnZSZ975Z7FYvNwaAAB8H6Gqk3KOVHHnHwAA5iBUdVKsUQUAgLl8JlQ9++yzGjZsmEJCQtStWzePjjEMQ1OmTFF8fLyCg4M1cuRIbd++3a1MZWWlfv3rXysqKkqhoaG69dZbdeDAgVY4g/alvIpH1AAAYCafCVVVVVUaM2aMfvnLX3p8zPPPP6+ZM2fqlVde0bp16xQbG6tRo0aprKzMVWby5MlatGiR5s+fr6+++krHjx/XzTffrJqamkZq9n0nTl3+Y6QKAABz+MwwxdSpUyVJc+fO9ai8YRiaNWuWnnrqKd15552SpLfeeksxMTF677339P/+3/9TSUmJ5syZo3feeUfXXXedJGnevHlKSEjQ0qVLdf3117fKubQHrpEq1qgCAMAUPjNS1VR79+5VQUGB0tPTXftsNptGjBihVatWSZI2bNggu93uViY+Pl6DBw92lemoXCNVXP4DAMAUHfYvakFBgSQpJibGbX9MTIx++OEHV5nAwEB17969Thnn8fWprKxUZWWla7ukpESSdPToUdntdlPaL0mV5WVyVJa7to8cOWJa3YeLjspRWS5L1QlT6+0o7Ha7ysvLdeTIEVmtVm83p9Og39sefe4d9Lt3nN3vzulAhmGYUr9XQ9WUKVNcl/Uasm7dOqWmpjb7Pc5eg8kwjHOuy3SuMjNmzKi33YmJic1rpIeiZplf58unPgAA6KzKysoUHh7e4nq8GqoeeeQR3X333Y2W6d27d7Pqjo2NlVQ7GhUXF+faX1hY6Bq9io2NVVVVlYqLi91GqwoLCzVs2LAG687KylJmZqZr2+Fw6OjRo4qMjDR9Ic3S0lIlJCRo//79CgsLM7VuNIx+9w76ve3R595Bv3vH2f1uGIbKysoUHx9vSv1eDVVRUVGKiopqlboTExMVGxurnJwcJScnS6q9g3DFihV67rnnJEkpKSmyWq3KycnR2LFjJUn5+fnatm2bnn/++Qbrttlsstlsbvs8XeahucLCwvjB8wL63Tvo97ZHn3sH/e4dZ/a7GSNUTj4zpyovL09Hjx5VXl6eampqtHnzZknSBRdcoC5dukiSBgwYoBkzZuiOO+6QxWLR5MmTNX36dPXr10/9+vXT9OnTFRISonvuuUdSbUeOHz9ev/nNbxQZGamIiAj99re/VVJSkutuQAAAAE/4TKj6wx/+oLfeesu17Rx9WrZsmUaOHClJ2rlzp2vSuCQ98cQTOnnypCZNmqTi4mINGTJES5YsUdeuXV1lXnjhBQUEBGjs2LE6efKkrr32Ws2dO1f+/qzfBAAAPOczoWru3LnnXKPq7Nn7FotFU6ZM0ZQpUxo8JigoSC+//LJefrl9Tte22Wz64x//WOdyI1oX/e4d9Hvbo8+9g373jtbud4th1n2EAAAAnViHXfwTAACgLRGqAAAATECoAgAAMAGhqh177bXXlJiYqKCgIKWkpGjlypXeblKHMmPGDF1++eXq2rWroqOjdfvtt2vnzp1uZQzD0JQpUxQfH6/g4GCNHDlS27dv91KLO54ZM2a4lj9xos9bz8GDB3XfffcpMjJSISEhuvTSS7VhwwbX6/S9uaqrq/X0008rMTFRwcHB6tOnj6ZNmyaHw+EqQ5+33JdffqlbbrlF8fHxslgs+ve//+32uid9XFlZqV//+teKiopSaGiobr31Vh04cKDpjTHQLs2fP9+wWq3GG2+8YezYscN47LHHjNDQUOOHH37wdtM6jOuvv9548803jW3bthmbN282brrpJqNnz57G8ePHXWX+/Oc/G127djUWLlxofPPNN8a4ceOMuLg4o7S01Ist7xjWrl1r9O7d27j44ouNxx57zLWfPm8dR48eNXr16mX8/Oc/N9asWWPs3bvXWLp0qfHdd9+5ytD35vrf//1fIzIy0vjPf/5j7N271/jXv/5ldOnSxZg1a5arDH3ectnZ2cZTTz1lLFy40JBkLFq0yO11T/p44sSJRo8ePYycnBxj48aNxjXXXGNccsklRnV1dZPaQqhqp6644gpj4sSJbvsGDBhgPPnkk15qUcdXWFhoSDJWrFhhGIZhOBwOIzY21vjzn//sKlNRUWGEh4cbf/vb37zVzA6hrKzM6Nevn5GTk2OMGDHCFaro89bzu9/9zrjyyisbfJ2+N99NN91kPPjgg2777rzzTuO+++4zDIM+bw1nhypP+vjYsWOG1Wo15s+f7ypz8OBBw8/Pz1i8eHGT3p/Lf+1QVVWVNmzYoPT0dLf96enpWrVqlZda1fE5F46NiIiQJO3du1cFBQVuXwebzaYRI0bwdWihX/3qV7rpppvqPLmAPm89H330kVJTUzVmzBhFR0crOTlZb7zxhut1+t58V155pT7//HPt2rVLkrRlyxZ99dVXGj16tCT6vC140scbNmyQ3W53KxMfH6/Bgwc3+evgM4t/diZFRUWqqalxPfjZKSYmRgUFBV5qVcdmGIYyMzN15ZVXavDgwZLk6uv6vg4//PBDm7exo5g/f742btyodevW1XmNPm8933//vV5//XVlZmbq97//vdauXatHH31UNptN999/P33fCn73u9+ppKREAwYMkL+/v2pqavTss8/qZz/7mSS+39uCJ31cUFCgwMBAde/evU6Zpv7NJVS1YxaLxW3bMIw6+2CORx55RFu3btVXX31V5zW+DubZv3+/HnvsMS1ZskRBQUENlqPPzedwOJSamqrp06dLqn3U1/bt2/X666/r/vvvd5Wj782zYMECzZs3T++9954GDRqkzZs3a/LkyYqPj9cDDzzgKkeft77m9HFzvg5c/muHoqKi5O/vXychFxYW1knbaLlf//rX+uijj7Rs2TKdf/75rv2xsbGSxNfBRBs2bFBhYaFSUlIUEBCggIAArVixQi+99JICAgJc/Uqfmy8uLk4XXXSR276BAwcqLy9PEt/vreF//ud/9OSTT+ruu+9WUlKSMjIy9Pjjj2vGjBmS6PO24Ekfx8bGqqqqSsXFxQ2W8RShqh0KDAxUSkqKcnJy3Pbn5ORo2LBhXmpVx2MYhh555BF98MEH+uKLL5SYmOj2emJiomJjY92+DlVVVVqxYgVfh2a69tpr9c0332jz5s2uj9TUVN17773avHmz+vTpQ5+3kuHDh9dZMmTXrl3q1auXJL7fW0N5ebn8/Nz/zPr7+7uWVKDPW58nfZySkiKr1epWJj8/X9u2bWv616FZ0+vR6pxLKsyZM8fYsWOHMXnyZCM0NNTYt2+ft5vWYfzyl780wsPDjeXLlxv5+fmuj/LycleZP//5z0Z4eLjxwQcfGN98843xs5/9jNudTXbm3X+GQZ+3lrVr1xoBAQHGs88+a+zevdt49913jZCQEGPevHmuMvS9uR544AGjR48eriUVPvjgAyMqKsp44oknXGXo85YrKyszNm3aZGzatMmQZMycOdPYtGmTawkiT/p44sSJxvnnn28sXbrU2Lhxo/GTn/yEJRU6mldffdXo1auXERgYaFx22WWuW/1hDkn1frz55puuMg6Hw/jjH/9oxMbGGjabzbj66quNb775xnuN7oDODlX0eev5+OOPjcGDBxs2m80YMGCAMXv2bLfX6XtzlZaWGo899pjRs2dPIygoyOjTp4/x1FNPGZWVla4y9HnLLVu2rN7f5Q888IBhGJ718cmTJ41HHnnEiIiIMIKDg42bb77ZyMvLa3JbLIZhGM0eVwMAAIAk5lQBAACYglAFAABgAkIVAACACQhVAAAAJiBUAQAAmIBQBQAAYAJCFQAAgAkIVQAAACYgVAEAAJiAUAUAAGACQhUAeKiqqsrbTQDQjgV4uwEA0F6NHDlSgwcPVmBgoN5++20NGjRIK1as8HazALRThCoAaMRbb72lX/7yl/r666/F8+cBNMZi8FsCAOo1cuRIlZSUaNOmTd5uCgAfwJwqAGhEamqqt5sAwEcQqgCgEaGhod5uAgAfQagCAAAwAaEKAADABIQqAAAAE3D3HwAAgAkYqQIAADABoQoAAMAEhCoAAAATEKoAAABMQKgCAAAwAaEKAADABIQqAAAAExCqAAAATECoAgAAMAGhCgAAwASEKgAAABMQqgAAAEzw/wPqpfmMBDWE2gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check the Hamiltonian constraint initially satisfied\n",
    "Ham = get_Ham_diagnostic(initial_state, np.array([0]), my_grid)\n",
    "\n",
    "# plot the profile for Ham\n",
    "plt.plot(r, Ham[0])\n",
    "\n",
    "plt.xlabel('r')\n",
    "#plt.xlim(-0.5,R+2)\n",
    "plt.ylim(-1e-6,1e-6)\n",
    "plt.ylabel('Ham value')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d97150a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 83%|██████████████████████████████████▏      | 834/1000 [03:27<00:41,  4.03‰/s]"
     ]
    }
   ],
   "source": [
    "# for control of time integrator and spatial grid\n",
    "T = 80.0 # Maximum evolution time\n",
    "num_points_t = 256 # time resolution (only for outputs, not for integration, which is decided by python)\n",
    "\n",
    "# Work out dt and time spacing of outputs\n",
    "dt = T/num_points_t\n",
    "t = np.linspace(0, T-dt, num_points_t)\n",
    "eta = 2.0 # the 1+log slicing damping coefficient - of order 1/M_adm of spacetime\n",
    "\n",
    "# Solve for the solution using RK45 integration of the ODE\n",
    "# to make like (older) python odeint method use method='LSODA' instead\n",
    "# use tqdm package to track progress\n",
    "with tqdm(total=1000, unit=\"‰\") as progress_bar:\n",
    "    dense_solution = solve_ivp(get_rhs, [0,T], initial_state, \n",
    "                               args=(my_grid, eta, progress_bar, [0, T/1000]),\n",
    "                        #atol=1e-5, rtol=1e-5,\n",
    "                        max_step=(0.3*my_grid.base_dx), #for stability\n",
    "                        method='RK45', dense_output=True)\n",
    "\n",
    "# Interpolate the solution at the time points requested\n",
    "solution = dense_solution.sol(t).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45fcb464",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot a single point versus time\n",
    "var1 = idx_u\n",
    "var2 = idx_lapse\n",
    "\n",
    "idx = num_ghosts+2 # Choose an inner point\n",
    "r_i = np.round(r[idx],2)\n",
    "var1_of_t = solution[:, var1 * num_points_r + idx]\n",
    "plt.plot(t, var1_of_t, 'b-', label=variable_names[var1])\n",
    "var2_of_t = solution[:, var2 * num_points_r + idx]\n",
    "plt.plot(t, var2_of_t, 'g-', label=variable_names[var2])\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('value at r is '+str(r_i))\n",
    "plt.legend(loc='best')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61ee31bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the profile for some variable at a selection of times\n",
    "#var = idx_lapse # I suggest looking at shiftr, lapse, K or phi to see the gauge evolution\n",
    "var = idx_u\n",
    "\n",
    "for i, t_i in enumerate(t) :\n",
    "    if (i < num_points_t) and (i % 10 == 0) and (t_i > 0.0):\n",
    "        labelt = \"t=\"+str(round(t_i,2))\n",
    "        f_t = solution[i, var * num_points_r: (var + 1) * num_points_r]\n",
    "        plt.plot(r, f_t, label=labelt)\n",
    "\n",
    "#plt.legend(loc=4)\n",
    "plt.xlabel('r')\n",
    "plt.xlim(0.0,60.0)\n",
    "#plt.ylim(-0.04,0.04)\n",
    "plt.ylabel('value over time of ' + variable_names[var])\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94800712",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the diagnostics, just the Hamiltonian constraint for now\n",
    "Ham = get_Ham_diagnostic(solution, t, my_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cd0795a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the profile for Ham at a selection of times\n",
    "\n",
    "for i, t_i in enumerate(t) :\n",
    "    if (i < num_points_t) and (i % 20 == 0) :\n",
    "        labelt = \"t=\"+str(round(t_i,2))\n",
    "        Ham_t = Ham[i]\n",
    "        plt.plot(r, Ham_t, label=labelt)\n",
    "\n",
    "plt.legend(loc=4)\n",
    "plt.xlabel('r')\n",
    "plt.xlim(-1,10)\n",
    "plt.ylim(-0.05,0.05)\n",
    "plt.ylabel('Ham value over time')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94615ed5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the diagnostic for the flux\n",
    "Flux_diagnostic = get_Flux_diagnostic(solution, t, my_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0edd0cbb-db91-4cae-9b15-5730bedce408",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the profile for Flux at a selection of times\n",
    "\n",
    "for i, t_i in enumerate(t) :\n",
    "    if (i < num_points_t) and (i % 5 == 0) :\n",
    "        labelt = \"t=\"+str(round(t_i,2))\n",
    "        Flux_t = Flux_diagnostic[i]\n",
    "        plt.plot(r, Flux_t, label=labelt)\n",
    "\n",
    "#plt.legend(loc=4)\n",
    "plt.xlabel('r')\n",
    "plt.xlim(-1,60)\n",
    "#plt.ylim(-0.05,0.05)\n",
    "plt.ylabel('Flux values over time')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d1d2bbd-f566-402e-88ae-576a1d9ba7eb",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}