{ "cells": [ { "cell_type": "markdown", "id": "5d7dd3f0", "metadata": {}, "source": [ "# Analytical benchmarks\n", "\n", "Two analytical solutions from `pytes.validation` check the packed-bed model: Schumann's\n", "solution for a bed of a sensible material, and the melting of a sphere (the spherical\n", "Stefan problem). Both models are built from the packed bed's own submodels with `pytes.ComposedModel` (see [composing models](composing_models.ipynb)), swapping only what each benchmark needs." ] }, { "cell_type": "code", "execution_count": 1, "id": "3ee69c76", "metadata": { "execution": { "iopub.execute_input": "2026-11-01T11:17:13.301193Z", "iopub.status.busy": "2026-11-01T11:17:13.300720Z", "iopub.status.idle": "2026-11-01T11:17:14.333554Z", "shell.execute_reply": "2026-11-01T11:17:14.332218Z" } }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pybamm\n", "\n", "import pytes\n", "from pytes import validation\n", "from pytes.materials import SensibleMaterial\n", "from pytes.submodels import (\n", " BaseSubModel,\n", " DeviceEnergy,\n", " EnthalpyStorage,\n", " WallCoupling,\n", ")\n", "from pytes.submodels.storage import INITIAL_TEMPERATURE\n", "\n", "base = pytes.PackedBedModel()\n", "param = base.param\n", "geometry = base.submodels[\"storage\"].storage_geometry" ] }, { "cell_type": "markdown", "id": "820fed7d", "metadata": {}, "source": [ "## Schumann: a sensible packed bed\n", "\n", "The capsules are a sensible material, conductive enough to be lumped (Biot number 0.008),\n", "and the fluid enters at 70 °C into a bed at 32 °C. In Schumann's variables, `xi` is the\n", "number of transfer units from the inlet and `eta` the time since the first fluid reached\n", "the outlet, in units of the particles' time constant." ] }, { "cell_type": "code", "execution_count": 2, "id": "23fa3c8a", "metadata": { "execution": { "iopub.execute_input": "2026-11-01T11:17:14.335623Z", "iopub.status.busy": "2026-11-01T11:17:14.335239Z", "iopub.status.idle": "2026-11-01T11:17:15.005885Z", "shell.execute_reply": "2026-11-01T11:17:15.004946Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhgAAAFzCAYAAAB8X3AUAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAX91JREFUeJzt3XlcVNXfB/DPDPu+iuybAm64oOKCu6JgmgvuS6mZpmk/Syt5MlOjUMuy1BazKDW13MisXDMVTdTELUVkERGURZB9nbnPH+TUCKPDMHhZPu/Xw+9hzj33zmdGYr7ce+45EkEQBBARERFpkVTsAERERNT4sMAgIiIirWOBQURERFrHAoOIiIi0jgUGERERaR0LDCIiItI6FhhERESkdbpiBxCDXC5HWloazMzMIJFIxI5DRETUYAiCgPz8fDg6OkIqVX2eokkWGGlpaXBxcRE7BhERUYOVkpICZ2dnldubZIFhZmYGoPLNMTc3FzkNERFRw5GXlwcXFxfFZ6kqTbLAeHhZxNzcnAUGERGRBp40xICDPImIiEjrWGAQERGR1rHAICIiIq1rkmMw1CGXy1FWViZ2DHoK9PT0oKOjI3YMIqJGhQVGNcrKypCUlAS5XC52FHpKLC0tYW9vz3lRiIi0hAXGIwRBwN27d6GjowMXF5fHTiJCDZ8gCCgqKkJGRgYAwMHBQeRERESNQ70oMO7evYuEhAT4+vrCwsJCrX1SUlKQnp4Ob29vrd5qWlFRgaKiIjg6OsLY2Fhrx6X6y8jICACQkZEBOzs7Xi4hItICUf88P3/+PMaMGYMOHTqgd+/eiImJeeI+JSUlCAkJgY+PD6ZOnQp7e3usW7dOa5lkMhkAQF9fX2vHpPrvYTFZXl4uchIiosZB1ALj6tWrGD9+PKKjo9XeZ/ny5Th79iwSEhJw/fp1bNu2Da+88kqNjqEOXotvWvjvTUSkXaJeIpk2bRoA4M6dO2rvExERgTlz5iiulY8cORLt2rVDREQEunXrVhcxiYiIqIbqxRgMdaWlpSE9PR2dO3dWavf393/s5ZXS0lKUlpYqHufl5dVZRqpq5MiRcHZ2xvr166t9TNSYCIIAmawMZWUFKCzOR2rabWRkZSC/IBd5hbkoLCxAUUkhyspLUV5ejtbtPCBIZJAJMsTfuI3M9PuoqKiATFaBcpkMsooKyOQyyGQV8AtoAameFHKhAvFXU5GWnA0BqPySV34nCAIAwK+fJwyN9QEISLyejjvxWf/kq8wIQUBlTwEd+7rDxNwAkAhIvp6FW9czFa9FgAA83AdAx35uMLMxBADcictG4uVMle+Fb29nWDU3AQCkJeTgZky6yr5tezrBxtEUAJCenIsb5+6p7NvK3wF2rpVj77JS83HtzzSVfb397GHvUTm2L/teAa5Epars26KDHZy9rAAAuZlFuHg8RWVf93Y2cG1lAwDIzylBzNFklX1dW9nAvZ0tAKAorxTnDt5S0VOAs5cVWnS0AwCUFJUj+pdElcd18LSEd+fmAIDyUhlO74uv0sejnS2GDxyN54b+n8rj1JUGVWBkZ2cDAGxsbJTabWxsFNuqEx4ejuXLl9dpNrGFhYXh7bffVmqzsLDAgwcPlNo2b96MsLAwJCcnw8vLC+Hh4Rg+fPhTTEpUP8kqKnDvbgJuJl1BXFIsktNu4W56GrLuZyGvMB9dhrmiVF6CUqEMZ/Ym4158HspKZSgvlaOiVA5ZqRzyCgFChYC237RVXHa7vf428s6r/qOm9eetoWNUObD4zp47eHDygcq+l9umQ89SDwCQdi4N2UdV/9671aUIBsYGAIB7V+4h69cslX3TOlfA0KqyaMiIz0DG/gyVfTM7A8bNK8csZd3Owr0DqguBjPZymDpVFg3ZqdlIO6C6EEj3kcH8n6Ih514OUg+oLgTuupfD0tMSAJCbmYuUA6oLgVSHUlh7WwMA8nPykXxQdSGQalsCmzaFAIDC/EIkHUxS2TfFrAjNfEsAAMVFxUhQWTQAyfqFiOtU+UduaWkpbh5S3feWpAAJXSsAAGUVZYg7pDqvVb883OpeOW6worQCsdX0TTErhse90yqPUZcaVIGhp1f5H1ZJSYlSe3Fx8WMHZYaGhuK1115TPH64Elxj061bN5w5c0bl9gMHDuCFF17A119/jeHDhyMiIgKjR4/GmTNnqpwVImpMZBUV+PvvKPwe9Quuxl9DSuodPCjMQ7uRdshBMXKkMpz/MBkFsYXVH0AKZD2rA4nOP0VDRiHybqroC0CoECDRq+wr1ZdAqi+BRFcCqa4EUl1p5f+XAhKpBC7lgL4OIIUEclsD6LQwhkQqgUQqgVRS2Uf6z+MOFUYwKtGDRCJFgosVUv0lkKDyeSSSyv+R/POgm8QSxqX6kECCBJcK3AnQ+acTUNlVAlT+H7ro28K0zBASAMkuUtzqr6cokJT7StDRxBEW5UaQQII7zoZICDRQ+T60N3eCVfk/ZzAcjHEzUE9l37ZWTrAtryxG7tmZIG6w6o+n1s0cYFdeWYxk2pjh2mDVwwm97ZvDodwSAJBtWYgrg1V2RQtHOziXVxYjD8yKYD1YUNnXw9UWruX/nMEwKoHlYJnKvq6e1vAobwYAKDIohfngysHk1Y38cmppjZbl/5zB0C2H+RDVA88dPS3gVW4PACiHDKZDSqv08XCzha99L5XHqEsNqsB4OC9FaqpyZZuamgpXV1eV+xkYGMDAQPV/BE3Fhx9+iOHDh+O5554DALz22mvYvn071q5diy1btlS7jyAIWLt2LT7//HOkpKTA19cXH330EXr16qXY/vHHH+Ozzz7DnTt34Onpiddeew0zZ85UO9fBgwfxf//3f4iNjYWjoyNefPFFLFy4kLeLUo3JZTKkpZ1H9PXfcCPzb+z+/gwSrmQhP6MU8rJHPiykQOEYE0ikEgBSSE0qf96kBlLomerAwEQHRiZ6MDbRg4mpIQaUOsNI1xgGusa4G+yB8j6AuZkFLMysYWVpC2srO5iZWsDMxAIuzm4wNjKDqbEFjCcawUBPV72BxC/W4MXOqkFf0kxoDfq+UYO+rz25i8IrNej7cg36PgX1vsCIj49Hfn4+OnXqBGNjY/Ts2RP79u3DlClTAACFhYU4cuQIli1bVifPLwgCistVV6Z1yUhPp0Z3N1y6dAlmZmYwNDSEv78/Vq5cCV9fXwCVr+PMmTN4//33lfYZOHAgdu7cqfKYS5cuxRdffIFNmzahf//+uHHjBr755htFgfHuu+/ihx9+wLZt29CuXTucO3cOY8eOhampKSZMmPDEzAUFBRg1ahRWr16NGTNmICMjAxs3bsSVK1fQsWNHtV87NU0xMUfx3Y+f4UxMDJISM5GTUQLvj33+KRqAlJx85N7554ynDmBkow8zG31YWBrB0tocA8o7wMHaA3aWnjBeYwt3p1Zws3eDubE+pFLeWURUG6IWGBkZGYiLi0NmZuUgoStXrkBXVxeurq6KMxIrV67EmTNncPXqVQCVYw0CAwMRGhqKHj16YN26dbCzs8OsWXVTzheXy9Bm6cE6OfaTXFsxBMb66v0T2dvb45tvvkFQUBAePHiAt956C7169cLVq1fh4uKC/Px8FBYWolmzZkr72dnZ4d696q+hFhUVYc2aNfj4448xYsQIAEDXrl3RtWtXAJWXqlavXo3IyEj4+/sDAPr27Yt58+bh66+/VqvAyMjIQHFxMZ599lkYGxvD3d29ShFE9FBBQTo2fPMOfoj8FTf+TkdRRtX1girulqKlrSGayczRrk9HGPa2RYf2vdCvxwi0dnOFgS7PjBE9DaIWGH/99Rfee+89AEBAQAB++OEH/PDDD5gxYwZmzJgBAPDy8lJadKxv3744duwYNmzYgLNnz8LX1xdbtmyBqampKK+hvvjvJQkrKyt8++238PDwwFdffYUVK1ao3E8QBJVnSWJjY1FcXIw+ffpUu/369esoLCxEUFCQ4lgPvzw9PdXK7eHhgeDgYPTo0QMTJkxA//79MXDgQMXsmkS//7EDF+/9guiCq7ihW460IxnIPPbPnQtSwNLNCI5ulvDwbIEufgMwPGgG2ru7QE+H0/wTiUnUAiM4OBjBwcGP7fPmm29WaQsICEBAQEBdxVJipKeDayuGPJXnqu65NaWvrw9vb2/cvHkTAGBmZgYTExPF2aKHMjMz0bx582qP8fC2NFUFyMPF4C5fvow2bdpolFMikeCXX37BiRMncPjwYbz11luYPXs2Dh8+rPExqeG7du1PvL36VRw5cgl5qSVw/Z8rzDuZA5DA09cUzfJ10aFDF4SMWoCg7n1gYlDvr/YSNTn8r/IJJBKJ2pcp6pOysjLExcXBz88PQOXr6N69O44dO4ZXXvl31NDvv/+Onj17VnuM1q1bw8jICCdOnECrVq1Ubv/tt99qVQxIJBL07dsXffv2RVhYGLp27YpvvvkGH374ocbHpIZHVlGBT798E2s3RuD2lRz8M0kDJDoSWKaUI7CVLzq0mIABIwLhZMkzXET1Hc8hNhJjx47FiRMnkJ+fj6SkJEybNg0PHjxQGpuyaNEi7N+/H5s3b0ZOTg4+/vhjxMTEYMGCBdUe09jYGK+++iqWLFmCn3/+Gfn5+Th//jzmzJmj2P7GG29gxYoV2L59O/Ly8pCcnIzPP/8cYWFhauU+c+YMZs+ejUuXLqGkpAQXL17EnTt30KJFi1q/J9QwlJcX4/Ptc2HhbILX5n2E25criwvrFsYYOKkzNv60Dye2JuPDl7ZhauCzLC6IGoiG96c5VWvBggVYsWIFoqOjYWxsDH9/f0RHR8PHx0fRJygoCF9//TXeffddvPjii/Dy8sLu3bsfOwdGWFgYLC0tsWDBAqSlpaFDhw5KZxaWLVsGOzs7hIWFYdq0aXB0dMSwYcOqTPqlSteuXXH58mU8//zziIuLQ/PmzTFr1izMnj1b8zeDGoTS0gL8+Ptb2Jx2FHd1AJmhFDpGUvh0a45nR76EuVMXwcWaKxoTNVQS4eGF9iYkLy8PFhYWyM3NrbLUe0lJCZKSkuDh4QFDQ0OREtLTxn/3p6e8vAyvvTUO3/1wEM6hHtAx0oFVhRwuyfYY3OstTBjQl3d6ENVjj/sM/S+ewSCipyZiy7tYEPo+8lIr56YoOXIfwX37YGjACgRPb8W5J4gaERYYRFTnbiVfx6hpg3Hxj8qVk3WMpega6IVXF3yPMT39WFgQNUIc5ElEdWrVxy/D29dXUVx4BTRD+Dc7cGzn3xjXqzOLC6JGimcwiKhOlJbmY/XeKfhgyzGU58tg4miAYZPG4IPQjRy8SdQEsMAgIq1LSTuHVw7MRLyOHE7TnWBoV4DXQ/diYp/ONVpfh4gaLl4iISKten/NHHSbFoSbUhmsKuQIxCDs23YDk/p2YXFB1ITwDAYRaYVcLsfIqf74edtfAAAXTxMMH7wWr04bCV2uC0LU5LDAIKJaKysrRY8gb1w4dhsA0Lpvc7w57ReM7u4rcjIiEgsLDCKqlQcPMuE3oBWSYrIBCdA7xBeffXwc7ZytxI5GRCLieUuqF1asWIHPPvtM5WOqn7Jz0tG6ewskxWRDoidB0PT+2PFVNIsLIuIZjMairKwM33//PY4fP47y8nJ07doVs2fPhpGR8sJQN2/exLp165CcnAwvLy+8+uqrcHJyqnEfbbtw4QKcnZ1VPqb6p6ysEJPC+uFeXD6khlKMnD0Wm97bDCsTfbGjEVE9wDMYjUTPnj0RFRWFfv36ISgoCBEREQgICEBxcbGiz82bN+Hv74/79+9jwoQJiI2Nhb+/PzIyMmrUh6isrBDztwcitb0uPF50QsjcydiyeiuLCyJS4BmMRuLAgQOwtbVVPB40aBCcnJzw66+/IiQkBADw7rvvwsvLC1u3boVEIsGYMWPQsmVLrFmzBqtWrVK7T3UuX76MTZs2ISUlBb6+vli4cCEsLCwU2y9evIhNmzbhzp078PT0xJw5c+Dl5aX267t79y7Wr1+P2NhYODo6Yvr06fDz86vp20RaUFZWilciBuO0YT4M5AL6t5+E8Jnvw1ifv06I6F88g/EkggCUFYrzVYOFbv9bXACAlZUVdHV1UVBQoGg7ePAgRowYoZiLQE9PD8OHD8eBAwdq1OdR+/btg7+/P8rLyzFx4kTo6elh0qRJiu2//vorBg4cCDs7Ozz33HMwMDCAn58fLl26pNZrq6ioQO/evXHjxg1MnjwZbdu2xbx58xAbG6vW/qQ9giCg97DW+HbFX5BnlqJ7wRC8P/N9mBiwuCAiZfyt8CTlRcD7juI89/+lAfomGu366aefQiqVYuDAgQCAoqIiZGRkwMXFRamfs7Mzbt26pXafR8nlcsyZMwevvPIKVq9erWjPy8sDUPmBNGfOHHzwwQeYMWMGAGD06NHIzc1FWFgYdu7c+cTXcuvWLSQkJODPP/9Es2bNAAAvvvgiSkpKnvxGkFaNnd4DZw8nARLA7boT3lu3EqYsLoioGvzN0AgdOHAAb731FjZs2KAYKFlWVgYAVQZ9GhsbK7ap0+dRsbGxSEtLw4QJE5Tazc3NAQA3btzA7du3sXnzZuzZsweCIEAQBCQnJ0NQ8wyNi4sLnJycMGvWLMydOxc9e/aEiYkJTEw0K75IM++8PwO7v4sGAPQc2RobP/wFFsZ6IqciovqKBcaT6BlXnkkQ67lr6Pfff8fo0aMRFhaGWbNmKdpNTU2hq6uL7Oxspf7379+HlZWV2n0e9fASjKrtubm5AIDnnnsOjo7KZ4IeLWRUMTAwQHR0NNavX4/Q0FBcu3YNI0aMwJdffqkoZKhuHTqyFWErvgMAtOlnjy8+OQVnKy5YRkSqscB4EolE48sUT9uxY8cwfPhwvP3223jzzTeVtunq6qJt27aIiYlRao+JiUGHDh3U7vMod3d3AMC1a9fg4eGhcruZmRmCgoI0eVkAACcnJ4SHhyM8PBx37txBr169sHbtWixdulTjY5J60tISMXryC5CXytHMxwQffPgHfF04zwURPR4HeTYSx48fx7Bhw7BkyRKEhoZW2+f555/Hjz/+iKSkJACVc00cPHgQzz//fI36/JednR2GDx+OZcuW4f79+wAAmUyG77//HgDQvHlzjBgxAkuXLsWdO3cU+928eRN79uxR67XduHEDe/fuVTy2t7eHubk5ZDKZWvuT5gS5HIPG90RhRhkMbPQw763PMbSzj9ixiKgBYIHRSIwaNQpSqRTHjx9HUFCQ4mvr1q2KPvPnz0dwcDB8fX3RrVs39OrVC7Nnz1YaP6FOn0d98803MDc3h6enJwICAuDm5oaUlBTF9oiICHh5ecHHxwc9e/ZE69atMXLkSFhbW6v12mxsbLBlyxY4OjqiX79+8PDwgJGREebPn6/BO0U1sf3IG5CHWMK0tQlCZk3HW5OniB2JiBoIiaDuSLtGJC8vDxYWFsjNza1yDb+kpARJSUnw8PCAoaGhSAlr7vDhw9X+Rd+yZUu0bNlSqS0hIQG3b99Gy5Ytq9wxUpM+j7p58ybS0tLQtm3bKrfNAsDt27eRkJAAR0dHeHl5QSr9t76NiYmBoaEhWrduXe1jALh37x7i4uLQvHlz+Pho96/ohvrvXpcSk09g3O9zUSqVoGeOG1bN2gNLY06kRdTUPe4z9L9YYDSSAoNqh//uyvJy72Pocn/kdjBG6yIpXh9yCF09m4sdi4jqAXULDF4iIaIqhoz3x6mPE/Fg9z30dHqXxQUR1RgLDCJSsunb5ThzMBEA0Nm5O+Y9M0zkRETUELHAICKFvNxsLHjzfQCAd0AzfPLuDujq8NcEEdUcf3MQkULIjD4ozCiDvqUuXl74NdxsGsYcMERU/7DAICIAwM+/fIUjkX8DAALH9MXLI3hphIg0xwKDiCCrqMD0+QsAOeDSyRIfLN8FHalE7FhE1ICxwCAi/BIVDqMh1jB0MMCMl1ejtaOl2JGIqIHjWiRETVxx8QN8kvADLLtbItinNRY/N0PsSETUCPAMBtULP//8M/744w+Vj6nurNs3Gxm6EjQvl2NE749hqKcjdiQiagR4BqMRunv3Ln744Qd4eXnhmWeeUdpWUVGBQ4cOITk5GV5eXhg4cCAkEkmN+2jb119/DWdnZ/Tr16/ax1Q3Tp6KxJK5e2E3yg4DOz6DoX5txI5ERI0EC4xGRi6XY8qUKYiOjkZQUJBSgVFQUIBBgwYhMzMTvXr1wnvvvQdfX1/s27cPenp6avehxmP6/2aiPLscFWdyMWv5qjovJImo6WCB0ci8//77MDc3r/Yv/5UrVyI1NRWXL1+GlZUVUlJS0LZtW2zcuBEvv/yy2n2qU1ZWhiNHjiAlJQW+vr7o2bOn0vbS0lIcPnwYd+7cgaenJwYMGABd3Zr9+J08eRKxsbFwdHTEgAEDYGRkVKP9SdmPu9ci4a/7gAQYO2Uu2jhZih2JiBoRjsFoRE6dOoWNGzfiq6++qnb7jz/+iPHjx8PKygoA4OLigmHDhuHHH3+sUZ9HJSQkoG3btnjttdcQHR2N119/HbNnz1ZsT0xMRLt27fD+++/jwoULCA0NRZcuXZCdna32awsJCcGUKVNw7tw5fPXVV+jRowfS0tLU3p+qevXtdwAAnv42eGfeCpHTEFFjwzMYTyAIAoorikV5biNdI7VPWefk5GDy5MnYtGlTtUull5WVIT4+Hq1atVJqb9WqFQ4fPqx2n+pMmTIFLVq0wL59+6CvX7mc9+nTpxXbp02bhhEjRuDDDz8EUPmeBgcHY9myZfj000+f+Npu376NPXv2IDExER4eHgAqixaeztfcZ1+9hbTreZDoSDDx+cWwNTUQOxIRNTIsMJ6guKIY3bZ1E+W5oydFw1jPWK2+L7zwAkaOHInBgwdXu72wsBCCIMDS0lKp3crKCvn5+Wr3edTt27dx5swZHDt2TFFcAFBcIklNTcXJkyfRo0cPrF+/HoIgQBAEGBkZ4dSpU2q9NjMzMxgYGODnn3/G7NmzYWBgAE9PT7X2parkcjnefn8tAKBVr+Z4Y/r/xA1ERI2SxpdI/vrrL8yfP19pEOF3332HgoICrQQj9R09ehS//PILHBwcsHbtWqxduxZJSUmIj4/H2rVrkZeXB2PjykLl0ULhv9vU6fOou3fvAgDc3d2r3X7nzh0AQFZWFuLj45GQkIDExES4ublhxIgRar0+Kysr/Pjjj4iIiIC1tTUGDRqE7777DoIgqLU/Kdvw9SJk3yqC1ECKabNWwdyQg3eJSPs0OoPxyy+/YOzYsRg5ciR+/fVXRXtycjLWrl2LJUuWaC2g2Ix0jRA9KVq051aHnZ0d5syZo/iwB4Di4mKUlZXh1q1bqKiogLm5OVxdXZGQkKC0b0JCAry8vAAABgYGT+zzKBsbGwDAvXv3qi0yHm6fOnVqrW45ffbZZ/Hss88iMzMT+/fvx7x58yCTyTBjBieFqglBEBAl/QOeSz3heksfr4ydLHYkImqsBA34+fkJO3fuFITKPyEV7bGxsYK7u7smh3yqcnNzBQBCbm5ulW3FxcXCtWvXhOLiYhGSac8zzzwjhISEKLXNnz9f8Pb2FkpKSgRBEIScnBzB1tZWeP/992vU51GtWrUSpk6dqtQWHx+v+L5NmzbC8OHDBZlMpmirqKgQrl27png8YsQI4eWXX672cXp6upCenq50/L59+woLFy58/JtQA43l3/1JTl+IENp9207w+6at8M2h38SOQ0QN0OM+Q/9LozMY169fx9ChQwFAaaCdk5MTUlNTa3y8lJQUpKenw9vbG+bm5k/sLwgCkpKSkJOTAxcXF9jZ2dX4OZuiJUuWYP/+/ejfvz8CAwMRGRkJJycnzJ8/v0Z9HvXtt98iKCgIgwYNQt++fXHjxg2UlJRg165dAIAtW7YgKCgIAQEBCA4ORk5ODo4cOYI5c+agdevWT8ydm5uL4OBg9O7dGz4+PoiNjcW5c+ewZs2a2r8pTcxHx9cCVkDHAktM7BcodhwiasQ0GoNhaWmJlJQUAMoFxunTp+Hs7Kz2cUpKShASEgIfHx9MnToV9vb2WLdu3WP3uXz5Mtq0aYOePXti1qxZ8PDwQEhICIqLxbnTo74aMWIEhg1TXm7bzs4OMTExmDx5MkpKSjB//nz8+eefMDU1rVGfR3Xr1g03btzAiBEjUFRUhGHDhuGHH35QbPfz80NcXBxmzJiBwsJCuLu7Y/fu3Zg7d66iz/Dhw9G/f/9qH3t5eeHixYvo27cv8vLy0KlTJ8TFxaFz5861fp+akk3fLcfu1y8jfftd+Hu+yinBiahuaXJ6ZOHChcKgQYOElJQUQSqVCiUlJcL+/fsFFxcX4e2331b7OIsXLxacnZ2FtLQ0QRAEYe/evQIA4cyZMyr36dmzpzBkyBChvLxcEARBSE5OFiwsLIQPPvhA7edtCpdIqGaawr+7YxtzAYDg3ctOKCgpFzsOETVQ6l4i0egMRlhYGKytreHi4gK5XA5TU1MMGzYMPXv2rNEAz4iICMycORMODg4AgJEjR6Jdu3aIiIhQuU9mZia6du2qmAXS1dUVzs7OyMzM1OSlEDUJu3/agLRreYAUmDj5DZgY8A51IqpbGv2WMTQ0xA8//IB3330Xf/31F+RyOfz8/NS6nv5QWloa0tPTq5zm9vf3R0xMjMr9VqxYgUWLFsHd3R1ubm44dOgQioqKlE63P6q0tBSlpaWKx3l5eWrnJGoMlqx6FwDg7meFV59/ReQ0RNQUaFRgODg44O7du/D29oa3t7dGT/xwmuiHtzE+ZGNj89gppAcMGAB/f3+EhobC2dkZiYmJWLx4MVxdXVXuEx4ejuXLl2uUk6ihu3I1CrFn0wEAw0ZOhoUR570gorqn0SWSsrIy5OTk1OqJH67MWVJSotReXFysNCPkfwn/TDENVE7gdOHCBVy5cgWffPIJ3nvvPZXPFRoaitzcXMXXwwGqRE3B/96eBcgAm5bGeHPOSrHjEFEToVGBMWnSJHzwwQe1mknRxcUFUqm0ym2tqampKs9GpKWl4cKFC3jxxRcVRYiLiwueffZZ7Nu3T+VzGRgYwNzcXOmLqCkoLMzBqWM3AQB9AgPgbG0iciIiaio0ukQSFxeHQ4cOYceOHWjVqlWVMw6RkZFPPIaxsTF69uyJffv2YcqUKQAq18I4cuQIli1bpugXHx+P/Px8dOrUCdbW1pBKpYrppx9KSUlBs2bNNHkpRI1aZFQ4PJa3QPnJHLy5YJPYcYioCdGowPDy8lI5dXRNhIWFITAwEKGhoejRowfWrVsHOzs7zJo1S9Fn5cqVOHPmDK5evQojIyO8+OKL+L//+z/IZDJ4enri0KFDOHDgAH755Zda5yFqTAS5HD/cOQA9Sz0M7t8d3bxVj1MiItI2jQqM9evXa+XJ+/bti2PHjmHDhg04e/YsfH19sWXLFqVJnby8vFBWVqZ4vGHDBnTr1g0HDx7E7t274ebmhtOnT6N79+5ayUTUWJy68B2SdAUYyAUM7PKW2HGIqImRCLUZSNFA5eXlwcLCArm5uVXGY5SUlCApKQkeHh4wNDQUKaE4UlJSYGho2CQvNzXGf3cv/2a4W1iMAc94Yk/4RejqaLx4MhGRwuM+Q/9LozMYvXr1euz2qKgoTQ5LWlBSUoKMjAw4OTlBR6dmU0FPnz4dXbp0wcqVvNOgoYu9cRbxF7IAGdB6RjCLCyJ66jQqMAYNGqT0WC6X4+bNm9i5c6fS+Al6epKTkzFnzhz8/vvvaN68OXJycjB8+HCsWLECLVq0EDsePWWvr5gDyAArTyMsnLlM7DhE1ARpVGD89y6P/woODsZvv/1WmzykoQkTJsDU1BR3796FlZUVysvLsWfPHvz5559VCozi4mJkZGQobhWuTk5ODuRyeZWJ0G7dugVzc3NYW1sr2tLS0iCVSmFvbw9A+VJLTk4OysrK0Lx5c0X/O3fuVHu7cHp6umJ+FRsbm2ov1Tx67OoyqtOnMSsvL8PRQ1cAAD36dIKdhZHIiYioSdLmAih5eXlCs2bNtHnIOtHYFjsrLi4WJBKJsGXLlsf2KyoqEmbNmiUYGBgIDg4Ogq2trbB582bF9oEDBwqjRo0SunbtKjg5OQmGhoZCYGCgUFhYqOjToUOHKgvLjR8/XnjhhReUjjN06FChbdu2gpOTk6Cnpyc888wzwpkzZ4TWrVsr2ubPn690nHfeeUfw8fERfHx8BCsrK8HT01M4fvy4Uh91MqrTp7r3sKH9u6uyau08AYCgYyIVDpz9S+w4RNTI1OliZ6pcuXJFafn2xqSwsFDl16OzkT6u76PLyqvqVxOGhoZwdHTE1q1bq0xc9l/z5s3D0aNHcenSJaSlpeHmzZu4e/euUp8DBw7gk08+wZ07d3Dr1i38/fff+OKLL2qUBwD++OMPRERE4M6dO7h8+TIOHTqE4OBgbNu2DXfu3EF0dDQ2bNiA6OhoxT7Lli1DbGwsYmNjcf/+fcyZMwcTJ06s8p6pk1Fbr6Mh+uzb7wEArf0dENi5k8hpiKip0ugSycyZM6u05eTk4ODBg3jppZdqHao++u+ts48aOnSo0jwcdnZ2KCoqqrZv37598ccffygeu7u7Iysrq0o/oYY392zfvh3PP/88nJ2d4e3tjYCAAIwZMwZDhw4FADx48ADfffcdtm3bBh8fHwCApaUl3njjDaXjhISEoEePHgCA5s2bIzAw8LGLz6kSEhKCrl27AgBatWqFtm3bomvXrujYsSMAoFOnTvD09ERMTAy6deumtG92djaysrIQFBSExYsX49q1a0qL4qmTUVuvo6G5HnsWyZcrLzONGjULUmnjLPiJqP7T6AxGQUFBlS87Ozts3LgRq1ev1nZGUkPv3r2RkJCAv/76C/Pnz0dGRgaeeeYZvP766wCAmzdvQiaTwc/P77HHcXJyUnpsamqK/Pz8GudxdHRUemxiYlJtW0FBgeLxkSNH4O3tDTc3NwQHB2PMmDGQy+VVzsqok1Fbr6OhOXjpSzhNd4J7X2ssmP6m2HGIqAnT6AzGhAkTMHLkyGq3RUZGqtzWkP33g/BRj94OmpGRobLvo4Mqb926Vatc/yWRSODn5wc/Pz/MmzcPixcvxpo1a/Dee+/BwMAAAFSeWanJczyqoqKiVscEAJlMhrFjx2Lx4sVYtGgRdHR0UF5eDiMjI8jl8lofv6n4I/88rHpbIahLL1ibGogdh4iaMI0KjFGjRqk8hf+4bQ2ZiYn6i0TVVV9VBEFARUWFYoXahx4uGieTydC6dWtYWVnh119/Rfv27RV9KioqoKur/o9Bs2bNlMZtyOVyXL58GX369KnVa7h79y4ePHiAkJAQRcF2/PhxyGSyWh23Kfk77lck6MmhKwjo3/lVseMQUROnUYGhSmZmJiwsLLR5SFJDeXk5WrVqhZkzZ8Lf3x82Nja4dOkSli9fjkmTJsHIqPI2xfDwcLz66qvQ0dFBv379cOPGDezZswd79uxR+7mGDh2Kd999F71794aDgwO++OILJCUl1brAcHR0hJubG5YuXYqFCxciMTERCxcubLSDhuvC/5YtQhZy0Ktjcwye2lXsOETUxNWowAgKCqr2e6DyL9kbN27U+oOGak5fXx9RUVH44osvsGrVKty/fx8ODg5Yvnw5ZsyYoeg3e/ZsODg44Msvv8S2bdvQoUMHfPzxx4rtrq6usLOzUzr2w7ktHpo3bx4KCwuxatUqmJiYYNSoUbC0tFQaBFvdcdzc3GBra6vU5uHhoZifQiqV4tdff8Xbb7+NGTNmwMHBAZ9//jmWLl0KMzOzGmVUp09jU1DwAMf3xaGiUAZTp96cuZOIRFejtUgWLVoEAFizZg0WLlyotE1PTw/u7u6YNGmS0gdCfcS1SOhRDf3f/a2w5/H+25uhb6mL8xeT4evm+OSdiIg0UCdrkXz44YcAAFtbWyxevLh2CYlIa77f/TMAoE1XRxYXRFQvaHQelcUFUf1x+3Ysbl95AAAY8ex0ccMQEf1D40Ge9+/fx9GjR3H79u0qtymyACF6epavWQBBJsDU0QDzn+d/e0RUP2hUYPz555945plnYG1tjYSEBLRt2xbx8fEoLS1Fhw4dWGAQPUW/Hj4NAPDt4gYbs4Y3foSIGieNLpG8/vrrCA0NRXx8PADg6tWrSE1NRXBwMIYMGaLVgESkWkZmHIpNJJDoSjB+7CtixyEiUtDoDMalS5cUa29IpVKUlpbCxsYGn332GQICArBq1SqthhRDY5wsjFRrqP/ee06thcs8V7g+kGHm2NlixyEiUtB4LZKHE2rZ2dkhOTkZAGBhYYHs7GztpRPBw1kky8rKRE5CT9PDKdQfnQ21vjuW9ScAoI1hG5gYaHXePCKiWqn1b6QBAwZg0aJFmDVrFr7//nt06tSwl4fW1dWFsbExMjMzoaenV2XtEGpcBEFAUVERMjIyYGlpWWVdmfrs0pXfEXO/EIbN9dHHd47YcYiIlGhUYPx39sfVq1dj2rRpmDBhAlq1aoWIiAithRODRCKBg4MDkpKSFGdmqPGztLRscLN9/t/K13Fz20206GODIUf7ix2HiEhJjWbybCzUmYVMLpfzMkkToaen16DOXDxk7W6MnORiBE31x2+bo8WOQ0RNRJ3M5PmQg4OD0oqajZFUKm2QU0ZT03D23G/ISS4GJMC0yW+JHYeIqAqNBhiUlZUhJydH21mISE3hG5YDAGxbmmD0oGEipyEiqkqjAmPSpEn44IMPGuytfUQNXVTUVQBAx87e0OPKqURUD2l0iSQuLg6HDh3Cjh070KpVK+jr6yttj4yM1EY2IqrG1asnkZVYCACYMmGRyGmIiKqnUYHh5eUFLy8vbWchIjWErVsCCICVhxEmPDNe7DhERNXSqMBYv369tnMQkbp8c+EyzwVd4QwD3YZ39wsRNQ21nmirvLy8wc1+SNRQ5eTcQqxpBSy6WGCyxxtixyEiUkmj0WHl5eVYtmwZXFxclMZfzJ07V7EAGhFp396odZBJJHArFTC0x3Cx4xARqaRRgREeHo6dO3dizZo1Su19+vTBihUrtBKMiKpau2E70vemwz7bDoZ6vDxCRPWXRgXGt99+i+3bt2PcuHFK7f369cP+/fu1EoyIlD14kI6/j91D5k+ZcDfvJXYcIqLH0qjASE1NVdxFIpFIFO26urooLi7WTjIiUvLRF29BXibAwFoXL08NFTsOEdFjaVRgeHt749SpUwCUC4ytW7eiffv22klGREr2HvgNANDStxmszIxETkNE9Hga3UUSGhqKKVOmYMmSJQCA3bt348CBA4iIiMDOnTu1GpCIAFlFBeIuZgAA+vQaIHIaIqIn06jAmDRpEiQSCcLCwiCXyzFmzBh4e3tjy5YtGDVqlLYzEjV5W35Yg7LcCkgNpZg/M0zsOERET1Tr5dqLi4shl8thYmKirUx1Tt2lZonqiz4jWuHkvhtw87PCrb+yxY5DRE1YnS7X/lBOTg5u3rwJoHL6cCsrq9ocjohUSJdlQ2osRdcuncSOQkSkFo0GeRYUFGDmzJmws7NDt27d0K1bN9jZ2WHWrFkoLCzUdkaiJu3WnQvQH9scbT9phflz3hc7DhGRWjQqMObMmYNTp05h7969uHfvHu7du4e9e/fixIkTmDt3rrYzEjVp+6M3AQC8ZDro3cFf5DREROrR6BLJ3r17cebMGbRr107RNmzYMLi5uaFHjx747rvvtBaQqKk7kfAnBDsBLaQtlG4LJyKqzzQ6g2FpaQl7e/sq7Q4ODhyHQaRFeXlZ2P32JcQtioOHWR+x4xARqU2jAmPEiBF4++23UVpaqmgrLS3FkiVLMGLECK2FI2rqPv1qGeRlAiQVAiY/y8uPRNRwaHSJJD4+HocOHcKuXbvQvn17CIKAK1euICsrC4MHD8bIkSMVfSMjI7UUlajpiTxUubaPe1trmBsbipyGiEh9GhUYXl5eirVIHmrTpo3GIVJSUpCeng5vb2+156WQyWS4fv06jI2N4enpqfFzE9VnsZfvAgC6+3NwJxE1LBoVGOvXr9fKk5eUlGDy5Mn47bff4ObmhuTkZKxatQrz589/7H579uzByy+/DGNjYxgbG8Pe3h7bt2+Hra2tVnIR1QdHj+9C4b0yQAeYPW2p2HGIiGqkVhNt1dby5ctx9uxZJCQkwMHBAZGRkRg1ahT8/f3RrVu3avc5ceIExo4di6+++gozZswAAPzxxx9IT09ngUGNyobvPgYA2HqaoEf7LiKnISKqGY2mChcEAfv27UNUVBRycnKqbN+0aZNax7G3t8ecOXPwzjvvKNp8fX0REBCAL774otp9BgwYAD09PRw8eLCmsRU4VTg1BM6+lki9mot+Ie1wbNcVseMQEQGo46nCX3/9dXz22Wfo06cPLC0tNQqYlpaG9PR0dO7cWand398fMTEx1e5TWlqKqKgofPTRR8jPz0dcXBwcHR3h4ODw2OcqLS1VuuMlLy9Po8xET0tpaT6MB1vCxkEXo5+dKnYcIqIa06jAiIiIwLFjx1RexlBHdnblgk02NjZK7TY2Noptj8rKykJ5eTkuX74MHx8f2NvbIz4+Hj169MC2bduqHOuh8PBwLF++XOOsRE/boejvYNDeDG3bmGDm+AVixyEiqjGN5sGQy+Vo3bp1rZ5YT08PQOVAz/8qLi6Gvr7+Y/c5fPgwLl68iAsXLiApKQlJSUlYtGiRyucKDQ1Fbm6u4islJaVW2Ynq2unEAwAA73ILGBlU/98DEVF9plGBMXr0aLXHWaji4uICqVSK1NRUpfbU1FS4urpWu4+trS1MTEwwevRo2NnZAag84zF27FicPHlS5XMZGBjA3Nxc6YuoPvv5twsouFYAT2OunkpEDZNGl0jee+89tGnTBtu2bUOLFlXXR9ixY8cTj2FsbIyePXti3759mDJlCgCgsLAQR44cwbJlyxT94uPjkZ+fj06dOkEqlSIwMLBKUXLnzh00a9ZMk5dCVO9cuhqFv3+snP+iSxRnxiWihkmjAmPBggWQyWRwd3eHmZmZxk8eFhaGwMBAhIaGokePHli3bp1i2feHVq5ciTNnzuDq1asAgBUrViAgIADvvPMOAgICEB0djW3btuHHH3/UOAdRfbLum5UAAHMXQwzqGShyGiIizWhUYPz888+IiopCp061O33bt29fHDt2DBs2bMDZs2fh6+uLLVu2wNTUVNHHy8sLZWVlise+vr44ffo0Pv74Y6xevRqurq44fvw4evbsWassRPXFyTPnAADebey5eioRNVgazYPh5OSEa9euwcLCoi4y1TnOg0H1lVwmg5GNIcpyK/B62EysfusrsSMRESlR9zNUo0GeQ4YM0dp04UT0rx9++hxluRWQ6Esw74V3nrwDEVE9pdElktu3byMiIgLbtm1Dy5Ytq5zG5QqqRJrZsvtbAIC9lxlc7Z3FDUNEVAsaFRitWrVCq1attJ2FqMm7GpsAAPBt30LkJEREtSPqaqpE9K/SkgLYzXWEfqoNpnaaLXYcIqJa0WgMxn+Vl5drIwdRk3f03DaU6kph72yAMUOmix2HiKhWNCowysvLsWzZMri4uChN6z137lzEx8drLRxRUxL9z/TgLcvNYKhiunwiooZCowIjPDwcO3fuxJo1a5Ta+/TpgxUrVmglGFFTE7Hud9zZeAdW+S5iRyEiqjWNCoxvv/0W27dvx7hx45Ta+/Xrh/3792slGFFTcu9uElIu5uLB6Qfo7D1U7DhERLWmUYGRmpoKLy8vAFC6RVVXVxfFxcXaSUbUhKz/5j1ABhja6mHC0OfEjkNEVGsaFRje3t44deoUAOUCY+vWrWjfvr12khE1IUdO/A4AcPWxgq6ujshpiIhqT6PbVENDQzFlyhQsWbIEALB7924cOHAAERER2Llzp1YDEjUFN65Xrp7q16mjuEGIiLREowJj0qRJkEgkCAsLg1wux5gxY+Dt7Y0tW7Zg1KhR2s5I1KjF3YjBg5QSAMDU8a+InIaISDs0KjAiIyMxceJETJw4EcXFxZDL5TAxMVFsGzlypDYzEjVqGyLCAQAmTgYIDuAATyJqHDQqMEaNGoWHi7AaGRmp3EZET3brQSwMXQzh4W3L5dmJqNHQqMBQJTMzs8Eu4U4kFqmfHC17tMTzxsFiRyEi0poaFRhBQUHVfg8AcrkcN27cQJ8+fbSTjKgJSLsbi4R/Ju0c0nmyuGGIiLSoRgVGu3btAAAHDx5UfP+Qnp4eQkJCMGnSJO2lI2rkdh3dCFm5ABeJBL6eHcSOQ0SkNTUqMD788EMAgK2tLRYvXlwngYiakq+//hHXo9JgHdICmCl2GiIi7dFooi0WF0TakRiXBaFCQGs3X7GjEBFpVa2XaycizZyNOYGCtFIAwMxJC0VOQ0SkXSwwiESyafvHAAAzRwP4dwgQOQ0RkXZpVGBERkZqtI2I/nXmr3MAAE/vZiInISLSPo0KjMdNB86pwonUk3QjCwDQ07+byEmIiLRPq5dIONEWkXqiLvyBgtTK8RcvTeX4CyJqfDjRFpEI/riyG/bj7WF0rwzt2/UQOw4RkdZxoi0iEaSW/Q3bYFsEldmLHYWIqE5woi0iESRJMgEAbZt3FzkJEVHdqPVEW+Xl5VoLQ9QUnLtyBlejc1CWVYZBXXnGj4gaJ40KjPLycixbtgwuLi7Q19dXtM+dOxfx8fFaC0fUGH2z7SOkbkpF6oe34OzQWuw4RER1QqMCIzw8HDt37sSaNWuU2vv06YMVK1ZoJRhRYxV9PhoA4OFtK3ISIqK6o1GB8e2332L79u0YN26cUnu/fv2wf/9+rQQjaqwS4zIAAP5d/EROQkRUdzQqMFJTU+Hl5QUAkEgkinZdXV0UFxdrJxlRI/R37EXkppQAAJ6fME/kNEREdUejAsPb2xunTp0CoFxgbN26Fe3bt9dOMqJGaNO2DwABMG6ujz5dh4gdh4ioztToNtWHQkNDMWXKFCxZsgQAsHv3bhw4cAARERHYuXOnVgMSNSanok8DAFy9rUVOQkRUtzQqMCZNmgSJRIKwsDDI5XKMGTMG3t7e2LJlC9ciIXqM+BvpAAC/TjzTR0SNm0YFBgBMnDgREydORHFxMeRyOUxMTLSZi6jRybyfBudXXWF6swjTJ3L8BRE1bhoXGA8ZGRlpIwdRo3fs7FYIlnpw8zPHwG7DxI5DRFSn1C4wBg0apPZBjxw5olEYosbsSupJQA9oIbNQGhxNRNQYqV1gdOnSpS5zEDV6u3acQba5FPaDOogdhYiozqldYKxcuVLxfWRkJEaOHFltv8jIyNpmImp0srLSceP4fQgVAtpP7yt2HCKiOqfRPBiPu1OEd5EQVbVp6wcQKgTometibNDzYschIqpzGhUYqmRmZsLCwkKbhyRqFI6cPAgAsG9pBj09/Sf0JiJq+Gp0F0lQUFC13wOAXC7HjRs30KdPH+0kI2pEYmNvAwDatm0hchIioqejRgVGu3btAAAHDx5UfP+Qnp4eQkJCMGnSJO2lI2oEioqKcC+pAADwzMARIqchIno6alRgfPjhhwAAW1tbLF68WGshUlJSkJ6eDm9vb5ibm6u9X25uLq5cuQInJyd4eHhoLQ+RNm3+8VPIiuWQGkoxY9z/xI5DRPRUaDQGQ1vFRUlJCUJCQuDj44OpU6fC3t4e69atU2tfQRAwefJk9O3bF5988olW8hDVhVMxByE1lKJZC1MYG5mJHYeI6KnQaCbPXr16PXZ7VFSUWsdZvnw5zp49i4SEBDg4OCAyMhKjRo2Cv78/unXr9th9P/74Y8hksiqXaojqG+PWRWi9oTUGFrmJHYWI6KnRqMB4dFZPuVyOmzdvYufOnZg1a5bax4mIiMCcOXPg4OAAABg5ciTatWuHiIiIxxYYf/31Fz766COcP3++ymBTovqkokKGBN0CSHSk6NMuWOw4RERPjUYFxrJly6ptDw4Oxm+//abWMdLS0pCeno7OnTsrtfv7+yMmJkblfvn5+ZgwYQI2bNgAe3t7tZ6rtLQUpaWlisd5eXlq7UdUW9GXjiBdTwodQUBg14lixyEiemq0Og/GqFGjcPToUbX6ZmdnAwBsbGyU2m1sbBTbqjNnzhwMHDgQI0aoPxo/PDwcFhYWii8XFxe19yWqjXUbVyHuzTiU/5wNcxObJ+9ARNRIaLXAuHLlitqLOOnp6QGoHOj5X8XFxdDXr34ion379mH//v0YPXo0oqKiEBUVhcLCQqSlpT123EdoaChyc3MVXykpKWq+IqLauXT1b5Sll8GkhKsOE1HTotElkpkzZ1Zpy8nJwcGDB/HSSy+pdQwXFxdIpVKkpqYqtaempsLV1bXafSoqKtCuXTusWLFC0Xb37l0UFRUhLS0Nx48fh46OTpX9DAwMYGBgoFYuIm2RywXcjn8AAOjTPUDcMERET5lGBUZBQUGVNjs7O2zcuBETJkxQ6xjGxsbo2bMn9u3bhylTpgAACgsLceTIEaUxHvHx8cjPz0enTp0wevRojB49Wuk4HTt2RL9+/bB27VpNXgpRnfn91AEUZZQBEmDWlEVixyEieqo0KjB27NihlScPCwtDYGAgQkND0aNHD6xbtw52dnZKd6KsXLkSZ86cwdWrV7XynERPy7Y9nwMAzJwM4e3RXuQ0RERPl1bHYNRU3759cezYMSQnJ+OTTz5B27ZtERUVBVNTU0UfLy8v+Pn5qTxGp06d4Onp+TTiEtXIhUuVd0O5e9uKnISI6OmTCIIgaLLj77//jtWrV+P69esQBAFt2rTBG2+8gQEDBmg7o9bl5eXBwsICubm5NZqanEhdgiDAys0YuSklmPXGCHy5KlLsSEREWqHuZ6hGZzC2bt2KIUOGwNLSEq+++ioWLlwIS0tLDBkyBFu3btU4NFFjceNWPHS9jGDgaIDp4+aKHYeI6KnT6AxGq1at8Oabb2L69OlK7REREVi1ahViY2O1FrAu8AwG1bWNe97Fuvwf0bxCwJEXOH6IiBqPOj2DkZiYiJCQkCrtISEhSExM1OSQRI1KXNZpAEBLWIuchIhIHBoVGI6OjtVObHXixAk4OjrWOhRRQ3clNRHyCjnaWXUUOwoRkSg0uk31lVdewcSJEzF//nz4+/sDAKKjo7Fu3TqV65QQNRWxt5Nx+P2bgBRodeQdseMQEYlCowLjtddeg6WlJVavXo2VK1cCAFq2bIm1a9dixowZWg1I1NBs/eEjCBUC9Mx0MKBn1UuJRERNgUYFBgDMmDEDM2bMQHl5OYB/1xYhauqizlQu+OfY0hzSaqauJyJqCjQuMB5iYUGk7GZc5WJ67dt5iZyEiEg8os7kSdTY3M15gPTEyrV6Qoby8ggRNV0sMIi0aPOPn0JWJIfUQILxI+eJHYeISDQsMIi06NiJfQAAO09TGBoai5yGiEg8tR6DQUT/MmhRguZjmqOPPRfgI6KmTe0Co1evXmoftLpJuIgauwcFpchxE9CsRTO81PpFseMQEYlK7QIjKChI8f3du3fx2WefoW/fvujatSsA4Ny5czh+/DjmzuXCTtQ0/RHzK3J0pdATBPTsMEbsOEREolK7wFiyZIni+1GjRuGTTz7BK6+8otTnk08+wfHjx7WXjqgB2X8oAg+yHqCTpwUMDMzEjkNEJCqNVlO1tbVFYmJilVXUcnNz0aJFC2RlZWktYF3gaqpUF3x6NUfcqQz0GNoSp3+5KXYcIqI6UaerqQLA2bNnq22TSCSaHpKowSosrcCdhAcAgMDefcUNQ0RUD2h0F8lLL72EcePGYcGCBejatSsEQcD58+exdu1azJ8/X9sZieq9fb9HouheGSABZk5+Xew4RESi06jAePfdd+Hs7IyPP/4Yy5cvB1C52Fl4eDhmzZql1YBEDcG+X78GAJg7GcLFxUfkNERE4tOowJBIJHjppZfw0ksvoaysDBKJhGuSUJN25colAEDLVs1FTkJEVD/UeiZPFhfU1JWUy3A7IRsA0K9nD5HTEBHVDxoVGOXl5Vi2bBlcXFygr6+vaJ87dy7i4+O1Fo6oIfg9+gTyU0sBADMnvSpyGiKi+kGjAiM8PBw7d+7EmjVrlNr79OmDFStWaCUYUUNxPfkAfD7ygd/Lrmjt4y92HCKiekGjAuPbb7/F9u3bMW7cOKX2fv36Yf/+/VoJRtRQ3Mo7Dz0rPfTu7CV2FCKiekOjAiM1NRVeXpW/TP8774Wuri6Ki4u1k4yoASiXyZEqzQAAdHLoJnIaIqL6Q6MCw9vbG6dOnQKgXGBs3boV7du3104yogbg1MUYHP3sFjL3Z6JPuxCx4xAR1Rsa3aYaGhqKKVOmKNYn2b17Nw4cOICIiAjs3LlTqwGJ6rMfd69F3uUClKWVwt2pk9hxiIjqDY0KjEmTJkEikSAsLAxyuRxjxoyBt7c3tmzZglGjRmk7I1G99VdM5Zk8V28rTpNPRPQfGhUYADBx4kRMnDgRxcXFkMvlMDEx0WYuonqvXCZHYnw6AKB75w4ipyEiql80LjAeMjIy0kYOogbnzLXruH+rclDz82NeFDkNEVH9onaBMWjQILUPeuTIEY3CEDUkP+z6CEKFAH0LXfTrxQGeRET/pXaB0aVLl7rMQdTgnL9wAgDg4m0JqbTWs+4TETUqahcYK1eurMscRA1KWYUcOWW5kOhL0N2vndhxiIjqnVqPwSBqis4l3ILJRDu0HtcMS3otETsOEVG9o3GBcf/+fRw9ehS3b99GRUWF0rbFixfXOhhRfXb64g8ol0jQHAJ8vPqLHYeIqN7RqMD4888/8cwzz8Da2hoJCQlo27Yt4uPjUVpaig4dOrDAoEYvIfMEYAa0gRUkHH9BRFSFRr8ZX3/9dYSGhiqWZr969SpSU1MRHByMIUOGaDUgUX1TUi7DTxtP4+b/3YRRqo3YcYiI6iWNCoxLly5h1qxZlQeQSlFaWgobGxt89tln2Lp1q1YDEtU3p28kISO+AKVppejbIVjsOERE9ZJGBUZBQQEsLCwAAHZ2dkhOTgYAWFhYIDs7W3vpiOqhPb+sg7xYDh0DKYYEzxA7DhFRvVTri8cDBgzAokWLsH//fsydOxedOnHBJ2rcLpyrnEjOqaUZ9PT1RU5DRFQ/aVRgfPzxx4rvV69ejeLiYkyYMAE3b97El19+qbVwRPVNcZkMSXF3AADdu7QROQ0RUf2l0V0kCxYsUHzv5OSEw4cPaysPUb126no8MuILAABTRz8vchoiovqL99cR1cCen9dCXiyHrpEUQRx/QUSkkkYFxo8//ojJkydXaZ88eTJ27dpV61BE9VVWyVVY9bFC62520NXTEzsOEVG9pdElkmXLlmHPnj1V2t9++22MHTsWY8aMqdHxUlJSkJ6eDm9vb5ibmz+xv0wmQ1xcHHR1deHh4QFdXc54TnWvsLQC5c55cJrhhNdtB4odh4ioXtPoDEZiYiLs7e2rtDdv3lwx+ZY6SkpKEBISAh8fH0ydOhX29vZYt27dY/d577334OTkhJCQEAwePBju7u7Yv39/jV8DUU2dvpGAW4YyAED/9uNETkNEVL9pVGC0atWq2jMYu3btgpeXl9rHWb58Oc6ePYuEhARcv34d27ZtwyuvvILo6Ohq+8tkMhQXF+PatWu4du0akpKSMHPmTIwfPx737t3T5KUQqe3w8S+Qf6sEdmVyODl1FzsOEVG9JhEEQajpTtu2bcMLL7yAOXPmoE+fPhAEASdOnMDnn3+Ob775BpMmTVLrOPb29pgzZw7eeecdRZuvry8CAgLwxRdfqHWM9PR02Nvb45dffsHQoUPV2icvLw8WFhbIzc1V65IMEQD4D/fCuf3x8Oluh9g/08WOQ0QkCnU/QzUavDBp0iTIZDKEhYUp5sTw9vbGV199pXZxkZaWhvT0dHTu3Fmp3d/fHzExMWpnOXfuHACgRYsWKvuUlpaitLRU8TgvL0/t4xMBQF5JOW7dvAsA6O7nK3IaIqL6T+PRkVOnTsXUqVNRWFgIiUQCY2PjGu3/cEpxGxvlxaJsbGzUnm48KysL8+fPx7hx4+Dj46OyX3h4OJYvX16jfET/dfTCZWQlFgEAnh87S+Q0RET1X63nwTAxMUF6ejouXLiAsrIytffT++cWv5KSEqX24uJi6Ksx/XJubi6CgoJgb2+PTZs2PbZvaGgocnNzFV8pKSlq5yQCgMifPoBQLkDfQhf9+o4VOw4RUb1XowIjNTUVixYtUmpbuHAhPD090blzZ7Rp00ax8NmTuLi4QCqVIjU1tcpzuLq6PnbfvLw8DB48GDo6Ojhw4ADMzMwe29/AwADm5uZKX0Q1cfHCaQBAy3a2kEgkIqchIqr/alRgfPDBB/Dw8FA8Pn/+PD766COEh4fj5MmTcHV1VftShLGxMXr27Il9+/Yp2goLC3HkyBEEBgYq2uLj45XGZDwsLgDg0KFDilVdiepK8v1C3LqRCQAI6tdb5DRERA1DjcZg7N+/H0ePHlU8/vnnn9GlSxcsXrwYAPDpp59i2LBhah8vLCwMgYGBCA0NRY8ePbBu3TrY2dlh1qx/r3GvXLkSZ86cwdWrV1FeXo7g4GAkJSXhm2++wZUrVxT9vLy80Lx585q8HCK17Dv6E/JSKi/lzX7uTZHTEBE1DDUqMNLS0tCsWTPF4z///BP9+vVTPG7ZsiXS09W/fa9v3744duwYNmzYgLNnz8LX1xdbtmyBqampoo+Xl5dibEdRUREkEgm8vLwQHh6udKzFixfXqLghUldy1kF4hHrALKkU3t6dn7wDERHVrMBwdnbGmTNnMGDAAOTn5+P06dN4+eWXFdtTU1Ph5ORUowABAQEICAhQuf3NN//9i9HCwgJRUVE1Oj5RbVTI5MiQX4eJjwkmtGdxQUSkrhoVGFOmTMHkyZMxffp0nD59GiYmJhg0aJBi+4kTJzBgwACthyQSS0xyJm4ZFQKQYkDr4WLHISJqMGpUYISGhiI/Px87d+6EtbU1du3aBRMTE8X2PXv2YOXKlVoPSSSWvb98jr8j78HG1xR+E3l7KhGRujSaKryh41ThpK6+Y31xYtdVOLYxR+rfuWLHISISnbqfobWeaIuoscouLMPN2Mp5Xbp3bS1yGiKihoUFBpEKh/6KQfrNQgDAtLEzRU5DRNSwsMAgUmFf5CrIS+XQt9DFM8EzxI5DRNSgsMAgqka5TI5LF/4EALTuaA+plP+pEBHVBH9rElXjXGImkq7dBwCMGRoschoiooaHBQZRNQ4d/w4yqQQSXQlmT39H7DhERA1OjebBIGoqMoqOw3uVN9pl6aFZs5rNTktERDyDQVTF7ftFuKd/GwDwTNu+IqchImqYWGAQPeJg9Enc1JUBAAZ2mi5yGiKihokFBtEj9u1djWvzb6Do+3twaN5e7DhERA0SCwyi/8gpLMP1q5cgL5HDQc9S7DhERA0WCwyi//jtr2tIvvoAADA1ZIK4YYiIGjAWGET/8dOed1FRIIOuqQ6mT1kidhwiogaLBQbRP/JKynHprygAQLvODtDXNxA5ERFRw8UCg+gfhy4lIPFKNgDg+bHjRE5DRNSwscAg+sfeXe+iPLcCOkZSzJq2XOw4REQNGmfyJAJQWFqBMrNY2I22g5dgCmMTU7EjERE1aCwwiAD8fu020p1KYOdmh5Xt5okdh4ioweMlEiIAUee+RY6uFOYyAd06PC92HCKiBo8FBjV52YVl+DNqB3Kjc9GlwgZ6uoZiRyIiavB4iYSavMjovxH9cxLKHlSg+QdDxY5DRNQo8AwGNXmRu5ag7EEFdE10sGDOh2LHISJqFFhgUJN2Mz0fl8/9CQDo2suNd48QEWkJCwxq0nYcOYjblx8AABbNfV3cMEREjQgLDGqyZHIBR/a9B6FMgJmDAUYNny12JCKiRoMFBjVZJ+My8PeFOABAUJAfJBKJyImIiBoPFhjUZO05uhnFEgASYOmij8SOQ0TUqPA2VWqSkrIKkV+xF55vecIvxxDt2nQXOxIRUaPCMxjUJH3/RxSumucCAOb2ny9yGiKixocFBjU5eSXlOB+1AgX5crSQ6cC//VSxIxERNTosMKjJ+eHPGzj8w3ncWHgDHe634eBOIqI6wAKDmhSZXMDeba+h6F4ZpHoSzJ+xRuxIRESNEgsMalL2X0rFX6dOAwB6B3rD1sZB5ERERI0TCwxqMsoq5IjY8n/IuF4ISIA1yz8TOxIRUaPFAoOajB3RifjzwB4AQLsejvDr2F/kREREjRcLDGoSCksrsO3buci4XgiJjgTfrt8idiQiokaNE21Rk/D1H1eQUnwFEj0J+ga3QudOA8SORETUqLHAoEbv9v0iXLmyHNJAW3TtaI7NE/aJHYmIqNFjgUGNWrlMjre37Uai9U0AUizoMAwuzi3FjkVE1OhxDAY1ap8c+htn9r6G1PhitJDrYEz/98SORETUJPAMBjVaZxLvIzJiMuKOZUFy8j62nPoRerqGYsciImoS6sUZjJSUFJw/fx55eXl1ug81HWeTsvH2uqX486erAIAJ03ohwH+MyKmIiJoOUQuMkpIShISEwMfHB1OnToW9vT3WrVun9X2oafnl8l28vvEbnNn1NeRlAlzaWmLL57+LHYuIqEkR9RLJ8uXLcfbsWSQkJMDBwQGRkZEYNWoU/P390a1bN63tQ01D2oNibDyRiCNnjiH5RBgK7pRC31wXB3cfgo4urwYSET1Nov7WjYiIwJw5c+DgULkexMiRI9GuXTtERESoLBY02YcaJ0EQkPqgGBdTHuBYbCZ+vpiCdtIYFKV8hjuncwEJsOWbL9Dap6vYUYmImhzRCoy0tDSkp6ejc+fOSu3+/v6IiYnR2j4AUFpaitLSUsVjbY/biIrZj4jfVuHS8RSVfTx8m8G1lTUAID+7BBeOJqvs69bGBu5tbQEAhXmlOH/wlsq+zt5WaNGhGQCguLAcZ39NqtwgVO3r4GkB787NAQgoK6nAnz8nqjxuczdztPK3BwDIKuSI2huvsm8zZ1O06fHvomHHd95U2dfa3hi+vZ0Uj6P2xkNWUU1YAJbNjNChrzMAAQKA0/sSUVYqgyAIEARALgByoXJfEysdtA4ywU0DHdw5XgIAWLR4LsaFvKAyCxER1R3RCozs7GwAgI2NjVK7jY2NYps29gGA8PBwLF++vDZxHysh7SJ+z7uDW4dUFw13LEvQzLfyg6+osAiJj+l726gIcR3LAFSOOYl/TN9kaSESulQAAMrKyxD3mL7WA6xxq4ccAFBRUoHYw7dV9rUMsMTtXpXfywU5rj2mr3kXc9zp++9wnquP6Wva3hR3B/77Y/f30dsQyqovMIx9jJE+WE/x+PqJFMjyZdX2NfIwQtEIcwDAgBe90X/aSMyd/qbKHEREVLdEKzD09Co/OEpKSpTai4uLoa+vr7V9ACA0NBSvvfaa4nFeXh5cXFw0yl2dFo4dMdDCFZeGSFT28XRtBpeKysIoz7QY1kOq/1AFADdPW7hX/HMGw7AUlkMe/VD993lcvK3hKbMDABTrl8E8qELlcR1bWMJLVnlWokxHBtOgMpV9m7tZoJXMEQAggxzGQSUq+zZzMkcbmbPisWFQkcq+1vam8JW7Kh6bBJZAJpNX29eymQk6yN0AABKJBBYDy1FRWgGJRAKpRAKpFNCRSiGRAFY2Jpjdagq6tJsEG5PmKp+fiIieDtEKDBcXF0ilUqSmpiq1p6amwtXVVWv7AICBgQEMDAxqH1qFXp2GoVenYUBoDXZ6vQZ9X3tyF4VXatB3bg36zqpB3xk16Du9Bn2n1aAvERGJSrTbVI2NjdGzZ0/s2/fvuhCFhYU4cuQIAgMDFW3x8fGK8RXq7kNERETikgiCoPpcfR07fvw4AgMDsXDhQvTo0QPr1q1DUlISLl68CFNTUwDAzJkzcebMGVy9elXtfZ4kLy8PFhYWyM3Nhbm5eZ29PiIiosZG3c9QUSfa6tu3L44dO4bk5GR88sknaNu2LaKiopQKBS8vL/j5+dVoHyIiIhKXqGcwxMIzGERERJppEGcwiIiIqHFigUFERERaxwKDiIiItI4FBhEREWkdCwwiIiLSOhYYREREpHWiLtculod35mp7VVUiIqLG7uFn55NmuWiSBUZ+fj4AaHXBMyIioqYkPz8fFhYWKrc3yYm25HI50tLSYGZmBolE9QqoNfFwhdaUlBRO3qUFfD+1j++pdvH91D6+p9pVV++nIAjIz8+Ho6MjpFLVIy2a5BkMqVQKZ2fnJ3fUgLm5Of/D0CK+n9rH91S7+H5qH99T7aqL9/NxZy4e4iBPIiIi0joWGERERKR1LDC0xMDAAO+88w4MDAzEjtIo8P3UPr6n2sX3U/v4nmqX2O9nkxzkSURERHWLZzCIiIhI61hgEBERkdaxwCAiIiKta5LzYGhbfn4+bty4AVtbW7i7u4sdp8GpqKjA+fPnYWVlBR8fn2r7FBYWIjY2FlZWVvD09HzKCRuWiooKxMXFwcDAAO7u7tDR0am2361bt5CVlYVWrVrB1NT0KadsWAoLCxEXFwcbGxu4urpW26e0tBR///03TExMVP4ck7KCggJcvHgR9vb2aNmyZZXtqampuHv3Llq0aAErKysREjYMUVFRVdpatGgBBwcHpbaKigr8/fff0NXVRZs2bbQ20aRKAtXKV199JRgbGws+Pj6CiYmJEBwcLBQUFIgdq0EoKCgQ3n77bcHFxUUwMzMTxo8fX22/77//XjAzMxO8vb0FMzMzoX///sKDBw+ectr6r6KiQli6dKnQrFkzoU2bNoKLi4vg7u4uHD58WKlfQUGBEBwcLJiYmAg+Pj6CsbGxsGnTJpFS12/Z2dnCzJkzBWtra8HPz0+wtrYWOnbsKFy7dk2p36+//irY2NgInp6egpWVldC5c2chLS1NpNQNx7hx4wSpVCq88MILSu1lZWXCpEmTBENDQ6F169aCoaGhsHLlSpFS1n8AhHbt2gkBAQGKr927dyv1OX36tODk5CS4uLgIdnZ2QqtWrYQbN27Uba46PXojd/nyZUEqlQrbtm0TBEEQMjIyBHd3d2H+/PkiJ2sYkpKShGXLlgmpqanCM888U22BcfPmTUFPT0/YuHGjIAiCkJOTI/j4+AjTp09/2nHrvfz8fGH58uWK4ksulwuLFi0SzM3NhZycHEW/efPmCZ6enkJmZqYgCIKwZcsWQSqVCleuXBEjdr127do1Yfv27YJMJhMEQRBKSkqEwMBAISAgQNEnMzNTMDMzE1asWCEIgiAUFxcL3bt3F4KDg0XJ3FB8+eWXQu/evYUePXpUKTDCwsIEOzs74datW4IgCMKhQ4cEiUQiHD16VIyo9R6AKn9I/FdRUZHg6OgozJ07VxCEyj9Ghg4dKvj5+dVtrjo9eiP32muvCS1btlRqW7lypWBhYSFUVFSIlKphUlVgLF26VHBwcBDkcrmibf369YKhoaFQVFT0NCM2SPHx8QIA4fjx44IgCEJ5eblgbm4ufPjhh0r9PDw8hIULF4oRscFZunSp4O7urnj82WefCcbGxkJhYaGibdeuXYJEIuFZDBWuXr0qODg4CLdu3RICAgKqFBienp7CokWLlNq6d+8uTJ48+WnGbDAACJs3bxbOnz8v3L9/v8r2PXv2VPl5jIqKEgAIMTExdZaLgzxrISYmBp07d1Zq8/f3R25uLhITE0VK1bjExMTAz89P6Vqhv78/SkpKEBsbK2KyhuHcuXOQSCSKcSuJiYnIy8ur8nPbtWtXxMTEiBGxQbh8+TKOHz+OL774Al9++SWWL1+u2BYTE4PWrVvD2NhY0ebv7w9BEHDx4kUR0tZvxcXFGD9+PD788EO4ublV2Z6Xl4fExMRqf7fyZ1S1V199FTNmzICjoyNGjhyJzMxMxbaYmBg4Ojoqjcnw9/dXbKsrLDBqITs7GzY2NkptDx9nZ2eLEanR4XusubS0NLz22muYPn26YnG/h+9Zde8p30/VNmzYgDfeeAOhoaHo3LkzBg4cqNjGn9GaWbBgATp27IhJkyZVu50/ozW3ceNGZGZm4tKlS4iLi8ONGzcwc+ZMxfbqfkb19PRgZmZWp+8pC4xa0NPTQ0lJiVJbcXExAEBfX1+MSI0O32PNZGVlYciQIfDx8cH69esV7Xp6egBQ7XvK91O1L7/8EtHR0UhNTYWBgQEGDx4MuVwOgD+jNfH777/j+++/x8SJExEVFYWoqCjk5eUhPT0dUVFRkMlk/BnVwIsvvqg4y+vq6oq33noLP//8MwoLCwFU/zMKVL7Hdfme8jbVWnBzc0NqaqpS28PHqm5lo5pxc3NDfHy8Uhvf48e7f/8+Bg0aBGtra+zfvx9GRkaKbQ9PSaempqJr166K9tTUVL6fajA2NsbLL7+MQYMG4fbt23B3d4ebmxvOnTun1I8/o9UrLS1Fx44dER4ermhLTk5GRkYGFi9ejAMHDsDe3h4GBgbV/m7l+6me5s2bQxAE3L17Fy1btoSbmxvu3bsHuVwOqbTyvEJGRgbKy8vr9D3lGYxaCAwMxIkTJ5Cbm6to++mnn9CpU6cqp6NIM4GBgYiOjkZGRoai7aeffoKXl1e112+buuzsbAwaNAgWFhb49ddfYWJiorTd1tYWHTt2xL59+xRtOTk5OHHiBAIDA5923Hrv4V+A/xUfHw8dHR1YWloCqPwZTUhIwLVr1xR9fvrpJ1hbW8PPz+9pRW0QgoODFWcuHn75+vpi2LBhiIqKgqmpKXR0dNC/f3+ln9HS0lIcOHCAP6PVqO5n9NChQzAzM1MUD4GBgcjLy8Mff/yh6PPTTz9BX18fffr0qbtwdTZ8tAkoLi4W2rRpI/Tp00eIjIwUli5dKujo6Ai//fab2NEajFOnTgknT54UevbsKQwcOFA4efKkEB0drdheXl4u+Pn5Cd27dxf27NkjvPfee4KOjo6wa9cuEVPXT0VFRYKfn5/g5OQkHDhwQDh58qTi6+EtqYJQOWeDjo6O8M477wiRkZFC7969hTZt2gjFxcUipq+f3n33XWHatGnC9u3bhYMHDwrvvfeeYG5uXuWOm8GDBwtt2rQRdu7cKXzyySeCgYGB8Nlnn4mUumGp7i6Ss2fPCgYGBsKCBQuEffv2Cc8884zg7Oxc7R0STd2mTZuEkSNHCps3bxZ+/fVX4dVXXxV0dXWF9evXK/WbOnWq4OrqKnz//ffCV199JVhYWAhLly6t02xcTbWWsrKysGrVKsTExMDGxgYvvfQS+vfvL3asBmPgwIEoLS1VarOxscFPP/2kePzgwQOsWrUK586dg5WVFWbOnIkhQ4Y87aj13r179zBmzJhqty1fvlxpYOKxY8fwxRdf4P79++jUqRPefPNN2NraPq2oDcru3buxd+9eZGZmwtXVFePGjavyl3RRURE++ugjnDhxAsbGxpg8eTLGjh0rUuKGZe7cuXBzc8Obb76p1H7u3Dl8+umnSEtLQ+vWrbF48WLFYGVSdvDgQWzfvh337t2Dh4cHZsyYoXQJFADKy8uxfv16HDx4ELq6uggJCcG0adPqdDZPFhhERESkdRyDQURERFrHAoOIiIi0jgUGERERaR0LDCIiItI6FhhERESkdSwwiIiISOtYYBAREZHWcS0SInqse/fu4cSJExg3btxTf+7ExEScPXsWANC2bVv4+vpqfCxNXkdWVhaOHDkCAHBxcUFAQIDGz0/U1LDAIGqi7t+/j8OHDz+2T+fOnZGQkIBJkyaJUmD8/vvv+N///ofhw4dDKpXWqsC4ePFijV/H/fv3ERkZifPnz6Njx44sMIhqgAUGUROVk5ODyMhIxePz588jOzsbgwcPVrRZW1vDwcEB48ePFyFhpWbNmmHHjh21Po4mr8PHxwc7duzASy+9hKysrFpnIGpKWGAQNVEtW7ZU+uB+6aWXcP78+Sof5vfu3cOIESMUj+/cuYOzZ89i5MiRuHLlCpKTk9G+fXu4u7tDJpMhOjoa2dnZ6Nq1K5o3b17lee/evYvz58/DzMwMfn5+MDc3r1FuTZ+/efPm1b6OUaNG4cqVK7h9+zbatm0LDw+PGuUhouqxwCCix3r00sKZM2fw/PPPo127dtDX14dUKsWpU6ewdu1abN68GUZGRpDJZLhy5QqOHDmitOjSu+++izVr1qB79+4oKipCbGwstm3bhkGDBqmdR9PnV/U6evbsicLCQpiamuL48eP44osvMH36dC2+g0RNE+8iIaIaKyoqwvTp03Hy5EkcP34c06ZNw/z58zF37lwcP34cUVFRCA4ORnh4uGKfn3/+GevXr8fly5dx4MABnDhxAuHh4Zg6dSpKSkrq/PlVHWfkyJE4ffo0Dh06hPfffx+LFy/W6D0hImUsMIioxnR0dPDCCy8oHvfo0QOGhoZ47rnnlNri4uIUjyMiItC2bVucPXsWO3fuxI8//ghdXV3cu3cP169fr/PnV3WcWbNmKR7369cPGRkZePDgQY3yEFFVvERCRDVmamoKPT09xWMDAwNYWFhAKpUqtf33zMStW7dQWlqKXbt2KR1r/PjxSvvV1fOrexwANT6jQkRVscAgoqfC3Nwcrq6u2Lx5s9hRiOgp4CUSInoqgoKC8NNPP+Hu3btK7ampqSIlIqK6xDMYRPRU/O9//8P+/fvRrVs3zJ07F9bW1rhw4QJOnDiBa9euiR2PiLSMZzCICADQtWtXDBkypEr7oxNUubi4ICQkRKmPu7s7Ro4cqdTWsmVLDBs2TPHYyMgIx44dw/vvv4+kpCRcuHABXbp0wYULFx6bq7CwEDt27MCVK1dq9fzqvA5LS0uMHz8eRkZGACqnCt+xYwcSEhIem5GIqpIIgiCIHYKIqDrHjh3Dl19+CQAYPXr0U5+u/MaNG3jnnXcAAN27d8eCBQue6vMTNWQsMIiIiEjreImEiIiItI4FBhEREWkdCwwiIiLSOhYYREREpHUsMIiIiEjrWGAQERGR1rHAICIiIq1jgUFERERaxwKDiIiItI4FBhEREWkdCwwiIiLSuv8H68rH85hdcfQAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fluid = base.submodels[\"fluid\"]\n", "storage = EnthalpyStorage(param, SensibleMaterial(\"Capsule\"), geometry, label=\"Capsule\")\n", "sensible = pytes.ComposedModel(\n", " {\n", " \"fluid\": fluid,\n", " \"storage\": storage,\n", " \"coupling\": WallCoupling(param, param.fluid, \"Capsule\", geometry),\n", " \"device energy\": DeviceEnergy(param, fluid, [storage]),\n", " },\n", " param,\n", " name=\"Sensible packed bed\",\n", " defaults_from=base,\n", ")\n", "pv = base.default_parameter_values\n", "pv.update(\n", " {\n", " \"Capsule conductivity [W.m-1.K-1]\": 1e3,\n", " \"Capsule density [kg.m-3]\": 1000,\n", " \"Capsule specific heat capacity [J.kg-1.K-1]\": 2000,\n", " \"Heat transfer coefficient [W.m-2.K-1]\": 300,\n", " },\n", " check_already_exists=False,\n", ")\n", "h, R, eps = 300, pv[\"Capsule radius [m]\"], pv[\"Bed porosity\"]\n", "Z, A = pv[\"Bed length [m]\"], pv[\"Bed cross-sectional area [m2]\"]\n", "rho_f = pv[\"Heat transfer fluid density [kg.m-3]\"]\n", "c_f = pv[\"Heat transfer fluid specific heat capacity [J.kg-1.K-1]\"]\n", "a = 3 * (1 - eps) / R\n", "u = pv[\"Mass flow rate [kg.s-1]\"] / (rho_f * eps * A)\n", "xi = h * a * Z / (eps * rho_f * c_f * u)\n", "rate = h * a / ((1 - eps) * 1000 * 2000)\n", "\n", "t = np.linspace(0, 3000, 151)\n", "exact = [validation.schumann(xi, rate * (ti - Z / u))[0] for ti in t]\n", "fig, ax = plt.subplots(figsize=(6, 4))\n", "for n_z in (50, 200, 400):\n", " solution = pybamm.Simulation(\n", " sensible, parameter_values=pv, var_pts={base.r: 5, base.z: n_z}\n", " ).solve([0, t[-1]], t_interp=t)\n", " theta = (solution[\"Outlet temperature [K]\"].entries - 305.15) / 38\n", " ax.plot(t / 60, theta, label=f\"{n_z} cells\")\n", "ax.plot(t / 60, exact, \"k--\", label=\"Schumann\")\n", "ax.set_xlabel(\"Time [min]\")\n", "ax.set_ylabel(\"Scaled outlet temperature\")\n", "ax.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "d3ac1e2b", "metadata": {}, "source": [ "With the default, second-order scheme the model is within 1 % of Schumann's solution on 50\n", "cells. On finer meshes about 0.3 % is left: the difference between a capsule that conducts\n", "heat and Schumann's lumped one. (With first-order upwinding, numerical diffusion leaves\n", "1.3 % on 400 cells: see the [advection scheme](advection_scheme.ipynb) notebook.)" ] }, { "cell_type": "markdown", "id": "c25ab236", "metadata": {}, "source": [ "## Spherical Stefan: a melting capsule\n", "\n", "One capsule, solid at its melting temperature, with its surface held 1 K above it (Stefan\n", "number 0.011), against the quasi-steady solution, exact as the Stefan number tends to 0." ] }, { "cell_type": "code", "execution_count": 3, "id": "4d4965fc", "metadata": { "execution": { "iopub.execute_input": "2026-11-01T11:17:15.007442Z", "iopub.status.busy": "2026-11-01T11:17:15.007301Z", "iopub.status.idle": "2026-11-01T11:17:15.796604Z", "shell.execute_reply": "2026-11-01T11:17:15.795658Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhgAAAFzCAYAAAB8X3AUAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZilJREFUeJzt3Xd4FNXbxvHvbnoHAiGF0ENvoQTpCNJEpIs/BVFEpajYldcGCIKKWMCORrGigPQuiESFEAgSIAESQkiBdNKzKXvePyKrkQSSsMmkPJ/r4jI7O2fmmWFl78ycOUenlFIIIYQQQpiRXusChBBCCFH7SMAQQgghhNlJwBBCCCGE2UnAEEIIIYTZScAQQgghhNlJwBBCCCGE2UnAEEIIIYTZWWpdgBaMRiNxcXE4OTmh0+m0LkcIIYSoMZRSZGRk4OnpiV5f+nWKOhkw4uLi8Pb21roMIYQQosaKjo6mSZMmpb5fJwOGk5MTUHRynJ2dNa5GCCGEqDnS09Px9vY2fZeWpk4GjKu3RZydnSVgCCGEEBVwoy4G0slTCCGEEGYnAUMIIYQQZicBQwghhBBmVyf7YJSF0WgkLy9P6zJENWZlZYWFhYXWZQghRLUkAaMEeXl5REZGYjQatS5FVHP16tXD3d1dxlMRQoj/kIDxH0opLl26hIWFBd7e3tcdRETUXUopsrOzSUhIAMDDw0PjioQQonqpFgHj0qVLRERE0LlzZ1xcXMrUJjo6mvj4eNq0aWPWR00LCgrIzs7G09MTe3t7s21X1D52dnYAJCQk4ObmJrdLhBDiXzT99TwoKIhJkybRtWtXBgwYQHBw8A3b5ObmMnHiRNq2bcu0adNwd3dn5cqVZqupsLAQAGtra7NtU9ReV0Nofn6+xpUIIUT1omnAOHnyJFOmTOHw4cNlbrNw4UICAwOJiIggNDSU7777jscff7xc2ygLuacuykI+J0IIUTJNb5Hcf//9AMTExJS5jb+/P7Nnzzbd8x43bhydOnXC39+f3r17V0aZQgghhCinatEHo6zi4uKIj4+nR48exZb7+fld9/aKwWDAYDCYXqenp1dajeJa48aNo0mTJqxatarE10IIUR0ppcgrNJKTV0hWXiHZhoKi/+YVkG0oJCuvgJy8QvILjeQXKgqNisICA/q8TPR5megKsrDMz8bCaEBvNKDLz8FGr7BUBiyN+ZwOO0dGeho52dnkZGeTlZPL0B6taeHRAFDoUFcrQafU9Uq9Lrt2Q+l662SznJPyqFEBIyUlBQBXV9diy11dXU3vlWTp0qUsXLiwUmvT2rFjx3jjjTc4cOAA+fn59OrVi6VLl+Lr61tsvTVr1rB48WKioqLw8fFh6dKljBkzRqOqhRCichmNivTcfFKz87mSnceV7Hyu5OSRmpXPlZx/lqVm55GWnYcypGOZl4513hVs8tNxVJm4kIlD/hWs8lKxNKSjM2RiNGRRmJeLi1UBd7SxxIFc7MnluZ2ZXEw3kparSDOov/8LabmKjm56gh9xNNU2b3Um4SnFh0O4zeYkfXVWZj0Hf1ragASM67OyKjrpubm5xZbn5ORct1Pm/Pnzeeqpp0yvr84EV5u88sorTJ8+nffffx+9Xs///d//MXToUEJCQvDy8gJg586dPPjgg3z++eeMGTMGf39/JkyYwKFDh665KiSEENVRQaGR5Kw84tNziU83EJ+eS3JmXlFAyCkKClf+DhNp2QZ0OSnUy0ugvkrFyy6fRro0GurSCA8NxZiTjjEvm4LcXGzyDDjl5pFhUDSvr+edkbbw9/e8x9sZXM4s+QqCr7ueWW3/CQ3bzuUTkVryuil5VkQ5dadQb02BzhqfZsE41TPgYGdj+pPTqg2BjRsDoHQ6oKifl6Li/b0cWg+ocNubUaMCxtVxKWJjY4stj42NpWnTpqW2s7GxwcbGprLL09TWrVuLvV61ahVffvkle/bsMfV1Wb58OWPGjOG+++4D4KmnnuL777/n3Xff5euvvy5xu0op3n33XT766COio6Pp3LkzK1asoH///qb333nnHT788ENiYmJo2bIlTz31FDNnzixz7bt27eL//u//CAsLw9PTk4ceeoinn35aHvsUog4xGhWp2XlFoSEjl4R/BYj4dAMJGblcTsslKdOAtTGbBtlRuBYm0aqejka6KzTmCns2HiAzI4tcQx45hnwyDEbSDQqjgsHNLVg/3cG0v4Y7MkjOKSU0eFhS4OCOsq2Hzr4+Dg6/Q2Ymer0OFwdbXBztcXFyxNnJgXatW8DMBWDtCDaOvOS5mdwChUv9Bri4uFzz59/DKmx/opJPqsaqfcAIDw8nIyMDX19f7O3t6du3L5s3b2bq1KkAZGVlsXfvXhYsWFAp+1dKkZNfWCnbvhE7K4sKP6WQmppKYWGhaVwRpRSHDh3i9ddfL7be0KFD+emnn0rdziuvvMLHH3/M6tWrufXWWzlz5gxffPGFKWC89tprrF27lu+++45OnTpx5MgRJk+ejKOjI3ffffcN68zMzGT8+PG8+eabzJgxg4SEBD799FNCQkLo1q1bhY5dCFH95OYXEpOaQ3RKNhdTsolOySb2Sg6X03NJ+DtA5BcacSYbd10KHroU3Ejm/NE/yM9IIis9jcL0LAozDJxPL+SMEW5racG70/4JDa9eyCAxu+TQkGm0heYDwNENHBoxot8esvJ1uNSrh0u9BjjXb4BLAzdcGjbGq2kLLEePNrX9Y0wC9vb2ODg43PDf5PsfnmOeE1YLaBowEhISOHv2LImJiQCEhIRgaWlJ06ZNTVckli1bxqFDhzh58iQAixcvZtiwYcyfP58+ffqwcuVK3NzcePjhhyulxpz8Qjq8sqtStn0jpxeNwN66Yn9F8+bNw9vbmxEjRgCQkZFBVlYWjRo1Kraem5sbly9fLnEb2dnZvP3227zzzjuMHTsWgF69etGrVy+g6FbVm2++ycaNG/Hz8wNg0KBBPProo3z++edlChgJCQnk5ORw5513Ym9vT/Pmza8JQUKI6s9oVCRkGIhOzeZi8t8hIjXbFCji0w2AwpksXBL/wiH9AtaZsRjTE7DPSKVBeiaJ6bk0bahj6z3/DHLoFlByaNDrIE9vBy0GgWNjcHRjxVNRWDnWp76bFy5uTXBp3BQX9+a4NGhYNDDev8LBt6PeKPOxubm53dS5qas0DRhHjx5lyZIlAPTr14+1a9eydu1aZsyYwYwZMwDw8fEpNunYoEGD2L9/Px988AGBgYF07tyZr7/+GkdHxxL3URe9+OKLbN++nX379t1wNFKlVKmJPCwsjJycHAYOHFji+6GhoWRlZTFy5EjTtq7+admyZZlqbdGiBaNGjaJPnz7cfffd3HrrrQwdOtQ0SqYQovrINBRwMbl4cDD9NzWHvAIjFioft6xzuFwJwzo1koLUODrY5PHlAAea6BJx0uXg9kPpVxqs9Xqwqw/OXuDsyX3DIyiwsMfbuylNmrfEu3VHvNt0xaNFGywti3+FTR1RFWdBlJWmAWPUqFGMGjXquus8//zz1yzr168f/fr1q6yyirGzsuD0Im0+tXZW5e+DsGDBAt5//3127NhhutIA4OTkhIODg+lq0VWJiYk0/rtD0X+pvx+LKi2AXJ0M7sSJE3To0KHctV7d9rZt2/jtt9/Ys2cPL774Io888gh79uyp8DaFEBVjNCpir/xzG+Piv8JDdEo2KVlFv+xZqnzc86Lwscugue4yg3WXWbv1IJGXrhCVksf5guLb9Wmgp/2gf34J7OZlS1KuBd6N6+Pt2bgoPLTwwdunE83a9wCf9qZ1l99bJYcuKkG174OhNZ1OV+HbFFVt4cKFvP3222zfvt3UR+IqnU7HLbfcwv79+3n88cdNy/ft20ffvn1L3F779u2xs7Pjt99+o127dqW+v2PHjpsKAzqdjkGDBjFo0CAWL15Mr169+OKLL1i+fHmFtymEKN3VIHEuIYOz8Zmcjc/gXHwm4QmZ/+pzpnDOvEjDKyewSw3HIzUap9RkElMziErJw8pFT8Dcf0LDqsRMwhKKfumw0EPzBjb4eDWgdTMv2rdrA/feD/Wagos3uxfIPE91Qc345hQ3tGTJEpYvX8727dsZMKDkR5KeeeYZxowZw5o1axgzZgxffvklwcHBfPTRRyWub29vz5NPPslLL72Eh4cHgwcP5syZM3z++ed89NFH2Nvb89xzz7Fo0SI8PT0ZPXo0qampbN++neTkZF566aUb1n3o0CH8/f2ZM2cObdu2JSwsjJiYGFq1anVT50MI8U+QCE8oChFn4zM5l5BBeEIm2Xn/BIn6WVG4phynQ9ZFJnV2pKN1HC1UDP2/jGN/vLHEbV9MUxgbdUTfsBU0aMlrDbMwOjTGp1tfmnfshZXM51TnScCoJV577TUMBsM1/SVeffVV0xM2I0eO5PPPP+e1117joYcewsfHh/Xr1193DIzFixdTr149nnjiCeLi4ujatWuxKwsLFizAzc2NxYsXc//99+Pp6ckdd9zByy+/XKa6e/XqxYkTJ5g+fTpnz56lcePGPPzwwzzyyCPlPwlC1GHJmQbOXM4g7HJG0X/jMzgXn1EsSHiSTIP4QFpdCiE/OYaUpCQiE7KIyi4KEbaW8GMXJywKi26LdmlsQWaBJT5eDfBp4U3r1m3w6dgNn259aNahJ/p/Pf4/clhVH7Go7nRK3cT4ozVUeno6Li4upKWlXTPVe25uLpGRkbRo0QJbW1uNKhQ1hXxeRFXLzS/kXHwmYZfTOXM5gzPxRaEiMeOf6RAc8lNxTwrELvEkutSLvDnKhRYqFluVw8Qfs9kQWryThA5o2dCGjs0a8uVL91G/pS80aotybY3OWm5niOKu9x36b3IFQwghqqksQwEnYtI4Hn2FkNgrhF3K4EJyFkbTr4WKJrpEmsYfwvPiMXISY7h4KYWwpDxO/+tXxw9vccTWRQ96SwZ28MRgV0jHNi3p2LkrHXsOoP0tw7B3rn/N/mWuYHEzJGAIIUQ1UGhUnEvI4PjFKxyPLvpzNj7DFCZ0qhDPzFA6JxzBKukMz/a3p6tVNLbGLJ4Ny2X5n3nFtufmaEG3Fo3o3K4V+jH3Q+cB0KAl8yysmFf1hyfqIAkYQgihgfj0XIJNYSKVkJg0sv7uL6HDiFfOOTrEBGCReIaUy3GciUvjj6x/LkvMaWqPbQtLsLBmkG8TzqscfDt3xLf3AHwHj8ajZYcKjwQshDlIwBBCiEqWnVdAyN+3Oq7+uZR2ddJGRWNDDD4JvzOocRYjGibSqiCcL0KuMGd78YkdLfTQ3sMR37bNcBl5DwwZCw3bcIeFFXdU/WEJcV0SMIQQwowKjYrwhEz+ir5C8L9udRT+fa/DxphNs6Q/8Lh8lLz480TGJBCUmI9Rwe132NLBuejxzlua2dPV25JeHVrSs2cvfPsPp/OA27FzkFGLRc0gAUMIIW5CQkauqd9E8MUrhMSmkWm4+pSGohmXuFMfQT/7SKwuneB/q8M4W3Dtdpo2sKGwaT8Y+wB4dsO3YVuOfyj/RIuaSz69QghRRkajIvRyOkciUzh28QpHo1KJvZJjel+vCmiZfoyGl/4g79IZwi5cpn8bHUuG2oIRUp0VOQXgbKunV2s3evt2pvfA2/AbPgn3pmWbv0eImkIChhBCXEdMajYB55IICE/ij4hk03wcALYY8FNh2JzeTvLFc/wVmcQvGcVHvnS1tQTv3tCkF/Wb9OTc3Q1p2aUveovyzzUkRE0iAUMIIf4lLTufP88ncfBcEr+HJ3EhOdv0npMxg54pv9FJhTOrQxZumWHojPk0+iWT5JyiPhaWeujZ0pV+PTvTf8go+t1xD3g0MW2jdZUfkRDakIAhqoVFixbRsGFD5syZU+JrISqLoaCQo1GpBPwdKEJi00xjT9iqbDolHaRR/GFiz5/l8PkrnDTAYWcdrz7hWPQYqEsTHh9ji0WDpvQfPo5ewyZh7ygdMYWQgFGLXLhwgU8//ZRjx47x3HPPMWTIkGvWMRqN/PDDD2zfvh2dTsc999zDqFGjiq1z7tw5Vq5cSVRUFD4+Pjz55JN4eXlVau3Hjh2jSZMmpb4Wwlyu9qO4etvjyIUUcvOLbmtYk09PXTh3OJ/j993b+CkwjjM5xWdTcLGzwLdtU7JGvoJju1uhXlNekfEmhLiGBIxawt/fn9dee40ZM2awa9cupk6des06RqOR8ePH89dff/HCCy/QqFEjPv74YwoKChgzZgxQFC78/Py4/fbbufvuu/n222/x8/MjODgYNze3qj4sIcwiJjWb38OLbnv8ux+FDiOts07gGbefjIunWDOyEEcLA+TBsdxcUnMUjjZ6BnZqwuABA7h17L34DhiOhfSfEOKGJGDUEqNHj2b69Ono9fpSZzL99NNP2b17N6dOnaJly6Ie6xMnTiQ9Pd20zmuvvYaPjw/ffPMNOp2OSZMm0bp1a95++23eeOONUvd/4sQJVq9eTXR0NJ07d+bpp5/GxcXF9P7x48dZvXo1MTExtGzZktmzZ+Pj41Pm47t06RKrVq0iLCwMT09PHnjgAbp3717m9qJuudqPIiA8iYBzxftRuObH0yd+N/roo4ScjWLv5XzTe4fb2zO0sye0GMicLm34n2Mbetx6B5ZWVlochhA1mgSMG1EK8rNvvF5lsLKHMl56LcvVBX9/f+68805TuLjq37Ph7dq1i8cff9w0xLCVlRVjxoxh586dpQaMzZs3c9ddd/HAAw/wv//9jzNnznDPPfewbds2ALZv3860adOYN28eQ4YM4ciRI3Tv3p2AgAC6du16w7oLCgoYMGAA3bp149577yUhIYFHH32UL774gnbt2t2wvaj9DAWFHIu6QkB4IgHhyYTEXDH1o7Amn/4WZ5nkcoaEv/Yz78dIjv1rHAqdDro3b8DwgX40nTUPeo8AnY422hyKELWGBIwbyc+G1z212ff/xYG1g9k2FxISwsSJE3n//ffZv38/bm5uTJw4keHDhwOQnZ1NQkIC3t7exdo1adKECxculLhNo9HI7Nmzefzxx3nzzTdNy69eFVFKMXv2bN566y1mzJgBwIQJE0hLS2Px4sX89NNPN6z7woULRERE8Oeff9KoUSMAHnroIXJzc2/QUtRmmYYCfj2TwM6Tl/n1TOK/BrcCz4KLdE7cQ0b4Ue5umc6drYyQDcedCsktAK961gzv3ZHho8cxdPJMGrlr9P+4ELWYBIw6QilFbm4ub7/9NqNHj2bq1KmcPn2asWPHsnz5cubOnUteXtF9aTs7u2Jt7e3tTe/9V1hYGHFxcdx9993Fll+9KnLmzBkuXrzImjVr2LBhA0oplFJERUWhlCppk9fw9vbGy8uLhx9+mDlz5tC3b18cHBxwcDBf+BI1Q0pWHntD49l18jIHw5PIKyjqnKnDiG/eX3hG7+ZsyF8cDE/nz6J5w9BnWHFnl2bQeihdxt3KqZlNaN+jr0wEJkQlk4BxI1b2RVcStNq3meh0OlxcXGjRogVffPEFUNT/IjMzk2XLljF37lwcHR2xtLQkJSWlWNvk5GTq169f4nYzMzMBSn0/LS0NgPvuuw9Pz+K/Jf43yJTGxsaGw4cPs2rVKubPn28KRp988kmx2zuidrqUlsPuU/HsPHmZwAsppjk9rMnnLpcIxtodZ+bSH9mYUDwEN3W1ZfSA7kz83/0w6UHQ69EDHar+EISokyRg3IhOZ9bbFFry9fW9pq9GixYtSEpKQimFpaUlHTt2JDg4uNg6wcHBpfaVaN68OQCnT5+mRYsWpb7v5OTEyJEjK1y7l5cXS5cuZenSpcTExNC/f3/effddXnnllQpvU1Rf5xMz2XUqnp2nLvNX9BXTcsf8ZHok7aRJ5ikWdU/CwpAJBrDVF2Chh/7tGjN6xG3cfu8cOnTvI1cphNCQBIw65MEHH+SJJ54gNjYWLy8vDAYDa9eupX///qZ/iKdPn85rr73G888/T4sWLTh27Bi7du3iq6++KnGbbm5ujBkzhgULFnDLLbfg6upKYWEhP/zwA/feey+NGzdm7NixvPLKK/Tp08c0tsW5c+cICQlhwoQJN6z7zJkznD59mvHjxwPg7u6Os7MzhYWFZjozojo4F5/B9pDLbA+5xJn4DNPyxjnn6Ri/g7hTQfx2JpFT+WBjAfPbOuHYwAPajuLr1R1o0nss9Ro21vAIhBD/JgGjljh+/DgvvPCC6fVbb73FN998w8iRI3niiScAuPfeezl69CgdOnSga9euRERE4OHhwbp160ztHnvsMYKCgujcuTMdO3YkJCSERx555Jo+Fv/2xRdfMGXKFFq2bEmnTp2Iiori0UcfNb3v7+/P9OnTadu2LV27diU1NRW9Xs8HH3xQpmNzdXXl66+/Zu7cubRp04aIiAg8PT157LHHynmWRHWilOLMv0JFeELm1XfoYhFFu9gt/H7wDw5FZhD4r+46TetbM26wL9l3vYZjl6Gg19NJkyMQQlyPTpW1p10tkp6ejouLC2lpadfcw8/NzSUyMpIWLVpga2urUYXll5yczJEjR65Z7u3tTceOHYsti42N5dy5c3h4eNCmTZsSLyNHRERw8eJFWrdufc1TJaU5d+4ccXFxdOzYkYYNG17z/sWLF03hwMfHB71eb3ovODgYW1tb2rdvX+JrgMuXL3P27FkaN25M27Zty1RTZaupnxetKKUIvZTB9pBLbD95ifOJWQBYUkC3rENMd7/AEMtgHHIu8e4hA0/uMgDQrakT44b1Z+y0OXQdOFpufQihoet9h/6bBIxaEjCENuTzcmNKKU7FpbM95BI7Tl4mMqkoVFhQSNfMP3C/sIugv05zJCaPz++0ZYavNVjZE1v/Fn4878j4GU/RvIOvxkchhLiqrAFDbpEIIcxOKUVIbBrbQy6z4+Qlov4eSVOHkW7ZQXhGbSXw2El+jvnnyQ+dDsL07eF/S6DlYLys7HhSqwMQQtw0CRhCCLM5n5jJhmOxbPorluiUnL+XKnpZnucR179on/orzd86z9/DV6DXweDOTZg0aTLjH3wad8/KnVRPCFF1JGAIIW5KWnY+W07Esf5YDMEXr5iWtzBG0+PSBvRxwbw5IA/SAD0MamFDrlV9/nfXJCbPegE3DwkVQtRGEjCEEOWmlOLQ+RS+C7zIrpOXySssuiThTCa3pW8l+sgeth2P59cC0AFP9nLDo9cY6DiBbc/2x8ZBBkgToraTgCGEKLOUrDzWHY3m+8BoU2dNPUZGGA9hDN7AtsNneSfDaFq/dWMH7hk3CsvHloNXMwBsNKlcCFHVJGAIIa5LKcWf55P5PjC62NWKTtaXecb9GP0y9/JVwEUe2ls0+VwDewvuHnEL9z36An63yiOlQtRVEjCEECUq6WqFszGN25J+5mLQPu5sls3gHtYATO7ekF0pdtx7/yPcPu1RrG3kOoUQdZ0EDCGESUlXKywpoG/Gr9if2caOwAg++vsWSFyiBQ//byx0+x8ubUby0wIJFUKIf0jAEGa1ZcsWnJycGDx4sNallMmuXbuwsLDgtttu07oUTZV0taKVLpauF74h8HAg318wmNZt6GjJtNEDmPHUq+A3SKuShRDVnAQMYVaff/45TZo0qXDA2Lx5My4uLgwaVDVfXF9//TW2trZ1MmCUdLXCmnwmWR9ljtNBWmYFc2dwNr9fKECvg5E9WzHj4TmMue9RrK2ttS5fCFHNScAQZjVmzBjq1atX4faffvoprVu3rrKAUReVdLWiSWEMvWK+59jhw8wdY0kLvR50ep66awB+GU24/8mFNGnWQuPKhRA1iQSMWuTSpUvs2LEDJycn+vfvz4kTJ4pd/v/2229p37493bt3N7XZsWMHVlZWpnX27dvHiRMngKJZTHv16kW7du2u2dfBgwcJCwvD09OTIUOGYGdnBxRN3+7k5HTDWktqv3v3biIjI8nIyODdd98FiqaPr1+/PgaDgT179hATE0PLli0ZMmQIlpb/fHzLWvd/z9G/HT58mLNnzzJt2rRrao2Njb3ujLLVnVKKw5EpfHf4Ijv/vlphQSG3Zu/DKmQjPwde5PfsommJVp9yZcmLT0P3aQx2acJgbUsXQtRQEjBqicDAQG677TZ8fX1p2rQpzzzzDHZ2dvTv398UHt566y2mTp1aLGB89dVXODo6mtZJTk7mwoULQNEX7pw5c1i4cCFPPfWUqc3EiRMJCgpixIgRJCQk8OKLL7J9+3Y8PT3LdIuktPZJSUnk5OSQnp5uqiE/P5/z588zYsQIGjVqRKdOnfjss88oLCxk3759NGjQoMx1X+8cAeh0Ou6//34GDx5cbAbZOXPmMG7cuAr9vWjtSnYe64/F8t3hKCL+nrm0nkpjbPKPhP6xl69Cs7g622GTBnbMvGc8Dz6zBJo116xmIUQtoeqgtLQ0Bai0tLRr3svJyVGnT59WOTk5xZZnZmaW+qc862ZnZ5dp3fLy8/NT999/v+n1kSNHlE6nUw8++KBpWdeuXdVbb71VrN2UKVOKrfNfv//+u7KxsVGJiYlKKaWioqIUoM6fP29aJyIiQsXFxSmllBo7dqyaO3duqdu7UfvRo0erefPmFWszYMAA9fTTT5teG41GNWLECPXYY4+VuW6lynaOunTpohYtWmR6ffjwYaXT6VRERESJ+ynt86K16JQs9eqmk6rtS9tVs+e3qmbPb1UTXv5AHXl3irryoqtyskZB0Z9h3Vuqjd98qvLz87UuWwhRA1zvO/Tf5ApGGTk6Opb63u233862bdtMr93c3MjOzi5x3UGDBvHrr7+aXjdv3pykpKRr1lNKXbOsNJcvXyYwMJCVK1ealvXs2ZPevXuXeRv/du7cOY4cOUJSUhJGo5GCggJCQ0MZMGAATk5O2NjYsGXLFh555BFsbGxo2bJlmbdd3vaxsbEcPHiQPn36sGrVKpRSKKWws7Pj999/L3PdZT1HM2fO5N133+Wll15Cp9PxxRdfMGjQoHIdo5bOXM7gkwMRbP4rjgKjwpICJqoDOJzeyMLuSehSdWAJjw7xJrNeB+b83xu069xV67KFELWQBIxaIDY2FgAvr+KTRjVp0qTc21q4cCFvv/02Q4cOxdPTEysrK3Q6HSkpKQDUr1+fH3/8kVdffZX58+fTp08fpk2bxn333XfNiI0RERFs2bLF9Hr48OF06NChzO0BYmJiAEhKSsJg+OdRyWbNmuHr61vmust6jqZOncpzzz3Hr7/+yi233MIPP/xQLJRUV0ejUvjo1wj2hiYA4KQyGZ+5jr9+3cE7p4tugwx3c6L/qEnQ+xFef7VX0fzoQghRSSRglFFmZmap71lYWBR7nZCQUOq6er2+2Our/QZuRuPGjQFITEws9gWakJCAi4uL6bWlpSWFhYXF2mZlZZmuzmRnZ7No0SJ27tzJsGHDgKLjfu+994pdUbnzzju58847SUxMZOvWrTz66KMUFhYyY8aMYtvOzs4udnxXz2FZ20NRh02AadOmldqvoyx1l/Uc1a9fnwkTJvDFF1+YQsmkSZNK3K/WlFIEhCex8pdwAi8UBSlP4yUGxa5h1/4/ee1SgWndkX5tsJv+PgwcoVW5Qog6RgJGGTk4OGi+bmm8vLxo27Yt3333Hd26dQPg/Pnz/PHHH/j4+JjW8/b25tSpU6bXV65c4c8//zR1YMzIyMBoNNKoUSPTOl999VWxfV0NT25ubjRq1IgHHniAr776itOnT19TV+fOnU1Pg5S1vbOzc7HbS61bt6ZDhw6sWLGCgQMHmgJaYWEhZ8+epX379mWqu6znCOChhx7i9ttv5+zZs9x9992mJ2SqC6UUv55N5P1fzpmmR+9pGcH9xs3M+nAfyzKKQpWtlY77xgxm3qsr6NClm3YFCyHqpGoRMKKjo4mPj6dNmzY4O994GmelFJGRkaSmpuLt7Y2bm1sVVFl96XQ6VqxYwdixY4mPj6dZs2Z8++23pt/ar5o1axZjxozB3t4eDw8P1q9fj62tren9xo0bc9ttt3H33Xczbdo0zp8/z/bt24s9DpqWlsaoUaMYMGAAbdu2JSwsjCNHjvD222+XqdYbte/fvz+vvvoqTZs2xdHRkenTp/P1118zcuRI+vXrx6hRo0hNTWXv3r3Mnj2b9u3bl6nusp4jgMGDB+Pl5XVNnw2tKaX4JTSB9/ed40RMGjqMDNcH8VKD/TTN/ItCo8LWEtxdbHh0xv94ZP6bNPxX6BJCiKqkacDIzc3l3nvvZceOHTRr1oyoqCjeeOMNHnvssVLbnDhxgilTppCamoqXlxdhYWGMHDmSb775ptr9plmVbr/9do4cOcL69etxdHRk/fr1LF++vNg6I0aMICAggB07duDg4MDatWs5cuQINv+amGrr1q2sWbOG8PBwOnfuzOuvv86KFSto1aoVAD4+Phw/fpx169Zx9uxZfH19WbJkiem2w40G2rpR+1mzZuHm5kZwcDBJSUnk5+fTvXt3zp49y08//UR4eDjNmzdn/fr1tGnTpsx1l3aOTp48iZWV1TV1Dh8+HBsbG/z8/Mr/l2FmRqNi9+l43v/lHKcvpWNBIbelbyf3yE/sCU9l5VxHsLbGoutktm2+nea9hhf7OxVCCC3oVHkeVzCz+fPn88033xAYGIiHhwcbN25k/PjxHDp0qNQnIPr164eTkxNbt27F0tKSixcv0qVLF1566SWeeeaZMu03PT0dFxcX0tLSrrlikpubS2RkJC1atCj2231NNHXqVGxtbVm9erXWpdQoBQUFtGjRgueee+66YRcq9/NiNCp2nLzMyn3nCLucgQ15DElZz6U/N7P5ZIZp/IqfX72LcU+9C84eZt2/EEKU5Hrfof+m6RUMf39/Zs+ejYdH0T+M48aNo1OnTvj7+5caMBITE4uN4ti0aVOaNGlCYmJildUtaq/Vq1ezY8cOdDpdiZ1Oq0KhUbH1RByr9oVzLiETR7IZk/A9Z3/fycdnc0zrjenbgRcWr6DvrdJxUwhR/WgWMOLi4oiPj6dHjx7Flvv5+REcHFxqu0WLFvHMM8/QvHlzmjVrxu7du8nOzmbOnDmltjEYDMUecUxPT7/5A6gBRo4cWeLlf1G6mJgYOnXqxFtvvWWWDrjlUVBoZMuJOFbuC+d8YhbOZPGc7W76p26hq39R51i9DqYM68X8pavo3F372zdCCFEazQLG1fEJrj6GeJWrq6vpvZIMGTIEPz8/5s+fT5MmTTh//jwvvPACTZs2LbXN0qVLWbhwoXkKr0GmTp2qdQk1zoIFC6p8n/mFRjYGx/LB/nAuJGfjpDK5P2c98xv9jk1hJtSHOzvXo1Grrjz/+gf4tO9Y5TUKIUR5aRYwrv5mnZubW2x5Tk5OqVNBK6UYNWoUzZo1IyYmBmtra6Kjo/Hz86OgoICXXnqpxHbz588vNidFenp6sbkmhNBCXoGRDcdi+ODXcKJTcnBUWYy7/BVBB/ayLCaPhx53xKtVRxj0HD+/dCd6y2rx0JcQQpSJZv9ieXt7o9frTYMZXRUbG1vq1Yi4uDiOHTvG4sWLTSHE29ubO++8k82bN5caMGxsbKRXvag2DAWF/BQUw0e/RhB7JQcnspkQv4bA/bt5LyoPAGtLHX96z2LS7DdAr0d/g20KIUR1o1nAsLe3p2/fvmzevNl0KT8rK4u9e/cWu0wdHh5ORkYGvr6+NGjQAL1ebxo++qro6OhigyyZg4YP14gaxGg0lnnd3PxCfgi8yMcHznM5PRdbDExO/Y7gX7byTkRRHyFrSx0PTx7JC298gpdcZRNC1GCaXnNdvHgxw4YNM81JsXLlStzc3Hj44YdN6yxbtoxDhw5x8uRJ7OzseOihh/i///s/CgsLadmyJbt372bnzp3FJhu7GVfnsEhMTKRRo0Ylzo8hhFKKvLw8EhMT0ev1pd7Wu7ru1hOXWLo9lLi0XCwpYI5DAPfmrcPniwvkFIClXsdDk4fz4lufSbAQQtQKmgaMQYMGsX//fj744AMCAwPp3LkzX3/9dbGZS318fMjLyzO9/uCDD+jduze7du1i/fr1NGvWjD/++INbbrnFLDVZWFjQpEkTYmJizDJPiKjd7O3tadq06TVzzFx1Oi6dBVtOERiZgg4jkzjAq/V34JQTAxbw+MCGXHbowKvvfE6LVq2ruHohhKg8mg60pZWyDBJSWFhIfn5+FVcmahILCwssLS1LvMqVmpXH23vO8N3hixiVYlD+H9ge/pwvDyWwf7o9fdt5wqDnUL73obOS/kFCiJqjRgy0VZ1ZWFhcM0uqEDdSUGjk28MXWbHnLGk5+bTPD6XdyVV8/msUGX9fiPs5uxd9520DawfkBpwQoraSgCGEmfwZkcyCzac4E5+BuzGR4eHv8e2ev9iZWXSR0Le1O2+sWMWwMRM1rlQIISqfBAwhblJCRi6vbwtl4/E4HMjhJdut/Pz99yw/W3TJooWbE0tef50pD8wpta+GEELUNhIwhKigQqPim0NRLN91hmyDgXst9jPf7mccC1KxaKfn92hLXn7mMR59camMwyKEqHMkYAhRAcejr/DSxhBOxqbTIy+IBodX0cUjC8dOVuDamvtef5Uxq/vQ0MzjswghRE0hAUOIcriSncebu87wfeBF3AsTmBC2nG/3hZKQpdjhbMHYx9/Ctt/DWFhY0VDrYoUQQkMSMIQoA6UU647GsHRHGFlZmdxx6Qt+37OTdy4VANC2SQPeXfUxtgMna1ypEEJUDxIwhLiBsMvpvLzxJEcupNA/5wAFv33KquPpADjbWbLghad4dP5i0wR+QgghJGAIUaosQwHv/XKOzwMiaaIu8bXNVxRePsqtx7MBmDF+CEs/+g63xo01rlQIIaofCRhC/IdSip0nL7Nwy2mS0zO5L28d/+eyAyuVD63seXWaHyNmvkyfgUO0LlUIIaotCRhC/MuFpCxe3XyKA2cT8c0Lpv0f77DieAozZjvi7TsERq9ggWsrrcsUQohqTwKGEBRNpf7xgQg+/DUCu/wrTD6/gu92HSUuo2gUzp+tJ/H4tC9AZtcVQogykYAh6rz9ZxJYsPkUUclZDM/aTuwef5afyQWglbsLH6/257bR4zWuUgghahYJGKLOSsnK46WNIWwPuYwnSdwTuoQPdp0j3QBWFjqenzOd/3vjQ+zs7LQuVQghahwJGKJO+iM8iSfWHichI5eplvt42fp7Xk1OId0AvTs05fPvN9GxSzetyxRCiBpLAoaoU/ILjbyz5ywfHYjAU13mcz5jqGUoGOHVe/rRasIQZsybj4WFhdalCiFEjSYBQ9QZ0SnZPPZ9MH9Fp3BH2k8Eb/+BFwuNDH7YFYthr2LX+xEe0kuwEEIIc5CAIeqETcdjeennkzgbYpj012t8sDeKnAJwtLEgZOBquvUZo3WJQghRq0jAELVabn4hr246xdqgi4zI3MKZbf4sv5APwNBe7Vn9w1aat2ypcZVCCFH7SMAQtVZkUhazvzlK/KUYJoUt4bNdp0k3gL21BctfX8Csp15EJ+NaCCFEpZCAIWql7SGXeG7dCXrlH+FLm08ZfSSWdAP06dScNet30rpNW61LFEKIWk0ChqhV8gqMLN0Ryg+/h/Gi5bdMtf4FgG8e7MTGwsE8u+htLC3lYy+EEJVN/qUVtUbslRwe/e4YOReCGBCwkHDrLBhsC30epf2Ql2lvZat1iUIIUWdIwBC1wq9nEnjih2D6Xf6OQ5u+ZdelQvQ6mPrKF7Qecq/W5QkhRJ0jAUPUaIVGxbt7z7Jm/3FGhr3Gh9tPkZkHDRytWf3ZZxIuhBBCIxUKGAUFBXz11Vf8/vvvpKSkXPP+xo0bb7YuIW4oMcPAvB+CyTxzgA77lvBmcAYAA7u15tvN+2ji7a1xhUIIUXdVKGDMmzePb7/9ltGjR9OkSRNz1yTEDR0+n8zj3x1lbNZPbPjen6Nxheh08PK8h3j5rQ+lI6cQQmisQv8K//jjj+zdu5eePXuaux4hrstoVHx68Dyf7jrKcosPGGJzHLceVrycbc23369lyEgZkVMIIaqDCgUMvV5Pu3btzF2LENeVlp3P0z8dJ/rk77yX8w4DGl5BWdry4CvvMbn1OFzq1dO6RCGEEH/TV6TRsGHD2LBhg7lrEaJUJ2KuMHrlQTjiT+r3zzL1m2iSrLzQPbgbXc/7JVwIIUQ1U6ErGLa2tsyYMYNNmzbRunXra4ZbXrZsmVmKE0IpxTeHL7JsywnuvLicL34+SEKWor6DFef6vkNDj65alyiEEKIEFQoYFy5cYPDgwaSlpXH06FFz1yQEAFmGAuZvCOHQ8RBGHn+Rt/bEUKigayt3Nuz8jZatfbQuUQghRCkqFDD27t1r7jqEKCb0UjpzvzuG0+VAmu9eyIoT2QDcc8dgPlu7DXt7e40rFEIIcT0V6oMhRGX6IzyJiR/9QdfknbQIeJWfTmRjoYd3Fr/IN5v3SbgQQogaoMKDBURGRvLee+8RGhqKUooOHTowb948WrRoYc76RB2z53Q8j34XxDx+YI71ZlIHWnEoyYEl73/J0FF3aF2eEEKIMqrQFYx9+/bRvn17Dhw4QIsWLWjVqhUHDhygffv27Nu3z9w1ijpi0/FYnvrmd+bFL2CO5WYA6g9/hj/D4iVcCCFEDaNTSqnyNurVqxejR49mwYIFxZYvWLCAbdu2ceTIEXPVVynS09NxcXEhLS0NZ2dnrcsRwDeHovho4z56HHqBlQcSWHWHI3OXfA5d7tK6NCGEEP9S1u/QCgUMGxsb4uPjqfefsQdSU1Nxd3fHYDCUu+CqJAGjevno1wh2bPsZ612vsj6kqDPnc7Om8cZHazSuTAghxH+V9Tu0QrdI6tevz9mzZ69ZfubMGRo0aFCRTYo6SCnFmzvD+GPTZ6T+OJ/1IdlY6nV8+t4bEi6EEKKGq1Anz6lTp3LXXXexePFi/Pz8ADh8+DAvvvgiU6dONWuBonYyGhULtpwiedfbHPlxDWeSjTjbWbJ+3QZuu13mExFCiJquQgHj9ddfR6/XM3PmTNPtEBsbGx5//HEWL15s1gJF7VNoVDz703G8jizjs69/4nKmwruhI9v2/Ebnbr5alyeEEMIMKtQH46rc3FwiIiLQ6XS0bNkSW1tbc9ZWaaQPhnYKCo08szaIgaELmWARwIo/DXx9vgHbfj2Cp5eX1uUJIYS4gUrt5FnTScDQRl6Bkee/C+CO0OcZanMao84C/Z3vk9t+Uo0Jp0IIUdeV9Tu0zLdIHn30UQBWrVpl+rk0q1atKutmRR1hKCjk+TX7sV57H0//Fc8vMxviOv1b8LkNiRZCCFH7lDlgxMTElPizOURHRxMfH0+bNm3KfEWhsLCQ0NBQ7O3tadmypVnrEeaVm1/I8/47Kfh6Ou8FJAPwk+NMZvncpnFlQgghKkuZA8bGjRtNPy9btox27dqVuF5YWFiZd56bm8u9997Ljh07aNasGVFRUbzxxhs89thj1223YcMG5s6di729Pfb29ri7u/P999/TsGHDMu9bVI2cvELmr95Eqv8DfB2UDsDrLz7FI89LZ2AhhKjNKvQUSfv27Smt68b13vuvhQsXEhgYSEREBB4eHmzcuJHx48fj5+dH7969S2zz22+/MXnyZD777DNmzJgBwK+//kp8fLwEjGomy1DAS5+sJXb1w6wLyUangw+XL2bWUy9qXZoQQohKVqFOnjqdrsQQkZmZibu7O5mZmWXajru7O7Nnz+bVV181LevcuTP9+vXj448/LrHNkCFDsLKyYteuXeUt20Q6eVa+jNx8Fny4hnOrH2VLaC4Wevj6sw/534zZWpcmhBDiJpi9kyfAE088UeLPAEajkePHj9O9e/cybSsuLo74+Hh69OhRbLmfnx/BwcEltjEYDAQEBLBixQoyMjI4e/Ysnp6eeHh4XHdfBoOh2PDl6enpZapRVExaTj5LP/FnZvSLDIsxYG2p48dv1zD2LhmETQgh6opyBYzw8PASfwawsrKiR48ezJs3r0zbSklJAcDV1bXYcldXV9N7/5WUlER+fj4nTpygbdu2uLu7Ex4eTp8+ffjuu++u2dZVS5cuZeHChWWqS9ycK9l5LPnYn1fTXsHRJY/NT/clueezDBs9VuvShBBCVKFyBYytW7cCcP/99/Pll1/e1I6trKyAoo6e/5aTk4O1tfV12+zZs4fjx4/j5uZGcnIyffr04ZlnnsHf37/EdvPnz+epp54yvU5PT8fb2/um6hfXSs40sGDlp4w99yqOTQvJ8upP9+k/gbW91qUJIYSoYhWa7OxmwwWAt7c3er2e2NjYYstjY2Np2rRpiW0aNmyIg4MDEyZMwM3NDSi64jF58mQOHjxY6r5sbGxwdnYu9keYV0JGLi+9+yknP3uGO9aksjG1DQ4SLoQQos6qUMAICAhg9uxrO+vNnj2bgICAMm3D3t6evn37snnzZtOyrKws9u7dy7Bhw0zLwsPDTX0y9Ho9w4YNuyaUxMTE0KhRo4ocijCD+PRcXn5/NWf8n+HXyDysLC1oOG6phAshhKjLVAX06tVLHTly5JrlgYGBqnfv3mXezq+//qqsrKzUCy+8oDZt2qRuu+021apVK5WRkWFa58EHH1QdO3Y0vT5x4oRycnJSr7zyitq1a5datGiRsrS0VBs2bCjzftPS0hSg0tLSytxGlCw2NVvNfG2lGtbaWgHKwcZCHdy3W+uyhBBCVJKyfodW6ArG1U6W/9W2bVtOnDhR5u0MGjSI/fv3ExUVxXvvvUfHjh0JCAjA0dHRtI6Pj0+xJ1M6d+7MH3/8QUxMDG+++SaRkZEcOHCA8ePHV+RQxE2IvZLDS6u+JO7rZ9gTnoedlZ5tW7fS/9ZhN24shBCiVqvQOBg+Pj6sWLGCMWPGFFu+adMmnn766WueMKluZByMm5dpKODx977nypqH+flULrZWerZu2sTQUXdoXZoQQohKVNbv0ApdwZg1axYPPvggn376KWFhYYSGhvLJJ58wc+ZMZs2aVeGiRc1QaFQs+nonC9JfxlZXgJWFjvU/rZVwIYQQwqRCQ4U/9dRTpKSk8MQTT5CTkwOAnZ0dTz75JE8//bRZCxTVz6otfzDr4tM0tUzlswc682SvN+k1QCYuE0II8Y8K3SK5Kjs7m9DQUHQ6He3atcPevmY8NSC3SCru5z9Pk/TBaB5slUyuvRf2s/aCi5fWZQkhhKgilTJU+H/Z29tfM9S3qL2OhccStGQ0r2+7wJ4uDny/Z5OECyGEECWqcMAoKCjg9OnTXLx4kYKCgmLvjRs37mbrEtXMpdQMdr40kmXbLwDQauAULNzaaFuUEEKIaqtCAePcuXOMGzeOM2fOUFhYiJWVFfn5+QA4ODiUeTZVUTMY8gv4/v8m8Nr6kxgVPHDXHbz+/mqtyxJCCFGNVegpkieeeIJBgwaRlZUFFM0ncvToUbp168aSJUvMWqDQ3levz2HJV7+QWwDD+nXn029/RqfTaV2WEEKIaqxCnTxdXV0JDQ3Fzc0NnU5HXl4eVlZWnDhxgvHjxxMREVEZtZqNdPIsu61r3ubxec8SeUXRxcebP4JDcXBw0LosIYQQGqnUcTBSUlJMk401bNiQy5cvA9CqVatr5gkRNVf40b2oXQuJzVA0aejEnoNHJFwIIYQokwoFjH/z8/Nj2bJlhIeHs2TJElq3bm2OuoTGUmPCcN3yAGN8dHz8kB+7f/0Tt8aNtS5LCCFEDVGhTp6PPPKI6edly5Zxxx138OGHH1K/fn1+/PFHsxUntFGQkUTKp+NppU/njL4VE9/YhbNzPa3LEkIIUYPc1EBbVymliI2NpXHjxlhZWZmjrkolfTCuo7CAd+/ryjvbQlk9xZMWz/1K61ZyVUoIIUSRSu2DceuttxZ7rdPpaNKkSY0IF+L6fn7tHl5Ye5qLaYo1ObdKuBBCCFEhFQoYR48elbEuaqG/fl7B3LfXYSiEnr6d+erLr7QuSQghRA1VoYAxatQovv/+e3PXIjSUEhbAI48/z6VMhXfj+uzedxC9/qb7AAshhKijKtTJ08nJiUceeYT169fToUMHrK2ti72/bNkysxQnqobKiGfe1Ns5HFOAo60l63b+Rv16LlqXJYQQogarUMC4cOECQ4YMoaCggBMnTpi7JlGVCvP5YPZgvjmagQ5488PP8evWSeuqhBBC1HBlDhhBQUH07NkTgL1791ZaQaJqXV73LPd6xbKltTWNBkxj9gP3aV2SEEKIWqDMN9l79epl+lkG06od0o5twD3Un/p2OiY88xZfffaJ1iUJIYSoJcocMOrVq2caBry6zzUibiwv8Tx7334IpRQ/Wo9n6v1zsLCw0LosIYQQtUSZb5GMGzeOzp0706JFCwDT7ZKSBAUF3XxlovIU5LFg+q0s3ZHCxB4NWbZlBQ42FeqOI4QQQpSozN8qq1evZuzYsYSHh3Ps2DEmTZpUmXWJSrT+tXt5c+dFAJoNf4TWHg00rkgIIURtU6GhwidNmsS6desqo54qUZeHCj+7dw23Tbyf6HRFvz5+HPz9EDqdTuuyhBBC1BBl/Q41y1wkNU1dDRg5iZHcNbADW8Ny8WzoxImz0bjWl/EuhBBClF2lzkUiaiBjIe/MGs7WsFysLHR8sXaLhAshhBCVRgJGHXFy7QJe3xIOwJx5TzJiyCCNKxJCCFGbScCoA4yXQmh75kM23m3PwFt8WfHWW1qXJIQQopaTgFHbFRhI/24GVhSgb3ELP+36XSYxE0IIUenK/JjqzJkzy7zR1atXV6gYYX5HPn2MhhfDKKhXj7gByxjibKd1SUIIIeqAMgeM3Nxc08+ZmZls2rQJLy8vfH19AQgODiY2NpaxY8eav0pRIdlh+7n35c+5lGlk/PS7WT24h9YlCSGEqCPKHDC++eYb088PPfQQ8+bNY/ny5VhaFm2ioKCAZ555huzsbPNXKcrPkMn/PTKZcylGGjjZMuXh57G2lFsjQgghqkaFxofetm0bJ0+eNIULAEtLS15++WW6dOlituJExR18/yHe/y0ZgP4zFzC6p0xQJ4QQoupU6FfajIwMoqKirll+8eJFMjIybroocXNyzv3Gg8t+RAGtu/Zg5YuPaV2SEEKIOqZCAeOuu+5i8uTJ/Pjjj0RGRnL+/Hl+/PFHJk+ezJQpU8xdoyiPAgML507hXIqR+o42TH/lI5q62mtdlRBCiDqmQrdIVq1axfz587nvvvswGAwA2NraMmvWLJYuXWrWAkX5/PXlsyzfexmAZmOf5Kk7fDWuSAghRF10U3OR5OTkcP78eQBatmyJnV3NeASy1s5FkhBG9sp+vLY/k/UprXnj212M922idVVCCCFqkbJ+h1boCsZVdnZ2dOzY8WY2IcxFKdj+DPYWBQy5tS8nGy1kbFcvrasSQghRR1X4ucUffviBYcOG0apVK9OyxYsXk5CQYJbCRPlkBH6HMfI3cpUVCwum8+IdHdHrZRp2IYQQ2qhQwPjss894/PHHGTBggOkWCUCDBg14/fXXzVacKKO8LB6aNYc+n2fxclx/2rfvxC0tXbWuSgghRB1WoYCxYsUK1q1bxyuvvFJs+e23387atWvNUpgou1/em8va4+kExRnZYBzE/FHttC5JCCFEHVehPhiRkZH06tULAJ3un8vw9evXJyUlxTyViTLJuxTK3DeLRllt392PSWOH07KRo8ZVCSGEqOsqdAXDy8uLU6dOAcUDxrZt22jdWkaMrEorn/kfZ5IKqe9ghfG2Z3l8qI/WJQkhhBAVCxizZ8/mgQceYO/eveh0Ok6dOsXbb7/N7NmzmTt3rrlrFKWI/2Mti9b/BUDjgf9j3ihfGjhYa1yVEEIIUcFbJE8//TQZGRmMHTuWwsJCOnXqhJ2dHc899xxz5swxd42iJAUGXnx6DukGaObRgAYDpjG9b3OtqxJCCCGACl7BOHPmDAsXLiQ5OZljx44RFBREUlISCxYsICwsrNzbi46OJigoiPT09HK1S0tLIyAggMjIyHLvs6bL++1d/rp4BQD9kCd5/vYO2FpZaFuUEEII8bcKBYz27dsDRcOD+/r60qNHD+zt7Yu9Vxa5ublMnDiRtm3bMm3aNNzd3Vm5cmWZ2iqluPfeexk0aBDvvfde+Q+iJku/hPWf73LoQQcm3PM/uvS5lTFdPLSuSgghhDC5qZE8/yszMxMHB4cyr79w4UICAwOJiIjAw8ODjRs3Mn78ePz8/Ojdu/d1277zzjum2zN1zm9vQn42x2nDUe97WDe6fbHOtkIIIYTWyhUwnnjiiRJ/BjAajRw/fpzu3buXeXv+/v7Mnj0bD4+i377HjRtHp06d8Pf3v27AOHr0KCtWrCAoKIiRI0eW5xBqvMKEs3zy0Sfc382SZcb/MaqTBz2bN9C6LCGEEKKYcgWM8PDwEn8GsLKyokePHsybN69M24qLiyM+Pp4ePXoUW+7n50dwcHCp7TIyMrj77rv54IMPcHd3L9O+DAaDadZXoNx9PaqTNQsfZO72HJYfsUc/sy2/yKBaQgghqqFyBYytW7cCcP/99/Pll1/e1I6vDsjl6lp8SGtXV9frDtY1e/Zshg4dytixY8u8r6VLl7Jw4cKKFVqN5EYe4ZVvfgfAutMI7u3bkmauZb8lJYQQQlSVCnXyvNlwAUVXPKCoo+e/5eTkYG1d8lgOmzdvZuvWrUyYMIGAgAACAgLIysoiLi6OgICAUvc1f/580tLSTH+io6Nvun4tfPzyTGLSFa5Otlj1ncFjQ2RQMyGEENVTma9gPProo2Xe6KpVq264jre3N3q9ntjY2GLLY2Njadq0aYltCgoK6NSpE4sWLTItu3TpEtnZ2cTFxXHgwAEsLK59VNPGxgYbG5sy118dZZzczes/hwBQv99dzBvRgXr2MqiWEEKI6qnMASMmJsasO7a3t6dv375s3ryZqVOnApCVlcXevXtZsGCBab3w8HAyMjLw9fVlwoQJTJgwodh2unXrxuDBg3n33XfNWl+1ohQrX5lLYraicX0H3AdN5b4+zbWuSgghhChVmQPGxo0bzb7zxYsXM2zYMObPn0+fPn1YuXIlbm5uPPzww6Z1li1bxqFDhzh58qTZ919TpJ/YwfIdRZ1qHfpN47lRHbG2rNDdLSGEEKJKaPotNWjQIPbv309UVBTvvfceHTt2JCAgAEfHf2YD9fHxue6jr76+vrRs2bIqytVM5i9vM7SlJZ4NnWncZwK3d5ZBtYQQQlRvOqWUqkjDH374gc8//5zz588TEREBFF2RePjhh3FzczNrkeaWnp6Oi4sLaWlpODs7a13O9V0IgC9Hk48lfTOWMWPMEGYPbqV1VUIIIeqosn6HVugKxmeffcbjjz/OgAEDOH/+vGl5gwYNeP311yuySVGaA28A8EPBYNLsmjCll7fGBQkhhBA3VqGAsWLFCtatW8crr7xSbPntt9/O2rVrzVKYgOwzB3jio92cTdXxUcGdjOniKdOxCyGEqBEqFDAiIyPp1asXQLE5MOrXr3/dQbJE+Xyx7GneO5zHoG8KiFWu3O0nVy+EEELUDBUKGF5eXpw6dQooHjC2bdtG69Yy+JM55F8OY/nPxwBw6TmGlo0c6dmsvsZVCSGEEGVToYAxe/ZsHnjgAfbu3YtOp+PUqVO8/fbbzJ49m7lz55q7xjrpx+VPEpWmaOBgTXanyUzq0URmTBVCCFFjVGi69qeffpqMjAzGjh1rmjLdzs6O5557jjlz5pi7xjpHZSay/Nu9AHj2uJUsaxsmdm+icVVCCCFE2VUoYOh0OhYuXMj8+fMJDQ3FaDTSvn177O3tzV1fnbTvk//j+OUC7Kz0XOl2H0PbNMLdxVbrsoQQQogyq1DAuMrW1hZfX19z1SIA8nNYsfpbADr49iDJzoW7ekrnTiGEEDVLuQLGyJEjy7Tezp07K1SMAPXXD3RxLeSwvQUJ3R6kgb0VQ9tX74HLhBBCiP8qV8DYtWsXzZo1o3///pVVT92mFLrDn7D0Nlu6THqE+ZebMK6bFzaW184QK4QQQlRn5QoYzz77LGvWrCEwMJAZM2Ywffp0PDxkXgyzufgnJIairOxZnugHwOSe0rlTCCFEzVOux1TffPNNoqOjWbZsGb/99hvNmzfnzjvvZNOmTRQUFFRWjXXG9k8WsCeigLMNh5NcaEdHT2c6erpoXZYQQghRbuUeB8PKyooJEyawfft2IiIi6NWrF0888QRNmshv2jdDZSTw9Gf7GP5NNo8dKLqwNLmHnFMhhBA1001N156WlsaVK1fIysrCxUV+074Zv325kLCkQhxtLAhrPBJLvY4xXT21LksIIYSokHIHjPT0dD799FNuueUWevfuTUpKChs2bODMmTOVUV/dYCzkky+LHk0d0Lcneht7BrVphKujjcaFCSGEEBVTrk6e06dPZ926dXTu3JkHH3yQu+++Gycnp8qqrc5ICtrI+r/SAEjscC8A43y9tCxJCCGEuCnlChhr1qyhWbNm1KtXj/Xr17N+/foS15NxMMrn6w/fJK8QurRoRKJjSxxtLBnWobHWZQkhhBAVVu4rGMK8VHYKn+84CkDbgWMJBEZ1csfWSsa+EEIIUXOVK2B8+eWXlVRG3ZV48CvyCxW2VjpOuY8AYHx3uT0ihBCiZrupp0jEzXO7uJWwuQ5sevdJsrDDw8WWW1q4al2WEEIIcVMkYGgp8QzEBqGzsGKXXdE8L3d280Sv12lcmBBCCHFzJGBoKO6XT8ktUOS3vI3N4fkATPCVwbWEEELUfBIwtKIUjy9dTePlGSz5qx75hYr2Hs60dZfHfoUQQtR85erkKcznSuhvbDmVTl4hHLbsAQrG+8rInUIIIWoHuYKhkZ8/X05eIbTzcuG00ROdDsZ2k6dHhBBC1A4SMLSgFGu3/QpAj74D0el09GvVkMbOttrWJYQQQpiJBAwNJJ/8hb1nMwGIaj4WgPEyNLgQQohaRAKGBjZ+8Q6FCjo1rU+03h1bKz0jOrlrXZYQQghhNhIwqprRyPqdvwHQoVd/AIZ3cMfRRvrbCiGEqD3kW62qxQbx9hBFHy8ndrUcD8CdXeXpESGEELWLXMGoamd20L6RBc8+OJ4YvRs6Hfi1bKB1VUIIIYRZScCoaud2AxDm3BeAdu7OONtaaVmREEIIYXZyi6QK5SWe56FPjnB7a0vOd2oHZNKreX2tyxJCCCHMTgJGFTr404es+SufXZE6eg4omtCsZ3O5PSKEEKL2kVskVWjrls0AjOzXhdOXMwDkCoYQQohaSQJGVcnPYduR8wB0GjQWowKvenZ4uNhpXJgQQghhfhIwqkjEbz9yLrkQSz3ktxkFgF8LuT0ihBCidpKAUUV2rv8agH4dvAhJygegp9weEUIIUUtJwKgiuw4GATD8tiEEX7wCQC/p4CmEEKKWkoBRBVRaLClpRZObtR48hZz8QlzsrGjdyFHjyoQQQojKIY+pVgHdhYMEzHDgskNHNtdrC5yhZ7P66PU6rUsTQgghKoVcwagKkUWTm7l3u42gqCsA9JIOnkIIIWoxCRiVTSkKzv1a9GOLQQRFpQLS/0IIIUTtJgGjkuVeOkOjV0K57esc/spvRkpWHrZWejp7uWhdmhBCCFFpqkXAiI6OJigoiPT09DKtX1hYSGhoKOfOnaOgoKCSq7s5gTu/50ounErWcSqtaFKzbt71sLasFqdeCCGEqBSafsvl5uYyceJE2rZty7Rp03B3d2flypXXbbNkyRK8vLyYOHEiw4cPp3nz5mzdurWKKi6/3/bvAWBg11am2yN+cntECCFELadpwFi4cCGBgYFEREQQGhrKd999x+OPP87hw4dLXL+wsJCcnBxOnz7N6dOniYyMZObMmUyZMoXLly9XcfVl81vQaQAGDhhAYGQKIB08hRBC1H6aBgx/f39mzpyJh4cHAOPGjaNTp074+/uXuL6FhQWLFy+mQYN/vqBnz55NdnY2x44dq5Kay6MwK5VD59MAaDdoHLFXcrDQ6+jeVEbwFEIIUbtpNg5GXFwc8fHx9OjRo9hyPz8/goODy7ydI0eOANCqVatS1zEYDBgMBtPrsvb1uFkn9/1ERh442ehJa9ARCKGjpzMONjL8iBBCiNpNsysYKSlFtwtcXV2LLXd1dTW9dyNJSUk89thj3HXXXbRt27bU9ZYuXYqLi4vpj7e3d8ULL4c/920H4JZ2HhyKLOp/0VtujwghhKgDNAsYVlZFT1Tk5uYWW56Tk4O1tfUN26elpTFy5Ejc3d1ZvXr1ddedP38+aWlppj/R0dEVL7wcPNUlRvtYMnxQXwLOJQHQ36dRlexbCCGE0JJm1+q9vb3R6/XExsYWWx4bG0vTpk2v2zY9PZ3hw4djYWHBzp07cXJyuu76NjY22NjY3HTN5WI0cqdbLHfeY0/MlKdZ9VUS1hZ6eYJECCFEnaDZFQx7e3v69u3L5s2bTcuysrLYu3cvw4YNMy0LDw8v1ifjargA2L17Ny4u1XTAqpQIyMsASzt+TW0IQI9m9bGzttC4MCGEEKLyadrbcPHixQwbNoz58+fTp08fVq5ciZubGw8//LBpnWXLlnHo0CFOnjxJfn4+o0aNIjIyki+++IKQkBDTej4+PjRu3FiLwyhR0umD5GUY8WzXkd/Ci/pf9PdpqHFVQgghRNXQNGAMGjSI/fv388EHHxAYGEjnzp35+uuvcXT8ZxpzHx8f8vLyAMjOzkan0+Hj48PSpUuLbeuFF17gjjvuqNL6r8f/6+94bnUmD9x2mSN9kwHo31oChhBCiLpB8+cl+/XrR79+/Up9//nnnzf97OLiQkBAQFWUddOOhYQBUN+rJRmGAtycbGT+ESGEEHWGTIhRGZQi+HwCAHmNuwAwtH1j9HqdllUJIYQQVUYCRiXIuhTO2cR8AIIsOgIwrIObliUJIYQQVUoCRiU4+ft2FODmZMklowv21hb0bSX9L4QQQtQdEjAqQcjRPwFo6VE05sXwDo2xtZLHU4UQQtQdEjAqwYmTRTOoOrgW3RYZ372JluUIIYQQVU7zp0hqo1Et9VjeYs3Zpl1wc7KhXyvXGzcSQgghahEJGOamFKM8khk1wpZhhoGM7+6FpYVcKBJCCFG3yDefuWVcAkM6BUpPtM6D+/o017oiIYQQosrJFQwziw/9k/CLBdg08GRo9yZ41bPTuiQhhBCiyknAMLOf1q/nMf9sujVPZu2zPlqXI4QQQmhCbpGY2fGQoidIPDw8aNP4+tPICyGEELWVBAwzuxBzCYA2beTqhRBCiLpLAoYZKaW4lJQGQNdu3bQtRgghhNCQBAwzik7OJDqlaGr5Hn6lzxArhBBC1HYSMMwo+PgxMvIUAG269tG4GiGEEEI7EjDMKPyvQwC4O1li6+CocTVCCCGEduQxVTOqb5HFylG2JNs21boUIYQQQlMSMMzI1TKLGX7WnPXopXUpQgghhKbkFokZWWXFA6B38tC4EiGEEEJbEjDMKDL8HH9GF2CwaaB1KUIIIYSmJGCYiVKKj7efpO8X2RyLvKJ1OUIIIYSmJGCYSYahgKQMAwCt2rTXuBohhBBCWxIwzORKZg6JmYUANG3dUeNqhBBCCG1JwDCT2PNnKFSgA9xbddK6HCGEEEJTEjDMJO78KQBc7fVY2thpXI0QQgihLQkYZhJ/MRIAV0drjSsRQgghtCcBw0wS42MBqO9oo3ElQgghhPZkJE8zaeNVj/dG2mBo2ELrUoQQQgjNScAwk2YNrPlfbxuC6rfTuhQhhBBCc3KLxEx0uVcAKLB20bYQIYQQohqQKxhmci78PIUpBWQ0sNK6FCGEEEJzcgXDTD7acpQB/tn8GRavdSlCCCGE5iRgmElWdtEw4Q71XDWuRAghhNCeBAwzycrNA8CpvpvGlQghhBDak4BhJtmGAgAc6jfSuBIhhBBCexIwzCTDYATAuV59jSsRQgghtCcBwwyUUmTnKQCcXBpoXI0QQgihPQkYZpCTk4P6+2cXuYIhhBBCyDgY5qDX61k+3Ja8AoWzSz2tyxFCCCE0JwHDDGytLHm6T9EsqlF2DhpXI4QQQmhPbpGYQ6HB9KONrZ2GhQghhBDVg1zBMIMrKUmExRTgYKWjia291uUIIYQQmpMrGGZwIiSEPp9nc8f32dja2GhdjhBCCKE5CRhmkJWVCYClXoeNpZxSIYQQolp8G0ZHRxMUFER6enqltqksOX8HDCsLPXq9TuNqhBBCCO1pGjByc3OZOHEibdu2Zdq0abi7u7Ny5Uqzt6lsuYZcAAkXQgghxN807eS5cOFCAgMDiYiIwMPDg40bNzJ+/Hj8/Pzo3bu32dpUNkNODiABQwghhLhK0ysY/v7+zJw5Ew8PDwDGjRtHp06d8Pf3N2ubylaQV3QFQ6eTgCGEEEKAhlcw4uLiiI+Pp0ePHsWW+/n5ERwcbLY2AAaDAYPhn7EqzN1vIz+/aKp2dNWiS4sQQgihOc2+EVNSUgBwdXUtttzV1dX0njnaACxduhQXFxfTH29v75sp/RpejRuwYJANd/WUqdqFEEII0PAKhpWVFVDUafPfcnJysLa2NlsbgPnz5/PUU0+ZXqenp5s1ZHTzG0hjw0Pg5G62bQohhBA1mWYBw9vbG71eT2xsbLHlsbGxNG3a1GxtAGxsbLCpxAGwmrXvQbP2PW68ohBCCFFHaHaLxN7enr59+7J582bTsqysLPbu3cuwYcNMy8LDw039K8raRgghhBDa0imllFY7P3DgAMOGDePpp5+mT58+rFy5ksjISI4fP46joyMAM2fO5NChQ5w8ebLMbW4kPT0dFxcX0tLScHZ2rrTjE0IIIWqbsn6HavrYw6BBg9i/fz9RUVG89957dOzYkYCAgGJBwcfHh+7du5erjRBCCCG0pekVDK3IFQwhhBCiYmrEFQwhhBBC1E4SMIQQQghhdhIwhBBCCGF2EjCEEEIIYXYSMIQQQghhdhIwhBBCCGF2mg0VrqWrT+aae1ZVIYQQora7+t15o1Eu6mTAyMjIADD7rKpCCCFEXZGRkYGLi0up79fJgbaMRiNxcXE4OTmh0+nMss2rM7RGR0fL4F0VIOfv5sk5vDly/m6enMObU1POn1KKjIwMPD090etL72lRJ69g6PV6mjRpUinbdnZ2rtYfjOpOzt/Nk3N4c+T83Tw5hzenJpy/6125uEo6eQohhBDC7CRgCCGEEMLsJGCYiY2NDa+++io2NjZal1Ijyfm7eXIOb46cv5sn5/Dm1LbzVyc7eQohhBCicskVDCGEEEKYnQQMIYQQQpidBAwhhBBCmJ0EjHIIDw/n6NGj5OTkVGqb2iozM5OgoCAuXLhQ5jaXL1/m+PHjMqz7386fP8/Ro0fJysoqV7vExEQCAgKIiYmppMpqhuzsbI4ePUpERES52sXFxREcHIzBYKikymqOCxcuEBQURGZmZpnWLywsJDw8nGPHjpGcnFzJ1VV/hYWFBAYGcvr06TK3ycnJ4ejRo4SHh1diZZVAiRtKTExUffv2VS4uLqp169bKxcVFbdiwwextajN/f3/l4OCg2rZtqxwcHNSIESNURkZGqev/+uuvqnfv3qpx48aqa9euys7OTj322GOqsLCwCquuPq5cuaJuvfVW5eTkpNq0aaOcnJzUt99+W6a2eXl5qnfv3kqv16ulS5dWcqXV1w8//KCcnZ2Vj4+PcnJyUoMGDVIpKSnXbRMfH69GjBihnJycVM+ePVXz5s3Vxo0bq6ji6iUjI0ONGDGi2P/H/v7+123z559/qlatWilPT0/l6+ur7Ozs1LRp01ReXl7VFF2N5OTkqEWLFqlmzZopZ2dnNXr06DK127BhQ7Hvkb59+6rExMRKrtY8JGCUwaRJk1T37t1VZmamUkqpt956S9nZ2amYmBiztqmtTp06pSwsLNSaNWuUUkolJSWpVq1aqdmzZ5faZvXq1SowMND0+sSJE8rR0VG99957lV5vdXT//ferDh06qCtXriillProo4+UlZWVCg8Pv2HbZ599Vk2fPl01bty4zgaM8+fPK2tra/Xhhx8qpZRKS0tTHTp0UNOmTSu1TUFBgerZs6caPHiwSktLU0oplZqaqr755psqqbm6mTVrlvLx8VHJyclKqaJfGiwsLNTp06dLbdOhQwc1efJk0y8GYWFhysbGRn322WdVUnN1cvnyZfXiiy+qqKgoNWXKlDIFjOjoaGVnZ6dWrFihlCoKeV27dlWTJ0+u7HLNQgLGDaSkpCgLC4ti/6gYDAbl4uKi3nrrLbO1qc2ee+451bx582LLli9frhwdHcv1m8zIkSNrzP9Y5pSdna1sbW3Vxx9/bFpWWFioGjdurF599dXrtt25c6dq1aqVSk9Pr9MBY9GiRcrNza3YFbCPP/5Y2djYmH4J+K/169crQJ05c6aqyqy28vLylKOjo3r33XeLLW/atKl6/vnnS23XqFEjtXz58mLLvLy81Ouvv14pddYUZQ0Yb775pqpfv77Kz883Lfvyyy+VpaWlSk1NrcQKzUP6YNxASEgIhYWF9OjRw7TM2tqarl27EhwcbLY2tVlwcHCxcwHg5+dHZmZmme8pGgwGQkJCaN26dWWUWK2FhoaSm5tb7Bzq9Xp69Ohx3c9TfHw8M2bMYM2aNTg5OVVFqdVWcHAwvr6+xSZm8vPzw2AwlHov/JdffqFt27a0adOG06dPExoaSl5eXlWVXK2cO3eOzMzMa/4/7tmz53U/g0uWLGHlypWsWbOGPXv2MHfuXJydnZkxY0Zll1wrBAcH06VLFywt/5k2zM/Pj4KCAkJCQjSsrGzq5GRn5ZGSkgKAq6trseWurq6m98zRpjZLSUmhRYsWxZZdPTdlPR9PP/00OTk5zJ071+z1VXfX+zxFRkaW2EYpxbRp03jooYfo27dvpddY3aWkpODl5VVs2Y0+g3FxcTg7OzNw4ECSkpLIzc0lKyuLTz75hHHjxlV2ydXK9T6DFy9eLLXdyJEj2bBhA88//zzu7u5ERUWxcOFCGjduXKn11hYpKSklnvOr71V3cgXjBqysrADIzc0ttjwnJwdra2uztanNrKysSjwXQJnOx5IlS/D392fDhg3XfEnUBRX5PH322WecPn2awYMHExAQQEBAAPn5+URFRXH48OFKr7m6qchn0MrKiiNHjvDAAw9w+vRpzp8/z6xZs5g6dSoJCQmVXnN1UpHPYEFBAUOHDqVx48ZER0cTHBzM4cOHeeWVV1i1alWl11wb3Oy/nVqTgHEDzZo1AyA2NrbY8tjYWJo2bWq2NrVZs2bNSjwXwA3Px7Jly1iyZAlbtmxh0KBBlVZjdVaRz5NOp6N58+a89NJLvPDCC7zwwgtkZGSwY8cOXnvttUqvubqpyGewefPm2Nracv/995uWzZo1i6ysLI4dO1ZptVZHFfkMnjlzhnPnzvHII4+YLvH7+PgwfPhwNm/eXLkF1xI3829ndSAB4wY6dOiAp6dnsf8hzp8/T0hICMOGDTMtCwkJITQ0tFxt6ophw4YREBBAamqqadmmTZvo3Lmz6VJpeno6AQEBxZ6tf/PNN1m0aBFbtmxhyJAhVV53ddG8eXNat25d7PN06dIlAgMDi32eTp8+zcmTJwF46KGHTFcurv5p0KABs2bNYuvWrVV+DFobNmwYQUFBXLp0ybRs06ZNtGjRglatWgGQlZVFQEAAaWlpAIwYMQKDwUBiYqKpzdVxRBo1alSF1WvP3d2dTp06FfsMJicn8/vvvxf7DIaFhZn6Blw9R/8deyU6OrrOnb+yysnJISAggCtXrgBFn9sTJ04QFRVlWmfTpk14eXnRvn17jaosB617mdYEX331lbKyslLLly9X69evV926dVP9+/cv1iN90KBBauzYseVqU1fk5uaqTp06qf79+6uff/5ZLViwQFlYWKgtW7aY1jl48KAC1JEjR5RSSn3wwQcKUAsXLlQHDx40/Tlx4oRWh6GpdevWKUtLS7VkyRK1YcMG5efnp7p3716sd/no0aPV0KFDS91GXX6KpKCgQPXq1Uv16tVLbdiwQS1dulRZWlqqtWvXmtYJDg5WgNq/f79p2e2336769eunNm3apNauXas6dOighg0bpoxGowZHoa0tW7YoCwsLtXDhQvXzzz+rfv36qc6dOyuDwWBaZ+LEiapfv36m11OmTFGenp5q9erVaufOnWr27NnKwsJCHTx4UItD0NyhQ4fUwYMH1dChQ1Xfvn3VwYMH1R9//GF6PzQ0VAFqx44dSqmip8X69eunfH191fr169Vbb72lLC0t1VdffaXVIZSLdPIsg/vuu4969erx1VdfsXPnTsaMGcOzzz5brEd6586di/XUL0ubusLGxoYDBw6wbNkyVq5ciaurK7t27WLo0KGmdVxcXOjXr5/pHEZFRdGvXz92797N7t27Tet16dKFDz/8sMqPQWsTJ05k69atrF69mn379jFkyBCef/75Yr3LO3bseN2nHHr37o23t3dVlFvtWFhYsGfPHt58800++OAD6tWrx+bNmxk1apRpHUdHR/r164eLi4tp2fr163nvvff44IMPsLe3Z+bMmcyZMwedTqfFYWjqjjvuYNeuXXzyySccOHCAvn378sILLxTrC9C+fXvc3NxMr7/++ms+//xzdu3axZUrV2jZsiVBQUF069ZNgyPQ3muvvWa6OqHT6XjhhRewt7c3/Rtnb29Pv379qF+/PlD0tNiOHTt46623+Oijj3B2dmb9+vXceeedWh1Cuch07UIIIYQwu7r367QQQgghKp0EDCGEEEKYnQQMIYQQQpidBAwhhBBCmJ0EDCGEEEKYnQQMIYQQQpidBAwhhBBCmJ0EDCGq2OXLl/nxxx+1LuOGjh49yk8//cSRI0fKtP5/j0vL4yzLvqt7fULUdDLQlhBmkpyczJ49e667To8ePYiIiOCOO+6goKCgiiq7lsFg4Oeff2bChAklzso4d+5cNm3aRN++fbntttt4+OGHb7jNnTt3Fjuu/76uSv/dd1xcHL///juTJ08udZ3KouW+hdCSDBUuhJmkpqayceNG0+ugoCBSUlIYPny4aVmDBg3w8PBgypQpGlT4j3379jF//nzuvvvua94zGo18/vnnbNmypdZMznfs2DGmTZtW7Eu+qv4etNy3EFqSgCGEmbRu3ZoffvjB9HrWrFkEBQUVWwZFl8fHjh1reh0TE0NgYCDjxo0jJCSEqKgounTpQvPmzSksLOTw4cOkpKTQq1cv0+yz/3bp0iWCgoJwcnKie/fuODs737DWrVu3cscdd1yzPD4+nm3btmEwGDhy5AjJycnceuutZGZmEhYWxujRo03rJiQksG/fPqZMmVKhuTmq6riTk5M5ePAgRqPR9HfRtm1bPDw8Svx7GD9+PCEhIVy8eJGOHTvSokWLa7YZGBhIYmIiHTt2xNHRkb1795Z4Hsq776r8DAhR6bSda02I2uuRRx5RPXr0uGb5jh07lIWFhen1Tz/9pOzt7ZWfn5/q37+/GjhwoLKwsFArV65UvXr1UgMHDlT9+vVTzs7OKjAwsNi2Fi1apFxcXNSIESPUgAEDVKNGjdSePXtuWFuzZs3Uzp07r1keEhKiJkyYoAA1fPhwNWXKFBUSEqI++ugj1apVq2Lr7t+/XwGmGV3/e1z/ff1flXnc/973uXPn1IABA5Rer1dTpkxRU6ZMUV999VWpfw+33Xab6tOnjxo2bJiytrZWX3zxhWmdwsJCNX78eOXi4qJGjhypmjVrpkaOHFnsPPxbefddlZ8BISqbBAwhKkl5AgagPvroI9OyBx98UAHK39/ftGzKlClq/PjxptebN29Wbm5uKioqyrRs9erVyt3dXeXk5JRa14kTJ5SDg4PKzc0t8f3U1FQFqCNHjpiWVVbAqKzj/u++t2zZomxsbIrtv7S/h1WrVpmWLV++XLm5uZlef/fdd8rBwUGFh4crpZTKzMxUvr6+pQaM8u67qj4DQlQFeYpEiGrAwsKCBx980PS6T58+2Nract999xVbdvbsWdNrf39/OnbsSGBgID/99BM//vgjlpaWXL58mdDQ0FL3tXXrVoYNG4aNjU3lHEw5VOVxl7Wef3doHTx4MAkJCaYpttetW8f48eNp1aoVAA4ODsydO/em9vnvfVencyHEzZI+GEJUA46OjlhZWZle29jY4OLigl6vL7YsNzfX9PrChQsYDAbWrVtXbFtTpkwp1u6/tm7dyowZM8xYfcVV5XFXtB7AtP/o6GhGjBhRrE3z5s1vap/X27eW50KImyUBQ4gaytnZmaZNm7JmzZoyt0lOTiYwMJD169eXa196vR6j0Vhs2b+/6KpSRY7bXBo0aGC6mnFVampqlddxlZbnQogbkYgrRA01cuRINm3axKVLl4otj42NLbXN9u3b6datG+7u7uXal5eXF5cuXSIrK8u0bP/+/eUr2EwqctyOjo7k5+dfE5LKq1+/fmzbtq3Y+BWbNm26bhtz7bskFTkXQlQVuYIhRA01b948tm7dSu/evZkzZw4NGjTg2LFj/Pbbb5w+fbrENqU9nnojt956Kw0bNmTy5MlMmjSJ48ePX3NZvqpU5Lg7duyItbU1zz//PD169KBt27YV2vdjjz3GJ598wogRI5gyZQpBQUHs3LkToNRHdc2175JU5FwIUVXkCoYQlaRXr17X3K+HawdZ8vb2ZuLEicXWad68OePGjSu2rHXr1sXCgZ2dHfv37+f1118nMjKSY8eO0bNnT44dO1ZiPQUFBezevfuGAcPa2popU6bQoEED0zJ7e3sOHTpEt27dOHToEO3atWP79u3F7vX/97huNJhUZR73f/fdqFEj9u7di8FgYNOmTYSEhJTp76FevXpMmTIFOzs70+vAwED8/PwIDAykY8eOrFy5Ejs7OywsLEo8zoruuzI+A0JUJRkqXIg6Yv/+/dx7773ExsZWaGAsUSQlJaVY+Jo1axanTp3i4MGDGlYlRPUjt0iEqCNSUlJ44YUXJFzcpMmTJ9OnTx+aN2/OwYMH+emnn4oNES+EKCJXMIQQohxSUlJYvXo1Z86coUmTJtxzzz1m7VchRG0hAUMIIYQQZiedPIUQQghhdhIwhBBCCGF2EjCEEEIIYXYSMIQQQghhdhIwhBBCCGF2EjCEEEIIYXYSMIQQQghhdhIwhBBCCGF2EjCEEEIIYXb/DxDxvfVz5i0iAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "class SurfaceAt(BaseSubModel):\n", " \"\"\"Holds the PCM surface at a temperature, and sets the initial temperature\"\"\"\n", "\n", " def __init__(self, param, T_surface):\n", " self.T_surface = T_surface\n", " super().__init__(param, \"Surface\")\n", "\n", " def build(self):\n", " self.variables = {\n", " \"Phase-change material boundary temperature [K]\": pybamm.PrimaryBroadcast(\n", " self.T_surface, \"flow path\"\n", " ),\n", " INITIAL_TEMPERATURE: self.param.T_init,\n", " }\n", "\n", "\n", "capsule = EnthalpyStorage(param, param.pcm, geometry, inner_condition=\"temperature\")\n", "melting = pytes.ComposedModel(\n", " {\"storage\": capsule, \"surface\": SurfaceAt(param, param.T_m + 1)},\n", " param,\n", " name=\"Melting capsule\",\n", " defaults_from=base,\n", ")\n", "pv = base.default_parameter_values\n", "pv[\"PCM smoothing width [K]\"] = 0\n", "pv[\"Initial temperature [K]\"] = pv[\"PCM melting temperature [K]\"]\n", "tau = 861 * 213e3 * R**2 / (6 * 0.15 * 1)\n", "t = np.linspace(0, 1.05 * tau, 2001)\n", "fig, ax = plt.subplots(figsize=(6, 4))\n", "for n_r in (20, 160):\n", " solution = pybamm.Simulation(\n", " melting, parameter_values=pv, var_pts={base.r: n_r, base.z: 1}\n", " ).solve([0, t[-1]], t_interp=t)\n", " ax.plot(t / tau, solution[\"Average liquid fraction\"].entries, label=f\"{n_r} cells\")\n", "melted = np.linspace(0, 1, 101)\n", "ax.plot(validation.sphere_melting_time(melted), melted, \"k--\", label=\"quasi-steady\")\n", "ax.set_xlabel(\"Time / full melting time\")\n", "ax.set_ylabel(\"Melted fraction\")\n", "ax.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "a40ee977", "metadata": {}, "source": [ "Both meshes follow the quasi-steady solution. Early on, a coarse mesh deviates slightly:\n", "the melting front crosses the large outer cells of the sphere one at a time (the enthalpy\n", "method with a sharp melting point). Refined, the model is within 1 % of the quasi-steady\n", "solution from half melted onwards." ] } ], "metadata": { "kernelspec": { "display_name": "pytes (3.14.4)", "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.14.4" } }, "nbformat": 4, "nbformat_minor": 5 }