{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# General Usage" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Import pygcc" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.5.2\n" ] } ], "source": [ "from importlib.metadata import version\n", "print(version(\"pygcc\"))\n", "\n", "from pygcc.pygcc_utils import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Read database by specifying the direct- or sequential- access and source database" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Using the default sequential-access database - speq21.dat" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "ps = db_reader(sourcedb = 'thermo.2021', sourceformat = 'gwb')\n", "# ps.dbaccessdic, ps.sourcedic, ps.specielist" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Using the user-specified sequential-access database" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Duplicate found for species \"AlOH++\" in slop07.dat\n", "Duplicate found for species \"MgCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"CaCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"SrCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"BaCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"BaF+\" in slop07.dat\n", "Duplicate found for species \"Sr(Succ)(aq)\" in slop07.dat\n", "Duplicate found for species \"Sc(Glut)+(aq)\" in slop07.dat\n", "Duplicate found for species \"AMP2-\" in slop07.dat\n", "Duplicate found for species \"HAMP-\" in slop07.dat\n", "Duplicate found for species \"+H2AMP-\" in slop07.dat\n" ] } ], "source": [ "ps = db_reader(dbaccess = './database/slop07.dat', sourcedb = 'thermo.2021', sourceformat = 'gwb')\n", "# ps.dbaccessdic, ps.sourcedic, ps.specielist" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate water properties" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*With IAPWS95*" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Density: [1009.73580931 1005.83998999 991.7058825 967.43838406 927.69054221\n", " 877.9652082 816.08878805 734.71208475]\n", "Gibbs Energy: [-56194.26867312 -56592.33889817 -57211.87993461 -58000.04190115\n", " -59093.00868798 -60294.40876454 -61596.26593847 -62995.11966287]\n", "Enthalpy: [-68680.33499723 -68236.29564253 -67613.18053963 -66897.34619192\n", " -65991.9279271 -65062.66635227 -64087.83220987 -63021.29280159]\n", "Entropy: [15.14005087 16.69550855 18.67153539 20.7005717 22.97720284 25.05223509\n", " 27.00976077 28.9549847 ]\n" ] } ], "source": [ "water = iapws95(T = np.array([ 0.01, 25, 60, 100, 150, 200, 250, 300]), P = 200)\n", "print('Density:', water.rho)\n", "print('Gibbs Energy:', water.G)\n", "print('Enthalpy:', water.H)\n", "print('Entropy:', water.S)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*With ZhangDuan*" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[175.8286937 314.74510733] [-88767.03099 -89082.94223142]\n" ] } ], "source": [ "water = ZhangDuan(T= np.array([1000, 1050]), P = np.array([1000, 2000]))\n", "print(water.rho, water.G)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate water dielectric constants" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([1.19637079, 2.52660344]),\n", " array([0.17582869, 0.31474511]),\n", " array([12.8720893 , 5.29642074]),\n", " array([0.5403405 , 0.48797969]))" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dielect = water_dielec(T= np.array([1000, 1050]), P = np.array([1000, 2000]), Dielec_method = 'DEW')\n", "dielect.E, dielect.rhohat, dielect.Ah, dielect.Bh" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([55.83017195, 44.79388023]),\n", " array([0.96293375, 0.92769054]),\n", " array([0.59551378, 0.67352582]),\n", " array([0.34191435, 0.35183609]))" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dielect = water_dielec(T= np.array([100, 150]), P = np.array([100, 200]), Dielec_method = 'JN91')\n", "dielect.E, dielect.rhohat, dielect.Ah, dielect.Bh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate quartz properties" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*With Maier-Kelly powerlaw*" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([-205319.43164294, -206728.65136032, -210399.56629393,\n", " -215044.47020067, -220518.04896025, -223492.57970966]),\n", " array([12.34074489, 13.89377788, 16.14397572, 16.103911 , 16.491911 ,\n", " 16.685911 ]))" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Temp = np.array([100, 200, 400, 600, 800, 900])\n", "ps = db_reader()\n", "sup = heatcap( T = Temp, P = 1, Species = 'Quartz', Species_ppt = ps.dbaccessdic['Quartz'], \n", " method = 'SUPCRT')\n", "sup.dG, sup.dCp" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*With Holland and Power's formulation*" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Duplicate found for species \"Fluorphlogopite\" in supcrtbl.dat\n", "Duplicate found for species \"Arsenic\" in supcrtbl.dat\n", "Duplicate found for species \"Arsenolite\" in supcrtbl.dat\n" ] }, { "data": { "text/plain": [ "(array([-205491.44716714, -206900.38310661, -210575.79120014,\n", " -215228.32290096, -220682.44526921, -223649.25812766]),\n", " array([12.4075447 , 13.99685509, 16.16426364, 16.05341748, 16.6660167 ,\n", " 16.90252028]))" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Temp = np.array([100, 200, 400, 600, 800, 900])\n", "pshp = db_reader(dbHP_dir = './database/supcrtbl.dat')\n", "hp = heatcap( T = Temp, P = 10, Species = 'Quartz', Species_ppt = pshp.dbaccessdic['Quartz'], \n", " method = 'HP11')\n", "hp.dG, hp.dCp" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*With Berman's formulation*" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Duplicate found for species \"Cordierite\" in berman.dat\n" ] }, { "data": { "text/plain": [ "(array([-205498.95550508, -206907.11144614, -210575.29239247,\n", " -215221.44900321, -220654.77893829, -223612.45803826]),\n", " array([12.323813 , 13.87414095, 16.29997608, 16.24449507, 16.72930639,\n", " 16.90352372]))" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Temp = np.array([100, 200, 400, 600, 800, 900])\n", "ps = db_reader(dbBerman_dir = './database/berman.dat')\n", "bm = heatcap( T = Temp, P = 1, Species = 'Quartz', Species_ppt = ps.dbaccessdic['Quartz'], \n", " method = 'berman88')\n", "bm.dG, bm.dCp" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAobxJREFUeJzt3QV4U2cXB/B/S43iXtzd3Z3hzrcNdzaGuzvDbcBwZ8CAbQwfDHd3d3d3Smnv95z3kiwtBdrS9Eb+v+e5NE3S5E1o3p77yjkumqZpICIiIiK752p0A4iIiIgofDCwIyIiInIQDOyIiIiIHAQDOyIiIiIHwcCOiIiIyEEwsCMiIiJyEAzsiIiIiBwEAzsiIiIiB+FmdAPsTUBAAG7fvo1o0aLBxcXF6OYQOTXJr/7ixQskSpQIrq48Tw1P7OuI7LOvY2AXStLRJU2a1OhmEJGFGzduIEmSJEY3w6GwryOyz76OgV0oydmr6c2NHj260c0hcmrPnz9XwYfpc0nhh30dkX32dQzsQsk0JSEdHTs7ItvAqcLwx76OyD77Oi5KISIiInIQDOyIiIiIHAQDOyIiIiIHwTV2REQUZv7+/vDz8zO6GURW5e7ujkiRIsEeMLAjIqIw5dW6e/cunj59anRTiCJEzJgx4ePjY/ObtRjYERFRqJmCuvjx48Pb29vm/9gRfc1JzOvXr3H//n31fcKECWHLGNgREVGop19NQV2cOHGMbg6R1UWOHFl9leBOfu9teVqWmyeIiChUTGvqZKSOyFl4f/h9t/U1pQzsiIgoTDj9Ss7ExU5+3xnYEREZ6PHjx6hXr56q7iCLs5s1a4aXL19+9mdKlCih/shYHi1btgx0n+vXr6NSpUpqlEGmjrp27Yr3799b+dUQkdG4xo6IyEAS1N25cwcbNmxQUzxNmjTBDz/8gEWLFn3251q0aIFBgwaZv7ecFpU1cBLUyQ6+3bt3q8dv2LChStkwdOhQq74eIjIWR+yIiAxy5swZrFu3DjNnzkT+/PlRpEgRTJw4EYsXL8bt27c/+7MSyEngZjos67n++++/OH36NBYsWIAcOXKgQoUKGDx4MCZNmoR3797BWT148AA//fQTkiVLBk9PT/W+lStXDrt27TLfR0Y/ly9f/tHPNm7cGNWrVw921NTLywuZMmXC5MmTA/2MvNcjR45E9uzZ1f9X3LhxUbhwYcyZM8e8Tkse1/Q4EninTJkS3bp1w9u3bzF37tyPRmaDHlevXg3z+5EiRQr88ssvH10/YMAA9Xtj+b3p+dzc3NTPdezYMdDIcrt27ZA7d271vlr+rIm8HnmtWbNmVY9h+V5S+GJgR0RkkD179qjp1zx58pivK1OmDFxdXbFv377P/uzChQtVoJAlSxb07NlTpWOwfFz5A5ogQQLzdRLAPH/+HKdOnYKzqlWrFo4cOYJ58+bh/PnzWLlypQrQHj16FKbHk1FTGQ2VIPq7775D69at8fvvv5uDOnnPhw8frkZgZeR0//796j4SvFv+P5QvX149zuXLlzFu3DhMmzYN/fv3x/fff6+uNx0FCxY0P6fpSJo0KSJC5syZ1fNJIDlixAhMnz4dnTt3DnSfpk2bqjYHR0aRZWepBIDyO07Ww6lYIiIDc8HJ+jdLMpoRO3Zsddun1K1bF8mTJ0eiRIlw/PhxdO/eHefOncOyZcvMj2sZ1AnT9596XF9fX3WYSBDoSCQ9y44dO7B161YUL15cXSfvYb58+cL8mKZRU9OolkyfS7BYp04dNRK2fft2HDx4EDlz5jT/TKpUqfDtt98GGjk1jR4KCdQk8JGpeQmgTGk2hIeHR6DnjEjye2l6XgneNm3apF6rBKFiwoQJ5lFR+Z0MKkqUKJgyZYq6LCOkTGxtPQzsrOXtWxljBxo0ACpVMro1RBSBevToof4of2kaNqxkBMhERuYkYWrp0qVx6dIlpE6dOkyPOWzYMAwcODDMbYKmARajhhFG1haGYLdi1KhR1SHTrAUKFFDBVHiTIMwUsMmIqgRolkGdiUy5yhGckydPqtE9CTptmeVrpbA5fRqQbqJyZeDbbxFuGNhZi5y9LFkCrFwJbNwIFCpkdIuIKILIFJWsJ/ocGbmRERBTNnsT2bkqO2VDMyoj6/PExYsXVWAnPyvTfpbu3bunvn7qcWU6t1OnToFG7EI1zSdBXdSoiHCyzitKlBCNOMmaNZnKnDp1KnLlyqVG7mrXro1s2bJ9VRNkmlGmYGWkyhR0X7hwQU3zhsTq1atV0Cn/9zJqKlPxv/76KyKCjPb26dMn0HUSsMmawU85dOiQGp0sVapUBLTQca1bB8yfDzx8yMDOPnTsCGzbBqxdq4fjO3bIIgWjW0VEESBevHjq+BJZMyVTUvKHUhaei82bNyMgIMAcrIXE0aNHA5U6kscdMmSIOUu+kKk92WDxqT/YMoJljVEsW1tjJ7uFZUp27969+Oeff9TmBtm88qVAPDiyWUJ+VgIhqUQgGwpkc4apDFVIlSxZUk1Tvnr1Sq2xkyBU2vo16+GuXbumLhctWlS9zk+RNDhBX7tMq8o0sqUTJ06o4FOCWHm98j5GVPDpqDZu1L+G95JDBnbW4u6OvxuMRYkHbxDrwBZZuQzs3g0kS2Z0y4jIRmTMmFEtnDeNIslOyTZt2qhRJFk/J27duqWmWefPn6/Wg8l0q4yWVKxYUZXzklEiCSiKFStmHnkqW7asCuAaNGigAhdZVyejMrJw32rBm0yJfiH/ntWeNxRkB+s333yjjr59+6J58+Zqo4IpuIkWLRqePXv20c9JAB4jRoyPUtX07t1bTUtKUC0jbSbp0qXD2bNnQ9QmWX+WJk0adXn27NlqF+2sWbNUTsOwWLt2rXnXreUaveDIBhzTc5vIGs+g0qdPr9bUSdApv5uy3o/CTmaxZezHGoGd3eyKlR1M1apVU7+EctYpaQG2bNny0f1kqF06N/nwypmqdGSWpBOUMxi5XaYZpNOzhjFjDqBmnSLIfD0Z3qbPJr2z9Lb6mCsR0QeyFitDhgwqeJNgTfo22XFoIn+gZWOEader/EHduHGjCt7k52TaV0Z3Vq1aZf4ZGT2S6T35KqN39evXV3nsLPPehTtZ5yZTohF9fGU1AAmAZaTMMoCREVRLMkp17NgxFaxZkkBPgqLEiRMHCupMG1zk/0l24QYl/6eWz2lJHqdXr14qEH/z5k2YXpOsz5N2mdoWHuT3Th5PUp0wqPt6e/fqqxdkQD1LFjjniF3lypWRNm1aNU0hZyCy40iuk7NX05qRsWPHYsyYMRg1apSaxpAPjmWOH1kzIp2hLGiVs2MZWpbt2ZJuwHIxcnhwcZGdZ09w59485IjaCaeSPEOkc+f0jRSbNhmzFoWIbI6MjnwuGbH8IbWc1pMT0m2mU/0v/HGXkRvSSUoT2Y0qfb6c/MvInOxYlZN7GTQwkXWGMlImQbOM6snfEUlP8uTJEzW6F1IdOnTAmjVrVMAuOQQlYDc9p2yskRG54PK9CWmnTJFK3sEuXbrAHsj6TslrJ6PDEpCalgdI4GwKBCUtjEzjyhrSFy9emO/zqffBWaZhXcK7UplmBx48eCC9mrZ9+3bzdc+fP1fXbdiwQX3/+PFjLXLkyNrGjRs/+TiTJ0/WYsWKpfn6+pqv6969u5Y+ffoQt+XZs2fqeeXrl3TqNEvdV45COQdrAbHjSPesaWXLappFG4gobELzeaTwe2/fvHmjnT59Wn21F2/fvtV69Oih5cqVS4sRI4bm7e2t+v4+ffpor1+/DnTfhQsXarlz59aiRYumJUiQQKtYsaJ27NixQPcpXry41r59+y8+57Bhw7SsWbNqXl5eWuzYsbXChQtrc+fO1fz8/NR9GjVqpFWrVu2jn5Wfixcvnvby5ctQPWdoJE+eXBs3btxH1/fv31/Lnj37J78PjrTN9PfO8rhy5Uqg5wvuPvbiTTj+3hcooIcDs2eHf19nF+9oQECA+gA2b95c/ZLLB2LUqFFa/PjxVUAnlixZonl6emrz5s3TMmTIoCVOnFj79ttvtevXr5sfp0GDBh99gDZv3qzeLNPjBPfBlDfSdNy4cSNUf0jq1Blu/uWtWeIXTfP21v8369TRNH//r3pfiJwdAzvrcbTAjuhrhdfv/dOnmubqqocC166Ff19nF2vspIyJaa2CDGXL+jiZdpVSPLFixVL3kYzdspNM6iDKNO2ff/6phntlKN2UaycsSTslt5OsozAdoc3yvXBhN5Qs2VFdXra1C9qUHi377gHJTi47Z0Oxc4qIiIjs29atQECAbLCxzn5KV6OTeH6pDp7sKpKRRdkEIZshZJu65GeSOnNVqlRRJU6EBHWyIFW2aUsZF0lAKXmFJJdQcJssQkpyO8kOKdNx48aNMASlo5E1a33JUIVJqzpjSI0J/+W6GzYszG0jIiIi+7LRSmlObGLzREiTeMqGCdnhJYtXTYWuJX+Q5GWSmn8SIJryN1nmaJI8UrKL9vr16+p72WRhStIZ0qSd4ZHbSXY5HTw4G2nSPMKNG/+gzx+9kbD2ODRd3BHo3VvfFhOKRblERERknzY6cmAX0iSepm3+QbeTy/cyUicKFy6svkpagCRJkqjLMhX78OFDc2kW2fYvOYdkZM9UzkWCQ9nebprStRYPD3ecPPkHUqYsg8eP96L5ktFI8O1AVPqjP/Djj5JMCKhe3aptICIiIuPcvAlIekMJZ0qWtM5z2MUaOwnIJPBq1KiRyiUkOe1kK/iVK1dU9msh+YVky3r79u1VnT2ptyf3ly3rktXblFdItl3LVvZTp05hyZIlGD9+fKAyOtYUPXoUHD++Gt7eGaFpt1Bt2e/YX7W9Ptleu/Z/2QqJiIjI4WzapH/NmxeIGdOJAzuZTpWNEpIjR2rT5cmTBzt37sSKFStUhm4Tycwu+esk2JMagDIqJz9nGp2TzQ///vuvCgilfI9MBffr1y/cc9h9TuLEcXDgwHq4uyeBv/9ZFFu/BxdL1wN8fYGqVYFjxyKsLUREROQ407DCRbbGWu/hHY8kOZYAUTZSmNb7hcWOHWdQokQRBAQ8RrSoZXExswvi71svi/2AXbtkcWG4tpvIEYXX55FC996+fftWnSCnTJlSZSkgcgZvv/L3XqItqRQoSThkT2eJEtbp6+xixM4RFS2aEX/8IVnhvfHi5b/IcikWXmXOo/+PS+mxIJs8iIiIyH6dPq3/iZfyvQULWu95GNgZqGbN/Pj117/UHpYHDxcj28t88EuWCrh0CahQQUJ0o5tIRERE4TgNW6yYZNyA1TCwM1jr1uXRq9dcdfnytckoFL0etLjxACkcLbtk3741uolERA7hwYMH+Omnn5AsWTKVxkrSXEne012y/MUi9+jy5cs/+llJzSX5U01KlChhzrcq03KSakvScFmS5PhSi1bWgnt7e6v14pLBYc6cOSo7g+lxTY8j68Flmq9bt25q2m/u3LlfzPVqWQ89tKQOselxIkWKhESJEqnNhZJazN4tXLjQ/L5LOjSpESz1gi1JMQPJiiH156X4QMeOHdX7bi0bNlh/fZ1gYGcDhgyph6ZNf1GXD54cjMrpuwFRo+qT8PXrA/7+RjeRiMju1apVS1Uwkvynkl1h5cqVKkAL+gc/pFq0aKGS5Etx+++++04l0pfE+KagToLG4cOHqw16kq1BkuvLfSZOnKgyM5iUL19ePY5UUBo3bhymTZuG/v374/vvv1fXmw7JEGF6TtMR2mpIQQ0aNEg9juR7lWBo+/btaNeu3Vc9pqnak1F27dqFhg0bmjNg/PHHH+q9l/fOZNGiRSoHrrzPZ86cwaxZs1SmjF69elmlTRLHS8WJiAjs7KJWrLPUpqxQoeeHurKuWtPyv2qah4deTK5lSymYG+7PR2TvWCvWehytVuyTJ0/U69m6detn7yf3+fvvvz+6vlGjRoFqjUvR+/bt2we6T9q0abXatWuryyNGjNBcXV21w4cPf/RY7969U3XPg3tcUbNmTS1nzpwf/Vxwz/k1kidPro0bNy7QdYMHD9YyZcoU6LodO3ZoRYoU0by8vLQkSZJobdu2Nbff9DiDBg1S9dijRYumXtOcOXO0GDFiaKtWrdLSpUunRY4cWatVq5b26tUrbe7cuepnYsaMqR7r/fv35seaP3++ljt3bi1q1KhaggQJtDp16mj37t0z375lyxb1f7Rx40Z1P3ncggULamfPnjXfZ9SoUVqqVKkCvYYJEyaoGvImrVu31kqVKhXoPp06ddIKFy78yffra37vd+zQ/5zHjRu2MvEOVyvWWaxZMwS5czeTAmmYva4z+lYeJ/MCwNSpwMCBRjePiMhuRY0aVR0yzeor6aWsQKb0TKNVMvpVpkwZ5MyZ86P7yZRrlChRgn0MycEqo3uSczWi3bp1C6tWrVJpw0wuXbqkRhRltPP48eNqVEvSjbVp0ybQz44ePVpNfcqIaN++fc3FBaTM5+LFi1Xqsa1bt6JGjRpYu3atOn777Tc1Oim13U1kinrw4MEqZ638X8lUc3AVqqTYwJgxY3Dw4EG4ubmpqVaTggULqvKf8hwSq0uFKXmOihUrmu9TqFAhHDp0SI3kCRktlftb3sca6+tKl9aTE1tV6ONG52btEYJ37/y0VKmqfxi5i65N+n68HubLMWmSVZ6TyF5xxM52RuxkUkEGcSL6CM1kxp9//qnFihVLjTwVKlRI69mzp3bs2LGvHrGTEafffvtN/eyvv/6qrpORpHbt2n2xTfK4kSJF0qJEiaJ5enqqx5CRPmlrRIzYeXh4qOeW90SeO3/+/Gp006RZs2baDz/88NEInrTR9P8vj1O9evVA95ERO3m8ixcvmq/78ccfNW9vb+3Fixfm68qVK6eu/5QDBw6oxzH9jOWIncmaNWvUdZa/j0uXLlWjfm5ubuq2KlWqqJFSS+PHj9fc3d3N92kps2Of8TUjdjIQKH/GZ8zQwoQjdnbM3d0Nx48vQrx4RSVzDdosHYpl3w7Rb5QzpKVLjW4iEdFHpPKjLA2O6ONDxckQkVGn27dvq7V1MgolI0i5cuVSmxTCQjZLyCigjNTJ+i1ZfC+bM0RoUsRKdaSjR49i3759qmJSkyZNVFvDKnPmzOYRygqSYeEzpIqTPLeMxm36UBZBkvz7f1jbLSNn8v6YHk8OWTso5Twlp5uJFA4ISjYupE6d2vx9ggQJ1IYNeQzL6+7fv2/+XkbRqlSpoja4RIsWTRUbEKaa7ybZsmUzXzbVijc9zunTp1UVKilAII8no4Uy8teyZUvzz8j//dChQ9X/4eHDh7Fs2TKsWbNGjRaGN0lwsXdvBK2vM7pWLAUvSpTIOHlyJVKnLo6XL4/j279mY1v1LiiyfLS+mSJ27Ij57SAicjCyg/Wbb75Rh0wZNm/eXC2gN033STAhSWCDevr0qUoQa6levXpqSlACOwkuLOuZS5nLs1IUNARkWjZNmjTq8uzZs9WUpizml8X/YSFTiqZdt9K2z5GduqbnTps2rdopKlOZW7ZsUVPJUvHpxx9/DHZDhQRflq8hKFPVJxPTzt+g15lqvr969UoFjXLIVLbUkpeATr4PuiHD8nHkMYTpcYYNG6Z2H0vQagoCpX1FixbFzz//rP6v5P++QYMG6v9fZM2aVT2/bHSR/9Ogtem/xvbt+h5IeZtTpIDVccTORsWPHxNHjqyDp2dKBARcQqk1m3Dmm8b61poaNYCDB41uIhGRmbc38PJlxB/yvF9D0pTIH3QTSX8hozyWZPRKRq4kWLMkgZ4ERYkTJ/4oEJDa5Bs3blRrzoKSoMvyOS3J48jOzD59+uDNmzdhek3JkydX7TK1LTQk7YkwPbeMaMoImOnxLI/wXgcogbDsUJadxBKESa13y9G8kJK1fUH/P0yvyzSSGpL72FMZMUsM7GxYmjQJsW3bekSKFB9+fkeQZ9c13C5cWe/NZIHnhQtGN5GISJFBExm0iejjw2DNF0nAILXGFyxYoKYdZRpR0mBInrlq1aqZ79epUyfMnDlTTdFduHBBTVPKKI7kdjON7oREhw4d1KhR6dKlMWnSJBUYygL9pUuXokCBAuqxP+Xbb79VQYb8nLW9ePECd+/eVSlPZCOBjHLJSJlsLhDdu3dXmzlks4S8F9JuqdMedPNEeJARQAkWJR2MvFcyZR6WqdEqVaqoqdUpU6aox5H0JzLimC9fPpWrz3QfuV02dsjvwoYNG9QonlxvCvDsNbDjVKyNy58/LVas+AdVqpTA69dbkPV0DVzNXhDRju3RS49JYs0Pv6hERBQ8Wdcluz0lT5zs9JRRM8kBJ2vjLHOX1alTR43YjB07VuU5k3ViuXPnVvndZD1YSEkCZAkWTHnpunTpoh4rY8aMKsjIkiXLJ39WdnlK4CRBp6zZ+9QO2vAg69DkEBLQ5c2bF//++y/ixIljnsbctm2bmp6UUTR5b2TdnOTYC2/y/LKeT/4/ZDetjBbKbtuqVauG6nEaN26sAtZff/0VnTt3RsyYMVVQP2LECPN9ZERUpnDlq+wGlueWoG7IkA9r2sPJnTuApCyUE5CSJREhXGQHRcQ8lWMwquj4rFmb0by5LIJ9h6SJm+Oi5w54XD4nCwP0CfyYMSOsLUTO/nl09vf2a4uhE9mjt2H4vf/tN6BhQ9lcAhw4EDF9Hadi7USzZqUwePBCicVx49ZM5HGvhoAECYETJ2RMWRZEGN1EIiIiMnAaVjCwsyN9+vwPrVrptQhPnBuJb5K3BSRy37kTCMXaDyIiIrIumQ9lYEdfNGlSS9SooVeh2Ly/N+rlGyxbeaTwnX4QERGR4STbze3bkmIHKFw44p6XgZ0d+uuvvihUqLUqTrFoYxd0zddHv6FVK8niaHTziIiInN7GD6N1RYrowV1EYWBnh2Qnz/bt45E+/XeSEQmj94zG+OR1AEmqKUk2PyRpJCIiImMYMQ0rGNjZKcmzc+TIfPj4lJZ83eh0fTfOeqUAtmwBxo0zunlE5ASYVIGciRaK3/f37/U/x+KbbxChGNjZsciRPXHixDJ4eKRCgHYNpTzzQP3aSU6m48eNbh4ROShTOSfJ3k/kLF5/+H0PWhYtOJLa5MULvQJojhyIUExQbOfixo2OOXN+Q716RXHn2Z9onqwFZl2fIUUM9d8s5pgiIivMGEjSV1O5J0m8a6rXSeSII3WvX79Wv+/yex+SyhSmadjSpaVMHCIUAzsHULduIcyf3xPr1w/B7Ot/okHMPChx8iDQuzcwZozRzSMiB+Tj46O+hqWWJ5E9ihkzpvn33lbX1wkGdg5i5cr+iBdvHZ4/P4SqvtFwH27wGjsWqFQJKFXK6OYRkYOREbqECRMifnypZe1ndHOIrMrd3T3ENWSlnPuePfplBnYUZh4e7li1agGKF8+FF2+2oEriH7Dh1nSgUSN9vV2sWEY3kYgckPyxC++i6UT2bPt2QM51UqYEUqWK+Ofn5gkHUqxYBjRvPlpd3nhrHn6LXw64eRNoLTnviIiIyNqMnIYVDOwczPTpPyFx4goAfNH88V08dIkB/P67fhAREZFVMbCjcF/3snnzbLi6xsW798dQKl4N/YaffgJu3DC6eURERA7r7l3gxAn9slHL2xnYOaB06XwwePAMdfnE/XkYnPhDVQpZb8eqFERERFaxebP+NWdOSUcGQzCwc1C9elVH1qxNVT3ZAbd34RyrUhARETn0NKxgYOfAtmz55UNViuso5ZWLVSmIiIisRCqOmQK7iC4jZomBnQOLEycaZs36Tf033366DC2SNgXevQPq1wfevjW6eURERA7jwgV9KbunJ1CkiHHtYGDn4OrXL4SyZXupy7Nu/I3tMXLpKzv79DG6aURERA5j44fRusKFpZa7ce1gYOcEVq7sh+jR8wB4gsrvosMXkQCpSiFr7oiIiMgh1tcJBnZOwNNTr0oBRMaLN1tRNXFTfTGA7JJ9+tTo5hEREdm19+//2xHLwI4iRLFi6dG06Rh1+d9b87Ew/jf6YgBWpSAiIvoqhw7pWcVixgRy5YKhGNg5kZkzWyJRooqqKkXTJ/fwyCUasGgRsHix0U0jIiKyW6ZpWElKbHTpZAZ2TleVYpZelcLvOKtSEBEROdD6OrsK7M6fP49q1aohbty4iB49OooUKYItFov/586dqwKX4I779++b77d161bkypULnp6eSJMmjfo5Z5I+vQ8GDtSrUhy//xuGJv5OX2fHqhRERESh9uoVsHu3fpmBXShUrlwZ79+/x+bNm3Ho0CFkz55dXXdXCrMB+P7773Hnzp1AR7ly5VC8eHHEjx9f3efKlSuoVKkSSpYsiaNHj6JDhw5o3rw51q9fD2fSp49UpWimqlL0vbMXF7yS6Ttkf/nF6KYRERHZlZ079RSxyZIBadIY3Ro7CewePnyICxcuoEePHsiWLRvSpk2L4cOH4/Xr1zh58qS6T+TIkeHj42M+IkWKpILAZs0kgNFNnToVKVOmxJgxY5AxY0a0adMG//vf/zDOCctsbd48Tq9KEXAdJb1y61Upevb8r3oxERERhWoa1sUFhrOLwC5OnDhInz495s+fj1evXqmRu2nTpqmRuNy5cwf7M3Jfb29vFbiZ7NmzB2WCjJPKqJ5c/ym+vr54/vx5oMMRxI0rVSkkBYorbj39Gz8mbaKfctSrJy/a6OYRERHZhY02UEbM7gI7WSe3ceNGHDlyBNGiRYOXlxfGjh2LdevWIVasWMH+zKxZs1C3bl01kmci07YJEiQIdD/5XoK1N2/eBPs4w4YNQ4wYMcxH0qRJ4Sjq1y+IMmX0qhQzbizHjhg5WJWCiIgohB48AI4e/W9HLJw9sJOp1U9teDAdZ8+ehaZpaN26tRqh27FjB/bv34/q1aujSpUqai1dUDICd+bMmUDTsGHVs2dPPHv2zHzccLDdo6tX90O0aB+qUvjFwDv5lRgzRnaZGN00IiIim7b5Q1Li7NmBD8v5Dedm5JN37twZjRs3/ux9UqVKpdbKrV69Gk+ePFE7YsXkyZOxYcMGzJs3TwWIlmbOnIkcOXJ8NE0ra+/u3bsX6Dr5Xh7TcmTPkuyelcORq1KsXLkAJUvmxPPX21A1STOsuzkLaNgQOH5cz7ZIRERENp3mxCYCu3jx4qnjS2SThHB1DTzAKN8HBEnR8fLlSyxdulRNoQZVsGBBrF27NtB1EhzK9c6sRIn0aNJkDObMaYX1NxdgUYIyqHtjI9CmDbBA1uERERGRJanMuWGD7QV2drHGTgIvWUvXqFEjHDt2TOW069q1qzl9iaUlS5aozRX169f/6HFatmyJy5cvo1u3bmqKV0b9JAjs2LEjnN2sWS2RMKFelaLJ43t47BIVWLiQVSmIiIiCcekScO0a4O4OFC0Km2EXgZ0kJZaNEjIaV6pUKeTJkwc7d+7EihUrVD67oJsmatasiZjBTCFKqpM1a9aoUTr5OUl7ItO2sjPW2QWuSnECpeJX129gVQoiIqJPTsMWKgREiQKbYehUbGhIMBeSRMK7TemfP6FEiRJqdy19LEMGH/TvPxP9+1fHsXsLMSzxt+h56w9A1kHKeHOQqXAiIiJntdEG19cJ/qWmQPr1q4YsWfSqFH3u7MNFr6T6tp/x441uGhERkU3w9/9vRywDO7J5W7b8Ag+P1B9XpfhQ5YOIws/jx49Rr149tTtflpBImiZZdvKlmYegqaFkDbGl4NJHLeaaWaJwIRN/T54AMWLIjCJsCgM7+kjcuFExY8Zv6tfj5tPlaJm0oV6NglUpiMKdBHWnTp1Sa38lrdP27dvxww8/fPHnWrRoEag29siRIz+6z5w5cwLdR/J/ElH4TcOWLAm42diiNgZ2FKyGDQuidOne6vL0GyuxM2Z2Pa9d375GN43IYUgiddkYJpu48ufPjyJFimDixIlqZO327duf/VkpmWhZH9uU49OSjABa3keq9hCR466vEwzs6JNWr+6LaNHyAniKyu8+VKUYPZpVKYjCiVTJkeBLNoeZSD1rydG5b9++z/7swoULVcaALFmyqAo5pnyflqRij9wnX758mD17tqriQ0RfRyqQ7txpu4GdjQ0gki3x8nLHihULUKpUTjx7vR3VkzTF2puzWZWCKJxI/WoplWjJzc0NsWPHVrd9itTBTp48ORIlSoTjx4+je/fuOHfuHJYtW2a+z6BBg1R6KBnZ+/fff9GqVSu1dq9du3bBPqavr686TKSGNhF9bNcufVVSkiRAunSwOQzs6LNKlkyHxo3HYO7cn/DPzYVYHL8kat/YwqoURJ8hZQ5HjBjxxWnYsLJcg5c1a1YkTJgQpUuXxqVLl5A6dWp1fV+LZRM5c+bEq1evMGrUqE8GdlKtZ+DAgWFuE5EzTsO6uMDmcCqWvmjWrB+RMKFU+PBF4ycP8QTeelWKJUuMbhqRTZI62BK4fe6QOtiy7u3+/fuBflYq58hOWbktpGR9nrh48eJn73Pz5s1Ao3KWZDr32bNn5uMGE5MTBcsWy4hZ4ogdfZGrqws2bZqJLFmywleqUiSohyP3FkqNNqBwYX08mohCXQdbyiU+ffoUhw4dQu7cudV1mzdvVjWwTcFaSBw9elR9lZG7z91HSjN6enoGe7tc/6nbiEj38KGe6kSULg2bxBE7CpGMGX3Qr99MdfnovUUYkbgm8PSpXpUiIMDo5hHZpYwZM6J8+fIqdcn+/fuxa9cutGnTBrVr11br58StW7eQIUMGdbuQ6dbBgwerYPDq1atYuXIlGjZsiGLFiiFbtmzqPqtWrVI7bU+ePKlG8aZMmYKhQ4eibdu2hr5eInu3ZQsge5CyZAFCMageoRjYUYj1718NmTM3V1Upet05gEueiYBNm4AJE4xuGpHdkt2tErjJGrmKFSuqlCfTp0833+7n56c2Rph2vXp4eGDjxo0oW7as+jmZ9q1Vq5YK5kzc3d0xadIkNSKYI0cOTJs2DWPHjkX//v0NeY1EjmKjDac5MXHRuP89VGSnWIwYMdQalODyRjm6Bw9eInHiHPDzu4Sksarh2pMVcJHpm4MH9VMYogjk7J9Ha+J7S/Qx2Zt0+bKkAwMqydJzG/w8csSOQiVePKlKIbthI+HGkxVolayBvu+7fn1WpSAiIod1+bJ+SKWJYsVgsxjYUag1alQApUrpVSmmXl+F3TGyAMeOAf36Gd00IiIiq5CVR6JgQSBaNNgsBnYUJmvW9EHUqHpViop+seAHF2DUKGDbNqObRkRE5JTr6wQDO/qqqhSAN5693oHqSRrqW4WkKsWzZ0Y3j4iIKNxI8gfTiB0DO3JYpUqlQ6NGY9XltTcXY0n84sD163pVCiIiIgdx7Bjw6JE+BZtXJqtsGAM7+iqzZ/8AH5/KqipFoycP8RSR9VJjS5ca3TQiIqJwrTZRooSkE4JNY2BH4VKVwtU1Hnz9TqGUTw39BqlKceuW0c0jIiJymvV1goEdfbVMmRKgb1+9KsWRu79jVOLqwJMnrEpBRER27+1bYMcO/TIDO3IaAwZURaZMLVRVih53DuGyZ0L9FGfiRKObRkREFGa7d+vBnZRizpgRNo+BHYWbLVvGwt09DQICbqCEd16okibduwOnThndNCIioq+ehnVxgc1jYEfhJn78qJg+/bcPVSlWonWyeno1inofvhIREdmZjXa0vk4wsKNw1bhxAZQsqVelmHJ9DXbHyMyqFEREZJeePNFLoQsGduTkVSnyqaoUlfxisyoFERHZpS1b9Nz7mTIBiRLBLjCwo3AXObI7li/Xq1I8fb0DNZI0YFUKIiKyOxvtbBpWMLAjqyhdOi0aNNCrUqy5uRh/xC+qV6Vo29bophEREYUIAzsiC3Pn/oAECaQqxTs0ePIYz+AF/PYb8McfRjeNiIjos65dAy5cACJFAooXh91gYEdWr0rh4qJXpShtqkrx44+sSkFERHYxWpc/PxA9OuwGAzuyqsyZE6BPn1nq8qG7izE6STV9m1GTJqxKQURENmujHU7DCgZ2ZHWDBlVBxowfqlLcPowrnj56RWVWpSAiIhsUEABs2qRfZmBH9JmqFP5SlSIKq1IQEZHtOnECePAAiBJFn4q1JwzsKEIkSBAVU6dKCpRIuP54FdokraNXo6hfH3j3zujmkROKHTt2qI44ceLgmqymJiKnmYYtXhzw8IBdcTO6AeQ8mjbNj99+64OtWwdi8o21aBA9AwocPapXpRg+3OjmkZN5+vQpfvnlF8SIEeOL99U0Da1atYK/v3+EtI2IjLXRTtfXCRdNeiwKsefPn6s/BM+ePUN0e9omYyPevPFDvHhF8OrVfsSMUgQPXu2Em1RV3roVKFbM6OaRE30eXV1dcffuXcSPHz9E948WLRqOHTuGVKlSwRmwryNn5esrI/rA69f6lGyWLPb1eeRULBlXleLVTtRMyqoUZIyAgIAQB3XixYsXThPUETmzvXv1oC5BAsnsALvDwI4iXJkyaVG//jh1edWNJfgzfmE9E2S7dkY3jYiInNyuXf+tr5MJJXtjN2vszp8/j65du2LXrl149+4dsmXLhsGDB6NkyZLm+xw4cAA9evTAoUOH4OLignz58mHkyJHInj27+T7Hjx9H69at1X3jxYuHtm3bolu3bga9Kuc1b14LbNiwCvfurUaDp09RBl6IOX8+UKUK8L//Gd08cgIrV64M8X2rVq1q1bYQke3Ys0f/WrAg7JLdBHaVK1dG2rRpsXnzZkSOHFktepbrLl26BB8fH7x8+RLly5dXHfDkyZPx/v179O/fH+XKlcONGzfg7u6u5qjLli2LMmXKYOrUqThx4gSaNm2KmDFj4ocffjD6JTpdVYqNG2ciW7asePvuFMr41MbBu4v1qhTyaUqc2OgmkoOrXr16iO4nJ4ncNEHkHDRNn4q158BOdnvZvAcPHsgGD2379u3m654/f66u27Bhg/r+wIED6vvr16+b73P8+HF13YULF9T3kydP1mLFiqX5+vqa79O9e3ctffr0IW7Ls2fP1GPKV/p6ffqsVO+n7OMZnaSSfKY07ZtvNM3f3+imkR3g59F6+N6SMzp/Xv8z5OmpaRahgl19Hu1ijZ3kj0qfPj3mz5+PV69eqdG4adOmqYXPuXPnVveR2+V+s2bNUlO1b968UZczZsyIFClSqPvs2bMHxYoVg4dFUhoZ0Tt37hyeSJmrYPj6+qqRPsuDws/gwVWQIYOMlmrofvsYrnrG16tSTJ9udNOIiMhJp2Fz57a//HUmdhHYyVTIxo0bceTIEZVywMvLC2PHjsW6desQK1YsdR+5fuvWrViwYIGaqo0aNaq6/Z9//oGbmz7jLKkNEsg2Fwum7+W24AwbNkxtMTYdSZMmtfrrdTZbtoz5UJXiJkp459WvlNx2DKIpAm3btg1VqlRBmjRp1CHLOnbs2GF0s4goAu2x8/V1hgd2stFBgrbPHWfPnlXJQWXDg4zQSUe7f/9+tT5GOuE7d+6ox5IRumbNmqFw4cLYu3ev2mSRJUsWVKpUSd0WVj179lR5Y0yHrNej8OXjExVTpuhVKa49WYO2MSvotVxGjDC6aeQk5IRQ1t56e3ujXbt26pATxNKlS2PRokVGN4+IIsgeBwjsDE1Q/ODBAzx69Oiz95G8URLMyaYHmS61TMwnmykkmJMAUaZde/XqpQI9STwqZEpWRvTkttq1a6Nhw4ZqKnX58uXmx9iyZQtKlSqFx48fm0f/PodJO62nRImB2LZtAIAY2IPIKOD1VLZDAxwlJSt/HmXJhmyg6tixY6DrZWZgxowZOHPmDJwN+zpyNi9eADFjSo5L4NYtIFEi2OXn0dBdsZJuRI4veS2ZAj9kirck30uSUdN95HsZ5bO8Xb433adgwYLo3bs3/Pz81C5ZsWHDBrU+LyRBHVnXP//0Rrx4/+DVq3341iM7brzdDvTtC8yda3TTyMFdvnxZzQAEJdOxcsJIRI7vwAE9qEuWzLaCOodcYycBmQRejRo1UiV9TDntrly5oqZaxTfffKNG9GTKVs6uT506hSZNmqj1daZcd3Xr1lUbJ2SUT25fsmQJxo8fj06dOhn8CklEjuyGefNmql/Lm++2YwYyAZLb7sgRo5tGDk7Wzm7atOmj62VtL9fVEjmHPQ4wDWs3eezixo2rNkLIaJtMm8qIW+bMmbFixQpz8uEMGTJg1apVGDhwoAoEZbQuZ86c6ucSJkyo7iPDmP/++68K/mQ3rTxuv379mMPOhtSqlQVp0zbBhQuz0NnNA83ea3Dt2lXfKWuPKcDJLnTu3Fmtqzt69CgKFSqkrpN1unPnzlUnf0Tk+PY4SGBn6Bo7e8R1J9Z36NBt5MmTVibY0cW1AEYF7AXWrgUqVDC6aeTAn8e///4bY8aMMa+nk3V3MjNQrVo1OCP2deRMNE2WhwGy7H/fPiBfPtjt55GBXSixs4sYxYr1x44dg+DumgJPAm4gSuYMwNGjwIfUNUSCn0fr4XtLzuT8ecmHC3h5Ac+e2V4Ou9B8Hu1ijR05nyVLusLVNQH8Aq6ikXtx4NQpbqKgCCHlCZmUnMi57HGAxMQmDOzIJiVMGBW1aw9Sl5f5HcU1RNN3yL58aXTTyAGZNmJFiRJFnRXLZi05pI40d8wTOb49DrK+TnBei2zWrFlN8eef4/Hu3Wl87/kN9t7dAIwZA/Tvb3TTyMHUr19fJUKfPXu2qkZjmTaJiBzfHmcL7FauXBnqB5b0I5K5nSisvLzc0L37SAweXBn7fLdjN3xQaORIQHYxf9jpTBQeJI3SoUOHVE5LInIustri5EknC+ykfFdoyNnuhQsXVNUIoq8xYEBFTJxYCk+fbkZ9zwy4/HqrXkd2xgyjm0YOJG/evKpcIAM7Iuezf7+emDh5cscYMwjxVOzdu3dVrdaQiBYt2te0icjM1dUFEyeOQoMGuXHFdysWIh3qzZ4NtG8PZMlidPPIQcycORMtW7bErVu3VI1pU2Uak2zZshnWNiKyrj0ONA0b4sBOKj6EZlpV1qtwezyFl/r1c6Fv3wa4evU3tHGPjjp+AXDt1k3PbUcUTnWrL126pKrVWM48yLo7+erv729o+4jIevY4WGDHPHahxNxOxti16zqKFJFpsrfo75IHA7SDejWKMmWMbho5wOcxU6ZMKiFxt27dgt08kVzmaJwM+zpyBgEBUt0KePJEn5LNmxd2/3nkrliyC4ULJ0O+fB2wf/9wDHN5ih6aK7yk1NihQzJfa3TzyM5du3ZNbRJLkyaN0U0hoghOTPzkiZ6Y+EOFUrsXosCuZs2aIX7AZcuWfU17iD5p6dIeSJlyJt4FXERz9xJYcHQrsGAB0LCh0U0jOyc1qGVnLAM7Iuechs2Tx/4TE4cqsJPhPyKjJU8eAzVq9MeyZW3x+/tTGAVvJOzdG/j2W4CpdegrVKlSBR07dsSJEyeQNWvWjzZPVK1a1bC2EZH1A7tCheAwuMYulLjuxFivXvkhVqzM8PO7gGKeZbDNdyMwdCjQs6fRTSM7/jy6fmY631k3T7CvI2eQNauew+7vvyW1G2xWhNSKlV1kO3fuVIdcJooIUaK4o127Eerydt9dOIK4wLBhwP37RjeN7FhAQMAnD2cM6oicwbNnehlyR9oRG6bA7tWrV2jatCkSJkyIYsWKqSNRokRo1qwZXr9+bZ1WElkYObI6okUrAuANantmA168AAYONLpZZIcaNmyIv/76S/VrRORc9u8HZM4yZUogQQI4b2DXqVMnbNu2DatWrcLTp0/VsWLFCnVd586drdNKoiBJi0ePHq0un/fdgr+QCpg2DTh3zuimkZ2RzRJDhw5F3LhxUaFCBUyZMkUlKSYix7fHwfLXhTmwk7PbWbNmqU5Q5nnlqFixImbMmIE///zTOq0kCuKHH/IjceLvAWho6REPkOmy7t2NbhbZmX79+qkasVICUTZQLF++HKlTp0bu3LkxaNAgHD161OgmEpGV7GFgp5PpVkngGZSUG+NULEWkefOGAnDHw3f7MMIlJ7BiBbB9u9HNIjuUJEkStGrVCuvXr1drhrt3745z586pNCiSnLhNmzY4ZVqMQ0R2TxIT792rX3b6wK5gwYLo378/3r59a77uzZs3GDhwoLqNKKKULp0KOXK0VZcHuL6Gn1zo0kX/xBKFkdS6/u6777Bw4UIV5M2ePRuRIkXCHtPpPRHZvXPngKdP9UxZjlYKOtSVJ8aPH49y5cqpM9zsH9I0S2JPLy8vdbZLFJGWLOmN9Oln463/ObRyK4oZB3ZIJmOgdm2jm0YOQAK60qVLq4OIHMeeD+dpUkIsSNpK5wvssmTJotajyNns2bNn1XV16tRBvXr1EJlJYimCpUsXGxUr9sXatZ0xx/8ChsETcXv00BMSSY0Yos/ImTPnR3Vhg3Jzc4OPjw+++eYb/Pjjj/BwlPT0RE5sj4OurwtzrVhvb2+0aNEi/FtDFAYLF7ZG3Li/wt//Cup6lMa/1zYBv/6qT8sSfUb1EGQklVx29+/fx88//4wzZ85g8uTJEdI2IrKePQ4c2IWp8sTt27dVYmLp7KTTs9SuXTs4MmZjt02tWy/B5Mky/RoFJ+GBzDE14OJFIE4co5tGDvJ53L59u1p7d/fu3XB7zMePH6Nt27YqfZRUv6hVq5Za7hI1atTP/pys9+vduzf27dunpotz5MihlsKYZk3C+riW2NeRIycmjhVLz2F3755s/oTNs2rliblz5yJlypQqIbHkEhs3bpz5+OWXX76m3URhNmHCd/D2zicptFHbM6e+Kvbnn41uFtmJ33///ZO3de3aVX3NlSsX6tatG67PK0tYZLfthg0bsHr1ahU8/vDDD18M6sqXL4+yZcti//79OHDggNq1a1kWLSyPS+Qs9u3Tg7pUqewjqAs1LZSSJEmi/fzzz5q/v7/mjJ49eyYjnOor2ZZx47ar/xvAVVuLZJrm7q5pFy8a3Syyg89jjBgxtLVr1350fYcOHTQfHx/NGk6fPq3afuDAAfN1//zzj+bi4qLdunXrkz+XP39+rU+fPuH+uEGxryNHNWCAhHWaVq+eZjdC83kMUx672rVrf7ZoNpEROnQoigQJZM1UAJp6JAH8/ICePY1uFtkB2Qwmm8BkiYmJTGUuXboUW7ZsscpzyshbzJgxkSdPHvN1ZcqUUX2rTLEGR5a/yG2SN7RQoUIqp2jx4sUDtTssj0vkTPY48Po6EeroTKZg//jjD+u0hugrzZw5Qu0JuvtuNyYiMyC/q8w/Rl9QqVIltSmiatWqqhKFJCtetmyZCuoyZMhgleeUtXoSoAXdgRs7duxPruO7fPmy+jpgwAC1gW3dunVqiljSsUi2grA+rvD19VXreCwPIkcT4MCJicO8K3bYsGGoXLmy6lCyZs0K9yAJYMaOHRue7SMKlcqV0yFTph9x+vQk9HQDWr0HIsnuWBnR+EJaC3Jusn5Oal8XLlwY8eLFU/WvpZZsaPXo0QMjRsgJxqfJ7tqwMG1Wk7QrTZo0Mads2bRpk0qkLP1zWMnPSqJ5Ikd29qy+ecLb2/ESE39VYCe7r9KnT6++t8wB9aV8UEQRYcmS/siadT5evT+FjpEKYsLu3cDffwM1axrdNLIhnTp1CvZ6CepkFMwyrUloTlg7d+6Mxo0bf/Y+qVKlUrnxZGrV0vv379WOVrktOAkTJlRfM2XKFOj6jBkz4vr16+pyWB5X9OzZM9B7IiN2SZMm/ezrILLnxMRuYUr4ZvtC/bLGjBmjzgy/1HERGSVLlngoVaoXNm/uiSnaDQyGG2J07y7DeQCTy9IHR44cCfZ6GaWToMZ0e2hPWCUwlONLpASjjBDK1G/u3LnVdZs3b1ajcvnz5w/2Z1KkSIFEiRKpOraWzp8/jwoVKoT5cYWnp6c6iBzZHgdfX6eEdmdGggQJtPPnz2vOijvF7MP9+681V9ek6v+qskdJfQvU+PFGN4vCmb1/HsuXL6/lzJlT27dvn7Zz504tbdq0Wp06dcy337x5U0ufPr263WTcuHFa9OjRtT/++EO7cOGC2iHr5eWlXbTYAf6lx3WG95YoOJky6X8OVqzQ7IpVd8W2b98eEydOtE6USRRO4sWLjEaNhqjLa94dwiVEAwYN0vPbEdnQblzZnCGbHypWrIgiRYpg+vTp5tv9/PzU6JxkIzDp0KGDmjbt2LGjqtct6+skX13q1KlD/LhEzujpU+D0af1ygQJwWKGuPFGjRg01rB8nThxkzpz5o80TspPMkTEbu/149y4A0aPnga/vEeTyKoVDbzcD3boBX1jYTs7xeaxZs6ZKuB7Sn5Okv5KIPeiOU0fFvo4czfr1QPnygJwDSWEie2LVyhOSH0k6RMmdFDduXPVElgeRrfDwcEW/fqPV5cNvt2MrEgLjxwPXrhndNLIBK1aswIMHDz5K8RHcIZ2plOd6+fKl0c0mojDa4wzr60KzeUKmAry9vTFnzhzrtogoHPXsWQpjxlTE48dr0dAzNa773gF69wYWLDC6aWQwmaxIly6d0c0gogiyh4FdYDI6V6pUKZXAs1q1airjOZGtkw2NU6aMxPffr8MN352YhXRotnAh0LEj8GHHIDmnsFSUSJw4sVXaQkTWFRCg14h1hsAuxGvsJEeSTF3IIeVrZNGuBHlySKJiZ8F1J/YpbdofcPHiDER3z44nfsfgWqKE5IBg0mI7x8+j9fC9JUdy6pSkwgKiRNE3UdhbDjurrLFLliyZqp24ceNG3Lt3T+3MOnHiBIoWLaqSbcr3sqnC398f1iB5mmSkUEYO5UXJLq+gZ9yyO0zqJ0aLFk0l4uzevbtKzGnp+PHjqs1eXl4q+ebIkSOt0l6yLYsWSUb9KHjudwy9XPMCW7cCa9YY3SwiIooAe5wgMXGYN08IiRqlYPbixYvV4uOpU6eqgE5K3EhiTtlqH96kjJkEaRI8SuJNGTGU60y1D48dO6a29ZcvX14lFl2yZAlWrlypyvtYRrxly5ZF8uTJ1WOMGjVK1VxkGgDHlzdvQhQu3FVdHov7eCW/+l27Skp+o5tGRERWtsdJ1tcp4Z1E7/Dhw9r+/fvD9TEfPHigEvNt377dfN3z58/VdRs2bFDf9+zZU8uTJ0+gn1u5cqVK3Cn3FZMnT9ZixYql+fr6mu/TvXt3lQA0pJi0037dvPlSc3FJqP7//udRXM9SOWWK0c2ir8DPo/XwvSVHkjGj3uWvXKk5/OcxxAOSMoX5JW5ubmo0LHbs2AhPkjNPatPOnz9f1XCUsjfTpk1T+aRMJXN8fX3V9KqlyJEj4+3bt2p0rkSJEtizZw+KFSsGD4uyUuXKlVMFu588eYJYsWKFa7vJtiROHAXffz8Iixe3wF/vjuM6oiBZ//6SoAyIFs3o5hERkRU8eQKcOeP4iYlDPRWbI0cO5MyZU3391JElSxY1FSvB1smTJ8OtkVKrUdb2yRSrrJ+TAE6Kcq9bt84cjEmAtnv3bvz+++9qWvjWrVsYJJUGANy5c0d9lWnboLt5Td+bpnSDkoAxaE4rsl+zZzeBh0cWaHiC773yA1IsnessnV7//v1xjfkNiRzSvg+7YdOkkapEcHghDuyuXLmCy5cvq6+fOi5duoRdu3apQtU//fTTFx9T1r9J0Pa54+zZsyrfVOvWrdUI3Y4dO7B//35Ur14dVapUMQdtsnZO1sy1bNlSjehJfipZc6depGuYlhIqw4YNC5SAWTZckP2KHDkSunTRA7m9b3diL+IBY8YAt24Z3TQykOz2l5JcUoJr0aJF6oSOiBzDHmdaXxeWkmJSuzBoGTGThw8fql2rFy9eVJsbXr169dnHko0Xjx49+ux9ZMetBHMSuMl0qeU237Rp06JZs2aBNkjIy5FgT0byrl69ikyZMqlAMG/evGjYsKEacVu+fLn5/rKzVvLzPX78ONipWOngLTt5+XkJ7pgCwH4FBGiIHbssnj3biFSeRXHJdwfQpIkM5xndNDIwJYfMCEgCdhn1l41atWvXRtOmTVXf4YyY7oQcRdmywIYNwOTJQAjGnJyvpJh0dsHFgpICRdaxiZQpU6pp0S+RaVspVP25Q9bDmQpgBx15k+8DJOugBRnlS5QokVpfJx20BGGyLk8ULFgQ27dvV8GpiRTPlvV7n1pfJ6N/8iZaHmTfXF1d8Msvo+S3BZd9d2ARUgFz58pCUqObRgaSpSYTJkzA7du3MWvWLNy8eROFCxdGtmzZMH78eNWhEpF9CXCixMRhDuwkUXHz5s0DXScjZBLUSSAmIkWKpEbswosEZBJ4NWrUSKU1kZx2Xbt2VdO/lSpVMt9PpmIlt96pU6cwePBgDB8+XHXU0h5Rt25dFSjKKJ/cR1KiSIfdqVOncGsr2YfGjXMgefKG6nIb9xj6yYqkPyGnJ78LcvL37t07dVn6nl9//VWdJEqfQUT24/RpGe3SExNLgmKnENott/fv39cyZMigdezYUX1/69YtLV26dNq3336r+fv7a9Zy4MABrWzZslrs2LG1aNGiaQUKFNDWrl0b6D4lS5bUYsSIoVKc5M+f/6PbxbFjx7QiRYponp6eWuLEibXhw4eHqh1MAeA4tm+/oQFe6v9zoGtOfS/8unVGN4sM+jwePHhQa926tepjEiZMqFIhXbhwwXz7hAkTtPjx42vOgn0dOYLp0/WuvWRJzWk+j6FeYydu3LihKj/UqlULq1evVlOdkpTYNDLmyLjuxLHkzdsbBw8OhWekVHjqfwVeWbPIYisZdja6aRSBn0cpiygbtWQtb4sWLdTGrKD9mawhlg1cQZd/OCr2deQImjYF5swBevUChgyB3bLqGjshUxKyNk2CuXz58qm1bM4Q1JHjWbq0O1xc4sHX/zJ+cC8MnDgBzJtndLMogn333Xdqs9WaNWvUjvvg+jPZGOYsQR2Ro9jjZDtiRYhG7GSNiWxKCEo2NcjmAstOUHaXOjKexTqeGjUmY/ny1nB1iYvb2gskSBRHihPrizLIpvHzaD18b8nePX4sBQ70yw8eyMkZnOLzGKLKE7/88kt4tY3I5syf3wKxY0/A+/fnUNurFLbc3gyMHQv07Wt00yiCfGoDlZzQSkL0NGnSoFq1auFeVYeIrGffh92wadPad1AXWiEK7GQ3KpGjihbNHW3bjsC4cdWx9e1uHEUs5BgxAmjRAvDxMbp5FAEkh93hw4dV1RpJfyRk973MRshu/8mTJ6Nz587YuXOnyo1JRLZvjxNOw4Z4jV1oy2i9ePEirO0hMsSoUVURNWoxAG9RxysbIMm1BwwwulkUQWQ0rkyZMiqHndSWlkPy2H3zzTeoU6eOKlEodaY7duxodFOJKIT2MLD7/Bq7+1JTM4QSJ06syo8R2YtIkVwwcuRodfns2+1YgWTAzJl6EiRyeJIDU3JfWq5dkfUsAwYMwMiRI+Ht7Y1+/fqpgI+IbJ+/v/MlJg7VVKzsr5g5cyaiRo0aoge1rOxAZC9++ikvfv65Dm7f/h0tPOKj6rvrcOneHVi1yuimkZXJgmQ5eQ06zSplD00zFjFjxlRJi4nI9p0+LbOHgIQtTpOYODSBXbJkyTBjxowQP6iPj88n68kS2bK5c4eibNm/8ODdQYx2zYKuq1dLQWGgZEmjm0ZWnoqVurBjxowx14Y9cOAAunTpotKfCKk5nS5dOoNbSkShmYbNl8/50pKGKLCT/E5EzuCbb1IgW7Z2OH58NPq7vEMHAO5dushfeSkya3TzyEqmTZum1s9JLez379+r69zc3NTGsXHjxqnvZROFzFwQke1z1vV1gn+piIJYsqSXrCzFG//zaONeEDh8GFi0yOhmkZXITljZEStr6R49eqR2yMohl6dPn44oH/IZ5siRQx1EZPv2MLAjIpMMGWKhfPl+6vKs95fxCO5A797AmzdGN42sQFKaSCmxp0+fqnXE2bJlU0dI1xQTke0lJj53Tr9coACcDgM7omAsWtQKkSKlgr92D/W8CgPXrwMTJhjdLLKSLFmycCc/kYPYu1f/KktiTZUnnAkDO6JgxIrlgRYthqvL69/uxxlEB4YOlUrwRjeNrODnn39WGyVWr16NO3fuqJ2wlgcR2Y89TjwNKxjYEX3CxIn/Q+TIMo7/Gt975ZRM3cCgQUY3i6ygYsWKOHbsGKpWrYokSZKo3J1ySIoT+UpE9mOPkwd2IdoVaylFihQqLUDjxo1VGhQiR+Xm5oJBg0aja9ciOPF2B9YjEcpNmQK0aaOP8ZPD2CIpbYjI7vk7cWJiExdNsg+Hwi+//IK5c+fi5MmTKFmyJJo1a4YaNWrA09MTzkCmZSQjvSQ0tcxST45JPh0JEtTCgwfLkMgzP2757gNq1gT++svophE/j1bF95bs0fHjQPbsUgMcePLEcXLYhebzGOqp2A4dOuDo0aMqWWfGjBnRtm1bJEyYEG3atFEpA4gciYsLMGOGrLVzw23ffZiE9MCyZcDOnUY3jcLZjh07UL9+fRQqVEjVhhW//fYbdvL/mshu7HHixMRfvcYuV65cmDBhgiqa3b9/f5W4UzK2S56n2bNnqzJkRI6gWrW0yJDhJ3W5u5sr/OWCJC3m77jD+Ouvv1CuXDlEjhxZnaD6+vqq6+XseKhsmiEiu7B7t3NPw35VYCf1YJcuXaoWG3fu3Bl58uRRwV2tWrXQq1cv1KtXL3xbSmSgxYslr110vHp/Bp3d8uqLOP74w+hmUTjuip06daoqnWhZDrFw4cKciSCyI3ucfONEmDZPSCc3Z84c/P7773B1dUXDhg1VyR0pt2Mia+5M9RaJHEH27HFRokQvbN3aA5MCbmEQIiF6z54ynAc4yfpSR3bu3DkUK1bso+tlTYskLiYi2yfZqC5ccN7ExGEesZOA7cKFC5gyZYpahzJ69OhAQZ1ImTKlqrlI5Eh+/70dXF2T4X3AbTTwLAxIQtvJk41uFoUDHx8fXLx48aPrZX1dqlSpDGkTEYUtMXH69EDs2HBaoQ7sJDv7unXr8O233waasrAktRVlVI/Ikfj4REb9+kPU5VW+R3AZUYDBg/WtV2TXWrRogfbt22Pfvn1wcXFRa4cXLlyokhb/9JO+vpKIbBunYcMY2CVPnjy0P0LkMKZPrwtPz1zQ8ALfeeXVg7oherBH9qtHjx6oW7cuSpcujZcvX6pp2ebNm+PHH39UO/+JyPYxsAtjHjvJwi5ntEHJdV5eXkiTJo1KXtykSRM4IuZ2okGDtqB//1Jqieo2xEYxj6fA2bOyBsHopjmd8P48vnv3Tk3JSnCXKVMmRI0aFc6KfR3Zk/fvgZgxgVev9Fx2WbPCoVg1j12/fv3UpolKlSph4MCB6pDLcl3r1q2RLl06NXUhu8uIHFHfviURK1Zl6UrQwCu1RANAr15GN4vCgYeHhwro8uXL59RBHZG9OXlSD+qiRQMyZYJTC/WuWFlMLKkBWrZsGej6adOm4d9//1X5oLJly6Zy3Mm6FSJHIwPWkyePRJ06/+D62z2Yg9Rosngx0LGjnhWT7M6rV68wfPhwbNq0Cffv30dAQMBHa4uJyPanYfPnd97ExGEO7NavX48RI0Z8dL2sTZF8dqaC2rJmhchR1a6dEb17N8fly9PQwd0bjf0AF0lavG2bHvmRXZH1dNu2bUODBg1UJZ3glpsQke3i+rqvCOxix46NVatWoaOMTliQ6+Q209lvNBkPJXJgixYNQIECC/Hc7wR6R8qFoTt2ACtWANWrG900CqV//vkHa9asUQmJicj+MLD7isCub9++ag3dli1b1DoUceDAAaxdu1ZlbhcbNmxA8eLFQ/vQRHYlf34fFCzYDXv29MNo7QH6wAXe3bsDlSoBn0gFRLZJNoWZTkyJyL48eACY0lAWcOLExGHePCHr5mTKQnLVLVu2TB3e3t7qumbNmqn7yJTskiVLrNFeIpuyZEknuLgkgl/ADTSVpMXnz0tOFKObRaE0ePBgtTHs9evXRjeFiMKYmDhjRjlJM7o1djZiJ/VhJa+TjNpJSTEiZ5c0aRR8++1gLF3aDEt9T2I0vJBkwACgfn2pR2V08yiExowZg0uXLiFBggRIkSLFR8nXWS+WyHZxGvYrAjvp7GTXqwR2RKSbM6cR/v77F/j5ncD3XsWx6+E2QDYYDR1qdNMohKpzXSSR3WJg95UJihs1aoQcOXJ8tHnCWTBpJwWnRw/ZLV5eTn+wD9GRz+uVPi2bNKnRTXNo/DxaD99bspfExDI5IqsoJJdd5syAs38eQ715Im3atBg0aBB27dqF3Llzq7V2ltq1axf6FhPZuaFDy2HKlLJ4/vxf1PXKgItvdwF9+gDz5hndNPqM/fv3q34s0icSX/n6+mLFihX47rvvIrxtRPRlJ07oQZ0Ed7LGjsIwYpfyM2WTJPeToyfy5FksfcqcOcfRtGkOABp+RzLUdrkBHDoE5MxpdNMc1td+HiWgu3PnDuLHj6++l8c4evQoUqVKpb6/d+8eEiVKBH9/fzgb9nVkDyZPBlq3BsqWlTy7RrfGTkfsrly58jVtI3JYTZpkQ//+jXHjxhy09oiF799d15MWb9zIpMU2Kuh5bXDnuaE89yWiCMT1deGQ7sSyWPa5c+fwXia4iUj57bfBACLj8btjGBIpG7B5s2S/NbpZ9BVYhYLIdjGwC4fATvI8Sb46yV2XOXNmXL9+XV3ftm1bVWvRWiTdwDfffIOYMWMiTpw4+OGHH/Dy5ctA95G2VKpUSbVNpla6du36UeC5detW5MqVC56enkiTJg3mzp1rtTaT8ylePDFy59ZL6/2MF/CVC1276it8iYgo3Ny/D1y69F+NWApjYNezZ08cO3ZMBUheXl7m68uUKWO1pMS3b99Wjy+B2L59+7Bu3TqcOnUKjRs3Nt9H1sBIUCcjibt378a8efNU0CZJRy2nkeU+JUuWVOtoOnTooGpESv1bovCyZEk3uLjEh6//FfzoURA4fVrqjxndLPqE06dP4/jx4+qQadezZ8+av5d+hohsOzFxpkxAzJhGt8aGaKGULFkybc+ePepy1KhRtUuXLqnLFy5c0KJFi6ZZw7Rp07T48eNr/v7+5uuOHz8uC1/U84q1a9dqrq6u2t27d833mTJlihY9enTN19dXfd+tWzctc+bMgR77+++/18qVKxfitjx79kw9r3wl+pSqVaeo3xNXl9jaXXhoWqZMmmbx+0vh42s/jy4uLqrfkK9BD9P18tUZsa8jW9ejhyyA1bRmzTSH9ywUn8dQj9g9ePDAvIPM0qtXr6y2FkVSDnh4eMDV9b/mRo4cWX3duXOn+rpnzx5kzZpVZY43KVeunNpJYjrrlvvIyJ8luY9cTxSefvutOdzcMiBAe4w6kfLoo3arVxvdLApCRvFlJ798DXqYrnf0nf5E9orr64IX6sAuT548WLNmjfl7UzA3c+ZMFLTSu1uqVCncvXsXo0aNUlOtT548QY8ePdRtkqpAyO2WQZ0wfS+3fe4+Evy9efPmk0Gl3G55EH1J9OhuaNVqpLq8xf8QjiEGMGyYnFwa3TSykDx58hAdRGRbZNnygQP6ZQZ2XxnYDR06FL169cJPP/2kNiaMHz8eZcuWxZw5czBkyJBQPZYEZxIYfu6Q9S6ySUPWzEk9R9kY4ePjo/LpSVBmOYpnDcOGDVO5Y0xHUlYSoBAaM6YyokQpLqcHaOGSUV8Qsn270c0iIrJ7x4/riYllbV2GDEa3xraEOioqUqSI2nggQZ1Mff77779qalamMyWDe2h07twZZ86c+exhShRat25dNeJ269YtPHr0CAMGDFDTwqbbJdiTZKKWTN/LbZ+7jyT7M03tBrdZRBICmo4bN26E6jWS83Jzc0Hnznpd5QPaUVyGN2DFneNkfx4/fox69eqpPkh2/EvGgaC7/YMj/a3MZEjlH/nZYsWKBZp1SJEixUcnydbMWkBk1DSs7Ia18viO3Ql1gmKROnVqzJgx46ufPF68eOoIDdNU6uzZs9WuXEmBImQaWEYM79+/b14DuGHDBtXpZZItMx/us3bt2kCPJ/f53BSypEWRgygs+vYthREj8sLX9wB+QH5sXLcOOHoUyCEVKsjZSVAny0mkH/Lz80OTJk1UKqdFn9lFLUFd+fLl1UnnxIkT4ebmpjIVBJ29kNKPLVq0MH8fLVo0q74WoojE9XWfEZbdGbI79dy5c9qOHTu0bdu2BTqsZeLEidqhQ4fU8/76669a5MiRtfHjx5tvf//+vZYlSxatbNmy2tGjR7V169Zp8eLF03r27Gm+z+XLlzVvb2+ta9eu2pkzZ7RJkyZpkSJFUvcNKe4Uo9Bq1WqZ+p1xQXTtPtw1rXZto5vkMOz583j69GnV9gMHDpiv++eff9RO3Fu3bn3y5/Lnz6/16dPns4+dPHlybdy4cU773pLjS5VK3xG7fr3mFJ6F4vMY6sBOUp2kTJky2BQB1kwL0KBBAy127Niah4eHli1bNm3+/Pkf3efq1atahQoVVNAXN25crXPnzpqfn1+g+2zZskXLkSOHepxUqVJpc+bMCVU72NlRaL1546+5uWVUvzc1kUfT5HNy8aLRzXII4fV5fP36tfbq1atAfYkERuut+Fdj1qxZWsyYMQNdJ/2VnGwuW7Ys2J+5d++eer0TJkzQChYsqNJAFStWTJ1kBw3sEiRIoPpM6e9Gjhz5UV/4JezryFbdu6cHdS4umvb0qeYUnoXi8xjqqdiWLVuad8YmTJgwwsrtzJ8//4v3kd1rQadagypRogSOHDkSji0j+jwvL1d8+213/P57Y6zANbwM0BB19GhgyhSjm0YfVKtWDTVr1lT929OnT5E/f364u7vj4cOHGDt2rNosFt5kzXDQ1FEyrRo7dmzzTv6gTKlXZI3x6NGjkSNHDtU3li5dGidPnkTatGnV7e3atVMVduSxJGG7TNvKlK+8lk+RDABymDADANkq0zSsrLKKEcPo1tieUC85vHDhgtoZmzFjRrXY13LHqBxE9LFJk+rC1TUZ/PEAHZEDmDNH/rIb3SyyKFlYtGhRdfnPP/9Ua3mvXbumgqYJEyZYZbd/WAQEBKivP/74o1qPlzNnTowbNw7p06dX645NOnXqpE5is2XLpoJVySgg6/EsA7egmAGA7AXX14VzYCdnshcvXgztjxE5tVix3FGuXFd1eb7LA/jJH9hffjG6WWRRA9u0uUB2+svonWxGKFCggArwrLHbX3bpy2YvS5JtQHbKmnbyByWzJMK0IcxETrRNdbs/1W/LY1+9evWT92EGALIXDOw+L9RTsW3btlUdl0wVSLoTma6wJGeIRPSxadOaIlmyQXin3URfZMFwmYrt2ZNzCTZA6lAvX74cNWrUULWjO3bsqK6XwEt21ltjt7/sxpdp30OHDplTRW3evFmNykkgFhxJY5IoUSKcO3cu0PXnz59HhQoVPvlckqJKAtXgqgaZMAMA2QM/PyYm/qLQLuBz9pqKXFBMX6NQoSHq9yeKS2rNXz5+w4YZ3SS7Fl6fxz/++ENzd3dXfViZMmXM1w8dOlQrX768Zi3y2Dlz5tT27dun7dy5U0ubNq1Wp04d8+03b97U0qdPr243kU0dUgNb2iy1smWHrJeXl3bxw4ac3bt3q/tIdgCp5b1gwQKVIaBhw4ahahv7OrJFBw/qGydk35Ezld9+Zs3NE1I7kYjCZtq0VsiadTheaZcwBunQVaZj27eX4sdGN82p/e9//1PJ12WDQfbs2c3Xy6YEGcWzloULF6JNmzbqeWRErVatWoHW9EluOxmdk6likw4dOuDt27dqVFGmbaW9kgdP8osKGXVbvHix2mAha+qkSo/cV9bdETnKNGyBAkxM/CkuEt198lb6iOwUk4XFsgYltFM0RCJbth44cWIEYrlkwCPtLFxkSrZlS6ObZZes8Xk0rS1z9s0D7OvIFtWrB0j+7oEDgX79YP8ePwZixw7Xz2OI491WrVoFKnXz+++/49WrV+bvZa1IxYoVQ/pwRE5r0qQOMq6CJ9pZzEMyYNQovaI1GUY2FvTt21d1nLKOTQ653KdPHzVqRkS2waE2Thw9CiRJoq+1/rDjPTyEOLCbNm1aoOkA2W5vWXdVhvxl0TERfV7Roj5ImbKputzLJaYkJ5McG0Y3y6nJprDp06dj5MiRKs+lHHJ51qxZKiccERlPQg5ZDSbpcz+xv8h+yGSpLMORGs/yosJxXjnEjxR0xpYzuERhN3aspD6JhDvacayGDyAF2vmZMozUZp07d646YZWd/XLIZQnsPle3lYgifrQuc2bA7lcH/PEHsH27vr565MhwfWguPSQyQPXqKeHjU1td7uiSEDh2DFi3zuhmOS3ZcCDTr0HJxgMPDw9D2kREDjoN+/o10KWLfrl7dyBZsnB9eAZ2RAYZPLiH+npRO4o9iKmP2pEhZGfq4MGDA1VmkMtDhgxRtxGR8RwmsBs1SnZp6QFdVz1xfXgKVbqTfv36wdvbW11+9+6d6vRMZcQs198R0Zc1a5YFXbpUwbNnq/ATUuOoDMvv3g0UKmR005yCVJewtHHjRiRJksSc7uTYsWOqn5NUJERkLNnDdPCgAwR2168DI0bol6Vm+IeYypDArlixYoGynRcqVMhckNryPkQUMrIAuHv3nujVaxWO4TjOIAoyygd+xQqjm+YUgta2lhxylpw93QmRLZHVKrLPIFYsIF062K+uXfUXUry4JNC0ylMwj10oMbcThSfZ4R41agm8ebMNRZEb23EIOHlSXx1MX8TPo/XwvSVbMnEiIBvUpXLe2rWwT9u2ASVK6DtgDx8GLJKhG5LHjojCn3y+f/qpp7q8E2dwC+7hvkOKQk5qw+7YsUMdcpmIbIPdr6/z99fTm4gffghVUBdaDOyIDDZ0aFm4u+eEhtdoiSx6WvVr14xullORs+EGDRogceLEKF68uDrkcv369dUZMhEZy+4Du5kz9fnkmDFl55xVn4qBHZHBPD1dULeuPmq3Fpfx9L2/JLozullOpUWLFti3bx9Wr16tqujIIZcPHjyo8tkRkXHu3gWuXtXXJefLB/vz5AnQu7d+WWqhxY1r1adjYEdkAyZMqAlX13QIwDO0RTZgxgzgwQOjm+U0JIibPXs2ypUrp9avyCGXZ8yYgVWrVhndPCKnZhqty5LFThMTDxwIPHoEZMoka2+s/nQM7IhsQPTokVC5cjd1eQnu4K3smpLVwhQh4sSJ89EuWSHXxZJteERkGLuehj19Gvj1V/3y+PGAu7vtBXbr1q3Dzp07zd9PmjQJOXLkQN26dfFEhhuJKEymTm0AF5fE8MN9dEcmvTN48cLoZjmFPn36oFOnTrgrcz4fyOWuXbuib9++hraNyNnZbWCnaUCHDvrGiWrVgDJlIuRpQx3YSUcnC43FiRMn0LlzZ1SsWBFXrlxRHSMRhU3ChB4oVqyzujzD5Tn85URJpmTJKnLmzIlcuXKpY+rUqdi7dy+SJUuGNGnSqEMu7969G9OmTTO6qURO6907O05MvGoVsGEDIGUJx4yJsKcNVeUJIQFcJpknBvDXX3+hcuXKGDp0KA4fPqwCPCIKu2nTWiBDhp/xRruJoUiDvtIZtG4tOyyMbprDqV69utFNIKIvkI2kb98CsWPbWWJiX1/ANNjVuTOQOrXtBnZSENtUPkxK8DRs2FBdjh07tnkkj4jCJn36qMiZsx2OHBmA0S4B6HP7NlwWLgSaNjW6aQ6nf//+RjeBiEI4DVuggL4r1m6MGwdcuiRTMUBPPetBRAn1VGyRIkXUlKsUzN6/fz8qVaqkrj9//ryqs0hEX2fKlLYAouC5dhlTkUxPWCxrNIiInIxdrq+7fRv4+Wf9spSJjBbNtgO7X3/9FW5ubvjzzz8xZcoUlcRT/PPPPyhfvrw12kjkVPLnj420afXcaQNcvAGp0cz6seFOZhkePnyoLsvOV/n+UwcRGcMuA7uePYFXr/Rhxnr1IvzpWSs2lFg/kSLCP//cQsWKKQH44Q8kwP/yJAX277ezuQjb/jzOmzcPtWvXhqenp7r8OY0aNYKzYV9HRrtzB0iUSC+9+PRphA98hc2+fXpAZ7ocThmVQ/N5DPUaO+Hv74+///4bZ86cUd9nzJhRLUSWkTwi+noVKiRG4sSNcOvWTHRBPPxPtoVt3gyULm100xyGZbDmjIEbkT0lJraLoC4gAGjXTr/cuLFhZTJCHYmdOnUKVapUwb1795A+fXp13YgRIxAvXjyVoT2L/A8Q0VcbOrQbGjWajWs4iS2IhZLDhzOwszLp3+TE1SRSpEjInDmzoW0iclZ2Nw3722/6zIpEocOGGdaMUK+xa968uQrebt68qVKcyHHjxg1ky5YNP/zwg3VaSeSEGjRIi9ix/6cut0FS2Yb+X0InChc7duxA3rx5zd8XKFBA5beTpOtySL8mu/+JKOLt3m1Hgd2LF0CPHvplSWru42M/gd3Ro0cxbNiwQGV25PKQIUNw5MiR8G4fkdOS5XS9eukdxWmcxFFE0XdYUbiZPHkyGjRoEOi6LVu2qHydly9fRvv27dUmMSKK+MTEhw7ZUWA3ZIiUqwHSpPlvOtZeArt06dKpadig7t+/r7K1E1H46dgxJ6JEkd3mAWiJNJIVXHILGd0sh3Hw4EGUKlUq0HWStil58uRIkSKFCvr2mOaDiCjCyDiR5PiNEwdImxa27eJFPW+dkK8GJ5R3DeluDNMho3Xt2rVT6U5kOlYOudyhQwe11o6Iwo/sBmvbVk9uuQ9ncEVzB0aNMrpZDkP6L9lpZiK7Y30splAk1cmjR48Mah2R87KrxMSdOulDjJLy7UNuXyOFaPNEzJgx4WLxzkqGlO+++858nSljimyqsFx4TERfb8CAohg7thDevduNH5EN/0pqjgEDgA85JCnsokWLhkuXLiFp0qTq+5o1awa6XaZkmeqDKOLZzcaJ9ev1mrCSFURG62wgCg1RYCdrTojIGJ6eLmjUqCdmzKiCjbiEB37+iCcdyOjRRjfN7uXPnx/z589HiRIlgr197ty56j5EFLHsIrDz85P1Mvrltm2BDBlgC0IU2BUvXtz6LSGiTxo7thJmz84Kf/8TaI0sWDptmuys0CtjU5hJecQyZcogTpw46Nq1K+LHj29eMyxLSxYsWIB///3X6GYSOZVbt4AbN/SlKAalgguZSZMAyecbLx7Qrx9sRZgyCj958gSzZs0yJyjOlCkTmjRpwtI7RFYSNaoLatTogT//rIdluIVXL18iyuTJQJ8+RjfNrpUsWRITJ05Ex44dMXbsWDXtKktMJLu7JFz/5ZdfPtpcQUQRM1qXNav0fbBNDx7oS2JMO2JjxoStCHVJse3bt6u1dLLgOE+ePOq6Q4cO4enTpypBcbFixeDIWGaHjHL//nv4+KSHpl3Gj8iAqXEfAteuAd7ecFbh9XmUXJyyCezChQvq+7Rp0+J///ufee2dM2JfR0bp3FlmKYCWLQGbzTb044/A9OlAzpzAgQOSzdxmPo+hDuyyZs2KggULqtxOkpVdyIaJVq1aYffu3Thx4gQcGTs7MtI330zFxo0/wRMJ8Bz34DFxItCmDZwVP4/Ww/eWjFKokD5qJ/vEGjaEbeZiyZ1bdo5KlnOgSBGb+jyGOo/dxYsX0blzZ3NQJ+SyrFWR26xFKlx88803aoeurIeRKhcvX74MdB9Jw5I7d25V1Fuyxgfn+PHjKFq0KLy8vNTZ+MiRI63WZqLwNnVqYwA+8MU9DJS8drKBQhbwEhE5AMldZ9OJiTUNaN9e/1q7doQEdaEV6sAuV65c5rV1luS67Nmzwxpu376tFjhLAuR9+/Zh3bp1qqZjYymyG0TTpk3x/ffffzLiLVu2rEo+KtPHo0aNwoABAzBdhlOJ7EDq1F7Im1ffhTUBvgiQqdglS4xuFhFRuA2GSUq4uHH1Ig42Z+lSfZQucmTARgeGQrR5Qka5LEfFpMyOjM5JXUWxd+9eTJo0CcOlSLkVrF69Gu7u7uo5XGWbjBq5mKrqOEo7TBUvJkyYoL4+ePAgUJtNFi5ciHfv3mH27Nnw8PBQxb2lRJosmmadW7IXU6e2RO7cw/ASNzAeSdFRPnd16+pbyIiI7JhNJyZ+/Rro2lW/3LMnYKNrcEMU2Mm0puwUs1yO161bt4/uV7du3U+Oln0NX19fFYiZgjoRWaJlADt37gxxKTMpDSSbO+SxTMqVK6fSGshOX8v6t5bPLYflqB+RkXLlio6MGVvjzJkh+NnFHR1OnYLLmjWSIdzophEROW7+upEj9TwsyZMDXbrAVoXoFN9UEFu+fu6Q+1iDpBu4e/eumjqVETcJwnr00Iuj37lzJ8SPI4+RIEGCQNeZvpfbgiMl1GTBoulw5l1yZDsmTmwvpzd4rF3GAiSQX1R9zQcRkR2z2cDu2jXAVDZV1jZ/GFyy28BO1qSF9AgNCc5kJPBzx9mzZ9WUqdRwHDNmDLy9vVUtx5QpU6qgzHIUzxp69uypdqGYDkmLQGS00qXjIVmy5upyT8TUe8OdO41ult2Snf2jR49Gvnz5VP8iOTktDyKyvps39UP+rOfNC9vStSvw9i0gVWpq1YItC9FU7MqVK1GhQgW1zk0uf07VqlVD/OSyuza4DRCWUqVKZZ7mlePevXuIEiWKCvpkbZzp9pCQDlt+3pLpe8vC35Zkh60cRLZm5MguqF17Cm7hHP5BLFSQtXZFixrdLLs0cOBAzJw5U/VJffr0Qe/evXH16lUsX74c/WwoozyRM4zWZctmY4mJt24F/vhDjzjHj7fBxX9hCOyqV6+upiql3I5c/hQJtuTMN6TixYunjtAwTZ3KBghJWSIpUEJK8u9Jh+3n56eCVLFhwwakT58+2PV1RLbsu++SoU2benj4cB46ICEqrF0rO530XpFCRTZWzZgxA5UqVVI75evUqYPUqVOrDVqyOUw2jRGRE07D+vvr6U1MSYntoH8N0TxmQECAuYaiXP7UEZqgLrR+/fVXlcvu/PnzandsmzZt1Po3yWtnIjtkZZerBKFv3rxRl+WQdXlCRvxk40SzZs1UupQlS5Zg/PjxKgcfkb2Rk8b+/bvLJZzHaexDlP/WgFCoSJ8hyddF1KhR1bILUblyZayRjSlE5JyB3YwZ+gmzDP4MHgy7oNmJBg0aaLFjx9Y8PDy0bNmyafPnz//oPsWLF5fV4x8dV65cMd/n2LFjWpEiRTRPT08tceLE2vDhw0PVjmfPnqnHlK9ERvP317Ro0Wqo38lcyKxprq6adumS5izC6/OYLl06be/evepy4cKFtWHDhqnLixcv1uLFi6c5I/Z1FJHevtU0Dw/ZAaZpFy5otuHxY02LE0dv1MSJdvN5DNFUrJARsE2bNqkzWNOmAss0IFJ9YvDgwWp61Brmz5//xftslXnwL5CplR2SXJDIAciSj06demLgwL9xGOdwLsAd6WXH1uTJRjfNrtSoUUP1b/nz50fbtm1Rv359zJo1C9evX0fHjnpCaCKynsOH9cTEsjordWrYhgEDgEePgMyZ9cK1diLEtWIlIbBMSaxatUp9Hy1aNLVb1ZRPTnavSm47R+8EWT+RbI1UFIsWrQx8fTehBDJji+dFfWt+kNQ+jshan0dZVye1r9OmTYsqTpofkH0dRaSxY2VDpWzABFasMLo1AE6dAqSaliwx27hRUhEY2hyr1IqVxcVBqzMsWrQIW7ZsUYfkmFsqpTaIKELJPqDmzXuqy9twEbd9/fSdWxRi27dvx/v3783fS1UdWXsr2QDkNiJyovV1mgZ06KAHdTVqGB7UhVaIAzvZmGBaXCxkytUyh5zkfzp9+nT4t5CIvmjEiFJwc8sLDb5ohQzApEnAhw0A9GUlS5bE48ePP7pezo7lNiKyLpsK7Nav10fpJNWZLG2xMyEO7J4+fRpoTZ3UY02RIoX5e9kVa3k7EUWcKFFc8O23+qjdKtzAMyl9N22a0c2yG7IiRdI1BfXo0SOVN5OIrEfy/t+6JWv1gTx5bOJMGUqrVpJMF/YmxJsnkiRJgpMnT6qcb8E5fvy4ug8RGWPixGpYvDgjArQz6IB0mDNuHCD516y0ockR1KxZU32VoE6SpVsmI5f0TdKvFSpUyMAWEjnPaJ0saTP8PGr/fj0hsZub7EyDPQrxiF3FihVVBva3UlIjmB2zkrldknsSkTHixHFFuXKS1w5YiId4I/WPQ7Cb3JmZakDLiJ1sCLOsCy3VaGRd8YIFC4xuJpFDs6lp2FGj9K/16smIFuxRiHfFSumtHDlyqAS/khw4Xbp06vpz586p5MGy8PjIkSPmyhCOijvFyJZdu+aHFCnSALiOTkiFMald5EOqz3E4oPD6PMqJaZcuXTjtaoF9HUWUAgWAffsAOYeSeMowFy8CEttIWHTiBJAlC2yFVXbFSsAm2/8zZsyIHj16qLxPckg+u0yZMmHnzp0OH9QR2brkyd1RsGAXdXkKXuP9pUvAn38a3Syb179/fwZ1RAaQSUDJYWcTI3ZjxuhBncw+2lBQZ7URO0uye0x2yYo0adIgduzYcBY8iyVbd+LEa2TLJhubHmAokqBnjrh6z2njhasj+vOYM2fOYDdMBEfKGTob9nUUEXbvBgoXBqRqqaweMaybundPzowB2QS6bRtQrBjs9fMY4s0TliSQk/QmRGR7smb1RpYs7XHyZB+MgAu6Hz0K13//BcqVM7ppNqV69epGN4HI6VmurzP03PPXX/WgLn9+oGhR2LMwBXZEZNsmTWqN4sVH4BluYBYSoMXw4Qzsgpl+JSJj2cTGiZcv9dyfols3u5/dCPEaOyKyH8WKxUTKlD+py33hDU227+/da3SzbJrk6pw5c6ZaN2xKVixTsLckwRYRhTtZCGYTgd2sWcCTJ0DatEC1arB3DOyIHNTo0VK32RP3cAV/I9Z/STfpI5KvTnb6jxgxAqNHj1ZBnli2bJkK9IjIOomJb9/WU8YZlpjYz08vVCu6dHGIDAIM7IgcVI0aPkiQoIm63AVxgeXLgTNnjG6WTZK6sJKg+MKFC6pcomX+TtaKJbJ+YmJvb4MasXQpcP26vnujYUM4AgZ2RA5KlokMGtRVfcyv4AK2IQowcqTRzbJJBw4cwI8//vjR9YkTJ8Zd2apHROHO8GlYTfuvT2zf3mGq9DCwI3JgzZunQowYtdXlNkiiZwCVs1MKREqJSTqBoM6fP4948eIZ0iYiR2d4YPfvv7IOQ69j9pO+JtkRMLAjcmCurrLJq4e6fBLncey963/rScisatWqGDRoEPxkvc2H2rHXr19H9+7dUatWLas+t2zUqFevnspNFTNmTDRr1gwvZZfeJ1y9elW1L7jjjz/+MN9P2i9lHr29vRE/fnx07dpVVQgispXExEeOGBzYjfwwWteiBRArFhwFAzsiB9e1a1ZEjlxZ5h3wE1IBM2YADx8a3SybMmbMGBVMSQAkta+LFy+ukq9L/dghQ4ZY9bklqDt16hQ2bNiA1atXqzV9UqP2U5ImTYo7d+4EOqQkWtSoUVGhQgV1H39/fxXUvXv3TlUMmjdvHubOnavqfRPZgkOH9H0LUrAqheRTj2gHDwKbN+ubJTrKRjMHIpUnKOSePXsmlTrUVyJ70bHjLvV7C7hrV+Cqaf37a44gvD+PO3bs0CZNmqSNGDFC27Bhg2Ztp0+fVu0/cOCA+bp//vlHc3Fx0W7duhXix8mRI4fWtGlT8/dr167VXF1dtbt375qvmzJlihY9enTN19c3RI/Jvo6sadQoWeCmadWrG9SA777TG1C/vmYPQvN55IgdkRP4+edCcHeXEjl++AlpgYkT9aScFEiRIkXQqlUrdOvWDWXKlLH68+3Zs0dNv+axyPUgz+vq6op9UhU9BA4dOoSjR4+qKVzLx82aNWug+t3lypVT6whldDA4vr6+6nbLg8gh19ddsqih3VU2mDkWBnZETkBSCdStq+djW4/reCAJeGfONLpZNiEgIACzZ89G5cqVkSVLFhUQyZq7+fPny4yGVZ9bdtzK9K8lNzc3VbYxpLtxZ82ahYwZM6JQoUKBHtcyqBOm7z/1uMOGDVO1KE2HTPkSOWRi4rFj5YMPlC8PZMsGR8PAjshJjBtXDq6uOaHhDdrJqN2YMcC7d3BmErhJENe8eXNVYUKCusyZM+PatWsqr12NGjXC9Lg9evT45AYH03H27Nmvbr+sB1y0aFGg0bqwkkTMUmDcdNyQ7LFEViAb8+/cMSgx8YMHwOzZ/5UPc0CsFUvkJGLFckHlyj2wcuX3+BP3MP3mc0RbtAho3BjOSjYUyGaFTZs2oWTJkoFu27x5M6pXr65G7hqGMnFp586dVWD4OalSpYKPjw/u378f6HrZuSo7ZeW2L/nzzz/x+vXrj9onP7t///5A1927d89826dSvshBZG2m0bocOYDIkSP4yX/9Vd+SKxFliRJwRByxI3IikybVgotLWrzHc/RACr3MmExJOKnff/8dvXr1+iioE6VKlVIjbwsXLgz140ruuwwZMnz28PDwQMGCBVX5MlknZxlQyvRw/vz5QzQNKyOOQXPtyeOeOHEiUNAou24lpUqmTJlC/XqIwpNh07CvXumBnWm0TrK4OyAGdkROJEmSSChaVJ9+mIUX8JXpwBUr4Mw1YsvLOptPkPQhx44ds9rzy9o4ef4WLVqoEbZdu3ahTZs2qF27NhIlSqTuI1PEEggGHYG7ePGiGm2UaeSgypYtqwK4Bg0aqPavX78effr0QevWrTkqR84b2M2eLYkjZbgcqFkTjoqBHZGTmTKlAYBE8MUjDEFiYPhwfTWzE5Ipz6CbDCzJbU+ePLFqG2REUAK30qVLq9q0sjN3+vTp5tslafK5c+fUlKsl2fCRJEkSFcQFFSlSJJUTT77K6F39+vXVdK0kYSYy0ps3BiUmfv/+v+TsXbro+esclIvkPDG6EfZEUgDIjjFZXCzTGkT2KGfOsTh6tDOiwgdPcReRJFFnMNORjv55lMBHdol+qmyYrEuTkTNJ+Ots2NeRNezcCRQtKms9gdu3I3A2dPFioE4dWScBXLtmwOK+iPs8cvMEkROaNOkHFC48BC9xF5MQH+1k1M4OA7uvJee1ssnhU9OTktuNiKwzDRthQZ2m/Vc+rG1buwvqQotTsUROqFChqEiTpq26PAge0KQY9uHDcDaNGjVSeeQs87dZHnJbaHfEEtGnGbK+btMmff5XEnq2agVHxxE7Iic1blxbVKkyGo9wE4sQC/Vkh+ySJXAmc+bMMboJRE7DsMTEIz+M1slGozhx4Og4YkfkpCpVioOECX9Ul3sghl5i58IFo5tFRA5KlrZJ4RNJTJw7dwQ96ZEjkutH3yzRsSOcAQM7Iicl61uGDesEwB03cRXrAiIDo0YZ3SwiclCm0bqcOSNwmduoD33a998DKVLAGTCwI3JiDRokRqxYjdTlDvAB5s3Tt6oREYWzCJ+GvXIFWLpUv9y1K5wFAzsiJ+bqCvTuLQmLXXEOl7DvnQvwyy9GN4uIHFCEB3bjxgGSqkhyPUr9MifBwI7IybVvnxbe3v9Tl1sjmWQwBqyclJeInC8x8dGjERjYPXwIzJzpdKN1goEdkZOThczt2vVQlw/hEs6+fKMHd0RE4eTgQb34Q8KEQLJkEfCEkyfr0aQs6CtdGs6EgR0RoV+/nPDwkJqpAWiNlPp0rHSKRET2lpj49Wtg4kT9crduEZgJ2TbYTWB3+PBhfPPNN4gZMybixImDH374AS9fvjTfLoWu69Spg6RJkyJy5MiquPb48eM/epytW7ciV65cKtN8mjRpMHfu3Ah+JUS2R3aoNWrUU13egmu49eCBJHkzullE5CAidH3d3Ln6VKzsgv2fvszEmdhFYHf79m2UKVNGBWL79u3DunXrcOrUKVUKyOTQoUMqS/yCBQvUbb1790bPnj3x66+/mu9z5coVVKpUCSVLlsTRo0fRoUMHNG/eHOvXrzfolRHZjpEjiyJSpELQ4Ie2SKWnCZC5EyKir0xMvHt3BAV2798DY8bolzt31teaOBkXTYol2rjp06ejb9++uHPnDlxlGx+AEydOIFu2bLhw4YIK+ILTunVrnDlzBpulwDmA7t27Y82aNTh58qT5PrVr18bTp09VsBgSLIxNjux//1uNv/6qAld44wFeI/aCBUC9erBV/DxaD99bCi+XLwOpUwPu7vJ7BXh5WfHJli7Vc9ZJhQnJiBwlCpzt82gXI3ZSiNvDw8Mc1AmZbhU7d+785M/JGxA7dmzz93v27FEjf5bKlSunrv/cc8sbankQOaqJEyvBxSUrAvAa3ZAcGD5cP90mIgqHxMRWDeqkrxr5oXxYmzYOE9SFll0EdqVKlcLdu3cxatQovHv3Dk+ePEGPHvouPhnFC87u3buxZMkStRbPRB4jQYIEge4n30uw9uYTC8WHDRsWqCi4rOEjclQJE7qgVCn9szUfT/FaRrfXrjW6WURkxyJsfd2WLbIuS1803Lo1nJWhgZ0EZy4uLp89zp49i8yZM2PevHkYM2YMvL294ePjg5QpU6qgzHIUz0SmWqtVq4b+/fujrCQm/AqyTk9G/kzHjRs3vurxiGzd5MnfAUgFPzzDACTWR+2IiGw9sDON1jVtCsSLB2dl6KrCzp07B9oAEZxUqVKpr3Xr1lXHvXv3ECVKFBX0jR071ny7yenTp1G6dGk1UtenT59At0lAKD9vSb6X+WrT1G5QsntWDiJnkS6dG/Lk6YqDB3/CJLzFzzt3wkOWPBQpYnTTiMjOvHolWSsiILCTJ5GNkK6uQCepge28DA3s4sWLp47QME2lzp49G15eXioFionshpVp20aNGmHIkCEf/WzBggWxNsi00oYNG9T1RPSfyZMbI1++gXiNuxiHeOg+YgQDOwo/AQGApKt6+lSvciJfLY+g15m+l5Psli2lyLG+Ep/sIjGxVPVKlAiw6kom2cUvvv1WRoTgzOxmH7CkLSlUqBCiRo2qgrGuXbti+PDhKq+dafpVgjrZDNGpUye1nk5EihTJHDy2bNlSPU63bt3QtGlTtVt26dKlaqcsEf0nb14vZMjQEWfPdsdwRELX1avheuIEkDWr0U0jWyOF1m/fDl2A9uyZHtyFxYEDgJy4y4xM/foM8GxchCQmlt2vixc7Zfkwuw7s9u/fr9bMSVLiDBkyYNq0aWggZ20f/Pnnn3jw4IHKYyeHSfLkyXH16lV1WdblSRDXsWNHlbw4SZIkmDlzpgoGiSiwX35pifLlh+Ep7mI2YqK5jNpZfLaIlC5dgGXLwvazHh5ArFiAnKDLYXk5uOvk5ELWUUn+DFlHZRngOWG+MnsQIevrxo3ThwWldFju3HB2dpHHzpYwtxM5k2TJ+uDGjSHwQRLcdr0Nl4sX5QwJtoKfRxt4bzt2BFavDnlwZnldWHJfyKKtqVMBOdGQCilCkqQxwLM5El3I6in5b9q1CyhUyApP8vixXnxWfi9kjV3Zr9swaatC09cxsAsl/iEhZ7Jo0QPUq5ccwBssQxTUaN1Y1kXAVvDz6MTvrfwhnzJFH8GzDPD69tWTajPAM9ylS4DUD7BqYmLTqG327MCRIw5bF9bhEhQTkTHq1ImHOHGaq8tdEBeYNQu4f9/oZhHpyWdlGljW+MnCeVlLLZGEZFrIkAGYN48l8WxkGjZXLisFdZJ/dsIE/XK3bg4b1IUWAzsi+iTpJ/v166KW417GNWx/q/3XkRLZWoAno3dx4/4X4GXMCMyfzwDPUdfXyf+tnGgmT67vhiWFgR0RfVarVskQNapeL7atJCyWqdgXL4xuFtHHAZ7siJQAT9bfSYAna0IbNQIyZQJ++40BniMFdrJZYvRo/bLkrePuaDMGdkT0WbJUqWPH7jJ+h+O4jGPPXurTXES2KGpUfVrOMsC7cAFo2JABXgQvgTx+3IqB3fLleuAu9eCbNbPCE9gvBnZE9EW9emWEp2d1dbk1kuvTsWHNQ0YU0QGelMWLE+e/AC9zZj11j4z6kFVIukF5exMntkJiYtnzKUG7aNVKH60lMwZ2RPRFsvC5efOe6vIuXMNV+QMpqQWI7CHA694dkHympgDv/Hm9eoWM4C1cyADPCqw6Dbt9ux45SiWStm2t8AT2jYEdEYXIkCF5ESlSSVncgq5IBYwfb3STiEIf4MkI3rBh+hSeBHiS+05G8Bjg2U9gJ5tkRJMmQPz4VngC+8bAjohCJEYMoFKl9uryctzHKxmxO3vW6GYRhU60aECPHvoI3tCheoB37tx/Ad6iRQzwwmGm1GqBnVQfkZrvsmW/c+dwfnDHwMCOiEJszJjKUpwP7/ESo5EAmDjR6CYRhT3A69nz4wBPkhszwPsqkm3m4UO9YpzksAtXpp2wtWrp2Y/pIwzsiCjE0qSJhCxZ9DUt4+ECbe5cvag7kb0HeDJFK1UMLAO8LFmA339ngPcViYllGVy4uXFDD7iFbIyhYDGwI6JQGTq0qSxYwhPcxd+vXYDZs41uEtHXkzJNvXrpAd7PP+u1bGWpQd26QNaswOLFDPBCyGrTsL/8oqeqKVECyJs3nB/ccTCwI6JQqVw5BuLEaawu9zFNx/IPHjlSgNe7tz5FawrwzpyR+noM8IwM7J48AaZP1y9ztO6zGNgRUajImuVOndqpy2dwGWeuXgdWrza6WUTWCfBkBG/wYCBmzP8CvGzZgCVLGOAF4+VLKyUmnjpVf3CZHi9fPhwf2PEwsCOiUOvYMS3c3Supy12QgqlPyLG3g/fpo4/gmQK806eB2rX/C/CYrNtM0svJ25EkiX6Ei7dv/+tjZLROzi7pkxjYEVGoRY4M/O9/euqT9biDZ1u26GkIiJwhwBs06OMAb+lSBnjWmoaVMnD37uklLOT9ps9iYEdEYTJiRBkAmeCPNxiMRHqZMSJnCPD69tWnaAcO1L8/dQr4/ns9wPvjD6cO8MI9sJPpblOKk44dAXf3cHpgx8XAjojCJGlSF+TOra+1m4Z38Jez6kePjG4WUcSQEbt+/fQRPMsA77vvnDbAk8TEe/eGc2C3cqVeIUTe7+bNw+lBHRsDOyIKs1GjGgCIhZd4iIW+HsCMGUY3ici4AG/AgMABXvbswJ9/Ok2Ad/GinphYctflzBlOkeKIEfrlVq30nIP0RQzsiCjMSpTwho/PD+ryAMQGJk3S80wROWOA17+/HuDJV9lVe/Ik8O23QI4cwF9/qQDv+XN9OV6DBkDZskDLlvpM4/Ll+jLV169ht0zTsLlzh1Ni4p07gX379AdrqydGpy9zC8F9iIiCJZvTevZsjfbtR+MKruHQzUjI/fff+h8zImcN8GTkrn17fSfnuHG4duIZVv1vC1ZGTYitbwvA7/3nx1QSJdKrZaVOrX81HfK9DAg6zfq6kSP1r40aAT4+4fSgjs9F02Ssk0Lq+fPniBEjBp49e4bockZG5OTevZM/Nt/h7ds/UAIpsaVwIv1MOwLw82g9fG/DTv6qHj6sLw9b+fd7HD0ReAwlvecVVC3vh4xV0uDKNVc1hSnHhQtfrtAXN+7HwZ7pcpw4xmYCkYHJY8f02Wcp5fpVZMex1OyVFyQVQNKlgzN7HorPI0fsiOirSKHvevXaY9asP7ANt/Bg1xXEO3RIn48hchKSak2y/kgwt2oVcOuW6RY3uLoChfP7oWqUzaiypxfSvzoMrABwM7ceBaVIYX6cx4+BS5f0QM/ykOsk44esYZPDtEnBkozmBTfKJ18TJrRu0PfixX8Zj8JlxM60E7ZGDacP6kKLI3ahxLNYoo/dvashUaK80LRD+BFJMLVhKWDePKs/Lz+P1sP39sskwFqzRg/m1q8HXr3677YoUfQCCVWrAhUr6iNt5shNap7KIdFQ/PjAihVAgQJffD65u2XQZ3n55s3P/6y3tx7kBQ385JBEwpEifd17sXkzULo0kCwZcO3a1z2WiopTpgT8/PQINn9+OLvnHLEjoojk4+OCwoXbY+fOhpiPV5j4++9wl/UxCRIY3TSicCWZNyQOk2Bu9+7AG14TJ9YDOTmkTr2XVzAPEDu2nuD4hx+AKlWAo0f1O8uJkOTC+wzZFCrTnXIE9eaNnlovuJE+2c8hmzJkRC24POIy6i5xVHDTuzKYGJLUceG6vk7WJkpQV7Qog7owYGBHROFi1KjvULBgN7zBXUz3i4nW06bpaSCI7Jjkx5WgRa2XWwmcOxf4dgmyTMFcrlyhmO6UYbIdO4C6dfW5W6moIJFYr15hmjOVajCZMulHcOtgZRQt6CifHJcv67fL6wr62oSM5MkoXHDTu6lS6c8broHds2d6XVhT+TAKNU7FhhKnJ4g+LXnyQbh+vT8SIwlu+rzX/5rIcICV8PNoPc783kqt+X//1QO51asD592W0auSJfVATgbcJOj56sixa1e1e1Zp2BCYPj2c8oWE7OllGje4kT75KiOBX4pPJdA7eFCfipbsJPnyfUWDZKS/e3c9QpXhRVmgSAjN55GBXSg5c2dH9CUzZ95Dixbyl+4dtsIdxRfMkZ0VVns+fh6tx9neW1nWJQNnEsxt2qSPYpnEigVUqqQHc+XK6Snqwp2MUrVpo0daMgUpaYNkm6uBJDq4cyf4kT45JCefJVlXKEsIw3wu5+urzwnLk86ZAzRuHB4vwyEwsLMiZ+vsiEJDchPHjNkYr17NQz6kxL588fRTeCvh59F6HP29lb98kprDNMUqG7ktyShUtWp6MFe4MOAWEQuXZJhQckBKxCRznTJcmD49bPX9k5FMyyldeZ9kA0WYzZ4NNGumJ/KTBYNWHO23N9w8QUSGkD9+zZq1x4QJ87AfN3Bz/xUkkV1tIdjxR2RtMgq3bdt/wdz16//dJsva5NfUtF4uY0YDcsJJKQrZkVG5sh4tyYK1Zcv0zRU2Rt4b2ekrR7h8vGUXyqhR+uWOHRnUfQVOXhNRuOrfPydcXYvK+B16IikwYYLRTSInJlODCxfqG04lCJHY6ddf9aBOFv7LqNysWfrsn8RUPXroy7sMS/QrSXlllFuipSdP9AbLtKSjk9FJSUQso1GyY5jCjIEdEYUryeZQunQHdXkpnuKNFMb8L1srkdXJejDZiyCbHCRNXP36en1WyQMnlalatNDX08lUotRobdrUxjLzSKMlMZxEo5L2Qxoou2Utc6s4GlP5sJ9+stIiRufBqVgiCnejR1dD9uzJ8Q7XMME/NrrLwvDBg41uFjkoiXf27/8vv5xUo7KUJYs+vSqjc3ny2MlGSxlOXLRIr7ogn51hw/SaY5LvTrINO5Jdu/RDpl/btTO6NXaPgR0Rhbts2SIhbdo2uHChK0bDA92mToVL796fyNhKFHqScHfjxv9KeN2/Hzj3WvHi/6UkkXxrdkkiUElmLBspmjfXy49JCiF50TL06ChMa+saNNA3TtBXsYfzFiKyQ4MGNZNCRniIu1j78AXw++9GN4ns3N27klJHD9gkE4hpfZwEdTJ7Jzl+ZZDrwQM9ZUn79nYc1FmS3HYSxco6hwMH9GoMwZWQsEeyrk6GWkWXLka3xiEwsCMiq/juu1iIHr2RutwLCfVNFMyuRKEgSzMXL9bTu2XPrg/mmNbHvX0rCbH1mTuJeSSYk3OHOnX0vHMOp1gxfVOFTM3Kzg/JLfLPP7B7o0frXyVKz5DB6NY4BE7FEpHVZpHatGmHoUOn4Diu4cLRq0grJZTkDxRREBLzS0kr+RXZuVP/KqnMgsqb97+UJFmzGrh71QgyJSu1u/73P2DLFj0titRVlcjXHt2+Dfz2m36Z5cOcb8Tu8OHD+OabbxAzZkzEiRMHP/zwA15K3ZcPHj16hPLlyyNRokTw9PRE0qRJ0aZNG5XUz9LWrVuRK1cudZ80adJg7ty5BrwaIufQvXsGRIpUXv5sozuSMfVJMB4/fox69eqppKPSvzVr1ixQ3xbU1atX4eLiEuzxxx9/mO8X3O2LZfjLhpJZy6zi2LFAjRr6RlDJHSeZLubP14M6OTmQ+qsyKicvTaZiZZNEnz6yjtPJgjoTmY5dt07fKSu7Rtq21d8gqVhhb6Q/kOSCMvpYqJDRrXEcmh24deuWFitWLK1ly5ba2bNntf3792uFChXSatWqZb7P48ePtcmTJ2sHDhzQrl69qm3cuFFLnz69VqdOHfN9Ll++rHl7e2udOnXSTp8+rU2cOFGLFCmStm7duhC35dmzZzKXpL4S0ZdVrfqP+sy4IbL2zMVF065eDbfHdoTPY/ny5bXs2bNre/fu1Xbs2KGlSZMmUL8V1Pv377U7d+4EOgYOHKhFjRpVe/Hihfl+8r7MmTMn0P3evHlj2Hv78qWmbdqkaQMGaFrp0poWJYq0MfDh5aVpxYtrWp8+mibdsh3/t1pfQICmDR/+35tXqZKmPX+u2Q35z40eXW/7ihVGt8bmhebzaBeB3bRp07T48eNr/v7+5uuOHz+uXuSFCxc++XPjx4/XkiRJYv6+W7duWubMmQPd5/vvv9fKlSvnVH9IiCLSuXP+GpBefW56I558EMPtse398ygnmNJ+OSE1+eeffzQXFxd1QhtSOXLk0Jo2bRroOnncv//+27D39uFDTVu+XNM6d9a0fPk0zc3t40AuZkxNq1xZ00aM0LRduzTt7dswN9d5/fmnpkWOrL+h2bJp2rVrml0YNUpvc4YMmmbxt52+/vNoF1Oxvr6+8PDwgKtF8qHIkuMHshZjZ7A/c/v2bSxbtgzFZc/7B3v27EGZMmUC3a9cuXLqeiKyjnTpXJEli56b6lfJOTZ9OvDqldHNsgnS98j0ax5JrvaB9FHS1+0LYY3dQ4cO4ejRo2oKN6jWrVsjbty4yJcvH2bPni0n8p/tZ2XpiuURUvKwkoVjwQLgxx/1yg1S5aF6dWDMGH36VKZekyTRNzdMnqxv6pQEwbIRQpZXyUycp2eIn5JMatXS66RJ+pPjx/UdszLHbctk+lUySIuuXe0ksaD9sIt3s1SpUrh79y5GjRqFd+/e4cmTJ+ghdV8gZWDuBLpvnTp14O3tjcSJE6s1KzNlb/wH8hgJgqQXl++lA3vz5k24d3ZEpBs6tCGAGHiGB/jjqa9e44lUnxRfFpdZcHNzQ+zYsdVtITFr1ixkzJgRhYKsURo0aBCWLl2KDRs2oFatWmjVqhUmTpz4yccZNmyYKjJuOmSdckhIUJYsGZAihZ6GTOL2M2f020xr5mR9/NWr+mZOSUcixQUkaTD/nocT2VEiJwKym0R+b2RAQ2rM2ir5JZCNEwkTAvXqGd0ah2Pox0qCs08tAjYdZ8+eRebMmTFv3jyMGTNGBW0+Pj5ImTKlCsosR/HEuHHj1EaLFStW4NKlS+jUqdNXtTGsnR0R/ady5aiIHbu5utwP8Rw+9UlI+7avJSekixYtCna0rm/fvihcuDBy5syJ7t27o1u3burk+FN69uyJZ8+emY8bN26EqA0PHwI3b0pACuTLB3TurJfpkvQjUgFi2jS9pJekJnHKzQ4RRaJrmcGqUEF+MfSRvBEjbO9zJhs+TL+HHTpwmNYKXGQ+FgZ58OCB2s36OalSpVLTsCb37t1DlChRVMcoI3Kyy+vbb78N9mdlmrZo0aJqWjZhwoQoVqyY2hH7yy+/mO8zZ84cdOjQQXVknxqxk8NERuwkuJP7y/MTUcgMGXIVffqklp4dx+CCbBs3SFHZr3pM+TzKCZetfR5D2rctWLAAnTt3VrMQJu/fv4eXl5fa4VpDtot+xm+//aaCulu3biFevHifve+aNWtQuXJlvH37VmUFCK/39sgRvVa9zABGifLFhyVrkzlvGdAwjc7K7tkpU/RyXbZg9Wq9HEi0aICcPMSIYXSL7EJo+jpD89hJR/Slzigo01SqrBeRzk9SoHxKwIeCyabArGDBgli7dm2g+8g0hVz/KdIBhqQTJKLP69gxBQYOrA4/v2XohqRYJ/m3vjKws1Uh7duk73n69KlaJ5c7d2513ebNm1XflV8ipRBMw1atWjVEzyXr8GLFihXu/VnOnOH6cPS1ZOhURsQlkbGU3pg9W88d89dftpG52TRaJ4sxGdRZh2YnJDXJoUOHtHPnzmm//vqrFjlyZLXr1WTNmjXa7NmztRMnTmhXrlzRVq9erWXMmFErXLjwR+lOunbtqp05c0abNGkS050QRaA6dbapz48rPLWH0v1cvPhVj+co6U5y5syp7du3T9u5c6eWNm3aQOlObt68qVI3ye2WJCOA7J6VXbRBrVy5UpsxY4bqD+V+kgpK+r5+/fo51Xvr9Nas0bSoUfXdp+nTyy+Nse3Zs0dvi7u7pt24YWxb7IzDpTsRDRo00GLHjq15eHho2bJl0+bPnx/o9s2bN2sFCxbUYsSIoXl5eanOsXv37tqTJ08C3W/Lli0qNYA8TqpUqVSep9BgZ0cUdlevBmhADvUZao8Emtahw1c9niN8Hh89eqQCOclDFz16dK1JkyaB8tHJiaq8Rum7LPXs2VNLmjRpoDRQJhLsST8njxklShSVJ2/q1KnB3teR31vSNO3YMU1LmlQPqOLE0bTt241rS82aejsaNzauDXYqNJ9HQ9fY2SNbXdNDZC/y5JmLQ4eaIApi4Um093CXgqCy3iYM+Hm0Hr63DkR2ykoNNkmDImvtZs3Sd7REpPPn9VqwEnKcPAlkzhyxz+9En0duNieiCDV8eG1ZhYZXeIL5LwKAefOMbhKRY5Mcd1u36jVmJYec5KXp1y9id8xKQkN5Pqlvy6DOqhjYEVGEKl3aCz4+P6nLgxFTX+j9YaMTEVmJtzewZInk4tG/HzwYqFsXePs2YkYMTSdwkviQrIqBHRFFKMll1qOHBHbuuIZb2HvhMrB+vdHNInJ8kvd12DB9p6zsnl28WCoAAPfvW/d5JfWKZKcoUAAoUsS6z0UM7Igo4rVs6QNPz+/V5W5IDEjqEyKKGE2aAP/+q6c/kZKaklrn1CnrPNeLF3oNOdNoHbNUWx0DOyKKcJJKrX799uryTtzGHRmxC4dKDEQUQiVL6kFd6tR6vTcpSbdhQ/g/j5T1fPpUz6snGzjI6hjYEZEhBg/OAxeXQtDwHn2R8L9M+UQUMdKnB/buBYoWlW2XejkyqQEXXvz8gLFj9ctdugCRIoXfY9MnMbAjIkNI/e/ChTuoywvwEm/nztXP7Iko4sSNq4/UNWwI+PvLOgm94K9c/lqyhk8KCUvFKNmJSxGCgR0RGWbkSKmFmhS+eIFpr930Rd1EFPFrI+TE6uef9e9llK1mTeDly7A/pqQ2GTlSvyylzby8wqet9EUM7IjIMAULuiFZstbq8nB4Q5Pp2PAYKSCi0JFNDb176ylRJNBbuVKfopURt7BYt05PRBw1qj4KSBGGgR0RGapfvxYAIuMu7mLT1ZvA6tVGN4nIeX33nZ7MOF484OhRfcfs4cOhfxzTaN0PP+i7bynCMLAjIkM1bBgb3t76+psesomCqU+IjCX55vbtAzJlAm7f1kfuVqwI+c/v368Hh5Irr4O+jpYiDgM7IjKUuzvQvHk7dfkQbuHKli3AiRNGN4vIuaVMCezeDZQtC7x+DdSooa+9C0kZslGj9K9S2SJpUqs3lQJjYEdEhuvbNzNcXb8BEIBeSKSXGSMiY8WIAaxZo6+Rk4BOdsv+9JOexuRTLl4E/vrrvxQnFOEY2BGRTWRcKF1aT1j8F57i5W+/AY8eGd0sIpLpVKkcMW6cvsFC8txVqvTp1ERjxuhBYMWKQNasEd1aYmBHRLZi5MgKANLCD68x3jcyMGOG0U0iIiEBnayVk3V2UaLoee+kUsWVK4Hvd+8eMGeOfrlrV0OaSgzsiMhG5MjhijRp2qrLY+GOgF9/Bd6/N7pZRGRSpQqwcyeQODFw5oy+Y1bW4ZnIZ9bXF8ibFyhe3MiWOjUGdkRkMwYNagwgOh7jAVY8eqnnwSIi25Ejh77rNVcu4MEDoFQp4Pff9WTGkybp9+nWTR/lI0MwsCMim/Hdd9EQPXpTdblPirz6HxEisi2JEgHbtwPVq+sjdLL7VdbUPXkCpE6t76AlwzCwIyKbITXC27SR6VgXnD67EWfOnDW6SUQUHFlrJ7tfTTtfd+zQv8r38kEmwzCwIyKb0rVrKri5yQ7ZSXj4MInRzSGiT3F11XPWyU5ZCeZSpAAaNTK6VU7PzegGEBFZihkT+O23cciWTU98T0Q2TsqGVa2qZxuPHNno1jg9BnZEZHNq1za6BUQUKj4+RreAPuBULBEREZGDYGBHRERE5CAY2BERERE5CAZ2RERERA6CgR0RERGRg2BgR0REROQgGNgREREROQgGdkREREQOgoEdERERkYNgYEdERETkIBjYERERETkIBnZEREREDoKBHREREZGDYGBHRERE5CAY2BERERE5CDejG2BvNE1TX58/f250U4icnulzaPpcUvhhX0dkn30dA7tQevHihfqaNGlSo5tCRBafyxgxYhjdDIfCvo7IPvs6F42nuqESEBCA27dvI1q0aHBxcYGtRPLS+d64cQPRo0eHPXOk1yL4eqxLui/p6BIlSgRXV64sceS+ztZ+974WX49te27HfR1H7EJJ3tAkSZLAFskvny38AoYHR3otgq/HejhS51x9nS397oUHvh7bFt0O+zqe4hIRERE5CAZ2RERERA6CgZ0D8PT0RP/+/dVXe+dIr0Xw9RCFD0f73ePrsW2edvx6uHmCiIiIyEFwxI6IiIjIQTCwIyIiInIQDOyIiIiIHAQDOxs1bNgw5M2bVyUHjR8/PqpXr45z584Fus/bt2/RunVrxIkTB1GjRkWtWrVw7969QPe5fv06KlWqBG9vb/U4Xbt2xfv372Gk4cOHq4SnHTp0sNvXcuvWLdSvX1+1N3LkyMiaNSsOHjxovl2Wrvbr1w8JEyZUt5cpUwYXLlwI9BiPHz9GvXr1VI6kmDFjolmzZnj58mWEvxZ/f3/07dsXKVOmVG1NnTo1Bg8eHKh0jT29HrIfjtzPCfZ1ttU3+DtLXyebJ8j2lCtXTpszZ4528uRJ7ejRo1rFihW1ZMmSaS9fvjTfp2XLllrSpEm1TZs2aQcPHtQKFCigFSpUyHz7+/fvtSxZsmhlypTRjhw5oq1du1aLGzeu1rNnT4Nelabt379fS5EihZYtWzatffv2dvlaHj9+rCVPnlxr3Lixtm/fPu3y5cva+vXrtYsXL5rvM3z4cC1GjBja8uXLtWPHjmlVq1bVUqZMqb1588Z8n/Lly2vZs2fX9u7dq+3YsUNLkyaNVqdOnQh/PUOGDNHixImjrV69Wrty5Yr2xx9/aFGjRtXGjx9vl6+H7Iej9nOCfZ3t9Q1DnKSvY2BnJ+7fvy+nFNq2bdvU90+fPtXc3d3VL6bJmTNn1H327NmjvpcOwdXVVbt79675PlOmTNGiR4+u+fr6RvhrePHihZY2bVptw4YNWvHixc2dnb29lu7du2tFihT55O0BAQGaj4+PNmrUKPN18ho9PT2133//XX1/+vRp9foOHDhgvs8///yjubi4aLdu3dIiUqVKlbSmTZsGuq5mzZpavXr17PL1kP1yhH5OsK+zzb6hkpP0dZyKtRPPnj1TX2PHjq2+Hjp0CH5+fmqY2CRDhgxIliwZ9uzZo76XrzJsniBBAvN9ypUrp2rgnTp1KsJfg0w/yPSCZZvt8bWsXLkSefLkwbfffqumSXLmzIkZM2aYb79y5Qru3r0b6PVIKZj8+fMHej0yhC+PYyL3lzJO+/bti9DXU6hQIWzatAnnz59X3x87dgw7d+5EhQoV7PL1kP1yhH5OsK+zzb6hkJP0dawVayfFuGWNRuHChZElSxZ1nfzyeXh4qF8wS9IZyG2m+1h2DqbbTbdFpMWLF+Pw4cM4cODAR7fZ22u5fPkypkyZgk6dOqFXr17qNbVr1069hkaNGpnbE1x7LV+PdJSW3Nzc1B+0iH49PXr0UH805A9MpEiR1DqUIUOGqDUkpraa2m8Pr4fskyP0c4J9ne32DT2cpK9jYGcH5Ozv5MmT6szCHt24cQPt27fHhg0b4OXlBUf4AyRna0OHDlXfy1ms/P9MnTpVdXb2ZunSpVi4cCEWLVqEzJkz4+jRo+oPbKJEiezy9ZB9svd+TrCvs21LnaSv41SsjWvTpg1Wr16NLVu2IEmSJObrfXx88O7dOzx9+jTQ/WV3ldxmuk/Q3Vam7033iQgy/XD//n3kypVLndnIsW3bNkyYMEFdlrMhe3ktQnZLZcqUKdB1GTNmVDvZLNsTXHstX4+8J5Zk15vstoro1yM77uRMtnbt2moKqEGDBujYsaPasWhqq6n99vB6yP44Qj8n2NfZdt/Q1Un6OgZ2Nko2tkhn9/fff2Pz5s1qe7al3Llzw93dXa0XMJE0AfKBK1iwoPpevp44cSLQL6GcScoW7aAfVmsqXbq0aoecHZkOOQuU4W/TZXt5LUKmioKmZJA1G8mTJ1eX5f9KPuCWr0eG/2X9heXrkc5d/hCYyP+znCHLeo6I9Pr1a7U+xJJMU0hb7PH1kP1wpH5OsK+z7b7htbP0dUbv3qDg/fTTT2rL9datW7U7d+6Yj9evXwfaNi+pATZv3qy2zRcsWFAdQbfNly1bVqUSWLdunRYvXjzD0wAIy51i9vZaJI2Bm5ub2jp/4cIFbeHChZq3t7e2YMGCQFvmY8aMqa1YsUI7fvy4Vq1atWC3zOfMmVOlEdi5c6faRWfElvlGjRppiRMnNqcAWLZsmUqv0K1bN7t8PWQ/HL2fE+zrbKdvaOQkfR0DOxslMXdwh+R8MpFftFatWmmxYsVSH7YaNWqoTtHS1atXtQoVKmiRI0dWv8CdO3fW/Pz8NFvr7OzttaxatUp1vrINPkOGDNr06dMD3S7b5vv27aslSJBA3ad06dLauXPnAt3n0aNHqjOQPEqSyqBJkyYqTUJEe/78ufq/kD82Xl5eWqpUqbTevXsHSq1gT6+H7Iej93OCfZ3t9A3PnaSvc5F/jB41JCIiIqKvxzV2RERERA6CgR0RERGRg2BgR0REROQgGNgREREROQgGdkREREQOgoEdERERkYNgYEdERETkIBjYERERETkIBnZEBtq6dStcXFzUUb169XB97BQpUpgfO2jRcSKiiMS+LuIwsCMz0wfjU8eAAQPgaKRD+OWXX4xuhiq0PXfuXPP3r169Qu3atZEwYULUqVNHFa+2dPfuXbRt2xapUqWCp6cnkiZNiipVqgQqXn3gwAH89ddfEfo6iOwB+zrjsK+zPgZ2ZHbnzh3zIR1A9OjRA13XpUsX2AOpkvf+/fsIfc5379591c/Hjx8fMWPGNH8v73/UqFHx77//InLkyIE65KtXryJ37tzYvHkzRo0ahRMnTmDdunUoWbIkWrdubb5fvHjxEDt27K9qF5EjYl8XduzrbB8DOzLz8fExHzFixFBnrpbXLV68GBkzZoSXlxcyZMiAyZMnB/oAyv2XLl2KokWLqg9o3rx5cf78eXU2lSdPHvXhrVChAh48eGD+ucaNG6th+YEDB6oPp3SwLVu2DNR5BAQEYNiwYUiZMqV63OzZs+PPP//8aIj/n3/+UZ2AnNXt3LkTly5dQrVq1ZAgQQL13NKejRs3mn+uRIkSuHbtGjp27Gg+Uxdytp4jR45A7410NnLGG7TdQ4YMQaJEiZA+fXp1/Y0bN/Ddd9+pjks6Gnl+eW9C68mTJ0iXLh2yZs2q3mvL6YVWrVqptu7fvx+1atVS98ucOTM6deqEvXv3hvq5iJwN+zr2dQ5NIwrGnDlztBgxYpi/X7BggZYwYULtr7/+0i5fvqy+xo4dW5s7d666/cqVK5r8OmXIkEFbt26ddvr0aa1AgQJa7ty5tRIlSmg7d+7UDh8+rKVJk0Zr2bKl+XEbNWqkRY0aVfv++++1kydPaqtXr9bixYun9erVy3yfn3/+2fy4ly5dUm3z9PTUtm7dqm7fsmWLeu5s2bJp//77r3bx4kXt0aNH2tGjR7WpU6dqJ06c0M6fP6/16dNH8/Ly0q5du6Z+Tu6TJEkSbdCgQdqdO3fUIfr3769lz5490Psxbtw4LXny5B+1u0GDBqrdcrx7907LmDGj1rRpU+348ePqPahbt66WPn16zdfXN9j32dT2J0+eBLpe3s9UqVJpbm5u6rXfvHnT3GYXFxdt6NChIfp//NTjE5GOfR37OkfDwI5C1NmlTp1aW7RoUaD7DB48WCtYsGCgzm7mzJnm23///Xd13aZNm8zXDRs2TH34LTsN6TRfvXplvm7KlCmqI/H399fevn2reXt7a7t37w703M2aNdPq1KkT6AO9fPnyL76uzJkzaxMnTjR/Lx2YdGSWQtrZJUiQIFAn9ttvv6nXFhAQYL5Obo8cObK2fv36UHdG8vqlA7Z8vH379qn7L1u27Iuv9UuPT0Ts69jXOR43o0cMyfbJ4lYZ6m/WrBlatGhhvl7Wdsg0hqVs2bKZL8u0gJAhdsvr7t+/H+hnZLrB29vb/H3BggXx8uVLNdQvX2Ux7TfffBPoZ2T6ImfOnIGukykQS/KzMtWwZs0atW5G2vvmzRtcv34d4UFel4eHh/n7Y8eO4eLFi4gWLVqg+719+1a9f6Hl6uqqpoUsyckYEVkH+7rgsa+zLwzs6Iuk0xAzZsxA/vz5A90WKVKkQN+7u7ubL5vWcQS9TtaRhPa5pcNKnDhxoNtkfYmlKFGiBPpeFkBv2LABo0ePRpo0adSalf/9739fXPwrnUzQTsXPz++j+wV9PmmrrHtZuHDhR/eVNTXhIW3atOo9PHv2bLg8HhH9h30d+zpHwMCOvkjOPGXR7OXLl1GvXr1wf3w5+5OzS+mMhCyKlQXAsq1dFuVKpyZnnsWLFw/V4+7atUst/K1Ro4a5Mwq6uFfOQv39/T/qmGSLvXR4pg776NGjX3y+XLlyYcmSJWrXlyyMtgZ5P8qVK4dJkyahXbt2H3W4svDYcscZEYUc+zr2dY6Au2IpRGQnl+zWmjBhgtr9JdvO58yZg7Fjx371Y8tZpUx9nD59GmvXrkX//v3Rpk0bdTYpQ/1yNiq7uebNm6eG+Q8fPoyJEyeq7790xrds2TLVUUmHWrdu3Y/OoGX31/bt23Hr1i08fPjQvINMdrONHDlSPZ90LLIL7UvkD0HcuHHV7rAdO3bgypUrahebdEo3b95EeJH2SAedL18+lbvpwoULOHPmjPq/kakdIgo79nXs6+wdAzsKkebNm2PmzJmqg5P1FnJGKUkmZVv+1ypdurTqmIoVK4bvv/8eVatWDZQgdPDgwejbt6/qbCUFQfny5dV0xZeeWzriWLFioVChQiqhpZz9yZmmpUGDBqkz29SpU5unEOQ5JL2BdCqyJka22ockr5WsnZGOM1myZKhZs6Z6HOnEZd1JeJ7VSqJO6fAll1Pnzp2RJUsWtS5HEnZOmTIl3J6HyBmxr2NfZ+9cZAeF0Y0g5yXTBzKkvnz5cjgjOcuVTktyOVljWsHaj09EIcO+jn1dROGIHZENSJIkiSqnE54kkackSSUishXs66yPmyeIDCQ772TdiJBF1OFJ1vCYdrhZa4EzEVFIsK+LOJyKJSIiInIQnIolIiIichAM7IiIiIgcBAM7IiIiIgfBwI6IiIjIQTCwIyIiInIQDOyIiIiIHAQDOyIiIiIHwcCOiIiIyEEwsCMiIiKCY/g/LsVETk0gnS4AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "fig, (ax1, ax2) = plt.subplots(1, 2)\n", "ax1.plot(Temp, sup.dG*4.182/1e3, 'r', Temp, hp.dG*4.182/1e3, 'b', Temp, bm.dG*4.182/1e3, 'k')\n", "ax1.set_xlabel('Temperature [°C]'); ax1.set_ylabel('Gibbs Energy [kJ/mol] ')\n", "plt.legend(['SUPCRT', 'HP11', 'Berman88'])\n", "ax2.plot(Temp, hp.dG*4.182/1e3 - sup.dG*4.182/1e3, 'r', Temp, bm.dG*4.182/1e3 - sup.dG*4.182/1e3, 'b')\n", "ax2.set_xlabel('Temperature [°C]'); ax2.set_ylabel('Delta Gibbs Energy [kJ/mol] ')\n", "plt.legend(['SUPCRT - HP11', 'SUPCRT - Berman88'])\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate olivine solid solutions" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([ 9.2646654 , 23.54844166, 19.7311314 ]),\n", " {'type': 'ol',\n", " 'name': 'Fo85',\n", " 'formula': 'Mg1.70Fe0.30Si1O4',\n", " 'MW': 150.1557,\n", " 'min': ['Mg1.70Fe0.30Si1O4',\n", " ' R&H95, Stef2001',\n", " np.float64(-466862.12264868076),\n", " nan,\n", " np.float64(26.21037219328013),\n", " 44.049,\n", " 24.058078393881452,\n", " 0.017393236137667304,\n", " -890893.8814531548,\n", " 171.38145315487571,\n", " -3.6586998087954105e-06],\n", " 'spec': ['H+', 'Mg++', 'Fe++', 'SiO2(aq)', 'H2O'],\n", " 'coeff': [-4, 1.7, 0.30000000000000004, 1, 2],\n", " 'nSpec': 5,\n", " 'V': 44.049,\n", " 'source': ' R&H95, Stef2001',\n", " 'elements': ['1.7000', 'Mg', '0.3000', 'Fe', '1.0000', 'Si', '4.0000', 'O']})" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "calc = calcRxnlogK( X = 0.85,T = np.array([300, 400, 450]), P = np.array([200, 200, 200]),\n", " Specie = 'olivine', dbaccessdic = ps.dbaccessdic, densityextrap = True)\n", "calc.logK, calc.Rxn" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate glauconite mineral properties" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([ 7.33889672, 6.08136856, 5.34485756, 4.59351803, 3.11784876,\n", " 1.73061272, 0.45134995, -2.30826164]),\n", " {'type': 'Smectite',\n", " 'name': 'Glauconite',\n", " 'formula': 'K0.68Na0.21Fe0.08Ca0.04Mg0.34Al0.69FeIII1.18Si3.654O10(OH)2',\n", " 'MW': 426.25271512,\n", " 'min': ['K0.68Na0.21Fe0.08Ca0.04Mg0.34Al0.69FeIII1.18Si3.654O10(OH)2',\n", " 'B2015, B2021',\n", " np.float64(-1173001.1408482206),\n", " np.float64(-1256479.2241755025),\n", " np.float64(90.78782850627151),\n", " np.float64(141.34069593680297),\n", " np.float64(77.70879420984171),\n", " np.float64(83.36446524554509),\n", " np.float64(-17.841879278113826)],\n", " 'struct_formula': '(K0.679Na0.212Ca0.036)(Mg0.337Fe0.079Al0.341FeIII1.175)(Si3.654Al0.346)O10(OH)2',\n", " 'spec': ['H+',\n", " 'K+',\n", " 'Na+',\n", " 'Fe++',\n", " 'Ca++',\n", " 'Mg++',\n", " 'Al+++',\n", " 'Fe+++',\n", " 'SiO2(aq)',\n", " 'H2O'],\n", " 'coeff': [np.float64(-7.381),\n", " 0.679,\n", " 0.212,\n", " 0.079,\n", " 0.0359,\n", " 0.337,\n", " 0.687,\n", " 1.175,\n", " 3.654,\n", " 4.692],\n", " 'nSpec': 10,\n", " 'V': np.float64(141.34069593680297),\n", " 'dG': np.float64(-4907.836773308955),\n", " 'dHf': np.float64(-5257.109073950302),\n", " 'S': np.float64(379.85627447024),\n", " 'Cp': 345.17439398884756,\n", " 'source': 'B2015, B2021',\n", " 'elements': ['0.0359',\n", " 'Ca',\n", " '0.6790',\n", " 'K',\n", " '0.2120',\n", " 'Na',\n", " '1.2540',\n", " 'Fe',\n", " '0.3370',\n", " 'Mg',\n", " '0.6870',\n", " 'Al',\n", " '3.6540',\n", " 'Si',\n", " '2.0000',\n", " 'H',\n", " '12.0000',\n", " 'O']})" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Temp = np.array([0.01, 20, 30, 40, 60, 80, 100, 150])\n", "calc = calcRxnlogK(T = Temp, P = 'T', Specie = 'clay', dbaccessdic = ps.dbaccessdic, densityextrap = True,\n", " elem = ['Glauconite', '3.654', '.687', '1.175', '0.079', '0.337', '0.679', '0.212', '0.0359', '0'])\n", "calc.logK, calc.Rxn" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Create new reactions and calculate equilibrium constants" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### An example with pyrite, pyrrhotite, magnetite (PPM) reaction\n", "Pyrite + 4 H2O + 2 Pyrrhotite → 4 H2S(aq) + Magnetite\n", "\n", "Then we can include the reaction in sourcedic, one of the output of db_reader, which is a dictionary of list of reaction coefficients and species. An example with the format for sourcedic is as follows:\n", "\n", "> ps.sourcedic['Name'] = ['formula', number of reactants in the reaction, 'coefficient of specie 1', 'specie 1', 'coefficient of specie 2', 'specie 2']\n", "\n", "> AB → 0.5 A2(aq) + B\n", "\n", "> ps.sourcedic['AB'] = ['AB', 2, '0.5', 'A2(aq), '1', 'B']" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-9.61842122, -8.38742387, -7.21969753, -4.97772897, -4.29319566,\n", " -3.81890541, -3.32234015, -2.52289322])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ps = db_reader(sourcedb = './database/thermo.2021.dat', sourceformat = 'gwb')\n", "ps.sourcedic['Pyrite'] = ['', 4, '4', 'H2S(aq)', '1', 'Magnetite', '-2', 'Pyrrhotite', '-4', 'H2O']\n", "Temp = np.array([300.0000, 325, 350.0000, 400.0000, 415, 425.0000, 435, 450.0000])\n", "Press = 500*np.ones(np.size(Temp))\n", "log_K_PPM = calcRxnlogK( T = Temp, P = Press, Specie = 'Pyrite', dbaccessdic = ps.dbaccessdic,\n", " sourcedic = ps.sourcedic, specielist = ps.specielist).logK\n", "log_K_PPM" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Generate GWB thermodynamic database" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "_Note: GWB requires exactly 8 temperature and pressure pairs_" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "# Vectors for Temperature (C) and Pressure (bar) inputs\n", "T = np.array([ 0.010, 25.0000 , 60.0000, 100.0000, 120.0000, 150.0000, 250.0000, 300.0000])\n", "P = 350*np.ones(np.size(T))\n", "nCa = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write GWB using default sourced database, with inclusion of solid_solution and clay thermo properties" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new GWB database is ready for download\n", "CPU times: user 2.95 s, sys: 49.2 ms, total: 2.99 s\n", "Wall time: 3.08 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = T, P = P, cpx_Ca = nCa, solid_solution = 'Yes', clay_thermo = 'Yes', \n", " dataset = 'GWB')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write GWB using user-specified sourced database, with inclusion of solid_solution and clay thermo properties" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new GWB database is ready for download\n", "CPU times: user 2.49 s, sys: 38.6 ms, total: 2.53 s\n", "Wall time: 2.67 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "src_db = os.path.join(os.getcwd(), 'database', 'thermo.29Sep15.dat')\n", "write_database(T = T, P = 175, cpx_Ca = nCa, solid_solution = 'Yes', clay_thermo = 'Yes',\n", " sourcedb = src_db, dataset = 'GWB')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write GWB using Jan2020/Apr20 formatted sourced database" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new GWB database is ready for download\n", "CPU times: user 2.84 s, sys: 31.2 ms, total: 2.87 s\n", "Wall time: 2.9 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "src_db = os.path.join(os.getcwd(), 'database', 'thermo.com.tdat')\n", "write_database(T = T, P = 125, cpx_Ca = nCa, solid_solution = True, clay_thermo = True,\n", " sourcedb = src_db, dataset = 'GWB')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write GWB using Jan2020/Apr20 formatted sourced database with logK as polynomial coefficients, using Tmax and Tmin" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new GWB database is ready for download\n", "CPU times: user 2.63 s, sys: 28.2 ms, total: 2.65 s\n", "Wall time: 2.69 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "src_db = os.path.join(os.getcwd(), 'database', 'thermo.com.tdat')\n", "write_database(T = [0, 300], P = 150, clay_thermo = True, logK_form = 'polycoeffs', \n", " sourcedb = src_db, dataset = 'GWB')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write GWB using Mar21 formatted sourced database with logK as polynomial coefficients" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new GWB database is ready for download\n", "CPU times: user 2.67 s, sys: 33.6 ms, total: 2.7 s\n", "Wall time: 2.76 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "src_db = os.path.join(os.getcwd(), 'database', 'thermo_latest.tdat')\n", "write_database(T = [0, 300], P = 145, clay_thermo = True, logK_form = 'polycoeffs', \n", " sourcedb = src_db, dataset = 'GWB')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write GWB using user-specified sourced database and direct-access database (slop07) and FGL97 dielectric constant" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Duplicate found for species \"AlOH++\" in slop07.dat\n", "Duplicate found for species \"MgCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"CaCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"SrCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"BaCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"BaF+\" in slop07.dat\n", "Duplicate found for species \"Sr(Succ)(aq)\" in slop07.dat\n", "Duplicate found for species \"Sc(Glut)+(aq)\" in slop07.dat\n", "Duplicate found for species \"AMP2-\" in slop07.dat\n", "Duplicate found for species \"HAMP-\" in slop07.dat\n", "Duplicate found for species \"+H2AMP-\" in slop07.dat\n", "Success, your new GWB database is ready for download\n", "CPU times: user 2.41 s, sys: 50 ms, total: 2.46 s\n", "Wall time: 2.47 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "db_db = os.path.join(os.getcwd(), 'database', 'slop07.dat')\n", "src_db = os.path.join(os.getcwd(), 'database', 'thermo.29Sep15.dat')\n", "write_database(T = [0, 400], P = 300, cpx_Ca = 0.5, solid_solution = 'Yes', Dielec_method = 'FGL97',\n", " dbaccess = db_db, sourcedb = src_db, dataset = 'GWB')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write GWB using default sourced database and direct-access database (slop07 with Berman mineral data) and FGL97 dielectric constant" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Duplicate found for species \"AlOH++\" in slop07.dat\n", "Duplicate found for species \"MgCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"CaCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"SrCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"BaCO3(aq)\" in slop07.dat\n", "Duplicate found for species \"BaF+\" in slop07.dat\n", "Duplicate found for species \"Sr(Succ)(aq)\" in slop07.dat\n", "Duplicate found for species \"Sc(Glut)+(aq)\" in slop07.dat\n", "Duplicate found for species \"AMP2-\" in slop07.dat\n", "Duplicate found for species \"HAMP-\" in slop07.dat\n", "Duplicate found for species \"+H2AMP-\" in slop07.dat\n", "Duplicate found for species \"Cordierite\" in berman.dat\n", "Success, your new GWB database is ready for download\n", "CPU times: user 1.9 s, sys: 55.3 ms, total: 1.95 s\n", "Wall time: 2.01 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "db_berman = os.path.join(os.getcwd(), 'database', 'berman.dat')\n", "db_db = os.path.join(os.getcwd(), 'database', 'slop07.dat')\n", "write_database(T = [0, 340], P = 150, Dielec_method = 'FGL97', dbaccess = db_db,\n", " dbBerman_dir = db_berman, dataset = 'GWB')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write GWB using Mar21 formatted sourced database and direct-access database (speq21 with HP mineral data) and DEW dielectric constant" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Duplicate found for species \"Fluorphlogopite\" in supcrtbl.dat\n", "Duplicate found for species \"Arsenic\" in supcrtbl.dat\n", "Duplicate found for species \"Arsenolite\" in supcrtbl.dat\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/noahawolayo/Documents/pygcc_build/src/pygcc/pygcc_utils.py:1924: UserWarning: Warning: Using SiO2(aq) data from Sverjensky et al. (2014) GCA. This thermodynamic model requires that you also add data for Si2O4(aq) to achieve accurate solubility calculations.Please ensure thermodynamic data for Si2O4(aq) are added to both your source GWB or EQ3/6 database AND your source direct-access database.\n", " warnings.warn(\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Success, your new GWB database is ready for download\n", "CPU times: user 1.48 s, sys: 74.2 ms, total: 1.55 s\n", "Wall time: 1.63 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 275], P = 1100, Dielec_method = 'DEW', dbaccess = os.path.join(os.getcwd(), 'database', 'speq21.dat'),\n", " dbHP_dir = os.path.join(os.getcwd(), 'database', 'supcrtbl.dat'), \n", " sourcedb = os.path.join(os.getcwd(), 'database', 'thermo_latest.tdat'), dataset = 'GWB')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write GWB using user-specified sourced Pitzer database and default direct-access database along the saturation curve" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new GWB database is ready for download\n", "CPU times: user 869 ms, sys: 22.2 ms, total: 892 ms\n", "Wall time: 931 ms\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 350], P = 'T', dataset = 'GWB', sourcedb = os.path.join(os.getcwd(), 'database', 'thermo_hmw.tdat') )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write GWB using user-specified sourced EQ3/6 Pitzer database and default direct-access database" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new GWB database is ready for download\n", "CPU times: user 749 ms, sys: 10.5 ms, total: 759 ms\n", "Wall time: 773 ms\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 350], P = 250, dataset = 'GWB', sourcedb = os.path.join(os.getcwd(), 'database', 'data0.hmw'),\n", " sourceformat = 'EQ36')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Generate EQ3/6 thermodynamic database" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write EQ3/6 using default sourced database" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new EQ3/6 database is ready for download\n", "CPU times: user 3.03 s, sys: 128 ms, total: 3.16 s\n", "Wall time: 5.09 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = T, P = P, cpx_Ca = 1, solid_solution = 'Yes', clay_thermo = 'Yes', \n", " dataset = 'EQ36')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write EQ3/6 user-specified sourced database" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new EQ3/6 database is ready for download\n", "CPU times: user 2.84 s, sys: 90.1 ms, total: 2.93 s\n", "Wall time: 3.09 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 400], P = 350, cpx_Ca = 0.5, solid_solution = 'Yes', clay_thermo = 'Yes',\n", " sourcedb = os.path.join(os.getcwd(), 'database', 'data0.geo'), dataset = 'EQ36', sourcedb_codecs = 'latin-1')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write EQ3/6 using user-specified sourced Pitzer database" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new EQ3/6 database is ready for download\n", "CPU times: user 717 ms, sys: 22 ms, total: 739 ms\n", "Wall time: 746 ms\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 350], P = 200, sourcedb = os.path.join(os.getcwd(), 'database', 'data0.hmw'), dataset = 'EQ36')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write EQ3/6 user-specified sourced database using FGL97 dielectric constant" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new EQ3/6 database is ready for download\n", "CPU times: user 2.76 s, sys: 31.9 ms, total: 2.8 s\n", "Wall time: 2.88 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 400], P = 300, cpx_Ca = 0.1, sourcedb = os.path.join(os.getcwd(), 'database', 'data0.geo'), \n", " dataset = 'EQ36', solid_solution = 'Yes', Dielec_method = 'FGL97', clay_thermo = 'Yes')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write EQ3/6 user-specified sourced database using DEW model" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new EQ3/6 database is ready for download\n", "CPU times: user 1.39 s, sys: 29.5 ms, total: 1.42 s\n", "Wall time: 1.45 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "Temp = np.array([50, 100, 150, 300, 450, 500, 600, 700])\n", "write_database(T = Temp, P = 1500, sourcedb = os.path.join(os.getcwd(), 'database', 'data0.geo'), dataset = 'EQ36', \n", " Dielec_method = 'DEW')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write EQ3/6 using default sourced GWB database and direct-access database" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new EQ3/6 database is ready for download\n", "Database for EQ36 generated successfully using JN91 dielectric constant, thermo.2021.dat gwb source database, solid solution included with cpx excluded, clay thermodynamics included\n", "CPU times: user 11.3 s, sys: 172 ms, total: 11.5 s\n", "Wall time: 11.8 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = np.array([0.010, 25, 60, 100, 150, 200, 250, 300]), P = 250, dataset = 'EQ36', \n", " sourcedb = 'thermo.2021', sourceformat = 'gwb', solid_solution = True, clay_thermo = True, \n", " print_msg = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Generate PHREEQC thermodynamic database" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write PHREEQC using default sourced database" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new PHREEQC database is ready for download\n", "CPU times: user 1.9 s, sys: 17.7 ms, total: 1.92 s\n", "Wall time: 1.92 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = T, P = P, cpx_Ca = 1, solid_solution = 'Yes', clay_thermo = 'Yes', dataset = 'PHREEQC')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write PHREEQC user-specified sourced database" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/noahawolayo/Documents/pygcc_build/src/pygcc/pygcc_utils.py:4891: UserWarning: Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied\n", " warnings.warn('Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied')\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Success, your new PHREEQC database is ready for download\n", "CPU times: user 1min 9s, sys: 725 ms, total: 1min 10s\n", "Wall time: 1min 13s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 600], P = 350, cpx_Ca = 0.5, solid_solution = True, clay_thermo = True, sourceformat = 'PHREEQC', \n", " sourcedb = os.path.join(os.getcwd(), 'database', 'PHREEQC_ThermoddemV1.10_15Dec2020.dat'), dataset = 'PHREEQC')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write PHREEQC using user-specified sourced Pitzer database" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new PHREEQC database is ready for download\n", "CPU times: user 693 ms, sys: 16.7 ms, total: 710 ms\n", "Wall time: 811 ms\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 350], P = 200, sourcedb = os.path.join(os.getcwd(), 'database', 'pitzer.dat'), sourceformat = 'PHREEQC', dataset = 'PHREEQC')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write PHREEQC user-specified sourced database using FGL97 dielectric constant" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new PHREEQC database is ready for download\n", "CPU times: user 1min 4s, sys: 417 ms, total: 1min 5s\n", "Wall time: 1min 6s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 600], P = 300, cpx_Ca = 0.1, sourcedb = os.path.join(os.getcwd(), 'database', 'PHREEQC_ThermoddemV1.10_15Dec2020.dat'), dataset = 'PHREEQC',\n", " solid_solution = True, Dielec_method = 'FGL97', clay_thermo = True, sourceformat = 'PHREEQC')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write PHREEQC user-specified sourced database using DEW model" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new PHREEQC database is ready for download\n", "CPU times: user 3.3 s, sys: 50 ms, total: 3.35 s\n", "Wall time: 3.64 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "T = np.array([50, 100, 150, 300, 450, 500, 600, 700])\n", "write_database(T = T, P = 1500, sourcedb = os.path.join(os.getcwd(), 'database', 'PHREEQC_ThermoddemV1.10_15Dec2020.dat'), \n", " sourceformat = 'PHREEQC', dataset = 'PHREEQC', Dielec_method = 'DEW')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Generate ToughReact thermodynamic database" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write ToughReact using user-specified EQ3/6 database" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/noahawolayo/Documents/pygcc_build/src/pygcc/pygcc_utils.py:5750: UserWarning: Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied\n", " warnings.warn('Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied')\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Success, your new ToughReact database is ready for download\n", "CPU times: user 43.4 s, sys: 437 ms, total: 43.8 s\n", "Wall time: 44.6 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = T, P = P, cpx_Ca = nCa, solid_solution = 'Yes', sourcedb = os.path.join(os.getcwd(), 'database', 'data0.dat'),\n", " dataset = 'ToughReact', sourceformat = 'EQ36')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write ToughReact using user-specified GWB database" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new ToughReact database is ready for download\n", "CPU times: user 1.62 s, sys: 20.2 ms, total: 1.64 s\n", "Wall time: 1.67 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 350], P = 250, cpx_Ca = 0.25, clay_thermo = 'Yes', dataset = 'ToughReact',\n", " sourcedb = os.path.join(os.getcwd(), 'database', 'thermo.com.tdat'), sourceformat = 'GWB')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write ToughReact using EQ3/6 user-specified Ptizer sourced database and JN91 dielectric constant" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new ToughReact database is ready for download\n", "CPU times: user 865 ms, sys: 5.8 ms, total: 871 ms\n", "Wall time: 876 ms\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 300], P = 200, sourceformat = 'EQ36', sourcedb = os.path.join(os.getcwd(), 'database', 'data0.fmt'),\n", " dataset = 'ToughReact', Dielec_method = 'JN91')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Generate Pflotran thermodynamic database" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " write Pflotran using user-specified EQ3/6 database" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/noahawolayo/Documents/pygcc_build/src/pygcc/pygcc_utils.py:5313: UserWarning: Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied\n", " warnings.warn('Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied')\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Success, your new Pflotran database is ready for download\n", "CPU times: user 37.4 s, sys: 455 ms, total: 37.8 s\n", "Wall time: 38.4 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = T, P = P, clay_thermo = 'Yes', sourcedb = os.path.join(os.getcwd(), 'database', 'data0.dat'),\n", " dataset = 'Pflotran', sourceformat = 'EQ36')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "write Pflotran using user-specified GWB database" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Success, your new Pflotran database is ready for download\n", "CPU times: user 2.64 s, sys: 20 ms, total: 2.66 s\n", "Wall time: 2.69 s\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "write_database(T = [0, 350], P = 250, cpx_Ca = 0.1, solid_solution = True, clay_thermo = True,\n", " sourcedb = os.path.join(os.getcwd(), 'database', 'thermo.com.tdat'), dataset = 'Pflotran', sourceformat = 'GWB')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate clay mineral thermodynamics" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "#%% specify the direct access thermodynamic database\n", "db_dic = db_reader(dbaccess = os.path.join(os.getcwd(), 'database', 'speq21.dat'))\n", "db_dic = db_dic.dbaccessdic\n", "\n", "folder_to_save = 'output'\n", "if os.path.exists(os.path.join(os.getcwd(), folder_to_save)) == False:\n", " os.makedirs(os.path.join(os.getcwd(), folder_to_save)) \n", "fid = open('./output/logK_test.txt', 'w')\n", "\n", "T = np.array([ 0.01, 25, 60, 100, 150, 200, 250, 300])\n", "TK_range = T + 273.15\n", "\n", "logKRxn = calcRxnlogK(T = T, P = 'T', Specie = 'Clay', dbaccessdic = db_dic,\n", " elem = ['Clinochlore', '3', '2', '0', '0', '5', '0', '0', '0', '0', '0'],\n", " densityextrap = True)\n", "logK, Rxn = logKRxn.logK, logKRxn.Rxn\n", "\n", "# output in EQ36 format\n", "outputfmt(fid, logK, Rxn, dataset = 'EQ36')\n", "# output in GWB format\n", "outputfmt(fid, logK, Rxn, dataset = 'GWB')\n", "# output in PHREEQC format\n", "outputfmt(fid, logK, Rxn, *TK_range, dataset = 'PHREEQC')\n", "# output in Pflotran format\n", "outputfmt(fid, logK, Rxn, dataset = 'Pflotran')\n", "# output in ToughReact format\n", "outputfmt(fid, logK, Rxn, dataset = 'ToughReact')\n", "fid.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate plagioclase solid-solution thermodynamics" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [], "source": [ "#%% specify the direct access thermodynamic database\n", "db_dic = db_reader(dbaccess = os.path.join(os.getcwd(), 'database', 'speq21.dat'))\n", "db_dic = db_dic.dbaccessdic\n", "\n", "folder_to_save = 'output'\n", "if os.path.exists(os.path.join(os.getcwd(), folder_to_save)) == False:\n", " os.makedirs(os.path.join(os.getcwd(), folder_to_save)) \n", "fid = open('./output/logK_test.txt', 'w')\n", "\n", "T = np.array([ 0.01, 25, 60, 100, 150, 200, 250, 300])\n", "logKRxn = calcRxnlogK(T = T, dbaccessdic = db_dic, P = 'T', X = 0.634, Specie = 'Plagioclase')#densityextrap = True)\n", "logK, Rxn = logKRxn.logK, logKRxn.Rxn\n", "\n", "logKRxnph = calcRxnlogK(T = T, P = 'T', X = 0.634, Specie = 'Plagioclase', \n", " dbaccessdic = db_dic, densityextrap = True, Al_Si = 'Arnórsson_Stefánsson')\n", "\n", "# output in EQ36 format\n", "outputfmt(fid, logK, Rxn, dataset = 'EQ36')\n", "# output in GWB format\n", "outputfmt(fid, logK, Rxn, dataset = 'GWB')\n", "# output in PHREEQC format\n", "outputfmt(fid, logKRxnph.logK, logKRxnph.Rxn, *TK_range, dataset = 'PHREEQC')\n", "# output in Pflotran format\n", "outputfmt(fid, logK, Rxn, dataset = 'Pflotran')\n", "# output in ToughReact format\n", "outputfmt(fid, logK, Rxn, dataset = 'ToughReact')\n", "fid.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate new reaction equilibrium constants and export to specific formats" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "#%% specify the direct access thermodynamic database\n", "db_dic = db_reader(dbaccess = os.path.join(os.getcwd(), 'database', 'speq21.dat'), \n", " sourcedb = os.path.join(os.getcwd(), 'database', 'thermo_latest.tdat'), sourceformat = 'GWB')\n", "\n", "Temp = np.array([25, 80, 140, 190, 240, 300, 350, 400])\n", "TK_range = Temp + 273.15\n", "\n", "folder_to_save = 'output'\n", "if os.path.exists(os.path.join(os.getcwd(), folder_to_save)) == False:\n", " os.makedirs(os.path.join(os.getcwd(), folder_to_save)) \n", "fid = open('./output/logK_test.txt', 'w')\n", "\n", "db_dic.sourcedic['HWO4-'] = ['HWO4-', 2, '1.000', 'H+', '1.000', 'WO4--']\n", "\n", "logK = calcRxnlogK(T = Temp, P = 250, Specie = 'HWO4-', dbaccessdic = db_dic.dbaccessdic, densityextrap = True,\n", " sourcedic = db_dic.sourcedic, specielist = db_dic.specielist, Specie_class = 'aqueous').logK\n", "\n", "Rxn = {}\n", "Rxn['type'] = 'tungsten'\n", "Rxn['name'] = 'HWO4-'\n", "Rxn['formula'] = ps.dbaccessdic['HWO4-'][0]\n", "Rxn['MW'] = calc_elem_count_molewt(ps.dbaccessdic['HWO4-'][0])[1]\n", "Rxn['min'] = ps.dbaccessdic['HWO4-']\n", "Rxn['V'] = Rxn['min'][5]\n", "Rxn['source'] = Rxn['min'][1]\n", "Rxn['spec'] = ['H+', 'WO4--']\n", "Rxn['coeff'] = [1.000, 1.000]\n", "Rxn['nSpec'] = len(Rxn['coeff'])\n", "d = calc_elem_count_molewt(ps.dbaccessdic['HWO4-'][0])[0]\n", "Rxn['elements'] = [item for sublist in [['%s' % v,k] for k,v in d.items()] for item in sublist]\n", "\n", "# output in EQ36 format\n", "outputfmt(fid, logK, Rxn, dataset = 'EQ36')\n", "# output in GWB format\n", "outputfmt(fid, logK, Rxn, *TK_range, logK_form = 'polycoeffs', dataset = 'GWB')\n", "# output in PHREEQC format\n", "outputfmt(fid, logK, Rxn, *TK_range, dataset = 'PHREEQC')\n", "# output in Pflotran format\n", "outputfmt(fid, logK, Rxn, dataset = 'Pflotran')\n", "# output in ToughReact format\n", "outputfmt(fid, logK, Rxn, dataset = 'ToughReact')\n", "fid.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate CO$_2$ activity and molality" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "T = np.array([ 0.010, 25 , 60, 100, 150, 175, 200, 250])\n", "P = 250*np.ones(np.size(T))\n", "\n", "TK = convert_temperature(T, Out_Unit = 'K')\n", "\n", "#%% Calculate CO2 activity and molality at ionic strength of 0.5M\n", "# with Duan_Sun\n", "log10_co2_gamma, mco2 = Henry_duan_sun(TK, P, 0.5)\n", "co2_activity = 10**log10_co2_gamma\n", "\n", "# with Drummond\n", "log10_co2_gamma = drummondgamma(TK, 0.5)\n", "co2_activityD = 10**log10_co2_gamma" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fluid Temperature [C]: 0.01\n", "Fluid Pressure [bar]: 250.0\n", "Molality of CO2 in aqueous phase: 1.9318448998874207\n", "Activity of CO2 in aqueous phase (Duan_Sun): 1.1291811186477334\n", "Activity of CO2 in aqueous phase (Drummond): 1.1340282640704162\n", "\n", "\n", "Fluid Temperature [C]: 25.0\n", "Fluid Pressure [bar]: 250.0\n", "Molality of CO2 in aqueous phase: 1.4446373503231484\n", "Activity of CO2 in aqueous phase (Duan_Sun): 1.117606348338049\n", "Activity of CO2 in aqueous phase (Drummond): 1.1228777554452154\n", "\n", "\n", "Fluid Temperature [C]: 60.0\n", "Fluid Pressure [bar]: 250.0\n", "Molality of CO2 in aqueous phase: 1.1670235810670366\n", "Activity of CO2 in aqueous phase (Duan_Sun): 1.1097025985737983\n", "Activity of CO2 in aqueous phase (Drummond): 1.1184645553679928\n", "\n", "\n", "Fluid Temperature [C]: 100.0\n", "Fluid Pressure [bar]: 250.0\n", "Molality of CO2 in aqueous phase: 1.0886091317050277\n", "Activity of CO2 in aqueous phase (Duan_Sun): 1.109458812741491\n", "Activity of CO2 in aqueous phase (Drummond): 1.1250331593722136\n", "\n", "\n", "Fluid Temperature [C]: 150.0\n", "Fluid Pressure [bar]: 250.0\n", "Molality of CO2 in aqueous phase: 1.1736136778845458\n", "Activity of CO2 in aqueous phase (Duan_Sun): 1.1191485747389718\n", "Activity of CO2 in aqueous phase (Drummond): 1.1457708279040804\n", "\n", "\n", "Fluid Temperature [C]: 175.0\n", "Fluid Pressure [bar]: 250.0\n", "Molality of CO2 in aqueous phase: 1.2783771607594785\n", "Activity of CO2 in aqueous phase (Duan_Sun): 1.1276219378379773\n", "Activity of CO2 in aqueous phase (Drummond): 1.1602093887224851\n", "\n", "\n", "Fluid Temperature [C]: 200.0\n", "Fluid Pressure [bar]: 250.0\n", "Molality of CO2 in aqueous phase: 1.4208468773780791\n", "Activity of CO2 in aqueous phase (Duan_Sun): 1.1385618727912967\n", "Activity of CO2 in aqueous phase (Drummond): 1.1769259206651324\n", "\n", "\n", "Fluid Temperature [C]: 250.0\n", "Fluid Pressure [bar]: 250.0\n", "Molality of CO2 in aqueous phase: 1.77261489074432\n", "Activity of CO2 in aqueous phase (Duan_Sun): 1.1700221965116315\n", "Activity of CO2 in aqueous phase (Drummond): 1.2163347543456382\n", "\n", "\n" ] } ], "source": [ "for i in range(len(T)):\n", " print('Fluid Temperature [C]: ', T[i])\n", " print('Fluid Pressure [bar]: ', P[i])\n", " print('Molality of CO2 in aqueous phase: ', mco2.ravel()[i])\n", " print('Activity of CO2 in aqueous phase (Duan_Sun): ', co2_activity.ravel()[i])\n", " print('Activity of CO2 in aqueous phase (Drummond): ', co2_activityD.ravel()[i])\n", " print('\\n')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate water activity" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fluid Temperature [C]: 0.01\n", "Fluid Pressure [bar]: 250.0\n", "Water activity: 0.9832827062135284\n", "Water osmotic coefficient: 0.9357947724582291\n", "NaCl mean activity coefficient: 0.7028482566450444\n", "\n", "\n", "Fluid Temperature [C]: 25.0\n", "Fluid Pressure [bar]: 250.0\n", "Water activity: 0.9834680138300679\n", "Water osmotic coefficient: 0.925334742394683\n", "NaCl mean activity coefficient: 0.6831227245014528\n", "\n", "\n", "Fluid Temperature [C]: 60.0\n", "Fluid Pressure [bar]: 250.0\n", "Water activity: 0.9834674900419101\n", "Water osmotic coefficient: 0.925364305805263\n", "NaCl mean activity coefficient: 0.6731415110785025\n", "\n", "\n", "Fluid Temperature [C]: 100.0\n", "Fluid Pressure [bar]: 250.0\n", "Water activity: 0.9836234739653248\n", "Water osmotic coefficient: 0.9165610286207091\n", "NaCl mean activity coefficient: 0.6470401472320056\n", "\n", "\n", "Fluid Temperature [C]: 150.0\n", "Fluid Pressure [bar]: 250.0\n", "Water activity: 0.9839617696261588\n", "Water osmotic coefficient: 0.8974734053857932\n", "NaCl mean activity coefficient: 0.601472728601265\n", "\n", "\n", "Fluid Temperature [C]: 175.0\n", "Fluid Pressure [bar]: 250.0\n", "Water activity: 0.9841776700481395\n", "Water osmotic coefficient: 0.8852951070804821\n", "NaCl mean activity coefficient: 0.5748272882994092\n", "\n", "\n", "Fluid Temperature [C]: 200.0\n", "Fluid Pressure [bar]: 250.0\n", "Water activity: 0.9844315601276349\n", "Water osmotic coefficient: 0.870977343196605\n", "NaCl mean activity coefficient: 0.5454072718625895\n", "\n", "\n", "Fluid Temperature [C]: 250.0\n", "Fluid Pressure [bar]: 250.0\n", "Water activity: 0.9850902029390286\n", "Water osmotic coefficient: 0.8338513425870471\n", "NaCl mean activity coefficient: 0.4767054502731833\n", "\n", "\n" ] } ], "source": [ "#%% Calculate Water activity, osmotic coefficient and NaCl mean activity coefficient at an ionic strength of 0.5M\n", "aw, phi, mean_act = Helgeson_activity(T, P, 0.5, Dielec_method = 'JN91')\n", "for i in range(len(T)):\n", " print('Fluid Temperature [C]: ', T[i])\n", " print('Fluid Pressure [bar]: ', P[i])\n", " print('Water activity: ', aw.ravel()[i])\n", " print('Water osmotic coefficient: ', phi.ravel()[i])\n", " print('NaCl mean activity coefficient: ', mean_act.ravel()[i])\n", " print('\\n')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Calculate salt water properties" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Pressure of vapor + liquid + halite coexistence [bar]: 11.151719135466353\n", "Composition of halite-saturated liquid (halite liquidus) [-]: 0.12609895126323659\n", "Composition of halite-saturated vapor [-]: nan\n", "Composition of liquid at vapor + liquid coexistence [-]: [-0.30583232]\n", "Composition of vapor at vapor+liquid coexistence [-]: [nan]\n", "Liquid NaCl density [kg/m3]: 0.918850325064366\n", "molar volume [mol/m3]: 20.010808061663727\n", "Specific enthalpy of an H2O–NaCl solution [J/kg]: 856.75482704851\n", "viscosity [Pa-s]: 0.00014914486904138454\n", "Isobaric heat capacity [J/kg/K]: 4161.461344041091\n" ] } ], "source": [ "from pygcc.pygcc_utils import Driesner_NaCl, concentration_converter\n", "\n", "water_salt = Driesner_NaCl(T = 200, P = 500, xNaCl = concentration_converter(val = 0.5, In_Unit='x_mol_kgsol', Out_Unit='x_mol'))\n", "\n", "print('Pressure of vapor + liquid + halite coexistence [bar]: ', water_salt.PVLH)\n", "print('Composition of halite-saturated liquid (halite liquidus) [-]: ', water_salt.xL_NaCl)\n", "print('Composition of halite-saturated vapor [-]: ', water_salt.xV_NaCl)\n", "print('Composition of liquid at vapor + liquid coexistence [-]: ', water_salt.xVL_Liq)\n", "print('Composition of vapor at vapor+liquid coexistence [-]: ', water_salt.xVL_Vap)\n", "print('Liquid NaCl density [kg/m3]: ', water_salt.rho)\n", "print('molar volume [mol/m3]: ', water_salt.vm)\n", "print('Specific enthalpy of an H2O–NaCl solution [J/kg]: ', water_salt.H)\n", "print('viscosity [Pa-s]: ', water_salt.mu)\n", "print('Isobaric heat capacity [J/kg/K]: ', water_salt.Cp)" ] } ], "metadata": { "kernelspec": { "display_name": "py313", "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.13.7" } }, "nbformat": 4, "nbformat_minor": 4 }