Negativity_noninteracting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import block_diag\n",
    "import random\n",
    "from math import pi, sqrt\n",
    "import sys\n",
    "import os.path\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generates single particle states\n",
    "def wf_gen(N,N_pt,BC,t1,t2):\n",
    "    h1=np.ones(N-1)\n",
    "    h1[0::2]=0\n",
    "    h2=np.ones(N-1)\n",
    "    h2[1::2]=0\n",
    "#     print(h1,h2)\n",
    "    hop= np.diag(t1*h1+t2*h2+0j,1)\n",
    "    hop[N-1,0]= t1*BC\n",
    "    H_t= -(hop+ np.matrix(hop).H)/2 \n",
    "    energies, evecs= np.linalg.eigh(H_t)\n",
    "    return evecs[:,:N_pt]\n",
    "\n",
    "def exact_renyi_calc(r,GA,epsilon=1e-9):\n",
    "    chi0, _ =np.linalg.eigh(GA)\n",
    "    chi1=chi0[np.nonzero(np.abs(chi0)>epsilon)]\n",
    "    chi2=chi1[np.nonzero(np.abs(chi1-1)>epsilon)]\n",
    "#     return -np.sum((1-chi2)*np.log(1-chi2)+chi2*np.log(chi2))\n",
    "    return np.sum(np.log((1-chi2)**r+chi2**r))/(1-r)\n",
    "\n",
    "\n",
    "def exact_log_neg_calc(Na,Nb,GA,epsilon=1e-9):\n",
    "    \"\"\"computes log negativity of 1d free fermions\"\"\"\n",
    "    Nt=Na+Nb\n",
    "\n",
    "    Gm=np.eye(Nt)-2*GA\n",
    "    Sp=np.array([[0,1],[0,0]])\n",
    "    Gm_1=block_diag(-Gm[:Na,:Na],Gm[Na:,Na:])\\\n",
    "        +np.kron(Sp,-1j*Gm[:Na,Na:])+np.kron(Sp.T,-1j*Gm[Na:,:Na])\n",
    "    Gm_2=block_diag(-Gm[:Na,:Na],Gm[Na:,Na:])\\\n",
    "        +np.kron(Sp,1j*Gm[:Na,Na:])+np.kron(Sp.T,1j*Gm[Na:,:Na])    \n",
    "    Gx=np.dot(np.eye(Nt)-Gm_2,np.dot(np.linalg.inv(np.eye(Nt)+np.dot(Gm_1,Gm_2)),np.eye(Nt)-Gm_1))/2\n",
    "    Gxh=(Gx+np.transpose(Gx).conj())/2\n",
    "\n",
    "    chi0 =np.linalg.eigvalsh(Gxh)\n",
    "    i0=np.nonzero(np.abs(chi0)>epsilon)[0]\n",
    "    chi1=chi0[i0]\n",
    "    i1=np.nonzero(np.abs(chi1-1)>epsilon)[0]\n",
    "    nu=chi1[i1]\n",
    "    eA=np.sum(np.log(nu**0.5+(1-nu)**0.5))\n",
    "    \n",
    "    chi,VR=np.linalg.eigh(GA)\n",
    "    sA=np.sum(np.log(chi**2+(1-chi)**2))/2    \n",
    "\n",
    "    return eA+sA\n",
    "\n",
    "def exact_Renyi_neg_calc2(r,Na,Nb,GA,epsilon=1e-9):\n",
    "    \"\"\"computes log negativity of 1d free fermions only works for even r\"\"\"\n",
    "    Nt=Na+Nb\n",
    "\n",
    "    Gm=np.eye(Nt)-2*GA\n",
    "    Sp=np.array([[0,1],[0,0]])\n",
    "    Gm_1=block_diag(-Gm[:Na,:Na],Gm[Na:,Na:])\\\n",
    "        +np.kron(Sp,-1j*Gm[:Na,Na:])+np.kron(Sp.T,-1j*Gm[Na:,:Na])\n",
    "    Gm_2=block_diag(-Gm[:Na,:Na],Gm[Na:,Na:])\\\n",
    "        +np.kron(Sp,1j*Gm[:Na,Na:])+np.kron(Sp.T,1j*Gm[Na:,:Na])    \n",
    "    Gx=np.dot(np.eye(Nt)-Gm_2,np.dot(np.linalg.inv(np.eye(Nt)+np.dot(Gm_1,Gm_2)),np.eye(Nt)-Gm_1))/2\n",
    "    Gxh=(Gx+np.transpose(Gx).conj())/2\n",
    "\n",
    "    chi0 =np.linalg.eigvalsh(Gxh)\n",
    "    i0=np.nonzero(np.abs(chi0)>epsilon)[0]\n",
    "    chi1=chi0[i0]\n",
    "    i1=np.nonzero(np.abs(chi1-1)>epsilon)[0]\n",
    "    nu=chi1[i1]\n",
    "    eA=np.sum(np.log(nu**(r/2)+(1-nu)**(r/2)))\n",
    "    \n",
    "    chi,VR=np.linalg.eigh(GA)\n",
    "    sA=r*np.sum(np.log(chi**2+(1-chi)**2))/2    \n",
    "\n",
    "    return eA+sA\n",
    "\n",
    "def exact_Renyi_neg_calc(r,N,Na,Nb,V1):\n",
    "    '''computes moments of partial transpose tr(\\rho^T \\rho^T\\dag ... )\n",
    "    for 1d free fermions ''' \n",
    "    k_sw=np.arange((-(r-1)/2),(r-1)/2+0.1,1)\n",
    "    Zk=0\n",
    "    if (r % 2) == 0:\n",
    "        delta=pi\n",
    "    else:\n",
    "        delta= pi*(r-1)/r\n",
    "#     delta=pi\n",
    "\n",
    "    Nflip=Na+Nb\n",
    "    N1=int((N-Nflip)/2)\n",
    "\n",
    "    for i_k in range(len(k_sw)):\n",
    "#         op_k=block_diag(np.eye(N1),np.exp(1j*2*pi/r*k_sw[i_k])*np.eye(Na),np.exp(1j*delta-1j*2*pi/r*k_sw[i_k])*np.eye(Nb),np.eye(N1))\n",
    "        op_k=block_diag(np.exp(1j*2*pi/r*k_sw[i_k])*np.eye(Na),np.exp(1j*delta-1j*2*pi/r*k_sw[i_k])*np.eye(Nb),np.eye(N-Nflip))\n",
    "        Zk +=np.real(np.log( np.linalg.det( np.dot(np.dot(np.matrix(V1).H,op_k),V1) )))\n",
    "    return Zk\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "N=100 # system size\n",
    "Lsub_list=np.arange(0,int(N/2))\n",
    "# Lsub_list=np.arange(2,7,1)\n",
    "# Lsub_list=[6]\n",
    "N_pt = int(N/2)\n",
    "\n",
    "# reference slater determinant\n",
    "t= 0\n",
    "# hopping amplitudes\n",
    "t1= 1-t\n",
    "t2= 1+t\n",
    "BC=np.exp(1j*pi) # boundary condition on a chain, you can put BC=0 for open chain\n",
    "# BC=1 periodic boundary condition and BC=-1 is anti-periodic\n",
    "# do not put BC=1 since the gs is not unique in that case\n",
    "V1=wf_gen(N,N_pt,BC,t1,t2) # eigenvectors in \n",
    "\n",
    "\n",
    "def ent_plotter(r):\n",
    "    '''r=2 # Renyi/ Negativity index'''\n",
    "\n",
    "    Rr_ex=np.zeros(len(Lsub_list)) # Renyi entropy\n",
    "    RNr_ex=np.zeros(len(Lsub_list)) # Renyi negativity\n",
    "    RNr2_ex=np.zeros(len(Lsub_list)) # Renyi negativity 2\n",
    "    LN_ex=np.zeros(len(Lsub_list)) # Log negativity\n",
    "\n",
    "    Gmat=np.dot(V1,np.matrix(V1).H)\n",
    "\n",
    "    for i_L in range(len(Lsub_list)):\n",
    "        Lsub=Lsub_list[i_L]\n",
    "        inds_A= np.arange(0,Lsub_list[i_L])\n",
    "        inds_AB= np.arange(0,int(2*Lsub_list[i_L]))\n",
    "\n",
    "        Rr_ex[i_L]=exact_renyi_calc(r,Gmat[np.ix_(inds_AB,inds_AB)])\n",
    "        RNr_ex[i_L]=exact_Renyi_neg_calc(r,N,Lsub,Lsub,V1)\n",
    "        RNr2_ex[i_L]=exact_Renyi_neg_calc2(r,Lsub,Lsub,Gmat[np.ix_(inds_AB,inds_AB)])\n",
    "    plt.figure()\n",
    "#     plt.plot(Lsub_list,-Rr_ex,'o',label='(-1) Renyi Entropy')\n",
    "    plt.plot(Lsub_list,(RNr_ex),'v',label='Renyi Neg')\n",
    "    if r % 2 ==0:\n",
    "        plt.plot(Lsub_list,(RNr2_ex),'x',label='Renyi Neg 2')\n",
    "    plt.title('r=%d Renyi index' % (r))\n",
    "    plt.xlabel('$\\ell$')\n",
    "    plt.legend()\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEYCAYAAABfgk2GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAi2UlEQVR4nO3de5hU1Z3u8e/LRdGA2ERQLiJtbBPFCEKbiPE2CpEgBNQQEpWAJofkxCgkJtEEz0xmHpnHjE6CRGYcRgwmesQrI0EmGdBh0KNRG0ISEB2NeGnsADKKGjWi/M4fVZCmqeru6l3d1VX7/TxPPb1r16paa3d3vbVq7b3XVkRgZmaVr0upG2BmZh3DgW9mlhIOfDOzlHDgm5mlhAPfzCwlHPhmZinhwDdrJ5K+L+nmPI+dKumZBK8dko5qe+ssjRz4VnYk7S9poaQXJb0p6TeSPtNM+emSPpD0lqQ3JP1W0vj2bmdE/H1EfCXPYw9HxEfbuw1mjTnwrVOQ1K2A4t2Al4HTgd7A/wHukjSkmec8FhE9gYOBfwIWSzq4TY01K1MOfCsZSS9IulLS74A/tTb0I+JPEfGDiHghInZFxDJgEzCyFc/dBfwc+BBQk23H/pKul/SSpC2SbpJ0QPaxMyTVS7pC0lZJDZIuzj52Yrb8nnZLOl/SuuzyDyTdlmfbz5BU3+R38W1Jv5O0Q9Kdkno0evw72bpfkXRJk9dqrv3LJf1jo7J3Srqlpd+TVSYHvpXaF4FzgIMj4n1JyyS9nue2LNcLSDoUOBrY0FJlkroCFwM7gRezq3+Yff5w4ChgIPDXjZ52GJlvEgOBLwPzJVVFxJPAdmBMo7IXkflAaYvPA2OBauB4YHq2zWOBb2frqQFGN3lec+2/BJgq6UxJFwInAjPb2D4rdxHhm28luQEvAJckfI3uwErgX5opMx14H3idTNC/A3w++5iAPwEfaVR+FLApu3xGtny3Ro9vBU7KLl8J3J5d7gO8DfTP3v8BcFueNp0B1Df5XVzU6P4/ADdll28Brm302NFAkAn3ZtufvX8emSGwV4FTSv139610N/fwrdRebusTJXUh05t+D/hGC8V/HREHA1XAUuDU7Pq+wIHAmt3fJIBfZtfvtj0i3m90/22gZ3b5NmCCpJ5keugPR0RDGzfpj3nqGMDev6cXGy23pv3LgK7AMxHxSBvbZhXAgW+lttd0rZL+PXs0Ta7bvzcqJ2AhcChwfkTsbFVlEW8BXyczzHECmV7vO8DQiDg4e+sdmR28rXm9zcBjwLnAVNo+nNOcBuDwRvcHN1puTfvnABuB/pK+2A7tszLhwLdOJSI+ExE989waH3r5z8AxwISIeKfAOrYDNwN/HZmduP8K/FhSPwBJAyWdXcBL/gz4LvBxYEkhbWmlu4Dpko6VdCDwN7sfaKn9kk4js8/iS9nbTyQNbIc2Whlw4FvZkXQE8FUyOyn/2OgbwIUFvMxcYJyk48mMwz8H/FrSG2T2CRRyjPwS4AhgSUT8qYDntUpE/DuZ9j6UbedDTYrkbL+kg8h8GH0jIjZnh3MWAj/NfkOylFGEL4BilpSkPwBfjYiVpW6LWT7u4ZslJOl8Mvsimva8zTqVQs5uNLMmJK0CjgWmZsfTzTotD+mYmaWEh3TMzFKiUw/pHHLIITFkyJBSN8PMrGysWbPm1Yjom+uxTh34Q4YMoa6urtTNMDMrG5JezPeYh3TMzFLCgW9mlhIOfDOzlOjUY/hm1rnt3LmT+vp63n333VI3JXV69OjBoEGD6N69e6uf48A3szarr6+nV69eDBkyBE/P03Eigu3bt1NfX091dXWrn1eUwM9ekecGMnNu3xwR1zZ5XNnHx5GZ53t6RKwtRt2N3fLDmazYMZDHdg3ds25Ulw2M6b2ZS668odjVmaXeu+++67AvAUl8+MMfZtu2bQU9L/EYfvaScfOBz5A5xfyLko5tUuwzZC7NVgPMIDO1bdHFgBHM7z6PUV0yV7ob1WUD87vPIwaMaI/qzAwc9iXSlt97MXbafgJ4LiKej4j3gMXAxCZlJgI/i4xfAwdL6l+EuvcyYeIUZu2axY3d5/HNbndnfu6axYRJU4pdlZlZ2SlG4A9k78uv1WfXFVoGAEkzJNVJqiv060q/g3oweMTZ3LFrDDO7LeGOXWM4fORY+vXqUdDrmFn56Nq1K8OHD+e4445jwoQJvP7660V9/bq6Oi6//PJ91q9atQpJ/OIXv9izbvz48axataqo9RdTMcbwc32vaDojW2vKZFZGLAAWANTW1hY8s9sVNVuIdSu44f1zmdp1Baq5EDiu0JcxsyIbd8PDPNXwxj7rj+1/EMtnnprjGa1zwAEHsG7dOgCmTZvG/PnzmT17dptfr6na2lpqa2tzPjZo0CDmzJnDhAkTilZfeypGD7+eva+3OQh4pQ1lktu0mqrlM1hy1BzmfjCZJUfNoeqBGbBpddGrMrPCjBh8MN277t33695VjDiiqmh1jBo1is2bNwPwhz/8gbFjxzJy5EhOPfVUnn76aQCmT5/O5Zdfzsknn8yRRx7JPffcA8DUqVO5//7797zWhRdeyNKlS1m1ahXjx4/PWd+wYcPo3bs3K1as2OexNWvWcPrppzNy5EjOPvtsGhoy17Z/8sknOf744xk1ahTf+c53OO64juuQFiPwnwRqJFVL2g/4ArC0SZmlwJeUcRKwIyIailD33javhcmLmDBxCicO6ZMZu5+8KLPezErq8rNq6NJkR2NXicvPOqoor//BBx/w4IMP8tnPfhaAGTNm8JOf/IQ1a9Zw/fXX8/Wvf31P2YaGBh555BGWLVvGVVddBcBXvvIVfvrTnwKwY8cOHn30UcaNG9divVdffTXXXHPNXut27tzJZZddxj333MOaNWu45JJL9nzruPjii7npppt47LHH6Nq1a1G2vbUSD+lExPuSvgH8isxhmbdExAZJX8s+fhOwnMwhmc+ROSzz4qT15nTKLAD6AXd9dVRmXa/ToPq0dqnOzFqv30E9mDxyEHfWvczOD4LuXcXnag9PvI/tnXfeYfjw4bzwwguMHDmSMWPG8NZbb/Hoo48yefLkPeX+/Oc/71meNGkSXbp04dhjj2XLli0AnH766Vx66aVs3bqV++67j/PPP59u3VqOyFNPzQxHPfzww3vWPfPMM6xfv54xY8YAmQ+j/v378/rrr/Pmm29y8sknA3DBBRewbNmyRNtfiKIchx8Ry8mEeuN1NzVaDuDSYtRlZuXr8rNquHtNPRBF693vHsPfsWMH48ePZ/78+UyfPp2DDz54z9h+U/vvv/+e5cYXgZo6dSq33347ixcv5pZbbml1G2bPns2cOXP2fEBEBEOHDuWxxx7bq9xrr71WwJYVn+fSMbMOs7uXL1GU3n1jvXv3Zt68eVx//fUccMABVFdXc/fddwOZAP7tb3/b4mtMnz6duXPnAjB06NDmCzfy6U9/mtdee21PHR/96EfZtm3bnsDfuXMnGzZsoKqqil69evHrX/8agMWLFxeyiYk58M2sQ11+Vg0nDulTtLH7xk444QSGDRvG4sWLuf3221m4cCHDhg1j6NChe+2QzefQQw/lmGOO4eKLCx91nj17NvX19QDst99+3HPPPVx55ZUMGzaM4cOH8+ijjwKwcOFCZsyYwahRo4gIevfuXXBdbdWpr2lbW1sbvgCKWee1ceNGjjnmmFI3o2jefvttPv7xj7N27dp2C+K33nqLnj17AnDttdfS0NDADTe0beqXXL9/SWsiIudxpO7hm5kBK1eu5GMf+xiXXXZZu/a6H3jggT0nij388MNcffXV7VZXU54t08wMGD16NC+99FK71zNlyhSmTCnNdC/u4ZuZpYQD38wsJRz4ZmYp4cA3M0sJB76ZlbVKnB553bp1jBo1iqFDh3L88cdz5513Jn5NcOCbWUd5ZO6+M9duWp1Zn8DuqRXWr19Pnz59mD9/fqLXa6q2tpZ58+blfGz39MjFduCBB/Kzn/2MDRs28Mtf/pJZs2YV5YPMgW9mHWPgCLh7+l9Cf9PqzP2BxbsEaaVMj3z00UdTU1MDwIABA+jXr1/B16/NxYFvZh2j+rTMdOV3T4eH5mR+Tl5UtNlsK3V65CeeeIL33nuPj3zkI637RTTDJ16ZWcepPg1qvwyr/wFO+25Rwr6Sp0duaGhg6tSp3HrrrXTpkrx/7sA3s46zaTXULcyEfd1CqD41cehX6vTIb7zxBueccw7XXHMNJ510Uquf1xwP6ZhZx9g9Zj95EZw5+y/DO0W6BGklTY/83nvvce655/KlL31pr28pSTnwzaxjZC9BuqdHv3tMv4iXIK2U6ZHvuusuVq9ezaJFixg+fDjDhw/P+22lEJ4e2czazNMjF87TI5uZlZinR26BpD7AncAQ4AXg8xHxWpMyhwM/Aw4DdgELIqJtH2dmZu3E0yO37CrgwYioAR7M3m/qfeCKiDgGOAm4VNKxCes1s06iMw8LV7K2/N6TBv5E4Nbs8q3ApKYFIqIhItZml98ENgIDE9ZrZp1Ajx492L59u0O/g0UE27dvp0ePwi4Cn/Q4/EMjoiHbgAZJ/ZorLGkIcALweDNlZgAzAAYPHpyweWbWngYNGkR9fX1RTvu3wvTo0YNBgwYV9JwWA1/SSjLj703NLqQiST2Be4FZEfFGvnIRsQBYAJmjdAqpw8w6Vvfu3amuri51M6yVWgz8iBid7zFJWyT1z/bu+wNb85TrTibsb4+I+9rcWjMza7OkY/hLgWnZ5WnAPmc2SBKwENgYET9KWJ+ZmbVR0sC/Fhgj6VlgTPY+kgZIWp4t8ylgKnCmpHXZW8tT0JmZWVEl2mkbEduBs3KsfwUYl11+BFCSeszMLDmfaWtmlhIOfDOzlHDgm5mlhAPfzCwlHPhmZinhwDczSwkHvplZSjjwzcxSwoFvZpYSDnwzs5Rw4JuZpYQD38wsJRz4ZmYp4cA3M0sJB76ZWUo48M3MUsKBb2aWEg58M7OUSBT4kvpIWiHp2ezPqmbKdpX0G0nLktRpZmZtk7SHfxXwYETUAA9m7+czE9iYsD4zM2ujpIE/Ebg1u3wrMClXIUmDgHOAmxPWZ2ZmbZQ08A+NiAaA7M9+ecrNBb4L7GrpBSXNkFQnqW7btm0Jm2dmZrt1a6mApJXAYTkemt2aCiSNB7ZGxBpJZ7RUPiIWAAsAamtrozV1mJlZy1oM/IgYne8xSVsk9Y+IBkn9ga05in0K+KykcUAP4CBJt0XERW1utZmZFSzpkM5SYFp2eRpwf9MCEfG9iBgUEUOALwAPOezNzDpe0sC/Fhgj6VlgTPY+kgZIWp60cWZmVjwtDuk0JyK2A2flWP8KMC7H+lXAqiR1mplZ2/hMWzOzlHDgm5mlhAPfzCwlHPhmZinhwDczSwkHvplZSjjwzcxSwoFvZpYSDnwzs5Rw4JuZpYQD38wsJRz4ZmYp4cA3M0sJB76ZWUo48M3MUsKBb2aWEg58M7OUcOCbmaWEA9/MLCUSBb6kPpJWSHo2+7MqT7mDJd0j6WlJGyWNSlKvmZkVLmkP/yrgwYioAR7M3s/lBuCXEfExYBiwMWG9ZmZWoKSBPxG4Nbt8KzCpaQFJBwGnAQsBIuK9iHg9Yb1mZlagpIF/aEQ0AGR/9stR5khgG/BTSb+RdLOkD+V7QUkzJNVJqtu2bVvC5pmZ2W4tBr6klZLW57hNbGUd3YARwD9HxAnAn8g/9ENELIiI2oio7du3byurMDOzlnRrqUBEjM73mKQtkvpHRIOk/sDWHMXqgfqIeDx7/x6aCXwzM2sfSYd0lgLTssvTgPubFoiIPwIvS/podtVZwFMJ6zUzswIlDfxrgTGSngXGZO8jaYCk5Y3KXQbcLul3wHDg7xPWa2ZmBWpxSKc5EbGdTI+96fpXgHGN7q8DapPUZWZmyfhMWzOzlHDgm5mlhAPfzCwlHPhmZinhwDczS4lER+mUi1t+OJMVOwby2K6he9aN6rKBMb03c8mVN5SwZWZmGR2RU6no4ceAEczvPo9RXTYAmV/i/O7ziAEjStwyM7OMjsipVPTwJ0ycwqzrtnFj97nc9sFoLuq6km/umsV1k6aUumlmZkDH5FQqevj9DurB4BFnc8euMczstoQ7do3h8JFj6derR6mbZmYGdExOpaKHD3BFzRZi3QpueP9cpnZdgWouBI4rdbPMzPZo75xKRQ+fTaupWj6DJUfNYe4Hk1ly1ByqHpgBm1aXumVmZhkdkFPp6OFvXguTFzHhw5/gV2//hgmTzoJXqzPrq08rdevMzDokpxQRRXmh9lBbWxt1dXWlboaZWdmQtCYick5WmY4hHTMzc+CbmaWFA9/MLCUc+GZmKeHANzNLiUSBL6mPpBWSns3+rMpT7puSNkhaL+kOST7F1cysgyXt4V8FPBgRNcCD2ft7kTQQuByojYjjgK7AFxLWa2ZmBUoa+BOBW7PLtwKT8pTrBhwgqRtwIPBKwnrNzKxASQP/0IhoAMj+7Ne0QERsBq4HXgIagB0R8R/5XlDSDEl1kuq2bduWsHlmZrZbi4EvaWV27L3pbWJrKsiO608EqoEBwIckXZSvfEQsiIjaiKjt27dva7fDzMxa0OJcOhExOt9jkrZI6h8RDZL6A1tzFBsNbIqIbdnn3AecDNzWxjabmVkbJB3SWQpMyy5PA+7PUeYl4CRJB0oScBawMWG9ZmZWoKSBfy0wRtKzwJjsfSQNkLQcICIeB+4B1gK/z9a5IGG9ZmZWIM+WaWZWQTxbppmZOfDNzNLCgW9mlhIOfDOzlHDgm5mlRDouYp7HLT+cyYodA3ls19A960Z12cCY3pu55MobStgyM6tUpcydVPfwY8AI5nefx6guG4DML31+93nEgBElbpmZVapS5k6qe/gTJk5h1nXbuLH7XG77YDQXdV3JN3fN4rpJU0rdNDOrUKXMnVT38Psd1IPBI87mjl1jmNltCXfsGsPhI8fSr5evz2Jm7aOUuZPqHj7AFTVbiHUruOH9c5nadQWquRA4rtTNMrMKVqrcSXUPn02rqVo+gyVHzWHuB5NZctQcqh6YAZtWl7plZlapSpg76e7hb14Lkxcx4cOf4Fdv/4YJk86CV6sz66tPK3XrzKwSlTB3PHmamVkF8eRpZmbmwDczSwsHvplZSjjwzcxSwoFvZpYSiQJf0mRJGyTtkpRzr3C23FhJz0h6TtJVSeo0M7O2SXoc/nrgPOBf8hWQ1BWYT+Yi5/XAk5KWRsRTCetuN55F08yS6ow5kqiHHxEbI+KZFop9AnguIp6PiPeAxcDEJPW2N8+iaWZJdcYc6YgzbQcCLze6Xw98sgPqbTPPomlmSXXGHGmxhy9ppaT1OW6t7aUrx7q8p/dKmiGpTlLdtm3bWllFcXkWTTNLqjPmSIs9/IgYnbCOeuDwRvcHAa80U98CYAFkplZIWHebeRZNM0uqs+VIRxyW+SRQI6la0n7AF4ClHVBv23kWTTNLqhPmSKIxfEnnAj8B+gIPSFoXEWdLGgDcHBHjIuJ9Sd8AfgV0BW6JiA2JW96ePIummSXVCXPEs2WamVUQz5ZpZmYOfDOztHDgm5mlRLovcdgGnfF0aTMrnXLKBPfwC9QZT5c2s9Ipp0xwD79AnfF0aTMrnXLKBPfwC9QZT5c2s9Ipp0xwD78NOtvp0mZWWuWSCe7hF6oTni5tZiVURpngHn6hOuHp0mZWQmWUCZ5awcysgjQ3tYJ7+EVSTsfimllhKuX97TH8IimnY3HNrDCV8v52D79IyulYXDMrTKW8v93DL5JyOhbXzApTKe9v9/CLqFyOxTWzwlXC+9s9/GIpo2NxzaxAFfL+dg+/WPIci/vk7T/gR2+PLfu9+2ZpkO9onG8d+EtOvLA8jrVvjnv4xXLKLKg+jX4H9eCur47KjO1Vn8bvjphWEXv3zdIg39E4vztiWs73N6fMKm2DC5Qo8CVNlrRB0i5JOQ/0l3S4pP+UtDFbdmaSOsvNhIlTmLVrFjd2n8c3u92d+blrFhPKbO++WRpU+vs1aQ9/PXAe0NxA1vvAFRFxDHAScKmkYxPWWzYqZe++WRpU+vs10Rh+RGwEkNRcmQagIbv8pqSNwEDgqSR1l5NK2LtvlhaV/H7t0J22koYAJwCPN1NmBjADYPDgwR3TsPaU3bu/8Kg5zH3qEHp+9K/48gMzWH7/ifz8rZHemWtWAvl2zk7tWce4bnX7vF85cFFZ7ZzNp8UhHUkrJa3PcZtYSEWSegL3ArMi4o185SJiQUTURkRt3759C6mic9p99M7EKZw4pE9mLHDyIg7puZ935pqVSL6ds4f03D/n+5XNa0vb4CIpymyZklYB346InFNbSuoOLAN+FRE/au3rVvJsmVvfeJdvX3cjP+7yl1O1v7VrFtd99xsVM15o1llV8vuvudky2/2wTGUG+BcCGwsJ+0pX6TuHzDqztL7/Eo3hSzoX+AnQF3hA0rqIOFvSAODmiBgHfAqYCvxe0rrsU78fEcuT1F0Jcu0cuuWH/1oR07CadQbNTWt8xTkTKnbnbD6JevgRsSQiBkXE/hFxaEScnV3/SjbsiYhHIkIRcXxEDM/eUh/2+U7VPqyqp8f2zYok31j9YVU9K2KqhEJ5aoVSyTMVw2nPP8HXX55V9tOwmnUG+aY1nv8R4MjynyqhUL7EYSd09ZLfc9jaH/GNrvdx4wfn8ceRV3DNpMr+qmnWXtL2fvIlDstMvhM/KuUya2bF1NL7opJPpCqUJ0/rbJqZhrVSLrNmVkzNvi8qZFrjYnEPv7PJM7bP5rVMmPi1irjMmlkxNXv5wd/elPf9VMlj9fl4DL/M5BqPPPIZH8ppla+5oZvnj/5Kqsbpm1PSE6+suK6o2cIFXTLjkRd0WcEVNX/0UI+lQnP/57neF7YvB345yTMeed4JAyt6Dm8zyD9X/XkjBnqcvpU8hl9O8ozvV21emzlNfO0GZnbLfKU9fORYlv3TlR7qsbLT3NDN4BFf2ef/vOq1lR6nbyX38MtJnssocsosD/VYxSh46KaZ94XtzYFfCTzUYxXEQzftx0M6lcBDPVaGPHTT8dzDrwQe6rEy5KGbjucefiXLc3nF88YtYOZ/556gzdM3WLHl+586s1d9duhm7//DuSMGUvXAvv+3lXKZwVJyD7+S5bm8YtXr6/Ne/MG9fyu2fP9TXQaNzPl/WPXa+oq+zGAp+UzblHpt/Uri7un8/IPRTO26En1+EVVDR+e99NspH3qJh94Y5J6/5dTcN8Px//uHeS8n2P3FR3L+H1rb+Uxb21szE0rlu/Rbl4Ej3fO3vJr7Zpj3coKvPuGjbjqYx/DTqJkJ2qg+Led0sjsPzz1B1Smb1/DF729xzz8lCh2P3z2xX84pijev91E3HSzRkI6kycAPgGOAT0RE3vEXSV2BOmBzRIxvzet7SKcENq2Gu6ezsP/fcM1Th3D1sa/y5Ya/hcmLuHpd1T4TVFW/uYZzn5vNpTsv57FdQ/f07B7SJ7n3vU/6g6BM5Qv28/d7nDPj8X3+3kuOmsOmniNyT2DWzP+Ug7342vMCKOuB84B/aUXZmcBG4KCEdVp7aqb3f0XNca3u+X/yyD7Mf+Gv9/0geOeTfPH7/+APgk4if7DXM7/7vfv8/RYf8XfMen5Uzp589xcfyX2hkRa+UVrHSRT4EbERQFKz5SQNAs4B5gDfSlKntbPsMc39gLu+OiqzrlfmTVl19/R9D5WbvCjnyV2fO/MoZl33P63+IFgyYI4PCW1HRQv2z13AKyuf3efv3e/VJyDHIcAcuCj//5TDvsN11Bj+XOC7QK8Oqs+KrcCef9VBxxX0QXDdpCn84t9g/ts5hof8raBVmvvAjAEjcv5uCw72Xj08Hl/GWgx8SSuBw3I8NDsi7m/F88cDWyNijaQzWlF+BjADYPDgwS0Vt45SaM//wEVcUfNeqz8I+vXqkffKRR4e2luhvfUlA+bkvypUocG+6X9ynsy3ezzevfjOrcXAj4ikB8V+CvispHFAD+AgSbdFxEV56lsALIDMTtuEdVt7y9fz//29VD39i1Z/EMBxfzl8L+HwUL4Pgq91e4Cb3j+n031A5AvwfO1ty/h6v165f7cFB/sxE9yTL2PtPqQTEd8DvgeQ7eF/O1/YWxnK1/NvwwdBvkNCCx0eyvdBcGePycx/d15B+w/yhW6h68f03gxQUM88X3vbMgwDeQ6NbGOwuydfnhKdeCXpXEn1wCjgAUm/yq4fIGl5MRpoZSrfRFd9qvOfNt/MCWG5JtPKd0LP586/IOf0uqdP+0He6aLznTj01AHFWR8DRuStY+sR5xTU3s997oK8U2PkvdRfvt/t+ntz/z2qqj1RWQVKepTOEmBJjvWvAONyrF8FrEpSp5W55o7YeGRuUYaH8n0jOHZA74L3H1w57X8x658HJl5/3aQpEBQ0jt5cewvqrR+4qMVDI91jTwefaWudRxGHh/LtJ8gZlM3sP8gXuoWu3z2sUlCA52tvW4ZhfGik4cC3clCsD4Id36LqkR8VtP+guQ+JQtdDgePo+drr8XVrIwe+la9CPwj+37z8QQmFhW6h6w9clL+OfAHeXHsd7NYWEdFpbyNHjgyzDvHwjyOe/6/YsuOdmHzTo7HljXcinv+viJ+fX5z1D/84fx0P/7i0224VBaiLPJnq+fDNzCqI58M3MzMHvplZWjjwzcxSwoFvZpYSDnwzs5To1EfpSNoGvNjGpx8CvFrE5pQLb3e6eLvTpTXbfURE9M31QKcO/CQk1eU7NKmSebvTxdudLkm320M6ZmYp4cA3M0uJSg78BaVuQIl4u9PF250uiba7Ysfwzcxsb5Xcwzczs0Yc+GZmKVFxgS9prKRnJD0n6apSt6c9SbpF0lZJ6xut6yNphaRnsz+rStnGYpN0uKT/lLRR0gZJM7PrK327e0h6QtJvs9v9t9n1Fb3du0nqKuk3kpZl76dlu1+Q9HtJ6yTVZde1edsrKvAldQXmA58BjgW+KOnY0raqXS0CxjZZdxXwYETUAA9m71eS94ErIuIY4CTg0uzfuNK3+8/AmRExDBgOjJV0EpW/3bvNBDY2up+W7Qb4q4gY3uj4+zZve0UFPvAJ4LmIeD4i3gMWAxNL3KZ2ExGrgf9psnoicGt2+VZgUke2qb1FRENErM0uv0kmBAZS+dsdEfFW9m737C2o8O0GkDQIOAe4udHqit/uZrR52yst8AcCLze6X59dlyaHRkQDZMKRzBUAK5KkIcAJwOOkYLuzwxrrgK3AiohIxXYDc4HvArsarUvDdkPmQ/0/JK2RNCO7rs3bXmnXtFWOdT7utAJJ6gncC8yKiDekXH/6yhIRHwDDJR0MLJF0XImb1O4kjQe2RsQaSWeUuDml8KmIeEVSP2CFpKeTvFil9fDrgcMb3R8EvFKitpTKFkn9AbI/t5a4PUUnqTuZsL89Iu7Lrq747d4tIl4HVpHZf1Pp2/0p4LOSXiAzRHumpNuo/O0GICJeyf7cCiwhM2zd5m2vtMB/EqiRVC1pP+ALwNISt6mjLQWmZZenAfeXsC1Fp0xXfiGwMSJ+1OihSt/uvtmePZIOAEYDT1Ph2x0R34uIQRExhMz7+aGIuIgK324ASR+S1Gv3MvBpYD0Jtr3izrSVNI7MmF9X4JaImFPaFrUfSXcAZ5CZMnUL8DfAvwF3AYOBl4DJEdF0x27ZknQK8DDwe/4ypvt9MuP4lbzdx5PZQdeVTEftroj4O0kfpoK3u7HskM63I2J8GrZb0pFkevWQGX7/vxExJ8m2V1zgm5lZbpU2pGNmZnk48M3MUsKBb2aWEg58M7OUcOCbmaWEA9/MLCUc+GZmKeHANytQ9gzIG7PTE5uVDQe+WeG+BuwPnFLqhpgVwoFvVrixwH8D60rcDrOCOPDNCiCpB5n5bEYA/1Xi5pgVxIFvVpgaMoH/dETsLHVjzApRaRdAMWtvfYGjqeBLZ1rlcg/frDADyFx8pYukqlI3xqwQDnyzVpLUjczY/WHATcAHpW2RWWE8H76ZWUq4h29mlhIOfDOzlHDgm5mlhAPfzCwlHPhmZinhwDczSwkHvplZSvx/RlUpgc/BzQEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEYCAYAAABfgk2GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAesklEQVR4nO3df5RcZZ3n8feHJhCEkIASCDQxcWiBBEkkDWtAlh8JGDKBqIiAigTWE12FMOeMDnHC0dkZM8uu7K6omWWygOOscSOwZolJFBLcLOEAC90uaJCAiBGa9JKY5Vfkh0n47h91O1Q6Vd1161Z1ddX9vM6pk6pbT93neaor33rqe5/7XEUEZmbW+vZrdAPMzGxoOOCbmeWEA76ZWU444JuZ5YQDvplZTjjgm5nlhAO+WZ1I+mtJt5Z57kxJT2XYd0g6rvrWWR454FtTkvQDSb2SXpX0tKTPDVB2nqTdknYk5R+XNKfebYyIv4+Iku2KiA0RcXy922BWzAHfhgVJ+6d8yb8FJkTEocBFwDckTRug/EMRcQgwBvgHYLmkMdW01axZOeBbw0jaLOl6Sb8E/pgm6EfEExHxVt/D5PZnFbzubeC/AgcDHUk7DpR0k6TnJL0o6RZJByXPnS2pR9JfStqa/Kq4Knnu1KT8nnZLuljSY8n9v5H0gzJ9P1tST7/34suSfinpFUk/kjSy6PmvJHVvkXR1v30N1P41kv5DUdkfSbp9sPfJWpMDvjXa5cCfA2MiYpekVZJeLnNbVfxCSf8g6XVgE9ALrBmsMkltwFXATuD3yeZ/B7wfmAocBxwDfK3oZUcBo5Pt/wpYIumwiHgU2A6cV1T2MxS+UKrxSWAWMBE4GZiXtHkW8OWkng5gZr/XDdT+q4ErJJ0r6dPAqcB1VbbPml1E+OZbQ27AZuDqjPtoAz4M3ACMKFNmHrALeJlCoH8D+GTynIA/An9WVH468Lvk/tlJ+f2Lnt8KfCi5fz2wLLl/OPA6MC55/DfAD8q06Wygp9978Zmix/8euCW5fztwY9Fz76fwi+a4wdqfPP448DzwB+DDjf67+9a4m0f41mjPZ3lxROyOiAeAduBfD1D04YgYAxwGrATOTLYfAbwL6O77JQH8LNneZ3tE7Cp6/DpwSHL/B8CFkg6hMELfEBG9VXbn/5ap42j2fp9+X3S/kvavovDF+FTyXllOOeBbo+21XKuknyazaUrdfjrAfvanshz+DuCLFNIcH6Qw6n0DmBwRY5Lb6Cgc4B288REvAA8BHwOuoPp0zkB6gWOLHo8vul9J+xcDTwLjJF1eh/ZZk3DAt2ElIi6IiEPK3C4AkDRW0mWSDpHUJukjFI4F/LzCOrYDtwJfi8JB3P8C/CdJY5P9H5Pss1L/DPwV8AFgRYrXVeoOYJ6kSZLeBXy974nB2i/pX1I4ZvHZ5PYdScfUoY3WBBzwrRkFhfRND/AScBPwFxFxd4p9fAuYLelkCnn4Z4CHJb0KrAPSzJFfAbwXWBERf0zxuopExE8ptPfnSTv7f7GVbL+kQyl8GV0TES8k6ZzbgO9JUq3bacOfInwBFLOsJP0W+HxErGt0W8zK8QjfLCNJF1P41VFRSsmsUdKe3WhmRSStByYBVyT5dLNhyykdM7OccErHzCwnhnVK5z3veU9MmDCh0c0wM2sa3d3df4iII0o9N6wD/oQJE+jq6mp0M8zMmoak35d7zikdM7OccMA3M8sJB3wzs5wY1jl8M2sNO3fupKenhzfffLPRTWkZI0eOpL29nREjRlT8Ggd8M6u7np4eRo0axYQJE/AyPtlFBNu3b6enp4eJEydW/LqaBPzkijw3U1hz+9aIuLHf80qen01hne95EfGLWtRdbPbNG/h176v7bJ807lDWXHdmiVeY2VB48803HexrSBLvfve72bZtW6rXZc7hJ5eMWwJcQOEU88slTepX7AIKl2brAOYD/zlrvaWcMn4MI9r2/kCNaBOnvPewelRnZik42NdWNe9nLQ7angY8ExHPRsSfgOXA3H5l5gL/HAUPA2MkjatB3XtZMKOD/fq9CW0SC2YcV+uqzMyaTi0C/jHsffm1nmRb2jIASJovqUtSV9qfK2MPHckl09r3jPJHtIlPdB7L2FEjU+3HzFpPW1sbU6dO5aSTTuLCCy/k5Zdfrun+u7q6WLBgwT7b169fjyR+8pOf7Nk2Z84c1q9fX9P6K1GLHH6p3xX9V2SrpExhY8RSYClAZ2dn6pXdFszo4M7uHiA8ujdrQvU6FnfQQQfx2GOPAXDllVeyZMkSFi1aVPX++uvs7KSzs7Pkc+3t7SxevJgLL7ywZvVVoxYj/B72vt5mO7ClijI10TfKl/Do3qwJDcWxuOnTp/PCCy8A8Nvf/pZZs2Yxbdo0zjzzTDZt2gTAvHnzWLBgAaeffjrve9/7uOuuuwC44ooruPvudy6u9ulPf5qVK1eyfv165syZU7K+KVOmMHr0aNauXbvPc93d3Zx11llMmzaNj3zkI/T29gLw6KOPcvLJJzN9+nS+8pWvcNJJJ2Xudy0C/qNAh6SJkg4ALgNW9iuzEvisCj4EvBIRvTWou6QFMzo4dcLhHt2bNaF6H4vbvXs39913HxdddBEA8+fP5zvf+Q7d3d3cdNNNfPGLX9xTtre3lwceeIBVq1axcOFCAD73uc/xve99D4BXXnmFBx98kNmzZw9a7w033MA3vvGNvbbt3LmTa6+9lrvuuovu7m6uvvrqPb86rrrqKm655RYeeugh2traatL3zCmdiNgl6RrgHgrTMm+PiCckfSF5/hZgDYUpmc9QmJZ5VdZ6BzL20JHc8fnp9azCzOqk71f6j7qeZ+fuqNmxuDfeeIOpU6eyefNmpk2bxnnnnceOHTt48MEHueSSS/aUe+utt/bc/+hHP8p+++3HpEmTePHFFwE466yz+NKXvsTWrVv58Y9/zMUXX8z++w8eSs88s5CO2rBhw55tTz31FBs3buS8884DCl9G48aN4+WXX+a1117j9NNPB+BTn/oUq1atytR/qNE8/IhYQyGoF2+7peh+AF+qRV1m1vrqcSyuL4f/yiuvMGfOHJYsWcK8efMYM2bMntx+fwceeOCe+8UXi7riiitYtmwZy5cv5/bbb6+4DYsWLWLx4sV7viAigsmTJ/PQQw/tVe6ll15K0bPKeS0dMxt26nksbvTo0Xz729/mpptu4qCDDmLixInceeedQCEAP/7444PuY968eXzrW98CYPLkyRXXff755/PSSy/tqeP4449n27ZtewL+zp07eeKJJzjssMMYNWoUDz/8MADLly9P08WyHPDNbFiq57G4D37wg0yZMoXly5ezbNkybrvtNqZMmcLkyZP3OiBbzpFHHsmJJ57IVVelz04vWrSInp4eAA444ADuuusurr/+eqZMmcLUqVN58MEHAbjtttuYP38+06dPJyIYPXp06rr6G9bXtO3s7AxfAMWs+T355JOceOKJjW5Gzbz++ut84AMf4Be/+EVNAnEpO3bs4JBDDgHgxhtvpLe3l5tvvnmvMqXeV0ndEVFyfqhH+GZmKaxbt44TTjiBa6+9tm7BHmD16tV7ThTbsGEDN9xwQ+Z9erVMM7MUZs6cyXPPPVf3ei699FIuvfTSmu7TI3wzGxLDOX3cjKp5Px3wzazuRo4cyfbt2x30a6RvPfyRI9PNXnJKx8zqrr29nZ6entTrt1t5fVe8SsMB38zqbsSIEamuzGT14ZSOmVlOOOCbmeWEA76ZWU444JuZ5YQDvplZTjjgm5nlhAO+mVlOOOCbmeVEphOvJB0O/AiYAGwGPhkR+1yqRdJm4DVgN7Cr3NKdZmZWP1lH+AuB+yKiA7gveVzOOREx1cHezKwxsgb8ucD3k/vfBz6acX9mZlYnWQP+kRHRC5D8O7ZMuQDuldQtaf5AO5Q0X1KXpC4vtGRmVjuD5vAlrQOOKvHUohT1nBERWySNBdZK2hQR95cqGBFLgaVQuMRhijrMzGwAgwb8iJhZ7jlJL0oaFxG9ksYBW8vsY0vy71ZJK4DTgJIB38zM6iNrSmclcGVy/0pgn8u9SzpY0qi++8D5wMaM9ZqZWUpZA/6NwHmSfgOclzxG0tGS1iRljgQekPQ48AiwOiJ+lrFeMzNLKdM8/IjYDswosX0LMDu5/ywwJUs9ZmaWnc+0NTPLCQd8M7OccMA3M8sJB3wzs5xwwDczywkHfDOznHDANzPLCQd8M7OccMA3M8sJB3wzs5xwwDczywkHfDOznHDANzPLCQd8M7OccMA3M8sJB3wzs5xwwDczy4lMAV/SJZKekPS2pM4Bys2S9JSkZyQtzFKnmZlVJ+sIfyPwceD+cgUktQFLgAuAScDlkiZlrNfMzFLKek3bJwEkDVTsNOCZ5Nq2SFoOzAV+naVuMzNLZyhy+McAzxc97km2lSRpvqQuSV3btm2re+PMzPJi0BG+pHXAUSWeWhQRd1dQR6nhf5QrHBFLgaUAnZ2dZcuZmVk6gwb8iJiZsY4e4Niix+3Aloz7NDOzlIYipfMo0CFpoqQDgMuAlUNQr5mZFck6LfNjknqA6cBqSfck24+WtAYgInYB1wD3AE8Cd0TEE9mabWZmaWWdpbMCWFFi+xZgdtHjNcCaLHWZmVk2PtPWzCwnHPDNzHLCAd/MLCcc8M3McsIB38wsJxzwzcxywgHfzCwnHPDNzHLCAd/MLCcc8M3McsIB38wsJxzwzcxywgHfzCwnHPDNzHLCAd/MLCcc8M3McsIB38wsJ7Je4vASSU9IeltS5wDlNkv6laTHJHVlqdPMzKqT6RKHwEbg48A/VlD2nIj4Q8b6zMysSlmvafskgKTatMbMzOpmqHL4AdwrqVvS/IEKSpovqUtS17Zt24aoeWZmrW/QEb6kdcBRJZ5aFBF3V1jPGRGxRdJYYK2kTRFxf6mCEbEUWArQ2dkZFe7fzMwGMWjAj4iZWSuJiC3Jv1slrQBOA0oGfDMzq4+6p3QkHSxpVN994HwKB3vNzGwIZZ2W+TFJPcB0YLWke5LtR0takxQ7EnhA0uPAI8DqiPhZlnrNzCy9rLN0VgArSmzfAsxO7j8LTMlSj5mZZeczbc3McsIB38wsJxzwzcxywgHfzCwnHPDNzHLCAd/MLCcc8M3McsIB38wsJxzwzcxywgHfzCwnHPDNzHLCAd/MLCcc8M3McsIB38wsJxzwzcxyItN6+M1u9s0b+HXvq/tsnzTuUNZcd2YDWmRmVj+5CPjlAvth7xrBiDaxc/c710of0SZOee9hQ9k8M7MhkfUSh9+UtEnSLyWtkDSmTLlZkp6S9IykhVnqrMYp48cwok17bRvRJs45YSz7ae/tbRILZhw3lM0zMxsSWXP4a4GTIuJk4Gngq/0LSGoDlgAXAJOAyyVNylhvKgtmdJQM7AsvOIFLprXv+TIY0SY+0XksY0eNHMrmmZkNiUwBPyLujYhdycOHgfYSxU4DnomIZyPiT8ByYG6WetMae+jIsoG9+MvAo3sza2W1nKVzNfDTEtuPAZ4vetyTbBtS5QJ735eBhEf3ZtbSBg34ktZJ2ljiNreozCJgF7Cs1C5KbIsS2/r2NV9Sl6Subdu2VdKHigwU2BfM6ODUCYd7dG9mLW3QWToRMXOg5yVdCcwBZkREqUDeAxxb9Lgd2DJAfUuBpQCdnZ1lvxiqsWBGB09v3bFPYB976Eju+Pz0WlZlZjbsZJqWKWkWcD1wVkS8XqbYo0CHpInAC8BlwKey1FstB3Yzy7OsOfzvAqOAtZIek3QLgKSjJa0BSA7qXgPcAzwJ3BERT2Ss18zMUso0wo+IkknviNgCzC56vAZYk6UuMzPLxmvpmJnlhAO+mVlOOOCbmeWEA76ZWU444JuZ5YQDvplZTjjgm5nlhAO+mVlOOOCbmeWEA76ZWU444JuZ5YQDvplZTjjgm5nlhAO+mVlOOOCbmeWEA76ZWU444JuZ5UTWa9p+E7gQ+BPwW+CqiHi5RLnNwGvAbmBXRHRmqbeRZt+8gV/3vrrP9knjDmXNdWc2oEVmZpXJOsJfC5wUEScDTwNfHaDsORExtZmDPcAp48cwok17bRvRJk5572ENapGZWWUyBfyIuDe5SDnAw0B79iYNbwtmdLCf9g74bRILZpS8vK+Z2bBRyxz+1cBPyzwXwL2SuiXNr2GdQ27soSO5ZFr7nlH+iDbxic5jGTtqZINbZmY2sEFz+JLWAUeVeGpRRNydlFkE7AKWldnNGRGxRdJYYK2kTRFxf5n65gPzAcaPH19BF2pvsDz9ghkd3NndA4RH92bWNAYd4UfEzIg4qcStL9hfCcwBPh0RUWYfW5J/twIrgNMGqG9pRHRGROcRRxxRTZ8yGyxP3zfKl/Do3syaRqaUjqRZwPXARRHxepkyB0sa1XcfOB/YmKXeeqskT79gRgenTjjco3szaxpZc/jfBUZRSNM8JukWAElHS1qTlDkSeEDS48AjwOqI+FnGeuuqkjz92ENHcsfnp3t0b2ZNI9M8/IgoObxNUjizk/vPAlOy1NMIztObWavxmbZlOE9vZq0m0wi/1S2Y0cHTW3d4dG9mLcEBfwB9eXozs1bglI6ZWU444JuZ5YQDvplZTjiHX2deTtnMhguP8OvMyymb2XDhgF9nXk7ZzIYLB/w683LKZjZcOOAPgeJRvkf3ZtYoPmhbI4MdnL1kWjvLHnnOo3szaxiP8GtksIOzXk7ZzBrNAb9GBjs46+WUzazRHPBrxAdnzWy4c8CvIR+cNbPhzAG/hryGvpkNZ56lU2NZ19D3UgxmVi9ZL2L+d5J+mVzP9l5JR5cpN0vSU5KekbQwS53DXdaDs16KwczqJWtK55sRcXJETAVWAV/rX0BSG7AEuACYBFwuaVLGeluWl2Iws3rJFPAjojj3cDAQJYqdBjwTEc9GxJ+A5cDcLPW2Ms/2MbN6yZzDl7QY+CzwCnBOiSLHAM8XPe4B/sUA+5sPzAcYP3581uYNWwPl6v/pqlO5s7sHCI/uzaxmBh3hS1onaWOJ21yAiFgUEccCy4BrSu2ixLZSvwRI9rc0IjojovOII46otB9NZ6BcvWf7mFk9DDrCj4iZFe7rh8Bq4Ov9tvcAxxY9bge2VLjPlrVgRseeUXyf4tF81tk+Zmb9ZZ2l01H08CJgU4lijwIdkiZKOgC4DFiZpd5WMFiu3ksxmFmtZc3h3yjpeOBt4PfAFwCS6Zm3RsTsiNgl6RrgHqANuD0inshYb0soHuVnydV77r6ZVSJTwI+Ii8ts3wLMLnq8BliTpa5W1DfKz7ps8injx/Cbra+xc/c76SHP3Tez/ry0QoPVYtlkz903s0p4aYUG68vVZ93HJdPa+VHX8+zcHZ67b9aEhiI164DfRDx336x1DUVq1imdJuK5+2atayhSsx7hN5Fazd33rB6z4WcoUrMe4TeRWs3d94qcZsNTvS+i5IDfZGrxgfCsHrPhqd6pWad0mkwt5u57Vo9Z4wyWUq3nsioO+E2oFrl6z+oxa4zBZuPUYqp2OQ74TajSD8RAH6y0vxR8oNesNgabfFFPzuG3sMFy9WnO8vWBXrPaaORFjhzwW1gtV+T0gV6z2qn3bJxyHPBbXK0+WL70olntNOpESefwW1wtc/VpDvQ65282sEZc5Mgj/ByoVa4+zajEOX+zgTXiIkcO+DlQy1x9pV8ezvmbDT9O6dheBjspq9IpoT65y/KiXPrywP33461db++zvZFpzUwBX9LfAXMpXOJwKzAvudpV/3KbgdeA3cCuiOjMUq/VV5pLL9by5C7n/a0ZlTvfZfzh72Lz9j8OqyvRZU3pfDMiTo6IqcAq4GsDlD0nIqY62A9/tcrVp52J4Ly/NaNy6cubL5s67NKamQJ+RBQPxw6m+NQxa2q1ytWnOWDsvL81o3JTlicdPXrYTWXOnMOXtBj4LPAKcE6ZYgHcKymAf4yIpQPsbz4wH2D8+PFZm2dVqlWuPs26IGnz/k4B2XBRLg2aJj06FAYN+JLWAUeVeGpRRNwdEYuARZK+ClwDfL1E2TMiYoukscBaSZsi4v5S9SVfBksBOjs7/YuhCdQq59+3UmCl+xqKS8JZvlQ7iCh3vkstVretpUEDfkTMrHBfPwRWUyLg9x3IjYitklYApwElA741nzQf6kpWCqx0X2kXofIvAhtMuc9n7ytvMGHh6n3KF392yp1I1YgTrMrJOkunIyJ+kzy8CNhUoszBwH4R8Vpy/3zgb7PUa8NPpR/qSoJ0/30NFKjTpID8i8AGU+7zee4JY1n5+JYBPzvl0pf1XO44raw5/BslHU9hWubvgS8ASDoauDUiZgNHAitUOBi3P/DDiPhZxnptmKnl/Pz++xooUC8497iKU0DVLEvrXwXNLe3fr9znc8G5x7Hql700YknjWso6S+fiiDgpmZp5YUS8kGzfkgR7IuLZiJiS3CZHxOJaNNyaV9oF3QaavZNm6mc1C8B5qmhzq+bvV+rz2SqLB/pMWxtyaQ9kDfarIE0KKO3JYD5OMHxU896W+/s98rvtA+bkS30+h9uMm2o44FtDpD2QNdB/tjQpoFp/2fRXzXGCvH5JpO13Ne9tub8fEfzuD+XPgi31+RxuM26q4YBvDZH2QFYtZ++U+s9cq18F1RwnSBvIqvmCSPuaepeH9P2udrRecrAQDPh3Kvf5HE4zbqrh1TKtaVR61m41V/pKu0TE7Js3MGHh6n1u8773aOpcb9ozjKvJS6d9Tb3LQ/p+l/u7njbh8AHrLvX3qzYn34gljWvJI3xrGml+FaTNt6b9VVDNzKFaTS+t5ldE2tfUuzyUT7fMu/3RdL+2Bhmt97Wv/8i8FXLyaTngW0uqda6+/5fNQAFu7KjSdaf9kqjVF0Ql/atn+bTpsm+v+03qYzBpp/oWt7mZc/JpOeBby6rlgeH+0s4c6r//PgN9SdTyC6JcLjttMK7mWEfag+jVHIOpdrTe7Dn5tJzDt5aVNt+adjnngc4nKFX3YHnj/sco0p5/MFgePe1ralF+sH6U6nc1x2CqvSh4s+fk03LANyuSZjnnaoJMmi+JWn5BlOtf2mBczXLY1QTwtCfnlavb9uaAb1Yk7YgvbZCp56+ISmaepH1N1vKV9KMW79NAdds7HPDNMqgmyNTzV0S1I+Nql7qotI5qfw15xF5bDvhmQ6yevyKqHRmneU21+fJqfg15xF5bnqVjNsylPSu5mpkn1cxoSlvHcFomOK8UMXwvKtXZ2RldXV2NboaZWdOQ1B0RnaWec0rHzCwnHPDNzHLCAd/MLCcc8M3McsIB38wsJ4b1LB1J2yhcHL0a7wH+UMPmNAv3O1/c73yppN/vjYgjSj0xrAN+FpK6yk1NamXud7643/mStd9O6ZiZ5YQDvplZTrRywF/a6AY0iPudL+53vmTqd8vm8M3MbG+tPMI3M7MiDvhmZjnRcgFf0ixJT0l6RtLCRrenniTdLmmrpI1F2w6XtFbSb5J/D2tkG2tN0rGS/qekJyU9Iem6ZHur93ukpEckPZ70+98k21u6330ktUn6P5JWJY/z0u/Nkn4l6TFJXcm2qvveUgFfUhuwBLgAmARcLmlSY1tVV/8EzOq3bSFwX0R0APclj1vJLuAvI+JE4EPAl5K/cav3+y3g3IiYAkwFZkn6EK3f7z7XAU8WPc5LvwHOiYipRfPvq+57SwV84DTgmYh4NiL+BCwH5ja4TXUTEfcD/6/f5rnA95P73wc+OpRtqreI6I2IXyT3X6MQBI6h9fsdEbEjeTgiuQUt3m8ASe3AnwO3Fm1u+X4PoOq+t1rAPwZ4vuhxT7ItT46MiF4oBEdgbIPbUzeSJgAfBP43Oeh3ktZ4DNgKrI2IXPQb+BbwV8DbRdvy0G8ofKnfK6lb0vxkW9V9b7VLHKrENs87bUGSDgH+O/AXEfGqVOpP31oiYjcwVdIYYIWkkxrcpLqTNAfYGhHdks5ucHMa4YyI2CJpLLBW0qYsO2u1EX4PcGzR43ZgS4Pa0igvShoHkPy7tcHtqTlJIygE+2UR8eNkc8v3u09EvAysp3D8ptX7fQZwkaTNFFK050r6Aa3fbwAiYkvy71ZgBYW0ddV9b7WA/yjQIWmipAOAy4CVDW7TUFsJXJncvxK4u4FtqTkVhvK3AU9GxH8seqrV+31EMrJH0kHATGATLd7viPhqRLRHxAQK/59/HhGfocX7DSDpYEmj+u4D5wMbydD3ljvTVtJsCjm/NuD2iFjc2BbVj6T/BpxNYcnUF4GvA/8DuAMYDzwHXBIR/Q/sNi1JHwY2AL/inZzuX1PI47dyv0+mcICujcJA7Y6I+FtJ76aF+10sSel8OSLm5KHfkt5HYVQPhfT7DyNicZa+t1zANzOz0lotpWNmZmU44JuZ5YQDvplZTjjgm5nlhAO+mVlOOOCbmeWEA76ZWU444JullJwB+d1keWKzpuGAb5beF4ADgQ83uiFmaTjgm6U3C3gaeKzB7TBLxQHfLAVJIymsZ3MK8L8a3ByzVBzwzdLpoBDwN0XEzkY3xiyNVrsAilm9HQG8nxa+dKa1Lo/wzdI5msLFV/aTdFijG2OWhgO+WYUk7U8hd38UcAuwu7EtMkvH6+GbmeWER/hmZjnhgG9mlhMO+GZmOeGAb2aWEw74ZmY54YBvZpYTDvhmZjnx/wF7we5iEYS7hwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEYCAYAAABBS301AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAds0lEQVR4nO3dfZRU1Znv8e9Dg7YGxPaFICDSSmsEhBZaI45ioqCooBgHiRoUiNNm+QKdwUQjriSTyCwzehMk9FwXCYgaIwLKqEBeQIcLLvFqNwGFqLka3xoRkPAigy8Iz/2jqtsGqqrrVJ3q6lP1+6zVi6pTdc7Zu6l+ap9n77O3uTsiIhJd7fJdABERyY4CuYhIxCmQi4hEnAK5iEjEKZCLiEScArmISMQpkIsEZGZ3mtlvk7x2rpm9kcWx3cx6Z146KUYK5NImmVmFmX1qZr9L8Z5xZrbXzHaZ2U4zW2tmI3JdNnf/d3e/IclrK939lFyXQaQ5BXLJKTNrn+GutcDLabxvlbt3BI4E/hOYa2ZHZnhOkUhSIJfQmdk7Zna7mb0C/E/QYG5m3wa2A8+mu4+77wMeAb4CVMSPc6iZ3Wdm75nZJjN7wMwOi7/2DTNrMLPJZrbZzDaa2fj4a2fE399UbjO70szWxB//NNmVQuNxD/hd3GZmr5jZDjN73MxKm73+g/i5PzCzCQccK1X5l5jZ/2r23sfNbHa6vy8pLArkkitXA5cCR7r7F2a2yMy2J/lZ1LiTmR0B/AyYHORkZlYCjAf2AO/GN/8COBmoBHoD3YEfN9utK9A5vv27QK2Zlbn7y8BWYFiz936H2BdFJq4ChgPlQH9gXLzMw4Hb4uepAIYesF+q8k8AxprZ+WZ2LXAGMCnD8knEZXrZK9KS6e7+fuMTd083d/1zYJa7v29m6bz/LDPbTqwl/gXwHXffbLGd/wXo7+7/ADCzfwd+D/wovu8e4Gfu/gWwxMx2AacALwIPEQvefzCzo4CLgJvSrMOBprv7B/EyPEMsMEMswD/o7uvir/2U2BcgLZXf3T80s+/Fy3kYMMrdP86wfBJxCuSSK++3/Jb9mVklsVbp6QF2e9HdzzGzjsAs4FxgHnAscDhQ3+wLwYCSZvtujQfxRruBjvHHvwNeix/3KmClu28MVqMmHx5wjm7xx92A+mavvdvscTrlXwTMAN5w9+czLJsUAAVyyZX9ptU0sz8QC7KJrHT3i4FvAL2A9+LBqyNQYmZ93H1gypO57zKzm4C34rnitcAnQF933xC48O4bzGwVcAUwFvjfQY+Rho3A8c2e92z2+CNaLv9U4DWg3MyudvfHclBGiQDlyKVVuPvF7t4xyc/F8bfNBE4ilnqoBB4AFhNLa6Rzjq3Ab4Efxzs/fwP8ysy6AJhZdzNL61hxDwM/BE4DFgbYL13zgHFm1sfMDgd+0vhCS+U3syHE+gSui//82sy656CMEgEK5NJmuPtud/+w8QfYBXzq7lsCHGYacImZ9QduB94EXjSzncAyYjnwdC0ETgAWuvv/BNgvLe7+B2LlfS5ezucOeEvC8sc7hB8GbnH3DfG0yizgQUuzY0EKi2lhCZHkzOwt4EZ3X5bvsogkoxa5SBJmdiWxXP+BLWWRNkWdnSIJmNlyoA8wNp6vFmmzlFoREYk4pVZERCIuL6mVY445xnv16pWPU4uIRFZ9ff1H7n7sgdvzEsh79epFXV1dPk4tIhJZZvZuou1KrYiIRJwCuYhIxCmQi4hEnMaRi0hCe/bsoaGhgU8//TTfRSk6paWl9OjRgw4dOqT1fgVyEUmooaGBTp060atXLzSFS+txd7Zu3UpDQwPl5eVp7RNKII+vdHI/sbmSf+vu94Rx3EazfzGJpTu6s2pf36Ztg9utZ1jnDUy4/f4wTyUicZ9++qmCeB6YGUcffTRbtqQ/V1zWOfL4Elu1wMXEbmm+2sz6ZHvc5rzbQGo7TGdwu/VALIjXdpiOd0s5RbWIZElBPD+C/t7D6Ow8E3jT3f/u7p8Dc4HLQzhuk5GXj6FmXw0zOkzn++3nx/7dV8PIUWPCPI2ISCSFEci7s/+yXg3xbfsxs2ozqzOzuiCXDABdjiil58CLeGzfMCa1X8hj+4Zx/KDhdOlU2vLOIhJJJSUlVFZW0q9fP0aOHMn27dtDPX5dXR0TJ048aPvy5csxM5555pmmbSNGjGD58uWhnj9MYeTIE10DHDQTl7vPJLYCDFVVVYFn6ppcsQlfs5T7v7iCsSVLsYprgX6BCysi4bvk/pX8dePOg7b3Oe4IlkxKtsJfaocddhhr1qwB4Prrr6e2tpYpU6ZkU8z9VFVVUVVVlfC1Hj16MHXqVEaOHBna+XIpjBZ5A/uvO9gD+CCE437p7RWULalmYe+pTNs7moW9p1K2uBreXhHqaUQkMwN7HkmHkv3bdB1KjIEnlIVy/MGDB7NhQ2zp0rfeeovhw4czaNAgzj33XF5//XUAxo0bx8SJEzn77LM58cQTWbBgAQBjx47lqaeeajrWtddey9NPP83y5csZMWJEwvMNGDCAzp07s3Tp0oNeq6+v57zzzmPQoEFcdNFFbNwYW5P75Zdfpn///gwePJgf/OAH9OvXeg3NMAL5y0CFmZWb2SHAt4GnQzjulzashtFzGHn5GM7odVQsNz56Tmy7iOTdxAsqaHdAB12JGRMv6J31sffu3cuzzz7LZZddBkB1dTW//vWvqa+v57777uOmm25qeu/GjRt5/vnnWbRoEXfccQcAN9xwAw8++CAAO3bs4IUXXuCSSy5p8bx33XUXd999937b9uzZw6233sqCBQuor69nwoQJTVcJ48eP54EHHmDVqlWUlJRkXe8gsk6tuPsXZnYL8Cdiww9nu/v6rEvW3Dk1AHQB5t04OLat0xAoHxLqaUQkM12OKGX0oB48Xvc+e/Y6HUqMf646Pqt+rE8++YTKykreeecdBg0axLBhw9i1axcvvPACo0ePbnrfZ5991vR41KhRtGvXjj59+rBp0yYAzjvvPG6++WY2b97Mk08+yZVXXkn79i2HvnPPjaWEVq5c2bTtjTfeYN26dQwbNgyIfckcd9xxbN++nY8//pizzz4bgGuuuYZFixZlXPegQhlH7u5LgCVhHEtEomniBRXMr28APJTWeGOOfMeOHYwYMYLa2lrGjRvHkUce2ZQ7P9Chhx7a9Lj5ojljx47l0UcfZe7cucyePTvtMkyZMoWpU6c2BX53p2/fvqxatWq/923bti1AzcKnuVZEJBSNrXIzsm6NN9e5c2emT5/Offfdx2GHHUZ5eTnz588HYoF17dq1LR5j3LhxTJs2DYC+ffumfnMzF154Idu2bWs6xymnnMKWLVuaAvmePXtYv349ZWVldOrUiRdffBGAuXPnBqli1hTIRSQ0Ey+o4IxeR4WSG2/u9NNPZ8CAAcydO5dHH32UWbNmMWDAAPr27btfR2YyX/3qVzn11FMZP3584HNPmTKFhoYGAA455BAWLFjA7bffzoABA6isrOSFF14AYNasWVRXVzN48GDcnc6dOwc+V6bysmZnVVWVa2EJkbbttdde49RTT813MUKxe/duTjvtNFavXp2zALtr1y46duwIwD333MPGjRu5//7MpxBJ9Ps3s3p3P2jMpFrkIlLQli1bxte+9jVuvfXWnLaSFy9e3HQD08qVK7nrrrtydq4DafZDESloQ4cO5b333sv5ecaMGcOYMfmZNkQtchGRiFMgFxGJOAVyEZGIUyAXEYk4BXIRaZMKcRrbNWvWMHjwYPr27Uv//v15/PHHsz4mKJCLSBien3bwbKRvr4htz1DjLfrr1q3jqKOOora2NqsiHqiqqorp06cnfK1xGtuwHX744Tz88MOsX7+eP/7xj9TU1ITyBaVALiLZ6z4Q5o/7Mpi/vSL2vHs4yzEWyjS2J598MhUVFQB069aNLl26BFqbMxkFchHJXvmQ2NTS88fBc1Nj/46eE8oMpYU6je1LL73E559/zkknnZTeLyIF3RAkIuEoHwJV34UV/wFDfph1EC/kaWw3btzI2LFjeeihh2jXLvv2tAK5iITj7RVQNysWxOtmQfm5WQXzQp3GdufOnVx66aXcfffdnHXWWWnvl4pSKyKSvcac+Og5cP6UL9MsISzHWEjT2H7++edcccUVXHfddftdVWRLgVxEshdfjrGpBd6YMw9pOcZCmcZ23rx5rFixgjlz5lBZWUllZWXSq4sgNI2tiCSkaWyD0TS2IiI5omlsRUQiTtPYikhRy0fqVYL/3hXIRSSh0tJStm7dqmDeytydrVu3Ulqa/uLVSq2ISEI9evSgoaEhlFvIJZjS0lJ69OiR9vsVyEUkoQ4dOlBeXp7vYkgalFoREYk4BXIRkYhTIBcRiTgFchGRiFMgFxGJOAVyEZGIUyAXEYm4rAK5mY02s/Vmts/MDpqRS0REci/bFvk64FtA9rPHi4hIRrK6s9PdXwMws3BKIyIigbVajtzMqs2szszqNHeDiEh4WmyRm9kyoGuCl6a4e8trLMW5+0xgJsRWCEq7hCIiklKLgdzdh7ZGQUREJDMafigiEnHZDj+8wswagMHAYjP7UzjFEhGRdGU7amUhsDCksoiISAaUWhERiTgFchGRiFMgFxGJOAVyEZGIUyAXEYk4BXIRkYhTIBcRiTgFchGRiFMgFxGJOAVyEZGIUyAXEYk4BXIRkYhTIBcRiTgFchGRiFMgFxGJOAVyEZGIUyAXEYk4BXIRkYhTIBcRiTgFchGRiFMgFxGJOAVyEZGIUyAXEYk4BXIRkYhTIBcRiTgFchGRiFMgFxGJOAVyEZGIUyAXEYk4BXIRkYjLKpCb2b1m9rqZvWJmC83syJDKJSIiacq2Rb4U6Ofu/YG/AT/KvkgiIhJEVoHc3f/s7l/En74I9Mi+SCIiEkSYOfIJwB+SvWhm1WZWZ2Z1W7ZsCfG0IiLFrX1LbzCzZUDXBC9Ncfen4u+ZAnwBPJrsOO4+E5gJUFVV5RmVVkREDtJiIHf3oaleN7PrgRHABe6uAC0i0spaDOSpmNlw4HbgPHffHU6RREQkiGxz5DOATsBSM1tjZg+EUCYREQkgqxa5u/cOqyAiIpIZ3dkpIhJxCuQiIhGnQC4iEnEK5CIiEadALiIScQrkIiIRp0AuIhJxCuQiIhGnQC4iEnEK5CIiEadALiIScQrkIiIRp0AuIhJxCuQiIhGnQC4iEnFZzUfels3+xSSW7ujOqn19m7YNbreeYZ03MOH2+/NYMhGRcBVsi9y7DaS2w3QGt1sPxIJ4bYfpeLeBeS6ZiEi4CrZFPvLyMdTcu4UZHabxu71D+U7JMr6/r4Z7R43Jd9FEREJVsC3yLkeU0nPgRTy2bxiT2i/ksX3DOH7QcLp0Ks130UREQlWwLXKAyRWb8DVLuf+LKxhbshSruBbol+9iiYiEqmBb5Ly9grIl1SzsPZVpe0ezsPdUyhZXw9sr8l0yEZFQFW6LfMNqGD2HkUefyZ92/4WRoy6Aj8pj28uH5Lt0IiKhMXdv9ZNWVVV5XV1dq59XRCTKzKze3asO3B7pFrnGiouIRDxHrrHiIiIRb5FrrLiISMRb5BorLiIS8RY5aKy4iEikW+QaKy4iEvUWucaKi4hoHLmISFQkG0eeVWrFzH5uZq+Y2Roz+7OZdcvmeCIiEly2OfJ73b2/u1cCi4AfZ18kEREJIqtA7u47mz39CtD6eRoRkSKXdWenmU0FrgN2AN9M8b5qoBqgZ8+e2Z5WRETiWuzsNLNlQNcEL01x96eave9HQKm7/6Slk6qzU0QkuIwnzXL3oWme4/fAYqDFQC4iIuHJdtRKRbOnlwGvZ1ccEREJKtsc+T1mdgqwD3gX+F72RRIRkSCyCuTufmVYBRERkcxEe64VERFRIBcRiToFchGRiIv27Ich0vqfIhJVapHHaf1PEYkqtcjjtP6niESVWuRxWv9TRKJKLfJmtP6niESRWuSNtP6niESUWuSNtP6niESU1uwUEYmIjKexLTQaLy4ihabocuQaLy4ihaboWuQaLy4ihaboWuQaLy4ihaboWuSg8eIiUliKrkWu8eIiUmiKr0Wu8eIiUmA0jlxEJCI0jjxkGo8uIm1F8eXIQ6Lx6CLSVqhFniGNRxeRtkIt8gxpPLqItBVqkWdB49FFpC1QizxTGo8uIm2EWuQtSDY65V8P/yNnXKvx6CKSf2qRtyDZ6JRXTrgeyofQ5YhS5t04OJYbLx8C59Tkt8AiUnTUIm+BRqeISFunFnkLNDpFRNo6tcjToNEpItKWhRLIzew24F7gWHf/KIxjthnx0Smzek9l2l+PoeMp3+S7i6vh8DmBOzV1W79Iccr1337WqRUzOx4YBryXdWnaosbZEi8fwxm9jmLkqDEwek5se0C6rV+kOOX6bz/r2Q/NbAHwc+ApoCqdFnmxzn64eeen3HbvDH7V7suO03/dV8O9P7xFOXeRAhbW336y2Q+zapGb2WXABndfm81xioU6TkWKU67/9lvMkZvZMqBrgpemAHcCF6ZzIjOrBqoBevbsGaCIhUUdpyLFKZd/+y0Gcncfmmi7mZ0GlANrzQygB7DazM509w8THGcmMBNiqZVsCt3WJevYGNuxjkva14XScSoiERLioIlEMk6tuPur7t7F3Xu5ey+gARiYKIgXm2QdG8d0PDS0jlMRiZAQB00kEtpSb2b2DursBNSpKSK5kfOl3uKtcqFZx8bq9Uxq/yQz9n4ro44NjTsXkXToFv0cmVyxiWvaxTo2rmm3lMkVwTNOGncuIunQLfq5EFLHhibsEomOfF5Bq0WeCyF1bGjcuUh05PMKOrTOziAKvbMzTNvWLcPnj+ORvUMZW7IMu2oOZX0TjggVkTxqjUEOOe/slMxp3LlI9IU1yCETCuRtgHcbSO3uKdy8ZyKr9vVtuiR7s+M3YWSw5eQ00kUkf/J157Zy5G3AyMvHULOvhhkdpvP99vNj/+6rodf43wReTk4jXUTyJI8LsqtF3gaEeUmmkS4iedI4yCEPC7IrkLcRYV2S5TNPJ1LU4lfKXYB5Nw6Obes0pFX6sxTI24KQJ9TRDIsiudMW+6EUyNuCgJdkKT9IV10Z6EuhLX4oRdqyZIMTFnabmrcyqbOzLTinJlCnZsoOzYA3I6lzVCSYZIMTRuaxH0ot8ghK2aEZz4Wnm6dT56hIYqmuVnsOvKFN9UOpRR5BYd66r2kARBJLdbUaxqR4YVKLPKKCdmimal1MvnRkaMdSXl0KRbKr1WkDu1O2OHer/WRCLfIoyuDGg2Sti65lHUM7lvLqUkiSXa2WbVvX5lb6Uos8ijK48SBZ66L2JODEcI6lvLoUmoRXvn1rgPyMF09GgTyKMrjxINmNQh0v6BfasZRXl7YqVToQiPykdQrkRSRIXr2lPHiYxxLJtZbGfoc1aV2+KEdeLALm1VPmwcM8lkgrSDX2O8xJ6/JFLfJiETCvnjIPvvaB8I6VgFrwEraW0oFRTxUqkBeLgHn1lB/8JMeaPe8Jli75j6xvoGiLt0BL25LJl32qdGDU5ydSIJekgn64UwXgIMfSqBhpSbLP2nOffJ2r7zy4MZGy4xJCnbQuHxTIJbEMZmQM6waKTEbFKB1TXJJ91r5+4lHUvvPjYB2XkLd5xMOiQC6JZTBWPVkALtu2LOGxXn70p/xy9/BQ7jYNmo5R4G9bgv5/JPus/fP5vam59x8HX82NvwU6laZMLbalceFBKZBLYhlOkh/kBopXTrie2jcPDr4vll0f2tVAsnSM8vBtS9BUSbIv+7Ij+kW+4zITCuQSnoDpmEzuNp0974lQZqTLJA+vVnx6Mvk9BU2VpPqyn1zxeaQ7LjOhQC7hCZiOyeRuU3/7oVA6VDPJw7dG+ibXXxatUaZMWtcTbr8/UKok6Zf9q09Q9vozke64zIQCuYQng3RM0JExQTtUlzx1Bo/sGhRKHr410je5/rJojTIFbV03HifRZyFZqiTpl30eF0DOJwVyyZ8MRsYE7VA95pnfUPvZ9LQvzVMF/mStxjDTN7n+sgizTOdsqOfqOzelneZK2hE5akzKz0KgVEkeF0DOJwVyyZ8MW09BOlR7jTuTmntnpH1pnizwp2o1zv7Fb0JbSSZoyidoYE51/KAr4rT7uJ7aT9JPc6XsiFyb5LNQpKmSoLIK5Gb2U+BfgC3xTXe6+5JsCyVFIpPWU8BWfNA8fLLAn6rV2LXsamp3Jw7+QQN/sgnJgn5ZLPrP2wMvJBL0hq49xwe/byBp6zrZZ6FIUyVBhTFp1q/cvTL+oyAuuRVwcWkg0LJcKZe+S3LuISeVJZx06VsDuyecXKxrWcfAE5Kl2idR/TJZSCTZ5FHJ6tFl60vBFl549YnAi5gEXZi8WCm1ItEStBWfQR4+aQdsknN3LB9Cz52vpp23H/L3l7jp/ZpAE5Il2ydZ6/dbl8xk0t8Ofn+qoZ1dyocE6n9gw2omV/RLf+EFta5zxtw9851jqZVxwE6gDpjs7tuSvLcaqAbo2bPnoHfffTfj84qk7flp0H0gm48+k1se+wszrjmdLh+9FAseiVp1b6+A+eOYddxPuPuvx3BXn4/47sZ/i7UoUwSbbeuW4fPH8cjeoYwtWYZdNYeyvkOTvv+uha/SdfUvuaUkFjA/HDSZu0elHuuccJ9jliWt311bLgh8jkD1yPB3JZkzs3p3rzpoe0uB3MyWAV0TvDQFeBH4CHDg58Bx7j6hpcJUVVV5XV1dOuUWaV1BAz9kFNCCBv5M9gl8jqD1yOR3JVnJOJAHOEEvYJG7t3gLlQK5FJTWaPUH3SeTcygwt3k5CeRmdpy7b4w//j7wdXf/dkv7KZBLUcskYAbdR0G5IOUqkD8CVBJLrbwD3NgY2FNRIBcRCS5ZIM9q1Iq7j81mfxERyZ4WXxYRiTgFchGRiFMgFxGJOAVyEZGIC20ceaCTmm0BMr218xhiNyEVG9W7+BRr3VXv5E5w92MP3JiXQJ4NM6tLNPym0KnexadY6656B6fUiohIxCmQi4hEXBQD+cx8FyBPVO/iU6x1V70DilyOXERE9hfFFrmIiDSjQC4iEnGRCuRmNtzM3jCzN83sjnyXJ1fMbLaZbTazdc22HWVmS83s/8X/LctnGXPBzI43s/82s9fMbL2ZTYpvL+i6m1mpmb1kZmvj9f63+PaCrncjMysxs7+Y2aL484Kvt5m9Y2avmtkaM6uLb8u43pEJ5GZWAtQCFwN9gKvNrE9+S5Uzc4DhB2y7A3jW3SuAZ+PPC80XxJYLPBU4C7g5/n9c6HX/DDjf3QcQmxZ6uJmdReHXu9Ek4LVmz4ul3t+ML1rfOHY843pHJpADZwJvuvvf3f1zYC5weZ7LlBPuvgL4xwGbLwceij9+CBjVmmVqDe6+0d1Xxx9/TOyPuzsFXneP2RV/2iH+4xR4vQHMrAdwKfDbZpsLvt5JZFzvKAXy7sD7zZ43xLcVi682LtoR/7dLnsuTU/GlA08H/i9FUPd4emENsBlY6u5FUW9gGvBDYF+zbcVQbwf+bGb18YXpIYt6Z7WwRCuzBNs0drIAmVlH4Amgxt13miX6ry8s7r4XqDSzI4GFZtbi2rdRZ2YjgM3uXm9m38hzcVrbP7n7B2bWBVhqZq9nc7AotcgbgOObPe8BfJCnsuTDJjM7DmJrpRJruRUcM+tALIg/6u5PxjcXRd0B3H07sJxYH0mh1/ufgMvM7B1iqdLzzex3FH69cfcP4v9uBhYSSx1nXO8oBfKXgQozKzezQ4BvA0/nuUyt6Wng+vjj64Gn8liWnLBY03sW8Jq7/7LZSwVddzM7Nt4Sx8wOA4YCr1Pg9Xb3H7l7D3fvRezv+Tl3/w4FXm8z+4qZdWp8DFwIrCOLekfqzk4zu4RYTq0EmO3uU/Nbotwws8eAbxCb1nIT8BPgv4B5QE/gPWC0ux/YIRppZnYOsBJ4lS9zpncSy5MXbN3NrD+xzq0SYo2ree7+MzM7mgKud3Px1Mpt7j6i0OttZicSa4VDLL39e3efmk29IxXIRUTkYFFKrYiISAIK5CIiEadALiIScQrkIiIRp0AuIhJxCuQiIhGnQC4iEnEK5CLNxO+6mxGfRlYkEhTIRfb3PeBQ4Jx8F0QkXQrkIvsbDvwNWJPncoikTYFcJM7MSonNdzIQ+D95Lo5I2hTIRb5UQSyQv+7ue/JdGJF0RWlhCZFcOxY4mQJdQlAKl1rkIl/qRmxRi3aFuHK7FC4FchHAzNoTy413BR4A9ua3RCLp03zkIiIRpxa5iEjEKZCLiEScArmISMQpkIuIRJwCuYhIxCmQi4hEnAK5iEjE/X8ITZSh+PC48gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# print(8*np.log(2))\n",
    "ent_plotter(r=2)\n",
    "ent_plotter(r=3)\n",
    "ent_plotter(r=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-15-8f879df0b7f2>:13: RuntimeWarning: divide by zero encountered in log\n",
      "  plt.plot(np.log(Lsub_list),LN_ex,'s')\n",
      "<ipython-input-15-8f879df0b7f2>:14: RuntimeWarning: divide by zero encountered in log\n",
      "  xplot=np.log(np.tan(pi*Lsub_list/N))\n",
      "<ipython-input-15-8f879df0b7f2>:15: RuntimeWarning: invalid value encountered in subtract\n",
      "  yplot=(xplot-xplot[0])/4+ 0.47\n",
      "<ipython-input-15-8f879df0b7f2>:16: RuntimeWarning: divide by zero encountered in log\n",
      "  plt.plot(np.log(Lsub_list),yplot,label='CFT')\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEMCAYAAAAxoErWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZWklEQVR4nO3df4xlZZ3n8feHpntRmtVduv0BTQtLmDjgSI/WNpo2I5jRNIwOcceYZvw1RtPByEYzxh00E1HXzZpMYnZVlGmdHtbsADEqLmtAILuuzGhwunFbaQRMixrKVkFUfow/2/nuH/cUe62+p/rWj1PnVtX7ldz0Pc/znHu/fSjq29/nPOecVBWSJI1yXN8BSJIml0lCktTKJCFJamWSkCS1MklIklqZJCRJrTpLEklOS/KFJHcnuSvJW0aMSZIPJjmU5OtJnjPUtzPJvU3f5V3FKUlq12UlcQR4W1X9LvA84M1Jzp415kLgrOa1G/goQJJ1wJVN/9nAJSP2lSR1rLMkUVXfr6qvNu8fBe4GTp017GLgEzVwO/DkJE8HtgOHquq+qvoVcF0zVpK0jI5fji9Jcjrw+8BXZnWdCtw/tD3dtI1qP6/ls3czqEI48cQTn/vMZz5zaYKWpDXgjjvu+FFVbW7r7zxJJNkIfBp4a1U9Mrt7xC41R/vRjVV7gD0AU1NTtX///kVEK0lrS5LvztXfaZJIsp5Bgvi7qvrMiCHTwGlD21uAw8CGlnZJ0jLqcnVTgL8B7q6qD7QMuwF4bbPK6XnAw1X1fWAfcFaSM5JsAHY1YyVJy6jLSmIH8BrgziQHmrZ3AlsBquoq4EbgIuAQ8DPg9U3fkSSXATcD64C9VXVXh7FKkkboLElU1T8w+tzC8JgC3tzSdyODJLIov/71r5menuYXv/jFYj9qYpxwwgls2bKF9evX9x2KpFVuWVY39Wl6epqTTjqJ008/ncEM2MpWVTz00ENMT09zxhln9B2OpFVu1d+W4xe/+AUnn3zyqkgQAEk4+eSTV1VlJGlyrfokAayaBDFjtf19JE2uNZEkJEkLY5JYBj/4wQ/YtWsXZ555JmeffTYXXXQR3/zmN3nCE57Atm3bHn+95z3vefz9unXrHn//wQ9+sO+/gqQ1atWfuO5bVfHyl7+c173udVx33XUAHDhwgB/+8IeceeaZHDhw4LfGX3HFFQBs3LjxqD5JAph636386LFfHdW+aeMG9v/li5f0u6wkOvaFL3yB9evXc+mllz7etm3bNk477bQ59pKkdqMSxFzti7G2KombLocf3Lm0n/m034ML39/affDgQZ773OeO7PvWt77Ftm3bANixYwdXXnnl0sYmSYu0tpLEhBk13SRJk2RtJYk5/sXflXPOOYdPfepTy/69krQUPCfRsRe96EX88pe/5GMf+9jjbfv27eO7353z7rySNBFMEh1LwvXXX8+tt97KmWeeyTnnnMO73/1uTjnllL5Dk7RCbdq4YV7ti7G2ppt6csopp/DJT37yqPaDBw+27vPYY491GZKkFWypl7nOxUpCktTKJCFJarUmksTgsRWrx2r7+0iaXKv+nMQJJ5zAQw89tGpuFz7zPIkTTjih71AkLbPlvB3HjFWfJLZs2cL09DQPPvhg36EsmZkn00laW5bzdhwzOksSSfYCLwUeqKpnjeh/O/CqoTh+F9hcVT9O8h3gUeA3wJGqmlpoHOvXr/cJbpK0QF2ek7ga2NnWWVV/VVXbqmob8A7gi1X146EhFzT9C04QkqTF6SxJVNVtwI+POXDgEuDarmKRJC1M76ubkjyRQcXx6aHmAm5JckeS3f1EJkmahBPXLwO+NGuqaUdVHU7yFODWJPc0lclRmiSyG2Dr1q3dRytJPdm0cUPr6qauTEKS2MWsqaaqOtz8+UCS64HtwMgkUVV7gD0AU1NTXkAgadVazttxzOh1uinJk4AXAv9jqO3EJCfNvAdeArTf5EiS1Jkul8BeC5wPbEoyDVwBrAeoqquaYS8Hbqmqfxra9anA9c2Fb8cD11TV57uKU5ImXR8X0c3oLElU1SVjjLmawVLZ4bb7gHO7iUqSVp4+LqKb0fvqJknS5DJJSJJamSQkSa1MEpKkVpNwnYQkaYS2VU0zuryIboZJQpIm1FwJ4jvv/6NlicHpJklSK5OEJKmVSUKS1MokIUkTaOp9t/YdAmCSkKSJ1PeqphkmCUlaYZbzluEmCUlSK5OEJKmVF9NJ0gQ51lXWy81KQpImyLESxHKetAYrCUlaMZbrVhzDrCQkSa2sJCRpAkzauYgZnVUSSfYmeSDJwZb+85M8nORA83rXUN/OJPcmOZTk8q5ilKRJMYkJArqdbroa2HmMMX9fVdua13sBkqwDrgQuBM4GLklydodxStLEW+4T1jM6m26qqtuSnL6AXbcDh6rqPoAk1wEXA99YwvAkaSKMM83UxwnrGX2fuH5+kq8luSnJOU3bqcD9Q2Omm7aRkuxOsj/J/gcffLDLWCVpyU3qNNOMPpPEV4FnVNW5wIeAzzbtGTG22j6kqvZU1VRVTW3evHnpo5SkNay31U1V9cjQ+xuTfCTJJgaVw2lDQ7cAh5c7Pknq0rirmfo6FzGjtySR5GnAD6uqkmxnUNU8BPwUOCvJGcD3gF3An/YVpyR1YZwE0ee5iBmdJYkk1wLnA5uSTANXAOsBquoq4BXAm5IcAX4O7KqqAo4kuQy4GVgH7K2qu7qKU5LUrsvVTZcco//DwIdb+m4EbuwiLknq00qZZprhFdeStIxWyjTTjL6XwEqSJpiVhCQtg5U2zTTDJCFJy2ClTTPNMElIUocm9e6u4/KchCR1aNwEMWnTTDOsJCSpA/OpICZxmmmGlYQkdWAlTzENs5KQpCU033MQkzrNNMMkIUlLaD4JYpKnmWY43SRJPZj0CmKGlYQkLYHVcqJ6NisJSVoCK32paxsrCUlaoPmepF5JFcQMk4QkzdNCrqJeaRXEDJOEJM3TfBPESqwgZpgkJGkMK/0eTAvliWtJGsNCE8RKnWaa0eUzrvcCLwUeqKpnjeh/FfAXzeZjwJuq6mtN33eAR4HfAEeqaqqrOCVpLgutIFbyFNOwLiuJq4Gdc/R/G3hhVT0b+I/Anln9F1TVNhOEpD4tJEGs9OphWGeVRFXdluT0Ofq/PLR5O7Clq1gkaT4WWj1s2riB/X/54g4i6s+knLh+A3DT0HYBtyQp4K+ranaV8bgku4HdAFu3bu00SElrw1pavXQsvSeJJBcwSBIvGGreUVWHkzwFuDXJPVV126j9mwSyB2Bqaqo6D1jSqrSY6mE16zVJJHk28HHgwqp6aKa9qg43fz6Q5HpgOzAySUjSYixmaetqriBm9LYENslW4DPAa6rqm0PtJyY5aeY98BLgYD9RSlrt1urS1nF1uQT2WuB8YFOSaeAKYD1AVV0FvAs4GfhIEvj/S12fClzftB0PXFNVn+8qTklrj9XD+Lpc3XTJMfrfCLxxRPt9wLldxSVJVg/j6/3EtSQtB6uHhfG2HJLWBKuHhbGSkLQqLfaGfKvxwriFOGaSSPIog4vbhj0M7Afe1pxDkKTeLcWdWtfy1NIo41QSHwAOA9cAAXYBTwPuBfYyWMEkSb1bbIJY61NLo4yTJHZW1XlD23uS3F5V703yzq4Ck6RxWD10a5wT1/+c5JVJjmterxzq8zYYknpl9dCtcSqJVwH/FfgIg6RwO/DqJE8ALuswNkk6ipXD8jpmkmhOTL+spfsfljYcSZqblcPyGmd10+8AHwWeWlXPam7K98dV9b7Oo5O05i1F5eBy1oUbZ7rpY8Dbgb8GqKqvJ7kGMElIWnJLkRRmOK20eOOcuH5iVf3jrLYjXQQjSUuVIJxWWhrjVBI/SnImzUqmJK8Avt9pVJLWDCuHyTZOkngzgye/PTPJ94BvA6/uNCpJa4aVw2Qbd3XTHzYPADquqh7tPixJq9FSVg1g5bAcWpNEkte2tANQVZ/oKCZJq9RSJggrh+UxVyXxb0e0hcE1E6cCJglJc1rKysFlrP1oTRJV9e9n3mdQPrwK+AsGV1z/p+5Dk7TSeWX0yjfnOYkkxwN/BrwN+Arwiqq6d5wPTrIXeCnwQFU9a0R/GNzu4yLgZ8CfVdVXm76dTd864ONV9f5x/0KS+rHU5xucTpoMc52TeDPwFuB/MbgT7Hfn+dlXAx+mfVrqQuCs5nUeg6u6z0uyDrgSeDEwDexLckNVfWOe3y9pGVk1rE5zVRIfAh4AXgD8z5kT1gzOS1RVPXuuD66q25KcPseQi4FPVFUBtyd5cpKnA6cDh2YeZpTkumasSUKaAEtdMYBVwySbK0mc0fF3nwrcP7Q93bSNah9+nsVvSbIb2A2wdevWpY9S0m/xwre1Za4T1/OdXpqvjGirOdpHqqo9DC72Y2pqyudbSEuki4phhpXDyjHOFdddmQZOG9rewuAxqRta2iUtIy96E4x3g7+u3AC8NgPPAx6uqu8D+4CzkpyRZAODZ2rf0GOckhbBqmFl66ySSHItcD6wKck0cAWwHqCqrgJuZLD89RCDJbCvb/qOJLkMuJnBEti9VXVXV3FKa5032NNcxnno0J0cfU7gYWA/8L6qemjUflV1yVyf26xqenNL340MkoikjnmDPc1lnEriJuA3wDXN9q7mz0cYXAvR9mhTSROirVpY6C92K4a1Y5wksaOqdgxt35nkS1W1I4m3DJdWgLZqYSFVhBXD2jJOktiY5Lyq+gpAku3AxqbPJ9RJq5gVg8ZJEm8E9ibZyOAahkeANzTPl/jPXQYnaXxLPaUkwXgPHdoH/F6SJwGpqp8OdX+yq8Akzc9Cp5Q2bdxgclGrcVY3PYnB8tU/aLa/CLy3qh7uODZJy8BnNGgu40w37QUOAq9stl8D/C3w77oKStLR5ppOOtYveqsFLdQ4SeLMqvqToe33JDnQUTySWixmhZLVghZqnNty/DzJC2Y2kuwAft5dSJKkSTFOJXEp8Inm3ATAT4DXdReStLYsZhpp9ninlLTUxlnd9DXg3CT/stl+JMlbga93HJu0JizVhW5OKakLY98FtqoeqapHms0/7ygeSdIEWeitwkc9GEhSh9qmjZxOUpcWeqtwnwAnHcNSnWuY4XSS+tCaJJI8yuhkEOAJnUUkrRJLeVM9qS9zPeP6pOUMRFqrXJWkSdbnM64l4TSSJptJQhrTUp9jkFaCha5uGkuSnUnuTXIoyeUj+t+e5EDzOpjkN0n+ddP3nSR3Nn37u4xTGofnGLQWdZYkkqwDrgQuBM4GLkly9vCYqvqrqtpWVduAdwBfrKofDw25oOmf6ipOqSsuWdVq0OV003bgUFXdB5DkOuBi4Bst4y8Bru0wHmlZOQWl1aDL6aZTgfuHtqebtqMkeSKwE/j0UHMBtyS5I8nuti9JsjvJ/iT7H3zwwSUIW5I0o8tKYtRV2W0X4b0M+NKsqaYdVXU4yVOAW5PcU1W3HfWBVXuAPQBTU1Ne5Kdj8gS0NL4uK4lp4LSh7S3A4Zaxu5g11VRVh5s/HwCuZzB9JS3aYh7zOZ92aTXospLYB5yV5AzgewwSwZ/OHtTcgvyFwKuH2k4EjquqR5v3LwHe22Gs0jFZZWgt6ixJVNWRJJcBNwPrgL1VdVeSS5v+q5qhLwduqap/Gtr9qcD1SWZivKaqPt9VrJKk0Tq9mK6qbgRunNV21aztq4GrZ7XdB5zbZWySpGPr9GI6SdLK5m05tGL4mE9p+ZkktGL4mE9p+TndJElqZZKQJLUySUiSWpkkJEmtTBJaMbwthrT8XN2kFcNVSdLyM0loyXmXVWn1cLpJS87HfEqrh0lCktTKJCFJamWSkCS1MklIklqZJLTkvJ5BWj1cAqsl5zJXafWwkpAkteo0SSTZmeTeJIeSXD6i//wkDyc50LzeNe6+kqTudTbdlGQdcCXwYmAa2Jfkhqr6xqyhf19VL13gvhrBK54lLZUuK4ntwKGquq+qfgVcB1y8DPuueV7xLGmpdJkkTgXuH9qebtpme36SryW5Kck589yXJLuT7E+y/8EHH1yKuCVJjS6TREa01aztrwLPqKpzgQ8Bn53HvoPGqj1VNVVVU5s3b15orJKkEbpMEtPAaUPbW4DDwwOq6pGqeqx5fyOwPsmmcfaVJHWvyySxDzgryRlJNgC7gBuGByR5WpI077c38Tw0zr6SpO51trqpqo4kuQy4GVgH7K2qu5Jc2vRfBbwCeFOSI8DPgV1VVcDIfbuKdbXZtHFD6+omSZqPDH4nrw5TU1O1f//+vsOQpBUjyR1VNdXW7xXXkqRWJglJUiuThCSplUlCktTKJCFJamWSkCS1MklIklqZJCRJrUwSkqRWa/oZ1z6cR5LmtqYrCR/OI0lzW9NJQpI0N5OEJKmVSUKS1MokIUlqtaaTRNtDeHw4jyQNrOklsC5zlaS5relKQpI0t06TRJKdSe5NcijJ5SP6X5Xk683ry0nOHer7TpI7kxxI4jNJJakHnU03JVkHXAm8GJgG9iW5oaq+MTTs28ALq+onSS4E9gDnDfVfUFU/6ipGSdLcuqwktgOHquq+qvoVcB1w8fCAqvpyVf2k2bwd2NJhPJKkeeoySZwK3D+0Pd20tXkDcNPQdgG3JLkjye4O4pMkHUOXq5syoq1GDkwuYJAkXjDUvKOqDid5CnBrknuq6rYR++4GdgNs3bp18VFLkh7XZSUxDZw2tL0FODx7UJJnAx8HLq6qh2baq+pw8+cDwPUMpq+OUlV7qmqqqqY2b968hOFLkrpMEvuAs5KckWQDsAu4YXhAkq3AZ4DXVNU3h9pPTHLSzHvgJcDBDmOVJI3Q2XRTVR1JchlwM7AO2FtVdyW5tOm/CngXcDLwkSQAR6pqCngqcH3TdjxwTVV9vqtYJUmjpWrkaYIVaWpqqvbv95IKSRpXkjuaf5yP5BXXkqRWJglJUiuThCSplUlCktTKJCFJamWSkCS1MklIklqZJCRJrUwSkqRWJglJUiuThCSplUlCktTKJCFJamWSkCS1MklIklqZJCRJrUwSkqRWJglJUiuThCSpVadJIsnOJPcmOZTk8hH9SfLBpv/rSZ4z7r6SpO51liSSrAOuBC4EzgYuSXL2rGEXAmc1r93AR+exrySpY11WEtuBQ1V1X1X9CrgOuHjWmIuBT9TA7cCTkzx9zH0lSR07vsPPPhW4f2h7GjhvjDGnjrkvAEl2M6hCAB5Lcu8C490E/GiB+3bN2BbG2BbG2BZukuNri+0Zc+3UZZLIiLYac8w4+w4aq/YAe+YX2tGS7K+qqcV+TheMbWGMbWGMbeEmOb6FxtZlkpgGThva3gIcHnPMhjH2lSR1rMtzEvuAs5KckWQDsAu4YdaYG4DXNqucngc8XFXfH3NfSVLHOqskqupIksuAm4F1wN6quivJpU3/VcCNwEXAIeBnwOvn2rerWBuLnrLqkLEtjLEtjLEt3CTHt6DYUjVyql+SJK+4liS1M0lIklqtqSSxmNuETEh85yd5OMmB5vWuZYprb5IHkhxs6e/tuI0RWy/HrPnu05J8IcndSe5K8pYRY3o5dmPG1tfP2wlJ/jHJ15rY3jNiTF/HbZzYevuZa75/XZL/m+RzI/rmf9yqak28GJwA/xbwbxgssf0acPasMRcBNzG4TuN5wFcmLL7zgc/1cOz+AHgOcLClv8/jdqzYejlmzXc/HXhO8/4k4JuT8jM3Zmx9/bwF2Ni8Xw98BXjehBy3cWLr7Weu+f4/B64ZFcNCjttaqiQWc5uQSYmvF1V1G/DjOYb0dtzGiK03VfX9qvpq8/5R4G4GdxMY1suxGzO2XjTH4rFmc33zmr3Cpq/jNk5svUmyBfgj4OMtQ+Z93NZSkmi7Bch8x3Rl3O9+flPq3pTknOUJ7Zj6PG7j6P2YJTkd+H0G//Ic1vuxmyM26OnYNVMmB4AHgFuramKO2xixQX8/c/8F+A/AP7f0z/u4raUksZjbhCyHcb77q8Azqupc4EPAZ7sOakx9Hrdj6f2YJdkIfBp4a1U9Mrt7xC7LduyOEVtvx66qflNV2xjcbWF7kmfNGtLbcRsjtl6OW5KXAg9U1R1zDRvRNudxW0tJYjG3CVkOx/zuqnpkptStqhuB9Uk2LVN8c+nzuM2p72OWZD2DX8J/V1WfGTGkt2N3rNj6PnbN9/4U+D/Azlldvf/MtcXW43HbAfxxku8wmK5+UZL/PmvMvI/bWkoSi7lNyETEl+RpSdK8387gv99DyxTfXPo8bnPq85g13/s3wN1V9YGWYb0cu3Fi6+vYJdmc5MnN+ycAfwjcM2tYX8ftmLH1ddyq6h1VtaWqTmfw++N/V9WrZw2b93Hr8gZ/E6UWcZuQCYrvFcCbkhwBfg7sqmbJQpeSXMtgxcamJNPAFQxO2PV+3MaIrZdj1tgBvAa4s5nDBngnsHUovr6O3Tix9XXsng78twwePnYc8Mmq+tyE/L86Tmx9/swdZbHHzdtySJJaraXpJknSPJkkJEmtTBKSpFYmCUlSK5OEJKmVSUKS1MokIUlqZZKQlliSx8YYc2KSDzdXvUoTyyQh9eNS4F8AL+g7EGkuJgmpA0lOz+Cpbx/L4AlmtzT3+pmxk8GDfg70E6E0HpOE1J2zgCur6hzgp8CfwOARmAzuz/Uc4Iu9RSeNwSQhdefbVXWgeX8HcHrz/iwGSeKeqvp1D3FJY1szd4GVevDLofe/AWammzYDv8OEPJ5WmouVhLT8TmHwsJ/jkvyrvoOR5mKSkJZRkuMZnIt4GnAVgwpDmlg+T0KS1MpKQpLUyiQhSWplkpAktTJJSJJamSQkSa1MEpKkViYJSVKr/wd/H5Welb+ehQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plots log neg\n",
    "LN_ex=np.zeros(len(Lsub_list)) # Log negativity\n",
    "\n",
    "Gmat=np.dot(V1,np.matrix(V1).H)\n",
    "\n",
    "for i_L in range(len(Lsub_list)):\n",
    "    Lsub=Lsub_list[i_L]\n",
    "    inds_A= np.arange(0,Lsub_list[i_L])\n",
    "    inds_AB= np.arange(0,2*Lsub_list[i_L])\n",
    "    LN_ex[i_L]= exact_log_neg_calc(Lsub,Lsub,Gmat[np.ix_(inds_AB,inds_AB)])\n",
    "    \n",
    "plt.figure()\n",
    "plt.plot(np.log(Lsub_list),LN_ex,'s')\n",
    "xplot=np.log(np.tan(pi*Lsub_list/N))\n",
    "yplot=(xplot-xplot[0])/4+ 0.47\n",
    "plt.plot(np.log(Lsub_list),yplot,label='CFT')\n",
    "plt.legend()\n",
    "plt.ylim([0,2])\n",
    "plt.xlabel('$\\ln \\ell$')\n",
    "plt.ylabel('Log Neg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}