SquareWell.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-fd9a8a256ce96e0b",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-360083120ef5ab52",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# 1D square well\n",
    "\n",
    "Consider a quantum particle, described by the wave function $\\psi(x)$ moving in an external potential $V(x)$.\n",
    "For simplicity, only consider a one-dimensional real variable $x$.\n",
    "\n",
    "The wave function of a particle obeys the Schoedinger equation\n",
    "\n",
    "\n",
    "$$\n",
    "-\\frac{\\hbar^2}{2m} \\psi''(x)  + V(x)\\psi(x) = E\\psi(x)\\;.\n",
    "$$\n",
    "\n",
    "Here $\\hbar$ is the Planck's constant, $m$ is the mass of the particle (both are known), and $E$ is the energy of the particle.\n",
    "\n",
    "Given boundary conditions, the non-zero solution for $\\psi(x)$ _only exists for some values of the energies_. Physisists say that the energy of a particle is _quantized_. Mathematicians say that the problem is nothing but the eigenvalue problem of the Sturm-Liouville operator\n",
    "\n",
    "\n",
    "$$\n",
    "L = -\\frac{\\hbar^2}{2m}\\frac{d^2}{dx^2} + V(x) \\;.\n",
    "$$\n",
    "\n",
    "\n",
    "Here we consider the so-called square well potential, \n",
    "\n",
    "$$\n",
    "V(x) = %\n",
    "\\begin{cases}\n",
    "&0, &0 < x < a \\\\\n",
    "&\\infty, &\\text{otherwise.} \n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Define\n",
    "$$\n",
    "E = \\frac{\\hbar^2}{2m} k^2\n",
    "$$\n",
    "\n",
    "The Schroedinger equation for $0 < x < a$ becomes\n",
    "\n",
    "$$\n",
    "\\psi''(x) = -k^2 \\psi(x)\n",
    "$$\n",
    "\n",
    "with the boundary conditions $\\psi(0) = \\psi(a) = 0$.\n",
    "\n",
    "Solutions are\n",
    "\n",
    "$$\\psi_n(x) = \\sqrt{\\frac{2}{a}} \\sin{k_n x}$$\n",
    "\n",
    "\n",
    "with\n",
    "\n",
    "$$\n",
    "k_n a = \\pi n\\;, \\qquad n = 1, 2, \\cdots\n",
    "$$\n",
    "\n",
    "So that the eigenenergies are\n",
    "\n",
    "$$\n",
    "E_n = \\frac{\\pi^2\\hbar^2}{2ma^2} n^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-3be39ca2ce355116",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Eigenfunctions, numerically\n",
    "\n",
    "Define a mesh\n",
    "\n",
    "$$x_j = a \\frac{j}{N}, \\qquad j=0, \\cdots, N\n",
    "$$\n",
    "\n",
    "and the mesh function $\\psi_j$.\n",
    "\n",
    "The Schroedinger equation becomes\n",
    "\n",
    "$$\n",
    "\\psi_{j-1} + \\psi_{j+1} - 2\\psi_{j} = -k^2 h^2 \\psi_j \\; , \\qquad j = 1, \\cdots, N-1\n",
    "$$\n",
    "\n",
    "where $h=a/N$ is the grid spacing, and the boundary conditions are $\\psi_0 = \\psi_N = 0$.\n",
    "\n",
    "Note that the last equation defines the eigenvalue problem in the $N-1$ dimensional space spanned by vectors $(\\psi_1, \\psi_2, \\cdots, \\psi_{N-1})^T$. Here the matrix, whose eigenvalues we are looking for, is tridiagonal, and the allowed values of $(k h)^2$ are its eigenvalues."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-eab0189cf5a42ef7",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Write the function that returns the matrix for the eigenvalue problem above.\n",
    "\n",
    "The numpy function `np.fill_diagonal` may be handy. Also note that to fill sub- and superdiagonals, you may slice the matrix, like so: `a[1:, :]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-aafe722b76d0a6b5",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def hamiltonian(npts):\n",
    "    \"\"\"Construct the eigenvalue matrix for the discretized 1D Schroedinger equation.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    npts : int\n",
    "        The number of nodes in the mesh: $x_j = a*j /npts$, with $j=0, ..., npts$\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    ham : ndarray, shape(npts-1, npts-1)\n",
    "        The matrix of the mesh version of the Sturm-Liouville operator.\n",
    "    \"\"\"\n",
    "    ### BEGIN SOLUTION\n",
    "    ham = np.zeros((npts-1, npts-1))\n",
    "    np.fill_diagonal(ham, 2.0)\n",
    "    np.fill_diagonal(ham[1:, :], -1)\n",
    "    np.fill_diagonal(ham[:, 1:], -1)\n",
    "    return ham\n",
    "    ### END SOLUTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-b55c38e4cbf73e1f",
     "locked": true,
     "points": 2,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "N = 20\n",
    "ham = hamiltonian(N)\n",
    "assert ham.shape == (N-1, N-1)\n",
    "\n",
    "### BEGIN HIDDEN TESTS\n",
    "from numpy.testing import assert_allclose\n",
    "N = 10\n",
    "ham = hamiltonian(N)\n",
    "assert_allclose(np.diag(ham),\n",
    "                2, atol=1e-15)\n",
    "assert_allclose(np.diag(ham, 1),\n",
    "                -1, atol=1e-15)\n",
    "assert_allclose(np.diag(ham, -1),\n",
    "                -1, atol=1e-15)\n",
    "\n",
    "for s in range(2, 8):\n",
    "    assert_allclose(np.diag(ham, s),\n",
    "                    0)\n",
    "    assert_allclose(np.diag(ham, -s),\n",
    "                    0)\n",
    "### END HIDDEN TESTS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-e3c64f6e9d4d3f08",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Solve the eigenvalue problem\n",
    "\n",
    "Write a function which returns the eigenvalues and eigenvectors of the discretized hamiltonian given the number of the mesh points,."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-c338c12b8b0ff65f",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def eigensystem(N):\n",
    "    \"\"\"Construct the Hamiltonian on a mesh and compute its eigenvalues and eigenvectors.\n",
    "    \n",
    "    Return eigenvalues in the sorted order.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    N : int\n",
    "        The number of nodes in the mesh: $x_j = a*j /N$, with $j=0, ..., N$.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    vals : ndarray, shape(N-1,)\n",
    "        Eigenvalues, in the ascending order.\n",
    "        \n",
    "    vects : ndarray, shape(N-1, N-1)\n",
    "        Eigenvectors arranged in columns: `vect[:, j]` is the \n",
    "        normalized eigenvector corresponding to `vals[j]`.\n",
    "    \"\"\"\n",
    "    ### BEGIN SOLUTION\n",
    "    vals, wavefunctions = np.linalg.eigh(hamiltonian(N))\n",
    "\n",
    "    # eigvals are not necessarily sorted, so sort them manually\n",
    "    # XXX: should also sort the wavefunctions\n",
    "    idx = np.argsort(vals)\n",
    "    vals = vals[idx]\n",
    "    \n",
    "    return vals, wavefunctions\n",
    "    ### END SOLUTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-79dc8ca40332b3cb",
     "locked": true,
     "points": 3,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "N = 10\n",
    "ham = hamiltonian(10)\n",
    "vals, vects = eigensystem(10)\n",
    "\n",
    "assert (vals[1:] >= vals[:-1]).all()\n",
    "\n",
    "for j in range(vals.shape[0]):\n",
    "    assert_allclose(ham @ vects[:, j],\n",
    "                    vals[j] * vects[:, j], atol=1e-10)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-af0f96f6ffa29419",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Long-wavelength limit of the FD problem\n",
    "\n",
    "We expect that in the limit of $N\\to\\infty$ solutions of the discretized system tends to the solution of the continuous problem. \n",
    "\n",
    "Plot the spectrum of the FD problem for several values of $N$. Note that the spectrum is bounded, i.e. we can only hope that low-lying eigenstates match: for large enough $N$, several first eigenstates of the FD problem are close to the corresponding eigenstates of the continuous problem. (See e.g. A.A. Samarski and A.V. Gulin, Sec 1.4.4 for details)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the spectrum of the FD problem for several values of N. \n",
    "# Plot the spectrum of the continuous problem on the same plot.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-79c6a06c54de6b6d",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Make the function in the cell below return a value of $N$ such that the first 10 eigenvalues of the FD problem and the continuous problem match within the relative tolerance of $10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-a0dccf755c881b94",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def get_large_N():\n",
    "    ### BEGIN SOLUTION\n",
    "    NN = 500\n",
    "    return NN \n",
    "    ### END SOLUTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-d85b44c5d92936e0",
     "locked": true,
     "points": 6,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "NN = get_large_N()\n",
    "\n",
    "vals, _ = eigensystem(NN)\n",
    "\n",
    "# the rest of this cell is a test which compares first 10 eigenvalues\n",
    "\n",
    "### BEGIN HIDDEN TESTS\n",
    "\n",
    "# Experimentally, NN=500 is enough\n",
    "continuous = np.pi**2 * np.arange(1, NN)**2 / NN**2\n",
    "assert_allclose(vals[:10],\n",
    "                continuous[:10], rtol=1e-3)\n",
    "### END HIDDEN TESTS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-4a4d614a4c80bf28",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Wave functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-fd83c17823e931ad",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Plot the wave functions (i.e. eigenvectors) for several lowest eigenstates. Compare to expected form. How many zeros does the eigenfunction number $k$ have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the wave functions here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-d0a5a2108321c342",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# Perturbation\n",
    "\n",
    "Consider now the modified potential: add to the square well potential the additional term \n",
    "\n",
    "$$\n",
    "V(x) = -\\frac{\\hbar^2}{2m}\\, \\gamma\\, x(x-a)\n",
    "$$\n",
    "\n",
    "for $0 < x < a$, where $\\gamma$ is a constant.\n",
    "\n",
    "An analytic solution is no longer possible, but the numeric solution is a straightforward extension of what we did in the previous section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-65d49ff5520f5952",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Numeric solution\n",
    "\n",
    "Again, define a grid function $\\psi_j$. The Schroedinger equation then becomes\n",
    "\n",
    "$$\n",
    "-\\left( \\psi_{j-1} + \\psi_{j+1} - 2\\psi_{j} \\right) - \\gamma h^2 x_j (x_j - a) \\psi_j = k^2 h^2 \\psi_j \\; , \\qquad j = 1, \\cdots, N-1\n",
    "$$\n",
    "\n",
    "so that the diagonal elements of the Hamiltonian matrix are no longer constants.\n",
    "\n",
    "Construct the hamiltonian matrix, solve for the eigenvalues and plot the ground state wave functions for a range of $\\gamma$. Discuss the structure of the wave function for $\\gamma > 500$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-f7dc2f611e4b19af",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def hamiltonian_perturbed(gamma, npts):\n",
    "    \"\"\"Construct the eigenvalue matrix for the discretized 1D Schroedinger equation.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    gamma : float\n",
    "        The extra potential is $gamma * x * (x-1)$\n",
    "    npts : int\n",
    "        The number of nodes in the mesh: $x_j = a*j /npts$, with $j=0, ..., npts$\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    ham : ndarray, shape(npts-1, npts-1)\n",
    "        The matrix of the mesh version of the Sturm-Liouville operator.\n",
    "    \"\"\"\n",
    "    ### BEGIN SOLUTION\n",
    "    ham = np.zeros((npts-1, npts-1))\n",
    "    \n",
    "    x = np.arange(1, npts) / npts\n",
    "    v = x * (x - 1) / npts**2\n",
    "    \n",
    "    np.fill_diagonal(ham, 2.0 - gamma*v)\n",
    "    np.fill_diagonal(ham[1:, :], -1)\n",
    "    np.fill_diagonal(ham[:, 1:], -1)\n",
    "    return ham\n",
    "    ### END SOLUTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-ad4764121a24f8dc",
     "locked": true,
     "points": 4,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "gamma, N = 590, 11\n",
    "h = hamiltonian_perturbed(gamma, N)\n",
    "\n",
    "assert h.shape == (N-1, N-1)\n",
    "\n",
    "### BEGIN HIDDEN TESTS\n",
    "x = np.arange(1, N) / N\n",
    "v = x * (x - 1) / N**2\n",
    "\n",
    "from numpy.testing import assert_allclose\n",
    "assert_allclose(np.diagonal(h),\n",
    "                2.0 - gamma*v, atol=1e-15)\n",
    "\n",
    "assert_allclose(np.diagonal(h, 1), -1)\n",
    "assert_allclose(np.diagonal(h, -1), -1)\n",
    "\n",
    "# check the rest is zeros\n",
    "np.fill_diagonal(h, 0)\n",
    "np.fill_diagonal(h[1:], 0)\n",
    "np.fill_diagonal(h[:, 1:], 0)\n",
    "assert_allclose(h, 0)\n",
    "### END HIDDEN TESTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-c6b6aa782bac0d62",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f44eea7fb20>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABVhElEQVR4nO2dd3xUVdqAn3OnJITee5UqFpCAIBZsiBV7XRUVcXd1dde1ILoqYl37tzYExAaiogJKB+mQSqJA6CUkQHohyZBMuef7485MJmECKTOZlPP488dtM/fcyb3nvW8XUkoUCoVC0XjRQj0AhUKhUIQWJQgUCoWikaMEgUKhUDRylCBQKBSKRo4SBAqFQtHIMYd6ANWhXbt2slevXqEehkKhUNQr4uPjs6SU7ctvD4ggEEKMAz4ETMBMKeWb5fZfDHwAnAPcKaWc77PPBWxzrx6WUt5wuvP16tWLuLi4QAxdoVAoGg1CiGR/22ssCIQQJuBj4EogFYgVQiySUib5HHYYmAA85ecrTkgph9R0HAqFQqGoHoHQCEYA+6SUBwCEEPOA8YBXEEgpD7n36QE4n0KhUCgCSCCcxV2BFJ/1VPe2yhIuhIgTQkQJIW6s6CAhxCT3cXGZmZnVHKpCoVAoylMXooZ6SikjgbuBD4QQZ/g7SEr5uZQyUkoZ2b79Sb4OhUKhUFSTQAiCI0B3n/Vu7m2VQkp5xP3vAWAtMDQAY1IoFApFJQmEIIgF+gkhegshrMCdwKLKfFAI0VoIEeZebgeMxse3oFAoFIrgU2NnsZTSKYR4DFiOET76hZRyhxDiFSBOSrlICDEc+AVoDVwvhJgqpRwMDAKmu53IGvBmuWgjhaLOEZ+cS9SBbEb2aQtQpeVhPVuHZtAKxSkQ9bEMdWRkpFR5BIpgUNEkX+JwsXZPJhYNZmw4iMMl0TQAga5LNOFelu5lUbpdCIFLl1hMGu/cdg7Nw83sOHqcUWe0K3MO32UlMBTBQAgR7/bJlqFeZhYrFIEkPjmXLfuzEMAHq/fidEmEAAmc6j3JpYNxFLhkuWV58rLdpfP4vESfb9iDcH/K5BYYupRYTRovXj+YXJtdCQVFraAEgaJR4Pum79J1fvvzGBFWE/szi1i1M/2kCb8iASAATRNIKTFpAgG4dGMZBC5dL7OsGaoCLl1i1jSG927D5n1ZeL7e86+vwCh26kz5xUi2N2uCe87vQfNwC5cO7AAorUEReJQgUDRYPJN/U6uJ15fswmG8wuM7x1tMosykb9IEuCd5hMDl0sssW8waL15X+rYOVfMRGOPKweHUTzqHxBAYAO5/cOqSr7YYVQE+WrMP9/CwmjVeUlqDIkAoQaBoUBiTfxY2u4vp6w7g1Cu27WgCbovsxs9bj+BwnjzJN82IJzfpd1qfeRmAd3lgpzQ4tAG0iwAYZq7CcvcRLLjB4vd7izoMI+pANq0jrLzy2w4cTt3rX/BchedySny0BotJ8NzVAznh0JVQUFQL5SxW1Hvik3PZtC+T7CI730Yd9r5V+2IShjXepAmGiL0MZwfxYjBPT7yvzIQ/sH0T2LsCNBNs/h+4HCDcUdZSB4RhH/IsGztOvSyEcbwwQc/RcHgz6C7jez37TBa4/GVwnIDeF7Er7Ti5Sb9T0HEkX2w6yDC5g1jOJFEOwKkbmo0/GWfWBE+O7Y+UKKGgOImKnMVKECjqJcbkn0VmQQlzopP9TorDTXsZIZKIF4N5cHRvmqdH0aVDe7rHvg66A4QJrcdIY2KWtVQGS5hAuipzICBBmNARxmdMVg6PeIljaUfKCIho/Uzi9H4nfYMSCoryKEGgqPd4zD4FJ5zM2HjwpDd/AUS6J/9i0YQplrlougMQGBZ4P/e6ZjLezk9ClGoCmsn4V3e5lwXozqovm6ww7k1YNhlcdmOflKWCQVYwRr9oSCGQUkfXrEx13EsLWUC0PJM418lCwWISfDFhOBFWs3I2N2JU+KiiXuJx+Jo1wdvLd3tt/ueJPYw07SRKH0QTzcl1YjMttSLGabEIjLf70lq35SZXYTK2maww7q2yE3P5SftENvQy7Pwc2lDz5e4joOOZJ+9r0tb/OLxCyOm+FM9F6QjpjmLSS3jFNNPYLUz84LyEfBnBShmJrsNIzfidJnxR6muwlvOHKKHQuFEagaLOEn8oh7tmRGN3R/ucJ/YwUttJsTQz2fI9ZozJUYhSq3wZhHsSPdUE330EpMRUPGnXJhWNw7N8krDAj4AoRQK6NLQhBxaeaTqNlNxir2BIpD9gCIU5E0cqYdAIUKYhRb1h874svtx8iKgD2fQtSeJCbRutxHHuM63GhI7wmfV9XbOggaYZJpZTvdHX9gQfSPwJC18BIYTb1GU8176/j0uzGs5vwIGZlx330UYUEqUPovNZF3Nml5ZKO2jgKEGgqNPEJ+eyemc6SUfzOb53M+NMMXQmm2tMMWiUZvp6MnFFZd/2GwseAVGhiUmD1r2QWXtLY5qk8Vs6MDPN8RdaiBPEicHceN2N5J5wKKHQAFGCQFFn+e3Po3w573tu0DbSnQwuMf3prt3jO/kLhMe52lDf9gPFqUxMX91gCAkkUuonCQUnZl5z3EVTYWerZoTXKmHQcFCCQFGn2BW7iiMJK9h/oikdMqO4wbTZmzULhoXDKQWaZkLDz+SvJvzq4Udz0AF0V6nw9REK37aYyMC2JtoMvoyBw68I4cAVgUAJAkXI2RW7ityk33HpkvMPTceMCyHAhUC4q3Z6Jn/jzd+CdvVbavIPFj5CQV/6rDt5jgqFwtYzn0UW5xmJd0oo1EuUIFCEBvdkk1xoonPUNCw4DXu/NN76XVKwu9N1DMxaaUxEavIPDX6EgkQipH6SULBjIW30VHo1KVZ/o3qGyiNQ1D4pMeizr0PoJfQAb7SPS4KOhpCGo1KLnIDW6XFl5w8l3Ud4f3fNnedw6EQ4HTe9jEUaIaomDKEQjoOem6cYGc+mMLQJv6q/WT1HCQJF4EneAnGzKNm7jjC9BDDq4rjcnVEdmPnzrMnotuyyZgY1mdQN3EKhF7Cr9QByk35Hi2jLOdvfxCKdCCQaEk1IpLOYol/+SdNBV8LAa9XfsJ6iTEOKwJASA/tWQsZu5M6FhqnZ/eYPxuS/tNs/6WwtUjbmeorHx3O0JIJrjnyIFQcgvaYjhIa48F9gbao0uzqKMg0pgsfupfD9X3yyXDHC19H43nUpR2lnVPocq0IR6zMDh18Bw68gPjmXB2Y2Z5jcQWeyuNP0O2YhkboOG941DjZZ4f7foMf5oR20olIoQaCoPjsWwKb/g2MJ3ph0pwSJCSElDswcH3gbEV2H87RKTmowDOvZmqcn3kfUgWwKj8Ti2LsBpNNouelxLrvs8MN9MGIiIKD3xUpDqMMoQaCoGikxkLQIUmOQKdEA6Ahc0oyGjgMzr7ruoxUFhhZw0TglABogw3q2ZljP1sQnt+WB3S8wTO4gj+a8YPoai3SiI9DtdsJ/f9Wo+6pZ0K7+LxTnKrNRHUT5CBSVZ+di+OFekC5vY3cj9l9jnutSjsp2xMhB9B9+BV1bNVElChoJngqxR/NOsCt2FecLo6jdBaYk/mX6EZPwrXukgTkM7l+khEEIUD4CRfXZuwo2vA0pMV4TkI5AR0N3m4AWyovZqvfDYtZ47rxuSgA0Ikq1g1zu2TqQRGd/zCaN9hFh2EsWlIafCh3Qkc4TiA3vQffhSjuoIyiNQOGflBjYvwbStsGuXwGQaDik5jUBveK8lzaikHgxmBuuu0nVtld4tYORfdoC8OaMrxguk8iRzXjJ8g1WHEb4qSfSyBwG4/6rEghrCZVZrKg8h6Phq2u9JYs9hd+caMxzGiagaDmIAcoEpDgN8cm5rN+TydLtx2iWsZWR2k66iwzuMK9FQ7rvLXdfZ5MyGQWboJqGhBDjgA8BEzBTSvlmuf0XAx8A5wB3Sinn++y7H3jBvfqqlPKrQIxJUUU8JQaadYa1b/gIAYFTau7mJmYW6BeTIA0T0BRlAlKcBo/Z6OL+7blrho2tzv6cJ/Zwo2mTUW7EnZgmpARnMSTMURnmIaDGGoEwCsPvAa4EUoFY4C4pZZLPMb2AFsBTwCKPIBBCtAHigEiMF894YJiUMvdU51QaQYBJiYGvrgdnCSDBEmHUmpE6JdLEy3bDBKS0AEVNiE/OZePeTFbsSCcsLY6R2k5yacY06zeYpd19lDD+Vw7loBBMjWAEsE9KecB9onnAeMArCKSUh9z7yvfTuwpYKaXMce9fCYwDvgvAuBSnIyUGDm6AA2uNtzEABIx6lIW2wRyIWcZGx0D+EAOQulRagKJGeLSDC/u1564ZhWx1Gq0ytdaDuabZPs5ybqfl0fWAWzuI+0JpB7VEIARBVyDFZz0VqGw6ob/Pdg3AmBSnIyUGvrrOrQUAbjutNIXx1r7ufHbACtyA1azxyvWqybkicAzr2ZrvHh7Jxr2ZxB7KYc4+mEMnRpg7MS8sBs3l1kz/+A5DOwhX2kGQqTfho0KIScAkgB49eoR4NPWYlBg4sA72LCsjBDL730ViQXN+zOrFigPtvC0hXS6dXJudRy/tG8JBKxoaHu3g4zX72Lw/G11CjLMvz7SZxs1tDjGQZNoc+g1DOzgB0Z8p7SCIBEIQHAG6+6x3c2+r7GfHlPvsWn8HSik/Bz4Hw0dQ1UEqOFkLEBogcGkW/r5jALGufgA8dGEv5kQfxuHUsZg1byigQhFoRvZpi9Ws4XDqSAnzM7oyP6Mr55s78l3YKjSXHdBh+0+w/WelHQSJQAiCWKCfEKI3xsR+J3B3JT+7HHhdCOGxN4wFngvAmBS+eLWApT5agEZmvztJsrXk62PdiXX1AsAkoE3TMOZMHOmNB1fmIEWwGNaztfdeO5p3grnRh5FAtLMvUzu9wbjm++gvjtL2wC8o7SB4BCSPQAhxDUZ4qAn4Qkr5mhDiFSBOSrlICDEc+AVoDRQDaVLKwe7PPghMcX/Va1LK2ac7n4oaqgIVaAG6ZuGOkueIdRpagFkTSGk4hOdMHKkmf0WtE5+cyz0zo7A7dXSfaWmkeR9zw14v1Q4A5TuoHiqhrDEiJXx/D+xa7N6gkdn/TvbaW/NjVi9+yTL88pqAO0f0UGGhipDjyUxOzbUxLyYFz+z0aN9sLrTs8tEO3Fz4JFzxUkjGWh9RtYYaG3tWwIoXIGt3GS3gsaSBRDsNx68mjKhti1njFhUWqqgD+NYt+iXhiFc7+HhfWz5hNOdb9jHXurQ0sih6Ohw/CsMfUppBDVAaQUMk+nNY+rSxrFnAXf73w/0deX+XMdkrLUBR1/FoB/szC/l5a2n8yVsjTnBH+2RwOmDdm4AEzQT3LYJeF4ZuwPUApRE0BpK3GOUhDq4r3SZ1Uo4d4aljVxB9MMfI2xRgVVqAoo7jqx0s2XaMEoeOBN7d2ZI91qt5UP5MV6GBdIHugp8egqvegNyDypFcRZQgaCjsWQ7f3QlSBzQwmUF34dIs/DOqGfF6DiZN8OJ1gygscSktQFFv8I0syiosZvamZGZtPMif5mZ8H2ZB0zE0ghP5MP8BlCO56ihBUN9JiYHEOfDH924hAAhBZt9b2ZrfnC+PdCNeP8PYLiWFJS6VHKaod/gmoGkCdAmxzn78u9U0bmt3iDaDL2Ng3gbY9AHeEhUJ36ow00qiBEF95nA0fHkt6EalUEwW0HV0zcLft5cmiFlMAt1dK0glhynqM2US0IBfsrryS1ZXwvY5WDj+IgaaPystnrj1a1XeupIoQVBfsdtg8ZOlQkCYYOhfcDbvxuT4lsQWdQGMBLHbIrsrp7CiQVBRAlqJU2fOkU5Mu/9XQws4uN4opigluEqMbUoQVIgSBPWRpIWw5GkoTAfNbNzsJisxLcbxxCYLx/KLyySIKaewoiHh60T+aWsqdnd5iu9iD9Mq4gzCLeO5fNC5DDy8xdAOpA4ZSbD+XeitzET+UOGj9Q3f0FCTFa5+G05ks8LWj0lrTIBhCpp6w1mqYqiiweMJMT2zc3PeXbGH7UePI4Awi8aCGywMLIiGnb9C5k7jA+ZwuP/XRisMVPhofSclxnCE7VpSuk13ceRYKq/mX83yHWmlm3WpKoYqGgUe7QBg+9HjbD96HAkUO3RmHmpH73Z3cGN3na6Zu/A6keNnKydyOZQgqA8kbzF6COsuQBiagDs09ImoZsS50hCA1SRwKaewopFywRnt+Niyz5tvMD/+CJqADeYWzLWGldYqSpxrZNsrJ7IXJQjqOiUFsOgxtxDAuIGH3kNx0678O6Y5cUVGBXBNwK3KKaxoxPg6kv9IyWNFUrq3z8HPwz7l1rYH4chW2PWb4TdQTmQvShDUZXYuhsX/gsIMo1SE1MFkZVPTK/nHRit5RfYyoaHKKaxo7Pg6ktftyaTEXatobVFP0luezeVnnMvAfasME5HUIXMPbHi30ZuJlLO4rpI4DxY8Yiz7OIWXFfbjr+sMp7DVpPHyDaqNpELhj/jkXDbsyWT5jjR2phUAEO5xIhdGw46FRlHGRpSJrJzF9YkD6+DXx0vX3U7h145fzbLtpU5hl67aSCoUFeHRDswmwa60Aq8TeVZyO3q1vYMbe7jomrUbrxP5wLoGLwgqQgmCukRKDER9auQJtOxm5Am4HN56QbFep7CGS1dtJBWKyjDqjHaE+TqR444gBGwwt2SuNby0pPW2+aA7oe/ljU4gKNNQXSElBmZfY2QKCw3umQ9hzXEeWM/TsS28TWRMAu5Q5aMViirhyTfYejiX1TszACPA4r/nlxhO5Mw98Oc84+AGnGugTEN1GSlh1dTSchEIOJbI2g738sKWQlLzTqhMYYWiBvg6kTfuzfI6kRNkP9Kd53CjZR5d0QDdMBPt/K1BCoKKUIIg1CRvgeVT4OhWo14QGJFBzoE8MDsWicoUVigCxbCerZn78Eh+35XOosSjzIk+jAA2WJoz12pF0+1GNFH8bHDZ4aybG4VAUIIglBzaZDSWl7pRM+iad+BEDnuaDOHhhQ4kRu6AyhRWKAKHRzuwmjTeX7UXCUQ7+vLzsE+4te0hcBTD+v9C9KcQ9wVMaPjagRbqATRaHCeMyCBPDwEpOXLsCP86ejnX/VJC83AzYWYNk0A5hRWKIHBhv/aEWzQEIIEvktvzYcn1HCmS7j7fGElnCd+Gcpi1gtIIQsH+tfDrE5B3yJsopnsjg4yIhjduOZuWTaxEHchW5iCFIgj4ZiInHc1n8bY0ko4VsNncnO/CrGi6w3hRS5wDLidETmiwmoGKGqpt9q2GObe4S0db4Op3kLZsXvyzFd+kdgKMyKAnxw5QpiCFopb4eM0+3lm+G89sOHVoEfd3SQHNCiv/Axil3pmwuF4Lg4qihpRpqDax5cDCRw0hAKDrpBxN4frEEXyT2gmTJpQpSKEIASP7tCXMoqEJY/3tpBa8mHMVR3LyfcxEdoieHrpBBhFlGqot9qwwiscVZZWpHvrPqGZs149j1gSv3DCY3BMOZQpSKGoZXzOR3anz4eq9fL0lmV3m5swLs6DpABK2zzfMRSP/Vq81g/IoQVAb7FkOc+/Ao14eOv9lDqem8kVqV+L13gBIKck94VDmIIUiRHiiiT5esw9N4K1cOrXTG4xrvo8uHdrTM+YV2PEz7PoVJixpMMIgIKYhIcQ4IcRuIcQ+IcRkP/vDhBDfu/dHCyF6ubf3EkKcEEIkuv//LBDjqVMUZsLCx8BtfZS6i582/Ml9ey5kra03ZmUOUijqFCP7tMXqjtgD+Cq1I3fvHM1PW3Z5rbq4HLDl45CNMdDUWCMQQpiAj4ErgVQgVgixSEqZ5HPYQ0CulLKvEOJO4C3gDve+/VLKITUdR51kz3JDCJzI8ZqDnJjZ5BwIGCnutw9XPQQUirqEr5nocE4R38emIoGNjoH8o4kFi3QCEpIWwA/3w6hH671mEAjT0Ahgn5TyAIAQYh4wHvAVBOOBl93L84GPhBAiAOeuu+xZeZI5aPeBZL5I7cpW2RdNgFWVi1Ao6iS+JSkWJh6lxKGzVfbn5VZvcF2L/XTp6DYTJS2A3YvrvZkoEIKgK5Dis54KnF/RMVJKpxAiH/DYQXoLIRKA48ALUsoN/k4ihJgETALo0aNHAIYdRGw5sOhRypqD/uB/jvEI4Omx/UEIpQUoFHWcUu0gi5iDOczZA3OOduLxA4v4lwYCDDNR9GeNXhDUhGNADyllthBiGLBACDFYSnm8/IFSys+Bz8HII6jlcVaefauNEFGf6CAHZjY5BwGGOQghlFO4mjgcDlJTUykuLg71UOol4eHhdOvWDYvFEuqh1Bs82gHsY/2eLP9mou2/QMezjOV62O0sEILgCNDdZ72be5u/Y1KFEGagJZAtjWy2EgApZbwQYj/QH6if2WL71/okixnmoPid+5iT1oNE+iuncABITU2lefPm9OrVi4ZuXQw0Ukqys7NJTU2ld+/eoR5OvcOTa+AxEz1qfpm/dEqh46ALGLD1NVg9FdDAHFbvup0FQhDEAv2EEL0xJvw7gbvLHbMIuB/YAtwK/C6llEKI9kCOlNIlhOgD9AMOBGBMtY/jhFE2Qpaag+av/4OPnOPRBLx0/ZkUlriUOaiGFBcXKyFQTYQQtG3blszMzFAPpV7i60ROyy/mmyhYcbwnYQc01gy7ki6ZOwHdqE90aEO9EgQ1Dh+VUjqBx4DlwE7gBynlDiHEK0KIG9yHzQLaCiH2AU8CnhDTi4E/hRCJGE7kv0opc2o6plrn0Cb49EKjdpDJghQm7NLEZpdhDhJAYYmLRy/tq4RAAFBCoPqo365mDOvZmkcv7UunluHeLOQSp84K+1lGQxswEs6K82HDu0bDqXpAQHwEUsolwJJy2170WS4GbvPzuZ+AnwIxhpCRvMWnlLSFg8NfZv0fu1mY14c/xABMSGUOUigaGJ5cA7u7wc37u1rTa/Qsmh1exTnpi7Bu+tAoTWGqH2aiUDuL6ze6Dsue9ZaSllJn/sY/+dg5HrMmmKZKRigUDRJfM1G3Vk2Y+msSE1YJNHElky1ZPKz9ipC6UZ+oHpiJVNG56nI4GmaPg2N/gGZGChMl0sQWtznIt2SEEgINj23bttGzZ08+/fTTGn3PsmXLGDBgAH379uXNN98M0OgUtYHHTDR+aFdujewGGGUpljuH4RJW4yCpG1HkddxMpARBdUiJgS+vgZRo0MzsGfYy0013cnfJFLZpA1R0UCPg7LPPZt68eXz99dfV/g6Xy8Wjjz7K0qVLSUpK4rvvviMpKen0H1TUOa4a3IkwszGdxuv9+aLv//ijyx04TVb4fRr8/ip8dUOdFQbKNFQdNrwHuhMAXUoWbNnGJ87xqrdwHSU+OTcoDX46dOjAjh07qv35mJgY+vbtS58+fQC48847WbhwIWeeeWaghqioJTy9kNfsyuCXhFRe3waC8bxhzeIObTVCyjptJlKCoKpsmw97loLQkAjs0kSU2xykegvXLlN/3UHS0ZNyD8tQUOxgV1oBujSS+QZ2ak7z8IqTqc7s0oKXrh9cqfNPnjyZkpISkpOT6dmzZ5l9F110EQUFBSd95p133uGKK64A4MiRI3TvXpqC061bN6Kjoyt1bkXdw5N4ZtIEH642eiH/4LiI28LXY5IOQEK3uicEQAmCqrH5Y1j5AnQ6m+j+TxG9YTkbnQPYrg3ApOvKHFQHOV7sRPf0AZLG+qkEQWVZunQpRUVFXHvttezYseMkQbBhg99KKYpGwMX92zN93X6KnUbi2ftd32N8ya/0zVyBWPMaHL4c+lxSpzQDJQgqS/xXsGIKAK6MPbyVfJCt8jqsZo2XbxiszEEhoDJv7vHJudwzMwqH0xDUH945tMZ/o+LiYp599lkWLVrE7Nmz2b59O9dcc02ZYyqjEXTt2pWUlNIyXampqXTt2rVGY1OEnmE9WzPn4ZFs2pfJsu1pfLQPPmYC/7GaefDwEji8xXAe16GwUiUIKkPWXlj6rHdVuhyM1Hay1dUfl0tX5qA6jG+YX6AE9auvvsp9991Hr169OPvss1m0aNFJx1RGIxg+fDh79+7l4MGDdO3alXnz5jF37twaj08RejxmIiEESccKkECWqzlSEwgkOIvrlL9ACYJTkRIDu5dAwhwwh+HSXUiXEwdm4sVgFR1UTygtGlZzdu/ezcqVK9m0aRNgRA+9/vrr1fous9nMRx99xFVXXYXL5eLBBx9k8ODK+ScU9YMLzmjHx5Z9lDh0ovRB2LFipcQQBmk7DM2gDhSpE1LW3UKeFREZGSnj4oJcly4lBr663pDcwJazX+HtOJ2R2k7ixWDGX3+TMgeFgJ07dzJo0KBQD6Neo37D2sWIWsti895sThzcwihtB9eaYxgkkhEIozRFLZmJhBDxUsrI8tuVRlARB9Z6hYBEY1PiDrbK8Wx1GVVEL1bmIIVCUQl8NdJ3DvRnq6s/GpJBlmRA1okidSqhzB+6bvxhAB1BsTSzt8kQwtx9TJU5SKFQVJWRfdoRZtEQwCb9LOxYkRilaWgR2iABpRGUJyUGVvwHUqL42jWWNL01Ufog/nHXLUxqYg1KYpJCoWj4+AYu7Djaibu2w5VaHHeY19F8yWQs6Ukw6LqQaAZKEPiSEgOzrwHdgQsTC5wXsFUapqCkYwWqbpBCoagRHjPRx2v2sXSbYSbaLbvzvvgUNn8IMdPh/l9rXRgo05Av0Z+B7gCMonEjTTuVKUihUAQcT7czAXQROehG92Oks9jwT9YySiMAQxNI+BZ9+wKkFEgEDswMHnUNT0YMUKYghUIRUHzNRDLlYuz7f8GKA5OQFP25iKZCg94X15pmoARBSozRWMZZAgheddxDE+EgRg7i0oizVGSQQqEICqVmIrhn5xRGajs5S+znmuw4+H17rYaVKtPQ7qVIZwkAuhREaA6m6+PZbhqozEEKhSLojOzTliTzQD7Tx7Nd9sHd5sowEx2qnZpVjVsQlBRy4s8FIMEpDXPQWaOv5cmxA5gzcaQyBykqJFCNaR588EE6dOjAWWedVWb7qRrWqGY2DQuPmejfYwfQevDllEirUShRSlLzHbUyhsYrCJK3wOdjCDt+kNedd/Oe83budUzhQPhgFR2kOC2BaEwDMGHCBJYtW1Zm26ka1qhmNg0TT7cze+dI/uKYwgfOW9gnu9Ax/l1Y/O+gN7RpnILgcDTyy2shey9OqZEgByhzUEMmJSYorQJr2pgG4OKLL6ZNmzZltvk2rLFard6GNafbp6j/jOzTlh2mgXyk38IbjrsxSTsydib6l9cFVRg0Smdx/vLXaCldAGhInh+czeYu16nooPrG0smQtu3Ux5Qch/TtRu9YoUHHsyCsRcXHdzobrq6cuaWmjWkq4lQNa1Qzm4aNbzTR4P0bkCkCTUhwlnA07le6BMlx3LgEQUoMMno6LY+swyU1JODAzP6mQ1V0UEOlON8QAmD8W5x/akFQSVRjGkWw8EQTzc8fgT1lFlbpQENi3fULrGsOfcYEPJKo8QiClBj02dcidDsuKZjqup/m2IgXg3l66KWhHp2iOlTmzT0lxmga7rKDyQq3zKzxQxSoxjQVcaqGNaqZTeOh99BLeSD+BYbJHXQnjTtL1iHXvIZc/y7ahMBmHzcaQXAkYTldXHaEAIng5kHN2dTlMZ5W5qCGTfcRRiz2oQ0Bq/seqMY0FXGqhjWqmU3jYVjP1jw98T6iDmTT9eAs9OR1aAJ0ZwlHElfQNYCCICDOYiHEOCHEbiHEPiHEZD/7w4QQ37v3Rwshevnse869fbcQ4qpAjMcfW/TBFGPBKTUcmNnnNgcpIdAI6D4CLvp3QISApzHNP//5T8CIHtq+fXu1v++uu+5i1KhR7N69m27dujFr1qwyDWsGDRrE7bff7m1Yc6p9ioaHJ5oovc1wSrDilAI7Fra4zgzoeWrcmEYIYQL2AFcCqUAscJeUMsnnmL8D50gp/yqEuBO4SUp5hxDiTOA7YATQBVgF9JfS7cmtgOo0polPzuXtmV8zTO4wzEET71NCoB6imqrUHPUb1j8CNX8FszHNCGCflPKA+0TzgPGAb3DzeOBl9/J84CMhhHBvnyelLAEOCiH2ub9vSwDGVQZfNUuZgxQKRX0i2PNXIARBVyDFZz0VOL+iY6SUTiFEPtDWvT2q3Gf9er6EEJOASQA9evSo1kAD2btWoVAoapNgzl/1JqFMSvm5lDJSShnZvn37UA9HoVAoGgyBEARHgO4+693c2/weI4QwAy2B7Ep+VqFQKBRBJBCCIBboJ4ToLYSwAncC5ePpFgH3u5dvBX6Xhpd6EXCnO6qoN9APCG5RDYVCoVCUocY+ArfN/zFgOWACvpBS7hBCvALESSkXAbOAb9zO4BwMYYH7uB8wHMtO4NHTRQwpFAqFIrAEJKFMSrkEWFJu24s+y8XAbRV89jXgtUCM43RsTd/K+tT1jOk+hiEdhtTGKRUKhSIgBHP+ajSZxYkZiTyw7AF0dObsnMOMsTOUMFAoFPWCxIxEJq6YiEN38E3SN8y6alZA5696EzVUU+LS45AYyXN2l5249KolpCkUvgSqMY2nTMWQIUOIjCzN81GNaRS+xKbF4tCNJjUu6Qr4/NVoBEFkx0isJquxIox1haK6BKoxDcCaNWtITEzEky2vGtMoytPU0hQAgcCiWQI+fzUaQTCkwxBmdXmGf2zryhkpLiIsEaEekqKWSMxIZOa2mSRmJAb0ewPRmMYfqjGNojxrU9YyJL0JbxwczswuTysfQXWxJSQQ9q/XuchhZ4Qm+bXru/x7wvRQD0tRA96KeYtdObtOeUyhvZDdubuRSASCAa0H0MzarMLjB7YZyLMjnq3U+QPRmEYIwdixYxFC8MgjjzBp0iTVmEZRhqTsJLJiN/PKPIHJuQWxIB7b7H5EDB0asHM0HkEQE4ssKQHAgiB/yybeO+s9LutxmXIaN2AKHAVe35BEUuAoOKUgqCyBakyzceNGunbtSkZGBldeeSUDBw6s8dgUDYfEjERe2fIK56aaMDkdICXS7sAWE6sEQXWIGDEczGZwOhFmM9t76OzdMZvvdn2nIojqKZV5c0/MSOThFQ/j0B1YNAtvXvRmjf/WgWxM42kq06FDB2666SZiYmIYPXq0akyjIDEjkYdWPITdZcfUXXCbEAgpEVarMZ8FkMYjCIYOpfO0aRx77jlSrj+Pvd3iAbDrRgSREgQNkyEdhjBj7Azi0uOI7BgZkL9zoBrTFBUVoes6zZs3p6ioiBUrVvDiiy+qxjQKwIh0tLvsAOztAq5wM017nkGnF/8TUG0AGpEgAGg5/gYy3nmHrvlmrJoVu25HIFQEUQNnSIchARP0nsY0mzZtAozooddff71a35Wens5NN90EgNPp5O6772bcuHEA3uYzLpeLBx980G9jmvL7FA2Lvq1K+6j3yzRjtp2g7QMTAi4EoJEJAqFpNB19AWzcxKwXZvBy1CscKzxW5gdXKE7FgAEDyjhnBwwYwNatW6v1XX369OGPP/7wu++aa645ydxUmX2KhkN8ejwCwX1n3sdV2YXA9zS94IKgnKvRhI96aDZ6NK6cHAbmhPPqha9ic9n4Yc8PoR6WQqFQeMkvyeeH3T8wrvc4nhr+FC0SDxI2aBDmdu2Ccr5GJwgiRo0CoHDTJs5qdxajOo/i6x1fU+wsDvHIFAqFwmDurrnYnDYeOush9KIibAkJNBsdHG0AGqEgsHToQNiAARRtNGy8E8+eSHZxNs+seybgCUcKhUJRVaKPRjNr2yyGtB/CgDYDKIqJAYeDpqNHB+2cjU4QADS9cDS2rVvRbTasJisCwZrUNUxcMVEJA4VCETISMxL566q/UuIqISk7icSMRIo2bUaEh9Nk2LCgnbdRCoJmo0eDw4EtNrZM8SZVjE6hUISSqGNROKUTKC0uV7RpExEjhqNZrUE7b6MUBE2GDQOLhcxPP2N4ZgvCTGHefUM7BD40S6FQKCpDgd1IQtTQsGgWhu/WsR88iLVnr6Cet1EKguKkJHC5KE5MJOzJN5jZ5Wmu7X0tEsmxomOhHp5CoWiE2F12lh9aTv/W/Xls6GPM7PI01pf/B0De999jS0gI2rkbpSCwxcSCdNefsdvpujeP1y96nb6t+jLjzxm4dNUtU6FQ1C4L9y8k3ZbOv4f9m4fPeZiue/PAaZiJpMtlzFtBolEKgogRwxEWi7GiaYb9TWg8cu4jHMg/wHMbn1NOY8UpCVRjmgcffJAOHTpw1llnldlencY0qmFN/SUuPY4P4j+gT8s+jOpihLhHnOc2UwuBsFgCXl/Il8YpCIYOpceXX6K1aEGTc8/1pmx3aNIBgWDpwaUqgkhxSgLVmGbChAksW7aszLbqNKZRDWvqL57CiMftx0kpSOGPTCPbXLqtFi2uu5Yes78ISmkJD41SEIAhbVtcew3FO3ei243CTlszSksFqAiihoMtIYGs6Z8H3MYaiMY0F198MW3atCmzrTqNaVTDmvpL9LFonLphAtKl7p13ijZsALOZTi+9FFQhAI2s1lB5ml10MXnfzeNEfDxNR43ytrMscRl9C87rcF6IR6g4FWmvv07JzlM3pnEVFlKya5fhExKCsIEDMTWruB9B2KCBdJoypVLnD0RjGn9UpzGNalhTfznhPAGURgp5imAWrt9AxHnnnfJ+DRSNWhA0HXk+wmKhcP0Gmo4axZAOQ5g5diZzd85l6aGlpBWlhXqIihqiHz/uDQxASvTjxwPyYAWqMY2iceNwOVhycAl9Wvbhuj7XMbzTcIZ0GIIjPZ2S3bvp8NS/a2UcjVoQaBERRAyPpHDDejo++wxglCw+p/057M3by6d/fMrYXmMxa436Z6qzVObN3ZaQwOEHHkQ6HAiLhS7vvF1jNTuQjWn8carmMxXtUw1r6ie/7PuFY0XH+OyKzxjdtbSERJH7RaLpRRfXyjga/QzX9KKLyXjrLRxHj2Lp0gUATWg8OuRR/rX2Xzy7/lnuPfNe1bimnhIxdCg9Zn9htPYbMbzaQsBls6EXFaE1bcorU6dyz8030719ewb37cuiRYtw5uUZGseJE2hNmrDm11+RxcVozZqBpnk/a4qIKPNdpoiIk85VncY0AwYMUA1r6hmxabG8F/8e/Vr144IuZQvKFa7fgLljR8L696uVsdRIEAgh2gDfA72AQ8DtUspcP8fdD7zgXn1VSvmVe/taoDNwwr1vrJQyoyZjqirNLr6IjLfeonDDRlrfcbt3e9smbREIViSvYF3qOmaOnamEQT0lYujQagkAl82GXlgIusSZnQVSsufgQVYuXcrvX3+N/eBBBrRqxfYdO3CkppZ+zvdLsrLKfKcWEYF+4oTXZ/HAiy+ybuNGsnJy6NatG1OnTuWhhx6qVmMa1bCm/pCYkciklZNw6k4OHT/EH5l/eOcX6XBQtHkzLa4ehxCiVsZTU41gMrBaSvmmEGKye71MI1m3sHgJiAQkEC+EWOQjMO6RUoYsPMfapw+mdu3I+fprwvr3804Y8enx3mM8EURKEDR8XAUFuPLykE4nelHRSfv79+7Nep837f69e7Plh8r3s/AKAQApmT11qrEshFFrXghcNlu1GtOohjX1h81HN58UKeSZX/Lmz0cvLMTcrVutjaem4aPjga/cy18BN/o55ipgpZQyxz35rwTG1fC8AeNEYiKu3Fzs+/dz+IEHvCGGkR0jvTWIJJKz250dymEqgoSUEmduLiUHD1K8ezf25GRc+fl+hIAw/hfu/6u1rGHp3Nm9ftJAcGZm4szIwH7gAPYjR3CkZ+Cy2YL8CyhCQYbNMHyUjxSyJSSQ9uprAGR/8mlQy0r4UlONoKOU0lOcJw3o6OeYrkCKz3qqe5uH2UIIF/AThtlI+juREGISMAmgR48eNRx2KbaYWNB1AKTdYdiShw71Nj1fuG8h8/fOZ3vWds7vfH7AzqsIHS6bDT0/H6nr6IWFSIfD/4FCGDqsEFg6d0K6XGhNmwJ4bfxVXTZFRCDCw9GLihAmE45jaaUaAqW3vivXrTBnZmBq0wZhNqM1a+bXp6CoX+QW57L80HKGdRjGhd0uJLJjpFcbsMXEgsswLkqn0zsfBZvTCgIhxCqgk59dz/uuSCmlEMLvJH4K7pFSHhFCNMcQBPcCflM1pZSfA58DREZGVvU8FRIxYjjCakWWlIAQZdK4PU3P023pzNo+i1v730rLsJaBOrWiFvHY+6XDUTrJAsJc7hGoYPIvPwH7rldn2bPuVygIfIQDuHJyjIWMTMwdOwBU6GhW1H1mbpuJzWnjP6P+wxmtziizz9q7t7FQC2UlfDmtIJBSVhjrJoRIF0J0llIeE0J0Bvw5eo8AY3zWuwFr3d99xP1vgRBiLjCCCgRBsDDKTczm2JTncRUU0GTIkJOOeeK8J7jt19t4I+YN+rbqW0aCK2ofKWWlnGie6BwhBI709DKTqweteQtc7oif003+waBqQkHiTE83FoXA0rlzlcdagcKtqCVWJ6/m253fcmGXC08SAoD379vm/vtpftXYWtEGoOamoUXA/cCb7n/95bQvB14XQrR2r48FnhNCmIFWUsosIYQFuA5YVcPxVIuIoUNpM2ECaS+9RMnevYT3719m/4A2A7igywUsPrAYDQ2rycqMsTOUMAgB4eHhZGdn07Zt21MKA5fNhv3gQb+Tv+9bv6l1K0ytW50ynLO2OK1QALzmIylxHD1qLAsNa+9epx27lJLs7GzCw8ODdAWKU5GYkciT655ElzrRadEkZiSeNIcUrl2LtXdvOk5+1v+XBImaCoI3gR+EEA8BycDtAEKISOCvUsqJUsocIcQ0wFND9RX3tqbAcrcQMGEIgRk1HE+1aTZmDACFa9aeJAgA+rXux6ajm9DRcegOFUUUIrp160ZqaiqZmZl+9+t2O7K4GL242FvCtwxCYGrRAqnriLAwtOTk0n3lQj3rCrrTgSwpQWgaLt9MaR9Ebg7CajWu6RSdrMLDw+lWi9EoilKWHlyKLg1/pFN3njSHuAqLsMXE0Pree2t9bDUSBFLKbOByP9vjgIk+618AX5Q7pggIXhPOKmLp2IHwwYMpXLOGdo9MOmn/5T0u55ukb3BJFyZh8nr5FbWLxWKht8eOWo6i6GgOPzQRnE5MAJoRFCesVjo+NxlXXn6NksrqAraEBGwxsZhatST9jTcN35aPYBBWKz2++rJeX2NDREpJQoYRAaSJspFCHoo2bUI6HDS/dEytj6/RZxb70uzSS8n6+GOc2dmY27Yts29IhyF8esWnPLHmCbo378657c8N0SgV5bFtTSB3zhwK1qwp1QI0jVa33YalS5d6P/n74pscF9a/P7aYWOypKeTP/wmkRNrtHHvhBZpfeSXNLrmkwVx3fWdNyhp25uzk/jPvp1V4K79+xsI1a9BatqRJCP5mjbYMtT+ajRkDUlK4br3f/aO6jOKpyKfYk7uHVYdD4s5QlOP4kqUk/+UvHF+8GGmzgdkMJhPCaqXljeNp98ikBjsZRgwdSrtHJtHq5psRYWFg0kDTsO8/QPZn00m+9z5sW7ee/osUQcXhcvBu3Luc0fIM/jnsn0w8e+JJQkC6XBSuW0eziy46OZKtFlAagQ/hg8/E1Lo12V98gbV3L78TyM39bua7Xd/xRvQb7M/bz8jOI5WvIATYYmPJ+vRTiqKivXkgmEy0uuWWBqcFnA7fekqOo0fJ+/FH4zdxOkl9/AlaXHsNLcaNazS/R13jnbh3OFxwmKcjn66wgGXej/Nx5eZi7d2rdgfnRtTHcLLIyEgZFxf4qhS2hASS/3IvuFyI8DB6zJ7t9+H5Jukb/hv7XwSCMFOYiiCqZfJ+WcCxKVO8IZ/CbDacvxZL0Ds51XV8q60C3uQkzGZ6fvUlEcPqjFuuUbA+ZT2P/v4oAOGmcL9zRZl5JyyMHl/6n3cCgRAiXkp5koNTmYZ88Jdl7A9P4xqJVJ3MahFbXByHH3yoVAgAaBotb76Z9o8/3uiFAJRqB+0ff5xWt97qdZjjdHLkqac5vnRpULq1KfzzUeJH3mVPtGF5/GUT1zbKNORDxIjhiDArsvjkLGNfPHWISlwlIFARRLXA8RUrOfLEE6VagMXi1QJa3ji+0QsAXzwOZVtCAvkLFxragabhzMnhyL+eNH6/sDAlOINMQkYCO3N2YhImAL+RQgCWrkb5+9rOJvZFCQIfjLep2aS9PBXH0aM0Odt/oTlPJ7NPEj9hy7EtFDlOrlKpCAy2hARy587l+LLlJ2kBjc0XUFXK92Io+H0NOTNmGNFFxcXkzf+pxn0aFP5x6k5ejXqVTk078eoFr7Ite1uFFQkc7oZCbSZOpPlll4bkb6F8BH44vnwFR554gh5ffknTkRUXmrO77Nyy6BaKncXc0v8W5TgOMLbYWJInPOBVm321APU2W3W8/gPf3ANNM3IP1O8ZUN6KeYtvd37LE0OfYOI5E0957IGbb0YLC6fXd8FvJKR8BFWg2UUXIsLCKFi58pTHWU1Wbh9wO2m2ND5J/ISHVzxMYkZi7QyygeNIz+DIM8+WOjtNJuULqCFe/8E//0mzyy8zNuo6sqQkJHbphsqaw2v4due3AEz/c/op5wR7SgolSTtpfuWVtTQ6/yhB4ActIoKmF11IwapVSE9oYgUox3HgyZ03j/3XXI0zKwthsRh5AW5fQEPOC6gNPLkHbSdONHIPAKSkaMtmsj6tvfr3DRUpJe/Hv+9dr8hB7KFgpZGP1HxsaAWB8hFUQIuxYylctZriP//0W5HUg6/jWCI5r8N5tTfIBkjGe++R/blRckpYrXR8fkqDKA1R1/BU3bVFRVEUHYMtKgpbVDTis+lBDV9s6Cw7tIyDxw9iFmYkskIHsYeClSsJO3MQ1hDXf1KCoAKajRkDFgvHV648pSDwOI6/3fktyw8tZ3vWds7rqIRBVbHFx5Px3vuciC9tESpdLlx5+X5rPylqjrdchWbCFh1tOJFLSsj55hvlRK4GOcU5vBH9Bme3O5unIp9ia8bWU5asd6RncCIhgfZPPF67A/WDEgQVYGrRgvAzzyRv/k80v/wKIs6r+IEY0mEI57Y/l2JnMR9u/ZAMWwZX9LxCOY4rSdGWLUaxOF03nJdmM9LlClkoXWPDCJsOQ9rtoOsULFlKwdJlKsS0ikxeP5l8ez7PDHyG8zqed9oXwuwvvwTA0r17LYzu1CgfQQXYEhIo3rEDPT+fwxMmnNZ2KoTgln63YNftfJX0FRNXTFSO40rgzMri6LOTS8tECEHLm25STuFaxOtEfuIJWlx3nbHRrR0oJ3LlmPnnTLYc24KUkqlbpp722bclJJD7ldHu/dgL/wm5b0YJggook2XsqDjL2Jf9+fsRGM1SlOP41NgSEkh/8y0O3HwLzrw85RQOMR4ncut77i7jRC7evp2s6dNDPlHVZTJsGXz252eAETRyOgcxQMHva6o8vwQTZRqqAG8v4+JikJKI4afPHvY4jotdxUgkbcPbnvYzjRFbQgKH759gmCKATtNeIaxvX2WXrgN4nMhFW6IoWL2KgpUrKVi5EhEerjQ0P+hS54WNLyCRWDUrLuk6rYMYMOYVMEyhdcAEqgRBBXjU5Zyvv6Zg6TL/LQ/LMaTDEGaMncHGIxv5YfcPfJL4CceKjnFBlwuUv8CHvJ9+9goBNA1XTm6ZOvuK0OL5WwhNI3NHEmBMXLaoaPU38iExI5EZf85gy7Et/Gfkf+jfuj9x6XGV6ml+IjERa69etLzppjrx8qMEwSmIGDqU8P792bNmLfm//Vapyo1DOgxhSIchtAlvwxsxb/DpH58ye/tsVaEUt130u+84vniJ0TdYCITVGvK3IYV/Is4fgQgP92Yi5y9bhpQ6TUeNCvnEFWoSMxJ5aPlD2HU7mtDo37q/99k/HfZDhyjeto0OzzxD2wcfCP5gK4ESBKdBa9qU5pddRsHSZXSaMsWwZVcCm9OGQCCRlLhKGn2PY1tCAsn33Q8OBwhB+2eeBrujTrwNKfzjW6uoeM9uChYvIWv3brKnf97ocw02HtmIXbd716vyfOcvXgxC0OKaq4M0uqqjnMWVoMV11+HKy6No8+ZKfyayYyRWk9FEXCLp0KRDsIZXL8iZ/aUhBMAojWx3KIdwPcDjRA7vP8DQ4gBZUkLRxo0hHlno0KVOzLEYADQ0rJq10hWIpZQc/20xEcOHY+nUKZjDrBJKI6gEzS4cjdayJfm/LabZJZdU6jOeRLN1qev4afdPfJjwIYcLDnNh1wsblWZgS0gge+YsClevNgRACEvtKqpP+VyDvAULQUqaXnRRoxLmiRmJfPrHpyRkJnDvoHtp06RNpXwCHoqTkrAfPEibByYEdZxVRQmCSiCsVlpcdRV5Cxdi7dmDpqNHV+rm99gMuzTtwitRrzD9z+l8ueNLZo6d2SiEgS0hgeR77zMaymsaHZ9/Hr2wUJmD6iFl2mGmHSPvu3lkffIp2bO+aDRmosSMRB5c/iAO3YEmNMb2Glvl5zj7i9mgaVg6dw7OIKuJMg1VkrABA6CkhKyPP+HwAw9WKa46357fKPMLsj+fYQgBACHQCwuVOage4zETWTp1LmsmWr8hxCOrHZYdWoZDN8ybAlHl59gWF0fB4sWg66T+4/E6lZuhBEEl0QsKjAUpq5wA4skv8DiPc07kBGmUdQNbQgIpf/87hWvWGOYgd6KYMgc1DDxmIk8bzLyFC8n8v/+rUxNboMm0ZbLs4DLA8AtUJlegPLk//uhdrgtJZL4o01AliTh/BJjN4HQizOYqTWqe/IKYYzFsPLqRb3d+S15JHrcPuL3BmYjKm4M6vfA8roKGZQ5KzEj0xosDVVpuCH/vMmaijHTy5sxt0Gai6KPRTNk0hQJHAdNGTyPrRFa1/paOI0eNhTr4YlQjQSCEaAN8D/QCDgG3Sylz/Ry3DBgJbJRSXuezvTcwD2gLxAP3Sint5T9fF4gYOpSuH37Ikcceo/m4cVW+2T3+gnPbn8vDKx/m1wO/svzQcmZdNatBTA4ecmbPLmMOchUU1qvqof4m+SEdhlBQUkDUsSjMmpnvdn2HU3eiCeONWJc6GhpSSKSUCARCCO92hHGMWTPz/PnP09TSlP15+xnddbT3HPVNYHiSzrKmf26Yidy1iYo2bmpQgiAuLY5JqyahSx2LZqFXi17c2PfGKn+PI8OoNNri+usI69uvzr0Y1VQjmAysllK+KYSY7F5/1s9xbwMRwCPltr8FvC+lnCeE+Ax4CPi0hmMKGi0uv4z8Sy7Btnkz0q0ZVJVt2du8JiK7bue3/b/V+Ye+MtgSEows7BUr63x0kO9kP6TDEOLS4lh9eDV2l52f9/6MUzq9Ph1JxRnlLukqXcaF51CJ9C77bnfoDl7e8rL3M576NIBXeEhp1LD/x9B/4NAdDO9k/H51VUCUjybK/+03EIKmoy+oUxNddXDoDqZFTUOXRk0gXerVzgfKX7gQXC7a/f3vhPXuHeCR1pyaCoLxwBj38lfAWvwIAinlaiHEGN9tQggBXAbc7fP5l6nDggCg1W23krp2LYXrN9D8skur/HlPfoHdZUciWbB/ASbNxNW9r65zD3llsSUkcPi++5EOhxEdNOU59CJbnXrrScxIJDYtFiR8+uenOHUnQghaWVuRU3Kyz6YiASAQaEJDSolJMwHGBGESxrJLusose7QGl+7CrJkZ0XkEm45sOun7JYY2AWDX7bwb/653n+d8Fs3CvyP/jc1pqzNCwddMZE85TP78n8j66COyZ8yo12ai+PR4Xot6jQP5BzBrZu/vX1W/ABi5A3nz5xMRGVknhQDUXBB0lFIecy+nAR2r8Nm2QJ6U0m1HIBXoWtHBQohJwCSAHj16VGOogaHZxRdjat+OvPnzqyUIPP6CuPQ4pJT8X8L/MXfXXObvmV9vzUR53/9gCAEwooOKbCEzB/m+7RfYC1h6cCm5xblsProZnbJtRz0Tr0dD80zyQIUTu0Wz8MzwZ8i351fLRwCGucGhO/yew6k7vaYkD55lu27njZg3vOMb13scLa0tuarXVWhCC5nWUJGZyBYVVS8FQVxaHA+teMgw5wkzU0ZM8f69q/Pb2mJjcSQfpv3f/x74wQaI0woCIcQqwF8K3PO+K1JKKYQ4fWW2aiKl/Bz4HCAyMjJo5zkdwmKh1Y03kT1rFunvvkfzyy6ttr9g5raZaGjo6Nh1wyxR3wRBUXSMN2U+VJUUPZM/Ej754xNviF9FeCZdi2bhsaGP8d/Y/+LQHZWe5Mv/jXzXK7PseRGo6BwtrS29Y6pIQLiki8UHFgMwd9dcrynLoll4evjTFDoKa10oeM1E7tpERTGxtJ00CWEy1doYaorNYePlzS97f2eJJN+ez8SzJ1b7O7M+/xxhtWLuWJX35NpFyEpU1azww0LsBsZIKY8JIToDa6WUAyo4dgzwlMdZ7DYNZQKdpJROIcQo4GUp5VWnO29kZKSMiwtdLP7xJUs58uSThh28Bl2cEjMSeXjFw14zkUBwVe+ruHvg3fVCIOT+8ANp017F3KEDnV56kZKdu2rNHJSYkcj61PXkFOewYN+CMvb68mhohnnF3UPWd7If0mHIST6DuoA/p7WvgBAIXNJ1Sh+GWZh54KwHaGJuwvBOw2vl2mwJCdhiYnGmp5E79zuaXjqGiHOHEHH+iDqvHWxM3cjULVNJs6WVMQfVpGBkwfr1pE4yXKN1oZS3ECJeSnmSfaumpqFFwP3Am+5/F1b2g24NYg1wK0bkUJU+H0rsKSnGgk9OQXX+uL5mojAtjLfj3mbpwaWsPLSSL8Z9wdAOdffByVu4kLQXXwLAlZWFqXnzoJuD4tLi+O3Ab6QVpbH56Ga/k6DnDfpU5hx/b/R1RQB4KD8mz3K/1v0q1hrcVlbP26xTOpmxbQZg+BkeOushzml/Dvvy9gVN6PmWE3edOMHxXxZQtGZtnZgET8Vv+39jysYp3peF50Y8VyNzkIdcdztKoEZzRbCpqSB4E/hBCPEQkAzcDiCEiAT+KqWc6F7fAAwEmgkhUoGHpJTLMRzL84QQrwIJwKwajqdWiBgxHGGxeJ2jNTGF+JqJPLZqp3Ty8uaXmTxiMjuyd9SpN1UwyuimTX3Fuy6dzqDd4GtS1rBw30KOFBxhV+6uk/YLBCZh8vu2D3U32qa6+AoIj1CoSGsAQyhIJLrUvUIBDEF535n30czSjBGdRwTl9wnr1cu77Gl7WdcmwcSMRObumsuKQyu8Lxa61GtsDgLQ7XZOJCXV+Sg6qKFpKFSE2jQEYNuaQOo/HkNr1pwzli1FuFPuq4vHTOSr9oMx0VlN1jrTz+D4qlUce24K0uUCp9PbZL6mb3seU8iwjsM4mH+QxQcWc/j4YdJsaScde6rJvy78RqHG81v6CgWzZubibhezKnnVSZqUJjRu6nsT7Zq048KuFwKBEaC2hAQOP/CgtxtX83HjCB80qM5Ek8WmxfLwioe9z5pFs3jzBQLxvOUtWMCxyc/RYfJkZElJnbjuikxDShDUgNwffyTtPy/S48svaTry/Bp/n69deMG+Bfy09yfAmPhu7ncz3Zp3C+lkV7BuHal//RtIibBa6fj8FFx5+dW6wX2vtcRZwt9W/+2UTt7T2fkV/inva/C8bECptlAekzAhpcRqstb4d7YlJGCLiub4ihWU7NxZY79aIEjMSGTxwcUs2b+E447jgHHNt/S7hc7NOgfknpJScvDmW5AOO31+/bXGL4qBIlg+gkZNy+uvJ/O998n56quACILyduHfDvxGiasEieTnvT+HVDtwHT/OsRf+423ZKV0uXHn5VfIL+L6pvhX7FnaX3XhAJCeFdoKa/ANB+XvK45OqyIQEpYlyxa5ipkVNA8BqsvLs8Ger/Pt7fAZSSkMQSIm020NmJoo+Fs0jKx/xXqNZmL331/VnXB+w+8oWE0vJzp10emVqnRECp0IJghqghYfT+q47yfrkU9LffIvmV40N2M3t6WcQlx5H9LFooo5FIZEUu4r5YfcPtWr7LtoSxdEpU3BlZxu+EV2vtL3TM/m3sLTgrdi3cOiOMm+hUkr6terHoeOHTkrMUpN/4PHnYyjveJZII1yV0sS6EldJjYRC01EjyZ4+3Qgt1XWcWZlkTf+81swlCRkJfJv0LetT13uFgCY0bu53c8C0AF8y/+//EOHhWHv2DNh3BhNlGqohBWvWkPo3I1EkWJERiRmJTFwx0RtmCqW+g+q8pVWFothYDt8/AdyTf8cXnj+tOcgz+ZuEif8l/M+vyUdzF771aDhQP2vuNBQqClf1aG5wcra1QGDRLDw7/FmOO46f9u9lS0igaMMGcn/+BVdampF3YrUGzUzkySZPK0rjxz0/escfqNDQisj/9TeOPv10nTCDlUeZhoJEyZ693uVgqby+2sGenD0sPbTU2wt5WtS0oJmMpMNB2stTQXcn1+h6heagxIxENh3ZRL49nx92/+A3rt+TtWvV/NueK0rAUgSfyoSrvhX7Fg6X25SE7q2XNS3a0BTMmplHznkEkzD5zVvwmolcOtnTp4OuB+2ZiToaxd9W/82r2XgwCRM39w2OFuAha8bnxkINw8trEyUIakiZbEohghYe5nlQEzMSWZOyxqsdeP4rdhXzTuw7DOs0jEu7G6UvavJWbdu6lbRXpmHfv98ovy0lwmLhSL9WLNg2k8iOkeSV5LHs4DKyTmQRkxZzyrh+Zeapn5zKlORp3u7NW9CdfJz4MWAI/XsG3kP/Nv1JL0rn/M6GDy0uPY7h53YhzPPM6Dq6vaTaZiJfTebMtmfy/a7vWbh/Ifvy9pV5GQmWL6A8xbt3Y9+zF9zZ1HU5ZNQXZRoKALaEBDLefpsTiX9wxpLFWH3ip4OBP6erv5BAT+SHr/kITl8Pp3+q5NBf/oLQdaTJhP2f93Hg8J8c7d+GWXKDN3HJHyq0s3HgL0QVKo5EgrLF855veisttiTRZfUOsJ0wSpRYLTjef57Y9scrda9GmCJ4N/5dI+RaGHWiPBqAhoamlZ6vtu7DlMcewxYdQ9f336N4R1KdCBn1RYWPBhlnVhb7rhxL8yuuoOvb/62183oeyKOFR/lpz09+o288+JY6LlNL32fZhMY788LoctAIq3MJ+P5ijQUXVNzMTkX3NG78CQWLZuHyHpez5OCSU5bBuG+Vi2tjJQLjXvvxEhO/jDLuJ08+je/96Um69PednZp2Ir0oHYkMeDhoZTixfQeHbr2Vdv94jPaPPhr081UH5SMIMuZ27Whzz91kz5yFqXVrWlxd9eY11cHXZPTr/l/LRH64dFeZQmW+pY7L1NL3LEvJHWuddDlYgssd8eY0wY4exorv276K7lF4OFW28+rDq09ZPG/LQI2xW11YXKBJONzWmOT93p+c7LDWhAbSCDqYdPakMoIomCYgf6S9+ioiLIwmdUgDqCxKIwggoS4wdbpCZb41abwTuW68cfU7onPbWgdnJUuyrx7Oxx220T/Zwa6eZvZ1M522BLOa/BX+qMw92TfVxajtDi75U6coHDYM1kjob2ZfdxMu3VVh34e6VEAwZ85c0qdNq5ORQr4o01AtkDX9czI/+MBIutI02j/xRJ1o03i6HrvDM1tg+ccrCJcLqWn0+vYb9nTTVDinImj4uydHROdjfmemcYDVguODF07rL6gL96N0ONg7ZgyubHeDI5OJ9o8/Xiee/fIo01AtYEQQWZHFRj32JkPODfWQgIpDAwHObX8uhyc8gM3lrm0kBLbYOIacN0mFcyqChr97Mmvz52RqmhGubHfQKSGVif9+sswx/pZDTc7X3+DKzqlysmVdQgmCAGK07ZtN3k8/kz9/PrboaJqeX/PSE8HCtjWBzA8+wBYTU+/C3RQNj4gRwxFWq7f/cd5P8wFJ88suq5NmFgBHWhqZH39Ms0svpe2kh42cgToWKVQZlGkoSBx5+hmOL1lCm3vvDWjpiUBh25pA8r33gssFJhOd/vMfXPnVKyCnUAQKT2MbV0EBOTMNM5EIC6uT/Y9tCQkce/El7IcOccbSJVi7dQv1kE5LRaahimMCFTWixTVXg8tFzpdfcviBB7ElJIR6SF6k00naa68ZQsCNK9/IGK5rD5uicRExdCjtHpmEqXlzo44/Ri+D40uWhnhkZbElJHD4vvux790LUuLMzAz1kGqEEgRBomTPXiNJhtKmHHWBothYDtx4EyU7dhgZwyaTMgcp6hweM5FHGOT++CPHXp5aZ16oClatMhpTAUhZZ57v6qJ8BEHCW3qiuBikRGvWLNRDonDDBlIe+avhjDOb6fSfF6rdT0ChCCaGv+0LbDGxSCRZH3xI3rx55P/0Ez2+/iqk96t0OChcv8FY0bQG8SKlBEGQ8NzIRRs3kffzz2T97384MzJoNuaSWr+JbQkJFKxYSd78+d4CckhZ5X4CCkVt4ilSlzX9c0O7dhdxy/jvf2k2ZgwRI0aE5FnKfPc97Hv30u6JxxGaqUG8SClBEEQ8N7KlR3eOPTuZ7OnTyfnqS3rMrj3Hly0hgcP3TzAiMaBMAbn6/hajaBx4o4kcDtB1TiQkciIhsdaTNj1+AelwgMlE05Ej670A8KAEQS3gTEsvfaMpLsEWFVUrN5CUkqxPPysVAppGq1tuwdKlS4N4i1E0DnzNRI6jR8j74Uf3s1RMwe9rau0+zvvhx1K/ANSL8tKVRQmCWqBMqWopOf77GqNk9fnnB+1GKty0ifQ33sC+b7/hcBMCYbHQ8sbxDebmVTQePNq1LSGB/IWLvM9S7pw5uHKyaXXrrUG7r20JCUZu0KJFxgtdA/EL+KLyCGoJT3z0ie3bKFy5CgARHhYUM1H2l1+R8eabxorZTKcXXlA5AooGg+dZkk4HWR99bJR0MZnoMXMGTUeNCvi5vKZVIWj/zNNgd9TbZ0mVmAgxvo6vwlWrvWaigtWrA3ZDFW7cSMZ771OSlFS6UUpvjoBC0RDw50TG5SL1iX/S/rFH0U8UB2SilrpO1ieflJpWhQC7o0E+S0oQ1DJeM5E7jT732zm4srJpdcft1bpxbQkJFG3aRMnuPRSsXGls1DSE2Yx0uRqcCqtQePB1IguTCWE2k/76GwAIq5WOz0+pdnh04YYNpL/xJvYDB8qYVhvqs6RMQyHAq9q6XGT973/eaqVtHngAU4sWp71xPZ/HYibzvffBWa5jmMlEq1tvVU5hRYPH8yxEjBhO0ebNpaYiMN7ghTCEwnOTTysUbAkJFK3fwImkJIrWrTM2NjDTalBMQ0KINsD3QC/gEHC7lDLXz3HLgJHARinldT7bvwQuAfLdmyZIKRNrMqb6gF/VVtfJmTXLOMBspvO0aVh79sAWG+d9CynaEoV02MmeOQt8ohcA43tMJm9oqHIKKxoDnmfJQ/aMmUZkj9RBl97oorSprxjPhtVKxynPeYUCQOG6dThzcsn/6acyZVeARmNarZFGIIT4L5AjpXxTCDEZaC2lfNbPcZcDEcAjfgTBb1LK+VU5b33XCDzYEhI4/MCDpSFp5W/CyuBTNbQybz0KRUPGoyGYWrUk/fU3Su37VZnnyr1U1dUmM9UhWM7i8cAY9/JXwFrgJEEgpVwthBhTfntjxzc+2tSqJelvvGnYO81mmgwdgi0q2v8H1eSvUPjFV0MI69+/rFDwvHDpfvp6l5v8G9tzVVNB0FFKecy9nAZ0rMZ3vCaEeBFYDUyWUpbUcEz1Cn83rkdlPZz4oDeLUYDX+dvYblKFojpUKBTcL1zquSrltKYhIcQqoJOfXc8DX0kpW/kcmyulbF3B94wBnipnGuqMIUCswOfAfinlKxV8fhIwCaBHjx7DkpOTTznuhoCvIwyot00vFIq6RGN+roLSs1gIsRsYI6U85p7U10opB1Rw7BjKCYKq7PelofgIFAqFojYJVmOaRcD97uX7gYVVHFRn978CuBHYXsPxKBQKhaKK1FQQvAlcKYTYC1zhXkcIESmEmOk5SAixAfgRuFwIkSqEuMq9a44QYhuwDWgHvFrD8SgUCoWiitTIWSylzAYu97M9Dpjos35RBZ+/rCbnVygUCkXNUa0qFQqFopGjBIFCoVA0cpQgUCgUikZOvSw6J4TIBKqbSNAOyArgcOoD6pobB+qaGz41vd6eUsr25TfWS0FQE4QQcf7iaBsy6pobB+qaGz7Bul5lGlIoFIpGjhIECoVC0chpjILg81APIASoa24cqGtu+ATlehudj0ChUCgUZWmMGoFCoVAofFCCQKFQKBo5DVYQCCHGCSF2CyH2udtolt8fJoT43r0/WgjRKwTDDCiVuOYnhRBJQog/hRCrhRA9QzHOQHK6a/Y57hYhhBRC1OtQw8pcrxDidvffeYcQYm5tjzHQVOK+7iGEWCOESHDf29eEYpyBRAjxhRAiQwjhtyKzMPg/92/ypxDivBqdUErZ4P4HTMB+oA9G05s/gDPLHfN34DP38p3A96Eedy1c86VAhHv5b43hmt3HNQfWA1FAZKjHHeS/cT8gAaN/OECHUI+7Fq75c+Bv7uUzgUOhHncArvti4DxgewX7rwGWAgIYCUTX5HwNVSMYAeyTUh6QUtqBeRj9lX0Zj9FnGWA+RolsUYtjDDSnvWYp5Roppc29GgV0q+UxBprK/J0BpgFvAcW1ObggUJnrfRj4WEqZCyClzKjlMQaaylyzBFq4l1sCR2txfEFBSrkeyDnFIeOBr6VBFNDK09+lOjRUQdAVSPFZT3Vv83uMlNIJ5ANta2V0waEy1+zLQxhvFPWZ016zW2XuLqVcXJsDCxKV+Rv3B/oLITYJIaKEEONqbXTBoTLX/DLwFyFEKrAE+EftDC2kVPV5PyU1bV6vqIcIIf4CRAKXhHoswUQIoQHvARNCPJTaxIxhHhqDofGtF0KcLaXMC+WggsxdwJdSyneFEKOAb4QQZ0kp9VAPrL7QUDWCI0B3n/Vu7m1+jxFCmDFUyuxaGV1wqMw1I4S4AngeuEFKWVJLYwsWp7vm5sBZwFohxCEMW+qieuwwrszfOBVYJKV0SCkPAnswBEN9pTLX/BDwA4CUcgsQjlGcrSFTqee9sjRUQRAL9BNC9BZCWDGcwYvKHePbb/lW4Hfp9sLUU057zUKIocB0DCFQ323HcJprllLmSynbSSl7SSl7YfhFbpBGB736SGXu6wUY2gBCiHYYpqIDtTjGQFOZaz6Mu1OiEGIQhiDIrNVR1j6LgPvc0UMjgXwp5bHqflmDNA1JKZ1CiMeA5RhRB19IKXcIIV4B4qSUi4BZGCrkPgynzJ2hG3HNqeQ1vw00A350+8UPSylvCNmga0glr7nBUMnrXQ6MFUIkAS7gaWm0lK2XVPKa/w3MEEL8C8NxPKGev9QhhPgOQ6C3c/s+XgIsAFLKzzB8IdcA+wAb8ECNzlfPfy+FQqFQ1JCGahpSKBQKRSVRgkChUCgaOUoQKBQKRSNHCQKFQqFo5ChBoFAoFI0cJQgUCoWikaMEgUKhUDRy/h+p9A9NA1SpLwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "### BEGIN SOLUTION\n",
    "N = 100\n",
    "\n",
    "for lam in [0, 100, 500, 1000]:\n",
    "    ham = hamiltonian_perturbed(lam, N)\n",
    "    vals, psi = np.linalg.eigh(ham)\n",
    "\n",
    "    grid = np.linspace(0, 1, N-1)\n",
    "    plt.plot(grid, psi[:, 0], '.-', label=r'$\\lambda = %s$' % lam)\n",
    "plt.legend(loc=\"best\")\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-3cfcd31b5b76539e",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Bonus: discuss the structure of the spectrum at large $\\gamma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Create Assignment",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}