{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "195dfb1f-429b-4672-bbfe-53e6c77926c6",
   "metadata": {},
   "source": [
    "# Complete Example\n",
    "\n",
    "The subcomponents are all joined together into the main flasher class, which combines the superancillary to first check the phase, followed by K-D tree lookup for the best starting point, and finally Newton iteration to get the correct temperature and density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "09dd9f4c-c7fd-4dbf-8492-abf5cbf08044",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-06T11:32:38.895160Z",
     "iopub.status.busy": "2025-01-06T11:32:38.894745Z",
     "iopub.status.idle": "2025-01-06T11:32:39.185857Z",
     "shell.execute_reply": "2025-01-06T11:32:39.185578Z"
    }
   },
   "outputs": [],
   "source": [
    "import timeit, json, functools\n",
    "\n",
    "import teqpflsh, teqp \n",
    "import CoolProp\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1aeec3d9-8574-47f1-9a58-d116c82d2e4f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-06T11:32:39.187406Z",
     "iopub.status.busy": "2025-01-06T11:32:39.187304Z",
     "iopub.status.idle": "2025-01-06T11:32:40.507482Z",
     "shell.execute_reply": "2025-01-06T11:32:40.507229Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10bae2210>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGdCAYAAADJ6dNTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAX6VJREFUeJztnQd8VVW2/3+3l/TkJkBIoYdOABECSpFmQGzYpVgRu+hDxjfqKIzlMfOc58xzfM78bTPi6DijjiKKoIING4o0B6WDQCCUtHuT287/s09IOAnJzb0hnHLz+84czz5r/e7Ovockd2WvtfcxSZIkgRBCCCFEJ5i1HgAhhBBCiBIGJ4QQQgjRFQxOCCGEEKIrGJwQQgghRFcwOCGEEEKIrmBwQgghhBBdweCEEEIIIbqCwQkhhBBCdIUVBiQcDmPfvn1ISkqCyWTSejiEEEIIiQKx72tFRQWys7NhNpvjKzgRgUlubq7WwyCEEEJIK9izZw9ycnLiKzgRMyZ1by45OVnr4RBCCCEkCsrLy+XJhbrP8bgKTupSOSIwYXBCCCGEGIuWSjJYEEsIIYQQXcHghBBCCCG6gsEJIYQQQnQFgxNCCCGE6AoGJ4QQQgjRFQxOCCGEEKIrGJwQQgghRFcwOCGEEEKIrmBwQgghhJD4CU4ef/xxeZe3u+66q942duxY2aY85s6d2+B1u3fvxtSpU+F2u5GVlYX58+cjGAyeylAIIYQQEie0evv6r7/+Gs888wwGDhx4ku/GG2/EwoUL669FEFJHKBSSA5OOHTvi888/x/79+zFr1izYbDY8+uijrR0OIYQQQtpzcFJZWYmrr74af/7zn/HrX//6JL8IRkTw0RTvv/8+Nm/ejJUrV6JDhw4oLCzEokWLsGDBAjz00EOw2+3Q+nHOvkBI0zEQQghp37hslhafPxPPtCo4ufXWW+XZjwkTJjQZnCxZsgQvvfSSHKBMmzYNDzzwQP3syZo1azBgwAA5MKlj8uTJuPnmm7Fp0yYMHjz4pP5qamrkQ/lUw9MVmFzyf2vw9db9MNucp+VrEEIIIS1R2MmFN+4Y124DlJiDk1deeQXffvutnNZpiquuugr5+fnIzs7G+vXr5RmRLVu24PXXX5f9Bw4caBCYCOquha8pHnvsMTz88MM43YgZExGY7Hniktpb48mDNSMHVmcirO5UWBPTZJ3ZXhtomZ1158ST+qqzhasrG/gt7mSEvOUnXYuzoLm2kqbsVncKgt6yZt9bY3+k65bakWxteRbE4mvJFk27OX+s/bT0mlh8Wly3VqOWTW1tJLtafenNZwS/UTXhQDW+3noMh8sq4Ek9+TOgPWCSxHRBlOzZswdnnHEGVqxYUV9rIgpgRWrmf/7nf5p8zYcffojx48dj69at6N69O+bMmYNdu3Zh+fLl9Rqv14uEhAQsW7YMxcXFUc2c5ObmoqysDMnJbfcP5/UHkeCwNeO1A45kwO4C7DbxqQzYbIAzEWa7C2bbiXSUmHWxOmoDl7AkwWwy1WocbphNZoSlMMI13vpr2J0wm2trk8Ph8EltZXAjgp3GdnEW13K/vsqTRt7YH+kaEiK2Tc5ESGIMTdja8qy8B9H4WrJF027udXVBS7T9tPSaWHxaXAtao1HLprY2kl2tvvTmM4LfqJo6XajGhy9/NRV52U2XSBgV8fmdkpLS4ud3TKt11q5di4MHD2LIkCGwWq3ysXr1avz+97+X26LYtTHDhw+XzyI4EYhUT0lJSQNN3XVzdSoOh0N+E8pDffxATSVQWQ6IWQNvKeA9BpSXIlx5BEHF4a8oRbDqmHz4y0pq25VHEaw8jLDfK58bXFeUIlztlY+m25WKo2m7oKGuMqI/0nXEtq+8eVsbn2P1tWSLpt2cvzV9tpVPi+vWatSyqa2NZFerL735jOA3qkYgBU/8Qd4eiSk4ETMgGzZswLp16+oPMZMiimNF22KxnPQaYRd06tRJPhcVFcl9iCCnDjETIwKOvn37nvo7IoQQQkj7qTlJSkpC//79G9hEOiYjI0O2b9u2DS+//DKmTJki20TNybx58zB69Oj6NNCkSZPkIGTmzJlYvHixXGdy//33y0W2YoZE31TLqQxUS4BVAuoKlfxWIKyYNQqFELbVvpew34ewoqBJstoRrqmWc4r110E/JOsJfeO2FAyceH0ocJJdCvprzwF/A239axr5I13Xv6aZthirxeZs0taW50h9N+VryRZNO9LrYuknmtfE4tPiurUatWxqayPZ1epLbz4j+I2qiUUXr7R6n5OmEMuAxRJhUX9SVVUl14VMnz5dDj7qELMrS5culVfniFkUEdzMnj27wb4o+g9QAkDg+N4tThGkmIGwYgl0OIxgtRcQMUkwgKAk6kskWEwWmLwVCAUDCAWq669FgBKqqZL14bB0ctty4p8p5K86yW6S3HKKRYLU5Dd0Y3+ka/H/SG2T3dWsrU3PgdopTZPFFp2vJVs07QivE+eo+4niNbH41L5urUYtm9raSHa1+tKbzwh+o2qUOp/Pi/bKKQcnq1atqm+LYETUoLSEWM0jil+NmQVzAFZREGsT0RhgdQI2URyruJVWO6zHV/IE4ZWLY81WB8w2G6zuJEjecpgkZ/11SPzV7kiQ9eFg4KS22aoIfOwJJ9nFUjNR1GoKVDfUHqexP9K1/C4jtMWMT3O2Nj0fn3kKecui87Vki6Yd4XXiHHU/UbwmFp/a163VqGVTWxvJrlZfevMZwW9UjVLncp3YwLS90aYzJ/GPHTDZawMROUBx1h52B6BYrSOCExGMCMxitY4cmDhhtjvlaFic5ZUudddmi3yW9XbXSW2T+FrHMQVtJ9lNdUGKCDSsTcycNPJHupbTVi21I9na8ByrryVbNO3m/K3ps618Wly3VqOWTW1tJLtafenNZwS/UTWx6OIVBicxwZoT1pxE324rnxbXrdWoZVNbG8muVl968xnBb1SNgDUnJEZYc8KaE9acaGlTWxvJrlZfevMZwW9UjVLna8c1J6f0VOL2ebtcgDUBcNoBl0jnuGprThzOBoeoOZE3YrPaats2JyzHa0zks+LabDbJNSQWe0KTbVEXUnc0Za+vGRFnhbY5f6TrltoIBZq1tenZJlJhjuh9LdmiaTfzOvlfPpZ+WnhNLD4trlurUcumtjaSXa2+9OYzgt+oGqXOxZoTEh2sOWHNSfTttvJpcd1ajVo2tbWR7Gr1pTefEfxG1QhYc0JigDUnrDmJvt1WPi2uW6tRy6a2NpJdrb705jOC36gaAWtOSIyw5oQ1J6w50dKmtjaSXa2+9OYzgt+oGqXOx5oTEh2sOWHNCWtOmtOoZVNbG8muVl968xnBb1SNUudizQmJDtacsOYk+nZb+bS4bq1GLZva2kh2tfrSm88IfqNqBKw5ITHAmhPWnETfbiufFtet1ahlU1sbya5WX3rzGcFvVI2ANSckRlhzwpoT1pxoaVNbG8muVl968xnBb1SNUudjzQmJDtacsOaENSfNadSyqa2NZFerL735jOA3qkapc7HmhEQHa05YcxJ9u618Wly3VqOWTW1tJLtafenNZwS/UTUC1pyQGGDNCWtOom+3lU+L69Zq1LKprY1kV6svvfmM4DeqRsCaExIjrDlhzQlrTrS0qa2NZFerL735jOA3qkap87HmhEQHa05Yc8Kak+Y0atnU1kayq9WX3nxG8BtVo9S5WHNCooM1J6w5ib7dVj4trlurUcumtjaSXa2+9OYzgt+oGgFrTkgMsOaENSfRt9vKp8V1azVq2dTWRrKr1ZfefEbwG1UjYM0JiRHWnLDmhDUnWtrU1kayq9WX3nxG8BtVo9T5WHNCooM1J6w5Yc1Jcxq1bGprI9nV6ktvPiP4japR6lysOSHRERZJDkBkcEzHCy7qUF42cp1kb3xuK6QY/ZGum2lLEWxtfY7V15Itmnakf8ZY+2wrnxbXrdWoZVNbG8muVl968xnBb1TN6fh4MBoMTmKGaR2mdZjW0dKmtjaSXa2+9OYzgt+oGqXOx7QOiQ6mdZjWYVqnOY1aNrW1kexq9aU3nxH8RtUodS6mdUh0cCkxlxJH324rnxbXrdWoZVNbG8muVl968xnBb1SNgEuJSQxwKTGXEkffbiufFtet1ahlU1sbya5WX3rzGcFvVI2AS4lJjLDmhDUnrDnR0qa2NpJdrb705jOC36gapc7HmhMSHaw5Yc0Ja06a06hlU1sbya5WX3rzGcFvVI1S52LNCYkOLiXmUuLo223l0+K6tRq1bGprI9nV6ktvPiP4jaoRcCkxiRGmdZjWYVpHS5va2kh2tfrSm88IfqNqlDof0zokOpjWYVqHaZ3mNGrZ1NZGsqvVl958RvAbVaPUuZjWIdHBtA7TOtG328qnxXVrNWrZ1NZGsqvVl958RvAbVSNgWofECNM6TOswraOlTW1tJLtafenNZwS/UTVKnY9pHRIdTOswrcO0TnMatWxqayPZ1epLbz4j+I2qUepcTOuQ6GBah2md6Ntt5dPiurUatWxqayPZ1epLbz4j+I2qETCtQ2KEaR2mdZjW0dKmtjaSXa2+9OYzgt+oGqXOx7QOiQ6mdZjWYVqnOY1aNrW1kexq9aU3nxH8RtUodS6mdUh0MK3DtE707bbyaXHdWo1aNrW1kexq9aU3nxH8RtUImNYhMcK0DtM6TOtoaVNbG8muVl968xnBb1SNUudjWodEB9M6TOswrdOcRi2b2tpIdrX60pvPCH6japQ6F9M6JDqY1mFaJ/p2W/m0uG6tRi2b2tpIdrX60pvPCH6jagRM65AYYVqHaR2mdbS0qa2NZFerL735jOA3qkap8zGtQ2K7ZccDEZsFMJkAs7U2QKk7hMpqhVkED+FgbVukQaQQLA6XfFZeI1hTm4Kx2JtsmxT/a8pen5IIBhpom/NHum6pHa6patbWpufj6Y6ofS3Zomk38zr53zOWflp4TSw+La5bq1HLprY2kl2tvvTmM4LfqBqlzsW0DokttSNmJ6xAIARYLHIAUheU1KuCwiZHKbVtiw0mkwWhGh8kkwVhSaq/htVRO+Mh9E20xV/z9f2G/Cfbj//VL+pb5BmNxjT2R7iumw1prm12JDRra9Pz8b8sova1ZIumHeF14hxrn5FeE4tP7evWatSyqa2NZFerL735jOA3qkap83HmhEQHC2JZEMuC2OY0atnU1kayq9WX3nxG8BtVo9S5OHNCooMFsSyIjb7dVj4trlurUcumtjaSXa2+9OYzgt+oGgELYkmMsCCWBbEsiNXSprY2kl2tvvTmM4LfqBqlzse0DokOpnWY1mFapzmNWja1tZHsavWlN58R/EbVKHUupnVIdDCtw7RO9O228mlx3VqNWja1tZHsavWlN58R/EbVCJjWITHCtA7TOkzraGlTWxvJrlZfevMZwW9UjVLnY1qHRA/3OeE+J9znREub2tpIdrX60pvPCH6japQ6F9M6JHq4zwn3OeE+J1ra1NZGsqvVl958RvAbVaPU+ThzQqKDBbEsiGVBbHMatWxqayPZ1epLbz4j+I2qUepcnDkh0cGCWBbERt9uK58W163VqGVTWxvJrlZfevMZwW9UjYAFsSRGWBDLglgWxGppU1sbya5WX3rzGcFvVI1S52Nah0QPC2JZEMuCWC1tamsj2dXqS28+I/iNqlHqXEzrkOhhQSwLYlkQq6VNbW0ku1p96c1nBL9RNUqdjzMnJDpYEMuCWBbENqdRy6a2NpJdrb705jOC36gapc7FmRMSHSyIZUFs9O228mlx3VqNWja1tZHsavWlN58R/EbVCFgQS2KEBbEsiGVBrJY2tbWR7Gr1pTefEfxG1Sh1PqZ1SPSwIJYFsSyI1dKmtjaSXa2+9OYzgt+oGqXOxbQOiR4WxLIglgWxWtrU1kayq9WX3nxG8BtVo9T5OHNCooczJ5w54cyJlja1tZHsavWlN58R/EbVKHUuzpyQ6OHMCWdOOHOipU1tbSS7Wn3pzWcEv1E1Sp2PMyckejhzwpkTzpxoaVNbG8muVl968xnBb1SNUufizAmJHs6ccOaEMyda2tTWRrKr1ZfefEbwG1Wj1Pk4c0Kig5uwcRM2bsLWnEYtm9raSHa1+tKbzwh+o2qUOhdnTkh0cBM2bsIWfbutfFpct1ajlk1tbSS7Wn3pzWcEv1E1Am7CRmKEm7BxEzZuwqalTW1tJLtafenNZwS/UTVKnY9pHRI9LIhlQSwLYrW0qa2NZFerL735jOA3qkapc7XjtM4pBSePP/64XHdw11131duqq6tx6623IiMjA4mJiZg+fTpKSkoavG737t2YOnUq3G43srKyMH/+fARF0aihCmJRWxAriemDICCFTxzHC2LDoWB9QawogJUaFcTWXdcVuCqLXZVtSfG/puyNC1wb/y+mgtgW2g0KYhvZ2rogtr54LBpfS7Zo2s28Tv73jKWfFl4Ti0+L69Zq1LKprY1kV6svvfmM4DeqRqnzceYkdr7++ms888wzGDhwYAP7vHnz8Pbbb+O1117D6tWrsW/fPlx88cX1/lAoJAcmfr8fn3/+OV588UW88MILePDBB2EMOHPCmRPOnGhpU1sbya5WX3rzGcFvVI1S5+LMSWxUVlbi6quvxp///GekpaXV28vKyvDss8/iiSeewDnnnIOhQ4fi+eefl4OQL774Qta8//772Lx5M1566SUUFhaiuLgYixYtwlNPPSUHLPqHMyecOeHMiZY2tbWR7Gr1pTefEfxG1Sh1Ps6cxIZI24jZjwkTJjSwr127FoFAoIG9d+/eyMvLw5o1a+RrcR4wYAA6dOhQr5k8eTLKy8uxadOmJr9eTU2N7Fce2sGZE86ccOZES5va2kh2tfrSm88IfqNqlDpXO545iXm1ziuvvIJvv/1WTus05sCBA7Db7UhNTW1gF4GI8NVplIFJnb/O1xSPPfYYHn74YegDbsLGTdi4CZuWNrW1kexq9aU3nxH8RtUodT7OnETHnj17cOedd2LJkiVwOp1Qi/vuu09OGdUdYhzawZkTzpxw5kRLm9raSHa1+tKbzwh+o2qUOhdnTqJDpG0OHjyIIUOGNChw/fjjj/G///u/WL58uVw3cuzYsQazJ2K1TseOHeW2OH/11VcN+q1bzVOnaYzD4ZAPfcCZE86ccOZES5va2kh2tfrSm88IfqNqlDofZ06iY/z48diwYQPWrVtXf5xxxhlycWxd22az4YMPPqh/zZYtW+Slw0VFRfK1OIs+RJBTx4oVK5CcnIy+fftC/3DmhDMnnDnR0qa2NpJdrb705jOC36gapc7FmZPoSEpKQv/+/RvYEhIS5D1N6uzXX3897r77bqSnp8sBx+233y4HJCNGjJD9kyZNkoOQmTNnYvHixXKdyf333y8X2epndiQSnDnhzAlnTrS0qa2NZFerL735jOA3qkap83HmpO343e9+h/POO0/efG306NFyqub111+v91ssFixdulQ+i6BlxowZmDVrFhYuXAhjwJkTzpxw5kRLm9raSHa1+tKbzwh+o2qUOhdnTlrPqlWrGlyLQlmxZ4k4miM/Px/Lli2DMeHMCWdOOHOipU1tbSS7Wn3pzWcEv1E1Sp2PMyckejhzwpkTzpxoaVNbG8muVl968xnBb1SNUufizAmJHs6ccOaEMyda2tTWRrKr1ZfefEbwG1Wj1Pk4c0KihzMnnDnhzImWNrW1kexq9aU3nxH8RtUodS7OnJDo4cwJZ044c6KlTW1tJLtafenNZwS/UTVKnY8zJyR6OHPCmRPOnGhpU1sbya5WX3rzGcFvVI1S5+LMCYkezpxw5oQzJ1ra1NZGsqvVl958RvAbVaPU+ThzQqKHMyecOeHMiZY2tbWR7Gr1pTefEfxG1Sh1Ls6ckOjhzAlnTjhzoqVNbW0ku1p96c1nBL9RNUqdjzMnJHo4c8KZE86caGlTWxvJrlZfevMZwW9UjVLn4swJiW3mxFd763wmwC4B5irAovhmk+wIV1dBnv4IBBA2VSHskGAxmRH0ViAUCiAUqKm/DpvMCNXURsjhsHRyWwQ5xxHXje0mySX/1S9JUu3MSiMa+yNdy183QrtulqUpW5ue6/6yiNbXki2adoTXyecY+4z0mlh8al+3VqOWTW1tJLtafenNZwS/UTVKnY8zJyR2rIo7KHIwjTHVmsXMSpP+pl/SZLs5TSy+pvyRrptpmyLY2vocq68lWzTtSP8EsfbZVj4trlurUcumtjaSXa2+9OYzgt+omlh08QqDk1bdMlHrYQNcDsBmA2xuwOpUHHaYnW6YHcJuldtWmwMWiw1Wd1LtWXFtNptgcbhhsbubbItUTt3RlN1kMskpibq0T+OjsT/SdUtt6Xgqpylbm56Ppzui9rVki6bdzOvkf/VY+mnhNbH4tLhurUYtm9raSHa1+tKbzwh+o2qUOhfTOiR6WBDLglgWxGppU1sbya5WX3rzGcFvVI1S52Nah0QPC2JZEMuCWC1tamsj2dXqS28+I/iNqlHqXJw5IdHDmRPOnHDmREub2tpIdrX60pvPCH6japQ6H2dOSPRw5oQzJ5w50dKmtjaSXa2+9OYzgt+oGqXOxZkTEj2cOeHMCWdOtLSprY1kV6svvfmM4DeqRqnzceaERA9nTjhzwpkTLW1qayPZ1epLbz4j+I2qUepcnDkh0cNN2LgJGzdh09KmtjaSXa2+9OYzgt+oGqXOx5kTEjvchK1xk5uwcRM2tWxqayPZ1epLbz4j+I2qEXATNhIj3ISNm7BxEzYtbWprI9nV6ktvPiP4japR6lxM65DoYUEsC2JZEKulTW1tJLtafenNZwS/UTVKnY9pHRI9LIhlQSwLYrW0qa2NZFerL735jOA3qkapc3HmhEQPC2JZEMuCWC1tamsj2dXqS28+I/iNqlHqfJw5IbHDgtjGTRbEsiBWLZva2kh2tfrSm88IfqNqBCyIJTHCglgWxLIgVkub2tpIdrX60pvPCH6japQ6F9M6JHpYEMuCWBbEamlTWxvJrlZfevMZwW9UjVLnY1qHRA8LYlkQy4JYLW1qayPZ1epLbz4j+I2qUepcnDkh0cOCWBbEsiBWS5va2kh2tfrSm88IfqNqlDofZ05I7LAgtnGTBbEsiFXLprY2kl2tvvTmM4LfqBoBC2JJjLAglgWxLIjV0qa2NpJdrb705jOC36gapc7FtA6JHqZ1mNZhWkdLm9raSHa1+tKbzwh+o2qUOh/TOiR2mNZp3GRah2kdtWxqayPZ1epLbz4j+I2qETCtQ2KEaR2mdZjW0dKmtjaSXa2+9OYzgt+oGqXOxbQOiR6mdZjWYVpHS5va2kh2tfrSm88IfqNqlDof0zokdpjWadxkWodpHbVsamsj2dXqS28+I/iNqhEwrUNihGkdpnWY1tHSprY2kl2tvvTmM4LfqBqlzsW0DokepnWY1mFaR0ub2tpIdrX60pvPCH6japQ6H9M6JHaY1mncZFqHaR21bGprI9nV6ktvPiP4jaoRMK1DYoRpHaZ1mNbR0qa2NpJdrb705jOC36gapc7FtA6JHqZ1mNaJPa1jggSHRUJq4DBSTX5Iwf0IWY8iIeiDGWHUWMvhDpYj2ZQKb+gIaqxlSAyVItlce11tLUNSaD+Szcnwho7Cay1HUmgn0sxJqAgfg89WAXd4K5LMaTgWKke51Qt3aA8SzKkoC1Wg0loJZ+gAXOZ0VITK5Wt76CgSzOkoDx5DpbUKrmAlHCYPQjDDGzyKMqsPCFghwYRjgUrUWKxM6zCtYxi/UTVKna8dp3UYnLSaOE3rSJHbDdI6UvtK67hRg2xTCVJQCZupBJ0sJqTAC5PlELKkIDLt1XLAkWEKo6/djCOOCgQlHxLNIaTZA6hEEIlOK/IcFux2hHDIH0Cmwyb3XSkFkeg47rOHEEoJwWKxIM8exbUthFByCP6wGQV2a+yvd4QQsh6/dljk8dSPz2lDntOC3TUhefwWRwJcZifMNie22U3Y5xehSyJSbJk4JiVihxXYHg7CaukEyZyNMiRin8mPEoRRiZQG95NpHWP7jOA3qkbQ3tM6DE5iRkQjYpbEWpvWsVpr0zpmxa201KZyBOFwqDbFI1I9x9M44VAAJlu4/jpY45XTNIJQje+ktki91GN3n2SvT794yxtqj9PYH+la/EREagd95bC4kpu0teW5btpTTG1G42vJ1nLbjw4uE9J929HfaUKOqRRB7270sdWg0O6D31GK6mAFPA4r3DYTDiGIBJulth0OwpNggS8IVDlDcjshZEbnBGB3IIjcFOEzQZIs8LjNgNUFd6INGS4L3E473A4bpKoQPMkOwGFDXqYFu4/6kZeRAJjNyMs2Y/cRP/IyEwGTGXkm4a9BnicBMFmQZ7Zgy/4K9OqYBG9NAHndTdh9qAJ56U5ACiMvHMLu0irkpdtrrzuGsPtwFfLSHEA4hLwMcV2NvFSL7BeHCMQyrFVwW4IAAvK4xfgT7H5k2ALwBSvgdEjICtS+94TwDvle7LaGEEqvDXQ8JrNsK00MY3cggI6JTlSZ01BuzcB6mxNbYIXFmY8Kc2fslTKx1ebG/pCfaR0D+IzgN6pGqXMxrUOih2kdI6d1LBYr8oO70MO0Dx77VgxCOTzBEqQ4S+APBNAlwwR/wAyXzYRKdxAJZgsSqs2QbBJKbSa4XImAOx3utERkZHaEyZ2BTHMi3KlZyEzviNJqC9wpHrhT0gGbC3kBkxy8elIy4A4C7uQMOeDwiMGUlsLjkVvI9Hrhdp/4RZTX6N+wpesCAF5FH7G+vvF1g/GFQ3AHfMgsK4XbBnjLjsBtk+D2e4GSA3Cb/XCbfID3KPJ8R7F79y54bNVwBcvgLT8ChEuRk3IE1cEgergP47CvFB5zGGNTw/AHvkZBihXegIQ94RD86Q7sCGbCj1wcsGZjfdiJLVYXDgQSUaOYeWFaR1ufEfxG1Sh1PqZ1SOwwraP3tI4NIfQ0/Yx880YMsx1BbnAP3Ka9yDbXqkLOECxhC5x2oNQfREaiHcfsHdCroDd87s6Q7BlI8OTDk98bXmsaTJZkuNOy5ADA7a39pSG3FV+/NtQ4gdLX+G+gusCkrp9TpS36aHJ8ZgvgSIQ7K7H266R1OaHpevLr6gKdumBJjKq0ZD/yHEF4D+6AJ3AE7kO7sGfbv9ErxY/SfbvgqtqLdP9e7C3zoSjlZ7jN+2CSTNiDEELu2pmYBEdHfF+TjR+lHHxvzsBGUyIOIqnB12ZaRz2fEfxG1QiY1iExwrSOXtM66ahA1/A6nJF4EJ6KnzDRcRh2UwClgSAS6oOQECRLMjLyCnDYmYfcfsOB9O4wSSnwdBso/7Ui8BwPQOo+8N2NA402DgTiEeU98nToVGtLza09i9mesbUBjLjX4pwQ8KHg6G7g6A7g6E7gyA7kHNmOH3/YiF6uI/D5DmKs9SDGYh32msKosoeQ7HChxNQF68PdsNaRhy/KHThmy2owDqZ1mNYxkkapczGtQ6KHaR29pHXclhD6+ddisH07hvh2ob9jP0JJIewql1CQbkKZ14T01FS4uvdFeq8RMHXsB1NiPjzdBsmplQTFB2jjtIbcPwOQ00598Cef3UBKBtBlcL1fBC0iZQXfUXgqd8G761t496xH+sEtCO/cDBu86BX8CYW2rSjyhXB3chjloRQcQU98F+6Jby0D8E0gmWkdpnUMo1HqfEzrkNgJAlIQCJqBUODkSbhg4PjZL38QwuJH2GyBZLUjLGxBxXU4DMlWO+MRDlSf1BYRdB1SKHiSXQqeOCu19a9p5I903fg1J7fF+7E0Y2vDczN9dwrtxyTLtzjLsgHu8h9gs5pgl8JIt1kgSSZUJ/fAGWedhZ+tXVBQNBXI6Aavz1ebXmCwYTjq/83EOaMz3Pkj5VkXEbSYqqrgri4B9n0H7/Yv0Xnfd9j94wZ0clSie2gdxti+hy/4D+w12/ApOuFnyyBsMQ3DeqkbxHeY8vtMSax2qNSX3nxG8BtVE4suXmFw0mqcgMkK2K2ASKUo0ymibbUfr8uQatsirWO1w2S2yueQSPfUXVvMMIlZEJEmsTlPapuOpxoEJrGpWyO7yXr8bLPB5D+hPfGahv5I1/V1JM20xZhFgWRTtjY9H8dstaGHtAdnOzajf3AjuloPwQMbfNUSbEmAPzEXrm7D4e55FjyDiuE1J8gfaIWKQMSdIOZISDxxItjMBHL7wz18phyw5JcfhbtiJ/DzWni3fgrX/i/hrDqEcxK2wYKdyDUvRZXkwD8CufjS3gvrcSaOoFEaSPH9F429Na9pTV968xnBb1RNLLp4hcFJzLDmRI2ak+72ozjX8hW6lK9CN2cZvP4Q0l0WhCQbpM5nomfRRfjJ3B0FZ4yT36OyNoS084ClY2eg5yi4x94Bb2Ul8g5twZ6vlyLHtxlVP30Kb9UxFCf8G72ObERP51LsCebgG0shPjOfgS9tHRFgzUmLPiP4japR6lysOSHRw5qT01VzkmgJYkzgA0xxfoNe1SXwBSUkJ4Wx3+dE9uDRsPScCFfvc+Dp3E3+2srZEUKawp2YCCQORUHXobVFt+EwEko2Ys+a19F972eo2vM9err2oFd4Ly7yL0VJyI1l5m7YEDob3+MsBI7/imTNCQznN6pGqfOx5oTEDpcSt81SYgl9TLsx2foVCgPrYEYQySYJIYsd3YeMxbbkESg8Zwa8ITPy8poqWyUkxvqVxBEo6D4CpaWlyEI1vJuWw/vDSqSXfIK9u0pxScJ3uCC4Hkmhv+ALy1C8FzoDnyBXVJk1C5cS69NvVI2AS4lJjDCt0xZpHafLjZE1n2Jm8jfA0V1w2y1IcoZh9vSEqd+FyBh+OcwZ2ShS7AVCyOnYx8U95npgzPUoPViCwtL1+HHli+hc+ikk7yGcVb0ak2wfY0vIiWXBHlhvL8ZGUeei+OhgWkeffqNqlDoX0zokepjWOZW0TkL1QRTbvsFY/xp0tJRhT4mEM3ISUJE/Eaa+FyJn6CS5gJUra4jaeLI6AFkTUdh3IryVFfD++DEStryLqo3vwC3tx9XJ32Fa5fdw2DrhffPZWGoZi71SFtM6OvUbVaPU+ZjWIbHDpcQn25o/p5m9OF9ajZHhtUgKVqMmbIKUkoFBE2Zhj+csDB89ngEJ0Q3uxCS4h0wFhkyFu/JxJGz5CD+teA7dyj6Dz3sAV5r+gatM/8QOZ288FyrE+0GxVW7qSf1wKTGXErdGI+BSYtJKuJQ4mqXEqeEyTHd8hjNrvoAt7IPTZoYzqyc6n3UDqnPPRmaXrujD1A3ReVGte+g0eIZOk7faN69/G9KmN+HY/wW6+H7AXdiEy20ufOYfhk9Ml2O31KH+tVxKzKXErdEIuJSYxAhrTqKpOUl12XCefxmm2b/G/qMVyPPY4EvrD/eoOXAXnIPcvDzOlBDD4cnOB7Jvg3f0dfDu34qEjf8EvloC7749uDPtU8yq+BRbHP3wpnkiPrD1abNlyXrzGcFvVI1S52LNCYke1pxEqjmRqstxnm0tLvB/gJD3KLw2Mwb2H4hDva5C99GXwZOZyaCExMeeKt0HAt0Hwn3OPUjY8A4Of/ocUg+uwcDARgwybcL5Fcl40zQcnwbORbCFJyobzWcEv1E1Sp2PNSckdlhz0tg2JLQRM8zLkV59AOV+YECvXjjUZzaCAydj3OATz0shJJ5wJyXDPfJKeEZeidKf1sL81V8QXv8a8qxHcEP4HUyv+hjbkyfhZdNUlCBdd7UjrDnRn0bAmhPSSlhzUmfLsnoxx/wWcivXw2Q2Izk1DZ3OnoND2eMwbvyEtr/1hOgUT8+hQM+h8E59AOZPnsPW957Bmc796Fn5Fs6zvYuPLcPxrGkcfrQ6DVFXonXNCGtO2i8MTmKGNSd1NSdmUxiXJ3yLMWVvIRCoQW6aHYG+l8I8ZCby+g+u30eCkPaGO9WDvGn3wjPuFuz5+CW4Nv4N0s9fY6rpM/Sq+RRvV3XDV2mzsb3xc310VFeidc0Ia04CrDkhscCaE3HOwUHMxlvIOHAQWSmAu3shDva7AcMmXcKghBDFSp+CKXPhHTsL3q1rcHjlE8DhlZiRsQVjSu9HecfheN50Ef5t7aW7uhKta0ZYc2JjzQlpDe2z5sQUrMFFtjXo730HFoTQr1sGDhbMRNakORjap88p31VC4raAduB4eAaOR8bmT/H1X3+FIclfwhf4BoOC32CtZSD+Yr8cG5Cqq7oSrWtGWHPSfmFwQqKmg+kIbre9BvuxHchOl1DmOROBCx7E2KJxXIFDSJR4+p6F4sc+QOkPn6NmxRNwbHsPo7AeQ33r8S/fAPzVNBZHGqV7CGlvMDhpNe2rIHaEaRNuM32AY2UVGJSfjl09ZmPApBtR0Lv36brBhMQ1nj4j5aN0y5c4/N5iOLe/i4m279Hbtw7f+8/C27gOhxrtOsuC2Nj8RtUIWBBLYqR9FcQ6nU7M9L+Kfr61gMOKQUOKcHDw7Rg7/nzOlhDSBngKhsNT8E+UblyFY28uRM7OzzAoYQ0mHvkWHyVOwhL7paiCiwWxMfqNqlHqXCyIJdHTfgpiO5oOY37Ni6ip2I2MJDNSRsxG+dBrMG7UqNN9kwlpd3j6j5WP3Z/9A0fefRzdrJuR6nsLY2pWY4n9EvwrMIQFsTH4japR6nwsiCWxE98Fsf3DP+LimrdgCvoxpFsmtvW8EYUX3oS8vLw2vYuEkIbkjboEeSOnY8uy/4PpsyfRuWo37g0/h77H3sGS8ATsShjZ5OtYEBs/GgELYglpgIQLLZ+gr+8tJFqBnn2G4vCwebhgykVM4xCiFiYTCqbeDO/ZV6J0xf+g9MOnUJSwD10qnsNm3/d4M+VW7AOX7JP4LqCImqeffhoDBw5EcnKyfBQVFeHdd9+t948dO1auX1Aec+fObdDH7t27MXXqVPmDLisrC/Pnz0cwGERcFMQqD1HMJGZARG6xiYLYBtd1Ba7KYtdGha/1R1N2ZYGrUtucv5lrm82COa4VGHnsHfT1mOHqMwEJlz6N4kuuZmBCiAa4k1ORN/0h9P7V17D0vxBpbhsuS9mIpyvuwoyav8MqZnCPw4LY+NEIWBAbAzk5OXj88cfRs2dPuSbhxRdfxAUXXIDvvvsO/fr1kzU33ngjFi5cWP8a5YdaKBSSA5OOHTvi888/x/79+zFr1izYbDY8+uijMAbxWRCb4jRjbuAVhA9vR78uLpT0vRZDp9zO1TiE6AB3VhcU3vUqdn/9HrzLf42UQ9/ihvBrGHr4Y7yQchPW2vJYEBsnGqXOxYLY6Jg2bVqD60ceeUSeTfniiy/qgxMRjIjgoynef/99bN68GStXrkSHDh1QWFiIRYsWYcGCBXjooYdgtxshUoy/gth0KYiFoRexvWQ3RnVNxu4h83DGedexvoQQnZE37Fx4+56N0o+fRem7j2Jkcgm6lD+Mf1QNxiuWqfCxINbwGqXO144LYmNK6ygRsyCvvPIKqqqq5PROHUuWLJG3L+/fvz/uu+8+eL0nbu6aNWswYMAAOTCpY/LkySgvL8emTZua/Vo1NTWyRnnopyA2WFsQK4IC5SGKS48Xvta1RSGsFA6eKIituxYFrqHgiWLXk9oBxdGEvXGBa+MjQgGsJ3QId4X/gn17d+OSYdnY0PNWnHvNvQxMCNEp7oQE5BXfgd4PfoP9udPgtAIXu9Ziftl/YVzNKrlurDEsiDWWRsCC2BjZsGGDHIxUV1cjMTERb7zxBvr27Sv7rrrqKuTn5yM7Oxvr16+XZ0S2bNmC119/XfYfOHCgQWAiqLsWvuZ47LHH8PDDD8c6VNIC+aYDuML0L1RX+DDlzAKsyZ+D226+W+thEUKiwO3JkVM9pd+9iw3Pz8MZnu040/wqxpZ+j+c8d2M/MrQeIiHqBScFBQVYt24dysrK8I9//AOzZ8/G6tWr5QBlzpw59ToxQ9KpUyeMHz8e27ZtQ/fu3Vs9SDEDc/fdJz40xcxJbm5uq/sjQC/TXlxjeg0lAR+GDuyL3Wc+iNkXX671sAghMeIZXIxxA8Zj3V9+gfBXz+GK1J+Qve9OrOw4C2+YJ2o9PELUSeuIupAePXpg6NCh8ozGoEGD8OSTTzapHT58uHzeunWrfBa1KCUlJQ00ddfN1akIHA5H/QqhukN7jLtap4ftIO4KvYAjFV6cMbgQlWf/EsUMTAgxLlY7Cq97Ar3vW41vAj1xdqcaXFr6DH5d9gA64ihX6xhMI2jvq3VaXXNSh9hATNSENIWYYRGIGRSBSAeJtNDBgwfrNStWrJCDjbrUkDFumUssvK1drSOeQSNW61idisNeu0JHrLoRwYTTDavNAcvx1TnyWXFtNpvk1TcWu7vJtlhRU3c0Za9ffQNTA21T/m7Ww7g5/HeUHivHmKIzsX/AXFx06VVa31RCSBvgzh+Ecb9Zi03dbkJqggsT3Vvwy7KFONv3EczHN3dUwtU6+tModS6u1ok+vVJcXCwXS1ZUVODll1/GqlWrsHz5cjl1I66nTJmCjIwMueZk3rx5GD16tLw3imDSpElyEDJz5kwsXrxYrjO5//77ceutt8qzI8bAuKt18rEfV/j/iQp/DUaNGIGtvW/C7KuvOf23jBCiHmYLim76HUq3XIGfnr8FXdwbMdf3DoYc3I5/dvkljuDEzDNX6+hPo9T5uFonOsSMh9iXRNSdiFqSr7/+Wg5MJk6cKKd7xBJhEYD07t0b99xzD6ZPn4633367/vUWiwVLly6Vz2IWZcaMGXJ/yn1RjIOxVut0Nh3CnNCrCPhr0K13IQ6deS+uYGBCSFw/ULDo0a9gOetO1EhO3NB5O3598E6cFfqqgY6rdfSnEXC1Tgw8++yzzfpEgaoojG0JsZpn2bJlsXxZcoqkowIP2/6BDWUVGDawD6pGz8e0aRdpPSxCyOnGbEHhNYuRM2wafnruJpydtQupB/8LgxPOwZ8SbkT7/vgjcV1zQvSNG9WYZ3sNm3YexKWjCrCz1w2sMSGkneHpdzaKHl+LjVkXolOiCVNCH2FhyZ0oMO3VemiENAmDkzjGihCusbyH0kOHcMmwfHyUdS3m3M59TAhpl9hcKLp7CfznPYUjplQUZ5XihgOPYpz/4yY3biNESxicxC0SLrGshtN7AEXdUvBNtzsw584FWg+KEKIxBZOuw+CHvsRa65nonR7GLPO7uPHo75GI9lt8SfQHg5M43eek2L4OXWp+RN8sK44MuR3nz2z4dGhCSPvFnZmHcb/+CJaz7oLJ4sAtnX/EfT/fgT6mXbrap4T7nLRfGJzE4T4nw2zbMcz/JfJTAMvwORh10Y0Nng5NCCEwmVB4zX9hxC/+hc1VGbg8vxx3/bwA59cslWdetd6nhPucBLjPCYmffU6yQ/swvvpdWM2Ao/+F6D7pRj7EjxDSLJ5BkzDuie/x0a8vRGHmN8gsfx5dKtbjKdtM1Oh4HxO97U/CfU7aFs6cxNE+J45gBS4JvYMDZWF0LhiChJE3oLCwUOsbRQjRO+50jPv1KuzvdwskWHCR+1vcfGAhMkMNHzeit31M9LY/Cfc5aTsYnMQNEm60LcOeg+UoHtQBJd0vR/HUqVoPihBiFMxmFF73W2Re+1ds9yWjuMNh3HD4DxhU863WIyPtEAYnccJlllUwH/oJk3rY8FXGhbjt7nu1HhIhxIDkjbwYRQs/w4/ojrOz/bg78CdcdPQvXG5MVIXBSRxQYN6DsdUr0cENlHS7FNfc/p9aD4kQYmDcnXrgoj98i8ouxchJMmOO5S3cVfoIXKjWemikncDgxOC4UIMZeAc7D/vR7+zzkDfqEhbAEkJOHZsT4x74F8zj7oM3aMEF7nW49+e7kY1SrUdG2gEMTgyNhOmWj1FeVobRA7pjR/Z5uOgiPjOHENJGmEwouPRBZF7zF7kOZXr2Idy++x70D27WemQkzmFwYmCGmbcgp2YbembYsaf/XFx6xUyth0QIidc6lIc/wWeVeZia68VtJQ9jRFXLD3olpLUwODEoqaZynFXzCVIdgLnwCpwx8TJutEYIOW24s3uh+ImvsdY6DAVpIVwfWIKxR99ioSw5LTA4MSQSzsUXSLSEkNV1ALLPmcM6E0LI6ceRhHGLPkRZryvlyzsyv8CMQ3+AHbV7LRHSVjA4MSADzTuQYz6ETmkJSJhwD4pGjdJ6SISQ9oLZgsLbnsOAqx7BoWob7snfgnk/34tUVGg9MhJHMDgxGMmmKow0b0IgCHQadz36nTlO6yERQtohedPuxZj5L2PjURcuytyDO3bNQ47pkNbDInECgxNDIWGKZQ0qq4HCgX3h7zqJ6RxCiGZ4hl2MQfPfxtqyNFySV4Zbdv0HCoI/aj0sEgcwODEQZ5h/hLOyBEM6m3Gw6yXcnp4QojmePiMx7pFP8FlFHs7L9eKOkocxoOJTrYdFDA6Dk1bjBExWwG4FLPaTD6sdsInDUdu22mG22mEyW+Vzg2ubEyaLFSZrM22LDW5LCGOCa9AxAajudi7Ov+p6rW8AIYTIuDt2R/F/f4G10gD0TavBf1Q9iZFVq2p/10XgdPuNqolFF68wOGnVLXOJ7RMBlwOw2QCbG7A6FYcdZqcbZoewW+W21eaAxWKD1Z1Ue1Zcm80mWBxuWOzuJttmix0XOL5BlbcaBd3y0GPCtUznEEL0hSsV4x79GAeyxiIvWcJ1wb9hXMUymMUfaE0g7OKp6qfLb1SNUudytd/tIRicxEwYgK/28NUANQEgUAUEqxWHH+HqKoSrvUAgILeDgRqEQgEEvRW1Z8V1OCwhVOOVj6ba2aG96Hj0O+QnAyVdLsTk4ila3wRCCDkZqwNF//kODnSZjhQ7cIP7Q1xY8kyTe6GEAzWAxVZ7boJT9RtVo9T5fF60VxictJogIAWBYBAIBYCQv+ERDMhBSu1R2w4H/ZDCQfnc4DpQDSkUlI+T2wEUh1bhUCVg73Im+o+9kJutEUL0i9mMonl/RcKYu+APm3FL+peYvX8xzPIfdg2RxO/CCJyq36iaWHTxCoMTnTPQsh2+IyUY3c2FivxzMW4clw4TQnSOyYTCaxejxyWLcLQauCbla1y792FYxR91hEQBgxMdY0MQhTXfIi8F2NthPGbddLvWQyKEkKgpmP4L5F71JLaXmXFj1ibctPcBONC+ZwRIdDA40TGDzdsQ8vuRnp6BflNuhMfj0XpIhBASE3kTb8Kgm1/ApiN2zOqwFbft/U8kyHV7hDQPg5NToqa23iQQaqLmpK4e5Xg74JdziKLGRM4lKq7DwZoTdSbH2+5QJXJ8PyLLDYT6XoRhI0Zq/WYJIaRVeEZcikF3LMHGo05cmrUTc3ctQAqqIAVaqCk5Rb9RNbHo4hUGJ6dECBBTlIEA4A+eOMQ3lbxqpwYIiyCl+rgtgGCNt7ZAVnEdliSE/D75qGsPCa2DVTxnK7UjCsZdyVkTQoih8Qw5D4PmvYa1hxMwO28/Lt6+CMlWP6RmVq7Idqut1X6japQ6n6/9zjAxODklxCY5DsBhA+yKQ2y+Ju934hA76dS27cJmg1Xe+8TW4NpsMsFid8mHaHscAfQ370KaA0gddjmGjxih9RslhJBTxjNgPMY9tBwr9yfiup7HcMG+p5Bha/qJxiaxF0goUHtuhd+oGqXO5RJ7arVPGJycEmJ3WHFY5DXpJw4ReIigxKpo22E6viOsOCuvzVbH8Z1ga9tjbZtQHQQGDh6AHiMmc+kwISRucHcdiov/+1N8cigZt/WvxLUHHkUGyprUyr8rI9CS36iaWHTxCoMTnZFhqkCmbzf6eoDK/MkoLi7WekiEENK2ZPbCNU99jq+PpeLGXmWYu+s/4GkmQCHtEwYnOuNM02aU+4GkbsNQNPlirYdDCCGnB09PXPD4h/jkUAouzzuGa7fdg0wc03pURCcwONERmaajSPcfRM90oDxvIoqKirQeEiGEnDbcuf0x7uEVWHkgBXN6lWHGlruRhaNaD4voAAYnOmK0ZQOCYcCZMxBnTblU6+EQQogqAcrFv/sM/9iZjFv6VuDqLZxBIQxOdDVr0qFmD/KTAKnHZBQWFmo9JEIIUYf0brjyD2vwj121Acrsbf+BNFRqPSqiIQxOdEKR+QeU1QCdehdi8BgWwRJC2hnpXXHl7z/H23uTMatbGa7a+ygS0H6fytveYXCiA1LgRWrlDgzMBCo7n8WH+xFC2ifp3XDuwuX4pCQRk7MO47JDzyAB1VqPimgAgxMdMMi8DTVis9nMXhgz7Sqth0MIIZrhzh+E4fP/iR8rE3FpbilmHfgNnOJRIaRdweBEY5zwo2PNTnRJAaTuEzBo0CCth0QIIZri6Xc2Jt3/GrZWJuLGrvtx3c5f8mnG7QwGJxrT37wTJqn2GTpnFF/N3WAJIeT4VveT/vNVrD3kxKzsnZi18wHYENR6WEQlGJxoiAUhdA38iGzxuJ2e52AQV+gQQkg9nsLJGH7Xi1h70Ibrs7dh1t5FMCGs9bCICjA40ZBuphK4TCEkpiai7/grOGtCCCGN8Aw9H8Nvfw6bDpsxI30Trtr/hHhur9bDIqcZBienRI385EgEQkDI3/AIBmuPunbADynohxQOymdxnV2zHU4zkDngHPQbwFoTQghpCs/wS9Dz6t9iX4WEq91rcP7+P8p2KdByHYoRNbHo4hWr1gMwNmKJjQg0rIAkNbyrVrH8zQSERZBSDZhM8hOMgzVeIBhAUuAo3KYaWExA2sBi5OXlaflGCCFE1+RNuhm+skMoXfYI5ro/xK59IXzqHgcp0PxKHtlntRlKo9T5fD60VxicnBLikdZ2wGGTA496LDbA6qxtS6Hatt0uf7NZHW74Q0F0dx4RoQvy+g1G78EjNHsHhBBiFAoufRDwHkXpqj/ieteHOFIqYVvqlc3qTTaHPLstnw2iUepcLhfaK0zrnBLW2qDEYqkNSOoPEYhYa4/6th0mcZitcFuBHNNR9M0AbD3O4Vb1hBASJQWznkDC4MtggYSZtpUYUPVFRL34vdsSetPEootXGJxoQE/zzzCbAEd6Fs6adrXWwyGEEONgMqHwthdh7joWfTLNuMa3BH1867QeFWljGJyojoTswE6k2IGU3qORk5ur9YAIIcRYmM0oXrQMpSmDMS43hKv2/QZ5wZ1aj4q0IQxOVKYTjiLJFJYzPTnDpsDj8Wg9JEIIMR5WOy76zUdY4+2K87vW4NJtDyEDZVqPirQRDE5UXkrcPbwVx3xAlz6FyO/eW+s3QAghxsWRiIt+uxrvHkjH3L5VuHTLArgaPShQb8uEuZQ4OhictMlS4gDgD544xDeVWD4crDmxlDjghylYDVv5fhSkAccSe6GoqEjrN0AIIcYmMRNXPLEa/9jmwj19D+OyrQ/AfHwXWb0tE+ZS4uhhcHJKiGpqR+1SYrvisIkVOk7A6gDM1vqlxN3sx9ApEfDbEzDt6pu0HjwhhMQHnh64cvEyvL3dipvzduDSPYvl+j69LRPmUuLoYXCi4lLirtgv3/CBo8Yju3NnrQdPCCHxQ95wjL7tKWw6FMbs5K8xef+fdblMmEuJo4PBiUokwQt3sBIdxEP+upyFgoICrYdECCFxhWfUTORMXYCj1RLm2N7DkMPvaT0k0koYnKhEnvmgfHZ26IYzx0zUejiEEBKXFFyxEAlDr5Lbs8L/Qn7Zd1oPibQCBieqrNYJIMm7D1luwN3tTPTq1UvrgRNCSNxSeMuz8OeNRYYzhNnSW+hQvSOinqt19AeDExVW69hDx5BmBfxhoO/Y6XC73VoPnBBC4hezGUX3vQl/Wm8UdQrgmmPPIBmVTUq5WkefMDhRYbVOrs0rflbQt28fdMrvofWgCSEk/rG5UPzoCmytzsA1Pctx+dYHYEXwJBlX6+gTBicqrNbpaDpSm9LpPoIpHUIIUYukjrhg4b+wco8D1+ftwYU7HpOXGDeGq3X0B4OT00wCquAyA247kDN0ElM6hBCiIu7uwzHm9qex9kAYt2Stw5g9tUuMib5hcHKaC2JTQ6UIhYHc3oPQtQe3qyeEELXxjLwSA6bfhx3HwrjVtRyFR1Y08LMgVn8wODmtBbHVSAx54bICUoeBKCws1HrAhBDSLsm7+FfIGHIhEm3AjIrn0T24TbazIFafMDg5jQWxiTYJbmttSmfQORdpPVhCCGm/mEwovP0l7HL2w8RcPy7atgjpKGdBrE5hcHIaC2KzbWVItAK5XbuiU25XrQdLCCHtG6sdxY+swLsHMnBT3ypc+O/7YUWIBbE6hMHJaaw5SfAfgsMMpHQp5CodQgjRA+50XPHbFfjHVjvu7LUP43b8L2tOdAiDk9NUc2IOlcMlATYLkFEwnKt0CCFEL2T1wQUL/h8+2BnEjIyNGFzxMWtOdAaDk9NUc5Jt88m1JpmdPOg3+EytB0oIIUSBe8glOOPiefCHTZhhX41eodoC2aZgzYn6MDg5TTUnWaajsJmBngPPRHZ2ttYDJYQQ0oi8yx5Hep+x6J8JzK58DumBkma1rDlRFwYnp6XmpApmfwDJoi62Uz8UFBRoPVBCCCFNPYNnwevYY8rD9K5eTN++EHYEmpSy5kTHwcnTTz+NgQMHIjk5WT6Kiorw7rvv1vurq6tx6623IiMjA4mJiZg+fTpKShpGort378bUqVPlGoysrCzMnz8fQbFhWRzVnDiCXtlrdjpRMLhI60ESQghpDkcSLvqv5Vi2y4XbepZg6k+PnCRhzYnOg5OcnBw8/vjjWLt2Lb755hucc845uOCCC7Bp0ybZP2/ePLz99tt47bXXsHr1auzbtw8XX3xx/etDoZAcmPj9fnz++ed48cUX8cILL+DBBx9EPNWcpFm9cFmAHn0Gom/fvloPkhBCSCQyuuOCX/4FK3eGcHv2Bozc81wDN2tOdB6cTJs2DVOmTEHPnj3lpbGPPPKIPEPyxRdfoKysDM8++yyeeOIJOWgZOnQonn/+eTkIEX7B+++/j82bN+Oll16Sd0stLi7GokWL8NRTT8kBizFnTiQgjOMPk6o9XKGg/KC/1G6DmNIhhBAD4B54HoZdtkDe4v4O57vocfiThoKTnxd4Mm2liUUXp7S65kTMgrzyyiuoqqqS0ztiNiUQCGDChAn1mt69eyMvLw9r1qyRr8V5wIAB6NChQ71m8uTJKC8vr599aYqamhpZozz0m9aphtMEVIWAPiOnaD1AQgghUZI3fSEyBkxCsj2Mayr/H1IDh2Q70zoGCE42bNggz5Y4HA7MnTsXb7zxhpy6OHDgAOx2O1JTUxvoRSAifAJxVgYmdf46X3M89thjSElJqT9yc3Oh17ROij0gNodFj26d0Sm3i9YDJIQQEssW93f+DbukHEzLrcIF234NM8JM6xghOBFpinXr1uHLL7/EzTffjNmzZ8upmtPJfffdJ6eN6o49e/ZAr2kdl+SFW2y81mWAXKNDCCHEQDiSULzoHSzb6cDcrnsxaccTtXamdfQdnIjZkR49esg1JWJGY9CgQXjyySfRsWNHuW7k2LFjDfRitY7wCcS58eqduus6TVOIWZq6FUJ1hz7TOgE4woDFBCR3GQSPx6P1AAkhhMRKZgEm3v47fHsghBuTPkfvkveY1jHaPifhcFiuCRHBis1mwwcffFDv27Jli7x0WNSkCMRZpIUOHjxYr1mxYoUcbBhzVUvDtI7VHkaiHUhPBvqcMVrrwRFCCGklnrOvRedRV8vtq8PLkOHdybSOilhjTa+IFTaiyLWiogIvv/wyVq1aheXLl8u1INdffz3uvvtupKenywHH7bffLgckI0aMkF8/adIkOQiZOXMmFi9eLNeZ3H///fLeKGJ2xHgo0joWCcmokq09exagE3eFJYQQQ1M490/46BfrMUjaiOm7/4lXMCRytoVpHW2CEzHjMWvWLOzfv18ORsSGbCIwmThxouz/3e9+B7PZLG++JmZTxEqcP/7xj/Wvt1gsWLp0qVyrIoKWhIQEuWZl4cKFMCZ1aR0rIEkwhQNwJAL+xM58CjEhhBgdqx3jHlyKv80dhKu6HoN33x/xr5Rfi7mNk6RM62gYnIh9TCLhdDrlPUvE0Rz5+flYtmwZ4gOR1rHXpnUsFqRKgMkCDBw5gU8hJoSQeCA5GxfcvwRv3n8e5nT6AYf3v4RPO808Sca0TtvCZ+u0UVrHhmrYzUBmMpCZ31vrgRFCCGkj3H0nYej0u3GsWsKM0NvoWv5t00KmddoMBidttFrH5S+DPwR0yul60l4vhBBCjE3BVY/B33kkeqWFMPPIU3AFGq5MZVqnbWFw0kardcSzdDLcQELH7hg3bpzWAyOEENKWmEwY98vX8W9vBqZkl2H6z//dYHqDaZ22hcFJm6R1JCSagZAEnDGWW9YTQkhc4k7H2HtfxtoDwJWpmzBq/5KGfqZ12gwGJ22S1imXv5EcorzYna71oAghhJwm3AVj0GPqncfrT95CXvk62c60TtvC4OSUMcNtCcj1JgU9cpGVlaX1gAghhJxGCq58FP5Ow1GQFsLlB/8AR6BMTsVI1VUtpnVa0ih1LqZ1SOsJw+wHOriBGqcHY8aM0XpAhBBCTidmM8bd/wbWlafi4txjmP7z7+TZDpMzocWZk5Y0Sp2PMyfkVApi09xAQOyAe96F3N+EEELaAwkenPuLJfi2BJiRtgHDy1ewILYNYXByyjUnNbCI50S5AUdmd60HRAghRCXcfcajR/Ht8AUknOtdhg7Vu1t+EQtio4LBySkHJwEEQkBCWjIcLs6aEEJIe6Lg6sdxMHkg+mUEMLXin7AFap+x1hQsiI0eBieniJg1sVuA7Lxu8pOZCSGEtCPMFhQvfAfbfKm4tEs5ig+90KyUBbHRw+DkFBHfYm4LkNShu/y0ZkIIIe2MxCxc9sAL2FYG3N5pAwYdfLtJGQtio4fBSRsEJxYLMGLC+VoPhRBCiEa4B5yHIcXXoToIXOn7O9IDB07ScOYkehictEFaJy1RPJI5SeuhEEII0ZCCa/+Ag45uKM714bxtj8MkngqrgDMn0cPg5FQxA52zO6Jr165aj4QQQoiWWB0ofuhNLNthw81d92DUjqcbuDlzEj0MTk4RmxlwpHdCTk6O1kMhhBCiNZkFmHjTInx7IIQbnB8h//AX9S7OnEQPg5NTQMS0LhvQpdcAeDwerYdDCCFEB3jG34H03mOQ6QrjsqN/lre3F3ATtuhhcHKK+8MmOYD8fmdqPRRCCCF6wWRC0X+8ii2+NHl7+wv2/v6Ej5uwRQWDk1PEYgeSM7O1HgYhhBA9kZCBc+99Ad+WSLgmfR2GHnyLm7DFAIOTU0DUYRfkd0RqaqrWQyGEEKIz3P3ORY8J18vb21/mfRUZwUMsiI0SBiengPj28tlTMWDAAK2HQgghRIcUzPofHHT2wOhO1Ri9508wsyA2KhicnAImAIOHDmMxLCGEkKax2uXlxav22nBt/gEU7Po7Z06igMHJKWCzAYkduL8JIYSQCHh6YuIND2HtgSCuSvkemcc2RJRLYmbFwpoT0sobZ7cB2V17az0UQgghOsczcR7yew9DzzQJ0yv/Dov8VPsWpubbMQxOTmHb+i4dXUhhMSwhhJCWMJtR/Ks3sK8mBTf2KMXYHX9oVmoSaZ8g9zkhrcAmjoRUdOrUSeuhEEIIMQIpnXHeHb+Vd4+d5foEXQ9/3qRMYkEsg5PWImbccrt0Q2FhodZDIYQQYhA8o69Der8J8DglXHL0WdgCFSdpTCyIZXDSWqwiCM7M1XoYhBBCjLZ77N1LsMWXKu8ee97ek9M7EmdOGJy0FqcF6DZgqNbDIIQQYjTc6Tj3nmfxXUkYM5K+Qb8jHzRwmzhzwuCktbicYvpE/IcQQgiJDffA85A9/FI4rcDFR19EYuBIvU/izAmDk9bWm6SnOdG9e3eth0IIIcSgFN78J/wY6IQpOVWYsuuJeruJMycMTlqDeFhkalo6n6lDCCGk9dgTMHnBX/H5zxJmezZjwL7XZbPETdgYnLR25sSWkIaioiKth0IIIcTAuHudjf7n3ohj1RJmh99EovfnWgc3YSOxIp6KMLJohNbDIIQQEgcUzPwt/MndMCStChce/D9uwsbgpHXYAfhtKVoPgxBCSDxgc2L47c/ih8PAlambUPDzUhbEaj0AI+KyAXndC7QeBiGEkDjB3XMUeky4QU7vXBheCXvZHs6ckNiw2gCHI/IjrwkhhJBY0zvHnPkYklGFib73OHNCYiM9NQFZWVlaD4MQQkg8YXNi/C9ewtajJszvuRvl65aivcLgpBV06tQRY8aM0XoYhBBC4oxw9lD4Bl6BypAVZ+//f4D3xOZs7QkGJ60gMc0Dt9ut9TAIIYTEIU8GL0ZG564wVR0E3r0X7REGJ60gNauz1kMghBASp9TAjvmBmyCZzMCG14Af2l96h8FJjIja6fwe/bQeBiGEkDjme6kHgiNur71YOq/dpXcYnMSITfzHYtV6GIQQQuKcwNn3Apm9AZHeee8XaE8wOIkRpwPo2LGj1sMghBAS71idwAV/BER6Z/2rwE8r0V5gcBIj+Z2SkJubq/UwCCGEtAdyhgLDb65tL70LqKlAe4DBSYyELC7068eaE0IIISpxzi+B1DygbA/wwSK0BxicxEh6ugd5eXlaD4MQQkh7wZ4AnPc/te2v/gTs+QrxDoOTGMlmSocQQoja9BgPDLoKgAT86zYgWIN4hsFJDJgA5HbrrfUwCCGEtEcmPwIkZAKlW4BP/hvxDIOTGLDIs2upWg+DEEJIe8SdDhQvrm1/8gRQshnxCoOTGBC7m6R7PFoPgxBCSHul30VAwRQgHADeuh0IhxCPMDiJAbcV8DA4IYQQohUmEzD1vwFHMvDzN8CXzyAeYXASAzkdE5Gfn6/1MAghhLRnkrOBiQ/Xtj9cBBzdiXiDwUkM+CUbevbsqfUwCCGEtHeGXAPkjwICXuDtuwBJQjzB4CQGMjLSmdYhhBCiPWYzMO33gMUBbP8I+P4VxBMMTmIgJT1d6yEQQgghtXh6AGMX1Lbf/2VcPbmYwUkMJCYzOCGEEKIjim4HMvsA3sPAigcRLzA4iYGefQdqPQRCCCHkBFY7MO341vbf/RXY+RniAQYnMRAIiz1iCSGEEB2RNwIYek1te+k8IOiH0WFwEgMdOnTQegiEEELIyUx46MTW9p8/CaPD4CRKbPJqnQyth0EIIYScjCsNmPxYbXv1b4DD22BkGJxEicsG5OXlaT0MQgghpGkGXAJ0GweEaoB37jb03icMTqLE5XZgwIABWg+DEEIIiby1vbz3ySpgw2swKgxOoqRrbiduwEYIIUTfZHQHxsyvbS//T8PufcLgJEr8Yd4qQgghBmDknYCnAKg6BKx8CEaEn7hR4nInaj0EQgghJLa9T759Edi1BkaDwUmUuBKTtB4CIYQQEh35I4HBMw279wmDkyhJTmO9CSGEEAMxcSHg9gCHfgDW/C+MBIOTKElN4x4nhBBCDIQ7HZj069r2x78Bju1BXAYnjz32GIYNG4akpCRkZWXhwgsvxJYtWxpoxo4dC5PJ1OCYO3duA83u3bsxdepUuN1uuZ/58+cjGAxCz3ClDiGEEMMx6AogrwgIeIHl9yEug5PVq1fj1ltvxRdffIEVK1YgEAhg0qRJqKqqaqC78cYbsX///vpj8eLF9b5QKCQHJn6/H59//jlefPFFvPDCC3jwQX0/TTEtLU3rIRBCCCGt2/vEZAF+eBv4aSWMgDUW8XvvvdfgWgQVYuZj7dq1GD16dL1dzIh07NixyT7ef/99bN68GStXrpSfVVNYWIhFixZhwYIFeOihh2C326FHOnXqpPUQCCGEkNjp0A8YPhf44ilg2X8At3wB2JyI25qTsrIy+Zyent7AvmTJEjkN0r9/f9x3333wer31vjVr1sg7rSofojd58mSUl5dj06ZNTX6dmpoa2a881MQByEEYIYQQYkjG/gJI7Agc3QF8pv8HA7Y6OAmHw7jrrrswatQoOQip46qrrsJLL72Ejz76SA5M/vrXv2LGjBn1/gMHDpz0dN+6a+FrrtYlJSWl/sjNzYWa2B1AcXGxql+TEEIIaTOcycDkR2rbnz4BHNmBuEnrKBG1Jxs3bsSnn37awD5nzpz6tpghEemQ8ePHY9u2bejevXurvpYIcu6+++76azFzomaAkpSUrNrXIoQQQk4L/afXbsq242PgvV8AV72KuJo5ue2227B06VJ5diQnJyeidvjw4fJ569at8lnUopSUlDTQ1F03V6ficDiQnJzc4FATt1PfuTlCCCEkquLYKb8FzDbgx/eAfy9DXAQnkiTJgckbb7yBDz/8EF27dm3xNevWrWtQUFpUVIQNGzbg4MGD9Rqx8kcEHH379oUesYq8DiGEEGJ0MguAoltr2+8uAPwnakING5yIVI6oJ3n55ZflvU5EjYg4fD6f7BepG7HyRqze2blzJ9566y3MmjVLXskzcOBAWSOWHosgZObMmfj++++xfPly3H///XLfYoZEjzj5XB1CCCHxwph7geQcoGw38Ml/w/DBydNPPy2v0BEbrYmZkLrj1Vdr81ZiGbBYIiwCkN69e+Oee+7B9OnT8fbbb9f3YbFY5JSQOItZFFEsKwKYhQsXQq8k8Lk6hBBC4gV7AnDuY7Xtz38PlNaWXRi2IFakdSIhilTFRm0tkZ+fj2XL9JvrakxiCjdgI4QQEkf0mQb0mABsXQm8Ox+Y8XptTYpO4LN1oiAjg8/VIYQQEkeYTEDxYsBiB7Z9CGz+F/QEg5MoEPU1hBBCSFyR0R0YdVdte/kvAX/DR9FoCYOTKOBD/wghhMQlZ98NpOQB5XuBT56AXmBwEgWNt+cnhBBC4gKb68TOsaI49sh26AEGJ1GQnZ2t9RAIIYSQ01cc220cEPID7/0n9ACDkygQy6IJIYSQuC6ONVuBH98Ffnxf6xExOImGwsJCrYdACCGEnD4yewHD59a231sABGtgyAf/xSNNbeNiAeD1B7UYDiGEkHaG1x/S7ouPWQBseK227mTNU7XFshrB4ESBLxBqcmrplYVXazIeQggh7Y8Hj38y21asBswqJziSOgGVJcDHvwUGXg6kdIYWMDhpATuA66zvaT0MQggh7Y2vNfzagarazdmGzNTkyzM4UZCRYMehimpkJjnl6wQAx169BoHUfK2HRgghpB1hNZu13U3ekQz0u1CzL8/gRIHJZIIn0dHiM4QIIYQQcvrgah1CCCGE6AoGJ4QQQgjRFQxOCCGEEKIrGJwQQgghRFcwOCGEEEKIrmBwQgghhBBdweCEEEIIIbqCwQkhhBBCdAWDE0IIIYToCgYnhBBCCNEVDE4IIYQQoisYnBBCCCFEVzA4IYQQQoiuMORTieueGlxeXq71UAghhBASJXWf23Wf43EVnFRUVMjn3NxcrYdCCCGEkFZ8jqekpDTrN0kthS86JBwOY9++fUhKSoLJZGrzqE4EPXv27EFycnKb9m1UeE+ahvelaXhfTob3pGl4X9rffZEkSQ5MsrOzYTab42vmRLyhnJyc0/o1xDdEvH1TnCq8J03D+9I0vC8nw3vSNLwv7eu+pESYMamDBbGEEEII0RUMTgghhBCiKxicNMLhcOBXv/qVfCa18J40De9L0/C+nAzvSdPwvjSNg/fFmAWxhBBCCIlfOHNCCCGEEF3B4IQQQgghuoLBCSGEEEJ0BYMTQgghhOgKBicKnnrqKXTp0gVOpxPDhw/HV199hXjhsccew7Bhw+RddbOysnDhhRdiy5YtDTTV1dW49dZbkZGRgcTEREyfPh0lJSUNNLt378bUqVPhdrvlfubPn49gMNhAs2rVKgwZMkSuNO/RowdeeOEFGIHHH39c3nH4rrvuQnu/Jz///DNmzJghv2+Xy4UBAwbgm2++qfeLOvoHH3wQnTp1kv0TJkzATz/91KCPI0eO4Oqrr5Y3kUpNTcX111+PysrKBpr169fj7LPPln/mxI6Yixcvhl4JhUJ44IEH0LVrV/k9d+/eHYsWLWrwjJD2cF8+/vhjTJs2Td7hU/y8vPnmmw38at6D1157Db1795Y14nt02bJl0Ns9CQQCWLBggTy+hIQEWTNr1ix5l/N4vienjFitQyTplVdekex2u/Tcc89JmzZtkm688UYpNTVVKikpkeKByZMnS88//7y0ceNGad26ddKUKVOkvLw8qbKysl4zd+5cKTc3V/rggw+kb775RhoxYoQ0cuTIen8wGJT69+8vTZgwQfruu++kZcuWSR6PR7rvvvvqNdu3b5fcbrd09913S5s3b5b+8Ic/SBaLRXrvvfckPfPVV19JXbp0kQYOHCjdeeed7fqeHDlyRMrPz5euueYa6csvv5THv3z5cmnr1q31mscff1xKSUmR3nzzTen777+Xzj//fKlr166Sz+er15x77rnSoEGDpC+++EL65JNPpB49ekhXXnllvb+srEzq0KGDdPXVV8vfl3/7298kl8slPfPMM5IeeeSRR6SMjAxp6dKl0o4dO6TXXntNSkxMlJ588sl2dV/E9/gvf/lL6fXXXxdRmfTGG2808Kt1Dz777DP552jx4sXyz9X9998v2Ww2acOGDZKe7smxY8fk3w+vvvqq9O9//1tas2aNdOaZZ0pDhw5t0Ee83ZNThcHJccQ3y6233lp/HQqFpOzsbOmxxx6T4pGDBw/KP0SrV6+u/wES38TiF24dP/zwg6wRP0x1P4Bms1k6cOBAvebpp5+WkpOTpZqaGvn63nvvlfr169fga11++eVycKRXKioqpJ49e0orVqyQxowZUx+ctNd7smDBAumss85q1h8Oh6WOHTtKv/nNb+pt4l45HA75F6ZA/GIU9+nrr7+u17z77ruSyWSSfv75Z/n6j3/8o5SWllZ/n+q+dkFBgaRHpk6dKl133XUNbBdffLH8YdFe70vjD2I178Fll10m/5soGT58uHTTTTdJWtJUwNbUH0NCt2vXrnZxT1oD0zoA/H4/1q5dK08/Kp/fI67XrFmDeKSsrEw+p6eny2fx/sX0o/IeiKnBvLy8+nsgzmKasEOHDvWayZMnyw+p2rRpU71G2UedRs/3UaRtRFqm8bjb6z156623cMYZZ+DSSy+V01SDBw/Gn//853r/jh07cODAgQbvSTwrQ6RClfdFTE2LfuoQevFz9eWXX9ZrRo8eDbvd3uC+iHTj0aNHoTdGjhyJDz74AD/++KN8/f333+PTTz9FcXFxu74vStS8B0b7uWr8+1ekf8R9EPCenAyDEwClpaVyPln5ASMQ1+IHLd4QT3UWdRWjRo1C//79ZZt4n+Kbvu6Hpal7IM5N3aM6XySN+LD2+XzQG6+88gq+/fZbuSanMe31nmzfvh1PP/00evbsieXLl+Pmm2/GHXfcgRdffLHB+4r08yLOIrBRYrVa5WA4lnunJ37xi1/giiuukANUm80mB23i50jUCbTn+6JEzXvQnEbv90jUsYkalCuvvLL+oX7t/Z7EzVOJyanPFGzcuFH+q689Ix5Hfuedd2LFihVy8Rg5EbyKv+AeffRR+Vp8CIvvl//7v//D7Nmz0V75+9//jiVLluDll19Gv379sG7dOjk4EQWO7fm+kOgRM7GXXXaZXDQs/gAgzcOZEwAejwcWi+WkVRjiumPHjognbrvtNixduhQfffQRcnJy6u3ifYr01rFjx5q9B+Lc1D2q80XSiL8QROW+nhBpm4MHD8qraMRfKeJYvXo1fv/738tt8RdHe7snArHKom/fvg1sffr0kVclKd9XpJ8XcRb3VolYwSRWJMRy7/SEWIVVN3siUnkzZ87EvHnz6mfd2ut9UaLmPWhOo9d7VBeY7Nq1S/6DqG7WpD3fk0gwOAHkqfuhQ4fK+WTlX4/iuqioCPGAiNRFYPLGG2/gww8/lJdDKhHvX0xVK++ByGWKD6S6eyDOGzZsaPBDVPdDVvdhJjTKPuo0eryP48ePl9+P+Au47hAzBmKavq7d3u6JQKT7Gi8zF3UW+fn5clt874hfdsr3JFJUIjeuvC8iqBMBYB3i+078XIn6gzqNWIIpfmkr70tBQQHS0tKgN7xer1wDoET8USPeU3u+L0rUvAdG+rmqC0zEkuqVK1fKS/SVtMd70iKtKqON06XEoqL8hRdekCun58yZIy8lVq7CMDI333yzvLxv1apV0v79++sPr9fbYNmsWF784Ycfystmi4qK5KPxstlJkybJy5HFUtjMzMwml83Onz9fXtny1FNP6XrZbGOUq3Xa6z0RKwmsVqu8dPann36SlixZIo//pZdearBcVPx8/Otf/5LWr18vXXDBBU0uFx08eLC8HPnTTz+VV0Qpl0aKVRxiaeTMmTPlpZHiZ1B8Hb0smW3M7Nmzpc6dO9cvJRbLRsWycbEaqz3dF7G6TSybF4f4CHniiSfkdt3KE7XugVg2K75Pf/vb38o/V7/61a80WzYb6Z74/X55OXVOTo78O0L5+1e58ibe7smpwuBEgdh/QnwQif1OxNJisd48XhA/ME0dYu+TOsQvj1tuuUVeria+6S+66CL5B0jJzp07peLiYnl9vfjFfM8990iBQKCB5qOPPpIKCwvl+9itW7cGX8NowUl7vSdvv/22HHSJgL13797Sn/70pwZ+sWT0gQcekH9ZCs348eOlLVu2NNAcPnxY/uUq9gIRS6uvvfZa+Ze4ErEPhli2LPoQH/zig02vlJeXy98b4neE0+mU/x3F3hbKD5j2cF/E93JTv0tE8Kb2Pfj73/8u9erVS/65Esv133nnHUlv90QEss39/hWvi9d7cqqYxH9anl8hhBBCCFEH1pwQQgghRFcwOCGEEEKIrmBwQgghhBBdweCEEEIIIbqCwQkhhBBCdAWDE0IIIYToCgYnhBBCCNEVDE4IIYQQoisYnBBCCCFEVzA4IYQQQoiuYHBCCCGEEF3B4IQQQggh0BP/HzwcoCuf2O8dAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "name = \"n-Propane\"\n",
    "\n",
    "j = json.load(open('PROPANE_exps.json'))\n",
    "sa = teqpflsh.SuperAncillary(json.dumps(j))\n",
    "import CoolProp.CoolProp as CP\n",
    "AS = CP.AbstractState('HEOS', 'n-Propane')\n",
    "def calc(T, rho, AS, key):\n",
    "    AS.specify_phase(CP.iphase_gas)\n",
    "    AS.update(CP.DmolarT_INPUTS, rho, T)\n",
    "    val = AS.keyed_output(key)\n",
    "    AS.unspecify_phase()\n",
    "    return val\n",
    "sa.add_variable(k='H', caller=functools.partial(calc, AS=AS, key=CP.iHmolar))\n",
    "sa.add_variable(k='S', caller=functools.partial(calc, AS=AS, key=CP.iSmolar))\n",
    "sa.add_variable(k='U', caller=functools.partial(calc, AS=AS, key=CP.iUmolar))\n",
    "Tcrit = sa.get_approx1d(k='D', q=1).xmax\n",
    "\n",
    "path = f'{teqp.get_datapath()}/dev/fluids/{name}.json'\n",
    "jresid = {\"kind\": \"multifluid\", \"model\": {\"components\": [name], \"root\": teqp.get_datapath()}}\n",
    "jidealgas = {\"kind\": \"IdealHelmholtz\", \"model\": [teqp.convert_CoolProp_idealgas(path, 0)]}\n",
    "rf = teqpflsh.RegionedFlasher(\n",
    "    ideal_gas=json.dumps(jidealgas), \n",
    "    resid=json.dumps(jresid), \n",
    "    mole_fractions=np.array([1.0])\n",
    ")\n",
    "\n",
    "# Now we make a region with rectangular shape in T, rho coordinates with only its 5 corners\n",
    "Tmin = 240 # K\n",
    "Tmax = 450 # K\n",
    "rhomin = 1e-6 # mol/m³\n",
    "rhomax = 12000 # mol/m³\n",
    "Tpoly = np.array([Tmin, Tmin, Tmax, Tmax, Tmin])\n",
    "rhopoly = np.array([rhomin, rhomax, rhomax, rhomin, rhomin])\n",
    "\n",
    "def makeVLE(Trange):\n",
    "    x, y = [], []\n",
    "    for T in Trange:\n",
    "        x.append(sa.get_yval(T=float(T), q=1.0, k='D')); y.append(T)\n",
    "    for T in reversed(Trange):\n",
    "        x.append(sa.get_yval(T=float(T), q=0.0, k='D')); y.append(T)\n",
    "    return np.array(x), np.array(y)\n",
    "\n",
    "ptr = teqpflsh.GeometryFactoryHolder()\n",
    "box = ptr.makeclosedpolygon(rhopoly, Tpoly)\n",
    "VLE = ptr.makeclosedpolygon(*makeVLE(np.linspace(Tmin, Tcrit, 1000)))\n",
    "rhoreg, Treg = box.difference(VLE).getXY()\n",
    "\n",
    "NT = 1000\n",
    "Nrho = 1000\n",
    "rf.add_region(T=Treg, rho=rhoreg, NT=NT, Nrho=Nrho)\n",
    "for reg in rf.get_regions_rw():\n",
    "    reg.add_pair(proppair=teqpflsh.PropertyPairs.PS, Nsplit=5)\n",
    "    reg.add_pair(proppair=teqpflsh.PropertyPairs.ST, Nsplit=5)\n",
    "\n",
    "# Build the helper class holding the models for ideal gas and residual Helmholtz energy\n",
    "helm = teqpflsh.teqpHelmholtzInterface(ideal_gas=json.dumps(jidealgas), residual=json.dumps(jresid))\n",
    "\n",
    "# Build the main flasher which includes all the parts\n",
    "mf = teqpflsh.MainFlasher(regions=rf, superancillary=sa, helm=helm)\n",
    "\n",
    "# Plot the region that is being mapped\n",
    "plt.plot(*box.getXY())\n",
    "plt.plot(*VLE.getXY())\n",
    "plt.fill(rhoreg, Treg)\n",
    "reg = mf.regioned_flasher.get_regions_ro()[0]\n",
    "propset = reg.propset_Trhogrid\n",
    "plt.plot(propset.rho, propset.T, 'k.', ms=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "77b589c2-ae22-49f3-8eac-36025f4e57d8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-06T11:32:40.508766Z",
     "iopub.status.busy": "2025-01-06T11:32:40.508680Z",
     "iopub.status.idle": "2025-01-06T11:32:43.653295Z",
     "shell.execute_reply": "2025-01-06T11:32:43.653011Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "907865\n",
      "3.4448237535615216 μs per flash call\n"
     ]
    }
   ],
   "source": [
    "reg = mf.regioned_flasher.get_regions_ro()[0]\n",
    "propset = reg.propset_Trhogrid\n",
    "\n",
    "val1 = propset.p\n",
    "val2 = propset.s + np.random.random(val1.shape)\n",
    "T = np.zeros_like(val1)\n",
    "rho = np.zeros_like(val1)\n",
    "q = np.zeros_like(val1)\n",
    "print(len(val1))\n",
    "\n",
    "tic = timeit.default_timer()\n",
    "mf.flash_many(teqpflsh.PropertyPairs.PS, val1, val2, T, rho, q)\n",
    "toc = timeit.default_timer()\n",
    "print((toc-tic)/len(val1)*1e6, 'μs per flash call')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "32439ed1-daf9-4b85-8e56-01acb2b53422",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-06T11:32:43.654543Z",
     "iopub.status.busy": "2025-01-06T11:32:43.654425Z",
     "iopub.status.idle": "2025-01-06T11:32:45.746390Z",
     "shell.execute_reply": "2025-01-06T11:32:45.746111Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.295959200958003 μs per flash call\n"
     ]
    }
   ],
   "source": [
    "reg = mf.regioned_flasher.get_regions_ro()[0]\n",
    "propset = reg.propset_Trhogrid\n",
    "\n",
    "val1 = propset.s + np.random.random(val1.shape)\n",
    "val2 = propset.T\n",
    "T = np.zeros_like(val1)\n",
    "rho = np.zeros_like(val1)\n",
    "q = np.zeros_like(val1)\n",
    "\n",
    "tic = timeit.default_timer()\n",
    "mf.flash_many(teqpflsh.PropertyPairs.ST, val1, val2, T, rho, q)\n",
    "toc = timeit.default_timer()\n",
    "print((toc-tic)/len(val1)*1e6, 'μs per flash call')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}