{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "source file for the documentation in the [jupyter notebook](SystemMatrix.ipynb) and converted to markdown by `jupyter nbconvert`\n",
    "\n",
    "Extracting System Matrices\n",
    "=======================\n",
    "\n",
    "We illustrate the system matrix reading in an example of a acoustic cavity coupled to an elastic plate with a piezo-patch.\n",
    "Them mesh is generated with the trelis using the [coupled.jou](coupled.jou), which produces the mesh file [coupled.cdb](coupled.cdb).\n",
    "It is intentionally coarse to better show the structure of the system matrix.\n",
    "We have an acoustic domain (yellow), plate structure (green), and piezo patch (purple).\n",
    "\n",
    "![mesh](mesh.png)\n",
    "\n",
    "The XML input [coupled.xml](coupled.xml) with matrial file [mat.xml](mat.xml) defines a coupled acousric-mechanic-piezoelectric analysis.\n",
    "The input is targeted at anaysing the system matrix structure, and not intended as a pysical model.\n",
    "\n",
    "## Defining Matrix Export\n",
    "\n",
    "Ths system matrix can be easily exported by cfs using\n",
    "```XML\n",
    "<linearSystems>\n",
    "  <system>\n",
    "    <solutionStrategy>\n",
    "      <standard>\n",
    "        <setup idbcHandling=\"elimination\"/>\n",
    "        <exportLinSys format=\"matrix-market\" system=\"true\" rhs=\"true\" baseName=\"coupled\"/>\n",
    "      </standard>\n",
    "    </solutionStrategy> \n",
    "  </system>\n",
    "</linearSystems>\n",
    "```\n",
    "in the sequence step.\n",
    "The result will be written in ASCII files in the working directory, called `baseName_step_<n>_sys_0_0.mtx`, where `<n>` is the frequency or time step.\n",
    "Depending on the analysis type you can also request the output of mass, stiffness and damping matrix (e.g. for eigenfrequency steps).\n",
    "The right hand side force vactor is exported in a similarly called file (`*.vec`).\n",
    "One can read in the matrices with python or matlab.\n",
    "\n",
    "\n",
    "## Analysing the System Matrix Structure\n",
    "We use targeted right-hand side load definitions in each PDE to *mark* certain DoFs, siuch that we can recognise them in the exported linear system by analysing the right hand side load vecor.\n",
    "This is done by using a harmic analysis with 4 frequency steps\n",
    "```XML\n",
    "\t\t\t<harmonic>\n",
    "\t\t\t\t<frequencyList>\n",
    " \t\t\t\t\t<freq value=\"1\" />\n",
    "\t\t\t\t\t<freq value=\"100\" />\n",
    "\t\t\t\t\t<freq value=\"200\" />\n",
    "\t\t\t\t\t<freq value=\"300\" />\n",
    "\t\t\t\t</frequencyList>\n",
    "\t\t\t</harmonic>\n",
    "```\n",
    "Then we define math parser conditionals in the `<bcsAndLoads>` section to switch certain load on for certain frequency steps:\n",
    "\n",
    "In the acoustic PDE we use a right hand side load density on the complete volume region\n",
    "```XML\n",
    "<rhsDensity name=\"V_air\" value=\"(f lt 50) ? on : ( (f lt 150) ? 0 : 0 )\"/>\n",
    "```\n",
    "which will be active in step 1 and produce calues in the right hand side load vector of the linear system.\n",
    "Adittionally we mark the interface DoFs in frequency step 4 by\n",
    "```XML\n",
    "<rhsDensity name=\"S_NC-air\" value=\"(f lt 250) ? 0 : on\"/>\n",
    "```\n",
    "\n",
    "In the mechanic PDE we use the same techniqe (with `<forceDensity>`) active in step 2\n",
    "```XML\n",
    "  <forceDensity name=\"V_plate\">\n",
    "    <comp dof=\"x\" value=\"(f lt 50) ? 0 : ( (f lt 150) ? on : 0 )\" phase=\"10\" />\n",
    "\t<comp dof=\"y\" value=\"(f lt 50) ? 0 : ( (f lt 150) ? on : 0 )\" phase=\"20\" />\n",
    "\t<comp dof=\"z\" value=\"(f lt 50) ? 0 : ( (f lt 150) ? on : 0 )\" phase=\"30\" />\n",
    "  </forceDensity>\n",
    "```\n",
    "where additionally the *phase* is used to mark x, y and z DoFs.\n",
    "Note that this must berepeated for all mechanic regions (`V_plate` and `V_piezo`).\n",
    "\n",
    "All electrostatic DoFs are marked in step 3 using\n",
    "```XML\n",
    "  <chargeDensity name=\"V_piezo\" value=\"(f lt 50) ? 0 : ( (f lt 150) ? 0 : on )\" phase=\"90\"/>  \n",
    "```\n",
    "for the volume regions and, in step 4 for the electode\n",
    "```XML\n",
    "  <charge name=\"S_piezo-top\" value=\"(f lt 250) ? 0 : on\"/>\n",
    "```  \n",
    "\n",
    "The value of the variable `on` was defined in the `<domain>` section\n",
    "```XML\n",
    "\t\t<variableList>\n",
    "\t\t\t<var name=\"on\" value=\"1e+3\"/>\n",
    "\t\t</variableList>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import and plot the matrixes\n",
    "We use python to read the exported data.\n",
    "For the matrix-market files use `scipy.io.mmread`,\n",
    "for the right hand side load vectors use `numpy.loadtxt`.\n",
    "Let's import some useful stuff first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.io import mminfo,mmread,mmwrite\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's import the system matrixes for the first 3 steps and print the size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(312, 312)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S1 = mmread(\"coupled_step_1_sys_0_0.mtx\").toarray()\n",
    "S2 = mmread(\"coupled_step_2_sys_0_0.mtx\").toarray()\n",
    "S3 = mmread(\"coupled_step_3_sys_0_0.mtx\").toarray()\n",
    "S1.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a harmoic anaysis we can't directly export mass, stiffness and damping matrix.\n",
    "However, having computed several frequenc steps, we can compute it from the frequenncy dependent system matrices.\n",
    "They are related by\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "S_{ij}(\\omega_1)\\\\\n",
    "S_{ij}(\\omega_2)\\\\\n",
    "S_{ij}(\\omega_3)\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "1 & j\\omega_1 & -\\omega_1^2\\\\\n",
    "1 & j\\omega_2 & -\\omega_2^2\\\\\n",
    "1 & j\\omega_3 & -\\omega_3^2\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "K_{ij}\\\\\n",
    "C_{ij}\\\\\n",
    "M_{ij}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Thus, we solve a linear system to obtain the matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = 2*np.pi*np.array([1,100,200]) # omega values\n",
    "A = np.vstack([np.ones(3),1j*w,-w**2]).T # coefficient matrix\n",
    "b = np.vstack([np.reshape(S1,(1,-1)),np.reshape(S2,(1,-1)),np.reshape(S3,(1,-1))]) # rhs vector\n",
    "kij,cij,mij = np.linalg.solve(A,b)\n",
    "\n",
    "K = kij.reshape(S1.shape) # stiffness\n",
    "C = cij.reshape(S1.shape) # damping\n",
    "M = mij.reshape(S1.shape) # mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's define a conceinance function for plotting the matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matplot(A,limfrac=1,cscale='lin'):\n",
    "    lim = np.max(np.abs(A))*limfrac\n",
    "    vmin = -lim\n",
    "    vmax = lim\n",
    "    fig,axs = plt.subplots(ncols=2)\n",
    "    if cscale == 'lin':\n",
    "        Ar = A.real\n",
    "        Ai = A.imag\n",
    "    elif cscale=='log':\n",
    "        Ar = np.log(np.abs(A.real))\n",
    "        Ai = np.log(np.abs(A.imag))\n",
    "        vmin = 0\n",
    "    else:\n",
    "        print(cscale,'not known')\n",
    "    caxR = axs[0].matshow(Ar,cmap=plt.cm.RdBu_r,vmin=vmin,vmax=vmax)\n",
    "    caxI = axs[1].matshow(Ai,cmap=plt.cm.RdBu_r,vmin=vmin,vmax=vmax)\n",
    "    fig.colorbar(caxR,ax=axs,extend='both',orientation='horizontal')\n",
    "    for ax in axs[:]:\n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "    return fig,ax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mass matrix looks like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD3CAYAAAAjdY4DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATXUlEQVR4nO3dW6wkx13H8d9/do2JL8RRdp372hFG9majmF0DApYHZJAwiIASFAEPCEGeECCMAAHyAxKPREKWwk2RQvISEW6JQEEQAkSKtQIEXifIq81NoCzGAcdYdiyvTOwzfx6me6anT3dP9XT3VHX39yOtTp+e6p7qOef8trqmqsbcXQCAw1vErgAAzBUBDACREMAAEAkBDACREMAAEAkBDACREMAAEAkBDACREMAAEMnJNoW/3k74rQ2HnDl/Ttceu7L1FQj1pWvX9PTTT9uhn/fUqVN+x5kzh35azETT73WrAD53/l5962PPrL9/+PpVSdKDN51dfX/pkh686ezWVyDUxYsXozzvHWfO6BK/qxhI0+81XRAAEEmrAL722BVJq5bvw9ev6sGbzq5bv7m8VZx/RbXi69bX9pBCnjNk/6HqC4wBLWAAiMTaLEd52m70H9Hrah/PW8XFr0Coixcv6tHLlw/+Jtx9Fy44fcAYStPvNS1gAIik1SiIM+fPNY5syPv36OcDgN1aBbC0Ha7lYWgAgHCtA7iqX5cgBoD2aAEDQCS9toABAOEYBQEAkRDAABAJAQwAkRDAABAJAQwAkRDAABAJAQwAkRDAABBJ9AAuz6Lrsgh52/0AEFP0AM7XDi5+32a7fK42+/vApz0A2Ff0AJaOh/CYFOs+5usAcHhJBLCkys+YaxtmsbogQlrnAFAWPYBDw5FbfQBTEz2Ay63H8i19XbmiqvJN+4fS9IYiAJRFD+AqbftS68ofuk+26g1FQhhAnegBXO5aqHtDa1cXROwQrqs3IQygTquPpY/x8d3Fj7fno+6njY+lxxSN+mPpGeYFYKqiB3DITLimEO4ycw4AYooewOUuhfJoh7oQriu/a/tQGBEBYJfoAbzLWLsgqv5jGUvdARxG9ACeYhdE3TXxBiKAougBvKsLom57VzjXdVkcQtM1AUAuegDvq20IA0BqRhPAbbomCGEAY5BkANfNaCtup9wFAQAhkgzgkFYrrV8AYxc9gLusBUEIAxiz6AFctRxl3WPF7TF0QTAbD0CT6AG8rzG0flOqC4D0jCaAx9QFUbccJSEMoCjJAB7rKIjy86XUHQIgPUkG8FhHQTQtxk7rF0BZ9ACeyiiIkE/EIIQBFCX/iRh9KX+yhkR3QGr4RAxM0ag/EaMv5ZYorVEAsc0mgKXq7oG+7BPmKS6lCeBwZhPA5X7YvkNunzBP9dM8ABzGbAI4hqm3Xrtc39RfGyDEbAK4qQ94qP7gqfczd7m+qb82QIjZBLDU/DltQ93q73PeMfUNd3nd6F7B3M0qgHNN4437fp59zKFvOJX/QICYZhnAdEHEx2sDzDSAc3XrS5R1CYqxt1SHxGuDuZt1AJfVTX8eS2utbd/xGK4JmDICuCT2mhL72HdtjJSvCZgDArhC12UtDx1qbd+0O9SbeYQ70IwALgkJ3pClMmMI6cPuUqYt+niBZgRwScgtemq37m3q26UMgH4RwAUh/aQphlRIl0JfZQD0hwAuqOt22OdNrBjh3FdXQ9fuiFT+YwJSRwCXtL1Fbyp/aH11NXTtjkjp7gBI2Ww+EaOL8qdpcHs+DD4RA1PEJ2LsaSpjZ7tM0BjbtQJjQgA3mMoEhi4TNMZ2rcCYEMA7hLwxV5ZiYPUxWSPF6wLGLIkATv0PO3Q88NBrC/ct5HUvlkl1JiAwVkkE8Bhuc6c4gSGkvoyCAIaTRABL1a2rvv6I+3pTqWkCQ1Prt21Ls+2xXYS0aneV6TIMj6DGnCURwOU/wr5u5ev6bvvU1Pptcx0hs9NSs+9aGU3HAnOSRACXw6uvW9hD3Ao3dUFM/Va8y1oZU39tgBBJBLB0vCXUdZRB137ZfVtxUw7hNmFLCAO7JRPAVepubdve0u9zmxtyWx3SJ1z1/SF06XeuE9pN0vbnRAhjrpIM4Lo/yJA/6NA33Nq+Mdf2uWMHSx9vru3S188JmKskA7hLeIUOCxti+FhTyzd2IA/hED8nYMqSDGCpWx9wmxCu2t5HyBtxYw+atn29bX9O1x670r2SwIgkG8DSfn2LdesZDK3qjbhYdRlKn33AVa/NmfPn+qgmMBpJB3BZ23fWY7Y4H75+NZm6HNqYfk5ATKMKYKn9O+ux/7i71mVMC/4U7fNzogsCczO6AA415Oy3tprGCe/aHuMMuTaK10EXBOZmsgGcsrpWcUotdwDDm10ApxJqhDCA2QVwCrfubYI3hfoWtZ3AMtY+bOAQZhfAKRhz67dtHbusFQFM3ewCeKg/+n3WVRhzCFdth5QP2Q/MxewCeKg/+jbnbdsF0fa2H8A4zC6AU9A2eMfQKgbQ3uwCuBxeMWfKEcLAvM0ugMtdBTH6IfsI3lj9p4yCAPozuwBOwZhbv4yCAPozuwAeqgsi9Dx9tX5jh3DVdkh5qX6aNWtBYG5mF8BDdUGEnqcuvMbSBdGHurqzFgTmZnYBnLIxtH4B9Gf2AZxamBHCwHzMPoBTupWfYhcEoyCAerMP4JRMsfXLKAig3uwD+NAhsO9H9Yw9hIsYBQGszD6AD30r33bYVsioiVzI7X4Kwc0oCGBl9gE8ZnkI72oZj7n1DEwZAVwhtRYjgGkigCu0nekV08PXrwYtX0krGEgPAVxS1/pNMbCKdRrbdGYABPAxYwqspj5gQhhIHwFcIWQUQiqauiBC15roU+p3DUBKCOCSvmZuHTp8Qlq2Q7Z+xzxbD4iFAC7pa+ZWjPCJGcJ0bQDtEcAV6sIz5Rbdrm6HtmX2kXp3DZAaArgkJIxSad3VjYIoit01AaAeAVzSNrBi6nMmHCEMHB4BXCHkVjqVW+w+JmIwRA2IgwAuCZmIkUow9TURg3HCQBwEcMmYgqnPiRipXyswRQRwhdDlIFPQ50SMXWVCgjjFOwUgVQRwSds1dVMZKdHXG21NZZr+02EiBtAeAVzSNpiayh/akCG86z8Uui2A9szdgwvfd+GCX7p0acDqoG8P3nR2qwuh6j+GkDKHcPHiRT16+bId+nn5vcaQmn6vaQE3SKV7oQvGAAPpIoAbpD4RIxQhDKSJAN4h9YkYuzARA0gXAdwg9YkYIRgDDKSLAG4wlWAihIE0EcA7pD4RI9ShJmIACEcAB5pyd0RomZCJGOVtAPUI4EBTuUVnIgaQDgK4hSl1R/RRBkA3BHAHoRM1UmwNhqzXUCxTN3St7lzl805hUgvQNwK4g9AxwqF9pyH7+7JPKzjkjbyyujIs2AMQwNGNbbZd277esV0fcEgEcALGMttu3zchCWGg2mwDeIg//H36fsfUN0oXBNCv2QbwEK2vqnGzU2sd0voF+jPbAJb6bX01TV4ICamh69cHuiCAfs06gPvUNAQrRpB26f6oK9M0DnrXTDhCGDiOAJ6oLi3TqjJ9zIQjhIFtBPAEdXlzrEo+ESOkXMiaw/vWA5gaAjhQyGiFptvwQ7bw+mr9SuH1brPS2j71AKaIAA7UdvWwpllkh0AXBJA+ArhGXVBU2dXSO/RtNl0QwDgQwDW6BMKQa0GEtBSHXO0sZLRDcTskkGn9Yq4I4MjGdrveZQxwyDYwJwRwAtpOxIh1ux4SnrtauU3b1x67MvxFAAkhgBEsZCJGl+0z588NfxFAQghg9KKvEAbmhABO3NhWS6vbpgsCOG62Abzr3fyu5wx9jl37Uxu21fZ1owsCqDfbAB7iFrhpJbS+RjvEvl1nFATQn9kGsNS8gllbdWESGsK76heyf2iMggD6NesA7lNVn2fIFOUxYRQE0C8CeCD59N253GYzCgJojwAeWNegYRQEMF3m7uGFzb4i6UvDVQczd4e7nz70k/J7jYHV/l63CmAAQH/oggCASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIiEAAaASAhgAIjk5NBPYGbnJD0k6W5JPvTzAcAAnpb0W5I+6e695Zj1eK7tE5uds9vufFxHX9PiNfdqcctrZYsT2/9OnJAkLbb2L9bbktbl8n2L4v7SP0lbZVflLStvMjMtFrbaXpgWVnhskT8mWba/XNYK51pk57PFplx23evvT2T/JK23TyxMJ4/tX2w9Lkkn7Hj5qnPV7rfN9Z2wVR3zffn3kjbfZ8+3dcx636o+6/KL/LhC+WyfmbRQ9tXy34XtfdnuQvnsddsqv9lnKp6ruG91vCSZu+RLmS8lX0r57/XWvuzfcpkdUyibl8uOyf9ZobyOPUe2f3kkXy6l5dFq++hoc8zySL5cfe/Z46v9y6zc0eb4rHzl+bJj8vP5cik/Wh3jR8v193601DIrX9zny6WWxfLlY5bFx462HltVq1jeN+daupZHXtjvm2OWvi7rR6vtzbm2j1vt96yenpVfba/q5Vq668hdRy4dZT/ffLv4tXn/Zt9S9Y+Vn+NZvaTLek4v6GVd0Cv1MT216COIe28B58Frt75ei1NntbjltX0/BQAc1G26Qffr1DqIX68bl2+313QO4t4CmOAFMHV9B3HnACZ4AcxNX0G8dwATvADmrmsQtw5gghcAtu0bxMEBbGZ32W13foHgBYBqVUFsZt/l7peqyreZiPE2f+FpLV51l+zm23uqLgBMzyt1UnfrZr2opSTdX1cuuAXs7h8xs9v9xWeeOvr8FZ04fU72qjfLjMl0ACBJLteTelGX9Zxu1kk9o5fucffP1ZVv1Qfs7l+RZGZ2miAGgJVy8D6p/7vH/cXa4M3tNQqCIAaA/YM312kcMEEMYI66Bm+ul5lwBDGAOegreHO9rgVBEAOYor6DNzfYamiSZGanF7e/9anl81+WveLVq1XGzFZhnG1n5bIVyBaSafN49lhe1rRYr1S22WeSVewvPi5l582+M2WPqfAc2VOajh2Tn3Zdp8Ixm6fMz7Z9vkXpmM2KYLZ1bF62WKf8KfJVv9artNn2+Wr3q/jcWu8rXMr6NTFtX8dm36YeKpVblSlex6aMCs+hrNyxx/PXuLJ8dm7b1EWlclY42NwlebbgaeF3Ot8vL6yQ5tnx5X35cZvzWHF/fu6tsspWRsu/utwL583LV+1fuqR8NbbiuVxeOG6rjvl51ser8P1q26vKuxee2uXL1TG+9KyKm2PcS3VevzzZvmWxvLbPXziXZ9d1/DnyembnXC8255tjXFvPr/wllmtZfLl8ta90qvVL6tl5j+/fHJOfa1l8rLhf0nN6SbfopL6gFxpHNbQ1aACvn8Ts1ZLuHPyJhvcOSR+NXYmeTOlaJK4ndWO/nqfc/T/7PulBAngqzOxf3f1bYtejD1O6FonrSd3UrqcvdMwCQCQEMABEQgC3877YFejRlK5F4npSN7Xr6QV9wAAQCS1gAIiEAAaASAjgFszsPWb2WTP7NzP7qJndFrtOXZjZu8zsipktzWy0Q4TM7AEz+5yZfdHMfi12fbowsz80s6fM7PHYdenKzN5kZp80s6vZ79kvxK5Tagjgdj4h6a3u/jZJn5f065Hr09Xjkt4p6VOxK7IvMzsh6Xclfb+kt0j6cTN7S9xadfJBSQ/ErkRPXpb0S+5+VtK3S/rZkf9sekcAt+Duf+vuL2ff/pOkN8asT1fufrXPaZWRfJukL7r7v7v71yR9WNIPR67T3tz9U5KeiV2PPrj7l939crb9vKSrkt4Qt1ZpIYD399OS/jp2JaA3SCpOEX1C/JEnx8zulHRe0j9HrkpSel0NbQrM7O8kVX3i6EPu/hdZmYe0ur360CHrto+Q6xk5q9jH2MqEmNktkv5c0oPu/tXY9UkJAVzi7t/b9LiZ/aSkH5T0PU0fN52KXdczAU9IelPh+zdKejJSXVBiZjdoFb4fcvePxK5PauiCaMHMHpD0q5J+yN2vx64PJEn/IumbzOzNZvZ1kn5M0l9GrhO0Whhc0vslXXX3345dnxQRwO38jqRbJX3CzD5tZn8Qu0JdmNk7zOwJSd8h6a/M7OOx69RW9qboz0n6uFZv8vyJu1+JW6v9mdkfSfpHSXeb2RNm9u7YdergoqSfkHR/9vfyaTP7gdiVSglTkQEgElrAABAJAQwAkRDAABAJAQwAkRDAACat7wWOzOxvzOxZM/tYaf/7zewz2WJdf5ZNQGlEAAOYug+q3wWO3qPV8LqyX3T3e7PFuq5pNTyyEQEMYNKqFjgys2/MWrKPmtkjZnZPi/P9vaTnK/Z/NTu3SXqFAqbEE8AA5uh9kn7e3e+T9MuSfq+Pk5rZByT9t6R7JL13V3nWggAwK1nf7HdK+tNVY1WSdGP22Dsl/WbFYf/l7t+369zu/lPZGtXvlfSjkj7QVJ4ABjA3C0nPuvs3lx/IFgzqtGiQux+Z2R9L+hXtCGC6IADMStZX+x9m9i5p1WdrZvd2OWd2jrvybUlvl/TZncexFgSAKcsWOPpuSack/Y+k35D0D5J+X9LrJN0g6cPuXtX1UHW+R7Tq471F0v9KerdWH1f2iKRv0GqN6s9I+pld6x8TwAAQCV0QABAJAQwAkRDAABAJAQwAkRDAABAJAQwAkRDAABDJ/wMa++EhUoKbdQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "matplot(M,limfrac=1e-12);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sftiffness matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADqCAYAAACcPZ9GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAg40lEQVR4nO2da4wkV3XH/3d2/RrbO+snrGGXXQLYgwEPC0pgJsmHVZQHEnLsKBKSpWAmEkKaSN6J8oHIX5D4FKRo1pFGRI40mEQofAkEEiXKa/OQdoAEllmwNRgDCzbxBntt02N27PUyffOh61Sfvn1vTfXU41Z3/X/SaKpv36q6t3vmX6dOnXOusdaCEEJI/UzFHgAhhLQVCjAhhESCAkwIIZGgABNCSCQowIQQEgkKMCGERIICTAghkaAAE0JIJCjAhBASif2jdDb7r7W37li8cN2NQ+8dfe0SDr5rFhvfeRpzdx1JfxOSl6effhoXL140dZ/31ltvtW86wr9VUg0/yvi7HkmA3/2ut+NxO4ePbpwGAKxsbwIAZuaX8GMAT5xZxRcOH8f9Z85gZn4JZ86sFhw6aRMLCwtRzvumI0dw5syZKOcmk0/W33XpLoiPHH4/AOAzz3yl7EMTQshEMZIFvPGdp7H/zrnU8l2enu29MXci7dNZ71m99z9ztqQhEkLIZDKSBXzP7VdXNQ5CCGkdI1nA5557Dftv6vl8AaSW72LiEwZ673XWV9PfhBBC/IzsA6ZvlxBCymEkC3juriO4/8wZ3B94f3l6Fpg7kf4mhBASZiQBBnouhs888xV85PD7UxeDdkksbpzGWvKbEEJImJEFOI1y8LQtT8+m4rs2dwIrpQyREEImkz1ZwEDPFyyhZrSACSFkdEZ+CNdZX0VnfXUgzlfaOuurWNneTH8DgOTf6Tw8vS3iLTHFqZg7hNpH7VM2y9OzA3OYmV/qx0erPoQQ4lJpMR4DwDq/kWzL+4sbp1PXhYSuaVET3FA3jYieG/b22KG5ykVZLjh6Wy4+ug8hhLiM7IIYBVd8tfBa9XplexNIRFi/DuEK7U3H7sEaMORzfvDCBs5VbH1qvzeAdA6d9VU8dmgODx1bYDw0IcRLpQLss4AtBq3htbkTWJtfAuZODCRwrGxvDj3Ug3qtCwHpaAz9QHBlezM9TlVW6Mr25oDw6+0HL2zgwUrOSgiZBCp1QfjE12DYAu6sr3of2oloastYi28Wuk+VLoCZ+aXUzXH0gUdTd4j2+y5Pz9biDiGEjBfGWrt7r4Tjx4/bImX7XDcE1OsQjx2aw7nOZQB9IQaAR86fwbnO5dTCBXpCKyLXWV/F0QcexUvnz9EFMCYsLCzg7NmztdcDfk/Bv2tCslhYWMA3An/XlT+E0799LE/Pppaja0kCwEPH+rU0Hznf/yfR7drC1X1++LmPFhg9IYRUS+UP4eR3KBLC9aECg8Kp44rloZr2D6+gZyXLtvvgLUY8slw8bjp2T2qBS1teFwohZPKpfE243e4nxX3ghp5J+xfvW/bup9sfvLBRaIxF0H5fbcU/cv7M0IWEEEI0lfqAZ+aXvOFZy8qK1YgPNyuqwW13zzdKO2kW9AGTSSSaD1gnKOjXYr3uNVKhSdakznzzRUBIOyGEuFQqwOKn1W4GAHjp/LmB1z5rGBhOuBBG9aFWbf3q6Iy1uRMDr4FmXTAIIc2h0odwK9ubWJtfwlqSZOGi21MRTlKSBdeK9uFbm07ateujCvRxy7pgEELaQaUC7PpzgX7hHtkeuD13qqhpSzLkG9b9XEtzIIuOEEIaRqUCrB+4acE8+sCjaYSAbtdifMqpB5HH0nRZ3Dg9dJyq8IWZ6YeQtIIJIS6Vh6FpH7Dg+oBdTm1v4mQJ9RtWSjpOHtxoD2nz+YQJIQSoQYB1YR3dBvTDzrTluzZ3IhVNvY/uE4oqcIVuZn4p+ICvLGT84s/eUm4VFqYnhGRRqQtCI1XPHjl/Bg8dW8Dixml88b7lTNfCgHg6D+eE5elZdBIL1y2A06nB8nXHbFVbeqGpfBSEkHGkcgtYEEvwXOdyemsurohRj0EIIZNApQKsExN8sbGyGoZOXnBdDeJDDa0158bg7rZNCCFNofIoiJTErzuTuCFC9Rtcl4REEuhC7GtKpLUg59muAx1yB9QTj0wIGT9qc0HoaIiHji14U3az9hXxddd+C1m6dRVkD+GKPt0nhBCXygU4VIZRlgtyVxEORTvIop15oxryRE2Ugb6Q6G1fejJFmBCiqXVFDADpQpUSBbHXoula6HwxuC689W8+rIZGJpFo1dB8FImCCOEmPxBCyDhQWxxw+mAqcUPkjY11Ld1UbJ36D26kRF7LmBBCYlG5AM/ML+GmY/cMFVZ/7NDcUCSEL1zNFeyVQN+t9VVYII2QEBHurK8yEYIQ0kgqd0EsbpwecjUsT8/ioWMLhR+OichK7Qgdsibi66Y6E0JIU6i+HGWy6GbIBbFbOco8LggLvwvCAIP7EEJIg6gtEUPX8tWJGFnRCSuebS3cqbXrJGucSpI1Tk7Pptsx8FneXBWZECLU9hBO89CxhXTp+CJi5D5oE3E7OT2LreS4J2sqyrMbtMQJIS61CbCvMHve/YBBgXUTNLQIr80vYSsRYaBfUa1Oq1Os80U1xpXtTZwCcGB+iQ8FCSEAahBgiXYQAdSrYbhIZIMWWl/Yml7KSEc7GCAVXsGgnlRkX72HmfklbK2vpsK7uHF6aN06Qkh7qVyAxd0gIvjS+XPBpXq0D3cU14Eu2qNX0xBBjuVz3VpfxQHnQsFEEUKIULkAu/G/rgvCrRy2NncirZymaz5kuSBW1PZaYjXL+ap2QaRjmTuBrWQOYu2ubG+mwruiokHogiCEABFqQcRiWYm6fhgWs0aEvjDobb3qc93+65iwFgSZRLJqQUSJgoiBFt8mugG0Je/GQtNiJmQyqb0YTyy0+EqmXJlhYaH16rIInX9cwtWKZBkuT8+idlOXkIbRGgvYV6Tn1PZm+pAMKOaO8O27V9dBaDWPprkminxeK9ubyO/8ImQyaY0F7KNMAchrAeu2Iq6QJljJZVvAz3zzieKDImSMaI0FrDmVWF/L07NpskQV+CIwRvXpnlKWYsx17spmJbn70Fb0C9fdGHFEhNRPayzgzvpqGgZ3oIIKaaHbcddNoPvp90Lr3B1wQu5827HorK+W6raZu+tI0SERMla0RoB3owp/qrveHbD7gzkgbN1OkgXsc1/QBUHaBgUYGEr6qJKyLNdTExgbTBcEaRutScTIg1vPwVffoYmIn3lGZeC57wODF5rYERQ+mIhBJhEmYuRkKN1ZhayFxK0JyIM9PT4ttHIRWQEGtt00cEJIvdAF4SBW4jj6WH11MgghzYUWsCKtqqZeA4Orb+QJIavbstTjFNaS4kA2YxxVXWR0DWda14SEoQWs0MLhPpiTEDZ5ep/10K5oeNao+JI7Oon4AsNJIjI2qRYnfaRf0YVM5bgUX0KyoQArxAesaxWPgztCXA2hZY90AXu9TQiJCwU4gGsFymsR6LW5E40WMm2lu/PYUi4LnyWvL0LL07O1hegR0jboA1a4Sx3p23NZbcOq1/o9HR0RywccKrWpx7G4cXogtVnG7cvEc6Mm8hYA0mFxdEMQEoYWsIO2cDWy1NEox6gbPe6QO0JHSrhRE2WhY5IJIWGYiJGBzxoel+QMjev/9cUz6/ZYMc9MxCCTSFYiBi3gnIgPWG7zt8bsYZbPGs2K5vBFRxBCyoUCnIF2R2i3grgjxuUWO+SaCG0zUoKQeqAAZ6DjZAUpZyn4LMSmRQ3kicnVfdz5bAXaQ+wWJ00I6VGpAMs/oS8cCgCOPvBo2u775w79w4/6z13Umuusrw6sISe1hV3LOBSy1SRrUgttVtSEFuzdHj4WmV+TPhtC6qZSARbBCtUkeOn8OQDjUds2VmRDlYQiJUJ9tnKmNOc5bmhfQtpEpQIsVpROj9XojDP93m63sHmfzvss0r2m24bcEdIu0RGSxiznH+W2fZT2MhjVNeHGy+jU7VGPK98xCwaRNlOLC0IQl4Mg/3w+Mcwq5JLXBSGCqC03OWZZ9RpE5LU7QgvT4sbpXOcJ9WlSqJvvO/Klao/iA27S/Aipm1pdEOJycPHdwgIYWjWXNIu0QpxzZzCKa4PfMWkzlQpwSFiB7FtrvXimRsfi7oakwfqOkwd3fLut7+ZacnJ+SWyoOzIgj/sjjxsmq4+uEKeRu5e83zGjJkhbqbQWhK/OgOaR82fSfr46u7798j7gqWMBS9fyc2tHuOcfdUn6JtDxuG+y+gh5fdfyfa7NncA7sLGnMRIyrtSSivzYoTk8eGFjoO3oA4/ih5/7qLe/LoA+qo9Qpw/rouSCFIrRhWXKKJ7jS1F2z5E1p1Chm1B7EwrdaJEddSzL07NDxY0eve5n6G5fZCoymSiipyI/dGxhqO2l8+eCVpLOPivCVnKLq29vZTVhXXC9jBAzNyZYnyPP7fWodRfGPXzr1PZm6n6Q7/uWV16OPSxCaqWWKAhtEQo6CiHkA3RdDSJmeW5v9b56W1vDK0qM91p4RvtI5Ti++Wwl/uBRCe0TY3FQ93OXC06ozsRu35PsK/5gLktP2kYtPmDhnplr0m1dccu9fQ35hPNaxXLcAzVU9cq69ZaVlbfWV3FShaqNK77vyYd2A4U4kHw2nfXV1M1CC5i0jVprQbiuiL0E4ecRMQlvGiUjq0zcAj4W/XG3KfEgz1wNqqtLTEjTqXVFDHdlhiL7h5DjioW1uHHa+8BKHpqJBVbmCr5iwYuorCk/56jzDj1sa8JDuBCjjEvuUtCiCxMhAguyR0BHYgCD1rFbIhKI4++NAQuyk0mksiiIK13g0pVukUO0El07Ahj/iAZCyN4o5IK4agq4amo8Sgq7hWNiWpfLzgM511eq2ynOhEwurVkVOQ05izwOQQur9kMDgxeJpmXPaX95nrXjstagA1iMh7SbQubrt59/BV968sWyxtIatNCKO+LUGPl55eLRtqgOQsqmkAX8ztuuwztvu66ssVRKqPZsWRZYKNJC0m012vqVPienZ8ciEkA+w5Vk2031dlncOB38DNrycJGQEIyCiEBWCraOjhBhbkt0BKMgyCSSFQVRyAJ+dcdiCsDV+1jVdRRCD9Z0dEQnyZ4jhEwuhQT42jES3plAMoR+XZdF6VsKHgCgrF13lQ3p746ZEDK+tCYKQvyzbkRBkyIMgH5EhN72uSBiRUfo6IU8URC+cpoyP7f9mW8+UfJoCWk2hQT4a89ewhtvvBpvuPGqssbTCkKWrOuCEPFaa1CdBLfQ/G4XArf2sm5vcjo1IXVQSIB/6Y7ryxpH5fgyz6QmrVBEDEJREMCgpTczv4TFZNvngpD23aIjYrkgRo2CkM9di/Ta3AlgehZbnugIQtoEoyAahLuah+BGRwiMgigHRkGQKqksCuLSlS6u3T+FMXoW12hC69ixdgQhk0khAb7+qvGoA+HiWppNiavNEx0h/SjChIw/rYmC0GhfbdOiIPKgIyWajq/mg47s0Bc9RkGQtlFIgP/5Bz/F/OEDuGFMLeGmkWXh6hWVQyslN5GQpd7EQkOE1E0h5fz1Nx8cS/GVhT317yrwHTvP0jsiTr4ylW4x96bjm4e0bTEEjbQcRkGMMbrYz0wgnM5dINN3698UGAVBJpHSoiB2HK2+0gX2GWCKURBR0MV6QrhF33nrT0hzGEmA3XCzMfQ+ZKJDvdwi6U3M2NKrLQPwuibSW31VYQ1j4r4gZNJpZRRECFnNOPS6abj+5Jn5JXQ8S7wfUNl3YjHrlZqbeHEhpA1MmA1bDP3QbGZ+KX1d9cO6oqxsb6aWsNykLG6cHnjIJT5fsYhFfPkgjJB40AJWyMMpX0pwk32ny9Oz6GxvwqBn4S4CA3UkdNzwlqozLDWHm/hAjpA2QAFWiBCFEjWaJr7eaIfk9YHEHSFjlgvLmnJH+Arl7BW96rS+cNG9QUgYuiActMtB0O6IUY5RN7L+muuOcNHuiLIEUsf70q9MSD5oAXtwyyzqUK481qIu2Vgn4nbQ7gjXveC6IMrOqtOV2zp0bRCSCRMxxhi9zJLGV3fBxe0DDMcK123FMhGDTCJZiRh0QSjE9SDuB/e3bDctIkJHMriiKe4IaXejI6SPCC5D0wipDwpwBlvJumcuPqszBjKOk9Ozqcj6LgziEhAfsbvasoiyjpRo0gWGkEmFPmCFm3hhMWxRhhb3bAI6hlmjl0E6ML+ELWcppkUwWYOQGNACVmgXhOuKCCVoxETcC6dUVMMp9eBLux10dIQglm/60JDJGoTUCi3gDGTpdYmIkIdeQLxIB03qCkksVolwWISqC5GI8slkEUxgtGQN2WayBiHlQwFW+Go/dFQ4WtNW0tDjXZtfGvDtLqtQM5v0PeBJ1thykjWWp2e97ggR+ZArIk9EBiFkELogHHxRDqHXoRoR2n1RF2mqtHIj+KzWkDtC8LkjOombImsdOl1j2BdpQQgZhgKcg1CNCDftNiZ6uaI1x6WgWRF3xPqqNyICGIyC0MfJckPI56MjLSjChGTDRIwJRy+ACfQt3FGTNeqAiRhkEiltRYy24VvmR1vB4/Rg6pTyAYeSNaSPbOe/NBNC9gIFeERGrQsRG3FHuNER7sVD3AY6OmKcLjCEjCMU4Ax8t+AiSuMgvsKyJzrCdTsc2CVZg9EMhJQPH8JloCMiQgkaTUceErrJGiK+WzmTNQgh5UMLeETciIimo10QvmQNScTwlbKkO4KQaqEAZ5CVeDEut+Rp7YrtzVzJGr6VNbYSEc6TiKGPMS6fESGxoAsiB74ylfrHl4jRRHzuiLzJGqPGOtNtQcjuUIALoAvWjAPaBZGVrHEqkKyRZdG668D5EjwIIYOMlIhxz7uP26+u9wPW/+67L+Kpi5fwh/OHqxgbqRAdBRESVt3HjXuemV8Kxkb7XBC6v0a3MxGDTCKlJWLscw7xwbfdDLzt5r2PbIzwpeOK0OhqaUIe8YmJuBeyHq65Ap21flxaiwLjl6RCSCz4EC4nPkERYfJVS/P186HXY3NFu8pi6HmO6fZxkzj0hUVXZnMvUosbp9OHerodQGaBH0ImnUI+4C89+SI+/fVnyxpLKxkX/7GglzeSeskUUUL2RiEL+N472+F+qJLQEkdNW/pIu1cWPW2+Kmx0QxCSTeuiINxQsrLw1QbOs4KyzrTL0x6LzvrqYEyxWm1ZwtkG3BE5xLdpfnFC6oY+4JJY2d5MV40YxWrV667ZHO2x2VpfTWtEiDuikyG2OjxNfy5iKbPqGmkzrRNgsbrKsr70w7JF9NeRy3t8eXjlilCoPRYDLojE5zszv4SOrD2HwaLwekUNn49Y9mWaM2kzrRPgsvH5assW+VHJE+OrU4dDDEQqOGFry9OzQ8fWi3mmfTwREa67Qs7zjtojgAmJSyEf8Pd/+hp+dqVb1lhIjeQpKrS4cbpv3TsZf+62zpjzHVdeu1EfTfFxExKDQhbwLxy8uqxxkBLR2WtZfdZ2CSFzax/7rHzZHiry4wmvC51rnNK5CSmT1kVB5MWNXPBFMrhRDm7tYF/fOhCh2823KkV3QrgFiHyF3fW2Pp9bjEeiJVyknRERpI1QgHPiu63Wt+jy2icmdd9mZxXa0YRWRdbo44TEWsLOZuaXBtwRur8kb4REmK4I0kb4EC6AG07my/gaukVPIhd8YlL3bbZEGWQhtX5DkRZumrHu51rXy0kBdwCp+J70REO4n6F8rp3tTSwsLOScHSGTAZelj0xolWUJ5XIL4FRdIyLveIXdxuLWjtAP9VaS1TckrvjnT34Z3e2LrIZGJorKlqU/+3/bOHrwGtx87b4ih2k1IX+tiG/TVl92LdhQnK8g8wAAzJ3oJ2UkLgsJcSOkjRQS4OOvny5rHK0l60GZrjAmxKoR4S5VlNVH/MauZewr2yn9DYDpW75c8qgJaTb0AZNchFwMciHQrpSZZIl72YZaiWNle7PfJ3FDiJvilldermcyhDSEQlEQz7x8BTtNyZUlUXGXsndjkWXRT7GMB2pEJO6IF667sfZxExKTQhbw4RuvKmscxCFUqD3UHhvXXeJaxu4DOPmtC7b/DWgBk3bRujhgETBf2cgi+EpPLk/PwgAw6v1QgoYvZldu1zWxM8bcFaGlLYSUrJSECzcKYkuVs6QFTNpG6wRYkgNO7VIHYVTcpAw5x4H5pVwrBPuiIGbml4aSH3SZytjIw7YsS1zm4QtBk4uOJGLM3XWkrqET0ggYB1wQNy5X4lp1sRs35rVJroO8hNZw0wkXOgpC5qorpOl4YLmA6GSN6285zDhgMnFkxQFTgCvArZ2rn/RnrSw8zmiB1qnQuiym/gy0D1vaj0wfxPP2MgWYTBRZAtw6F0QdiJXoe9If24dbNSKsM/NLAxcaAwxZybItonz43XfHGDIh0WhtHLC7iGSZUQW6JoT7pH+U8enj7dYeG/ezc1+LW0ay3mTVEO2iOTJ9sK7hEtIIRnJBGGOeB/Cj6oZDWs6brLW31X1S/l2Tign+XY8kwIQQQsqDPmBCCIkEBZgQQiJBASaEkEhQgAkhJBIUYEIIiQQFmBBCIkEBJoSQSFCACSEkEhRgQgiJBAWYEEIiQQEmhJBIUIAJISQSFGBCCIkEBZgQQiJBASaEkEhQgAkhJBIUYEIIiQQFmBBCIkEBJoSQSFCACSEkEhRgQgiJBAWYEEIiQQEmhJBIUIAJISQSFGBCCIkEBZgQQiJBASaEkEhQgAkhJBIUYEIIiQQFmBBCIrG/6hMYY+4G8DCAOwHYqs9HCCEVcBHApwD8u7W2NB0zJR5r8MDG3G0OHn0cO69h6nX3YOqG18NM7Rv82bcPADA10D6VbgNI+0nblG53fgAM9O31N0l/A2MMpqZMb3vKYMqo96bkPcAk7W5fo441lRzPTPX7JfNOX+9LfgCk2/umDPYPtU8NvA8A+8xwf9+xgu2mP799pjdGaZPXAPqvk/MN7JO29caT9p+S/VT/pM0YYArJbyN/C4NtSbPqn3xuA/37bQb6WLqttz8AGGsB24WxXcB2Afm7HmhLfrrdZB/VV/ol+8iPUf0xdI6kvbsD2+0C3Z3e9s5Of5/uDmy399om7/fau0m/nf7+SX/v8ZJ95Hi224Xd6e1jd7rpa7vTRTfpr9tst4uu7u/u09Xv7Qy81xuW7m/7x+padHesarf9fbo27Wt3etv9Yw3u12u3yTht0r+33RuXRdda7FiLHQvsJN+vbOvf2e39ti7C77nn+Cmu4Cw6uISf4zhm8Pd4bqoMIS7dAhbhNTfegalbZzF1w+vLPgUhhNTKQVyFE7g1FeI7cE33g+Z1hYW4NAGm8BJCJp2yhbiwAFN4CSFtoywh3rMAU3gJIW2nqBCPLMAUXkIIGWSvQpxbgI0xbzEHjz5F4SWEED8+ITbG/LK19oyv/yiJGO+yly5i6qa3wFx/e0nDJYSQyWMG+3Enrser6ALAiVC/3BawtfYLxpjb7asvPrfz3Sew77a7YW46BmOYTEcIIQBgYfEsXsVZdHA99uNFXLnLWvtkqP9IPmBr7fMAjDHmNgoxIYT0cIX3WVy+y9pXg8Ir7CkKgkJMCCF7F16hUBwwhZgQ0kaKCq9QSiYchZgQ0gbKEl6h1FoQFGJCyCRStvAKlVVDAwBjzG1Tt7/jue7LF2Cuu6VXZcyYnhgn20m/pALZFGDQfz95T/oaTKWVyvptBjCedv0+kBw3eWWQvAd1juSUBkP7yGHTMal9+qeUow0eb8rZp18RzAzsK331mOQUUvUrrdJmBo8XbIc+N9I2NZX0MzEYnEe/rT8OOP16ffQ8+n2gzoGk39D78hl7+yfHNv2xwOln1M7GWgA2KXiq/qalHVZVSLPJ/m6b7Nc/jtHtcuyBvkgqo8lvC2vVcaW/r71rAUg1Nn0sC6v2GxijHCfdH+p1b9v6+lurTm1hu719bNcmQ+zvY60z5vTjSdq6uj8Gj6+OZZN5DZ9DxpkcMy02Z/v7WAycH/IRw6KrPy7ba3MOlX6kNjnucHt/HzlWV7+n2wF0cAU3YD+ewqXMqIZRqVSA05MYcwuAo5WfKMx9AL4Y8fwxaNucOd/JJ+acn7PWPlP2QWsR4NgYY75urX1v7HHUSdvmzPlOPpM4ZzpmCSEkEhRgQgiJRFsE+NHYA4hA2+bM+U4+EzfnVviACSGkibTFAiaEkMZBASaEkEiMtQAbY242xvyLMeap5PdNgX6/aYx50hjzPWPMx1X77xpjnjDGdI0x73X2+eOk/5PGmN+oei55KGG+3v2NMUeNMa8YYzaSnz+va04+QuNX7xtjzJ8l73/LGHN8t33zfnaxqGjOnzDG/K/6Xj9Q13x2o+B814wxzxljHnf2afR37EWyT8bxB8CnAHw82f44gD/x9NkH4PsA3gzgagDnALw9eW8WwJ0A/gPAe9U+b0/6XQPgWLL/vgmYr3d/9JJkHo89v93Gr/p8AMA/opcH9z4AX9vr3JvwU+GcPwHgj2LPr8z5Ju/9KoDj7t9sk7/j0M9YW8AA7gXw2WT7swB+29PnFwF8z1r7A2vtawA+n+wHa+2m9acV3gvg89bay9ba8wC+lxwnNoXmm3P/2GSNX7gXwF/aHl8FcNAYc2iXfZs896rm3FSKzBfW2v8C8KLnuE3+jr2MuwC/zlp7AQCS3761kt4AQKcQ/jhpy2Iv+9RB0flm7X/MGPNNY8x/GmN+pfyh5ybPZx/qs9e5x6aqOQPAHyS38GsNuiUvMt8smvwdeym1GloVGGP+FYBvBdCH8x7C07Zb7N1e9imFSPO9AOCItfYFY8x7APytMeZua+1WznOWSZ7xh/pE+94KUtWcPw3gk8nrTwL4UwCLexxjmRSZ70TReAG21v5a6D1jzE+MMYestReS25PnPN1+DOCwev1GAM/uctq97FMKFc/Xu7+19jKAy8n2N4wx3wfwNgBfLz6jkcnz2Yf6XJ2xb57PLhaVzNla+xNpNMb8BYC/L2/IhSgy3yya/B17GXcXxJcBfDjZ/jCAL3n6/A+AtxpjjhljrgbwoWS/3Y77IWPMNcaYYwDeCuC/SxpzEYrO17u/MeY2Y8y+ZPvN6M33B5XMYHfyfF9fBvB7yZPy9wHoJLecI8+9IVQyZ/GZJtwH4HE0gyLzzaLJ37Gf2E8Bi/wAuAXAvwF4Kvl9c9J+B4B/UP0+AOC76D15fVi134felfYygJ8A+Cf13sNJ/ycB/FbsuZY039D+vwPgCfSeRp8F8MHI8xwaP4CPAfhYsm0ArCbvfxuDESwjzb0pPxXN+a+Svt9CT5wOxZ5nSfP9a/TcZleS/9/fH4fv2PfDVGRCCInEuLsgCCFkbKEAE0JIJCjAhBASCQowIYREggJMCCGRoAATQkgkKMCEEBKJ/wdW97iBVXk/yAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "matplot(K,limfrac=1e-12);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The damping matrix is almost zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD3CAYAAAAjdY4DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXiUlEQVR4nO3dfawsd13H8c93du85t6e3D9felpaH24IILa0FERAElQdRCIrRBMUEJZpoMIqAolCIGv1HAaOSKCYEqn+AEAVRg8qDgaBUIEB5KlYoYmmlhdvS5156H3a//jEzuzNzZs7O7JnZ3+zO+5Wc3N3ZefjunnM/+9vf/uY35u4CAKxeFLoAABgqAhgAAiGAASAQAhgAAiGAASAQAhgAAiGAASAQAhgAAiGAASCQcZOVjxw54hcfPdpVLRi4r910k26//XZb9XGPHDniR5O/a5OUPTe0WIyXLAP2stffdaMAvvjoUV1zzTXtVAUUPPWpTw1y3IuPHtVHr7mmNGyleSgTvFjGXn/XdEEAmofroplRmDkFbSKAAe1u7RZZ4V8sVufNKvQbmlfcbrrtsghgDF7xPxJdDu3Yz+u3qmBu63e8bL2N+oCBTVb25Rsh3K2q1zfE6970mFZxuwlawIDKRzuU/Qu0iRYwBm+v1m6xT5iuCbSJFjCQcO1u6e7nSxpgEQIYSJh2t26zy2j5om0EMKC9+3qzjxHCaBN9wIAW9wMDXaAFDACBEMAAEAgBDACBEMAAEAgBDACBEMAAEAgBDACBEMAAEAgBDKxA7gw7ny5ep8ZyrD8CGOiInX5gfjv3QPl/uz7NjUvorwYBDHTExwc1uufW0GUsxSRNMilMIHeDAAY6NDn7Ikn1uiCqhAq/UabpzXwY3SCAgQ4U5xEeV7SE64RrVfjRKl1/BDDQgeL1wiZnXzRv+Wb6gJdpWWav4Iz1RgADq2KR7jqx/+4Hk2Qn7mulpGXroPXdDgIY6EK2n9ens/vnbkc6dnxSvknJsuI16Wbrbh/KjbLoTFJ32ZVCsH9MyA50ITvUrDDs7IKdUe5+eqWNpsPQfHxw2erqqxgyh3bw6gKBpN0IxSFfGA4CGFiBRd0II6sed1sVzivpgkCnCGBgBYrdCLMLfY4PzvpZq8bdjir6IFbSBVE85sqPuNkIYCCAXLeDRWtzxlzZewGhvDwCGOjAXqGUPjay+e3J2RfpwI2fLN9gcrJiR82GtO1XsSskvc+IiOURwEAH9gql4kkaqVOXPDG33izvRlvlgb7iEQrFrpCqrhHURwADPZMdHZG2fpkcZzMRwEBAVaMjZt0Lo61ZII9MOp352E8Irz8CGAioanSELJqFsG8fmj0+NvpeNwkBDKzQolZrttuhODoi++XdaZq/G4EABjpQnI4yVezLzY0sSIM30+0wOfsija774K5tx6bq0RFYGwQw0IGqkQ5SfvRAbiRBZrRDttthcsWzc2fMZUdHYL0RwECPVI128PFBnUgeMIXvgqhq4aMZAhjomeJohzSQt0em6PidkuIuCHvgnjAFKj8Kw6SVnxSyKQhgoIeyox2y3Q7TncOzvl8/eHaY4srmCGbayqXwqgErVOdKFmWjHUyatX412tLrPnpTJ/UtqikuJtL13zpR/hgaIYCBFUhbs759KN91kB3JkGlZpqE2tvny6c7h2bavetpRff62b3dbdIYpP/3lZedt587Sw3IIYGAFsqMdcl0H2ZEMVRfrzCzPbnvl+WestPW5a/pLRmHsGwEMrLHciRtYOwQwsO46bIl21cKm3zhGAAOo1FX/Lv3GMQIY2CAnMmdx0MrsPwIY2CDbI8udMUcIL2c/r1uTK1wTwMCG2c4MueCj/nL287o1uVIIAQxsqBNNmmIIggAGNkwau9sjk6ang9YyRE3e9ghgYMPkPgFHY/qBV6xJ9wUBDGy4TeoHrvNmsk5vOAQwMBCh5xBuKlduyTwZVeuv0+gPAhgYiOwUl1LzkFp1qFXNh1HVot/rKiRtavN1IICBARllZldrGlLr3pXRVnC2+ToQwMDQWLQW3RG1SqxxJY7clTt6hgAGBiQ3z/CCdUKrFZg1rsTRx+BNEcDAgOwKo5IW5F6BFSKc2+qr3m/tXTx3AhgYMotyV7rIKgucEK3JuqMasl0NZV0s+x0d0cVzJ4CBgfPxwcYt4VWrU0t2naoulj49J4kABiBJFunkml1Z3mvc7jsCGIAkaSvSWoVwcdxv1YkY2dt3Z55gH4KaAAYGqiyAtqJm89n2SdWJGNnb52wtPqFjlQhgoCN9z7GqABqZdHwdBgrvwzLPjlEQwBoxrUdrsqzEnbHluiOqRkqsq2VGRHQxxwQBDHSo9OoINc7eaqqL4VVbkXTfqbhWHx/cxxH6aZkuiLa7LQhgYNVqnL3VeJct7y9tuR86EEmTk7nH1qBRvzYIYKBDxbBqK7y6DsGRZbodRlsa3/7V2WMmddKK75NVvckQwECHii3TvVqqTfqLTcv1yzYJlmy3w+kjj9D42JczBWxOdDQ946/NSeE351UE1lT6n7XJ1XSl5fpli4doEsinL3jUrnHCq+6O6OJ4dV/2JrOq1d0nAQwE1rT/ts0RYnWOnT3cVqSg3Q+rGLtb9fIyFwSw5vaTndmpJK//1ok2yqllV59vZgKfTewP3u9wsyZnExLAwArtNTZ4UR9wdtvLztveNTqhU4W5Inx8UKM7vjZ7bNPsp7W7FdEHDPRWVV9vnT7g3DqjrVbqqSsNljRcJt9xsUZ33rzSGtYFfcDAGmn6zXqosbimfLfD5PDDFB2/M1A1/bfo90QAAz1Q+5v1JS+o2bpMt8N05/BanHIdwqL+ZAIYWCc97W9tOoRuSPZ6afr52wSwnjJfDNIoXowABtCe0VZ/uknWAAEMoF0W0fqtiQAG0JqqywKVrdNn+/lSscmmBDCwAfoSasVLAZVdWaOLic3btp8vFZtsSgADG6Cv/a07FdeH72u9q0YAA+jchl9ibmkEMIDOja1eCA8tpwlgAK0rC9Kx1ZtwaEgIYACtq5r1bWTSiYoU7tPpzG1MG1oHAQxsgB5l18zIyicQ2h6Z7jqxe9Lc4voh7aclzigIYGD6+tG9OCwtde52VBrOfX0eXSGAAQSR7aYYWvCmCGAAwYxqjo7YVAQwgCCy17i779RmXVeuLgIYwEqV9fceOhBJ09MhymkdoyAA9M7CiXqiceUQtXXCKAgAvVM1IiJre2RLtYTXNbYJYAD9Eo2lycn8l3Ne3ke87sPXCGBgKDIh1vsW42hLp6aZKiuuhbeuwZsigIGhyITYOgTXGSOrPD25928gNRHAwID1PchGptyFPlO73kAquij6jgAGhqIkpNbh6hQabS2eqMeitQxhAhgYijXuRx3tMZXlbHHF82tbm29YBDAwAG2FRsjW8siqrzG3Sm0ejwAGBqCtrobQreWdsZWG8Kq0fWQCGBiI0OG5jLLA2xlb4yF1bQVn268hAQxAUj+/jKsMPItyJ2Es6h/u65sPAQwMUcWIiL4o5mlZvprmU1mOTGs5yoMABoZoRSMGllV8M6gK0rEtHgVR3FefArnfvwUAK9OXYKpq/RaDtOnVNPrYHUEAAwOU7TNdOE3kitVt/RavpvGN4+WzqOWCt2cnaxDAwACNMilXZ5rI0KrqGmceuHBnXDqVZW7bJbteyi4g2gYCGBiAVQ7VCioqD+H9yraeOREDQCN1uhf60gWxb9FY0fE7299vZuhbWwhgYCDqtNz62gXR1HTncCf9vZyIAQCq0VovtFj72LongIEBqgqjPoZUlb3OgMuuYyePz273DQEMDFBVGPUxpMqkuTuqUbBv7fT2jYUABgam7x/L6yh+YbioRd/XNxYCGBiY4rjfbHjdfbJfJyrsZdf45RrzW9x076n5nQ6GqzVFAAMDlw2pc7bWLxKaXBHj6FkH5v3G0Ti/fQDr92oD6Nw6dU007V4YmZRt6IfsniCAAcz0bV6ItqUjIrYi6YY7TwZ/jgQwgJl1mBdiP3xrZ3b7kYe3mm/fZjEigAEUbMIY4b1kW/km6d5T9b94bDr15SIEMICcdR8jvEjxeZx1oHkMtjXsjQAGAEmanJzdvG9Bq3ivN6kmnxQIYACQpNGWTiRj1A4diJaezKfJJwUCGEAjm9IXXGY7N1N9dTxWvQaL5qYoIoCBDdNFQGb3Wfdj9qJLxffZtxck6V6XSaILAhiwLr4sK7tO2yJVE+X0+cu8NDzPGJmie7+557pt9AMTwAA2yn5a2NlQnZ71oKX3xSgIAIPUZgu77pWUlw1qAhjAxuikf7nGJD/Lhj4BDAxY07Pe+v4FWtv9y7nnmxkn3BYCGBiwpl8k9eELtDpvAnXfKOpcKXpmtNU4hBftnwAGBqRJC7YYtqFbv3VO8216KnDjN5TRVqM+4UX7J4CBAdlPCzZ067fO8duosc7VlheNcWYUBAAsoc443pHlT9bIBu+ikziyCGAAKKjTgj0jc6aJFZZzIgYArEDVl5V1QpgABoA91Jr7t8YVmcsQwACwh1oT1Fu01NhpAhjA4LUxxM6Uv9pydnkVAhjYMF5xu011vulfp7Psin22dWopC9utSLNJ3evshwAGNkxXVzbOhkmdb/qbXlsu1DjjsrG7e32Jli7fisrX2R6Z7OTxhfuRCGAANRXDJPSJGW3ZzxtF1Tq+tTM7bZkuCACt2JTQXVad7pPZ7dEWXRAAutdW/23o+SYWqTN5UdXt0u3c6z9lM7tN0tdqbwA0c7G7n7/qg/J3jY5V/l03CmAAQHvoggCAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAiEAAaAQAhgAAhk3PUBzOxySa+V9GhJ3vXxAKADt0t6vaQPu3trOWYt7iu/Y7PL7dxLrtPkpKIHPVbRoQtl0Sj/MxpJkqLc8mh2W9JsvXRZlF1e+JGUWzde35L1TWamKLL4dmSKLPNYlD4mWbK8uK5l9hUl+7Novl7yvGf3R8mPpNntUWQa71oe5R6XpJHtXr9sX5XLbf78RhbXmC5L70ua30+Ol9tmtiyuZ7Z+lG6XWT9ZZiZFSv619G8hvyxZnFk/ed1y68+XmbL7yi6Lt5ckc5d8KvOp5FMp/bvOLUt+ptNkm8y66XrJNumPZdbXrmMky6cT+XQqTSfx7clkvs10Ip/G9z15PF4+TdabzLdP1i/dX7JNuj+fTuWTeBufTGf3fTLVNFk/u8ynU02z6xe3mWYfm+Qei8vKru/zfU1d04lnlvt8m6nP1vVJfHu+r/x28XJP6vRk/fh2XJdr6q6JuyYuTZLfb3o7++/ey+fLpqp+rHiMu3RK1+pu3a/TerzO0Xt1LGojiFtvAafBa2c9WNGRyxQdurDtQwDASp2rA3qmjsyC+MHanv64PWjfQdxaABO8ADZd20G87wAmeAEMTVtBvHQAE7wAhm6/Qdw4gAleAMhbNohrB7CZPdLOveQGghcAypUFsZk9zd2vKVu/yYkYV/r9tys6/EjZmRe0VC4AbJ5zNNajdaYe0FSSnlm1Xu0WsLv/vZld4A/ccWzy5S9qdP7lssMPlxkn0wGAJLlct+gBXau7dabGukOnLnX3L1Wt36gP2N1vk2Rmdj5BDACxYvDeohOXuj9QGbyppUZBEMQAsHzwpvY1DpggBjBE+w3eVCtnwhHEAIagreBNtToXBEEMYBO1HbypzmZDkyQzOz+64Ipj03tvlZ1xXjzLmFkcxsntZL1kBrJIMs0fTx5L1zVFs5nK5stMspLl2celZL/JPVPymDLHSA5p2rVNuttZTZlt5odM95bfX1TYZj4jmOW2TdfN1pQeIp31azZLm+X3V7lc2WNrtizzVGaviSn/PObL5nWosF68TvZ5zNdR5hhK1tv1ePoal66f7NvmtaiwnmU2NndJnkx4mvmbTpfLMzOkebJ9cVm63Xw/ll2e7ju3rpKZ0dJ/Xe6Z/abrly2fuqR0Nrbsvlye2S5XY7qf2fbK3I9ve9n67plDu3wab+NTT0qcb+NeqHn28iTLptn1ld9/Zl+ePK/dx0jrTPY5m2zO59u4csdX+hLLNc2+XB4vK+xq9pJ6st/dy+fbpPuaZh/LLpd0t07pkMa6QffvOaqhqU4DeHYQs/MkXVJY/JOS3tP5wZuhpnqoqR5qqmcdajrm7je3fZCVBHDpgc0+5e5PCHLwCtRUDzXVQ031DLkmOmYBIBACGAACCRnAbw547CrUVA811UNN9Qy2pmB9wAAwdHRBAEAgBDAABBI0gM3scWb2cTP7rJl9ysyeFLKelJm91My+ZGZfNLPXh64nZWavNDM3syM9qOUNZvbfZvZ5M3uPmZ0bsJbnJL+vr5jZq0PVkannYWb2YTO7PvkbelnomiTJzEZm9hkze2/oWlJmdq6ZvSv5W7rezJ7Sg5pekfzerjOzd5jZwa6OFboF/HpJv+/uj5P0u8n9oMzsGZJ+QtKV7n65pD8OXJKk+D+1pGdLuil0LYkPSrrC3a+U9GVJV4UowsxGkv5C0nMlPUbSz5rZY0LUknFa0m+6+2WSnizpV3tQkyS9TNL1oYsoeKOk97n7pZIeq8D1mdlDJP26pCe4+xWSRpJe2NXxQgewSzo7uX2OpFsC1pL6FUl/5O4nJMndjwWuJ/Wnkn5buXNgw3H3D7j76eTuxyU9NFApT5L0FXf/qruflPROxW+gwbj7re5+bXL7XsWh8pCQNZnZQyU9T9JbQtaRZWZnS/pBSW+VJHc/6e53BS0qNpZ0hpmNJe2ow1wKHcAvl/QGM7tZcUszSCuq4FGSfsDMPmFmHzGzJ4YuyMyeL+nr7v650LVU+EVJ/xro2A+RlD1F9P8UOOyyzOwSSd8j6ROBS/kzxW/g0wXrrdIjJN0m6a+SrpG3mNmZIQty968rzqKbJN0q6W53/0BXx2t1NrQyZvZvksqu4PlaSc+S9Ap3f7eZ/bTid8IfDlzTWNJhxR8dnyjpb83sEXUuMd1hTa+R9CNdHr/MXjW5+z8m67xW8Ufut6+ytgwrWdaLTwlmdkjSuyW93N3vCVjHjymey+DTZvb0UHWUGEt6vKSXuvsnzOyNkl4t6XdCFWRmhxV/gnq4pLsk/Z2Zvcjd39bJAdMZiEL8SLpb87HIJumekPUkdbxP0tMz9/9H0vkB6/luScck3Zj8nFb87nxhD16rF0v6mKSdgDU8RdL7M/evknRVD16bA5LeL+k3elDLHyr+ZHCjpG9IOi7pbT2o60JJN2bu/4Ckfw5c0wskvTVz/+clvamr44XugrhF0g8lt58p6YaAtaT+QclVTM3sUZK2JN0eqhh3/4K7X+Dul7j7JYr/Iz3e3b8RqiYpHnkg6VWSnu/uxwOW8klJ32VmDzezLcVfmPxTwHpk8Zyeb5V0vbv/SchaJMndr3L3hyZ/Py+U9CF3f1HgspT8Dd9sZo9OFj1L0n8FLEmKGzdPNrOd5Pf4LHX4xWDnXRAL/JKkNyad3Q9I+uXA9UjS1ZKuNrPrJJ2U9GJP3gqR8+eStiV9MJlD+OPu/pJVF+Hup83s1xS3NkeSrnb3L666joKnSvo5SV8ws88my17j7v8SrqTeeqmktydvnl+V9Ashi/G4K+Rdkq5V/GnzM+rwtGRORQaAQEJ3QQDAYBHAABAIAQwAgRDAABAIAQxgkMzsajM7lox4amN/r0sm8LnOzH6mzjYEMICh+mtJz2ljR2b2PMVn9T1O0vdJ+q1kros9EcAABsnd/13SHdllZvadZvY+M/u0mf2HmV1ac3ePkfQRdz/t7vdL+pxqhDsBDABzb1Y8N8X3SnqlpDfV3O5zkp6bnEF3RNIzJD1s0Uahz4QDgF5IJk/6fsUT8KSLt5PHfkrSH5Rs9nV3/1F3/0Ayc+J/Kp7h7WOKz6Tb+5icCQdgqJLpQt/r7lckfbZfcveLWtjv3yie8GjP08/pggAASR5PGfq/ZvYCKZ5UycweW2fb5HJP5yW3r5R0paSF8wjTAgYwSGb2DklPl3RE0jcl/Z6kD0n6S0kXKZ5S9J3uXtb1UNzXQcUT+EjSPZJe4u6fXbgdAQwAYdAFAQCBEMAAEAgBDACBEMAAEAgBDACBEMAAEAgBDACB/D+brN6f+6mq0AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "matplot(C,limfrac=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determining the system structure\n",
    "\n",
    "The acoustic DoFs are marked in step 1.\n",
    "Read it and extraxt the indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
       "       17, 18, 19, 20, 21, 22, 23])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rhs1 = np.loadtxt(\"coupled_step_1_rhs.vec\")\n",
    "r1 = rhs1[:,0] + 1j*rhs1[:,1]\n",
    "Ia = np.argwhere(np.abs(r1)>1e-12).ravel()\n",
    "Ia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mechanic DoFs are marked in step 2.\n",
    "Read the rhight hand sinde load vector and also determine it's pahse to infer x, y, and z DoF.\n",
    "One can see that CFS assembles alternatingly dor each node (x, then y, then z)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 42,  45,  48,  51,  54,  57,  60,  63,  66,  69,  72,  75,  78,\n",
       "         81,  84,  87,  90,  93,  96,  99, 102, 105, 108, 111, 114, 117,\n",
       "        120, 123, 126, 129, 132, 135, 138, 141, 144, 147, 150, 153, 156,\n",
       "        159, 162, 165, 168, 171, 174, 177, 180, 183, 186, 189, 192, 195,\n",
       "        198, 201, 204, 207, 210, 213, 216, 219, 222, 225, 228, 231, 234,\n",
       "        237, 240, 243, 246, 249, 252, 255, 258, 261, 264, 267, 270, 273,\n",
       "        276, 279, 282, 285, 288, 291, 294, 297, 300, 303, 306, 309]),\n",
       " array([ 43,  46,  49,  52,  55,  58,  61,  64,  67,  70,  73,  76,  79,\n",
       "         82,  85,  88,  91,  94,  97, 100, 103, 106, 109, 112, 115, 118,\n",
       "        121, 124, 127, 130, 133, 136, 139, 142, 145, 148, 151, 154, 157,\n",
       "        160, 163, 166, 169, 172, 175, 178, 181, 184, 187, 190, 193, 196,\n",
       "        199, 202, 205, 208, 211, 214, 217, 220, 223, 226, 229, 232, 235,\n",
       "        238, 241, 244, 247, 250, 253, 256, 259, 262, 265, 268, 271, 274,\n",
       "        277, 280, 283, 286, 289, 292, 295, 298, 301, 304, 307, 310]),\n",
       " array([ 44,  47,  50,  53,  56,  59,  62,  65,  68,  71,  74,  77,  80,\n",
       "         83,  86,  89,  92,  95,  98, 101, 104, 107, 110, 113, 116, 119,\n",
       "        122, 125, 128, 131, 134, 137, 140, 143, 146, 149, 152, 155, 158,\n",
       "        161, 164, 167, 170, 173, 176, 179, 182, 185, 188, 191, 194, 197,\n",
       "        200, 203, 206, 209, 212, 215, 218, 221, 224, 227, 230, 233, 236,\n",
       "        239, 242, 245, 248, 251, 254, 257, 260, 263, 266, 269, 272, 275,\n",
       "        278, 281, 284, 287, 290, 293, 296, 299, 302, 305, 308, 311]),\n",
       " array([ 42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,\n",
       "         55,  56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,\n",
       "         68,  69,  70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,\n",
       "         81,  82,  83,  84,  85,  86,  87,  88,  89,  90,  91,  92,  93,\n",
       "         94,  95,  96,  97,  98,  99, 100, 101, 102, 103, 104, 105, 106,\n",
       "        107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119,\n",
       "        120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132,\n",
       "        133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145,\n",
       "        146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158,\n",
       "        159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171,\n",
       "        172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184,\n",
       "        185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197,\n",
       "        198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210,\n",
       "        211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223,\n",
       "        224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236,\n",
       "        237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249,\n",
       "        250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262,\n",
       "        263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275,\n",
       "        276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288,\n",
       "        289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301,\n",
       "        302, 303, 304, 305, 306, 307, 308, 309, 310, 311]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# mech is active in step 2, we use the phase to distinguish components\n",
    "rhs2 = np.loadtxt(\"coupled_step_2_rhs.vec\")\n",
    "r2 = rhs2[:,0] + 1j*rhs2[:,1]\n",
    "Im = np.argwhere(np.abs(r2)>1e-12).ravel()\n",
    "Ix = np.argwhere(np.abs(np.angle(r2)*180/np.pi - 10)<1e-3).ravel()\n",
    "Iy = np.argwhere(np.abs(np.angle(r2)*180/np.pi - 20)<1e-3).ravel()\n",
    "Iz = np.argwhere(np.abs(np.angle(r2)*180/np.pi - 30)<1e-3).ravel()\n",
    "Ix,Iy,Iz,Im"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In step 3 the electrostatic nodes get non-zero right hand side loads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,\n",
       "       41])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rhs3 = np.loadtxt(\"coupled_step_3_rhs.vec\")\n",
    "r3 = rhs3[:,0] + 1j*rhs3[:,1]\n",
    "Ip = np.argwhere(np.abs(r3)>1e-12).ravel()\n",
    "Ip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we look at the right hand side load vector from step 4 to get the interafce and electode indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs4 = np.loadtxt(\"coupled_step_4_rhs.vec\")\n",
    "r4 = rhs4[:,0] + 1j*rhs4[:,1]\n",
    "\n",
    "# the electrode\n",
    "I = np.argwhere(np.abs(r4[Ip])>1e-1).ravel()\n",
    "Ipe = Ip[I]\n",
    "\n",
    "# the acou-iterface nodes\n",
    "I = np.argwhere(np.abs(r4[Ia])>1e-6).ravel()\n",
    "Iai = Ia[I]\n",
    "\n",
    "# the mech interface dofs\n",
    "I = np.argwhere(np.abs(r4[Im])>1e-6).ravel()\n",
    "Imi = Im[I]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're rady to polt the system matrix and mark the DoFs in it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ff96a87bf10>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATMAAAD8CAYAAAAbkUOLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABaF0lEQVR4nO2deZwU1bX4v6d7FrZh30REwCAKsi8SzU9ccM2ikhh9SdTEROCXmDyzvKfRpzGLiS+/aIxmIWh4Sp4ENG4kMRoxohJkFRQRUQRmGJZhmxlmYWZ6ub8/uquprq6tu6una4b6fj7zma6qW7du3a4+de65554jSikCAgICOjqhYjcgICAgwAsCYRYQENApCIRZQEBApyAQZgEBAZ2CQJgFBAR0CgJhFhAQ0CkoujATkctEZJuIbBeR24tw/V0isllENonI+uS+viLysoh8mPzfp0DXXigiB0TkXd0+y2uLyPeT/bRNRC5th7bcIyJ7kn2zSUSuKHRbROQUEXlVRLaKyBYR+ffk/nbvF5u2FKNfuojIWhF5O9mWHyb3F+V58SVKqaL9AWHgI2AkUAa8DYxp5zbsAvob9v0cuD35+Xbgvwt07fOAycC7TtcGxiT7pxwYkey3cIHbcg/wPZOyBWsLcBIwOfm5Avggeb127xebthSjXwTokfxcCqwBZhTrefHjX7E1s+nAdqXUDqVUG7AEuLLIbYJEGx5Pfn4cuKoQF1FKvQ4ccXntK4ElSqlWpdROYDuJ/itkW6woWFuUUvuUUm8lPzcAW4GTKUK/2LTFikK2RSmlGpObpck/RZGeFz9SbGF2MrBbt12N/cNSCBTwDxHZICJzkvsGKaX2QeKBBga2Y3usrl2svrpFRN5JDkO1IUy7tEVEhgOTSGghRe0XQ1ugCP0iImER2QQcAF5WShW9X/xEsYWZmOxr7/VV5yqlJgOXA98QkfPa+fpuKUZf/Q44DZgI7APub6+2iEgP4GngVqXUUbuiRWhLUfpFKRVTSk0EhgLTReQsu2YXsi1+pNjCrBo4Rbc9FNjbng1QSu1N/j8APEtCFa8RkZMAkv8PtGOTrK7d7n2llKpJ/oDiwCMcH6YUtC0iUkpCeDyhlHomubso/WLWlmL1i4ZSqg5YAVyGj56XYlNsYbYOGCUiI0SkDLgOWNZeFxeR7iJSoX0GLgHeTbbhxmSxG4Hn26tNNtdeBlwnIuUiMgIYBawtZEO0H0mSq0n0TUHbIiIC/AHYqpR6QHeo3fvFqi1F6pcBItI7+bkrMAt4H5N+EZGFwPeAu71oi4i8KCJ1IvJXw/4RIrImOZO6NPkbLh7FnoEAriAxS/QRcGc7X3skiRmft4Et2vWBfsArwIfJ/30LdP0/kRimREi8Sb9qd23gzmQ/bQMub4e2/BHYDLxD4kdzUqHbAnyCxHDoHWBT8u+KYvSLTVuK0S/jgY3Ja74L3G31rHJ8ZrrGi7YAFwGfBv5q2P8kcF3y83zg/xbid+L2T5INCQgI6EQkJyz+qpQ6K7l9GvAbYADQDNyslHo/i/rOJ+GO8qnktgAHgcFKqaiIfBy4RylVNH+2kmJdOCAgoF1ZAMxTSn0oImcDvwUuzKO+fkCdUiqa3C76bGkgzAICOjnJ2dhzgKcSChWQcKZFRGYDPzI5bY+DluW72dJAmAUEdH5CJLSoicYDKjFD+0zGGc4cAnqLSElSOyv6bGmxZzMDAgIKjEr4xu0UkWsgYe8SkQl51qmAV4HPJXe196x/BsEEQEBAJ0NE/gScD/QnMaP5A+CfJJx9TyKxFGqJUspseGlW3xvAGUAP4DDwVaXUSyIyksQSxL4kZlq/pJRq9fZusqBQ06QkHPq2kVgTdruL8nOKOa3rx7b4pR1BW4K2dIS/ggwzRSRMYhr4chKr9/9NRMY4nDbH4Xh74pe2+KUdELTFiqAtPqFQNjO/RsMICAjwCBEJ+SlOWkFsZiLyOeAypdTXktvXA2crpW7RlZlD8k3SrytThvdOl6uxOIRtRG1Dm6Kiohd07U1lZSWnDkx8plu/vNt/8OBBBgwYkHc9naUdELTFio7Ulg0bNhxSSnnWWBFJCQ+llJmrRrtSKNcMRx8UpdQCEo58TB0SVuvn9MBwPFFR0i9GKYXORya5HQeOoFSP1Gc+9QOY+mXPbiQgoLMgIpUe1uU7T4hCNciTFft64WVFhma5cVG2lwkICMiemH5Dr6UVi0IJs7yjYRgFmd122rFYW7ZtDQgIyAI/amVQoGGmSiw8vQV4iUSc/4VKqS2WJ4RKMoaRJnW60tSIBsIsIKDAxMx2iogqpu2sYMuZlFIvAC+4KhyPklwqZlaP7akZAq7/KFeXDCgMkUiE6upqWlpaTI9XVlZy6qmnut52gx/qdHMNrYxVWafjbuvs0qULQ4cOpbS01FUd2eBXrQx8tDbTSutypY1phErg3H/3qEUBuVBdXU1FRQXDhw83/e6ampo488wzXW+7wVjeabsQdbq5hnZvVveo7cumD4x1KqU4fPgw1dXVjBgxwlUdWWKqlWkUUzvzxXIms9lMJ5raYrREod/AIWzbtYfRI4fBuM/DxT8sUCsD3LB161bOOOOM7F5CAZ6ilOL999/PEIgiskEpNTXXepNama0wS16/KF++b1VGJyHbrTRE364haKzh9H5hOLoX/vUgrH+sXdoXYE0gyIpLAfv/KpfXP7tQDbDDt8LMDaZfWuCaEVAEHnzwQZqbm1PbV1xxBXV1dcVrUGH4jstytxe0FRb4Vpjl/HYJXDMCioBRmL3wwgv07t27eA0qDF1dlitKxFnfCjMnLP3MAteMDseGylp+8+p2NlTWelbnVVddxZQpUxg7diwLFiwA4MUXX2Ty5MlMmDCBiy66CIAjR45w1VVXMX78eGbMmME777wDwD333MMvfvGLVH1nnXUWu3btoqmpiU9+8pNMmDCBs846i6VLl/LQQw+xd+9eLrjgAi644AIAhg8fzqFDhwBYtGgR48ePZ8KECVx//fWe3WMxEBHbv2Lim9lMMxuZfimT2X5TSoqb7SogOzZU1vLFR1fTFo1TVhLiia/NYMqpfZxPdGDhwoX07duXY8eOMW3aNK688kpuvvlmXn/9dUaMGMGRI0cA+MEPfsCkSZN47rnn+Oc//8kNN9zApk2bLOt98cUXGTJkCH/7298AqK+vp1evXjzwwAO8+uqr9O/fP638li1buPfee/nXv/5F//79U9ftqOS0Kqed8Idm1qW3o5R3Lf0DzaxDsXrHYdqiceIKItE4q3cc9qTehx56iAkTJjBjxgx2797NggULOO+881LuCn379gVg5cqVKW3pwgsv5PDhw9TX11vWO27cOJYvX85tt93GG2+8Qa9evWzb8c9//pPPfe5zKSGnXbejEgqFHP+K1raiXVlPj4FoTTEEmwOytJ917+9cJsA3zBjZj7KSEGGB0pIQM0bmH/VkxYoVLF++nDfffJO3336bSZMmMWHCBNPnyGpEUFJSQjweT+3TnIBPP/10NmzYwLhx4/j+97/Pj35kH6zV9cqVDoLTMLOY9+qbYWYkGqMkbH7MaZi5qbyMv/TozqFwiH4c4cxtT/H+kfdRKM7seyb1bfX0KutFfVs9Uwcl3GzW16xn6qCpTBw4sRC3E+CSKaf24YmvzWD1jsPMGNnPkyFmfX09ffr0oVu3brz//vusXr2a1tZWXnvtNXbu3JkaZvbt25fzzjuPJ554grvuuosVK1bQv39/evbsyfDhw/nrXxMJvN966y127twJwN69e+nbty9f+tKX6NGjB4899hgAFRUVNDQ0ZAwzL7roIq6++mq+/e1v069fv9R1Oyp+Fsz+EGZtjZSWWHeSXQduKi/jK4MHEk2VaUS9+UPT0EFKKcKSkJixeAwVVzzx6SeYOHAiImIadkjfBrvjTrg5R3+NYmFsp7HvrV4s+bR7yql9PBFiGpdddhnz589n/PjxjB49mhkzZjBgwAAWLFjA7NmzicfjDBw4kJdffpl77rmHr3zlK4wfP55u3brx+OOPA/DZz36WRYsWMXHiRKZNm8bpp58OwObNm/mP//gPQqEQpaWl/O53vwNgzpw5XH755Zx00km8+uqrqbaMHTuWO++8k5kzZxIOh5k0aVJKAHZEnIaRbp8DEVkIfAo4oJKJig3HvwjcltxsJJEt/W3bOov94wGYOv5Mtf6zNaCsnItDQNz0yI/79eHJih6QEess80doJuCmDJzCY5c/Bji/dbLtKzfDGidh0dHYunVrTsuHArzF7HvwYAXAhq5du062K6OUoqWlZZ1SarpDXeeREFKLLITZOcBWpVStiFxOIlu6rTOuPzSzsu5w04uw/AewfwuEQhAuT0SOPfv/wqAxCe/+Q9sTs5V11dDWBCIcCmfeQmYQR0FM4kUqpVizY42t4DPDSnuxE0SdzXYScGLi1TOslHpdRIbbHF+l21xNIiaiLf4QZgCnTIev/N36+HWLzff/899h9z8tT9OEmEKlPusFW49ePSwFjbbfTXiiXHF7jfbCSbC7HX4GdE5cPqMDRGS9bntBMrJ0rnwVsBEOCfwjzHKkX1f72S9NiJkNMUWEuIqn9qedl9zW/3cKGOlkczM7V5u59YMgA3Oh5HZfQOfHpc3soNMw0y0icgEJYfYJx7Z5ccFicmZfc/uMlQZhZHD3wTlf25i3L9tz7bb9htX0u1+m5QPah/Z0zRCR8cCjwJVKKUcHxA6vmb1/5P3MncomPhrpQ6YBXRPJatzMZurP66w/XDt7oN1QPODEoL2+axEZBjwDXK+U+sDNOR1emClMNBqb/k6znYmwr2lfYr9OcFkNq/TDU/1/0+uYZpMSS2HhF+w0RivtUb9/69athWlYgC/w0DXjT8D5QH8RqQZ+AJQm65gP3A30A36b/I1EnWZiO/ww8zOnfYYwCd8xpRQoZ18uPfWt9WnH7P6Mw0q7YabRZ8zOZqb/H+BP9AvH86UjhwfyapiplPo3pdRJSqlSpdRQpdQflFLzk4IMpdTXlFJ9lFITk3+OLiUdXjODpOCIJzUuwdQNwwylFEebj7p2zbAacln5tdnVGQzPTlxeeMFdagw/4udntsNrZutr1hNTsUQnZ9nPIkLXLl1NNSw7zcs4E6l9NqvfSEfVwNy+kSsrK1m/fr1DbQZ2r4U37k/894Bdu3Zxxhln8LWvfY2zzjqLL37xiyxfvpxzzz2XUaNGsXbtWpqamrjpppuYNm0akyZN4vnnnwcgFovxve99j3HjxjF+/HgefvjhVL0PP/wwkydPZty4cbz/fsJWu3btWs455xwmTZrEOeecw7Zt2wB47LHHmD17NpdddhmjRo3iP//zP1P1dOTwQH5eaN7hNbOpg6YSljBRFT1uh8K9YGuLxrj8wdcpKwnxyiM/IXKoku5nXUi4W2/m3HAdv7rnP/jWPf+PRxYt4Y5n3mH25KFs29/A0nVVzFm0nrkzT0tt97/qDo6ufYa2ve+n+Y6VDTmDvpd+g/5X3cGdz26mbMgZaW0oG3IGdzy7OXV89uShni7vyQU3fmRG/zgR4e9/d3QHSmf3Wnj8M4mgmuEyuHFZwucwT7Zv385TTz3FggULmDZtGosXL2blypUsW7aMn/70p4wZM4YLL7yQhQsXUldXx/Tp05k1axaLFi1i586dbNy4kZKSkrSQPf379+ett97it7/9Lb/4xS949NFHOeOMM3j99dcpKSlh+fLl3HHHHTz99NMAbNq0iY0bN1JeXs7o0aP55je/ySmnHM+N3RHDA/lZM+vwwkwjl9k2pSAeD/PevqMA9L/8m4QEYvHED/eJNVX0v/ybPLG6km6nf5wn1lTxp7VVqeMiwsvv1RBP/tC7jz6Hnmeew5Nzz0kJo/W7jnDt71cRiR2v85Qv38+GylqUUmyorOWUG3+Rdo2l63ezdM7HiyrQ3GiQZlpp1suZdr2REGQqlvi/6w1PhNmIESMYN24ckFgfedFFFyEijBs3jl27dlFdXc2yZctSARhbWlqoqqpi+fLlzJs3j5KSxE9Dvyh89uzZAEyZMoVnnnkGSCxqv/HGG/nwww8RESKRSKr8RRddlAoRNGbMGCorK9OEWUcMD+RnYda5hpkmOHZ+tHfaMCmuSNvW6tAf19erSN+OxeHpt6pT5z79VjXReHod0ZhKlXE63ukZ/n8SGpmEE/+H/x9Pqi0vP56HNRQKpbZDoRDRaEKLf/rpp9m0aRObNm2iqqoqlarN6pnR6giHw0SjUQDuuusuLrjgAt59913+8pe/pOUL1bdBf45GR7SbtqefWbZ0eM1s6qCplIZKaYsngjLaa2L6xeaJfW2153repkMNraafzco4HW9PFq+p4j8ffxUp785nzz6N6SP68b0FfyHUox9SUo6KtnHGqYO56dwRfOHsYSxeU5U4Xt6DMz42nJvOHcGknlle9JTpiaHlrjcSgswDrcwNl156KQ8//DAPP/wwIsLGjRuZNGkSl1xyCfPnz+f8889PDTPtNKb6+npOPjkR8j7baBgdMTyQV64ZhSAvYSYiu4AGErn0okqpqSLSF1gKDAd2AZ9XSnkX3N3AxIET+cOlf+AvH/2Fj+o+YuVbK/nE5MTKh9rWWkpDpUTiEUpDpWzavIlPTP4E7+09SkPbMSL104jWeZ8Vq665zfSzWRmn4+3F4jVV3PHsZsIVieVhz23ay3Ob9lLaX5dZu6SU7QcauePZzazdeTjtuLb/z9cNy/7ip0xvNyGmcdddd3Hrrbcyfvx4lFKp+GVf+9rX+OCDDxg/fjylpaXcfPPN3HLLLZb1/Od//ic33ngjDzzwABdeeGFWbeiI4YH8rEnmFQIoKcymKqUO6fb9HDiilLpPRG4H+iilbrOqA2Dq1Kkq6xmwzLZkhPnR79d/nnX/Cj480JjzF+PkZvGxAd1Z/t3zWbymiv96bnNqaKrnYwO6c9MnRvKDZe+m7GnG48u/e35O7cuFK3+9krer650LJulaGuJYJDMs0+OzhzBz+iQvmxaQA4UKATR06FDHEEB79uxxDAFUCAphM7sSeDz5+XFcJg7NF7sJAEh/o7y98mXLsm7WTFo5wGqMHNCDxWuq+P4z7xCLm6/bbGyJcsezm2mLxk2Pl4bb15w5qGeXrMr3KDdX6sM+fnMH5I+fXTPyvbIC/iEiG0RkTnLfIKXUPoDk/4F5XsNzjq59xnS/mQEzW+1NgLkzT2PpuirLOgGOJIeRVscjMfNglIVi7szTCBmaYXXn4ZDw7YtHY5S34ZBQ0aW0IO0L8AedeQLgXKXUXhEZCLwsIiarvs1JCr85AMOG5WBnyawPMI/kanQhKO1/qmsnVzOsFqSLCArYtr+Btmh+wqi9NbMpp/bhqXnnMP+1jzhwtIVrpw1j9OAKnn6rmu01DRxpaqNv9zJGDapI+cFpxw81tDKgopzZk4dS1ry/Xdsd0L742WaWlzBTSu1N/j8gIs8C04EaETlJKbVPRE4CDlicuwBYAAmbWT7t0GMl1LRjSikqxl+cVjaf65htL11X5ahZOd1we2tmkBBoj9wwlQ2Vtcx/7SMWrtxB3+5l9O6WLsQ2VNZyx7Ob04TcgIpy5wsEdHg65WymiHQHQkqphuTnS4AfAcuAG4H7kv+f96KhTpgZ/42fNT598Uz+8V5NztdyGooO6tnFcTayPBwiErPKeQB9uxcnmfGGylqu+d0qYiYa7pJ1Vfz4ynHc9dxmooZZjbW7almyropn/m14O7Y2oL3xs2aWz1hmELBSRN4G1gJ/U0q9SEKIXSwiHwIXJ7fbFacOf/z712e7jNP9tUnYn5woCfvvodhQWcvc/11PHPOgi7E4/PAvW4gp6+MNLRGbKwR0dPxsM8tZmCmldiilJiT/xiql7k3uP6yUukgpNSr5v90XnDmpuqtWrSKW9AK3CuHjVL9R+zOef6TJXjNrM3HHyOZ8r9E0soNHW23jlrU62AKNGptfyDWEz4oVK1i1apVzQQ/rMpabP38+ixYt8qQN+eJnYdbhVwBoGJcfmfmZaVx26y/odvrHbTveaamJ1VBTAff9fSv761tMzjpOa8R6iAntPwEw/7WPUhqZFW4eVP/pm/mxYsUKevTowTnnnJNxLBqNptZw5luXXbl58+Zl1+gC4mebWYdfm6lh1LL0bwmjcDvnossdf5h2S6LstgHW7jxMQ2vU9hpRC/8zjfaeADhw1F74uiYHabbpwCYe3fwomw5s8qQJ//u//8v06dOZOHEic+fOJWawTVodf/HFF5k8eTITJkzgoosuYteuXcyfP59f/vKXTJw4kTfeeIMvf/nLfOc73+GCCy7gtttuY9OmTcyYMYPx48dz9dVXU1ubWOzy0EMPMWbMGMaPH891111nWtdf/vIXzj77bCZNmsSsWbOoqakxLXfPPfekFsRv376dWbNmMWHCBCZPnsxHH33kSZ+5JdDMfMa104bxdvXmnM5154fm/KU6HW9vzSyfPtHTvSy7R2rTgU3c/I+baYu1URYu45FLHmHiwIk5X3/r1q0sXbqUf/3rX5SWlvL1r3+dJ554wvH45Zdfzs0338zrr7/OiBEjUusk582bR48ePfje974HwB/+8Ac++OADli9fTjgcTsU8mzlzJnfffTc//OEPefDBB7nvvvvYuXMn5eXl1NXV0bt374y6amtrWb16NSLCo48+ys9//nPuv//+jHKvvPJKqv1f/OIXuf3227n66qtpaWkhHm/fl56fJwA6jTAzLmfSMM5yigijB1egYlEIhTPq0MrZbbcH7a2ZfeHshK/fAy9vo645QlggphLLluJAayRO19IQg3p15awhPXl3Tz37jrYQiSr6di9lUM8uXDttGN3Lm7K67vqa9bTF2ogTJxKPsL5mfV7C7JVXXmHDhg1MmzYNgGPHjjFw4EDH46tXr+a8885jxIgRgH04nmuuuYZwOEx9fT11dXXMnDkTgBtvvJFrrrkGgPHjx/PFL36Rq666iquuusq0nurqaq699lr27dtHW1tb6tpWNDQ0sGfPHq6++moAunTJbtWGF/h5mNlphJkVRlsaJEP0hMKWwinfVQBeUAzXjC+cPSwl1PRsqKxNc57dW3eMs0f2Y+yQXry7t55DDa3UNbexcOUO7jmvN02tUbpbLHcyMnXQVMrCZalgAFMH5bx0EEj8mG688UZ+9rOfpe3XFnBbHV+2bJnr77l79+6OZf72t7/x+uuvs2zZMn784x+zZcuWjDLf/OY3+c53vsNnPvMZVqxYwT333GNbpx+iFPtZM+tUNjPtv9lnfblDDa2+/lL8xIbKWq79/SqeWF3Jmp1H+PBAI2t2HuGJZJSNJ1ZX8tKW/aljxyJxdhxspKk16lw5iagnj1zyCLdMuiXvISYkwur8+c9/5sCBhK/2kSNHqKysdDz+8Y9/nNdee42dO3em9gNUVFTQ0NBgeq1evXrRp08f3njjDQD++Mc/MnPmTOLxOLt37+aCCy7g5z//OXV1dTQ2NmbUpQ8f9Pjjj6f2W12zZ8+eDB06lOeeew6A1tZWmpubc+qnXAlsZj6jI3irt7drBpDy/NcvZ7rtz2+ngkeaYbZfAbXNba61s4kDJ+YtxDTGjBnDT37yEy655BLi8TilpaX85je/cTw+Y8YMFixYwOzZs4nH4wwcOJCXX36ZT3/603zuc5/j+eefT8sHoPH4448zb948mpubGTlyJP/zP/9DLBbjS1/6EvX19Sil+Pa3v03v3r0z6rrnnnu45pprOPnkk5kxY0ZKkNpd849//CNz587l7rvvprS0lKeeeoqRI0d60ndu8LMSkFcIIK/wIgSQht52ZuWeoUW0MDsX8reZeRFB9MzBFfz91vPyqiMbNlTW8vn5q1J+YiKSFkLc7f088pmTGDRsJN3LSjhtYI+CtTfAnkKFAJo4caJjCKC3337bMQSQiCwEPgUcUEqdZXJcgF8BVwDNwJeVUm/Z1dlphplu0H6Q7+6tT22bqcdO2+1Be08APP1WdZpnP5iHEHeLX51nA/LDw2HmY8BlNscvB0Yl/+YAv3OqsNMJM6vUb3qfs0cWLSmYcDKz0Tn5pvkhnpnXYbpLjPGEAjoFXgkzpdTrgN3qoCuBRSrBaqC3JAJXWHJC2sxizXUFq9vNTKibMu2tmXltR+xS2unekwG0q2vGycBu3XZ1ct8+qxM67RNn96ZoevefGYEFAUKC6X4N4yoDsz+vaG/NbPbkoRn3notypVCgoHe34kT9CCisC4dLzWyAiKzX/c1xqtfsUib7bG/shNDMzCYEnpybCES482BjWtBBIC3g4NghvVi4cgfbDzZ1aqfZKaf2SfWJVXDG0nCISCye6i+jn9mRpjaOxUL0K2mlW1mvdm1/QAKlFIcPHy6YQ63L38BBD3IAVAOn6LaHAnvtTuhUwkzf0fq3k9XM5iM3HJ/YWbymik//aPHxtGoqjkgIpeKEwqVIKJSVMMt3RrMYTrNacEbjvmyIRCJUV1fz/vuugw4HeEyXLl0YOnRoQepux2HmMuAWEVkCnA3Ua+H4rehUwswO4+yk0VXjjmc3p6dV8/B6VnjhwuE3SktLHZflBHRcvHpeReRPwPlAfxGpBn4AlAIopeYDL5Bwy9hOwjXjK051njDCzMjiNVWpdYjRLH2prMhWONmVLYbTbECAE14JM6XUvzkcV8A3sqmz0wgzOwdXo82sx4TLbJ1m7WgvJ+P2ngAICHCDn0cSnUaYadi5Qmj/B3/pF5Zlc6m/EBQjoYlbFq+pYuHKHew72pKKpqFF1tCibfTtXsqts0abLlwP6LgEUTPaAbtONGppcxatzyuhSXtQrIQmTmj2RX1/6wWvlgHgQEMbdzybiI8WCLTOQ6CZ+QgRoWzIGZx8wy/wesWNlwZ9P/lp6UMAbdpdB7h/qJeuqwqEWSciEGZFQG8jM9rM4HiEiJ0HG9P8p4A0n6rScIijLRH21LX4Ip5Ue2NcgK7h9qEuLwlsf52JYJhZJOyWDSml0nyq9NqHnoouJRxsajWtr5A45d0sJHpBf6ixNbUAPSDAz89BpxJm+hlNK6dZM6y0D+O57flFFss1wy4JcC4ELiadi0CYtTNmwsv4JazfdSSlfeyvb7HVPtozjplGMVwzjEmAvSBwMelcBMKsnTGzkRntZVc//DoSNk9o4lR3e9DerhlGjcysL4zarhv87GISkD2BzawdMC5XstoH0P+qOxyTALu5TiGpOtLMrPtXZExOtEbjXDvNPPFIPniVBNhIJKa4edH6tMXrZovZ9Qv7tcXrAyrKmT15aNbrQwMKR4fWzMQkvK2I9AWWAsOBXcDnlVK1yWPfB74KxIBvKaVeKkjLDWTzRrjy1yt5u7q+gK3JD6UUrVHFhwcaM46JSCq/pZcCzbMkwAZ2HWpk1+FE+rm3q+sRIJ78rt6urjcNza23US5dv5ulcz4eCDSf4Gdh5sag8RiZ4W1vB15RSo0CXkluIyJjgOuAsclzfisiYXzGtdP86/ekz8ZuF8Vz4codnl63UH2iRRtJCSqcQ3Prt6MxlUgN2InYUFnLzYvWc+WvV7J4TZUndS5eU8WVv17JnEXr2VBZ60mdZoRCIce/YuGomSmlXheR4YbdV5JY8Q7wOLACuC25f4lSqhXYKSLbgenAmx611xPcJLzV9hn/25VxOhYKC5GoQsUVkbjKKJNN3Pz99d5qUlqf/PeLW6k/Fk3TjsrDkmpnNK6Iqcw1sF1LQsQUlJcIrdE4vbuV0q97Oe/tO5r329zoLtORMc6cayOEfLRs46qMV96v4cm55xREm/WzZparzWyQFltIKbVPRLSU0ScDq3XltFC3vsMq4W0xGXXnC7RF3Qm0xrYYlz/4OmUlIc9saF84exjv7q1n8ZqqtIf2lL7dWP7d8wG449nNGccBrp4ylJ9ePS5t35xF63lv39G827XjUBNX/nolg3p2Ye7M09i2v4Gl66pS2x1pCKpPHKOxcOWOvL4/TUvX6ozFE9cJhFl+uA51mwylOwdg2DB/CZVi4fYx0d7AmqDw4u2uYZbYpP5YxPa41X6vHH8PNrRyMJm4+eX3alI2NyisFlIInPo3F1oiMVfX8YLOOJtZIyInJbWyk4ADyf2uQ90qpRYACyCRNzPHdnQq3HaC2dsx37e7hllik166daJWiU+0/YvXVKW0pj0eDoX1Njf9diG1kELg1L+5MGZIL6rr0vtau44+sXOoW+/+eV2IzqmZLQNuBO5L/n9et3+xiDwADCGR825tvo08UcjnMfHKhjZ2SCJ2v94mdtaQnrbHtf365Mrt+dB3JJuaU//mwsj+3TPqHDukV4Z9rqRn/qGUO7Qwswhvex/wpIh8FagCrgFQSm0RkSeB94Ao8A2lVKYOHGBKPuppq0fOqfoEyal9e+ptj2v71+44bHqs0HSkJVNO/ZsLy7fWZNa5t55399Z7vq62Qwszm/C2F1mUvxe4N59GnaiUh0NEYrmJNK8esXxsZvXHIraL+wtFR1oyVQibmdn5J6LNrOM8BScAg3t3zflcrx6hfGxm+dp+cqUjLZkqhM3M7PwBFeWeJ3YG7zKaF4JAmPmIm87NzGrk9k1X7pF2YpYMWN+u2ZOHZiQHDkliv1n7jWUL8aj7NSqvGU79mwvG8yV5ndmTh3re334WZp1mbWZnQJuNXLquirZonK373Ru2S8LePERmyYD1s6Tb9jcQNzjNxlVi/9qdhzMmBub8n5F8dKjJdG1mXXOEyiPNebfZT1F5nXDq31xYuzNhq0y9+ETYlnx2FN4O/ZyGmfF48bRk8UP01KlTp6r169cXuxm+I+FEG3f9tutWFqYsLPSv6MJZQ3ry7p76tKQjg3p15aZzR+T147Fa1zphaC8+qGmguS2W1t4JQ3vx/C2fMK1r1v0r2H6wyfGaxtBKxu0zB1dwsLGVptYYl44dxPQR/Vi4cgd1LRF6dy3LuGctIYv++Nqdh3lpy34qupR0uEQsZ971d45F0oXIhKGJWVP9d7Xv8Vtp3fdhzm89Edlw3XXXTbYrE4/HefLJJ9c5ZTQXkcuAXwFh4FGl1H2G472A/wWGkVC6fqGU+h+7OgPNzMf07lrKwUbnmTrthdTUGqUJqDsWZfuBxoykIw2tjXknGRnUswuQKcwG9ezC3rpjGT+qRHlzRg7o4UqYOU0qbN3fkLrX5zbt5blNe1PbhwyJVYxLf7Tj2vaxSLzDJWLpUV7CsUj6c3K8370NqOD0YnXz4k2u1/4NcDEJ39R1IrJMKfWertg3gPeUUp8WkQHANhF5Qill+YMIbGY+5tsXj3ZVzspmYbX/ruffzXmB89yZp2XYfMIhYe7M0zLaG5JE+WzqMuJWjbBbrA5wz7J3mXrvy9yzbIvpcW2fxm9e/ZCpP3mZM+96kVuXbHTZiuJg1e9zZ56WYbPMF49sZtOB7UqpHUnhtITEum49CqiQRIU9gCMk3L0sCTQzH6O3oWmJQfTJVo40tXGkOWIa4dYq6q1SilicnLUPzeajj0GmxRzTvPDdrps0q8ssnplxLaZ+e9v+Bld2t9ZonINHW10P2fUJbJ7blFjE8uB1k1yd297onxNjvz8177h9bvfRQ5X5XsujqBgnA7t129XA2YYyvybhhL8XqACuVUrZGuQCm1kHJt+4bHb2LD36ZC9Hmtro272MUYMqfBE40a3dzQmnkOddS0Ns/fHleV/HDYXqbxHZoJSa6lzS+vzrr7/e0Wb2xBNP7AIO63YvSC5f1Oq5BrhUKfW15Pb1wHSl1Dd1ZT4HnAt8BzgNeBmYoJSyjFwQaGYdGDP7VTZ5CNykgbNK9rJ2Vy1L1lUVfZG3W7tbvoTbyeXA7/3t0mZ20GECwM0a7q8A96mEtrVdRHYCZ2CzPDKwmXVgzGwi2fj5bK6uZ9b9K2wD+s1/7aPUkhijXSQWTxwvJnNnnlYQ3zUjsbhKC6i4obKWO57dzJxF67nz2c2eBUT0e397FJxxHTBKREaISBmJgK7LDGWqSK4yEpFBwGjANiJpoJl1AuxS4dnZzpojMT480Mj2g00s31rDU/My3/pO4bQLFW47W7JNB5ht+WPROP/Ysh/IDP8tIp5pTX7vby9mM5VSURG5BXiJhGvGwuS67nnJ4/OBHwOPichmEvNAtymlDtnVGwizDsz81z5KhZ22wq2Aiyv4yv+s5fbLz0ybFLh22rBUzgEzih2CfP5rH6VCcWdDLp7qab5uhu1YHL61ZCNjT+qZVyIWv/e3F8IMQCn1AvCCYd983ee9wCXZtC0QZh2YnQczE55o2GlkVuXrj0UyZjm1/wtX7vDcAdcL7PpATy6JnLM9p/pIM9VHmvPS1Pze38WM8e9EIMw6MHbRIqx+gG723/X8u8DxH9bowRVMH9kvY3Zt9OCKXJvuGW4jZuSriWVbPhaH/3p2M2UloazDe/u5v73SzApBIMw6MGUuZiOzRSlFNKZSGtrowRVc+/tVGaGJ1u6q9UUauEL0gVdoYc1FxHV47w2Vtb7ubz8LM/8+CQGO5Gs/MRty6mfPFq7cwdNvVRONm8+u+SENXLFtSE7oZyLd9JXf+zuImhFQELRh4E/+9h7NbeYBfbO1nempPxZxDPJXqCCAbtGnDWxqjdGrawn7jzq3qT1saMZybvrK7/3tZ5uZf1sW4IovnD2MqydZZ/PL1namp1e3MscAf4UIAJgtXzh7GOv/62K2/vgyLjpzkKtzctEi8tU83PSV3/vbz5pZIMw6AbMnD3Udz0wpZaph6Pdr3HTuCNu6S8PC7MlDc2x1YcimL5ww65Ncy4dD7vrK7/3tZ2EWDDM7AVNO7cPSOR9PrefbuLvOMpeAG01NgCsnDkkN4X70mbNYuHIHB5taicdhaO+uTD61jy/WZhrR94XZ4vU3dxymocU2+EKKfH3XTu3bjQMNLakYaW76asqpfXzd334eZgbCrJOgj1rx+fmrWLsrt+U1SikUiUgR00f0Y/TgCu5+fnNKOIoIHx5s5CdXjyv6D8sKfV8YyadvsmXX4cSa0WOROHc9v5nRgytczWb6ub/9PJt5wgszfZJUsxDGWgQDY7gbPfrEt9n4E2VzjWwYNagi5x+s/mH8zasf0rU0nJpd09DPqnnZ7vYgn77JBeNspr5/zJ49/WymhtbfbvrW6Xn26n5yPV5ITmhhpkUo0EZk2jIS7QHQfH6iuihKRl8fLXJpgnrX/kT6NjhdI1tmTx7K0vW7iTqkrXOKsLGnznod4MbKWp5cV+Vpu9sDt31TCPTJiq2ePauExm4SHTs9z17g52Gmf1uWJWk2H4vPRhJvwfSHeuHKHYbj6ecYfX0WrtyRZvR160+UzTWyRbMbfWxgD9ty+bxFq+uaPW93e6D1zRfOHsYlYwbxxbOHMbR3F8e+CIfg4jGDmDC0F9OH5yasW3Ud9vRb1Rjl6cKVOywTGrtJdGxVp5cEEwA+xS4h6+I1VTy3cY/jefkkYLW7xj/fq+HMt6pTxmM47kt16dhBjlFPp5zah//+7HiuXfBmSgux8znL9iFsbDH3a3O698VrqmzvwyzZiNdDJaNNLV27NqdrSZhHbkjENdxQWZvWr275+Mh+qes9tX53xvH6YxFG9O8OJvHZnNLp3bpkI8s2GUOC5Z9g2EgwzGwH7BZQGxERlFIMqCjP6Pxe3coykl4Yy9Q1t6WVNyYdceML5HSNfclQL1qCDX05t2Gc9TN7z2/cQ5PBsXb68D5sP9jIkabMB97JQTTuom+MGO/ZeB9WyUagsMlFjGGnAf7xXk1aGX2CZuOM6Y5DTWw/kLng3dg3FV1Lbb/3Y20xS5veqEHW6zJvXbIxLYmLsU4v8bMwcxxmishCETkgIu/q9t0jIntEZFPy7wrdse+LyHYR2SYilxaq4V5gliT1pnNHsHRdItmHldqsV/mtErA64XQN7ZjV9kvJ2FpOTDm1Dz+9ehx3fnJM2v5wSLjt8jN55IZppklFnIYMbvrGiNk96+/Dqk+0/YXkC2cP4/lbPsGCG6aaJloxfs9avy64YSr//dnxGb5hZomSZ4zsZ/u9t8bitkmWrdD60KpOL/EoOGNBcKOZPUYiucAiw/5fKqV+od8hImNIRI0cCwwBlovI6Uopb18PJmjaFqRrY/r9RhIJbdPLrt15mDajMciAPlKDMfGtStbrZAR3uoadIINEerFssEt68eTcc/ivZzenpWzLdThqF8XCLMy3/j7s0ti1J9km6jXzbdt9pJnXPjgIJPpSS5Rs972XiNgmWbZ6pmJx6+GuV5nuNTq0ZqaUep1Emic3XAksUUq1KqV2AttJpJVqV9x2qJWmEHF4m2nHNfXeeE0nTWLxmiq2u4zDBQltz3hHVmnoFq+pYtb9K5h678vMeuC1tJRya3ce5oOaBjbtrk1lvIbED+WDAw2p+8hnCdS2mgbLNHZmYb7192GXxi4f7PrEim37G9hYVcsHNY2pjOFO5dfuOMxbu2tZs/MIb+5InGN8LuyerZISsdTcfviXLZbttrPcNbbFPE2V11knAG4RkRuA9cB3lVK1JFJIrdaVqU7uKwpmWpl+n9lbMhZXjkkySsMhbl2ykWctjPd2b18nW5nG4IpyBvXqktKiAEdNwc7mtHbn4VR79YluAb7/zDupz24eRruJhDhiaeeacmqftNRnxvuwS2OXK7nY4fTnuLFRml3DjLZo3PbZKg2FTJ8dpRStUevkxFbJoo22SS/ws2tGrsLsdyRidKvk//uBm8hUIMDixSEic4A5AMOGFc64azfMNHtLKhIawivv12D1Eo3E4ml2Cjf1aujfvHb06FKSkQZOm03Ltu6l66r4oKYh45heg8zmjepGc1u6rspUWEw5tY/tfdh57+eCXZ9YCTOzc+xslG6/00gsztyZp7F8aw1mI8PG1ih9TZ4dp3799sWjTWdjc7GxOtGhh5lmKKVqlFKxZFLORzg+lHSTQkqrY4FSaqpSauqAAQNyaYaxPtPPdowc0COj83t3LU1pCCf3NrfVjBzQwzb1mN00upv0bk51ZFt3eUnItL3lySiohXgA3d5nobHrk2zOsfu+rWx6xn4dOaBHSjsd0CPz++3VtdTxezdr2xfOHsZPrx5Ht7Kw5XlepcrrdMNMETlJKbUvuXk1oM10LgMWi8gDJCYARmGT567YmL0l9TacMUN6ZXjBhyRx3qrttoliTNlQWcuOw+5yPNYcbeXKX6+0HGaOHlzhybIVO00hHz6oaWRDZW2GluWH5WNeY6bJC+lDEu250Rjev3vG0PDbF4/muY35OB1bf4nNkRihbr3751E54G/NzFGYicifgPOB/iJSDfwAOF9EJpLovV3AXIBkuqgngfeAKPCN9pjJNGmzq/AtVjYcq0Ss+i8qZjOz12pi9zDWqT/XzAa163ATHE6Ue/m9GmLRKBIOIyKpZSra/WlZzc2uq+03a29rNJ7WBzsPNhKJKSqPNJvaGo3Y2c7qjkW4Zv6qtPR1GyprueZ3q4gZ2q1fPmbsI2NikMVrqtJsfFYp8vT3mM1+q2Mxm2fJyta3bX9DhtA19oGG1o927bI67mSHVUoRU1DSs/+ptpW7wM82Mzezmf+mlDpJKVWqlBqqlPqDUup6pdQ4pdR4pdRndFoaSql7lVKnKaVGK6X+Xtjm54bT28MuEWtcJY5fOnaw5flmoZyNdWbTTgUpQWY8ru17YPk2yxDS104bZtperbxmx/rvz02gKRJNqzuX2U1jX2nMf+0j4ibt1h839pEx8a1WXitjvIbVPbrdb3XM7vuGdL+ze20iXBj7wO13aNc2Y58Y8XL41+mGmX7EaeZSj7ZvQ2UtVz/8OhJO2Bo0TcEpfdnOg43c9ImRGbOZ2hdplkXHrk6rB9CpjJ6mlmhamjLjciDt2Etb9qeWSOmHeEaNwUrjctseDf19myWwbdLFFrNKcKvf39QSzbi2XWJcuz5xOieb5WNG9BqkiKQ0SLvn4GhzhNGDK9ISDGuICIL5s1Xf7O2SJTu8GmaKyGXAr0gkAX5UKXWfSZnzgQeBUuCQUmqmXZ2dRphZYSfk5r/2UYbG88DybQzobr8cqW/3srS3oZH5r32UMWPnNiVarmiag15wGbH7Qeo1Bityeevq79sswa1e47FKgKvXRi4dOzjD1cBJm7HrEy/P0WN8PjQN0u456Nu9zDapscL82WpPXciLYaaIhIHfABeTmDRcJyLLlFLv6cr0Bn4LXKaUqhKRgU71dhphlsts5oGjLRkPTVNLlJN7dbU4I8GoQRVs3F1nedzs7evkiGuke1mYuFKWC82nj+iXl+ZgxE67yQf9fTtpPG60KK28lYbpF6w0SLvUeIN6dnE1KjDSnsGMPNLMpgPblVI7kucsIeFw/56uzBeAZ5RSVQBKqQNOlXYaYeYWvdHaSlOYPqKfqYYAx2O5P79xj6VzrNn0et/uZabREKwYO6QnT847J7W9obKWScP6cOBoSyoCrH47X6y0onwx9oWTxqMlwNUM6WbDqukj+rHzUBODenbxRWJcM+w0SKt+vnbasMRsps1zYvZslYdDRGLtM8/mUpgNEJH1ut0LlFILdNsnA/qwIdXA2YaqTgdKRWQFUAH8SillXFKZRqcWZlYRMzScNIWFK3ew72gLrZE4FV3CTBveLzUr1bUsTKNFRILe3bL3EbNDb9vTZjP1dhXjrGAu6PuirjlCWCCmjs+EtQfGBLgiYhoMUz+bmW0wzPbCSYN84OVtHD0WpW/3Ugb17JKaSV+xzV4B8frZyhaXw8yDSim7ZYxunOtLgCnARUBX4E0RWa2U+sCq0k4tzNxgpSk4aRBmoX80zEIAuQmuZ1X+6beqM2x7RrvKwpU78h5umd3znc9u5gkXaxmteLu6ns/PX8WRpjZKwyEisXjqv6ZlaMe0gI/6+zKGjNaCDepnO92GlG5vrIb+ds9WLqnmvI6MYYdHw0w3zvXVJIz+TUCTiLwOTAACYaaRTeowO246d4TpEhKrlGLZDjP1w4lDDa2OD4nXQfg08gkzrZSiJRJjzU63cQoSGO9VHzK6/lgk43ixE+N6iV1/W6Wa8+Ew04l1wCgRGQHsIRFp5wuGMs8DvxaREqCMxDD0l3aVnnDCzCv0IXXaonEiscQiYiuP9GwTaWgrANqicfbUHXMs36tAww9jGrv39zdwtECp2qzQa6m5BsPsKBj7+0hTG327lzFqUIXlovvBvbuaBocsBF4IM6VUVERuAV4i4ZqxMOlwPy95fL5SaquIvAi8A8RJuG+8a11rIMzyIpvp+2yTjGgrAIxYPSzG4IFeol/87TZkdC6huK3QuzMYNWK3yXU7EtkutrcaJRQCr1YAKKVeAF4w7Jtv2P5/wP9zW2cgzNoJszeu0YZUXXuMY5HjwwU7YXBy7y6MGdKrYCnFrLDSHOC47WurLlaaF5i5d/hpbaYTxrWk4BzOSY/TetUvnD2MtTsP89KW/XQpDTG4Z9e056o0HKKsJMTuo4cq872XDr02M8A7nN64037yckqYOT0UbdG4Y0igQuF0H6ff+QKtUe+mQI2Opvk6tLYnxlSEy7cmcgtoy36d0sGZpSJ8akM1f7p5RtoMr+YGciwS52hLg+kMr3yzLvvoCAb8LMz8u2r0BKRXtzLXD0OhbGRe4LUnR7YOx37CmIowrsiIUGKXDs4sFWFbNM7qHcdtEMbzs013mA0dPQdA0UgttPZoBtLvZGP7KKSNzG06uGPROGNP6pkx1CsPm0dMzZVcYrv5BTezzHZlrGZq++heZvmkO8yWQDPLETdhfDoTWnIUuz992UKgDYsONrTS3BbluU1702LIa8c/PNBI9ZFm/vFeDZ///So2VLqfqT2RvlM3qd7syuhT9+n77VWdc63Z+XYp//IhiJoR4Aq7UNxWZb3GKWS02XGj46qTE2e2D/zaXbV87I4XUqsStP99u5f6dm2mhhuHVrsyercUfb/p12iaxVrL1knbLYFmFuCKHuUlrh+GbFPNucUsBHRmOrhMxOJzvmgaaSQWpyXpz6f9rznayh3PbnaVbalY9O5a6vid9u5aanls5IAejvvNYq1ZnZcvfraZ+VqYFVttbW+s0seZcbipjdF3vsCZd73oaSqxXNLBGT3TvRxEuhnStEeS4FwxfqchgfNG9bcto8dN+r0Hr5uUVmdpOP/0fFYEw8wcOZFsK5CeUNiJmNKSv8YdU6Flgxfp4Nw+zl451voleYoZVn5xi9dU8fd393H5WSc5Jhh2k35v0VfPZkNlLat3HGbGyH4F873z8zDT18LsRCMbm5nVuV6Qbzo4t68grx78YkeSyJYNlbW8uu0AR7NcT2smyDSHXP2Suhkjj4eE0ieQ8SKhiZ9zAATCzEf0KC/hWCQ3l4ZC2dBywWvXjI6MMWTRy+9pTrPuwjcZk7w8saaKeeeN5PYrzsxIZALw4YHGVIhuIC0UuhcJTQLNLEdOND8zq2Subs/1inzSwS1eU0WTC3cELymUG4IXGMNna0+yMVS7lTDTJ3lJ7Xt9B8P6dbcM3a5P8uIUCj1bAmGWIyeKENOTyz17+QDlkw5u2/6GNC2kvSiUG4IXuEk2YlfGKoz20nVVOZ2XL34eZvq3ZScgxpRhbv/05+ZLPungnFKeFQo/rxBw0wt2ZawSoLRF447nFaJf/DybGQgzH2GWBCObc70gn3Rw+bQ/H/w8ARAOOf/Aw0ZfGB1WCVDKSkKO540a5H1+hECYBbjCKdFsoc7V45QAN9tEw+2Bn21mbvok24TS2n6n82ZPHprhM5gvfhZmvraZnWjok2DE44pITJkuVTEiIp5kaIL80sHpEw3H44lEKF1LQ8SB1kg8tQwpagwbkSd+tpmZJTYB9wmGs+1v4xKvp+adw31/38ruI81UN9fn7b/jZ5uZozATkVOARcBgEpMjC5RSvxKRvsBSYDiwC/i8Uqo2ec73ga8CMeBbSqmXCtL6Toj+wZ51/wo+PNDo6m3nRUITDad4YbkmGtb4/PxVWYUQd8LPNjMw75Nsviu79Hta2j29n5kx/V7vbmW0ReOoaCTvUBodfTYzCnxXKfWWiFQAG0TkZeDLwCtKqftE5HbgduA2ERlDIkHBWGAIsFxETldKte98fSfALHmHXdmAzodd+r1t+xtM/cy09HtA2szzCe9nppTaB+xLfm4Qka0kknheCZyfLPY4sAK4Lbl/iVKqFdgpIttJZDB+0+vGd3bs0tmZlW0vnPzMnEJaez0s9PMwM1+04Ix6IaGl31ubDNBoFCD64IxGH7V86dDDTD0iMhyYBKwBBiUFHUqpfSIyMFnsZGC17rTq5L6ALDlrSE8+rHGOpy8inDWkZzu0yFlT0PzMRMQyQa+Vu0GueF2fn9hu8f1vr2lgf735zLLdefnSoTUzDRHpATwN3KqUOmrTaLMDGRZfEZkDzAEYNsy/8aiKybt7Eg6rbh4QrWyhyUZTsErQ63UY7O0HG/n8/FX07lZGXXNbKslK725laZqjplFaJWJxShdYDKy0ziNNba7joHlJhxdmIlJKQpA9oZR6Jrm7RkROSmplJwFa6Es32YpRSi0AFgBMnTr1xHP1d0E2drBj7bQW0ioc86GGVtfhm7NNiGzEGG0jElOZEwq6+peu382PPnMWdz+/OT2evkkbth9sstQoi4FVX/XtXsb++hbL5L+pSZE8+tkMr4SZiFwG/IpE3sxHlVL3WZSbRmKkd61S6s92dbqZzRTgD8BWpdQDukPLgBuB+5L/n9ftXywiD5CYABgFrHW6TkAm2djMmloizLp/RSph7NghvXh3b73rRLJusUq4O6Ci3HWC3mwTIhvJ9u0fjSkWrtyRkRjECiuNshhY9dWoQRUcaY5YJv/VHGa9nDUGZ5uZm+V4IhIGfgNcTEL5WSciy5RS75mU+28SyYKd2+aizLnA9cCFIrIp+XcFCSF2sYh8mGzUfcmb2QI8CbwHvAh8w2kmc8OGDab7i6my+oFskpbUHYuy/WATa3fV8kQymsLiNVWs3VWbtv/aBW9mFa/fyOzJQykJp38vWnBGY3utEvSa1VFo7OxLZhTK5pQt2fS3htbvsycPzQjsmC8eOc1OB7YrpXYopdqAJSQmDo18k8SI8IDJsQzczGauxHr52EUW59wL3OumAQBTpkyxurbbKjol+sB+WgBCvX2n/ljEteamodm3ctU69EmAjbOZWp1Os5lOCZHLSkKM6N+dnYeaTO+7NBxi+8HG1CSEG9zE4tfjlxlSt/2t9zPT9/uTc48H2vRREuCTgd267WrgbEM9JwNXAxcC01y1zQ8CY+rUqWr9+vXFbkaHY86i9fwjGR8rGy4ZM4gFRUog7BW53rtbQkDX8jDxOI5e+h0FEdmglMr5ixeRDc8+++xkuzKxWIzPfe5zuwB9+rAFSRu5Vs81wKVKqa8lt68Hpiulvqkr8xRwv1JqtYg8Bvw1b5tZgH+xsl8V6jw/MbJ/d8Bd6G39C9uN5qCUIgY0JhfYexmWvKPjUjM7qJSablPMzSThVGBJsr7+wBUiElVKPWdVqS8cdAKbWW7kupD4b5v3+TqjkRvetHAYNSPbRdBmNqBCpfbraHhkM1sHjBKRESJSRmLF0DJ9AaXUCKXUcKXUcODPwNftBBn4RDMLbGa5sW1/A3GVfT/VNUdSEW39nHPSDjdhub1KmAIQDl6sgDezmUqpqIjcQmKWMgwsVEptEZF5yePzc2mbL4RZQG6YJeTN9vyOKszcON56qdl3KQt7VpeRxWuqeODlbdQ1R0yjjHQtDTGoV9dUtIxi4pWfmVLqBeAFwz5TIaaU+rKbOgNh1oFJJOTN3fPfKqFvR2DkgB5sz9IhNFtNTV++pUB5DYxJSfRCWnNBjsTiNLQ2+kKb7vArAAqNnc0sGGpaM3fmabzyfg25rA6S5PluMFsG5IUDbj7MnXkay7fWkE1oNLMfmp2A0+/P1rXDLdlq18XWpgNh5kBgM8sNLUHs/Nc+4uX3ajL6S3uwrPa7wZjARGPtrtpUIpNiCDR9suKdBxvZU3ssqyVdPbuU0NwWJeJS4Sov0GL2bBMYFzvhsRc2s0LhC2EWkDtpmcRdaBgaikQSErtkv2Ce6kxDS2TiVIcficYV8SyUrea2GJc/+DqRWDxtgboftNT2JNDMAgqKVZIRL85zKpPrtfPFmBJPw60fWVOr+wQwmt/Ze/uOmh7PR0ttzTJAQLblvcbPwswXfmYBubOhstYyg48TNUdbmHX/CuYsWm+5XtMqoYbb44XCmBIvXz+ybMqbna9pqdmSbf8Vq781QqGQ41+xCDSzDoyVPcvJVqbt35dcfL39YBPLt9bw1LxMzUKfUGPf0RZaI3FfuAoUKsmtHsEkEJ8NubRJn0BGc80IhYVIVKHiikhc+aK/NfysmQXCrANjZ88C+wfLeCyurO1fTglOioFX0WWtZjP7dislphT1x9wPR6uONPP5+auytqH5sX+t8LMwC4aZHZhc7VVWD1yx7F+5kE20WqVU6s+IlbvGkeYIdc3ug2MqpWiNxlmz84gnoZb8ip+HmYEw68B4bT853NjaYdZsZpNezs4+ZiXgsnWwNV5DC7XU2fBobWZBCIaZHZjRgysICcQ8Sqq7u/aYL7zM25Nchuhu8UuARy/x8zAzEGYdmKffqibucSox8DahcKHwS/BEO96urmfW/Ssck6ToU/cBqaQsWjBKP/myBcIsoCBYJRfJl46QUDjfpChGvIywodXXEonx4YFG2yQp+tR9xuV7+u21u2pTKf2KKdCCFQABBaFQQRaPtkRZvOb4GkA/rs3MNymKEa81CmN9VklSjKn7jOfpt/MNee4FgWbmQLDQPDdmTx7K0vW7iWYRC98NrdF4ynY2enAF1/5+VUaKtmJrCrMnD2XJuqqcFtm7wWtNDcw16WztasW2w/lZmPliNjNYaJ4bWrKLjw3o7qq8sT+d+nfhyh0pzcGMYs7YaYvsLx4ziI8N6E5ZjtmerPqgvX6U2dr+im0rDFwzAgrGlFP7MHJAD1dls/2B1h+LWKbl0iiU3c4NU07twyM3TGX5d89n4im9c6qj2KHZs3ExyaW81/jZNSMQZp2AXG1nTg9er25ljnkG/JIcRUt662esEiJnQ7HvMxBmDgQJTfKjUEl1bzp3BFNO7cNPrhpnqqFpyWj9QCES3npJWUnIMiGy23b7ob/9LMyCvJmdBP2M445DTRxqbMs6WOOAHmX071FOWUmIa6elrxfcUFmbCoToh9lMM/RtHDmgByP7d+fNHYdpi8bZcagpq/A5PZIx/xuzCJfdJZm8WO8bNnZIL2qb25gxsp+tn5m+b+F4rLTe3RLb+uS/uSIe5M186623HPNmTps2bZ1DqrmCEAizTsiVv17J29XWuQGsZuomDO3F87d8opBNKxqz7l+Rdc6AbKNmAFSUh30T4cKIF8Js48aNjsJs6tSpRRFmjgquiJwiIq+KyFYR2SIi/57cf4+I7BGRTcm/K3TnfF9EtovINhG5tJA3EJCJU6ISq6FAR05w4oRZlA2rF7m2KD1usTjd7pyjLVG2H0gkH+ko61yzwc/DTDej9SjwXaXUmcAM4BsiMiZ57JdKqYnJvxcAkseuA8YClwG/FRHbPF2Bzcxb5s48zXEW0ozXPjjIrUs2et4eP2AWZcNuXWa2P06z8lqyks6EV64ZInJZUtnZLiK3mxz/ooi8k/xbJSITnOp0dJpVSu0D9iU/N4jIVuBkm1OuBJYopVqBnSKyHZgOvGl1QuBn5i3b9jegyL7/WiIxntu0F4AHr5tUgJYVD6/in2WDm0TFHQ0vnGaTys1vgIuBamCdiCxTSr2nK7YTmKmUqhWRy4EFwNl29Wa1AkBEhgOTgDXAucAtInIDsJ6E9lZLQtCt1p1WjYnwE5E5wByAYcP8ZVvo6OSbHPilLfu9bI4vyCb+mVfsONTkeqH5/Nc+4sDRFq6dNozRgytSC8+9MPx7iUcrAKYD25VSO5LnLCGhBKWEmVJqla78asBxGte1MBORHsDTwK1KqaMi8jvgxyRspD8G7gduAtMRToaKoJRaQELaMnXq1EAF85B8kwP3KPfFKjdPySVpsBVuljoZF5pbhSU3JmZ5u7oeAeLJbREpako/Ix4tND8Z2K3brsZe6/oq8HfHtrm5soiUkhBkTyilngFQStUopWJKqTjwCAlpqzXsFN3pQ4G9dvUHNjNvmTvztLx8rspLQmnG6w2Vtdzx7GbmLFrPnc9u7pARVK36JBdThsuhVJoNTQtLbsSYmAUSb379dq7JUgqBywmAASKyXvc3x1iNSdWmX4SIXEBCmN3m1DbHV7AkWvcHYKtS6gHd/pOS9jSAq4F3k5+XAYtF5AFgCDAKWGt3jcBm5i25JgfW0AdpHD24Ii1pit80BbdofaKP/pGtplYWFtryWNRvlvDEbRKU9kjg4gaXw8yDDq4ZrhQeERkPPApcrpQ67NQ2N+OJc4Hrgc0isim57w7g30RkIgmJuguYC6CU2iIiT5IY/0aBbyil3HseBniCm+TATsceWL6NSaf0yUia0lGT/045tU9av5x+5wu0Rt0Jp77dShnUswtb9+cetaKxJcqVv17JoJ5dmDvzNAD217vLu1AaDqUFcbSypRntb177unlkM1sHjBKREcAeEt4PXzDUMwx4BrheKfWBm0rdzGauxFwtfMHmnHuBe900IKBwOCUocXrwmlqilnV0pOQnVoRD7l0vIvF4zvlJNfY3tLIv2W8vv1cDpNvG7DjaEnHUkM3sb+BtCHQvbGZKqaiI3AK8BISBhUklaF7y+HzgbqAfCdcugKiTw6+PV7MF5Eu+CU+mDe9r+QO+dtqwDm9Lu3Ts4Ix9Vj/GaFRl3Z9WyVJEJOE6Q7ptzI765khKQ7aypZnZ3x5Yvi2rNjvhldOsUuoFpdTpSqnTksoPSqn5SUGGUuprSqk+Oj9Wx2FA55u2CkhhTODb2JLIAak9cOeN6s/qHYdN1yyKCG98eCilOej3a49rR7elab50L23ZTzyuaLWxh8XJ7M+mVnvriZcTWK0WriV6W5qZXa0+i3R5bvBomFkQAmHWydESzE6992UaW6JpD1t5aZiZowfyjy37TR9CTXMw2//A8m1Z2dL0i6lLw6GMxezFCs2tdw62s6F1Ty481yfsvfzB1/OyoWVDzELQas7Ai9dUsbv2WMZxr0VLkAMgoOg0GQQZHLd75fI2bWoxz/RtZkvbUFmbpsVpaDYdLTR3xPCDbe/Q3HY2NLPhdjY2NGWwjRm3nYhZlI/E4ixeU8Udz242rdNr0eJnzSywmZ0gmNmHrp02LGe7mll9Wp1Gnn6rOs3eo7et6ENzmx1vz9DcVvcEcNXEzBV8bvtOc7LV/9CzsS/Zle/bvYyFK3dYlin3eBmXVzazQuBrzczJHyrAPXr7UEWXEm6dNTptluuBl7dR1xwhLImhxKVjBzF9RL+Ufag1EqdraSgtvI12vK4lQu+uZZZhb+xCa9cfiziG3m6v0Nz6PupSGqJnlzJaozGumngyt19xZkZ5ow2tNRLP0D5zIdvQQ3buHV4H7fSzZuZrYXYiC7HFa6pYuq4qzScpX/8hq8XjejsQHLdvLVy5g77dyxg5oEeaX5M2i6llChrZrzujBlUwerB5SGe70Nq9upU5ht5uz9Dc2S6w1/edVdifbH/g2Tz1R5raLCcHgLycfM0IbGYOOC1nOtGEmmYDSVDP8q2aT1Jiz9vViWOFCP6n2bdSvwGdl/zS9bv50WfO4u7nN2eVem7skF6W15t1xkCG9bPPLmV3vp8oRnKX0nDI1sjvtR3Jz5qZL2xmdsuZTjRBBolhi/6+4wpihuGLZifxkg2VtXxryUasXubRmOJnf9+adeq5d/daL3r/6FCT7XGn8/1EMZK7RGJxW00uZJeNJgcCm1lAVtQfc/YNclMmGzI0MgsaLGYxNcyS1NppLO/tO8pYetrWWcx0dtlgpUEWIqGwhlEzM16rxGNh5udhpi80s4B0enUry3j4jdu9unmbP1GbccwXsyS1dhpL19Kwr2xm+WClQRZSWykrCaVpZsZrNbbav3yyJdDMArLipnNH6GxmpPJW6keaN507wtNreqX9mCWpnT15KEvWVWFmp77p3BGMHlzB0vW7iZpI03zSqy1eU2U7G2vGrUs28tKW/ZSEhDgQj8OlYwdZTgzor+G0IqAQfFjTYGsza4spbl2y0bPIwX62mQXCzIdoPzSvZzPt8Er7MUtSqw9JZLUCYOmcj3u6AsDoSAoJ+1JDa2PqRWHsw1uXbOTZjXsy6rIKJW52DcjvB23q+GoxTFVK0dQWs7yeVpfWfi/w8zDTF8LMbjbzRJwAgEx3CaCgIXfstCe3WCW6hYRAs2u/MTxPvjiFDl+6riqjf5dvrbEsbxZKPN/w5GaY1WVVfzZakleh0APNzIEgOGPxsdKePj6yHx8dakpL/jt2SC/e3VvP9poGWqNxPj6yHxVdS20T3bY3TslEzI6HsNaCjLPJkH94ci8ZXFHO0dYIzW3m923W/lwIhJkDgZ+Zd2g2HDOvfM0eZFwBsHhNlasVALuPNPNWVR1dS/cwqFdXzhrSk3f31PPnjdUca4vxy5c/oG/30ozVBcXAKYGJ2fEWGwFo9gTOnXlaIpJvto0rAD26lNCjS4ll9Fyv2hgIMwcCzcwbjDacQw1tKfvQ2p2HU/agY5F42gTD9595J/U54fAR57lNe3lu0950m5P2P2l72n6gMeM7OqC7ZjEFmlNquWxTz5n9RN2m9DPav/T7sl1wboXT/XglYgKbWUC7YGXDWbquig+S/l/6Y/oktV7YZYzXLKYwy0Uzs8PsJ5qLzczJ5SZX2iuVnp81s8DPrBNRbhGSprwkRNjkISsvCTGoZ5eCPIBWbWkvzFxE9FqD2XGzfRpm0ScSNrP2xUrz6du9zLb9dseyIfAzC/Atc2eexvKtNZjZh7ON3mCH2+CMen83Y9IObfF9WzROJBZ3TK6bcT8OP7Rjkez8xObOPI1X3q/JeQa4S2mIlkh2J1vdw+bqesTm/XGgoZVQt979s7qYCcEwM6BdMAt/re2PmczStUbjTDm1D0/NM/cBGz24gpsXrTf16s+mLcZEGxpOwRn1obi37W9Is+0BfHig0TK5rlVf2B1vi8YthYWx7WCevq40HCISi6dmfptao6a+awDH2o4LTzcajZ2/WXNSEFsdjyko6dn/VMeLOODnYWYgzDoR104bloqoYdy/dufhDOdJLbignQ/Y9y4ZnTZZkE1bNPSJNow8sHwbl4wZnArOaEQLxb1xdyJZirGMllzX2H6rvjBrn8alYwdbOphaBW40+sfpNcwPaxpYX2ne7lzQTxoYAz1qmAm8QvvBFepa2RIIs06EPlig0TVDO2YVnNFNncZlQZprRl1LhGNtMVojcVPXDLu0dE0tUceZtgNHW0zDftvVb9duq+VM+uCMbpcz6bEKD+4Gq0xOZrOddgKj0MIkEGYB7YbZygGNXNfn2dXpBjst6dKxgx1XH1hplvrjZuTS7nzWMM5/7aOMJC9ucTOb7Cb6RiEjdEBgMws4wdEEipljriY8tNUHa3cepqk1lqFF6TXLeDxhA/KLg65GoRMjuxFSHUUzE5HLgF+RSAL8qFLqPsNxSR6/AmgGvqyUesuuzkCYBbQLbrSkARXljB5UkbbQXB+Ke/qIfuw81JQ2m2kVqrsYONnpOgNeCDMRCQO/AS4GqoF1IrJMKfWertjlwKjk39nA75L/LXEUZiLSBXgdKE+W/7NS6gci0hdYCgwHdgGfV0rVJs/5PvBVIAZ8Syn1kt01guVMJzYbKmu59verbENxb9vfkDERsf1gE6+8X+ObxMNGm6VmRwwLxBSp/11LQ8SB1kicWFz5YjmUWzwaZk4HtiuldgCIyBLgSkAvzK4EFqlEhatFpLeInKSU2mdVqRvNrBW4UCnVKCKlwEoR+TswG3hFKXWfiNwO3A7cJiJjgOuAscAQYLmInK6UsnTiCZYzZeKU0GT04IqChQTS+4T17V5G72TSEX1CE2O4HiDNNcHoA2asUx/eR0s1Z4YWinvtjsOm9qBYPBFY0g/CDMw1UKN/nL6/IjFF5ZHmYjQ1F1R9fT29e/e2LLB7926AHiKyXrd7gVJqgW77ZGC3bruaTK3LrMzJQO7CLCkZtbzvpck/RUJynp/c/ziwArgtuX+JUqoV2Cki20lI4jedrmXkRNXM3CQ00Tu0epngJNeEJkb0WhOQUefaXbUpHzKzUNtpddU02KZTczq/mKR/l0ksFoOHQ9CtNJzS2sICrR5nV8qTRx566KEpd999t2WBn/3sZwD/oZT6m009ZmNR4426KZOGK5tZcoy7AfgY8Bul1BoRGaSpfEqpfSIyMFn8ZGC17nRNohrrnAPMARg2zPxHeKIJMQ0toYkmzOMqc5bK2DMLV+7IW5i5SWhy97J3LbUoI5rWBJjWqfmQOTnlOqVTy8Wpt1Do/czqmtt4a3et6xnGPt3KmDysD3XNbSktd/vBxgxn4iLyP08//fT8b33rW6ba2a5du9iyZQvACw71VAOn6LaHAsapajdl0nC1gE4pFVNKTUxWOF1EzrIp7kqiKqUWKKWmKqWmDhgwwE0zThhySVaSb4ITTSPbU3vMtpxZaGs7ttc0OCY0cVo32Ld7mW1mbq/WHeaLZvtbvKaKf7xXw9pdtUSzWCF1qLEtdd72g01s3d/gJ0GGUqrt61//Og899JDp8Z/+9Kfccccd2mjOjnXAKBEZISJlJMxSywxllgE3SIIZQL2dvQyyXGiulKojMZy8DKgRkZMAkv8PJItlLVED0nGT0MTsnHzwKqGJkSNNbY4JTcxCbesZNaiCwb27WvaB0/nthZXtr5iOpF4zb9688qeffpq6urq0/ZpW9qlPfcpRpiilosAtwEvAVuBJpdQWEZknIvOSxV4AdgDbgUeArzvV63hhERkgIr2Tn7sCs4D3SUjOG5PFbgSeT35eBlwnIuUiMoLE1Opap+sEHMeYrCQkx5OaaBh/HvkmOClUOrfScCKUtpViddO5I5g9eSglYfMfvJbQxOr+wqHcE554jds+9MJ8UiwTjJV2loVWptXzglLqdKXUaUqpe5P75iul5ic/K6XUN5LHxyml1tvX6M5mdhLweNJuFiIhRf8qIm8CT4rIV4Eq4JpkI7aIyJMkplmjwDfsZjIDMnGT0MTr2cxCpXOLxOKeJDTRZivziZpRaNz2oZfrNIvBvHnzysePH9+q2c40reyRRx4J+XoFgFLqHSBjjYdS6jBwkcU59wL35t26ExjjFP/iNVVsrKqlqTXG2p2HAdK28xVmXiQ0MUOzZzklNNm2v4G1Ow6n1lEebmqlt2HovHbnYT6oaUitADjaEmHb/gbfCLPZk4dapszrTCil2n7/+9/z0EMPcffdd6e0sk9+8pNFvXHxw4zh1KlT1fr1jlrkCcviNVWp8DdWIZivmjgk79yIxphjkVic+mMRDjbmPls4fXgfnpx3jm0Zs5Rt2n2GQ4mlTn98c5dpKB0R4adXj/PNkibjbKYxLJAmoLVjhxrbqPM4O70V+x6/ldZ9H3qi0olI2fjx41sXLVrE17/+dVatWhVyO8QsFMFypg6AVXhm/bYXqcTMtKcNlbVcu+DNnLUNN24TduGnNfcO7f7MyhQ7RLeebFPmbais5fO/X+W5RlxoNO1s1qxZPPbYY0XXysAnmpmIHASagEPFbkuS/vijLf2BQyV9hpwWKu/W266giscikQM737ErkytS1rV7qEvFUCktQ0rKuoqEwm7Pjbc210Vr935kV8bp/uLHGg5KebfeEgqX5nqNApL3syJlXbuHu/cZLOHSLkrFogAi4RJQccIl5dn0tx3R+gPEmus9M7Yl3Sq+Dvyq2FoZ+ESYAYjIeqVU4bLcZoFf2uKXdkDQFiuCtviHIKFJQEBApyAQZgEBAZ0CPwmzBc5F2g2/tMUv7YCgLVYEbfEJvrGZBQQEBOSDnzSzgICAgJwJhFlAQECnIBBmAQEBnYJAmAUEBHQKAmEWEBDQKfj/b4Gf4nUlRHwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "lim = np.max(np.abs(S1))/1e+20\n",
    "fig,ax = plt.subplots()\n",
    "cax = ax.matshow(np.abs(S1.real),cmap=plt.cm.gray_r,vmin=0,vmax=lim)\n",
    "fig.colorbar(cax,extend='both')\n",
    "ax.plot(*Im[np.argwhere(np.abs(S1[np.ix_(Im,Im)].real)>1e-12)].T,'.',label='acoustic')\n",
    "ax.plot(*Ia[np.argwhere(np.abs(S1[np.ix_(Ia,Ia)].real)>1e-12)].T,'.',label='mechanic')\n",
    "ax.plot(*Ip[np.argwhere(np.abs(S1[np.ix_(Ip,Ip)].real)>1e-12)].T,'.',label='electrostatic')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, we can plot a vector of degrees of freedom marking the different regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD+CAYAAACZd9ZDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVoUlEQVR4nO3df6zddX3H8eebcqEttAIC5dLWFWapgy5CaIpbE4MylSmRuQwCy5DFzroEgmQuo+AyjI6sZoqaSIzXQYRMKFUhNAZ1QEYIC78KGrGtUwLILq0UpPIjAtr63h/ni9zAveece3u/3/O59/t8JM2953s/53u+3DSvvvicz/l8IzORJDVjv0FfgCS1iaErSQ0ydCWpQYauJDXI0JWkBhm6ktQgQ3eSIuL0iPjfiHgkItYP+npmo4i4JiJ2RcSPB30ts1VELI2I/46I7RGxNSI+PuhraotwnW7/ImIO8FPgPcAo8ABwbmZuG+iFzTIR8U7gReC6zFw56OuZjSJiGBjOzIciYgHwIPAX/l2un013clYDj2Tmo5n5G2AjcOaAr2nWycy7gGcHfR2zWWbuzMyHqu9fALYDiwd7Ve1g6E7OYuD/xjwexb+omuEiYhlwEnDfgC+lFQzdyYlxjjk/oxkrIg4Gvg1cnJnPD/p62sDQnZxRYOmYx0uAHQO6FmmfRMQQncD9RmbeNOjraQtDd3IeAJZHxDERcQBwDrB5wNckTVpEBHA1sD0zrxz09bSJoTsJmbkHuBD4Pp03HjZl5tbBXtXsExE3APcAKyJiNCLWDvqaZqE1wHnAuyPih9Wf9w/6otrAJWOS1CCbriQ1yNCVpAYZupLUIENXkhpk6E5RRKwb9DXMdv6Om+Hved9FxOMR8XC1CmRLt7G1hW4LduPyL2r9/B03w9/z9HhXZp6Ymau6DaoldKvduK4C/hw4Hjg3Io6v47UkaSapZZ1uRPwJ8KnMfF/1+FKAzPy38cbPW3hoLjhyZu0b89LzzzJv4WGDvoy+HTl30Fcwebt37+bQQw8d9GVMyot7Z96M3Yu/epaDD5k5f5d//pOHn8nMI/blHAuXr849v36ur7Ev7fjpVuDlMYdGMnNk7JiIeAzYTWcvlq++/udj7T/5y+3LeLtxnTJ2QDWPtA7g4COGOfvfv1nTpQjgohPmDPoSpGmxYsWKn+/rOfb8+jne9vdf6WvsD/7ltJd7TRkAazJzR0QcCdwWET+ptih9g7r+We65G1dmjmTmqsxcNZMaoyS9XmbuqL7uAm6ms/f2uOoKXXfjktQKEXFQdfcNIuIg4L3AhLeaqmt64fe7cQFP0tmN669rei1JGqRFwM2djdvYH7g+M7830eBaQjcz90TEq7txzQGucTcuSbNRZj4KvL3f8XU1XTLzVuDWus4vleh/npuBy0TUqJm3vkWSZrDamq7URmve9HLvQWo1m64kNcjQlaQGGbqS1CDndCW1zsHzhvjTlcN9jf3BNL+2TVeSGmToSlKDDF1JapChK0kNKuKNtOFnt/PJG07pPVBT9sK/dr1tk6SGFBG60myx4J977XWttnN6QZIaZNOVppHTOA345opBX8E+MXSlaeT0gnoxdKVpZNNtgE1XkmaW+UP7cfKSgwby2oauNI2cXlAvrl6QpAYZupLUIKcXpGnkG2kNmOFvpNl0JalBRTTdnx3yJj7wwdMGfRmz2sZBX4AkoJDQXXjAAZx2zLGDvgxJql0Rofv8b37DHY89OujLmNU+esSgr0AS7GPoRsTjwAvAXmBPZq6KiMOAG4FlwOPA2Zm5u9t5bLqS2mI6mu67MvOZMY/XA3dk5oaIWF89vqTbCebtP5c/PmL5NFyKJPU2b/85vO3QhQN57TqmF84ETq2+vxa4kx6h+9Kel3n46Z/VcCl61TuOeOegL0ES+x66CfxXRCTw1cwcARZl5k6AzNwZEUf2OolNV1Jb7GvorsnMHVWw3hYRP+n3iRGxDlgHsGDRQptuzWy6Uhn2KXQzc0f1dVdE3AysBp6KiOGq5Q4DuyZ47ggwAnDMymVp05XUBlMO3Yg4CNgvM1+ovn8v8GlgM3A+sKH6ekuvczmnWz+brlSGfWm6i4CbI+LV81yfmd+LiAeATRGxFngCOKvXiZzTldQWUw7dzHwUePs4x38JTOozvU/9+nm+eP/3p3op6sPG1TZdqQRFfCJt0fyFXLz6fYO+DEmakoiYA2wBnszMM7qNLSJ0bbr1s+lKtfo4sB3o+YmLIkLXpitppoqIJcAHgCuAf+g1vojQtenWz6YrveaAJ3/Essve0u/wwyNi7O70I9WS11d9EfgnYEE/JysidG26kgr2TGaOe8fRiDgD2JWZD0bEqf2crIjQtenWz6Yr1WIN8MGIeD8wF1gYEf+ZmX8z0ROKCF23dpQ0E2XmpcClAFXT/cdugQuFhK6bmNfPTcylMhQRujZdSTNdZt5JZyvbrooIXZtu/Wy6UhmKCF2brqS2KCJ0bbr1s+lKZSgidG26ktqiiNC16dbPpiuVoYjQtelKatLQQfNZ9I437Ew7vm/ePa2vXUTouom5pLYoInS9XU/9vF2PVIYiQtemK6ktighdm279bLpSGYoIXZuupLYoInRtuvWz6UplKCJ0bbqS2qKI0HUT8/q5iblUhiJC19v1SGqLIkLXpls/m65UhiJC16YrqVEHzoNj+vwYMA1/DDgirgFevePlyurYYcCNwDLgceDszNxd/exSYC2wF7goM3tWWJtu/Wy6Uhn6abpfB74MXDfm2HrgjszcEBHrq8eXRMTxwDnACcDRwO0RcVxm7u32AjZdSW3RM3Qz866IWPa6w2cCp1bfX0vnvkCXVMc3ZuYrwGMR8QiwGrin22vYdOtn05XKMNU53UWZuRMgM3dGxJHV8cXAvWPGjVbH3iAi1gHrABYsWujWjpJaYbrfSItxjuV4AzNzBBgBmH/sYekm5vVyE3OpDFMN3aciYrhqucPArur4KLB0zLglwI5eJ3MTc0ltMdXQ3QycD2yovt4y5vj1EXElnTfSlgP39zqZt+upn01XKkM/S8ZuoPOm2eERMQpcTidsN0XEWuAJ4CyAzNwaEZuAbcAe4IJeKxfApiupPfpZvXDuBD86bYLxVwBXTOYibLr1s+lKZSjiE2k2XUltUUTourWjpEYdeBAce0qfg6+a1pcuInTdxLx+bmIulaGI0LXpSmqLIkLXpls/m65UhiJC16YrqS2KCF2bbv1sulIZighdm66ktigidN3asX5u7SjVIyLmAncBB9LJ1G9l5uUTjS8idN3EXNIM9grw7sx8MSKGgLsj4ruZee94g4sIXZtu/Wy6Uj0yM4EXq4dD1Z9xt7SFQkLXpiupYIdHxJYxj0eq/cB/LyLmAA8CbwWuysz7JjpZEaFr062fTVd6zbaXnmPVw5v7Hf5MZq7qNqDaTfHEiDgEuDkiVmbmj8cbW0To2nQlzQaZ+auIuBM4HSg3dG269bPpSvWIiCOA31aBOw/4M+CzE40vInTd2lHSDDYMXFvN6+4HbMrM70w0uIjQdRPz+rmJuVSPzPwRcFK/44sIXZuupLYoInRtuvWz6UplKCJ0bbqS2qKI0LXp1s+mK5WhiNC16UpqiyJC160dJbVFEaHrJub1cxNz6TWT+b/rH0zzaxcRujZdSW1RROjadOtn05XK0DN0I+Ia4AxgV2aurI59Cvgo8HQ17LLMvLX62aXAWmAvcFFm9txUwaYrqS36abpfB74MXPe641/IzM+NPRARxwPnACcARwO3R8Rx1bZnE7Lp1s+mK5WhZ+hm5l0RsazP850JbMzMV4DHIuIRYDVwT7cn2XQltcW+zOleGBEfBrYAn8jM3cBiYOx9gUarY28QEeuAdQBDb57v1o41c2tHqQxTDd2vAJ+hcx+gzwCfBz4CxDhjx71XUHW7ixGAY1YuSzcxl9QGUwrdzHzq1e8j4mvAq3tHjgJLxwxdAuzodT43Ma+fTVcqw5RCNyKGM3Nn9fBDvHZbis3A9RFxJZ030pYD9/c6n7frkdQW/SwZuwE4lc4dMUeBy4FTI+JEOlMHjwMfA8jMrRGxCdgG7AEu6LVyAWy6TbDpSmXoZ/XCueMcvrrL+CuAKyZzETZdSU0a5IqpIj6RZtOtn01XKkMRoWvTldQWRYSuTbd+Nl2pDEWErpuYS2qLIkLX2/XUz9v1SGUoInRtupLaoojQtenWz6bbjHPu/+SgL0GFKyJ0bbqS2qKI0LXp1s+m24yNqyf1uSBNwQq+NehL2CdFhK5NV1KT5u5/IMcd0uJPpLmJuaS2KCJ0vV1P/bxdTzN8I029FBG6Nl1JbVFE6Np062fTbYZvpNXPN9KmgU1XUlsUEbo23frZdJvhnK56KSJ0bbqS2qKI0HVrx/q5tWMznNOtX2lzuhGxFLgOOAr4HTCSmV+aaHwRoesm5pJmsD3AJzLzoYhYADwYEbdl5rbxBhcRujbd+tl0m+GcbvtUd0bfWX3/QkRsBxbTuUHvGxQRujZdSU3avWMvN336uX6HHx4RW8Y8HsnMkfEGRsQy4CTgvolOVkTo2nTrZ9NthnO69RvAnO4zmbmq16CIOBj4NnBxZj4/0bgiQtemK2kmi4ghOoH7jcy8qdvYIkLXpls/m24znNNtn4gI4Gpge2Ze2Wt8EaHr1o6SZrA1wHnAwxHxw+rYZZl563iDiwhdNzGvn5uYN8M53fqVtk43M+8Got/xPUN3ooW/EXEYcCOwDHgcODszd1fPuRRYC+wFLsrMrnMHNl1JbdFP0x134S/wt8AdmbkhItYD64FLIuJ44BzgBOBo4PaIOC4z9070Ajbd+tl0m+GcrnrpGbpdFv6eCZxaDbsWuBO4pDq+MTNfAR6LiEeA1cA9E72GTVdSW0xqTvd1C38XVYFMZu6MiCOrYYuBe8c8bbQ69vpzrQPWAQy9eb5Nt2Y23WY4p1u/0uZ0J6vv0H39wt/OKonxh45zLN9woPOJjhGAo952VNp0JbVBX6E7wcLfpyJiuGq5w8Cu6vgosHTM05cAO7qd360dJTVp7vwh/ujtiwby2v2sXpho4e9m4HxgQ/X1ljHHr4+IK+m8kbYcuL/ba7iJef3cxLwZvpGmXvppuuMu/KUTtpsiYi3wBHAWQGZujYhNdHbY2QNc0G3lAth0JbVHP6sXui38PW2C51wB9P2Ogk23fjbdZvhGWv1a80ZanWy6ktqiiNC16dbPptsM53TVSxGha9OV1BZFhK5bO9bPrR2b4Zxu/ZzTnQZuYi6pLYoIXZtu/Wy6zXBOV70UEbo2XUltUUTo2nTrZ9NthnO69ZuOOd0DD9yft/zh4dNwNZNXROjadCW1RRGha9Otn023Gc7pqpciQtdNzCW1RRGh6+166ucm5s1wTrd+rtOdBjZdSW1RROjadOtn022Gc7rqpYjQtelKaosiQtemWz+bbjOc062fc7rTwKYrqS2KCF2bbv1sus1wTle9FBG6Nl1JTTrggDm8ZembBvLaRYSum5hLaosiQtfb9dTP2/U0wzfS6ucbadPApiupLYoIXZtu/Wy6zfCNNPVSROjadCW1RRGha9Otn023Gc7p1s853Wlg05U0U0XENcAZwK7MXNlrfM/QjYilwHXAUcDvgJHM/FJEfAr4KPB0NfSyzLy1es6lwFpgL3BRZnbdodxNzOvnJubNcE63lb4OfJlOTvbUT9PdA3wiMx+KiAXAgxFxW/WzL2Tm58YOjojjgXOAE4Cjgdsj4rjM3DvRC3i7HkkzVWbeFRHL+h3fM3Qzcyews/r+hYjYDizu8pQzgY2Z+QrwWEQ8AqwG7pnoCTbd+tl0m+Gcbv0GMKd7eERsGfN4JDNHpnqySc3pVml+EnAfsAa4MCI+DGyh04Z30wnke8c8bZRxQjoi1gHrAN589GE2XUmN+cUvnuSzn13f7/BnMnPVdL1236EbEQcD3wYuzsznI+IrwGeArL5+HvgIEOM8Pd9woPMvxQjA/GMPS5tuvWy6zXBOV730FboRMUQncL+RmTcBZOZTY37+NeA71cNRYOmYpy8BdnQ7v3O6ktqin9ULAVwNbM/MK8ccH67mewE+BPy4+n4zcH1EXEnnjbTlwP3dXsM53frZdJvhnG79SlunGxE3AKfSmfsdBS7PzKsnGt9P010DnAc8HBE/rI5dBpwbESfSmTp4HPgYQGZujYhNwDY6Kx8u6LZyAdzaUdLMlZnnTmZ8P6sX7mb8edpbuzznCqDvf/LdxLx+bmLeDOd01UsRn0iz6UpqiyJC16ZbP5tuM5zTrV9pc7qTVUTo2nQltUURoWvTrZ9NtxnO6aqXIkLXpiupLYoIXbd2lNSk+fPnc/LJJ/c19uqrJ1xyOyVFhK6bmNfPTcylMhQRujZdSW1RROjadOtn022GS8bq55KxaWDTldQWRYSuTbd+Nl2pDEWErk1XUlsUEbpu7Vg/t3Zshh+OUC9FhK6bmEtqiyJC16ZbP5uuVIYiQtemK6ktighdm279bLrNcJ1u/aZjne7cuXNZsWLFNFzN5BURujZdSW1RROjadOtn05XKUETourWjpLYoInTdxLx+bmLeDNfpqpciQtemK6ktighdm279bLpSGYoIXZuupLYoInRtuvWz6UplKCJ0bbqS2qJn6EbEXOAu4MBq/Lcy8/KIOAy4EVgGPA6cnZm7q+dcCqwF9gIXZWbXRbhu7SipSUNDQyxZsmQgr91P030FeHdmvhgRQ8DdEfFd4C+BOzJzQ0SsB9YDl0TE8cA5wAnA0cDtEXFcZu6d6AXcxLx+bmIulaFn6GZmAi9WD4eqPwmcCZxaHb8WuBO4pDq+MTNfAR6LiEeA1cA9E72GTVdSW/Q1pxsRc4AHgbcCV2XmfRGxKDN3AmTmzog4shq+GLh3zNNHq2OvP+c6YB3AgkULbbo1s+k2ww1v6teKG1NWUwMnRsQhwM0RsbLL8BjvFOOccwQYAThm5bK06Upqg0mtXsjMX0XEncDpwFMRMVy13GFgVzVsFFg65mlLgB3dzuucbv1sulIZ+lm9cATw2ypw5wF/BnwW2AycD2yovt5SPWUzcH1EXEnnjbTlwP3dXsM5XUlt0U/THQaureZ19wM2ZeZ3IuIeYFNErAWeAM4CyMytEbEJ2AbsAS7otnIB3NqxCW7tKJWhn9ULPwJOGuf4L4HTJnjOFUDf7yi4ibmkmSoiTge+BMwB/iMzN3QbX8Qn0my69bPpStOvmgG4CngPnfezHoiIzZm5baLnFBG6Nl1JM9Rq4JHMfBQgIjbS+azChKEbnc8+DFZEPA38fNDXMUmHA88M+iJmOX/HzZhpv+c/yMx92sIpIr5H57+7H3OBl8c8HqmWvBIRfwWcnpl/Vz0+DzglMy+c6GRFNN19/QUOQkRsycxVg76O2czfcTPa+HvOzNOn6VR9fS5hrP2m6YUlqY0m/bkEQ1eSpu4BYHlEHBMRB9DZ7GtztycUMb0wQ40M+gJawN9xM/w9T1Fm7omIC4Hv01kydk1mbu32nCLeSJOktnB6QZIaZOhKUoMMXUlqkKErSQ0ydCWpQYauJDXI0JWkBv0/JUV8i1ZpQCMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "D = np.zeros([S1.shape[0],3])\n",
    "\n",
    "# in column 1\n",
    "D[Ia,0] = 4.6 # acoustics in blue\n",
    "D[Ip,0] = 3.8 # electrostatics in orange\n",
    "D[Ix,0] = 2.8 # mechanics in green\n",
    "D[Iy,0] = 2.6 # bit lighter for y\n",
    "D[Iz,0] = 2.4 # and z\n",
    "\n",
    "# in column 3\n",
    "D[Iai,2] = 4.4 # interface (acou)\n",
    "D[Ipe,2] = 3.8 # piezo electrode\n",
    "D[Imi,2] = 2.8 # interface (mech)\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "cax = ax.matshow(D,cmap=plt.cm.tab20c_r,vmin=0,vmax=5)\n",
    "fig.colorbar(cax)\n",
    "ax.set_aspect('auto')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also extract certain sub-matrix-blocks and plot them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADqCAYAAACcPZ9GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAR6ElEQVR4nO3df4zk9V3H8dd7Zu+Wu9vjFjigBXocCqEtSGNTDQkNVWITwR/VxD80WvuHiWnaakFJrWlibPjHn4CNv+KPP2ps0hhFJU2N0dpiQoOxomkktOXEFhuIcNJd7oDb253v2z/mOzuzu5/393bmZua9M/t8JBtmP/P5fr+fz+zea4ad93w+5u4CAExfK3sAALBfEcAAkIQABoAkBDAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIsDNP5+PIxv/6aq3becclSeExHVr4j+ABeO+i+l1VBe9Oz20Yw/4WG+Xfm6DEr+cZzz+n06dNTn83x48f9xIkT074s9onnGn6vhwrg66+5Sk986hM72qsbbw+PeVWLxfZO8BHooweSX5R7FKeSrDy2tSAZFxuScWWtfJ3lxXj+Z9fLxyxlP2Zjcscdd6Rc98SJE3r88cdTro351/R7PR//cgFgBhHAAJCEAAaAJEP9DViXLBX/3ts69UR4SOfb7yy2t21vvnO0VsXjWmwH7a1oSc/4XAsjPPVtNPx5GsDs4RUwACQhgAEgCQEMAEkIYABIQgADQBICGACSDFWG1pEVP1oclZpJ0vmH7iu2X/mhXy22VwuXhed6LVhA4XDDAgpngo/vRh95fj1apEHxR4vPbpT7Lx0IT6Wl9vArSEyjcm+Uj0iPovRzida6AOYVr4ABIAkBDABJCGAASEIAA0ASAhgAkgy3GI+XF1JvWlgnqnZ46XfL7Vd85HfCc0WLuDctejPsoj/TWiRndaN8oWMH42OOLYyy98Zw2lN6Su4UpkIRBPYbXgEDQBICGACSEMAAkIQABoAkBDAAJBmqCqJtw28bH63tEFU7tL/yhfj6N76r2G4b58JjOh6VFZSrI0bZ4n2UdRKOHRzhua81XNFKIy9XVJSqE0a11rC4Q6naZG9uUgVMDq+AASAJAQwASQhgAEhCAANAEgIYAJIQwACQZIx1TWXRNkLRwjpRqZkkvf7gveVz3fdQeMyx9dVie3Ug3vpoGkbZ+mfY7ZUar3++3L58YHwL/qxXcRlaqUKNxXiw3/AKGACSEMAAkIQABoAkBDAAJCGAASDJxKsgDi9ES6yU2xsX1gmqHZb/65/DY84EVRWHwyOmY5QFfEapdogcCn4uZzvlayyNcOmlht+uqAoD2E94BQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQTL0OLFpBpW7kMKt7DLV5YJyo1k+IFfA4He9JNyygL64yygE8kWiinYRu3oTWVmpUW/QkrFoE5xStgAEhCAANAEgIYAJIQwACQhAAGgCQTr4IYfgGZ+K3waBuhpoV1omqH9le+UGzv3HxneC7rlN/WP2eL5XM1VBR0gp1/ouoISVoY49Pl0hgX9ok0V2fw3A/wrwAAkhDAAJCEAAaAJAQwACQhgAEgyfBVEL7zXfq1Kq5ceH2jXAoQvaM/jXfnpbja4eUH7w+POXLvbxfbD62tlM/VPhZf34PHpaEKZDVYC2JajxmA8eJfLgAkIYABIAkBDABJCGAASEIAA0ASAhgAkgxfhmY7M3uxHXdfbO/NfWaihXWiUjNJOvzMY8X2qKStvRFfPyo3awVbNUnSEcrNgLnCv2gASEIAA0ASAhgAkhDAAJCEAAaAJENVQVSS1gr77Cy24r13zgaVABvBzjvN29iMT7SNULSwjjT8Aj4Hf/63wnNFC+s0VTpEC/gAmE28AgaAJAQwACQhgAEgCQEMAEkIYABIMlQVREvR2g7x+gVLB4Yc0ZQUijkkNW8jFK3tEFU7rD10X3iua+8vrznROrcaHvOthfLYXgu2fTq8MPw6HCtBdUa74an66AhrVJxZ33md6GcCzCteAQNAEgIYAJIQwACQhAAGgCQEMAAkIYABIMlQZWgbXi5TWmiI8aV2uaxpdaN80LGD8cmiEqmmBXxK5U6S1AkWA2pa8CbaRihaWCcqNZOke5e+o9j+8ZWnwmOuWH222H760m8rto9ShhbNv91QajiKVwqPGWVo2G94BQwASQhgAEhCAANAEgIYAJIQwACQZKgqiAUbZcugqNphyNNotO2KooViouqIqNJBklpWvi/aRqhpYZ2o2mHj4XgBn42gqqLdCQ8Z2sHiYktSO5j7qEqPWWu8lwD2PF4BA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgyVBlaB2XzhbKtzaChW0kKapeOrYQHNSKhxSVjjXtSRYt4BMtIBQtrCPF5WbRAjbRHm5SvLBOVGomxQv4PLAaL+AzrNIiOVI8d2m8i/4A+wmvgAEgCQEMAEkIYABIQgADQBICGACSDFUF0TZpqeHd8OGMb2GdJsMu4DO++UmvbcTv9EfbCDUtrBNVO5x/qLyAz2u/8HB4rvPB/j/XHi73b1pYaHW9XO3RVBuxwMo7AK+AASALAQwASQhgAEhCAANAEgIYAJIMVQWB4TStkTDK+gmRqNrh6KnHwmM6N99ZbF85X35Obh+I17W41M4X28/ZYniMqp1VGNRFYL/hFTAAJCGAASAJAQwASQhgAEhCAANAEgIYAJLsvTI0j7cEWilXO+lQQ0nXeqHcSRrvojvZooV1olIzSXr5wfuL7f7B3yy2LzQUiUXlZofWVsJjTldHd7QF0wDm1vykEADMGAIYAJIQwACQhAAGgCQEMAAkmXgVxMpauaqhHUR/Jy6C0PKB8p1nO/HzSOY769HcJanj5YEdbMfVBq8E54u2EYoW1pHiaocHLr+12P7xlfJ2SJKKC+tI5UqHnutf//qOtoPVWnwNYA7xChgAkhDAAJCEAAaAJAQwACQhgAEgycSrIJYXx5nx5XMt7dGnkajSQ5LawdoKbYurII4E61e0zq2Wz9WwjVC0tkNU7bDx8H3huQ4EWyJF45WkjcMnC4Nq2MIImEN7NLoAYP4RwACQhAAGgCQEMAAkIYABIAkBDABJ9t6WRHPk6Ji3PTocbL20ul4uN7vUgj2cFG8jFC2sE5WaSdLRU48V25u2RCotFLTBlkTYZ3gFDABJCGAASEIAA0ASAhgAkhDAAJBkz1VBrDXsIbQevEO/1DCLlaAQYLyLBOWKlu8JKx0kHVpbKbZH2wg1LawTVTu8/OD94THRlkjAfjI/KQQAM4YABoAkBDAAJCGAASAJAQwASQhgAEgy8TK0M+tVsb1TbtZCw1NCVKEWlZpJ0vKB4EJTeO6J5i5Jr6yV72ss9/LyA7DQCgrRgrI9KS43u/71rxfbi3u41UoL60jNpWYPXH7rjrbn9ULYH5hHvAIGgCQEMAAkIYABIAkBDABJCGAASGIevLNe7Gz2kqRvTG442Oeud/crp31Rfq8xYeHv9VABDAAYH/4EAQBJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkIQABoAkBDAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAuTvoCZ3SLpY5JuluSTvh4ATMBpSb8h6fPuPrYcszGea+uJzW6x5ZP/qc55ta5+m1pLb5C12lu/2m1JUmtLe2vztqTNfr221mD7ti9JW/p2+1vd32RmarWse7tlatnAfa3efZLV7dv72sC5WvX5rNXvV8978/t2/SVp83a7ZVrY0d7acr8ktW1n/9K5wnbrz69t3TH22nrfS+p/X19vyzGbbd3xbPZv9Y4b6F+3mUkt1f+13u/C1ra6eaB//bht6d9vMw2ea7Cte7wkmbvklcwrySup93u9pa3+qqr6mIG+vX71Mb0vG+ivHdeo26uOvKqkqtO93en0j6k68qr7vdf3d9urul+nf3zdv3i++pje+byq5J3uMd6pNr/3TqWq7j/Y5lWlarD/9mOqwfs6W+7rDmuwv/fPVbmqjg+0e/+Yyjf7eqd7u3+urcd1270ep9f9u7e743JV7uq4q+NSp/759m4P/re5vd9WKb5v+zVWtK4ntapXtaG365g+oxdb4wjisb8C7gWvHb1GreNvUWvpDeO+BABM1bIO6C4d3wzia7RY/ZBdfdFBPLYAJngBzLtxB/FFBzDBC2C/GVcQjxzABC+A/e5ig3joACZ4AWCrUYN41wFsZjfa8slnCF4AKCsFsZm9090fL/Uf5oMYt/mrp9W67EbZkavGNFwAmD/HtKCbdUTnVEnSXVG/Xb8CdvdHzOwqP/fyi52vPaX2lbfILrtBZnyYDgAkyeV6Xuf0pFZ1RAt6WetvdvevRv2H+huwu78kyczsSoIYALq2B+/zWnuz+7kweHtGqoIgiAFg9ODtuag6YIIYwH50scHbM5ZPwhHEAPaDcQVvz1jXgiCIAcyjcQdvz8RWQ5MkM7uyddWtL1ZnXpAduqK7yphZN4zr23W/egWylmTq31/f1+tram2uVNZvM8kK7YP3S/V56+9M9X0auEZ9SdOOY3qn3RzTwDH9S/bOtvV8rW3H9FcEsy3H9voOjql3id6qX5urtNnW84XtGry2NtsGprL5mJi2zqPf1h+HtvXr9hmcR7+PBq6hut+O+3uPcbF/fW7rj0Xb+tnAweYuyesFTwd+p3vt8oEV0rw+fntb77j+eWywvXfuLX1Vr4zW+6/LfeC8vf6l9sol9VZjGzyXyweO2zLG3nk2j9fA993bXurvPnBpl1fdY7zyeoj9Y9y3jXnz4anbqsH+2nr+gXN5Pa+d1+iNsz7n5mJz3j/GteX66j3EclWDD5d327adavMh9fq8O9v7x/TOVQ3eN9guaVXrWtKCntGrjVUNw5poAG9exOwKSScnfJkflfTXE77GXsb8mf9+nr802cfgRXf/n3GfdCoBPA1m9iV3f0f2OLIwf+a/n+cvzeZjwB9mASAJAQwASeYpgP8oewDJmP/+tt/nL83gYzA3fwMGgFkzT6+AAWCmEMAAkGRmA9jMLjezfzCzZ+r/XtbQt21m/25mn5nmGCdpN/M3szeZ2efN7Gkze8rMPpwx1nEys+83s6+a2Skz+2jhfjOzT9T3f9nM3p4xzknZxfx/sp73l83si2b2toxxTsqF5j/Q77vMrGNmPzbN8Q1rZgNY0kclfc7db5L0ufr7yIclPT2VUU3Pbua/IekX3f0tkm6X9EEze+sUxzhWZtaW9HuS7pb0Vkk/UZjP3ZJuqr9+VtIfTHWQE7TL+f+3pHe5+22SHtAMvjEV2eX8e/1+XdLfT3eEw5vlAH6PpE/Wtz8p6UdKnczsOkk/IOlPpjOsqbng/N39BXd/sr59Rt0noWunNcAJ+G5Jp9z9WXc/L+nT6j4Og94j6c+86wlJy2b2xmkPdEIuOH93/6K7f6v+9glJ1015jJO0m5+/JP2cpL+S9OI0BzeKWQ7gq939BakbNJKifZIelvQR9T8mPi92O39JkpmdlPSdkv5l8kObmGslDX4c9Jva+YSymz6zati5/Yykv5voiKbrgvM3s2vV/UjyH05xXCMb62po42Zm/yiptPvnx3Z5/A+q+xnufzOz7xnj0KbiYuc/cJ4ldV8R3Ovur4xjbEms0La9jnI3fWbVrudmZt+rbgC/c6Ijmq7dzP9hSb/k7p3BBbL2qj0dwO7+fdF9Zva/ZvZGd3+h/l/M0v9u3CHph83sHkmXSLrUzP7c3X9qQkMeqzHMX2Z2QN3w/ZS7PzKhoU7LNyW9aeD76yQ9P0KfWbWruZnZber+ye1ud/+/KY1tGnYz/3dI+nQdvscl3WNmG+7+N1MZ4ZBm+U8Qj0p6X337fZL+dnsHd/9ld7/O3U9K+nFJ/zQr4bsLF5y/dX8L/1TS0+7+4BTHNin/KukmM7vBzA6q+zN9dFufRyX9dF0Ncbuk1d6faubABedvZickPSLpve7+tYQxTtIF5+/uN7j7yfrf/F9K+sBeDV9ptgP41yS928yekfTu+nuZ2TVm9tnUkU3HbuZ/h6T3SrrLzP6j/ronZ7gXz903JH1I3Xe3n5b0F+7+lJm938zeX3f7rKRnJZ2S9MeSPpAy2AnY5fx/RdIVkn6//nl/KWm4Y7fL+c8UPooMAElm+RUwAMw0AhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEn+HzOZsGEAQhroAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD3CAYAAAAjdY4DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAASFklEQVR4nO3db4xc11nH8d8zM+vZjb3+010nzl8nobXTJEpLoEFqyj9TmhCgpUgVIF4g0Td9AW1fFBU1EhVEBRUKQkICUakRIPpHoLRSVSAKIaWUQJq0VhIS0qQhNImbGnudrB3vxuuduQ8v5s7OzHr3nBnvrJ/NzPcjrTJzz7nnPjNe//ZmfM5Zc3cBAC68SnQBADCuCGAACEIAA0AQAhgAghDAABCEAAaAIAQwAAQhgAEgCAEMAEFqg3R+w8yMX3Hl/nXbJ9TID1JJX7Loo47sTw1v5gexarK5jxGUHkGS9/FqbPN/Bvaz1tE2vQopV8nzz7+guRMnLkwpXWZnZ33/VVdd6MtiTDz/wguam5tb8/t6oAC+4sr9+vK/fG3d9suKV7JjFDv2JtuXmvm4qFfTf0ftzKnsGD65M9m+2MjXcVEtU8fS6Xwd9R3ZPhvVx1uqzFs6FNY4k2x/+48f2vwi1rD/qqv04IMPhlwbo+/WW29dt42PIAAgCAEMAEEIYAAIQgADQBACGACCDDQLYkKN5EyHb7x6UXaMt2X+0X9y+dXsGF5Nz2CYr+RnFuzKtG9vLuTrqKWvM+dT2TFmsj02rtZMzz6QJK9NbnodS5X0NYoLNBkO2Cq4AwaAIAQwAAQhgAEgCAEMAEEIYAAIQgADQBACGACCEMAAEGSghRiq1JLbSeYWWUjS0cX0nsF76tPZMaZOH0+2767Vs2O40os5hrFN5MxkdsfgochtnVmv5hdZDKPS3Pab9Wr6263S187FwOjgDhgAghDAABCEAAaAIAQwAAQhgAEgCAEMAEEIYAAIMtA84ELpXxvfz2bquXm+jx9bzI5xy570PN/a/JHsGMv7rk+2D+NXued+Dbs0nI3Q65lC5pea2TGGMmc5M8/XXjuZPr/I1wmMEu6AASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCAEMAAEGWghRkXpSf9eTW9yLuU3U88tspCkh19JLxo4MHNddoxdmfbcIot+DGORRT9qZ9Mboc/W8n/MPoQt2ecaE8n26e2XJNuLymC/HwB4veMOGACCEMAAEIQABoAgBDAABCGAASAIAQwAQQhgAAgy2MRLb8rOnFq3eb6yIzvE7trGN1PPzfO96/5ns2N86o4DyfZhbKZ+tsgOoW3D+BG40Y3QJfn0xucsT29LzyX+7smzyfaz/eyCD4wQ7oABIAgBDABBCGAACEIAA0AQAhgAghDAABCEAAaAIAQwAAQZbCGGVeWT62+6ntvkXJJc6U3bl/ddnx0jd53cIgtJOrrYSLbv3JbfGH5HZnP5+sRUdgyvpxevLPWxOKGeWRDS2JFfZJHbjj21AKdtW+J7Q5IOTqdXpkxufE944HWFO2AACEIAA0AQAhgAghDAABCEAAaAIAQwAAQhgAEgyEDzgJuSFhvrz0vd3lzIjpGb99rPntxVS7f3s5l6bp7vk8cXs2Pc8ob0PN/aKy9kx8jNe65VMi+2D6eX8zvD78rtDF/btuE6Kgsn0h2a6bnZwKjhDhgAghDAABCEAAaAIAQwAAQhgAEgCAEMAEEIYAAIQgADQJCBFmJUJV1UW39hgNfSiyz6usbG1x3IMxuUS/nN1HOLLCTp4ZfTxV6752B2jJlMe20pvxF6apN8SdpdOZsfQ+n3bL7IL8TIbZT/2vTlyfaiOpG9BjBKuAMGgCAEMAAEIYABIAgBDABBCGAACEIAA0AQAhgAggw0D1heyJZOr9s85/m5szOT1WR7P5up5+b5ns3vP676xMY3U8/N8/2DB/4nO8an7jiQ7nAhNkKX1NyVnqM7Vdv4z+qjC8vJ9uWij934gRHCHTAABCGAASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCCDLcSwiry+/qbruc3F+9HPZuo52/r4sZJ6HZK0vO/67Bi515tdZCHp6GIj2b6nXs+O0WykFzDUd6YXWUjSRGaD+notX4cs/cbvn0p/u9W5HcCY4VseAIIQwAAQhAAGgCAEMAAEIYABIAgBDABBCGAACDLYPOARstRMz52tVSw7Rm3pVKZDfjP13Dzfx48tZse4ed/2ZPv8UjM7xmxmnm9t/kh2jMbstcl2e+1keoAiXycwSrgDBoAgBDAABCGAASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQ5IIvxFjMbR5e7WMBxNnT6Q7V/MuqD2Hjd5/cueExcpup5xZZSNLhowvJ9gMzU/lCPP2zOLfIQpJqc88l24upXZka0u8FMGq4AwaAIAQwAAQhgAEgCAEMAEEIYAAIQgADQBACGACCDDQP2CWl9jGvNc9kx6hX0/Nv+9s8PF12duNvSY0d6TpOLxfZMXZXzibbKwsnsmPUd16ebO/n/cjN873r/mezY/zxu65Ktvfznubm+TYeuTfZ7ov5awCjhDtgAAhCAANAEAIYAIIQwAAQhAAGgCAEMAAEIYABIAgBDABBBlqIYZJS+6V7H5ucVzPtM5O5HpJnRvHpjdexa1v+Z5MrfZ3mrvQiC0maOH082T5bq2fHyG2mnltkIUkPH08v+Dgwszc7RqNIb6g+fej9yXafvjt7DWCUcAcMAEEIYAAIQgADQBACGACCEMAAEIQABoAgBDAABBloHvAosTOn0h1q27JjzBfpPlO1/M+3emaeb23+SHaMxuy1yfZ+NlPPzfPtZ1P337/9QLL9uyfTG9ifTe32D4wg7oABIAgBDABBCGAACEIAA0AQAhgAghDAABCEAAaAIAQwAAS54AsxbOl0ukM1X9JcYyLZPr0tv6n7tsmd2T45uzY8giRL/wzMLbKQpNrcc8n2YipfaW4z9dwiC0n62L3PJNs/+KPXJNstsdk/MIq4AwaAIAQwAAQhgAEgCAEMAEEIYAAIQgADQBACGACCDDgP2GWNM+u2LlUmsyPUM/N8+9k8fHr7Jcn23MbfknRwuki2VxZOZMd4bfryZPvRheXsGPunNv5+5Ob5Nh65NzvG9KH3J9v7eU9z83xvvO3DyfbGMy9mrwGMEu6AASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCAEMAAEMff0Rtw9nc2OS3p+88rBmNvv7nsv9EX5vsYmW/f7eqAABgAMDx9BAEAQAhgAghDAABCEAAaAIAQwAAQhgAEgCAEMAEEIYAAIQgADQBACGACCEMAAEIQABoAgBDAABCGAASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCAEMAAEIYABIAgBDABBCGAACEIAA0AQAhgAghDAABCEAAaAIAQwAAQhgAEgSG2zL2BmN0i6U9JBSb7Z1wOATTAn6Q8lfdXdh5ZjNsSxegc2u8F2X/2EmmdVueQtquzYJ6tUe7+qVUlSped4ZeWxpJV+7WOV7uOrviT19G31t7K/ycxUqVjrccVUsa62SrtNsvL46r7WNValHM8qnX7l6155Xi2/JK08rlZMtXOOV3raJalq5/Zfa6x1j1vn9VWtVWP7WPu5pM7z8no956wca9Wz0r/SPq+rf3nMTKqo/K+1vxd6j5WHu/qX71tP/84xU/dY3cda50uSuUteyLyQvJDa39c9x8qvoijP6erb7lee0/6yrv465xrl8aIpLwqpaLYeN5udc4qmvGg997K9dbwo+zU755f91xyvPKc9nheFvNk6x5vFynNvFirK/t3HvChUdPdffU7R3dbsaWuV1d3fO2MVrqLpXce9c07hK3292XrcGav3vNZxL+v0sn/rcasuV+GupruaLjXLP9/24+7/po93jhVav231Nea1rMM6qQU1dLN26Ss6VhlGEA/9DrgdvDZ9mSqzb1Zlx75hXwIALqjdmtAhza4E8WWqFz9vl2w4iIcWwAQvgFE37CDecAATvADGzbCC+LwDmOAFMO42GsQDBzDBCwC9zjeI+w5gM3uj7b76OwQvAKxtrSA2s3e4+4Nr9R9kIcZNvjCnyp43yrZfPKRyAWD07FJNB7VdZ1RI0qH1+vV9B+zuXzSzi/3My8eazzyp6t4bZHuukRmL6QBAklyul3RGh3VS21XTy1q+zt2fXq//QJ8Bu/txSWZmewliAGhZHbwvaek69zPrBm/bec2CIIgB4PyDt21D84AJYgDjaKPB2zaUlXAEMYBxMKzgbRvqXhAEMYBRNOzgbdu03dAkycz2Vi6+8Vjx6vdlUzOtXcbMWmFcPi77lTuQVSRTp71sa/c1VVZ2KuscM8nWON7dLpXjls9MZZu6rlFe0nTOOe1hV2rqOqdzyfZoveNVVp3T2RHMes5t9+2uqX2J9q5fK7u0We946x5X97W1cqzrpay8J6be19E51qlDq/q1+nS/jk4fdV1DZb9z2tvv8Zr9y7GtU4tW9bOuk81dkpcbnnZ9T7ePy7t2SPPy/NXH2ud1xrHu4+2xe/qq3Bmt/V+Xe9e47f5rHS9cUns3tu6xXN51Xk+N7XFWzlfX89ZjX6u/e9elXV60zvHCyxI757ivqnnl7SmPFd391Tt+11hevq5zr9GusxxzZbM575zj6rm+2m+xXEX32+WtY6uGWnlLvRz33OOdc9pjFd1t3cclndSydqim72ghOathUJsawCsXMZuRdPWmX6h/75X0pegi1kFt54fazg+19eeYu7847EEvSABvNWb2TXf/4eg61kJt54fazg+1xeKDWQAIQgADQJBxDeBPRxeQQG3nh9rOD7UFGsvPgAFgKxjXO2AACEcAA0CQsQ1gM7vLzB43s0fN7D4zuyy6pjYz+yMz+3ZZ35fMbHd0TW1m9j4ze9LMCjPbElOEzOx2M3vazJ41s9+OrqfNzO42s2Nm9kR0LauZ2ZVm9lUze6r88/xQdE2SZGaTZvawmT1W1vW70TVtprH9DNjMdrr7qfLxByVd7+4fCC5LkmRm75L0gLs3zOyTkuTuHw0uS5JkZm9Wa3HQX0r6iLt/M7ieqqRnJP20pCOSHpH0K+7+35F1SZKZ/Zik05L+xt1vjK6nm5ldKulSdz9sZtOSviXpF6LfN2stxdzu7qfNbELSv0v6kLs/FFnXZhnbO+B2+Ja2q2dtaSx3v8/dG+XThyRdEVlPN3d/aphLMYfgFknPuvtz7n5W0hckvSe4JkmSu/+bpJej61iLu3/f3Q+Xj1+V9JSky2OrkrzldPl0ovzaMn83h21sA1iSzOwTZvaipF+V9DvR9azj1yX9U3QRW9jlkrqXiB7RFgiS1xMzu1rSD0r6RnApklr/V2Nmj0o6Jumf3X1L1LUZRjqAzex+M3tija/3SJK73+nuV0r6rKTf2Eq1lX3ulNQo69tStW0htsaxkb1jGjYz2yHpHkkfXvV/hWHcvenub1Xr//xuMbMt9fHNMA11O8qtxt3f2WfXz0n6B0kf38RyeuRqM7Nfk/Rzkn4q9WutN8MA79tWcETSlV3Pr5D0UlAtryvlZ6z3SPqsu38xup7V3H3ezP5V0u2Sttw/ZA7DSN8Bp5jZm7qevlvSt6NqWc3Mbpf0UUnvdvfF6Hq2uEckvcnMrjGzbZJ+WdKXg2va8sp/7PqMpKfc/U+i62kzs73tWT9mNiXpndpCfzeHbZxnQdwj6aBa/6L/vKQPuPv3YqtqMbNnJdUlnSgPPbSFZmi8V9KfSdoraV7So+5+W3BNd0j6U0lVSXe7+yci62kzs89L+glJs5L+T9LH3f0zoUWVzOwdkr4u6b/U2Ub3Y+7+j3FVSWZ2k6S/VuvPsiLp79z99yJr2kxjG8AAEG1sP4IAgGgEMAAEIYABIAgBDABBCGAAY2nYmyWZ2Se7Fi39Uj/nEMAAxtVfqbXIY8PM7Gcl3SzprZJ+RNJvmdnO3HkEMICxtNZmSWb2A2Z2r5l9y8y+bmbX9Tnc9ZK+5u4Nd1+Q9Jj6CHcCGAA6Pi3pN939hyR9RNKf93neY5J+xswuMrNZST+p3iXyaxrpvSAAoF/lxkRvl/T3rZXaklorUmVmvyhprRV533P329z9PjN7m6T/kHRc0n+qtZFW+pqshAMwrsqtOL/i7jeWn9k+7e6XDmHcz0n629zSbj6CAACt/JKG/zWz90mtDYvM7C39nFvuYTxTPr5J0k2S7suexx0wgHG01mZJkh6Q9BeSLlXrt3F8oZ/NgMxsUtLh8ukptTb3ejR7HgEMADH4CAIAghDAABCEAAaAIAQwAAQhgAEgCAEMAEEIYAAI8v/0/mhEFZyd2AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD3CAYAAAAjdY4DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAk8UlEQVR4nO2dXYwk13Xf/+dWTc/scLnLJZfkkit+iLFCSfSHICg2IwaGIDiKICBQDERAgACRAQGBHxIgiV8ECLGN6CGIgBhBghiIHwI7D3GAAHEiGLYVObDihLAC2wwpS9YHaYUSv7+5y+XsTE9VnTzUvdW3qm9VV/d0d3VV/X/A7HTfunXrVvfOv06de84pUVUQQgjZPqbrCRBCyFihABNCSEdQgAkhpCMowIQQ0hEUYEII6QgKMCGEdAQFmBBCOoICTAghHUEBJoSQjoiX6Xz58mV94OoVaDSB2DYFIJoCEuWv7b/VbQCQKhBJYOBZ7+CWDLxSjIEf/PCHeP3112v/h2yKy5cv6wP337/tw5KR8PQzz+Cdd975iqp+srptKQF+4P778T//9+PYv/4C0otXi/aTVLEfCaAZEhjE9k8otVnOkTSLr9q+cWB7UtO+MbIEMPnHMs2APQOIZoDMLgEKQLx+ZD089thjnRz3gfvvx+OPP97Jscnweeyxx/CnTzwxJ77ACoZlLEB68SpkepQLrGaITW73pp74ArngOtFtEt/Q9kQBZEk+nmbLTnN1PFGdGGuTS/ljkko/QghZhZXv7JO9Q8Qn16FiEGmCTAEjudW4DJIlEM1yK9MjFiCVGKkCCR0QhJABsrSyFW4FZMgOLiC68RpuapSLb6qY+CNqVli4zlr2UeQiCzFzIptqfoxIgL3To2WnSQghO89SAqwA4pPrOEkVKQwkS/Duucs4/90/gEmOEZmyHyGFgaRTALmYTnX+Vj5CBmQJ9t9+rrTN7/9SMlnytAghZPdZSoAFgB5cwH4kMNZFcBAJsr/8UZh334DA+m4tkQCIrHiKKVvHlhQGMDHSi/fMbZuYfLwrh9v1t6aaX2xC7hR3B5CwjDIh5IwsbQErAGg28/mmitPoANMLV3Hwyp/PRSw4wQrpVWrHOEkVNzUqH8sufG0zAsLN1XenVF0nbr50ixBCzsrSPmCTHCOBQSTAaQbsR4JYcoE9vfJBxN//etE3VSDSBIAN26qOJUCm+RjnKmEQhRjeeG3ZKa6ME1d3XtMst+J9/7Sbr04OtzYvQsgwWd4FER8UVqnvUnBtyUOPQqZHQJbACAAxUACJxHO37YJc4PxmJ7xuvOz8nUXbphHk4urOy/2uhtYRQsg62Eh8l8YTRG+/gFQBOblRCOhx4jlVqxaxW6wTQJLjUr9tCTAhhGyTpQU4UZQFUbPS+1QBNTFOLj2A/esv4Hp0HoLcirxVprN+MltYEwAn2MuHA5BEB8W2oyzCwcm1ZadJCCE7z9KLcHsn18ttYhBlubA6n69ohr2T63j78B7c9hd/WPT9v2/MlDpCOcTg8MaLAPK0X+c3BoDzp9fxot6yzDQJIaQXrBSG5vtBBShCzSIBYGKoGJzuX8BhbDB9+GOYXH8B0AwfuvtcsV81JvjkQl5bIoVB4lnH04OLWw9Dm8NZ+b7bJEvymhC+u4QQQpZgaQu46pP1w7Sq210tiPTCPdh75Tulsfa8I0+z2UJXJZdjq2Foic7ie905nKQKiClddKYZABMjU0AjJokQQlZjaQv4RloWIz9jzS9Sc5xkkOQYT75yExCD0ysfxJv/8p8U+/kiPjHAy0dJMUZVdP3Qtk1TXSjcjwRPvnITqebnDszmG8ruI4SQtix9b39+ryI4YjAJWK3n9wwUB/jQ3bP223/hVyDTozxKwqsmpphlu7nkDN+KPnnvoyinaWwGN293XFhXiHOd+Ofu5ls9d0IIacvWzTc/RK3kV02nQJYgslloQG6NptWoC0IIGQibE2DnC9YM08yGl1VC1CJkSCUuFvKm1iAXybPrzke5OyBUQ4IQQvrOmaWtmsVWvDexvY03RWHzWPLfgnxhzrz7Rt7HJmHs2doLE2MXv0y83adhEELIFllZgF3EQ2atXOe7zVq4CyJNkMIgO38nomsvANEEkhwj09nji/YjwUmqOEoU75xu8YkYhBCyJZYPQ/N+F+FZmhW+W1fboarDrrSjIndDOMs3vXgVey//OTQ+KIrhRPbJGvuR4DAW3Fpd+COEkAGwvLKl09KiWKJAdO2l3Bq2bQJA0mkuptdeQKKz0C23DUBh+boQNckSHKQ3AZRD0774tWdXPkFCCNlVlo4DRjQpfLmQ/CGc6aX7ABOXS0pGE+xHgvTi1cKPW2S0iZktzEUHUAC3/+MvIXrnFejeQVFv4sphXkHtn37swTOf6Lpwfm6XsHGSKqM0CCErsZF7+1IRdusfLj0bzsSYGBQ+32mqkGRaZMzFMsuI63oRzoXBObeKm+9pprmFb4QlKgkhK7ERAXYlJd0CnaL8tImTVIti527BTSeHSGCKou6ugPs2C7KHMK7YfJYLr5vvuUiKYvTLPgmaEEKAszyW3q+ZoFkhtK5d44OiOE8sM0s2VSAyUoSaRZ6ACezt/YMfwd6r3wOQF2TvAmf1uvC5icmF15+vO3euERJCVmFl6Yg1KR4dP1UDiAnWcQhxnGRFnYUT++y1l4/yLDhXT8KvHbHNWhAAiuJCId+uJMdztSBEaQITQpZnpWpoifXjAigSLZzlW7gMsgRHicLcvDb3pOSDeFZy0omcW3Dzx7j9F34FR7/+y0geenSFU1seZ8GnKBcceuc092MnCryL/WK+7tzT7Wd0E0IGwNJREFPY7DTrcpDpEV4+Suae4TZFXBRTd9vi738dmeaLbzOXRK66X/zas6UxXP/Dn/tl/KPDD5zxNNvh5hZpgpupFvN59zQ/10zzAvFuvu7c/QLyhBDSlqWroRV1GcQgBqCTQ1yplMR1dYCzcxdxxWtPHsqrmvnWpUuyCIWaOcv3Xx19O3cL2EfVO9/sunHzBmKc89pn4XP5Obn5KjB37tX51b0OvSeEjIv+3DuLKT19IpRtty7OMm6T4Pqkil76jtfxuTNsmpCc/ggw8sgK8+4bRUGfXbQeBTN/eHV+LsoDyBcv++o7Pmviif85EDJmeqMAzvLKbrmj9Gw2rfysg7MIu0vWcK+rcyqs48pjjvx+uyxObSNd2oxDyNjpjQCX/mBNDEmOi1v8XbOImyrCnVjzUVC2JP2500IkZBz0RoCrZPEBxD6ZeN0W8FnGEaA2NVmRJ3M4mizgXbmYVFnX5xwagxcdMjZ6JcDVxS01MSRL1m79rmus6rx8y7YqZLsquFXW9VmHxujLZ0DIuuiVAAcxcZG5tut/wE0WXlWMaQ0SMnz6L8DAXIjaJjmre8Lhkj5C20Lvh8QmQwgJ6RO9F+Dilj4+gBxf73QO6+QoGZ9Eje+MydhZOhNuF5DKb4ceXMhTpG3GXDXrLLTPqsde1LZo/0U1hA9trFd13k0LdLu8eOdTN8c+zJ2QddJ7C7iKiilqCQOz2123eDQWK8u/zR/LORPSNwYnwMAsOiIkPEMS4abz8KMVhnTOhAyJwQlwYemaGJJOS0LkW8NtqQsb2xahsK86d0rV5eLPfRu39+uOxyZk6AxOgAFPAOxTl9fBrmba7eLjkHblcyJk1xmkAPtofFDULgbqrcQ6uhaT0Bz9BbyJKff1qV2sXMvM5tmlixQhfWBwAuwvuBVtYgp3RFO/XaSpotqi/eqqljVl5BFCtsfgBNinJC7RBPBqR4T69oG2VdOaqrJVs+5C4XqEkM0zOAEOWYjFItyChblFgtaVOFXrSbSJo617XR2jqSbFqudMESekHYMTYGAmpkGxWnFhblfcFb4oNrkZgPpyl5ukL8kghOwCgxRgoHnhyS3MVfv2QTiqF4Km0pdN21Y5Vtt9CCHtGKwA+wStYTFAZWGu7yxTUc0PX2sqIE8I2RyjEGCfUrJANAHS6eAjAULnt+d983WF4Qkhm2UUAtxUQEejCSSd1u7bJ3EOhdmFtrV1LSx77n35nAjZFUYhwI66hTknwn0S20X4LohF7oimhTzRHUy1I2QgrFWAU63/Q0909rj2ZVCc/THoDidEoThYjSaIbl6bVVLzs+cqItQXkW5TwrLpGXZ5h/B/kU09F46QMbGaAHvlHn0iAaYNanmcLG9Npbo+AQaab7+z/Vuw98p3Zu/tUyvk5EZgYvVui11jmWpo1ad01BWGF83mvpdlq66xShsZO6sJsJmv4+6sX5GwvMUCnN9b8nCaQQBMkAt+U0bXWREAMDFO7/lRmBuvIYVBqvlFRQ8uzDqmU4hm0GiyePoNc9y28LSN9qhaw64wfJW05r9OW9+y/5uQsbI2F4Sx1u+kZsRlLdlU8z/ySICbGgEIP2V4rWgGaIbTW+7E5J2XMDHzIpGaSf4Mupq7AJ+2WWu7SugC4t5HMi/WywrqriS3ENIVKwlwqATiaQbsR1JbHjESINLFolVMTPI/6ESBc4ue37MmVExh+aYX7kH82jO5tRucYC+f5rQ0wRhhe6GqQpcCIcuxkgCHPAkTk1utdRbwSaqtRKtwZSB3W9TcAbeibVRDtZ+zfJPLDyG6/tLMStMst/qs/7drsdlEiFi1n6l8/gLkC3Niap84sohtlcckZNdZSYDrQpMiZOFIB82wHwmefOVmq/F9V8U7pxlePmpvOfsse4srmiE+uV7sq2KQXrwK89TvAVkyW4yLJoiuvWAn291i3C5anE2RMCFYO4KMmdXuo2tCkyAmOGAKgwjAj911buHQRnJ3hvurPIwNzsWruapDf9yhtuK9mNKCm2vPfuxnYI7fgU7s/DVDevGqnXB3rohVrdomlvGxh7Yv4y3yq9BRhMkY2UoiRqQJkCW42TIMTXUmGYqOaxVoBpkeQSfnEL39Yt5m60jUhqiNlGUTWSi+ZOysTYDdrWcw0sHEgIkXhqGlan2/1vEoyTGmqeLg5NrC44f++JtSkJsnMp1FbVirWOOD3Cf85g8AAFOZIEKG7OBCoxtik/V027pY1hltoJWfsxzHD0fbNVcKIdtgJQEO/bG4JIy6W9A2D480ki/WRZL31/gAh7EgO3dx4b7rEJlCCEw8C7Ny4mr93untD2DvxW/ki402HC019THBq8xpVyzCJmEMzXFZl4hgvRcHQvrGaotwgRjYRIGDtGaRLUswMVi4mJYq4HTa7//Frz27yjSXRoC8JoSYWbafW3CzbodpBpze++N47Z//QyQKRNdemosUGCuhVO+mfkBZjGkFk7Gxtky4WACdHAa7q4mRKnDlXPPhYplNyO//ix+9a6VprkQ0KbLiABuOdvFq7l4xE0wkv0Rc+vy/wf71F3DjwnsGW7AmZJ36bXXhZIuuR6HkDoovGSNbWYQzyTEiZJhq8+FOUkVkzUm//0vJ4rTfTeGsYhei5qzj/befw43z9+L8t3+/Ni13rKySEUfIGNmKcmh8AIipTdJwxEYQO/+vGEwkL/hy5aDlcc4+1TDRpFiMk3SKVGIkl+7DOaOYPvLXsXfzrU0duVcs69MNuSEIGRM7lU/rFvAKoZY8fhjSbppb+QOOJvmcgGJ+2eGlfJHOZodVy0D6vxF4PTbGeM6EhOC987oQU1RvA+ZF96whW4SQ4dErAd7luFEFikgJwIqrJ8jV37s2f0LI9umNADuL0tWa2LUn+RbWbDSBOXqrSOJw21ySiXvSBq1fQkhvBNgJViwtHqPTEe6acHruEqJsVjFNYQu7I19c9Nv7wFwR/IGG3RGybXojwMC8C2JTArZq+rCzcCNBbgnfvDYXI+yuG4VF3DPELjYuoukz7MuFh5BN0xsB3pXb9qbU3MIPbF9nB7fWFuvZhXNZBW0hvkBzZltfz52QddMbAa6WSdxkFEHTuHXbQtEOKgbZwYUiddtfiHPpy6GFxV1zTyxTorJuP0LIPL0R4F3HF1fnXiguEiaGJMdQ2CeDeP0IIeOlNwIcsgrrYmy7wJ9LseDmtWfxASSdYj+S0oMtfeH2kzZoPRIyfHojwL5AAbklWU102AUhdgtx7rUT0mmqRYha0Tc5BkBLmJCx0hsBdjhB27dxaILN+4Tb4B8/FGrmLN/08FLhE9b4YM5iblvScZvMhaERQtZCrwQ4VRSJDA4nVsUTLNbAWZ5i4UdrpJonjFQLkMPEMDfnn/Lh+vY1RG0RxecaqCdNYSdjpDcCrMhLVKYwG1/AOov4+fvGMvPzuvk6t0N27iKit54r+X0j5KnLsSAoUrvAMp9NqP6FAnmt5cpjnIZ4wSFkETtVDa0JAWbV0TDLhJPK+64IFS53KLywM+t2AID00n354+5dgXvPdQET74worTqP0H5+yjYhY6c3FnCf8f3DfohaIbQ2RM23GHet1gUhZP30UoCr4We7rlW+uFYX3AAbopYlZXfErpi/hJCN0UsBrivx2Ccyu6BYYC3hPp5LiKZIjmXbCRkqvRTgvuG7HJyFawSlcpUK++imHV182wZDufgQ0pZeCnCXyRerhqj5GW4u1Ky6DUAeouYla4wGlrgkI6SXAlxN2XWieJRsXo7XZaX5oWbVBbfs8BKit55b05F2i9qqdmJAJwQZG70UYCBcoesw3s2b2GqWnm/t1i24pZfuK5JO2lZK26XFybrMxOZvaDe/P0I2RW8FuI5Qsfa+2lUqpkhb9qnKVFWg/TsEQsjuMjgBrmPX6iu0RU0cFGGfaj0Mii8h/WBwAhwqzLMLxXrOQkiEqy4J/6f6NIo+XngIGQODE+AhMVfUfXoEwayouyJfwKuzgFMNCzIhZDcYnACHLMJdWpxahrnzmBxCkuOiFCcQLv7ubyOE7C6DE+Cqu8FZiH0R3Sr++ZykCo0PYG68NjtHzzUx5/+16c1Af90vhAyZwQmwjx/i5YtyX8SoWlFtYou6Z+fvBLKkVDWt2rco++i1bYJQxElfL3aEbJvBCnBIBFYRhy6z7nzmQs9MDFPzyPsqmyzwXnoSCPpzcSNkFxisADshmGbltmWLue+K+8JfcHOuBd0/X5u27AvhNnzBvvhShAlpx2AF2DEJnOGygtRlGJsv/sWCm+dayA4vFXUUfBFsugPYRIgaRZeQ5Rm8AI8CMUA6nYt9rlr6VXcBQ9QI6ZbBCXCdlRfy4Xbt111E3dxCqciIJrOFOdtWZ+nX1qYghGyVwQnwMplwu54dV1fMJpT5BpQz5qquiFBctAS2E0K2x+AEeEhUC+vM1Q72+hVCakW4rZ+3rrQnIWTzDFKAQ1af3970epeoim+dFV/9rSaGmR7N+jQUO3chanXiTgjZHIMU4CpuQcqJS0h0dlGEq6FdxaKaDUOrbiuJ8OSwyJiDmFqfbyQY9WOQCOmSQQqwL0Tutb8g1ZfsuOqcjMxnuDncOflWfyljDmELuppNt2zdjOodBn3JhLRnkAJct7C0KDuuSTi6FpaQCyKUUFJcUFxssM2YC/UNRYqE3B5tEc128kJGyK4ySAFu8pmGRLQqOtX+u2DR+S4ToD7MrNjuPXE5mZxHnB4vHL/kwsDydwXumISQdgz2L6YuhKvNIlZon10IWatzmdRdIFx7JMgfee9lzLU51qpz6vpzIqQvDFaASQAvY87VyFi2NgYhZH2MQoDrwtB2xb2wDupcK3N9bMacq5HRlDFHCNksoxDgkG9XAJzWh8f2ijZZf27bNEMe9XByY9YvnW58joSQeUYhwEBYlCZmWP7KJl+w+9kz9vX+eSCdFlZx27AzQsj6GI0A19VEGIIrIpRK7Py6Icu/2C+alDLm+vwZENJHRiPAIfwnCmc9Vh9feN1rl+FWjeIoMubsYpxODhG99RwAr+A7IWQrjEaAq9lxgnB2XF/xq58VMcM2Y86/thTnad0OCiC5dB8knRYF390+ba5JoYVNQkg7RiPAZB6nxZkCiCYwN6/NLeglVFRCNsaoBLgpFdn9VGNi+2bVhcpL1hn3pUQNANm5i5Dj6yUrOF7yzqDnNxKEbJVRCXBTdlzILeFv7xNVK7buAhLKYssOLhRF3X1CdTVC4h5K5a6jbxc3QtbNqAR4jLQtt1mqM2FimGS+doR/d5Dp/MJlHy9WhHTJ6AQ4ZHWFQtGaQtT6ZLktU1jH75PFB6UQNaB8d2Bk/m5hlZRmCjYZM6MT4DpcKcUjb9Wpz6FpjmXEN9GyO0Ynh4iuvVBs95+sIQCgWTjCojJukwtkAB8xISszOgGurYhmSyke2lUn5w+u679I1JoSPLaZ/LGMhRlacEsvXoUkx3m4mvffxVm7ixb46rbV+ZAJGRPzj1Yga6FJWHZddHzxFOTuCMkSxCYuRU5Ur99+MkjofZVV6w4TMhRGZwGTdhT+XM3y1yaGWJ9wpgjWFm6qt1wHxZeMmV4JcGihbFvHqra3iSqo+kgXjb0rCHJ3hHPNxC5OeHIIc3w9jxv2nrjRROg8/QeIEjJmeiXAwPaKh1cXnZbeXzNUn0ZcGnvlkbsjU0APLiB68wfFIlwdpcW8yrZYELSgCRkbvRJgZ5n579eJs1oXiXwrAe3x89FcRqBfWU0ARMgt+vT2B2DefQMQs7AEJlCTrNHSgiZkyPRXJTaAAEj9j8QKTJ/jgM+Cs3IjuOcXmaItu+WOXIQr+KJbXYzzxx3D50fIIhgFUSGUitz0vivaRg9UF8na9nOfg3pulOKY9sKU3nJHnrbsRUeEIhua2trMjZChQguYACj7bKvtda/zSLS4qDscoinet+6YhIwFCnBPqbu9D/XbpMgpUBLhkKVb9AvsOxZ3DiEhBi3Abf+462pBhPpVf3cpHm3cEE3hcNV+1c/BX4xclNnm4oRDC25NiRi0fsmYGbQA1xFaFKqGTPmvq1ERu7CItFQGWU04XLCr93rRU0LcHNznoZPDIlnDH2/RXCnCZKwMXoCXsciq/d1rI+GxuhaOtheBNpa6b43W+WtDbVUfr04O5x5zX3fBqrvAETIWBi3AbYXkLOO0jSxoQ8gN0DSnNueyDh9waF6hcLOCaBIUYT+22LUBsydyUIfJ2Bi0AHfNIldFXwRnJfGOJoB9soYTVyPlsapiTsjYoABvmGV8n4L1WKxnpc4KrwtJqxVSEwO2zrIfDVHNkOvalUNIV1CAyRzVRUnfx+tb9dWokarFX6QcV2pG1I1ByNigAJNGmpIzqhEOIRHOO5raZA1awWTMUIA7YiiWny+g/msX2VC4HEwMsQ/6DAq1dVUQMiYowB0xlAW6OkILbu7JGlUUs/rChIyJXv2vXyRY6xKtRaFg1Uy4Rf2CBG7J/fF2yRr055V4lm0TwboRQP7I+5MbpT6+xUzImOiVAAcXeTZ0nKZIhLaZcI0iauJgrGzXERCLcMXUl/kM/IuKAND985Dj6/PjEjIyeiXAQIvqXGegaWW+LgHhTMeMJmfZe2OELPvSuXrF1OsiIvzf1YgKANCDC4VPuLo/IWOhdwK8SZpSZkPW99Czt/woh2XuPKqRDXVZcxofFHcBu1Bfg5BtQwGuMCcSCBeUGXL4VKguRN2dh0+1tkM14cJvK4gmhSU8xM+SkCYowAsICdA6/NCN+53hYaDbJhSG5mo7AJW7BHte1YU8jQ9gjt7a8EwJ2T16LcBVEdvGSrrvz/TbVhmnSjGmmLkFur7jsuIE4YW87PASkKVdTY+QThDV9qolIq8B+MHmpkNGzgOqeue2D8r/12TDvA/AH6nqJ6sblhJgQggh66PXLghCCOkzFGBCCOkICjAhhHQEBZgQQjqCAkwIIR1BASaEkI6gABNCSEdQgAkhpCMowIQQ0hEUYEII6QgKMCGEdAQFmBBCOoICTAghHUEBJoSQjqAAE0JIR1CACSGkIyjAhBDSERRgQgjpCAowIYR0BAWYEEI6ggJMCCEdQQEmhJCOoAATQkhHUIAJIaQjKMCEENIRFGBCCOkICjAhhHQEBZgQQjqCAkwIIR1BASaEkI6IN30AEXkEwBcAPAxAN308QgjZAK8D+BKAP1DVtemYrHGs8sAij8htD34T6RTm7p+AOX8FYqLyTxQBAEyp3RSvART9XJvx2ys/AEp98/5i+wtEBMZI/toIjHjbjNsGiG2v9hVvLGPHEzPrZ8+7eB/ZHwDF68gI4rl2U9oOAJHM9w+NVdsus/OLJJ+ja3PvAcze2+OV9ina8vkU/Y3bz+tv20QAA/tb3P+Fcptt9vrbz63Uf9Ym8Mfy2/L9AUBUAc0gmgGaAe7/danN/mSZ3cfr6/rZfdyPeP0xdwzbnqXQLAOyNH+dprN9shSa5e/Vbs/bM9svne1v+wfHs/u48TTLoGm+j6ZZ8V7TDJnt77dpliHz+1f3yfxtaWlbPi2/v87GyhRZql67zvbJtOiraf56NlZ5v7xd7TzV9s9f5/NSZKpIVZEqkNrv1732fze3z9oy1G+rHuNtnOIJXMO7SPBhXMRv41WzDiFeuwXshFduvRfm8gdgzl9Z9yEIIWSr3IY9fByXCyG+F/vZ35S7zyzEaxNgCi8hZOisW4jPLMAUXkLI2FiXEK8swBReQsjYOasQLy3AFF5CCCmzqhC3FmAR+RG57cGnKbyEEBImJMQi8tdU9fFQ/2USMX5c330d5tKPQG65a03TJYSQ4XERMR7GLThGBgAfr+vX2gJW1f8iInfp8Zuvpt/7FqI7H4Fcei9EmExHCCEAoFC8iGM8gWu4BTHexOn7VfW7df2X8gGr6msARETupBATQkhOVXhfxMn7VY9rhdexUhQEhZgQQlYXXseZ4oApxISQMXJW4XWsJROOQkwIGQPrEl7HWmtBUIgJIUNk3cLr2Fg1NAAQkTvNXT/6avbOS5Bzd+RVxkRyMbavbT9bgcwAgtl2u831FZiiUtmsTQAJtPvbATuufSew2+Adwx5SMLePG7aYk7fP7JButPJ4prLPrCKYlPZ1ff05uUO4ql9FlTYpj1fbDv/YKNq8Uyk+E0H5PGZts3mg0i/v45/HrA+8Y8D2m9vuPuNgfzu2zOaCSj/xdhZVAGoLnnr/p1071KuQpnb/apvbbzaO+O1u7FJf2Mpo7rdC1RvX9Q+1ZwrAVWPzx1Kot19pjm6cYn947/PXGuqv6h1aoVm+j2ZqpzjbR7Uy5+LjsW2Z3x/l8b2x1J7X/DHcPO2YRbE5ne2jKB0f7iOGIvM/Ls3bKkMVH6nacefbZ/u4sTJ/m98O4BpOcR4xnsa7jVENy7JRAS4OInIHgAc3fqB6fhbAb3V4/C4Y2znzfIdPl+f8qqo+t+5BtyLAXSMif6KqH+l6HttkbOfM8x0+QzxnOmYJIaQjKMCEENIRYxHgX+t6Ah0wtnPm+Q6fwZ3zKHzAhBCyi4zFAiaEkJ2DAkwIIR0xSAEWkc+IyLdEJBOR2rAVEfmkiHxXRJ4Rkc9vc47rRkRuF5GvisjT9velmn7PisificiTIvIn257nWVn0nUnOv7bbvyEiH+5inuuixfl+TESu2e/zSRH5xS7muS5E5N+LyKsi8s2a7YP6fovskyH9APgAgIcBfA3AR2r6RAD+AsBDACYAngLwwa7nfoZz/hKAz9vXnwfwL2r6PQvgctfzXfEcF35nAD4F4HeR58g9CuD/dD3vDZ/vxwD8dtdzXeM5/zSADwP4Zs32wXy/qjpMC1hVv62L0wV/EsAzqvp9VZ0C+E8APr352W2MTwP4Dfv6NwD8re6msjHafGefBvAfNOfrAG4TkXu2PdE1MbT/owtR1T8E8GZDlyF9v8MU4JZcBeCnFj5v2/rK3ar6EgDY33XPjVIA/11E/lRE/v7WZrce2nxnQ/pe257LXxWRp0Tkd0Xkke1MrTOG9P2utxraNhGR3wcQejLoF1T1v7UZItC20zF5Tee8xDCPqeqLInIXgK+KyHes1dEH2nxnvfteG2hzLk8AeEBVb4jIpwD8VwDv2/TEOmRI329/BVhVf+aMQzwP4D7v/XsAvHjGMTdK0zmLyCsico+qvmRvyV6tGeNF+/tVEfkt5Le5fRHgNt9Z777XBhaei6pe917/joj8qohcVtXXtzTHbTOk73fULog/BvA+EXmviEwA/B0AX+54TmfhywA+a19/FsDcXYCI3CIit7rXAD4BILjavKO0+c6+DODv2dXyRwFcc66ZHrLwfEXkitgaoiLyk8j/pt/Y+ky3x5C+38FGQfws8ivlCYBXAHzFtt8L4He8fp8C8D3kK81f6HreZzznOwD8DwBP29+3V88Z+Wr6U/bnW30859B3BuDnAfy8fS0A/q3d/meoiYLpy0+L8/0H9rt8CsDXAXy06zmf8Xx/E8BLAE7t3/Dnhvz9MhWZEEI6YswuCEII6RQKMCGEdAQFmBBCOoICTAghHUEBJoSMmkUFgCp9f1pEnhCRRET+dmXbZ20xrKdF5LN1Y/hQgAkhY+fXAXyyZd8fAvg5AP/RbxSR2wH8EoCfQp7c9Et1FQl9KMCEkFGjgQJAIvKXROT3bM2U/yUi77d9n1XVbwDIKsP8DQBfVdU3VfUtAF9FC1HvbSoyIYRskF9DnvzxtIj8FIBfBfDxhv4rFQmiABNCiIeInAfwUQD/2WZ5A8D+ot0CbQuz3CjAhBBSxgB4W1U/tMQ+zyMvju94D/IHQiw8ECGEEIvmFeb+n4h8Bigeg/QTC3b7CoBPiMglu/j2CdvWCAWYEDJqROQ3AfwRgIdF5HkR+RyAvwvgcyLiCld92vb9KyLyPIDPAPh3IvItAFDVNwF8EXkFuz8G8M9sW/OxWYyHEEK6gRYwIYR0BAWYEEI6ggJMCCEdQQEmhJCOoAATQkhHUIAJIaQjKMCEENIR/x+gzjiIRmio8QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAABsCAYAAABU1GiDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAANJElEQVR4nO3dfaxk9V3H8fdn7mXXsosschcWCru0tvJkLC7aYG8bG4rabrS1jTX4hxprUjG2gSZoHzC1sfpHoT5rNai10fTBGEvbEPpgpYkVpQHW7QNCgbZLgV3YbpEFli7sznz94/zOzJm5c+bO3Dszvzv3fl7JZM/8zu/he39z9nt/c2bOuYoIzMxs+hq5AzAz26icgM3MMnECNjPLxAnYzCwTJ2Azs0ycgM3MMnECNjPLxAnYzCwTJ2Azs0zmR6m8sLAQu3bunFQsY3Pi0Yfa2/M7zl11f+W1gqoZp3eMo8dbAGw5afjfbw8deRaAc0/d3LevUfsbl3L4UYaum6/qNZe9+wAe/Pa3OXz4cL9dEzUrx7XNpkHH9UgJeNfOndx2223jiWqCvnv91e3t03/nz1bdX11CKcfpHeOOg0cB+PGztgw9xrW33AfA+/f8UN++Ru1vXB595gQAO04e/lBZaQJeXFwcJbSxmZXj2mbToON6pAQ8KwYl3brk+HRa6m3ts9S7s6ZN3TiXbT4MQJPhE+aNf1D09f49f9VV/sJ/elfnSc94dYlunEZJvMuZ+tLWbI3zOWAzs0xWtbwZ9Fb7+1/2WwA8+V/dK7pTX/5WAI785190lde9BR/U10rUvY3vt/Jdrk2d5g/sGqk+1P9sg1bzdSvKuf+9tYjjostHjqPXoFX23P8V58Cbp63+PPsdB492ne822wi8AjYzy2RVK+BBK8M3/+7Vfcuf6Fn5lvqtfEv3fX71H6SVVvIB2XNpYbZpDL+udOIYADH/fUPFVb4zgKVzVPeu4W37zynKL1o6fl2bupXuoPO2b/vv76W+usvrzpkPcsUb38Hxbz4ydH2z9cArYDOzTJyAzcwy0Sh/kujS3btj2O9L1n3YVmfQh3ArOW2wEazka2jT+OraSr4H/FwLfvIVi/zP3r1T/7baKMe12agWFxe5q+a49grYzCyTiV2IMezKt3TDgA/hLlx43kh9NStLrbme3zt1F1ys5IqvOoNW83Xj17W55uQL29t/+sw9XfvqloqD3n3UtalbtQ66QKUu5sdq5rL6LbPeDzQ3NbwasI3Hx7yZWSYjnQOW9B3gwcmFYxvcrojYPu1BfVzbhNUe1yMlYDMzGx+fgjAzy8QJ2MwsEydgM7NMnIDNzDJxAjYzy8QJ2MwsEydgM7NMnIDNzDJxAjYzy8QJ2MwsEydgM7NMnIDNzDJxAjYzy8QJ2MwsEydgM7NMnIDNzDJxAjYzy8QJ2MwsEydgM7NMnIDNzDJxAjYzy8QJ2MwsEydgM7NMnIDNzDJxAjYzy8QJ2MwsEydgM7NMnIDNzDJxAjYzy8QJ2Mwsk/lJDyDpYuA64HwgJj2emdkEHAauB74QEWPLYxpjX90dSxdr23lfo/kcjTNfQmPrDtSY637MzQHQ6CpvtLeBdr2yrFEt73kAXXWL+kr1hSQaDRXbDdFQZV+j3AdK5b11VemrkfpTo1Mv/dzt53PpAbS35xpifkl5o2s/wJyW1u/XV225Oj/fnIoYy7LyOdB5nsbratMuK+Jp12+U7Sr1U5kEDdK/Ko+F7rJUXKmf5q2rfqdMVPuqlhXtARQB0ULRgmhBeVx3laVHq5XaVOqW9VKb8qFKfZaMkcpbTaLVglaz2G42O21aTaJVPI+0vyhvpXrNTvtUv29/qU3ZX7RaRLNoE81W+3k0W7RS/WpZtFq0qvV727Sq+5pd+4qwqvWj01craDWjUh6dNq1o141msd3pq7tdUR4pzkj1i+0irqAVQTOCZkAzvb7ldvXfweWdshb1+3rHeILj7OUIRznBbk7lZg41xpGIx74CLhOvTjmbxsKFNLbuGPcQZmZTtY2TuJyFdiI+m82tn9OZq07EY0vATrxmtt6NOxGvOgE78ZrZRjOuRLziBOzEa2Yb3WoT8cgJ2InXzKzbShPx0AlY0ou07bz7nXjNzPrrl4glvTwibutXf5QLMX4kjh6mcdqL0JYzxhSumdn6cyrznM8WjtECuLyu3tAr4Ij4uKQz4tjjh5r33c3c9ovRaS9A8sV0ZmYAQXCAY+zlCFuY53GOXxARX6+rP9I54Ij4DiBJ252IzcwKvYn3AM9eEHGsNvGWVvQtCCdiM7OVJ97Sqr4H7ERsZhvRahNvaSxXwjkRm9lGMK7EWxrrvSCciM1sPRp34i1N7G5oAJK2N8744UOtpw6i551e3GVMKpJx2k710h3IGiA6+9O+sq5otO9U1ikTqE95dT+kftMzkfZRGSMNKZa0Kbttx1Rp0xmy7K27v0ZPm84dwdTVtqxbjakcorzrV/suberur7ac6ti0yyo/SntORPfP0SnrxEFPvaJO9efo1KEyBqnekv3lHPetn/pWJxZ66qnSWBFApBueVo7pspyo3CEtUvvesrJdpx9Vy8u+u+qS7oxW/htEVPot6/crbwVQ3o2t2lcQlXZdMZb9tNtTeV5sR7/6EZWhg2gVbaIVKcROm4iemNvTk8pa1fp091/pK9LPtXSMMs7UZ/tmc9FpE3SNTznFBK3qdEVR1tNVe0oj9bu0vNOm7KtV3VctB45wnK3Mcz9HB36rYVQTTcDtQaTTgfMmPlDH64GbpjjeSs1CnLMQIzjOcZqFGGG6cR6KiIfG3elUEvC0SbozIn4sdxzLmYU4ZyFGcJzjNAsxwuzEOYhPzJqZZeIEbGaWyXpNwDfmDmBIsxDnLMQIjnOcZiFGmJ04a63Lc8BmZrNgva6AzczWvHWRgCX9s6R96bFf0r6aevslfTXVu3PKYSLpPZIeqcS6p6beqyV9XdIDkt4x5RhvkHSvpK9IuknStpp6WeZyublR4c/T/q9I2j2t2NL450r6gqR7JN0t6eo+dV4p6UjlOHj3NGOsxDHwNcw9lymG8yvztE/Sk5Ku6amzJuZzRcovP6+XB/BHwLtr9u0HFjLG9h7g2mXqzAHfAF4IbAK+DFw0xRh/GphP2+8D3rdW5nKYuQH2AJ+muETjMuBLU47xLGB32j4FuK9PjK8Ebp728Tfqa5h7Lmte/0eBXWtxPlfyWBcr4JKKy7h+Efho7lhW4aXAAxHxzYh4DvgY8LppDR4Rn4uIE+np7cA50xp7CMPMzeuAf4zC7cA2SWdNK8CIOBgRe9P2U8A9wPOnNf6YZZ3LPl4FfCMiHswYw1itqwQMvAJ4LCLur9kfwOck3SXpzVOMq+ot6e3cByWd1mf/84HqFTcPk+8/8JsoVkD95JjLYeZmzcyfpPOAHwW+1Gf3T0j6sqRPS7p4upG1Lfcarpm5TK6kfnG1FuZzZGO9Gc8kSfo80O8P0V0XEZ9M27/E4NXvYkQckHQG8G+S7o2I/5hWnMBfA++lOPDfS3G65E29XfRpO9avqgwzl5KuA04AH67pZuJz2ccwczPx+RuGpK3AvwLXRMSTPbv3UryNfjp9DvAJ4MVTDhGWfw3XxFwCSNoEvBZ4Z5/da2U+RzYzCTgirhi0X9I88Abg0gF9HEj/HpJ0E8Vb2rEmjeXiLEn6W+DmPrseBs6tPD8HODCG0NqGmMtfBX4WeFWkk2x9+pj4XPYxzNxMfP6WI+kkiuT74Yj4eO/+akKOiFskfUDSQkQcnmacQ7yG2eey4jXA3oh4rHfHWpnPlVhPpyCuAO6NiIf77ZS0RdIp5TbFh01fm2J89Jw/e33N+HcAL5b0gvRb/0rgU9OID4pvGQBvB14bEc/U1Mk1l8PMzaeAX0mf4F8GHImIg1OIDWh/DvH3wD0R8cc1dXakekh6KcX/w+9OK8Y07jCvYda57FH77nYtzOdKzcwKeAhLzg9JOhv4u4jYA5wJ3JRep3ngIxHxmSnHeL2kSyjexu0HfqM3zog4IektwGcpPvX9YETcPcUY/xLYTPGWFOD2iLhqLcxl3dxIuirt/xvgFopP7x8AngF+bdJx9VgEfhn4qjpfh3wXsLMS4y8AvynpBPA94Mq6dxoT1Pc1XGNzCYCkk4GfIv1/SWXVONfCfK6Ir4QzM8tkPZ2CMDObKU7AZmaZOAGbmWXiBGxmlokTsJmta+mq00OSxvJVSUmfkfSEpJt7yj8k6VuVmwJdslxfTsBmtt59CHj1GPu7geKrhv38dkRckh77luvICdjM1rV0efXj1TJJP5hWsndJ+qKkC0bo79+Bp8YRmxOwmW1ENwJvjYhLgWuBD4yp3z9MN9v6E0mbl6u8nq6EMzNbVrpR0suAf0lXAkJx9SeS3gD8fp9mj0TEzyzT9Tsp7le8iSLBv72mrzYnYDPbaBrAExFxSe+OdPOkJTdQGkblPhnPSvoHipX1soGYmW0Y6e5p35L0Rmj/6aWXrLbf8mZb6cZAP88QN6jyvSDMbF2T9FGKP1u0ADwG/B5wK8X9uc8CTgI+FhEDTxdU+vsicAGwleKua78eEZ+VdCuwneI+yvuAqyLi6YF9OQGbmeXhUxBmZpk4AZuZZeIEbGaWiROwmVkmTsBmZpk4AZuZZeIEbGaWiROwmVkm/w8q0QeUAgpMzQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADqCAYAAAABMJbOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQ80lEQVR4nO3dfZBddX3H8c/33s1uNhDYJBtCNCSQlPCQCGhFUUQLAgW0Dg/W6rSFKdNOdfrEUIogtUadTg10KOPItDqOTC32Aa1Wp4VmBC2WtIolEMIyPD9GIbIBAgnZZPeeX/84D/fcvXeX3ZNzvnuXfb9m7uy9v/M753x/59587i/n3j1rIQQBAHzUZroAAJhLCF0AcEToAoAjQhcAHBG6AOCI0AUAR4QuADgidAHAEaELAI56ZrqA2WhwcDCsWrlyWutsuffe4RDC0opKAkox3df20888o+HhYauwpDccQreAVStXavPmzdNap3/BgqcrKgcozXRf26eeemqF1bwxcXoBABwx0y3BbY+/LEk6beUhWdt7N9whSbr1mtNnoiQAXYqZLgA4YqZbgnPXDLS1bfnLs/wLAdD1mOkCgCNCFwAcEboA4IjQBQBHfJBWgt/91pAkafVhB2dtX7jqWknSSRd8dEZqAtCdmOkCgCNmuiX46ofXtbV96n9ubHnc/xmvagB0M2a6AOCI0AUAR4QuADjinG4J0m8vXPvBY7O2G3/8rCTpj989vevuAnhjY6YLAI6Y6Zag07cXPv0rR/oXAqDrMdMFAEeELgA4InQBwBHndEtw+vXxH/I7bvXirG1Nch2GE5cf0nEdAHMTM10AcMRMtwQ/vJw/Qw1gapjpAoAjQhcAHBG6AOCI0AUAR3yQVoJ3bPiBJOnfr3xf1vb8nv2SpGMWz5+RmgB0J2a6AOCImW4J7t5wRlvbYQv6Z6ASAN2OmS4AOCJ0AcARoQsAjjinW4IjfuNLkqT1Z7wza9t44XpJ0uL5HGIATcx0AcAR07ASPPsvfzjTJQCYJZjpAoAjQhcAHBG6AOCI0C3gyZf26ndueUBP7hrVk7tGs/bdo1F2S42F+AYAEqELAK4IXQBwxFfGCjhqUb9u+sj6tvaD57W/h/WYR0UAZgtmugDgiNAFAEeELgA4InQLGNr+itZffpuGhkc0NDyStX/vkRezW+rlfZFe3hd12gyAOYjQBQBHfHuhgHUrDtHm689ta//Q2sVtbQN9vK8BaCIRAMARoQsAjghdAHDEOd0ChveO6Wv37dCFxy2V1Dxve/E/b8v6fP2jb5EkvbC34V8ggK7FTBcAHDHTLWCwv0eXnrSsrT2d3eYt7a97lARglmCmCwCOCF0AcEToAoAjQhcAHPFBWgHDe8f09/fv0MG98eG76NglkqS9jeYfQ+uvx1cv39fgD6QBaGKmCwCOmOkWMNjfo0tOaP/KWDq7zevr0AZg7mKmCwCOCF0AcEToAoAjzukW8PjOvfr1r9+nU9cOSpIuO2WFJOmnz+3J+py8/CBJ0pO7Rv0LBNC1mOkCgCNmugWsWdKvb158Ult7OrvNO+rQeQ4VAZgtmOkCgCNCFwAcEboA4IhzugXs2tfQpide1rGD8TncVYfE52033vVM1ueT71kpSRoaHvEvEEDXYqYLAI6Y6RZwaF9dv7p6oK09nd3mrRuc71ARgNmCmS4AOCJ0AcARoQsAjghdAHDEB2kF7B0LGhoe0YZbH5Kk7FeC90fNPr3J29mro5EAIMVMFwAcMdMtoL/HtG5wfttFb3o7vIUtnMf7GoAmEgEAHBG6AOCI0AUAR5zTLWDHnlHd8OPt+sJf3SxJev67V0mSHnlpf9Zn7aJeSdLX7tvhXyCArsVMFwAcMdMtYNlB83TZKSt0WTLDTaWz27xLT1omSfoDl8oAdDtmugDgiNAFAEeELgA44pxuAS+NNPRvD7+o7279uSTppo+slyT9/rcfzPp8+cLjJUnHX3arf4EAuhYzXQBwROgCgCNOLxSwaH5d5x+zWOcfs7ilPT2lkPfgDedJkvq/4lIagC7HTBcAHBG6AOCI0AUAR5zTLWA0kp5/bUx7R4Mk6chD50mSzrh+c9bnB5efKkla9bG/9S8QQNdipgsAjpjpFjCvJh2+oP3Q/TCZ3eY980+fkCT1L/jTyusC0P2Y6QKAI0IXABwRugDgiHO6JRgaHpEkfeKmn2Ztf/Jr8W+nXfnX/zkjNQHoTsx0AcARM90SrBucL0n60Z+d1rbsoq/+piSp/x9/z7UmAN2JmS4AOCJ0AcARoQsAjghdAHDEB2kFhOSWaiQPrtn0aNb2mTN/SZJ09nX/7VcYgK7HTBcAHDHTLcCSW6onebDxnKPb+t519XslSf2fr74uAN2PmS4AOCJ0AcARoQsAjjinW8BYkHbtj9TfE79n9SZvXf+w7RdZn99+y2GSpJ/tHnOvD0D3YqYLAI6Y6RbQY9Khve3vV+nsNu/NB3OIATQx0wUAR4QuADgidAHAESccCxgL0s6Rhgb66pKkWvIbaX9x++NZn8+duUaS9Mr+yL0+AN2LmS4AOGKmW0CPSUvm19vaP5/MbvM6fcsBwNxFIgCAI0IXABwRugDgiNAFAEd8kFaC9KI2X7zryazt0++PP1R7jgveAMhhpgsAjpjpliC9qE2nP9dz9KJe73IAdDFmugDgiNAFAEeELgA44pxuASG5jffQzn3Z/WOX9EmSdo9ywRsATcx0AcARM90CLLmNd1wyu81bOI/3NQBNJAIAOCJ0AcARoQsAjjinW8D4by80kgfff+LlrO2cNQOSpP2NTt9zADBXMdMFAEeELgA44vRCAeO/MtaTPDg3OaWQ11fv9OUyAHMVM10AcEToAoAjQhcAHHFOtwRDwyOSpE9978Gs7ZZL3iqJC94AaMVMFwAcWQh8eX+6zOwFSU9Pc7VVIYSlVdQDlKXAa5vX9TQRugDgiNMLAOCI0AUAR4QuADgidAHAEaELAI4IXQBwROgCgCNCFwAcEboA4IjQBQBHhC4AOCJ0AcARoQsAjghdAHBE6AKAI0IXABwRugDgiNAFAEeELgA4InQBwBGhCwCOCF0AcEToAoAjQhcAHBG6AOCI0AUAR4QuADgidAHAEaELAI4IXQBw1FP2Bs2sT9Klki6pYvsA4MAk/YekG0IIL5a64RBCORsy66uteNdINPyQagNHqrb0eNXm9ctqdVmtHvep1WX1umpJW7O9lj1O+2T9a536t/Ydv6xWM1nNZGaq1Sxpj9tqZtl9SUlfZX1bluX61pLtxdvShPuod7hJUk9LWy1rz9pssv7W3r/TMovHV7e4vro1x5G2pfez9lxb+jheltteLV7Wso+kzdJtKb5fM8U/lRyTtE3N/ul/r8ws11/Z8W32b/aRWrdhIUghkoUoXhgiKd/W0h5J0fj2cX1z/S3pr477iKSooRDFPyXFjxuNeJ2ooZC1R/H9KF4W98mtk24jvZ+0j99evOu4T2gkt+S+JIVGpKjRujxeJ1LUaF0n7R8akaKsfyPrny1r20dI2oKiRhjXHpLtBYVGSIYRKTRC23rN/SfrRVGuPdlOul6STY0Q1AitP+P2fFvz/mTrRB2WdeyvoEe1Ww/oVa3SAm3RriVlhe8Bz0TTsNX8RdLYiHqOPk9W7y2jNgCYEXWZjtVCHa2D9ah2a4nm7fxlGyglfAuHLmEL4I2uivCddugStgDmmjLDd8qhS9gCmOvKCN/XDV3CFgBaHUj4Thq69SPeHQhbAOhskvA9KITwWqd1Jg3d8OJjsvmHqrZoDYELABOoy3SkFmiXxvSY9kjSYkkFQve1F2r1w0+KGs/8SOobUH3ZCbK+hRWUDACz04ga2qZX9bRe03ot1Csa6w0hjE7Uf/LQjX9zwszM6stOJHwBIDE+bHdqtPe/ws4JwzY1pW8vEL4AECsatqlpfU+X8AUwVx1o2KYK/UYa4QtgrigrbFMHdO0FwhfAG1XZYZsq7SpjUpy+9dVnRdGOrVLPAllPn2Qms1raIXmcXrErbVfcJ21PL11lJlOt2ZZe5UsmWWu7taxjyTYtudZVvA8lj+P1miXF3ZvrNMfTeZ3mUNIraTUvg2XJ1brSdbKrcCX7SNdrtrf2TXeRXU0rvRJX8xi3tTfrSPbRPNTZ8Wq2NY/J+Lb2bbX2ydbJxpE7Hvljlj5HLY+TPq2Hcdzy3HHssM1su+lVxhSk7OWb3EnbWx7HP02h+Xh839x2WredXyfdVpSs1rxiWQihub3cPlvao6D4OlfJvtNthKCQr3f8TYq3k2yjeV/ZspDrH1r2n/yMmsckRCEps71/+jhdL9t/lPZXbkz5IYdsefs+1NxmulLUYZ201nR5dtSDovRwBWXHKjQ31Tyc2S7CuPbWdYKUbTPdR9YuKVLQTu3Xei3UnXpx0m8jTFepoZttNP4Xuk5SX4fFF0j6Tuk79cUYugNjmHmzvX6p8xiCpG1lhm2qktCddIdm/xdCeLvrTkvGGLoDY5h5s71+yX8M/LkeAHBE6AKAo5kI3a/MwD7Lxhi6A2OYebO9fsl5DO7ndAFgLuP0AgA4InQBwFEpoWtmi83s+2b2aPJz0QT9zjGzh83sMTO7airrm9kJZva/ZjZkZtvMbH4ZNXuOIVm+0sx2m9kVVdRf5RjM7Cwzuyc5/veY2RkV1N6xptxyM7MvJsvvN7O3FR1PVSoaw3Vm9lDS/ztmNjDbxpBbfoWZBTMbnI1jMLM/SpYNmdm1hQsM2W+gFL9JulbSVcn9qyRt7NCnLulxSasl9UraKun4ydZX/GvK90s6MXm8RFK9jJq9xpBb918lfVPSFVXUX/Hz8FZJb0rur5f0s5LrnrCmXJ/zJN2m+BfVTpH0kwN9TmbJGM6W1JPc3zgbx5AsP0LSJklPSxqcbWOQdLqk2yX1JY8PK1xjSQN9WNLy5P5ySQ936PMuSZtyj6+WdPVk6ycH5+aqniCPMSSPz5d0naQNqjZ0KxtDrr9J2pm++Eqqe8Kacm1flvSx8WM90PF0+xjGrX+BpG/MxjFI+pakEyU9pWpDt6rX0i2SziyjxrLO6S4LITwnScnPwzr0ebOkZ3OPtydtk62/VlIws01mtsXMriyp3k4qGYOZHSTpk5I+W1HdeVU9D3kXSbo3hLCvtKonr+n1+hzoeMpS1RjyLlU8Q6tKJWMwsw8p/t/R1rIL7qCq52GtpNPM7CdmdqeZnVy0wOn8CfbbJR3eYdE1U91Eh7bX+75aj6T3SDpZ8d8busPM7gkh3DHFfbYWMDNj+Kykvwkh7DbrtPr0zNAY0n2vU/xf3LOnuK+pmkpNE/UpPJ6SVToGM7tG0pikbxSqbmpKH4OZLVD82iz7NTORqp6HHkmLFJ+OOFnSLWa2OiTT4OmYcuiGEM6caJmZ7TCz5SGE58xsuaRfdOi2XfF5ndQKST9P7k+0/nZJd4YQhpP93CrpbZIKhe4MjeGdkj6cnHgfkBSZ2UgI4UuzaAwysxWKLwpycQjh8SK1T2Kyml6vT+8k607leJSlqjHIzC6R9EFJ7y/yj3waqhjDGklHSdqaTDpWSNpiZu8IITxfavWT1zeVPpM9D9slfTs5/nebWSRpUNIL066wpPMo16n1A4trO/TpkfSE4icgPUm9brL1Fb+zbJG0IFn/dkkfqOhcUCVjGLf+BlV7Treq52Eg6XdRRXVPWFOuzwfU+uHH3WU8J7NgDOdIelDS0qpqr3oM49Z/StWe063qefi4pM8l99cqPg1hhWosaaBLFM8+H01+Lk7a3yTp1ly/8yQ9ovgTwmteb/1k2W9JGpL0QMX/aCobQ67PBlUbupWMQdKfS9oj6b7crfCntxPU3lZT8kL/eHLfJN2YLN8m6e1lPCezYAyPJf/A0+P+d7NtDOO2/5QqDN0Kn4deSTcrzqEtks4oWh+/BgwAjviNNABwROgCgCNCFwAcEboA4IjQBQBHhC4AOCJ0AcDR/wPw9SREbazgRQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for A in [ K[np.ix_(Ia,Ia)], K[np.ix_(Ip,Ip)], K[np.ix_(Im,Im)],K[np.ix_(Iai,Imi)],K[np.ix_(Imi,Iai)]  ]:\n",
    "    matplot(A)"
   ]
  }
 ],
 "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": 4
}