{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d6633e83",
   "metadata": {},
   "source": [
    "# Buchdahl equation of state example for O$_2$sclpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a9e763c",
   "metadata": {},
   "source": [
    "See the O$_2$sclpy documentation at https://awsteiner.org/code/o2sclpy for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4e1c97bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import o2sclpy\n",
    "import matplotlib.pyplot as plot\n",
    "import ctypes\n",
    "import numpy\n",
    "import sys\n",
    "\n",
    "plots=True\n",
    "if 'pytest' in sys.modules:\n",
    "    plots=False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "366d67d6",
   "metadata": {},
   "source": [
    "Link the O$_2$scl library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dae07403",
   "metadata": {},
   "outputs": [],
   "source": [
    "link=o2sclpy.linker()\n",
    "link.link_o2scl()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3461505",
   "metadata": {},
   "source": [
    "Get a copy (a pointer to) the O$_2$scl unit conversion object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6da48214",
   "metadata": {},
   "outputs": [],
   "source": [
    "cu=link.o2scl_settings.get_convert_units()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "227d3d42",
   "metadata": {},
   "source": [
    "Create the Buchdahl EOS object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6ca1dcae",
   "metadata": {},
   "outputs": [],
   "source": [
    "b=o2sclpy.eos_tov_buchdahl(link)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c659e14",
   "metadata": {},
   "source": [
    "Create the TOV solve object, set the EOS and compute the M-R curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "44f497e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.400000e+00 solar mass mode.\n",
      "Central P: 5.2000e-05 (Msun/km^3), M: 2.2703e+00 (Msun), R: 1.5779e+01 (km)\n",
      "Central P: 5.2000e-05 (Msun/km^3), M: 2.2703e+00 (Msun), R: 1.5779e+01 (km)\n",
      "Central P: 5.2000e-05 (Msun/km^3), M: 2.2703e+00 (Msun), R: 1.5779e+01 (km)\n",
      "Central P: 1.5746e-05 (Msun/km^3), M: 1.2137e+00 (Msun), R: 1.5329e+01 (km)\n",
      "Central P: 2.2140e-05 (Msun/km^3), M: 1.4453e+00 (Msun), R: 1.5395e+01 (km)\n",
      "Central P: 2.0889e-05 (Msun/km^3), M: 1.4027e+00 (Msun), R: 1.5381e+01 (km)\n",
      "Central P: 2.0811e-05 (Msun/km^3), M: 1.4000e+00 (Msun), R: 1.5380e+01 (km)\n",
      "Central P: 2.0812e-05 (Msun/km^3), M: 1.4000e+00 (Msun), R: 1.5380e+01 (km)\n",
      "Central P: 2.0812e-05 (Msun/km^3), M: 1.4000e+00 (Msun), R: 1.5380e+01 (km)\n",
      "Gravitational mass is: 1.400000e+00\n",
      "Radius is: 1.538043e+01\n",
      "Exact radius is 1.538043e+01, computed radius is 1.538043e+01.\n",
      "Relative difference 5.976744e-10.\n"
     ]
    }
   ],
   "source": [
    "ts=o2sclpy.tov_solve(link)\n",
    "ts.set_eos(b);\n",
    "ts.fixed(1.4,1.0e-4)\n",
    "print('Exact radius is %7.6e, computed radius is %7.6e.' %\n",
    "      (b.rad_from_gm(1.4),ts.rad))\n",
    "print('Relative difference %7.6e.' %\n",
    "      (abs(b.rad_from_gm(1.4)-ts.rad)/ts.rad))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bfcc1ee",
   "metadata": {},
   "source": [
    "Get the table for the TOV results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66d67093",
   "metadata": {},
   "outputs": [],
   "source": [
    "tov_table=ts.get_results()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3774bba0",
   "metadata": {},
   "source": [
    "The compactness of a 1.4 solar mass NS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8d34b357",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta=ts.mass*b.G_km_Msun/ts.rad"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d08055ac",
   "metadata": {},
   "source": [
    "Construct two lists, a radius grid and a list containing the\n",
    "relative difference of the exact and calculated enclosed\n",
    "gravitational mass:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f8a89a4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "radial_grid=[]\n",
    "rel_diff=[]\n",
    "for i in range(1,tov_table.get_nlines()):\n",
    "    r=tov_table['r'][i]\n",
    "    radial_grid.append(r)\n",
    "    enc_mass=r*(1.0-1.0/b.exp2lam_from_r_gm(tov_table['r'][i],\n",
    "                                            beta))/2.0/b.G_km_Msun\n",
    "    enc_mass2=tov_table['gm'][i]\n",
    "    rel_diff.append(abs(enc_mass-enc_mass2)/enc_mass)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ebb56ce",
   "metadata": {},
   "source": [
    "Initialize the plotting object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e5e11ac7",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plots:\n",
    "    pl=o2sclpy.plotter()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9940eeb",
   "metadata": {},
   "source": [
    "Plot the enclosed gravitational mass as a function of\n",
    "radius for a 1.4 solar mass neutron star:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "033e8f4d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAGoCAYAAADiuSpNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA39ElEQVR4nO3de3Qb530n/C8o6n7hEJQli7YscWDJtzi2AMqxkzTeiKDlt47bNALtpkna7b4Wab+7yXZbSzB3322os3tMglJPt3upBcjddhu7jkRYdds3aWIOlasvEoGh5CS2ZQsX3aiLhcGQulG8YN4/VCAECZIYEsAMgO/nHBwJwwH4k8TBV8/MM7/HommaBiIiIhOpMLoAIiKiiRhORERkOgwnIiIyHYYTERGZDsOJiIhMh+FERESmU2l0AYW0cuVKrF+/3ugyiIjKVjQaxaVLl2bcr6zCaf369QgEAkaXQURUturr67Paj6f1iIjIdBhORERkOgwnIiIynZINp7a2NlgslrRHf3+/0WUREVEWSjqcNE1Le9TW1hpdFhERZaFkw4mIiIoXw4mIiEyH4URERKbDcCIiItNhOBERkekUNJxUVYXf74fb7c5q387OzgJURUREZlPQ3nrJvnaqqs64b3t7e56rISIisyroyMnpdEIUxRn3k2UZNputABUREZEZmfKak6IoWYUYERGVJtMtmSFJEpxOJyRJMroUIqKyNTQyhg/ODeJk7BrOxK9hZExDhcWCf+/cUJDvb6pwUlUVVqs16/3b2tqwa9eurPdfs2bNbMoiIioLkUtX8cNfnceVoVEsqKzAfbUr8Onbq/B/3X8rFlbOK2gtpgonn88HURQRDofR29sLWZZTI6lM2tra0NbWlvX7Z7vIFRFRubgxOoZ/ONqP0CdXIK5cit/dvBbCkgVGl2WOcFJVFYIgYOfOnZO+NlUwERHR7A2NjOHVw6egXL2B33rgNjxVv9boktIUdEKEJEnwer2QJAl+vz+1vampCbIsp57Lsoz9+/dDkiReeyIiyqFEQkNX4DT+m/QxHrt3NXZsvRt33brc6LImsWiaphldRKHU19en7rUiIio3H124jFffPYltjtvx6dsFQ2rI9nPYFKf1iIgof8YSGv73zyOYP8+C//yle1E5z5R3EaVhOBERlbBzA9ex98chfO3hddi42nyn76bCcCIiKlFvnbiEt0OX0Pqb92DR/MJOBZ8rhhMRUQl67cgpAMCOrXcbXMnsMJyIiErIyFgC/6PnYzws1uCzd640upxZYzgREZWIoZExdPzzh/iDz65H3cqlRpczJ+afsjFLbW1tsFgsaY/+/n6jyyIiyosrN0bx4vc/wLOP2oo+mIASDydN09IetbW1RpdFRJRz8avD6PjnD/BHzo24tWqR0eXkBE/rEREVMfXaMP6s+zh2Pn43Viyab3Q5OVOyIyciolJ3eWgEe94svWACGE5EREXp+vAYOn9wHH/ceFfJBRPAcCIiKjrDowl4fvAhvrnlTliXGr+8RT4wnIiIioimafizN4/j//58HVatKI3JD5kwnIiIioj3p2H85v1rsNa6xOhS8orhRERUJA7KZ2C7ZRkeWCsYXUreMZyIiIrA2ycu4drwGBrvXW10KQXBcCIiMrnTyjW8FbqErz+8zuhSCobhRERkYteHx+D7aRjfathgdCkFVbLhxN56RFTsNE3Dn0sf4d9tuRMLK4trPaa5KulwYm89Iipmr7x7Eo/duxqrS3jK+FRKNpyIiIrZO6EY5s+rQP16q9GlGILhRERkMpeu3MCPjl/E05vXGl2KYRhOREQmkkho+J+HTuBbDRtgsViMLscwDCciIhP523ei2Ga/HcsWlveKRgwnIiKTOHZaxbwKC+6/vcroUgzHcCIiMoHBoRG8cfRsWd1oOx2GExGRCfyvQyfwzS3lfZ1pPIYTEZHBvvfeOXz2zpUluzbTbDCciIgMdHFwCB+eH8SjG28xuhRTYTgRERlE0zS89JMQnn3UZnQpplOy4cTeekRkdv7gGTxx/xosLfNp45mUdDixtx4RmdVp5RouDA6VbXuimZRsOBERmVUioeGvfh7B9i+IRpdiWgwnIqICOxA4jd/ZdFvZLYOhB8OJiKiAzg8M4cLgDTywVjC6FFNjOBERFYimadj3szC2f6HO6FJMj+FERFQg3//FeXzxrlVYsoCz82ZS0L8hVVUhSRJ6e3vh8Xgy7iNJEsLhMEKhEGw2G5qbmwtZIhFRXsSvDuOX/QNwP3630aUUhYKGUyAQAHAzpDIJh8MAkAokm80GURThdDoLUh8RUb74fhbGs1/gzbbZKuhpPafTCVGceuqkLMvwer1p+8uyXIjSiIjy5u3QJdxXuwJVS+YbXUrRMNU1J5fLhX379qWeh8PhacOMiMjshkcT6H7/Ap64f43RpRQVU4UTAAiCAODmqT9FUeByuYwtiIhoDv72nSh+/5H1XApDJ9OFU5Lb7UZPT8+0+2Tqnzfdg731iKiQzsSvYWRMQ93KpUaXUnRMGU6dnZ1wu92pUdRUMvXPm+7B3npEVEh/81YUf/i59UaXUZRMEU7jZ+9JkgSXy5W61iRJkkFVERHNnvT+BXxuw0osms8WRbNR0HCSJAlerxeSJMHv96e2NzU1QZZlyLKMpqYmOBwOVFdXo7q6upDlERHlxPXhMRyOxPDFu1YZXUrRsmiaphldRKHU19en7rUiIsqXv/zxCXz5wdtQKyw2uhTTyfZz2BSn9YiISkXk0lUsrJzHYJojhhMRUQ698u5JfP3hO4wuo+gxnIiIcuQnH32CR8QartOUAwwnIqIcGBlL4CfHP0HDPZwEkQsMJyKiHPjukVP46kNr2QkiRxhORERzpFwdxqUrw9iwernRpZQMhhMR0Rz99VsRdoLIsZINp0x999hbj4hy7YNzg7i1ahGEJQuMLqWklHQ4sbceEeWTpmnY33saT9evNbqUklOy4URElG8//NV5NN67GpXz+FGaa/wbJSKahRujYwhE4/jcnSuNLqUkMZyIiGbh7w6fwtceXmd0GSWL4UREpJN6bRjxayNcRDCPGE5ERDr97Tsn8QePcNSUTwwnIiIdTivXsHj+PNQsW2h0KSWN4UREpMMrh0/i67zWlHcMJyKiLL13RsWGVcuxeAG7jucbw4mIKAuapuGgfBa/s+k2o0spCyUbTmxfRES51PPBRXzx7lWYV8Gu44VQ0uHE9kVElAujYwm8FbqERzfeYnQpZaNkw4mIKFf8wTNocrB/XiExnIiIpnH1xihOKtdwb+0Ko0spKwwnIqJpfOfdk/gGp44XHMOJiGgKsSs3MDyaQK2w2OhSyk6l3hccPXoU4XAY4XAYgiDAarVCFEU8+OCDeSiPiMg433n3JP7gkfVGl1GWsgqnaDSKjo4ORCIRiKIIURQhCAI0TUMoFMKbb76JcDgMm80Gt9uN9evX57lsIqL8Oqtex5IF81C9lCvcGmHGcNq9ezcURYHH40FVVdW0+w4MDMDn86G6uhrPPPNMzookIiq0V989iX/7xTuNLqNsTRtOu3fvhsvlQl1dXVZvVlVVhR07diASiWDPnj14/vnnc1IkEVEhnbh4GWuExVi6UPeVD8oRi6ZpmtFFFEp9fT0CgYDRZRCRyb34/Q/wJ49txMJK9tDLtWw/h3M2Wy8ajebqrYiIDPPeGRV3rV7OYDLYrMNpcHAw7eHxeHJZ15yxtx4RzcYbff34Mpu7Gk53OO3btw8VFRWorq6GIAipX30+Xz7qmzX21iMivd4+cQmfEa1s7moCusMpFAohHo9jbGwMiUQi9euOHTvyUR8RUUFomobuDy7gsXtXG10KYRbh1NjYmHFKeWtra04KIiIywpvvX0DjPathsXDUZAa6w8lisWSc/NDV1ZWLeoiICm4soeFwWMFn71xpdCn0L3RP4t+7dy/6+voAAKIoAgBisRgikQhvvCWiovRG31n89oO8Jm0musNJVVV0dHRAEITUNk3T0NnZmcu6iIgK4sboGI5fuIxtjtuNLoXG0R1OHo8HmzZtmrS9pqZmxteqqgpJktDb2zvl1HNZllNNZVVVhcvl0lsiEVHWbi4kyGAyG93XnDIFEwBUV1fP+NrkXcGqqk65T3t7O1wuF5xOJ3p7exEOh/WWSESUlaGRMZxSrmHD6uVGl0IT6A6niTffJh9ut3vG1zqdztR1qkwkSYLVak09t9lskCRJb4lERFnpCpzGU/Vcft2MdJ/WEwQBFosF41vyWSyWtGtQs5U8nZdktVrR3d095/clIppoaGQMZ9Uh2G5ZZnQplIHukVNzc3PqxtvkIxAI5GSEo6rqpGtXiqLM+X2JiCb67pFT+N3NHDWZle5wyjSRYdOmTYhEInMuRhAExGKxtG3jT/NNlKl/3nQP9tYjIgC4NjyKi5dvYP3KpUaXQlPQHU5TLTiYixGOKIppkyUURYHD4Zhy/0z986Z7sLceEQHAa0dO46sP3WF0GTQN3dectm7dOmlbOBye05RvVVUhCAKcTie8Xm9qezAYzGqiBRFRtq7eGIVy9QbWWpcYXQpNQ3c4aZo2KTBEUcxqtVxJktDV1QVJkuD3+1OB1tTUBI/HA7vdjtbWVvj9fgiCgMbGxmln9xER6fXakVMcNRUB3Svh9vX1TXmvk9lxJVyi8nZ5aAS+n4bxJ4/dZXQpZStvK+FWV1fj6NGjAG7e87Rnzx7s2bNHd4FERIX2d4dP4fc+w1FTMdAdTh0dHamuDQ0NDQiFQmhoaGBAEZGpDVwfwdXhMaypWmx0KZQF3decGhsb8ZWvfAWRSAShUAi9vb0AwDZDRGRqrx05ha9x1FQ0ZnVaD7g5ucHpdKa2c4EuIjKrgWsjuD48htUrFhldCmVJ98gpGAwiHo/D4/HA5/MBAHp6etjJgYhM65XDJ/G1hzlqKia6R047duyAoijwer3YsmULenp6IMtyPmojIpqz+NVhjI5pWLWco6ZiMuPI6dChQ9iyZUvatu3bt6d+39DQgIaGBhw6dCj31RERzdGrh0/i6xw1FZ0ZR05dXV1ZvVG2+xVKpr577K1HVF5iV24AAGqWLTS4EtJrxptwKyoqZpzsoGkaLBYLxsbGclpcrvEmXKLy8t97PsY3Hl6H6qULjC6F/kW2n8MzntZL9rqrr6+fspVQLBZDZ2enzhKJiPLnk8s3MK/CwmAqUjOGU/L6Ul9fH4LBIABMugZVVVWFpqamPJRHRDQ7rx4+iX/92fVGl0GzlPVU8vH99Hp6emCxWCCKItavXw/g5sQIIiIzuDg4hIWV8yAs4aipWOm+zwn4dRANDAzg0KFDUFUVTqcTK1asyGlxRESz8erhU/g3n595pQQyL933OY3X09ODjo4OuFwurrtERKZwbuA6liyYh6rF840uheZAdzgdPXoUzz77LCoqKtDe3o6mpiYkEgm89NJL+aiPiEgXdh4vDVmd1hscHITX64XX60U8Hsf27dsRCoWyWmCQiKhQzqrXsWLRfCxfxFFTsZtx5PTYY4+huroagUAAXq8XsVgMHR0dk4KptbU1b0USEWXjtcOn8FWOmkrCjCOnQCCAjo4OCIKAcDiMSCSS+lry5tt4PA6fz4f29va8FktENJXTyjVUL12AZQtnNc+LTGbGf8Xm5mbs2LFjxjcKhUI5KShX2trasGvXrrRta9asMagaIsq37/aewr/94p1Gl0E5MuNpvZaWlqzeyGyz9dra2qBpWtqjtrbW6LKIKA9Oxq7ilmULsWQBR02lYsZwynbSAydHEJFRvtt7Gk9v5rWmUjJtOO3evRvRaFT3m0YiEezZs2e2NRERZS38yRWsqVqExQvmGV0K5dC04bRjxw50d3ejtbU1q5AaHBzECy+8gNdffx3PP/98rmokIprSgcAZPFW/1ugyKMeyavwaiUSwd+9e9PX1QRRFCIIAm80GVVURi8WgqipCoRBsNht27tzJU3xEVBAnLl7GWutiLJrPUVOpyerqYV1dHTo6OgDcPGUXDocRDodRVVWFuro6iKKY1hiWiKgQuoJn8MeNG40ug/JA99SWuro61NXVsQs5ERnq+PnLWF+zFAsrOWoqRXNq/EpEZJTX5TPYZr/d6DIoTxhORFR03u8fxJ23LMOCSn6ElSr+yxJR0Xnj6Fn8jv02o8ugPGI4EVFR+eXZAdy1ejnmz+PHVykr2X/dtrY2WCyWtEd/f7/RZRHRHP3D0bP47QfZiqzU6Q6naDSKo0ePArh50+2ePXtM2Q2CvfWISs+x0yrurV2BSo6aSp7uf+GOjg6Ew2EAQENDA0KhEBoaGkwZUERUWv7xWD9+6wFeayoHuu9zamxsxFe+8hVEIhGEQiH09vYCQCqwiIjyIXgyjgfWCphXYTG6FCoA3SOn6upqAIAkSXA6nantFgt/YIgof/75F+fwxP1ck61c6B45BYNBxONxeDwe+Hw+AEBPTw8URcl5cUREANAbVWBfV81RUxnRPXLasWMHFEWB1+vFli1b0NPTA1mW81EbEREA4Ae/PI/H77vV6DKogHSPnKLRKDZv3owHH3wQg4ODkGUZFosFzzzzTD7qI6Iy904ohs3rrajgqKmszHm2Xjgc1jVbT5Zl+P1+SJIEv98/5T4+ny/1IKLypGkaut+/gK33rTa6FCqwgs/Wa29vR1dXFwDA7XbDbrdDFMW0fQKBAJqbmwEAfr8fsizDbrfrLZWIitzboRgesdVwwlUZKuhsPUmSYLVaU89tNhskSUrbR1XVVHgBgKIoEARBb5lEVOQ0TUPPBxfhvGeV0aWQAQo6Wy8cDqcFjdVqRXd3d9o+giBAFEXYbDZ4PB5YrdZJIysiKn0/+egTfGHjSo6aylRBZ+upqoqampq0bZlCzePxwG63Y/v27dOGXqb+edM92FuPqDhomoaffnQJj268xehSyCCzalC1ffv21Eq49fX1EEUR9fX1M75OEATEYrG0beNP8wE3A8zn86GrqwvBYBBer3fKiROZ+udN92BvPaLi0PPBRTTcs4qjpjI26+6Jg4ODiEajiMfjcDgc8Hq9M75GFEWoqpp6rigKHA5H2j4HDhxIXcsSRRE9PT2TTv0RUenSNA1vh2L4rK1m5p2pZOkOp56eHlitVtTV1cHhcMBut8PhcKCxsXHG1zqdzrTTdMFgMBVEydASRXHSqbxs3puISsMPf3Uej923mqOmMqc7nCRJgqIoiMViOHDgABRFmTTRYTqtra2p+5waGxtTkx2ampogyzKcTifC4TD8fj98Ph8OHDgAl8ult0wiKkKJhIbeaBwPixw1lTvds/XGX1uKRCIAgKqqqqxfb7fbM96zNP7UXfIeJyIqL9/7xTn85v1sU0SzGDkJgoCDBw8CAOLxOI4dOwYA7K9HRHMyltBw9LQKxzrrzDtTydMdTlarFS+++CKi0Siam5uxbds21NTUIBQK5aM+IioT/3SsH7/1AGfU0k26T+tt2rQJgUAg9fzEiRPo6+vDpk2bcloYEZWP0bEEfnl2AF/exFVu6aZZTyUfb9OmTTh06FAu3oqIytDf951lMFEa3SMnADh48OCk6d5erzfVBJaIKFvDowl8fPEKmurXGl0KmYjucHr22WdTU8eT3R0URUm7uZaIKFsH5TPYZr/d6DLIZHSHk8PhwN69eydt37dvX04KypW2tjbs2rUrbduaNWsMqoaIMrkxOobIpav43YfuMLoUMplZzdbLxGxdHDL13WNvPSJz6Qqc4ek8ykh3ONntdhw6dAjRaBSDg4Oph8fjyUd9RFSihkbGcCZ+HXeuWmZ0KWRCuk/rSZKElpaW1HOLxQJN02CxWPDSSy/ltDgiKl37e0/j6c0cNVFmukdOqqoiHo8jkUggkUhgbGwMiUQCO3bsyEd9RFSCrg2P4sLgEOpWLjW6FDKpWZ3Wy9RLr7W1NScFEVHpe+3IaXyVkyBoGrrDyWKxIBqNTtre1dWVi3qIqMRduTGK+NVhrLUuMboUMjHd15z27t2Lvr4+AEgtdxGLxRCJRPDMM8/ktjoiKjmvvnsSv/cZjppoerrDSVVVdHR0pK3fpGkaOjs7c1kXEZWg+NVhXB8ZQ62w2OhSyOR0h5PH48nY5LWmhouDEdH0Xnn3JL7x8Dqjy6AioPua01Tdx9mVnIimc2FwCPPmWVCzbKHRpVARyElXciKimbzy7kl8naMmylLJhlNbWxssFkvao7+/3+iyiMrSqdg1VC2ejxWL5htdChWJkg4n9tYjMoe/O3IKX/sMR02UvZINJyIyh48uXMbt1YuxeME8o0uhIsJwIqK86gqcxlPsPE465SycXn755Vy9FRGViGOnVdx96wosqOT/g0mfGe9z2rp164xvomkagsEgO0QQUZp/ONqP//TEPUaXQUVoxnDSNA1ut3vGN/J6vTkpiIhKwzuhGB6qq8a8CovRpVARmjGcpuoIMdFUK+QSUfnRNA1vvn8ef/qle40uhYrUjCeCZwqmgYEBvP7667BY+L8jIrrpR8cv4ot3reLnAs3arK9SDg4OIhqNIh6Pw+Fw8LQeEQEAEgkNP/v4En5jw0qjS6Eiprvxa09PD5qamlL/I0ou0b5v376cF0dExed7vziHJ+5fw1ETzYnukZMkSVAUBbFYDAcOHICiKAiHw2lLaJgB2xcRFd7IWALHTquoX89r0DQ3usOpvr4+9ftIJAIAGZdtNxrbFxEV3oHAaTy1mTfc0tzpDidBEHDw4EEAQDwex7FjxwAAsizntjIiKipXbozibPw6Nq5ebnQpVAJ0h5PVasWLL76IaDSK5uZmbNu2DTU1NQiFQvmoj4iKBJfEoFzSPSFi06ZNCAQCqecnTpxAX18fFxskKmOfXL6BkdEEl1+nnMlJw6tNmzYhGo3m4q2IqAh9592T+P1H1htdBpWQOd3nNP7h8XhyWRcRFYnopasQFs9H1RIuJEi5ozuc9u3bh4qKClRXV0MQhNSvPp8vH/URkcm9duQUfu8zdxhdBpUY3decQqEQ4vH4pOnjL7zwQlavl2U5dV+UqqpwuVwZ9/P5fBBFcdp9iMhY751RsWH1ciyaz4UEKbd0j5waGxsz3tfU2tqa1evb29vhcrngdDrR29uLcDg8aZ+WlhY4nc7UPpymTmQ+mqbhjb5+/M6m24wuhUqQ7nCyWCwZJz90dXXN+FpJktK6l9tsNkiSlLaPqqoIh8MQRRHAza7odrtdb5lElGc//fgSPr+hhktiUF7oPq23d+9e9PX1AUAqQGKxGCKRyIyLDU5sc2S1WtHd3Z22jyRJEAQBfr8fAKAoCpqbm/WWSUR5lEho+NGHF/HtJ7kkBuWH7nBSVRUdHR1pIaNpGjo7O7N6bU1NTdo2RVHSnofDYYTD4dR1ppaWFkiSBKfTOen92trasGvXrqxrX7NmTdb7EtHU/vFYP37rwVo2d6W80R1OUy0+ODF0MhEEYVIniYmLFAqCkNa/z2azobu7e8pwamtry7Ly9L6ARDQ7N0bH8P65QXyZ15ooj3Rfc0oG0+DgII4ePYrBwcG07dNJzr5LUhQFDodj0j5EZF5/d/gUvvoQp45Tfs3qJtxnn30WgiBgy5YtqK6uxtNPP53V65xOZ9ppvGAwmBoRJUNr4j6hUAiNjY2zKZOIcky5OoyB6yOoW7nU6FKoxOkOp927d6OxsRGJRAKKomBsbAxPPfUU9uzZk9XrW1tb4ff7IUkSGhsbUyOlpqam1JTx1tZWdHZ2wu/3w+FwZDylR0SF9zdvR/GvP7ve6DKoDOi+5iSKIrZt25a2bdu2bVmvhGu32zNODR8/a2+qfYjIOKFPrsC6ZD6EJQuMLoXKwKzuc8okmwkRRFS8Xjt8Cl9lmyIqEN3hFAqFUpMgkqLRKI4cOZKzoojIXN4JxeBYV42FlWxTRIWh+7Rec3MztmzZAovFAqvVCkVRoKoqgsFgPuojIoMlEhrefP88/vRLvOGWCkd3OFVVVSEQCMDv9yMSiWS8BkVEpeONo2fx2w/exhtuqaB0h1PSxE7hL7/88ozti4iouFwfHsPxC5fxFfvtRpdCZWbGcHruuefQ1NSELVu2AAC2bt06aR9N0xAMBhlORCXmO+9G8Y2H1xldBpWhGcNJ07RJz91u96T9vF5v7qrKgUx999hbjyh7FweHMDyawO3VS4wuhcqQRZuYPjPo6+vL2Kpoqu1mUl9fj0AgYHQZREVh9w8/RMujNqxYxOXXKXey/RyedW+9pOTaTmYPJiLK3i/ODOAO6xIGExlGdzi9/PLLac81TUNPT8+k7URUnDRNw+vyGbgca40uhcqY7nAa31UcAOrq6tDQ0JCreojIYP94rB9PfHoNV7glQ2U1lXz37t3o7e2FxWKBLMvo7e1N+7qqqhBFkbP1iIrcteFRvH9uEL/9INdqImNlFU47duwAALzwwgvQNA1PPfVU2tdFUeQ1J6ISwK7jZBa6bsLt6OhAT08PT+MRlaDTyjVUVliwpmqx0aUQ6b/mNDGYBgYG8Prrr+Po0aO5qomIDPC370Tx+4+sN7oMIgCzXAkXuLlMezQaRTweh8PhMN1NuESUvbdPXIL9jmosms+u42QOunvr9fT0oKmpKdUEUtM0WCyWrBcbJCJzGR1LoPuDC+w6Tqaie+QkSRIURUEsFsOBAwegKArC4TAEQchDeUSUb9/tPY2nN69l13EyFd3hVF9fn/p9JBIBcHMZDbNpa2uDxWJJe/T39xtdFpGpxK8O48LgEO6+dYXRpRCl0R1OgiDg4MGDAIB4PI5jx44BAGRZzm1lc9TW1gZN09IetbW1RpdFZCp/9fMI/s3n6owug2gS3eFktVrx4osvIhqNorm5Gdu2bUNNTQ1CoVA+6iOiPDl6WsUdNUtQvXSB0aUQTaJ7QsSmTZvSOsqeOHGiKDqSE9GvjSU0/L18Bt9+8j6jSyHKaNZTyZOi0SiDiajI7O89jac2r0UF++eRSbErOVGZuXTlBs4PDuG+WvNNZCJKYldyojLzVz+P4Jnf4CQIMjd2JScqI71RBXffupyLCJLpsSs5UZkYHUvge++dw7efZCcIMj92JScqE68ePoWvfeYOdoKgojDnruRJnBBBZF7nB4YweH0EG1YvN7oUoqzMOHJ67rnn0NTUhC1btgAAtm7dOmkfTdMQDAZNdc2pra0Nu3btStu2Zs0ag6ohMtbLPwvjjx/baHQZRFmbMZw0TZv03O12T9rPbEtmtLW1oa2tLW3b+L6AROXizV+dx+c2rMSSBbrvuScyzIw/rXv37k177vF4Mk5+sFqtuauKiHLi8tAIgifjaP3Ne4wuhUgX3decfD5fxu2crUdkPvt+FsEzvyEaXQaRbrrDaf/+/Xj55ZcxODiYj3qIKEfkU3GsrV6MW5YvNLoUIt10n4Tet28ftm3bhp6eHoTDYdhsttRkCSIyh5GxBP7xaD/vaaKipXvktG3bNgA3p5Rv374ddXV12Lp1K/bs2ZPz4ohodv7P21H8/iPreE8TFS3d4XT06NHUr88++ywcDgc0TYPdbs91bUQ0C5FLV5HQNIi3LDO6FKJZ031ar6mpCYIgIBwOo7W1FZFIRNcy7bIsIxwOQxAEqKoKl8s15b6qqsLn82Hnzp16yyQqS5qm4f+8HcV/5Ow8KnK6w0nTNHR0dMy6hVF7ezu6uroAAG63G3a7HaKYeTZRe3v7rL4HUbnqCp7Bkw/UYkHlnJdqIzKU7p9gj8eTFkwDAwN4/fXXU6f7piNJUtr9UDabDZIkZdxXlmXYbDa95RGVrX71Os4PDMGxrtroUojmbNYTIgYHBxGNRhGPx+FwOLLqEJE8nZdktVoRDAYz7qsoypQjKiJKp2ka9v0sjOYv8Jih0qA7nHp6emC1WlFXVweHwwG73Q6Hw4HGxsYZX6uqKmpqatK2KYoyaT9JkuB0Omd8v7a2Nlgslqwf/f392f9BiYrI6/JZPHH/GiyaP8/oUohyQnc4SZIERVEQi8Vw4MABKIoyaUQ0FUEQEIvF0rZNbHukqmrWrZDa2tqgaVrWj9ra2qz/nETF4tzAdZyNX0f9erYQo9KhO5zGN0+NRCIAkPVsPVEU05Z5VxQFDocjbR+fz4dwOAy/34/u7m7IsjzldSmicqdpGvb9NMLTeVRydIeTIAg4ePAgACAej+PYsWMAbk5gmInT6Uw7jRcMBlOn75KhtXPnTrhcLrhcLmzevBl2uz2rU3xE5eiNo2fx+KduxeIFPJ1HpUV3OFmtVrz44ouIRqNobm7Gtm3bUFNTg1AolNXrW1tb4ff7IUkSGhsbU5Mempqa0gJOlmXs378fkiRx5ESUwYXBIUQvXcNDdTydR6XHok1csGkW+vr6iqIreX19PQKBgNFlEM2Zpmn4L//fB3h+60au00RFJdvPYd0jp4MHD066p6kYgomolPiDZ/D4p25lMFHJ0h1O3/3udzNu5xIaRIVxKnYNFy/f4Ok8Kmm6w+npp5/OeHPsVIsQElHujCU0/NXPw9jOBQSpxOk+J9Dd3Y0XXngBAFKdyDVNQ19fH55//vncVkdEaf76rQh+7zPr2DuPSp7ucAoEAti5c2fGm2eJKH9+eXYA8+dV4K5blxtdClHe6Q6niY1fk7LpEEFEszM0MgZ/8Az+9Etc2ZbKg+5zA8muEBM1NDTg6NGjWXUnL4RMfffYW4+K1d6fhND8BREVFVzZlsqDrnAaGBhAXV1dxq/t27cPzzzzDLxeLw4dOpST4uYiU9899tajYvTj4xdhu2UZaoXFRpdCVDC6wqmqqgrhcBj19fXYsGFDqo0RcHO23qFDh/DSSy9l1cqIiGZ28fIQDkcUPPkA/2NF5UX3ab1QKASPx4M333wTR44cSd3fFA6HsWLFCgDZN4IloqklEhr+8kchfHPLnUaXQlRwusPpoYceQkNDA+rq6tDR0ZGxDYXFwvPiRHP1N29H0VR/O7tAUFnSHU5HjhzByZMnMTg4iD179kAURQwMDGB8iz5OKyeam75TccyvrMB9tTwLQeVJ93/JWlpa0NDQgHg8jubmZnR3dyMYDMLj8aTCiqf1iGZvcGgE/3TsHP7zl+4xuhQiw+gOp7q6Opw4cSJt2/bt2wHc7E4eCARSz4lIH03T8D8PncC/23InT49TWcvpyexNmzaxQznRHBwInMaWu1fBunSB0aUQGYoNuohM4r0zKq7eGMPDYo3RpRAZjuFEZALqtWH8fd9Z/OHn1htdCpEpMJyIDJZIaPiLno/xR86NvM5E9C9KNpzYW4+Kxf9+K4Jt9ttRtXi+0aUQmUZJhxN765HZvR26hBWL5uNTt/H2C6LxSjaciMzu3MB1/OSjT9BUf7vRpRCZDsOJyABDI2P4yx+F8B94nYkoI4YTUYFpmoY/lz7Cc//KhkXz5xldDpEpMZyICuyVd0/Cec9qrs9ENA2GE1EBvXXiEhZUVmDzeqvRpRCZGsOJqEBOxa7hnVAMT2++w+hSiEyP4URUAFdujOLln4fxrYYNRpdCVBQYTkR5NjqWwJ4fHse/b9iABZU85IiywSOFKM/++6ET+PrD61CzbKHRpRAVDYYTUR69evgkHq6z4s5Vy4wuhaiolGw4sbceGe3QhxewqHIePnvnSqNLISo6JR1O7K1HRvnl2QF8cO4ytjnYmohoNko2nIiMcla9jjf6zuK5R21Gl0JUtBhORDkUvzoM709C2PH4XaioYM88otliOBHlyLXhUfxZ93Hs2HoXFlayZx7RXDCciHJgZCyBzh8cx7caNmD5Ii4aSDRXDCeiOUokNOz54XH84efWY9XyRUaXQ1QSKgv9DWVZRjgchiAIUFUVLpdr0j6SJCEcDiMUCsFms6G5ubnQZRJlRdM0/K8fncCTD9RiXc1So8shKhkFD6f29nZ0dXUBANxuN+x2O0RRTH09HA4DQCqQbDYbRFGE0+ksdKlEM3r5ZxHUr7dymXWiHCvoaT1JkmC1/nqpAJvNBkmS0vaRZRlerzf13Ol0QpblgtVIlK2/eSuCu25djkdsNUaXQlRyChpOydN5SVarFcFgMG0fl8uFffv2pb1m/MiKyAxePXwSa61L8IWNtxhdClFJKmg4qaqKmpr0/2UqijJpv2SAqaoKRVEyXpcCMrcomu7B9kWUCwcCp1GzdCEa7lltdClEJaug4SQIAmKxWNq28af5JnK73ejp6Zny65laFE33YPsimqs3+s5i8fx5ePxTtxpdClFJK2g4iaIIVVVTzxVFgcPhyLhvZ2cn3G532mlAIiO90XcWAPDkA/xPDlG+FTScnE5n2mm8YDCYmoU3PrQkSYLL5Upda5o4aYKo0A70nsa8Cgu+vOk2o0shKgsFn0re2toKv98PQRDQ2NiYCqCmpiZ4PJ7U78dLTj0nMsIr757EymULeSqPqIAKHk52ux12u33S9u7u7tTv4/F4IUsimtJfvxXBupol2HI3Jz8QFVLBw4moGGiaBt9Pw7ivtgqf38DFAokKjb31iCZIJDT8ufQxNt1RzWAiMghHTkTj3Bgdw+4fHEdT/Vrcdetyo8shKlsMJ6J/MTg0gj0/PI6WR224TVhsdDlEZY3hRATg4uAQ/sehE/iTxzZCWLLA6HKIyh7Dicrexxcu47Ujp/GfnrgHi+ZzBVsiMyjZCRGZ+u6xtx5N9OPjF/GDX55nMBGZTEmHE3vr0VQ0TcN33oni0pVhfLNhA+ZVWIwuiYjG4Wk9KjvDown8Rc9H+MKGW/AZkWsxEZkRw4nKSuzKDfxFz8d45vMi7qhZYnQ5RDQFhhOVDflUHN977xx2bL0LyxfNN7ocIpoGw4lKnqZpePXwKWiahv/3iXtgsfD6EpHZMZyopF29MYr/Jn2Erffdivr1Uy9sSUTmwnCikvXh+UG8dvgUvtmwASuXLTS6HCLSgeFEJSeR0PDK4ZMYHk3gT5+8j9PEiYoQw4lKyoXBIbz04xC+Yr8Nn75dMLocIpolhhOVjDd/dR7HzqjY+fhdWLKAP9pExaxkO0RQ+YhfHUbHP38Ii8WCHVvvZjARlYCSDSf21it9mqbhe++dg+9nYTz3qA2N93IpdaJSUdLhxN56pevC4BD+6/c+wIrFlXA/fjeqlvCmWqJSwvMfVFRGxxLYHziNTy7fwB83bsTShfwRJipFPLKpaByJKPj+L87hqfq1uLd2hdHlEFEeMZzI9M4PDOGv34rgvtuq8O0n72X7IaIywHAi07pyYxTfeeckEpqGbzVs4Ck8ojLCo51M58boGA70nsa5gSF845F1WFO12OiSiKjAGE5kGmMJDf90rB+/6h/A05vX4s5Vy40uiYgMwnAiw42MJfAPR/vx4blBPPHpNfjyptuMLomIDMZwIsMMjYzBHzyD08o1PPlALVyO240uiYhMguFEBRe7cgOvy2cQvzaCbfbb8PWH1xldEhGZDMOJCuaXZwfw/V+cw/JF87HNcRtWLV9kdElEZFIlG05tbW3YtWtX2rY1a9YYVE35uj48hh/+6jzePzeIe9eswB85N2JBZcl2zSKiHLFomqYZXUSh1NfXIxAIGF1GydM0DX2nVRz64CLmVVjw+KduxT1r2NGBiLL/HC7ZkRMVXviTK3jz/QuIXxvGprUCvtWwgaMkIpoVhhPNSTKQBq6PoG7lUvzu5rUQliwwuiwiKnIMJ9JleDSBQFTB4YiC0UQC62oYSESUewwnmpamafj44hUcjig4p17H/HkVeKjOiv/nizYsrJxndHlEVKIYTpRmLKHh+PnLOBKJ4dKVYVgswJ2rlsF5zyr2uCOigil4OMmyjHA4DEEQoKoqXC7XrPahuRtLaIhcuoL3zgwgGrsGaBoqKizYsGo5nvh0LW5ZvtDoEomoTBU8nNrb29HV1QUAcLvdsNvtEEVR9z6UvbGEhtPKNZy4eAUnPrmCazdGAQAWiwXiLUvx4FoBX37wNlRUcJ0kIjKHgoaTJEmwWq2p5zabDZIkobm5Wdc+lG5oZAyxq8M4p17Hmfh1nFWv48bI2M0vWiyosABrq5fgzlXL8LCtBsu4LhIRmVxBP6WSp+qSrFYruru7de9TasYSGoZHE7gxOoahkQSu3BjB4NAoLg+N4vLQSOrXqzfGkOmO6YWVFahZugC1wmLcf3sVHv/UrVg0n5MViKh4FTScVFVFTU1N2jZFUXTvY4SffPQJ5JPxrPfPtJJ4sheHxfLr3wPAvAoLFlRWYGFlBRZWzsOyRZVYvqgSKxZVYk3VIixfVInli+Zj6YJ5XKKciMpCQcNJEASEQqG0beNP4WW7T1Km/nnTmUtvvUc33oJHN94y69cTEVH2CtpbRhRFqKqaeq4oChwOh+59ktra2qBpWtaP2trafPyxiIgoxwoaTk6nM+0UXTAYhNPpBIBUIE23DxERlYeCdyWf6h6mxsZGeDwe2O32vN3nxK7kRETGyvZzmEtmEBFRwWT7Ocz1DIiIyHQYTkREZDoMJyIiMh2GExERmQ7DiYiITIfhREREplOy4dTW1gaLxZL26O/vN7osIiLKQkmHE9sXEREVp5INJyIiKl4MJyIiMh2GExERmU5Z9dZbuXIl1q9fn7atv78/62tR2e7L9zT39+Z7lud7ltqfp1jfMxqN4tKlSzO+pqzCKROLxYJs/wqy3Zfvae7vzfcsz/cstT9PKb7neDytR0REpsNwIiIi02E4ERGR6ZR9OH3729829D2z3Zd15vZ7s87cfu9iqJPHUG73zUed45X9hAg9ZnNRzwjFUGcx1AiwzlxjnblVynWW/ciJiIjMh+FERESmw3AiIiLTYTgREZHpMJyIiMh0GE465HvqZK4UQ53FUCPAOnONdeZWKdfJqeRERGQ6HDkREZHpVBpdQDGQZRnhcBiCIEBVVbhcLqNLmkSSJITDYYRCIdhsNjQ3Nxtd0rRUVYXP58POnTuNLmVKPp8Poiia9t8cuPmzGQgEUs/N8u+uqiokSUJvby88Hk9qu9mOpanqNNvxNFWd479u9PE0XY2zOZYYTllob29HV1cXAMDtdsNut0MURYOr+rVwOAzg1x9MNpsNoijC6XQaWda02tvbjS5hWi0tLXC73RBFMfWr3W43uqxJAoFA6t/d7/dDlmVT1JkMTFVV07ab7VjKVKcZj6ep/j6TzHA8TVXjbI8lntabgSRJsFqtqec2mw2SJBlY0WSyLMPr9aaeO51OyLJsYEXTk2UZNpvN6DKmpKoqwuFw6kPT4/GY4gN/IlVVUx/0AKAoCgRBMK6gcZxO56TQMeOxlKlOMx5PmepMMsvxlKnGuRxLDKcZJE9BJFmtVgSDQeMKysDlcmHfvn2p5+N/GMxIURRT1ydJEgRBgN/vh9/vh8/nM7qkjARBgCiKsNls8Pv9sFqtpv57LYZjCeDxlEtzOZYYTjNQVRU1NTVp2xRFMaiaqSUPelVVoSiK4efypyJJkqlPNwI3P4zC4TBcLhdcLheCwaDh/8OfSvJ/otu3bzflz+V4xXIsATyecmUuxxLDaQaCICAWi6VtG39qwmzcbjd6enqMLiMjVVVN/XeXJAgC6uvrU89tNhu6u7sNrCiz5EXwrq4uBINBeL1e+P1+o8uaUrEdSwCPp7may7HEcJpBcoZJkqIocDgcxhU0jc7OTrjdbtNcd5jI5/MhHA7D7/eju7sbsiybckRi1lMkEx04cCD1v2ZRFNHT02PKEE0qpmMJ4PGUC3M5lhhOM3A6nWmnHoLBoCmH0ZIkweVypX4YzPZDCgA7d+5MDe83b94Mu91uyr/Lif/moVAIjY2NBlaUmSiKk06LmbHOpGI5lgAeT7kyl2OJHSKyYLZ7MyaSZRkNDQ1p27q6ukz3g5okyzLa29sRDofh8XhMWWfyf6HJADD6Ppep+Hw+WK3W1AeAWeqUJAldXV2QJAkejyd1zJjtWMpUpxmPp6n+PgHzHE/T/ZvP5lhiOBERkenwtB4REZkOw4mIiEyH4URERKbDcCIiItNhOBERkekwnIiIyHQYTkQlqrOzs6jfn8obw4koD1RVRVNTEywWS9r2lpYWtLS05P37+3y+1E2QkiTB4XCgqakpp9/D5XIxoChvGE5EeSAIQtpaS0lNTU05D4mJkiu4JlvvOJ3OvASiKIqIxWKpxfmIconhRFRATqcz7+1lvF7vpDDKV/fq1tbWjMuGE80Vw4moxMiyXLDO6oIgcOREeVFpdAFEZiJJEtxuN+rr69HY2Ij9+/fj6aefTjUEVRQFqqqit7cXjY2Nk0ZBbrcbNpsNVqt10lILsizD7XYDQGqJg/HPVVWF2+2Gz+dDPB5Pvb6zsxN2ux2qqqK7uxtut3vK8Mlm1VZVVVFXV4f6+nq43W5YrdbU0hDJEZcsyxAEIbVEuaIoqTWjJhJFEbIsm3IpeypiGhGl8Xq9miiKWjwe14LBoBYMBjVN0zS73a51dXWl9hMEIe11Tqczta+maVooFNImHmLd3d2a0+mc8rmmaRoALR6Pp2oZ/z27urrSvsdE3d3d2s6dOydt7+rq0lwuV6our9c76euiKGqhUEjTNE2Lx+MaAK27uzvtzze+liSPx5NxO9FccORENEHy+owgCGmjga6urkmjElVVIQhCaimI8fvn4tSaKIpoaWmBoihwOp0zLjERDodhs9mm/Hpy+YKdO3embU+O0pI1T3ye3JZpWXWe2qN84DUnogwynaKyWq3o7OyEz+dLLT6X/LAOBAJ5WTHV6XTC6/Wiu7sbDocDDocjbTXZTKb6uizLCAQCaG9vz7hPpjAdP5FiukkVZl0tlooXw4koSw6HA06nE83NzWnXmlRVhSiKeRk9SJIEp9OJrq4uxONx1NfXw+fzTbm/1WpFLBbL+LVk7a2trdi+fXvOalRVNW+zAal8MZyIspCcFJAcUY0feSQDJDkxYPxrZjIx1Ca+JnkaLmmm+5VEUZxxZLVz507Isgy/35+2PdMpu2y+Pv6eKqJcYTgRjSNJErxeL2RZRmdnZ+qD3m6346mnnkJnZyckSUIgEMC+ffvg8XhSH8w9PT3Yv38//H4/JElKfZi3tLRAVVXIsgyv14tAIJAa/YiiCJfLBZ/PlwoLQRCwffv21HLm4XAYfr8/9b4TrxeNZ7fbJ43gkn8mSZLSvkdTUxM6OzshyzI8Hg/C4TB8Pl9q1iBwc/Zhcnuy/vFhCWDStTaiXOAy7UQlpqWlBR6PpyDXgZJBlmmKOdFccOREVGLcbjfa29sL8r18Pl9qlEWUSwwnohIjiiJqamryPr07edqR15soHxhORCVo586dkyY85Jrf70dzc3NevweVL15zIiIi0+HIiYiITIfhREREpsNwIiIi02E4ERGR6TCciIjIdBhORERkOv8/ulRohNYIwbMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if plots:\n",
    "    pl.canvas()\n",
    "    plot.plot(tov_table['r'][0:tov_table.get_nlines()],\n",
    "              tov_table['gm'][0:tov_table.get_nlines()])\n",
    "    pl.xtitle('radius (km)')\n",
    "    pl.ytitle('gravitational mass (Msun)')\n",
    "    plot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbd3bded",
   "metadata": {},
   "source": [
    "For the enclosed gravitational mass, plot the relative\n",
    "difference of the exact results and that computed from\n",
    "the tov_solve class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "69c45734",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAGyCAYAAABEAduNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4yUlEQVR4nO3df1wb950n/pfA2MZ2zCCcH3Z+mVEcx+5uGwT2tr22SUHEaZI+7i4W8d1+7/b2u18bnO1jt9ddB5n0l9y91hbxttu77/ccIHvXx6a3lwTZu3dN2iQauG02bRqQZKdNnNiBAccO8S+kAfyTH5rvH0QKAoEYGM2MpNfz8dAj0TAa3jaWXnxmPvP+2FRVVUFERGQBBWYXQEREFMdQIiIiy2AoERGRZTCUiIjIMhhKRERkGQwlHbS2tkKSpFmfExHR/DCUdFBVVQVFUWZ9TkRE88NQAuD3+1FbWztjmyRJaG1tNakqIqL8w1AC4Ha7k577/X4AgMvlAoDEqTi/35/04GiIiEhfeRFK4XAY4XA48by1tXXOQOnu7oYoigAAURQTr3W73UkPQRAATIZWd3d34pjTnxMR0fwsMbsAIzidzsToJxgMwuVyJQIllelhMjg4OOfxGxsb53xORETzkxcjJWBylNPS0gIAiVHQbARBQCQSMaIsIiKaIm9Cye/3o6GhAQAgy/Kc+27ZsiUxWpJlecYkCCIiyoy8CKVwOAxBEOB0OlFfXw9JkpJO0UmShGAwmDjF53a7IctyYr/4hAciIsosG5euICIiq8iLkRIREWWHvJh9t2bNGqxfv97sMoiI8lZ/fz8uXryYdr+cCyWv14t9+/YlbVu7di2CwaBJFRERUVVV1bz2y7nTd16vF6qqJj3WrVtndllERDQPORdKRESUvRhKRERkGQwlIiKyDIYSERFZBkOJiIgsg6FERESWwVAiIiLLYCgREZFlMJSIiMgyci6UvF4vbDZb0mNgYMDssoiIaB5yMpTYZoiIKDvlXCgREVH2YigREZFlMJSIiCilnvOX8MNXT+CjoauGfc+cW0+JiIgW58TZERwJn8EdZSvwOccavPvRMNaWFBvyvRlKREQEAOi7eBnPdX2A8jUr8ZcPbMTSJQU4N3wNgePnDKuBoURElOfOj1zDs2+cQtnKpfhG7d1YXlSY+NqaVctw8dJ1w2phKBER5alrYxP46W9O4eroBHZ+UURJcdGMfQoLbIipxtXEUCIiyjOqquLV4+fQ1RfBH33uTtxZtjLt/qqqwmazZbw2zr4jIsoj758bwfdePI7iokJ8+5HNaQMJAF4InsZvzwwZUF0OhhLbDBERzXRtbAL/3//pwes9F/HkQ5vwpbtvnPdrVy5bAvHG9OGlB5uqqgaeLTRHVVUVgsGg2WUQEZniN/IgXn77LP6fL5TjdvsKza+fiKkoLFjcqbv5fg7zmhIRUY4aujKGltd6sfGWG/Ddr25e8DWhxQaSFgwlIqIc9PLbH+HY6SHUf0mEfeVSs8uZt5y7pkRElM+UK6P4wc/fxfKiQuz9yj1ZFUgAR0pERDnjtZMX8Kuei/ja/XehZMXMe46yAUdKRERZ7sroOH746gmMXBtH00ObsjaQAANHSn6/H6IoQpZldHd3w+fzzdgnHA5DlmUIggBFUeB2u40qj4goK4U/iOJnbw1g930O3Lx6udnlLJohIyVFUeDxeOB0OuF2uyFJEsLh8Iz99u/fD7fbDZfLhe7ubsiybER5RERZJxZT8cw/y3jvoxF855HNORFIgEGhJAgCQqFQ4rmiKHA6nUn7SJIEu92eeO5wOCBJkhHlERFllfMj1/C9F4/j8441+MM/uMOQ9j9GMez0nSAIAIDW1taUp+7ip+3i7HY7AoGAQdUREWWH105ewBvyIBof3IgVS3Nvrpqhf6JwOJw0GppKURSUlZUlbYtEIkaURURkeWMTMbS+JuO20mJ4HrzH7HIyxtDZd/FrSi0tLWhtbU36miAIGBwcTNo2W4Cl6m8314O974gom32oXMVfvXgcX/m9W/Av773V7HIyypBQ8vv9aGhoSDwXRTHpGlN8m6IoieeRSASVlZUpj+f1ehOt1OfzWLduXUb+XEREmfbayQv4+zdP4cmHNkG8cZXZ5WScYRMdamtrE89lWUZdXR0AJILI5XIlna4LhUJwuVxGlEdEZDnx2XXnR67jiW33JK0Gm8sM6xLu9/sBTI6AFEVBY2MjAKC2thY+nw9OpzNj9ymxSzgRZZPha2P4UeAktjtvw+/dWmJ2ObqY7+cwl64gIrKQE2dH8D+7PsDXazagNMv61s2FS1cQEWWZn701gDPRq/j2I5sNXS7CStj7jojIZGMTMfyNdBLFRYV4/H5H3gYSwJESEZGphq6M4YeBE/gPn1+fF7Pr0mEoERGZpOf8JTz7Rj/+4oGNKCnO3s7eemIoERGZ4LWTFxA6FcW3H9mMJYW8khLHUCIiMpCqqvjpmx9g2ZICfKP2brPLsZyci+dULYjYZoiIrGB0PIanXjmBDTetwmNVt5tdjiXlZCixzRARWU3k8ij+00vH8W+23IHPimXpX5CnePqOiCjDTp4bwd+/+QGe2LYRNyznhIa5MJSIiDKo871zeOfD4by+IVYLhhIRUYY8+5tTWLm0EH9Ws8HsUrIGQ4mISGcTMRU/7ngfnxXt+LxjjdnlZBWGEhGRjq6MjqP55RP4d5+9A3fddIPZ5WQdhhIRkU7OD1/Df+nswZ/XbMCNNywzu5ysxFAiItLByXMjeK7rNJ58aBOKl+bHgnyZwFAiIlqk19+/iNCpKL758CbOsFskhhIR0SIcCZ/BtbEYvu7iDDs95GRHB7YZIqJMU1UVz/yzjFXLluAP/+AOs8vJGTkZSmwzRESZNBFT8aPASXzmdgEPfOoWs8vJKTx9R0SkwfXxCTz18gnUVd2OjbdwyrfeGEpERPM0fG0MB185gYb7HLhVKDa7nJzEUCIimof4PUh/UXs3SlcuNbucnMVQIiJKQ75wCX/3xineg2QAhhIR0RzeOq3g5XfO4lsPb+Ky5QZgKBERzeKXJy/gnYEhNG7bCJuNN8UagaFERJTC/zr2IS5fn8Cf3n+X2aXkFY5FiYim+cmv+rC0sIA3xZqAIyUioo/FYir+S2cPtpSXch0kk+TcSIlthohoIUbHYzjw8ntwbb6JgWSinAwlthkiIi0uXR/H9186jn//2TvxqXUlZpeT13j6jojy2sVL1/Fj6X38R9cGlK3iwnxmYygRUd76YPAK/vZ1GXu/cg9WLuPHoRVo+in09/dDURTce++9GB4eRmtrKwBgz549GSmOiChT3v5wCD97awDfemQzinhTrGVo+kkcOHAAsiwDAGpqatDb24uamhocPHgwI8UREWXCr3su4v+8dx6eB+9hIFmMppFSbW0tHn30UfT19aG3txfd3d0AkAgqIiKre+m3H2Hw8nX8WQ1XirUiTb8ilJaWAgAkSYLL5UpsZ/sNIsoG/+PNU5hQVfzR59abXQrNQtNIKRQKIRqNwufzJa4ndXR0IBKJZKQ4IiI9qKqK//pPvfj9W0vwpbtvNLscmoOmkdITTzyBSCSClpYWVFdXo6OjA+FwOFO1EREt2thEDL6XT+CLG9YwkLKA5tl3W7ZsScy+C4fDsNls2LlzZ6bqIyJasKujE/C9/B7++PPrsX7NSrPLoXlY1Ow7WZYtN/uObYaICACil0fxg5+/iz/9soOBlEU0hdL02XeHDh1CRUUFysvLM1WfZmwzRERnolfwI+kkGh/ciJtuWG52OaSBptN3nH1HRFb33tlh+INn8M2HN2HZEi5dnm04+46Icsab8iB+3TuIJx/ahIIC/rKcjTj7johywivvnMU7A8P4j64NDKQsprkD4a5duxL/X1VVBUVR4HA4dC2KiEiL57s/QFFhAf7kC9a5vk0Ls6CmT8PDw+jv70c0GkVlZSVaWlr0rouIKC1VVdHyy17cdMNyPOq8zexySAeaRkodHR2oq6tLTGxQVRU2mw1tbW1pXytJEmRZRm9vLxwOB+rr62fs09DQgLq6OlRVVcHj8TDsiGhWEzEVPwycQM2mm+G8o9TsckgnmkJJkqTEpIaOjg7U1NRgaGgIoVBoztfF722KB5HD4YAoikkz+AAgEokkQomBRESzuTY2geaXT+AP/+AO3HXTKrPLIR1pOn1XVVWV+P++vj4AQElJ+qWDw+FwUsi4XK6UEyR27NiBaDSKQCAAURS1lEZEeWLoyhi+/9K7qP+SyEDKQZpCSRAEHDlyBAAQjUbx1ltvAUDaGXhutzvpFJ8syylDR5Zl+P1++P3+xJRzIqK4s0PX8NeBE9izbSNuKeFNsblI0+k7u90Oj8cDp9OJ+vp6VFZWIhqN4rHHHkv7WkEQAACKoiASicDtds/Yp76+PrFfXV0dXC4XR0xEBADoOT+Cv3/zNJ58aBOWF/Gm2FylaaRUUVGBYDCI9evXo6SkBD09PZAkCYcOHZr3MTweDzo6OlJ+LR5IACCKIvx+f8r9UvW3m+vB3ndE2S10Kop/PDqAbz7MQMp1i14HuKKiAp2dnfPat7m5GR6PJyl84iRJQkNDw7yOk6q/3VwP9r4jyl6d751D6FQEf/nA3SjkTbE5T/PNs0eOHJnRVqilpSWxNPpsJEmC2+1OnI6L989TFAWCIEAURdTV1SX2nz45gojyz5HwGYzHVNR/iTfo5wtNobR7927IsgxBEGC32wFMTuNWFGXO14XD4aTAAYD29nYAk9eOfD4fnE4nwuEw/H4/ZFlGQ0MDrycR5SlVVfG3r/fhDvsKPPCpW8wuhwykKZQqKyvx9NNPz9ie7uZZp9OJaDSa8muBQCDx/6kmPxBRfpmIqfgb6STuu/tGVK23m10OGUzz7LtUamtrdSmGiPLbJzfF3o67brrB7HLIBJomOjidTnR2dqK/vx/Dw8OJh8/ny1R9RJQnopdH8f2X3sXu+0QGUh7T3GZo6gw5m82W6H+nZVo4EdFUpyNX8Lev96HxwY24YXmR2eWQiTSFkqIoiEajM1oL7d27V9eiiCh/vP3hEH722wE8+dAmLF2y6LtUKMtpCiWn05my111TU5NuBRFR/vjlyQt4+8Mh7H3wnsTqA5TfNP1aYrPZ0N/fP2N7fHo3EdF8/cPRMxhQruJrX76LgUQJmkZKTz/9NI4ePQoAiXuIBgcH0dfXh507d+pf3QJ4vV7s27cvadvatWtNqoaIplNVFW3/LGN92Ur86wreg0TJNF9TOnDgQFKbIFVV0dzcrHddC+b1euH1epO2TV1yg4jMMzoew18HTuArv7cW994umF0OWZCmUPL5fKioqJixvaysTLeCiCg3DV0Zw18HTmDXF0Xcbl9hdjlkUZpCKVUgzbWdiAiYnPLd9s8y/rJ2I0pWcMo3zU5zQ1YiIi2OnVbwi7c/wrce3swp35QWQ4mIMuaVd87i1OBlTvmmeeOvLUSkO1VV8Xdv9OPa2ATqv+RgING8caRERLoan4jhP3e8j8851uBzDk6CIm0YSkSkm8vXx/HUKyfw7z57B5uq0oIwlIhIF+eGr+H/7ezB110bsGbVMrPLoSzFUCKiRXvv7DBe6D6DJx/ahOKlhWaXQ1ls0RMdDh8+jG3btulRiy68Xi9sNlvSY2BgwOyyiHLWaycv4NV3zuFbDzOQaPEWHUpOpxMul0uPWnTh9XqhqmrSY926dWaXRZRzVFXFT39zCmeHruHPazagoIAz7GjxFn36rry8HE888YQetRBRlhibiOHH0vv4F3dxhh3pS9NI6fHHH89UHUSUJaKXR/FXLx5HXdVtDCTSnaZQev755/HMM89geHg4U/UQkYX1nB/Bjzvex55tG3Fn2Uqzy6EcpOn0XVtbG7Zv346Ojg7IsgyHw4Hq6upM1UZEFvJPJ87jrdND+PYjm1HI60eUIZpGStu3bwcA1NTUYNeuXSgvL8e2bdtw8ODBjBRHROZTVRXPvtGPCyPX8XXXBgYSZZSmUDp27Fjiv7t370ZlZSVUVYXT6cxEbURkstHxGJ565QTuvvkG1FXdbnY5lAc0nb6rq6uDIAiQZRlNTU3o6+tDSUlJpmojIhNdvHQd/7njfS7KR4bSFEqqquLAgQOoqanJVD1EZAG/PaPgfx0bQOOD92DVMjZ+IeNoXg59aiANDQ1BkiQ4HA7ce++9etdGRCb4x6Mf4uKl6/jWw5u45AQZbkETHYaHh9Hf349oNIrKykq0tLRkpLiFYJshooUZm4jhR4GTKCkuws4vigwkMoWmUOro6IDdbkd5eTkqKyvhdDpRWVmJ2traTNWnGdsMEWl38dJ1/NWLx/GvKm7Fl++5yexyKI9pOn0nSRIikQiAyYCqqanB0NAQQqFQRoojosz73Zkh/MPRD3n9iCxB07/AqqqqxP/39fUBAGffEWWxqdeP2FCVrEDT6TtBEHDkyBEAQDQaxVtvvQUACIfD+ldGRBkz/foRA4msQlMo2e12/OAHP0B/fz/q6+uxfft2lJWVobe3N1P1EZHOLozw+hFZl6bTdxUVFQgGg4nnPT09OHr0KCoqKnQvjIj019UXgfTuOXgevAcref2ILEjTSOnIkSOJVkNxDCQi61NVFT/5VR96L1xC01cYSGRdmkLpueeeS7k915eyeP/cCFRVNbsMogUZujqG7714HBV3lOLfbr2D9x+RpWkKpR07dkAUxRnbW1tbdSvIin7+u7MYjzGUKPu8/eEQfhQ4iT+v3oDP3C6YXQ5RWprG8IFAAHv37gWARGdwVVVx9OhR7NmzR//qLEKFiomYiqJCsyshmr/24GkoV8bwnUc2c3YdZQ1NoRQMBtHY2Ai73Z60XVEUPWtaFK/Xi3379iVtW7t27aKOqaqTD6JscHV0Aj/ueB9f3LCGy01Q1llUQ9Y4QRD0qmfRvF4vvF5v0rapN/0uhApggqlEWeD9cyN49jen8Kf334VbSpabXQ6RZppCKd7FYbqamprErLyc7BauqogxlMjCVFWFP3QGg5dH8Z1HNmNJoabLxUSWMe9/uUNDQygvL0/5tba2NuzcuRMtLS3o7OzUrTirUAHEONGBLOrS9XHs/8V7uLW0GLvvczCQKKvN+19vSUkJZFlGVVUVNmzYkGg3BEzOvuvs7MShQ4dysuWQqgLMJLKitz8cwsFXTqD+SyI+71hjdjlEi6bpV6re3l74fD68+uqr6OrqStyfJMsyVq9eDSA3G7TGZ98RWYWqqvgfb55CV18E33lkM9asWmZ2SUS60BRKW7duRU1NDcrLy3HgwIGklkNxuXhj3uTsO4YSWcOFkevY97PjuOeW1fiTL5RzujflFE0THbq6ulBZWYnS0lK0trbC7XZjaGgo6QN7tunhkiRBlmX09vbC4XCgvr5+xj7hcBiyLEMQBCiKArfbre1PkyGcfUdW0fneOXT3R7Fn20aufUQ5SdO/6oaGBtTU1CAajaK+vh6BQAChUAg+nw8HDx6EKIopT9/JsgwAiSByOBwQRREulytpv/3796O9vR0A4PF44HQ6U3aQMBqvKZHZro5O4L/+Uw82r10Nz4P3mF0OUcZoCqXy8nL09PQkbdu1axcA4OjRowgGg4nnU4XDYTz//POJEHK5XAiHw0mhJElS0k25DocDkiSlHFEZTYXK2Xdkmrc/HII/dAa773Pw3iPKebqN/ysqKmbtGO52u5MCSJZl1NbWJu0TP20XZ7fbEQgE9CpvcVTwPiUy3ERMxU9+3Y+iQhtbBVHeMOykdDxwFEVBJBKZcb1IURSUlZUlbYtEIkaVNycV4Ow7MtSpwcv4b6/34d9svQOb1q42uxwiwxh+l53H40FHR8eM7YIgYHBwMGnb9B57cV6vFzabbd6PgYGBRdWsqiqvKZEhYjEVP/3NKfz8d2fxzYc3M5Ao7xgaSs3NzfB4PCl75YmimDRzLxKJoLKyMuVxvF4vVFWd92PdunWLqptTwskIHwxewb6fvYPP3Cbg8fsdWLqEnRko/xj2r16SJLjd7sRsOkmSAHwyhdzlciWdrguFQjNm55mFU8Ipk+Kjoxd/N4AnH96E378t925AJ5ovXUIp3cqz4XAYdXV1iXucSktLE1+rq6tLtCZqamqC3++HJEmora21xHRw4OMp4TGzq6Bc9MHgFXzvxeP49G0l+NP778KyJVy0i/KbLhMdPB4PDh06NOvXnU4notFoyq9NnWHndDoTiwdaiQp2CSd9xWIq/mf3Bxi6Ooamh+5hGBF9bM5QKigomHfboLlCKdupnBJOOjoduYK/fb0PjzpvxadvE8wuh8hS5gwlt9uNF154IfG8o6MDoigmLWERDodnHQXlEk4Jp8WKxVQ8130a0SujHB0RzWLOa0ptbW1Jz1OtqeR0OnOyCetUnBJOi3Xy3Ai+9+JxfGrdanzty7x2RDSbOUdK0/vYzXYz62xNWHOFCp6+o4W5MjqO//Z6H1YXF+FbD2/iAnxEaWh6h/T09CSWPY87duwYuru79azJciZn3zGUSJvA8XP4G+l9uCtvxx99bj0DiWgeNM2+O3DgAB544AH09fVBEATIsgxRFFN2aMglKlTep0TzdiZ6BT/5VT++sGENnnxok9nlEGUVzVPCX3311URHcFEUUVNTk4m6Fszr9WLfvn1J29auXbuoY052dFjUISgPjE3E8HdvnMLoeAx7tm3E8iJeNyLSSvP5hOHhYYRCIdTW1qKmpgadnZ2ZqGvBUrUgWnSbIXD2Hc3tTXkQ33/pXdy/8UY8fr+DgUS0QJpCqaOjA9XV1QiFQomF+8rLyy0XTHrjfUo0m8FL1+F7+T2cHb6G7351Mxw3rjK7JKKspun0XSAQQDAYBIDEdaTy8vJEm6DcxY4OlCwWU+EPncGZ6BXsvs+BkuIis0siygmaQmnr1q0pt+f+fUrsfUefOHZawT8e/RD/uuJWPLbldrPLIcopmk7fdXV1YWRkBMAnQdTf34+uri79K7MQnr4jADg7dA0HfvEe5AuX8J1HNuMztwtml0SUczSNlJqamlBRUZHo8q0oCgRByIsp4Qyl/HX5+jh++ptTGI+p+LPqu7BymWELNhPlHU3vrpKSEvT09ODw4cOJe5S2b9+eqdosY3KkZHYVZLTR8RheCJ7GgHIV/9dn78StQrHZJRHlPE2h1N/fD0VRsH37dgwNDaGtrQ0HDx7Enj17MlWfJXBKeH6ZiKn43299iOMDw3is6nZsuPkGs0siyhuarikdOHAgMRXc5XKht7cXNTU1OHjwYEaKswpeU8oPqqqi491z+P5L72J92Up88+HNDCQig2kaKdXW1uLRRx9FX18fent7Ez3v4kGVq3hNKfd190fw8999hOp7bsK3H9mU8zNKiaxK00gpPsFBkiS4XK7Ediu9gb1eL2w2W9JjYGBgcQfllPCc9e5Hw/irF4/j/PB1fPvhzfjihhst9e+ZKN9oCqVQKITDhw/D5/Nh9+7dACZvop1tSQszZKLNkM1mY0PWHPPB4BXs/8W7+O0ZBXu/cg8e/vRaFBQwjIjMpun03RNPPIG2tja0tLSguroaHR0dCIfDiRFUriosmLzeQNnv/Mg1PPvGKQgrluIbrrvZo47IYjTfcLFr1y4MDw/j2LFj2LJli+W6hGdCgc2GCZ6+y2rnhq/hua7TWFJow84vimwLRGRRmkNp9+7daG1thSAIGBoagtvtxvPPP5+J2iyjoMDGiQ5Zqv/iZbSHTmP18iL88efXo2QFw4jIyjRdU3rqqadQW1uLWCyGSCSCiYkJPPbYYzk/JbzAxinh2eadgSHs//m7eO39C/jal+9Cw30OBhJRFtA0UkrVwWH79u1oa2vTtSirKbTZuBx6lgh/EMUvfvcR7r75BvzlAxuxdAmXICfKJppCabapsmVlZboUY1U2m41thixsIqYicPwsgv1R3HuHgL1f2YRCzqQjykqafo3s7e3F8PBw0rZ86BJeYOM1JSu6fH0cz77Rj+aX38OaVcvwzYc34ZFPr2MgEWUxTSOl+vp6VFdXw2azwW63IxKJQFEUhEKhTNVnCYUFvKZkJacjV3Ak/CFiqopHnbfizrKVZpdERDrR3CU8GAzC7/ejr68vb7qEc0q4+VRVRVdfBJ0nzmNdSTH+5AvrccNyTlwgyjULWhjG7XYnPX/mmWewc+dOXQpaLK/Xi3379iVtW7t27aKOySnh5jk3fA3/+9gABi+PourOUjRuu4en54hymE2do1XBtm3b0h5AVVWEQiEMDg7qWpieqqqqEAwGF/z6p155D8uXFOLPajZgIqbiyug4f0vPoLGJGDrePY+jH0Rx8+rl+Opn1uHGG5aZXRYRLcJ8P4fnHCmpqgqPx5P2IC0tLfOvLAsVTpl99w9HP8SV0XH80efWm1pTLpIvXMJLv/0IV8cmULPpZuz9yj1sjkqUZ+YMJZ/Ph4qKirQHsdvtuhVkSVNm37394RDKVi41uaDccX7kGn721kc4P3INjjWr8B/+xXqs5iiUKG/NGUqpAunYsWO49957k57PJ7iy2dTf1UuKizDOm5YWZejKGALvnkPP+UtYs2opvvqZdbh59XKzyyIiC9A00eHw4cPweDzw+/2JYFJVFZ2dnaiurs5EfZQjhq6O4ZV3zqLv4mWUFBfBtelmuCtvM7ssIrIYTaEkCAJ6enqStlVUVODIkSO6FmVVsZgKm21yeXRK78roOKR3z+P4wDBWFy/Btk/dgseqbje7LCKyME2hNDQ0lHK7lRb5y6TxmIqiwgKMjvOmpdlcH5/AL09cwNHTClYUFaJm08346qfXcsICEc2LplDq6uqCKIpJ15Q6OzsRCoUsc59SJo3HYrxHJoWxiRje6B3Em32DWFJQgPs23ojGbRsZRESkmaZQampqQk1NTaKbgyzLEEURHR0dmarPUsZjKpYwlABMTlb4+dsf4Uz0CpYUFOAPRDv+onYjQ5uIFmVBbYYkSUoEUz6sPBs3MTEZSvl6SenU4GUEjp+DcmUMxUsL8dVPr8MdZSvMLouIcsiC2gy5XK6k5/39/Vi/fr0e9SxaJtoMxY3FYigsLMCyJQW4Pj6BZUsKdTmuVV0ZHcevegbxuzMKJlQV64RibHfehlLep0VEGbKgUJq+fIXP58OhQ4d0KWixvF4vvF5v0raqqipdjj0RU1FUYMPq4iIMXR3DTTfkViipqoqe85fwi7fPYnQ8huVFBficYw2+7rqbp+WIyBCaQqmtrQ0NDQ2w2WxQVTXpv1YJpUwan1BRWGDDimVLMHx1DDfdkN03fKqqit4Ll/Fm3yA+Uq7BZgPWl63Eri+KKF6aW4FLRNlBUyj19vYiGo2ipKQkafvevXt1Lcqq4lPC15UsR//FK7jrphvMLkmz8yPX0PHueXwYvQoAEG9cifs33oRbhWKTKyMi0hhKtbW1MwIJmJyVlw/GJyanhFfeWYrvv/QuXJtvNrukOQ1fG0OXHMHR01EU2GwYj6m4cdUyuDbdzAkKRGRJmkLJZrOlnNTQ3t6eJ/cpqSgqtMFms2GtUIzBS9dRtsoaSypELo/itZMXIF+8jOKiQly6PobiokJ8/q41+IbrbgDAksICk6skIpqbplB6+umncfToUQCAKIoAgMHBQfT19eVFKE3EVBQWTH6w37x6GS5eGjUslMYmYlhSMBmI/Rcvoz10GrYprWJLiotw38Yb8S/vXYfhq+MoWcFO20SUfTSFkqIoOHDgAARBSGxTVRXNzc1pXydJErq7u+Hz+VLu09DQgLq6OlRVVcHj8VhyjaZ4MABA2cplGLx8HcAn15V+3XMRv3dbyYKWXlBVFdfHJ4+/pLAAoVNR/PLkBaiqirEJFdfHJ7C8qBATMRV32Ffga1++CyuWpv7xMZCIKFtpCqXZ1lcqKyub83Xx1QYVRZl1n0gkkgglKwYSMDlSWlL4cSitWooTZ0cwdHUM7cHT+PRtAn4jD6L34mX8+8/eOeO1o+MxFBXaMDah4uroBH7VexErly3BHfYVCJ2K4te9F1FcVAj7yqW4MjqBT99Wgq/XbEBhgS0v7okiIgI0hlJFRQWGh4fxwgsvwOVyYf369fNatsLlciEcDs+5z44dO9De3q6lHMONfTwlHADsK5cicnkUP/lVP84OX8PItXF8rfou/PDVk7g+PoFn3ziFbZ+6Bb88eQHXxiYgf7xkw5Xr4yheugSOG1di5bIlCJ2KoqjQhv2P/j5UFVheNDN8GEhElC80hVJHRwc8Hg+2bNkCURSxfv16lJeX67KekizL8Pv9ACZHTfX19Ys6XiZMxFQsK5q8plS6YimiV0YRi6n4rGjHL09cwLIlhfhXFbei4dkQvl6zAa/3XETV+lIsLSzA/122kjegEhGloSmUAoFA4lRcvAlreXl52lHQfNTX1yeuVdXV1cHlciUmU1jFWCyGFQWTo5bJ02oxFBcV4r67b0xc39m0djX++x9vgc1mQ8UdpWaWS0SUdTTNEd66dWvK7XosUTB18oQoiolRUyperxc2m23ej4GBgUXXB0w2ZC0q+OSv7HTkCu4sWwFhxVLUTrlniUs2EBEtjKZQ6urqwsjICIBPPnj7+/vR1dW1qCIkSUJDQ8O89/d6vVBVdd6PdevWLao+AFAxeZ/S1FNwg5dGIa5ZtehjExHRJM3rKVVUVKC0dPK0lKIoEARhwespxV8viiLq6uoS28PhsCVn4I3HJmfQxdlXLoV440oTKyIiyi2a11Pq6enB4cOHEwv8bd++Pe3rJElCe3s7JEmC3++H2+0GMHntyOfzwel0IhwOw+/3Q5ZlNDQ0WO56kg3xm2c/CaUdW27HymULarROREQp2FRVzfk166qqqhITNBbiR4GTAIDb7Suwdb2dfeOIiDSa7+cwm6FpMBGLJW6eJSIi/TGUNBj7eDl0IiLKDIaSBpNthvhXRkSUKfyE1WD6lHAiItKXLqHU39+vx2Esb3wieUo4ERHpa0GhNDw8nPSYbTmKXMOREhFRZmkKpba2NhQUFKC0tBSCICT+29ramqn6NEvVgkivNkPj09oMERGRvjR9wvb29iIajWJiYgKxWCzx3yeeeCJT9WmWqgWRHm2GACCmqijgSImIKGM0hVJtbS1KSkpmbG9qatKtICIiyl+aQslms6Wc1GD1xfmIiCg7aGrc9vTTT+Po0aMAkOhNNzg4iL6+PuzcuVP/6oiIKK9oCiVFUXDgwIGktY9UVUVzc7PedRERUR7SFEo+nw8VFRUztpeVlelWkFWpmOwUTkREmaPpmlKqQAKAUCikSzFWxkAiIsq8OUdKjz/+OOrq6lBdXQ0A2LZt24x9VFVFKBTK+WtKBTYbxmMxs8sgIsppc4bS9KWWVFWFx+OZsZ8VV4nVW0nxEvRdvGx2GUREOW3OUHr66aeTns92Tclut+tblQUJK5ZCuaqYXQYRUU7T5ZrSbNvNkKk2QyUriqBcGdOhQiIimk3ONXLLVJshobgIQ1cZSkREmZRzoZQpwoqlUNPvRkREi8BQmiehuAhFbMZKRJRRDKV5Wl1chCIuhU5ElFGaPmWPHDmCY8eOZagUaysssMG+cqnZZRAR5TRNofTcc8+l3D48PKxLMVZXtoqhRESUSZpCaceOHYnu4FNZaeXZTPryPTeZXQIRUU7T1JA1EAhg7969AACn0wlgssvD0aNHsWfPHv2rs5gvb2QoERFlkqZQCgaDaGxsnNHBQVEUPWsiIqI8pXnpipqamhnbp66vREREtFCarinFA2l4eBjHjh1LTHBIFVRmyVSbISIiyjzNN97s3r0bgiCguroapaWl2LFjRybqWrBMtRkiIqLM0xRKTz31FGpraxGLxRCJRDAxMYHHHnsMBw8ezFR9RESURzSFkiiK2L59e9K27du3o6SkRNeiiIgoP2kKJZstde+3srIyXYohIqL8pimUent7Z3Rv6O/vR1dXl65FERFRftI0Jby+vh7V1dWw2Wyw2+2IRCJQFAWhUChT9RERUR7RFEolJSUIBoM4fPgwZFlOeY2JiIhooTSF0pEjRxhERESUMewSTkRElsEu4UREZBnsEk5ERJaRc13CvV4v9u3bl7Rt7dq1JlVDRERa5FyXcK/XC6/Xm7StqqrKnGKIiEgTTdeUhoaGcOzYsRnbrdQlnIiIshdn3xERkWVw9h0REVkGZ98REZFlGDL7TlEUSJKE7u5u+Hy+lPuEw2HIsgxBEKAoCtxut5bSiIgoBxgy+y4YDAKYO7z279+P9vZ2AIDH44HT6Ux5qpCIiHKXpmtKNTU1GB4exjPPPIP+/n4AQGdnZ9rZdy6Xa86AkSQpafTlcDggSZKW0oiIKAdoCqWOjg5UV1cjFApBlmUAQHl5OTo7OxdVRPy0XZzdbudyGEREeUjzRIf4qbiOjg4Ak6EUDocXVYSiKDNWr41EIos6JhERZR9NI6WtW7em3D7bMunzJQgCBgcHk7ZNn0wxldfrhc1mm/djYGBgUfUREZExNIVSV1cXRkZGAHwSRHoshy6KYtIkiEgkgsrKyln393q9UFV13o9169Ytqj4iIjKGptN3TU1NqKioQGlpKYDJ026CICRO5WkVf73L5UJLS0tieygUgsfjWdAxiYgoe2leDr2np0fzcuiSJKG9vR2SJMHv9yfuQaqrq4PP54PT6URTUxP8fj8EQUBtbS2ngxMR5SGbqqqq2UVkWlVVVWKCxkL8KHAS36i9W8eKiIjyy3w/hzVdU8pHeZDZRESWwVBKQ1WBRU4uJCKieWIopaECsIGpRERkBIZSGqqqcqRERGQQhlIakyMlIiIyAkMpDV5TIiIyTs6FUqoWRItpM6RCXXQbJSIimp+cDCU92wxxRjgRkXEWHUqHDx/Gtm3b9KjFsjhQIiIyxqJDyel0wuVy6VGLJakqp4QTERll0aFUXl6OJ554Qo9aLGnympLZVRAR5Qddrik988wzehzGkiZHSkREZIQ5u4TP51qRqqoIhULYuXOnbkVZiQpeUyIiMsqcoaSq6rzWNZq6FlKuUVWV15SIiAwyZyj5fD5UVFSkPchcS5dnO46UiIiMM+c1pVSBNDw8jGeeeQb9/f0AgM7OznkFV7bifUpERMbRNNGho6MD1dXVCIVCkGUZwOTsu87OzowUZwkq2NGBiMggmkIpEAggGAzi0KFDicXvysvLEY1GM1LcQmSkzZCO9RER0ew0hdLWrVtTbrfSSCITbYYs9McjIsppmkKpq6sLIyMjAD4Jov7+fnR1delfmUVw6QoiIuPMOftuuqamJlRUVKC0tBQAoCgKBEFAR0dHRoqzgslF/hhLRERG0BRKJSUl6OnpweHDhyHLMkRRxPbt2zNVmyVwSjgRkXE0hdKRI0fyIoimYpshIiLjaLqm9Nxzz6XcPjw8rEsxVqSCMx2IiIyiKZR27NgBURRnbG9tbdWtIMvhSImIyDCaTt8FAgHs3bsXwOQ6SsDkRICjR49iz549+ldnAbymRERkHE2hFAwG0djYOKPXnaIoetZkKVzkj4jIOJpCyefzoaamZsZ2QRD0qsdyuMgfEZFxNF1TShVIc203g+5thnhNiYjIMLqsPGslurcZAq8pEREZJedCSW9c5I+IyDgMpTRUNr8jIjIMQ2kemElERMZgKKWhcpE/IiLDMJTS4CJ/RETGYSilwUX+iIiMw1BKg1PCiYiMw1BKI8Yp4UREhmEopcHTd0RExmEopaWaXQARUd7IuVDKSO87DpWIiAyRk6Gkd++7AmYSEZEhci6U9Mb1lIiIjMNQSoPrKRERGYehlAbXUyIiMg5DKQ1OCSciMo5hoRQOh+H3+yFJEvx+f8p9GhoaIEkSFEVBQ0ODUaXNSQXXriAiMophobR//3643W64XC50d3dDluUZ+0QiEdTV1aGurg4ej8eo0ubEkRIRkXEMCSVJkmC32xPPHQ4HJEmasd+OHTsQjUYRCAQgiqIRpc0LM4mIyBhLjPgmsixDEITEc7vdjkAgkHK/+Km9SCSC+vp6I8qbE2+eJSIyjiGhpCgKysrKkrZFIpEZ+9XX1yfCq66uDi6Xy/QRE9dTIiIyjiGn7wRBwODgYNK2qafzpu4XJ4rirBMiUrUSmuux+DZDC345ERFpYEgoiaIIRVESzyORCCorK5P2kSRp3jPuUrUSmuux2DZDDCUiImMYEkoulyvpdF0oFILL5QKARFiJooi6urrEPuFwGG6324jy5qRyPSUiIsMYck0JAJqamuD3+yEIAmpraxPXiurq6uDz+eB0OhP3MsmyjIaGBtOvJwEfL1zBTCIiMoRhoeR0OuF0OmdsnzoLzwojo+nYZoiIyDhsM5SWyinhREQGYSilwZESEZFxGEppcPYdEZFxGEppcJE/IiLjMJTSUFUu8kdEZBSGUhpcuIKIyDg5F0qpWhAtts0QU4mIyBg5GUr6thliRwciIqPkXCjpjg1ZiYgMw1BKg2fviIiMw1BKg4v8EREZh6GUhgpOCSciMgpDKQ22GSIiMg5DKQ22GSIiMg5DKQ2VNyoRERmGoZQGR0pERMZhKKXDa0pERIbJuVDSvc0QF/kjIjJMToaSrm2GOFIiIjJMzoWS3lS2GSIiMgxDKY3JuXdMJSIiIzCU0uAif0RExmEopaGaXQARUR5hKKXBa0pERMZhKKXFRf6IiIzCUEqDIyUiIuMwlNIoWVGE4qJCs8sgIsoLS8wuwOo+71hjdglERHkj50ZKercZIiIi4+RkKOnZZoiIiIyTc6FERETZi6FERESWwVAiIiLLYCgREZFlMJSIiMgyGEpERGQZDCUiIrIMhhIREVkGQ4mIiCyDoURERJaRc6HE3ndERNkrJ0OJve+IiLJTzoUSERFlL5uqqqrZRWTamjVrsH79+sTzgYGBeY2e5rsfj2neMc3+/jwmj8ljzm/f/v5+XLx4Me1r8iKUprPZbJjPH3u++/GY5h3T7O/PY/KYPKY++8bx9B0REVkGQ4mIiCyDoURERJaRl6H03e9+17RjavnerFPffVmnvvuyTn33zec6p8rLiQ7ztZCLdGZgnfrLllpZp75Yp7440YGIiLIaQ4mIiCyDoURERJbBUCIiIstgKBERkWUwlOaQ6amPemGd+suWWlmnvlinvhZSJ6eEExGRZXCkRERElrHE7AKsKBwOQ5ZlCIIARVHgdrvNLiklSZIgyzJ6e3vhcDhQX19vdklzUhQFra2taGxsNLuUWbW2tkIURUv/3IHJf6PBYDDx3Co/e0VRIEkSuru74fP5Etut9p6arU6rvadmq3Pq163wnpqrTq3vKYZSCvv370d7ezsAwOPxwOl0QhRFk6tKJssygE8+jBwOB0RRhMvlMrOsOe3fv9/sEubU0NAAj8cDURQT/3U6nWaXlVIwGEz87P1+P8LhsCVqjQeloihJ2632nkpVpxXfU7P9fcZZ5T01W50LeU/x9N00kiTBbrcnnjscDkiSZGJFqYXDYbS0tCSeu1wuhMNhEyuaWzgchsPhMLuMWSmKAlmWEx+UPp/PEh/yqSiKkviAB4BIJAJBEMwraAqXyzUjbKz4nkpVpxXfU6nqjLPSeypVnQt9TzGUpomfYoiz2+0IhULmFTQLt9uNtra2xPOpP3wrikQilq5PkiQIggC/3w+/34/W1lazS5qVIAgQRREOhwN+vx92u93Sf7d8T2VGrr6nGErTKIqCsrKypG2RSMSkauYWf6MrioJIJGL6efrZSJJk6dOKwOQHkCzLcLvdcLvdCIVCpv82P5f4b527du2y7L/POL6n9JfL7ymG0jSCIGBwcDBp29RTD1bk8XjQ0dFhdhkpKYpi+b8/YPLnXlVVlXjucDgQCARMrGh28Yvb7e3tCIVCaGlpgd/vN7usWfE9pa9cf08xlKaJzxKJi0QiqKysNK+gNJqbm+HxeCxzTWG61tZWyLIMv9+PQCCAcDhsyRGIlU+DTPfCCy8kfksWRREdHR2WDVCA7ym95fp7iqE0jcvlSjq1EAqFLDtMliQJbrc78cO34j/MxsbGxPB9y5YtcDqdlvz7nP5z7+3tRW1trYkVzU4UxRmnv6xaK8D3lN5y/T3Fjg4pWO2eilTC4TBqamqStrW3t1vyHycwWe/+/fshyzJ8Pp8l64z/xhn/0Df7HpW5tLa2wm63J970VqlVkiS0t7dDkiT4fL7Ee8dq76lUdVrxPTXb3ydgrffUXD93re8phhIREVkGT98REZFlMJSIiMgyGEpERGQZDCUiIrIMhhIREVkGQ4mIiCyDoUSUY5qbm7P6+JTfGEpEOlIUBXV1dbDZbEnbGxoa0NDQkPHv39ramrhxUZIkVFZWoq6uTtfv4Xa7GUyUMQwlIh0JgpC01lFcXV2d7uEwXXzF1HiLHJfLlZEgFEURg4ODiUXxiPTEUCIygMvlyngbmJaWlhkhlKlu0k1NTSmX5yZaLIYSUY4Ih8OGdTsXBIEjJcqIJWYXQGQFkiTB4/GgqqoKtbW1eP7557Fjx45Eo85IJAJFUdDd3Y3a2toZox6PxwOHwwG73T5jyYNwOAyPxwMAiaUGpj5XFAUejwetra2IRqOJ1zc3N8PpdEJRFAQCAXg8nllDZz6rpCqKgvLyclRVVcHj8cButyeWaIiPsMLhMARBSCwFHolEEms2TSeKIsLhsGWXjacspRKRqqqq2tLSooqiqEajUTUUCqmhUEhVVVV1Op1qe3t7Yj9BEJJe53K5Evuqqqr29vaq099agUBAdblcsz5XVVUFoEaj0UQtU79ne3t70veYLhAIqI2NjTO2t7e3q263O1FXS0vLjK+Loqj29vaqqqqq0WhUBaAGAoGkP9/UWuJ8Pl/K7USLwZES0cfi118EQUj67b+9vX3GKERRFAiCkFiSYer+epxCE0URDQ0NiEQicLlcaZd6kGUZDodj1q/HlxBobGxM2h4flcVrnv48vi3V8uU8hUeZwGtKRFOkOhVlt9vR3NyM1tbWxKJv8Q/pYDCYkRVKXS4XWlpaEAgEUFlZicrKyqTVW1OZ7evhcBjBYBD79+9PuU+qEJ06QWKuyRJWXZ2VshdDiSiNyspKuFwu1NfXJ11LUhQFoihmZLQgSRJcLhfa29sRjUZRVVWF1tbWWfe32+0YHBxM+bV47U1NTdi1a5duNSqKkrHZfZS/GEpEc4hf7I+PoKaONOLBEb/gP/U16UwPs+mviZ9ui0t3v5EoimlHUo2NjQiHw/D7/UnbU52am8/Xp94TRaQXhhIRJgOmpaUF4XAYzc3NiQ94p9OJxx57DM3NzZAkCcFgEG1tbfD5fIkP5I6ODjz//PPw+/2QJCnxId7Q0ABFURAOh9HS0oJgMJgY7YiiCLfbjdbW1kRICIKAXbt2JZYNl2UZfr8/cdzp14OmcjqdM0Zs8T+TJElJ36Ourg7Nzc0Ih8Pw+XyQZRmtra2JWYDA5GzC+PZ4/VNDEsCMa2lEeuBy6EQ5oqGhAT6fz5DrPPEASzVVnGgxOFIiyhEejwf79+835Hu1trYmRlVEemIoEeUIURRRVlaW8Wna8dOLvJ5EmcBQIsohjY2NMyYy6M3v96O+vj6j34PyF68pERGRZXCkRERElsFQIiIiy2AoERGRZTCUiIjIMhhKRERkGQwlIiKyjP8fh31cSebA9KAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if plots:\n",
    "    pl.canvas_flag=False\n",
    "    pl.canvas()\n",
    "    plot.plot(radial_grid,rel_diff)\n",
    "    pl.xtitle('radius (km)')\n",
    "    pl.ytitle('rel. error in enclosed grav. mass')\n",
    "    plot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f17c764",
   "metadata": {},
   "source": [
    "For testing using ``pytest``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bf8fffcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_fun():\n",
    "    assert numpy.allclose(b.rad_from_gm(1.4),ts.rad,rtol=1.0e-9,atol=0)\n",
    "    for i in range(0,len(rel_diff)):\n",
    "        assert numpy.allclose(rel_diff[i],0.0,atol=5.0e-11)\n",
    "    return"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}