0_Synthetic data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Synthetic dataset\n",
    "\n",
    "This notebook generates a synthetic dataset that is based on the following conceptual framework: \n",
    "- Log interval from 500 - 1000 mMD in a vertical well. \n",
    "- Two lithologies: 500 - 750 and 750 - 1000 mMD. \n",
    "- Two fault zones with associated increased fracture frequency, centered at 600 and 900 mMD.\n",
    "\n",
    "A CSV file is generated that replicates an output from BHI analysis done with log analysis software."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import mplstereonet as mps\n",
    "\n",
    "from utilities import *\n",
    "from fractoolbox import strike2dipaz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "strike_cluster1 = random_strikes(100, 90, 10) # n, mean, std_dev\n",
    "dip_cluster1 = random_dips(100, 65, 10)\n",
    "depth_cluster1 = random_depths_normal(100, 600, 80, 500, 1000) # n, mean, std_dev, lower_bound, upper_bound\n",
    "aperture_cluster1 = random_aperture_lognormal(2.5, 3, 0.01, 20, 100)  #(0.02, 3, 0.01, 20, 100) \n",
    "attribute_cluster1 = 100 * ['Conductive, HC']\n",
    "\n",
    "strike_cluster2 = random_strikes(100, 280, 10)\n",
    "dip_cluster2 = random_dips(100, 65, 10)\n",
    "depth_cluster2 = random_depths_normal(100, 900, 50, 500, 1000)\n",
    "aperture_cluster2 = random_aperture_lognormal(0.1, 3, 0.01, 10, 100)\n",
    "attribute_cluster2 = 100 * ['Conductive, HC']\n",
    "\n",
    "strike_background = np.random.randint(0, 360, 50)\n",
    "dip_background = np.random.randint(0, 90, 50)\n",
    "depth_background = random_depths_normal(50, 750, 250, 500, 1000)\n",
    "aperture_background = random_aperture_lognormal(0.02, 5, 0.01, 5, 50)\n",
    "attribute_background = 50 * ['Conductive, LC']\n",
    "\n",
    "strike_bedding1 = random_strikes(100, 170, 15)\n",
    "dip_bedding1 = random_dips(100, 15, 5)\n",
    "depth_bedding1 = random_depths_uniform(100, 500, 750)\n",
    "aperture_bedding1 = 100 * [0]\n",
    "attribute_bedding1 = 100 * ['Bed boundary']\n",
    "\n",
    "strike_bedding2 = random_strikes(100, 220, 30)\n",
    "dip_bedding2 = random_dips(100, 20, 15)\n",
    "depth_bedding2 = random_depths_uniform(100, 750, 1000)\n",
    "aperture_bedding2 = 100 * [0]\n",
    "attribute_bedding2 = 100 * ['Bed boundary']\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAGjCAYAAACBlXr0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHOUlEQVR4nO3de1yUZd4/8M+AgiAH8ZQKhoopBh4QsFTygAmGSXbW1VxL1HQ3i3R/uatmPmVPz2OluU8aiHnMyMcVgyTRTfHxsK2AaEhCNikopoYHBkRAmfv3BzssIIeZYWaue+77896XL5TDzBdz58P3ur73dWskSZJAREQkkIPoAoiIiBhGREQkHMOIiIiEYxgREZFwDCMiIhKOYURERMIxjIiISDiGERERCccwIiIi4RhGREQkHMOIiIiEYxgREZFwDCMiIhKOYURERMIxjIiISDiGERERCccwIiIi4RhGREQkHMOIiIiEYxgREZFwDCMiIhKOYURERMIxjIiISDiGERERCccwIiIi4RhGREQkHMOIiIiEYxgREZFwDCMiIhKOYURERMIxjIiISDiGERERCccwIiIi4RhGREQkHMOIiIiEYxgREZFwDCMiIhKOYURERMIxjIiISDiGERERCccwIiIi4RhGREQkHMOIiIiEYxgREZFwDCMiIhKOYURERMIxjIiISDiGERERCccwIiIi4RhGREQkHMOIiIiEYxgREZFwDCMiIhKOYURERMIxjIiISLg2ogsgsgfl5eW4fPkydDodSktLUVpairKysnq/v3PnDvR6fe0vAHBwcKj95eLiAjc3N7i7u8Pd3b3e7z08PNCjRw+4uroK/k6JxNBIkiSJLoJIFEmScPPmTRQVFeHSpUv3vTX8/ubNm2jbti08PT3rhUjdUHFxcYGjo2Nt+ACoDabq6mrcuXPnvhAz/LmkpAR3796Fl5cXfHx84O3t3eRbLy8vaDQawX9zRJbFMCLVqKqqQm5uLjIzM5GVlYWsrCz8+OOPKC8vh4eHB7y7eMKn3R149xkAn8Gj7wuCzp07Wy0EJElCcXFxo4FY961Op4OrqysefvhhBAcHIyQkBMHBwQgICICTk5NVaiOyBYYRKVJVVRXOnDlTGzpZWVn44Ycf4OLigqFDh9a+kA8cOBA9e/aEu7s7sDoQKLkIePYEYs+I/hYaVVpaiosXL+KHH36o/b5OnjyJO3fuYNCgQQgODq793hhQZE8YRmT3jA2e4OBg9OnTp3YJ7T4ZG4Gjq4GwWCB0lm2/iVbQ6/X45Zdf6nV8jQVUcHAwAgMDGVAkSwwjsktFRUVISUlBcnIyDh48iHbt2pkWPArXVEBVVFQgPDwc0dHRmDRpEry9vUWXSgSAYUR2QpIknDp1qjaATp06hREjRiA6OhoTJ05E//79ZRk8NxMTURy/AZ3nzIbXlClCa9Hr9cjPz8fevXuRnJyM48ePY8iQIYiOjkZ0dDQGDx7MwQgShmFEslVZWYn09HQkJycjJSUFN2/exIQJExAdHY2oqCh06tRJdIktOhc+DvcuX0abHj3w0MHvRJdTT3FxMVJTU5GSkoJ9+/bBy8sLkyZNQnR0NMaMGQNnZ2fRJZKKMIxIVhp7gTQsKdnjC6ScOqPmKCH4yb4xjEg4nU6H7du3IzExkUtHMiBJEk6fPo3k5OR6S6JTp07F9OnTayYPiSyMYUTC5OTkYP369di2bRsGDBiAl19+GdHR0dxUl5lLly4hJSUFn3/+OfLy8jBjxgzMmzcPgYGBoksjBZHfji8pWlVVFRITEzFq1CgMGzYMFRUVOHjwIE6cOIF58+YxiGTIx8cH8+bNQ0ZGBg4ePIjy8nKEhoZi1KhRSExMRFVVlegSSQHYGZFNFBYWIj4+Hhs2bICHhwfmzZuHmTNnomPHjqJLszp72TcyxfXr17F582asX78eZWVlmD17NubMmYOePXuKLo3sFMOIrEav1+Pvf/871q1bh9TUVERFRWH+/Pl4/PHHZTmGbS1ynqhrLb1ejwMHDmDdunX49ttvMXHiRMyfPx/jxo1T1X9jaj3+ayGLu3HjBj7++GP0798fL730EgIDA/Hzzz9jz549iIiIUNSL1M3ERJwLH4ebiYlNfk7nObPRpkcPdJ4z24aV2YaDgwMiIyPx9ddf49y5c3j44Ycxbdo0+Pv7Y/Xq1bh586boEslOsDMii7l27Rref/99xMXFISQkBPPnz8czzzxjd+PYplBy12OuyspK7N69G+vWrUNWVhbmzJmDJUuWoEuXLqJLIxlTzo+oJExpaSneeecd+Pn5QavV4vjx4zhy5AimTp2q6CACbNv1GNOFGfM51ubs7IypU6fiyJEjOHbsGLRaLfr06YMVK1agtLRUWF0kb+yMyGyVlZWIi4vDe++9h379+uGDDz5AWFiY6LJkw9KDC8Z0YXLt1I4cOYLFixfj3LlzWLp0KebOnav4H1TINOyMyGTV1dXYtm0b/P39sWHDBmzcuBFHjhyxyyCyRidxMzER+Y88iisr/gP3Ll9GcfwGizyuMV1Y3c+RQ5dk8Nhjj+Ho0aPYuHEj4uPj4e/vj+3bt9feEZeInREZTZIkpKam4s9//jN0Oh3+4z/+A9OmTYOjo6Po0sxyMzERV959D6iutmgnYehOAACOjui2bGmTnZE1x77l2iVVV1fjiy++wNtvvw1PT0+8//77iIqK4kkbKsfOiIxy/PhxjBo1CjNnzsSsWbOQn5+PGTNm2G0QAajpWKqrAUdHi+75dJ4zGw6ennDw9Gw2iADg2uo1uHf5Mq6tXmOx569bhxyn+BwdHTFjxgzk5+fjlVdewcyZMzF69GgcP35cdGkkkkTUjJycHCk6Olpyc3OTli9fLul0OtElme3Gl19KP40Nl258+WWjfxYhb9gj0o/9/aW8YY/Ue78carOVkpIS6e2335bc3Nyk6Oho6cyZM6JLIgHYGVGjSkpKMHfuXISGhqJXr17QarV455137PqQzOL4DfX2cLymTMFDB78TeipC19g30KZHD3SNfaPe+xvWqmQeHh5YsWIFtFotfH19ERISgrlz56KkpER0aWRDDCO6T1paGgIDA3H+/Hnk5ubik08+QdeuXUWX1Wq2WrZqbnCg4ceaCkTDUp/+9m1ZDCDYQteuXbF27Vrk5ubi/PnzCAwMRFpamuiyyFZEt2YkH7du3ZJmzZolubu7S/Hx8ZJerxddktWZshxm7Of+NDZc+rG/v/TT2HCTPtaaz1UavV4vxcXFSe7u7lJMTIx069Yt0SWRlbEzIgD/7oYKCwtx5swZzJ49WxXTTcYuhxkm74z53OY6MFO6s+Y+V05j29ag0WgwZ84c5OTk4MKFC+yS1EB0GpJYauyG6jK12/nx4QCjBiCs0XE1Vk+9rulEgiR9HFDzVkHYJakDrzNSsbS0NMTExGDAgAFISEjAgw8+KLok2WrseqC61/F0njO73sdNucbH2M+tWwNQ09W5BgWhPDu75nl/fQ8ouQh49gRiz1jum5eJgoICxMTEIC8vDwkJCYiMjBRdElkQl+lUqKSkBDExMXj++eexbNkypKWlMYha4DVlSm3gGJbG6i6jGZb7rrz7Hm4mJlpsOa6uhkuK+tu3oUtN/ff7wmJrgigstvXfsAz5+vpi//79WLZsGZ5//nnMnj2bE3cKws5IZQzdkL+/PxISEuDr6yu6JLvRXAdTtHARdHv3AoBR3ZA5Jy/U/RpDMAFo8ZQHUx7XXm7+xy5JedgZqURZWVm9bmj//v0MojqaGgio+/6GHUzdj5VnZ9d8QTOnOdT9fHOuI6o7Bm7KKQ8tscdrmgxd0tKlS2u7pLKyMtFlUSswjFRAq9Vi+PDh+Pnnn5GTk4M5c+YodlLO3Cmzpl6Qm3uhrvsxQ1DVDYaGtRg+/9rqNdDfvg0HT0+zr3nymjIF/f/5Pfr/8/tWdzNyPTaoJRqNBnPnzkVOTg7OnTuH4cOHQ6vVii6LzMQwUrj9+/cjNDQU48aNw4EDBxTfDZn7U35TL8iN7QsZHrvuxxq7eLXhuXOGzwcAfUkJHNq3l8WymBxOomgNX19fHDhwAOPGjUNoaCgOHDgguiQyh9hhPrIWvV4vrVq1SnJ1dZU2b94suhybseaZbqY+dmvOnVPT2XSWtGnTJsnV1VVatWqV6i5TsHccYFCg8vJyzJ49G4cPH0ZSUhJCQ0NFl6RKzQ0GtDQ0INfbP9iDEydO4JlnnsHo0aOxYcMGuLq6ii6JjMBlOoUpKChAWFgYCgsLkZWVxSAygbn7TU19XXPLXy0tJxqW9FyDglp10oIxgxlKM2zYMGRmZtb+f6GgoEB0SWQEhpGCpKenIyQkBI888gi+++47PPDAA6JLkqWmXojN3W+q+3XGvsi3NDRgCLLy7OxWTbqZM5ihhKDq1q0bDh48iGHDhiEkJATp6emiS6KWiF4npNbT6/XS2rVrJVdXVykuLk50ObJn0XsInUiQ7r3bS/p1ch/pR/8B0tkhQWYfbtrY87d276ipr2/ucZV2QOtnn30mubq6Sn/961+5jyRj3DOyc3q9HgsXLsSOHTuQlJSEESNGiC7Jpsy5YPPKc/3RsfdvuHG+C7rtym9dAasDgZKLqLrtCG3KA4BGgzbduzdZT3P1ymWfyB4vgm3JsWPH8Mwzz2DatGn48MMP4eDARSG5YRjZsaqqKsycORMnTpxAWloa/Pz8RJdkc+a8gFe/1xuO926guk1HOC4937oCMjYCR1fj5tUHcSWlEB4TJsD7ow/NqleJISAnWq0WERERePTRR7Fp0yY4OTmJLonqYBjZqdLSUjz77LO4fv06UlNTVbs/ZNYL+L8CBGGxQOgsq9fT2AGnDBwxrl69iieeeAJdunTB3/72N7i5uYkuif6FYWSHrl27hokTJ8LT0xNJSUl2fStwpWms85HL8pul2Wsnp9Pp8Mwzz0Cn02Hv3r3o0qWL6JIInKazO+fPn0dYWBj8/Pywd+9eBpENmDJd1tiUXEtj2jcTE5H/yKPIf+RRu5pgs8cz7QDAw8MDe/fuRe/evTFy5EhcuHBBdEkEdkZ25aeffsLYsWPx7LPPYs2aNdyEtRFLdTZNPY7h/YBxJ37Lhb12RgZ6vR6vv/46kpKScPDgQfTr1090SarGVzM7cfbsWYwePRovvfQSPvnkEwaRhTXX/TTW7ZhzLU5z598ZTuA297BSEdcG2fuZdg4ODli7di2mTZuGMWPGIC8vT3RJ6iZqppyMd+bMGalr167SkiVLeJ2ElZh6bU1T1yoZmHN9kDnXBJlbP/2bXq+X/vKXv0hdu3aVzpw5I7oc1eKP1zJ39uxZhIeH49VXX8W7776r2Fs/iGap2ygYOhTDid2m7KeYc1qCgb3eBkIONBoN3nvvPcydOxfjxo1jhyQIw0jG8vPzER4ejtmzZ+Odd95hELVCS8tYJi05ZWxE30lX0Sm4DbrGvlHvQ3XvwOrg6Qn97dtGL50ZcxuLptj7kploGo0GK1aswKxZsxAeHo6ffvpJdEmqwwEGmTJMzb300kv4z//8TwZRK1l0vPpfpy7AsycQe6beh24mJtbevwiouW+R4TmLFi6Cbt++Fi+MJXEkScLixYvxxRdf4OjRo+jVq5foklSDnZEM/fbbb4iMjMSzzz7LILIQiy5jhcXWBFFY7H0f8poyBQ7t20NfUgIAtc95MzERur17gepq6Pbta/bhlXBQqb3SaDT44IMP8PTTTyMyMhK//fab6JJUg52RzJSVlWHs2LHw8/PDjh07ODVnhxobea47vu0xcaLZRwaRbej1ekydOhXnz5/HwYMHeVKDDTCMZKSqqgqTJk1CdXU19u7dC2dnZ9ElkYWYck2OvV+/oxSVlZWIiopCmzZtkJKSwrPsrIxhJBN6vR4vvfQS8vLycOjQIXh4eIguya7J5QVdLnWQeXQ6HcaMGYOHH34YW7du5UqFFfFvVgYkScLChQvxz3/+E6mpqQwiCzBMtV1bvQbXpg1E9Xu9aw5ItYLm9njs9cgcquHh4YFvv/0W33//PRYtWgT+7G49DCMZ+J//+R/s2LEDaWlpqj1929IMAwsA0KFHERzv3ag5qdsKmgscS53eQOI88MADSEtLwxdffIFPP/1UdDmKxTASLD09HYsXL0ZSUpIq70dkLYbrbrrGvoFbl71R3aZjo9NvDVnymJ+GS3QNL4i9tnoNQ8lO+Pn5Yffu3Xjrrbd4C3Mr4Z6RQAUFBQgJCcH777+P2bN55bwcNJxka82eT8PHMvzZwdOzZvz79u161yGR/MXHx2PJkiXIysrCgw8+KLocRWFnJEh5eTmefvppPPfccwwiGWnY5bRmz6fhYxn+3DX2jdqujUf42Jc5c+bg2WefxeTJk1FeXi66HEVhZySAJEmYNm0aCgsLcfDgQY6MylhjS22cjlO3qqoqhIeHw9fXF9u3b+dF6RbCMBLgww8/xJo1a5CZmYlu3bqJLodMwAtSCQCuXLmCkJAQxMbGYuHChaLLUQQu09nY/v37sXz5ciQlJTGI7JCxxwpxYk7ZunXrhqSkJLz99ts4cOCA6HIUgZ2RDWm1WoSGhmLNmjWYMWOG6HLIimqP/3F0RLdlS7mkp1CbN2/Gm2++iczMTPTp00d0OXaNnZGNlJWVYfLkyZgxYwaDSCFaujssHB2B6mpe8KpgM2fOxEsvvYTJkyejrKxMdDl2jZ2RjcTExODnn3/GgQMH0LZtW9HlkAW0tH/EYQd1uHv3LsaPH4+HHnoIGzbwBw9zMYxsIC0tDc8//zxycnLg6+sruhyyEIYNGRQUFGDgwIHYtWsXIiIiRJdjl9qILkDpSkpKEBMTgw8//JBBpDCGADIswzGQ1MvX1xerVq1CTEwMcnJy4OnpKboku8POyMpiYmJQUFCA/fv383oEBeKoNxlIkoTx48ejd+/eXK4zAwcYrCgtLQ07d+5EQkICg0ihLHoHWbJrGo0GGzduxFdffYX9+/eLLsfusDOykpKSEgQGBmLZsmWYM2eO6HLIBriHRAAQFxeHlStXcrnORAwjK4mJiUFhYSHS0tLYFakEl+wIqFmui4iIQK9evbhcZwIu01kBl+fUiUt2BNQs1yUkJHC5zkTsjCyMy3NEBHC5zlQMIwvj8hwZcA9J3bhcZxqGkQWlpaXhhRdeQE5ODm+8RdxDIl4MawLuGVlI3YtbGUQEcA+Jai6G/fDDDxETE4OSkhLR5cgaOyMLmTt3Ls6fP8/lOSKqx7Bc16dPH8TFxYkuR7YYRhZw5swZhIaGIjc3l8fIE9F9fvnlFwQEBCAzMxMBAQGiy5ElhpEFPPXUU/D19cXatWtFl0JEMvXaa6/h4sWL2LNnj+hSZIlh1ErHjh3DhAkToNVq0bVrV9HlEJFMXb16FX379kVaWhpGjBghuhzZYRi1giRJGDVqFMaNG4d33nlHdDlEJHPLly/HoUOHcPjwYe4tN8BpulZITU1FXl4e3nzzTdGlkJGauzsrkbUtXLgQZ8+exbfffiu6FNlhZ2Sm6upqBAUFYdasWXj99ddFl0NG4rU/JNqaNWuwadMmZGdnw8GB/YAB/ybMtGPHDuh0Orz66quiSyET8NofEu3VV1/FrVu3sGPHDtGlyAo7IzNUVlbC398fK1aswIwZM0SXQ0R2ZsuWLXjnnXeQn58PJycn0eXIAjsjM8TFxcHNzQ3Tpk0TXQoR2aHp06ejffv2vAi2DnZGJiotLYWfnx82btyISZMmiS6HiOxUcnIyYmJioNVq4e7uLroc4dgZmejjjz9Gv3798OSTT4ouhYjs2KRJk/DQQw9h9erVokuRBXZGJrh27Rr8/Pzw7bffIiwsTHQ5RGQNGRuBo6uBsFggdJZVn+rIkSOYOHEitFotunTpYtXnkjt2RiZ4//33MWbMGAYRWV7GRmB1YM1bEuvoaqDkYs1bK3vssccwatQorFy50urPJXfsjIx048YNeHt74/jx4wgKChJdDinN6sCaF0DPnkDsGdHVWIYNOwyLsnHd2dnZGDlyJIqKiuDl5WX155MrdkZG2rx5M0JCQhhEZB1hsTVBFBYruhLLsWGHYVGhs2p+ILBRgAYFBWHo0KHYvHmzTZ5PrhhGRtDr9Vi/fj3mz58vuhRSKhu/ANqEEgPWSubPn4/169dDr9eLLkUYLtMZYf/+/XjppZdQWFgIZ2dn0eUQkT1qZvmvsrISPXv2xBdffIHx48cLKlAsdkZGWLduHWJiYhhERGS+ZpYtnZ2dERMTg3Xr1gkoTB7YGbWgsLAQffv2xblz5+Dr6yu6HCKyVy0MRhQUFOChhx6CVqtFz549BRQoFsOoBUuXLkVOTg6+/vpr0aUQkcJFR0dj8ODBePfdd0WXYnMMo2ZUVVWhZ8+e2Lp1KyIjI0WXQ0QKt2/fPsycOROFhYWqO0CVe0bN2L17Nzw8PFS7oUhEthUREQE3NzckJSWJLsXmGEbNWLduHebNm8cbYBGRTTg4OGDevHmqHGTgMl0TcnJyMGzYMBQVFaFjx46iy2mZvV7tTkT1XL9+HT4+PsjIyEBgYKDocmyGP/I3Yf369Zg6dap9BBFgv1e7E1E9nTp1wpQpU7B+/XrRpdgUO6NG6HQ6eHt749ChQwgJCRFdjnHYGREpRkZGBsLDw3H58mXV3OuInVEjtm/fjgEDBthPEAHKPE6GLIsng9uN0NBQ+Pv7Y/v27aJLsRmGUSMSExPx8ssviy6DyLK4lGtXXn75ZSQmJoouw2YYRg0UFxfj+PHjvKU4yYMluxkeXGpXoqOjcezYMVy/fl10KTbBMGogNTUVQ4YMgY+Pj+hSiCzbzXAp1674+Phg8ODBSE1NFV2KTTCMGkhJSUF0dLToMohqsJtRtejoaKSkpIguwyY4TVdHZWUlOnfujCNHjmDIkCGiyyEilcvOzsbo0aPx22+/Kf6uAeyM6khPT4eXlxcGDx4suhSSI06jkY0NGTIEnp6eOHz4sOhSrI5hVEdycjImTZoEjUYjuhSSI06jkY1pNBpMmjQJycnJokuxOobRv0iSxP0iah73b0iA6OhoJCcnQ+k7Ktwz+pdTp07hscceQ3FxseLXZonIfhj2so8eParoLQR2Rv+SnJyMCRMmMIio9bi3RBbk7OyMyMhIxS/VMYz+JTk5mUt0ZBncWyILMyzVKRmX6QBcunQJvXr1wtWrV9GpUyfR5ZC946G1ZGHFxcXo1q0bCgoK4O3tLbocq2BnBOCbb77BiBEjGERkGTzpgCysc+fOGD58OL755hvRpVgNwwhcoiMi+VP6Up3ql+mqqqrg4eGB7OxsDBgwQHQ5RESNOnv2LIKCgqDT6eDk5CS6HItTfWd05swZtGvXDv7+/qJLISJqUv/+/eHs7Izc3FzRpViF6sMoKysLQ4cO5akLRCRrDg4OGDp0KLKyskSXYhUMo6wsBAcHiy6DiKhFwcHBDCOlysrKsq/bixORaoWEhCg2jFQ9wFBVVQV3d3fk5uaib9++osshImrWuXPnEBgYiNLSUsUNMai6Mzpz5gxcXFzg5+cnuhQiohb5+fmhXbt2ihxiUHUYcXjBTvHsN1IpJQ8xqD6MOLxgh3j2G6mYUocYVB9GHF6wQ7yvEKmYUocYVDvAwOEFIss4fSAVJ/bswrDJz2Hw+CjR5SjeuXPnMHDgQJSWlqJt27aiy7EY1XZGubm5HF4gsoATe3ZBV3wNJ/bsEl2KKvj5+SnyJAbVhlFmZiaHF4gsYNjk5+DRuSuGTX5OdCmqYBhiyMzMFF2KRak2jDi8QGq3M38nInZFYGf+TqPe35TB46Mw+9PPuURnQ0ocYlB1GHF4gdRqZ/5OrPznSvx6+1ck5CTU+1hCTkKj7yf5UOIQgyrDSJIk5ObmYuDAgaJLITKLqd1LQwk5CdBLejhoHBAzMKbex2IGxqB7++73vZ/kY9CgQcjNzYWS5s9UOU1348YNdOrUCTqdDu7u7qLLITJZxK4I/Hr7V3Rv3x37n9tv8tfvzN+JhJwExAyMwQv9X7BChWRNOp0Onp6euHHjBry8vESXYxGq7IyKiorg4eHBICK71dru5YX+L2D/c/sZRHbK8PpVVFQkuhSLUWUYXbp0Cd7e3qLLIDJbU2HS2uU7sh/e3t64dOmS6DIsRpVhVFRUBB8fH9FlEBnFlIAxDB+s/OdKBpLC+fj4sDOyd+yMyGYscKirKdNtMQNj4KBxgF7St/j57KLsGzsjBWBnRDZjyqGuTQSXKftDL/R/AUseWdLo5zcMH45w2zd2RgrAzohsxpRDXZsIrob7Q3VDxZwlPEP41Iac+wDeksMOsTNSAHZGZDOhs4DYMzVvW2JkcNXdF1p7cu19e0RNdTwNO6zakPvxO96Sww6xM1IAdkYkS0YGV919IV2VDhpo6u0RNbWs1+QE3sPjEObbEyM7t+P+kR1RWmekuotey8vL0b59e1y7dg1dunQRXQ6R0epeqAoAK/+5EnpJD08nTwBAZXUlAMDZ0RkLhi5o8Roiw+PdvnsbuiodAJh9ES3Z3rVr1/DAAw+gvLwcLi4uostptTaiC7C1y5cvo23btujcubPoUoiatDN/Jz45+Qk00NQGS93lN0NgGMIpIScBJVUlAICK6gqs/OdKALgvkN5Keg77SvIwwKUbzlb+Vhtmnk6ekCDxCCA70rlzZ7Rp0waXL19WxK1wVLdMZzhGg7eOIDlLyEmArkqHkqqS+5bfgroGIWJXBADULrvFDIyBp5Mn2jm2AwDoJT1WZay6b7hhX0ke9BoNcu9cqT2bbsHQBTg69SiOTT3GExnsiIODAzw9PaHT6USXYhGqC6PS0lK4ubmJLoOoWTEDY+Dh5AFPJ8/7Bg6yr2XXG1AwLLctGLoAGdMz4OHkAaCmQ/r19q/45OQntY87wdMfDpKEAJdu6N6+O5Y8soQBZMfc3NxQWloqugyLUGUY8Uw6EsmYcewX+r+AY1OP4ejUo/eFRcMBhYbTc68PfR0Omn//X1sDTe3zZlfrsGT420h88e88m04B3N3dGUb2qqysjGFEQhlzsWlzgfVC/xewyvMNlKzdh9MHUhHUNQgaaHCz4ibe+r+3kJCTgAm9JtR2VguGLjD6ecm+uLu7o6ysTHQZFqG6AQZ2RiSaYeCguWGButcSZV3NQva1bAR1DUL2tWzEDIxByZ590BVfw9df/BXfji2CBAkV1RXYd2Ef9JIe5UXlaN+2fb1bRBjzvGRf2BnZMe4ZkWiGvZ/+hW7Y8IdXcPpA6n2fU/daon0X9uHX27/Wvl17ci3+8eBllLvo8UOfW5AgQQMN2jm2gwMcoIEGPloHjExxwL49n9d2WVlXm74zKM+ps0/cM7JjXKYjABY5wLS1Dv3vFuiKr+HQ/26pfZ8hFABgySNL4OHkAScHJ3g6eWJCrwno3r47JEjI6FaEfZEl+PUhB7RzbAd3J3c4OTrhnnQPEiQMOOcKt4o2GKT1rO2yDGHW2DIdl/Dsk5KW6VQXRnfu3FHEBWLUgpbCprkDTG0UVD/4laCs3T384FdS+z5DKHxy8hMk5CRAAw0qqiugq9LhaNFRxAyMwetDX6+90HVW9QQ89V1XdD+nR1V1Ve3j5PiVoMzlHsY+//vagQdDmBlGw+t2QbzVuH1ydXVFeXm56DIsQnVhpNfr4ejoKLoMsraWTstu7hw4U07aboUJk1/BsUl6TJj8Su37DKGggQa/3v619kJWCRJ0VTq89/17taPauiodLn/3PVzKNRhyviOcHZ1rHyfftwzlrwxG/oNltftE/zXqvxodDQd451d75eDgAKUcoqPKMHJwUN23rT4tHTra3Dlwppy0bWBGN9VcAIz0HllvPBuoGdE2hJKuSgcPJw/0GPcoPDp3RY9xj0KChDaaNtBAg6jeUfivUf9Ve5Dq2pNrax+HXZByODg4QK/Xiy7DIlQ3TccwUonQWcadlG2pr63bTYXOws6017Hv+xMYVPAAxk6dg8Hjo4x6GMMyHa7V7Bm99/17kCDB08kTI71H4tvz30L61/8A4JXpy4DpQMSuCOhu6+CgccDSR5fWBpzhvDrDW6AmBOsG4OkDqTixZxeGTX7O6DpJHhwcHFBdXS26DIvgqzKRJfR8BNA41rwFkFD0HXx/8UT17Xs4sWeX0Q9Tt2t5of8LWProUnRv3x0Lhi5A9rXs2hACULtH9Nb/vYUrt68AwH13eK27dNfUtNyJPbugK75mUp3mOn0gtckJQlI31YWRktpakpGL/wSk6pq3AGK8x6GgTwk0zg6ounPH6Bffhkt3hnPnEnISENQ1CN3bd689f87J0QkAsO/Cvtrx7obLbwuGLkD39t3h5OjU5LTcsMnPwaNzVwyb/Fyr/gqMYcvgUwMl7YEzjIgsocE+0wuRn+Dz5f+Au3tnVNwubdWLr2HpLvtaNvY/tx/hD4ZDAw3u3LuDkV+OxICOA6CBBs6OzvUucgWA9mkXMOFv7TBdG9TkPtHg8VGY/ennNlmiMzn4ZDCCL2dK2nZQxndhAiWtsZKMNDEQYXjx7dHf3+jlqYZLWQ0HDgzLdXf1d6Gr0uFGxQ10a9+t9tYRhqW4nfk78ePxw5D0etzJLbTttFwTIWJy8NlostFe6fV6xdyBQHVh5OLigjt37ogug1TC8OJ7OT/P6OUpw1LW0S+3YcMfXkH/Qrfapbqd+TtrT/Ru59gO7RzbofxuOYK6BtWe2GBYikvIScCFbreh10joPzzM2t9qfWaESKP7SeZMNqpIeXk5XF1dRZdhEaoLIyUdn0H2o+7yVEub+IbPhUaqDbC6JyQYTvTOmJ4Br3ZeKKkqQfa1bCx5ZEm9DipmYAzOhTnBd8XvMXHB/zOqTosNGJgRIo3uJxl5K3a1UtLxZqoLIyUdn0H2o+7yVGMvunVDwPC5YVNm1AZYU9cGNZy+azj8YPizsWfPWWzAwIwQseUghVIo6Xgz1V1npKRTbsk+DZv8XO11PQZ1Q8CwnzJ4fNS/f4/7byEO3H/NUFMadlam1GYrdb9fuZPLtVlKuguB6jojLtOR6Amtxjbxje0KzD1d29hTFwy15T9YxlO8myGXEXUlLdNpJKUcbGSkw4cPY+bMmTh//rzoUkiU1YE1m+uePWuWkuxIxK4I/Hr7V3Rv3x2rPN+w2k/ndZ9n/3P7LfrYSiCXzqhXr17YunUrRo0aJawGS1FdZ8Q9I7KHCa2mBgnqdjiN/XRuaud0+kAq1sx5Ea/81wSLnuKt9JMWbHltVnOUtGekujDy8PBASUmJYk66JTPIbEKrsRfuppaB6t5yvEd///uW9hrel6ilUDixZxeqS27D90dY9BRvuSxjKZler0dJSQk8PDxEl2IRqgujHj164O7duyguLhZdChEA4MTOTTUv3Ds31b6vsT0kQ7Ac/XIbdMXXcDk/776fzmPcB6B7tYQY9wE1j91CKAyb/BwcPduj4GGY1QWdPpCKT1+Zik9nTakXeEIm41R2WkNxcTHu3buHHj16iC7FIlS3ZwQAHTt2xHfffYegoCDRpZDCGbO3cHrxcJy42B7Det7G4A/+0eRjbfjDK9AVX0M7Nzc4tXNt/DEb7IdZe2/DUBMAeHTuitmffm7x5zCaHe8FmuPkyZMYP348rl+/LroUi1BdZwQAPj4+KCoqEl0GqYAxy1WDn43B7JArGPxs852JodsImzKj6f2KBvth1t7bGDb5ObRr7452bm7irw+yg71ASyoqKoKPj4/oMixGlZ3RE088gaeeegqvvvqq6FJI4eQydSUn/DuxjPXr1yMlJQWpqcoYEmFnRGRFhs4EQKODBE3tubSaBfdPLD0Zx+EGy1BaZ6TKMPL29salS5dEl0FyY8UN8KZegE/s2YWK26WoKCu772OtCipLHVTaTO3m4rE/lnHp0iV4e3uLLsNiVBlG7IyoUVa8XYFhb6Wqovy+qbOm9lyaC6oWWeqgUpgfHk2Fm1yu0bF3SuuMVHc2HcDOiJoQFlsTRFbYAG94QGpj58815Di8N+7tvYo2eg169Pc37QlDZ5l8HVVT59KZe2ZcY98vWQ47IwVgZ0SNsvLFsM12GI0sESa2O4IKp5obQV7Oz2vdkwu4BofLcdaltM5IldN0N27cQKdOnaDT6RRzlAbZuUaukdmZvxP79nyOQVpPjH3+963rLoy4BsdwzZDw64WoRTqdDp6enrhx4wa8vLxEl2MRquyMvLy84OLigosXL4ouhahGI3s8L/R/AZ+/tQ9vxH+FweOjWjfVZsQeklw6GaWfa2cJly5dgqurKzp06CC6FItRZRhpNBoEBAQgJydHdClENYxYImxpqq2pF/HTB1KxYfMxnA7872YfXy6DBRz9btkPP/yAgIAAaDQa0aVYjCrDCACCg4ORmZkpugyiWsbejrypzqW58fHG3i/XDqRHf39oHBxMH9pQkczMTAQHB4suw6JUHUZZWVmiyyCqZQiN7z7/rNGAaKlzaSqsmnq/XDuQy/l5kPT61g9tKFhWVhbDSClCQkJw8uRJ3kqCZGPY5OegcXCApNfXBoQp3UtTYdXU+y21R9RSjaZ2YML2ruzk1G+9Xo+TJ08iJCREdCkWpdowCggIwJ07d6DValv3QHbyD5haycj/zq1Z+ho8PgrjXnm13guxsd2Lsc9b9/PqhlRr6m6pRlM7MGF7V1a86NmStFotKisrERAQILoUi1JtGDk5OWHQoEGtX6qzk3/A1EpG/neu+8Jr6gt8YweIGtslGPuCb+q+kjFaqrGlj8tm78pOTv3OysrCoEGD0LZtW9GlWJRqwwiw0BCDnfwDplYy8r9z3RdeU1/gG/t8Y7sEY0PL1H0lY7RUY0sfl83elczuANwUJQ4vACq96NVgw4YN+PLLL3Hw4EHRpZACmXqrBLXeWkGt37e5xo4di2nTpiEmxvQ788qZqsPo5MmTCA8Px82bNxU1r09EyqTX6+Hl5YX09HTF3ala1ct0gYGBlhliIBKouT0Xi+3HcFBHFrRaLSoqKhQ3vACoPIwsNsRAJFBzey4W24/hoI4sZGVlYfDgwXBychJdisWpOowAnsRA9q+54QOLXbPDQR1ZUOrwAqDyPSOAQwzUiIyN/76vkcDJKsPGfo/+/ricn8cNflLs8ALAzgjBwcE8iYHqk8mSlGGJLf8fR81aapPN9TtkEYaTF5TaGak+jAIDA1FRUYG8PJ6DZY+s8oLbYElK1Iu6YYmt//Aws5baZHP9DllEfn6+Ik9eMFD9Mh0AREVFITw8HIsWLRJdCpnIFjeEs6ebztW9ZgcAr99RkFWrViE9PR179+4VXYpVqL4zAoDo6GgkJyeLLoPMYItDNUUc3GluN1a3G5LL/YnIMpKTkxEdHS26DKthZ4Sauyb26tULV69eRadOnUSXI5sNdBKnYTdm7CkFPM1AmYqLi9GtWzcUFBTA29tbdDlWwc4IgI+PD4YMGYLUVJls9MpkA53EadiNyXX/h0MStpGamoqgoCDFBhHAMKolq6U6XtOheg2X2JpbKqwbCLYOLbmGpNIofYkO4DJdrVOnTuGxxx5DcXExnJ2dRZdDZLS6S3qG08LNWaYzZ4mvtcuCXFZsWUVFBTp37oxjx45h8ODBosuxGobRv0iSBF9fX2zYsAGRkZGiyyEymqVe0BubGjTlsVua5Dt9IBVHv9wGaCSETZmBweOj7GpSUZR9+/Zhzpw5KCgoUPSBzlym+xeNRoNJkybJZ6mOyEiGJT0Ardq/GTb5ObRr746qivLaxzBlGa7u5zb2dSf27ELF7VJUlJXVvl/YLcbtiGGJTslBBDCM6omOjkZKSgpPYyC71Nr9m8Hjo+Dk4mJ2WNT93Ma+zhB27dzcat/P8fPmSZKElJQUxe8XAVymq6eyshKdO3fGkSNHMGTIENHlEJnEEst1llry416QZWRnZ2P06NH47bffFL+XzTBq4Pnnn8fAgQPx9ttviy6FqFkiXvBbek7Dx6vu3EHF7VLuBbXSihUrkJubi507d4ouxeq4TNcA943IVlp7jY6IseqWntPwcWgk7gVZQHJyMiZNmiS6DJtgGDUQFRWFU6dO4dKlS6JLIYVrbZiI2Pxv6TkNHw+bMoN7Qa106dIlnD59GlFR6vg75DJdI0aNGoWpU6di3rx5okshBZPjMhvJx7p16/DVV1/h8OHDokuxCXZGjZgyZQo2bdokugxSOGtMkrW09McTE+zHpk2bMGXKFNFl2AzDqBHTp0/H2bNnkZGRIboUIpO0FDa8rsc+ZGRkIC8vD9OnTxddis1wma4J8+fPR0VFBT7/nJNAZD/LW7aq017+PuzVyy+/DFdXV3z66aeiS7EZhlETcnJyMGzYMBQVFaFjx46iyyHBlHRsjSWCREl/H3Jz/fp1+Pj4ICMjA4GBgaLLsRku0zVh4MCBCA0NxebNm0WXQjKgpOUtS+wbKenvQ242b96MYcOGqSqIAHZGzUpMTMTSpUvx008/wcGhlbmtlBvmKeX7UDG7WmJT2b83vV6Pfv36YeXKlXjxxRdFl2NTDKNmVFVVoWfPnti6dWvrT/JeHVhzwzzPnkDsGcsUKIJSvg+yDyr797Zv3z7MnDkThYWFcHJyEl2OTXGZrhlOTk6YPXs21q1b1/oHU8oN85TyfZB9UNm/t3Xr1mH27NmqCyKAnVGLCgsL0bdvX5w7dw6+vr6iyyEihbpw4QL69esHrVaLnj17ii7H5tgZteDBBx9EVFQU4uPjRZdCMtTa8+WIDOLj4zFx4kRVBhHAMDLK/PnzkZCQgMrKStGlkMzwRAOyhMrKSiQkJGD+/PmiSxGGYWSExx9/HB4eHti9e7foUkhmOOJMlvC3v/0NHTp0wLhx40SXIgz3jIz08ccfIykpCUeOHBFdChEpTFhYGJ599lnExqpjUKMxDCMj3bhxA97e3jh+/DiCgoJEl0NkEru6tkhlsrOzMXLkSBQVFcHLy0t0OcJwmc5IHTt2xNy5c3kHWGqSnIcZuLclX8uWLcOcOXNUHUQAw8gkf/nLX5Ceno6jR4+KLoVswNRwEfWCb0yd3NuSpyNHjuD//u//sGTJEtGlCMcwMkHXrl2xaNEiLF68GFzdVD5Tw6WxF3xbdEvG1GmNeydR60iShMWLF2PRokXo0qWL6HKEYxiZ6M0338RPP/2Eb775RnQpZGWmdhONveDbolti12OfUlJScO7cOVUPLdTFAQYzrF27Fhs2bMCpU6fg6OgouhySMQ4OUGOqq6sxePBgzJ07F6+99procmSBYWSGyspK+Pv7Y8WKFZgxY4bocojIzmzZsgXvvPMO8vPzVXkOXWMYRmbatm0bli1bhvz8fDg7O4suh4jsREVFBfr374+VK1eq6rbiLeGekZl+97vfwcPDA5999pnoUojIjnz22Wfo0KEDfve734kuRVbYGbXC3r17MXPmTGi1Wnh4eIguh1SG+1H2R6fTwc/PD1u2bEFUFP+b1cXOqBWioqLg7++Pjz/+WHQppEK8kNX+fPTRRxgwYACeeOIJ0aXIDjujVjp+/DgiIyOh1WrRtWtX0eWQirAzsi9Xr15F3759sX//fgwfPlx0ObLDMLKAp556Cr6+vli7dq3oUohIpl577TVcvHgRe/bsEV2KLDGMLCA3NxchISHIzc1Fnz59RJdDRDLzyy+/ICAgAJmZmQgICBBdjiwxjCxk7ty5OH/+PNLS0qDRaESXQ0QyIUkSIiIi0KdPH8TFxYkuR7Y4wGAh//3f/42zZ88iISFBdClEJCMbNmxAfn4+Vq1aJboUWWNnZEFpaWl4/vnncebMGTz44IOiyyEiwQoKCjBw4EDs2rULERERosuRNYaRhcXExKCwsJDLdU3Y/n0B1qdrMW+MH6Y/6iu6HCKrMSzP9erVCxs2bBBdjuxxmc7CPvroI5w9e5b/+JqwPl2Lolt3sD5dK7oUIquKj49Hfn4+PvroI9Gl2AWGkYV5enoiISEBixYtQkFBgehyZGfeGD94d3DBvDF+okshspqCggL86U9/QkJCAk9nMRKX6ayEy3VkDVzmlD8uz5mHnZGVcLmOTGHsHWG5zCl/XJ4zD8PISrhcZx+2f1+AkR8cxPbvzf9v1NJjGPMcxp4zx2VOebtw4QKX58zEZTori4mJQUFBAfbv38/lOsEaW+Ia+cFBFN26A+8OLji2ONysx23pMYx5Dp4zZ/8kScL48ePRu3dvroiYgZ2RlX300UfIy8vjP04ZaGyJyxKdRkuPYcxzDB4fhdmffs4gsmPx8fH46aefuDxnJnZGNmC4GDYnJwe+vtx0FoWb/2QtFy5cwKBBg3hxayswjGwkJiYGP//8Mw4cOIC2bduKLoeILOTu3bt4/PHH0a9fP66AtAKX6WxkzZo1uH79Ov70pz+JLkX1LDG0oIQayDIWLVqEmzdvYvXq1aJLsWsMIxtxc3PDnj17sHXrVmzZskV0Oaomh/FoOdRArbd582Zs374de/bsgZubm+hy7BrDyIb8/PyQmJiI+fPnIyMjQ3Q5qmWp8ejWdDcc0bZ/J06cwB/+8AckJibyPmYWwD0jAT788EOsWbMGmZmZ6Natm+hyyAjWGgsn+3TlyhWEhIQgNjYWCxcuFF2OIjCMBJAkCdOmTUNhYSEOHjwIJycn0SVRCwzB46gBVjwViOmP+nI6T6WqqqowduxY9O7dG9u2beP1gxbCZToBNBoNEhISUF5ejtdff110OcLY0yb+vDF+cNQA1VLNfg+DSL0WLFiAiooKbNiwgUFkQQwjQVxdXZGUlIRdu3YhPj5edDlCCNvEz9gIrA6seWuk6Y/6YsVTgbX7PBxAUKe4uDjs3r0bSUlJcHFxEV2OonCZTrD09HRMnDgRBw4cwIgRI0SXY1OW6i5MfpzVgUDJRcCzJxB7plXPGezrhayCm+yQVODYsWOIiIhAamoqRo8eLbocxWEYycBf//pXvPfeezh+/Dj8/DhdZaqGgwSNhVO99zn+HTi6GgiLBUJnWfS5SZm0Wi1GjBiBZcuW4Y9//KPochSJy3Qy8Mc//hHTpk1DZGQkrl69Krocu9NwTLqxJTTD+z5My8fI73pj+/C9tUFkkxFtM5YGSR6uXr2KiIgITJ8+HX/4wx9El6NY7IxkQq/X46WXXsLZs2eRnp6u+OPnrTUAsP37AnyYlg8AWBTZHwDqLafdrryHW3fuQgPA06UtFkX2rw2q5rqbVtdrgaVBsj2dTocxY8bg4YcfxtatW+HgwJ/frYV/szLh4OCATZs2oUuXLnjmmWdQWVkpuqT7WHL6zZwBAGOef326Frfu3EV75zaY/qhv7fNkFdzEscXhWBTZH44aQAJw685dfJiWj3lj/NDBpS1uV95r8rFbPbAQFlsTRGGx5n092VxlZSWefvppdO3aFZ9//jmDyMr4tysjTk5O+Nvf/gadTocZM2ZAr9eLLqkeS06QmXMCQWPP3zCgGj5uwz8bpuLqDuROf9QX7Z3b4Nadu01+b60+MSF0Vk1H1Mo9KqvhMmI9er0eM2bMQGlpKXbt2sVrAW2Ay3Qy9Ntvv2HkyJGYMGECPvnkE9lcy2DppTVTH8+SpyA0fKwWhx6UPinHZcRakiRhwYIF2L9/P44ePYouXbqILkkVGEYydeHCBYwcORLTp0/HBx98IJtAsiRLBAkAqwWGqiblMjZabMLQnkmShLfeegs7duzA0aNH0atXL9ElqQaX6WSqV69eOHToELZu3YqlS5dCiT8zGLP01dg+Ud3luumP+uLY4vB6QWSpva2m6jPl8Zv7XFmdQCH3ZUQbkCQJS5Yswfbt23Hw4EEGkY2xM5K5vLw8jB49Gq+++ireeecdRXZITdn+fQGWf30G1RLqdSctLZ8ZOpoOLm3R3rmNxbsmUzqm5j5XVZ2XzEmShOXLlyM+Ph7p6enw9/cXXZLqsDOSOX9/fxw8eBCfffYZli1bpsgOqSnr07WolgBHDep1J411QwYLvszG5Vt34NK25p92SwMP5jB2mGH79wW4XXkPHVzaNvq5hscJ9vWST4ekQpIkYenSpYiPj8d3333HIBKEYWQHAgICkJ6ejo0bN+LPf/6zYgOpqck4wynZxtj7w2VIAKru6bEosv99odFwIs+aS2UNx8wbPpchVLMKbvKcO0EkScLixYuxadMmpKenIyAgQHRJqsUwshMDBgzA4cOHsX37dixYsEB2Y9+W0DAomuuAmjJxUA84amreNvb1xpzWYGqdTTH2uXijPTGqq6uxYMECfPHFFzh8+DA7IsG4Z2Rnzp8/j8jISAQFBWHr1q1wdnYWXZLFNLcXZM0TG0wd6Ta3FlWNistcZWUlZsyYgVOnTiEtLY3DCjLAMLJDv/32G6KiouDp6Yndu3cr/uggwLqb/Q1DwhLP1fBYIoaPfOh0Ojz99NMoLS3F3r17eR2RTHCZzg516dIFhw4dgqOjI8aMGaOKw1WtuZRVd/mspaEDUx7z1p27tUcO2cV4twpcuXIFY8aMQZs2bXDw4EEGkYwwjOyUm5sbUlJSMGDAAIwcORJarbI3vxvu/xhexBd8mW3R6biGQwctaSpMDOfddXBpi8p71bUnhjfEm/TZjlarxciRI/Hwww8jJSUFbm5uokuiOhhGdszJyQnbtm3DpEmTMGLECBw7dkx0STZjeBHf+8PlVr+Y1w06Yw5NbayOhs8//VFfnFoegVPLI+DcxrHJr+fwgm0cO3YMI0aMwFNPPYWtW7fyrDkZYhjZOQcHB3z88cdYtmwZIiIiVHMLc8OL+MRBPSz6Yl730NTlX59pMZCMCS/DiPmofl3u66KmO/4dx5wX1NzwrxFcxmu9uLg4REREYNmyZfjoo494+rZMcYBBQdLT0/H888/j2Wefxdq1a/nTnxGamqZr7OSHpgxZsR+37txFB5e2OLU8osnPa3QwooUDSnlKg/mqqqrw2muvYffu3di1axdvFS5z/BFBQcaMGYPMzEycOHEC4eHhuHLliuiSzGLLbqCxZTbDbSYsvXzW6JJcC/c5ErqMZ8e3lbhy5QrCw8ORkZGBrKwsBpEdYGekQOXl5Zg9ezYOHz6MpKQkhIaGii7JJKZ2A625LsgS1/4o9rojO72txIkTJ/DMM89gzJgxiI+Ph6urq+iSyAjsjBTI1dUV27dvR2xsLMaMGYMtW7aILskkpnYDLU2kNfdxw9CCYazbHOY+hjE3CzSWVbpJO7w77ebNmzF27FjExsZi27ZtDCI7wjBSKI1Gg4ULFyIpKQmxsbF4/fXXcffuXdFlGcXUY4BaCq+WPt5UWC34Mht+f96LBV9mt1iDOSPajdVl7qi3VUbE7ei2Enfv3sXrr7+ON998E3v27MHChQtVdcK9EnCZTgW0Wi0mT56MTp06YcuWLfD1leGSkEBNLZf1XrwXEgANgPMfTDTpMWy9dCf7JT8rKigowIwZM3Dz5k3s2bMHffr0EV0SmYFhpBJlZWWIjY3FV199hVWrVmHOnDmy/MnRGi+q5j7mgGXf4s5dPVzaOuDsu0+Y9JycgrM+SZIQHx+PP/3pT3jxxRexevVqXshqx7hMpxJubm7YsGED/vd//xfvvfceIiIiUFAgv2tXrLHcZO5jLpn4MLw7uGDJxIdNfk5ezGpdBQUFGD9+PFauXIldu3Zhw4YNDCI7xzBSmcjISJw5cwa9evXCwIEDERcXJ6v7I1njRdzcxzTnFhaW+FpqmiRJiIuLw8CBA9G7d2/k5OQgIqLpa7vIfnCZTsXS0tIQExMDf39/JCQkcC/JQoxZFlTzHo+5CgoKEBMTg/z8fCQkJDCEFIadkYo17JLi4+Nl1SXZqw/T8ps8GNWAB6Qar2431KtXL3ZDCsUwUjlPT8/avaR3331XtntJxmjpWpvmPm7rM+C4p2ScgoICRERE1Nsb8vT0FF0WWQHDiAAoo0tqzcWvluxUDAejLors3+TncE+peeyG1IdhRLUadkmRkZH45ZdfRJdltNZc/GqpTmXBl9lY/vUZBPt6MWjM9Msvv7AbUiGGEd3H0CX17t0bAQEBWLBggV3cTbalY3ma60Ys1ans/eEyqqWatwa8DYRxrl69itdeew0BAQHo06cPuyGVYRhRozw9PREXF4fMzEwUFhaib9++WL58OXQ6nejSmmXrwYCGQTNxUA84amreGj6+bM+ZFgca1Eyn0+Htt9+Gn58fLl68iMzMTMTFxbEbUhmGETUrICAAe/bsQVpaGg4dOgQ/Pz+sWbMGlZWVoktrlKl3ajVH3TPrGobf2qlB0P7nRKydGgSgJhzta+fNdiorK7FmzRr4+fkhPT0d+/fvx549exAQECC6NBKAYURGGTFiBA4fPowtW7Zg06ZN6NevH7Zs2YLq6mrRpdVT906t1uqO6i7FGbNP1cGlLTq4tG12oEFNqqursWXLFvTr1w+bNm3Cli1bcPjwYYwYMUJ0aSQQL3olk+n1euzYsQPLli1D+/bt8f7772PSpEmyOevO2heULvgyG3t/uIyJg3rUdkDUMkmSkJKSgr/85S8oLy/Hu+++i6lTp/I24ASAYUStUFlZifj4eLz77rvo168fPvjgA4SFhYkui2ToyJEjWLx4Mc6dO4dly5Zh7ty5cHJyEl0WyQh/JCGzOTs747XXXoNWq8X48ePxxBNP4Mknn0R2dsv3/yF1yM7OxpNPPomJEyciMjISWq0Wr732GoOI7sMwolZzd3fH8uXL8csvv6Bv374YOXIkwsLCsGPHDtkOOpiK49nGq6ysxI4dOxAWFoaRI0eib9++0Gq1ePvtt+Hu7i66PJIpLtORxd24cQNbtmzB+vXrUVJSgpiYGMyZM8euD2Ll/YlaVlBQgLi4OCQkJKBDhw6YP38+fv/738PLy0t0aWQH2BmRxXXs2BGxsbHIy8vD9u3b8eOPP+Khhx7CU089hbS0NOj1etElmoxnyTVOr9dj3759iI6OxkMPPYS8vDzs2LEDeXl5eOONNxhEZDR2RmQTFy9eRHx8fO1N0ObNm4eZM2eiU6dOoksjM1y/fh2bNm3CZ599hrKyMsyZMwezZ89Gz549RZdGdophRDZVVVWFpKQkrFu3DidOnMCUKVMwf/58hIaGii6NjHDixAmsW7cOX331FYYNG4b58+fj6aef5kACtRqX6cimnJyc8OKLL+Lw4cPIyMiAq6srwsPDERoaivXr1+PSpUuiS6QGLl26hHXr1iE0NBSPP/443NzckJGRgcOHD+PFF19kEJFFsDMi4UpLS7F9+3YkJibi2LFjGDx4MKKjoxEdHY0hQ4ZAo9Hwzqg2JEkSTp06heTkZCQnJ+P06dMYOXIkpkyZgunTp3MijqyCYUSycv36daSmpiIlJQXffvstOnTogOjoaBwq74nbnfzh09mD02xWUFlZiUOHDiE5ORkpKSkoKSnBhAkTEB0djSeeeIJ7e2R1DCOSrcrKShw+fBjJyclI2vM1im/cwNDhYzDv9y8iKioKnTt3Fl2iXSsuLsbevXuRkpKCtLQ0dOzYEZMmTUJ0dDRGjx4NZ2dn0SWSijCMyC5IkoTTp08jJSUFycnJyM7OxvDhwxEdHY0nn3wS/fv35xlnLdDr9cjPz8c333yD5ORk/OMf/0BQUFDtkuigQYNkc74gqQ/DiOxSUVFR7Yvqd999B2dnZwwdOhQhISEIDg5GcHAw/Pz8VBtQer0eWq0WWVlZyMzMRFZWFk6ePInKykqMGzeuNsS9vb1Fl0oEgGFEClBVVYXc3FxkZWXV/jp9+jTatWunioBqKngqKiowePDg2u89ODgYAQEBnH4jWWIYkSK1FFDBwcEICQnBoEGD4OPjAw8PD9ElG0Wn0+HSpUv44Ycf7ut4Bg0aVBs6ISEhCAgIQNu2bUWXTGQUhhGpxt27d5Gbm1v7Ip6VlYXc3FyUl5fD3d0d3t7e8PHxafJt586drdZV6fV6FBcX49KlSygqKmrybWlpKVxdXREQEMDgIUVhGJGqSZKEW7duNRkAht/fuHEDbdu2hYeHB9zc3ODu7l77y/BnV1dXODg41P4CakLG8Ku8vBylpaUoKytDaWlp7a+ysjKUlJTg3r176NixY7OB6O3tjQ4dOnDQgBSHYURkhPLycvz666/Q6XT1QqTu78vLyyFJEvR6fe3t2B0dHeHg4ACNRgNXV9d6QVb39x4eHujevTtcXV0Ff6dEYjCMiIhIOGWNFRERkV1iGBERkXAMIyIiEo5hREREwjGMiIhIOIYREREJxzAiIiLhGEZERCQcw4iIiIRjGBERkXAMIyIiEo5hREREwjGMiIhIOIYREREJxzAiIiLhGEZERCQcw4iIiIRjGBERkXAMIyIiEo5hREREwjGMiIhIOIYREREJxzAiIiLhGEZERCQcw4iIiIRjGBERkXAMIyIiEo5hREREwjGMiIhIOIYREREJxzAiIiLhGEZERCQcw4iIiIRjGBERkXAMIyIiEo5hREREwjGMiIhIOIYREREJxzAiIiLhGEZERCQcw4iIiIRjGBERkXAMIyIiEo5hREREwjGMiIhIOIYREREJxzAiIiLhGEZERCQcw4iIiIRjGBERkXAMIyIiEo5hREREwv1/BDEaX3ZGBYwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(4,4))\n",
    "ax = fig.add_subplot(111, projection='stereonet')\n",
    "\n",
    "ax.pole(strike_cluster1, dip_cluster1, marker = 'o', color = 'tab:red', markersize=1)\n",
    "ax.pole(strike_cluster2, dip_cluster2, marker = 'o', color = 'tab:blue', markersize=1)\n",
    "ax.pole(strike_background, dip_background, marker = 'o', color = 'tab:orange', markersize=1)\n",
    "ax.pole(strike_bedding1, dip_bedding1, marker = 'o', color = 'tab:green', markersize=1)\n",
    "ax.pole(strike_bedding2, dip_bedding2, marker = 'o', color = 'tab:brown', markersize=1)\n",
    "\n",
    "ax.set_azimuth_ticks([])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Frequency')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/YAAAFzCAYAAACU+pIWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAArEElEQVR4nO3deZyVBb0/8O+wjWAjGgoDymZKqJiiVFc0Ew0yxOvyurihkFpXLy6omUZ4fy5XGZcbUvIClwr1uma5dV9ikZoblMnmUoEhAiKEJc0IyiDM8/vDy4mRYWQOZ+ach+f9fr3OK8+zfp8z32/6meecM2VJkiQBAAAApFKrYhcAAAAA5E+wBwAAgBQT7AEAACDFBHsAAABIMcEeAAAAUkywBwAAgBQT7AEAACDFBHsAAABIsTbFLqC51dXVxTvvvBMVFRVRVlZW7HLImCRJ4v33349u3bpFq1Yt/3s0/U+xmQGyTP+TdWaALGvp/t/ug/0777wT3bt3L3YZZNzSpUtjjz32aPHz6n9KhRkgy/Q/WWcGyLKW6v/tPthXVFRExMcv6E477VTkasiampqa6N69e64PW5r+p9jMAFmm/8k6M0CWtXT/b/fBfuPbbnbaaScDTdEU6+1f+p9SYQbIMv1P1pkBsqyl+t+X5wEAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApFhRg/1zzz0Xxx57bHTr1i3Kysri0Ucfrbc+SZK46qqrolu3btG+ffs44ogj4vXXXy9OsQAAAFCCihrs16xZEwcccEBMmjSpwfU33nhjTJgwISZNmhR/+MMforKyMgYPHhzvv/9+C1cKAAAApalNMU/+jW98I77xjW80uC5Jkpg4cWKMGzcuTjzxxIiIuOuuu6JLly5x3333xTnnnNOSpQIAAEBJKmqwb8yiRYtixYoVMWTIkNyy8vLy+OpXvxozZszYYrCvra2N2tra3POamppmrxVKhf4n68wAWab/yTozQJaV7JfnrVixIiIiunTpUm95ly5dcusaUlVVFR07dsw9unfv3qx1bq13b5mUe0BzKdX+h5ZiBsgy/U/WmQGyrGSD/UZlZWX1nidJstmyTY0dOzaqq6tzj6VLlzZ3iVAy9D9ZZwbIMv1P1pkBsqxk34pfWVkZER/fue/atWtu+cqVKze7i7+p8vLyKC8vb/b6oBTpf7LODJBl+p+sMwNkWcnese/du3dUVlbG9OnTc8vWrVsXzz77bAwcOLCIlQEAAEDpKOod+9WrV8df/vKX3PNFixbF3Llz47Of/Wz06NEjLrroohg/fnzsvffesffee8f48eOjQ4cOcdpppxWxagAAACgdRQ32L7/8cgwaNCj3/JJLLomIiFGjRsWdd94Zl112WXz44YcxevToWLVqVXz5y1+OX//611FRUVGskgEAAKCkFDXYH3HEEZEkyRbXl5WVxVVXXRVXXXVVyxUFAAAAKVKyn7EHAAAAPp1gDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYiUd7NevXx9XXHFF9O7dO9q3bx977rlnXHPNNVFXV1fs0gAAAKAktCl2AY254YYb4tZbb4277ror9ttvv3j55ZfjzDPPjI4dO8aYMWOKXR4AAAAUXUkH+5kzZ8Zxxx0XxxxzTERE9OrVK+6///54+eWXi1wZAAAAlIaSDvaHHXZY3HrrrbFgwYLo06dPzJs3L1544YWYOHHiFvepra2N2tra3POampoWqBRKg/4n68wAWab/yTozQJaV9GfsL7/88jj11FOjb9++0bZt2+jfv39cdNFFceqpp25xn6qqqujYsWPu0b179xasuHm8e8uk3AMasz32P4V18/QFucf2yAyQZfqfrDMDZFlJB/sHH3ww7rnnnrjvvvti9uzZcdddd8V///d/x1133bXFfcaOHRvV1dW5x9KlS1uwYigu/U/WmQGyTP+TdWaALCvpt+J/97vfje9973txyimnRETE/vvvH4sXL46qqqoYNWpUg/uUl5dHeXl5S5YJJUP/k3VmgCzT/2SdGSDLSvqO/QcffBCtWtUvsXXr1v7cHQAAAPyfkr5jf+yxx8Z1110XPXr0iP322y/mzJkTEyZMiLPOOqvYpQEAAEBJKOlgf8stt8R//ud/xujRo2PlypXRrVu3OOecc+L//b//V+zSAAAAoCSUdLCvqKiIiRMnNvrn7QAAACDLSvoz9gAAAEDjBHsAAABIMcEeAAAAUkywBwAAgBQT7AEAACDFBHsAAABIMcEeAAAAUkywBwAAgBQT7AEAACDFBHsAAABIMcEeAAAAUkywBwAAgBQT7AEAACDFBHsAAABIMcEeAAAAUkywBwAAgBQT7AEAACDFBHsAAABIMcEeAAAAUkywBwAAgBQT7AEAACDFBHsAAABIsTbFLiCL3r1lUkRE7HbB+UWuBMiym6cviIiIiwf3KXIlJeKZqq3fdtDY5qsDKJjJcydHRMToA0cXuRJgezfjoXsjImLg8BFFOb879gAAAJBigj0AAACkmGAPAAAAKSbYAwAAQIoJ9gAAAJBigj0AAACkmGAPAAAAKSbYAwAAQIoJ9gAAAJBigj0AAACkmGAPAAAAKSbYAwAAQIoJ9gAAAJBigj0AAACkmGAPAAAAKSbYAwAAQIoJ9gAAAJBigj0AAACkWF7BftGiRYWuAwAAAMhDXsF+r732ikGDBsU999wTa9euLXRNAAAAwFbKK9jPmzcv+vfvH9/5zneisrIyzjnnnHjppZcKXRsAAADwKfIK9v369YsJEybEsmXLYurUqbFixYo47LDDYr/99osJEybEu+++W+g6AQAAgAZs05fntWnTJk444YT42c9+FjfccEMsXLgwLr300thjjz1i5MiRsXz58kLVCQAAADRgm4L9yy+/HKNHj46uXbvGhAkT4tJLL42FCxfG008/HcuWLYvjjjuuUHUCAAAADWiTz04TJkyIqVOnxvz582Po0KFx9913x9ChQ6NVq49/T9C7d++47bbbom/fvgUtFgAAAKgvrzv2U6ZMidNOOy2WLFkSjz76aAwbNiwX6jfq0aNH/OQnP9nmApctWxann356dOrUKTp06BAHHnhgzJo1a5uPCwAAANuDvO7Yv/HGG5+6Tbt27WLUqFH5HD5n1apVceihh8agQYNi2rRp0blz51i4cGHsvPPO23RcAAAA2F7kFeynTp0an/nMZ2L48OH1lj/00EPxwQcfbHOg3+iGG26I7t27x9SpU3PLevXqVZBjAwAAwPYgr7fiX3/99bHrrrtutrxz584xfvz4bS5qo8cffzwGDBgQw4cPj86dO0f//v3jjjvuaHSf2traqKmpqfeArND/ZJ0ZIMv0P1lnBsiyvIL94sWLo3fv3pst79mzZyxZsmSbi9rozTffjClTpsTee+8dv/rVr+Lcc8+NCy+8MO6+++4t7lNVVRUdO3bMPbp3716weqDU6f9suXn6grh5+oJG1zW2TVOOlxZmgCzbUv/f8UrjN0Vge+HfAZSCGQ/dGzMeurfFz5tXsO/cuXO88sormy2fN29edOrUaZuL2qiuri4OOuigGD9+fPTv3z/OOeec+Pa3vx1TpkzZ4j5jx46N6urq3GPp0qUFqwdKnf4n68wAWab/yTozQJbl9Rn7U045JS688MKoqKiIww8/PCIinn322RgzZkyccsopBSuua9euse+++9Zbts8++8QvfvGLLe5TXl4e5eXlBasB0kT/k3VmgCzT/2SdGSDL8gr21157bSxevDiOOuqoaNPm40PU1dXFyJEjC/oZ+0MPPTTmz59fb9mCBQuiZ8+eBTsHAAAApFlewb5du3bx4IMPxn/913/FvHnzon379rH//vsXPHBffPHFMXDgwBg/fnycdNJJ8dJLL8Xtt98et99+e0HPAwAAAGmVV7DfqE+fPtGnT59C1bKZL37xi/HII4/E2LFj45prronevXvHxIkTY8SIEc12TgAAAEiTvIL9hg0b4s4774ynnnoqVq5cGXV1dfXWP/300wUpLiJi2LBhMWzYsIIdDwAAALYneQX7MWPGxJ133hnHHHNM9OvXL8rKygpdFwAAALAV8gr2DzzwQPzsZz+LoUOHFroeAAAAoAny+jv27dq1i7322qvQtQAAAABNlFew/853vhM//OEPI0mSQtcDAAAANEFeb8V/4YUX4plnnolp06bFfvvtF23btq23/uGHHy5IcQAAAEDj8gr2O++8c5xwwgmFrgUAAABooryC/dSpUwtdBwAAAJCHvD5jHxGxfv36+M1vfhO33XZbvP/++xER8c4778Tq1asLVhwAAADQuLzu2C9evDiOPvroWLJkSdTW1sbgwYOjoqIibrzxxli7dm3ceuutha4TAAAAaEBed+zHjBkTAwYMiFWrVkX79u1zy0844YR46qmnClYcAAAA0Li8vxX/xRdfjHbt2tVb3rNnz1i2bFlBCgMAAAA+XV537Ovq6mLDhg2bLX/77bejoqJim4sCAAAAtk5ewX7w4MExceLE3POysrJYvXp1XHnllTF06NBC1QYAAAB8irzein/zzTfHoEGDYt999421a9fGaaedFm+88Ubsuuuucf/99xe6RgAAAGAL8gr23bp1i7lz58b9998fs2fPjrq6ujj77LNjxIgR9b5MDwAAAGheeQX7iIj27dvHWWedFWeddVYh6wEAAACaIK9gf/fddze6fuTIkXkVAwAAADRNXsF+zJgx9Z5/9NFH8cEHH0S7du2iQ4cOgj0AAAC0kLy+FX/VqlX1HqtXr4758+fHYYcd5svzAAAAoAXlFewbsvfee8f111+/2d18AAAAoPkULNhHRLRu3TreeeedQh4SAAAAaERen7F//PHH6z1PkiSWL18ekyZNikMPPbQghQEAAACfLq9gf/zxx9d7XlZWFrvttlsceeSR8YMf/KAQdaXOu7dMKnYJW7Sxtt0uOL/IlQDN4ebpC4q6P5/imaqt33bQ2OarAwpo8tzJEREx+sDRRa6k8Lbna4NtNeOhe3P/PHD4iG3eduM2n3YsPl1ewb6urq7QdQAAAAB5KOhn7AEAAICWldcd+0suuWSrt50wYUI+pwAAAAC2Ql7Bfs6cOTF79uxYv359fP7zn4+IiAULFkTr1q3joIMOym1XVlZWmCoBAACABuUV7I899tioqKiIu+66K3bZZZeIiFi1alWceeaZ8ZWvfCW+853vFLRIAAAAoGF5fcb+Bz/4QVRVVeVCfUTELrvsEtdee21mvxUfAAAAiiGvYF9TUxN//etfN1u+cuXKeP/997e5KAAAAGDr5BXsTzjhhDjzzDPj5z//ebz99tvx9ttvx89//vM4++yz48QTTyx0jQAAAMAW5PUZ+1tvvTUuvfTSOP300+Ojjz76+EBt2sTZZ58dN910U0ELBAAAALYsr2DfoUOHmDx5ctx0002xcOHCSJIk9tprr9hxxx0LXR8AAADQiLzeir/R8uXLY/ny5dGnT5/YcccdI0mSQtUFAAAAbIW8gv3f//73OOqoo6JPnz4xdOjQWL58eUREfOtb3/Kn7gAAAKAF5RXsL7744mjbtm0sWbIkOnTokFt+8sknx5NPPlmw4gAAAIDG5fUZ+1//+tfxq1/9KvbYY496y/fee+9YvHhxQQoDAAAAPl1ed+zXrFlT7079Rn/729+ivLx8m4sCAAAAtk5ewf7www+Pu+++O/e8rKws6urq4qabbopBgwYVrDgAAACgcXm9Ff+mm26KI444Il5++eVYt25dXHbZZfH666/He++9Fy+++GKhawQAAAC2IK879vvuu2+88sor8aUvfSkGDx4ca9asiRNPPDHmzJkTn/vc5wpdIwAAALAFTb5j/9FHH8WQIUPitttui6uvvro5agIAAAC2UpPv2Ldt2zZee+21KCsra456AAAAgCbI6634I0eOjJ/85CeFrgUAAABoory+PG/dunXx4x//OKZPnx4DBgyIHXfcsd76CRMmFKQ4AAAAoHFNCvZvvvlm9OrVK1577bU46KCDIiJiwYIF9bbxFn0AAABoOU0K9nvvvXcsX748nnnmmYiIOPnkk+NHP/pRdOnSpVmKAwAAABrXpM/YJ0lS7/m0adNizZo1BS2oMVVVVVFWVhYXXXRRi50TAAAASlleX5630SeDfnP6wx/+ELfffnt84QtfaLFzAgAAQKlrUrAvKyvb7DP0LfGZ+tWrV8eIESPijjvuiF122aXZzwcAAABp0aTP2CdJEt/85jejvLw8IiLWrl0b55577mbfiv/www8XrsKIOO+88+KYY46Jr33ta3Httdc2um1tbW3U1tbmntfU1BS0Fihl+p+sMwNkmf4n68wAWdakYD9q1Kh6z08//fSCFtOQBx54IGbPnh1/+MMftmr7qqqquPrqqz91u3dvmRQREbtdcP5WLdt0eUPLSkVD10B2bG3/03xunv7xXwq5eHCfgh5n4/N8j73p/k3Zdluvo6W12Aw8U9X850iD5nodBo1tnuOmrd4m2tr+nzx3ckREjD5wdIPPG9p2o4a2aep58j1ePuckW/x3UOHNeOjeFtmnUDaee+DwEdu0zdbsk89xmlOTgv3UqVObq44GLV26NMaMGRO//vWvY4cddtiqfcaOHRuXXHJJ7nlNTU107969uUqEkqL/yTozQJbpf7LODJBlTQr2LW3WrFmxcuXKOPjgg3PLNmzYEM8991xMmjQpamtro3Xr1vX2KS8vz31UALJG/5N1ZoAs0/9knRkgy0o62B911FHx6quv1lt25plnRt++fePyyy/fLNQDAABA1pR0sK+oqIh+/frVW7bjjjtGp06dNlsOAAAAWbRNf8ceAAAAKK6SvmPfkN/+9rfFLgEAAABKhjv2AAAAkGKCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAAp1qbYBaTJu7dMavbj7XbB+c12PmD7cPP0BSV1zosH92nBSsjLM1Vbv+2gsc1XR3NoyrVl2OS5k5t929EHjm7y/ls658Zj5euTx2uopm09B1DfjIfuzWvdtpxr4PAR23yMTW08XqHrbQnu2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAACkm2AMAAECKCfYAAACQYoI9AAAApJhgDwAAAClW0sG+qqoqvvjFL0ZFRUV07tw5jj/++Jg/f36xywIAAICSUdLB/tlnn43zzjsvfve738X06dNj/fr1MWTIkFizZk2xSwMAAICS0KbYBTTmySefrPd86tSp0blz55g1a1YcfvjhRaoKAAAASkdJB/tPqq6ujoiIz372s1vcpra2Nmpra3PPa2pqmr0uKBX6n6wzA2SZ/ifrzABZlppgnyRJXHLJJXHYYYdFv379trhdVVVVXH311Q2ue/eWSXkvKxUN1bbbBedv1XZbWt/Q/tuiOY9N4xrrf0rHzdMXRETExYP71Hu+6bJSsWltaZCpGXimqtgVsNGWfhZr1rZoGaXa/5PnTm7yuk2Xjz5wdJOP25S6Nh5/W49H8ZXqDJCfGQ/dW+wSUqWkP2O/qfPPPz9eeeWVuP/++xvdbuzYsVFdXZ17LF26tIUqhOLT/2SdGSDL9D9ZZwbIslTcsb/gggvi8ccfj+eeey722GOPRrctLy+P8vLyFqoMSov+J+vMAFmm/8k6M0CWlXSwT5IkLrjggnjkkUfit7/9bfTu3bvYJQEAAEBJKelgf95558V9990Xjz32WFRUVMSKFSsiIqJjx47Rvn37IlcHAAAAxVfSn7GfMmVKVFdXxxFHHBFdu3bNPR588MFilwYAAAAloaTv2CdJUuwSAAAAoKSV9B17AAAAoHGCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKRYm2IXkHbv3jKp2CVscw0b99/tgvO36pibblfIOqCU3Dx9QUREXDy4T5PW5XOOxpY1tE2p2VKN2/r6FNRzP4jYcYdiV/HpnqkqdgVqyLDJcycXdf+tPf7oA0c3y3GBppnx0L1FP+6Wts23tk/br6H1zfU6NJU79gAAAJBigj0AAACkmGAPAAAAKSbYAwAAQIoJ9gAAAJBigj0AAACkmGAPAAAAKSbYAwAAQIoJ9gAAAJBigj0AAACkmGAPAAAAKSbYAwAAQIoJ9gAAAJBigj0AAACkmGAPAAAAKSbYAwAAQIoJ9gAAAJBigj0AAACkmGAPAAAAKSbYAwAAQIoJ9gAAAJBigj0AAACkmGAPAAAAKSbYAwAAQIoJ9gAAAJBiqQj2kydPjt69e8cOO+wQBx98cDz//PPFLgkAAABKQskH+wcffDAuuuiiGDduXMyZMye+8pWvxDe+8Y1YsmRJsUsDAACAoiv5YD9hwoQ4++yz41vf+lbss88+MXHixOjevXtMmTKl2KUBAABA0bUpdgGNWbduXcyaNSu+973v1Vs+ZMiQmDFjRoP71NbWRm1tbe55dXV1RETU1NTEug8/bL5iC6S8pib3z+8XuN5PO/bWnnvT7Tb1acfMopr/u/4kSVrkfI31P02zds3qiGj4tct33Se32d7V1NSUzgx8ULulXaDZbOy7Yvf/2jVrG9x+43x+uHrr/3sjn3221SfP2VI1+HfntiuZfwf4WeZtzQcftMh5Nv0ZtdQ5m1tL938kJWzZsmVJRCQvvvhiveXXXXdd0qdPnwb3ufLKK5OI8PAoqcfChQtbYmT0v0fJPsyAR5Yf+t8j6w8z4JHlR0v1f1mStNSvEJrunXfeid133z1mzJgRhxxySG75ddddF//zP/8Tf/7znzfb55O/qaurq4v33nsvOnXqFGVlZS1S96Zqamqie/fusXTp0thpp51a/PzFlvXrr66ujh49esSqVati5513bvbzlVr/R+iBrF9/1mcg6z//rF+//s/2zz/r1x9hBrLeA1m//pbu/5J+K/6uu+4arVu3jhUrVtRbvnLlyujSpUuD+5SXl0d5eXm9ZS3xQn6anXbaKZMNvVHWr79Vq5b5OotS7f8IPZD168/6DGT955/169f/2f75Z/36I8xA1nsg69ffUv1f0l+e165duzj44INj+vTp9ZZPnz49Bg4cWKSqAAAAoHSU9B37iIhLLrkkzjjjjBgwYEAccsghcfvtt8eSJUvi3HPPLXZpAAAAUHQlH+xPPvnk+Pvf/x7XXHNNLF++PPr16xdPPPFE9OzZs9ilbZXy8vK48sorN3tbUFa4/mxff4TXwPW7ftfv+l2/68+qrL8Grt/1t+T1l/SX5wEAAACNK+nP2AMAAACNE+wBAAAgxQR7AAAASDHBHgAAAFJMsM/TsmXL4vTTT49OnTpFhw4d4sADD4xZs2bl1idJEldddVV069Yt2rdvH0cccUS8/vrr9Y5RW1sbF1xwQey6666x4447xr/+67/G22+/3dKX0mTr16+PK664Inr37h3t27ePPffcM6655pqoq6vLbbM9Xf9zzz0Xxx57bHTr1i3Kysri0Ucfrbe+UNe6atWqOOOMM6Jjx47RsWPHOOOMM+If//hHM19dfvR/dvo/wgw0xAxkZwb0/+b0v/7fSP9nq/8jzEBJz0BCk7333ntJz549k29+85vJ73//+2TRokXJb37zm+Qvf/lLbpvrr78+qaioSH7xi18kr776anLyyScnXbt2TWpqanLbnHvuucnuu++eTJ8+PZk9e3YyaNCg5IADDkjWr19fjMvaatdee23SqVOn5H//93+TRYsWJQ899FDymc98Jpk4cWJum+3p+p944olk3LhxyS9+8YskIpJHHnmk3vpCXevRRx+d9OvXL5kxY0YyY8aMpF+/fsmwYcNa6jK3mv7PVv8niRn4JDOQrRnQ//Xpf/2/Kf2frf5PEjNQyjMg2Ofh8ssvTw477LAtrq+rq0sqKyuT66+/Prds7dq1SceOHZNbb701SZIk+cc//pG0bds2eeCBB3LbLFu2LGnVqlXy5JNPNl/xBXDMMcckZ511Vr1lJ554YnL66acnSbJ9X/8nB7pQ1/rHP/4xiYjkd7/7XW6bmTNnJhGR/PnPf27mq2oa/Z/d/k8SM5AkZiDLM6D/9b/+fyT3XP9vbnv++W9kBh7JPS+1GfBW/Dw8/vjjMWDAgBg+fHh07tw5+vfvH3fccUdu/aJFi2LFihUxZMiQ3LLy8vL46le/GjNmzIiIiFmzZsVHH31Ub5tu3bpFv379ctuUqsMOOyyeeuqpWLBgQUREzJs3L1544YUYOnRoRGz/17+pQl3rzJkzo2PHjvHlL385t82//Mu/RMeOHUvu9dD/+n9TZsAMZHkG9L/+1//6P0v9H2EGNlVqM9BmWy8oi958882YMmVKXHLJJfH9738/XnrppbjwwgujvLw8Ro4cGStWrIiIiC5dutTbr0uXLrF48eKIiFixYkW0a9cudtlll8222bh/qbr88sujuro6+vbtG61bt44NGzbEddddF6eeempExHZ//Zsq1LWuWLEiOnfuvNnxO3fuXHKvh/7X/5syA2YgyzOg//W//tf/Wer/CDOwqVKbAcE+D3V1dTFgwIAYP358RET0798/Xn/99ZgyZUqMHDkyt11ZWVm9/ZIk2WzZJ23NNsX24IMPxj333BP33Xdf7LfffjF37ty46KKLolu3bjFq1Kjcdtvr9TekENfa0Pal+Hrof/3fEDNgBrI8A/pf/+v/f9L/2/fP3wxsrlRmwFvx89C1a9fYd9996y3bZ599YsmSJRERUVlZGRGx2W9YVq5cmfuNTmVlZaxbty5WrVq1xW1K1Xe/+9343ve+F6ecckrsv//+ccYZZ8TFF18cVVVVEbH9X/+mCnWtlZWV8de//nWz47/77rsl93rof/2/KTPwMTOQzRnQ/x/T//p/U/p/+/75m4F/KrUZEOzzcOihh8b8+fPrLVuwYEH07NkzIiJ69+4dlZWVMX369Nz6devWxbPPPhsDBw6MiIiDDz442rZtW2+b5cuXx2uvvZbbplR98MEH0apV/dZp3bp17s9cbO/Xv6lCXeshhxwS1dXV8dJLL+W2+f3vfx/V1dUl93rof/2/KTPwMTOQzRnQ/x/T//p/I/2//f/8zcA/ldwMbPXX7JHz0ksvJW3atEmuu+665I033kjuvffepEOHDsk999yT2+b6669POnbsmDz88MPJq6++mpx66qkN/umDPfbYI/nNb36TzJ49OznyyCNL8s88fNKoUaOS3XffPfdnLh5++OFk1113TS677LLcNtvT9b///vvJnDlzkjlz5iQRkUyYMCGZM2dOsnjx4iRJCnetRx99dPKFL3whmTlzZjJz5sxk//33L8k/9aL/s9X/SWIGPskMZGsG9H99+l//6//s9n+SmIFSngHBPk+//OUvk379+iXl5eVJ3759k9tvv73e+rq6uuTKK69MKisrk/Ly8uTwww9PXn311XrbfPjhh8n555+ffPazn03at2+fDBs2LFmyZElLXkZeampqkjFjxiQ9evRIdthhh2TPPfdMxo0bl9TW1ua22Z6u/5lnnkkiYrPHqFGjkiQp3LX+/e9/T0aMGJFUVFQkFRUVyYgRI5JVq1a10FU2jf7PTv8niRloiBnIzgzo/83pf/2v/7PZ/0liBkp5BsqSJEm2/v4+AAAAUEp8xh4AAABSTLAHAACAFBPsAQAAIMUEewAAAEgxwR4AAABSTLAHAACAFBPsAQAAIMUEexr11ltvRVlZWcydO7fJ+z799NPRt2/fqKury/v8tbW10aNHj5g1a1bex4B86X+yTP+TdWaALNP/6SPYl7BvfvObUVZWFmVlZdG2bdvo0qVLDB48OH76059u06A0dr7jjz++YMe77LLLYty4cdGqVf5tVl5eHpdeemlcfvnlBauLdND/+j/L9L/+zzozYAayTP/r/3wI9iXu6KOPjuXLl8dbb70V06ZNi0GDBsWYMWNi2LBhsX79+mKXt0UzZsyIN954I4YPH77NxxoxYkQ8//zz8ac//akAlZEm+l//Z5n+1/9ZZwbMQJbpf/3fVIJ9iSsvL4/KysrYfffd46CDDorvf//78dhjj8W0adPizjvvzG1XXV0d//7v/x6dO3eOnXbaKY488siYN29ebv1VV10VBx54YNx2223RvXv36NChQwwfPjz+8Y9/5Nbfdddd8dhjj+V+Q/jb3/42t/+bb74ZgwYNig4dOsQBBxwQM2fObLTuBx54IIYMGRI77LDDZjX89Kc/jR49esRnPvOZ+I//+I/YsGFD3HjjjVFZWRmdO3eO6667rt6xOnXqFAMHDoz7778//xeSVNL/+j/L9L/+zzozYAayTP/r/6YS7FPoyCOPjAMOOCAefvjhiIhIkiSOOeaYWLFiRTzxxBMxa9asOOigg+Koo46K9957L7ffX/7yl/jZz34Wv/zlL+PJJ5+MuXPnxnnnnRcREZdeemmcdNJJud8OLl++PAYOHJjbd9y4cXHppZfG3Llzo0+fPnHqqac2+tvC5557LgYMGLDZ8oULF8a0adPiySefjPvvvz9++tOfxjHHHBNvv/12PPvss3HDDTfEFVdcEb/73e/q7felL30pnn/++W163dg+6H+yTP+TdWaALNP/NCqhZI0aNSo57rjjGlx38sknJ/vss0+SJEny1FNPJTvttFOydu3aett87nOfS2677bYkSZLkyiuvTFq3bp0sXbo0t37atGlJq1atkuXLl2/xfIsWLUoiIvnxj3+cW/b6668nEZH86U9/2mLtHTt2TO6+++56y6688sqkQ4cOSU1NTW7Z17/+9aRXr17Jhg0bcss+//nPJ1VVVfX2/eEPf5j06tVri+dj+6P//0n/Z4/+/yf9n01m4J/MQPbo/3/S/1uvTbF+ocC2SZIkysrKIiJi1qxZsXr16ujUqVO9bT788MNYuHBh7nmPHj1ijz32yD0/5JBDoq6uLubPnx+VlZWNnu8LX/hC7p+7du0aERErV66Mvn37Nrj9hx9+WO8tOBv16tUrKioqcs+7dOkSrVu3rvflGl26dImVK1fW2699+/bxwQcfNFoj2aH/yTL9T9aZAbJM/7Mlgn1K/elPf4revXtHRERdXV107dq13udhNtp55523eIyN/6ew8X8b07Zt2832a+xbOXfddddYtWpVo8fZeKyGln3y2O+9917stttun1on2aD/yTL9T9aZAbJM/7Mlgn0KPf300/Hqq6/GxRdfHBERBx10UKxYsSLatGkTvXr12uJ+S5YsiXfeeSe6desWEREzZ86MVq1aRZ8+fSIiol27drFhw4aC1Ni/f//44x//WJBjRUS89tpr0b9//4Idj/TS/2SZ/ifrzABZpv9pjC/PK3G1tbWxYsWKWLZsWcyePTvGjx8fxx13XAwbNixGjhwZERFf+9rX4pBDDonjjz8+fvWrX8Vbb70VM2bMiCuuuCJefvnl3LF22GGHGDVqVMybNy+ef/75uPDCC+Okk07KvQWnV69e8corr8T8+fPjb3/7W3z00Ud51/31r389XnjhhW27+E08//zzMWTIkIIdj3TQ/x/T/9mk/z+m/7PLDHzMDGST/v+Y/t96gn2Je/LJJ6Nr167Rq1evOProo+OZZ56JH/3oR/HYY49F69atI+Ljt6088cQTcfjhh8dZZ50Vffr0iVNOOSXeeuut6NKlS+5Ye+21V5x44okxdOjQGDJkSPTr1y8mT56cW//tb387Pv/5z8eAAQNit912ixdffDHvuk8//fT44x//GPPnz8//4v/PzJkzo7q6Ov7t3/5tm49Fuuh//Z9l+l//Z50ZMANZpv/1f1OVJUmSFLsImt9VV10Vjz76aMydO7fFznnZZZdFdXV13Hbbbdt0nOHDh0f//v3j+9//foEqI2v0P1mm/8k6M0CW6f/scMeeZjNu3Ljo2bPnNn1mp7a2Ng444IDcZ4kgLfQ/Wab/yTozQJbp/+Jwxz4jivHbOigV+p8s0/9knRkgy/R/dgj2AAAAkGLeig8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAApJtgDAABAign2AAAAkGKCPQAAAKSYYA8AAAAp9v8Bkkq6M9IWhgcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1200x400 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make a set of five histograms using plt.subplots and defining explicitly the axes\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=5, figsize=(12, 4), sharey=True)\n",
    "\n",
    "axes[0].hist(depth_cluster1, bins=36, color = 'tab:red', alpha = 0.5)\n",
    "axes[1].hist(depth_cluster2, bins=36, color = 'tab:blue', alpha = 0.5)\n",
    "axes[2].hist(depth_background, bins=36, color = 'tab:orange', alpha = 0.5)\n",
    "axes[3].hist(depth_bedding1, bins=36, color = 'tab:green', alpha = 0.5)\n",
    "axes[4].hist(depth_bedding2, bins=36, color = 'tab:brown', alpha = 0.5)\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_xlabel('Depth (m)')\n",
    "    ax.set_xlim(500, 1000)\n",
    "\n",
    "axes[0].set_ylabel('Frequency')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Frequency')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+0AAAF4CAYAAAAhcxjKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsOElEQVR4nO3dfZRVdb0/8M8RZBxwIBWYgeSpBC1JF0qppIIWlHZdInWXpSU+1NVQgsxFGbfFVAaGq5GKm9q9RHhviNXSrrcHLlMiGegvRAilVPKSoDLiw8STMAjs3x9cz2WYQZnDzJwvnNdrrb2c89377P2ZfT5f6j37nH1yWZZlAQAAACTniGIXAAAAADRPaAcAAIBECe0AAACQKKEdAAAAEiW0AwAAQKKEdgAAAEiU0A4AAACJEtoBAAAgUUI7AAAAJEpoBwAAgEQVNbRXV1dHLpdrtFRVVeXXZ1kW1dXV0bt37ygvL48RI0bEqlWrilgxAAAAtJ+iX2k/+eSTY/369fnliSeeyK+bPn161NTUxMyZM2Pp0qVRVVUVI0eOjM2bNxexYgAAAGgfRQ/tHTt2jKqqqvzSo0ePiNhzlX3GjBkxefLkGDNmTAwePDjmzJkTr7/+esydO7fIVQMAAEDb61jsAlavXh29e/eOsrKyOOOMM2Lq1Knxrne9K9asWRN1dXUxatSo/LZlZWUxfPjwWLJkSVx77bXN7q+hoSEaGhryj3fv3h2vvfZaHHfccZHL5dr894G9ZVkWmzdvjt69e8cRR7T938j0P6kxByhl+p9SZw5Qylq1/7Mi+vWvf539/Oc/z1auXJnV1tZmw4cPzyorK7NXXnklW7x4cRYR2QsvvNDoOZ/73OeyUaNG7XefU6ZMySLCYklqWbduXVtPJ/1vSXoxByylvOh/S6kv5oCllJfW6P9clmVZJGLr1q3x7ne/OyZNmhRnnnlmfPCDH4wXX3wxevXqld/mc5/7XKxbty7mz5/f7D72/Qvbxo0bo2/fvrFu3bro2rVrm/8OsLdNmzZFnz594u9//3t069atzY+n/0mNOUApS73//3Xlv0ZExOdO+Vyb10ZpSn0O/L/7740zLrm0zeuiNLVm/xf97fF769KlS7zvfe+L1atXx+jRoyMioq6urlFo37BhQ1RWVu53H2VlZVFWVtZkvGvXrv4PG0XTXm/J0v+kyhyglKXa/+VHl+fXQ1tKdQ506dxZ/9PmWqP/i34jur01NDTEX/7yl+jVq1cMGDAgqqqqora2Nr9+x44dsWjRohg2bFgRqwQAAID2UdQr7TfddFNcdNFF0bdv39iwYUPccsstsWnTphg7dmzkcrmYOHFiTJ06NQYOHBgDBw6MqVOnRufOneOyyy4rZtkAAADQLooa2p9//vn41Kc+Fa+88kr06NEjzjzzzHj00UejX79+ERExadKk2LZtW4wbNy7q6+vjjDPOiAULFkRFRUUxywYAAIB2UdTQPm/evLdcn8vlorq6Oqqrq9unIAAAAEhIUp9pBwAAAP6P0A4AAACJEtoBAAAgUUI7AAAAJEpoBwAAgEQJ7QAAAJAooR0AAAASJbQDAABAooR2AAAASJTQDgAAAIkS2gEAACBRQjsAAAAkSmgHAACARAntAAAAkCihHQAAABIltAMAAECihHYAAABIlNAOAAAAiRLaAQAAIFFCOwAAACRKaAcAAIBECe0AAACQKKEdAAAAEiW0AwAAQKKEdgAAAEiU0A4AAACJEtoBAAAgUUI7AAAAJEpoBwAAgEQJ7QAAAJAooR0AAAASJbQDAABAooR2AAAASJTQDgAAAIkS2gEAACBRQjsAAAAkSmgHAACARAntAAAAkCihHQAAABIltAMAAECihHYAAABIlNAOAAAAiRLaAQAAIFFCOwAAACRKaAcAAIBECe0AAACQKKEdAAAAEiW0AwAAQKKEdgAAAEiU0A4AAACJEtoBAAAgUcmE9mnTpkUul4uJEyfmx7Isi+rq6ujdu3eUl5fHiBEjYtWqVcUrEgAAANpREqF96dKl8cMf/jBOOeWURuPTp0+PmpqamDlzZixdujSqqqpi5MiRsXnz5iJVCgAAAO2n6KF9y5Ytcfnll8e//uu/xjHHHJMfz7IsZsyYEZMnT44xY8bE4MGDY86cOfH666/H3Llz97u/hoaG2LRpU6MFSoX+p9SZA5Qy/U+pMwc4XBU9tF9//fXxsY99LD784Q83Gl+zZk3U1dXFqFGj8mNlZWUxfPjwWLJkyX73N23atOjWrVt+6dOnT5vVDqnR/5Q6c4BSpv8pdeYAh6uihvZ58+bF448/HtOmTWuyrq6uLiIiKisrG41XVlbm1zXn5ptvjo0bN+aXdevWtW7RkDD9T6kzByhl+p9SZw5wuOpYrAOvW7cuJkyYEAsWLIijjjpqv9vlcrlGj7MsazK2t7KysigrK2u1OuFQov8pdeYApUz/U+rMAQ5XRbvSvmzZstiwYUOcfvrp0bFjx+jYsWMsWrQovve970XHjh3zV9j3vaq+YcOGJlffAQAA4HBUtND+oQ99KJ544olYsWJFfhk6dGhcfvnlsWLFinjXu94VVVVVUVtbm3/Ojh07YtGiRTFs2LBilQ0AAADtpmhvj6+oqIjBgwc3GuvSpUscd9xx+fGJEyfG1KlTY+DAgTFw4MCYOnVqdO7cOS677LJilAwAAADtqmih/UBMmjQptm3bFuPGjYv6+vo444wzYsGCBVFRUVHs0gAAAKDNJRXaH3rooUaPc7lcVFdXR3V1dVHqAQAAgGIq+ve0AwAAAM0T2gEAACBRQjsAAAAkSmgHAACARAntAAAAkCihHQAAABIltAMAAECihHYAAABIlNAOAAAAiRLaAQAAIFFCOwAAACRKaAcAAIBECe0AAACQKKEdAAAAEiW0AwAAQKKEdgAAAEiU0A4AAACJEtoBAAAgUUI7AAAAJEpoBwAAgEQJ7QAAAJAooR0AAAASJbQDAABAooR2KHG31z4Tt9c+U+wyAACAZgjtAAAAkCihHQAAABIltAMAAECihHYAAABIlNAOAAAAiRLaAQAAIFFCOwAAACRKaAcAAIBECe0AAACQKKEdAAAAEiW0AwAAQKKEdgAAAEiU0A4AAACJEtoBAAAgUUI7AAAAJEpoBwAAgEQJ7QAAAJAooR0AAAASJbQDAABAooR2AAAASJTQDgAAAIkS2gEAACBRQjsAAAAkSmgHAACARAntAAAAkCihHQAAABJV1NB+xx13xCmnnBJdu3aNrl27xllnnRW/+c1v8uuzLIvq6uro3bt3lJeXx4gRI2LVqlVFrBgAAADaT1FD+/HHHx+33nprPPbYY/HYY4/F+eefHxdffHE+mE+fPj1qampi5syZsXTp0qiqqoqRI0fG5s2bi1k2AAAAtIuCQvuaNWta5eAXXXRRXHjhhTFo0KAYNGhQfOtb34qjjz46Hn300ciyLGbMmBGTJ0+OMWPGxODBg2POnDnx+uuvx9y5c1vl+AAAAJCygkL7CSecEOedd178x3/8R2zfvr1VCtm1a1fMmzcvtm7dGmeddVasWbMm6urqYtSoUfltysrKYvjw4bFkyZL97qehoSE2bdrUaIFSof8pdeYApUz/U+rMAQ5XBYX2P/3pTzFkyJD40pe+FFVVVXHttdfGH//4x4IKeOKJJ+Loo4+OsrKyuO666+L++++P9773vVFXVxcREZWVlY22r6yszK9rzrRp06Jbt275pU+fPgXVBYci/U+pMwcoZfqfUmcOcLgqKLQPHjw4ampq4oUXXojZs2dHXV1dnH322XHyySdHTU1NvPzyywe8rxNPPDFWrFgRjz76aHz+85+PsWPHxp///Of8+lwu12j7LMuajO3t5ptvjo0bN+aXdevWtfwXhEOU/qfUmQOUMv1PqTMHOFwd1I3oOnbsGJdcckn89Kc/jW9/+9vx7LPPxk033RTHH398XHHFFbF+/fq33UenTp3ihBNOiKFDh8a0adPi1FNPje9+97tRVVUVEdHkqvqGDRuaXH3fW1lZWf5u9G8uUCr0P6XOHKCU6X9KnTnA4eqgQvtjjz0W48aNi169ekVNTU3cdNNN8eyzz8aDDz4YL7zwQlx88cUt3meWZdHQ0BADBgyIqqqqqK2tza/bsWNHLFq0KIYNG3YwZQMAAMAhoWMhT6qpqYnZs2fH008/HRdeeGHcfffdceGFF8YRR+z5G8CAAQPirrvuipNOOukt9/PVr341LrjggujTp09s3rw55s2bFw899FDMnz8/crlcTJw4MaZOnRoDBw6MgQMHxtSpU6Nz585x2WWXFVI2AAAAHFIKCu133HFHXH311XHVVVfl38a+r759+8asWbPecj8vvfRSfOYzn4n169dHt27d4pRTTon58+fHyJEjIyJi0qRJsW3bthg3blzU19fHGWecEQsWLIiKiopCygYAAIBDSkGhffXq1W+7TadOnWLs2LFvuc3bhfpcLhfV1dVRXV3dkvIAAADgsFDQZ9pnz54dP/vZz5qM/+xnP4s5c+YcdFEAAABAgaH91ltvje7duzcZ79mzZ0ydOvWgiwIAAAAKDO3PPfdcDBgwoMl4v379Yu3atQddFAAAAFBgaO/Zs2esXLmyyfif/vSnOO644w66KAAAAKDA0P7JT34yvvCFL8TChQtj165dsWvXrnjwwQdjwoQJ8clPfrK1awQAAICSVNDd42+55ZZ47rnn4kMf+lB07LhnF7t3744rrrjCZ9oBAACglRQU2jt16hT33ntvfPOb34w//elPUV5eHu973/uiX79+rV0fAAAAlKyCQvubBg0aFIMGDWqtWgAAAIC9FBTad+3aFT/+8Y/jd7/7XWzYsCF2797daP2DDz7YKsUBAABAKSsotE+YMCF+/OMfx8c+9rEYPHhw5HK51q4LAAAASl5BoX3evHnx05/+NC688MLWrgcAAAD4XwV95VunTp3ihBNOaO1aAAAAgL0UFNq/9KUvxXe/+93Isqy16wEAAAD+V0Fvj//DH/4QCxcujN/85jdx8sknx5FHHtlo/X333dcqxQEAAEApKyi0v+Md74hLLrmktWsBAAAA9lJQaJ89e3Zr1wEAAADso6DPtEdE7Ny5M37729/GXXfdFZs3b46IiBdffDG2bNnSasUBAABAKSvoSvtzzz0XH/3oR2Pt2rXR0NAQI0eOjIqKipg+fXps37497rzzztauEwAAAEpOQVfaJ0yYEEOHDo36+vooLy/Pj19yySXxu9/9rtWKAwAAgFJW8N3jFy9eHJ06dWo03q9fv3jhhRdapTAAAAAodQVdad+9e3fs2rWryfjzzz8fFRUVB10UAAAAUGBoHzlyZMyYMSP/OJfLxZYtW2LKlClx4YUXtlZtAAAAUNIKenv87bffHuedd168973vje3bt8dll10Wq1evju7du8c999zT2jUCAABASSootPfu3TtWrFgR99xzTzz++OOxe/fuuOaaa+Lyyy9vdGM6AAAAoHAFhfaIiPLy8rj66qvj6quvbs16AAAAgP9VUGi/++6733L9FVdcUVAxAAAAwP8pKLRPmDCh0eM33ngjXn/99ejUqVN07txZaAcAAIBWUNDd4+vr6xstW7ZsiaeffjrOPvtsN6IDAACAVlJQaG/OwIED49Zbb21yFR4AAAAoTKuF9oiIDh06xIsvvtiauwQAAICSVdBn2h944IFGj7Msi/Xr18fMmTPjgx/8YKsUBgAAAKWuoNA+evToRo9zuVz06NEjzj///PjOd77TGnUBAABAySsotO/evbu16wAAAAD20aqfaQcAAABaT0FX2m+88cYD3rampqaQQwAAAEDJKyi0L1++PB5//PHYuXNnnHjiiRER8cwzz0SHDh3itNNOy2+Xy+Vap0oAAAAoQQWF9osuuigqKipizpw5ccwxx0RERH19fVx11VVxzjnnxJe+9KVWLRIAAABKUUGfaf/Od74T06ZNywf2iIhjjjkmbrnlFnePBwAAgFZSUGjftGlTvPTSS03GN2zYEJs3bz7oogAAAIACQ/sll1wSV111Vfz85z+P559/Pp5//vn4+c9/Htdcc02MGTOmtWsEAACAklTQZ9rvvPPOuOmmm+LTn/50vPHGG3t21LFjXHPNNXHbbbe1aoEAAABQqgoK7Z07d44f/OAHcdttt8Wzzz4bWZbFCSecEF26dGnt+gAAAKBkFfT2+DetX78+1q9fH4MGDYouXbpElmWtVRcAtI+F04pdAQDAfhUU2l999dX40Ic+FIMGDYoLL7ww1q9fHxERn/3sZ33dGwAAALSSgkL7F7/4xTjyyCNj7dq10blz5/z4pZdeGvPnz2+14gAAAKCUFfSZ9gULFsR///d/x/HHH99ofODAgfHcc8+1SmEAAABQ6gq60r5169ZGV9jf9Morr0RZWdlBFwUAAAAUGNrPPffcuPvuu/OPc7lc7N69O2677bY477zzWq04AAAAKGUFvT3+tttuixEjRsRjjz0WO3bsiEmTJsWqVavitddei8WLF7d2jW3u5e/PbPS4x/gbilQJAAAA/J+CrrS/973vjZUrV8YHPvCBGDlyZGzdujXGjBkTy5cvj3e/+92tXSMAAACUpBZfaX/jjTdi1KhRcdddd8XXv/71tqgJAAAAiAKutB955JHx5JNPRi6XO+iDT5s2Ld7//vdHRUVF9OzZM0aPHh1PP/10o22yLIvq6uro3bt3lJeXx4gRI2LVqlUHfWwAAABIXUFvj7/iiiti1qxZB33wRYsWxfXXXx+PPvpo1NbWxs6dO2PUqFGxdevW/DbTp0+PmpqamDlzZixdujSqqqpi5MiRsXnz5oM+PgAAAKSsoBvR7dixI/7t3/4tamtrY+jQodGlS5dG62tqag5oP/Pnz2/0ePbs2dGzZ89YtmxZnHvuuZFlWcyYMSMmT54cY8aMiYiIOXPmRGVlZcydOzeuvfbaQsoHAACAQ0KLQvv//M//RP/+/ePJJ5+M0047LSIinnnmmUbbHMzb5jdu3BgREccee2xERKxZsybq6upi1KhR+W3Kyspi+PDhsWTJkmZDe0NDQzQ0NOQfb9q0qeB64FCj/yl15gClTP9T6swBDlctenv8wIED45VXXomFCxfGwoULo2fPnjFv3rz844ULF8aDDz5YUCFZlsWNN94YZ599dgwePDgiIurq6iIiorKystG2lZWV+XX7mjZtWnTr1i2/9OnTp6B64FCk/yl15gClTP9T6swBDlctCu1ZljV6/Jvf/KbR588Pxg033BArV66Me+65p8m6fa/eZ1m23yv6N998c2zcuDG/rFu3rlXqg0OB/qfUmQOUMv1PqTMHOFwV9Jn2N+0b4gs1fvz4eOCBB+L3v/99HH/88fnxqqqqiNhzxb1Xr1758Q0bNjS5+v6msrKyKCsra5W64FCj/yl15gClTP9T6swBDlctutKey+WaXOE+mM+wZ1kWN9xwQ9x3333x4IMPxoABAxqtHzBgQFRVVUVtbW1+bMeOHbFo0aIYNmxYwccFAACAQ0GLrrRnWRZXXnll/i9Y27dvj+uuu67J3ePvu+++A9rf9ddfH3Pnzo3//M//jIqKivzn1Lt16xbl5eWRy+Vi4sSJMXXq1Bg4cGAMHDgwpk6dGp07d47LLrusJaUDAADAIadFoX3s2LGNHn/6058+qIPfcccdERExYsSIRuOzZ8+OK6+8MiIiJk2aFNu2bYtx48ZFfX19nHHGGbFgwYKoqKg4qGMDAABA6loU2mfPnt2qBz+Qz8Tncrmorq6O6urqVj02AAAApK5Fn2kHAAAA2o/QDgAAAIkS2gEAACBRQjsAAAAkSmgHAACARAntAAAAkCihHQAAABIltAMAAECihHYAAABIlNAOAAAAiRLaAQAAIFFCOwAAACRKaAcAAIBECe0AAACQKKEdAAAAEiW0AwAAQKKEdgAAAEiU0A4AAACJEtoBAAAgUUI7AAAAJEpoBwAAgEQJ7QAAAJAooR0AAAASJbQDAABAooR2AAAASJTQDgAAAIkS2gEAACBRQjsAAAAkSmgHAACARAntAAAAkCihHQAAABIltAMAAECihHYAAABIlNAOAAAAiRLaAQAAIFFCOwAAACRKaAcAAIBECe0AAACQKKEdAAAAEiW0AwAAQKKEdgAAAEiU0A4AAACJ6ljsAlL08vdnNhnrMf6GIlQCAABAKXOlHQAAABIltAMAAECihHYAAABIlNAOAAAAiRLaAQAAIFHuHg9ERMTttc/kf/7iyEFFrAQAAHiTK+0AAACQqKKG9t///vdx0UUXRe/evSOXy8UvfvGLRuuzLIvq6uro3bt3lJeXx4gRI2LVqlXFKRYAAADaWVFD+9atW+PUU0+NmTNnNrt++vTpUVNTEzNnzoylS5dGVVVVjBw5MjZv3tzOlQIAAED7K+pn2i+44IK44IILml2XZVnMmDEjJk+eHGPGjImIiDlz5kRlZWXMnTs3rr322vYsFQAAANpdsjeiW7NmTdTV1cWoUaPyY2VlZTF8+PBYsmTJfkN7Q0NDNDQ05B9v2rSpzWuFVOh/Sp05QCnT/5Q6c4DDVbI3oqurq4uIiMrKykbjlZWV+XXNmTZtWnTr1i2/9OnTp03rhJS0dv/fXvtMo7vKQ+r8bwClTP9T6swBDlfJhvY35XK5Ro+zLGsytrebb745Nm7cmF/WrVvX1iVCMvQ/pc4coJTpf0qdOcDhKtm3x1dVVUXEnivuvXr1yo9v2LChydX3vZWVlUVZWVmb1wcp0v+UOnOAUqb/KXXmAIerZK+0DxgwIKqqqqK2tjY/tmPHjli0aFEMGzasiJUBAABA+yjqlfYtW7bEX//61/zjNWvWxIoVK+LYY4+Nvn37xsSJE2Pq1KkxcODAGDhwYEydOjU6d+4cl112WRGrBgAAgPZR1ND+2GOPxXnnnZd/fOONN0ZExNixY+PHP/5xTJo0KbZt2xbjxo2L+vr6OOOMM2LBggVRUVFRrJIBAACg3RQ1tI8YMSKyLNvv+lwuF9XV1VFdXd1+RR2El78/s8lYj/E3FKESAAAADgfJfqYdAAAASp3QDgAAAIkS2gEAACBRQjsAAAAkSmgHAACARBX17vGHEneGp5TcXvtMsUuA9rVwWsR5Nxe7CgCAJlxpBwAAgEQJ7QAAAJAooR0AAAASJbQDAABAooR2AAAASJTQDgAAAInylW8HobmvgXu7bXxNHAAAAAfKlXYAAABIlNAOAAAAiRLaAQAAIFFCOwAAACRKaAcAAIBECe0AAACQKKEdAAAAEiW0AwAAQKKEdgAAAEiU0A4AAACJEtoBAAAgUR2LXQARL39/ZpOxHuNvKEIlAAAApMSVdgAAAEiU0A4AAACJEtoBAAAgUUI7AAAAJEpoBwAAgES5ezzQIrfXPhMREV8cOajZ8b3tuw0AANAyrrQDAABAooR2AAAASJTQDgAAAIkS2gEAACBRQjsAAAAkyt3j29nL359Z0HY9xt/QFuUko9R+XyBBC6dFnHdzsasAAGjElXYAAABIlNAOAAAAiRLaAQAAIFFCOwAAACRKaAcAAIBEuXv8IeJA7zq/r7a8C/uB1uRO8Ie+22ufabV1Xxw56IC3eXPdvo+b209z+3yrbeCALJy257/7u6u8O84DAG3MlXYAAABIlNAOAAAAiRLaAQAAIFFCOwAAACRKaAcAAIBECe0AAACQKF/5dphrza+KK3Rf+z7vQL4CrtBjFepAfl9fXdfYW33VWzGO1dw2B/LVbwfytXSUkDe/4u1Atmvuq97eHG9u/b5je2+7N18hBwDsxZV2AAAASNQhEdp/8IMfxIABA+Koo46K008/PR5++OFilwQAAABtLvnQfu+998bEiRNj8uTJsXz58jjnnHPiggsuiLVr1xa7NAAAAGhTyYf2mpqauOaaa+Kzn/1svOc974kZM2ZEnz594o477ih2aQAAANCmkr4R3Y4dO2LZsmXxla98pdH4qFGjYsmSJc0+p6GhIRoaGvKPN27cGBERmzZt2u9xNm/b1grVHl7KmjlfrXWe2nLfhTqQmprb5u282XdZlhVWWAsV0v/bt25p87r2tm8tzR3/zW3eqrYD2c/+tj3Q43PwDoU5EFu3v/2O933+1u17xt7879uN772uuW3f6lgcslLv/21btr3lejhYqc+Bra+/rv9pM63a/1nCXnjhhSwissWLFzca/9a3vpUNGjSo2edMmTIliwiLJanl2WefbY8po/8tyS7mgKWUF/1vKfXFHLCU8tIa/Z/Lsnb601cBXnzxxXjnO98ZS5YsibPOOis//q1vfSv+/d//PZ566qkmz9n3L2y7d++O1157LY477rjI5XJNtt+0aVP06dMn1q1bF127dm2bX+Qw4ny1zMaNG6Nv375RX18f73jHO9r8ePq/bTlfLZfyHPB6tpxz1jIp93+E17OlnK+WS3kOeD1bzjlrmdbs/6TfHt+9e/fo0KFD1NXVNRrfsGFDVFZWNvucsrKyKCsrazR2ICepa9eumq8FnK+WOeKI9rl9hP5vH85Xy6U8B7yeLeectUzK/R/h9Wwp56vlUp4DXs+Wc85apjX6P+kb0XXq1ClOP/30qK2tbTReW1sbw4YNK1JVAAAA0D6SvtIeEXHjjTfGZz7zmRg6dGicddZZ8cMf/jDWrl0b1113XbFLAwAAgDaVfGi/9NJL49VXX41vfOMbsX79+hg8eHD8+te/jn79+rXK/svKymLKlClN3kpD85yvlkn9fKVeX2qcr5ZL+ZylXFuqnLOWSf18pV5fapyvlkv5nKVcW6qcs5ZpzfOV9I3oAAAAoJQl/Zl2AAAAKGVCOwAAACRKaAcAAIBECe0AAACQqJIO7T/4wQ9iwIABcdRRR8Xpp58eDz/8cLFLSsbvf//7uOiii6J3796Ry+XiF7/4RaP1WZZFdXV19O7dO8rLy2PEiBGxatWq4hRbZNOmTYv3v//9UVFRET179ozRo0fH008/3WibVM+XOdA8/d8yh+oc0P/7Zw4cOP1/+NH/LWMOHF70f8u0V/+XbGi/9957Y+LEiTF58uRYvnx5nHPOOXHBBRfE2rVri11aErZu3RqnnnpqzJw5s9n106dPj5qampg5c2YsXbo0qqqqYuTIkbF58+Z2rrT4Fi1aFNdff308+uijUVtbGzt37oxRo0bF1q1b89ukeL7Mgf3T/y1zKM4B/f/WzIEDp/8PP/q/ZcyBw4v+b5l26/+sRH3gAx/IrrvuukZjJ510UvaVr3ylSBWlKyKy+++/P/949+7dWVVVVXbrrbfmx7Zv355169Ytu/POO4tQYVo2bNiQRUS2aNGiLMvSPV/mwIHR/y13KMwB/X/gzIGW0f+HF/3fcubA4UP/t1xb9X9JXmnfsWNHLFu2LEaNGtVofNSoUbFkyZIiVXXoWLNmTdTV1TU6f2VlZTF8+HDnLyI2btwYERHHHntsRKR5vsyBwqX4eqYm9Tmg/w9Oaq9navT/4S211zNF5sDhK7XXMkVt1f8lGdpfeeWV2LVrV1RWVjYar6ysjLq6uiJVdeh48xw5f01lWRY33nhjnH322TF48OCISPN8mQOFS/H1TMmhMAf0/8FJ7fVMif4//KX2eqbGHDi8pfZapqYt+79j65V56Mnlco0eZ1nWZIz9c/6auuGGG2LlypXxhz/8ocm6FM9XijUdKpy75h1KcyC1eg41zl9T+r90OH/NMwdKg3PXvLbs/5K80t69e/fo0KFDk79ubNiwoclfQWiqqqoqIsL528f48ePjgQceiIULF8bxxx+fH0/xfJkDhUvx9UzFoTIH9P/BSe31TIX+Lw2pvZ4pMQcOf6m9lilp6/4vydDeqVOnOP3006O2trbReG1tbQwbNqxIVR06BgwYEFVVVY3O344dO2LRokUlef6yLIsbbrgh7rvvvnjwwQdjwIABjdaneL7MgcKl+HoW26E2B/T/wUnt9Sw2/V9aUns9U2AOlI7UXssUtFv/F3pnvEPdvHnzsiOPPDKbNWtW9uc//zmbOHFi1qVLl+xvf/tbsUtLwubNm7Ply5dny5cvzyIiq6mpyZYvX54999xzWZZl2a233pp169Ytu++++7Innngi+9SnPpX16tUr27RpU5Erb3+f//zns27dumUPPfRQtn79+vzy+uuv57dJ8XyZA/un/1vmUJwD+v+tmQMHTv8ffvR/y5gDhxf93zLt1f8lG9qzLMv+5V/+JevXr1/WqVOn7LTTTsvfmp8sW7hwYRYRTZaxY8dmWbbn6wumTJmSVVVVZWVlZdm5556bPfHEE8UtukiaO08Rkc2ePTu/Tarnyxxonv5vmUN1Duj//TMHDpz+P/zo/5YxBw4v+r9l2qv/c/97MAAAACAxJfmZdgAAADgUCO0AAACQKKEdAAAAEiW0AwAAQKKEdgAAAEiU0A4AAACJEtoBAAAgUUI7AAAAJEpop12ce+65MXfu3HY/bkNDQ/Tt2zeWLVvW7seGN+l/Spn+p9SZA5Qy/d86hPZ2sGTJkujQoUN89KMfbdfj/u1vf4tcLhcrVqxo1+Pu65e//GXU1dXFJz/5yXY/dllZWdx0003x5S9/ud2PzR76X/+XMv2v/0udOWAOlDL9r/9bi9DeDn70ox/F+PHj4w9/+EOsXbu2XY65Y8eOZPb3ve99L6666qo44ojitNvll18eDz/8cPzlL38pyvFLnf7X/6VM/+v/UmcOmAOlTP/r/1aT0aa2bNmSVVRUZE899VR26aWXZl//+tcbrV+4cGEWEdkvf/nL7JRTTsnKysqyD3zgA9nKlSsbbbd48eLsnHPOyY466qjs+OOPz8aPH59t2bIlv75fv37ZN7/5zWzs2LFZ165dsyuuuCKLiEbL8OHDsyzLsuHDh2cTJkxotP+LL744Gzt27Fvu70Dq2NfLL7+c5XK57Mknn2w0HhHZnXfemX3sYx/LysvLs5NOOilbsmRJtnr16mz48OFZ586dszPPPDP761//mn/OlClTslNPPTWbNWtW1qdPn6xLly7Zddddl+3cuTP79re/nVVWVmY9evTIbrnlliZ1jBgxIvva17623zppG/pf/5cy/a//S505YA6UMv2v/1uT0N7GZs2alQ0dOjTLsiz7r//6r6x///7Z7t278+vfnLDvec97sgULFmQrV67M/uEf/iHr379/tmPHjizLsmzlypXZ0Ucfnd1+++3ZM888ky1evDgbMmRIduWVV+b3069fv6xr167Zbbfdlq1evTpbvXp19sc//jGLiOy3v/1ttn79+uzVV1/NsuzAJ+y++zuQOvZ1//33Z126dMl27drVaDwisne+853Zvffemz399NPZ6NGjs/79+2fnn39+Nn/+/OzPf/5zduaZZ2Yf/ehH88+ZMmVKdvTRR2ef+MQnslWrVmUPPPBA1qlTp+wjH/lINn78+Oypp57KfvSjH2URkT3yyCONjjdp0qRsxIgRB/CK0Zr0v/4vZfpf/5c6c8AcKGX6X/+3JqG9jQ0bNiybMWNGlmVZ9sYbb2Tdu3fPamtr8+vfnLDz5s3Lj7366qtZeXl5du+992ZZlmWf+cxnsn/6p39qtN+HH344O+KII7Jt27ZlWbZngo0ePbrRNmvWrMkiIlu+fHmj8QOdsPvu70Dq2Nftt9+evetd72oyHhHZP//zP+cfP/LII1lEZLNmzcqP3XPPPdlRRx2VfzxlypSsc+fO2aZNm/JjH/nIR7L+/fs3+gfhxBNPzKZNm9boeN/97nez/v37N1sjbUf/6/9Spv/1f6kzB8yBUqb/9X9r6ngwb63nrT399NPxxz/+Me67776IiOjYsWNceuml8aMf/Sg+/OEPN9r2rLPOyv987LHHxoknnpj//MWyZcvir3/9a/zkJz/Jb5NlWezevTvWrFkT73nPeyIiYujQoa1a/777O9A69rZt27Y46qijmt3/Kaeckv+5srIyIiLe9773NRrbvn17bNq0Kbp27RoREf3794+KiopG23To0KHRZ2UqKytjw4YNjY5VXl4er7/++tv+zrQe/a//S5n+1/+lzhwwB0qZ/tf/rU1ob0OzZs2KnTt3xjvf+c78WJZlceSRR0Z9fX0cc8wxb/n8XC4XERG7d++Oa6+9Nr7whS802aZv3775n7t06XJAdR1xxBGRZVmjsTfeeKPJdvvu70Dr2Fv37t2jvr6+2XVHHnlk/uc3f9fmxnbv3t3sc97cprmxvZ8TEfHaa69Fjx49mq2DtqH/9X8p0//6v9SZA+ZAKdP/+r+1Ce1tZOfOnXH33XfHd77znRg1alSjdR//+MfjJz/5Sdxwww35sUcffTTf9PX19fHMM8/ESSedFBERp512WqxatSpOOOGEFtXQqVOniIjYtWtXo/EePXrE+vXr84937doVTz75ZJx33nlvub9C6hgyZEjU1dUd0D9QbenJJ5+MIUOGFO34pUb/76H/S5P+30P/ly5zYA9zoDTp/z30f+vylW9t5Je//GXU19fHNddcE4MHD260fOITn4hZs2Y12v4b3/hG/O53v4snn3wyrrzyyujevXuMHj06IiK+/OUvxyOPPBLXX399rFixIlavXh0PPPBAjB8//i1r6NmzZ5SXl8f8+fPjpZdeio0bN0ZExPnnnx+/+tWv4le/+lU89dRTMW7cuPj73//+tr9TIXUMGTIkevToEYsXL37b/belhx9+uMk/nLQd/b+H/i9N+n8P/V+6zIE9zIHSpP/30P+tS2hvI7NmzYoPf/jD0a1btybrPv7xj8eKFSvi8ccfz4/deuutMWHChDj99NNj/fr18cADD+T/SnbKKafEokWLYvXq1XHOOefEkCFD4mtf+1r06tXrLWvo2LFjfO9734u77rorevfuHRdffHFERFx99dUxduzYuOKKK2L48OExYMCAt/0LW6F1dOjQIa6++upGn4Fpb4888khs3LgxPvGJTxSthlKj//fQ/6VJ/++h/0uXObCHOVCa9P8e+r915bJ9P9hAu3rooYfivPPOi/r6+njHO95R7HLaxEsvvRQnn3xyLFu2LPr169fux//Hf/zHGDJkSHz1q19t92Pz1vR/29P/6dL/bU//p80caHvmQLr0f9s7nPrflXbaXGVlZcyaNSvWrl3b7sduaGiIU089Nb74xS+2+7EhQv9T2vQ/pc4coJTp/9bjSnuRlcJf2WB/9D+lTP9T6swBSpn+pyWEdgAAAEiUt8cDAABAooR2AAAASJTQDgAAAIkS2gEAACBRQjsAAAAkSmgHAACARAntAAAAkCihHQAAABL1/wGHc8wV74vAyQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1200x400 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make a set of five histograms using plt.subplots and defining explicitly the axes\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=5, figsize=(12, 4), sharey=True)\n",
    "\n",
    "axes[0].hist(aperture_cluster1, bins=30, color = 'tab:red', alpha = 0.5)\n",
    "axes[1].hist(aperture_cluster2, bins=30, color = 'tab:blue', alpha = 0.5)\n",
    "axes[2].hist(aperture_background, bins=30, color = 'tab:orange', alpha = 0.5)\n",
    "axes[3].hist(aperture_bedding1, bins=5, color = 'tab:green', alpha = 0.5)\n",
    "axes[4].hist(aperture_bedding2, bins=5, color = 'tab:brown', alpha = 0.5)\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_xlabel('Aperture (mm)')\n",
    "    ax.set_xlim(-1, 20)\n",
    "    ax.set_ylim(0,50)\n",
    "\n",
    "axes[0].set_ylabel('Frequency')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# concatenate lists\n",
    "\n",
    "strike = np.concatenate((strike_cluster1, strike_cluster2, strike_background, strike_bedding1, strike_bedding2), axis=None)\n",
    "dip = np.concatenate((dip_cluster1, dip_cluster2, dip_background, dip_bedding1, dip_bedding2), axis=None)\n",
    "depth = np.concatenate((depth_cluster1, depth_cluster2, depth_background, depth_bedding1, depth_bedding2), axis=None)\n",
    "aperture = np.concatenate((aperture_cluster1, aperture_cluster2, aperture_background, aperture_bedding1, aperture_bedding2), axis=None)\n",
    "attribute = np.concatenate((attribute_cluster1, attribute_cluster2, attribute_background, attribute_bedding1, attribute_bedding2), axis=None)\n",
    "\n",
    "# convert to dataframe\n",
    "\n",
    "df = pd.DataFrame({'strike': strike, 'dip': dip, 'depth_m': depth, 'aperture_mm': aperture, 'type': attribute})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>strike</th>\n",
       "      <th>dip</th>\n",
       "      <th>depth_m</th>\n",
       "      <th>aperture_mm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>450.000000</td>\n",
       "      <td>450.000000</td>\n",
       "      <td>450.000000</td>\n",
       "      <td>450.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>188.856626</td>\n",
       "      <td>42.696792</td>\n",
       "      <td>752.122350</td>\n",
       "      <td>1.395988</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>75.794837</td>\n",
       "      <td>24.645009</td>\n",
       "      <td>157.483882</td>\n",
       "      <td>2.913713</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>5.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>157.653506</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>108.798179</td>\n",
       "      <td>17.358540</td>\n",
       "      <td>610.541891</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>186.981233</td>\n",
       "      <td>44.734488</td>\n",
       "      <td>752.270637</td>\n",
       "      <td>0.049566</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>265.714827</td>\n",
       "      <td>65.568442</td>\n",
       "      <td>895.314876</td>\n",
       "      <td>1.289191</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>348.000000</td>\n",
       "      <td>89.000000</td>\n",
       "      <td>1386.520969</td>\n",
       "      <td>18.159196</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           strike         dip      depth_m  aperture_mm\n",
       "count  450.000000  450.000000   450.000000   450.000000\n",
       "mean   188.856626   42.696792   752.122350     1.395988\n",
       "std     75.794837   24.645009   157.483882     2.913713\n",
       "min      5.000000    0.000000   157.653506     0.000000\n",
       "25%    108.798179   17.358540   610.541891     0.000000\n",
       "50%    186.981233   44.734488   752.270637     0.049566\n",
       "75%    265.714827   65.568442   895.314876     1.289191\n",
       "max    348.000000   89.000000  1386.520969    18.159196"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate dip azimuths\n",
    "\n",
    "dipazs = []\n",
    "\n",
    "for strike in df.strike.values:\n",
    "    dipazs.append(strike2dipaz(strike))\n",
    "\n",
    "df['dip_az'] = dipazs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# subselect a dataframe for export\n",
    "\n",
    "export_df = df[['depth_m','dip_az','dip', 'aperture_mm', 'type']].copy()\n",
    "export_df = export_df[(export_df.depth_m > 500) & (export_df.depth_m < 1000)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>depth_m</th>\n",
       "      <th>dip_az</th>\n",
       "      <th>dip</th>\n",
       "      <th>aperture_mm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>446.000000</td>\n",
       "      <td>446.000000</td>\n",
       "      <td>446.000000</td>\n",
       "      <td>446.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>752.315736</td>\n",
       "      <td>201.734263</td>\n",
       "      <td>42.568513</td>\n",
       "      <td>1.395895</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>150.197112</td>\n",
       "      <td>111.337939</td>\n",
       "      <td>24.712021</td>\n",
       "      <td>2.920963</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>500.034701</td>\n",
       "      <td>0.013227</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>611.794292</td>\n",
       "      <td>168.171203</td>\n",
       "      <td>17.283704</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>752.270637</td>\n",
       "      <td>239.087762</td>\n",
       "      <td>44.379844</td>\n",
       "      <td>0.049566</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>895.062074</td>\n",
       "      <td>284.395986</td>\n",
       "      <td>65.568442</td>\n",
       "      <td>1.289191</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>998.242393</td>\n",
       "      <td>359.885757</td>\n",
       "      <td>89.000000</td>\n",
       "      <td>18.159196</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          depth_m      dip_az         dip  aperture_mm\n",
       "count  446.000000  446.000000  446.000000   446.000000\n",
       "mean   752.315736  201.734263   42.568513     1.395895\n",
       "std    150.197112  111.337939   24.712021     2.920963\n",
       "min    500.034701    0.013227    0.000000     0.000000\n",
       "25%    611.794292  168.171203   17.283704     0.000000\n",
       "50%    752.270637  239.087762   44.379844     0.049566\n",
       "75%    895.062074  284.395986   65.568442     1.289191\n",
       "max    998.242393  359.885757   89.000000    18.159196"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "export_df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "export_df.to_csv('0_Synthetic_data.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bhi_env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}