{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "40a514e6",
   "metadata": {},
   "source": [
    "# DSH example for O$_2$sclpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "009e5479",
   "metadata": {},
   "source": [
    "See the O$_2$sclpy documentation at\n",
    "https://awsteiner.org/code/o2sclpy for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5c599678",
   "metadata": {},
   "outputs": [],
   "source": [
    "import o2sclpy\n",
    "import matplotlib.pyplot as plot\n",
    "import numpy\n",
    "import sys\n",
    "\n",
    "plots=True\n",
    "if 'pytest' in sys.modules:\n",
    "    plots=False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de741f15",
   "metadata": {},
   "source": [
    "Link the O$_2$scl library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c4fd81d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "link=o2sclpy.linker()\n",
    "link.link_o2scl()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "133c1b4b",
   "metadata": {},
   "source": [
    "Get a copy (a pointer to) the O$_2$scl unit conversion object, which\n",
    "also allows access to the constant library, then get ħc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2a7e1f0c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ħc = 1.973270e+02\n",
      "\n"
     ]
    }
   ],
   "source": [
    "cu=link.o2scl_settings.get_convert_units()\n",
    "ħc=cu.find_unique('hbarc','MeV*fm')\n",
    "print('ħc = %7.6e\\n' % (ħc))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c25ddf88",
   "metadata": {},
   "source": [
    "Use the cloud_file object to download the EOS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "33f65524",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Function cloud_file::get_file_hash() succeeded to obtain file named 'dsh.o2'.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cf=o2sclpy.cloud_file(link)\n",
    "cf.verbose=1\n",
    "cf.get_file('dsh.o2','https://isospin.roam.utk.edu/public_data'+\n",
    "            '/eos_tables/du21/fid_3_5_22.o2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e147343",
   "metadata": {},
   "source": [
    "Read the tensor which stores the average mass number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4c2e0f74",
   "metadata": {},
   "outputs": [],
   "source": [
    "hf=o2sclpy.hdf_file(link)\n",
    "tg_A=o2sclpy.tensor_grid(link)\n",
    "hf.open('dsh.o2')\n",
    "o2sclpy.hdf_input_tensor_grid(link,hf,tg_A,'A')\n",
    "hf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddb996f6",
   "metadata": {},
   "source": [
    "In order to make a plot at fixed Ye, we first need to construct a\n",
    "tensor index object. We want to include all values of nB (index 0 in\n",
    "the tensor object) and all values of T (index 2 in the tensor\n",
    "object), but for Ye, we select the value in the grid which is\n",
    "closest to Ye=0.4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "079cb4c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ix=o2sclpy.std_vector_size_t(link)\n",
    "ix.resize(3)\n",
    "ix[1]=tg_A.lookup_grid(1,0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83fc56d6",
   "metadata": {},
   "source": [
    "Create a table3d object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f8aa2de3",
   "metadata": {},
   "outputs": [],
   "source": [
    "t3d=o2sclpy.table3d(link)\n",
    "tg_A.copy_table3d_align_setxy(0,2,ix,t3d,'nB','T','A')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00950172",
   "metadata": {},
   "source": [
    "Now plot the results. Raw matplotlib works, but o2sclpy has\n",
    "a couple functions which make it easier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4c6d7495",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAGqCAYAAAC4UxuPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABBF0lEQVR4nO3dz3PbZrov+O8LkJRkxxZFxffOOfdXDN5TUzVVU9WhlD53ZjUTQ+duZjFVh4z/gTZVvW+b0eoom1EkZzlVHSn9DyRip2Z5JoSyvVMtkZ3VbM4R7LlVc2bhFkW7O3Fsi3xnAQICSJAEKZAEwO8nxZIJgQRIMXr0vO/zPhBSSgkiIqIYUOZ9AkREREExaBERUWwwaBERUWwwaBERUWwwaBERUWwwaBERUWyk5n0C85ARS1jG7XmfBhEtgMy6gj/96U83eo7//D/fxkWzHcr5vP9X/xP+8R//MZTnmofEB63d3V189tlnnm0ZLONvFX1OZ0REkTCjJaqXHzy78XNcNNv4w//570M4G+CX/8vNAui8JX54cHd3F1JKz20JK/M+LSKiwCSATkj/xV3iMy0ioviTaMv4B5wwMGgREUWclWmx4x6wAMODRESUHMy0iIhiIAnzUWFg0CIiijgJiTYvyAGAw4NERBQjzLSIiGKAhRgWBi0iooiTANoMWgAYtIiIYoGZloVzWkREFBvMtIiIIk4CrB7sYtAiIooBrtKyJH54cHd3F0IIz+0NXs/7tIiIaAILEbTY5Z2I4kxCoh3SLe44PEhEFHUSaMc/3oQi8ZkWERElBzMtIqKIsy5NQgCDFhFRDAi0IeZ9EpHAoEVEFHESQIdzWgA4p0VERDHCTIuIKAY4PGiZadBqtVowDAOnp6fY3993thuGAdM0cX5+jnw+j3K5DABoNBowTRPZbBatVgvFYnGWp0tEFAlWl3cGLWDGQevs7AyAFbxspmkCgBOo8vk8NE2DruvY29vD8fExAKBSqaBQKEDTtFmeMhERRchM57R0Xe8LOo1GA4eHh559Go0GDMNALpdztufzeRiGMbNzJSKKko4Uodzibu5zWsViEbquO/dN08TW1pYzLGjL5XKo1WpzOEMiovni8OC1SFQP2sGp1Wqh2WyiWCyi1WphfX3ds1+z2ZzD2RERzZeEQBtKKLe4i9QrqFQqODk5AWAFsouLC8/33cOFbn6d3Ifd2OWdiCieIhO0Dg4OUKlUnKxL0zRPwUaz2cTGxobvY/06uQ+7scs7EcUN57QskQhahmGgWCw6RRqGYUDXdc9wYL1e98x9EREtCntOK4xb3M00aBmGgcPDQxiGgWq1CsCqHiyVStjY2MDa2hrW1tac/Xd2dlCtVmEYBra2tljuTkS04GZaPajrel+2VCgUcHl56bt/oVBAoVCYxakREUWYQFtGYmBs7uZe8k5ERMNZlyZh0AIiMqdFREQUBDMtIqIYmFURRbVahaZpME3T0yd2UC/YWfeIZaZFRBRxUlpzWmHchmm1Wk6f12KxCMMw0Gg0AAB7e3tOB6PT01Onb+yg7dPCoEVEFAMdiFBuw2SzWdTrded+q9VCoVAY2At2Hj1iGbSIiMhhN3g4Ojpyhgb9esHW6/WB26eJc1pERBFnLS4OJ8d48eIFNjc3nfvlctm5NJSt0Wh4MqhBvWDn0SOWQYuIKPLCW6d1794959qGg9hrZLe2ttBsNpHNZnF+fu7ZJ5fLDdw+TRweJCIiAFbl4Pb2tnNf0zTU6/WBvWDH6REblsQHLb8O8OzyTkRxYi8uDuM2TDabxdbWlnPfNE2USqWBvWDn0SNWSCnlVI8QQXdFDn+rsPku0UKb0a++y8KzkcNxo/zH//4WDv6P/zaU8/nfHqpDz8fuC2vPWT158gRAdNZpcU6LiIgcg4LOoF6ws+4Ry6BFRBRx9pWLiUGLiCgWOqF1eY/3jBCDFhFRxIW5Tgtoh/Q888F8k4iIYoOZFhFRxEkItOVsurxHHYMWEVEM8CKQFr4LREQUG8y0iIgiTkqE1nsw7hi0iIgib/S1sBYFgxYRUcRJMNOy8V0gIqLYSHzQYpd3IkqCNpRQbnEX/1cwwu7uLqSUntsSVuZ9WkREgUkIdGQ4t7hLfNAiIqLkYCEGEVEMJGFoLwwMWkS0OIQCyM68z2JsEmF2eY83vgtERBQbzLSIiCJPoM3FxQAYtIhoUYj4DixxePAagxYRUQww07IwdBNR8rmzLKHEOutadMy0iIgiTkrB4cEuBi0iSraEZFVsmGvhu0BERLHBTIuIkishWZYEeD2trmT8RIdgl3ciij+BtlRCucVd/F/BCOzyTkRxZ63TYpd3YAGCFhEtqIQMDZIX57SIiGKAXd4tDFpElDwJy7Lsi0AShweJiChGmGkRUbIkLMuydZhjAGDQIiKKPCmBNocHAXB4kIiSJKFZFl2baabVarVgGAZOT0+xv7/vbG80GjBNE9lsFq1WC8Viceh2IqJFw0IMy0z/LDk7OwNgBS+3vb09FItF6LqO09NTmKY5dDsR0SKxqgeVUG5xN9NXoOs6NE3zbDMMA7lczrmfz+dhGMbA7UREvhI+NNiGCOUWd3P/KdvDf7ZcLod6vT5wOxERLa65Vw+2Wi2sr697tjWbzYHbiYj6JDzLsnsPUgQyrWw2i4uLC8+2XC43cLsfv07uw27s8k5E8cI5LdvcX4GmaZ7CjGaziY2NjYHb/fh1ch92Y5d3ogRJeJZFXnP/aeu67hn2q9fr0HV94HYiokXUgQjlFnczndMyDAPHx8cwDAPVatVZd7Wzs4NqtYpsNoutrS2nwnDQdiIiAAuTZbEjxrWZBi07g+pVKBRQKBQCbyciWjRJmI8Kw9yrB4mIKDoMw4Bpmjg/P0c+n0e5XAYAbG9vo1QqYXNzE5VKBYeHhwBm37mIQYuI4mlBhgaB2V1Py+46ZAeqfD4PTdOcGgM7aNkBC7A6Fx0fHwMAKpUKCoXCVKdyFuenTkQUY7MoxGg0Gp6ApOs6Go0GAODhw4e4vLxErVZzgtI8Ohcx0yIiIgBwer3aTNPE1taW8+9qtQrAWoJULpd9OxfVarWpniODFhFRxIXZEePFixfY3Nx07pfLZWc4EIAThFqtFprNpjNHVS6Xne+VSiXouj6XzkUMWkQUPws0n2ULq3rw3r17zhU3hqlUKjg5OXHuuzMqTdOc5Ujn5+eexw3qXBSWxfvJExHRUAcHB6hUKk6gMgwD29vbffuN07koLAxaRERRJ63qwTBuoxiGgWKx6Cm20DQNpVLJ2afRaDjzX7PuXMThQSKKlwUcGpTATFowNRoNT3ACgOPjY2iahkajgWq1CtM0sb29PbfORYkPWru7u/jss8882zJYntPZEBFNZhbrtAqFAi4vL32/N2jR8Kw7FyX+Txa/DvDs8k4UQ0JZyCyLvBKfaRERxR0vAnmNQYuIoo8ZFoNWFz8JREQUG8y0iCjamGXNrGFuHDBoERHFQBKuOhwGBi0ioqiTnNOyMe8mouji0CD1YKZFRBRxLHm/xqBFRNHELMuDQcvCTwUREcUGMy0iih5mWR4seb+W+E/G7u4uhBCe2xu8nvdpERGNRUoRyi3uFiJosWEuUYwwy6IhODxIRBQDXFxsYdAiIoo4ycXFDgYtIooODg0OlIT5qDDwE0JERLHBTIuIooFZ1hAsebcxaBERxQCHBy0MWkQ0X8ywaAwMWkREEceGudcYtIiIok5aZe/EoEVE88ShwcDCWlwc93yNnxgiIooNZlpENB/MsgKTCK96kJlWxLHLOxHFn7VOK4xb3C1E0GKXd6KICSnLEkr8fwnTeDg8SEQUA6wetDBoEVEs2VmWUARkJ/m/0dkRw5L44UEiIkoOZlpERBEnJTMtG4MWEc1OCAUYi1p8kYTKvzAwaBERxQALMSyc0yKi2WCWRSFgpkVEFAOc07JEJmg1Gg2cnZ0598vlsrPdNE1ks1m0Wi0Ui8V5nSIRTWpaWVb3eYUKyHb7xseIKgnBoNUVmeHBs7MzlMtllMtl5HI5NBoNAMDe3h6KxSJ0Xcfp6SlM05zzmRIR0bxEImi1Wi0cHx8795vNJrLZLAzDQC6Xc7bn83kYhjGPUySiSQjlxlmWUMTQLGtRyJBucReJn3o2m4Wmacjn86hWq8jlctA0zRkWtOVyOdTr9fmdKBHRPHTXaYVxi7tIBC0A2N/fR6FQwKNHj9BsNgFYGdj6+rpnP/t7bn6d3Ifd2OWdiGKHqRaAiAStVquFo6MjHB8fo16v4/DwENVqFdlsFhcXF5593cOFNr9O7sNu7PJONAOzKnFXhHWjhRCJoPXNN99A13UAgKZpODk5Qa1Wg6ZpaLVazn7NZhMbGxtzOksiovnh8KAlEiXvmqb1DfttbW1B13UcHh462+r1OiqVyqxPj4jGMe1rZbmf3+70LkQSRr6GmlVHDMMwYJomzs/Pkc/nRy4/mvWypEgELV3XcXR0hGq16gQv+43a2dlxhgq3tragado8T5WIKLHsJUX27998Pg9N06DrOvb29pwq70qlgkKhAE3TBm6flkgELeD6TepVKBRQKBRmfDZENJFpz2MNyLKgRGKmY2okZtMRo9Fo4Ouvv3ama3Rdd9bM+i0/0jTNd/ug3+dhiEzQIiKiASSAGQQtu5GDzTRNbG1t+S4/qtVqADBw+7Qk+88TIpqNEBYRW08zZB6rJ8sSQkCoKmDfKBR2EGq1Wmg2mygWiwOXHwVdlhQmBi0iohiwLgR589uLFy+wubnp3I6OjnyPV6lUcHJyAgADlx8FXZYUJg4PEhHFQUjVg/fu3fM0J/dzcHCASqXiZF2apnm6EdnLjwZtnyZmWkR0MyENCwYuce8ODUJRAFWFUFWITBrq+hh/4cvOzU545sJZoxWkmMMwDBSLRacC0DAM6LruGfar1+vQdX3g9mlipkVERACs6sFSqeTZZpezD1p+NOtlSWMHrR9++AGmaTrVJHZz21/84hdTOD0iiqw5LCIGrkvcRSplbU+lrH+nAv46i12W1TWDxcWFQgGXl5cDv+e3/GjWy5IC/ZSfP3+Ozz//HM+ePYOmadA0DdlsFlJKnJ+f47vvvoNpmsjn86hUKvjggw+mfNpERAtE8srFtpFB6+nTp2g2m9jf38fq6urQfV++fImjoyOsra3hV7/6VWgneRO7u7v47LPPPNsyWJ7T2RAlxLTL2916FxHbX9MpQLHms5BSITNpKP/d36Dzf//Tjc+NomvoJ+/p06coFovY29sbGbAAYHV1FY8fP8aDBw/wxRdfhHaSN+HXAZ5d3okodnhpEgAjMq3Hjx9P9KT379/Hb37zm4keS0REfjg8CLB6kIjGNauhwQHDgkIIq9Q9nbaKLzJpyKU0ZCaF9q1Mcn+1JyBLCsPIT9/z589ncBpERESjjQxavH4VEQGYbn9B93O7FhC7ewsKVYVQuwUYmTSwlIFcWYJcWULn9hKu7izh3d00fvpfN298jpHEOS0AAYJWrVbDzs4Ovv3221mcDxER9bK7vIdxi7mRc1r7+/t49OgRXr58ia+++govX75EsVjkWiyiRTHNRcSjro9lZ1uKNY8FRYVIqUAqBXlrGTKTQmcljfZKCle3VLy7reLtewK37Od0LyQO6XXQfI0MWo8ePQJglbPb/z45OcHx8XGk1mMRESWZTMDQXhgmqh588OABHjx4gJcvX+J3v/sd6vU6tre32cqJKEmilmGpilUtmFIh0yl0bi+hvZxC+1YK726reHdL4N17Am/vCPw///C3+A+7/8V7TNmJd7bFoAXgBiXvP/zwAw4PD/H111+j1Wohm80yaBER0VSNDFqvXr3C3bt3nX9/8803ODw8RL1eR6FQwFdffYW///u/n/qJEhEttAQUUYRhZNCqVCoolUo4PDxEtVrF6uoqyuUyvvnmG9y/f38W50hEszSr7u2jhgWFAqjda2a5hgWRTkEupXF1J3NdfHHbGhp89x7w7q7E1Z2YdnIfQnB4EECAoHV4eIijoyM8ePAA3333HR48eDCL8yIiIltC1liFYeSfVPaVKeMasHZ3d61Fiq7bG7ye92kRRU9Ii4etpxq9gNhZPNxt0SRSKYh0yiq86LZnEuk0RCbTbdWUgVzuLiR+L4M32TR+zqr4eU3gTU7gTQ54s97Bu7U2lNybUF4HRc/IT+j29nagDu9RxS7vRBR/IS0sTsC82Mig5Vdk8erVK/zud79z+hJ+//33oZ8YEc1IyBlWX5bl16IpaIa1vAS53G3XdGsJnfcyuLqTwbu7afycVfBzTsHPOYGf1yXerLfRzl1hKfcaudUf/U8wrlctBtjGqWvsT+rJyQk+/vhj1Ot1mKYJwLoUCQMXERFN29jrtGq1Gs7OzgBYAQywglaj0Qj3zIhoekJeZDv2AmJV7f67p0qw+z2kVMiUalUKplXIpTTaSyrayyqubqt4d1vBm5zAuzvAu/c6aL/XgXrnHW7d/hnZlZ/x/spfsP5/3cJ//U8/hfo65yoBWVIYxg5av/zlL323CxH/sVIioshi0AIwwfDgH/7wB/z5z38GcB2onj9/jj/84Q/hnhkREVGPsTOtnZ0dfPjhh1hbWwMAp4WTPVRIRBEX4tDg2MOC9tWHVWX4sKBqLSSWSyo6mRQ6SyquVlRc3VLw7raCt3es4ourOx3I966w/N5b3Ln1M3IrPyG39BPeX/oL1tM/4r9iKbTXOlf2pUlo/KC1urqKf/7nf8bvf/97mKYJTdPYxomIaMqS3hHjiy++wG9+85uR+40MWt9//z0+/vjjvu29gWrQfkQUAXMovHDmud0Zll18kUoNz7DSKmRaQWcphfaSgvaKine3FKuT+x2Bt6vAu/UrqHfe4fatN1i79Rq55R+xvvQTcpkfkUv9iPfUn4GkZFpArOe0Xr16hWazCdM0YZomWq0Wzs/PnW2tVgumaYYTtI6PjwMFo6D7ERHRYnj27Bk2Njbw8uVLAICUEtlsFpqmIZfLodls4uHDh1hdXUW1Wg30nIF7Dw4jpYQQAr/97W8DHZSIZmQK14/yXTxsG5Vhqaq1oDid7suwrIa41q2TUdHJKLhaVtFeVnC1Iq6vl7UKvM12sJz7Gdn3fsLa8mu8v/wX5DI/YS31E1ZTP+GO8jNuK28AvB/666fg7t+/D9M08emnn2Jra6tvhO7p06d4/PgxACCfzwd6zkBBCwA2NzehaZrvPhcXFzg4OAh0QCIiGl9c57Sy2Sy+/PJL/PGPf8QXX3yBcrnsXO7KvVQqaG/bkUHr0aNHAIA//vGPqNfrANA3DLi6uopSqRTsFRAR0cL58MMP8eGHH+Krr76CEAK/+tWvJnqewNWDH374ofPvk5MTCCGgaRo++OADAMGj5Kzt7u7is88+82zLYHlOZ0M0A1O6pPzAYUF3WTvQPyxoXx9Lsb4KRQGWMt5hwZTiDAvKtGJ1v8hYw4JXKwre3QKublkdMN6udSBzb/HfZF/h3spf8P7Sj8imfsJq6jVW1Z9wS3mD28rb7vBggiSk5P3Ro0d4+fIlnj59imazOfbjJ/p0P3jwAB9//DHW1tbw/fff49tvv8WrV68meaqpY5d3Ioq9sJrlRmSIcXV1FY8fP8Ynn3yCL774Aq9evcK3334b6LFjr9NyOzk5wZdffgnDMLC9vc1CDKJ5mVJ2ZT31gPL2AVce9iwe7mZYoruI2Oo7qEBm0lZfQTvDSlvZViejoJ1R0FkSuFpW8G5F4OqWwNUt4OoW8Ha1A7n2Ftnsj/jgThO59I9YS/+EVfU17qivcUt5g2XxDsviHTKiPbX3ZC4iEnDCZA8ZPn36FJ9++ina7dE/s7GD1g8//IAvv/wSR0dH2NjYQLlcxnfffTfRCRMRET1+/Binp6eB9g0UtF69eoXDw0McHh7i8vISjx49wvn5Oe7fv3+jEyWiG5hidmU9/ZAFxHZpu19Zu92eSSj9GZa9gHg5bZW3pxR00io6aSvL6qQF2hmB9pLwZFntFeDqlkTn7hXurr7Gv119iX+z3HIWEd9W3uB2N8tKiyukRRsZJCvTimv1YFD7+/uB9hv5qf+7v/s7rK2t4ezsDIeHh7i4uMDnn3/eF7B2dnYmO1MiIhotQXNafoImQSMzrbOzM3z++efIZrMwTRPPnj1zvmcvKr68vMTR0RH29vYmP2MiCmbWGZbf4uEhC4eda2XZGVZKhbQzrm6lYGcljU5KsSoFMwpkSlhzWWmBdga4WrayrHe3gfaKRHsZuHqvjcydt/iru6/w729d4q8yLdxRXuO28hbLyttultWGik73a4R/Q9PERgatcrnsrFge5vz8PJQTIiIiH4zBAAIEre3t7UBPVKlUbnwyRETUT8jkz2kFNTJoBR1nZFEG0ZTNo6wd8F887BkCFNfDgnbXdlUBFKVvWNAucW+vpKzii5RAJyPQTgvr32mB9hKsXoO3reKL9jLQWe4AK22svvca/+52C/9uuYl19S9OiXtGtJEWV1AhnaFBRUg8ffZf8Pj+/zC1941mb+j/BU+fPsXz58/HftJnz57hiy++GPtxR0dHMAzD0+230WigWq32bSciWihShHOLuaFB6/Hjx6jVatjZ2QkUvF69eoVPP/0Uv//97wNdF8Vte3sbuq5D13Wcnp6i0WgAAPb29lAsFp3tpmmO9bxEsSUU720qhxD+hRfOwmCr8EKo3ZZL3ZtIpYB0yrouVioFkU5bndvtRcOZtLWAeCkDuZxGZzllFV8sp9FeTuFqWcXVimLdlhW0l4RVfLFsF2FcF2B0VjqQSx2kl6/w/q0f8dfLLbyf+jOy6k9WmbvyzinCcGdZKmSyijESXj0YVKCGuc+ePXO69Gqahmw2i3w+j1arhYuLC+eCXvl8Hk+ePBl7qNC+AJjdRd6u1zcMA7lcztkvn8/DMAyUy+Wxnp+IKO44p2UJtLj4/v37+PzzzwFYQ3/21SdXV1dx//59aJrmaag7LsMwkM1mneG/ZrOJcrkM0zSRzWad/XK5HGq12sTHIYq8KZezXx/mBtfEEoqVXfUsHHbmsFLWNmcOS7VbNSmQqlXa3l6x5rPaaVhzWSlApmAtLF6Gc+usdCDTHYhMB0vL7/CvV/6Mf5V+hfXUXzwLie2sSun+ZrczLCUJqQV5jN3G6f79+7h//36oXd3tIFgsFgFYQ4WGYaDVamF9fd2zr19XYL9O7sOwyzsRxQ7jL4AJu7yHLZvNYnNz07mfz+dRq9WQzWZxcXHh2dc9XGjz6+Q+7MYu7xQpU5638h6qZw6rd/5KVa4XD9vzV2q34a09h5VJQ2TSQLo7f7VkzV/JpTTkUgZyqTuHtZRCZzmNzkoK7WUV7WUVVyvdKxEvCVwtWS2bOhmgk4bz7/YScLUs0V62sixkOlDSHaxk3uH9zF+sqkHxxrWYWPbNYynd20AyZhFAXpe93/Q2TKvVQrVa7VvC5E4k3Mug5lEoF4mgNeiKyJqmodVqOfebzSY2NjZmdFZERIvl7OwMADy/dwHrd2+pVEKpVPIEtHkUykUiaOm67hn2Oz8/x9bWVt/2er0OXdfncYpERPM1g+pBXdd9k4iHDx/i8vIStVrN+f6gQrlpu9H1tMK0s7ODg4MDaJqGjY0NJzjt7OygWq0im81ia2trYFZGFCszKri4PlyAqw7bRReA7zWxnKsOq6o1LGgXYqgKpKIAqZ7ii5RAR7X6CcruwmGr6MLq4t5RAalahRhStQsxgPayRGdJQmYkkJYQqoSidHA78xbvp/+CO+prLCvvPOXtQH/xhZq0crs5vhzTNCNTKBeZoFUoFFAoFAJvJyJaJPOMweVy2QlQpVIJuq4HLpQL28ig9Td/8zfQdR0bGxvY3NzEL37xi6mfFFEizTi7sg45OMMaWdbee00sIa7L2dOpbvd2AamqgCr6MizZzaqk6s6yuiXuaXQzLQGpWFlWJ929ZYDOUgdIdyBSHSiqhKJ2cDv9FmupH3FbWB3d3cUXtsRmWSF68eKFp/CtXC6PXPvqzqg0TXNGv3obpfsVyoVtZNB68OABfvvb3079RIiIaPru3bvnFFwEYRgGjo+PcXh46NmuaRrq9bpzf1aFciP/9Mvn81M/CaJEmkEbJv/DiqGl7UJVrLZMvWXt7tZMmTREKuUpbUfGVdq+nEYn0y13X1Kt8vYlFZ2MgvaSivaSdZ2s9pKCq+Xu1Ygz3cXDaXQXE3fnslQr6+qoQCcj0V6SQEpCpDoQioQQEooi8V76De4qr53GuO4sy13i7s6y2MbJ5zYmTdNQKpWc+41Gw6kYnEeh3MhMyxlCICKiRLOzKnvdVbFYhKZpznos0zSxvb3tFMTNo1BuZND68ssv0Ww2oes6Njc3cffuXd/9nj9/jg8++CDs8yOKlznMW10fevjlRZw/QLtZlmf+yq9C0N2ayb4CsV0lmFbRSVmPkykBqQh0nHktWHNYqnAqAzupbrVgqptRpay5LCi4/tqd15KZ7nyWKj2Z1m31LW4pbzwLigFvq6YkZ1mzmKazm5b3srsV9ZpHodzI/8M0TcPZ2RmKxSLW1tbw0UcfYWdnB99//z1evXrl7MfLhhARTRG7vAMIELS2trbw3Xffodls4p/+6Z9QLpdxfn7uCWK//vWv2ciWiIimbuTwoHuiTdM0aJqGR48eAbAWnJ2cnKBWq81kJTRR5MxxONA6fMArDivdbfa1sRRxXcbuuQqx4j8sqIrrzu2KQCfTLXO3hwVdN2to0Cpn76TQHSaEU3QhFe/QoFTsIUOrEAMpeT00CADCWlx8J/2zp9QdGD00GImWP2FJQJYUhpE/U/tijH7sAPbNN9/g8ePHoZ5YWHZ3d62L2Llub/B63qdFRBSYwGwa5sbByKB1//59/PrXv/bMX/mJamm8Xwd4dnmnic2pjN17CkNK2t1l7T1XHBapFES6W9aeznivOLyUsUraezu2u0rbO8spdDKqc2unrbL2TkZBJ62gk+mWtqetbKttLyhOW0UW0l5MrOI6yxLo/kYGOmkJmZaAXeoOQAhpJYNC4j3Vukpx0AIMBYCapOJnzmkBCFg9CAAnJydYX18f2BHDHjIkIiKalsC9B8O86CNRrMx53so2yfyVU9puz2MJYS0iHjR/pVrzVtY2YbVlUrul7Up33kqxMimpCEDBdasmex7L89VePCw8c1juLEsKOA1zpXs+yzWeJYTEHfVnZNDua93k1ptlReMnF4KEDO2FITINc4mIaAgGLQAJ+kOEiIiSj5kWUa+IDAcCYwwJAsG6tQsBZNL9Q4J2p3ZFue7Y3h0StEvbpeh2ZFevu2DYQ3/u4cCO6u0paO8nFVjDgT0FGFJxl7tfF2G4qYrEe+rPUESn7+3o7ejuHhrMJKkNHTMtAAxaRESxwDktC4MWLbYIZVXAgMzq+pvW1yDXwurtJai4sqp0qn/R8ICiCysT6hZgdP8NYfcOdGVaSvAsy86w3KQqrSKMlLcIQwirCENVOritvPH0HAT6S917s6xo/XRviEELAOe0iIgoRphp0WKJU2Z1vZP11c6wXNkVgO5cleKUtI9qyyTTKWfeysmyBsxfQfFmWM5Xew7LJ9salWVJ+yXb5e6KfX0t6TufBQCq0nFaOAHehcUAeq5eDHTXLyMdsZ/3xBKyMDgMDFpERDHAOS0LgxYlW4T+0g6UVV3vbH3tXTQ8YMGwJ7sSwn/RsF0hmFa7WdCA+atuJgUheqr8utmUcr2A2Gp2682y7H2guB8Lb5blvE50r6Xlqhz0aZKXVjpYVt55tg2qGlQBpIVAWihQffM2ijMGLSKiOGCmBWABCjHY5Z2IkoBd3i0LEbTY5X1B9HZgj0AXdt+O7IMf2F0MrF53au92Z3c6tadTEJm0tUDY7tK+tOTt1L6c6XZozwDdTu1yKQ25nEJnqXvLqJDuDu0p4Xy9XliMbr9A4SqqcC0w7h0adK6N5e0z6AwRugow3KN2UgU6aQDdnoPO2+Haxy7EUHx+6/YVYQgBFdZNGfQrLkmLjhcMhweJiOIgAVlSGBi0KJ6ieMXgyZ7I+upXcAH0F10oqlXiPqTgwi5j7/3qZEzdTMouXbeKLbxFFwAghV2cgb5CDE+3dvumegsuZG+GBXgzrO73OqrrasWA7xjWknqFtLga+DbaRRgZIZCGVYShQIE6KKOSMYsALHl3MGgREUWcTxORhcWgRfEwx/mpkJ/w+t9B2zG5S9pTqYmyKwjhug+n1B3oaV5rn4+dYcG9sLg3y+rJuIR3G3qumTXot651VWMJxS53R/+UU0Zp+1+xuPtvVQBp65BQhBieZVGsMWgREcUBhwcBMGhR1Mw4owo9kxp8IOvrTZvdpnpbM43OrtyLhp3FxH7zV3Ddd+amXFmWa66qN6Nyt3jyjGUNWlCM6+NIFdbC4gE/CiGA5dQ738uSoHv43qrB6+8NqR6M2bxWEsrVw5D4knciIkoOZlpERHHATAsAgxbNS1KHAb0Htb66y9lVtbvNNSRol6+7ry48oKTd6dDevcKwt8BiyJDgsLJ2wDss2L3vDA36DAm6Fw7bndoDlbo7x72+30lJQB3esiGjXEGF//Cg3W/Q7upuF2EMHBoEYjc0CIBBq4tBi4go6hLSgikMDFo0fTPKquaSTXlPwPrqzqwAb8HFsGtf+S0W7m6zMytPh3a7+EG1r391nRn5ZldAX1m7Tbp/RL1ZluItsnA6wTv3r/f3tGzyK8JAz3GU/uto+RVkrKj9Hd7tKxXbRRgjFxR7jh2/QgyyMGgREcUBYyyABQhau7u7+OyzzzzbMlie09kk3AwyqrlnU26TLBS22zANy66U/jJ2p5xdVZxrVdlXGA6SXQHBMixnP1e2BHfZu6sBrrvFU+98lvd96vnqIlUr03IvLPY8tLttRX3naYzrpgJIBy11dw4cvwjA4UFL4kve2eWdiCg5Ep9p0RQsWkZlG5VZdf/tVxXoLBS22zAJ4d+GSVxXBloZVU9lYGrC7Mq1DejJsHz2lz3Z0/Xtej7resHxBPNZ3WNIVQLK6I/ULeUtAPRlW6qw5rPcrZsGZllSXv+c4jinFbPTnRYGLSKiGODwoCXxw4NERJQczLRosCkOA0Zy+M+P6z0Qql3SPmA4UHSH+tzDgc5NXC8O7i1jV5RuA73r3oGwS9i75eR+178CMNaQIOAzLOizj2cBcM/Qn19vwUmLMACrw7t7YbFftboQErfUt33bFddNhRhe6m4PDbq/xgmvp+Vg0CIiigMGLQAMWjTloorYZFS2YQuEuy2YPFcS7m2/5OrI7im2sLOrtNqfWXUzNSngLBTuXeBrf9++BtY42RUQLMPyK3UfVIThuf6WT/Y1rAjD/T2pAFDloETMsaS869vmXKm4u6gYGFDq7i7AsO/HjADntGyc0yIiotiIXKbVarVwdHSEJ0+eAAAajQZM00Q2m0Wr1UKxWJzzGcYQs6nB/DIrpbutd4Gwsxi423rJb4Fwd+7KybR65q1kWvVmVu4GtwNK2QFvg1nf7AoInmH57NubqfVmVL035/sI3iTX7xw63RZOwwghsSyu0O55Ivd81kgxzK76JOAlhCFymdbe3h4uLi4894vFInRdx+npKUzTnOPZERHNh5AylFvcRSpoNRoN5PN5575hGMjlcs79fD4PwzDmcWpERBQBkRoebDab0DQN5+fnAOAMC9pyuRxqtdqczi7iWJ4+2qgu7N1/DyxhT6U8Q4GeQgt3v0DXsKBzvSt7SFAZMRxonZh3OND19UZDggMe4+5J6D4H63veIgx7W1/xheu5AxdhCKv3IEb0HQSAtLhy/m13eFdhDQ0O64ShBBk6jAOWvDsiE7QMw4Cu655MqtVqYX193bNfs9mc9akREc1dWNWDw56m1WrBMAycnp5if3/f2T6otmAeNQeRGB5stVqeYUBbNpv1zG8B8N1vd3cXQojAtzd4PbXXMjVCGX678dOLgbfYcr8/qtVdXaRTEKpqZU3pNEQmA5Hpfl1egljKACvLEMvLELdWrK/LSxAry5ArS5C3liGXlyBXltC5tQy5koFcyaBzK4POShqd5bT1dSWFzlIK7SUVHfuWUaxbSqCTViBTAp2UgEwpTn9Ba3ExrAzG00FdeLu0+2Q5Thf2oe/J6O/7FX04x/AUX3jL3cctwoBd7t79bTxobbAQQFq0AcBTjGH3HBwkMVmWTYZ0G+Ls7AyA9TvZbVBtwTxqDiIRtI6OjmCaJqrVKmq1GhqNBgzDgKZpnjev2WxiY2Oj7/F+ndyH3djlnYion67r0DTNs21QbcG8ag4iMTxol7e76boOADg8PHS21et1VCqVmZ3XzHFeanI9HdiBnvJ1u0zdfW2rYeXr9nyV/RyqAplylbHbi4MFnKsI97ZeGquM3f0VweaugACZ1YjnGDSfJQX6My6fDG9g+ya/c/XLDNVu66YRY19pcYV2z9/Y7vZNvfNZicuyMJvhQT/DagvmUXMQiaBlazQa+Prrr2GapjPHtbOzg2q1imw2i62trb6/AoiIFsKcCjEG1RbMq+YgUkGrUCjg+Pi4b1uhUJjTGYWIC3zDM0kVYCp1fQ2sUVmVXyWg0s2o3Ne2srOmIZkVgP7sym/uB8GzK+DmGdbA5+2rVvRWDjrfm6BysO9Y3Tm7YX1uAauKUO1JM+wsa+QViqnPixcvsLm56dwvl8sol8sD989ms05Fty2Xyw3cPm2RClpERORDhjc8eO/ePafgIghN01Cv1537dm3BoO3Txj9TiIjiYAbVg350XfcM+9Xrdei6PnD7tDHTCguLKKbHLpyAz1Bg7zCgq/O655pWKdXZ3+kP6Bke9A4FXhdVuLquq4rvUCAwYjgQGFgCPs6QIDDmsOCI5+srwnC/FvuxfkUX6Pm+3/aeY/sNG0pVAsroFMIuie/0vPhB5e5JLMKYFcMwcHx8DMMwUK1WnXVXg2oL5lFzwKBFRBRxAQosQ2FnUL0G1RbMo+aAQSsIFlHMxqgCC7+sSlH9Cyxc2ZRUrzMtbyY1OKvyK7Rwrh5s/7h6Wi8BAzIO17/7Miu/fX2EmWENf9yA19GTVQ18HWOQCqwWTr6ncf0b2mrb1LH+jW4LJzG43D2xEtDsNgwMWkREMcCLQFoYtDgXNXuTNq51LwS256rcGZXfXJUQnutZ2QuCR2VVED5zPhNkVjbfzCTAx2PsDCvg8/oey6fcfdgxehcWD3ovBj2P7F4Qy69ZrudQPt8f9LZwPiv5GLSIiKKOXd4dixu0bphhMYsKaETln7N91DyV+3Ig7ozKXgjczZJ656qclkTuVks+WRUwohoQrsc6r633tXrvzjS7GuP5gcGVg3779WVKvpNQ4x4f3Z/R8Ga5AKB2r2zcf+XiBZnL6hKdeZ9BNCT+p+7XAT6WXd6JiGgxgha7vBNR7M1pcXHULO7w4AAc9pvApB3WexcAO8OB1/92egJ2hwfdw4BOocWgAgvrRKzt9tWDAf+Fwa5z7y228G5zv+7+t2JoKXhEhgaHHn/AMJ9vN3dXCfxEx1LkyL6DAKAq3nExq++g/W/FtT3Z/++yetDCoEVEFHUSXKfVtZhBSzCjGtuobGpYSyU707IX/vYu/rX3U9xZmBhZXOFbtu6cS2+xwYSZVffxvv929r1ZdgXMP8Oynqf3/Rjy3MOO6fresCxMqrDShxEphJ1ptV1vkiqshcW0eBYzaBERxQyHBy0MWuTVk1EFzqbs+Sm/jEpVr7Op3jkrvzkqwL9kHegvW+/u25dV2fvCW77u3Ifr8fDJCAKUec89wxrjOKOMWkjst2/fwuKg52I/RgGEMvo3cbpnTktd1ASLQQvAAlQPEhFRcjDTIiKKOAEOD9oYtBbROL3/VNUzBCjcl6XvGerz7VbRW0RhHdA7LOhsw3VxBeAMBQLoPq+1eeRQoG1IoUXfvq7nGHjfedzNhwWBEIYGb2pYkYTfeyF6hgZHPM/IUviAv4XdJe9Kd4ys961Lerk7pGT1YNe8/7chIiIKjJlWUg3q+WdtCF5Q4ddJ3ZVZDSpRt4/TW2gBYHhG1VNgYW0TgbOqvu3jFFuM2B65DGvC5EL2vicBsqbQCNdXxf94vV3d00rbc19xvi7W39wcHrQwaBERxQGDFgAGrfgKMC8lnPmoMVsnOYt2hX95unVAV3Y0aDvGz6i6+3myGjF+VtX3Pfv4fobO7YxIPcbMTOY+jwVASBnKlYf9nzzYbnaGF+Q0Mup1pqV20w016XNYPphpWSLwv9B0+XZ5l+zyTkQUR4nPtHZ3d7G7u+vZdlfJzedkxjEqk7LnpOx9/K5F5Z6T6qnwkz0Vft2D+GZNfQt+AW/G1D2H3jmnUYuArcd7MyrPc7pef3/W1P+Xtm+12gTZVd95TfB43+dM4J+IA69KPOpnIayGuUHSh4xyNfT7ia8cBKyhwQ5TLWABghYRUSIwZgFYgOFBIiJKDmZa8xJ0+G9Qx/QRC3yl/f3eknRgwDCftyjCb9GvZ7/eIUD3YwIsAra+73o/BhVXOO/XiCILZ78xtzvPF3CIacKRqNgMDUpYBRJyzOtkjfu+iOvHjHrrM92S9073g6UCUKZVSBJhLMSwMGgREcUBO2IAYNCaDr8sym9hr72PXxGFEEAq5V9E4X6+YZnUoCzKPg/X1/5rTIn+EnXn/OF9Xk/2NCSj6rmGWZDiCt/9ep7fV1iZVdDjDTpOXDKsaet5/2T/x2SgJfW6EENd4IkdZloW/i9FRESxwUxrXEGbzarKRIt6nYzJ/T34LOB1ffVkDQP3sb/vMzfll011n2tgifqguSnXMXy/N0425TxmyPeCfB+zy6yc48X4z8Gx57OmbEV9N+9TmD8JVg92MWgREUWcgNXJhBi0vFxZ1NDLdfi0Ruq7ZEdK7cugrhfvuuapMEEW5bOfta+9X/exSu9z9TxW6Xn+AdnUoCq/vu/1nluPG2VTQfdxjjXbzAqId3blS/rc/HZzPlc994exfgsHOg2/xcV2G6eFWFhMHgxaRERx0Bm9yyJg0CIiigEOD1oWMmgJCIhMpntH9A/7uYoonPJztTvc11uC3h0WvC4rt4cE7aFGXB/H/RU3HPq7aSGFZ1vPGzTmQt+h+3oeN+L74+6HCYYBJzjGyHNI2rAg4CwyDsWQ5+m9dpafpRG9B2mxJPF/Nw+/Lu8/s8s7EcWJ3/zipLeYS3ym5dflfTX9r6DcvXNdONHb7qjvGlI9hRTAZBlUz34jM6kRWZTnWEGKKFzP33fcQec4bL++xwXYZ5z9PMe/4Z/9Ic/XJzO7khAQQ3/WQs7+915aXF9PK4lvezCSHTG6Eh+0iIiSgB0xLIsZtFIq5NpdT5Y0rJksMOAv/QBZ1NC5qN79/LKoCeeifL8/6VyU+3yDmjCzCeWKulOqgk5kduVnyHzW2IuOB+0vgz2ZEMCywsXFdG0xgxYRUdxweBDAogYtISCX0gAmz6Csx9r7j6joCzoP5dnmfe6+Yw46v1H7+plRFgWElEm5zWBtadIyrNAqLm/w3g8a6vKrJnTPaQGAGvZnKA4kILhOC8Aiz2sSEVHsLGamRUQUNxweBLCgQUsKQKqif1jthh3Re/cN+2q9Q/fte2yAfcJ4jI/QhwCBmQwDukV+SPCm74ewh65ndLxeYwx1LfkUYihQFq/vIGMWgAgFLcMwYJomzs/Pkc/nUS6XAQCNRgOmaSKbzaLVaqFYLM75TImIZo9tnCyRCFqmaQKAE6jy+Tw0TYOu69jb28Px8TEAoFKpoFAoQNO0mx1QCMi06lt+PlbhhF/peUKzJ2BKGVSvOf7xHPnsaoYmXUTs+RyPWKRsrZcVECOOtMhXK6Z+kfjftNFo4PDw0Lmv6zoajQYMw0Aul3O25/N5GIYxj1MkIpovKcO5jbC9vQ3DMNBqtbC9ve1sbzQaqFarMAwD1Wp1mq90qEhkWsViEbquO/dN08TW1pYzLGjL5XKo1Wo3P6AA2ssp599A+NnT1DKnmzxugJlkULYITUPENrOa0XvoZFvTON4YffDS4rphrhqhz89MSczs0iTNZhOlUgmbm5ueZGIqo14TiMz/tnZwarVaaDabKBaLaLVaWF9f9+zXbDbncHZERIvh4cOHuLy8RK1Wc4JSlEa9IhO0bJVKBScnJwCsQHZxceH5vvuNs/l1ch92e/P2zzN5LUREYRCQEDKc2yimaaJaraJareLo6MjZ1jvqVa/Xp/Vyh4rE8KDt4OAAlUrFeXM0TfO8Mc1mExsbG32P8+vkPsydu/8Wncx1vA7Soy8OhRJuMx3y6xXhIZzYDgnaovLeip6vkzxFwJYtCjvFWkKqHnzx4gU2Nzed++Vy2SmCs+/bv4NLpRJ0XY/UqFdkgpZhGCgWi550VNd1z5hqvV5HpVKZ1ykSEcXevXv3cHZ2NvD77oxK0zRUq1Vks1mcn5979vMb9ZqFSPzd2Wg0UCqVsLGxgbW1NaytrTnf29nZcSpWtra2wpn4E0BHFc5NKr039N0gpnQLQHa70I9zm4oQX9MseX6OBMC7sHjszu22YY8L8pw+iYP0ORl1VhUIUTeD6kHDMDwVgzZN09BqtZz7g0a9ZiESmVahUMDl5eXA7xUKhRmfERFRhIRZPagO/pamaSiVSs59ezmSpmmRGfWKRNCauVEZwRSzhbnONQ0T0dOaFLOqIabxs57kObt/9Es5cM39jakiOR+EWXTE0DTNWY9lmia2t7ed0S171CubzYY36jWBxQxaRETka1CrvKiMei1u0LrBX3aRzZZ6xeQ0w7AQmdUsFxX3HGvovNdN/l8K8HB1ggtJJSnDcrD3IIBFDlpERLERrAXTIkjgnyNERJRUC5tpxWaIzxaz0522hRgOdJvXz1+4hgZHnMPEpfPwFmNIKSB6FhQrPUOE6pD/f5M5NAhmWl0LG7SIiGKFy9UAMGjNB7OmiSxcdmWbYol6X3YURv150IdPIXFIZJZFHon/Cfs1033z86t5nxYR0Vhm1TA36hYiaEkpPbel5bs3f+Ipt24iy8K3YJrC5+Umc0/u55AinOdynnPA79N2wB9+4rOsGV0EMuo4PEhEFHUSQCf+AScMixu0mPFE1sJmVb2i+hkN8byi+hIpuhY3aBERxUYyhvbCwKBFRBQHDFoAGLQoIjgk2GPe42Yx+/2Y+CIMcjBoERHFATMtAAxaNEfMrnzMMMMKUq4uZ7xMY5K1zQuRZbF60MGgRUQUeRKQ7OMEMGjRjDG7GmJWGc0kx5njHFs7+T0QaAwMWkREccA5LQAMWkRE0cc5LQeDFk0dhwQDmHeJu4+hhRphnG8Iz7EQRRjkkfifOLu8E1EisGEugAUJWlPp8k6BMMsaIWmd/2f0WpTk/+rqx6AFYAGCFhERJQfntGgqmGEFMMcMK5TrYIXwHELE/y//2UhGlhQGBi0ioqiTADpcXAwwaFHImGEFFJd5rN7zHHDeYV7BuFdnwJMvXOUgMy0AnNMiIqIYYaZFRBQHzLQAMGhRSDgsGNCchwVlkBbqPa3WRw79TfCaxqm/UFisAUCyI0YXf9UQEVFsMNOiG2GGRTQDEpC8NAkABi0ionjg8CAABi26AWZZY4hLiTtFFwsxAHBOi4iIYiTxQYtd3sMnFWZZY0liljXD16SCczmQ0uqIEcYt5hL/q4dd3okoEdjlHcACBC0iIkoOFmLQWDgsOKYIDg2O0ydwmj0FaTwyAUN7YWDQIiKKvGQM7YWBQYsCYYY1pphnKFHqnKQKZhh0jUGLiCjqJLi4uItBi0ZiljWmmGdZY7nJaxUYmtIF6e27UNjGCQCDFhFR5EkAkpkWAAYtGoFZ1hjikBnE4Rx7KFxcTC4MWkREUSflzIYHG40GTNNENptFq9VCsVicyXGD4t/RREQxIDsylNsoe3t7KBaL0HUdp6enME1zBq8uOAYtojDEcNgtzpSJLpcsvF+pj2EYyOVyzv18Pg/DMOZ4Rv0YtIiI4kB2wrkNYQ8L2nK5HOr1+pRf2HgWck6r/baJ1//v12M/7l/+5V/w13/91zN5HI81v8fxWNcEgP/P53FBfnGMOt6tgI/5351/LQP4N+Mfq/BsovMbZNzHPX/+fOxj9Pof//Mv8ac/+b+Ocb1+/Rqbm5vO/XK5jHK5DABotVpYX1/37N9sNkM5bliElIvXG0QIgUle9iwfx2PN73E81vwel9RjxcXR0RHOz8+xv78PAKhWq6jVajg8PJzzmV3j8CAREQEANE1Dq9Vy7jebTWxsbMzvhHwwaBEREQBA13XPcGC9Xoeu63M8o34cHozo43is+T2Ox5rf45J6rDiJ+jqthSzEmNQ//MM/zOxxszzWpOLwuvh+3PxYk+L7EU+FQgGFQmHepzEQMy3i+9GD74cX3w8vvh/zxTktIiKKDQYtIiKKDQYtIiKKDQYtIiKKDQYtIiKKjYUMWix19eL74cX3w4vvhxffj/layJJ3IiKKp4XMtIiIKJ4S2xEjSCuSqLcrCVOQ17q9vY1SqYTNzU1UKpVIdXYOU6vVgmEYOD09dbpZ91qkz0aQ92NRPhs2wzBgmibOz8+Rz+edS3e4LdJnJFJkQhWLReffT548kefn5xPtkxRB349sNit1XU/0e1Gr1eTx8bEsl8sD91mkz0bQ92MRPhtSSnl+fi5rtZpzX9M0z33bIn1GoiSRw4NBLhkdh8tKhyXoa3348CEuLy9Rq9WgadosT3GmdF0f+voW6bMBjH4/gMX5bABWBuXOJHVdR6PR8OyzaJ+RKEnk8KDfJaNrtdrY+yRF0Ndqmiaq1SoA6zo6fkMii2CRPhtBLdJno1gsei7HYZomtra2PPvwMzI/iQxaQS4ZHYfLSocl6Gstl8vO/4ilUinQX+BJtEifjaAW7bNhv9ZWq4Vms9k3X8XPyPwkcngwm83i4uLCs82dygfdJymCvlb3X46apjl/WS+aRfpsBLWon41KpYKTk5O+7fyMzE8ig1aQS0bH4bLSYQnyWg3DwPb29ozPLJoW6bMRxKJ+Ng4ODlCpVDwB28bPyPwkMmgNu2S0/UGLw2WlwxLk/dA0DaVSydmn0WgsXAnvIn42hlnkz4ZhGCgWi84QqF1kwc/I/CW2I8agNRRbW1vY399HoVBYqHUWQd4Pe8jHNE1ompbY98MwDBwfH8MwDOzv7y/8ZyPI+7Eonw3A+n/lwYMHnm3Hx8fQdX1hPyNRktigRUREyZPI4UEiIkomBi0iIooNBi0iIooNBi0iIooNBi0iIoqNRLZxIoojuwmrYRjQdR2FQmHep0QUOcy0iCLANE1n/U+hUMDe3t68T4kokhi0iCJA0zSnS7hpmvjoo4/mfEZE0cSgRdTj4OBgbseuVqs4Pz/HkydPInE+RFHDoEXkcnR01NeOp1qtolKpoFQqwTTNqR6/WCxifX3dc0HBYrHIwEXUxUIMoi7TNHF+fu65wGGr1cKjR49weXl5o8txHB0debqC2+yCC/t72WwWuq6jVCrh/PwcgDV0eHFx4fT9I1pk7D1I1FWpVLC9ve0JDHYT2Wlflfbg4AAXFxfY39+HaZrY2NjA5eWl8/1Wq4VKpeK5DDzRImKmRdTVaDTmlsmUy2UYhuF0XD8+PvZ8P5vNTn1okigOGLQoFo6OjlCv17G/v+/M95yenuLhw4ehrGfyG3prNBo4Pj6GaZo4ODhwrpdkXxjQvjBio9FwhvUajQaazSbq9fpYWVE2m3Xm0gZdl0nTNDQaDa7fooXGoEWRZxgGPvnkE9TrdZRKJc9Q3eHhoSc4NBoNlEolZ5gvaACxr4vkVigUnOILdzXf9vY2KpUKNE2DpmnY3NzE2toaarWa51pU1Wo11Gss5fN5mKbJoEULjUGLIk/TNGSzWZydneGrr75ytg8KNAA8QWZra8vpMjGIaZrI5/OBzsc+pp2Z9d63t7mvbBsGDhESseSdYsAOBr1Zxtdff42trS3Pvn7DfKZpYnNzc+Rx/Kr7Rp2TWy6X8/13mHqDNNGiYaZFsWAYhifwtFotmKYJXdc9gcowDGSzWedS6KZpolarjfxln8vlnBLzqGq1Wix5p4XHTItioVarebKqb775Bp988gkAeBbi1ut17OzsoFAooFgsOuubRtE0baxMa9TQX9hDgwBwfn7OoEULj0GLYsE0TU9Rg67ryGazODo68iwGNgzDM4T40UcfBariKxQKfcGt0Whgf38fZ2dnODg4QKvVcraZpuksGK5UKgCsqkJ7e6PRwOHhoSeg3hSLMIi4uJgSpNVq4cGDB6jX6862UqmEra0tT2AbZHt7G/v7+5GcN+LiYiILgxYlgn1pD9M0nfVTp6enyOfzgQKW/RyHh4fY39+f5qlO5ODgAMVikcODtPBYiEGJoGlaXxYy7hopTdOwvr4euR5/dml/lM6JaF44p0Xk8uTJkxs1xp2GarUaOFskSjoODxIRUWww0yIiothg0CIiothg0CIiothg0CIiothg0CIiothg0CIiotj4/wGdo3ojLBsxgAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if plots:\n",
    "    pl=o2sclpy.plotter()\n",
    "    pl.colbar=True\n",
    "    pl.xtitle(r'$ n_B~(\\mathrm{fm}^{-3}) $')\n",
    "    pl.ytitle(r'$ T~(\\mathrm{MeV}) $')\n",
    "    pl.ttext(1.25,0.5,u'$ A $',rotation=90)\n",
    "    pl.den_plot(t3d,'A')\n",
    "    plot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5869c8e6",
   "metadata": {},
   "source": [
    "For testing purposes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "cfee5b0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_fun():\n",
    "    assert numpy.allclose(t3d.get(0,0,'A'),81,rtol=0.1)\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "154fe708",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}