diff --git a/examples/BHEvolutionWithSF.ipynb b/examples/BHEvolutionWithSF.ipynb
new file mode 100644
index 0000000..3b211d5
--- /dev/null
+++ b/examples/BHEvolutionWithSF.ipynb
@@ -0,0 +1,340 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ff3f7338",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<script>Jupyter.notebook.kernel.restart()</script>"
+      ],
+      "text/plain": [
+       "<IPython.core.display.HTML object>"
+      ]
+     },
+     "execution_count": 9,
+     "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": 1,
+   "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.bhinitialconditions import *       # go here to change the initial conditions\n",
+    "from source.hamdiagnostic import *             # go here to change the Ham constraint diagnostic\n",
+    "from source.Grid import *"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "e0b8dc67",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Base dx is  0.03366381638499971\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEGCAYAAAB1iW6ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAApUUlEQVR4nO3de3xU9Z3/8ddnJpNMLhAugYgkEC9YREjVgNVCK9a2q1Sxtdbq1rbWdtnaumv7q7/V9vdYd9v97cXHr7/aZallrbba37pQd6312nVbm9QbioCAAqKICAHknsuQTDKX7++PmYRcJskMuczt/Xw85jHnfM/3nHy+TPjMN9/zPeeYcw4REcl+nnQHICIiI0MJXUQkRyihi4jkCCV0EZEcoYQuIpIjCtL1gysqKlxNTU1K+xw/fpzS0tLRCSiDqd35Re3OL6m2e/369Yedc1MSbUtbQq+pqWHdunUp7dPQ0MDixYtHJ6AMpnbnF7U7v6TabjN7b6BtGnIREckRQyZ0M/Ob2Voz22RmW8zs+wnqLDazZjPbGH/dOTrhiojIQJIZcukAPuacC5iZD3jBzH7rnHu5T73nnXNXjHyIIiKSjCETuovdGyAQX/XFX7pfgIiMiVAoRGNjI8FgMN2hjIry8nK2bdvWr9zv91NVVYXP50v6WJbMvVzMzAusB84EfuKcu73P9sXAI0AjsA+4zTm3JcFxlgHLACorK+tWr16ddKAAgUCAsrKylPbJBWp3flG7eysrK6OyspLy8nLMLA2Rja5IJILX6+1V5pyjubmZAwcOEAgEem275JJL1jvn5ic8mHMu6RcwAagH5vYpHw+UxZeXAG8Pday6ujqXqvr6+pT3yQVqd35Ru3vbunWri0ajYxvMGGppaUlYHo1G3datW/uVA+vcAHk1pVkuzrkmoAG4rE95i3MuEF9+GvCZWUUqxxYRGUgu9syHcjJtHnIM3cymACHnXJOZFQMfB+7qU+cU4IBzzpnZBcRmzxxJORrpFok6fvHiu7S0h2jc08kHF3QysbQw3WGJSAZLZpbLNODB+Di6B3jYOfekmX0dwDm3ErgGuNnMwkA7cF38TwM5Sdv2t/C/nzpxomTBlve5/oIZaYxIRPrqukCyoqL3gMTjjz/O1q1bueOOO8Y0nmRmuWwGzktQvrLH8gpgxciGlt/aOiMA3POF8/nGQxs43hFOc0QikqylS5eydOnSMf+5ulI0QwVDsYQ+oTg2ZakjHE1nOCJ5bdeuXcyePZsvf/nL1NbWcs0119DW1gbAv/zLv3D++eczb9483nzzTQAeeOABbrnlljGPM233cpHBdSX08cU+rMe6SD77/hNb2LqvZUSPOefU8fzNlecMWW/79u3cf//9LFy4kJtuuol77rkHgIqKCjZs2MA999zDD3/4Q+67774RjS8V6qFnqGC8R+73efB5lNBF0q26upqFCxcCcMMNN/DCCy8AcPXVVwNQV1fHrl270hUeoB56xupK4EUFXnxeCIY05CKSTE96tPSdRti1XlRUBIDX6yUcTu+5LvXQM1RHPKH7fV4KPaYeukia7d69mzVr1gCwatUqFi1alOaI+lNCz1AdPYZcCr0nhmBEJD3OPvtsHnzwQWprazl69Cg333xzukPqR0MuGSrYo4euMXSR9PN4PKxcubJXWc8x8/nz59PQ0ADAjTfeyI033jh2wcWph56hgqEoHoMCj1Ho1ZCLiAxNCT1DBUMR/D4vZkahFzp0UlQkbWpqanjjjTfSHcaQlNAzVDAcS+hA7KRoWD10ERmcEnqGCoai+AtiH09s2qISuogMTgk9Q3UNuQAUeqBdCV1EhqCEnqE6wlGK4gnd5zVdWCQiQ1JCz1CxHnrs4ynUtEWRtNq1axdz585NdxhDUkLPUB2hKEXxMfRCr2mWi0iyNj8Md8+Fv50Qe9/88Jj82L6X/afjNgC6sChDBcMRJsWfUOTzQGckSiTq8Hry71FcIknb/DA88ZcQao+tN++JrQPUXjusQ0ciEf7sz/6Ml156ienTp/PYY49x+eWX8+EPf5gXX3yRpUuX8sQTT/Ra/853vjPMBqVGCT1DBUMR/AXxk6LxB4J3hCOUFOojkzz22zvg/dcH3t74KkQ6epeF2uGxW2D9g4n3OWUeXP5PQ/7ot99+m1WrVvGzn/2Ma6+9lkceeQSApqYm/vjHPwLwxBNP9FofaxpyyVDBULTHGLp1l4nIIPom86HKU3Daaadx7rnnAr1vlfv5z3++V72+62NJ3b0M1WvaovdEmUheG6onfffc2DBLX+XV8JWnhvWju26TC7Fb5ba3x4Z1SktLe9Xruz6W1EPPUB3haHdC93m7euhK6CKDuvRO8BX3LvMVx8rzwJAJ3cz8ZrbWzDaZ2RYz+36COmZmy81sh5ltNrPzRyfc/BEMRSjqMW0xVqYhF5FB1V4LVy6P9cix2PuVy4d9QjRbJDPk0gF8zDkXMDMf8IKZ/dY593KPOpcDs+KvDwE/jb/LSXDOxS4s6nNSVPdzEUlC7bUjnsD73pzrtttuS1iv6/a56TJkD93FBOKrvvjL9al2FfDLeN2XgQlmNm1kQ80fPR9uAeDzaMhFRIaW1ElRM/MC64EzgZ84517pU2U60PNMRGO8bH+f4ywDlgFUVlam/G0WCATS/g04FgKdse/Lxl3v0tDQSKSzHTBe3bCRzj35cx47Xz7vvtTu3srLy2ltbR37gMZIJBIZsH3BYDCl34WksoNzLgKca2YTgEfNbK5zrufNgRNd7dK3F49z7l7gXoD58+e7xYsXJx0oxP6cSXWfbPR+cxD+8Cxzz/4Aiz80gz1P/AFoZ9bsc1g8L3/+8MmXz7svtbu3bdu2MW7cuLEPaIy0trYO2D6/3895552X9LFSmuXinGsCGoDL+mxqBKp7rFcB+1I5tpxw4vFzXUMuvctFRBJJZpbLlHjPHDMrBj4OvNmn2uPAl+KzXS4Emp1z+5GTcmIMve88dM1yEZGBJTPkMg14MD6O7gEeds49aWZfB3DOrQSeBpYAO4A24CujFG9e6NtDL9RJURFJwpAJ3Tm3Geg3iBNP5F3LDvjmyIaWv7oSt6YtimSGsrIyAoHA0BXTTFeKZqBgv2mL8XINuYjIIJTQM1DfHrqZUVTgoUNDLiJpFQgEuPTSSzn//POZN28ejz32GBB7AMbs2bP58pe/TG1tLddccw1tbW0A3HHHHcyZM4fa2truC5IOHTrEZz/7WRYsWMDFF1/Miy++OCLx5c+k5ixyYgzd213m93k1hi557661d/Hm0b5zMoZn9qTZ3H7B7UnV9fv9PProo4wfP57Dhw9z4YUXsnTpUgC2b9/O/fffz8KFC7npppu45557uOmmm3j00Ud58803MTOampoAuPXWW/n2t7/NokWL2Lp1K5/97GfZtm3bsNuihJ6Bup5O1DXk0rWsIReR9HLO8b3vfY/nnnsOj8fD3r17OXDgAADV1dUsXLgQgBtuuIHly5fzrW99C7/fz9e+9jU+9alPccUVVwDw+9//nq1btwIQjUZpaWkZdD56spTQM1BHeIAeuk6KSp5Ltic9Wh566CEOHTrE+vXr8fl81NTUEAwGgdjQaE9mRkFBAWvXruXZZ59l9erVrFixgj/84Q9Eo1HWrFlDcXHxiCTyLhpDz0DBUO956AD+Ag25iKRbc3MzU6dOxefzUV9fz3vvvde9bffu3axZswaAVatWsWjRIgKBAM3NzSxZsoQf//jHbNy4EYBPfvKTrFixonvfrvLhUkLPQCdOimrIRSSTfOELX2DdunXMnz+fhx56iNmzZ3dvO/vss3nwwQepra3l6NGj3HzzzbS2tnLFFVdQW1vLxRdfzN133w3A8uXLWbduHbW1tSxYsICVK1cO9CNToiGXDBQMR/B6DJ/3REIv8nlpVw9dJC265qBXVFR098J72rVrFx6Pp19iLikpYe3atf3qV1RU8Ktf/QoY/F4uqVIPPQMFQ1H8Bb0/Gr/Pq2mLIjIoJfQM1PN5ol2KNeQikrH6PgAjXZTQM1AwFO2X0DXLRfJZ7O4i+eVk2qyEnoE6wieeJ9pFs1wkX/n9fo4cOZJXSd05x5EjR/D7/Sntp5OiGSg2ht63h64hF8lPVVVVNDY2cujQoXSHMiqCwWDCxO33+6mqqkrpWEroGShhD12X/kue8vl8nHbaaekOY9Q0NDSk9FSiwWjIJQMFQ5F+PfQin5eOcDSv/uwUkdQooWeg2EnRvj302HrX04xERPpSQs9AiaYtdvXYNewiIgNRQs9AwXCChO7rSujqoYtIYkroGahjkCEX9dBFZCBDJnQzqzazejPbZmZbzOzWBHUWm1mzmW2Mv+4cnXDzQzAU6X5aUZfuHrouLhKRASQzbTEMfMc5t8HMxgHrzex3zrmtfeo975y7YuRDzD/BcP8rRYvj6+2dSugiktiQPXTn3H7n3Ib4ciuwDZg+2oHlq2jU0RmO9rp1LtA9L11j6CIykJTG0M2sBjgPeCXB5ovMbJOZ/dbMzhmJ4PJR17TEAU+KashFRAaQ9JWiZlYGPAJ8yznX0mfzBmCmcy5gZkuA3wCzEhxjGbAMoLKykoaGhpSCDQQCKe+TbQKdsQuHGnftpKFhT6wsEGD3xg0ArH9tM7Y/Py7wzYfPOxG1O7+MaLudc0O+AB/wDPA/kqy/C6gYrE5dXZ1LVX19fcr7ZJt9TW1u5u1Pun9/5b3usvr6evfOwVY38/Yn3aMbGtMY3djKh887EbU7v6TabmCdGyCvJjPLxYD7gW3OuR8NUOeUeD3M7AJiQzlHRuD7Ju90dD9PtP+9XEDTFkVkYMn87b4Q+CLwupltjJd9D5gB4JxbCVwD3GxmYaAduC7+TSIp6hoj73+3RSV0ERnckAndOfcCYEPUWQGsGKyOJOd4RzyhF/a/fS7EpjSKiCSiK0UzTFNbJwATSwp7leteLiIyFCX0DHOsLQTAxBJfr3KPxyj06iEXIjIwJfQM09VDn9Cnhw5QXOilrTM81iGJSJZQQs8wTW0hvB5jvL//6Y3JpYUcPd6ZhqhEJBsooWeYY22dTCj2EZ8F2suk0kKOBJTQRSQxJfQM09QWYkKf8fMuk8vUQxeRgSmhZ5hjbZ39Zrh0mVRaxJHjHWMckYhkCyX0DHOsLZTwhCjExtCPtYWIRnXNloj0p4SeYZraOgcdcolEHc3toTGOSkSygRJ6hokNuSRO6JNKYz13DbuISCJK6BkkGIoQDEUHGXIpAtBMFxFJSAk9gxwb4LL/LpPLYuWa6SIiiSihZ5BjxxNf9t9lcnzI5bASuogkoISeQQa77B9gYjyhH9WQi4gkoISeQbpvzFWauIfu83ooL/bppKiIJKSEnkGGGkOH2LDLEQ25iEgCSugZ5MSQS+IeOsSmLmrIRUQSUULPIMfaQpQUeinq8/i5niaXFWrIRUQSUkLPIIPdx6XLpNIiTVsUkYSU0DNI8yB3WuxSEb/jou7nIiJ9DZnQzazazOrNbJuZbTGzWxPUMTNbbmY7zGyzmZ0/OuHmtmOD3Mely6TSQqIOmnQ/FxHpI5keehj4jnPubOBC4JtmNqdPncuBWfHXMuCnIxplnmga5E6LXbru53JU4+gi0seQCd05t985tyG+3ApsA6b3qXYV8EsX8zIwwcymjXi0OW6wG3N1qSiL3c/lsGa6iEgf/R9cOQgzqwHOA17ps2k6sKfHemO8bH+f/ZcR68FTWVlJQ0NDSsEGAoGU98kWUedoagvRfGg/DQ1Hem3r2e49rVEAnl/7GsHdKX18WSeXP+/BqN35ZUTb7ZxL6gWUAeuBqxNsewpY1GP9WaBusOPV1dW5VNXX16e8T7Y4drzDzbz9SXff8zv7bevZ7gPN7W7m7U+6X67ZNYbRpUcuf96DUbvzS6rtBta5AfJqUrNczMwHPAI85Jz7dYIqjUB1j/UqYN9Jfsfkpe7L/ocYcum6n8uRgMbQRaS3If9mt9jj5+8HtjnnfjRAtceBW8xsNfAhoNk5t3+AupJA9cRinv+rSygfIqF33c9Fc9FFpK9kBmEXAl8EXjezjfGy7wEzAJxzK4GngSXADqAN+MqIR5rjCrweqieVJFV3cmmhHnIhIv0MmdCdcy8ANkQdB3xzpIKSwenyfxFJRFeKZqGKsiIOtiqhi0hvSuhZaMakEhqPtevyfxHpRQk9C82YXEJnOMr7LcF0hyIiGUQJPQvNnFQKwK4jx9MciYhkEiX0LDRzcmw2zO4jbWmOREQyiRJ6Fjp1QjE+r/HeUSV0ETlBCT0LeT1G1cQS3tOQi4j0oISepWZOLuE9DbmISA9K6Flq5qRYQo9d0yUiooSetWZMLiXQEdY9XUSkmxJ6lqqJz3TRiVER6aKEnqW6pi7qxKiIdFFCz1JVE0swQydGRaSbEnqW8vu8TBvv18VFItJNCT2LzZhcosv/RaSbEnoWq5lcym6dFBWROCX0LDZjcgmHA50EOsLpDkVEMoASehY7bXLsros7DwXSHImIZAIl9Cx2zqnlALyxtyXNkYhIJlBCz2LVk4opL/bx+t6mdIciIhlgyIRuZj83s4Nm9sYA2xebWbOZbYy/7hz5MCURM2Pe9HJe39uc7lBEJAMk00N/ALhsiDrPO+fOjb9+MPywJFlzp5ez/f1WOsKRdIciImk2ZEJ3zj0HHB2DWOQkzJteTiji2P5+a7pDEZE0s2Ruv2pmNcCTzrm5CbYtBh4BGoF9wG3OuS0DHGcZsAygsrKybvXq1SkFGwgEKCsrS2mfXDBYuw+1Rfmfz7Xz5TmFXDLDN8aRjS593vlF7U7OJZdcst45Nz/hRufckC+gBnhjgG3jgbL48hLg7WSOWVdX51JVX1+f8j65YLB2R6NRV/u3z7g7Htk0dgGNEX3e+UXtTg6wzg2QV4c9y8U51+KcC8SXnwZ8ZlYx3ONKcrpOjG5u1IlRkXw37IRuZqeYmcWXL4gf88hwjyvJm1dVzlsHdGJUJN8VDFXBzFYBi4EKM2sE/gbwATjnVgLXADebWRhoB66L/1kgY6TnidHaqgnpDkdE0mTIhO6cu36I7SuAFSMWkaRs3vTYFaObG5uV0EXymK4UzQFVE4uZWOJj456mdIciImmkhJ4DzIwLT5/MSzsOo9EukfylhJ4jPjJrCvuag7xzSA+8EMlXSug54iOzYjNFn3/7UJojEZF0UULPEdWTSqiZXMLzbx9OdygikiZK6DnkI7Om8PLOI3SGo+kORUTSQAk9h3xkVgVtnRE27D6W7lBEJA2U0HPIRWdMxusxjaOL5Ckl9Bwyzu/j/BkTNI4ukqeU0HPMR2ZN4fW9zRwJdKQ7FBEZY0roOebSs6fiHPz2jffTHYqIjDEl9BwzZ9p4Zk0t47GNe9MdioiMMSX0HGNmXHXuqby66xiNx9rSHY6IjCEl9Bx01bnTAXh80740RyIiY0kJPQdVTyqhbuZEHntNCV0knyih56hPn3sq2w+0sm1/S7pDEZExooSeo5bMm4bXY/xGJ0dF8oYSeo6aXFbE4rOm8Mj6vXrWqEieUELPYTcurOFwoIPHN2osXSQfKKHnsEVnVjD7lHHc/8K7epKRSB4YMqGb2c/N7KCZvTHAdjOz5Wa2w8w2m9n5Ix+mnAwz46uLTuPN91t1fxeRPJBMD/0B4LJBtl8OzIq/lgE/HX5YMlKWnnsqU8YVcd8L76Y7FBEZZQVDVXDOPWdmNYNUuQr4pYv9Tf+ymU0ws2nOuf0jFWS+uGvtXbx59M1+5U1NTTz4Xw+e9HEnntHOq0fb+PzjEygp9A4nxDE13HZnK7U7v4xku4dM6EmYDuzpsd4YL+uX0M1sGbFePJWVlTQ0NKT0gwKBQMr7ZJPGo400dTb1K49EIjQ19S9Plt+Bx2DXwWZOLcue0ybDbXe2Urvzy4i22zk35AuoAd4YYNtTwKIe688CdUMds66uzqWqvr4+5X1ywUi0+4fPvOlm3v6ke233sWEfa6zo884vandygHVugLw6Et21RqC6x3oVoHlyGebPLz6DirJC/uGpbZrxIpKjRiKhPw58KT7b5UKg2Wn8POOUFRXw7U+cxdpdR/nvrQfSHY6IjIJkpi2uAtYAHzCzRjP7qpl93cy+Hq/yNLAT2AH8DPjGqEUrw/L5+dWcObWMu377JqFINN3hiMgIS2aWy/VDbHfAN0csIhk1BV4P3718Nl99cB2/ePFdln30jHSHJCIjKHumPMiI+NjsqXxiTiX/97/fYsfBQLrDEZERpISeZ8yMv//MXIoLvdz2H5uIRHWCVCRXKKHnoanj/Pzgqrls3NPEvc/tTHc4IjJClNDz1JW107h87inc/bu39BAMkRyhhJ6nzIy/+/Rcykt8/Pn/W09TW2e6QxKRYVJCz2MVZUWsvKGO95uD/MWq1whrKqNIVlNCz3N1Myfyg6vO4fm3D/N/ntme7nBEZBhG4uZckuWuu2AGW/a18K/P7eSMqWVcO7966J1EJOMooQsAf33FHN49fJw7HtlMaWEBn6qdlu6QRCRFGnIRAAoLPNz7pTrqZk7k1tWv8ew23e9FJNsooUu3ksIC7r9xAXNOHc/ND23g+bcPpTskEUmBErr0Mt7v48GvXMDpFaXc9MCrPL5Jd0IWyRZK6NLPxNJCfrXsIs6bMZG/XPUaP3tup+6hLpIFlNAlofISH7+86QKWzDuFv396G99/YqvmqYtkOCV0GZDf52XF9efzlYU1PPDSLv70vlc42BJMd1giMgAldBmUx2P8zZXncPfnP8jrjc0sWf4CL71zON1hiUgCSuiSlM+cV8VjtyykvLiAG+57hR8+s52OcCTdYYlID0rokrSzKsfx+C2LuPr8KlbU7+CK5S+wcU9TusMSkTgldElJaVEBP/zcB/nFVxYQ6Ahz9T0v8g9PbyPQEU53aCJ5TwldTsolH5jKM9/+KJ9fUM29z+3kkh828PC6PUT1BCSRtEkqoZvZZWa23cx2mNkdCbYvNrNmM9sYf9058qFKphnv9/GPV9fy6Dc+TNXEYv7qPzez9CcvsOadI+kOTSQvDZnQzcwL/AS4HJgDXG9mcxJUfd45d2789YMRjlMy2HkzJvLrmz/MP193LkcCnVz/s5e5/t6XeXmnErvIWEqmh34BsMM5t9M51wmsBq4a3bAk25gZV507nfrbFvPXV8xhx6EA1937Mtfdu4YXdxzWlaYiY8CG+o9mZtcAlznnvhZf/yLwIefcLT3qLAYeARqBfcBtzrktCY61DFgGUFlZWbd69eqUgg0EApSVlaW0Ty7IxnZ3RhwNe8I89W6I5g5HVZnxiRofF00roNBrSR0jG9s9EtTu/JJquy+55JL1zrn5ibYlcz/0RP/7+n4LbABmOucCZrYE+A0wq99Ozt0L3Aswf/58t3jx4iR+/AkNDQ2kuk8uyNZ2fxK4MxTh8U37+PkL7/KLN1p57F24bkE1186vpqaidND9s7Xdw6V255eRbHcyCb0R6PkImypivfBuzrmWHstPm9k9ZlbhnNMlhXnO7/Ny7fxqPldXxcs7j/LzF99l5R/f4Z6Gd7jgtElcO7+aJfNOoaRQz1oRGa5k/he9Cswys9OAvcB1wJ/2rGBmpwAHnHPOzC4gNjavM2LSzcy46IzJXHTGZA60BHlkQyP/sa6R2/5jE3c+9gYfP7uSK2qn8dGzpuD3edMdrkhWGjKhO+fCZnYL8AzgBX7unNtiZl+Pb18JXAPcbGZhoB24zuksmAygcryfbyw+k5svPoN17x3j1xsa+a833ufxTfsoKyrgE3Mq+ZNzKiGsXyGRVCT1d65z7mng6T5lK3ssrwBWjGxokuvMjAU1k1hQM4kfXDWXNe8c4cnN+3hmywEefW0vBQYf3r2Wj589lUs+MJXqSSXpDlkko2ngUjKCz+vho2dN4aNnTeHvPxNl3a5jPPC7dbx9tI07H9sCbOG0ilIWnjmZRWdWcNHpFZSX+NIdtkhGUUKXjOPzerjojMl07Cli8eLF7DwUoGH7IV565zCPbtjLv728GzM459Tx3T38+TUTmTrOn+7QRdJKCV0y3ulTyjh9Shk3LTqNUCTK5sYmnn/7MK/sPMqqtbv5xYu7AKiZXML5MydyXvUEPlg9gdmnjKewQLcrkvyhhC5Zxef1UDdzEnUzJwHQGY6yZV8zr+46ytp3j/HcW4f49Ya9ABR6Pcw5dTy1VeWcc+p45kwrZ1ZlmWbRSM5SQpesVljg4bwZEzlvxkSWfRScc+xrDrJpTxOb9jTx2p4mfr1hL79c8x4AXo9x5pQyzp42jg+cMp6zKss4q3Ic0ycU4/EkdwWrSKZSQpecYmZMn1DM9AnFLJk3DYBo1LHnWBtb97WwZV8LW/e38PLOo/xm44nr40oKvcyaWsYZU8s4Y0oZp1eUcvqUMmZOLlGPXrKGErrkPI/HmDm5lJmTS7k8nuQBmttD7DjYylsHArx1oJW3DrTy0o4j3UM2AGYwfUIxp1WUMmNSCTMmlTBzcgkzJpUyc3IJpUX6LySZQ7+NkrfKi329xuO7HO8I8+7h47xzKMC7h4+z89Bx3jtynKde309TW6hX3Umlhd1/EVRNLGb6xNjy9InFTCsvZmKJDzMN5cjYUEIX6aO0qIC508uZO72837bm9hC7j7Sx+2gb7x09TuOxdhqPtfP2wVbqtx+kIxztVb+owMO0cj+nlPs5tbyYU8r9VI73Uzm+iCnjut6LKCrQsI4MnxK6SArKi33MqypnXlX/ZO+c48jxThqPtbOvqZ33m4Psb25nf3OQ95uDvPLuUQ60BAkneEzfxBIfU8f5qRhXSOR4kOcDW6koK6KirJCKcUVMKStiUmkhk0oLNaYvA1JCFxkhZhZPwkWcWz0hYZ1o1HG0rZMDLUEOtnZwsCXIwZYODrTG3g8HOtjTFGXTK7tpD0USHqO00MukskImlRYxOZ7kJ5b4mFBSyMSSHsulPiaWFDKhxKe/APKEErrIGPJ4TiT9cwao03V/7OMdYQ4HYkn+UGsnR493cvR4B0eOd3LseCdHjse+GLbtb+FYWyfBUHSAI4Lf56G82NfrNb7nst/HOH8B43ssd5WX+QvwakpnVlBCF8lQpUUFlBYVMHPy4A8C6RIMRTjW1smx4yGa2jo51hbiWFsnze0hmttjZV3Le5uCbNvfSnN7iEBHeOhYCr2U+QsoKyqgzO9jXFFseZy/oLu8K97SQi+lRSfKyoq8lBTGlksKvfi8unp3tCihi+QIv8/LtPLY7JpUhCNRAh1hWtrDtARDsVd7mNZgiJZg7D0QDBPoCNPaESYQLzvYGowtd8S2JXvD7EKvh5IiLyU+LyXxL4DiQi/trUH+c98GSgpjXwB+nze+HNte7OtaLqDYF1v3+zz4fbHt/nhZPv81oYQukucKvB4mlBQyoaTwpI/hnKM9FCHQEeZ4R4TjHeHYqzNMoMd6W2ck/orVaw/F3ts6wzR1OFr2tdAe394eihCKpH5P/EKvhyKfJ57wTyR9f4GXoq5lnxd/Qayev+BEvaJ4naKC+HK8TvdygZfCgtj2nu+FBR4KvZ60T1FVQheRYTMzSgoLYo8SHHdyx0j0bM1QJEp7KBJP8rH39lCY9s5YeTAU6X4PhiK9yzsjdIQjBENRguFYWWswzKHWDjrD0dg+Xe+hCAkmH6WssMBDUfwLpdDbI9nHE77P2/uLwOf1cPRQB7879np3na7y2Mso6rk+xM3mlNBFJGN1JbLx/tG/9304EqUjnuA7wtH4K0JHqPdyMByhMxyNvSJROkJd7xE6IrHyjvCJ91C8Xlf91mCYI/HlUCRK6/EIW5vej9WPl53s896U0EVEiA09FXg9Y347h0R/mUSijlD8CyYUT/KhsKMzEmHWXQMfSwldRCTDeD2G1+NN+SIyzR8SEckRSfXQzewy4J8BL3Cfc+6f+my3+PYlQBtwo3NuwwjHmts2PwzP/gCa94B5wUW63z+KBxqi/cq734vjN5dqP5qg7BgUTxyZ5fIqmPVJePu/oblx8G0ns37pnbHjPfsDaG7kwqIKCCwdtM6IlNVe2+Pff5TKe33Gg+8z9dTPAYsH36ff700KPy+ZfZOtMxLHiNe5cM33oOHwwHWSPdZo1BvNuidTfwDmhhh9NzMv8BbwCaAReBW43jm3tUedJcBfEEvoHwL+2Tn3ocGOO3/+fLdu3bqUgk001pQTNj8MT/wlhNrTHUl6eXyx+9VGOlOrM5wyXzF88E9h07/3/vcfyfIrl8eW+37GA+wT8RTh/fSKgfe5cvmJL4KBtg9n355fCKN9jGTrpLPeaNYFtv7q+8zZ8dOk608p9bx76Hj09H4bSC6hXwT8rXPuT+Lr3wVwzv1jjzr/CjQ451bF17cDi51z+wc6rhJ6D3fPjfXMJTd5i2LvkY6R2cdbBFULoPHVgbcPZ9+qBbHloeqMxDGSrZPOeqNZF4jsfgWvCyVd/4Pf/a/OTe9HivrvkNyQy3SgZ7ZpJNYLH6rOdKBXQjezZcAygMrKShoaGpL48ScEAoGU98kGFzc3kr/XtqWfg4T//iNWHv/PPWL7RDpobmqiPNIxKvs2NzUBDFlnJI6RbJ101hvNugDliZL5IPV9Hga8AiyZhD7Q72GqdXDO3QvcC7Eeeqq97Zztob9WpR56GlnXeYfRKi+vji0k+IxPap/yaiZ8+8UB/7Ib7r4Tvv1ibGWoOiNxjGTrpLPeaNYFgv94Jv6OQ0nXD93lHXBMMplZLo1AdY/1KmDfSdSRgVx6Z2zMLN95fOAd4vLzRHWGU+Yrhrob+//7j2T5pXcm/owH2CfiKRp8n66Tu4NtH86+XcbiGMnWSWe90awL7Dz9iynV39fq9ibcQHIJ/VVglpmdZmaFwHXA433qPA58yWIuBJoHGz+XPmqvjZ0A6epZmbfXe7TrY+pT3v1ePOnErJZ+ZTZyy+XVMP+r8TiH2HYy65++B676SXdZsGjKkHWGXXblcrjiRz3+/UehvPbaPp/x4Pts/8A3B9+n60TZYNuHs2/C38tROkaPOsGiKQPXSfZYo1FvNOsCBysvTqn+4TZ3NOEGiN1UZ6gXsdkrbwHvAP8rXvZ14OvxZQN+Et/+OjB/qGPW1dW5VNXX16e8Ty5Qu/OL2p1fUm03sM4NkFeTmofunHsaeLpP2coeyw74ZjLHEhGR0aErRUVEcoQSuohIjlBCFxHJEUroIiI5YshL/0ftB5sdAt5LcbcK4PAohJPp1O78onbnl1TbPdM5NyXRhrQl9JNhZuucc/PTHcdYU7vzi9qdX0ay3RpyERHJEUroIiI5ItsS+r3pDiBN1O78onbnlxFrd1aNoYuIyMCyrYcuIiIDUEIXEckRWZXQzexzZrbFzKJmlvPTm8zsMjPbbmY7zOyOdMczVszs52Z20MzeSHcsY8XMqs2s3sy2xX/Hb013TGPBzPxmttbMNsXb/f10xzSWzMxrZq+Z2ZMjcbysSujAG8DVwHPpDmS0xR/O/RPgcmAOcL2ZzUlvVGPmAeCydAcxxsLAd5xzZwMXAt/Mk8+7A/iYc+6DwLnAZfFnKuSLW4FtI3WwrErozrltzrnt6Y5jjFwA7HDO7XTOdQKrgavSHNOYcM49Bwx8E/8c5Jzb75zbEF9uJfaffHp6oxp98Vt8B+KrvvgrL2ZqmFkV8CngvpE6ZlYl9Dwz0IO3JceZWQ1wHvBKmkMZE/Fhh43AQeB3zrm8aDfwY+CvgOhIHTDjErqZ/d7M3kjwyoveaQ9JPXhbcouZlQGPAN9yzrWkO56x4JyLOOfOJfYs4gvMbG6aQxp1ZnYFcNA5t34kj5vUE4vGknPu4+mOIUPowdt5xsx8xJL5Q865X6c7nrHmnGsyswZi509y/YT4QmCpmS0B/MB4M/s359wNwzloxvXQpVsyD+eWHGFmBtwPbHPO/Sjd8YwVM5tiZhPiy8XAx4E30xrUGHDOfdc5V+WcqyH2f/sPw03mkGUJ3cw+Y2aNwEXAU2b2TLpjGi3OuTBwC/AMsRNkDzvntqQ3qrFhZquANcAHzKzRzL6a7pjGwELgi8DHzGxj/LUk3UGNgWlAvZltJtaJ+Z1zbkSm8OUjXfovIpIjsqqHLiIiA1NCFxHJEUroIiI5QgldRCRHKKGLiOQIJXQRkRyhhC4yAIvR/xHJGvplFenBzGri9yS/B9hA79sviGQ0XVgk0kP8Toc7gQ87515OczgiKVEPXaS/95TMJRspoYv0dzzdAYicDCV0EZEcoYQuIpIjdFJURCRHqIcuIpIjlNBFRHKEErqISI5QQhcRyRFK6CIiOUIJXUQkRyihi4jkiP8Psb7ICiImlaEAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "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 = 124 # 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": 3,
+   "id": "5515ca2e",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEGCAYAAACpXNjrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAATSElEQVR4nO3de5BkZXnH8e/DLCsKEkCWYbnoblJbRKAMyJSiSaUGgYSLCiaFgqVuhBRJlcRLYRmUP4SQVFEmmotBrFWJa8WwodSEDYUaJY6m4g0WNgIixYaLrCwsIoIDcll48keflXad2XnPzNt7Zqa/n6qu6XP6vN3PM8z2j3Pe06cjM5Ekqabdui5AkrT4GC6SpOoMF0lSdYaLJKk6w0WSVN2SrgvYlfbff/9csWLFrMY+9thj7LnnnnULWkCGuX97H87eYbj77+99w4YNP87MZW3GD1W4rFixghtuuGFWYycmJhgfH69b0AIyzP3b+3jXZXRmmPvv7z0i7mk73sNikqTqDBdJUnWGiySpOsNFklSd4SJJqs5wkSRVZ7hIkqozXApdf/82Hpp8susyJGlBMFwKPPzYU1y28UnOWTu7D2BK0rAxXArcsXUSgP97cLLjSiRpYTBcCnxw/a0A/OyJbR1XIkkLg+FS4IV7DNUl2CRpzgyXArtF1xVI0sJiuEiSqjNcCmR2XYEkLSyGiySpOsNFklSd4VLAo2KS1I7hUsJ0kaRWDJcCabpIUiuGSwHPFpOkdgyXAs+aLpLUiuFS4FmzRZJaMVwKpHsuktSK4VLAaJGkdgyXAu64SFI7hksBT0WWpHYMlwLuuUhSO4ZLgWc8XUySWjFcCrjnIkntGC4F/BClJLVjuBQwWiSpHcOlgHsuktROp+ESESdFxO0RsSkiLpji8YiIf2ge/15EvLx0bFVmiyS10lm4RMQIcBlwMnA4cFZEHL7DZicDq5rbucDlLcZW84x7LpLUypIOX/sVwKbMvBMgItYBpwHf79vmNOAz2bu417cjYp+IWA6sKBhbTf9hsXf8y41MPrGNp595lmeeTTKf+5DlzjJoocfTI4/8nI/e9s2uy+iEvQ9n77C4+r/w1Jfy8hfvu8ter8twORi4t295M/DKgm0OLhwLQEScS2+vh9HRUSYmJloX+vjjT/zi/oZN9/P83YMlAbsFRED80utN/zw7eWjeG+EZnph8pOsyOmHvw9k7LK7+N950I4/eOVK8/eTk5KzeL7frMlymeq/d8X/wp9umZGxvZeYaYA3A2NhYjo+Ptyix53nfug6eeIKTjzyQy99yTOvxi8HExASz+d0tBvY+3nUZnRnm/ufae5fhshk4tG/5EOC+wm2WFoytxg/oS1I7XZ4tdj2wKiJWRsRS4Exg/Q7brAfe1pw1dizwSGZuKRxbzfY5F+f1JalMZ3sumbktIs4DvgyMAFdk5q0R8afN4x8HrgVOATYBjwNv39nYgdU6qCeWpEWqy8NiZOa19AKkf93H++4n8I7SsYPiN1FKUjt+Qr/Ahae+FNj5mWCSpOcYLgWet6T89D1JkuEiSRoAw0WSVJ3h0oLz+pJUxnCRJFVnuLTg2WKSVMZwkSRVZ7hIkqozXCRJ1RkukqTqDJcWPBVZksoYLpKk6gyXFjwVWZLKGC6SpOoMF0lSdYaLJKk6w6UFzxaTpDKGiySpOsOlBc8Wk6QyhoskqTrDRZJUneEiSarOcJEkVWe4tOCpyJJUxnCRJFVnuLTgqciSVMZwkSRVZ7hIkqozXCRJ1RkukqTqDJcWPBVZksoYLpKk6gyXFjwVWZLKGC6SpOoMF0lSdYaLJKk6w0WSVJ3hIkmqrpNwiYj9IuIrEXFH83PfabY7KSJuj4hNEXFB3/qLIuJHEbGxuZ2y66qXJM2kKFwi4vkRcVjF170AuC4zVwHXNcs7vuYIcBlwMnA4cFZEHN63yd9m5lHN7dqKtUmS5mjGcImI1wEbgS81y0dFxPo5vu5pwNrm/lrg9Cm2eQWwKTPvzMyngHXNOEnSPLekYJuL6L3RTwBk5saIWDHH1x3NzC3N822JiAOm2OZg4N6+5c3AK/uWz4uItwE3AOdn5sNTvVBEnAucCzA6OsrExETrYm+9fxsADz744KzGLwaTk5P2PoSGuXcY7v7n2ntJuGzLzEei5cfTI+KrwIFTPHRh6VNMsW771b0uBy5pli8BPgycPdWTZOYaYA3A2NhYjo+PF778cx6/eQtsvJFly5YxPn5M6/GLwcTEBLP53S0G9j7edRmdGeb+59p7SbjcEhFvBkYiYhXwTuCbMw3KzBOmeywiHoiI5c1ey3Jg6xSbbQYO7Vs+BLivee4H+p7rE8A1BX1IknaRkgn9PwOOAJ4ErgQeBd49x9ddD6xu7q8Grp5im+uBVRGxMiKWAmc242gCabs3ALfMsR5JUkUz7rlk5uP0DmWVHs4qcSlwVUScA/wQOAMgIg4CPpmZp2Tmtog4D/gyMAJckZm3NuM/FBFH0TssdjfwJxVrkyTN0YzhEhFf47m5jl/IzNfM9kUz8yHg+CnW3wec0rd8LfArpxln5ltn+9qSpMErmXN5b9/9PYA/BLYNphxJ0mJQclhsww6r/icivj6geiRJi0DJYbH9+hZ3A45h6lOMJUkCyg6LbaA35xL0DofdBZwzyKLmm/yVGSdJ0s6UHBZbuSsKkSQtHtOGS0T8wc4GZuYX6pczP7W8OIEkDb2d7bm8biePJTA04SJJamfacMnMt+/KQiRJi0fJhD4RcSq9S8DssX1dZv7FoIqSJC1sJd/n8nHgTfSuMRb0LtXykgHXJUlawEouXPnqzHwb8HBmXgy8il++WvGi56nIktROSbj8vPn5eHNhyacBT0+WJE2rZM7lmojYB/hr4EZ6Z4p9YpBFzTeeiixJ7ZR8iPKS5u7nI+IaYI/MfGSwZUmSFrKSCf3/jYgPRMRvZOaTBoskaSYlcy6vp3dNsasi4vqIeG9EvHjAdUmSFrAZwyUz78nMD2XmMcCbgZfRu3ilJElTKv0Q5QrgjfQ+7/IM8L4B1jTveCqyJLVT8n0u3wF2B64CzsjMOwdelSRpQSvZc1mdmT8YeCXzmKciS1I7JXMuQx0skqT2Ss4WkySpFcNFklRdyYT+CHAqsKJ/+8z8yODKml88W0yS2imZ0P8P4AngZuDZwZYzvzmxL0llSsLlkMx82cArWQDcg5GkMiVzLl+MiN8beCXzmHssktROyZ7Lt4F/i4jd6H2XSwCZmXsPtDJJ0oJVEi4fpvftkzdnemBIkjSzksNidwC3GCySpFIley5bgImI+CLw5PaVnoosSZpOSbjc1dyWNreh5cS+JJUp+Zrji3dFIQuBezCSVKbkE/rL6H1/yxHAHtvXZ+ZrBljXvOIeiyS1UzKh/1ngB8BK4GLgbuD6AdYkSVrgSsLlRZn5KeDpzPx6Zp4NHDvguiRJC1jJhP7Tzc8tEXEqcB9wyOBKmn+ca5GkdkrC5S8j4teA84GPAnsD7xloVfOUcy+SVKbkbLFrmruPAMfVeNGI2A/4V3qX8b8beGNmPjzFdlcArwW2ZuaRbcdLkroxbbhExEeBaQ8IZeY75/C6FwDXZealEXFBs/znU2z3aeAfgc/McnxVHh6TpDI7m9C/AdjQ3F7fd3/7bS5OA9Y299cCp0+1UWZ+A/jJbMfX4uEwSWonSi4ZFhE3ZebR1V404qeZuU/f8sOZue80264ArtnhsFib8ecC5wKMjo4es27dutb1Xn//Ni7b+CRjoyOcd/QeMw9YhCYnJ9lrr726LqMT9j6cvcNw99/f+3HHHbchM8fajC+Z0IedHB6bTkR8FThwiocubPtcc5GZa4A1AGNjYzk+Pt76OR6/eQtsvJFly5YxPn5M5QoXhomJCWbzu1sM7H286zI6M8z9z7X30nBpLTNPmO6xiHggIpZn5paIWA5sbfn0cx3finMtktTOtHMuEfGziHg0Ih4FXrb9/vb1c3zd9cDq5v5q4OpdPH5WnHuRpDLThktmvjAz925uS/ruv7DCt1BeCpwYEXcAJzbLRMRBEXHt9o0i4krgW8BhEbE5Is7Z2XhJ0vwwsMNiO5OZDwHHT7H+PuCUvuWz2owfNA+PSVKZkmuLDT0Ph0lSO4aLJKk6w0WSVJ3hUsC5Fklqx3BpwbkXSSpjuEiSqjNcWvDwmCSVMVwKeDhMktoxXCRJ1RkuBTwcJkntGC4teHhMksoYLpKk6gwXSVJ1hksLzr1IUhnDpYBzLZLUjuEiSarOcCng4TBJasdwacHDY5JUxnCRJFVnuEiSqjNcWnDuRZLKGC4FnGuRpHYMlwLusUhSO4ZLC+7BSFIZw0WSVJ3hIkmqznCRJFVnuLTgxL4klTFcCjiRL0ntGC4F3GORpHYMlxbcg5GkMoaLJKk6w0WSVJ3hIkmqznBpwYl9SSpjuEiSqjNcWvBsMUkqY7hIkqozXCRJ1XUSLhGxX0R8JSLuaH7uO812V0TE1oi4ZYf1F0XEjyJiY3M7ZddULkkq0dWeywXAdZm5CriuWZ7Kp4GTpnnsbzPzqOZ27QBqlCTNUlfhchqwtrm/Fjh9qo0y8xvAT3ZRTTPyVGRJKhPZwTtmRPw0M/fpW344M6c7NLYCuCYzj+xbdxHwR8CjwA3A+Zn58DTjzwXOBRgdHT1m3bp1rev97v3b+NjGJxkbHeG8o/doPX4xmJycZK+99uq6jE7Y+3D2DsPdf3/vxx133IbMHGszfslAqgIi4qvAgVM8dGGFp78cuATI5ueHgbOn2jAz1wBrAMbGxnJ8fLz1iz32vS2w8UYOOGAZ4+PHzLbmBW1iYoLZ/O4WA3sf77qMzgxz/3PtfWDhkpknTPdYRDwQEcszc0tELAe2tnzuB/qe6xPANbOvVJJUW1dzLuuB1c391cDVbQY3gbTdG4BbpttWkrTrdRUulwInRsQdwInNMhFxUET84syviLgS+BZwWERsjohzmoc+FBE3R8T3gOOA9+za8iVJOzOww2I7k5kPAcdPsf4+4JS+5bOmGf/WwVUnSZorP6HfgqciS1IZw0WSVJ3h0oJXRZakMoaLJKk6w0WSVJ3hIkmqznCRJFVnuEiSqjNcJEnVGS6SpOoMF0lSdYaLJKk6w0WSVJ3h0oIXrpSkMoaLJKk6w6UFL1wpSWUMF0lSdYaLJKk6w0WSVJ3hIkmqznApMLJbbyZ/6Yi/Lkkq4btlgRMPH+XUlbvzwdcd0XUpkrQgGC4FRnYLzjhsKfvuubTrUiRpQTBcJEnVGS6SpOoMF0lSdYaLJKk6w0WSVJ3hIkmqznCRJFVnuEiSqoscoq9XjIgHgXtmOXx/4McVy1lohrl/ex9ew9x/f+8vycxlbQYPVbjMRUTckJljXdfRlWHu396Hs3cY7v7n2ruHxSRJ1RkukqTqDJdya7ouoGPD3L+9D69h7n9OvTvnIkmqzj0XSVJ1hoskqTrDpUBEnBQRt0fEpoi4oOt6BikiDo2Ir0XEbRFxa0S8q1m/X0R8JSLuaH7u23WtgxIRIxFxU0Rc0ywPU+/7RMTnIuIHzd/Aq4al/4h4T/M3f0tEXBkReyzm3iPiiojYGhG39K2btt+IeH/zHnh7RPz+TM9vuMwgIkaAy4CTgcOBsyLi8G6rGqhtwPmZ+VLgWOAdTb8XANdl5irgumZ5sXoXcFvf8jD1/vfAlzLzN4Hfovd7WPT9R8TBwDuBscw8EhgBzmRx9/5p4KQd1k3Zb/MecCZwRDPmY81747QMl5m9AtiUmXdm5lPAOuC0jmsamMzckpk3Nvd/Ru/N5WB6Pa9tNlsLnN5JgQMWEYcApwKf7Fs9LL3vDfwu8CmAzHwqM3/KkPQPLAGeHxFLgBcA97GIe8/MbwA/2WH1dP2eBqzLzCcz8y5gE733xmkZLjM7GLi3b3lzs27Ri4gVwNHAd4DRzNwCvQACDuiwtEH6O+B9wLN964al918HHgT+qTks+MmI2JMh6D8zfwT8DfBDYAvwSGb+J0PQ+w6m67f1+6DhMrOYYt2iP387IvYCPg+8OzMf7bqeXSEiXgtszcwNXdfSkSXAy4HLM/No4DEW12GgaTVzC6cBK4GDgD0j4i3dVjWvtH4fNFxmthk4tG/5EHq7y4tWROxOL1g+m5lfaFY/EBHLm8eXA1u7qm+Afht4fUTcTe/w52si4p8Zjt6h97e+OTO/0yx/jl7YDEP/JwB3ZeaDmfk08AXg1QxH7/2m67f1+6DhMrPrgVURsTIiltKb1FrfcU0DExFB75j7bZn5kb6H1gOrm/urgat3dW2Dlpnvz8xDMnMFvf/O/5WZb2EIegfIzPuBeyPisGbV8cD3GY7+fwgcGxEvaP4NHE9vvnEYeu83Xb/rgTMj4nkRsRJYBXx3Z0/kJ/QLRMQp9I7FjwBXZOZfdVvR4ETE7wD/DdzMc/MOH6A373IV8GJ6/xDPyMwdJwMXjYgYB96bma+NiBcxJL1HxFH0TmZYCtwJvJ3e/4Qu+v4j4mLgTfTOmLwJ+GNgLxZp7xFxJTBO79L6DwAfBP6dafqNiAuBs+n9ft6dmV/c6fMbLpKk2jwsJkmqznCRJFVnuEiSqjNcJEnVGS6SpOoMF0lSdYaLNI9Ej/8uteD5Ryx1LCJWNN+d8jHgRn75MhvSguSHKKWONVefvhN4dWZ+u+NypCrcc5Hmh3sMFi0mhos0PzzWdQFSTYaLJKk6w0WSVJ0T+pKk6txzkSRVZ7hIkqozXCRJ1RkukqTqDBdJUnWGiySpOsNFklTd/wODZhbTm1EtBwAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "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(-0.1,0.1)\n",
+    "plt.ylabel('Ham value')\n",
+    "plt.grid()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "d97150a0",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|████████████████████████████████████████▉| 999/1000 [01:35<00:00, 10.41‰/s]\n"
+     ]
+    }
+   ],
+   "source": [
+    "# for control of time integrator and spatial grid\n",
+    "T = 48.0 # Maximum evolution time\n",
+    "num_points_t = 128 # 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": 5,
+   "id": "45fcb464",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyNUlEQVR4nO3deXxU5bnA8d+TBRJIgiRAQIKCgiJVBAGlaitaV7RivXDdd6TYuqC1hXqrdV+wFXetVau3oqhVKygVW0zcrlbAlUWUyhZ2QlgChpDkuX+8MzCELJPJOXNmMs/Xz/nMzJkz5zzvIPPwrkdUFWOMMaa50oIOwBhjTHKyBGKMMSYmlkCMMcbExBKIMcaYmFgCMcYYE5OMoAPwQ6dOnbRnz54xfXbr1q20b9/e24CSiJXfyp+q5U/lsgPMmTNnvap2bs5nWmUC6dmzJ7Nnz47psyUlJQwbNszbgJKIld/Kn6rlT+WyA4jI0uZ+xpqwjDHGxMQSiDHGmJhYAjHGGBOTVtkHYowx0dixYwelpaVUVlbSoUMHFixYEHRIvsvKyqKoqIjMzMwWn8sSiDEmZZWWlpKbm0vPnj2pqKggNzc36JB8paqUlZVRWlpKr169Wny+QJuwRORpEVkrInMbeF9E5EERWSQiX4rIYfGO0RjTelVWVlJQUICIBB1KXIgIBQUFVFZWenK+oPtAngFObuT9U4A+oW0M8FgcYjLGpJBUSR5hXpY30CYsVX1PRHo2csgI4H/VrTn/sYjsJSLdVHWVH/Hc9u5tbF61mfYr2nNwl4PJzsz24zLGGNMqJHofSHdgecTr0tC+PRKIiIzB1VIoLCykpKSkWReqrq1m4kcTqaiu4A/f/IE00ihqV0Tv9r05uMPBDNxrIPu227fV/2uloqKi2d9da2LlT63yd+jQgS1btgBQU1Oz83k8devWjVWr3E/ajBkzGD9+PNOmTaNHjx6+XbOystKTP+dETyD1/VrXewcsVX0CeAJg8ODBGsuM0k3HbuKFt14gu2c2n6/+nC/WfMGclXN4Z907AAzeezB//dlf6dupb7PPnSxSfTaulT+1yr9gwYKdHedbtmwJrBM9NzeXmTNnMn78eN5++232339/X6+XlZXFwIEDW3yeRE8gpUBkGi4CVvp1sTRJo3t2d4YdNIwzDzoTcKMWFm9czIxFM7ix+EYO+9Nh3HfSfYwdPNavMIwxKeb999/n8ssvZ/r06b4nDy8legKZClwpIlOAI4BNfvV/NERE2K/jflwx5ArO6HsGl7x+CVe8eQVHdD+Cgd1ansGNMYlh/Pi2zJ/v7TkHDID772/8mO3btzNixAhKSkro2ze5WjeCHsb7AvARcKCIlIrIZSIyVkTC/7yfDnwHLAL+DPwioFAB6JbbjclnTiZd0nl5/stBhmKMaSUyMzM58sgjeeqpp4IOpdmCHoV1ThPvK/DLOIUTlYJ2BRzX6zhenv8ydxx3R6vvVDcmVdxzz3Zyc9vE/bppaWm89NJLHH/88dx5553ccMMNcY8hVkHPA0lKI/uNZNGGRXy19qugQzHGtALt2rXjjTfeYPLkyUlVE7EEEoMz+p5BmqTx8jxrxjLGeCM/P5+33nqL22+/nddffz3ocKJiCSQGXdp34Zh9j+Hl+S/jWtmMMSY2FRUVO5/36NGDxYsXM2LEiAAjip4lkBiN7DeShWULmb/O42EbxhiTJCyBxOjMg85EEF5fmBxVTWOM8ZolkBh1zelKt9xuLNqwKOhQjDEmEJZAWqBHXg+Wb17e9IHGGNMKWQJpgR4delC6uTToMIwxJhCWQFqgR14Plm9abiOxjDEpyRJICxTlFbF1x1Y2Vm4MOhRjTJLKyckJOoSYWQJpgR55bqFg6wcxxqQiSyAt0KODSyDWD2KMaamKigp+8pOfcNhhh3HIIYfsnI2+ZMkS+vbty0UXXUT//v0ZOXIk27ZtA2DChAn069eP/v37c/311wOwbt06/uu//oshQ4YwZMgQPvzwQ99iTvTl3BNaUV4RAMs3WQ3EmGQ3vng88zd4OzF4QNcB3H/y/VEdm5WVxWuvvUZeXh7r169n6NChnH766QAsXLiQp556iqOOOopLL72URx99lEsvvZTXXnuNr7/+GhFh48aNAFxzzTVce+21HH300SxbtoyTTjqJBQsWeFquMEsgLdAtpxvpkm5NWMaYFlNVbrjhBt577z3S0tJYsWIFa9asAdwSJ0cddRQA559/Pg8++CDjxo0jKyuL0aNHc+qpp3LaaacB8K9//Yv5ETc22bx5s293W7QE0gLpaensnbu3JRBjWoF7jr0nsFvaAkyePJl169YxZ84cMjMz6dmzJ5WVlQB73DZCRMjIyOCTTz5h5syZTJkyhYcffph33nmH2tpaPvroI7Kzs32P2fpAWsjmghhjvLBp0ya6dOlCZmYmxcXFLF26dOd7y5Yt46OPPgLghRde4Oijj6aiooJNmzYxfPhw7r//fj7//HMATjzxRB5++OGdnw3v94MlkBYqyiuyPhBjTIudd955zJ49m8GDBzN58uTdbm970EEH8eyzz9K/f382bNjAFVdcwZYtWzjttNPo378/xxxzDJMmTQLgwQcfZPbs2fTv359+/frx+OOP+xazNWG1UI+8HkxdOBVVtbsTGmOaLbyce6dOnXbWMiItWbKEtLS0PRJBu3bt+OSTT/Y4vlOnTrz44ov+BFuH1UBaqEdeDyqrKyn7vizoUIwxJq4sgbSQzQUxxvipZ8+ezJ07N+gw6mUJpIVsLogxyS3V1rLzsryBJhAROVlEForIIhGZUM/7HURkmoh8ISLzROSSIOJsjC1nYkzyysrKoqysLGWSiKpSVlZGVlaWJ+cLrBNdRNKBR4ATgFJglohMVdXIqaC/BOar6k9FpDOwUEQmq2pVACHXqzCnkMy0TKuBGJOEioqKKC0tZd26dVRWVnr2w5rIsrKyKCoq8uRcQY7COhxYpKrfAYjIFGAEEJlAFMgVN7wpB9gAVMc70MakSRrd87pTuqX+PpC334ZzzoEBA2DYMPjlLyE/P64hGmMakJmZSa9evQAoKSlh4MCBAUeUXIJMIN2ByH+2lwJH1DnmYWAqsBLIBc5S1dr6TiYiY4AxAIWFhZSUlMQUVEVFRbM/m1uby1dLv6r3c5MmHci2bV1YuvR7fv/79vzrX+u55ZZ5McUWD7GUvzWx8qdu+VO57DFT1UA2YBTwZMTrC4CH6hwzEpgECNAbWAzkNXXuQYMGaayKi4ub/Zlz/naO7vfAfnvsr61V7d5ddeRI9/rGG1VB9YsvYg7Pd7GUvzWx8hcHHUJgUrnsqqrAbG3m73iQneilQI+I10W4mkakS4BXQ+VbhEsgfUkw3XK6saZizR77v/4aVqyAE05wr6+9FvLy4NZb4xygMcb4IMgEMgvoIyK9RKQNcDauuSrSMuAnACJSCBwIfBfXKKOQn53P1h1b2V69fbf9b7/tHsMJpGNHGDcOXnkFvvwyvjEaY4zXAksgqloNXAnMABYAL6nqPBEZKyJjQ4fdBhwpIl8BM4Hxqro+mIgbVtCuAIDyyvLd9v/zn9C7N4T66ACXQPLy4JZb4higMcb4INC1sFR1OjC9zr7HI56vBE6Md1zNlZ/thlWVbSuja05XAKqqoKQELrxw92M7doQrr4S77oJly2CffeIcrDHGeMRmonsgnEA2fL9h576PP4atW+HEetLf6NGgCs88E6cAjTHGB5ZAPFBfAnn7bUhPh2OP3fP4Xr3g+OPhqaegpiZeURpjjLcsgXigINv1gUQmkHffhcGDoUOH+j8zerRrwpo5Mx4RGmOM9yyBeGBnH0jEku7LlsEBBzT8mTPOcDPSn3rK5+CMMcYnlkA8kNMmh4y0jJ01EFVYtQr23rvhz7RtCxdcAK+9BusTblyZMcY0zRKIB0SE/Oz8nQmkrAx27IBu3Rr/3GWXueP++tc4BGmMMR6zBOKRguyCnQlk1Sq3r6kEcsghcMQR8OSTrtZijDHJxBKIR/Kz83f2gawMLcjSWBNW2OjRMH8+/PvfPgZnjDE+sATikcgmrGhrIABnnQXt27taiDHGJBNLIB6JNYHk5rokMmUKbNniY4DGGOMxSyAeiewDWbnSzf9o1y66z44e7Watv/iijwEaY4zHLIF4JD87n4qqCqpqqli1KrraR9jQodCvn80JMcYkF0sgHolczqS5CUTE1UI+/hjmzvUpQGOM8ZglEI9EJpCVK6MbgRXpggsgM9NqIcaY5GEJxCPhe4KUbWt+DQSgUye3vMn//i9s397k4cYYEzhLIB4J10CWrStj+/bmJxBwzVgbNsDrr3scnDHG+MASiEfCCWTxGjcSq7lNWOCWeN93X5sTYoxJDpZAPBJe0n35epdAYqmBpKXBpZe6W+EuWeJhcMYY4wNLIB4Jr8i7aqNbziSWBAJw8cVuVNZf/uJdbMYY4wdLIB4Jr8i7dkvsNRBw90g/6SR4+mm7W6ExJrFZAvFQeDmTnBy3REmsRo+G0lJ3W1xjjElUgSYQETlZRBaKyCIRmdDAMcNE5HMRmSci78Y7xuYoyC5gU9WGmGsfYT/9KXTubJ3pxpjEFlgCEZF04BHgFKAfcI6I9KtzzF7Ao8DpqvoDYFS842yO/Ox8tmpZTCOwIrVpAxdeCFOnwtq13sRmjDFeazSBiMixIvKwiLwuIq+IyN0i0tujax8OLFLV71S1CpgCjKhzzLnAq6q6DEBVE/rnND87n+3S8hoIuLsVVldbZ7oxJnFlNPSGiNwNFAIzga7AYuA/wMsicqeqvtzCa3cHlke8LgWOqHPMAUCmiJQAucADqvq/DcQ7BhgDUFhYSElJSUxBVVRUxPzZreu3Up25gerq5ZSU/Cemc0QaOPBQ/vjHbAYN+jcZGfG5ZWFLyt8aWPlTt/ypXPaYqWq9G/BVxPMM4MPQ847A3IY+F+2Ga456MuL1BcBDdY55GPgYaA90Ar4FDmjq3IMGDdJYFRcXx/zZ3719u3IzetfE7TGfI9Kbb6qC6nPPeXK6qLSk/K2Blb846BACk8plV1UFZmszf8cba8KqFZH80PO9gfRQwikHxIPcVQr0iHhdBKys55i3VHWrqq4H3gMO9eDavkjb7r6u3C4bPDnfySfDQQfBH/9o90w3xiSexhLIncBnIvI28AFwG4CIdAa+8ODas4A+ItJLRNoAZwNT6xzzOvAjEckQkXa4Jq4FHlzbF+mhBJKZ500CSUuDa6+Fzz4Dq1kbYxJNgwlEVV8EBgI3AP1V9c3Q/nWqem5LL6yq1cCVwAxcUnhJVeeJyFgRGRs6ZgHwFvAl8AmuySth75ghoQRCtjcJBNwy7507wx/+4NkpjTHGEw12ooeU45qujhMRxTUxfRJqL2sxVZ0OTK+z7/E6r+8F7vXier7b5tbDqm3rXQLJyoKrroKbboLPP4cBAzw7tTHGtEiDNRARORHXaX0zMBw4FbgF+Db0nqmjpsLVQKoyyjw971VXQV4e3Hmnp6c1xpgWaawG8gBwvKouidwpIr1wtYaDfIwrKVVtcgmkUryrgQDstRdceSXcdRcsWOA61o0xJmiNdaJn4EZB1bUCyPQnnOS2bWMu1GSwabu3CQRcZ3p2ttVCjDGJo7EE8jQwS0TGi8i5oW088G/A7txdj43lQlqVW1DRa506wRVXwPPPw6JFnp/eGGOarbFRWHfhlhIR4IfAkaHn54XeM3WUl0PmjnzKvve2DyTs+uuhbVu45RZfTm+MMc3S6Cis0DDahJ13kWg2boS26k8NBKBrV7j6apg4EcaPh4MP9uUyxhgTlZhW4xWRmz2Oo1UoL4ds/EsgAL/5jbvXyI03+nYJY4yJSqzLuc/xNIpWorwcctILfE0g+fnw61/D3/8On3zi22WMMaZJMSUQVZ3mdSCtQXk5dMj0rw8kbNw4Nzt9wgRbI8sYE5ym7gdykog8JiJTQ/cEeUxETo5XcMlkxw7YuhX2ysqnoqqCqpoq366VkwM33wzFxfDGG75dxhhjGtXYTPT7gWuAd4GJuOVE3gWuFpEH4hJdEtm40T0WtHOTCcu/L/f1emPGQN++bmTWjh2+XsoYY+rVWA1kuKoOV9UpqvpBaJuCW9JkeJziSxrloXzRJceth+VnPwhARgbcey988w386U++XsoYY+rVWAKpFJHD69k/BKj0KZ6kFU4gXfNcDcTvfhCAU0+F445zzVkb/M1Xxhizh8YSyMXAQyIyX0TeDm0LgIdC75kI4QSyd75LIH7XQABEYNIk13z2P//j++WMMWY3DU4kVNVPgSNEpCvu/uUClKrq6ngFl0zCfSA9CuKXQAD693eTC++/Hy65BA6vr85ojDE+aHIYr6quVtU5qjrbkkfDwjWQnoXx6QOJdPPN0K2bWyurpiZulzXGpLhYJxKaOsIJpKhzLumSTtk2//tAwvLyXFPWp5/C4483fbwxxnjBEohHysvd3QOzs4X8bH+XM6nPqFFwwgmuL2S11RONMXHQZAIRkf1FpG3o+TARuVpE9vI9siRTXg4dO7rn+dn5bKiMbwIRgYcfhu+/d0udGGOM36KpgbwC1IhIb9x9QHoBz/saVRLauHFXAilo5+96WA054AC3Su9zz0FJSdwvb4xJMdEkkFpVrQZ+BtyvqtcC3by4uIicLCILRWSRiExo5LghIlIjIiO9uK4f6tZA4tkHEum3v4VevWDsWKi02TrGGB9Fk0B2iMg5wEVAeOWlFt/SVkTSgUeAU4B+wDki0q+B4+4BZrT0mn4qL3f3LgcC6QMJy852HekLF7rRWcYY45doEsgluDsS3qGqi0WkF/CcB9c+HFikqt+pahUwBRhRz3FX4ZrR1npwTd/sVgPJCi6BAJx4Iowe7ZY6sSXfjTF+iWYeyHxVvVpVXwi9Xqyqd3tw7e7A8ojXpaF9O4lId1zTWcIPTo1MIAXtCthStYUdNcGtcviHP8Dee7vJhdaUZYzxQ4Mz0UXkJVX9bxH5CtjjrhOq2r+F15Z69tW9zv3AeFWtEanv8IiTiYwBxgAUFhZSEmMvckVFRbM/W1sLmzcfw6ZNSykpWcL6FesBmDZzGvlt8mOKwwtXXdWR8eMP5bLLlnL55Yuj+kws5W9NrPypW/5ULnvMVLXeDegWety3vq2hz0W74ZrFZkS8/i3w2zrHLAaWhLYKXDPWGU2de9CgQRqr4uLiZn9mwwZVUJ00yb1+4asXlJvR+WvnxxyHVy67TDUtTfWTT6I7PpbytyZW/uKgQwhMKpddVRWYrc38HW9sLaxVocelnmctZxbQJ9SnsgI4Gzi3Tgy9ws9F5BngDVX9u0/xxCw8Cz3ciV6QHf/lTBryxz/CjBmuKWvOHGjbNuiIjDGtRWAz0dUNDb4SN7pqAfCSqs4TkbEiMjaouGIRTiCRw3ghMRJIhw7wxBMwbx7cckvQ0RhjWpMGayDxoKrTgel19tXbYa6qF8cjplg0lEDicU+QaJxyClx2GdxzDwwfDkcfHXRExpjWoKl7oqeLiBdDdlu18FLuiVgDCbv/fjfB8PzzYdOmoKMxxrQGjSYQVa0BOotImzjFk5Tq1kDy2uaRLukJlUByctwSJ6WlcNVVQUdjjGkNomnCWgJ8KCJTga3hnap6n19BJZu6negiwazI25ShQ+HGG90M9VNPhbPOCjoiY0wyi6YTfSVuCZM0IDdiMyHl5ZCRAe3b79qXn52fMH0gkf7nf1wiGTsWli9v+nhjjGlIkzUQVbWxO00Iz0KPnOvYqV0n1m1dF1xQDcjIcE1ZAwbAhRfCzJmQZneFMcbEwH46PFBWBp067b6ve153VmxZEUxATdh/f3jwQbfk+x//GHQ0xphkZQnEA+vXQ0HB7vuKcoso3VwanlGfcC6+GM480zVpffZZ0NEYY5JRNHckPCqafamsvhpIUV4R23ZsY2PlxkBiaoqIm2DYuTOcd567k6ExxjRHNDWQh6Lcl7LqrYHkFQFQurk0gIiiU1AAzzwDCxbAb34TdDTGmGTT2Gq8PwSOxM0DuS7irTwg3e/AkoVqw30gACu2rOCQwkMCiCw6J5wA48a5iYbDh7sbUhljTDQaq4G0AXJwSSZy+O5mIGFvLRtvFRVQVZWcNZCwu+6CQw5xCy6Wl7f4ZpPGmBTR2Gq87wLvisgzPq7Im/TKQlM96tZAuuV0Q5CkSCBZWTB5MgwZAvfeeyBnnLH7kGRjjKlPNH0g20TkXhGZLiLvhDffI0sS6929o/aogWSmZ1KYU5gUCQRcDeSee+Cjjzrxpz8FHY0xJhlEk0AmA18DvYBbcEubzPIxpqTSUA0EXDNWos4Fqc9VV8GQIRu47jr4+uugozHGJLpoEkiBqj4F7FDVd1X1UmCoz3EljYZqIOASSLLUQMDNSB8//mvat3dDe6uqgo7IGJPIokkgO0KPq0TkVBEZCBT5GFNSaawG0j23e1IlEICCgiqefBI+/RRuuinoaIwxiSyaBHK7iHQAfgVcDzwJXOtrVElk/XrX4RxeiTdSUV4RGys3srVq655vJrARI2DMGJg40S13Yowx9WkygajqG6q6SVXnquqxqjpIVafGI7hkUFYG+fmQXs/MmPBQ3mTqBwm77z7o0wcuuGDXcvXGGBPJ1sJqofpmoYcl01yQutq3d0N7V692S78n6JJexpgAWQJpofpmoYd1z3Wz0ZMxgQAMHgy33AIvveQ2Y4yJZAmkhRqrgexczmRz8jVhhf3mN3DEEfCLX7jaiDHGhEWzGm+hiDwlIv8Ive4nIpd5cXEROVlEForIIhGZUM/754nIl6Ht/0TkUC+u66XGaiDtMtuRn52ftDUQcDegeuYZ2LrVmrKMMbuLpgbyDDAD2Dv0+htgXEsvLCLpwCPAKUA/4BwR6VfnsMXAMaraH7gNeKKl1/VaWVnDNRAIDeXdkrwJBKBvX7jjDnj9dXj++aCjMcYkimgSSCdVfQmoBVDVaqDGg2sfDixS1e9UtQqYAoyIPEBV/09Vw2OAPibB5p9s2+buo9FYAinKK0rqJqywcePgyCPdbPWVK4OOxhiTCJq8JzqwVUQKAAUQkaHAJg+u3R1YHvG6FDiikeMvA/7R0JsiMgYYA1BYWEhJjBMYKioqov7s2rVtgR9SVraQkpJV9R6TtjWN79Z/F3M88dZY+ceOzWb06MGMHFnOHXfMbZULLjbnz781SuXyp3LZY6aqjW7AYcCHuKTxIa4Jq39Tn4vivKOAJyNeXwA81MCxxwILcMuqNHnuQYMGaayKi4ujPvbTT1VB9dVXGz7m5uKblZvR7dXbY44pnpoq/6RJrszPPBOXcOKuOX/+rVEqlz+Vy66qCszWZv6ORzOR8FPgGNzNpX4O/EBVv/Qgd5UCPSJeFwF7NI6ISH/c7PcRqlrmwXU909gyJmH7dNgHgKUbW8eK+FdfDT/6EVxzjY3KMibVRTMK60LgXGAQrjZyTmhfS80C+ohILxFpA5wN7DbDXUT2AV4FLlDVbzy4pqcaW0gxrH9hfwA+X/25/wHFQVoa/PnPru9n3LigozHGBCmaTvQhEduPgJuB01t6YXWd8VfiRngtAF5S1XkiMlZExoYOuwkoAB4Vkc9FZHZLr+ulaGogB3c5mIy0jFaTQAAOPBB+9zt48UV4882gozHGBKXJTnRVvSrydWhhxb96cXFVnQ5Mr7Pv8Yjno4HRXlzLD+EaSH5+w8e0zWhLv879+Gz1Z/EJKk7Gj4cpU9wEw3nzICcn6IiMMfEWy0z0bUAfrwNJRmVlbhXejCbS8ICuA1pdAmnTBp54ApYtg9tuCzoaY0wQoukDmSYiU0PbG8BC4HX/Q0t8jS1jEmlg14GsrljN6orW1et81FFwySVu5d4FC4KOxhgTb9HMA/lDxPNqYKmqJvfUao80toxJpIFdBwLw2arPOKXPKT5HFV933w2vveYmGP7zn7TKuSHGmPpFM4z33YjtQ0seu0RbAxnQdQBAq2vGAujSBW6/HWbOhL/9LehojDHx1GACEZEtIrK5nm2LiGyOZ5CJKtoaSIesDuzXcb9WNRIr0tixMGAAXH+9G95rjEkNDSYQVc1V1bx6tlxVzYtnkIlINfoaCLhmrNZYAwF3N8b77nMd6g88EHQ0xph4iXoUloh0EZF9wpufQSWDNWvcEuc9e0Z3/MCuA1m0YRGbt7fOytuxx8Lpp8Odd7rvxhjT+kUzCut0EfkWt7T6u8ASGlnUMFXMm+cef/CD6I4P94N8sfoLfwJKABMnuias3/8+6EiMMfEQTQ3kNmAo8I2q9gJ+gltUMaU1N4EM7OZGYn266lOfIgregQfCFVe4pU6+/jroaIwxfosmgewILWKYJiJpqloMDPA3rMQ3fz507AiFhdEd3y2nG73zezP1m6lNH5zEfvc7yM6Gm28OOhJjjN+iSSAbRSQHeA+YLCIP4OaDpLR581ztI9p5DyLCuQefS/Hi4lZxg6mGdOniFll88UX4ovW21hljiC6BjMAtX3It8BbwH+CnfgaV6FR3JZDmOK//eSjKlLlT/AksQfzqV9ChA9x4Y9CRGGP8FE0CGQPsrarVqvqsqj6YaPfliLc1a6C8HPrVvYN7Ew4oOIAhew/hua+e8yewBNGxI/z61zBtGvz730FHY4zxSzQJJA+YISLvi8gvRSTKVv/Wq7kd6JHO738+n6/+nHlr53kbVIK55ho3R+aOO4KOxBjjl2iWMrlFVX8A/BLYG3hXRP7le2QJLJxAmlsDATjrB2eRLulM/mqyt0ElmJwcd/fCadPgq6+CjsYY44fmLOe+FlgNlAFd/AknOYRHYHXt2vzPFuYUcvx+x/PsF8+yZfsW74NLIFde6RLJ3XcHHYkxxg/RTCS8QkRKgJlAJ+ByVe3vd2CJrLkjsOq66ZibWLVlFeP/Nd7bwBJMfj78/OfuxlPffRd0NMYYr0VTA9kXGKeqP1DV36vqfL+DSmThEVixNF+FHdnjSK4dei2PzX6Mdxa/411wCei669wNt+69N+hIjDFei6YPZIKqfh6HWJJCeARWLB3okW4/7nYOKDiAS1+/lI2VGz2JLRHtvTdceCE888yue8gbY1qHWG5pm9JaMgIrUnZmNs+MeIYVW1Yw9MmhfFP2TcuDS1BXXw2VlfDUU0FHYozxkiWQZqipcavNZme7+1+01A97/JCZF86k7PsyDv/z4UyZO4Wa2pqWnzjBHHIIDBsGjz7qvkNjTOsQaAIRkZNFZKGILBKRCfW8LyLyYOj9L0XksCDiDLvnHnjnHXjooejvA9KUH+/7Y2ZdPoteHXtxzivn0PeRvjzyySMs3bjUmwskiKuugqVL3bBeY0zrEM090X0hIunAI8AJQCkwS0Sm1umkPwXoE9qOAB4LPcZVTY273/dNN8FZZ8Gll3p7/p579WT25bN5dcGrTPy/iVz5jyu58h9X0ju/N4P3Hkzfgr7sn78/ndt1plO7TnRu7x7bZbbzNhAfnX469OgBDz4IZ5wRdDQm2akqiu58Dux8Heu+yppKtu3Y5tn5mrvPC2mSRqd2Udwm1SMSLki8icgPgZtV9aTQ698CqOpdEcf8CShR1RdCrxcCw1R1VWPnHjx4sM6ePbvZMY0aBcuXl9GhQwHV1VBdDdu2uXkf27ZBr17w2WdunSe/qCrz1s1j5nczeWfJO3y15iuWbFxS7/9oaZJG2/S2tM1ou9tjVkYWmemZCIKIkCZpCKFHkQafp0kaG8s3kp+fv/MzXvrPd/DNQvjRj9z8kES0YcMG8vPzfb+Oojt/BBt6BJo8JvIx2vM2dv5t27aRlZ3lyfm9iqm+Y0z9CtsXsvr61TF9VkTmqOrg5nwmsBoI0B1YHvG6lD1rF/Ud0x3YI4GIyBjcul0UFhZSUlLS7IAWLBjI9u3pbNmyifR0JSPDbaecso3evbcwdOgGPvtsR7PPG4tDOZRDux4KXWF7zXbWbl/Lph2b2LhjI5t2bGLTjk18X/M9O2p3UKVV7rHWPe6o2UF1dfVufyFrQv/Vau3O/fU9r66pZuvarb78Ra1tK9Auh7nfVVFYuN3z83uhpqaGTWs2+XoNRXcmZwn/J/W/juaYlnwm9GLnPyZqsmrIyMjYdTzi3o84vu6+pl43+Bl2XTfydd346jtmt/jrfC7yO6irseOqqqpo06ZNi89X3/HRHtdSbdPbxvTbF6sgE0h931rdX61ojnE7VZ8AngBXAxk2bFizA5o7F0pKStjzs/7/izRR1F9+74wcCe++Bl+sgNDf1YTid/kTXSqX34+yq7qWjJoaqK3dfatvX31bNMeFj8nIgB/+0NMiNCrIBFIK9Ih4XQSsjOEYk0QuvRReeQXeeAPOPDPoaOKvthYWLXL3Spk71w0sKC11c2QqKtwtgUUgLQ3y8qBTJzdgo1Mn6NwZ9tsP+vRxW5cusa+GYHaprnZzu5Yvz+bjj2HDBreVl7vHTZtg69Y9t23b3J9XVRXs2FH/Fu9Rh4WFsDq2FqyYBJlAZgF9RKQXsAI4Gzi3zjFTgStFZAqueWtTU/0fJrGddBJ07+7mhKRKAtmwAV57DWbMcKP4whMq09LcRMuiIrfl5roh4uB+eDZtcsfOnw/r17vntbW7zpuX524jPHAgHHaY2w45BLKy4l/GRFVVBYsXu6S9YgWsXOm2Vat2PV+zxtUUGhqfk5MD7du7rV27Xc8LCtyfV2Zm01tGhvvzrrulp9e/P9bj4v1nH1gCUdVqEbkSmAGkA0+r6jwRGRt6/3FgOjAcWIS7qdUlQcVrvJGeDhdfDHfd5f5Cd+8edET+UHUJ47HH4B//cP8a7d4dTjsNfvxjN4+oX7/m/YWvrnY1lm++gW+/dY/z58NLL8ETT7hjMjLcJNdwQhk82F2rNSeVyCTx7bduCz9funT3pCviam7durnkfdhh7nnnzrB69XyOProf+fns3Dp0cN+pqV+gX42qTsclich9j0c8V9wy8qYVueQSd5+QZ5+FG24IOhpv1da6xSPvuss1UXXt6ubAnHuu+7FqSZNTRgbsv7/bTjll135VWLIE5syBTz9127Rp8Je/7PrcIYfAkCEuoQwZ4pJMZmaLihpXO3Y0niQim4o6dHBNfEOHwvnnu+e9e7th5IWFDZe7pGQtw4a1YJG7FGS51cTd/vu7f4U/9xz89retpx3//ffd4pGzZ7sf7GefhbPP9n+wgIgbYt6rlxukAC6plJa6WGbNcltkTSUry9VMhgxxie3gg+Ggg/yNsynbt7tk8J//7EoS4USxZMnuSSIvzyWGIUNccg73C/Xu7fqLWsv/U4nOEogJxLnnwtix8OWXcOihQUfTMhUVcP318Kc/ub6Mv/7VlS8twHUeRNy/uHv0gJ/9zO1TdT/Os2btSixPP+1WVgjr1u0IBg1yP8Y9e8K+++567NAh9h/m2lrXKb1mjduWLXM1iu++c4+LF7v+iMhpaTk5Lo5Bg1wiDieIPn1ck5MlieBZAjGBGDnS3XDq+eeTO4F8/LFLFkuWwK9+Bbfe6jpaE5GI+wHu3RvOOcftq6lxSWXePLe9884WlizJZuZMN8IoUmamu5Fafr573Gsvty/cwZue7hLAtm27j1YqL4d161wfTt14und3Nafjj99ViwqPNLNRZonPEogJREGBG5H1wguuvyDIf63H6umn4Yor3I/gu++6GfbJJj0dDjjAbT/7GRx99HyGDeuCqhv5tWSJa1ZautS9Dg9t3bAB1q51SSE8DyHcxBQ5SmmffVzNpbBw19ali6sZ7bMPtG0baPFNC1kCMYE591x480348MPk+vGtrXV9HQ88ACec4DrN47D6SVyJuGaizp1dP4Mx9UnCf/eZ1uL0011zzwsvBB1J9Kqr4aKLXPIYNw6mT299ycOYaFkCMYHJyYERI9zooB3xWWKsRaqq3GrMzz3nhiFPmmRzBExqswRiAnXWWW6GdXFx0JE0rrYWLrgAXn0V7r+/9c1fMSYWlkBMoE46ydVE/va3oCNpmKprrnrpJZg4Ea65JuiIjEkMlkBMoLKy4Kc/dWtF1R3mmSgmTnRzJa67zs33MMY4lkBM4EaNckNE43gbg6hNn+5my599Ntx7r81LMCaSJRATuJNPdvMGXn456Eh2t2iRG2p86KFu9eBknKtijJ/sr4QJXHa2W6U2kZqxtm1zE+vS011ciTq73JggWQIxCWHUKLfcxXvvBR2J86tfuaU9pkxxa0EZY/ZkCcQkhFNOcf/KT4RmrGnT4PHHXYf5CScEHY0xicsSiEkI7dq5ZqxXX43/bUAjbdjQhksvdUud33ZbcHEYkwwsgZiEMXKkW6Dv/feDub4qTJrUh4oKmDzZFvozpimWQEzCGD7cdagH1Yz1yivwwQedufVWd7tZY0zjLIGYhNG+PZx6qvshj3czVnm5uz9Jnz5buPba+F7bmGRlCcQklFGj3B3rPvwwvtf99a/dZMZf/3qhLZBoTJQsgZiEMny4W94kns1YH3/sJgpedx306VMRvwsbk+QCSSAiki8i/xSRb0OPHes5poeIFIvIAhGZJyK2hF0KyMlxSeSVV9wKuH6rrYWrr4Zu3eDGG/2/njGtSVA1kAnATFXtA8wMva6rGviVqh4EDAV+KSLWtZkCRo2CVavi04z17LMwaxbccw/k5vp/PWNak6ASyAjg2dDzZ4Ez6h6gqqtU9dPQ8y3AAqB7vAI0wTn1VDeE1u8l3jdvdgslDh0K553n77WMaY1EVeN/UZGNqrpXxOtyVd2jGSvi/Z7Ae8DBqrq5gWPGAGMACgsLB02ZMiWm2CoqKsjJyYnps61BopT/xht/wNdf5/Hiix/5tojh44/vx4sv7sNjj82hb98tQOKUPyipXP5ULjvAscceO0dVBzfrQ6rqywb8C5hbzzYC2Fjn2PJGzpMDzAHOjPbagwYN0lgVFxfH/NnWIFHKP3myKqh+8IE/51+4UDUzU/Xii3ffnyjlD0oqlz+Vy66qCszWZv7O+zZgUVWPb+g9EVkjIt1UdZWIdAPWNnBcJvAKMFlVX/UpVJOATjvNNWO9/DIcdZT357/uOjfa6667vD+3MakiqD6QqcBFoecXAa/XPUBEBHgKWKCq98UxNpMA8vLc7W79GI31j3/Am2+6UVddu3p7bmNSSVAJ5G7gBBH5Fjgh9BoR2VtEpoeOOQq4ADhORD4PbcODCdcEYdQoKC2FDz7w7pxVVXDttdCnj93b3JiWCmTOraqWAT+pZ/9KYHjo+QeA3UA0hf3sZ64m8uST8OMfe3POhx+GhQvhjTegTRtvzmlMqrKZ6CZhtW/vbin78sturaqWWrMGbrnF3UJ3uNVljWkxSyAmoV1+OVRWwvPPt/xcv/udu1XtpEkgVrc1psUsgZiEdthhbvvzn939OmI1a5Zb7+rqq6FvX+/iMyaVWQIxCe/yy+GLL2D27Ng+X1MDY8e6EVe//723sRmTyiyBmIR37rnulrcPPRTb5x99FD791DVd5eV5G5sxqcwSiEl4eXnwi1+428wuWNC8z65a5fo+TjwR/vu//YnPmFRlCcQkhfHj3aism26K/jOq8POfw/bt8Mgj1nFujNcsgZik0KmTmwD4t7+55qhoPPQQTJsGEydC797+xmdMKrIEYpLGdddBx45www1Nj8j67DN3m9rTToOrropPfMakGksgJml06ODWr5oxw9UqGrJ8uVsGpVMn+MtfrOnKGL8EspSJMbG65hr45BOYMAH23x9Gjtz9/f/8B37yEzdz/a23XBIxxvjDEohJKmlprlaxbBlccIFLGBdc4JZmf/llt1RJVRW88w4MGhR0tMa0btaEZZJOVhb8/e/uVrQTJkCPHm6S4NixkJ8P775rycOYeLAaiElKnTtDcTF88w0895wbqnv22TBggPV5GBMvlkBMUjvgALj11qCjMCY1WROWMcaYmFgCMcYYExNLIMYYY2JiCcQYY0xMLIEYY4yJiSUQY4wxMbEEYowxJiaWQIwxxsREtKl1sZOQiKwDlsb48U7Aeg/DSTZWfit/qpY/lcsOcKCq5jbnA61yJrqqdo71syIyW1UHexlPMrHyW/lTtfypXHZw5W/uZ6wJyxhjTEwsgRhjjImJJZA9PRF0AAGz8qe2VC5/KpcdYih/q+xEN8YY4z+rgRhjjImJJRBjjDExsQQSIiIni8hCEVkkIhOCjsdvIvK0iKwVkbkR+/JF5J8i8m3osWOQMfpJRHqISLGILBCReSJyTWh/SnwHIpIlIp+IyBeh8t8S2p8S5QcQkXQR+UxE3gi9TpmyA4jIEhH5SkQ+Dw/hbe53YAkE9z8S8AhwCtAPOEdE+gUble+eAU6us28CMFNV+wAzQ69bq2rgV6p6EDAU+GXozzxVvoPtwHGqeigwADhZRIaSOuUHuAZYEPE6lcoedqyqDoiY/9Ks78ASiHM4sEhVv1PVKmAKMCLgmHylqu8BG+rsHgE8G3r+LHBGPGOKJ1Vdpaqfhp5vwf2QdCdFvgN1KkIvM0ObkiLlF5Ei4FTgyYjdKVH2JjTrO7AE4nQHlke8Lg3tSzWFqroK3A8s0CXgeOJCRHoCA4F/k0LfQagJ53NgLfBPVU2l8t8P/AaojdiXKmUPU+BtEZkjImNC+5r1HbTKpUxiIPXss/HNKUBEcoBXgHGqulmkvv8VWidVrQEGiMhewGsicnDAIcWFiJwGrFXVOSIyLOBwgnSUqq4UkS7AP0Xk6+aewGogTinQI+J1EbAyoFiCtEZEugGEHtcGHI+vRCQTlzwmq+qrod0p9R0AqOpGoATXJ5YK5T8KOF1EluCaq48TkedIjbLvpKorQ49rgddwTfnN+g4sgTizgD4i0ktE2gBnA1MDjikIU4GLQs8vAl4PMBZfiatqPAUsUNX7It5Kie9ARDqHah6ISDZwPPA1KVB+Vf2tqhapak/c3/V3VPV8UqDsYSLSXkRyw8+BE4G5NPM7sJnoISIyHNcumg48rap3BBuRv0TkBWAYbgnrNcDvgb8DLwH7AMuAUapat6O9VRCRo4H3ga/Y1Q5+A64fpNV/ByLSH9dJmo77h+RLqnqriBSQAuUPCzVhXa+qp6VS2UVkP1ytA1xXxvOqekdzvwNLIMYYY2JiTVjGGGNiYgnEGGNMTCyBGGOMiYklEGOMMTGxBGKMMSYmlkCMiQMR2UtEfhF0HMZ4yRKIMfGxF2AJxLQqlkCMiY+7gf1D9164N+hgjPGCTSQ0Jg5CK/6+oaopsWChSQ1WAzHGGBMTSyDGGGNiYgnEmPjYAuQGHYQxXrIEYkwcqGoZ8KGIzLVOdNNaWCe6McaYmFgNxBhjTEwsgRhjjImJJRBjjDExsQRijDEmJpZAjDHGxMQSiDHGmJhYAjHGGBOT/weglw4TcmHpHgAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# Plot a single point versus time\n",
+    "var1 = idx_K\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": 6,
+   "id": "61ee31bd",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEGCAYAAACQO2mwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+g0lEQVR4nO3deXzcdZ348dd7jmRyn216pKV3S+ltoRUstqBcKsWTIipUkAXEKusPQdEVuu4CW0RFcJHl8IQuVEFwCyhHEBEpVw9KW3qmTdq0uZt7rvfvj+80JGmSTpKZTJt5Px+P6Vzf4/1JmnnP5/yKqmKMMcZ0x5XoAIwxxhy/LEkYY4zpkSUJY4wxPbIkYYwxpkeWJIwxxvTIk+gA+qqwsFDHjRvX7/2bmprIyMiIXUAnECt7cpYdkrv8Vnan7G+99VaVqg7r6zFOuCQxbtw43nzzzX7vX1JSwuLFi2MX0AnEyr440WEkTDKX38q+GAARKe3PMay5yRhjTI8sSRhjjOmRJQljjDE9siRhjDGmR5YkjDHG9ChuSUJEHhKRQyLybg/vi4jcLSI7RGSjiMyLVyzGGGP6J541iV8B5/Xy/vnA5MjtKuC/4xiLMcaYfojbPAlV/ZuIjOtlk6XAb9RZq/yfIpIrIiNV9UC8YjpeqCqNbUFqmvwcbgnS2BakJRCkNRDGHwzjD4UJhZVgWAmFwoTU2ScUVhQIq3JkhXeNPNb2Y0fuOXoJ+L27mtja9E9SQ414Qy14wq14wq24wwHcGsAVDuDWIKIhXBpCCCMaxhW5hzCiiqCIKuA8JnKT9lM6DyRyHyBIMwFaNYCfMAFCBMS5D6GECBMiTFiVMJGbaOSoH9yHO5RMu5RQj/oXCCtefxhvW4hwQwtvv3obrqDiDinuoOIKKa6wE7eEFVfY+QFKOBK74rynkaMqkXLToaydH7f/AvCgpABeVL0oXue5eFHcON/PXKCu9seKC0SAbm7t55Ae7jv+xI/8vwg5N8KENcy2X/3xg5+ghjv/FLXzT7PLT/JoR7YXetjy+LkEgaqy+aHfdfPG4McyUFmj8/jK7asG9ZyJnEw3GtjX4XlZ5LWjkoSIXIVT26CoqIiSkpJ+n7SxsXFA+/dVKKzsqAuzrTZE6eEwBxrDVLUq/lB8zpdGK7NkNzNcu5ggFZwkFYyUGoZJPdnSDPtjf85mEbaleNmSksI+r4cDHg8VHjeH3G7qXW78Ljn2Qfoho0WZUKEU1cGwemVYPRTWK4UNkN4K6f64nBYFWn0FNGaMpimjCH9KDv6UHNpSsvGnZOFPySbkSYvPuVVBGwmHatFwLRqqRcOHUW0FbUXDzj0E43J+k1i6L9ynz69YfN4lMkl098nRbW5X1fuB+wHmz5+vA5k9OVizLw81tPLwq3t47I19VDc5n1YTCjOYOS6T4rx0RmT7yMtIIdvnIdPnIT3Fg8/rIsXtwhu5uV3i3ERwucAlgkvE+bKJ86VTEORwGa4tT+Ha8hSUv4mok4E0LR/yJ0DOAsgczq5DjUw4eS74siElA7zp4E0DTyq4U8DlBbcXXB5wuUHckXtXh5sQCId4u2ojJeV/5x8Vr7P7cGn7d890TzojM0YwIr2IqelF5PpyyfJmkZ2aTYY3A587DZ8nlRR3KqmuFDxuLx6XB6/LuXe53LjFjUtcH9wiraKhPaX4172Ff9Nm/Js2E9zTYQKpx4NnRBGe0aPwjBiBOzsbV3YW7qwsXFlZbC0tZca8eYg3BVdqKpLiRVJSEI8HcXvA7Wq/d7k84BJwuQDhcI2f/bsbqT7QSnVFMzUHWvC3fpDlU3xu0rNTyMj2Miw7hfTsFNIyvXh9brwpbjwpLjwpHufe68LtdeNyCy6XIG7B5XE5j4/cIr/jQGsLB/ds59Cu96nat4e6gweoP7ifoP+D7OdJSSVn+HDSsnLwZYzGl5lFakYmvszM9sfelBQ2v7eFufPm4fZ6cXm8eLxeXB4PLrcHt8uFuATEhYjgcjn3uFwIglPBcWo0zp04r0ceu8T1wX/ID+7o+Ccu0uXPvZu/fun2I6HrRn3/wtHb3/xRcQ0xsfi8S2SSKAPGdHheTFy+5w4uVeU3r5Wy6rltNPuDfHx6EUvnjOaMiYXkpHtje7L978Crd8N7TzrNByNmwkeuhzELYPQ8JKOw0+Z7S0qYsGBxv0/3XvV7PLr1UV4ofYGGQAOp7lROHXEq542/gJMLTubk/JMZnj48Zn94Gg7T+u67NPz1rzT89Xn8e/YA4C4sJG3WLNI+/RnSZs0kZfx4PMOGIW53j8cKlJSQE+UfS/NhP+Xbatm3tYqyLbU01LQCTjIoKM5kyoIRFBZnUjA6k/xRGaT4Bv5npOEw1eX7OLB9Gwe2b2X/+1upLt/X3qyTN3IUeSNHM27WHPJGjnZuo0aRmVcQ1c+7vLmNk2bNGXCcJyIn6Q7tZBBPiUwSTwHXichqYAFQf6L3RwRCYb792Aae2rCfj04Zxi0XnsL4wjgsLNZYCX/5PmxcDanZ8OHr4EOXQ8HEmJ8qrGH+UvoXfvfe79hQuYE0TxrnnHQOZ409i4UjF5LuTY/5Of2lpdT98Qnq//QnghUV4HaTseA08r7yZTLPPBPv6NEx/aMPh5WDuw9T+m4Vpe9WU7WvEYDUdA+jp+Yx79yxjJ6aR25ReszO23y4nood7zsJYfs2Kna8j7+l2TlvRgYjJ09j6ocXMXLyVEZMmoIvIzMm5zWmr+KWJETkUWAxUCgiZcAPAS+Aqt4HrAUuAHYAzcDyeMUyGMJhZcWj7/DMuxXccO5Url08MT7fXt57Cp76Bvib4CP/6tQcfNkxP42q8kr5K9z99t1sq93G2KyxfOfU77B00lKyU2J/vnBLC/XPPkvVH/7A4Y0bCHg8uGfPxPvZpXgmT6LW7eJgayuhda8SDoXQcIhwKEQ4HEZcLlwuNy73B7cUXxpenw+vz0dKqo/GinKqy/Y6TTEZmbQ1K2Xbail9t5q971XT1hREXMKICdksWDqBMSfnM2xsFq4B9qf4W5qpPbCfmgPl1O4vp2Z/GQd3baeuwvk+JOKi8KRxTDvjTEZOnsaoKdPIGzEKcdkUJnN8iOfopkuO8b4CX4/X+QfbT59/n2fereD7nziZKxdNiP0JQkF44Rb4x89h1Dz49H0wbGrsz4PTrHTHujt4+9DbjMkaw+2Lbuf88ed/0PbcRwF/G4011c6ttoam2poP7veX07i/nOamJvxucdqcTz7J2bG5Bv7+vHPrQsSFy+1CXG5Uw4SDIVTDvcax7YlHOzxzg/hwuX34MrPIKcghtyiPjJxs2hoy2LvRx4FtvvZE4031IXToNFNnZFA4FKa1saH91tJwmNbGRprr66it2E9TbU3HoMkuHMbwcROZeda5jJw8laIJk0jxxaeT25hYOOGWCj8evVVay89f2sHnPlTMFR8ZH/sTBFphzVdh2//BqVfCuf/pdDbHWF1rHT9/5+c8/v7j5KXmcfOsf+Oj+R+jrT7EttcqaGkM0NYUoK0lRNDv3PytbQRa6wm2RW7+ekKBBoL+egJt9QTbDhP0Nx11LpfLTWpISWlpwRcKUTByFDkzZpI1eQppWdn4MjNJzcj8oEaQ6iPF53M6Xl3ubr9pazhMa1Mb9ZVNVO6toWpfDdXlddQcqKO1sQm0DbfbT2a+i/SsMClpIaANf3MjLY0NHNh+gNbGxvZmn74SlwtfZhZpmVmkZWczbtY8py9h1GjyR44mZ8RIvCmx/70ZE0+WJAYoGArz3T9uZFROGrdceErsm5jaGuHRZbDnFTj/v2DBv8T08C2NfqrLGynZ9BrrtqwnrTmXr4Vvw9OUTtVLQdaESyDcgGoDGm4AbQRpgtBhwqEGwqGjP1Bd7jRcnmxc7ixc3iI87gxEMsGVibgyEckASUVFaAOCHiGY4aWlwUvtFg8paZGbz43bG8TlbsbtbsHlduYBhPxhgsEwIX+IgD9M82E/TXVtNNW1EWj7YNSRyyPkj8xnwtyx1LUd5MzzPkRhcSYud+81Ig2HCfr9+FtbCLS1EWhtIdDWiuqRwTVHRvEILreb1IxM0rKySPGlWTORGXIsSQzQH98u5/2Djfz3pfPITI3xjzPQCqsvgdJ/wGf+B2Z9od+HUlX8jcqOtw5RubeByn2HqdxbRXNdFRpuQMOHGa8teFJqcPMegcBhAi2H6ToqOTUjg6z8QrIKiskqGEZWQSFZhZH7gkIy8wvwpvoA8JeV0VjyMg0vvUTjm+8QUC/hEWNIXXQ27nkLCecOp7UxQGtT5NYYwN8aovmwn7qDzfhbg4SCTpNOOKiEw04sbu+RoaQuPClu0rO8FIzOYOwp+WTkppKV5yN/VAa5I9JxRxJCSckhhp8UXV+KuFztzUzGJDtLEgMQDIX52QvbmT0ml/NmjIjtwUNBWLMcdv8NPn1/VAlCQyHCLS1oayv+xlYO7muk9P0DVJRWUHPwIP6WajaHn0TD9aCH0XDn2WaelBQyModHPvCntieB7A6JICWt+9FMwdpa/Lt30/DXF2h+5x1a3lnvjEwCUsaNo/CLF5O5ZDHpH/pQr0NVey1fZGq5xGlynjHmaJYkBuCZdysor2vhh5+aHttmJlV45gbYthYuuBNmX9z+Vqi+ntbNm2nbsYO2nbsIlJfjP3CA6sZmKry51GTl0pjiIeBqdZIBHzS/iAppgTAZAT+uUBvgJ9OXztiCseTlDcOXk4s7MwNXWjqutDQk7EFqG+BwC+zYTaPfj7b5CTc2EqypJlRdQ/DQIfy7dxOqr28/j2fUSNLnzSNt7lwyPnIGqeNj008jIt1PwTTGxI0liQF4+NXdjC/M4GMnF8X2wK/dC28+BGd8k9D0S2h69jmaXnuN5nXraNu9mxavh5rMdCqHjaYu3UdrpgfNTAf8QCUeTzbZGXnk5xRTVJBFfmYGWSmpvLftbfa1baeqtoKRmsPMtElktkG4so7Qrn0cbmgg3NQE4d5HCQG4c3NxFxTgKSwk69xzSRk/npTx4/BNm4Z3RIxrVcaYhLEk0U87DjXw9t46vnfBtAGPpe984OcJ/98PaNAzqH+sgqZ/PZ0WwlQX5FE1ZgxV+TPxB490FodxuTPJyj+ZERMnMW7ONCafOhVfRucmod31u/nZ+l/wbOo/yPfls2Luv3PRpItwu45u9lFV1O8n3NyMBgIQCqHBIOJ2O8tYpKY6tQyP/dcxJhnYX3o/Pf5WGW6X8Om5xTE7ZttbJdSsvJr6XSOo9eznUHGYirnTaWxzZgATCOLyFJNVNIGxp0xnysJTKJ5a0OOyEDtqd/Dw5of5864/k+pO5Zzsc7j1E7eSmdLz7F0RcRJBqg3VNMZYkugXVeXPGw6waHIhw7IG/mHasmEDVffeQ8Xrr1OeP5z9M0fQEvIDgiuUhzf9Q4yYPJMpC6YzflYhOcN6XgojFA7x6v5XeXTro/y9/O+kedK49ORLuWLGFWx6fVOvCcIYY7qyJNEPm8rrKa9r4Zsfmzyg47Tt3EnFqjvZ+c4b7B5RSN20sYDgkpGkZE3jpJmnMnXhOMbNLCC1l8UBVZVttdv4y56/8PSup6loqiDfl8/X53ydZVOXkevLHVCcxpjkZUmiH57bXIHbJZwzvX8d1qHGRg789Cds/Mtadg4vxD9uBOLKwZM6k5FTFjLjo1OYOHdYr4mhuqWa9ZXreW3/a7xa/ipljWW4xMWHR36YG+bfwJIxS/C6Y7zqrDEm6ViS6IeXtlYy/6Q8ctNT+rxv/Qsv8Pqq29iamUZgZAHiGUV62inMKm7mlKu/Tu6IrPZtg+EgVS1VVDRVsK9hH3sb9rK9djtba7ZS3lgOQJonjVNHnMpXZ36Vs8eeTb4vP2blNMYYSxJ9dPBwK+8dOMyN503r037h5mbW/+BmXtu9g9YcH+IZRcbwM8gseI2mnP/h5TlLeXLbSmo21FDT6txqW2s7XaTTJS7GZo1lesF0Lp56MbOHzWZG4QxS3H1PVsYYEw1LEn30yvYqAD46ZVjU+1RueptnbllJpSeIpOYSyJnD6xM3sXPYnWjkAslZ+14kPy2fvNQ8xmSNYc7wORT4ChiePpwRGSMoziqmOLPYEoIxZlBZkuijf+ysIj8jhWkdmoW6o6psqNxAyf/cBesPE/KEEN88Kk5NpWB+mC81j2LsumcYMe+rDD/7VlLdNuTUGHP8sSTRB6rKP3dWs3BCfq8T6N6tepe73vwxY9aUk92QinhyGHbyuVz0zc+QnZMB1Tvhl2fCiHlw9o/Abb8GY8zxyT6d+qCstoX99a1cM6Gg2/fbQm389K2f8sjm3/Oll6bgakvB453A+dd8gylnTHI2CrbB45eD2wuffcAShDHmuGafUH3wZqlzlbH5444eQVTRVMGKF1fw/qGtXPHyHAJtNWSmzeaLq24ia1iHpqm//AAqNsKyRyEndrO1jTEmHixJ9MGbe2rJSvUwpahzf8Se+j1c8ZcraGlp5spXTqWt7SDDMz/EJffejMfXoaN5y59h3S9hwTUw7YJBjt4YY/rOLqPVB2+V1jJnbC7uDv0R5Y3lXPHcFQSDQb78yjzaWg4yOn0Wl/53lwRRtxf+dC2MnAMfv3XwgzfGmH6wJBGlZn+Q9w82MHdMbvtrjf5Grnn+GlpCLVy2bhHNh0sZ4ZnE5+/9Pq6UDgkiFIA1VzhLcH/+4bhcn9oYY+LBkkSUNu8/TFhhVnEu4Ix0+sGrP2Dv4b18vWwZtfvXk6Oj+cJPvoc7vcsCfC+shLJ1cOHdkD9h8IM3xph+siQRpQ376gCYNSYHgCd2PMHze5/nmtQrqXjzBVIo4As3XIl3+PDOO257Fv5xN8z/Ksz4zCBHbYwxA2NJIkqbyusZke1jeJaPQ82HWPXGKhbmfpiWp95A8HLBGQvIPvXUzjvV7YMnr4YRs+Dc2xITuDHGDIAliSi9t/8wM0ZnA/CTt36CP+Rn4SsnEQpWM8uTz4Sv/0vnHYJ+WLMcQkH4/K/A6xv8oI0xZoAsSUShNRBiZ2Uj00dms7l6M3/e9We+wpeo2/c62YFczvyP7yHuLpcCff4WKHsDlv4cCiYmJG5jjBkoSxJR2FbRQFhh+qgc7n3nXvI8+chf30fwcM6iU0kZN67zDlv+DP+8F079Gpzy6YTEbIwxsWBJIgpbKw4D4E0r55XyV/jC/qUE2vYxLpDGmGuv7rxx7R54MjIf4tz/GPRYjTEmlixJRGFrRQPpKW6eLVtNHgUE1m/Creks+eoluFI7zHk4si4TwBd+bfMhjDEnPEsSUXj/YAPjR4R5vvR5Ltp3LuFgJVNJJfcTn+i84V//Dfa/AxfdC3njEhKrMcbEUlyThIicJyLbRGSHiNzUzfs5IvK0iGwQkc0isjye8fTXtooGUnPfgIAL2bIDt6Zx+lWXIdJhufCt/wev3wcLroaTP5W4YI0xJobiliRExA3cC5wPTAcuEZHpXTb7OvCeqs4GFgM/FpHj6tJrtU1+qhrbOKR/5xOHziEcPMiEcArZZ5/9wUb1ZZF+iNnw8ZWJC9YYY2IsnjWJ04AdqrpLVf3AamBpl20UyBLnK3kmUAME4xhTn+2obMSdVkp9oIJh25oRvHz44os+qEWEQ/DHf3HWZ/qcrctkjBla4rlU+GhgX4fnZcCCLtvcAzwF7AeygItVNdz1QCJyFXAVQFFRESUlJf0OqrGxsU/7l+wL4MnewMSqaQRbdjOsNY3NubkQOcbY0jVMKP07W6euoGLTPjoX+fjS17IPJclcdkju8lvZSwZ0jKiShIikAWNVdVsfjt3d9T21y/NzgfXAWcBE4K8i8oqqHu60k+r9wP0A8+fP18WLF/chjM5KSkroy/4vP/UuXv+7nLH+o8AmTps/i2kf/7jzZsUm+NtqmL6UaZ9fyTTp+ZKmx4O+ln0oSeayQ3KX38q+eEDHOGZzk4h8CueD/NnI8zki8lQUxy4DxnR4XoxTY+hoOfBHdewAdgPTojj2oNlUvYlUDeCu3k9qKIMJl33ZeSMUgCevgbQ8+MRP4DhPEMYY0x/R9EncgtO/UAegquuBcVHs9wYwWUTGRzqjl+E0LXW0FzgbQESKgKnAriiOPWhKW97glLLpaKiacanppIwd67zx6k+dmsQnfwIZ3V/z2hhjTnTRNDcFVbVe+vhNWVWDInId8BzgBh5S1c0icnXk/fuAfwd+JSKbcJqnblTVqj6dKI78wTBN7nc5Zd9MoIZ5553lvFG9E15e5Sy5cfInExqjMcbEUzRJ4l0R+SLgFpHJwArgH9EcXFXXAmu7vHZfh8f7gXOiD3dwbTywjzR3Hd7mg6SFMii68EJQhbX/zxnFdN7tiQ7RGGPiKprmpm8ApwBtwKPAYeBbcYzpuPHCnleZemAKGq5jTEY27pwc2PYM7HwRltwMWSMSHaIxxsTVMWsSqtoM3AzcHJkgl6GqrXGP7DjwTuWbTN9XANQxa/FCp7P6rz+Awqlw6pWJDs8YY+IumtFNj4hItohkAJuBbSJyQ/xDS7zSxnfJbGrBG05n1Cc/Be/8Fqp3wMdvBXc8p5gYY8zxIZrmpumReQsX4fQvjAW+HM+gjgdVLVVk1aagwQoK3el48rKdzuoxC2DKeYkOzxhjBkU0ScIrIl6cJPEnVQ1w9KS4IWf9ofWcvG8UEGbKlPFOLaJhPyz+rs2JMMYkjWiSxC+BPUAG8DcROQmn83pI23BoA6OrfICLKed+DP7xcxg9HyYsTnRoxhgzaI6ZJFT1blUdraoXRGZGlwJLBiG2hHrnwCY8/kZ84Qyysg45V5w7Y4XVIowxSSWajusCEblbRN4WkbdE5GdAziDEljBhDVO7px4NVZLjSUfeeQhyxsA0mzhnjEku0TQ3rQYqgc8Cn4s8/t94BpVopYdLGV+WDyhjx+TD7r/B/OXgcic6NGOMGVTRJIl8Vf13Vd0duf0IyI1zXAm1pXoLo6vTADhl9GEQN8y5NMFRGWPM4IsmSbwkIstExBW5fQH4v3gHlkhba7eS2hrEo+nktf0NJp1ts6uNMUkpmiTxL8AjgD9yWw38q4g0iMiQHOW0c/8eNFSNL5yCq3k/zLo40SEZY0xCRLMsR9ZgBHI8aXq/EbSFPLeAx2eT54wxSSvaK9PlAZMB35HXVPVv8Qoqkerb6hm+1/mxTPHtgkkfg9TMBEdljDGJccwkISJXAt/EubLcemAh8BrOJUeHnF31uyg47AWEibmlMPX6RIdkjDEJE02fxDeBU4FSVV0CzMUZBjsk7azbidcfwqNpZOQEYPJxe7kLY4yJu2iSROuRpcFFJFVVt+JcZnRI2l29B8KHSQu78RfNhMxhiQ7JGGMSJpo+iTIRyQWeBP4qIrXA/ngGlUgH91QyMnyYfNy4Jg351UeMMaZX0Yxu+nTk4S0i8hLOkhzPxjWqRNrWAMBIbyUpkxYnNhZjjEmwHpOEiOR38/KmyH0mUBOXiBIoEA6QWREGYHTmARizMMERGWNMYvVWk3gL57oRHZc9PfJcgQlxjCshKhoryGxxAy78RYWQkp7okIwxJqF6TBKqOn4wAzke7G3YizsYxK1pHCyaw5REB2SMMQkWzeimpFFavxcNN+FTF4eHzU10OMYYk3BRzbhOFgfKD0L4MNl4CI2cl+hwjDEm4XqsSYhI0jU3NbxfCYTJdjeTVZR0xTfGmKP01ty0BkBEXhikWBLOVeoMfyUlwIictMQGY4wxx4HemptcIvJDYIqI/GvXN1X1rviFlRje2iAABzNzWJjtO8bWxhgz9PVWk1gGtOIkkqxubkNKU6CJlDYAN5uzJpOX7k10SMYYk3C9DYHdBtwhIhtV9ZlBjCkhKpoqcIVCuPBRkTMNETn2TsYYM8RFMwT2HyJyl4i8Gbn9WERy4h7ZINvfuB8Nt5KibtpyrNPaGGMguiTxENAAfCFyOww8HM3BReQ8EdkmIjtE5KYetlksIutFZLOIvBxt4LG2v+ogGm7AJ2Hys+0iQ8YYA9HNk5ioqp/t8PxWEVl/rJ1ExA3cC3wcKAPeEJGnVPW9DtvkAr8AzlPVvSIyvC/Bx1JlaQXgxytehmdZp7UxxkB0NYkWEfnIkScicgbQEsV+pwE7VHWXqvqB1cDSLtt8Efijqu4FUNVD0YUde63vVwAQcLsoyk5NVBjGGHNciaYmcTXwmw79ELXAZVHsNxrY1+F5GbCgyzZTAK+IlOCMmPqZqv6m64FE5CrgKoCioiJKSkqiOH33Ghsbu93fv6+GVKAmJY3G/XsoKSnv9zmOVz2VPRkkc9khuctvZS8Z0DGiuZ7EBmC2iGRHnh+O8tjdDQ/Sbs7/IeBsIA14TUT+qarvd4nhfuB+gPnz5+vixYujDOFoJSUldLf/m797DID300Zy8fxZLJ5W1O9zHK96KnsySOayQ3KX38q+eEDHiHqBP1U93IcEAU7NYUyH58UcfUW7MuBZVW1S1Srgb8DsPpwjZsSvgIsN6VMYlml9EsYYA/FdBfYNYLKIjBeRFJzJeU912eZPwCIR8YhIOk5z1JY4xtSt1mArrlAYNz4qPAUUZqUMdgjGGHNcitsqsKoaFJHrgOcAN/CQqm4Wkasj79+nqltE5FlgIxAGHlDVd+MVU08qWypBA7jxAEJhpnVcG2MMRJEkIt/wvw2MVdWvichkYKqq/vlY+6rqWmBtl9fu6/J8FbCqT1HHWGVTJaqteHCRm+7F67bLbBhjDETX3PQw0AZ8OPK8DPhR3CJKgIrqSjTchMdltQhjjOkomiQxUVX/CwgAqGoL3Y9cOmFVvr8TCCFuD4WZ1h9hjDFHRJMk/CKSRmT4qohMxKlZDBlNO8sAaHOnWk3CGGM6iKbj+ofAs8AYEfk9cAZweTyDGmzBQ40AHPJkWpIwxpgOoplM91cReRtYiNPM9M3InIYhI9zgXGxoq7uQj2RYc5MxxhwR7TCe0TjDWFOAM0XkM/ELKQFawwBsSx9DgdUkjDGmXTRDYB8CZgGbceYygNM/8cc4xjW4ggAe6lJzKLCOa2OMaRdNn8RCVZ0e90gSKRzGTSqIUGDNTcYY0y6a5qbXRGTIJomWYAuqIVyRfJlvScIYY9pFU5P4NU6iqMAZ+iqAquqsuEY2SGpbalD14xYnXxZkWJ+EMcYcEU2SeAj4MrCJD/okhozq2r2otuB2ZeJxCdlpcVvOyhhjTjjRfCLuVdWuq7cOGYfK9oC2Iu4c8jJSEBlSk8mNMWZAokkSW0XkEeBpOsy0VtUhMbqp+v3dAIQ9KdZpbYwxXUSTJNJwksM5HV4bMkNgmw7UANDqSbNOa2OM6SKaGdfLByOQRPHXtABQ6cklz5KEMcZ00mOSEJHvqOp/icjPOfra1KjqirhGNkiCTSEA9rryGJ9uScIYYzrqrSZx5DKibw5GIIkSbnXy305XJvOsJmGMMZ30mCRU9enIw2ZVfbzjeyLy+bhGNYjCztp+VPryyE/3JjYYY4w5zkQz4/q7Ub52YgqDkELQ7bU+CWOM6aK3PonzgQuA0SJyd4e3soksiTcUaFhx4SSHXOuTMMaYTnrrk9iP0x9xIfBWh9cbgOvjGdRgUg3iFjcAedbcZIwxnfTWJ7EB2CAij6hqYBBjGjStzTVAAOFIkrCahDHGdHTMPomhmiAA6uv3ElY/LpeTJHKtJmGMMZ1Ee2W6Iammugy0DZfLhcclZKba4n7GGNNRr0lCRNwismqwghlsVWX7gCDqdpObbov7GWNMV70mCVUNAR+SIfrpWbfvIAAhT6o1NRljTDeiaV95B/iTiDwONB15cSisAttUeRiAFk+ajWwy5gQTCAQoKyujtbW11+1ycnLYsmVLr9sMJT6fj+LiYrze2HymRZMk8oFq4KwOrw2JVWBb6p3/XLWeLHLSbGSTMSeSsrIysrKyGDduXK9NxQ0NDWRlZQ1iZImjqlRXV1NWVsb48eNjcsykXgU2EFncr8qVbTUJY04wra2tx0wQyUZEKCgooLKyMmbHPOboJhGZIiIviMi7keezROT7MYsggUJtztVY95NufRLGnIAsQRwt1j+TaIbA/g/OWk0BAFXdCCyL5uAicp6IbBORHSJyUy/bnSoiIRH5XDTHjZVwwPlhlrkzbUkOY4zpRjRJIl1V13V57ZhrN4mIG7gXOB+YDlwiItN72O4O4LkoYompcAhAqE/JJCfNahLGmL6pq6vjF7/4RdTbL1q0iDlz5jBnzhxGjRrFRRdd1O12bre7fbsLL7wwRtH2TzQd11UiMpHIhYci3/YPRLHfacAOVd0V2W81sBR4r8t23wD+AJwabdCxopEVYHG5LEkYY/rsSJK49tpro9r+lVdeaX/82c9+lqVLl3a7XVpaGuvXr49FiAMWTZL4OnA/ME1EyoHdwKVR7Dca2NfheRmwoOMGIjIa+DTOyKkek4SIXAVcBVBUVERJSUkUp+9eY2Nj+/6qYSTyIyjdvoWS2vf7fdwTQceyJ5tkLjsMzfLn5OTQ0NAAwB1/2cnWg43dbqeq/Wqnn1aUyY3nTOx1m29/+9vs3LmTWbNmsWTJEn70ox9FdeyGhgZefPFF7r777vYydLdNf7W2tlJSUhKT33s0SaJUVT8mIhmAS1Wjjby730rXy6D+FLhRVUO9/RJV9X6cRMX8+fN18eLFUYZwtJKSEhYvXgyqrP/lb3GJ8yM4c+GHmFWc2+/jngjay56EkrnsMDTLv2XLlvahrd4UL263u9vtQqFQj+/1xpviPebQ2R//+Mds27aNjRs30tDQwKJFi7rd7pFHHmH69A9a25944gnOPvtsRo8e3e32ra2tLFmyBI/Hw0033dRjs1RPfD4fc+fOjcnvPZoksVtEngX+F3ixD8cuA8Z0eF6Ms/x4R/OB1ZEEUQhcICJBVX2yD+fpn2ArSghXpFvGmpuMOXH98FOn9PjeYM2TyMrKirqJ6NFHH+XKK6/s8f29e/cyatQodu3axVlnncXMmTOZOLH3Wk28RJMkpgKfwml2elBE/gysVtW/H2O/N4DJIjIeKMcZEfXFjhuoavtsDxH5FfDnQUkQQLi5BiXQXpPItcl0xpgBiLYmUV1dzbp163jiiSd6PNaoUaMAmDBhAosXL+add945fpOEqrYAjwGPiUge8DPgZaDX+puqBkXkOpxRS27gIVXdLCJXR96/b6DBD0RDwwFUA4grFRHI9NkKsMaYvsnKymrvO4i2JvH444/zyU9+Ep/P1+37tbW1pKenk5qaSlVVFa+++irf+c53Yhl2n0S1VLiIfFREfgG8DfiAL0Szn6quVdUpqjpRVf8j8tp93SUIVb1cVdf0IfYBqa8/CNoGbhdZqR7cLpuUY4zpm4KCAs444wxmzJjBDTfcENU+q1ev5pJLLun02ptvvtne/LRlyxbmz5/P7NmzWbJkCTfddFOn/ozBdsyvzyKyG1iPU5u4QVWbet/jxFC1fz8QIuxyk2OzrY0x/fTII4/0afvuRhvNnz+fBx54AIDTTz+dTZs2xSK0mIimjWW2qh6OeySDrL6iCoCQy2Od1sYY04NompuyReQJETkkIgdF5A8iUhz3yOKssbIOgDZ3Ctk+SxLGGNOdaJLEw8BTwCicCXJPR147oTXXNwPQ5E6zmoQxxvQgmiQxTFUfVtVg5PYrYFic44q7tkZn+akaSbckYYwxPYgmSVSJyJci17t2i8iXcC5CdEILtB65lkQG2ZYkjDGmW9Ekia/iDHmtwFnY73OR105owTZnhZCDkkG2zZEwxphuHTNJqOpeVb1QVYep6nBVvUhVSwcjuHgKRRY7r/ZmWU3CGNMvfV0q/J577mHSpEmICFVVVe2vl5SUkJOT0748+MqVK7vdf/fu3SxYsIDJkydz8cUX4/f7B1yGY4lqMt1QpE5rEw0pGdYnYYzpl74miTPOOIPnn3+ek0466aj3Fi1axPr161m/fj3/9m//1u3+N954I9dffz3bt28nLy+PBx98sN+xRytp21nCCuBFXW4bAmvMie6Zm6Ci+wloaaEguPvxUTdiJpx/e6+b3HTTTezcuZM5c+bw8Y9/nFWrVvW6/dy5c/seR4Sq8uKLL7ZP3rvsssu45ZZbuOaaa/p9zGgkbZLQsCI4ySE7LWl/DMaYAbj99tt59913Wb9+PQ0NDcyZM6fb7bouFd6d1157jdmzZzNq1CjuvPNOTjml88q21dXV5Obm4vE4n1fFxcWUl5fHpBy9iWZZjiLgP4FRqnp+5BKkH1bV+Ndz4khVcYmzRqHVJIw5wfXyjb/lOFwqvKt58+ZRWlpKZmYma9eu5aKLLmL79u2dtlHtejke+nUxpb6Kpk/iVzgruY6KPH8f+Fac4hk0iiKRhWyzLEkYYwboSE2iu9t773W9anNn2dnZZGZmAnDBBRcQCAQ6dWwDFBYWUldXRzDojLopKytrX1I8nqJpZylU1cdE5LvQvgR4KM5xxVewzS44ZIwZsP4sFd6diooKioqKEBHWrVtHOBymoKCg0zYiwpIlS1izZg3Lli3j17/+dY/XyI6laGoSTSJSQOTSoyKyEKiPa1Tx1taIahARF1634PMm7SAvY8wA9HWp8Lvvvpvi4mLKysqYNWtW+/Lga9asYcaMGcyePZsVK1awevXq9qakCy64gP37nYt63nHHHdx1111MmjSJ6upqrrjiivgVLiKamsS/4qzdNFFEXsVZkuNzcY0qzkKttUAAEWdxv8Fo1zPGDE19WSp8xYoVrFix4qjXr7vuOq677rpu91m7dm374wkTJrBu3bq+BzkA0VyZ7m0R+SjOZUwF2KaqgbhHFkeNjVWoBkBcZNlsa2OM6VE0o5u+0uWleSKCqv4mTjHFXW3VQSCAusQ6rY0xphfRfI0+tcNjH3A2zmVMT9gkUbO/AoCQy21zJIwxphfRNDd9o+NzEckBfhu3iAbB4UpnEduAzbY2xphe9WdYTzMwOdaBDKbGWudqrK3isT4JY4zpRTR9Ek8TGf6Kk1SmA4/FM6h4az3cCkCDpFifhDHG9CKamsSdwI8jt9uAM1X1prhGFWdtTc7grFpJs+YmY0y/9XUV2EsvvZSpU6cyY8YMvvrVrxIIOJ9Fv//975k1axazZs3i9NNPZ8OGDd3uf/nllzN+/Pj2mdz9nbzXF9FcT+LlDrdXVbUs7lHFWaDNmTBe786w5iZjTL/1J0ls3bqVTZs20dLSwgMPPADA+PHjefnll9m4cSM/+MEPuOqqq3o8xqpVq9qXFO9pQcFY6vETUkQa+KCZqdNbgKpqdtyiirOgPwxAvdsuOGTMUHDHujvYWrO12/dCoRBut7vPx5yWP40bT7ux1236ulT4BRdc0P74tNNOo6zM+c59+umnt7++cOHC9tePBz0mCVWN/7KJCRKMTAWs9mZbTcIY02/9XSo8EAjw29/+lp/97GdHbfvggw9y/vnn93jOm2++mZUrV3L22Wdz++23k5qaOuBy9CbqT0gRGY4zTwJwLmsal4gGQTjkVJDqU7IsSRgzBPT2jb/hOFwq/Nprr+XMM89k0aJFnV5/6aWXePDBB/n73//e7X633XYbI0aMwO/3c9VVV3HHHXf0eBW7WIlmdNOFOJ3Wo4BDwEnAFuCU3vY7nmkYwE3A7bWOa2NMTDQ0NBz1oX9Ex5rErbfeSmVlJb/85S87bbNx40auvPJKnnnmmaNWgD1i5MiRAKSmprJ8+XLuvPPOGJage9F8jf53YCHwvKrOFZElwCXxDSu+wmFFIkW3JGGM6a++LhX+wAMP8Nxzz/HCCy/gcn0wbmjv3r185jOf4be//S1Tpkzpcf8DBw4wcuRIVJUnn3ySGTNmxKQcvYlmCGxAVasBl4i4VPUlYE58w4ovVRBxkoQ1Nxlj+quvS4VfffXVHDx4kA9/+MPMmTOHlStXArBy5Uqqq6u59tprmTNnDvPnz2/fp+NS4ZdeeikzZ85k5syZVFVV8f3vfz8+Besgmk/IOhHJBP4G/F5EDgHBaA4uIucBPwPcwAOqenuX9y8FjjQmNgLXqGr3A4RjSPWDq9JlWpIwxgxAX5YKP3JVua4eeOCB9uGwXXVcKvzFF1/sW3AxEE1NYinOUhzXA88CO4FPHWsnEXED9wLn48zSviRyfeyOdgMfVdVZOM1a90cfev8pYQQ3aV43XrddcMgYY3oSzdfoq4DHI5Poft2HY58G7FDVXQAishon4bRf7FVV/9Fh+38CxX04fr8duXSpNTUZY0zvovmUzAaeE5EaYDWwRlUPRrHfaGBfh+dlwIJetr8CeKa7N0TkKpxkRVFRESUlJVGcvnuNDQ2ohkC8uMP+AR3rRNPY2JhU5e0omcsOQ7P8OTk57Z3GvQmFQlFtN5S0trZSUlISk997NEuF3wrcKiKzgIuBl0WkTFU/doxdu7smaHczuImMmLoC+EgPMdxPpClq/vz5unjx4mOF3aO/vfAcEERIoSg/h8WLz+j3sU40JSUlDORndyJL5rLD0Cz/li1bopr/MFjzJI4nPp+PuXPnxuT33pf2lkNABVANDI9i+zJgTIfnxcD+rhtFks8DwPmRUVTxFWh2Ll3qEluSwxhjjuGYvbYico2IlAAvAIXA1yIdzcfyBjBZRMaLSAqwDHiqy7HHAn8Evqyq7/c1+P4ItNYDAVRs+KsxxhxLNEN7TgK+paqnqOoPVfW9Y+4BqGoQuA54DmeG9mOqullErhaRqyOb/RtQAPxCRNaLyJv9KEOftDTUA0pYhGxLEsaYAejrKrBXXHEFs2fPZtasWXzuc5+jsbERcIblr1ixgkmTJjFr1izefvvtbvffvXs3CxYsYPLkyVx88cX4/f6YlKM30SwVfpOqru/PwVV1rapOUdWJqvofkdfuU9X7Io+vVNU8VZ0Tuc3v/YgD56+vByAkYhccMsYMSF+TxE9+8hM2bNjAxo0bGTt2LPfccw8AzzzzDNu3b2f79u3cf//9XHPNNd3uf+ONN3L99dezfft28vLyePDBB2NSjt4k3VfpYEMTAH6Xy2oSxgwRFf/5n7Rt6X6p8GAoRE0/lgpPPXkaI773vV636etS4dnZzhUWVJWWlhZEnPE9f/rTn/jKV76CiLBw4ULq6ural+A4QlV58cUX2yfvXXbZZdxyyy09JpRYSbpPyUBzm3MvbqtJGGMGpD9LhS9fvpy1a9cyffp0fvzjHwNQXl7OmDEfjPMpLi6mvLy8U5Korq4mNzcXj8fTaZt4S74k0epcTKLFnWId18YMEb194z/elgp/+OGHCYVCfOMb3+B///d/Wb58OapHzw44Uss4Ippt4iHpPiXDkUuXNrlSrSZhjImZaJcKB3C73Vx88cWsWrWK5cuXU1xczL59H8w9LisrY9SoUZ2OUVhYSF1dHcFgEI/H0+028ZB0SSIYOHLp0nTrkzDGDEhflgpXVXbu3MmkSZNQVZ5++mmmTZsGwIUXXsg999zDsmXLeP3118nJyenU1AROrWHJkiWsWbOGZcuW8etf/5qlS5fGrWxHJN3qdhp0qmw1dn1rY8wA9WWpcFXlsssua1/q+8CBA+1XlbvggguYMGECkyZN4mtf+1qnEVMdlwq/4447uOuuu5g0aRLV1dVcccUV8StcRNJ9lQ5HVuqtdOdZn4QxZsCiXSrc5XLx6quvdvueiHDvvfd2+17HpcInTJjAunXr+h7kACRdTSLstDZR6821moQxxhxD0iUJDSngxu9NITPFahLGGNOb5EsSCoKHTJ8Xlyv+w8eMMeZEloRJQhHxkG3DX40x5piSM0ngsf4IY4yJQvIlCcIILnLSrD/CGGOOJWmThM22NsYMVKyWCl+1ahVz5sxhzpw5zJgxA7fbTU1NzVH7X3755YwfP75922iWARmopEsSaAgRIceam4wxAxSrpcJvuOEG1q9fz/r167ntttv46Ec/Sn5+frfHWLVqVfu2PS0oGEtJ1+aiBAGvJQljhpBXHnufqn2N3b4XCoVw92Op8MIxmSz6wpRet4nVUuEdPfroo1xyySV9jjdeki9JaAAkzUY3GWMGLFZLhR/R3NzMs88+217D6M7NN9/MypUrOfvss7n99ttJTU2NWXm6k1RJIhAIAAFEIC/DkoQxQ0Vv3/hPhKXCj3j66ac544wzemxquu222xgxYgR+v5+rrrqKO+64o339p3hJqj6J6qoqQFHBmpuMMTF1pCbR3e29997rtO2RpcL/8Ic/dHp99erVvTY1jRw5EhEhNTWV5cuXD8o6TklVk6gpc9ZrtyRhjImFWC0VDlBfX8/LL7/M7373ux6PceSSpqrKk08+yYwZM2JWlp4kVU2i7oCz3G7YreSlpyQ4GmPMiS5WS4UDPPHEE5xzzjlkZGR02q/jUuGXXnpp+/5VVVV8//vfj32hukiqmkRjZSUAIRFy060mYYwZuFgsFQ7OHIjLL7/8qNc7LhX+4osv9jm+gUqqmkRzbR0AQbeQl2E1CWOMOZakShJtDc446oDHQ1ZqUlWijDGmX5IqSfibWwEIpaZ2O4nFGGNMZ0mVJAKtbQBol44hY4wx3UuqJBFscy5w7cntfqKKMcaYzpIqSYQCIcBFWsGIRIdijDEnhKRKEuFQCPCSk1+U6FCMMUNAX1eBPeIb3/gGmZmZ7c9LSkrIyclpn6G9cuXKbvfbvXs3CxYsYPLkyVx88cX4/f5+xx6t5EoS4TAu8ZKfn5foUIwxQ0B/ksSbb75JXV3dUa8vWrSofQnwntZjuvHGG7n++uvZvn07eXl5PPjgg/0Ju0/iOg5URM4Dfga4gQdU9fYu70vk/QuAZuByVX07XvGEwyEED0WFhfE6hTEmAV761f0cKt3V7XuhYAi3p+9LhQ8/aQJLLr+q1236ulR4KBTihhtu4JFHHuGJJ57oUzyqyosvvtg+ee+yyy7jlltu4ZprrunTcfoqbklCRNzAvcDHgTLgDRF5SlU7rnR1PjA5clsA/HfkPi5UwwhuJo0ZGa9TGGOSSF+XCr/nnnu48MILGTny6M+g1157jdmzZzNq1CjuvPNOTjnllE7vV1dXk5ubi8fjfGwXFxdTXl4e8zJ1Fc+axGnADlXdBSAiq4GlQMcksRT4jaoq8E8RyRWRkap6INbBPPStbxMKV+OSLIpy0mN9eGNMAvX2jf94WSp8//79PP7445SUlBz13rx58ygtLSUzM5O1a9dy0UUXsX379k7bOB+TnQ3GfK94JonRwL4Oz8s4upbQ3TajgU5JQkSuAq4CKCoq6vaHfCya7nS/qGi/9h8KGhsbrexJaiiWPycnp30F1t6EQqGotuuPxsZGwuEwDQ0NNDQ0cN5553W73YMPPsiePXvYvn07EydOBJwLDE2YMIENGzYgIqgqDQ0NLFq0CL/fz549eygoKGg/RmpqKrW1tdTW1uLxeNi2bRvDhw/vtmytra2UlJTE5PcezyTRXYrrmgqj2QZVvR+4H2D+/Pm6ePHiPgczvqCANT/6Lt50L/3ZfygoKSmxsiepoVj+LVu2RFVDiGdNYuTIkTQ1NZGVlUVWVhYbN27scdtTTz2Vz3/+8+3PMzMz2bXL6UepqKigqKgIEWHdunWoKieddNJRNYWzzjqL5557jmXLlrFmzRo++9nPdls2n8/H3LlzY/J7j+fopjJgTIfnxcD+fmwTExm5zrVlU3PiX+00xiSHviwV3ps1a9YwY8YMZs+ezYoVK1i9enV7gui4VPgdd9zBXXfdxaRJk6iuruaKK66ISTl6E8+axBvAZBEZD5QDy4AvdtnmKeC6SH/FAqA+Hv0RAPmjipl73qdoy7bZ1saY2Il2qfCuGhsb2x9fd911XHfddd1u13Gp8AkTJgzK1eg6iltNQlWDwHXAc8AW4DFV3SwiV4vI1ZHN1gK7gB3A/wDXxisel9vNWcv/hbSCYfE6hTHGDDlxnSehqmtxEkHH1+7r8FiBr8czBmOMMf2XVDOujTFDS3fDQpNdrH8mliSMMSckn89HdXW1JYoOVJXq6mp8Pl/MjmmXZzPGnJCKi4spKyujMnLt+p60trbG9EPzeOfz+SguLo7Z8SxJGGNOSF6vl/Hjxx9zu5KSEubOnTsIEQ1N1txkjDGmR5YkjDHG9MiShDHGmB7JiTYyQEQqgdIBHKIQqIpROCcaK3vySubyW9kdJ6lqn2cTn3BJYqBE5E1VnZ/oOBLByp6cZYfkLr+VfWBlt+YmY4wxPbIkYYwxpkfJmCTuT3QACWRlT17JXH4r+wAkXZ+EMcaY6CVjTcIYY0yULEkYY4zpUdIkCRE5T0S2icgOEbkp0fHEk4iMEZGXRGSLiGwWkW9GXs8Xkb+KyPbIfV6iY40XEXGLyDsi8ufI82Qqe66IrBGRrZH/Ax9OlvKLyPWR//PvisijIuIbymUXkYdE5JCIvNvhtR7LKyLfjXwGbhORc6M5R1IkCRFxA/cC5wPTgUtEZHpio4qrIPBtVT0ZWAh8PVLem4AXVHUy8ELk+VD1TZwrIh6RTGX/GfCsqk4DZuP8HIZ8+UVkNLACmK+qMwA3zmWTh3LZfwWc1+W1bssb+QxYBpwS2ecXkc/GXiVFkgBOA3ao6i5V9QOrgaUJjiluVPWAqr4dedyA8yExGqfMv45s9mvgooQEGGciUgx8Anigw8vJUvZs4EzgQQBV9atqHUlSfpyVrdNExAOkA/sZwmVX1b8BNV1e7qm8S4HVqtqmqrtxLht92rHOkSxJYjSwr8PzsshrQ56IjAPmAq8DRap6AJxEAgxPYGjx9FPgO0C4w2vJUvYJQCXwcKS57QERySAJyq+q5cCdwF7gAFCvqn8hCcreRU/l7dfnYLIkCenmtSE/9ldEMoE/AN9S1cOJjmcwiMgngUOq+laiY0kQDzAP+G9VnQs0MbSaV3oUaXtfCowHRgEZIvKlxEZ1XOnX52CyJIkyYEyH58U41dAhS0S8OAni96r6x8jLB0VkZOT9kcChRMUXR2cAF4rIHpxmxbNE5HckR9nB+b9epqqvR56vwUkayVD+jwG7VbVSVQPAH4HTSY6yd9RTefv1OZgsSeINYLKIjBeRFJzOm6cSHFPciIjgtElvUdW7Orz1FHBZ5PFlwJ8GO7Z4U9Xvqmqxqo7D+T2/qKpfIgnKDqCqFcA+EZkaeels4D2So/x7gYUikh75Gzgbpz8uGcreUU/lfQpYJiKpIjIemAysO9bBkmbGtYhcgNNW7QYeUtX/SGxE8SMiHwFeATbxQbv893D6JR4DxuL8QX1eVbt2eg0ZIrIY+H+q+kkRKSBJyi4ic3A67VOAXcBynC+EQ778InIrcDHOCL93gCuBTIZo2UXkUWAxzpLgB4EfAk/SQ3lF5Gbgqzg/n2+p6jPHPEeyJAljjDF9lyzNTcYYY/rBkoQxxpgeWZIwxhjTI0sSxhhjemRJwhhjTI8sSRhjjOmRJQljYkwc9rdlhgT7j2xMDIjIuMi1G34BvE3n5Q+MOWHZZDpjYiCy2u4u4HRV/WeCwzEmZqwmYUzslFqCMEONJQljYqcp0QEYE2uWJIwxxvTIkoQxxpgeWce1McaYHllNwhhjTI8sSRhjjOmRJQljjDE9siRhjDGmR5YkjDHG9MiShDHGmB5ZkjDGGNOj/w8fIGqFmMH+2AAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "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",
+    "\n",
+    "for i, t_i in enumerate(t) :\n",
+    "    if (i < num_points_t) and (i % 20 == 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(1.0,9.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": 7,
+   "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": 8,
+   "id": "0cd0795a",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEGCAYAAACpXNjrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyK0lEQVR4nO3de4CUdb348fdnbnvf5SbLLgsCLl6Qm4RAEAURpVTgryzxZ4VI8UtDOp4yzUspv07CwSw72oWgIhMpMRX7oXQU18pMvCGK2AFBZLmzsLALe5uZz++PmV1md2eXZ3bnsjv7eZ0zzjz3z7NL89nv93m+n0dUFWOMMSaeXKkOwBhjTPqx5GKMMSbuLLkYY4yJO0suxhhj4s6SizHGmLjzpDqAZOrXr58OGTIkacc7ufcQma4c6r1VeOtzqA/W4REPbpcXAH9GkOx++Rw8dZDKukou7HNh0mJLtFOnTpGTk5PqMJLKzrln6Inn/Nprrx1V1XNi2aZHJZchQ4bw6quvJu14Gxf/hIuzP8Se4hco2juRPad20cfXh76ZAwA4cN4pLv3aZfzktZ+wettqXv1K8mJLtLKyMqZNm5bqMJLKzrln6InnLCJ7Yt3GusVSyF8fGmOU4cnAr378QX+KIzLGmPiw5JJCAX8ouWS6MwGoD9SnMhxjjIkbSy4p5G8IApDtyQbgtP90KsMxxpi4seSSdNL0KRDuBcv2hpNLgyUXY0x6sOSSQo3dYrneXACqG6pTGY4xxsSNJZcU8vtD3WK5vlByOdVwKpXhGGNM3FhySRpBpPmclt1i1fXWcjHGpAdLLknT+tEGLbvFTvmt5WKMSQ+WXJItovWi/tBEjjc02vdUvSUXY0x6sOSSAvWBWgBcgdC0XdA3xqQbSy4p0BCsA8AVup5PhjsDt7jtgr4xJm1YckkgaWN+fWNyCYSaLiJCjjfHkosxJm1YckmBhmCozIubIMFA+HZkb651ixlj0oYll2SQ5neKBTVAIOjHK0pDfSi55Pis5WKMSR89quR+VyAIgov6YC1egU2/207f4hxyPDnWcjHGpA1LLknScpRLQ7AOrwi73jjC6RN15IzL4WTdyZTEZowx8WbdYinSEKzFK6Eff01Vg11zMcaklZQmFxG5TET+JSI7ReTWKMtFRH4aXr5VRMa1WO4WkTdE5M/Jizo+/IE6vC43ADXVoeRigyiNMekiZclFRNzAg8DlwAjgahEZ0WK1y4Hh4ddC4Octln8T2J7gUBMi4K/FJ6FeyfoaP9luu+ZijEkfqWy5TAB2quouVa0H1gJzWqwzB/idhvwT6CUiRQAiUgJ8GliZzKDjJeCvwevyMjzDxRCfixx/Pqf9pwlqMNWhGWNMp6Xygv5AYG/EdDkw0cE6A4EDwE+A7wB57R1ERBYSavVQWFhIWVlZZ2KOSSAYShT19f5m8xUI+mvwunwM8wU4FfDz2u5KADY+v5EsV1bSYkyU6urqpP6suwI7556hJ55zR6QyuUQbwN7ypqqo64jIZ4DDqvqaiExr7yCqugJYATB+/HidNq3d1ePqL3/aAoA3w9tqmdRX4xY3brebhmAdpYXnw364ZOIlDMgZkLQYE6WsrIxk/qy7AjvnnqEnnnNHpLJbrBwYFDFdAux3uM4UYLaIvE+oO+3jIvL7xIUaL2dypdafeaRxpjsDX20oAdlASmNMOkhlcnkFGC4iQ0XEB8wF1rdYZz3wlfBdY5OAE6p6QFW/q6olqjokvN0mVf1SUqOPidCsUSYgEc9u8bp8eKpCXWh2Ud8Ykw5S1i2mqn4RWQRsBNzAr1V1m4h8Pbz8F8AGYBawEzgNzE9VvPHmCjRvoXiPBSDPWi7GmPSQ0hH6qrqBUAKJnPeLiM8KfOMs+ygDyhIQXkJ5glXNpl0n6yy5GGPSho3QT4oW9yUoeIKh7q/K+sMAuKtDd5RV11u3mDGm+7PkklBtPdEFvMHjBDXIsVM7QmvWha7JWMvFGJMOLLkkQ5Qck8UJnt33Wyrf3wiAOxBa6WS9Fa80xnR/VhU5RbK9DVz8rzfJPlVLw5g6PEEX+b58KusqUx2aMcZ0miWXJJKI/7p9Sn5VDQC1gdN48dA7ozeVtZWpCs8YY+LGusVSQFFcvjM1xOr9p8lweejr6s/xuuMpjMwYY+LDkkuKuL1nkkuD/zQZbi/9dIB1ixlj0oIll2SKGKTv9p2ZaKivJsOVQS9/b47XWsvFGNP92TWXFHGFWy7efCVYf5IMVyb5NdlUuipRVUTavo3ZGGO6Omu5JIFEuRfZHb7mklXogboTiAj5JzKpC9RR469JdojGGBNXllySrDHNuPNyAfD1ycBNqEJy7onQUruob4zp7iy5pIinaCD5554mb0QffBmh6y+Zp0LJxW5HNsZ0d5ZcEqhVZ1jEDOk9iIEfriRzcD8yC7IB8NWFfh3WcjHGdHeWXJJEWz5ks6Ak9O7LJXdAH4IawBcIJxe7Y8wY081ZckmVxuSSkUvOwBJqAqfwaehplDbWxRjT3VlySShte1FB+OnNGXn4BhZT5z9FpsuLL5hpLRdjTLdnySUZog1ZyR8Yevfl4S0upq6hmky3j0IZaC0XY0y3Z8klVXoNDr1nFuApKqKh/iRZ7iz6BwstuRhjuj0boZ8qBQPhyt/AsGm4s3qhdSfxuTLoX9uXI7UHUh2dMcZ0irVcUmnk5yC7DyKCS0OPNy6u6mUtF2NMt2fJJanarhfm8/kB6FOdxbHaY8kKyBhjEsKSSxK8fDAQMRU9wWT2DpWDyavJ4HjtcRoCDUmIzBhjEsOSSxLsOhFsNt1qQCWQW3IOANkNPhTlcM3hpMRmjDGJYMmli8gfei61gdNkqQ+AQ6cOpTgiY4zpOEsuXUTGoBJqGk6S7fLhDno4dNqSizGm+7Lk0kV4i4qoqT9BtjuTnLpe1nIxxnRrllwSSKNUfxGJfknf078//rpKcjw59G8ospaLMaZbs+SSDNruJADidqP+E7jFw/DaIRw8dTApoRljTCJYckmitke5hGR4awE4t6qPtVyMMd2aJZcuJLtfaKxL/9psu+ZijOnWHCcXEclJZCBpSc7WVmku/7zQM17y6zI5UnOEhqANpDTGdE9nTS4iMllE3gG2h6fHiMjPEh5ZGtAW72dTcMFwavzV5GhoIGVFTUWiQjPGmIRy0nL5MfApoAJAVd8EPprIoNKG06wSlnHuudT4q8hyZeIKuuyivjGm23LULaaqe1vMCkRd0XSKt7iY2roT5LizyKvry8HTllyMMd2Tk+SyV0QmAyoiPhH5NuEuss4SkctE5F8islNEbo2yXETkp+HlW0VkXHj+IBF5XkS2i8g2EflmPOJJNfF68TccJ9uTS5/ThRystuRijOmenCSXrwPfAAYC5cDY8HSniIgbeBC4HBgBXC0iI1qsdjkwPPxaCPw8PN8PfEtVLwImAd+Ism235NYTuMTF6Kph7Knak+pwjDGmQ876JEpVPQpck4BjTwB2quouABFZC8wB3olYZw7wO1VV4J8i0ktEilT1AHAgHF+ViGwnlPzeoUs7+91j2QWhdUpPFfLXk1sSHI8xxiTGWZOLiAwFbgSGRK6vqrM7eeyBQOS1nHJgooN1BhJOLOH4hgCXAC93Mp4E0IhPjZ/bTzC9LxoG70FhXQ57TljLxRjTPZ01uQBPAKuAp4Bg+6vGJNq3bMv7q9pdR0RygceAf1PVk1EPIrKQUJcahYWFlJWVdSjYjggGQz8ubXEairYZh0f9nOM/Sb5mcbjmMM9seoZMV2aiQ4276urqpP6suwI7556hJ55zRzhJLrWq+tMEHLscGBQxXQLsd7qOiHgJJZaHVfVPbR1EVVcAKwDGjx+v06ZN63TgTj2zbkuby9qKwz9mDFu/9//I8+biDWQweMxgRvTtfpeTysrK2jzHdGXn3DP0xHPuCCcX9O8Xke+LyIdFZFzjKw7HfgUYLiJDRcQHzAXWt1hnPfCV8F1jk4ATqnpARIRQa2q7qt4Xh1i6DE/v3tTVHSPPk0fB6XPYc9K6xowx3Y+Tlsso4MvAxznTLabh6Q5TVb+ILAI2Am7g16q6TUS+Hl7+C2ADMAvYCZwG5oc3nxKO6S0R2RKed5uqbuhMTIkSYxUYCFbidfkYVjWE90+8n4iQjDEmoZwkl/8FDFPV+ngfPJwMNrSY94uIz0qU255V9e84ufWqi3BFtA+dBJ2d6wdg7Mlzef/k+wmJyRhjEslJt9ibQK8Ex5HW3C0zylnKwvQ+fyAAg0/3sW4xY0y35KTlUgi8KyKvAHWNM+NwK3KP4XbF1sjqM+lSDuw6TF/N5f2T76OqSMx9a8YYkzpOksv3Ex5Fuooc2hJDEcus8y/gZO3b9PL1InjKxYFTByjOLU5EhMYYkxBORui/kIxA0pkrxuTizs2hrvYg5+Scy8CT57KtYpslF2NMt9LmNRcR+Xv4vUpETka8qkQk6oBF01xjPnF3oEsrw1uJS1x86Nhwth3dFt/AjDEmwdpsuajqR8LveckLJz01v+TiLNH0vbAQDsDFp4p5uuL1hMRljDGJ4uRJlA85mWfa1jqdnL2PrPATM6jxV3NOMI9tFdsI3ZVtjDHdg5NbkS+OnBARD/ChxISTXjqTDjIvuJCTtUfo7cmj/lSA8uryuMVljDGJ1t41l++KSBUwOvJ6C3AIeDJpEfZQ4vFQX3eQPG8B5x27kG0Vdt3FGNN9tJlcVPWe8PWW5aqaH37lqWpfVf1uEmNMG7Fe1s/JPo5LXHz08MW8U9HFH1VjjDERztotZomk80JJJfZOsoEzLsUfbODCuv68dui1eIdljDEJ4+Sai+kwjfhvy7ln13vm5Ryr2U9/b292HNjNyXq7A9wY0z1YcunCXFlZ1NXupZevNxcduYhXDryS6pCMMcaRdpOLiLhE5O1kBZOuOlMVrKBfNQDTDl/EP/b/Iz4BGWNMgrWbXFQ1CLwpIoOTFE9aaRya0vIxx7EY/PnLqQ/UUtowgJf2vxSnyIwxJrGcdIsVAdtE5DkRWd/4SnRg6aQzLZecD03h6KndFGcW0rDfzd6Te+MWlzHGJIqTqsh3JzyKNCVR0kqsiUZE8AR3kOG+iM9+MImy8jK+POLL8QnQGGMSxMmtyC8A7wPe8OdXACt25UDjXWHalFE61oY5739dSn2glkvqzuXP7/05HqEZY0xCOakt9jVgHfDL8KyBwBMJjCntxFhxv5X8mVdypHoXAzMHcHz3KXYe3xmv0IwxJiGcXHP5BjAFOAmgqjuA/okMyjQnLheZnp34XBl8cffHeGrXU6kOyRhj2uUkudSpan3jRLhwpZXojUE8HlA8/JrpnGo4wYd0KH/513MEgoE47NUYYxLDSXJ5QURuA7JEZCbwKGB/OndYx1JN9sTLOV65mXMyC5m8dTTPffBcnOMyxpj4cZJcbgWOAG8B/wfYANyRyKDSRpvPYOlAw0+ECz+RT32glhmnL2bFll8R1GCnwjPGmERxklzmAL9T1S+o6pWq+iu1J1fFJB7dYgC9vnAjh09spSSrhNLXBvH83ufjtGdjjIkvJ8llNvA/IvKQiHw6fM3FONBeBg4GO5CfPRmcN/oQ/mA9X6iaxMqXf23XXowxXZKTcS7zgVJC11r+N/CeiKxMdGDdnT/QfpfV4aq6Du33nAXf5cixTRRmFvGpv4/l4e0Pd2g/xhiTSI6qIqtqA/A0sBZ4jVBXmWlHnT8iuUjrjrEPjp3u2I6zenPJlwZTWXeYKd5RvLB+E++feL9j+zLGmARxMojyMhH5LbATuBJYSajemGlHs+QSKZxn9nY0uQC+j/0fsngUr8vHokOf4T82/JC6QMdaQsYYkwhOWi7XEhqRf76qzlPVDarqT2hUaaC2oe1rIQrsPd7x5ILLxXm33MGxisc4J7OQL704mds23mHXX4wxXYaTay5zgZeAmSLyGRGx0fkOtNlyCetwt1ijASO55N+ms+/YPzgv9zyueHokP9z4Q/xBy/vGmNRz0i32BWAz8AXgi8DLInJlogPr7ur8Adq7Cbn8WE2njyGjrmT852vZf+yvDM0ZxmXPlLL0wTvsccjGmJRz0i12B3BpuEvsK8AE4M7EhtX91TVEXtBvvkzoZLdYBPcnv8eEL2dy8NhT9M8cwNwPprPh33/OX999IS77N8aYjnCSXFyqejhiusLhdj1aW9dcGvPMwZO14dZNJ4nA1G8x/hvjcNX/gBp/NR/JnUzvXx5j3TeX88bb/+z8MYwxJkZOBkQ+IyIbgUfC01cRKgFj2hH1mks4syihyjDvHT7FiOL8+BxwzFwG3fURip64iTf/MZCc3BlMyphEzepKyk7/iur804y+egaDR4+Mz/GMMaYdZ00uqnqziHwO+Aihr8cVqvp4wiPr5lonl+Yj8r1u4Y+v7uWu2RfH76AFJXjmPcqHpr9M4G/38+aL2fjdUxmaOxy3uvH//jBv1P+Rkw2V1GfUk1vci0ETRlA04WLcmRnxi8OYLq6hoYHy8nJqa2tj3ragoIDt27cnIKrUy8zMpKSkBK/X2+l9OSrloqp/Av7U6aO1ICKXAfcDbmClqi5tsVzCy2cBp4FrVfV1J9umWltdXkooQ396VBHrXivn25+6gNyMOFfUGTwR9zVrGDf7ILz7/6h542e89nZfgoFRFPj6c27OcDwuDxwFNpzi/af+Sk3gFLWBGuqDtfi1gaAriMvnwpubiS8/m8zeOWSd04vcgX3JHVxEZu8cXC7rHTXdU3l5OXl5eQwZMgSJMsi5PVVVVeTl5SUostRRVSoqKigvL2fo0KGd3l/K6oSJiBt4EJgJlAOviMh6VX0nYrXLgeHh10Tg58BEh9umVG1DkMx2ll87ZShPbNnPnU+8zbLPj8bnScAXdd4AuHQBWZcu4CPBIBz9Fxzaxsl9L7Bl6yGqDuXhrutHFvlkuXLJcmeT7y0gw5WF2xX+p1Ebfh0G/uWnnkMc4xCqil/r8Qf9+LWBgPoJBAOhdwJkBgO89KffEiSAiqKuYOhKnUsQj+DyCOL14PZ5cGd4cft8eDIz8WRl4s7y4c4MvTzZGXgyM3Bn+HBnefFkeHF7Pbh9bjxeDy6fB4SYvyBMz1ZbW9uhxJLORIS+ffty5MiRuOwvlUUoJwA7VXUXgIisJVRWJjJBNFZkVuCfItJLRIqAIQ62Tak6f6Dd5DJ2UC/+feb53Pff/8MbHxznUxcPYNg5ORRkefG6Q4kmEFQCQaU+EMQfUPzBIPUBJRAI4g8q/vDyZi9VgqqohopjBpXwtKJAUC9A9QJ0OGgpgKKBIJmBE2TWH0WDuwg0HMRVU0nWCT/Zpzxk1Pnw+TPwBTPwBjPw4sOLF4/48IgHt8sTehc3bvHgc2Xg9pyZdosHjyuima1AQ/gV9aa5AFAD1DRbtT2h8wuG38Ofw/NBm/4v9P+N88Lnr02fzrxrqzlNwTevCX5mnT7A1icfa7b2meO13qLVOTTbruU20bdoe3nbhVFb1zSPrYhq5NrZwOY/taxvl5yi6W39HJ3IuOpCTu091qFtBag+XtHhY7etaxSbrztWzcs3/q7T+3GUXEQkCxisqv/q9BHPGAjsjZguJ9Q6Ods6Ax1uC4CILAQWAhQWFlJWVtapoJ3a9n4Dk8OfNeK5K4KgQFlZGaPdcOMlGTy7p5Zf/W0XHSmUHEkAl4RuIHPJmWki38MrSmM04c8ChH7NwxCGndlpZugl0uqO6tAeBJQg0IBHasmQE/hc1bipxSu1ZOhpMhpOk9XgJ88fJNuvZAVc+IJCZtCDN+jCo27c6sKjLjx4EFy41Y0bNy5cCBJ+d4Xe5cw8EFwSmgqdjyAq4b9IJeI8JTwlZ+ZK6JM2LhEBlfDPp3FLzuynxV+50vQTOfOp8UPUZZErtPgZNouLNn7YzbaWaDOjHrF9Z4suVu3vIX4Nhc7tKIALt3Slbt2u04Jy4SLXW9Dp/Zw1uYjIZ4F7AR8wVETGAktUdXYnjx3tp9ny67WtdZxsG5qpugJYATB+/HidNm1aDCF23Pay9+DlbQCIq+WPWWmMYxrwLaAhEOTgiVqq6/zUh28GcLsEt0vwugWv24XH7cLrEjxuFx634Akvd0vovSs18cvKykjWz7qrsHPuPrZv307WoN4d2jYe11wqKytZs2YNN9xwg6P1d+/ezdy5czl27Bjjxo3joYcewufztVpv9erV/OAHPwDgjjvuYN68eTHH5q3O5qL7WtQm/nHMu3E0XuUuQl1YlQCquoVQt1RnlQODIqZLgP0O13GybUq1V1ssGq/bxaA+2VxUlM+YQb0YM6gXIwcWcFFRPqX98zi3bw4De2XRPz+TPjk+8jO9ZPs8ZHjceNyuLpVYjDHtq6ys5Gc/+5nj9W+55RZuuukmduzYQe/evVm1alWrdY4dO8bdd9/Nyy+/zObNm7n77rs5fvx4PMOOiZNuMb+qnkjAl9crwHARGQrsA+YSel5MpPXAovA1lYnACVU9ICJHHGybUnX+YETHiDGmq7r7qW28s995yaRAIIDb7W53nRHF+Xz/s20PM7j11lt57733GDt2LDNnzmT58uVtrquqbNq0iTVr1gAwb9487rrrLq6//vpm623cuJGZM2fSp08fAGbOnMkzzzzD1Vdf7fTU4spJcnlbRP434BaR4cBi4B+dPbCq+kVkEbCR0O3Ev1bVbSLy9fDyXxAarDmLULn/08D89rbtbEzxFO1W5Nj6wo0x6Wrp0qW8/fbbbNmyhaqqKsaOHRt1vTVr1tC/f3969eqFxxP6ui4pKWHfvn2t1t23bx+DBp3p0GlrvWRxklxuBG4H6giN0t8I/N94HFxVN9BitH84qTR+VuAbTrftSloOouwa94EYY1pqr4URTbzHueTl5bFly5Y2l0e7NThaT5K2vg0wpd3lTkbonyaUXG5PfDjpI7JwpbVXjDFtqaqqYurUqVGXrVmzhosuuojKykr8fj8ej4fy8nKKi4tbrVtSUtLsbtjy8vKU3mzh5G6x54nyh7eqfjwhEaWJWn8gjrddGmPSSV5eHlVVVU2f22u5AEyfPp1169Yxd+5cVq9ezZw5rZ80/6lPfYrbbrut6SL+X/7yF+655564x+6Uk7vFvg3cHH7dCWwBXk1gTGkhsuWiLmu9GGPO6Nu3L1OmTGHkyJHcfPPNZ11/2bJl3HfffZSWllJRUcGCBQsAePXVV/nqV78KQJ8+fbjzzju59NJLufTSS/ne977XdHE/FZx0i73WYtaLImIPCzmLuJTTN8akrca7v5wYNmwYmzdvbjV//PjxrFy5smn6uuuu47rrrotLfJ3lpFssMvW5gA8BAxIWUZqo8wfbaK5YG8YYk/6c3C32GmdGxfuB3cCCRAaVDuoaAjbOxRjTYznpFut87WUDWJIxxvQcbSaX8APC2hR+xotpw5OLPsKzN74enpJOVXA1xpjupr2Wy2fbWaYk4OFh6UqtyWKM6WHaTC6qOj+ZgRhjjEkfTp/n8mngYjjz/CtVXZKooNKNtPHZGNMzxVpyf+rUqU2DLg8fPsyECRN44oknWq3ndrsZNWoUAIMHD2b9+vVxizlWTm5F/gWhB85NB1YCVwKtb7g2MbDrL8b0ZI0l950ml7/97W9Nnz//+c9HHaEPkJWVddbR/snipOUyWVVHi8hWVb1bRH6EXW9xxFKIMd3A07fCwbccr54V8IP7LF+dA0bB5UvbXBxLyf1IVVVVbNq0id/85jeO400VJ8mlJvx+WkSKgQrAbk+OgXWFGWMixVJyf8SIEU3Tjz/+ODNmzCA/Pz/q+rW1tYwfPx6Px8Ott97KFVdckYDonXGSXP4sIr2A5cDrhP4g/1Uig0pflmaM6XLaaWFEU5PkkvuRHnnkkaZaYtF88MEHFBcXs2vXLj7+8Y8zatQozjvvvDhFGhsngygbn93ymIj8GchU1ROJDcsYY3qGs5Xcb2y5VFRUsHnzZh5//PE299VYin/YsGFMmzaNN954o+smFxF5E/gD8AdVfY/QQ8NMDLRF7X27FmNMzxZryX2ARx99lM985jNkZmZGXX78+HGys7PJyMjg6NGjvPjii3znO9+JZ9gxcVJyfzahmmJ/FJFXROTbIjI4wXEZY0zairXkPsDatWu5+uqrm82LLLm/fft2xo8fz5gxY5g+fTq33nprs+s1yeakW2wP8J/Af4rIcELPdFlG6Nn1xgG70mKMaSmWkvtAs6dMNoosuT958mTeesv5XW+J5nQQ5RDgi8BVQABIXVurG1JLL8aYHsbJNZeXAS/wR+ALqror4VEZY4zp1py0XOap6rsJjySNSTtTxhiTjs56Qd8SSxy0zCd2u5gxJs05uVvMdJIIllCMMT2KJRdjjDFxd9bkIiJuEZktIotF5N8bX8kIzhhj0lFjVWSnHnjgAUpLSxERjh492jS/rKyMgoICxo4dy9ixY1myJPqTUHbv3s3EiRMZPnw4V111FfX19Z0+h7Nx0nJ5CrgW6AvkRbyMQ82e52LX843p8WJNLlOmTOHZZ5/l3HPPbbVs6tSpbNmyhS1btvC9730v6va33HILN910Ezt27KB3796sWrWqw7E75eRusRJVHZ3wSHoUuwBjTFexbPMy3j3m/L6lQCCA293+GPIL+1zILRNuaXN5rCX3L7nkEsfxtaSqbNq0qWnQ5rx587jrrru4/vrrO7xPJ5wkl6dF5JOq+peERmKMMT1ER0vuR/PSSy8xZswYiouLuffee7n44oubLa+oqKBXr154PKGv+5KSEvbt2xeX82iPk+TyT+BxEXEBDYR6eVRVoz9QwDhgfWPGdBXttTCiqUphyf2Wxo0bx549e8jNzWXDhg1cccUV7Nixo9k6qq17SiQJ/fNOksuPgA8Db2m0KE2bxJKIMeYsnJbcjybyoWGzZs3ihhtu4OjRo/Tr169pfr9+/aisrMTv9+PxeCgvL28qzZ9ITpLLDuBtSyydIII2u85iP0pjerKOlNyP5uDBgxQWFiIibN68mWAwSN++fZutIyJMnz6ddevWMXfuXFavXs2cOXM6ewpn5eRusQNAmYh8125FNsaYzou15P5Pf/pTSkpKKC8vZ/To0U1l9tetW8fIkSMZM2YMixcvZu3atU1dXrNmzWL//v0ALFu2jPvuu4/S0lIqKipYsGBB4k4uzEnLZXf45Qu/TCdZu8UYE0vJ/cWLF7N48eJW8xctWsSiRYuibrNhw4amz8OGDWPz5s2xB9kJTp7ncne8DyoifQg93XII8D7wRVU9HmW9y4D7CT07ZqWqLg3PXw58FqgH3gPmq2plvOPsvFAasbEtxpiexskI/XNEZLmIbBCRTY2vTh73VuA5VR0OPBeebnlcN/AgcDkwArhaRBqvbP03MDI8/uZ/gO92Mp6Esue5GGN6GifXXB4G3gWGAncTamm80snjzgFWhz+vBq6Iss4EYKeq7lLVemBteDtU9S+q6g+v90+gpJPxJFTL1GKpxhiT7pwkl76qugpoUNUXVPU6YFInj1uoqgcAwu/9o6wzENgbMV0entfSdcDTnYwniSy1GGPSn5ML+g3h9wMi8mlgPw5aCiLyLDAgyqLbHcYW7Vu42bVwEbkd8BNqXbUVx0JgIUBhYWHU51AnSjAYCtcfCDWyGse9KJrUOFKhuro67c+xJTvn7qOgoKDpVuBYBQKBDm/bHdTW1sbld+okufxARAqAbwH/BeQDN51tI1X9RFvLROSQiBSp6gERKQIOR1mtHBgUMV1CKLE17mMe8BlgRntjcFR1BbACYPz48Tpt2rSzhR43zz62FaCp7EKkZMaRCmVlZWl/ji3ZOXcf27dv7/Ao+3iP0O9qMjMzO1XLrJGTJ1H+WVVPqOrbqjpdVT+kqus7edz1wLzw53nAk1HWeQUYLiJDRcQHzA1v13gX2S3AbFU93clYUsC6xozpyWKtinzNNddwwQUXMHLkSK677joaGkIdSg8//DCjR49m9OjRTJ48mTfffDPq9tdeey1Dhw5tKs3f0UGbsWgzuYjIf4nIT9t6dfK4S4GZIrIDmBmeRkSKRWQDQPiC/SJgI7Ad+KOqbgtv/wChsv//LSJbROQXnYzHGGOSpiPJ5d133+Wtt96ipqaGlStXAjB06FBeeOEFtm7dyp133snChQvb3Mfy5cubSvO3VSgzntrrFns14vPdwPfjdVBVrQBmRJm/H5gVMb0B2BBlvdJ4xZIMrdspNozSmK7i4A9/SN125yX3/YEAx85Scj/jogsZcNttbS6PteT+rFlNX4tMmDCB8vJyACZPntw0f9KkSU3zu4I2k4uqNt4qjIj8W+S0iU3kOBfrEDPGdLTkfkNDAw899BD3339/q3VXrVrF5Zdf3uYxb7/9dpYsWcKMGTNYunQpGRkZnT6P9ji5oA/2p3YHScR/jTFdUXstjGhSWXL/hhtu4KMf/WirKsrPP/88q1at4u9//3vU7e655x4GDBhAfX09CxcuZNmyZW0+tTJenCYX0wmWmY0xbXFacv/uu+/myJEj/PKXv2y2ztatW/nqV7/K008/3aoicqOioiIAMjIymD9/Pvfee28czyC6NpOLiFRx5nsxW0RONi7CHhZmjDEdFmvJ/ZUrV7Jx40aee+45XK4z92F98MEHfO5zn+Ohhx7i/PPPb3P7AwcOUFRUhKryxBNPMHLkyLicR3vavFtMVfNUNT/88kR8zrPEEhvrFjPGRIq15P7Xv/51Dh06xIc//GHGjh3LkiVLAFiyZAkVFRXccMMNjB07lvHjxzdtE1ly/5prrmHUqFGMGjWKo0ePcscddyTmxCJYt1giRe0Ps1RjjImt5L7f7486f+XKlU23JbcUWXJ/06bO1hqOnZPaYqaTXJZQjDE9jCWXZLDcYozpYSy5GGOMiTtLLklgDRdjTE9jycUYY0zcWXJJAbVhlcaYNGfJJSnEEooxpkmsVZEXLFjAmDFjGD16NFdeeSXV1dUAqCqLFy+mtLSU0aNH8/rrr0fdfvfu3UycOJHhw4dz1VVXUV9fH5fzaI8ll5SwqzDG9GSxJpcf//jHvPnmm2zdupXBgwfzwAMPAPD000+zY8cOduzYwYoVK7j++uujbn/LLbdw0003sWPHDnr37s2qVavich7tsUGUyWC5xJgu629//B+O7q12vH4gEMB9lpL7/QblMvWLbZdjibXkfn5+qCiKqlJTU4NI6EvlySef5Ctf+QoiwqRJk6isrGwq9dJIVdm0aVPToM158+Zx1113tZmI4sWSizHGJFlHSu7Pnz+fDRs2MGLECH70ox8BsG/fPgYNOvM0+JKSEvbt29csuVRUVNCrV6+mx603rpNollwSKUqLxRoxxnQt7bUwoklVyf3f/OY3BAIBbrzxRv7whz8wf/58VFtfy21s1TRysk4iWHJJipa/XLu4b4wJcVpyH8DtdnPVVVexfPly5s+fT0lJCXv37m1aXl5eTnFxcbN99OvXj8rKSvx+Px6PJ+o6iWDJJQmstWKMiRRLyX1V5b333qO0tBRV5amnnuLCCy8EYPbs2TzwwAPMnTuXl19+mYKCgmZdYhBqpUyfPp1169Yxd+5cVq9ezZw5cxJ2bo3sbrEksiRjjIHYSu6rKvPmzWsqmX/gwIGmp0jOmjWLYcOGUVpayte+9rVmd6BFltxftmwZ9913H6WlpVRUVLBgwYLEnVyYtVySwVK4MaYFpyX3XS4XL774YtRlIsKDDz4YdVlkyf1hw4axefPm2IPsBPvaSyrB2i/GmJ7Akosxxpi4s+SSFNZaMcb0LJZcksKSizGmZ7HkkgRJGK9kjDFdiiUXY4wxcWfJJQmkzQljTE8Ur5L7y5cvZ+zYsYwdO5aRI0fidrs5duxYq+2vvfZahg4d2rSuk3IznWXJJWm06c2KvxjTs8Wr5P7NN9/Mli1b2LJlC/fccw8f+9jH6NOnT9R9LF++vGndtgplxpMNokwBa7wY03U8/9sVHN6zy/H6AX8At6f9kvv9zx3G9GsXtrk8XiX3Iz3yyCNcffXVDs4gOSy5GGNMksWr5H6j06dP88wzzzS1aKK5/fbbWbJkCTNmzGDp0qVkZGTE7XyiseSSDNZUMabLaq+FEU1XKrnf6KmnnmLKlCltdondc889DBgwgPr6ehYuXMiyZcua6pMlil1zSbJkPEfBGNN9NLZcor3eeeedZus2ltx/7LHHms1fu3Ztu11iRUVFiAgZGRnMnz8/KXXGrOWSdIpd0jemZ4tXyX2AEydO8MILL/D73/++zX00PvpYVXniiScYOXJk3M6lLdZyMcaYJItXyX2Axx9/nE9+8pPk5OQ02y6y5P4111zTtP3Ro0e544474n9SLaSk5SIifYA/AEOA94EvqurxKOtdBtwPuIGVqrq0xfJvA8uBc1T1aILD7jBp9cEY09PFo+Q+hMawXHvtta3mR5bc37RpU8zxdVaqWi63As+p6nDgufB0MyLiBh4ELgdGAFeLyIiI5YOAmcAHSYm4E1SsI8wY07OkKrnMAVaHP68GroiyzgRgp6ruUtV6YG14u0Y/Br6DfW8bY0yXk6oL+oWqegBAVQ+ISP8o6wwE9kZMlwMTAURkNrBPVd88291XIrIQWAhQWFhIWVlZ56N3SIOhvOdvCERGhEJS40iF6urqtD/Hluycu4+CgoKmC+qxCgQCHd62O6itrY3L7zRhyUVEngUGRFl0u9NdRJmnIpId3scnnexEVVcAKwDGjx+v06ZNc3j4znv2sbcA8Hlb/5iTGUcqlJWVpf05tmTn3H1s3769w2NV4j3OpavJzMzkkksu6fR+EpZcVPUTbS0TkUMiUhRutRQBh6OsVg4MipguAfYD5wFDgcZWSwnwuohMUNWDcTsBY4wxHZaqay7rgXnhz/OAJ6Os8wowXESGiogPmAusV9W3VLW/qg5R1SGEktA4SyzGGNN1pCq5LAVmisgOQnd8LQUQkWIR2QCgqn5gEbAR2A78UVW3pSheY4yJm1irIje68cYbyc3NbZouKyujoKCgaUT/kiVLom63e/duJk6cyPDhw7nqqquor6/vcOxOpSS5qGqFqs5Q1eHh92Ph+ftVdVbEehtU9XxVPU9V/6ONfQ3pymNcgKarR2IDXYwxdCy5vPrqq1RWVraaP3Xq1KZS+m3VC7vlllu46aab2LFjB71792bVqlUdCTsmVv4lWbTNCWNMClU+9R71+085Xj8Q8FPjbv+r01ecQ6/Pntfm8lhL7gcCAW6++WbWrFnD448/7jhWCI3w37RpU9OgzXnz5nHXXXdx/fXXx7SfWFlySShrqRhjWou15P4DDzzA7NmzKSoqarXOSy+9xJgxYyguLubee+/l4osvbra8oqKCXr164fGEvu5LSkrYt29f3M+pJUsuSWFJxpiuqr0WRjTJLrm/f/9+Hn300ahjT8aNG8eePXvIzc1lw4YNXHHFFezYsaPZOqqte0qSUZ3dkosxxqRQVVUVU6dOjbpszZo17N69m507d1JaWgqEHgxWWlrKzp07m55QCaFClTfccANHjx6lX79+TfP79etHZWUlfr8fj8dDeXk5xcXFiT0pLLkkhz3DxRgTIZaS+yNGjODgwTMjLXJzc9m5cycABw8epLCwEBFh8+bNBINB+vbt22x7EWH69OmsW7eOuXPnsnr1aubMmUOiWcn9JLM7xowxsZTcb8+6desYOXIkY8aMYfHixaxdu7apyyuy5P6yZcu47777KC0tpaKiggULFsTlPNpjLRdjjEkBpyX3W6qurm76vGjRIhYtWhR1vciS+8OGDUvK0ycjWcslCaytYozpaSy5JI2NbTHG9ByWXIwxPVK0W3R7unj+TCy5GGN6nMzMTCoqKizBRFBVKioqyMzMjMv+7IJ+Ctg/Z2NSq6SkhPLyco4cORLztrW1tXH7Au5qMjMzKSkpicu+LLkkgV3QN6Zr8Xq9DB06tEPblpWVxeVhWunOusWSwJKLMaanseRijDEm7iy5JIM1XYwxPYz0pLslROQIsCcFh+4HdO0HmsWfnXPPYOfcM1ygqjGVgu5RF/RV9ZxUHFdEXlXV8ak4dqrYOfcMds49g4i8Gus21i1mjDEm7iy5GGOMiTtLLsmxItUBpICdc89g59wzxHzOPeqCvjHGmOSwlosxxpi4s+RijDEm7iy5JJiIXCYi/xKRnSJya6rjSTQRGSQiz4vIdhHZJiLfTHVMySAibhF5Q0T+nOpYkkVEeonIOhF5N/z7/nCqY0okEbkp/G/6bRF5RETSrnqliPxaRA6LyNsR8/qIyH+LyI7we28n+7LkkkAi4gYeBC4HRgBXi8iI1EaVcH7gW6p6ETAJ+EYPOGeAbwLbUx1Ekt0PPKOqFwJjSOPzF5GBwGJgvKqOBNzA3NRGlRC/BS5rMe9W4DlVHQ48F54+K0suiTUB2Kmqu1S1HlgLzElxTAmlqgdU9fXw5ypCXzgDUxtVYolICfBpYGWqY0kWEckHPgqsAlDVelWtTGlQiecBskTEA2QD+1McT9yp6l+BYy1mzwFWhz+vBq5wsi9LLok1ENgbMV1Omn/RRhKRIcAlwMspDiXRfgJ8BwimOI5kGgYcAX4T7g5cKSI5qQ4qUVR1H3Av8AFwADihqn9JbVRJU6iqByD0xyPQ38lGllwSK1rJyh5x77eI5AKPAf+mqidTHU+iiMhngMOq+lqqY0kyDzAO+LmqXgKcwmF3SXcUvs4wBxgKFAM5IvKl1EbVtVlySaxyYFDEdAlp2JRuSUS8hBLLw6r6p1THk2BTgNki8j6hbs+Pi8jvUxtSUpQD5ara2CpdRyjZpKtPALtV9YiqNgB/AianOKZkOSQiRQDh98NONrLkklivAMNFZKiI+AhdAFyf4pgSSkSEUD/8dlW9L9XxJJqqfldVS1R1CKHf7yZVTfu/aFX1ILBXRC4Iz5oBvJPCkBLtA2CSiGSH/43PII1vYGhhPTAv/Hke8KSTjXpUVeRkU1W/iCwCNhK6u+TXqrotxWEl2hTgy8BbIrIlPO82Vd2QupBMgtwIPBz+w2kXMD/F8SSMqr4sIuuA1wndEfkGaVgGRkQeAaYB/USkHPg+sBT4o4gsIJRkv+BoX1b+xRhjTLxZt5gxxpi4s+RijDEm7iy5GGOMiTtLLsYYY+LOkosxxpi4s+RijDEm7iy5GNOFSIj979J0e/aP2JgUE5Eh4eeh/IzQIL1BZ9vGmK7OBlEak2Lh6tG7gMmq+s8Uh2NMXFjLxZiuYY8lFpNOLLkY0zWcSnUAxsSTJRdjjDFxZ8nFGGNM3NkFfWOMMXFnLRdjjDFxZ8nFGGNM3FlyMcYYE3eWXIwxxsSdJRdjjDFxZ8nFGGNM3FlyMcYYE3f/H2F/7plnUjmsAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "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": []
+  }
+ ],
+ "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.9.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/source/bhinitialconditionswithSF.py b/source/bhinitialconditionswithSF.py
new file mode 100644
index 0000000..f98ee22
--- /dev/null
+++ b/source/bhinitialconditionswithSF.py
@@ -0,0 +1,66 @@
+# bhinitialconditions.py
+
+# set the initial conditions for all the variables for an isotropic Schwarzschild BH
+# see further details in https://github.com/GRChombo/engrenage/wiki/Running-the-black-hole-example
+
+from source.uservariables import *
+from source.tensoralgebra import *
+from source.Grid import *
+import numpy as np
+from scipy.interpolate import interp1d
+
+def get_initial_state(a_grid) :
+    
+    # For readability
+    r = a_grid.r_vector
+    N = a_grid.num_points_r
+                     
+    initial_state = np.zeros(NUM_VARS * N)
+    [u,v,phi,hrr,htt,hpp,K,arr,att,app,lambdar,shiftr,br,lapse] = np.array_split(initial_state, NUM_VARS)
+
+    # Set BH length scale
+    GM = 1.0
+    
+    # set non zero metric values
+    grr = (1 + 0.5 * GM/r)**4.0
+    gtt_over_r2 = grr
+    gpp_over_r2sintheta = gtt_over_r2
+    phys_gamma_over_r4sin2theta = grr * gtt_over_r2 * gpp_over_r2sintheta
+    # Note sign error in Baumgarte eqn (2), conformal factor
+    phi[:] = 1.0/12.0 * np.log(phys_gamma_over_r4sin2theta)
+    # Cap the phi value in the centre to stop unphysically large numbers at singularity
+    phi[:] = np.clip(phi, None, 10.0)
+    em4phi = np.exp(-4.0*phi)
+    hrr[:] = em4phi * grr - 1.0
+    htt[:] = em4phi * gtt_over_r2 - 1.0
+    hpp[:] = em4phi * gpp_over_r2sintheta - 1.0    
+    lapse.fill(1.0)
+    #lapse[:] = em4phi # optional, to pre collapse the lapse
+    
+    # overwrite inner cells using parity under r -> - r
+    a_grid.fill_inner_boundary(initial_state)
+             
+    dhrrdx     = np.dot(a_grid.derivatives.d1_matrix, hrr)
+    dhttdx     = np.dot(a_grid.derivatives.d1_matrix, htt)
+    dhppdx     = np.dot(a_grid.derivatives.d1_matrix, hpp)
+
+    # assign lambdar values
+    h_tensor = np.array([hrr, htt, hpp])
+    a_tensor = np.array([arr, att, app])
+    dhdr   = np.array([dhrrdx, dhttdx, dhppdx])
+        
+    # (unscaled) \bar\gamma_ij and \bar\gamma^ij
+    bar_gamma_LL = get_metric(r, h_tensor)
+    bar_gamma_UU = get_inverse_metric(r, h_tensor)
+        
+    # The connections Delta^i, Delta^i_jk and Delta_ijk
+    Delta_U, Delta_ULL, Delta_LLL  = get_connection(r, bar_gamma_UU, bar_gamma_LL, h_tensor, dhdr)
+    lambdar[:]   = Delta_U[i_r]
+
+    # Fill boundary cells for lambdar
+    a_grid.fill_outer_boundary_ivar(initial_state, idx_lambdar)
+
+    # overwrite inner cells using parity under r -> - r
+    a_grid.fill_inner_boundary_ivar(initial_state, idx_lambdar)
+            
+    return initial_state
diff --git a/source/fluxdiagnostic.py b/source/fluxdiagnostic.py
new file mode 100644
index 0000000..37a8d4d
--- /dev/null
+++ b/source/fluxdiagnostic.py
@@ -0,0 +1,126 @@
+#hamdiagnostics.py
+
+# python modules
+import numpy as np
+import time
+
+# homemade code
+from source.uservariables import *
+from source.Grid import *
+from source.tensoralgebra import *
+from source.mymatter import *
+
+# The diagnostic function returns the Hamiltonian constraint over the grid
+# it takes in the solution of the evolution, which is the state vector at every
+# time step, and returns the spatial profile Ham(r) at each time step
+def get_Ham_diagnostic(solutions_over_time, t, my_grid) :
+
+    start = time.time()
+    
+    # For readability
+    r = my_grid.r_vector
+    N = my_grid.num_points_r
+    
+    Ham = []
+    num_times = int(np.size(solutions_over_time) / (NUM_VARS * N))
+    
+    # unpack the vectors at each time
+    for i in range(num_times) :
+        
+        t_i = t[i]
+        
+        if(num_times == 1):
+            solution = solutions_over_time
+        else :
+            solution = solutions_over_time[i]
+
+        # Assign the variables to parts of the solution
+        (u, v , phi, hrr, htt, hpp, K, 
+         arr, att, app, lambdar, shiftr, br, lapse) = np.array_split(solution, NUM_VARS)
+        
+        ################################################################################################
+        
+        # first derivatives
+        dudx       = np.dot(my_grid.derivatives.d1_matrix, u      )
+        dvdx       = np.dot(my_grid.derivatives.d1_matrix, v      )
+        dphidx     = np.dot(my_grid.derivatives.d1_matrix, phi    )
+        dhrrdx     = np.dot(my_grid.derivatives.d1_matrix, hrr    )
+        dhttdx     = np.dot(my_grid.derivatives.d1_matrix, htt    )
+        dhppdx     = np.dot(my_grid.derivatives.d1_matrix, hpp    )
+        darrdx     = np.dot(my_grid.derivatives.d1_matrix, arr    )
+        dattdx     = np.dot(my_grid.derivatives.d1_matrix, att    )
+        dappdx     = np.dot(my_grid.derivatives.d1_matrix, app    )
+        dKdx       = np.dot(my_grid.derivatives.d1_matrix, K      )
+        dlambdardx = np.dot(my_grid.derivatives.d1_matrix, lambdar)            
+
+        # second derivatives
+        d2udx2      = np.dot(my_grid.derivatives.d2_matrix, u     )
+        d2phidx2    = np.dot(my_grid.derivatives.d2_matrix, phi   )
+        d2hrrdx2    = np.dot(my_grid.derivatives.d2_matrix, hrr   )
+        d2httdx2    = np.dot(my_grid.derivatives.d2_matrix, htt   )
+        d2hppdx2    = np.dot(my_grid.derivatives.d2_matrix, hpp   )      
+            
+        ##############################################################################################
+
+        h = np.array([hrr, htt, hpp])
+        a = np.array([arr, att, app])
+        em4phi = np.exp(-4.0*phi)
+        dhdr   = np.array([dhrrdx, dhttdx, dhppdx])
+        d2hdr2 = np.array([d2hrrdx2, d2httdx2, d2hppdx2])
+              
+        # Calculate some useful quantities
+        ########################################################
+        
+        # \hat \Gamma^i_jk
+        flat_chris = get_flat_spherical_chris(r)
+        
+        # rescaled \bar\gamma_ij
+        r_gamma_LL = get_rescaled_metric(h)
+        r_gamma_UU = get_rescaled_inverse_metric(h)
+        
+        # (unscaled) \bar\gamma_ij and \bar\gamma^ij
+        bar_gamma_LL = get_metric(r, h)
+        bar_gamma_UU = get_inverse_metric(r, h)
+        
+        # \bar A_ij, \bar A^ij and the trace A_i^i, then Asquared = \bar A_ij \bar A^ij
+        bar_A_LL = get_A_LL(r, a)
+        bar_A_UU = get_A_UU(a, r_gamma_UU, r)
+        traceA   = get_trace_A(a, r_gamma_UU)
+        Asquared = get_Asquared(a, r_gamma_UU)
+        
+        # The connections Delta^i, Delta^i_jk and Delta_ijk
+        Delta_U, Delta_ULL, Delta_LLL  = get_connection(r, bar_gamma_UU, bar_gamma_LL, h, dhdr)
+        bar_Rij = get_ricci_tensor(r, h, dhdr, d2hdr2, lambdar, dlambdardx, 
+                                       Delta_U, Delta_ULL, Delta_LLL, bar_gamma_UU, bar_gamma_LL)
+        bar_Rij_diag = np.array([bar_Rij[i_r][i_r],bar_Rij[i_t][i_t],bar_Rij[i_p][i_p]])
+        bar_R   = get_trace(bar_Rij_diag, bar_gamma_UU)
+        
+        # Matter sources
+        matter_rho            = get_rho(u, dudx, v, bar_gamma_UU, em4phi )
+
+        # End of: Calculate some useful quantities, now start diagnostic
+        #################################################################
+
+        # Get the Ham constraint eqn (13) of Baumgarte https://arxiv.org/abs/1211.6632
+        Ham_i = (  two_thirds * K * K - Asquared
+                         + em4phi * ( bar_R
+                                      - 8.0 * bar_gamma_UU[i_r][:] * (dphidx * dphidx 
+                                                                        + d2phidx2)
+               # These terms come from \bar\Gamma^r d_r \phi from the \bar D^2 \phi term
+                                      + 8.0 * bar_gamma_UU[i_t][:] 
+                                            * flat_chris[i_r][i_t][i_t] * dphidx
+                                      + 8.0 * bar_gamma_UU[i_p][:] 
+                                            * flat_chris[i_r][i_p][i_p] * dphidx
+                                      + 8.0 * Delta_U[i_r] * dphidx)
+                         - 2.0 * eight_pi_G * matter_rho ) 
+        
+        # Add the Ham value to the output
+        Ham.append(Ham_i)
+        
+    # end of iteration over time  
+    #########################################################################
+    
+    end = time.time()
+    #print("time at t= ", t_i, " is, ", end-start)
+    
+    return Ham