From 5e9e3c9d2316bbe5976f328d0c99b1c01c4b30a0 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Fri, 27 Mar 2020 15:41:43 +0000 Subject: [PATCH] #729 add experimental data --- examples/notebooks/compare-ecker-data.ipynb | 180 +++++++++++++++++- pybamm/expression_tree/symbol.py | 2 +- pybamm/input/discharge_data/Ecker_1C.csv | 62 +++--- pybamm/input/discharge_data/Ecker_5C.csv | 66 +++---- .../cells/kokam_Ecker2015/parameters.csv | 2 +- .../electrolyte_diffusivity_Ecker2015.py | 31 +-- pybamm/parameters_cli.py | 3 +- pybamm/simulation.py | 2 + 8 files changed, 256 insertions(+), 92 deletions(-) diff --git a/examples/notebooks/compare-ecker-data.ipynb b/examples/notebooks/compare-ecker-data.ipynb index 09f57a6039..23fdc414d3 100644 --- a/examples/notebooks/compare-ecker-data.ipynb +++ b/examples/notebooks/compare-ecker-data.ipynb @@ -6,7 +6,7 @@ "source": [ "# Comparing with Experimental Data\n", "\n", - "In this notebook we show how to compare results generated in PyBaMM with data. We compare the results of the DFN model (see the [DFN notebook](./models/DFN.ipynb)) with the results from Ecker et. al. [1]. Results are compared for a constant current discharge at 1C and at 5C." + "In this notebook we show how to compare results generated in PyBaMM with experimental data. We compare the results of the DFN model (see the [DFN notebook](./models/DFN.ipynb)) with the experimental data from Ecker et. al. [1]. Results are compared for a constant current discharge at 1C and at 5C." ] }, { @@ -23,18 +23,192 @@ "outputs": [], "source": [ "import pybamm\n", - "import numpy as np\n", "import os\n", + "import pandas as pd\n", + "import numpy as np\n", "import matplotlib.pyplot as plt\n", "os.chdir(pybamm.__path__[0]+'/..')" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We then load the Ecker data in from the `.csv` files using `pandas`" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "voltage_data_1C = pd.read_csv(\"pybamm/input/discharge_data/Ecker_1C.csv\", header=None).to_numpy()\n", + "voltage_data_5C = pd.read_csv(\"pybamm/input/discharge_data/Ecker_5C.csv\", header=None).to_numpy()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the data is Time [s] vs Voltage [V]." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We load the DFN model and select the parameter set from the Ecker paper [1]. We update the C-rate an `InputParameter` so that we can re-run the same model at different C-rates without the need to rebuild the model. This is done by passing the flag `[input]`." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# choose DFN\n", + "model = pybamm.lithium_ion.DFN()\n", + "\n", + "# pick parameters, keeping C-rate as an input to be changed for each solve\n", + "chemistry = pybamm.parameter_sets.Ecker2015\n", + "parameter_values = pybamm.ParameterValues(chemistry=chemistry)\n", + "parameter_values.update({\"C-rate\": \"[input]\"})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For this comparison we choose a fine mesh of 1 finite volume per micron in the electrodes and separator and 1 finite volume per 0.1 micron in the particles" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "var = pybamm.standard_spatial_vars\n", + "var_pts = {\n", + " var.x_n: int(parameter_values.evaluate(pybamm.geometric_parameters.L_n / 1e-6)),\n", + " var.x_s: int(parameter_values.evaluate(pybamm.geometric_parameters.L_s / 1e-6)),\n", + " var.x_p: int(parameter_values.evaluate(pybamm.geometric_parameters.L_p / 1e-6)),\n", + " var.r_n: int(parameter_values.evaluate(pybamm.geometric_parameters.R_n / 1e-7)),\n", + " var.r_p: int(parameter_values.evaluate(pybamm.geometric_parameters.R_p / 1e-7)),\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We create a simulation using our model, parameters and number of grid points" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "sim = pybamm.Simulation(model, parameter_values=parameter_values, var_pts=var_pts)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can then solve the model for a 1C and 5C discharge " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "C_rates = [1, 5] # C-rates to solve for\n", + "t_evals = [\n", + " np.linspace(0, 3800, 100), \n", + " np.linspace(0, 720, 100)\n", + "] # times to return the solution at\n", + "solutions = [None] * len(C_rates) # empty list that will hold solutions\n", + "\n", + "# loop over C-rates\n", + "for i, C_rate in enumerate(C_rates):\n", + " sim.solve(t_eval=t_evals[i], solver=pybamm.CasadiSolver(mode=\"fast\"),inputs={\"C-rate\": C_rate})\n", + " solutions[i] = sim.solution" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally we plot the numerical solution against the experimental data" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAEYCAYAAABCw5uAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd5xU1f3/8dfZ2TLb+7KVXqUrIoKggqgxgjERosbEEjUajcbEtF8SBEJ6okZNURNL7DVfsSsIVgTpvdctsJXtfc/vjxnWZdlZdoed2fZ+Ph7z2Jk79849Q8w987nncz7HWGsRERERERER8bWAzm6AiIiIiIiI9A4KQEVERERERMQvFICKiIiIiIiIXygAFREREREREb9QACoiIiIiIiJ+oQBURERERERE/EIBqIiIiIiIiPiFAlCRHsIYc7sxZrUxptoY80Sz96KMMfcbYw4aY8qMMXvcrxM6qbkiIiJdjjFmuTGmyt1XlhljdjR5L8UY8x9jTI4xptQYs90Ys8AYE96ZbRbpbhSAivQc2cAi4LGmG40xwcBSYCRwMRAFnA0UABP93EYREZGu7nZrbYT7MQzAGBMHrABCgbOttZHATCAGGNR5TRXpfgI7uwEi0jGsta8CGGMmAOlN3voO0Bc431pb5t6WC/zGvy0UERHptn4ElALXWGsbAKy1h4A7O7VVIt2QRkBFer4LgHeaBJ8iIiLi2e+NMfnGmE+NMee5t10AvHos+BQR7ykAFen54oGczm6EiIhIN/AzYCCQBjwCvG6MGYT6UpEOowBUpOcrAFI6uxEiIiJdnbV2pbW21Fpbba19EvgUuAT1pSIdRgGoSM+3BLhIVfpERETazQIGV196uTFGv51FTpH+TyTSQxhjAo0xTsABOIwxTmNMIPAUcAh4xRgz3BgTYIyJN8b8P2PMJZ3aaBERkS7CGBNjjLnoWP9pjPkWMA14B7gXVxX5J40x/dz7pxlj7jXGjOnEZot0OwpARXqOXwGVwM+Ba9zPf2WtrcZVPGE78D5QAqwCEoCVndNUERGRLicI13JmeUA+8APga9bandbaQmAyUAusNMaU4lrirBjY3UntFemWjLW2s9sgIiIiIiIivYBGQEVERERERMQvFICKiIiIiIiIXygAFREREREREb9QACoiIiIiIiJ+EdjZDWivhIQE279//85uhoiI9HBr1qzJt9YmdnY7Opr6URER8QdP/Wi3C0D79+/P6tWrO7sZIiLSwxljDnR2G3xB/aiIiPiDp35UKbgiIiIiIiLiFwpARURERERExC8UgIqIiIiIiIhfdLs5oCIiPVltbS2ZmZlUVVV1dlN6DafTSXp6OkFBQZ3dFBEROUXqR/2vvf2oAlARkS4kMzOTyMhI+vfvjzGms5vT41lrKSgoIDMzkwEDBnR2c0RE5BSpH/Uvb/pRpeCKiHQhVVVVxMfHq9P0E2MM8fHxulMuItJDqB/1L2/6UQWgIiJdjDpN/9K/t4hIz6Lrun+199+71wagb23KIfPDJ+G+UTA/xvV344ud3SwREZEuz1rLMysPUFhe09lNERGRbqZXBqD1DZZVr/2LuA/uhuJDgHX9ff0OBaEi0us5HA7GjRvHyJEjGTt2LH/9619paGgAYPny5URHRzNu3DjGjRvHBRdcAMD8+fMJCwsjNze38XMiIiI6pf3ie/krnua8t6YT8+ckGu4dqb5TRKQJ9aOt65UBqCPAMC/sZcJMszu3tZXUv7+gcxolItJFhIaGsn79erZs2cL777/P22+/zYIFX14bp06dyvr161m/fj1Llixp3J6QkMBf//rXzmhyt2eMcRhj1hlj3mjhvRBjzAvGmN3GmJXGmP5N3vuFe/sOY8xFfmnsxhdJXPYT0kw+AVgCSjKxuoErItJI/WjremUAChBQktXidlOSxWUPfcLfluxic1Yx1lo/t0xEpOtISkrikUce4aGHHjrp9fCGG27ghRdeoLCw0E+t61HuBLZ5eO+7QJG1djBwH/BHAGPMacCVwEjgYuAfxhiHz1u6dCHUVh63ydRWYpcu9PmpRUS6G/WjJ+q9y7BEp7vTb49X5uyDI8Bw/9Kd3LdkJynRTqYPT+LqsJWctvU+THGW69gZ82DM3E5ouIj0Fgte38LW7JIO/czTUqO4Z9bIdh0zcOBA6uvrG9OCPv74Y8aNGwfAnDlz+OUvfwm4UoVuuOEG/va3vx13p1daZ4xJB74K/Bb4UQu7XAbMdz9/GXjIuCo+XAY8b62tBvYZY3YDE4EVPm1wcWb7touIdBL1o11T7w1AZ8xzzflsehc3KJSor/6GV8dMIb+smmXbc1my7Qg1655ngHkEcyxlt/gQdvEdGFAQKiK9ztSpU3njjRMyRQG44447GDduHHfffbefW9Wt3Q/8FIj08H4acAjAWltnjCkG4t3bP2+yX6Z72wmMMTcDNwP07dv31Frr4QZuVkM8GzflcMnolFP7fBGRHq6396O9NwA9FjguXei6a9tsVDMhIoQ5EzKYMyGDhnuvJaDk+Pmipq6S/Nd+ycf1k7lgRB8inUH+/gYi0sO19w6rr+zduxeHw0FSUhLbtnnKEnWJiYnh6quv5u9//7ufWte9GWMuBXKttWuMMef56jzW2keARwAmTJhwanNLWriBazFsCD2LX7yykTHp0aTHhp3SKUREOoL60a6p9wag4Ao22zCC6Wm+aFx9Hne9sIHgwABmntaHb53Vl7MHauFbEek58vLyuOWWW7j99tvbfG370Y9+xJlnnkldXZ2PW9cjTAFmG2MuAZxAlDHmaWvtNU32yQIygExjTCAQDRQ02X5Munubb42ZCwc/h9WPAa5Y1mD5Sv0HfGgHcufzkbxw8yQCHb22zISISCP1oydS79AW0ektbjbR6bxy62SuntiXT3blc/WjK5nx1w/54MWHXGXptb6oiHRDlZWVjeXjL7jgAi688ELuueeeNh+fkJDA5ZdfTnV1tQ9b2TNYa39hrU231vbHVVDog2bBJ8Bi4Fr38yvc+1j39ivdVXIHAEOAVX5p+K73OBZ8HhNQV8n8sJdZc6CIB5bu8kszRES6IvWjrTPdrcrrhAkT7OrVq/170o0vtjhflFkPNI6gVtXW8+bGHA5++ATfK/7b8Uu8NNtXRMSTbdu2MWLEiM5uRq/T0r+7MWaNtXaCv9rgTsG921p7qTFmIbDaWrvYGOMEngLGA4XAldbave5jfgncANQBP7TWvn2y83RIPzo/huYBqPtbcPfID3llbSbP3jiJswfFn9p5RETaSf1o52hPP+rzEVBv1zbrUsbMdQWQ0RmAcf1tFlA6gxx844x07jLPt7i+aNW792hJFxER8chau9xae6n7+Txr7WL38ypr7Rxr7WBr7cRjwaf7vd9aawdZa4e1JfjsMB4yg4hOZ8HskfSPD+eO59exJbvYb00SEZHuwR8puO1e26xLGjMX7toM84+6/noazfRQhj64LIcL7v2Qxz/dR3FlrefzbHzRlbar9F0REemqZsxzZfc0FRQKM+YRHhLIw98+g8AAwxX/XME7mw93ThtFRKRL8mkA2mRts3972OUy4En385eBGaa7V/DxcFe4MiyZCGcQC17fypm/XcLN/13Na+uzKKtuMrn4WKpv8SHAuv6+foeCUBER6VpOkhk0tE8kr902haHJkdzy9Br+vmy3soBERATwfRVcb9c2y2+6U4euX+ZrHtYXDf/KQl4bM4VNmcW8ui6Ttzbl8N7WI4QEBjB1SALThiZy1WfzCWp6HLg+Z+lCzR8VEZGupXkl+WMZPO6lzZJmzOOFm7/Bz17ZyJ/f3cG6g0f5yUXDGJbs6SeBiIj0Bj4LQDtybbMOXb/M106yvujo9GhGp0fz66+exuoDRby5MZtlO/JYsi2Xa0KyoKXxXw9pvSIiIl1C82J97gweJ3D/N+cwMjWKB5bu5qL7P+Iro5K5Y8YQRqREdWqTRUSkc/hyBPRU1jbr3tqwvmhAgGHigDgmDohjAXCgoJyKR1OIqMo5Yd+CwET+9/FeTu8Xy8jUKEICHcfvsPFFjwGviIiIzy1deHzmDzRm8Jgxc7l52iDmTsjgP5/s44lP9/P25sOcOzSRayb1Y/rwJBwB3Xv2jYiItJ3PAlBr7S+AX8BxpeU9rW22guPXNut1+sWHwyULT0jfrTYhPMjVPPGmq45TsCOA01KjGJcRw7iMGCZXfEDi8p9gmt11BhSEiohXHA4Ho0ePbnx95ZVX8vOf/9xn51u8eDFbt2716TmWL19OcHAwkydP9tk5ejVPmTpNtseEBfPjC4dx4zkDeeKz/Ty76gA3/Xc1qdFOrpzYlzkT0kmJDm35c0REuhH1o63z9RzQEzRd2wz4D/CUMWY37rXN/N2eLqWF9N2QGfOYP2Yut5ZUse5gEWsPHmX9oaO88MUhnvhsP58Ez8cEnHjX2brvOotID+eDDIjQ0FDWr1/fQQ1sXV1dHbNnz2b27Nk+Pc/y5cuJiIhQAOor0enuAnotbG++KSyIOy8YwvfPH8TSbUd4ZuVB7n1/J/cv2ck5QxKZc0Y6M0/rgzPIceLniYh0NPWjbdKR/ajpbgOOHbKAdg9QV9/Artwyhj/cF9PCYuANGK7NeJdRadGMSo1mdFo0GXGhHFdkWKm7Il1OuxbQbj7vDlxLYTRbp7i9IiIiKCsrO25bcXExEydOZPHixQwbNoyrrrqK6dOnc9NNNxEREcFNN93Ee++9R3JyMs8//zyJiYns2bOH2267jby8PMLCwnj00UcZPnw41113HU6nk3Xr1jFlyhTGjBnD6tWreeihh7juuusIDQ1l3bp15Obm8thjj/Hf//6XFStWcNZZZ/HEE08A8N5773HPPfdQXV3NoEGDePzxx4mIiKB///5ce+21vP7669TW1vLSSy/hdDqZNGkSDoeDxMREHnzwQaZOnXrSf3dPC2h3dz7pR0/xv8UDBeW8siaTl9dkkl1cRZQzkEvHpvKN09M4vW8s3b1Avoj4j/rRrt+P+mMdUPGBQEcAI1KiMB6WfSkOSqKgrIZHP9rLbc+uZdqflzF2wXtc+cgKFr2xlS8W/4uGxVryRaRba2Xe3amorKxk3LhxjY8XXniB6Ojoxo7t+eefp6ioiJtuugmA8vJyJkyYwJYtWzj33HNZsGABADfffDMPPvgga9as4S9/+Qvf//73G8+RmZnJZ599xr333nvC+YuKilixYgX33Xcfs2fP5q677mLLli1s2rSJ9evXk5+fz6JFi1iyZAlr165lwoQJx31OQkICa9eu5dZbb+Uvf/kL/fv355ZbbuGuu+5i/fr1J3Sa0gFaW5alDetb94sP50cXDuPjn03nqe9O5PzhSby6NpNv/HMF5/1lOfcv2cmBgnL/fy8R6dnUj3ZKP+r3FFzpYB6WfYmdtYi3xkyluq6enYfL2JxdzKasYrZkl/DU5we4LuDPBLSQulvz7nzqh3+D0GClPol0eW2Yd+cNT6lDM2fO5KWXXuK2225jw4YNjdsDAgL45je/CcA111zD17/+dcrKyvjss8+YM2dO437V1dWNz+fMmYPD0fJ1ZtasWRhjGD16NH369GmcRzNy5Ej2799PZmYmW7duZcqUKQDU1NRw9tlnNx7/9a9/HYAzzjiDV1991dt/BmmvlgrweaiO27h/M44Aw9QhiUwdkkhZdR1vb8rh1bVZ/G3pLu5fsovxfWO4fHwal45JJS482MdfSER6PPWjgP/7UQWg3d1Jln0JCXQ0Lv1ylfuQuvoGHL9pudhwYFk2w+95h4GJEZyWEsVpqVGNfxMiQr7cUem7Ip2vHfPuOkJDQwPbtm0jLCyMoqIi0tNbPo8xhoaGBmJiYjzOgQkPD/d4npAQ17UmICCg8fmx13V1dTgcDmbOnMlzzz3X6vEOh4O6uro2fTfxkdZGF07SZ0SEBDJnQgZzJmSQfbSSxRuy+d/aLOa9toWFr29l6pAEvjY+jZmn9SEsWD9nRMQL6kdbPd5X/ahScHuCMXPhrs0w/6jr70k69UBHgMfU3erwFG6fPoT+8WGsOVDEH97eznceW8WERUs487dL+M5jq3jtqfupe+0HSt8V6Wwz5rnmqjQVFOra7gP33XcfI0aM4Nlnn+X666+ntrYWcHWoL7/8MgDPPvss55xzDlFRUQwYMICXXnoJAGvtcXd7T8WkSZP49NNP2b17N+BKXdq5c2erx0RGRlJaWtoh55d26KDRhdSYUG45dxDv3jWNt++cynenDmD74VLufH49Z/xmCXc8t46l245QU9fQAY0WkV5D/Sjg/35Utwx7Kw+pu6EXL+BHY4Y2bjpaUcPWnBK2ZpewLaeUrTklTCh8kEBTdfzn1VZS9tY8NoXPYERKJDFhSo0S8bmTZEB469jclWMuvvhirr/+ev7973+zatUqIiMjmTZtGosWLWLBggWEh4ezatUqFi1aRFJSEi+88AIAzzzzDLfeeiuLFi2itraWK6+8krFjx55S2wASExN54oknuOqqqxrTkRYtWsTQoUM9HjNr1iyuuOIKXnvttRaLJ4iP+GB0YURKFCNSovjZRcP5Yn8h/7c+m7c357B4QzYxYUF8ZVQys8amctaAeK0vKiKtUz/aKf2oquD2Zl6m0dr5MS1X3rWGgdXPAJAS7WR4ciTDU6Jcf5OjGHT4LQKX/UZpuyKtaFf1vi6ipWp/3Y2q4PqIjypMNldT18DHu/JYvCGb97YcobK2nqTIEC4dk8qssSmMy4hRJV2RXkL9aOdoTz+qEdDerKWCEW1gPNzRtlFp/HfWRLbllLD9cCnbckr4ZHc+tfWW2QGf8IegfxNoalw7Fx+i/rUfcLS8hrhJ39IPAxGRnqg9owunUFsgODCAGSP6MGNEHypq6li6LZfFG7J5+vMDPPbpPjLiQpk1JpVZY1MZnhypPkdEpBMpAJX285C+65h5D9OGJjJtaGLj5pq6Bvbml9H3yR8TVllz3Mc46quofPseRr+TyJA+EQzrE8nQxkcEiZEh+pEg0g1097u24mNtudnZzmq5rQkLDmTWWFewWVxZy7tbDvPGxhwe/mgv/1i+h8FJEe5gNIWBiRFefCERkY7V2/pRBaDSfu24ox0cGMDw5CioPNziR6UFFPD18WnsPFLKu1sO8/wXX46sRocGMbRPBN8IWsGleY8SXnWY+sg0HDPvwSh1V3owa61uvvhRd5uK0iOdQrXc1kSHBjF3QgZzJ2RQUFbNW5sP88aGbO5fupP7luxkZGoUs8amcumYFNJjw07xS4hIV6F+1L/a248qABXvtDd910ParolOZ+FlowDXf7x5ZdXsPlLGziOl7Mwto8++17is5G+E4ho9DSzNpPKV2/jH+zvJzpjFoKRwBidGMCgpgr5xYQQ5VNhZujen00lBQQHx8fHqPP3AWktBQQFOp7Ozm9K7+WgtvqbiI0L49qR+fHtSPw4XV/HmJlfhoj+8vZ0/vL2dM/rFMmtMCl8dk0piZMjJP1BEuiT1o/7lTT+qIkTiH94WorhvVIuBa54jia8G/JPc0i8X5A0MMPSND2NgQjgDEyM4r3o5Z+x+gODyHIhOx6jokXQDtbW1ZGZmUlVVdfKdpUM4nU7S09MJCgo6bruKEPmRh2s90Rmu5cV86GBBBa9vzOb1DdlsP1xKgIEpgxOYPTaVi0YlE+UMOvmHiEiXoX7U/9rbjyoAFf/xpsDE/BhooeIuGJh/lJKqWvbmlbM3r4zduWXsyy9nb145IwvfZVHAI4SZL+edVhHCc8l3kz9gNv3iw+kXF0b/hHCSNNdURFqgANSP/FQt92R2HSll8YZsXlufzcHCCoIDA5g+LInLxqVy/vAknEEOv7VFRKS7UwAq3ZOXd8XtfaMwLRx32CQypfoB6hu+/O/eGRRA37gw+saF0y8+jHOrlnHm3odwVuRgo9IIuOAejZyK9EL+CECNMU7gIyAE17SYl6219zTb5z7gfPfLMCDJWhvjfq8e2OR+76C1dvbJztll+9G23qQ8hWq5bWWtZUNmMa+tz+L1DTnkl1UTGRLIhSOTmT0ulcmD4jXlQ0TkJBSASvfk7V3xVkZOa39dSPbRSvbll3OwsIIDBa7HwcJyRhe9x2/M8SOnlYTwSNSdHEy/lIy4UDJiw+gbH0ZGbBhJkSEEaKFzkR7JTwGoAcKttWXGmCDgE+BOa+3nHvb/ATDeWnuD+3WZtbZdpVy7dT/aCSOl9Q2WFXsKeG19Fu9sOUxpVR1x4cF8ZVQyl45JZeKAOBzqB0RETqB1QKV7as8ack15KHpEdDpBjgBXCm58+Alv2/tuwxQfv1xMKNVcXfYEs3ZP4khpFU3v2QQHBpAeE0p6XBgZsaGkx4aREef6mx4bSnx4sNJ7RcQj67oLfKz+fpD70dqd4auAe1p5v2fzUbXc1jgCDOcMSeCcIQksunwUH+7IY/GGbF5Zm8kzKw+SEBHCV0Ylc8noFAWjIiJtoABUur72VtwFj2uVMmNeq4cZDxUXExvy+Pz/zaC6rp6sokoOFlZwqKiSzKIKMgsrOVRUwabMo0ytWsaswBdJNflk2wTu5krWR88kLTaMtJhQ0mNDSYsJJSXaSWpMKMnRTqVxifRyxhgHsAYYDPzdWrvSw379gAHAB002O40xq4E64A/W2v/zcOzNwM0Affv27cDW+5k31XI7MGU3JNDBhSOTuXBkMuXVdSzbkcubG3N4ac0hnvr8AAkRwcw8rQ8XjUxm8qAEggN1fRcRaU4BqPRMPhg5BdePj4GJES0vXr7xRezixzF1rqA33eTz+4B/81RoFP9XPoXNWcUUlh8/umoMXBO6kh/wLIkN+RQHJbFq8A8oHXw5SVEh9IlykhgRQnRokFJ9RXooa209MM4YEwP8zxgzylrb0iT3K3HNEa1vsq2ftTbLGDMQ+MAYs8lau6eFczwCPAKuFFwffA3/OMk1+gTNU3aLD7lewymPmIaHBHLpmFQuHZPaGIy+u+UIi9dn89yqQ0SGBHLe8CRmntaH84YlqpquiIib5oCKNHUq84vaUDCpoqaO7KOVZB+tIqe4koid/2Pm7t8SbL9cTqbCBvPz2htZ3HBO47bAAEN8RDDx4SHERwQTFx7MtKplzMx5mMjqI1SGpnDo9LupPe0KokODiAoNIjIkUEGryCnojCq4xph5QIW19i8tvLcOuM1a+5mHY58A3rDWvtzaObp1P9rea3QnLO9SVVvPZ3vyeWfzYZZuy6WgvIYgh+GsAfHMGJHEBSP6kBEX5pNzi4h0JZoDKtIW3o6cQptSw8KCAxmcFMngpEjXhk8ehSbBJ0CYqeHe+Nf41uV3k1taTV5pNfllxx41FJTXMPjwW1xS8w9C3cWSwiqzyfjk5/x82e7GwDXAQERIIJHOICKdgUS5/55bvYxZ+f8mpjaX0pA+rBl8B7n9ZxEeEkh4SCARIYGEBTuIcL8ODw7EGRSguawiPmCMSQRqrbVHjTGhwEzgjy3sNxyIBVY02RaLK1itNsYkAFOAP/mn5Z2kvdfo9qbsdkC6rjPIwfThfZg+vA/1DZZ1B4t4f+sR3t92hAWvb2XB61sZkhTB9BFJTB+WxOn9YjUVQ0R6FY2AinQUb+60n2Sd0/aeqzIsldenv0dJZS0llbUUV9ZSWlVHSVUdpVW1TCh5nx+UP4ST1kdcm/ta4Kf8NPAFkikg1yTwTPi1rIme2RiwRoQEEuF0/Y1yup5HhrgC3khnEFGhgUSFBhERrFFZ6T78VAV3DPAk4AACgBettQuNMQuB1dbaxe795gNOa+3Pmxw7GXgYaHAfe7+19j8nO2ev6kfbc132Q4Xd/fnlLN2ey9JtR1i1r5C6BkukM5BpQxM5b2gi5w5LJCnS2SHnEhHpbFqGRcTXvPnx4m16WAcHrnWR6ez51ueUVddRUVNHeXUdZdX1VNTUkXLwdc7bsYighqrG/atNCP+MvIN3AqZRVl3nelTVUddwYptmB3zCT48VZiKBh8zVrAifTkxYMDGhQcSGBREbHkxcWDBxEcHEhweTEBHiekSGEB7s0OirdIrOSMH1h17Vj7bnuuzndN3Sqlo+3Z3PB9tzWbYjj7xS143B01KiOG9YIucOTdToqIh0a0rBFfE1b9J3vazW2+5CHMd4SDsLLM1iWHJky8d8/jA0CT4BQmw1PzTP88Mf/qpxm7WW6roGyqrrKHWPuIZse4VBnz9GYL3r+HTyWWge5rnISJYEnUtRRQ1788soKq+lrLruuHMcC1zDTAFHAhJ4Iep6didfQkq0s7GKcGq0q7JwTFiQglQROVF7rst+TteNdAZx8agULh6VgrWWrTklLN+Rx4c78nj4o738Y/keIkICmTwonmlDXQGp5o6KSE+gAFSkI7V3yRhv55z6M3Bt448yYwzOIAfOIAcJESGujS/dB/XHB6/BtpprK//Ltbf89Ljt1XX1FJXXUlBejdn0MkNXfRm4Jts8vlf8N/5UWcsTFWdRU9dw3LGRIYGkx4XRNy6U/gnhDEwIZ0BCBAMSwkmI0FqsIr1aW6/L7bk+dnB1XWMMI1OjGZkazW3nD6akqpbPdhfw4c48PtqZx3tbjwDQPz6MqUMSOWdIAmcPildlXRHplpSCK9JdeXP33Z9pwtDhqcJEZ2B/uInC8hqyj1aRddS9Fqt7bdaDhRUcLKigpv7LADUmLIihSZEM7hPBsD6RnJYaxfDkSCL1w01OQim4vUwXTde11rInr5xPduXx8a58VuwtoKKmHkeAYWx6NOcMSWTqkATGZcQoXVdEuhSl4Ir0NO0dbT12DPgnTRg6PFWY4kyMMcRHhBAfEcLo9OgTdmnY8CINSxbgKM2i3JnMG4k38Urt2by5MYdnKw827tc3LoxRaVGMTothbHo0o9KjNZog0pv5Ml0XvE7ZNcYwOCmCwUkRXDdlADV1Daw9WMQnu/L5eHc+D32wiweW7iI82MHEAXFMGZzA5EEJDE+OVNE3EemSNAIqIifn7Vwnb6tKeju60Mr57Og5HCmpZltOCVtzStiSXcymrGIOFX6575CkCM7oF8vp/WI5o18sAxPClb7bi2kEVDxq7zXKhxV2iytqWbE3n093F/Dp7nz25pcDEB8ezKRB8UwZ5ErX7R8fpuuZiOzuBxQAACAASURBVPiVquCKSOfwV6oweBW4FpbXsCmrmA2HjrLuYBFrDhRRUuUqiJQQEczEAXGcNSCeSQPjGdonQj/gehEFoOJRe69Rp5Ky285raNbRSj7bnc+KPQV8uiefIyWu6rop0U4mDYzn7EHxnD0wnvTYUF3PRMSnlIIrIp3DX6nC4FVaXFx4MOcOTeTcqmWwcSGWTOoSUvli0A94qeZsVu4t4K1NhwFXQHr2oASmDIpnyuAEVaQU6a3ae43yJmUXvCp2lBYTypwJGcyZkIG1lr355azYU8CKPQV8tDOP/63LatzvrIFxTBoYz1kD4ugbpxFSEfEPjYCKSM/hg9RdxszlUGEFK/YW8NnufD7dU9C4Xt/AhHCmDklg6pBEzh4UT3iI7un1JBoBlQ7j7bWpg4sdNTRYduWWsXJfAZ/vLeDzvYUUltcAkBzlZOKAOM4cEMfE/nEMSYrQHFIROSUaARWRns/bgklLFx5/DLheL10IY+aSERdGRlwYc90jCrtzy/h4Vz4f78rjxdWZPLniAEEOw5n949wLyCcpXVdEvuTttcnbkVMPAgIMw5IjGZYcyXfO7t94PVu5r5CV+wr5fG8BizdkAxAdGsSEfrGc0T+WCf3iGJMejTPI4dV5RUSaUgAqIj2HH1J3jTEM6RPJkD6R3HDOAKrr6lmzv4gPd+axfEcev3trO797azup0U7OG57E+cOSmKzRUZHezdtrk7eVxI85yfzRptezayb1w1rLwcIKvthfxBf7CvniQCFLt+cCEOQwjEqL5vS+sYzvG8PpfWNJjQltWztERJpQCq6ISAcWCCk6++e8EzCN5Tty+WRXPuU19QQ7Apg4wDU6ev7wJFXX7SaUgiud7lSq53ZQ5d3C8hrWHChi9YFC1uwvYlNWMdV1rrWW+0SFMC4jhrEZMYzLiGF0WrTWWBaRRqqCKyLiibc/1E5yXE1dA6v3F7JsRy7Ld+SxK7cMgIy4UM4bmsT5wxM5e2ACocFKa+uKFIBKl+DtMlgdPH/0mJq6BrYfLmHtgSLWHTrKhkNH2V9Q0fj+wMRwRqdFMzotmlFp0ZyWGqU1lkV6KQWgIiKt8eZHXjt/4B0qrGD5zjw+3JHLp7sLqKytJzgwgLMGxHHu0ETOG5bIoETNHe0qFIBKtzY/BmjpN56B+Uc79FRF5TVsyDzKpsxiNmYVszmrmJziqsb3+8aFcVpKFCNSohieEsmI5CjSY0NV5Eikh1MAKiLS0U7hB15VbT1f7C9k+Y48PtyZx2736GhaTChThyQwbWgiUwYlEB2mkYPOogBUujUfjYC2VV5pNVuyi9mSXcLWnBK2Zpewv6CcYz87w4IdDEmKcM1BTYpgSJ8IBiVGkB4bhkOBqUiPoCq4IiIdzdsCIRtfxLl0IVOLM5kanc6vL5hHZsalfLQznw935vLmxhye/+IQAQbGZsQwdXAC5wxJZHzfGIIcAb75LiLSs3hbebeDJEaGcN6wJM4bltS4raKmjp1Hyqhc8xwjttxHdH4uh/Pj+f26ufy+4RwAggMDGBAfzoCEcPonhDMgIYz+8eH0iw8nKTJEo6YiPYDPRkCNMU7gIyAEV6D7srX2nmb7XAf8Gchyb3rIWvvv1j5Xd25FpMvwZu5oG46pq29gQ+ZRPtzpWuplw6GjNFgID3YwaWA85wxJ4JzBCQxOUrquL2kEVLo9b+eP+rpNza6BNjCUfZN/x+rImezJK2NPXhl788s5VFhBbf2Xv1NDAgPIiAujb1wY6bGhZMS6/qbHhpEa4yQuPFjXRJEuxO8puMZ1BQi31pYZY4KAT4A7rbWfN9nnOmCCtfb2tn6uOk4R6VLa+wPPi7S44spaVuzJ5+Nd+XyyO58D7oIffaJCmDLYFYxOGZxAnyhnR3wjcVMAKuID7bgG1tU3kH20in0F5RwsrOBQYQUHCso5WFhJZmEFpdV1x+3vDAogNSaU1OhQUqKdpMSEkur+mxLtJDnaqYJIIn7k9xRc64psy9wvg9yP7jXhVETkZMbMbd+IghcLy0eHBnHxqBQubvgY9i/EOjOpcCbzYswNPLhjPK+udSWRDEmK4JwhCUwbkshZA+MIC9YsCxHpYtpxDQx0BNA3Poy+8WEtH1JZS2ZRBZlFlWQfrSSrqJKso5VkF1exa1ceuaXVNB9nCQ920CfaSUq0kz5RTpKjXIFp0+cJESGahyriQz79dWKMcQBrgMHA3621K1vY7RvGmGnATuAua+0Jt8WMMTcDNwP07dvXhy0WEfGxU5g3eixtzQDhVTlcX3Av185+gK0JF/Hpbtfo6LMrD/L4p/sJchjO6BfLtKGJnDs0kdNSopSaJiKdz9trYEsfFRpEdGg0I1OjW3y/tr6BIyVVHC6uIqfY9Te7uLJx2+d7CjhSWk19w/FRqiPAkBgRQp9oJ8lRIe7ANJTk6BD6RDlJiQ4lOcqpJbREvOSXKrjGmBjgf8APrLWbm2yPB8qstdXGmO8B37TWTm/ts5Q6JCLdmrdrjrYxba2qtp7V+4v4eFceH+3KZ1tOCeAqCDJtSCIzRiQxdUiCFotvA6XgiviAt9dAX7Rj6UJscSYNkWlknfETdiR9hcMlVRwprnL9dQeqh0uqKK2qO+EjopyBrmDUPaKaHO10pf/GuILU1BinMlGkV+vUKrjW2qPGmGXAxcDmJtsLmuz2b+BP/miPiEinOfYDq72FQdqYtuYMcriKFA1J4BdAbkkVH+3K56OdeSzdfoRX1mYS5DBMHBDHjOF9uGhUMmkxoaf+vURE2sLba2BHapZR4ijNpO+nP6fvrDCY1HI7yqvrXAFpk6A05+iXr7dkl5BfVn3CcTFhQaREh5IW4yQ1JpS0mFDX39hQ0mNDSYwIUXaK9Dq+LEKUCNS6g89Q4D3gj9baN5rsk2KtzXE/vxz4mbV2Umufqzu3ItIrncqafk3u9NeEp/JO8s08lDeeXe61R8emR3PxqBS+OjrF41yr3shfI6CnWjXeGHMt8Cv39kXW2idbO5/6Uen1fLRGanVdPbkl1WQfrSTHnfabfdQ9P9X9t6TZSGpIYEBjJd+MuFD6xoWRERtGRlwY/eLDlK0i3VpnjICmAE+654EGAC9aa98wxiwEVltrFwN3GGNmA3VAIXCdD9sjItJ9ebumX7M7/SHlWVx28A9cNusB9qV+lbc35/DO5sP88Z3t/PGd7ZzRL5avjU/j0tEpxIYH+/QrSaNqYHrTqvHGmLebVo13e6F51XhjTBxwDzABV6G/NcaYxdbaIr+0XKQ78qIYXFuEBDrIiHMFj42aVUqvnPVLDqZdSmZRBVlHK8ksqiSzqIJDhZVsyDzK0Yra4z4zLjyYfvGutVD7x4fTPyGMgQkRDEgMJyJE6b3SPfllDmhH0p1bEem1vFnTr413+jOLKnh9Qw7/W5fJziNlBDkMM0/rwzWT+nH2wPhemSLWGXNAjTFhuJYtu7Vp4T5Py5YZY64CzrPWfs/9+mFgubX2OU/nUD8qvZ6PRkBP4MV815KqWg4VVnCwoIIDhRUcKHAtPXOgwBWwNtUnKoSBCREMSgpncGIEQ/pEMjgpgqRIpfVK19Cpc0BFRKQDtHfJF2jznf702DBuPW8Qt5w7kK05Jby6NouX12Ty1qbDDEoM55pJ/ZgzIUN33H3kFKrGpwFNf0lnurc1/3xVkxc5xtuMkvZauvD4c4Dr9dKFx1/Lm9xcjIpOZ+SMeYxs4VpfVVvPgYIK9uWXsSevnL155ezNL2Px+uzjUnujnIEMS450PfpEMjwliuHJkUrnlS5DvyRERHqydi55YIxhZKprWYOfXDSMNzbm8PTnB1jw+lbuX7KL66f057rJ/YkJU3puR7LW1gPjjlWNN8aMalo1HngdeK5J1fgngVarxjf7/EeAR8A1AtqBTRfpfvxVCKktNwCbj5IWH3K9btpON2eQozGwbMpaS15ZNbuPlLErt4ydR0rZcbiU19ZnH1e9NyMulBHJUYxMjWZ0ehSjUqNJinKe8tcUaS8FoCIiPZk3d/rdd+OdxZlcEZ3OFTPmsS5mJn9ftof7l+zi0Y/28p3J/bnl3EFEh+qOekfyomp8FnBek/fSgeW+baVID+BNRkl7teUGYFtHSY9pYSqGGTOXpEgnSZFOJg9OaNzVWktOcRXbD5ewLaeUrTklbMsu4f1tRzg2Ay8pMoQx6TGMy4hmbEYMY9JjdF0Xn1MAKiLSk7X3Tr+Hu/HjZz3Av6+dy7acEv6xfA//+nAPL3xxiB9fOJQrz+yLI0DzjbzVQtX4mcAfm+3TWDUemA1scz9/F/idMSbW/fpC4Bd+aLaInExbbgC2pyBSO0ZLwZXRkupe9mX68D6N28uq69iaXcKmrGI2ZxWzIfMoS7YdaXx/SFIEZ/SL5fR+sZzRL5aBCeGaUyodSkWIRETkS20szrE5q5iFr29l1f5ChidHMn/2SCYNjPdjQ33Pj8uwjMGVUtu0avzCplXjjTG/xxV4Hqsaf6u1drv7+BuA/+f+uN9aax9v7XzqR0X86GTF49pTEMnb4kltKGBXXFnLpsxi1h0sYu3BItYePEpxpasib0JEMBMHxDGxfxyTBsUzNCmSAN10lDbw1I8qABURkS/Nj8G1mkdzBuYfPW6LtZa3Nx/md29tI7OokhumDOCnFw/DGeTwS1N9rTOq4PqD+lGRLqQ9lXLbcX326vObaGiw7M0vY/X+IlbtK2TlvsLGKrzx4cFMGhTP5EHxTB2cqPWjxSNVwRURkZNrR9EiYwyXjE7h/GFJ/P7tbTz26T4+3pXH/VeOY2RqtB8aKyLSzbVnmkQ7i8o1fm575pi6BQQYBidFMjgpkisnuipnZxZV8PneQj7bk89nuwt4c6NrVkC/+DCmDUlk2tBEpgyOJyxY4YW0TiOgIiLyJS/vlgMs35HLT1/eSFFFDf/vkhFcN7l/t543pBFQEelSvLk+eztqepKA2FrL3vxyPt6Zx8e78lmxt4CKmnqCAwOYNDCeGcOTmD48iYw4jY72ZhoBFRGRk/N2eYKNL3Le0oWsrM2kwJnIwjevYF/+t5h36WkEOgJ8324RkZ7Om+tze0dN21joyBjDoMQIBiVGcN2UAVTX1bNmfxEfbM/lg+253LN4C/cs3sKIlCguGtmHi0YmMzw5slvflJSOoxFQERE5NS3cla8JcHJ31Q2UDrmcB68+nYiQ7ne/UyOgItLttXfU1NtCR83syy9nydYjvLf1MKsPFGGtK1X3ktEpfHV0CiNToxSM9gIqQiQiIr7h4QdLmTOFsSX3MrRPJE9ef2a3W/BcAaiI9AhtSKlt1N6U3TZ8dl5pNUu2HeHtzYf5dHc+9Q2WfvFhzBqTytfGpzI4KfKUv6J0Te0OQI0xX2/D51ZZa9861ca1hzpOEZEuppUfLB9evYvvP72GtNhQXvze2cSEBfu7dV5rawDaVftLT9SPiohH7RkB9WJOamF5De9tOcwbG3P4bE8+DRZGpkbxtXFpXDY+laTI7nWjUlrnTQBaALwGtDY+Ps1aO6hjmtg26jhFRLqYk/xg+Wx3Ptc9/gUj06J45sazuk2FxHYEoF2yv/RE/aiIeNSeoPIU03VzS6t4Y0MOr63PYkNmMY4Aw7lDE7nijHRmjEgiJLBnLOnVm3lThOhta+0NJ/nQp0+5ZSIi0r3NmNfyD5YZ8wCYPDiBB68ez61Pr+F7T63h39dO6Gk/LNRfikjP0J5CR8WZLX+Gp+3NJEU6ueGcAdxwzgD25JXxyppMXl2bxfe3ryUmLIivj0/nqokZDOmjFN2eprUR0CBrba2f23NSunMrItIFtWEe0IurD/HTlzdyyehkHrzqdBwBXbsARTtGQLtkf+mJ+lER6RAdUbCoWd/RMH0en4SezwtfHOK9rYeprbdM6BfLVRP78tUxKTiDetTNyx7PmxTcXGAx8Bzwge0i1YrUcYqIdF+PfrSX3761jbsvHMrt04d0dnNa1Y4AtEv2l56oHxWRDnEK60a35fj8smpeXZvJc6sOsS+/nNiwIL55Zl++dVZfrS/aTXjqR1tbnG0E8AXwK+CQMeZvxphJvmqgiIj0fDdOHcCssanct2QXq/cXdnZzOor6SxHpfcbMdQWL0RmAcf1ta/AJrpHPpsEnuF4vXQhAQkQIN08bxAc/PpdnbjyLiQPieOSjPUz78zJufHI1n+3Jp4vf7xMP2rQMizEmFZgDXAkkAc9ba3/p47a1SHduRUS6oSZpVg1RaSyqnMO7jmm8dcdUosOCOrt1LfJmGZau1F96on5URLoEL5Z8qXt/AY7SLHKI5w81c9nV5xJumNKf2eNSe1ptgR7BmxHQRtbabOA/wD+BUuDGjm2eiIj0WMfSrIoPAZaAkkx+1fAvJpYt4WevbOxRd7DVX4qItFF0etu3u/uRwNJMDJZU8rk39DHOrVrGT17eyNQ/LuNfH+6hpKrbTMfv1VoNQI0xTmPMHGPMq8BuYDrwcyDVH40TEZEeoIU0q4D6ShaGv8o7Ww7zzMqDndSwjqP+UkSknWbMc835bKpJBfXjtNCPBNZX8bPgF3nquxMZlhzJH97ezuTff8Dv39pGbkmVDxsup8rjMizGmGeBC4APgWeAq621+l9TRETax0NJ/ojqw5w7NJFFb27l/OFJpMWEtrhfV6f+UkTECx2w5IspzmTqkESmDklkc1Yxj3y0l0c/3svjn+3nyjMz+N65g7pt39KTtbYO6DvA96y1pf5qjIiI9EDR6S2W6jfR6fz28lHM+OuH/Omd7fztyvGd0LgOof5SRMQbY+a2rWiRh36kabruqLRoHrhqPD++cCj/XL6H0i+ehbUvYE0B9ZFpBM68p+0FksSnWkvBLTxZZ2qMubSD2yMiIj1NK2lW6bFh3DxtIK+tz2bNgaLOad+pU38pIuJL7UjX7Rcfzh+GbOevzsdIM/kYLIGlmdT873aKPn/aTw2W1rQ2AvpnY0wW0NpK4b8D3ujYJomISI9ykjSrW84dxAtfHGLhG1v5362TCQhordvpktRfioj4UnvSdd37BdQdP2c02FaT9/Y9PJA3ntvPH0x8RIiPGy2etBaAHgHuPcnxuzqwLSIi0lO1kmYVHhLITy8ezt0vbeC1DVlcPt5DZcSuS/2liIivtTVdFzzOGU01Bfx3xQFeWp3JTVMHcuPUAYSHtBYOiS94/Be31p7nx3aIiEgv9vXxafx3xX7++PYOLhqZTFhw9/lBoP5SRKSLaaX2wLvfmsZf3t3BfUt28tTn+7nzgqFcdWYGgY42rU4pHUD/0iIi0nk2vgj3jSJgYSwvVd7MxLIl/OvDvZ3dKr9xL9+yyhizwRizxRizoIV9fmSM2WqM2WiMWWqM6dfkvXpjzHr3Y7F/Wy8i0kW1Mmd0cFIE//r2Gbz6/ckMTIzg1/+3mT//5TdU/WkEzI+B+0a5+ibxGQWgIiLSOdwLi7vuUltCyrP4c8h/OPLpfymvruvs1vlLNTDdWjsWGAdcbIyZ1GyfdcAEa+0Y4GXgT03eq7TWjnM/ZvunySIiXdyYuTDrAYjOAIzr76wHjkvhPb1vLC/cPInXpmVxV+XfcVZkA9bVJ71+h4JQH+o+OU4iItKztLCweIit5gcNz/Hymhu5dnL/zmmXH1lrLVDmfhnkfthm+yxr8vJz4Br/tE5EpBtrw5xRYwxjdzyA615gE7WVNCxZQICWbfGJk46AGmPCjDG/NsY86n49ROXkRUTklHkqEhFQwOOf7qOhwbb4flflbX9pjHEYY9YDucD71tqVrez+XeDtJq+dxpjVxpjPjTFfa+UcN7v3W52Xl9fGbyQi0gt46IsoyeK9LYf925Zeoi0puI/jui1wtvt1FrDIZy0SEZHeIbrlardVoSnsL6hg6fZcPzfolHnVX1pr662144B0YKIxZlRL+xljrgEmAH9usrmftXYCcDVwvzFmkIdzPGKtnWCtnZCYmNjmLyQi0uN56IvyTAI3P7WGW55aQ25JlWuju26B5oqemrYEoIOstX8CagGstRW0vtaZiIjIyXkoEhFy0XxSo53855NuV4zolPpLa+1RYBlwcfP3jDEXAL8EZltrq5sck+X+uxdYDow/hfaLiPQ+Hvqi+Mt+y88uHs6yHbnMuPdDPvvfP7FN6hZorqj32hKA1hhjQnHPSXHfXa1u/RAREZGT8FAkwjHum1w7uT+f7y1kS3ZxZ7eyPdrdXxpjEo0xMe7nocBMYHuzfcYDD+MKPnObbI81xoS4nycAU4CtHfd1RER6AQ99UeC4b3LreYN454fTGJkaRd91f8E0q1tAbaWrnoG0S1uKEN0DvANkGGOewdXBXefLRomISC/hoUjElRP78relu/jPJ/u4d+64TmiYV7zpL1OAJ40xDlw3hV+01r5hjFkIrLbWLsaVchsBvGSMATjorng7AnjYGNPgPvYP1loFoCIi7dVKwaIBCeE8e+MkzMKClo/1NIdUPDppAGqtfd8YsxaYhCuV6E5rbb7PWyYiIr3Pxhdh6UKiizP5LCSJBRuvIPcrw0mKdHZ2y07Km/7SWruRFtJmrbXzmjy/wMOxnwGjT6nRIiJyUgEBxjVXtPjQiW96mEMqnrWlCu7pQD8gB8gG+hpjBhljtISLiIh0nGbrgsbUHuG3jkdZ8/ojnd2yNlF/KSLSg7UwV7SSELaPvOvEfVWsqFVtmQP6D1zrjj0CPAqsAF4CdhhjLvR0kDHGaYxZZYzZYIzZYoxZ0MI+IcaYF4wxu40xK40x/b36FiIi0v21sC5omKnh9F0PdFKD2s2r/lJERLqBZnNFayPSuM95G19Zlsxf3t1BXX2Da79mN1NVrOhEbQlAs4Hx7vLtZ+BKFdqLq1DCn1o5rhqYbq0dC4wDLjbGTGq2z3eBImvtYOA+4I/t/QIiItJDeJhHk9iQz84jpX5ujFe87S9FRKQ7GDMX7toM848SdPdWfvijXzLnjHQeWrabqx79nJziyhZvpqpY0fHaEoAOtdZuOfbCXeBguLvku0fWpcz9Msj9aL6q+GXAk+7nLwMzjLvCgoiI9DIe5tFkE8+bG3P83BiveNVfiohI9xQWHMifrhjL364cx9bsEmY9+AnWU1EiFStq1JYAdIsx5p/GmHPdj38AW92l32tbO9AY4zDGrAdygfettSub7ZIGHAKw1tYBxUB8C59zszFmtTFmdV5eXhuaLCIi3Y6HtdhejbmBtzZ1iwDU6/5SRES6r8vGpfHa7VOIcgaRbU8IZVxUrKhRWwLQ64DdwA/dj73ubbXA+a0daK2tt9aOA9KBicaYUd400lr7iDulaUJiYqI3HyEiIl2dh7XYYiZ9i125Zd0hDfc6vOwvRUSkexucFMn/3T6FNxJvosIGH/9mUKjrJqsAbVuGpRL4q/vRXFkL21r6jKPGmGXAxcDmJm9lARlAprtKYDTgYZEdERHp8VpYi+3i0iruWbyFNzfmMHRmZCc17OQ6or8UEZHuK8oZxE3f/xlvPxfK2B0PkBpQgI1Kw3HBPR7XGe2N2rIMyxBjzMvGmK3GmL3HHm04LtEYE+N+HoqrCMP2ZrstBq51P78C+MBa23yeqIiI9GJJkU4m9o/r8mm43vaXIiLScwQEGL76rTtZ8/WPGF77HDPt3zmYdumJO/bipVrakoL7OPBPoA5XCtF/gafbcFwKsMwYsxH4Atcc0DeMMQuNMbPd+/wHiDfG7AZ+BPy8vV9ARER6vkvHpHSHNFxv+0sREelhLhuXxlPfnUhBWQ1f/+enrD909Ms3e/lSLW0JQEOttUsBY609YK2dD3z1ZAdZazdaa8dba8dYa0dZaxe6t8+z1i52P6+y1s6x1g621k5UpUAREWnJRaOSMYauXg3Xq/5SRER6prMGxvPKrZMJDXZw5SMr+Ginu5hqL1+qpS0BaLUxJgDYZYy53RhzORDh43aJiEhv1yQ9KenfE/hRn/VdPQ1X/aWIiBxncFIEr946hQEJEdz45Gre3XLY85IsvWSplrYEoHcCYcAdwBnANcB3fNkoERHp5VpIT7q15AFG5L/TldNw1V+KiMgJEiNDeP6mSZyWGsX3n1lLRWhyyzv2kqVa2hKA9rfWlllrM62111trvwH09XXDRESkF2shPSmwoYqfBr7Iij1dtli6+ksREWlRdFgQT994FhP6xfKLksupcziP36EXLdXSlgD0F23cJiIi0jE8pCGlBhQcX8iha1F/KSIiHkWEBPLE9RM5OuhyflR5A+WhKTRd97pxqZYeXiHX4zqgxpivAJcAacaYB5q8FYWrwp+IiIhvRKe702+PVxSY2OUCUPWXIiLSVqHBDh7+9hl898kGRu85h4euPp1LRqd8ucOxKSjHsoCOVciFHrOWaGsjoNnAGqDK/ffYYzFwke+bJiIivdaMea50pKaCQlkz5A725ZdztKKmc9rVMvWXIiLSZs4gB49+ZwKn943ljufW8cH2I1++2Qsq5HocAbXWbgA2GGOettbqDq6IiPjPsbu8Sxe60nGj02HGPCLCp8O6law/dJTzhiV1bhvd1F+KiEh7hQUH8tj1Z/KtR1dyy9NreeK6M5k8OKFXVMhtLQV3E2Ddz09431o7xnfNEhGRXm/M3BPSjcZU12EMrDvYdQJQ9ZciIuKNKGcQ/71hIt98ZAXfe2oNL986mWEepqD0pAq5HgNQ4FK/tUJERKQNIkICGdYnsqvNA1V/KSIiXokND+bx6ydy+d8/5frHV/H29F8QveTHx6fh9rAKuR7ngFprDxx74JrXMtr9qHRvExER8btxGTFsyDyKtbazmwKovxQRkVOTFhPKY9edydHKWq5Z1Y/qr9zvqozbUoXcHuCky7AYY+YCq4A5wFxgpTHmCl83TEREpCXjMmI4WlHL/oKKzm7KcdRfioiIt0alRfPgVePZkl3MbZsGUX/nJph/FO7afHzw2QOWaGnLOqC/BM601l5rrf0OMBH4tW+bJSIi0rJxfWMAWHewcW5KQQAAG7xJREFUqJNbcgKv+ktjjNMYs8oYs8EYs8UYs6CFfUKMMS8YY3YbY1YaY/o3ee8X7u07jDGquisi0k3NGNGHBbNHsmRbLn96d/uJOxxboqX4EGC/XKKlmwWhbQlAA6y1uU1eF7TxOBERkY7R5I7vsOcmMyf4s642DxS87y+rgenW2rHAOOBiY8ykZvt8Fyiy1g4G7gP+CGCMOQ24EhgJXAz8wxjjOLWvISIineXbZ/fnW2f15eEP9/LO5sPHv9lDlmhpS8f4jjHmXWPMdcaY64A3gbd82ywRERG3Znd8TfEhFjkeJWrX/zq7Zc151V9alzL3yyD3o/kE18uAJ93PXwZmGFfJ3cuA56211dbafcBuXCOvIiLSTc2bdRpjM2K4+6UN7M0r+/KNHrJEy0kDUGvtT4CHgTHuxyPW2p/5umEiIiJAi3d8Q2w1V5U+QVVtfSc16kSn0l8aYxzGmPVALvC+tXZls13SgEPu89QBxUB80+1ume5tzT//ZmPMamPM6ry8/9/e3UfJVZcJHv8+/VYJ3SGBpEVIAiRMwMEIgY0MLMirIjgKegbnxHHUGXfl6OIKC84M6iyLnHHPOOPqinrksKCO4ysrqOCgwiqDOLMiLxMghLcQYAEDdECTdAh5ffaPuh2KprvzQlfVvd3fzzl16ta9t6qeJ7dSt5/6/e7vN7BriUmSWqrW1cmX330kPV0dfPAbd/D8pmKK6dGmYqnYFC2jFqAR8aWIOBYgM6/JzPOLW+l+cpYkTWCj/LK7L89y72/WtjiYlxuP82Vmbs3MRcAc4KiIWDieMWbm5Zm5ODMX9/f3j+dLS5KaYL8ZU7l0yRGseGaQC6++pz7y+ykX1adkaVTBKVrGagF9EPhMRDwaEX8XEYtaFZQkSduN8svub3JmWQYiGrfzZWb+DriJ+vWcjZ4E5gJERBcwnfo1ptvXF+YU6yRJFXfcgllccOohXHvXb7jq9sfro+G+7dLKT9Ey1jygn8/MY4ATqJ/kvhoR90fEf4uIg1sWoSRpchvlF98rev60FAMRvdLzZUT0R8SMYnkq8CZg+PCH1wLvK5bPAn6e9YlQrwWWFKPkzgMWUJ8KRpI0AXzohIP49wfN5JPXLeexZ9fXi83/smzkKVoqYmeuAX0sMz+dmUcA7wLeDtzX9MgkSYJRf/EdmHdmKQrQIa/gfLkvcFNE3A3cRv0a0B9FxCURcUaxz5XAzIhYAZwPXFi8573AVcBy4CfAOZlZngtjJUmvSEdH8Jl3Hk5nR3D+VXexdduwMeoqOC9o1452KLr6nE59mPdTgH8GLm5qVJIkNTrsj1/2K+/rfvsw/3TPKtZs2Mz0qd1tCuxFu3u+zMy7gSNGWH9Rw/ILwDtHef6ngE/tTsySpPLbb8ZU/ubtCzn3O0u57OaHOeek36tvGBolfmigvqF5QaHULaNjDUL0poj4CvUR9T5AfTj5gzJzSWb+sFUBSpI0kr17ewBYu2FzW+PwfClJarYzDt+Ptx62L5+78UGWPbmmvrKi84KO1QX3Y8C/Ar+fmWdk5rcyc32L4pIkaUzTavVOPOuHhqdvH8+XkqSmigj+5u0LmdVX47zvLmXjlq2VnRd0rEGITs7MKzKzFEMMSpLUqHeoAN3Y3gLU86UkqRVm7NHD3/7R61jxzCCX37yysvOC7nAQIkmSymioAF33QttbQCVJaokTD3kVf3jYvnzhphUMHPVXlZwX1AJUklRJ06YMtYA66KskafK46K2H0tPZwfn3H0xWcF5QC1BJUiUNtYAObmzvIESSJLXSPntO4aOnHswtD63mnzjupfOCQumnZbEAlSRVUl/PUAFqC6gkaXJ5zzEH8rrZ0/nkdctZ+0LxQ+zQtCxrHgfyxWlZSlaEWoBKkiqpt9YJtH8QIkmSWq2zI/jv73gdzw5u5LM3PFhfWZFpWSxAJUmV1NXZwZTuDgYtQCVJk9Dr5kznT/5gf77xq8d4ZPX6ykzLYgEqSaqsvlq3BagkadI695SDqXV18Pc/vb8y07JYgEqSKquv1smg07BIkiap/mk1PnD8fK6/5ykeXXRBJaZlsQCVJFVWb63La0AlSZPaB94wn1l9Nf7ywddUYlqWrnYHIEnS7uqrddkFV5I0qfXWujj3jQv4rz9Yxs+PP4FThqZjKSlbQCVJlWUBKkkSLHn9XObP6uXTP7mfrduy3eGMyQJUklRZfVPsgitJUndnB3/x5kN48OlBrr6zXKPeDmcBKkmqrN5aF4Mbt7Y7DEmS2u60ha9m4ew9+fI/P/xiK+jdV8HnFsLFM+r3d1/V3iBpYgEaEXMj4qaIWB4R90bEuSPsc2JErImIpcWtXEM0SZJKrd4Fd3O7w5Akqe0igg+d8Hs8sno9N9z7VL3YvO4jsOZxIOv3132k7UVoMwch2gJckJl3RsQ04I6IuDEzlw/b75bMfGsT45AkTVB9tS5e2LyNLVu30dVppx5J0uR22sJXc+DMPbjs5oc5bfMlxOYNL91h8wb42SVtHRm3aWfrzFyVmXcWy+uA+4DZzXo/SdLk01ur/4663m64kiTR2RF84Pj53PXEGlgzyrWgo61vkZb8XBwRBwJHALeOsPmYiLgrIn4cEa8d5flnR8TtEXH7wMBAEyOVJFVJX60TgMFNDkQkSRLAHx05h1l9NZ7t7B95h+lzWhvQME0vQCOiD7gaOC8z1w7bfCdwQGYeDnwB+MFIr5GZl2fm4sxc3N8/yj+kJGnS6at1A1R2JNydHC/hLxrGSlgWEVsjYu9i26MRcU+x7fbWZyBJKpsp3Z28/7gDuWTDWWzrmvrSjd1T4ZT2DrvT1AI0IrqpF5/fzMxrhm/PzLWZOVgsXw90R8SsZsYkSZo4eosW0HUvVLMA5cXxEg4FjgbOiYhDG3fIzL/PzEWZuQj4GHBzZj7XsMtJxfbFrQtbklRmf3r0AdzUfSJfn3U+TJ8LRP3+bZe29fpPaOIgRBERwJXAfZn52VH2eTXwdGZmRBxFvSB+tlkxSZImlmlThq4BrWYBmpmrgFXF8rqIGBovYfiAfUPeBXy7ReFJkipqzynd/MnR+3PJL7Zw0kdv5YCZve0OabtmtoAeC7wHOLmh69BbIuKDEfHBYp+zgGURcRdwKbAkM7OJMUmSJpChQYgGK1qANtrBeAlExB7AadR7Fg1J4IaIuCMizh7jtR1LQZImmfcfO4+I4Du3Pd7uUF6iaS2gmflLIHawzxeBLzYrBknSxNbbMzEK0B2MlzDkbcC/DOt+e1xmPhkRrwJujIj7M/MXw5+YmZcDlwMsXrzYH3olaRLYZ88pnHhwP9fc+QQXvOng0kxXVo4oJEnaDVXvggs7Hi+hwRKGdb/NzCeL+2eA7wNHNStOSVL1vHPxHJ5eu5FbHloNd18Fn1sIF8+o3999VVtisgCVJFXW9i64FR2EaGfGSyj2mw6cAPywYV1vREwbWgZOBZY1N2JJUpWc/Jp92Lu3h0du+ipc9xFY8ziQ9fvrPtKWIrRpXXAlSWq27s4Oero6qjwP6NB4CfdExNJi3ceB/QEy87Ji3TuAGzJzfcNz9wG+X69h6QK+lZk/aUnUkqRK6Onq4B1HzObNt30QYsNLN27eAD+7pOWj4lqASpIqbVqtq7JdcHdmvIRiv68BXxu2biVweFMCkyRNGO9cPId9b1s98sY1T7Q2GOyCK0mquN5aV2W74EqS1GyvefWerO7sH3nj9DmtDQYLUElSxfXVuhjcuLXdYUiSVFoPLDyf57PnpSu7p8IpF7U8FgtQSVKl1QvQze0OQ5Kk0jrstA/w19vO5nfd+wAB0+fC2y5t+fWf4DWgkqSK6611snpwU7vDkCSptKbv0c2WQ8/ihAdP4Nd/fQq1rs62xWILqCSp0vqmdFd2ECJJklrlbYfvx5oNm7njsd+2NQ4LUElSpfXVOllnASpJ0piOOWgmXR3BLx4cZUTcFrEAlSRVWl+Fp2GRJKlV+mpd/LsD9uKWhwbaGocFqCSp0nprXTy/aStbt2W7Q5EkqdSOP7ife3+zltWDG9sWgwWoJKnS+mr18fTWb7IVVJKksbxhwSwA/mVF+7rhWoBKkiptewFqN1xJksb02v2ms9ce3dz8YPu64VqASpIqrbcoQAdfsACVJGksnR3BcQv6ueWh1WS259IVC1BJUqUNtYAO2gIqSdIOvWHBLAbWbeSBp9e15f0tQCVJldY3ZagL7tY2RyJJUvkNXQd6S5umY7EAlSRVWm/PUAvo5jZHIklS+e07fSoLXtXHL9o0HYsFqCSp0qZNGSpAbQGVJGlnHH9wP7c+8hwvbG79udMCVJJUaS8OQmQLqCRJO+MNC2axacs2fv3Icy1/bwtQSVKl9dY6AVi/yRZQSZJ2xh/Mm0lPZwe3tKEbrgWoJKnSal2d9HR2OAquJEk7aWpPJ6+ftxe/XPFsy9/bAlSSVHm9tU7nAZUkaRcsnD2dh58ZZOu21s4HagEqSaq83loX620BlSRppx00q49NW7fx5G83tPR9LUAlSZXXV+uqZBfciJgbETdFxPKIuDcizh1hnxMjYk1ELC1uFzVsOy0iHoiIFRFxYWujlyRV2bz+XgAeXj3Y0vftaum7SZLUBFUtQIEtwAWZeWdETAPuiIgbM3P5sP1uycy3Nq6IiE7gS8CbgCeA2yLi2hGeK0nSy8yfVS9AHxlYz0mHtO59bQGVJFVe35RqdsHNzFWZeWexvA64D5i9k08/CliRmSszcxPwHeDM5kQqSZpo9u7tYc8pXaxscQuoBagkqfJ6a12sq2AB2igiDgSOAG4dYfMxEXFXRPw4Il5brJsNPN6wzxPsfPEqSZrkIoL5/X2sHFjf0ve1AJUkVV5fTzVbQIdERB9wNXBeZq4dtvlO4IDMPBz4AvCD3Xj9syPi9oi4fWCg9XO+SZLKaf6sXh5ZbQEqSdIuqXfB3druMHZLRHRTLz6/mZnXDN+emWszc7BYvh7ojohZwJPA3IZd5xTrXiYzL8/MxZm5uL+/f9xzkCRV0/z+XlateYHnN7XuR1wLUElS5fUWgxBta/FcZq9URARwJXBfZn52lH1eXexHRBxF/dz9LHAbsCAi5kVED7AEuLY1kUuSJoL5/X0ALW0FdRRcSVLlTavVT2fPb95KX61Sp7ZjgfcA90TE0mLdx4H9ATLzMuAs4EMRsQXYACzJzAS2RMSHgZ8CncBXMvPeVicgSaquecVIuCsH1vPa/aa35D0rdZaWJGkkvUXROfjClkoVoJn5SyB2sM8XgS+Osu164PomhCZJmgSGCtBWtoDaBVeSVHm9tU6Aqs4FKklSW0zp7mT2jKmsHGjdVCwWoJKkyps2pd7qWeWRcCVJaof5/a0dCdcCVJJUeb09RRdcC1BJknbJvFm9rBxYT314geZrWgEaEXMj4qaIWB4R90bEuSPsExFxaUSsiIi7I+LIZsUjSZq4tl8DagEqSdIumT+rl3Ubt7B6cFNL3q+ZLaBbgAsy81DgaOCciDh02D6nAwuK29nAl5sYjyRpgrILriRJu2doKpZWXQfatAI0M1dl5p3F8jrgPmD2sN3OBL6edb8CZkTEvs2KSZI0MdkCKknS7tk+FUuLrgNtyTWgEXEgcARw67BNs4HHGx4/wcuLVCLi7Ii4PSJuHxgYaFaYkqSK6rMAlSRpt8yeMZWero6WDUTU9AI0IvqAq4HzMnPt7rxGZl6emYszc3F/f//4BihJqrxaVwddHcHgCxagkiTtio6OYN7M3up3wQWIiG7qxec3M/OaEXZ5Epjb8HhOsU6SpJ0WEfTWurwGVJKk3TC/v7f6XXAjIoArgfsy87Oj7HYt8N5iNNyjgTWZuapZMUmSJq6+WheDG7e2OwxJkipn3qxe/t+zz7N567amv1dXE1/7WOA9wD0RsbRY93Fgf4DMvAy4HngLsAJ4HvjzJsYjSZrA6gXo5naHIUlS5czv72PLtuSJ327YPihRszStAM3MXwKxg30SOKdZMUiSJo/eWifrbQGVJGmXze8vRsIdGGx6AdqSUXAlSWq2vindrPMaUEmSdtn8ouhsxUi4FqCSpAmhr9bpIESSJO2GGXv0sHdvDw8PNL8AbeY1oJIktcxJh7yKQ/bZs91hSJJUSe895gDm7rVH09/HAlSSNCG8c/HcHe8kSZJGdN4bD27J+9gFV5IkSZLUEhagkiRJkqSWsACVJEmSJLWEBagkSZIkqSUsQCVJkiRJLWEBKkmSJElqCQtQSZIkSVJLWIBKkiRJkloiMrPdMeySiBgAHhunl5sFrB6n12oXcyiPiZCHOZTDRMgBqp/HAZnZ3+4gxpvn0TGZT7mZT7mZT3m1K5cRz6OVK0DHU0TcnpmL2x3HK2EO5TER8jCHcpgIOcDEyUOjm2jH2HzKzXzKzXzKq2y52AVXkiRJktQSFqCSJEmSpJaY7AXo5e0OYByYQ3lMhDzMoRwmQg4wcfLQ6CbaMTafcjOfcjOf8ipVLpP6GlBJkiRJUutM9hZQSZIkSVKLWIBKkiRJklpiUhagEXFaRDwQESsi4sJ2x7MjEfFoRNwTEUsj4vZi3d4RcWNEPFTc71Wsj4i4tMjt7og4sk0xfyUinomIZQ3rdjnmiHhfsf9DEfG+EuRwcUQ8WRyLpRHxloZtHytyeCAi3tywvm2ft4iYGxE3RcTyiLg3Is4t1lfmWIyRQ9WOxZSI+HVE3FXk8cli/byIuLWI6bsR0VOsrxWPVxTbD9xRfm3M4WsR8UjDsVhUrC/d50njo53/l3bXeJ2XymI8v9/LYDy/I8skIjoj4t8i4kfF48rmExX8e3QsETEjIr4XEfdHxH0RcUxV84mIQxrOwUsjYm1EnFfafDJzUt2ATuBhYD7QA9wFHNruuHYQ86PArGHr/g64sFi+EPh0sfwW4MdAAEcDt7Yp5uOBI4FluxszsDewsrjfq1jeq805XAx8dIR9Dy0+SzVgXvEZ62z35w3YFziyWJ4GPFjEWpljMUYOVTsWAfQVy93ArcW/8VXAkmL9ZcCHiuX/BFxWLC8BvjtWfm3O4WvAWSPsX7rPk7dx+RxU7jxaxP2Kz0tluo3X93tZbuP1HVm2G3A+8C3gR8XjyuZDBf8e3UE+/wD8x2K5B5hR5Xwa8uoEngIOKGs+k7EF9ChgRWauzMxNwHeAM9sc0+44k/p/HIr7tzes/3rW/QqYERH7tjq4zPwF8Nyw1bsa85uBGzPzucz8LXAjcFrzo68bJYfRnAl8JzM3ZuYjwArqn7W2ft4yc1Vm3lksrwPuA2ZToWMxRg6jKeuxyMwcLB52F7cETga+V6wffiyGjtH3gFMiIhg9v6YbI4fRlO7zpHFRyfPoOJ2XSmMcv99LYRy/I0sjIuYAfwhcUTwOKpzPKCr5eYuI6dR/lLoSIDM3ZebvqGg+w5wCPJyZj1HSfCZjATobeLzh8ROM/cdsGSRwQ0TcERFnF+v2ycxVxfJTwD7Fcpnz29WYy5rLh4vuCl8Z6spABXIouvMcQf1X5Uoei2E5QMWORdEVaynwDPWi62Hgd5m5ZYSYtsdbbF8DzKTNeQzPITOHjsWnimPxuYioFetKeyz0ikyk41fFc+nLvMLv99IYp+/IMvmfwF8C24rHM6l2PhPl71Go9yAaAL5adJG+IiJ6qW4+jZYA3y6WS5nPZCxAq+i4zDwSOB04JyKOb9yYmcnYrRClU8WYC18GDgIWAauA/9HecHZORPQBVwPnZebaxm1VORYj5FC5Y5GZWzNzETCHeivSa9oc0i4bnkNELAQ+Rj2X11PvVvtXbQxR2i1V+S4cbiJ8vw+ZCN+RQyLircAzmXlHu2MZRxPp79Eu6l3yv5yZRwDrqXdR3a5i+QBQXFN8BvC/h28rUz6TsQB9Epjb8HhOsa60MvPJ4v4Z4PvUv5SfHmoqL+6fKXYvc367GnPpcsnMp4sT5Dbgf/Fi18fS5hAR3dT/OPlmZl5TrK7UsRgphyoeiyFFN5+bgGOod3vpGiGm7fEW26cDz1KSPBpyOK3oCpiZuRH4KhU6FtotE+n4VfFcut04fb+Xziv8jiyLY4EzIuJR6t3UTwY+T3XzmUh/j0K9xe+Jhl4836NekFY1nyGnA3dm5tPF41LmMxkL0NuABVEfhayHejP1tW2OaVQR0RsR04aWgVOBZdRjHho58n3AD4vla4H3FqNbHQ2saWh6b7ddjfmnwKkRsVfRvfLUYl3bDOsf/w7qxwLqOSyJ+ih284AFwK9p8+etuH7kSuC+zPxsw6bKHIvRcqjgseiPiBnF8lTgTdSv2boJOKvYbfixGDpGZwE/L369HC2/duVwf8PJLahfX9J4LEr1edK4qNR5dAeqeC4FxvX7vRTG8TuyFDLzY5k5JzMPpP5/5OeZ+W4qms8E+3uUzHwKeDwiDilWnQIsp6L5NHgXL3a/hbLmkyUYranVN+ojPz1I/dqCT7Q7nh3EOp/6CIN3AfcOxUv9uoCfAQ8B/wfYu1gfwJeK3O4BFrcp7m9T7xa5mfqvTP9hd2IG3k99kJUVwJ+XIId/LGK8m/p/3n0b9v9EkcMDwOll+LwBx1HvbnE3sLS4vaVKx2KMHKp2LA4D/q2IdxlwUbF+PvUCcgX1LjO1Yv2U4vGKYvv8HeXXxhx+XhyLZcA3eHEky9J9nryN22ehMufRhpjH5bxUltt4fr+X4Tae35FluwEn8uIouJXMh4r+PbqDnBYBtxefuR9QH5W9yvn0Um81n96wrpT5RBGEJEmSJElNNRm74EqSJEmS2sACVJIkSZLUEhagkiRJkqSWsACVJEmSJLWEBagkSZIkqSUsQCVJkiRJLWEBKpVMRMyMiKXF7amIeLLh8b824f3+LCIGIuKKMfaZWrz/poiYNd4xSJI0XjyPSuXW1e4AJL1UZj5LfXJkIuJiYDAzP9Pkt/1uZn54jJg2AIsi4tEmxyFJ0ivieVQqN1tApQqJiMHi/sSIuDkifhgRKyPibyPi3RHx64i4JyIOKvbrj4irI+K24nbsTrzHa4vXWRoRd0fEgmbnJUlSK3geldrPFlCpug4Hfh94DlgJXJGZR0XEucB/Bs4DPg98LjN/GRH7Az8tnjOWDwKfz8xvRkQP0Nm0DCRJah/Po1IbWIBK1XVbZq4CiIiHgRuK9fcAJxXLbwQOjYih5+wZEX2ZOTjG6/5f4BMRMQe4JjMfGv/QJUlqO8+jUhvYBVeqro0Ny9saHm/jxR+XOoCjM3NRcZu9g5Mmmfkt4AxgA3B9RJw8znFLklQGnkelNrAAlSa2G6h3IwIgIhbt6AkRMR9YmZmXAj8EDmteeJIklZrnUWmcWYBKE9tHgMXFIAjLqV+XsiN/DCyLiKXAQuDrzQxQkqQS8zwqjbPIzHbHIKmNIuLPgMVjDR/fsO+jxb6rmx2XJElV4HlU2jW2gEraAJy+MxNoA93Ur42RJEl1nkelXWALqCRJkiSpJWwBlSRJkiS1hAWoJEmSJKklLEAlSZIkSS1hASpJkiRJaon/D2rJJIuC9yMHAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n", + "\n", + "# plot the 1C results\n", + "t_sol = solutions[0].t\n", + "ax1.plot(solutions[0][\"Time [s]\"](t_sol), solutions[0][\"Terminal voltage [V]\"](t_sol))\n", + "ax1.plot(voltage_data_1C[:,0], voltage_data_1C[:,1], \"o\")\n", + "ax1.set_xlabel(\"Time [s]\")\n", + "ax1.set_ylabel(\"Voltage [V]\")\n", + "ax1.set_title(\"1C\")\n", + "ax1.legend([\"DFN\", \"Experiment\"], loc=\"best\")\n", + "\n", + "# plot the 5C results\n", + "t_sol = solutions[1].t\n", + "ax2.plot(solutions[1][\"Time [s]\"](t_sol), solutions[1][\"Terminal voltage [V]\"](t_sol))\n", + "ax2.plot(voltage_data_5C[:,0], voltage_data_5C[:,1], \"o\")\n", + "ax2.set_xlabel(\"Time [s]\")\n", + "ax2.set_ylabel(\"Voltage [V]\")\n", + "ax2.set_title(\"5C\")\n", + "ax2.legend([\"DFN\", \"Experiment\"], loc=\"best\")\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For a 1C discharge we observe an excellent agreement between the model and experiment, both in terms of the overall shape of the curve and the capacity. The agreement between model and experiment is less good at 5C, but in line with other implementations of the DFN (e.g. [2]). " + ] + }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", - "[1] Ecker, Madeleine, et al. \"Parameterization of a physico-chemical model of a lithium-ion battery II. Model validation.\" Journal of The Electrochemical Society 162.9 (2015): A1849-A1857." + "\n", + "[1] Ecker, Madeleine, et al. \"Parameterization of a physico-chemical model of a lithium-ion battery II. Model validation.\" Journal of The Electrochemical Society 162.9 (2015): A1849-A1857.\n", + "\n", + "[2] Richardson, Giles, et. al. \"Generalised single particle models for high-rate operation of graded lithium-ion electrodes: Systematic derivation and validation.\" Electrochemica Acta 339 (2020): 135862" ] }, { diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py index b5032ab881..5b0f895121 100644 --- a/pybamm/expression_tree/symbol.py +++ b/pybamm/expression_tree/symbol.py @@ -737,7 +737,7 @@ def shape(self): ) # Pick a y that won't cause RuntimeWarnings y = np.linspace(0.1, 0.9, min_y_size) - evaluated_self = self.evaluate(0, y, y) + evaluated_self = self.evaluate(0, y, y, inputs="shape test") # Return shape of evaluated object if isinstance(evaluated_self, numbers.Number): diff --git a/pybamm/input/discharge_data/Ecker_1C.csv b/pybamm/input/discharge_data/Ecker_1C.csv index 96cee89957..dbd1215881 100644 --- a/pybamm/input/discharge_data/Ecker_1C.csv +++ b/pybamm/input/discharge_data/Ecker_1C.csv @@ -1,31 +1,31 @@ -0.042382325553453,4.10984760218981 -0.286442281697757,4.06170981679401 -0.516421855756812,4.02086563524606 -0.821496800937191,3.98148017446767 -1.10779790333724,3.94063599291972 -1.39879246643237,3.91146157752832 -1.71325433300291,3.87499355828908 -1.99486197470788,3.84581914289769 -2.29054999849809,3.81226856519758 -2.59562494367847,3.79038775365404 -2.8866195067736,3.76121333826264 -3.19169445195398,3.73641508517996 -3.47799555435403,3.7247453190234 -3.77368357814424,3.70578194901899 -4.06467814123938,3.69848834517114 -4.35567270433451,3.68827729978416 -4.6560541888198,3.67514881285803 -4.96582259469527,3.66056160516233 -5.26151061848548,3.63576335207965 -5.54781172088553,3.60075405360997 -5.84819320537082,3.55407498898374 -6.13449430777087,3.51031336589665 -6.43018233156108,3.47384534665741 -6.73056381604638,3.43008372357031 -6.97931723288577,3.38194593817451 -7.2421510318104,3.30755117892646 -7.41580907752846,3.21273432890442 -7.54722597699078,3.11500003734325 -7.62232134811211,3.01872446655165 -7.69741671923343,2.89910936344693 -7.75373824757442,2.76636577341609 +20.3084233101775,4.10984760218981 +137.255118370379,4.06170981679401 +247.454888715569,4.02086563524606 +393.638257540821,3.98148017446767 +530.825726746058,3.94063599291972 +670.262170856298,3.91146157752832 +820.943489491555,3.87499355828908 +955.88198379179,3.84581914289769 +1097.56740280703,3.81226856519758 +1243.75077163228,3.79038775365404 +1383.18721574253,3.76121333826264 +1529.37058456778,3.73641508517996 +1666.55805377301,3.7247453190234 +1808.24347278826,3.70578194901899 +1947.6799168985,3.69848834517114 +2087.11636100874,3.68827729978416 +2231.05075492899,3.67514881285803 +2379.48309865925,3.66056160516233 +2521.16851767449,3.63576335207965 +2658.35598687973,3.60075405360997 +2802.29038079997,3.55407498898374 +2939.47785000521,3.51031336589665 +3081.16326902045,3.47384534665741 +3225.0976629407,3.43008372357031 +3344.29333290591,3.38194593817451 +3470.23592758612,3.30755117892646 +3553.44799907126,3.21273432890442 +3616.41929641137,3.11500003734325 +3652.40289489144,3.01872446655165 +3688.3864933715,2.89910936344693 +3715.37419223154,2.76636577341609 diff --git a/pybamm/input/discharge_data/Ecker_5C.csv b/pybamm/input/discharge_data/Ecker_5C.csv index 87e7c364f5..621281ddff 100644 --- a/pybamm/input/discharge_data/Ecker_5C.csv +++ b/pybamm/input/discharge_data/Ecker_5C.csv @@ -1,33 +1,33 @@ --0.005659174238801,4.00772806063782 -0.022016756115274,3.9630519063978 -0.099509361106684,3.90466605332464 -0.232353826806242,3.85148732048119 -0.424855297935696,3.80418814691602 -0.636268654807101,3.75781321732131 -0.86305752854188,3.72697310974463 -1.08215864384497,3.69851401469471 -1.32047915522728,3.67723564231538 -1.55879966660959,3.65476333775538 -1.82787121171865,3.63351237271449 -2.13538154898615,3.60393814158189 -2.89262575450736,3.55566182788309 -3.56146073806417,3.51924603989306 -4.08422831141892,3.4910575923102 -4.26873451377942,3.47331305363065 -4.49167950829835,3.45679670630482 -4.7415316573282,3.43791647603876 -4.98369604792635,3.41067186867334 -5.21432880087697,3.38938664445941 -5.533370775792,3.35265909799465 -5.75631577031094,3.32539736104271 -5.99079240247741,3.29098230875862 -6.20220575934881,3.25535276879001 -6.41746299543606,3.214950926016 -6.61734471465993,3.17572931175344 -6.78647540015705,3.11737737526159 -6.93254281035912,3.06736240853066 -7.07092246212949,3.00540126815833 -7.15164392566221,2.94219480684576 -7.24005314762661,2.87660733300644 -7.31693073194349,2.79787632742773 -7.36690116174946,2.71195774734383 +0,4.00772806063782 +2.10996257174233,3.9630519063978 +9.53641973294314,3.90466605332464 +22.2674891521443,3.85148732048119 +40.7157517827462,3.80418814691602 +60.9764236014983,3.75781321732131 +82.7105988252504,3.72697310974463 +103.708022346502,3.69851401469471 +126.547325124005,3.67723564231538 +149.386627901507,3.65476333775538 +175.172937489009,3.63351237271449 +204.643005589013,3.60393814158189 +277.21304828527,3.55566182788309 +341.310446402776,3.51924603989306 +391.409562172782,3.4910575923102 +409.091603032784,3.47331305363065 +430.457402405285,3.45679670630482 +454.401832736538,3.43791647603876 +477.60951136529,3.41067186867334 +499.712062440292,3.38938664445941 +530.287258094045,3.35265909799465 +551.653057466548,3.32539736104271 +574.123984392801,3.29098230875862 +594.384656211552,3.25535276879001 +615.013703881554,3.214950926016 +634.169248146556,3.17572931175344 +650.377785601557,3.11737737526159 +664.37606794906,3.06736240853066 +677.63759859406,3.00540126815833 +685.373491470311,2.94219480684576 +693.846136049062,2.87660733300644 +701.213653074063,2.79787632742773 +706.002539140314,2.71195774734383 diff --git a/pybamm/input/parameters/lithium-ion/cells/kokam_Ecker2015/parameters.csv b/pybamm/input/parameters/lithium-ion/cells/kokam_Ecker2015/parameters.csv index 6d1c423d59..d6f636959e 100644 --- a/pybamm/input/parameters/lithium-ion/cells/kokam_Ecker2015/parameters.csv +++ b/pybamm/input/parameters/lithium-ion/cells/kokam_Ecker2015/parameters.csv @@ -11,5 +11,5 @@ Electrode height [m],1.01E-01,, Electrode width [m],8.50E-02,, ,,, # Electrical,,, -Cell capacity [A.h], 0.15625,, +Cell capacity [A.h], 0.15625, 7.5/48 (parameter set for a single layer cell), Typical current [A], 0.15652,, diff --git a/pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Ecker2015/electrolyte_diffusivity_Ecker2015.py b/pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Ecker2015/electrolyte_diffusivity_Ecker2015.py index 0671a99e18..702404696d 100644 --- a/pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Ecker2015/electrolyte_diffusivity_Ecker2015.py +++ b/pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Ecker2015/electrolyte_diffusivity_Ecker2015.py @@ -1,5 +1,5 @@ import pybamm -from pybamm import exp +from scipy import constants def electrolyte_diffusivity_Ecker2015(c_e, T, T_inf, E_D_e, R_g): @@ -37,29 +37,16 @@ def electrolyte_diffusivity_Ecker2015(c_e, T, T_inf, E_D_e, R_g): Solid diffusivity """ - # Depends on electrolyte conductivity. Have just hard coded in now for - # convinience, but should be able to call the conductivity directly - - # mol/m^3 to mol/l - cm = 1e-3 * c_e - - # value at T = 296K - sigma_e_296 = 0.2667 * cm ** 3 - 1.2983 * cm ** 2 + 1.7919 * cm + 0.1726 - - # add temperature dependence - C = 296 * exp(E_D_e / (R_g * 296)) - sigma_e = C * sigma_e_296 * exp(-E_D_e / (R_g * T)) / T - - ## Depends on the electrolyte conductivity - # E_k_e = pybamm.Parameter("Electrolyte conductivity activation energy [J.mol-1]") - # sigma_e = pybamm.FunctionParameter( - # "Electrolyte conductivity [S.m-1]", c_e, T, T_inf, E_k_e, R_g - # ) + # The diffusivity epends on the electrolyte conductivity + E_k_e = pybamm.Parameter("Electrolyte conductivity activation energy [J.mol-1]") + sigma_e = pybamm.FunctionParameter( + "Electrolyte conductivity [S.m-1]", c_e, T, T_inf, E_k_e, R_g + ) # constants - k_b = 1.38 * 1e-23 - F = 96487 - q_e = 1.602 * 1e-19 + k_b = constants.physical_constants["Boltzmann constant"][0] + F = constants.physical_constants["Faraday constant"][0] + q_e = constants.physical_constants["electron volt"][0] D_c_e = (k_b / (F * q_e)) * sigma_e * T / c_e diff --git a/pybamm/parameters_cli.py b/pybamm/parameters_cli.py index 39c6290fc3..8f2197eaef 100644 --- a/pybamm/parameters_cli.py +++ b/pybamm/parameters_cli.py @@ -33,7 +33,7 @@ def get_parser(description): """ parser = argparse.ArgumentParser(description=description) parser.add_argument( - "parameter_dir", type=str, help="Name of the parameter directory", + "parameter_dir", type=str, help="Name of the parameter directory" ) parser.add_argument("battery_type", choices=["lithium-ion", "lead-acid"]) parser.add_argument( @@ -150,6 +150,7 @@ def list_parameters(arguments=None): Available package parameters: * graphite_Chen2020 * graphite_mcmb2528_Marquis2019 + * graphit_Ecker2015 * graphite_Kim2011 Available local parameters: """ diff --git a/pybamm/simulation.py b/pybamm/simulation.py index 1ae9621a4c..6e85bf8566 100644 --- a/pybamm/simulation.py +++ b/pybamm/simulation.py @@ -387,6 +387,8 @@ def solve( # to correspond to a single discharge elif t_eval is None: C_rate = self._parameter_values["C-rate"] + if isinstance(C_rate, pybamm.InputParameter): + C_rate = inputs["C-rate"] try: t_end = 3600 / C_rate except TypeError: