CODE_2_PES2MP_2D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "11183bfd",
   "metadata": {},
   "source": [
    "Contact : dhilip@iitrpr.ac.in <br>\n",
    "\n",
    "    - Apoorv Kushwaha & Dr. T.J. Dhilip Kumar\n",
    "        Quantum Dynamics Lab (410), IIT Ropar."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fc97f14",
   "metadata": {},
   "source": [
    "### 2D multipole expansion code for a rigid rotor-atom collision PES:\n",
    "*Uses Legendre functions from scipy library of python*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "46a1938f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import important libraries\n",
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.special import legendre\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f14eaf97",
   "metadata": {},
   "source": [
    "#### read input file (PES)  give separation (remove header such as r, theta, phi, etc)\n",
    "*The code assumes first column to be R (Radial Coordinate), 2nd to be theta (Angular coordinate) and 3rd column to be E(Potentials)*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d78f36e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_inp = pd.read_csv('2D_PES.dat',header=None,sep='\\s+')  # import file\n",
    "lm = 4  # Difine lambda max "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f8920115",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_inp.sort_values(by = [ 0,1], inplace=True, ascending = True)    # sort by (R, theta) in ascending order\n",
    "E_inf = -78.70465827                                      # define E_infinity (Asymptotic Energy)\n",
    "df_inp[2] = (df_inp[2] - E_inf)*219474.63                     # convert to cm-1\n",
    "df_inp.reset_index(inplace=True, drop = True)                 # Resetting index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9731ae8f",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2645.156525</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.5</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2640.813122</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2.5</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2631.509593</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2.5</td>\n",
       "      <td>3.0</td>\n",
       "      <td>2615.856662</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2.5</td>\n",
       "      <td>4.0</td>\n",
       "      <td>2595.561843</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9095</th>\n",
       "      <td>12.4</td>\n",
       "      <td>86.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9096</th>\n",
       "      <td>12.4</td>\n",
       "      <td>87.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9097</th>\n",
       "      <td>12.4</td>\n",
       "      <td>88.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9098</th>\n",
       "      <td>12.4</td>\n",
       "      <td>89.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9099</th>\n",
       "      <td>12.4</td>\n",
       "      <td>90.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>9100 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         0     1            2\n",
       "0      2.5   0.0  2645.156525\n",
       "1      2.5   1.0  2640.813122\n",
       "2      2.5   2.0  2631.509593\n",
       "3      2.5   3.0  2615.856662\n",
       "4      2.5   4.0  2595.561843\n",
       "...    ...   ...          ...\n",
       "9095  12.4  86.0     0.000000\n",
       "9096  12.4  87.0     0.000000\n",
       "9097  12.4  88.0     0.000000\n",
       "9098  12.4  89.0     0.000000\n",
       "9099  12.4  90.0     0.000000\n",
       "\n",
       "[9100 rows x 3 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_inp # preview"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8388613b",
   "metadata": {},
   "source": [
    "#### Declare Variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ff9d89f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "nc=100                      # number of Radial coordinates (Must be same for all angles)\n",
    "ngm = 91                    # number of angular coordinates (91 for symmetric molecule, 181 for non-symm.)\n",
    "px = np.zeros((ngm,lm))     # Matrix to stores legendre coeffinients\n",
    "f = np.zeros(ngm)           # Array to store part of ab initio energy \n",
    "R = np.zeros(nc)            # Array for distance R\n",
    "E = np.zeros(nc)            # Array to store multipole expanded potentials\n",
    "df_out = pd.DataFrame()     # dataframe that stores final V lambda (Radial coefficients)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d9991b69",
   "metadata": {},
   "outputs": [],
   "source": [
    "V_nf= np.zeros((nc,lm))     # Numpy 2D array to store V_lambdas as they are calculated for each radial term\n",
    "V_n= np.zeros(lm)           # Stores V_lambda for one radial term (Depreciated part of code no longer used!)\n",
    "symmetric = True            # Verify if rigid rotor is symmetric (else put False)\n",
    "if symmetric:\n",
    "    sym = 2\n",
    "else:\n",
    "    sym = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba4955ce",
   "metadata": {},
   "source": [
    "#### Calculate legendre coeff and take pseudo-inverse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "abd2108b",
   "metadata": {},
   "outputs": [],
   "source": [
    "for j2 in range (0,ngm,1):      # loop over anglular terms (considering 0-90 with 1 degree interval)\n",
    "    for j3 in range (0,lm,sym):      # loop over legendre terms\n",
    "        pxc = legendre(j3)         # Uses j3*2 for symmetric molecule (only even V_lambdas); and *1 otherwise\n",
    "        ang = math.radians(j2)         # convert angles to radians\n",
    "        px[j2,j3]= pxc(math.cos(ang))  # store legendre coefficient for corrosponding angle and lambda (2D)\n",
    "\n",
    "np.save(\"2D_L_coeff.npy\", px)    # save Legendre coefficients to numpy readable file for future use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a5bb24a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "px = np.load(\"2D_L_coeff.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "eaf2d75b",
   "metadata": {},
   "outputs": [],
   "source": [
    "A_inv = np.linalg.pinv(px)             # take pseudo-inverse of px matrix (synonymus to least squares fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9b9648a",
   "metadata": {},
   "source": [
    "#### Calculate V_lambda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d4376f74",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>2</th>\n",
       "      <th>4</th>\n",
       "      <th>6</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>27.670772</td>\n",
       "      <td>49.251762</td>\n",
       "      <td>6.366982</td>\n",
       "      <td>-0.044231</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>8.548032</td>\n",
       "      <td>28.952351</td>\n",
       "      <td>3.792775</td>\n",
       "      <td>0.006748</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>-3.300328</td>\n",
       "      <td>15.707746</td>\n",
       "      <td>2.112041</td>\n",
       "      <td>0.205769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>-10.223749</td>\n",
       "      <td>7.352154</td>\n",
       "      <td>1.123517</td>\n",
       "      <td>0.238140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>-13.847281</td>\n",
       "      <td>2.430609</td>\n",
       "      <td>0.634301</td>\n",
       "      <td>-0.028438</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>-0.005858</td>\n",
       "      <td>-0.003273</td>\n",
       "      <td>0.000129</td>\n",
       "      <td>0.000010</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>96</th>\n",
       "      <td>-0.005047</td>\n",
       "      <td>-0.001974</td>\n",
       "      <td>-0.000416</td>\n",
       "      <td>0.001024</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>97</th>\n",
       "      <td>-0.003658</td>\n",
       "      <td>-0.003119</td>\n",
       "      <td>0.000211</td>\n",
       "      <td>-0.000308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>98</th>\n",
       "      <td>-0.002881</td>\n",
       "      <td>-0.002007</td>\n",
       "      <td>-0.000323</td>\n",
       "      <td>0.001012</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99</th>\n",
       "      <td>-0.001723</td>\n",
       "      <td>-0.002852</td>\n",
       "      <td>0.000546</td>\n",
       "      <td>-0.000797</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>92 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            0          2         4         6\n",
       "8   27.670772  49.251762  6.366982 -0.044231\n",
       "9    8.548032  28.952351  3.792775  0.006748\n",
       "10  -3.300328  15.707746  2.112041  0.205769\n",
       "11 -10.223749   7.352154  1.123517  0.238140\n",
       "12 -13.847281   2.430609  0.634301 -0.028438\n",
       "..        ...        ...       ...       ...\n",
       "95  -0.005858  -0.003273  0.000129  0.000010\n",
       "96  -0.005047  -0.001974 -0.000416  0.001024\n",
       "97  -0.003658  -0.003119  0.000211 -0.000308\n",
       "98  -0.002881  -0.002007 -0.000323  0.001012\n",
       "99  -0.001723  -0.002852  0.000546 -0.000797\n",
       "\n",
       "[92 rows x 4 columns]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for i in range (nc):           # loop over all R   \n",
    "    ct = i*ngm                 # extract start point. Since input dataframe is sorted by R and theta, \n",
    "    f = df_inp[2][ct:ct+ngm]   # potentials (V) are extracting for each R value at a time \n",
    "    V_n1 = A_inv.dot(f)        # A-inv * V gives Radial coefficients\n",
    "    V_nf[i,:] = V_n1           # radial coefficients stored in 2D matrix\n",
    "a12 = np.arange(lm)            # creates header for lambda terms\n",
    "df_Vnf = pd.DataFrame(V_nf, columns = a12*sym) # saves final matrix into dataframe with appropriate header\n",
    "df_Vnf[8:]                     # prints first 8 terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "33903c83",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-15.826502499355303"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(df_Vnf[0])   # prints minima for isotropic term"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "578172b6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_dummy = np.arange(2.5,12.5,0.1) # Radial coordinates for plotting data\n",
    "len(x_dummy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2797c66b",
   "metadata": {},
   "source": [
    "#### Plots "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "6c31e6c7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVdUlEQVR4nO3de5CddX3H8c9nN7KEALNAwi0hAQpKg9MgHLBUpbUJcpGRwogDTpEKNTjjhTIdW5DOaKvWWpkBZ+xFFFSmGEUwQqXlkrQVHRtlo8QSInKNJiBshADKzkLYb/84z8LJ5pzdc3v2PM/5vV8zO3vOc55znt/Znfnsb7+/y3FECACQloFeNwAAMPsIfwBIEOEPAAki/AEgQYQ/ACSI8AeABBH+AJAgwh+lZfsrtj+Z02uH7SOaPPfQ7Pw5ebQFyAPhj56yfYftv6tz/EzbvyJQG7N9afYzetb2dbaHet0mlAfhj177iqTzbXvK8fMl3RARO2a/ScVn+xRJl0laLulQSYdL+ttetgnlQvij174taV9Jb5k8YHsfSWdIur7ZF7G9j+3v2B61/Ux2e1HN4/9j+5O2f2D7N7b/3fZ+tm+w/Zzte2wfOuVlT7f9iO1ttj9reyB7rUHbV2bHH5H09iltea/tTbafz55/cas/lCZcIOnaiNgYEc9I+oSkP8vhOuhThD96KiLGJN0o6T01h98l6WcRsaGFlxqQ9GVJSyQtljQm6fNTzjlX1f8oFkr6HUn/mz1nX0mbJH1syvlnSapIOlbSmZIuzI6/T9U/Tm/IHn/nlOc9lT2+t6T3SrrK9rH1Gm37zba3T/P15gbv92hJtT+fDZIOsL1fg/OBnRD+KIKvSjrH9tzs/nuyY02LiF9HxM0R8UJEPC/pU5L+cMppX46IhyPiWUn/KenhiFiTlZa+qWqY1/pMRDwdEb+QdLWk87Lj75J0dUT8MiKelvTpKW25LbtORMR3Jd2pmv9sppz7/YgYnubr+w3e8p6Snq25P3l7rwbnAzsh/NFzWcCNSjrT9uGSjpf0tVZew/Yetr9ge7Pt5yTdLWnY9mDNaU/W3B6rc3/PKS/7y5rbmyUdnN0+uM5jtW05zfY620/b3i7pdEnzW3k/TfiNqv9ZTJq8/XyXr4M+RfijKK5Xtcd/vqQ7I+LJGc6f6i8lvU7SGyNib0knZcenDiS34pCa24slPZ7dfqLOY9WLVWfc3CzpSkkHRMSwpP9o1A7bb8nGIBp91f2PQdJGSctq7i+T9GRE/LqF94eEEf4oiuslrVC1nt5SySezl6q99+2299Wu9ft2fCQbSD5E0iWSvpEdv1HSh20vyganL6t5zm6ShlT9T2aH7dMkva3RBSLiexGx5zRf32vw1OslXWR7adaGv1F15hTQFMIfhRARj0n6gaR5km5t4yWuljRX0jZJ6yTd3oVm3SJpvaR7Jd0m6drs+Bcl3aHqIOuPJX1r8gnZeMOHVf0D8Yykd6u99zOtiLhd0j9K+m9Vy06b1Z0/eEiE+SQvAEhP7qsnbT+m6iDUy5J2REQl72sCAKY3W0vn3xoR22bpWugjtjeqOnd/qosj4obZbg/QL9g3BYUWEUf3ug1AP8q95m/7UVUHvkLSFyLimimPr5S0UpLmzZt33FFHHZVrewCg36xfv35bRCxo5TmzEf4HR8TjtveXdJekD0XE3fXOrVQqMTIykmt7AKDf2F7f6nhq7lM9I+Lx7PtTklZLOiHvawIAppdr+NueZ3uvyduqLna5L89rAgBmlveA7wGSVmdbtc+R9LVscQoAoIdyDf+IeEQ77z8CACgAtncAgAQR/gCQIMIfABJE+ANAggh/AEgQ4Q8ACSL8ASBBhD8AJIjwB4AEEf4AkCDCHwASRPgDQIIIfwBIEOEPAAki/AEgQYQ/ACSI8AeABBH+AJAgwh8AEkT4A0CCCH8ASBDhDwAJIvwBIEGEPwAkiPAHgAQR/gCQIMIfABJE+ANAggh/AEgQ4Q8ACSL8ASBBhD8AJIjwB4AEEf4AkCDCHwASlHv42z7V9gO2H7J9Wd7XA4AiGx0f1z3bt2t0fLxrx7T77nu02o457b6BZtgelPRPkk6WtEXSPbZvjYj787wugPIYHR/XY2NjOnTuXC0YGurrY6u2btVFGzZot4EBvTgxoWuXLZOkjo/p4INf2+rP3RHR6nOaf3H7REkfj4hTsvuXS1JEfLre+ZVKJUZGRnJrD4DG+ikMi3jsqqVLden992tsYuKVn/nutmx3fuzSSxUPPuhWft95h/87JZ0aEX+e3T9f0hsj4oM156yUtFKSFi9efNzmzZtzaw/Q78oUwrmGYQGPDdnabWBAz7/88ivH5g0OSpJ+2+mxAob/OZJOmRL+J0TEh+qdT88fmP0A71UI5xqGBTy25+CgXorQeEF6/rnW/FWt8x9Sc3+RpMdzviZQSLMV4JMBceG9974SDtMdu2TjRu02sPPcj8GBXeeCdPvYawYG9OKUzufLExOy3Z/HIvS57Pf0moEBvTTld9nJsbGJiZfVorx7/nMk/VzScklbJd0j6d0RsbHe+fT8UUbdCvVOeuCd9KJz7ZFOc2zuwMAr77nbYVjUY+ctXJjLOMr+w8ObYmxsqVqQa/hLku3TJV0taVDSdRHxqUbnEv4oujxDvVcB3ssQzisMi3wsD7bXR0SlpefkHf6tIPzRK+303rsd6r0M8F6GMDpH+ANN6Fbvvduh3usAR3kR/sAUefbe8wh1AhztIPyBGs0EfSe997xCHWgV4Y+k1QapJC1Zuzb33juhjiJoJ/zznucP5GKmcs5HjzhCuw0M7BTgdeeVtzD3+ryFC3X2QQftEuoLhoZ2Cfh6x4AiIfxReK3U7SfD/lMPPtj0IptGob5i/vymgh4oI8IfhdZM0Ndbobrb4KA+cvjh+vuHHqL3DtRB+KNQptbts6Xr0wZ9vXLOSxMTunjJEl28ZAm9d6AOwh+F0e26fW3YA9gZ4Y+emamX3426PYD6CH/0RDO9/G7U7QHUxzx/zIqZ5uA3mlu/eflySaJHD0yDef4opE56+dTtgXwQ/ui6dmr5083OAdB9hD+6il4+UA6EP7pmdHycXj5QErt+sCbQgtHxcd2zffsrpZ56K20/esQRmjswoL3nzNHcgYFXevkLhoZ0/PAwwQ/0AD1/tK3e1gsv1pR3JHr5QFHR80dbaks8z+7YobGJCV16//26aulSevlACdDzR9NqZ/FMlnimbr1w7PCwNi9fTi8fKLjShT8fktEbzZZ4Jn8v/G6AYitV2WfV1q1asnatTl63TkvWrtWqrVt73aQktFriAVB8pen515tGeNGGDVoxfz6BkzNKPED/KU34Nwqgx8bGCJ0c1JbXDp07lxIP0GdKE/7TBRC6a2p9/9ply3TtsmW7fKYtoQ+UV2nCf8HQEAE0CxqV1zYvX06JB+gjpQl/STpv4cK6H6qN7pmuvMY8faB/lCr8JT6wIw/N1vcB9I/ShT+6i/o+kCbCP2HU94F0Ef4Jo74PpKtUK3zRXdT3gXQR/omp3X9/cvosWzQA6aHsk5B6g7tMnwXSlFvP3/bHbW+1fW/2dXpe18LM6m3OdtGGDa/8B0CNH0hL3j3/qyLiypyvgSawNxKAWtT8E8HgLoBaeYf/B23/1PZ1tvepd4LtlbZHbI+Mjo7m3Jy0MLgLoBFHRPtPttdIOrDOQ1dIWidpm6SQ9AlJB0XEhdO9XqVSiZGRkbbbg1c1Gtzlk9CA/mN7fURUWnlORzX/iFjRzHm2vyjpO51cC82b6YNvCH0Aec72Oajm7lmS7svrWtjZ5OBurcnBXQCQ8p3t84+2j1G17POYpItzvBZqMLgLYCa5hX9EnJ/Xa2N6fPANgJmwwrdPTB3IZeUugOkQ/n2g0cweBncBNMIir5KbbtsGAGiE8C85ZvYAaAfhX3LM7AHQDsK/5Ni2AUA7GPDtA8zsAdAqwr+E6u3Pw8weAK0g/Eum0bROAGgFNf8SYVongG4h/EuEaZ0AuoXwLxGmdQLoFsK/RJjWCaBbGPAtGaZ1AugGwr/gmNYJIA+Ef4ExrRNAXqj5FxTTOgHkifAvKKZ1AsgT4V9QTOsEkCfCv6CY1gkgTwz4FhjTOgHkhfAvOKZ1AsgDZR8ASBDhXyCj4+O6Z/t2pnMCyB1ln4JgQReA2UTPvwBY0AVgthH+BcCCLgCzjfAvABZ0AZhthH8BsKALwGxjwLcgWNAFYDYR/gXCgi4As4WyDwAkiPDvERZ0Aeglyj49wIIuAL1Gz3+WsaALQBF0FP62z7G90faE7cqUxy63/ZDtB2yf0lkz+wcLugAUQadln/sknS3pC7UHbS+VdK6koyUdLGmN7ddGxMsdXq/0WNAFoAg66vlHxKaIeKDOQ2dK+npEjEfEo5IeknRCJ9eaTpkGT1nQBaAI8hrwXShpXc39LdmxXdheKWmlJC1evLjlC5Vx8JQFXQB6bcbwt71G0oF1HroiIm5p9LQ6x6LeiRFxjaRrJKlSqdQ9p5HawdOxrJRy0YYNWjF/fuEDlQVdAHppxvCPiBVtvO4WSYfU3F8k6fE2Xmdak4OnYzU19MnBU4IVABrLa6rnrZLOtT1k+zBJR0r6UbcvwuApALSn06meZ9neIulESbfZvkOSImKjpBsl3S/pdkkfyGOmT1kGT8s0IA0gDY5oqcyeq0qlEiMjIy0/b3R8vLCDp2UckAZQLrbXR0Rl5jNf1RcrfBcMDen44eHCBT+reQEUVV+Ef1GxmhdAURH+OWJAGkBREf45KsuANID0sKVzzljNC6CICP9ZwGpeAEVD2QcAEkT4A0CCCH8ASBDh32Vs5QCgDBjw7SK2cgBQFvT8u4StHACUCeHfJWzlAKBMCP8uYSsHAGVC+HcJWzkAKBMGfLuIrRwAlAXh32Vs5QCgDCj7AECCCH8ASBDhDwAJIvw7wFYOAMqKAd82sZUDgDKj598GtnIAUHaEfxvYygFA2fVt+OdZj2crBwBl15fhv2rrVi1Zu1Ynr1unJWvXatXWrV19fbZyAFB2fTfgW1uPH8t65xdt2KAV8+d3NZzZygFAmfVd+E/W48dqyjKT9fhuBzRbOQAoq74r+1CPB4CZ9V34U48HgJn1XdlHyqcePzo+Tn0fQN/oy/CXuluPZzUvgH7Td2WfbmM1L4B+RPjPgNW8APpRR+Fv+xzbG21P2K7UHD/U9pjte7Ovf+28qb3B7CEA/ajTnv99ks6WdHedxx6OiGOyr/d3eJ2eYfYQgH7U0YBvRGySJNvdaU1BsZoXQL/Js+Z/mO2f2P6u7bc0Osn2StsjtkdGR0dzbE5nFgwN6fjhYYIfQF+Ysedve42kA+s8dEVE3NLgaU9IWhwRv7Z9nKRv2z46Ip6bemJEXCPpGkmqVCrRfNNb1+xcfeb0A+h3M4Z/RKxo9UUjYlzSeHZ7ve2HJb1W0kjLLeySZufqM6cfQApyKfvYXmB7MLt9uKQjJT2Sx7Wa0excfeb0A0hFp1M9z7K9RdKJkm6zfUf20EmSfmp7g6SbJL0/Ip7urKnta3auPnP6AaSi09k+qyWtrnP8Zkk3d/La3TTdXP3a+j5z+gGkIokVvo3m6q/Ztm2nT/xas20bc/oBJMERuU6waUmlUomRkfzGhGt7+ZK0ZO3anT70Ze7AgDYvXy5JzPYBUBq210dEZeYzX9W3u3rWU7vT5z3btzf8xC/m8wPod0mUfeqhvg8gZcmGP3v2AEhZUmWfqdizB0Cqkg5/qbuf+AUAZZFs2QcAUkb4A0CCCH8ASBDhDwAJIvwBIEGEPwAkiPAHgAQR/gCQIMIfABJE+ANAggh/AEgQ4Q8ACSL8ASBBhD8AJIjwB4AEEf4AkCDCHwASRPgDQIIIfwBIEOEPAAki/AEgQYQ/ACSI8AeABBH+AJAgwh8AEkT4A0CCCH8ASBDhDwAJ6ij8bX/W9s9s/9T2atvDNY9dbvsh2w/YPqXjlgIAuqbTnv9dkl4fEb8n6eeSLpck20slnSvpaEmnSvpn24MdXgsA0CUdhX9E3BkRO7K76yQtym6fKenrETEeEY9KekjSCZ1cCwDQPXO6+FoXSvpGdnuhqn8MJm3Jju3C9kpJK7O7v7H9QBfb1In5krb1uhFdwPsoFt5HsfTL+3hdq0+YMfxtr5F0YJ2HroiIW7JzrpC0Q9INk0+rc37Ue/2IuEbSNU21dhbZHomISq/b0SneR7HwPoqln95Hq8+ZMfwjYsUMF71A0hmSlkfEZMBvkXRIzWmLJD3eauMAAPnodLbPqZL+WtI7IuKFmodulXSu7SHbh0k6UtKPOrkWAKB7Oq35f17SkKS7bEvSuoh4f0RstH2jpPtVLQd9ICJe7vBas61wpag28T6KhfdRLMm+D79aqQEApIIVvgCQIMIfABJE+Ddge9D2T2x/p9dtaZftx2z/n+1725kKVhS2h23flG0lssn2ib1uU6tsvy77PUx+PWf7L3rdrnbYvtT2Rtv32V5le/det6kdti/J3sPGMv0ubF9n+ynb99Uc29f2XbYfzL7vM9PrEP6NXSJpU68b0QVvjYhjSj6X+XOSbo+IoyQtUwl/LxHxQPZ7OEbScZJekLS6t61qne2Fkj4sqRIRr5c0qOpWLqVi+/WS3qfqzgPLJJ1h+8jetqppX1F125xal0laGxFHSlqb3Z8W4V+H7UWS3i7pS71uS+ps7y3pJEnXSlJEvBgR23vaqM4tl/RwRGzudUPaNEfSXNtzJO2hcq7h+V1VZye+kG1R811JZ/W4TU2JiLslPT3l8JmSvprd/qqkP5npdQj/+q6W9FeSJnrcjk6FpDttr8+20SijwyWNSvpyVob7ku15vW5Uh86VtKrXjWhHRGyVdKWkX0h6QtKzEXFnb1vVlvsknWR7P9t7SDpdOy9MLZsDIuIJScq+7z/TEwj/KWyfIempiFjf67Z0wZsi4lhJp0n6gO2Tet2gNsyRdKykf4mIN0j6rZr4l7aobO8m6R2SvtnrtrQjqyWfKekwSQdLmmf7T3vbqtZFxCZJn1F1Z+LbJW1QdU1SMgj/Xb1J0jtsPybp65L+2Pa/9bZJ7YmIx7PvT6laXy7jzqpbJG2JiB9m929S9Y9BWZ0m6ccR8WSvG9KmFZIejYjRiHhJ0rck/UGP29SWiLg2Io6NiJNULaM82Os2deBJ2wdJUvb9qZmeQPhPERGXR8SiiDhU1X/P/ysiStezsT3P9l6TtyW9TdV/dUslIn4l6Ze2J3ctXK7qyvGyOk8lLflkfiHp923v4eqy/uUq4QC8JNneP/u+WNLZKvfv5VZJF2S3L5B0y0xP6OaWziiWAyStzrbdmCPpaxFxe2+b1LYPSbohK5k8Ium9PW5PW7La8smSLu51W9oVET+0fZOkH6taJvmJyrtFws2295P0kqpb0DzT6wY1w/YqSX8kab7tLZI+JukfJN1o+yJV/0CfM+PrsL0DAKSHsg8AJIjwB4AEEf4AkCDCHwASRPgDQIIIfwBIEOEPAAn6f+VkMSHqNYOMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUB0lEQVR4nO3dfZBldX3n8fenZ2QciOyooMIgIAvqopUx2D6tkU0yGIG4EqiYgtQSS9kMqTKRULvZlbBbSe2GzZO1YlV2s2LwqRZxiWQCqwkibq0mlUWZUcYwQwgPDjpgpAkOGu1qGPq7f/RtvPTcme7bt2/f7v69X1Vdfe/vPH1P35nPPed3fvfcVBWSpLaMjboASdLyM/wlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/rVpJPprkt4e07kpy6gLnPbkz//ph1CINg+GvkUry2ST/qUf7eUn+3kDtLck7kuxM8t0k+5L8vn8r9cPw16h9FLg4Sea0XwxcV1UHlr+kVeFI4NeAY4DXAVuBfzvKgrS6GP4atT8Dnge8abYhyXOBtwIfX+hKkjw3yaeTTCT5TufxCV3T/2+S307y10n+Mcn/TvL8JNd1jp7vSHLynNWem+SBJI8m+YMkY511rUvyvk77A8DPzKnlnUnuTvK9zvKX9vtHmU9V/VFV/WVVPVFVDwHXAW9c6u1o7TL8NVJVNQncAPxiV/PPA39bVbv6WNUY8BHgJOBEYBL4wznzXMjMGcVm4J8C/6+zzPOAu4HfnDP/+cA4cAZwHvCuTvsvMfPm9GOd6T83Z7lHOtOPBt4JvD/JGb2KTvLjSfYf5ufHF7j/ZwK7FzivZPhrRfgY8PYkGzvPf7HTtmBV9Q9VdWNV/aCqvgdcBfyLObN9pKrur6rHgb8A7q+q2zpdS3/CTJh3+72qeqyqvgFcDVzUaf954Oqq+mZVPQb8zpxaPtPZTlXVF4Bb6TqzmTPvX1XVpsP8/NV8+57kncy8Cb1vvnmlWYa/Rq4TcBPAeUlOAV4DfKKfdSQ5MskHkzyY5LvAF4FNSdZ1zfbtrseTPZ7/yJzVfrPr8YPA8Z3Hx/eY1l3LOUluT/JYkv3Aucz0zS+5JD8L/C5wTlU9OoxtaG0y/LVSfJyZI/6LgVur6tvzzD/XvwFeBryuqo5mphsEYO6F5H68uOvxicDDncff6jFtZmPJBuBGZo7CX1hVm4A/P1QdSd7UuQZxqJ+eZwydZc8GPgT8y6r6m0Xsnxrm0DCtFB8H/gPwo8Dli1j+Ocwcve9P8jwO7r9fjF9P8iVmzgguA/5rp/0G4D1JPg18H3hv1zJHABuYOZM5kOQc4KeBu3ptoKr+koPPOOaV5KeYuch7flV9ud/lJY/8tSJU1V7gr4GjgJsXsYqrgY3Ao8DtwC1LUNZNwE7gTuAzwLWd9g8BnwV2AV8B/nR2gc71hvcw8wbxHeAXWNz+zOc/Av8E+POus4S/GMJ2tEbFb/KSpPYMvdsnyV7ge8BTwIGqGh/2NiVJh7dcff4/6UgELUaS3cyM3Z/r0qq6brnrkdYKL/hqRauqV4y6BmktGnqff5KvM3Phq4APVtU1c6ZvA7YBHHXUUa9++ctfPtR6JGmt2blz56NVdWw/yyxH+B9fVQ8neQHwOeBXq+qLveYdHx+vHTt2DLUeSVprkuzs93rq0Id6VtXDnd+PANuB1w57m5Kkwxtq+Cc5KslzZh9zmA+7SJKWz7Av+L4Q2N65Vft64BNVtRQfvpEkDWCo4V9VDwBbhrkNSVL/vL2DJDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWrQmg3/iakp7ti/n4mpqVGXIkkrzrC/zGUkrn/oIS7ZtYsjxsZ4Ynqaa7ds4aLNm0ddliStGGvuyH9iaopLdu1icnqaxw8cYHJ6mkt27fIMQJK6rLnw3zs5yRFjz9ytZ42NsXdyckQVSdLKs+bC/+SNG3lievoZbU9OT3Pyxo0jqkiSVp41F/7HbtjAtVu2sHFsjKPXr2fj2BjXbtnCsRs2jLo0SVox1uQF34s2b+asY45h7+QkJ2/caPBL0hxrMvxh5gzA0Jek3tZct48kaX6GvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaNPTwT3J2knuS3JfkvcPentSyiakp7ti//xlfW2rb2m/j2c8+kj4N9ZbOSdYB/w14M7APuCPJzVW1Z5jblZbCxNTUQd8JsZLbrn/oIS7ZtYsjxsZ4Ynqaa7dsAbCtgTaOP/6l8/xzPkiqqt9lFr7y5A3Ab1XVWzrPrwCoqt/pNf/4+Hjt2LFjaPVobWs5SN9/+ulcvmcPk11fYfrshCS2tdB2+eXUvfeGPgw7/H8OOLuq/nXn+cXA66rqV7rm2QZsAzjxxBNf/eCDDw6tHq1Oyx3WqzFINyQcMTbG95566um2o9atA+D7tq39tkWE/7C/yatXMc94t6mqa4BrYObIf5jF9AoRjc5ShXp3WM8G4rvuvPPpgOy37bLduzli7JmXw9aNHXx5bCW1PWtsjCfmHMg9NT1NEtsaauvHsMN/H/DirucnAA8PeZs99QqRizZvHkUpa95yh/pSh/WqDNIqPtD5ez1rbIwn5/xNbVvbbZPT0z88JVigYXf7rAf+DtgKPATcAfxCVe3uNf+w+vwnpqY46fOff8Zp8saxMR7cutUzgAENGuqzBunK+JF163iyiqkl6kLZODb2dI0r6T/4fG0Xbd68Yi4+27a8bS/YtOnumpw8nT4MNfwBkpwLXA2sAz5cVVcdat5hhf8d+/fz5ttv5/EDB55uO3r9em57/et5zaZNS769tWq+oB9VqA8jrFdrkKpNSXZW1Xhfyww7/Pvhkf/KsZigH2WoDyOspdXC8D+M2fCaGxg62GKDftShLrXK8J+HgdFb998FOOgsaaFBb6hLo7GY8B/2aJ8V5dgNG5oPlPm6c37j1FM5YmzsGeHfc/TLIUaXXLR5Mxccd9xBAX7WMccc1Nbr9fA1kpZHU+HfusN158yG/VX33rvgYYSHCnpDXVr5DP81bG53Tmc88GHHxx+xbh2/fsop/Jf77lt00Eta+Qz/NWqx3TlPTk9z6UkncelJJxn00hpm+K8R8x3l99Od0x32ktYmw38NWMhRfr/dOZLWtqaGeq4V8w3NPNTY+ge3bgUw6KU1xqGeDRjkKN/uHEmzDP9VZGJqakF9+Ye7aCtJ4Be4r3jd39m5d3Ky59DM3zj1VDaOjXH0+vVsHBt7+ij/2A0beM2mTQa/pIN45L+C9fpQ1hNd3TvgUb6kxWk+/FfSvWTmG655+Z49Pe+dY1++pH41Hf4r6du9FvqhrDM2beLBrVtXzBuWpNWp2fDvdfH0kl27OOuYY5Y9UPu5kDsb+Ia+pEE0e8G318XTZ42NsXdyclm2P8iFXEkaVLNH/idv3Njz4ulsf/sweSFX0qg1e+R/7IYNXLtly7IfWXd38Tx+4ACT09NPX8h1uKak5dLskT/ARZs39/ySkaXWPYpntovHC7mSRqnp8Ifh36Z4oV08XsiVtJya7fY5nO6LsYOup58uHklaLs0f+c816Nh/u3gkrQaGf5dBx/7bxSNptbDbp0u/Y/+7u4fs4pG0mnjk3+VwY//n3gPI2zFIWs0M/y6zY/8v2bXrGTdOu+3RRw/qzrl8zx5vxyBp1TL855g79h9++DWJs0F/2e7dPW/HcLhvz5KklcTw76H7SP2O/ft7duc8Mee7j70dg6TVxAu+8+h1HeCpKj7g7RgkrWIe+c/jUNcBLtq8mQuOO86jfEmrkuG/AIe6B5AXciWtVob/Ahn0ktYS+/wlqUFDC/8kv5XkoSR3dn7OHda2JEn9GXa3z/ur6n1D3oYkqU92+0hSg4Yd/r+S5GtJPpzkub1mSLItyY4kOyYmJoZcjiQJIDXnk6p9LZzcBryox6QrgduBR4EC/jNwXFW963DrGx8frx07diy6HklqUZKdVTXezzID9flX1VkLmS/Jh4BPD7ItSdLSGeZon+O6np4P3DWsbUmS+jPM0T6/n+RVzHT77AUuHeK2JEl9GFr4V9XFw1q3JGkwDvWUpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaNFD4J3l7kt1JppOMz5l2RZL7ktyT5C2DlSlJWkrrB1z+LuAC4IPdjUlOBy4EXgEcD9yW5KVV9dSA25MkLYGBjvyr6u6quqfHpPOAT1bVVFV9HbgPeO0g25IkLZ1h9flvBr7Z9Xxfp+0gSbYl2ZFkx8TExJDKkSR1m7fbJ8ltwIt6TLqyqm461GI92qrXjFV1DXANwPj4eM95JElLa97wr6qzFrHefcCLu56fADy8iPVIkoZgWN0+NwMXJtmQ5CXAacCXh7QtSVKfBh3qeX6SfcAbgM8k+SxAVe0GbgD2ALcA73akjyStHAMN9ayq7cD2Q0y7CrhqkPVLkobDT/hKUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JatBA4Z/k7Ul2J5lOMt7VfnKSySR3dn7+x+ClSpKWyvoBl78LuAD4YI9p91fVqwZcvyRpCAYK/6q6GyDJ0lQjSVoWw+zzf0mSryb5QpI3HWqmJNuS7EiyY2JiYojlSJJmzXvkn+Q24EU9Jl1ZVTcdYrFvASdW1T8keTXwZ0leUVXfnTtjVV0DXAMwPj5eCy9dkrRY84Z/VZ3V70qragqY6jzemeR+4KXAjr4rlCQtuaF0+yQ5Nsm6zuNTgNOAB4axLUlS/wYd6nl+kn3AG4DPJPlsZ9KZwNeS7AI+BfxyVT02WKmSpKUy6Gif7cD2Hu03AjcOsm5J0vD4CV9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBA4V/kj9I8rdJvpZke5JNXdOuSHJfknuSvGXgSiVJS2bQI//PAa+sqh8F/g64AiDJ6cCFwCuAs4H/nmTdgNuSJC2RgcK/qm6tqgOdp7cDJ3Qenwd8sqqmqurrwH3AawfZliRp6axfwnW9C/hfncebmXkzmLWv03aQJNuAbZ2n/5jkniWsaRDHAI+Ouogl4H6sLO7HyrJW9uNl/S4wb/gnuQ14UY9JV1bVTZ15rgQOANfNLtZj/uq1/qq6BrhmQdUuoyQ7qmp81HUMyv1YWdyPlWUt7Ue/y8wb/lV11jwbfQfwVmBrVc0G/D7gxV2znQA83G9xkqThGHS0z9nAvwfeVlU/6Jp0M3Bhkg1JXgKcBnx5kG1JkpbOoH3+fwhsAD6XBOD2qvrlqtqd5AZgDzPdQe+uqqcG3NZyW3FdUYvkfqws7sfK0ux+5Ic9NZKkVvgJX0lqkOEvSQ0y/A8hybokX03y6VHXslhJ9ib5myR3LmYo2EqRZFOST3VuJXJ3kjeMuqZ+JXlZ53WY/flukl8bdV2LkeTyJLuT3JXk+iTPHnVNi5Hkss4+7F5Nr0WSDyd5JMldXW3PS/K5JPd2fj93vvUY/od2GXD3qItYAj9ZVa9a5WOZPwDcUlUvB7awCl+Xqrqn8zq8Cng18ANg+2ir6l+SzcB7gPGqeiWwjplbuawqSV4J/BIzdx7YArw1yWmjrWrBPsrMbXO6vRf4fFWdBny+8/ywDP8ekpwA/Azwx6OupXVJjgbOBK4FqKonqmr/SIsa3Fbg/qp6cNSFLNJ6YGOS9cCRrM7P8PwzZkYn/qBzi5ovAOePuKYFqaovAo/NaT4P+Fjn8ceAn51vPYZ/b1cD/w6YHnEdgyrg1iQ7O7fRWI1OASaAj3S64f44yVGjLmpAFwLXj7qIxaiqh4D3Ad8AvgU8XlW3jraqRbkLODPJ85McCZzLMz+Yutq8sKq+BdD5/YL5FjD850jyVuCRqto56lqWwBur6gzgHODdSc4cdUGLsB44A/ijqvox4Pss4JR2pUpyBPA24E9GXctidPqSzwNeAhwPHJXkX422qv5V1d3A7zFzZ+JbgF3MfCapGYb/wd4IvC3JXuCTwE8l+Z+jLWlxqurhzu9HmOlfXo13Vt0H7KuqL3Wef4qZN4PV6hzgK1X17VEXskhnAV+vqomqehL4U+Cfj7imRamqa6vqjKo6k5lulHtHXdMAvp3kOIDO70fmW8Dwn6OqrqiqE6rqZGZOz/9PVa26I5skRyV5zuxj4KeZOdVdVarq74FvJpm9a+FWZj45vlpdxCrt8un4BvD6JEdm5mP9W1mFF+ABkryg8/tE4AJW9+tyM/COzuN3ADfNt8BS3tJZK8sLge2d226sBz5RVbeMtqRF+1Xguk6XyQPAO0dcz6J0+pbfDFw66loWq6q+lORTwFeY6Sb5Kqv3Fgk3Jnk+8CQzt6D5zqgLWogk1wM/ARyTZB/wm8DvAjckuYSZN+i3z7seb+8gSe2x20eSGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAb9f4TBsjAyQEigAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAShUlEQVR4nO3df7DldX3f8efr7sq6bKSo+ItFQCtq0ckavDVaI226GJU4UpyagbTGUSZLZmykTKeplM4k00bTNEwl0zQpm+CvKWKIZAPVFhEy1WZSlLvKWhZCXJCVBYFLEDR6Z2G57/5xzpXL3bN799xzzp577uf5mLlzz/l8f3ze373M65zv5/P9fklVIUlqy9S4C5AkHX2GvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4a+JleSTSX5zRPuuJK84wnVP7a6/fhS1SKNg+Gusknwxyb/v0X5OkgcN1OUl+XM/fNQvw1/j9kngvUmypP29wFVVdeDolzQ5kvwzwNBX3wx/jdufAc8D3rLQkOS5wDuBTx/pTpI8N8nnk8wm+V739UmLlv/vJL+Z5C+T/G2S/5Hk+UmuSvL9JLcmOXXJbs9Ock+SR5L8TpKp7r7WJbms234P8PNLanl/kjuT/KC7/YX9/qMc4TH/HeDXgV8bxf61thn+GquqmgOuAX5pUfMvAH9VVbv62NUU8AngFOBkYA74vSXrnEfnjGIz8HeB/9vd5nnAnXSCdLFzgWngDOAc4APd9l+m8+H0U93l/3TJdg93lx8HvB/4WJIzehWd5GeSPHaYn585zDF/FPgD4MHDrCP1ZPhrNfgU8J4kG7vvf6nbdsSq6m+q6tqq+lFV/QD4CPAPl6z2iaq6u6oeB/4XcHdV3dQdWvoTOmG+2G9X1aNV9R3gcuD8bvsvAJdX1X1V9SjwW0tq+UK3n6qqLwM3sujMZsm6f1FVxx/m5y96bZdkGngz8F+O7F9IeibDX2PXDbhZ4JwkLwf+PvCZfvaR5NgkVyTZm+T7wFeA45OsW7TaQ4tez/V4/xNLdnvfotd7gRO7r0/ssWxxLe9IckuSR5M8BpwNnNDP8RxOd/jp94GLnBPRShn+Wi0+Tecb/3uBG6vqoWXWX+pfAa8CfrqqjgPO7LYvnUjux0sXvT4ZeKD7+rs9lnU6SzYA1wKXAS+qquOB/3moOpK8pTsHcaifXmcMx9EZbvrjJA8Ct3bb9x1ifekgXiWg1eLTwL8DfhK4eAXbP4fOt/fHkjyPg8fvV+JfJ/kqnTOCi4D/3G2/BvhQks8DPwQ+vGibY4ANdM5kDiR5B/BzwO29Oqiq/8PBZxzLeZynz0Kg80H0NeD13X6lZfnNX6tCVd0L/CWwCbh+Bbu4HNgIPALcAtwwhLKuA3YCtwFfAK7stv8h8EVgF/B14E8XNujON3yIzgfE94BfZGXHc0jduYQHF354OvAfqqonhtmX1q74f/KSpPaMfNgnyb3AD4CngANVNT3qPiVJh3e0xvx/tqoeOUp9aQ1JspvOtftLXVhVVx3teqS1wglfrWpV9Zpx1yCtRSMf80/ybToTXwVcUVXblyzfBmwD2LRp0+tf/epXj7QeSVprdu7c+UhVvaCfbY5G+J9YVQ8keSHwJeBXq+orvdadnp6umZmZkdYjSWtNkp39zqeO/FLPqnqg+/thYAfwhlH3KUk6vJGGf5JNSZ6z8JrD3OwiSTp6Rj3h+yJgR/dR7euBz1TVMG6+kSQNYKThX1X3AFtG2YckqX8+3kGSGmT4S1KDVnX4z+7fz62PPcbs/v3jLkWS1pRVe4fv1fffzwW7dnHM1BRPzM9z5ZYtnL9587jLkqQ1YVV+85/dv58Ldu1ibn6exw8cYG5+ngt27fIMQJKGZFWG/71zcxwz9czSnjU1xb1zc2OqSJLWllUZ/qdu3MgT8/PPaHtyfp5TN248xBaSpH6syvB/wYYNXLllCxunpjhu/Xo2Tk1x5ZYtvGDDhnGXJklrwqqd8D1/82bOOuEE7p2b49SNGw1+SRqiVRv+0DkDMPQlafhW5bCPJGm0DH9JapDhL0kNMvwlqUGGvyQ1yPCXpAZNXPj7pE9JGtyqvs5/KZ/0KUnDMTHf/H3SpyQNz8SEv0/6lKThmZjw90mfkjQ8ExP+PulTkoZnoiZ8fdKnJA3HRIU/+KRPSRqGiRn2kSQNz5oIf2/8kqT+TNywz1Le+CVJ/Zvob/7e+CVJKzPR4e+NX5K0MhMd/oe78ct5AEk6tIkO/0Pd+HXTI49wys0389ZbbuGUm2/m6vvvH3epkrSqTPyE79IbvwBOuflm5ubnmeueFVywaxdnnXACgDeISRJrIPzhmTd+3frYYxwzNfXj4IfOPMAVe/fy0T17vCpIkpjwYZ9ees0DPPHUU3x0z56eVwU5NyCpRWsu/HvNA1x62mk9rwq6Yu9e5wYkNWnk4Z/k7UnuSrInyYdH3R905gH2bt3KTW98I3u3buXCU07p62xgEvU6g7HNNtvaaOPZzz6WPo10zD/JOuC/Am8F9gG3Jrm+qu4YZb9w8APgrtyyhQt27eJZU1M8OT/Pv33FK7jsnnsOmhtYuEdg1BPDs/v3H9THStt63eUM2GabbY20ceKJr6RPqap+tznynSdvAn6jqt7WfX8JQFX9Vq/1p6ena2ZmZmT1LA5NePqqoAUbp6b42Omnc/Edd6x4Yvhoh/VCvYuP49kJSWyzzbZW2i6+mPrWt0IfRn21z2bgvkXv9wE/vXiFJNuAbQAnn3zySItZ7mxgcZAeyWWiKwn1Xn184LbbfvxH7Lftot27D5rPWDd18GiebbbZtvbb+jHq8O/1SfSMU42q2g5sh843/xHX8wxL7xFYeFzEkVwmCisL9WGH9bOmpnhiydnbU/PzJLHNNtsaa+vHqMN/H/DSRe9PAh4YcZ99WXo2cLiJ4WGE+tDDuorf7X7wLJzBLP5wss0229Z+29z8/FP0adRj/uuBvwa2AvcDtwK/WFW7e60/6jH/I7EwdLN0YvjxAwd+vM6mdesA+OFTT/97/8S6dTxZxf75w4/TLZ5XGNZ/AOdv3jzUCWTbbLNtstpeePzxd9bc3On0YaThD5DkbOByYB3w8ar6yKHWXQ3hD8tPDA8a6qMIa0ntSrKzqqb72mbU4d+P1RL+Sy09GxhGqEvSsBj+I2SoS1qtVhL+a+LBbkfD0onhQ7VJ0iRYc8/2kSQtz/CXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBIwv/JL+R5P4kt3V/zh5VX5Kk/qwf8f4/VlWXjbgPSVKfHPaRpAaNOvz/RZJvJvl4kuf2WiHJtiQzSWZmZ2dHXI4kCSBVtfKNk5uAF/dYdClwC/AIUMB/AF5SVR843P6mp6drZmZmxfVIUouS7Kyq6X62GWjMv6rOOpL1kvwh8PlB+pIkDc8or/Z5yaK35wK3j6ovSVJ/Rnm1z39K8jo6wz73AheOsC9JUh9GFv5V9d5R7VuSNBgv9ZSkBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkho0UPgneU+S3Unmk0wvWXZJkj1J7krytsHKlCQN0/oBt78deDdwxeLGJKcD5wGvAU4Ebkryyqp6asD+JElDMNA3/6q6s6ru6rHoHOCzVbW/qr4N7AHeMEhfkqThGdWY/2bgvkXv93XbDpJkW5KZJDOzs7MjKkeStNiywz5JbgJe3GPRpVV13aE269FWvVasqu3AdoDp6eme60iShmvZ8K+qs1aw333ASxe9Pwl4YAX7kSSNwKiGfa4HzkuyIcnLgNOAr42oL0lSnwa91PPcJPuANwFfSPJFgKraDVwD3AHcAHzQK30kafUY6FLPqtoB7DjEso8AHxlk/5Kk0fAOX0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNGij8k7wnye4k80mmF7WfmmQuyW3dn/82eKmSpGFZP+D2twPvBq7osezuqnrdgPuXJI3AQOFfVXcCJBlONZKko2KUY/4vS/KNJF9O8pZDrZRkW5KZJDOzs7MjLEeStGDZb/5JbgJe3GPRpVV13SE2+y5wclX9TZLXA3+W5DVV9f2lK1bVdmA7wPT0dB156ZKklVo2/KvqrH53WlX7gf3d1zuT3A28Epjpu0JJ0tCNZNgnyQuSrOu+fjlwGnDPKPqSJPVv0Es9z02yD3gT8IUkX+wuOhP4ZpJdwOeAX6mqRwcrVZI0LINe7bMD2NGj/Vrg2kH2LUkaHe/wlaQGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGjRQ+Cf5nSR/leSbSXYkOX7RskuS7ElyV5K3DVypJGloBv3m/yXgtVX1k8BfA5cAJDkdOA94DfB24PeTrBuwL0nSkAwU/lV1Y1Ud6L69BTip+/oc4LNVtb+qvg3sAd4wSF+SpOFZP8R9fQD44+7rzXQ+DBbs67YdJMk2YFv37d8muWuINQ3iBOCRcRcxBB7H6uJxrC5r5The1e8Gy4Z/kpuAF/dYdGlVXddd51LgAHDVwmY91q9e+6+q7cD2I6r2KEoyU1XT465jUB7H6uJxrC5r6Tj63WbZ8K+qs5bp9H3AO4GtVbUQ8PuAly5a7STggX6LkySNxqBX+7wd+DfAu6rqR4sWXQ+cl2RDkpcBpwFfG6QvSdLwDDrm/3vABuBLSQBuqapfqardSa4B7qAzHPTBqnpqwL6OtlU3FLVCHsfq4nGsLs0eR54eqZEktcI7fCWpQYa/JDXI8D+EJOuSfCPJ58ddy0oluTfJ/0ty20ouBVstkhyf5HPdR4ncmeRN466pX0le1f07LPx8P8m/HHddK5Hk4iS7k9ye5Ookzx53TSuR5KLuMeyepL9Fko8neTjJ7YvanpfkS0m+1f393OX2Y/gf2kXAneMuYgh+tqpeN+HXMv8ucENVvRrYwgT+Xarqru7f4XXA64EfATvGW1X/kmwGPgRMV9VrgXV0HuUyUZK8FvhlOk8e2AK8M8lp463qiH2SzmNzFvswcHNVnQbc3H1/WIZ/D0lOAn4e+KNx19K6JMcBZwJXAlTVE1X12FiLGtxW4O6q2jvuQlZoPbAxyXrgWCbzHp6/R+fqxB91H1HzZeDcMdd0RKrqK8CjS5rPAT7Vff0p4J8stx/Dv7fLgV8D5sdcx6AKuDHJzu5jNCbRy4FZ4BPdYbg/SrJp3EUN6Dzg6nEXsRJVdT9wGfAd4LvA41V143irWpHbgTOTPD/JscDZPPPG1Enzoqr6LkD39wuX28DwXyLJO4GHq2rnuGsZgjdX1RnAO4APJjlz3AWtwHrgDOAPquqngB9yBKe0q1WSY4B3AX8y7lpWojuWfA7wMuBEYFOSfz7eqvpXVXcCv03nycQ3ALvo3JPUDMP/YG8G3pXkXuCzwD9O8t/HW9LKVNUD3d8P0xlfnsQnq+4D9lXVV7vvP0fnw2BSvQP4elU9NO5CVugs4NtVNVtVTwJ/CvyDMde0IlV1ZVWdUVVn0hlG+da4axrAQ0leAtD9/fByGxj+S1TVJVV1UlWdSuf0/M+rauK+2STZlOQ5C6+Bn6NzqjtRqupB4L4kC08t3ErnzvFJdT4TOuTT9R3gjUmOTee2/q1M4AQ8QJIXdn+fDLybyf67XA+8r/v6fcB1y20wzEc6a3V5EbCj+9iN9cBnquqG8Za0Yr8KXNUdMrkHeP+Y61mR7tjyW4ELx13LSlXVV5N8Dvg6nWGSbzC5j0i4NsnzgSfpPILme+Mu6EgkuRr4R8AJSfYBvw78R+CaJBfQ+YB+z7L78fEOktQeh30kqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWrQ/wcw9UFi8+zjrgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAASfUlEQVR4nO3df7DddX3n8efrJhID1YkKiiQYZMEf6G4Uby1uK91uYlXqyOKsHejUuuo0dMat1tnprpQ/2tktdtsyhe5025IWf42IpdIU1h+IpFOt46YlUWAJlBIwkYCVSzFqMRsMee8f53vZm5tzc3PvuSfn3nyej5kz55z399f73My8zvf7Od/vN6kqJEltGRt1A5KkY8/wl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/LVkJflokt8c0roryVlHOe8Z3fzLh9GLNAyGv0YqyReS/Nc+9QuT/KOBOrMkZyb5TJLvJ3ksye+MuictHYa/Ru2jwDuSZFr9HcB1VXXg2Le0+CU5Afgi8FfAqcAa4BMjbUpLiuGvUftL4LnA6ycLSZ4DvAX4+NGuJMlzur3giSTf6V6vmTL9r5P8ZpKvJvnnJP8ryfOSXJfke0luT3LGtNVekOTBbq/6d5OMdetaluTKrv4g8DPTenlXknu7PfIHk1w61z/KUfgPwCNV9XtV9URV/d+qumsI29FxyvDXSFXVPuAG4BemlH8W+PuqunMOqxoDPgKsBV4E7AP+YNo8F9M7olgN/Avgf3fLPBe4F/j1afNfBIwD5wIXAu/u6r9I78vp1d30fz9tuUe76c8G3gVcleTcfk0n+Ykke4/w+IkZPu95wK4kn+++hP46yb+cYV7pMIa/FoOPAW9PsrJ7/wtd7ahV1T9V1Y1V9YOq+j5wBfCT02b7SFU9UFXfBT4PPFBVt3VDS39OL8yn+u2qeryqvglcDVzS1X8WuLqqHqqqx4HfmtbLZ7vtVFV9CbiVKUc20+b9SlWtOsLjKzN85DX0vsz+B3Aa8Fngpm44SJqV4a+R6wJuArgwyZnAjwKfnMs6kpyY5Joku5N8D/gysCrJsimzfXvK63193v/ItNU+NOX1bnohS/c8fdrUXt6cZGuSx5PsBS4ATp7L5zkK+4CvVNXnq+pJ4ErgecDLF3g7Ok4Z/losPk5vj/8dwK1V9e1Z5p/uPwEvBX6sqp4NnN/Vp/+QPBenT3n9IuCR7vW3+kzrbSxZAdxIL4xfUFWrgM/N1EeS13e/Qcz06HvEANwFeD92zZvhr8Xi48AGeuPpcxry6TyL3t7w3iTP5fDx+/n41e6H5NOB9wN/1tVvAN6XZE334/QHpyxzArCC3pHMgSRvBn56pg1U1d9U1Y8c4fE3Myz6CeC8JBu6o5tfAR6j99uFNCvDX4tCVe0CvgqcBNw8j1VcDaykF4BbgVsWoK2bgO3AHfTG1K/t6n8CfAG4E/ga8BeTC3S/N7yP3hfEd4CfY36f54iq6j7g54E/7rZzIfDWbghImlX8n7wkqT1Dv3oyyS7g+8BTwIGqGh/2NiVJR3asLp3/qap67BhtS8eRJDvonbs/3aVVdd2x7kc6XnjfFC1qVfWKUfcgHY+GPuaf5Bv0fpAq4Jqq2jRt+kZgI8BJJ530mpe97GVD7UeSjjfbt29/rKpOmcsyxyL8T6uqR5I8n96NqH65qr7cb97x8fHatm3bUPuRpONNku1z/T116Kd6VtUj3fOjwGbgtcPepiTpyIYa/klOSvKsydf0Lna5e5jblCTNbtg/+L4A2Nzdqn058MmqWoiLbyRJAxhq+FfVg8C6YW5DkjR33t5Bkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDv4+J/fu5fe9eJvbvH3UrkjQUy0fdwFxN7N/Prn37OGPlSk5ZsWLG2nxd//DDvOfOOzlhbIwnDx7k2nXruGT16oVoXZIWjSUV/v2CGRgorKd+cUyua9/Bg+w7ePDp9xtOPnngLxVJWkyGPuyT5E1J7kuyM8kHjzTvE089dchQy9Thl4n9+58O5u8eOMC+gwd59x13HFZ7z513Pj3/bEM31z/8MGu3bOENW7eydssWrtm9mxPGDv2TPGNsjF379vVdvt82rFmzZu1Y13jmM09kjoa6559kGfA/gTcAe4Dbk9xcVff0m//+J55g7ZYtfffof+2sszhhbOzpPXKAZWOHf3c9Y2yMa3bv5kM7dx7xaGDql8nkOq+4/36SHDLfDw8e5IyVKw8bWjraoxBr1qxZG3aN0057CXOUqprrMke/8uR1wG9U1Ru795cBVNVv9Z3/7LOLq67imQlJDgn6QWorx8bYvX79IUM3t+/dyxu2buW7Bw48XXv28uX86pln8qGdO3nG2Bg/nOGPftU55/CBe+5ZsP6sWbNmbaDaBz5A3X//oXuusxj2mP9q4KEp7/cAPzZ1hiQbgY0AnHIK0H+P/oRly44YzJO1XzvrLK588MFD/khTh24m997PWLmSJ6fMA729/EvXruXStWsP+R1g7ZYthxwhvH/HjsOGh/r1bM2aNWvHsjYXww7/ft9EhxxqVNUmYBN0e/7AUwcP9h1+mR7Mk3vyG04++ZCw/tDOnYct+7W9e/nJr371kEOma9etO+SL49p1655e5+Tz7Xv3Hjbc9IyxMZ6cdsTUr2dr1qxZO5a1uRh2+O8BTp/yfg3wyEwzjyWsGBvru0ffL5gnnbJixSG16aE+dZhm6lk8u9evZ/f69Uc8TbTfEcJTVfx+t84jHYVYs2bN2rGo7Tt48CnmaNhj/suBfwDWAw8DtwM/V1U7+s3/8le/ur68deuCnL8/ddld+/b1Hd+/7bzz+NFVq2Zd1+SPu1P/6JesXt23P2vWrFk71rXnr1p1b+3bdw5zMNTwB0hyAXA1sAz4cFVdMdO84+PjtW3btgXvYWL//qfH7Sf1+xF4tnUs1IVkkrSQkmyvqvG5LDP0i7yq6nPA54a9nSM5ZcWKI47vH+06DH1Jx4sldYXvIC5ZvfqQH4YNckktayb8wb13SZrkXT0lqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDRpa+Cf5jSQPJ7mje1wwrG1JkuZm+ZDXf1VVXTnkbUiS5shhH0lq0LDD/z8muSvJh5M8p98MSTYm2ZZk28TExJDbkSQBpKrmv3ByG3Bqn0mXA1uBx4AC/hvwwqp695HWNz4+Xtu2bZt3P5LUoiTbq2p8LssMNOZfVRuOZr4kfwJ8ZpBtSZIWzjDP9nnhlLcXAXcPa1uSpLkZ5tk+v5PkVfSGfXYBlw5xW5KkORha+FfVO4a1bknSYDzVU5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JatBA4Z/k7Ul2JDmYZHzatMuS7ExyX5I3DtamJGkhLR9w+buBtwHXTC0mOQe4GHgFcBpwW5KXVNVTA25PkrQABtrzr6p7q+q+PpMuBD5VVfur6hvATuC1g2xLkrRwhjXmvxp4aMr7PV3tMEk2JtmWZNvExMSQ2pEkTTXrsE+S24BT+0y6vKpummmxPrXqN2NVbQI2AYyPj/edR5K0sGYN/6raMI/17gFOn/J+DfDIPNYjSRqCYQ373AxcnGRFkhcDZwN/N6RtSZLmaNBTPS9Ksgd4HfDZJF8AqKodwA3APcAtwHs900eSFo+BTvWsqs3A5hmmXQFcMcj6JUnD4RW+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkho0UPgneXuSHUkOJhmfUj8jyb4kd3SPPx68VUnSQlk+4PJ3A28Drukz7YGqetWA65ckDcFA4V9V9wIkWZhuJEnHxDDH/F+c5OtJvpTk9TPNlGRjkm1Jtk1MTAyxHUnSpFn3/JPcBpzaZ9LlVXXTDIt9C3hRVf1TktcAf5nkFVX1vekzVtUmYBPA+Ph4HX3rkqT5mjX8q2rDXFdaVfuB/d3r7UkeAF4CbJtzh5KkBTeUYZ8kpyRZ1r0+EzgbeHAY25Ikzd2gp3pelGQP8Drgs0m+0E06H7gryZ3Ap4FfqqrHB2tVkrRQBj3bZzOwuU/9RuDGQdYtSRoer/CVpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaNFD4J/ndJH+f5K4km5OsmjLtsiQ7k9yX5I0DdypJWjCD7vl/EXhlVf0r4B+AywCSnANcDLwCeBPwh0mWDbgtSdICGSj8q+rWqjrQvd0KrOleXwh8qqr2V9U3gJ3AawfZliRp4SxfwHW9G/iz7vVqel8Gk/Z0tcMk2Qhs7N7+c5L7FrCnQZwMPDbqJhaAn2Nx8XMsLsfL53jpXBeYNfyT3Aac2mfS5VV1UzfP5cAB4LrJxfrMX/3WX1WbgE1H1e0xlGRbVY2Puo9B+TkWFz/H4nI8fY65LjNr+FfVhlk2+k7gLcD6qpoM+D3A6VNmWwM8MtfmJEnDMejZPm8C/gvw1qr6wZRJNwMXJ1mR5MXA2cDfDbItSdLCGXTM/w+AFcAXkwBsrapfqqodSW4A7qE3HPTeqnpqwG0da4tuKGqe/ByLi59jcWn2c+T/j9RIklrhFb6S1CDDX5IaZPjPIMmyJF9P8plR9zJfSXYl+T9J7pjPqWCLRZJVST7d3Urk3iSvG3VPc5Xkpd2/w+Tje0l+ZdR9zUeSDyTZkeTuJNcneeaoe5qPJO/vPsOOpfRvkeTDSR5NcveU2nOTfDHJ/d3zc2Zbj+E/s/cD9466iQXwU1X1qiV+LvPvA7dU1cuAdSzBf5equq/7d3gV8BrgB8Dm0XY1d0lWA+8DxqvqlcAyerdyWVKSvBL4RXp3HlgHvCXJ2aPt6qh9lN5tc6b6ILClqs4GtnTvj8jw7yPJGuBngD8ddS+tS/Js4HzgWoCqerKq9o60qcGtBx6oqt2jbmSelgMrkywHTmRpXsPzcnpnJ/6gu0XNl4CLRtzTUamqLwOPTytfCHyse/0x4N/Nth7Dv7+rgf8MHBxxH4Mq4NYk27vbaCxFZwITwEe6Ybg/TXLSqJsa0MXA9aNuYj6q6mHgSuCbwLeA71bVraPtal7uBs5P8rwkJwIXcOiFqUvNC6rqWwDd8/NnW8DwnybJW4BHq2r7qHtZAD9eVecCbwbem+T8UTc0D8uBc4E/qqpXA09wFIe0i1WSE4C3An8+6l7moxtLvhB4MXAacFKSnx9tV3NXVfcCv03vzsS3AHfSuyapGYb/4X4ceGuSXcCngH+b5BOjbWl+quqR7vlReuPLS/HOqnuAPVX1t937T9P7Mliq3gx8raq+PepG5mkD8I2qmqiqHwJ/AfzrEfc0L1V1bVWdW1Xn0xtGuX/UPQ3g20leCNA9PzrbAob/NFV1WVWtqaoz6B2e/1VVLbk9myQnJXnW5Gvgp+kd6i4pVfWPwENJJu9auJ7eleNL1SUs0SGfzjeB85KcmN5l/etZgj/AAyR5fvf8IuBtLO1/l5uBd3av3wncNNsCC3lLZy0uLwA2d7fdWA58sqpuGW1L8/bLwHXdkMmDwLtG3M+8dGPLbwAuHXUv81VVf5vk08DX6A2TfJ2le4uEG5M8D/ghvVvQfGfUDR2NJNcD/wY4Ocke4NeB/w7ckOQ99L6g3z7rery9gyS1x2EfSWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5Ia9P8A5PoRRtMOsf8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot raw data\n",
    "for i in range(0,lm):\n",
    "    y_dummy = df_Vnf[i*2]\n",
    "    # Plot the noisy exponential data\n",
    "    plt.scatter(x_dummy, y_dummy,s=20, color='#00b3b3', label='Data')\n",
    "    plt.title(\"V_lambda = %d\" %(i*2))\n",
    "    plt.ylim(-20, 5)\n",
    "    plt.xlim(3.3, 10)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff114e3c",
   "metadata": {},
   "source": [
    "#### save data for further use \n",
    "    1) either directly into MOLSCAT\n",
    "    2) OR curve fit into appropriate functions and then use in MOLSCAT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "69806fca",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_Vnf.to_csv('2D_Vlam.dat', index=None, header=True,\n",
    "                sep=',')    # save V_lam coefficients to file separated by comma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c925a01",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}