{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3fc8183b",
   "metadata": {},
   "source": [
    "# Tolman-Oppenheimer-Volkov example for O$_2$sclpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1de243df",
   "metadata": {},
   "source": [
    "See the O$_2$sclpy documentation at https://awsteiner.org/code/o2sclpy for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8c46b427",
   "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": "cd8c039f",
   "metadata": {},
   "source": [
    "Link the O$_2$scl library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b1f114c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "link=o2sclpy.linker()\n",
    "link.link_o2scl()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38301d17",
   "metadata": {},
   "source": [
    "This is the APR EOS as computed in ``o2scl/examples/ex_eos_had_apr.cpp``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0778e337",
   "metadata": {},
   "outputs": [],
   "source": [
    "eos=[[4.500000e-02,2.156413e-01,6.998033e-04],\n",
    "     [8.500000e-02,4.084332e-01,2.456061e-03],\n",
    "     [1.250000e-01,6.025862e-01,7.007571e-03],\n",
    "     [1.650000e-01,7.987918e-01,1.572468e-02],\n",
    "     [2.040000e-01,9.927169e-01,2.946785e-02],\n",
    "     [2.440000e-01,1.193221e+00,3.249203e-02],\n",
    "     [2.840000e-01,1.395525e+00,5.124024e-02],\n",
    "     [3.240000e-01,1.600975e+00,7.771385e-02],\n",
    "     [3.640000e-01,1.810261e+00,1.139703e-01],\n",
    "     [4.040000e-01,2.024147e+00,1.620335e-01],\n",
    "     [4.440000e-01,2.243442e+00,2.238295e-01],\n",
    "     [4.840000e-01,2.468981e+00,3.012283e-01],\n",
    "     [5.240000e-01,2.701614e+00,3.960900e-01],\n",
    "     [5.640000e-01,2.942202e+00,5.102769e-01],\n",
    "     [6.040000e-01,3.191619e+00,6.456452e-01],\n",
    "     [6.440000e-01,3.450743e+00,8.040340e-01],\n",
    "     [6.840000e-01,3.720458e+00,9.872584e-01],\n",
    "     [7.240000e-01,4.001649e+00,1.197107e+00],\n",
    "     [7.640000e-01,4.295201e+00,1.435340e+00],\n",
    "     [8.040000e-01,4.602000e+00,1.703690e+00],\n",
    "     [8.440000e-01,4.922926e+00,2.003863e+00],\n",
    "     [8.840000e-01,5.258860e+00,2.337541e+00],\n",
    "     [9.240000e-01,5.610674e+00,2.706378e+00],\n",
    "     [9.640000e-01,5.979240e+00,3.112009e+00],\n",
    "     [1.004000e+00,6.365422e+00,3.556046e+00],\n",
    "     [1.044000e+00,6.770081e+00,4.040080e+00],\n",
    "     [1.084000e+00,7.194071e+00,4.565686e+00],\n",
    "     [1.124000e+00,7.638242e+00,5.134417e+00],\n",
    "     [1.164000e+00,8.103438e+00,5.747813e+00],\n",
    "     [1.204000e+00,8.590497e+00,6.407393e+00],\n",
    "     [1.244000e+00,9.100252e+00,7.114660e+00],\n",
    "     [1.284000e+00,9.633531e+00,7.871098e+00],\n",
    "     [1.324000e+00,1.019116e+01,8.678169e+00],\n",
    "     [1.364000e+00,1.077394e+01,9.537317e+00],\n",
    "     [1.404000e+00,1.138270e+01,1.044996e+01],\n",
    "     [1.444000e+00,1.201824e+01,1.141747e+01],\n",
    "     [1.484000e+00,1.268134e+01,1.244122e+01],\n",
    "     [1.524000e+00,1.337281e+01,1.352252e+01],\n",
    "     [1.564000e+00,1.409343e+01,1.466262e+01],\n",
    "     [1.604000e+00,1.484397e+01,1.586275e+01],\n",
    "     [1.644000e+00,1.562519e+01,1.712404e+01],\n",
    "     [1.684000e+00,1.643786e+01,1.844757e+01],\n",
    "     [1.724000e+00,1.728270e+01,1.983432e+01],\n",
    "     [1.764000e+00,1.816046e+01,2.128515e+01],\n",
    "     [1.804000e+00,1.907185e+01,2.280082e+01],\n",
    "     [1.844000e+00,2.001758e+01,2.438194e+01],\n",
    "     [1.884000e+00,2.099830e+01,2.602899e+01],\n",
    "     [1.924000e+00,2.201470e+01,2.774223e+01],\n",
    "     [1.964000e+00,2.306741e+01,2.952175e+01]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63a59491",
   "metadata": {},
   "source": [
    "Create a table, the columns, and the units for each column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "43669964",
   "metadata": {},
   "outputs": [],
   "source": [
    "tab=o2sclpy.table_units(link)\n",
    "tab.line_of_names(\"nb ed pr\")\n",
    "tab.line_of_units(\"1/fm^3 1/fm^4 1/fm^4\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2731bf36",
   "metadata": {},
   "source": [
    "Fill the table with the Python data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "578f4521",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(0,len(eos)):\n",
    "    tab.line_of_data([eos[i][j] for j in range(0,3)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a66c63c",
   "metadata": {},
   "source": [
    "Create a new column in the table for the chemical\n",
    "potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fc86f1e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "tab.function_column('(ed+pr)/nb','mu')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e08d4c99",
   "metadata": {},
   "source": [
    "Create the object which interpolates the EOS for the TOV\n",
    "solver. Use the default crust EOS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ce77424f",
   "metadata": {},
   "outputs": [],
   "source": [
    "eti=o2sclpy.eos_tov_interp(link)\n",
    "eti.default_low_dens_eos()\n",
    "eti.read_table(tab,'ed','pr','nb')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52852053",
   "metadata": {},
   "source": [
    "Specify the EOS and determine the M-R curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "dc26b3ab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mass versus radius mode.\n",
      "Central P: 7.0000e-07 (Msun/km^3), M: 1.0722e-01 (Msun), R: 5.7745e+01 (km)\n",
      "Central P: 7.7000e-07 (Msun/km^3), M: 1.1205e-01 (Msun), R: 4.6451e+01 (km)\n",
      "Central P: 8.4700e-07 (Msun/km^3), M: 1.1754e-01 (Msun), R: 3.8926e+01 (km)\n",
      "Central P: 9.3170e-07 (Msun/km^3), M: 1.2361e-01 (Msun), R: 3.3646e+01 (km)\n",
      "Central P: 1.0249e-06 (Msun/km^3), M: 1.3020e-01 (Msun), R: 2.9792e+01 (km)\n",
      "Central P: 1.1274e-06 (Msun/km^3), M: 1.3722e-01 (Msun), R: 2.6891e+01 (km)\n",
      "Central P: 1.2401e-06 (Msun/km^3), M: 1.4460e-01 (Msun), R: 2.4649e+01 (km)\n",
      "Central P: 1.3641e-06 (Msun/km^3), M: 1.5264e-01 (Msun), R: 2.2817e+01 (km)\n",
      "Central P: 1.5005e-06 (Msun/km^3), M: 1.6165e-01 (Msun), R: 2.1265e+01 (km)\n",
      "Central P: 1.6506e-06 (Msun/km^3), M: 1.7164e-01 (Msun), R: 1.9961e+01 (km)\n",
      "Central P: 1.8156e-06 (Msun/km^3), M: 1.8259e-01 (Msun), R: 1.8868e+01 (km)\n",
      "Central P: 1.9972e-06 (Msun/km^3), M: 1.9444e-01 (Msun), R: 1.7954e+01 (km)\n",
      "Central P: 2.1969e-06 (Msun/km^3), M: 2.0714e-01 (Msun), R: 1.7189e+01 (km)\n",
      "Central P: 2.4166e-06 (Msun/km^3), M: 2.2062e-01 (Msun), R: 1.6546e+01 (km)\n",
      "Central P: 2.6582e-06 (Msun/km^3), M: 2.3480e-01 (Msun), R: 1.6004e+01 (km)\n",
      "Central P: 2.9241e-06 (Msun/km^3), M: 2.4970e-01 (Msun), R: 1.5542e+01 (km)\n",
      "Central P: 3.2165e-06 (Msun/km^3), M: 2.6593e-01 (Msun), R: 1.5134e+01 (km)\n",
      "Central P: 3.5381e-06 (Msun/km^3), M: 2.8354e-01 (Msun), R: 1.4777e+01 (km)\n",
      "Central P: 3.8919e-06 (Msun/km^3), M: 3.0246e-01 (Msun), R: 1.4467e+01 (km)\n",
      "Central P: 4.2811e-06 (Msun/km^3), M: 3.2260e-01 (Msun), R: 1.4201e+01 (km)\n",
      "Central P: 4.7092e-06 (Msun/km^3), M: 3.4383e-01 (Msun), R: 1.3972e+01 (km)\n",
      "Central P: 5.1802e-06 (Msun/km^3), M: 3.6604e-01 (Msun), R: 1.3776e+01 (km)\n",
      "Central P: 5.6982e-06 (Msun/km^3), M: 3.8428e-01 (Msun), R: 1.3630e+01 (km)\n",
      "Central P: 6.2680e-06 (Msun/km^3), M: 3.9798e-01 (Msun), R: 1.3505e+01 (km)\n",
      "Central P: 6.8948e-06 (Msun/km^3), M: 4.1319e-01 (Msun), R: 1.3369e+01 (km)\n",
      "Central P: 7.5843e-06 (Msun/km^3), M: 4.2997e-01 (Msun), R: 1.3230e+01 (km)\n",
      "Central P: 8.3427e-06 (Msun/km^3), M: 4.4825e-01 (Msun), R: 1.3091e+01 (km)\n",
      "Central P: 9.1770e-06 (Msun/km^3), M: 4.6789e-01 (Msun), R: 1.2956e+01 (km)\n",
      "Central P: 1.0095e-05 (Msun/km^3), M: 4.8929e-01 (Msun), R: 1.2826e+01 (km)\n",
      "Central P: 1.1104e-05 (Msun/km^3), M: 5.1268e-01 (Msun), R: 1.2702e+01 (km)\n",
      "Central P: 1.2215e-05 (Msun/km^3), M: 5.3792e-01 (Msun), R: 1.2586e+01 (km)\n",
      "Central P: 1.3436e-05 (Msun/km^3), M: 5.6485e-01 (Msun), R: 1.2479e+01 (km)\n",
      "Central P: 1.4780e-05 (Msun/km^3), M: 5.9351e-01 (Msun), R: 1.2380e+01 (km)\n",
      "Central P: 1.6258e-05 (Msun/km^3), M: 6.2438e-01 (Msun), R: 1.2289e+01 (km)\n",
      "Central P: 1.7883e-05 (Msun/km^3), M: 6.5730e-01 (Msun), R: 1.2206e+01 (km)\n",
      "Central P: 1.9672e-05 (Msun/km^3), M: 6.9200e-01 (Msun), R: 1.2133e+01 (km)\n",
      "Central P: 2.1639e-05 (Msun/km^3), M: 7.2847e-01 (Msun), R: 1.2066e+01 (km)\n",
      "Central P: 2.3803e-05 (Msun/km^3), M: 7.6716e-01 (Msun), R: 1.2007e+01 (km)\n",
      "Central P: 2.6183e-05 (Msun/km^3), M: 8.0780e-01 (Msun), R: 1.1954e+01 (km)\n",
      "Central P: 2.8801e-05 (Msun/km^3), M: 8.5003e-01 (Msun), R: 1.1908e+01 (km)\n",
      "Central P: 3.1681e-05 (Msun/km^3), M: 8.9404e-01 (Msun), R: 1.1866e+01 (km)\n",
      "Central P: 3.4850e-05 (Msun/km^3), M: 9.3994e-01 (Msun), R: 1.1829e+01 (km)\n",
      "Central P: 3.8335e-05 (Msun/km^3), M: 9.8731e-01 (Msun), R: 1.1796e+01 (km)\n",
      "Central P: 4.2168e-05 (Msun/km^3), M: 1.0359e+00 (Msun), R: 1.1765e+01 (km)\n",
      "Central P: 4.6385e-05 (Msun/km^3), M: 1.0860e+00 (Msun), R: 1.1737e+01 (km)\n",
      "Central P: 5.1023e-05 (Msun/km^3), M: 1.1372e+00 (Msun), R: 1.1711e+01 (km)\n",
      "Central P: 5.6126e-05 (Msun/km^3), M: 1.1891e+00 (Msun), R: 1.1686e+01 (km)\n",
      "Central P: 6.1738e-05 (Msun/km^3), M: 1.2420e+00 (Msun), R: 1.1661e+01 (km)\n",
      "Central P: 6.7912e-05 (Msun/km^3), M: 1.2953e+00 (Msun), R: 1.1637e+01 (km)\n",
      "Central P: 7.4703e-05 (Msun/km^3), M: 1.3487e+00 (Msun), R: 1.1612e+01 (km)\n",
      "Central P: 8.2174e-05 (Msun/km^3), M: 1.4023e+00 (Msun), R: 1.1585e+01 (km)\n",
      "Central P: 9.0391e-05 (Msun/km^3), M: 1.4556e+00 (Msun), R: 1.1557e+01 (km)\n",
      "Central P: 9.9430e-05 (Msun/km^3), M: 1.5085e+00 (Msun), R: 1.1528e+01 (km)\n",
      "Central P: 1.0937e-04 (Msun/km^3), M: 1.5607e+00 (Msun), R: 1.1496e+01 (km)\n",
      "Central P: 1.2031e-04 (Msun/km^3), M: 1.6119e+00 (Msun), R: 1.1461e+01 (km)\n",
      "Central P: 1.3234e-04 (Msun/km^3), M: 1.6619e+00 (Msun), R: 1.1424e+01 (km)\n",
      "Central P: 1.4558e-04 (Msun/km^3), M: 1.7105e+00 (Msun), R: 1.1384e+01 (km)\n",
      "Central P: 1.6013e-04 (Msun/km^3), M: 1.7575e+00 (Msun), R: 1.1341e+01 (km)\n",
      "Central P: 1.7615e-04 (Msun/km^3), M: 1.8026e+00 (Msun), R: 1.1295e+01 (km)\n",
      "Central P: 1.9376e-04 (Msun/km^3), M: 1.8457e+00 (Msun), R: 1.1245e+01 (km)\n",
      "Central P: 2.1314e-04 (Msun/km^3), M: 1.8866e+00 (Msun), R: 1.1192e+01 (km)\n",
      "Central P: 2.3445e-04 (Msun/km^3), M: 1.9252e+00 (Msun), R: 1.1136e+01 (km)\n",
      "Central P: 2.5790e-04 (Msun/km^3), M: 1.9613e+00 (Msun), R: 1.1076e+01 (km)\n",
      "Central P: 2.8369e-04 (Msun/km^3), M: 1.9949e+00 (Msun), R: 1.1014e+01 (km)\n",
      "Central P: 3.1205e-04 (Msun/km^3), M: 2.0259e+00 (Msun), R: 1.0948e+01 (km)\n",
      "Central P: 3.4326e-04 (Msun/km^3), M: 2.0541e+00 (Msun), R: 1.0879e+01 (km)\n",
      "Central P: 3.7759e-04 (Msun/km^3), M: 2.0796e+00 (Msun), R: 1.0808e+01 (km)\n",
      "Central P: 4.1534e-04 (Msun/km^3), M: 2.1024e+00 (Msun), R: 1.0734e+01 (km)\n",
      "Central P: 4.5688e-04 (Msun/km^3), M: 2.1225e+00 (Msun), R: 1.0657e+01 (km)\n",
      "Central P: 5.0257e-04 (Msun/km^3), M: 2.1399e+00 (Msun), R: 1.0579e+01 (km)\n",
      "Central P: 5.5282e-04 (Msun/km^3), M: 2.1547e+00 (Msun), R: 1.0498e+01 (km)\n",
      "Central P: 6.0811e-04 (Msun/km^3), M: 2.1669e+00 (Msun), R: 1.0416e+01 (km)\n",
      "Central P: 6.6892e-04 (Msun/km^3), M: 2.1766e+00 (Msun), R: 1.0332e+01 (km)\n",
      "Central P: 7.3581e-04 (Msun/km^3), M: 2.1838e+00 (Msun), R: 1.0247e+01 (km)\n",
      "Central P: 8.0939e-04 (Msun/km^3), M: 2.1888e+00 (Msun), R: 1.0161e+01 (km)\n",
      "Central P: 8.9033e-04 (Msun/km^3), M: 2.1915e+00 (Msun), R: 1.0074e+01 (km)\n",
      "Central P: 9.7936e-04 (Msun/km^3), M: 2.1922e+00 (Msun), R: 9.9865e+00 (km)\n",
      "Central P: 1.0773e-03 (Msun/km^3), M: 2.1908e+00 (Msun), R: 9.8987e+00 (km)\n",
      "Central P: 1.1850e-03 (Msun/km^3), M: 2.1877e+00 (Msun), R: 9.8107e+00 (km)\n",
      "Central P: 1.3035e-03 (Msun/km^3), M: 2.1827e+00 (Msun), R: 9.7228e+00 (km)\n",
      "Central P: 1.4339e-03 (Msun/km^3), M: 2.1762e+00 (Msun), R: 9.6353e+00 (km)\n",
      "Central P: 1.5773e-03 (Msun/km^3), M: 2.1682e+00 (Msun), R: 9.5483e+00 (km)\n",
      "Central P: 1.7350e-03 (Msun/km^3), M: 2.1589e+00 (Msun), R: 9.4621e+00 (km)\n",
      "Central P: 1.9085e-03 (Msun/km^3), M: 2.1484e+00 (Msun), R: 9.3768e+00 (km)\n",
      "Central P: 2.0993e-03 (Msun/km^3), M: 2.1367e+00 (Msun), R: 9.2927e+00 (km)\n",
      "Central P: 2.3093e-03 (Msun/km^3), M: 2.1240e+00 (Msun), R: 9.2099e+00 (km)\n",
      "Central P: 2.5402e-03 (Msun/km^3), M: 2.1105e+00 (Msun), R: 9.1286e+00 (km)\n",
      "Central P: 2.7942e-03 (Msun/km^3), M: 2.0962e+00 (Msun), R: 9.0490e+00 (km)\n",
      "Central P: 3.0736e-03 (Msun/km^3), M: 2.0813e+00 (Msun), R: 8.9712e+00 (km)\n",
      "Central P: 3.3810e-03 (Msun/km^3), M: 2.0658e+00 (Msun), R: 8.8954e+00 (km)\n",
      "Central P: 3.7191e-03 (Msun/km^3), M: 2.0499e+00 (Msun), R: 8.8217e+00 (km)\n",
      "Central P: 4.0910e-03 (Msun/km^3), M: 2.0336e+00 (Msun), R: 8.7503e+00 (km)\n",
      "Central P: 4.5001e-03 (Msun/km^3), M: 2.0170e+00 (Msun), R: 8.6812e+00 (km)\n",
      "Central P: 4.9501e-03 (Msun/km^3), M: 2.0003e+00 (Msun), R: 8.6147e+00 (km)\n",
      "Central P: 5.4452e-03 (Msun/km^3), M: 1.9835e+00 (Msun), R: 8.5508e+00 (km)\n",
      "Central P: 5.9897e-03 (Msun/km^3), M: 1.9667e+00 (Msun), R: 8.4896e+00 (km)\n",
      "Central P: 6.5886e-03 (Msun/km^3), M: 1.9499e+00 (Msun), R: 8.4312e+00 (km)\n",
      "Central P: 7.2475e-03 (Msun/km^3), M: 1.9332e+00 (Msun), R: 8.3758e+00 (km)\n",
      "Central P: 7.9723e-03 (Msun/km^3), M: 1.9168e+00 (Msun), R: 8.3234e+00 (km)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ts=o2sclpy.tov_solve(link)\n",
    "ts.set_eos(eti)\n",
    "ts.mvsr()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8111ba1f",
   "metadata": {},
   "source": [
    "Obtain the table of results and the maximum mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3b870437",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M_max: 2.192160e+00\n"
     ]
    }
   ],
   "source": [
    "mvsr_table=ts.get_results()\n",
    "print('M_max: %7.6e' % mvsr_table.max(\"gm\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7e756aa",
   "metadata": {},
   "source": [
    "Delete rows in the table which correspond to unstable configurations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ac512034",
   "metadata": {},
   "outputs": [],
   "source": [
    "mvsr_table.set_nlines(mvsr_table.lookup(\"gm\",mvsr_table.max(\"gm\"))+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1f0cd4f",
   "metadata": {},
   "source": [
    "Plot the mass-radius curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6b891deb",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAGoCAYAAAAXR6H9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAp30lEQVR4nO3dfXRb933f8c8FCBJ8BkCKepbIC8m27NiVKcaJ2+bJgpyHeu2pA8Zbt279oyaTtWc97SLxqDtN1fQkEmV3Xbu1Mels67rT7sykndPVXRIDUpqmTZaIhJTWkR8J0g+SrAeCECVSJEUC+0MmLJpPIAniXly8X+fgSLgAwa+uyfvx/d3f/X2NdDqdFgAANuWyugAAAJZCUAEAbI2gAgDYGkEFALA1ggoAYGsEFQDA1kqsLiCf6uvr1djYaHUZAFC0hoaGdOXKlRV9TVEFVWNjo/r6+qwuAwCKVktLy4q/hqE/AICtEVQAAFsjqAAAtkZQAQBsjaACANgaQQUAsDWCCgBgawQVAMDWCCoAgK0RVAAAWyOoAAC2RlABAGyNoAIA2Jpjg+rIkSMyDGPO4/z581aXBQBYISOdTqetLiJfWlpaVt3m4+2RcfX0vS1JSkvylXvUtKFSZn2ltvsr5HIZOawUAJxpNcfhoupHtRbb/BX6zQN3ZJ6PjE0pfmVMPxpMqLf/baXSae0MVOrhezbKV1FqYaUA4CwE1Sr5K0u1r7JU+3b6M9uGrozpL374pq5PTmv/XQ1qaQxYWCEAOINjr1FZobG+Ur/2iV369wfu0OVrk/ryX5/Vi+euWl0WABQ0gmodlLhd+vS9m/Uffm6PXrowqs5vvayRsSmrywKAgsTQ3zpyuwy1tmzX1Rs39ad/+7oeNOv08TsbrC4LAAoKZ1R5UFvu0eFP79H41Iz+KPqaZlJFM9ESANaMoMqjz9y7WT933yb9/vNnNTY5bXU5AFAQCKo829VQrV9/aJeOffNlXRqdsLocALA9gsoC9VVl+u3P7NF/+c7rGr4+aXU5AGBrBJVFykvd6vjUXfpP0dd0dfym1eUAgG0RVBaqLCvRwU/dqSdfeEU3pmasLgcAbImgsliN16Nff2iX/mPkFRXRsosAkDWCygY21nj1mXs368++P2R1KQBgOwSVTdy/wy9fhUffffWy1aUAgK04NqgKsR/VL96/TT8YGGYmIADcxtFBlU6n5zy2bNlidVnL+refCOpP/3aA61UA8C7HBlWhqvF6FNqzUd84fc7qUgDAFggqG3owWKc3hsdZuQIARFDZ1uMfNfVf/2HQ6jIAwHIElU1VlZXo7s01+mF82OpSAMBSBJWN/fxPbdE3X3xH0zMpq0sBAMsQVDZmGIb++QPb9b/73rK6FACwDEFlc3dtqtFbiRusBQigaBFUBeCXHtih//WjN60uAwAsQVAVgB11FUqMTenaBO1AABQfgqpA/MsP79Bf/pCzKgDFpyTf3zAajSoej2tgYEDBYFBtbW3z3hOLxRSPx+Xz+ZRMJhUOh/Ndpu1sri3X9clpjU9Nq6I07//ZAMAyeT3ixeNxScqEUzAYlGmaCoVCc9539OhR9fT0SJI6OjrU3Nws0zTzWaotPdq8Tc/FzulffXin1aUAQN7kdegvFoupq6sr8zwUCikWi815TzQaVSAQyDwPBoOKRqN5q9HOmuor9VZiXDMpFqwFUDzyGlThcFhPP/105nk8Hp93pjQ75DcrEAiov78/XyXa3sP3bFTk7DtWlwEAeZP3yRSzIZRMJpVIJOZdf0omk6qrq5uzLZFILPhZC/WcWuph935U2di3M6D+N0asLgMA8sayWX8dHR06ceLEvO0+n0/Dw3PXt7t9KPB2C/WcWupRCP2osvGBrbX6yfmrVpcBAHlhSVAdP35cHR0dc4b4ZpmmqWQymXmeSCS0b9++/BVXAD55zyZ9+0WG/wAUh7wHVTQaVTgczlybmp0oMRtOoVBozlBff3//vFmBxc7rccvlMjQ+NW11KQCw7vI+66+1tVX79u2T3++X3+/PvNba2pqZAXj48GH19vYqGo3qwIEDTE1fwCP3bdbz/3jB6jIAYN3l9T6q5uZmjYwsPBEgEonMeV9zc3O+yipIuxqq1dP3ttVlAMC6YwmlAraroUqvX7pmdRkAsK4IqgL2yQ9s0rd/ctHqMgBgXRFUBazG69G1iWml06xUAcC5CKoCt3d7rf7xbe6pAuBcBFWB+/idDTr58iWrywCAdUNQFTivx63pVIqFagE4FkHlAC07A4q9yfp/AJyJoHKAB4N1+ofXr1hdBgCsC4LKAbwetyanU8z+A+BIBJVD3LGxSq9dum51GQCQc44NqoV6VTmhH9ViHrpzI7P/ADiSo4PKqf2oFlJb4dHVGzetLgMAcs6xQVWMNlSV6fK1SavLAICcIqgc5CO765n9B8BxCCoH2dVQpddYTR2AwxBUDmIYhiQpxSoVAByEoHKYPZtr9PI7nFUBcA6CymF+Jsh1KgDOQlA5jL+yVMkbU1aXAQA5Q1A5kMft0tR0yuoyACAnCCoHum9brf7pHM0UATgDQeVALY0B9Q0lrC4DAHKCoHKgGq9H1yenrS4DAHKCoHIol2FoeobrVAAKH0HlUPdsqdFPzo9aXQYArBlB5VAfbAzoFNepADiAY4Oq2PpRvZ+/slSjtP0A4ACODqpi6ke1EFb8A+AEjg0qSA01Xr1zdcLqMgBgTQgqB7t/u09n3hqxugwAWBOCysHu2lStly6wkjqAwkZQOViJ26VUmitVAAobQeVw3PgLoNARVA63q6FKA5fHrC4DAFaNoHK4e7fW6kVWUgdQwAgqh9sRqNAbiXGrywCAVSOoHM7lMpRmQgWAAkZQFQFDUipFWAEoTARVEWjaUKnBYSZUAChMBFURYEIFgEJGUBWBxrpKDV7hjApAYSKoikCJ28U1KgAFy7FBVez9qADAKRwdVMXej+p21V6PrtJIEUABcmxQYa47N1Xr1YuspA6g8BBUReKuTdV6+R2CCkDhIaiKxIbqMl2+Nml1GQCwYgRVkTAMQ2IpJQAFiKACANgaQVVE/JWlSoxNWV0GAKwIQVVEdjVU6fVL160uAwBWhKAqIsENVRq4TFABKCwEVRHZXOvVhasTVpcBACtCUBURZv4BKEQEVbExDKsrAIAVIaiKTFmJSxM3Z6wuAwCyRlAVmZ11FXpjeNzqMgAga44NKtp8LKyxrlJDtKUHUEAcHVS0+Zivsb5SbxBUAAqIY4MKC6sqK9H1Sa5RASgcBBUAwNYIKgCArRFURcjrcenGFMN/AAoDQVWEdgYq9WaCKeoACgNBVYR2BCoIKgAFg6AqQtsD5XqLoAJQIAiqIlRb7lHyxk2rywCArBBURchgYVoABYSgAgDYWl6DKplMqre3Vx0dHYu+p729XdFoVMlkUu3t7Xmsrvik6U0FoACU5POb9fX1SboVWItJJBJqbW1VS0uLurq68lRZ8amrLNXw2JTqq8qsLgUAlpTXoAqFQorFYku+57HHHlNPT0+eKipes1PUCSoAdpfXoMpGPB5Xb2+vpFtnV21tbRZX5Exb/eV65Z1rat7ht7oUAFiS7SZTtLW1KRwOKxwOKxKJKB6PL/rehXpOLfWgH9V7tvrKdS55w+oyAGBZtgsqn8+X+btpmpmzq4Us1HNqqQf9qN5TWVaicdb7A1AAbBVU0WiUmX4AgDlsEVSzswBN01Rra2tmeywWUzgctqiqIsD0dAAFIK9BFY1G1dXVpWg0OmdIr7W1VbFYTKZpZu61On78uNrb22WaZj5LBADYTN6np4dCoXnbI5FI5u+cQeVPtdej0YmbqvF6rC4FABZli6E/WGOrv1znRpj5B8DeCKoitrnWq/NMUQdgcwRVEdviK9eFqxNWlwEASyKoilh9VZmuXJ+0ugwAWBJBVcTcLkOpFFPUAdgbQQUAsDWCqtjR7ReAzRFUoIEiAFsjqIqcr9yjqzduWl0GACyKoCpyt+6lYoo6APtybFAt1KuKflTzbaz16uI1ggqAfTk6qOhHtbxNNV5d5KZfADbm2KBCdjZUl+nSNW76BWBfBFWR87hdmuamXwA2RlABAGyNoAIA2BpBBVrSA7A1ggryuF2anJ6xugwAWBBBBW2s8erSKDP/ANgTQQU11DBFHYB9EVR494yKm34B2BNBBW2s8eoiQQXApkpW+gVnzpxRPB5XPB6Xz+dTIBCQaZrau3fvOpSHfPBXeJQYZwV1APaUVVANDQ3p2LFjGhwclGmaMk1TPp9P6XRaAwMDeuGFFxSPxxUMBtXR0aHGxsZ1Lhu5ZNA8EYCNLRtUTzzxhBKJhDo7O1VbW7vke69evaru7m75/X796q/+as6KxPojqgDY1ZJB9cQTTygcDqupqSmrD6utrdXBgwc1ODioJ598Ul/84hdzUiQAoHgtOZni4MGDWYfU7ZqamiwPKfpRrQxrUwCwK8fO+qMf1cqUlbg0cZPVKQDYz4pn/b3fiRMnZBiG0um0DMOQaZpMpihAG6rLdPnapLYHKqwuBQDmWFNQnTx5Uvv375+z7fTp05JEWBWYjTVeXbo2QVABsJ01Df35/f552+6//37F4/G1fCws0FBdxnp/AGxp1UF15swZ3X///ZKkZ599Vrt37868lqZtRMFpoCU9AJtadVA1NTVpaGhIkhQKhfTUU09lXuMG0sLjryhVYmzK6jIAYJ5VB1Vtba0GBgYyf7/9WhVnVIXH5TKYog7AltZ0jaqlpUUnT57MPB8cHNRzzz03b4IFCgT/gwHAhtY066+2tlYPPfSQTp8+rZGREZmmqUcffTRXtSHfGLIFYEPLnlHdfsa0mNlJFUxJBwDk2rJB1dPTk9UHZfs+2JfHZWhqOmV1GQAwx7JB1dXVJbfbveTD5XKpu7s7H/ViHW2oLtPl60xRB2Avy16j6urqknRr4oRpmgu+Z3h4WMePH89tZci72Zb0W33lVpcCABnLBtXjjz8u6dbSSP39/ZKkhx56aM57amtr1draug7lIZ82VJfpXPKG1WUAwBxZz/qbnTAhvbcQ7e0L0DIlvfBtrPHq9JsjVpcBAHOs6j6q/fv366GHHpLf79fJkyf13HPPaXR0NNe1rQn9qFaurrJUibGbVpcBAHOs6YbfEydO6NixYwqHw+ro6MhVTTlBP6qVc7kMzXDTLwCbWXFQnTlzRp///Oflcrl09OhRtba2KpVK6Wtf+9p61Ic845ZfAHaT1TWq0dFRdXV1qaurSyMjI3r88cc1MDCwqjb1AACsxLJnVA8//LD8fr/6+vrU1dWl4eFhHTt2bF5IHT58eN2KRP4w8AfAbpY9o+rr69OxY8fk8/kUj8c1ODiYeW22/fzIyIi6u7t19OjRdS0W66+sxKXJ6RmVlbitLgUAJGURVG1tbTp48OCyHzTb8gOFbbbTLy3pAdjFskN/7e3tWX2Q3Wb9YXU21nh1cXTC6jIAIGPZoMp2wgQTK5xhU61XF0dZ7w+AfSwZVE888USm3fxKDA4O6sknn1xtTbDQxmqvLlxlGSUA9rHkNaqDBw/q6aefVjweV3t7+7L9pkZHR/XVr35V9fX1+uIXv5jLOpEnNeUlGp2YtroMAMjIalHawcFBPfXUUzp9+rRM05TP51MwGFQymdTw8LCSyaQGBgYUDAZ16NAhhgELmGEY3PQLwFayuuG3qalJx44dk3RrWC8ejysej6u2tlZNTU0yTXPOorUAAORK1qunz2pqalJTUxOrpQMA8mJNi9LCmQxDSqVYowKAPTg2qGjzsXp1laUaHpuyugwAkOTwoKLNx+psqi3XO1e56ReAPTg2qLB6m2u9Os+9VABsgqDCPFt85bqQJKgA2ANBhXn8FR4lxmlJD8Aelg0qlkIqPobBLb8A7GPZoHrqqaf0jW98Q6Ojo/moBwCAObJaPf0rX/mK/H6/6urq9Nhjj+nrX//6vMVqz5w5s04lAgCK2bIrU0QiEZ04cUKBQECSFI1G9cwzz6itrU1+v1+hUEgHDhxQf3+/vva1r617wciPUrdBp18AtpDVZIr9+/fLNE0NDg6qvb1dL7zwglKplCKRiFpaWvTMM8+ou7t7vWtFHm2uLdeFJPdSAbBe1rP+amtr9eijj+rUqVM6efKkJKm5uVkHDx7UCy+8kFm0Fs6w1V+uc0xRB2ADK56evn//fu3bt0/PPvvsnOtUoVAol3XBYtsDFXorMW51GQCwuvuoamtr9dnPflYjIyOZs6ts2nwkk0n19vaqo6Nj0ffEYjH19vYqGo2qt7d3NeUhBzbVePXOKEN/AKy36ht+R0dHZRiGRkZG9MlPfjKrWX99fX2SbgXWYo4ePapwOKxQKKRTp04pHo+vtkSsgdtliAXUAdjBsrP+Tp48qXg8roGBAcVisUzTREkyTTPTOLGvr0979+5d8rNCoZBisdiir0ej0czsQkkKBoOKRqNqa2vL8p8DAHCaZYMqFAqptbVVTU1NCofDMk0zE1C5Fo/H5fP5Ms8DgYAikUjOvw8AoHAsO/QXDof1uc99TsFgUJJUV1e3LiEl3RoSrKurm7MtkUgs+v6Fek4t9aAf1cpUl5VodII1/wBYa9kzqsOHD8+ZKDE4OKhnn31WiURChmGopaVFe/fu1Re+8IU13/Dr8/k0MDAwZ9vtQ4Hvd+TIER05ciTrz29paVltaUVpR92tmX/3bKm1uhQARWzZoHr/bL6mpqY5Z1SDg4N64okn9Mwzz6w5qEzTVH9/f+Z5IpHQvn371vSZWL2ddRWKXx4jqABYas1tPpqamnTw4EG1trau+jNmZwGGQqE5Q339/f3cn2WhxrpKDV4Zs7oMAEUuZ/2oOjs7l31PNBpVV1fXvHukWltbM7MBDx8+nLmP6sCBAzJNM1clYoW8Hrcmp1NWlwGgyC079Jet2trlh4dCodCCZ0i3z+xrbm5Wc3NzrsoCABQ4OvwCAGyNoMKS/BUejYxNWV0GgCJGUGFJTfWVijOhAoCFCCosyayvUvzydavLAFDECCosaZu/XG+P0JcKgHUIKizJ5TKUTrOMOgDrEFQAAFsjqLCsaq9HV2+wOC0AaxBUWNYdm6r16sVrVpcBoEgRVFjWXZuq9co7BBUAazg2qBbqVUU/qtVpqC7TxdEJq8sAUKQcHVTpdHrOY8uWLVaXVZAMw7C6BABFzLFBhdxyGYamZ1hJHUD+EVTIyq6GKg1cZiklAPlHUCEr926t1T+du2p1GQCKEEGFrOysq9Abw5xRAcg/ggpZYUIFAKsQVMia22XoJhMqAOQZQYWs3bWpRi9f4MZfAPlFUCFrH2z069RQwuoyABQZggpZq6sqU4K29ADyjKDCihiG6E8FIK8IKqzIzrpKDQ2PW10GgCJCUGFFuE4FIN8IKqzIjkCF3uSMCkAeEVRYEcMwlFaa61QA8saxQUU/qvWzq6FKr1+6bnUZAIqEo4OKflTr4yO7N+h7r12xugwARcKxQYX1U19VpuGxSavLAFAkCCqsSqnbrYmbM1aXAaAIEFRYlQ82+dU3NGJ1GQCKAEGFVWnZGdCPuJ8KQB4QVFiV0hKXZlIppVJMUwewvggqrNre7X6deTtpdRkAHI6gwqp9ZHe9/u7Vy1aXAcDhCCqsmtfj1tR0ilUqAKwrggprcu/WWr14btTqMgA4GEGFNfn4nQ068fJFq8sA4GAEFdakvNSt6Zm0pqZTVpcCwKEIKqzZ/j0NOvESZ1UA1gdBhTXbu93HNHUA64agwpoZhqFNNV6dT96wuhQADuTYoKIfVX79wt6t+qsz7F8AuefooKIfVf4EKks1OnGTSRUAcs6xQYX8e+S+zfrrH3NWBSC3CCrkzD1bavXShVFWqgCQUwQVcuojd9CmHkBuEVTIqY+yUC2AHCOokFOGYejuLTX6p7evWl0KAIcgqJBzP/9TW/SN0+esLgOAQxBUyLkSt0sfMgP6/utcqwKwdgQV1sXDd29U5KWLzAAEsGYEFdaFYRh6+O5N+vZP3rG6FAAFjqDCunkwWKdTQyOsVgFgTQgqrKtf+tAO/fkPhqwuA0ABI6iwroIbqjSTSuuN4TGrSwFQoAgqrLtf+ZlG/fd/GGJiBYBVIaiw7spK3PrMvZtpAwJgVRwbVPSjspcHmgJ6/dJ1Xb42aXUpAAqMo4OKflT28oWPB/Un33ldqRRDgACy59iggv1UlpWotWWb/uz7Q1aXAqCAEFTIq3u21KrM41LszRGrSwFQIAgq5N0vPbBDz//4gq7euGl1KQAKAEGFvDMMQ/9u/y79YeRVTc+wagWApRFUsISvolT/+sGd+qMTr3F/FYAlEVSwjLmhSg8G6/QXP3zT6lIA2BhBBUv9dLBe5R63Tr580epSANgUQQXLfXbfNr3yznX1v8FMQADzEVSwhc9/zNT3XrusM28lrS4FgM3kPahisZh6e3sVjUbV29u74Hva29sVjUaVTCbV3t6e5wphBcMw9Bv7d+vky5f04rmrVpcDwEbyHlRHjx5VOBxWKBTSqVOnFI/H570nkUiotbVVra2t6ujoyHeJsIhhGPrN0G5968V3dPb8qNXlALCJvAZVNBpVIBDIPA8Gg4pGo/Pe99hjj2lkZESRSESmaeazRFjMMAz91oE79M0XL3DNCoAkqSSf3ywej8vn82WeBwIBRSKRBd83OyyYSCTU1taWrxJhAy7XrbB6+ntxDV+f1MP3bLK6JAAWyusZVTKZVF1d3ZxtiURi3vva2toUDocVDocViUQWHB6UFm7lsdSDNh+FwzAMtX00qOuT0/pL7rMCilpeg8rn82l4eHjOttuHAm9/3yzTNBeddLFQK4+lHrT5KDyPNm/TVn+5/vOJ12gPAhSpvAaVaZpKJpOZ54lEQvv27Zvznmg0ykw/zPGxOzZo/56N+vLzZ5Ucn7K6HAB5ltegCoVCc4b6+vv7FQqFJCkTYKZpqrW1NfOeWCymcDiczzJhQ3dvqdFvHrhDf3TiNf3j20mrywGQR0Y6zyuCxmKxzKSKZDKZCaEDBw6os7NTzc3NmaG+eDwu0zRzFlQtLS3q6+vLyWfBGqlUWn/2/SFVlrn1uZbtMgzD6pIArMBqjsN5DyorEVTO8YOBYX3nlUv6tY/vUm2Fx+pyAGRpNcfhvE5PB3LlwWCd9myu1p/87ev62V31+ugdG6wuCcA6Ya0/FCxfRal++zN7NDpxU3/wwisan5q2uiQA64CgQsF75L4t+uUP79ST335V33/9itXlAMgxggqO0FDj1e88skcj4zd19Jsv6fK1SatLApAjXKOCYxiGoZ+7b7M+cke9vv69QW2p9epzLdvlcjEzEChknFHBcWq8Hv3WgTt095Yaffn5syxuCxQ4ggqOdd82n770yN16MzGmr/7fl/RWYtzqkgCsAkN/cDSXy9Av3r9Nn7pnRv/z/w1pajqlf/PTjar2cu8VUCgIKhSF8lK32j4a1KXRCXV9N666qlL9iwd2yOtxW10agGUQVCgqDTVeffGTd2roypj+MPqqGusqFd63TR43o+CAXTn2t3OhXlX0o8KsxvpKHf70Ht2/w6fj33pZ3zj9tmZoIwLYEmv9AZLOvJXUX//4vPbt9OtT92xiSjuwTlZzHHbsGRWwEnu3+/Q7j9yt+qoy/f7fnNXfv8YKF4BdEFTAbR5oCuhLj9yt65PT+srfnGVKO2ADTKYA3scwDH3qA5v08Ts36M9/MKSZlPQrP92o8lJmCAJWIKiARXg9t6a0n0ve0B9GX9Vdm6r1C3u3ys31KyCvGPoDlrHVV67f/swemRuq9JW/eUmRsxdVRHOQAMsRVECW9m736Uv/7G6Ve9z68vNn9cP4sNUlAUWBoAJW6Gd31+tLj9ytxNiUfv/5s3rx3FWrSwIcjWtUwCoYhqFP37tZB+7eqP/z4/P6qzPn9EBTnR66q4FrWECOEVTAGpS4XXq0eZvS6bR+NJjQ8W+9rC2+cj3avJWFb4EcIaiAHDAMQx8y6/Qhs05vJcb13/5+SKl0Wo82b9XOukqrywMKGkEF5Nj2QIV+I7RbY5PTeu70Ob39oze1Z1ONHr5noypK+ZUDVorfGmCdVJaV6Jc/vFOS9NKFUX39e4O6cXNGDzQG9LO761mxHcgSQQXkwZ7NNdqzuUapVFqnhhL64xOvyZD0sTsb1LzDJ8NgAgawGIIKyCOX671rWZPTM/ruK5d1/NuvyFvi1ofNgJp3+jnTAt7HsW0+jhw5ot/7vd+bs23z5s30pIItTdyc0Q8HE+p/Y0TTMyk11lXqY3du0MYar9WlATm1mjYfjg2qhdCPCoVi6MqYvvvqZV26NqGyErc+bNapeYdPJZxtocCt5jjM0B9gQ431lWqsvzWtfeLmjH4QH9Yfn3xdM6mUGqq92rfTrz2ba7i5GEWBoAJszutx6xN3NugTdzZIki6NTqj/jRG9cPaiUqm0ykvd+qltPt2/w6fKMn6l4Tz8VAMFpqHGq0/fu1mfvnezJOnG1IzOvJXU//jBkMYnZ+QypN0bq9XS6Nfm2nKLqwXWjqACClx5qVsPBuv0YLBOkpRKpfXqpWs68dIlXRydkCGpzOPW7oYq3b2lRlt95UyHR0EhqACHcbkM3bWpRndtqslsm7g5o1cvXtPfv3ZF55M3Mtu3Byp095Ya7W6oVmkJEzVgTwQVUAS8Hrfu2+bTfdt8mW2pVFpvjYzrpQujOvHSJd2cSUmSass92rO5RuaGSm2q8XL2BcsRVECRcrkM7ayr1M66Sn3qA+9tT45P6aUL13Ty5Uu6ODopvXsHS4nbpZ11FTLrq9S0oVJVTNxAnvCTBmAOX0XpnGtes6amU3ozMa745ev6/sAVjU1OS++ebdV4S2RuqJRZX6Vt/nLu90JOEVQAslJa4tKuhirtaqia99rV8ZuKX7mu2Jsj+qsz55W6bR2B0hKXtvi82uqr0DZ/uTbWeLn/CytCUAFYs9oKj+7f4df9O/zzXpu4OaMLVyd0buSG/u7Vy7o4OqmZ24KsxGVoU61X2/zl2uar0KZaLxM7MAdBBWBdeT1uNdVXqql+4QaSN2dSeufqhN4euaEfDSV0IXlDN1NzV3arLitRQ02ZNlSXqaHaq401ZaoqK2GiR5EgqABYyuN2aXugQtsDFQu+nk6ndX1yWpeuTerS6KR+cv6qvvPypK5PTs99n6Ryj1sN1WVqqLkVaA3VZfJVeAi0AkdQAbA1wzBU7fWo2utRcMP862O3uzE1o0vXJnTp2qRev3Rr0kdy/KbSkmajKq1bw43+Co/8laUKVJTe+rOyVL4Kj8pK3Ov9T8IKEVQAHKO81J2Zcr+U6ZmUkjduamRsSomxKb0xPKbTbyY1Mj6lqemUbh94nA04j9uYF2z+ilvhRg+x9eXYoFqsHxUAlLhdqq8qU31VWdZfMzWdUnJ8SonxW+H22sXrSoxPKTk2Ne+a2qxS962zwZryElWXeVRT7lG1t0Q15R7VeEtUWVoiFzMgl0U/KgBYJ1PTKV2buKnRielbf9549893/359clpLHYANSVVlJbeCzutRzWzoeW8FXbXXU3AzJOlHBQA2UlriUl1VmepWcOZ2u1QqrbGp6TlBd/napAYuX8+E3tTMe1FnSAsOW3o9blWWuVVRWqLKUrcqyt79s7Tkve1lbpV73LaceEJQAYBNuVzvTSSRVteyJZ1Oa3I6pbHJaY1PzWhsalpjkzMan5rWletTGp+a1tjUjMYnp3Xj5oyyHWPzuI15QXcrCEtUUebO/Fnhca95pRKCCgAczDAMeT1ueT1u1S3/9qxNTad0493gG383/MampnVxdOLdbTMam7y1ffYG79WeqxFUAIAVKy1xqbTEpdoKz4q+7i8Or/x7FdZVOABA0SGoAAC2RlABAGyNoAIA2BpBBQCwNYIKAGBrBBUAwNYIKgCArRFUAABbI6gAALbm2KA6cuSIDMOY8zh//rzVZQEAVsjRQZVOp+c8tmzZYnVZAIAVcmxQAQCcgaACANgaQQUAsDWCCgBgawQVAMDWCCoAgK0RVAAAWyOoAAC2RlABAGyNoAIA2BpBBQCwNYIKAGBrBBUAwNaMdDqdtrqIfKmvr1djY+Oqv/78+fNrXoGdz3DmZ9ihBj4j959hhxqc9hlnz57V+Pj4ir6mqIJqrQzD0Fp3F5/hzM+wQw18Ru4/ww418BkM/QEAbI6gAgDYGkEFALA1girPfvd3f9cxn5ELdvm32GF/2OXfYYd9Idnj38K+yP1nrAaTKVYgFxcSnYT98R72xVzsj/ewL+ZiMgUAwHEIKgCArRFUAABbI6gAALZGUAEAbI2gWgG7TFW1C/bHe9gXc7E/3sO+mGs1+4Pp6QAAW+OMCgBgayVWF2BHyWRS0WhUp06dUmdnZ2Z7LBZTPB6Xz+dTMplUOBy2sMr8WWx/RKNRxeNxDQwMKBgMqq2tzcIq82OxfXH7693d3Tp06JAF1eXfUvuju7tbpmkWze/KUseNvr6+zPNi+D2ZtdgxYqXHUoJqAbM/VMlkcs72o0ePqqenR5LU0dGh5uZmmaaZ7/LybqH9EY/HJb33SxcMBmWapkKhUN7ry6fFfjZmHT16NI/VWG+x/dHe3q6Ojg6Zppn5s7m52YIK82exfdHX15f5Pent7VUsFnP8vpCWPkas9FjK0N8CQqHQvJ0WjUYVCAQyz4PBoKLRaL5Ls8RC+yMWi6mrq2vOe2KxWL5Ly7uF9sWsWCymYDCY54qstdD+SCaTisfjme2dnZ1FcWBebF/MHpAlKZFIyOfz5bkyayx2jFjNsZSgytLsaeqsQCCg/v5+6wqyWDgc1tNPP515fvuBqVglEomi3wfSrf+p8/l86u3tVW9vr7q7u60uyTI+n0+maSoYDKq3t1eBQKBofkYWO0as5lhKUGUpmUyqrq5uzrZEImFRNfYw+8OWTCaVSCSK4jrEYqLRqOOHPbMVj8cVj8cVDocVDofV399fNKMPC5k9o3z88ceL7pix0DFiNcdSgipLPp9Pw8PDc7bdfvpazDo6OnTixAmry7BMMpnkZ+E2Pp9PLS0tmefBYFCRSMTCiqwzO7mmp6dH/f396urqUm9vr9Vl5d3tx4jVHEsJqizNzl6alUgktG/fPusKsonjx4+ro6OjaMbdF9Ld3a14PK7e3l5FIpHMOHyxKpahrWw888wzmTNt0zR14sSJogvt9x8jVnMsJaiyFAqF5pye9vf3F/1QTzQaVTgczhyYivXgfOjQocww1wc/+EE1NzcX9c/G+39XBgYGdODAAQsrso5pmvOGtYppXyx0jFjVsTSNeSKRSLqtrS1tmma6p6cns72/vz/d09OTjkQic7Y73UL7o7+/P+3z+eY8IpGIxZWuv8V+NtLpW/skHA6nm5ubi2JfpNNL/650dname3p60l1dXRZWmD+L7Yuurq7MfiiWfZFOL32MWOmxlCWUAAC2xtAfAMDWCCoAgK0RVAAAWyOoAAC2RlABAGyNoAIA2BpBBdjU8ePHC/rzgVwhqIAci0ajam1tlWEY6ujoyCwh09ramvU6b93d3ZlFfqPRqPbt26fW1tac1hkOhwkrFAQaJwI5Ntt3J5lMzut66/f7JWnJleZnO6LONpwLhUJqb2/P+RpxpmlqeHiYFi2wPc6ogHUQiUQWXNMtEAjo1KlTS35tV1eX2tvb533dejh8+PC8MAXshqAC1sFC/ami0agSicS8EHq/WCyWtzMcn8+XaRkO2BVDf0COxWIx+Xy+TPv1eDyuaDSqSCSi/v7+JUMom2G4ZDKppqYmtbS0qKOjQ4FAINNGYTYEZ2uYHYZMJBKZfkjvZ5qmYrFYUbSLR2EiqIAci0ajMk0zM3Hi1KlTSiaT6unpWfZr39+meyGJREKdnZ2Za1iS1N7ero6ODpmmKdM01dLSIr/fr0gkkrkeduDAAfX29s67PhYMBhWPxwkq2BZDf0CORSIRtbe3Z3pUdXZ2KpFIZDXDLh6PKxgMLvp6LBZTb2/vnJCSNKcp3ULPZ7ct1PKb4T/YHUEF5NhC16eSyaQGBgay+vrbu5/eLhaLqa+vT0ePHl3wPQsNGd4+CWOpCRnF3KEZ9kdQATkUjUbl8/nmhUZfX9+SZ0qzAoGAhoeHF3wtFAqpra1Nhw8f1uOPP56TeqVbwbheswqBXCCogByKRCJqaWmZtz2ZTGbCa6mbfk3TXPSMatahQ4cyQ4C3W2hYL5vXBwYGuI8KtkZQATmQTCZ1/PjxTCB1d3fPeb2zs1OnTp1Sb2/vkpMWmpub510vikaj6urqUjQazYSTz+dTa2urjh8/rlgsps7OTsXjcXV3dyuZTKqjo0OS1NHRkdkei8Uyn3M7JlLA7mhFD9hMe3u7Ojs783LdaDbUFpq2DtgFZ1SAzXR0dOjo0aN5+V7d3d2Zsy/ArggqwGZM01RdXd26TxmfvWeL61OwO4IKsKFDhw5lvdL6ai10PxZgR1yjAgDYGmdUAABbI6gAALZGUAEAbI2gAgDYGkEFALA1ggoAYGv/H4mW9noa1zzsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if plots:\n",
    "    p=o2sclpy.plot_base()\n",
    "    p.xlimits(9.5,20)\n",
    "    p.plot([mvsr_table,'r','gm'])\n",
    "    p.xtitle(r'$ R~(\\mathrm{km}) $')\n",
    "    p.ytitle(r'$ M~(\\mathrm{M}_{\\odot}) $')\n",
    "    p.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34a79ec3",
   "metadata": {},
   "source": [
    "For testing using ``pytest``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6eb9c29c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_fun():\n",
    "    assert numpy.allclose(mvsr_table.max(\"gm\"),2.192,rtol=4.0e-4)\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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}