{ "cells": [ { "cell_type": "markdown", "id": "b3d34331-5842-4e5b-9b78-d985fcf67860", "metadata": {}, "source": [ "# Regioned Flasher\n", "\n", "The regioned flasher breaks up the temperature-density plane into multiple region. In each region, a K-D tree is constructed for points distributed within the region in T-D coordinates (because they are the coordinates of the EOS and require no iteration to obtain them). Thus distributing the points cannot fail (like iterative calculations in P-H coordinates might). To start, lets begin with a rectangular region in the supercritical region as a demonstration. The real regions are much more complex as they need to handle the complete fluid domain, deal with solid-liquid phase equilibria, etc.." ] }, { "cell_type": "code", "execution_count": 1, "id": "8a1f9562-40f5-4ffb-a0d4-d9478d3ad12d", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:47.391096Z", "iopub.status.busy": "2025-01-06T11:32:47.390565Z", "iopub.status.idle": "2025-01-06T11:32:47.658092Z", "shell.execute_reply": "2025-01-06T11:32:47.657745Z" } }, "outputs": [], "source": [ "import timeit, json\n", "\n", "import teqpflsh, teqp \n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "id": "dc129d47-71e5-4025-a813-20248654879f", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:47.659680Z", "iopub.status.busy": "2025-01-06T11:32:47.659577Z", "iopub.status.idle": "2025-01-06T11:32:48.280376Z", "shell.execute_reply": "2025-01-06T11:32:48.280068Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# of regions: 0\n", "# of regions: 1\n" ] } ], "source": [ "name = \"n-Propane\"\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", "\n", "rf = teqpflsh.RegionedFlasher(\n", " ideal_gas=json.dumps(jidealgas), \n", " resid=json.dumps(jresid), \n", " mole_fractions=np.array([1.0])\n", ")\n", "# To start off there are no regions in the regioned flasher\n", "print('# of regions:', len(rf.get_regions_rw()))\n", "\n", "# Now we make a region with rectangular shape in T, rho coordinates\n", "# As we will see, a rectangular shape with only the corners defined doesn't work so well when transformed into \n", "# other coordinates\n", "Tmin = 400 # K\n", "Tmax = 450 # K\n", "rhomin = 1e-6 # mol/m³\n", "rhomax = 6000 # mol/m³\n", "Tpoly = np.array([Tmin, Tmin, Tmax, Tmax, Tmin])\n", "rhopoly = np.array([rhomin, rhomax, rhomax, rhomin, rhomin])\n", "NT = 1000\n", "Nrho = 1000\n", "\n", "rf.add_region(T=Tpoly, rho=rhopoly, NT=NT, Nrho=Nrho)\n", "print('# of regions:', len(rf.get_regions_ro()))" ] }, { "cell_type": "code", "execution_count": 3, "id": "cdf82ace-ed39-4d0b-824c-1479fb9cd72e", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:48.281665Z", "iopub.status.busy": "2025-01-06T11:32:48.281551Z", "iopub.status.idle": "2025-01-06T11:32:49.136805Z", "shell.execute_reply": "2025-01-06T11:32:49.136540Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAG4CAYAAABBxj3JAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAANPhJREFUeJzt3Ql8FdW9wPF/9hAIW1gCj1VBkQdhrYhVyyuURR7Fltf6lCLlWVFKFbQfpFiKAio8nq9aCyqmT6EFpdW6IGWtyKIgIItisCiKkgohFmQJgazzPv9jc82F3EDI3Jk7M7/vp7c3d/7/OZmMyb1/zpxzJs6yLEsAAAB8LN7tAwAAAIg2Ch4AAOB7FDwAAMD3KHgAAIDvUfAAAADfo+ABAAC+R8EDAAB8j4IHAAD4XqLbBxArysvL5eDBg5Keni5xcXFuHw4AALgAun7yyZMnpWXLlhIfH7kfh4Lnn7TYad26tduHAQAALkJubq60atUqYpyC55+0Z6fihNWvX9/twwEAABfgxIkTpsOi4nM8Egqef6q4jKXFDgUPAADecr7hKAxaBgAAvkfBAwAAfI+CBwAA+B4FDwAA8D0KHgAA4HsUPAAAwPcoeAAAgO9R8AAAAN9j4cEoKiu3ZOv+o5J/8ow0S0+VK9s3loR47tMFAAiOshj5LIypHp7Zs2eblRInTpwY2tavXz+zrfLjjjvuCNvvwIEDMnToUElLS5NmzZrJpEmTpLS0VNy08v1Dcs1/r5Wbst+WCUt2mWd9rdsBAAiClTH0WRgzBc+2bdtk/vz5kpWVdU7stttuk0OHDoUec+bMCcXKyspMsVNcXCybNm2ShQsXyoIFC2TatGniFv0POW7RDjl0/IxkyhHpG59jnvOOnzHbKXoAAH63MsY+C2PiklZBQYGMHDlSsrOz5cEHHzwnrj03mZmZVe67evVq2bNnj/z1r3+V5s2bS/fu3WXmzJkyefJkeeCBByQ5OVmc7rqb/toesUTkhwlvyH8nZove3sOyvn7IiyJlL15M6xVdgAki8QnhofhEKU9Kk/jS0+abmK+tUimPS/zqOaWBSGp9kTMnJL7o+NevK1S1vW5TKc3sJomH3hUp/OLcwzk7Xt3rOKn+69NfSmmdRlVvu5BYdc8dBpjDTdz315rlXMw2B762Og6WuI9WRjVmx2vzGxulnNrsd6Ht2bF/LG73Uiya+9oRd6oNO3Ocaks/C+9fmhP6LJyd+DuJj7Ok3IqTX5T+RF4o+zfzWfmdzpmOXd6KsyzzEeyq0aNHS+PGjeXRRx81l7C0aHnsscdMTF/n5OSIHqYWPcOGDZNf/epXpghS2pOzdOlS2bVrV6i9/fv3yyWXXCI7duyQHj16VPk9i4qKzOPsu60eP368VjcP3fzxEdNlp1Xs5pQ7TbHjZfrbUd3PcHa8utfn+7q6bRcSq+5ZXUzOxWyL9tcFVorUiyuKaqy2r98p72ja7x3/ke05tdnvByXT5YWk+8/bnh37R8p1c7vySiya+9oRd6oNO3PcaKuqz0LN7Vv0W8mTDHn+tquk76UZUhv6+d2gQYPzfn673sOzZMkSU5joJa2q3HzzzdK2bVtp2bKlvPfee6bnZu/evfLSSy+ZeF5enunZqazitcYimTVrlkyffu5/wNrSQVlqfbz3ix11vp/h7Hh1r8/3dXXbLiR2vueLzbmYbdH8up589cEezVhtX+ub4NnnyK6c2ux3l/wp9AZd0+9R0/0j5bq1vZ9sN6+9EFsnvcxzNPa1I66caMPOHDfaUjfEbwh7v6nIvSF+ozxVfkPoM9P3PTy5ubnSu3dvWbNmTWjsztk9PGdbu3at9O/fX/bt2yeXXnqpjB07Vj777DNZtWpVKKewsFDq1q0ry5cvlyFDhrjSw/Nx0s2ScNZVJwDuKkvNkIQzRxzZP1KuW9tLLh1g/mmd9MnrMR8rufGPkrTkh1HZ1464cqINO3OcbktnZf342W3yh8QZcm3i387J21h6hYwq/VVweni2b98u+fn50rNnz7BByBs2bJC5c+eagiThrKqhT58+5rmi4NHLXFu3bg3LOXz4sHmONO5HpaSkmIfddLpdiwapIoW2Nw2glhJ0jJtD+0fKdWt7UtHJr69lxngsKTlRRHOisK8tceVEG3bmONzWtR2bms/ClFMlVeYlS4mJ62dmIGZpaU/N7t27zfibiof2+OgAZv367GJHVYzVadGihXnu27evaUMLpwraY6RVXufOncVpOvjq/mHOf18AFyAh2bn9I+W6tT2pjkhymjdiKlr72hF3qg07cxxuq+Kz8JSVWmVaoZVi4k6ux+NqD096erp06dIlbJteisrIyDDbP/74Y3nuuefk+uuvN9t0DM/dd98t1113XegS2MCBA01hM2rUKDNdXcftTJ06VcaPHx+VHpwLMbhLi4uchQUgqsqLnds/Uq5b20u058fySEyfC6O0rw1xp9qwM8eFtvSz8EizJJEvz03LapYkGV2+6rhwiuuDlqujU8p1urmO5zl16pQZYzNixAhT0FTQXqBly5bJuHHjTG+PFkw662vGjBmuHjvDd4AYVF7m3P6Rct3aXvr1mEVPxEqLo7dvbeNOtWFnjhttiUhGhH6HSNsDVfCsW7cu9LUWOOvXrz/vPjqLSwcoAwAAeKLgAYCoiUtwbv9IuW5tT0yJONg35mLVjUWq7b52xJ1qw84cN9qqSV6Qbi0BAFHHoGVvxBSDlu3NcaOtmuQ5gB4eAMHBoGWPxBi0bHuOG23VJM8BFDwAgoNBy96JMWjZ3hyXBi1f1H/HKOGSFgAA8D0KHgAA4Htc0gIQHMzS8kZMMUvL3hzFLC0ACAhmaXkjppilZW+OYpYWAAQEs7Q8EmOWlu05illaABAQzNLyToxZWvbmKGZpAQAA+Bs9PACCg0HL3ogpBi3bm6MYtAwAAcGgZW/EFIOW7c1RDFoGgIBg0LJHYgxatj1HMWgZAAKCQcveiTFo2d4cxaBlAAAAf6PgAQAAvsclLQDBwSwtb8QUs7TszVHM0gKAgGCWljdiilla9uYoZmkBQEAwS8sjMWZp2Z6jmKUFAAHBLC3vxJilZW+OYpYWAACAv1HwAAAA3+OSFoDgYJaWN2KKWVr25ihmaQFAQDBLyxsxxSwte3MUs7QAICCYpeWRGLO0bM9RzNICgIBglpZ3YszSsjdHMUsLAADA3yh4AACA73FJC0BwMEvLGzHFLC17cxSztAAgIJil5Y2YYpaWvTmKWVoAEBDM0vJIjFlatucoZmkBQEAwS8s7MWZp2ZujmKUFAADgbxQ8AADA97ikBSA4mKXljZhilpa9OYpZWgAQEMzS8kZMMUvL3hzFLC0ACAhmaXkkxiwt23MUs7QAICCYpeWdGLO07M1RzNICAADwNwoeAADgexQ8AADA9xjDAyA4mJbujZhiWrq9OYpp6QAQEExL90ZMMS3d3hzFtHQACAimpXskxrR023MU09IBICCYlu6dGNPS7c1RTEsHAADwNwoeAADgexQ8AADA9xjDAyA4mJbujZhiWrq9OYpp6QAQEExL90ZMMS3d3hzFtHQACAimpXskxrR023MU09IBICCYlu6dGNPS7c1RTEsHAADwNwoeAADgexQ8AADA9xjDAyA4mJbujZhiWrq9OYpp6bFj9uzZEhcXJxMnTjwnZlmWDBkyxMRfeeWVsNiBAwdk6NChkpaWJs2aNZNJkyZJaWmpg0cOwBOYlu6NmGJaur05imnpsWHbtm0yf/58ycrKqjL+2GOPmWLnbGVlZabYyczMlE2bNsmhQ4fklltukaSkJHn44YcdOHIAnsG0dI/EmJZue44K+LT0mOjhKSgokJEjR0p2drY0atTonPiuXbvkf//3f+WZZ545J7Z69WrZs2ePLFq0SLp37256gWbOnCnz5s2T4uJavrkB8JegT0uvbopwLMXMc5T2tSPuVBt25iimpbtv/PjxppdmwIAB58QKCwvl5ptvNgWM9uKcbfPmzdK1a1dp3rx5aNugQYPkxIkTkpOTE/F7FhUVmZzKDwAA4E+uX9JasmSJ7Nixw1zSqsrdd98tV199tQwfPrzKeF5eXlixoypeayySWbNmyfTp02t17AAAwBtcLXhyc3NlwoQJsmbNGklNTT0nvnTpUlm7dq3s3LnT9u89ZcoUueeee0KvtYendevWtn8fAADgPlcvaW3fvl3y8/OlZ8+ekpiYaB7r16+Xxx9/3HythdDHH38sDRs2DMXViBEjpF+/fuZrvcx1+PDhsHYrXld1CaxCSkqK1K9fP+wBAAD8ydUenv79+8vu3bvDto0ZM0Y6deokkydPliZNmsjtt98eFtfxOo8++qgMGzbMvO7bt6889NBDpnDSKelKCyUtYDp37uzgTwMg5rEOjzdiinV47M1RAV+Hx9WCJz09Xbp06RK2rW7dupKRkRHaXlUvTZs2baR9+/bm64EDB5rCZtSoUTJnzhwzbmfq1KlmILT24gBA2JtvaaEz+0fKdWt7deuexGLsYtZvYR2eyDmKdXi8LSEhQZYtWybjxo0zvT1aMI0ePVpmzJjh9qEBiDWsw+ORGOvw2J6jAr4OT8wVPOvWras2risun61t27ayfPnyKB4VAF8I+jo8kcRi7GLWb2Ednsg5inV4AAAA/I2CBwAA+B4FDwAA8D0KHgAA4HsxN2gZAKKGdXi8EVOsw2Nvjgr4Ojz08AAIjkhvvtHYP1KuW9t13ZPq1kSJpZiK1r52xJ1qw84cxTo8ABAQrMPjkRjr8Nieo1iHBwACgnV4vBNjHR57cxTr8AAAAPgbBQ8AAPA9Ch4AAOB7FDwAAMD3GLQMIDhYh8cbMcU6PPbmKNbhAYCAYB0eb8QU6/DYm6NYhwcAAoJ1eDwSYx0e23MU6/AAQECwDo93YqzDY2+OYh0eAAAAf6PgAQAAvkfBAwAAfI+CBwAA+B4FDwAA8D0KHgAA4HtMSwcQHKy07I2YYqVle3MUKy0DQECw0rI3YoqVlu3NUay0DAABwUrLHomx0rLtOYqVlgEgIFhp2TsxVlq2N0ex0jIAAIC/UfAAAADfo+ABAAC+R8EDAAB8j4IHAAD4HgUPAADwPaalAwgOVlr2Rkyx0rK9OYqVlgEgIFhp2RsxxUrL9uYoVloGgIBgpWWPxFhp2fYcxUrLABAQrLTsnRgrLdubo1hpGQAAwN8oeAAAgO9R8AAAAN+j4AEAAL5HwQMAAHyPggcAAPgeBQ8AAPA9Ch4AAOB7LDwIIDi4l5Y3Yop7admbo7iXFgAEBPfS8kZMcS8te3MU99ICgIDgXloeiXEvLdtzFPfSAoCA4F5a3olxLy17cxT30gIAAPA3Ch4AAOB7FDwAAMD3KHgAAIDvUfAAAADfo+ABAAC+R8EDAAB8j4IHAAD4XkwVPLNnz5a4uDiZOHFiaNvtt98ul156qdSpU0eaNm0qw4cPl7/97W9h+x04cECGDh0qaWlp0qxZM5k0aZKUlpa68BMAAIBYFDMFz7Zt22T+/PmSlZUVtr1Xr17y7LPPygcffCCrVq0Sy7Jk4MCBUlb21Wqi+qzFTnFxsWzatEkWLlwoCxYskGnTprn0kwAAgFgTEwVPQUGBjBw5UrKzs6VRo0ZhsbFjx8p1110n7dq1k549e8qDDz4oubm58umnn5r46tWrZc+ePbJo0SLp3r27DBkyRGbOnCnz5s0zRRAAhAT9bunV3bk6lmIqWvvaEXeqDTtzFHdLd9/48eNNL82AAQOqzTt16pTp7Wnfvr20bt3abNu8ebN07dpVmjdvHsobNGiQnDhxQnJyciK2VVRUZHIqPwD4HHdL90ZMcbd0e3NUwO+W7nrBs2TJEtmxY4fMmjUrYs4TTzwh9erVM48VK1bImjVrJDn5qz/qvLy8sGJHVbzWWCT6/Ro0aBB6VBRQAHws6HdLr+7O1bEUM89R2teOuFNt2JmjAn63dFcLHr00NWHCBFm8eLGkpqZGzNPLXTt37pT169fLZZddJj/84Q/lzJkztfreU6ZMkePHj4ceeiwAfC7od0uv7s7VsRQzz1Ha1464U23YmaMCfrf0RHHR9u3bJT8/34zNqaCDkDds2CBz5841l50SEhJCvTAdO3aUq666yozzefnll+Wmm26SzMxM2bp1a1i7hw8fNs8aiyQlJcU8AACA/7naw9O/f3/ZvXu37Nq1K/To3bu36dHRr7XYOZvO0tKHFkOqb9++pg0tnCroJa/69etL586dHf15AABAbHK1hyc9PV26dOkStq1u3bqSkZFhtn/yySfyxz/+0UxD1zV4/v73v5u1enRNnuuvv97ka0wLm1GjRsmcOXPMuJ2pU6eagdD04AAAgJgYtFwdHdezceNGU9x06NBBbrzxRlMk6Xo7usCg0l6gZcuWmWft7fnRj34kt9xyi8yYMcPtwwcAADHC1R6eqqxbty70dcuWLWX58uXn3adt27YXlAcAAIIppnt4AAAA7EDBAwAAfI+CBwAA+B4FDwAA8D0KHgAA4HsUPAAAwPcoeAAAgO9R8AAAAN+LuYUHASBq4hKc2z9SrlvbE1P0ZoTeiKmE5Ojsa0fcqTbszHGjrZrkOYAeHgDBEenNNxr7R8p1a3tSHZHkNG/EVLT2tSPuVBt25rjRVk3yHEAPD4DgKC92bv9IuW5tLzktIpZHYvpcGKV9bYg71YadOW60VZM8B1DwAAiO8jLn9o+U69b20qKqt8dqrLQ4evvWNu5UG3bmuNFWTfK8eEmroKDA7iYBAACcK3geffTRauMnT56UQYMG1e6IAAAA3Cx47rvvPvn9739fZezUqVMyePBgOXLkiF3HBgAA4HzB84c//EFuv/12Wbp06TnFjvbsfPHFF/LGG2/Yc2QAAABuDFr+j//4Dzl27JjcdNNN8pe//EX69esX6tk5fPiwrF+/Xlq0aGHXsQEAALgzS+snP/mJHD16VIYPHy6vvvqqTJs2TQ4ePGiKnZYtW9pzVAAAADa6qGnp9957ryl6+vfvL+3atZN169ZJq1at7DwuAAAAdwqe73//+2Gvk5KSpEmTJjJhwoSw7S+99JI9RwcAAOB0wdOgQYOw1zqWBwAAwFcFz7PPPhu9IwEAAIgSbh4KAAB8j4IHAAD4HgUPAADwPQoeAADgezUqeHSRwe3bt0fvaAAAANwueP7+97/LkCFDzCKD48aNkxUrVkhxcXE0jgsAAMCdgueZZ56RvLw8ef755yU9PV0mTpxoFh4cMWKEuYu6rr4MAADg+TE88fHxcu2118qcOXNk7969smXLFunTp4/Mnz/f3Evruuuuk0ceeUQ+//zz6BwxAACA04OWr7jiCnNvrbfeektyc3Nl9OjRsnHjRtMLBAAA4Nmbh0bStGlTufXWW80DAAAgVjAtHQAA+F6NCp6//vWvYllW9I4GAADA7YJn0KBB8sUXX0TjOAAAAGKj4KF3BwAAeBFjeAAAgO/VuOB58skn5fXXX5cvv/wyOkcEAADg9rT0uXPnyvTp0yUuLk5at24tPXv2DHtkZmbafYwAAADOFjw5OTlSWloqO3fulB07dphHdna2WXRQiyAteFhlGQAAeLbg0YJG6S0k9DF06NBQ7MiRI+ZO6rt27bL/KAHADnEJzu0fKdet7YkpOvPEGzGVkBydfe2IO9WGnTlutFWTvFgreKqbpZWRkSEDBw40DwCISfrmW1rozP6Rct3anlQn8rHGYiw5LXr71jbuVBt25rjRVk3yYq3gWblypTRo0CB6RwMA0VRe7Nz+kXLd2l5yWv/Z6pGYPhdGaV8b4k61YWeOG23VJC/WCh56bwB4WnmZc/tHynVre2lR1dtjNVZaHL19axt3qg07c9xoqyZ5DmAdHgAA4HsUPAAAwPcoeAAAgO9R8AAAAN+j4AEAAL5HwQMAAHyPggcAAPgeBQ8AAPA9Ch4AAOB7FDwAAMD3KHgAAIDvUfAAAADfi6mCZ/bs2RIXFycTJ040r48ePSp33nmnXH755VKnTh1p06aN3HXXXXL8+PGw/Q4cOCBDhw6VtLQ0adasmUyaNElKS0td+ikAAECsqdHd0qNp27ZtMn/+fMnKygptO3jwoHk88sgj0rlzZ/nss8/kjjvuMNtefPFFk1NWVmaKnczMTNm0aZMcOnRIbrnlFklKSpKHH37YxZ8IAADEipjo4SkoKJCRI0dKdna2NGrUKLS9S5cu8uc//1mGDRsml156qXz729+Whx56SF577bVQD87q1atlz549smjRIunevbsMGTJEZs6cKfPmzZPi4gi3pQcAAIESEwXP+PHjTS/NgAEDzpurl7Pq168viYlfdU5t3rxZunbtKs2bNw/lDBo0SE6cOCE5OTkR2ykqKjI5lR8AAMCfXL+ktWTJEtmxY4e5pHU+//jHP0zvzdixY0Pb8vLywoodVfFaY5HMmjVLpk+fXqtjBwAA3uBqD09ubq5MmDBBFi9eLKmpqdXmag+M9gLpWJ4HHnig1t97ypQppreo4qHHAgAA/MnVHp7t27dLfn6+9OzZM7RNByFv2LBB5s6day47JSQkyMmTJ2Xw4MGSnp4uL7/8shmQXEEHK2/dujWs3cOHD4dikaSkpJgHAADwP1d7ePr37y+7d++WXbt2hR69e/c2A5j1ay12tGdn4MCBkpycLEuXLj2nJ6hv376mDS2cKqxZs8aM89HeIAAAAFd7eLTHRmdiVVa3bl3JyMgw2yuKncLCQjMLq/Lg4qZNm5qCSONa2IwaNUrmzJljxu1MnTrVDISmBwcAAMTEoOXq6GDmLVu2mK87dOgQFtu/f7+0a9fOFD3Lli2TcePGmd4eLZhGjx4tM2bMcOmoAQBArIm5gmfdunWhr/v16yeWZZ13n7Zt28ry5cujfGQAAMCrYmIdHgAAgED18ABA1JSXO7d/pFy3tpdVs/J8LMYi5dixb23jTrVhZ44bbdUkzwEUPACCwyp3bv9IuW5tLy2qenusxiLl2LFvbeNOtWFnjhtt1STPARQ8AIIjLl7EKnNm/0i5bm1PrGbWaizGIuXYsW9t4061YWeOG23VJM8BFDwAgiM+XqTMof0j5bq1PSE58rHGYixSjh371jbuVBt25rjRVk3yHMCgZQAA4HsUPAAAwPcoeAAAgO9R8AAAAN+j4AEAAL5HwQMAAHyPggcAAPgeBQ8AAPA9Ch4AAOB7FDwAAMD3KHgAAIDvUfAAAADfo+ABAAC+R8EDAAB8j4IHAAD4HgUPAADwPQoeAADgexQ8AADA9yh4AACA71HwAAAA36PgAQAAvkfBAwAAfI+CBwAA+B4FDwAA8D0KHgAA4HsUPAAAwPcoeAAAgO8lun0AAOCY8nLn9o+U69b2suKqt8dqLFKOHfvWNu5UG3bmuNFWTfIcQMEDIDiscuf2j5Tr1vbSoqq3x2osUo4d+9Y27lQbdua40VZN8hxAwQMgOOLiRawyZ/aPlOvW9sSUyMcai7FIOXbsW9u4U23YmeNGWzXJcwAFD4DgiI8XKXNo/0i5bm1PSI58rLEYi5Rjx761jTvVhp05brRVkzwHMGgZAAD4HgUPAADwPQoeAADgexQ8AADA9yh4AACA71HwAAAA36PgAQAAvkfBAwAAfI+CBwAA+B4FDwAA8D0KHgAA4HsUPAAAwPcoeAAAgO9R8AAAAN9LdPsAAMAx5eXO7R8p163tZcVVb4/VWKQcO/atbdypNuzMcaOtmuQ5gIIHQHBY5c7tHynXre2lRVVvj9VYpBw79q1t3Kk27Mxxo62a5DmAggdAcMTFi1hlzuwfKdet7YkpkY81FmORcuzYt7Zxp9qwM8eNtmqS5wAKHgDBER8vUubQ/pFy3dqekBz5WGMxFinHjn1rG3eqDTtz3GirJnkOYNAyAADwPQoeAADgezFV8MyePVvi4uJk4sSJoW1PP/209OvXT+rXr29ix44dO2e/o0ePysiRI01Ow4YN5dZbb5WCggKHjx4AAMSqmCl4tm3bJvPnz5esrKyw7YWFhTJ48GC57777Iu6rxU5OTo6sWbNGli1bJhs2bJCxY8c6cNQAAMALYmLQsvbGaNGSnZ0tDz74YFisordn3bp1Ve77wQcfyMqVK03B1Lt3b7Ptt7/9rVx//fXyyCOPSMuWLR34CQAAQCyLiR6e8ePHy9ChQ2XAgAE13nfz5s3mMlZFsaO0nfj4eNmyZUvE/YqKiuTEiRNhDwAA4E+u9/AsWbJEduzYYXpoLkZeXp40a9YsbFtiYqI0btzYxCKZNWuWTJ8+/aK+JwAA8BZXe3hyc3NlwoQJsnjxYklNTXX0e0+ZMkWOHz8eeuixAAAAf3K1h2f79u2Sn58vPXv2DG0rKyszg47nzp1rLjslJCRU20ZmZqZpo7LS0lIzc0tjkaSkpJgHAADwP1cLnv79+8vu3bvDto0ZM0Y6deokkydPPm+xo/r27Wumqmvx1KtXL7Nt7dq1Ul5eLn369InasQPwIG4e6p0YNw+1N0dx81D3pKenS5cuXcK21a1bVzIyMkLbdRyOPvbt22dea4Gk+7Vp08aM07niiivMtPXbbrtNnnrqKSkpKZGf/exn8p//+Z/M0AIQjpuHeifGzUPtzVEBv3loTMzSqo4WMT169DAFjbruuuvM66VLl4ZydAyQ9gppj5FOR7/mmmvMgoUAcM6NNZ3aP1KuW9v1Zo3V3cgxlmKVn+3e1464U23YmaO4eWhsOXu9nQceeMA8qqM9Pc8991yUjwyA53Hz0KrFYoybh9qbo7h5KAAAgL9R8AAAAN+j4AEAAL5HwQMAAHyPggcAAPgeBQ8AAPA9Ch4AAOB7MbcODwBEDbeW8E6MW0vYm6O4tQQABAS3lvBOjFtL2JujAn5rCQoeAMGht12wypzZP1KuW9urW8o/FmMXc0sCbi0ROUdxawkACAhuLVG1WIxxawl7cxS3lgAAAPA3Ch4AAOB7FDwAAMD3KHgAAIDvUfAAAADfo+ABAAC+x7R0AMHBSsveibHSsr05ipWWASAgWGnZOzFWWrY3R7HSMgAEBCstVy0WY6y0bG+OYqVlAAgIVlquWizGWGnZ3hzFSssAAAD+RsEDAAB8j4IHAAD4HgUPAADwPQoeAADge8zSAhAcLDzonRgLD9qbo1h4EAACgoUHvRNj4UF7cxQLDwJAQLDwYNViMcbCg/bmKBYeBICAYOHBqsVijIUH7c1RLDwIAADgbxQ8AADA9yh4AACA71HwAAAA32PQMoDgYB0e78RYh8feHMU6PAAQEKzD450Y6/DYm6NYhwcAAoJ1eKoWizHW4bE3R7EODwAEBOvwVC0WY6zDY2+OYh0eAAAAf6PgAQAAvkfBAwAAfI+CBwAA+B6DlgEEB+vweCfGOjz25ijW4QGAgGAdHu/EWIfH3hzFOjwAEBCsw1O1WIyxDo+9OYp1eAAgIFiHp2qxGGMdHntzFOvwAAAA+BsFDwAA8D0KHgAA4HuM4QEQHExL906Maen25iimpQNAQDAt3TsxpqXbm6OYlg4AAcG09KrFYoxp6fbmKKalA0BAMC29arEYY1q6vTmKaekAAAD+RsEDAAB8j4IHAAD4XkwVPLNnz5a4uDiZOHFiaNuZM2dk/PjxkpGRIfXq1ZMRI0bI4cOHw/Y7cOCADB06VNLS0qRZs2YyadIkKS0tFTfVZpgAgCgJ+rT06qYIx1Ks8rPd+9oRd6oNO3NcmpZ+7NSpKtMibQ9EwbNt2zaZP3++ZGVlhW2/++675bXXXpMXXnhB1q9fLwcPHpTvf//7oXhZWZkpdoqLi2XTpk2ycOFCWbBggUybNk3cMmv5HioeIBYFfVp6dVOEYylW+dnufe2IO9WGnTkutLXy/UPyxbETVaZ9ceykiQeu4CkoKJCRI0dKdna2NGrUKLT9+PHj8n//93/y61//Wr797W9Lr1695NlnnzWFzdtvv21yVq9eLXv27JFFixZJ9+7dZciQITJz5kyZN2+eKYKcVlxaLtkb9zv+fQGcX1kt3/Jqsn+kXLe2lyQmS0lCkidihcWlUdvXjrhTbdiZU+hwWyfPlMj9S3OkyKp6Mrhun/7aHikrtyRQBY9estJemgEDBoRt3759u5SUlIRt79Spk7Rp00Y2b95sXutz165dpXnz5qGcQYMGyYkTJyQnJyfi9ywqKjI5lR92+MPmT8XB/34AauBkWZxj+0fKdWv79kMlsjOvxBOxztNWRW1fO+JOtWFnTmeH2+r6wGo5fKJISiKsflMsiXLo+BnZuv+oBKbgWbJkiezYsUNmzZp1TiwvL0+Sk5OlYcOGYdu1uNFYRU7lYqciXhGLRL9fgwYNQo/WrVvb8vN8drTQlnYA2O8LaeTY/pFy3dr+idVK9lmtPRFT0drXjrhTbdiZ40Zb6oiEf35XOCoNzHP+yTMSiIUHc3NzZcKECbJmzRpJTU119HtPmTJF7rnnntBr7eGxo+hp2zjNPBeXiKTGi8TV7h+UwDm0AzEuyjE7XksUcy52v1bf/aVYf7nror9HTfaPlOvW9u/ddt9XXy9YF/OxG1r2kriDTaKyrx1x5UQbduY43Zb23Pz42W2yqHSA9E94N+yz0LLEbFfN0lOD0cOjl6zy8/OlZ8+ekpiYaB46MPnxxx83X2tPjY7DOXbsWNh+OksrMzPTfK3PZ8/aqnhdkVOVlJQUqV+/ftjDDqP6tpP4OJEr5DnzH1UfgJ3iGrWPeqzWr7vdbB5RyanFfnW+MfrC2rNh/0i5rm1v18c8vBBLS06M2r52xJ1qw86cNIfburZjU2nRIFXWSy95p7xj6LNQn/W1btf4le0bi1PiLMu9j+STJ0/KZ599FrZtzJgxZpzO5MmTTY9L06ZN5fnnnzfT0dXevXtNXMfuXHXVVbJixQr593//dzl06JCZkq6efvppMzVdiyktbC6E9vDopS0dKF3b4kdnac3f8NXA5Q/kZkk+e2xXgvlfNbQUTvxqifgK8bpT0le/Lbok9zn3K0kVSWkgUnT8q9f6tVUiEqf7lIjUaSKS1lik8IjI6SNfv65Q1fb05iKZXUTydouczD/3MM+OV/daf6bqvj5zTCS1YdXbLiRW3XPrq7463ty3a5ZzMduc+LpVL5G/b49uzI7XKlo5tdnvQtuzY/9Y3O6lWDT3tSPuVBt25jjYls7CGrdoh/n6W7Jd/i3xXXmjtJspdtSTP+opg7u0kNq60M9vVwueqvTr18/MtnrsscfM63Hjxsny5cvNVHP9Qe68806zXWdqVUxL1/yWLVvKnDlzzLidUaNGyU9+8hN5+OGHL/j72lnwVBQ9Olur8gBm7fm57dr2MuX6zrVuHwCAWLfy/UNmNpYOUK6gPTv3D+tsS7FTk8/vmL956KOPPirx8fGmh0dnVukMrCeeeCIUT0hIkGXLlpnCqG/fvlK3bl0ZPXq0zJgxw9Xj1qLm5wM7mVlbOpBZx/bo5a7kRNfHiQMA4Agtar7TOdOM6dEByjpmRy9jJWgPgMNirofHLXb38AAAgNj5/Ka7AQAA+B4FDwAA8D0KHgAA4HsUPAAAwPcoeAAAgO9R8AAAAN+j4AEAAL5HwQMAAHyPggcAAPhezN9awikVC07rio0AAMAbKj63z3fjCAqeSnduV3qHdgAA4L3Pcb3FRCTcS+ufysvL5eDBg5Keni5xcXG2Vp5aROXm5nKPrvPgXNUM5+vCca4uHOfqwnGuYuNcaRmjxU7Lli3NzcYjoYfnn/QktWrVKmrt639g/iAuDOeqZjhfF45zdeE4VxeOc+X+uaquZ6cCg5YBAIDvUfAAAADfo+CJspSUFLn//vvNM6rHuaoZzteF41xdOM7VheNceetcMWgZAAD4Hj08AADA9yh4AACA71HwAAAA36PgAQAAvkfBE2Xz5s2Tdu3aSWpqqvTp00e2bt0qfrZhwwYZNmyYWfFSV6x+5ZVXwuI6Rn7atGnSokULqVOnjgwYMEA++uijsJyjR4/KyJEjzeJUDRs2lFtvvVUKCgrCct577z259tprzXnV1TvnzJkjXjNr1iz5xje+YVb3btasmdxwww2yd+/esJwzZ87I+PHjJSMjQ+rVqycjRoyQw4cPh+UcOHBAhg4dKmlpaaadSZMmSWlpaVjOunXrpGfPnmaGRIcOHWTBggXiJU8++aRkZWWFFi3r27evrFixIhTnPEU2e/Zs87c4ceLE0DbO19ceeOABc34qPzp16hSKc67Cff755/KjH/3InA99D+/atau888473niP11laiI4lS5ZYycnJ1jPPPGPl5ORYt912m9WwYUPr8OHDll8tX77c+uUvf2m99NJLOvvPevnll8Pis2fPtho0aGC98sor1rvvvmt997vftdq3b2+dPn06lDN48GCrW7du1ttvv21t3LjR6tChg3XTTTeF4sePH7eaN29ujRw50nr//fet559/3qpTp441f/58y0sGDRpkPfvss+Zn2LVrl3X99ddbbdq0sQoKCkI5d9xxh9W6dWvr9ddft9555x3rqquusq6++upQvLS01OrSpYs1YMAAa+fOneb8N2nSxJoyZUoo55NPPrHS0tKse+65x9qzZ4/129/+1kpISLBWrlxpecXSpUutv/zlL9aHH35o7d2717rvvvuspKQkc+4U56lqW7dutdq1a2dlZWVZEyZMCG3nfH3t/vvvt/71X//VOnToUOjxxRdfhOKcq68dPXrUatu2rfXjH//Y2rJli/m5Vq1aZe3bt88T7/EUPFF05ZVXWuPHjw+9Lisrs1q2bGnNmjXLCoKzC57y8nIrMzPT+p//+Z/QtmPHjlkpKSnmF1rpm4Hut23btlDOihUrrLi4OOvzzz83r5944gmrUaNGVlFRUShn8uTJ1uWXX255WX5+vvnZ169fHzo3+qH+wgsvhHI++OADk7N582bzWt9c4+Pjrby8vFDOk08+adWvXz90fu69917zhl7ZjTfeaAouL9Pfgd/97necpwhOnjxpdezY0VqzZo31rW99K1TwcL7OLXj0w7cqnKtw+j57zTXXWJHE+ns8l7SipLi4WLZv32668yrfr0tfb968WYJo//79kpeXF3ZO9P4neqmv4pzos3Zx9u7dO5Sj+XrutmzZEsq57rrrJDk5OZQzaNAgcznoyy+/FK86fvy4eW7cuLF51t+fkpKSsPOlXe1t2rQJO1/apdy8efOwc6E36svJyQnlVG6jIserv4dlZWWyZMkSOXXqlLm0xXmqml6G0cssZ/9MnK9z6SUXvQx/ySWXmEsteolKca7CLV261Lw3/+AHPzCX7nr06CHZ2dmeeY+n4ImSf/zjH+aNufIfgdLX+gsRRBU/d3XnRJ/1D6myxMREUwRUzqmqjcrfw2vKy8vNGItvfvOb0qVLl9DPon/w+uZQ3fk637mIlKNvyKdPnxav2L17txlDoWMg7rjjDnn55Zelc+fOnKcqaEG4Y8cOM07sbJyvcPphrONpVq5cacaK6Ye2jh3Ru29zrsJ98skn5hx17NhRVq1aJePGjZO77rpLFi5c6In3eO6WDsTIv8bff/99efPNN90+lJh1+eWXy65du0xP2IsvviijR4+W9evXu31YMSc3N1cmTJgga9asMQM+Ub0hQ4aEvtaB8VoAtW3bVv70pz+ZQbcI/4eZ9sw8/PDD5rX28Oj71lNPPWX+HmMdPTxR0qRJE0lISDhnNL++zszMlCCq+LmrOyf6nJ+fHxbX2Q46qr9yTlVtVP4eXvKzn/1Mli1bJm+88Ya0atUqtF1/Fr00euzYsWrP1/nORaQcnSHhpTd0/Ze2zm7p1auX6bno1q2b/OY3v+E8nUUvw+jfkM4I0n8560MLw8cff9x8rf9S5nxFpr05l112mezbt4/frbPozCvtVa3siiuuCF0CjPX3eAqeKL456xvz66+/HlYd62sddxBE7du3N7+slc+JdunqdduKc6LP+uaib9oV1q5da86d/surIkenv+u19Qr6r1ntAWjUqJF4hY7r1mJHL83oz6jnpzL9/UlKSgo7X3oNW99cKp8vvdRT+Q1Ez4W+kVa8MWlO5TYqcrz+e6i/E0VFRZyns/Tv39/8rNobVvHQf5Xr2JSKrzlfken06I8//th8uPO7FU4vuZ+9dMaHH35oesQ88R5fqyHPOO+0dB2dvmDBAjMyfezYsWZaeuXR/H6jM0N0aqY+9Nfr17/+tfn6s88+C01Z1HPw6quvWu+99541fPjwKqcs9ujRw0x7fPPNN81Mk8pTFnXUv05ZHDVqlJmyqOdZp3x6bVr6uHHjzPTNdevWhU2JLSwsDJsSq1PV165da6bE9u3b1zzOnhI7cOBAM7Vdp7k2bdq0yimxkyZNMjNM5s2b57kpsb/4xS/M7LX9+/eb3xt9rbM6Vq9ebeKcp+pVnqWlOF9f+/nPf27+BvV366233jLTy3Vauc6aVJyr8GUOEhMTrYceesj66KOPrMWLF5ufa9GiRaGcWH6Pp+CJMl1vQf9YdD0enaau6w742RtvvGEKnbMfo0ePDk1b/NWvfmV+mbUY7N+/v1lXpbIjR46YX/569eqZqZ1jxowxhVRlur6DTo/UNv7lX/7F/JF5TVXnSR+6Nk8FfZP46U9/aqZo6h/89773PVMUVfbpp59aQ4YMMetU6Bu1voGXlJSc89+le/fu5vfwkksuCfseXvBf//VfZv0PPX79MNHfm4piR3GealbwcL7Cp4e3aNHC/Az6XqKvK68rw7kK99prr5kCT997O3XqZD399NNh8Vh+j4/T/7v4/iEAAIDYxxgeAADgexQ8AADA9yh4AACA71HwAAAA36PgAQAAvkfBAwAAfI+CBwAA+B4FDwAA8D0KHgAA4HsUPAAQRXqjRL1hZ/fu3aVLly6SnZ3t9iEBgcStJQAgisrKysxd3dPS0uTUqVOm6HnnnXckIyPD7UMDAoUeHgD4p379+snEiRNtbTMhIcEUO0oLn3/etNnW7wHg/Ch4AMSkMWPGyNSpU8UPx6KXtbp16yatWrWSSZMmSZMmTWw9PgDnR8EDICYvAy1btky++93v+uJYGjZsKO+++67s379fnnvuOTl8+LCtxwjg/Ch4AETdm2++KVdeeaWkpqaa3o3f/OY31eZv2rRJkpKS5Bvf+EbES0933nmnufzUqFEjad68uRkMrGNktDcmPT1dOnToICtWrAjto5eT7rrrLmnWrJk5jmuuuUa2bdt23mOvfCwX830r03zt6dm4ceN5vy8Ae1HwAIiq5cuXy/e+9z356U9/Ku+9957cfvvtcvfdd8unn34acZ+lS5fKsGHDJC4uLmLOwoULTfG0detWU4SMGzdOfvCDH8jVV18tO3bskIEDB8qoUaOksLDQ5N97773y5z//2eyncS1MBg0aJEePHq32+M8+lpp+X+3NOXnypPn6+PHjsmHDBrn88ssv6lwCqAWdpQUA0XD69GmrVatW1uLFi0PbSktLrXr16lkLFy6MuF/Hjh2tZcuWRYx/61vfsq655pqwNuvWrWuNGjUqtO3QoUM6MtjavHmzVVBQYCUlJYUdR3FxsdWyZUtrzpw5Ye1OmDAh4rHU9PuqLVu2WN26dbOysrKsrl27Wk899dR5zhqAaEisTbEEANVZu3atnD59Wm688cawWUvaW5KSklLlPh988IEcPHhQ+vfvX23bWVlZYW3qNO+uXbuGXT5S+fn58vHHH0tJSYl885vfDMX1MpVeZtPvF0lVx1KT76v0e+zatavanwVA9HFJC0DUvPHGG2bBPS0MKuzbt89c4unRo0fES0jf+c53zDib6mjBUpkWUZW3VVyCKi8vv+jjr+pYnPi+AOxHwQMganbu3CnFxcVh25544gnp1auXXHbZZVXu8+qrr8rw4cNtPY5LL71UkpOT5a233gpt0x4fHbTcuXPniPtF41gAuINLWgCiWvDoInu///3vpU+fPvLCCy/Ik08+aWY+VUUvA+kqxNqzYqe6deuawcW6Bk7jxo2lTZs2MmfOHDOw+NZbb3X0WAC4g4IHQFQcOHDAzIDSNWx+8YtfyIcffmjGv6xcuTLi5azXXnvNjHmJxsJ8s2fPNpeZdAaVXlLT+1utWrXKTC93+lgAOI97aQGICu0Z0bVpjhw5csH76OJ+uj6OTiF3WywdC4DaYwwPgKhdzqo8e+lCaIFx0003SSyIpWMBUHv08ACIihtuuMGMlXn88cfdPhQAoOABAAD+xyUtAADgexQ8AADA9yh4AACA71HwAAAA36PgAQAAvkfBAwAAfI+CBwAA+B4FDwAA8D0KHgAA4HsUPAAAwPcoeAAAgO9R8AAAAN+j4AEAAOJ3/w8OBCO/LXWvowAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Here is the bounding region and the points used for construction of the K-D tree\n", "# This all looks nice\n", "for reg in rf.get_regions_rw():\n", " reg.add_pair(proppair=teqpflsh.PropertyPairs.DT, Nsplit=5)\n", " \n", " pset = reg.propset_bounding\n", " plt.plot(pset.rho, pset.T, 'o-')\n", "\n", " pset = reg.propset_Trhogrid\n", " plt.plot(pset.rho, pset.T, '.')\n", "plt.gca().set(xlabel=r'$\\rho$ / mol/m$^3$', ylabel='$T$ / K');" ] }, { "cell_type": "markdown", "id": "673b70e9-dde8-4710-a560-73fa14ca07bc", "metadata": {}, "source": [ "But when you shift to another variable pair, here density and entropy, the rectangular box and its sampled (in $T$, $\\rho$) points do not all map together" ] }, { "cell_type": "code", "execution_count": 4, "id": "dc9f7ae0-dace-420e-9826-b5d9b8c6d458", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:49.138115Z", "iopub.status.busy": "2025-01-06T11:32:49.138026Z", "iopub.status.idle": "2025-01-06T11:32:49.893696Z", "shell.execute_reply": "2025-01-06T11:32:49.893435Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAG4CAYAAACq8YbKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcVBJREFUeJzt3QdYlWUfBvCbjWwBQRH3QhQB98IyTXHvlZmaZpmLLC0bZlM/LXOm2VArzZF7oeYIVNziFrfiQFBkyz7f9TwEgoKLBziHc/+u63yn95yX87ycDw83z/obaDQaDYiIiIj0hGFRXwARERFRYWL4ISIiIr3C8ENERER6heGHiIiI9ArDDxEREekVhh8iIiLSKww/REREpFcYfoiIiEivGBf1BWij9PR03Lp1C9bW1jAwMCjqyyEiIqJnIPZtjo2NhYuLCwwN8+7fYfjJhQg+5cqVK+rLICIiohcQGhoKV1fXPJ9n+MmF6PHJfPNsbGyK+nKIiIjoGcTExMjOi8zf43lh+MlF5lCXCD4MP0RERLrlaVNWOOGZiIiI9ArDDxEREekVhh8iIiLSKww/REREpFcYfoiIiEivMPwQERGRXmH4ISIiIr3C8ENERER6hZscFpK0dA0OXolEeGwinKzN0bCSPYwMWTeMiIiosDH8FAL/U7fxxYYzuB2dmPVYGVtzfN7JHb61yxTptREREekbDnsVQvAZ/udRGXzmYRpOmAyW92HRifJx8TwREREVHoafAh7qEj0+GgCXTF+Dr9kx2BglyfuLpq/Jc8Tz4jwiIiIqHAw/BUjM8cns8RHTezLrrIl7cfwjpsnnxXlERERUOBh+CpCY3Cy0MTmWFXwyiWPxePbziIiIqOAx/BQgsapL5XlERESUfww/BUgsZxerup7EKNtwGBEREelZ+Jk3bx7q1KkDGxsbeWvSpAm2bNmS9XxiYiJGjBgBBwcHWFlZoUePHrhz506O17h+/To6dOgACwsLODk5Ydy4cUhNTS2C7wZyHx+xnP1J0jTAaz/vx7St55CSll5o10ZERKSvtCr8uLq6YsqUKThy5AgOHz6MV155BV26dMHp06fl8++99x42bNiAlStX4t9//8WtW7fQvXv3rK9PS0uTwSc5ORn79u3D4sWLsWjRIkycOLHIvie5j49RHk8aAT3qukIs9pq76xJ6zNuHSxFxhXyFRERE+sVAo9Fo9Tpre3t7TJs2DT179kSpUqWwdOlS+d/CuXPnULNmTQQFBaFx48ayl6hjx44yFDk7O8tz5s+fjw8//BAREREwNTV9pjZjYmJga2uL6Oho2QOVb5Nsn/BcNDaduI2P15xE9IMUmJsY4tMO7ujfqDwMOB5GRET0zJ7197dW9fxkJ3pxli1bhvj4eDn8JXqDUlJS0Lp166xz3NzcUL58eRl+BHHv4eGRFXyEtm3byjcjs/coN0lJSfKc7LfC1KFOGWz1a4FmVR2QmJKOT9eewtDFh3E3LqlQr4OIiEgfaF34OXnypJzPY2ZmhnfeeQdr1qyBu7s7wsLCZM+NnZ1djvNF0BHPCeI+e/DJfD7zubxMnjxZJsXMW7ly5VDYStua4483G+HTDjVhamSIHefC4TsjADvO5pzTRERERMUs/NSoUQPBwcE4cOAAhg8fjoEDB+LMmTMF2uaECRNkF1nmLTQ0FEXB0NAAQ30qY/2oZnArbY27cckYsvgwPllzEgnJRTNpm4iIqLjRuvAjeneqVq2KevXqyR4ZT09PzJw5E6VLl5YTmaOionKcL1Z7iecEcf/o6q/M48xzciN6mTJXmGXeipJbaRusHdEMQ5pXksdLDlxHx1l7cOJGzu+diIiIikH4eVR6erqckyPCkImJCXbs2JH1XEhIiFzaLuYECeJeDJuFh4dnnbN9+3YZZsTQmS4xNzHCZx3d8eeQRnC2McPlu/Ho/uM+zN11kbXAiIiIikv4EcNPAQEBuHr1qgwx4nj37t3o37+/nIszZMgQjB07Frt27ZIToAcPHiwDj1jpJbRp00aGnAEDBuD48ePYunUrPv30U7k3kOjd0UXNqznKydDtPUojNV2DaVtD0OenIIRGJhT1pREREekkrQo/osfmjTfekPN+WrVqhUOHDskA8+qrr8rnf/jhB7mUXWxu2KJFCzmUtXr16qyvNzIywsaNG+W9CEWvv/66fL0vv/wSuszOwhRzX6uL73t5wsrMGIev3Ue7mYFYdeQGtHynAiIiIq2j9fv8FAXV+/ykTbLNdZ/DNBHYJkU/12uJHp/3lgfLAJS5TP6brrVlQCIiItJnMbq+z09xMXnzmYyUk5u0/55/DuXsLbBsWGN80KY6jA0N5AaJvjMCsffiXSXXS0REVNwx/BSg5NR0/Bx45YnniOfFec/D2MgQI1+phlXDm6KyoyXCYhLR/5cD+HrjGSSm5JW0iIiISGD4KUB/BF2VdbueRDwvznsRnuXssHF0c1kKQ/hlzxV0nbsX58IKd4dqIiIiXcLwU4CuPeOKrGc9LzcWpsb4ppsHfnmjPhwsTXEuLBad5+zFr3uuIJ1L4omIiB7D8FOAKthbKD3vSVq7O8PfrwVecXOSw2hfbTyDN347iLDoxHy/NhERUXHC8FOABjSpCMNnKMzeqLKDkvZKWZvh14H18XXX2rI6/J6Ld9F2RgA2n7yt5PWJiIiKA4afAmRqbIi3fDJKVDxJj3n7sGjvFSV79hgYGOD1xhWwabQPPMraIvpBCt5dchQfrDyO2MSUfL8+ERGRrmP4KWAT2rsj101+BCPgpeqlkJSajkkbzmDQwkMIj1UzTFWllJVcDTaiZRUYGAB/H7mB9rMCcfhqpJLXJyIi0lUMP4XgCdkHiwY3wBeda8HM2BD/no+Qe/ZsP5OzOGt+ep7GtXXD8mFNUNauBEIjH6D3T0H4flsIUtKeb3k9ERFRccHwU8TEMNXAphWxYVRz1Cxjg8j4ZLz1+2FMWH0SCcmpStpoWMkeW/x80L1uWbm0fvbOi+g5bx8uR8QpeX0iIiJdwvCjJao7W2PtiKYY1qKyPP7r4HV0nLUHJ25EKXl9G3MTTO/thTmvecPG3BjHb0Sjw6w9WHrgOuuDERGRXmH40SJmxkb4uH1NLB3aCKVtzHH5bjy6/7gPc3ddRJqiPXs61nHB1vdaoGkVBzxIScPHa07Knqa7cUlKXp+IiEjbMfxooaZVHeHv54MOHmWQmq7BtK0h6LdgP27cf/HNELMrY1sCfw5phE871ISpkSH+ORsO3xkB2HlOzVwjIiIibcbwo6VElXYxRPVdL09Ymhrh4NVItJsRiLXHbip5fUNDAwz1qYy1I5qhurMV7sYl481Fh/Hp2pN4kMz6YEREVHwx/Gj5ZOie9VyxeYwP6pa3Q2xSKvyWB2P0X8fk/j0quLvYYP3I5nizWcZ+RH/uv44OswNx8ka0ktcnIiLSNgw/OqCCgyVWvN0Efq2rwcjQAOuP30L7mYE4cPmektc3NzHCxE7u+GNIQzjbmOFyRDy6/bhX6VwjIiIibcHwoyOMjQzh17q6DEHl7S1wM+oB+v68H//zPydreangU60U/Me0QLvapXPMNQrNR+FVIiIibcPwo2PqVSgph8F61XOFWKE+b/clWR7jkqI9e0pamuLH/nUxrWedrLlGopdpzbEbXBJPRETFAsOPDrIyM8a0Xp4ypNiWMMHJm2LPnkD8uf+asvpgveqXw5YxLWTYEnON3lt+HKPEXKME1gcjIiLdxvCjw9p7lJFL4ptVdUBiSjo+XXtK7tlzT9GePeUdLLB8WGOMfbW6nGu08cRt+M4MwL6Ld5W8PhERUVFg+NFxYs+eP95shE/aP9yzp+2MQOwKCVc212h0q2qySGpFBwvcjk5E/18P4NvNZ5GUyiXxRESkexh+igGxZ89bLTL27KnmJPbsScLghYfw+bpTSExRE1C8ytlh02gf9GtYXs41WhBwGV3m7EVIWKyS1yciIiosDD/FiNizRxRIHdS0ojxeHHQNnWbvwelbavbssTQzxuTuHvj5jfqwtzTFubBYdJqzB7/tuYJ0LoknIiIdwfBTzIg9eyZ1roVFgxvA0coMF8Lj0G3uPiwIuKQsoLzq7iznGr1co5RcZv/lxjMYuPAg7sQkKnl9IiKigsTwU0y9XMMJW/180LqmM5LT0vHt5nN4/dcDuB39QMnrO1mbY+GgBviqSy2YGRsi8MJdtJ0RAP9Tt5W8PhERUUFh+CnGHKzM8PMb9fBtNw+UMDHCvkv34DsjEJtPqgkoYkn8gCYVsWl0c9Qua4OohBS88+dRjFt5HHFJqUraICIiUo3hp5gTAeW1RuVlQKnjaitrgr275Cg+UBhQqjpZY/XwZhj+chUYGAArj9yQGyMeuXZfyesTERGpxPCjJyqXspLL1Ue0zAgofysOKKbGhvjQ1w3L3mqMsnYlcD0yAb3m78P07eeRkqam/AYREZEKDD96xMTIEOPaumH5sCZZAaX3T0H4Yft5pCoKKI0qO2CLnw+6eZeFmF89a8cF9JwfhCt345W8PhERUX4x/OihhpXsZUDp6uUiq7bP3HEBvX4KwrV7agKKjbkJfujjhVn9vGFjbozjoVGy/Mayg9dZH4yIiIocw4+eEgFlRl9vzOzrBWtzYxy7HiWHwVYcDlUWUDp7usDfrwWaVHZAQnIaPlp9EsP+OKKs/AYREdGLYPjRc128ymLLGB/ZGxSfnIbxf5+QE6LvxycreX0XuxJYMrQRPm7vBhMjA2w/c0dp+Q0iIqLnxfBDcC1pgb/eaozxvjVgbGiALafCZAHTPRfuKiu/MaxFlcfKb0xcdwoPklkfjIiIChfDD0miavu7L1fF6neborKjJe7EJMlNEb/eeEZZAdNaLray/MbgZhnlN34X5Tfm7MGpm2rKbxARET0Lhh/KoY6rHTaObi73BhJ+2XNFFjA9fydWWfmNzzvVwuI3G6KUtRkuivIbP+7FvN2X5ORrIiKigsbwQ4+xMDWWu0JnL2DacfYeLNx7Rdlk6Jeql8JWvxZoW8sZKWka/M//HPr9vB837icoeX0iIqK8MPzQUwuYiqAiCph+sUEUMD2EcEUFTEWwmv96PUztWQeWpkY4eCUS7WYGYl3wTSWvT0RElBuGH3pqAVNRIf6LzhkFTAPOR8B3ZiC2nQ5TVn6jd/1y2DzGB97l7RCbmIoxy4Ix+q9jshQHERGRagw/9EwBZWDTinKycs0yNoiMT5b79UxYfQIJyWrqg1VwsMTKt5vgvdbV5eTr9cdvod2MAARduqfk9YmIiDIx/NAzq+5sjbUjmmJYi8ry+K+Doeg4aw9O3IhS8vrGRoYY07oa/n6nCSo6WOBWdCJe+2U/Jm8+q2zFGREREcMPPRczYyN83L4mlg5thNI25rh8Nx7df9yHubsuKlut5V2+JDaN9kHfBuUg5lf/FHAZ3ebuwwVFK86IiEi/MfzQC2la1VFOhu7gUQap6RpM2xqCfgvUrdayNDPGlB518NOAeihpYYIzt2PkirNFClecERGRfmL4oRdmZ2GKOa9547tenhmrta5Got2MQKw9pm61VttapeWSeLHiLCk1HZMUrzgjIiL9w/BD+Z4M3bOeq1ytVVes1kpKhd9ytau1nGweX3HWdkYA/E+pWXFGRET6heGHlK3WWvF2E/i1rpa1WktUiT9w+Z7SFWcbRzWHexkb3E9IwTt/HsGHf59AfJKaFWdERKQfGH5IGbFay691dax8pwnK21vgZtQD9P15v9y9WWySqEI1ueKsGd55qQoMDIDlh0PRflYgjl6/r+T1iYio+GP4IeXqli8ph8F61XOVq7VE3a4e8/bhUkScktc3NTbER+3cZCV6F1tzXLuXgF7zgzDjn/NITVMTsoiIqPhi+KECYWVmjGm9PDGvf13YljDByZvR6DArEH/uv6ZstVbjyg7Y4tcCXbxc5DL7Gf9cQM/5Qbh6N17J6xMRUfHE8EMFqp1HGblaq1lVBySmpOPTtafw1u+HcS8uScnri2A1s683Zvb1grW5MYJDo+Qw2PJD17kknoiIcsXwQwWutK05/nizET5pXxOmRob452w42s4IxK6QcGVtdPEqC3+/FmhUyR4JyWn4cNVJOSFalOIgIiLKjuGHCoWhoQHealFZTlau5mSFu3FJGLzwED5fdwqJKWpKV5S1K4GlbzWW84FMjAyw9fQduST+3/MRSl6fiIiKB4YfKlTuLjayQOqgphXl8eKga+g0ew9O34pW8vpimb1YCbbm3Wao6mSFiNgkDPztICatP60sZBERkW5j+KFCZ25ihEmda8mNCx2tzHAhPA5d5+7FgoBLSFdUH6x2WVu5J1BmyFq076rSkEVERLqL4YeKzMs1nLDVzwetazojJU2Dbzefw+u/HsDt6AdKQ9bCwQ1QyvphyPrp30vKirASEZHuYfihIuVgZYaf36iHb7t5oISJEfZdugffGYHYdOK2sjZa1nCC/xgfvOqeEbImbzmH/r/sl5swEhGR/mH4oSInSle81qg8No1ujjqutrIm2IilR/H+iuOIU1S6QoSsBQPq4X89PGBhaoT9lyPhOyNAluEgIiL9wvBDWqNyKSusGt4UI1pmlK5YdfSGrA925Np9ZSGrT4Py2DzaB17l7BCbmCoLsPotU1eElYiItB/DD2kVEyNDjGvrhuXDmsil69cjE9D7pyD8sF1d6YqKjpay/tiYVtVgaACsDc4owrpfURFWIiLSbgw/pJUaVrLHFj8fdP2vdMXMHRfQ66cgXLsXryxkvfeqKMLaNKsIa7+f92PKFnVFWImISDsx/JDWsjE3wYxspSuOXY+SPTQrDocqK11Rr0JGEdY+9cvJIqzz/72Ebj/uxcXwWCWvT0RE2ofhh7SeKF2xZYyP7A2KT07D+L9P4N0lR3FfUekKUYT1fz3rYP7rdWFnYYLTt2LQYdYe/B50lfXBiIiKIYYf0gmuJS3w11uNMd63BowNDbDlVBh8ZwZgz4W7ytrwrZ1RhNWnmiOSUtMxcd1pDF50COGxicraICKiosfwQzpDlK549+WqsnRF5VKWuBOTJDdF/HrjGSSlqild4WxjjsWDG2JSJ3eYGhtid0iE3Hdo2+kwJa9PRERFj+GHdI6Ha0bpiv6NysvjX/ZcQZc5e3H+TqyyIqyDmlWSbdQsYyMrww/74wg+WnUC8Yr2HSIioqLD8EM6ycLUGN9088Avb9SHg6UpzoXFouPsPVi494qyeTrVna2xdkRTvN2istx3aNmhUHSYFYhj19XsO0REREWD4Yd0Wmt3Z7kk/qXqpeQS9S82nMHAhYcQHqNmno6ZsREmtK+JJUMboYytOa7eS0DP+UGY+c8FZfsOERFR4WL4IZ3nZG0uK8R/0bkWzIwNEXA+Ar4z1c7TaVrFEf5jWqCTZ8a+Qz/8c15uvqhq3yEiIio8DD9ULIjSFQObVsSGR+bpTFh9AgnJaubp2FqYYFZfL8zo4wVrM2McLYB9h4iIqOAx/FCxkjlPZ1iLyvL4r4Nins4eHA+NUhayunqXlUNt2fcdGv6nun2HiIhIj8LP5MmT0aBBA1hbW8PJyQldu3ZFSEhIjnNefvll+Qso++2dd97Jcc7169fRoUMHWFhYyNcZN24cUlO5SkdfiHk6H7eviaVDG6G0jTmu3I1Hj3n7MGfnBTlkVRD7DvmfDkPbGQFyyI2IiLSbVoWff//9FyNGjMD+/fuxfft2pKSkoE2bNoiPzzmv4q233sLt27ezblOnTs16Li0tTQaf5ORk7Nu3D4sXL8aiRYswceLEIviOqCg1reoIfz8fdPAog9R0Db7bdh59FwQhNDJB6b5Da0c0Q5VSlgiPTcIbvx3EFxtOIzFFzb5DRESknoFGiycrREREyJ4bEYpatGiR1fPj5eWFGTNm5Po1W7ZsQceOHXHr1i04OzvLx+bPn48PP/xQvp6pqelT242JiYGtrS2io6NhY2OT/29kku0TnovO/+vTE4kf8VVHb+LzdafkMJWYr/NV19py+EqVB8lpmLzlLH4PuiaPqztbYUYfb7i7KPj5ISKiZ/Ksv7+1qufnUeLiBXt7+xyPL1myBI6OjqhduzYmTJiAhISHf8kHBQXBw8MjK/gIbdu2lW/I6dOnc20nKSlJPp/9RsWHGBrtWc9VFjCtW94OsUmp8FsejNF/HUP0gxQlbZQwNcKXXWpj4aAGcLQyw/k7ceg6dy8WBFxCuqKhNiIiUkNrw096ejr8/PzQrFkzGXIyvfbaa/jzzz+xa9cuGXz++OMPvP7661nPh4WF5Qg+QuaxeC6vuUYiKWbeypUrV2DfFxWdCg6WWPF2E/i1riaHrNYfvyVXax24fE9ZGy3dnORQW+uazkhOS8e3m8+h/y8HcCvqgbI2iIiomA57DR8+XA5h7dmzB66urnmet3PnTrRq1QoXL15ElSpVMGzYMFy7dg1bt27NOkf0DFlaWmLz5s1o165drj0/4pZJ9PyIAMRhr+Lr6PX78FsWjOuRCXL35ndeqoL3WleX9bxUEP+sxI7QX244gwcpabAxz9iRWuwTREREBUOnh71GjhyJjRs3yt6dJwUfoVGjRvJehB+hdOnSuHPnTo5zMo/Fc7kxMzOTb1L2GxVvdcuXlMNgveq5QsT/ebsvyRVhlyLilA219WtYXrbh6WqLmMRUjPrrGN5bHoyYRDVDbURE9GK0KvyIv5ZF8FmzZo3s0alUqdJTvyY4OFjelylTRt43adIEJ0+eRHh4eNY5YuWYCDTu7u4FePWka6zMjDGtlyfm9a8L2xImOHkzWtbu+nP/NWWbFlZytMTfw5ti9CtVYWgArDl2E+1mBOLglUglr09ERDo+7PXuu+9i6dKlWLduHWrUqJH1uOjCKlGiBC5duiSfb9++PRwcHHDixAm89957sndIrAjLXOouVoO5uLjIJfBins+AAQMwdOhQfPvtt890HVztpX/CohPx/spg7L2YMf+ndU0n/K9HHThYmSlr48i1SDnROjTygRxqG/5SFfgpHGojItJ3Mc/4+1urwo8YKsjNwoULMWjQIISGhsrJzadOnZJ7/4h5Od26dcOnn36a45sUc37EnKHdu3fLuT4DBw7ElClTYGxs/EzXwfCjn8SqrN/2XsFU/xA5WVms2prWqw5a1nBS1kZcUiq+WH8aK4/ckMe1y9rIJfFVnayUtUFEpK9idDH8aAuGH/125lYM/JYfk8vVhYFNKsjK7uYmRsra2HzyNj5ecxJRCSkwNzHEJx3c8Xqj8nn+AUBERMV8wjNRURIbE64f2RyDmlaUx4uDrqHT7D04fUtdUG3vUUZWifep5ojElHR8tvYU3lx0CBGxD1cdEhFRwWD4IcqF6OWZ1LkWFr/ZEKWszXAhXP2mhaVtzbF4cENM7Ogu5/3sComA74wA/HMm52pFIiJSi+GH6Aleql4K/mN88Kq7M1LSNHLTwtd/PYDb0Wo2LTQ0NMCbzSthw8jmcCttjXvxyRj6+2E5JJaQzGK8REQFgeGH6CnEiq8FA+rh224eKGFihH2X7sF3RiA2nbitrI0apa2xbmQzDGtRWR4vPXAdHWbtwfHQKGVtEBFRBoYfomcgJiK/1qg8No1ujjqutrIm2IilR/H+iuOIVbRpoZmxET5uXxNLhzZCaRtzXLkbj+7z9mH2jgtITUtX0gYRETH8ED2XyqWssGp4U4xoWUXu1bPq6A20nxUo9/BRpWlVR1kfrEOdMkhL1+D77efRZ8F+XL/3sIAvERG9OIYfoudkYmSIcW3dsHxYE5S1KyE3Lew1PwjTt59X1kNjZ2GKOf288UMfT1ibGePItfsyZP195Iay3aeJiPQVww/RC2pYyR5b/HzQ1csFYgHYrB0X0HN+EK7ejVc21NbN21XWB2tQsaTcIPGDlcflcNv9+GQlbRAR6SOGH6J8sDE3wYy+3pjZ1wvW5sYIDo2SPTQrDoUq66EpZ2+BZcOaYFzbGjA2NMDmk2HwnRmAwAsRSl6fiEjfMPwQKdDFqyy2jPGRvUEJyWkYv+oE3l2irofGyNAAI1pWxZp3m6FyKUvciUnCgF8P4ssNZ5CYkqakDSIifcHwQ6SIa0kL/PVWY4z3zeih2XIqo4dmz4W7ytrwcLXFplE+eL1xeXksapF1mbMXZ2/HKGuDiKi4Y/ghUkj00Lz7cs4eGrEp4tcbzyApVU0PTQlTI3zd1QO/DaoPRytThNyJlQHol8DLynafJiIqzhh+iAqA6KHZOKo5+jfK6KH5ZU9GD835O7HK2njFzRn+fi3QuqaTrEL/9aazSnefJiIqrhh+iAqIhakxvunmgV/eqA8HS1OcC4tFx9l7sHDvFWWToR2tzPDzG/XxTbfasjp8Qew+TURU3DD8EBWw1u7Ockn8yzVKITk1HV9sOIOBCw8hPCZR2ZL4/o0qYNNonxy7T49dEaxs92kiouKE4YeoEDhZm2PhoAb4skstmBkbIuB8BHxnBmLb6TBlbVT5b/fpUa9UhaEBsProTbSbGYhDV9XtPk1EVBww/BAVEtFD80aTinIuUM0yNoiMT8awP45gwuoTyiq4i92n329TAyvebgLXkiVw4/4D9PkpCNO2nkMK64MREUkMP0SFrJqzNdaOaCoruIv6YH8dDFVewb1+RXu571CPuq5y9+m5uy6hx7x9uBQRp6wNIiJdxfBDVAQyK7gvGfKwgrsIJ3N2XpDFTFWwNjfB9709Mfe1urAtYYITN6LRYVYg/tx/jfXBiEivMfwQFaGsCu4eZZCarsF3286j74IghEaqq+AuqsNv9WuBZlUdkJiSjk/XnsLQxYdxNy5JWRtERLqE4YeoiMkK7q9547tenrA0NcKhq/fRfmYg1h67qayN0rbm+OPNRvi0Q02YGhlix7lw+M4IwI6zd5S1QUSkKxh+iLRkMnTPehkV3OuWt0NsUir8lgdj9F/H5NJ1FQwNDTDUpzLWj2oGt9LWuBuXjCGLD+OTNSeVTbgmItIFDD9EWqSCg6VcqeXXuposlbH++C3ZC7T/8j1lbbiVtsHaEc0wpHklebzkwHV0nLUHJ26om3BNRKTNGH6ItIyxkSH8WlfHyneaoLy9BW5GPUC/n/fjf/7n5CaJKpibGOGzju74c0gjONuY4fLdeHT/cR/m7rqobMI1EZG2Yvgh0lJ1y5eUw2C96rlCLM6at/sSus/bi4vh6parN6/mKCdDt/coLSdcT9saIvcFUjnhmohI2zD8EGkxKzNjTOvliXn9M5arn7oZg46z1S5XFxOuxXL473t5yvYOX7svd4ZedeQGl8QTUbHE8EOkA9p5PL5c/a3f1S1XFxOue9RzlRsj1q9QEnFJqXh/5XGM/OsYohKSlbRBRKQtGH6IdMSjy9X/OSuWqwdiV0i4sjbK2Vtg2bDG+KBNdRgbGsjq8KKNvRfvKmuDiKioMfwQ6ZDM5epitVZ1ZyvZ8zN44SF8vu4UElPSlE24HvlKNVkktbKjJcJiEtH/lwP4euMZZW0QERUlhh8iHeTuYoP1I5tjUNOK8nhx0DV0mr0Hp29FK2vDs5wdNo5ujv6NysvjX/ZcQde5e3EuLEZZG0RERYHhh0hHieXqkzrXwuI3G6KUtRkuhMfJcLIg4BLSFS1XtzA1xjfdPPDLG/XhYGmKc2Gx6DxnL37dc0VZG0REhY3hh0jHvVS9FPzH+OBVd2ekpGnw7eZzeP3XA7gd/UBZG63dneHv1wKvuDnJvYa+2ngGb/x2EGHRicraICIqLAw/RMWAg5UZFgyoh8ndPVDCxAj7Lt2TE5XFhGVVRO/SrwPr4+uutWFuYog9F++i7YwAbD6prg0iosLA8ENUTIjl6v0alsem0c1Rx9VW1gQbsfQo3l9xHLGJKcraeL1xBWwa7QOPshltvLvkKD5Yqa4NIqKCxvBDVMxULmUlV2qNbFkVhgbAqqM30H5WII5ci1TWRpX/2hjRsgoMDIC/j2S0cfiqujaIiAoKww9RMWRiZIgP2tbAsmFNUNauBEIjH6DX/CBM334eqWlq6oOZGhtiXFs3LM/WRu+fgvD9thCkKGqDiKggMPwQFWMNK9lji58Punq5QCzOmrXjAnrOD8LVu/HK2+het6xsY/bOi+g5bx8uR6irQUZEpBLDD1ExZ2Nughl9vTGzrxeszY0RHBolh6hWHApVVrtLtDG9txfmvOYNG3NjHL8RjQ6z9mDpgeusD0ZEWofhh0hPdPEqK2t3iZ6ahOQ0jF91AsP/PIr78epqd3Ws44Kt77VA0yoOeJCSho/XnFRag4yISAWGHyI94lrSAn+91RjjfWvI2l3+p8PgOzMAey6oq91VxrYE/hzSCJ+0z16DLAC7zqmrQUZEpHXhJy6OY/1E2srI0ADvvlwVa95thsqlLHEnJkluiviVwtpdogbZWy2y1yBLxuBFh/DZ2lN4kMz6YESkY+Hnhx9+eOLzsbGxaNu2bX6uiYgKgYerLTaOeli7S5SsEOUxQsJildcge7NZJXn8x/5r6DA7ECdvqKtBRkRU4OHn448/xu+//57rc/Hx8fD19cW9e/ee+0KIqPDlVrur05w9+E1h7S5Rg2xiJ3f8MaQhnKzNcDkiHt1+3Iu5uy4ijfXBiEgXws8ff/yBt99+G+vXr38s+Igen4iICOzatUvlNRJRARO1u8Ry9ZdrlJK1u77ceAaDFh1CeIy62l0+1Uphq18L+NYqjdR0DaZtDUG/BfsRGpmgrA0iogIJPz179sTs2bPRr18/7N69O0ePz507d+RjZcqUed6XJaIi5mRtjoWDGuDLLrVgZmyIgPMR8J0ZiG2nw5S1UdLSFPNer4tpPevA0tQIB69Gov3MQKw5doNL4olIuyc8Dx06FJ9//jm6dOkiw067du1w69Yt2ePj4uKi/iqJqFCI2l1vNKko5wLVLGODyPhkDPvjCCasPoGE5FRlbfSqXw5bxrRA3fJ2iE1KxXvLj2PUX8cQncD6YESkxau9xo8fj+HDh6NVq1a4efOmDEGurq5qr46IikQ1Z2usHdEUw1pUlrW7/joYKjctPB4apayN8g4WWPF2E4x9tbpcgbbxxG257H7fJXXL7omIcmOgec6+5u7du+c43rx5Mzw9PVG2bNkcj69evRq6KiYmBra2toiOjoaNjU3+X3CS7ROe46oX0m77Lt7F2BXHERaTKPcG8mtdDcNfrioDiypi12m/Zcdw9V6CDFtv+VTG+22qw8zYSFkbRFT8xTzj7+/n7vkRL5r9Jub+uLu7P/Y4ERUPTas6wt/PBx08ysiJyt9tO4++C4KUTlT2KmeHTaN90K9hOYg/xxYEXEaXOXtx/o66ZfdERC/c86MP2PND9DjxUbHq6E18vu4U4pPTYG1mjK+61kZX75y9vvklJlh/tPqknG8kKsdPaOeGgU0qyo0TiYiKpOdn4sSJOHLkyPN+GRHpODFRuWc91xwTlf2WB2O0mKj8QN1E5Ta1Ssuepsxl919sOIOBCw/ijsJl90Sk3547/Ny4cUOu7hKTm8WE5y1btiA5WV1hRCLSbpkTld9rnTFRef3xW3K5+v7L95Qvu//qv2X3gRfuou2MAPifuq2sDSLSXy807JWeno69e/diw4YNWLduHW7fvo1XX31VLn3v2LEj7O3tocs47EX0bI5ev4/3lgfj2n8Tld95qYoMRWK4SpWL4bGyh+nUzRh53Lu+KyZ2qgUrM2NlbRBR8fCsv7+VzPk5e/ZsVhASQ2INGzZE586d5WToR1eB6QKGH6JnF5eUii83nMaKwzfkce2yNpjRxxtVnayUtSGGv3745zzm/3tJTogub2+BH/p4oV6FksraICLdV6jhJ7vw8HAZhET5Cx8fH3zwwQfQNQw/RM9vy8nbcqKymP9jbmKITzu4y6KpYq6QKgcu35PL7m9GPYCY/zzylWoY9UpVmBip62kiIt1VZOGnOGD4IXoxYdGJeH9lMPZezJj/07qmE6b0qANHKzNlbcQkpuDzdaex5tjNrGXyoheokqOlsjaISDcVSPgZO3bsM1/A9OnToasYfohenKgG/9veK5jqH4LktHQ4WpliWk9PtHRzUtqOmGj96ZqTiElMhYWpESZ2dEefBuWU9jQRkW4pkPDTsmXLZzpPfPjs3LkTuorhhyj/ztyKgd/yYzh/J04ev9GkAj5uXxPmJup2bb4V9QBjVwRj/+VIefyquzOmdPeAg8KeJiLSHRz2ygeGHyI1ElPSMGXLOSzad1Uei0nQM/t6oZaLrdKepl/2XMa0rSFISdPIIbZpveqgZQ21PU1EpP0KJfxERUXh119/lau9hFq1auHNN9/U+fIWDD9Eav17PgIfrDyOiNgkmBgZ4IM2NWT9LpW7Np++FQ2/ZcG4EP6wp2lCu5ooYcr6YET6IqagdnjOdPjwYVSpUgU//PADIiMj5U3M8xGPHT169EVfloiKoZeql4L/GB85LCV6ZyZvOYfXfz2A29EPlLUhepM2jGqOQU0ryuPfg66h05w9OHWTf2AQkaKeH7GMvWrVqvj5559hbJyx2VhqaiqGDh2Ky5cvIyAgALqKPT9EBUN83Cw7FIovN5zBg5Q02JYwwbfdPNChTpkC7Wka+2oNDGtRWWkleiLSw2GvEiVK4NixY3Bzc8vx+JkzZ1C/fn0kJKir+FzYGH6ICtbliDi5a/OJGxk//z3qumJSZ3dYm5soa0MURp2w+gS2nr4jjxtWssf03p5wLWmhrA0i0rNhL/Gi169ff+zx0NBQWFtbv+jLEpEeqFzKCquGN8XIllXlZoWrjt5A+1mBOHItY9WWCvaWppj/ej1M7VFHLoU/eCUS7WYGYl1wxv5ARKS/Xjj89OnTB0OGDMHy5ctl4BG3ZcuWyWEvUdaCiOhJxK7MH7StgWXDmqCsXQmERj5Ar/lBmL79PFLT0pW0Ibbd6N2gHLaM8YG3qESfmIoxy9RXoici3fLCw16ikvu4ceMwf/58OddHMDExkZXep0yZAjMz3d1ng8NeRIUrt12bZ/TxQkWFuzaLQDV31yXM2nkBaekauNia4/veXmhSxUFZG0SkJ/v8iLk9ly5dkv8tVnpZWOj+eDrDD1HRELs2f7LmpOyhEUNVkzrVQq/6rkp3bX60Ev0wn8oY26Y6zIy5JJ5I13GTw3xg+CEqOqJo6djlwThwJWP+j2+t0pjc3QMlLU2VtRGflIqvNp6RK88E9zI2cvPFas6cr0ikywol/CQmJuLEiROyknt6es4x+s6dO0NXMfwQFS0xLLUg4DK+3xaC1HQNnG3M8H0vLzSv5qi0na2nw/DRqhO4n5ACM2NDTGjnhoFNK7I+GJGOKvDw4+/vjzfeeAN37959/EUNDJCWlgZdxfBDpB1O3ojGmOXHcDkiXh4PaV4J49rWUFofLDwmEeP+PiH3BhJaVC+F73rWgZONubI2iKiYLHUfNWoUevXqhdu3b8ten+y3Fw0+kydPRoMGDeRSeScnJ3Tt2hUhISGP9TaNGDECDg4OsLKyQo8ePXDnTsY+HpnEEvwOHTrI+UfidcTE7MxJ2USkOzxcbbFplA/6Nyovj3/dcwVd5+5FSFissjZEyFk0uAG+6FxL9v4EnI9A2xkB8D8VpqwNItIuLxx+ROAYO3YsnJ2dlV3Mv//+K4PN/v37sX37dqSkpKBNmzaIj8/4q0947733sGHDBqxcuVKef+vWLXTv3j3reRG8RPARq9H27duHxYsXY9GiRZg4caKy6ySiwiNqc33TzQO/vFEfDpamOBcWK8tW/LbniixqqoLorRbDXRtHNZfzf8Qw2Dt/HsGHf5+Q84OIqHh54WEvUcC0WbNmcq+fghIRESF7bkTIadGihezGKlWqFJYuXYqePXvKc86dO4eaNWsiKCgIjRs3xpYtW9CxY0cZijKDmViO/+GHH8rXMzV9+qRJDnsRaafw2ESM//sEdocU3BBVcmo6vt8eIucciU/HCg4W+KGPF+qWL6msDSLS0Tk/Yom7GPYSYcTDw0Pu8ZPd6NGjkV8XL15EtWrVcPLkSdSuXRs7d+5Eq1atcP/+fdjZ2WWdV6FCBfj5+cleIdHDs379egQHB2c9f+XKFVSuXFkWXPX29n6snaSkJHnL/uaVK1eO4YdIC4mPrD/2X8M3m84iKTUdJS1MMKVHHbStVVppO0GX7uH9FcG4FZ0oa4KNeqWq3JHa2OiFO8yJSEvCT0ZF0hfw119/Ydu2bTA3N8fu3btzrI4Q/53f8CPmDolAI3qXRPARwsLCZM9N9uAjiB4e8VzmOY8OxWUeZ56T21yjL774Il/XS0SFQ3y+vNGkIppUdsDoZcE4ezsGb/9xBP0alsNnHd1hYfrCH2s5iM0Pt/i1wMR1p7Au+BZm/HNB9jip3nyRiArfC/8J88knn8jAINLV1atXZe9K5k1Udc8vMffn1KlTsmRGQZswYYL8PjJvolQHEWk3sSfP2hFNZbV28bfXXwdD0WHWHhwPjVLWhqg6P7Ovt9wDyNrcGMGhUbIG2fJD12UPFBHpWfgRE4pFfS9DQ/VdwCNHjsTGjRuxa9cuuLq6Zj1eunRp2W5UVNRjk6/Fc5nnPLr6K/M485xHiVIconss+42ItJ/Ylfnj9jWxZEgjlLYxx5W78egxbx/m/FfCQpUuXmVlfbBGleyRkJyGD1edlBOiReV4ItI9L5xcBg4cKIuaqiT+khLBZ82aNXJ+T6VKlXI8X69ePTm3aMeOHVmPiaXwYml7kyZN5LG4F3OExMaLmcTKMRFo3N3dlV4vEWmHplUd4e/ngw4eZeSmiN9tO4++C4IQGpmgrA3XkhZY+lZjfNTODSZGBth6+o5cEp+5PxAR6Y4XnvAs5vT8/vvv8PT0RJ06dR6b8Dx9+vTnfs13331XruRat24datSokfW4mLxUokQJ+d+icOrmzZvl8nURaMR+Q4JY1p651N3LywsuLi6YOnWqnOczYMAAWW3+22+/fabr4GovIt0kPs5WHb2Jz9edQnxyGqzNjPFV19ro6l1WaTunbkbDb3kwLobHyeNBTSvKUKRy80Ui0sLVXi1btsz7RQ0MZM/N88prS/mFCxdi0KBBWZscvv/++3LCtVih1bZtW/z44485hrSuXbsmQ5KYiG1paSl7qUSleWPjZ5sIyfBDpNuu30uA3/JjOHo9Y4i8s6eLDEFiDo8qiSlpmLz5LBYHXZPH1ZysMKOvF2q5POHfOxEVKBY2zQeGHyLdl5qWjrm7LmHWf/N/ytqVwPe9PdG4soPSdnaFhMu9hyJik+Rw2AdtamCoT2W5PJ6Iikl5C7GPzpEjR/J7fUREBUrsxzOmdTWsfKeJ3KhQVIvv9/N+/M//nNzIUJWWNZzgP8YHr7o7IyVNg8lbzqH/L/tle0SknZ47/Ny4cQPt2rWTq7DE0JLYUVmswCIi0kZiZ+ZNo33Qu76r3LF53u5L6D5vb9Z8HRUcrMywYEA9/K+HByxMjbD/ciR8ZwRg/fFbytogoiIMP7/99pucRCzm3IgCpGIjQkdHR1lgVEyAjoyMVHh5RET5Z2VmjKk9PTH/9bqwszDBqZsx6Dg7EH/uv6Zsvx4xZ7FPg/LYPNoHXuXsEJuYitF/HYPfsmOIfpCipA0iUkPJnJ+zZ8/KYqNilZYYEmvYsCE6d+6Mfv36oWxZtassCgPn/BAVX2HRifhg5XHsuXhXHreu6STLYzhamSlrIyUtHXN2XsTsnRcgthsqqPlGRKQlE57F/joiCIn6Wj4+Pvjggw+gaxh+iIo3UQ3+t71XMNU/BMlp6XC0MsW0np5o6eaktJ0j1+7jveXBuB6ZIHehfrtFFYx9tTpMjVkfjKggcLVXPjD8EOkHURfMb1kwQu7EyuM3mlSQO0ar3K8nLikVX204g+WHM8rm1HKxkeUyqjpZK2uDiAo4/HTv3v2p54j9dMS+O6+++io6deoEXcPwQ6Q/xH49YgXYwr1X5XFVJysZTlTv1+N/6jY+Wn0SUQkpMDM2xCcdamJA4wp57m9GRFq01F286NNuYjfmCxcuyNpfYmk8EZG2Er08n3eqhcVvNkQpazO5Cqzr3L346d9LcnhMFd/aZbDVrwV8qjkiKTUdE9edxuBFhxAem6isDSJ6NgU67CWKk4qSFaL2li5hzw+RfroXlyR7Z7afySiG3KSyg5yo7GKXUV5HBRGofg+6im+3ZOw3ZG9piindPdCmVu6Fl4lIC3p+nkfz5s1Rv379gmyCiEj5fj2Tu3ughIkRgi7fk/v1bDyhbr8eQ0MDDGpWCRtHNUfNMjayMvywP47go1UnEJ+UqqwdIsobJzzngj0/RHQ5Ik4WLz1xI+PfaPe6ZfFF51qwNldXHywpNQ3Tt53HgsDLcgPGig4W+KGPF7zLl1TWBpE+idGGnh8iIl1VuZQVVg1vipEtq0KU6Vp99CbazwrEkWvqNnI1MzbChPY1sWRoI5SxNcfVewnoOT8IM/+5IGuTEVHBYPghIsqDiZEhPmhbA8uGNZEbFYZGPkCv+UGYvv280nDStIoj/Me0QCdPF1mE9Yd/zqP3T0G4di9eWRtE9BALmxIRPUXDSvbY4ueDbt5l5Y7Ns3ZckD00V++qCye2FiaY1dcLM/p4wdrMGEevR6H9zECsOByqrAQHEWVgYVMiomdgY24i5+PM6ucNa3NjBIdGyWGwFYfUhROx509X77IyaInAFZ+chvF/n8DwP4/ifjw/Z4lUYWFTIqLn0NnTBf5+LdCokj0SRDhZpT6cuJa0wF9vNcZ43xowNjSA/+kwtJ0RgIDzEcraINJnLGyaC672IqKnEXNzFgRcxvfbQpCaroGzjRm+7+WF5tUclbZz6mY0xiw7hksRGUNsg5tVxIe+bkpLcBAVF0VW2ysiIkIWNWVh02wYfoiKrZM3ojFm+TFc/i+cDGleCePa1lAaTh4kp2HylrP4PeiaPK7ubIUZfbzh7qLg84moGGFh03xg+CGi5w0nX286gyUHMnazdyttjZl9vVGjtNripbvOhWPc3ydwNy4JpnIlWnUMbV5ZbpxIROA+P0REhaWEqRG+6eaBX96oDwdLU5wLi0WnOXvw254rSuuDtXRzgr+fD1rXdEZyWjq+3XwO/X85gFtRD5S1QaQPGH6IiBRp7e4sJ0O3rFFK1u36cuMZDBLFS2PUFS91tDLDz288XoJjw3F1JTiIijuGHyIihURl+N8GNcCXXWrBzNhQrtASK7W2ng5T1oZYEt+vYXlsHuMDT1dbxCSmYtRfx/De8mDEJKYoa4eouHrh8PPmm29i0aJFWcfXrl2Te/6IcTYiIn0mwskbTSrK4qXuZWxwPyEFb/9xBBNWn0BCsrripZUcLfH38KYY/UpGCY41x26i3YxAHLzCLUeICiT8bN68GW5ubvK/o6KiUK9ePXTt2hXu7u4ICQl50ZclIio2qjlbY82Ipni7RWUYGAB/HQxFh1l7cDw0SmkJjrFtamDlO01Qzr4EbkY9QJ8FQZjqf04OvRGRwvAjengy9/BZtWoVSpcuLWdZ9+nTBxMmTHjRlyUiKlayipcOaYTSNua4cjcePebtw5ydF+ReQarUq2CPzaN90Kueq6wQ/+PuS+g+by8uhscpa4MI+h5+ypUrhytXrsj/XrlyJQYNGgQzMzO888472Lt3r8prJCLSeU2rOsqVWh08yshNEb/bdh59FwQhNDJBWRvW5iaY1ssTP/avCzsLE5y6GYOOswPxx/5rrA9GpCL8iLAzevRofPbZZ9ixY4cc8hLS09MRF8e/NIiIHmVnYYo5r3nju16esDQ1wqGr92Xx0jXHbigNJ+09ysgq8c2rOiIxJR2frT2FNxcdQkRskrI2iPQy/IihrV69eiEgIABTpkxB1apV5eOHDh1C+fLlVV4jEVGxmgzds54rtoxpgbrl7RCblIr3lh/H6GXBiH6gbqVWaVtz/P5mQ0zs6A5TY0PsComQS+L/OXNHWRtEukr5Ds/Tpk1DYmKi7BHSVYW6w7N8nivkiPRRalo65u66hFn/zf9xsTXH9D5eaFzZQWk7IWGxsj6Y2HxReK1ReXzaoSYsTI2VtkNU1FjeIh8YfoioMB29fl/u0XPtXoJcFfbOS1XwXuvqssdGlaTUNHy3NQQ/B17JWiY/o48XPMvZKWuDqKixvAURkY6oW74kNo32Qe/6GSu15hXASi2x6uyTDu5YMvThqrPu8/Zh9o4LsgeKSJ8w/BARaQErM2NM7emJ+a/nXKn1p+KVWs0yV53VKSOH2r7ffh59FuzH9XvqVp0RaTuGHyIiLeJbO+dKrU/XnsJbvx+WldyVrjrr543pvT1l6Dpy7T7azwrE30fUrjoj0lYMP0REWiZzpZaYlGxqZIh/zobLlVq7zoUrXXXWva5YdeaDBhVLIi4pFR+sPI4RS4/ifnyysnaItBHDDxGRFjI0NMBQn8pYN7IZajhb425cMgYvOoSJ604hMSVNWTvl7C2wbFgTjGtbA8aGBth8Mgy+MwMQeCFCWRtE2obhh4hIi9UsYyMD0OBmFeXx70HX0HH2Hpy+pW6VqJGhAUa0rIrV7zZFZUdL3IlJwoBfD+LLDWeUBi0ibcHwQ0Sk5cxNjPB5p1pY/GZDlLI2k6vAus7di5/+vYR0hfXB6rjaYePo5ni9ccZGtb/tvYIuc/bi7O0YZW0QaQOGHyIiHfFS9VLY6tcCbdydkZKmweQt59D/lwO4FfVAWRti48Ovu3rgt0H14WhlipA7sTIA/RJ4WWnQIipKDD9ERDrE3tIUPw2ohyndPVDCxAhBl+/JydAbT9xS2s4rbs7w92uBVm5OSE5Lx9ebzmLAbwdwO1pd0CIqKgw/REQ6RqzU6tuwPDaP8YGnqy1iElMxcukxjF0RjNhEdfXBHK3M8MvA+vimW22Ymxhi70URtAKx6cRtZW0QFQWGHyIiHSVKVPw9vClGtqwKQwNg9dGbcr+eI9cilQat/o0qyB2o67jayuKrYjm86qBFVJgYfoiIdJiJkSE+aFtDLlcva1cCoZEP0Gt+EKZvC0GKwrIVVUpZYdUjQavdzEAcuqouaBEVFoYfIqJioGEle2zx80E377IQ85Jn7bwoQ9DVu/HKg9byt5vAtWQJ3Lj/AH1+CpIFU1UGLaKCxvBDRFRM2Jib4Ic+XpjVzxvW5sYIDo2Sw2DLD11XWraiQUV7uTN0j7quMmjN2XURPebtw6UIdYVYiQoSww8RUTHT2dNFrtRqVMkeCclp+HDVSQz/U23ZCmtzE3zf2xNzX6sL2xImOHEjGh1mqS/ESlQQGH6IiIohMf9n6VuN8aGvmyxb4X+6YMpWiOrwokp8s6oOWYVYhy5WW4iVSDWGHyKiYkqUrRj+chWsebcZKpd6WLbiq41qy1aUsS2BP95slFWIdce5jEKsO87eUdYGkUoMP4XBvnpRXwER6TEPV1tsGuWD/o0yylb8uueKLI8REharvBDr+lEPC7EOWXwYn6w5iQfJrA9G2sVAw8HZx8TExMDW1hbR0dGwsbFR86KTbJ/yvLoihUREefnnzB18uOoE7sUnw9TYEB/5umFQ04oyvKgiepWmbQ2RIUsQvU4z+njJ2mFE2vD7mz0/RER6pLV7RtmKljVKITk1HV9uPINBiw4hPCZRaSHWzzq6488hjeBsY4bLEfHo/uM+zN11EWmsD0ZagOGHiEjPiMrwvw1qgC+71IKZsSECzkeg7YwAbD0dprSd5tUcZSHW9h6lkZqukb1BYl+g0MgEpe0QPS+GHyIiPSTKVrzRpCI2jmoO9zI2uJ+Qgrf/OIIJq08gITlVWTt2FqZyOfz3vTxhZWaMw9fuy52hVx+9wSXxVGQYfoiI9Fg1Z2usGdEUb7eoDAMD4K+Doegwaw+Oh0YpDVo96rnKjRHrVyiJuKRUjF1xHCP/OoaoBHV7DxE9K4YfIiI9Z2ZshAnta2LJ0EYoY2uOK3fj5Y7Nc3ZeUDpHp5y9BZYNa4wP2lSXew+J6vCiSvzei3eVtUH0LBh+iIhIalrFEf5jWsiNC8Ucne+2nUffBWrn6BgbGWLkK9VkkdTKjpYIi0lE/18O4JtNZ5CUyiXxVDgYfrTFDK+ivgIiIthamGBOP++sOTqHrt5H+5mBWHNM7Rwdz3J22Di6edbeQz8HXkGXOXtxLixGWRtEeWH40RZRGfthEBEVtexzdOpVKInYpFS8t/w4Ri8LRvSDFGXtWJga45tuHvjljfpwsDTFubBYdJ6zV+4PlM4l8VSAGH6IiCjPOTrLhzXG2Fery1IZG47fQrsZAdh/+V6B7D30ipuT3HtIlN9447eDCItWt/cQUXYMP0RE9MQ5OqNbVcPKd5qggoMFbkUnot/P+zFlyzkZVFTuPfTrwPr4umttmJsYYs/Fu3Lvoc0nbytrgygTww8RET1V3fIlsWm0D3rXd4WY+jP/30voPm8vLobHKR1ue71xBWwc5YPaZW3kENu7S47ig5XHEZuobriNiOFHm0xyKOorICLKk5gAPbWnJ+a/Xhd2FiY4dTMGHWcH4o/915ROhq7qZIXVw5thRMsqcu+hv4/cQPtZgTh8NVJZG6TfGH60irpdVYmICopv7TJySXzzqo5ITEnHZ2tPYejiw7gbl6SsDVF0dVxbNywf1gRl7UogNPIBev8UhO+3hSAlTd1wG+knhh8iInpupW3N8fubDfFph5owNTLEjnPh8J0RgF3nwpW207CSPbb4+aC7d1mIBWCzd15Ez3n7cDlC3XAb6R+GHyIieiGGhgYY6lMZ60Y2Qw1na9yNS8bgRYcwcd0pJKao27DQxtwE0/t4YXY/b9iYG+P4jWhZgmPpgeusD0YvhOFH20yyLeorICJ6LjXL2MgANLhZRXn8e9A1dJy9B6dvRSttp5OnC7a+1wJNqzjgQUoaPl5zEm/9rna4jfQDww8REeWbuYkRPu9UC4vfbCiXrYtVYF3n7sVP/15SumFhGdsS+HNII3zSPmO47Z+zBTPcRsUbww8RESnzUvVS2OrXAm3cnZGSpsHkLedk7a5bUQ+UDre91aIy1o5ohurOVlnDbWLi9YNk1gejpzPQcMD0MTExMbC1tUV0dDRsbGyKZkhrktruYiKiwiR+tSw/FIovNpyRQ1Rirs633T3QsY6L0nbE3KKp/iH4bW9GiaDKpSwxs483PFw5hUAfxTzj72/2/BARkXJiw8K+Dctj8xgfeLraIiYxFSOXHsPYFcFKNywUw20TO7njjyEN4WRthssR8ej2417M3XURaawPRroQfgICAtCpUye4uLjIfzhr167N8fygQYPk49lvvr6+Oc6JjIxE//79ZeKzs7PDkCFDEBeng0siOfGZiIqBSo6W+Ht4U4xsWRWGBsDqozflhoVHrqndsNCnWsZwm2+t0khN12Da1hD0W7AfoZEJStuh4kGrwk98fDw8PT0xd+7cPM8RYef27dtZt7/++ivH8yL4nD59Gtu3b8fGjRtloBo2bFghXD0REeXGxMgQH7StgeVvP9ywsNf8IExXvGFhSUtTzHu9Lqb1rANLUyMcvBqJ9jMDsebYDS6JJ92Y8yN6ddasWYOuXbvm6PmJiop6rEco09mzZ+Hu7o5Dhw6hfv368jF/f3+0b98eN27ckD1KOjPnRzAwBT6PUNc+EVERi0lMwefrTmPNsZvy2KucHWb08UJFR0ul7Vy/lwC/5cdw9HqUPO5Ypwy+6eoBWwsTpe2Qdim2c352794NJycn1KhRA8OHD8e9e/eyngsKCpJDXZnBR2jdujUMDQ1x4MCBPF8zKSlJvmHZb1pBk1zUV0BEpJTYsPCHPl6Y1c8b1ubGCA6NksNgyw+p3bCwvIMFVrzdBGNfrQ4jQwNsPHEbvjMDsO/SXWVtkO7SqfAjhrx+//137NixA//73//w77//ol27dkhLy1jaGBYWJoNRdsbGxrC3t5fP5WXy5MkyKWbeypUrB63BuT9EVAx19nSBv18LNKpkj4TkNHy46iSG/3kU9+PV/dFnbGSI0a2qYdXwpqjoYIHb0Yly2f23m88iKZVL4vWZToWfvn37onPnzvDw8JDDYWJOjxjiEr1B+TFhwgTZRZZ5Cw0NhVaZXruor4CISDkx/2fpW43xoa8bjA0N4H86TPbOBF5QO9wvhtY2jfZBv4blIDqXFgRcRpc5e3H+TqzSdkh36FT4eVTlypXh6OiIixcvyuPSpUsjPDznLp+pqalyBZh4Li9mZmZybDD7rUC86N49MVoWxoiIFBFDUsNfroI17zaTe/TciUnCgF8P4quNZ5TWB7M0M8bk7nWwYEA92Fua4lxYrCzBsXDvFaU7UJNu0OnwIyYxizk/ZcqUkcdNmjSRE6KPHDmSdc7OnTuRnp6ORo0aQadx+IuIijGxKeGmUT7o36i8PP51zxVZHiMkTG3vTJtapeHv54OXa5RCcmq63IRx4MKDuBOTqLQd0m5atdpL7MeT2Yvj7e2N6dOno2XLlnLOjrh98cUX6NGjh+zFuXTpEsaPH4/Y2FicPHlS9t4IYg7QnTt3MH/+fKSkpGDw4MFyAvTSpUuf+ToKbLWXihDDnZ+JqJj758wdfLjqBO7FJ8PU2BAf+bphUNOKsqyFKuJX35/7r+HrTWL+TzrsLEwwpbsHfGtn/DFNuulZf39rVfgRc3dE2HnUwIEDMW/ePDnP59ixY7J3Ryxbb9OmDb766is4OztnnSuGuEaOHIkNGzbIVV4iLM2aNQtWVlbFI/zI12AAIqLiLSI2CeP/Po5dIRnzf3yqOeL7Xp5wsjFX2s7F8Fj4LQ/GqZsZq3x713fFxE61YGVmrLQdKhw6GX60RYGGn6+cgTQF3asMQERUzIlfT3/sv4Zv/uudKSl6Z3rUQdtaec/hfBFi+OuHf85j/r+X5ITo8vYWcjl+vQollbZDBY/hR1vDj8r5OwxARKQHLtyJxZhlwThzO6N3pm+Dcviso7ucxKzSgcv3MHbFcdyMeiBLcYx8pRpGvVJV7lBNuqHYbnJI2XASNBHpgWrO1lgzoineblEZBgbAskOhcqXW8dCM3ZtVaVTZAVv8fNDNuyzEArBZOy7IMhxX7sYrbYeKHsOPrmMAIiI9YGZshAnta2LJ0EYoY2suA0mPefswZ+cFpdXbs+9AbfPfDtQdZgVi2UG1O1BT0eKwV1EMe0XfBH5wV/uaVmWAD86pfU0iIi0UnZCCj9eexKYTt+Vxg4olMb23F8rZWyht51bUA4xdEYz9lzMq0L/q7ixXhDlYZawuJu3DYS9tZltW/WvG3c7oBbrxcI8jIqLiSBQnndPPW67+EquyDl29XyDV213EDtRDG+Pj9m4wMTLA9jN30HZGIHaF5NxMl3QPe36KouenwIerDIBJasfCiYi0UWikqN4ejCPX7svjTp4u+LprbdiWUFu9/fStaPgtC8aF8Dh5/EaTCpjQriZKmBopbYfyhz0/2q5AV2ppMsLVzm8KsA0ioqInhrqWD2ucVb19w/FbaDcjAPsv31PaTi0XW2wY1Vxutij8HnQNnebswambXHWri9jzU1Q9P4U5WbnRcKDdlMJpi4ioiBy7fl/2Al27lyBXhb3doooMRWKXaJX+PR+BD1YelxsxiuGwsa/WwLAWlWX4oqLFfX50IfzMagBEnkeh8RoAdJ1TeO0RERWy+KRUfLnhDJYfzigIXbusDWb08UZVp2ff5f9ZRMYnY8LqE9h6+o48bljJHtN7e8K1pNpJ1/R8GH50IfwU1VJ19+5A74WF3y4RUSHxP3UbH60+iaiEFJibGOKTDu54vVF5GIguIUXEr8+Vh29g0obTSEhOg7W5sZxv1MWrABa10DNh+NGZ8OMAIBVFopov0H950bRNRFTARKX291ccx56Ld+VxKzcn/K9nHTgqXqp+7V68HG47dj1joUlnTxd8VQCTrunpGH50Jfxoy0aFLJVBRMVQeroGv+29gqn+IUhOS4ejlSmm9fRESzcnpe2kpqVjzq6LmL3zotx00cXWHN/39kKTKuIPXCosDD+6FH60JQAJ7A0iomLo7O0YuVQ95E6sPB7QuAI+bq9+qfrR6/fxXrZJ18N8KmNsm+pyh2oqeAw/+aDX4Sc7zg0iomIkMSUN//M/h4V7r8pjMQl6Rh8v1C5rq3zS9Vcbz8gaZIJ7GRvM7Osla5RRwWL40bXwo60BSOAqMSIqRh5dqv5+mxqyh8ZQ8VL1rafD8NGqE7ifkAIzY0NMaOeGgU0rKp10TTkx/Ohi+AnxB/7qA61mUw6o2ZH7BhGRThNL1UUw2XYmY6l648piqbqXLGmhUnhMIj74+wQCzkfI4xbVS+G7nnXgZGOutB3KwPCji+FHmFRSTNGDzmgxHnjlk6K+CiKi5yZ+/S0/FIovNpzBg5Q0WcX92+4e6FjHRXk7YkfobzefRVJqOkpamGBy9zrwrV1aaTsEhh+dDT/aPPz1RCZAvz+BGr5FfSFERM/lyt14+C07huM3Mla9dq9bFl90rgVrc7VL1S/cicWYZcE4cztGHvepXw4TO7nD0sxYaTv6jOFHl8OPzgagbExsgJH7C6aCPRGRYilp6Zj5zwX8uPsi0jWiZlgJORm6XgV7pe0kpaZh+vbzWBBwGeK3bwUHC/zQxwt1y4tef8ovhh9dDz/FIQA9Oldo7Kmivgoioic6dDVSLom/GfUAYv7zyJZVMapVNZgYqa0PFnTpHt5fEYxb0YmyJtioV6rKtowVt6NvYhh+ikH4KW4BKLvS3sA7u4v6KoiIHhOTmILP153GmmM35bFXOTvZC1TR0VJpO9EPUvDZ2lNYf/xWgbajT2IYfopJ+CnOASg7M1ugSivuK0REWkOEkk/WnERsYiosTI3weSd39K5fTvlS9XXBN/Hp2lMF3o4+iGH4KUbhR18CUA4mQPXWQM1OgHf/or4YItJTYvhr7PJgHLgSKY99a5XG5O4eKGlpqrSdG/cTZB2yzHba1nKWK8LsFbdT3MUw/BSz8KOXAegRFk7A+AtFfRVEpGdErS4xQXn69hCkpGngbGOG73p5wqdaKeXt/Bx4Gd9vy2inlHVGOy9VV9tOccbwUxzDjzC1GpAQXtRXoT04d4iICsmpm9EYvewYLkfEy+MhzSthXNsaMDcxUt6OqBJ/MTxOHg9qWhEftXNT3k5xxPBTXMNPJn3vBXoSRzfAvTM3XyQi5R4kp+GbzWfw5/7r8tittDVm9vVGjdLWyuuQTd58FouDrsnjaqIOWV8v1HLhZ/+TMPwU9/AjMAA9O3MH4KPLRX0VRFRM7Dh7B+P/PoF78ckwNTbER75usodGdX2wXSHhGLfyBO7GZdQh+6BNDQz1qSyXx9PjGH70IfwIv7QFbuwv6qvQMeJD478fe0MTwKUeMHRrUV8UEekYURh1/N/HsSsko26XTzVHfN/LU3ndrntxSfho9Ulsz1aH7PveXiiruA5ZccDwoy/hJxN7gdQxtgB6LQQubAOqtWHJDiLKk/gV+sf+a/hm08O6XVN61EHbWqWVt7PicEYdsoTkNFibG+Obbh7o7Km2DpmuY/jRt/AjrBgMnFld1FdRfJWqBUScBUrVBEbsK+qrISIt8mjdrr4NyuGzjurrdsk6ZMuDcTw0Sh539XLBF11qw7aE2jpkuorhRx/DT6a1I4HgP4r6KvSDoRmQngTACKjehvsSEekxWbdr23ksCMyo21XRwQIz+nrLnZtV1yGbs/MiZu+8IOuQieGv73t7onFlB+i7GIYfPQ4/mTgUVvTsqwORFwD7asDoQ0V9NURUCPZduis3LLz9X90uv1bV8G7LqsonKR+5dh/vLQ/G9cgEiM2g325RBWNfrS4nYOurGIafF1dswk+mSeKvDv7frFWMzIG0xIz73os5v4iomIlOSMHHa09i04nb8rh+hZKyens5ewul7cQlpeLLDaex4vANeVzLxQYz+3qhqpPapfe6guEnH4pd+MnEEKRbXBsDt45wNRqRjhK/XlcfvYnP15+WIcXKzBhfda2Frl5lldft8j91W64Ii0pIgZmxIT7pUBMDGlfQu/pgMQw/L67Yhp9MO78BAqYW9VXQizKzA5KiMu4nZGyARkTaKzQyQU5SFsNUQidPF3zdVf0k5Tsxifhg5XEEXrgrj1+uUQpTe9aBk7XapffajOEnH4p9+Mm05SPgwLyivgpSSqwsSc24n3QvI+ieWc8dr4mKWGpaOn7cfQkzd1yQNbxcbM0xvY+X8knK6ekaLA66islbziE5NV0WRp3S3QNtFC+911YMP/mgN+EnO64Q0y+TooH5LwN3TgDOdVgfjaiQHLt+X/YCXbtXsJOUz/+39P7sf0vv+zUsh087qF96r20YfvJBL8NPdlwlpt9EMJpSGUi8x7IgRAUgXk5SPoPlh0Plce2yNpjRxxtVnawKfOm9mHTtXb4kiiuGn3zQ+/CTiaUzKK9w9KVTxv5GYp+jieFAiD9XrBHlY5KyuYmYpOyO1xuVVz5J+dGl96NfqYYRLavA2Kj4LYln+MkHhp88zG0KRJwu6qsgbefaCLhx8L+VhQbApKiHQ6vn/YHqvkDXOUV9lURa4dFJyq3cnPC/nnXgaGWmfOn9J2tPYuN/S+/rlreTvUAVHCxRnDD85APDzzP4zg2Iy/hHRPRUBsaAJjXvY9GblNnbyOX9pGfEJOWF+67if2KSclo6HK1MMa2nJ1q6OSltR/y6Xxd8C5+tPYXYpFRYmhrh88610Kuea7FZEs/wkw8MP89pcoWMpddEBSkzIM3wAqKuAHaVAL/gor4qImXE5GS/ZcEIuRMrj8U+PR+3r4kSpkZK27lxPwFjlx/HwauR8ti3VmlM7u6Bkpam0HUMP/nA8JMPx5YA694t6qsgfQtF2Sfpi+Op1YCEcMDCCRh/oSivjui5JKak4X/+57Bw71V5LCZBz+jjhdpl1S5ESUvX4KeAS3JCdGq6Bk7WZviulydaVC8FXcbwkw8MPwqx0jxpa0DK9G05IDkGMLUBPs5YfUNU1ALOR+D9lccREZsEEyMDvN+mBob5VIah4vpgp25GY8yyY7gUES+PBzeriA993WBuora3qbAw/OQDw08BmtUgo9Any2yQtgWipwWlb8oCKXGAiRXwyc3Cu07SW5Hxyfho1QlsO3NHHjeubI/pvb3gYldCaTsPktPw7eaz+GN/xo7x1Z1Fb5M33F107/cfw08+MPwUMs4ZIl2UGYpyC0pflwFSEwBjC+BTLgygFyd+RS8/FIovNpzBg5Q02Jgb45tuHrJEhmq7zoVj3N/HcTcuGaZGhhjXtgaGNK+kvLepIDH85APDTxES5RjObgAizhb1lRA9megBEj1B+Q1KAsMSPcWVu/HwW3YMx29k/Nx09y6LL7rUgrW52vpgd+OS8NGqk/jnbEZvU5PKDvi+t6fy3qaCwvCTDww/2ri/kAhD6UV9JUQvRoQaEW6eZ+hN3jMs0UMpaemY+c8F/Lj7ItI1QDn7EvihtxfqV7RX2o5Go8GyQ6FyF+qC7m1SjeEnHxh+dEDmJFWi4kh1WOKk7mLl0NVIuST+ZtQDiBGpkS2rYlSrajBRvGPz5Yg4vLc8OEdv06QutWCjuLdJJYaffGD40dFVZZd3AomcO0R66EXD0pMCk8AtA7RWTGIKPl93GmuOZUy+9yxnh5l9vFDR0VJ5b9PsHRcwZ1dGb1NZuxJyZ+iGldT2NqnC8JMPDD/FSOaHNxE9lNkD9Ky9SI8+97TQlLmq074aMPqQssumx60/fgufrDmJ2MRUWIgdmzu5o3f9csp3bD5yLVJWow+NfCCr0Q9/qQr8WquvRp9fDD/5wPBTzG35CDi3CYi+XtRXQqSdRE9Pbn80ZPYAFVRoWtQFCN0PlGsMDFr34tevZ8Tw19jlwThwpWB3bI5NTJHzgFYeuSGPPcrayl4g1dXo84PhJx8YfvRU9i5+9hYRPT/76kDk+dwfFz1A+QlNj57z6HP75gBn1gHuXYCmI6FvxI7NCwIuY/r2EKSkaeBsk7Fjs0819Ts2bz55Gx+vKfhq9C+C4ScfGH4oi6hEfmErUK0tELyEK86ICkLFl4Gru3N/PLMH6EnhaGoVICGjKrpk4QiMv5Tz3C9KAZpkwMAU+DwCxZXYsXn0smO4/N+OzWKfHrFfj+odm8OiM6rR77mY8b63rFEKU3t6opS12mr0z4vhJx8YfuiZcYNGooKV16TtTG2+AbZ9kvvjmT1A+e1V0rHJ32LH5m82n8Gf+zOG9t1KW2NmX2/UKG2tvBr9on1XMcX/HJJT0+FgaYr/9aiD1u7OKCoMP/nA8EP58ktb4NYRwKUecGN/UV8NUfEOP0/72swen0dl7wF6nnlKuT2f23m5PS9WpV7aAVRpBfReiIK24+wdjP/7BO7FJ8uJyR/5umFQ04rKd2wOCYuV9cHOhWVUo3+tUXl82qEmLEyNUdgYfvKB4YcKRIg/cHE7UPVVYNmAh13wuX0wE5Ga8POkr32W5/NaMfpoD9BTA1TJR4bNDYFJ93P5GrGEPA2AETApYwJzfojCqOP/Po5dIRlBz6eaI77v5QknG3OolJSahu+2huDnwCvyuJKjpaxGL5bga+Pvb+1ao0ZUnNXwBTp8n3Ev/uIUH4yZ95k38YEqZN4TUdHKa/FD9sfzClCZj4sen8fmC6b/9/ij54vgI6Q9+XUzb3n5zk0+X+onT/w2qAG+6lILZsaGCLxwF21nBGDr6TCoZGZsJCc+LxnaCKVtzGU5jh7z9sk9glLTtG+uJMMPkTaRy4ijH95n3lwbA4YmGffZiYmdRKTdzqx++uOyxycXjz7+aODJtcfJFoj7r+xJ3G0YfGGHAU0qYtPo5qjlYoP7CSl4+48jsmL8g2cJUseWAEv7Ztw/RbOqjvD380GHOmWQmq7B99vPo8+C/bh+LyFrVVrQpXtYF3xT3ovjolD4A3JE9PyGbs37ObHE9+x6oGZnYNvEnF3mLzpcQESFLO3pjz+xJyj6YY9Pbr5zQ9UPzmHNu83w/fYQuSz+2+M+MDDM43UyzfQC7mcMZeH8FiBgGjAmOPc2/iujYmdqgzkTrqOVmxMmrjuNI9fuo/2sQPSoWxZbz9yRK8UylbE1lxsz+tYug8LEnh8iXSdWtAzZlnEvA0/0w7kC2XuPRAkEQdznNhmTiHRfZo9PHo+Lic8T2tXEGbPXZPDJPvVZ82jAEj09mcEnkzjOrQdIfF1mvcXkGNnb1L2uK7aM8UGDiiURl5SKxUHXZPA5i9dwyeQ1eS+Oh/95FP6nCrdAL8MPkb4Q1b9F6MmsAp7XXCMGI6Jir4RY8fTIYwaP9j+tezf3L370cdHjk5tvy6GcvQWWDG0MK7OMgabLpq/B3AwwMoK8v2T6mnz8iw1nCnUIjMNeRPT4viWPBqC5TYGIs0CpmsCIfRxOIyquNMDqozfQzbvsY+EoT5k9Pnk8Loa9RM+P6OkRm0BnbgSdeX8Gr6Fm9FIcvBKJJlUcUBgYfojo6UTgye5pNZke3e/kKycgLalwrpWIXlw6MHbFcew8F46ZGbMH8y08NmOOj6nJw8CTSRyLx5Hy8LzCwPBDRPn3aBHKR8PRZ+EPC8q6dQDaTXk8ID1aooCICp8RYJRmgI0nbmOmoq2AnKzNlZ6nAuf8EFHhEIHnvZMZ94/OORJELSZRkqBco4x7+Xjmn4kG/y35b5TzNR89JqJ8MQLw9ztNUMHBIu8FaM+pYSV7uarrScTz4rzCwp4fItIeYsVa9orckx6pmzZ0W86dssWGkY/uiuvZBzi+9OHXeL6W8zjrtXYCv7xSUN8Jkc7yLl8Sm0f7AJPVvJ6RoYFczo7leZ8jnhfnFRaGHyLSLSLwZIaeTI+WAWgw9OEcJNd6QLd5uZcNeDQY5RWUciU+qFkdiIony/9WZ+VG/NQ/b0wR+/ik/Z3Hk0YZzxcmhh8iKn5E4BG37HKrkyRC0bMEpTXDHw9J8jyueiP9ky62+4lLgqOV2XN9XV6Tp1VMqn5enPNDRPpNBJ4mI3KGpUc3ixRBRwyTtf02414cy/Ois310G+U9LymvvZPE3CYiXZMG+M4IwK5zedQ90wFaFX4CAgLQqVMnuLi4wMDAAGvXrs3xvChAP3HiRJQpUwYlSpRA69atceFCzv1JIiMj0b9/f1nN1c7ODkOGDEFcXFwhfydEpBchKbegJOYl9Vue0aMk7sWxPC+XCdxiftOj9dnEcV5hyWuA6u+K6PkZAXfjkjF40SF8tvaUqnnR+ht+4uPj4enpiblz5+b6/NSpUzFr1izMnz8fBw4cgKWlJdq2bYvExId7A4jgc/r0aWzfvh0bN26UgWrYsGGF+F0Qkd4Tc5I6fJ/L3KSo/4JStoncj65yE8fy3EcCkDjuOgcweGS2gjjOKyxVfPn5Hid6RoObVZT3f+y/Bl2kVXN+2rVrJ2+5Eb0+M2bMwKeffoouXbrIx37//Xc4OzvLHqK+ffvi7Nmz8Pf3x6FDh1C/fn15zuzZs9G+fXt89913skeJiEjrV7llyi3UfH4PWDsSuLAVqNY2IxBlnvvo3kny3jbvfZlyrQgeDcxqAESef/w5++q5P056xUiuzqqFljWc8P7K40BiEU3cKS49P09y5coVhIWFyaGuTLa2tmjUqBGCgoLksbgXQ12ZwUcQ5xsaGsqeorwkJSUhJiYmx42ISGuJwDPuwsPgk+nRvZMyHxM9PUbmGfePPvfo1wujD+Xerng8r14m8bh5HqUJ8nqcdFqL6qWw1a+FzgUfnQo/IvgIoqcnO3Gc+Zy4d3L6r0Djf4yNjWFvb591Tm4mT54sg1TmrVy5PIq0ERHpItHT89mdx3fiziswZT4uenrE/CRx/yyh6aPLubcvHn9SaDLMY9VQXo+T1rC3NNWdIJGNLl6zchMmTEB0dHTWLTQ0tKgviYio6Mmenqjce4KeFJoye3rE/bOEpol5rBrKfFyssMtNXo9ToTLI43ExETo2MQXaSGfCT+nSpeX9nTt3cjwujjOfE/fh4Tn/EaWmpsoVYJnn5MbMzEyuDst+IyKiF5TZ05NbT9CTQlNmT4+4z/68WGEn9lbKThyLx5/UoyS4d8/9+bweJ3XSgPazAnH4ai57bBUxnQk/lSpVkgFmx44dWY+JuTliLk+TJk3ksbiPiorCkSNHss7ZuXMn0tPT5dwgIiLSYqKnR4SW3HqC8tpr6Uk9SkLvhbn8qjP87/E8JpVnf9wmj2kQeT1ODxkBoZEP0PunIEzfFpLnkvg0fQ8/Yj+e4OBgecuc5Cz++/r163LfHz8/P3z99ddYv349Tp48iTfeeEOu4Oratas8v2bNmvD19cVbb72FgwcPYu/evRg5cqRcCcaVXkRExXSvpSf1KMnn7mf09JjbZdyL40e/Nq/jsadyv5bsjz8tQBmY5v58Xo8XI929yyJdA8zaeTHvlJMG+J+6rb/h5/Dhw/D29pY3YezYsfK/xcaGwvjx4zFq1Ci5b0+DBg1kWBJL283NH1aLXbJkCdzc3NCqVSu5xL158+ZYsGBBkX1PRESkBURPz0fXHvb4PFd4in7Y0yPu8zonr+PPI3JvM/vjee32rcO7gBsBmN7HCzP6eD21FtgXG84gTaQkfdzn5+WXX5b7+eRF9P58+eWX8pYXsbJr6dJnLUxIRET0DPLqAcourx6gzOe+KAVokjN6fB4NRGKfpz0/AAl3c+72nX3/p0f3csqtXbtKQNSVx88RjxcRZxvzp5YAvh2diINXItGkioP+9fwQEREVWyLwiKCSV09QXrt9P2sPk+CXMW3kMY8+PukpQ3VP7B95vn6T8NhEpeepwPBDRESkLURPz5Btue/4/SxDdJnPZ/b0iPtnCjrI5fheHl/3yONPWVHnZJ0xNSU5RVRryHmKOBaPZz+vMBhonjTOpKfEKjKx2aHY84fL3omISK9NEkNRqRk9PnkGopIA0rM9YJg1sVzM5Wn+v50Ii07EJdPXYGAgprFkBB9xq5K8FKVtzbHnw1dgZPi02UFqfn+z54eIiIjyJgKP7Gm694Rz8l5RJwLN553c5X+LoJOYBKSlQd6LY0E8n9/g8zzY85ML9vwQERGpJZazi1VdYnJzpjK25jL4+NYuU6i/v7VqtRcREREVT761y+BV99JyVZeY3Czm+DSsZF+oPT6ZGH6IiIioUIigU1jL2Z+Ec36IiIhIrzD8EBERkV5h+CEiIiK9wvBDREREeoXhh4iIiPQKww8RERHpFYYfIiIi0isMP0RERKRXGH6IiIhIr3CH51xkljsTNUKIiIhIN2T+3n5a2VKGn1zExsbK+3LlyhX1pRAREdEL/B4XBU7zwqruuUhPT8etW7dgbW0NAwMDpYlUBKrQ0FBWi38KvlfPh+/Xs+N79ez4Xj07vlfa8V6JSCOCj4uLCwwN857Zw56fXIg3zNXVtcBeX/yfzX8cz4bv1fPh+/Xs+F49O75Xz47vVdG/V0/q8cnECc9ERESkVxh+iIiISK8w/BQiMzMzfP755/Kenozv1fPh+/Xs+F49O75Xz47vlW69V5zwTERERHqFPT9ERESkVxh+iIiISK8w/BAREZFeYfghIiIivcLwU4jmzp2LihUrwtzcHI0aNcLBgwdRnAUEBKBTp05yp02xU/batWtzPC/m2k+cOBFlypRBiRIl0Lp1a1y4cCHHOZGRkejfv7/cCMvOzg5DhgxBXFxcjnNOnDgBHx8f+b6KXUOnTp0KXTN58mQ0aNBA7iru5OSErl27IiQkJMc5iYmJGDFiBBwcHGBlZYUePXrgzp07Oc65fv06OnToAAsLC/k648aNQ2pqao5zdu/ejbp168qVFlWrVsWiRYugS+bNm4c6depkbZDWpEkTbNmyJet5vk95mzJlivy36Ofnl/UY36+HJk2aJN+f7Dc3N7es5/le5XTz5k28/vrr8v0Qn+EeHh44fPiwbnzGi9VeVPCWLVumMTU11fz222+a06dPa9566y2NnZ2d5s6dO5riavPmzZpPPvlEs3r1arGiULNmzZocz0+ZMkVja2urWbt2reb48eOazp07aypVqqR58OBB1jm+vr4aT09Pzf79+zWBgYGaqlWravr165f1fHR0tMbZ2VnTv39/zalTpzR//fWXpkSJEpqffvpJo0vatm2rWbhwofwegoODNe3bt9eUL19eExcXl3XOO++8oylXrpxmx44dmsOHD2saN26sadq0adbzqampmtq1a2tat26tOXbsmHz/HR0dNRMmTMg65/LlyxoLCwvN2LFjNWfOnNHMnj1bY2RkpPH399foivXr12s2bdqkOX/+vCYkJETz8ccfa0xMTOR7J/B9yt3Bgwc1FStW1NSpU0czZsyYrMf5fj30+eefa2rVqqW5fft21i0iIiLreb5XD0VGRmoqVKigGTRokObAgQPy+9q6davm4sWLOvEZz/BTSBo2bKgZMWJE1nFaWprGxcVFM3nyZI0+eDT8pKena0qXLq2ZNm1a1mNRUVEaMzMz+cMtiA8G8XWHDh3KOmfLli0aAwMDzc2bN+Xxjz/+qClZsqQmKSkp65wPP/xQU6NGDY0uCw8Pl9/7v//+m/XeiF/wK1euzDrn7Nmz8pygoCB5LD5oDQ0NNWFhYVnnzJs3T2NjY5P1/owfP15+uGfXp08fGb50mfgZ+OWXX/g+5SE2NlZTrVo1zfbt2zUvvfRSVvjh+/V4+BG/iHPD9yon8TnbvHlzTV60/TOew16FIDk5GUeOHJFdftnrh4njoKAg6KMrV64gLCwsx3si6rGI4cDM90Tci27Q+vXrZ50jzhfv3YEDB7LOadGiBUxNTbPOadu2rRwyun//PnRVdHS0vLe3t5f34ucnJSUlx/sluuPLly+f4/0S3c7Ozs453gtRRPD06dNZ52R/jcxzdPXnMC0tDcuWLUN8fLwc/uL7lDsxVCOGYh79nvh+PU4My4ih+sqVK8vhGDGMJfC9ymn9+vXys7lXr15yeM/b2xs///yzznzGM/wUgrt378oP6ez/IARxLH449FHm9/2k90Tci39U2RkbG8tAkP2c3F4jexu6Jj09Xc7JaNasGWrXrp31vYh//OKD4knv19Pei7zOER/ODx48gK44efKknHMh5ky88847WLNmDdzd3fk+5UKEw6NHj8p5ZY/i+5WT+MUs5t/4+/vLuWXiF7iYayKqhPO9yuny5cvyPapWrRq2bt2K4cOHY/To0Vi8eLFOfMazqjuRFv6VfurUKezZs6eoL0Vr1ahRA8HBwbKH7O+//8bAgQPx77//FvVlaZ3Q0FCMGTMG27dvl5NF6cnatWuX9d9iUr0IQxUqVMCKFSvkhF3K+Uea6LH59ttv5bHo+RGfW/Pnz5f/HrUde34KgaOjI4yMjB5bFSCOS5cuDX2U+X0/6T0R9+Hh4TmeF6smxOqA7Ofk9hrZ29AlI0eOxMaNG7Fr1y64urpmPS6+FzF8GhUV9cT362nvRV7niJUWuvThLv4CF6tk6tWrJ3s0PD09MXPmTL5PjxBDNeLfkFhZJP6iFjcREmfNmiX/W/wFzfcrb6KXp3r16rh48SJ/th4hVnCJ3tbsatasmTVMqO2f8Qw/hfRBLT6kd+zYkSM1i2MxT0EfVapUSf7gZn9PRLevGOfNfE/EvfigER/gmXbu3CnfO/EXWeY5Ykm9GIvPJP7KFT0DJUuWhK4Qc8JF8BHDN+J7FO9PduLnx8TEJMf7Jca8xQdN9vdLDAdl/zAR74X4UM38kBLnZH+NzHN0/edQ/EwkJSXxfXpEq1at5Pcqeskyb+KvdTGXJfO/+X7lTSy5vnTpkvxFz5+tnMSw/KPbcZw/f172lOnEZ3y+pkvTcy11F7PcFy1aJGe4Dxs2TC51z74qoLgRK0zEck9xEz9q06dPl/997dq1rGWQ4j1Yt26d5sSJE5ouXbrkugzS29tbLqXcs2ePXLGSfRmkWD0glkEOGDBALoMU77NYRqprS92HDx8ul4Tu3r07xzLbhISEHMtsxfL3nTt3ymW2TZo0kbdHl9m2adNGLpcXS2dLlSqV6zLbcePGyZUqc+fO1bllth999JFcBXflyhX5cyOOxeqQbdu2yef5Pj1Z9tVeAt+vh95//335b1D8bO3du1cuWRdL1cXqS4HvVc6tE4yNjTXffPON5sKFC5olS5bI7+vPP//MOkebP+MZfgqR2M9B/MMR+/2Ipe9iX4PibNeuXTL0PHobOHBg1lLIzz77TP5gi2DYqlUruW9Ldvfu3ZP/EKysrORy0cGDB8tQlZ3YP0IsuRSvUbZsWfkPTtfk9j6Jm9j7J5P4wHj33Xflsk/xj79bt24yIGV39epVTbt27eQ+GOJDW3yYp6SkPPb/i5eXl/w5rFy5co42dMGbb74p9xcR1y9+sYifm8zgI/B9er7ww/cr55LzMmXKyO9BfJaI4+z71vC9ymnDhg0y7InPXjc3N82CBQtyPK/Nn/EG4n9evN+IiIiISLdwzg8RERHpFYYfIiIi0isMP0RERKRXGH6IiIhIrzD8EBERkV5h+CEiIiK9wvBDREREeoXhh4iIiPQKww8RERHpFYYfIqJCIoo4imKiXl5eqF27Nn7++eeiviQivcTyFkREhSQtLU1Wn7ewsEB8fLwMQIcPH4aDg0NRXxqRXmHPDxFRLl5++WX4+fkpfU0jIyMZfAQRgv4rLq20DSJ6OoYfItJ6gwcPxqefforicC1i6MvT0xOurq4YN24cHB0dlV4fET0dww8Raf1Q0caNG9G5c+dicS12dnY4fvw4rly5gqVLl+LOnTtKr5GIno7hh4gK1Z49e9CwYUOYm5vLXo+ZM2c+8fx9+/bBxMQEDRo0yHN4atSoUXKIqmTJknB2dpYTicWcGtFLY21tjapVq2LLli1ZXyOGnEaPHg0nJyd5Hc2bN8ehQ4eeeu3Zr+VF2s1OnC96gAIDA5/aLhGpxfBDRIVm8+bN6NatG959912cOHECb7/9Nt577z1cvXo1z69Zv349OnXqBAMDgzzPWbx4sQxSBw8elIFk+PDh6NWrF5o2bYqjR4+iTZs2GDBgABISEuT548ePx6pVq+TXiedFSGnbti0iIyOfeP2PXsvztit6eWJjY+V/R0dHIyAgADVq1Hih95KI8kGs9iIiKmgPHjzQuLq6apYsWZL1WGpqqsbKykqzePHiPL+uWrVqmo0bN+b5/EsvvaRp3rx5jte0tLTUDBgwIOux27dvi1nFmqCgIE1cXJzGxMQkx3UkJydrXFxcNFOnTs3xumPGjMnzWp63XeHAgQMaT09PTZ06dTQeHh6a+fPnP+VdI6KCYJyf4ERE9Kx27tyJBw8eoE+fPjlWP4leFDMzs1y/5uzZs7h16xZatWr1xNeuU6dOjtcUS8c9PDxyDDEJ4eHhuHTpElJSUtCsWbOs58VQlhiKE+3lJbdreZ52BdFGcHDwE78XIip4HPYiokKxa9cuubmfCAmZLl68KIeBvL298xxmevXVV+W8nCcR4SU7EaiyP5Y5TJWenv7C15/btRRGu0SkHsMPERWKY8eOITk5OcdjP/74I+rVq4fq1avn+jXr1q1Dly5dlF5HlSpVYGpqir1792Y9JnqCxIRnd3f3PL+uIK6FiIoGh72IqNDCj9jQ7/fff0ejRo2wcuVKzJs3T66gyo0YKhK7H4seF5UsLS3lxGSxx469vT3Kly+PqVOnyknJQ4YMKdRrIaKiwfBDRAXu+vXrciWV2CPno48+wvnz5+V8GX9//zyHvDZs2CDnyBTEJoBTpkyRQ1FiJZYYdhP1trZu3SqXrBf2tRBR4WNtLyIqcKLHROx9c+/evWf+GrGRoNh/RyxLL2radC1ElH+c80NEhTLklX0V1LMQYaNfv37QBtp0LUSUf+z5IaIC17VrVzm3ZtasWUV9KUREDD9ERESkXzjsRURERHqF4YeIiIj0CsMPERER6RWGHyIiItIrDD9ERESkVxh+iIiISK8w/BAREZFeYfghIiIivcLwQ0RERHqF4YeIiIj0CsMPERER6RWGHyIiItIrDD9EREQEffJ/XK/F4FMDxFEAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "for reg in rf.get_regions_rw():\n", " reg.add_pair(proppair=teqpflsh.PropertyPairs.DS, Nsplit=5)\n", " \n", " pset = reg.propset_bounding\n", " plt.plot(pset.rho, pset.s, 'o-')\n", "\n", " pset = reg.propset_Trhogrid\n", " plt.plot(pset.rho, pset.s, '.')\n", "plt.gca().set(xlabel=r'$\\rho$ / mol/m$^3$', ylabel='$s$ / J/mol/K');" ] }, { "cell_type": "markdown", "id": "d124dc3c-7f1c-49c8-9ab5-d5cc90d065f7", "metadata": {}, "source": [ "In order to get around this problem, the box needs to be sampled at more than 5 points\n", "\n", "Exercise for reader: build a more dense polygon defining the boundary of the box, with many more points along each side" ] }, { "cell_type": "markdown", "id": "cb086490-211a-4019-80c5-16305059ba94", "metadata": {}, "source": [ "To better understand the timing of each step, it can be useful to profile each step independently. The ``_many`` methods has been written for this purpose. The overhead in nanobind with pre-allocated buffers being passed to the function is functionally zero." ] }, { "cell_type": "code", "execution_count": 5, "id": "5fff9393-b345-4b60-a3d9-a8443c4bba68", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:49.895109Z", "iopub.status.busy": "2025-01-06T11:32:49.895026Z", "iopub.status.idle": "2025-01-06T11:32:50.541839Z", "shell.execute_reply": "2025-01-06T11:32:50.541574Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.19432550005149096 μs to look up a point from the K-D tree\n", "The K-D tree consumes 32.41525650024414 MiB\n", "0.1970805840101093 μs to look up a point from the K-D tree and return the (T, rho) point and its distance. This should be a smidge slower than the above calculation\n" ] } ], "source": [ "N = 500000\n", "\n", "for reg in rf.get_regions_rw():\n", " reg.add_pair(proppair=teqpflsh.PropertyPairs.DP, Nsplit=5)\n", " \n", " tree = reg.get_kdtree(teqpflsh.PropertyPairs.DP) \n", " X = np.linspace(2000, 2001, N)\n", " Y = np.linspace(0.25e7, 0.251e7, N)\n", " idx = np.zeros_like(X, dtype=int)\n", " d2 = np.zeros_like(Y)\n", " tic = timeit.default_timer()\n", " tree.get_nearest_indexd2_many(X, Y, idx, d2)\n", " toc = timeit.default_timer()\n", " print((toc-tic)/N*1e6, 'μs to look up a point from the K-D tree')\n", " print(f'The K-D tree consumes', tree.get_used_bytes()/1024**2, \"MiB\")\n", "\n", " TT = np.zeros_like(X)\n", " DD = np.zeros_like(X)\n", " tic = timeit.default_timer()\n", " reg.get_starting_Trho_many(teqpflsh.PropertyPairs.DP, X, Y, TT, DD, d2)\n", " toc = timeit.default_timer()\n", " print((toc-tic)/N*1e6, 'μs to look up a point from the K-D tree and return the (T, rho) point and its '\n", " 'distance. This should be a smidge slower than the above calculation')" ] }, { "cell_type": "markdown", "id": "5b5502e6-4283-4a21-9d8a-ba76e0cdfabe", "metadata": {}, "source": [ "Putting it all together, here is an example of using the entire flash calculation. The steps are:\n", "1. Find the nearest point in the K-D tree to get a starting value for $T$ and $\\rho$ for further iteration\n", "2. Do the iteration to find the right $T$, $\\rho$ satisfying the problem statement (trivial in this case because $T$, $\\rho$ are input variables)" ] }, { "cell_type": "code", "execution_count": 6, "id": "741a2939-dd1e-4e25-89d4-813dcb3ae040", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:50.543097Z", "iopub.status.busy": "2025-01-06T11:32:50.543014Z", "iopub.status.idle": "2025-01-06T11:32:50.937454Z", "shell.execute_reply": "2025-01-06T11:32:50.937239Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.3062583201099187 0.9456310399999999 0.25357912 0.0\n" ] } ], "source": [ "for reg in rf.get_regions_rw():\n", " reg.add_pair(proppair=teqpflsh.PropertyPairs.DT, Nsplit=5)\n", "N = 50000\n", "X = np.linspace(2000, 2001, N)\n", "Y = np.linspace(250, 251, N)\n", "TT = np.zeros_like(X)\n", "DD = np.zeros_like(X)\n", "steps = np.zeros_like(X, dtype=int)\n", "maxabsr = np.zeros_like(Y)\n", "newtontime = np.zeros_like(Y)\n", "candtime = np.zeros_like(Y)\n", "tic = timeit.default_timer()\n", "rf.flash_many(teqpflsh.PropertyPairs.DT, X, Y, TT, DD, steps, maxabsr, newtontime, candtime)\n", "toc = timeit.default_timer()\n", "print((toc-tic)/N*1e6, np.mean(newtontime), np.mean(candtime), np.mean(steps))" ] }, { "cell_type": "markdown", "id": "dec3beab-33e6-4157-a989-983b260e8e60", "metadata": {}, "source": [ "The timing is carried out at a fairly granular level. The ``candtime`` argument is the time required (in μs) to do preparation of the candidates from the K-D tree values. The ``newtontime`` is the time spent (in μs) preparing the iteration object and actually doing the iteration. In this case the inputs do not require any iteration, but the newton iterator is still constructed." ] }, { "cell_type": "code", "execution_count": 7, "id": "89a67529-6345-4ed4-a4df-39e7373fe7b6", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:50.938661Z", "iopub.status.busy": "2025-01-06T11:32:50.938561Z", "iopub.status.idle": "2025-01-06T11:32:51.263240Z", "shell.execute_reply": "2025-01-06T11:32:51.262956Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "400.04716315487264 1.0000000000000004e-06 5.041 16.0 10.834 1\n" ] } ], "source": [ "# Here we iterate for two variables, it is much slower\n", "for reg in rf.get_regions_rw():\n", " reg.add_pair(proppair=teqpflsh.PropertyPairs.PS, Nsplit=5)\n", " \n", "# Take the points in the K-D tree to do calculations\n", "# They trivially satisfy the stopping conditions!\n", "propset = reg.propset_Trhogrid\n", "X = propset.p\n", "Y = propset.s\n", "o = rf.flash(teqpflsh.PropertyPairs.PS, X[0], Y[0])\n", "print(o.T, o.rho, o.candidate_duration_us, o.total_duration_us, o.newton_duration_us, o.step_count)" ] }, { "cell_type": "code", "execution_count": 8, "id": "3040cf17-0640-4549-b378-784d2dcb22c8", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:51.264534Z", "iopub.status.busy": "2025-01-06T11:32:51.264446Z", "iopub.status.idle": "2025-01-06T11:32:54.482868Z", "shell.execute_reply": "2025-01-06T11:32:54.482589Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.2074755386728637 2.643995782586261 0.46512329597547125 3.0077775161447087 1.4765034340778192e-10 0\n" ] } ], "source": [ "# Now p, s inputs but with a bit of noise in entropy to force the \n", "# Newton iterator to actually do something\n", "\n", "# Input variables\n", "X = propset.p\n", "Y = propset.s + np.random.random(X.shape)\n", "\n", "# Output buffers\n", "TT = np.zeros_like(X)\n", "DD = np.zeros_like(X)\n", "steps = np.zeros_like(X)\n", "maxabsr = np.zeros_like(Y)\n", "newtontime = np.zeros_like(Y)\n", "candtime = np.zeros_like(Y)\n", "tic = timeit.default_timer()\n", "\n", "rf.flash_many(teqpflsh.PropertyPairs.PS, X, Y, TT, DD, steps, maxabsr, newtontime, candtime)\n", "toc = timeit.default_timer()\n", "print((toc-tic)/len(X)*1e6, np.mean(newtontime), np.mean(candtime), np.mean(steps), np.mean(maxabsr), np.sum(TT<0))" ] }, { "cell_type": "markdown", "id": "a6ff2bd2-e001-45ca-b15d-289c0980e9b2", "metadata": {}, "source": [ "Ok, that's good. Iteration was carried out, and the deviations between the specified and iterated entropy and pressure were good. It took approximately 1.1 μs per iteration step of the Newton iterator, which isn't too bad, but we did start not too far from the actual solution." ] } ], "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 }