FDTDEquations2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "source": [
    "from sympy import *\r\n",
    "import sympy\r\n",
    "import pyperclip\r\n",
    "\r\n",
    "def easycopy(sympyObj):\r\n",
    "    pyperclip.copy((latex(sympyObj)).replace(\"{\\\\left(x,y,z,t \\\\right)}\",\"\"))\r\n",
    "\r\n",
    "\r\n",
    "def gradient(scalar_function, variables):\r\n",
    "    matrix_scalar_function = Matrix([scalar_function])\r\n",
    "    return matrix_scalar_function.jacobian(variables)\r\n",
    "def tensorDiv(mat, axes):\r\n",
    "    res = [0,0,0]\r\n",
    "    for i in range(3):\r\n",
    "        for j in range(3):\r\n",
    "            res[i] += diff(mat[i,j],axes[j])\r\n",
    "    return Matrix(res)\r\n",
    "\r\n",
    "mu, lam = symbols(\"\\mu \\lambda\")\r\n",
    "\r\n",
    "c = Matrix([\r\n",
    "[2*mu+lam,lam,lam,0,0,0],\r\n",
    "[lam,2*mu+lam,lam,0,0,0],\r\n",
    "[lam,lam,2*mu+lam,0,0,0],\r\n",
    "[0,0,0,mu,0,0],\r\n",
    "[0,0,0,0,mu,0],\r\n",
    "[0,0,0,0,0,mu],\r\n",
    "])\r\n",
    "\r\n",
    "\r\n",
    "x,y,z,t = symbols(\"x y z t\")\r\n",
    "dims = [x,y,z]\r\n",
    "us = [Function(\"U_x\")(x,y,z,t),Function(\"U_y\")(x,y,z,t),Function(\"U_z\")(x,y,z,t)]\r\n",
    "\r\n",
    "U=Matrix(us)\r\n",
    "easycopy(c)\r\n",
    "U"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\operatorname{U_{x}}{\\left(x,y,z,t \\right)}\\\\\\operatorname{U_{y}}{\\left(x,y,z,t \\right)}\\\\\\operatorname{U_{z}}{\\left(x,y,z,t \\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[U_x(x, y, z, t)],\n",
       "[U_y(x, y, z, t)],\n",
       "[U_z(x, y, z, t)]])"
      ]
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "source": [
    "gradientMat = gradient(U,symbols(\"x y z\"))\r\n",
    "gradientTranspose = gradient(U.transpose() ,symbols(\"x y z\"))\r\n",
    "uTerm = (1/2)*(gradientMat + gradientTranspose)\r\n",
    "easycopy(uTerm)\r\n",
    "uTerm"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1.0 \\frac{\\partial}{\\partial x} \\operatorname{U_{x}}{\\left(x,y,z,t \\right)} & 1.0 \\frac{\\partial}{\\partial y} \\operatorname{U_{x}}{\\left(x,y,z,t \\right)} & 1.0 \\frac{\\partial}{\\partial z} \\operatorname{U_{x}}{\\left(x,y,z,t \\right)}\\\\1.0 \\frac{\\partial}{\\partial x} \\operatorname{U_{y}}{\\left(x,y,z,t \\right)} & 1.0 \\frac{\\partial}{\\partial y} \\operatorname{U_{y}}{\\left(x,y,z,t \\right)} & 1.0 \\frac{\\partial}{\\partial z} \\operatorname{U_{y}}{\\left(x,y,z,t \\right)}\\\\1.0 \\frac{\\partial}{\\partial x} \\operatorname{U_{z}}{\\left(x,y,z,t \\right)} & 1.0 \\frac{\\partial}{\\partial y} \\operatorname{U_{z}}{\\left(x,y,z,t \\right)} & 1.0 \\frac{\\partial}{\\partial z} \\operatorname{U_{z}}{\\left(x,y,z,t \\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1.0*Derivative(U_x(x, y, z, t), x), 1.0*Derivative(U_x(x, y, z, t), y), 1.0*Derivative(U_x(x, y, z, t), z)],\n",
       "[1.0*Derivative(U_y(x, y, z, t), x), 1.0*Derivative(U_y(x, y, z, t), y), 1.0*Derivative(U_y(x, y, z, t), z)],\n",
       "[1.0*Derivative(U_z(x, y, z, t), x), 1.0*Derivative(U_z(x, y, z, t), y), 1.0*Derivative(U_z(x, y, z, t), z)]])"
      ]
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "source": [
    "result = tensorDiv(uTerm*c[0:3,0:3],dims)\r\n",
    "result"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2.0 \\lambda \\frac{\\partial^{2}}{\\partial y\\partial x} \\operatorname{U_{x}}{\\left(x,y,z,t \\right)} + 2.0 \\lambda \\frac{\\partial^{2}}{\\partial z\\partial x} \\operatorname{U_{x}}{\\left(x,y,z,t \\right)} + 2.0 \\lambda \\frac{\\partial^{2}}{\\partial z\\partial y} \\operatorname{U_{x}}{\\left(x,y,z,t \\right)} + \\left(1.0 \\lambda + 2.0 \\mu\\right) \\frac{\\partial^{2}}{\\partial x^{2}} \\operatorname{U_{x}}{\\left(x,y,z,t \\right)} + \\left(1.0 \\lambda + 2.0 \\mu\\right) \\frac{\\partial^{2}}{\\partial y^{2}} \\operatorname{U_{x}}{\\left(x,y,z,t \\right)} + \\left(1.0 \\lambda + 2.0 \\mu\\right) \\frac{\\partial^{2}}{\\partial z^{2}} \\operatorname{U_{x}}{\\left(x,y,z,t \\right)}\\\\2.0 \\lambda \\frac{\\partial^{2}}{\\partial y\\partial x} \\operatorname{U_{y}}{\\left(x,y,z,t \\right)} + 2.0 \\lambda \\frac{\\partial^{2}}{\\partial z\\partial x} \\operatorname{U_{y}}{\\left(x,y,z,t \\right)} + 2.0 \\lambda \\frac{\\partial^{2}}{\\partial z\\partial y} \\operatorname{U_{y}}{\\left(x,y,z,t \\right)} + \\left(1.0 \\lambda + 2.0 \\mu\\right) \\frac{\\partial^{2}}{\\partial x^{2}} \\operatorname{U_{y}}{\\left(x,y,z,t \\right)} + \\left(1.0 \\lambda + 2.0 \\mu\\right) \\frac{\\partial^{2}}{\\partial y^{2}} \\operatorname{U_{y}}{\\left(x,y,z,t \\right)} + \\left(1.0 \\lambda + 2.0 \\mu\\right) \\frac{\\partial^{2}}{\\partial z^{2}} \\operatorname{U_{y}}{\\left(x,y,z,t \\right)}\\\\2.0 \\lambda \\frac{\\partial^{2}}{\\partial y\\partial x} \\operatorname{U_{z}}{\\left(x,y,z,t \\right)} + 2.0 \\lambda \\frac{\\partial^{2}}{\\partial z\\partial x} \\operatorname{U_{z}}{\\left(x,y,z,t \\right)} + 2.0 \\lambda \\frac{\\partial^{2}}{\\partial z\\partial y} \\operatorname{U_{z}}{\\left(x,y,z,t \\right)} + \\left(1.0 \\lambda + 2.0 \\mu\\right) \\frac{\\partial^{2}}{\\partial x^{2}} \\operatorname{U_{z}}{\\left(x,y,z,t \\right)} + \\left(1.0 \\lambda + 2.0 \\mu\\right) \\frac{\\partial^{2}}{\\partial y^{2}} \\operatorname{U_{z}}{\\left(x,y,z,t \\right)} + \\left(1.0 \\lambda + 2.0 \\mu\\right) \\frac{\\partial^{2}}{\\partial z^{2}} \\operatorname{U_{z}}{\\left(x,y,z,t \\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[2.0*\\lambda*Derivative(U_x(x, y, z, t), x, y) + 2.0*\\lambda*Derivative(U_x(x, y, z, t), x, z) + 2.0*\\lambda*Derivative(U_x(x, y, z, t), y, z) + (1.0*\\lambda + 2.0*\\mu)*Derivative(U_x(x, y, z, t), (x, 2)) + (1.0*\\lambda + 2.0*\\mu)*Derivative(U_x(x, y, z, t), (y, 2)) + (1.0*\\lambda + 2.0*\\mu)*Derivative(U_x(x, y, z, t), (z, 2))],\n",
       "[2.0*\\lambda*Derivative(U_y(x, y, z, t), x, y) + 2.0*\\lambda*Derivative(U_y(x, y, z, t), x, z) + 2.0*\\lambda*Derivative(U_y(x, y, z, t), y, z) + (1.0*\\lambda + 2.0*\\mu)*Derivative(U_y(x, y, z, t), (x, 2)) + (1.0*\\lambda + 2.0*\\mu)*Derivative(U_y(x, y, z, t), (y, 2)) + (1.0*\\lambda + 2.0*\\mu)*Derivative(U_y(x, y, z, t), (z, 2))],\n",
       "[2.0*\\lambda*Derivative(U_z(x, y, z, t), x, y) + 2.0*\\lambda*Derivative(U_z(x, y, z, t), x, z) + 2.0*\\lambda*Derivative(U_z(x, y, z, t), y, z) + (1.0*\\lambda + 2.0*\\mu)*Derivative(U_z(x, y, z, t), (x, 2)) + (1.0*\\lambda + 2.0*\\mu)*Derivative(U_z(x, y, z, t), (y, 2)) + (1.0*\\lambda + 2.0*\\mu)*Derivative(U_z(x, y, z, t), (z, 2))]])"
      ]
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "source": [
    "import numpy as np\r\n",
    "import matplotlib.pyplot as plt\r\n",
    "\r\n",
    "fig = plt.figure()\r\n",
    "ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])\r\n",
    "\r\n",
    "f0 = 10000\r\n",
    "t0 = 1/f0\r\n",
    "\r\n",
    "tRange = np.linspace(0,t0*4, 1000)\r\n",
    "vRange = np.exp(-(2*(tRange-2*t0)/t0)**2)*np.sin(2*np.pi*f0*tRange)\r\n",
    "envelope = np.exp(-(2*(tRange-2*t0)/t0)**2)\r\n",
    "\r\n",
    "ax.plot(tRange, vRange)\r\n",
    "\r\n",
    "ax.set_xticks([])\r\n",
    "ax.set_yticks([])\r\n",
    "\r\n",
    "ax.set_xlabel(\"Time, 0 to 4t0\")\r\n",
    "ax.set_ylabel(\"Displacement\")\r\n"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Displacement')"
      ]
     },
     "metadata": {},
     "execution_count": 25
    },
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAECCAYAAAAB2kexAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkEElEQVR4nO3deXRkZ3km8OetRUtJJZXW0m5JbUnd6sVqL70YrywGGzBOxxlowuZkIEDAeCDmhHBmcAYTmHgSTCeBYGwMxmCY8cSxjTcwOG5v3U43vS/qRepFW2lXlZYq1fLNH1WlLstqqSTVrXvr6vmdoyNVqere97ilx1fffb/vE6UUiIjIHCx6F0BERKnDUCciMhGGOhGRiTDUiYhMhKFORGQiNi0OWlpaqurr67U4NBERAdi7d++gUqps9vOahHp9fT327NmjxaGJiAiAiJyd63kOvxARmQhDnYjIRBjqREQmwlAnIjIRhjoRkYkw1ImITIShTkRkIgx1IiIT0WTyEdFKMTYZxDefOYoD50dx8/pKfOldTbBaRO+yaAXjlTrREoUjCn/x6B48ub8bhbl27PjdSXzz10f1LotWOIY60RI9daAbuzqG8a0/Wo/HP3c17nhHPX7y+hnsPTuid2m0gjHUiZYgHFH43osn0VpZgNsvrwEA/NVNLSjOy8L3fndS5+poJWOoEy3B66cHcWZoEp+7YRUssTH0vGwbPn1tI3aeGEB7n0/nCmmlYqgTLcG/7+uBM9uG97S63/L8f7myBjaL4P/9oUunymilY6gTLVIgFMYLR/rwvnUVyLFb3/K9kvxs3Li6HE/s60Y4onSqkFYyhjrRIu09O4LxQAg3ra2Y8/sfaqvCgC+AP5zjDVNKP4Y60SK9cnIQNotgS2PxnN+/rrkMNovgxWOeNFdGxFAnWrRXTw7i8roiOHPsc36/IMeOzY3FePEoQ53Sj6FOtAg+fxCHe8Zw9aUl877uxpZynB6YQM/oVJoqI4piqBMtwsGuMSgFbKwrmvd1Wxqjob+7cygdZRHNYKgTLcL+86MAgLYa17yva60sQGGuHW+cZqhTejHUiRZh37lRNJblodAx93h6nMUi2NRQjF0dw2mqjCiKoU6UJKUU9p8fQVutK6nXb2kswbnhSXRzXJ3SiKFOlKSukSkMjk8vOJ4et7kh2vLIBb4onRjqREk60uMFAKyvLkzq9S0VTmTbLDgQG4cnSgeGOlGS2vt8EAGa3flJvd5utWBddeHMzVWidGCoEyWp3eNFXbEDjqzkNwxrq3XhcPcYguGIhpURXcBQJ0rS8T4fWtzORb2nrdaFQCjCpXgpbRjqREnwB8M4MziB1RWLD3UAHIKhtGGoEyXhVP84IgpoqShY1PtqinJRkpfFm6WUNgx1oiQcjw2ftCzySl1EcFmtCwe6RjWoiujtGOpESTjVPw67VVBf4lj0e9dWFeD0wAT8wbAGlRG9FUOdKAmdg+OoK3bAZl38r8yaygKEIwonPeMaVEb0Vgx1oiR0Dk6goTS5/vTZ1lRGx+GP9XpTWRLRnBjqRAuIRBTODE2isSxvSe+/pNgBR5YVRxnqlAYMdaIF9IxNYToUQX3J0kLdYhGsrnAy1CktGOpEC+gcnAAANJQuLdQBoLWqAMd6vVBKpaosojkx1IkWcCYW6ksdfgGi4+o+fwhdI1yGl7TFUCdaQMfgBBxZVpQ7s5d8jNbYzVIOwZDWGOpEC+gcnEB9SR5EZMnHaKlwQoQdMKQ9hjrRAjoHJ9CwjKEXAHBk2dBQmoejPQx10hZDnWgewXAEXSNTaFhi50uiFrcTJ/s5AYm0xVAnmkfvqB/hiELdEpYHmK3J7cTZIS4XQNpiqBPNo2tkEkB0tcXlanbnI6Ki68gQaYWhTjSPeAtibdHyr9TjG2yc7OeGGaQdhjrRPLpGJmERoKIwZ9nHqi/Ng90qOMGFvUhDDHWieXSNTKGyMBf2JazOOJvdakFDaR5OenilTtphqBPNo2tkCtUpGE+Pa3I7eaVOmmKoE82ja2QyJTdJ45rLnTg3PInJ6VDKjkmUiKFOdBHToQj6vH7UpOAmaVyzO7omOztgSCsMdaKL6BvzI6JS084Y1xzb45RDMKQVhjrRRaSyRz3ukmIHsqwW3iwlzTDUiS4ilT3qcTarBY1leTjBUCeNMNSJLiKVPeqJmtkBQxpiqBNdRNdo6nrUEzW789E9OoXxADtgKPUY6kQXkeoe9bim+HIBHIIhDTDUiS6ie2QqpTdJ45pnQp1DMJR6DHWiOYQjCn1ePypTPJ4OAHXFDmTbLFzYizTBUCeaw/DENMIRhYqC1Ie61SJYVZbPm6WkCYY60Rw8Xj8AoFyDUAeiN0s5pk5aYKgTzaHfFw11t0ah3uR2omfMD58/qMnxaeViqBPNweMNAADcBdmaHL+pnGvAkDYY6kRz8Hj9EAFK87UJdXbAkFYY6kRz8HgDKMnLTvnEo7jaWAcMlwugVGOoE82h3+tHuVObq3QgoQOGwy+UYgx1ojl4fH7NxtPjmt35OMUrdUoxhjrRHDzegGadL3HsgCEtMNSJZgmFIxgcD2jWox4X74A5ySEYSiGGOtEsg+PTUEq7dsa4eAfMKXbAUAox1Ilmic8mdTu1vVJnBwxpgaFONEu/Lz7xSNtQt1oEl5azA4ZSi6FONMvMlbrGwy9AdFyda8BQKjHUiWbp9/phEaBEo9mkiZrcTvSyA4ZSiKFONIvHG0CZMxtWi2h+rpnlAjgEQynCUCeaJTrxSNvx9Lhmd6ytkUMwlCIMdaJZPN6ApksEJKopinfA8EqdUoOhTjRLv9ev+cSjuHgHDIdfKFUY6kQJpkMRDE1Ma96jnqjZ7eTwC6UMQ50owcC4tptjzKXJnY/eMT+87IChFGCoEyW40KOeviv1pvLYcgEcgqEUYKgTJeif2XA6fVfq7IChVGKoEyW4sDdp+q7Ua4scyLGzA4ZSI6lQF5EvJfMcUabr9/lhswiKHVlpO6clvgYMr9QpBZK9Uv/kHM99KoV1EBlCvEfdkobZpImayp0cU6eUsM33TRHZDuCjABpE5KmEbzkBDGtZGJEePGnsUU/U5M7HE/u64fUHUZBjT/v5yTzmDXUArwPoBVAK4B8SnvcBOKhVUUR66fcGUF/qSPt5m2MdMCc947jikqK0n5/MY95QV0qdBXAWwNb0lEOkL4/Pj00NxWk/78wuSP0+hjotS7I3SreJyEkRGRMRr4j4RMSrdXFE6eQPhjE6GUzrxKO4mqJc5NqtaO/juDotz0LDL3F/D+CDSqljWhZDpKeB2I5HeoypWyyClgonjvaOpf3cZC7Jdr94GOhkdnrMJk3UWlWAoz1eKKV0OT+ZQ7KhvkdEfiUi22NDMdtEZJumlRGl2YWJR+kffgGA1soCeP0hdI9O6XJ+Modkh18KAEwCuCnhOQXg31JeEZFOZq7U07hCY6LWqgIAwJEeL2qK0t+BQ+aQVKgrpe7QuhAivXl8fmRZLXA59OkTX1NRAIsAR3u8eO/aCl1qoMyXVKiLSDOAHwBwK6XWicgGALcqpe7VtDoyJK8/iNdPDaFndAouhx3XNJWiXKer21Tq9wZQXpANkfTOJo3LzbKioTQPR3vZWEZLl+zwy48A3A3ghwCglDooIr8AwFBfQcYDIXzvxRP4+e5zmJwOzzxvswg+eXU97n5vC3LsVh0rXJ7+NO5NejGtVYX4w9kRXWugzJZsqDuUUm/OuoIJaVAPGdR/nhnGXb/cj56xKdzWVo3tm+rQVB7d3OFnu87goVc7cbTHi4c+dSUcWcn+WBmLxxuYWQZXL2urCvD0gR6MTk7DlcZFxcg8ku1+GRSRVYjeHIWI3I7o8gG0Ajyxrwsf/dEu2KyCxz+7Fd/9cBs2NRSjKC8LrVUF+Pa2Dbj/w23Y3TmEux8/mLEteR6vX/dhpNbK6M1SDsHQUiV7SfWXAB4AsFpEugF0AviYZlWRYTyw8zT+7tnj2NJYjB9+7EoUXuQm4m0bq9Hn9eM7zx3HtZeW4iOb6tJc6fJMTofg84fSujnGXNbEQ73Hi6tXlepaC2WmpK7UlVIdSql3AygDsFopdY1S6oymlZHuHnylA3/37HF8YEMlHvmzzRcN9Li/uK4RmxuK8e3njmMottdnpuiP96jrfKVe5sxGuTMbR3t4pU5Lk+zaLy4RuRPANwF8S0R2iMgObUsjPT30aifufeYY3r++Evd/uA1ZtoV/VEQE9962Dj5/EA/s7EhDlamj92zSROuqC3Gom8sF0NIkO6b+LIB6AIcA7E34IBP66etn8M1fH8XN6ypw/0faYLMmv+thk9uJWy+rwiNvnMVgBl2te3z6ziZNdFmNC6cGxuHzB/UuhTJQsr+tOUqpLyulHlZK/TT+oWllpIvH3jyHbzx1BDe1urFj+0bYFxHocV94ZxP8oTB+8tqZ1BeokQsbTut/pd5W54JSwMEuXq3T4iX7G/szEfm0iFSKSHH8Q9PKKO2e2NeFv3niEG5oKcM/fXRpgQ4Al5bn48aWcvxqz3kEw5EUV6kNj9ePHLsFBTn6t2O21bgAAPvOsV+dFi/Z39ppAPcBeAMXhl72aFUUpd+zh3rxlf9zAFsbS/CvH7sC2bblTSL60811GPAF8OJRT4oq1JbHG4C7IEe32aSJCh12NJblYf/5Ub1LoQyUbKh/BcClSql6pVRD7KNRy8IofV440oc7H9uHy+uK8OAnr0zJrNAbWspRVZiDx/7zfAoq1J7H69e98yVRW60L+8+PZmzPP+kn2VA/hegqjWQyj+/twuce3Yv1NYV4+I6rUjYb1GoR3LaxGq+dGsTwxHRKjqmlfl9A9x71RBtrXRgcn0bXCJfhpcVJNtQnAOwXkR/G2xnZ0pjZlFJ46NVO/NX/PYCrV5Xi0T/fDGeKd7F//4ZKhCMKLxzpS+lxtdDv1X/dl0RttdF9SjkEQ4uV7GXZv8c+yASmQxF846nDeOzN83jf2gp8b3vbssfQ59JaWYD6EgeePdSL7QaeYToeCGFiOmyIdsa41ZVOZNss2H9+FB+8rErvciiDJLue+k9FJBdAnVKqXeOaSEPnhydx16/2Y+/ZEfzljavw5fe0wGrR5uagiOCW9ZX44c4OQy9QFZ94pPe6L4nsVgsuq3Fhz5lhvUuhDJPsjNIPAtgP4PnY4zYReUrDuijFlFL4xe5zeN/9O9He58M/f3Qj7n7vas0CPe5da9wIRxReOTmo6XmWYybUDXSlDgBbGotxqHuMk5BoUZIdU78HwCYAowCglNoPgN0vGWJXxxBu+5fX8DdPHEJbnQvP33UtPrAhPX/St9W64HLY8VJ7f1rOtxQz674YaEwdALY0liCigD1n2K9OyUt2TD2olBqb1cObGbNKVih/MIwXjvTh4dfOYP/5UVQU5OC+2zfg9itq0tqLbbUIrm8uw8vtA4hEFCwa/2WwFEZa9yXRxroi2K2CXZ1DuHF1ud7lUIZINtSPiMhHAVhFpAnAnQBe16Kg88OTmI7NQhRgJoDiURDPI4k9MzufRJJ/z8xbE74/+3sLHWvmk8z/HosI7FaLZsMdgVAYx3t9ONg1ildPDWLniUFMBcNoKM3DNz7Yiu2b6nTbleiGljI8ub8Hh3vGsCE2W9JIPN4A8rKsyM/WfzZpotwsK9pqXdjVYcxx9UhEoXt0Csd6vTg3PIkBXwADvgAmpkOYDkUQDCuEIhHYrRbYLNGff7vNArtFkGWzRD+sVmTbLciyRh9nxz7e8v3ExzYLBNHfp+jvXOzzHF9bJP77Km/5/TQKu9WC2uLUbzCe7E/xFwF8HUAAwC8AvACNtrL75MNvomNgQotDG4LdKsiyWpBtt8Y+X/ghzrZZZ3647VYLsmyxX4TYR5ZVICIIhMKYmg5jKhjGgC+A3jE/+n0BhCPRiSoVBTn44yuq8b61lbh6VYnuV8fXXFoGAHj99JAxQ90A29hdzOaGEvzg5dPw+YMpbzldiqHxAJ4/0oedJwbw+ukh+PwXNkDLslpQmp8FZ459JoCtIpgIhRAMKwTDEQTDEUyHIwiGFKbDEUyHIgiEwgiGV94kq1VlefjdV25I+XGT7X6ZRDTUv57yCmb5+i1rMB4Ixc4bO390w6ULj2een6nvwuOLvQez3zvH9xOP85bzXeT5Weef7/hhpRAMKQRCYQRCF36YL3wdexyMYDIYRjAUmfklCIYVArHHkYhCTpYVuXYrcuwWlOZnY+uqElS7crGmsgAbagpR7co1xHT3uDJnNlaV5WF3xxA+e/0qvct5m36v33A3SeO2rirBP790Crs7hvHuVrdudew5M4yHXzuD3xztQzCsUO3KxfvXV2JDjQurK51oLM1DYa59yT93kUgs5MMRBIKRmcCP/57MfB2OACr6+xWJfgml1IXPKv5c7DUKhp2Vq9VfhkkdVUR+C+BPlFKjscdFAH6plHpvqgt61xr9fnBJO1saS/Dk/h6EwpFFLeWbDh5vABvrXHqXMacr64vgyLLi9+39uoT64e4x3PdCO14+MYAihx2f2FqP26+oweoKZ0ovHCwWQY7FGh0iNOYfTRkj2f9VlMYDHQCUUiMiwjs3lLQtjSX4+e5zONLjxWW1Lr3LmaGUiq77YtDhl2ybFddcWoqXjvdDKZW2v8D8wTC++9sT+NErHSjIteNrN6/GJ7bWIzdLn/sylLxkQz0iInVKqXMAICKX4MLoA9GCNjdGV2re3TlkqFD3+kMIhCIodxpz+AUA3rm6HL856sHxPt/MHqZaau/z4XM/34uOgQls31SLv755DQpz9R/Pp+QkG+pfB/CqiLyM6E3kawF8RrOqyHTKnTloLMvDro5hfOY644yr9xu0nTFRvJ3x98f7NQ/1pw/04KuPH0R+jg2P/vlmXNPEza8zTbIbTz8P4HIAvwLwSwBXKKVe0LIwMp/NDcXYc2YYkYhx/sjzxCYeGflK3V2Qg3XVBfjdMW3Xpn/wlQ588bF9WFtVgGe+eA0DPUMt5o5VGEA/AC+AVhG5TpuSyKw21hbB6w+hY9A4LatGnXg0283rKvGHc6PoGtFmBewHX+nAvc8cw83rKvCLT28xxLZ+tDTJrv3yXwHsRLQ//W9jn+/Rriwyo3iHiZGWk/X4jLnuy2y3xlZqfPpAb8qPnRjoO7ZvRJbNWN1JtDjJ/ut9CcBVAM4qpW4EsBGxdWCIkrWqLB/ObJuh9t7s9wbgzLGlbHMQrdQWO3B5nQtP7u9O6XFnB/pS96Ul40j2X9CvlPIDgIhkK6WOA2jRriwyI4tF0Fbnwr5zo3qXMsPI7Yyz3baxGsf7oktBpAID3ZyS/VfsEhEXohtl/FZEngRwVquiyLzaal1o9/gwOR1a+MVpEA11Yw+9xN22sRqOLCseeWP5v3oMdPNKtvvlj5RSo0qpewD8dwAPAbhNw7rIpDbWuRCOKBzqGtO7FADR7hcjbTg9n4IcO7ZdXo2nDvQsa99XBrq5zfuvKSLFsz8AHALwKoD8tFRIpmKkvTeVUuj3+TOq0+OTW+sRDEfw4CsdS3o/A938Fro7tBfRmaNzzU1W4EYZtEjFeVmoduXiSI9X71IwMhlEMKwyZvgFAJrcTtx6WRV+/FonPnV1/aL+h/SjnR341rMMdLOb919VKdWglGqMfZ79wUCnJWmtKsDhHv2HXzKlR322L7+nGaGwwneeO57U65VSuP/FE/jWs8dwy3oGutkl/S8rIttE5B9F5B9E5DYNayKTW1dViM7BCUwE9L1Z2u+Lb2OXOVfqAHBJSR4+f8Mq/Nu+bjx3aP6+9XBE4d5njuH+F0/ijy+vwY6PMNDNLtnJR98H8FlEx9MPA/isiPyLloWRea2tKoBSwLFefYdgZjaczpAbpYm++K4mbKgpxN2PH7xoi+PQeACfevhNPPRqdKjmvts3GG7ZY0q9ZP+F3wngvUqph5VSDwO4JfYc0aKtqy4EEF2rW0/xxbzKDLzuy8XYrRY88PEr4XLY8ZEHduHnu89iOhTdBtLrD+KhVzvxrn98Gbs7hvGdbetxz61rdd8Bi9Ij2Wl0pwDU4UJvem3sOaJFcxdkozQ/S/ebpR5vAC6HXbe9W5erojAHj3/2atz5y334+hOHce+vj6E4Lwt9Xj/CEYUtjcX4nx9ah2a3U+9SKY2SDXUngGMi8mbs8VUA9ojIUwCglLpVi+LInEQErVWFOKx7qPszpkf9YioKc/Crz2zBzpOD+I/2foxNBlFdlIubWiuwvqZQ7/JIB8mG+v/QtApacdZVFeCBnR0IhMLItulzpezxBQy/kFcyRATXN5fh+uYyvUshA0g21PcAmFJKRUSkGcBqAM8ppYLalUZmtraqEKGIwom+cd2uKPu9fjSVc81wMpdkb5TuBJAjItUAfgPg4wB+olVRZH6tVdEdfI726nOzNBJR6PcFMq6dkWghyYa6KKUmAWwD8H2l1J8AWKddWWR2dcUO5NgtON7n0+X8QxPTCEdUxk08IlpI0qEuIlsB/CmAZxb5XqK3sVoEzW4n2nUK9UzuUSeaT7LBfBeArwF4Qil1REQaAbykWVW0IrToGOr9vvgSARx+IXNJdundl5VStyql/lfscYdS6k5tSyOza6lwYmhiGgOx6frpNLPhNIdfyGTm7X4RkfuVUneJyNOIrsr4FuxPp+VYXRG9Wdre50v7rM4Lwy+8UidzWail8Wexz/9b60Jo5WmpiM50bPf4cE1TelsLPV4/SvOzuLgVmc68oa6U2hv7/LKIlMW+HkhHYWR+Zc5slORlob0v/TNL+8YyZ29SosVY8DJFRO4RkUEA7QBOiMiAiHCGKaVES4U+N0v7vAFUMNTJhBbazu7LAN4B4CqlVLFSqgjAZgDvEJH/lo4CydxaKpw44RlHJPK2Wzaa8nj9cBcy1Ml8FrpS/ziA7UqpzvgTSqkOAB8D8AktC6OVYXWFE1PBMM4NT6btnIFQGMMT07xSJ1NaKNTtSqnB2U/GxtXt2pREK0lLrAMmnTNL+72ZueMRUTIWCvXpJX6PKCnN7nyIIK3j6pm6NylRMhZqabxMROZqTRAA/I2gZXNk2VBb5MAJT/pCvS8W6hUcUycTWqilMTO3hKGM0ux2pjfUx2Khzit1MiHOvCDdNbvz0Tk4MbPHptY8Xj+ybRYU5vK2EJkPQ51011LhRCii0Dk4kZbzebwBuAtyIMKNmMl8GOqku6byC8sFpEOf18+hFzIthjrprrEsD1aL4GSaQp0Tj8jMGOqkuxy7FfUljrS0NSql0DfmRwV71MmkGOpkCOnqgBmbCiIQirBHnUyLoU6G0Ox24uzwJPzBsKbn6ePEIzI5hjoZQrPbCaWAU/3jmp4nvuMRJx6RWTHUyRBaKvIBQPMhGA8nHpHJMdTJEC4pyYPdKpq3NcaHX8p5o5RMiqFOhmC3WrCqLB8nPdoOv/R5/Shy2JFt4woYZE4MdTKMZrf2uyD1jflRWZir6TmI9MRQJ8Nodueje3QK44GQZufoGZ1ClYuhTubFUCfDaHZHlwvQcmZp9+gUql28SUrmxVAnw4iHulYdMF5/ED5/iFfqZGoMdTKM2mIHcuwWnNDoZmnvaLTzhaFOZsZQJ8OwWgSXludrdqXeMzoFAKji8AuZGEOdDEXLDpjumVDnlTqZF0OdDKXF7US/L4DRydTva947NgWrRVDu5JU6mRdDnQzlws3S1I+r94xGN8ewWrjjEZkXQ50MpblCuw6YaDsjh17I3BjqZChVhTnIz7ZpEuo9o1Oo5E1SMjmGOhmKiKDJnfoOmHAkuuMRb5KS2THUyXCay6MdMEqplB1zwBdAKKIY6mR6DHUynOYKJ0YmgxgcT10HTLydkUsEkNkx1MlwWjRYA6Z3jD3qtDIw1Mlwmt3RXZBSuWFGfDYpl90ls2Ook+GUObPhcthT2qveNTIFZ7YNBTm2lB2TyIgY6mQ4IoLmcmdKO2DODk2irsQBEU48InNjqJMhNVdE2xpT1QFzfngSdcWOlByLyMgY6mRIzW4nfP7QzEbRyxGOKHSNTKGuhKFO5sdQJ0OKrwGTihUbPV4/psMRXqnTisBQJ0O6sLXd8m+WnhueBACGOq0IDHUypOK8LJTmZ6ekrZGhTisJQ50Mq6UiPyUTkM4NTcJqEU48ohWBoU6G1VTuxAnPOCKR5XXAnBueRJUrB3Yrf9zJ/PhTTobVUuHEVDA8s27LUp1jOyOtIAx1MqyZ5QKW2QHDHnVaSRjqZFhN8bbGZYyrjwdCGJqYRi1DnVYIhjoZVkGOHVWFOcu6Wdo5MAEAaCzNS1VZRIbGUCdDa3I70b6MXvWOweh7G8vyU1USkaEx1MnQWiqcOD0wjlA4sqT3n+4fh0WAS7hEAK0QDHUytKbyfEyHIjgbm0C0WKcHJlBX7EC2zZriyoiMiaFOhtZSsbxdkE4PjHPohVYUhjoZ2qXl8bbGxY+rhyMKnYMTWFXGm6S0cjDUydAcWTbUFTvQ7vEu+r09o1MIhCJYxSt1WkEY6mR466oLcLBrbNHvOz3AzhdaeRjqZHhttS50jUxhcDywqPfFZ6I2lTPUaeVgqJPhXVbjAgAcOD+6qPcd7fWisjAHRXlZqS+KyKAY6mR462sKYZHFh/qxXi/WVBZoUxSRQTHUyfAcWTY0u53Yt4hQ9wfDOD0wgTWVTu0KIzIghjplhLZaFw6cH4VSya2tfqp/HOGI4pU6rTgMdcoIbbUueP0hdA5OJPX6o73RFkiGOq00DHXKCJdfUgQA2HNmJKnXH+3xItduRX0JJx7RysJQp4zQVJ6Pkrws7OoYSur1+8+PYl11AawW0bgyImNhqFNGEBFsaSzBro6hBcfV/cEwjvSMzVzdE60kDHXKGFsai9Ez5sf54fn3LD3cPYZgWOHyOoY6rTwMdcoYW1eVAADe6Bic93V/OBcdd2eo00rEUKeMsaosH+6CbPxH+8C8r3vj9BAaSvNQ5sxOU2VExsFQp4whInhPqxsvnxiAPxie8zWBUBi7OoZxbVNpmqsjMgaGOmWUm1orMDkdxuun5x6C2Xt2BFPBMK5tKktzZUTGwFCnjLKlsQTObBteOOyZ8/u/P9YPu1WwpbE4zZURGQNDnTJKls2Cd7e68ezhXkxNv3UIJhxRePpgD65vLoczx65ThUT6YqhTxtm+qQ4+fwhPH+h5y/O7O4bg8QbwobYqnSoj0h9DnTLOVfVFaHbn44FXOhAKR2ae//FrZ+By2PHuNW4dqyPSF0OdMo6I4Cs3teBU/zgeeeMsAGBXxxBePObBHVc3IDfLqnOFRPqx6V0A0VLc1OrG9c1l+PZzx3BueBJPHehBXbEDn76uQe/SiHTFK3XKSCKCHds34ppLS/GT18+g3JmNn/7ZJjiyeJ1CKxt/AyhjFeba8fAdm+APhpFts0CEKzISMdQp4+XYOYZOFMfhFyIiE2GoExGZCEOdiMhEGOpERCbCUCciMhGGOhGRiTDUiYhMhKFORGQiopRK/UFFBgCcTfmBiYgo7hKl1Nu2+NIk1ImISB8cfiEiMhGGOhGRiTDUiYhMhKFOuhGREhHZH/voE5Hu2NfjIvJ9jc8tIrJDRE6JyEERuXyO17hE5PNLPL5VRPaJyK8TnrtLRBwJj68QkUOxGnYI1w6mFGCok26UUkNKqTalVBuAfwXw3djjfKXUksJ0EW4G0BT7+AyAH8zxGheApdbxJQDHZj13FwBHwuMfAPh0Qh3vW+K5iGYw1MlwROSG+BWuiNwjIj8VkVdE5KyIbBORv49d4T4vIvbY664QkZdFZK+IvCAilQuc5kMAHlFRuwC45njPdwCsiv31cF/s6v4+ETkcO/+HL1J/DYD3A3gw4bk7AVQBeElEXoqdq0AptUtFW9AeAXDbov9jEc3CUKdMsArAOwHcCuBRAC8ppdYDmALw/liw/xOA25VSVwD4MYBvLXDMagDnEx53xZ5L9NcATsf+ergbwDYAbQAuA/BuAPdd5H8e9wP4KoBI/Aml1A4APQBuVErdGDtX1wLnJ1o07nxEmeA5pVRQRA4BsAJ4Pvb8IQD1AFoArAPw29iwtBVArwZ1XAPgMaVUGIBHRF4GcBWAp+IvEJEPAOhXSu0VkRs0qIFoXgx1ygQBAFBKRUQkqC7MmIsg+jMsAI4opbYu4pjdAGoTHtfEnluudwC4VURuAZADoEBEHlVKfWyO89docH5a4Tj8QmbQDqBMRLYCgIjYRWRt7OsviMgX5njPUwA+ERsn3wJgTCk1++reB8CZ8PgVAB+OdbaUAbgOwJuJb1BKfU0pVaOUqgfwEQC/Twj0mePFzuUVkS2xrpdPAHhyqf8BiOJ4pU4ZTyk1LSK3A9ghIoWI/lzfD+AIgNUAXpvjbc8CuAXAKQCTAO6Y47hDIvKaiBwG8Byi4+RbARwAoAB8VSnVt4hSHwDwvIj0xMbVPw/gJwByY8d/bhHHIpoT134hU4t10WxTSk3rXQtROjDUiYhMhGPqREQmwlAnIjIRhjoRkYkw1ImITIShTkRkIgx1IiIT+f/64x3XKOZimAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {}
    }
   ],
   "metadata": {}
  }
 ],
 "metadata": {
  "orig_nbformat": 4,
  "language_info": {
   "name": "python",
   "version": "3.9.7",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.9.7 64-bit"
  },
  "interpreter": {
   "hash": "f373b9c250ed3b0a998e280c91e2688b2dd6d73ddffc7cc226348105802fd29b"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}