{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating bond orientational order parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example illustrates the calculation of bond orientational order parameters. Bond order parameters, $q_l$ and their averaged versions, $\\bar{q}_l$ are widely used to  identify atoms belong to different crystal structures. In this example, we will consider MD snapshots for bcc, fcc, hcp and liquid, and calculate the $q_4$ and $q_6$ parameters and their averaged versions which are widely used in literature. More details can be found [here](https://pyscal.readthedocs.io/en/latest/steinhardtparameters.html). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyscal.core as pc\n",
    "import pyscal.crystal_structures as pcs\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we analyse MD configurations, first a set of perfect bcc, fcc and hcp structures and another set with thermal vibrations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perfect structures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create atoms and box for perfect structures, the :mod:`~pyscal.crystal_structures` module is used. The created atoms and boxes are then assigned to `System` objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcc_atoms, bcc_box = pcs.make_crystal('bcc', lattice_constant=3.147, repetitions=[4,4,4])\n",
    "bcc = pc.System()\n",
    "bcc.atoms = bcc_atoms\n",
    "bcc.box = bcc_box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcc_atoms, fcc_box = pcs.make_crystal('fcc', lattice_constant=3.147, repetitions=[4,4,4])\n",
    "fcc = pc.System()\n",
    "fcc.atoms = fcc_atoms\n",
    "fcc.box = fcc_box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "hcp_atoms, hcp_box = pcs.make_crystal('hcp', lattice_constant=3.147, repetitions=[4,4,4])\n",
    "hcp = pc.System()\n",
    "hcp.atoms = hcp_atoms\n",
    "hcp.box = hcp_box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next step is calculation of nearest neighbors. There are two ways to calculate neighbors, by using a cutoff distance or by using the voronoi cells. In this example, we will use the cutoff method and provide a cutoff distance for each structure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Finding the cutoff distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cutoff distance is normally calculated in a such a way that the atoms within the first shell is incorporated in this distance. The :func:`pyscal.core.System.calculate_rdf` function can be used to find this cutoff distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "bccrdf = bcc.calculate_rdf()\n",
    "fccrdf = fcc.calculate_rdf()\n",
    "hcprdf = hcp.calculate_rdf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the calculated rdf is plotted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7f5cfc5f11d0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApEAAAEWCAYAAAA+dsw2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO29e7QcZ3nm+zx93RfJbAlko5tj42NMDCvYoDgkJAzBEAzDYDhnyJiVgCGe4yQL5sAkmRyY5JzJZMUzzCQhMyQHGC4OZgI4HgjBQwzBOBiSjMHIIAy+YfkCliVbwpJsaW+pu6v7PX9Ufd21e1d3V/W1qvr5rbVXd1dX93631F9/bz3vjWYGIYQQQgghklCYtQFCCCGEECJ7yIkUQgghhBCJkRMphBBCCCESIydSCCGEEEIkRk6kEEIIIYRIjJxIIYQQQgiRGDmROYHkwyRfPms7hMg7JC8g+W2SJ0j+X7O2R4g8oD0sm5RmbYAQQmSM3wZwq5ldPGtDhBBilkiJFEKIZPwYgLtmbYQQQswaOZH54idJ3k3yGMk/J7kAACQvJ7mP5FMkHyB5WXB8a3DeweA1fz1b84VINyT/DsDPA/gzkidJPp/kH5P8AcknSf4DycXg3J8l+b9IHif5CMm3zNR4IdLPRSTvDNbSX8bYw24l+R9J3h685nMkt872T5gv5ETmi18C8EoA5wF4NoDfJXkJgI8D+DcAVgC8BMDDwfn/HcASgOcCOBPAn0zZXiEyhZm9DMDfA3i7mW0CcDWAFwL4GQBb4Ye6WyTPBvAFAH8KYBuAiwDsm4nRQmSHXwRwGYBzAfwEgLcM2MMA4M0AfgXADgAegPdN0d65h5qdnQ9IPgzgPWb2weDxq+FvYF8GsGZm/7rr/O0AHgXwdDM7NmVzhcgsJG8F8BcArgWwCuBFZvadrnPeDeASM3v99C0UInsEe9jvmtlfBI//M4AzABARe1hwzq0Avm5m7woeXwj/Ym3RzJpTMn2ukRKZLx4J3f8B/Cuz3QAeiDh3N4CjciCFGJpnAFhA7/UVdVwI0ZvHQvfXAGzC4LXUve+V4a9NMQXkROaL3aH7ZwM4CH+BnRdx7iMAtpJcmYZhQuSQHwE4jd7rK+q4ECIZg9ZS977XgL82xRSQE5kv3kZyV5BY/G8B/CWAjwJ4K8lLSRZI7iT5HDM7BD9n6/0kt5Ask3zJLI0XIkuYWQt+SPu9JHeQLJL8aZJVAJ8A8HKSv0iyRPLpJC+arcVCZJLIPSz0/C+TvJDkEoDfB/BphbKnh5zIfPFJAF8C8GDw8wdmdjuAt8IvmnkSwFfhtygBgDfBv2q7F8BhAO+ctsFCZJzfAvBdAN8EcBTAfwJQMLMfAng1gN8Mju8D8PxZGSlEVhmwhwF+gejH4IfCFwBoAMAUUWGNEEIIITKHK3Izs4/M2pZ5RUqkEEIIIYRIjJxIIYQQQgiRGIWzhRBCCCFEYqRECiGEEEKIxJRmbQAAPOMZz7Bzzjln1maIWXLfff7tBRfM1o4xcccdd/zIzLbN2g6H1tiMydnne9akbX0BWmMCuVrncddYKpzIc845B3v37p21GWKWvPSl/u2tt87SirFB8geztiGM1tiMydnne9akbX0BWmMCuVrncdeYwtlCCCGEECIxciKFEEIIIURi5EQKIYQQQojEyIkUQgghhBCJkRMphBBCCCESIydSiClB8lqSh0l+L3TsL0nuC34eJrkvOH4OyVOh5z4Yes0LSX6X5H6S7yPJWfw9Qggh5ptUtPgRYk74GIA/A/Bxd8DM/oW7T/KPATwZOv8BM7so4n0+AOBqAF8HcBOAywB8YQL2CiGEED2REinElDCzrwE4GvVcoCb+IoBP9XsPktsBnGFmt5k/s/TjAF43bluFz+fvPIhjq/VZmyGEmBD3PXYC33w48mtZxEBOpBDp4OcAPG5m94eOnUvy2yS/SvLngmM7ARwInXMgOLYBkleT3Ety75EjRyZjdY558lQDb//kt/G5fY/O2hQhxIR43y3343c/+73BJ4pI5EQKkQ7eiPUq5CEAZ5vZxQB+A8AnSZ4BICr/0aLe0Mw+ZGZ7zGzPtm2pmhCXCWpeM7htzdgSIcSkWKt7ONVoztqMzKKcSCFmDMkSgP8dwAvdMTOrAagF9+8g+QCAZ8NXHneFXr4LwMHpWTs/NJoW3MqJFCKv1LxW+4JRJEdKpBCz5+UA7jWzdpia5DaSxeD+swCcD+BBMzsE4ATJFwV5lG8G8LlZGJ13vMB5dM6kECJ/+E6kLhSHRU6kEFOC5KcA3AbgApIHSF4VPHUFNhbUvATAnSS/A+DTAH7NzFz2968D+AiA/QAegCqzJ4JzHr2WNhgh8krda6EuJ3JoFM4WYkqY2Rt7HH9LxLHPAPhMj/P3AnjeWI0TG3DOoyclUojcUvOaUiJHQEqkEEJE4LVzIuVECpFX6l4LzZa101dEMuRECiFEBK6gRuFsIfKLUyHrciKHQk6kEEJE4LWkRAqRd5wTWWvIiRwGOZFCCBFBW4mUQiFEbqlLiRwJOZFCCBGB167OlhIpRF5pDxWQEjkUciKFECIClwupZuNC5JNWy9rpKvWmGo4Pg5xIIYSIoO4FSqRyIoXIJeEQ9mkpkUMRy4kkuULy0yTvJXkPyZ8muZXkzSTvD263hM5/N8n9JO8j+crJmS+EEJNBSqQQ+SYcwlavyOGIq0T+VwBfNLPnAHg+gHsAvAvALWZ2PoBbgscgeSH8CRzPBXAZgPe78W1CCJEV2n0ilRMpRC6phULYmlozHAOdSJJnwB/B9lEAMLO6mR0HcDmA64LTrgPwuuD+5QCuN7OamT0EfzTbJeM2XAghJomqs7MFyWtJHib5vdCx3yP5KMl9wc+rQ89FRsxIvpDkd4Pn3hfMqBc5ZL0SqZzIYYijRD4LwBEAf07y2yQ/QnIZwFlmdggAgtszg/N3Angk9PoDwbF1kLya5F6Se48cOTLSHyGEEOPGVWUrJzIzfAx+9KubPzGzi4Kfm4CBEbMPALgawPnBT9R7ihwQzomUEjkccZzIEoAXAPiAmV0MYBVB6LoHUVdtG76FzexDZrbHzPZs27YtlrFCCDEtnALZ0MSaTGBmXwNwNObpkREzktsBnGFmt5mZAfg4OlE2kTOUEzk6cZzIAwAOmNk3gsefhu9UPh4sOAS3h0Pn7w69fheAg+MxVwghpoNr/SElMvO8neSdQbjbFYD2ipjtDO53H49EEbVsEw5hy4kcjoFOpJk9BuARkhcEhy4FcDeAGwFcGRy7EsDngvs3AriCZJXkufDDAbeP1WohhJgwqs7OBR8AcB6AiwAcAvDHwfFeEbNYkbT2E4qoZZpwCFvh7OEoxTzvXwH4BMkKgAcBvBW+A3oDyasA/BDAGwDAzO4ieQN8R9MD8DYzU8aqECJTNDSxJvOY2ePuPskPA/h88LBXxOxAcL/7uMghYfVRhTXDEcuJNLN9APZEPHVpj/OvAXDNCHYJIcRMaY89lBKZWUhudwWgAF4PwFVu3wjgkyTfC2AHgoiZmTVJniD5IgDfAPBmAH86bbvFdJASOTpxlUghhJgrOuFsKZFZgOSnALwUwDNIHgDw7wC8lORF8EPSDwP4VWBgxOzX4Vd6LwL4QvAjcsh6JVJO5DDIiRRCiAhc+w9P1dmZwMzeGHH4o33Oj4yYmdleAM8bo2kipawvrFE4exg0O1sIISJoT6yREilELlE4e3TkRAoxJTRRI1u0+0QqJ1KIXKJw9ujIiRRienwMmqiRGRqaWCNErnHq4+ZqSUrkkMiJFGJKaKJGtvCUEylErnF5kJsXSlIih0ROpBCzZyITNTRNYzTCOZG+vy6EyBM1r4UCgcVKUYU1QyInUojZMrGJGpqmMRqNUJPxphqOC5E76l4L1VIRC+WiwtlDIidSiBliZo+bWdPMWgA+DOCS4ClN1Jgx4SbjmlojRP6oeS1USgVUSgWFs4dETqQQMyTIcXR0T9TYMIM+mL5xguSLgqrsN6Mzt16MkXBrH1VoC5E/al4L1VIBVTmRQ6Nm40JMCU3UyBbhghpVaAuRP2peM1Aii3jyVGPW5mQSOZFCTAlN1MgWYfWxoQptIXLHOiWyocKaYVA4WwghIgiHs6VECpE/XGFNtVRojzkVyZATKYQQEYQLa5QTKUT+WFdY09AaHwY5kUIIEUG4Ilvzs4XIH3WvGYSzi1Iih0ROpBBCRLAunK2cSCFyh1MilRM5PHIihRAiAq/ZQrnI4L6USCHyRq3RyYlUi5/hkBMphBAReC3DQrkIQDmRQuSRerOFarnQLqzReNPkyIkUQogIGs0WFgMnUhNrhMgfNa+JatEvrDFT7vMwyIkUQogIvKZhsSIlUoi8UvecEumvcxXXJEdOpBBCROC1QkqkFAohckfNa6ESKJEAVFwzBLGcSJIPk/wuyX0k9wbHtpK8meT9we2W0PnvJrmf5H0kXzkp44UQYlLUvVY7J1LV2ULkj1qjhWrZL6wBoOKaIUiiRP68mV1kZnuCx+8CcIuZnQ/gluAxSF4I4AoAzwVwGYD3kyyO0WYhhJg4Xsuw1A5nS4kUIm/Um8HYw7LvCtXlRCZmlHD25QCuC+5fB+B1oePXm1nNzB4CsB/AJSP8HiGEmDpeM+xEanNJOySvJXmY5PdCx/6Q5L0k7yT5WZIrwfFzSJ4Komv7SH4w9JoXBpG3/STfR5Kz+HvEZPGaLTRb5oezi/46lxKZnLhOpAH4Esk7SF4dHDvLzA4BQHB7ZnB8J4BHQq89EBxbB8mrSe4luffIkSPDWS+EEBOi0QqFs6VEZoGPwY9+hbkZwPPM7CcAfB/Au0PPPRBE1y4ys18LHf8AgKsBnB/8dL+nyAGuiMa1+AGkRA5DXCfyxWb2AgCvAvA2ki/pc27UVduGb2Az+5CZ7TGzPdu2bYtphhBCTJ5my2CGdmGNlMj0Y2ZfA3C069iXzMwLHn4dwK5+70FyO4AzzOw285sGfhydKJvIEW5W9rrCGk+FNUmJ5USa2cHg9jCAz8IPTz8eLDi38A4Hpx8AsDv08l0ADo7LYCGEmDTOaXQtftQnMhf8CoAvhB6fS/LbJL9K8ueCYzvh72GOyEiaQxG17OJC1yqsGY2BTiTJZZKb3X0AvwDgewBuBHBlcNqVAD4X3L8RwBUkqyTPhR8OuH3chgshxKRwTmOnxY82lyxD8ncAeAA+ERw6BOBsM7sYwG8A+CTJMxAzktZ+QhG1zOJC15ViAdVgnSucnZxSjHPOAvDZILe4BOCTZvZFkt8EcAPJqwD8EMAbAMDM7iJ5A4C74S/at5mZNGIhRGZwTmNn7KGUyKxC8koArwFwaRCihpnVANSC+3eQfADAs+Erj+GQtyJpOcWFrqvlAipFhbOHZaATaWYPAnh+xPEnAFza4zXXALhmZOuEEGIGOKexE86WQpFFSF4G4P8G8E/MbC10fBuAo2bWJPks+BGzB83sKMkTJF8E4BsA3gzgT2dhu5gs7XB2qdhu8aNwdnLiKJFCCDFXOKdxUUpkZiD5KQAvBfAMkgcA/Dv41dhVADcH0bSvB5XYLwHw+yQ9AE0Av2Zmrijn1+FXei/Cz6EM51GKnOAcxkoprETKiUyKnEghpgTJa+GH1Q6b2fOCY38I4J8BqAN4AMBbzew4yXMA3APgvuDlbvMDyReis8ndBOAdLkwnxkPD686J1D9v2jGzN0Yc/miPcz8D4DM9ntsL4HljNE2kkHY4O9RsXE5kcjQ7W4jp8TGoj10maLQ6PeQKVDhbiLxRDymR1ZIKa4ZFTqQQU0J97LKDUx5LhQJKxYLC2ULkjE5OZCHU4keFNUmREylEehh7HzsxHK5PZKlIlAtUs3EhckY9VFjjciKlRCZHOZFCpIA+feyeCHIg/5rkc5Ggj10wovRqADj77LPHb3SOcX0iy0WiVCyoT6QQOSOsRBYKRLlI5UQOgZRIIWZMqI/dL4X72AVttGBmd8AvuknUx06NkIfHOY2lQgHlItHQxBohckW4sMa/LbZHIYr4yIkUYoaE+ti9truPHclicD/cx+4QgBMkX0S/Z8mb0ZkWJcaEy4EsFYlSQUqkEHkjXFgD+M5kvamcyKTIiRRDceDYGp79O1/AvY89NWtTMkPQx+42ABeQPBBMe/ozAJvh97HbR/KDwekvAXAnye8A+DQ29rH7CID98BVK9bEbM64au1wsoFSkWvwIkTPCzcYB35mUEpkc5USKoXj02CnUmy388Ik1POeZZ8zanMT8/v+8G5c975m45NytU/ud6mOXHTrV2US5WFA4W4icEa1EyolMipTIOeY/3HQP/nH/j4Z6rbuKy2Iispnh2n98CLfc+/isTREpxVVjl4sFlApUOFuInFHzmigViGLBr1WUEjkcUiLnmD//x4dwfK2OF/9vz0j8WncVl8WWCO5qM4u2i+ngciL9cLb6RAqRN2qNVruoBggKa9QnMjFSIueUVsvQaBoeOXpqqNe3HbEMKjRZVlHFdHA5kaWi3/pDE2uEyBf1ZqsdygYUzh4WOZFzilssjxxbG3BmNO6KrdbI3pVbllVUMR3aSmShgHKxkPvCGjPD+265H/sPn5i1KUJMBV+JLLYfK5w9HHIi5xSnwh168vRQ+V5tRyyDV25SIsUgvNDEmtIcTKw51WjivTd/Hzd997FZmyLEVKg3W6iWpUSOipzIOcUpic2W4dCTpxO/PstqXsf27KmoYjq4amw/nF3IvRPpFJgsrmchhqHmNdvjDgE1Gx8WOZFzSnizOHAseV5kltW8LDvAYjo4JbJcCPpE5rzFT5ZznIUYhlpjvRJZKRVUWDMEciLnlLDzN0xeZC3Djlg7nzODtovp4HVNrMl7dbZTYLKY4yzEMNSbrS4lspDJ/WzWyImcU9YpkUeTO5F1KZEixzRCE2vKxfz3idSFlZg3Igtr9PlPjJzIOaU2x+HsLNsupkN4Yk2pWMh9ODvLkQUhhqG2obCmqM//EMiJnFPcYilwuHB2ltW8LNsupoNTHosFojwH1dm6sBLzRq3RVVhTlhI5DLGdSJJFkt8m+fng8VaSN5O8P7jdEjr33ST3k7yP5CsnYbgYDRe+2r11aSglst5sBrfZW3Rt1SWDtovp0GgZKsUCSPqFNXnPiWyHs5UTKeaDutdCtRwKZxf9Fj+tnEcdxk0SJfIdAO4JPX4XgFvM7HwAtwSPQfJCAFcAeC6AywC8n2QRIlU4Fe68bZvw2FOnE28eWU7Ez3KjdDEdGl4LpaI/U9cPZ+f7giPLOc5CDEPNa21QIgGJC0mJ5USS3AXgnwL4SOjw5QCuC+5fB+B1oePXm1nNzB4CsB/AJeMxV4yLWtuJXIYZcPB4sl6RWW4JkuVG6WI6eC1DqeA7kX44O9/qRB7C2SSvJXmY5PdCxxJHzEi+kOR3g+feR5LT/lvE5Kl5XS1+Aocyy2tgFsRVIv8LgN8GEP7XPcvMDgFAcHtmcHwngEdC5x0IjokUEVYiAeBAwrzILOcVtjdMNZYVPWg0WygHm4o/9jDfn5U8OJEAPgY/+hVmmIjZBwBcDeD84Kf7PUUOqHtNVMOzs4PQdhb3tFky0Ikk+RoAh83sjpjvGXXVtuEynuTVJPeS3HvkyJGYby3GhQvpnnem70Q+cjRZXmSWN5126C7njoEYHq9p68LZjZznSWX5otBhZl8DcLTrcKKIGcntAM4ws9vMzAB8PPQakSNqXguV0vo+kf5xpTklIY4S+WIAryX5MIDrAbyM5F8AeDxYcAhuDwfnHwCwO/T6XQAOdr+pmX3IzPaY2Z5t27aN8CeIYXCbxdlbl1AqcK6UyHYo3mvB3yeEWE+j1UKp4JTIeajOzm1hTdKI2c7gfvfxSCSGZBMz88PZoT6RHScy32t93Ax0Is3s3Wa2y8zOgS///52Z/TKAGwFcGZx2JYDPBfdvBHAFySrJc+GHA24fu+ViJNxCWawUsWNlEY8krNDOshMZDmMrL1JE4TUNZadEFgow8+fM55VOodzcrIdeEbNYkbT2ExJDMonLca5GKJFZ3NNmSWmE174HwA0krwLwQwBvAAAzu4vkDQDuBuABeJuZ5e7yNus4J7JSLGD31kU8knBqTZaVC9eeCAjaPJTUPECsx2u1UApyIl1Yu9FsoVjI52cly4VyA3ic5HYzOxQzYnYguN99XOQIt2+FnciKlMihSNRs3MxuNbPXBPefMLNLzez84PZo6LxrzOw8M7vAzL4wbqPF6LirrWqpgN1bkveKzPKEi7Daoi8MEUWjGarODpzIPE+tyXLLrgEkipgFIe8TJF8UVGW/OfQakRPC+5/DiQlZ3NNmiSbWzCmuRxZJ7NqyiB+drOFUPf4GkmXlImzzNL8w1IIkO3ih6myXG5nnCu0sDw9wkPwUgNsAXEDyQBAlew+AV5C8H8Argscws7sAuIjZF7E+Yvbr8NvZ7QfwAAAJITmjHYlTYc3IyImcU/wwrv/fv3vrEgDg0ePxQ9pZzqGaoRL5MagFSSbwWp2cyHI7nD0HSmSGi83M7I1mtt3MykEe/0eHiZiZ2V4ze17w3Nstq/8goie1thIZmljjnMgM7mmzRE7knFLzmu1Fs2vLIoBkbX6cYpHFNjmzUiLVgiQ71L1wTmSgROZ4ao3bVM3y7SwLAXS+9ytR4ewM7mmzRE7knLJOidziK5GPJGjzE67OztqFethxTEH+y8RakKj9yPCElUiXG5nn+dnr1oQ2UZFz+hfWKJydBDmRc0q40eq2zVVUS4VExTXhhZa1TSdse4q/MEZuQaL2I8PjNcN9Iv3bPPeKXLcm8ldcI8Q6osLZVYWzh0JO5JwSbm1DEju3JGvzU/daWChns69WLV1KZNKm/WpBMgUazXBOpAtn51eJDK8JdSwQeSc6nB3sZzm+WJwEciLnlHBOJIBEbX68ZgstAzYvlAGkwhFLRC0Uyk/BhqkWJCnEC02scX0is/Y5T0LKUjyEmCh9w9lSIhMhJ3JOqTfXzw3dtWUxdk6kc7w2L5TWPc4Kda/VdoCnabtakGSH8OzsuegTKSVSzBEqrBkfo0ysERmm1mituwrbvXUJx9caOHG60XaweuEW4OZqad3jrFD3WjhjoYQfnaxN9QvDzN7Y46lLe5x/DYBrIo7vBfC8MZomumi05qtPZEbyhIUYC52cyM4e6C4WlROcDCmRc0q3EukqtOOEtJ3j1Q5nZ2xzrXnNjoqqLwwRgReaWFOagz6Rda/V/nuzdlEoRFJcyLpa7hTWkES1VJASnxA5kXNKtxLZ6RU5OKTtFmDHEcvWoqs3W9gU2J41B1hMh0bT2v0hy3PSJzKr6SlCJMX1N64U17tAciKTIydyTvGVyM5VmJta80gsJdJX7za5cHYzW2perdHC5mq5fV+IbrxWa+76RGa1UE6IpLgIVLW83gWqlIpyIhMiJ3JOqTWa65TIlUV/A3nqVGPwa7314eysLbp6s6O6SIkUUTS8Tk7kfPSJbOGMRadEZuuiUIik1PsokbqISoacyDmlOyeyUPDzQU7HyBHMQ3V2O5ydMdvFdGi0OtXZpbmozm521HmtCZFz2jmRpahwti6ikiAnck7pzokEgMVKMZYTWe9yIrPmiNW8FpYqRRQL1BeGiMRrtlAurK/OzrsSmdWLQiGS4ia2+a12O1SUE5kYOZFzSq1LiQSAxXIRp3LuRHrNFpotQ7VUVOhCRNJqGVrWUSBdyCvPOZG1GfVOFWIW1L0WqsWN7k+1XNSekBA5kXOIma0be+hYKBdxKkahSZZzItu5MKWCrjpFJI2gCrvdJ7Idzs7nZ8V9H6jtlZgXal5zQ1ENAFSLCmcnRU7kHOIcqe5w9kI5WTh7UwabjbcnFRQLqBSlRIqNOMWxu09kPadKpPs+OGMxm31fhUhK3WttKKoB/Gpt7QnJkBM5h7hFsiEnshyvsMa19OmEs7Nz5daeVFAuoFqWEik20nYiXXV2zifWtCML1Wz2fRUiKTWvta7RuKNS1J6QFDmRc4hbJBtyIitFnKrHqM5uZDicLSVSDKATzu6qzs6pEunW80K5oE1UzAU1r7lBRAEgYWEI5ETOIb2UyIVSEadjqIqdsYfZC2d3lMgiqmosKyLohLO7+kTmNCeyk96iYjMxH9S9jYWlgL8G9PlPhpzIOaSXErkQU4l0i2yxUkSB2cqhcknTlaIrrMlOKF5MB9fKpzQnE2tcIU2n2ExrQuSbmrexxR3gwtn6/CdhoBNJcoHk7SS/Q/Iukv8+OL6V5M0k7w9ut4Re826S+0neR/KVk/wDRHI6Id31OSGL5SJOJ6jOrmawwrneZbuuOkU3ja5pFsW2E5nPz0p4PWt2sJgHeiqRKqxJTBwlsgbgZWb2fAAXAbiM5IsAvAvALWZ2PoBbgscgeSGAKwA8F8BlAN5PcmMGq5gZ7kprY3V2IVafyFqG8wq1YYpBuMk0TokkiXKRaOR0Yk37wqpcyGWfPJIXkNwX+nmK5DtJ/h7JR0PHXx16jYSQHFOLaHEHqLBmGAY6keZzMnhYDn4MwOUArguOXwfgdcH9ywFcb2Y1M3sIwH4Al4zVajES9R7h7MUELX5ct/9qOVt5heG/XflfIop2OLvQWR+lQiH3SmSlWMxlOM/M7jOzi8zsIgAvBLAG4LPB03/injOzmwAJIfOACmvGR6ycSJJFkvsAHAZws5l9A8BZZnYIAILbM4PTdwJ4JPTyA8Gx7ve8muReknuPHDkyyt8gEhJW48K4iTVm/RWXmtdsd/vP2qbTCWf7hTVZyucU08HlPrrqbMBXJRs5zYlcr0Tm/sLqUgAPmNkP+pwjISTn9CusabYstxeMkyCWE2lmzeAqbheAS0g+r8/pjDi24dvXzD5kZnvMbM+2bdviWSvGQi8lslouwmxwy57wAsyamhcuKlIRgYjCTaYphZoRV4qF3E6sCRebzUGKxxUAPhV6/HaSd5K8NpTXH0sIASSGZJWehTXBMYkL8UlUnW1mxwHcCl/if5zkdgAIbg8Hpx0AsDv0sl0ADo5sqRgbnZzIjYU1wOBmw/XQAsxacYprlF5VOFv0wCmO5cJ6JTK31dkhJTJrhXJJIFkB8FoA/yM49AEA58HP9T8E4I/dqREvj/zPlxiSTXorkTwpGGsAACAASURBVIX28yIecaqzt5FcCe4vAng5gHsB3AjgyuC0KwF8Lrh/I4ArSFZJngvgfAC3j9twMTz9mo0DGFhcU+tSIrO06TgHOU2zs5X4ny66J9YAfk5kXtWJDSkeKVgTE+JVAL5lZo8DgJk9HkTZWgA+jE7IWkJIzulZWBPsa2nYF7JCKcY52wFcFyQWFwDcYGafJ3kbgBtIXgXghwDeAABmdhfJGwDcDcAD8DYzU8wwRfTKiVwIBtIPciLDV3HZUyJDTmRKKsvN7D74agiCdfYo/MT/t8JP/P+j8Pldif87AHyZ5LO1zsZDox3O7ghS5VwrkZ0+kdV8p3i8EaFQNsntLq8fwOsBfC+4fyOAT5J8L/z1JSEkZ/QsrCnFi8aJDgOdSDO7E8DFEcefgJ+kHPWaawBcM7J1YiL0np3tL6BBFdr1ZucqrloqxmoLlBbCf3tKiwjaif9kVFQNQCjxH8BDJF3i/21TsjHXtAtrwtXZuc6JzG7f17iQXALwCgC/Gjr8n0leBD9U/bB7TkJIvmm1DI2m9Q9nN/XfHZc4SqTIGbVQ+CrMQjluOLu5Tok8fqo+ASsnw7rCmmIRXsvQbFm7oXQKiEr8fzOAvQB+08yOwU/y/3ronJ4dEABcDQBnn332xAzOG15zoxJZKsxBdXaO84TNbA3A07uOvanP+RJCckp4zGc3bl+LM3RD+Gjs4RzSqzrbOZGnB4w+rHut9jSPtISE4xJulF4tpyuJetyJ/0r6Hw63yZRDOZHl4hz0icyxEimEo1dNABBWIrUG4iIncg7p12wcAE4PyImqe622A5bSkHBPnIpKsu0IpygHTIn/KaBXn0gvpxNr2rOziwVUS8X2YyHySK+JbUCosEZKZGzkRM4hNa+JUoEbQrjt6ux6/wVU61Iis6Rc1L1Wp1F6+to5bEj8Dz3XnfivDggTIqpPZLlQaE+yyRu1pt+yi6Qfzk7h32lmA4cgCBEH5yD2LaxJj7CQepQTOYeE+zyGWSjFy4nMdHV2WEVNUTsHJf6nh159IrP0OU9CrbF+PTeaqcsTxh/8zT24//BJfPxXNDhGjEa4Q0c36hOZHDmRc0i4z2OYhYpLKh7cJzJcnZ0GJywu61TUFDmRSvxPD53CmvXV2asDcoWzSvd6BvxN1EUm0sADR07igcMnZ22GyAEdJXLj5ztNwkJWUDh7Dqn3aLQat8VPLfNK5MYNUwiHy30MV2dXisxtYU04MpFWJWat1sRa3Zu1GSIHdKqze4ez0/b5TzNyIueQcIueMO0WPwOrs5vrNp16s5WZfKWa12wrkZ2rznwqTGI4GlF9IguFXDcbD48xdcfSxMmah9VaumwS2cQVjvUtrJETGRs5kXNIvRkdzi4XCygVGGvs4cZNJxuLrjuf0x0TwhHZJ7LI3BbWhNdEWsN5q3UP9WZLa1WMTHhWfDcSFpIjJ3IOqTWiC2sAP6Tdr9Gqma1zQrPWV6sWEbpL24YpZkvDhbML4bGHhfY4xLyxbk2Ui+1jacKpkAppi1Fpt7grRuREpqx3cBaQEzmH9FIiAWCh0n+ModcymCH1OVS9kBIpBuE1WygViPDYyVIh37OzXS5YCnunAgBWa77zeLImJ1KMRj8lsvP5154QFzmRc0g/JXKhXOhbWNPd7T9z4exmWIksto8J4fBatm5aDeBXZ+d57GF4eIA7lhaaLWtf2K7ltEJeTA83F7tS3LgHlooFFAv5bec1CeREziG1ZguViOpswIWze39Rd0IB2VTzunviAelTXcRsqXutdfmQgD+9xstxOHtjsVl6/tZwCFtKpBiVdoufCCUS8NeA9oT4yImcQ2qNZt+cyH7h7Ho7FJDNNjm+Eulsz5YDLKaD12ptVCJzXZ2dzgb8jrD6uKYKbTEi7XB2DyFF8+OTISdyDumXE1ktF/u2+HFXaOGxh+HjaafWaGY2FC+mg9e0dUU1gK9E5ro6u9iV4pGiNRFWH6VEilGpd6VkdVPNWO/jWSMncg4ZWJ3dZwF1L8CshbPDDnTWbBfTodGMyolkuwl53lhXWJPCFI+w+rgqJ1KMiPts99oDpUQmQ07kHBIuLulmsVzE6b5K5Ppu/1kLCavFjxiE14rKiSyg2bLMNNVPQmQ4u0+br2kTVh/V4keMSt1rgcSGaIOjWipmZj9LA3Ii5xA/JzI6H2ShXOibE5n16ux1IxvVzkFEEB3O9j8reazQjgxnpyh0v7ounJ0ehVRkEyckhFt4hVFhTTLkRM4h/XIiFysxq7Mz6ESa2bq54SQzN/tbTJ5GM6qwhu3n8kZYiWyv5wFTq6bJal1KpBgftdAeEIXC2cmQEzlnmNm6kG43C4Oqs5vrK9vSqFz0omN752+vFnXVKdbjtWxDOLsUOJV5q9D2mi00W7ahY0GaNtHwzGwV1ohRCUejoqjKiUyEnMg5w02ciWq0CvhOZN9m413D66spVC560d3jEoCUSLGBRrOFUmH9+igHTmXeRh+6C6sNY0xTtCac+rhcKarFjxgZv5Cst+tTKRXlRCZgoBNJcjfJr5C8h+RdJN8RHN9K8maS9we3W0KveTfJ/STvI/nKSf4BIhmdPo+9C2saTYPXQ1ns3nTaFc5ZUCIj/na1cxDdeE3bcJHlnMq8KZHtxsvBOi4VCygwXUqkUx+3ba7ipMLZYkTqMZRI7QnxiaNEegB+08x+HMCLALyN5IUA3gXgFjM7H8AtwWMEz10B4LkALgPwfpK9ExDEVKlFqHFhFoMm4r3a/LQdsRQrF72I+tuV/yK6aTQ3Vme7x3nLiexOT3H303RRuFZvYqFcwBmLZawNGc4m+TDJ75LcR3JvcExCyBwyKCdShTXJGOhEmtkhM/tWcP8EgHsA7ARwOYDrgtOuA/C64P7lAK43s5qZPQRgP4BLxm24GI5OYUzv6mwAPRuOZ7k6O1qJVDsHsZ5Gy9o5kA4Xzs5br0inRIaVmWq5kKr0lJM1D5uqJSxXSuvyI4fg583sIjPbEzyWEDKH9KsJAAJhIUUtrtJOopxIkucAuBjANwCcZWaHAN/RBHBmcNpOAI+EXnYgONb9XleT3Ety75EjR5JbLoZiUKPVBadE9thENszOLmZRiezsBxVddYouvGYL5a4WP51wdvo/50mI+j6oFNOlzq/WPCxVSliuFsddWCMhZA6pe80B4ex0KfFpJ7YTSXITgM8AeKeZPdXv1IhjGy7fzexDZrbHzPZs27YtrhliRAaNfFqsxHMi3ezsUrGAYoGZcCK7Q/Hufhq+MBRuSw9ec2N1druwJm85kRHfB9VyunLCVmtNLFdLWK6WRmnxYwC+RPIOklcHx0YSQgCJIVlkkBJZLaVLiU87sZxIkmX4DuQnzOyvgsOPk9wePL8dwOHg+AEAu0Mv3wXg4HjMFaPSPXGmm4UgzN2rzU/37Gx3PwtqXr0Z2F5anxOZog1T4bYU0Gi1NoSz20pkzqqzo74PqimrTl2teViuFLFUKY3SbPzFZvYCAK+Cn9f/kj7nxhJCAIkhWaQew4lMg7CQFeJUZxPARwHcY2bvDT11I4Arg/tXAvhc6PgVJKskzwVwPoDbx2eyGIUo5SGMUyJ75UQ6h6scUmpS5oj1JCr/K+WFNQq3zQCvaRvC2eVSPifWdMLZoRSPlF0UrtU9LFdL2FQtDq1EmtnB4PYwgM/CXy8SQuaQeIU1rVyOOJ0EcZTIFwN4E4CXBaG2fSRfDeA9AF5B8n4Arwgew8zuAnADgLsBfBHA28wsPd9Ic07UphFmYUB1dq25cWRUVq7calHNxtPjAI893KZQ23B4zY1KpHMq85YTGZXeUi2n68LqZM3DctVXItfqTbQSFjeRXCa52d0H8AsAvgcJIXNJbUBOZKVUgFn+LhgnRWnQCWb2D4iW9wHg0h6vuQbANSPYJSbEwJzIcn8lstbY2GMrK9Vs0UpkakJ3LzazgyTPBHAzyXv7nBs77xjAhwBgz549+kaMSaNl65R2oDOxJm8bS3Q4O11O5Fq9ieVKCZuq/na1ljxf7SwAnw0ufEsAPmlmXyT5TQA3kLwKwA8BvAHwhRCSTgjxICEkVwwOZ3emsPVzNoXPQCdS5IuBOZFB+5uehTXNjQuwUiq0Vb40E90TLx1KZDjcRnJduM3MDincNj28iIk1pZxOrHHfBwvl9RdWT55qzMqkDfhKZAlLVX/drias0DazBwE8P+L4E5AQMnfEafED+FPY3IWL6I3c7DkjqkI5TJzq7O5QeFZ6LXaPbATS0eJH4bZ04TUN5Q3h7HxOrOm07ErfhRUAmJmvRFaL7Q1d87PFKMSZWANkYwpbGpCbPWd0lMjonMh2OLtndXaPcHZKNp1+dI9sBFLTE0/hthRRb7Yiwtn5zIls50iX03Vh5ah5LTRb5rf4qQThbM3PFiMwsLCm7JTIfK31SSEncs4YlBO5MMCJrHvNDSMTq8UC6inZdPoR2ScyBT3xFG5LF16rT5/IKU6s+eQ3fohb7nkcH33LT07sd/TsnTrimjjdaOJv73oMr33+jnVFeElxquNypRPOlhIphsVr+hclfQtrAlU+BeJCJlA4e84YNLHGHT/d4yqs7rXWqRZAOhyxOEQ2Vi6qnYPoYGZotmxjTuQMJtbc/tAT+Nr9Ryb62YxcE2MoNvvC9w7hHdfvw/7DJ0d6H6c6+i1+AiVy+IbjYs6pR3To6KYdzs7AnvatHx7D7Q8dnakNciLnjEFKJEkslos9cyJrXmuDEpmSkPBAukc2Ap3JO3mruhXD4T4HvcPZ0/ucHF1roNG0iSpv7Y4FxfVK5KgTO544WQcA/Ci4HZaOEum3+AkfEyIp7vMeq7AmA9G1//A39+Df/8+7ZmqDnMg5Y1B1NuBXavZrNh6VE5mFq7aa10SxwHU9AN3mmYUvDDF53ESaDX0iXYufKVZnH1+rB7eTq5SuNzeuiXH0fT266tt+bG00J9KpjuuVSK1VMRydvPj+zcaBbCiRjx4/hUePn5qpDXIi54y610KBGzfJMIvlYu+cyIgWP1lpNl6PUlEz9IUhJo9TIktdE2tKhRkokYEj5m4nQa0RvZ5HTfFwzuOotreVyGpx6BY/QjjiKJEuOpX26Fqj2cLjT53G8bXGTNeEnMg5Y1C3fgBYqPQJZ2e42XhkPqfaOYgQLudxQ4uf9tjD6X1OjjknckQ1rx9RPfPGMbGjrUSO6EQ61TFcnb2q6mwxJC7i1L+wJhvRqceePA1X53foydmpkXIi54yoPo/dLJR6O5G+EhnRJzIDTlhkPme7sWz67ReTxwu+lTdUZ7vCmilVZ9e8JlYDB2pUR6wfUekp4Ykdw3Js1Q/Bj+oAh6uziwVioVzAqgprxJDESedqt/hJuRJ5MBTGPnBMTqSYEnFGOS1W+oSzM5wTGa1Ejr5hivzglMZyr4k1U/qch/Mgj00wJ7LmNTdeFLb75A2vxDjncWQlstbJiQSATdWSCmvE0LSdyHJvIaWjRKZ7TzgYUh8PHj89MzvkRM4ZUTlQ3fjV2dELKCocnqbmxP2QEikG4XIey6XonMhp9YkM5xJOUomMDGePYRPthOJHc4BX2+Fsf9NfqpTajqUQSYnq0NGNu4hKuzDiHMcC16uS00bNxueMWgwlcqFc6JkQH7XpVEsFNJqGVstQKAzfWHjSRE0q6OREpt8JFpPHKZHdfSJJolTg1PpEhh3HSeZERoazR9xEWy1rF9aM6gCfrHkoFdje9JerJZxUTqQYkqgJTd24PSLtSuSBY6fw9OUKFspFOZFievhK5ICcyD59InuFswE/JLxQ6P/esyQqlN/pCZbuLwwxHXr1iQT8kPa0ciJdCLvATqufSRB9UTjaJvrU6UY74X/U6uy1moflaqk99WZTtahm42JoanGUyIz0iTx4/BR2rCxioVzAgRk6kQpnzxmxciJ7OJFm5m86Ec3GgfQ7YrVGdCgeSL/tYjq0+0QWNq6RcqEwtepspz6evXVpsi1+InIiR61OdQ7w05crIzvAJ2tNLFc69i1VSmrxI4bGqesLfZRI9/lPfzj7FHauLGLnyuJMlUg5kXNGrdEcmBO50KNPpFNpupOS3eO0L7pePS6B9NsupkO7T2QvJXJKfSJdGPhZ2za1K50nwSTC2Ufbti9jtd7sGdWIw1rdaxfVAH5hzaqajYsh6VRn946YFQpEuchUCwtm1lYid6ws4rEnT6M5pShJN3Ii54woR6qbXtXZ7W7/XUqkUybTXuHcq7EyICVS+PTqEwn4Dfq9KU2sObZWx+ZqCWedUR156ks/JlFY4xzg87ZtAjDaxJ2TNQ9LISdyqVKUEimGZtDYX0e1VEy1sPDUKQ+r9SZ2rCxgx8oivJbhyInaTGyREzlnxKnOXgiqs7snVriWHz1DwiPO2500vXpcAlIihU+7T2REgVi5wKnNWD+2WsfKchkrSxUcW6uPND2mH9FK5GhrwoXinRM5Sjh+rd7EpmpnzS6rxY8YgXZhzUAnMt0dRw4cXwOAdjgbwMzGH8qJnDPi5EQu9Gi26pTGniHhtCuRPdoTueeEaFdn91Iip5YT2cDWpQq2LlXQaNrEHKd+HQuGzol0SuSZy/7jEZTU1ZrXnlQD+K1+pESKYanFVCLTPoXNtffZsbKInVvkRIopEpVI381ioESc6so96hUKyMr86cjZ2RlJohbTwetTnV0ucmp9Io+t1rFluYIty5Xg8WTyImtec0O7k1GLzY6u1VEpFbBzZcl/PIISudqVE7lcLWFGqV8iB8TpEwn4F1JpFkVcIc3OLYvY/rSFdcemjZzIOSPKkerGOZGnu5SIXldxWalwjpxYk5HGsmI69K3OnqISeWytji1LFWxZKrcfT4KoBvyj5gkfW637KqpzgEdSIpvtRuMA1qmSQiSl5jVRKjAy0hAm/UrkKVRKBTx9uYLNC2WcsVBKrxNJ8lqSh0l+L3RsK8mbSd4f3G4JPfdukvtJ3kfylZMyXAxHLcKR6max0l+JzGpeYeTEmoy0JxLTweU8VkrR1dnTzIncstRRIifVcDzq+2DUPpFHVxvYslzBSuAAj6JEntwQzpYTKYYnKgc4imqpmGol8kDQ3sf1T90xwzY/cZTIjwG4rOvYuwDcYmbnA7gleAySFwK4AsBzg9e8n2R6u0/PIXGUSLeJdFdoD1Ii0+5ERimRpWIBxQJTb7uYDr0m1rhj0+gTWfOaWK03sXW5jK1LLpw9fifSzPw10d0ncsRCuWNrdWxdLqNcLOCMhdLQtjeaLdS91vpwdiX5dkJyN8mvkLyH5F0k3xEc/z2Sj5LcF/y8OvQaiSE5JKobQRRpL6zx2/sstB/v2rKIA8dS6kSa2dcAHO06fDmA64L71wF4Xej49WZWM7OHAOwHcMmYbBVjoBbjSswpkd393dyi6q3mpXfRNVsGr2WoFDduQpXibL8wtMmlB69Pn8jylPpEupY4K0sVbHFO5IgzqKOYVKHcsbU6VgK7ty5XhrZ9LRhvuFRZX509BB6A3zSzHwfwIgBvCwQPAPgTM7so+LkJkBiSZ2qNeEpkFsLZriobmK0SOWxs4CwzOwQAZnaI5JnB8Z0Avh4670BwTKQAr9lCs2WxC2tOdy2idji7R15hmkPCvWx3x2asRLpN7lskNwO4g+TNwXN/YmZ/FD65a5PbAeDLJJ9tZun14jNCo9WnT2RhOn0iXfh363IFmxdKKBY4ESWy03i5R07kkJuoy4kEgC3LlaFzIleD8YabqqOFs4O9yu1XJ0jeg/77UlsMAfAQSSeG3Jb4l4tUEdXmLYpqqZDaVlJ1r4XDJ2rY0eVEPnXaw4nTDWxeKE/VnnEX1my8fAciL91JXk1yL8m9R44cGbMZIop2s/CYLX56VmdnsMK5X1VepTjbSjwzO2Rm3wrunwAQe5OT4j9e2kpkRJ/IaeVEOodxy1IFhQKxZak8kZxI5yR2O5EkA3U++ZpotgzHTzXauZxblypD50S6Vj5L1fUtfkaB5DkALgbwjeDQ20neGeT+u9z+nQAeCb2spxiifSxb+N1Jsq1EPvbkaZhhgxMJAIeePD11e4Z1Ih8nuR0AgtvDwfEDAHaHztsF4GDUG5jZh8xsj5nt2bZt25BmiCTUeygP3bRb/PTIiczi1Jd2KD7ib0/TF8a4NzmRjH59IstTmljjHMYty76isLJUmYgS2e+isloaTp1/8lQDZsDWpdFtd+MNN42pOpvkJgCfAfBOM3sKwAcAnAfgIvhK5R+7UyNeHnn1oH0sW+ShsMb1gwyHs2fZcHxYJ/JGAFcG968E8LnQ8StIVkmeC+B8ALePZqIYF3FHPi2Uo3Mis1yd3csBdsdqKfjCGPcmJ5UkOW5iTVSfyFJhOjmRLofQhYS3Lg0fEu6HK5yJCu9Vy8PlCTvVsa1ELg+vojolchzV2STL8NfWJ8zsrwDAzB43s6aZtQB8GB01P7YYIrJFosKalE5gc7mPO6KcyBkU18Rp8fMp+LkgF5A8QPIqAO8B8AqS9wN4RfAYZnYXgBsA3A3giwDepjyt9NDPkQrT04nsoVxURkzEnwb9JhVUSsWZK5GT2OSkkiTH61OdXS5OpzrbKXcr7bzC8kSajfddE0OGs52z63pEblmu4HSjtSE1Jg5tJ3LEcDb9PigfBXCPmb03dHx76LTXA3Bt7CSG5JREhTUpFUWcE+majAPAts1VlAqcSXHNwMs6M3tjj6cu7XH+NQCuGcUoMRnijnxq94nsDmcPnJ2dzkUH9FZR/WOzzYnst8m5AjZs3OQ+SfK98AtrtMmNicaAiTXeFMalHF2tY3O11F5XW5cr+NYPj4/99/RLb6mWi0NFFo6G8jmBjpp6dK2OnZXFnq+LwhXWhJ3ISrEQma86gBcDeBOA75LcFxz7twDeSPIi+Cr+wwB+FfDFEJJODPEgMSQ31JotrFQGF55US8N9/qfBo8dP4Rmbqm2xBwCKBeKZT1tIpxMp8kNn+Hz/q/mFYFPZUJ3doyVIsUCUCkS9md7v2V62A74TXJ9teyJtcinBa7VQLLDdxDeMPzt7Gi1+6lhZ7mx0K0sVHF+rw8wi7RqW2oALq2HC2cdWNyqR7ng4hysOJ4MWP+HekCQTh7TN7B8QnQJyU5/XSAzJIbVGE9XN1YHnpVmJfPT4KewM9Yh0+G1+pl9YIydyjohbWFMqFlAucoMS2S+nsjJkIv606KWiArNv56BNLj14TYtUIQFfiZyGYn10rdFW8ABfzWs0DSdr3ljbd0xiPbeLgkJ9IoHhptasRYSzgeEajgsB+GJCvMIaPzo17gu3cXDw+Ck8+6zNG47vWlnENx7qbuk9eTQ7e46ImxMJ+HmR3XlMNa8FMrr9SZqv3ID+SuSwlagifzSahnJEPiQQ9ImcUk6kU/CAsJo33rzITmSiR7HZEGvi+FoDC+VCOyWm0yw9uRPpciIXy+udRo0+FMNSa8TsE5nS3sdmhoPHT68rqnHsWFnEY0+dnsp3VBg5kXNE3OpswP/i7g5n1YPKtqgrs7Q7Yi5fs2eLnxTbLqZHo9mKnFYD+H0ip1OdXW87XwCwJWiXM+4K7fZFZUQD/mHXxNFQo3Ggo0QO0+Zntd7EcqWIQtdF65KcSDEkcSa2AeEpbOnaF46tNXCq0YxMDdmxsohmy3D4RG2qNsmJnCP65UB100uJ7DV3O+3h7I4SGZX/ld4kajFdvFYrskckEFRnT6FP5LHVLidyuVOcMk76NeAfdk10q6hPWyyD9EP0SVmteZEO46YRG46L+aUes9l4tZzOtnVR7X0cbpb2tItr5ETOEUmVyKhm45UeDuiwLUGmRd/8r2K6HWAxPfxwdg8lcgp9Ik83mlitN7E1VFjjlL1xNxxvh7PL4yusObpWb6uPgF90t7JYHlqJ3BThRI7ScFzMNzWvFam8d9MZoJGuesWoRuOOXVtm03BcTuQc0S8HqpuFSjFydnav11ZLxVQ7kX3zv4ZsrNwLs8mHPMVk8Jq9lchSsQCvZRP9/z0eKHYr68LZwxen9KNfjvSw4exuFRXwldRhVNTVmhfZF1I5kWIYzMyfnd1jfYdxayJt4kJHidxYnb39aXIixYRJokQulAoRSmTvUEBlSOViWkxTiZzGfGUxGRot65kT6RTKSfaK7G7WDQCbF0ooFth2MMdFvzUxbI7z0dX1SiQQTNwZRomseViKUB1HnZ8t5pNG02AWrbx3k9ZRvo8eO4WFcmHDGgP8i6uVpbLC2WJyxG02DvgNx6PGHvZ6bdpzIvtPrBlvKH4a85XFZPCard7V2YGCMcmQ9rGuZt0AUCgQW5aGHx/Yi/6jQJNHFhrNFp467W1QIleWKkOpqKt1T+FsMTacyNErrz9MJaVO5MEnT2HHymLPtkM7njb9XpFyIueIuH0igSAnsquwpt7sF86e7dSXQQzaML2WoTkmhUlKZHbxmn2UyOD4JItrjkYokYDvVI49J7LRe1OtDDE72Cml4XxO93i4Fj9NLEX0hFQ4WwxDe/+LlROZzsKaR4+f7tu0f+eWRSmRYnK0r8Ri9ok83RWe7jd31B9Yn64FF6bWpxK1Mub8l2n36RLjww9n967OBiasRAaOmGvr49iyVBl/i5+g8XLPll0JP8fOvi3dDvByBcdWG4lzSVdr0UpklGMpxCCS9ElOa2HNweOnsONpfZzIlUU8ekxOpJgQ/Vp6dOO3+Nk49rBvODvFzpMLxffaMN0542Aa85XFZPCaLVT69Il050wKpzaubChOKY+/2Xijf6Fco5lMne+em+3YulRBvdnCaj3ZhtwrJzLKsRRiEElqAsYtLIyD040mjpyoRbb3cexYWcCJmoenTo/3u6IfciLniFofR6qbxXJ0TmSvHpNp77VY93pX5bXzX8Y0+7uRYmda9MdrGko9ciJdruQkL5aOrtaxuVrasNFtHbLCuR+1Pt0WhtlEj69FO5Fbhmg43moZ1hrNyJ6QajYuhiFJn2R3TppyIh970s913LmlnxPpLo84GwAADrZJREFUPzfNkLacyDmi36bRzWKlsMGJrHnN3s3Gi+muzq55zb4qKoCxheOnMdVETIb6gIk1wGT/f4+v1bGyvHE+9spSBcfX6mNtL9T/ojC5E3l01eVEblQigWQTd041mjCLdhjVbFwMw3CFNenZ0/q193HsHJMTmUQIkRM5RyRxIheCYpPwh6nep1FrtZzu6uz+PS7HqzCpOju7eK1WO/exm3Z19kQLaxrrxgY6ti5V0GgaTgbzpMdBv5ZdndnB8TdR5ySudOdzLifvc7la9//OqCIaVWeLYUhWWJO+cHa/RuMO99yjI1Zo37jvYOxz5UTOEf2Uh24Wg+T1cK/Ier+xhymf+tJvZmp1zEpk3ZMSmVX8cHb/PpGTrL7vHhvo6ISEx5fr1Ldl1xCzg4+u1rFcKWKhqw9fe352AiVyteZ/70SpjqrOFsOQLJydvhY/rnXPM5/WW4l8xqYqykWOVFzTahn+29ceiH2+nMg5ol9hTDduIwiHtPs5YuPutThuBuVzAlIihR/GGahETtCJPBox8QXoVGuPMy+yX2TCNWROsqZ7OcBb2xN34jvAq4HiGt1sXE6kSE7WC2sePb6GbZurfZ3gQoHY/rTR2vx85b7D+P7jJ2OfLydyjqg14g2fB0JOZL0rnN3HEfNahlZKK5MHVZYD4/vCUJ/I7OL1mVhTmkKfyONrPZzIIdS8Qfjh7Oj13FEi44ezu+dmO9zEnSSFNc6JjG42rpxIkZwkY3/TWFhzcECPSMeOlYWRnMgPfvWBWL/HISdyjkiiRC6WN4aza3EcsZRWJscqrBlTErX6RGaXONXZk1IiTzeaWK03NzTrBkLFKWNsON4vnO3yxpJcWEXNzQaGm7jjciLVbFyMiyQT29xggaQN9yfJweOnYjl3O1eWhnYi7/jBUXzz4WO46mfPjf0aOZFzRL++cN0sVvzzXDjbzAaOPXS/I43EKqxRn8i5xw9nz6ZPpJv40i8ncpjxgb3oG84eIieslxIJ+NXlyZRIlxO50WFcjDH7WIhukjQbJ+kP0EiJIGBmePT4qb6V2Y6dKwt47KnTQ7Wa+8CtD2JlqYwrLtkd+zVyIueIRDmRpfVKpFMYB246Y+q1OG4G5XO6c8aB+kRml37h7M7Yw8lcJBzr0WcRAM4IQsLO0RwHsYrNEimRjQ2V2Y6tCednt3MiI5zIQo/CJyH6kaSwxj8vPVPYnlito+a1+jYad+xYWUTLgMefSlahff/jJ/Dlex7Hm3/6nMhc5F7IiZwj+uVAdbPQVZ09aO52GhORw/RXIsc7J1V9IrNLv8KaztjDyXzGj/WY+AL4ykjSkPAgxrkm6l4LJ2teZHsiwJ+4k8QBdtNtNqmdjxgTSQpr/POKqUnPOhijvY+j03A8mRP53772IBbKBbzlZ85J9LqJOZEkLyN5H8n9JN81qd8j4tOvRU83LmTkckIG5ZOksSVCGD8U36OIYOzh7Mn/G2h9TQavab2rs4OcyEkVTjkHsVdIeEvCkPAg+l1UJp0dfLzH3GxH0ok7HSVydqFrrbF8kaSwxp2XFiWy02g8Rk7kluQNxw89eQqf2/co/sWe3T2/f3oxESeSZBHA/wfgVQAuBPBGkhdO4neJ+PQLX3WzUE6mRKaxOWuYePlf4wnF1yesRGp9TY5Gs9W7T6QLZ09KiXQ5kT1CwluWKmOuzh5fjnNcBzjuxJ3VuodKqdDToZ80WmP5w32W4wop1VJ6prC55uGxlMinuYbj8Z3Ij/79Q2gZ8C9/7lmJbZtUrOASAPvN7EEAIHk9gMsB3B118vcfP4FXvPerEzJFOB49dgo/ec7WWOc6JfI/3nQv3v+VB9ob56BN5+r/vredT5mEPzpwHADwWxP6HDz21OmBtv/XW/bj47f9YOTf9dTp8eWt9SDR+gK0xuLi50T27xP5B39zN953y/2J3jfO57sz8aV3SPgr9x0Z2//jk6cafSIL/hr+T1+8Fx/86uDGw+5is2dO5HIFXsvw8vd+FQUOzmk8fKI261Y+WmM540cnaygXGTuntlIq4NaE621S+9gTq3Uslos911eYxUoRW5cr+PDfP4i//vajsd7/B0+s4TU/sR27ty4ltm1STuROAI+EHh8A8FPhE0heDeBqADhjx7Nw/lmbJmSKcDz7rM345y/cFevcMzdX8dYXn7MuOffis7fgp5/1jMjzX3D2FvwfL9iFU43hxrI5p3VSn4Nnn7UZr794Z+Rzm6sl/Oo/eRYeObo2tt93+9jeKZKB6wvQGhuGC565GZc995mRz+3esog3vejH8MRqLfH7xv18n3/m5p6O3S/91I+hOMaikmc/czNe+/wdkc9Frf9B/Mx5T8dFu1cin3v5j5+FOw88GTvV4/yzNuEFZ2+J/bsngNZYzjj/rE24cPsZsc+/6mfPxVfuO5zod0xqHzsfwMW7t4AxLsAA4B2Xno9vPPRE7Pf/8e1n4Dde8eyhbGPc8EKiNyXfAOCVZvYvg8dvAnCJmf2rqPP37Nlje/fuHbsdIkO89KX+7a23ztKKsUHyDjPbM6H3TrS+AK2xmZOzz/esmeT6Ct5fa0wkJ0frPO4am1TCyQEA4UZDuwDEn+gthOiH1pcQk0VrTIgYTMqJ/CaA80meS7IC4AoAN07odwkxb2h9CTFZtMaEiMFEciLNzCP5dgB/C6AI4Fozu2sSv0uIeUPrS4jJojUmRDwm1snVzG4CcNOk3l+IeUbrS4jJojUmxGA0sUYIIYQQQiRGTqQQQgghhEiMnEghhBBCCJEYOZFCCCGEECIxE2k2ntgI8gSA+2ZtR4hnAPjRrI3oIm02yZ7+XGBmm2dthENrbCCypz9psydV6wvQGouB7OlP2uyJtcYmVp2dkPsmOX0gKST3pskeIH02yZ7+kEzb6AqtsT7Inv6k0Z5Z2xCB1lgfZE9/0mhPnPMUzhZCCCGEEImREymEEEIIIRKTFifyQ7M2oIu02QOkzybZ0x/Z0x/Z0x/Z05+02QOkzybZ0x/Z059Y9qSisEYIIYQQQmSLtCiRQgghhBAiQ8iJFEIIIYQQiZm5E0nyMpL3kdxP8l0ztuVakodJfm+WdjhI7ib5FZL3kLyL5DtmbM8CydtJfiew59/P0h4HySLJb5P8fApseZjkd0nuS0MbkjStr8AerbH+9miNDbZFa6y/PalZY1pf8UjT+gKSrbGZ5kSSLAL4PoBXADgA4JsA3mhmd8/InpcAOAng42b2vFnY0GXPdgDbzexbJDcDuAPA62b470MAy2Z2kmQZwD8AeIeZfX0W9oTs+g0AewCcYWavmbEtDwPYY2YzbxqbtvUV2KQ11t8erbHBtjwMrbF+NqVmjWl9xbYrNesrsOdhxFxjs1YiLwGw38weNLM6gOsBXD4rY8zsawCOzur3d2Nmh8zsW8H9EwDuAbBzhvaYmZ0MHpaDn5lWZpHcBeCfAvjILO1IKalaX4DWWAx7tMayhdZYH7S+BpP19TVrJ3IngEdCjw9ghh+wNEPyHAAXA/jGjO0oktwH4DCAm81spvYA+C8AfhtAa8Z2OAzAl0jeQfLqGdui9ZUArbGeaI31RmssJlpfPUnb+gISrLFZO5GMOKaeQ12Q3ATgMwDeaWZPzdIWM2ua2UUAdgG4hOTMwiUkXwPgsJndMSsbInixmb0AwKsAvC0ILc0Kra+YaI1FozU2EK2xGGh9RZPS9QUkWGOzdiIPANgderwLwMEZ2ZJKgryNzwD4hJn91aztcZjZcQC3Arhshma8GMBrg/yN6wG8jORfzNAemNnB4PYwgM/CD3fNCq2vGGiN9UVrrD9aYwPQ+upL6tYXkGyNzdqJ/CaA80meS7IC4AoAN87YptQQJAF/FMA9ZvbeFNizjeRKcH8RwMsB3Dsre8zs3Wa2y8zOgf/Z+Tsz++VZ2UNyOUgeB8llAL8AYJYVklpfA9Aa64/W2EC0xvqg9dWftK0vIPkam6kTaWYegLcD+Fv4Cbc3mNlds7KH5KcA3AbgApIHSF41K1sCXgzgTfCvTvYFP6+eoT3bAXyF5J3wvzxvNrNUtCRICWcB+AeS3wFwO4C/MbMvzsqYtK0vQGssBlpj/dEaG0DK1pjWV/ZItMY09lAIIYQQQiRm1uFsIYQQQgiRQeRECiGEEEKIxMiJFEIIIYQQiZETKYQQQgghEiMnUgghhBBCJKY0awOEECINkPw9ACcBnAHga2b25R7nvQ7A983s7imaJ4QQqUNKpBBChDCz/7eXAxnwOgAXTsseIYRIK3IihRBzC8nfIXkfyS8DuCA49jGS/zy4/x6Sd5O8k+QfkfwZAK8F8IdB4+TzSP6fJL9J8jskP0NyKfQ+7yP5v0g+6N4zeO63SX43eM17gmPnkfwiyTtI/j3J50z9H0QIIRKgcLYQYi4h+UL4o8Yuhv9d+C0Ad4Se3wrg9QCeY2ZGcsXMjpO8EcDnzezTwXnHzezDwf0/AHAVgD8N3mY7gJ8F8Bz44/A+TfJV8NXMnzKzteD3AMCHAPyamd1P8qcAvB/Ayyb4TyCEECMhJ1IIMa/8HIDPmtkaAATOYZinAJwG8BGSfwOg13i05wXO4wqATfBH4Dn+2sxaAO4meVZw7OUA/tz9XjM7SnITgJ8B8D/8ccMAgOpIf50QQkwYOZFCiHmm59xXM/NIXgLgUviK5dsRrQx+DMDrzOw7JN8C4KWh52qh+wzddv/eAoDjZnZREuOFEGKWKCdSCDGvfA3A60kuktwM4J+FnwzUwaeZ2U0A3gnAOXgnAGwOnboZwCGSZQC/FOP3fgnAr4RyJ7ea2VMAHiL5huAYST5/hL9NCCEmjpxIIcRcYmbfAvCXAPYB+AyAv+86ZTOAz5O8E8BXAfzr4Pj1AP4NyW+TPA/A/wPgGwBuBnBvjN/7Rfj5kXtJ7gPwW8FTvwTgKpLfAXAXgMtH+POEEGLi0KxnNEcIIYQQQohIpEQKIYQQQojEyIkUQgghhBCJkRMphBBCCCESIydSCCGEEEIkRk6kEEIIIYRIjJxIIYQQQgiRGDmRQgghhBAiMf8/YYeulhvHgC4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 792x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(11,4))\n",
    "ax1.plot(bccrdf[1], bccrdf[0])\n",
    "ax2.plot(fccrdf[1], fccrdf[0])\n",
    "ax3.plot(hcprdf[1], hcprdf[0])\n",
    "ax1.set_xlim(0,5)\n",
    "ax2.set_xlim(0,5)\n",
    "ax3.set_xlim(0,5)\n",
    "ax1.set_title('bcc')\n",
    "ax2.set_title('fcc')\n",
    "ax3.set_title('hcp')\n",
    "ax2.set_xlabel(\"distance\")\n",
    "ax1.axvline(3.6, color='red')\n",
    "ax2.axvline(2.7, color='red')\n",
    "ax3.axvline(3.6, color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The selected cutoff distances are marked in red in the above plot. For bcc, since the first two shells are close to each other, for this example, we will take the cutoff in such a way that both shells are included."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Steinhardt's parameters - cutoff neighbor method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcc.find_neighbors(method='cutoff', cutoff=3.6)\n",
    "fcc.find_neighbors(method='cutoff', cutoff=2.7)\n",
    "hcp.find_neighbors(method='cutoff', cutoff=3.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have used a cutoff of 3 here, but this is a parameter that has to be tuned. Using a different cutoff for each structure is possible, but it would complicate the method if the system has a mix of structures. Now we can calculate the $q_4$ and $q_6$ distributions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcc.calculate_q([4,6])\n",
    "fcc.calculate_q([4,6])\n",
    "hcp.calculate_q([4,6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thats it! Now lets gather the results and plot them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "bccq = bcc.get_qvals([4, 6])\n",
    "fccq = fcc.get_qvals([4, 6])\n",
    "hcpq = hcp.get_qvals([4, 6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f5cfc4edc10>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEPCAYAAACDTflkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAbtUlEQVR4nO3dfXRV9Z3v8feXx8RYdC4QiQkTLCMRndUqHnNFqWZU1AuMT5cRxtYW5zqgrYMPtJbaK3i9Klxv60ir1rLAB1qV8cZqEbE+DLqWC71KCGqNEcvFAQPRE0IhNYbykO/945ykIZyE8EvO2Sfk81rrrJyzf3uffAhr55P9cPY2d0dERCREv6gDiIhI76USERGRYCoREREJphIREZFgKhEREQk2IOoAmTZs2DAfNWpU1DFERHqNdevWbXf34anG+lyJjBo1ioqKiqhjiIj0Gma2uaMx7c4SEZFgKhEREQmmEhERkWAqERERCdbnDqyLiPQZK63jsSk9c91EbYmIiByJOiuQrox3kUpERESCqURERCSYSkRERIKpREREJJhKREREgqlEREQkmEpERORIdKjPgfTQ50T0YUMRkSNVDxVFZ7QlIiIiwVQiIiISTCUiIiLBVCIiIhJMJSIiIsFUIiIiEkwlIiIiwVQiIiISTCUiIiLBVCIiIhJMJSIiIsFUIiIiEkwlIiIiwVQiIiISTCUiIiLBVCIiIhJMJSIiIsFUIiIiEkwlIiIiwVQiIiISTCUiIiLBVCIiIhJMJSIiIsEiLxEzu9jMNpjZRjObm2K8zMx2mdm7yce8NmPHmlm5mX1kZtVmNj6z6UVE+rYBUX5zM+sPPAhMBGqAtWa2wt0/bDfrG+4+JcVbLAJ+5+5TzWwQcFR6E4uISFtRb4mUAhvdfZO77wGWA5d2ZUEzGwKcAywFcPc97r4zbUlFROQgUZdIIfBpm9c1yWntjTez98zsRTM7JTntq0Ad8KiZrTezJWaWl+qbmNlMM6sws4q6uroe/QeIiPRlUZeIpZjm7V5XAsXu/nXg58BzyekDgHHAL9z9NKAROOiYCoC7L3b3mLvHhg8f3jPJRUQk8hKpAUa2eV0EbGs7g7s3uPsXyeergIFmNiy5bI27v52ctZxEqYiISIZEXSJrgRPN7ITkgfHpwIq2M5jZCDOz5PNSEpnr3f0z4FMzK0nOej7Q/oC8iIikUaRnZ7n7PjO7AXgJ6A884u5VZnZdcvxhYCpwvZntA5qA6e7essvrX4AnkgW0Cbgm4/8IEZE+zP7y+7hviMViXlFREXUMEZFew8zWuXss1VjUu7NERKQXU4mIiEgwlYiIiARTiYiISDCViIiIBFOJiIhIMJWIiIgEU4mIiEgwlYiIiARTiYiISDCViIiIBFOJiIhIMJWIiIgEU4mIiEgwlYiIiARTiYiISDCViIiIBFOJiIhIMJWIiIgEU4mIiEgwlYiIiARTiYiISDCViIiIBFOJiIhIMJWIiIgEU4mIiEgwlYiIiARTiYiISDCViIiIBFOJiIhIMJWIiIgEU4mIiEgwlYiIiARTiYiISDCViIiIBIu8RMzsYjPbYGYbzWxuivEyM9tlZu8mH/Pajfc3s/VmtjJzqUVEBGBAlN/czPoDDwITgRpgrZmtcPcP2836hrtP6eBtbgSqgSHpSyoiIqlEvSVSCmx0903uvgdYDlza1YXNrAiYDCxJUz4REelE1CVSCHza5nVNclp7483sPTN70cxOaTP9fuBWoLmzb2JmM82swswq6urquh1aREQSoi4RSzHN272uBIrd/evAz4HnAMxsChB393WH+ibuvtjdY+4eGz58eHczi4hIUtQlUgOMbPO6CNjWdgZ3b3D3L5LPVwEDzWwYcDZwiZn9B4ndYOeZ2a8zklpERIDoS2QtcKKZnWBmg4DpwIq2M5jZCDOz5PNSEpnr3f1H7l7k7qOSy612929lNr6ISN8W6dlZ7r7PzG4AXgL6A4+4e5WZXZccfxiYClxvZvuAJmC6u7ff5SUiIhGwvvb7OBaLeUVFRdQxRER6DTNb5+6xVGNR784SEZFeTCUiIiLBVCIiIhJMJSIiIsFUIiIiEkwlIiIiwVQiIiISTCUiIiLBVCIiIhLssEvEzL5hZiVtXn/XzKrM7E9m9kHLJUtEROTIF7Il8gvgeAAz+x7wv0hcnv164FlgYXK6iIgc4UIuwDga2JR8fi1wg7s/3jJoZu8D/5PEbW9FROQIFrIl8idgWPL58cC77cYrgb/uTigREekdQkpkFXBD8vlrwJXtxqcBH3cnlIiI9A4hu7PmAmvM7A3gbeBmMzsHqAZKgDOBy3ouooiIZKvD3hJx98+AccAbwH8hcZ/0UuBCEre7PdvdX+zJkCIikp2C7mzo7ruA25IPERHpow5ZImb2AbCexAHz9cD6ZImIiEgf15UtkZHAycBVLRPM7BMOLJZKd4+nJaGIiGStrpTICKCcxPGPTUA9UAj8V+CKlpnMrJaDi2VLTwcWEZHs0ZUSuRc4FzjX3d9omWhmfwvMA6YCHwL5wOTkA6C5i+8vIiK9VFd+yU8FlrctEAB3/wC40szmAd8C/gYYApxG4uytU3s4q4iIZJmulMixwPaOBt39TjO7DLjN3eeSOM33+R7KJyIiWawrnxPZAPzdIeZ5FfiH7scREZHepCsl8ghwhpnd3Mk8/4nklX1FRKTv6EqJPAi8AvzEzFaZ2ZltB83sXOAfgW1pyCciIlnskMdE3H2/mU0hUSbXAheZWT2wBRhK4oq9Bvw8nUFFRCT7dOnaWe6+191nAv8Z+DWwn8QZWEUkjpn8k7vfn7aUIiKSlQ7rcxzuvhb4DoCZDQb2u/u+dAQTEZHsF/xhQHf/c08GERGR3ifkplQiIiKASkRERLpBJSIiIsFUIiIiEkwlIiIiwSIvETO72Mw2mNlGM5ubYrzMzHaZ2bvJx7zk9JFm9pqZVZtZlZndmPn0IiJ9W6T3+zCz/iQ+CT+RxNV/15rZCnf/sN2sb7j7lHbT9gFz3L3SzL4CrDOzV1IsKyIiaRL1lkgpsNHdN7n7HmA5cGlXFnT3WnevTD7/E1BN4o6LIiKSIVGXSCHwaZvXNaQugvFm9p6ZvWhmp7QfNLNRJG6G9Xaqb2JmM82swswq6urqup9aRESA6EvEUkzzdq8rgWJ3/zqJizw+d8AbmB0NPAPc5O4Nqb6Juy9295i7x4YPH94DsUVEBKIvkRpgZJvXRbS7pLy7N7j7F8nnq4CBZjYMwMwGkiiQJ9z9N5mJLCIiLaIukbXAiWZ2gpkNAqYDK9rOYGYjzMySz0tJZK5PTlsKVLv7fRnOLSIiRHx2lrvvM7MbgJeA/sAj7l5lZtclxx8GpgLXm9k+oAmY7u5uZhOAq4Hfm9m7ybe8Lbm1IiIiGWDu7Q9BHNlisZhXVFREHUNEpNcws3XuHks1FvXuLBER6cVUIiIiEizSYyIiIj2hoaGBeDzO3r17o47S6wwcOJD8/HyGDBkStLxKRER6tYaGBj7//HMKCwvJzc0leTKndIG709TUxNatWwGCikS7s0SkV4vH4xQWFnLUUUepQA6TmXHUUUdRWFhIPB4Peg+ViIj0anv37iU3NzfqGL1abm5u8K5AlYiI9HraAume7vz8VCIiIhJMJSIiIsFUIiIiWWDGjBnEYik/FJ7VVCIiIhJMnxM5hLemToWqqoMHTjmF8eXlmQ8kIpJFtCXSiQ4LBKCqKjEuIkcEb26m7vnneX/qVComTOD9qVOpe/55vLk5ozmee+45TjrpJHJycpgwYQIffvhh69j+/ftZsGABY8aMYfDgwRQVFTFjxowDln/22WcpLS0lNzeXoUOHMmnSJDZv3py2vCqRznRUIF0dF5FewZub2TB7Npvmz6exqoq99fU0VlWxaf58Ntx4Y8aKZPPmzdxyyy3cfvvtPPnkk+zatYuLLrqI3bt3AzBr1izmz5/PlVdeycqVK/npT39KY2Nj6/K/+tWvuOKKKxg9ejRPP/00jz76KGPGjCGdtwXX7iwR6fO2v/ACu958k+ampgOmNzc1sWvNGravWsXwKVPSn2P7dn77299y1llnAXD66aczevRoHnvsMcrKyli6dCmLFi1i9uzZrctMmzYtkbW5mblz53L55Zfz1FNPtY5fcsklac2sEhGRPq/28ccPKpAWzU1N1D72WEZKJD8/v7VAAIqLizn99NN55513aLn3U/vdVy02bNjAtm3buOaaa9Kesy3tzhKRPm/PZ591a7yn5Ofnp5xWW1tLfX09eXl5HV4ksb6+HoCCgoK0ZmxPJSIifd6gESO6Nd5TUl0EMR6PU1BQwNChQ2lsbKShoSHlskOHDgWgtrY2rRnbU4l05pRTujcuIr1CwXe+Q78OLuLYLzeXgg52IfW0eDzOm2++2fp6y5YtVFZWUlpaynnnnQfAsmXLUi5bUlJCYWEhjz/+eEaytlCJdGJ8eXnHRaHPiYgcMYZNnswxZ511UJH0y83lmLPPZtikSZnJMWwYV199NU8++STPPvsskydPJj8/nxkzZlBSUsLMmTOZM2cO8+bN49VXX6W8vJzp06cnsvbrx7333sszzzzDN7/5TVauXMkLL7zAnDlzqKioSFtmHVg/BBWFyJHP+vWj5Gc/Y/uqVdQ+9hh7PvuMQSNGUDBjBsMmTcL6Zebv7eLiYm677Tbmzp3L5s2bicViPPXUU+Tk5ADw0EMPUVxczJIlS1i4cCH5+flMnDixdfmrrrqKnJwc7r77bqZOnUpeXh5nnnkmw4cPT1tmazni31fEYjFPZyuLSGZVV1czduzYqGP0ep39HM1snbunvLCXdmeJiEgwlYiIiARTiYiISDCViIiIBFOJiIhIMJWIiIgEU4mIiEgwlYiIiARTiYiISDCViIiIBFOJiIhkiTvvvJPCwkL69evX4c2nso0uwCgikgUqKiqYP38+99xzD2VlZSlvUJWNVCIiIlngo48+AuB73/teh3cvzEaR784ys4vNbIOZbTSzuSnGy8xsl5m9m3zM6+qyIiJd5s1Q8wS8EYOXj0t8rXkiMT3NZsyYwdVXXw3AMcccg5nx+uuvU19fz6xZsygoKCAnJ4eSkhLuv//+1uX279/PggULGDNmDIMHD6aoqCjju8Ei3RIxs/7Ag8BEoAZYa2Yr3P3DdrO+4e5TApcVEemcN0PFFbD9VdjfmJi2Jw6/nwW15RB7Bix9f3PffvvtjBw5krvuuovVq1eTm5vL2LFjmTBhAvF4nPnz53PSSSexceNGNm7c2LrcrFmzWLZsGbfeeivnnnsuO3bsoDzD90CKendWKbDR3TcBmNly4FKgK0XQnWVFRP5i61MHFkiL/Y2w/RXYthwKr0rbtx89ejSjR48G4IwzzuDoo4/ml7/8JVVVVVRWVnLqqacCtN4iFxK7v5YuXcqiRYuYPXt26/Rp06alLWcqUe/OKgQ+bfO6JjmtvfFm9p6ZvWhmLfer7eqymNlMM6sws4q6urqeyC0iR5JP/vXgAmmxvxE23ZfZPMDq1as57bTTWgukvddeew0g8rO4oi4RSzGt/a0WK4Fid/868HPgucNYNjHRfbG7x9w9ls7bRIpIL9X06SHGazKTo436+noKCgo6Hc/Ly4v8IHzUJVIDjGzzugjY1nYGd29w9y+Sz1cBA81sWFeWFRHpktyRhxgvykyONoYOHUptbW2n442NjTQ0NGQw1cGiLpG1wIlmdoKZDQKmAyvazmBmI8zMks9LSWSu78qyIiJdcsLN0D8v9Vj/PPjqLZnNA5x//vmsX7+e999/P+V4y/GRZcuWZTLWQSI9sO7u+8zsBuAloD/wiLtXmdl1yfGHganA9Wa2D2gCpru7AymXjeQfIiK9W+E/Qu3/Ofjgev88GDYRjp+e8Ujf/va3efDBB7nwwgu54447KCkp4ZNPPuHjjz9m4cKFlJSUMHPmTObMmUM8Huecc85h586dlJeXs3z58ozljPrsrJZdVKvaTXu4zfMHgAe6uqyIyGGzfhD7TeIsrE33JY6B5BYltkCOn57W03s7kpOTw+rVq5k7dy7z5s2joaGBUaNG8d3vfrd1noceeoji4mKWLFnCwoULyc/PZ+LEiRnNaYk/6vuOWCzmFRUVUccQkR5SXV3N2LFjo47R63X2czSzde4eSzUW9TERERHpxVQiIiISTCUiIiLBVCIiIhJMJSIiIsFUIiIiEkwlIiIiwVQiIiISTCUiIiLBVCIiIllgxowZxGIpPxSe1VQiIiISTCUiIiLBVCIiIkCzN/NS9ctc88S1TPrlJVzzxLW8VP0yzd6c0RyvvPIKX/va18jLy2PChAlUVf3lDhf79+9nwYIFjBkzhsGDB1NUVHTA7XHLysqYOnUqixcvZtSoUeTm5jJ58mS2bt2atryRXwpeRCRqzd7Mj57/76zdspamvbsB+OOXf2Thq/+b1X94nQV/fxf9MnA5+C1btvCDH/yAH//4x+Tm5vL973+fK6+8kg8++AAzY9asWSxbtoxbb72Vc889lx07dlBeXn7Ae7z11lts2LCB++67j927d/PDH/6Qyy67jLVr16Yls0pERPq8Vz569YACabF7327e2byWVzb8OxedlP77dOzYsYM1a9Zw4oknAtDc3Mzll1/Ohg0bAFi6dCmLFi1i9uzZrctMmzbtgPeIx+O8+eabFBcXA1BcXMyECRP43e9+x8UXX9zjmbU7S0T6vOWVTx9UIC1279vN8nX/lpEco0aNai0QgJNPPhmAmpoaXnvtNYADdl+lMm7cuNYCATj77LPJz8/nnXfe6fnAqERERPj8i3i3xnvKsccee8DrQYMGAbB7927q6+vJy8tjyJAhnb5Hfn5+ymm1tbU9F7QNlYiI9HnHHX3wL97DGc+EoUOH0tjYSENDQ6fzxeMHF148HqegoCAtuVQiItLnTR93JTkDclKO5QzIYfrp01KOZdJ5550HwLJlyzqdr7Kyki1btrS+XrNmDfF4nNLS0rTk0oF1EenzJp50Aav/8DrvbF7L7n1/OTaSMyCH0uIzmFhyfoTpEkpKSpg5cyZz5swhHo9zzjnnsHPnTsrLy1m+fHnrfPn5+UyZMoU77rij9eyscePGpeWgOqhEREToZ/1Y8Pd38cqGf2f5un/j8y/iHHd0PtNPn8bEkvMzcnpvVzz00EMUFxezZMkSFi5cSH5+PhMnHnjW2Pjx47ngggu46aabqKuro6ysjMWLF6ctk7l72t48G8ViMa+oqIg6hoj0kOrqasaOHRt1jKxQVlbGsGHDDvrsSFd09nM0s3XunvLCXtlRryIi0itpd5b0aeP/9Rsdjr118xsZTCLSO6lEpM/qrEBaxlUk0pu8/vrrGf+e2p0lIiLBVCIi0uv1tROEelp3fn4qERHp1QYOHEhTU1PUMXq1pqYmBg4cGLSsSkREerX8/Hy2bt3Kl19+qS2Sw+TufPnll2zdujXlNbe6QgfWRaRXa7kg4bZt29i7d2/EaXqfgQMHctxxxx3ywo4dUYmISK83ZMiQ4F+C0j3anSV91qFO39XpvSKHpi0R6dNUFCLdoy0REREJphIREZFgKhEREQmmEhERkWB97n4iZlYHbE7DWw8DtqfhfXtCNmeD7M6nbOGyOV82Z4Psy1fs7sNTDfS5EkkXM6vo6KYtUcvmbJDd+ZQtXDbny+ZskP352tLuLBERCaYSERGRYCqRnrM46gCdyOZskN35lC1cNufL5myQ/fla6ZiIiIgE05aIiIgEU4mIiEgwlcghmNnFZrbBzDaa2dwU42ZmP0uOv29m49qN9zez9Wa2MtvymdmxZlZuZh+ZWbWZjc+ibDebWZWZfWBmT5lZTk9m62K+k8zsLTP7s5l9/3CWjSqbmY00s9eS/59VZnZjtmRrMx71OtHZ/2vU60Rn2dK+TgRxdz06eAD9gf8HfBUYBLwHnNxunknAi4ABZwJvtxu/BXgSWJlt+YDHgWuTzwcBx2ZDNqAQ+ATITb5+GpgRwc8uHzgDuBv4/uEsG2G2AmBc8vlXgI+zJVsWrRMd5suCdaKj/9e0rxOhD22JdK4U2Ojum9x9D7AcuLTdPJcCyzzh/wLHmlkBgJkVAZOBJdmWz8yGAOcASwHcfY+778yGbMmxAUCumQ0AjgK29WC2LuVz97i7rwXa3y6vK/+2SLK5e627Vyaf/wmoJvELKPJskB3rREf5smGd6OxnR/rXiSAqkc4VAp+2eV3DwStkZ/PcD9wKNGdhvq8CdcCjyV0LS8wsLxuyuftW4CfAFqAW2OXuL/dgtq7mS8eyGXt/MxsFnAa83SOpErqbLRvWiY5kwzqRUobWiSAqkc5Zimntz4lOOY+ZTQHi7r6u52N1/r27OM8AYBzwC3c/DWgEenLffnd+dn9F4i+0E4DjgTwz+1YPZutqvnQsm5H3N7OjgWeAm9y9oUdSJd86xbQuZcuidaIj2bBOpF4wM+tEEJVI52qAkW1eF3HwJmRH85wNXGJm/0Fis/U8M/t1FuWrAWrcveWv1HISK1A2ZLsA+MTd69x9L/Ab4KwezNbVfOlYNu3vb2YDSRTIE+7+mx7M1d1s2bJOdLZs1OtERzKxTgRRiXRuLXCimZ1gZoOA6cCKdvOsAL6dPNPoTBKbmbXu/iN3L3L3UcnlVrt7T//l0J18nwGfmllJcr7zgQ+zIRuJTfYzzewoM7NktuoezNbVfOlYNq3vn/x5LQWq3f2+HszU7WxZtE50lC8b1omOZGKdCBP1kf1sf5A4g+hjEmdV/Dg57TrguuRzAx5Mjv8eiKV4jzLScCZKd/MBpwIVwPvAc8BfZVG2/wF8BHwA/AoYHMHPbgSJvx4bgJ3J50M6WjYbsgETSOwieR94N/mYlA3Zsmid6Oz/Nep1orNsaV8nQh667ImIiATT7iwREQmmEhERkWAqERERCaYSERGRYCoREREJphIREZFgKhEREQmmEhHJIDMbYGazzew9M9ttZp+Z2QPJTyLvMrOe/IS0SNoNiDqASF+RvNTF88CFJD4V/TNgGPBPJK4gOwRIy42aRNJFJSKSOQ+QKJAfuPtPWiaa2ePA68mXlRHkEgmmy56IZICZnQG8A7zs7helGG+549357r460/lEQumYiEhm3JD8emcH4/XJr+vbD5jZbWbmZvZAWpKJdINKRCQzLgLq3X1NB+OFJO4X8ce2E5OXyP9nEleVFck6KhGRNDOzHOA4EveESDX+tyTuVre+3fRjgCeA/wb8McWiIpFTiYik3/7kY2gH4/OSX9sfVF8MlOsYiWQzlYhImnnidqZ/AP7azP6uZXryjo7zgH9ITlrfZuyfgb8Bbs9kVpHDpVN8RTLjXuAR4AUzewrYQeK+2V8hcQvWk0luiSRvz3oP8A133xNNXJGu0Sm+IhliZrcA/0Li+Mc2oBy4m8RWyl53Pz453wzgURK7wFr0J3Hb22Ygz93/nLnkIh1TiYhEyMxGkjjgvsrdJyenHQsUtZv1URJlcw9Q5VpxJUtod5ZItE5Lfm09qO7uO4GdbWcys0Zgh7t/kMFsIoekA+si0WopkYM+ZCjSG2hLRCRaB22JpOLuZemPInL4dExERESCaXeWiIgEU4mIiEgwlYiIiARTiYiISDCViIiIBFOJiIhIMJWIiIgE+//apPQ7FE5rhgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(bccq[0], bccq[1], s=60, label='bcc', color='#C62828')\n",
    "plt.scatter(fccq[0], fccq[1], s=60, label='fcc', color='#FFB300')\n",
    "plt.scatter(hcpq[0], hcpq[1], s=60, label='hcp', color='#388E3C')\n",
    "plt.xlabel(\"$q_4$\", fontsize=20)\n",
    "plt.ylabel(\"$q_6$\", fontsize=20)\n",
    "plt.legend(loc=4, fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, we can see that Steinhardt parameter values of all the atoms fall on one specific point which is due to the absence of thermal vibrations. Next, all the points are well separated and show good distinction. However, at finite temperatures, the atomic positions are affected by thermal vibrations and hence show a spread in the distribution. We will show the effect of thermal vibrations in the next example. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structures with thermal vibrations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again, we create the reqd structures using the :mod:`~pyscal.crystal_structures` module. Noise can be applied to atomic positions using the `noise` keyword as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcc_atoms, bcc_box = pcs.make_crystal('bcc', lattice_constant=3.147, repetitions=[10,10,10], noise=0.1)\n",
    "bcc = pc.System()\n",
    "bcc.atoms = bcc_atoms\n",
    "bcc.box = bcc_box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcc_atoms, fcc_box = pcs.make_crystal('fcc', lattice_constant=3.147, repetitions=[10,10,10], noise=0.1)\n",
    "fcc = pc.System()\n",
    "fcc.atoms = fcc_atoms\n",
    "fcc.box = fcc_box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "hcp_atoms, hcp_box = pcs.make_crystal('hcp', lattice_constant=3.147, repetitions=[10,10,10], noise=0.1)\n",
    "hcp = pc.System()\n",
    "hcp.atoms = hcp_atoms\n",
    "hcp.box = hcp_box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### cutoff method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcc.find_neighbors(method='cutoff', cutoff=3.6)\n",
    "fcc.find_neighbors(method='cutoff', cutoff=2.7)\n",
    "hcp.find_neighbors(method='cutoff', cutoff=3.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now, calculate $q_4$, $q_6$ parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcc.calculate_q([4,6])\n",
    "fcc.calculate_q([4,6])\n",
    "hcp.calculate_q([4,6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gather the q vales and plot them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "bccq = bcc.get_qvals([4, 6])\n",
    "fccq = fcc.get_qvals([4, 6])\n",
    "hcpq = hcp.get_qvals([4, 6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f5cfc743a50>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEPCAYAAAB/WNKuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOy9e3xU5b3/+37mkoQkEAi3AEJABEQoKoLRBqw7rbTbtl6qVWtrW3f9qW3d3XS3p/W3PVIPblvbo0fa/bNb1GqtrYp300qrtalboBpAUCDcA4RbbiSQK0nm8pw/nlkza61ZM5lMMplJeN6vF6+ZWbdZa4D1Wd+7kFKi0Wg0Gs1A40r3CWg0Go1meKIFRqPRaDQpQQuMRqPRaFKCFhiNRqPRpAQtMBqNRqNJCZ50n0CmMG7cODl9+vR0n4ZGo9EMKT788MMTUsrxTuu0wISYPn06mzdvTvdpaDQazZBCCFETa512kWk0Go0mJWiB0Wg0Gk1K0AKj0Wg0mpSgBUaj0Wg0KSHtAiOE+JwQYo8QYr8Q4m6H9ZcLIVqEEB+F/qwILZ9jWvaREKJVCLE8tO4+IcQx07orB/u6NBqN5kwnrVlkQgg38ChwBXAU2CSEKJdS7rRtuk5K+QXzAinlHuAC03GOAa+ZNnlESvlQyk5eo9FoNHFJtwVzMbBfSnlAStkDvABcncRxPg1USyljpstpNBpNr9SVw/a71Kum36RbYKYAR0yfj4aW2blUCPGxEOLPQoh5DutvAp63LbtLCLFNCPGUEGKM05cLIW4XQmwWQmxubGxM6gI0Gs0woa4ctn4Fah5Vr1pk+k26BUY4LLMPqNkCFEspzwf+C3jdcgAhsoCrgJdMi/8bmIlyodUCDzt9uZTycSnlIinlovHjHQtRNRrNmULj2xDoVO8Dneqzpl+kW2COAlNNn88Cjps3kFK2SinbQ+/XAl4hxDjTJv8MbJFS1pv2qZdSBqSUQeAJlCtOo9FoYjN+Gbhz1Xt3rvqs6RfpbhWzCZglhJiBCtLfBNxs3kAIUQTUSymlEOJilCg2mTb5Cjb3mBBikpSyNvTxWmBHis5fo9EMF4quggufV5bL+GXqs6ZfpFVgpJR+IcRdwFuAG3hKSlklhLgztP4x4Hrg20IIP3AauEmG5jwLIXJRGWh32A79CyHEBSh32yGH9RqNRhNN0VVaWAYQEbpXn/EsWrRI6maXGo0mirpyZdV4C8DXoq0bG0KID6WUi5zWpdtFptFoNJmFIShGDGbrVyLBf4CjTytXGmjh6QUtMBqNRmNgpCoHOpWQFF5uFRdQn3f8K/TUQ7A7stwQHi0yYbTAaDQajWG1dB60pip316PCwwHr9l2Ho49hpDZrgQmjBUaj0ZzZ1JXDli9DsAeVpOoCguDKgraPiRKXWLiydGqzDS0wGo3mzMSwWk5+EBIXgGDo1Q05U6GzOvHjBf3W4+qYjBYYjUZzBmKOtTjWmweg67jD8ngE4aNvgP+U+nh4NVz0yhktMumu5NdoNJrBx9wWhiCOXauCp/t+XENcAKQf9q5M5uyGDVpgNBrNmcf4ZSpmEiZF9YCnHZIBziC0wGg0mjMLI0Yy8hOp/y7vqNR/RwajBUaj0Zw5mFvyt2wFkeIwdOcB2H2P9fvPoHkzOsiv0WjODOrKYcf3rLEXCRRcBJ014DuRgi+VsP9BaPkICi6Ag6siRZxnQFGmFhiNRjN8MdxhPfVQ+wrRsZagsmTC6cmpIAiNa9UfgzOkKFMLjEajGZ5YUpHjkUpxicEZMm9Gx2A0Gs3wxJKKnEGIHJixfNhbL6AtGI1GM1zxFqT7DJyRXXDgYRhdoj4P447MWmA0Gs3QJV5bFl9LjJ0EKat7SZRgN3x8G/hPqoJMg2EW/NcuMo1GMzQxpxxv/Yo19beuXHVGdmU77Cghu2jQTjMmvkaruEAk+D9M0AKj0WiGJuYYi/nGbAhP41pAwvgrYdL1kZoXdy58YjWqDX+GITyZ69pLAu0i02g0Q5Pxy5RLKdBpzcoyC0+wB3JnwCf+jyp4rCs3uZ8SbMM/mEi/qpWBSEwGhmx3ZiFlmn2RGcKiRYvk5s2b030aGo2mLzjFYMzpye7cyHjjcMqyG7LGQ09d2k67d0JDzlzZgFRCaVxLhomMEOJDKeUip3Vpd5EJIT4nhNgjhNgvhLjbYf3lQogWIcRHoT8rTOsOCSG2h5ZvNi0vFEL8VQixL/Q6ZrCuR6PRDCJFVynrxHzTLbpK3YiLvxu5IVtSlgMZLi4Qtq6C3ZFZNUMwPpNWF5kQwg08ClwBHAU2CSHKpZQ7bZuuk1J+IcZh/klKae/xcDfwNynlgyHRuhv48UCeu0ajyWCKroqIzu574Oiz6T2fPhPKdLNbMEOsODPdMZiLgf1SygMAQogXgKsBu8D0lauBy0PvnwHeRQuMRnPmsfse2P/TdJ9FEkjABfnzYcJnh2yNTLoFZgpwxPT5KFDisN2lQoiPgePAD6WUVaHlEnhbCCGB1VLKx0PLJ0opawGklLVCiAlOXy6EuB24HWDatGn9vhiNRpNmjJiMt0A1mGz8a7rPqB8EofVDaN8OC19KXlzSOMI53QLjMEYuqgJqC1AspWwXQlwJvA7MCq0rlVIeDwnIX4UQu6WU7yX65SFBehxUkL/vp6/RaDKGunLYcoOKWwwngj1Qs9oqDomKhjnhIVYRZwoFKN1B/qPAVNPns1BWShgpZauUsj30fi3gFUKMC30+HnptAF5DudwA6oUQkwBCrw2pvAiNRpMmjPkqu+9RrfiHm7g4Ea/A1E6sWqFkjpUE6RaYTcAsIcQMIUQWcBNguUIhRJEQQoTeX4w65yYhRJ4QYmRoeR6wDNgR2q0c+Ebo/TeAN1J+JRqNZnAx3xz3/xS6atJ9RinCpWbJGNhFY+/K2EPMxi9TyQHgnCTQmwD1/8zTh5TSD9wFvAXsAl6UUlYJIe4UQtwZ2ux6YEcoBvMr4CapincmAutDyzcCb0op/xLa50HgCiHEPlSG2oODd1UajWZQyNRuyQNOEPb/PDIZc/wyawuc1g9jWyBOKdtmehOgfqILLUPoQkuNZohRVw5bvhypExn2uGHRq0okKj9vHWBmUPxdVRfUF/oZg8noQkuNRqNJHqc8oeFKIOLCKriAqF5qyVogTsWqA0S6s8g0Go0mORrfPjOC+gauLCUgdeWhfmUB1Rxz9MXga4WR50YEKEPqZbTAaDSaoYm52SUuLKOP3aMgezzkzQFPLtS+nK6zHEBC1po59iT9cGqjem0P5Thl0EwZ7SLTaDRDE3MAe/znrOsCrdBdq1xJbbvTc34DTbA7EisxAvOQ0TNltMBoNJqhixE/KL7DetMFdaOt/kXkyX7IIyKB+BnLI/Nt7GRQzzLtItNoNEMLp6wnw5qpWQ1N74Qyy9y2p/sMGJXcHyZdp16336WmdZqvLX+++g0yrGeZFhiNRjN0MLc+OfKkurFmT1QWjNFB2RgsNvJcqP+jKRFgCIvLiJmQNzty7a5sFfQ3uiyf+4B1Hs72uzJCaLTAaDSaoYNlWmW3KjIEaPobLHxRvT+4Sm3TsRtyiuF0dXrOdUBwgxDqGvb/HMucmPFXqmmdsYatZUCwXwuMRqNJD31p2GhsZ8kcM2EEwDsPWjOsTh9UsQp7IHzIEDAZXgHCbj53bsRqM+PU+iWNAqOD/BqNZvBJtMmifTtQT+WjLsJSZOnKVi36m96xHSAIo84H73iVuuzKS8XVDCIScKsgv5NwpLj1S1/RAqPRaAafRJssxtquYxfhm+2oi5R7rOWj6LYxwgMtW8DXqFKXgx2puJpBJqCC+U701ntskNECo9FoBp9En7SdtrM0uQzAmEvU2xM2kRoxE2SQIR3cd8KVpVyBsay+FLZ+6Ss6BqPRaAYf40m7txhMrO2MOIxZdOxxFvcILNX9Qx1Xtsqaa9+hGl02v5sRVko8tMBoNJr0YKQVxyNezYsxGrlmNbTvwdIuxpWltts/TIos8+erVOTGtyOZcxkQxO8NLTAajSYzcUq5BWtGmVO7fu84KL5d3ZA79g7hPmQhwbTXudittwxGC4xGo8lM7AH+Hd+DnjqVknz0aSi83HkWjO+EqoUZXQJTbhm6AuPOh/xZMGuFs/WWAYWUvaGD/BqNJjOxN3XsqolU5RvC48py3jfQqVxnGdL0MSkCrdDyIZyqtC7PoCB+b2iB0Wg0mYnxtJ4/P3qdO1d1Sh77GZUt5jR4rPEtaN0auynkUOHoH2KvM9rCxMooSzND/JfXaDTDGuMpPdyDK0uJSsEFkZYwuHFORQ7AyX8M4smmiK4jakSyuXK/rjzU2PNvEZdhBmaUpd2CEUJ8TgixRwixXwhxt8P6y4UQLUKIj0J/VoSWTxVC/F0IsUsIUSWE+DfTPvcJIY6Z9rlyMK9Jo9EMIObiwYUvQcmbqtDQXAtjHx88rAiqtGSj44GR/NC41uoyjOUOTKOVk1YLRgjhBh4FrgCOApuEEOVSyp22TddJKb9gW+YHfiCl3CKEGAl8KIT4q2nfR6SUD6X0AjQazeBgpDQbN0tvgXKTGdlUM5Yrl1jLFpytGRdM+hLUvka4YeRQwywi9l5ssTLK0tz8Mt0usouB/VLKAwBCiBeAqwG7wEQhpawFakPv24QQu4Apieyr0WiGIOabpSEqxvwTY33Mqv1QjMZbqLLMhlR1v6nBpXGtRqqy4TJ0anwJaW9+mW4X2RTgiOnz0dAyO5cKIT4WQvxZCDHPvlIIMR24EDCnW9wlhNgmhHhKCDHG6cuFELcLITYLITY3NjYmfREajSYB+uuqsd8sDXGpWQ1bbuhlcmVApSv7GskccUnw9uvKUq35DevDcBmOvzK+uEDam1+mW2AcUj+i/va3AMVSyvOB/wJetxxAiHzgFWC5lLI1tPi/gZnABSgr52GnL5dSPi6lXCSlXDR+/Pjkr0Kj0cQn0e7J8fAWRDLC3Lnqsz0WMaRIsI1NsBvaqpSQmn+35netsRkn0tz8Mt0usqPAVNPns4Dj5g1MooGUcq0Q4tdCiHFSyhNCCC9KXP4gpXzVtF298V4I8QTwp1RdgEajSYD+umrqylXWmPQTbldvCfQPc7pq1J/Gv8D4z6llif6eibTkSRHptmA2AbOEEDOEEFnATYBFioUQRUIIEXp/Meqcm0LLfgPsklL+f7Z9Jpk+XgsMk4ZEGs0Qpb+uGnsH5bpyZcFEZY85OUWGKMKjJnJaCGWUNf0tUmRq/J4ZWBOTVgtGSukXQtwFvIX6l/KUlLJKCHFnaP1jwPXAt4UQfuA0cJOUUgohlgC3ANuFEB+FDvkfUsq1wC+EEBeg3G2HgDsG9cI0Go2V/rY4sU+ybN8BBw/ApGuh7nVl2bhzYcKVQzhTzNSsE0DkwJjFUHuEKHeafWQyZNSoZAMhZaYEvNLLokWL5ObNm9N9GhrNmY25ezJYBamuHHbfYw3mF3830q7fCPg3rk3PufcH73jVd8yxMDSURWbGnWsVke13qfiWQfF3VTuZQUAI8aGUcpHTunTHYDQajUZhTkM+8gQg1JP64dUw80eqozBYU5UN8TEEqK3Kesy4N+504sZiZflOxTlHk7iMukgNWLNbgWF3YSCjuixrgdFoNJmBOc5i7pIs/bD/56o7cixXm1mczPhPQksmeibsniNf7E2FJ+ICnL0i2vVlJEAQUNvOWJ4R7jHQAqPRaDIFc5zFlQVBP5HYQyCSKeWUFWVJAjAh/ZlT9mKhlxTlrCLIGgcjzwV/6LoSKaaUfmj5KHqbNJHuLDLNEKK5ooID999Pc0VFuk9FMxwx12yc/UMQptuTKyu+28dbkPrzGyyEBxasVi7BhrWR8cixGL9MjVM2aHonYzLJtMBoEqK5ooJ9P/gB9c89x74f/ECLjCY1GLNOfC2hmpcQYz+j1u2+B979hHo142sZ3PMcEFw43oINYa1ZHV3r4kTRVTD205HPwZ6MmYOjBUaTEKc2bCDY1QVAsKuLUxs2pPmMNMMae91M8R1KVPb/VGWR7f8pfPjlSN2HfThZpuMeBaMuVAF7O8Ee2LsSTphEwpUd34IrviOtLWFioQVGkxCjS0tx5eQA4MrJYXRpaZrPSDMkSbQY0KnFiX2f2lcirWcg0p8rCkHGFWAGWqH1w9iZY61brBZc/vz4gfs0t4SJha6DCaHrYHqnuaKCUxs2MLq0lMKysnSfjmaoYe+G3NcboWHBOJE/P5LGvPPfobO6/+ebSbjzYcb3ItdorhdKs5joOhjNgFBYVqaFRZM8sfqRJXqzNN9cR56rAuDmyv7N14Y2TLCJZCZjNPU0rJhAe0RcR5dkZNW+E9pFptFoBgenfmR97bJ87gNw+Xa46CV1Y82fb1oZZFiICy5VWHrRK+AptK469Fjiwf8M6E2mBUaj0QwOTnECJ6vGzIdfhr+MVa9Oxzv3gcjT/rAhqLLiiq6C6XdaV/mbnRtdglVQBmI8wgAw3P5mNBpNJmMvkjQXV9qznzaURoLgtS8rkbnopejjzfwR7H+Q4WG9YP0dDLfgoceUuEB0o0vDzWh2mxVentZJlgZaYDQaTfqI1/rFnmHVaKu9MsduFr2mUntbP8Kpk/K69kIqO0dTknuKpfnNqbmWgSCnGOb/yioG5z5gjbsYadvmbeyWIKjtnIR7ENECM4g4ZWHpzCzNGU+s1i92xpv+f9if2C98Hi7bHBEdb0HYqlnXXsiKujl0STdvtk5kZdGezBWZkfOcEx96G3dgtwSL71B/0pxppgVmkDAq4YNdXdSvWcOU224jf8GC8LLGV19l1sMPxxUZLUaaMwZjPLKRRTXmkxH3mNG23/zEXrPa2qesrhxypkJXDZWdo+mSajBZl3RT2Tk6/QLjHgWBNqIapRVc4CyevQlELAFKc3aZFphBwlwJTyDAsSefZHRpqaU6vn7NGkcBaa6o4Oivf03Hzp0gJfUvvEDe3Lmc9Z3vaKHRDD/s45HP+bE1Rdmpa7LRf8sejwBKck/xZutEuqSbHBGgJPdU385HeEHG6XacDCOmWefaGDS8ZR0FbU586E100jgaORZaYAYJz8iR1gWBAK1btlgWnXrvPQDqn3sO3G5yiotxjxhBR5VtxkUwSEdVFXv+7d8Y/clPMvHGG6OExmztAI7vtThpMgrDLdR50Doe2dxnLFbXZKP/lj0zDVia38zKoj3Jx2AGWlwAvKOImmAJqrrfnW3btiB2DVGGoyv5Q6Sykr+5ooLqe+/F35was9yVkxN2rx1etYrGP/6Rnro6CAbB5QIp1R8hwO0Gv9+yj0aTdsxWhysbkEo07BX/deWw+UtEBfKN7SAU7N9KZmeVOUypDK/KBtkd+WxM7exPF4QUoiv508jhVas49uSTEEjdjPBgVxd7//3fAZDd3baVpv9kUoLfH96nfs0aQFs0mgzAMmzMIQ3XoOgqGHWBetI3M2O5et1yg9o/44nzYG8WF/PUznhB/gxFWzAhBtqCObxqFY3l5cqSyNTf2OVCeDzInh5t0WjSS1/6lDlZMcXfVa/mufRDDi+WyZZGf7UMF5N4FkzaK/mFEJ8TQuwRQuwXQtztsP5yIUSLEOKj0J8Vve0rhCgUQvxVCLEv9DpmsK4HQlbL6tX01NZmrrgABIPIHjWaNtjVRV3IounLYDE9hEwzICTaDdiI04wpCS9a117IQzv38tiuXTzUcA7r2gsddkyym7JndHL79UbeuQ4LA5HBYe7cISEuvZFWC0YI4Qb2AlcAR4FNwFeklDtN21wO/FBK+YVE9xVC/AJollI+GBKeMVLKH8c7l4G0YD789KfpOX58QI41mAiPh8nf+ha1zzyjstvcbqbcdhvTli933N6ceq0toDOEdHbxtWSHuYGApcZFuZ0EOSIYqnVpAlxQcCGM/yzs/xlm11RCxZfmVOlYuLKVtdH2ce/bGnjHg68xerlThb5TEWoGucoyOQZzMbBfSnkAQAjxAnA1sDPuXr3vezVweWi7Z4B3gbgCM1A0V1TQU19vXehyWWMhGYr0+6l78cWodOr8BQschcNpCJkWmGFMMvUZA4klOywAwmOpcTGslC7pojL/ZpYWB6034Y69quUMJF582ZtgiBzInwezV8CxZ6H2dSABkXESF3uFvtPvDUOmkzKk30U2BThi+nw0tMzOpUKIj4UQfxZCzEtg34lSylqA0OsEpy8XQtwuhNgshNjc2OjwF54EpzZssAb0hRgS4mIQOHnStiAQc3qlHkJ2htFbY8pUY+/GXHQNJXlt5Ajj/5uyTnJcUDL3BjV6GRw7CjsVXyaF7FIJB5uvDolXghaMnVEXRYuFvWvy7nsS76ScIaTbgnFyjNp9dluAYilluxDiSuB1YFaC+8ZFSvk48DgoF1lf9o3F6NJSGl5+WcU2hMjsGEyCGDU89k4ChWVlzHr44fCy9m3bOLxqFYVlZTHdapohTLzGlANFPPePPZOq8W2W5jWysihAZedo8l1+2oNZlEyawdK8ZqsFcORJS3ZZwsWXnsJIk8lUMuYS6/XWlaviUTPtO6Bzr3LJBbszajRyLNItMEeBqabPZwGW4IWUstX0fq0Q4tdCiHG97FsvhJgkpawVQkwCGlJy9jGQhgUzDMQFwN/W5tjqZtry5WGhMRIbAI7t28fpQ4fwjh2r05+HE6lOlY3lEqpR/67C7iPz9x59mqX5zSwd1aHiIO3bIbAftm6wdhS2pS4nVHyZOxOC/gEWGKP+xXAeBVXrfbtQNL6t6oDsBHtip3BnIOkWmE3ALCHEDOAYcBNws3kDIUQRUC+llEKIi1F/M03AqTj7lgPfAB4Mvb4xCNcCoDKxUljzkg5Gl5ZSv2ZNzNhMc0UFdS+8YNmn+a23ABLqsaYZQqSyHYndBVezWs0+McSh6R04+4eqst+pNqTx7Uh9jL2jsCs7dMOOPPQtzW+OU9UvVHJA3esDfJHS9Bp6Lx1c6GZr0Vx4igu666O7KWcoaRUYKaVfCHEX8BYqLeQpKWWVEOLO0PrHgOuBbwsh/MBp4CapUt8c9w0d+kHgRSHEt4DDgMO0otSQZDJkxuIZN472bdui4zCm2Ixh2ThhpD/rYk5Nr9hdcGC1PII9UP0LFXg3B7iNG+2pSuvxuutVAWbLR+pzxx7orE7wZCTUvoZT6/+BweTdkH4lprG6J3sLoPEtaAl1J2j9ELZ8GRa+lPEiowstQwxUmnJzRQV7/u3fwhXzwxaXi4IlSxBEeqjFxOPR7WlsrKteT2XNJkqKF7N05pJ0n86gE/P6zTEYsFXm23p3FX83EsgHqPw8NK61fpH56d+VpayFRFOJEyJOy5c+HSaUDm0vMt19T0RU7divP00MaKGlEGKpEGKO6fN3hBBVQog2IcQOw/o4UyksK2POL39JwWWXkT11au87DFWCQVree693cQFLexqjmNNguBdqrqtez0MVj7Cuer1l2Yq19/HKx6+yYu19lnVDgn7Oeo97/UVXqZum8RS/8EUVcxh/JZxztzWLzB636LaVB4ASJyOWEexRA70w0ppdhG+BSY1ddpG4z8K+ndv60RAQwzUI6vfd/3NncXGK22Qgyfyq/w38K7BHCPFdlDvqV8AuYA7woBDCLaUcyj0b+oXxhL7vBz9I85lkAC6XyqYLxaVaP/iA5oqKcOymL/NwhhrGjbTL382bVW+y8sr7WDpzCZU1m+jyq6fyLn83lTWbho4VMwC1MH26foeYz7o9a6nsnk1JRyFLjXNKNF339AEiFkdQCRdA+x44naj7zKAv5Qeh7xQeGLcM/KeiJ3YaNP1NWS515VhddG7Vhy174pCJwSRTBzMTOBB6fxtwl5TyHinl76WU9wL/CyVAZzSW+S/DCFdeHjlnn53Yxm43hVdcQdbEieFFsqcnHLtxKtQcTjjdSAFKiheT41EtQXI82ZQUL07bOfaZZGphbBZP0tdfV866bb9hxcERvHL4OCv+dA/ryi9WfclqHlVpvCFLZF3HeB5qOt+hbYyq9geUFVBwgUoe6LO4JIn0qwwwX2vsbYLdyi1mmRfjhknXqnTmISIukJzAtAHjQu8nAx/Z1m8BpvXnpIYDo0tLEV5vuk9jwAl2dtJ14ECv27kLCphy222c+p//sbTNMRdkDvdCzfzsPNxCuULMN9KlM5ew8sr7uO78L4WtmiGDvdixNzeNYfHUPKpe68qTv/6a1VS250UKJINBKk80EX7KD3bDuGWsy/0aK+rn8UrzSFbUzYkhMgBCDfhySgeOQmB1cwmst88+uNu8BTEEwuS6M7vF8ueroWsNay2/41AgGRfZWuAu4Fbg78ANwMem9TeieoSd0RSWlVFw6aWJxSiGEgkmhRTddBP+tjaLFZc1aRK5s2aFP9sLNQfKPZYJAfR11etZs+VFAjKAW7i4ceENADxU8Uj4vIaUsBj0tRYmxqCsZK8/boGk8EDxHVTurKYrcAgIVen7Z7OUD6IPFuyGnhMJfKsLRJaq2jfImQpdh62fz/qqSqH2FqiJnE6D0QAOPAxjPw1ZRdBTF1luuL9OvB1JXXZlq6aXQ3TgWDICczewQQixDqgEvi+EuIxIDOYS4JqBO8WhS97cuZxat27YFFwmSu68eUxbvpzmiopIVwOPB19TE6dqa2nduDEcbzH+DBSx4h79OV4iYrWuej2vb3sDCVy74GqLeywgg1QequT3m/5AQAYH5LzS2vCwL7UwA1n9X3wHS5veYSV7qOwspCSvRVXsA+CCmT+Coqso6VjPm1Vv0uXvVpbj5OnQ6SAwAGMWQ92x2JllntEqXmIWF3cunPU1q4h01cCBh1TqMKgiz+56ta89NTrYHcp2M6yeULHl7BXq79R8Llkh9/JgdFFIAX0WGCllnRBiIap55NUoW/FiVFX9BqBUSpma0ZBDiOaKCmqfeeaMExfh9TL1O9+JLDCuPxBAmrLJUtUYcyAD6LHEyi4mAPe8uQJfQM3y2Hx4Mzcv+go5nmy6/N14XR72Nu4nGHoq7fJ389q2N8Lna7jOnvrgaU6ebuFzc5dxZ+ntsU/MFGhft+tVKnPfoWTuDZlpEQ1k9X/RVbDwJZY2vs1Sb4Gqb+mujwp6Gy648INBXjNsfTUkBra04qyJcNErKnOrZYvVogAlEGZENky4Um+xgqAAACAASURBVFkqE64MN88ElKtt6y1KjCxuN9X5OZog4IKcaUqwjN/GEBJQVpJR8zIEB44lVWgppWwB/iP0R+PAcA3y94qI+KlPbdiA9IUGKEmpxjUHAimNt5QULw4/vXpdHmpbjrOuen2v1kdlzSbys/No7+4IWyuxgvRmMfngUCU53pzwZwBf0E/loUoWnnVhOKT8j0PWJ+iNhzex8fAmAsEAr3z8KgKBDN34ntn4LEBskQm5S1RH4Ol0ya28eXBn5sZziq5iXUchlTs3UdIR/+8ikWMB1uFks1dE3XCjXHAzlithHnmuimUEOpXV0HlQrS95UyUi9DawTHZbRcVOwCl4H69YM6hE5OAqGF0SEeQd/xpxwQV7lACWvDlkhMWgV4ERQuwAtqKC91uBrSGB0RDdANJgdGkpja++Gp6VMvpTnwq3TxnOGFlihWVlUb/BpG98A39bW7/jLbF+c4g8vb627Q02H/mQfxz6gC1Ht3LjwhvC4lFVt5P11RtYMrOUeUXnha0Ug/Lt5SyZuYQtR7dav7ezmde3vWERE4nktO901DnublBhSIFg/qR5uIWbgIzcaAJB601H2or11ldviBKYsBD6RtDeOJvaHnck4J3B6c4D7bbsczyirjzizjp9IFLd3/Q35apqfEvFP/JmYC3mNAL5qajmt1lS5oy8xrche7w1xjNEScSCmQqch6lHmBDiIFbR2SKlHNSGkplAvDoOewAboPmddyJ9ylwu8s8/n/bt24dd1b/RfTkVQfxEameWzlxiEYIuf3c4/vH69jfCN/fqpgNMKZhsERdQFsjf970b9d1Oy3pDItleu6P3DW1MK7QmYppv0ooJeAW4hYuADOJ1ezM23bkvbsuEYl6JxCPqyiMFi2AVJF+LShVuNH7LgGq/YvQxCyPVOtyQNws6did6yYALJn0J2nbb0o1R5xwWuXeUheLOVckBhmXmyopU97uylQvQfn1DwF2WiMAUAS8D/4yqf2lCzV25DviSsZEQopZo0Rn6EhyH3gZumQPYu+64w9oEMxiko6pqeIiLbSxB7TPPhBthxgrix7NC4mH/zTe8/Dhrjj5DblYunT2dYauk8vCm8D6u0E0Yoi2HYy2ZOXl0/YENXP/UjfiDAc4rmsvh5sPRQihBhBsmJh/rS3XWndltGbPupa6cdTtfYsWOo3QF/PEtnd7iOnXl1hYzwqNu2MaN3BCkw6sTbBsT6KO4AARVfOfcW6yiMfYz0UPFzM06w92fe9SMmGB3dFJFuge/9YFEBOYXwKeAT0kp1xkLhRDzgRWoZpQ7UUO9Ph/6A8rOTHe35pQS5QaLEVdorqig5f33o5bLnkTy7zOfwmXLOH3gAKf37QN6D+L3p4J/dGkpda+8jOjuwecRvOE9QHVTdnh9ddMB8rJyLUKS5c6iyz+04mG+gC8sfvVtDi1QQhiuNV/Qn5SLbMDdVw5EBd3txw/dMCtri+gKTAaSq/AP0/i2tUmm9MM4h1HEo85Xgf2oXmKhivn2HVFt/uPjAuGKWB2GkOXNVenQZ31VpRzHu45wgN8VGcF88EAkPmNc3xBJWU6k0PJ64AWzuABIKXdIKW8A7gOygHNQBZZXA/8P8MeBPdXMw3ABTbz55rg3SUuwGyyB8KFA4Wc/y+jLLou5bs6qVUxbvtxSNOkZOTJmj7H+VPBXFWfxzD/lsW5eDr8ry6dqenbUNh091vqDoSMuyVshLuGyWAZOPdCciJXIMBCYz2HpzCX8sOz7zoIRumGW5J4KT6fsV4eD8ctCTS5DuLKU1WD0ODMsgJYPAZeqR8E0dtk7BiZ8VvVByym2HVzApOtVo8lJ10eKKoVHucSEcUuVqrvzli8r11tXjap/iVcgWXSVcp0JD2BqymnvmNDXYtc0koiFMRqIWY0kpVwphLgG+A8p5d2oQWDDXlwMEqnjGF1aSsNLL0VExu3GM3Ikfvt44kzCo/5peEaPpmPnTnxNTVGbTLnjjvDkSnO8xTNyJLXPPGOxUICwSyxRy8+M4cY53nKcbdPcbJuWP4AXmwmYWpj0gtvljnL1BWWQqrqdLJ25hMc2PB6OOZVvL2fxtEXMmjDLkiFnkJD7Kgn6ZBmFYipL85tZOfkQlbnX9y/t2miSaR9UZrB3pakIMgA9japSvvEtJTq+E7D/p0pA7IxaCBe9FPlsd3GZG2se/b01XTnY3bu14WuJdtvZRSTVg98GkEQEZg/wT71s8w5q5srd/T6jYUhUVb/fTzDT3WOh2JD/xInoKeNud3iipRlDbA/cf7/FQqlfs4bWjRstgtOX4L/5ZuV2ueNuO3RJTFzmTJjDvsZ9juvWV29gXtF5/H7Tc+GYky/o5x+HPginSdvreSprNnHjwhvY27BvQGcZ9bmhZeiGuXT8MpYmesPsbbwyKJExhKboKtVEMiqYH1A39oDNHVb7Char0pWtUqLt32P+7iNPKFERHui2xffMbrNY2AeNjf20c++xVA5+G0ASEZingFVCiO9LKR+JsU0hqi/ZGU+s4PXEG28M32Rxuwl2dKTxLJPHXVDAOT/9abgbciIp2hKiXGJn33svhWVl6iYXap8CWPz0j214nPXVG/C6PZGq+GAqUkaHDkdPHQkXbNpZMrOUyppNlnRoO2Y3mCHaAoHL5SIQDLDl6NYBicP02TLq6w2zt0B3XTl8eF3EGmh6RxUrHv2D8/Fat6r995szvkzikj9fxU96PceQTMsglm7LOcUw/1e97z+ErJNE6HXgmBDCjeo/9hnU9MiVUsoPTOs/FVpfJ6WcmcJzTSkDMXDMHLx2espvrqigbs0aWtavh2BfWn1nDoZbzHytToPEzOIDOG5rduW4hUsJkQyS48nm0hmXWtKCzYWI8Zh3qJvZx3zsneJ1jM8MR0Z4R3DDhdczr+i8cP2PL+DDJVwIISyi7HV5WDxtERJ4/5Bz+5Trzv8SPyz7fr/PK6XZafaiyESGj426KBI4d+KcUN24vSBTeFQbGnuA3kxdOez4noq1GJiHiM1Ybh31PIyIN3CsVwtGShkQQnwBeBTVnv+zQogm1CjisajAvgD+a+BOeWhiqd53mFt/asMG9XwzRMXFM24c+QsWAH1L0QaiXGLrqtdbXDkB01N5l7+bd/f9j+W7ExWXWyrayPJDyZ4uni3jjBAZQ1zu+dO9+IJ+3MLFJ6dfwjWhNjZGl4K9DfvCxadetxcXgqDtdx3IOExKG3om05urNTRyOBZ15XD59oiQ7L4nMvDLXGnvtJ9l8ibKvTXxi3BykxKaAw+F4jI2aytV9SwZUieTUBqxlNIH3C6EeALVSXkZsBBVhbQH+IWU8repOsmhwujSUurXrInUuzjMrRdeL8LjCfflGkr4T5xg7/e/z+xHHrEkLgivN+HWLzvrdvJBxcdU1e6I68pJRFDszD7mIyv0s2b51eczQWDWV2/gg0OV+ILq4gNSyYZxczdeH6p4BF/IavEFfEwpmMzxllokEq/by6KpF3HtgqtTIgpJWTO9xVicXEnGPgUXhLoSm/+fBYk74tiVrUTFsDR8LYSr+OOlA9vTogE8Y6D+j9HLzcdJRT2LUWDa9Df13Wmuk+lTnYqUchPwDQAhRDYQkHJAB1wPaQrLyphy220ce/JJS88t89O+JV05hPB4ENnZ6YvLeL3gcF5OGK1gRpeWYrhXndyshsXmGTmSjl27aHn/faTPR8Aj2F2Wz+4U3Pj3TvFSsqeLLD/0eNTnM4Hqpuj5PCc7my2jAcDWp83tpaGtAYnELVzcfNFN8Rts9oOkam0SufnGK0DEDbnTbZ2M3TCmxDpJMu9cCJyG7mPWav6jTyu3lju3dyvJW0CUcPXU4yhkRv8zQwjj1bP01QqxXD+xjzuIJF0IKaXsSwVSTIQQnwN+iUpEf1JK+WCM7RYDHwA3SilfFkLMAcwD3s8GVkgpVwkh7kNN1mwMrfsPKaXNIZsapi1fTv6CBVHBbyPobWfErFlMW76c9m3bOLZ6ddT6QcHnU2nJIavKO3Eivnrn4j6RlaWKHdesiXQh8PupW7OG9m3bOPqXP9Hq9jGq5gQiEO2OyPLLmJZFf+MnVdOzebaMMy4GY2DEqdwuN9VNB9ndsNdyQzcXPB5vOR6OwQRkkPbu1D3cJNXhOpliQvM+BNSN3JWtnuSNOIqvxSowuWfD6cPQZXtONlrKXPi8teWMHaPPWZSYSCx9zYQHRp6vijcb10Lzu/EFLBnrxnL9IdJcJ5PWSvtQAsGjwBWo+plNQohyKeVOh+1+jkoyAEBKuQe4wLT+GPCaabdHpJQPpfYKnLHHH4wakcOrVoWr3QFwu5m2fDmFZWXpHxfs94PLBcFgTHEBmHzrrRSWlSlXoImugwc59d57CGAUsZNuY1kWAxU/qZqefcYJi8HsCbOZP2keO2qr2NOwB1A39Ne3vWFxTxkpyluPbh3w+hcnkqq1SSbGMn4Z1DxGpDllUKX52iv4w2nAWRFXkh3zdza/q7Zvflfd6E9VquMUXRVKbzbf1EOWjLnfGKhU48a3IxaSWcCcrJRkBNaS4uzQliYNpLuVy8XAfinlAQAhxAuoTgA7bdv9K/AKEOtf5qeBaillTYz1accQHHuWmbHcqZBx0Ekg+cDf1sZjGx7n6IRa/tnjxuUPILKy6Go9GRYVs7gY5YM+lxKWD+bmcGzuBDhtbch9psZPBgqv28u3LvkmoLpBG7hdbjYd3owv6He0Zoy5Nqmk11YxTiSTrlt0lSqYrP5FJHvLfoM1H7fzoDXTLH9+RDSM79x+l/VGv3dlRCT274DcmdY+Z71li9lFM1Z6djICm4EpzukWmCnAEdPno0CJeQMhxBTgWqCM2AJzE/C8bdldQoivA5uBH0gpo8rmhRC3A7cDTJs2zb56wHHqLmzEKtq3b0/59/eKrWmlHQlsa9nPMxvfhkKoK8vlat/ZnD/7Eo488XhYWAxRCQjYPj2LtlyXxWUlTrdGpR2fqfGTgUAgWHK2qoGpqt0RDvYDjMzO51RIzJ3cU1tCVszWAap/iUVSGWXJFBOe+4DK9op3kzWOu/ueSCKAO9e5zmX8MjjypLJyXNnRI5Y7q0MFlFc6Wwv2OEqiApCsWGRYAWa6BcbJk2K/w60CfhxKl44+gBBZwFXA/zYt/m/g/tCx7gceBv4l6oukfBx4HFQdTBLn32fM7rPmigr2Ll+uAv+uRNrCpRbPmDFkTZrE6T17LFluhmAIYOxfNjHv0/lhV1TX2Gzm1bThCkZ+vhMjXeyelhUzDuKUIXamx0/6g0Ty3v73LKneAC4ErV1t4c9el8finhro6Z/xLJRUd2wGrDfzT/wf9Xn7XbE7Lh9cFcoycyvLI+aNWUZexyyGWpujJNit3HD2AH2sbK7eRMggw8QiGdItMEdR82YMzgLs/dMXAS+ExGUccKUQwi+lfD20/p9RowHCgQPz+1Bq9Z9ScO79pm7NmkhWWTBI1tSp4PfTU1cX25IwBeMHGn9zM8HOTiZ/61t07NrFqX/8A/x+y1OAJ2AN0k8rnMa7zSc4x+vC6wvS44HyS/N0/GQQceGKEhdA1biY/h2dPW5mzD5kbuEmPzsvqe/vLUtsMDo2RwXFZyyPDBk78kR0PMKeEOBriRzHfLO39xfLmqgKMo/+AbqOAMHoFjB9yeYaQq33kyHdj82bgFlCiBkhS+QmwNJuVEo5Q0o5XUo5HTWX5jsmcQH4Cjb3mBBikunjtUDfJz4NAnZ7LHfGDC6qqGDK7beTNXky2VOn4h4zxrpRiutngl1ddOzaFXMQmvR6yF68kFHZo7iytYjC377JPw6+z7OhDsfPlo2kano2uVm5jB5RwLxD3Vy7oZ15hwYk6VBjYoR3BHMmzOaWi7+K1927S7Ew1/pvaenMJdy48IbQ0LIAa7a82Gv3ZSd668icyo7NYexB8bpy62yVxrXqRm50MzZ3XDYEwrjZ1zwa2dapc/G5D6i2Ly7j+dz2MNiXbC6nYP4wIq0CE6qhuQuVHbYLeFFKWSWEuFMIcWdv+wshclEZaK/aVv1CCLFdCLEN1aiz/30vUkDe3LkW15hrxAiaKyqofeYZeo4fp6euLi0zY06tX2/p9CyBznEjCS6aT/DuO3g9t4apexq57OUdlFad5paKNiTwWmmkff7pntNM2dXALRVtLK3q4paKNi0yA8xp32kONB2k8lAlgUD8Bw+vy8O1ocp+M+3dHWHrJ9mbf0nxYnI86u/dKUust/UxMVxcphb3ljEE5vV2ISi6KvLZIOoGLq2vsTK3LnxetaIxWxd2yyZWO31XlorPxLJMhlDr/WRIt4uMUH3KWtuyx2Js+03b505Uuxr7drcM4CmmBENIzJlbzW+9RffRo5aiTKfCzFQisrOR3VYhEMCMZV/k7Hvv5aGKR+iq7+4160sidWbYIOAL+NjdsDfuNjPHns0dpf/L0S01EO36e8sSSyqLzMF1tK6jMOJq21HOyqI9LM2tjbiW7EHx0SXWOIj5Bu4kELEyt5xiIfGyvPoTzIfYcaMhSNoF5kzF0rfMhO/kSVw5OY7r+kv21KnIQICe47HHBIcTDswpyx5PuBWMcUPaO6U7nPXlc0Fha4B5h7otAqIzw9KPQMQUF0jy5h/jOPH27XMWmYM1Udk4M+JqC/ipbB8BwUIqO0dTkv0SS8uedQ6SOwXRnQRiILO8zJaO+XOsY6WqdUya6bWb8pnCQHRT7guWzssmptxxBwDHnnhiwJtiuvLycI0Ygf9EzPlxAGRNnUqwvZ1gTw8jpk/nrO98x1I4amQETd5VR/Zf32fC/hN4ApIeD+EYjMGZ2N04k5gzYTa//epv0n0afaeuXE2DDPYoN9PClywWTI7bw40FR1hzcgJd0k2O28PKz9/fNxFLphVLotvvvsdaj5OIWPTWITpD6Vc3ZU3yxJqXAtaaGF9TE6cPHKCwrIz8BQs4vGpVlAUxEMH9YEdHQv3Oeo6o0iRXTk6UuBgcbzlOeeBDvuBpY3JAPaQ4ucF0ZtjA4dT92MwI7whO+05bll0yvSRlKcKpTz22lu7ara3KXS/S1bwVCFk08dKsncShL2nATtYFOAtOXbnqwmxulFmzum+V+MMkHqMFJkWYLRRjiqOTyNhnqDhZNanOHItFsKuLXWtfptRmvRhPkaDdYINJPHEBosQFYG/DPp7b/Dy+oJ/y7eU88IX7AfotDClPPTZ3KDaNGra42k5u5M3qj+gKyvjxo0RcT71ZJ3aXXc3qSAsZ+zEb3ybSriZE0zuR9jJO52d8d4ZV4veXdKcpD1uc5qX0ZR8AsrJSdXoJERDw++DHltTV17e9ERYXg32TvFRN9Ua5x5wQCApzC1NyvhorOZ5sTnY2hyv7fUE/v/ngaVasvY9XPn6VFWvvSyotGZJMPXbICotJb9lVdeUsbbyXlROruG50AysvvTK2wPWWCuyUnmzHW2D93LEn9jHN525gzzSL9d2g3GLDQFxAC0zKGF1aiisnByDctr8v+wAQJ0XZlZ8f+0AuF55x46IW25IyVWsY234GQVSbl4+musI3j3XV69l0OBKnml/jU+nHR3zMqk0s221E1giaO5sT2lYTG7dwMz5vHF6Xl+Ix0yxpwN+4+BauO/9LrLzyvigxbzndMiA1KX1OPU7kJm4mVnqwQUg0luY388Pxe1maHT2yIEwssTIEr2Z177UoPmvvPDoPRd67sqzHbHxbFXqOv1KtM3+vXWSHeR2MdpGlCKe+Y4nuE9V1GSyZXcLjYdJXv8rxp592rJPxFBQQaLH+hwieU8y7o5qY2NhNwWnJ5OnnMvWiUmqfeSbcfLPwM5/hxP7dyOoa3MB5R3pYcNhPyVXq5lFZs8nS52ruMX9UGjLA3OMBdk12O1oznT2dUcs0zriEC7dwWX5zUIH7b11ya1S1fCy316bQCGWv28tn5y5jzZYX+91Fuc/ZZ8l0B44XI+lLvMIp48vsNnNlWxtWOh3L/H24sbjAxn5GvVZ+XrnCjONc+Hyki7JxTLurbhjGXczoLLIQg51FZsc+w37P974XmYwJjL7ssnBPsIk33ghA/Zo1tO/Ygb/ZahFkTZ0aDtQD5M2bx76xUNV6gE/tOE2WHwJZHs575JdqDo1pQNqJmeMprIrsu/+SYm55+i9AdPzF3GLfaGx53pGecDzG7DLLcmfRExj8otGhilG7AvDatjeorNlIUAZxCxc/++IDvd7QzYIDWN4bHZRTNb3SEfMNPdGsqkSOmWy8wp6xNf5K1U/MWxC7G7Lxfd6CSBsao4Oy8dmMPQssVpZYhow3TpZ4WWRaYEKkQ2DMUx8NS8KVk8Oshx+OuvGbkwTMyQAiKwukDBdk5l94IR1VVRHLxuVCuN1qmqQAt+mvu+Cyy8g56yzqn3suvCx4TjH+Q4fJ8kt8LnAvnM/cW78d/u511etZveEJqpsOMO9QN1ds6WTKiYB6prMdf/28HF4tzcclXAQdemVprKiWLUpESoov5pqQADxU8QivfBxpVnHp9EuYXDA5bDnYrRfzg0COJzscgI+1fNBI0420tbWVhoYGfOai5UBnqDNy6LEtK+RSti+zx1LM+we6wJ2jXgNttg0c9nf6zljHzwC8Xi8TJkxg1KhRcbfTApMAAy0w8VKUjfXm2TB2ayX7rLPwjByJv60t6hgH7r/fIgrm7ZsrKqLdayYMKwgAl4vcuXPp3LMnnKkmsrLwX/MZTlVtp3BvLcLnjxK4xzY8zuYXnwxbL2YMkenxwB8+U8D2aTqrLBE8Lg+BoN+SJ2aIABAWBq/LA0LgC/jI8WRz48IbLC4vw23lJEjmKZYAn5x+CQ9f+/8O0hWmh9bWVurr65kyZQojRozA0pHddwp8reAdBd7R0HkYehoi67MmQG6cMR7G/tIHPtM0EHcu5ExWx4y1j/GdGYqUktOnT3Ps2DEmTpwYV2R0Hcwgk0iKsiVjLBAIi4zIygrPr7ff2A3sw8lcI0YwurTUOcXZhiWsHwzSWVVlCfbLnh4me0Yz+fyl1FcpETOy4IzzaO/usLSBMejxwP/MH0GOT6rCSi0uCeMPRqeiG9MoJxVM5saFN9De3UFty3H+ERKJLn8366s3RAXtze1fvC4Pm498iO/QB3jdXtwuN4GgepjZdHgz66rXJ2fFDBG3TkNDA1OmTCE318FS8I623uS9o8B3AmQQhEt9joXvFHQeUNvaEd7Y4mH/TvsxM0R8hBDk5uYyZcoUjh8/3qsVEwudRZYCEklRtmeZTbntNibefDOjLrkk7O6Kte/pA9aMmfbt22O2nklozoyUSLcrfC6jS0vjZsGVFC/m4LRcekKPJwG3i+Z5Uzl251W8c+loXivNRwiXYxflbHe27rDcB94/9AGvfPwqz21+ntqW48yaMMuSvbVkZmlUNpcRgL/u/C+xeNoifAH178kX8DEuL5Jd6Av6k8si62tGWD+wNLfsC6FsLV9XCyNGjEhsH+9oyD07ZLmcHf8m72t1FhcAfyu071OCkSiGYPU0qNe+7JtCRowYYXUt9hFtwaSA0aWlNL76ajim4pSiHCvLrLmigraNGx33NdxuI84+2+IG66mvxzNyZLiHmfB6kVIqt5chMPa2M0IghUAE1QyXdQtGMC1rHGdf8YXwucTKgls6cwl89+fsm/YSc476OO/K6yksK1M3gT9tZN6hbr4Wcp+V7Oni2TJV0e8SLuYc7ORmh3XDCY/L42iRJIPhMvMF/fzj0AdU1mzisnOWUphbGBaTeUXnRWVzGQWJ66rXh6dW5niy+dxAZJElkxGWBEkXc5oTCiZ/EeFvSdwiiGdhWLazWTuuHFOQX4K/RcVlehMqA7NgyWDIkkm/C81pyGNf0AKTAgzxqF+zJm7ttb2S37yvk/AYLjBXTo41UywQwN/WxqyHH6ZuzRpO798faWjp9+PKzSXYac1wmXL77bzDQU6sf8/UJ6yDnKbXWVl9AUtnLqGwrIwXs/ez/48PMufp+5n7z9ex5Oa7gNANbHnkP3tzRQUNa37N7JyOmF2UpZTMPHJ62HdYTlZc8oSPDunBedCrIiADvLf/Pb62+Kth68O46Vo+h1xYS8cvi0ondhKkPjFIqbVJT9y0zGOR0TfrgXBFGdaOcRyIdpn1RSj64p4bQmiBSSGtIUukbeNGx1hKX7C73XJnzMDf2GixdNq3baNl/XqLteLKyWH00qU0v/WW5Xj+tjZm3fxlnnVXWSrzzf+R7cH8no8eZT1wXtF5FgE0xG9yVxe3eAT/Mz+HHg9R7WMkkv1nZQ94a5lMb6hpZIfFwyuC/KRoH7iy+U3nJ6k/3U1bd6vjfgEZ5PebniMgA7xZ9aYl0P9m1Zuqqr3x3nC9xdILn2dpWWQkklNn4z71FUt2XnwfSXqUgKVmRVhv1ubYie9E4haGeX+zOJn3zT0buhsh0KomifZFKOyClQHWy0CgBSZFOMVhEhGYWAkCdrfbxBtvZOKNN4Zv9O3btnFs9WrLsbImT2bGPfdQWFbGnuXLaX77bZAyLEhnh3z1r297g02HN+ML+i3/kddXb2CBzRo59NzvcR09jejuCZ+f+Vqz/JJRAS/PlnnCN/2d03MAidfloXnBFP4gDnHOke4BEQRzLU5vLrdEbvQDzZwJswHBnoY9Ueu8Lg8zx51NYW4h10ydxNLsA6zrPpua99eqAL3byzljz+FA08FwHAUIT6CEGIH+g++xdIR6gl/XkkPluqcoWVgYUziSckUNwrz4pEcJmAVQjLNZL/1wRfUmTobgJGshxXHPrVy5ktWrV1NbW8vXv/51fvvb3yZ+3DSiBSZFJBKHcSKWMMVynRkWxLEnn4w6Vu4554S3m7NqlWPq9LyaHqbsGMXxmdext2Efc2p8zK3pgZmwZGYpm6fsslgcHT2diG5rEoL5Wns8gh2TpK2LsnIU+oMBjrUc49g074ClL/dlqFlJ8cVU1e+i5XSL4/pUcMrhu6YUTOGS6SVRN811ATnH+gAAIABJREFU1etZvfWJsFj4Aj7mT5rPty65lcqaTeRn57G3YR8nO5upDomOEeg/fPIwvqAfr8tDyYzLoPGvrGvJYUXdHLpkE2/W3RdTOJJ2RQ0CfZ4jY2AI4K5d1uX9cUUlKk6JxnESZPPmzfzkJz/hpz/9KZdffjkTJkwYsGOnGi0wKSKZVjEQX5icYjagRMlcRwOqnYxR8R9rf7O15Ha5mBtyre372/vMevhh7iy7nW8equRZtoetEYBZteqmbpyfca1/X/Nr/pxzOHyDH5U9itbu1vD3yV66AQtEr9sUjZxIU0dTuH1Kot2c3S43NScPJ/w9yVI8Zlr4ewAa2xujikwb2uodxcXcJQGsWWGGy3JjzSYCMoDX5eGT0y/hGmMMshGMFQLGXAxnPU/luqfokiqlPZ5wDMRUyyFDf1xRaYqT7N69G4Dvfve7SacLpwudppxCCsvKOPvee/sUezFu1hNvvjnhuI2lSabLRd68ecz+5S8BVZTZXFHhuJ8ltdkUtzGnR3/rklupPmcUr5Xmh62SF64oJPj5f7KcX2FZGY1fX8buGSolNMeTzbXnX42rD//E5kyYhdvljrvN3KK5PPCF+7nu/C/xiUnzqZqezbNlI1k3LyduN+dAMMCxluO0nG5BIplSMAWvOzErSiCYUjAZj8v5eSzLnYUIBebr2+r5xKT54XVOHQyc0oPNVgSoVjFmi2Nd9Xp+t/H3YdeYL+hnUsFkls5conrEmVKRK2s2QdFVlCz8YUINKc1pzYNe3Z8OvKNVAWVfrYy+pDEPEN/85je55RY1Ab6goAAhBO+++y5NTU3ccccdTJo0iZycHObMmcOqVavC+wUCAX72s58xe/ZssrOzOeuss/jmN7+Z8vO1k3YLRgjxOeCXqA5yT0opH4yx3WLgA+BGKeXLoWWHgDZU5zm/UU0qhCgE1gDTgUPADVLKk9FHzUxiWSrxtrdbS2brpH7NGqbcdhvTli+37Ge2liwIEbaczH7w/Ow82rs7KLlqMaUOQeI1W14MtTpxc+PCG5hXdB4jskbQ0dP7kDOv28u/XHIrD/99FfVt9TG3+/u+d5k2ZiolxYt57ePXgeSGmrV1tXHzRTfR3t1BVe2OmHPtJ46cyMyxM8j25nCsJXrUdI4nm4VnXWgpfjxy8kjUdvZ97Dd7uxVhH3P82rY3LFaXS7jCx4hlgfQlhpG0K+pMY4DdX71x7733MnXqVP7zP/+TiooKRowYwdy5c1myZAkNDQ385Cc/4dxzz2X//v3s378/vN8dd9zB7373O370ox/xqU99iubmZl5++eVBO2+DtLaKEUK4gb3AFcBRYBPwFSnlToft/gp0AU/ZBGaRlPKEbftfAM1SygeFEHcDY6SUP453LuludhmP3trOOGFvJ4PbzZxf/Spq/+aKCo7++td0VFWFlxV+9rPMMT0NJYK9X9a5E2azp2FflCvK7XJTNHIitS11BFFP91MKpjB7wiwONx8mNyuX7bU74n7XqOxRTC4oiikKiWJuxXLPn+6N6locq5uxgdHVGLD0+Mr25jjGebxuL4umXhSzyWS8TK4fvPZ/hUUM1O/7tGkUcuqnSw5Ndu3axdy5c/t/ICNwf3I9NP9jUDsY/Pa3v+XWW2+lra2N/Px8Vq9ezbe//W22bNnCBRdcELX97t27mTt3Lr/85S/53ve+1+/v7+03zORWMRcD+6WUBwCEEC8AVwM7bdv9K/AKkKhz+Grg8tD7Z4B3gbgCk6kk0nbGidGlpdSvWROJzQQCjplshrV0OJQEUFhWFmXpJIK9Pcnexv2OcY5AMECOJycsLqDGLx9rORb+/IlJ8+ns6SQ3K5eddbvCbiGD1u5WWhta6S9GXOKHZd/ngS/cz2ubnmBz3X580oWbIEW5BRzrcE4ImFIw2dIy32wpVNXt5JmNz4a3/adZl1sKI2MRz4q4ZsHVlrb7/xIStkT21fQTI3vsxN9hz/8Nwa7YkzEHgYqKCi688EJHcQH4+9//DpAWl5iddAvMFMDsTzgKlJg3EEJMAa4FyogWGAm8LYSQwGop5eOh5ROllLUAUspaIYRj2oUQ4nbgdoBp0+I0tUsjiaY7262cwrIyptx2m6Ujc7xMtmnLlyclLAZmd4y5X5YdI+vpUPOhcMqwXYga2ht5/TZlzpu7N6eCHbU7Iv24alZDy0kEMCu7g+dORcePXMKFAI61HOd///Eevrb4q9xZervlBm+8rq/ewJKZpdxZenu/z3PpzCU88PmVvLbtjThlmJqBoK2rjY6eTvKychmZMzKSPXaqUokLqBqb2ldh7GW9u8ySTVs29rONAWhqamLSpEkxd2tqaiIvLy8jEgLSLTBO/1fsj72rgB9LKQMObQtKpZTHQwLyVyHEbinle4l+eUiQHgflIuvDeQ8aiaQ7x7Jypi1fTv6CBX12ryWLU3sSl3Axa/w5LDs1juxt1Uz81OUsCd1wjWJBe1bXmBEFlmMCURlWZkbnFHCqK9rScAsX54yfRV1bXczU5D0Ne1mx9j5VrLjjKF2BseSIAAg3PtO/iDkTZjN/0nxLR2Kj4HFe0XlR1sOdpbc7Ckt/XVlbQ7/rlqNbkw/ID5FGlemgrauN4y3HCSJpOX2KyUxmpJE9NroE6v+oRMaVA6M+oSybeAH/ZAs7zft1N1hWjR071hJvsTN27Fg6OjpobW1Nu8ikW2COAlNNn88C7JHURcALIXEZB1wphPBLKV+XUh4HkFI2CCFeQ7nc3gPqhRCTQtbLJKCBIUoi6c7xrJy+JgwMBPbg8ryaHvY9FGpzs7Ga5qLzuLPs9nDLkvzsPH6/+TkCwQBul9vi/jFuyDcuvIEPDm1kb8PeKIsnLzuPLn8XXf5u3MIddqkFZJDC3DFcMv1ii8vKTrhYMaBiLV3SjRy1EG/XkbBLynCHPbbhcUvL+4AMJFw38tiGx/n9pj8QkMG+9dYKMSD1KuY+XWl082QqHT2dBEP/voJIOno6GTlqohKGSeMgZxI0vK3EZeynei/WTLaw09JM0/rv/dOf/jQvvfQS27ZtY8GCBVG7loX+v//ud7/jrrvuSui6U0W6BWYTMEsIMQM4BtwE3GzeQEo5w3gvhPgt8Ccp5etCiDzAJaVsC71fBqwMbVoOfAN4MPT6RqovJJX0JhLJFnWmErPL6MBz9zsKoNniKZm2OGrKork2xOv24g/4HeM6dW31nDNuJmNzC5k1YRbPbX4+HJjfeHgTB5sOOp6jQP3X9bo8LJlZyrGWY+FAfU7+NAK1h9SGpkSY9u4O2zEEx1uOR7W9dxoCpiw2ddPoTSCcLJ0BqVcZpEaVQ5W8rFxaTp8iiMSFIC8r1ObfyB7LvRkmXBmxLnqrh3GqnUnEZWbez+bo+frXv86jjz7KsmXLuO+++5gzZw4HDx5k7969PPjgg8yZM4fbb7+dH/zgBzQ0NHDZZZdx6tQpXn75ZV544YWB+aESJK0CI6X0CyHuAt5CpSk/JaWsEkLcGVr/WJzdJwKvhSwbD/CclPIvoXUPAi8KIb4FHAa+nKpryASSLepMBHtsJ5mMtuMzxxDwCLL8kh6P4PjMMZwdWmcWEbdwMXvCLEvzRnNVeywCwQB7GvaS48nmmgVXs3jaonAMKBAMUOeQ9vyJSfPZWb+LQDBAEMm8ovMs6di/2/j7sJgZdStLZy6hpHgx5Tv+iC/gQyBwCcH7hz5gq8ll5WSpVIYKJA3cpjRjO7FatyTdOsXMMJ8B319G5oxkMpOtMRg7fSnWjNcUM57LzLxftjWEnJOTQ0VFBXfffTcrVqygtbWV6dOn853vfCe8za9//WuKi4t58sknefDBB5kwYQJXXHFFMj9Jv9ATLUNkcppyumiuqGDv97+P7OlBZGUx+dZbo0Y7JyIyD1U8wu43ng93A+i6aG64zsOe3uwWbn72xf+MGvHrdXtBSnxBP27h4rJzLuP9g+9HxWWuO/9L5GfnxXWJXTr9EgRYEhHM0x3tKcEu4eLBLz4QPicjpdmFCLtTQBVHLplZGo4tmc+ppHixSUjha3Mv4c7PWqdJGlaLffLkded/iR+amlX2mzMoBjNgacoDRV+nZmYAQzlNWTPAJGNhxKJ+zRpkTw+gJl02/vGPSTXwLClezJvnvEnV9JAYNB1gxVrVG6ukeDGvb3s97DoyxzTsT+xAlNvJqVFnvCFaOZ5srl1wNb/54GnLchnjPcCkUUVha+G1bW+E3W9BpCXmU910wJIdBxFLZenMJay89Eoqt68mX3TRfryBdZsnsXTRvwPR7kCvyxPVfNSg3zUvg9Co8oylN/fXMG3LHwstMMOIZGtmYh6v86Sl0Yt3zBj8zc19jvUYQmFONzbXoHxt8VfDT/1u4SY/O8+yr/kman9vCI39hmuuyVk8bRGzJsxSXQiKF/PnXW+xx1Sk6RYurjV6eqHiQBtrNoaFoqG9MTxRcfORDyP7udycM24mp063hLsPGJ0MjGv52uKbI+nL2QdgxIlQA0o3b254nZVjLg63ezG7Az85/RImFUyO27MsmUQBTQpJJGNsmLblj4UWmGFEojUziTwBr6tez3Pjj/EVF3iDID0ezgr5eO0WUiLHs6cbm5/M7wynLf+BgAywZsuLjqm/sb7HSYTssQpj3z/veou/73vXctyS4ouj9i8pvjjsJvMFfKze8AQTR06wxIKCwSB7GvZG9U+zT5wMM34Zlaf/SpdU23cFpSW2Yw7gXxOj2j+TOx+nlWTcfgMxeMxyvPR0W85ktMAMIxLJJrM/Ad8/9homV5+McqlV1mzio6kufJ8ZyexjPsYtuYxPmlKfYx0v3hN1vCB1e3dH3Ayr3r7HLj5m0XHqVGxm1oRZjsvMcZjqpgMcPnkYr9uLL+CzpkMHrZ0GCnMLnWMmRVdR8on9vLnhdbqCMqm+YWdU5+NESSb1ur+Dx5w4w9xfiaAFZhiRSDaZ+Ql45v5WeHI19T3+KJdaSfFiyreXUzUd9p6dxwNfcE7E6+sTdayWJr3dOO3fs3rDE+Hj9SY+9k7FdpwsJns6MqhsskunX8LkgsnkZ///7d1/kBTlmcDx7wMsLtlI0IWBDeqiHGLIXQywEnWRkA2oIRaKMRCgPMmlROtCOCIYOS05SgpDuEBCTk1CQBRLoTxU9IyXxApnVQQMIomIUoSQcwmwsgiueOAqsM/90d2TZuj5PT3dA8+namt2+sfMM729/cz7o9+3JjmLZFWnLiCSnJ8l00X/qoY7uO+cYTmVxAL3L0VPstNNIV2vs5U2Uks3OXUtPrOqv3JhCeY0k+2eGf+FfNC+Djp/7DZYB1Wp+ecYyeH1ivlGne3C6X8fcEoUXkeBbEkuNcYrLryC13ZvSc5Vk20fj9dBwNvOP7c9kPNFv9hxw2zcsRSFdL3OVNpILd107Q0f78+ttOOv/ip1FVwFsgRzhvFfyIf1rqLTn5cHVqkFzTESdFEr5TfqTBfOTB0F/MnA6yTgVZl5UwxMGDI+2cifWuoJSoyB0xTk0PZTKjY6ch78UyTn2gaTqbSRWro53pb/3fhhVMFVIEswZ4igNgqAQ30GBVap5VMyKcc3ai/+1Dvuvc8zYcj4ZCeBJ15bnbxnxuMNze8fkDJbYiz158o1aVhPsQIU0vU6XWN7aummSw/o+Ci/tpVCh4g5zViCqWC53vOS6YKVrkotTnX9qaUNrzTyybNqkve8+DsJBN313378I9ZufTZjR4BSfNZ0r5NP0oiqp5iVmlxBpZsuNflVd1mDP2AJpmLlc89LoResuNT1p8a/s3UninNPyrETx/jlm79kwpDxVHc565S7/j1VnauSN2SmXuBLVWLI9Dr5/A1Omlunc1XgWGelZqUml7/dxH+Hfb5di63BHyCPCdNNrATd85LOF+ovy2ludr9D69bxl3nzOLRuXWkCLoI/fi9RbHz7lWRJpf34R/zfR0eS88rP/+p9zL9uHl+79EZuGXYzX7v0RhrOH5pMON4F3hN08S9EptfJ52/glR6v7Hc5qLLx7VeY88JcfrZ+KT9c96PkTZ+lVKpjUNG8dpOPW53HY23FvV5VDydJFZlcpkyZQkND4EgssWclmAqVzwjK+VZ3lXpEgGJ5bSwv71pPl85dTroLHzipLSZdo/vvdr2cnEsl9QKfb0+4dFVJmV4n37+Bd3e/PykWM9R/NnZ/DdZuEgJLMBXIa3upu+UWjn/wQU7jjuVT3ZXriABh+t2ul5OzNw5IDPjbPSe+cbq8YWDS3fXul+kCn8/FP1NVUrbXybfKMbV3nHdjZxjtMnFqc4tM3NpNvOq6jo+jjaMIlmAqjL90kc+IxvmIen4Z/4jFAK80b6LD14CfbpyubLJ1g87ltbK1pZSy3Sq1q7SXZMMqYcSlzS0ycWo38XdzPv4B6AlefPFFZs6cya5duxg8eDA///nP+exnPwvAiRMnWLhwIStWrKC5uZlevXoxatQoHnnkEQBGjhxJz549ufrqq7n//vvZv38/TU1NLF26lL59+4b2MSzBVJhylC7CnF8mF/6qIYAO7aCzdOKEdmQcpysXxfaUKndVkv+i77+x84xOBGGKyzhhKTNa7t79V+68807uueceunXrxqxZsxg/fjzbtm1DRLjttttYuXIl3/ve9/jiF7/IoUOHWLNmzUkvuXHjRnbs2MHixYtpb2/nrrvu4oYbbuDVV8Nrb7MEU2GCShelHKLfE8VUyx5vmBovyVR1rmL4RY3sPrSb4f0bi0ouxfaUirIq6YwvYcRA2bpyp8xoeei9NtZv2MiAAc64eR0dHYwbN44dO3YAsHz5cpYsWcL06dOTLzFhwoSTXrK1tZUNGzZQX18PQH19PcOHD+dXv/oV1157bSgfwxJMhUktXQCxapAvhav6D2f+dfMC22D2vr83cKTlXJTq/pJCLvRxvcckrnHFUVm7cvur67qcTb9+/ZLJBWDQoEEA7Nmzh507dwJOb7NMhgwZkkwuAI2NjSQSCTZt2hRagrFuyhXo3KYmLrr3Xs5tasqru3Iluar/cBaP+3du+Nz1vLxrfUm60BbSXTsv7zwHb0xzHn28C9NTrz/NnBfmhtLNuBBxjSuuyt6V2+vm3KkrPXqcXG3XtWtXJ472dg4ePEhNTQ3du2fulJBIJAKXtbS0lC7mFJZgKlyPxkY6VVcDRNIgXyq/2/XyKfd4eBdAb+wxKC4xeNVbX7v0xtJ/+/SGjG9+0Hn0JZm43mMS17jiKvQvKAWqra3lyJEjHD58OON2ra2tgcvq6urCCs0STKXzqsx6T5pUsdVj6b5Jpw6z37/2oqITw1X9hzOr6bulr9oIGjLeFdcLU1zjiqtQv6AUocn9n1+5cmXG7bZs2cLu3buTz9evX09rayvDhg0LLbbI22BE5FpgCdAZWKaqC9JsdxnwCjBBVdeIyPnASqAP0AEsVdUl7rZzgVuBA+7ud6vqC6F+kAhF2SBfCunaRlJ7bN3WeGts/qlPkWHI+LjeYxLXuOIsjh0tBg4cyNSpU5k5cyatra2MGDGCtrY21qxZw+rVq5PbJRIJrrvuOubOnZvsRTZkyJDQ2l8g4gQjIp2BB4HRwB7gVRF5TlXfCtjuB8CvfYuPAzNVdYuInA28JiIv+vb9kar+MPxPYYqVrutvpVwAnYbyXXyh1zyuOusvgUPGx/HCBPGNy+TnoYceor6+nmXLlrFgwQISiQSjR48+aZsrrriCUaNGMWPGDA4cOMDIkSNZunRpqHGJqob6BhnfXOQKYK6qXuM+/1cAVf1+ynYzgGPAZcDzqrom4LWeBR5Q1RfdEsz/5ZNgGhoadPPmzQV/FlOcSu3NlDrSc5yqTgxs376dz3zmM1GH8TcRTULm3WiZem9MLrIdQxF5TVUDB0uLug2mL/BX3/M97rIkEekLjAN+lu5FRKQfMBj4vW/xNBHZKiIPi8g5afabKiKbRWTzgQMHgjYxZRJa20jIrKHc5KzUg2lWgKgTTNBcvKlFqh8Dd6m6AzGlvoDIJ4GngBmq6nWj+CnQH/g80AIsCtpXVZeqaoOqNvTq1auQ+E0ZBPUwi8v7WEO5yVnQYJqnuagb+fcA5/uenwfsS9mmAVgtzrzwPYExInJcVdeKSBVOcnlcVZ/2dlDV/d7vIvIL4PmQ4jchK9fNbYW+T6W0E5kYiHAwzZdeeqls7+UXdYJ5FRggIhcCe4FvAJP8G6jqhd7vIvIIThvMWnEyznJgu6ou9u8jInWq6t09NA7YFt5HMGEq1+yOxbyPNZSbnMRpMM0yibSKTFWPA9NweodtB55U1TdF5HYRuT3L7o3AzUCTiPzR/RnjrlsoIm+IyFbgS8B3w/oMlapc1U7FKlcVlFV1mbIo0SRklSLSXmRxcib1Iqu0nk/l6mFWET3Z3nnOuYkzoCu0OVXsepFVoGJ6kUVdRWYiUK5qp1IpVxVU7Ku6vOFoThx1buocvMqSTD4i6iJ8Jou6F5mJgFUHVagMw9GYLM7ALsJxYCWYM5D1fKpQGYajMVkEdRG2UkzoLMGcoWJfHWRO1WesUy1mbTD5i7CLcLGmTJnCtm3bCGwjjnm1nyUYYypJn7GWWApxOnYR9qr9tMNJnp+4KHafyxKMMebMUNUjdhfgolRAtZ818htjTAVYu3Ytl1xyCdXV1QwfPpy3du5zqvuAEx3K9xcv4+KLL+ass87ivPPOO2UK5WeeeYZhw4bRrVs3amtrGTNmDM3NzaHGbAnGGGNirrm5mTvuuIN7772XJ554gvfff59rrvs67Z0+DV0T3DbrP/i3+xYwfvx4nn/+eRYtWsSRI0eS+z/22GPceOON9O/fnyeffJIVK1Zw8cUXE/Ygv1ZFZowxOTq0bh1t69fTo7GxrJP8vfvuuzz77LNceeWVAAwdOpT+/fvzyONrGTlyJMtXPMaSJUuYPn16cp8JEyYA0NHRwezZsxk3bhyrVq1Krh87Nvy2PEswxhiTg0Pr1rFz5kw62ts58PTTZZ2iPJFIJJMLQH19PUOHDmXTpk14o7GkVol5duzYwb59+/jmN79ZjlBPYlVkxhiTg7b16+lobwego72dtvXry/beiUQicFlLSwsHDx6kpqaG7t2Du14fPHgQgLq6ulBjDGIJxhhjctCjsZFO1dUAdKqupkdjY9neu7W1NXBZXV0dtbW1HDlyhMOHg+eXqa2tBaClpSVwfZgswRhjTA7ObWpiwKJF9J40qazVY+Akkw0bNiSf7969my1btjBs2DCa3DhWrlwZuO/AgQPp27cvjz76aFli9bM2GGOMydG5TU1lTSyenj17cvPNNzNv3jy6devGnDlzSCQSTJkyherqaqZOncrMmTNpbW1lxIgRtLW1sWbNGlavXk2nTp1YuHAhkydPZvLkyUycOBERYd26dUycOJGGhsCBkEvCEowxxsRcfX09d999N7Nnz6a5uZmGhgZWrVpFtVtl99BDD1FfX8+yZctYsGABiUSC0aNHJ/efNGkS1dXVzJ8/n5tuuomamhouv/xywp4q3uaDcZ1J88EYc6aw+WCKV8x8MNYGY0yU3nkO3pjmPBpzmrEEY0xUvAnEmh90Hi3JmNOMJRhjomITiJnTnCUYY6LS62pn4jCwCcSysarEihR5ghGRa0Vkh4j8WURmZ9juMhE5ISI3ZdtXRM4VkRdFZKf7eE7Yn8OYvHkTiNV/23m0eV6CFVmVaB2ZClfssYs0wYhIZ+BB4CvAIGCiiAxKs90PgF/nuO9s4LeqOgD4rfvcmPjpMxb+4QFLLpkUUZVYVVXFhx9+GFJgp78PP/yQqqqqgvePugQzDPizqv5FVT8GVgPXB2z3HeApoDXHfa8HvNtWHwVuCCN4Y0wZFFGVmEgk2Lt3L0ePHrWSTB5UlaNHj7J3797AcdByFfWNln2Bv/qe7wG+4N9ARPoC44Am4LIc9+2tqi0AqtoiIoFHSESmAlMBLrjggsI/hTEmPF5V4oHfOMklj9KeNwDkvn37OHbsWFgRnpaqqqro3bt32kE0cxF1gpGAZalfM34M3KWqJ0RO2jyXfTNS1aXAUnButMxnX2NMGfUZW3A1Yvfu3Yu6SJrCRZ1g9gDn+56fB+xL2aYBWO0ml57AGBE5nmXf/SJS55Ze6ji5as0YY0wZRN0G8yowQEQuFJGuwDeAk7qIqOqFqtpPVfsBa4B/VtW1WfZ9DrjF/f0W4NnwP4oxxhi/SEswqnpcRKbh9A7rDDysqm+KyO3u+p/lu6+7egHwpIh8C9gNfD3Mz2GMMeZUNtilywa7NMaY/Nlgl8YYY8rOEowxxphQWBWZS0QOAM1Rx5FBT+DdqIPIwOIrjsVXnLjHB/GPsdD46lU1cOYySzAVQkQ2p6vnjAOLrzgWX3HiHh/EP8Yw4rMqMmOMMaGwBGOMMSYUlmAqx9KoA8jC4iuOxVecuMcH8Y+x5PFZG4wxxphQWAnGGGNMKCzBGGOMCYUlmAhkmyZaHD9x128VkSG+dW+LyBsi8kcR2exbXrJpoguNT0QGunF5P4dFZIa7bq6I7PWtGxNifJeIyEYR+UhEZuWyb5mPX2B8InK+iPyPiGwXkTdF5F9860p2/IqJ0V0Xh3Mw3TGMyzk42f3f2CoiG0Tk0mz7lvn4BcZX8nNQVe2njD84A3PuAi4CugKvA4NSthkD/DfOnDeXA7/3rXsb6BnwuguB2e7vs4EfRBFfyuu8g3MTFsBcYFaZjl8CZ3K6+f73zLRvmY9fuvjqgCHu72cDf/LFV5LjV2yMMToH08YXk3PwSuAc9/eveP8jMToH08VX0nPQSjDll8s00dcDK9XxCtBDnHltMinVNNGliu/LwC5VLfXoCFnjU9VWVX0VSJ3CsBzTbBccn6q2qOoW9/cPgO04M7eWWjHHMJPIj2GKKM/BDar6nvv0FZz5qrLtW87jFxhfqc9BSzDlFzTVc+ofMNM2CvxGRF4TZ8pnz0nTRON8w4siPs83gFUpy6a5RfKHiyj+5/LeheyEaU8aAAAEK0lEQVRbzuOXlYj0AwYDv/ctLsXxK0WMcTgHcxGXc/BbOCX+bPtGdfz88SWV4hy0BFN+uUz1nGmbRlUdglOs/baIjChlcFneO6dtxJkAbizwn771PwX6A58HWoBFIcYXxr65Kvo9ROSTwFPADFU97C4u1fErRYxxOAczv0BMzkER+RLOBfyufPctQjHxectLcg5agim/XKaJTruNqnqPrcAzOMVhcKeJBpDipokuKj7XV4AtqrrfW6Cq+1X1hKp2AL/wxR1GfIXsW87jl5aIVOH8Yz+uqk97y0t4/IqOMSbnYDaRn4Mi8jlgGXC9qh7MYd+yHr808ZX0HLQEU35Zp4l2n/+jOC4H3lfVFhGpEZGzAUSkBrga2ObbpxTTRBccn2/9RFKqJlLaaMb54g4jvkL2LefxCyQiAiwHtqvq4pR1pTp+xcYYl3Mwm0jPQRG5AHgauFlV/5TjvmU7funiK/k5WEgvBfspuhfKGJzeGbuAe9xltwO3u78L8KC7/g2gwV1+EU6PkNeBN7193XW1wG+Bne7jueWOz133CeAg8KmU13zM3Xare7LXhRhfH5xvcYeBNvf37un2jeD4BcYHDMepytgK/NH9GVPq41dkjHE5BzP9jeNwDi4D3vP9HTdn2jeC4xcYX6nPQRsqxhhjTCisiswYY0woLMEYY4wJhSUYY4wxobAEY4wxJhSWYIwxxoTCEowxxphQWIIxxhgTCkswxsSAiHQRkeki8rqItIvIOyLygIh8QkTeF5G3oo7RmHx1iToAY8507nAe/4Uz7Mpm4CdAT+CfcO6c7w48H1mAxhTIEowx0XsAJ7ncqao/9BaKyKPAS+7TLRHEZUxRbKgYYyIkIpcBm4DfqOo1Aeu9mQm/rKrryh2fMcWwNhhjojXNfbwvzXpvGPU/pK4QkbtFREXkgVAiM6ZIlmCMidY1wEFVXZ9mfV/gf/Vv09sC4E6TcCvOyLbGxJIlGGMiIiLVQG9gd5r1fw98mpTSi4h8CngcZybC9wJ2NSYWLMEYE50T7k9tmvVz3MfUBv6lwBprkzFxZwnGmIio6jGcyaUucOdGB5xZBUVkDvB1d9EffOtuBf4OuLecsRpTCOumbEy0FgIPA78UkVXAIWAUcDbwFjAItwQjIgOB+4GrVPXjaMI1JnfWTdmYiInIHcB3cNpb9gFrgPk4pZtjqvppd7spwAqcajVPZ5wpbjuAGlX9qHyRG5OZJRhjYkhEzsdp/H9BVb/qLusBnJey6QqcRHQ/8KbaP7SJEasiMyaeBruPyQZ+VW0D2vwbicgR4JCqbitjbMbkxBr5jYknL8GccoOlMZXCSjDGxNMpJZggqjoy/FCMKYy1wRhjjAmFVZEZY4wJhSUYY4wxobAEY4wxJhSWYIwxxoTCEowxxphQWIIxxhgTCkswxhhjQvH/uU7uCNjoTJYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(fccq[0], fccq[1], s=10, label='fcc', color='#FFB300')\n",
    "plt.scatter(hcpq[0], hcpq[1], s=10, label='hcp', color='#388E3C')\n",
    "plt.scatter(bccq[0], bccq[1], s=10, label='bcc', color='#C62828')\n",
    "plt.xlabel(\"$q_4$\", fontsize=20)\n",
    "plt.ylabel(\"$q_6$\", fontsize=20)\n",
    "plt.legend(loc=4, fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The thermal vibrations cause the distributions to spread, but it still very good. Lechner and Dellago proposed using the averaged distributions, $\\bar{q}_4-\\bar{q}_6$ to better distinguish the distributions. Lets try that. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcc.calculate_q([4,6], averaged=True)\n",
    "fcc.calculate_q([4,6], averaged=True)\n",
    "hcp.calculate_q([4,6], averaged=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "bccaq = bcc.get_qvals([4, 6], averaged=True)\n",
    "fccaq = fcc.get_qvals([4, 6], averaged=True)\n",
    "hcpaq = hcp.get_qvals([4, 6], averaged=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets see if these distributions are better.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f5cfc125650>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEPCAYAAACDTflkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3de3jU9Zn//+edyZFTkMjJoAEpoNZaxYBoUFlU9KteqLvWU2vLfn8t+GtpF6VrrS7o166r61VbbasrVNYVrWJLrVJ1q/v9RbcQVyRgpSKgHAoGgaSJ4RxymPfvj89MmJlMTpP5zEzC63Fduch8TnODY+68T/fbnHOIiIgkIivdAYiISO+lJCIiIglTEhERkYQpiYiISMKUREREJGHZ6Q4g1U488UQ3evTodIchItJrrF279q/OuaHxzh13SWT06NFUVlamOwwRkV7DzHa0d07dWSIikjAlERERSZiSiIiIJExJREREEqYkIiIiCVMSERGRhCmJiIj0dXtWwJ/nen8mmZKIiEhfEpsw9qyA92+GHY97fyY5kSiJiIj0FZvugcrrvISx7itewqh5E1oOe+dbDnuvk0hJRESkL9izArY8BAS918FG2LEIhs6AQD/vWKCf9zqJjruyJyIifUq4tXF4O60JJNKImXDOC941Q2d4r5NISUREpLcKj3e0HIasXLBscM3Hzh/d610zYmbSk0eYurNERHqryPGOYCMM/DIMOvfY+f1rj42N+ERJRESktxo6w2uBhB34AJrro68JNiZ9MD2SkoiISG81YiYM/NKx164ZDu9se12SB9MjKYmIiPRmucNjDjTFvA7Armd9e3slERGR3qzw7LbHsgoiXrTA7uWw9iu+vL1mZ4mI9Dbhab1DZ0DTvrbng0faHqsp9yUUJRERkd4kclpv1dMwKE5LJJ6h030JR91ZIiK9SWwZk33rO7+n/2lw7m98CUdJRESkN4ksY0IWBBs6v+fQJq+ulg+UREREMlW8Eu4jZsKYeYDhlTlpbufmGFseUil4EZHjRnsl3PesgJ2/BFw3Hxj0CjImmQbWRUQyUXsl3MOD6olo3Juc2CKoJSIikonilXDfsSjxBAJxFib2nFoiIiKZKLaEe/1qqPlD4s/LyoWSOcmLL/zYpD+xm8zsCjPbbGZbzOyuOOenmdk+M/tT6GthzPmAmb1vZq+mLmoRkRQYMRO+9Avv+y3/Stz9Qrpi0Lkw8Te+lINPa0vEzALA48BlQBWwxsxWOOc+irl0pXPu6nYe8w/ARmCQf5GKiKTRjkVAS+L3H/ggaaHESndLZDKwxTm3zTnXCCwDrunqzWY2CrgKeMqn+ERE0mvTPVDzRs+e4Zp9Kwef7iRSDHwa8boqdCzW+Wb2gZn9p5l9MeL4o8CddNLGM7PZZlZpZpU1NTU9DlpEJCX2rAh1Y/WgFRKWU9jzZ8SR7oF1i3MsdvLzOqDEOXfQzK4EXgbGmdnVQLVzbq2ZTevoTZxzi4HFAKWlpd2dXC0iklrhAoufv0tSEgjEL9SYBOlOIlXAyRGvRwGfRV7gnNsf8f3rZvaEmZ0IlAEzQ4klHxhkZs85576WgrhFRPwRWWAxmXxqiaS7O2sNXqtijJnlAjcBUevyzWyEmVno+8l4Mdc6537onBvlnBsduq9cCUREer3IRYbJ5FNLJK1JxDnXDMwF3sCbYfVr59wGM7vNzG4LXXY98KGZfQD8DLjJOacuKRHpm6IKLCaRD6vVAex4+3lcWlrqKisr0x2GiEj7Nt2TvAH1sJyhcHl1Qrea2VrnXGm8c+keExERkbDwgPrh7SQ1gfhISUREJBPsWQHrvgLBRrBs78t1scx7V5R8K3nPiqAkIiKSCXYs8hIIJDd5hA0+L/nPJP2zs0REJBU23OHLY5VEREQyQckcyMrz7/lHtvqyRa6SiIhIJhgxEyb+Gkq+AyOv9+c9fNgeV2MiIiKZYsTMY+XaK8rg83eS//wkUxIREckUe1bAx/fDkZ3QVJvcZxeeC6c9kNxnoiQiIpIZ9qyAtX/nz8wsAjBuYeeXJUBjIiIi6bZnhTfo7UsCAUZe50tXFqglIiKSXn5V7Y10aLtvj1ZLREQknfyq2hvpyE7fHq0kIiKSTkNn+Ls+BKD5c1+m94KSiIhI+uUO9/f5Pu6xrjEREZF02bMC1t0AwaM+v1HAa/H4QElERCSVwuXeh84IFV30O4EA2QM1O0tEpNeLnIlV9TT0Pz0179tc7723D4lEYyIiIqkSOROr5TDkDYes3NS8945FvjxWSUREJFVyCqNfF54Np36f3vyjWN1ZIiKp0rQv+vW+P4XWcAT9fd+sXK/UvB+P9uWpIiLS1tAZEOjnfZ+VB399Ew5+6O97DjoXJv5GA+siIr1S5Gws8AbTj+z0ZmUF9/v//vv/5OvjlURERPwSORvr06fAtfhXZLFdLV4SU0tERKSXiZyNlYr1IHH5t9AQNCYiIuKfoTNSN4W3PV/4gW+tEFASERHxV+7INL73CF92M4ykJCIi4odN90Dl30LDjvTFkF/sjcv8ea6q+IqI9Bp7VsDWh4GWdEfiDezveNz704dEoiQiIpJsNW92MgsrACOvP7ZmxC9HdkaXWfGhHLxmZ4mIJNvQGV6BxfZ2LAz0h+bDMGaet4o9pxBq3oB9a5MbR1NNxHv282WWVtpbImZ2hZltNrMtZnZXnPPTzGyfmf0p9LUwdPxkM3vLzDaa2QYz+4fURy8ix7X2xhtGzIRzXoCS73gtDiz6fMt+qHnd6/Kqfxeq3/B1C1vyT/Hi8WGWVlpbImYWAB4HLgOqgDVmtsI591HMpSudc1fHHGsG5jvn1pnZQGCtmf1XnHtFRJIvtqx77A/p8Pfv3wy4+M9wzclvfcQz8EzfpvmmuyUyGdjinNvmnGsElgHXdOVG59xu59y60PcHgI1AsW+RiohEii3rHm+8IfKatAn4VnwR0p9EioFPI15XET8RnG9mH5jZf5rZF2NPmtlo4BxgtR9Bioi0EVlMMXK8IbKLa+gMsDQPPffxxYYW51hsu28dUOKc+zLwc+DlqAeYDQB+C8xzzsWtZmZms82s0swqa2pq4l0iItI9keMe4a6scBdXeEotwNg7Sf+PWv+k+29WBZwc8XoU8FnkBc65/c65g6HvXwdyzOxEADPLwUsgv3LOvdTemzjnFjvnSp1zpUOHDk3230FEjlcjZsKXfnEsgWy6p20X12kPwNAr0hejT4sMw9KdRNYA48xsjJnlAjcBUX9jMxthZhb6fjJezLWhY0uAjc65n6Q4bhGRY8ItkMi9QSK7uErmkLYftz52ZUGak4hzrhmYC7yBNzD+a+fcBjO7zcxuC112PfChmX0A/Ay4yTnngDLgVmB6xPTfK9Pw1xCR40m8ab07FkUPoA84M7qLq+ZNKBiT+lhTUDvLvJ/Hx4/S0lJXWVmZ7jBEpDeKnNYb6OclCoB1X4Fgo/d9Vh5M/HX0GEnLYe94sJF2p/tm9YfgoeTGO/RKOO+1Hj/GzNY650rjndOKdRGRrmpvWm84gQAUXXKsCyl2P5HCc7191ePV1Ep2AsF8ndoblu4xERGR3iPetN7IYwSg8Oz2rx+3EEpfgoKxKQ3bT0oiIiJdFW9a74iZXg0sywZaYPujx8ZL2rt+WCpmazlfCi7GUneWiEh3hBNBpKZ9x6r2hru5wtfEu76zAo3JkJXn67a4rW/j+zuIiPR17a1eD4ud0RVboDE3mevXAt7YS3hw32dqiYiI9FQ4KdS86SWQyB/e7RVqjCzQmNQWSQsc+LDzy5JELRERkZ4ItzLg2Or1SB0VavSrQGPwqLd2JQWUREREYnV1X/LYWlmR1+9ZAauv8vYLycr1jsV2dUXN7OqdlERERCJ1lBhitdfK2LPCW4BY87q3X4gLegv/4u05MiyZhTZCP9KzclOyRiTiHUVEBOjaPiFh7Q2o17wZvQDRNUO/MW3HSv48Fz5/L3mxD73CG6yf+JuUDKqDBtZFRKJFTr/tbF/y2AF18BJDTqHXGogshRL5nMjB9mT+Lt9c743LpJCSiIhIpI5mWrV3PXgD2bX/nzeoHegHp34/VOIEr2sp8jlRA+rB5MW+b72XxLoSd5IoiYiIxIq3QLA9Ua2KkJbD3gLE9oofDp0Bnz7lJZxkCh70xnLi7fnuE42JiIj0RLxpul3pBiu6xL+YOhvLSSIlERGRnogcXM/KjT8LK56SORHTe5P9ozgrJSVPQN1ZIiI9090xFDi2UdWYeV631+fvwv61yYup8ByNiYiI9BqJjqGEN7YaOgMqr6XdDau6w7K9kvMpou4sERG/Ra6Aj7cOZcRM+MIPScqP5EFfTlkrBJRERET8FbsCPqcw/gLF0x7wFgv21IE/d16uJYmURERE/BTb8mjad6wM/Jh53vnwD/2SOUAg5gGxrzsRbEzZzCxQEhER8Ve80igjZnp/bn80ukbXiJkw8rro+/uP6977WXbKZmaBBtZFRPzV3uyt9mp0HdgUff+hmNcdCsDYO1M6JqIkIiLit3izt3IK8bqqWrwWSk4hrLshZhV76HxHCsZC8Y1eN1kKy52EKYmIiKTanhVeVxYtXvfTmHlena3YMigFo+HIdjqsrzXsCm9QPvLZ3Vmz0kMaExERSbXIrizX7LUiDm5ue92RrXRaoLH+3WMD893ZCyVJup1EzOxCM5sQ8frbZrbBzA6Y2YdmdltyQxQR6WNiB9tzCuHIXxJ71r61XjfY6qu8SsJd3QslSRLpzvo34LvAZjP7DvAQ8DNgIzABeMjMAs65x5MXpohIHxIebA/vg77vT3Q69tGR4FFvF8WsvGP7mHRWBDJJEkkiY4Ftoe+/Ccx1zj0TPmlm64EfAUoiIiIdqXvbazFk5XoJIHgUyALL8rq5CEDuUGjc0/be/FOgcW/0OErwqFcAst+YlI2JJJJEDgAnAjuAk4A/xZxfB5zSw7hERPq2yHGRYGP0D38IbXL1f70EYtnePu3h8ZGsXDjz59HXhVsfsRtg+SyRgfXXgbmh798Cbog5fyPwcU+CEhHp8yLHRQhA4dnesfA4Rr8xx7bXdc1EDbDnn3xs2vB5r3l7qpd8J2UbUUVKpCVyF1BhZiuB1cDtZnYRx8ZEpgDXJi9EEZE+aMRMb2rv1oe9JLHtEcB5iaPqae9coN+x7q5gC63jJoe3wqZ7YPB5x6bzpnhv9bBut0Scc3uAicBK4H8BBkwGZgBVQJlz7j+7+jwzu8LMNpvZFjO7K875aWa2z8z+FPpa2NV7RUQyWtO+UCsDbzwj3PKIrbE18TeQXxx9b9VzKZ/OG09Ciw2dc/uAu0NfCTOzAN4A/GV4CWiNma1wzn0Uc+lK59zVCd4rIpKZhs7wWh0th6PHPSJrbIW7p+pXw5Z/OXZvy5H4JeVTrNOWSGjtx7NmdnuoVVCYxPefDGxxzm1zzjUCy4BrUnCviEj6haf6Dr0SLAAEj61gj00Ipz0AX7gbcoYCBk01x86laDpvPF1piZwMnAHcEj5gZtuB9/FmYr0PrHPOVSfw/sXApxGvq4Dz4lx3vpl9AHwGfN85t6Eb92Jms4HZAKecooljIpJBRsz0WhE1oam64RXssfas8NaTNNURtQPigDO9BJOGVgh0LYmMAJbjjX9sA2rxfoD/HfC34YvMbDdtE8vOTp5tcY7F7g+5Dihxzh00syuBl4FxXbzXO+jcYmAxQGlpaRL2nxQRSaLIbq14rYrILXWjBNKaQKBrSeRh4GLgYufcyvBBMzsTWAhcD3wEDAOuCn2BNx+ts+dX4bV0wkbhtTZaOef2R3z/upk9YWYnduVeEZFeob1y8WGRa0rCLDvlZd/j6UoSuR5YFplAAJxzHwI3hGZLfQ34AjAIOAdv9tbZXXj2GmCcmY0BdgE3EdFtBmBmI4C9zjlnZpPxxnFqgfrO7k23uvJy6isqGFxWxpDp09Mdjohksnjl4sMiWypZuVB0afxFhSmu4AtdSyKDgb+2d9I5d7+ZXQvc7Zy7C6+F8PuuvLlzrtnM5gJv4BXO/3fn3IZwEUfn3JN4Sez/NbNm4Ahwk3POAXHv7cr7pkJdeTmfzJ9PsKGBmpdeYtwjjyiRiEhiOmupQHSXV9XTKVt42JUkshn4m06u+b/AV/AWInaLc+51vFXwkceejPj+F0DcVTTx7s0U9RUVBBsaAAg2NFBfUaEkIiKJ66ilAvF3SsyQ/UT+HZhkZrd3cM0QvDpaEjK4rIys/HwAsvLzGVxW1nqurrycbT/6EXXl5ekKT0T6mnh7uaeAeT1DHVzgLep7HbgUr+vofufcuxHnLw6d3+OcG+tjrElRWlrqKisrfXl27BhIvDGRyG6urPx8dXOJSOe6Otbh05iIma11zpXGPddZEgk9IAdvdfg38abR1gI7gSK8ir0G3OGcezRZQfvFryQSmRwAyM4me+BABk2ezNGqKprq6hg6cybNBw6w9/nnW+8bfsstnLpgQdLjEZE+InKsI9AvLUUWO0oiXaqd5Zxrcs7NxlvM9xxeFbCJeNNqNwP/uzckkGSK7ZKKHAMBoLmZ5s8/p+6NNzi0YQONu3eza9Eimmpr2+3mEhFpI95YRwbpVu0s59wa4BsAZpYHtDgXrh52/IhsdVQvX86gKVMYcPrpXbr3yLZtjHvkEU39FZGu6WwhYpolVIARwDl3tPOr+qbIVodrbGTfH//IvlWrunRvwamnMmT6dCUPEemarkzvTaNENqU67g0uK8Nyc6MPBoPxL46RU1TkQ0Qi0qeNmOntF5JhCQSURBLmupg0YmUPHJjkSERE0kdJJAF7X3wRmhMbCtq1ZInWh4hIn6EkkoDG2trEb25uZus//ZMSiYj0CUoiCYhXg747mj//nE/mz1ciEZFeT0kkAdlJGBwP19MSEenNEp7iezwKlzEZcPrp7KuogJaWhJ+lhYYi0hcoiXRRXXk5H8+bh2tqwnJyyDvpJI5++mnnN8ZheXmqmSUifYK6s7poz4sv4pqaAHBNTTTV1SX8LHf0uF2nKSJ9jJJIF8UOpmfFLjbspj0vvtij+0VEMoGSSBcNv/HG1lXqlpvL8Btu6NHzGrZvT0ZYIiJppTGRLhoyfTrjf/rTqMKJny1dijtyJKHntRw8mOQIRURST0kkAQfXr6e+ooLsQYNoSjSJRJaNFxHppZREuqjNplM95I4c4f2rruKc115LyvNERNJBYyJd1GbTqSRo2LaNzfPmJfWZIiKppCTSRYPLylp3JEymfe++2/lFIiIZSt1ZXXRw/XqyCgoIHj0KXdiXvqsKp0xJ2rNEjlf79++nurqaptBaLum6nJwchg0bxqBBgxK6X0mkC3Y++ii7Fi1K+nOtoIAJjx5XW9OLJN3+/fvZu3cvxcXFFBQUYNbTEqnHD+ccR44cYdeuXQAJJRJ1Z3VBze9/78tzXUODKvmK9FB1dTXFxcX069dPCaSbzIx+/fpRXFxMdXV1Qs9QEumC7BNOiHqd1b8/ZCXhn845b4MrEUlYU1MTBQUF6Q6jVysoKEi4K1BJpAtO/va3IRDwXgQCjHv4YYq/9S1Iwm899RUVao2I9JBaID3Tk38/jYl0wZDp05nws59FrVYHvNZID8rBA9DSQn1FhSr6ikivpCTSRUOmT4/6Qf/pE0/0PIEABALaV0REei11ZyWgrrycwxs3JuVZUS0bETmu3X///RQXF5OVlcWsWbPSHU6XpD2JmNkVZrbZzLaY2V0dXDfJzFrM7PqIY7eb2QYz+9DMXjCz5K8GjKO+ogKCwaQ8a/iNNyblOSLSu1VWVnLvvfcyd+5cKioqWLBgQbpD6pK0JhEzCwCPA/8LOAO42czOaOe6fwXeiDhWDHwPKHXOnQkEgJtSEXcyV68fXL8+Kc8Rkd5t06ZNAHznO9/h/PPPZ+zYsWmOqGvS3RKZDGxxzm1zzjUCy4Br4lz3XeC3QOxE5mygwMyygX7AZ34GGzZk+nTGPfII/b/4xR4/a/ezz2p2lshxbtasWdx6660AFBYWYma8/fbb1NbWMmfOHEaOHEl+fj4TJkzg0YgFyi0tLTz44IOMHz+evLw8Ro0alfJusHQPrBcDkRuVVwHnRV4QanFcB0wHJoWPO+d2mdmPgZ3AEeBN59yb8d7EzGYDswFOOeWUpAV/KAnjIsHDh/n49tsZ/9OfamxE5Di1YMECTj75ZP75n/+Z8vJyCgoKOP3005k6dSrV1dXce++9nHbaaWzZsoUtW7a03jdnzhyWLl3KnXfeycUXX0xdXR3Lly9PaezpTiLxJifHFqZ6FPiBc64lci6zmZ2A12oZA9QDvzGzrznnnmvzQOcWA4sBSktLk1L4KpnjIq6xUdN8RTLFnhVQ8yYMnQEjZqbkLceOHdvafTVp0iQGDBjAokWL2LBhA+vWrePss88GYHrEz4hNmzaxZMkSHnvsMb73ve+1Hr8xxeOs6U4iVcDJEa9H0bZLqhRYFkogJwJXmlkzkANsd87VAJjZS8AFQJsk4ofsgQOT9izLztY0X5FMsGcFvH8ztByGqqfhnBdSlkhilZeXc84557QmkFhvvfUWQNpncaV7TGQNMM7MxphZLt7A+IrIC5xzY5xzo51zo4HlwLedcy/jdWNNMbN+5mWYS4DkzLvtguYDB5L2rEEXXKBWyHFg5dZV/Lj8p6zcuirdoUh7at70Egh4f9bE7SFPidraWkaOHNnh+f79+ydcfTdZ0ppEnHPNwFy8WVcbgV875zaY2W1mdlsn967GSyrrgD/j/V0W+xxyq8FlZVhubs8fFAgwQtN8+7yVW1ex8PX7+O0HL7Hw9fuUSDLV0BkQ6Od9H+jnvU6ToqIidu/e3eH5Q4cOsX///hRG1Va6WyI45153zo13zo11zj0QOvakc+7JONfOcs4tj3h9r3PuNOfcmc65W51zR1MZe3hfEcvJOVZbq5uKv/lNtUIySLzWQldbEB1dt3rHGhqavY9nQ/NRVu9Yk9zAJTlGzPS6sEq+k9auLIBLLrmE999/n/XtLAMIj48sXbo0lWG1ke4xkV6rvqICF6p66RKsfjnk8ss5RdvjZoxwa6Gh+SivbXiN+6+8D6DNsQvHTo177z2vLqAp2MyKP6/gltKbOXj0EOeVTOLCsVM5r2QSr214jYbmo+Rn53FeyaQ2z5AMMWJmWpNH2Ne//nUef/xxZsyYwX333ceECRPYvn07H3/8MQ899BATJkxg9uzZzJ8/n+rqai666CLq6+tZvnw5y5YtS1mcSiLdUFde3lqEcXBZGTUvvdSjfde10DCztNdaiD0WL4n8bv0rNAWbAWgKNvPse78iSDAq8dx/5X2s3rGmNbGIdCQ/P5/y8nLuuusuFi5cyP79+xk9ejTf/va3W6954oknKCkp4amnnuKhhx5i2LBhXHbZZSmNM+3dWb1FXXk5n8yfz97nn+fj229nz4svMvjiiwkUFra9uIt7jTTu3auFhhnkvJJJ5GfnAZCfnceAvP7s3vcZOVnZrcfaa0HEzlUP4k3/VteVdNWsWbNwzjFgwIDWY0VFRfzyl7+kurqahoYGNm3aFDWdNxAIcPfdd7Nt2zYaGxupqqri6aefTmncaol0UX1FRWurwzU2su+Pf2z32v6nn07T55/T+FknC+iDQa0PySDh1sLv1r/C54freH7tMppamghkBRicX8jwQcMBr+tq9Y41DMjrz8GjhxiQ158dnx9bM5tlWQTdsTVEA/L6x+0qU2tE+gIlkS4aXFZG9W9+0+n4h+Xmcvjjj73rAoFOy8Unc72JJMf7Ve+3dmEBtARbqG/YR33DPn6w4m6yswKtXVfx9Mvtx8GjB1tff1L9Cau2VnSpW0ykt1ES6aIh06dTeP751MdpgWTl5zPyG9+g+cABGqqqjrVSurDfyMEklZSX5IgcF4nH4TpMIEBUAskJ5LBmZ2XUPRpYl75ESaQbht94I/vfe8/r1goEGHLppeQUFbXuCVJXXs7eF1/EcnJwTU3eny0tHZZH0aaeqbdy6ypeXv8KDrjurGvYsOcj3tj4JgePHuJg48FO7++OLIyjwWOt1xEDh3Nq0ZikvodIOimJdEO4em+bbXKBnY8+yq6nnoKWFiw3l8KLLmpdRPjxHXfgjrb97dZyc7WfiM/C4xfhGVErt67intcW0tTi/WBf/ZfVBNuUa0ueoy2Nrd8HLIvaw3XsObCX1TvW8LVJt3Bb2Wzf3lskFZREuil2m1zwZm6FEwh4A+/5o0a1XnfSrFnsWrQo6h7Ly2PQeVEFiyXJnqxYzHNrfkWLOzbVdvWONa0JBPA1gcRqcUFaWoKh71t4bs2v+OKIMzQ2Ir2apvj2UF15OTsffTR6/CNm3/RT5s1jwDnnRN/Y0sK+P/6RT+bP1zRfH6zcuopnQwkEvMHsRRW/ZEBefwJZiVUXSLYWF9T0X+n11BLpgfDakagFh4FAVCmT8ALF4m9+k5oVK9j37rtkDxrE0U+9KaHBhgZN8/XBy+tfiZpmC7C1dhs76z/FudS1PiIZYBHTf3OysjXALr2ekkgPRK4dASgYN45T5s2LSiDhJLM3XIYgGKTl4EEsNxfX2EhWfr7KwPugvTQR2ZWVauOHTWBz9ebW15NOKVVXlvR66s7qgci91rPy86MSCMQkmWDw2CytlhYKxo1j+C23MO6RR9QK8cF1Z13TutI8U5xUODJqRfy1Z8XbCVqkd8ms/8t6mY5ma4GXZPa++GLc9SK5RUWcumBBqkI97lw4dioPXP0jfrf+Fd7buYaWYOdrdvz29if/zbRxFzOk3xDVz5I2Zs2axYcffkhlZWW6Q+kWtUR6aMj06Zy6YEHc1sSQ6dMZcumlbY5raq//Vm5dxb+/+zTbarczYuDwdIcDeAsV3/7kv5VApE9RS8RnOUVFUa9zR46k37hxaYrm+LBy6yp++Oo/ZUTrI5bDseTd/1A1X+kz1BLxWeS4ieXk0FRbS72m9vpq9Y7M6L5qz8fVH2uHQ2nXf/3Xf3HWWWfRv39/pk6dyoYNG1rPtbS08OCDDzJ+/Hjy8vIYNWpU1B7r06ZN4/rrr2fx4sWMHj2agoICrrrqKnbt2uVbvEoiPguPmwy/5RYKzz8f1+itYA5P7ZXky/Rpsy40d6yh+Si/W/9Kmo4uysoAABJ4SURBVKORTLJz507+8R//kXvuuYcXXniB6upqbrjhhtZp6XPmzOHee+/lhhtu4NVXX+WRRx7h0KFDUc/4n//5H37+85/zk5/8hCVLlrB+/XquvfZa32JWd1YKhFe515WXt9be0tRe/1w4dirfmHwrz7z3bLpDacOArKxAa0up8tO1rNy6St1aGSS2VE4q1dXVUVFRwbhQl3cwGOS6665j82ZvaviSJUt47LHHovYUuTFmfLW6upp33nmHkpISAEpKSpg6dSp/+MMfuOKKK5Ies1oiKRTZKtHUXn/dVjabkhNOSXcYbWRlBRiYe6z8f1NLk1atZ5Dwvi/p6m4cPXp0awIBOOOMMwCoqqrirbfeAojqvopn4sSJrQkEoKysjGHDhvHee+8lP2DUEkm5eLW3JPlWbl3FzoiNojKFtzdJfetrlYXPLPG2SE5la2Tw4MFRr3Nzc71YGhqora2lf//+DBo0qMNnDBs2LO6x3bt3Jy/QCEoi0mes3LqqdVfC7XV/aR17yFRji05lTtm31JWVQc4rmcRrG16jofloxiX4oqIiDh06xP79+ztMJNXV1XGPjRw50pe4lESkT1i5dRX3vLqg0w2j0i0QGg/Jz85TAslA4S2SM3EK9vRQD8bSpUuZO3duu9etW7eOnTt3csopXnduRUUF1dXVTJ482Ze4lESkT1i9Y03GJxCAyadM4qTCkzLuB5Qcc+HYqRn532bChAnMnj2b+fPnU11dzUUXXUR9fT3Lly9nWbg2H17X1dVXX819991HQ0MDP/jBD5g4caIvg+qgJCJ9xIC8/ukOoUsMlEAkYU888QQlJSU89dRTPPTQQwwbNozLLrss6przzz+fSy+9lHnz5lFTU8O0adNYvHixbzFZuspip0tpaanrbbVppGPhGTUd7Y2eCbJCZeDzs/O4/8r7lEiSZOPGjZx++unpDiMjTJs2jRNPPJHly5d3+96O/h3NbK1zrjTeOU3xlV4vckZNJsrJyqG48KTWfUTCs35E+gIlEen1ziuZRMAy96PcFGxi9/49rTsq5mfnMSCvPz8u/6nKnkivpzER6fUuHDuVM0acwZ93f5juUNoVdEHMGReMnsK4YeN4cd2vaWg+2rr3u7q2JBnefvvtlL9n5v76JtINNQdr0h1CpxyOkYUncfDooTYL2kR6q7QnETO7wsw2m9kWM7urg+smmVmLmV0fcWywmS03s01mttHMzk9N1JJpwl1FmSy8eO28kklROxxm0oI2ke5Ka3eWmQWAx4HLgCpgjZmtcM59FOe6fwXeiHnEY8AfnHPXm1ku0C8FYUuGWbl1FXsO7E13GHFNGDaBKaMnc/DooaipvZm6oE2ku9I9JjIZ2OKc2wZgZsuAa4CPYq77LvBboPVXNjMbBFwEzAJwzjUCjf6HLJkmU/cPyc/Op6jfCXxxxBltEkWmLmgT6a50d2cVA5FV8qpCx1qZWTFwHfBkzL2nAjXA02b2vpk9ZWZxV5yZ2WwzqzSzypqazO87l65ZuXUVPy7/KXWH69IdSlwNzQ2885d3uee1hZqFJX1WupOIxTkWu/rxUeAHzrnYXzWzgYnAvznnzgEOAXHHVJxzi51zpc650qFDh/Y0ZkmxcLKI/EEcWbL7vz/57zRG1zmVe5e+LN1JpAo4OeL1KOCzmGtKgWVm9hfgeuAJM7s2dG+Vc2516LrleElF+pD29neIXGAYzPBqvYGsgAbPpVOzZs2itDTuovCMlu4ksgYYZ2ZjQgPjNwErIi9wzo1xzo12zo3GSxTfds697JzbA3xqZhNCl15C27EU6eXi7e8ARM1wClgWFrdRmz5fGnkmWaH/vbIyLDaRZEprEnHONQNz8WZdbQR+7ZzbYGa3mdltXXjEd4Ffmdl64GzgX/yLVtKho+mwE0edw4Rh48nKCuBwGJbWlevFhcWMLTqVb0y+lfHDxhPEK3PSFGxWd5b0WeluieCce905N945N9Y590Do2JPOudiBdJxzs5xzyyNe/yk01nGWc+5a59znqYxd/Bfe3+Hvvvy3rSu7w11c7/zlXbbUbKWppQnwFvONG/oFxhadSskJp1CQXUB+dn5K4jSMf7h4Ls99/RluK5uttSCSsJdffpnTTjuN/Px8pk6dykcfHetgaWlp4cEHH2T8+PHk5eUxatSoNtvl/u53v2Py5MkUFBRQVFTElVdeyY4dO3yLN91TfEU6FTsdNrKLq8W1ELAsWlyQnKxsttZub00qqZRlx7qsVm5dxeoda7hx4g1t1oeIdGTHjh3ccccd/OhHP6KgoIB7772Xyy+/nE8++YT8/HzmzJnD0qVLufPOO7n44oupq6uLqtj77LPP8vWvf52bbrqJBQsW4JyjvLycmpqaqH3Xk0lJRHqd2C1Mwz+sd+/7jHf+8m5aYmpxwdYuq3BZepV8753qysupr6hgcFkZQ0K7CabKX//6V1555RUuuOACAM4991zGjh3Lf/zHfzBt2jSWLFnCY489xve+973We2688UYAgsEgd911F9dddx0vvPBC6/mZM2f6GrOSiPQ68bYwXbl1FS+vf4WcrOyU7nAYud3teSWT4k4EUBLpPerKy/lk/nyCDQ3UvPQS4x55JKWJZNiwYa0JBKCkpIRzzz2X9957j/DeT7HdV2GbN2/ms88+4+///u9TEWorJRHplSK7uCI3pcoJ5DBi4PCklUHJIqt1gDyegbkDyM/J5/LTZ7TGE9lK0lhI71JfUUGwoQGAYEMD9RUVKU8i8Y7t3r2b2tpa+vfvz6BBg+LeW1tbC8DIkSN9jTFW2gfWRXpi5dZVLKr4Zetv/00tTZxaNKZ1UDtSbiC328+PTSAjBg4nJ5DT+rq+YR97Duzl+coXWLl1VdyJANJ7DC4rIyvfm4yRlZ/P4LKylL5/dXV13GMjR46kqKiIQ4cOsX///rj3FhUVAbB7925fY4ylJCK91pMVi/nh7+9ha+221mP52Xlce9Y13H/lfVwwego5WV5jOxDamrYn8rPzuONv5vHAVfcztujUqHOR03gvHDuV70+/XQmkFxoyfTrjHnmE4bfckvKuLPASxjvvvNP6eufOnaxbt47JkyczPRTL0qVL4947YcIEiouLeeaZZ1ISa5i6s6RXWrl1Fc+teZ6WiMQwYuBw7vibea0/vMNjJUvefZrN1R9DD5JIYUEh135pZlRiuOe1ha0zwXKystV11UcMmT495ckj7MQTT+TWW29tnZ21cOFChg0bxqxZs8jPz2f27NnMnz+f6upqLrroIurr61m+fDnLli0jKyuLhx9+mK9+9at89atf5eabb8bMKC8v5+abb/ZtNbySiPRKq3esoSWmnFptnEKMF46dyqKKX/b4/fYd2cfzlS9EVeQtPflc6g5/zpB+J3DdWdeo5SE9VlJSwt13381dd93Fjh07KC0t5YUXXiA/1MX2xBNPUFJSwlNPPcVDDz3EsGHDuOyyy1rvv+WWW8jPz+eBBx7g+uuvp3///kyZMgU/awZaeMT/eFFaWuoqKyvTHYb0UORgeqS/+/Lf8v3pt0cde7JiMc+892zr67zsXI42x981oLjwJC6dcAnPr10Wd73J+aOncN1Z12gabwbZuHEjp59+errD6PU6+nc0s7XOubhNGY2JSK8UHsCOHPdobzbUbWWz+cbkW1tLktw08cZ2nztl9BRuK5vNA1fdzwWjp1BYUBh1flvtdl5e/4q2txUJUXeW9Frhab7hFeLtrQxfuXUVB48eYk7Zt6LOv7HxTfYeqMZFVAEekNe/zbMjxz72HthL3aHa1vUhkWMhncUh0hepJSK9XkezodorJX9b2Wx+983ljB82Lur6j6s/afPs2NlYTcHm1oVfhMqdtPc+In2dkoj0ae2Vkg8r6jck6nW8ou0Xjp3KnLJvRZWeD08XDm841dn7iPRVSiLSp3VWTffas65pXTyYE8jh2rOuifucyEWEF33hoqhzA/L6q2pvmh1vE4SSrSf/fhoTkT4tXp2t2PMPXHV/l8YywuMkPy7/adTxg0cPdfo+4p+cnByOHDlCv3790h1Kr3XkyBFycnI6vzAOTfEV6abI6cWa4pt++/fvZ+/evRQXF1NQUICZdpLsKuccR44cYdeuXQwfPrzdulwdTfFVS0Skm9TqyCzhH3yfffYZTU2p30umt8vJyekwgXRGSUQkAbEbZUl6DRo0KOEfgtIzGlgXEZGEKYmIiEjClERERCRhSiIiIpIwJREREUmYkoiIiCTsuFtsaGY1wI40vf2JwF/T9N6dydTYMjUuyNzYFFf3ZWpsmRJXiXMu7s5Wx10SSSczq2xv1We6ZWpsmRoXZG5siqv7MjW2TI0rkrqzREQkYUoiIiKSMCWR1Fqc7gA6kKmxZWpckLmxKa7uy9TYMjWuVhoTERGRhKklIiIiCVMSERGRhCmJJImZXWFmm81si5ndFee8mdnPQufXm9nEmPMBM3vfzF7NlLjMbLCZLTezTWa20czOz6DYbjezDWb2oZm9YGb5KYzrNDP7HzM7ambf78696YjLzE42s7dC/w03mNk/JDOunsQWcT5dn/+O/lum+/PfUWy+ff67zTmnrx5+AQFgK3AqkAt8AJwRc82VwH8CBkwBVsecvwN4Hng1U+ICngG+Gfo+FxicCbEBxcB2oCD0+tfArBTGNQyYBDwAfL8796YprpHAxND3A4GPkxVXT2PLgM9/u3FlwOe/vf+evn3+E/lSSyQ5JgNbnHPbnHONwDLgmphrrgGWOs+7wGAzGwlgZqOAq4CnMiUuMxsEXAQsAXDONTrn6jMhttC5bKDAzLKBfsBnqYrLOVftnFsDxG6j15W/U8rjcs7tds6tC31/ANiI94MoWXryb5bWz397cWXC57+jfzP8+/x3m5JIchQDn0a8rqLt/6QdXfMocCcQzKC4TgVqgKdD3QxPmVn/TIjNObcL+DGwE9gN7HPOvZnCuPy4NyXPNrPRwDnA6qRE5elpbOn8/LcnEz7/cfn8+e82JZHksDjHYudOx73GzK4Gqp1za5MfVuJx4f2mMxH4N+fcOcAhIJl9/D35NzsB77e2McBJQH8z+1oK4/LjXt+fbWYDgN8C85xz+5MSVejRcY51KbYM+Py3JxM+//Fv9Pfz321KIslRBZwc8XoUbZuX7V1TBsw0s7/gNWmnm9lzGRBXFVDlnAv/xroc73+qZOlJbJcC251zNc65JuAl4IIUxuXHvb4+28xy8BLIr5xzLyUppmTElu7Pf0f3pvvz3x4/P//dpiSSHGuAcWY2xsxygZuAFTHXrAC+HppxNAWvCbrbOfdD59wo59zo0H3lzrlk/VbRk7j2AJ+a2YTQdZcAHyUprh7FhteMn2Jm/czMQrFtTGFcftzr27ND/0ZLgI3OuZ8kKZ6kxJYBn//24sqEz397/Pz8d1+6RvT72hfeTKKP8WZc3BM6dhtwW+h7Ax4Pnf8zUBrnGdNI4uyUnsYFnA1UAuuBl4ETMii2/wNsAj4EngXyUhjXCLzfJPcD9aHvB7V3b7rjAqbidZWsB/4U+royE2LLgM9/R/8t0/357yg23z7/3f1S2RMREUmYurNERCRhSiIiIpIwJREREUmYkoiIiCRMSURERBKmJCIiIglTEhERkYQpiYikkJllm9n3zOwDM2swsz1m9ovQ6uN9ZpbMVdEivstOdwAix4tQeYvfAzPwVkL/DDgR+N94VWMHAUndlEnEb0oiIqnzC7wE8o/OuR+HD5rZM8DboZfr0hCXSMJU9kQkBcxsEvAe8KZz7vI458O73F3inCtPdXwiidKYiEhqzA39eX8752tDf74fe8LM7jYzZ2a/8CUykR5QEhFJjcuBWudcRTvni/H2iPg88mCoBP638CrJimQcJRERn5lZPjAcbx+IeOfPxNuh7v2Y44XAr4D/B/g8zq0iaackIuK/ltBXUTvnF4b+jB1UXwws1xiJZDIlERGfOW8L00+AU8zsb8LHQzs2LgS+Ejr0fsS5bwFfABakMlaR7tIUX5HUeBj4d+A1M3sBqMPbK3sg3rarZxBqiYS2ZP0X4ELnXGN6whXpGk3xFUkRM7sD+C7e+MdnwHLgAbxWSpNz7qTQdbOAp/G6wMICeFvcBoH+zrmjqYtcpH1KIiJpZGYn4w24v+6cuyp0bDAwKubSp/GSzb8AG5z+x5UMoe4skfQ6J/Rn66C6c64eqI+8yMwOAXXOuQ9TGJtIpzSwLpJe4STSZpGhSG+glohIerVpicTjnJvmfygi3acxERERSZi6s0REJGFKIiIikjAlERERSZiSiIiIJExJREREEqYkIiIiCVMSERGRhP3/jhBmuhMikrQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(fccaq[0], fccaq[1], s=10, label='fcc', color='#FFB300')\n",
    "plt.scatter(hcpaq[0], hcpaq[1], s=10, label='hcp', color='#388E3C')\n",
    "plt.scatter(bccaq[0], bccaq[1], s=10, label='bcc', color='#C62828')\n",
    "plt.xlabel(\"$q_4$\", fontsize=20)\n",
    "plt.ylabel(\"$q_6$\", fontsize=20)\n",
    "plt.legend(loc=4, fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks much better! We can see that the resolution is much better than the non averaged versions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is also the possibility to calculate structures using Voronoi based neighbor identification too. Let's try that now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcc.find_neighbors(method='voronoi')\n",
    "fcc.find_neighbors(method='voronoi')\n",
    "hcp.find_neighbors(method='voronoi')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcc.calculate_q([4,6], averaged=True)\n",
    "fcc.calculate_q([4,6], averaged=True)\n",
    "hcp.calculate_q([4,6], averaged=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "bccaq = bcc.get_qvals([4, 6], averaged=True)\n",
    "fccaq = fcc.get_qvals([4, 6], averaged=True)\n",
    "hcpaq = hcp.get_qvals([4, 6], averaged=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the calculated points.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f5cfbe5a510>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEPCAYAAACDTflkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3de3iU9bnv//edSUggnAQJIGhQinhYCysG1AaVldZD0WJt3eKhtnStFlitddFi1erPw8/War1qN65dXUq1KrUVLVWk1VXt3qlbjYpEbFEEBLRgEEwkciaQzHz3H8/MMJlMkslknnkm4fO6rlyZeQ7z3NImd76n+2vOOURERDJREHQAIiLScymJiIhIxpREREQkY0oiIiKSMSURERHJWGHQAeTa4Ycf7saMGRN0GCIiPcabb775iXNuWKpzh1wSGTNmDLW1tUGHISLSY5jZxvbOqTtLREQypiQiIiIZUxIREZGMKYmIiEjGlERERCRjSiIiIpKxQ26Kr4jIoaKxuprtNTUMrqxkSFWVL89QS0REpBdqrK5m3bx5fPy737Fu3jwaq6t9eY6SiIhIL7S9poZIUxMAkaYmttfU+PIcJRERkV5ocGUlBSUlABSUlDC4stKX52hMRESkFxpSVcW4u+/WmIiISLu2LoW3r/K+H6Iaq6t5/8c/bjPmsWn+fDbNn0/hgAG+JRBQS0REeqqtS+GtyyC8F+oehpMfhxHTg44qJ2KzrgoHDGDLo48SaWqi4amnGPmNb9CyaxfN27bR+PzzAGxetw6Ao+bO9SUWJRER6ZkaXvASCHjfG144JJJIbNZVpKkJCgogEgG8wfPNv/pV/H2ihj/+0bckou4sEemZhp0DoX7e61A/730vlNxdlTjrqk3CSJFAAA7U1/s2xVctERHJb1uXeq2MYee0bmmMmO51YaU610sktjo+fuIJ+k+YQNPGdrf2aF9LC9tranwZG1ESEZH81dm4x4jpvTZ5bK+pYX9d3cFWRzjM7rfeyvgzCwcMyFJ0SZ/ry6eKiGRD8rjHO9+DdbdBn+FQPrvXJJDE8iQA733/+7gDB7L6jN2rV2f182KUREQkP21dCns/gII+EIn+Qm3a5H0BbPvfMPH3PT6RNFZX897cubjmZup//3v6HXts1hMIgGX9Ez1KIiKSfxK7sQqKoeSog8kjJnIA1tzove7BiWTrE0/gmpsBcM3N7I1Oyc2qwkKGz5iR/c9FSURE8lFiN1ZkPxQPa5tEAHa/A7UXwciveF1cPXCAvWXbtlbv/WiF9Ckry/pnxmiKr4jkn+Tpuy7pfEFJwpsIbFkMG+/1Wi89YPX6pvnz+dv06aydO5c9a9b4/rwDH33Ee3Pn+jLNVy0REck/idN3D3wMW/7Q+nzRYNi/te19PWDR4duXXx6fZbXPj66rdrjmZrY+8UTWp/kqiYhI/qr/M+zb0Pb4/vrU1+f5osO1c+d2a5pud/kxuK4kIiL5Z+tSb6yD1CuwWx0vKYfRV0DzjrweE2msro7XswpK6fHHZ/0zlUREJL9sXQpvXUn7CSTJ6CvguNt9DSkb/NoUqiu2PPoo/SdMyGqXlgbWRSR/rLnRa4GEd6Z/T/MO/+LJosGVlVhRUaAx+LHDYeBJxMzOM7O1ZrbezK5PcX6qme0ws79Fv25OOh8ys7fM7E+5i1pEsm7rUlj/M9JugcQc+NiXcLJtSFUVg04/PdggQqGs73AYaHeWmYWAe4GzgTpguZktdc69m3Tpy865C9r5mP8AVgMD/YtURHwTK7D46etAuOv3b1nstWB6QJdWQd++AT68gFHf+lavm501GVjvnHsfwMwWARcCyUkkJTMbDZwP3A78wK8gRcQniSvTu/s5eZ5ENs2fn/uBdTNKTziBwqFDGTFjRq+s4jsK+DDhfR1waorrTjezvwMfAdc451ZFj88HrgU6LE9pZrOAWQBHHXVUd2MWkWxJXJneHXk8IytWWNGv/TzaFQrRp6yMwVOm+LYhFQSfRFJNW05em7oCKHfO7TazacASYJyZXQDUO+feNLOpHT3EObcAWABQUVGR/PkiEpSiQXi/BrrxY9l3bF62QhL3Aml46in6+TC9tkPhMAe2bGHzAw8A/m2PG/TAeh1wZML70XitjTjn3E7n3O7o6+eAIjM7HKgEppvZP4BFQJWZPZaTqEWk+7YuhQ130a0EQgG4loOFGPNI4g6EkaYmDmxNscI+R/xsBQWdRJbjtSqONrM+wKVAq8I3ZjbCzCz6ejJezNuccz9yzo12zo2J3lftnPtabsMXkYxtfMBLAN0SgaaNsP6neZdIBldWUlDi1fgqKCmhaMiQwGLxYywkJtDuLOdci5ldBTwPhIBfO+dWmdmc6Pn7gYuBfzezFmAfcKlzTl1SIj3Z1qXefiDZ/sw86tYaUlXFuLvvZntNDYUDBvDJc88FEke/E0/s1WMisS6q55KO3Z/w+pfALzv5jBeBF30IT0SybetSeOfqgxtNZUueDa7HBtULBwxgy6OPHtzmNpfMOPI73/H1EYEnERE5BMTWghQNgvfv9vYIySqDwakmdgYjcVA9SP0/+1lfu7JASURE/NZqLUiIjBYUdsrlVQn4xEH1IO1euZLG6mpfE0nQA+si0tu1WgviRwKJKhrk32d30eDKSgiFgg4DwmFvkWMvnp0lIr1d4i6FBX3w7ddOHhViHFJVxahvfQsKgv8Vu2/dOtbNm+dbIgn+v1BEerfYLoXDpkFoEF0usJiuPGqJgLe4b9S3vx3Y8/uMHBl/7Uf13hglERHx3+bfQMNz0Nzg3zPevzvv9lffs3p17h9aWMio2bPpN25cvPR8QUlJ1qv3xh/ny6eKiMSsudGrtOu3yP68GlyH7q3Fz1S/8ePjU4qtTx8GnXmmb8UXQS0REfFbrloHBcV5s796Y3U17//4x/Q//nh/x0WSP7uwkD5Dh8ZnhrkDBygZPbr3rlgXkUPAgONg9zv+P2f4l/KiFZK4RsT69KFwyBBaPvnEn4dFIlhhIa6lxdsv5N/+jf4TJrDzjTeINDX52o0VoyQiIv7qMzw3z9m6xGv1BJxIEteIuAMH/EsgUa4lWn8sEqFl165W5VYGV1b6vthQ3Vki4q9h55CTXzWuxRsTCdjgykqsT5+cPzex1TGkqopjbrrJ9wQCSiIi4rcR06HIxwq2Fu1QCfXL6ZhIbNwj1fqLeOsgRwr69WPkN76Rk6TR5tk5f6KI9B5bl8LbV3U+eF4+y5/nF/SHU/4A5d/11qLkqCsrNu7x8e9+12YhX91990HEp7Uw7Yjs3cuWRx/N/e6JKImISKZiNbE23ut9T0wkycllz3v+xBDZDRt+5s9ndyB5w6nEhXzNn36a83hSxZErGlgXkcwk1sQK7z24RiOx4OLG/8Lb/tbHmlmfvup91T2cs9bI4MpKGp56Kj4Dq6mujsbqanavXBlo4cXCAQNy/ky1REQkM4k1sRLHI1oVXIzgawJJFEtkORCbATX4zDPBOXa89BJrr76azQ88QEtjY05iSKVl166cP1MtERHJTKwmVsMLXgKJtQCCqmGV44H1IVVVbK+pwTU3ewfCOUqW7cjFmpCUz835E0Wk9xgxHf75l627kIKoptv/n3I6sB6TuI96kKXfS088kXF3363ZWSLSCww7J1ryPVcKvL3VA1hkGOvWGn755Yz/z/+k/8kn5zwGgP4nnRRIAgF1Z4lItiRugWslQJb3UG/PZ64PdJX6kKqqVr/A1151FTifSi8WFkLSGpSgurFilEREpPtabYGbS/m1t/r2mpqOE0hxMezvxv7yCQnEiooYdPrpDPexQm861J0lIt238YEAEgiA856dJzrbFre4rCwrz+k7bhzHzp/P8Q88EGgCASUREemurUth2/8JOoq8EN8W1yzleQcUDe9iQcpQiNITT4zX4yooKeGouXMDTx4x6s4Ske5peMHbECqmoBQie3L08AIon52jZ6XnqLlz2b16NTteeqnNuQMfftilz+p34okc+Z3vMKSqisbq6pxV5u0KJRER6Z7kdSE5SyAAocBLv6cyYsYMdkX39OiOAQmzrpIH8POFurNEpHt2/C24Zxf2De7ZHUic+jvk3HPTvq/0xBPjYypBz7pKl1oiIpK5oMdDxlwV3LM7EWs5bJo/P63rC0pKGP2d7wBk3G0VRJdX4EnEzM4D7gFCwIPOuTuTzk8FngE+iB56yjl3m5kdCSwERuAV6FngnLsnZ4GLSNvxkFwqPc5bZJjHGqur2fzgg51e13fcuDaD5bGKvOkmg8RteRueeipnK9gDTSJmFgLuBc4G6oDlZrbUOfdu0qUvO+cuSDrWAsxzzq0wswHAm2b2lxT3iki2JS4sLOgDkdjCQsObg+SzwsHwL6v9f043ba+paVtTy6zVWpLk2VaZJoNU5el7fRIBJgPrnXPvA5jZIuBCoNNE4JzbAmyJvt5lZquBUencKyLdkLiwMNQP+v8z7HwzejIHCQS83QyXne+9Lp+dl4Pr4K0b+fh3v2t1rM/IkRw2dSqFAwbQsmtXm66nTJNBYnn6XI6nBJ1ERgGJc97qgFTLT083s78DHwHXOOdWJZ40szHAycAyf8IUkbjkfUSKh3vJJJZUCg+D/Zv9jaG5ERqei8bzZ6/0SR52bQ2pqmLIuefS+Pzz8WOFhx3W4ZhF8p4g6e4REhvMz/WYSNCzs1KtyEn+U2YFUO6cOwn4X8CSVh9g1h/4AzDXObcz5UPMZplZrZnVNjQ0ZCFskUNY8j4ihf2gaBgMPMWrpGu5qGYbaf16w12db9EbkPHz5zNq9mz6HHEEhELsXbWqzZa6cHDP9t2rW3fTdWWPkCFVVRxz0005nQocdBKpA45MeD8ar7UR55zb6ZzbHX39HFBkZocDmFkRXgL5rXPuqfYe4pxb4JyrcM5VDBs2LNv/DSKHltg+IuXfhbJpsGUxNG30urQ2/Aya6nIfk2vJ2YZUHYklguQEcdTcuRw2dWp8fCR5K9vEPdt3vv46VlQE9IxpvkEnkeXAODM72sz6AJcCrf6cMLMRZl4NATObjBfztuixh4DVzrlf5DhukUNL8p7psX1Edq1pfd2nr9K6leCTwz4H/ca2PhbUZlhRiYkgVUsjce+R5OSQOA7iDhzwCitefnlge4R0RaBjIs65FjO7Cngeb4rvr51zq8xsTvT8/cDFwL+bWQuwD7jUOefMbApwJfC2mcVWO90Qba2ISLYkDqR/+CsY+oWDg9lFA4OJKbwfhp0HG+89eCwHm2F1tA6jswHxjsYskgfFg67M2xXm/Kp7n6cqKipcbW1t0GGI9BxvX9X6lzV4YyFHz/XGIlxL6vt8FYLPXAcfzD84oO/zzoaJU28LSkratBI6O5/O5+djbSwAM3vTOVeR6lzQs7NEJN8NOwfqHm5d6j2812uhBJJAAMJeyyPVHu8+6U5LIx35WhurM0oiItKx2ED6xge8EieR/VBQ7H1ZYTCJJNTvYOLI0RqRdNZh9NRE0B3qzhKR9G1dGk0m/zu6Sr2AnAykxxV4YyEBLTBM7nLK5y6obFJ3lohkx4jpXvdRQ6zMSYTcJhIHu1bB9mWBJJHElkZQtaryTdBTfEWkp0lebNjv6Bw+3HlrUtb/FNbcmMPntpVqjORQpCQiIl2TuNjw5MehdHwwcdT9NpjnRnW07uNQou4sEem65AHthj+T27ERvPpcW5cGVnwxqFpV+UZJRES6Z8R0rwDi+p/m9rmxUicBVvA9FGdjJVN3loh03+Yncvcsi/7tG5vmK4FSS0REumfrUti3IXfPO/wcbzA/BwsMpXNKIiLSdbGdDYedk9vquaF+eb0J1aFI3Vki0jWxgowb7/W+Fw06OOXXT4NO8b0+lnSdWiIi0jXJOxsm1rAqGuS93/467Hiz48/pipEXwym/z97nSdYoiYhI1yQWZExVw2rrUvj0dbzdHcLZeWbL3s6vkUAoiYhI18QWG6aqnvvm//B2Osy2T14IdE2ItE9jIiLSNbEijHs/aH18zY3+JBDw1oRsfMCfz5ZuUUtERNK3dSmsuMQrBw9eNd+Jv/daCFuXdnyv9EpqiYhI+hpeOJhAwCsHH5viO+A4f5896LP+fr5kRElERNI37BxvM6qYgj4HV40nd29lWw72UJeuU3eWiKRvxHSY+OTB8Yny2d73ZefDjhX+PTcxWUleURIRka5Jns771mWt91/Ptj4jYMIDmpmVp9SdJSKZS1x46BcX9p6jgfu8pCQiIpkbds7BqropFQBF3XtGc8PBEitKJHlHSUREMjdiOoy9Fm91eioRoDk7zwrvzW2xR0mLxkREpHuOu937/o/7oaXRv+do/5C8pJaIiHTP1qXwwXx/E0hJuSr45iklERHpnlwMru/f7O/nS8a6nETM7AwzG5/w/jtmtsrMdpnZO2Y2p4ufd56ZrTWz9WZ2fYrzU81sh5n9Lfp1c7r3ikgODDvHh/1EksZYVDsrb2XSEvkv4AgAM/su8DNgCfDvwNPAndHjnTKzEHAv8EXgBOAyMzshxaUvO+c+G/26rYv3ioifYlV9h03rZKZWF4RKoWhYdj5LfJXJ/+Jjgfejr78FXOWcezR20sxWAj/G+wXfmcnAeufc+9F7FwEXAu/6fK+IZEtsq9zy2d7Xxgdg/8ew6+9eCyIT4Z3gir2k5Fq8Feux1fGSVzJpiewCDo++PgL4W9L5FcBRaX7WKODDhPd10WPJTjezv5vZf5vZiV28V0T8krxVLsCpz8KZtXB4N2dSRfZ7n1H+3YOVgiXvZJJEngOuir7+K3BJ0vkZwHtpfpalOOaS3q8Ayp1zJwH/C6/rLN17vQvNZplZrZnVNjQ0pBmaiHQqeavcxHUc5bO7MFaS4leRFXqf8c+/VALJY5kkkeuBM83sZby//r9vZi+b2QIz+7/ALcCP0vysOuDIhPejgY8SL3DO7XTO7Y6+fg4oMrPD07k34TMWOOcqnHMVw4apn1UkaxIH1ZPXcSSOlRT08Y4V9Gln3CSS9N68RYxKHnmvy2MizrmtZjYRuA5vDMLwxieOBGqASudcbZoftxwYZ2ZHA5uBS4HLEy8wsxHAx845Z2aT8RLfNmB7Z/eKiM862io3dj62YVXDC1A0COqfh32boGggNNV53VYFxYDz9ichBJ+57uAiRslrGU2lcM7tAG6IfmXMOddiZlcBz+PN6fu1c25VbJqwc+5+4GLg382sBdgHXOqcc0DKe7sTj4hkILGqb6JY4khMLm9+9eBge3gnHDPP2yck1oJpLxlJ3jLv93EHF5i9A7yFNzbxFvBWNIn0SBUVFa62Nt2GkohkJLFEfKif11rZ+AA0PNf6uvLvemMektfM7E3nXEWqc+m0RI7EW4cR7yoysw9onVhWOOfqsxCriPQGHQ24x5m3G+LWpa33J1FrpEdJJ4mMABbjLep7H288YhTwVeArsYvMbAttE8umbAcsIj3AsHOg7uGDLZFYd9W2/xPdo70ArMBrmTS+6LVU4GDrpe5h1crqIdJJIncBZwFnOedejh00s38CbsYbs3gXKAPOj36BN91CVYJFDkXtDbhPfNI7tveDg11biS2V5NaLkkjeS+eX/MXAosQEAuCcewe4JFrL6mvAZ4CBwMnAROCzWY5VRHqSVAPuibO1Gl9s21JJ1XqRvJZOEhkMfNLeSefcbWb2ZeAG59z1eOs3/pil+ESkN2qvpdLRdGHJS+nMzvobsN85d2oH19wFfNU5NzbL8WWdZmeJiHRNR7Oz0lmx/mtgkpl9v4NrhhCt7CsiIoeOdJLIvcBfgJ+b2XNmdlriSTM7C7iMdkqOiIhI79XpmIhzLmxmF+Alk28B55rZNmATMBSvYq/hFUcUEZFDSFoFGJ1zzc65WcCpwGNAGG8G1mhgLfCvzrn5vkUpIiJ5qUvrOJxzy4FvAJhZMRB2LtNdZ0REOqDV6z1CxosBnXP7sxmIiEhcYu0trV7Pa5nsJyIi4q+0am9JPlASEZH809FmV5JXVNtKRPJPZ5tdSd5QEhGR/NTeZleSV9SdJSIiGVMSERGRjCmJiIhIxpREREQkY0oiIiKSMSURERHJmJKIiIhkTElEREQypiQiIiIZUxIREZGMKYmIiEjGlERERCRjgRdgNLPzgHuAEPCgc+7Odq6bBLwOzHDOLY4e+z7evu8OeBv4pnOuKSeBi0je2LlzJ/X19TQ3NwcdSo9TVFREWVkZAwcOzOj+QJOImYWAe4GzgTpguZktdc69m+K6nwHPJxwbBVwNnOCc22dmTwKXAo/kKHwRyQM7d+7k448/ZtSoUfTt2xczCzqkHsM5x759+9i8eTNARokk6O6sycB659z7zrkDwCLgwhTXfQ/4A1CfdLwQ6GtmhUA/4CM/gxWR/FNfX8+oUaPo16+fEkgXmRn9+vVj1KhR1Ncn/3pNT9BJZBTwYcL7uuixuGiL4yLg/sTjzrnNwM+BTcAWYIdzLuUemmY2y8xqzay2oaEhi+FLspc3vMLPq/8nL294JehQ5BDR3NxM3759gw6jR+vbt2/GXYFBJ5FUfza4pPfzgeucc+FWN5odhtdqORo4Aig1s6+leohzboFzrsI5VzFs2LAshH3oSSc5vLzhFW5+7lb+8PenuPm5W5VIJGfUAume7vz7BT2wXgccmfB+NG27pCqARdH/yMOBaWbWAhQBHzjnGgDM7Cngc8Bjfgd9qIklh6aW/Ty76llum3YrZ4yd0ua6ZRuX09SyH4Cmlv0s27g85XUi0nsE3RJZDowzs6PNrA/ewPjSxAucc0c758Y458YAi4HvOOeW4HVjnWZm/czLMJ8HVuc2/N4rseWRKjmkcmr5JEoKiwEoKSymf3GpurZEerlAk4hzrgW4Cm/W1WrgSefcKjObY2ZzOrl3GV5SWYE3vbcAWOBzyIeExG6pH/3x/6Nxb2Or5HBq+aSU950xdgq3TbuVr570FWZMvIQnVjypri2RLrjtttsYNWoUBQUFzJw5M+hw0hJ0dxbOueeA55KO3d/OtTOT3t8C3OJbcL1ArCVxavmktLuWElseYRfmpfUv8bVJV7B7/542nxP7/P7FpbxXvw4DvjzhQnVtiXRRbW0tt9xyCz/96U+ZOnUqZWVlQYeUlsCTiPgn3bGMZKeWT+Lpvy8hQgSAsIvwXv06jhh0RLufn2j5h28y5ZhKQhYi7MIdtl5ExLNmzRoAvvvd72a88C8IQY+JiE9e3vAKD9T8Kq2xjJSSJmu8sWl5vHvrhj/dxM+r/ydLVj7TJoEANIebeWn9S4RdmJAVMGPiJWqFiHRg5syZXHnllQAMGjQIM+PFF19k27ZtzJ49m5EjR1JSUsL48eOZP39+/L5wOMwdd9zBscceS3FxMaNHj855N5haIr1QqhZCUUFhq9bAyxte4emVz8S7nxJ/yS9Z+QwRF2n1meGIN8M67ML8dd2L8c8sChXRHG47vzzsDrZidu/fk63/NJFe6aabbuLII4/kJz/5CdXV1fTt25fjjz+eKVOmUF9fzy233MJxxx3H+vXrWb9+ffy+2bNns3DhQq699lrOOussGhsbWbx4cU5jVxLphRLHI+Ki88BjyWPZxjfiiWLZpuV85vDPMKTfYZQUlfDqP15P6znNkRZGDTqCXU27OazfYDZ9+iEuaZmPurKkx9m6FBpegGHnwIjpOXnk2LFjGTt2LACTJk2if//+PPDAA6xatYoVK1bw2c9+FoCqqqr4PWvWrOGhhx7innvu4eqrr44fnzFjRk5ijlES6YVOLZ/E0reX0hxpiR9rDjezZOUzrKh7q02CCUfCrK1fm9GzNu/wlvXs3L+zzbmxQ49hduW31ZUlPcfWpfDWZRDeC3UPw8mP5yyRJKuurubkk0+OJ5Bkf/3rXwECn8WlMZFeJtbSaIm0WuBPUUEhDlKOYfihqKBQCUR6noYXvAQC3veGlJWUcmLbtm2MHDmyw/OlpaWBD8KrJdKDJU7fBW8sY/mm2lYtkJjmSAslRSW+xzRiwHCOHno0F0XHWTKZYiwSmGHneC2Q8F4I9fPeB2To0KGtxj9Snd+zZw87d+4MNJEoifRQiYPnS1YuwUGbwfBksQFxv4QKQm0SyMEYn+Frky5nTuUsX2MQ6ZYR070urByPiaTy+c9/nt///vesXLmSCRMmtDkfGx9ZuHAhV111Va7Di1MS6aFaLwjsOHnkSjgS5rV/vM5bdW9x27Rb2yxafGz5bzlxxAlqkUh+GzE90OQR8/Wvf517772Xc845h1tvvZXx48fzwQcf8N5773HnnXcyfvx4Zs2axbx586ivr+fMM89k+/btLF68mEWLFuUsTo2J9FD9i0uDDqFdsTUpp5ZPImSh+PGwi7S7VkUl5EVaKykpobq6mi996UvcfPPNfPGLX+Suu+7iiCMOLvq97777uOWWW3jssceYNm0ac+fOzXlZfHMuufJ671ZRUeFqa2uDDqPbfvD0D3ktzam4uRYqCHHHBT/hjLFTuL9mAY8t/y1hF6GksDjlqvnEbq/2rhFpz+rVqzn++OODDqPH6+jf0czedM5VpDqn7qw80t4g9MsbXmHJymfYtreRof2G8OUJF6bciCVfhCNhVm19lzPGTmFO5SxOHHFCh4PrqrMl0nMpieSJ9upcvbzhFW589uZWq8Jf+8eyvN+E55UNNfFB9DPGTukwKZxaPolnVz0bb4locaJIz6EkkmPttTba+2v8168/3KasiMOR792QU8ZWpj29N1ZCXlOBRXoeJZEcSm5tzJh4Sby8eqq/xu+vWcCa+veCDrtLQhbizM+cwYkjTuhSBeHOWisikp80OyuHklsbjy3/XXzTJoDbpt3K+LJjGdx3MKu2vsszby/t4NPyU9iFee2D13g6ocJvlysIi0iPoZZIDiW2NkJWQNh5pUliv2Qb9zayNtryePSN3wQZarc0tezH8IovapxDpHdTEsmhxL7//sWlPLHiyWhCCdG4t9H3FeV++ueR/8S6hnXxpPHlCRfGdzjUOIdI76UkkmPJff+/eeMxwi7M/133UoBRdd+xZcdy5aQr2iSN2HfV0BLpnTQmEpCXN7zCn1e/QCS6/0ZsK9qeKGQF8eRwTdX3AVqtPo9NKIiN/2hVukjvoZZIjsVKtdd++PKZGUwAABPFSURBVGbKHQF7oq9NuqJViyN5VpYWE4r0XkoiOZRq29qebsSA4cypnBXvrtqy46M2CUOLCUU6N3PmTN555x16WlkmJZEcSrltbQ83rP+wVskxcd/1WMLQYkKR3ktJJIf6F5dGp/b23PGPZG9veYeHXn84nhybIy2cPuY0jhh0RJsBdiUPkd5HA+s58vKGV3hixZOEXYSQhejfp3/QIWXNjn07KCksBry1IRdNuJBrqr6vpCGSgb/85S9MmDCB0tJSpkyZwqpVq+LnwuEwd9xxB8ceeyzFxcWMHj261R7rU6dO5eKLL2bBggWMGTOGvn37cv7557N582bf4lVLJEeSN2jafWB3wBFlz7nHn9NppV4R6dymTZv44Q9/yI033kjfvn255ppruOSSS3jnnXcwM2bPns3ChQu59tprOeuss2hsbGTx4sWtPuO1115j7dq1/OIXv6CpqYnrrruOL3/5yyxf7k/VCCWRHHh5wyus2vIOhuHI78KJXfUv46a2qtYr0tMFuaapsbGRmpoaxo0bB0AkEuGiiy5i7dq1ADz00EPcc889XH311fF7ZsyY0eoz6uvrefXVVykvLwegvLycKVOm8Oc//5nzzjsv6zEH3p1lZueZ2VozW29m13dw3SQzC5vZxQnHBpvZYjNbY2arzez03ESdvlgp9zX17/WaBFJc2IfPjTmNu6bfwU8v+HHQ4YhkTdBrmsaMGRNPIAAnnHACAHV1dfz1r38FaNV9lcrEiRPjCQSgsrKSsrIy3njjjewHTMBJxMxCwL3AF4ETgMvM7IR2rvsZ8HzSqXuAPzvnjgNOAlb7G3F6Erd6XbZxea9ZDxKzv+UAyzf1rGmIIulItaYplwYPHtzqfZ8+fbxYmprYtm0bpaWlDBw4sMPPKCsrS3lsy5Yt2Qs0QdAtkcnAeufc+865A8Ai4MIU130P+ANQHztgZgOBM4GHAJxzB5xz2/0PuWPJf8n0Ly4lVBDq/MYepjnSosq80uucWj6p1SSRfFrTNHToUPbs2cPOnTs7vK6+vj7lsZEjR/oSV9BJZBTwYcL7uuixODMbBVwE3J907zFAA/Cwmb1lZg+aWWmqh5jZLDOrNbPahoaG7EWfQvJfMs+8vZRwJOzrM4NQVFBI/+LSVuVNRHq62Jqmr570lU73wMm1qqoqABYuXNjhdStWrGDTpk3x9zU1NdTX1zN58mRf4gp6YD3VHq/JAwfzgeucc+GkLWELgYnA95xzy8zsHuB64KY2H+jcAmABQEVFha8DE4mrswG279vh5+Nyqk+oD0cPPZoh/Q7j2LJx8SrE6Ww6JdJT5OuapvHjxzNr1izmzZtHfX09Z555Jtu3b2fx4sUsWrQofl1ZWRkXXHABt956a3x21sSJE30ZVIfgk0gdcGTC+9HAR0nXVACLognkcGCambUArwN1zrll0esW4yWRQMX+knmg5lds2PZ+0OFkrLRPP4oLS2jc2xg/Nqz/MB654kHAK7CoelgiuXXfffdRXl7Ogw8+yJ133klZWRlnn312q2tOP/10vvCFLzB37lwaGhqYOnUqCxYs8C2moJPIcmCcmR0NbAYuBS5PvMA5d3TstZk9AvzJObck+v5DMxvvnFsLfB54N1eBdyT2y/T6pTfEq/T2NPuamxg9eHSrJHLUYQfzvephiWTXI4880ubYmDFjcO7g75BQKMQNN9zADTfc0OFnzZkzhzlz5mQ7xJQCTSLOuRYzuwpv1lUI+LVzbpWZzYmeTx4HSfY94Ldm1gd4H/imrwF3wRljp3Dl5K/12B0KIy7C0H5DKCoopDnSQlFBIRdNODjnQfWwRASCb4ngnHsOeC7pWMrk4ZybmfT+b3jdXXlpTuUsXv/HG6ytXxt0KCnF6ngVYG1aTOnsTpivfccikjuBJ5He7t9Om8mNz96cd2tFjis7ln897ZtttuotKihk0lEVfHnChW12JxSR/Pbiiy/m/JlKIj47Y+wUbj//NpasfIbX/rEsb1atjxx0RPz1iSNOUNeUiGRESSQHzhg7hSUrn8mbBALw7tbVvPbBa62m6Ma2thURSVfQiw17vZc3vMK8p3/Iso3+1K3J1O79uwIt7yAivYNaIj7qaDvcUEGISCQSWOtkz4G98deaoisimVJLxEftbYdrGCWh4pQJpKigKBehtTJj4iUaBxGRjCiJZFFi9V5oXcytwA7+Uzsce5r3trk/VBBiytjKrMUzdugx3DX9Du6afgefG3MaIUv9P/fu/Xuy9kwRObQoiWRJqn0IEou5XTnpinhCaY+LRPhoR/fKNVu0HFlJYTFTxlbGxzpGDjoi5d7u6soSyQ8zZ86koiJvl721S2Miaepst7NU+xDEFuMlXv+b5b8lkuKXOUAExweffNCtOI8YNJLTxpzWau3Hs6ueZcbESygpLI6XKZkx8RJ279+jKb0i0i1KImlIHCBvr2Jte7WkEpPP7v172k0gMQciB+KvDeOwfoe1ql8FXrfXgOL+7Gra1aZ1Ub+7gVPLJ7VJarv379FaEBHJOnVnpSGd3c5S7UOQaoOqzrq0YoYPGE5hQYjGvY2ECkIMHzCcfxk3la+e9BXuuOAn/PecP3HHl27nqyd9hfFl4+P3NYeb44miqMD7G6GooDCeOK6p+r4SiEgeW7JkCccddxwlJSVMmTKFd989WFc2HA5zxx13cOyxx1JcXMzo0aPbbJf79NNPM3nyZPr27cvQoUOZNm0aGzdu9C1etUTSkG7F2uSuq45aA/2LS+PdSQBLVj7D8k21NEdaKCksZuzQo/l418cAhCNhphxT2WYxYOx5sX3cm8PNFIWKDsYX23/FUm3bIiL5ZuPGjfzgBz/gxz/+MX379uWWW27h3HPPZd26dZSUlDB79mwWLlzItddey1lnnUVjYyOLFy+O3/+b3/yGr3/961x66aXcdNNNOOeorq6moaGh1b7r2aQkkoZMK9amSj7tFS2MJYPYM1ZtfZdX//F6/Hz/4pSbNh4UKxcd/Z64t3usdaIWiEjnGqur2V5Tw+DKSoZEdxPMlU8++YRnnnmGz33ucwCccsopjB07lkceeYSpU6fy0EMPcc8993D11VfH75kxYwYAkUiE66+/nosuuojHH388fn769Om+xqwkkqZMKtZ2NfkkPiO5y6yjabjLNi6nOdICHNz7XPt9iHRdY3U16+bNI9LURMNTTzHu7rtzmkjKysriCQSgvLycU045hTfeeCO+r0hy91XM2rVr+eijj/jmN3O7I4aSiM8yLZfelSSQfG3/4lKWbVyuGVgiXbS9poZIUxMAkaYmttfU5DyJpDq2ZcsWtm3bRmlpKQMHDkx577Zt2wAYOXKkrzEmUxLJI8nTiNNtxZwxdgozJl7CKxtqOGrIUfGpvSWFxdr7XKQLBldW0vDUU0SamigoKWFwZfYW/6ajvr4+5bETTzyRoUOHsmfPHnbu3JkykQwdOhSALVu6t9asqzQ7K0+0t1gxndlUL294hSdWPMmGbe/z0vqXVFhRJENDqqoYd/fdDL/88px3ZYGXMF599dX4+02bNrFixQomT55MVTSWhQsXprx3/PjxjBo1ikcffTQnscaoJZIn2lus2NV7wy5CyEKEXVhjISIZGFJVlfPkEXP44Ydz5ZVXxmdn3XzzzZSVlTFz5kxKSkqYNWsW8+bNo76+njPPPJPt27ezePFiFi1aREFBAXfddRdXXHEFV1xxBZdddhlmRnV1NZdddplvq+GVRPJEdwbCk+/VWIhIz1ReXs4NN9zA9ddfz8aNG6moqODxxx+npKQEgPvuu4/y8nIefPBB7rzzTsrKyjj77LPj919++eWUlJRw++23c/HFF1NaWsppp53GsGHDfIvZYiP+h4qKigpXW1sbdBgpdVZapaPzsXOJ60+UQORQsHr1ao4//vigw+jxOvp3NLM3nXMpmzJqieSRjmZydVZ6Jfa6s/IsIiLZpIH1HiKd0ivpXCMikk1KIj1E4t4k7Y2ZpHONiEg2qTurB5k4+mQccNGEC9stnaJKvSKSS0oiPUDieEhJYTEXTbiw3WszXSEv0pM55zAVGs1YdyZYqTurB9BYh0j7ioqK2LdvX9Bh9Gj79u2jqKgoo3uVRHoAjXWItK+srIzNmzezd+/ebv1FfShyzrF37142b96csm5XOgLvzjKz84B7gBDwoHPuznaumwS8Dsxwzi1OOB4CaoHNzrkLchByzmmsQ6R9sTpSH330Ec3NzQFH0/MUFRUxfPjwdgs7dibQJBJNAPcCZwN1wHIzW+qcezfFdT8Dnk/xMf8BrAYy+xfoITTWIdK+gQMHZvxLULon6O6sycB659z7zrkDwCIg1ajx94A/AK1KXJrZaOB84EG/AxURkbaCTiKjgA8T3tdFj8WZ2SjgIuD+FPfPB64FIh09xMxmmVmtmdU2NDR0L2IREYkLOomkmpOXPDI2H7jOORdudaPZBUC9c+7Nzh7inFvgnKtwzlX4WYhMRORQE/TAeh1wZML70cBHSddUAIuic8APB6aZWQtwKjDdzKYBJcBAM3vMOfc1/8MWEREIPoksB8aZ2dHAZuBS4PLEC5xzR8dem9kjwJ+cc0uAJcCPosenAtcogYiI5FagScQ512JmV+HNugoBv3bOrTKzOdHzqcZBREQkTxxy+4mYWQOw0efHHA584vMzuivfY8z3+EAxZkO+xweKEaDcOZdyQPmQSyK5YGa17W3gki/yPcZ8jw8UYzbke3ygGDsT9OwsERHpwZREREQkY0oi/lgQdABpyPcY8z0+UIzZkO/xgWLskMZEREQkY2qJiIhIxpREREQkY0oiXWBm55nZWjNbb2bXpzhvZvaf0fMrzWxiwrnvm9kqM3vHzB43s5KAYjzOzF4zs/1mdk1X7g06RjM70sz+amaro/+W/5FP8SWcD5nZW2b2Jz/i626MZjbYzBab2Zrov+XpeRhjvvy8XBH9WV5pZq+a2Unp3htkfLn6WQG8na301fkX3or6DcAxQB/g78AJSddMA/4br7DkacCy6PFRwAdA3+j7J4GZAcVYBkwCbscrFZP2vXkQ40hgYvT1AOC9bMfYnfgSzv8A+B1eiZ6g/r/YbozAo8C3oq/7AIPzKcY8+3n5HHBY9PUXE36mff956WZ8vv+sxL7UEklfOnufXAgsdJ7XgcFmNjJ6rhDoa2aFQD/aFprMSYzOuXrn3HIgeQu4dPd2CSxG59wW59yK6OtdeJuRtdo6IMj4IGd73GQco5kNBM4EHoped8A5tz2fYozKl5+XV51zn0bfvo5XJDate4OML0c/K4C6s7qi071P2rvGObcZ+DmwCdgC7HDOvRBQjH7c2xVZeY6ZjQFOBpZlJaqDuhtfWnvcdFN3YjwGaAAejna5PWhmpdkOkG7EmMc/L/+G19OQyb2Z6E58cT7+rABKIl2Rzt4nKa8xs8Pw/oI4GjgCKDUzPyoOpxOjH/d2RbefY2b98Xa6nOuc25mVqBI+PsWxtOKzLuxx003d+TcsBCYC/+WcOxnYA/jRn9+df8e8+3kxs3/B+yV9XVfv7YbuxBc77ufPCqAk0hXp7H3S3jVfAD5wzjU455qBp/D6MoOI0Y97u6JbzzGzIrwfit86557KcmzQvfgq8fa4+Qde10OVmT2W3fCA7v/vXOeci/1VuhgvqWRbd2LMq58XM5uA1z15oXNuW1fuDTC+XPysAEoiXRHf+8TM+uDtfbI06ZqlwNfNcxpeM3wLXrP8NDPrZ2YGfB6vjzKIGP24NycxRv/tHgJWO+d+4UNs3YrPOfcj59xo59yY6H3Vzp89broT41bgQzMbHz30eeDdfIqRPPp5MbOj8JLYlc6597pyb5Dx5ehnxePHaH1v/cKbffUe3oyJG6PH5gBzoq8NuDd6/m2gIuHe/x9YA7wD/AYoDijGEXh/4ewEtkdfD2zv3nyKEZiC15xfCfwt+jUtX+JL+oyp+DQ7Kwv/O38WqI3+Oy4hOrsnz2LMl5+XB4FPE/7/VtvRvfkSX65+VpxzKnsiIiKZU3eWiIhkTElEREQypiQiIiIZUxIREZGMKYmIiEjGlERERCRjSiIiIpIxJRGRHDKzQjO72sz+bmZNZrbVzH4ZXZ29w8z8WD0u4pvCoAMQOVRES1f8ETgHb8X4fwKHA/+KV113IODbRlYiflASEcmdX+IlkB86534eO2hmjwIvRt+uCCAukYyp7IlIDpjZJOAN4AXn3Lkpzsd2sPu8c6461/GJZEpjIiK5cVX0+23tnI+V8H4r+YSZ3WBmzsx+6UtkIt2gJCKSG+cC25xzNe2cH4W3h8aniQejWwp8G68aq0jeURIR8ZmZlQDD8fbJSHX+n/B28Hsr6fgg4Ld4O9Z9muJWkcApiYj4Lxz9GtrO+Zuj35MH1RcAizVGIvlMSUTEZ87b4nUdcFR0L2zA233OzG4G/kf00FsJ574NfAa4KZexinSVpviK5MZdwK+BZ83scaARby/xAXjb055AtCUS3br2p8AZzrkDwYQrkh5N8RXJETP7AfA9vPGPj4DFwO14rZRm59wR0etmAg/jdYHFhPC2O40Apc65/bmLXKR9SiIiATKzI/EG3J9zzp0fPTYYGJ106cN4yeanwCqnH1zJE+rOEgnWydHv8UF159x2YHviRWa2B2h0zr2Tw9hEOqWBdZFgxZJIm0WGIj2BWiIiwWrTEknFOTfV/1BEuk5jIiIikjF1Z4mISMaUREREJGNKIiIikjElERERyZiSiIiIZExJREREMqYkIiIiGft/RBRNsIe4PFsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(fccaq[0], fccaq[1], s=10, label='fcc', color='#FFB300')\n",
    "plt.scatter(hcpaq[0], hcpaq[1], s=10, label='hcp', color='#388E3C')\n",
    "plt.scatter(bccaq[0], bccaq[1], s=10, label='bcc', color='#C62828')\n",
    "plt.xlabel(\"$q_4$\", fontsize=20)\n",
    "plt.ylabel(\"$q_6$\", fontsize=20)\n",
    "plt.legend(loc=4, fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voronoi based method also provides good resolution,the major difference being that the location of bcc distribution is different."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "lammps",
   "language": "python",
   "name": "lammps"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}