freefermion_python.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def renyi_calc(r,GA,epsilon=1e-9):\n",
    "    \"\"\"computes von Renyi entropy of 1d free fermions\"\"\"\n",
    "    \"\"\"GA: correlation matrix\"\"\"\n",
    "   \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",
    "def EE_calc(GA,epsilon=1e-9):\n",
    "    \"\"\"computes von Neumann entanglement entropy of 1d free fermions\"\"\"\n",
    "    \"\"\"GA: correlation matrix\"\"\"\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",
    "\n",
    "\n",
    "def mut_inf_calc(Na,Nb,GA,epsilon=1e-9):\n",
    "    \"\"\"computes mutual information of 1d free fermions\"\"\"\n",
    "    \"\"\"GA: correlation matrix\"\"\"\n",
    "    Nt=Na+Nb\n",
    "    S12 = EE_calc(GA)\n",
    "    S1 = EE_calc(GA[:Na,:Na])\n",
    "    S2 = EE_calc(GA[Na:,Na:])\n",
    "    return S1+S2-S12\n",
    "\n",
    "def log_neg_calc(Na,Nb,Gm,epsilon=1e-9):\n",
    "    \"\"\"computes log negativity of 1d free fermions\"\"\"\n",
    "    \"\"\"Gm: covariance matrix of A and B\"\"\"\n",
    "    \"\"\"Na, Nb: number of complex fermions in A and B\"\"\"\n",
    "    Nt=Na+Nb\n",
    "\n",
    "    Gm_1= np.block([\n",
    "        [-Gm[:2*Na,:2*Na], -1j*Gm[:2*Na,2*Na:]],\n",
    "        [-1j*Gm[2*Na:,:2*Na], Gm[2*Na:,2*Na:]]\n",
    "    ])\n",
    "    Gm_2= np.block([\n",
    "        [-Gm[:2*Na,:2*Na], 1j*Gm[:2*Na,2*Na:]],\n",
    "        [1j*Gm[2*Na:,:2*Na], Gm[2*Na:,2*Na:]]\n",
    "    ])\n",
    "    Gx=np.eye(2*Nt)- np.dot(np.eye(2*Nt)+1j*Gm_2,np.dot(np.linalg.inv(np.eye(2*Nt)-np.dot(Gm_1,Gm_2)),np.eye(2*Nt)+1j*Gm_1))\n",
    "    Gx=(Gx+np.transpose(Gx).conj())/2 # make sure it is hermitian\n",
    "\n",
    "    chi0 =np.linalg.eigvalsh(Gx)\n",
    "    i0=np.nonzero(np.abs(chi0-1)>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(((1+nu+1j*0)/2)**0.5+((1-nu+1j*0)/2)**0.5))/2\n",
    "    \n",
    "    chi =np.linalg.eigvalsh(1j*Gm)\n",
    "    sA=np.sum(np.log(((1+chi)/2)**2+((1-chi)/2)**2))/4    \n",
    "\n",
    "    return np.real(eA+sA)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### SSH chain sanity checks\n",
    "\n",
    "N = 100\n",
    "Lsub_list=np.arange(2,int(N/2)+1)\n",
    "N_pt = int(N/2)\n",
    "r=0.5\n",
    "\n",
    "t=1.0\n",
    "BC = -1\n",
    "hop= np.diag(np.ones(N-1),1)\n",
    "hop[N-1,0]= t*BC\n",
    "H_t= -(hop+ np.matrix(hop).H)/2 \n",
    "energies, evecs= np.linalg.eigh(H_t)\n",
    "V1 =  evecs[:,:N_pt]\n",
    "Gmat=np.dot(V1,np.matrix(V1).H)\n",
    "\n",
    "Rr=np.zeros(len(Lsub_list)) # Renyi entropy\n",
    "LN=np.zeros(len(Lsub_list)) # log neg\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[i_L]=renyi_calc(r,Gmat[np.ix_(inds_A,inds_A)])\n",
    "    Γ = np.kron(np.eye(2*Lsub)-2*Gmat[np.ix_(inds_AB,inds_AB)],[[0,1],[0,0]]) \\\n",
    "        -np.kron(np.eye(2*Lsub)-2*Gmat[np.ix_(inds_AB,inds_AB)],[[0,0],[1,0]])\n",
    "    LN[i_L]=log_neg_calc(Lsub,Lsub,Γ)\n",
    "#     LN[i_L]=log_neg_calc(Lsub,N-Lsub,Gmat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEYCAYAAABMVQ1yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhoklEQVR4nO3de5RU5Znv8e9D2yMaELyAGBBpl5gAIgoNahTESxQZjBjkRMeQ4OUQL1noSWJ0YpKZZMVZ6zieHESJxomXmDBqjGAcg+dETTyAmmjDgJdg1ATUVoZucYniBWl4zh+1q6kuatetq2pX7f37rNXLrr137Xo3yNNPP+/N3B0REWl8faJugIiIVIYCuohITCigi4jEhAK6iEhMKKCLiMSEArqISEwooIuUwcy+Y2Y/Czk32cz+0ot7u5kdVn7rJKkU0KXumNl+ZrbUzD4ws9fM7B/yXDvXzHaY2VYze8/M1prZjGq30d3/xd0vDjm3wt0/U+02iGRTQJeqM7M9SnzLIuAT4EDgfOAWMxuT5/qn3b0fMBD4CXCvmQ0so6kiDU0BXarCzDaY2dVm9hzwQbFB3cw+BcwCvufuW919JfAQMKfQe919J/AL4FPAyOB+e5rZDWb2upltMrNbzWyv4NxUM2s3s2+aWYeZbTSzC4JzE4Pru9ttZrPMbE3w/T+b2S9DnmGqmbVn/Vl8y8yeM7MtZnafmfXNOH9V8NlvmdmFWffK1/5lZva/Mq69z8zuKPTnJPGlgC7VdB7w98BAd+8ys4fN7N2Qr4eD9xwO7HD3lzPusxbIl6EDYGZNwAXAduC14PD/DO55FHAYMBT4fsbbhgADguMXAYvMbF93fxbYDHw+49ovk/qBUY7/BkwDWoAjgblBm6cB3wo+ZyRwatb78rX/QmCOmZ1sZucDE4ErymyfxECpvwqLlGKhu7+RfuHuxdS2+wFbso5tAfrnec+xZvYuqcy8C/iyu3eYmQH/HTjS3d8BMLN/Af4d+MfgvduBH7p7F7DMzLYCnwH+CPycVBB/xMz2A04HLiviGXJZ6O5vBW34D1IBGlKB/k53fyE498+kfhBSqP3u/l9mdknQzr2Ame7+fpntkxhQhi7V9EbhS3azFdgn69g+QL5A9Ud3HwjsS6o8Mzk4PgjYG1iV/k0A+D/B8bTNQTBP+5DUDxWAXwJnmlk/UoF3hbtvLPmJUv4r5DM+Tc8/p9cyvi+m/Q8DTcBfgvKUJJgCulRTj6U8zeyRYDRKrq9HgsteBvYws5EZbx0HvFjww9y3ksqg55jZ0cDbwEfAGHcfGHwNCDpQCzfe/U3gaeBsUjX8csst+WwEDs54PTzj+2Lafx2wDjjIzM6rQvukgSigS824+xnu3i/k64zgmg+AJcAPzexTZnY8cBZFBlN33wz8DPh+0En6b8D/NrPBAGY21MxOL6HZdwPfBsYCS0t4X7F+Bcw1s9FmtjfwT+kThdpvZlNI9Rl8Jfi6ycyGVqGN0iAU0KUeXUaqJtwB3ANc6u4FM/QMC4DpZnYkcDXwKvBHM3sPeIxUjbxYS4FDgKXBD5uKcvdHSLX390E7f591Sc72m9k+pH7YfN3d3wzKLbcDdwa1d0kg0wYXIvmZ2V+Br7n7Y1G3RSQfZegieZjZLFJ9AdmZs0jd0bBFkRBm9gQwGpgT1LNF6ppKLiIiMaGSi4hITERWcjnggAN8xIgRUX28iEhDWrVq1dvuPijXucgC+ogRI2hra4vq40VEGpKZvRZ2TiUXEZGYUEAXEYkJBXQRkZioq3Ho27dvp729nY8//jjqptSlvn37MmzYMJqbm6NuiojUoboK6O3t7fTv358RI0ag5Sh6cnc2b95Me3s7LS0tUTdHROpQXQX0jz/+WME8hJmx//7709nZGXVTRBJj+o0r+PPG93Y7vucefdjWtfvk4bDjow/ah3M+/jWPbhnK0zt3bb51XJ8X+fyAN7nw6hsr0t66CuiAgnke+rMRSQkLtKMPSu2NUqkgPH74QF7peJ/tO3bNqG9uMobvtzcbNn9Q9PHxh+yLvz+eRR9ey+Xb5/P0zjEc1+dFFjUvZOmnryvt4fOou4AuIslTaia8797NNDdZzsCJe8WC8PyTD+P+Ve1k7tXSZMaN5x7F2T95qujj8085DPwwrvzXTm5uXsAvd5zKl5se43/svJJ/nfmlUv+4QimgZ2lqamLs2LF0dXXR0tLCL37xCwYOHFix+7e1tXH33XezcOHCit1TpJ6UGpzLyYRP+uxgfvvcRnIHTioWhAf378vsCcO4r+0Ntu9wmpuMc1oPZvSnB5R0fHD/vgAMH38696x+kSv2WMLNO77IwROmdZ+rhIYN6Pl+5Vp2xeQc7yjOXnvtxZo1awD46le/yqJFi7j22mvLvl+21tZWWltbK3Y/kWoqp7QRlj1XMhO+5ozPsndzU2jgrGQQnn/KyO62df/QKOM4wDdHbsLXPMqNXWczp+lRbOT5wBGl/JXk1bDj0McPH0hzU8+acvevXBVy3HHH8eabbwLw17/+lWnTpjFhwgQmT57MSy+9BMDcuXOZP38+n/vc5zj00EP59a9/DcCcOXP4zW9+032v888/n4ceeognnniCGTNmVKyNIqWYfuMKRlzz292+PvPdR3Ie37jlo9B/Z2H/Bk/67GD6ZPX3pINzruPzTzmMwfukMuH0/bKDbfbxwf37Mv+Ukd33yw6cYedKPQ50t82MHoG+1OOsX86+y+ax9LDrWLBjNksPu459fzsP1i8v7i+vCA0b0DP/AtKy/yJ6Y8eOHTz++ON84QtfAGDevHncdNNNrFq1ihtuuIHLLrus+9qNGzeycuVKHn74Ya655hoALr74Yu68804AtmzZwlNPPcX06dMr0jaRfMKC9vQbV4QG4eH77V1ScJ5/ymGh/wavOeOzJQdnKD3YhgbOPOdKDsKB+aeMZOKI/XaLLyUdf3M1zL6LM8/6EhNH7MeZM78Es+9KHa+Qhi25pP8Cwn5NKtdHH33EUUcdxYYNG5gwYQKf//zn2bp1K0899RSzZ8/uvm7btm3d38+cOZM+ffowevRoNm3aBMCJJ57I5ZdfTkdHB0uWLGHWrFnssUfD/nFLHQorh+TrMKxFaSOdPZdajkj/m178zOs5g2328fT9Xu7YmjORCztX6vF0G371teN6d/yEK1PnYNe5/lOgZcpu7y9Xw2bokP/XpHKla+ivvfYan3zyCYsWLWLnzp0MHDiQNWvWdH+tW7eu+z177rln9/eZG4bMmTOHxYsXc+edd3LBBRf0um2STGEZd1g5JF9WXavSRq0y4XTgzJXIhZ0r9XgjaeiAXuh/jt4YMGAACxcu5IYbbmCvvfaipaWF+++/H0gF7bVr1xa8x9y5c1mwYAEAY8aMyX+xJFo5ZZKwwB1W8qhlaSN9v1KCc/p+cQ22tdDwNYB8vyb11tFHH824ceO49957Wbx4MZdeeik/+tGP2L59O+eeey7jxo3L+/4DDzyQUaNGMXPmzIq3TRpTJcsk+coh9VDaKLVMIb0X2Z6ira2tnr3Bxbp16xg1alQk7amGDz/8kLFjx7J69WoGDBhQkXvG7c8orvIF7q3bunYL3GeO+zS/fW5jj3Hafffow/KrT2Jw/758d+nzPQL3lyYO50czj6DjvY+ZfP0f2Na1s8f1AN9d+jyLn3md8485hB/N7Dk0ruO9j/n6Pf/Jzf9wdI/AHXZc6oeZrXL3nGOfGz5Dr1ePPfYYF154Id/4xjcqFsylvuQbox02USZsQkyhzsewjLvcrFrZczwpoFfJqaeeyuuvvx51M6QC6rlMAqWXPCS+GrpTVKQWSu2UzDeaJN+okXJHgKjDUNKUoYsEwjLxkYP7BQG4PsskImkFA7qZHQzcDQwBdgK3ufuNWddMBX4DrA8OLXH3H1a0pSIVUE7d+5hD9+eYlv1UJpG6V0yG3gV8091Xm1l/YJWZPeruf866boW7a5ESqQuVrHtnr+BXbLYNCtxSWwVr6O6+0d1XB9+/D6wDhla7YVHp169f1E2QCqh03bucCTSg+rbUVkk1dDMbARwN/CnH6ePMbC3wFvAtd38xx/vnAfMAhg8fXnJje1i5AIaO77kOwvrlqYVugjUTJP5qVfdOn1O2LfWs6FEuZtYPeAC40t2z/wWtBg5x93HATcCDue7h7re5e6u7tw4aNKjMJgeGjof75+5aenL98tTroeN7d9+Au3PVVVdxxBFHMHbsWO677z4Adu7cyWWXXcaYMWOYMWMG06dP714yN9PUqVO5+uqrmTRpEocffjgrVqwAUqs4XnXVVUycOJEjjzySn/70pyXdV3oKy8SPOXT/stYmKTTNXdm21LOiMnQzayYVzBe7+5Ls85kB3t2XmdlPzOwAd3+7ck3N0jIltfTk/XOh9SJouz31ukIrly1ZsoQ1a9awdu1a3n77bSZOnMiUKVN48skn2bBhA88//zwdHR2MGjWKCy+8MOc9urq6eOaZZ1i2bBk/+MEPeOyxx7j99tsZMGAAzz77LNu2beP444/ntNNOY9WqVUXfN2nydWTedcHEmtW9RepdMaNcDLgdWOfuPw65Zgiwyd3dzCaRyvw3V7SlubRMSQXz5dfDlG9XdBnKlStXct5559HU1MSBBx7IiSeeyLPPPsvKlSuZPXs2ffr0YciQIZx00kmh9/jiF78IwIQJE9iwYQMAv/vd73juuee6s+8tW7bwyiuvlHTfuCqnI7PQMsoaHihJUkzJ5XhgDnCyma0Jvqab2SVmdklwzTnAC0ENfSFwrtdikZj1y1OZ+ZRvp/5bwZ0/wppfymOll9Vtamqiq6ur+/033XRT9zK869ev57TTTivpvnFVTkcm5F9GWZNxJEmKGeWy0t3N3Y9096OCr2Xufqu73xpcc7O7j3H3ce5+rLs/VfWWp2vms++Ck6/dVX6pUFCfMmUK9913Hzt27KCzs5Ply5czadIkTjjhBB544AF27tzJpk2beOKJJ0q67+mnn84tt9zC9u3bAXj55Zf54IMPen3fRhK2VOyf1r9T1nKwqnuLpDTuTNFgO6fuMku6pv7m6oqUXs4++2yefvppxo0bh5lx/fXXM2TIEGbNmsXjjz/OEUccweGHH84xxxxT0uJbF198MRs2bGD8+PG4O4MGDeLBBx/s9X0bSaUn8IDq3iKg5XPLsnXrVvr168fmzZuZNGkSTz75JEOGDKnJfRvlz6hQR2Z6yde09NKvOGUtByuSFFo+t8JmzJjBu+++yyeffML3vve9igTzat43CmFZeLU6MkVEAb0s1apvx6lunlkiSSt2s2BN4BEpT90FdHfHsjrGJKUeR8LkK63ky8LzjQVX4BYpT12th963b182b95cl4Erau7O5s2b6du3vkZrhA01HH/IvnmHE0L+NVBEpHR1laEPGzaM9vZ2Ojs7o25KXerbty/Dhg2L5LNLXTNl/imHMbh//hmZysRFKquuAnpzczMtLS1RN0NyKGeoIagjU6SW6iqgS7QqvWYKKAsXqaW6qqFLtPLVw8tdK1xEakcBXbpldmKm9XbNFBGpHZVcEkhDDUXiSRl6AmmooUg8KaAnUL7SSqF6uFYvFKlfKrnEWLmlFQ01FGlMytBjrNzSirJwkcakgB5jvSmtiEjjUcmlweUrqyy7YrJKKyIJogy9weUrq0D+seMqrYjEiwJ6gys0GUilFZHkUEBvcPmm5Kdp7LhIMqiG3iCKXTgr12QgzeIUSQZl6A2imIWzVFYRSTYF9AZRzMJZKquIJJsCeoMoVCvXiBURUQ29jhQaUz7/lJF5a+UikmzK0OtIoTHlqpWLSD4K6HWkUJ08fY1q5SKSiwJ6HSlmTLlq5SISRjX0CPRmTLmISBhl6BHQmHIRqQYF9AhoTLlITK1cAOuX9zy2fnnqeA0ooEdAY8pFYmroeLh/7q6gvn556vXQ8TX5eAX0iBTajFlEGlDLFJh9VyqI//661H9n35U6XgPqFK2iYjafWPzM66qVi8RJyxRovQiWXw9Tvl2zYA7K0KuqmM0nVCsXiZn1y6Ht9lQwb7t995p6FRUM6GZ2sJn9wczWmdmLZnZFjmvMzBaa2atm9pyZ1aZgVOeK2XxCtXKRGEnXzGffBSdfu6v8UqOgXkyG3gV8091HAccCl5vZ6KxrzgBGBl/zgFsq2soGVcxEIRGJkTdX96yZp2vqb66uyccXDOjuvtHdVwffvw+sA4ZmXXYWcLen/BEYaGYHVby1DUidnyIJcsKVu9fMW6akjtdASZ2iZjYCOBr4U9apocAbGa/bg2Mbe9O4RlCo4zOdpavzU0SqrehOUTPrBzwAXOnu2RHMcrzFc9xjnpm1mVlbZ2dnaS2tU4U6PkGdnyJSG0UFdDNrJhXMF7v7khyXtAMHZ7weBryVfZG73+bure7eOmjQoHLaW3eKWSFRnZ8iUgvFjHIx4HZgnbv/OOSyh4CvBKNdjgW2uHvsyy2gjk8RqR/F1NCPB+YAz5vZmuDYd4DhAO5+K7AMmA68CnwIXFDxltYx7SQkIvWgYEB395XkrpFnXuPA5ZVqVKNRx6eI1ANN/S9SMft9vtyxVdm5iERGU/+LVMx+n+r4FJEoKaAXqZjRLCIiUVJAL5JGs4hIt4g3sgijgF4CTeMXESDyjSzCKKCXQPt9iggQ+UYWYTTKpUQazSIiQKQbWYRRhl4ijWYRESDSjSzCKEPPUmi8uYhIj40sWqZAy+S6KLsoQ89SzOqJIpJwEW9kEUYZepbMdVnSNKJFRHrItWFFy5TI6+jK0LNovLmINCoF9Bw03lxEgLqdQBRGAT0HjTcXEaBuJxCFUQ09hMabi0iPCUStF6WGJ9bBBKIwCugh0uPNRSTh6nACURiVXERE8qnDCURhEpmha/KQiBSlTicQhUlkhq7JQyJSlDqdQBQmkRm6Jg+JSLeVC1KjVjIz7vXLU0G7TicQhUlkhq7JQyLSrcGGJuaTyIAOmjwkIoE6Xdu8HIkN6Jo8JCLdMocmtl7UkMEcEhzQIZWlTxyxn7JzkaRroKGJ+SSyUzRNk4dEpNGGJuaT6AxdRBIkbKGtJxc21NDEfBTQRSQZwkazHD9/90y8ZUruIYt1LtElFxFJkAZbaKscsQ7omuIvIj000EJb5Yh1yUVT/EUSKN+mFDEZzRIm1gE9c/JQmiYRicRcWK28zx67Rq+cfO2u8kuMgnqsA7qm+IskUNjMz51dsRnNEibWNXTouRCXsnORhMhVK89VL6/jhbbKEesMHTTFXySWCm3eHPNaeZjYB3TQFH+R2Mm3QmLmzM+Y1srDmLsXvqoKWltbva2tLZLPFpEYSAfu7DHlpa5v3mDMbJW7t+Y6F/sauog0sELBOdeY8gbblKKSCpZczOwOM+swsxdCzk81sy1mtib4+n7lmykiiVSotJLAOnk+xWTodwE3A3fnuWaFu8+oSItERNLCputDbFZIrKSCGbq7LwfeqUFbRCSp8o1aybX5RINt3lwrlRrlcpyZrTWzR8xsTNhFZjbPzNrMrK2zs7NCHy0iDa/U0soJV8ZmhcRKqkRAXw0c4u7jgJuAB8MudPfb3L3V3VsHDRpUgY8WkYZRKAvPNbsTEjsEsRy9Duju/p67bw2+XwY0m9kBvW6ZiMRLviwcVFqpgF4PWzSzIcAmd3czm0Tqh8TmXrdMROKl0Hrk2aWVlsmJHoJYjoIB3czuAaYCB5hZO/BPQDOAu98KnANcamZdwEfAuR7VbCURiV45Y8djtK9nlAoGdHc/r8D5m0kNa4yMNrIQqSPp0ko6GGcG61xZeKHSigJ60WIxU3T88IG80vE+23fs+sVAG1mIVFm+TLzUseMqrVRELBbn0kYWIhHI18mpDs5IxCJDTy+Re1/bG2zf4drIQqRSCtXDwzo51cEZiVhk6NAzS1d2LlIh5Qw1TPDytVGLTUDXRhYivRA26SezHp454SdsqGHme1RaqbnYBHTQRhYiZSu1Hp55TXYmnl2iAU3Lr5FYBfTB+/TlV187Ttm5SJhyMvGwZWqVidedWAV0ESH/mimlZuL56uFaIKvuKKCLxE2hoF1KJq4svKHEYtiiSCKVOrEnMyhnT78vZeq9hhrWLWXoIo2qnI7M9HXKxGNJGbpIvSsnEw9bM0WZeKwpQxepB7XqyFQmHmsK6CL1oFYdmRqZEmsquYjUkjoypYqUoYtUQ1gJ5Z316siUqlFAF6mGsBLK2Fmlz8jUFHspkkouIuUqd2lZKK18ot18pEjK0EUKKad8AvkXtVJHplSBMnSRQvLtkTl2Vmm72IM6MqVqFNBF0sodgVLKLvajzlT5RKpGJRdJlkpO4ElfU8rSsvu2qHwiVaMMXeIpLNt+Zz08tTB3+SRzAk8xU+lBu9hLXWmYgD79xhX8eeN7ux0ffdA+LLticgQtkrpWbt27lAk8Kp9InWmYgD5++EBe6Xif7Tu8+1hzkzH+kH0jbJVErpJ17/R7szPxUoYNKguXCDVMDX3+KSPpY9bjWJOZ9g9NglrVvTWBRxpcwwT0wfv0ZfaEYTQ3pYJ6c5NxTuvB2j80CSq5cJVWIpQYM3cvfFUVtLa2eltbW0nv6XjvYyZf/we2de2k7x59WH71SQrocZKvfJIO6rnKJ5AK5ukSysnX7l73Tr8edSYcMSt8dqdInTOzVe7emutcw2TosCtLN0PZeSOLcuEqDRuUGGuogA6pWvrEEfupdt7ItHCVSFU0VMlFGkihhavSgThXCaXY8km6vp3vc0RiJjYlF6lDWrhKpG40zDh0qVNauEqkbiigS3G0cJVI3VPJRYqjhatE6p4ydNml3B14tHCVSF0omKGb2R1m1mFmL4ScNzNbaGavmtlzZja+8s2UiqpkR2bYsMEXHtCsS5EaK6bkchcwLc/5M4CRwdc84JbeN0uqqpxx4OnrNIFHpG4VLLm4+3IzG5HnkrOAuz01oP2PZjbQzA5y942VaqSUqRYdmRqBIlI3KtEpOhR4I+N1e3BsN2Y2z8zazKyts7OzAh8tedWiI1MlFJG6UYlOUctxLOf0U3e/DbgNUjNFK/DZoo5MEQlUIkNvBw7OeD0MeKsC95Vi5MvCQR2ZIglSiYD+EPCVYLTLscAW1c+rIGxkSmY9XB2ZIolWzLDFe4Cngc+YWbuZXWRml5jZJcEly4C/Aa8C/wZcVrXWJlm59XCtRCiSGMWMcjmvwHkHLq9Yi5KskvXwliml7YUpIg1PU//rSSXr4euXayVCkYRRQI9CLerh6sgUSRwF9CioHi4iVaCAHoVSd6oHZeIiUpBWW6ymQp2c2VPs802v18QeESlAGXo15SutqB4uIhWmDL2aMksrmUMNQQtdiUjFKUPvrbARKysXpL7P1cmpTFxEqkABvbcKjR3PVVrR+HARqQKVXHorrKxSqJNTZRQRqTBl6MXKV1oJGzuu0oqI1JACerFKHbECKq2ISE2p5FKsSoxYERGpImXo2UotraisIiJ1QgE9W6mlFZVVRKROqOSSTaUVEWlQytBzUWlFRBpQMjP0Qotm5doBSItjiUidS2aGXqhOHrYDkIhIHUtmQM+3HrlKKyLSoJJZcoGedfL0euSg0oqINKxkZugQPrtTRKRBxTtDD+v8fP4BeOk/NARRRGIl3hl6WOenoTq5iMROvDP0fEvb5rpW2bmINLB4Z+gQvrStiEjMxD+gq/NTRBIi3gFdk4REJEHiEdDDlrx9cqE6P0UkMeIR0MNGsxw/X0vbikhixGOUSymjWUREYioeGTpoNIuIJF58ArpGs4hIwsUjoGs0i4hITAK6lrwVEYlJp6iWvBURKS5DN7NpZvYXM3vVzK7JcX6qmW0xszXB1/cr31QREcmnYEA3syZgEXAGMBo4z8xG57h0hbsfFXz9sMLtDJ88tHJBxT9KRKQRFZOhTwJedfe/ufsnwL3AWdVtVg759gEVEZGiAvpQ4I2M1+3BsWzHmdlaM3vEzMZUpHWZ8u0DKiIiRQV0y3HMs16vBg5x93HATcCDOW9kNs/M2sysrbOzs6SGApo8JCKSRzEBvR04OOP1MOCtzAvc/T133xp8vwxoNrMDsm/k7re5e6u7tw4aNKj01mrykIhIqGIC+rPASDNrMbO/A84FHsq8wMyGmJkF308K7ru5oi3V5CERkbwKjkN39y4z+zrwf4Em4A53f9HMLgnO3wqcA1xqZl3AR8C57p5dlumdfJOHVHoREcEqHXeL1dra6m1tbZF8tohIozKzVe7emutcPKb+i4iIArqISFwooIuIxIQCuohITCigi4jERGSjXMysE3gteHkA8HYkDYlekp8dkv38evbk6s3zH+LuOWdmRhbQezTCrC1sGE7cJfnZIdnPr2dP5rND9Z5fJRcRkZhQQBcRiYl6Cei3Rd2ACCX52SHZz69nT66qPH9d1NBFRKT36iVDFxGRXlJAFxGJiUgDuplNM7O/mNmrZnZNlG2pBTO7w8w6zOyFjGP7mdmjZvZK8N99o2xjtZjZwWb2BzNbZ2YvmtkVwfHYP7+Z9TWzZ4ItGl80sx8Ex2P/7Glm1mRm/2lmDwevk/TsG8zseTNbY2ZtwbGqPH9kAd3MmoBFwBnAaOA8MxsdVXtq5C5gWtaxa4DH3X0k8HjwOo66gG+6+yjgWODy4O87Cc+/DTg52KLxKGCamR1LMp497QpgXcbrJD07wEnuflTG2POqPH+UGfok4FV3/5u7fwLcC5wVYXuqzt2XA+9kHT4L+Hnw/c+BmbVsU624+0Z3Xx18/z6pf9xDScDze8rW4GVz8OUk4NkBzGwY8PfAzzIOJ+LZ86jK80cZ0IcCb2S8bg+OJc2B7r4RUkEPGBxxe6rOzEYARwN/IiHPH5Qc1gAdwKPunphnBxYA3wZ2ZhxLyrND6of378xslZnNC45V5fkLbkFXRZbjmMZQxpyZ9QMeAK509/eCrWhjz913AEeZ2UBgqZkdEXGTasLMZgAd7r7KzKZG3JyoHO/ub5nZYOBRM3upWh8UZYbeDhyc8XoY8FZEbYnSJjM7CCD4b0fE7akaM2smFcwXu/uS4HBinh/A3d8FniDVl5KEZz8e+IKZbSBVVj3ZzH5JMp4dAHd/K/hvB7CUVLm5Ks8fZUB/FhhpZi1m9nfAucBDEbYnKg8BXw2+/yrwmwjbUjWWSsVvB9a5+48zTsX++c1sUJCZY2Z7AacCL5GAZ3f3f3T3Ye4+gtS/8d+7+5dJwLMDmNmnzKx/+nvgNOAFqvT8kc4UNbPppOprTcAd7n5dZI2pATO7B5hKaunMTcA/AQ8CvwKGA68Ds909u+O04ZnZCcAK4Hl21VK/Q6qOHuvnN7MjSXV8NZFKon7l7j80s/2J+bNnCkou33L3GUl5djM7lFRWDqkS97+7+3XVen5N/RcRiQnNFBURiQkFdBGRmFBAFxGJCQV0EZGYUEAXEYkJBXQRkZhQQBcRiQkFdJEswey+m4MlbkUahgK6yO4uAfYEToi6ISKlUEAX2d004GVgTcTtECmJArpIBjPrS2rNlfHA/4u4OSIlUUAX6WkkqYD+krtvj7oxIqWIcoMLkXo0CDicmG+HKPGkDF2kp0+T2oSjT5x3opd4UkAXCZjZHqRq50OAW4Ed0bZIpDRaD11EJCaUoYuIxIQCuohITCigi4jEhAK6iEhMKKCLiMSEArqISEwooIuIxMT/B9Z9Rn0gq+9NAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "#     plt.plot(Lsub_list,-Rr_ex,'o',label='(-1) Renyi Entropy')\n",
    "plt.plot(Lsub_list,(Rr),'v',label='Renyi')\n",
    "plt.plot(Lsub_list,(LN),'x',label='log neg')\n",
    "plt.title('r=%d Renyi index' % (r))\n",
    "plt.xlabel('$\\ell$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Kitaev chain Hamiltonian\n",
    "\n",
    "def BdG_Hamiltonian(N,t,Δ,μ, Tm,BC):\n",
    "    hopping = np.diag(np.ones(N-1),-1)\n",
    "    hopping[0,N-1]=BC\n",
    "    H = np.block([[-t*(hopping+hopping.T)-μ*np.eye(N), Δ*(hopping-hopping.T)],\\\n",
    "        [-Δ*(hopping-hopping.T),t*(hopping+hopping.T)+μ*np.eye(N)]])\n",
    "    \n",
    "    En, Vn = np.linalg.eigh(H/2)\n",
    "    nF = np.diag(1/(1+np.exp(En/Tm)))\n",
    "#     Uvec = np.matrix(Vn[:N,:N])\n",
    "#     Vvec = np.matrix(Vn[N:,:N])\n",
    "    Uvec = np.matrix(Vn[:N,:])\n",
    "    Vvec = np.matrix(Vn[N:,:])\n",
    "    Fmat = Uvec.dot(nF.dot(Vvec.H))\n",
    "    Cmat = Vvec.dot(nF.dot(Vvec.H))\n",
    "    Γ = 1j*np.kron(Fmat.H+Fmat+Cmat-Cmat.T,[[1,0],[0,0]])\\\n",
    "        -1j*np.kron(Fmat.H+Fmat-Cmat+Cmat.T,[[0,0],[0,1]])\\\n",
    "        +np.kron(np.eye(N)+Fmat.H-Fmat-Cmat-Cmat.T,[[0,1],[0,0]])\\\n",
    "        -np.kron(np.eye(N)+Fmat.H-Fmat-Cmat-Cmat.T,[[0,1],[0,0]]).T\n",
    "    return Γ\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adjacent intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "δ = 0.\n",
    "Tm = 5e-3\n",
    "# from Haining's notes\n",
    "t = -(1+δ)/2\n",
    "Δ = t\n",
    "μ = (1-δ)\n",
    "Nc = 100\n",
    "N = 2*Nc # number of majoranas\n",
    "BC = -1\n",
    "\n",
    "Gm = BdG_Hamiltonian(Nc,t,Δ,μ,Tm,BC)\n",
    "Lsub_list=np.arange(2,int(N/3),2)\n",
    "LN=np.zeros(len(Lsub_list)) # log neg\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",
    "    LN[i_L]=log_neg_calc(int(Lsub/2),int(Lsub/2),Gm[np.ix_(inds_AB,inds_AB)])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVkklEQVR4nO3df5BV9Znn8fcjaJioQzJ0b3BoHRj8HYOCLRi0NCETQlxrjDGpKMVmdt2URWUYk3XLUWsrW5vNH1ub2qQyCZbEis5sFXE0hcaxjDWQ6KhRUoQmIJEASiuDPWZCw5ZMacUI+Owf90IuzW36Xujm3vvl/arq6j7f8z3nPg8/Phy+5/bpyEwkSeU6qdUFSJLGlkEvSYUz6CWpcAa9JBXOoJekwo1vdQH1dHV15dSpU1tdhiR1jHXr1u3KzO56+9oy6KdOnUpfX1+ry5CkjhER/zzcPpduJKlwBr0kFc6gl6TCGfSSVDiDXpIKZ9B3mGXP9LO6f9chY6v7d7Hsmf6m543muVo1r51rs9cye2hVr8fCoB9jo/2bPaNnIkseWH9w7ur+XSx5YD0zeiY2PW80z9Wqee1cm72W2UOrej0mmdl2H5deemmW4vltgznzf67K57cN1t1udl7tvm+s3DLsnEbnjea5WjWvnWuz1zJ7aFWvRwL05TCZ2vJQr/dRUtBnjs1v9jdWbsk/uePx/MbKLUd87Ubmjea5WjWvnWsb7XntXFuj89q5ttGeN9qvOZwjBb1LN8fB3OldLJpzFt9+ahuL5pzF3OldxzRvdf8ulq/Zwa3zzmb5mh2HLfk0M280z9Wqee1cm72W2UOrej1qw/0L0MqPTriiv+fpbXWXX+55etthc8fiv4KjsRQ02stKrZjXzrXZa5k9tKrXkXCsSzfAAmArsA248wjzLgP2A5+pbp8J/BOwGdgEfKmR1+uEoG/VH7BG/4FpZN5onqtV89q5Nnsts4dW9TqSIwV9VPYPLyLGAS8BHwcGgLXATZn5qzrzfgy8DdyfmSsi4gzgjMz8RUScDqwDPjX02KF6e3uzEx5qduDu+KI5Z7F8zQ6WLpx52HLLsmf6mdEz8ZDx1f272Diwh8VXT296niTVExHrMrO37r4Ggv7DwP/IzE9Ut+8CyMz/NWTel4G9VK7qH8/MFXXO9Q/A0sz88ZFes1OCHuCbq7by7ae2ceu8s7lt/nmtLkfSCepIQd/IzdgpwGs12wPVsdoXmAJcDyw7QhFTgZnAmmH23xIRfRHRNzg42EBZrTfmN1AkaRQ0EvRRZ2zofwO+BdyRmfvrniDiNOBh4MuZ+W/15mTmvZnZm5m93d11n53fVg4s2yxdOJPb5p/H0oUzD/mmB0lqF4384JEBKjdVD+gBXh8ypxd4MCIAuoBrImJfZj4aESdTCfnvZ+Yjo1BzW9g4sOeQNfm507tYunAmGwf2DPu2SElqhUbW6MdTuRn7MeBfqNyMXZiZm4aZ/3dU1+ijkvz/F/h/mfnlRovqpDV6SWoHx7RGn5n7gCXASipvk/xBZm6KiMURsXiEw68A/gMwLyI2VD+uabJ+SdIxaOhnxmbmE8ATQ8bq3njNzP9Y8/Vz1F/jlyQdJz4CQZIKZ9APcTyeDS1Jx5NBP8RxeTa0JB1HDa3Rn0gOvE1ypEcbSFKn8Iq+jkYfFyxJncCgr8NHG0gqiUE/hI82kFQag36IIz3aQJI60YiPQGgFH4EgSc051scUS5I6mEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqXENBHxELImJrRGyLiDuPMO+yiNgfEZ9p9lhJ0tgYMegjYhxwN/BJ4ELgpoi4cJh5/xtY2eyxkqSx08gV/WxgW2a+kpnvAA8C19WZ91fAw8DOozhWkjRGGgn6KcBrNdsD1bGDImIKcD2wrNlja85xS0T0RUTf4OBgA2VJkhrRSNBHnbEcsv0t4I7M3H8Ux1YGM+/NzN7M7O3u7m6gLElSI8Y3MGcAOLNmuwd4fcicXuDBiADoAq6JiH0NHitJGkONXNGvBc6JiGkRcQpwI/BY7YTMnJaZUzNzKrAC+GJmPtrIscfTsmf6Wd2/65Cx1f27WPZMf4sqkqSxN2LQZ+Y+YAmVd9NsBn6QmZsiYnFELD6aY4+97KMzo2ciSx5YfzDsV/fvYskD65nRM7FVJUnSmIvMukvmLdXb25t9fX1jcu4D4b5ozlksX7ODpQtnMnd615i8liQdLxGxLjN76+074b4zdu70LhbNOYtvP7WNRXPOMuQlFe+EC/rV/btYvmYHt847m+Vrdhy2Zi9JpTmhgv7Ass3ShTO5bf55LF0485A1e0kq0QkV9BsH9hyyJj93ehdLF85k48CeFlcmSWPnhLsZK0kl8masJJ3ADHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMI1FPQRsSAitkbEtoi4s87+6yJiY0RsiIi+iLiyZt9/iYhNEfFiRPx9REwYzQYkSUc2YtBHxDjgbuCTwIXATRFx4ZBpTwIXZ+YlwM3A96rHTgFuBXoz8yJgHHDjqFUvSRpRI1f0s4FtmflKZr4DPAhcVzshM9/MzKxungpkze7xwB9ExHjgvcDrx162JKlRjQT9FOC1mu2B6tghIuL6iNgC/IjKVT2Z+S/A/wF2AL8G9mTmqnovEhG3VJd9+gYHB5vrQpI0rEaCPuqM5WEDmT/MzPOBTwFfA4iI91O5+p8G/DFwakQsqvcimXlvZvZmZm93d3eD5UuSRtJI0A8AZ9Zs93CE5ZfMfBaYHhFdwJ8Br2bmYGbuBR4B5h5DvZKkJjUS9GuBcyJiWkScQuVm6mO1EyLi7IiI6tezgFOA3VSWbC6PiPdW938M2DyaDUiSjmz8SBMyc19ELAFWUnnXzP2ZuSkiFlf3LwNuAD4fEXuB3wKfq96cXRMRK4BfAPuA9cC9Y9OKJKme+P2bZdpHb29v9vX1tboMSeoYEbEuM3vr7fM7YyWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVrqGgj4gFEbE1IrZFxJ119l8XERsjYkNE9EXElTX73hcRKyJiS0RsjogPj2YDkqQjGz/ShIgYB9wNfBwYANZGxGOZ+auaaU8Cj2VmRsQM4AfA+dV9fwP8Y2Z+JiJOAd47qh1Iko6okSv62cC2zHwlM98BHgSuq52QmW9mZlY3TwUSICL+ELgKuK86753MfGOUapckNaCRoJ8CvFazPVAdO0REXB8RW4AfATdXh/8UGAT+NiLWR8T3IuLUei8SEbdUl336BgcHm2pCkjS8RoI+6ozlYQOZP8zM84FPAV+rDo8HZgH3ZOZM4C3gsDX+6vH3ZmZvZvZ2d3c3UrskqQGNBP0AcGbNdg/w+nCTM/NZYHpEdFWPHcjMNdXdK6gEvyTpOGkk6NcC50TEtOrN1BuBx2onRMTZERHVr2cBpwC7M/Nfgdci4rzq1I8BtTdxJUljbMR33WTmvohYAqwExgH3Z+amiFhc3b8MuAH4fETsBX4LfK7m5uxfAd+v/iPxCvCfxqAPSdIw4vd53D56e3uzr6+v1WVIUseIiHWZ2Vtvn98ZK0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFa6hoI+IBRGxNSK2RcSddfZfFxEbI2JDRPRFxJVD9o+LiPUR8fhoFS5JasyIQR8R44C7gU8CFwI3RcSFQ6Y9CVycmZcANwPfG7L/S8DmY65WktS0Rq7oZwPbMvOVzHwHeBC4rnZCZr6ZmVndPBU48DUR0QP8ew4Pf0nScdBI0E8BXqvZHqiOHSIiro+ILcCPqFzVH/At4K+Bd4/0IhFxS3XZp29wcLCBsiRJjWgk6KPOWB42kPnDzDwf+BTwNYCIuBbYmZnrRnqRzLw3M3szs7e7u7uBsiRJjWgk6AeAM2u2e4DXh5ucmc8C0yOiC7gC+POI2E5lyWdeRCw/+nIlSc1qJOjXAudExLSIOAW4EXisdkJEnB0RUf16FnAKsDsz78rMnsycWj3uqcxcNKodSJKOaPxIEzJzX0QsAVYC44D7M3NTRCyu7l8G3AB8PiL2Ar8FPldzc1aS1ELRjnnc29ubfX19rS5DUgvs3buXgYEB3n777VaX0pYmTJhAT08PJ5988iHjEbEuM3vrHTPiFb0kHU8DAwOcfvrpTJ06leqKsKoyk927dzMwMMC0adMaPs5HIEhqK2+//TaTJk0y5OuICCZNmtT0/3YMekltx5Af3tH82hj0klQ4g15Sx1r2TD+r+3cdMra6fxfLnuk/pvOedtppx3R8uzHoJXWsGT0TWfLA+oNhv7p/F0seWM+Mnoktrqy9GPSSOtbc6V0sXTiTJQ+s55urtrLkgfUsXTiTudO7RuX8mcntt9/ORRddxIc+9CEeeughAN59912++MUv8sEPfpBrr72Wa665hhUrVhx2/Ec+8hHuuOMOZs+ezbnnnstPf/pTAPbv38/tt9/OZZddxowZM/jud7/b1Hmb5dsrJXW0udO7WDTnLL791DZunXf2qIU8wCOPPMKGDRt44YUX2LVrF5dddhlXXXUVzz//PNu3b+eXv/wlO3fu5IILLuDmm2+ue459+/bx85//nCeeeIKvfvWr/OQnP+G+++5j4sSJrF27lt/97ndcccUVzJ8/n3Xr1jV83mYY9JI62ur+XSxfs4Nb553N8jU7uHz6pFEL++eee46bbrqJcePG8YEPfICrr76atWvX8txzz/HZz36Wk046icmTJ/PRj3502HN8+tOfBuDSSy9l+/btAKxatYqNGzcevFrfs2cPL7/8clPnbUYRSzdjdUNGUns7sCa/dOFMbpt/3sFlnKF5cLSGe3JAM08UeM973gPAuHHj2Ldv38Hjv/Od77BhwwY2bNjAq6++yvz585s6bzOKCHpvyEgnpo0Dew5Zkz+wZr9xYM+onP+qq67ioYceYv/+/QwODvLss88ye/ZsrrzySh5++GHeffddfvOb3/D00083dd5PfOIT3HPPPezduxeAl156ibfeeuuYzzucIpZuam/ILJpzFsvX7BjVGzKS2tPiq6cfNjZ3eteo/d2//vrr+dnPfsbFF19MRPD1r3+dyZMnc8MNN/Dkk09y0UUXce655zJnzhwmTmz8wvILX/gC27dvZ9asWWQm3d3dPProo8d83uEU9VCzb67aevCGzG3zzxuDyiSNtc2bN3PBBRe0uowRvfnmm5x22mns3r2b2bNn8/zzzzN58uTjct56v0YnxEPNxvKGjCQNde211/LGG2/wzjvv8JWvfGVUQn6szltE0NfekJk7vYvLp08a9ffTSlKt0Vo/Px7nLeJm7FjfkJF0fLXjknK7OJpfmyKu6Mf6hoyk42fChAns3r3bRxXXceB59BMmTGjquCKCXlI5enp6GBgYYHBwsNWltKUDP2GqGQa9pLZy8sknN/XTkzSyItboJUnDM+glqXAGvSQVri2/MzYiBoF/HjLcBYzOk4paxx7agz20hxJ6gPbp408ys7vejrYM+noiom+4b+/tFPbQHuyhPZTQA3RGHy7dSFLhDHpJKlwnBf29rS5gFNhDe7CH9lBCD9ABfXTMGr0k6eh00hW9JOkoGPSSVLi2D/qIWBARWyNiW0Tc2ep6GhUR90fEzoh4sWbsjyLixxHxcvXz+1tZ40gi4syI+KeI2BwRmyLiS9XxjukjIiZExM8j4oVqD1+tjndMDwARMS4i1kfE49XtjqofICK2R8QvI2JDRPRVxzqqj4h4X0SsiIgt1b8XH+6EHto66CNiHHA38EngQuCmiLiwtVU17O+ABUPG7gSezMxzgCer2+1sH/BfM/MC4HLgL6u//p3Ux++AeZl5MXAJsCAiLqezegD4ErC5ZrvT6j/go5l5Sc37zjutj78B/jEzzwcupvJ70v49ZGbbfgAfBlbWbN8F3NXqupqofyrwYs32VuCM6tdnAFtbXWOT/fwD8PFO7QN4L/ALYE4n9QD0UAmQecDjnfpnCdgOdA0Z65g+gD8EXqX6JpZO6qGtr+iBKcBrNdsD1bFO9YHM/DVA9fO/a3E9DYuIqcBMYA0d1kd12WMDsBP4cWZ2Wg/fAv4aeLdmrJPqPyCBVRGxLiJuqY51Uh9/CgwCf1tdRvteRJxKB/TQ7kFf78fL+H7Q4ywiTgMeBr6cmf/W6nqalZn7M/MSKlfGsyPiohaX1LCIuBbYmZnrWl3LKLgiM2dRWYr9y4i4qtUFNWk8MAu4JzNnAm/Rjss0dbR70A8AZ9Zs9wCvt6iW0fCbiDgDoPp5Z4vrGVFEnEwl5L+fmY9UhzuuD4DMfAN4msq9k07p4QrgzyNiO/AgMC8iltM59R+Uma9XP+8EfgjMprP6GAAGqv8jBFhBJfjbvod2D/q1wDkRMS0iTgFuBB5rcU3H4jHgL6pf/wWVNe+2FZUf2HkfsDkzv1mzq2P6iIjuiHhf9es/AP4M2EKH9JCZd2VmT2ZOpfLn/6nMXESH1H9ARJwaEacf+BqYD7xIB/WRmf8KvBYR51WHPgb8ik7oodU3CRq4AXIN8BLQD/y3VtfTRN1/D/wa2EvlSuA/A5Oo3FR7ufr5j1pd5wg9XEllqWwjsKH6cU0n9QHMANZXe3gR+O/V8Y7poaaXj/D7m7EdVT+V9e0Xqh+bDvxd7sA+LgH6qn+eHgXe3wk9+AgESSpcuy/dSJKOkUEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQSw2qfnfn0upjjqWOYdBLjVsMvIfKdwxLHcOglxq3gMrjODa0uA6pKQa91ICImACMo/K0wmdaXI7UFINeasw5VIJ+S2bubXUxUjPGt7oAqUN0A+cC17W6EKlZXtFLjfljKj+A5aSIeH+ri5GaYdBLI4iIAz9CbjKwDNjf2oqk5vg8ekkqnFf0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQV7v8DzLiD9Tk9vWkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(Lsub_list,(LN)/np.log(2),'x',label='log neg')\n",
    "plt.xlabel('$\\ell$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEMCAYAAADUEk3/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAq2klEQVR4nO3de5yN5f7/8dfHmBxKYpAY4zAoyaDGYFQOlWSLrXYR2p0k1WR3PrerX7V1buecitpNNp2IUkoOYYQRKYc0I4fRJqMoymFmrt8fa2a+c541rJl1mPfz8fB4WOu+1z0fWb27uu7r/lzmnENERIJfFX8XICIivqFAFxEJEQp0EZEQoUAXEQkRCnQRkRChQBcRCRFVvTnJzPoALwNhwGvOuacLHK8DTAGigUPA9c6570q6Zr169VyzZs2OpWYRkUpr9erV6c65+kUdKzXQzSwMGA9cBKQBq8xstnNuQ57THgTWOucGmtkZ2edfUNJ1mzVrRnJysrd/BhERAcxsW3HHvJlyiQNSnHNbnHNHgOnAgALnnAl8AeCc2wQ0M7NTj7FeERE5Bt4EemNgR57Xadnv5fUNcBmAmcUBTYFIXxQoIiLe8SbQrYj3CvYLeBqoY2ZrgduANUBGoQuZjTCzZDNL3rNnT1lrFRGREnhzUzQNaJLndSTwU94TnHO/AdcBmJkBP2b/osB5k4HJALGxsYWayBw9epS0tDQOHTrkbf1yDKpXr05kZCTh4eH+LkVEfMibQF8FtDKz5sBOYDAwJO8JZnYK8Ef2HPtw4MvskC+TtLQ0atWqRbNmzfD8d0F8zTnH3r17SUtLo3nz5v4uR0R8qNQpF+dcBpAAzAM2Au8459ab2UgzG5l9WhtgvZltAi4B/nEsxRw6dIiIiAiFeTkyMyIiIvR/QSIhyKt16M65ucDcAu9NyvP75UArXxSkMC9/+mcs4kcpX0C91nBKk9LPLSM9KVrASSed5O8SRCQU/fELzLoFEi+DpS+Wy48I2kCftDiVpNT0fO8lpaYzaXGqnyoSESnGhtkwvjN8Mx3OuwsuHl0uPyZoAz0msjYJ09bkhnpSajoJ09YQE1nbJ9d3znHPPfdw1lln0a5dO2bMmAFAVlYWt9xyC23btqVfv3707duX9957r9Dne/TowX333UdcXBytW7dmyZIlAGRmZnLPPffQqVMnYmJieOWVV8p0XREJIr/vhhlXwztXQ62GMGIRXPBPCK9eLj/Oqzn0QBQfXY9xQzqSMG0NwzpHkbhiO+OGdCQ+up5Prv/BBx+wdu1avvnmG9LT0+nUqRPnn38+y5YtY+vWrXz77bf8/PPPtGnThuuvv77Ia2RkZLBy5Urmzp3L448/zvz583n99depXbs2q1at4vDhw3Tr1o3evXuzevVqr68rIgHOOc9o/NP74eifcMGjEH8bhJXvUuGgDXTwhPqwzlGMWZDCqF4tfRbmAEuXLuWqq64iLCyMU089le7du7Nq1SqWLl3KFVdcQZUqVWjYsCE9e/Ys9hqXXXYZAOeccw5bt24F4LPPPmPdunW5o+/9+/fzww8/lOm6IhLA9m2HObdD6hfQpAv0Hwv1W1fIjw7qQE9KTSdxxXZG9WpJ4ortdImO8FmoF7d5dlk21a5WrRoAYWFhZGRk5H5+7NixXHzxxfnO/fjjj4+xUhEJCFlZkPw6zH/MM0K/5DnoNByqVNzMdtDOoefMmY8b0pE7e5+eO/1S8EbpsTr//POZMWMGmZmZ7Nmzhy+//JK4uDjOPfdc3n//fbKysti9ezeLFi0q03UvvvhiJk6cyNGjRwHYvHkzBw8ePO7riogfpf8Ab/SFuXdDk85w61fQeUSFhjkE8Qh9Xdr+fHPmOXPq69L2+2SUPnDgQJYvX0779u0xM5599lkaNmzI5ZdfzhdffMFZZ51F69at6dy5M7Vre38jdvjw4WzdupWzzz4b5xz169dn1qxZx31dEfGDzKOQNBYWPQ3hNeCvE6H9VeCnZz2sLFMIvhQbG+sK9kPfuHEjbdq08Us9ZXHgwAFOOukk9u7dS1xcHMuWLaNhw4YBe92iBMs/a5GA9b918OGtsGsdtOkPfZ+HWuXfNdzMVjvnYos6FrQjdH/q168f+/bt48iRIzzyyCM+C93yuq6I+NDRQ/Dls7D031AzAq78D5xZcIsI/1CgH4Pymt/WvLlIgNu+AmYnQPpmaD8ELn4Katb1d1W5FOgiIqU5fAAWPAErXoHakTDsfWh5ob+rKkSBLiJSktQFMOcfnvXlcSM8T3pWq+XvqoqkQBcRKcqfv8K8h2FtIkS0hOs+haZd/V1ViRToIiIFbZwDH98FB9Ph3Duh+33l1n/Fl4L2waLysmvXLgYPHkx0dDRnnnkmffv2ZfPmzWzdupUaNWrQoUOH3F+PP/547u/DwsJyfz9mzBif17Vp0ya6du1KtWrVeP755wsdv+mmm1i2bBn33HMPZ5xxBjExMQwcOJB9+/b5vBaRkHXgZ3jnGpgxDE5qADcugAsfDYowB43Q83HOMXDgQK655hqmT58OwNq1a9m9ezdNmjQhOjqatWvX5vvMo48+Cnj6qBc85kt169ZlzJgxzJo1q8jjK1asYMKECRw8eJDRo0dTtWpV7rvvPkaPHs0zzzxTbnWJhIRCzbT+CfGjyr2Zlq9phJ7HwoULCQ8PZ+TIkbnvdejQgfPOO89nPyMzM5O7776bdu3aERMTw9ixY736XIMGDejUqVORGztv3LiR1q1bExYWRu/evala1fPf6S5dupCWluaz2kVC0r4d8PbfYNZIqH86jFzq6VkeZGEOgTxC/+R+2PWtb6/ZsB1c8nSxh7/77jvOOeecYo+npqbSoUMHALp168b48ePLXMLkyZP58ccfWbNmDVWrVuWXX34B4I477mDhwoWFzh88eDD3339/idf85JNP6NOnT6H3p0yZwqBBg8pco0ilEADNtHwtcAM9ABU15VJW8+fPZ+TIkbmj6Lp1PQ8lvPTSS8d8zXnz5jF16tR87z311FNUrVqVoUOHHnuxIqEqPcXzgND25dCiJ1z6MtRp6u+qjlvgBnoJI+ny0rZt23LfJcg5V+Qmzcc6Qv/jjz/Yt28fjRo1yn3vzTff5KOPPuKLL77QhtAieWVmwPKxsHC050bngAnQYYjfmmn5mleBbmZ9gJeBMOA159zTBY7XBhKBqOxrPu+cm1roQgGuV69ePPjgg7z66qvceOONAKxatYo//viDpk3L9l/vmTNnsnLlSkaPzr93YO/evZk0aRI9evTInXKpW7fuMY/QFy5cmG8zjE8//ZRnnnmGxYsXU7NmzWO6pkhI2vWtp5nW/76BNpdmN9MKrX5JpU4WmVkYMB64BDgTuMrMzixw2q3ABudce6AH8IKZneDjWsudmTFz5kw+//xzoqOjadu2LY899li+0a+3UlNTOfnkkwu9P3z4cKKiooiJiaF9+/ZMmzbNq+vt2rWLyMhIXnzxRZ588kkiIyP57bffCs2fJyQk8Pvvv3PRRRfRoUOHfDd4RSqljMPwxRMwuQf89j9PM61BiSEX5uBF+1wz6wo85py7OPv1AwDOudF5znkAaIIn2JsBnwOtnXNZxV03mNvnemPYsGG89NJL1K9fv1x/ztlnn82KFSuKXP1SklD6Zy1SrB0r4cMESP8+IJtpHYvjbZ/bGNiR53Ua0LnAOeOA2cBPQC1gUFFhbmYjgBEAUVFRXvzo4JWYmFghP+frr7+ukJ8jElSOHPSMyldM8jTTGvo+tAq8Zlq+5k2gF3W3oOCw/mJgLdALiAY+N7Mlzrnf8n3IucnAZPCM0MtcrYhIaVIXwpxRnmZanYbDhY8FbDMtX/Mm0NPwTKfkiMQzEs/rOuBp55m/STGzH4EzgJU+qVJEpDR/7oPPHoI1iVA3Gq77BJrG+7uqCuVNoK8CWplZc2AnMBgYUuCc7cAFwBIzOxU4HdhyLAUVt6xPfMdf2w6KlJtNH8NHd8LBPXDuHdnNtGr4u6oKV2qgO+cyzCwBmIdn2eIU59x6MxuZfXwS8ATwhpl9i2eK5j7nXHpZi6levTp79+4lIiJCoV5OnHPs3buX6tWDo9mQSIkO7IFP7oH1M+HUdjBkOjTq6O+q/CagNok+evQoaWlpHDp0yC81VRbVq1cnMjKyzCtjRAKGc7DuHfj0Ps8N0PPvhXNvD8r+K2UVNJtEh4eH07x5c3+XISKBbH8azLkdUj6HyDgYMM7TVEsCK9BFRIqVlQWrp8Dnj4HLhD7PQNyNUCXM35UFDAW6iAS+vakw+zbYtgxa9MhuptXM31UFHAW6iASuzAxYPg4WjYaq1aD/OOg4LGSaafmaAl1EAtOubz2P7f9vLZzRD/7yQkj2X/ElBbqIBJaMw/Dlc7D0JahRB654E84coFG5FxToIhI4dqzybDyxZxO0vwou/lfQN9OqSAp0EfG/IwdhwZPw1UQ4uTEMfQ9aXeTvqoKOAl1E/GvLIpg9CvZtq3TNtHxNgS4i/vHnPvj8Efj6P55mWtfOhWbd/F1VUFOgi0jF2zQXPr4TDvwM3W6HHvdXymZavqZAF5GKc2APfHIvrP8ATj0LBk+Dxmf7u6qQoUAXkfLnHHz7LnxyHxw5AD0frjTNtCqSAl1Eytf+NE+v8h/mQWQnz9OeDc7wd1UhSYEuIuUjKwtWT4XPH81upvU0xI1QM61ypEAXEd/bm+pZirhtKTQ/Hy4dA3XVGru8KdBFxHcyM+Cr8bDwXxBWDfqPhY5X67H9CqJAFxHf2PWd57H9n9bA6X/xNNM6+TR/V1WpKNBF5PhkHIYvn4elL0L1U+BvU6HtQI3K/UCBLiLHLm8zrXZXem58nhjh76oqLQW6iJTdkYOw4Cn4agKc3AiGvAute/u7qkrPq0A3sz7Ay0AY8Jpz7ukCx+8Bhua5ZhugvnPuFx/WKiKBYMtimDMKft0KsTd4mmlVP9nfVQleBLqZhQHjgYuANGCVmc12zm3IOcc59xzwXPb5lwJ3KMxFQsyf++Dzf8LXb0LdFmqmFYC8GaHHASnOuS0AZjYdGABsKOb8q4D/+qY8EQkIuc20dkP8KOj5oJppBaAqXpzTGNiR53Va9nuFmFlNoA/wfjHHR5hZspkl79mzp6y1ikhFO5gO710P069ib9aJfHPJ+9D7idwwT0pNZ9LiVD8XKTm8CfSi1h65Ys69FFhW3HSLc26ycy7WORdbv359b2sUkYrmHKx7F8Z1gg2zoedDbB4wh+vmZZKUmg54wjxh2hpiImv7uVjJ4c2USxrQJM/rSOCnYs4djKZbRILb/p2e6ZXNn/J7vQ7Uum4SNGhDV+DmHn9ywxvJ3HhecxJXbGfckI7ER9fzd8WSzZsR+iqglZk1N7MT8IT27IInmVltoDvwoW9LFJEKkZUFyVNgfGfYspgfz3mInr88QNLvnv+bTkpNZ+KiLfQ5qyFjFqQwrHOUwjzAlDpCd85lmFkCMA/PssUpzrn1ZjYy+/ik7FMHAp855w6WW7UiUj72prLzrRE03pec20yred3m3HRKar4R+c09WjBx0RZG9WpJ4ortdImOUKgHEHOuuOnw8hUbG+uSk5P98rNFJFtWpufhoAVPkWFh/CtjGBcOvZv4lvVz58i7t67PzDU7GdixEYs3p+dOs+Qc17RLxTKz1c652KKO6UlRkcpq9wZ+ThxOg9/XQ+tLqNrvRS7ccwI3vbWado1rs2nX7/lG5K8u+ZE7e7fKDe/46HqMG9KRdWn7FegBQoEuUtlkHIElL8CSF6gTXov77Xb6xyUQf3J92JPO0cwsklL3MrBjIyYu2pI7Au8SHUHCtDW0bVQ7X6grzAOHAl2kEvlgzof0SXmCmvt/gHZXEt7naaK/3scNb67mxvOaMzVpK+FhVRhxXguNyIOQAl2kMjjyByx8ioFfT2C3q8O2nq/RpvsVhVauVA+vwpRrO2lEHqS8WbYoIkFq0uJUvls2ByZ2heXjsHOu5dPuM7lsfi1e/Ox7Eqat4eYeLVi8eQ/doiMID/u/SMg7IpfgoBG6SKg6tJ+/7niGhikz+LNWU2pc+zFJmWcwZtqa3BF5wXnygitXNCIPLgp0kRAzaXEq3V0ybVY/SsMDu9l55o38dcP5tPq8Kpt2rdHKlRCmQBcJcpMWpxITmT3PfTCdy398lPpb5/C/6i04bfh8th1qyu/frtLKlUpAc+giQS4msjYJb3/N5vlTYHwcEds/ZTxX0vfQk7y4/iRuems14WFVGNWrJZ9+t5ube7QockQuwU8jdJEgk29EDsTXP8LMumNounQJ/zupLaN4iDuuHsDVqXu1cqWS0QhdJMjERNYmYdoaklL2wOo3yBjbiQbpK5jV4Fa6pT9A1y7nApC4YrtWrlQy6uUiEgQKjspXr11Nxqzb6Mx6VnEWqV2f4tmVRxnWOYqpSVsBeOXqc9RzJQSV1MtFI3SRIJA7Kv9hNySNo8OcvpzptnD/0eFMO2Msz648yrghHbmz9+n0izkt32c1Kq88NIcuEoAKzZNH1+PB2CxOTOwLlsISzuZJu5G+3c8ptPRw9GUxXNq+Ub6lh5onrxwU6CIBKGdEPm5IR+Kbnsz22U8yYN14Dlc9kVF/JvBZWDemXBunG52Sj6ZcRALApMWpuXt1gieQb+7RgpfemE76i12IWvcy2069iP7uRfY2v5TwsLB852pKRUCBLhIQcufIs0P9q+93UG3BY8wIe4SjB3/ltch/cWX6DTw5tAdv39iFV64+J9/58dH1GNk92p9/BAkAmnIR8ZO88+Q5o+yb3lrNlfW2ck36i3RhFx/YhfwU9yDjl+/hzt5FPxCkqRXJoRG6iJ8UHJWHHfmdh7Im80j6vdQIr8KIKo/RcNgrJPQ9h9evjWXioi2FpmU0Kpe8FOgiFaS4efIb3khm5vTXaTajF1dU+YLkRkPpfehpOvXsr0f0pUwU6CIVpOCIPCk1nf8uXENinVcZuOlO9rsTWX/J+8SOmMD4a7tpRC5l5lWgm1kfM/vezFLM7P5izulhZmvNbL2ZLfZtmSLBp6QR+YvzNjE7cSwfh91FzG+LeL/WMAbZMxyo1z73XI3IpaxKvSlqZmHAeOAiIA1YZWaznXMb8pxzCjAB6OOc225mDcqpXpGgkW8tefYj+BMXbWHQ6WG0W3oLF4Wt5tujLcnsN43LY7txmjaXkOPkzQg9Dkhxzm1xzh0BpgMDCpwzBPjAObcdwDn3s2/LFAl8JY7IP/uehLe/5t+t1nJXyt/pEf4dz2QNY0Wv/9Ihtlvu+RqVy/HwJtAbAzvyvE7Lfi+v1kAdM1tkZqvN7O9FXcjMRphZspkl79mz59gqFglQRc2R52zAPHPhMqZXH835m54k69R2hN+6nPOueZwJi7dpnlx8xpt16FbEewVbNFYFzgEuAGoAy83sK+fc5nwfcm4yMBk83RbLXq5I4MoZYSdMW8OwzlEkrtjOLd2bsn/hOJ6tMZ3DB6uwpM1DnHfl3VClCvERaC25+JQ3I/Q0oEme15HAT0Wc86lz7qBzLh34EmjvmxJFAlPBKZYcZzSsxZgFKQxq+judFwzhLvcm4S17sunyz/nHDx1J+vGX3HM1Ihdf8ibQVwGtzKy5mZ0ADAZmFzjnQ+A8M6tqZjWBzsBG35YqEliKmmK56a3VbNy5l2mtFnFH6nBahu+By1+Hq6YTG9NOc+RSrkqdcnHOZZhZAjAPCAOmOOfWm9nI7OOTnHMbzexTYB2QBbzmnPuuPAsXqUgF29nmuLjtqblTLFOTtnImqUw55Q1O3PE9e5r3Y9D2gTxZswfx5pm51MoVKU9e9XJxzs0F5hZ4b1KB188Bz/muNJHAUdQSxJzX9U+qxuQF63m+3lz+cvB9LPNUGPxf6p/RlydT0zVHLhVGzblEvFDUDc9xQzoCsOGrT1h+ymvUObCDXS0H0fBvz0H12rmfU5hLRdGj/yJFKO2G57DOUYQd+Z2diTfzWtaj1KkRxncXvkXfH68gaedRP1QsokAXKVJxNzy/3bmfUb1a8uNXszj9g978jfnQNQFuTuKsc/vrpqf4laZcRIpQcIplatJWAF6/ojlx3z8PWTNIzYpke593ielyYb7PaYpF/EUjdKn0ippeScq+mTmscxRjFqQQ0+hk3jtvF3EfXwLfvQ/d7+fnIZ+RdLi5n6oWKUyBLpVeUdMrCdPWEFYFElds54FzazP8p0c4fckoOKUJ3PQl9HyArq0b66EgCSiacpFKr6gVLDf3aMHEham8E5dCyzVPkxl2iBez/k6Xng8Rf2pDf5csUiSN0KVS8WZ6ZVjnKGoeTGNBgxdp+dUD0LAdYbcsp8uwR1n300E/VS5SOgW6VCqlTa/8o2dzspZPYPCqQZzyy3fQ7yW4Zg5ERKvvigQ8TblIpVLs9MqiLbze9yQ6rhkFbhVLMs+mxoAxxMa083fJIl7TCF1CUnFTK5MWpxIfXS/f9IrLOMqss5bRcW5/2JsKl71G2LB3SP61pp+qFzk2CnQJScVNrcRE1iYpNZ3EFdsZ1aslq79ayNB11xD1zUtwRj+4dSXEXEF8y/qaXpGgoykXCUkl9V5JmLaGCVe2ocv2yTg3lvRfa7O9xyu06THYz1WLHB+N0CXoebtyJT66HuvS9vOfC47SZV5/WPYy1nEYWwYtYLF18lP1Ir6jQJegV9rKlVG9WpK4YjsrNv7IyN/Hc9ZnV0FWBvz9Q+g/ls5tWmh6RUKCplwk6JW0ciWnf3mf6t9SZ8aFOPZiXW6BXg/DCSf6u3QRn9IIXUJCwZUrmVmeDZjjT6sCH9zEmQtuoHbtOsw8eyr0Ga0wl5CkQJegUNIyxJzf551eiWl8MvGHlsD4OPjuPTj/XmretozL+g/0R/kiFUJTLhIUStoCLu/v46Prcd5pmfyeOARYCad1gKtnQcOz/PwnECl/5pzzyw+OjY11ycnJfvnZEpxygjvvMsT46Hr/t4FziwhYkwjzHiIz4xArm42k65B/QpjGLRI6zGy1cy62qGP6pkvQyDtPPqpXy9yNJEZ2j4Zft8JbN8CWRdC0G2H9x9I1QitXpHLxag7dzPqY2fdmlmJm9xdxvIeZ7Teztdm//un7UqWyKzhPnpSaDlmZ8NVEmNAV0lbDX16Eaz4ChblUQqWO0M0sDBgPXASkAavMbLZzbkOBU5c45/qVQ41SSeROneTZwi3nAaGCc+hdoiN44e05nBXxJifvWQOtens6I9aO9OOfQMS/vBmhxwEpzrktzrkjwHRgQPmWJZVRSf1X1qXtzw1zMo8Sv3Mq73Iv4ft+hMtehSHvKMyl0vNmDr0xsCPP6zSgcxHndTWzb4CfgLudc+sLnmBmI4ARAFFRUWWvVkJacf1X8m28/NMa+PA22P0tVdoOpMYlz8FJ9f1buEiA8CbQrYj3Ci6N+Rpo6pw7YGZ9gVlAq0Ifcm4yMBk8q1zKVqpUBsXd+OTon7DoaUgaCyfWh0FvQxvN8Ink5U2gpwFN8ryOxDMKz+Wc+y3P7+ea2QQzq+ecy/8kiEgpCt747BIdQXzVH2B2AuxNgY5XQ+8nocYp/i5VJOB4E+irgFZm1hzYCQwGhuQ9wcwaArudc87M4vDMze/1dbES3Eq66Tmye3ShB4Tio6qxLfFWcJ/CKU09zbRa9PDfH0AkwJV6U9Q5lwEkAPOAjcA7zrn1ZjbSzEZmn/Y34LvsOfQxwGDnryeWJGCVdNMTyH/j84f5dPmkH1e6eayLHAK3LFeYi5RCT4pKhSruac9cf/wC8x6Eb/4L9U6HAeOgSZz/ChYJMHpSVAJGsTc9AdbPgrl3w5+/wnl3Q/d7oWo1v9UqEmwU6FKhirzp2SDDE+Qb58Bp7eHqmdCwnb9LFQk6CnTxidJueOa8zve0Z4u6fJb4AnHhiVTNPAwXPgZdb1MzLZFjpH7o4hOl3fCEAjc9f91GfNKNPOYm8HONaLg5Cc69Q2Euchx0U1R8ptQbngBZWbDqVZj/OJjBRY/DOddDFY0tRLyhm6JSIUq84Qmw53uYfRvsWAEtL/I00zqlSdEXE5EyU6CLzxR5wzO7mRbLXobFz3j28hz4CsQM8ozQRcRnFOjiE4VueEZHkDBtDVMvDqf91w/Drm+h7UC45Fk4qYG/yxUJSQp08Yl8NzyB+KiTmH3GfBp98iqcWE/NtEQqgAJdiuXNUsQc+V5vWw6zE4jcmwIdh2U306pTUWWLVFpaWiDF8mYpYj6Hf4eP74apfSDziOcBoQHjFeYiFUQjdClWSRtOFJIyH+bcDvvToPPN0OthqHZShdcsUpkp0KVEpS5F/OMXmPcQfDMN6rWG6+dBVFEbWolIeVOgS4mKXYoIsGE2fHwX/PmLp5nW+fdAeHX/FixSiSnQpVjFLUV8ZWAkndb/CzbO9jTTGvY+nBbj73JFKj0FuhSr0FLEFhHM6JxKk1nXg1MzLZFAo38TpVj5liL+ug0+up1WqQsgqiv0Hwv1Cu0DLiJ+pGWLlcSkxam5yw9zJKWmM2lxaskfzMqCFa/AhK6wYyX0fR6unaswFwlACvRKosxrygH2bIapl8An90JUF8++nnE3qjOiSIDSlEslUaY15XmbaYXXhL9OgvaD1UxLJMAp0CuRUteUA/zvG/jwVk8zrTMHeKZY1ExLJCh49f/OZtbHzL43sxQzu7+E8zqZWaaZ/c13JYqvFFxTnm9O/eghz6YTk3vCgZ/hyrfgyv8ozEWCSKkjdDMLA8YDFwFpwCozm+2c21DEec8A88qjUDk+xa0pHzekI/HhKfBhAuz9AToMg4vVTEskGHkz5RIHpDjntgCY2XRgALChwHm3Ae8DnXxaofhEoTXl0fWYcEVrTlzwEOycAbWbwLAPoOUFfq5URI6VN4HeGNiR53UakK9Zh5k1BgYCvSgh0M1sBDACICoqqqy1ynEo2O6WlC/o8sntsH8HxI2AC/6pZloiQc6bQC9qaUPBnaX/DdznnMu0ElZCOOcmA5PBs0m0lzWKL/35q6eZ1tq3IaIVXP+pZ0miiAQ9bwI9Dci7k28k8FOBc2KB6dlhXg/oa2YZzrlZvihSfGTDbJh7NxxMh/PugvPvVTMtkRDiTaCvAlqZWXNgJzAYGJL3BOdc85zfm9kbwEcKc98qy+5Bhfy+2xPkG2dDw3Yw9F1PUy0RCSmlLlt0zmUACXhWr2wE3nHOrTezkWY2srwLFI9jetLTOVj7XxgfB5vnQa9H4MaFCnOREGXO+WcqOzY21iUnJ/vlZwernBAv9UlPgH3bPTsIpX4BTTpD/3FQv3WF1isivmdmq51zsUUd05OiQcSrJz2zsiD5dZj/mGeEfsmz0En9V0QqAwV6EClx9yCA9B9g9m2wfTlE94J+/4Y6Tf1Wr4hULAV6kCjxSc9mp0DSGFj0NITXgL9OhPZXqZmWSCWjQA8SRT3pOW5IR3ZuXAnzn/Y01WpzKfR9AWqd6udqRcQfFOhBotDSxKOHiN86AVb/G2pGwBVvQtu/+qM0EQkQCvRgtH0FzE6A9M3QYSj0fhJq1vV3VSLiZwr0YHL4ACx4wrMlXO0mMOx9aHmhv6sSkQChQA8WqQtgzj9gn5ppiUjRFOiB7s9fYd7DsDbR00zruk+gaVd/VyUiAUiBHsg2zoGP7/I00zr3Tuh+n5ppiUixFOiB6MDPMPce2DDL00xryDvQqIO/qxKRAKdADyTOwboZ8On9cOQPzzx5/CgIC/d3ZSISBNTgw0cmLU7Nv+kynqc7Jy1O9e4C+3bA21fAzJugXmsYudTTs1xhLiJeUqD7yDG1twVPM62Vr8KELrAtCfo847nxqc6IIlJGmnLxkZxH8b1ubwuQnpLdTCsJWvSES19WMy0ROWYKdB/yqr0tQGYGLB8LC0d7Vq0MmAAdhqiZlogcFwW6D5Xa3hZg17fw4a2eZlpn9IO/vAC1GvqnYBEJKQp0HymxvW10Pcg4DIufhWX/hhp14Io34My/alQuIj6jQPeR4trbrkvbT/wJW+DDBEj/3tOn/OJ/qZmWiPicAt1HCrW3BeIjqxO/+Xl4fRLUjoSh70MrNdMSkfLh1bJFM+tjZt+bWYqZ3V/E8QFmts7M1ppZspmd6/tSg0zqApjYFVZMhE7D4ZblCnMRKVeljtDNLAwYD1wEpAGrzGy2c25DntO+AGY755yZxQDvAGeUR8EB78998NlDsCYRIlpmN9OK93dVIlIJeDPlEgekOOe2AJjZdGAAkBvozrkDec4/EXC+LDJobPwou5nWHjj3juxmWjX8XZWIVBLeBHpjYEee12lA54InmdlAYDTQAPiLT6oLFnmbaZ3aDoZMh0Yd/V2ViFQy3gR6UevqCo3AnXMzgZlmdj7wBFBowtjMRgAjAKKiospWaSDK10zrIPR6GLrdrv4rIuIX3gR6GtAkz+tI4KfiTnbOfWlm0WZWzzmXXuDYZGAyQGxsbHBPy+zbAR/dASmfQ2QcDBgH9U/3d1UiUol5E+irgFZm1hzYCQwGhuQ9wcxaAqnZN0XPBk4A9vq62ICQlQWrp8Dnj4LLgj5Pe7aEqxLm78pEpJIrNdCdcxlmlgDMA8KAKc659WY2Mvv4JOBy4O9mdhT4ExjknAvuEXhR9qZ6mmltWwYtemQ302rm76pERAAwf+VubGysS05O9svPLrPMDFg+DhaNhqrVPE96dhiqx/ZFpMKZ2WrnXGxRx/SkaGl2fQezE+CnNWqmJSIBTYFenIzD8OXzsPRFNdMSkaCgQC/KjlWeUfmeTWqmJSJBQ4Ge15GDsOBJ+GoinNwYhr4HrS7yd1UiIl5RoOfYsghmj4J926DTjXDho1Ctlr+rEhHxmgL9z33w+SPw9X+gbjRcOxeadfN3VSIiZVa5A33Tx/DRnXDwZ+j2D+jxgJppiUjQ8qofeiCYtDiVpNR8nQRISk1n0uLUsl/swB549zqYPgRqRsDwL+Ci/6cwF5GgFjSBHhNZm4Rpa3JDPWcPz5jI2t5fxDlY9w6Mj4NNH0HPh2HEImh8dvkULSJSgYJmyiVnj86EaWsY1jmKxBXb8+3hWar9aZ7plR/mQWQn6D8OGlTOPThEJDQFTaCDJ9SHdY5izIIURvVq6V2YZ2XB6qnZzbQy1UxLREJWUAV6Umo6iSu2M6pXSxJXbKdLdETJob431bMUcdtSaN7d00yrbvOKK1hEpAIFTaDnzJnnTLN0iY7I9zqfzAz4agIsfArCqnmmVzoO02P7IhLSgibQ16XtzxfeOXPq69L25w/0vM20Tv+Lp5nWyaf5qWoRkYoTNIE+snt0offio+v9X5hnHIYlL3h+VT8F/jYV2g7UqFxEKo2gCfQSpSXDh7d6mmnFDPLc+FQzLRGpZII70I8chAVPeebLT24EQ96F1r39XZWIiF8Eb6BvWQxzRsGvWyH2BrjwMah+sr+rEhHxm+AL9MMHYN4D2c20WqiZlohItuAL9CpVPRtQqJmWiEg+wRfo4dXhpsWezZpFRCSXV825zKyPmX1vZilmdn8Rx4ea2brsX0lm1t73peahMBcRKaTUQDezMGA8cAlwJnCVmZ1Z4LQfge7OuRjgCWCyrwsVEZGSeTNCjwNSnHNbnHNHgOnAgLwnOOeSnHO/Zr/8Coj0bZkiIlIabwK9MbAjz+u07PeKcwPwyfEUJSIiZefNTdGinp13RZ5o1hNPoJ9bzPERwAiAqKgoL0sUERFveDNCTwOa5HkdCfxU8CQziwFeAwY45/YWdSHn3GTnXKxzLrZ+/frHUq+IiBTDm0BfBbQys+ZmdgIwGJid9wQziwI+AK52zm32fZkiIlKaUqdcnHMZZpYAzAPCgCnOufVmNjL7+CTgn0AEMME83Q0znHOx5Ve2iIgUZM4VOR1e/j/YbA+wrYwfqw3sP44feyyfL8tnSjv3eI4Xd6wekO5Vdf5xvH9n5Xltf38fSjtH34eKv3ZZr+GP70NT51zRc9bOuaD5BUyu6M+X5TOlnXs8x4s7BiT7+++lPP/OQvn7cKx/5/o+lN+1y3qNQPk+5Pzy6knRADLHD58vy2dKO/d4jh/vn91fyrPuYP8+lHaOvg8Vf+2yXiNQvg+AH6dcxDfMLNnpfoVk0/ehcgu2EboUpjYLkpe+D5WYRugiIiFCI3QRkRChQBcRCREKdBGREKFAD2Fm1sLMXjez9/xdi/iHmZ1oZm+a2atmNtTf9Uj5UqAHKDObYmY/m9l3Bd4vcfeovJynh/0N5VupVLQyfjcuA95zzt0I9K/wYqVCKdAD1xtAn7xvFLd7lJm1M7OPCvxqUPElSwV5Ay+/G3i6o+bsZ5BZgTWKHwTfJtGVhHPuSzNrVuDt3N2jAMxsOp52xaOBfhVcovhJWb4beNpfRwJr0QAu5OkvOLiUafcoM4sws0lARzN7oLyLE78q7rvxAXC5mU0keNsFiJc0Qg8uXu8eBeA8G42MLL9yJIAU+d1wzh0ErqvoYsQ/NEIPLl7tHiWVkr4bokAPMqXuHiWVlr4bokAPVGb2X2A5cLqZpZnZDc65DCBn96iNwDvOufX+rFMqnr4bUhw15xIRCREaoYuIhAgFuohIiFCgi4iECAW6iEiIUKCLiIQIBbqISIhQoIuIhAgFuohIiFCgixSQvcvPODPr4u9aRMpCgS5S2EigGnCuvwsRKQsFukhhfYDNeDaFEAkaCnSRPMysOhAGnA0s9nM5ImWiQBfJrxWeQN/knDvq72JEykI7FonkVx9ojWc/TpGgohG6SH6NgPeBKmZWx9/FiJSFAl0km5lVxTN33hCYBGT6tyKRstEGFyIiIUIjdBGREKFAFxEJEQp0EZEQoUAXEQkRCnQRkRChQBcRCREKdBGREKFAFxEJEf8fptiQu7UqIAwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "# plt.plot(Lsub_list,(LN)/np.log(2),'x',label='log neg')\n",
    "plt.plot(np.tan(pi*Lsub_list/N),(LN),'x',label='log neg')\n",
    "xx = np.logspace(-1.5,0.8)\n",
    "plt.plot(xx,np.log(xx)/8+0.65, label=\"CFT, c=1/2\")\n",
    "plt.xlabel('$\\ell$')\n",
    "plt.xscale(\"log\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEMCAYAAADUEk3/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAArV0lEQVR4nO3deXhU1f3H8feXsIkIyi6EAAYRQZElLAYQpLKIIMWlIqLWqogacakLalu1anFpXRA04l4pRX+uCCgqq4IiIIgsggQQAiIEBBRkSXJ+f9wknSSTZBJmMks+r+fJ8zAzNzffMeOHw7nnfo855xARkehXKdwFiIhIcCjQRURihAJdRCRGKNBFRGKEAl1EJEYo0EVEYkTlQA4yswHA00Ac8KJz7pECr58AvAwkAgeBPznnVhZ3znr16rnmzZuXpWYRkQpr6dKlGc65+v5eKzHQzSwOmAD0BdKBxWY21Tm32uewe4DlzrmhZtY65/jfFXfe5s2bs2TJkkDfg4iIAGb2Q1GvBTLl0gVY75zb4Jw7DEwBhhQ4pg0wC8A59x3Q3MwalrFeEREpg0ACvQmwxedxes5zvr4BLgAwsy5AMyA+GAWKiEhgAgl08/NcwX4BjwAnmNly4CZgGZBZ6ERmI81siZkt2blzZ2lrFRGRYgRyUTQdaOrzOB7Y5nuAc24fcBWAmRmwMeeLAsdNBCYCJCUlFWoic+TIEdLT0zl48GCg9UsZVK9enfj4eKpUqRLuUkQkiAIJ9MXAyWbWAtgKDAOG+x5gZscDB3Lm2K8B5ueEfKmkp6dz3HHH0bx5c7y/FyTYnHPs2rWL9PR0WrRoEe5yRCSISpxycc5lAinATGAN8KZzbpWZjTKzUTmHnQqsMrPvgHOBm8tSzMGDB6lbt67CPITMjLp16+pfQSIxKKB16M65GcCMAs+l+vz5C+DkYBSkMA89/TcWCaP1s6BeKzi+acnHlpLuFC2gZs2a4S5BRGLRgd3w3g0w6QL4/ImQ/IioDfTUeWksTMvI99zCtAxS56WFqSIRkSKsngoTusI3U6Dnn6H/2JD8mKgN9HbxtUmZvCwv1BemZZAyeRnt4msH5fzOOe644w5OO+00Tj/9dN544w0AsrOzueGGG2jbti2DBg1i4MCBvPXWW4W+v3fv3tx111106dKFVq1a8dlnnwGQlZXFHXfcQefOnWnXrh3PP/98qc4rIlHkl5/gjcvhzcvhuIYwcg787m9QpXpIflxAc+iRKDmxHuOHdyBl8jJGdE1g0qLNjB/egeTEekE5/zvvvMPy5cv55ptvyMjIoHPnzpx11lksWLCATZs28e2337Jjxw5OPfVU/vSnP/k9R2ZmJl999RUzZszggQce4NNPP+Wll16idu3aLF68mEOHDtG9e3f69evH0qVLAz6viEQ457zR+Edj4MhvXognj4a40C4VjtpABy/UR3RNYNzs9Yzu0zJoYQ7w+eefc+mllxIXF0fDhg3p1asXixcv5vPPP+fiiy+mUqVKNGrUiLPPPrvIc1xwwQUAdOrUiU2bNgHw8ccfs2LFirzR9969e/n+++9LdV4RiWB7NsMHt0DaLGjaFc4fD/VblcuPjupAX5iWwaRFmxndpyWTFm2mW2LdoIV6UZtnl2ZT7WrVqgEQFxdHZmZm3vc/88wz9O/fP9+x06dPL2OlIhIRsrNhyUvw6f3eCP3cx6DztVCp/Ga2o3YOPXfOfPzwDtzW75S86ZeCF0rL6qyzzuKNN94gKyuLnTt3Mn/+fLp06UKPHj14++23yc7O5qeffmLu3LmlOm///v157rnnOHLkCADr1q1j//79R31eEQmjjPXw6kCYcTs07QI3fAFdryvXMIcoHqGvSN+bb848d059RfreoIzShw4dyhdffMEZZ5yBmfHYY4/RqFEjLrzwQmbNmsVpp51Gq1at6Nq1K7VrB34h9pprrmHTpk107NgR5xz169fnvffeO+rzikgYZGXCF8/AnLFQ5Rj4/XNwxqUQpns9rDRTCMGUlJTkCvZDX7NmDaeeempY6imNX3/9lZo1a7Jr1y66dOnCggULaNSoUcSe159o+W8tErG2fwvv3wg/fgOnng8D/+mtZAkxM1vqnEvy91rUjtDDadCgQezZs4fDhw/z17/+NWihG6rzikgQHTkI8x+HBU/BMXXgD/+GNgW3iAgPBXoZhGp+W/PmIhFu8yKYmgIZ6+CM4dD/YahRJ9xV5VGgi4iU5NCvMPtBWPQ81I6HEW9Dy3PCXVUhCnQRkeKkzYYPbvbWl3e+Fs65D6odF+6q/FKgi4j489vP8PFfYNkkqNsSrvoQmiWHu6piKdBFRApaMw2m3wb7M6DHrdBrTMj6rwRT1N5YFCrbt29n2LBhJCYm0qZNGwYOHMi6devYtGkTxxxzDO3bt8/7euCBB/L+HBcXl/fncePGBb2u7777jjPPPJNq1arxz3/+s9Dr1113HQsWLOCOO+6gdevWtGvXjqFDh7Jnz56g1yISs37dAW9eCW9cBjUbwLWz4Zz7oyLMQSP0fJxzDB06lCuvvJIpU6YAsHz5cn766SeaNm1KYmIiy5cvz/c99913H+D1US/4WjDVqVOHcePG8d577/l9fdGiRTz77LPs37+fsWPHUrlyZe666y7Gjh3Lo48+GrK6RGKCc7DiDa+Z1uH90Oev0P3mkDfTCjaN0H3MmTOHKlWqMGrUqLzn2rdvT8+ePYP2M7Kysrj99ts5/fTTadeuHc8880xA39egQQM6d+7sd2PnNWvW0KpVK+Li4ujXrx+VK3t/T3fr1o309PSg1S4Sk/Zsgf9cDO9e5+0kNOpzOOv2qAtziOQR+odjvDuxgqnR6XDuI0W+vHLlSjp16lTk62lpabRv3x6A7t27M2HChFKXMHHiRDZu3MiyZcuoXLkyu3fvBuDWW29lzpw5hY4fNmwYY8aMKfacH374IQMGDCj0/Msvv8wll1xS6hpFKoQIaKYVbJEb6BHI35RLaX366aeMGjUqbxRdp453U8KTTz5Z5nPOnDmTV155Jd9zDz/8MJUrV+ayyy4re7EisSpjPUy9CTYvhJPOhsFPwwnNwl3VUYvcQC9mJB0qbdu2DfkuQc45v5s0l3WEfuDAAfbs2UPjxo3znnvttdeYNm0as2bN0obQIr7yNdOqDkOehfbDw9ZMK9gCCnQzGwA8DcQBLzrnHinwem1gEpCQc85/OudeKXSiCNenTx/uueceXnjhBa699loAFi9ezIEDB2jWrHR/e7/77rt89dVXjB2bf+/Afv36kZqaSu/evfOmXOrUqVPmEfqcOXPybYbx0Ucf8eijjzJv3jxq1KhRpnOKxCTfZlqtB8F5/4LjYqtfUomTRWYWB0wAzgXaAJeaWZsCh90IrHbOnQH0Bv5lZlWDXGvImRnvvvsun3zyCYmJibRt25b7778/3+g3UGlpadSqVavQ89dccw0JCQm0a9eOM844g8mTJwd0vu3btxMfH88TTzzBQw89RHx8PPv27Ss0f56SksIvv/xC3759ad++fb4LvCIVUuYhmP0QTOwN+7bBxa/BsP/EXJhDAO1zzexM4H7nXP+cx3cDOOfG+hxzN9AUL9ibA58ArZxz2UWdN5rb5wZixIgRPPnkk9SvXz+kP6djx44sWrTI7+qX4sTSf2uRIm35Ct5PgYy1Xp/y/v+IqGZaZXG07XObAFt8HqcDXQscMx6YCmwDjgMu8RfmZjYSGAmQkJAQwI+OXpMmTSqXn/P111+Xy88RiSqH93uj8i+fg1pN4LK34OS+4a4q5AIJdH9XCwoO6/sDy4E+QCLwiZl95pzbl++bnJsITARvhF7qakVESrJhLkwdDXt+gM7XeHd6RmgzrWALJNDT8aZTcsXjjcR9XQU84rz5m/VmthFoDXwVlCpFREry256cZlqvQ53EqGimFWyBBPpi4GQzawFsBYYBwwscsxn4HfCZmTUETgE2lKWgopb1SfCEa9tBkZD5bjpMuw3274Tut0DvMd4enxVMiYHunMs0sxRgJt6yxZedc6vMbFTO66nAg8CrZvYt3hTNXc65jNIWU716dXbt2kXdunUV6iHinGPXrl1Urx4dzYZEivXrTvjwTlj1DjQ8HYZPgcYdwl1V2ETUJtFHjhwhPT2dgwcPhqWmiqJ69erEx8eXemWMSMRwDla8CR/d5V0APetO6HFLVPZfKa2o2SS6SpUqtGjRItxliEgk25sO026F7z+G+M5w/nho0DrcVUWEiAp0EZEiZWfD0lfgk/vAZcGAR6DLSKgUF+7KIoYCXUQi3640byniD59Di15eM606+td8QQp0EYlcWZnw5QSY8w+IqwbnPwMdLo+ZZlrBpkAXkci0fSVMTYFty+CU87xmWrVODHdVEU2BLiKRJfMQfPYv76v68XDRK9B2qEblAVCgi0jkSF/iNdPauQbaXQL9x8KxdcNdVdRQoItI+B3eD7Mfhi+fhVqNYfib0Kp/uKuKOgp0EQmvDfPgg9Hw8yZIutprplW98F4CUjIFuoiEx8G98PFf4evXoM5J8Mfp0LxHuKuKagp0ESl/az/07vb89SdIHg1n31Mhm2kFmwJdRMrP/gyvmdbKt6FBW28ruCadwl1VzFCgi0joOQffvuWF+aFf4Ox7vTa3laNu6+GIpkAXkdDauxWm3wbrPoImSTBkPDTQfrahoEAXkdDIzvYueH7yN8g64m3Q3HWUmmmFkAJdRIJvVxp8cDNs+gya94Tzx3krWSSkFOgiEjzZWd7NQbMf9jabGDwOOl6h2/bLiQJdRILjp9Xw/o2w7WtodS4MesK761PKjQJdRI5O5mGfZlq14aKXoe0FGpWHgQJdRMoufak3Kt+5Bk6/GAY8qmZaYaRAF5HSO3wA5uQ006rZSM20IkRAgW5mA4CngTjgRefcIwVevwO4zOecpwL1nXO7g1iriESCjfNh6k1eM61OV0HfB7ypFgm7EgPdzOKACUBfIB1YbGZTnXOrc49xzj0OPJ5z/GDgVoW5SIw5uNdbU770VTihBVw5DVr0DHdV4iOQEXoXYL1zbgOAmU0BhgCrizj+UuC/wSlPRCLC2o9ymmlthzNTvFv3q9YId1VSQKUAjmkCbPF5nJ7zXCFmVgMYALxdxOsjzWyJmS3ZuXNnaWsViVmp89JYmJaR77mFaRmkzksLU0U59mfAW1fDfy+BY46Hqz+F/g8rzCNUIIHub+2RK+LYwcCCoqZbnHMTnXNJzrmk+vXrB1qjSMxrF1+blMnL8kJ9YVoGKZOX0S4+THPTuc20JnSB1e9D77th5DyIV2fESBbIlEs60NTncTywrYhjh6HpFpFSSZ2XRrv42owf3oGUycsY0TWBVxZuYlC7E0lOrFf+Be3bBtNug3Ufeq1tzx8PDduUfx1SaoGM0BcDJ5tZCzOrihfaUwseZGa1gV7A+8EtUSS25Y7OAUZ0TWDc7PUcycpm8BnlfJelc94FzwldYcNc6PcQXP2JwjyKlDhCd85lmlkKMBNv2eLLzrlVZjYq5/XUnEOHAh875/aHrFqRGJI7Mk9OrMf44R247vWlHMrMpnIlo0pcIGOtINq9AaaO/l8zrcFPQ93E8q1BjlpA69CdczOAGQWeSy3w+FXg1WAVJhLrckfm44d3AOBQZjaHM7MZ2qEJFyfF570W0mmX7Cz48jmY/RBUqgyDnoKOV0Klcv4LRYJCvzWRMMhdvZI7b/7w9NUczsymbeNazFu3M++1Fel7Q1fEjjXwUj/4+F5ocRbcuAiSrlKYRzH95kTCwHfevFereqza9gtV4ox7zzs1L+QBRvUKwbRH5mGY+yik9vSmWi54EYa/AbX9rkaWKKJeLiLlzHdVy3WvL+W3w1lUrgSVcroT5s6pr0jfG/zplq1L4f2bYMcqOO0iOPdRODYMK2kkJDRCFylnuaPzVdv2ciQrm8xsR+W4Stzev1XeWvTkxHrBHZ0fPgAf/wVePAd+2w2XToGLXlKYxxgFukg58p07f3zmOrKyHVXjDOegbePaoZk33/gZpHaHhc9Ah8u9ufJTzg3uz5CIoEAXKUe+o/NKBkeyHJUqWd7oHII4b35wH3xwC7w2CFw2XDHV29tTnRFjlubQRcpBwTXnf3xlMdnZLm/Nue/oPCjz5utmes20fvlRzbQqEI3QRUIsdV4acZXI16slO9uRme0YfEZjnr+8U/BG5/t3wdvXwuQ/QLVa3p2eaqZVYWiELhJiudMs1/c+iZTJy6hVvTKZ2Y4eLesyb91OLk6KP/rRuXOw6h2YcafXt7zXGOj5Z6hcNbhvRiKaAl0kRApOs6RMXsaxVePYtOsAPVrWY9I1XfO6Ko4f3qHso/N9P8L022DtDGjcEYaMh4Ztg/tmJCpoykUkRHxb4iYn1qPNibXY8vNvND3hGFb/uC/v+TKvbHEOvv6310wrbbZPMy2FeUWlEbpIkPkbmTeuXZ2V2/bRo2U9Vv+4L2/6JbdXS6mnWnZvhA9uho3zoFkPb/WKmmlVeBqhiwSZv5H5ym37OK1xLSZd05Xxwzvw3NwNXN/7pNKPzLOz4Itn4blk2Po1DHoSrvxAYS6ARugiQVVws4pererx+foMmtetwba9BwtNs5Rq3nzHdzA1BdIXw8n9vDCvHR+6NyNRR4EuEkS+LXF7tarHu8u2USXO+McFpwOUbZol6wh8/hTMfwyq1oQLXoDTLwbztzukVGQKdJEg8LdZxYHDWVQyqF4lDihj061ty+D9FPhpJZx2IQx4FGpqP17xT4EuEgT+NqvIynZ+N6sIKMyP/AZzx3r9V45tAMP+C60HhvhdSLRToIuUke+o3HdkXqNKXN7OQ2W6cWjTAph6E+xO85pp9XsIjjk+5O9Hop9WuYiUke9qllwHj2Tx0y+HGNqhMU9e0r50m1Uc3AfT/wyvDoTsTLjife8mIYW5BEgjdJFS8rfOvFer+kxf8SNmxtAOJzJvXUahFS3Fjs6//8TrjLhvK3S7EfrcC1WPLbf3JLFBgS5SSr7z5cmJ9ejVqj7vLttK1cqVePWqziQn1st3S3+x8+YHdsNHd8OKKVC/tXenZ9PO5fuGJGYEFOhmNgB4GogDXnTOPeLnmN7AU0AVIMM51ytoVYpEAH8j8zYn1uLz9Rm0bVyLzbsP5B1b4sjcOVj1Lsy4Aw7ugbPuhLNuh8rVyu8NScwpMdDNLA6YAPQF0oHFZjbVObfa55jjgWeBAc65zWbWIET1ioRNwZF5bpj7a7SVOyr3G+a/bPfmyr+bBie29+bKG51W7u9HYk8gI/QuwHrn3AYAM5sCDAFW+xwzHHjHObcZwDm3I9iFioRDUStZEurUYJVPb5aA5sudg2WTYOa9kHUI+v7dmy+P08ynBEcgq1yaAFt8HqfnPOerFXCCmc01s6VmdoW/E5nZSDNbYmZLdu7cWbaKRcpRUStZVm3bx9AOjfN6s5S4ufPPm+D1od6t+41Og1ELoPvNCnMJqkA+Tf7uL3Z+ztMJ+B1wDPCFmX3pnFuX75ucmwhMBEhKSip4DpGI4ztfPqJrAq8s3FS6lSzZWfDVCzDrAbBKcN6/oNOfoJJWDEvwBRLo6UBTn8fxwDY/x2Q45/YD+81sPnAGsA6RKOM7zQLkrWQZN3s91auUYiXLzrXebfvpX0HLvjD4KTXTkpAKZJiwGDjZzFqYWVVgGDC1wDHvAz3NrLKZ1QC6AmuCW6pI+Sg4zfLCZ2m8t2wrbRvXokrc//6XKXJziqwjMP9xSO0Bu76HoRPhsv9TmEvIlThCd85lmlkKMBNv2eLLzrlVZjYq5/VU59waM/sIWAFk4y1tXBnKwkWCqbiLn6u37eOe81pzbc/EkleybFvuzZNv/xbaDoVzH1czLSk3AV2Rcc7NAGYUeC61wOPHgceDV5pI+Sm4JBHgSFZ23sXPa3t6FzqLnC8/8hvMexQWjINj68El/4FTB4XjrUgFpkvsIvi/+FklrhIje57EpEWb8y5+5h6bL8x/+MIble9aDx1G5DTTOiFM70QqMl1qlwopdV5avqWIuVo3Oo5xs9dzJCub5y/vxG39Tsm3LDGfQ7/A9NvhlQGQdRgufw+GTFCYS9go0KVCKnjhc2FaBte9vpRvt+6le2Ldki9+rv8Unj0TFr8IXa+H67+AxLPL+22I5GPOhWc5eFJSkluyZElYfrYIkHeBM3eKBeD5yzv5XZKY58Bu707PbyZDvVZw/nhI6BqeNyAVkpktdc4l+XtNI3SpMApOs/iuL2/XpHZemOe+VmhUvvp9mNAVvn0Tet4O132mMJeIokCXCqOo9eVDOzRhzfZfCh2fdxv/L9vhjRHw5hVQ60S4dg787q9QpXp5vwWRYmmVi1QY+TekqMd7y7YVub4c8JppLZ8MM++GIwfhnPvhzJvUf0UilkboEpP8rWJZmJbBivS9jOiawLvLtvH7ItaXA/DzD14zrfdvgAZt4PqF0ONWhblENAW6xCR/q1hSJi8jrhJMWrSZ0X1a5jXXypWcWI9RPVvAoue9FSzpi2HgP+GPM6Bey3C9FZGAabghMangjUKTFm3m+t4n8dzcDXnTKt0S6+afZtm5DqbeBFu+hJbnwKAn4fiEcL8VkYBphC4xoaQbhUZ0TSArm3xz5Lmhv3JzBnz2L6+Z1s7v4PepcNlbCnOJOgp0iQnF3Sg0uk9LJi3anK8lbq7kGlsZufZamPV3OGUApCyG9peC+dsGQCSyacpFYoK/XizwvxuFCk2vHDmY00zraa+Z1h9ehzbnh/dNiBwlBbpEnYIbUED+FSzjZq+ne2JdbuzT0u+NQslV1nsbT+z6HtqPgP5qpiWxQVMuEnUCWcHi90ah+OqM2p8KLw+AzEMw4h34vZppSezQCF2iTplWsKyfBR/cAnu3QJeR8Lu/QbWa4X4rIkGl5lwStZ74eC3jZq9ndJ+W1KhW2e80zHcbt/CnX1+A5f+BuifDkPGQ0C2MVYscneKac2mELhGruLnydvG186ZXJi3aXLgrIpB8aCHJy2+H/RnQ889w1p3qvyIxTXPoErGKmyvPnU7xuwHFrzvgjcvhzcuhZkMYOdebYlGYS4zTlItENN+e5bkj8dwReqGR+5Y9jDr+K/jobm+Pz95jIPkmiKsSxncgElyacpGolZxYL28p4uicZYgFp1YAkuseIHnhnTB3FjTt5s2V1zs5DBWLhE9AUy5mNsDM1prZejMb4+f13ma218yW53z9LfilSkW0MC0j31x5odv7s7Phqxe8Zlqbv4RzH4erPlSYS4VU4gjdzOKACUBfIB1YbGZTnXOrCxz6mXNuUAhqlBhX1MXPD77ZxsxVPxW9FDHje6+Z1uYvILEPDH5a/VekQgtkhN4FWO+c2+CcOwxMAYaEtiypSIq6+AnFNdN6Ap7rDjvWwO+f824SUphLBRfIHHoTYIvP43TA30aKZ5rZN8A24Hbn3KqCB5jZSGAkQEKC/ucTj78bhfwtQwRIrrGN5LU3wvYVcOpgGPgvOK5hGKoWiTyBBLq/tnMFl8Z8DTRzzv1qZgOB94BCk5jOuYnARPBWuZSuVIll/i5+5nPkIMx/DD5/CmrUhT/8G9roH4oivgKZckkHmvo8jscbhedxzu1zzv2a8+cZQBUzKzy8EilCsRc/Ny+C53t6PcvbXQI3LlKYi/gRyAh9MXCymbUAtgLDgOG+B5hZI+An55wzsy54f1HsCnaxEp2Ku+NzVK/CGzTnXvx89uJWdNv4rLclXO14GPG2t5OQiPhV4gjdOZcJpAAzgTXAm865VWY2ysxG5Rx2EbAyZw59HDDMheuOJYk4RV30bBdfG4AV6XsLXfx8vfd+2rw3ABalQpdr4YYvFOYiJdCdolIu/N3x6e+iJ7/9DB//BZZNgrot4fzx0OzM8i9YJELpTlEJuxIvegKsmQbTb/OaafW4FXqNUf8VkVJQoEu5KHjRs1ti3f+F+q87YMYdsPo9aHQ6DH8TGrcPZ7kiUUmBLketrBc9x1/anuT9s+CjMXB4P/T5K3S/Wc20RMpI7XPlqJXloufEIQ2Jn3ElvHsd1GsFoxbAWbcrzEWOgi6KSlAEfNEzOxuWvgyf3AfOwTn3QedroZLGFiKB0EVRCbmALnpmrM9pprUQTjrba6Z1QrPyL1YkRinQJSiKveiZlQlfjIe5Y6FyNRgyAdpfBuavq4SIlJUCXY5akRc9h3cg+dgf4f0b4cdvoPUgOO9fcFyjcJcsEpMU6HLU/F30nHBJG6rOHwtbXoVjToCLX4O2vw9rnSKxToEufpW0FNFXwcds+YozP06BjLVwxqXQ/x9Qo055lC1SoWlpgfhV0lJEvw7vhw/HwEv9vD9f9hYMTVWYi5QTjdDFr9JsOgHAhrkwdTTs+cFbhnjOfVDtuHKtWaSiU6BLkQJaivjbnpxmWq9DnURvg+ZmyeVeq4go0KUYxS5FBPhuOky7DfbvhO63QO8xUOWYsNUrUtEp0MWvYpciNnTw4Z2w6h1oeBoMnwKNO4S7ZJEKT4Eufvlbijj+0vb8tuS/sOkJOPwrnP0X6HGL+q+IRAgFuvhVaCni3nSSF90G38+E+M7exhMNWoenOBHxS8sWK4jUeWn5N17Gm1ZJnZdW/DdmZ8Pil2BCN9j0GQx4BP40U2EuEoEU6BVEmdaV70qD1wZ7uwg16QjXL4Ru10OluHKqWkRKQ1MuFUSp1pVnZcKXE2DOPyCumje90mGEmmmJRDgFegUS0Lry7SthagpsWwannOc106p1YvkXKyKlFtCUi5kNMLO1ZrbezMYUc1xnM8sys4uCV6IES8F15fnm1DMPeSPyib1gbzpc/CoM+4/CXCSKlDhCN7M4YALQF0gHFpvZVOfcaj/HPQrMDEWhcnSKXVdedaM3Kt/5HbQbBgPGqv+KSBQKZMqlC7DeObcBwMymAEOA1QWOuwl4G+gc1AolKPytK3/24lM4Zs7fIH0y1GoMw/8PWvULc6UiUlaBBHoTYIvP43Sgq+8BZtYEGAr0oZhAN7ORwEiAhISE0tYqR6HQuvIN8+g2czT8vAmSroZz7ofqtcJRmogESSCB7m9pQ8GdpZ8C7nLOZVkxKyGccxOBieBtEh1gjRJMB/fCx3+Fr1+DOifBH2dA8+7hrkpEgiCQQE8Hmvo8jge2FTgmCZiSE+b1gIFmlumcey8YRUqQrP0Qpt0Kv/4EyaPh7HvUTEskhgQS6IuBk82sBbAVGAYM9z3AOdci989m9iowTWEePKXZPciv/RleM62Vb0ODtjBssnejkIjElBKXLTrnMoEUvNUra4A3nXOrzGyUmY0KdYFSxrs8AZyDFf8H4zvD6qlw9r0wcq7CXCRGmXPhmcpOSkpyS5YsCcvPjka5IR7Q7kEAe7d6t+yv+wiaJMGQ8dDg1PIrWERCwsyWOueS/L2mO0WjREB3eYLXTOvr1+CTv0HWEW+D5q6j1H9FpAJQoEeJEncPAti9wdvXc9Nn0OIsGDwO6rTwf0IRiTkK9ChQ7F2eifUgOwu+fA5mP+RtNjF4HHS8Qs20RCoYBXoU8Lt70PAOrEjfS3LNHd5t+1uXQqtzYdAT3l2fIlLhKNCjgL+licnNapG85UV4/p/eHZ4XvgSnXahRuUgFpkCPRulLvVH5jtVw+sUw4FE4tm64qxKRMFOgR5PDB2DOw/Dls1CzEVz6BpwyINxViUiEUKBHi43zYepNXjOtTldB3wegegk3FolIhaJAj3QH93prype+Cie0gCunQYue4a5KRCKQAj2Srf0op5nWdki+CXrfA1VrhLsqEYlQCvRItD8DPrwLVr4FDdrAJZMgvlO4qxKRCKdAjyTOeR0RP7wTDu6D3ndDj9ugctVwVyYiUSCgTaKlZKnz0vJvuox3h2fqvLTATrBvG/z3Unj7aji+GVw3H3qPUZiLSMAU6EFyVC1ul74KE7rChrnQ7yG45lNo2CbkNYtIbNGUS5Dk3o5fqha3vs20mveEwU9D3QA2rBAR8UOBHkSBt7jNgkWpMOtBqFQZBj0FHa+ESvoHk4iUnQI9iAJqcbtjDbyfAluXwMn9YdCTULtJeAoWkZiiQA+SElvcZh6Gz5+E+Y+rmZaIhIQCPUiKbXFb/Qd4/ybYscoL8XMfg2OLmVsXESkDBXqQ+G1x27QGyWlPw4vjoWZDGPZfaD0wDNWJSEUQ0FU4MxtgZmvNbL2ZjfHz+hAzW2Fmy81siZn1CH6pUWbT55DaHRaOgw4j4IYvFeYiElIljtDNLA6YAPQF0oHFZjbVObfa57BZwFTnnDOzdsCbQOtQFBzxDu6DT++DJS/DCc3hiqlwUq9wVyUiFUAgUy5dgPXOuQ0AZjYFGALkBbpz7lef448FXDCLjBrrZnrNtH75Ec5MgbPvVTMtESk3gQR6E2CLz+N0oGvBg8xsKDAWaACcF5TqosX+XfDRGPj2Tah/Kvzh3xCfFO6qRKSCCSTQ/a2rKzQCd869C7xrZmcBDwLnFDqR2UhgJEBCQkLpKo1EzsGqd2DGnV7f8l5joOef1X9FRMIikEBPB5r6PI4HthV1sHNuvpklmlk951xGgdcmAhMBkpKSontaZt+PMP02WDsDGneEIeOhYdtwVyUiFVgggb4YONnMWgBbgWHAcN8DzKwlkJZzUbQjUBXYFexiI4JzsOx1mPkXyDrkNdPqej3EaQWoiIRXiSnknMs0sxRgJhAHvOycW2Vmo3JeTwUuBK4wsyPAb8AlzrnoHoH7s3sjfHAzbJwHzXrA+ePUTEtEIoaFK3eTkpLckiVLwvKzSy07CxY9D7MfBIuDfn+Hjn9UMy0RKXdmttQ553fVheYJSrLjO5iaAumL1UxLRCKaAr0oWUf+10yrak244EU4/SI10xKRiKVA92fbMq/F7U8rvWZaAx6FmvXDXZWISLEU6L6O/AZzx8LCZ+DYBjBsMrSuWPdIiUj0UqDn2rQApt4Eu9Og4xXQ90E45vhwVyUiEjAF+qFf4NP7YfGLcHwzuOJ9OKl3uKsSESm1ih3o338CH9wC+7ZCtxugz1+g6rHhrkpEpEyiZiF16rw0Fqbl6yTAwrQMUuellf5kB3bDO9fBfy7yAvzqT2DAWIW5iES1qAn0dvG1SZm8LC/Uc/fwbBdfO/CTOAer3oUJXWDlW3DWnTDqM2jaOURVi4iUn6iZcsndozNl8jJGdE1g0qLN+fbwLNEv22H6n+G7aXBie7j8PWh0WihLFhEpV1ET6OCF+oiuCYybvZ7RfVoGFubOwbJJMPNer5lW379DtxvVTEtEYk5UpdrCtAwmLdrM6D4tmbRoM90S6xYf6j9v8pppbZgLzbrD4HFQr2V5lSsiUq6iJtBz58xzp1m6JdbN9zif7Cz46gWY9YDXTOu8J6DTVWqmJSIxLWoCfUX63nzhnTunviJ9b/5A37nWu20//Sto2RcGPwW148NTtIhIOYqaQB/Vq3Df8eTEev8L86wjsOApmPeYt/xw6ERo9wc10xKRCiNqAr1Y25bnNNP6Ftr8HgY+DjUbhLsqEZFyFd2BfuQ3mPtITjOt+nDJf+DUQeGuSkQkLKI30H/4wtt4Ytd66HA59HsQjjkh3FWJiIRN9AX64QPw6X3eKpbjm3o3CCWeHe6qRETCLvoCPa4KbP4Sul2vZloiIj6iM9Cv+RQqVwt3JSIiESWgO23MbICZrTWz9WY2xs/rl5nZipyvhWZ2RvBL9aEwFxEppMRAN7M4YAJwLtAGuNTM2hQ4bCPQyznXDngQmBjsQkVEpHiBjNC7AOudcxucc4eBKcAQ3wOccwudcz/nPPwS0K2ZIiLlLJBAbwJs8XmcnvNcUa4GPjyaokREpPQCuSjq79555/dAs7PxAr1HEa+PBEYCJCQkBFiiiIgEIpARejrQ1OdxPLCt4EFm1g54ERjinNvl70TOuYnOuSTnXFL9+vXLUq+IiBQhkEBfDJxsZi3MrCowDJjqe4CZJQDvAJc759YFv0wRESlJiVMuzrlMM0sBZgJxwMvOuVVmNirn9VTgb0Bd4FnzuhtmOueSQle2iIgUZM75nQ4P/Q822wn8UMpvqw3sPYofW5bvL833lHTs0bxe1Gv1gIyAqguPo/2dhfLc4f48lHSMPg/lf+7SniMcn4dmzjn/c9bOuaj5AiaW9/eX5ntKOvZoXi/qNWBJuH8vofydxfLnoay/c30eQnfu0p4jUj4PuV/RtifbB2H4/tJ8T0nHHs3rR/vewyWUdUf756GkY/R5KP9zl/YckfJ5AMI45SLBYWZLnK5XSA59Hiq2aBuhS2FqsyC+9HmowDRCFxGJERqhi4jECAW6iEiMUKCLiMQIBXoMM7OTzOwlM3sr3LVIeJjZsWb2mpm9YGaXhbseCS0FeoQys5fNbIeZrSzwfLG7R/lyXg/7q0NbqZS3Un42LgDecs5dC5xf7sVKuVKgR65XgQG+TxS1e5SZnW5m0wp8NSj/kqWcvEqAnw287qi5+xlklWONEgbRt0l0BeGcm29mzQs8nbd7FICZTcFrVzwWGFTOJUqYlOazgdf+Oh5YjgZwMU+/4OhSqt2jzKyumaUCHczs7lAXJ2FV1GfjHeBCM3uO6G0XIAHSCD26BLx7FIDzNhoZFbpyJIL4/Ww45/YDV5V3MRIeGqFHl4B2j5IKSZ8NUaBHmRJ3j5IKS58NUaBHKjP7L/AFcIqZpZvZ1c65TCB396g1wJvOuVXhrFPKnz4bUhQ15xIRiREaoYuIxAgFuohIjFCgi4jECAW6iEiMUKCLiMQIBbqISIxQoIuIxAgFuohIjFCgixSQs8vPeDPrFu5aREpDgS5S2CigGtAj3IWIlIYCXaSwAcA6vE0hRKKGAl3Eh5lVB+KAjsC8MJcjUioKdJH8TsYL9O+cc0fCXYxIaWjHIpH86gOt8PbjFIkqGqGL5NcYeBuoZGYnhLsYkdJQoIvkMLPKeHPnjYBUICu8FYmUjja4EBGJERqhi4jECAW6iEiMUKCLiMQIBbqISIxQoIuIxAgFuohIjFCgi4jECAW6iEiM+H8/csEGkYQ8+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "# plt.plot(Lsub_list,(LN)/np.log(2),'x',label='log neg')\n",
    "plt.plot(np.tanh(pi*Lsub_list*Tm),(LN),'x',label='log neg')\n",
    "xx = np.logspace(-1.5,0.8)\n",
    "plt.plot(xx,np.log(xx)/8+0.65, label=\"CFT, c=1/2\")\n",
    "plt.xlabel('$\\ell$')\n",
    "plt.xscale(\"log\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bipartite negativity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "δlist = np.linspace(-1,1,201)\n",
    "Tm = 5e-3\n",
    "Nc = 64\n",
    "N = 2*Nc # number of majoranas\n",
    "BC = -1\n",
    "\n",
    "Lsub_list=np.array([8,16,32,64])\n",
    "LN=np.zeros((len(δlist),len(Lsub_list))) # log neg\n",
    "\n",
    "for i_d in range(len(δlist)):\n",
    "    δ = δlist[i_d]\n",
    "    t = -(1+δ)/2\n",
    "    Δ = t\n",
    "    μ = (1-δ)\n",
    "    Gm = BdG_Hamiltonian(Nc,t,Δ,μ,Tm,BC)\n",
    "    \n",
    "    for i_L in range(len(Lsub_list)):\n",
    "        Lsub=Lsub_list[i_L]\n",
    "        inds_A= np.arange(0,Lsub)\n",
    "#         inds_AB= np.arange(0,int(2*Lsub_list[i_L]))\n",
    "        inds_AB= np.arange(0,N)\n",
    "\n",
    "        LN[i_d,i_L]=log_neg_calc(int(Lsub/2),Nc-int(Lsub/2),Gm[np.ix_(inds_AB,inds_AB)])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+VUlEQVR4nO3dd3hUVfrA8e+bRgIEQkhCCyHUEEpooYkFVBCwYEGFxYKirLsWdFddlV3Lqj/LqruKWFARexdBehORToAEQgkEAiEESGihpCfn98cMOMIkmSRzM0l4P88zT2buOefOO5PJvDn3nnuOGGNQSimlzuXl6QCUUkpVT5oglFJKOaUJQimllFOaIJRSSjmlCUIppZRTPp4OwJ1CQkJMZGSkp8NQSqkaY/369YeNMaHOympVgoiMjCQuLs7TYSilVI0hIntLKtNDTEoppZzSBKGUUsopTRBKKaWcqlXnIJRSF5aCggLS0tLIzc31dCjVnr+/P+Hh4fj6+rrcRhOEUqrGSktLIzAwkMjISETE0+FUW8YYjhw5QlpaGq1bt3a5nWWHmESkpYj8IiLbRGSLiExwUkdE5C0RSRaRTSLS06FsqIgk2cuesCpOpVTNlZubS+PGjTU5lEFEaNy4cbl7WlaegygE/m6MiQb6AfeLSKdz6gwD2ttv44F3AUTEG5hsL+8EjHbSVimlNDm4qCLvk2UJwhhzwBizwX7/JLANaHFOtRHAp8ZmNRAkIs2APkCyMWa3MSYf+NpeV6kaJzsujpz4eE+HoVS5Vck5CBGJBHoAa84pagHsc3icZt/mbHvfEvY9Hlvvg4iICPcErJSbGGPY/+hjALRbuAApxwlCpTzN8mGuIlIf+AF42Bhz4txiJ01MKdvP32jMFGNMrDEmNjTU6dXiSnlM7patFB48SOHBg5xcvNjT4ShVLpYmCBHxxZYcvjDG/OikShrQ0uFxOJBeynalapRTSxaDlxc+TZty9LPPPR2OssjAgQPZs2fPH7Y98MADtGrVqsx6rvrvf/9L586d6dKlC6NHj66Sob1WjmIS4CNgmzHmjRKqzQTusI9m6gdkGWMOAOuA9iLSWkT8gFH2ukrVKCcXLyGgZw+Cx95Jzvr15O7Y4emQVBVISUlh6dKl5Ofnc/LkyUrvb//+/bz11lvExcWRmJhIUVERX3/9tRsiLZ2V5yAGALcDm0Uk3r7tKSACwBjzHjAHGA4kA9nAXfayQhF5AJgPeANTjTFbLIxVKbfLT9tPXlISYY89RuDAgWS8/Aq5mxPx79DB06HVSs/9vIWt6ecexa6cTs0b8My1ncvd7plnnuGf//wnH3zwAVu2bKFfv36VjqWwsJCcnBx8fX3Jzs6mefPmld5nWSxLEMaY5Tg/l+BYxwD3l1A2B1sCUapGOjNyqd7FA/Bt2RKpU4e8nTs9G5Sy3JYtW0hMTOSTTz5h+fLlZSaISy65xGkv47XXXuPKK68EoEWLFjz66KNEREQQEBDAkCFDGDJkiGWv4Qy9klopi+Sn7gUR/CIjEW9v/Nq0IS852dNh1VoV+U/fChMnTuT5559HRIiOjiYxMbHU+r/99luZ+zx27BgzZswgJSWFoKAgbr75Zj7//HNuu+02d4XtlCYIpSxSsDcVnyZN8KpTB4A67dqRvW6dh6NSVlqzZg3z588nPj6e+++/n9zcXGJiYkpt40oPYtGiRbRu3ZozIzVvvPFGVq5cqQlCqZoqPzUVP4drc+q0b8+Jn3+m6ORJvAMDPRiZsspTTz3FrFmzuOKKKwA4dOgQPXr0KLWNKz2IiIgIVq9eTXZ2NgEBASxevJjY2Fi3xFwane5bKYvkp6bi18ohQbRrB6CHmWqphQsXkpeXdzY5ADRp0oTTp09z9OjRSu27b9++jBw5kp49e9K1a1eKi4sZP358ZUMukyYIpSxQdOoURUeO4PuHHoQmiNps8ODBLF++/LztWVlZBAcHV3r/zz33HNu3bycxMZHPPvuMOvZDl1bSBKGUBQpSUwHwi2hFbmEuhcWF+LZogQQEkK8JQtUQmiCUskC+PUH4RrTk3gX3ctuc2ygwhdRp04a8nZogapuxY8cSFBTktnrVhZ6kVsoC+XttCWKb/zHiM+MBmLRxEqPbteX0Wh3JVNuMHTvWrfWqC+1BKGWB/NS9eIeG8GXqjwT6BXJd2+uYtmUaJ5oGUnjgAMWnT3s6RKXKpAlCKQsU7E1FWjRj0d5F3NjuRh7p9QgASYGnAMir4IRtSlUlTRBKWSA/NZXDjX0pMkVc3+56QgJCaNOwDev8bJMS56fs8WyASrlAE4RSblacnU1hRgbpQcUE+ATQJqgNAL2b9mYpSeDlRf7u3R6OUqmyaYJQys3y96UBsKPeSaIaReEltj+z2KaxZJlsTLNQ8lI0QajqTxOEUm6Wn7oXgI1+B+kY3PHs9tgmtqkRjjepp4eYVI2gCUIpNztzkdzewFyiG0ef3R4SEELrhq1JCconPyUFU1zsqRCVcokmCKXcLH9vKkUN65HtL3/oQQD0DOvJ5rpHMXl5FKQf8FCEyt2qYsnRu+++m7CwMLp06XJe2fHjxxk5ciQdO3YkOjqaVatWVeg5zmXlkqNTRSRDRJxOhi4ij4lIvP2WKCJFIhJsL9sjIpvtZXFWxaiUFfJTUzkZWg8f8aFdULs/lHUL7cauoDxbvd27PBGeqgLuXnIUbBfZzZs3z2nZhAkTGDp0KNu3bychIYHo6Gin9crLyiuppwFvA586KzTG/Af4D4CIXAs8YoxxnPJwkDHmsIXxKWWJ/NS9HGgBbYPa4uft94eybmHdeCXEdj9v507qX3qpByKspeY+AQc3u3efTbvCsJfL3cyKJUcvvfRSp72PEydOsGzZMqZNmwaAn58ffn5+59WrCMt6EMaYZYCrc9yOBr6yKhalqkpxXh6FBw6SXO80UcFR55VHNojEu2FDshsFkLdjhwciVFY7s+TorbfeSnR0NFu2bCm1/iWXXEL37t3Puy1atMil59u9ezehoaHcdddd9OjRg3vuuYfTbrpS3+NzMYlIXWAo8IDDZgMsEBEDvG+MmVJK+/HAeLAtqqGUJxWkpYExpATmMCD4/G6+l3gRExrDvrB1NErSBOFWFfhP3wpWLDlamsLCQjZs2MCkSZPo27cvEyZM4OWXX+b555+v1H6hGiQI4FpgxTmHlwYYY9JFJAxYKCLb7T2S89iTxxSA2NhYY324SpXszCR9hxqdf4L6jG6h3djeaBlR63dhCgoQX9+qDFFZyKolR0sTHh5OeHg4ffv2BWDkyJG8/LJ7kmV1SBCjOOfwkjEm3f4zQ0SmA30ApwlCqeokPyUFgPRgSk0Qn4YCBQXk79lDnfbtqzBCZSWrlhwtTdOmTWnZsiVJSUlERUWxePFiOnXqVKl9nuHRYa4i0hC4DJjhsK2eiASeuQ8MAUrvoylVTeTvSSEn0I/gsAjq+9V3WqdrSFf2hdn+9HL1PEStYeWSowCjR4+mf//+JCUlER4ezkcffXS2bNKkSYwZM4aYmBji4+N56qmnKv18YGEPQkS+AgYCISKSBjwD+AIYY96zV7sBWGCMcTyj0gSYLiJn4vvSGON8bJdS1UxeSgoHgks+vARQ368+AW3bUeSVRF7SDrj66iqMUFll8ODB3Hvvvedtz8rKcsv+v/qq5HE83bt3Jy7O/VcEWJYgjDGjXagzDdtwWMdtu4Fu1kSllLXydqeQ0jKfaCcnqB11adad9JCdNNiRVEWRKVV+eiW1Um5SlJVF8dGjpJfRgwDbeYjdYcWcTtyMMTq2oqarrUuOaoJQyk3OnqBuzB/mYHKme2h3kpsLHD5K4cGDVRGespAmCKVUqfLsM7TmNA8mJCDk94LN38OKNyEr7eymVg1acbBVoK1+QkJVhqmUyzRBKOUm+SkpFHlBSFuHIYYpy+DHe2Hh0/BmN9izHAARIbhLTwp8hJyETR6KWKnSaYJQyk1ydydzKEjoEGpPEKcPw/fjILgt3LcCGjSHOY9BUSEAseH92NXEkLVhnQejVqpkmiCUcpNTyTvY39jhArm4j+F0Btz8MTTtAkNehIytEGcbv96veT92thAKtm3HFBR4MHKlnNMEoZQbFOfnQ2o6+0KxDXE1BuK/gMhLbDOCAkRfa3v82+tQmE/7oPYcbBWIV34huTovk6qGNEEo5Qb5yclIcTEHm/kTHhgOqavgWAp0H/N7JREY8DCcOgRbpiMiNIy1zZ+TvU4PM6nqRxOEUm5wpgfg074dXuJl6z341YdO1/2xYtvLIaQDrH4HjCGm0yDSgyHjt8UeiFqp0mmCUMoNsrdvJd8HmnbsAUUFsO1niL4O/Or9saKXF/T9MxyIh31ruKj5RWyOFIrWJ+h5iBrM6iVHc3Nz6dOnD926daNz584888wzZ8v27dvHoEGDiI6OpnPnzrz55psVeQlOaYJQyg2ytiSQGgIxTbrD3hWQmwXR1ziv3G00+DeE1e/SpF4TjnVpiXdeATmbdLhrbeHuJUfr1KnDkiVLSEhIID4+nnnz5rF69WoAfHx8eP3119m2bRurV69m8uTJbN26tdLPCdVjum+lajRjDIU7ktnbWhgR0gV+ewt8AqDNIOcN/OpBzzth1WQ4vo/wy4ZS/OkUMn9bQqtevao2+FrklbWvsP3odrfus2NwR/7R5x/lbufuJUdFhPr1bbMDFxQUUFBQgH1CU5o1a0azZs0ACAwMJDo6mv3797tlym/tQShVSYWZmficyCazWV3C67WA7XOg7SDwq1tyoz73AgbWfcAl0cPY3QwO/6rnIWoDq5YcLSoqonv37oSFhTF48OCzCwQ52rNnDxs3bnRaVhHag1CqkvKSbDOy+nXsgBzaDCfSYOATpTcKioCO18D6T4i69HHmRjWgza97KTx6FJ/g4CqIuvapyH/6VrBqyVFvb2/i4+M5fvw4N9xwA4mJiXTp0uVs+alTp7jpppv43//+R4MGDSr1Gs7QBKFUJWUlbAAgLKavrfcgXhA1rOyG/f4K22Yim74h4PLL8Fr6M4cXzqXprWPKbquqpapYcjQoKIiBAwcyb968swmioKCAm266iTFjxnDjjTe658Vg7YJBU4FrgAxjTBcn5QOxrSSXYt/0ozHm3/ayocCbgDfwoTGmeqxGrpQTRzasJiMYOrfqDT8/Di37Qr2QshtG9INm3WDN+1x08xQyGv7MqVnfa4KowaxacjQzMxNfX1+CgoLIyclh0aJF/OMfth6TMYZx48YRHR3N3/72t8q/CAdWnoOYBgwto85vxpju9tuZ5OANTAaGAZ2A0SLingVWlXIzYwxmyw6Smwtd/ILg4GaIGu5aYxHo+xc4nETXrAy2dw7Ef2MSRadOl91WVTtWLjl64MABBg0aRExMDL1792bw4MFcc41tlNyKFSv47LPPWLJkydnzF3PmzKnU851h5Ypyy0QksgJN+wDJ9pXlEJGvgRGAe8ZtKeVGhenp+GVlc2JQUxrsWmbb2LEcS4h2uQmWvIAs+w91r7gcn5UzOLB4NuEjbrEmYGUZK5ccjYmJYePGjU7LLr74YssWnfL0KKb+IpIgInNFpLN9Wwtgn0OdNPs2paqdUwnxAAR27wlJsyEkChq3dX0HPn5w8cOQtpaLusRwJBBSv/3UkliVKi9PJogNQCtjTDdgEvCTfbs4qVtiehSR8SISJyJxmZmZ7o9SqVLsX7uUfG9o26UP7FlRvt7DGT1uh8BmtN3wOUm9m9Bg4y4K9LNco+iKcm5mjDlhjDllvz8H8BWREGw9hpYOVcOB9FL2M8UYE2uMiQ0NDbU0ZqXOdTo+npSm0PP0UTBF0GlE+Xfi6w8XPwJ7V9DsiovwLoYtX73n/mCVZTRBuJmINBX7pYAi0sceyxFgHdBeRFqLiB8wCpjpqTiVKklxfj51d6WTEdmQkJ2LIaiVbVRSRfS8A+o3YdDJeHY39yb7p58tO66slKssSxAi8hWwCogSkTQRGSci94nIffYqI4FEEUkA3gJGGZtC4AFgPrAN+NYYU/qliEp5QHZCAj4FxUhMFOxeaus9iLMjpC7wDYABDxOwZwVZl3WgUfpJDixfVHY7pSxk5Sim0WWUvw28XULZHMA947SUssieX2cjQETbJhBXULHDS456jYWVkxjgn0GGPxz64H80v2SwO0JVqkI8PYpJqRrrxJpVpIZB7xMp0CAcWlRyoj2/unDls7TK3Exy38aErNvNiX273ROsUhWgCUKpCjD5+dRLSuNgu2CCd/1qWxiolMNLxcWGzJN5Ze+4683QvCfRLY8AsOGtf7srZKXKTROEUhVwLGE9fvnF+EU1haI82+JATuTkF/GP7zcR89wCer+4iGsm/cbsTQdK3rGXFwx9me6F6WzvGkDQvLWcPlTiID6lLKUJQqkK2L1kBgBtGudB/aa2+ZfOceRUHiPfW8m36/cxvGtTHrsqivzCYu7/cgPPzEikoKjY+c4j+kKXm2jdPhPfAsOa/0608qUoVSJNEEpVQM6KVaQ086LHoXUQfa3tP38HxcWGv32bwM6MU3x0ZyyvjuzG/YPaMeehS7jn4tZ8smovf/5sPTn5Rc6f4Mpn6RVQwI5OvjSas5qs9D3WvyhVYVYvOQpw/PhxRo4cSceOHYmOjmbVqlV/KC8qKqJHjx5n52hyB00QSpVT4YkTNErO4EhUEP4FuRBz63l1pq5I4dcdmfzrmk5c3rHJ2e0+3l7885pOvHhDF35JyuCOqWvIynGyFnVQBFzxL9p3OIRPIax5wb2zdCpruXvJUYAJEyYwdOhQtm/fTkJCAtHR0X8of/PNN8/bVlm6HoRS5ZSyeAbexdCoRTE0bAPhsX8oP5iVy2sLkrgyOozb+kY43ceYvq0ICvDj4W82cuv7q/h0XB/CAv3/WKnvfXTZ8hM/dt1H1C/bSN20koiYi6x6WTXewf/7P/K2uXfJ0TrRHWn61FPlbufuJUdPnDjBsmXLmDZtGgB+fn74+fmdLU9LS2P27NlMnDiRN954o1LP5Uh7EEqV0/7Fs8iuA72Kd9h6D+eMXnpjYRJFxYanr+l8dt1gZ66OacbUsb1JPZrNyHdXkXok+48VvLxhxGT6RWWR4wdJ/3pMr66uAaxYcnT37t2EhoZy11130aNHD+655x5On/59WviHH36YV199FS8v936law9CqXIwxhAQt52UNnXoZQoh5o/Tcm8/eILv1qdx94DWRDQuZU1qu0vah/LFPX25a9o6bnpvJZ/e3YfoZg7LRYZ2oPmQx4lPf4vWy4+yatorXHRXGcuZXqAq8p++FaxYcrSwsJANGzYwadIk+vbty4QJE3j55Zd5/vnnmTVrFmFhYfTq1YulS5e66VXYlJggRMSVdety7Vc9K3VByEhYQ4Pj+UgfX9vIpeA2fyh/ac52Auv48ODl7VzeZ4+IRnz35/7c/tFabnl/FR+P7U1spMO61P0f5Mrts1m6az9Bkz7l6NBRBDeLdNMrUu5k1ZKj4eHhhIeH07evbbTcyJEjefll20KbK1asYObMmcyZM4fc3FxOnDjBbbfdxueff17p11NaD+IDbEuClja5zKXolBjqApL006c0FugUdABiJvyhbPnOw/y6I5OnhnckqK5fCXtwrn2TQL7/S3/u+Ggtt320hnfH9GJQxzBbobcPfjd/Qos9l1E03Zc1f7uLoV8uKfXwlfIMq5Ycbdq0KS1btiQpKYmoqCgWL15Mp062hTZfeuklXnrpJQCWLl3Ka6+95pbkAKWfg5hrjLnbGHNXSTdsk/EpdcGQZWtJaelNRx8DnW84u7242PB/c7bRIiiAO/pHVmjf4Y3q8u19/WkXVp97P43jp437fy9s2ILOY98nNTafyI0H+e2D5yr5SpS7WbnkKMCkSZMYM2YMMTExxMfH81QVHFIrMUEYY24rq7ErdZSqLY4kbyEk/TT5rYqQDldB3d8PA03fuJ+tB07w+NAo/H29K/wcIfXr8NW9/YiNbMTD38Tz8YqU3wvbXcFVt9/J3vBiGrz1DUlr5lXm5Sg3Gzx4MMuXLz9ve1ZWFsHBwU5alE/37t2Ji4tj06ZN/PTTTzRq1Oi8OgMHDmTWrFmVfq4zSj3lLSINROS89RNFpPSDakrVQlumfwxAu9Cj0OP3/41yC4p4fUESXVs05NqY5pV+nkB/X6bd1YchnZrw3M9beWNB0tnRS76XT6THyA5k+0PmhEc5lpFa6edTqiQlJggRuQXYDvwgIltEpLdD8TSrA1OquileuIzUpkKXkMbQ7vdpuKeuSCE9K5enhkfj5eWe8wL+vt68M6YnN/cK560lyfxrRiJFxQa8vGly11fUuzaQBieLWHPvzRQW5LvlOVXFXYgryj0F9DLGdAfuAj5zGNmkZ8fUBeVQUgJNUk+S3ToPr553gLdtfMeRU3m8+8surugYRv+2jd36nD7eXrw6MoY/X9qGz1enMv7TONtV13Xq0+1vP3PoMi9aJZ1g/iMjKS4uYV4nVSUuxAThbYw5AGCMWQsMAiaKyENAmVfriMhUEckQEaeDgEVkjIhsst9Wikg3h7I9IrJZROJFJK58L0kp99vyjW2N6M7NT0GP289un7QkmdP5hTwxrKMlzysiPDk8mueu68yvOzK5fvIKkg6ehMAmXPnsjyR3M7RZtJP5/9TTgcr9SksQJx3PP9iTxUBgBNDZhX1PA4aWUp4CXGaMiQGeB6acUz7IGNPdGBN7flOlqo4xBp8lq0kJhw6dB0HDFgDszjzF56v3cmvvCNo3CbQ0hjsviuSr8f04lVfIDe+sYNamdCQsiuGvfc6uaEPkjxtZ8OI9lsZQXenV5a6pyPtUWoL4C+ccSjLGnMT2pX+3C8EsA0oc22WMWWmMOWZ/uBoILzNapTxgz7olhB7MpTgyF2LvOrv95bnbqePjxSOD21dJHL0jg5n14MVEN2vAA19u5MXZWylq1pOr3vyc3e2g5WcrWPyfv1ZJLNWFv78/R44c0SRRBmMMR44cwd/fv+zKDkq8UM4Yk1DC9gLgi/KFV6ZxwFzHpwEWiIgB3jfGnNu7OEtExgPjASIinE+MplRl7PjifZp6Q2z7utDOdmXr6t1HWLD1EI8O6XD+JHsWatLAn6/u7ccLs7fywW8prNx1hP/d2p0r3/6EJffdSeRHvzDvyM0M+b9v3D4vT3UUHh5OWloamZmZng6l2vP39yc8vHz/h0tZmVdETnL+OYcsIA74uzGmxEVzRSQSmGWM6VJKnUHAO8DFxpgj9m3NjTHpIhIGLAQetPdIShUbG2vi4vSUhXKfgtwcEvr3Yn/LIkb8fRxc9jjFxYYRk1dw+FQeS/4+kAC/il/3UBnztxzkyR83czqvkCeHdWRUq5MsenAUbbcZdl/SjKHvLcDbW6dbU6UTkfUlHcp35V+MN4DHgBbYDgM9im0ajq+BqZUMLAb4EBhxJjkAGGPS7T8zgOlAn8o8j1IVtWH6FOrlGIJb50Ev2+GlGQn72bw/i8euivJYcgC4qnNT5j18CRe1bcyzP2/l3vk59Hx7Hrt616HNbweYP7IvJw7rcqWq4lxJEEONMe8bY04aY07YD/cMN8Z8A5x/KZ+LRCQC+BG43Rizw2F7PREJPHMfGAKUPh2iUhbJ+uYbjtWHvgOHQ/1QcvKLeHWe7aK467u38HR4hAX6M3Vsb164vgvr9hxlyJRtHBn/DalXhRCxPZuN111J8qoZng5T1VCuJIhiEblFRLzsN8f5jUs8PiUiX2GbqylKRNJEZJyI3Cci99mrPA00Bt45ZzhrE2C5iCQAa4HZxhidU0BVudTE1bTcfowjnfLxG/AgAB8t382BrFwmXu2+i+IqS0S4rV8r5j98Kb1aNeKZ2cm81eZFMu67FP9cw8nxT7D0f/frtRKq3Fw5B9EGeBPob9+0CngE2I/tQrrzJx/xED0Hodxp7kM3Eb5oK40fiKT5X+eScTKXQf9ZyoB2IUy5o3qOvjbGMGfzQZ77eQuZp/KY0Go3nb57h5b7YXenAPq9+QWNW7p3WUpVs1XqHIQxZrcx5lpjTIj9dq0xJtkYk1OdkoNS7pR94ihhv24lpX0RzYc8CsDr83eQV1hs2UVx7iAiXB3TjMV/v4yxF0Xy1t42TLzoBTZf2oiI7TnsvP5GVrz7iPYmlEvKTBAiEi4i0+1XRR8SkR9ERK9ZULXamo9foW4etOhWH9peTvy+43wTt4+7L25Nm9D6ng6vTIH+vjxzbWcWPHIpPdq34vHgiXx8zWCy60Lwm/OYf1139qxz36yfqnZy5RzEx8BMoDm2kUw/27cpVSsVFxcjP8xlfxNDj+smUGzg6RmJhAXWKddKcdVBu7BAptwRyw9/6U9q61E82P8ZlvdvTJPUArLGPsbsvwwk6+AuT4epqilXEkSoMeZjY0yh/TYNCLU4LqU8ZvOCL2mSUYDp6oVX99F8G7ePTWlZPDU8mkB/X0+HVyG9WgXz3X39eXPsIGb1epl/Dr6D5LY+tPnlEDuGXsO8x4aTffyQp8NU1YwrV9EcFpHbgK/sj0cDR0qpr1SNlj7lbRoHQP9b/0JWvvDq/CR6RzZiRPfKr/XgSSLClZ2acEV0GCt3deSdpYOYv+kzrk36hXY/p7BpyUCOXdqSfn//H43CO3k6XFUNuJIg7gbeBv6LbVjrSlyYi0mpmmjnuoVEbs1iVz/o0+9enpmTxPHsfJ67rm+tWQNaRBjQLoQB7UKI39eRKb8m479mEoOS1tFu7j5SFt/EithG9HjoOVp0H1z2DlWtVeYw15pEh7mqypoz+mKabjlC5GtjOdjlL1wz6Tdu79eK50aUOFtMrXDoRC5frU0leeH7XLptAZ12FgKwp70vIVdfQa87X8DXv56Ho1RWKG2Ya4kJQkQmUcqFcMaYh9wTnvtoglCVsXfzKk7dfDe7e8E1H6/n1qkJJGee4pe/D6Rh3Zp57qG8CoqKWbj1EL8u+J7otZ8QvfM0DbLheH3I7BlCzN2PENnvxrJ3pGqM0hJEaYeY9JtWXVASX3+KcB/oM+oWvk04zNo9R3nlpq4XTHIA8PX2YnjXZgzv+iCHT43n57XbODrzRdpt30q73w6Ts2wii8Imkt29OTF/ul+TRS2nh5iUAtKT4zly3Wj2dIV+767kirfW0rFZA76+t1+1mVLDk3ZnnmLJgumYJR8SsSuTiIO27430UMjqGELEkGvoft1D+NQJ8HCkqrwqeojpWWPMs2XsuMw6VUkThKqo2fdcTsTKAwS98ideOXk1C7YcYu7Dl9C2BlwUV9UOZOWwaM735C+ZSovkQ7Q8aPAycMof9kf64d0tii7X3EGrnsPw8vbcbLfKNRVNEGnYpvoucb/AvcaYajPvgCYIVRGpW9eQNXIse7t4Ue/ZBdz1aTx/G9yBh66ompXiarLcgiJWrlnO3p8nU39HEhF782mQbSs7FgiHwv0w7SNoe8UNdLn8Nrx9/TwbsDpPRRPEMy7s+5Qx5vXKBOdOmiBURcwe3Z8Wm48T9MZfGbOhM3X9vJn90CX4+dT+Fdnc7dDxk6yYMYVT6xZQf286zdILCTptKzvlDxlNvDjdrAF12nWgw8DriOp/vfYyPKxCCaIm0gShymv7yp8puvtxUvr6sv6mL5m2ai/f39efXq2CPR1arXAqJ4/VC77h4G8/4rdnD8EZeTQ5DD72uQJP14GMMC9OhdaF5k0JiupOx0uuI6Jjb88GfgHRBKFUCeZe25PQ1Bx8X3uS65c05ra+rXj++tp9zYOnpR9IZ+P8Tzi+aQW++/YTlJFH6FGDf8HvdU7UhaPBwumgOhQ1boBvs3CC2nWhXeyVhLfvqb0ON6roMFelarWNM94hcmcOey5vyH8TWxEWWMBjQ6M8HVat17xZc5qPffIP23LzcklcPZ/UuAXk7U7C99AR6h/NIzwll4aJuUAGsIHTfMpGPzjWEE438CG/vh/FQfXxCg7Bv2kEjVtH07LzRTRv1RFvH/16qyzLehAiMhW4Bsgwxpz3L5nY5i14ExgOZANjjTEb7GVD7WXewIfGmJddeU7tQShXFRUV8suV3al3soj1j77Ea/G+fHxXbwZFhXk6NOXAGMP+fbvZuW4BR3dsID99D96Hj+F/PJfAE0UEnoZ6uee3y/OBE/XhdD0hL8CbwgAfiur6Y+rVxathQ3yDwqgb2pyGzdsQFhFFs1ZRBNRrUPUvsBpwew9CRK4xxpQ1mfw0bHM4fVpC+TCgvf3WF3gX6Csi3sBkYDCQBqwTkZnGmK0ViVUpZ359/a+0OFDE7pta83qCL6P7RGhyqIZEhPCItoRH/MVpuTGGzEP7Sdm8nMO7N5OdnkLh4Qy8jp/A71Qu/qcKaXiskIADhdTLycWn+DiQDmw7u488YA+2pJLrB3l1IN9PKLDfCv28KarjQ3EdP0wdP8TPD/Grg9Spg1edALwD6uETUA/fgHr41WtInfoNCQgMpm7DxtRr0Ji69QPx969HnXoN8PGtY/2b5kYV7YP1BkpNEMaYZSISWUqVEcCnxtaFWS0iQSLSDIgEko0xuwFE5Gt7XU0Qyi2OZ+yj3te/sa8FvNXgr7QggIlX6zKcNZGIENY0nLCmo4BRpdYtLCjk0P5kDu7azLH0XZzOTCP/aAZFJ45Ddg7k5uOdV4BPfiE++cX45hdT/0QxdfKK8CvIJyAv++zJdVcY4JT9dkaxQJHX77fiMz+9f39c7AVG7PMcnf0p5zy21cF+P8/fixGzE10PzkUVShDGGFeGwJalBbDP4XGafZuz7X1L2omIjAfGA0RERLghLFXbLX/qdlpnw7ZbBpF0TPh6fDfq19Hj1bWdj68PLSI70iKyYpduFRcVcfrEUY4dPsjprExOZx0h99Qx8k+fIO9UFoU5JynMyaYo9zTFeTmYvDxMfj6mqAiKiqCoGCkuQoqK7feLkWJj+1lk7PcNXvafwNnZ8ATAGOTMGQFzZpvtp7FopuEy/ypExNlkK1nAZmNMRiWe29krMqVsd8oYMwWYArZzEJWIR10Adq74gcgVh9jVzZdnjw/n7gGt6demsafDUjWAl7c3gY1CCWx04ayX5sq/TeOA/sAv9scDgdVABxH5tzHmswo+dxrQ0uFxOLaDg34lbFeqUoqLi9n93DM0rgPfdhxHmwb1eVxHLSlVIlcuFS0Goo0xNxljbgI6YTuv0xf4RyWeeyZwh9j0A7KMMQeAdUB7EWktIn7YDizOrMTzKAXALy/fRURqEVsGNGdVQXveuKU7/r46nl6pkrjSg4g0xjguVpsBdDDGHBWRgpIaichX2HobIfZ5nZ4BfAGMMe8Bc7ANcU3GNsz1LntZoYg8AMzHNsx1qjFmS3lfmFKODu3aQNA3a9nbUnix/v38fUgHurcM8nRYSlVrriSI30RkFvCd/fFIYJmI1AOOl9TIGDO6tJ3aRy/dX0LZHGwJRCm3iHt0HOGFML3bUPq0b8F9l7b1dEhKVXuuJIj7gRuBi7GdQP4E+MH+BT/IwtiUcosV7z5Cm225rOnXgLiGw5h7Szdd40EpF5SZIIwxRkTisJ0jWCQidYH6wEnLo1Oqko6kJuL94TzSQ+E/YY/w9i3dCGvg7+mwlKoRyjxJLSL3At8D79s3tQB+sjAmpdyiuLiY1RNup14OfN3zCkZfGqNXSytVDq6MYrofGACcADDG7AT0r0xVe8v+e5/t0FLvBhxsN4rHh1abta2UqhFcOQeRZ4zJF/uVeiLiQykXrilVHaQnLiPw09/Y2xymRD7JD2N66gJASpWTK38xv4rIU0CAiAzGNprpZ2vDUqriCvNy2fToX/Eqho973Mirf7qIlsF1PR2WUjWOKwniCSAT2Az8Gdvw039aGZRSlbHwiRG02lPE4j4tGDB8LAP1vINSFeLKKKZi4AP7TalqLf6H/xAxL5Ut7b3ZMuAFPhrUztMhKVVjuTJZ3wDgWaCVvb5gG/3axtrQlCqfo3s3kf3yVAobwpd9H+PjUb30egelKsGVk9QfAY8A64Eia8NRqmKKCvJZff8YWp6GyVddwQt330DDur6eDkupGs2VBJFljJlreSRKVcL8vw2ldXIh8/o1ZcTdE+nU/MJcPlIpd3IlQfwiIv8BfsQ2iysAZ9aPVsrTVr3/N1otPMCmKF98xrzLsK7NPB2SUrWCKwnizGpujotaG+By94ejVPnsXj0dv8lzSQ+DuBte49UrdX0HpdzFlVFMOiGfqpZOZuwl9e8TqesNs64axytjrkAsWnpRqQuRXlqqaqTiwkKW3XsdoUcNcwb2YuJDD+jiP0q5mSYIVSPNeehK2iTls7RvKLc89g5hgTpDq1LuZmmCEJGhIpIkIski8oST8sdEJN5+SxSRIhEJtpftEZHN9rI4K+NUNcvil26n7ZJDJHT0o9uT3+mIJaUs4sqFcjc62ZwFbDbGZJTSzhuYDAwG0oB1IjLTGLP1TB1jzH+A/9jrXws8Yow56rCbQcaYwy69EnVBiPvi34R9FseucMHviS+5OKqJp0NSqtZyZRTTOKA/8Iv98UBgNdBBRP5tjPmshHZ9gGRjzG4AEfkaGAFsLaH+aOArF+NWF6Dk5d/Cq19xOAgOPPBf7u3X2dMhKVWruXKIqRiINsbcZIy5CeiE7XqIvsA/SmnXAtjn8DjNvu089lXqhgI/OGw2wAIRWS8i40t6EhEZLyJxIhKXmZnpwstRNVHmzrUc+PszFHlDwu0PcM+IIZ4OSalaz5UEEWmMOeTwOAPoYD8UVFBKO2fjDUtaR+JaYMU5h5cGGGN6AsOA+0XkUmcNjTFTjDGxxpjY0NDQUsJRNdWpw3vZOH4sgadh2Y3X8OCf/6rDWZWqAq4cYvpNRGZhWwcCYCSwTETqAcdLaZcGtHR4HA6kl1B3FOccXjLGpNt/ZojIdGyHrJa5EK+qRfJOHmHZmOG0OmhYcE0vJjzxCt46AZ9SVcLVJUc/BroDPYBPgPuNMafLuIhuHdBeRFqLiB+2JDDz3Eoi0hC4DJjhsK2eiASeuQ8MARJdekWq1ijMy2H+bYNovbeYpZe35S//94muCqdUFXLlSmojIsuBfGyHiNYaY8pcctQYUygiDwDzAW9gqjFmi4jcZy9/z171BmCBMea0Q/MmwHT7YQQf4EtjzLxyvC5VwxUXFTHrzkuISipg+cVh3PH6dL0QTqkqJmV914vILdiGoi7Fdl7hEuAxY8z3lkdXTrGxsSYuTi+ZqPGMYfq4AXRceYy1sYFc9/5yGtXz83RUStVKIrLeGBPrrMyVcxATgd5nrnkQkVBgEVDtEoSqHX66fxAdVx4jvqs/V739iyYHpTzElQO6XudcEHfExXZKldv0+y8naskhEqP8GPj+rzQNqufpkJS6YLnSg5gnIvP5fZTRrcAc60JSF6rpDw6m4+IDJEb5MuCjpTQL1ik0lPIkV05SPyYiNwEDsJ2DmGKMmW55ZOqCMn3CEDouTGNrlA8DPlpK85BGng5JqQueKz0IjDE/8MernJVym+kPXknHhfvZ1sGH/lOW0Dwk2NMhKaUoJUGIyEmcX/ks2Ea/av9fVUpxURHT7x9Ep6WZbG3vy4CPltA0NMTTYSml7EpMEMaYwKoMRF1YiouKmD5uAJ1WZ5HY0Y9BU38lJDjI02EppRy4dIhJKXcqKshnxh396bQxm4SuAQybtpyG9ep6Oiyl1Dl0uKqqUvk5J5kxqg/RG7PZ0CuQ6z5brclBqWpKE4SqMiePHGDuzRcRvSWPuIsac9PHK6nrrxfBKVVdaYJQVWL/zg2suPkK2iUXsv6KcG776Df8/fQIp1LVmSYIZbltK2ey884xND9kSLilF7dNXqjrOShVA+i/cMpSq6dPghfeoX4hJP/5Ov700CueDkkp5SJNEMoy8yf/jdD355JTB07980FuvPmvng5JKVUOmiCUJb7/581E/ZhIZjA0eOl1Blwy3NMhKaXKydJzECIyVESSRCRZRJ5wUj5QRLJEJN5+e9rVtqp6KizI49s7+9H5+0RSW3jR+uPviNXkoFSNZFkPQkS8gcnAYGzrU68TkZnGmK3nVP3NGHNNBduqaiQzLZlVf76erruK2BwTwLAPlhDYMMjTYSmlKsjKHkQfINkYs9sYkw98DYyogrbKAxKXz2DzqGtpu7uIjYMjuOnLtZoclKrhrEwQLYB9Do/T7NvO1V9EEkRkroh0LmdbRGS8iMSJSFxmZqY74lbltGjq05x+6AkanoSke67kT5Pm4+2jp7eUqums/Ct2NtD93NlhNwCtjDGnRGQ48BPQ3sW2to3GTAGmgG1N6gpHqyrk2yduJHrmNo42BJ58nBuvu8vTISml3MTKBJEGtHR4HA6kO1YwxpxwuD9HRN4RkRBX2irPOnU8k7l/voquCTnsjvCi09tf0LpDd0+HpZRyIysPMa0D2otIaxHxA0YBMx0riEhTsV9SKyJ97PEccaWt8pytq2ez8vrL6JKQQ0Kfhlw+fZUmB6VqIct6EMaYQhF5AJgPeANTjTFbROQ+e/l7wEjgLyJSCOQAo4wxBnDa1qpYlevmTn6Uxh/OJqQItozpy6h/TfN0SEopi4jt+7h2iI2NNXFxcZ4Oo1YqKizk+weH0GXpATKCwe/JfzDgmrGeDkspVUkist4YE+usTIeaqDKlJiew4eHbiEkuZFsHX/q9/SPNI9p5OiyllMU0QahSLZz2AgHvfEHbUxA/uCW3/HeODmFV6gKhf+nKqcKCPL5/aBidfz3A8UA4+MTtjL7jKU+HpZSqQpog1HlSNq9ky6P30m1vMdvb+9L7za8Jb9PJ02EppaqYJgj1B7MmP0bI1Fm0zIX44W255dWf9JCSUhco/ctXAGQdzWDehGvpsu4EmcGQ+8T9jL75AU+HpZTyIE0Qit9+nEzBG28Tcxi2xARw2f++J7R5G0+HpZTyME0QF7CC/Dx+eORqOi3dT3YdSBo3kJGPvevpsJRS1YQmiAvU5uUz2ffsk3RLK2ZHWx86v/w+/bpe5OmwlFLViCaIC0xRYSE/PPMn2v68mabAphEdGfnid3oiWil1Hv1WuIAkbVjKtqcfpGtyIXtbCI0n/ptbLx/p6bCUUtWUJogLwJleQ+vZm2ldAAmDmnLDG7OpE1DX06EppaoxTRC13JbV89j13KN0TSlib3Mh6NHHGDVcF/VRSpVNE0QtVVRYyPdP3kS7+TuIMBA/JIIbX5muvQallMs0QdRCG3/5gfSXniYmtZiUll6EPfEso6+42dNhKaVqGE0QtUj26RPMePxGOi7bTzOxTZUx8uUf8PWr4+nQlFI1kJVLjiIiQ0UkSUSSReQJJ+VjRGST/bZSRLo5lO0Rkc0iEi8iugpQGRZMe4FVQ/vSffF+9oV74/Puq4x+Y5YmB6VUhVnWgxARb2AyMBhIA9aJyExjzFaHainAZcaYYyIyDJgC9HUoH2SMOWxVjLVB2u5EVj15J50SssmqD1tv78cNT3yIl7e3p0NTStVwVh5i6gMkG2N2A4jI18AI4GyCMMasdKi/Ggi3MJ5apaiwkJ+ev4MWszbSKRu29KzPpS9/yUUR7T0dmlKqlrDyEFMLYJ/D4zT7tpKMA+Y6PDbAAhFZLyLjS2okIuNFJE5E4jIzMysVcE2xduHXzL26G52+2cjJ+sKhp8dxy5fraKrJQSnlRlb2IMTJNuO0osggbAniYofNA4wx6SISBiwUke3GmGXn7dCYKdgOTREbG+t0/7XF4QMpLHpqDNFrj9HcG+KvasUNL/2Af916ng5NKVULWZkg0oCWDo/DgfRzK4lIDPAhMMwYc+TMdmNMuv1nhohMx3bI6rwEcSEoKixk+gtjaTZ7Pd1OwrZoPzpNfIPRsVd4OjSlVC1m5SGmdUB7EWktIn7AKGCmYwURiQB+BG43xuxw2F5PRALP3AeGAIkWxlptLZ/+LvOHdaPz1+vJ8YfUR2/hxukJdNTkoJSymGU9CGNMoYg8AMwHvIGpxpgtInKfvfw94GmgMfCOiAAUGmNigSbAdPs2H+BLY8w8q2KtjvbtSmTVP8fSKf40depAwvA23PDid3oltFKqyogxteewfWxsrImLq9mXTOScPsGMp8fQ+pdk6ufYVnjr/ewUWkfHejo0pVQtJCLr7f+Yn0evpK4mjDHM/N8EAn9YSLfDsKeFkPOP8dxy68OeDk0pdYHSBFENLJ/xAcfe/R8d9hRzuCFsu70fI/7xgS7io5TyKP0G8qDkzSvZ+MIDRG/KIcAPEq4M5+rnvySwUainQ1NKKU0QnnAscz/zn76d9isOEF0AW7vXo/e/3mFUpz6eDk0ppc7SBFGFck6fYOazt9Pylx10OwU7W3sT8sDfueVqXcBHKVX9aIKoAkWFhfz00jgaz11LzFFIbSYcGzeC6/7ykqdDU0qpEmmCsNistx/F77s5dDpkOBgM2+64iBGPv68noJVS1Z5+S1lkyZevkT3tY9qmFnM0EDZd35Hrn/lCL3RTStUYmiDcbOXsjzn0/n/puKOAkwEQf2ULhj33GQMaN/N0aEopVS6aINxk7YLPSXv3P0RtzyfABxL6B3PJ0x/Qp3UnT4emlFIVogmiktYv/o4977xI1NY82nhDYmwDYh//L6O6XuTp0JRSqlI0QVTQxl+ns2vSv+m4JZd23rYV3bo/9jq3dr/U06EppZRbaIIop00rZ5P0v3/RMTGH9mK7yC3m0Ve5tdflng5NKaXcShOEizb+8gPJ775I1JYcOhrY3rUu0Y+8wC39hnk6NKWUsoQmiDKsmvUR6R+9RdT2fKK8IKlLAFEPPcfNA671dGhKKWUpTRAl+OXrN8j64mPa7ywkwNd2jiHm4Ze4ufeVng5NKaWqhKUJQkSGAm9iW1HuQ2PMy+eUi718OJANjDXGbHClrVXmT/03+d99S7uUIhr4web+jej9t9e4VUclKaUuMJYlCBHxBiYDg4E0YJ2IzDTGbHWoNgxob7/1Bd4F+rrY1m2KCguZ+94TeM2YS+t9xbYL3C4N5eJ/vEOvtl2seEqllKr2rOxB9AGSjTG7AUTka2AE4PglPwL41NjWPV0tIkEi0gyIdKGtW+xP2Uri3SNpe8BwvJ5tTYYrJ06hT7PW7n4qpZSqUaxMEC2AfQ6P07D1Esqq08LFtgCIyHhgPEBERET5g2zdiTUNfUnoEc7QiR/SX6fEUEopwNoEIU62GRfruNLWttGYKcAUgNjYWKd1ynLjTwkVaaaUUrWalQkiDWjp8DgcSHexjp8LbZVSSlnIy8J9rwPai0hrEfEDRgEzz6kzE7hDbPoBWcaYAy62VUopZSHLehDGmEIReQCYj22o6lRjzBYRuc9e/h4wB9sQ12Rsw1zvKq2tVbEqpZQ6n9gGENUOsbGxJi4uztNhKKVUjSEi640xsc7KrDzEpJRSqgbTBKGUUsopTRBKKaWc0gShlFLKqVp1klpEMoG9FWweAhx2YzjuonGVj8ZVPhpX+dTGuFoZY0KdFdSqBFEZIhJX0pl8T9K4ykfjKh+Nq3wutLj0EJNSSimnNEEopZRyShPE76Z4OoASaFzlo3GVj8ZVPhdUXHoOQimllFPag1BKKeWUJgillFJOXVAJQkRuFpEtIlIsIiUOCRORoSKSJCLJIvKEw/ZgEVkoIjvtPxu5Ka4y9ysiUSIS73A7ISIP28ueFZH9DmXDqyoue709IrLZ/txx5W1vRVwi0lJEfhGRbfbf+QSHMre9XyV9VhzKRUTespdvEpGerratDBfiGmOPZ5OIrBSRbg5lTn+fVRjbQBHJcvj9PO1qW4vjeswhpkQRKRKRYHuZJe+ZiEwVkQwRSSyh3NrPlzHmgrkB0UAUsBSILaGON7ALaINt4aIEoJO97FXgCfv9J4BX3BRXufZrj/EgtgtcAJ4FHrXg/XIpLmAPEFLZ1+XOuIBmQE/7/UBgh8Pv0S3vV2mfFYc6w4G52FZJ7AescbWtxXFdBDSy3x92Jq7Sfp9VGNtAYFZF2loZ1zn1rwWWWP2eAZcCPYHEEsot/XxdUD0IY8w2Y0xSGdX6AMnGmN3GmHzga2CEvWwE8In9/ifA9W4Krbz7vQLYZYyp6FXjrqrs6/XY+2WMOWCM2WC/fxLYhm2tc3cq7bPiGOunxmY1ECQizVxsa1lcxpiVxphj9oersa3aWBUq87o9+p6dYzTwlZueu0TGmGXA0VKqWPr5uqAShItaAPscHqfx+xdLE2Nb8Q77zzA3PWd59zuK8z+cD9i7mFPddSinHHEZYIGIrBeR8RVob1VcAIhIJNADWOOw2R3vV2mflbLquNK2osq773HY/gs9o6TfZ1XG1l9EEkRkroh0LmdbK+NCROoCQ4EfHDZb+Z6VxtLPl5VrUnuEiCwCmjopmmiMmeHKLpxsq/RY4NLiKud+/IDrgCcdNr8LPI8tzueB14G7qzCuAcaYdBEJAxaKyHb7fz4V5sb3qz62P+SHjTEn7Jsr/H6du3sn2879rJRUx5LPWRnPeX5FkUHYEsTFDpvd/vssZ2wbsB0+PWU/P/QT0N7FtlbGdca1wApjjON/9la+Z6Wx9PNV6xKEMebKSu4iDWjp8DgcSLffPyQizYwxB+zduAx3xCUi5dnvMGCDMeaQw77P3heRD4BZVRmXMSbd/jNDRKZj694uw8Pvl4j4YksOXxhjfnTYd4Xfr3OU9lkpq46fC20rypW4EJEY4ENgmDHmyJntpfw+qyQ2h0SOMWaOiLwjIiGutLUyLgfn9eAtfs9KY+nnSw8xnW8d0F5EWtv/Wx8FzLSXzQTutN+/E3ClR+KK8uz3vGOf9i/JM24AnI54sCIuEaknIoFn7gNDHJ7fY++XiAjwEbDNGPPGOWXuer9K+6w4xnqHfbRJPyDLfljMlbYVVea+RSQC+BG43Rizw2F7ab/Pqoqtqf33h4j0wfY9dcSVtlbGZY+nIXAZDp+5KnjPSmPt58vdZ92r8w3bl0EakAccAubbtzcH5jjUG45t1MsubIemzmxvDCwGdtp/BrspLqf7dRJXXWx/KA3Paf8ZsBnYZP8QNKuquLCNkkiw37ZUl/cL2yETY39P4u234e5+v5x9VoD7gPvs9wWYbC/fjMPouZI+Z256j8qK60PgmMN7E1fW77MKY3vA/twJ2E6gX1Qd3jP747HA1+e0s+w9w/bP4AGgANt317iq/HzpVBtKKaWc0kNMSimlnNIEoZRSyilNEEoppZzSBKGUUsopTRBKKaWc0gShlIVE5HmxzSa7UUQ+9nQ8SpWHDnNVykIichjb2PQ9IhJgjMnxdExKuUp7EEpZ6x1glogsxHaFrVI1Rq2bi0mp6kJEWmCboz8GCAKWi8hOY8xWjwamlIu0B6GUda7HNutnsbHN/PkDtrU8lKoRNEEoZR0f/thLF5xPw6xUtaQJQinr/ApcLyJ17bN83kDVTAGtlFvoOQilLGKMiReRKdhWshPgfWNMvGejUsp1OsxVKaWUU3qISSmllFOaIJRSSjmlCUIppZRTmiCUUko5pQlCKaWUU5oglFJKOaUJQimllFP/D+xBEB5ykrBaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i_L in range(len(Lsub_list)):\n",
    "    plt.plot(δlist, LN[:,i_L]/np.log(2),label=r\"$ |A| = %d$\" % (Lsub_list[i_L]))\n",
    "plt.legend()\n",
    "plt.xlabel(\"δ\")\n",
    "plt.ylabel(\"log neg. [log2]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}