Example_3_Evaluation.ipynb 105 KB
Newer Older
Maximilian Schanner's avatar
Maximilian Schanner committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluation\n",
    "The third and last step in the `CORBASS` algorithm consists of evaluating the output from the integration step. Thus again this example makes use of the output from the previous example. A standard `run` of `CORBASS` concludes by calculating mean, variance and percentiles of the Gauss coefficient (compound) posterior. Therefore we first load the results from  the `Integration` step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sys\n",
    "import os\n",
    "import numpy as np\n",
    "\n",
    "# relative import\n",
    "sys.path.append(os.path.abspath('') + '/../')\n",
    "from corbass.utils import load\n",
    "from corbass.integration import IntegrationResult\n",
    "\n",
    "pars = load('./Example_Parfile.py')\n",
    "# This example focuses again on the interval 1400-1500 A.D.\n",
    "year = 1450\n",
    "\n",
33
    "with np.load(f'{pars.bin_fname(year)}{IntegrationResult.suffix_large}') as fh:\n",
Maximilian Schanner's avatar
Maximilian Schanner committed
34
35
36
37
38
39
40
41
42
43
44
45
46
47
    "    posterior = fh['posterior']\n",
    "    mu_coeffs = fh['mu_coeffs']\n",
    "    cov_coeffs = fh['cov_coeffs']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then pass the results to the routine `evaluation.coeffs`:"
   ]
  },
  {
   "cell_type": "code",
48
   "execution_count": 2,
Maximilian Schanner's avatar
Maximilian Schanner committed
49
   "metadata": {},
50
51
52
53
54
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
55
      "/home/arthus/Dokumente/CORBASS/examples/../corbass/evaluation.py:217: RuntimeWarning: covariance is not positive-semidefinite.\n",
56
57
58
59
      "  for it in par_samps]\n"
     ]
    }
   ],
