tutorial/pyfdam_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial\n",
    "\n",
    "In this Jupyter notebook, we will do a walkthrough of using PyFDAM to fit a mock impedance data and using it to validate a galvanostatic discharge data. Before you run this tutorial, make sure you've installed PyFDAM!\n",
    "\n",
    "First, let's import the library. We will also import `numpy` as it will be useful later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyfdam\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll define the **model** and its **parameters**. In PyFDAM, a model is a formula $ Z(s) $ for the impedance data, where $ s $ is the frequency argument. In Python, the model is written as a function followed by parameters as an additional argument:\n",
    "\n",
    "```python\n",
    "def Zmodel(s,parameters):\n",
    "    # Define your model which takes in 'parameters' as argument\n",
    "```\n",
    "\n",
    "The argument ` parameters ` cannot be defined/managed arbitrarily! PyFDAM provides an object called ` Params ` that will manage the parameters, how they're labeled, and whether they should be designated as fitting vs. fixed parameters. \n",
    "\n",
    "\n",
    "For this tutorial, our model is defined as a symmetric cell with a porous electrode of thickness $ L $ and specific interfacial area $ a$. The contact surface area of the cell is denoted as $A$. Also, the electrolyte's conductivity is denoted by $ \\kappa $ and the current collector's conductivity is denoted by $ \\sigma $. The cell is connected to an external resistor $ R_\\mathrm{ext} $, which model's separator resistance + other external resistance contributions, as well as an inductance element $ L_\\mathrm{ind} $ in series.  Finally, double layer capacitance at the electrode/electrolyte interface is modelled as a simple capacitance $ C_\\mathrm{dl} $.\n",
    "\n",
    "In total we have 8 parameters ($ L, a, A, \\kappa, \\sigma, R_\\mathrm{ext}, L_\\mathrm{ind}, C_\\mathrm{dl} $). For this tutorial, we assume that $ L, a, A, \\sigma $ will be known and is therefore fixed while the rest will be treated as fitting parameters. To do this, we will create the `Params` object and add a list of parameters to it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Step 1 -  Define the Parameters List\n",
    "myparams = pyfdam.Params()\n",
    "\n",
    "#Step 2 - Write the fixed parameters\n",
    "#Electrode thickness (cm)\n",
    "myparams.add('L',0.0242,vary=False)\n",
    "#Conducitivity of current collector (S/cm)\n",
    "myparams.add('sig',8,vary=False)\n",
    "#The interfacial area (cm^2/cm^3)\n",
    "myparams.add('a',6*10**6,vary=False)\n",
    "#The cell area (cm^2).\n",
    "myparams.add('A',71.43,vary=False)\n",
    "\n",
    "#Step 3 - Write the fitting parameters\n",
    "#Seperator/solution resistance, (ohm)\n",
    "myparams.add('Rext',1.0,vary=True,min=0)\n",
    "#External inductance (F/cm2)\n",
    "myparams.add('Lind',1e-5,vary=True,min=0)\n",
    "#double layer capacitance (F/cm2)\n",
    "myparams.add('Cdl',1e-5,vary=True,min=0)\n",
    "#Conductivity of electrolyte (S/cm)\n",
    "myparams.add('kap',1.0,vary=True,min=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that each parameter can be accessed via `.value` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "0.0242\n"
     ]
    }
   ],
   "source": [
    "#As an example, we can print the conducitivity of the electrolyte\n",
    "print(myparams['kap'].value)\n",
    "\n",
    "#Or print the electrode thickness \n",
    "print(myparams['L'].value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the parameters have been defined, we can actually write the model's function. Let $ Z(s) $ be the cell's total impedance, then the equation for the $ Z(s) $ can be written as\n",
    "\n",
    "\\begin{gather}\n",
    "Z(s) = \\frac{2 L}{(\\kappa + \\sigma)} \\left[1+\\frac{2+\\left(\\frac{\\kappa}{\\sigma}+\\frac{\\sigma}{\\kappa}\\right)\\cosh \\nu}{\\nu \\sinh \\nu}\\right] +R_\\mathrm{ext} +s \\cdot L_\\mathrm{ind}\n",
    "\\\\\n",
    "\\nu = \\sqrt{a \\left(s \\cdot C_\\mathrm{dl} \\right) L^2 \\left(\\frac{1}{\\kappa}+\\frac{1}{\\sigma}\\right)}\n",
    "\\end{gather}\n",
    "where the first term is the porous electrode contribution.\n",
    "\n",
    "The corresponding Python function is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Zmodel(s, params):\n",
    "    #Access the raw numbers from the Params object\n",
    "    Rext = params['Rext'].value\n",
    "    Cdl = params['Cdl'].value\n",
    "    Lind = params['Lind'].value\n",
    "    K = params['kap'].value\n",
    "    sig = params['sig'].value\n",
    "    L = params['L'].value\n",
    "    a = params['a'].value\n",
    "    A = params['A'].value\n",
    "    \n",
    "    res = Cdl*s\n",
    "    nu = (a*res*L**2*(1/sig+1/K))**0.5\n",
    "    imp = (1+(2+(sig/K+K/sig)*np.cosh(nu))/(nu*np.sinh(nu)))\n",
    "\n",
    "    return (2*L*imp/(A*(K+sig))+Rext+Lind*(s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we've defined the parameters and model, we can now use them to do impedance fitting. To do this we construct an `ImpedanceFitter` object that will take our defined model, parameters, and impedance data filename as an input.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "MyFit = pyfdam.ImpedanceFitter(model=Zmodel,parameters=myparams,namedata=\"EIS_Mock_Data.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the datafile must be formatted with three columns and no headers! The first column is frequency, the second column is the real part of the impedance, and the third column is the imaginary part of the impedance. In this tutorial, the experimental data is `EIS_Mock_Data.txt` so check that file for reference!. \n",
    "\n",
    "Finally, we can run the fitting process using the `.run()` method. The fitting run is done by selecting a suitable numerical method as an initial error minimization step. This is then followed by some number of cycles, denoted by `ncycles`, of Leveberg-Marquadt algorithm to further converge the parameters and obtain their uncertainty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Variables]]\n",
      "    L:     0.0242 (fixed)\n",
      "    sig:   8 (fixed)\n",
      "    a:     6000000 (fixed)\n",
      "    A:     71.43 (fixed)\n",
      "    Rext:  0.01531600 +/- 1.1009e-04 (0.72%) (init = 0.02352427)\n",
      "    Lind:  4.8185e-06 +/- 1.6489e-07 (3.42%) (init = 3.985611e-06)\n",
      "    Cdl:   7.6876e-06 +/- 7.8261e-09 (0.10%) (init = 7.660804e-06)\n",
      "    kap:   0.01425053 +/- 1.6015e-04 (1.12%) (init = 3.835463)\n",
      "[[Correlations]] (unreported correlations are < 0.100)\n",
      "    C(Rext, kap) =  0.690\n"
     ]
    }
   ],
   "source": [
    "MyFit.run(method='nelder',ncycles=1) #In here, we've chosen Nelder-Mead ad the initial method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A list of numerical methods for the initial fitting can be found in [the LMFIT package's documentation](https://lmfit.github.io/lmfit-py/fitting.html).\n",
    "\n",
    "Next, we can take a look at the result of the fitting. To do this, we need to grab the output result as well as the experimental impedance data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get both the frequency (expfreq) and impedance data (expdata) of from the loaded textfile!\n",
    "expfreq, expdata = MyFit.get_impedanceexp()\n",
    "\n",
    "#Get the impedance output from the model using .get_impedancemodel() method. It takes a list of frequencies as an argument\n",
    "modeldata = MyFit.get_impedancemodel(expfreq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Afterwards, we can plot the data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 0.07)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT0AAAFzCAYAAABMyd6BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAk1klEQVR4nO3dfXBd9X3n8ff3PulKsh78iIytYEMcYmMoGEMItHlOB5MAnW2nTdpN2my3lC2h7XZnu8lkdic73ZnttjPdbCgblibpljbbNE1pi1M2hKak7YaBYAIBGWEwjkHGli3b6Mp6uLpP3/3jHIFwhC18r8/vyufzmtFY9zzofpHER7/fOb/f75i7IyKSFpnQBYiIJEmhJyKpotATkVRR6IlIqij0RCRVFHoikiq50AWciVWrVvmGDRtClyFvxvRxwKFrJc8cnKC/K8/5/Z2hq5Il5PHHHz/q7qub/TpLMvQ2bNjArl27Qpchb8bQvTB9DK7+FS777AP8i23r+exNl4SuSpYQM3uxFV9H3VtJRq4I1RkANBxeQlLoSTLyRajNvvrSLGAtkmoKPUlGrgiNGtSraupJUAo9SUauGP1bKwNgqKknYSj0JBlzoVcth61DUk+hJ8nIv9bSU+9WQlLoSTJO7t6qdyuBKPQkGfNCT2s4SkgKPUnGSdf01NCTUBR6koyTurcioSj0JBmZDOQKupEhwSn0JDm5Tt3IkOAUepKcXAdUy+g+hoSk0JPk5Oe39NTUkzAUepKcXIduZEhwCj1JTnxNz3UrQwJS6Ely4mt6oHF6Eo5CT5KT74RGjYzXQ1ciKabQk+TkOgAoUlFTT4JR6ElyctGDgApUAxciaabQk+TELb0OKoELkTRT6Ely8lFLr4OKVk6WYBR6kpx40YGiq3sr4Sj0JDlx6BWsorm3EoxCT5Iz19LTNT0JSKEnyYmXlyoo9CQghZ4kK1ekSFW3MSQYhZ4kK1fUkBUJSqEnyYpDTzcyJJRc6AIkRUaH4Llv8nP2DFMjDRj9BAxsDV2VpExLWnpmdr2Z7TGzvWb2qQX2m5l9Pt7/lJltm7ev38y+bmbPmtmwmb2zFTVJmxkdgofvgNosk1akWDsRvR4dCl2ZpEzToWdmWeBOYAewBfiomW056bAdwKb44xbgC/P2/Q/gm+7+duDHgOFma5I2NLwTiv1Q7CNLg0quJ3o9vDN0ZZIyrWjpXQ3sdfd97l4BvgrcfNIxNwP3eOQRoN/M1ppZL/Au4EsA7l5x9/EW1CTtpjQCxV7IZMnQAPfodWkkdGWSMq0IvXXA/N/cA/G2xRxzITAG/LGZPWFmXzSz7hbUJO2mbxDKE2BZgCj4yhPRdpEEtSL0FroPd/J64G90TA7YBnzB3a8ApoAfuSYIYGa3mNkuM9s1NjbWTL0SwuYboTwO1WnAKdZK0evNNwYuTNKmFaF3AJj/53o9cHCRxxwADrj7o/H2rxOF4I9w97vdfbu7b1+9enULypZEDWyFa2+HYj89XqaS7Y5e6+6tJKwVofcYsMnMNppZAfgIcN9Jx9wHfDy+i3sNUHL3Q+4+CoyY2cXxce8HnmlBTdKOBrbCNf+GB/1Knlj38wo8CaLpcXruXjOzTwIPAFngy+6+28xujfffBdwP3ADsBaaBT8z7ErcDX4kDc99J++Qc49k8ANmGlpeSMFoyONnd7ycKtvnb7pr3uQO3vcG5TwLbW1GHLAGZ6FcupyXjJRBNQ5NEeaYAQMZrgSuRtFLoSaI8G4VetqHQkzAUepIoV/dWAlPoSbLU0pPAFHqSKM/kcDeyejiQBKLQk0Q5RpWcQk+CUehJohynQo6s7t5KIAo9SZQ7aulJUAo9SVyVnMbpSTAKPUlclaxaehKMQk8S5Q5Vz2nIigSj0JNEOa5rehKUQk8Sp+6thKTQk0S5Q4UcGXVvJRCFniTK0ZAVCUuhJ4lyd6qeI+MNaNRDlyMppNCTxFXJYQbUK6FLkRRS6EminOiaHgB1dXEleQo9SdTcNDRAoSdBKPQkWQ41ogd+q3srISj0JFFzq6wYKPQkCIWeJK7qOTDUvZUgFHqSqNdf01NLT5Kn0JNEzQ1OVvdWQlHoSaLcPR6yYureShAKPUmcBidLSAo9SZQDDTI0LAMNtfQkeQo9SZR7/K/l1b2VIBR6kignSr1GJq/urQSh0JMgFHoSikJPkhV3b6PQU/dWkpcLXYCkiwMX24usG/tnmJiBYy/A5hthYGvo0iQl1NKTRGWPPMMt2fvJNGahoxdmxuHhO2B0KHRpkhIKPUlUx96/o0QX9dwyoAGd/VDsh+GdgSuTtFDoSaKyEwc4QRevXtwDKPZCaSRYTZIuCj1JVLVnHT1Mk/E6ZPLRxvIE9A2GLUxSQ6EniZq56Ab6mCZfOwGZXHRNrzwe3cwQSYBCTxJVXbWFu+s30LA8zJaia3rX3q67t5IYDVmRRLnDHn8Lx1deydqr3wcXvS90SZIyaulJohru5KlH1/RynaHLkRRS6EmiHOhkNlpaKq/Qk+Qp9CRRDXc6iOfc5rvCFiOppNCTRLlDp1UwM7X0JAiFniTK3elkNnqhlp4EoNCTREXX9CrRg4HU0pMAFHqSqKh7OwtmkCuGLkdSSKEniWq4U6SCZzsgo18/SZ5+6yRR7tGQlbrG6EkgCj1J1FxLr6GurQSi0JPEdVoFz+nOrYSh0JNERS29WRpZtfQkDIWeJMobTicVXN1bCUShJ4nyeoWc1WnoRoYEotCTRFltBkB3byUYhZ4kqxZNQXOFngSi0JNEZarTgEJPwlHoSaKsHnVvFXoSikJPklWJQ0+LDUggCj1JVKY+g7vhuY7QpUhKKfQkUVadYZY8ZtnQpUhKKfQkOaNDrH7hr/jJzGOs/8HnYHQodEWSQgo9ScboEDx8B9nKCY7RS74yAQ/foeCTxCn0JBnDO6HYTyOTo0aORkcfFPuj7SIJUuhJMkojUOwFdxoYGNHr0kjoyiRlWhJ6Zna9me0xs71m9qkF9puZfT7e/5SZbTtpf9bMnjCzb7SiHmlDfYNQngAMw6Nt5Ylou0iCmg49i27D3QnsALYAHzWzLScdtgPYFH/cAnzhpP2/AQw3W4u0sc03Qnkcq5cxGuRmS1Aej7aLJKgVLb2rgb3uvs/dK8BXgZtPOuZm4B6PPAL0m9laADNbD3wI+GILapF2NbAVrr2dWq6LFUxSL/bDtbdH20US1IrQWwfMvzBzIN622GM+B/w20DjVm5jZLWa2y8x2jY2NNVWwBDKwlSPrfpKvN97N0St/S4EnQbQi9GyBbb6YY8zsw8ARd3/8dG/i7ne7+3Z337569eozqVPaQMOMDK4HoUkwrfjVOwDMvxq9Hji4yGOuA24ys/1E3eL3mdmftaAmaVMNz5ChQcYW+jsocva1IvQeAzaZ2UYzKwAfAe476Zj7gI/Hd3GvAUrufsjdP+3u6919Q3zeP7j7v2xBTdKm6pYhS4NsRqEnYeSa/QLuXjOzTwIPAFngy+6+28xujfffBdwP3ADsBaaBTzT7vrI01dXSk8CaDj0Ad7+fKNjmb7tr3ucO3Haar/Ed4DutqEfal6OWnoSly8mSqDpzLb3QlUhaKfQkUXWMrLq3EpBCTxLVIEMWV/dWglHoSaLqnsXU0pOAFHqSqLrF3Vv95kkg+tWTRDXcyFoDLRYvoSj0JFH1+Fcu8yMzFUWSodCTRNU9Dj075foSImeNQk8SVY/XnsieelEdkbNGoSeJqsdX87Lq3kogCj1JVM2jlp5RD1yJpJVCTxLV8LnurVp6EoZCTxL1WvdW1/QkDIWeJKoW38gwV/dWwlDoSaLqnsEAc7X0JAyFniSqhmFmoNCTQBR6kqi6ZzAD1L2VQFqycrLIYvVPv8hFPAMPfAbWbIke9q1HQUqC1NKT5IwOcfnYfRStAsvWwMw4PHwHjA6FrkxSRKEnyRneyUy2m4oVAIPOfij2w/DOwIVJmij0JDmlEWboxOY/+73YC6WRcDVJ6ij0JDl9g3TUJqPIy8Qr6pUnoG/wVGeJtJRCT5Kz+UaK9Qk6qIJlomt65fHoZoZIQhR6kpyBrXxv2QeoZgoweSS6pnft7bp7K4nSkBVJ1NH8Gp7OXcoVN/0BFLpDlyMppJaeJKtWIZMxyHaErkRSSqEnibJ6BctkIatOhoSh0JNEZeoV6pYPXYakmEJPEmX1Werq2kpACj1JVKZRxTNq6Uk4Cj1JVKYxSyNTCF2GpJhCTxKVaVRpZBV6Eo5CTxKVbVTwjK7pSTgKPUlUtlHB1dKTgDRYSpIzOsT26uN0je+Ghya0gKgEoZaeJGN0CL77eTq8zGR+pRYQlWAUepKM4Z3Q0cMsBcjktYCoBKPQk2SURqDQTQOwTPxrpwVEJQCFniSjbxBmJ3D36BGQoAVEJQiFniRj841QLlGgEq2crAVEJRCFniRjYCuN7b/CTCNPT+MVLSAqwWjIiiRmdsXbeNS3sPbCn2XTe3eELkdSSi09ScxspQpARy4buBJJM4WeJKZcrQHQkVfoSTgKPUlMuVoHFHoSlkJPEjPX0svndClZwlHoSWJm45ZeUS09CUihJ4l5rXurlp6Eo9CTxMzO3cjI6ddOwtFvnyRmVi09aQMKPUlM4fhzvMOeYc33fhce+q9aVkqCUOhJMkaHuGD/X1C0CvQMaD09CUahJ8kY3kk508UsBXLZnNbTk2AUepKM0gjVRjRUJVfojLZpPT0JQKEnyegbpFg5Sp0smUIx2qb19CQAhZ4kY/ONFCvj0arJ7lpPT4JR6Eky1mxhf8cmSrYcJl7WenoSjAZMSTJmjvOK9/DN7vew46d+KXQ1kmJq6UkyJg9Trtapd68JXYmknEJPkjF5mOmak+leFboSSTmFniTjxGEOVnvo6y6GrkRSTqEnZ587PnmYlyo99HXlQ1cjKafQk7Nv9gS18hSH6r30dxZCVyMp15LQM7PrzWyPme01s08tsN/M7PPx/qfMbFu8fdDMHjKzYTPbbWa/0Yp6pM3ENzHGvJ9+tfQksKZDz8yywJ3ADmAL8FEz23LSYTuATfHHLcAX4u014N+5+2bgGuC2Bc6VpW7yMDM1Z4x+liv0JLBWtPSuBva6+z53rwBfBW4+6ZibgXs88gjQb2Zr3f2Qu38fwN1PAMPAuhbUJO1k8jCTmWVUydGn7q0E1orQWwfMnzV+gB8NrtMeY2YbgCuAR1tQk7STE4cpZVYAqHsrwbUi9GyBbf5mjjGzZcBfAb/p7hMLvonZLWa2y8x2jY2NnXGxkrDqDJRLHGU5oNCT8FoxDe0AMH+pjPXAwcUeY2Z5osD7irvf+0Zv4u53A3cDbN++/eRQlXY0OgTf/1N48busyf6Qi+1y3b2V4FrR0nsM2GRmG82sAHwEuO+kY+4DPh7fxb0GKLn7ITMz4EvAsLv/QQtqkXYxOhStjDw5Ch09TJdn+WTHN+k8Phy6Mkm5pkPP3WvAJ4EHiG5EfM3dd5vZrWZ2a3zY/cA+YC/wR8CvxduvAz4GvM/Mnow/bmi2JmkDwzujlZHNINfBoeoyvKNPKyVLcC1ZZcXd7ycKtvnb7pr3uQO3LXDe/2Ph632y1JVGoHcdvPJDKHRTKlc5b1mfVkqW4DQjQ86OvkEol6A6TSPXyYlyjTWFilZKluAUenJ2bL4RpsagMs10I8cyn2RlblorJUtwCj05Owa2wqU/C/kiM1PjlLyb0o/dqpWSJTiFnpwdo0Pw7N9BeYIjufV8s3EVK9+6LXRVIgo9OQvmhqtMHYGulcxMTfCruftZW34hdGUiCj05C+YPVyl0c6RapF7oJbfnG6ErE1HoyVlQGoke5F2rQK6DY5MV8t39Gq4ibUFPQ5PW6xuEo3th/EXqR5/jrdPLWN5zAfRdEboyEbX05CxYvRle3gXVWco1o9em2VB+JtouEphCT1pvbBjWXwXFHuqVaSbpJn/BO6LtIoEp9KT1SiOwYiOcfzkvspaXeraRXXmhrulJW1DoSev1DUJ5gmq9zmSlzrrlnVCe0BQ0aQsKPWm9zTdCeZzxV47TcLiguwLlcU1Bk7ag0JPWG9gK197OoekMPZRZvXoArr1dU9CkLWjIipwV0yvezhdfuZz3nHcRl33gM6HLEXmVWnpyVvz1Ey8zW6myfeOq0KWIvI5CT1rO3bnn4Re5YHkHg6t6Qpcj8joKPWm5R/YdZ8/hE7z3bSuwjK6gSHtR6EnL/fF3f8jyrjxXDvZCJhu6HJHX0Z9haalv7R7lW88c5r9cYxT23AfTx2HiUDRcRXdvpQ2opSctc3Rylk/f+zQ3rDnGR+t/A5Up6FwBM+PR+nqjQ6FLFFHoSWu4O5++92lOzNb4nbc+R7ZRBSzq3nb2R+vr6fGP0gbUvZUzMzoUhVhpBPoG+VZ9O7ufOcidV9dYefAfobAMCt3QuzY6vtirubfSFhR68ubNLQdf7Ifu1Yy98CQDL/0ln11+Me9feSU0NkGuE/rW8epjjTX3VtqEQk8WZ37L7pX9UOiD2gylo6PsOzxBPr+Md71tFZnrfh2OPheF4kwpauGVJ6K5t9s+Fvq/QkShJ4sw17Lr6IVMDo4+D95gong+3znaS6X4dj585YUUZ0YhV3h17u387i/bPqa7t9IWFHpyesM7IV+EY3uhPguFZUxXG+w5VqFUXM/PXLmezvrk67uvA1sVctKWFHryxua6tE/9BTjQtQLWb2dP/XwaLz3C6twUW7adT1d9Ut1XWTIUerKwV7u0fWBZmC3RyOZ4/IVRvns4y9W9m9i0skyufFjdV1lSFHqysLln184ch45l1L3BgYka9fIe3rn2Mq5au4aM1siTJUihJwsrjURj7SZeZrJnI/8wNsOayn4290zRt+VCTSuTJUuhJwvrXQcvPszxWeMvXzScHrZveQd9562F9346dHUiZ0zT0GRBvvwCjh09zEMvw7JCnp//sT7OL8zqORey5Kmll3bzBx1nOwColyf44f59PDR+IQOrz+OmwRr5/lWw+RPq0sqSp9BLs1fv0PaA5WD/P1N3Z3iii/LsNNcNrOTtP30LmfMvDV2pSMso9NJseCdkCzD2LLyyn1q9wZHJKn2Ncbo2vIsL166BPd8AhZ6cQxR6aTU6BM9+A6aPQbaDSqXKyEyODBkGevJ0btwEuFZGkXOOQi+N5rq1tTI0nKlskZmZ4zTo5fzVA3R2dYNZtGCAVkaRc4zu3qbR8E6oV6CwjBmyPP9KlRJ9bOyYpLMxBSs3Rasdl8d1t1bOOWrppcX8u7QvPQL5bkaz5/HQxAouKRzhLcuL5LwIG34iWlSgc0BTy+ScpNBLg/mLfk4dh8nD1Go1jtQm6O/azMXbd5CvnoiWddfAYznHKfTSYG4e7fQxOLiLGetmpj7J+sw4l/TsJzveDdmcVkmRVNA1vTQojUQrGE8cYLqRZ6zcoJRbTX9XniwNmDwULfqprqykgFp6adA3CDPjTJRe4fAU9OWMlSu6yHStgwuug4mXFXiSGgq9c93oEEweYXrP31ObOMbKTJa+/vPImsGqi/XAHkkddW/PZfENjIOTNR6b6MMyWfqzs2QzOVh/VTQbQ8NSJGXU0juXDe/k2MQ0Ez98mK2ZCXpWrSfTdz7UZqBR1bAUSSWF3jmstO/7jL30LH02TV//SnKFTjhxEPo3wE/9z9DliQSh0Fvq5g867ht8dUXj5146RG3kOVbYFKt6Osnl89ETzWqzMFsKXbVIMAq9pWxu0HG9BpMH4cD34OmvUVr7Ezw0XOYnqbCqK0euozsap1ctgzeg2Be6cpFgFHpL2fDOKPAOPQmNOlSnacxOUDj+57yNS1mx6VpyHRZ1acsT0Vi9FRfCyotCVy4SjEJvKSuNQOklmHkFamUatTJVNyrkuXp1jW4mYdbhvEujwCtP6G6tpJ6GrCxlfYNw4jDgeK1CzR3zBsusRnd9CpZviB7w09kfDUDu7NfMC0k9tfSWkrmbFqNPQbkEGEwdwclR8wZ1z1DIZsjmi9FKKcXeKOy0iIDIqxR6S8XcTYtGHY7vB8tAfRa3LPXqJLiTzxbJFpdFx+eKmm0hsgB1b5eKuZVSDg9BZRKqU1CZZKYOhxv9VAoryBWKkMlD1yro6NX1O5EFqKW3VJRGoutzk0ei2RTAdCPPsWqBfFc/3T15eMs7YfRpmDkOa6+Fq35Z1+9ETqLQWyr6BuHo3ijwGjXKuT4my5P05LP0rrsQ6vHUsk0ffHWAsoj8KIXeUrH5Rrj3Fij205gYpVE9Sp9BvtBN5pV98HN/qqATWQSFXjsbHYLHvgQHHotel0s0skVq3iDjkCt0kcnFd2pFZFEUeu1qdAi+/Ttw/AUo9IAB1WlmTxzneGMZy7u7yHYWobg8On54p1p6Ioug0GtXwzthauy1Ne9mJ6jPTlFozLIyU6PYuRKyHdGza9dfpYdyiyySQq9dlUaiFVGmxqAyiVuGqUaeDqtRyAAYdC2PVj/OFqK18UTktDROr131DUKuAzCwDLVanaxXyFqWjGWi8XgXXKfVj0XepJaEnpldb2Z7zGyvmX1qgf1mZp+P9z9lZtsWe25qbb4RuldDdRpvGNaYpWBVsoVO6FkLJw5F1/00n1bkTWk69MwsC9wJ7AC2AB81sy0nHbYD2BR/3AJ84U2cm04DW+H9/xFWXkS9UaXmBlbE+tZB14poeai1l0XzahV4IovWipbe1cBed9/n7hXgq8DNJx1zM3CPRx4B+s1s7SLPTa+BrXDznYznVvASA2SXrYRMLrp5MXCpbl6InIFWhN46YP7/fQfibYs5ZjHnAmBmt5jZLjPbNTY21nTRS0Vt9Ra+XbmE5XknQx06emDwHdGCAlpMQORNa8XdW1tgmy/ymMWcG210vxu4G2D79u0LHnNOiZeReuXlvdRqNYpdBTh/C5y39bXFQLd9LHSVIktOK0LvADC/ybEeOLjIYwqLODd95paRKvYxWqqxLfM83bXxaEkpJ7qWp0c3ipyRVnRvHwM2mdlGMysAHwHuO+mY+4CPx3dxrwFK7n5okeemz9wyUqWXGTz2HfqyVTIdvZDvjBYG1YICImes6dBz9xrwSeABYBj4mrvvNrNbzezW+LD7gX3AXuCPgF871bnN1rTklUag2Iu/8kPKDaPbqtFA5eP74Ojz0XxcETkjLZmR4e73EwXb/G13zfvcgdsWe27q9Q3CzDi1WhW8ThdVoBMK3YDDvoeiLrBaeyJvmmZktKPNN0J5nHKtQb9Pkm1UoTIF9WjxUDpXRF1gEXnTFHrtaGArvG0H9WqFnHl0i9syMFuCmXGN0RNpgkKvXY0N82z3O3iZ1VAoRqGXyUU3MjRGT+SMaZWVdlUa4Ui1g1L+Ei5Y9hKYQXUWqjMaoyfSBLX02lXfIJQnmO5aDxvf/dr0s+5VWmBApAlq6bWj0SGYPMIVs49St+Xg10Df+uju7Y7fU+CJNEEtvXYTz8bwTJ5/ql9CIZuB5x+EI3ugVonu2o4Oha5SZMlS6LWbeDZGZXaaPqaoda6CTDaapdz/luju7cN3KPhEzpBCr93EszEq5WlWWYlVE7th+li0bPyxvVCvRFPUNE5P5Iwo9NpNfANjonMdZc/TUZuMhqt09EVzbw98L7qhoXF6ImdEoddu4tkY9alxilTwbAEa1WgdvUJ39AS00ac1Tk/kDCn02s3AVrj2diYzy1jDOLW+DVHgVafBPfqYOa4HAYmcIYVeOxrYymMbfpW/bVyHD74Dzrs0ejLa7IlokPKF79WwFZEzpNBrUyfKVb7ZuIpC7UTUulv1tuiRj6s2wVW/HLo8kSVLodemTpRr7M9uJHvdr0M2B5VJPe5RpAU0I6NNTVVqdHfkooDb+BOw6mK4+PrQZYkseWrptanpSp3OfBYajWiRgUJX6JJEzgkKvTY1U6nTVchCdSq6plfoDl2SyDlB3ds2NV2pszkzAt/5Fuz7R5g+Go3R0/U8kaaopdemVk0+x8/M3huFXUcPVMuacyvSAgq9NrVt5rvM5nog1xmNzetaqTm3Ii2g0GtTK2uHqWS74cQhyOSjwcnFXs25FWmSQq9NHfRVrKsfiGZhrLwoWnSgPKE5tyJNUui1qUd9C6srL0fd264V0Tp65XHNuRVpkkKvHTUarKu/zIs922DtZTBxULMxRFpEQ1ba0cgjrGgcY/fAT3PtBz4UuhqRc4pCr12MDkV3Zo8+j08eplQfpLP3raGrEjnnKPTaQfwwIDp6YfYEXq/y48sOMpN5Cbg4dHUi5xSFXjuIHwbE7ARUp8gMXMq7z89BdhfwwdDViZxTdCOjHZRGoNgTDU/pXh090Ftj8kTOCrX02kHfYDQk5bwt4I1om8bkiZwVaum1g/hhQMyUANOYPJGzSKHXDuKHAdHZDxMva0yeyFmk7m27GNiqkBNJgFp6IpIqCj0RSRWFnoikikJPRFJFoSciqaLQE5FUUeiJSKoo9EQkVRR6IpIqCj0RSRWFnoikikJPRFJFoSciqaLQE5FUUeiJSKoo9EQkVRR6IpIqCj0RSRWFnoikikJPRFJFoSciqaLQE5FUUeiJSKoo9EQkVRR6IpIqCj0RSRWFnoikikJPRFJFoSciqdJU6JnZCjN70Myej/9d/gbHXW9me8xsr5l9at723zezZ83sKTP7azPrb6YeEZHTabal9yng2+6+Cfh2/Pp1zCwL3AnsALYAHzWzLfHuB4Gt7n4Z8Bzw6SbrERE5pWZD72bgT+LP/wT4qQWOuRrY6+773L0CfDU+D3f/lrvX4uMeAdY3WY+IyCk1G3rnufshgPjfNQscsw4Ymff6QLztZP8K+L9N1iMickq50x1gZn8PDCyw6zOLfA9bYJuf9B6fAWrAV05Rxy3ALfHLWTMbWuT7J2EVcDR0EfO0Wz3QfjWpnlNrt3oALm7FFzlt6Ln7B95on5kdNrO17n7IzNYCRxY47AAwOO/1euDgvK/xi8CHgfe7u/MG3P1u4O74nF3uvv10tSdF9Zxeu9Wkek6t3eqBqKZWfJ1mu7f3Ab8Yf/6LwN8ucMxjwCYz22hmBeAj8XmY2fXAfwBucvfpJmsRETmtZkPvd4EPmtnzwAfj15jZ+WZ2P0B8o+KTwAPAMPA1d98dn/+HQA/woJk9aWZ3NVmPiMgpnbZ7eyrufgx4/wLbDwI3zHt9P3D/Ase99Qzf+u4zPO9sUT2n1241qZ5Ta7d6oEU12Skuo4mInHM0DU1EUqWtQu+NpqvN229m9vl4/1Nmtm2x5waq6ctmdqSVw2vOtB4zGzSzh8xs2Mx2m9lvBK6naGbfM7MfxPX855D1zNufNbMnzOwbrain2ZrMbL+ZPR1f827J3csm6+k3s69bNH102MzeGaoeM7s4/r7MfUyY2W+e9g3dvS0+gCzwAnAhUAB+AGw56ZgbiAYwG3AN8Ohiz026pnjfu4BtwFAbfI/WAtviz3uIpv019T1qsh4DlsWf54FHgWtC/rzi/b8F/B/gG6F/ZvG+/cCqdvj/LN73J8C/jj8vAP2hf2bzvs4ocMHp3rOdWnpvOF1tnpuBezzyCNBv0fjAxZybdE24+z8Bx1tQR9P1uPshd/9+XNcJojvpC82MSaoed/fJ+Jh8/NHsBeamfl5mth74EPDFJutoWU1nwRnXY2a9RH/IvwTg7hV3Hw9Vz0nHvB94wd1fPN0btlPoLWa62hsds9ipbknWdDa0pB4z2wBcQdS6ClZP3JV8kmhQ+4PuHrQe4HPAbwONJutoZU0OfMvMHrdoVlLIei4ExoA/ji8BfNHMugPWM99HgD9fzBu2U+iddrraKY5ZzLlnopmazoam6zGzZcBfAb/p7hMh63H3urtfTjRL52oz2xqqHjP7MHDE3R9vsoaW1RT/e527byNapeg2M3tXwHpyRJdrvuDuVwBTLLCyUoL1RDujSQ83AX+5mDdsp9A75XS10xyzmHOTrulsaKoeM8sTBd5X3P3e0PXMibtI3wGuD1jPdcBNZrafqIv1PjP7sybrabYmPBrzirsfAf6aqDsYqp4DwIF5LfKvE4VgqHrm7AC+7+6HF/WOzVyEbOUH0V+RfcBGXrugeclJx3yI11/Q/N5iz026pnn7N9C6GxnNfI8MuAf4XJv8zFYTXwQHOoF/Bj4c+ucVH/MeWncjo5nvUTfQM+/zh4HrQ36P4p/TxfHnnwV+P/TPjOiP1CcW/Z6t+h+gRb8gNxDdVXwB+Ey87Vbg1vhzI1qQ9AXgaWD7qc5tg5r+HDgEVIn+Wv1yqHqAHyfqEjwFPBl/3BCwnsuAJ+J6hoD/FPrnNe9rvIcWhV6T36ML4xD4AbC7Vb/XTf5OXw7sin9ufwMsD1xPF3AM6Fvs+2lGhoikSjtd0xMROesUeiKSKgo9EUkVhZ6IpIpCT0RSRaEnbcPMfsnM/vAMzx08acWNuVU3/lur65SlramVk0XahbuPEI0hA8DMLiVarftzgUqSNqXQk7ZkZv8bmAHeDlwAfILo4VPvJFpa6JdOcW6R6HGit3n8XGaROereSjtbDrwP+LfATuC/A5cAl5rZ5ac47/eA77r7fWe9Qlly1NKTdrbT3d3MngYOu/vTAGa2m2hO85Mnn2BmO4APAG31zFZpH2rpSTBmdtu8mw7nL3DIbPxvY97nc69/5A+2ma0G/hfwC67nKMsbUEtPgnH3O4kmkgNgttCyaW/Kl4E73P2JZr+QnLsUenJOiB9Q82HgLWb2C/N2Peju/z5QWdKGtMqKiKSKrumJSKoo9EQkVRR6IpIqCj0RSRWFnoikikJPRFJFoSciqaLQE5FU+f9Xd/BIgSaJzwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(6,6))\n",
    "#Plotting model data\n",
    "plt.plot(modeldata.real,-modeldata.imag)\n",
    "\n",
    "#Plotting experimental data\n",
    "plt.plot(expdata.real,-expdata.imag,'o-',alpha=0.5)\n",
    "\n",
    "#Impedance plot should be equal axes\n",
    "plt.gca().set_aspect('equal', adjustable='box')\n",
    "\n",
    "#Some labelling\n",
    "plt.xlabel('Re Z')\n",
    "plt.xlabel('-Im Z')\n",
    "#Adjust the range as necessarry!\n",
    "plt.ylim([-0.02,0.07])\n",
    "plt.xlim([0,0.07])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also save the impedance data for later plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We use the list of frequencies from the experimental data, although we could use any other list of frequencies\n",
    "MyFit.save_impedance(expfreq,'EIS_Model_Data.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What's next? Well we can use this set of parameters to run a discharge experiment. The data that we will compare  with starts  at a fully charged cell at 2.4 V, which was discharging at a current of 200 A. To setup the galvanostatic discharge simulation, we need to create a `Galvanostatic` object and feed it wth out model and the fitted parameters, which we can get from the previous `ImpedanceFitter` object via the `.get_params()` method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "MyGalvanostatic = pyfdam.Galvanostatic(model=Zmodel,parameters=MyFit.get_params())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run the simulation, we can use the `.run()` method which will use the Gaver-Stehfest (GV) method to perform numerical inverse Laplace transform. The `.run()` method takes in the following arguments: \n",
    "1. `time` ($ t $) is an array of time which starts from **intial time larger than zero** (the GV method cannot handle t=0) and up to some final time. \n",
    "2. `DeltaV0` ($ \\Delta V_0 $)the initial total difference in potential across cell\n",
    "3. `Iapp` ($ I_\\mathrm{app} $) is the applied current, with (-)/(+) sign indicating discharging/charging experiments respectively,\n",
    "4. `A` ($ A $) is an optional parameter that is divided against the impedance function $ Z(s) $. This is used only if the original model is not described in the correct Ohmic units, e.g., the model impedance is multiplied by surface area so we need `A` to be the cell's area. By default, it is set to one. For our tutorial, `A` will be the surface area parameter $ A $ since the original model function $ Z(s) $ is expressed in terms of Ohm cm^2. \n",
    "5. `Nterms`($ N_\\mathrm{terms} $) is a parameter for the Gaver-Stehfest method, which should be around 6 to 20 to obtain good convergence of the time-domain profile. By default, it is set to 10.\n",
    "\n",
    "Under the hood, the code is inverting the following function in Laplace space:\n",
    "\n",
    "\\begin{equation}\n",
    "V(s) = \\frac{\\Delta V_0}{s}+\\frac{I_\\mathrm{app} Z(s)}{s \\cdot A}\n",
    "\\end{equation}\n",
    "\n",
    "by using the following series expansion:\n",
    "\n",
    "\\begin{gather}\n",
    "V(t) \\approx \\Delta V_0+\\frac{\\ln 2}{t} \\sum_{n=1}^{N_\\mathrm{terms}} c_{n} \\frac{I_\\mathrm{app} t}{A(n \\ln 2)} \\cdot Z\\left(\\frac{n \\ln 2}{t}\\right)\n",
    "\\\\\n",
    "c_{n}=(-1)^{n+\\frac{N_\\mathrm{terms}}{2}} \\sum_{k=\\frac{n+1}{2}}^{\\min (n, N_\\mathrm{terms} / 2)} \\frac{k^{N_\\mathrm{terms} / 2}(2 k) !}{(N_\\mathrm{terms} / 2-k) ! k !(k-1) !(n-k) !(2 k-n) !}\n",
    "\\end{gather}\n",
    "\n",
    "PyFDAM automates the whole series expansion computation so all we need to do is just feed in the right parameters!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load the discharge data\n",
    "expdata = np.loadtxt(\"Discharge_Mock_Data.txt\")\n",
    "#The time array is obtained from the mock data. Note that we also exclude the initial time!\n",
    "MyGalvanostatic.run(time = expdata[1:,0], V0=2.4,Iapp=-200,A=myparams['A'].value)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to the impedance run, we can get the data, which is done through the `.get_output()` method, save the data, which is done through the `.save_output` method, and plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAE9CAYAAABtDit8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA200lEQVR4nO3dd3iTZffA8e9pm7RhFqEiQwQXorJRQFy4UBRElKUsEREQBMFKUWQL1QJuX8UBDjZiRXkV9Qc4wAVSQAVUFLWFVxBoGS1N2t6/P5LUtrSkQJMnTc7nunLZJE/yHKEe73luMcaglFKqZBFWB6CUUsFOE6VSSvmgiVIppXzQRKmUUj5oolRKKR80USqllA9RVgdwomrUqGHq169vdRhKqRCzYcOGf4wxccW9V+4SZf369Vm/fr3VYSilQoyI/FHSe9r1VkopHzRRKqWUD5oolVLKh3I3RqnUqXK5XKSmpnL06FGrQ1EWiImJoW7duthstlJ/RhOlCjupqalUrlyZ+vXrIyJWh6MCyBjDvn37SE1NpUGDBqX+nHa9Vdg5evQo1atX1yQZhkSE6tWrn3BvQhOlCkuaJMPXyfzdh3SiTN6YRrvEVTRIWEG7xFUkb0yzOiSlAIiMjKRZs2b5j8TERL/eb/ny5X6/x5o1a1i3bp3P6+bOncuwYcOKfT0uLo7mzZtz3nnn0aFDh1J9X3JyMj/99NNJxVxaITtGmbwxjbHLtpDlygUgLT2Lscu2ANCleR0rQ1MKh8NBSkpKQO6Vk5ND586d6dy5s1/vs2bNGipVqsRll1120t/Ro0cPnn/+eQBWr15N165dWb16NY0aNSrxM8nJydxyyy1ceOGFJ31fX0K2RZm0cnt+kvTKcuWStHK7RRGp8ipQPZOMjAwaNmzI9u3u39FevXrxyiuvAFCpUiVGjx5NixYtuPbaa9m7dy8AO3bs4MYbb6Rly5ZcccUVbNu2DYD+/fszatQo2rdvz5gxYwq14vr378+QIUNo3749Z599Np999hkDBgygUaNG9O/fPz+ejz/+mLZt29KiRQu6devG4cOHAffuuAkTJtCiRQsaN27Mtm3b2LlzJy+99BJPPfUUzZo144svvuD999+ndevWNG/enOuuu46///77hP482rdvz6BBg5g9ezYAr7zyCpdccglNmzbl9ttvJzMzk3Xr1rF8+XLi4+Np1qwZO3bsKPa6UxWyiXJXetYJva5Ucbw9k7T0LAz/9kxONVlmZWUV6novWrSIqlWr8vzzz9O/f38WLlzIgQMHuPfeewE4cuQILVq04Pvvv+eqq65i0qRJAAwaNIjnnnuODRs2MGPGDIYOHZp/j59//plPP/2UmTNnHnP/AwcOsGrVKp566ik6derEgw8+yI8//siWLVtISUnhn3/+YerUqXz66ad8//33tGrVilmzZuV/vkaNGnz//fcMGTKEGTNmUL9+fQYPHsyDDz5ISkoKV1xxBZdffjlff/01GzdupGfPnjz55JMn/OfUokWL/OTftWtXvvvuOzZt2kSjRo147bXXuOyyy+jcuTNJSUmkpKRwzjnnFHvdqQrZrnftWAdpxSTF2rEOC6JR5dXxeianMoRTUtf7+uuvZ8mSJdx///1s2rQp//WIiAh69OgBQO/evenatSuHDx9m3bp1dOvWLf+67Ozs/J+7detGZGRksffv1KkTIkLjxo2pWbMmjRs3BuCiiy5i586dpKam8tNPP9GuXTsAnE4nbdu2zf98165dAWjZsiXLli0r9h6pqan06NGD3bt343Q6T2g5jlfBM71++OEHxo0bR3p6OocPH6ZDhw7Ffqa0152IkE2U8R0aFhqjBHDYIonv0NDCqFR5E+ieSV5eHlu3bsXhcLB//37q1q1b7HUiQl5eHrGxsSWOdVasWLHE+0RHRwPuBOz92fs8JyeHyMhIrr/+ehYsWHDcz0dGRpKTk1PsNcOHD2fUqFF07tyZNWvWMHHixBLjKcnGjRvzxyf79+9PcnIyTZs2Ze7cuaxZs6bYz5T2uhMRsl3vLs3rML1rY+rEOhCgTqyD6bddrBM56oSU1APxV8/kqaeeolGjRixYsIABAwbgcrkAdwJdunQpAPPnz+fyyy+nSpUqNGjQgCVLlgDu1lfBVuipaNOmDWvXruXXX38FIDMzk59//vm4n6lcuTKHDh3Kf56RkUGdOu7/3t54440TjuGzzz5j9uzZ+cMPhw4dolatWrhcLubNm1fifUu67lSEbKIEd7Jcm3ANvyfezNp+NeiyoR8c3GV1WKocie/QEIetcPe1LHomRccoExIS+Pnnn3n11VeZOXMmV1xxBVdeeSVTp04F3K3DH3/8kZYtW7Jq1SrGjx8PwLx583jttddo2rQpF110Ee+9994pxeUVFxfH3Llz6dWrF02aNKFNmzb5Y4Ul6dSpE++++27+ZM7EiRPp1q0bV1xxBTVq1CjVfRctWkSzZs04//zzmTZtGu+8805+i3LKlCm0bt2a66+/ngsuuCD/Mz179iQpKYnmzZuzY8eOEq87FVLezvVu1aqVOal6lH//BK9eB6dfAP3/C7aYsg9OlQtbt2497nKTopI3ppG0cju70rOoHesgvkPDgPdMKlWqlD/rrE5dcb8DIrLBGNOquOtDdozyGDUvhK4vw6LesGIU3PoC6O4MVQpdmtfRIZswF9Jd72M06sS2hkMgZR6Txg3X3Tqq3NDWpLXCp0WJuwv1yE9X8jTf8ljU2/x9sBpjlzkB3a2jlCpZWLUok1ZuJ9NleMA1jA3mPJ62vUDL3BTdraOUOq6wSpTetW9HieYeZzw7TB1ets3i9IzNFkemlApmYZUoC659O0hF+jrHsNfE8lb0E5C2wcLIlFLBLKwSZdE1cXupxt1mPFLhNHjrNti10cLoVDgREfr06ZP/PCcnh7i4OG655ZYT+p769evzzz//nPA1/fv35+WXXy70WnJyMh07dizxe/r375+/6H3gwIHFljYrqYRaQUXLsb300ku8+eabx/2M1cIqURa3W2dE1/Z82W4Ou4/aSX/5ZgZMe0VnwpXfVaxYkR9++IGsLPdw0CeffJK/iyUQevXqxcKFCwu9tnDhQnr16lWqz7/66qsnXdasaKIcPHgwffv2PanvCpSwSpRQZLdOwjUAjPxoP92yx3HIVODp7PHMX/auJkvldzfddBMrVqwAYMGCBYWS1P79++nSpUv+rpjNm93j6Pv27eOGG26gefPm3HfffYWKRrz99ttceumlNGvWjPvuu4/c3MLFPAq67rrr2LZtG7t37wbcWxQ//fRTunTpwuTJk7nkkku4+OKLGTRoEMVtSrn66qvxbvyYM2cO559/PldddRVr167Nv6a4MmvFlWObOHEiM2bMACAlJYU2bdrQpEkTbrvtNg4cOJB/vzFjxnDppZdy/vnn88UXX5zUn/nJCrtEWZS3OkyqiaOncxwHTCVejZjKsveWaXV05Vc9e/Zk4cKFHD16lM2bN9O6dev89yZMmEDz5s3ZvHkz06ZNy29xTZo0icsvv5yNGzfSuXNn/vzzT8C902TRokWsXbuWlJQUIiMjj7vPOTIykq5du7J48WLAXQG9ffv2VK5cmWHDhvHdd9/lt3g/+OCDEr9n9+7dTJgwgbVr1/LJJ58U6o4XV2atuHJsBfXt25cnnniCzZs307hx4/xycuAenvj22295+umnC70eCGG1jrI4BavApBFHT+djzLdP5UWmMiAjnjQaaXX0UPZhAvxvS9l+5xmN4Sbfxy40adKEnTt3smDBgmPGBr/88kveeecdAK655hr27dtHRkYGn3/+eX5Zs5tvvplq1aoB8H//939s2LCBSy65BHDvJT/99NOPe/9evXoRHx/PiBEjWLhwYX4yXr16NU8++SSZmZns37+fiy66iE6dOhX7Hd988w1XX301cXFxgLtCubd4xomWWcvIyCA9PZ2rrroKgH79+hUqIVewtNvOnTuP+11lLexblEWrwOymOj2c4/mfOY037E/QLsL9H5FWR1f+0LlzZx566KFjxgaL6+56D8Uq7nAsYwz9+vUjJSWFlJQUtm/f7rOsWbt27di9ezebNm1i3bp1dOzYkaNHjzJ06FCWLl3Kli1buPfee32eWFjSYV3Dhw9n2LBhbNmyhZdffvmUz1EvTWk3fwn7FmVxdSv3UI0ezsd42z6N120zuM81kjV5zbU6eigqRcvPnwYMGEDVqlVp3LhxobqJV155JfPmzeOxxx5jzZo11KhRgypVquS/Pm7cOD788MP8Mbxrr72WW2+9lQcffJDTTz+d/fv3c+jQIc4666wS7y0idO/enX79+tGxY0diYmJIT08H3BXMDx8+zNKlS7njjjtK/I7WrVszYsQI9u3bR5UqVViyZAlNmzYFSi6zVrlyZQ4ePHjMd1WtWpVq1arxxRdfcMUVV/DWW2/lty6tFvYtyuJmwqtVsLGPqvRyjmO7qcts2yw6RHyr1dFVmatbty4jRow45vWJEyeyfv16mjRpQkJCQn6imTBhAp9//jktWrTg448/pl69egBceOGFTJ06lRtuuIEmTZpw/fXX50/UHE+vXr3YtGkTPXv2BCA2NpZ7772Xxo0b06VLl/yufElq1arFxIkTadu2Lddddx0tWrQo9O9QXJm1ouXYCnrjjTeIj4+nSZMmpKSk5JeTs1r4lFk7AQVPcKzCEebYn6Sp7CCl5XRadR7s13sr/zvRMmsq9JxombWwb1EWp2Ar8xAVGeOYxIG4S2j1fQJsOPFKzUqp8i3sxyhLckwNQtc17lqW7z/AUx9u4tnD11pWxFUpFVjaoiwtm4PlF8zgE3MJD+a8xuDI5WV2dKlSKrj5LVGKyJkislpEtorIjyJyzIi1iNwlIps9j3Ui0tRf8ZSFJz79nSHZw3kv9zLG2BYyKmoxWa4cXTZUDpW3sXlVdk7m796fXe8cYLQx5nsRqQxsEJFPjDEFd9L/DlxljDkgIjcBs4HWxX1ZMNiVnoUhigddQ8kydh6ISqYC2Tye3tvq0NQJiImJYd++fVSvXr3ENYAqNBlj2LdvHzExJ3Zmlt8SpTFmN7Db8/MhEdkK1AF+KnDNugIf+Roo/hDjIFE71kFaehZ5RDA2ZyBZRDMw6kNio1xcPj2GtIxsHbcsB+rWrUtqaip79+61OhRlgZiYmBLPSy9JQCZzRKQ+0Bz45jiX3QN8GIh4TlbBxemGCCbl9OWoxDAk8j0iMo8Sz3263bEcsNlsPrfTKVWQ3xOliFQC3gFGGmOOXY7vvqY97kR5eQnvDwIGAfkLbK3gTXz/Hl1agdnO3hzKjuZh22JicDLCNYwsl/saTZRKhQa/LjgXERvwAbDSGDOrhGuaAO8CNxljfvb1nYFYcH4iGiSswAB3R37IBNtbrMptxhDXSJzY+T3xZqvDU0qVkiULzsU9Sv4asPU4SbIesAzoU5okGYy82xrn5N5EgmsgV0dsYo7tSc6panFgSqky4891lO2APsA1IpLieXQUkcEi4t0HOB6oDrzoeT94moqlVPB4iYW51/CgawiXRmxjUYUnISvd2uCUUmVC93qXgeSNaQXGLR3MavwXl6wfxS/Uo1fWwzhia+pMuFJB7nhdb93CWAaKbndM3pjGfTnxPBcxg4X2KdyV/ghjlznzr1VKlS+6hdEPklZu5xNXY/o5E6gt+1hin8xprv/pDh6lyilNlH7gLfD7jWlEH+dYqskhFkdPxp7xm8WRKaVOhiZKPyhY4HejOY87neOIwcni6Cn0njZXDyxTqpzRROkHBWfCAX409entGk+egWezH+VC+V0rDylVjmii9IPijpf4X0x9ujknkEkMC+xTaSE/64FlSpUTujwoQLw7eGqxj3n2x6kpBxjoeoiv8y7SHTxKBQE9CiIIeMctvcfh/mVOZ67tSW6J2Uy7xFU6bqlUENNEGSAFxy33EktP5zh+NnWZaZJofPAzDOi4pVJBShNlgBQdt6wYezpDIyewyZzDC7ZnuC3CfWynjlsqFXx0Z04AFd3B0yBhBX1J4BXbTGbaXsKR42R+7rX56zCVUsFBE6WF3BXT4R5XPC/YnmGa7TUcZPOO/VbaJa7K3zuu+8SVspZ2vS3kHbfMxs4Q14N8kNuax2xv0z9nMWnpmTpuqVSQ0ERpoYLjljlE8WSFeJZzFSMjlzAmaiHgXrql45ZKWUu73hYrOm55dsJRMqJsDIl6HwfZTMrpiyFCxy2VspAmyiBTK7Yij6XfTSbR3Be1ggpkk5BzL7ViK1odmlJhSxNlkPGe9DjddSeZJoYHbe9QKcKJ6/qXrA5NqbCliTLIFDzp8dn027E5KjEs5w3YnsByM40nPt2ps+FKBZgmyiBUeNzyZviuCawYTfXtaezPHokhRs8PVyqAdNa7PLhkIFNtw2nDFt6wP0FlMgGdDVcqUDRRlhOvHWrLcNdwmsuvvG2fRlUOA+hsuFIBoImynKgd6+C/eW24z/UgF8hfLLRPoQYZhaqpK6X8QxNlOeHdxbMqrwV3u+I5S/awJHoyt55ttEybUn6mibKcKLiL56u8ixkVPYFakQe588f7iMjYqdsdlfIjrXBejg2Y9gozsyeSjY27nI+ww7hnv+vEOlibcI3F0SlVvmiF8xC1+mBtejgfI5I8Ftmn0Ej+AHSCR6mypomyHKsd6+BncybdneNxEsVC+xSaya86waNUGdNEWY55J3h+N7Xo7pxAuqnEPPs07q6bphM8SpUhTZTlWMEJnjQTxwOO6RyJOYO7fhnFOQe/1gkepcqITuaEmJunv8uTWRM4V9IY7hrOx3mXADrBo5QvOpkTRn7KsNPL+Sg/mfq8aHuGzhHrAJ3gUepUaKIMMbVjHRykEr2dY1lvGvK07QW6Ra7RCR6lToEmyhDjneA5goP+zof5Iq8xSbbZjK3+mU7wKHWSNFGGmIITPNlEM6HCOH6qegW3pD1N50OLdIJHqZOg9ShDUNFzeK6cLozKzWWMbSEV5Cgzc7rll2jTWpZK+aaJMgz8leFiFEPJMtEMj0qmIkeZnNNHJ3iUKiVNlGGgdqzD3d3OGUgW0QyI+ggH2SRF3Ue7xFV6tIRSPugYZRjwTvCAMDmnD8/n3EqvqNVMzHuO/6Uf1nFLpXzQRBkGCk7wCMKCSv15Tu6kc8RaXrA9ix0XoEdLKFUSv3W9ReRM4E3gDCAPmG2MeabINQI8A3QEMoH+xpjv/RVTOCs6wdMgIYv9kZFMsL3FbGYx2DWSo0TruKVSxfBnizIHGG2MaQS0Ae4XkQuLXHMTcJ7nMQj4jx/jUQXUjnUwJ/cmxrju5cqIzcyxJVGRLF2YrlQx/JYojTG7va1DY8whYCtQdKbgVuBN4/Y1ECsitfwVk/qXd9xyUW57RrqGcknENuZFJ/JI+zOsDk2poBOQWW8RqQ80B74p8lYd4K8Cz1M9r+0ORFzhzNsNT1q5nffT21HRXonHc2bSbOMg/mv+w+Nr9upsuFIefk+UIlIJeAcYaYw5WPTtYj5yTDkjERmEu2tOvXr1yjzGcFV43PJm+PVScuffyfkf9sCV/QiGavmz4d7rlQpHfp31FhEb7iQ5zxizrJhLUoEzCzyvC+wqepExZrYxppUxplVcXJx/glVw7nU8EDWOM9jHEvsk6speQGfDlfJbovTMaL8GbDXGzCrhsuVAX3FrA2QYY7TbbaH/HjyH3s5HiJXDLLJPpoG4/zp0NlyFM3+2KNsBfYBrRCTF8+goIoNFZLDnmv8CvwG/Aq8AQ/0YjyqF2rEOUsy59HKOIwYni+2TOV/+0tlwFda0wrkqJHljGmOXbSHLlcs5ksZ8++NEk8PGq+fQvv0NVoenlN9ohXNVagV38fxm6jA8ehr2ClW4fO3dDJ72gtazVGFJW5TKp5Xr1tNwZW9OZz8DXaNZl3cxDlsk07s21plwFTK0RalOyeTPD9It+zH+NKczx5ZE+4iNOhOuwoomSuXTrvQs9hJLT+c4tpkzmW2bRceIr3UmXIUNTZTKJ++MdzqV6e18hI3mXJ6zPceAyl9bHJlSgaGJUvn0bz1LOEQF+jnH8A0X85jrWfjuNYujU8r/NFEqnwrXs4TTYqvxT6c3+V/Nq2HFKB5/dKjOhKuQpkdBqFIpWs8yeWMa1++6l+kc5VHbPByHshm77I78a5UKJZoo1UlJWrmdQy7hAYaRZeyMsi2lQs5Rkj4aoIlShRxNlOqkeGe884jg4ZxBZBHN4KgPqHgkG/LaQ4SO6qjQob/N6qQU3PttiGB8Tn9eyulEn6hP4L37ITfHwuiUKluaKNVJKTgT7iY8I3fxf2fcA5vms2JiR66avlIneFRI0K63OikFK6R7K6G3vyCOYRtu4K68o4yzzSM6czqjl41i/R/7Wb1NK6ar8kv3eqsy0y5xFWmescu7Ij/lcdvrfJl7EYNco8kkJv863SeugtEp7fUWkRgRuUNEnhGRJSLypog8LCIXlX2oqjwruKVxXu51jHIOpm3ET7xhT6Qymfnv6T5xVd4cN1GKyERgHdAW98FgLwOLcR9Fmygin4hIE38HqcqHosV9l+VdyTDXAzSTHbxtn0Ysh/Lf033iqjzx1aL8zhjTwhgz2hgz3xjzqTHmA2PMLGNMJ+AuwB6AOFU5cOwED3yU15pBrlFcIH+x0D6VONKBY5OqUsHMV6J0iEh0SW8aY/YYY3TAUAHHbnWsE+vgrjb1+DqyFXe74qkne1hon0ID2wHiOzS0OlylSu24kzki8i7us28+AhYAHxtjcgMUW7F0Mqf8Sd6YRtLK7dTKSGFudBIRjli+bPc6k77M0plwFTSON5njc9ZbRKoAtwE9gabAe8ACY8znZR1oaWiiLOd2bcQ551YOOCO40/kIO4w7OepMuLLaKc16G2MOGmPeMMbcBDQGUoDnROSvsg1ThYXazRkgk4ggj0X2KVwoOwGdCVfBrdQ7c0SkGtAV6AGcBrzjr6BUaFt78HS6OceTjY0F9qk0k18BnQlXwcvX8qDKItJHRP4LbAUuAaYC9YwxIwMQnwpBtWMd7DS16J49ngOmMm/bp9FatupMuApavlqUvwM3Av8BzjTGDDLGrDLlbTuPCireZURpxNHdOZ7dpjpv2BNJar7X6tCUKpavvd71jDGZPq5R6oQU3icOIx2PMy8mkcu+HQb1KkKjW6wNUKkifC0PWg7MBlYaY1xF3jsb6A/sNMa87s8gC9JZ7xCVdQDevgN2bYTbXoYm3ayOSIWZU5n1HgRcCWwTke9E5L8iskpEfsO9nXFDIJOkCmGOatA3mb3VW5K37F4SHn1Iz+FRQeO4XW9jzP+Ah4GHRaQ+UAvIAn7WLrkqa8k/HWTi38N5mhkk2l5h4qFsxi67GdBzeJS1Sl2P0hizE9jpt0hU2EtauZ10VySDGMWzPM9E25tUcGWTtNKuiVJZSgv3qqDhXUfpxMb9rgeYwUs8bFtExcNHwbQHEYsjVOFKE6UKGrVjHfmFf3OJZLRrCFnGzv1R77F4am/GHOlF7dgKui9cBZyemaOCRtEybXlEMCHvXubm3kj33A+YFvUqu9OPMHbZFp3kUQF13BaliGwBils/JIAxxmjRXlVmijuHJ9OZw8TMPhwyMQyPSsYh2TzkGkzSyu3aqlQB46vrrSt/VUB1aV6nUAJskLACEGbmdCfLRPOwbREOnDyQPty6IFXY8bU86I9ABaJUcQqOW76YeyuZRDPR9iavR87imunR/J6Rp/Usld+VaoxSRNp4FpwfFhGniOSKyEF/B6dU0XHLubk3MjZnEG3NZqZlTaICWaSlZ+m4pfKr0k7mPA/0An4BHMBA4Dl/BaWUV3HHS3xkv56RrvtpJdt52z6dKhzWepbKr05kwfmvIhLpOQpijois82NcSuUrbtxyOZeR5bLzvO1ZFtofp48zgV3p1sWoQltpW5SZImIHUkTkSRF5EKh4vA+IyOsiskdEfijh/aoi8r6IbBKRH0Xk7hOMXYUpb93KT/Jaca9rNGfLLhbZp3BezEHaJa6iQcIK3SeuylRpE2UfIBIYBhwBzgRu9/GZubhrWZbkfuAnY0xT4GpgpicZK3VcBcctP89rSj9nAmfIfl4145GMPzGg45aqTJUqURpj/jDGZHnOz5lkjBlljPnVx2c+B/Yf7xKgsogIUMlzbU5pA1fhq+i4ZWrVFgyJGE9VDrPYPokGshvQc3hU2fFVj3KxMaZ7SQvPfS0491Qc+sAYc3Ex71UGlgMXAJWBHsaYFb4C1nqUqjgNElZwgfzBW/bpGITezrFsN/UA9wSQHourfDmVepQjPP+8BehUzONUdMB9omNtoBnwvOdo3GOIyCARWS8i6/fu1eMC1LFqxzrYas6ih/MxcolgoX0qjeU3BHc3XLvj6lQcN1EaY3Z7fhzq6X7nP4Chp3jvu4Flxu1X3OfzXFBCHLONMa2MMa3i4uJO8bYqFHnHLXeYOnRzjuewcTDf/jgtZVuh67Q7rk5GaSdzri/mtZtO8d5/AtcCiEhNoCHw2yl+pwpTBcctU01NhjumscfE8qb9CdpFbCl0rR6Lq06Ur6IYQ3C3HM8Wkc0F3qoMrPXx2QW4Z7NriEgqMAGwARhjXgKmAHM9458CjDHG/HOS/x5KHbPe8pbpUSRljed12wyGuEawKq8FgB6Lq06Yr8mcqkA1YDqQUOCtQ8aY481o+41O5qjSSt6YRuKyr3hZHudC+YORrvtZFdmO21vWYfW2vTrBowo53mSOr6IYGUAG0EtEIoGans9UEpFKxpg/yzxapcqIO/m15aGPpjAtazLP2p8nuV5lxm1wj1XCvxM8/16v1LGO26LMv0hkGDAR+BvI87xsST1KbVGqk+I8Agvvgt9WM851N2/nFh52rxPrYG3CNRYFp4LBSbcoCxgJNDTG7CuzqJQKJHtF6LWQTyZ3ZKptDg6yeSX333KrOsGjjqe0s95/4e6CK1V+2WKYUiGB93Pb8KhtPiOjluLdR6ETPOp4Stui/A1YIyIrgGzvi8aYWX6JSik/GXXjRSQsG8HRHDsjo5bhIJunpQ/xHRpaHZoKYqVNlH96HnbPQ6lyyTthM+OjEWQdiea+qBVcX78yZzftQPLGtELn9ehsuPIq1WRO/sUiFY0xR/wYj086maPKjDHw6QRY+wx/ntmFjju7c9j179sOWyTTuzbWZBkmTmWvt/cL2orIT8BWz/OmIvJiGcaoVOCJwHWToP2j1PsrmUSewVaggJVud1RepZ3MeRp3EYt9AMaYTcCVfopJqcARgaseZqrrLm6J/Ib/2J4iGmf+2zobrqD0iRJjzF9FXsot41iUssyHle/gEdc9XBORwmu2JBwcBXQ2XLmVenmQiFwGGBGxi8hDeLrhSoWC+A4NeTfiBka7BtM24ifetCcSZ8um/QVxeryEKvWs92DgGaAOkAp8zKmXWVMqaHgnbJJW2hl2MNq93dGRyB0b4tntqgDodsdwVtotjO2MMWt9vRYIOuutAuLnj8mefxe/59Wkj3Mse4nNf0u3O4amU571pvgzvPVcbxW6zr+Bu53xnCl7WGSfTC3+3b2rEzzh57iJ0rMsaDQQJyKjCjwm4j6VUamQ9UeVVvRxjqWGHGRJ9CTqyd+ATvCEI18tSjvuExKjcBfr9T4OAnf4NzSlrBXfoSFboxrRy/koFTjKEvskLrbt1u2OYai0Y5RnGWP+8JycaIwxh/0fWvF0jFIFkndbY8WMX5gfPZ3Kdlh32auM+zpCtzqGmOONUZY2UV4MvAWc5nnpH6CfMeaHMouylDRRKsvs20HmqzfjyjxIf+fDbDTnAbrVMVSUxWTObGCUMeYsY8xZwGjPa0qFj+rn0DtvIgdMJd6yT6dNxE+AbnUMB6VNlBWNMau9T4wxa4CKfolIqSC2MaMy3Zzj2WWqM9f2BFdFbAJ0JjzUlTZR/iYij4lIfc9jHO5zuJUKK7VjHeylGj2dj/GrqcMrthl0iPhOZ8JDXGkT5QAgDljmedQA7vZXUEoFq/gODXHYItlPFe50PsoPpgEv2J5hRM2NutUxhPlaRxkjIiNxn8H9I9DaGNPCGDPSGHMgEAEqFUy6NK/D9K6NqRPr4BAVGeOYzJ+Vm3LHzilcfui/GP7d6qjJMnT42uv9BuACvgBuAhrhPmhMqbDVpXmdQjPc7adHMD5vOk/YXqEiR3k996b8CR6dCQ8NvhLlhcaYxgAi8hrwrf9DUqp82ZmRx32M4hmeZ7ztLRxk80JuF53gCSG+EmV+YXxjTI6I+Dkcpcqf2rEO0tJhmOsBkniZeNtiKshRXonqTbvEVbowPQT4msxpKiIHPY9DQBPvzyJyMBABKhXsvBM8uUQy2jWY+TnXcH/Uckblvsau9CM6bhkCjtuiNMZo4QulfPi3lqX7BMcXKg4jL9tB38gV2E02j+QMJI8IHbcsx0pbuFcpdRxFJ3gaJGSxP8rGA1HJOMTJaNdgcojScctyShOlUn5QO7YCs9K7k2ViGGNbSAxOhruGExdbxerQ1Eko9eFiSqnS845b/ie3MxNc/egQuZ7Xo2eScG09q0NTJ0FblEr5QcFxyzfTOxDtqMTYnBeRHx7gA/M001el6mx4OaKJUik/KTxueTNsaUHeskGcubMXh7IfxlBJDywrJ7TrrVSgNL6DsVHxXMBOFtofpzoZgJZpKw80USoVQIsPNWGg6yEayG4W2adQk/2AlmkLdpoolQqg2rEOvshrQl9nAjXlAEvsk6gre6jqsGn1oSCmiVKpAPLOhn9nLqC3cyxVJJOl9smc7vyTtPQs3cUTpDRRKhVABcu0bTbn8kD0VGySy7yoSVwgf+Zfp+OWwcVvs94i8jpwC7DHGHNxCddcDTwN2IB/jDFX+SsepYJF0V08147N4W37dBbap9DXmcBmcw6g45bBxJ8tyrnAjSW9KSKxwItAZ2PMRUA3P8aiVNA6WvVcujnHc9BUYJ59Gq1kG4AeLxFE/JYojTGfg2dKr3h3AsuMMX96rt/jr1iUCmbxHRqyL6oW3Z3j2WNiecueyLW2H4nv0NDq0JSHlWOU5wPVRGSNiGwQkb4WxqKUZbzjlpGxdenpHE9qRG1m256ki2OT1aEpDyt35kQBLYFrAQfwlYh8bYz5ueiFIjIIGARQr57ulVWhp9C4ZeaN8Pbt5C3qw8SoEbx1qKVudbSYlS3KVOAjY8wRY8w/wOdA0+IuNMbMNsa0Msa0iouLC2iQSgVchdP4oNlLfJ93LhOcs7g98jNdMmQxKxPle8AVIhIlIhWA1sBWC+NRKmhMX5VGn+yHWZt3MTNsL9M78hNdMmQhfy4PWgBcDdQQkVRgAu5lQBhjXjLGbBWRj4DNQB7wqjHmB3/Fo1R5sis9C0MM97pG8zzPMdU2h4ocZXZ6J6tDC0t+S5TGmF6luCYJSPJXDEqVV+4Dy7LIxs4Q1wie4kXG2hZwms1Fu+kx7Mo4quOWAaQ7c5QKQt6tjgA5RDHCNYwluVdzH0vpd+Q1DEbHLQNIE6VSQajgVkcBasVWZHrUEObm3MCgqBVMjXodIU/HLQNEC/cqFaSOPbBsBRPpRxbRDIl6H4c4edg1SLc6BoAmSqXKCe+45RM5PTliYnjItoQYsnmywkNWhxbytOutVDnx77il8HzubUxx9ebmyG9ZFPsCuI5aHV5I00SpVDlRdNzyo8q3s7HpRGr+/Tnrp1/HRQnvaNFfP9Gut1LlSNFxy+SNDVmw8R+mmxd4w57I3ekP62FlfqAtSqXKsaSV21nsvIz7XSNoIjuYb59KjOuAzoSXMU2USpVj3hnvj/IuZZBrFOdJGgvtU3Gl77I4stCiiVKpcqxgcd81ec252/UwdWUv7zimQvpfFkYWWjRRKlWOFdzBA/BV3kUMzBvHGZGHYc5NsG+HhdGFDk2USpVjRWfC68Q66NH1dr5sN4cDGRn8/ew13DXtDZ0JP0U6661UOXfsTHgaQ1flUidnHPPs03g2+1HuXZYJdNaZ8JOkiVKpEJO0cjtZrlx+pS7dneOZZ5/G3IjJPPQhQGeSVm5nV3qWVh86Adr1VirEFNz7/Yc5g+7Z49lnKvNU9gSWLVtIWnoWBrT60AnQRKlUiCl6zO0uatDdOZ40E8fsiOlcHZGS/55WHyodTZRKhZiiM+EAh2016OEcxy+mDrNtM+kQ8W3+e1p9yDdNlEqFmOJmwqd3bUyF2Jrc6RzHFnM2L9iepUvEl8CxLVB1LJ3MUSoEFZ0J9xq7bAt9nGN51TaDWbb/UMW4aNHhQQsiLF+0RalUmPC2NKvFVmOA62G+jmzB5IhX6HI02erQgp4YY6yO4YS0atXKrF+/3uowlCr/cpzwzj2wdTkvR91J4uFbwnrJkIhsMMa0Ku497XorFa6i7Lx37lTkxwPclzOf3KjDPJneQ8u0FUMTpVJh7MlPdrDLeR+HouwMjVqOg2wmu/qQtHK7JsoCNFEqFcZ2pWdhiODRnAFkYWdg1IdUIJtH0gdaHVpQ0USpVBjzHlgGwtSc3hwhhhFR71IlysWV06P5K8MV1uOWXjrrrVQYK7w4XXgqpxtJOb24ibWMy0zEhku3OqKJUqmwVtzi9Pn2rox39eOGyA28aptBDNlhv9VRlwcppQppkLACA3SLXENi1CusNw25x/kQh6lAnVhHyFYeOt7yIG1RKqUK8W5pXJJ7NSNcw2ghv/C2fRqxHA7bykOaKJVShRQct/wgry1DXCNpJH+ywD6V6mTkXxdO3XFNlEqpQoqOW26tcjn3uOKpL/9jsX0yNdmff224VB7S5UFKqWMULarRLhH6ZiTwuj2JJfZJ3Ol6lFRzethUHtIWpVLKp/gODfkh6iLucj5CFclkiX0yjWx/0/6CONolrqJBwgraJa4K2TFLTZRKKZ+83fF9VS+ml/MxoiWXJdGT2bJhXVhM8OjyIKXUifvnF/Y83wG7OUpfZwKbzTn5b9WJdbA24RoLgzs5ujxIKVW2apzH7dmPcdBUYJ59GpfItvy3QnGCRxOlUuqk5FU9i27OCewxsbxpT+TyCHd5tlCc4NFEqZQ6KfEdGnLQFkd353h2mjN4zZZER9tG4js0tDq0MqeJUil1UrwTPDGxZ9DLOY5fIxrwfNQsukR9ZXVoZc5v6yhF5HXgFmCPMebi41x3CfA10MMYs9Rf8Silyl6h9ZZHb4IFPeGdgeDKghZ9rA2uDPmzRTkXuPF4F4hIJPAEsNKPcSilAiGmCty1FM5pD8uHwTezrY6ozPgtURpjPocCe52KNxx4B9jjrziUUgFkr8DyRrP4POJS+DCe/0y9PyTWVVo2RikidYDbgJesikEpVbaSN6Yx5r3tDMgcxvLctgzJeZtd7z5G8vepVod2Sqzc6/00MMYYkysix71QRAYBgwDq1avn/8iUUiclaeV2sly5QBQjXfeTZaIZGvUOC1a4SGY6SR//XC7rWVqZKFsBCz1JsgbQUURyjDHJRS80xswGZoN7Z04gg1RKlV7BxeZ5RJCQM5BMormb5SxMzmSXsz+GiPztjlA+jsW1rOttjGlgjKlvjKkPLAWGFpcklVLlR9HF5oYIJuX05T85nekZ8SkzbS8RSS5QvupZ+i1RisgC4CugoYikisg9IjJYRAb7655KKWsVPqzMzWGL4omcniS5utM18kuesz2HjRyg/Gx39FvX2xjT6wSu7e+vOJRSgePtRiet3F5oLDJp5XZeSO9CFtGMt71FDLMY4hqJw1GRdomrgn7cUgv3KqXKVNGiv15jl23hdddNZGHn8ajXmWt/ksHOeNKyogGCetxStzAqpfyu4PESC3OvZYp9BJdEbOP1yGlU4Uj+dcE6bqktSqVUQBRuad7M4EfgWdtzzLc/Th9nAgeoAgTnuKW2KJVSlthS5SoGuUZzrqSx0D6VOA4AwVmmTROlUsoS8R0a8k1kS/q7xlBX9rLYPpmzbQeCskybJkqllCW845Z/VWlJb+cj1Ig4xAeVptKl3lGrQzuGnpmjlAoOuzfBW7dxNFcYaB5j7cG4gC4Z0jNzlFLBr1ZT/q/1HA4ezeXZ7Ee5UH4PmpMdNVEqpYLG+K9yucM5nkxiWGB/nBbyc1AsGdJEqZQKGrvSs/jT1KR79nj+MVV4yz6dthE/kpaeRbvEVTRIWEG7xFUBb2FqolRKBQ3v0qBd1KCHczypJo45tie5OiKFtPQsDFjSHddEqZQKGgWLauwllp7Ocfxi6jDbNpMbI77Nvy7Q3XFNlEqpoFFwq6MAFWJrcqdzHJvNObxge4bbIr7IvzaQO3h0C6NSKqgULarRLnEVfdITeMU2k5m2l3DkOJmfe21Ad/Boi1IpFdTiOzQEW0XuccWzOq8Z02yvMdj+UUB38GiiVEoFNW93vEZsVYa4HmRVRFsSIt6ky8F5EKANM9r1VkoFvULd8dyb3eeGr34cnEfguong44DCU6WJUilVvkRGwa0vgs0Ba58GVybc+ARE+K+DrF1vpVT5ExFBcu3RzI+8Fb6dzfuP30Hy93/673Z++2allPKT5I1pjH33Bx450p2nc7rSKff/sCcP4r0NO/1yP02USqlyJ2nldrJcuYDwdM4dTHf1omPEV5y2YiC4yr5Mm45RKqXKnaKLzV/O7UQm0UyxzWXP7NvoefABfs/IK7MybdqiVEqVO8UtNn8r9wbGM5Tqe75ietZEKpJZZvvCNVEqpcqdgnvCvRy2SN6PaM8DruE0kx20jPgFKJt94ZoolVLlTtE94XViHUzv2pj0TBcr8tpwRfbTfJbXNP/6U90XrmOUSqlyqeiecHBP8qSlZ7GHaoVeP9V94dqiVEqFjJK65Ke6L1xblEqpkOFtYSat3M6u9Kwym/XWRKmUCinFdclPlXa9lVLKB02USinlgyZKpZTyQROlUkr5oIlSKaV80ESplFI+aKJUSikfNFEqpZQPYgJ0illZEZG9wB+luLQG8I+fwyktjaV4GsuxgiUOCL9YzjLGxBX3RrlLlKUlIuuNMa2sjgM0lpJoLMEbB2gsBWnXWymlfNBEqZRSPoRyopxtdQAFaCzF01iOFSxxgMaSL2THKJVSqqyEcotSKaXKREgmShG5UUS2i8ivIpJgYRxnishqEdkqIj+KyAirYvHEEykiG0XkA4vjiBWRpSKyzfNn09bCWB70/N38ICILRCQmgPd+XUT2iMgPBV47TUQ+EZFfPP+sdrzv8HMsSZ6/o80i8q6IxFoVS4H3HhIRIyI1AhGLV8glShGJBF4AbgIuBHqJyIUWhZMDjDbGNALaAPdbGAvACGCrhff3egb4yBhzAdAUi2ISkTrAA0ArY8zFQCTQM4AhzAVuLPJaAvB/xpjzgP/zPLcqlk+Ai40xTYCfgbEWxoKInAlcD/wZoDjyhVyiBC4FfjXG/GaMcQILgVutCMQYs9sY873n50O4E0LZll4uJRGpC9wMvGrF/QvEUQW4EngNwBjjNMakWxhSFOAQkSigArArUDc2xnwO7C/y8q3AG56f3wC6WBWLMeZjY0yO5+nXQF2rYvF4CngYCPjESigmyjrAXwWep2JRcipIROoDzYFvLArhady/ZHkW3d/rbGAvMMczDPCqiFS0IhBjTBowA3cLZTeQYYz52IpYCqhpjNkN7v/RAqdbHI/XAOBDq24uIp2BNGPMJivuH4qJUop5zdKpfRGpBLwDjDTGHLTg/rcAe4wxGwJ972JEAS2A/xhjmgNHCFz3shDP+N+tQAOgNlBRRHpbEUswE5FHcQ8jzbPo/hWAR4HxVtwfQjNRpgJnFnhelwB2p4oSERvuJDnPGLPMojDaAZ1FZCfuoYhrRORti2JJBVKNMd6W9VLcidMK1wG/G2P2GmNcwDLgMoti8fpbRGoBeP65x8pgRKQfcAtwl7FuLeE5uP9ntsnzO1wX+F5EzghUAKGYKL8DzhORBiJixz04v9yKQEREcI/FbTXGzLIiBgBjzFhjTF1jTH3cfx6rjDGWtJyMMf8D/hIR70HL1wI/WREL7i53GxGp4Pm7uhbrJ7uWA/08P/cD3rMqEBG5ERgDdDbGZFoVhzFmizHmdGNMfc/vcCrQwvO7FBAhlyg9g8/DgJW4f+kXG2N+tCicdkAf3C24FM+jo0WxBJPhwDwR2Qw0A6ZZEYSnVbsU+B7Ygvu/h4DtABGRBcBXQEMRSRWRe4BE4HoR+QX3DG+ihbE8D1QGPvH87r5kYSyW0p05SinlQ8i1KJVSqqxpolRKKR80USqllA+aKJVSygdNlEop5YMmSqWU8kETpQoqnhJsQws8ry0iS/10ry4iUuK2OBFpLCJz/XFvVb7oOkoVVDzFQz7wlD3z973W4d51UuIxqCLyKTDAGBPw0l4qeGiLUgWbROAcz06QJBGp7y3gKiL9RSRZRN4Xkd9FZJiIjPJUIfpaRE7zXHeOiHwkIhtE5AsRuaDoTUTkfCDbmyRFpJuneO8mEfm8wKXvE9galSoIaaJUwSYB2GGMaWaMiS/m/YuBO3HXHX0cyPRUIfoK6Ou5ZjYw3BjTEngIeLGY72mHe+ui13iggzGmKdC5wOvrgStO4d9HhYAoqwNQ6gSt9hRBPiQiGbhbfODeq93EU9LuMmCJu84FANHFfE8t3HUxvdYCc0VkMe4qQl57cJdgU2FME6Uqb7IL/JxX4Hke7t/nCCDdGNPMx/dkAVW9T4wxg0WkNe4q8Cki0swYsw+I8Vyrwph2vVWwOYS7Ys1J8RRG/l1EuoG71J2INC3m0q3Aud4nInKOMeYbY8x44B/+rWl6PnDMIVcqvGiiVEHF04pb65lYSTrJr7kLuEdENgE/UvyZSZ8DzeXf/nmSiGzxTBx9DniPHGgPrDjJOFSI0OVBKmyJyDPA+8aYT0t4Pxr4DLi8wCFbKgxpi1KFs2m4T14sST0gQZOk0halUkr5oC1KpZTyQROlUkr5oIlSKaV80ESplFI+aKJUSikf/h8L2DaX3SQwkwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Grabbing the output as well as time used\n",
    "t, modeldata = MyGalvanostatic.get_output()\n",
    "#Saving the output\n",
    "MyGalvanostatic.save_output(\"Discharge_Data_Results.txt\")\n",
    "\n",
    "plt.figure(figsize=(5,5))\n",
    "plt.plot(expdata[::2,0],expdata[::2,1],'o', label=\"Experimental Data\")\n",
    "plt.plot(t,modeldata,label=\"Model Validation\")\n",
    "plt.xlabel(\"time (s)\")\n",
    "plt.ylabel(\"Potential (V)\")\n",
    "plt.legend(loc=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extra Notes\n",
    "\n",
    "The module only consists of one file, which is just an `__init__.py` file. All codes are implemented in there. PyFDAM also do not implement the fitting routines by itself, and delegate it to the LMFIT package. Any issues pertaining to the fitting process can be traced directly to the LMFIT package. I highly recommend looking at LMFIT package for more troubleshooting because non-linear curve fitting, just by the nature of the mathematical problem, can be difficult! "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}