notebooks/sphere-ED.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## *This notebook demonstrates basic usage of the code*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load QHED code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "qhed_dir = \"../\"\n",
    "# CG_dir = joinpath(qhed_dir, \"CG\") # specify directory for precomputed Clebsch-Gordan coefficients\n",
    "\n",
    "include(joinpath(qhed_dir, \"src\", \"QHED.jl\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set matplotlib defaults"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using PyPlot\n",
    "\n",
    "PyPlot.matplotlib[:rc](\"font\", size=14)\n",
    "PyPlot.matplotlib[:rc](\"text\", usetex=true)\n",
    "PyPlot.matplotlib[:rc](\"xtick\", direction=\"in\", bottom=\"on\", top=\"on\")\n",
    "PyPlot.matplotlib[:rc](\"xtick.minor\", visible=\"on\")\n",
    "PyPlot.matplotlib[:rc](\"ytick\", direction=\"in\", left=\"on\", right=\"on\")\n",
    "PyPlot.matplotlib[:rc](\"ytick.minor\", visible=\"on\")\n",
    "\n",
    "# bblue = \"#0C2340\"\n",
    "# borange = \"#FC4C02\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set parameters and build setup\n",
    "$N$ = number of electrons\n",
    "\n",
    "$N_\\phi$ = number of flux quanta penetrating sphere\n",
    "\n",
    "For a LLL problem at filling factor $\\nu$, we will in general have $N_\\phi = \\nu^{-1}N - S$, where $S$ is the \"shift\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Building states list for N = 12, Nphi = 21, Lz = 0 ...\n",
      "Finished building states list of size 16660\n",
      "elapsed time: 0.313629542 seconds\n",
      "\n",
      "Building L^2 matrix of size 16660 x 16660 ...\n",
      "Working on LpLm, column j = 1: 1/16660 = 0.006002400960384154%\n",
      "Working on LpLm, column j = 833: 833/16660 = 5.0%\n",
      "Working on LpLm, column j = 1666: 1666/16660 = 10.0%\n",
      "Working on LpLm, column j = 2499: 2499/16660 = 15.0%\n",
      "Working on LpLm, column j = 3332: 3332/16660 = 20.0%\n",
      "Working on LpLm, column j = 4165: 4165/16660 = 25.0%\n",
      "Working on LpLm, column j = 4998: 4998/16660 = 30.0%\n",
      "Working on LpLm, column j = 5831: 5831/16660 = 35.0%\n",
      "Working on LpLm, column j = 6664: 6664/16660 = 40.0%\n",
      "Working on LpLm, column j = 7497: 7497/16660 = 45.0%\n",
      "Working on LpLm, column j = 8330: 8330/16660 = 50.0%\n",
      "Working on LpLm, column j = 9163: 9163/16660 = 55.0%\n",
      "Working on LpLm, column j = 9996: 9996/16660 = 60.0%\n",
      "Working on LpLm, column j = 10829: 10829/16660 = 65.0%\n",
      "Working on LpLm, column j = 11662: 11662/16660 = 70.0%\n",
      "Working on LpLm, column j = 12495: 12495/16660 = 75.0%\n",
      "Working on LpLm, column j = 13328: 13328/16660 = 80.0%\n",
      "Working on LpLm, column j = 14161: 14161/16660 = 85.0%\n",
      "Working on LpLm, column j = 14994: 14994/16660 = 90.0%\n",
      "Working on LpLm, column j = 15827: 15827/16660 = 95.0%\n",
      "Working on LpLm, column j = 16660: 16660/16660 = 100.0%\n",
      "Working on Lz^2, column j = 1: 1/16660 = 0.006002400960384154%\n",
      "Working on Lz^2, column j = 1666: 1666/16660 = 10.0%\n",
      "Working on Lz^2, column j = 3332: 3332/16660 = 20.0%\n",
      "Working on Lz^2, column j = 4998: 4998/16660 = 30.0%\n",
      "Working on Lz^2, column j = 6664: 6664/16660 = 40.0%\n",
      "Working on Lz^2, column j = 8330: 8330/16660 = 50.0%\n",
      "Working on Lz^2, column j = 9996: 9996/16660 = 60.0%\n",
      "Working on Lz^2, column j = 11662: 11662/16660 = 70.0%\n",
      "Working on Lz^2, column j = 13328: 13328/16660 = 80.0%\n",
      "Working on Lz^2, column j = 14994: 14994/16660 = 90.0%\n",
      "Working on Lz^2, column j = 16660: 16660/16660 = 100.0%\n",
      "elapsed time: 1.667097743 seconds\n",
      "\n"
     ]
    }
   ],
   "source": [
    "N = 12\n",
    "Nphi = 2N - 3\n",
    "\n",
    "setup = HaldaneSphereSetupLLL(N, Nphi, 0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate and plot Coulomb pseudopotentials for the lowest and first excited Landau levels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SphericalPseudopotentials(10.5, 9.5, [21.0, 20.0, 19.0, 18.0, 17.0, 16.0, 15.0, 14.0, 13.0, 12.0  …  9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0], [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0], [0.647293, 0.448627, 0.50082, 0.352973, 0.298615, 0.26634, 0.244242, 0.227975, 0.215456, 0.205535  …  0.185507, 0.180997, 0.177255, 0.174166, 0.171644, 0.169624, 0.168059, 0.166914, 0.166163, 0.165791])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_LLL = spherical_Coulomb_pseudopotentials(get_ell(setup), 0)\n",
    "pp_SLL = spherical_Coulomb_pseudopotentials(get_ell(setup), 1)\n",
    "\n",
    "# pp = deepcopy(pp_LLL)\n",
    "pp = deepcopy(pp_SLL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGyCAYAAAAYveVYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X9sm/dh5/GPpGiKvUqiJcNxfiiryTROuLRNSNNDi8Mtnck2WHA4IJXsBUH+ulqsK+xwyLXitOWuLS4dS7c73OXAGKKzv4L8kYgIDndLlxtpJAMOOZxlcSmW1XVqPl5Pzs9Z1EMqcKKpEu8PlaxpiRIlPY/48OH7BRgtnx9ff6Mn4vPJ92dHuVwuCwAAwCU6m10BAAAAKxFuAACAqxBuAACAqxBuAACAqxBuAACAqxBuAACAqxBuAACAqxBuAACAq9zS7Ao0w8rKit577z319vaqo6Oj2dUBAAANKJfLWlhY0B133KHOzvrtM20Zbt577z0NDQ01uxoAAGAbZmdnddddd9U935bhpre3V5Lk8/mUy+UsK7dUKmloaEizs7Pq6+uzrFxJCoVCmp6ednyZdpVrR5l2Pa92/7naUS6/W/aVy7Nq32cltd73YCAQUD6fr77H62nLcFPpiurq6rL8l0+S+vr6LC/Xjrra9c/fSnWVrH9e/Fz53WqlnyvPimcltdb3oKRNh5S09YDikydPNrsKDRsbG2uJMu0q16662oGfK8+rlX6uPCuelR3sqmuj7+2OdtwVvFQqqb+/X8Vi0dJkaVe5sAfPq3XwrFoHz6q1tNrzarS+bd1yEwqF5Pf7lUwmLSmvp6dH3/3ud9XT02NJebAXz6t18KxaB8+qtbTK80omk/L7/QqFQg1dT8tNCyRVAABAyw0AAGhTbTlbCgDgLktLS1peXm52NbBFXV1d6u7utrxcwo1FllfKOn+loI8WPtWB3lt19NCAujpZ/RgA7FQqlXTt2jUtLi42uyrYpp6eHu3fv9/SYSKEGwu89vb7eubVi7o6/0n12F379ujpR+/XIw/c3sSaAYB7lUolvfvuu/rMZz6j/fv3q7u7my11Wki5XNbS0pKKxaLeffddSbIs4BBudui1t9/XqRdzOnbfAT37+EM6fFuvLn24oOdev6xTL+Z05okAAQcAbHDt2jV95jOf0V133UWoaVF79uxRb2+vrl69qmvXrlkWbhhQvAPLK2U98+pFHbvvgFJPHlHg7n367Z5bFLh7n1JPHtGx+w7oBz+5qOWVtpuQBgC2Wlpa0uLiovr7+wk2La6jo0P9/f1aXFzU0tKSJWUSbnbg/JWCrs5/om995R513jS+prOzQ6cevkezhU90/kqhSTUEAHeqDB62YzAqdl/lOVo1KJxwswMfLXwqSTp82/obeB0+2FtzHQDAWrTauIPVz5FwswMHem+VJF36cGHd85c+WKi5DgAA2I9wswNHDw3orn179Nzrl7Vy07ialZWyzrxxWUMDe3T00ECTaggAQPthttQOdHV26OlH79epF3MafeGCTj18jw4f7NWlDxZ05o3LOvfzj3TmiQDr3QAA2oJhGEokEvL5fPJ4PJKk0dHRXa8H4WaHHnngdp15IqBnXr2or595s3p8aGAP08ABAG3DNE1FIhHNzMxUg000GlUqldr1gEO4scAjD9yuiP8gKxQDgMuw+nzj4vG4wuFwNdhIUiwWUzAYJNy0qq7ODn3JN9jsagAALMLq81uTTqcVi8Vqjnm9XpmmqVwup0AgsGt1YUAxAAA3qaw+f9/BXr3yrS/rH77/Nb3yrS/rvoO9OvViTq+9/X6zq+g4hmHI6/WuOe7xeJTNZne1LrTcAABwg5tXn68s0lpZfX70hQv6wU8uKuI/2LQuqnQ6LcMwNDc3p0QioVQqJUnK5/OSpEQisav1MU2z7rmBgQHNzc3tYm0INwAA1KisPv/s4w/VXX3+62fe1PkrhaYNRygUChofH68ufjcxMVEd67Jv3z6dOHFiw26gaDSqCxcuNPz3TUxMaHh4eMP6bGSj8GOHtg43oVBIXV1dGhsb09jYWLOrAwBwAKevPp/NZnX8+HEZhiFJOnHiRM0g3kaCxOTkpG31s0MymVQymWx4e4a2DjfT09OW7UAKAHCHG1efD9y9b835Zq8+Hw6HJa2GHK/XW9NCUwk8uzl490brBatCoVATvraj0ghRKpXU39+/6fUMKAYA4Aatsvp8JpOpBp2KdDrdlGAzMLD6s1ive8o0TQ0O7m73XVu33AAAcLNWWX0+m82uGTj80ksv6cSJE5vea/WYG4/HU532vZ6bQ5jdCDcAANzE6avPm6YpwzBqQkNlPZmpqSlJ2nBlYDvG3AwPD1dna1XkcjlJu99NRrgBAGAdTl59vjLe5saxLIZhVFtQUqnUrreWTExMKBgMyjTNar0mJyerYWs3EW4AAKjDqavPG4axppsoEAjoyJEjSqfT8nq96y6oZyePx6NMJqN4PC6fz6d8Pq9gMLhhd5ZdOsrlcnnzy9ylMtq6WCwyWwoAWtCnn36qK1eu6NChQ7r11ubMWoJ1Gn2ejb6/mS0FAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABchXADAABcpa3DTSgUkt/vVzKZbHZVAABwjVwup1gsZll5yWRSfr9foVCooevbOtxMT0/rZz/7mcbGxppdFQAAXCGXy+nYsWOWljk2Nqaf/exnmp6ebuj6Wyz92wEAQFsyDEOxWExer1cDAwNNrQvhBgCAm5mz0vW5+uf3Dkqeod2rTwvwer2ampqSJGWz2abWhXADAMCNzFkpeVRaul7/mu690th5Ao5DEW4AALjR9bnVYPPYWWn/vWvPX3tHeuXk6nVNCjfpdFqGYWhubk6JREKpVEqSlM/nJUmJRKIp9XIKwg0AAOvZf690x4PNrsW6CoWCxsfH1dHRIUmamJiQx+ORJO3bt08nTpxQIBCoe380GtWFCxca/vsmJiY0PDy8s0rvIsINAAAtJJvN6vjx4zIMQ5J04sSJarCRJNM0Ny1jcnLStvo5QVtPBQcAoNWEw2F5PB7lcjl5vd6aFppK4Nmo1aYdEG4AAGhBmUxG4XC45lg6nW77YCPRLQUAQEvKZrNrBg6/9NJLOnHixKb3MuYGAIB2dO2drR3fRaZpyjCMmpYb0zSVy+Wqa82kUimNjo6ue7/bx9wQbgAAuNHewdV1bF45Wf+a7r2r1zVJNpuV1+utGUhsGIY8Ho+8Xq9SqdSaLqvdZJpmQwOb7eKocGMYhhKJhHw+X/WB1UudN983OTmpwcFBzc3NaXBwUOPj43ZXFwDgRp6h1QX6HLxCsWEYa7qJAoGAjhw5onQ6La/XK6/Xu6t1Mk1T8XhchmHIMAylUikVCgWFQqFdfyd3lMvl8q7+jXWYpqlgMKiZmZlqsIlGowoGgxsGnMoP88Z+x2w2q6mpqbrNbqVSSf39/SoWi+rr67P2HwQAYLtPP/1UV65c0aFDh3Trrbc2uzrYoUafZ6Pvb8fMlorH49XpbRWxWGzTLdPj8bii0WjNsXA4vKWBUgAAwD0cE27S6bSCwWDNMa/XWx0gVY9hGOtu0NXsHUkBAEBzOCbcGIaxbv+gx+PZcHfRUCikaDSqdDpdPZbL5WpagAAAQPtwxIDijUZUDwwMaG6u/qCu8fFxvfTSSxoZGdHw8LCi0agymUx1KtxGSqVSzeeenh719PQ0XnEAAGCbxcVFLS4uVj/f/N6uxxEtN4VCYcPzm00nm5mZUTgcVjqdViQSaWgBI0kaGhpSf39/9U88Hm+4zgAAwF7xeLzmPT001NgMNUeEm52KxWIaGRmpzo4KBoPV7d83Mjs7q2KxWP0zMTFhd1UBAECDJiYmat7Ts7OzDd3niG6pivVaaAqFwobjZ2KxmHw+X3W6eDgc1sjIiKLRqI4cObLhHht9fX1MBQcAwKG2O1zEEeGmMrNpve4p0zQ1OFh/Fch0Oq18Pl/97PV6NTMzo2AwqJdeeokNxADAxRyyVBt2yOrn6Ihuqcpy0fXG1tRbQrreDCtJdDEBgIt1dXVJkpaWlppcE1ih8hwrz3WnHBFuJGl4eLimBUZSdX2beq0vXq9XhmGse66y5DMAwH26u7vV09OjYrFI602LK5fLKhaL6unpUXd3tyVlOqJbSlptaQkGgzJNszrGZnJysmZKt2EYikQimpycrLbmRKNRxWKxmu0XDMNoeDo4AKA17d+/X++++66uXr2q/v5+dXd3q6Ojo9nVQoPK5bKWlpZULBb18ccf684777SsbMeEG4/Ho0wmo3g8Lp/Pp3w+r2AwuGZjsJvH5YyPjyudTisajVZD0eDgIMEGAFyuMiHk2rVrevfdd5tcG2xXT0+P7rzzTksn+Dhm48zdxMaZAOAuS0tLWl5ebnY1sEVdXV1b6opq9P3tmJYbAAC2q7u727LxGmh9jhlQDAAAYAXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcBXCDQAAcJW2DjehUEh+v1/JZLLZVQEAAHUkk0n5/X6FQqGGrmfjTDbOBACgJTT6/m7rlhsAAOA+hBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqhBsAAOAqbR1uQqGQ/H6/kslks6sCAADqSCaT8vv9CoVCDV3fUS6XyzbXyXFKpZL6+/tVLBbV19fX7OoAAIAGNPr+buuWGwAA4D6EGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4CqEGwAA4Cq3NLsCLc+cla7P1T+/d1DyDO1efQAAaHOEm50wZ6XkUWnpev1ruvdKY+cJOAAA7BLCzU5cn1sNNo+dlfbfu/b8tXekV06uXke4AQBgVxBurLD/XumOB5tdCwAAoDYfUBwKheT3+5VMJptdFQAAUEcymZTf71coFGro+rZuuZmenlZfX1+zqwEAADYwNjamsbExlUol9ff3b3p9W7fcAAAA9yHcAAAAV2nrbinLXHtna8cBAIBtCDc7sXdwdR2bV07Wv6Z77+p1AABgVxBudsIztLpAHysUAwDgGISbnfIMEV4AAHAQBhQDAABXIdwAAABXIdwAAABXIdwAAABXcdSAYsMwlEgk5PP55PF4JEmjo6MN3zs5OanBwUHNzc0pFAppeHjYzuoCAAAHcky4MU1TkUhEMzMz1WATjUaVSqU2DTjZbFaJREJTU1PyeDwyDEORSEThcLhaFgAAaA+OCTfxeHxNGInFYgoGg5uGm5GREZ07d656r2EYKhQKttYXAAA4k2PG3KTTaQWDwZpjXq9Xpmkql8vVve/06dPyer0KBALVY+FwWPPz87TaAADQhhwTbgzDkNfrXXPc4/Eom83WvW9ycnLd+wAAQHtyRLeUaZp1zw0MDGhurv72BoZhaHh4WKlUqqas8fFxaysJAABagiPCzWbjY+qFn8pxwzAUjUarLTinT5/WyMiIpqamNiy3VCrVfO7p6VFPT0+j1QYAADZaXFzU4uJi9fPN7+16HNMttR03hqIbu6ZGR0eVTqc3HKsjSUNDQ+rv76/+icfjttUVAABsTTwer3lPDw01tpejI1puKtZroSkUCnUHBlcCTSgUqjleuT6bzdYMNL7Z7Oys+vr6qp9ptQEAwDkmJib01FNPVT+XSqWGAo4jws3AwICk9bunTNPU4ODgtsrN5/Mbnu/r66sJNwAAwDm2O1zEEd1SHo+nOu17PeFwuO69gUCg7oBjn89nSf0AAEDrcES4kaTh4eE1LS2VMTMbdS1Fo9E1Y2sMw6iWCQAA2otjws3ExISy2WxN683k5GTNjCfDMOTz+WrWvRkdHZVhGDUBJxaLaXR0lPVvAABoQ44YcyOtdk1lMhnF43H5fD7l83kFg8E1rS/rjcuZmZlRLBaTx+ORaZoKhUKscwMAQJvqKJfL5WZXYreVSiX19/erWCwyoBgAgBbR6PvbMd1SAAAAViDcAAAAVyHcAAAAVyHcAAAAVyHcAAAAVyHcAAAAVyHcAAAAVyHcAAAAVyHcAAAAVyHcAAAAVyHcAAAAV2nrcBMKheT3+5VMJptdFQAAUEcymZTf71coFGroejbOZONMAABaAhtnAgCAtkS4AQAArkK4AQAArkK4AQAArkK4AQAArkK4AQAArkK4AQAArkK4AQAArkK4AQAArkK4AQAArkK4AQAArkK4AQAArnJLsyuAXWTOStfn6p/fOyh5hnavPgAA2IBw0y7MWSl5VFq6Xv+a7r3S2HkCDgCgpdkSbs6ePatIJKLPfvazdhSP7bg+txpsHjsr7b937flr70ivnFy9jnADAGhhloSbs2fPKp1OKxAI6MSJEzp58qSef/55feMb37CieFhp/73SHQ82uxYAANjGkgHFpmnqhz/8obxer8bHxzU4OKhMJmNF0QAAAFtiScuN1+vVQw89pIceekgnT560oshdEQqF1NXVpbGxMY2NjTW7OgAAYB3JZFLJZFLLy8sNXW9ZuHnrrbf04IOt1d0xPT2tvr6+ZlcDAABsoNIIUSqV1N/fv+n1lnRLZbNZBQIBfe1rX9OPf/xjvfXWW1YUCwAAsGWWtNx4PB6trKzo7/7u75TNZjU+Pq75+XlNT09bUTysdO2drR0HAKDFWBJuBgYGJKk67uY73/mOFcXCSnsHV9exeWWDMVHde1evAwCghdUNN1sZQ+P1epn67XSeodUF+lihGADgch3lcrm83onPfe5z+sUvftFQId/85jeVzWY1Pz+vcDisEydOKBwOO3awbmVAUrFYdGwdAQBArUbf33UHFOfzef3pn/5pQ3+Zz+fT5cuXZRiGjh8/rr/5m7/RyMjI1msNAACwQ3Vbbo4cOaKJiQlduHBB0Wh0w60UisWizp07p8cee8yuelqqlVpullfKOn+loI8WPtWB3lt19NCAujo7ml0tAAB2XaPv77rhplgsVueSnz17VsViUd/+9rftqe0ua5Vw89rb7+uZVy/q6vwn1WN37dujpx+9X488cHsTawYAwO7bcbi5WbFYVDwe11e/+lX9wR/8gWUVbYZWCDevvf2+Tr2Y07H7DuhbX7lHh2/r1aUPF/Tc65d17ucf6cwTAQIOAKCtWB5uKs6dO6dsNrtpV5WTOT3cLK+U9fs/el33HexV6skj6ryhG2plpazRFy7o0ocLeuPbX6GLCgDQNnY8oLieY8eOKR6PK5PJ6Mc//vGOKon1nb9S0NX5T/Str9xTE2wkqbOzQ6cevkezhU90/kqhSTUEAMC5tr2I38mTJ1UsFvUnf/InruiqcpKPFj6VJB2+rXfd84cP9tZcBwAAfqNuy83zzz+/6c3z8/OKRCJKJBL6oz/6I5VKJUsr164O9N4qSbr04cK65y99sFBzHQAA+I26LTeTk5OKRCLK5/MyDKP6v5U/pmlWr60M28nn8+wnZYGjhwYU7P9Yf/XaX+vBR+9XZ8cNY27KZb36vy4q6OnS0UMDTawlAADOVHdAcWdnpzp+/VKtXOLxeOT1eqt/fD5f9f8fOnRo92q9Q04fUCxzVr/6byHdsvxJ3Ut+1bVHt/zxNNslAADaRqPv77otNx6PR4lEoiXDS8u7Pqdblj/RT4/+SP/1px36sLRYPXWwv0f/9gtlffH8d1b3iSLcAABQo264GR0d1cmTG+wgDdt98cGjOvvIF9euUPzBT6Xzza4dAADOVDfc/PCHP9zNeqCOrs4Ofck32OxqAADQMra8zo2bhEIh+f1+JZPJZlcFAADUkUwm5ff7FQqFGrp+2+vcuMH09LQzBxQDAICqsbExjY2NVQcUb6atW24AAID7tHXLjeNde2drxwEAAOHGkfYOSt17pVc2mK3WvXf1OgAAUINw40SeIWns/Oo6NvXsHWSNGwAA1kG4cSrPEOEFAIBtYEAxAABwFcINAABwFcINAABwFcINAABwFUcNKDYMQ4lEQj6fTx6PR9LqBp5bFYlElMlkrK4eAABoAY4JN6ZpKhKJaGZmphpsotGoUqnUlgJOKpVSNpu1q5oAAMDhHNMtFY/HFQ6Hq8FGkmKxmGKxWMNlmKapqakpO6oHAABahGPCTTqdVjAYrDnm9XplmqZyuVxDZaRSKUWjUTuqBwAAWoRjwo1hGPJ6vWuOezyehrqZ6t0PAADaiyPCjWmadc8NDAxobm6DbQh+LZ1Oa3h42MpqAQCAFuSIcFMoFDY8v1H4kaRsNkuwAQAAkhw0W2oncrmcwuHwlu8rlUo1n3t6etTT02NVtQAAwA4sLi5qcXGx+vnm93Y9jmi5qVivhaZQKNTMoLrZVqeK32hoaEj9/f3VP/F4fFvlAAAA68Xj8Zr39NBQYxtKO6LlZmBgQNL63VOmaWpwcHDd+yphaKPws5HZ2Vn19fVVP9NqAwCAc0xMTOipp56qfi6VSg0FHEeEG4/HU532vZ56XU7ZbFYzMzM1078Nw5C0ugCgx+NRIpGo+/f29fXVhBsAAOAc2x0u4ohwI0nDw8PK5/M1xyrr2wQCgbr33DyQuLJC8eTkpD0VBQAAjuaYMTcTExPKZrM1rTeTk5M1Kw4bhiGfz7fhujebzawCAADu5piWG4/Ho0wmo3g8Lp/Pp3w+r2AwuKZlpt608cqmm5XgE4lENDIysu3BxgAAoDV1lMvlcrMrsdtKpZL6+/tVLBYZcwMAQIto9P3tmG4pAAAAKxBuAACAqxBuAACAqxBuAACAqxBuAACAqzhmKjh2z/JKWeevFPTRwqc60Hurjh4aUFdnR7OrBQCAJQg3bea1t9/XM69e1NX5T6rH7tq3R08/er8eeeD27RVqzkrX5+qf3zsoeRrb7AwAgJ0i3LSR195+X6dezOnYfQf07OMP6fBtvbr04YKee/2yTr2Y05knAlsPOOaslDwqLV2vf033XmnsPAEHALArCDdtYnmlrGdevahj9x1Q6skj6vx1N1Tg7n1KPXlEoy9c0A9+clER/8GtdVFdn1sNNo+dlfbfu/b8tXekV06uXke4AQDsAsJNmzh/paCr85/o2ccfqgabis7ODp16+B59/cybOn+loC/5Brf+F+y/V7rjQYtqCwDA9jFbqk18tPCpJOnwbb3rnj98sLfmOgAAWhXhpk0c6L1VknTpw4V1z1/6YKHmOgAAWlVbh5tQKCS/369kMtnsqtju6KEB3bVvj557/bJWVmr3Sl1ZKevMG5c1NLBHRw8NNKmGAACsL5lMyu/3KxQKNXR9W4+5mZ6ebptdwbs6O/T0o/fr1Is5jb5wQacevkeHD/bq0gcLOvPGZZ37+Uc680SA9W4AAI4zNjamsbGx6q7gm2nrcNNuHnngdp15IqBnXr2or595s3p8aGDP9qaB3+jaO1s7DgCATQg3beaRB25XxH/QuhWK9w6urmPzysn613TvXb0OAIBdQLhpQ12dHdub7r0ez9DqAn2sUAwAcAjCDXbOM0R4AQA4RlvPlgIAAO5DuAEAAK5CuAEAAK5CuAEAAK5CuAEAAK5CuAEAAK5CuAEAAK5CuAEAAK7CIn6wzPJK2bptHQAA2CbCDSzx2tvv65lXL+rq/CfVY3ft26OnH71/extymrNs6QAA2BbCDXbstbff16kXczp23wE9+/hDOnxbry59uKDnXr+sUy/mtr7juDkrJY9KS9frX9O9d3VPKwIOAOAmhBvsyPJKWc+8elHH7jug1JNH1PnrbqjA3fuUevKIRl+4oB/85KIi/oONd1Fdn1sNNo+dlfbfu/b8tXdWdyG/Pke4AQCsQbjBjpy/UtDV+U/07OMPVYNNRWdnh049fI++fuZNnb9S2PpO5Pvvle540MLaAgDaQVvPlgqFQvL7/Uomk82uSsv6aOFTSdLh23rXPX/4YG/NdQAAbFUymZTf71coFGro+rZuuZmenlZfX1+zq9HSDvTeKkm69OGCAnfvW3P+0gcLNdcBALBVY2NjGhsbU6lUUn9//6bXt3XLDXbu6KEB3bVvj557/bJWVso151ZWyjrzxmUNDezR0UMDTaohAKDdEG6wI12dHXr60ft17ucfafSFC5r55bw+XvyVZn45r9EXLujczz/Sn/3h/ax3AwDYNW3dLQVrPPLA7TrzREDPvHpRXz/zZvX40MCerU8Dv9G1d7Z2HAAAEW5gkUceuF0R/0FrVijeO7i6js0rJ+tf07139ToAAG5CuIFlujo7tj7dez2eodUF+lihGACwDYQbOJNniPACANgWwg0czfLNONmzCgBcj3ADx7JlM072rAIA1yPcwJEs34xTYs8qAGgThBs4ji2bcd6IPasAwNVYxA+OU9mM81tfuafuZpyzhU90/kqhSTUEADgZ4QaOw2acAICdINzAcW7cjHM9bMYJANgI4QaOw2acAICdYEAxHKeyGeepF3MafeGCTj18jw4f7NWlDxZ05o3LOvfzj3TmicD217uxes8q1s4BAEch3MCRbNmM0449q1g7BwAch3ADx7J0M06pumfVm39/Sc//7yv6sLRYPXVbX4++8S8O6cufP7y1EMLaOQDgOG0dbkKhkLq6ujQ2NqaxsbFmVwfrsGwzzl977eotOvWTRR2770H9pyfuqVkc8ImffKQz+27RI55tFMzaOQBgm2QyqWQyqeXl5Yau7yiXy+XNL3OXUqmk/v5+FYtF9fX1Nbs62CXLK2X9/o9e130He2sWB5RWByqPvnBBlz5c0Bvf/krjrUPvvSWlfl8a/dv1w81m5wEADWv0/c1sKbQNFgcEgPZAuEHbYHFAAGgPbT3mBu3lxsUBA3fvW3PecYsDMsUcALaFcIO2cePigOuNudnR4oDX3tFyuax/eLekwvV/1sDe39Lv3tmnrrlfbK+yTDEHgG0j3KBt2LI44A1r53RJ+sJ612x17RyJKeYAsAOEG7QVyxcH9Azpjchf6cf//f/o6GcHNHJkSJ8d3Kt/nLuuqQuzOv+PBX37D7+kh7c9r01oAAARPklEQVQbQJhiDgBbRrhB27FyccDllbKefr2o+w7/np6+oavr/t+Rnn5wdXr5f3ijqDeOlLe/+CAAYEsIN2hLVi0OWJle/uzjD9WdXv71M2/q/JWCpYsRbhuDlAG0AcINsAO7Mb18eaW8tpVpOwUxSBlAmyDcADtg9/TyN/PXNP7C67o6/0n12F379uj0l1f05a0WxiBlAG2CcAPsgK3TyyW9/No5ffV3Pq+RSO1A5Zdf+3t9uXublbZjkDLdXQAcxFHhxjAMJRIJ+Xw+eTyruxeOjo5uel82m1Umk5FpmjIMQyMjIw3dB+yULdPLJS3vGdA/q0f/pfs56T1J/2P1+P2S/qMkdUufqEe/tWebXVRWorsLgMM4JtyYpqlIJKKZmZlqsIlGo0qlUhsGlVwup1wup0QiUS3n0KFDmpmZ0eTk5K7UHe3N8unlks4Xflv//tMf6S+Pe3X/wbXjeX72/oK+MWXoLwq/rS+t7Q3blGXjeCS6uwA4jmPCTTweVzgcrgYbSYrFYgoGgxuGm8nJyZoQ4/F4lEgkFI1GFYvF5PV6ba03IFk7vVxaHYD8nvbr7t/9ktSz9tf07sFf6b2p0rYGKls6judGVnd30dUFYJscE27S6bRisVjNMa/XK9M0lcvlFAgE1r3v5Zdfls/n0/j4ePXYkSNHJK12V9E9hd1i1fRyyd6Byn/+1z/XfYd/T88+/pAO39arSx8u6LnXL+vP//r/6q9+a2f1tgxdXQB2wDHhxjCMdVtZPB6Pstls3XAzMDCgubkN/usOaEF2DFReLpfVJelf3/mx/k34FnV2/FKakwK3SKnwLfrLjz+W/uk3122HZd1ddHUB2AFHhBvTNOue2yy85PP5NccMw5D0mxacekqlUs3nnp4e9fT0bHgPsBvsGKj81lyX7i/36OQ/xaWz8ZpznZJOSrpe7tHFuS4F79x6nW3p7mJmF9DWFhcXtbi4WP1883u7HkeEm0KhsOH5jcLPeiYnJxUOh+u29lQMDdV+gX33u9/V9773vS39XYBdrB6ofHVlUH+8+COdO/WA9nSvbU+5vrSs8Jm3FVsZVHAb9aW7C4DV4vG4vv/972/5PkeEGyul02kZhqGZmZlNr52dnVVfX1/1M602cBorByof6L1V72m/LnZ4Fbhj7Tiei7+c13v6YMvjeOzu7mqJmV20BgG2mJiY0FNPPVX9XCqV1jRMrMdR4Wa9FppCoVAzg2qz+2OxmDKZTEP39PX11YQbwImsGqhs14KDdnZ3tcTMLlqDANtsd7iII8LNwMDql+l63VOmaWpwsLEv9pGREWUyGaZ/A+uwa8FBO7u77OzqaonBz7QIAdviiHDj8Xiq077XEw6HNy0jGo0qkUjUBJuNppAD7ciOBQft6O6yu6urJQY/29UiRGBCG3BEuJGk4eHhNTOfcrmcJG0aUFKplEZGRmquMwxDhmEQboCbWL3goB3dXXbP7LKrRcjx44MITGgTjgk3ExMTCgaDMk2zOl5mcnJSU1NT1WsMw1AkEqnOhpJWF+qbmppSJBKphiFJymQy1S0ZANSycsFBO7q77OrqqrT0HP3sgJ6+IYgF7t6n1JNH9Mzzv5De216LkJ3jg5YPftG5XWh2jjkiNGGbHBNuPB6PMpmM4vG4fD6f8vm8gsGghoeHa667eVzOyMiITNNUNptdUyatNsDusLq7y66ZXf/wbklfkPTk5xbV+cFPa851SnrynkXpvV9ft8UWIbtagxzfhWbnDDRambBNjgk30up2Cxu1tni9Xs3Pz9ccu/kzgOawsrvLrpldH/xqr+4p9+jQ3/476W/Xnj+k1e6uD361V19osEw7W4OkFulCk6wfc0Qrk31BrA0CnqPCDYDWZlV3l10zu3pv8yq82MBu67c1PuPSrtagVuxCszwwSe3bymRny1UbtIgRbgA4kh0zu44eGlDnviH9xd/fqtSDX1zTIvSfMxfUNTC0pRYhO1qDpN+EppEjQzX1lKTOzg4NB4cc1YVm25gj2dPKZNk4Jsme0GRXEGu1FrFtItwAcCyrZ3bZ0SJkR2uQJBWu/7Mk6VD5qvTeW2vOe3W15rpG2Nka1Epjjlpi8LedZVpdrgM3uiXcAHA0K2d2Sda3CNnRGiRJvQO36Xq5R3v/5zfXPX+rVluEegdua7hMO7rQWnHMEUHMvoHqtgSxbSDcAGg7VrYI2TU+6MEHPq/jr/43fXFwWU8/er86O24ITeWynnn1on5a6NLLD3y+4TLt6EJrpTFHBDF7y7Wza3KrCDcA2pKVLUJ2jA/q6uzQyX/1L3XqxZz+X/ZXa0PTlX5HdKHZPebIytBk1zimVgpirdg1uR2EGwCwgNXjgyplOr0Lza4xR3aEJjvGMUn2hCa7gpgd5drdIrYdhBsAsIjV44Mk53eh2TbmyIbQZMc4Jsme0GRXELOjXDtn9m0X4QYAHM7JXWh2jTmyIzTZMY5Jsic02RXE7CjXriC2E4QbAGgzVneh2TXmyOrQZMc4Jsme0GRXELOjXLuC2E50lMvl8q79bQ5RKpXU39+vYrGovr6+ZlcHAFxh3cX2djDmSJJee/t9PfPqxZoZOEMDe/Rnf3j/tkKTnWWeejGnY/cdqBvEtlq2HWXaUe7ySlnHEy9vHpjGj+/434dG399tHW7uvfdedXV1aWxsTGNjY82uFgBgHXaEpnYOYnaUa1cQq0gmk0omk1peXtY777xDuFkPLTcAADu0ShCzo1y7gtiNaLnZAOEGAADr2RXEKhp9fzOgGAAAWMKO5RC2o7PZFQAAALAS4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALgK4QYAALhKW4ebUCgkv9+vZDLZ7KoAAIA6ksmk/H6/QqFQQ9d3lMvlss11cpxSqaT+/n4Vi0X19fU1uzoAAKABjb6/27rlBgAAuA/hBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuArhBgAAuEpbh5tQKCS/369kMtnsqgAAgDqSyaT8fr9CoVBD13eUy+WyzXVynFKppP7+fhWLRfX19TW7OgAAoAGNvr/buuUGAAC4D+EGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4CuEGAAC4yi3NrsCNDMNQIpGQz+eTx+ORJI2Ojtp2HwAAcB/H7C1lmqaCwaBmZmaqASUajSoYDG4YVLZzH3tLAQDQelpub6l4PK5wOFwNKJIUi8UUi8Vsuc8Oi4uL+t73vqfFxcVd/7uxdTyv1sGzah08q9bi1uflmJYbn8+nWCy2prWlo6NDMzMzCgQClt1nV8sNLUKthefVOnhWrYNn1Vpa7Xm1XMuNYRjyer1rjns8HmWzWcvvk6RUKrX1ijZJMplsiTLtKteuutqBnyvPq5V+rjwrnpUd7Kpro+9tR4Qb0zTrnhsYGNDc3Jyl91WcPXu2sQo6QCv9orRSXe3Az5Xn1Uo/V54Vz8oOdtW10fe2I2ZLFQqFDc/XCzHbva/SE7e0tKSrV69Wj/f09Kinp2fDMjdSKpVq/tdKy8vLlpdrR5l2lWtHmXY9r3b/udpRLr9b9pXLs2rfZyU5/3twcXGxZjzQ0tKSpN+8x+txRLjZbQsLC5KkK1euaGhoyPLy7ShTkvr7+1uiTLvKtauudjwvfq72lMvvVuvUlWfVWnVtpe9BafU9vlH5jgo367W0FAqFmplQVtx3xx13KJ/Pq7u7Wx0dHdXjO225AQAA1rm55aZcLmtpaUl33HHHhvc5ItwMDAxIWr+byTRNDQ4OWnpfZ2fnuoOQAQBA63PEgGKPxyOv11t3jEw4HLb0PgAA4F6OCDeSNDw8rHw+X3Msl8tJUt01bnZyHwBnyuVydRfhNAxD0WhUp0+fViqVaqnlHNxoo2cVDAaVzWZlmqZM01QqldLp06d3uYZoV47olpKkiYkJBYNBmaZZHSszOTmpqamp6jWGYSgSiWhycrLaKtPIfbuB/a1aRzAYVCKR0JEjRyRJL7/8skzT1Pj4eJNrhlwup2PHjq37u2OapiKRyJqtVlKpFL9rTbDRs6qcj0Qi1c/hcFiZTGa3qocbZLNZZTIZmaYpwzA0MjKy5rm57R3mmHDj8XiUyWQUj8fl8/mUz+cVDAY1PDxcc93N42savc9OfOm2Fr50nccwDMViMXm93upYupvV22pls/3nYK1GnpW0+mKs/IdnOBymJb1JcrmccrmcEomEpNX31aFDhzQzM6PJycnqMbe9wxyz/UIri8ViMk2z+i+KtPoFEAwGNT8/38SaYT2VjVX50nWmYDCocDhc/TKu2O4WLbBPvWclSadPn6Y11AGi0WjNu0laXeU3Go0qn8/L6/W68h3mmDE3rSydTisYDNYcqwx0roz/gXP4fD6Njo5qfHycF2IL2clWK0C7evnll9eMdap0yVd+b9z4DnNMt1Qr2+xLlxcosDM73WoFu29ubk6nT5+Wx+OpPj9acnZfI78fbnyHEW52iC/d1sOXbuvZ7lYraJ7KANWKaDSqWCy2bhcW7HPzbGJp9dlIqy04bn2H0S21Q3zpth7DMDQ+Pl7tmsrn83WnswLYnptnrI6MjOj06dN8JzpAZcZxIBBw7TuMcIO2w5du69ruFi1ovpvHeaA50um0DMPY9eVSdhvhxiJ86bYuvnSdb7tbraA5otFo3d+nSpcIdp9pmorFYspkMmveTW57hxFudogv3dbCl25rYquV1pJKpdb8PlW+I1txcKpbjIyMKJPJ1Awedus7jHCzQ3zptha+dFsXW620jsqYthtls1l5PJ5qSyl2VzQaVSKRqAk2uVzOte8wwo0F+NJtHXzpOl9lL6KbTUxMVPcqqmjGViv4jXrPKhKJKJ1O11yXSCR09uzZlu3maGWpVEojIyM17yPDMKr/oefGdxgrFFvANE0Fg8E1S1dHIpFd3QYCm6u8HCvPpfLsEokEz6qJTNNUPB6XYRjVl+Lw8LBCoVDNNH3DMDQ5OVndaqWyICN2T6PP6sb9jAqFgqLRaMu2ArSybDarRCJRs+WMJGUyGSUSCQUCAVe+wwg3FuFLt3XwpQugXezbt69ul9ONr3+3vcMINwAAwFUYcwMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFyFcAMAAFzllmZXAAC2I51OyzAMzc3NKZFIKJVKSZLy+bxM09Tk5KSy2axyuZwkaXp6WlNTU82sMoBdQrgB0JIKhYLGx8fV0dEhSZqYmJDH45Ek7du3Tz6fT+FwWOPj45KkSCSi06dPVz8DcC+6pQC0nGw2q+PHj8swDEnSiRMnqsFGkkzT1NzcnAKBQLOqCKCJOsrlcrnZlQCA7Uin04rFYsrn89VjhmHI5/Npfn6+JvB0dHQok8koHA43o6oAdhEtNwBa1nphJZfLyev11gSbyrgbgg3QHgg3AFpWNptVJBKpObZe4HnppZcINkAbIdwAaEmmacowjDWhZb3Ak06nFY1GJUmnT5/etToCaA7CDYCWlM1m13Q/1Qs8hmEoEAjIMAx5vd7driqAXcZUcAAtyTAMDQ8PrzkWDodrAo8kjY+PK51Oy+v1rrkHgPswWwoAALgK3VIAAMBVCDcAAMBVCDcAAMBVCDcAAMBVCDcAAMBVCDcAAMBVCDcAAMBVCDcAAMBVCDcAAMBVCDcAAMBVCDcAAMBVCDcAAMBV/j/IacmgAQNDcgAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x131e09c10>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots()\n",
    "ax[:plot](pp_LLL.m, pp_LLL.V, ls=\"none\", marker=\"o\", mfc=\"none\", label=\"\\$n=0\\$\")\n",
    "ax[:plot](pp_SLL.m, pp_SLL.V, ls=\"none\", marker=\"s\", mfc=\"none\", label=\"\\$n=1\\$\")\n",
    "ax[:set_xlabel](\"\\$m\\$\")\n",
    "ax[:set_ylabel](\"\\$V_m\\$\")\n",
    "ax[:legend](loc=\"best\")\n",
    "ax[:set_ylim](bottom=0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Building CG table for ell1 = 10.5, ell2 = 10.5:\n",
      "Working on alpha1 = 1: 1/22\n",
      "Working on alpha1 = 2: 2/22\n",
      "Working on alpha1 = 3: 3/22\n",
      "Working on alpha1 = 4: 4/22\n",
      "Working on alpha1 = 5: 5/22\n",
      "Working on alpha1 = 6: 6/22\n",
      "Working on alpha1 = 7: 7/22\n",
      "Working on alpha1 = 8: 8/22\n",
      "Working on alpha1 = 9: 9/22\n",
      "Working on alpha1 = 10: 10/22\n",
      "Working on alpha1 = 11: 11/22\n",
      "Working on alpha1 = 12: 12/22\n",
      "Working on alpha1 = 13: 13/22\n",
      "Working on alpha1 = 14: 14/22\n",
      "Working on alpha1 = 15: 15/22\n",
      "Working on alpha1 = 16: 16/22\n",
      "Working on alpha1 = 17: 17/22\n",
      "Working on alpha1 = 18: 18/22\n",
      "Working on alpha1 = 19: 19/22\n",
      "Working on alpha1 = 20: 20/22\n",
      "Working on alpha1 = 21: 21/22\n",
      "Working on alpha1 = 22: 22/22\n",
      "elapsed time: 96.792481813 seconds\n",
      "\n",
      "Building 2-body interaction matrix ...\n",
      "Working on alpha1 = 1: 1/22\n",
      "Working on alpha1 = 2: 2/22\n",
      "Working on alpha1 = 3: 3/22\n",
      "Working on alpha1 = 4: 4/22\n",
      "Working on alpha1 = 5: 5/22\n",
      "Working on alpha1 = 6: 6/22\n",
      "Working on alpha1 = 7: 7/22\n",
      "Working on alpha1 = 8: 8/22\n",
      "Working on alpha1 = 9: 9/22\n",
      "Working on alpha1 = 10: 10/22\n",
      "Working on alpha1 = 11: 11/22\n",
      "Working on alpha1 = 12: 12/22\n",
      "Working on alpha1 = 13: 13/22\n",
      "Working on alpha1 = 14: 14/22\n",
      "Working on alpha1 = 15: 15/22\n",
      "Working on alpha1 = 16: 16/22\n",
      "Working on alpha1 = 17: 17/22\n",
      "Working on alpha1 = 18: 18/22\n",
      "Working on alpha1 = 19: 19/22\n",
      "Working on alpha1 = 20: 20/22\n",
      "Working on alpha1 = 21: 21/22\n",
      "Working on alpha1 = 22: 22/22\n",
      "elapsed time: 5.209750639 seconds\n",
      "\n",
      "Building Hami matrix of size 16660 x 16660 ...\n",
      "Working on H, column j = 1: 1/16660 = 0.006002400960384154%\n",
      "Working on H, column j = 833: 833/16660 = 5.0%\n",
      "Working on H, column j = 1666: 1666/16660 = 10.0%\n",
      "Working on H, column j = 2499: 2499/16660 = 15.0%\n",
      "Working on H, column j = 3332: 3332/16660 = 20.0%\n",
      "Working on H, column j = 4165: 4165/16660 = 25.0%\n",
      "Working on H, column j = 4998: 4998/16660 = 30.0%\n",
      "Working on H, column j = 5831: 5831/16660 = 35.0%\n",
      "Working on H, column j = 6664: 6664/16660 = 40.0%\n",
      "Working on H, column j = 7497: 7497/16660 = 45.0%\n",
      "Working on H, column j = 8330: 8330/16660 = 50.0%\n",
      "Working on H, column j = 9163: 9163/16660 = 55.0%\n",
      "Working on H, column j = 9996: 9996/16660 = 60.0%\n",
      "Working on H, column j = 10829: 10829/16660 = 65.0%\n",
      "Working on H, column j = 11662: 11662/16660 = 70.0%\n",
      "Working on H, column j = 12495: 12495/16660 = 75.0%\n",
      "Working on H, column j = 13328: 13328/16660 = 80.0%\n",
      "Working on H, column j = 14161: 14161/16660 = 85.0%\n",
      "Working on H, column j = 14994: 14994/16660 = 90.0%\n",
      "Working on H, column j = 15827: 15827/16660 = 95.0%\n",
      "Working on H, column j = 16660: 16660/16660 = 100.0%\n",
      "elapsed time: 4.520362586 seconds\n",
      "\n",
      "106.749168 seconds (273.53 M allocations: 8.444 GiB, 1.05% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time Hami = HaldaneSphereHamiLLL(setup, pp, build_CG_table(get_ell(setup)));\n",
    "# @time Hami = HaldaneSphereHamiLLL(setup, pp, load_CG_table(get_ell(setup), CG_dir));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Diagonalize Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diagonalizing Hami of size 16660 x 16660 with eigs, nev = 50 ...\n",
      "elapsed time: 7.080410053 seconds\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# eig!(Hami);\n",
    "eigs!(Hami, 50);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Organize energy spectrum according to total angular momentum quantum number $L$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Organizing spectrum according to quantum number L ...\n",
      "elapsed time: 0.035451185 seconds\n",
      "\n"
     ]
    }
   ],
   "source": [
    "organize_spectrum!(Hami);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot energy spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHMCAYAAAA067dyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3b9vI/eB9/GPskEEBPBqRF4VnIBd0udCzcWcZfOUt2RcpMuS8T8QkWeoW+REs4ljN1xufE9HxOT+AyeRzlUpHM467QNIpF25MI6zxiPAnciRFgigPNnoKRacLCWSGv6e4bxfgAqS8x1+NSNyPvr+mo2rq6srAQAAhNgPVl0BAACAVSMQAQCA0CMQAQCA0CMQAQCA0CMQAQCA0CMQAQCA0CMQAfCldrutbDY7dptaraZsNivTNGWapp4+fbqk2q2PaY6hl3MDBA2BCFgj/QvbxsaGNjY25DjOyO22t7e1sbGheDyufD6/5JqO5jiOarWaTNMcWX9JKhQKikQiqtfrarVaqtfrKpVKisfjc6vLOhzPcSY9hl7PDRBEBCJgjdTrdT179kypVEqSVCqVxm6XyWTU6XRUrVaXWc2Rstms2/JgGMbI7drttiQpk8m4z8ViMdXrddm2PbdAEvTjOc6kx9DruQGCaoOVqoH18vTpU2UyGZmmKUnq9XpDt6vVanrw4IESicQyq+fZ9va2Hjx4oGazeeO1fD6vcrk89MK8vb0tx3E0r6+2dTme181yDMedGyCoaCEC1kyn01EsFtMvf/lLOY4jy7JGbheUi/d1Jycnun///tBum1gsJkmybXsu77Wux3OZxxAIAgIRsKYKhYIkqVwuD309yGNAIpGIHMdRt9td2nuu2/FcxTEE/OyHq64AgPlpt9tu104sFlMikZBlWXIcZ6BrxLZtd7sg6o9z6bdkvKnfqjHstUmt8/Fc1jEEgoIWImCNWJblDgCWpGKxKOnmYODr2wWNYRhDu6f6YeXg4GAu77POx3NZxxAICgIRsEb64136+jOIarXa2O28KBQK2t7enupn1LibeSsUCorFYm5wmdUij+dtVnW8530MgaCgywxYc7lcTrVabaAVY5rxLuVyeeT4GT94+vSp2u22Op3OQqeFT3s8G43GwBT326zieC/rGAJ+RAsRsCZGjWO5PhjYtu25Ll7oB7Ztq1Qqqdlszq2lZt7H0++DlxdxDIEgIRABa2LUOJbrg4Hb7XbgxrvcJpvNql6vz/X3CtvxXMQxBIKEQASsiVarNfI/+/6qw7VaTcfHx4FaL+c22WxWxWJx7hfyeR7PUbO5/GJRxxAIEgIRsCbGjfnI5XKSXs+Omna9HD8Oqi4UCkqn0zfG5tRqtZkXFZzH8azVaqrVaiqXy7JtW4VCwXO9lnW8F3kMgSBhUDWwBryMY+kPBp52/JDfBlU3Gg1Fo1E3nLyp2WzeeP762kHjzON4ptNpFQoFpVIp1Wo15XI5OY6jhw8fqtVq3VqHZRzvSY8hsM4IRMAaKJfL7o03R8nn86rVaoHpFhm3irJt29rb21MsFtPh4eGNcteZpunOnvLSdTXr8czn80qn0zdeMwzDHXe06m7LSY/h9df9PkgcmNgVgMDK5XJXsVjsStKVYRhXmUxm7PapVGpJNZvOwcHBVSqVujIM40rSlaSrWCx2lUqlrur1urtdIpFwXx/2c/33zGQyV7FY7KparY59/3kdzze/Wlut1lWr1Rp4rdfrjd3vMkx6DL2eGyCouNs9gNDod10tUn+6fq/Xu/GetVpN9Xqdu8QDPsSgagCh0G63lzLTKxaLuTdOvf7+9Xpd9Xp94XUAMDlaiACEQqFQWNqgcNu2Va1WlUwm1Ww23QUeGaQM+BeBCMDaazQaSiQSS18LqN1uyzAMX69BBOA1uswArL1MJrOSUHJyckIYAgKCQAQAAEKPLjMAABB6tBABAIDQIxABAIDQIxABAIDQIxABAIDQ4+auHv3973/X999/r7feeksbGxurrg4AAPDg6upKL1++1E9+8hP94Aej24EIRB59//332tnZWXU1AADAFE5PT/XP//zPI18nEHn01ltvSXp9QO/evTu3/V5cXGhnZ2fu+x0nmUzq+Ph4Ke+1qvdc5vuF4Ryu+98M5zDY77eK8yet9zFd9vst8hz2992/jo9CIPKo30129+7dhXzgFrXfYe7cubPUL41VvOcqfsd1Podh+JuROIdBfj9puedPWv9jum7n8LbhLgyqDqH9/f21f89V/I7LFIbjyTnk/fxu3Y9pGM7hm1ip2qOLiwttbW3pnXfe0Z07d7S/vz+XP5b+fs/Pz5eexDEfnMPg4xwGG+cv+BZxDiuViiqVil69eqVvv/321n3TZTah4+PjuX7gNjc39dFHH2lzc3Nu+8RycQ6Dj3MYbJy/4FvEOew3XPTD1m1oIfKI/0AAAAger9dv37YQ2batcrmseDwuwzAkSblcbi7lbNtWtVpVNBrV2dmZotGoDg4O5v9LAACAQPBlC5HjODJNU61Wyw01+XxepmmODUVeyjmOo1KppHK57JazLEv1el3VanXkvmkhAgAgeLxev30ZiAqFghzHGQgotm3LNE31er2ZyhUKBeXzecVisYGy/SA1CoEIAIDg8Xr99uW0+0ajIdM0B56LxWJyHEftdnumcrZty7KsG2Ujkcgcag4AAILIl4HItu0bLTiSZBjG0DAzSblkMql8Pq9Go+G+3m633S42AAAQPr4bVO04zsjXIpGIzs7OZip3cHCgw8NDZbNZZTIZ5fN5NZtN1et1T/W7uLgYeLy5uclUTwAAfOLy8lKXl5fu4+vX7VF810LU7XbHvj4q+ExSrtVqKZVKqdFoKJ1O6/333/dcv52dHW1tbbk/pVLJc1kAALBYpVJp4Drt9cbsvgtEy1AoFJTNZt3B16ZpqlareSp7enqq8/Nz96dYLC6yqgAAYALFYnHgOn16euqpnO+6zPqGtQR1u91bx/rcVq5QKCgej7vT8FOplLLZrPL5vB48eKBEIjF2/8u+eSAAQPrLX/82U/kf/8i3lzvM2bRDWXz3F9Kf7TWsC8xxHEWj0ZnKNRoNdTod97VYLKZWqyXTNHV4eHhrIAIAP1r3wLD7my9mKv/dk5/PqSZYV777BBiG4U6VHyaVSk1dbtQsNOl1E9vx8fF0lQaAFSMwALPxXSCSpEwmM9CKI8ldR2hcC46XcrZtDy3b7XaVTCanrjMAYHG++eS9VVcBa86XgahYLMo0TTmO4479qVarA1PjbdtWOp1WtVp1W428lMvn8yoUCgO37rBte6Kp9wDgN+seGPzepYfg8+VfmGEYajabKpVKisfj6nQ6Mk1TmUxmYLvr44W8lDs4OFCj0VA+n3dDUzQaJQwBCDQCAzAbX97LzI+4lxkAAMHj9frNvxQTSiaTunPnjvb397W/v7/q6gBzse4zlACET6VSUaVS0atXrzxtTwuRR7QQYZ3d+/CPM5VnhhIAvwr03e4BAACWiXZuzAVdLsG27jOUAOA2XIWWIAxhgUXhgi0If2MAsEh8Cy4BYQF+F4bQDgDj8C2GuaDLJdgI7cDq8Y/JanH0liAMYYEPIgDMhn9MVour2BIQFuB3YQjtCDZaT7Bo/IVMiIUZw2ndv4z9Xj8gDK0n/GMyXyzMuCAszBhuLFwIrBafQUyLW3cAwBvWvZVv3dF6gkXjEw54wJdx8IWhy2WdEUiDz+//lPAXBnjAlzEAzMbv/5TwLQ8gFGjlAzAOgQhAKNDKB6yW3/8p4RsCAAAsnN//KfnBqisQNMlkUru7u6pUKquuCgAAGKFSqWh3d1fJZNLT9qxD5BHrEAEAEDxer9+0EAEAgNDzd4ceAMATv6/xAvgdnwAAXEzXgN/XeMHt+ByuFkcPABdTwAf4HK4WgQgA1oDf13gB/I5ZZh4xywzrjKZ6YPX4HC4Gd7sH4BlfpMDq8TlcLabdT4iFGQEA8D8WZlwQuswAAJjeqroE6TID5oi+fQCYjd9n0fn2W9q2bZXLZcXjcRmGIUnK5XIzl6vVanIcR6lUyn39TbFYbE6/AdaJ3z/IAIDZ+DIQOY6jdDqtVqvlhpZ8Pq9arTY2FHkp12q1VKvVRu6j1+sNDUoAAGB6fl8awpdjiAqFghzHUbVadZ+zbVumaarX681ULpvNKp1OKxKJDJQ9Pj5WMplUJpMZum/GEIUbXWYAEEyBHkPUaDRUKBQGnovFYnIcR+12W4lEYupyyWRyaCvT8fHxyDAEEGgAYL35ctq9bdtDx/IYhiHLsmYqNywMFQoFFYvFGWoMwO/+8te/zfQDYL357t9ex3FGvhaJRHR2djZTuevjg9rttpLJpOdxQxcXFwOPNzc3tbm56aksgNVhYDwQDpeXl7q8vHQfX79uj+K7FqJutzv29VHBZ9py1Wp1oq6ynZ0dbW1tuT+lUslzWQAAsFilUmngOr2zs+OpnO9aiJbJsqyJZ5Sdnp4ODMqidQgIBr/PcAEwH8ViUY8fP3YfX1xceApFvg1Ew1p0ut3urQFmknLlcvnGIOzb3L17l1lmWDthmEUXhDoCmN20Q1l89w3Rnw4/rAvMcRxFo9G5lbMsS+VyeZbqAmuB8TUAws53gcgwDHeq/DCpVGou5fqzzliZGgD8LwytmOvO7+fQl38hmUxGnU5n4Ll2uy1JI9cgmrRc/3lWpQYYXwP/oxUz+Px+Dn0ZiIrFokzTlOM4bmCpVquq1+vuNrZtK51Oq1qtuq0/Xsr1jZq+D4RRGP579vt/pwBWy5efcMMw1Gw2VSqVFI/H1el0ZJrmjenx18cLeS0nSdFolO4yIET8/t/prNY98NGKGXx+P4e+vJeZH3EvMyDY7n34x5nK+z0QrfvvB0wr0Pcy87NkMqk7d+5of39f+/v7q64OAI/8/t8pgPmqVCqqVCp69eqVp+1pIfKIFqJwW/fuCAQfLUTAcLQQAXO07uNPACDsCEQAsAboEgRmQyACPOBiA7+jWxaYDZ8gwAMuNgCw3n6w6goAAACsGoEIAACEHoEIAACEHgMjJsTCjEAwsZYUEC4szLggLMwYblxMg4+FC4FwYmFGYI5YmBEA1huBCEAosJYUgHEIRIAHXEyDj27LYKPbGovGXwjgAV+mwGrRbY1FY9o9AAAIPf7tBYA1sO5dSnRbY9H8/QnwIdYhAuBH696l5PfABv9hHaIFYR0iAH7GOkvAcKxDBAAhQpcSMBsCEQCsAbqUgNnwCQIQCus+6BjAbPiEAwiFdR90DGA2rEMEAABCjxYiAKHAoGMA4xCIAIQCY4AAjEOX2YSSyaR2d3dVqVRWXRUAADBCpVLR7u6uksmkp+1ZmNEjFmYEACwSMyEXg4UZAQAIEGZCrpZvA5Ft2yqXy4rH4zIMQ5KUy+XmVs62bVWrVUWjUZ2dnSmZTCqTycz3lwAAAIHgyy4zx3FkmqZarZYbavL5vEzTHBuKvJazLEvlcln1el2GYci2baXT6YFy19FlBgBYJLrMFsPr9duXgahQKMhxHFWrVfc527ZlmqZ6vd7M5ba3t/X8+XMlEglJrwNSNpvVixcvCEQAAKwRr9dvX84yazQaMk1z4LlYLCbHcdRut2cq9/TpU8ViMTcMSVIqlVKv1xsZhgAAwHrzZSCybVuxWOzG84ZhyLKsmcpVq9Wh2wAAgPDyXYej4zgjX4tEIjo7O5upnG3bymQyqtVqA+UODg481e/i4mLg8ebmpjY3Nz2VBQAAi3V5eanLy0v38fXr9ii+ayHqdrtjXx8VfLyU65e1bVupVEq5XM4NQtls1lP9dnZ2tLW15f6USiVP5QAAwOKVSqWB6/TOzo6ncr5rIVqkN0PTm91muVxOhUJB7XZ7YGzRMKenpwODsmgdeo3ZEQAAPygWi3r8+LH7+OLiwlMo8u1VaFhLULfbvXXg87hy/RB0fRnv/j4ty7o1EN29e5dZZkOwoBgAwA+mHcriuy6zSCQiaXgXmOM4ikajcy33pk6nM0lVAQDAmvBdC1G/JWfUWKFUKjVTuUQiMXJgdjwen6LGkKRvPnlv1VUAAGBqvgtEkpTJZG601vTXERrXpeWlXD6fV71eH9jGtm23PKbDGCAAQJD5rstMej0gyrKsgdaearU6EGRs21Y8Hh9Yl8hLuVwuJ9u2BxZ4LBQKyuVyrE8EAEBI+fLfesMw1Gw2VSqVFI/H1el0ZJrmjRac6+OFvJZrtVoqFAoyDEOO4yiZTHpehwgAAKwfX97LzI+4lxkAAMET6HuZ+VkymdTu7q4qlcqqqwIAAEaoVCra3d29sdTOKLQQeUQLEQAAwUMLEQAAgEe+HFQNAJgMt88BZsMnAADWALfPAWZDIAIA+B4tYFg0/kIAYA2s++1zaAHDohGIAGAN0AICzIZPEADA99a9BQyrx7T7CbEwIwAA/sfCjAvCwowAsDr3PvzjTOUZQxReLMwIAADgEWOIAAC+xxgiLBqBCPCANVCA1eIzhEXjLwzwgDVQAGC9MYYIAACEHi1EgAeMXwCA9UYgmlAymdSdO3e0v7+v/f39VVcHS8L4BQAIlkqlokqlolevXnnannWIPGIdIgAAgod1iAAAADwiEAEAgNAjEAEAgNAjEAEAgNBj6gwAVuIGEHp8iwFgJW4AoUeXGQAACD1aiCbEwoxYR6zEDWDdsDDjgrAwIwAAwcPCjAAAAB75tsvMtm2Vy2XF43EZhiFJyuVycylnmqbK5bIePHggSTo6OpLjODo4OJjzbwEAAILAl11mjuPINE21Wi031OTzeZmmOTYUeS23sbExUC6VSqnZbI6tE11mAAAEj9frty8DUaFQkOM4qlar7nO2bcs0TfV6vZnL9UOS4zhKpVJKJBK31olABAQbay0B4eT1+u3LT3ij0VChUBh4LhaLyXEctdvtkQHGa7l4PO6p+w3A+lj3tZYIfMBsPH0CksmkIpHIzG92dXWl7777Tt9+++3Y7WzbViwWu/G8YRiyLGtkIJq2HAAE3boHPmDRPAWi7e1tffHFbB+2vv5A5lEcxxn5WiQS0dnZ2czlzs7O9PTpUxmG4ZbzOqD64uJi4PHm5qY2Nzc9lQWwOqy1BITD5eWlLi8v3cfXr9ujeApE2Wx2uloN8f777499vdvtjn19VPCZpFx/JlpfPp9XoVAYeG6UnZ2dgccfffSRfvvb395aDsBqrXuXEIEPeK1UKunjjz+euJynb4i9vb2JdzzKf/zHf8xtX9Oq1+sDj7PZrNLptIrFojs7bZTT09OBQVm0DgHwg3UPfIBXxWJRjx8/dh9fXFzcaMwYZu6foO+++07tdlvdbleRSESJREL37t2beD/DWoK63e6tgWWacv1uPMuylMlkxu7/7t27zDIDAMCnph3KMrdA9OLFC2WzWZ2fn2tra0vSP8LJxsaGms2mp2DUH7w9rAvMcRxFo9GZyuXzeWWzWaVSqRvb2bZ9a/0AAMD6mVsg+vzzz/X8+XM3DL3JcRyVy2WVSqVb92MYhjtVfphhQWaScrVaTaZpDrzWD1HMQgMAIJzmdi+z+/fvDw1D0uuwctvssjdlMhl1Op2B59rttqTxocVLuYODgxtrEFmWNXEdAQDA+phbILqtu+nk5MTzvorFoizLGmjtqVarA4OhbdtWPB6XZVkTlUun02o0Gu7jfuvVs2fPbh2fhPD6y1//NtMPAMDf5tZllslk9Pbbb2tjY0OGYSgSiajb7brh5PrMrnEMw1Cz2VSpVFI8Hlen05FpmjcGPF8fL+SlXCqVkmVZ7m0+ut2uqtXqyK44QGLROwBYd3O/l9lXX32lk5MTOY7jdkO9++6783yLleBeZuF278M/zlSeQAQAqzH3m7teXFzo8PBQ5+fnisVi+sUvfjG3ygZB/4C+8847unPnjvb397W/v7/qamFJuE8UAARLpVJRpVLRq1ev9O23384vEH366af69a9/LUk6Pz/X0dGRNjY29Ktf/Wro9l9++aUePHiwNq0ptBABABA8Xq/fngdVP3r0SB988IFevnypra0t7e3tKZvN6ne/+52+/PLLG9v/27/9m6dp9gCwDAyMBzDORGOIbNvWL3/5SyWTSf37v/+7/vVf/1XS6zWIer3ejdair776ShsbG/rpT38631qvAC1EWGdh6BJkHBgQTl6v3xN9i8ViMZ2cnOjzzz/Xo0eP1Ov1lEqllE6nJb2e9v5mq1Cv15uy+gCWiVl0AMJupllm7XZbR0dHajQa7jpE29vbisVikl5PcV+XbjNaiLDOwtB6EoZWMAA3zX2WmRcvXrxQr9fT0dGRnjx5Mq/d+gKBCOuMsABgXa0kEL3p66+/XouxQ30EIgAAgmfus8wmtU5hCAAArLeZAtGnn3468Pji4mLg8ddffz3L7n0pmUxqd3dXlUpl1VUBAAAjVCoV7e7uKplMetp+pkB0dXU1EIKuD6Ce5IauQXF8fKxvvvmGVaoBAPCx/f19ffPNNzo+Pva0/UyBKJFI6OjoyH18fZr99vb2LLsHAABYipkC0cOHD3V0dKQ///nPkqSNjY2B11mHCAAABMHMc2WfPHmiVCqleDwuSfrzn/+sq6srNZtNRaPRmSsIAACwaDPPMkskErJtW/fu3VOn09HDhw+VzWYVjUbdm8ECAAD42VxWUzMMQ/V6XZJ0fn6ura2teewWAABgKTy1EF2fTj/ObWFokn0BAAAsg6dAdHh4qC+//HLmN/vyyy8HZqUFEesQAQDgf5OuQ+T51h3Pnj2TbdvK5/O6d+/eRJX67rvv9OTJE7399tuBHVfErTsAAAger9dvz2OI9vb29OLFC+VyOW1sbCidTiuVSikWi914g4uLC52cnKjdbutPf/qTNjY29Nlnn+n+/fvT/0YBxo0zAQDwt6lu7np+fq6joyPV63XZti3btgfWIDIMQw8ePFA2m1U2m12LQdaztBDd+/CPM733d09+PlN5AADCau4tRG/a2trS3t6e9vb23OfOz8/d1wAAAIJkbn0xBKHRvvnkvVVXAQAAjHFrIHrx4oU2NjYmHkiNfwjDGCDGSQEAgszTVejJkydqtVpKpVJTzTLD+tv9zRczlWecFABglW4NRPfv39dnn30mSfrqq68IRwAAYO1MNctMkp4/f656vR6acNQfpf7OO+/ozp072t/f1/7+/qqr5Rt0mQEA/KRSqahSqejVq1f69ttvb51lNnUgetMiwpFt2yqXy4rH4zIMQ5KUy+UWUi6dTqvZbI7dhoUZAQAIHq/X77kEojf1w9GLFy+USCSmCkeO48g0TbVaLTfU5PN5maY5NtxMU65Wqymfz+u2w0AgAgAgeFYWiN40bctRoVCQ4ziqVqvuc7ZtyzRN9Xq9uZVzHEfZbFaWZRGIAABYQ16v355u7jqthw8f6rPPPtPx8bFSqZSePHmi999//9ZyjUZDpmkOPBeLxeQ4jtrt9tzK9VuHAABAuC1tJOvDhw/18OFDT9vatq1YLHbjecMwZFmWEonEzOVGbQsAAMJnoS1E03AcZ+RrkUhEZ2dncynXaDSUyWSmqyQAAFgrngPR119/rYuLC/fx+fn5wON56Xa7Y18fFXwmKWdZFmEIAAC4bg1EL168UCQSUSKR0Pb2tv7zP//Tfa1arerOnTsLreAitNvtqbvLLi4uBn4uLy/nXDsAi/CXv/5tph8AwXB5eXnjWu3FrWOICoWCnj17pkePHsm2bX344Yf69NNP9etf/1p7e3sqFAozV36YYS1B3W7XnU4/bblareZpPaNRdnZ2Bh5/9NFH+u1vfzv1/gAsB7eXAcKhVCrp448/nrjcrS1EyWRSjx49kvR6xtbR0ZG2trb0hz/8QRsbG9rY2Ji8tmNEIhFJw7vAHMdRNBqdulw/LN0WqsY5PT3V+fm5+1MsFqfeFwAAmK9isThwnT49PfVUbqpZZnt7e3r+/LmOjo6mKT6WYRjuVPlhUqnU1OUsy1Kr1RqYam/btqTXCzgahqFyuTy2fnfv3mUdIiCAvvnkvVVXAcASbG5uanNzc+JytwaiVCqlDz/8UL/73e/U6XTchRUfPnyoq6urWxc0nEYmk1Gn0xl4rr+O0Kgp917KJRKJG4Opa7WaLMsaWMwRwPrhfnkAxrm1y+zdd99VsVjUycnJjVWmU6nUjQAyD8ViUZZlDbT2VKtV1et197Ft24rH47Isa6Jy142brg8AAMLB079MW1tbevfdd4e+dv/+/blWSHrd/dVsNlUqlRSPx9XpdGSa5o3WnevjhbyWk/5xE9h+oEqn08pmszMNuAYAAME0073M+vcq+9nPfqZf/OIX7vPPnj1TMpnUT3/607lU0g+4lxkAAMGzlHuZWZalzz77TFdXV/r666/d5/f29nR8fLyQhRsBAADmbaZA1F/c8NGjRzo5ORl4bW9vT7VabZbdAwAALMVM0y7eHJDc6/VuvL61tTXL7n0pmUzqzp072t/f1/7+/qqrAwChMOtq4cwyDJ9KpaJKpaJXr1552n6mv5BUKqUPPvhAv//977W9vX3j9Xkv2ugHx8fHjCECgCVjpXFMqt9w0R9DdJuZAtG7776r+/fvKxqN6sGDB4pGo+6ss5OTE7VaLf3qV7+a5S0AAAAWbqZZZn3tdlu5XM5dBHFjY0O5XE6///3vZ66gXzDLDABWhy4zTMvr9XsufyGJRMIdVH1+fr6WY4cAAKtDoMGieZpl9sEHH3je4W1haJJ9AQAALIOnyN2/Aeo8XJ+eDwQBzfUAsN4837rjvfdmv1P01dXVQm71ASwaM1wAYL15CkRHR0eLrgcAAMDK0I4/IRZmDKdvPpm9hRQAsDyTLsw4l2n3YcC0ewAAgmcpN3cFAABYBwQiAAAQegQiAAAQegQiAAAQegQiAAAQep6n3X/66aeKxWKKxWIyDEP37t1bYLUAAACWx3ML0cHBgQ4PD9XpdNRqtfTdd98NvP7ixQt9/fXX866f7ySTSe3u7qpSqay6KgAAYIRKpaLd3V0lk0lP23teh+jtt9/W//zP/4zd5vnz52o2m/qnf/on5XK5tVqvh3WIsM64VxuAdeX1+u35WywWi926zcOHD/Xw4UM9f/5c9+7d089+9jP913/9l9e3ALAi3KsNQNh57jIzDGPg8bNnz/T3W7qXAAATwElEQVSHP/zhRteZ9DoYHR0dqV6vz1xBAACARfPcQrSxsTHweG9vT59//rni8bhSqZTS6fRAN1kqleLO9kBAcK82AGE3U8f/o0eP9O677+qLL4Y3tycSiVl2D2BJGAMEIOw8d5nZtq2XL1/eeP7Bgwcjy0Sj0elqBQAAsESeA1Gr1ZJhGPqXf/kXffDBB/rv//5vXVxc3OhKe5PHCWwAAAAr5bmdPJFIKJ/P609/+pMODw9VrVYHwpBpmkqlUgMLNo4LSwAAAH7hORClUint7e1pb29P0uuFGFutlizLUrPZHAhIqVRKiURCJycni6n1CiWTSd25c0f7+/va399fdXUAAMAQlUpFlUpFr1698rS954UZb3N+fq7j42M1m009f/5c7XZbGxsbnividyzMCABA8Mx9YcbbbG1tKZVKKZVKSZIcxxk74Po2tm2rXC4rHo+7ayDlcrm5lOu3akWjUXU6HRmGoXK5PHVdAQBAsC1srq1hGG44mpTjOEqn0+5AbknK5/Oq1WpjQ5GXcpZluaGpLx6Py7ZtFpIEACCk5tZlNk+FQkGO46harbrP2bYt0zTV6/VmKpfP52VZljqdjrtNPzSNOxR0mQEAEDxer9+ep90vU6PRkGmaA8/FYjE5jqN2uz1TuXQ6zYKRAABggC8DkW3bQ28maxiGLMuaqVwmk7nRNWZZljKZzIy1BgAAQeW7QOQ4zsjXIpGIzs7O5lquVqtJen2zWgAAEE6+u4FRt9sd+/qo4DNpOcuyVK/XdXJyomaz6Q7Cvs3FxcXA483NTW1ubnoqCwAAFuvy8lKXl5fu4+vX7VF810K0LKlUStVqVc+fP1c6nXZbim6zs7Ojra0t96dUKi24pgAAwKtSqTRwnd7Z2fFUznctRH3DWoK63e6tLTmTljMMQ4VCQfl8Xg8ePLh1wPXp6enAKHVahwAA8I9isajHjx+7jy8uLjyFIt8FokgkIml4F5jjOIpGozOVa7fbisViAwGpv4Dk4eHhrYHo7t27TLsPob/89W8zlf/xj3z3UQOAtTTtUBbffUsbhuFOlR9m1GKPXso5jiPTNJXL5QbWKgJus/ubL2Yq/92Tn8+pJgCARfBdIJJeT41/c+FESe46QuNacLyUi8ViyufzA9vYti3p9RpFAAAgfHwZiIrFokzTlOM4btdWtVodWD/Itm2l02lVq1W31chLuethSHo9AOvN+7AB133zyXurrgIAYIF8GYgMw1Cz2VSpVFI8Hlen05FpmjcWT7w+XshLuYODAzUaDVWrVRmGIdu29f777+vg4GApvxuCiTFAALDefHkvMz/iXmYAAARPoO9lBgAAsEwEogklk0nt7u6qUqmsuioAAGCESqWi3d1dJZNJT9vTZeYRXWYAAAQPXWYAAAAeEYgAAEDoEYgAAEDoEYgAAEDoEYgAAEDoEYgAAEDoEYgmxDpEAAD4H+sQLQjrEAEAEDysQwQAAOARgQgAAIQegQgAAITeD1ddAQBYhr/89W8zlf/xj/i6BNYZn3AAobD7my9mKv/dk5/PqSYA/IguMwAAEHq0EAEIhW8+eW/VVQDgYwSiCSWTSd25c0f7+/va399fdXUAeMQYICBcKpWKKpWKXr165Wl7Fmb0iIUZAQAIHhZmBAAA8IhABAAAQo9ABAAAQo9ABAAAQo9ABAAAQo9ABAAAQo9ABAAAQo9ANKFkMqnd3V1VKpVVVwUAAIxQqVS0u7urZDLpaXsWZvSIhRkBAAger9dv365lb9u2yuWy4vG4DMOQJOVyubmUsyxLzWZTjuPItm1ls1lP+wYAAOvJl4HIcRyl02m1Wi031OTzedVqtbHBxUu5drutdrutcrnslrl//75arZaq1eqCfzMAAOBHvuwyKxQKchxnIKDYti3TNNXr9WYql8/nbwSfWq2mfD6vTqejWCw2dN90mQEAEDyBvpdZo9GQaZoDz8ViMTmOo3a7PVO5o6MjPX36dGCbBw8eSHrdlQYAAMLHl4HItu2hLTWGYYwNLV7KRSIRnZ2dza+yAAAg8Hw3hshxnJGvjQszXst1Op0br9u2LekfLUUAACBcfBeIut3u2NdHBZ9py0lStVpVKpVSIpG4tX4XFxcDjzc3N7W5uXlrOQAAsHiXl5e6vLx0H1+/bo/iyy6zZWo0GrJtW/V63dP2Ozs72tracn9KpdKCawgAALwqlUoD1+mdnR1P5XzXQtQ3rEWn2+260+nnUc5xHBUKBTWbzVv323d6ejowSp3WIQAA/KNYLOrx48fu44uLC0+hyHeBKBKJSBreBeY4jqLR6NzKZbNZNZvNkVPth7l79y7T7gEA8Klph7L4rsvMMAx3qvwwqVRqLuXy+bzK5fJAGBo3pR8AAKwv3wUiScpkMjdmg/XDyriBz17L1Wo1ZbPZgeds23ZnmwEAgHDxXZeZ9Lr/zzRNOY7jju2pVqsDA59t21Y6nXZniHktZ1mW6vW60un0QItQs9l0b+cBAADCxZeByDAMNZtNlUolxeNxdTodmaapTCYzsN318UJeymWzWTmOM3SBRy/T7gEAwPrx5b3M/Ih7mQEAEDyBvpcZAADAMhGIJpRMJrW7u6tKpbLqqgAAgBEqlYp2d3eVTCY9bU+XmUd0mQEAEDx0mQEAAHhEIAIAAKFHIAIAAKFHIAIAAKFHIAIAAKFHIAIAAKFHIJoQ6xABAOB/rEO0IKxDBABA8LAOEQAAgEcEIgAAEHoEIgAAEHoEIgAAEHoEIgAAEHoEIgAAEHoEIgAAEHoEogmxMCMAAP7HwowLwsKMAAAEDwszAgAAeEQgAgAAoUcgAgAAoUcgAgAAoUcgAgAAoUcgAgAAoUcgAgAAoUcgmhALMwIA4H9rszCjbdsql8uKx+MyDEOSlMvl5lqu3W7r8PBQ5XL51v2yMCMAAMHj9fr9wyXWyTPHcZROp9VqtdxQk8/nVavVxoaiScq12209fPjQU8gCAADrzZddZqVSSalUyg01klQoFFQoFGYuZ9u2stmsDg8PFYlE5l95AAAQOL4MRI1GQ6ZpDjwXi8XkOI7a7fZM5WKxmOr1usrl8kBwAgAA4eXLQGTbtmKx2I3nDcOQZVlzLwcAAMLNd4HIcZyRr0UiEZ2dnc21HAAAgO8CUbfbHfv6qOAzbTkAAABfzjLzs4uLi4HHm5ub2tzcXFFtAADAmy4vL3V5eek+vn7dHsV3LUR9w1p0ut3urQOhpy3n1c7Ojra2ttyfUqk0l/0CAIDZlUqlgev0zs6Op3K+ayHqT4Uf1gXmOI6i0ehcy03q9PR0YGEnWocAAPCPYrGox48fu48vLi48hSLfBSLDMNyp8sOkUqm5lpvU3bt3WakaAACfmnYoiy+7zDKZjDqdzsBz/XWEEonE3MsBAIBw82UgKhaLsixroLWnWq2qXq+7j23bVjweH1hfyEu5NzmOw+wzAADgvy4z6XX3V7PZVKlUUjweV6fTkWmaymQyA9tdHy/kpZzjOCqVSrJtW7Ztq1arqdvtKplM6uDgYCm/HwAA8Bff3u3eb7jbPQAAweP1+u3LLjMAAIBlIhABAIDQIxBNKJlMand3V5VKZdVVAQAAI1QqFe3u7iqZTHranjFEHjGGCACA4GEMEQAAgEcEIgAAEHoEIgAAEHoEIgAAEHoEIgAAEHoEIgAAEHoEogmxDhEAAP7HOkQLwjpEAAAED+sQAQAAeEQgAgAAoUcgAgAAoUcgAgAAoUcgAgAAoUcgAgAAoUcgAgAAoUcgmhALMwIA4H8szLggLMwIAEDwsDAjAACARwQiAAAQegQiAAAQegQiAAAQegQiAAAQegQiAAAQegQiAAAQegSiCbEwIwAA/sfCjAvCwowAAASP1+v3D5dYp4nYtq1yuax4PC7DMCRJuVxuLuWm3TcAAFhPvmwhchxHpmmq1Wq5gSWfz8s0zbHBxUu5afdNCxEAAMHj9frty0BUKBTkOI6q1ar7nG3bMk1TvV5vpnLT7ptABABA8AT6XmaNRkOmaQ48F4vF5DiO2u32TOWm3TcAAFhfvgxEtm0rFovdeN4wDFmWNVO5afcNAADWl+8GVTuOM/K1SCSis7OzqctNu+9Z/eWvf5up/I9/5LvTBADAWvHdlbbb7Y59fVSo8VJu2n2/6eLiYuDx5uamNjc3x5bZ/c0Xt+53nO+e/Hym8lisy8tLlUolFYvFW/8W4E+cw2Dj/AXfPM/h5eWlLi8v3cfXr9uj+LLLzM92dna0tbXl/pRKpVVXCSt2eXmpjz/+eOADiGDhHAYb5y/45nkOS6XSwHV6Z2fHUznfBqJhrTXdbtedKj9LuWn3LUmnp6c6Pz93f4rF4q1lvvnkvZE//+fX/0v/938/0v/59f8auc28rWKV7WW/57qvJB6G48k55P38bt2PaVDPYbFYHLhOn56eeirnu0AUiUQkDe8CcxxH0Wh06nLT7vtNd+/eHfjx0rT34x/9cMzPHV39v0v9+Ed3Rm4zb2G4uAX1g+xVGI4n55D387t1P6ZBPYebm5s3rtVe+G4MkWEY7jT4YVKp1NTlpt23JPWXa/LaF+lVf3/z3u84r169Wur7reI9l/l+YTiH6/43wzkM9vut4vxJ631Ml/1+izyH/X3etuyi7wKRJGUyGXU6nYHn+msEJRKJmcpNu++XL19Kkue+yEktar+jbG1tLfX9VvGey36/dT+HYfib4RwG+/2Wff6k9T+m63QOX758Ofb38eVK1aNur5FOp5XJZCS9Xk8onU6rWq26LTteynnZZpi///3v+v777/XWW29pY2NjYb87AACYn6urK718+VI/+clP9IMfjB4p5MtAJL0OPNVqVfF4XJ1OR/F4fOBeY/3bbdTr9YGurtvKed0GAACEh28DEQAAwLL4bpYZAKxKOp1edRUArIgvB1WHhW3bKpfLisfj7ngmuu6Cw7IsNZtNOY4j27aVzWY5fwFWq9W4n2FA9YdBRKNRnZ2dKZlMjh0TCv+4fu6i0agODg5WU5krrESv17uKxWJXvV7PfS6Xy11Vq9UV1gpetVqtq3K57D7u9XpXhmFc5XK5FdYK0+r1elepVOqKr8TgaTabV6lUyv0u7XQ6N75b4U+9Xu/q4OBg4Llms7my71HGEK1IoVCQ4ziqVqvuc/2B4r1eb4U1gxf5fH7g3EmvWxjy+bw6nY5isdiKaoZpPH36VLFYTNls9ta1SuAv29vbev78ubtsimVZymazevHihae7D2B1CoWC8vn8je/L/kzwZWMM0Yo0Gg2ZpjnwXH/RyP66SPCvo6MjPX36dOC5Bw8eSBLdLgFj2zYBNqD6QfbNNeRSqZR6vR5hKABs2x76fdm/q8SyEYhWZNSXsGEYXFADIBKJ6OzsbNXVwBw0Gg3GmwRUtVolzAZYMplUPp9Xo9Fwn2u32ysLswyqXoFRtw6RuNAGxfXVzqXXIVf6R0sR/M+yLMJQgNm2rUwmo1qtJukf360rG5SLiRwcHOjw8FDZbFaZTEb5fF7NZlP1en0l9aGFaAWG3Vz2TeMCE/yrv2r6uFvAwF/a7TYtDAHV/560bVupVEq5XM4NQtlsdpVVwwRarZZSqZQajYbS6bTef//9ldWFQATMQaPRkG3bK/vPBpOr1WoskxBgb/5j+WaozeVyajQajMUMiEKhoGw2605SMU3TbfFbNgLRCg1rCep2uwwGDBjHcVQoFNRsNjl3AdH/7HG+gqsfgpLJ5MDz/XPKWEz/KxQK7q2zcrmcOp2OEomE8vn8SgItY4hWoD+CfljXmeM4ikajy64SZpDNZtVsNul6CRDLstRqtZTP593n+mPA8vm8DMNQuVxeVfUwB8PG+cFfGo3GwHmKxWJqtVoyTVOHh4dLH35AIFoBwzDcKfbDvHmzWvhbPp9XuVweCEPtdptxRD6XyWRuDKbur1R9fX0p+FcikRg5CSUejy+5NpjEuOUuisWijo+Pl1wjusxWJpPJ3PgPpt9EyMU0GGq1mrLZ7MD5sm3bbWlAsDCZIXiGda30P3/MHvS3WCw28ruy2+3e6ApdBlaqXhHHcdzVOPt93vl8Xul0mg9yAFiWpXK5fONmoM1mU+VymVAbIP17ClqW5c5Y4r50wRGPx1Wv193PXDabVSQSoaUvAJ4+faqzs7OB7mnbtlUoFFYyQYVAtEL9m9rF43F1Oh13cBn8b3t7e2SLAh8pYHn6kxoMw5DjOIrH46xDFCCNRmNgQsoqb+5KIAIAAKHHGCIAABB6BCIAABB6BCIAABB6BCIAABB6BCIAABB6BCIAABB6BCIAABB6BCIAABB6BCIAoVUoFGSapjY2NhSPx5XNZtVoNFZdLQArwErVAEKtVqspn8+r0+mMvPs2gPVHCxGAUGs2m4rFYoQhIOQIRABCzbIs907pAMKLQAQgtBzHkeM4SqfTq64KgBUjEAEILcuyJEmpVGrFNQGwagQiAKHVbDZlGAbjhwAQiACEl2VZtA4BkEQgAhBSjuPItm0lk8lVVwWADxCIAITSycmJpNHjh2q12jKrA2DFCEQAQqnZbErSyCn3nU5nmdUBsGIEIgChNG78UK1WoysNCBkCEYDQcRxH7XZ7aOtQrVZToVBQJpNZQc0ArMoPV10BAFimfD7vrj9kWZby+bwkybZtnZycyHEc5XK5VVYRwApwc1cAABB6dJkBAIDQIxABAIDQIxABAIDQIxABAIDQIxABAIDQIxABAIDQIxABAIDQIxABAIDQIxABAIDQIxABAIDQIxABAIDQIxABAIDQ+/9rHniv00u7VAAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x13e6b6650>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots()\n",
    "ax[:plot](Hami.spectrum[:, 1], Hami.spectrum[:, 2] - minimum(Hami.spectrum[:, 2]), ls=\"none\", marker=\"_\", ms=\"12\", mew=\"1.5\")\n",
    "ax[:set_xlabel](\"\\$L\\$\")\n",
    "ax[:set_ylabel](\"\\$E~[e^2/(\\\\epsilon\\\\ell_0]\\$\")\n",
    "ax[:set_title](\"\\$N=$(N), N_\\\\phi=$(Nphi)\\$\")\n",
    "ax[:set_ylim](bottom=-0.001);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define ground state wave function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi = Hami.eigenstates[:, 1];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate and plot entanglement spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating entanglement spectrum for NA = 6, LzA = [18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 31.0, 32.0, 33.0, 34.0, 35.0, 36.0, 37.0, 38.0, 39.0, 40.0, 41.0, 42.0, 43.0, 44.0, 45.0, 46.0, 47.0, 48.0] ...\n",
      "elapsed time: 0.530126771 seconds\n",
      "\n"
     ]
    }
   ],
   "source": [
    "subsystemA = collect(Int(get_Norb(setup)/2+1):get_Norb(setup))\n",
    "NA = Int(N/2)\n",
    "LzAvec = collect(sum(setup.states.possible_mz[subsystemA][1:NA]):sum(setup.states.possible_mz[subsystemA][end-NA+1:end]))\n",
    "\n",
    "ent_spec = entanglement_spectrum(psi, setup.states, subsystemA, NA, LzAvec);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAHUCAYAAAAgOcJbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3c9vI+l95/GPRnZkGB6JrZ7LTtLYniJgIDosArIJ7CmHiLQPAfaQkPZf0OSBt9lAnA4G7p5gMEoJey0ErP4LuknvbQ4J2ZdcFkg3iRwWPjhgwUAj3osllThGYMWRaw8KaVEiJVJkseqper+ABppUVenRU1XP863nV20EQRAIAADAQB9EnQAAAID7IpABAADGIpABAADGIpABAMSC7/vyfT/qZMAwBDIAgMh1u13t7+9rf38/6qTAMAQyAIDIFYtF2bYddTJgIAIZAABgLAIZJE6/31elUrl1G9d1ValUlM/nlc/ndXR0tKbUJcd98nCec4P5jPJy9K/RaESdpFRa5Jrm+g8HgQwkaVwhbWxsaGNjY+aAu0qlogcPHmhjY0PZbFa1Wm3NKZ3N9325rqt8Pn/rgMFGo6Hd3V21Wi31ej21Wi0dHh4qm82uLC1JyM/bLJqH856b+0p6fl/nuq6ePn0q27bVarXUarVUKBRWEsysIi9931ej0VC73Va73Vaj0Rgfp9vtKp/Pq1Qqqd1uq1ar3TgPruuq2+3Kdd3YBmiLXNNhX/+pFwD/qdfrBcViMZAUHBwczNyu1WoF5XJ5jSm7W7lcDorFYtBsNoNMJhMUi8Wp2/V6val/W6fTCSQF1Wp1ZWkyOT9vs2gezntuVpGuJOb3db1eL8hkMsHp6enE95ZlBblcbmW/Y5m8zOVyQa/XG38eDAaBZVkT+1mWFQwGg6DX6wW2bQdB8Ptr6Orf1mq1Qrtm7muRa3pd13+aEchgzLbtYDAYBJlMJshkMjO3azabE4VU3NxWWFSr1RsVwNX9VhnbJyU/r1smD8MsyJOa39dZljU14J71/X0sk5ejIOW6XC4XNJvNIAguA5Zpx+10OlP3lRTbc7bINU0gEw66ljA2GAxkWZZ+9KMfyfd9dbvdmdvlcrk1p2413r17p08++WRq865lWZIkz/NW8ruSmp/rzMNFJDW/r+p2u/I8b+o4i8FgoGazuZLfs0xevn37VplM5sa2u7u7GgwG48+ja2VeUVxTMAOBDG4Y9UnPmgppch/v7u6ufN/XycnJ2n5n0vIzijxcRNLy+6pWqyXpcqryOtwnLwuFwtSgw/O8pcahLRr4ID0IZCDpcjR9Pp+XdFlg5HI5dbvdGwWV53nj7Uw0Gpw6rVAcFb6rKDCTnJ9h5mGtVrtXZZfk/L7qasvI0dGRjo6O1Gg0Vjogdtm8LJfL2t3dVb/fn9hWkqrV6p2/3/O8id/luq6KxaKxrWgIH4EMJF0WkFef8p49eyZJOjw8vHU702QymakF4qigPjg4WMnvSXJ+hpmHo66TRbsRkpzfV3mep0wmI9d1Va1WdXBwINu2lc1mlc1mV9LatIq87PV6evXqlVzXleu6ajab6vV64/1s21a/31ej0bhxrpvNprrdrtrttlzX1WAwUKfTWfrvQoJFPUgH8TBtkKCkGwPybpvBMMvBwcF40OCi/zqdzsK/7z4D6nK5XGBZ1sxBrIsKMz/vsu78HpknD+86N6NZLIuKMr+vHz/MvB/9TaNBs1dZlrWSvy8ueXlfcS5vGOwbDgIZBEEwvfCqVquBpImbe5XTk8OyaGFh23YgKRgMBitLwyrzs9VqrSxdYZk3D8MqyMO4fqvV6lKBXRgk3ZiePFIul2f+bBFJKgvWgUAmenQtYWZf9/WBfssO1osjz/N0eHioTqezssGEq87PuA6qHQkjDxf9/au+fkeLscV1YPC0WUGrmDGW5rIA5iKQwcy+7usD/fr9vtHjC6apVCpqtVor/bvSlp9h5OEiwsjv0ZiMuE35zWQyU4MYSXr48KGky+nx95W2axfJQCCDmTNQJI2XDnddV2/fvk3UzIFKpaJnz56tvEBeZX56nhfraadh5eEiVn39Hh0dqVarybKsiXVP4uDJkyczW4mOj48lLTfrLqyyoFarzVyLBlgWgQxmPuFJv58ueXh4eO9m9kajoQcPHtzrX1iFX6PRUKlUUrlcnvjedd2ln8JXkZ+j2R62bcvzvKmzO2ZZV36HmYeLWOX1e3X6uGVZC3frhZ33pVJJ0vT1W0bfLRPIhFEWrLubLo7lDUIW9SAdRGswGEydAXHVaKDf6H0ocXfXgLpWqzXzb5n23phFBk+uIj+LxeJ4UOXoWKenpyt7j84qLJqHI6uetbTq6/fqANZyuRyrPA+Cy79X1wbdjuRyuanpnff6DassODg4MKr8WBSDfaNHi0zK2bZ95xPcqEnZlD7x21ad9TxPT58+1atXr5TP5yf+ZbPZiUW8JCmfz+vBgwdztzAsm5+1Wk2lUunGzzKZzHhsQtQWzcOr7loRuFQqLfSG4FVev67rTiz9b1lW7MbIWJalarV6Y6Xdfr+vfr+vly9fTny/yPUbRlkQ5266VVlkles4r4httKgjKUSjWq0GlmWN14e4623AcX+KODg4CIrF4vilhZICy7KCYrE4MX05l8uNfz7t3/W/s1wuB5ZlzfWkuor8vHpL9nq9idYJrWBq7SosmofznpsguMzHedbzWPX1O2rxOjg4GP8bvf05jorFYlCtVoNerxc0m83AsqypeTbP9RtWWTAYDMatMMVi0eg3jl+3yDW9yLa4n3jepUDM3BXIrMLobcPTfmez2Yx9MLkqvV5v7eu3THujd6vVik3wOM0oiGm1WnemcR3X73Vx76ZDcnxrHa0+gMn6/f5aZg5ZljV+IePVQZf9fl+tVmv8wsCkOzk5WetMLdd1VSqVbgx0vbouSxxn6+VyubnSta7r96pp3XQMpEVYGCMD3OHVq1drGx/U6XR0eHiodrutXq8n13X17t07dTqdW2eUJMm6Kl7P81SpVFSr1W68y6fdbo8XgWs0GkZXwuu8fqXLcSDNZlOdTmf8Qst+vx/bxQVhvo0gCIKoEwHEVbvdVi6XW/sTbb/fVyaTifUaMmHwfV+u667s5Z1pF8X1W6vVZNv2RODdbrdVqVR0enqamoAc60OLDHCLcrkcSTDx7t271AUxkvT69WuCmBVa9/U7TzcdsGoEMgBiY7ToGsySlm46xBNdSwAAwFixa5Hp9/vj6P26brerRqMxXjTMdd01pw4AAMRJrKZf9/t97e/vT21eHq1cOVrR0vd9ffLJJ+r1emo2m+tOKgAAiIFYtMiM+ldfvXql3d3dqds0m82JQYCZTEa2ba/9BXUAACA+YhHIWJalVqt1Y8reVa9fv9bR0dHEd0+ePJEkBo8BAJBSsQhk5rG7u6vj4+OokwEAAGIkVmNkbjPtzamjLqVRywwAAEgXY1pkpmk2myoWi7F8DwoAAAifMS0y17XbbXmep16vd+t2v/vd7/SLX/xC3/72t7WxsTH+fmtrS1tbW2EnEwAAzOH8/Fzn5+fjz0EQ6Le//a0eP36sDz6Y3e5iZCDj+74ajcZcL9L75S9/qWw2u6aUAQCAVXr//r3+6I/+aObPjQxkKpWKOp3OXO8Q+fDDDyVJP/vZz8b/l1bXIlMoFPT27duljzPNcDjUo0eP9P79e21vb6/8+GGm3fTjm3xeJXPzZtnj/9u//8etP//TP/1T/eM//uPMn3/3D+5fJHJezT3+bdfNN998o+9///v6+c9/PlGHXHXXdRPmdRnlNS+t9rxeb5H55ptvtLe3NzPfR4wLZEZvVr0axPT7/ZnjZEbdSX/4h38YSuGyubkZWqE1sr29bWTa43z8u27+D/7gO/rWd7478+fL3vxSeOdVCjfv43xe/9tnX9++wf+w9d//1/+Z+eNf/O2f3+v3XsV5Ne/4d103//V//m+Vmv935s/vum7CvC6jvubDPK/D4VCSJoaFTGNUIOO6riqVykTQ4nmePM+LbMBvvV6P5PeuQthpj/Px937y97dv8Odf3brNKiq8MIWZ93E+r0nHeY3u+JguDvkeu5dGZrNZFYvFG68d6Ha7sm1bpVJp4vtOpyPbtmcGMsPhUDs7Ozo7Owu95WTVTE573D2+6ynmDssEMpzX8NzV0naXZbuWOK9muu26GQ6/0ccf/xf98pf/T9vb4XQt3WWZrqVljh21ee+pWPwFvu/r8PBw3Lriuq5OTk5UKBTGryWoVCryfX/qKr5JnX69tbWl58+fp3Z2VZg36M/+5odLHfsut6X9YuNb+uuffKGLjW/N3C7OhUucRZlvab9fTXbbdbP54Xf1k7/+TA8+/K627nl9hXldUlbEsEVm1XhKMleUrSbLMjntmM3kp1+T0450MqpFBgBMcOfYqjtEGaCanHbgNgQyuLewn/DC7P4xOe1ILlpNgMXRtYR7M7n7xOS0h43KdLaw8ybM65LzCtPQtQTgXuiCmM3kytzktAO34cpGKtH1gzgy+bqkxQdR4cpBKpleaFJpJJPJ54WWPETF3LsmIUxeKMnkp0fTUWkgbUwuKxEucj9iYVZIYVd23LzJRICaTKY/2JhcViJcqamJCoWCNjc3Va/XY/FuCGAZYVYaBKjJxIMNTOE4jhzH0cXFxVzbM/06YjSXzmZ6+jEd5zUapi85QFmZPky/NgTv4JiN5t5kCvO8UiHNZnqXIWUlZuHsAUgMgt/ZqKyRVFzZiC3TnyAxHec1GrRWIam4MhFbFJzJFOZ5JUiajdYqJBU1BVKJp9PZTM6bJJ8XANNx1yOVeDqdzeS8MTkICxutVUiq5N61AFLH5CAsbEkO0pBuXNlIJZ5OZ3v3+X7USQCAuRHIIJV4Op3tyZdvlto/ylYNAlQgfSjNASQGASqQPtz1ACbQqgHAJAQyACbQqgHAJJRYgIGYZpxMnFdgcam56guFgjY3N1Wv11Wv16NODrAUphknE+cVkBzHkeM4uri4mGv71AQyb9++vfU14AAAIHqjBofhcKidnZ07t98IgiBYQ7oiM8qIs7MzAhkkBl0QyWTyeTU57YineetvrhzAQBT60Qi7sjb5vNIthqiYe9cAwJpRWQPxQyADGIhmfCwq7GuG9YcQFUozwEC0DEQj7Mo6zGAj7GuG4BhR4coDMIHWntnC/tsIUIHFJbfEARIszJYBKtNkousHSUUgA4SA2S24jzCDDa4ZJBVXNmLL5C4Ok1s1eHKPDsEGsDjuGsSWycGAyahMAZiEEgsIAa0awGqZ3EKLcHFmEVsmBwMUmsBq0UKLWShtEVsEA8nEkzWAVaJEALBWPFlHw/QA0uQWWoQrNYFMoVDQ5ubm+PXgAJAmpgeQUQdSWB/HceQ4ji4uLubafiMIgiDkNEVq3teAA6tk+tNvmMibaDz+7Oul9o86kEH6zFt/UyIAITD96RfJQ9cMkopABsBaEeRFg5YsJBVXNhACnn4BYD0IZBBbJo+lMPnpN+x8f/f5/lLHB4CrzC1tkXh0QcwWZrARdr4/+fJNqMcHkC4EMoCBCPIA4FLsApl+v69Xr17Jtu0bP/M8T7ZtK5vNKpPJSJKq1eq6k4g1YZxJNMLOd85rMpncFQyzxerK6ff72t/fnxqc+L6vUqmkXq83DmJqtZpc1yWYSSgKttnCDAbCzvcwj09lGh1aCRGVWNy1nuep0WjIsizt7u5O3ebw8FDFYnEcxEhSo9FQPp8nkMHCTK/wov79cUVlCqRPLEpDy7LUarUkSd1ud+o27XZbjUbjxn6+76vf7yuXy4WeTiRH2BWe6YESsCi6DBEVY0pLz/NkWdaN7zOZjLrdLoEMYoWWgWhQmUaH4BtRMeLK831/5s92d3d1fHy8xtQgCajwkonKFEgfI+76k5OTW39+W6AzMhwOJz5vbW1pa2trqXTBXGFXeARKs9HtBmCa8/NznZ+fjz9fr7dnSU2J8OjRo4nPz58/14sXL6JJDBKPynY2ut0ATHN4eKgvvvhi4f2MKm2ntbycnJxMzGSa5f379xOvAac1Znk8WUeHvAeQNM+ePdOnn346/jwcDm80QkxjRGk2mpI9rYvJ9309fPjwzmNsb29PBDJYHk/W0TE57+l2AzDNfYd8GBHIZDKZ8VTraYrF4ppTBOC+aA0CsErGlCjlclmDwWDiu36/L0lMvY4IT9bRIe+jQZceED+xu6t835/a8vLs2TPl83n5vj8eE9NsNscL6WH9KJSjQ95Hw+QuPSCpYlEa+r6vw8NDeZ4nz/Pkuq5OTk5UKBR0cHAg6bJ7qdPp6PDwUNlsVoPBQPl8XuVyOeLUAzfx5B4N8h1In40gCIKoExGm4XConZ0dnZ2dMdgXa/P4s6+X2p8n9/sJO98JlID1mbf+5q4CgDkRiADxQ4sMEAKe3KNBvgPJQYsMECEqxGiQ70D6fBB1AgAAAO6LxxcAmNOvfv2bpfb/6HvfWVFKAIwQyADAnJ58+Wap/ZmNBqweXUsAAMBYtMgAwJzefb4f6vGZdQUsLjVXfaFQ0Obmpur1uur1etTJAWCgsMe48AoEQHIcR47j6OLiYq7tWUcGAGLC5BWhaU3CqrGODAAYxuS3mtOahKgQyACYwJN1dMg7YHHcNUAITA4GeLLGfZjcmgSzEcgglcIONAgGkDZhB98mPxwgXJxZpBKBxmw8WUeHynq2MO9Z8t1s5D4QApODAZMLZdMrJALsaJDvZjO3xAKWEHagEXWFmFZUSMll8sMBwkVpi1Qi0EAcUVnPFuY9S76bjdIcQGKYXiERYEeDfDcbZw9AYlAhAenDXQ9ggukDZgGkCyUOgAkMmAVgkg+iTgAAAMB9paZFplAoaHNzU/V6XfV6PerkAEuh+wdAUjmOI8dxdHFxMdf2G0EQBCGnKVLzvgYcWKWwA43Hn3291PGjXOWUIAzAPOatvykRgBCYPM4k7ECBpeYBrBJ3LWAg09dLCUvYASSBEhA/3FVACAg0ZjM5b0xuaQOSikAGsWXy06/J3TNhY6l5AKtEIIPYMrmyxmxhBqhhB5AESkD8EMgglX71698stf9H3/vOilJyPyZXqCYHqIxxAeKHuxKxFWZl/eTLN0vtn+TWHpO79ACkDyUOYosKcbYwWzXCbjExuTUJQPxQUyC2wmwZePf5/lLHxv2FGaDSmgSkD3ctYivMloGox7jEmcktJiaPvwFwPwQygIHCDDZolQBgEt61hHvjnTyIG9OvGdPTD6wS71pC6MJuxqdQxqJMv2Z4DxWwuNRcmYVCQZubm6rX66rX61EnZy0ouBBHXJfRYPwQTOE4jhzH0cXFxVzb07WUYI8/+3qp/XmBHsIQ9nVpsjDvKfIdpqFrCaEjEAFWi/dQAYujJkowCi4AIzx4IKm4shOMggsAkHTUdADWipZCAKtEIANgrWgpnI0B9MDiuOoBICaYIg0sjkAGCAFP1gCwHsaVlp7nqdls6uHDhzo+PtbDhw91cHAQdbKACTxZ4z5MHj9E8I6oGHXl+L6vZrMp27bH33W7XdVqNTWbzQhTBiANwq6sw6zMw047wTuiYlQgc3h4qFqtNvFdsVhUo9GIKEUwVdiF+rvP95c6Pu6Hyno2k9MO3MaoQMbzPHW7XVWr1Ynvd3d3I0oRTBV2of7kyzehHh/TUVlHx+RuMZjNqECmUCioVqtpd3dX5XJZktTv95XJZCJOGYB5MZYiGmEHGpwXRMW4l0bm83n1+32Vy2XVajV1Op2JMTPXpfmlkZgt7MqUynq2MF9eyHkFkiOxL43s9XoqlUpqt9tqt9vq9Xpz7TccDic+b21taWtrK4wkwgBhVyhUWNHgvALmOj8/1/n5+fjz9Xp7FuNaZBqNhrLZrCSNB/42m80b42ZGRhHddc+fP9eLFy9CSyeWx9NvMnFeo0G+I+5evHihL7744sb3d7XIGBXIjIKYUdDieZ4qlYr6/b56vZ5yudyNfUaBzPv37ycyghaZ+AuzCwJIG+4nxN20FplHjx4lq2up3W5rMBiMP1uWpV6vp3w+r1evXk0NZEa2t7cZI7NiPOHhPrhuAExz3wYGY0oEz/NkWdbUnz179kxv375dc4oQ9lRXpnMmE1Oko8H9hKQyJpCxLEue50392cnJiQqFwppThLDx5A2sDvcTksqoMTJHR0c6Pj6emG7teZ4ajYZardbUfZh+HR66CHAfYV43pl+TpqcfWKV562+jAhnpcpxMp9MZL4J310sjCWSQRFR405k+oNX09AOrlNh1ZMrl8nhVXyCtGGcCAJeMC2QAYBbTB7Sann4gCgQygIGo8KYzvcvM9PQDUeCuAQxEhQcAlz6IOgEAAAD3xWMdgMRgNld0yHtEhSsHQGIwmys6Yec9gRJm4cwCAGKPIBWzpCaQKRQK2tzcVL1eV71ejzo5AELAbK7okPdYFcdx5DiOLi4u5treuJV9F8XKvsBiaMJHHHFdpk9iV/YFEC7GOiCOOO+YhSsDwFox1gHAKhHIAJjAWAcAJiGQATAh7CZ8kwMlusWA+OGuArBWJlfmdIsB8cMrCgAAgLHMfTQCgDULu1uMritgcVz1ADCnsAMFuq6AxdG1BAAAjEWLDIAJdG9Ex+QZXUBUKHEATKB7IzoEgcDiuGsAIAVoaUNScWUCQArQ0oakYrAvAAAw1kYQBEHUiQjT6DXg3//+97W5ual6va56vR51spBwJjfjm5x2zMZ5hSkcx5HjOLq4uNDPf/5znZ2daXt7e+b2qQlk7soIYJUef/b1UvvTjA8g7eatvwmxAawVLQMAVokSAQgB64HMxqBTAKtEIAOEgFYDAFgPSlsAa0VrFYBVIpABsFa0VgFYJdaRAQAAxiKQAQAAxqKNF0BiMLUbSB/uWgCJwdRuIH3oWgIAAMaiRQZAYjC1G0gfAhkAicEYFyB9uOsBJAaDfYH0Sc1dWygUtLm5qXq9rnq9HnVyAISAwb6A+RzHkeM4uri4mGv7jSAIgpDTFKl5XwMOwHyPP/t6qf0JZID4mLf+Tk2LDIDkY7AvkD4EMgASgzEuQPqwjgwAADAWjy8AMKdf/fo3S+3/0fe+s6KUABghkAGAOT358s1S+981mJjp48DiuOoBICaYPg4sjkAGMBBP7tOFnS/vPt9f6vgAVs/I0szzPDWbTT18+FDHx8cqFAoql8tRJwtYG57cpws7X8Ie48L0cWBxxgUy3W5Xtm2r1Wopk8nI8zyVSiUVi0VlMpmokwcA95bUljIgTHPfNf/8z/+sTCajx48fh5icu1UqFb1582YctHiep5OTk0jTBKwbT+7TkS9A+swdyHQ6HW1sbOiv/uqvwkzPrY6OjmRZlnK53Pi7YrGo09PTyNIERIEn9+nCzhfGJgHxs9Rd9dlnn+mjjz6aGdz89Kc/1Y9+9CNZliXbtvUXf/EXy/w6NZvNiSAGANaJsUlA/My1su/Lly/V6XRujEHJZrO3jkt5+vSpqtWq/uVf/kX/9E//tFxKddmNZFmWXNeV67o6OjrS0dHRXPsOh8OJf+fn50unBwAArMb5+fmNunoec7XIPH36VPv7+/rhD3+oH//4x/rwww8lSWdnZzo+Pp66z09/+lOdnZ2p0WhIugx6luH7vqTLYKZWq8myLEmX3U2VSkWtVuvW/R89ejTx+fnz53rx4sVSaQKQLozBAcJzeHioL774YuH95u5asixLf/Znf6bHjx/rxz/+sSzL0ldffaVSqTQzQcVicTw4eHd3d+HEXXV1QO8oiJGkarWqRqOhfr9/a7fT+/fvJ14DvrW1tVR6AKQPY1yA8Dx79kyffvrp+PNwOLzRCDHNQndls9lUJpPRy5cv5fu+bNtWLpfTs2fPdHh4ON7uzZs36vf76vf74+/+8i//cpFfdcMoeCkUChPfj7q2ut3urYHM9vb2RCADAHHDYGKk2dbW1r0aGRa+6m3blm3bN77/wQ9+oB/84AeSpK+++kq1Wk1/8id/snCC7mswGKztdwFAGBhMDCxursG+d9nf39c//MM/aGdnR7/61a/08uVL/d3f/d0qDj0hl8vNHJOz7BgcAABgnpW2Qz59+nSVh7uhVqvdGNTreZ4k8YoCAMZjMDGwuI0gCIKoE7GIbDarVqs1Hg9TqVS0u7urZrM5dfvhcKidnR2dnZ0xRgYAAEPMW38bNzKs1+up0Wgok8nI930VCgUdHBxEnSwAABAB41pkFkWLDAAwIwrmSWyLDABgccyIQlKtZNYSAABAFGiRAYAUYEYUkopABggB4xGSyeTzyjWFpOLKRiqFXSExHiGZOK9A/BDIIJWokBBHJrf4AFHhqgdCwHiEZAr7vBJgA4tLTSBTKBS0ubmper2uer0edXIQsbArJJ6Mk4nzCoTPcRw5jqOLi4u5tmdBPACICbqWgN9jQTwAMAyBCLA4FsQDAADGIpABAADGoh0TALA0xvcgKlw5AIClMXUcUaFrCQAAGIsWGcBANOMjblgEElGhNAMMRDM+4obgGFGhawkAABiLEBowEM34AHCJQAYwEM34AHCJriUAAGAsAhkAAGCs1AQyhUJBe3t7chwn6qQAAIAZHMfR3t6eCoXCXNtvBEEQhJymSM37GnAAABAf89bfqWmRAQAAyUMgAwAAjEUgAwAAjEUgAwAAjEUgAwAAjEUgAwAAjEUgAwAAjEUgAwAAjMWb5wBM+NWvf7PU/h997zsrSgkA3I1ABsCEJ1++WWr/X/ztn68oJQBwN7qWAACAsWiRATDh3ef7UScBAOZGIANgAmNcAJiEriUAAGCs1AQyhUJBe3t7chwn6qQAAIAZHMfR3t6eCoXCXNtvBEEQhJymSA2HQ+3s7Ojs7Ezb29tRJwdIvX/79/9Yav/v/gE94kAazFt/UyIAWKu9n/z9UvszvRvAVanpWgIAAMlDiwyAtfrZ3/ww6iQASBACGQBrxRgXAKtEiQIYyOQBsyanHUD8GF8ilEoldTqdqJMBrJXJA2ZNTjuA+DE6kHFdV91uN+pkADfQ6gAA62Fsaen7vlqtVtTJAKYKu9XB5AGzJqcdQPwYG8i4rqtarUaLDGAYWpsArJKRJYrnebIsK+pkADOF3eoQZosP3WIATGJkidNut3VwcKB2ux11UoCpTK40zz/gAAAPA0lEQVTMGYwLwCTGlbbdblflcjnqZACRYpwJAFwyLpDp9/sqFosL7zccDic+b21taWtra1XJgmHC7j4J+/hhtvgQJAGIwvn5uc7Pz8efr9fbsxgVyLiuq2q1eq99Hz16NPH5+fPnevHixQpShbCEGQyE3X1icveMyd1iAMx1eHioL774YuH9jCmxfN+XJGUymXvt//79+4nXgNMaE38mBwOYjcHEAKZ59uyZPv300/Hn4XB4oxFiGmNKhG63q16vp1qtNv7O8zxJUq1WUyaTkW3bM/ff3t6eCGSQbmF3n9A9MxsBKoBp7jvkYyMIgiCE9KzFaC2Z2/6E4XConZ0dnZ2dEcgYhif3ZHr82ddL7U8gA6TDvPW30SX9qLsJyUQgkky0VgFYJSNrCs/zZNv2eFXfUqmkSqVy74HAANaHABXAKhndtTQPupaAxZjcpWdy2gFMSkXXEoDVM3kwrslpB3A/H0SdAAAAgPuiRQbABAbjAjAJgQyACSaPEyEIA9LH3BILAK4xOQgDcD+MkQEAAMYikAEAAMYikAEAAMZKTSBTKBS0t7cnx3GiTgoAAJjBcRzt7e2pUCjMtT0r+wJYK1bfBTAPVvYFEEusvgtglVLTtQQAAJKHFhkAa8WidQBWiUAGwFoxxgXAKtG1BAAAjEUgAwAAjEUgAwAAjEUgAwAAjEUgAwAAjEUgAwAAjMU8SAATeIUAAJNQ4gCYwCsEAJiEQAYwEK0mAHCJ0gwwUJitJrxCAIBJUhPIFAoFbW5uql6vq16vR50cILZorQEQJcdx5DiOLi4u5tp+IwiCIOQ0RWo4HGpnZ0dnZ2fa3t6OOjnAStC1BCDp5q2/Kc0AAxGIAMAl1pEBAADGIpABAADGon0aAOYU9tgkxj4Bi+OqB4A5hb1YIIsRAoujawkAABiL6dcAMCe6loD1Yfo1AKxY2IECgQiwOLqWAACAsQhkAACAsQhkAACAsQhkAACAsRhZBiAxmPUDpE9q7tpCoaDNzU3V63XV6/WokwMgBCwoB5jPcRw5jqOLi4u5tmcdGQCJ8fizr5fan0AGiA/WkQGQOj/7mx9GnQQAa0YgAyAxGOMCpA+zlgAAgLEIZAAAgLEIZAAAgLEIZAAAgLEIZAAAgLGMG+Lf7XbV6XTk+748z1OlUlG1Wo06WQAAIAJGBTL9fl/9fl+2bUuSfN/XJ598ol6vp2azGXHqAADAuhnVtdRsNnVwcDD+nMlkZNu2XNeV53kRpgwAAETBqEDm9evXOjo6mvjuyZMnki67nAAAQLoYFcjs7u7q+Pg46mQAAICYMGqMzGAwuPHdqEtp1DIDAADSw6hAZppms6lisahcLnfrdsPhcOLz1taWtra2wkwaAACY0/n5uc7Pz8efr9fbsxjVtXRdu92W53lqtVp3bvvo0SPt7OyM/x0eHq4hhQAAYB6Hh4cT9fSjR4/m2m8jCIIg5LSFwvd95fN5dTodWZY1c7vhcKidnR29f/9e29vb4+9pkQEAID6mtcg8evRIZ2dnE/X3dcZ2LVUqlTuDmKu2t7dvzQgAABCd+zYwGNm1VKvVZNv2RBDT7/cjTBEAAIiCcYGM67qqVCoTg3s9z2NBPAAAUsiorqVut6tWq6VSqTTRAtPpdMavLQAAAOlh1GDfBw8eyPf9qT+b9WeMBvveNVgIAADEx7z1t1EtMqenp1EnAQAAxIhxY2QAAABGCGQAAICxCGQAAICxCGQAAICxCGQAAICxCGQAAICxCGQAAICxUhPIFAoF7e3tyXGcqJMCAABmcBxHe3t7KhQKc21v1Mq+98HKvgAAmGfe+js1LTIAACB5CGQAAICxCGQAAICxCGQAAICxCGQAAICxCGQAAICxCGQAAICxCGQAAICxCGQAAICxCGQAAICxCGQAAICxCGQAAICxCGQAAICxUhPIFAoF7e3tyXGcqJMCAABmcBxHe3t7KhQKc22/EQRBEHKaIjXva8ABAEB8zFt/p6ZFBgAAJA+BDAAAMBaBDAAAMBaBDAAAMBaBDAAAMBaBDAAAMBaBDAAAMBaBDAAAMBaBDAAAMBaBDAAAMBaBDAAAMBaBDAAAMBaBDAAAMBaBDAAAMFZqAplCoaC9vT05jhN1UgAAwAyO42hvb0+FQmGu7TeCIAhCTlOkhsOhdnZ2dHZ2pu3t7aiTAwAA5jBv/Z2aFhkAAJA8BDIAAMBYBDIAAMBYBDIAAMBYBDIAAMBYBDIAAMBYBDIAAMBY34o6AYvyPE+2bSubzSqTyUiSqtVqxKkCAABRMKpFxvd9lUol2batg4MDVatV9Xo9ua4bddJCcX5+rhcvXuj8/DzqpGCFOK/JxHlNJs5r/Bm1sm+j0ZDv+2o2m+PvPM9TPp/X6enp1H1MXtnX5LRjNs5rMnFek4nzGp1EruzbbreVz+cnvrMsS77vq9/vR5Imk9/dFHbaTT6+yedVMjtvTE572EzOG9OPHybO65ICg0gKOp3Oje8zmUxg2/bUfc7OzgJJwdnZWShp+uM//uNQjhsEZqfd9OObfF6DwNy8Cfv4nNdojm3y8Tmv0R1/3rw3ZrCv7/szf7a7u6vj4+OpPwv+s+fsX//1XzUcDsffb21taWtra+l0XVxcTBx3lUbHDev4Yabd9OObfF4lc/Mm7ONzXqM5tsnH57yu7/jn5+cTY5G++eYbSb+vx2cxJpA5OTm59eezAp1RRuzt7a08TSM7OzuhHVuSHj16FNqxw067ycc3+bxKZueNyWnnvCbz+JzX6I7/zTff3Po7jAlk7uvjjz/WYDDQt7/9bW1sbIy/X1WLDAAAWN71FpkgCPTb3/5WH3/88a37GRfITGt5OTk5Ga8pc90HH3wgy7LCThYAAIiAMbOWdnd3JU3vYvJ9Xw8fPlx3kgAAQMSMaZHJZDLjqdbTFIvFNacIAJBUpVJJnU5n4jtWlo8nYwIZSSqXyxoMBhPfjdaPyeVyUSRpad1uV51OR77vy/M8VSqVGzdGPp+Xbdt68uSJJOn169fyfV8HBwdRJBlzGJ3Xhw8fajAYKJPJyLbtiW0oFM00z7nlnjWb67rqdrsT341Wlu/1euP7tVaryXVd7tuohTYBPASnp6eBZVnB6enp+LtqtRq0Wq0IU3V/vV5vYv2b09PTIJPJBNVqdWI7SRP/isXiupOKBXQ6naDZbE58Z1lWUC6Xx59nXcvX90O8zHNug4B71mSnp6dBsVgMrlePBwcHN8rmwWAQZDKZdSYPUxj1igLp8im22Wwqm81qMBgom80aGw3XarWJ1y1Il08CtVpNg8FgPEi5Vqspn8/L930Vi0VjW5/SolarqdvtTrQejp7cRrfbfV63gejNc25H33HPmuno6EiWZalSqUyc02w2q0ajcaO+2djYUK/X4xxHyKiuJenylQTXm3FN9fr1a2Wz2Ynm5lFTdLfbHd8wJgdraVQqle5c96jdbqvRaEx8d/V1GxSK8TTPuZW4Z03led7MWa6zfpbJZNTtdrlnI2TMrKUkum1FYpirXC6r1WpNfNftdlUul8ef7yoUEU/znFuYq91uTz2X911ZHuthXItMklwfuCxdVnDS71tmJOn4+FhHR0fKZDLjG4pBg+ZwXVeS9PLlS0kUikly/dyOcM+a57aA9L4ry2M9CGRiptls3uhTH81uGanVamo0GonpYkuqbrerVquld+/eqdPpjGc6UCiab9a5HeGeNU+/32cZD0PRtRQj7XZbnufdaLq+/rlSqejo6IgKL+aKxaKazabevHmjUqk0fnqH+e46t9yzZpl3CvWiK8tjPQhkYsL3fTUajalPd9ddHRCM+MtkMmo0GqrVauN1jyQKxSSYdW6v456Nr9F9eNt9x8ry8UYgExOVSkWdTufGANDRdM9pRuNpEC/9fv9GkDKqyF69ekWhaLC7zq3EPWuabrerXq+nWq02/jdaFmHUJcjK8vHGGJkYqNVqsm17IogZTcF1XVf5fH5i+1EFyHS/+PF9X/l8XtVq9cYaQSMUimaa59xK4p41TLlcvjHId7Sy79XznMSV5ZOCFpmIua6rSqVyY3Dv6Mnt4ODgRt9tt9tVJpOZmNmEeBgFKbVabeL70fkslUqSKBRNNO+55Z4137SHjGfPnqnb7U78rNls3hgPhfUzbmXfJOl2u7Jte1wAjnQ6Hdm2rVwuN75xRk8Mo6dC27ZZuyKmjo6Obsw8y+fz2t3dHb+EbnQer7+3pVQqcV5jbJ5zyz1rrtFss263K8/zVCwWJ95/l6SV5ZOEQCZCDx48mNm9cPW0XH2x5MnJiWq1Gt0PMddut8cDtz3PU6FQuLGOCIWimeY5t9yzwPoQyAAAAGMxRgYAABiLQAYAABiLQAYAABiLQAYAABiLQAYAABiLQAZAol1fpwlAshDIAEiso6MjXtQIJByBDIBYazQayufz2tjYUDabVaVSUbvdvnO/fr8/XkmXlzUCyUUgAyDWbNsev9+o0+mo1WrNvdS/ZVnjFXgBJBOBDIDY63Q6sixr4g3xtxm9PV66DGYIZIDk+lbUCQCAu3S73bnfVeR5nrrdrt69ezf+7vqbxgEkB4EMgFjzfV++7889+8j3/YmXOA4GA1pkgASjawlArI1mHc3TIuO67rhLaSSbzRLIAAlGIAMg1jqdjjKZzK3jYzzPU6VSkW3bE0FLt9tVs9lUv9/X0dHROpILYM02giAIok4EAMySzWaVy+XUarWiTgqAGKJFBkBs+b4vz/NUKBSiTgqAmCKQARBbo5lHs8bHuK67zuQAiCECGQCx1el0JOnGAN4RplUDIJABEFu3rR/jui5dTgAIZADEk+/7Eyv0XuW6rhqNxtyvKgCQXCyIByB2arXaeP2Ybrc7fteS53l69+6dfN9XtVq9sV+j0ZDrunrz5o0sy9LJycncrzUAYCamXwNIBM/z5Hmenjx5otevX+vJkyczx9YASA4CGQAAYCzGyABIJN/3o04CgDUgkAGQCP1+X67ryvM8+b6vw8PDqJMEYA3oWgJgvNH7lXZ3d7W/vy9JarVaDPQFUoBABgAAGIuuJQAAYCwCGQAAYCwCGQAAYCwCGQAAYCwCGQAAYCwCGQAAYCwCGQAAYCwCGQAAYCwCGQAAYCwCGQAAYCwCGQAAYCwCGQAAYKz/D6hXOAaGdaHCAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x13e612f50>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots()\n",
    "ax[:plot](ent_spec[:, 1], ent_spec[:, 2], ls=\"none\", marker=\"_\", ms=\"12\", mew=\"1.5\")\n",
    "ax[:set_xlabel](\"\\$L_z^A\\$\")\n",
    "ax[:set_ylabel](\"\\$\\\\xi\\$\")\n",
    "ax[:set_title](\"\\$N=$(N), N_\\\\phi=$(Nphi); N_A=$(NA), N_A^\\\\mathrm{orb}=$(length(subsystemA))\\$\")\n",
    "ax[:set_xlim](23, 44)\n",
    "ax[:set_ylim](0, 12);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate and plot pair correlation function along equator of sphere (great circle)\n",
    "$g(\\mathbf{r}) \\equiv \\frac{1}{\\rho^2}\\langle\\hat\\psi^\\dagger(\\mathbf{r})\\,\\hat\\psi^\\dagger(0)\\,\\hat\\psi(0)\\,\\hat\\psi(\\mathbf{r})\\rangle $, where $\\hat\\psi(\\mathbf{r})$ is the electron field operator in the continuum projected into the LLL, and $\\rho$ is the 2D density of electrons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing 4-point functions ...\n",
      "Working on alpha2p = 1: 1/22\n",
      "Working on alpha2p = 2: 2/22\n",
      "Working on alpha2p = 3: 3/22\n",
      "Working on alpha2p = 4: 4/22\n",
      "Working on alpha2p = 5: 5/22\n",
      "Working on alpha2p = 6: 6/22\n",
      "Working on alpha2p = 7: 7/22\n",
      "Working on alpha2p = 8: 8/22\n",
      "Working on alpha2p = 9: 9/22\n",
      "Working on alpha2p = 10: 10/22\n",
      "Working on alpha2p = 11: 11/22\n",
      "Working on alpha2p = 12: 12/22\n",
      "Working on alpha2p = 13: 13/22\n",
      "Working on alpha2p = 14: 14/22\n",
      "Working on alpha2p = 15: 15/22\n",
      "Working on alpha2p = 16: 16/22\n",
      "Working on alpha2p = 17: 17/22\n",
      "Working on alpha2p = 18: 18/22\n",
      "Working on alpha2p = 19: 19/22\n",
      "Working on alpha2p = 20: 20/22\n",
      "Working on alpha2p = 21: 21/22\n",
      "Working on alpha2p = 22: 22/22\n",
      "elapsed time: 4.589514643 seconds\n",
      "\n",
      "Computing pair correlation function for 200 points ...\n",
      "Working on g[i], i = 1: 1/200 = 0.5%\n",
      "Working on g[i], i = 10: 10/200 = 5.0%\n",
      "Working on g[i], i = 20: 20/200 = 10.0%\n",
      "Working on g[i], i = 30: 30/200 = 15.0%\n",
      "Working on g[i], i = 40: 40/200 = 20.0%\n",
      "Working on g[i], i = 50: 50/200 = 25.0%\n",
      "Working on g[i], i = 60: 60/200 = 30.0%\n",
      "Working on g[i], i = 70: 70/200 = 35.0%\n",
      "Working on g[i], i = 80: 80/200 = 40.0%\n",
      "Working on g[i], i = 90: 90/200 = 45.0%\n",
      "Working on g[i], i = 100: 100/200 = 50.0%\n",
      "Working on g[i], i = 110: 110/200 = 55.0%\n",
      "Working on g[i], i = 120: 120/200 = 60.0%\n",
      "Working on g[i], i = 130: 130/200 = 65.0%\n",
      "Working on g[i], i = 140: 140/200 = 70.0%\n",
      "Working on g[i], i = 150: 150/200 = 75.0%\n",
      "Working on g[i], i = 160: 160/200 = 80.0%\n",
      "Working on g[i], i = 170: 170/200 = 85.0%\n",
      "Working on g[i], i = 180: 180/200 = 90.0%\n",
      "Working on g[i], i = 190: 190/200 = 95.0%\n",
      "Working on g[i], i = 200: 200/200 = 100.0%\n",
      "elapsed time: 101.029179439 seconds\n",
      "\n",
      "Throwing out imaginary part of g: maximum(abs, imag(g)) = 2.3603740944573712e-17\n"
     ]
    }
   ],
   "source": [
    "g = pair_correlation_function_equator_gc(psi, setup.states, pts=200, frac=0.5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAHOCAYAAACRo6NSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xdc03f+B/BX2JsAihOBxBnrIkSr3QK16zoU7LS1rUJbf9fr2RbKtXftXe+OYq83m1bSacdVhNrldQh2t1YlcW8JKA5QRhJkk3x/f1BSI8MAgW/G6/l4+GjJN99v3nwD+b74fD9DIgiCACIiIiI35CV2AURERESDhUGHiIiI3BaDDhEREbktBh0iIiJyWww6RERE5LYYdIiIiMhtMegQ0ZDR6XRIS0vr9TkajQZpaWlQKpVQKpVYtWrVEFXnPvpzDu15b4hcEYMOkZPrvGBJJBJIJBIYDIYenxcREQGJRAK5XI6MjIwhrrRnBoMBGo0GSqWyx/oBICsrC5GRkSgoKIBWq0VBQQFycnIgl8sdVos7nM/e9PUc2vveELkqBh0iJ1dQUIBXXnkFycnJAICcnJxen5eamorS0lLk5eUNZZk9SktLs7YUSKXSHp+n0+kAAKmpqdbHZDIZCgoKoNfrHRY0XP189qav59De94bIlUk4MzKR81u1ahVSU1OhVCoBAHV1dd0+T6PRIDExEQkJCUNZnt0iIiKQmJiIoqKiLtsyMjKQm5vb7QU3IiICBoMBjvq4cpfzeb6BnMPe3hsiV8YWHSIXUFpaCplMhsWLF8NgMKC4uLjH57nKRfl8JSUliI+P7/b2iUwmAwDo9XqHvJa7ns+hPIdEroJBh8iFZGVlAQByc3O73e7KfSwiIyNhMBhQW1s7ZK/pbudTjHNI5Ox8xC6AiHqn0+mst1hkMhkSEhJQXFwMg8Fgc4tCr9dbn+eKOvuRdLY8nKuzFaK7bX3lzudzqM4hkSthiw6RkysuLrZ2nAWA7OxsAF070Z7/PFcjlUq7vU3UGUIyMzMd8jrufD6H6hwSuRIGHSIn19mfpFPniBqNRtPr8+yRlZWFiIiIfv3rqV+Lo2VlZUEmk1kDyUAN5vm8ELHOt6PPIZEr4a0rIheUnp4OjUZj0+rQn/4kubm5PfZPcQarVq2CTqdDaWnpoA5/7u/5LCwstBnKfSFinO+hOodEzootOkROrKd+Iud3otXr9Q6dVM8Z6PV65OTkoKioyGEtK44+n87e6XcwziGRq2HQIXJiPfUTOb8TrU6nc7n+JBeSlpaGgoICh35fnnY+B+McErkaBh0iJ6bVanv8S7xzlluNRoNt27a51HwvF5KWlobs7GyHX6AdeT57Gt3kLAbrHBK5GgYdIifWW5+K9PR0AB2jhfo734szdkbOyspCSkpKl74vGo1mwJPdOeJ8ajQaaDQa5ObmQq/XIysry+66hup8D+Y5JHI17IxM5KTs6SfS2Ym2v/1znK0zcmFhIaKioqyh41xFRUVdHj9/7pveOOJ8pqSkICsrC8nJydBoNEhPT4fBYEBSUhK0Wu0FaxiK893Xc0jk7hh0iJxUbm6udcHFnmRkZECj0bjM7YneZu3V6/VYvnw5ZDIZ8vPzu+x3PqVSaR1NZM8tpIGez4yMDKSkpHTZJpVKrf16xL592NdzeP52Z+9cTdQvAhE5lfT0dEEmkwkABKlUKqSmpvb6/OTk5CGqrH8yMzOF5ORkQSqVCgAEAIJMJhOSk5OFgoIC6/MSEhKs27v7d/73mZqaKshkMiEvL6/X13fU+Tz341Kr1QpardZmW11dXa/HHQp9PYf2vjdEroyrlxORS+u8hTSYOoeld65yfu5rajQaFBQUcNVvIifFzshE5LJ0Ot2QjHySyWTWBTPPf/2CggIUFBQMeg1E1D9s0SEil5WVlTVknan1ej3y8vKgUqlQVFRknXiQnXuJnBuDDhG5pMLCQiQkJAz5XDY6nQ5SqdSp59Ahol/w1hURuaTU1FRRwkZJSQlDDpELYdAhIiIit8VbV0REROS22KJDREREbotBh4iIiNwWgw4RERG5LY9e68piseDkyZMIDQ2FRCIRuxwiIiKygyAIqK+vx+jRo+Hl1XubjdMGHZ1Oh/z8fLsnAysuLkZRUREMBgP0ej3S0tIuOJFXeXl5v1d9JiIiInFVVFRg7NixvT7HKYOOTqdDUlKS3TOO6nQ66HQ6aygyGAyIj4+HVqtFXl5ej/v5+voCAPbt24cxY8YMvPBuqFQqbNu2bVCOPVSvMZjHN5lMiImJQUVFBcLCwgblNQDXPkdD8RpD8T64+jka7OPzd8E5js/fBfGPb8970Pmc0NDQCx7PqYKOXq9HVlaWdV0Ze+Xl5dkEGqlUitzcXGRkZFiP153O21WhoaGD9gPt7e09qB9aQ/EaQ/E9hIWFufT34A7vMzC474M7nCNXfw8A1z9HQ/EeAPxdEPv4gH3vgT3dTpyqM7JMJkNBQQFyc3MhlUrt3m/dunVYtWqVzWOJiYkAOm5piWnFihUu/xpD8T0MNnc4R67+PrjDOXL19wBw/XPE98A5XsOV3gennTBQqVQiOTnZrj46crkcqampNs/V6XRQKpXIy8vr8RbY8ePHrc1jF7rHR4PDZDIhPDwcRqNxSP5Ko+7xfRAf3wPnwPdBfPa8B315n5zq1lV/lZaWdnlMr9cD+KVlpzv+/v42/6Wh5+/vj6effprvgcj4PoiP74Fz4PsgPke/B27RotOdlJQUAEBRUVGPz2FyJyIicj0e16JzvsLCQuj1emi1WruebzKZbL729/dnmiciInISLS0taGlpsX59/nW7N07VGdkRDAYDsrKyUFRUZHeH5piYGISHh1v/5eTkDHKVREREZK+cnByb63RMTIzd+7pdi05aWhqKiop6HFLenfPH6rM1h4iIyHlkZ2dj5cqV1q8759Gxh1sFnYyMDOTm5tqEHJ1Oh4SEhF73G+x5K4iIiKj/BtKlxG1uXWk0GqSlpdmEGr1ebx19RURERJ7HaVt0DAYDDAZDl8f1ej1SUlKQl5eH5ORkAB2TAhYUFCAlJQU6nc763KKion6P2iIiIiLX51RBx2AwICcnx9oSo9FoUFtbC5VKhczMTOvzamtrbfZLS0uDwWDodhbkC922IiIiIvfltPPoDAXOo0NEROR6+nL9dps+OkRERETnY9AhIiIil2G2CGhqNdv9fKfqo0NERESezdDYin0nTdhfWY+jNQ04WtOIKlMzahpaYWhsRZtZgKWl0e7jMegQERGRaE7XN+PbQ9XYoq/B1vJaHK2xP8TYg0GHiIiIhlR5dQM+2XkSRfursOu4scv2cZFBUIwKQ/zwYMRGBmG0NBCRwX6IDPZDgK83mhvPYsw/7XstBh0iIiIadI2t7fhox0ms3XoMO88LNzPGhmPe+GGYHR+JhHERCA/07fVYJrP98YVBh4iIiAZNWXUD3t58FAXaCtQ3twMAvL0kuGT8MFw/bSSumhyN6NCAQXt9Bh0iIiJyuC36Gqi/LsW3h85YH4uNCsKSi2Nx86wxGBYyNAtoM+gAUKlU8Pb2xooVK7BixQqxyyEiInJZ24/V4e9Fh/Dd4WoAgEQCXDUpGnfPjcXlE4bDy0vS72Or1Wqo1WqYzfYPL+fMyJwZmYiIaMD2nDDi70WH8OWB0wAAHy8JFqti8MDlcoyLCnLoa/Xl+s0WHSIiIuq30/XNyP3sIN7XHQfQ0f9m4awxeDhpAmIiHRtw+oNBh4iIiPqszWzBW5uP4p9Fh1Df0tHJ+KaZo/FI8kTEDwsWubpfMOgQERFRn2wurcEzH+/Fwap6AMD0seH4000XYWaMVOTKumLQISIiIruYmtvwlw37kV9SAQCICPJF5jWTsTgxBt4D6GQ8mBh0iIiI6IK+O3wGWYW7cNLYDIkEuGP2ODy+YBKkQX5il9YrBh0iIiLq0dmWdvzlf/vx3tZjADqWZ3g+dTrmyKJErsw+DDpERETUrW3ltXhk7Q6cMDQBAJbOi0PmNZMQ5Oc68cF1KiUiIqIhYbYIUH91BP8sPgSLAIyNCMTzqTMwV+4arTjnYtAhIiIiqypTMx5ZuwOb9TUAgIWzxuBPN1+EEH/XjAyuWTURERE53FcHTuPRgp2obWhFkJ83nr3pIixSjhW7rAFh0CEiIvJwZouAFzYexEtflwIAFKPC8OIdsyAbHiJyZQPHoENELsdsEXC0pgGlZxqgP3MWp4zNOHO2BXUNrWhqM6OlzQIvL8DP2wuBft6ICPJDVLAfRksDERsVDNnwYMiGBcPH20vsb4VIdHUNrXh47XbrIpz3zI1F9nVTEODrLXJljsGgQ0ROr81swfZjBnx3+Ay0R+uws8KAhlb7Vy/ujr+PFyaPCkPCOCnmxEdhTnwkIoKdez4QIkfbc8KIjLe1OGFoQoCvF3IXTcdNM8eIXZZDcfXy8HBMnDgR3t7eWLFiBVasWCF2WUQEoKXdjO8OVeOTXSfx5f7T1rV0OgX6ekMeHQzZsBCMiQhEdKg/IoP9EODrDX8fLwgA2totaGw1o7ahFdVnW3C8rglHaxpw5PTZLkFJIgESxkVg/uRoXHvRSLdosifqTaH2OJ78YDda2i2IjQrC6ruUmDKq95XAxaZWq6FWq2E2m3Ho0CG7Vi9n0LFzmXciGhpHTp/Ff7ccw/rtx2FobLM+HhHki8smDMfFsigkxEoxITq031POWywCjtY2YtdxA7aV12KLvhaHT5+1ec6scVIsTBiLX00f5fQzvxL1hdki4LnP9uOV78oAAPMnR+Mft85EeKCvyJXZry/XbwYdBh0i0QmCgM2lNXj5m1JrPwEAGBHmj+unjcb100dhZox0UNfSOWlowlcHT2Pj3ip8f6QaZkvHR6OftxeSpkTjjjnjcOn4YZBInHM9HyJ7NLa245G1O7BxXxUA4DdJE/CbpAnwctJ1qnrCoGMnBh0icQmCgO8OV+OFokPYWWEAAHhJgPmTR+DOi8fh8gnDRVko8HR9Mz7ecRLv605g/ymT9fHJI0ORfrkMN0wfDT8fdmQm13K6vhnL1pRg13Ej/Ly98Hya6/bHYdCxE4MOkXh2VhiQ89l+/KSvBdDROfg2VQyWXSZDTGSQyNX9Yt9JE9aVVGBdSQUaf+7XMzIsAPdeEoc75oxDaIDrNPeT5zpYWY/73tyGE4YmRAT5QnN3IlRxkWKX1W8MOnZi0CEaemfqW/D8FwewruQ4gI5bQ3ddHIuHrpJjWIi/yNX1zNDYine3HMObP5bjTH0LgI5+Qw9dOR5L5sa6zVBccj/fHjqDFe/qUN/SjvhhwXhjqQpxw4LFLmtAGHTsxKBDNHQEQcDabRX466f7Ud/cMYJq4awxeHTBJIyRBopcnf1a2s34aMdJrP6mFPozDQCAUeEB+E3SBKQqx3JuHnIq7209hqc+3AOzRcDs+Ejk3aV0i2kUGHTsxKBDNDSO1TQi8/2d1ttU08aE45kbp0IZGyFyZf3XbrZgve4E/ll8CCeNzQAA2fBgZF0zGVcrRrDTMonKYhGw6ouDWP1Nx0zHt8wag+cWTYO/j3u0PDLo2IlBh2jwfbj9BJ76cA/OtrQj0Ncbj149EfdeEi9KJ+PB0NxmxrtbjkH91RHUNrQCAC6bMAzP3DgVcs7FQyJobjNj5bod+HR3JQDgkeSOkVXuFL4ZdOzEoEM0eM62tOMPH+3Bet0JAMDsuEi8sHiGU3U0dqT65jas/qYUr3xbhlazBb7eEtx3STx+nTTBZVd9JtdTfbYFy9aUYEeFAb7eEqxKnY5bZrn2opzdYdCxE4MO0eDYfdyIh9duR1l1A7wkwG+SJuL/5o93m1ac3pRXN+DZDfuw6cBpAEB0qD9+f4MCN0wf5VZ/UZPzOVxVj3vf3IbjdU2QBvki7y4l5siixC5rUDDo2IlBh8jx3t5cjj9t2Ic2s4DR4QH41+2zXHoYa399eaAKf/pkH8prGgEASZOj8ezNF2G0C3W8Jtfxw5FqPPCOFvXN7YiLCsLrS1VuvYwJg46dGHSIHKfNbMEfP9mLd346BgC4ZupIPLdomkcvn9DSbsbqr/V48avDaDMLCPbzRta1k3HXnFiXm4mWnNe6bRX43Qe70W4RkBgbAc3diYh0g5FVvWHQsRODDpFjGBpb8dC7OvxYWgOJBHjimslIv1zGWzU/O1xVjyfW74b2aB0AQBkbgdxF0zA+OlTkysiVWSwCXig6CPVXHSOrbpwxGqtSp3vEnE4MOnZi0CEauCOnz2LZmm0or2lEsJ83/nXbLCQrRohdltOxWAS8u+UonvvsABpazfD38ULmNZNx77w4tu5QnzW3mfFYwU5s2HUKAPDr+eOxMmWix/xx0ZfrN2e2AqBSqaBQKKBWq8Uuhcil/HCkGre89APKaxoxNiIQ7z80jyGnB15eEiyZG4eilVfgionD0dJuwbMb9uGOV39CRW2j2OWRC6k524I7X92CDbtOwcdLgudTp+PRqyd5RMhRq9VQKBRQqVR278MWHbboEPXLp7tP4Tdrt6PNLEAVF4HVdykR5cRLODgTQRDw363H8Jf/7Udjqxkh/j74ww0KpCWO9YiLFfVf6ZmzuPeNbThW24iwAB+svkuJeeOHiV3WkOOtKzsx6BD1z3tbj+HJD3bDIgDXTxuFv986w21mXB1KR2sa8Oi6nSj5ue9O0uRo5CyahujQAJErI2e0ubQGD7yjhbGpDTGRgXhjqcpj+3nx1hURDZqXvy5F9vqOkHP77HH49+2zGHL6KTYqGPkZc5F97WT4eXth04HTWPCPb/Hp7lNil0ZO5n3tcdz9+hYYm9owa5wUHzx0iceGnL5iiw5bdIjsIggCnvvsAPK+1QMAHrpSjscXeEa/gKFwoNKElfk7se+UCQCwMGEM/njjVIQG+IpcGYlJEAT8o/gw/r3pMICOFtQXFs/wiJFVveGtKzsx6BDZx2IR8PuP9uDdLR1z5PzuuslIv1wuclXup7Xdgn9vOoyXvj4CiwCMjQjEP26d6ZETLlLHPEyZhbvw0Y6TADr+uHjs6kkcpQcGHbsx6BBdmCB0hJx3fjoGiQR4buE03KoaJ3ZZbq2kvBa/XbcDFbVN8JIAD105Hr9JngBfb/Y28BS1Da3IeLsE28rr4OMlwV9uuYi/d+dgHx0icghBEPD0x3utIef51Bn8sB0CiXGR+PThy7AoYSwsAvDiV0ew6OUfUXrmrNil0RAoPXMWC1/6AdvK6xDq74M3753N37sBcNoWHZ1Oh/z8fOTm5tr1fL1ej9zcXMjlckilUgBAenp6r/uwRYeoZ4Ig4I+f7MObP5ZDIgFWLZqOtMQYscvyOP/bdQq/+2A3jE1tCPD1wlPXK3DnnHHsG+Wmvj9cjYfe1cLU3I4x0kC8ca8KE0ew0/H5XL5FR6fTISkpye7nGwwGpKSkIDc3F5mZmUhPT4dWq4VGoxnEKonclyAI+NOGjpADALkLGXLEcv30Ufjikctx6fhhaG6z4KkP92DZmhJUn20RuzRysLc3l+OeN7bC1NyOhHFSfLjiEoYcB3CqoKPX65GWlob8/HxERtrf+S4nJwfJycnWlhwAyMrKQlZW1mCUSeT2nvv8AN74obzj/xdOw2IVQ46YRoYH4K37ZuP3Nyjg59MxDP2af36LTfurxC6NHKDdbMEfPtqD33+0F2aLgIWzxuC/yy/G8FBOwOkIThV0ZDIZCgoKkJubaxNaLqSwsBBKpbLLsQwGA3Q6naPLJHJrL39dirxvOoaQ//WWabhtNvsGOAMvLwnuvzQeH//fJZg8MhTVZ1tx/5oSPPnBbjS2totdHvWTsbEN9765DW9tPgqJBMi8ZhKHjzuYUwWd/tLr9ZDJZF0el0qlKC4uFqEiItf03y3HkPv5AQDAk9dNwR1zGHKczeSRYfhwxSVYdmk8AODdLcdww7+/x67jBpEro74qq27ALS//gO8OVyPQ1xur71LioSvHs/+Vg7l80DEYev7ljoyMRE1NzRBWQ+S6Nuw6iSc/3A0AWHGVHMsv7/rHAzmHAF9vPHWDAu8um4ORYQHQVzdg4Us/4sUvD8NsccrxJXSeH45U42b1D9CfacDo8AAUPjgXC6aOFLsst+QjdgEDVVtb2+v23oJQpyNHjiAkJMT6dVhYGEaOHIn29nacOXOmy/NHjRoFAKiurkZbW5vNNqlUisDAQDQ0NMBkMtls8/PzQ1RUFCwWC6qqut5bj46Ohre3N2pra9HSYtvRMDQ0FCEhIWhqauryPfn4+GD48OEAgFOnuk4dP2zYMPj6+sJgMKCpqclmW3BwMMLCwtDS0tLlXHp5eWHEiI6VqKuqqmCxWGy2R0ZGwt/fHyaTCQ0NDTbbAgMDIZVK0dbWhurq6i41dZ7DM2fOoL3dttm98xyePXsW9fX1Ntv8/f0RGRkJs9mM06dPdznuiBEj4OXlhZqaGrS2ttpsCwsLQ3BwcLfn0NfXF8OGdSyM1905HD58OHx8fFBXV4fm5mabbSEhIQgNDe32HHp7eyM6OhpA9+cwKioKfn5+3Z7DoKAghIeHd3sOJRIJRo7s+FDs7hxGREQgICCg23MYEBCAiIgIm3O4ubQGjxbsQJtZwJL5M/DY1ZO6PYfh4eEICgpCY2MjjEajzbbOn29BEFBZWdnlHHb+fHd3Djt/vpubm1FXV2ez7dyf78rKSpw/ULTz59toNKKx0XYV8M6f79bW1i5/9Jz783369GmYzWab7Z0/3/X19Th71nZYd+c5FPszQhYMvHnreDz32QF8qa/H3zYewqbdx/FkyjiMiQjq9hzyM0Lcz4iamhq8u+Uo/vNlx6SQ08ZGYM2KSxAdGuDUnxHnGjlyJCQSyZB+Rvj5+cHX1xfNzc0wGAxdfid74/JBxxHO79+zePFi5Ofnw2QyIS8vr8vzn3nmGQDAhx9+iOPHj9tsW7hwIaZPn469e/fi008/tdkml8uxZMkStLW1dXvcxx9/HMHBwfjiiy9w8OBBm20LFizA3LlzodfrUVBQYLNt1KhRyMjIAAC8+uqrXT6wH3roIURHR+Pbb7/t0mfp0ksvRXJyMk6dOoU333zTZltYWBhWrlwJAHj33Xe7fCgvXboUcXFx2Lp1K77//nubbQkJCbjxxhtRV1fX5Xv19vbG73//ewDA+vXru3xopKWlYerUqdi9eze++OILm22TJk3C7bffjubm5m7PYXZ2Nvz9/fHpp5+itLTUZtt1112H2bNn4/Dhw1i/fr3NtrFjx2LZsmUA0O1xH374YURGRuKrr77Crl27bLZdeeWVuPLKK1FRUYF33nnHZltkZCQefvhhAMCaNWu6XIjvv/9+xMTEYPPmzdi8ebPNNpVKheuvvx7V1dVdavL390d2djYAYN26dV0utrfffjsmTZqE7du3Y9OmTTbbFAoFFi9ejIaGBuTl5eGkoQnrdSfQbrFg4ohQPHPDbZBIJPjkk09QXl5us++NN96IhIQEHDhwAB9//LHNtri4OCxduhRms7nbc7hy5UqEhYWhqKgI+/bts9mWlJSEyy67DEePHsV7771ns2348OFYsWIFAOCNN97o8kdARkYGRo0ahe+//x7btm2z2TZ37lwsWLAAVVVVeO2112y2BQUFITMzEwCwdu3aLhfxu+66C+PHj4dWq8XXX39ts2369OlYuHCh03xGxAoCHpiWgP+W+2DLrv24bv1ruHJSNKaMCoVEIuFnxDnE/Iw4WFqG+5/6Bw6f7ggWU0aFIVk+GdGhNwBw3s+I8z311FPw8fEZ0s+Io0ePdvn5s5fTzqOjVCqRnJx8wXl09Ho95HI5CgoKkJqaarMtIiIC6enpPR6jcxy+Vqtliw7/WvPIFp2S/eW4f81WGJraMVcWhb+lzcC4saOH/K81tuh0GOhnRF2rFx5+ewu2He4IV8mTo/HEtZMRFRbEz4ififUZcaZZguVvbMbBY5Xw9ZZgZcokLEoYAx8fH6f+jHDmFh2lUunaS0DYG3QMBgMiIiKQl5fXZYJAiURinVunO5wwkDxZzdkWLHz5RxytacSMseF4L/1iBPmxkdfVmS0CVn9Tin8UHUK7RcCIMH+8kDYTl04YJnZpHuuLvZV4bN1O1Le0IzrUHy/flQBlLNcvGwiXnzCwL6RSqXUoeXeSk5OHuCIi59fUasb9a0pwtKYRMZGBePUeFUOOm/D2kmDFVeOx/qF5kA0LRpWpBXe9tgXPbtiH5jbzhQ9ADtPSbsYzH+9Fxtta1Le0QxUXgQ0PX8qQM8RcPugAQGpqapd7rZ33mRMSEsQoichpmS0CHsnfjh0VBoQH+uKNpbM5MZkbmj5Wig0PX4q7Lu6YIuC178tw04s/YO9J4wX2JEcoq27Aopd/tM4uvvyyePx3+cWIDg0QtzAP5LRBx2AwdNtK09kn59z5cbKzs1FcXGzz/Ly8vC6ddokI+PP/9uGLvVXw8/bCK3cnYnx0yIV3IpcU5OeDP988Da/dk4ioYD8crKrHjS/+gFWfH2DrziD6aMcJ3PDv77DnhAkRQb54fWkinrxewdXnReJUbdUGgwE5OTnQ6/XQ6/XQaDSora2FSqWy6Wdzfoc4qVSKoqIi5OTkQC6Xo7S0FEqlskvnZCJP99r3ZdalHV5YPAOz49mE7gmSpozA549cjqc/3oNPd1fipa9L8fmeSjy3aDp/BhzI1NyGP368D+/rOjqDz46PxL9um4lR4YEiV+bZnLYz8lBgZ2TyJJ/vOYUH39VBEIDsaycj4wq52CWRCD7fU4k/fLQHp+s7Rm3dOWccsq6djLAAX5Erc23fH65GZuFOnDQ2QyIBfj1/Ah6ePx4+bMUZFH25fjPoMOiQB9AercMdr/yElnYLllwciz/dNJXTzHswY1Mbcj7dj7XbKgAAw0L88eT1k3HzzDH8ueijxtZ2PPfZAby1+SgAIDYqCC+kzUBiHFvKBhODjp0YdMgTVNQ24ib1D6htaEXylGisvkvJvzIJAPBjaTWe/GAPyqo75meZHReJP908FZNH8vPQHlv0Nch8fxeO1nTM37Tk4lhkXzeZIxiHAIOOnRh0yN3VN7dh0cs/4lDVWUwbE478DM6VQ7Za2s149bsy/OfLw2hus8DbS4J75sbhkZQJvJ3Vg5qz3ib2AAAgAElEQVSzLfjrpwesfXFGhQdgVep0XDZhuMiVeQ4GHTsx6JA7M1sEpL9Vgk0HTiM61B8f/9+lGBnOoa3UvROGJjz7yT58vrdjxtqIIF88nDQBd86JhZ8PWwABwGIRsHZbBXI/PwBjU8eM17fPHocnrp2M8ECGwqHEoGMnBh1yZzmf7UfeN3r4+3ghP2MuZsZIxS6JXMDXB0/j2Q37UHqm43ZWbFQQMhdMxnXTRnp0/53dx434w8d7sP1YxzQmU0aF4S+3XISEcREiV+aZGHTsxKBD7qpQexyPFewEAPzrtpm4aeYYkSsiV9JutiC/pAL/KDqM6rMdo7Nmxkjx2NWTcMn4KI8KPOXVDfjbxoPYsKtjfatgP2+svHoS7pkby75uImLQsRODDrkj7dFa3K7ZglazBb+ePx6PXj1J7JLIRTW0tOOV7/TQfKtHY2vHBIMzY6R4OGk8rpoU7daB50x9C/7z5WH8d8sxtFsESCTATTNG44lrp/AWsBNg0LFT54maOHEivL29sWLFCqxYsULssoj67YShCTe9+D2qz7ZiwdQRePlOJby83PdiREPjdH0zXvqqFO9tPYaW9o7VtaeODsOv54/H1YqRbvUzVmVqxus/lOHtzUet4e6KicORec0kTB0dLnJ1pFaroVarYTabcejQIQadC2GLDrmThpZ2pK7ejP2nTJgyKgyFD8xFsD9HWJHjnKlvwavf6fH2T7+EgJjIQNw1JxaLE2MQEewncoX9d7iqHppv9fhwxwm0mTsuizPGhiPr2smYJ+fK786GLTp2YtAhd2GxCHjwXS2+2FuFYSF++Oj/LsUYKaedp8FR19CK138ow5ofy2FqbgcA+Pt44VczRuPuubGYPtY1Or63my347kg13t58FF8eOG19XBUXgfTL5Uie4t6351wZg46dGHTIXbyw8SD+8+UR+Hl74b30OVDGclZWGnxNrWZ8tOME3tp8FPtOmayPT4gOwQ3TR+NXM0ZBNtz5Fo09VFWP97XH8cH2E9alMCQSYIFiJNKvkHEklQtg0LETgw65g492nMBv1u4AAPwtbQZSlWNFrog8jSAI0B2rw1ubj+Kz3ZVoNVus2xSjwnDDjFG4cmI0Jo8MFaU/j8UiYN8pE745dAaf76nE7hNG67aIIF/cNHMM7pkXh/hhwUNeG/UPg46dGHTI1e2oMODWvM1oabcg43IZsq+bInZJ5OGMTW0o2leFT3aexPdHqmG2/HKJiQr2w8XyKFwiH4Z58ijERgUN2q2hSmMztpbX4puDZ/DNoTPWYfIA4OMlwfzJ0VikHIurJkVzQkQXxKBjJwYdcmWVxmbc+OL3OF3fgqTJ0dDcnQhvNxr9Qq6vtqEVn++pxBd7K7G1rBZNbWab7aH+Ppg4MhSTRoZi8shQTBwRijHSQIwIC7ArfFgsAqobWnC8rgkVtY04XHUWB6vqsfu4EZWmZpvnBvl5Y558GK6cNBzXXjQSUSH+Dv1eaWgx6NiJQYdcVVOrGbdqNmPXcSMmjgjB+w/OQyjXJSIn1tpuwY4KA34srcaPR2qwvaLOOrqpO+GBvggL9EGovy98fbzg4yWBRRDQ2m5BU5sZxsY2GJrabFqMzuUlASaPDMNlE4bhionDoYyLgL+P92B9ezTEGHTsxKBDrkgQBPz6ve3YsOsUIoJ88dGKSzEuKkjssoj6pLXdgrLqBhyoNOFgZT0OVtbjyJmzOGVsRmu75cIH+JmXBBgZFoCxEUGQR4dg4ogQKEaFYdrYcC5g68b6cv3mTwGRi/nPl0ewYdcp+HhJsPouJUMOuSQ/Hy9M+vm21bkEQUBdYxtqG1phbGpDfXNHq02bWYCXpGO/AF9vSIN8ERHkh8hgP/hyKQbqBYMOkQv5bPcp/L3oEADgzzdfhDmyKJErInIsiUSCyOCOAEPkCIzBRC5izwkjVq7rWKjz3kvicNvscSJXRETk/Bh0iFzA6fpmLH+rBE1tZlw+cTie5DByIiK7MOgQObnmNjMy3tbilLEZsuHB+M/ts+DDPglERHbhpyWRExMEAb9bvxvbjxkQHuiL1+5RITyQw8iJiOzFoANApVJBoVBArVaLXQqRjdXf6LF++wl4e0nw0p0JnKKeiDyaWq2GQqGASqWyex/Oo8N5dMhJfbG3Eg+8o4UgAM/eNBVL5saJXRIRkVPoy/WbLTpETmjvSSN+m78DggDcPTeWIYeIqJ8YdIiczOn6ZixfU4LGVjMuHT8Mf7hBIXZJREQui0GHyIl0jrA6aWyGbFgw1HckcIQVEdEA8BOUyEkIgoAn3t/1ywirpSqEB3GEFRHRQDDoEDmJl74uxYc7TsLHS4KXOcKKiMghGHSInMDne07h+S8OAgD+eNNUzBs/TOSKiIjcA4MOkcj2nDDit/kda1gtnReHO+fEilwREZH7YNAhEtFpUzOWrfllDaunrucaVkREjsSgQySS5jYzlr9VgkpTM8ZHh+DFO7iGFRGRo/FTlUgEgiDgsYKd2HncCGmQL167JxFhARxhRUTkaAw6RCL4Z/FhbNh16ucRVkrERnGEFRHRYGDQIRpi72uP41+bDgMA/nzzRZgrjxK5IiIi98WgQzSEfiytxhPrdwEAHrhCjttmjxO5IiIi98agA0ClUkGhUECtVotdCrmxI6frkfG2Fm1mAddPH4XMBZPELomIyKWo1WooFAqoVCq795EIgiAMYk1OrS/LvBMNxJn6Ftzy0g84XtcEZWwE3l02BwG+3mKXRUTkkvpy/WaLDtEga2o1Y9lbJThe14TYqCC8cnciQw4R0RBh0CEaRBaLgN/m78DOCgOkQb54Y6kKkcF+YpdFROQxGHSIBlHOZ/vx+d5K+Hl7QbMkEbLhIWKXRETkURh0iAbJ25vL8cp3ZQCA59OmY3Z8pLgFERF5IAYdokGwcW8lnv54LwDg0ZSJuGnmGJErIiLyTAw6RA6mPVqLX7+3HRYBWJw4Fv83f7zYJREReSwGHSIHOnL6LO5fU4KWdgvmT47GX2+ZBolEInZZREQey0fsArqj1+uRm5sLuVwOqVQKAEhPT7drv7y8PERFRaGmpgZRUVHIzMwc7HKJAABVpmbc8/pWGBrbMCNGytXIiYicgNMFHYPBgJSUFGi1WmvIycjIgEaj6TXsGAwG5OXlITc31/pYcXExMjIykJeXN+h1k2czNbfhnte34oShCfHDgvH6PYkI8nO6Xy8iIo/jdH9u5uTkIDk52RpyACArKwtZWVkX3C8jI8PmseTkZJSUlAxKnUSdWtrNyHhLiwOV9RgW4o+37puNqBB/scsiIiI4YdApLCyEUqm0eUwmk8FgMECn0/W4n16vR3FxcZfHIyM5pJcGj8Ui4LGCXdisr0GwnzfevFeFmMggscsiIqKfOV3Q0ev1kMlkXR6XSqXdBplOKpUKGRkZKCwstD6m0+lsWoaIHO2vn+7HJztPwsdLgtVLlLhoTLjYJRER0TmcqhOBwWDocVtkZCRqamp63J6ZmYn8/HykpaUhNTUVGRkZKCoqQkFBwQVf12Qy2Xzt7+8Pf3/eeqDe5X1Tile/75gQ8G9pM3DZhOEiV0RE5J5aWlrQ0tJi/fr863ZvnKpFp7a2ttftvQUhANBqtUhOTkZhYSFSUlJw66232vW6MTExCA8Pt/7Lycmxu2byTGu3HkPOZwcAAL+7bjJunsUJAYmIBktOTo7NdTomJsbufZ0q6AxUVlYW0tLSrKOslEolNBrNBferqKiA0Wi0/svOzh7sUsmF/W/XKWR/sBsA8OCVcqRfLhe5IiIi95adnW1zna6oqLB7X6e6ddWpu5ab2traXvvbZGVlQS6XW4egJycnIy0tDRkZGUhMTERCQkKP+4aFhSEsLGzghZPb+/bQGTySvx2CANw+exwyF0wSuyQiIrc3kC4lThV0OkdIdXcLy2AwICoqqsd9CwsLUVpaav1aJpNBq9VCqVQiPz+/16BDZA/t0VpkvK1Fm1nA9dNH4c83X8RZj4mInJxT3bqSSqXWoeTdSU5O7vbxnkZqAeBtKHKI/adMuPeNbWhqM+OKicPxj8Uz4e3FkENE5OycKugAQGpqqk3LDADr/Dk9tcrIZDLo9fput9XW1kKlUjm2SPIo5dUNWPLaVpia25EYG4GX70qAn4/T/eoQEVE3nO7TOjs7G8XFxTatOnl5eTbDxPV6PeRyuc28OhkZGV1mT9br9SgqKkJqaurgF05uqdLYjLte24Lqsy2YMioMry1VcWkHIiIX4nSf2FKpFEVFRcjJyYFcLkdpaSmUSmWXsHJ+P57MzEwUFhYiIyPD2mk5KirKrnl0iLpT19CKJa9twfG6JsRFBeGt+2YjPNBX7LKIiKgPJIIgCGIXIRaTyYTw8HAYjUaOuiIbZ1vaceerW7CzwoARYf4ofGAel3YgInISfbl+O92tKyKxNbeZkf5WCXZWGCAN8sU7989hyCEiclEMOkTnaDdb8PB72/FjaecinbMxYUSo2GUREVE/MegQ/cxiEfDE+t3YuK8Kfj5eeOXuRMyM4aKwRESujEGHCIAgCPjz//ajUHsc3l4SvHj7LMwbP0zssoiIaIAYdIgAvPjlEbz+Q8dK5LmLpuPqqSNFroiIiByBQYc83pofy/FC0SEAwB9uUCBVOVbkioiIyFEYdMijfbj9BJ7+eC8A4OGkCbjv0niRKyIiIkdi0AGgUqmgUCigVqvFLoWGUPG+KjxasBMAsHReHH6bPEHkioiIqDdqtRoKhaJPSztxwkBOGOiRftLX4J7Xt6Kl3YJbZo3BC2kz4MVFOomIXAInDCTqxZ4TRixbU4KWdguSp0RjVep0hhwiIjfFoEMepay6Afe8vhVnW9oxOz4SL96RAF9v/hoQEbkrfsKTx6gyNWPJa1tQ09AKxagwvHpPIgJ8vcUui4iIBhGDDnkEY2Mb7n5tq3Ul8jX3zUZYAFciJyJydww65PaaWs24b802HKyqR3SoP96+fw6Gh/qLXRYREQ0BBh1ya21mCx58Vwvt0TqEBfjgrftncyVyIiIPwqBDbstiEfB4wU58ffAMAny98Ma9KkweyWkEiIg8CYMOuSVBEPCnDfvw4Y6T8PGS4OU7lVDGRopdFhERDTEGHXJLq7/R480fywEAf0ubgasmR4tbEBERiYJBh9zOh9tPIPfzAwA6Fum8edYYkSsiIiKxMOiQW/nxSDUeL+xYv2r5ZfFcpJOIyMMx6JDbOFBpQsbbWrSZBVw/fRSyr50idklERCQyBh1yC6eMTVj6+jbUt7RjdlwkF+kkIiIADDoAAJVKBYVCAbVaLXYp1A+m5jYsfX0bKk3NGB8dAs3dSi7tQETkhtRqNRQKBVQqld37SARBEAaxJqfWl2XeyTm1my24b00Jvj10BsND/fHBQ/MwNoITAhIRubO+XL/ZokMu7c//249vD3VMCPj6PSqGHCIissGgQy7rnZ+OWufK+fvimZg2NlzcgoiIyOkw6JBL+uFINZ7+eC8A4LGrJ+K6aaNEroiIiJwRgw65HP2Zs3jwHS3MFgE3zxyNFVeNF7skIiJyUgw65FKMjW24f00JTM3tmDVOiucWTYdEwmHkRETUPQYdchkWi4BH8rejrLoBY6SB0CxJ5DByIiLqFYMOuYx/bTqMrw6egb+PF/KWKDE81F/skoiIyMkx6JBL+PJAFf616TAA4K+3TMNFYzjCioiILoxBh5ze0ZoGPLJ2BwBgycWxWKQcK3JFRETkKhh0yKk1tZqR8bYWpuZ2JIyT4vc3KMQuiYiIXIjPQA+wY8cO6PV61NbWAgBkMhlkMhni4uIGemjycIIgIHv9LhyorMewED+8dKcSfj7M5kREZL8+Bx2TyYS8vDzk5+ejrKwM8fHxiIyMhFQqBQAYDAbU1tairKwMiYmJSEtLw+LFi7mWFPXZf7cew4c7TsLbS4IX70jAyPAAsUsiIiIX06eg8/zzz2PdunXIyMhAQUEB4uPje31+WVkZiouLMX/+fDzwwANYtmzZgIolz3Gwsh5/+mQfACDrmkm4WBYlckVEROSK7LoPYDQasXjxYiiVSmzbtg3Lli27YMgBgPj4eCxfvhwlJSWIj4/Hgw8+CJPJNOCiHU2lUkGhUECtVotdCgFobjPj1+/p0NJuwRUTh2PZpTKxSyIiIiegVquhUCigUqns3kciCIJwoSe98sorWLx4McLDBzak12g0oqCgwGladvqyzDsNnSc/2I13txzDsBB/fP7IZRgWwvlyiIjoF325ftsVdNwVg47z+Wz3KTz4rg4A8Pb9s3HZhOEiV0RERM6mL9fvAQ9hccZbUeSaThiakPX+LgDAA1fIGXKIiGjABhR0HnjgAURERKC8vNxB5ZCnajdb8Mja7TA1t2NGjBSPXj1R7JKIiMgNDCjoKJVKbNy4kXPm0IDlfavHtvI6hPr74D+3zYKvN+fLISKigRvw1SQpKckRdZAHO1Bpwj+LDwEAnrlxKsZFBYlcERERuYsBBZ3ly5fjwQcfxM6dOx1VD3mYNrMFjxXsRJtZQPKUEViYMEbskoiIyI0MaAmIJ554Anl5edBoNJDJZEhOTkZKSgqSk5MHNIpJr9cjNzcXcrncOuNyenq63fvm5eUhKioKNTU1UKlUSE1N7XctNLhe/roUe06YIA3yxV8XXgSJRCJ2SURE5EYGFHSioqJgsVhQVlYGnU6HoqIiZGZmQq/XIyIiAunp6cjJyenTMQ0GA1JSUqDVaq0hJyMjAxqN5oJhp7i4GLm5uSgoKIBUKoVer7cGr85jkfPYe9KIf286DAD4441TER3KJR6IiMixBhR0ZDIZXn31VSxevBiLFi3CokWLAHRMDFhUVISysrI+HzMnJ6dLMMnKyoJSqbxg0ElLS8OmTZus+5672Cg5l9Z2Cx4r2IV2i4Brpo7EjTNGi10SERG5IYdMGLhp0yaHdUqWy+XIysrqEmokEgm0Wi0SEhK63W/VqlXIz8+HVqu1+7U4YaB4/l50CP/edBiRwX7Y+NvLOfsxERHZzeETBl5oUsC+hJwLHUuv10Mm67q2kVQqRXFxcY/75eXldbsfOZ+9J41Qf3UEAPDsTRcx5BAR0aCx69ZVfn4+5HI55s+fP6AX+/LLL6HX63tc68pgMPS4b2RkJGpqanrcrtfrkZqaCo1GY3OszMzMAVRMjma2CPjdB3tgtgi4btpIXD99lNglERGRG7OrRWf58uUoLS1FdnZ2v2ZBLi8vxwMPPACdTtfrgp4X6k/TUxDqfFyv1yM5ORnp6enWgJOWlnbB+kwmk82/lpaWC+5D/fPe1mPYWWFAiL8Pnv7VVLHLISIiF9DS0tLlWm0vuzsjL1++HGVlZUhPT4dEIrGOZpLJZF3uj5lMJpSUlECn02Hjxo2QSCRYvXo14uPj7f+u+uDcgHTu7av09HRkZWVBp9P12LcHAGJiYmy+fvrpp/HMM884vE5Pd6a+Bas+PwAAeOzqiRgRxlFWRER0YTk5OfjjH//Yr337NOoqPj4eGzduhNFoxLp166xDyfV6vc38J1KpFImJiUhLS0NBQQHCw8P7VFR3LTe1tbU9DhHvDDcqlcrm8c7nFxcX9xp0KioqbMKavz/7jAyGv366H6bmdlw0JgxL5saJXQ4REbmI7OxsrFy50vq1yWTq0kjRk34NLw8PD8fy5cuxfPly62NGo9G6rb8iIyMBdH8Ly2AwICoqql/HLS0t7XV7WFgYR10Nsh+PVOOD7ScgkQB/uXkavL04MSAREdnH39+/340QA5pH51wDCTidpFIpZDJZj31xkpOTe9w3ISGhx87Kcrl8wLVR/7W0m/HUh3sAAEsujsWMGE7eSEREQ8PplohOTU3t0gKj0+kAoNfbTxkZGdbnddLr9dZjkng03+ihr27AsBB/PHr1JLHLISIiDzKgoDN+/HioVCpkZ2dj/fr1feoF3ZPs7GwUFxfbtOrk5eWhoKDA+rVer4dcLreZVyc9PR16vd4m7HROPMj5dcRzrKYR//l5zpzf3zAF4YG+IldERESeZEC3rjIyMpCfn4+6ujprx2S5XI6EhATrqKy4uLg+HVMqlaKoqAg5OTmQy+UoLS2FUqns0irTXT8erVaLrKwsSKVSGAwGqFQqzqMjspzP9qO13YJ58igu80BERENuQEGntrYWJSUlNo9t2rQJRUVF1mHoQMeto1WrViE2Ntau48pkMuTm5va6va6ursvjUqkUeXl5ffgOaDBtLavFZ3sq4SUB/vArBVcmJyKiITegoNM5SupcSUlJSEpKQkpKCgwGA2QyGdauXYuEhASUl5cjNDR0IC9JLsJiEfDshn0AgFtV4zB5JEe1ERHR0BtwZ+T169d3+3hSUhLq6uowa9Ys5ObmYtu2bbyN5EE+3HECu08YEeLvg5UpE8Uuh4iIPNSAgs7jjz+O1atXY8GCBfjggw+6dEY+dyVxmUzW66gpch+Nre1Y9flBAMCKq8ZjeCgnYCQiInEMuEVn48aNSEpKwqJFixAREYGoqCioVCpERUVBqVTaPJd9NDyD5ls9Kk3NGBsRiHsviRO7HCIi8mASQRAERxzIaDRi27Zt2L59O4CODsjnrm0VGRmJW2+9FS+//LIjXs4hTCYTwsPDYTQaOTOyg1Qam3HV375GU5sZL94xCzdM50grIiJyrL5cvx06M3JycnKPsxdnZ2dzPhsP8PwXB9HUZoYyNgLXTxsldjlEROThHBZ0LuTxxx8fqpcikew5YcT7uuMAgN/fwOHkREQkPqdbAkIMKpUKCoUCarVa7FJc2gsbOzog3zhjNGZyPSsiInIwtVoNhUIBlUpl9z4O66PjithHx3G0R2ux6OXN8PaSYNPKKxA3LFjskoiIyE315frNFh0aMEEQ8PwXHa05acqxDDlEROQ0GHRowH44UoOf9LXw8/bCr5MmiF0OERGRFYMODYggCHj+5745d8wZhzHSQJErIiIi+gWDDg1I8f7T2FlhQKCvN1ZcNV7scoiIiGww6FC/WSyCdaTV0kviuNQDERE5HQYd6rcNu0/hQGU9Qv19kHE5J4MkIiLnw6BD/dJutuCfRYcAAMsvl0Ea5CdyRURERF0x6FC/fLzzJPTVDYgI8sV9l8ZfeAciIiIRMOhQn1ksAl76uhQAsOwyGUL8h2wlESIioj5h0KE++2JvJY6cPovQAB8smRsrdjlEREQ9YtChPhEEAS9+dQQAsHReHMICfEWuiIiIqGcMOtQnXx86g70nTQj09ca9l7BvDhEROTcGHbKbIAhQf9nRmnPnnHGIDOZIKyIicm4MOgBUKhUUCgXUarXYpTi1LWW1KDlaBz9vLyznvDlERDTE1Go1FAoFVCqV3ftIBEEQBrEmp9aXZd4JWPLaFnx3uBp3zhmHv9wyTexyiIjIQ/Xl+s0WHbLLzgoDvjtcDW8vCR64Qi52OURERHZh0CG7dI60umnmaMREBolcDRERkX0YdOiCDlXVo2hfFSQS4KEruUI5ERG5DgYduqBXvtUDABYoRmJ8dIjI1RAREdmPQYd6dbq+GR/tOAkAHGlFREQuh0GHevX25qNoNVswa5wUytgIscshIiLqEwYd6lFTqxnv/HQUALD8MrbmEBGR62HQoR69rzuOusY2xEQGYsHUkWKXQ0RE1GcMOtQti0XA69+XAQDuuyQe3l4SkSsiIiLqOwYd6tamA6ehr25AaIAP0hJjxC6HiIioXxh0qFuvfNcxpPyOOeMQ4u8jcjVERET9w6BDXew6bsDWslr4eEmwdF6c2OUQERH1G4MOdfHqdx19c341YzRGhQeKXA0REVH/MegAUKlUUCgUUKvVYpciupOGJvxv9ykAwP2XxotcDRER0S/UajUUCgVUKpXd+0gEQRAGsSan1pdl3j3F818cgPqrUlwsi8Ta9Llil0NERNRFX67fbNEhq+Y2M97bWgEA7JtDRERugUGHrP636xRqG1oxOjwAyVNGiF0OERHRgDHokNVbm8sBAHdeHAsfb/5oEBGR6+PVjAAAOyoM2HncCD9vL9ym4gSBRETkHhh0CADw1o/lAIAbZoxCVIi/uMUQERE5CIMOofpsCzbs6hhSfs/cOHGLISIiciCnnNtfr9cjNzcXcrkcUqkUAJCent7n46SkpKCoqMjR5bmd/G0VaDVbMCNGihkxUrHLISIichinCzoGgwEpKSnQarXWkJORkQGNRtOnsKPRaFBcXDxYZbqNdrMF7/x0FABwz9xYkashIiJyLKe7dZWTk4Pk5GRryAGArKwsZGVl2X0Mg8GAgoKCwSjP7RTvr8IpYzOigv1w3bRRYpdDRETkUE4XdAoLC6FUKm0ek8lkMBgM0Ol0dh1Do9EgIyNjMMpzO2t+7GjNuW12DAJ8vUWuhoiIyLGcLujo9XrIZLIuj0ulUrtuRfW0P3V1qKoem/U18JIAd87hbSsiInI/ThV0DAZDj9siIyNRU1NzwWMUFhYiNTXVkWW5rf9uOQYASJ4yAqOlXKWciIjcj1N1Rq6tre11e29BCACKi4v7FXJMJpPN1/7+/vD3d++5ZJpazXhfdxxAx0zIREREzqqlpQUtLS3Wr8+/bvfGqVp0Bkqn0/XrtlVMTAzCw8Ot/3JycgahOueyYddJ1De3IyYyEJeNHyZ2OURERD3KycmxuU7HxNg/g79Tteh06q7lpra21mYk1vn6Ovz8XBUVFTbLvLt7aw4A/Hdrx22r21Tj4OUlEbkaIiKinmVnZ2PlypXWr00mk91hx6mCTmRkJIDub2EZDAZERUV1u19nMOotCPUmLCzMJui4u/2nTNh+zAAfLwnSEseKXQ4REVGvBtKlxKmCjlQqtQ4l705ycnK3jxcXF0Or1doMKdfr9QA6JhuUSqXIzc11fMEuqrMT8tVTRyA6NEDkaoiIiAaPUwUdAEhNTUVpaanNY53z5yQkJPS4z/mdkDtnRs7LyxucQl1UY2s7Ptx+AgBwx2x2QiYiIjBLIBcAABf9SURBVPfmdJ2Rs7OzUVxcbNOqk5eXZzPTsV6vh1wu73VenQuN0PJUn+w8ifqWdsRFBWGevPtbgURERO7C6Vp0pFIpioqKkJOTA7lcjtLSUiiVyi4tNj0NRe9cELQzBKWkpCAtLa3fHZXdzbs/37a6fTY7IRMRkfuTCIIgiF2EWEwmE8LDw2E0Gj2iM/KeE0bc8J/v4efthc3Z8xEV4v6jy4iIyP305frtdLeuaPB0tuYsuGgkQw4REXkEBh0PcbalHR/v6OyEPE7kaoiIiIYGg46H2LDzJBpazZANC8bFskixyyEiIhoSDDoeYl1JBQBgsSoGEgk7IRMRkWdg0PEAR07XQ3fMAG8vCRYmjBG7HCIioiHDoOMB8rd1tObMnxzNmZCJiMijMOi4uTazBet1HZ2QFyfav9orERGRO2DQAaBSqaBQKKBWq8UuxeE27T+NmoZWDA/1x1WThotdDhERUb+p1WooFAqoVCq79+GEgW4+YeB9b27DlwdOI+MKGbKvnSJ2OURERAPGCQMJAFBlasbXB08D4G0rIiLyTAw6bqxQexwWAVDFRUA+PETscoiIiIYcg46bEgQBBT/PnZPG1hwiIvJQDDpuamtZLcprGhHs543rp40SuxwiIiJRMOi4qfyfW3N+NWM0gv19RK6GiIhIHAw6bqi+uQ2f7j4FgLetiIjIszHouKFPdp5Cc5sF46NDkDBOKnY5REREomHQcUOdt60WJ47lAp5EROTRGHTczMHKeuysMMDHS4KFCWPFLoeIiEhUDDpuZt3PrTlJU6IxLMRf5GqIiIjExaDjRlrbLfhgOxfwJCIi6sSg40Y27a9CbUMrokP9ccVELuBJRETEoONGOjshpyrHwsebby0RERGvhgBUKhUUCgXUarXYpfTbKWMTvj10BgDnziEiIvekVquhUCigUqns3kciCIIwiDU5tb4s8+7sXvzyMP628RBmx0diXcZcscshIiIaNH25frNFxw1YLALWlRwHwE7IRERE52LQcQNbympxrLYRIf4+uG7aSLHLISIichoMOm5g3TkLeAb5cQFPIiKiTgw6Ls7Y9MsCnosTORMyERHRuRh0XNwnO0+ipd2CiSNCMDOGC3gSERGdi0HHxa2zLuAZwwU8iYiIzsOg48L2nzJh13EjfL0luGXWGLHLISIicjoMOi6sszUnecoIRHEBTyIioi4YdFxUS7uZC3gSERFdAIOOiyraVwVDYxtGhgXgci7gSURE1C0GHRfVORPyIuUYeHuxEzIREVF3GHRc0PG6Rnx3uGMBT962IiIi6hmDzv+3d/cxjtz1Hcc/3stmLyS5m3iTS3PJ9nJ2HsjyFNbrAAnhIbGpKIIW1U6gpRENYFelIAqVraNSoVJh6xVClHaL7IoHQQnsrhsVENDiSYsoj3dnp4Q8QBIPSe6ScCHrndu7cNl7mv6xsYl37X062zO23y8pijzeGX9v5nb8ud/vN79fF8oXD8pxpFcFhrVr+Fy3ywEAwLMIOpLC4bBGR0c1NTXldilrOnXa0exz3VZvu47WHABA/5iamtLo6KjC4fC69/E5juO0sSZP28gy717xvQd/rds+t1fbtp6lvX8T0dbBLW6XBABAR23k+5sWnS4zvW9p7py3vvxSQg4AAGsg6HSRyjPH9Z37fyVJuiVMtxUAAGsh6HSRO0sHdeKUo5dcul0v2rnd7XIAAPA8gk6XcByntuTDrbTmAACwLgSdLnH3AVsPHjqqrYMDesu1O90uBwCArkDQ6RIzzw1C/v2XXKJtWwddrgYAgO5A0OkCzyye1Dd++oQk6VZmQgYAYN3OcruARizLUiaTUTAYlGEYkqREIrHmfqZpqlAoyLZtWZaleDy+rv287pv3PKlnjp/S7gvP1XW7/W6XAwBA1/Bc0LFtW9FoVMVisRZyksmkcrncqqGlVCqpVCopk8nUjrN7924Vi0Vls9mO1N4uX933mKSlda18PhbwBABgvTzXdTUxMaFIJFILOZKUTqeVTqdX3S+bzSqVStVeG4ahTCajXC4ny7LaVm+7PXToiEqP2doy4NMfhS51uxwAALqK54JOPp9XKBSq2xYIBGTbtkqlUtP9ZmZmNDk5WbdtfHxc0lKXVreqzoR80wt3aMf5W12uBgCA7uK5oGNZlgKBwIrthmGsGlj8fr/m5ubaWVrHHT95Wnfe/bgk6W3MnQMAwIZ5aoyObdtN31sryJTL5RXbql1W1ZadbmM+cEiVZ45rx/lDeu1VF7ldDgAAXcdTQadSqaz6/mpBqJFsNqtIJKKxsbFVf25hYaHu9dDQkIaGhjb0We1Q7baKj1+ms7Z4rvENAICOWFxc1OLiYu318u/t1fTst2c+n5dlWZqdnV3zZ0dGRrR9+/bafxMTEx2ocHWP28f0vYd+LWnpaSsAAPrVxMRE3ff0yMj6vxc91aJT1ajlplKp1D2Jtdb+6XRahUJhXfscOHBA27Ztq732QmvOV37ymBxHuj44rF3D57pdDgAArtmzZ48++MEP1l4vLCysO+x4Kuj4/UuT4TXqwrJtW8PDw+s6TjweV6FQaDiouZFt27bVBR23HT95Wl99rtvqHa/c5XI1AAC460yGlHiq68owjNqj5I1EIpE1j5FMJpXJZOpCzmqPpXtR4f5Devrooi46f0jR0YvdLgcAgK7lqaAjSbFYbMUTVNWgstag4lwup3g8XvdzlmV13YSB//bjRyUtPVI+yCBkAAA2zVNdV9JSP1woFJJt27XxNdlstm5QsWVZikajtaeqpKVJAWdnZxWNRutacAqFQm1ZiG7w8FNH9SNrTgM+6e3X/a7b5QAA0NU8F3QMw1ChUNDExISCwaDK5bJCoZBisVjdzy0fxxOPx2XbdsNJBddqCfKSO36ytK7VTS+8WDuNc1yuBgCA7uZzHMdxuwi3LCwsaPv27Tp8+LAnBiMfO35Kr/i4qYVnT+rzfxbW66/e4XZJAAB4zka+vxkA4iHfuOcJLTx7UiP+c/TaK5kJGQCAM0XQ8ZAvP9dt9cfX7dLAgM/lagAA6H4EHY+49/HD+ukBW4NbfIqPX+Z2OQAA9ASCjkd8+SdLj5S/8cWX6MLz3J+ZGQCAXkDQ8YCFZ0/oP+5+QpL0J6/gkXIAAFqFoOMB03sP6NiJU7rq4vN03W6/2+UAANAzCDouO3nqtL7ww0ckSbffsFs+H4OQAQBoFYKOpHA4rNHRUU1NTXX8s79z/yE9bh+T/9yz9Ycvv7Tjnw8AQLeYmprS6OiowuHwuvdhwkCXJwyMfeaH2v/ovN530xX60BuudqUGAAC6CRMGdomfHrC1/9F5DW7x6U9fucvtcgAA6DkEHRd97ge/lCS9+aU7tWPbVperAQCg9xB0XPKrw8/qm/c8KUm6/dW7Xa4GAIDeRNBxyRd/9IhOnnZ03W6/XnzpdrfLAQCgJxF0XHDs+CndsXdpXat30ZoDAEDbEHRccOfdB2X/5oRG/Ococs3FbpcDAEDPIuh02OnTjj7/g0ckSe+8fre2sEo5AABtQ9DpsP/++VN6+KmjOm/oLN3CKuUAALQVQaeDHMfR1HcfliS945W7dP7WQZcrAgCgtxF0OujHVkV3P2br7LMGdPurL3e7HAAAeh5Bp4P+5bnWnFvGL9OO85kgEACAdiPodMjPDh7W/z70tLYM+JR8TdDtcgAA6AsEnQ755/95SJL05pdeohH/C1yuBgCA/kDQkRQOhzU6Oqqpqam2HP++Jw7rv+47JJ9Peu/rr2jLZwAA0OumpqY0OjqqcDi87n18juM4bazJ0zayzPuZeM8X96tw/yG95WU79em3v7xtnwMAQD/YyPc3LTptdu/jh1W4/5AGfNL7b77S7XIAAOgrBJ02+5T5oCTpD669VFfsOM/lagAA6C8EnTYqPjov84GnNOCT3ncTY3MAAOg0gk6bOI6jf/j2A5KkeGhEgYtozQEAoNMIOm1iPvCU9j0yr62DA/qr6FVulwMAQF8i6LTByVOnlfnPn0uSbr9ht35nO7MgAwDgBoJOG8zsP6iHnzoq4wWD+vPXMQsyAABuIei02OHfnNAnvvMLSdL7brpS21ihHAAA1xB0WuyThV+o8sxxXbHjPN32ql1ulwMAQF8j6LTQA08u6Es/flSS9HdveZEGt3B6AQBwE9/ELeI4jj7ytft02pHe9JJLdMMVF7pdEgAAfY+g0yJf2XtAex+p6JzBLfrwm65xuxwAACCCTks8bh/Tx7+1NDngX//e1brUOMfligAAgETQkSSFw2GNjo5qampqw/s6jqMP3/kzHV08qdCuC/TO6y9vfYEAAEBTU1MaHR1VOBxe9z4+x3GcNtbkaRtZ5r2ZmX0HlPr3e3T2WQP61vtvZOFOAADabCPf37TonIGHDh3RR75+nyTpQ9GrCDkAAHgMQWeTjh0/pffeUdKxE6d045UX6j03BtwuCQAALEPQ2QTHcfS3X7tXDx46qovOH9Inb7lWAwM+t8sCAADLEHQ24bPf/6Vmiwc14JP+8W3X6qLzh9wuCQAANEDQ2aBv/+xJfey5R8n3vPEaXR9kYkAAALyKoLMB+x6p6APT/yfHkW571S69+8bdbpcEAABWQdBZpx88/LRu++xeLZ48rcg1O/SRN79IPh/jcgAA8LKz3C6gEcuylMlkFAwGZRiGJCmRSLRtv7WY9x/SX9xR0vGTp/Waqy7SP719TFsYfAwAgOd5rkXHtm1Fo1FlMhmlUiklEgkVi0XlcrmW77e4uFj3/+VOnXb06bseUuJL+3X85Gm9YfRi/ettIZ1z9pbN/wFRZ3FxUR/96EebXgN0BtfBfVwDb+A6uK/V18BzMyOn02nZtq1sNlvbZlmWQqGQ5ufnW7rfwYMHNTIyogMHDuiyyy6re+8J+5hS+Xv0/YefliTdOj6iv3/rizW4xXPZsKu1YnZqnDmug/u4Bt7AdXDfeq5BV8+MnM/nFQqF6rYFAgHZtq1SqdTy/Zazf3NcE996QK/7xHf1/Yef1jmDW/SJ+MuUib10UyFnM+tnee0zOvFnaLdeOEfdfh164Rx1+zWQuv8ccQ288RnddB08F3Qsy1IgsHKWYcMwZJpmy/eTpHufPKKZfQf0ri/s03Ufu0vZ71k6fvK0XrHbr6//5Q2KhS5bdf/V8BfaG3rhHHX7deiFc9Tt10Dq/nPENfDGZ3TTdfDUYGTbtpu+5/f7NTc319L9qr127/78Xg0MvaC2/aqLz9UHIlfrxisvlM/naGFhYT3lN3Tq1Kkz2t8Ln9HO41ePyzly9zM6cR26/Ry1+/j8Lnjj+PwuuH/8RtdgcXGxbszOkSNHJP32e3w1ngo6lUpl1febBZrN7nfixAlJ0uOfeWfd9gOS7tqz6iE3ZPv27a07mEuf0e7jj4yMtPX4Uvefo058RruvQy+co26/BlL3n6NO/D3id8H946/nGhw5cmTNOjwVdDrt8ssvV7lc1uDgYN2cOENDQxoaYlkHAAC8YHmLjuM4OnHihHbu3Lnmvp4MOo1aYCqVSm1unFbtNzAw0HBcDwAA6A2eGozs9/slNe6Ksm1bw8PDLd0PAAD0Nk8FHcMwao+ENxKJRFq6H4B60WjU7RIA9JFSqaR0Ot3wPcuylEwmNTk5qVwut+bEwc14rusqFoupXC7XbavOgzM2NtaS/dq1VAQ2xjRNFQoF2bYty7IUj8e5Di7K5XJrTsWA9rEsS9lsVsPDw5qbm1M4HFYsFnO7rL6x/PwPDw8rlUq5XVZPK5VKuvnmmxve96urHRSLxdr3dDKZVC6X2/j3hOMx8/PzTiAQcObn52vbEomEMzs7W3tdLpedQCDgFAqFDe232s9ls9l2/HHQRLFYdDKZTO31/Py8YxiGk0gkXKyqf83PzzuRSMTx4C2hLxQKBScSidTuS9V73PPvU2if+fl5J5VK1W0rFArcj9qkXC47sVjMSaVSTiAQWHHuHcdxUqnUivNfLpcdwzA2/HmevKuVy2UnlUo52Wy29v/l7xuGURd01rOf47T25GHzGt1AstmsI8kpl8suVNTfMpmMMzs7S9BxiWEYTrFYrL0uFAqOYRgEnQ5JpVIN7ztjY2MuVNNfxsbGGgadQCDQ8DtcUt3vynp4rutKWlq6IZPJrPp+o/Wr1tpPWloqYnl/4POXilitewytMzMzo2AwWNc0PD4+LmmpS4surM5pNqs4OmNyclKBQKDu3hOJRFZd2w+tZVlWw/tO9UEXdN5aqx1s5LvaU4ORO+FMlopA66w2YzU6K5/PMxbERdlslqDpsnA4rGQyqXw+X9tWKpXWnNIE7bHZ1Q6a8WSLTru0+uRh85YPHJeWQqj025YdtJ9pmoQcl1mWpVgsVnuipHqfYiBs56RSKU1PTysejysWiymZTKpQKGh2dtbt0vrSZlc7aKavWnRaffLQWtlsVpFIhO7DDiqVSrQmuKh6z7EsS5FIRIlEohZw4vG4m6X1nWKxqEgkonw+r2g0qltvvdXtktAifRV04F35fF6WZfEvqA7a1GOaaKnn/+Pr+YEzkUgon8/XpshA+6XTacXjcWWzWUlSKBTa9LwtaI3NrpKwXF8GnVadPLSGbdtKp9MqFApcgw6p/g5wvt1VDTfhcLhue/W6MG6wM9LptILBoBKJhBKJhMrlssbGxpRMJgmbLmj1agd9NUaHpSK8KR6Pq1Ao0IXSQaZpqlgsKplM1rZVx0glk0kZhrHmE4xov0Zj2dB6+Xy+7lwHAgEVi0WFQiFNT0/Tnd5hrV7toK+CDktFeE8ymVQmk6kLOTzm336xWGzFIOTqzMjVpnt0xtjYWNMHIYLBYIer6T+rTa+wZ88e7du3r8MVQdr8KgmN9F3XVStPHs5MLpdTPB6vO++WZdVaFtBZDMZ3R6PukervAE/EtV8gEGh6z6lUKiu6FdFatm03vPfs2bNHpmnWvZfNZjc1jtP33EyDfcO2bYVCoRXrZ0SjUW4qHWSapjKZzIpFJAuFgjKZDKGzg6prv5mmWXv6h3XHOisYDGp2drb29z4ej8vv99O61iGTk5Oam5ur6661LEvpdJoHJNrAtm1NTEzIsqza3EWxWEzhcLhuWoXq+mPBYFDlcrk2jmqj+i7oSK07edi8Cy64oGkLQh/+lUSfqw7INwxDtm2vmDUc7ZfP5+seiGBRz97Rl0EHAAD0h74bowMAAPoHQQcAAPQsgg4AAOhZBB0AANCzCDoAAKBnEXQA9CQmnwQgEXQA9CjW6gIgEXQA9Kj9+/c3XcPIsixNTk6y7AXQBwg6AHpOqVRqukhvdcmLVCqlmZmZDlcGoNP6avVyAP0hm80qnU43fC8ej+uuu+6SxEKmQD+gRQdAz2nWbZXP5zU+Pl5bU6pZ1xaA3kHQAdBTTNNs2m2VTqeVTCYlSTMzM4rFYp0sDYAL6LoC4Fn5fF6WZWlubk6ZTEa5XE6SVCwWlc1mG+4zOzvbsNuqVCqpUqlobGyMLiugjxB0AHhWpVJRKpWSz+eTbdvKZrMyTVPJZFKZTEaGYazYp1m31fT0dK2lJ5fLKZVKtb1+AO7zOY7juF0EACxnmqbGx8dVqVQUDAZVLBY1NjYmaenJqUZhxjRNFQqFhnPoBINBpdNp+f1+GYZR171lWZby+bwCgYAsy1IsFmP8DtAjaNEB4EnVIGKapgKBQC3kSGoaQpp1W0mqzZJcqVRWjM2JRqMql8uSlp7Euvnmm1UsFs/4zwDAfQxGBuBphUKh6eDi5Zp1W1VDTrlcViKRqHuvVCrVva4+kcXyEUBvIOgA8DTTNBWNRtf1c80CkW3bMgyjYZfW/v37V4z1MQxjRQAC0J0IOgA8q9qysp4WndnZ2dqj48tNT0837e6ybVt+v3/F9kqlsrFiAXgSQQeAZ5mmKcMwGj5dtVyzAcqlUqlhS45pmpKWWm+Wh5pm4QdA9yHoAPAsy7JWjKlpxDTNusHKz1fdHolEZJqmSqWS0um0xsfHJUnj4+MN59VpdjwA3YXHywF0vWQyqXQ6veYj4aZpyu/3rwgxwWCw7qmrUChUew2guxF0AHS9UCh0Ro+Dl0olTU9PKxwOa9++fUomk8yjA/QIgg6ArrbaJIEAwBgdAF1ttaetAICgA6CrVSoVupkANEXXFQAA6Fm06AAAgJ5F0AEAAD2LoAMAAHoWQQcAAPQsgg4AAOhZBB0AANCzCDoAAKBnEXQAAEDPIugAAICeRdABAAA96/8BuaWBIBk1qdgAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x121c3fb10>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots()\n",
    "ax[:plot](g[:, 1], g[:, 2])\n",
    "ax[:set_xlabel](\"\\$r/\\\\ell_0\\$\")\n",
    "ax[:set_ylabel](\"\\$g(r)\\$\")\n",
    "ax[:set_title](\"\\$N=$(N), N_\\\\phi=$(Nphi)\\$\")\n",
    "ax[:axhline](1, c=\"k\", ls=\"--\", lw=1, alpha=0.5)\n",
    "ax[:set_xlim](0, maximum(g[:, 1]))\n",
    "ax[:set_ylim](bottom=0);\n",
    "# fig[:savefig](\"g.pdf\", transparent=true, bbox_inches=\"tight\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}