Maximilian Schanner's avatar
Maximilian Schanner committed
60
61
62
63
64
65
66
67
68
69
70
71
72
73
   "source": [
    "from corbass.evaluation import coeffs\n",
    "ls, ms, mean, sd, err_16, err_84 = coeffs(posterior, mu_coeffs, cov_coeffs, r_ref=pars.r_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For illustration we show the dipole coefficients together with errorbars derived from the percentiles and one standard deviation:"
   ]
  },
  {
   "cell_type": "code",
74
   "execution_count": 3,
Maximilian Schanner's avatar
Maximilian Schanner committed
75
76
   "metadata": {},
   "outputs": [
77
78
    {
     "data": {
Maximilian Schanner's avatar
Maximilian Schanner committed
79
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD4CAYAAAD//dEpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAbjklEQVR4nO3df5BU5Z3v8ffnggLXCPIriWFAJoAmSkaIHTTc1bCRDOzdTTBZubJmS5KQIpiQzdb9x1jG5a6WVTHJlluGjYZ7dVFXrhqtLJTGa2D9QbYkxEHIgFHWISTSQiUjMxCMwjLke//oZ2LP0DNzmJ7pZmY+r6pTnvM95+l+upzhM+c8p8+jiMDMzKwn/6XaHTAzs4HBgWFmZpk4MMzMLBMHhpmZZeLAMDOzTIZXuwP9ZcKECTF16tRqd8PMbEDZtm3bGxExsdS+fg8MSb8CjgAngLaIyEkaBzwMTAV+BfyPiGhNx98ILEvH/01EPJXqlwBrgVHAj4CvRTf3BE+dOpWGhob++VBmZoOUpF93ta9Sl6T+NCJmRUQubX8d+LeImAH8W9pG0oXAEuAiYCHwPUnDUpu7gOXAjLQsrFDfzcyM6o1hLALuS+v3AVcV1R+KiGMRsRdoAuZIOhcYHRFb0lnF/UVtzMysAioRGAH8WNI2SctT7T0RcQAg/ffdqT4J2FfUNp9qk9J653oHkpZLapDU0Nzc3Mcfw8xsaKvEoPd/i4j9kt4NbJT0SjfHqkQtuql3LESsAdYA5HI5P/PEzKwP9fsZRkTsT//9LfBDYA7wm3SZifTf36bD88DkouY1wP5UrylRNzOzCunXwJB0lqSz29eBemAXsAFYmg5bCqxP6xuAJZJGSKqlMLj9s3TZ6oikyyQJuK6ojZmZVUB/X5J6D/DDwr/xDAfWRcT/k/QC8IikZcBrwGKAiHhJ0iPAL4A24CsRcSK91vW8c1vtk2kxM7MK0WB9vHkulwt/D8PM7NRI2lb0FYgO/GgQM7PT0DXf38I1399S7W504MAwMxtM/vnPC0s/cGCUcDomu5lZtTkwzMxOI3c/t4fn97zRofb8nje4+7k93Tf893+EvZs71vZuLtT7iAOjL/XjqaCZDQ11NWNYuW47h98+DhTCYuW67dTVjOm+4aQPww8+B0cPFbb3bi5sT/pwn/XNgVHkdE52Mxsa5k6bwOprZ9P02zfJt77FynXbWX3tbOZOm9B9w9orYPFaaH4FDv26EBaL1xbqfcSBUeR0TnYzGzrmTpvAe0aP4PVDR/nrS6f0HBbtaq+As8+Fw/sgt6xPwwIcGB2czsluZkPH83ve4De/O8akc0byL1tfO+nKR5f2boYjB2DMZGi45+QrH2VyYHRyuia7mQ0N7Vc2pr/7XdSM/a+svnY2K9dt7zk02q9sTPwAnHNe4Y/WH3yuT0PDgdHJ6ZrsZjY0NOYPs/ra2YwZdQbwzpWPxvzh7hu+/mIhJEaeU9huv/Lx+ot91jcHRpHTOdnNbGhY8bFpJ13ZmDttAis+Nq37hn/ytydf2ai9olDvIw6MIqdzspvZ0PLwlz7Kw1/66Kk3/PwThaUf+OGDJbR/y7tX/7PMzAYwP3zQzMzKVokpWgccn1mYmZ3MZxhmZpbJgAoMSQsl7ZbUJOnr1e6PmdlQMmACQ9Iw4J+APwMuBP5K0oXV7ZWZ2dAxYAIDmAM0RcQvI+I/gYeARVXuk5nZkDGQAmMSsK9oO59qfyRpuaQGSQ3Nzc0V7ZyZ2WA3kAJDJWodvkQSEWsiIhcRuYkTJ1aoW2ZmQ8NACow8MLlouwbYX6W+mJkNOQMpMF4AZkiqlXQmsATYUOU+mZkNGQPmi3sR0SZpJfAUMAy4NyJeqnK3zMyGjAETGAAR8SPgR9Xuh5nZUDSQLkmZmVkVOTDMzCwTB4aZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLxIFhZmaZODDMzCwTB4aZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJv0WGJL+l6TXJe1Iy38v2nejpCZJuyUtKKpfImln2nenJKX6CEkPp/pWSVP7q99mZlZaf59h3BERs9LyIwBJF1KYXvUiYCHwPUnD0vF3AcuBGWlZmOrLgNaImA7cAdzez/02M7NOqnFJahHwUEQci4i9QBMwR9K5wOiI2BIRAdwPXFXU5r60/ihwZfvZh5mZVUZ/B8ZKSY2S7pU0NtUmAfuKjsmn2qS03rneoU1EtAGHgfGd30zSckkNkhqam5v79pOYmQ1xZQWGpE2SdpVYFlG4vDQNmAUcAP6hvVmJl4pu6t216ViIWBMRuYjITZw48ZQ/j5mZdW14OY0jYn6W4yT9b+DxtJkHJhftrgH2p3pNiXpxm7yk4cAYoKX3PTczs1PVn3dJnVu0+WlgV1rfACxJdz7VUhjc/llEHACOSLosjU9cB6wvarM0rV8NPJ3GOczMrELKOsPowbckzaJw6ehXwJcAIuIlSY8AvwDagK9ExInU5npgLTAKeDItAPcAD0hqonBmsaQf+21mZiVosP6hnsvloqGhodrdMDMbUCRti4hcqX3+preZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLxIFhZmaZODDMzCwTB4aZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLxIFhZmaZlDun92JJL0n6g6Rcp303SmqStFvSgqL6JZJ2pn13ptn1SDPwPZzqWyVNLWqzVNKraVmKmZlVXLlnGLuAzwCbi4uSLqQwK95FwELge5KGpd13AcspTM06I+0HWAa0RsR04A7g9vRa44BVwKXAHGCVpLFl9tvMzE5RWYERES9HxO4SuxYBD0XEsYjYCzQBc9I836MjYkuak/t+4KqiNvel9UeBK9PZxwJgY0S0REQrsJF3QsbMzCqkv8YwJgH7irbzqTYprXeud2gTEW3AYWB8N691EknLJTVIamhubu6Dj2FmZu2G93SApE3Ae0vsuiki1nfVrEQtuqn3tk3HYsQaYA0U5vTuom9mZtYLPQZGRMzvxevmgclF2zXA/lSvKVEvbpOXNBwYA7Sk+rxObZ7tRZ/MzKwM/XVJagOwJN35VEthcPtnEXEAOCLpsjQ+cR2wvqhN+x1QVwNPp3GOp4B6SWPTYHd9qpmZWQX1eIbRHUmfBr4LTASekLQjIhZExEuSHgF+AbQBX4mIE6nZ9cBaYBTwZFoA7gEekNRE4cxiCUBEtEi6FXghHXdLRLSU028zMzt1KvwRP/jkcrloaGiodjfMzAYUSdsiIldqn7/pbWZmmTgwzMwsEweGmZll4sAwM7NMHBhmZpaJA8PMzDJxYJiZWSYODDMzy8SBYWZmmTgwzMwsEweGmZll4sAwM7NMHBhmZpaJA8PMzDJxYJiZWSYODDMzy6SswJC0WNJLkv4gKVdUnyrpbUk70nJ30b5LJO2U1CTpzjRVK2k614dTfaukqUVtlkp6NS1LMTOziiv3DGMX8Blgc4l9eyJiVlpWFNXvApZTmOd7BrAw1ZcBrRExHbgDuB1A0jhgFXApMAdYleb2NjOzCiorMCLi5YjYnfV4SecCoyNiSxTmhr0fuCrtXgTcl9YfBa5MZx8LgI0R0RIRrcBG3gkZMzOrkP4cw6iVtF3Sc5IuT7VJQL7omHyqte/bBxARbcBhYHxxvUSbDiQtl9QgqaG5ubnvPomZmTG8pwMkbQLeW2LXTRGxvotmB4ApEXFQ0iXAv0q6CFCJY6P9rbrY112bjsWINcAagFwuV/IYMzPrnR4DIyLmn+qLRsQx4Fha3yZpD3A+hbODmqJDa4D9aT0PTAbykoYDY4CWVJ/Xqc2zp9onMzMrT79ckpI0UdKwtP5+CoPbv4yIA8ARSZel8YnrgPazlA1A+x1QVwNPp3GOp4B6SWPTYHd9qpmZWQX1eIbRHUmfBr4LTASekLQjIhYAVwC3SGoDTgArIqIlNbseWAuMAp5MC8A9wAOSmiicWSwBiIgWSbcCL6Tjbil6LTMzqxAV/ogffHK5XDQ0NFS7G2ZmA4qkbRGRK7XP3/Q2M7NMHBhmZpaJA8PMzDJxYJiZWSYODDMzy8SBYWZmmTgwzMwsEweGmZll4sAwM7NMHBhmZpaJA8PMzDJxYJiZWSYODDMzy8SBYWZmmTgwzMwsk7ICQ9K3Jb0iqVHSDyWdU7TvRklNknZLWlBUv0TSzrTvzjTzHpJGSHo41bdKmlrUZqmkV9OyFDMzq7hyzzA2AjMjog74D+BGAEkXUpgx7yJgIfC99ilbgbuA5RSmbZ2R9gMsA1ojYjpwB3B7eq1xwCrgUmAOsCpN1WpmZhVUVmBExI8joi1t/hSoSeuLgIci4lhE7AWagDmSzgVGR8SWNF/3/cBVRW3uS+uPAlems48FwMaIaImIVgoh1R4yZmZWIX05hvEF3pmfexKwr2hfPtUmpfXO9Q5tUggdBsZ381onkbRcUoOkhubm5rI+jJmZdTS8pwMkbQLeW2LXTRGxPh1zE9AGPNjerMTx0U29t206FiPWAGugMKd3qWPMzKx3egyMiJjf3f40CP0XwJXpMhMUzgImFx1WA+xP9ZoS9eI2eUnDgTFAS6rP69Tm2Z76bWZmfavcu6QWAjcAn4qIt4p2bQCWpDufaikMbv8sIg4ARyRdlsYnrgPWF7VpvwPqauDpFEBPAfWSxqbB7vpUMzOzCurxDKMHq4ERwMZ0d+xPI2JFRLwk6RHgFxQuVX0lIk6kNtcDa4FRFMY82sc97gEekNRE4cxiCUBEtEi6FXghHXdLRLSU2W8zMztFeucq0uCSy+WioaGh2t0wMxtQJG2LiFypff6mt5mZZeLAMDOzTBwYZmaWiQPDzMwycWCYmVkmDgwzM8vEgWFmZpk4MMzMLBMHhpmZZeLAMDOzTBwYZmaWiQPDzMwycWCYmVkmDgwzM8vEgWFmZpmUO+PetyW9IqlR0g8lnZPqUyW9LWlHWu4uanOJpJ2SmiTdmWbeI83O93Cqb5U0tajNUkmvpmVp536YmVn/K/cMYyMwMyLqgP8AbizatyciZqVlRVH9LmA5hWlbZwALU30Z0BoR04E7gNsBJI0DVgGXAnOAVWmqVjMzq6CyAiMifhwRbWnzp0BNd8dLOhcYHRFb0nzd9wNXpd2LgPvS+qPAlensYwGwMSJaIqKVQkgtxMzMKqovxzC+wDvzcwPUStou6TlJl6faJCBfdEw+1dr37QNIIXQYGF9cL9HGzMwqZHhPB0jaBLy3xK6bImJ9OuYmoA14MO07AEyJiIOSLgH+VdJFgEq8Tvuk4l3t665N574up3C5iylTppT+QGZm1is9BkZEzO9ufxqE/gvgynSZiYg4BhxL69sk7QHOp3B2UHzZqgbYn9bzwGQgL2k4MAZoSfV5ndo820Vf1wBrAHK5XMlQMTOz3in3LqmFwA3ApyLiraL6REnD0vr7KQxu/zIiDgBHJF2WxieuA9anZhuA9jugrgaeTgH0FFAvaWwa7K5PNTMzq6AezzB6sBoYAWxMd8f+NN0RdQVwi6Q24ASwIiJaUpvrgbXAKApjHu3jHvcAD0hqonBmsQQgIlok3Qq8kI67pei1zMysQpSuIg06uVwuGhoaqt0NM7MBRdK2iMiV2udvepuZWSYODDMzy8SBYWZmmTgwzMwsEweGmZll4sAwM7NMHBhmZpaJA8PMzDJxYJiZWSYODDMzy8SBYWZmmTgwzMwsEweGmZll4sAwM7NMHBhmZpaJA8PMzDIpd4rWWyU1Stoh6ceS3le070ZJTZJ2S1pQVL9E0s607840VSuSRkh6ONW3Sppa1GappFfTshQzM6u4cs8wvh0RdRExC3gc+DsASRdSmGL1ImAh8L32Ob6Bu4DlFOb5npH2AywDWiNiOnAHcHt6rXHAKuBSYA6wKs3tbWZmFVRWYETE74o2zwLa53tdBDwUEcciYi/QBMyRdC4wOiK2RGFu2PuBq4ra3JfWHwWuTGcfC4CNEdESEa3ARt4JGTMzq5Dh5b6ApNuA64DDwJ+m8iTgp0WH5VPteFrvXG9vsw8gItokHQbGF9dLtOncl+UUzl6YMmVKrz+TmZmdrMczDEmbJO0qsSwCiIibImIy8CCwsr1ZiZeKbuq9bdOxGLEmInIRkZs4cWJ3H8vMzE5Rj2cYETE/42utA56gMN6QByYX7asB9qd6TYk6RW3ykoYDY4CWVJ/Xqc2zGftkZmZ9pNy7pGYUbX4KeCWtbwCWpDufaikMbv8sIg4ARyRdlsYnrgPWF7VpvwPqauDpNM7xFFAvaWwa7K5PNTMzq6ByxzC+KekC4A/Ar4EVABHxkqRHgF8AbcBXIuJEanM9sBYYBTyZFoB7gAckNVE4s1iSXqtF0q3AC+m4WyKipcx+m5nZKVLhj/jBJ5fLRUNDQ4fa8ePHyefzHD16tEq9qqyRI0dSU1PDGWecUe2umNkAIWlbRORK7Sv7LqmBJJ/Pc/bZZzN16lTS9wUHrYjg4MGD5PN5amtrq90dMxsEhtSjQY4ePcr48eMHfVgASGL8+PFD5mzKzPrfkAoM4JTD4prvb+Ga72/pp970r6EQjGZWOUMuMMzMrHccGF24+7k9PL/njQ615/e8wd3P7SnrdYcNG8asWbOYOXMmixcv5q233irr9czMKsWB0YW6mjGsXLedw28fBwphsXLddupqxpT1uqNGjWLHjh3s2rWLM888k7vvvjtTu7a2trLe18ysXA6MLsydNoHV186m6bdvkm99i5XrtrP62tnMnTahz97j8ssvp6mpid///vd84Qtf4CMf+QizZ89m/frCdxnXrl3L4sWL+eQnP0l9fT1vvvkmn//85/nQhz5EXV0djz32WJ/1xcysJ0PqttpTNXfaBN4zegSvHzrK33x8ep+GRVtbG08++SQLFy7ktttu4+Mf/zj33nsvhw4dYs6cOcyfX3giy5YtW2hsbGTcuHHccMMNjBkzhp07dwLQ2traZ/0xM+uJA6Mbz+95g9/87hiTzhnJv2x9jcumjS87NN5++21mzZoFFM4wli1bxty5c9mwYQPf+c53gMLtv6+99hoAn/jEJxg3bhwAmzZt4qGHHvrja40d62lBzKxyHBhdaB+zmP7udzFm1Bl8bf6MPrks1T6GUSwieOyxx7jgggs61Ldu3cpZZ53V4TjfKmtm1eIxjC405g+z+trZjBlVeKxG+5hGY/5wn7/XggUL+O53v0v7Y1q2b99e8rj6+npWr179x21fkjKzSnJgdGHFx6addCYxd9oEVnxsWp+/180338zx48epq6tj5syZ3HzzzSWP+8Y3vkFrayszZ87k4osv5plnngHgi1/8Ip2fm2Vm1teG1MMHX375ZT74wQ9WqUfVMRQ/s5n1XncPH/QZhpmZZeLAMDOzTMqdce9WSY2Sdkj6saT3pfpUSW+n+g5Jdxe1uUTSTklNku5MM++RZud7ONW3Sppa1GappFfTsrRzP8zMrP+Ve4bx7Yioi4hZwOPA3xXt2xMRs9Kyoqh+F7CcwrStM4CFqb4MaI2I6cAdwO0AksZRmCf8UmAOsCpN1Wp22ur1U47/+c8Li9lpqKzAiIjfFW2eBXQ7gi7pXGB0RGxJ83XfD1yVdi8C7kvrjwJXprOPBcDGiGiJiFZgI++ETP/zL7CZGdAHYxiSbpO0D/gsHc8waiVtl/ScpMtTbRKQLzomn2rt+/YBREQbcBgYX1wv0aZzX5ZLapDU0NzcXOYnMzt1vX7K8b//I+zd3LG2d3Ohbnaa6DEwJG2StKvEsgggIm6KiMnAg8DK1OwAMCUiZgP/E1gnaTRQ6mvK7WclXe3rrk3HYsSaiMhFRG7ixIk9fbTu9eMv8G233cZFF11EXV0ds2bNYuvWrcybN48LLriAuro6PvCBD7By5UoOHTpU9ntZZfX6KceTPgw/+BwcTf/P924ubE/6cL/21+xU9BgYETE/ImaWWNZ3OnQd8JepzbGIOJjWtwF7gPMpnB3UFLWpAfan9TwwGUDScGAM0FJcL9Gm//TTL/CWLVt4/PHHefHFF2lsbGTTpk1Mnlz4eA8++CCNjY00NjYyYsQIFi1aVN5nsIrr9VOOa6+AxWuh+RU49OvCz9ritYW62Wmi3LukZhRtfgp4JdUnShqW1t9PYXD7lxFxADgi6bI0PnEd0B48G4D2O6CuBp5O4xxPAfWSxqbB7vpU61/99At84MABJkyYwIgRIwCYMGEC73vf+zocc+aZZ/Ktb32L1157jZ///OdlvZ9VXvFTjv/60inZnz1WewWcfS4c3ge5ZQ4LO+2UO4bxzXR5qpHCP+RfS/UrgEZJP6cwgL0iIlrSvuuB/wM0UTjzeDLV7wHGS2qicBnr6wCp3a3AC2m5pei1+lc//ALX19ezb98+zj//fL785S/z3HPPlTxu2LBhXHzxxbzyyitlv6dVVuenHHce0+jS3s1w5ACMmQwN95x8SdSsysp6Wm1E/GUX9ceAkrP7REQDMLNE/SiwuIs29wL39r6nvdT5F7j28rJD413vehfbtm3jJz/5Cc888wzXXHMN3/zmN0seO1gf2zKY9fopx+2XPCd+AEaeA/Nu8GUpO+34m95dKf4FPue8wi/uDz7XJ3/1DRs2jHnz5vH3f//3rF69uuTMeSdOnGDnzp1+DtQA0+unHL/+YuFnbOQ5he32S6Kvv9iv/TU7FQ6MrvTTL/Du3bt59dVX/7i9Y8cOzjvvvA7HHD9+nBtvvJHJkydTV1dX1vtZZfX6Kcd/8rcnn0nUXlGom50mPIFSV9p/UZ+9/Z1a7RVlXx548803+epXv8qhQ4cYPnw406dPZ82aNVx99dV89rOfZcSIERw7doz58+f/cW5vG0I+/0S1e2DWJT/efJAbip/ZzHrPjzc3M7OyOTDMzCyTIRcYg/USXClD6bOaWf8bUoExcuRIDh48OCT+IY0IDh48yMiRI6vdFTMbJIbUXVI1NTXk83mGypNsR44cSU1NTc8HmpllMKQC44wzzqC2trba3TAzG5CG1CUpMzPrPQeGmZll4sAwM7NMBu03vSU1A7/up5efAGR8ZrVZSf4ZsnL118/QeRFRcsrSQRsY/UlSQ1dfnTfLwj9DVq5q/Az5kpSZmWXiwDAzs0wcGL2zptodsAHPP0NWror/DHkMw8zMMvEZhpmZZeLAMDOzTBwYvSRpsaSXJP1Bkm+PtEwkLZS0W1KTpK9Xuz828Ei6V9JvJe2q9Hs7MHpvF/AZYHO1O2IDg6RhwD8BfwZcCPyVpAur2ysbgNYCC6vxxg6MXoqIlyNid7X7YQPKHKApIn4ZEf8JPAQsqnKfbICJiM1ASzXe24FhVjmTgH1F2/lUMxsQhtR8GKdK0ibgvSV23RQR6yvdHxvwVKLm+9ptwHBgdCMi5le7Dzao5IHJRds1wP4q9cXslPmSlFnlvADMkFQr6UxgCbChyn0yy8yB0UuSPi0pD3wUeELSU9Xuk53eIqINWAk8BbwMPBIRL1W3VzbQSPq/wBbgAkl5Scsq9t5+NIiZmWXhMwwzM8vEgWFmZpk4MMzMLBMHhpmZZeLAMDOzTBwYZmaWiQPDzMwy+f/rtEoKpw2quwAAAABJRU5ErkJggg==\n",
Maximilian Schanner's avatar
Maximilian Schanner committed
80
81
82
83
84
85
86
87
88
89
90
91
92
93
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(ms[0:3]-0.05, mean[0:3], yerr=[mean[0:3]-err_16[0:3], err_84[0:3]-mean[0:3]],\n",
94
    "            marker='x', ls='', label='Perc.');\n",
Maximilian Schanner's avatar
Maximilian Schanner committed
95
96
    "\n",
    "ax.errorbar(ms[0:3]+0.05, mean[0:3], yerr=sd[0:3],\n",
97
98
99
    "            marker='x', ls='', label='SD');\n",
    "ax.set_xticks([-1, 0, 1]);\n",
    "ax.legend();"
Maximilian Schanner's avatar
Maximilian Schanner committed
100
101
102
103
104
105
106
107
108
109
110
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though this is the standard output from `CORBASS`, having the full information of the (compound) Gauss coeffient posterior available on the grid allows for calculation of other quantities of interest. For example we can calculate the full power spectrum:"
   ]
  },
  {
   "cell_type": "code",
111
   "execution_count": 4,
Maximilian Schanner's avatar
Maximilian Schanner committed
112
113
   "metadata": {},
   "outputs": [
114
115
    {
     "data": {
Maximilian Schanner's avatar
Maximilian Schanner committed
116
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD4CAYAAAANbUbJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXyV5Z338c+VfSUJ2SGsCSSEzSAKogKCpqJirbaPigt1t5badjrT6tPptH1m+kxnujx1xKlacamC1XFp3bDBDbUgW0DCDmENZE9IgJD1XM8fJ6HsW07OfXLu7/v18kVyguf+vg7ky53r3NfvNtZaREQkuIU4HUBERHqeyl5ExAVU9iIiLqCyFxFxAZW9iIgLhDkd4GRSUlLs4MGDnY4hItKrrFq1qsZam3qyrwVU2RtjZgIzc3JyWLlypdNxRER6FWPMrlN9LaCWcay1b1tr709ISHA6iohIUAmoshcRkZ6hshcRcQGVvYiIC6jsRURcQGUvIuICKnsRERcIqLI3xsw0xjzd0NDgdBQRkaASUGWv6+xFxM1ufmopNz+1tEeeO6DKXkREeobKXkTEBVT2IiIuoLIXEXEBlb2IiAuo7EVEXEBlLyLiAip7EREXCKiy1w5aEXFCT25mChQBVfbaQSsi0jMCquxFRKRnqOxFRFxAZS8i4gIqexERF1DZi4i4gMpeRMQFVPYiIi6gshcRcdiTi0tZUlpzzGNLSmt4cnGpz46hshcRcdiYrATmLFhNw+E2wFv0cxasZkyW7zaYquxFRBx28eC+PDw9hy2VB9hSeYA5C1Yzd1YBk7JTfHaMMJ89kw8YY2YCM3NycpyOIiJ+0DWP5pUHLnE4if9VNjazeEs1i7dU8/nWmiNn9fVNbTw8LcenRQ8BVvbW2reBt8ePH3+f01lERHyptd3Dql31LN5SzSebq9hUcQCA9D6RfGVkOpkJ0Tzx8TbS+0Ty0rLdTMxODt4zexGRYFJW3+Q9e99czZLSWg62tBMeahg/qC+PzMhjyvBU8jLiWbq9ljkLVpOTFkdCdDjfvXKYz5dyVPYiIj7S3NbB8h11R87eS6sPAdA/MZqvXtCPKcNTmZSTQlzksdW7tqyBubMKeOyDrQBMyk5h7qwC1pY1qOxFRALBjppDLN5cxeIt1SzdXktzm4eIsBAmDOnLrAmDmDI8lezUWIwxp3yOB6dkAxwpe/AWvpZxREQc0tTaztLS2iNvru6qbQJgSEost1w0kCm5qUwckkx0RKjDSY+lshcROQ1rLVurDrJ4s7fcl++oo7XDQ3R4KJOyk7n3siFMHp7KoORYp6OelspeRFzrycWlJ2xcWlJaw4oddeRmxB95c3VfQzMAw9PjmD1pEFNz0xg/OInIsMA6ez8dlb2IuFbXztW0uAjCQkP40Wtreb24DI+1eCzER4ZxaU4K35meypThqfRLjHY68nlT2YuIK3k8lsiwECYM7svC9RUArNvXyODkGK4ZncnU3DQKBiYSHhocgwZU9iLiGh6PZfWeet5ZW87CkgoqGpuJCA0hMiyElnYP91w2mJ9cN9LpmD1CZS8iQc3jsRTvrufdkmMLfvLwVB6ZkUd8VBgPvLiK/olRvLl6H9NHpPt8VEEgUNmLuFCwz6TpKvh31pbz/rrOgg8LYcrwVB4Zncf0EWnER4UfmS7ZkztXA4XKXkSCgsdjWbW7nndPUvCPjsljWp634I/mj52rgUJlLyK91tEFv3BdOZWNLUSEhTD1NAV/NH/sXA0UAVX2GnEsImfi8VhW7qrnvZITC/7aMZlnLHi3Cqiy14hjETmZroJ/d+0+Fq6roOqAt+CvyE3lmtGZTB+RfsJwMTmWXh0RCUgdHsvKnXWdZ/Dego8MC2GqCv686JUSEb871ZiCNXv2c+HAJO9lkusqqD6q4K8d049peWkq+POkV01E/O7ImIL4SPpEhTHv8+385/ubiQoPpeFwG5FhIVyRm8Y1nWvwKvju0ysoIn53ydBk/rFwOD/+8zoMsHxnPeEhhkuGJnPNmEym56URq4L3Kb2aIuI3lY3N/Hn1Xt4o3svmSu89WC1w9cgMfvO/xqrge5BeWRHpUYdbOyjaUMHrxXv5fGs1HgsFAxO569LBvLh0F+l9Ilm+s44vy/YH5fXtgUJlL+JnwT6qALyXSi7fWccbxWW8V1LBwZZ2+idG89DUHG4c15+KxmbXjCkIFCp7EfGZHTWHeLO4jDdW76Ws/jCxEaFcMzqTG8dlMWFIX0JCvPdhLdpQ6ZoxBYFCZS8i3dLQ1MY7Jft4fVUZxbv3YwxclpPCPxbmUjgynZiIE2vGTWMKAoXKXkTOWVuHh0+3VPNG8V4Wbayktd3DsLQ4HpmRxw0X9CcjIcrpiHIclb2InBVrLev3NfJ6cRlvrdlH7aFW+sZGMOvigdw0LotR/ftgjHE6ppyCyl5ETuv4yyUjQkOYPiKNG8dlMTU3NWhu2xfsVPYicoJTXS75rzeMYuaYTBJjIpyOKOdIZS/iIqeaSbO2rIH7Lx96ysslvzauP9mpcQ6lFl9Q2Yu4yNEzaRKivbfl+9ZLxUzLS2Xyrz4+crnkjNGZ3HTc5ZLSu6nsRVyk63r2O+ctJzYylDvnLafdY/nzmn1clpPCDwqH85WRGSe9XFJ6Xk9utNOfqIhL1B5s4a/rK1m4rpx2j6XhcDtJMeHcPzmbGwr6kZkQ7XRE6UEqe5EgVnOwhffXVfBeSTlfbK/FYyG9TyQhBpJjI+iwMHZAgoreBVT2IkGm6kAzf11Xwbsl5SzfUYfHwtCUWB6amkO/xCh+XbSF4enxmknjMgFV9rrhuMj5qWxs5v3Ogl+xsw5rITs1ljlX5HDNmExy0+MxxvDk4lLNpHGpgCp73XBc5OxVNDSzcF0575WUs3JXPdbCsLQ4Hp42jGvHZDIsLe6EHa2aSeNeAVX2InJ6+/YfZmHnGvyqXfUA5KbH873pw7lmdAbD0uMdTiiBSmUvEuDK6puOLNGs3r0fgBGZffjBVcOZMTqTnDRtdpIzU9mLBKA9dU28V1LOe+sq+HKPt+BH9uvDP30llxmjMhiq3axyjlT2IgFid20T75aUs3BdOWvLGgAY3T+BH12dx4xRGQxOiXU4ofRmQVX2brjdm/ReJ5tL82ZxGW+u2UvdoVbW7W0EYGxWAo/OyGPGqEwGJsc4EVWCUFCVvcjpOH0y0DWXJik6nJYOD5P/82N21zUBcMGARH58zQiuHpXBgL4qePE9lb2In2QlxpCbEc/S0loAwkIMt08YyLeuyKF/onawOskNqwFBUfanG9vadV2xiFPqD7Xy+EfbePGLnYSGGOIiQznY0sFDU7P5h8Jcp+OJSwTFLWa6fjxuONwGeIt+zoLVJ/wDIOJPzW0d/P6TUib/6mOeX7KDGwuy+M03xtLc5qF/YhQvLdvNktIap2OKSwTFmX3Xlu875i0nLMTw0EvF/Pft47QrUBzR4bG8UVzGbxdtobyhmel5afxoRh41B1uYs2A1OWlxmksjfhcUZ/bgLfzk2Aha2j2MHZCobx7xO2stn2yu4tr/+ox/em0tqfGRvHzfROZ98yKGp8eztqyBubMKSIgOB46dSyPS04LizB68Szd1h1oJDYG/ldawpLRGhS9+s25vA/++cCN/21bLgL7RPH5rAdeOzjzmLk+aSyNOCoqy71qjz0mLo+FwG5WNzXx7fjFP3KalHOlZe+qa+E3RZv68Zh9JMeH8dGY+t00YRERY0PzQLEEiKMq+68fjxz7YSoiB8oZmbpswUGNbpcfsb2pl7kfb+OPSXRgDD03N5sGp2fSJCnc6mshJBUXZH/3jcVxkGClxkeyoaeKJ28Y5nEyCTXNbBy8s2ckTH2/jQEs7Xx+XxT8UDtedniTgBUXZH80Yw1X56by1Zi8t7R1EhoU6HUmCgMdj+fOavfymaAt79x9mam4qj8zIIy+jj9PRRM5KUJV91y64jzdV8fLy3SwpreWK3DSHU0lv99nWav79vU1sKG9kVP8+/OrrY5iUo+VB6V2Cquy7XJKdTGxEKEXrK1X2ct7W72vglws38dnWGrKSonnslguYOabfMVfYSPe4YUxBoAjKso8KD2VqbhofbKzkF55R+uaUc1JW38Rvi7bw5pq99IkK55+vHcEdlwzSkqD0akFZ9gCFI9N5t6ScNWX7GTcwyek40gs0NLXx359s47klOwG4f/JQHpqSQ0KMrrCR3i9oy35qbhphIYai9ZUqezmtlvYOXly6i8c/2kZjcxs3FnivsNEkSgkmQVv2CdHhTByaTNGGCh6Zked0HAlAHo/lrS/38euizZTVH2by8FQeuTqP/H66wkaCT9CWPXiXcv7lL+vZVnVQN2V2sZONwH7601LmfbaDygMt5Gf24cV7RnP5sFSHEor0vKAu+ytHeMt+0YZKlb2LdY3ATouPJDzUcP3cz1lb1kBKbAS/u/kCrh+rK2wk+AX1AI9+idGMyUqgaEOF01HEQZOyU3jslgvYXHmAkr2NlJQ1cNuEgXz+yDRuKOivohdXCOqyB7hqRDqrd++nqrHZ6SjiEGst75WUY6338/suH8IvvjaaqHBdSinuEfRlXzgyA4BFGysdTiJOefyjbby8fA8hBvonRvFa8V7dIUpcJ+jLfnh6HIOSYyhar7J3o1dX7OG3i7YQERbC8LQ4spJimDurgDkLVqvwxVWCvuyNMRTmp7O0tJYDzW1OxxE/+nhzFY++WcLg5BjmzR5PQkwEoDtEiTsFfdkDXJWfQWuHh8Vbqp2OIn6ytmw/D71UTF5GPO88fPkJl1VOyk45MhpbxA1cUfYXDkqib2yElnJcYlftIe5+fgXJcRE8d9dFxEUG9RXGImfFFd8FoSGGK0eksbCkgtZ2j24ZF8RqD7Yw+9nltHssL9x9MWnxUU5HCkiaNuk+Pm89Y8xQY8w8Y8xrp3vM3wrzMzjQ0s6yHbVORZAedri1g3teWEl5QzPzZo8nO1Ub6US6nFXZG2OeNcZUGWPWHff41caYzcaYbcaYRwCstduttfcc/ftO9pi/XTYshejwUC3lBKn2Dg/febmYtWX7+a9bC7hwUF+nI4kElLM9s38euProB4wxocATwAwgH7jVGJPv03Q+FBUeyuThKSzaUInHY52OIz5kreUnf1nPBxur+PlXR/GVzr0VIvJ3Z1X21tpPgbrjHr4Y2NZ51t4K/An46vkGMcbcb4xZaYxZWV3dM1fNFOZnUNHYTMleXXIXTOZ+tI2Xl+/moanZ3DFxkNNxRAJSd9bs+wN7jvq8DOhvjEk2xjwJFBhjHgU42WPHs9Y+ba0db60dn5raM9MHp+WlERpiWLRBSznB4tWVe/jNoi3cWNCff/pKrtNxRAJWd67GOdn0KGutrQUePO7BEx5zQlJsBBcNTqJoQwX/qGLo9T7eXMWjb5Rw+bAUfnnTGIzRQDORU+nOmX0ZMOCoz7OAfd2L0/MK8zPYUnmQHTWHnI4i3bC2bD/fnl9Mbno8v7/9Ql1OK3IG3fkOWQEMM8YMMcZEALcAb/kmVs+5Kj8dgEUae9xr7a5t4u7nV5AUE8Hz2jQlclbO9tLLl4GlQK4xpswYc4+1th2YA/wV2Ai8aq1d350wxpiZxpinGxp67g3UAX1jyM/so3X7Xqr2YAuznztq01QfbZoSORtnezXOrdbaTGttuLU2y1o7r/Px96y1w6212dbaX3Q3jLX2bWvt/QkJCWf+zd1wVX46K3fVU3OwpUePI143P7WUm59a2u3n6do0tW//YZ65c7zuPiZyDly50Fk4Mh1r4UPNuO81ujZNfVm2n8duKWD84N67aeqVBy7RuALxO1eWfX5mH/onRms3bS9xzKap60dy9ShtmhI5V64se2MMhSPT+WxbDYda2p2OI2fQtWnqW1OzufOSwU7HEemVXFn24L0Es7Xdw2dbNeM+kP3PUZumfqi9ESLnzbVlf9HgJBJjwrWUE8A+2VzFI2+UcFmONk2JdFdAlb0/Lr3sEhYawrS8ND7cVEVbh6fHjyfnpqSsgYeObJoap01TIt0UUN9B/rr0skthfgYNh9tYsfP4GW/ipN21Tdz1/PIjm6bio8KdjiTS6wVU2fvb5OEpRIaFaCkngNQdamX2c8tp67C8cPdF2jQl4iOuLvuYiDAuH+adcW+tZtw7zbtpagX79h9m3uzx5KTFOx1JJGi4uuzBu5Szd/9h1u9rdDqKq3k3Ta1mzZ7ev2lKJBC5vuynj0gjxKBZOQ6y1vIvb63ng42V/GymNk2J9ATXl31yXCTjB/WlSGXvmCc+3saCZbt5cEo2sycNdjqOSFAKqLL356WXR7sqP52N5Y3sqWvy63HFu2nq10Vb+Jo2TYn0qIAqe39fetmla8a9zu79a/GWah7t3DT1HzeNISREm6ZEeoru+gAMToklNz2eRRsquOeyIU7HcYV1exv41kurGObHTVOaNCluFlBn9k4qHJnO8h111B9qdTpK0NtT18Q3n1uhTVMifqSy73RVfjoeCx9uqnI6SlCrO9TK7GeX09bh4YW7LyJdm6ZE/EJl32l0/wQy+kRRtF73pu0ph1s7uPeFFZRp05SI36nsO3XNuP90azWHWzucjhMUnlxcypLSGsB7Lf3Df1pN8e79XDs6Q5umRPxMZX+UwvwMmts8fL6txukoQWFMVgJzFqymoamVnbVNLNpQSUxEKN8YP8DpaCKuE1Bl79R19l0mDO1LfFSYlnJ8ZFJ2CnNnFbCl6iBVB1qICg/hmdnjmZSd4nQ0EdcJqLJ36jr7LuFHzbjv8Ggwmi9sqThA10t532VDVfQiDgmosg8EhfkZ1B1qZdWueqej9HqvrtjDz97egAH6JUQxf/nuI2v4IuJfKvvjTMlNJSI0REs53fSXNXv54etrCQ81DE+PY0DfGObOKmDOgtUqfBEHqOyPExcZxqScZIo04/68Fa2v4B9e/ZIBSdH84Y7xJMZEAH9fw19b5sx7MiJuprI/icL8DHbXNbG58oDTUXqdxVuqmbNgNaP7J7Dwe5OZmpd2zNcnZafw4JRsh9KJuJfK/iSuzE/DGHS7wnO0bHstD7y4kuy0OF6462LiIjV6SSRQqOxPIi0+ioIBibqhyTlYvbueu59fQf/EaF6852ISYjTvRiSQqOxP4ar8DEr2NrBv/2GnowS8Dfsamf3scpLjIpl/70RS4iKdjiQixwmosnd6U9XRCkd6Z9zr7P70tlUd4I55y4iLDGP+vRPISNBgM5FAFFBl7/SmqqNlp8aRnRpL0QZdgnkqu2oPcdszyzDG8NK9ExjQN8bpSCJyCgFV9oGmcGQGy7bX0dDU5nSUgLNv/2Fm/WEZLe0e5t87gaGpcU5HEpHTUNmfRmF+Ou0ey8ebNeP+aNUHWrj9mWU0Hm7jxbsnkJuhUcUigU5lfxpjsxJJi4/UUs5R6g+1cvszyyhvaOa5uy5idJbzS24icmYq+9MICTFcmZ/OJ5uraW7rfTPub35qKTc/tdRnz9fY3Madzy5nR+0hnpk9XjPpRXoRlf0ZFOan09TawdLSWqejOKqptZ27n1vBxvJGnrx9HJfmaHqlSG+isj+DS7KTiYsMc/VSTnNbB/f/cRXFu+t57JYCpuWlOx1JRM6Ryv4MIsNCmZKbyqINla6ccd/a7uHb84v5fFsNv/r6WK4dk+l0JBE5Dyr7s1CYn07NwVbW7HHXjPsOj+X7r6zhw01V/NsNo7jpwiynI4nIeVLZn4Ur8tIIDzUUuWg3rcdj+eFra3m3pJwfXzOC2ycOcjqSiHSDyv4s9IkKZ+LQZIrWu2PGvbWWf3lrHa8Xl/H9K4dz3+ShTkcSkW4KqLIPpNk4xyvMT2dHzSFKqw86HaVHWWv594WbeOmL3TwwZSgPT89xOpKI+EBAlX0gzcY53pX53itQ/hrkM+4f+3ArT3+6nTsvGcQjV+dhjHE6koj4QECVfSDLTIhmbFZCUE/BfPrTUn73wVa+fmEWP5s5UkUvEkRU9uegcGQGa/bsp7Kx2ekoPvfi0p383/c2cd2YTP7jpjGEhKjoRYKJyv4cXJUfnDPu/2flHn7yl/VcOSKd/3fzBYSq6EWCjsr+HAxLi2NwckxQXYL5ztp9/Oj1tVw+LIW5swoID9VfCZFgpO/sc2CMoXBkBktLazjQ3Ptn3H+woZLv/WkNFw5K4qk7LiQqPNTpSCLSQ1T256gwP522Dssnm6udjtItn2+t4aH5xYzs14dnv3kRMRFhTkcSkR6ksj9HBQOTSI6N6NVLOSt21nHfH1cyNDWWF+6+mPiocKcjiUgPU9mfo9AQw5Uj0vlkUxWt7R6n45yztWX7ueu5FWQmRvHiPRNIjIlwOpKI+IHK/jwUjkznQEs7X2zvXTPuN1U0cuezy0mKDWf+vRNIjY90OpKI+InK/jxcmpNCTERor5pxX1p9kNufWUZUWCgL7p1IZkK005FExI9U9uchKjyUycO8M+49vWDG/Z66Jm77wzIA5t83gQF9YxxOJCL+prI/T4Uj06lsbGHt3sAb2vbk4lKWlNYA3puPzHrmCw40t3H92H5kp8Y5nE5EnKCyP0/T8tIIDTEsCsClnDFZCcxZsJq6Qy1srGik+kCL943lfN1OUMStAqrsA3nE8fESYyKYMKQvRQE4BXPikGS+fmEWW6sO0dzmITw0hCfvuJBJ2c7cJPyVBy7hlQcuceTYIuIVUGUfyCOOT6YwP52tVQfZHkAz7teW7efG3y/h6U+3Ex7qnXFz16TBjhW9iASGgCr73ubKABqMVneolUffKOGrT/yNsvrDfGtqNtZC/8QoXlq2+8gavoi4k8q+G7KSYhjZr4+jZd/hsbz0xS6m/eYTXl25h3suHcIvbxrNKyv2kJMWR1ZSDHNnFTBnwWoVvoiLqey7qTA/g1W766k+0OL3Y6/aVcf1cz/nn/+8jhEZfVj43cv55+vy2VZ1kLmzCkiI9o5BmJTtnWi5tizw3wsRkZ6hsu+mwpHpWAsfbvz72f3NTy3l5qeW9tgxqw+08INXv+Sm3y+l9mArc2cVsOC+CQxPjwfgwSnZJ6zRT8pO4cEp2T2WSUQCm0YddlNeRjxZSdEUbajklosH9uix2jo8/HHpLn63aAvN7R08NDWbb1+RQ2yk/hhF5PTUEt1kjKEwP4OXlu3iUEt7jxXv0tJafvbWejZXHmDK8FR+OjOfodogJSJnScs4PlA4Mp3Wdg+fbvH9jPvyhsN85+XV3PqHLzjU2s7Td1zI83ddpKIXkXOiM3sfGD8oiaSYcIo2VDJjdKZPnrOlvYNnP9/J4x9tpcNj+d6Vw3hwSrbuJiUi50Vl7wNhoSFMy0tn0YYK2jq6P+N+8ZZqfv7WerbXHKIwP52fXJev4WUi0i0qex8pHJnO68VlLN9Rd97PsaeuiX99ZwNFGyoZkhLL83ddxNTcNB+mFBG3Utn7yORhqUSFh5zXBqvmtg6eXFzK7z8pJcQYfnh1LvdcNoTIMC3ZiIhvqOx9JDoilMuHpVK0voKspGiMMWf8f6y1LNpQyf95ZwNl9Ye5bkwmP752hG4sIiI+p7L3oavy01m0oZK+sRFnvARze/VBfv72BhZvqWZYWhwL7pugYWUi0mNU9j40PS+NEAP1Ta2nLPum1nYe/2gbz3y2naiwUH5yXT53XjKI8FBdBSsiPUdl70PJcZGMH9yXdXsbyEo69mvWWt4tKecX726kvKGZm8Zl8aMZuaTFRzkTVkRcRWXvY4X56SzfUUdzW8eRx7ZUHuCnf1nP0u215Gf24fFbCxg/uK+DKUXEbVT2PvTk4lIy+njP1Oub2mhsbuNHr63l/fUV9IkK519vGMWsiwcSGnLmN29FRHxJZe9DXfd+jQgNobKxmct++TGNzW1My0vj198YS9/YCKcjiohL6V1BH+qaG9/u8dDS7uFQazv/dsNInv3mRSp6EXGUyt7HJmWnHFnKeWhKNrdPHOxsIBERAqzsjTEzjTFPNzT03jsqLSmtoepAC/0To5i/XPd+FZHAEFBlb61921p7f0JCgtNRzsuS0hrmLFite7+KSMAJqLLv7daWNejeryISkFT2PqR7v4pIoFLZi4i4gMpeRMQFVPYiIi6gshcRcQGVvYiIC6jsRURcQGUvIuICKnsRERdQ2YuIuIDKXkTEBXTzkiD2ygOXOB1BRAKEzuxFRFxAZS8i4gIqexERF1DZi4i4gN6g7QF6Y1REAo3O7EVEXEBlLyLiAip7EREXUNmLiLiAyl5ExAVU9iIiLqCyFxFxAZW9iIgLqOxFRFzAWGudznACY0w1sOs8//cUoMaHcc6XchxLOY4VCDkCIQMox/G6k2OQtTb1ZF8IyLLvDmPMSmvteOVQDuUI/AzK4b8cWsYREXEBlb2IiAsEY9k/7XSATspxLOU4ViDkCIQMoBzH65EcQbdmLyIiJwrGM3sRETmOyl5ExAWCpuyNMc8aY6qMMesczjHAGPOxMWajMWa9Mea7DuWIMsYsN8Z82Znj507k6MwSaoxZbYx5x8EMO40xJcaYNcaYlQ7mSDTGvGaM2dT5d8TvtzUzxuR2vg5d/zUaY77n7xydWb7f+fdznTHmZWNMlAMZvtt5/PX+fh1O1lvGmL7GmEXGmK2dvyb54lhBU/bA88DVTocA2oEfWGtHABOBbxtj8h3I0QJMs9aOBS4ArjbGTHQgB8B3gY0OHftoV1hrL3D4WurHgPettXnAWBx4Xay1mztfhwuAC4Em4E1/5zDG9AceBsZba0cBocAtfs4wCrgPuBjvn8d1xphhfozwPCf21iPAh9baYcCHnZ93W9CUvbX2U6AuAHKUW2uLOz8+gPebub8DOay19mDnp+Gd//n93XhjTBZwLfCMv48daIwxfYDJwDwAa22rtXa/s6mYDpRaa893x3p3hQHRxpgwIAbY5+fjjwC+sNY2WWvbgcXA1/x18FP01leBFzo/fgG4wRfHCpqyD0TGmMFAAbDMoeOHGmPWAFXAImutEzl+B/wQ8Dhw7KNZoMgYs8oYc79DGYYC1cBznctazxhjYh3K0uUW4GUnDmyt3Qv8GtgNlAMN1toiP8dYB0w2xiQbY2KAa4ABfs5wvHRrbTl4Tx6BNF88qcq+hxhj4oDXge9Za5B8W34AAAIQSURBVBudyGCt7ej8UT0LuLjzR1a/McZcB1RZa1f587incKm1dhwwA+/S2mQHMoQB44DfW2sLgEP46Ef082GMiQCuB/7HoeMn4T2LHQL0A2KNMbf7M4O1diPwH8Ai4H3gS7xLsUFHZd8DjDHheIt+vrX2DafzdC4VfIL/39O4FLjeGLMT+BMwzRjzkp8zAGCt3df5axXe9emLHYhRBpQd9RPWa3jL3ykzgGJrbaVDx78S2GGtrbbWtgFvAJP8HcJaO89aO85aOxnvkspWf2c4TqUxJhOg89cqXzypyt7HjDEG75rsRmvtbx3MkWqMSez8OBrvN9Ymf2aw1j5qrc2y1g7Gu1zwkbXWr2duAMaYWGNMfNfHQCHeH9/9ylpbAewxxuR2PjQd2ODvHEe5FYeWcDrtBiYaY2I6v2+m48Ab1saYtM5fBwI34uxrAvAWMLvz49nAX3zxpGG+eJJAYIx5GZgKpBhjyoCfWmvnORDlUuAOoKRzvRzgf1tr3/NzjkzgBWNMKN5/1F+11jp26aPD0oE3vX1CGLDAWvu+Q1m+A8zvXELZDtzlRIjO9emrgAecOD6AtXaZMeY1oBjv0slqnBlZ8LoxJhloA75tra3314FP1lvAL4FXjTH34P0H8Rs+OZbGJYiIBD8t44iIuIDKXkTEBVT2IiIuoLIXEXEBlb2IiAuo7EVEXEBlLyLiAv8fYSU6MWRCnY0AAAAASUVORK5CYII=\n",
Maximilian Schanner's avatar
Maximilian Schanner committed
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from corbass.evaluation import spectrum\n",
    "ls, mu_ps, a16_ps, a84_ps = spectrum(posterior, mu_coeffs,cov_coeffs,\n",
    "                                     pars.r_ref, r_at=pars.r_ref)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(ls, mu_ps, yerr=[mu_ps-a16_ps, a84_ps-mu_ps],\n",
134
    "            marker='x');\n",
Maximilian Schanner's avatar
Maximilian Schanner committed
135
    "\n",
136
137
    "ax.set_yscale('log');\n",
    "ax.set_xticks(ls);"
Maximilian Schanner's avatar
Maximilian Schanner committed
138
139
140
141
142
143
144
145
146
147
148
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we analyse the dipole. Using routines readily available in `corbass.evaluation`, we can get the pdf of the dipole intensity via sampling:"
   ]
  },
  {
   "cell_type": "code",
149
   "execution_count": 5,
Maximilian Schanner's avatar
Maximilian Schanner committed
150
151
   "metadata": {},
   "outputs": [
152
153
    {
     "data": {
Maximilian Schanner's avatar
Maximilian Schanner committed
154
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD4CAYAAAAQP7oXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXQc1Z3o8e9PkjdijHdjbLBsLCDeMEYYQyAkECc2YWLyksyBxwQmwzseQpjJmzmZFzN5vDfvvSwksyTDDMEDTIhNQtgSgoJNjDE2qzcZ2/IqW5Y3ScaWF+RFXtRdv/dHXXW32lJ3aa1efp9z+nTd23Wrfl229FPdqrpXVBVjjDEmiIKwAzDGGJM9LGkYY4wJzJKGMcaYwCxpGGOMCcyShjHGmMCKwg6guw0dOlSLi4vDDsMYY7LKunXrDqvqsOT6nE8axcXFlJeXhx2GMcZkFRHZ21q9dU8ZY4wJzJKGMcaYwCxpGGOMCcyShjHGmMAsaRhjjAnMkoYxxpjALGkYY4wJzJKGMcaYwCxpGGOMCSznnwg3JlMVz1sUW97z6BdDjMSY4OxMwxhjTGB2pmFMSPpxhsmym+HyMWxpggGj4OLJ0Ktv2KEZ06ZASUNEZgH/ChQCT6vqo0mfi/v8dqAR+HNV/TBVWxEZDLwAFAN7gD9V1WMJ27wM2Ar8g6r+k6u7Fvgl0A9YDHxbbZJzk21OHoLyZ/hV75foRcSve+9d/72oD1wxC675M+g/PLwYjWlD2u4pESkEHgdmAxOAu0VkQtJqs4ES95oLPBGg7TxgmaqWAMtcOdFPgdeT6p5w22/e16z0X9GYDKEKW34Pz98DlYvjCSNR5CxsfRVevA8q/9jzMRqTRpAzjelAlapWA4jI88Ac/LOAZnOAhe6v/lUiMlBERuKfRbTVdg7wGdd+AbAC+K5b706gGjjVvAO3vQGqutKVFwJ3cn5iMSbzNJ2Bt38Mu95qUb1bR1LpjebNfR6fLNjLKDnMlFEXQVMjrPgR8xYu5fnoZwGxi+UmIwRJGqOA/QnlGuD6AOuMStN2hKoeAFDVAyIyHEBEPoGfPGYC30naR00r+ziPiMzFPyPhsssuS/3tjOlmE+b9lkeKnmViwR4/IQAMHsf39kxjk46LrxhVpsoufn/RGmjw/6v/18JlNGpfyrwbQ4jcmPMFuXtKWqlLvo7Q1jpB2ib7P8BPVfVkB+LwK1WfVNVSVS0dNuy8iaeM6TlnT/L9Xr9gYsGeeN0n/wS+/B8tEwYAwgYdz5UbvsqC/cOoqG0A4C+KFnOV7OuxkI1JJciZRg1waUJ5NFAXcJ3eKdoeFJGR7ixjJHDI1V8PfFVEfgIMBDwROQP81rVPFYcxGaF43iIKifK/ixYytaA2Vv/tvTfy6u5PwuKlbbY9S2++H/k6Pyx6misL9lOA8u2i30Lkfijq3RPhG9OmIGcaa4ESERkrIr2Bu4CypHXKgHvFNwNocF1PqdqWAfe55fuAVwFU9WZVLVbVYuBnwA9V9d/d9k6IyAx3t9a9zW2MyUT/rXAxUwuqYuXHI3N41bspUNsmivhJ5C5O0weAUXIYtv+hW+I0pj3SJg1VjQAPAUuAbcCLqrpFRB4QkQfcaovxL1xXAU8BD6Zq69o8CswUkZ341y9a3Mbbhm8CT7v97MIugpsMVSrb+WLhqlj5N9FbWeJNb9c26hnIryO3xSs+fBaaTndViMZ0iOT6Yw6lpaVaXl4edhgmn5xp4J0ffJGB4l+WW+lN4EeR/0rrl+VS60WEJ3v9M0PkOABPR26nzPuU3Ullup2IrFPV0uR6G0bEmK5W/kwsYRzTC3k8cicdSRjgd1P5t9z67ihcheB1RZTGdIglDWO60tHd/sN5zhPRL3GcT3Rqkyu8qZzSfgBcLEeZJjs7tT1jOsOShjFdac1ToP6ZQIU3jlXeJzu9ybP0Zql3baw8u3BNp7dpTEdZ0jCmi3zm4aep+GAxFbUNKMJ/Rm+no91Syf4YvS62XFqwA840dMl2jWkvSxrGdJGvFL4bW17pTWC3juyybdcxlB3qP/JUgAfVK7ps28a0hyUNY7rC8QPcXFARK74c/XSX72J5dGq8sPONLt++MUFY0jCmK2x/jQI3qs1G73KqdHSaBu33njcJr7m766PNcPxAl+/DmHQsaRjTWZ7X4i//xV7yeJ5do4H+rPdK4hXVy7tlP8akYknDmM6qW+9PrASc0AtY613Zbbt6x5sSLyQNs25MT7CkYUxn7YhPlvSON4VIN86ivNr7JE3N2z+8Ez7en7qBMV3MkoYxnXGuEXa/Eyu+6U3r1t010pd13hXxCjvbMD3MkoYxnVG9AiJnANinw9mll3T7Lt/1JscLljRMD7OkYUxnJHRNLYtOo6se5ktljXcVFPX1C8f2wNHqbt+nMc26r/PVmFx3/AAV5X7XlIewwru6R3Z7lt4w5gbY5e6e2vUWDE6eBdCY7mFnGsZ0VMJZxnqvhGMM6Ll9j4uPfMuu5ZDjUxyYzGFJw5iO8DzYsSRWXNbNF8DPc9kM6HWBv9xQ499JZUwPCJQ0RGSWiFSKSJWIzGvlcxGRx9znFSIyLV1bERksIktFZKd7H+Tqp4vIBvfaKCJfTmizwm2r+fPhnfv6xnTQRxVwwn8i+5T2Y3UXjGbbLkV9oDhh6li7IG56SNqkISKFwOPAbGACcLeITEhabTZQ4l5zgScCtJ0HLFPVEmCZKwNsBkpVdSowC/gPEUm89nKPqk51r0Pt/cLGdImEs4x3vcnxZyd60uUJXVTVK6yLyvSIIGca04EqVa1W1XPA88CcpHXmAAvVtwoYKCIj07SdAyxwywuAOwFUtdHNLQ7QF7CfBJNZmk63GGV2mXdNj4dQPG8R45+oZ2XNOb/ixAE4tK3H4zD5J0jSGAUkPnZa4+qCrJOq7QhVPQDg3mNdTSJyvYhsATYBDyQkEYBnXNfUIyLS6v2NIjJXRMpFpLy+vj7AVzSmHXa/C02NANTqUCrdkOU9LUIRK72Ek37rojI9IEjSaO0Xc/Jf/22tE6Tt+SuorlbVicB1wMMi4m5K5x5VnQzc7F5fb6P9k6paqqqlw4YNS7c7Y9on4a6pt7xr6IlnM9rS4kG/6uX+BXpjulGQpFEDJP4pNRqoC7hOqrYHXRcW7v286xOqug04BUxy5Vr3fgJ4Dr/7y5iec+Ig1H3oL4vwVrTnu6YSbdTLod9Av3DqMOxbGWo8JvcFSRprgRIRGSsivYG7gLKkdcqAe91dVDOABtfllKptGXCfW74PeBXArVvklscAVwJ7RKRIRIa6+l7AHfgXzY3pOTvfiF9wHnUtR7go1HA8CuCK2fGKrb8PLxiTF9Le8qGqERF5CFgCFAK/UNUtIvKA+3w+sBi4HagCGoFvpGrrNv0o8KKI3A/sA77m6m8C5olIE+ABD6rqYRH5BLDEJYxC4E3gqU4fAWOCUoUdf6Si1p+f+5/3Dgk5IN/03/fnqd7HEZQprPFHvh0YznUWk/sC3SeoqovxE0Ni3fyEZQW+FbStqz8C3NZK/bPAs63UnwKuDRKvMd3i4Bb/QTrgNH1Y1dPPZrThEINY613J9ILtfsXW38ONfxVuUCZn2RPhxgRVGf/b573oJH8MqAyxOJowW2Dl6/6Q7cZ0A0saxgRxrjE+QCCw1CsNMZjzrdfx1OlQv3DuVIsEZ0xXsqRhTBDVK2LPZtToMLaH9GxGW5QCyqI3xCs2/85uvzXdwpKGMUFULootLo1eS5jPZrRlmTcN+lzoF47Xwr4Pwg3I5CRLGsakc2wvfOTu7i4odA/0ZZ6z9Iar7ohXbHopvGBMzrKkYUw6la/Hly+7gQb6hxdLOpP+C4j7sa7bAIerwo3H5BxLGsakEo20GDaEq74YXiwBFH9/LY/vG01FbQMVtQ387Gc/oHjeovQNjQnIkoYxqexfDaeP+csXDIFLr0+9fgZ4Nfqp2PKnCyoYyIkQozG5xpKGMW0onreIp5/+t9gT4FwxCwoKww0qgB16Kdu9ywAoIsrthWtCjsjkEksaxrRhICcoLagEoKK2gRtf6ZM1XT1l3o2x5dkFqyFyLsRoTC6xpGFMG24rWE+BG8l/s1dMHUNDjii4D7yJ1Ks/+u1FcspuvzVdxpKGMa1R5fOF5bHiGxn2BHg6HgX+cxvNdrwRXjAmp1jSMKY1desZKUcAaKQvK72JIQfUfiuiV8cL+1fB6Y/DC8bkDEsaxrRm+2uxxeXRqRk1OGFQdQyNXRDHi9p0sKZLWNIwJtmZBtj9TqyYbV1TiZZ7U+OFXcvCC8TkDEsaxiTbuRSiTf6ijmK3jgw5oI5735uI1zxO1sEtcOpIuAGZrBcoaYjILBGpFJEqEZnXyuciIo+5zytEZFq6tiIyWESWishO9z7I1U8XkQ3utVFEvpzQ5loR2eS29ZiIZN6ocSa7qcK2P8SKb0Sz9ywD4Dj92eoV+wVV2PNuqPGY7Jc2aYhIIfA4MBuYANwtIhOSVpsNlLjXXOCJAG3nActUtQRY5srgz/tdqqpTgVnAfzTPGe62OzdhX7Pa+4WNSengFji2B4Az9OZt7+rU62eB971J8UJCt5sxHRHkTGM6UKWq1ap6DngemJO0zhxgofpWAQNFZGSatnOABW55AXAngKo2qmrE1fcF/0Z5t70BqrrSTS+7sLmNMV1me/zhvfeikzlDnxCD6RotpqWtW+9fszGmg4IkjVHA/oRyjasLsk6qtiNU9QCAex/evJKIXC8iW4BNwAMuiYxy7VPF0dx+roiUi0h5fX19gK9oDP7sfNXx2fmWZPEF8ERHuAhGuFuG1YM974cbkMlqQZJGa9cNNOA6Qdqev4LqalWdCFwHPCwifduzLVV9UlVLVbV02LBh6XZnjG/Pe9B02l8eNIbKDJudr1PGfjq+bE+Hm04IkjRqgMSfntFAXcB1UrU96LqcmrueDiXvWFW3AaeASW5bo9PEYUzH7Ux4arrk82Ti7HwdNiY+FhU15TYWlemwIEljLVAiImNFpDdwF1CWtE4ZcK+7i2oG0OC6nFK1LQPuc8v3Aa8CuHWL3PIY4Epgj9veCRGZ4e6aure5jTGdduoI1K6Ll8d/LrxYusNFl8IA15vbdBoObAw3HpO1itKtoKoREXkIWAIUAr9Q1S0i8oD7fD6wGLgdqAIagW+kaus2/SjwoojcD+wDvubqbwLmiUgT4AEPquph99k3gV8C/YDX3cuYTimet4g5Be9xf5E/b8Zmr5i//8G6NK2yjAiMuQE2veyX962ES68LNyaTldImDQBVXYyfGBLr5icsK/CtoG1d/RHgtlbqnwWebWNb5fhdVcZ0qc8WbogtL8/QOcA7o3jeIq4Wj//Xy79zasqFH8CNf+UnE2PawZ4IN3nvEg4zTg4A0EQRH2Th4IRBbNFiTjffQnziAHy8N9yATFaypGHy3qcKN8eW13lXcIp+IUbTfSIUsd4bH6/Ytyq8YEzWsqRh8t6Mgq2x5Q+85MEOcku5dyXgz0T43EvPZ81MhCZzWNIw+e3EQUqkFvAnLlrrXRVyQN1rrUsaABMK9tKfxhCjMdnIkobJbwkD+G3wLs/ZrqlmDfRnp/qPOxWgXFNQFXJEJttY0jD5bffbscVsnJ2vIxLPNq4r2B5iJCYbWdIw+avxKHy0CQBFWg7sl8MSk0ap7ADPCzEak20saZj8tX+1P8cEsNUbQwP9Qw6oZ+zSSziqAwDoL6fh0JY0LYyJs6Rh8tf+NbHF8oS/vnOftDjbsFtvTXtY0jD5yYtCzdpYcZ2WhBhMz2uRNPbaqLcmOEsaJj8d2gZnTwBwRAewRy8OOaCetVEvp6l5FKGj1XDio3ADMlnDkobJT/vjXTIfeiXk1DDoAZylNxXeuHjFvpXhBWOyiiUNk58Sr2doPl3PiGtxHWevJQ0TjCUNk38aj0J9pb9cUMhG7/Jw4wlJi+sadevjsxYak4IlDZN/EidbGjGJRvqGF0uIDjGIfTrCL0TPQe2H4QZksoIlDZN/En85ji4NL44M0PLWW+uiMukFShoiMktEKkWkSkTmtfK5iMhj7vMKEZmWrq2IDBaRpSKy070PcvUzRWSdiGxy77cmtFnhtrXBvYZ37uubvFS3Pr58ybS218sD5yUN97CjMW1JmzREpBB4HJgNTADuFpHk8aNnAyXuNRd4IkDbecAyVS0BlrkywGHgT1R1Mv7c4cmz+N2jqlPd61B7vqwxHD/gT0AE0KsfDMvtUW3T2a6X8UHNOSpqG6jYsQuO2ACGJrUgZxrTgSpVrVbVc8DzwJykdeYAC9W3ChgoIiPTtJ0DLHDLC4A7AVR1varWufotQF8R6dPB72dMS4lnGRdPgcJAMx7nLI8CyvWKeIV1UZk0giSNUcD+hHKNqwuyTqq2I1T1AIB7b62r6SvAelU9m1D3jOuaekSk9QmORWSuiJSLSHl9fX3qb2fyS13C9YxLcm8u8I5oMYeI3Xpr0giSNFr7xZzc8dnWOkHatr5TkYnAj4G/TKi+x3Vb3exeX2+trao+qaqlqlo6bNiwILsz+UA1dqZRUdvA5144bjPXAeu98XjNP6r12/xbko1pQ5CkUQNcmlAeDdQFXCdV24OuCwv3Hrs+ISKjgVeAe1V1V3O9qta69xPAc/jdX8YE8qmHF1CxYxcVtQ000pddeknYIWWEk1zAVm+MX1D1R/81pg1BksZaoERExopIb+AuoCxpnTLgXncX1QygwXU5pWpbhn+hG/f+KoCIDAQWAQ+r6vvNOxCRIhEZ6pZ7AXcAm9v9jU3emlywO7a82StG7Y7zmHIbwNAElPanRlUjwEPAEmAb8KKqbhGRB0TkAbfaYqAaqAKeAh5M1da1eRSYKSI7gZmujFt/PPBI0q21fYAlIlIBbABq3b6MCeTqgthJKxV5+hR4W1YnTkC1fw00nQkvGJPRAt06oqqL8RNDYt38hGUFvhW0ras/AtzWSv33ge+3Ecq1QeI15jyexySJn2lU6NgQg8k8tQxjnw5nCmchcsbvohp3S9hhmQxk5+cmPxzbzUVyCoDjegF7m4fPMDGrvITHr3a/E14gJqNZ0jD5IeH5jM06zq5ntOJ9b5L/kF9tA6vffg0i58IOyWQg+8kx+SEhaWxMnEfCxOzWizmogwHox9mWAzsa41jSMLnP8+DAxlhxk13PaIPwvjcxXrQuKtMKSxom9x2tjk3t+rH2p0btgc+2rEy8rrHnXYhGwgvGZCRLGib3HdgQW/TPMvJratf22KGjOaID/MLZEy3O0IwBSxomHyReBPesayoVpaDl2cbuFaHFYjKTJQ2T2+x6Rru9702KF3a/6x9DYxxLGia32fWMdtuqY/hY+/uF08fg4KZwAzIZxZKGyW0J1zM22/WMQJQCViUOK1L9dnjBmIxjScPktoTrGZvsekZgLbuo3rEuKhNjScPkLrue0WGbdSz0dXdRnar359kwBksaJpclXM+g3yC7ntEOUQr5lx3DY8OKWBeVaWZJw+SuhOsZXDIVu57RPu8ldFEtXfwSxfNes5kOjSUNk8MSk8bIqeHFkaUqdByN9AVghBzjckmesNPkI0saJjd5HtQln2mY9ohQxGrvqlj5UwU2UaYJmDREZJaIVIpIlYjMa+VzEZHH3OcVIjItXVsRGSwiS0Vkp3sf5Opnisg6Ednk3m9NaHOtq69y+7P+BtO6I1UtrmcwcEy48WSp96PxLqpPFWwBNLxgTEZImzREpBB4HJgNTADuFpEJSavNBkrcay7wRIC284BlqloCLHNlgMPAn6jqZPy5w59N2M8TbvvN+5rVni9r8kjN2vjy6FKwvy86ZL2WcJo+AIyUI1wqh0KOyIQtyJnGdKBKVatV9RzwPDAnaZ05wEL1rQIGisjING3nAAvc8gLgTgBVXa+qzZ2nW4C+ItLHbW+Aqq5008subG5jzHlqywGoqG3gL95Uu4DbQU0Usd4riZVLZUeI0ZhMECRpjAL2J5RrXF2QdVK1HaGqBwDc+/BW9v0VYL2qnnXtatLEYQw0nYGP4kNfbPDGhxhM9iv3rogtX1tgSSPfBUkarZ3XJ3dstrVOkLat71RkIvBj4C/bEUdz27kiUi4i5fX19UF2Z3LJgY0QbQJgnw7nKANCDii7fZhwpjGxYC+cOxViNCZsQZJGDXBpQnk0kHzvXVvrpGp70HU54d5jnaUiMhp4BbhXVXcl7GN0mjgAUNUnVbVUVUuHDbMHuvJOwvWM9XaW0WlHGUC1jgSgkKhNA5vngiSNtUCJiIwVkd7AXUBZ0jplwL3uLqoZQIPrckrVtgz/Qjfu/VUAERkILAIeVtX3m3fgtndCRGa4u6bubW5jTAvuegbQoj/edFy5d2W8sG91eIGY0KVNGqoaAR4ClgDbgBdVdYuIPCAiD7jVFgPVQBXwFPBgqrauzaPATBHZCcx0Zdz644FHRGSDezVf7/gm8LTbzy7g9Q5/c5ObTh2Go7v95cJebNHiUMPJFS2Sxv7VoHbrbb4qCrKSqi7GTwyJdfMTlhX4VtC2rv4IcFsr9d8Hvt/GtsqBSa19ZgwA+1bGly+ezNmq3uHFkkMqdTQntR/95bQ/gOHRahhyedhhmRDYE+Emp/znr38dG2Tvrz/oF3Y4OUMpYL0mXB+y6xp5y5KGyR1NZ5haUBUrrkkYAsN0XoWXcGZR+2F4gZhQWdIwuaN2Hb2IAP6tth8xJOSAcstGb1y8cGAjRCPhBWNCY0nD5I69sZvt7CyjG3zEYA7pQL/Q1Aj128MNyITCkobJDZ4Hez+IFdckznFtuogkdVHZdY18ZEnD5Ib6bXD6GAAN+gkqdXSaBqYjKjShi8qSRl6ypGFyw67lscW13lWo/dfuFi2uaxza6o/zZfKK/WSZ7Od5UL0iVkycptR0rWMMgEHFfiHaBAdtYqZ8Y0nDZL+Dm/0HzoATegEb1R46604/3Ngv9iyMdVHlH0saJvvteiu2+L43kSiFIQaT+yrUntfIZ5Y0THbzPNj9dqz4njc5xGDywyZvLF7zTAWHK+HM8XADMj3KkobJbgc2QONRf7nfIDbp2HDjyQOn6McudfOfqULd+nADMj3KkobJbrveivWv/6jyYrtrqodsSHxeo866qPKJ/YSZ7BU5l3TXlHVN9ZSNiUmjprztFU3OsaRhsteed+HsCQA+0sFs0TEhB5Q/tukYmppnVmiogZOHUjcwOcOShslelfE5uN70ptH6NPKmOzRRxFYvIUnbrbd5I1DSEJFZIlIpIlUiMq+Vz0VEHnOfV4jItHRtRWSwiCwVkZ3ufZCrHyIiy0XkpIj8e9J+VrhtJc/oZ/LNyUOxaV0V4a3otDQNTFfb4I2PXU/612dfCDsc00PSJg0RKQQeB2YDE4C7RWRC0mqzgRL3mgs8EaDtPGCZqpYAy1wZ4AzwCPCdNkK6R1WnupedE+epv/vRP1JR8zEVtQ1s9C7nMBeFHVLe2ZgwDtU1UmVTwOaJIGca04EqVa1W1XPA88CcpHXmAAvVtwoYKCIj07SdAyxwywuAOwFU9ZSqvoefPIw5nyozC+LdIX7XlOlpu/QSTqo/O+JAOQnH9oQbkOkRQZLGKGB/QrnG1QVZJ1XbEap6AMC9B+1qesZ1TT0iItaJnY8ObORi8Z/NaKQvq7zkE1/TE5QCG/U2DwVJGq39Yk4+D21rnSBt2+MeVZ0M3OxeX29tJRGZKyLlIlJeX1/fid2ZjLTjj7HFd6JTOEevEIPJbxttCti8EyRp1ACXJpRHA3UB10nV9qDrwsK9p70+oaq17v0E8Bx+91dr6z2pqqWqWjps2LB0mzXZ5Fxji2czlnrXhheLaZk06tbbFLB5IEjSWAuUiMhYEekN3AWUJa1TBtzr7qKaATS4LqdUbcuA+9zyfcCrqYIQkSIRGeqWewF3ADYuc76pXgFNpwF/HvCdmtxTanpSHUNaTgH7UUW4AZluV5RuBVWNiMhDwBKgEPiFqm4RkQfc5/OBxcDtQBXQCHwjVVu36UeBF0XkfmAf8LXmfYrIHmAA0FtE7gQ+D+wFlriEUQi8CTzVua9vss6O+LMZy6L2bEb4hLXeVXyxcJVf3LcSRtmNCbksbdIAUNXF+IkhsW5+wrIC3wra1tUfAW5ro01xG6FYX0Q+a6iBA/5fsh7CW941IQdkANZ6V8aTxt4P4IZWfxWYHGFPhJvskfAEeLl3JQ30DzEY02yTjuNs880IDTXw8f7UDUxWs6RhsoPnwY4lseKbdgE8YzRRxAZvfLxi38rwgjHdzpKGyQ615bEpXek3kLXeleHGY1po8e+x9/3wAjHdzpKGyQ4JXVOMn2lTumaYtd6V0Pys7YGK+MRYJudY0jCZ78xx2PNebHC8W/7QL+yITJJjDICLp/gF9Vo8S2NyiyUNk9GK5y3iwf/7j1Ts87umqnQUe/XikKMyrbr8s/Hl6uXhxWG6lSUNk/FuK4wPT7HMhkDPXGNvAXG/Uj7aBCdtCJ9cZEnDZLTL5CAlUgtAhELe9qaEHJFp0wWD4RL37IyqdVHlKEsaJqMlDoG+0pvASS4IMRqTSvG8Rcz9YAAVtQ1+xc4lqRuYrGRJw2SuaITPFGyIFa1rKvO9702Kzx1+eCfU7wg3INPlLGmYzLVvJRfJKQCO6ADW6/g0DUzYTtGP972J8YrKReEFY7qFJQ2TuSrjQ5Yt86ah9t81K7wRLY0Xdr4JkbPhBWO6nP0Umsx0sh72rYoVrWsqe2zWsTDADVl/7iRUvx1uQKZLWdIwmalysf+QGFDhjeMAQ0IOyAQn/N3mUbGHMdnyStgBmS5kScNkHs+D7fG+8CXedSEGYzrijWgpkeahXg5thYNbww3IdBlLGibz1KyFkwcBOKEXsMqbEHJApr0a6M87ic/UbP5teMGYLmVJw2Se7a/FFt/ypsZv4TRZ5Q/RG+KF6uVw6nB4wZguEyhpiMgsEakUkSoRmdfK5yIij7nPK0RkWrq2IjJYRJaKyE73PsjVDxGR5SJyUkT+PWk/14rIJretx0TE5vrMNY1HWwytbV1T2WuXjmKrN8a/trH/KN/74Q/CDsl0gbRJQ0QKgceB2cAE4G4RSe4vmA2UuNdc4IkAbecBy1S1BFjmygBngEeA77QSzhNu+837mra1UNQAABLuSURBVBXoW5rssa0MvKi/fPEkanR4uPGYTvmDFz/bmFWwBiLnQozGdIUgZxrTgSpVrVbVc8DzwJykdeYAC9W3ChgoIiPTtJ0DLHDLC4A7AVT1lKq+h588Ytz2BqjqSjcn+cLmNiZHRJtg66uxu27uWz0q7IhMJ63yJnBEBwD4D2pWLQ05ItNZQZLGKCBx0t8aVxdknVRtR6jqAQD3nu5PylGufao4ABCRuSJSLiLl9fU20mbW2LU8NnnPUR3Q8slik5WiFLa8trHhOf/uOJO1giSN1q4baMB1grQNKvC2VPVJVS1V1dJhw4Z1cHemR6nC5pdjxUXe9UTsAnhOeN2bTiN9/UJDDex5N9yATKcESRo1wKUJ5dFAXcB1UrU96LqcmrueDgWIY3SaOEy2+mgT1FcC0EQRS6J2ATxXnKYvi6PT4xUbnvP/SDBZKUjSWAuUiMhYEekN3AWUJa1TBtzr7qKaATS4LqdUbcuA+9zyfcCrqYJw2zshIjPcXVP3pmtjssj6X8UWV0Sv5jifCDEY09XKojfGb52u3w5168MNyHRY2qShqhHgIWAJsA14UVW3iMgDIvKAW20xUA1UAU8BD6Zq69o8CswUkZ3ATFcGQET2AP8C/LmI1CTccfVN4Gm3n13A6x383iaTHNoO+1f7yyK84t0cbjymy33MhSyLXhOv2Pib8IIxnSKa46eJpaWlWl5eHnYYJpUl34M97/nLl99K8Zulqdc3WWkkR1g57pnYmGJ85WkYWhJuUKZNIrJOVc/7YbQnwk246nfEEwbAtK+HF4vpVgcYAuNuiVd8uDC8YEyHWdIw4VGFVT8HoKK2gfn7R1P8k20hB2W6063LEka/3f0OHK4KOyTTTpY0THj2rYpdEPUQfhWdGXJAprtV6yWsTByAct0z4QVjOsSShgmHF4XV82PFJdHp1Kg9U5MPfhO9NV7Y855/I4TJGpY0TDi2vkrF5o1U1DawuvYMzyX+IjE5bY+O5D1vcrzCzjayiiUN0/NOHYY1T8WKL0c/TQP9QwzI9LTfRD8LzYNU71sFH20ONyATmCUN0/NWPg5NjQDU6lB+H70p5IBMT9uvI2D85+IVq56wp8SzhCUN07P2r4Fdb8WKP4/MsUmW8tT1b4xhfe1J/06qg5v9u6lMxrOkYXpO5Cy899NYcbk3lU06LsSATJgOMpjXotfHK9Y86Q+PbzKaJQ3Tc8qfoWLbNipqG1hZc45fRGaHHZEJ2QvRz3JK+/mFhhp/Ei6T0SxpmJ5xaBtUvBAr/jL6ebv4bTjJBbwYvSX2wN/7L/wznD0ZdlgmBUsapvtFzsGKR2NjDlV442zubxPzmncDh3QgABdKoz90uslYljRM91u/EI7tAeAsvfi3yJdpfU4tk4+aKOLZ6OfjFZteguMHwgvIpGRJw3Svwzth/a9jxQWRL3CQwSEGZDLRO95kdqqbvTl6rsVoASazWNIw3ScaadEtxcgpLPKuT93G5CWlgKcid8QrqldA3YbQ4jFts6Rhus+6Z+BIFRW1DayrPcWNH5Si9l/OtGG7Xsbb3tXxig/+DTwvvIBMqwL9BIvILBGpFJEqEZnXyuciIo+5zytEZFq6tiIyWESWishO9z4o4bOH3fqVIvKFhPoVrm6Dew3v+Fc33apuA2yId0stjMykjqEhBmSywYLIFyivbfTvpqpYB5WLwg7JJEmbNESkEHgcmA1MAO5OmH612WygxL3mAk8EaDsPWKaqJcAyV8Z9fhcwEZgF/Nxtp9k9qjrVvQ61/yubbnfmOCz/QWxYiI3e5ZR5N4YclMkGh7mIl6Ofjlesfdpuwc0wQc40pgNVqlqtqueA54E5SevMARaqbxUwUERGpmk7B1jglhcAdybUP6+qZ1V1N/584NM7+P1MT1OF9/6FisqdVNQ28EHNOX4W+Yp1S5nAXonexGG9yC+c/thm+MswQX6SRwH7E8o1ri7IOqnajlDVAwDuvbmrKd3+nnFdU4+ISKv3bYrIXBEpF5Hy+vr6dN/PdKXNv4Vdy2PFf4t+mSNcFGJAJtucpTfPRGfFKza/DB/vb7uB6VFBkkZrv5iTh6Nsa50gbduzv3tUdTJws3u1OqG0qj6pqqWqWjpsmE3s02PqNvgj2Dp/jF7HSm9iiAGZbPWuN5mt3hi/4EVb/L8y4QqSNGqASxPKo4G6gOukanvQdWHh3puvT7TZRlVr3fsJ4Dms2ypznKyHN/8hdnvtTh3FU9E7Urcxpk3CU9EvJsy5sRKq3w43JAMESxprgRIRGSsivfEvUiePKlYG3OvuopoBNLgup1Rty4D73PJ9wKsJ9XeJSB8RGYt/cX2NiBSJyFAAEekF3AHYzC2Z4OxJeP1/wOljfrnfQH7UdI8NeW46ZZeOgqu+GK94/2f+TRYmVGl/qlU1IiIPAUuAQuAXqrpFRB5wn88HFgO341+0bgS+kaqt2/SjwIsicj+wD/iaa7NFRF4EtgIR4FuqGhWRTwBLXMIoBN4E4tO/mXBEm2DpI3C0moraBjwK+J9NX+WwXccwXWDi4rH8vJfHIDnBlFH4kzV95rthh5XXRHN8tqzS0lItLy8PO4zc5Hmw4oewcykAFbUN/DTyVZZ714QcmMklMwq28vdFv2bKKPeHyBd+AMU222N3E5F1qlqaXG/3QZqOUYVVP6dixcuxYa2fjc60hGG63CpvAu97k+IVKx61AQ1DZEnDtJ+qP8vappdiVX+MXsdL0VtCDMrksscjc1hWI/4fKNU18Ob/hqbTYYeVlyxpmPZb98sWcx584E1kfvRL2HDnpruc5AJ+HLmLKG5wiPpK/269aCTUuPKRJQ0TnCqsmu8nDWeNdxX/FPlTPPuvZLrZDr2UJyMJd1PtW+VfU7PE0aPsnkgTTLQJ3vkn2PFHKmobAPjQK+HHkbuJ2H8j00Ne965nSPQ4U1jvV1Qtg3ON8Ll/gF59wwwtb9ifhya9EwfhD9+GHX+MVa3xruKHEXsWw/S8X0U/BxMShr/btxL+8Ndw0sYv7QmWNEzbPA+2lsFv74eDW2LVb0av5YeRezhHrxCDM/lL4Ka/gWv+DPBv9a7YsIZ3f/xlm7ipB9ifieZ8ngd73oEPn4UjVfF6KYDr/5LHdvfHLnqbMBU/vBi4hNsLPs3cokUU4HGRnIJFfwul98PVd0OB/U3cHSxpmLjIOah6Ezb+Bj7eF7t2AXBQB/Mvka+yrfrCEAM0pqXF3gz2No3gu0XPM1BO+oMbrnkSasvhM38P/W3A0q5mScP4Y0dt+4M/BPWpwy0+OkcvXonexEvRW6w7ymSkLTqWv216kO8W/YYpuLGpaj+El78BN/13uPy2+MCHptNsGJF8drLeTxRby6CpESB2dnGaPrwWnUFZ9EYa6B9mlMYEUkiUuwvf4muFbyNofNiRS6bCDQ/B0JJwA8wybQ0jYmca+ehoNWx8we+K8iItuqE+1v6UeTfyenQ6p+gXYpDGtE+UQn4Vncl6bzzfKXox/kHdBvjtf4NLr4dJX4HR19n1jk6wpJEvohH/1sQtv4PaD1skCoAaHcYr0ZtY4U2122hNVtuiY/lW07fZOvmAP5Okm+OF/av91wVD4PJbofhTMGIyFNr/9/awo5WrVP371g9u9p+c3bcSzp44b7Wt3hh+F72ZtXqlzeNtckYjfSkuG8tl8nXuKlzOjQWbKXATgE4ZhT9u2qaXoPcnYHQpjCqFkVNg4Bi7/pGGXdPIBV4UGmrg2B44thuO7PKfq2g80uKMItbHKwXM33cJv4/exHa9LJyYjelBl3CY2YWruaWgwr/Lyon9TDTrexFcPBlGXg3DroIh46H3BT0cbWZo65qGJY1s4HnQdMo/czjxEZw4AMfr3HKdnzCiTa02bZE0riyB8TNhwhyKf5Dlx8SYDijAY6pUcX3BdkoLKhkmH7f4/LwkIgIXjYYhJX4CGXCJ/7pwJPS5MKfPSjqVNERkFvCv+DPmPa2qjyZ9Lu7z2/Fn7vtzVf0wVVsRGQy8ABQDe4A/VdVj7rOHgfuBKPDXqrrE1V8L/BLohz9b4Lc1zRfotqShCl7E/2Ude29KeI8klN3n0XP+cuSsf7dSU6M/bk7TKX+Y5xbLp9w6p1sdAjr5mkTif/bEO6B2eKPZpmNY7V3FLr0EeyjPmGbKpXKIa2UnEwv2MEH2cqE0trrmeckEoKgvXDAY+g323y8YDP0G+eV+g6DXBdCrnz8mVlE/KOwNBYVQ2AsKiqCgV0ZfkO9w0hCRQmAHMBOowZ/3+25V3Zqwzu3AX+EnjeuBf1XV61O1FZGfAEdV9VERmQcMUtXvisgE4DfAdOAS/Gldr3BTvq4Bvg2swk8aj6nq66ni71DS+PBZ2P3O+b/0k5NBD0tOFImO6AD26XD3GsFOHcVeHWHXKYwJSPAYLfVMkL1MKNjH5VLHaDkUuxbSrLU/0JLrg+9UXPIoaiWhtFIuaOdl6IGX+c+qdEBnbrmdDlSparXb0PPAHPw5vJvNARa6v/pXichAERmJfxbRVts5wGdc+wXACuC7rv55VT0L7BaRKmC6iOwBBqjqSrethcCdQMqk0SGnDsHhHW1+3NYv77b+M3WF0/ThYx3CRzqIgzqIj3QwB/HfD+gQGrERPo3pDKWA/TqC/TqCJd50AHrTxBg56BJIPRfLUfbVHOViOUZvWv7h2FYCSfn7QtXvgYie64ZvRLdMVBUkaYwC9ieUa/DPJtKtMypN2xGqegBAVQ+IyPCEba1qZVtNbjm5/jwiMheY64onRaSyrS/XA4YCh9OulZ/s2LTNjk3bevTY7MTv7sgSScfmHWB+R7c1prXKIEmjtU7w5D6tttYJ0jbo/gJvS1WfBJ5Ms58eISLlrZ3iGTs2qdixaZsdm7b1xLEJ0uFdA1yaUB4N1AVcJ1Xbg64LC/fePBh+qm2NThOHMcaYbhQkaawFSkRkrIj0Bu4CypLWKQPuFd8MoMF1PaVqWwbc55bvA15NqL9LRPqIyFigBFjjtndCRGa4u7XuTWhjjDGmB6TtnlLViIg8BCzBv232F6q6RUQecJ/Px7+T6XagCv+W22+kaus2/SjwoojcD+wDvubabBGRF/EvlkeAb6lq1LX5JvFbbl+nOy6Cd72M6CbLUHZs2mbHpm12bNrW7ccm5x/uM8YY03XsJn5jjDGBWdIwxhgTmCWNAESkr4isEZGNIrJFRP6Pq/9HEdkuIhUi8oqIDExo87CIVIlIpYh8IaH+WhHZ5D57zF3Ux134f8HVrxaR4p7+nh3R1rFJ+Pw7IqIiMjShLu+PjYj8lfv+W9zoCM31eX1sRGSqiKwSkQ0iUi4i0xPa5MWxaSYihSKyXkRec+XBIrJURHa690EJ6/bcsVFVe6V54T8j0t8t9wJWAzOAzwNFrv7HwI/d8gRgI9AHGAvsAgrdZ2uAG9w2Xwdmu/oHgflu+S7ghbC/d2eOjStfin8TxF5gqB2b2P+bz+I/L9bHfTbcjk3s2LyR8N1uB1bk27FJOEZ/CzwHvObKPwHmueV5Yf2+sTONANTXPJ5yL/dSVX1DVSOufhXx50hiQ6Go6m78u8qmi/88ygBVXan+v1bzUCjNbRa45ZeB25r/KshkbR0bV/4p8D9o+RCmHRv/LsBH1R8qB1VtfkbJjo3/GuDqLyL+LFbeHBsAERkNfBF4OqE68fssoOX37LFjY0kjIHequAH/IcSlqro6aZW/IH4LcKphVdoaCiXWxiWiBmBIV36H7tLasRGRLwG1qroxafW8PzbAFcDNrlvgbRG5zq1uxwb+O/CPIrIf+CfgYbd6Xh0b4Gf4f3B5CXUthl4CEode6rFjY0kjIFWNqupU/LOJ6SIyqfkzEfke/jMlv26uam0TKepTtcl4rRybKcD3gP/Vyur5fmwm4T8fNQi/O+bv8J9XEuzYTMI/C/sbVb0U+BvgP93qeXNsROQO4JCqrgvapJW6bjs2ljTaSVU/xh+RdxaAiNwH3AHc404BoWNDocTaiEgR/qn50W75Et0k4djMwe9b3Sj+6MSjgQ9F5GLs2MzC/z6/c100a/D/mhyKHZtZ+KND/M599BL+KNuQX8fmU8CX3M/O88CtIvIrunbopQ4fG0saAYjIMHF3RolIP+BzwHbxJ5j6LvAlVU2cvaUjQ6EkDqvyVeCthCSUsdo4NutVdbiqFqtqMf5/0Gmq+hF2bLYDvwdudfVXAL3xRya1Y+P/UrvFrXYr/iCzkEfHRlUfVtXR7mfnLvy4/4yuHXqp48cm6BXzfH4BU4D1QAWwGfhfrr4Kv19wg3vNT2jzPfy7GCpxdyy4+lK3jV3AvxN/Kr8v/l9WVfh3PIwL+3t35tgkrbMHd/eUHRsFP0n8ytV9CNxqxyZ2bG4C1uHfDbQauDbfjk3ScfoM8bunhgDL8BPpMmBwGMfGhhExxhgTmHVPGWOMCcyShjHGmMAsaRhjjAnMkoYxxpjALGkYY4wJzJKGMcaYwCxpGGOMCez/A62ZtaTofpbPAAAAAElFTkSuQmCC\n",
Maximilian Schanner's avatar
Maximilian Schanner committed
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from corbass.evaluation import intensity\n",
    "arr, pdf, samps = intensity(posterior, mu_coeffs[:, 0:3], cov_coeffs[:, 0:3, 0:3],\n",
    "                            r_ref=pars.r_ref, ret_samps=True)\n",
    "fig, ax = plt.subplots()\n",
170
171
    "ax.hist(samps, bins=100, density=True);\n",
    "ax.plot(arr, pdf, lw=3, alpha=0.8);"
Maximilian Schanner's avatar
Maximilian Schanner committed
172
173
174
175
176
177
178
179
180
181
182
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can also get the pdf of the dipole location:"
   ]
  },
  {
   "cell_type": "code",
183
   "execution_count": 6,
Maximilian Schanner's avatar
Maximilian Schanner committed
184
185
   "metadata": {},
   "outputs": [
186
187
    {
     "data": {
188
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADnCAYAAADl9EEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd3gc1dX/P9sl7a7qqqyqJVnNam6S5YJNDS0FMCTBpiRvEkiDvEDyJu8vCZCENEoSAqRQEwI2AZIXCEkIkIBxL+qy1XtfSau22j4zvz9Wu9JKsi2tJNsk/j7PPivtzty5MzvfOeee+z3nyiRJ4jzO4zzOPcjPdgfO4zzOY26cJ+d5nMc5ivPkPI/zOEdxnpzncR7nKJRnuwPnsXjIZDINoAP0gAiMAxZJklxntWPnsSjIThWtFdzi+VDuWcDw8DANDQ00NTXS0NhIR3s7Y+NjWCwWxsctTExYGB8fZ3x8HIvFgiRJ6PV6dDodLpcbm82KxWJBqVT6Ptfr9Gh1OvR6HTqdjrDQMFLT0sjMyCAjI5OVK1ei0+nO9qn/R0KhlMvm+vw8Oc8SbDYbTU1NNDY20NDYSGPD5HtjAw6Hg8yMTDIyPcRJSU4mNDQUnV6PTqtDr9dPkU6vR6PR+Nrdu+8DLtiyFUmSsNvtPgL7yDxhwTI+zsjICK1tbZ7j1zfQ1NxEZGQkGRmZHsJmZpKRkUFmRiZpaWmoVKqzeLX+vXGenGcZVquVAwcPsGfPHvbs2UNFRTkpKSlkTFouLyEyMzKJjY1FJpvz9zotvORcKERRpLOzk8amRhobGicfGg00NjZiMpkoKSlh27YL2bZ1G+vXrz9P1iXEeXKeYcwkY2VlBYWFq9m2bRvbtm1jY8lGtFrtKduQJAlRFHC6XLhdblwuJ67p724XLpcLt8s1uY2LkdERoiKjUCgVKBRKlAoFCqUShULh+1utUqHV6dCGaFEoFKc9F7PZzL59e3l/zx4+2LOH5pbm82RdQpwn5zJDkiQqKyv58//92WcZ50tGl8s16XaOMTb57na7AVAqlChVKlTTX0olKpUalUrp+0yp9LwfPHSA4qINCIKA2+1GEAQEwY3b7XkXBAGH08mExcKEdQJRFNFogtDptOi0OrQ6HTqtDo1Gc1LrfTKyXnTRxVy//XrS0tKW7Tr/O+I8OZcJ7e3t7H5pN7t378JisXDDDZ/kkksuYdPGTXOS0WazMTQ0xPj4GOPj49gd9snATSh6vZ5QvR6dXo9apQ6oPwt1ayVJwuFwYLFYsExYmJh8dzgcyOVytFodOq0WnU5PZGSk3/jWCy9Z3377bf705z+xcmUGO3fs4Prrb8BgMAR0Hv9JOE/OJcTw8DCvvvoKu3bv5sSJ42y/bjs7duxk06ZNyOX+U8eiKDI8bMY0MMDg4ABqlZoog4FQfSj6UD1BmqCAx5dzIdAx51wQBIGJiQksk9Fhs3kIQRCIjIzCYDAQFRmJQuE/G+dyuXj7nbfZtetF/v73v7N161Z23LiDj33s4wQHBy9Jv/7dcJ6ci4Tdbudvf/sru3bv4l//+heXXXYZO3fs5PLLr5hlTex2OwMDJkwmE5aJCSIjIoiJiSEqyoBSubxTy0tJzrngdrsxm80MDg4wZDajVCoxGAwYDAbCw8L9HjTj4+O89tr/8eKuFzl27Bif+PgnuHHHDi668KJ5jXX/U3CenAGir6+PXzz6C5599hkKCgrZuWMH1123nbCwMN82kiQxPDyMacDE4OAACoWSmOgYYmJi0Ol0S2oZT4flJudMOBwOBgcHGBwcZGR0hMiISJKTkwkLC/fbrre3lz/+8Y/s2vUi/aZ+7rjjTm6/7Xb0ev0Z6+u5ivPkXCBaWlp45JGHefmVl9lx4w7uuutuVqxY4ftekiQGBkx0dXczPj5GeHgEMdExGAyGsxq5PNPknA5JkhgYHKCzswObzUZCfCIJCQmo1f7j58rKSh56+EHeeecdbr/9i9zx1TuIjo4+K30+F3Aycp6X781AdXU1Dz74U95+522+8IXbOF5zgpiYGN/3DoeDzq5Oenq6iYiIJD0tjdDQsDNqHc9VyGQyj8cQHYPT6aCru5vDRw6h1WpJSkrGEGVAJpNRWFjIC394kaamJh555GFW5eZw086buOuuu0lOTj7bp3HO4LzlnMSBAwf46YM/obS0lDvv/Bq333a7z3X1uK1m2trbsVonSEpMIiEhcdnHj4HgbFrOuSBJEiOjI3R2djAyMkJcnJGkxCS/4FBPTw+/ePQXPPfcs3z8Yx/nG9/4H7Kzs89ir88szru1J8G7/3yXH/3oR3R0tPP1e77Orbd+xnfjuFwuuru76ezqRKfTsSIlhfDwiHPaSp5r5JwOt9tNb18vXZ2dyBVykpKSiY2J9QWHzGYzv/r1r3jiicfZsmUL3/5/32H16tVnudfLj/PknIG2tjbuvucujh8/zr3fvZdPferTPksoCAKtba309HSTEJ9IUlLSrHHTuYpzmZzTYbFY6OzqxGQysSIlhaSkZN801MTEBE8//TQPPvRTrvnENfzgBw8QGRl5lnu8fDgZOf/j8jntdjsP/PABNpQUs359EZUVVezceRNKpRJJkujs7GD/gX3IZTI2b9pCenr6h4aYHybodDpysnPYtHETDoeDffv30tnViSRJaLVavva1r1FTfRyFQkF+QR6/ffK3iKJ4trt9RvEfZTn/+re/cvfdd5Gfn88jD/+MlJQUwDMuMplMNDQ2YDAYWJm+8kOrFf2wWM6ZcDqdNLc0Mzg4QHraSoxGo2/4UF5ezmUfuRRDlIEXXniR9evXn+XeLi3+o93a1tZW7r7nLmpra/nFLx7lisuv8H03PDJMXV0dwcHBZGVmnfMqFkmSEAQBl8sjene6nJ6/nS5cbhetrS2kp6+cErpPe6lUakJCQmapmM4l2O12mpqbGBkZJmNlBjExngydw4cPc8WVl6PRaLju2uv4wQ8eICoq6mx3d0nwHzmVYrfbeejhh3jssV9y111389LuP/rUPBMTE9TV1yIIArmrcgkNDT3LvfUQz2az+QTwExMWXE4PAQVBQIYMCQmlYlLwrlZNCuBVqFUqNGqPWF2hUCAIAk6nc1L47nm5nE6sNiuiKBKkCUKr0/q0s1qtjqCgpZUSBoKgoCDycvOw2Ww0NjbQ3NJMXl4+GzZs4Hv3f5/f//53yORy8gvy+P73vs9//dfnzumHzWLwb2s5jxw5wi233kx+fj4PP/SIz4V1OBw0NDYwNjZKVlY2hqizI8yWJAmLxYLZPOSXiRIcFIw+NJRQfSg6nRaVWo1apUIuV8yLOPNxa31i9wkLExMTngyViQnsdjvIICQkxEfaiIhIQkJCzhppx8bGqKquIjYmhrS0dG6+5SZUShX//d938ZWvfpnx8XFef+2ND3UmzH+MWyuKIr949Bc89NCDPPH4E1x33Xbf5y2tLfT0dLMyPcNvTHOmIAgCg0ODmEwmzGYzOq2WqCgDoaGejJSlGOcudswpiiJWm5UJi0fwbjYPYbXaCA8LI8pgwBBlICgoaNH9XGifmpubMJlMZGRk8NGPfpSrP/pR3nnnHVRKJXX1dTz+2ONs3379Ge3XUuE/gpyDg4P81+c+y+DgELte3OWT201MTFBZVYHBEM3K9JVn1A2y2WyYTP30m0zY7TYMUQZiYmKJjIxcln4sR0BIFEVGx0YZGhxkcGgQl8tFREQkRqORyIjIM/aQGxsbo7q6Co1Gw7XXXctVV17Fr3/9G1544QXu/Nod7Nx5E488/MgZf3gsFv/25Ny7by8333wTn/rUp3ngBw+gUqmQJImOjnY6OjvIzysgPDz89A0tEqIoMjIyTL/JxODgIGq1alIEH3vaygdLgTMRrRUEAfOwme7ubsbGRok3xpOQmEhw0PIH00RRpKWlme6ebgoLV/PO229zx5138LNHfsYbf3mDxoZGdu3aTVZW1rL3Zanwb0tOQRD46YM/5YknHuepp57mqiuvAjzBoKqqSkK0WnKyc5Y9RWl0dJS29jZGRoYJD48gNiYWg2H5U8Rm4kxPpbjdbnp6e+jq6kShUJKUlOSn+lkujI+PU1VdydNPP80tt9zKRRdehCRJPPXUk9x737088vAj7Nx507L2Yanwbxmt7evr49bP3ILD4eTwoSMkJiYC0NPbQ1NTIzk5q4g2LF+2gyAI9Pb20NHRgVqjZkXKCgryC856xPNMQqlUkpyUTHJSsk/109jYQFSUgeSkJEJDw07fSADQ6/VsLNnEe++9x9DQIKOjo4SFhXHbbbezYUMJO3beyHvvvcejj/7yjHgsy4EPbQz6/T3vU7yhiJKSjbz7zrskJibicrkoryijr6+PkpKNy0ZMq9VKbV0t+/bvZcJqZc2ataxfV4TBEP0fRcyZ8Kp+LtiylWhDNA2NDezbv5fWtlacTueSH08ul7Nt24U88cQT1Byvob6+DlEUKSws5PChIwiCQMnGDdTV1S35sc8EPpSW85VXXubOr93J88//gcsuvQwAm81KaVkpqStSSUhIXJbj2mw26hvqsFqtrEhJJSsz6992jm0xkMvlxMbGEhsbi8PhoLu7i8OHDxEZGUl6+solDdhs2bwFs9mM3W5HrlBw5Ohh1q5Zh06n47nnfsezzz7DJZdezKuv/ImNGzcu2XHPBD50Y84nnniCBx/6KX95400KCgoAGB0dobKqkvy8AiIiIpb8mC6Xi6bmJgYHB8nKzCQ6OuactZDnqnxPkiR6e3tobm4mymAgPS19zmJhgeDXv/417+95jz++9DL9/X00NDSwZs0adDpPlYW/v/V3PvvZz/DMM89y9VVXL8kxlxIf+oCQJEnce9+9vPrqK/ztr38nNTUVgP7+fhoa61m7dh3akKUdW4iiSHt7Gx2dnaSmppKUmHTOktKLc5WcXkiSRHd3Fy2tLcREx5CWnh5wpUEvxsbGSEtPpaqymvj4eMbGxqioLCcnO4foaE+i/KFDh9h+/XX8+Ec/5pZbbl2KU1kynIycivvvv/+kO0midPIvzyDcbjdf+tIXOXjgAO+8/S4JCQmAJ+2rs6uDovXFSxrGlySJ3r5eKirL0YZoKSgoJCLi3M7j9KKjo52U5JSz3Y2TQiaTERoaRlJSMjabnZrj1bjcbsLDwwO+vhqNhvb2dhqbGtm2dRsajYa4OCPHjx/H7XYRHh5OUlISV115Fbd/8TacTiebNm46Z35PuVz2vbk+P+ctp81mY8fOG7Hb7bzy8qvodDokSaK2rha73U5hQeGShu3N5iHq6urQh4aSmZG5ZK7XckGSJNxut08AX15exrq165DLFSgUcuRyBcrJiu/nIgRBoLmlmYEBE/l5BQFrnGtqarjsI5fylzfe9GWtCIJATU01MrmMvNx85HI5nZ2drF23hk9+8pM89svHz4mYwYfSrR0eHuaaa68hKSmRZ595DrVajSAIlFeUo9fpyMzMWrKnn8UyTm1dLTJkZGfnnPUVt1wuF2PjY0xYJjzEczpxulw4nZ6/3YLbt61SqUI9KYLv7+8jJiYWURAQRAFREHG5XYiiiEqlQqvVog3RodV69LMhISHnBHG96p/omJiAVVxvvPE6t3/xdl595U9s3rwZ8Dy8WlpbGDCZSEtL56t3fIWmxia0Oi0rUlbwzDPPnvV83Q8dOXt6erjqqiu5+JJLePihh5HL5TgcDkrLjpGYmERy0tIUghJFkYaGeobMQ+Rk5xAZeebTkARBwGw2Yx42MzY2htVqRaVUog8NRafVodaoUavUPgKq1WoUirmF8KcaczqdTo/QfWKCCavn3Tq5JINKpfYQV6slPCyM8PCIM05ar/65r6+XvLx8wsMWruj6x9v/4NZbb2HXrt1cfNHFvs/37HmfyqpK2tvb+eEDP0IURXbsvBGHw+HzyM4WPlTkNJvNbLtwKzfeuIP//db/IpPJsFgslJWXkpO9asnKKI6Pj1NZVYkxzkhaWtoZHYM4HA5MJhP9pj6sViuRkVFERUYRGhZKSHDgWSCBBIQkScLlck0S18LwyAgjw8Oo1CqiJqu7h4WFnzEX0GIZp6q6isgIz5KEC31I7PlgD5/61CfZunUra9euw2q18vTTT/GbX/8GQ7SB4qINKJVK3G43X/jC5+nr7+P11944axb0Q0NOq9XKRy7/CJs2beLBnz4IeMhac7ya1YVrliTvUpIk2tra6OruorCg8IzkckqSxPj4OP2mfgZMJpBBTEwssTGxS1p4eimjtXa7naGhQQaHhhgZGcZgiCYpMemMXa/Wtla6u7vIy81f8BRZV1cXe/d+QGlZGUODg3z/+z8gKSmJ7p5uOjs6KCoqQqHwEPT6G7YTqg/ld7/7/VkZg34oyOlyudh+/XVERkTy7LPPIZfLGRkdoaqqkqKipYnI2u12Kqsq0Ov0ZGVlL6vrJggCQ+Yh+vv7PSliOh2xMbHExESjVi9PoGm5plJEUWRgYIDOrg7sdjsJ8QmTBaOXN2A2MTFBVXUVYaGhZGVlzVqbJRB46w6vX1eEQqHAarVyxZVXUFRUxMMPPXzGo7jnPDklSeLzn/8c/aZ+/u/Pr6FSqbBYxikrL2P9uiJCQkIWfQzTgIm6ulpW5azCsAzSPu+1tFonaG9vZ2BwAIPBQGxM3LKliM3EmZjndDoddHf30N3TTZBGQ1JSEtHRMct2ft7sovaOdgoLVvsthREo2tvbMA2YWLd2PXK5HLPZzEUXX8jNN9/C1+/5+uI7vQCc88L3b3/n25yoreXdd95FpVJhs1kpKy9jzeq1S0LMjo4Ouro72VBcsizTI6IoMjA4QHtbG4IokJK8guzsnFk3rCRJ58z8WqBQqzWkpqaSmprK2NgYXV2d1NXXEW2IJnEZ3F6ZTEZKygqiDAbKy8tYmb4SozF+UW2mpKxAEEXKKzz3WGRkJH99829s3XYBMdHR54RQ4Zwg56OPPsrrr7/Gnvc/QKvV4nA4OFZ6jPy8gkUvdCNJEg0N9YxbLGwoLllyN9bj7nkq94WFhZOdnY1eP/vmnO6h/DsQ1IvQ0FBWrcpFFEVMAybqG+pxuVxkZmQQNbn8wlJBp9VRsmEj5eVljI+Pk5GRuaj201LTaBIEKqsqWF24hsTERP765t+49LJLiDIYzrrU76y7tS+9tJtv/e+32PP+B6SkpOByuThy5DCZWVmLzioRRZHKqkrUajWrclYF9EPOvD7eNiRJYnRslLq6OjQazYIr9y0XOc8F+Z7FYqGxqQG73U5GRuaS12mSJIna2hPY7DYKC1YvOme2vqEeu93uS/c7dOgQ1153Df/359coKSlZol6fHOfkmPOdd9/h1ltv4e1/vENeXh6CIHDk6GFWpKRiNBoX1bbL5aK07BgxMbGkpS68+NPJroskgd1uo6GhHpvdTk52TsBjoOUg6LlATi8slnEaGhtxOhzk5OTMWhZwsejs6qS9vZ11a9cSHBz40EeSJOrqanELAnm5echkMv7+1t/5/Oc/xz/f/deyr9tyzpGzs7OTDSXFvPTSH9l6wVYkSeJY6TFiY2MXLTCw2WwcKz3GypUrMcYFRvKZ10WSQBDcNLc0TxaayiQ2xj87xbvLQji31AQ9l8jpxdj4GDU11URERJIZwLzlqeCdZlu3dv2ikqolSeL4CU+F+ZzsHACeffYZfvHoLzh44NCyJmyfU+R0u91cfMnFXH311Xzzf74JQFNzEy6nk5ycVYtqe3R0lMqqiiVJH/NeG0mS6OjsoL29jeSkyXU9ppNqDoLN/OhUxD1Xl51fSvjmlrs6WZWbS9QSKrG8v/lio/qSJFFaeoz4+ATi4z0Bp8985lZUKhVPPfX0UnV3Fs6ptVLu/979aLUhfOPr3wA8Tz+TqZ+srMW5D6YBE1VVlaxbu37BxJQkyffy/g8wMDDA/gP7sNlslGzYRErKCn9iMrcLLEkgihKSNEVM7+enO3agWIo2lgsymYzU1FTWrVtPU1Mj1TXVuN3u0+84D4SFhVGQX8ix0qPYbLZF9bGwcDXNzU2Mj48D8PjjT7D/wH5efPGFJenrgvpzpi3n2++8zec//zmOHS0lJiYGp9PJocMHKVpftKhxQ2dnB51dnaxftz6gifGZ18HlclNbdwKn00nuqlyCpgsgZrq8k+++YJH3ssn8P5+JUxlM7z6CIDBhncDpcOJ0OnA4nTidTpwOB06nE4fTidvtQoYMZJ5gjE6rQ5rslUwm9xSJDgkhJERLiDYEbYj2rFZ3lySJrq4uWttayM7KJiYmdknaHR42U1NTQ1FR8aKqLYyNjVFZVcHGkk0olUoqKyu5/IqP8MGevWRmZi5JX6fjnHBre3t7Kd5QxPPP/8FXLe1Y6VGSkpKJi40LuN2u7i66u7t8io9AMP06jIyMUF1TTUpyCklJSXgXt/LdyzNcVO+eXlL6bvoZl9z3+TQfdyY/vHVuh0eGsVgsyGQytFodGo0Gtdojeteo1ajVU/8rlUpf2zPdWkEQsNlsTFgnsE5YsVonsFqt2O12JCSCNEGEaLWE6vVERRnOaHV3u91OzfEaFAoFuatyl0TbOjQ0xIna4xQXbVjUfHZXdxf9/f2sXbMWmUzGb3/7G5588kn27z+w5HVxzzo5BUHgiiuvYMuWLdx3730AtLQ0Y7PbyV2VG3C7JlM/TU1NFBdvCDikPt2VbW5uwjRgoiB/NVqtdtJVnNpWJmOKjTP2nw6ZTOZHTs82Mvx+BpnMo7m1jNHf38/AgAmlUklsTCyRkVHo9foFq24WMuaUJAm7w451wuopGj00hNVqJTQ0lIT4eAyG6DOiaurt7aGhsZGMlRm+sd5iMDg4SG1dLRuKixclL6yuqUar1ZKWmoYkSdy449NEG6J57LHHF93H6Tjr5Hzghw/w3nvv8fY/3kahUDA8MsyJE8fZWLIp4BvAbDZz/EQNG4pLAn7qes/fmyeq0+omXZeZ48qpv6d/Iwgesyqf4/oKooRC4cfGycZEzGYz/aa+Sc2tntjYWAyGaDQa//NYqBVbbEBIkiRGR0fp6u5iaGhwUvWTuGwlLr1wOp3UHK9BpVSSm5u36IeCacBEQ0M9xUUbAr43BEHg8OFDZGdnExkZxejoKEXF6/nxj368pEs/nFVy7vlgDzt37uDI4aPEx8fjcrk4eOjAoqJr3joxReuLA162z3vuTqeT0lJPnmhiYpLPWsr8LN+0v6ddlunXbzpBhWnbKBQyRFGk39RPf18flgkLkRGRxMbGERkxqbmVzW4Dzjw5p8Orfurq6sJqsxEfH09CfMKyLXcgSRKtrS30m0ysXbN20TLL/v5+mpobKS7aEPA6NDabjSNHj7CheANBQUEcPXqUj3/iY+zfd2DJFk86a+QcGxujoDCf3/zmt1xx+RWecHXZMRLiEwLWR1qtVo4eO8q6tesCTpKdLlIvLSsjKyuLmOgYRNHr4nq2mzlMFCctpQTIJ93S6RAEEbli6qnvcDjo6u7AZOon2hBNbKwRvU6PTC5D4bUOM36a6QQ9m+ScDqfLSW+PR+yuUChITEgiLi52SbJEZsKboLC6cPWiLXZfXy8trS2+HM6A+mMy0dHZwfp1nvInjz76KC+/8jJ7P9i7JG7/WSPnPV+/h7HRUd88UWtbKxaLhfy8/IDac7vdHDx0gPz8goAy5b2QJMkT+KmuoqCg0KfymUlOz2eThBSnNLHer70E9bq3XgiiQGtbM8PDZhITk4mLNaKcHqya/DkU06ymF15yStLc7vKpcCbmOScmJujq7qKvr4+42FhSU9OWPFHZYrFQXlFGVmY2MTExi2qrp6eHzq4Oios2BBzsKisvJSEhkdiYWERRZOu2rXz2M5/hc5/7/KL6BmeJnNXV1Xzk8suoqqwmOjqa0dERqmtq2FiyMeCoakVlBVGRkSQtUkXkrW+6du06n2stSZMu6+Slmkm46e6s16S63cLkd96PJXp6u+jq7iQxIYW4WKPfuSqV/mRUTLOy3p/C33J63+d3U51JEYIoinR1ddHW3kpMTCwr01cu6dowTqeTI0cPk5WZ5StxGSjqG+qRRJHsSfXPQmG32zly9DCbN21GoVBSXl7O1R+9iprq40RGRi6qb2dchCBJEnd+7U7u/e59REdHIwgCVVVVrF69OmBidnV1IYkiiYlJi+pbR0cHzS3NFBdv8CPmVOdBFGYLCKbDYXf7iOnFkHmAI8cOMTFhY83qIoxx8aclleAW5xAqzCVqOPfEBXK5nOTkZLZsvoDg4GAOHNxPf3//krWvVqspWl9MXX09g4MDi2orMyOTsbEx+vr6Ato/KCiIpKRkGpuaAFizZg3Xb7+e7373O4vq16mwbOR86aXdWMbHue222wBobmkmISEBnTawMaLFYqG1tYX8RSwUJEkS9Q31mAb6/aJ4kjjFDplc5nNtp/bzb8ft9phJwS0huCUmJixUVB6jt6+H7Mx8kpPSEEU5TpeI0yVO7iPgdgvY7S5EUfK9YGocO7OvkijNIu25qAKSy+WkJKdQXLSBru4uyspKcTgcS9K2RqOhuKiY2rpaBgcHA25HJpOxevUaGhrrsVgsAbWxImUFQ4ODPvXQ9773fV5/43VKS0sD7tepsCzkHBsb45vf+ia//OVjvjIQ/f19rFjhqdLe0NDAm399c96L2wiCQEVlBQUFhYuay6yqqsTlcrFu7Xqf9fa3mNIsazjtK2CKmOAJktTVH6eu7gRJSWlkrsxFo5kdyZzZpts1m4yiICK4BSRRQhQkv7nUk0n+zjUEBQWxbu06EhISOHzkEJ2dnUvSTw9BN1Bbd4KhoaGA21Gr1RQWrKa8ojwg6aBMJiM3N4+a4zVIkkRERAQP/OAB7rjzDl9cYimxLOS8/3v3c8GWC3y5cLV1J8jOyqG+vp7t11/Htgu3cs01n+BXv/rVvNqrraslIT4h4NQsSZKorqkmKDiY3FW503IyJzeQyRAlCdErRpipMsBr+URfe339XVRVl6IPjWDlygL0uhkJ1tNuSqtt9o3gdom4XSLiSazhXNZ75n1+LhIUIDY2jo0lmxgdHeHIkcNMTEwsuk2NRkPR+mJOnDiO2Rw4QcPCwkhdsYKq6sqArl94eDg6nY7u7i4AbrrpZlwuJ8/97rmA+3QyLHlAqPnBMisAACAASURBVKamhssuu5Rdu3cTHByE0RjPgMnE+vVFfPrGT5GYmMT37v8enZ2dXHjRNqoqq08Zjevr653UzBYF7M7WN9TjdDp9uXrgdRE9309/6s1FCq/lE9wSLpeTxuZaNGoNRmMqcrnHAk8PHlkmnD7Jn16v9hPKhwQr/cLvbkEkOGjKG5gZmlcop/0/zfWejunX5VzLShkeNlNzvIZ4YzypqWmLnnqw2W0cO3qUvLw8IiICD8RUVlZgiI4mIT5hwfu6XC4OHDxAVmYmjU1NDA4M8JWvfiXg4NAZidZKksQll17C9duv58tf/jIjI8McPHSQuNg4cnNzKSou4ve/e57Vq1cDnmkW68QEv/71b+Zsz2bzzGeWbNgYcKi+rb2NoaFB1q5Z50dMz7tnm1OR0+nwEFMmh2HzMI3NtaxITkep9reUSoUcy4RzWjtT34WFeibT7Q6PBQ0PC5q1r0rluWm9N6+3j0rl1P9OlxOHw4Fr8t3hcOBwOnA6PP87XU6f8F2lVhESEkJIcAjB3vfgYDQazRkXvIuiSFNzEyaTidWFqxddwNlms3H02JFFpQU6nU4OHjrIpo2bFixQmJiYoKy8FJfTRUnJRkJCQrjjjq8iSRKPP/7EgvtyRsj58st/5MEHH+Tw4SMoFApaW1twOp3o9Xoamxp5/PHH+d1zv/clro6MjJCbt4q/vvk3H2G9EEWRQ4cPkZ2VHXCoure3h/aOdorWF08bY05X90weSxB90xvTySm4RQTB8//g0ADtHc3kZBUgSkpc08aRw8N2AIKDpyzgzCGIRuMfoZ5OUOW06RS1WjG5v4jZPMTAkImxsTEkSUStUqPRaHwvtUaDRq0hKCjIJ4zft38vWzZfgNPpxGqzYrPaJt+tWG1WX6BGowlCp9MSGRFJlMGw6JW+5oOxsVEqKivIysomdpGZKN4Hd0FBYcDz3V1dXQwPm8nPL5jX9m63m8amRoaGBsnOyqa5uZnMzEwiIiIZHh4mLz93znv5dFh2coqiSFJyEh/5yGU8+8xzOJ1ODh855JsXGh8f5447v8pnP/tZcrJX+VzZp556kl27d/Ovf/7L74leV1eLQqkkY2XGvE9yOiwTFsrLyigp2ej3ZJxJTm9U1AuXyz94IwgSvX3d9PR2k5mRh1vwv45eYoLHMkaEB01eD//+zCSnatIiarVqHzndgpuhoUGGRwawWCYlfnFxhIeFo5oWCFusWytJEg6HA4tlnCGzmYEBE0GaIBISEoiJiV3WWr5Op5Oy8lIMUQbS01cuyopbrVaOlR4NuKKiJEkcPnKIrMzsU1pgb4pbS2sLK1JSSE5O8axCMGGhoqKcTRs3I5fLefTRR9l/YB8v//GVBfVj2cn5xhuvc9fddxEWGsbHPv5xrt++nSiDgfhpEr0LL7qQb37zm8TGxuJ2u1iVk0twcDDr16/jgR/+0FftzGweoqGxgQ3FJQH9eKIocvDgAXLz8mY9Vaefr08NNPkuiJLftIYkSbS1tzA8MsLK9FXI5Qqfa+rpp41xi5NQvb/bGhEe5COneXgq+dcYN+XOecnpdrtwOMYYGjLhcNqJiDAQFxNHeHgYMpkM+TThvHfsuhxjzrHxMbq7uzGZTISHhZGQkEhUVNSyuMCiKFJdXYUmKIisRS5GNTA4QFNTIxuKSwIaz1osFioqpwg2E8PDZk7UniA8PILMjMxZLnB9Qz1qlZrU1FQmJiZYmZHOe/96f0F1h5aVnJIksXnLZu65+262bt3Gxz/xMe66625uuP4Gvwv/wx/9kOHhYR5+6GHMZjO1tSeIiIigqqqK3z75JO/96z0kSWL/gX2sWb024LotJ04cJygoiLS09Dn7Cv7uqyRKPqG6TzsrSdQ31AIScXHpvrEfQE/vuO/vccvUONPrkgL0dY8SHTdjXKqUMTJsR5IkNBo7ojiCKApERhgwxhkJmVz8d/qx1NMsrven8qqKZLLZyqGlyEoxm81093QxPDxMdHQMCfEJhIaGLilRJUmiqrrKV7lwMW03NTfhsNvJzc0LaP+GhnoUCiXp6VP3i81mo66uFpfLxapVq3yrZM+E2+1m/4H9bNm8BYVCwQ8e+AHtbW08/fQz8z7+siqE3t/zPiMjw1xzzbVER0fz+GOP88Ybr/PLX/7Sb7tLL7mUd995B4DIyEg2bdqMXq/HaDQSH2/k4MGDdHS0ExVlCJiYJlM/FouF1FNU3DtVBF2ukCMIAtU1FZOLsKYjk8kQBAlBkGhpG8HumHJ99TrPWK27fZjWRs8keV/3KAD9PaN+bY8M2xHFCVyuNlwuCyq1kRBtOnFxyUiomLD6z/uK4pR6aC5ies5laadTZDIZUVFRFOQXcsGWrURGRNLU1MjBgwfo7+9fsuPJZDIK8gtwOBwead0i2k1PS8fusNPV1RXY/ukr6enpxmq1IggCjY0NHDt2lPiEBIqLN5yUmABKpRJjXJxvauUrX/4Kr7/xOp2dnQH1ZTqWxHJeceXlfOqTn+Kzn/0vBMHNvv37SElewaWXXcI3vv4Nbr/9i4DnKaMP1TE4MORHPpfLxcsv/xGny0lKSgpbL9gWkNjAq3881Rhk+hTKzPGmIHpW26qoLCUuNh6DwYjd7j9H2d45RbigSatWd3xKshYcMuX2TLfOITo5bnc/IEepjEEf6nFxB/rHKSpK9NvHG90FT/RWNc0iy2Qy5HL/CgpnYirFarXS2NiA3eEgPy9/Sarww9JZULfbzYGD+wNOIRwYHKCurtYTO0lMYsWK1Hm7yd5SOxds2YpMJuMb//MNBLebn/3s5/Paf9ksZ2lpKXV1dezceRMAnV1dJMQnsGLFCv7x1tv8+Cc/5vnnfw9AXV0dSUlJs6yiSqXimmuupbGxEbdboLKqEqvVuqB+SJJEZWUFq3JWzTs4IM149jgcdkrLjpCctAKDwVNSMyjo5A8Ju0Pg6MF2v88GTTOlYSJWWy8TlnYUiihUqgRkMhUD/eMM9HvcY5vdjc3unjWNMxOe+1ZaFjXK6RASEkJh4WoyVq6ktOwYTc1NS9KPpbKgSqWSVTmrqDleveA2RkdHaWpqwu12k5aaRlpa+oLGr2q1GkOUgd7eHgDu+u+7+MMLf1iU3BCWgJw/ffAn3PXfd6NWq/EuOJOcnAJAeno6b/39H9zz9Xvo7e31rbM5l/shiiI52Tl88MEHJCclU1p2jPqG+nnLrJpbmgkNCzvlAkUzl0SYDrfbRUVFGelpmURFRaNSyWeRd3TM4ZsC6e0aobdrBIDxMbvfduVHuyk/2kV7SxsOZztul4LhwQi62210tp5e3TI65kAmlyGTy5CQcDrdOJ1u37SOp/+nds+XC5GRUWzetAVJFDlwcD/Dw+ZFt7lUBDUYolEpVfMWtzscDqqqKzlRe4JVOavYtHEzrW2tAUn70tLSaWltQZIk4uPj2X7ddh57/LEFtzMdiyJnXV0de/fu5fOf9+S09ff3ExkZ5ScYyM7O5qadN/Gzn/+M3Nxc7rzza9x8y80Igv+URWVVJb/57W/4xMc/QXR0NJs3bUGtVnPg4H66u7tP+YMNDw97SmtmZp10m1MRUxQFSsvLSEpKmaU6sUw4GR6xMzwyRUAvKadjfMzO+JidzvZRtDqRrFwXIVqoLZfTWGNnep5Y/Qn/zI22Nk97Lrfoeznsbhz22TfJdIJ6zuWkp7xskMvlZGRksmb1WhoaGqiursLlci2qzaUi6KpVq2hobDhlf0RRpKWlmcNHDhFtiKZkQwlhYWFoNBpWpKRS31C/4ON6qkPIeOqpp7j4kot5f8/7/OLnP2dsbCyg84BFkvOhhx7kK1/5qs9NbWltITU1ddZ2F110Ec8++wySJPE/3/gfgoODuOwjl3H06FEAhkeGaWlpJndVLkVFRZ6OyeWkrkilZMNGzMNmDh46wMjIbFK4XC6qa6pYXbjmlK7IqcYyJ2qPY4iKJibGvwKgN1Kq13seNqIo0d9/qowGCWOim/hkN23NSnq61AiCjLHJcarV5p5TZzs+amPIPLcbPxdJZxL0dO7wckGr1VJcvIGIyEgOHDxAT0/PogI7PoLa7TQ3NwXUhlqtIS0tjbr6ujm/7zf1s3//PgRBYPOmLRiN/ml9SUlJmM1mbPb517/t7e1lw4ZivvKVL6PWqLjrrrt49pnnKCgo4MmnngzoPGAR5DSZTLz+xut85ctfATwEU6vVaEP8x5N/+MPz3Hb7bfzuud9NBjPkvPH6X7jx059m+/XXceOOT1NdVUV5WTnbtm2bdRy1Wk1+Xj55efnU1ddSUVmB3T5lxU7UniA9beWCAhSyafOFvX09uN0CKSkrTrvfwIBHwO1yzpW5IqEJGsIx5qDxhAqH3XNpvcSs/GcLjQc66Ju0kh1tw569JonV0znKwODJBeIqlWJy6sTzmpnpcrZE8DKZjMSERDaWbKTf1E9lVeUsr2ih7eXnFzA4NEi/KbDc0MSERCYmLAwPD/s+s1jGOXLksGfR3KKiky5nL5PJSE9Lo7WlZV7Hcjgc3PDJG7jiyispKy2noKCQTRs3sWnTJh577HEee+yXAV+PgMn58ssvc/VVV/uUFS0tLaRPm1d0u93c8/V7+OGPfsg/3/0XH/vYx33fKZVKvvCF26g9UceaNWs5VnqMP778x1NavlB9KBuKS4iLjePI0cM0NTd5nnA2a8DlFG02Ky0tzWTn5M5pWUemubJeYnrhcgo+kspkIiG6IUbMMlqPuxntnHJlhltHGG71t/jd9YN01w9y8B+NtLWYaWuZGredjKCCW0Rwz6jMIM10cc9eloparWZ14WoiwsM5fOTwotxcuVzO2jVrqa+vwzKx8NxLmUxGXl4+x0/U4HA4OH6ihsqqSjIyMlizeu1pV0g3GuMZHBrE6Tx9TurXvnYnRmMc9917HwqFgrTUNNo7PEHC1atXExsTy/t73l/wOcAi1ud8cdeL3HfffYAnzO5w2ImIiMDtdlNaWsp9998LwIH9B0+qjdVqtXz0ox9FLpczODhEYeGpNYkymYy4uDiio6NpaW3xSK8WuISDbxpFFKmqqiQ7OxeV0jP9oVR6nqRutzAleD9Ne2ERagShm95OOSNDcvRxOsb7LIx2jjJYNxWtG24dYWLIRqhxnLiC2brSthYz8UlhKORyzGYbZrON2Bgt2hB/RYrgFpHJJexeAbzdjt1ux+F04LB7atCeOHGcoOBggoOCJt/PjODdu8itSq3m6LEjFK0vDrjqnVqtobBgNRUVFWzauPDyqdoQLRq1hj0f7GFVTg6rTvIAngsymYzUFam0trae8v5qaGjgL2/+hbrael//wsLCsFjGEQQ3CoWSG2/cwa4XX+SSiy9ZUP8hwHnOxsZGLrxoG+1tHSgUCj7Yu4fGxibefPMv7N27l+SkZLZffz3f+ua3Tjlf6Z2b2rL5ggVf/J6eHvpN/chlMux2OzmrVhE6x6K10zE9G6WhsR4kSE1dOWs762R2idU29fRvmXRHJVHyiQv6+83EJVgZ7A+m/YQFtdZzI7a+3zarzYkhzxgm1OiZ3/QSNC490ifR27zVXzgRG+MZItgdNpz2EQbNA4iCgEwuJygoiOCgIDRBQQRpgjwZJ0EaykpLyc3Nw+6wY7fZsNlt2G0e8kqShEqlIjYmFqMxPuCSovNBX18vzS3NFBUVL0pUP5d653QYHBqkrq6WqMgoBgYHWbd23YJFLaIosm//XjaWnDxr5Tvf/Q4Oh4OHHnzI7/P6+jr0ej3x8Qn09PRQUJhPZ0fXSa/3ki47v2v3Lj55wydRKpXcc8/dZGRmUFFewY4bd/Dkb5+ad7W0js4OEhMSF0xMQRBoam70iQ1GRkaoqakmVB9KZmbWrPSymQ+g4REzZvMQ69cVI5fL/AIsgltEPqNUiWViiqQyuYzurlHUGgFjgpX+3hCcDo/FdU646K/xHyd5STkTfVWe7eIzp1bb6usbI25S8ud22+nrH2J0dAilUkVUZAyF+Wv9bhS1Rumz7F69rUKhwGA4+YrSDoeDvr5eKiorkCQRozEeY5xxyWvRxsUZkcnkHDlyhOKi4oBT/tLTV3Lg4H7i4uJOSzCr1cqJ2hOAxJo1a9GGaIk1D1FzvHrBlffkcjnJScm0t7exco7kC0EQeOGFP/CXN96c9V18fAJ19bW+1crWrVvHm2/+hRtu+OS8jw8BjDklSWL37l3s2LGTqqoqOjo72Lx5C7/5zW/55Cc/NW9iepT+nSQmJtHa2rqgCe3Ork6Mxnif2CA8PJyNJZuIiIjk0OGDtLZNtTd7gSIXx4/XUFi4etZDwTumU6kUyOUydFq1j5gx0VM3Rmp6MLFGG33dU8QEfMRUBiuZGLIxMWRjrHv+YyZJctLb28bAQC2jo10EaYLIyiwgK7MAgyEOhULpm/+UyWVIEiw0UKvRaEhJWcHGko2sXbMOJE/Zx0OHDtLe3javcdZ8ERsbS2ZGJkeOHg64ppBnHZU8qmuqTjqmdrvd1NXXUVpWyoqUFaxfV+QLTEZGRqFRawIqEJaUlEx3T/ec9+a/3vsXsTGx5OfPLvGq1+t9+bUAO3fsZNfuXQs+/oLJeeTIEeRyOevWreOer9/D1VdfjdPhpKKigvLycsrKyigtLaW0tPSUNYJaWlvo7uqmpGQDefm5fP8H35/X8b1Ch5RJoYMXMpmMhIQENm3cjNPpZP+BfQwMDMzat6ammvS0dIKDpqK7CoUMhUI2K9vD4RR8SdBeuFzjyBXD9PXoEESVr4D0TIs5HXMRVCaXEb/WSF/LMHK5SFT0BJJkApS43XEU5K/BYIhDqZw95pyJQGdSgoKCSE1NZdPGzRQWrvYsP3DkMHV1tQjC0izPFxMTQ1ZmNkeOHgmYoJGRkWi1Orq6/cUr3lSu/Qf2ExwUxOZNmzEYZi9xv3JlBk1NC5+aUSgUGKIMcxL7+d//nptvvuWk+8bFGenr9YghrrnmWvbs2bNgxdCC3dpdu17kxht3YLfbcblc2Gw2vvAFjwhBJpP5XIdxyzg5OTn86dU/zxmy/utf3+Tw4cM8+NBD5K7KZdPmjawuLOSaa6495fEHBkxEhEec1E1SKpVkZWaRlJhEbV0tbe2t5GSvQqvV0t3TjVwux2j0lKaQy+W+p6LXjZ1JUMBH0KGBQRSKQQR3HIIwlZlS+Zr/nJpMIWO0c9yvrbFuC5IgMdoxTliinrAkj5g6JgEMhlEGhzSEhSd4VhXTaXzzoSHBs38iL0G9AaylQHBwMGlp6aSmptHe3sb+A/vJzMgiNjZ20YGk6OhoZDIZR44epmh9YMvzZWdlc/DQAWKiY9BoNAyPDFNbe4LQ0DA2btx4ynGtTqdDExTE0NAgUVGzyXsqxCck0NbWOmuZwsGhQYzGk6+MF2+Mp6q6kuTkZEJDQ7niiit49dVX+OIXvzTvYy/IcrpcLl559RV23LiDkJAQXvjDH7j0kkspKyunrKyc0tIyjh0r5dixUqqrahgbG+fb3/n2rHaGh4eRJIlrPnENl15yKUajkV//+jf88Ic/PG0fWttafVX8ToWQkBDWrV1HWmo6FZXlHD9eQ0tz85xpRXNN4guCiHJaLqUguFGqhlCqEoBpqWEtZoxr/H+kjoMejaUkTpa3FDyv6VDJneQWgT5KRUtLCCPDSuxWF1qdvy54LtGCF7Jpr6WaRZHJZKxYkcqG4hL6+no5euwoE9bFF+gyGAzkZK+itOxYQPN+KpWKrMwsqqorqagop6Ghnvy8AvJy8+YVcMpYmUFjU+OCjxseFs74+PgsSd91117Hn/78p5PuFxISgiiKvjn5HTt2smv37gUde0HkfOfdd0hLS/dFznr7ejEajfT09DA8POx3Amq1mpd2v8SvfvUEH+z9wK+d1tYW1CoN7733nu8znVZHUHAwkiRxxZWXc+edd8xapdhTL1SGXn/yFJ6ZiIryaEHHLeO43G56e/1VLNOtgndFsOkJ1+GTGSJdXc3I5RHIZGrGx+yYTRP0tczWlY50jqE3zghcTDM8waFKNt0Yw9qro+jpCaa3R40oTm0w2D/OTFhtblQqBQ6n4PNOVCrF7Ir0SzjNqdFoWL16Delp6ZSVldLQ2LAocQF4CJqYkMjxE8cXvK8gCIyPjzM0NERwcAjFRRsWdB/o9XqUSpWfMGE+kMlkxMbGzRJEXHPNtbz11lunTNCIN8b7xPCXf+RyGhsbaG1tnfexF0TOV195hU9/6tOAZxBut9vp7u4hOSWJlRnphGiDMcbHcdHFF/GlL32R2794Gxs2lPgWgAFPWpfNbmPz5s288+47vpN77fXXWJmezp4P9kyuetzKvffd63f81rYW0uaQB54OExMWJFHigi0XYLFYOHT4wGkF2xrNlDvpdHiSomWyUAb6PAKD9Kwo1DqN7xWZHkFkegS66NlKJa8KKGW1nuvuX0l3rZW3f9XDUPsElr45yDitUJhCIfNLIZsLp6qMsFh4H24KhYL9B/YxOjpbQrkQJCenYLfb5j3+kiSJ3r5e9h/Yh0wuZ9vWC+k39QX0oEhLTaWjo/30G85AQnwC3d3dfp9FR0dTVFTEW2/9/aT7GY3x9Pb2Ah7Lv/267bzy6vxLmMybnJIk8fY7b7N27RrAk9QcExNLZmamL/vEYXdSeqyM7373u+TnF7Bu3Xpef+11P2ldf38fRmM8eXl5bNm8hTVrV/PlL3+JN954nf/93//H7557ji98/jaee/Z3vPTSbvbt2weA0+lgdHT0lFknJ0NDYwNZ2dmo1WpyclZRkF9Ic0sz5RVl2Gw2v7HyzOrrTqeD7p42IiNXkJAwd93c8R5/cbMuOgRjgaefXmKuLAmj+AYjbz3WxUCnG11sCD1lnh/OS9AQnYaQSbd2bMw+Y21PD+wOt1+lhJnWczkgl8tJT0tn3br1VFZVYjYHnokik8nIzyvgRO3x02Z/jI2NcfjIIUwmExuKSzyBvOBgUlJW0BSA9jYiIpLRsdEFE1un0+F0OmcFOK/ffj0PP/II5eXlc+6n0WhQKBW+YcFll13G+++/N+e2c2He5GxtbcVqtbL9+u3s2vUinZ2dGOOMyGQykpKSaW5uQiaTER8fz8UXXcyXv/xlvvPt78zSvHpWpopDLpfz3HO/4+GHH2FkdIR9e/eTmZmJXq/HarPS3d1NcHAwVVVVALR3dPgKKy0ENrsNm81G5LRsE51Ox/p1RSTEJ1BWXkpTcyOCIPjmOL2QJInG5lqSEtPR6YJwu0V0+ikrlp4VhSr45AoYLzFzL4mi4Mpo3v1tD46JyQCUW0ScDOwkZBkIC5tqNz5+Skyh0849nvI+UGQy2clX9l1iaEO0FBUVc/xEzaLWLgkODiYlOeWk2R9Op4Pq6iqOH68hJ3sVhQWFfjm6SYlJmEz9814xwAuZTEZsTGxA67nEG40+F9WLW265lSuuuIKSjRu4+eabTrKfR4gAcMEFWzl48OC8pY3zJueePe9z1VVX8fprb/D0M89w9NhRNm4q4brt19Lc3MS111532jacLicut9uPsB/76MfY9eJu3/zo9TfcwM9+9ghXXHk53/n2d/jSl76EKIr09PSQmJB4sqZPio72dlKSV/hZRy9iYmLZtHEzapWKQ4cP0N/fiyRJyOWeagPdvV0EaUIIC4tkfHzqRtDpNX4kjUyPIjI9ipnQx2lZ87EYsrcZ+OeTvbjs/lbOuNZIQtZU9LC5wd/VGx72H3MrFTKUitnnAZyx3LHgoGCKizZQV19Hf39giwKBx70dHxvzGwOKokhLawsHDx0iKiqKkpKNc1b5l8vlrEhZQWvr/MTp05GQkOArKbIQGOPjfSQDT+3aR372CE899SSXXHIJt91++5z7xcXF+fJLIyMjSUtNm/faKgsg5x62bd1GcXExr77yKhdsuYA/vvQy26/bzgt/eHFe1QdM/f3Exp66XumWzVu497v3UV1Vwy233IpMJqO3t4fYmJgFl2wUBIG+/j6MRuNJt5HL5aSmplFcVMzI2DBlFUcZGx+bzGroJzHRM8bVTlqwlFQPCS3js+fsZhJ0zdVRJBWEcfj1YQTXFHnGeycwrvX0qeO46ZTnMGS2+kjpxUkFG2eIoBqNhuLiDTQ3N/vdsAuBN/uk5ng1giBgMpnYv38fbpeLLZs3Ex+fcEovKTExib7+/gBcVP2ki7qwOdfgoGBkMhk2mydG8vzzz/Pmm2/y97+9xd/++ncu2HLBnPupVCpUSqVvjnfbtm3s+WDPvI45L3JKksSeD/awbduFAIyMDBMZGUleXh47d97E5s2b53Ww3r5ejHEnnxsCz8TvHXfc4TeZ3NbeRsqKFfM6xnT09PYQF+e/PqZMJpvzHlarNeTl5pO7KpfGpjpKy4+SnpaBRqNGpfLsr9WqfSQ9GXKvzSG2IJbVl4cTHh/M+8/2IAr4RXAzr/bX83oJmpEd7bOeISEqQiZF7645hAeCIOB0OhgbG6Xf1E97RzsNjQ3Y7XY6OtoxmUyMjY/hcrmWJVtFrVJTXLyB9o52Ojs7AmpDq9ViiDJ4goDdXaxfX0RmZta8VsuWy+XExcbSF4D19ljB3gXvFx+fQPfkw8hqs7Jp0yby8k5f8S8iIgLzZABy27Zt7Jlnlsq8yNna6indkJHh0RiOjI4SHr6wMvhutxu7zX7KSmZzYcI6gUqpOm2az0xIkkR7e9ssJRFM5UR6Md0S6fWhxMUaCQ8Lp6mpge6e9lmWKj1zKiiVnDw1PvQGcFbkKNBGh9DZrkYXO0VKvVGL3qglJCrE7wUeYnrf++aI4LrcIja7g/7+HsorSjl89ABVNZV0dncwOjaKTC4jLCwMhUKBKEqYzUM0NzdzrPQYH+zdQ1l5KX19vYueDpkOpVJJcVERPb29s9Q7p4PL5eLEieMMDQ2hUChYmZ6+YCF+QmJiQBX3EuIT6OntPv2GM2AwGPjnP9/l/u/dj91un3cRuojISIYng2gLGXfOq/U9e95n27Zt57mVBQAAIABJREFUPjfDYhlf8HoXQ0ODGKIXHmnt6+sj7jTWdi4MD5vRhmhnqVH85zhne4KiKNDV3UHROk9B66aWFqqqj5KQkIokaXFOVoQ3JnqKVXtLlniJGRwsEB4ho2KvG12sjn5OH3wo2JyMQiH3i7y6XCIqlRxRFBgY6GPIbMLtdhIbG0dOziqfdlSu8H++NjU3smKGlyFJEmNjY/T09tDQ2IBep8cYH09MdMyiFxZSKJSsW7uOAwf3ExUZSXDwqZPeJUmis7OD1rZW0lLTyMlZhcViobqmio0lmxYU8NNpdQiTU3oLUR1pNBoUCgUTExMLylZ59dVXsFqt1NfX8Ze//IW77rp7XvtFhEfQ2NgA+I87vavwnQzzJKdnvAlQVVWF0+lCFMUFjQFHRkaIWKC1Bc/Uy7q160+/4Qy0trXOWVR6JmZaxZ6eXqINMb6nYlZmJpERsXR2tWCzdRIWnohKNfWET0yJ9Lm9IKFQ9VK1d2qKwLjWSG9ZL/nbVyGKEoONQ1h6x9EZ9azMnnLd+7pHiZucqhEEcbK4swmrtZ9oQywZ6ZlotTq/MpnzhUzmsaphYWFkZ2UzOjpKT083jY0N5GTnBDQ9NR1KpZLcVblUVVedMvtjaGiI2roTREVGsWnjZl+GjV6vR6fVYTKZThuTmImEhAR6errn9Vv775dId083mRmZ89peEAR+8tOf8POf/4Lrrt3ORRdeREbm/PZVqVSIooQgCCgUCt+483TkPO1jc+Z48yc/+TGvv/4ahugoLrzoQu75+j3s3r1r1iTtTIyMjhAePveCM2/+9U3e+sdbsz73DKJlC14Hw2azYbfbZz0M5hp7zbQcXd2dJCQk+XxfwS3+f+reOzqu+lr//kxRmZFGvffeZVnFqi6ADcGATW8mAUywKblAKMFcSoAklITEQICEFloCobiQG0I3xbZkFUuyLduS1Xuvozr1vH+MZqSRRqM5IvmtvM9aWhjpfE/fZ/dn46ZUkJSQSlRkDMNDrQwPtyEIi81DZ9cxJBIPpifBPXDOski/PGX2WBL84n3xi/fF1X1p31Ui0dLTU4eEGVZn5BAVGYNS6WY5f6lUam2Xi4BEIsHLy4uUlFRysnNobWujsqpyUTWWWJiIwN1tJvmnpqaorKqkuaWJzNVZJCenLOqRjI2No6VVfPQ1ODiE7h7xQamgwCB6e3vm9fgKHPjmAC+99BL377qf+3fdb5WH/fTTf6HRaFiTk8Po6Ag7d97C2Wed7fDxvDw9GZ0t4HDU71xWOBf6m4GBgWy79jqaGpt5+OGHCQgI4IUXX+TmHT9dch+CINg1PR795S+59tpr+O777/j6wNfk5eWSlZXJc8+tjL2sra2VKAc4gcyyKpVKkUqlqNVqnJ2dF5ns5mohFxd3kpMycXF2w2jswGgcRRAEAkM8kclmkEp16LVunH15qmWtX/TS1sKcxjWfjwGdrg9n52FiouOJiIi3sDSAqdDd6mNi5kKSLPajHYFCoSQnO4fw8HAqjpbT9AO5aJOTkmhrb7NUfen1es7Un6Gy6igR4RGsycld0ox0czN9fMzRUEfh7OyMwlWBWj22/MbzIJfLUalUFoF5992/sWPHzZw+fQpfX19KSkp4+eW50ZRFRWuJjY1j23XbGFrB8N6V+J3LmrU1NScwGAxcceXlZGfnEBAYgFKpxMfHh00bN7Fp4ybqamvJyspech/T09NLEnDV1dXRP9DPnj17ueSSiwkMDOTpp54mJCSUkiPF/M///IyG+kaHfRFBEOjr7yfBBk3mwvmcC2FqRYsCTGMP5pNomVkRnJxkxERHMjrmT2trMzMzLTg5+eGsGGNm0heVl+k645L8aKwbJCrKZC2Y6S/NiJ4VWm8/N0YGJ5FIjGi1Hfj7BxMeFm51vebuE6tBulbXbHtuiqMI8A/A18eXhsYGqqoqyczMWtGkMZlMTlpqOsdPHCc8LJzmliYiwiMoKlzrkG8bFhZOR2enw6bm3DpTYCglRdzk88DAIAb6+xkZHuG+X9zH5599YRnft+WiLZyz8WyuuuoqAgIC8PHx4dN/fcojv3yE6upqlAq3Zc3S+fDx9uFk90nTv3188Pf3p6WlhQQ7pvGyd6y+oYELL7iQq6+6mrGxMUJCQggPD7faRi6X024nnD46NmpzhqIgCLz62qtccfkVbDxnIxXlRzlxvIZLL72M5ORkXFxceO6550S9dGr1GCqVyqGXYb620em0qMfH8PEx5SolmARUKjNN+vL1VVqVzUmlMmJi4omKSkKr7cXGpHornzIqyouCddEUbYixCKYZqzIDCQgeJzIyhsBAa6rG+SwMYP1hkdpmt1gRZDIZSYlJ+Pn5c7SyYsX9nFKplMnJCdrb28jPKxA11iA4KNjK1HQU/v4BDAwOiNb6vj6+DA0Ncf0N17Pr/ges5mqmpKRw3XU/5uF5XVVyuZynnnyKuNg4fv2bX3Hbbbeyf/8+m5StC6FQKJienrZcW3x8PPWzQaKlsOxda6ivZ82aNVx11dU89uijJCenLCLsevLJp3jvvXcpKyuzuY+x0VE8FwhnT08Pl152Cd99+y133fVzywm7uLgwPT3NDTdej8pdxUUXblnuFK3Q29dHUKC46K5EAu0dHYSGhCGVSJAAxgUvyPwmZ/W8AgSdToubmzuRUXG4eYwiCIMIgv2XJDTSm9BIk4CGhbnR23sGH58IVCrT7zTaOcFYqnxvPkyzU/49ghoVFUVoSCjl5eWiGPRmZmY4fvwYZ87UkZO9BoPRIDptI5fL8fTwZFik2SiVSvFdoinaHlxcXOjo7MDDQ8Vdd9216O+/fOSXfPb5Z5SXl1v9vqioiIcffoSEhERe/8vrREVHkpySRExsNIFBAUTHRC26dolEgquri6XkMCE+gYZ6+y1sDmlOc1RqZNR2UMff359ndz/Ljp0326x3NAWD5kyOb7/7luycLDIyVlNaWmYV+h8dHWXzBZuJioxi1SrHJg7Px+DAAP52Uja2vsqCINDd3UVoyOLyQNm8VIWHysUyixPAxVnOuLqfwMAIfLz9iIxKB5yADpycp/D2U+I+O74hOHyx7xkUpqSvrx4/v2gUCmtyMoWrHA+Vi4XPaL5JKwiCiblhBX6mIwgLCycqKoryirJl61dNfE6NlFeUERgYRG5uHl5eXkRHxayoAyQ8PJyOTvETusJCwywFAo6itbWV8vJynn/ueZva3dPTkyefeJI777rTSit7eXljNBi4++67+dcnn9Lb08fePfv45sC3nKw5hdFotJl/VbgqLIG3+PgES3plKSyvORvqLT7A2Kht8xRgy5attLa2zvZczkEQTJO7nJ1NL7VWq+XWW2/h1Vdf4/HHHrdiNOjt7WXjxnNYnZHBHXfcsWg0wnIwGo0YjUa7dIy2NEz/wADe3t44u8wyu88TYL3OgGAULCakTmcg0N8NF2e5afCtdgalQoV6QmNqVI6OIjwijaBQF5ydewkIcCIzZ7HQC4KRqcl2/P1jcHW1DkB5eSyds5NKJJYhuv9JBAeHEBMdy4ma40t+0Hp7e02tXEBR4VqCgoIs9zc4OJjevl7Rpqa3tw9qtVr0vBJPT0/UarUok/jAga8xGAx2rY7rrvsxcrmct99+y/I7uVyOwTinGV1dXUlJSSEqKgp/f3/L3JSFMJm2poBXfEI89Q32NafdgNDY2BhjY2McOVLCeef9iNHRUWJjF1NJgol2pKCgAF9f6/rShYne1157DQCjwcDevXvo7x+gpbWF1tYWSktL2bnzFh568CGOHz9md8amLUxMjOOuWr44wlTCN/cQu7o6iYk25ckkUgky5tj35PMiqvp5I+m9PF3o6OzF2zsAiUSywPx0RiqNRKedZmSknampfgT8kUpNHw1BgJmZfvz9A/H29kIul6JWm0zl+YKp0xmsIrr/L4RyPoKDg+nv7zNZFfOaDtTjampPn8bV1ZXcNXk2o/AymQxfXz8GBgZE5S4lEgnBQcF093QTER4hap1SqWBqasrhwoILL7yI9z94n33791FVVUWAfwABgYEEBgQQEBBAQEAggYGBPPnkk2zbdi2XXnqZxXKUyWTo9XqbVULxcXHs37ePDes3WGlks98Js2btMprTrnA2NDQQFBTEG2++wc07bubZZ5+ltLSM+IR4IsIjiIyMRCqVUlpWyjXXmpqwa2trSU5OtuxjdMza3zSnZd548w3kcjn+fv5ERUeTk53NQw8+TEZGBmDiIBLT6Q4wNqbG08OxiJ2lf9NoZGJiAi8vT6tqoYX0mGASVN2sgAqCwMBgL6tX5SAIUsbm+aGCYPIVJ4CIiBQmJkYYHGxHpfJBq/PEYJhBJp3B13eOctHDwwXvWRNYP+vfus3mQi3nZZbNBcpBEIT/GGF0SkoqR0pL8PX1QyqVUl9/BvW4mpTk1CXz1mZEhEdwpr5OdGFBWFg41ceqRQknmHKtQ0NDDgtnUFAQf3/vfQ4XH2JcPU5zczMC0NHeTl9/H319/QwM9NPX18f09DTvvvsuP/uZafyIQqFkenra5jv6+OO/4pprr2HrxVt4+613LApLoVAwNmZK+ZhnskxOTuLhafs9tyucjY0NrFmzhr+/9z7qcTXff/8dVZVVHPjma9rb22lrM9WdZmVl84tf3E9KcrIlH2rG2OgogfMCNHfddZdN53s+9Ho9UolUdGnZmHpMdKnf+Pi4ZaS6RLJgHP2sVJhfe4NRsJzTyMgQ7u4eyOVO6HQG/HxMVUODw4uT+e7u3ri5eTIy0otW04LRaCAiIsVKoDzn+bJy+cKBuQvKDCWLPy4zM9NoNKaG4IHBARSurri5uf9goXVyciIpMZnSslIkEoiLjSM1Nc2h/Xp4eKDVakWX15lqbAW0Wq0ovltPD49lx/8dP36c559/jsamJtraWhkYGGD37t2UlpUhl8sYHhnms0+tC2IEQWByctKq9lepUDA1PWVTOIODg/n6q6956OGHyM0zyU9ubi4KhZKpWc1pbmBvbGwkKzvT5rnaFc7a2lpaWlqoqakhLj6OwMAgXnjhRauTNpfxCYLAP//5f1RXV5OTk2N5eGPqMauc4/Hjx3n6t0/RUN/AwOAAY2NjBAYGEhYWTlpqKs8//0fUajUeHvbZ221BrR6zOwbQ5ppx62PNFwRbr59ZaLp7OwkPiwIECzufIICvt+kBanUGPFTO6HRzPpeHKpr2dj1T0xP09zfj7x9FVMScVtHPNnvbKtGTSOZSJ4Ig0N/fT09PN6Njo3h4eKBwVeDs4mLS6AMDTE9PMTExgZ+vH6GhYctquaUwMDjAmfo6ZDIpEeGRVuatIwgPi6Cjo514kblLkxYcJDjY8Tk4SqUbk0tw+oyOjvLYY4/y4Ucfcv8vdrH9ppuIiowiJCSEsrJSbv7pDgwGAwmJ8VRWVpKdPZe3l0gkiwpTFEql3aoqJycnfvfb31FYUMDFl2zlid88wY03brcawhWfEE9DQ/3KhLO+oZ7IyEg2bjqHzz/7AtcFZXQSiQSZTEZ3dzc7b9lBV2cXM5oZ9Ho9V111NY/+8lH0er0lQPPM75/hueeeZdf9D3Dfvb/A398fDw8P+vr6+OKLL/jzy38CTL6urSZbexAEwepYjmJcrcZ/ARG2WUClMuki2hIwdeprNTP4+HhbBFins502cHWRMTM7d8VgNDAxOUZaag5TUxN0dTUzPT1IbHQ8zs7OJo3pJAOJxGYgRRBMfvWp06dQKBREhEewalWGlRbr6uokJdlULmg0GhkcHKChoR7Z7ORnRzXY5OQkp2tPIZPKyM7KQS6XUV5RTkREhChtHBISQsmRYuLi4kWt8/DwYFytxk4r7iK4urqimbEeZGw0Gvnb3/7Kgw89yNYtW6k5cXJRXEShVDI1NYWHhwc/v+tu/rD797z3rn2mPKVC4RAP0oUXXkRGxstUVVVx000/tUqzxccn2A0K2RXO9vYOfvv0b1m/bgN/+vOf+MV99y3a5tSpU5x73iZuueVW9u97ELlcTkVFBevWr+W+e+9DgoRTp07R3NxEa2sLqalp3HnnnVb78PLy4sUXX2Dr1osBGBsbFR0MmpqaEjUG0Ay1Wr1kkAvmuj7mC+nAQD/BwSHI5vHeOs+W+JlnaTo7ydDqDDg5yyyasLenm8CAINzdXFC4OhHgn0P/QB/Vx48SHBRCtJnyUxDmmfRz9bTNLc309fVY/L3lXnapVEpAQCABAYH09fdRXlFGfFyC3eZznU5HY2MDQ8PDJCcn4+sz9yIrlW6MjolrYJDL5Xh5ejE8PCSKM9bNzY3eHnE9l/MrwCQSCceOHeOOO+9Ap9Wyf9/HltmvC6FUKpmaNgmnUqlwKFJs9jmXw1NPP4VUKuXZZ58DrANJMdHRS9YGwDKpFLV6DC8vL2699VbkMhknak4s2uYvb/yFnTtv4dFfPoqTkxOS2a9+Wmoabm5K9AY9a9cV8cKLL/Lmm2/S2trCvn17LS91XV0dV1x5OV9+9SX3/+J+AMYnJsQHg9RjDgeDzBAEAY1WY7OwfuF7L5VJLWmM8XE1Hh6eGI1Gy2gEM8zC6OQsw9V17tsnQUJ3bxchs4TWcrl0ltMmiJysPAwGA6VlJQwMLGRGMFGmnKmvQ60eIz+vYEUmamBAIAX5hbS2ttBpI48oCALtHe2UHClG6eZGUWGRlWACRERE0N4uvrHa28fHEghxFG5KtxXx5bq4ujIxMcH9u+5n8wXnc/1PfkJxccmSggmgVJg0pyAIvPbaa9x8845lj6NQKCz+oz1IpVLS01dZLDqFQsnM7GBeT09PxuzUBNsVzomJCSRSCU89/RQGo4GxUesidIPBwJ49c3SZZhw6fIh169Yxo9EwODDIlou28OzuZ4mOjub55//ILx/9JUHBgWzctJGzzt5AXl4+x6qPW8wNQTBaNEdDQwN//OMfueeeu7n0sktYu24tTzz5BB0d1i+YemwMD5HCqdFocHVZ2sxbLKCmUr6JyQnc3VQ22eGXKqmbmBjH2ckJpZsSmVxq1Ycpk8lIiE8gKzOb7t5uKqsqLHMppVIJjU2N6LRaMlbNzXdZSaDHycmJ3Nxcurq7rDo5hoeHKSkpZmJigsKCIiKXIFLz9fFlbGxU9OxNd3d3xifEzdk0NY0bxbM4CAK33nYrtbWnqTlxkh07di5bJ6xUKpiemqKiwkRRs2njJgfPb/kKqMLCQkpKSiz/P1+o3VUqJuzcF7tm7dDQEGedtYELNl/Ajh07SEuzHtpy8NBBAvwDSEqam2G4Z89HPPfcs7z99jtoNBpa21rZsGEDzS3NREVFsfn8zWw+fzPd3d1UVJSzfv0GywBemDVJkHDy5EmefOoJvvnmGy679DLi4uNZt24dHp6e7N+3j+ycLHJzc7nzzrs4d9O5jKnHHGKCnw/1uBrVgsDT/JfSZB6Z/z23jV6vx1Xhsmi69EKYC9alUgmNfV2Ehi2sSZ4XkZVKcHVVkJG+mpHREU6ePIGXlxeBAYEMDgyQn19gk6RMLGQyOVmZ2RwpLUGpUNLc0oxerycjY/WyDfQSicTE4drd5VDXjxnubu5MihROMHWc6HRaSwHLcjh8+DCvvvYKhYUF3Hbr7Q5H+81m7auvvsJPf3qzw+tkUpmlR3Mp5K7J5fjsNHbT2Ma5XKe7uzvj40vflyXPQiKRSLRaLWWl5bz99js4u7hYzL+uri7uve9errrqSktaRKfTse26a3n0sUfZu2cfmzZuQjOjoaOjk08/+5Sf/ex2ts6bbh0SEsLFF1+Ct7c3Wq2WycnJ2WoiLU7OTpx73iYyMlbTUN/In/70Z+65+x4uvfQyNp6zkRdffInWljYuv+xy7r33HgoK8hlXj4vu+1Sr1cvO9Jy7H6aXU6vV4jIb3pfJpMjlMmRSqdWPVCoxFc3PsviBiXfJ12eBzzUrZ+b5nGYT2dvLm7zcAlTuKiqOllv4lP5duUypVIJK5UFp2RHCQsPIXZPrMLNFUFAwgwPialidnJzQ6fWitaCb0o2JCcdMW6PRyCWXXsz6dRvYeM5GUWk4V1cFbW1tfPvdt9y0/SaH180vKlgK7u7uJCUmUVpWumiNyl1l96Nl7wpcJBIJr7z6Crfddit1dbXceecdXHX1lazONBUKHKs+bpm0dObMGcrLyzlaMUe/oNFouGn7dlTuKt56621uueVWqwP09fXx6GOPEhkVQWBQACoPdzZu2khXZxdDQ0Pcd+99S/qeSqWS7dtv4vixE+Tn59PR2cH27TeKGvE97mDKZr7GMlG0qCy/X7wtFsE0w8QWL7HSlPY6SqQyqUnwnZyICI9EbzBQcqSYoSFxU6oWwlxDXFxSjKenJz7ePqJ5e8xaRixcXJxF88y6uTnud0qlUp7d/Sxvv/OWQ10iZgiCwGOPP8bQ0BDffvOd3brshVA4eC9+8pPrueSSi8nKyuTU6VOWiLJKpWJ8YjFflBn2hFPl6uqKj48Pq1ebTJ7CwkIu3noxp0/V8off/4HQ0FDLxpOTk/j5+Vs9bL1BT3BwMG+99TbnnH2O1c53P7ub1LQUBgcG+P67g0yMT9LfN8ATTzxJdXU1rq6uDvUUSqVS7rrr5zg7OzM4NMh333237BozJiYn7WqMhWakRGISTpXK3WbR+SIfdVZzTk7OlZRJpRLk8nlaVWJ9LPPxJBLo7GgnMjKSpMRkMjOzaGlpobLyqN35HEthbGyM0tIjDA0NkZ+XT0x0DJGRkbR3iCtON52jdEWUlBN2XkRbcBNpDv/kJ9eTm5vH8RPHqaurW3I7tVrNZ59/xoMPPUhhUSH/+uQTkpOTCQoSV8nk7OyMzoEPzh133MFA/yBXXnUVH+/fb7EgTGbt0vfEns/p7uPjw677dwFQcqSYrVsuXpJxbGp6cSpDuqCG1QytVsvTTz9F6ZEy4uLm0hhKpXI2fO8jagpwaGgosbGxXLD5Qk6LGJIzP/C0/Lam6xifmLD6KEkkLMmybha06elJVO4qK21p/pu5u2ThbZqZmcFgNFosBzelGzk5axgcHKSyqpIAf39iY+OWZYDTaDTUnaljenqK1NQ0K0vB3z+AujN1osv/3Gc1mqMuAYDK3Z2JiQnR6ZTOLnEdKs/87hn2f7yPs885iysuv4JHHvmlhbD8tdde5fXXX6fuTB05OWtYt24dv/nNbygqLOLkqZPo9HqH/VtYXKNtD05OTly89WJuuXUnwmx6TPUDAkIqlRWN5dyJ6HQ61Gq1Jbp65MgRHti1a1G4WiKRLOqLBNO0sqSkZCvBNEMqMWkTc43tcjAajZSWldLa2sITT/6G66+/wfI3QTCVgC3li0pEzC6YM2snUC2g91wUtV1wyeMT40u+yBKJ9X/NrAbjE+NWgTIz/Pz8KCosor2jneKSYmJjYwm1QcBsNBppaW2hq6uT+PgEggKDFm0jlUpRKpTMzMyIMm/dVSomxidECae7uzs9y5TWLYTJP5tZfsN5kEqlBPgHcLLmFE88+QTpq9K45557iY6K4oknn+Cdd/5KXm7eondCKpVYxmc4CjHCCZCUZCJXM1s+CoUCrVaLRCKRC4KwKLlqT224u80z+cwPVqvVsuncTURFRxIcEkReXi7XXHs1d955J3/84wuLTt4Wg/O+vXu58oorbR5UIpHanJdpC4cOHyImNpq77/45zs4ufPvNdzz04EM8/MjDhIWHonRT4Oau5DdP/Mah/TkCR7TtQpPX1JnjPs9kXfqjYA48TYwvTT9qHkdQkF/A6OgoR0pLGBkdmT0/gb6+PoqLDyMYjRQVrrXMtLEFU5pDnLnp7uZmSfU4CicnZ9EpGNNw45Xx7Pr6+rL7D7spPlzC0aMVbL9pOx9+8BHr161fIq8txbhMk/yi81tC+Sy5vVTK6tWrGRwYnD2mpSzQ5oO2rzlV84XT9DW+447/wc/Pl7FRNb29vTQ2NZKVmWXzRVpKc0qlUpycbB9aKpUuyyQAMDg4yI9/fB0v/PEFtmzZytcHvmb/x/v5059e4qwNZ/HtN98RERHB2NgY6zesIzQkhO0LInGCLW6RZSCXyzEYDBZz0t6XUzLPn5y7Pse09fTMNB7LlDA6OzuTlprG+Pg4p2tPIZFImJycpKe3mzVrHJsg7Wily3wIiG9f0+v1ODlIwjwfK4pQS+aqhOLi4vjowz1MTEzYjS9IJSvTnGKEE2DNmlxa21otVUIqlQq1Wq0CFkWx7KkA6fwbo9Pp2HnLDqqrq3n7rXeQSqWEhIRgNBr53/99gOt+vI0LLtzMFVdebnnYS6n9rKxsqqqqlrzg5W6SIAjs2HkzV155FS6urmRlZdLQUE9Hezt7PtrL22+/Y6E8CQgI4M9/fplXZ/tIbe1LDOQyOfp5/Dr2taB50K0TBoPeJmvBUrlLmUyO0cGgi4uriyl4Mjk121vq7nCNscFoP09nC0v1MdqDTq9DLrLuWafTWbEPOgwbj3S5VJFUKl5zijVrAX583Y+RSCRcf8NP0Ov15mdv8yWyJ5zj83NMR44cISwsnEOHDltd6JtvvMH09DQXXXgRd9xxJzqdjt3P7p49eanNk1+7di379u/j8V89vmgcmyM36eWXX6ars4vExERuumk7Tz71FGlpabz44ks2y7Ty8/Lp7u5aNEBGKpWJ7tSXy+XoddbuwfxI68Koq3mNLcIs60iw9Tq5XL5sjafRaKS1rZUjR47g5enFWRvOQqlUIpFIKS45TI8DZFkr0WgrEU69TvxxtDotTiJaxswQEN/falxBT+xKhNPV1ZW83DxGRkbYvv1Gc0DIpo9gTzgnzGHs8fFxpqenefihhxaZSn39fVx51VVce+02Np+/meeefZ7nn3+Orq4u9AY9ba2tvPnmG9x998+56abtTE1NkZ6ezl/f+RvHqqtJSU1mx46bOXny5NwF29CcExMTfPb5Z+x6YBePPf4o52/ezG9/+zTfHPiWzedvtntDFAoFL77wIjt37qCpqcnye5lkegBXAAAgAElEQVRMZkU34QhkcmvN6Qic5HJ0s19JW8Jr85xdXZmcXDrHNzg4SMmRYmZmZigqLCIsLMyy39jYWPJy8+nv76esvNQu9+/U5CSuIufQmIRTnEbT63Xita1Wh7NIbbvSoU1ajQYXEZFasKmgl4cEJFIp+/bup6+/z1xzbFM47d2tcXOgoKura9aJXvxC9fVas91FR0ezc+ctZGVnsmaNaUDtzMwMqzIyKK+o4K9/fYedO2/hscceZWBwAE9PT5pbmvnR+eeRnp7O44/9ykJmNTk5ye5nd/PFF19w/PgxIiMj8fb24c033uIvb7zOrx7/FYmJibPXbP9l37JlK6draylaW0hhYSH33H0vzi5Opm4TEc/fEY3271jj5+dPVVXloj7Iqakpk3+JhKzM7CU7cVxcXMhYlcHY2CinTp3EXaUiMSHRqnlZEIRF/ayOQK9bgaDp9aILHnQ6regWQJMlIN4UXqoBwh60Go3oD5vE5BCjUCh4/+8fEBoWgnaJ6gy7mtOcIG1sakShUCz6KplrZyMjrSd5Pf7Y43z37fe88McXeOCBB3j77XfYumUrTU2NbNp0Lv/4x8cYjUYaG5rY89FeQkNCiY+L57LLLuPG7TfQ0tJCQ0MDhYUFfPDB+5w4cZyAgABSU1MZGOhHq9WwalUGp2trRd2YXffvoqmxmfz8Aq686gqT5lwBfaNBpKCZ6T7FrjEYDZYIp16vp66ulqrqSqKiosnOznGoRc7T04v8/AJ8fHwoLTtCS0uzxZQfU4/hofIQbc5Nz0yLfpH1Op14s1arE8WEYFqjtRC1iYHR6HjO2wyNZgUCrZ0z1Q0Gg93uq2V8TpO2PXnyJB4eHouE89DhQ6SkLOaSkclkJCcn4+TsjCCYBrj89Oaf8tCDD1NVVcldP7+Lx3/1K6RSKVlZWbz11tuoPFTUnq5l/76P6ezsJDsnC41GQ0REBHW1Z2hqbOaD9z8kMTERmUxGxqpVnDhxXNSNAZM2Li8vY/369ZbCZTGQy8RrQXNXv1iEhYbR0tJMZ2cHxSXFKJRKigrX4icikQ9zBetFhUXodDqKSw7TP9BPc3MT4SJ5egRBWDbyaQu6FZjCJs0pVjg1ogoJYOWm8MwKhHNGM2MhLZhYpjXSnnBOa7Va0xe7ttY0y2KBlf3VV19x3nnnLb3zWYf5pZdeoqbmBH/721955plneOONN7lg8wVz20ml/PWdv/HlV19y8SVbMRHwuvLjn/yET/75L0JC5qgq9Ho9MrmcVasyOHFicX+pPQwPD5O+Ko1//OMf6PV6ikuK+fvf/y5K2FxcXZieEZcYd3Z2RjAKovN8KpUHTc1NDI8MU1hQuGQrl6OQyeQkJCSSk72GlpZmBgcHcVU4zu0DZsI2T9HnYRIaccKp1elEr9GI5B0CrFJjoo6l0eDiKk44NTMaXGbjNuPj44sKWuZjyTMSBEHw9PRkfHycLVu2UFtXy8ZzNuHuNvfF9PPzo7WlZcmdSyQShoaGeOjhBwkKCuKee+7hyiuvsmk+eHl5UV11jN7eXs6cqePwoeJFcyTUajU1NTWEhoTi4eFhFeiwR1Vohre3N/v3fczQ8BDjajUCAsXFJfzo/B/x8f6PHWrw9vXxpaWlGRDHVeTnZ2Ikd4QTZ3pmmrq6WnRaHQkJCYyNqVf08iwFZ2dnDHoDSYlJHD92DG8fH+Lj4h3y7wYGBvD3C1h2u/kwa1s3N5HaVrsSzSleODUrCAbBrGYXaQ1oNDMoZ2erjo+P42Yv92pvR+bav8svvwIXFxdee906V3jZpZfx8T8+XtI0dHFxMaVW/rCbUydPc/XV19i1652dnYmIiMDF1cXmgJfyinJGR0c5ePAglVWVpKfPMcIrlUqmlulgkEgkFBYWsuWiLWzbdh3nnH0Od955J1GRkdx2+60OmTcr1YL+/gEMLNNqZTAYaGio52hFBSEhoeTm5hEbE4eLszMNjfYJiB2FIAgcP3Gc0LAwIiIiKSwswt3dnZIjJbS3ty97DwYHB/D3F2dWm1vzxGtb8YI2v6XPUWg0GpxFmqdmiL2mmXnadmJyApUdnmW7wjm/av7cTedSUlJMW9tcF0NsbCxBgUGUHCmxuV6hUOLs7MyOHTtF3WSp1HbXw6aNmzh08DCff/EZl112KavS55q/lW5uTE6K69bw9vZhZGSEF198icrKSg4eOujQOl8/PwZF+pBeXl6MjI4sPQ6ip5viksPIZDKKitYSGDDXIZGSksrY6ChNTU0r9o/AJPwnT9bg6uJCZIQpiCeRSIgIj6CwoJDJyQmKS4qXHHGn1WoRjIJon25goB9/f3HaFmByalI0L5RWI97n1K4gUmswGFY0FXxmZsbCvmGvRBOW1ZxzpmNAQCAFBQV8+umnVttcetll7Nu71+Z6V1dXZjTi/DOw38SakZHBvz75lM8++5w77pgjCnNTKkVzzri6uqLXmSKCN964nQ8/+MChdf5+/qIbjs0F2Z1d1jM01OoxSstKGRwYIC83n5iY2EUPXSKRkJ2dw8zMNEcrj4rW2mBKwZSWHcHN3Z3kWXa++XByciI5OYXVq1fT3NxEVVXlolmZTU2NRCyIzDuC/mXm19iCTqdbEXfxSnxOk1krPvAkVqBNx5qxrFOr1XZ9TrtXHhkZYRlGo1AoSEpK4osvrQl3L73kUv7xf/+wud5cQSH2a69ULN/Eun7deivT18vbm5GRYVHHAZMfOjIywtYtW9n/8X6HztXb29tSaC4GsbFxtLQ0YzAY0Gg0nKg5wanTp0hJTmHVqgy7D1sqlZKamkZYaCglR0o4U38GrVaz5PZmTE5NUnOyhsqqo6QkpxATHWPXFHN3c2dNTi7h4eEcrTzKmfozGAx6pmemGRwaJEwkb61Op8NoNIp+kQcGBywMEGIwOTnhMOO7GRo7nUtLwRSpFRdMA+vgU1t7O5FRUUtuazfKMJ9XU6FQ4O/vz8GDB63yO8tRfTg5OaPT63AW4dgrZ3lExWA+g5oYP8Bv1kT9aM9HhIWGObRWKpWiUCgWzYFZDs7OzgQFBVNVbRrzHhcXT3pauqjzDQ4OITAwiK6uTsorypHL5Hh7e+Ph6YlMKkOn09Hd3cXo6CijY6NIkBAdE0OagyztZvj7B+Dr60dbWxvFJcU4OTkTFyuOexZgcGhwRULW399HpAieIsBCCCa2Vlij0eAsMj21khznwg9/Q0M9F15w4ZLb29Wc8XFxNNSbhq24zLKJ+/n50draatnm8y8+50fnn7/kPpQKBdNT4hLwHp6ejImgmgCTll4YwXUEPr6+DA0Ncustt6IeV/Pqq684tM7fz3/ZAM9C9Pf30dvbw+jIKHm5+YQEh6woNSKVSgkPj2Bt0ToyM7Pw9PJiXK1meHgIg8HA5NQU/v4BrMnJpaCg0GYvp6PHiY6OJj1tFRMT4zS3NIuiAAGTvxkg0qQVBIGxsbElJ9otheVSE0thOd/PFkzMjeKEU6e3LuRvqG8gLi5+ye3tC+c8zWk2UXU6nVV97Reff875doTTUZ6V+fD0sM/nuRR8ffyWDGYsBWcnZwTBZKr+65NPefxXj1NTU7PsOv+AAHp7HSM9npiYoLyijK6uLtasySU+IWE2HfPD4eLiQlBgEAkJiSQlJePq6kp8XDz+/v6iS99sQRAEGpsayVydRXp6OnVnai1scsvBaDQyMjKClwgSajCNTvDyXJ40eyGGh4fxstGgbg+CIJgsu5X4qSKFc36OUxAE6hvqF80Wmg+7wpmQYBpTZlbHc8S4pnrC/v5+6hvqKSosWnIfs/1qoi5CKpXi7Ozs0AswH76+pjHiYuHr48Pw8DCxsbEUFhZSX39m2TVuSjecnJzsHk+n03Hq9CmOnzhGfFw8mZlZKFxNYxTG1Gor7tj/VjQ2NqBQKPDz88ND5UFebj6BgUGUV5TR1NRot8Kqs6uTgIBA0UGd/v4+C7WIGPT191pFuR3BSiLCYOKSUor1bedVB/X29qJQKGyyXZhh9675+voikUgsMyFamltmi89NO6w4WkFAQMAiHhS9Xs8777zNFVdezp6P9ogeIw7g4+Mrep051yk2AOXrN1de5zc7U9IRJCYmcqa+btHxBEGgrb2NkiMleHh4UFhQZDUIWCqVkpWZRVNTk2gz8f8luru7GRkdscxeAZMFFRQURFHhWgSguOQwvX29i+6BwWCgpaWFODujLpaCKRgkPrqr1xtEF9eb5vKIZ9C3RVezHGbmaU6T1rQ/3MmucEokEpNpW1/P7md3c+jwIf74/AsWDbpp4yYu2HwBWdmZfPvdt5Z1dXV13Hb7bVx4wYV8tOcjTp+utfJTHYGvjw9l5WV88eUXDq+RSCR4qDwYHxfpd3r7WMzhkNBQOru6HFrn7q7CTelGf//cCIXBoUGKSw4zPTVFUWER4WHhNs0zJycnsrOyOVFzfFHK4r8BIyMjNLc0kZWZbVPzyWQy4mLjyF2TR19vL+UVZVZMcq2trYSFhoo2raenp5HL5aLX9ff3r0jbrmRollanRS6Xiza7JycnLFq6ob6eBDsmLTgwdj4yMpL/ueNnvPrqK9x158+tyKxcXFzYvftZXnnlVW644Xr+98H/RavVkpqaiq+vL4WFRXz7zbc4OztxyaUX4+3jRVh4KPkF+ZSXly95zI6ODoaHRxgeHuaee+4WVZzuswLTVi6X4+riyvj4OJEREbSJ+JAkJCTQ0FjP5OQklZVHaWlpITMzi6Sk5GVL7pRKJelp6VRWVa4od/mfwtTUFDUnT5Cdlb3sNbi6upKRsZrEhCRqTp7g5MkaJicn6ezqRCwDP5hNWnGmKZhM2iCRQ3rB5KeKGcwEK9e2o2OjliaR+oYG4m1Uwc3HssKZlppKTHQMFeVHiY+Pt2mG/ei8H1F5tIq6ulrWri3iwDcHiIyM4pN/fYJMJuOsDWez56M9tLW2U15Wwd0//zmXXHqxaQbLPMHT6XT87pnfkbMmm/Ub1jE0NIS7uzv/+MfHDt+Alfqd4eERtLW1kpCQyJkzy/ucZjg5OYMApWWlREZGsiZnDW5Kx30Rb28f4uLiKS09Ito3/09gYKCfo5UVrFqVgULhuC/m5eVFQX4h3t7eFJccnuX2FR8h7u0T7zcajUbGxydQiWADBJienkIul4nW0mOzASsxEATBxHI42//Z0PBv0JwJCQlIZVJUKhUymQypVGLzK+/v78++vfu5+eabeWDXLpISE9mwfgMAPj4+jI6ahryGhIRw9dXXUFZazoEDB4iKjmT79ht55ZWXycvL5dtvv+FISSmtLW3krsnloQcf4umnn3bYjzRPpxLrdwYGBjIyOkpMTIyFy9UeBEGgs6uT4pJigoKDcXKS47NgKpejCAoMIjMzk+MnjtPRIX6K178DgiBwpv4Mzc3N5OXmi375wORW+Pr64erqilKhpLikmMFBx9NNk5OTCEZBdBHB0PCQJT4iBr29fQQFihgAOovRsVE8vcSZwgt5nRsaGpb1OZdtdUhJSbVqzfL28mZ0dMRmraREImHnzlvYufMWq98rZ0vr5hcIhIeH89WXX9HU1MTXX3/FwUMHue+++7j22m2WbVJSUhmYfbhZWZmcv3kzv/7Vr+2aWhKJBC8vb4aHhxcNSbUHiURCXGwcff19TE9PWyZ228LI6Ai1tafx8PCkIL9gtstDT0dHu+jEuRnu7ioK8gs4ebKGoeFh0tPSkMn+fZ0o9qDRaKg+VoW3tw+5uXk/qC2toaHeMgM0YmqK2rpaWlpbSUlOWVboGhsbRM9lBejr7SUwKGj5DRegt6+X1atXi143OTklyjqC2fTQrPk8OTlJe3u7Td7m+VhWcyYnJzM6Okpnp6km1Nvbh+ERcaVrEolklh1uctHv4+LiuPXW23j3b++xbdt1Vi+Gt7c3oyMjFBeX8Morr3LsWDW33bZ890hYaOiiGlZHEBQUxNjYKL6+vjaDIDMzMxw7Xk19/RnS01aRlppmyY/FxMTS3t6OegX5WTPkcjkZGavx8fGhuKSYzs4O0QRkYqDX62lubqK0rJTYmFgSExJ/kGD29/czOTVJ0KygKJVKsrOyiYmOofpYFbV1tUv2zqrVaqamp0UHdQRBYGh4eNEs0eWg1WoxGg0WM9NRTE5OolC4ir5Po/NM4ZIjJWRmZi1LXbrsp1kqlbJ+/XoOHvyebduuw8fHh7a2VlEnBibTdnh4WFQlhpOTk4VJLTc3l48+3MOmczdx22232mTZk0gkrF27jvj4eE6dPmU1nu3Kq65gcnKSVasyeGDXAzYH0JqmiOm4+eabrW6+KS3QTHdPN4kJSQQEBCx6OE5OTmRmZlFVXUnumjyHx7vbOoeI8AgCAwJpaW3h0OFDhM+2d/27ejo1Gg0tLc309fcTHh5OUWHRD963Wq2m7kwt+Xn5i+6Nr68vRYVrZ1nqDxMTHUPYgih2bd1pkpOSRb/0arUalcpddC61r79PtG8L0NvbQ1CQeFN4bGzUUnDw/fffs2HDhmXXOHRFG9afxfffm2glXVxc0Oq0or/oPj4+DC8oTJ+ZmeGtt97kllt2snp1BjfccP0i7err62cxbd3d3fnn//0TJycnSktLLT9lZWWUlZVx8OBBzjp7Axs3baS7u5uDhw7S22vKwancVbi4uDCuVi9K/czHG2/8hdWrV1vqdHt6eyguOYxEIqGocC2BgYF22dNTklOpqqoUTX+yEC4uLiQlJlkKPIpLiqmtq2VwcFA0TQqYuij6+vo4UXOc8ooy3N3dWbd2HTHRMT9YMM0WRVZm9pLtWhKJhMiISAoLihgfH6ekpJjhYdP7MDDQj5OT84omdnf3dK9IWFYqZL194qPCRqMRvV5vqS93VDgdeiobNmzgxZfmRi2Yc4liwskeKo9FJt9zzz/Hxx9/zE9+/BNuvnkHr7zyMkVFhezZs9dij4eGhFLfUE/w7I309/fnhRdeXPI4Wq2Wf/7z//jkX59QVFTEdddts3C1ODs7s3/fx3z2+WfccMP1hISEsjojg9WrV5OQmMgH77/PN998w+OP/4rTp0+hN+hRKJTk5eY7XKrl5+fH1NQkx48fIzMz6weZiWAydWNiYomKiqavr5e+PhNThFEw4unphbeXN0o3pYV7UK/XMzw8hCAIjE9MMDIywvi4Gie5E17e3gQHhZCetuoHn5cZBoPe1PGSkuqQVeTk5ERKSioTExOcrj2FrFXO5OQEOdlLj4Vf+tgG+vv7SEwQx0phMOiZmZ4RXU87PTONVCoT3S86MTFX8zs5OcmJE8cpyC9Ydp1DwpmammrxO8PCwvD28WF4eESUcM4vyXN1dUWv1/PKKy+zb+9+MjMzAcjJyeGVV15m3fq1vP76X7jwggvx8PBgenrKxP5tI+RdX1/PF198wdT0FFNTU3h5erF27Vpee/V1GhobKCsrx8vTi5aWFotW3nz+Zupqz1BdXc3x48eprq7mr3/7G2effTbVVcfo6+tlcGiQjIzVlo+CGERERDIxOUl9Q73oF2cpSKVSgoNDLDQnBoMBtXqMkZERenvmPno6nY6enh6QmEboxUTHoFKpVtQYvBwEQaD6WDUR4ZGiScfc3d3JXZPHqdMnGRoapLOrk9iYGFFBsK7uLoKDQ8SXB66gvxRMJXdBKwg8DQ4N4eNjqhArLilm9epMh0oGHboqs9/5/exgWh9vnxX1Tvr6+llC65988k/CwsItggkm0+fWW29j75593H77bfz6N7/GaDQSHBxCzxJ1qE8++QRffPk5o6OjyOVympoa2bHjZgIC/XnooQcpLT3CsWPHSExMpKDA9LU6evQol19xGQcPHeTKK6/klVde5dDBQ9x0003UnanDz8+PrMxsekVOxZqP5KRk1OqxFQWmHIFMJsPb24eYmFhSU9MsPwqFwvTvlDSiIqPw9PT8jwgmQG1dLe5u7oSHh69ovcGgZ3BwiA3rz0Iul3O4uJju7i6H0mCCINDW1mphdBCDFUd3Vyicvb29BM5yO3///fecddZZDq1z+KnN9ztNk6nETZkCCAkOoavbVBr30p/+xM9uv93mdoWFhZQeKeOrr77iwosuQDAKlnXzIQgC3x/8nj/8fjdPPfkUjzz8CC+++BLHjh2nob6RG2+4kYGBAd55520ioyLYt28vTz39FFu2XsTWLRfT3NREckoSjzzyMAcPfY9ep2NtUREhIaH4+/uj1WpWVNAApg9N5uos2lpbbZ77/59h7lSZnp4mMTFpxftpbm4mPDwcFxcXYqJjyM/LZ2hoiNLSI4yN2a857ujswNfXT3RniNFoZEwtvh1Nq9WsKLprbt4wBwgd9TdBjHBu2GCZNSKRSFAoXB2uCdVoNAiCKbms1+sZGBigtPQIl112+ZJrgoODOfD1Ac4662zWriuio6NjUXVSa2srOp3OJhmYr68vW7dezE9vupnrr7+ej/f/g1/c/wu+/vpryssquO2223j2uef4+ON/EBYexoMPPojBYLSYVRKJhPS0VZw6fRKtTty4dDPkcjm5uXn09vRw8mTNDw4S/TdAp9NRWVXJzPQ0maszV+y7Tk1N0dPbQ9S8vLCLiwvp6atITU3ldG0tx08cR6NZzPag1+tpaWkh3k4v5FLoH+jHz9dPfMFCn/VkA8fXzWnbiYkJampOkJ+X79Bah4UzJSWFsbExOjpMk4Z9ffwYGFie5Gp8fJycNdlcdfWVTE9PExIcinpcbeKIWcbccnJyYtf9uzh08DD//OcnvPzKy5aZKgBff/0VGzZssHujXVxc8PXzIyw8jKbGZg58fYCgoCBO157m+LFjpCQnc8vOW9lx8w7O2Xg2X371pWWtUqkkISGR6uqqFecbnZycyMrKRunmRllZ6YpGxv+3YGzMNAs0JCSEtLT0FZvLRqOR4yeOk5aaZnMfHh6e5OflE+AfQFl5KU3NTVb3v6mpkcjIyBX1q7a2tNilBlkKK47u9vRY4hbm/KajLWoO312z3/nNNwcAk2br7rFvrgmCwE9vvom83DycnZ05f/P5KBQKent78PX1dbg1KyEhgbfefIv4+Di2XXctV1x5OUnJiTz40INce+22ZdfHxsTS0mwaQ9A+28rl7u5OYWGRpeRu+/ab+OCDD7nxxhvomteVEhQYhI+3D7V14kY/zIdEIiEmOobk5GSOVlbQ39+3/KL/IgiCQGtbKzUna8jKyibEAe5de/s6deok/n5+dssdJRIJwcHBptY0o5Hi4sP09fUxNTVF/0A/ESKZ6sFEpiaTyay4lx2BTqdDo9GKLivUajUYjHNtbAcOHHDYpAURwglw+WWX88GHJoY6pVKJ0Wi02xC9+9ndtLe18+KLL/HXd/7GmjVrOPe8Teh0epKTkk1RRQchk8lYv249L7/8ChdecCF7PtpLX28/F1140bJrXVxcULq58d333zExOUlhQSER4RGLNO76deu5/fafceP2G62+1HFx8WhmZn5w3au3tw95ufm0tLZwxoH63f8G6PV6qo9VoVarKcgvFP1iL0RTs4neM9bBPk+ZTEZcXDxr1uTSM0sfGhEeuSKt3dLSsqLywP6BfgJX0I7WO2/Il9Fo5KOPPuRyO67cQoi6wi1btlJeXm6ZqWmKotoWsG+/+5Znn93Nhx9+hKurK1KplN8/83tu2v5T/vznP5GZlUVvn7hoqJ+fP4Jg5OqrryEtzbZJtBDT01NUVVWi02mRSiUkJSbZNYd23b+LY8eqreqJJRIJGRkZtLW3WxLnK4WLiwu5a0z1q6VlRxgRWQr5/wqCIDAwMMCRIyUEBgSyKn2VaOKshejq6mR4eJg0kaRmYAqoREZGoVQo6erq5NQpcbEAjUaDenxcVL21Gb29PQQFizdpe3q7LeuKS4rxUHmwatWqZVbNQZRwKpVKtly0hQ9m+V1DlkhxDA8Pc/31P+Gtt94mIsLa/Lj77rvx8vQiIiKcPpGpCnNxemNT47LbGgx6ztSf4WjlUcLDw8nPKyA4OITW1qXHRwDs37+P2JhY0tLS6O/v5+TJk9TV1SGTycnOyubkqRrRE8NsXUdCQiIpySk0NzdxpPQI/f19/xWa1Gg00tXVSXHJYbq7u8jKyiZUJB2mLQwODtLa1kpWZtaKtJ4gCNTWnWbVqgwKCgrx9PTkyJEjtLa1OnTfGpsaiI6KFv1R0Gg0TE1NiWY90Oq06HQ6S4H8e+++y7brrhO1D9F36dpt23jvvXcBkxaQyWSLghy//d1v2XLRFjZt3GRzH4888kuOHTuGxgHe1YUICAhErR5bUkAEwZR2KS4pxtnZmaLCtZYOmrjYOLq6u+wGZV548UVqTtbgrnJjVUY6l1x6MVsv3gKYuJPSUtNny/P09Pf38/WBr3nuueeoqqoSfS2enl5kZ+eQnpZOb28vxSWH6erq/I8Wuy8FvV5Pc0szh4sPoR4fJyd7DRkZq0X7WbagHldzuvYU2Vk5Ky4V7OnpRqXyQKVSIZFICAsz1QTPzMxQXHLYQqVjCxMT44yOjhIaGir6uG3tbVYDpC7aciExsdFERUcSHRPFRx99aHNdf1+fJbep0WjYu28v11x9jahji75T55x9Djd1befMmTMkJiYSEhJCd3eXheKvs7OTN998g2PVS4/nCwsLIz4+AR8fX9E8sxKJhNiYWJqam0hLTbP62+jYKKdPn0alUpGfV7CIUU0mk5GakkrNyRMW03Ih9u3dh8FgwM/PD7lczpdffckzzzwDmLTKm2+9SWdnB199/RUvvfQS6emriI+L43fP/JYvPv+S9PR0Tp06xWuvvcrGTSYal+XMQXd3d1atymBmZoaW1haamg8RHh5BRHj4f7xtTKPR0NLaQl9fH2FhYRTkF/5bWPvMmJmZ4dixarIc6MJYCgaDnsamRvJyrVMQcrmcpMQkpsIjOF17mtY2U2vawmjo6dOnSUlOEa01DQYDPT3drC1aZ/n/AwcOUHPiJC4uLrS3t3PFlZfjHxDAWXVBXSgAACAASURBVBvOslrb09tDcpKJe+nzzz8jLS19kRW5HEQ/eblczlVXXc17f3+Pxx97nKCgYErLSi3C+fHHH7N27VqrsX22sOv+XZyuPU1PTzchIeK+aEFBwTQ2NVpKAWdmZjhzpo7pmWnSUtPsTmr28fHFQ9VHQ2MDCTaaXReWddWfqScxIYHp6WluuPF6+vr6uf/++wnw9+fiiy8mY9VqpFIpZ7//d7ZsvYjzzj2Pf37yT7Zvv4knnniCe+65m/vuvY+dO29Z9uVwdXUlOSmZuNg4C0GYTCbDy8sbby9vvL29V/yCg8mqmJ6eZmR0hNGREUZnE/2REZEkrF33b68k0ul0HK2sIC01DfcV8Mmacfr0aSIjo5YsOFAqleRk5zA4OEhVdSV+vn7ExcUjl8vp7+9D7uRkRbDmKLq7TUX15o9rX18fPj4+lrrv8PBw3n33Pa655mp6unstz1ev1zMzM1e7++5777Lt2mtFH19iz1436I02/1hZWcm1267hTF09EomEiqMVJCUmoVKpGBwcJC09lW8OfEtKyuKZHPMxMzPD0coKigrXiv6qdXV1Mjo6iourK93dXSTEJ9rtGJkPQRAorygjMjJq2cTy3Xf/HF8/Pz777DOio6N4/bW/WASksamRocFBMjMzcXZ24bXXXqWpuZld9++yMBSWlZVx5513EJ8Qz6uvvCZ+MI9Oy+jIKCOjI4yMDKPRaFGpVHh7eePl5WmtWSVQVVVJdlaO6X8lpns8MjLC6OioqRtfocTL2xtvLy88Pb3+rVrS6ry1GiqOVhATHUvwCoIpZvT29dLZ0UF2do5Dz9ZoNNLR0U5rWxsx0dG0tLayJmeNaFY+QRA4XHyY3DW5lo9CWVkZd911J6WlZVbb+gf4UXu6zsJs39XdxcTEBIkJiYyOjhITG01TY/OSNJgyudTmha3IZsrKysLZ2ZnS0lIKCgqICA+nta2V9LR0/Pz8+N8HHuS+X9zLvz751O4NdXV1xdPDk/7+OfvcEQiCgFQqo629jZjoGIoK14qKJJpL60rLjuDu5mb3q56amsrd99zN6tWZ/PWdv1ldT1xsHB4qFaVlpaxalcGOHTsXrc/Ly+O7775n584dnHXWBvbs2SvKvHF2ciYgIMDShCwIAuPj44yMjtDR2YnROFt1JICAyUxtbGywjDl2cXHGy8ub0NAwFArFv60bxR5GR0c5UXOc5KTkFU0XM2N6eor6+jPk5xU4fN5SqZTIyCiCQ0I4erQCrVZj4u4RKZyDQ4N4eHhYaWt/f3/aO9oXzYENDwuno6PDIpwdHe2kpZom4O3bt5eNGzfa5add8lpEr8D0cl977VxgKCAgkNHRUUvO8/bbb6etrY1PP/vU3m4A03AfMaPtxsfHKa8oo7+/j/S0dKamp1cU4nd2dmZ1RibVx6rtMt/96EfnMz09zbnnnmvzBQkICCQ7K4eTJ2uWzIMqFAreeeevXH3NNRStLeTw4cOiz9cM89iJyIhIVqWvYnVGpulndSaZs90OmZlZZM3+pKakERoSilKp/I8LpiAINDc3cer0SbKzcn6QYGq1Wo5WHiU9bZVoNvbZk0Gv17MmZw0NDfVUH6tiesbxKHtzcxMxC3KiMTExhIWGLRoVGRoWSmenqXJubGy20MFi0r7HNgcKZWxhxU7Gtmu38dGej5ienkYikRAdHU3LbJrCycmJ3/3uGXbtun/ZxmClUolSqbSQOi8FrVbLyZM11Jw8QWJCEhkZqwkPj8BoNNA/0G937VLw8PAgLjaOY8eql/w4mKPC9io73NzcKMgvYHBwkJqaEzZraCUSCffecy+vv/4Xrrr6St544y8rOuf/Vmg0GiqOljMzM0NBfuEPivIaDAaOVlYQH5+ARqMRzfwPcKa+ntjYODw9vcjNzSMkJJSjFRU0NNQvW+M8MDCAk9zJ5qTzy6+4gr179lj9zqQ5Td1HzS1zQt3a2kpNzQk2b75A9PnDDxDO6Oho8vPzeeuttwBTznNgoN+SGL5g8wUEB4fwl7+8vuy+4uLiaGy0nbs0Go20trZwpPQI3t7eFOQXWnXMp6amUVdXi8Egnh0ATIUUKpVqyREM5pavvNw8u/uRyeSsXp2Ju7uKsvLSJb/SPzrvR3z/3UHu+8V9/7UFCGIxODg4Sw0aRUpK6g8KLBmNRqqqKzl86DCJiQmkpCZz0ZaLROWA1Wo14+PjVmWGgQGBFBWZ3J/iksP09HTb3KfRaKTuTK3NGaZgqpLb//F+q9+FhYfT0dnB9PQ0U5NTlrLE3bv/wE9/evOKg3g/KDy36/4H+MPu36PX6022fkSkhZBZIpHw26d/y69/82srJnBbcHdX4eTktKhHdGBwgOKSw2i0WooKiwi1MaLPPHvEPHBpJUhMTGJMPWaz2sl8PEdusNmCSExMoqKifMncW3x8POeeey779tkeOvz/F5he5DoaGxvIy81bESfPfAiCwMlTNXh5elF9rJo/Pv9H+vsGOH78mIVgzpF9nK49ZTN1IpVKiYmJJS83n4HBAcrKShexc7S2tRIYGLSkjzoxMbGITiU8LIyuzk5aW1uIijYVOvT19fH39//OXXfeJeIOWOMHCWdBQQEREZF88MH7AISFhdPd02MxZbOysth4zkZ+/4ffL7uvuLh4Gma15+TkJBVHy2lvbyc7K4fEhES7yevIyChGR0Z+cO9lY2MD6gWjHM4+62z0OnGtXr4+vuSuyaOpuZETNbbbnq695lr+/v77Kzrf/waYJmWXIpNKycvL/0EpHjMaGhuQSKT/X3vnHR5llfbhezLpycykk0YaSYD0HjpYV1dwBaVbVyzfunakqKuuuirsWkAsq+7qupIgdl3XgiCdFNJJSG+kJ5M2k2Qy7f3+mGSkBZKZCcXlvi6uCeHNzBny/uac85zn+T10dMj58ssv+eHHH9BqtfT394+6yLmxsQEHB8cz+hHZ2dkREx3L1KkRFJcUU1hUyODgIIODgzQ0HGNSyKQRf3bv3r3MmT3nhO/5T5xIS2srbe3txgqUza9vZumSpSYVZw9j9sHW2rVr2fjXjUafV38/f+qPC4w899zzvPHGllOaHZ2MTCZDEPTk5eeRl59LSPAkEhMSR3X0IBKJiIuPp7jkyJg2/cdjY2NDXFw8+fl5p5iMmYK9vT0pyal4uHuSlZVJcUnxCXuna665loKC/BMqYC4G1Go1pWWl5OTmMGXyFMLCwi0SaKqrr6O3tweJs4TFS24iLi6e7du38+WXXxAWGjaqIx9ln5KamhoiIyJH9ZoymYxpqdPx8PAgMyuD7MPZTAoJPWOAcd/+vcyePfuE7/l4+xAUFMhEf3+srKzo6enh3Xff4dFHV49qHCNhtjivvupqbG1s+c9/vgEMvVWOHfvFbzUgIIDp06fz9Qit6cGwFDl27Bh9ff10dXYyfdqMMScoO9gb7Dny8nJNLmqWSCTExsSRk3vYIgIViUT4+voya9ZsXF1dOZyTTVFRIf39/djb23PD724wVvlc6CiVCgqLCsnMzMDRwYGZM2YafXHMpaW1habGRuLjEtj8+mYWLVzEY6tX8+knn9LT04O1tTVq9ZmT3PV6PQX5+UTHxIwpRVAkEuHrY6hPVasHqaquPKEx1fEIgsC+ffuYM+fE4KCnpyexsbFMHCpje+vtt7jmN9cQHDz2XjHHY7Y4RSIRa9eu5aUNGxAEAWtra7y9vU/wzlmxfAVp6Wmn/fnOzk4OHjqAQtHL7Fmz8fPzo66+zqSxuLu54+PtS3HJkbNfPAIymcwo0LPtlUfCYFj8S37s8A0wc8YsvLy8yMvPJb8gnyVLl7ItPd3ksY43w5UpWVmZFBcX4z1hArNmzSYgwLSSrdPR0dFBRUUFiUlJiMVi0tPTePOtN7nl1ltwdXNj1aq7sHdwOOOHO0Bp6VF8fHxMaiMhCAJlpaUkJSaTnJRCQ2MDWdmZp6z2hlNDT/ZM6h/oJz8/H51Ox8DAAK+/vpk1a9aOeRwnY5H/4YULF9Hd3cXuIQOw4KBgamtrjNGwBQuu59ChQycUVw8MDJCbl0tlVQWxMXFERERiY2NDaGgYDUORL1MICgpCpzMUVZuKTCYjPi6B/Py8MdWcDnPb7bcSFh7K8395/oRlq0gkYsIEb2ZMn4mfrx9OTk5c+9trKSgYOQ/5fKDT6aivr2f/gX00NTcxZepUUlOn4eU1ugys0dLU3ERpWSnJScnY2tgyODhIe3s7Tz/9DP/8xz956MEH+frrrxCLxWdcara0tKBUKk2q1Rweh7NEglQqxcHBgYT4BEInhVFQmE9xSbHxHPzw4WxSU0+M2htc86vZt28f7e3tvP/++6SkpBAVFXW6lxoTFsmqFovFPLb6MTZseInL5l2Gra0t7u4eNLc04+vji5OTE3Z2dqhUKnQ6LVXV1bS2tDB58uRT2r2JxWLDRr34CElJY/cyFYlExETHcCjjEBKJ1KTMDDAscadNmz5UaNxD+BhaFeTn5zNpUig7d+5k06bXmDFjBqvuXMW11/7W2NfR09MTqVTKG2+8QVlZKQplLxMmeOPu5oZM5mJ27eRYEQSB3t5empubaG1rxcfHl5Tk1DEbaI32tWpra2lrayU1JdW4n7Szs0Op6MPe3n4o6nqUrWlbsbKyInwEi9HOTjkVlRWkpprW40Wr1VJVVcm01BN9ZN3c3JgxfSYNDcc4eOgAQUHBuLm7n3L8VVFZQVBgIBJnCY2Njbz8yt/Y+tHpV4ljxaTc2tOhVqsJnxzGZ59+TmJiIgMDA+TkHmbmjFk0NDSQOi2F7OzDVFVV4u8/kaDAoDMujfLycvHx9TXJVAkM0cTDOdljMoQ+HcPdtxS9vcTFxY8qMHHgwAG2b/+Yb/7zDYODg0RHR6NU9lFdXcXSJUtZsWIlLi4urFy5goDAAN595z2cnJxob2+js7OT7p5uxGJr3FxdcXNzx9XVddQ5sPv272X2rDlnvU6lUtHd001Pdzc9PT0MqAaQOEvw9vZhwoQJ4/bhoNFoKDpSiNhKTHR0jFnL497eXvIL8khOThmzK94wZWWl2NnZnbGXqEajobKqkoaGBl555WW+++/3gCEAlZ+fx8wZs5i/4DqmTo0gLy+PnT/tHNMYRsqttZg4AbZs2cJ/vv2G7/77PSKRiKKiQtzdDZGw9vZ2pk2bzuTwyaMSy+DgIJlZGWaVMLW3tw+VGqWavUdqam6iqqqS+KFEg9EgCAK79+zm5ptXUlVZzbFjx0jflk5a2laamprY8NJG7r333tN+4qvVajq7Ouns7KSrqxNBEAzVKa6uODk6YWdvh5XICpGVyPAoEiESidh/YB9zZv8SsNDr9QYXgN5eenq66e4x1MLa29shk7ngIpMhk7lgbz/25jxjpbNTzpHiI4ROCjtr1dLZGP7wTYhPHLNz+zB9fX3k5uUyc8bMUd0f3d3dvPLqyyxYcD2REZEcKS4idFIobm7u3HXXKr7+5mvS0tK54vIrxjSOcyJOjUZDUnIizzz9DAsXLjL0xDh0gI6ODhobm1jz2JoxDbq5uYmGxgaSEpNNvnEqqyoZHFQRGWH+HqBX0Ut+fh7h4ZPHNKMvuH4+N/zuBu68cxVarZbKykr8/PxOmx42Elqtlu5uQ3VKf38/g4MqBEFArxcQBL3xa4Wi11i1P9wEys7ODqlEikzmgkwmO2cJ8MPo9XoqKivolMuJjY0bc2XOyQwODpKVnUl0VIxJ/VWGx3Qo4yCREVFjeo6w8FA+2f4pnV1yxFZiZs6chbW1NQ899CDFJcXs+PGnMY/FolUpI2FjY8OmTZv5/e/vIDQsjK6uTnx9fcnOziYuLv7sT3ASPj6+dHV3U1l59kajIzEpZBK5eTkGd3ATe2cOI5VImZY63Wh4FRYaNqqb/JFHHuX++/9If/8Amza/Rk9PD+1tZ7cVPR5ra2s8PDyMlQ8jMdpl7bmiv7+f/II8PDw8SU2dZvYKRqvVcjjHUKJoqjDB4Fbv4+0z5ucIDTW4aUgkzvh4+xgdN9K3pZOZkWXyeE6HxX36582dR2pKKm+//RYzZ8wiMiIKOzs7pk6datLzTZk8BblcTruJye0ikYi42HhaWlvNiuAOY2trS0pyKjqtlqzszFH50M6bOw9vbx/27d/Ls39+1qyb6mKisalxSEhTCQ8LN1uYer2enNzDBAUGm1Xx0traSl9fn0nR3dBJoXR0dODr60dYWDjTp03nmaef5v4/3k+QCX64Z2JcPDBefvkV4hPiWHzTYl7f8jqDqkFmzJhBSHDImJdTVlZWxMcnkJmVgZOTs0lLIrFYTFJiEodzshGJRMbDYlMRiURMnRpBZ6ecnNzD+Pn6ERQUPOLNJxKJ+GmHYblz4MAB3MfY9OdiQ6vVcqT4CHq93mK2J4IgkF+Qj6enl0leQMMMDAxQVl562j6io2Hy5CnY2FgTPBRA2rlrJ3X1dXy1+muTxzQS49LhxtfXl7Vr1nLTYkMa1qeffoabm7vJyQV2dnZER8WYlf0zLNDGpkZj7Z25uLm5M2P6TDRa7ZA515mb8KjVatauW8uiRYuM39NoNBQWFvLhh/9i9WOr2fHTDouM7XwgCAKNjQ0cPHQAD3d3EuITLCbM4pJiHB0cTqmxHAt6vZ78gjwiI6LG3MZvmOCQIH7+eTdisZjBwUEefvghXn31NYvkFp/M+LSfAu6//wF8fHyIjorC3t6esNAw6uvrTpsEPhpcXV3x9/fnyJEiky0kxWJrkpOSaWxspKqq0iJWlGKxmMnhk0lOSqGnp5v9+/fRMIKD3iOPPIy7uxvhYWHcffddpKam4ObuyvIVy/jhxx+QyWT8/vd3sGHjhgvCJnO0CIKAXC7n0KGDdHV3My11Ov7+pnUeO91zlxwtRq/XmdU0CaCsvAx3dw+TvGsB5J1yJBIJ3377HwBefuVlIiIiuPaaa80a10hYNFp7Mj/v/pk77/w9R4qKcXR0pLW1haamJuLjE0x6PkEQKCgswNXV1aTWb8Po9XqKS4rRajXERMda9ExvcHCQquoqOjraCQkOwdfXDysrKz744H1W3bUKmUxGREQky5YuIyEhgZiYmBMKkxsaGrhp8U2EhASz9aO0MS+9zmVAaNidrq6uDgdHR8LDwswy8jrd8+cX5OHs5DymJJDTUV1TTXd3t8nNlwRB4MDB/cTFxhMcEsT33/3Ab6+7lsyMLLP3mufkKOV0rLx5BSEhk3ju2ecAyMnNwcfbx+RzLp1Ox6GMQ0RGRuLqYlr2zzB19XU0NBwjISHR5EPskRgcHKS6ppr29jaCg0JYsnQJCfHx3H33PcTGxp7xZ1UqFe4ebrQ0t47puAXOjThVKhV19XUGJ/QJ3gQGBll8WadWq8nJPYyvj6/ZUfaGhmM0NTeRlJhsclCqoqIcvSAwOXwyl11+GTqdjquvvponn3jSrLHBeRRnY2MjiUkJfPff74mPj0ej0ZCRcYi4uPgx33jDWCr7B345GI+OijbJPvFsqNWDVFdX09beRlBgMH5+vqPyorV3sKO3RzFm/5zxFGdXVxc1tTX09/cTFBiIj4/vuGQSGZIDcoyOiubQ2tpKVXUVqSkpJnsAy+VyysvLjEdBv7vhd2RnZ1FdVWORD6XzJk6AbdvS+fOzfyYrMxuJREKvopeCgnymT5thsgN4W1sbVdVVpCSnmH2DGPqp5BIQEGB2JHck1Go1tbU1NLe04CKT4e/vj5ub+2mXWDqdDnsHO9SDmvO+rNXr9TQ3N1FbV4uDvQNBQcG4urqOSxKDIAgcazhGXW0tMTExyEyoMDkeeaeco0dLSElJxdbGBJMwfslUS05KwcHBgZqaGmbMnM7XX31DcvLYc79Px3kVJ8A999xNf38/H374b0QiEQ2NDbS2tpIQn2DyL7qurpbWtjaSEpPMPkPT6bQUFBZgZ2vH1KkR49aqXRAEOrs6aWhooLu7C0dHJ9zd3XF3c0cqlVJbW8tdd9+FIAhjztEEy4hTo9HQ3t5GS2srCkUvE7wmEBgYZLTuEASBp595GgcHBwIDAggICCQ6OhqZTGbya6rVagqLCrG1tSFiaqTJH9rD9Pb2kF+QT0pyqsmzmyAIZGdnERgYxIQJE1Cr1cydN5elS5by0EMPmTW+4znv4uzv72f6jGk89OBD3HHH7wEoOlKEk5OTWeHxmpoa5PIOEhISzRaUIAhUVVXSIe8gIT7RNEvGMb5eX38fcrkcubyD2tpaiouL8fX1Y+HChbi6jH2GGqs41Ro1vT299PT20NvTg1KpxMrKCg9PT7wneCOVSo1jMKQI6uno6GDK1Mncc8+91NfXUVdXz9GjJfzmN7/h5ptv4eqrrh7TEUp7RztHj5YQHjbZLFuPYfr6+sjJPUxiQpJZLoCVVZWo1Woihsy+HlvzGOXlZXz5xVcWXTmcd3GCoVj18isuM7rB63Q6MjMzmDJlyhkbqZ6N6uoquoYicZaY8VpbWykrLyU6KsbkkrOxUFNTw11330V/fz9vv/U2Li4uyDs76O1V4OjgYEh2d3LC0dEJJyenM84qI4lTEARUg6oThNjX34+NjQ0yqRSpTIZMKsPZ2fm0/4c9PT1cffVVdMg7uPaaa8nKyiIrK9v473K5nE8+2c6WN7bg7e1jTLo4EzqdjtKyUpRKJbExsRbZv6lUKrKyM4mNiTNrJu/s7KS07CjTUqdjZWXFt//9lj/+8T4OZ+eYfBQzEheEOAE++OB9Xn3tVQ4dzMDR0ZGBgQGyD2eZtfwAqKysQKFQEGdiqPxk+vr6KCoqRCKVntVgzBz+8Y/3eOLJJ1i9+jEefujhE/bPgiDQP9BPd1c3ff199Pf30d/fj047QiKGyOAOJ5VIsbIyVKxoNBrjmamdnR0yqcwoxNEaTatUKq6bfx0RU6dyyy23svn1TUyY4M3Lf3v5lGsrKyv57XXXUl52ZjfEXkUvhYUFxuwqS/zONBoNmVkZTJ0y1awsLLVaTUbmIZISk3F0dDSWPH68bTuzZs0ye5wnc8GIUxAEbr/9Nuzt7fn7398BoKOjnYpK80u7ysvL6B8YIMbMOsHjx1p/rJ66oe5VHh6eZ/+hM6DRGAI8YrEYkUjEJ59sZ936dfznm29Nzj0+mb379jBzxiz0er3RNsac/wudTsfSZUuwtbXl3x9+dNbgW3V1NVf/5ioqK6pO+++GQusaGpsaiYmJRSoZuenUWNBqtWRnZxEcEmJyDTAYAmCHc7IJmBiAt7cPWq2WK6+6kmuuuYZ1a9dZZKwnc06qUkaDSCRiy5Y3SJ2WQnp6GsuXr8DDw5Ou7m5Ky0qN63tTCAsLp6qqkuzDWcTHJ5gcoTt+rIEBgXh5enHkSBENDQ2EhYcbG6KOha6uLgICJxrtV6ysrHB2dmbnT7ssJszhMZ/N1mMspKVtJSsri4ryylE9p1gsHjHFUqFQUFJSjEQiYfq0GRYbo0qlIifnMMHBwWYJUxAECosKcHNzx3vI4vLZ557F3t5uzOWOlmDc0vfOhEQiIW1rOg8/8jBlZQan9dBJofT1KU3y7BlGJBIRGhpGYEAQmRkZZ7XjHC0ODg4kJSXj6+tHfn4eBYUFo6pGOR5XV1cKC4q4+667cXFx4euvv6G9zdCh7ELmxhtvwsfbh1defWVU1yuVylPEOTg4SGFRIUVHigwdvSMiLSbM4d45U6ZMGXMryeMZThO0tbUjdJKhxd+On3bwwQfv868PPhy36P2ZOOcz5zBxcXG8+MKLLLh+Pnv37MPb25u42HgOZRzE2dnZ5AQFAG9vbxwdHcnNy2HKlKl4mVFeNIxIJMLLywtPT0/a2tvIy89FIpEQFhqGg8PoKmWCg4O5997/4/MvPidiasQ59wkaKwcPHuTlV/6GVCrlrbfexMHeAYVSga+PD3feueqU63f9vIvbbruVB4ZcznU6g/lVc0szYaHhREdFWzTK2SHvoKSkmPi4BLPuFzDELLRaHTHRBs/bgoICbr31FrZt+9jsRAhTOS8z5zB33PF7brvtdq677rf09PRgY2NDfFwCefm59PWb5xsrlUpJTZlGVVUl1dWj72J2NkQiERO8JjBj+kwmeHmTm5tLUVHhqNwC1Wo1d9xxOxte2jDmLsfnAkEQ2LtvryGyq1Jx1dVXctWVV5E6bRogYtvH2+jq6mLd+nW0tLQAhvfU3NzMc88/x2233cr773/AfX+4j6OlR9l/4ABia2tmzZyNj4+PRYXZ0NhAWWkpKcmpZguztraW3t5eYqJjEIlE1NTUcP3vFvD65teZe5JH7bnknAeETkYQBB588AGKS0r49j/fYm9vT09PDwWF+UNNT82ztNDr9RwpLkLQC0RFRVt8thIEgdbWViqrKpDJXAidFDpin42nnn6KwsICvvj8y3GzCTEnCeGvf/sr69evY/vH21m4cBEOjvbGFMK0tK088ugjPPbYGoMXUnoaarWagYEB3N3dSUxMZOOGv6JQ9qLVaA2Nib29Lb4cHLY86enpISE+wewoemNjg9EKRywW09bWxtx5c3jg/gf5v//7PwuN+sxcMNHa076OTsfKm1eg1+tJT9uGWCymq6uLoiOFZh+xwFCEsK6W5uZmEhMSx83usaW1haqqSlxkLgQFBZ1QoZGVlcXCRTeQczjXIgftI2GqONPT03jiySd4+qmnefGlF3nllVe5/voFtDS3Gq1RKioqeOjhB6mrq2PduvXMv24+jo6OhhrZY8eQyWQEBQUhlZp+vngm+vr6KCjMN7aVN1f4rW2tVFVWkpKSirW1NQqFgiuvvIJrrr2WPz/zZwuN+uxc0OIEQ9Bg/oL5hIeFsWXLG4hEIjo75RSXFFvMP3U4E8XcA+ozYZhJW6g/Vo9Go8Hfzx93dw9mzJzOE48/wdKly8bldYcxRZw7d+3klltuZsePPxEZGcnSZUtoamrm3nvuYeXKm0+4VhAEvvnmazZs3MCiRYuIjY0jJDgY/4kTzY6Oj4ShXUf9UPd0TomZ2AAAFHtJREFUyySGyDvllJSUkJqSiq2tLWq1mgXXzyc4KJi33nr7nBqgXfDiBIMP6eWXX8bvbriBPz35J8Bg119adpSU5FSLpNMp+5Tk5eUyKSTUbHvGs6FSqWhsbOBwTg4tLc3cuOhGPDw8x6VqfhhTxDlz1kzuvecebrnl1jNeNzg4SFtbG42NDQgIZGRksGHDBlavfoxHH3l0XG5oQ6S3ADs7eyKmRlgkGWR42zS8KtPr9dxy682oVCo+3rZ93BJORuKiECcYrPXnzJ3N6kdXc/fd9wCGCpSKynJSklMtYnthMDYuAgSjAdl4MZyyuHfPPqysRHTIO9CoNbi4uODu4YGHu7vJlhmnwxRxPvf8c7S1tvL661tO+L4gCCgUClrbWmlrax2KWE/Ax9vHmLOalraVx594nOqqGovvL1tbWyktK2XK5CkWi5gq+5Tk5uYY824FQeCRRx4mLz+f7/773YjxgvHkohEnGFLALrt8Hq+8/AqLFy8BDKKtqq4kKTHZYmJqaWmhrLyMsNBQfHx8Lf7Jr9frmTtvLsuXLecPf/iD8fs6nY7u7m465B3I5XJ0Oi1urm64e3jg7uZu1geQKeI8duwYKanJFBQUoh4yoO7t7UWtUePk5MQErwl4eXmd8iGi0WhITk7iqaeeYtGiG00e88lotVqOHi1BNagiJjrWYr9vpVJJbl4OcbHxSKWGzKS/vPAXPv3kE37+efd5c0W8qMQJhn4jC66fz5+e/JNxBu3oaKfkaAlxsXEWCzqoNWqOlpSg1miIjoq26JLz008/4eVXXuHA/gNnnFV0Oi2dXV3IOzqQd8oRBEMzJScnQ6K7k6MTjo6Oo5qZRiNOtVo9JMAeQ4t2pYKKigpcXFyIjY1DKpUilUjPKorXXnuN73/4zujwr9freeHFF1iyeAnh4ab5DHd1d3GkqIjAwEAmTgyw2AemXC6nuOSIUZh6vZ6169byww/f8/13P4z7FudMXHTihF+SqG+55VaefOJJRCKRcc8YFhpmTLGyBG3tbZSWHiUoKJiJ/hMtclMsXHQDixYuOute7mQ0Gg0KRS99ff309SmHkt4HYPh3JTKk/4mtxFiJxYjFVlhZGR5bWlqYODEAsZUVgiCg1mjQqNWoNWpjEryNtQ1SmRSpRIZUKsXZ2ZmvvvqSLW+8wa6du0Y1xpaWFmLjYti7Zx+TJ09GEAQeffQRPv/ic7w8vdi//8CYYgQ6nc7ovRQbE2dWqdfJNDY2UFNbS1JiEvb29mg0GlatupPqmhq++vIri/UZNZWLUpxguAnmz7+OadOmsWnTZsRiMRqNhty8HNzd3Jk0KdRin64ajYaKinLknXImTQrFx9v0g/POzk5CwyZRW1NnXEJZiuG6Sp1Oh06vQ6/TGx/z8nOJioxCp9MjEoGNjS22trbY2tpgY2M74vtRq9UEBQfyw/c/Eh0dfdYx7Ny1k2eeeYZ9e/cB8MKLL/DJ9u3s2vUzd/z+dvz9/FmzZi0TJ575g04QBJqaGqmqrsLP14/g4BCL7V0FQaCysoKu7m7jmahSqWTJ0sXY2NiQnrbN7NYQluCiFScYoriLblyEm5srH/7r38YIW8nREjQatcUd9FQqFZVVlXR3dxE6KYwJE8bel/K9997lp50/sS393HauNicJ4d133+HNN9/k4MFDZw2MtLS0EBMbTWtLGzqdDl8/H1595VVWrryZ9vZ2/nj/fezbt4+2tjaOFBUzZcqptpYdHR2UlZXi4uJCWFiYRQNjer2eoqJCrKysiIyMwsrKio6ODhZcv4DIiAjefvvv5zwqOxIjifO8pu+NFqlUyrf/+RaRSMR186+jp6cHKysroiKjcHN1IzMrA5VKZbHXs7e3JyoyisSEJNra2zh48ABtba1jSgFM37aN5cuWW2xM54JVq+4iMiqShx568KzXTpgwASsrK5qbm7G2tubzz75g9WOryc7OxtPTk4+3bedYfQM2NjY0NjWe8LPDyep19bXExcUTGWm6yfPp0Gg0ZGdn4SyREBUVjZWVFXV1dcydN4fLL7+cd99974IR5pm4KMQJhkLhtK3pREZEcPnllxmrVwIDgwgPm0xWdiY9PT0WfU0HBwdiomOIj0+gubmZQxkHae9oP6tIGxoaKCoq5JpxMhseL0QiEW+9+Tb7D+xn69aPznptVFT00JEUzJo1iz/84T4++ujfxms+//wzEhOTuPyyyw3G051ycnNzOHKkiLDQcLNtRE6HQqEgIzODgIAAJoVMGmpFWcTceXO45+57+cvzfzmnCQbmcNGIEwy1gps2bWbRjTcyZ+5sY7mZh4cHiQlJFBYV0NxiesnZSDg6OhIbG0dMTCwNx46RkXkIuVw+4vXbP9nODb+7YVzPT8cLiURCeto2Hl39KEePHj3jtdFRURw5csT49065nIBAg9m3Xq/nhRdeYP369dTV17H/wD6O1dcTHBzCtGnTLW7/IggC1TXV5BfkExsTg4+PIfq6Z+8efnPN1Wx4aQMPPPCARV9zvLmoxAmGT+wnHn+Cx9c/ztx5c0hPN7T4dnJyYlrqdOrr66morBiXdgbOTs7ExycQFRlNbV0NGZmHaG1tOaX1wrb0dJYtv7iWtMcTExPDC395gWXLl9LXN3J1kIurKw3Hfuk74+jkxNNPP8XUiCksW76MBdcvQOLsjFarJTUllbi4+HGx1RwY6DdubWZMn4FUKkOv1/PShpdYsWI5H37473FPmxwPLoqA0Ejk5+ezYuVyZs2cxWuvbcLR0RG9Xk9p6VEUSgUx0TFmV7WcCYVCQf2xOjo6OvD09GKi/0ScnZ2xtbNhoF91XvY1lvKtFQSBP/7xPn7c8SPPPfscS5YsPSGKOjAwQPjkML7+6htjwbher6exsZGy8jL6lEp8fHxJTja98fFoxtjY2Eh1TRWREVFG463W1lZuv+M2+vsH2PrRVvz9/cfl9S3FRR2tPRMKhYL77vsDBYUFhj1ppKFYVi7voLikmMDAIAIseJh9OvR6Pa2trRxrqEer0bLxrxt5/58fnBPnvpOxtKn07j27Wbd2LYhE7Pxpp3GPmJOTw+133EZR4ZET2jN4eXoREBhokpXLWFCr1RQdKcRabE1ERKQxq2rXz7u4/fbbuO2223n6qacvisDPBeMhZGkkEgn/+teH/OtfH3DFlZfz4gsvcvvtdxgqQabPpKys1NiifLzOtKysrPDx8cHHxweVSoWbmxt5+bk4Ojga8mc9PHCRuVw0gYjjmTd3Ht9++18mhYYYBSAIAn5+fkRHR3Pw4AEQQcDEAGbNnH1O3B3a2to4WnqUyeHhxkQUrVbLc88/xz//+Q/++c/3uerKq8Z9HOPNRS9OMOxDb7/9DlJSUlm+Yhm7du3izTffQiKREBkZhbxTTk7OYQKG3MnHUyT29vYcOnSItWvWERgYSEdHB3V1dRT2FOLo6ICHuyceHh44OztfNGLdu3cPs2fPpqu7i7a2Nrq6OnFydEIkEiGVyoiKijon4zg+53Za6i99chobG7n5lpuxsbEhO+vwuNbLnkt+FeIcJiIigkMHM3j44YdISU0mbWs68fHxuLu5M2PGDErLysjMyiAmOtais2htbS0fb/+YyooKkpKS0Ov19Pb2Ym9vj7+/P/7+/r+4u3d0UF5ehrLP4C/rMTSzjufeeKyo1Wp6Fb0oenvpVSgYUA2wePFiOjo68PH2ITIiEisrK1QqFaWlR8+JOLu6ujhyxJBzGzXxFy+i/373X+66axX33fdH1q5Ze8H7Mo2FX5U4wXDs8fe/v8O2belc+9trWLXqLtavW4+TkxOREZF0dnaSk3sYTw9PQkImmVUjaig63khlZQWLFi4iPj6BzMxMHB0cT/G1EYlEODs54+zkTGBgEIIg0NvbS0dHB4VFhQwODuI05Oju5OyMs5MTTk7O49oSQq/Xo1QqDUJUKOjt7UWlUmFra4NUIkUikRIYEIizkzPLli0lKDiId995zxgYio6OobCoiJtuWjxuY1SpVJRXlNPX10dCYqJxL9vW1sbjj69n566dpKdvY87sc9OT9Fxy0QeEzkRTUxNr1j5mcJH728vccMNCY/VEY2MjNbXVeHv7EBwUbFKZ1n33/QG1Ws2bb75ldp2pXq9nYGAAZZ+SPqWSvr4+lH19aDRqxGIxdrZ2Rj/aM/05WnqU2JhYrKzEaHVaNEOJ7xqNxpAEr9Gg0aiNhmTOzs4GIR5XiXK65bZKpWLNmsf4/ofvSU/bRmJiIp999ilb07by+WdfmPXeT4dWq6WqqpK2tjZCw8LwnuCNSCRCp9Px97//neeef5aVK2/mqT89ZfHc5XPNrzZaOxp279nNAw/cj5+fH5te22wsZ9Lr9UOO7nX4+/kRFBQ0ph6Ou37exfp168jMzBqvoQOGG1WtVhsT3XVaneHr4T/6X76ura3B19cXvU6P2NoaWxsbbGxssLG1/eVrG1vs7e1NSjD/9NNPuP+B+1m+fAXff/8dcXFxpG1Nt9h71ev11NXXUV9fT9BQ2djwOA8ePMgDD9yPVCZj86bN52yvO978T4sTDPmWW7Zs4aUNL7Jq1V08vv5x47GATqejrq6WYw0NhqDRxIBR7V20Wi2BQQHs2b2X0NDQ8X4Lo+JcdLaurKzk9dc3s3DRIubOmWuRwJYgCDQ3N1NZVYG3tw8hwSHGY5C2tjbWr1/Hjp92sHHDRpYuXXbRBNNGw0Wd+G4JbGxsePjhh8nLzae+ro7omCg+//wzBEFALBYTEjKJmTNmotVoOHBwP/X19adk/pyMtbU1ixYu4tPPPj1H7+LCIDQ0lE2bNjNv7jyzRaLX62lpaeHgoQN0dspJTZlGeFg41tbWaLVa3njjDWJio3Fzd+dIUTHLli3/VQnzTPzPzJwns3vPbh588AF8fHx4/rm/kJSUZPw3tUZtaBXf1sqkkEn4+vqNeEP8vPtn1q9bR0ZG5rka+hk5FzOnJVCpVNQfq6e5uRkPD3eCAoONKxlBEPhxx488vn49Lq6ubN602Zhc8mvkV5uEYCrz5s7jcHYO77zzDjctvpHJkyezdu06Lpt3GbY2tkyZPIXgoGCDY3xNNRP9J+Lr53eK/WOnXI67h2X7Nf5aEQQBuVxOXV0tqkHVUOLCTOM+X6fT8fnnn7Fx40bUGjVPPvEkN920+H9mpjyZ/9mZ83jUajVpaVvZ+NeNyGQurF+3jvnzFxgDEYODgzQ2NtDY1IiTkzMBEwNwd3dHJBJx9913ER0dw/3333+e34WBC3HmVGvUNDQ00NjQgFQmIzAwEBfZL2ZaarWajz76N3/9219xc3Nn3dq1XHfd/PPSPOh88D8fEBoNOp2OL7/8gg0bNqAaVLHmsTUsXbrshLS17u5u6o/V09PTjbe3DzfddCOffvKZyYZWluZCEmd3dzd1dYY+JIZkjIknHDkplUree+89Xtv0KlOnTmXt2nUWCzBdTFwS5xgQBIEdP+1gw4YN1NbW8Ogjq7njjjtOsO4YbtaamZXJjOkz8PP3x9PDw6IV/aZwvsWpUqlobWul4dgx7B0cCAwMxN3N/QTBdXZ2suWNLbz11pvMmTOHNY+tJTEx8byN+Xxzac85BkQiEVdfdTVXX3W1wdV840s8+9yfuenGm1ixYiXTp0/H2tqaAwcPUldby6o7V9HU3ExdXR16vQ53N3c8PDxxc3P7VaWTnQ5BEOjq6qKtvY2OjnasxdZ4enmROOR0N4xWq2XHTztI27qV777/jht+dwO7f97D5MmTz+PoL2wuzZyjpKamhvRt6aSlbWVwUI2Dgz3d3d288cabLJi/wHidVqtFLpfTIW+ns7MTG2ubofxZT2Qy2bgv2c7FzDk4OEh7extt7W0oFApcXFzx8vLCw93jhGWrIAhkZ2eTlraV7Z9sJygomJUrVrJ48WK8vMzvmfpr4dKy1kIIgkBOTg7P/PlpMjIyCAmZxMoVK1m6dOlpqyFUKhUdHR10yNvp6enF2cnJIFZPz3GpebSkOAVBYHBwEIVCgWIo/1ahVGBlZYWXpxdeXl5IJNJTPnAqKytJS08zulQsX76CFctXXDCJGhcal8Q5Duh0Onb9vIv0tDS+/uZrkpOTWbpkKVdccSUTJ0485XpBEFAqlUax9vcPYGdna0iId3YeSnh3xt7e3uQZ1lRxajQalEoFvQoFil6DEDVaDXZ2dkicJUikUiQSCRJnySlLdUEQKC8vZ8eOHaRvS6emppqlS5ayfPmKcXVC+LVwSZzjTH9/P9988zWff/E5e/bsQSaTMXfOXObOncucOXNHFKtarUapVBoT3pV9fahUA4hEVkOtGByxsbEZSmy3xtra8CgWi0/62vC4/8A+ozgFQUCn06HRqNFotMc9atBoNUOzYi8DAypsrK0N4pNIkEgMQhypImZYjHv27GbP3j3s2bMHGxsbLpt3GYuXLOHKK660SMOp/xUuifMcotfrKSkpYc+ePezZu5u9e/eOSqwnP0dfXx/9/X1otFpDYrtWi3YowV07/D2dFq3W8KjT6VAoFEicJQgIiBAhthYPJbsP/bH+5WtbWzskEslZZ+qRxDhv7jzmzp3L3LnzCAoKujRDmsglcZ5HTidWqVRKTEwMYWHhhIeFGR7Dw/H09DTrJjdnz6nX62lqaqK8opyK8nLKKyqoqCgnLy8Pa2vrS2IcJy6J8wJCr9dz9OhRjh4tMQhgSAjl5WXo9fpfBBtueAwICEQ6tOdzdnZGIpGMuGwcSZzDgR2lUolCoaCnp4eammqjACvKK6iorEAqlZ7y+lFR0QQHB18S4zhxSZwXCXK5nPLycoNgKiupqCintrYOpfIXYSmVSqytrY1ClThLcHJ2RiJxprOzE7HYmr4+5VCU1XC9IAjGPaWzszNSqYyQ4GCjAMPCwgkNDb3oC5cvRkwS5yUuTESGKcwekADOJz3qAQWgHHpUAEpBEAbPz2gvYSqXxHmJS1yg/G+k/V/iEhchl8R5iUtcoFwS5yUucYFySZyXuMQFyiVxXuISFyj/D4OB+hXGC8/3AAAAAElFTkSuQmCC\n",
Maximilian Schanner's avatar
Maximilian Schanner committed
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from corbass.evaluation import location\n",
    "lat, lon, dens = location(posterior, mu_coeffs[:, 0:3], cov_coeffs[:, 0:3, 0:3])\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "ax_proj = ccrs.Stereographic(central_latitude=90., scale_factor=30)\n",
    "fig, ax = plt.subplots(1, 1, subplot_kw={'projection': ax_proj})\n",
206
207
208
209
210
    "pltlat, pltlon, _ = ax_proj.transform_points(ccrs.Geodetic(), lon, lat).T\n",
    "ax.coastlines(zorder=1);\n",
    "ax.gridlines();\n",
    "ax.set_global();\n",
    "ax.tripcolor(pltlat, pltlon, dens, zorder=0, cmap='Purples');"
Maximilian Schanner's avatar
Maximilian Schanner committed
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
230
   "version": "3.7.6"
Maximilian Schanner's avatar
Maximilian Schanner committed
231
232
233
234
235
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}