e5bf977ca5
Rename module directories (module 6 → module_6, module 7 → module_7) to fix Google Colab "Open in Colab" link compatibility. Update Colab badge URL in Peng_Robinson_EOS.ipynb accordingly. Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
319 lines
73 KiB
Text
319 lines
73 KiB
Text
{
|
||
"cells": [
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "73040203-bdb6-4e2e-b631-64da6d08a77d",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Problem 4"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 9,
|
||
"id": "13c903f6-7b59-413c-88b3-7fa98f2c29ea",
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"Computing PR saturation curve for CO2...\n",
|
||
"T= 220.000 K Psat= 5.93635 bar ZL= 0.0117 ZV= 0.9175\n",
|
||
"T= 222.120 K Psat= 6.47773 bar ZL= 0.0128 ZV= 0.9119\n",
|
||
"T= 224.240 K Psat= 7.05598 bar ZL= 0.0139 ZV= 0.9061\n",
|
||
"T= 226.360 K Psat= 7.67271 bar ZL= 0.0151 ZV= 0.9001\n",
|
||
"T= 228.480 K Psat= 8.32954 bar ZL= 0.0164 ZV= 0.8938\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABEEAAAGZCAYAAABmNC8xAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAtbpJREFUeJzs3XlYVGUbx/HvgKwKKCriiriL+76V5m4qaba6b2lpZmqpbWq0aNqilWTZ65q5ZJmllWm554KGG2FaikuKoqKAIshy3j8m0AlQRoFh+X2u61xvc54zZ+6ZFzmHe57nvk2GYRiIiIiIiIiIiORzdrYOQEREREREREQkJygJIiIiIiIiIiIFgpIgIiIiIiIiIlIgKAkiIiIiIiIiIgWCkiAiIiIiIiIiUiAoCSIiIiIiIiIiBYKSICIiIiIiIiJSICgJIiIiIiIiIiIFgpIgIiIiIiIiIlIgKAkiBdLChQsxmUypW6FChShXrhyDBw/mzJkzqcdt3rzZ4jh7e3tKliyJv78/e/fuzfTr7d69m4cffpgKFSrg5OREqVKlaNGiBS+88MJdxf/jjz/y+uuv39Vzs1JoaCivv/46J06cSDM2aNAgKlasmOMxpYiOjubtt9+mcePGuLu74+TkRMWKFRkyZAjBwcE2iysrvf766xY/nxltDzzwQKbOd+LECYvn2dnZUbx4cbp27crOnTuz982IiEi2K4j3I/cq5VqbkYSEBEqVKkXz5s0zPCY5OZkKFSpQt27dTL9uyr1qdrwnkYJOSRAp0BYsWMDOnTvZsGEDw4YNY9myZdx///1cu3bN4ripU6eyc+dONm/ezKRJk9ixYwdt2rThr7/+uuNr/PDDD7Rs2ZLo6GhmzJjB+vXr+fDDD2nVqhUrVqy4q7h//PFHAgIC7uq5WSk0NJSAgIB0L9CTJk3i22+/zfmggGPHjtGgQQPeeecd2rZty7Jly1i/fj0BAQGcP3+eRo0aERUVZZPYstJTTz3Fzp07U7dVq1YB8Nxzz1ns/+STT6w6b8rzt23bxrRp0zhw4ABt27Zl37592fE2REQkBxTU+5Hs5uDgQP/+/dm9ezehoaHpHvPLL79w+vRphg4dmsPRiUh6Ctk6ABFbql27No0bNwagbdu2JCUl8eabb7J69Wr69u2belzVqlVTM/z3338/RYsWZeDAgSxZsuSOF/8ZM2bg6+vLzz//TKFCN//JPfnkk8yYMSMb3tXdi42NxdXVNUvOVbly5Sw5j7WSkpJ4+OGHuXjxIjt37qR27dqpY23atGHgwIH89NNPODg43PNrGYZBXFwcLi4u93yuu1GuXDnKlSuX+jjl5q9ChQq3/UbqTm59fqtWrahSpQrt27fnk08+4fPPP7+nmEVExDYK6v1IThg6dCjvv/8+8+fP57333kszPn/+fBwdHenXr58NohOR/9JMEJFbpPzhd/Lkydsel5I4OX/+/B3PeenSJUqUKGFxw5HCzs7yn+CKFSvo1KkTpUuXxsXFhZo1a/LSSy9ZzEwZNGgQgYGBABZLF06cOJG6nGHhwoVpXstkMllMWU2Z3hkcHMyjjz5KsWLFUhMXe/fu5cknn6RixYq4uLhQsWJFevfubfG5LFy4kMceewwwJ5BS4kh57fSWw8TFxfHyyy/j6+uLo6MjZcuW5dlnn+XKlSsWx1WsWJHu3buzbt06GjZsiIuLCzVq1GD+/Pm3/awBVq9ezaFDh3j55ZctEiC3evDBB1NvrjJatpPe9FeTycSoUaP49NNPqVmzJk5OTvzvf//Dy8uL/v37pznHlStXcHFxYdy4can7oqOjefHFFy0+gzFjxqSZfZSbZPbfhYiI5F4F9X4EzDMx2rdvj7u7O66urrRq1Ypff/01TWw//PAD9evXx8nJCV9f33QTGumpWbMmLVq04IsvviAxMdFi7MqVK3z33Xf06NGD4sWLA/D999/TokULXF1dcXNzo2PHjpladlqxYkUGDRqUZv8DDzxgsfQ1ZTn30qVLmThxIqVLl6ZIkSL4+/tz/vx5YmJiGD58OCVKlKBEiRIMHjyYq1evWpzTMAw++eQT6tevj4uLC8WKFePRRx/l+PHjFsft27eP7t274+XlhZOTE2XKlKFbt278888/mfrsRGxBSRCRW/z9998AlCxZ8rbHhYWFAVCtWrU7nrNFixbs3r2b0aNHs3v3bhISEjI89q+//qJr167MmzePdevWMWbMGL766iv8/f1Tj5k0aRKPPvoogMWSh9KlS98xlvT06tWLKlWqsHLlSj799FPAPKOgevXqzJo1i59//pnp06cTHh5OkyZNuHjxIgDdunVj6tSpAAQGBqbG0a1bt3RfxzAMevbsyXvvvUf//v354YcfGDduHIsWLaJdu3bEx8dbHH/gwAFeeOEFxo4dy3fffUfdunUZOnQoW7duve37Wb9+PQA9e/a8q8/jTlavXs2cOXOYPHkyP//8M+3ataNfv3588803REdHWxy7bNky4uLiGDx4MGD+ZqtNmzYsWrSI0aNH89NPPzFx4kQWLlzIQw89hGEYqc9NuSncvHlztrwPa2T234WIiOReBfV+ZMmSJXTq1Al3d3cWLVrEV199haenJ507d7ZIhPz666/06NEDNzc3li9fzrvvvstXX33FggULMhX/0KFDiYiI4IcffrDYv3TpUuLi4lKXwixdupQePXrg7u7OsmXLmDdvHpcvX+aBBx5g+/btd/XZZeSVV14hIiKChQsX8v7777N582Z69+7NI488goeHB8uWLWPChAl88cUXvPLKKxbPffrppxkzZgwdOnRg9erVfPLJJ/zxxx+0bNky9UvAa9eu0bFjR86fP09gYCAbNmxg1qxZVKhQgZiYmCx9LyJZyhApgBYsWGAAxq5du4yEhAQjJibGWLt2rVGyZEnDzc3NOHfunGEYhrFp0yYDMFasWGEkJCQYsbGxxm+//WZUr17d8PPzMy5fvnzH17p48aJx3333GYABGA4ODkbLli2NadOmGTExMRk+Lzk52UhISDC2bNliAMaBAwdSx5599lkjvX++YWFhBmAsWLAgzRhgTJkyJfXxlClTDMCYPHnyHd9DYmKicfXqVaNw4cLGhx9+mLp/5cqVBmBs2rQpzXMGDhxo+Pj4pD5et26dARgzZsywOG7FihUGYMydOzd1n4+Pj+Hs7GycPHkydd/169cNT09P4+mnn75trF26dDEAIy4u7o7vK704U6R8PrcCDA8PDyMyMtJi/8GDB9O8B8MwjKZNmxqNGjVKfTxt2jTDzs7O2LNnj8VxX3/9tQEYP/74Y+q+gIAAw97e3ti8eXOm3odh3Pz//9133830c9J7/vTp042EhAQjLi7O+P33340mTZoYgPHDDz/c1XlFRMT2CuL9yLVr1wxPT0/D39/fYn9SUpJRr149o2nTpqn7mjVrZpQpU8a4fv166r7o6GjD09Mz3ff4XzExMUaRIkWMhx56yGJ/o0aNjPLlyxtJSUlGUlKSUaZMGaNOnTpGUlKSxXO9vLyMli1bpu5LuVcNCwtL3efj42MMHDgwzWu3adPGaNOmTerjlPvX/77vMWPGGIAxevRoi/09e/Y0PD09Ux/v3LnTAIz333/f4rjTp08bLi4uxoQJEwzDMIy9e/cagLF69erbfzgiuYxmgkiB1rx5cxwcHHBzc6N79+54e3vz008/UapUKYvjnnjiCRwcHFKnUEZHR/PDDz9QtGjRO75G8eLF2bZtG3v27OGdd96hR48eHD16lJdffpk6deqkfpMBcPz4cfr06YO3tzf29vY4ODjQpk0bAA4fPpyl7z3FI488kmbf1atXmThxIlWqVKFQoUIUKlSIIkWKcO3atbuOY+PGjQBppnE+9thjFC5cOM201Pr161OhQoXUx87OzlSrVs3mSzLatWtHsWLFLPbVqVOHRo0aWXxbdPjwYYKCghgyZEjqvrVr11K7dm3q169PYmJi6ta5c+c0sz4mT55MYmJi6v//98owDIvXvHUzbpmBAjBx4kQcHBxwdnamUaNGnDp1is8++4yuXbtmSSwiIpLzCuL9yI4dO4iMjGTgwIEW173k5GS6dOnCnj17uHbtGteuXWPPnj306tULZ2fn1Oe7ublZzH65nSJFivD444/z448/ps6UCAkJ4ffff2fQoEHY2dlx5MgRzp49S//+/S2WIBUpUoRHHnmEXbt2ERsbm6nXy4zu3btbPK5ZsyZAmlm7NWvWJDIyMnVJzNq1azGZTPTr18/ic/P29qZevXqp9ytVqlShWLFiTJw4kU8//TTDwrAiuY0Ko0qBtnjxYmrWrEmhQoUoVapUhlM4p0+fTrt27YiNjWX9+vVMmzaNnj17snv3bpycnDL1Wo0bN06tJZKQkMDEiROZOXMmM2bMYMaMGVy9epX7778fZ2dn3nrrLapVq4arqyunT5+mV69eXL9+Pcve963Se899+vTh119/ZdKkSTRp0gR3d3dMJhNdu3a96zguXbpEoUKF0iypMJlMeHt7c+nSJYv9Ketmb+Xk5HTH109JnISFhVGjRo27ivV2MvoZGTJkCM8++yx//vknNWrUYMGCBTg5OdG7d+/UY86fP8/ff/+dYVHWW29As9qiRYtSl+X814IFCyySU88//zz9+vXDzs6OokWL4uvre9v2gCIikncUpPuRlGREyrKd9ERGRmIymUhOTsbb2zvNeHr7MjJ06FDmz5/PF198wYsvvsj8+fMxmUyp19+Ue5303muZMmVITk7m8uXLWVYU1tPT0+Kxo6PjbffHxcVRpEgRzp8/j2EYab4UTFGpUiUAPDw82LJlC2+//TavvPIKly9fpnTp0gwbNozXXnstS4rQi2QHJUGkQKtZs2bqjcDtVKpUKfW41q1b4+LiwmuvvcbHH3/Miy++aPXrOjg4MGXKFGbOnElISAhgnilx9uxZNm/ebPHt/3+Lht5OyrcX/62v8d8Ew63++8dtVFQUa9euZcqUKbz00kup++Pj44mMjMx0LP9VvHhxEhMTuXDhgkUixDAMzp07R5MmTe763Lfq3Lkzc+fOZfXq1RbxZ8TZ2TnN5wUZJyQySgb07t2bcePGsXDhQt5++22++OILevbsaTFrpESJEri4uGRY4LVEiRJ3jPdu+fv7s2fPnnTHfH19LR6XK1cuU/8uREQkb8vv9yMp19WPP/44w65ppUqVIiEhAZPJxLlz59KMp7cvIy1btqRmzZosWLCA559/niVLltCuXbvU62zKFzzh4eFpnnv27Fns7OzSzDa91e3uWbLyHqJEiRKYTCa2bduW7pd9t+6rU6cOy5cvxzAMDh48yMKFC3njjTdwcXHJ1H2YiC1oOYzIXZgwYQJVqlThnXfeuWPhp/QudHBzOmmZMmWAmxf//15sPvvsszTPTTnmv9+ClCpVCmdnZw4ePGix/7vvvrttjLcymUwYhpEmjv/9738kJSVlKo70tG/fHjAXKLvVN998w7Vr11LH71WPHj2oU6cO06ZNS72h+6+ff/45dbppxYoViYiIsOj0c+PGDX7++WerXrdYsWL07NmTxYsXs3btWs6dO2exFAbM01KPHTtG8eLFU7+Ju3VLr0tNVsnoNRs3bpzurBsREclfCuL9SKtWrShatCihoaEZXgMdHR0pXLgwTZs2ZdWqVcTFxaU+PyYmhjVr1mQ6ZjDPDA0NDeW1117jwoULFvcC1atXp2zZsixdutRiKeq1a9f45ptvUjvGZKRixYppPtOjR49y5MgRq2K8k+7du2MYBmfOnEn3M6tTp06a55hMJurVq8fMmTMpWrQowcHBWRqTSFbSTBCRu+Dg4MDUqVN5/PHH+fDDD3nttdcyPLZz586UK1cOf39/atSoQXJyMvv37+f999+nSJEiPP/884D524NixYrxzDPPMGXKFBwcHPjyyy85cOBAmnOmXHymT5/Ogw8+iL29PXXr1k3tQT9//nwqV65MvXr1CAoKYunSpZl+b+7u7rRu3Zp3332XEiVKULFiRbZs2cK8efPS1EBJaUE7d+5c3NzccHZ2xtfXN90/qjt27Ejnzp2ZOHEi0dHRtGrVioMHDzJlyhQaNGiQbovZu2Fvb8+3335Lp06daNGiBSNGjKBt27YULlyYkydP8vXXX7NmzRouX74MmOu9TJ48mSeffJLx48cTFxfHRx99lOYGKzOGDBnCihUrGDVqFOXKlaNDhw4W42PGjOGbb76hdevWjB07lrp165KcnMypU6dYv349L7zwAs2aNQPgjTfe4I033uDXX3/NsrogIiJScBXU+5GPP/6YgQMHEhkZyaOPPoqXlxcXLlzgwIEDXLhwgTlz5gDw5ptv0qVLFzp27MgLL7xAUlIS06dPp3DhwlbNhB0wYACvvPIK7777LkWLFqVXr16pY3Z2dsyYMYO+ffvSvXt3nn76aeLj43n33Xe5cuUK77zzzm3P3b9/f/r168fIkSN55JFHOHnyJDNmzMjy7m2tWrVi+PDhDB48mL1799K6dWsKFy5MeHg427dvp06dOowYMYK1a9fyySef0LNnTypVqoRhGKxatYorV67QsWPHLI1JJEvZqiKriC2lVNz+b5eO/0qprr1y5cp0x5s1a2YUK1bMuHLlSobnWLFihdGnTx+jatWqRpEiRQwHBwejQoUKRv/+/Y3Q0FCLY3fs2GG0aNHCcHV1NUqWLGk89dRTRnBwcJoK6/Hx8cZTTz1llCxZ0jCZTBbVw6OiooynnnrKKFWqlFG4cGHD39/fOHHiRIbV2C9cuJAm5n/++cd45JFHjGLFihlubm5Gly5djJCQkHSrks+aNcvw9fU17O3tLeJMr+vK9evXjYkTJxo+Pj6Gg4ODUbp0aWPEiBFpuuz4+PgY3bp1SxPXf6uf386VK1eMN99802jYsKHF596vXz/jt99+szj2xx9/NOrXr2+4uLgYlSpVMmbPnp1hd5hnn302w9dMSkoyypcvbwDGq6++mu4xV69eNV577TWjevXqhqOjo+Hh4WHUqVPHGDt2bGpXIsO4+f9Pep13MpJV3WHu9vkiIpJ7FdT7EcMwjC1bthjdunUzPD09DQcHB6Ns2bJGt27d0tzfff/990bdunUNR0dHo0KFCsY777yT7v3AnTz88MMGYIwcOTLd8dWrVxvNmjUznJ2djcKFCxvt27dPc2+SXneY5ORkY8aMGUalSpUMZ2dno3HjxsbGjRsz7A7z3/eX0f1vRv8fzJ8/32jWrJlRuHBhw8XFxahcubIxYMAAY+/evYZhGMaff/5p9O7d26hcubLh4uJieHh4GE2bNjUWLlxo1eclktNMhvGftgAiIiIiIiIiIvmQaoKIiIiIiIiISIGgJIiIiIiIiIiIFAhKgoiIiIiIiIhIgaAkiIiIiIiIiIgUCEqCiIiIiIiIiEiBUMjWAWS35ORkzp49i5ubGyaTydbhiIiI5BuGYRATE0OZMmWws9P3Kpmh+xIREZHskdn7knyfBDl79izly5e3dRgiIiL51unTpylXrpytw8gTdF8iIiKSve50X5LvkyBubm6A+YNwd3e3cTQiIiL5R3R0NOXLl0+91sqd6b5EREQke2T2vsSmSZCKFSty8uTJNPtHjhxJYGAghmEQEBDA3LlzuXz5Ms2aNSMwMJBatWpl+jVSppq6u7vrZkNERCQbaFlH5um+REREJHvd6b7Epgt49+zZQ3h4eOq2YcMGAB577DEAZsyYwQcffMDs2bPZs2cP3t7edOzYkZiYGFuGLSIiIiIiIiJ5kE2TICVLlsTb2zt1W7t2LZUrV6ZNmzYYhsGsWbN49dVX6dWrF7Vr12bRokXExsaydOlSW4YtIiIiIiIiInlQrinlfuPGDZYsWcKQIUMwmUyEhYVx7tw5OnXqlHqMk5MTbdq0YceOHRmeJz4+nujoaItNRERERERERCTXFEZdvXo1V65cYdCgQQCcO3cOgFKlSlkcV6pUqXTriKSYNm0aAQEB2RaniIiIiIhIQZeUlERCQoKtw5ACxMHBAXt7+3s+T65JgsybN48HH3yQMmXKWOz/b1ETwzBuW+jk5ZdfZty4camPUyrEioiIiIiIyL0xDINz585x5coVW4ciBVDRokXx9va+p6LsuSIJcvLkSX755RdWrVqVus/b2xswzwgpXbp06v6IiIg0s0Nu5eTkhJOTU/YFKyIiIiIiUkClJEC8vLxwdXVVhzDJEYZhEBsbS0REBIBFjsBauSIJsmDBAry8vOjWrVvqPl9fX7y9vdmwYQMNGjQAzHVDtmzZwvTp020VqoiIiIiISIGUlJSUmgApXry4rcORAsbFxQUwT4zw8vK666UxNk+CJCcns2DBAgYOHEihQjfDMZlMjBkzhqlTp1K1alWqVq3K1KlTcXV1pU+fPjaMWEREJO9ISjYICoskIiYOLzdnmvp6Ym+nb+1ERMR6KTVAXF1dbRyJFFQpP3sJCQl5Nwnyyy+/cOrUKYYMGZJmbMKECVy/fp2RI0dy+fJlmjVrxvr163Fzc7NBpCIiInnLupBwAtaEEh4Vl7qvtIczU/z96FL77qeRivUCAwMJDAwkKSnJ1qGIiNwzLYERW8mKnz2TYRhGFsSSa0VHR+Ph4UFUVBTu7u62DkdERCRHrAsJZ8SSYP57kU+5dZjTr+E9J0J0jbWePjMRycvi4uIICwvD19cXZ2dnW4cjBdDtfgYze421y+4gRUREJGclJRsErAlNkwABUvcFrAklKTlffw9ScMRFw+bpcP2yrSMRERHJ9ZQEERERyWeCwiItlsD8lwGER8URFBaZc0FJ9gn6DDZPhVl1YdM0uH7F1hGJiAgwb948OnXqZOswGDRoED179kx9/MADDzBmzJgsfY3NmzdjMplu2zp54cKFFC1aNPXx7Nmzeeihh7I0jsxQEkRERCSfiYjJOAFyN8dJLpZ4A3Z/Zv7v+GjY8g58WNc8MyQuyraxiYjkAv7+/nTo0CHdsZ07d2IymQgODr7tOU6cOIHJZGL//v2Zft34+HgmT57MpEmTrAn3tipWrMisWbPu+TyrVq3izTffvPeAbtGyZUvCw8Px8PDI9HOGDRvGnj172L59e5bGcidKgoiIiOQzXm6ZW6ed2eMkFyvkCE/9Cg0HgN2/9e7jom7ODNnyrnm5jIhIATV06FA2btzIyZMn04zNnz+f+vXr07Bhwyx/3W+++YYiRYpw//33Z/m575Wnp2eWNxtxdHTE29vbqsKlTk5O9OnTh48//jhLY7kTJUFERETymaa+nni7Z5zgMGHuEtPU1zPngpLsU8wHHvoYRu2FBv3A9G/LwLgrsOkt88yQre9BfIxNwxQRsYXu3bvj5eXFwoULLfbHxsayYsUKhg4dmi2vu3z58jRLPTZv3kzTpk0pXLgwRYsWpVWrVqnJmWPHjtGjRw9KlSpFkSJFaNKkCb/88kvqcx944AFOnjzJ2LFjMZlM99Ql5b/LYSIiIvD398fFxQVfX1++/PJLi1kn6c2EuXLlCiaTic2bN6e+t/8uh1m4cCEVKlTA1dWVhx9+mEuXLqWJ5aGHHmL16tVcv379rt+PtWzeIldERESylr2diXrlPTj3R9rlLim3TFP8/bC3U4vDfMXTF3oEwn3jzEmPg8vBSDYXTN34JuwMhFajockwcCpi62hFJJ/w/3g7F2Lic/x1S7o5sea5++54XKFChRgwYAALFy5k8uTJqcmDlStXcuPGDfr27Zst8W3bts3i3ImJifTs2ZNhw4axbNkybty4QVBQUGo8V69epWvXrrz11ls4OzuzaNEi/P39OXLkCBUqVGDVqlXUq1eP4cOHM2zYsCyNddCgQZw+fZqNGzfi6OjI6NGjiYiIuKdz7t69myFDhjB16lR69erFunXrmDJlSprjGjduTEJCAkFBQbRp0+aeXjOzlAQRERHJZzb9GcHPf5wHoKirA1diE1LHvD2cmeLvd8/tcSUXK14ZHp4DrV+ELTPg0Ff/JkMi4ZfXYcfH0HI0NB0GjoVtHa2I5HEXYuI5F527a0wNGTKEd999l82bN9O2bVvAvBSmV69eFCtWLMtf78qVK1y5coUyZcqk7ouOjiYqKoru3btTuXJlAGrWrJk6Xq9ePerVq5f6+K233uLbb7/l+++/Z9SoUXh6emJvb4+bmxve3t5ZFuvRo0f56aef2LVrF82aNQPMBV1vje1ufPjhh3Tu3JmXXnoJgGrVqrFjxw7WrVtncVzKrJgTJ04oCSIiIiLWO3vlOmO/2g/AgBY+TPGvRVBYJBExcXi5mZfAaAZIAVG8MvT67GYyJORrczIk9hL8MgV2fAQtn9PMEBG5JyXdnHL969aoUYOWLVsyf/582rZty7Fjx9i2bRvr16/PlthSlnY4O99cmurp6cmgQYPo3LkzHTt2pEOHDjz++OOULm3+UuLatWsEBASwdu1azp49S2JiItevX+fUqVPZEmOKw4cPU6hQIRo3bpy6r0aNGhZdXO72vA8//LDFvhYtWqRJggC4uLgQGxt7T69nDSVBRERE8omEpGSeW7aPK7EJ1CnrwavdamJvZ6JF5eK2Dk1sqURVeORzaD0etkyHkG8A499kyOvw20fQchQ0HQ5OWVsoT0Tyv8wsSckNhg4dyqhRowgMDGTBggX4+PjQvn37bHmt4sWLYzKZuHz5ssX+BQsWMHr0aNatW8eKFSt47bXX2LBhA82bN2f8+PH8/PPPvPfee1SpUgUXFxceffRRbty4kS0xpjAMA+C2NUbs7OwsjgVISEjI6PA0x95JZGQkJUuWzPTx90qFUUVERPKJ934+wu8nL+PmXIjAPg1xKmRv65AkNylZDR6dB8/uhtqPkloh5nok/PoGzKqjbjIikm89/vjj2Nvbs3TpUhYtWsTgwYPvqbjo7Tg6OuLn50doaGiasQYNGvDyyy+zY8cOateuzdKlSwFzDZFBgwbx8MMPU6dOHby9vTlx4kSa8yYlJWVprDVr1iQxMZG9e/em7jty5IhFgdOUBEV4eHjqvju1C/bz82PXrl0W+/77GMwFYePi4mjQoMFdRH93lAQRERHJB349fJ7Pth4H4N1H61GhuKuNI5Jcq2T1f5MhQVD3CTD9ezt4/bK5m8ys2rB5urnVrohIPlGkSBGeeOIJXnnlFc6ePcugQYOsPseRI0fYv3+/xZbRTI3OnTuzffv21MdhYWG8/PLL7Ny5k5MnT7J+/XqOHj2aWnujSpUqrFq1iv3793PgwAH69OlDcnKyxTkrVqzI1q1bOXPmDBcvXgTgzJkz1KhRg6CgIKvfD0D16tXp0qULw4YNY/fu3fz+++889dRTuLi4pB7j4uJC8+bNeeeddwgNDWXr1q289tprtz1vyoyXGTNmcPToUWbPnp3uUpht27ZRqVKl1DopOUFJEBERkTzun8uxjPvqAACDW1WkS+2sK5gm+VjJatBrLjy7B+r1vpkMiYuCzVNhZh3YNA2uX7FpmCIiWWXo0KFcvnyZDh06UKFCBauf/+STT9KgQQOL7ezZs+keO2zYMH788UeioswJZVdXV/78808eeeQRqlWrxvDhwxk1ahRPP/00ADNnzqRYsWK0bNkSf39/OnfuTMOGDS3O+cYbb3DixAkqV66cOjsjISGBI0eO3FNNjQULFlC+fHnatGlDr169GD58OF5eXhbHzJ8/n4SEBBo3bszzzz/PW2+9ddtzNm/enP/97398/PHH1K9fn/Xr16ebOFm2bFmWd7u5E5NhzWKdPCg6OhoPDw+ioqJwd3e3dTgiIiJZ6kZiMo9/tpP9p69Qr5wHK59piWOhnPmOQ9dY6+Xqz+zSMdj2PhxYDsYt062d3KHZ09B8JLh62i4+EbG5uLg4wsLC8PX1tSj6Kel7/PHHU5e/5DUVK1ZkzJgxjBkzJtteIyQkhPbt23P06FE8PDwy9Zzb/Qxm9hqrmSAiIiJ52PR1f7L/9BXcnQsxu0/DHEuASD5UvDL0/ASe2wsN+oHp35oy8dGw9V1zzZANk+HqBdvGKSKSR7z77rsUKaLuWxk5e/YsixcvznQCJKvoTklERCSP+vmPc8zbHgbAe4/Vo7yn6oBIFvCsBD0CYXQwNBwAdv82E7xxFX770JwMWfcyRIff/jwiInnEM888Q5EiRdLdnnnmmbs+r4+PD88991wWRpq/dOrUic6dO+f462o5jIiISB50OjKWrh9tIyYukafu8+W17n45HoOusdbLk5/ZldPw2ywIXgxJtxQAtHeChv2h1RgoWt5W0YlIDsqvy2EiIiKIjk6/M5a7u3ua+hhiO1mxHKZQdgcpIiIiWSs+MYlnlwYTE5dIgwpFmfhgDVuHJPlZ0fLQ7X24/0XY8RHsXQCJ1yEpHvb8D35fBPV7w33jwNPX1tGKiFjNy8tLiY4CRMthRERE8phpP/7JwX+iKOrqwOw+DXGw1+VccoB7aegyDcYchFbPg0Nh8/7kBPMskY8bwbfPwMW/bBuniIjIbeiuSUREJA/56VA4C3ecAOCDx+tRtqiLbQOSgqeIF3R8A8aGQOvx5u4xYO4oc2AZzG4CXw2E8IO2jVNERCQdSoKIiIjkEScvXWPC1+Y/LJ9uU4l2NUrZOCIp0Fw9od1rMOYQtH0VnIv+O2BA6Gr47H5Y+gSc3mPDIEVERCwpCSIiIpIHxCX8WwckPpHGPsV4sVN1W4ckYuZSFNpMMM8M6RAAhUveHDu6DuZ1gEX+cHwL5O96/CIikgcoCSIiIpIHvP3DYULORONZ2JGP+zRQHRDJfZzc4L4x5pkhD74L7uVujoVthcUPwbyOcGSdkiEiImIzuoMSERHJ5dYcOMsXu04C5jogpT1UB0RyMQcXaDYcRu+Dh2aDZ6WbY//sgWVPwKf3Q8gqSE6yXZwiIjmkf//+TJ061dZh5HpNmjRh1apV2f46SoKIiIjkYmEXr/HyqkMAPNu2Mg9UVws/ySMKOULD/jBqLzwyD7z8bo6dPwRfD4bAphD8BSTesF2cIpKv+fv706FDh3THdu7ciclkIjg4OMPn//7775hMJrZv357ueOfOnXnooYcyfP7Bgwf54YcfeO6556wLPAMnTpzAZDKxf//+LDlfbjJp0iReeuklkpOTs/V1lAQRERHJpeISkhj5ZTBX4xNp6uvJ2A7VbB2SiPXs7KHOo/DMb/DkUijT8ObYpb/h+1HwUX3YNQduXLNZmCKSPw0dOpSNGzdy8uTJNGPz58+nfv36NGzYMJ1nmjVq1Ih69eqxYMGCNGOnT5/ml19+YejQoRk+f/bs2Tz22GO4ubnd3RsoAG7cMCfCu3XrRlRUFD///HO2vp6SICIiIrlUwJpQDodHU7ywIx/3bkAh1QGRvMzODmp0g2Ebof+34HPfzbHoM7DuJZhZG7bMgOuXbReniOQr3bt3x8vLi4ULF1rsj42NZcWKFbdNYKQYOnQoX331FdeuWSZqFy5cSMmSJenWrVu6z0tOTmblypVpZop88sknVK1aFWdnZ0qVKsWjjz6aOrZu3Truu+8+ihYtSvHixenevTvHjh1LHff19QWgQYMGmEwmHnjggTvGD/Dyyy/TvHnzNPvr1q3LlClTANizZw8dO3akRIkSeHh40KZNmzSzZEwmE3PmzOHBBx/ExcUFX19fVq5caXHMoUOHaNeuHS4uLhQvXpzhw4dz9erV1PFBgwbRs2dPpk2bRpkyZahWzfwlj729PV27dmXZsmWZek93q1C2nl1ERETuynf7z7As6BQmE8x6sj6l3J1tHZLcg8DAQAIDA0lKUg0MTCao3M68nQ6CbR/A0Z/MY9cjYdPb8NuH0HgItHgW3LxtG6+I3N5nbeBqRM6/bhEveHrLHQ8rVKgQAwYMYOHChUyePBmTyQTAypUruXHjBn379r3jOfr27cv48eNZuXIlgwYNAsAwDBYuXMjAgQMpVCj9P6sPHjzIlStXaNy4ceq+vXv3Mnr0aL744gtatmxJZGQk27ZtSx2/du0a48aNo06dOly7do3Jkyfz8MMPs3//fuzs7AgKCqJp06b88ssv1KpVC0dHxzvGn/Ie3nnnHY4dO0blypUB+OOPPzh06BBff/01ADExMQwcOJCPPvoIgPfff5+uXbvy119/WcxkmTRpEu+88w4ffvghX3zxBb1796Z27drUrFmT2NhYunTpQvPmzdmzZw8RERE89dRTjBo1yiIR9euvv+Lu7s6GDRswbimW3bRpU2bMmJGp93S3TIaRv8tzR0dH4+HhQVRUFO7u7rYOR0REJF1JyQZBYZFExMSRkGQwafUhrickM7pdFcbl0na4usZaT59ZBs7/AdtnQsg3YNyyFtzeCer3gVbPg6ev7eITEQDi4uIICwvD19cXZ+d/k/Pv14SYszkfjFsZeOFwpg79888/qVmzJhs3bqRt27YAtGnThrJly7J06dJMnePJJ58kPDycLVvMiZdNmzbRrl07/vzzT6pXT/86vXr1ah599FESEhJSky+rVq1i8ODB/PPPP5laInPhwgW8vLw4dOgQtWvX5sSJE/j6+rJv3z7q16+fqdhT1KtXj0cffZRJkyYB8Morr/DLL78QFBSU7vFJSUkUK1aMpUuX0r17d8A8E+SZZ55hzpw5qcc1b96chg0b8sknn/D5558zceJETp8+TeHChQH48ccf8ff35+zZs5QqVYpBgwaxbt06Tp06lSaJ8/333/Pwww+TkJCAnV3aGbDp/gz+K7PXWM0EERERsbF1IeEErAklPCrOYn+1UkV4XnVApCAoVQse+R+0fRV2fAT7voSkePP2+wIIXgS1esF9Y8G7tq2jFZFbFbFRwW4rXrdGjRq0bNmS+fPn07ZtW44dO8a2bdtYv359ps8xdOhQOnXqxN9//02VKlWYP38+rVq1yjABAnD9+nWcnJxSEyAAHTt2xMfHh0qVKtGlSxe6dOnCww8/jKurKwDHjh1j0qRJ7Nq1i4sXL6YWCT116hS1a9/b77++ffsyf/58Jk2ahGEYLFu2jDFjxqSOR0REMHnyZDZu3Mj58+dJSkoiNjaWU6dOWZynRYsWaR6nFGo9fPgw9erVS02AALRq1Yrk5GSOHDlCqVKlAKhTp066s1hcXFxITk4mPj4eF5fs6YanJIiIiIgNrQsJZ8SSYNKblnn0/FU2hJ6jS+3SOR6XiE14+kL3mdBmIuz6BPbMhxsx5tkhIV+btyodzckQn5bmpTUiYluZWJKSGwwdOpRRo0YRGBjIggUL8PHxoX379pl+focOHfDx8WHhwoVMmDCBVatWMXv27Ns+p0SJEsTGxnLjxo3UP/jd3NwIDg5m8+bNrF+/nsmTJ/P666+zZ88eihYtir+/P+XLl+fzzz+nTJkyJCcnU7t27dTiofeiT58+vPTSSwQHB3P9+nVOnz7Nk08+mTo+aNAgLly4wKxZs/Dx8cHJyYkWLVpk6rVTEj2GYVgkfdI7BrBIktwqMjISV1fXbEuAgAqjioiI2ExSskHAmtB0EyAAJszFUZOS8/XKVZG03Lyh4xsw9hC0ew1ci98c+3sDLOwK8zrCnz9ANrdSFJH84fHHH8fe3p6lS5eyaNEiBg8enOEf6+kxmUwMHjyYRYsWsXTpUuzs7Hj88cdv+5yU5SqhoaEW+wsVKkSHDh2YMWMGBw8e5MSJE2zcuJFLly5x+PBhXnvtNdq3b0/NmjW5fNmyUHRKMuVuakyVK1eO1q1b8+WXX/Lll1/SoUOH1JkZANu2bWP06NF07dqVWrVq4eTkxMWLF9OcZ9euXWke16hRAwA/Pz/2799vUUT2t99+w87OLrUA6u2EhITctltPVlASRERExEaCwiLTLIG5lQGER8URFBaZc0GJ5CYuxaD1eBgTAg++Cx4Vbo79sweW94FPmpuXzyTe+7ekIpJ/FSlShCeeeIJXXnmFs2fPphY4tcbgwYM5e/Ysr7zyCk8++WSGsxlSlCxZkoYNG7J9+/bUfWvXruWjjz5i//79nDx5ksWLF5OcnEz16tUpVqwYxYsXZ+7cufz9999s3LiRcePGWZzTy8sLFxcX1q1bx/nz54mKigLg22+/TU1E3E7fvn1Zvnw5K1eupF+/fhZjVapU4YsvvuDw4cPs3r2bvn37pjsjY+XKlcyfP5+jR48yZcoUgoKCGDVqVOr5nZ2dGThwICEhIWzatInnnnuO/v37WyRcMrJt2zY6dep0x+PuhZIgIiIiNhIRk3EC5G6OE8m3HF2h2XAYHQy9Pgcvv5tjF4/AdyPho/qwMxDir2Z4GhEp2IYOHcrly5fp0KEDFSpUuPMT/qNChQp06NCBy5cvM2TIkEw9Z/jw4Xz55Zepj4sWLcqqVato164dNWvW5NNPP2XZsmXUqlULOzs7li9fzu+//07t2rUZO3Ys7777rsX5ChUqxEcffcRnn31GmTJl6NGjBwBRUVEcOXLkjvE89thjXLp0idjYWHr27GkxNn/+fC5fvkyDBg3o378/o0ePxssrbe2VgIAAli9fTt26dVm0aBFffvklfn7m38uurq78/PPPREZG0qRJEx599FHat29/x6VDAGfOnGHHjh0MHjz4jsfeC3WHERERsZGdxy7R+/Nddzxu2bDmtKhc/I7H5TRdY62nzyyLGAb8tcHcUebUDssx56LQdDg0exoKl7BJeCL51e06c0j64uLiqF69OsuXL09TUDQvMplMfPvtt2kSKFlh/PjxREVFMXfu3AyPyYruMJoJIiIiYiNNfT0p4pRxjXITUNrDmaa+njkXlEheYDJBtU4w5CcYsh6qd705FncFts6AmbXghxcg8rjNwhQRcXZ2ZvHixenW1hBLXl5evPnmm9n+OkqCiIiI2Mivh89zNT4x3bGUUm1T/P2wt1MHDJEMVWgGvZfByF1Qrw/Y/ZtYTIyDPf+DjxvBykFwJtimYYpI7vXMM89QpEiRdLdnnnnmns/fpk0b/P39syDS/G38+PGZqhtyr9QiV0RExAb+jrjKuK8OANC2ekn+PBdjUSTV28OZKf5+ao8rklleNeHhOdD2FXN73d8XQcI1c3vdP741b76todXzULm92uuKSKo33niDF198Md0xLV20lB+qadg8CXLmzBkmTpzITz/9xPXr16lWrRrz5s2jUaNGgPlDDggIYO7cuVy+fJlmzZoRGBhIrVq1bBy5iIjI3YmJS+DpL/ZyNT6Rpr6ezB3QGDuTiaCwSCJi4vByMy+B0QwQkbtQtDx0mQZtJsCeebD7U7h2wTwWttW8laoNLUdD7V5g72DbeEXE5ry8vNItACr5k02Xw1y+fJlWrVrh4ODATz/9RGhoKO+//z5FixZNPWbGjBl88MEHzJ49mz179uDt7U3Hjh2JiYmxXeAiIiJ3KTnZ4IWvDnDswjVKezgT2KchDvZ22NuZaFG5OD3ql6VF5eJKgIjcK5di0PpFc3vd7rPAs/LNsfMh8O1w+LA+7PxEHWVErJQfZgNI3pQVP3s2TYJMnz6d8uXLs2DBApo2bUrFihVp3749lSubL1KGYTBr1ixeffVVevXqRe3atVm0aBGxsbEsXbrUlqGLiIjclcBNf7M+9DyO9nbM6deIkm5Otg5JJH9zcIbGg2HUHnj8Cyjb+OZY9D/w88sw0w9+CYCYc7aLUyQPcHAwz5yKjY21cSRSUKX87KX8LN4Nmy6H+f777+ncuTOPPfYYW7ZsoWzZsowcOZJhw4YBEBYWxrlz5+jUqVPqc5ycnGjTpg07duzg6aefTnPO+Ph44uPjUx9HR0dn/xsRERHJhI1/nueDX44C8FbP2tQvX9S2AYkUJHb24PcQ1PSHUzvhtw/h6DrzWFwUbP8Ads6GOo9Dy+fAq4Zt4xXJhezt7SlatCgREREAuLq6YlJ9HckBhmEQGxtLREQERYsWxd7e/q7PZdMkyPHjx5kzZw7jxo3jlVdeISgoiNGjR+Pk5MSAAQM4d86cjf9vhdhSpUpx8uTJdM85bdo0AgICsj12ERGR20lKNixqfJR0c+L55fsxDOjXvAKPNylv6xBFCiaTCXxamreIw7DjYzj4FSQnQNIN2L/EvFXtZK4bUvE+FVEVuYW3tzdAaiJEJCcVLVo09WfwbpkMGy7ocnR0pHHjxuzYsSN13+jRo9mzZw87d+5kx44dtGrVirNnz1K69M3q+MOGDeP06dOsW7cuzTnTmwlSvnx5oqKiVNlXRERyxLqQcALWhFp0eylkZyIx2aCRTzGWDWuOY6G836U+OjoaDw8PXWOtoM8sl4oONxdQ3bsA4qMsx0rXN88M8esJ9jbvKSCSayQlJZGQkGDrMKQAcXBwuO0MkMxeY236m7x06dL4+flZ7KtZsybffPMNcDPLeO7cOYskSERERIb9g52cnHBy0vpqERGxjXUh4YxYEsx/v2FITDbvebxRuXyRABHJV9xLQ8cAcyHV4C/MLXajTpvHwvfDN0PNNUNajIQG/cGpiE3DFckN7O3t72lJgoit2PQurFWrVhw5csRi39GjR/Hx8QHA19cXb29vNmzYkDp+48YNtmzZQsuWLXM0VhERkTtJSjYIWBOaJgFyq1m//kVSsqrqi+RKTm7mRMfoffDIPPCue3Ms6hSse+nfIqqvm2ePiIhInmPTJMjYsWPZtWsXU6dO5e+//2bp0qXMnTuXZ599FgCTycSYMWOYOnUq3377LSEhIQwaNAhXV1f69Oljy9BFRETSCAqLtFgCk57wqDiCwiJzKCIRuSv2DlDnUXh6Kwz4Hqp0vDkWFwXbZ8KsOvDtM3DukO3iFBERq9l0OUyTJk349ttvefnll3njjTfw9fVl1qxZ9O3bN/WYCRMmcP36dUaOHMnly5dp1qwZ69evx83NzYaRi4iIpBURc/sEiLXHiYiNmUxQqY15Ox9q7h6TUkQ1OQEOLDNvlR6AFs9BlfYqoioiksvZtDBqTlABMhERySk7j12i9+e77njcsmHNaVG5eA5ElL10jbWePrN8IOYcBM2FPfMg7orlWMka0OJZc5tdB2ebhCciUlBl9hqrymwiIiJZpKmvJ97uGf/hYwJKezjT1Ncz54ISkazl5g3tJ8O4UOj6HhTzvTl24U/4/jmYVRu2zIBrl2wXp4iIpEtJEBERkSxib2fCr0z63zykTJCf4u+HvZ2my4vkeY6FoekweO53eGIJlG9+c+zaBdj0NsysBWvHwsW/bBeniIhYUBJEREQki6zYc4qNf0YAUNTVwWLM28OZOf0a0qV26fSeKiJ5lZ091PSHoT/DU7+CX08w/XuLnXgd9s6H2Y3hy8fh+BbI3yvRRURyPZsWRhUREckvgsIieW11CABjO1RjVLsqBIVFEhETh5ebeQmMZoCI5HPlGsPji+DyCdj1Kez7Am5cNY/99bN5K1XbXDek9iNQyMmm4YqIFEQqjCoiInKPTkfG0iPwNyKv3aBbndLM7tMAUwHoEKFrrPX0mRUw169A8GLY/RlE/2M5VqQUNBkGjYdA4bxfKFlExNZUGFVERCQHXItPZNjivUReu0Htsu6891i9ApEAEZFMcCkKrUbD8wfg0flQttHNsavnYdNbMNMP1oyBC0dtFaWISIGiJIiIiMhdSk42GLNiP3+ei6FEESfm9m+Mi6O9rcMSkdzGvpB5+ctTv8KQ9eDX45a6IXHw+wIIbAJLHoVjG1U3REQkG6kmiIiISCYlJRsWdT62/XWBDaHncSxkx9wBjShT1MXWIYpIbmYyQYVm5u3yCfMymeAv4EaMefzvDeatZE1oPgLqPg4O+r0iIpKVVBNEREQkE9aFhBOwJpTwqLg0Yx88Xo9eDcvZICrb0jXWevrMJI24KHMiZPdnEHXKcsy1uLlmSJOnwM3bNvGJiOQRqgkiIiKSRdaFhDNiSXC6CRAAVy2BEZG75ewBLUfB6H3w2CIo3/zmWOwl2PouzKwNq4bD2f02C1NEJL9QEkREROQ2kpINAtaEktG0SRMQsCaUpOR8PbFSRLKbfSGo1ROG/gzDNkKdx8Du35XryQlwcAXMbQPzH4TQ7yE5yabhiojkVUqCiIiI3EZQWGSGM0AADCA8Ko6gsMicC0pE8reyjeCR/8GYQ3DfOHApdnPs1A74qj981AB2zDYvpxERkUxTEkREROQ2ImIyToDczXEiIpnmXgY6TIGxodB9FpSofnPsyklY/yq8XxN+HA8X/7ZZmCIieYmSICIiIrfh5eacpceJiFjN0RUaD4aRu6DfN1C5/c2xhGsQNBdmN4IvH4O/f1WLXRGR21CLXBERkdto6uuJu7MD0XEJ6Y6bAG8PZ5r6euZsYCJS8NjZQZUO5u3CEdj9KRxYDgmx5vG/1pu3kjWg2dNQ90lzAkVERFJpJoiIiMht/Pb3Ra7GZ5wAAZji74e9nSndY0REskXJ6tB9JowLhY5vgEf5m2MX/oS1Y+GDmrBhCkT9Y7s4RURyGSVBREREMnDkXAzPfhlMsgHNfD3x9rBc8uLt4cycfg3pUru0jSIUkQLPpRi0eh5G7ze32K3Q4uZY3BX4bRbMqgtfDYSTO7VURkQKPC2HERERSUdETBxDFu4hJj6Rpr6eLB7alEJ2dgSFRRIRE4eXm3kJjGaAiEiukNJit1ZPOLsPdn0KId+Y2+saSRC62rx514XmI6BWL3BQLSMRKXhMhpG/08HR0dF4eHgQFRWFu7u7rcMREZE84PqNJJ6cu5MD/0ThW6Iwq0a0pFhhR1uHlevoGms9fWaSo2LOw975sHceXLtgOeZawlxstfFQcNdsNhHJ+zJ7jdVyGBERkVskJxuMWbGPA/9EUczVgQWDmigBIiJ5k1spaPsyjP0DHv4MSte/ORZ7Eba+C7Nqw9dD4HSQlsqISIFgVRIkKSmJLVu2cPny5eyKR0REJEclJRvsPHaJ7/afYeexS0z96TA//3EeR3s75g5oTMUShW0douQya9eupXr16lStWpX//e9/tg5H5M4KOUG9J2H4Zhi6AWo/Anb/ropPTjQvm5nXET5vBwdWQGK8TcMVEclOVi+HcXZ25vDhw/j6+mZXTFlK005FRCQj60LCCVgTSnhUXJqxD5+sT4/6ZW0QVd5REK+xiYmJ+Pn5sWnTJtzd3WnYsCG7d+/G0zNzLZIL4mcmuVT02X+Xyiwwzwq5VWGvf5fKDAE3b9vEJyJipWxbDlOnTh2OHz9+T8GJiIjY2rqQcEYsCU43AQLgVEgrRiWtoKAgatWqRdmyZXFzc6Nr1678/PPPtg5LxHruZaDda+alMj0+MRdMTXEtArZMh5m14OuhWiojIvmK1Xd4b7/9Ni+++CJr164lPDyc6Ohoi01ERCS3S0o2CFgTSka39CYgYE0oScm66c9vtm7dir+/P2XKlMFkMrF69eo0x3zyySf4+vri7OxMo0aN2LZtW+rY2bNnKVv25gyhcuXKcebMmZwIXSR7ODhDg77w9FYYvA78eoLJ3jyWnAghX5uXysx9APYvhYT0E8ciInmF1UmQLl26cODAAR566CHKlStHsWLFKFasGEWLFqVYsWLZEaOIiEiWCgqLzHAGCIABhEfFERQWmXNBSY64du0a9erVY/bs2emOr1ixgjFjxvDqq6+yb98+7r//fh588EFOnToFQHqriE2mjNskx8fH6wsjyRtMJvBpAY8vgjGH4P4XwbX4zfHw/bB6hHl2yK9vQpSSfyKSNxWy9gmbNm3KjjhERERyTERM5r7JzOxxknc8+OCDPPjggxmOf/DBBwwdOpSnnnoKgFmzZvHzzz8zZ84cpk2bRtmyZS1mfvzzzz80a9Ysw/NNmzaNgICArHsDIjnBoyy0nwStx8Mfq2D3Z+YkCJjrh2x7D7bPhJr+0OxpqNDCnEQREckDrE6CtGnTJjviEBERyTFebs5ZepzkDzdu3OD333/npZdestjfqVMnduzYAUDTpk0JCQnhzJkzuLu78+OPPzJ58uQMz/nyyy8zbty41MfR0dGUL18+e96ASFZzcIb6faBeb/hnD+z+FEK/My+TMZIgdLV5K1UHmg6DOo+Bo6utoxYRuS2rkyApYmNjOXXqFDdu3LDYX7du3QyeISIikjs09imGs4MdcQnJ6Y6bAG8PZ5r6Zq7jh+QPFy9eJCkpiVKlSlnsL1WqFOfOnQOgUKFCvP/++7Rt25bk5GQmTJhA8eLF0zsdAE5OTjg5OWVr3CLZzmSC8k3NW3Q4/L7A3FXmWoR5/PwhWDMaNkyGhv2h8VDwzBudJEWk4LE6CXLhwgUGDx7MTz/9lO54UlLSPQclIiKSnd7bcOS2CRCAKf5+2NtpendB9N8aH4ZhWOx76KGHeOihh3I6LJHcwb00tH0F7n8B/lgNQXPhzF7zWNwV2PEx7JgN1bqYZ4dUagt26rYlIrmH1b+RxowZw+XLl9m1axcuLi6sW7eORYsWUbVqVb7//vvsiFFERCTLLPgtjM+2mFu9D2zhQ2kPyyUv3h7OzOnXkC61S9siPLGhEiVKYG9vnzrrI0VERESa2SEiBV4hJ6j3BAz7FYZthHp9wN7x30EDjv4ES3pBYBNzTZE4FQUWkdzB6pkgGzdu5LvvvqNJkybY2dnh4+NDx44dcXd3Z9q0aXTr1i074hQREblnPx4K5421oQCM71ydZ9tWYbJ/LYLCIomIicPLzbwERjNACiZHR0caNWrEhg0bePjhh1P3b9iwgR49etgwMpFcrmwjeLgRdHoTghfBnnkQ/W8B4Ut/w08T4Nc3oN6T0HQ4lKxu23hFpECzOgly7do1vLy8APD09OTChQtUq1aNOnXqEBwcnOUBioiIZIXdxy8xZsV+DAP6N/dh5AOVAbC3M9GicsY1HSR/uXr1Kn///Xfq47CwMPbv34+npycVKlRg3Lhx9O/fn8aNG9OiRQvmzp3LqVOneOaZZ2wYtUgeUbiEeZlMy+fhyI/mpTIntpnHblyFPf8zb76tzcmQag+C/V2XKBQRuStW/9apXr06R44coWLFitSvX5/PPvuMihUr8umnn1K6tKYOi4hI7nP0fAzDFu/lRmIynfxK8fpDtdLUfZCCYe/evbRt2zb1cUrnloEDB7Jw4UKeeOIJLl26xBtvvEF4eDi1a9fmxx9/xMfH555eNzAwkMDAQNVOk4LBvhD4PWTezofCns/hwHJIiDWPh201b+7loPFgaDgQipS0bcwiUmCYDMMwrHnCl19+SUJCAoMGDWLfvn107tyZS5cu4ejomHrzkJtER0fj4eFBVFQU7u7utg5HRESyWVKyYbG8pbynC499upPwqDga+RTjy6ea4exgb+sw8wVdY62nz0wKrOtX4MAyCPocIo9Zjtk7Qq2HzbNDyjYyd6MREbFSZq+xVidB/is2NpY///yTChUqUKJEiXs5VbbQzYaISMGxLiScgDWhhEfFpe4rZGciMdmgcsnCfP1MS4oVdrzNGcQa93KNPXjwYKaPrVu3rrWh5Vq6L5ECLzkZjm+EoP/B0XXAf/4UKdMAmgyD2r3AwcUmIYpI3pQjSZCUp+bmKcW62RARKRjWhYQzYknwf2+nU73dszZ9m9/bkgaxdC/XWDs7O0wmExndhqSMmUymfLWERPclIre4fMJcRHXfF3D9suWYiyc07A+Nh0Ix/e4WkTvL1iTIvHnzmDlzJn/99RcAVatWZcyYMTz11FN3H3E20c2GiEj+l5RscN/0jRYzQP6rtIcz2ye2U+eXLHQv19iTJ09m+th7rceRm+i+RCQdCdch5BtzIdXwA/8ZNEG1zubZIZXbgZ2dTUIUkdwvs9dYqwujTpo0iZkzZ/Lcc8/RokULAHbu3MnYsWM5ceIEb731VqbP9frrrxMQEGCxr1SpUpw7dw4wzzQJCAhg7ty5XL58mWbNmhEYGEitWrWsDVtERPKxoLDI2yZAAMKj4ggKi1QnmFwiPyU2ROQeObhAg35Qvy/8s9ecDPnjW0hOAAzzspmj68CzknlmSIO+4FLM1lGLSB5ldRJkzpw5fP755/Tu3Tt130MPPUTdunV57rnnrEqCANSqVYtffvkl9bG9/c1idTNmzOCDDz5g4cKFVKtWjbfeeouOHTty5MgR3NzcrA1dRETyqYiY2ydArD1Oct6xY8eYNWsWhw8fxmQyUbNmTZ5//nkqV65s69BEJKeYTFC+iXnr/DYEL4K9CyD6jHk88jisfxU2vgV1HoWmw6B0PdvGLCJ5jtXzyZKSkmjcuHGa/Y0aNSIxMdHqAAoVKoS3t3fqVrKkuT2WYRjMmjWLV199lV69elG7dm0WLVpEbGwsS5cuzfB88fHxREdHW2wiIpK/ebk5Z+lxkrN+/vln/Pz8CAoKom7dutSuXZvdu3dTq1YtNmzYYOvwRMQWinhB6/Hw/EF4YglUeuDmWOJ1cx2Rz1rD/zrCwa8gMd5moYpI3mJ1EqRfv37MmTMnzf65c+fSt29fqwP466+/KFOmDL6+vjz55JMcP34cgLCwMM6dO0enTp1Sj3VycqJNmzbs2LEjw/NNmzYNDw+P1K18+fJWxyQiInlLU19Piro4ZDhuwlwTpKmvZ84FJZn20ksvMXbsWHbv3s0HH3zAzJkz2b17N2PGjGHixIm2Di9LBAYG4ufnR5MmTWwdikjeYl8IavrDgO/g2T3Q7BlwumWt/z9BsGoYzKwFv74BV07bLlYRyRMyVRh13Lhxqf+dmJjIwoULqVChAs2bNwdg165dnD59mgEDBvDxxx9n+sV/+uknYmNjqVatGufPn+ett97izz//5I8//uDIkSO0atWKM2fOUKZMmdTnDB8+nJMnT/Lzzz+ne874+Hji429mgqOjoylfvrwKkImI5GOb/oxg6KI9JKdzRUspgzqnX0O61C6do3Hld1lV5NPZ2ZlDhw5RtWpVi/1Hjx6lbt26xMXln2VMKowqkgXir8Khr8xtdiP+sBwz2UG1B6HJUKjUVoVURQqQLC2Mum/fPovHjRo1AszrdwFKlixJyZIl+eOPP9I893YefPDB1P+uU6cOLVq0oHLlyixatCg1wfLf9rsp7fIy4uTkhJOTk1VxiIhI3rXnRCTPLPmdZAOaVCzG6cjrnIu++Uezt4czU/z9lADJxUqWLMn+/fvTJEH279+Pl5eXjaISkVzLqQg0HgKNBsOpnRD0ORz+HpITwUiGIz+YN8/K5mRI/T4qpCoiqTKVBNm0aVN2xwFA4cKFqVOnDn/99Rc9e/YE4Ny5c5QuffPGNSIiglKlSuVIPCIikruFno1myMI9xCcm066GF5/1b4SdyURQWCQRMXF4uZmXwKgtbu42bNgwhg8fzvHjx2nZsiUmk4nt27czffp0XnjhBVuHJyK5lckEPi3NW8z5m4VUY86axyOPwc+vwK9vQp1HzG12y9S3acgiYntWd4fJTvHx8Rw+fJj7778fX19fvL292bBhAw0aNADgxo0bbNmyhenTp9s4UhERsbUTF68xYH4QMXGJNKlYjMA+DXGwN097VhvcvGXSpEm4ubnx/vvv8/LLLwNQpkwZXn/9dUaPHm3j6EQkT3ArBW0mwH3j4MiPsOd/ELbFPJZ4HfYtMW9lG0OTp6DWw+CgYtkiBVGmaoJklxdffBF/f38qVKhAREQEb731Flu2bOHQoUP4+Pgwffp0pk2bxoIFC6hatSpTp05l8+bNVrXI1dpbEZG8LynZsJjdUcHTlSfm7uSfy9epWdqd5cOb43GbwqiSPbLjGhsTEwOQ6et8XqP7EpEcdOEo7J0H+5dC/H86Rrp4QsP+5iU1nr62iU9EslSW1gTJLv/88w+9e/fm4sWLlCxZkubNm7Nr1y58fHwAmDBhAtevX2fkyJFcvnyZZs2asX79+nx7YyQiImmtCwknYE0o4VE363wUsjORmGxQsbgri4c0VQIkH9E1XkSyTMlq8OB0aDcJDq00zw45H2Ieux4Jv30Iv30EVTuaZ4dU6QB29raNWUSynU1nguQEfeMiIpJ3rQsJZ8SSYDK6UL3Vszb9mvvkaExyU1ZdYy9dusTkyZPZtGkTERERJCcnW4xHRkbea6i5hu5LRGzIMOD0bnMy5I/VkJxgOe5RARoPhoYDoHAJm4QoIncvT8wEERERyUhSskHAmtAMEyAAgZv+pnfTCip8msf169ePY8eOMXToUEqVKnXbLnAiInfNZIIKzc1b52mwb7G5kGrUafN41Cn4NQA2TwO/nubZIeWbmp8nIvlGliZBtm7dSr169fDw8MjK04qISAEUFBZpsQQmPeFRcQSFRaoQah63fft2tm/fTr169WwdSrYJDAwkMDCQpKQkW4ciIgBFSsL9L0CrMfDXBvPskL9/AQxIugGHvjJvpeqY2+zWeczcmldE8jy7rDzZAw88QKVKlXj//fez8rQiIlIARcTcPgFi7XGSe9WoUYPr16/bOoxs9eyzzxIaGsqePXtsHYqI3MrOHqp3gX5fw+hgaDkaXIrdHD9/CNaOgQ9qwo/jIeJPm4UqIlkjS5MgYWFhfPPNN1y8eDErTysiIgWQl1vmWhdm9jjJvT755BNeffVVtmzZwqVLl4iOjrbYRERyhGcl6PQmjDsMPT81t9NNER8NQXPhk2awsDuErILEG7aLVUTuWpYuh/Hx8cHHx4cHHnggK08rIiIFUJOKxXB1tCf2RvrLB0yAt4czTX09czYwyXJFixYlKiqKdu3aWew3DAOTyaQlJCKSsxxcoH5v83Z2v7nN7sGVkPjvjLUT28xbkVLQcCA0GgQeZW0ZsYhY4a6SIFeuXOHrr7/m2LFjjB8/Hk9PT4KDgylVqhRly+oXgIiI3BvDMHhv/dHbJkAApvj7qShqPtC3b18cHR1ZunSpCqOKSO5Spj489DF0fBMOLIM98+DSX+axq+dh6wzY9h5U7wqNh0CltmCXpZPtRSSLWZ0EOXjwIB06dMDDw4MTJ04wbNgwPD09+fbbbzl58iSLFy/OjjhFRKQA+ejXv/l0yzEAejctz+YjFyyKpHp7ODPF348utUvbKkTJQiEhIezbt4/q1avbOhQRkfS5FIXmI6DZMxC21VxI9c8fwEgCIxn+XGvePCuZkyH1+4KrZiqK5EZWJ0HGjRvHoEGDmDFjBm5ubqn7H3zwQfr06ZOlwYmISMHz2ZZjzPzlKACvdavJU/dXIinZICgskoiYOLzczEtgNAMk/2jcuDGnT59WEkREcj+TCSq1MW/RZyF4Mfy+EGLCzeORx2H9a7DxLajVy9xmt2xDtdkVyUWsToLs2bOHzz77LM3+smXLcu7cuSwJSkRECqZFO04w7Sdz5f0XO1XjqfsrAWBvZ1Ib3Hzsueee4/nnn2f8+PHUqVMHBwcHi/G6devaKDIRkdtwLwMPvGRutXvkJ/PskLAt5rHEODiw1LyVrgeNh0KdR8GxsG1jFhHrkyDOzs7pVmo/cuQIJUuWzJKgREQk//vv7I4Tl64y5fs/ABjVtgqj2lW1cYSSU5544gkAhgwZkrrPZDKpMKqI5A32DuD3kHm7+BfsnQ/7v4S4KPN4+AFYMxrWTzIXW208FEpWs23MIgWY1UmQHj168MYbb/DVV18B5puUU6dO8dJLL/HII49keYAiIpL/rAsJJ2BNqEWdjxRD7/PlhU66OSxIwsLCbB2CiEjWKFEVukyDdpMg5Bvz7JDw/eax+CjY/al5q3g/NBkKNbqbkygikmOsToK89957dO3aFS8vL65fv06bNm04d+4cLVq04O23386OGEVEJB9ZFxLOiCXBGBmMN/Yppu4gBcQrr7xCz549adq0qa1DERHJWo6u0LC/eTvzO+yZDyFfm5fJwC1tdr2h4QC12RXJQSbDMDK6D72tjRs3EhwcTHJyMg0bNqRDhw5ZHVuWiI6OxsPDg6ioKNzd3W0djohIgZaUbHDf9I3pzgABc+tbbw9ntk9sp8KnecC9XmMHDx7MDz/8gL29Pf7+/vTo0YMOHTrg5OSUDdHaVmBgIIGBgSQlJXH06FHdl4gURLGRN9vsRh6zHDPZQ/UH1WZX5B5k9r7EqiRIYmIizs7O7N+/n9q1a2dJoNlNSRARkdxj57FL9P581x2PWzasuQqh5gFZcY01DIPt27ezZs0avv/+e86cOUPHjh156KGH6N69OyVKlMjiqG1L9yUiQnIyhG02J0OO/GRus3srtdkVuSuZvcZalWIsVKgQPj4+KlAmIiJ3JSIm/Rkgd3uc5H0mk4n777+fGTNm8OeffxIUFETz5s35/PPPKVu2LK1bt+a9997jzJkztg5VRCRr2NlB5Xbw5Jcw5hC0mWheFpMipc3uBzXh2xHwz+9wd5P3RSQdVs+zeu2113j55ZeJjIzMjnhERCQf83JzztLjJP+pWbMmEyZM4LfffuOff/5h0KBBbNu2jWXLltk6NBGRrOdRFtq+AmND4PHF4Nv65lhKm93/tYO5beD3RXDjmu1iFcknrC6M+tFHH/H3339TpkwZfHx8KFzYstd1cHBwlgUnIiL5y/UbibcdT6kJ0tRX038FLl68yFtvvcXx48dtHYqISPaydwC/HubtwtF/2+wuNXeUAbXZFclCVidBevbsmQ1hiIhIfrf5SATPfHkzUW4Ciw4xKWVQp/j7qSiqAHDjxg1Onjxp6zBERHJWyWrw4DvQPqXN7rzbtNl9Cmp0U5tdEStYnQSZMmVKdsQhIiL52OYjEQz/4nduJCbTya8U/vXKMPXHwxZdYrw9nJni70eX2qVtGKmIiEgu4VjY3D634QBzXZC988xJkfTa7DYaCA0Hqs2uSCZYnQQRERGxxq0JkM61SvFx74Y4FrKja53SBIVFEhETh5ebeQmMZoCIiIiko1wj89bpLfMymb3zb7bZvXoOtkyHre+Z2+w2GQq+D6jNrkgGrE6C2NnZYTJlfJOqzjEiIgVXUrJhkdi4fiORZ74MTpMAAbC3M6kNroiIiDVcPaHlKGg+Mm2bXSMJ/lxr3jwrm5Mh9fuASzFbRy2Sq1idBPn2228tHickJLBv3z4WLVpEQEBAlgUmIiJ5y7qQcALWhFoscUnRuVYpZvdpiIO9vpWSm4oVK3bbL1YSE29fSFdEpMBKabNbuR1EnYHgRebuMVfPmccjj8HPr8Cvb0DtR80JkbINbRuzSC5hdRKkR48eafY9+uij1KpVixUrVjB06NAsCUxERPKOdSHhjFgSbFHo9Fb+9cooASJpzJo1y9YhiIjkfSltdluPhz9/gD3/M9cKAXP9kP1LzFuZBuZCqrV6gaOrbWMWsSGTYRgZ3bNa5dixY9StW5dr13JX7+ro6Gg8PDyIiorC3d3d1uGIiOQ7SckG903fmO4MELjZ9nb7xHaq+ZHP6BprPX1mIpIjLhz5t83uspttdlM4F4UG/aDxEChe2SbhiWSHzF5js+RruevXr/Pxxx9Trly5rDidiIjkIUFhkRkmQMDcBjc8Ko6gsMicC0oklwkMDMTPz48mTZrYOhQRKQhKVocHp8MLh8H/I/Cue3Ms7grsnA0fN4TFPeHwWkjS8kMpOKxeDvPf9buGYRATE4OrqytLlizJ0uBERCT3i4jJOAFyN8eJ5EfPPvsszz77bOq3VCIiOcKx8L/tcwfAmd/NS2VCVkFSvHn8+Cbz5l4WGg02H+dWyrYxi2Qzq5MgM2fOtEiC2NnZUbJkSZo1a0axYqo8LCJS0Hi5OWfpcSIiIpLFTCYo19i8dZ4K+5bA3nlw+YR5PPoMbHrL3GrX7yFz7ZAKLczPE8lnrE6CDBo0KBvCEBGRvOrYhau3HU+pCdLU1zNnAhIREZGMuXpCq9HQYhQc22ieHXJ0HWBAcgKEfGPevPzMXWXqPgFObraOWiTLWF0TZN26dWzfvj31cWBgIPXr16dPnz5cvnw5S4MTEZHcbf72MF5bHZL6+L/fF6U8nuLvp6KoIiIiuYmdHVTtAH2Ww/MH4L5x4Fri5nhEKPzwArxfE354ESIO2y5WkSxkdRJk/PjxREdHA3Do0CHGjRtH165dOX78OOPGjcvyAEVEJHcK3PQ3b6wNBeDp1pWY07ch3h6WS168PZyZ068hXWqXtkWIkg+88cYbbN261dZhiIjkb8V8oMMUGBcKvT6H8s1ujt2IgT2fwyfNYUG3f2uKJNguVpF7ZHWL3CJFihASEkLFihV5/fXXCQkJ4euvvyY4OJiuXbty7ty57Ir1rqgVnYhI1jIMg/fXH2X2pr8BGNOhKs+3r4rJZCIp2SAoLJKImDi83MxLYDQDJP/KiWusr68v58+fp3379qxZsyZbXiMn6b5ERPKM8IPmpTKHVkJCrOVYEW9oNMhcdNW9jE3CE/mvzF5jra4J4ujoSGys+R/BL7/8woABAwDw9PRMnSEiIiL5w3+TGk0qFmPaT38yb3sYAC8/WIOn21ROPd7ezkSLysVtFa7kQ2FhYcTFxbFlyxZbhyIiUrCUrgsPfQQd34ADy80JkUt/mceunoMt78DWd6Fmd2gyDCrep0KqkidYnQS57777GDduHK1atSIoKIgVK1YAcPToUcqVK5flAYqIiG2sCwknYE0o4VE3W9u6OtoTeyMJgICHajGwZUUbRScFibOzM507d7Z1GCIiBZNLUWj+DDR7GsK2QNDncORHMJLBSILQ78xbyZrmQqr1nlQhVcnVrK4JMnv2bAoVKsTXX3/NnDlzKFu2LAA//fQTXbp0yfIARUQk560LCWfEkmCLBAiQmgDp39xHCRDJUtu2baNfv360aNGCM2fOAPDFF19YFGMXEREbMpmg0gPw5Jcw5hC0Hg+FvW6OXzgMP74I79cwF1RVIVXJpayuCZLXaO2tiIh1kpIN7pu+MU0C5FalPZzZPrGd6n0UcFl1jf3mm2/o378/ffv25YsvviA0NJRKlSrxySefsHbtWn788ccsjNq2dF8iIvlK4g04/L15qcypnWnHK94PTZ6CGt3A3iHn45MCJbPXWKtnggQHB3Po0KHUx9999x09e/bklVde4caNG3cXrYiI5BpBYZG3TYAAhEfFERQWmUMRSX731ltv8emnn/L555/j4HDzJrlly5YEBwfbMDIREbmtQo5Q51EYsg6e2W4ulurgenP8xDZYORBm1YUt78LVCJuFKpLC6iTI008/zdGjRwE4fvw4Tz75JK6urqxcuZIJEyZkeYAiIpKzImJunwCx9jiROzly5AitW7dOs9/d3Z0rV67kfEAiImI97zrg/yGMOwxdpkPxKjfHYs7CprfgAz/45ik4HQT5e0GC5GJWJ0GOHj1K/fr1AVi5ciWtW7dm6dKlLFy4kG+++eauA5k2bRomk4kxY8ak7jMMg9dff50yZcrg4uLCAw88wB9//HHXryEiInfm5eacpceJ3Enp0qX5+++/0+zfvn07lSpVskFEIiJy11IKqT67B/p/C9W7genfPzuTE8wtd+d1hLltIPgLSLhu03Cl4LE6CWIYBsnJyYC5RW7Xrl0BKF++PBcvXryrIPbs2cPcuXOpW7euxf4ZM2bwwQcfMHv2bPbs2YO3tzcdO3YkJibmrl5HRETuzNvd+ba1PkyYa4I09fXMuaAkX3v66ad5/vnn2b17NyaTibNnz/Lll1/y4osvMnLkSFuHJyIid8PODiq3g95L4fkDcN9YcLnl3iH8AHw/Cj6oCetfg8gw28UqBYrVSZDGjRvz1ltv8cUXX7Blyxa6desGQFhYGKVKlbI6gKtXr9K3b18+//xzihUrlrrfMAxmzZrFq6++Sq9evahduzaLFi0iNjaWpUuXZni++Ph4oqOjLTYREcmcw+HRPD53J0nJ6U9RTUmNTPH3U1FUyTITJkygZ8+etG3blqtXr9K6dWueeuopnn76aUaNGmXr8LJEYGAgfn5+NGnSxNahiIjkvKIVoMPr5qUyPedAmQY3x65fhh0fw0cNYOkT8PevWioj2crqJMisWbMIDg5m1KhRvPrqq1SpYl7r9fXXX9OyZUurA3j22Wfp1q0bHTp0sNgfFhbGuXPn6NSpU+o+Jycn2rRpw44dOzI837Rp0/Dw8Ejdypcvb3VMIiIF0d4TkTzx2U4uxMRTw9uN6b3qUNrDcsmLt4czc/o1pEvt0jaKUvKrt99+m4sXLxIUFMSuXbu4cOECb775pq3DyjLPPvssoaGh7Nmzx9ahiIjYjoMz1O8DwzfDUxuh7pNg7/jvoAFH18GSXjC7CeyeC3H6QluyXpa1yI2Li8Pe3t6iqvudLF++nLfffps9e/bg7OzMAw88QP369Zk1axY7duygVatWnDlzhjJlyqQ+Z/jw4Zw8eZKff/453XPGx8cTHx+f+jg6Opry5curFZ2IyL+Skg2CwiKJiInDy828rGXr0QuM+PJ34hKSaeRTjPkDm+Dh6pDusZoBIinU7tV6+sxERP7j6gUIXgR750P0GcsxRzeo3xuaDocSVW0Tn+QZmb3GFrqbk1+5coWvv/6aY8eOMX78eDw9PQkNDaVUqVKULVs2U+c4ffo0zz//POvXr8fZOePieiaT5c22YRhp9t3KyckJJyenzL0REZECZl1IOAFrQi1a4BZ1cSA6LoFkAx6oXpI5fRvh4mgPgL2diRaVi9sqXCkg4uLi+Pjjj9m0aRMRERGptcdSqE2uiEg+VqQktH4RWo2BIz+YZ4Cc3G4euxEDQXPNW+V20PRpqNoR7OxtGrLkbVYnQQ4ePEj79u0pWrQoJ06cYNiwYXh6evLtt99y8uRJFi9enKnz/P7770RERNCoUaPUfUlJSWzdupXZs2dz5MgRAM6dO0fp0jenXUdERNxV7RERkYJuXUg4I5YE89/pf1euJwDQtGIxPh/QGAd7q1dKityTIUOGsGHDBh599FGaNm162y87REQkn7IvBH49zNu5ENjzORxYAYn/do85ttG8FasITZ6CBv3ApdhtTymSHquTIOPGjWPw4MHMmDEDNze31P0PPvggffr0yfR52rdvz6FDhyz2DR48mBo1ajBx4kQqVaqEt7c3GzZsoEEDc+GcGzdusGXLFqZPn25t2CIiBVpSskHAmtA0CZBbnb58HTv98Sk28MMPP/Djjz/SqlUrW4ciIiK5gXdt8P/QXEx13xII+hyunDSPXT5h7iazaSrUexKaPQMlq9syWsljrE6C7Nmzh88++yzN/rJly3Lu3LlMn8fNzY3atWtb7CtcuDDFixdP3T9mzBimTp1K1apVqVq1KlOnTsXV1dWqZIuIiEBQWKTFEpj0hEfFERQWqeUvkuPKli1r8cWKiIgIYJ7p0fI5aD4S/toAQZ+ZZ4MAJMSa64jsnQ+V2pqTIVU7mVvzityG1UkQZ2fndNvOHjlyhJIlS2ZJUCkmTJjA9evXGTlyJJcvX6ZZs2asX79eN0oiIlaKiLl9AsTa40Sy0vvvv8/EiRP59NNP8fHxsXU4IiKS29jZQ/Uu5u3CUXONkP1LIeGaefz4JvNWzBeaPQ31+4Kzik9L+qzuDjN8+HAuXLjAV199haenJwcPHsTe3p6ePXvSunVrZs2alU2h3h1VYRcRgZ3HLtH78113PG7ZsOaaCSKZllXX2AsXLvD444+zdetWXF1d03Sai4yMvNdQcw3dl4iIZJG4qH+Xysw1L5G5lWMRcyKk6XAoUcUm4UnOy7buMO+99x5du3bFy8uL69ev06ZNG86dO0eLFi14++237yloERHJHqU9nLG3M5GUnH7e2wR4e5hb4IrktN69e3PmzBmmTp1KqVKlVBhVRETuzNkDWjxrXgbz13rYNQfCtpjHblw1L50J+sy8RKbZ01C5Pej6ItxFEsTd3Z3t27ezceNGgoODSU5OpmHDhnTo0CE74hMRkXt06J8oBi/cc9sECMAUfz/s7XRzIDlvx44d7Ny5k3r16tk6FBERyWvs7KH6g+Yt4jDs/gwOLL/ZVeav9eatZA1zwqTuE+DoatuYxaasSoIkJibi7OzM/v37adeuHe3atcuuuEREJAtsPhLByC+Dib2RRM3S7gxs4cOHv/5lUSTV28OZKf5+dKld+jZnEsk+NWrU4Pr167YOQ0RE8jqvmuA/C9pPhn1fmLvKRJ02j134E9aOgV8DoNFgaDoM3MvYMlqxEauSIIUKFcLHx4ekpKTsikdERO5CUrJBUFgkETFxeLmZl7V88/s/vPztIZKSDe6rUoI5/Rri5uzAY43LpzlWM0DElt555x1eeOEF3n77berUqZOmJohqZ4iIiFVcPaHV89D8WTjyA+z6FE7tMI9dvwzbP4AdH0Gth6HZCCjXyLbxSo6yujDqggULWLlyJUuWLMHTM/evHVcBMhHJ79aFhBOwJtRidkcRp0JcjU8EoFeDsrzzSF0cC6llnGStrLrG2v3bzvC/tUAMw8BkMuWrL190XyIiYiNn95mTISHfQHKC5Vj5ZtB8BNTwB3urK0ZILpFthVE/+ugj/v77b8qUKYOPjw+FCxe2GA8ODrY+WhERuSvrQsIZsSSY/2azUxIgXWqX4v3H66nQpORqmzZtsnUIIiKS35VpAL0+g44BsOd/sHc+xF4yj53ebd7cy0Gz4dBwILgUtWm4kn2sToL06NFDN9MiIrlAUrJBwJrQNAmQWx04HUWyAfb6tS25WJs2bWwdgoiIFBRu3tDuNbj/BTi00txVJiLUPBb9D2yYDFtmQIP+0PwZKFbRpuFK1rN6OUxeo2mnIpJf7Tx2id6f77rjccuGNadF5eI5EJEUNFl1jd26dettx1u3bn3X585tdF8iIpLLGIa5te6uOXB0neWYyQ5q+kOLUVC+qW3ik0zL8uUwsbGxjB8/ntWrV5OQkECHDh346KOPKFGiRJYELCIi1omIibvzQVYcJ2IrDzzwQJp9t846zU81QUREJJcxmaDSA+bt4t+w6xPYv9TcYtdIhtDvzFu5ptDiWXNSxM7e1lHLPch0lbwpU6awcOFCunXrxpNPPsmGDRsYMWJEdsYmIiK34eXmnKXHidjK5cuXLbaIiAjWrVtHkyZNWL9+va3DExGRgqJEFej+AYz9w7xkpkipm2P/BMHKgfBRA/OskfgY28Up9yTTM0FWrVrFvHnzePLJJwHo168frVq1IikpCXt7ZcJERHLa4fCo246bAG8PcwtckdzMw8Mjzb6OHTvi5OTE2LFj+f33320QlYiIFFiFi0Pr8dByNBz6GnYGQsQf5rErJ2HdS7BpGjQaCM2eAY+yto1XrJLpmSCnT5/m/vvvT33ctGlTChUqxNmzZ7MlMBERSV9iUjKTvwvhjbWHU/f9t+5pyuMp/n7Y26kqquRNJUuW5MiRI7YOI0sEBgbi5+dHkyZNbB2KiIhkViEnaNAXRvwG/b+Fyu1vjsVHwY6P4MO68O0IOB9quzjFKpkujGpvb8+5c+coWbJk6j43NzcOHjyIr69vtgV4r1SATETyk5i4BEYt3ceWoxcwmeClLjWo4OnKG2tDCY+6WfujtIczU/z96FK7tA2jlfwuq66xBw8etHhsGAbh4eG88847JCQk8Ntvv91rqLmG7ktERPK486GwKxAOfgVJNyzHqnYyzx6peJ+51ojkqMxeYzOdBLGzs+PBBx/Eyckpdd+aNWto164dhQsXTt23atWqewg76+lmQ0TyoqRkg6CwSCJi4vByMy9pOXvlOkMX7eHo+as4O9gx64kGdKntneHxmgEi2S2rrrF2dnaYTCb+e0vSvHlz5s+fT40aNe411FxD9yUiIvnE1QjY/Rns+R/EXbEcK9MQWj2vIqo5LMuTIIMHD87UCy9YsCBzEeYQ3WyISF6zLiScgDWWMzuKF3bkRlIyMXGJeLk5MW9gE+qUS1tHQSQnZdU19uTJkxaP7ezsKFmyJM7O+a+or+5LRETymfirsO8Lc92QqNOWY8V8oeUoqN8XHFxsE18BkuVJkLxKNxsikpesCwlnxJJgMvrFXK6oCytHtKC0hy6kYnu6xlpPn5mISD6VlAB/rIbfPoTzhyzHXItD06eh6TBwVcH67JLZa2ymu8OIiEj2Sko2CFgTmmECBCAx2VDLW8kXPvroo0wfO3r06GyMREREJAvYO0Ddx6DOo3B8kzkZcnyzeSz2EmyeCr/NgkaDoMUodZSxISVBRERyiaCwSIslMOk5Fx1HUFgkLSoXz6GoRLLHzJkzM3WcyWRSEkRERPIOkwkqtzNv4Qfgt4/gj2/BSIKEWNj1CQR9DvWehPvGQvHKto64wFESREQkl4iIuX0CxNrjRHKzsLAwW4cgIiKSvUrXg0fnQfvJsHM2BC+GxDhITjDXEdn/Jfj1gPvGQem6to62wLCzdQAiImKW2WUuWg4j+ZlhGGm6xIiIiORpxXyg67swJsSc8HD6t16FkWyeJfLZ/bDkUTi5w7ZxFhBKgoiI5BLHL1697bgJKO1hbn8rkt8sXryYOnXq4OLigouLC3Xr1uWLL76wdVgiIiJZp0hJ6DAFxoZA+ylQuOTNsb83wIIHYX4XOLoe9IVAtrmr5TBHjx5l8+bNREREkJycbDE2efLkLAlMRKSgSEhK5s21oSzeebNNqAksCqSa/v3fKf5+2NuZEMlPPvjgAyZNmsSoUaNo1aoVhmHw22+/8cwzz3Dx4kXGjh1r6xBFRESyjrMH3D8Omo+AfUvMdUOiTpnHTu2EpY+Bdx1oPR5q+IOd5i5kJatb5H7++eeMGDGCEiVK4O3tjcl082bcZDIRHByc5UHeC7WiE5HcIinZICgskoiYOLzczDM6rsTeYOSXwewOiwTgxU7VqFyyCG+sDbUoklraw5kp/n50qV3aVuGLpJFV11hfX18CAgIYMGCAxf5Fixbx+uuv56v6IbovERGRNJIS4NDXsH0mXDxiOeblZ06G+PVUMuQOMnuNtToJ4uPjw8iRI5k4ceI9B5kTdLMhIrnBupBwAtZYJjZKFHEkORkiY29QxKkQM5+oT0e/UkD6CRPNAJHcJquusc7OzoSEhFClShWL/X/99Rd16tQhLi7/FAPWfYmIiGQoORmO/Ajb3oOz+yzHStYwJ0NqPQx29raJL5fL7DXW6lTS5cuXeeyxx+4pOBGRgmRdSDgjlgSnaX978eoNImNvULKII9+ObJmaAAGwtzPRonJxetQvS4vKxZUAkXytSpUqfPXVV2n2r1ixgqpVq9ogIhERERuws4Oa3WHYJuj7NZRtfHPswp/wzVAIbAYHVkBSou3izOOsrgny2GOPsX79ep555pnsiEdEJF9JSjYIWBPK7abc2dmZqFSySI7FJJLbBAQE8MQTT7B161ZatWqFyWRi+/bt/Prrr+kmR0RERPI1kwmqdoQqHeD4Jtg8HU7vMo9d+gu+HQ5b3oH7X4S6T4D9XZX6LLCs/rSqVKnCpEmT2LVrF3Xq1MHBwcFifPTo0VkWnIhIXhcUFplmBsh/nY+OJygskhaVi+dQVCK5yyOPPMLu3buZOXMmq1evxjAM/Pz8CAoKokGDBrYOT0RExDZMJqjcDiq1hbCtsGU6nPzNPBZ5HL4bCVtnwP0vQL3eYO9w+/MJcBc1QXx9fTM+mcnE8ePH7zmorKS1tyJiS9/tP8Pzy/ff8bgPn6xPj/plsz8gkSyka6z19JmJiMg9ObEdNr8DJ7ZZ7i/mC21fgdqPFNiaIZm9xlo9EyQ/VWgXEcluXm7OWXqcSH7Utm1b+vXrx6OPPoqHh4etwxEREcm9Kt4Hg9bCyR2wZYZ5uQzA5TBYNQy2fQDtXoUa3c0zSSSNe+qxYxgGVk4kEREpMG4kJrP24NnbHmPC3P62qa9nzgQlkgvVqVOH1157DW9vbx555BFWr17NjRs3bB2WiIhI7uXTEgashiHrwbf1zf0XDsOKfvB5W/j7F9Df62ncVRJk8eLF1KlTBxcXF1xcXKhbty5ffPFFVscmIpJnhUdd54m5O/ly96nUff/Nxac8nuLvp+4vUqB99NFHnDlzhu+++w43NzcGDhyIt7c3w4cPZ8uWLbYOT0REJPeq0AwGroEB30O5Jjf3n90HSx6BBV3Ns0YkldU1QT744AMmTZrEqFGjaNWqFYZh8NtvvxEYGMhbb73F2LFjsyvWu6K1tyKSnZKSDYLCIomIicPLzTyjY/fxSzy3bB+Xrt3A3bkQs56sz43EZALWhFoUSS3t4cwUfz+61C5tw3cgcvey6xobFxfHmjVrePvttzl06BBJSUlZdm5b032JiIhkG8OAoz/Dxrfg/CHLscrtod1rULahbWLLAZm9xt5VYdSAgAAGDBhgsX/RokW8/vrrua5miG42RCS7rAsJT5PYcHMuxNW4RAzAr7Q7n/ZrRIXirkD6CRPNAJG8LDuusefOnWP58uUsWbKE4OBgmjRpwu7du7Pk3LmB7ktERCTbJSdD6GrYNNXcUvdWNbpDu0ngVcMmoWWnbCuMGh4eTsuWLdPsb9myJeHh4daeTkQkT1oXEs6IJcH8N4scE5cIQHPf4iwc0gRnh5vVue3tTGqDK5KO6OhovvnmG5YuXcrmzZupVKkSffr0Yfny5VSpUsXW4YmIiOQtdnZQuxfUfAgOLofN0yHq3yXaf66FIz9Cw4HmbjJFvGwbqw1YXROkSpUqfPXVV2n2r1ixgqpVq2ZJUCIiuVlSskHAmtA0CZBbnYy8hoP9PdWeFikwSpUqxauvvkqtWrXYsWMHR44cYcqUKUqAiIiI3Av7QtCgHzz3O3R9D4p4m/cbyfD7AvioAWx9F27E2jbOHGb1TJCAgACeeOIJtm7dSqtWrTCZTGzfvp1ff/013eSIiEh+ExQWabEEJj3hUXEEhUVq5odIJnz33Xd06NABOzslDkVERLJcIUdoOgzq94Vdn8D2mXDjqnnb+BbsXWBeIlP3CfMsknzO6nf4yCOPsHv3bkqUKMHq1atZtWoVJUqUICgoiIcfftiqc82ZM4e6devi7u6Ou7s7LVq04KeffkodNwyD119/nTJlyuDi4sIDDzzAH3/8YW3IIiJZKiLm9gkQa48TKeg6deqkBIiIiEh2c3SF1i/C6H3QeAiY/l22HX0GVj8Dc9vA8fzflc3qmSAAjRo1YsmSJff84uXKleOdd95Jne66aNEievTowb59+6hVqxYzZszggw8+YOHChVSrVo233nqLjh07cuTIEdzc3O759UVE7saNxORMHefl5pzNkYiIiIiIWKmIF3SfCU2fhl+mwNF15v3nDsLih6BaF+j4BpSsbts4s0mmusNER0enVleNjo6+7bH3Wunc09OTd999lyFDhlCmTBnGjBnDxIkTAYiPj6dUqVJMnz6dp59+OlPnUxV2EckqhmGwcu8/TPruEPGJGf/qNAHeHs5sn9hO3V8kX9M1NvMCAwMJDAwkKSmJo0eP6jMTEZHc4/gWWP+aOQmSwmQPjQbBAy9DkZI2C80aWdoi197envDwcLy8vLCzs8NkSntTbxgGJpOJpKSkuwo4KSmJlStXMnDgQPbt24ezszOVK1cmODiYBg0apB7Xo0cPihYtyqJFi9I9T3x8PPHx8amPo6OjKV++vG42ROSeXItP5LXVIXy77wwANUu7cTg8BhNYFEhN+e04p19DutQundNhiuQoJUGsp89MRERypeRkOLgCfn0DYs7e3O/kAZ3fggb9IZ08QG6SpS1yN27ciKenJwCbNm3Kmgj/dejQIVq0aEFcXBxFihTh22+/xc/Pjx07dgDmivG3KlWqFCdPnszwfNOmTSMgICBLYxSRgiMp2SAoLJKImDi83Jxp6uvJ0fMxPLs0mOMXrmFvZ2Jcx2qMaFOZ9aHnCFgTalEk1dvDmSn+fkqAiGRSQkICnTp14rPPPqNatWq2DkdERKRgsrOD+r3BrwfsCoTts8yFU+Oj4PvnIGQVPPQRFK1g60jvWaaSIG3atEn9b19fX8qXL59mNohhGJw+fdrqAKpXr87+/fu5cuUK33zzDQMHDmTLlpvFWNJ7nfRmoqR4+eWXGTduXOrjlJkgIiJ3si4kPE1Sw93ZgdgbiSQmG3i7O/NxnwY0qWhOCnepXZqOft5pkiZaAiOSeQ4ODoSEhNz22i4iIiI5xNEVWo+HhgNhw2Q4sMy8//gm+KQFdHgdGg/N011krI7c19eXCxcupNkfGRmJr6+v1QE4OjpSpUoVGjduzLRp06hXrx4ffvgh3t7mHsbnzp2zOD4iIiLN7JBbOTk5pXabSdlERO5kXUg4I5YEp2l9Gx2XQGKyQe0y7vz4/P2pCZAU9nYmWlQuTo/6ZWlRubgSICJ3YcCAAcybN8/WYYiIiEiKIl7w8KfQ92twL2ved+Mq/PgiLOoOl47ZNr57YHV3mIxmYly9ehVn53vvhGAYBvHx8fj6+uLt7c2GDRtSa4LcuHGDLVu2MH369Ht+HRGRFEnJBgFrQrldgaRL127g4eKQYzGJFCQ3btzgf//7Hxs2bKBx48YULlzYYvyDDz6wUWQiIiIFXNWOMHKXeVbI7wvM+07+BnNaQbtXoflIsLO3bYxWynQSJGWJiclkYtKkSbi6uqaOJSUlsXv3burXr2/Vi7/yyis8+OCDlC9fnpiYGJYvX87mzZtZt24dJpOJMWPGMHXqVKpWrUrVqlWZOnUqrq6u9OnTx6rXERG5naCwyDQzQP4rPCqOoLBIWlQunkNRiRQcISEhNGzYEICjR49ajGmZjIiIiI05u4P/LKj1sLk+yJWTkHjd3FHmj9XQIxC8atg6ykzLdBJk3759gHmmxqFDh3B0dEwdc3R0pF69erz44otWvfj58+fp378/4eHheHh4ULduXdatW0fHjh0BmDBhAtevX2fkyJFcvnyZZs2asX79etzc3Kx6HRGR24mIuX0CxNrjRMQ6WV10XURERLJBpTYwcif8+ibs/hQw4Mxe+Ox+aDMBWo0B+9w/czpTLXJvNXjwYD788MM8U2tDrehE5E4+/OUvZv5y9I7HLRvWXDNBRG6ha6z19JmJiEi+cGoXfPcsXPr75j7vOtDjEyhd1yYhZWmL3FstWLDgngITEcktYm8k8ubaUJYF3b6zlQlz69umvp63PU5E7t6ePXtYuXIlp06d4saNGxZjq1atslFUIiIikq4KzeGZ7bD5HdjxERjJcO4QfN4W7hsHrV+EQk62jjJdVidBQDcqIpK3JCUbadrYHg6PZvTyfRy/cA2TCTrULMUvoecBLAqkplQjmOLvp84vItlk+fLlDBgwgE6dOrFhwwY6derEX3/9xblz53j44YdtHZ6IiIikx8EFOgaA30Pw3SiICIXkRNg6Aw6vMdcKKdfI1lGmYXUSRDcqIpKXrAsJJ2BNqEXhUzfnQsTeSCQpGbzdnfng8Xq0rFIi3WO9PZyZ4u9Hl9qlbRG+SIEwdepUZs6cybPPPoubmxsffvghvr6+PP3005QurX97IiIiuVrZRjB8C2x7H7a9Z06EXDgM8zpAi1HQ9hVzwiSXsLomSN26dXn66adTb1QOHDhgcaMSEBCQXbHeFa29FSm41oWEM2JJcIatb+uXL8qCQU0oVvhmoef0Zo1oBohI+rLqGlu4cGH++OMPKlasSIkSJdi0aRN16tTh8OHDtGvXjvDw8CyM2rZ0XyIiIvnauUPmWiHhB27u86xsnhXi0yJbXzqz11g7a0987NgxunXrBoCTkxPXrl3DZDIxduxY5s6de/cRi4hkoaRkg4A1oRkmQADOR8fh7mJZwdrezkSLysXpUb8sLSoXVwJEJAd4enoSExMDQNmyZQkJCQHgypUrxMbG2jI0ERERsYZ3HXhqI7SfAvb/1gSJPAYLHoQfJ0D8VdvGx10kQXSjIiJ5QVBYpMWylvSER8URFBaZQxGJSEbuv/9+NmzYAMDjjz/O888/z7Bhw+jduzft27e3cXQiIiJiFftCcP84eGYblGv6704Dgj6DOS3g+GZbRmd9TZCUG5U6deqk3qhs3LiRDRs26EZFRHKNiJjbJ0CsPU5Ess/s2bOJizP/W3z55ZdxcHBg+/bt9OrVi0mTJtk4OhEREbkrJavDkHWw+zP49Q1IvA5XTsHiHtBwIHR6E5w9cjwsq2uCREZGEhcXR5kyZUhOTua9995j+/btVKlShUmTJlGsWLHsivWuaO2tSMFzLT6R55fv45fDEXc8dtmw5rSoXDwHohLJf3SNtZ4+MxERKZAij8P3o+HEtpv73MpAv2+glF+WvERmr7FWzQRJTExkzZo1dO7cGQA7OzsmTJjAhAkT7i1aERErZVTAdM+JSF5ceYCTl26/PM+EufNLU1/PnAlYRCxER0dn+lglC0RERPI4z0ow4Hv4fQFsmAw3rpo7xnj65ngoViVBChUqxIgRIzh8+HB2xSMickfptrJ1d6JOuaL8cvg8hgFli7rwWKNyfPjrXwAWBVJTSp1O8fdT4VMRGylatCgm0+3//RmGgclkIikpKYeiEhERkWxjZwdNhkLVTrB2rLluiA1a51pdE6RZs2bs27cPHx+f7IhHROS2Mmp7ey46nnOh5wF4rFE5Jvn74e7sQI3SbmkTJh7OTPH3o0vt0jkYuYjcatOmTbYOQURERGyhaHno97XNXt7qJMjIkSN54YUX+Oeff2jUqBGFCxe2GK9bt26WBScicqvMtL0t5urAO4/UTZ3h0aV2aTr6eae7dEZEbKdNmza2DkFEREQKIKuTIE888QQAo0ePTt1nMpk0ZVVEsl1m2t5ejk0gKCzSotipvZ1JxU9FcrkrV64wb948Dh8+jMlkws/PjyFDhuDhkfNV40VERCT/sjoJEhYWlh1xiIjckdreiuRPe/fupXPnzri4uNC0aVMMw+CDDz7g7bffZv369TRs2NDWIYqIiEg+YXUSRLVARMRWYuISM3Wcl5tzNkciIllp7NixPPTQQ3z++ecUKmS+NUlMTOSpp55izJgxbN261cYRioiISH5hdRJk8eLFtx0fMGDAXQcjIpKeuIQkPthwlM+3Hr/tcWp7K5I37d271yIBAuaOdBMmTKBx48Y2jExERETyG6uTIM8//7zF44SEBGJjY3F0dMTV1VVJEBG5K0nJRrrFS3cfv8TEbw5y4lIsAE0rFmPPicuA2t6K5Bfu7u6cOnWKGjVqWOw/ffo0bm5uNopKRERE8iOrkyCXL19Os++vv/5ixIgRjB8/PkuCEpGCZV1IeJo2tqXcnaju7cbWoxcB8HZ3Zmqv2rSrUSrd49X2ViTveuKJJxg6dCjvvfceLVu2xGQysX37dsaPH0/v3r1tHZ6IiIjkI1YnQdJTtWpV3nnnHfr168eff/6ZFacUkQJiXUg4I5YEp2l7ez46nvPR8QD0aVaBlx6sgbuzA6C2tyL5zXvvvYfJZGLAgAEkJppr/zg4ODBixAjeeecdG0cnIiIi+UmWJEEA7O3tOXv2bFadTkQKgKRkg4A1oWkSILcqXtiRN3vUTpPgUNtbkfzD0dGRDz/8kGnTpnHs2DEMw6BKlSq4urraOjQRERHJZ6xOgnz//fcWjw3DIDw8nNmzZ9OqVassC0xE8r+gsEiLJS3puXTtBkFhkUp4iBQArq6u1KlTx9ZhiIiISD5mdRKkZ8+eFo9NJhMlS5akXbt2vP/++1kVl4gUABExt0+AWHuciOQtQ4YMydRx8+fPz+ZIREREpKCwOgmSnJycHXGISAETl5DE1qMXMnWsl5tzNkcjIrawcOFCfHx8aNCgAYZxu4VxIiIiIlnD6iTIG2+8wYsvvphmne7169d59913mTx5cpYFJyJ5W0Ztb3f8fZFXV4cQdvHabZ9vwtz1pamvZ84ELCI56plnnmH58uUcP36cIUOG0K9fPzw99e9dREREso/JsPKrF3t7e8LDw/Hy8rLYf+nSJby8vEhKSsrSAO9VdHQ0Hh4eREVF4e7ubutwRAqM9NrYerk5UbFEYYLCIlMf96hfhv9tCwOwKJCaUgZ1Tr+GansrkktlxTU2Pj6eVatWMX/+fHbs2EG3bt0YOnQonTp1wmTKfx2fdF8iIiKSPTJ7jbV6JohhGOnelBw4cEDf3ogIkHHb24iYeCJizG1vB7bw4YXO1XF3dqCRT7E0CRNvD2em+PspASKSzzk5OdG7d2969+7NyZMnWbhwISNHjiQhIYHQ0FCKFCli6xBFREQkH8l0EqRYsWKYTCZMJhPVqlWzSIQkJSVx9epVnnnmmWwJUkTyjsy0vS1RxJHJ/rVS2952qV2ajn7e6S6dEZGCI+U+wzAM1SATERGRbJHpJMisWbMwDIMhQ4YQEBCAh4dH6pijoyMVK1akRYsW2RKkiOQdmWl7e/Fq2ra39nYmtcEVKYBuXQ6zfft2unfvzuzZs+nSpQt2dna2Dk9ERETymUwnQQYOHAiAr68vLVu2xMHBIduCEpG863y02t6KSOaMHDmS5cuXU6FCBQYPHszy5cspXlzJUBEREck+VtcEadOmTep/X79+nYSEBItxFfkSyd8y6vgCEHImik82/52p86jtrYh8+umnVKhQAV9fX7Zs2cKWLVvSPW7VqlU5HJmIiIjkV1YnQWJjY5kwYQJfffUVly5dSjOe27rDiEjWSa/jS2kPZ17sVI19p6/w5e5TGIa5s0tGNUHU9lZEUgwYMCBfdoARERGR3MvqJMj48ePZtGkTn3zyCQMGDCAwMJAzZ87w2Wef8c4772RHjCKSC2TU8SU8Ko4XVh5MffxQvTK0qFycV1YdAtJvezvF309FT0WEhQsX2joEERERKWCsToKsWbOGxYsX88ADDzBkyBDuv/9+qlSpgo+PD19++SV9+/bNjjhFxIYy0/GlkJ2JxUOa0rJKCQCKuTqo7a2IiIiIiOQqVidBIiMj8fX1Bcz1PyIjIwG47777GDFiRNZGJyK5QmY6viQmGxbT2tX2VkTyq4cffpjNmzfTvn17vv76a1uHIyIiIlawuvdcpUqVOHHiBAB+fn589dVXgHmGSNGiRbMyNhHJJTLbyeW/x6W0ve1RvywtKhdXAkRE8oXRo0ezePFiW4chIiIid8HqJMjgwYM5cOAAAC+//DKffPIJTk5OjB07lvHjx2d5gCJiW4Zh8Nf5mEwdq44vIlIQtG3bFjc3N1uHISIiInfB6iTI2LFjGT16NGC+CTh8+DDLli37f3v3Ht90ff0P/JWmuTRtk7b0mt4o90u5g1CncpGrU0F8KN5hQ1EHOobbV2EqMDcuOhH8rjBABTaH4HcIQ0Sw27jtBwK2K4JAESiUS0vpLb3m0uT9+yNNIPSSFNrmk/T1fDzySD6XfPI+VOnh5P05b2RlZeGXv/xliw+QiFqX1SZw6Fwx/pF9BYfOFcNqu9H5I/NiKR5ZeRB/2nOuyWvIYF8lhiu+EJG37d+/Hw899BD0ej1kMhm2bdtW75yVK1ciJSUFarUagwYNwoEDB9p+oEREROQVze4Jcqvk5GQkJyff1nsXL16ML774AqdPn0ZQUBDuvvtuLF26FN27d3eeI4TAwoULsWbNGpSWlmLo0KFIT09H796973ToRO1eY0vezhrZBd/mluDLY1cBABqlHKN7RuPLY/kAuOILEUlXVVUV+vXrh5/97Gd49NFH6x3fvHkzZs+ejZUrV+InP/kJVq9ejQkTJuDkyZNISkoCAAwaNAgmk6nee7/55hvo9fpWj4GIiIhaj8czQQ4fPoyvv/7aZd9f/vIXpKSkIDo6GjNmzGgwYWjKvn37MHPmTHz77bfIyMhAbW0txo4di6qqKuc57777LpYtW4Y//elPOHr0KGJjYzFmzBhUVHg2PZ+IGuZY8vbWhqf5BiN+u+0Evjx2FTIZMGVwIvb+egQ+fHIgVj0zELE611teYnVqrHpmIFd8ISJJmDBhAn7/+99j8uTJDR5ftmwZpk+fjueffx49e/bE8uXLkZiYiFWrVjnPyczMxIkTJ+o9bqcAYjKZUF5e7vIgIiIi7/F4JsiCBQswYsQITJgwAQBw/PhxTJ8+HdOmTUPPnj3x3nvvQa/XY8GCBR5/+K5du1y2161bh+joaGRmZuK+++6DEALLly/Hb3/7W2cys2HDBsTExGDjxo148cUXPf4sIrrBkyVvlfIA/P3lNPRNCHPu44ovROTLzGYzMjMz8cYbb7jsHzt2LA4ePNgqn7l48WIsXLiwVa5NREREzefxTJDs7Gzcf//9zu1NmzZh6NChWLt2LebMmYMPP/zQuVLM7TIYDACAiAh7X4Hc3FwUFBRg7NixznNUKhWGDx/eaLLCb1yI3PNkyVuz1YYqk7Xefq74QkS+qqioCFarFTExMS77Y2JiUFBQ4PF1xo0bh8ceeww7d+5EQkICjh492ui5c+fOhcFgcD4uXbp02+MnIiKiO+fxTJDS0lKXpGHfvn0YP368c3vIkCF39ItdCIE5c+bgnnvuQWpqKgA4E5KGkpWLFy82eB1+40Lk3u0ueUtE5A9kMtfirRCi3r6m7N692+NzVSoVVCqVx+cTERFR6/J4JkhMTAxyc3MB2KeTZmVlIS0tzXm8oqICCoXitgcya9YsfP/99/jss8/qHWtOssJvXIiaXvHlbGEF/nqo4SLirbjkLRH5k8jISMjl8nqzPgoLC+t94UJERET+yeOZIOPHj8cbb7yBpUuXYtu2bdBoNLj33nudx7///nt07tz5tgbxyiuvYPv27di/fz8SEhKc+2NjYwHYZ4TExd1outhUssJvXKi9a2zFl1fv74pjl8rw+XeXYGuqGQjsK77EcslbIvIzSqUSgwYNQkZGBh555BHn/oyMDEycONGLIyMiIqK24nERxNFpffjw4QgJCcGGDRugVCqdxz/55BOX3h2eEELglVdewdatW7F3716kpKS4HE9JSUFsbCwyMjIwYMAAAPZZKPv27cPSpUub9VlE7YFjxZdbaxz5BiPmfnHcuT2mVwzSOnXAOztOAuCSt0TkPyorK3H27Fnndm5uLrKzsxEREYGkpCTMmTMHzz77LAYPHoy0tDSsWbMGeXl5eOmll7w4aiIiImorHhdBoqKicODAARgMBoSEhEAul7sc/7//+z+EhIQ068NnzpyJjRs34h//+AdCQ0Od01N1Oh2CgoIgk8kwe/ZsLFq0CF27dkXXrl2xaNEiaDQaPPXUU836LCJ/58mKLwq5DJ9OH4qhnToAAPRh6nqzRmJ1asx/qBeXvCUin/Tdd99h5MiRzu05c+YAAKZOnYr169djypQpKC4uxu9+9zvk5+cjNTUVO3fuRHJysreGTERERG1IJoRwMzG+FT+8kb4e69atw7Rp0wDYZ4ssXLgQq1evRmlpKYYOHYr09HRn81R3ysvLodPpYDAYoNVqW2roRJJz6Fwxnlz7rdvzPnthGNI6d3BuW22CS94S0W3h71jPpaenIz09HVarFWfOnOGfGRERUQvzNC/xahGkLTBBo/Zia9Zl/OrzY27PW/FEf0zsH98GIyIif8ffsc3HPzMiIqLW4envWI9vhyEi72psxobVJvDlsatYuivHo+twxRciIiIiImqvWAQh8gENrfgSq1XjoX5x2JtzHT8WVgIAZDKgsbldXPGFiIiIiIjaOxZBiCSusRVfCsqNWHsgFwCgVQdixn2dEB+uwZzN2QC44gsREREREdGtWAQhkjBPVnwJUQVi729GIiLYvmR1kCKAK74QERERERE1gEUQIgk7klviUsxoSKWpFjkFFc4VX8anxmFMr1iu+EJERERERHQLFkGIvKip5WltNoGMkwUeXaewwrVQIg+QuSyDS0RERERERCyCEHlNQ81O43RqvPXTXrDYbEjfcxZnrlV6dC2u+EJEJG3p6elIT0+H1Wr19lCIiIjaNZkQja0l4R88XSuYqC011uz0ViFKOSADqkzWBs91rPjyn9dH8XYXImpz/B3bfPwzIyKi9s5osWLn8XxEhqhwX7eoFruup79jOROEqI150uxUJgNm398V036SgkPnivDyp1mQgSu+EBERERGRbzp5tRybj+Zh63+voNxYi7s6RrRoEcRTLIIQtTFPmp0KAdyV0gG6IAXGp8Zh1TMDueILERERERH5lEpTLb48dhWbjuTh2GWDy7EjF0qQW1SFlMjgNh0TiyBELaipRqcO56571ufj5manXPGFiIiIiIh8gRAC2ZfKsOnIJXz5/VVUm137YakVAXiwrx5PDElExw6aNh8fiyBELaSxRqeO2RoXi6uw9sB5bD5yyaPr3drslCu+EBERERGRVJVVm7Htv1ew6eglnC6oqHe8t16LJ+5KwsT+emjVCi+M0I5FEKIW0Fij0wKDES99moWByWHIziuDre4EhVwGi7XhriCOZqd3pUS06piJiIiIiIjuhBACh3NLsOlIHnaeKIC51uZyPEQViIn99XjyriSkxuu8NEpXLIIQ3aGmGp069mVdLAMADO8WhZeGd0ZZtRm/+FuWyzkAm50SEREREZH0Xa8wYUvWZWw+egm5RVX1jg9MCsMTdyXhwb5x0CilVXaQ1miIfJAnjU4BYOmjfTBlSJJzm81OiYiIiIjIV1htAgd+vI7NRy8h4+Q11NpcvwYO0ygweUACnrgrEd1iQr00SvdYBCFyw12z05sbmDZFrZC7bLPZKRFR+5Geno709HRYrVb3JxMREUlIvqEGnx+9jM+/u4QrZTX1jt/duQOmDEnEuN6x9f7NI0UsghA1wV2z07OFFfhH9lWPrnVro1OAzU6JiNqLmTNnYubMmSgvL4dOJ417oomIiBpTZap1zvrYd+Y6bpn0gcgQFR4bnIApgxPRsY2XuL1TLIIQNcJds9NecVqczC93ex02OiUiIiIiIqk7f70Se3KuY29OIQ6fL4HZ6trkVCYDRnSLwpQhSbi/ZzQU8gAvjfTOsAhC1ABPmp06CiBje8UgNV6HDzLOuBwH2OiUiIiIiIikyWix4nBuCfacLsTenEJcKK5u8Dy9To3HhyTi8cGJ0IcFtfEoWx6LIEQN8LTZ6Yop/TFxQDwAoFtMCBudEhERERGRZF0urbbP9jhdiP93rghGi63B8+LDgjCiexTG9IrBvV2j/OoLXRZBqF1y1+z0WrlnzU5x098FbHRKRERERERSkm+oweHzJfj2fDG+PV/c6GyPwAAZBncMx8ju0RjZIxpdo0Mgk/nnv2NYBKF2p6lmp2mdIvH3rMtYu/+8R9e6tdkpG50SEREREZG3eFr0AICoUBVGdo/CyO7R+EnXSGjVijYcqfewCELtSmPNTvPrmp0q5DJYrPajMqDBniCOY2x2SkRERERE3iKEwJWyGnx3odSjoodSHoD+iWG4r1skRnSPRq84LQLa4ax1FkGo3Wiq2amDxSrQPSYEz93dEcHKQPxqczYANjslIiIiIiLvMlqs+P6yAf/NK0VWXimy8spwvcLU6PlKeQD6J4VhWKcOGNYpAgOTwqFWyNtwxNLEIgi1G542O13wcG+kdY4EAKgVAWx2SkREREREbUoIgculNfZix8VS/PdSGU5eLUetrfGvdBVyGQYkhmNYpwgM69QBA5LCEaRk0eNWLIKQX3DX6NRosWLn8XyPrlV4UzWVzU6JiKglpKenIz09HVar1dtDISIiCSosN+LEVQOOXy7H8SsGZF8qQ1Fl47M8ACBUFYj+SWEYkBjGokczsAhCPq+pRqcdI4Ox6cglbP3vFRhqLB5dj81OiYiopc2cORMzZ85EeXk5dDqdt4dDREReIoTAtXITjl8x4ETd4/gVg8sXsY3pGh2CAUlhGJgUjgFJ4egSHcIvZ28DiyDk09w1Or2ZXqdGubEWlabaBq/FZqdERERERNRSrDaBvJJqnM4vx6l8+wyP41fK3c7wAIBQdSAGJIVjYFIYBiSFo39iGHRB7WP1ltbGIgj5LE8anQLA+N4xeOKuJNzbNQoZJwvwcl1xhM1OiYiIiIioJRiqLThdUI7TBRU4XVCOk/kVOFNQgRqL+9sgtepApMbr0Cdeh951z8kRmna5cktbYBGEJMtdnw9PG51OvTvFeTvL+NQ4rHpmIJudEhERERFRsxktVpy/XoUfCytw5loFTudX4HRBBa6U1Xj0fl2QAn3idUiN1yE1Xos+8TokRWggk7Hg0VZYBCFJaqrPx4ju0dj9QwFW7T3n0bUKK1wLJWx2SkRERERETTFarDhbWImzhZV1BQ/764vFVWhigRYXyR006BEbih6xWvSMC0VvvQ4J4UEseHgZiyAkOe76fKgDA2CstXl8vVsbnQJsdkpEREREREBJlRnnr1fi/PUqnCuqxLnCSvxYWIm8kmoID4sdoapA9IgLRc84LXrEatEjLhTdY0IRrOI/t6WIPxWSFE/6fBhrbdDr1Jg8KAGbj1xCUaWpwfPZ6JSIiIiIiMy1NuSVVOHc9Sqcv16F89crce56Jc4XVaGs2rMVJAEgSCFHl+gQdI0OQdeYUHSNDkGPuFDEh3F2hy9hEYTaVEv1+fjjY/1wd5dIpOq1ePnTLMjARqdERERERO2VxWrD5dIaXCiuwoWiKlwsrkZuURUuFlfhUmkNrJ7ewwJAo5Sja3QIukSHomtMCLrFhKBrtL3YwWalvo9FEGozTfX5uK9bFDJOXsOa/ec9utb1umWl2OiUiIiIiKh9MFqsuFxag7ySKlwoqrYXPIqrcbG4CpebWegAgFitGp2igtE5KgSdooLRKSoEnSKDWezwcyyCUJtw1+dDKQ+A2Xp7fT7Y6JSIiIiIyPcJIXC9woS8kmqXx+WSGuSVVKOg3P2M8VtplHKkRN4ocDiKHimRwezZ0U7xp06tzpM+H2arDUkRQXi4nx6bjl5CcaW5WX0+2OiUiIiIiEjahBAoqTLjcmkNrpTV4HJpNS6X1uByaQ0ulVTjUmk1jBbPvxh1CFEFIrmDBh0jg9GxgwbJHYLRsUMwOkZqEBWiYr8OcuHVIsj+/fvx3nvvITMzE/n5+di6dSsmTZrkPC6EwMKFC7FmzRqUlpZi6NChSE9PR+/evb03aHLhrscHABw+X+xRn4+lj/ZFWudIpMbr2OeDiIiIiMjH2GwC1ytNuFJWg6tlNbhSV+C4udhRY7He1rUjQ5RIjNAgKUKDxPAbBY+OkcHoEKxkoYM85tUiSFVVFfr164ef/exnePTRR+sdf/fdd7Fs2TKsX78e3bp1w+9//3uMGTMGOTk5CA0N9cKI6WZN9fgY2ysW310sxc7j+fjiv5c9ul5hBft8EBGRf0pPT0d6ejqs1ttL/omIpKDSVGsvbpTVIL/MiKuOYkdZDa4aalBgMMJibV5fDgdlYACS6oocSREaJEZokBgehKQO9qIHb12hliITwtPVj1uXTCZzmQkihIBer8fs2bPx+uuvAwBMJhNiYmKwdOlSvPjiiw1ex2QywWQyObfLy8uRmJgIg8EArVbb6nG0F431+HDQqgNRbqxt1jU/e2GYyy0tnswyISIi7ykvL4dOp+Pv2GbgnxkRSVW50YICgxH5BiMKDDXINxiRX2ZEfvmN7Ypm5vc3UwUGID48CAnhGiSEB9U9bryODFaxGSndEU9/x0q2nJabm4uCggKMHTvWuU+lUmH48OE4ePBgo0WQxYsXY+HChW01TL/krvjgSY+PcmMtQlRyjO0diwm9Y/HWP07gWrmJfT6IiIiIiNqQ1SZQVGlCgcGIgnIjCsvtzwUGE645XxtRabr9AgcA6IIU0IcFIT5MDX1YEOJ0QS7FjsgQ3rJC0iDZIkhBQQEAICYmxmV/TEwMLl682Oj75s6dizlz5ji3HTNByDNN3eLiuA3lPz9e96jHx8qnB+K+btEAAKsQ7PNBRERERNRCbDaB0mozrpWbUFhhRGGFCYXl9ucCgxHXyo24Vm7C9UpTs5eOvZUqMABxOjVideq6QkcQ9HWP+DA14nRBvF2FfIbk/0u9tVoohGiygqhSqaBSqVp7WH6psVtcCgxGvPxpFp5LS0a+wYg9OYUeXa+02uJ8zT4fRERERETuWaw2FFWacL3ixqOwwl7osBc8TLheV+yovcPiBgCoFQHQ64IQ6yhy1L2Ou2k7TKPgLA7yG5ItgsTGxgKwzwiJi7vxD+TCwsJ6s0PIvTu5xcWxb8OhxmfgNCQ6VO2yPT41DmN6xbLPBxERERG1K9a6WRtFlSYUVdQ9V7oWOa5X2GdtlFSZW+QzZTIgMkSFWK0aMVo1YrR1r3Vq575YrRraoEAWOKhdkWwRJCUlBbGxscjIyMCAAQMAAGazGfv27cPSpUu9PDrf4sktLkdySzy6xeXRgfGYdncKXvjrd7hmMDarxwfAPh9ERERE5B9MtVaUVJlRXGkvajie7Q+zs8hRVGlGSZUJLTBpA4C9uNEhWInoUDWitSrE1D1Ha9WIDlUhOlSFGK0aUaEqKOQBLfOhRH7Eq0WQyspKnD171rmdm5uL7OxsREREICkpCbNnz8aiRYvQtWtXdO3aFYsWLYJGo8FTTz3lxVH7Fne3uLz/eD+oAuVYfzDXo+vd1y0KfRJ0WPBQL/b4ICIiIiK/UWu1oaTajJIqM0oqzSiusr8urjKjuK7IUVx1o8BxJyulNEQZGIDoUBWiQlWICql7DlUhOlSNyBAlYrT2YkdkCIsbRHfCq0WQ7777DiNHjnRuOxqaTp06FevXr8f//M//oKamBr/4xS9QWlqKoUOH4ptvvkFoaKi3huxTPLnFZc7nx5p1TcctLuzxQURERERSVm2uRUmVGaVVFpRUm1FaV9RwFDZKqkw3FTnMMNRY3F+0mRRyGSJDVHUPpf059MZ2dKjaWezQqnlbClFbkAkhWmhiljR5ulawr3HX4wMADp0rxpNrv3V7Lb1OjQf6xuGLrCsorTI3eYvLf14fVa+XCHt8EBG1T/76O7Y18c+M6PbUmK3OQkZptRml1RaU1c3aKKu22IsdddulVWaUVJthtNhaZSyh6kBEhqjQIViJDiFKdAhRITLY/tzBUegIsc/mYL8Norbj6e9YyfYEoca56/EhhEDOtQpsOprn0fVen9ADE/vHY3ByeLNvcWGPDyIiIiLyVK3VhrIaewGjrNriLGaUVVtQVmMvbhiqLc5Ch6PoYaptnYIGAISoAhERrEREsNJZ2IgIthc5Iuq2I+sKHBHBSqgC5a02FiJqfSyC+Jimeny89GkW0jp1wPmiSlwrN3l8Td7iQkRERETNYbRYUVZtgcFR0KixvzbUFTPsRQ3XbUO1BRWmlu2jcavAABnCg5WI0CgRHqxARLAS4Rqls8hx83aHEPtrtYJFDaL2hEUQCWmJZWwPnS8GYF/v+66OEfjvpbJGmzY1tIoLl7ElIiIi8n9CCBgtNhhqLCg32gsY5Y5CRo2juOG6r+ym1+ZWnJnhoJQHIExjL2SEaRQI1ygRHqxEuON1XaEjTHOjsMG+GkTkDosgEuHuFherTWDzkTyPlrGd90APPJfWEWqF3DlzBOAtLkRERET+xFxrQ7nRXqgoN9bWPVtQXlPr3G+oO+YocpTfVPSwWNumNaA8QIawIAV0QQpnMUNX9xwWpEBYsP05XGMvdtgfSgQr5SxoEFGLYxGklXnSOLSxW1zy625x6Z+oQ25Rtccdq2O0aue0Pt7iQkRERCQ9VptApakWFUYLKoy1dQ97gcKx7SxuOPfZixgVdftaq/FnYzRKOXR1xQxtkMKlsBGmUd54HWQvZji2Q1ScnUFE0sEiSCtyN7sDaPoWF4fsSwYAQJAiADUe/LJz9Phw4C0uRERE3pWeno709HRYrVZvD4VagKnW6ixUVNYVJypMjm17kaLSVIty481FjhvFjUqT/eENoapAaJ2FjED7s1rhUtxwKXTUFTO0agWUgQFeGTMRUUviErmtpLHZHY6yw8qnB6JrTAg2Hs7DJ//vgtvrvTOxNx4fnIgRf9yLAoOxWcvYEhERtQYu99p8/DPzHptNoMpciyqTFZUmCypNVlTeVJCoNFrsMzNMtTf2G2/ZrntttrbtDIybhagCoVUHIlR9o5ChVdsLFlp1YN2z6/5Qtb3YEapWMEckIr/FJXJbUUs0MJ25MQu2ZpSftEEKqBRyzH+oV7OXsSUiIiLyReZaG6rqig/2AkYtKk3WG/tMrvuq6ooYN/bfKFxUmb0/CydIIUeoOrDuYS9OaOueHfscBQ5H8cJxXBtkv62EeR4R0Z1hEaSZPLnF5UhuidsGpjYBKOQydIsJxQ9Xy91+LpexJSIiIikTQqDGYkWlqRbVprpns9VZvKg23XhdZXYULeqezbUu73Oc11aNO91RyGUIUQUiRB2IEFVd0UJlL1yE1BUtHDM0QtSBCFUp6s69UcQIUQdCIeftJERE3sYiSDM0dotLQV0D0xfuTYFNABknr3l0vaWT+2LigHjcs/Tfbm9x4TK2RERE1BKEEDBbbagxW1FltqL6pmJFtcmKaot9X5XZihpzrfOcKrMV1XW3k7g8O65hsUJKN1nLZECI0l58CFbZCxKOR3BdASNYJa97rUCo43hd8SJUfWNbFSj3djhERNRCWATxkCe3uKw9kNusa8aFBUEeILutW1y4jC0REREBwNfH83G6oAI1FntRotpsdRYzahzb5hvHasxW1Dbnntw2pFYEIEQVCI3SUbiwFymClTcKFjcXMhzPwSp5vSJHkEKOAH5BREREt2ARxEOe3OICAKN7RuOB1Dgs3nUaRRUmj2Z38BYXIiIiul3bj13F1ycK2vxzlYEBCFbK6woW9meN0lG0kENTV5DQKOUIVgZCU1eo0Cjtxx3Fi+Cb9gXydhEiImplLIJ4qLDCfQEEAB7qp8fE/vHQqOTNmt3BW1yIiIjodgQpm75VQyGXIUhhLzYEKeXQKOXQKOoXLoKUcmdRQ1NXxAhWyuv2uxYsNEo5+1sQEZFPYhHEQ47GpJ6edzuzO3iLCxERETXXz3+Sgof66aFR1BUwVDcKHUFKOZSBLFYQERE5sAjiobtSIhCnU7OBKREREUlKarzO20MgIiLyGSyCeIgNTImIiIiIiIh8G+dHNoPjFpdYneutMbE6NVY9M5ANTImIiIiIiIgkjDNBmom3uBARERERERH5JhZBbgNvcSEiIiIiIiLyPbwdhoiIiIiIiIjaBRZBiIiIiIiIiKhdYBGEiIiIiIiIiNoFFkGIiIiIiIiIqF1gEYSIiIiIiIiI2gUWQYiIiIiIiIioXWARhIiIiKiVpaeno1evXhgyZIi3h0JERNSuyYQQwtuDaE0GgwFhYWG4dOkStFqtt4dDRETkN8rLy5GYmIiysjLodDpvD8cnMC8hIiJqHZ7mJYFtOCavqKioAAAkJiZ6eSRERET+qaKigkUQDzEvISIial3u8hK/nwlis9lw9epVhIaGQiaTtdh1HVUmf/omhzH5Dn+MizH5Dn+MizHdHiEEKioqoNfrERDAO2w9wbzEc4zJd/hjXIzJd/hjXP4YE9D6cXmal/j9TJCAgAAkJCS02vW1Wq1f/YcJMCZf4o9xMSbf4Y9xMabm4wyQ5mFe0nyMyXf4Y1yMyXf4Y1z+GBPQunF5kpfwaxsiIiIiIiIiahdYBCEiIiIiIiKidoFFkNukUqkwf/58qFQqbw+lxTAm3+GPcTEm3+GPcTEm8nX++PNmTL7DH+NiTL7DH+Pyx5gA6cTl941RiYiIiIiIiIgAzgQhIiIiIiIionaCRRAiIiIiIiIiahdYBCEiIiIiIiKidoFFECIiIiIiIiJqF1gEqbN48WIMGTIEoaGhiI6OxqRJk5CTk+M8brFY8Prrr6NPnz4IDg6GXq/Hc889h6tXr7pcx2Qy4ZVXXkFkZCSCg4Px8MMP4/Lly20djpO7uABgwYIF6NGjB4KDgxEeHo7Ro0fj8OHDLudIKS5PYrrZiy++CJlMhuXLl7vs97WYpk2bBplM5vIYNmyYyzlSignw/Gd16tQpPPzww9DpdAgNDcWwYcOQl5fnPC6luDyJ6dafk+Px3nvvOc/xtZgqKysxa9YsJCQkICgoCD179sSqVatczpFSTIBncV27dg3Tpk2DXq+HRqPB+PHj8eOPP7qcI6W4Vq1ahb59+0Kr1UKr1SItLQ1ff/2187gQAgsWLIBer0dQUBBGjBiBH374weUaUoqHGse8hHmJ1GNiXiKNuPwxLwH8MzdhXiKhvESQEEKIcePGiXXr1okTJ06I7Oxs8dOf/lQkJSWJyspKIYQQZWVlYvTo0WLz5s3i9OnT4tChQ2Lo0KFi0KBBLtd56aWXRHx8vMjIyBBZWVli5MiRol+/fqK2ttYbYbmNSwgh/va3v4mMjAxx7tw5ceLECTF9+nSh1WpFYWGh8xwpxeVJTA5bt24V/fr1E3q9XnzwwQcux3wtpqlTp4rx48eL/Px856O4uNjlOlKKSQjP4jp79qyIiIgQv/nNb0RWVpY4d+6c2LFjh7h27ZrzHCnF5UlMN/+M8vPzxSeffCJkMpk4d+6cz8b0/PPPi86dO4s9e/aI3NxcsXr1aiGXy8W2bdskGZMQ7uOy2Wxi2LBh4t577xVHjhwRp0+fFjNmzKgXu5Ti2r59u/jqq69ETk6OyMnJEfPmzRMKhUKcOHFCCCHEkiVLRGhoqNiyZYs4fvy4mDJlioiLixPl5eWSjIcax7yEeYnUY2JeIo24/DEvEcI/cxPmJdLJS1gEaURhYaEAIPbt29foOUeOHBEAxMWLF4UQ9oREoVCITZs2Oc+5cuWKCAgIELt27Wr1MXvCk7gMBoMAIP75z38KIaQfV2MxXb58WcTHx4sTJ06I5ORkl2TDF2OaOnWqmDhxYqPvkXpMQjQc15QpU8QzzzzT6HukHpcn/09NnDhRjBo1yrntizH17t1b/O53v3M5b+DAgeLNN98UQkg/JiHqx5WTkyMAOH9RCyFEbW2tiIiIEGvXrhVC+EZc4eHh4qOPPhI2m03ExsaKJUuWOI8ZjUah0+nEn//8ZyGEb8RDDWNewrzEm5iX3CD1uPwxLxHCP3MT5iXei4e3wzTCYDAAACIiIpo8RyaTISwsDACQmZkJi8WCsWPHOs/R6/VITU3FwYMHW3W8nnIXl9lsxpo1a6DT6dCvXz8A0o+roZhsNhueffZZ/OY3v0Hv3r3rvccXYwKAvXv3Ijo6Gt26dcMLL7yAwsJC5zGpxwTUj8tms+Grr75Ct27dMG7cOERHR2Po0KHYtm2b8z1Sj8vd/1PXrl3DV199henTpzv3+WJM99xzD7Zv344rV65ACIE9e/bgzJkzGDduHADpxwTUj8tkMgEA1Gq18xy5XA6lUon//Oc/AKQdl9VqxaZNm1BVVYW0tDTk5uaioKDAZawqlQrDhw93jlXK8VDTmJcwL/Em5iXbnO+Relz+mJcA/pmbMC/xXjwsgjRACIE5c+bgnnvuQWpqaoPnGI1GvPHGG3jqqaeg1WoBAAUFBVAqlQgPD3c5NyYmBgUFBa0+bneaimvHjh0ICQmBWq3GBx98gIyMDERGRgKQdlyNxbR06VIEBgbi1VdfbfB9vhjThAkT8Le//Q3//ve/8f777+Po0aMYNWqU8y9MKccENBxXYWEhKisrsWTJEowfPx7ffPMNHnnkEUyePBn79u0DIO24PPm7YsOGDQgNDcXkyZOd+3wxpg8//BC9evVCQkIClEolxo8fj5UrV+Kee+4BIO2YgIbj6tGjB5KTkzF37lyUlpbCbDZjyZIlKCgoQH5+PgBpxnX8+HGEhIRApVLhpZdewtatW9GrVy/neGJiYlzOv3msUoyH3GNewrzEm5iXMC/xNn/MTZiXeDcvCWy1K/uwWbNm4fvvv3dW3G5lsVjwxBNPwGazYeXKlW6vJ4SATCZr6WE2W1NxjRw5EtnZ2SgqKsLatWvx+OOP4/Dhw4iOjm70elKIq6GYMjMzsWLFCmRlZTV7fFKNCQCmTJnifJ2amorBgwcjOTkZX331lcsvsltJISag4bhsNhsAYOLEifjVr34FAOjfvz8OHjyIP//5zxg+fHij15NCXO7+rgCATz75BE8//bRLVb8xUo7pww8/xLfffovt27cjOTkZ+/fvxy9+8QvExcVh9OjRjV5PCjEBDcelUCiwZcsWTJ8+HREREZDL5Rg9ejQmTJjg9nrejKt79+7Izs5GWVkZtmzZgqlTpzqTcwD1xuXJWKXyc6KGMS9hXuJNzEuYl3ibP+YmzEu8m5dwJsgtXnnlFWzfvh179uxBQkJCveMWiwWPP/44cnNzkZGR4fy2BQBiY2NhNptRWlrq8p7CwsJ6FbC25i6u4OBgdOnSBcOGDcPHH3+MwMBAfPzxxwCkG1djMR04cACFhYVISkpCYGAgAgMDcfHiRbz22mvo2LEjAN+LqSFxcXFITk52doyWakxA43FFRkYiMDAQvXr1cjm/Z8+ezi7sUo3Lk5/VgQMHkJOTg+eff95lv6/FVFNTg3nz5mHZsmV46KGH0LdvX8yaNQtTpkzBH//4RwDSjQlo+mc1aNAg5y/u/Px87Nq1C8XFxUhJSQEgzbiUSiW6dOmCwYMHY/HixejXrx9WrFiB2NhYAKj3zcnNY5ViPNQ05iXMS6QYU0OYl0j/Z+VreQngn7kJ8xIJ5CWt1m3Ex9hsNjFz5kyh1+vFmTNnGjzHbDaLSZMmid69e7t0KHdwNHbZvHmzc9/Vq1e92qjGk7ga0rlzZzF//nwhhPTichdTUVGROH78uMtDr9eL119/XZw+fVoI4XsxNaSoqEioVCqxYcMGIYT0YhLCs7jS0tLqNSCbNGmSePLJJ4UQ0ourOT+rqVOn1lupQQjfi8nRlHDnzp0u+2fMmCHGjBkjhJBeTELc3v9XZ86cEQEBAWL37t1CCGnGdatRo0aJqVOnOhuQLV261HnMZDI12IBMyvGQHfMSV8xL2hbzEuYlUv9Z+WJuwrxEOnkJiyB1Xn75ZaHT6cTevXtdlpCqrq4WQghhsVjEww8/LBISEkR2drbLOSaTyXmdl156SSQkJIh//vOfIisrS4waNcqry0u5i6uyslLMnTtXHDp0SFy4cEFkZmaK6dOnC5VK5dKZWEpxuYupIbd2YRfCt2KqqKgQr732mjh48KDIzc0Ve/bsEWlpaSI+Pr7eElNSicmTuIQQ4osvvhAKhUKsWbNG/Pjjj+J///d/hVwuFwcOHHCeI6W4PP3vz2AwCI1GI1atWtXgdXwtpuHDh4vevXuLPXv2iPPnz4t169YJtVotVq5cKcmYhPAsrs8//1zs2bNHnDt3Tmzbtk0kJyeLyZMnu1xHSnHNnTtX7N+/X+Tm5orvv/9ezJs3TwQEBIhvvvlGCGFfik6n04kvvvhCHD9+XDz55JMNLkUnlXioccxLmJdIOSbmJdKJyx/zEiH8MzdhXiKdvIRFkDoAGnysW7dOCCFEbm5uo+fs2bPHeZ2amhoxa9YsERERIYKCgsSDDz4o8vLyvBOUcB9XTU2NeOSRR4RerxdKpVLExcWJhx9+WBw5csTlOlKKy11MDWko2fClmKqrq8XYsWNFVFSUUCgUIikpSUydOrXeeKUUkxCe/6w+/vhj0aVLF6FWq0W/fv1c1ncXQlpxeRrT6tWrRVBQkCgrK2vwOr4WU35+vpg2bZrQ6/VCrVaL7t27i/fff1/YbDbnOVKKSQjP4lqxYoVISEhw/n/15ptvuvwDUghpxfXzn/9cJCcnC6VSKaKiosT999/vTDSEsH/LNH/+fBEbGytUKpW47777xPHjx12uIaV4qHHMS5iXSDkm5iXSicsf8xIh/DM3YV4inbxEJoQQDd0mQ0RERERERETkT9gYlYiIiIiIiIjaBRZBiIiIiIiIiKhdYBGEiIiIiIiIiNoFFkGIiIiIiIiIqF1gEYSIiIiIiIiI2gUWQYiIiIiIiIioXWARhIiIiIiIiIjaBRZBiIiIiIiIiKhdYBGEiIiIiIiIiNoFFkGIfJhMJmvyMW3aNG8PscWNGDECs2fP9vYwXEybNg2TJk1y2ff3v/8darUa7777rncGRURE5AXMTaSBuQlR4wK9PQAiun35+fnO15s3b8bbb7+NnJwc576goCBvDOu2WCwWKBQKv/i8jz76CDNnzkR6ejqef/75VvkMIiIiKWJuIs3PY25CdANnghD5sNjYWOdDp9NBJpO57Nu/fz8GDRoEtVqNTp06YeHChaitrXW+XyaTYfXq1XjwwQeh0WjQs2dPHDp0CGfPnsWIESMQHByMtLQ0nDt3zvmeBQsWoH///li9ejUSExOh0Wjw2GOPoayszGVs69atQ8+ePaFWq9GjRw+sXLnSeezChQuQyWT4/PPPMWLECKjVanz66acoLi7Gk08+iYSEBGg0GvTp0wefffaZ833Tpk3Dvn37sGLFCuc3ShcuXMD69esRFhbm8vnbtm2DTCarN+5PPvkEnTp1gkqlghACBoMBM2bMQHR0NLRaLUaNGoVjx47d9s/k3XffxaxZs7Bx40YmGURE1O4wN2FuQiR1LIIQ+andu3fjmWeewauvvoqTJ09i9erVWL9+Pf7whz+4nPfOO+/gueeeQ3Z2Nnr06IGnnnoKL774IubOnYvvvvsOADBr1iyX95w9exaff/45vvzyS+zatQvZ2dmYOXOm8/jatWvx29/+Fn/4wx9w6tQpLFq0CG+99RY2bNjgcp3XX38dr776Kk6dOoVx48bBaDRi0KBB2LFjB06cOIEZM2bg2WefxeHDhwEAK1asQFpaGl544QXk5+cjPz8fiYmJHv+ZOMa9ZcsWZGdnAwB++tOfoqCgADt37kRmZiYGDhyI+++/HyUlJR5f1+GNN97AO++8gx07duDRRx9t9vuJiIj8GXOT+pibEHmBICK/sG7dOqHT6Zzb9957r1i0aJHLOX/9619FXFyccxuAePPNN53bhw4dEgDExx9/7Nz32WefCbVa7dyeP3++kMvl4tKlS859X3/9tQgICBD5+flCCCESExPFxo0bXT77nXfeEWlpaUIIIXJzcwUAsXz5crdxPfDAA+K1115zbg8fPlz88pe/bDJ2IYTYunWruPmvuPnz5wuFQiEKCwud+/71r38JrVYrjEajy3s7d+4sVq9e7XZsDlOnThVKpVIAEP/61788fh8REZE/Y26ic9nH3IRIGtgThMhPZWZm4ujRoy7frlitVhiNRlRXV0Oj0QAA+vbt6zweExMDAOjTp4/LPqPRiPLycmi1WgBAUlISEhISnOekpaXBZrMhJycHcrkcly5dwvTp0/HCCy84z6mtrYVOp3MZ4+DBg122rVYrlixZgs2bN+PKlSswmUwwmUwIDg6+0z8OAEBycjKioqKc25mZmaisrESHDh1czqupqXGZZuuJvn37oqioCG+//TaGDBmC0NDQFhkzERGRv2BuUh9zE6K2xyIIkZ+y2WxYuHAhJk+eXO+YWq12vr65AZfjPtWG9tlstkY/y3GOTCZznrd27VoMHTrU5Ty5XO6yfWsC8f777+ODDz7A8uXL0adPHwQHB2P27Nkwm82NBwogICAAQgiXfRaLpd55t36ezWZDXFwc9u7dW+/cW+/jdSc+Ph5btmzByJEjMX78eOzatYvJBhER0U2YmzA3IZICFkGI/NTAgQORk5ODLl26tPi18/LycPXqVej1egDAoUOHEBAQgG7duiEmJgbx8fE4f/48nn766WZd98CBA5g4cSKeeeYZAPZE4Mcff0TPnj2d5yiVSlitVpf3RUVFoaKiAlVVVc5kwnFfbVMGDhyIgoICBAYGomPHjs0aa0OSkpKwb98+jBw5EmPHjsXu3bud31ARERG1d8xNst1+HnMTotbHxqhEfurtt9/GX/7yFyxYsAA//PADTp06hc2bN+PNN9+842ur1WpMnToVx44dw4EDB/Dqq6/i8ccfR2xsLAB7t/PFixdjxYoVOHPmDI4fP45169Zh2bJlTV63S5cuyMjIwMGDB3Hq1Cm8+OKLKCgocDmnY8eOOHz4MC5cuICioiLYbDYMHToUGo0G8+bNw9mzZ7Fx40asX7/ebRyjR49GWloaJk2ahN27d+PChQs4ePAg3nzzTWfjteZKSEjA3r17UVxcjLFjx8JgMNzWdYiIiPwNc5P1buNgbkLU+lgEIfJT48aNw44dO5CRkYEhQ4Zg2LBhWLZsGZKTk+/42l26dMHkyZPxwAMPYOzYsUhNTXVZZu7555/HRx99hPXr16NPnz4YPnw41q9fj5SUlCav+9Zbb2HgwIEYN24cRowYgdjYWEyaNMnlnF//+teQy+Xo1asXoqKikJeXh4iICHz66afYuXOnc+m6BQsWuI1DJpNh586duO+++/Dzn/8c3bp1wxNPPIELFy4470F2LJnX0LTUxsTHx2Pfvn0oKyvDmDFj6i3RR0RE1B4xN1ngNg7mJkStTyZuvVmNiKgJCxYswLZt2zya0ukP9u7di0ceeQTnz59HeHi4t4dDREREt2BuQkTNwZkgRERN2LVrF+bNm8ckg4iIiCSBuQnRnWFjVCKiJixZssTbQyAiIiJyYm5CdGd4OwwRERERERERtQu8HYaIiIiIiIiI2gUWQYiIiIiIiIioXWARhIiIiIiIiIjaBRZBiIiIiIiIiKhdYBGEiIiIiIiIiNoFFkGIiIiIiIiIqF1gEYSIiIiIiIiI2gUWQYiIiIiIiIioXfj/g5NG9jfyNwcAAAAASUVORK5CYII=",
|
||
"text/plain": [
|
||
"<Figure size 1100x420 with 2 Axes>"
|
||
]
|
||
},
|
||
"metadata": {},
|
||
"output_type": "display_data"
|
||
}
|
||
],
|
||
"source": [
|
||
"#!/usr/bin/env python3\n",
|
||
"# -*- coding: utf-8 -*-\n",
|
||
"\"\"\"\n",
|
||
"Pure-component saturation with Peng–Robinson EOS (phi^L = phi^V)\n",
|
||
"Outputs Psat(T) *and* saturated liquid/vapor Z and molar volumes.\n",
|
||
"Units: bar, L/mol, K\n",
|
||
"\"\"\"\n",
|
||
"import math\n",
|
||
"from typing import List, Optional, Tuple\n",
|
||
"\n",
|
||
"try:\n",
|
||
" import numpy as np\n",
|
||
"except Exception:\n",
|
||
" class _NP:\n",
|
||
" @staticmethod\n",
|
||
" def linspace(a, b, n):\n",
|
||
" step = (b - a) / (n - 1)\n",
|
||
" return [a + i * step for i in range(n)]\n",
|
||
" np = _NP()\n",
|
||
"\n",
|
||
"R = 0.08314472 # bar·L/(mol·K)\n",
|
||
"\n",
|
||
"def pr_params(Tc: float, Pc: float, omega: float) -> Tuple[float, float, float]:\n",
|
||
" a = 0.45724 * R**2 * Tc**2 / Pc\n",
|
||
" b = 0.07780 * R * Tc / Pc\n",
|
||
" kappa = 0.37464 + 1.54226 * omega - 0.26992 * omega**2\n",
|
||
" return a, b, kappa\n",
|
||
"\n",
|
||
"def alpha_pr(T: float, Tc: float, kappa: float) -> float:\n",
|
||
" Tr = T / Tc\n",
|
||
" return (1 + kappa * (1 - math.sqrt(Tr))) ** 2\n",
|
||
"\n",
|
||
"def A_B(a: float, b: float, alpha: float, T: float, P: float) -> Tuple[float, float]:\n",
|
||
" A = a * alpha * P / (R**2 * T**2)\n",
|
||
" B = b * P / (R * T)\n",
|
||
" return A, B\n",
|
||
"\n",
|
||
"def _real_cubic_roots(a3: float, a2: float, a1: float, a0: float) -> List[float]:\n",
|
||
" if abs(a3) < 1e-18:\n",
|
||
" if abs(a2) < 1e-18:\n",
|
||
" if abs(a1) < 1e-18:\n",
|
||
" return []\n",
|
||
" return [-a0 / a1]\n",
|
||
" disc = a1 * a1 - 4 * a2 * a0\n",
|
||
" if disc < 0:\n",
|
||
" return []\n",
|
||
" if abs(disc) < 1e-18:\n",
|
||
" return [-a1 / (2 * a2)]\n",
|
||
" rt = math.sqrt(disc)\n",
|
||
" return [(-a1 - rt) / (2 * a2), (-a1 + rt) / (2 * a2)]\n",
|
||
" p = a2 / a3\n",
|
||
" q = a1 / a3\n",
|
||
" r = a0 / a3\n",
|
||
" shift = p / 3.0\n",
|
||
" Acoef = q - p * p / 3.0\n",
|
||
" Bcoef = (2 * p**3) / 27.0 - p * q / 3.0 + r\n",
|
||
" disc = (Bcoef / 2.0) ** 2 + (Acoef / 3.0) ** 3\n",
|
||
" roots = []\n",
|
||
" if disc > 0:\n",
|
||
" sqrt_disc = math.sqrt(disc)\n",
|
||
" u = -Bcoef / 2.0 + sqrt_disc\n",
|
||
" v = -Bcoef / 2.0 - sqrt_disc\n",
|
||
" u_c = math.copysign(abs(u) ** (1 / 3.0), u) if u != 0 else 0.0\n",
|
||
" v_c = math.copysign(abs(v) ** (1 / 3.0), v) if v != 0 else 0.0\n",
|
||
" x = u_c + v_c\n",
|
||
" roots = [x - shift]\n",
|
||
" elif abs(disc) < 1e-18:\n",
|
||
" if abs(Bcoef) < 1e-18:\n",
|
||
" x = 0.0\n",
|
||
" roots = [x - shift]\n",
|
||
" else:\n",
|
||
" u = math.copysign(abs(-Bcoef / 2.0) ** (1 / 3.0), -Bcoef / 2.0)\n",
|
||
" roots = [2 * u - shift, -u - shift]\n",
|
||
" else:\n",
|
||
" rho = math.sqrt(-Acoef**3 / 27.0)\n",
|
||
" theta = math.acos(-Bcoef / (2.0 * rho))\n",
|
||
" m = 2 * math.sqrt(-Acoef / 3.0)\n",
|
||
" x1 = m * math.cos(theta / 3.0) - shift\n",
|
||
" x2 = m * math.cos((theta + 2 * math.pi) / 3.0) - shift\n",
|
||
" x3 = m * math.cos((theta + 4 * math.pi) / 3.0) - shift\n",
|
||
" roots = [x1, x2, x3]\n",
|
||
" roots = [float(r) for r in roots]\n",
|
||
" roots.sort()\n",
|
||
" return roots\n",
|
||
"\n",
|
||
"def pr_Z_roots(A: float, B: float) -> List[float]:\n",
|
||
" a3 = 1.0\n",
|
||
" a2 = -(1.0 - B)\n",
|
||
" a1 = A - 3 * B * B - 2 * B\n",
|
||
" a0 = -(A * B - B * B - B**3)\n",
|
||
" roots = _real_cubic_roots(a3, a2, a1, a0)\n",
|
||
" roots = [z for z in roots if z > -1e-10]\n",
|
||
" return roots\n",
|
||
"\n",
|
||
"def lnphi_pure(Z: float, A: float, B: float) -> float:\n",
|
||
" sqrt2 = math.sqrt(2.0)\n",
|
||
" term1 = Z - 1.0 - math.log(max(Z - B, 1e-300))\n",
|
||
" term2 = (A / (2.0 * sqrt2 * B)) * math.log(\n",
|
||
" max((Z + (1.0 + sqrt2) * B) / (Z + (1.0 - sqrt2) * B), 1e-300)\n",
|
||
" )\n",
|
||
" return term1 - term2\n",
|
||
"\n",
|
||
"def lnphi_L_minus_V_and_roots(T: float, P: float, Tc: float, Pc: float, omega: float):\n",
|
||
" a, b, kappa = pr_params(Tc, Pc, omega)\n",
|
||
" alpha = alpha_pr(T, Tc, kappa)\n",
|
||
" A, B = A_B(a, b, alpha, T, P)\n",
|
||
" roots = pr_Z_roots(A, B)\n",
|
||
" if len(roots) < 2:\n",
|
||
" return None, None, None, None, None\n",
|
||
" Z_L = min(roots)\n",
|
||
" Z_V = max(roots)\n",
|
||
" lnphi_L = lnphi_pure(Z_L, A, B)\n",
|
||
" lnphi_V = lnphi_pure(Z_V, A, B)\n",
|
||
" return lnphi_L - lnphi_V, Z_L, Z_V, A, B\n",
|
||
"\n",
|
||
"def find_psat_at_T(T: float, Tc: float, Pc: float, omega: float,\n",
|
||
" Pmin: float = 1e-6, Pmax: Optional[float] = None, nscan: int = 200):\n",
|
||
" if Pmax is None:\n",
|
||
" Pmax = 0.999 * Pc\n",
|
||
" Ps = [Pmin * (Pmax / Pmin) ** (i / (nscan - 1)) for i in range(nscan)]\n",
|
||
" vals = []\n",
|
||
" for P in Ps:\n",
|
||
" F, ZL, ZV, A, B = lnphi_L_minus_V_and_roots(T, P, Tc, Pc, omega)\n",
|
||
" if F is not None and math.isfinite(F):\n",
|
||
" vals.append((P, F))\n",
|
||
" if len(vals) < 2:\n",
|
||
" return None\n",
|
||
" bracket = None\n",
|
||
" for (P1, F1), (P2, F2) in zip(vals[:-1], vals[1:]):\n",
|
||
" if F1 == 0.0:\n",
|
||
" P = P1\n",
|
||
" _, ZL, ZV, _, _ = lnphi_L_minus_V_and_roots(T, P, Tc, Pc, omega)\n",
|
||
" VL = ZL * R * T / P\n",
|
||
" VV = ZV * R * T / P\n",
|
||
" return P, ZL, ZV, VL, VV\n",
|
||
" if F1 * F2 < 0.0:\n",
|
||
" bracket = (P1, F1, P2, F2)\n",
|
||
" break\n",
|
||
" if bracket is None:\n",
|
||
" vals_sorted = sorted(vals, key=lambda t: abs(t[1]))\n",
|
||
" if len(vals_sorted) >= 2:\n",
|
||
" P1, F1 = vals_sorted[0]\n",
|
||
" P2, F2 = vals_sorted[1]\n",
|
||
" for _ in range(50):\n",
|
||
" if abs(F2 - F1) < 1e-14:\n",
|
||
" break\n",
|
||
" P_new = P2 - F2 * (P2 - P1) / (F2 - F1)\n",
|
||
" P_new = max(Pmin, min(P_new, Pmax))\n",
|
||
" F_new, ZL, ZV, _, _ = lnphi_L_minus_V_and_roots(T, P_new, Tc, Pc, omega)\n",
|
||
" if F_new is None:\n",
|
||
" P_new = 0.5 * (P2 + P1)\n",
|
||
" F_new, ZL, ZV, _, _ = lnphi_L_minus_V_and_roots(T, P_new, Tc, Pc, omega)\n",
|
||
" if F_new is None:\n",
|
||
" break\n",
|
||
" if abs(F_new) < 1e-10:\n",
|
||
" VL = ZL * R * T / P_new\n",
|
||
" VV = ZV * R * T / P_new\n",
|
||
" return P_new, ZL, ZV, VL, VV\n",
|
||
" P1, F1, P2, F2 = P2, F2, P_new, F_new\n",
|
||
" return None\n",
|
||
" P1, F1, P2, F2 = bracket\n",
|
||
" ZL=ZV=VL=VV=None\n",
|
||
" for _ in range(80):\n",
|
||
" Pm = 0.5 * (P1 + P2)\n",
|
||
" Fm, ZL, ZV, _, _ = lnphi_L_minus_V_and_roots(T, Pm, Tc, Pc, omega)\n",
|
||
" if Fm is None:\n",
|
||
" Pm = math.sqrt(P1 * P2)\n",
|
||
" Fm, ZL, ZV, _, _ = lnphi_L_minus_V_and_roots(T, Pm, Tc, Pc, omega)\n",
|
||
" if Fm is None:\n",
|
||
" break\n",
|
||
" if abs(Fm) < 1e-12:\n",
|
||
" VL = ZL * R * T / Pm\n",
|
||
" VV = ZV * R * T / Pm\n",
|
||
" return Pm, ZL, ZV, VL, VV\n",
|
||
" if F1 * Fm < 0:\n",
|
||
" P2, F2 = Pm, Fm\n",
|
||
" else:\n",
|
||
" P1, F1 = Pm, Fm\n",
|
||
" if ZL is None or ZV is None:\n",
|
||
" return None\n",
|
||
" Pfin = 0.5 * (P1 + P2)\n",
|
||
" VL = ZL * R * T / Pfin\n",
|
||
" VV = ZV * R * T / Pfin\n",
|
||
" return Pfin, ZL, ZV, VL, VV\n",
|
||
"\n",
|
||
"def compute_psat_curve_pr(Tc: float, Pc: float, omega: float,\n",
|
||
" T_min: Optional[float] = None,\n",
|
||
" T_max: Optional[float] = None,\n",
|
||
" nT: int = 60):\n",
|
||
" if T_max is None:\n",
|
||
" T_max = 0.999 * Tc\n",
|
||
" if T_min is None:\n",
|
||
" T_min = 0.55 * Tc\n",
|
||
" Ts = np.linspace(T_min, T_max, nT)\n",
|
||
" rows = []\n",
|
||
" for T in Ts:\n",
|
||
" out = find_psat_at_T(T, Tc, Pc, omega)\n",
|
||
" if out is None:\n",
|
||
" rows.append((float(T), None, None, None, None, None))\n",
|
||
" else:\n",
|
||
" P, ZL, ZV, VL, VV = out\n",
|
||
" rows.append((float(T), float(P), float(ZL), float(ZV), float(VL), float(VV)))\n",
|
||
" return rows\n",
|
||
"\n",
|
||
"def _maybe_plot(results):\n",
|
||
" try:\n",
|
||
" import matplotlib.pyplot as plt\n",
|
||
" except Exception:\n",
|
||
" return\n",
|
||
" # T-P plot (single line)\n",
|
||
" T = [r[0] for r in results if r[1] is not None]\n",
|
||
" P = [r[1] for r in results if r[1] is not None]\n",
|
||
" if len(T) == 0:\n",
|
||
" return\n",
|
||
" fig, ax = plt.subplots(1, 2, figsize=(11,4.2))\n",
|
||
" ax[0].plot(T, P, marker=\"o\", lw=1.5)\n",
|
||
" ax[0].set_xlabel(\"Temperature, K\")\n",
|
||
" ax[0].set_ylabel(\"Saturation Pressure, bar\")\n",
|
||
" ax[0].set_title(\"PR Saturation Curve: T–P\")\n",
|
||
"\n",
|
||
" # T–V plot (two branches)\n",
|
||
" VL = [r[4] for r in results if r[1] is not None] # L molar volume\n",
|
||
" VV = [r[5] for r in results if r[1] is not None] # L molar volume\n",
|
||
" ax[1].plot(T, VL, label=\"V_L (sat. liquid)\", lw=2)\n",
|
||
" ax[1].plot(T, VV, label=\"V_V (sat. vapor)\", lw=2)\n",
|
||
" ax[1].set_xlabel(\"Temperature, K\")\n",
|
||
" ax[1].set_ylabel(\"Molar volume, L/mol\")\n",
|
||
" ax[1].set_yscale(\"log\")\n",
|
||
" ax[1].set_title(\"Saturated Volumes\")\n",
|
||
" ax[1].legend()\n",
|
||
" plt.tight_layout()\n",
|
||
" plt.show()\n",
|
||
"\n",
|
||
"if __name__ == \"__main__\":\n",
|
||
" # Example: CO2\n",
|
||
" Tc = 304.2 # K\n",
|
||
" Pc = 73.8 # bar\n",
|
||
" omega = 0.225\n",
|
||
" print(\"Computing PR saturation curve for CO2...\")\n",
|
||
" table = compute_psat_curve_pr(Tc, Pc, omega, T_min=220.0, T_max=0.995*Tc, nT=40)\n",
|
||
" for row in table[:5]:\n",
|
||
" T, P, ZL, ZV, VL, VV = row\n",
|
||
" print(f\"T={T:8.3f} K Psat={P if P is not None else float('nan'):9.5f} bar \"\n",
|
||
" f\"ZL={ZL if ZL is not None else float('nan'):7.4f} \"\n",
|
||
" f\"ZV={ZV if ZV is not None else float('nan'):7.4f}\")\n",
|
||
" _maybe_plot(table)\n"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "d9f12983-dcfa-4795-91c7-11eb45aa104c",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": []
|
||
}
|
||
],
|
||
"metadata": {
|
||
"kernelspec": {
|
||
"display_name": "Python [conda env:base] *",
|
||
"language": "python",
|
||
"name": "conda-base-py"
|
||
},
|
||
"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.9"
|
||
}
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 5
|
||
}
|