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:214: 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": {
79
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD4CAYAAAD//dEpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAbpUlEQVR4nO3df5BV5Z3n8fdnG/mxUVABE0KDdABj0GAbbxrDjoYpEJjMVDCJrIyZksxQhZjgZmr/SSyrwkaTSkxqy4nDRNOzZlAjK8ZUBkrjGhijZCIhNIqIUcYmRLlCxZYGQqIwNvnuH/dpvd3c7j707b6X7v68qk55zvec597nlt18+pzn3PMoIjAzM+vJf6l2B8zMbGBwYJiZWSYODDMzy8SBYWZmmTgwzMwsk2HV7kB/GTduXEyZMqXa3TAzG1C2b9/+RkSML7Wv3wND0m+Bo8AJoC0icpLOBdYBU4DfAv89Ig6l428GlqXj/0dEPJ7qlwFrgFHAT4AvRjf3BE+ZMoWmpqb++VBmZoOUpFe62lepS1J/HhH1EZFL218G/i0ipgP/lraRNANYAlwELAS+K6kmtbkLWA5MT8vCCvXdzMyo3hjGIuDetH4vcHVR/cGIOB4Re4FmoEHSBGB0RGxJZxX3FbUxM7MKqERgBPBTSdslLU+190bEAYD03/NSfSKwr6htPtUmpvXO9Q4kLZfUJKmppaWljz+GmdnQVolB7/8WEfslnQdslPRSN8eqRC26qXcsRDQCjQC5XM7PPDEz60P9foYREfvTf18Hfgw0AL9Ll5lI/309HZ4HJhU1rwX2p3ptibqZmVVIvwaGpPdIOqt9HZgP7AI2AEvTYUuB9Wl9A7BE0ghJdRQGt3+VLlsdlXS5JAHXF7UxM7MK6O9LUu8Fflz4N55hwNqI+H+StgEPSVoGvAosBoiIFyQ9BPwaaAO+EBEn0mvdyLu31T6WFjMzqxAN1seb53K58PcwzMxOjaTtRV+B6MCPBjEzOw1d+70tXPu9LdXuRgcODDOzweRf/rKw9AMHRgmnY7KbmVWbA8PM7DRy91N7eHrPGx1qT+95g7uf2tN9w3//B9i7uWNt7+ZCvY84MPpSP54KmtnQMLN2DCvXPsuRt94GCmGxcu2zzKwd033DiR+BH34Ojh0ubO/dXNie+JE+65sDo8jpnOxmNjTMnjqO1dddSvPrfyB/6E1Wrn2W1dddyuyp47pvWHclLF4DLS/B4VcKYbF4TaHeRxwYRU7nZDezoWP21HG8d/QIXjt8jL+ZNbnnsGhXdyWcNQGO7IPcsj4NC3BgdHA6J7uZDR1P73mD3/3+OBPPHskPtr560pWPLu3dDEcPwJhJ0HTPyVc+yuTA6OR0TXYzGxrar2xMO+9Mas/5r6y+7lJWrn2259Bov7Ix/kI4+/zCH60//FyfhoYDo5PTNdnNbGjYmT/C6usuZcyoM4B3r3zszB/pvuFrzxRCYuTZhe32Kx+vPdNnfXNgFDmdk93MhoYVH5960pWN2VPHseLjU7tv+Gd/f/KVjborC/U+4sAocjonu5kNLetu+BjrbvjYqTf820cLSz/wwwdLaP+Wd6/+Z5mZDWB++KCZmZWtElO0Djg+szAzO5nPMMzMLJMBFRiSFkraLalZ0per3R8zs6FkwASGpBrgn4C/AGYAfy1pRnV7ZWY2dAyYwAAagOaI+E1E/CfwILCoyn0yMxsyBlJgTAT2FW3nU+0dkpZLapLU1NLSUtHOmZkNdgMpMFSi1uFLJBHRGBG5iMiNHz++Qt0yMxsaBlJg5IFJRdu1wP4q9cXMbMgZSIGxDZguqU7ScGAJsKHKfTIzGzIGzBf3IqJN0krgcaAG+H5EvFDlbpmZDRkDJjAAIuInwE+q3Q8zs6FoIF2SMjOzKnJgmJlZJg4MMzPLxIFhZmaZODDMzCwTB4aZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLxIFhZmaZODDMzCwTB4aZmWXiwDAzs0z6LTAk/S9Jr0nakZZPFO27WVKzpN2SFhTVL5P0fNp3pySl+ghJ61J9q6Qp/dVvMzMrrb/PMO6IiPq0/ARA0gwK06teBCwEviupJh1/F7AcmJ6Wham+DDgUEdOAO4Db+7nfZmbWSTUuSS0CHoyI4xGxF2gGGiRNAEZHxJaICOA+4OqiNvem9YeBue1nH2ZmVhn9HRgrJe2U9H1J56TaRGBf0TH5VJuY1jvXO7SJiDbgCDC285tJWi6pSVJTS0tL334SM7MhrqzAkLRJ0q4SyyIKl5emAvXAAeB/tzcr8VLRTb27Nh0LEY0RkYuI3Pjx40/585iZWdeGldM4IuZlOU7SPwOPpM08MKlody2wP9VrS9SL2+QlDQPGAK2977mZmZ2q/rxLakLR5qeAXWl9A7Ak3flUR2Fw+1cRcQA4KunyND5xPbC+qM3StH4N8EQa5zAzswop6wyjB9+SVE/h0tFvgRsAIuIFSQ8BvwbagC9ExInU5kZgDTAKeCwtAPcA90tqpnBmsaQf+21mZiVosP6hnsvloqmpqdrdMDMbUCRtj4hcqX3+preZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLxIFhZmaZODDMzCwTB4aZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLxIFhZmaZlDun92JJL0j6k6Rcp303S2qWtFvSgqL6ZZKeT/vuTLPrkWbgW5fqWyVNKWqzVNLLaVmKmZlVXLlnGLuATwObi4uSZlCYFe8iYCHwXUk1afddwHIKU7NOT/sBlgGHImIacAdwe3qtc4FVwCygAVgl6Zwy+21mZqeorMCIiBcjYneJXYuAByPieETsBZqBhjTP9+iI2JLm5L4PuLqozb1p/WFgbjr7WABsjIjWiDgEbOTdkDEzswrprzGMicC+ou18qk1M653rHdpERBtwBBjbzWudRNJySU2SmlpaWvrgY5iZWbthPR0gaRPwvhK7bomI9V01K1GLbuq9bdOxGNEINEJhTu8u+mZmZr3QY2BExLxevG4emFS0XQvsT/XaEvXiNnlJw4AxQGuqz+nU5sle9MnMzMrQX5ekNgBL0p1PdRQGt38VEQeAo5IuT+MT1wPri9q03wF1DfBEGud4HJgv6Zw02D0/1czMrIJ6PMPojqRPAf8IjAcelbQjIhZExAuSHgJ+DbQBX4iIE6nZjcAaYBTwWFoA7gHul9RM4cxiCUBEtEq6DdiWjrs1IlrL6beZmZ06Ff6IH3xyuVw0NTVVuxtmZgOKpO0RkSu1z9/0NjOzTBwYZmaWiQPDzMwycWCYmVkmDgwzM8vEgWFmZpk4MMzMLBMHhpmZZeLAMDOzTBwYZmaWiQPDzMwycWCYmVkmDgwzM8vEgWFmZpk4MMzMLBMHhpmZZVJWYEhaLOkFSX+SlCuqT5H0lqQdabm7aN9lkp6X1CzpzjRVK2k613WpvlXSlKI2SyW9nJalmJlZxZV7hrEL+DSwucS+PRFRn5YVRfW7gOUU5vmeDixM9WXAoYiYBtwB3A4g6VxgFTALaABWpbm9zcysgsoKjIh4MSJ2Zz1e0gRgdERsicLcsPcBV6fdi4B70/rDwNx09rEA2BgRrRFxCNjIuyFjZmYV0p9jGHWSnpX0lKQrUm0ikC86Jp9q7fv2AUREG3AEGFtcL9GmA0nLJTVJamppaem7T2JmZgzr6QBJm4D3ldh1S0Ss76LZAWByRByUdBnwr5IuAlTi2Gh/qy72ddemYzGiEWgEyOVyJY8xM7Pe6TEwImLeqb5oRBwHjqf17ZL2ABdQODuoLTq0Ftif1vPAJCAvaRgwBmhN9Tmd2jx5qn0yM7Py9MslKUnjJdWk9Q9QGNz+TUQcAI5KujyNT1wPtJ+lbADa74C6BngijXM8DsyXdE4a7J6famZmVkE9nmF0R9KngH8ExgOPStoREQuAK4FbJbUBJ4AVEdGamt0IrAFGAY+lBeAe4H5JzRTOLJYARESrpNuAbem4W4tey8zMKkSFP+IHn1wuF01NTdXuhpnZgCJpe0TkSu3zN73NzCwTB4aZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLxIFhZmaZODDMzCwTB4aZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLpKzAkPRtSS9J2inpx5LOLtp3s6RmSbslLSiqXybp+bTvzjTzHpJGSFqX6lslTSlqs1TSy2lZipmZVVy5ZxgbgYsjYibwH8DNAJJmUJgx7yJgIfDd9ilbgbuA5RSmbZ2e9gMsAw5FxDTgDuD29FrnAquAWUADsCpN1WpmZhVUVmBExE8joi1t/hKoTeuLgAcj4nhE7AWagQZJE4DREbElzdd9H3B1UZt70/rDwNx09rEA2BgRrRFxiEJItYeMmZlVSF+OYfwd787PPRHYV7Qvn2oT03rneoc2KYSOAGO7ea2TSFouqUlSU0tLS1kfxszMOhrW0wGSNgHvK7HrlohYn465BWgDHmhvVuL46Kbe2zYdixGNQCMU5vQudYyZmfVOj4EREfO6258Gof8KmJsuM0HhLGBS0WG1wP5Ury1RL26TlzQMGAO0pvqcTm2e7KnfZmbWt8q9S2oh8CXgkxHxZtGuDcCSdOdTHYXB7V9FxAHgqKTL0/jE9cD6ojbtd0BdAzyRAuhxYL6kc9Jg9/xUMzOzCurxDKMHq4ERwMZ0d+wvI2JFRLwg6SHg1xQuVX0hIk6kNjcCa4BRFMY82sc97gHul9RM4cxiCUBEtEq6DdiWjrs1IlrL7LeZmZ0ivXsVaXDJ5XLR1NRU7W6YmQ0okrZHRK7UPn/T28zMMnFgmJlZJg4MMzPLxIFhZmaZODDMzCwTB4aZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLxIFhZmaZODDMzCwTB4aZmWXiwDAzs0zKnXHv25JekrRT0o8lnZ3qUyS9JWlHWu4uanOZpOclNUu6M828R5qdb12qb5U0pajNUkkvp2Vp536YmVn/K/cMYyNwcUTMBP4DuLlo356IqE/LiqL6XcByCtO2TgcWpvoy4FBETAPuAG4HkHQusAqYBTQAq9JUrWZmVkFlBUZE/DQi2tLmL4Ha7o6XNAEYHRFb0nzd9wFXp92LgHvT+sPA3HT2sQDYGBGtEXGIQkgtxMzMKqovxzD+jnfn5waok/SspKckXZFqE4F80TH5VGvftw8ghdARYGxxvUQbMzOrkGE9HSBpE/C+ErtuiYj16ZhbgDbggbTvADA5Ig5Kugz4V0kXASrxOu2Tine1r7s2nfu6nMLlLiZPnlz6A5mZWa/0GBgRMa+7/WkQ+q+AuekyExFxHDie1rdL2gNcQOHsoPiyVS2wP63ngUlAXtIwYAzQmupzOrV5sou+NgKNALlcrmSomJlZ75R7l9RC4EvAJyPizaL6eEk1af0DFAa3fxMRB4Cjki5P4xPXA+tTsw1A+x1Q1wBPpAB6HJgv6Zw02D0/1czMrIJ6PMPowWpgBLAx3R37y3RH1JXArZLagBPAiohoTW1uBNYAoyiMebSPe9wD3C+pmcKZxRKAiGiVdBuwLR13a9FrmZlZhShdRRp0crlcNDU1VbsbZmYDiqTtEZErtc/f9DYzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLxIFhZmaZODDMzCwTB4aZmWXiwDAzs0wcGGZmlokDw8zMMnFgmJlZJg4MMzPLxIFhZmaZODDMzCwTB4aZmWVS7hStt0naKWmHpJ9Ken/RvpslNUvaLWlBUf0ySc+nfXemqVqRNELSulTfKmlKUZulkl5Oy1LMzKziyj3D+HZEzIyIeuAR4CsAkmZQmGL1ImAh8N32Ob6Bu4DlFOb5np72AywDDkXENOAO4Pb0WucCq4BZQAOwKs3tbWZmFVRWYETE74s23wO0z/e6CHgwIo5HxF6gGWiQNAEYHRFbojA37H3A1UVt7k3rDwNz09nHAmBjRLRGxCFgI++GjJmZVciwcl9A0teB64EjwJ+n8kTgl0WH5VPt7bTeud7eZh9ARLRJOgKMLa6XaNO5L8spnL0wefLkXn8mMzM7WY9nGJI2SdpVYlkEEBG3RMQk4AFgZXuzEi8V3dR726ZjMaIxInIRkRs/fnx3H8vMzE5Rj2cYETEv42utBR6lMN6QByYV7asF9qd6bYk6RW3ykoYBY4DWVJ/Tqc2TGftkZmZ9pNy7pKYXbX4SeCmtbwCWpDuf6igMbv8qIg4ARyVdnsYnrgfWF7VpvwPqGuCJNM7xODBf0jlpsHt+qpmZWQWVO4bxTUkfBP4EvAKsAIiIFyQ9BPwaaAO+EBEnUpsbgTXAKOCxtADcA9wvqZnCmcWS9Fqtkm4DtqXjbo2I1jL7bWZmp0iFP+IHn1wuF01NTR1qb7/9Nvl8nmPHjlWpV5U1cuRIamtrOeOMM6rdFTMbICRtj4hcqX1l3yU1kOTzec466yymTJlC+r7goBURHDx4kHw+T11dXbW7Y2aDwJB6NMixY8cYO3bsoA8LAEmMHTt2yJxNmVn/G1KBAZxyWFz7vS1c+70t/dSb/jUUgtHMKmfIBYaZmfWOA6MLdz+1h6f3vNGh9vSeN7j7qT1lvW5NTQ319fVcfPHFLF68mDfffLOs1zMzqxQHRhdm1o5h5dpnOfLW20AhLFaufZaZtWPKet1Ro0axY8cOdu3axfDhw7n77rsztz1x4kTPB5mZ9RMHRhdmTx3H6usupfn1P5A/9CYr1z7L6usuZfbUcX32HldccQXNzc0A/OAHP6ChoYH6+npuuOGGd8LhzDPP5Ctf+QqzZs1iy5YtbNu2jdmzZ3PJJZfQ0NDA0aNH+6w/ZmbdcWB0Y/bUcbx39AheO3yMv5k1uU/Doq2tjccee4wPf/jDvPjii6xbt45f/OIX7Nixg5qaGh544AEA/vjHP3LxxRezdetWGhoauPbaa/nOd77Dc889x6ZNmxg1alSf9cnMrDtD6nsYp+rpPW/wu98fZ+LZI/nB1le5fOrYskPjrbfeor6+HiicYSxbtozGxka2b9/ORz/60XeOOe+884DCmMdnPvMZAHbv3s2ECRPeOW706NFl9cXM7FQ4MLrQPmYx7bwzGTPqDL44b3qfXJZqH8MoFhEsXbqUb3zjGycdP3LkSGpqat45zrfKmlm1+JJUF3bmj7D6uksZM6rwWI32MY2d+SN9/l5z587l4Ycf5vXXXwegtbWVV1555aTjLrzwQvbv38+2bYXHah09epS2trY+74+ZWSkOjC6s+PjUk84kZk8dx4qPT+3z95oxYwZf+9rXmD9/PjNnzuSqq67iwIEDJx03fPhw1q1bx0033cQll1zCVVddxbFjx9i/fz+f+MQn+rxfZmbFhtTDB1988UU+9KEPValH1TEUP7OZ9V53Dx/0GYaZmWXiwDAzs0zKnXHvNkk7Je2Q9FNJ70/1KZLeSvUdku4uanOZpOclNUu6M828R5qdb12qb5U0pajNUkkvp2Vp536YmVn/K/cM49sRMTMi6oFHgK8U7dsTEfVpWVFUvwtYTmHa1unAwlRfBhyKiGnAHcDtAJLOpTBP+CygAViVpmo1O231+inH//KXhcXsNFRWYETE74s23wN0O4IuaQIwOiK2pPm67wOuTrsXAfem9YeBuensYwGwMSJaI+IQsJF3Q6b/+RfYzAzogzEMSV+XtA/4LB3PMOokPSvpKUlXpNpEIF90TD7V2vftA4iINuAIMLa4XqJN574sl9QkqamlpaXMT2Z26nr9lON//wfYu7ljbe/mQt3sNNFjYEjaJGlXiWURQETcEhGTgAeAlanZAWByRFwK/E9graTRQKmvKbeflXS1r7s2HYsRjRGRi4jc+PHje/po3evHX+Cvf/3rXHTRRcycOZP6+nq2bt3KnDlz+OAHP8jMmTO58MILWblyJYcPHy77vayyev2U44kfgR9+Do6l/+d7Nxe2J36kX/trdip6DIyImBcRF5dY1nc6dC3wmdTmeEQcTOvbgT3ABRTODmqL2tQC+9N6HpgEIGkYMAZoLa6XaNN/+ukXeMuWLTzyyCM888wz7Ny5k02bNjFpUuHjPfDAA+zcuZOdO3cyYsQIFi1aVN5nsIrr9VOO666ExWug5SU4/ErhZ23xmkLd7DRR7l1S04s2Pwm8lOrjJdWk9Q9QGNz+TUQcAI5KujyNT1wPtAfPBqD9DqhrgCfSOMfjwHxJ56TB7vmp1r/66Rf4wIEDjBs3jhEjRgAwbtw43v/+93c4Zvjw4XzrW9/i1Vdf5bnnnivr/azyev2U47or4awJcGQf5JY5LOy0U+4YxjfT5amdFP4h/2KqXwnslPQchQHsFRHRmvbdCPwfoJnCmcdjqX4PMFZSM4XLWF8GSO1uA7al5dai1+pf/fALPH/+fPbt28cFF1zA5z//eZ566qmSx9XU1HDJJZfw0ksvlf2eVlmdn3LceUyjS3s3w9EDMGYSNN1z8iVRsyor62m1EfGZLuo/An7Uxb4m4OIS9WPA4i7afB/4fu972kudf4Hrrig7NM4880y2b9/Oz3/+c372s59x7bXX8s1vfrPksYP1sS2DWa+fctx+yXP8hTDybJjzJV+WstOOv+ndleJf4LPPL/zi/vBzffJXX01NDXPmzOGrX/0qq1ev5kc/OjlbT5w4wfPPP+/nQA0wvX7K8WvPFH7GRp5d2G6/JPraM/3aX7NT4cDoSj/9Au/evZuXX375ne0dO3Zw/vnndzjm7bff5uabb2bSpEnMnDmzrPezyur1U47/7O9PPpOou7JQNztNeAKlrrT/oj55+7u1uivLvjzwhz/8gZtuuonDhw8zbNgwpk2bRmNjI9dccw2f/exnGTFiBMePH2fevHmsX9/5RjQb9P720Wr3wKxLfrz5IDcUP7OZ9Z4fb25mZmVzYJiZWSZDLjAG6yW4UobSZzWz/jekAmPkyJEcPHhwSPxDGhEcPHiQkSNHVrsrZjZIDKm7pGpra8nn8wyVJ9mOHDmS2trang80M8tgSAXGGWecQV1dXbW7YWY2IA2pS1JmZtZ7DgwzM8vEgWFmZpkM2m96S2oBXumnlx8HZHxmtVlJ/hmycvXXz9D5EVFyytJBGxj9SVJTV1+dN8vCP0NWrmr8DPmSlJmZZeLAMDOzTBwYvdNY7Q7YgOefIStXxX+GPIZhZmaZ+AzDzMwycWCYmVkmDoxekrRY0guS/iTJt0daJpIWStotqVnSl6vdHxt4JH1f0uuSdlX6vR0YvbcL+DSwudodsYFBUg3wT8BfADOAv5Y0o7q9sgFoDbCwGm/swOiliHgxInZXux82oDQAzRHxm4j4T+BBYFGV+2QDTERsBlqr8d4ODLPKmQjsK9rOp5rZgDCk5sM4VZI2Ae8rseuWiFhf6f7YgKcSNd/XbgOGA6MbETGv2n2wQSUPTCrargX2V6kvZqfMl6TMKmcbMF1SnaThwBJgQ5X7ZJaZA6OXJH1KUh74GPCopMer3Sc7vUVEG7ASeBx4EXgoIl6obq9soJH0f4EtwAcl5SUtq9h7+9EgZmaWhc8wzMwsEweGmZll4sAwM7NMHBhmZpaJA8PMzDJxYJiZWSYODDMzy+T/A+JmUoKCrsW4AAAAAElFTkSuQmCC\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": {
116
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD4CAYAAAANbUbJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxV9Z3/8dc3O1khOxAEkgAhkX0RgoIViYJirU7V0WrdsVOcbr/OqO2002UerR1ta0unuC9VrFbajohaEAW1UNYoS9gS1kDIBgSy39z7/f2RwCCCstx7z8097+fjwYPkEu55P0Ly5uR7v+dzjLUWEREJbxFOBxARkcBT2YuIuIDKXkTEBVT2IiIuoLIXEXGBKKcDnEp6erodMGCA0zFERLqVtWvX1llrM071ZyFZ9gMGDGDNmjVOxxAR6VaMMbtP92chtYxjjJlpjHmioaHB6SgiImElpMreWrvAWntvSkqK01FERMJKSJW9iIgEhspeRMQFVPYiIi6gshcRcQGVvYiIC6jsRURcIKTKXvvsRcTNbnx8BTc+viIgzx1SZa999iIigRFSZS8iIoGhshcRcQGVvYiIC6jsRURcQGUvIuICKnsRERdQ2YuIuEBIlb0uqhIRCYyQKntdVCUiTgjklauhIqTKXkREAkNlLyLiAip7EREXUNmLiLiAyl5ExAVU9iIiLqCyFxFxAZW9iIgLqOxFRBw2d1kFyyvqPvHY8oo65i6r8NsxVPYi4hg3XLl6JobnpDB7Xin1jW20erwsr6hj9rxShuf4b5pAlN+eyQ+MMTOBmfn5+U5HEREJuPYOH2t3H+L9bXUkxERSXttEfHQEs+eVMufmURTnpfvtWCFV9tbaBcCCsWPH3uN0FhGRQKg81MyybbUs21rL8op6Gts6iIowjB3Qi4YWD0daO7j7kgv8WvQQYmUvIhJuWj1eVu86yNKttSzbVkt5TSMAfXv24JqRfZgyOIPivDQ27GvgtqdX0bdnHC+u3MOEvLTwPbMXEQkHu+qaOs/et9WyoqKeFo+XmMgILspN5aZx/bh0SAZ5GYkYYwCOr9HnZyaS0iOab1w+yO9LOSp7EZHz1NLu5R876lm6tYZl22rZVd8MQP+0eG4Ym8OUIRlMyE0jPubUlbu+soE5N4/isXe2A1Ccl86cm0exvrJBZS8i4hRrLRW1jceXZlbuPEh7h4+46Agm5qZxx6SBTBmcwYD0hDN6vvum5AEcL3voLHwt44iIBFljWwfLy+tY2vXi6r7DLQDkZyZy64T+XDokg3EDUomLjnQ46amp7EXEteYuq/jUXvblFXWsr2xg1uRcthw4yrJttSzdWsPa3YfweC0JMZFMyk/nX76Qx+RBGfRLjXco/dlR2YuIax27mCkzKZaUHtEsLqvmW698xNj+PXn27zupPtIGQEF2EnddnMuUwRmM6d+LmKjudz2qyl5EXKs4L50fXF3It1/9iMgIwz0vrAFg3Z7DXDIogymDM5gyJIOs5DiHk54/lb2IuM6e+mYWbqhi4Yb9bNx3BACf1zJuQC/+/coCRvbrSVRk9zt7/ywqexFxhVMV/Ih+Pbnlogv44+q9ZCfHUlHbRLvXF3ZFDyp7EVc6NnzslVkTHU4SWLvrm1i4oYo3N1R9ouC/N2Mo04dls+dgM7PnlTIogBczhQqVvYiElc8r+Jxe/7d75o31VQG/mClUqOxFpNs7m4I/UTAuZgoVIVX2GnEsImfqXAverUKq7DXiWEQ+y6kKfqQK/oyEVNmLiJzsWMEvXF/Fpv0q+HOlsheRoPusMQX3TclTwQeAyl5Egu7kMQXLK+r42ovruPLCbK76zQcq+ABQ2YtI0B3b4njr06tIiInk1qdX4fVZXlm9VwUfICp7EQmqo60e3tpwgNfWVeL1WY60dpCVHMvdF+eq4ANIZS8iAef1Wf5eXsf8dZX8bdMBWj0+spPjiDCQkRiLx2sp6pusog8glb1IkLllVAHAtuqjzF9byV8/2kf1kTaS46K4fnQOg7OTeOyd7QzOSgr7MQWhQmUvIn5V39jG6x/vZ/66SjbuO0JkhOELQzL44cwcLivIJC46krnLKlwzpiBUqOxF5Ly1dXh5d3MN89ftY+nWGjp8lqI+yfzg6kKuGdmH9MTYT3y8m8YUhAqVvYicE2stH+09zPx1lSz4uIqGFg+ZSbHcdfFAvjS6LwXZyU5HlBOo7EXkrOw73MJfS/cxf10lO2qbiI2K4IqibK4fk8OkvLSwnAUfDlT2IvK5mto6eGvjAf68rpIVO+qxFsYPTGXW5FxmDOtNUly00xHlc6jsReSUvD7LP3bUM39tJW9tPECLx0v/tHi+OXUw143uS79UbZPsTlT2IvIJ5TWNzF9XyV9L91HV0EpSXBTXjurL9aP7MqZ/L4wxTkeUc6CyF3GR0w0gW7njIGmJMcxft4+P9x4mMsIweVA637tqKJcPzSIuOtKhxOIvKnsRFzlxAFlSXBS/WbKd3767HZ+1eH0wtHcy379qKNeM7ENmUpzTccWPVPYiLjJ+QCr3Tc7lZ29vAQurdx0iOS6aG8bmcN3oHAr7aLtkuFLZi4Q5j9fHiop63txQxd82HeBQs+f4n10zog+/vGGEtkuGiECO0FDZi4Qhj9fH38vreHNDFYvKqjnc7CEhJpLLC7PIy0jkN0u2k5Ucy4fldazadVBXrrqAyl4kTLR3dBb8wg1VLC6rpqHFQ1JsFJcXZjH9wmwmD85g3Z5DzJ5XSn5mogaQuUxIlb0xZiYwMz8/3+koIt1CW4eXD7f/X8Efbe0gKTaKaYVZzBjWm0sGpxMb9X87adZXNmgAmUuFVNlbaxcAC8aOHXuP01lEQlWrx8sH2+t4a0MVizd3FnxyXBQlhdlcNTybSfmfLPgTaQCZe4VU2YvIqbV6vLy/rZY3N1TxzuYaGts6SOkRzZVF2cwY3ptJeenEROlFVjk9lb1IiGr1eFm6tbPgl2yupqndS8/4aK4a1psZw3tTnJdGtHbRyBlS2YuEkJZ2L0u31rBwQxXvbqmhud1Lr/horhnZh+kX9maiCl7OkcpeJEhON6pg7a5D5GYk8mZXwbd4vKQmxPDFkX25alhvJuSmah+8nLewKns33dtTzp7TXx8njipIjI3i0UVbeXxZBRbweC3piTFcN7qz4McPVMGLf4VV2YuEsuK8dL5bMoSH/rIBC6zZ3Tmq4Isj+zCjq+AjIzRR0gluOEFU2YsEwf7DLTy6aBt/Lq08/tj1o/vyi38aoYKXoAiLnxPnLqtgeUXdJx5bXlHH3GUVDiUS6dTQ4uHnb23h0keWsmD9fq4a1pvICEPfnnG8t7WWlTvrnY4oLhEWZX9sLbShuZ0jrR7e31bL7Hmln3oxTCRY2jq8PP3hTqb893s8/n4FVw/vzSNfHsHyinryMxPJ6RXPnJtHMXte6adOVEQCISzK/tgl39trGtlcdZSvz1unWR/iCJ/P8vrH+7n8l8v4yRtlDOubwhv3X8wvbxjJ/sMtzLl5FCk9Ou/XeuKoApFAC5s1++K8dDKT46hqaGVgWoKKXoJuRUU9P3trM+srGxjaO5kX7hzG5MEZx/9cowrESWFT9ssr6qg92kZcVAQb9zewvLyO4nx9E0ngbas+ys/f2sK7W2rokxLHo18ewbWj+uqFVwkpYVH2yyvqjo9tbe/wsaOuifteXMvcW8forEkC5kBDK79avI0/rd1LQmwUD0wv4PbiAbpfq4SksCj7E8e2erw+IgxMHZqpsa0SEEdbPTy+bAdPfbgDr89yx6SBzP5CPr0SYpyOJnJaYVH2J66FRkdGMG5AKmX7j/KrG0c5nEzCSXuHj5dX7eGxJds52NTONSP68N0rhtAvNd7paCKfKyzK/phjV8E99cEOfrpwM7vrm+ifluBwKunurLW8tfEAv3h7C7vqm5mQm8pDM4YyPKen09FEzlhYbL08WUlhNgCLy6odTiLd3epdB7nu98v5l5fWERMVwbO3j+Pleyao6KXbCasz+2MuSIunIDuJRWXV3H1JrtNxpBsqr2nk4be3sLismqzkWH5x/XCuH5OjHTZ+5oaZNKEiLMseoKQwiznvlVPf2EZaYqzTcaSbqDnaymPvbOePq/fSIzqS714xhDsnDaRHjHbYSPcWvmVflM1v3i1nyZYabhjbz+k4EuKa2jp48oMdPPH+Dto7fNw6oT/3X5avEwUJG2Fb9kV9kumTEsfismqVvZxWh9fHK2v28qvF26lrbOOqYb357hVDGJCuF/YlvIRt2RtjmFaYxStr9tLS7tWP4fIJ1loWl1Xz87e3sKO2iXEDevHEbWMYfUEvp6OJBETYlj10LuU8v2I372+v5YqibKfjiENOvh3guj2HeHD+erZWN5KXkcCTt43l8qGZGKMXXyV8hXXZjx+YSnJcFIvLqlX2LnZsBHavHtEcam7nuv9ZjgHuunggD04v0O3/xBXC+qs8OjKCywoyWbK5mg6vz+k44pDivHS+PW0QFXVNHGz20CM6kqdvH8t/XF2oohfXCPuv9JKibA41e1iz+5DTUcQhe+qb+fU724nsWqW555KBXFaQ5WwokSAL+7KfPDiDmKgIXU3rUvWNbXz12VW0eLxgOm8H+OLKPbo7lLhO2Jd9YmwUk/LSWFR2AGut03EkiFravdz1/BoqDzYTaQyDdDtAcbGwL3voXMrZe7CFLQeOOh1FgqTD6+P+l9fxceVhrh7Rh7m3jtHtAMXVXFH2U4dmYowGo7mFtZYfvL6JdzbX8KNrivjVjSM/dV+D4rz046OxRdzAFWWfmRTHqH49WVR2wOkoEgS/e6+ceSv38LVL87ht4gCn44iEBL+XvTEm1xjztDHmtc96LNhKirLZuO8I+w63OBVBguC1tZU8smgbXxrVl3+7YojTcULWK7MmauKky5xR2RtjnjHG1BhjNp70+JXGmK3GmHJjzAMA1tod1tq7Tvy4Uz0WbNMKO7favaOlnLD1/rZaHpi/novz03n4+uG6IlbkBGd6Zv8ccOWJDxhjIoHfAdOBQuCfjTGFfk3nR3kZieRlJGgpJ0xt3NfA115cy6CsJH7/ldHERLlihVLkjJ3Rd4S19n3g4EkPjwfKu87a24E/Al881yDGmHuNMWuMMWtqa2vP9Wk+U0lRNit3HKSh2ROQ5xdn7D3YzB3PraZnfAzP3TGOpLhopyOJhJzzOf3pC+w94f1KoK8xJs0YMxcYZYx5EOBUj53MWvuEtXastXZsRkbGecQ6vWmFWXT4LO9trQnI80vwHWpq56vPrqLN4+W5O8aRlRzndCSRkHQ+g9BOtSBqrbX1wH0nPfipx5wwMqcnGUmxLCo7wLWj+jodR85Tq8fLXc+vpvJQCy/edRGDspKcjiQSss7nzL4SOPGuIDnA/vOLE1gREZ0z7pdtraXV43U6jpwHr8/yry+XUrr3MI/dOJLxA1OdjiQS0s6n7FcDg4wxA40xMcBNwOv+iRU40wqzaGr3sqKi3ukoco6stfxowSYWlVXzg6sLmT6st9ORRELemW69fBlYAQwxxlQaY+6y1nYAs4G/AZuBV621mwIX1T+K89JIiInUrpxubO6yHbywYjf3Ts7ljkkDnY4j0i2c0Zq9tfafT/P4m8Cb/gpjjJkJzMzPz/fXU35KbFQklxZksrishv+61hIRob3YgXbj4ysA/HIRz19L9/Hw21uYOaIPD1xZcN7PJ+IWIbUZ2Vq7wFp7b0pKyud/8HkoKcyirrGN0r2HA3oc8a+/l9fx3dc+ZmJuGo98ebj+oxY5C2F9W8LTuXRIJlERhkVlBxjTXzeY7g7K9h9h1h/WkpeRyNxbxxAb1X1vIK8xBeKEkDqzD5aUHtFMzEvTFMxuovJQM7c/u4qkuCievWPc8VHFInLmXFn20LkrZ0dtE+U1jU5Hkc9wuLmd259dTYvHy/N3jqd3Sg+nI4l0S64t+8uHdg5G066c0NXq8XLvC2vZU9/Mk7eNZbAumhI5ZyFV9saYmcaYJxoaAn8HoT49ezA8J0VLOSHK57N8+9WPWLXrII/eMIIJuWlORxLp1kKq7IO1G+eYaUOzKN1zmJojrUE5npwZay0/WVjGmxsO8P2rhjJzRB+nI4l0eyFV9sFWUpQNwOLNOrsPJU99sJNn/76LOycN5O5Lcp2OIxIWXF32g7MS6Z8Wr6WcEPL6x/v5rzc3c9Xw3nz/qqFOxxEJG64ue2MM04Zmsby8nsa2DqfjuN7yijq+8+pHjB+YyqNfHqGLpkT8yNVlD51LOe1eH8u2BuaGKXJmthw4wqwX1jIwPYEnbx1LXHT3vWhKJBS5vuzH9O9FakKMtmA6aP/hFm5/ZjXxsZE8d8d4UuJ10ZSIv4VU2Qdz6+UxkRGGqQWZvLulBo/XF7TjSqeGFg+3P7uKprYOnrtjPH166qIpkUAIqbIP9tbLY0qKsjna2sHKHSffZlcCqa3Dy70vrGFnXROP3zqGob2TnY4kErZCquydcnF+OnHREVrKCSKfz/KdVz9m5c6DPPLlERTnpzsdSSSsqeyBHjGRTB6UweKyaqy1TsdxhZ+9tZk31lfx4PQCvjgyOPcDfmXWRE2cFNdS2XcpKcqmqqGVjfuOOB0l7D394U6e/GAntxcP4N7JumhKJBhU9l0uK8gkwmgwWqAtXF/FTxeWcWVRNv9xdSHGaC+9SDCo7LukJsQwbkAqizbpatpAWbmjnm+98hFjLujFr28aSaQumhIJGpX9CUqKstlafZTd9U1ORwkLc5dVsLyiDoDm9g7ueWEN6YkxTMpP10VTIkEWUmXvxD77E5UUds6416wc/xiek8LseaXUN7ax9UAjxhia2r1clJvqdDQR1wmpsndqn/0x/VLjKchO0lKOnxTnpfPz64ZRXttEu9eHtZbff2U0xXnaZikSbCFV9qGgpCibNbsPUt/Y5nSUbu9Iq4c575Uff//24gEqehGHqOxPUlKYhc/Cki01Tkfp1prbO7jrudVs3NdApIG+PeN4ceWe42v4IhJcKvuTFPVJpk9KnJZyzkOrx8usP6xlza5DxMdEMSgriZxe8cy5eRSz55Wq8EUcoLI/iTGGkqJsPiyvpaXd63Scbsfj9TF73jo+2F7HjGG9eeK2MaT06JxiWZyXzpybR7G+0pkX4EXcTGV/CtMKs2j1+Hh/u2bcnw2vz/KtVz7inc01/OTaC/ndLZ9+MbY4L537puQ5lFDEvVT2pzB+YCrJcVHagnkWfD7LA/PX88b6Kh6aUcCtE/o7HUlETqCyP4XoyAimDs1iyeZqOjTj/nNZa/nRgk38aW0l35g6iHsn68xdJNSo7E9jWmEWh5o9rNl9yOkoIc1ay8Nvb+X5Fbu5d3Iu37x8kNORROQUQqrsnb6C9kSTB2cQExWhpZzPMefdcuYuq+ArEy7gwekFGmwmEqJCquydvoL2RImxUVycn86isgOacX8aT32wg0cXb+O60X358TUXquhFQlhIlX2omVaYxd6DLWw5cNTpKCFn3so9/HThZmYMy+YX1w8nQhMsRUKayv4zTB2aiTEajHayv5RW8r2/buCygkx+feMooiL1ZSQS6vRd+hkyk+IY1a9nt72hyY2Pr+DGx1f49Tnf3ljF//vTeibmpvE/t4wmJkpfQiLdgb5TP0dJUTYb9x1h3+EWp6M47r2tNdz/cikj+/XkydvGaia9SDeisv8cx2bcv+PypZwVFfXc94e1DMlO4pnbx5EQG+V0JBE5Cyr7z5GbkUheRkK3Xcrxh3V7DnHX86u5IDWeF+686PisGxHpPlT2Z6CkKJt/7DhIQ7PH6ShBt3FfA199ZhWZSbG8dPdFpCbEOB1JRM6Byv4MlBRm4fVZ3tvqrhn35TVHue2ZVSTHRfPSPRPITI5zOpKInCOV/RkYkdOTzKRYVy3l7K5v4uYnVxIZYXjp7ovo27OH05FE5Dyo7M9ARITh8sIslm6tpdUT/jPu9x9u4eYnV+Lx+njp7osYkJ7gdCQROU8hVfahNBvnZCWFWTS3e1lRUe90lICqOdrKLU+t5EiLhz/cdRGDs5KcjiQifhBSZR9Ks3FONjEvjcTYqLBeyjnU1M6tT62i+kgrz905jgv7ht6/g4icm5Aq+1AWGxXJlCEZLC6rwecLv8FoR1o93PbMKnbWN/HUbWMZ0z/V6Ugi4kcq+7NQUphFXWMbpXsPOx3Fr5rbO7jz2dVsOXCEuV8ZTXF++uf/JRHpVlT2Z+HSIZlERZiwWspp9Xi594W1rNtziMduGsVlBVlORxKRAFDZn4WUHtFMzEsLmymYHq+P2fPW8WF5Hf/9TyOYMay305FEJEBU9meppDCLHbVNlNc0Oh3lvHh9lm++8hHvbK7hJ9deyPVjcpyOJCIBpLI/S5d3DUbrzks5Pp/l3+evZ+H6Kr43Yyi3TujvdCQRCTCV/VnqndKD4Tkp3XYpx1rLfy7YxGtrK/nm5YO4Z3Ku05FEJAhU9uegpDCL0j2HqTnS6nSUs2Kt5edvb+GFFbuZNTmXb0wd5HQkEQkSlf05mFaYDcDizd3r7P6375bz+LIdfGXCBTwwvUA3CBdxEZX9ORiclUj/tPhutZTz1Ac7+OXibVw/OocfX3Ohil7EZVT258AYQ0lhFsvL6znaGnoz7ucuq2B5Rd3x9+et3MNPF26mIDuJh68fRkSEil7EbVT252haYTbtXh/LttU6HeVThuekMHteKQ0tHuoa23joLxuIjjQ8NGMoUZH6JxdxI33nn6Mx/XuRmhATkks5xXnp/PamUWyrPkpFbRNREYYnbxvL5MEZTkcTEYeEVNmH8ojjk0VGGKYWZPLulhraO3xOx/mEjfsaeHTxVo7Na7v7koFcOiTT2VAi4qiQKvtQHnF8KiVF2Rxt7WDlztCYcX+4uZ3v/3UDM+d8SHlNIxEG+qTE8eqayk+s4QfbK7Mm8sqsiY4dX0RCrOy7m0sGpdMjOtLxpRyvz/Lyqj184ZGlvLxqL1cUZhEZYRiclUS/1Hjm3DyK2fNKHS18EXGWyv48xEVHcsmgdBaXVWOtMzPuS/cc4kv/83ce/PMGBmUl8cb9FzPygl787pbRpPSIBjrX8OfcPIr1laG/PCYigRHldIDurqQom0Vl1Wzcd4RhOZ3LTzc+vgIgoEsX9Y1tPPz2Fl5dU0lWciyP3TSSa0b0wRjD0N7JADz2zvbjH1+cl05xnubUi7iVyv48TS3IJMJ0DkY7VvaB1OH18dLKPTy6aCvN7V5mTc7l/qmDSIzVP6WInJ4a4jz1Sohh3IBUFm2q5jslQwJ6rNW7DvKD/93E5qojXJyfzn9eU0R+ZmJAjyki4UFl7wclRdn85I0ydtc30T8twe/PX3OklZ+9tYW/lO6jT0ocv79lNFdemK2RByJyxvQCrR+UdM249/euHI/Xx1Mf7OCyR5excH0V91+Wz5LvXMr0Yb1V9CJyVnRm7wf9UuMpyE5i0aZq7r7EP/Phl5fX8YPXN1Fe08gXhmTww5lFDEj3/08NIuIOKns/KSnKZs6726lvbDuv59l/uIX/enMzC9dXcUFqPE9/dSxTh+om4CJyflT2flJSmMVvlmxnyZaac/r7bR1envpgJ3PeLcdnLd+eNph7J+cSFx3p56Qi4kYqez8p6pNM3549WLTp7Nftl26t4UcLythZ18QVRVl8/6pC+qXGByCliLiVyt5PjDFMK8zij6v3cGGfFCLPYGb83oPN/PiNMhaXVZObnsDzd45niiZTikgAqOz9qKQwi+eW76KhxUNqQsxpP67V4+X3SyuYu6yCyAjDA9MLuHPSQGKitDlKRAJDZe9H4wamkhwXxaHm9lOWvbWWxWXV/PiNMioPtTBzRB8emlFA75QeDqQVETdR2ftRdGQEU4dmseDj/Z8ajLazrokfLdjE0q21DM5K5OV7JjAxL82hpCLiNip7P5q7rIL+qfF0+CxHWzsAeG9LNU+8v4O1uw8TGxXBf1xdyG0T+xOt2wOKSBCp7P1oeE4KX39pHQCHmtv51eJt/Pbd7fgsXD86h3+fPoTMpDiHU4qIG+n00o+K89L53S2jMcCBI208tmQ7/XrF89p9E3n0hhEqehFxjM7s/aw4L51e8dEcbPZw6eAMnr593BltwxQRCaSQOrPvTjccP53lFXUcae2gb8841u9rCJn704qIu4VU2Xe3G46fbHlFHbPnlZKfmUhOL937VURCR0iVfXe3vrKBOTeP0r1fRSTkqOz96L4peZ+6z2txXjr3TclzKJGISCeVvYiIC6jsRURcQGUvIuICKnsRERdQ2YuIuIDKXkTEBVT2IiIuoLIXEXEBlb2IiAuo7EVEXEBlLyLiAppnH8ZemTXR6QgiEiJ0Zi8i4gIqexERF1DZi4i4gNbsA0Br5SISanRmLyLiAip7EREXUNmLiLiAyl5ExAVU9iIiLqCyFxFxAZW9iIgLqOxFRFxAZS8i4gLGWut0hk8xxtQCu8/xr6cDdX6Mc66U45OU45NCIUcoZADlONn55Ohvrc041R+EZNmfD2PMGmvtWOVQDuUI/QzKEbwcWsYREXEBlb2IiAuEY9k/4XSALsrxScrxSaGQIxQygHKcLCA5wm7NXkREPi0cz+xFROQkKnsRERcIm7I3xjxjjKkxxmx0OEc/Y8x7xpjNxphNxphvOJQjzhizyhjzcVeOHzmRoytLpDGm1BjzhoMZdhljNhhjPjLGrHEwR09jzGvGmC1dXyNBv62ZMWZI1+fh2K8jxphvBjtHV5ZvdX19bjTGvGyMiXMgwze6jr8p2J+HU/WWMSbVGLPYGLO96/de/jhW2JQ98BxwpdMhgA7gO9baocAE4OvGmEIHcrQBl1lrRwAjgSuNMRMcyAHwDWCzQ8c+0RestSMd3kv9GPC2tbYAGIEDnxdr7dauz8NIYAzQDPwl2DmMMX2BfwXGWmsvBCKBm4Kc4ULgHmA8nf8eVxtjBgUxwnN8urceAJZYawcBS7reP29hU/bW2veBgyGQo8pau67r7aN0fjP3dSCHtdY2dr0b3fUr6K/GG2NygKuAp4J97FBjjEkGJgNPA1hr2621h51NxVSgwlp7rlesn68ooIcxJgqIB/YH+fhDgX9Ya5uttR3AMuBLwTr4aXrriwYwSMIAAAKZSURBVMDzXW8/D1zrj2OFTdmHImPMAGAUsNKh40caYz4CaoDF1loncvwa+DfA58CxT2SBRcaYtcaYex3KkAvUAs92LWs9ZYxJcCjLMTcBLztxYGvtPuARYA9QBTRYaxcFOcZGYLIxJs0YEw/MAPoFOcPJsqy1VdB58ghk+uNJVfYBYoxJBOYD37TWHnEig7XW2/Wjeg4wvutH1qAxxlwN1Fhr1wbzuKcxyVo7GphO59LaZAcyRAGjgd9ba0cBTfjpR/RzYYyJAa4B/uTQ8XvReRY7EOgDJBhjvhLMDNbazcDDwGLgbeBjOpdiw47KPgCMMdF0Fv1L1to/O52na6lgKcF/TWMScI0xZhfwR+AyY8yLQc4AgLV2f9fvNXSuT493IEYlUHnCT1iv0Vn+TpkOrLPWVjt0/MuBndbaWmutB/gzUBzsENbap621o621k+lcUtke7AwnqTbG9Abo+r3GH0+qsvczY4yhc012s7X2lw7myDDG9Ox6uwed31hbgpnBWvugtTbHWjuAzuWCd621QT1zAzDGJBhjko69DZTQ+eN7UFlrDwB7jTFDuh6aCpQFO8cJ/hmHlnC67AEmGGPiu75vpuLAC9bGmMyu3y8ArsPZzwnA68BXu97+KvC//njSKH88SSgwxrwMXAqkG2MqgR9aa592IMok4FZgQ9d6OcBD1to3g5yjN/C8MSaSzv/UX7XWOrb10WFZwF86+4QoYJ619m2HstwPvNS1hLIDuMOJEF3r09OAWU4cH8Bau9IY8xqwjs6lk1KcGVkw3xiTBniAr1trDwXrwKfqLeDnwKvGmLvo/A/xy345lsYliIiEPy3jiIi4gMpeRMQFVPYiIi6gshcRcQGVvYiIC6jsRURcQGUvIuIC/x9VmTPRC3312QAAAABJRU5ErkJggg==\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": {
154
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD4CAYAAAAQP7oXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3xV5Zno8d+ThJsockcEBNSARUEFimi1tlo70HoGe2bao8cpTMdPGat2pjOfmSlO65zpmU6P7cyZtp46UnWsaLWWXkkVReWioAYItyByC+GShFu4BcItyV7P+WO92Xtlk+y9cl3Zez/fz2d/9nrfvd61nr1IeLLetdb7iqpijDHGhJEXdQDGGGMyhyUNY4wxoVnSMMYYE5olDWOMMaFZ0jDGGBNaQdQBdLbBgwfrmDFjog7DGGMyyrp1646o6pDk+qxPGmPGjKGkpCTqMIwxJqOIyN7m6q17yhhjTGiWNIwxxoRmScMYY0xoljSMMcaEZknDGGNMaJY0jDHGhGZJwxhjTGiWNIwxxoRmScMYY0xoWf9EuDGZYMy81+LLex7/fISRGJOanWkYY4wJzc40jOlCdkZhMp2daRjTxfpwjlFyCI7vgfqzUYdjTKuEOtMQkRnAj4F84FlVfTzpc3Gffw44A/y5qq5P1VZEBgK/BMYAe4AvqerxwDavAD4C/llV/93VTQGeB/oAi4G/VlVtw/c2psuMmfcaI+Uwn81bxzM9PmKYHPM/WLgARGBQIXfnDeNNbyp19Ig2WGPSSHumISL5wJPATGACcJ+ITEhabSZQ6F5zgadCtJ0HLFXVQmCpKwf9EHg9qe4pt/3Gfc1I/xWNiVBNFX9f8ApP9niCe/JXJRJGI1U4soO5Ba/yZI8fc63sjiZOY0IK0z01DShT1XJVrQNeAWYlrTMLeEF9xUB/ERmepu0sYIFbXgDc07gxEbkHKAe2BOqGA/1U9QN3dvFCsI0x3YoXg02/hF9/hdvyNiMkTohj5LNfB0O/ESCJX8Fhcpx/7fEcbE/+W8mY7iNM99QIoCJQrgRuCrHOiDRth6nqAQBVPSAiQwFEpC/wTeAu4O+S9lHZzD4uICJz8c9IuOKKK1J/O2PaodkL2ycPwNLvwOGtTdYt8cbzWuwmNulVNFAAm6AvZ7kzbwP35i/jYjlLHh68830o6A1XfdounJtuJ0zSkGbqkq8jtLROmLbJvgP8UFVr/UslrYrDr1R9GngaYOrUqXbNw3SdfcWw7Ltw/lS8arcO58mGWezQUResfpo+FHm38L53LY/1eJGxcsDvslr+r3DpyK6M3JhQwiSNSiD40z4S2B9ynZ4p2h4SkeHuLGM4cNjV3wT8qYj8AOgPeCJyDviNa58qDmOis+X38N6P/P/0AfIKYMqf87e7LyVGfsqmR7iUb9f/BT/o8VMmUQ+xelj2L/Tkv9vFcdOthLmmsRYoFJGxItITuBcoSlqnCJgtvulAjet6StW2CJjjlucAiwBU9TZVHaOqY4AfAd9T1Z+47Z0Skenubq3ZjW2MidqsvFWw6ofxhLG8UplR/qeMWTgwbcJodIqL+Jf6L/tdUwDH9/Kl/BWdFLExbZM2aahqA/AIsATYCixU1S0i8qCIPOhWW4x/4boMeAZ4KFVb1+Zx4C4R2Yl//aLJbbwt+BrwrNvPLi68u8qYLndr3mYeKHid0qoaSqtq+E3lJXyj/mG2aeuvp+1nMNz8cLz8hfxVDOV4ihbGdK1Qz2mo6mL8xBCsmx9YVuDh5HYttXX1R4E70+z3n5PKJcB1YWI2piuMlQP8TcGv4+Ut3hi+0zCbc/Rq+0avuRu2vQbV2+hBA7ML3uTfG/5HB0RrTPvZE+HGtFWsnr8t+BU9aABgvw7mXxvub1/CAMjLg1u+Hi/emreZyznSvm0a00EsaRjTVusXMFoOAVBHD/6l/s+o5aKO2fZl18GoaQDkofxp/rsds11j2smShjFtUVMFG38RL77Q8FmqGNKx+7jxy/HFT+dvhDPHUqxsTNewpGFMW6x9Fjy/W2qbdwV/8KZ3/D6GT2K759+xnk/MnhQ33YIlDWNa6/A22LUsXvyv2Ey0k36VFnuBwRe2FoHndcp+jAnLkoYxrbXhxfji+951bG/DrbVhrfImUqt9/MKpg1CxutP2ZUwYNgmTMa1xfC/sWRUv/jyW8q7xNgmONwUFvOVN4Qv5bp873oDRN3f4Po0Jy840jGmN0l8mlkd/gkod2um7fNubnCjsfb/JuFbGdDVLGsaEdfoI7HwzUb7hvi7ZbYUOo1yH+4VYHZS/0yX7NaY5ljSMCWvrHyBWT2lVDQsrLmXMj/Z12a6XxW6MD1Py8kvPdtl+jUlmScOYEK6c9weWL/oZpVU1ABR5t3Tp/ld6E/Hc7ADX5e3xL4obEwFLGsaEMFW2M0hOAnBCL6bY+1iX7v84/djoXZ2o2PlWl+7fmEaWNIwJYUb+2vjyUm+yP/NeF1vh3ZAolL2dmLfDmC5kScOYdGqrmZy3I158MzY1kjCKvY9xvnFCpuN74OiuSOIwuc2ShjHplL1NnptZuNS7kgMMiiSMc/Si2JsQiMu6qEzXs6RhTDqB/5yXeTdGGAisiF2fKJQttWFFTJcLlTREZIaIbBeRMhGZ18znIiJPuM9LRWRyurYiMlBE3hKRne59gKufJiIb3WuTiHwh0GaF21bj553/ZJXJbUd3xbuB6ujBB8G/9COwUa/mpLrh109Xw4GNkcZjck/apCEi+cCTwExgAnCfiCT/5swECt1rLvBUiLbzgKWqWggsdWWAD4GpqnoDMAP4qYgErzrer6o3uNfh1n5hY1ql7O344mrvGs7SO8JgIEY+K71JiYpAfMZ0hTBnGtOAMlUtV9U64BVgVtI6s4AX1FcM9BeR4WnazgIWuOUFwD0AqnrGzS0O0BuwW0RMNDyvya2tK2I3pFi567zjBbqoyt+BhrrogjE5J0zSGAFUBMqVri7MOqnaDlPVAwDuPd7VJCI3icgWYDPwYCCJAPzMdU09JiLSXMAiMldESkSkpLq6OsRXNKYZBzf5XUDAKb2I9VoYcUC+bToK+l3uF+pque+ffsSYea8lDXRoTOcIkzSa+485+a//ltYJ0/bCFVRXq+q1wMeBR0WksU/gflWdCNzmXl9uof3TqjpVVacOGdLBs6mZ3LEz0fWz0ptIjPwIgwkSuDoxuu6n8uy6huk6YZJGJTAqUB4J7A+5Tqq2h1wXFu79gusTqroVOA1c58pV7v0U8DJ+95cxHa+hDnYnBgZs0iXUHVx9V3zx43nbuZgzEQZjckmYpLEWKBSRsSLSE7gXKEpapwiY7e6img7UuC6nVG2LgDlueQ6wCMCtW+CWRwPjgT0iUiAig119D+Bu/IvmxnS8yjWJIcgvGc7WTpxoqU0GjIbB4wAoIMYteVsiDsjkirRjIahqg4g8AiwB8oHnVHWLiDzoPp8PLAY+B5QBZ4CvpGrrNv04sFBEHgD2AV909bcC80SkHvCAh1T1iIj0BZa4hJEPvA080+4jYExzyt6OD064MHYjzfe0RqzwLjjiP6n+qbxNvOl9POKATC4INYCOqi7GTwzBuvmBZQUeDtvW1R8FLpj2TFVfBF5spv40MCVMvMa0S90Z2PNevPhO8BbX7uSqO6D4PwG4Nm8Pg6iJOCCTC+yJcGOS7VnlT3YE7NHLqNBhEQfUgr6D4XL/OVpBuT1vU8QBmVxgScOYZIEH5t6JdbML4Mmu/kx88fb80ggDMbnCkoYxQWeOQWViGPR3vYkRBhPC2E9S73qZx8oBf/RbYzqRJQ1jgspXgPqDAH7kjaaaAdHGk06vi1nrjU+Uy5ZGF4vJCV0/k4wx3VngP91u92xGQPDp71vyJiVuud21DKb+BTQ/WIIx7WZnGsY0OnkADrlHf/Lyec+7Ltp4QirxxnOWXn6hphKqt0cbkMlqljSMabQr0LUz8uOcpG90sbRCHT2SJmeykW9N57GkYUyj4H+2gbuSMsE7scCzJOXLbXIm02ksaRgD/kRLx3b7ywW9YPQnoo2nlTbpVdSoOzM6fcQfodeYTmBJwxiIXwAvrarhyd3DGPNPyyMOqHViJF2D2ZVZ8ZvMYUnDGM9rcj2j2z/Q14JVwaSxZ5V1UZlOYUnDmMNb4NRBAGq1T7eZbKm1PtIx0LufXzhzFKq3RhqPyU6WNIwJXABf5V1HQ4Y+vuSRB6NvTVTsXhldMCZrWdIwuS3W0KT/vzs/0BfGF5cUUFpV4w/rvmclaNqJMo1pFUsaJrdVrYNzbkjxvkP4SEdHG087bdSrmz7od3x3tAGZrGNJw+S2Js9m3Ilm+K9EPQWs88YlKqyLynSwUL8hIjJDRLaLSJmIzGvmcxGRJ9znpSIyOV1bERkoIm+JyE73PsDVTxORje61SUS+EGgzRUQ2u209IWID7Jh2qD8He1bGu3M+vahH1BF1iA+CT4fvfje6QExWSps0RCQfeBKYCUwA7hORCUmrzQQK3Wsu8FSItvOApapaCCx1ZfDn/Z6qqjcAM4CfNs4Z7rY7N7CvGa39wsbE7X0P6s8CUKFD2K2XRRxQxyjxxtFAvl84WuaPqWVMBwlzpjENKFPVclWtA14BZiWtMwt4QX3FQH8RGZ6m7SxggVteANwDoKpnVLXB1fcGFMBtr5+qfuCml32hsY0xbXLBiLbZceJ6lt5s8q5KVOx7P7pgTNYJkzRGABWBcqWrC7NOqrbDVPUAgHsf2riSiNwkIluAzcCDLomMcO1TxdHYfq6IlIhISXV1dYivaHLOuZNQsTpefDfWTecBb6PV3scShb2WNEzHCZM0mvvzK/k+vpbWCdP2whVUV6vqtcDHgUdFpHdrtqWqT6vqVFWdOmTIkHS7M7lo9zvg+Se0O3QUBxkUcUAdqyQ4MdP+jVB3OrpgTFYJkzQqgVGB8khgf8h1UrU95LqcGrueDifvWFW3AqeB69y2RqaJw5hwmswDnl1nGQBHuBQGXe0XvAaoLIk2IJM1wiSNtUChiIwVkZ7AvUBR0jpFwGx3F9V0oMZ1OaVqWwTMcctzgEUAbt0CtzwaGA/scds7JSLT3V1TsxvbGNMqtdVwwI0CK3ms7O7zgLfRtzdcknjQz7qoTAdJO16CqjaIyCPAEiAfeE5Vt4jIg+7z+cBi4HNAGXAG+Eqqtm7TjwMLReQBYB/wRVd/KzBPROoBD3hIVY+4z74GPA/0AV53L2NaZ/c7lFaeAKDUu5ITXBJxQJ1jtfcxvpS/wi9UFPsDGOZl9nMoJnqhBtlR1cX4iSFYNz+wrMDDYdu6+qPAnc3Uvwi82MK2SvC7qoxpu/LEsCHZepYBUKaXc0Ivpr/UwtkTcPgjuMx+fUz72J8dJrfUHoaD/jzgHtL0Qbgso+Q1vSC+74PogjFZw5KGyS3l78QXS72rOMnFEQbT+dboNYmCXdcwHcCShsktOdI11Wijd1Xi6fBj5fF5Q4xpK0saJnfUHoZD/n0YHnkUZ3HXVKNz9Gr6dLidbZh2sqRhcsfe9+KLpd6VnOKiCIPpOnZdw3QkSxomd+xN/IfZZJiNLLemydPhG6DuTHTBmIyXmfNaGtNadWdYX7yUAmJA0n+kWa6aATDwSv+aRqweqkpg7CejDstkKDvTMLmhqiSeMPbqMP8/0lwy+pbE8l7rojJtZ0nD5IYc7ZpqNON3DfEhRd5d+gf/6XBj2sCShsl+ntfkAvDaHOqaarRdR1KjfQH8J8SPbI84IpOpLGmY7Fe9Fc4eB+CEXswOHZmmQfZR8lingbnD7dZb00aWNEz2q1gTX1znjUNz9Md+jRd4OtxuvTVtlJu/PSa3VCXmktigV0cYSLQ2eFfjNf7KH9npDxFvTCtZ0jDZre40HPooXtwYfDo6x5ylN5u9sYkKO9swbWBJw2S3/RtA/TuFynV41g9QmE6TmwD2FUcXiMlYljRMdgtMc7rRy92uqUZrg9c1qkqg4Xx0wZiMFCppiMgMEdkuImUiMq+Zz0VEnnCfl4rI5HRtRWSgiLwlIjvd+wBXf5eIrBORze79jkCbFW5bG91raPu+vsl6lWvjixssaXCAQVTpYL/QcN4/EzOmFdImDRHJB54EZgITgPtEJHl40JlAoXvNBZ4K0XYesFRVC4GlrgxwBPhvqjoRf+7w5Fn87lfVG9zrcGu+rMkxpw5CTaW/nN+Tj3RMpOF0F03ONuzWW9NKYc40pgFlqlquqnXAK8CspHVmAS+orxjoLyLD07SdBSxwywuAewBUdYOq7nf1W4DeItKrjd/P5LJA1xTDr6fehloD/HG3Gp8OX/bmIlCNOiSTQcIkjRFARaBc6erCrJOq7TBVPQDg3pvravoTYIOqBjtef+a6ph4TEWkuYBGZKyIlIlJSXW23FeasqpL4f47feK9H1NF0G1t1NGfoDcBgqfEHMjQmpDBJo7n/mJP/NGlpnTBtm9+pyLXA94G/DFTf77qtbnOvLzfXVlWfVtWpqjp1yJAhYXZnso3nNb0Irrl7q22yGPms8woTFdZFZVohTNKoBEYFyiOB/SHXSdX2kOvCwr3Hr0+IyEjgd8BsVd3VWK+qVe79FPAyfveXMRc6sgPOnwL8oUP26GURB9S92K23pq3CJI21QKGIjBWRnsC9QFHSOkXAbHcX1XSgxnU5pWpbhH+hG/e+CEBE+gOvAY+qanyqNREpEJHBbrkHcDfwYau/sckNVevii/5ZRrM9mTlrvTcOr/GYHN4SH5vLmHTSJg1VbQAeAZYAW4GFqrpFRB4UkQfdaouBcqAMeAZ4KFVb1+Zx4C4R2Qnc5cq49a8GHku6tbYXsERESoGNQJXblzEXCiYNu9X2Aifpy3bvCr+g2mR8LmNSCXU7iaouxk8Mwbr5gWUFHg7b1tUfBe5spv67wHdbCGVKmHhNjmuog4Ob48VNOTx0SCprdTwfY69f2Ps+jPujaAMyGcGeCDfZ59BmiNUBUKWDOcqlEQfUPTUZ9bZyLcQaogvGZAxLGib7VK2PL9pZRsv26VAOa3+/UHcaDpZGG5DJCJY0TPYJXM/YZLfapiBNnw63UW9NCJY0THY5fwqq3VSmktd0KHBzgSa33trzGiYESxomu+zfGB8KncHjqOWiaOPp5jbrlVDgPx1OTSWcqEjdwOQ8SxomuwS6phhhN9ulU08BjJyaqLAH/UwaljRMVnl9yavx8aZm/f5s1OFkhituTizvsy4qk5oN+2myR201I+QI4P8FvVVHRxxQZpi84BTP96wBYFLeJv+6UK9LIo7KdFd2pmGyR6Br6iNvNHXYyLZhHKMfO9UffLq04hhf+c6PGTPvtYijMt2VJQ2TPQJJo1SvjDCQzPN+7Lr48m15m1OsaXKdJQ2THVSbJA2b2rV1VnmJpHFjXhl9setBpnmWNEx2OLEXzhwF4LT2YZdeHnFAmeUQAylzXVT5xJiWty3iiEx3ZUnDZIeKtfHFzToWtR/tVnsv0EV1a57NOmCaZ79ZJjtUJob2bjIrnQkt2EU1OW9nfBIrY4IsaZjMV3/OfxLcWe+NizCYzHWIgfFuvXxiNqyIaZYlDZP5DmyMD4W+T4dSTf+IA8pcq2ITE4XyFZHFYbqvUElDRGaIyHYRKRORec18LiLyhPu8VEQmp2srIgNF5C0R2eneB7j6u0RknYhsdu93BNpMcfVlbn82h6fhez99Pv4U+AbrmmqXYBcVlWuti8pcIG3SEJF84ElgJjABuE9EJiStNhModK+5wFMh2s4DlqpqIbDUlQGOAP9NVSfizx3+YmA/T7ntN+5rRmu+rMlOU/J2xpdLrGuqXYJdVMTqYa8Nl26aCnOmMQ0oU9VyVa0DXgFmJa0zC3hBfcVAfxEZnqbtLGCBW14A3AOgqhtUdb+r3wL0FpFebnv9VPUDN73sC41tTA6rqWS4+LfanqcHH+mYaOPJAu8FzzbKl0cXiOmWwiSNEUBwvORKVxdmnVRth6nqAQD3PrSZff8JsEFVz7t2lWniAEBE5opIiYiUVFdXp/hqJuNVJO6aKvWu8kdtNe2yKpbcRVUbXTCm2wmTNJq7bqAh1wnTtvmdilwLfB/4y1bE4VeqPq2qU1V16pAhQ8LszmSqitXxRbvVtmMcZFBSF5XdRWUSwiSNSmBUoDwS2B9ynVRtD7kuJ9z74caVRGQk8DtgtqruCuxjZJo4TC5pOA/7N8SL6+x6Rodp2kW1IrI4TPcTJmmsBQpFZKyI9ATuBYqS1ikCZru7qKYDNa7LKVXbIvwL3bj3RQAi0h94DXhUVd9r3IHb3ikRme7umprd2MbkqMoSP3EAVTqYQwyMOKDs0bSLao11UZm4tElDVRuAR4AlwFZgoapuEZEHReRBt9pioBwoA54BHkrV1rV5HLhLRHYCd7kybv2rgcdEZKN7NV7v+BrwrNvPLuD1Nn9zk/n2rIovFnvJN/SZ9jjIIBjsuvusi8oEhLpqqKqL8RNDsG5+YFmBh8O2dfVHgTubqf8u8N0WtlUCXNfcZybHeF6TWeaKvY9FGEx2+vv1A/hyvpucqXwFjPtstAGZbsGeCDeZ6dCHcPYEACf0YnboyDQNTGs1eTrcuqiMY0nDZKYmXVMfs1FtO8EBBlGuw/2CdVEZx37TTOZRtesZXcTuojLJLGmYzHN8N5ys8pd7XMRmm9q10zTtolrrjyhscpolDZN59ryXWL7iJnsKvBMdYBAV6h6QjdU1mVLX5CZLGibz7H43sTzm1ujiyBFrvWviowj/2zPPRx2OiZglDZNZairhyA5Kq2pYX1XLhGdroo4o660O3M48LW+bf7uzyVmWNExm2ZUYdXW9V8gZekcYTG7YpqM4qRcBMEBOwZEdEUdkomRJw2SWXcviiyu9SREGkjuUPNbp+ETF3vdaXtlkPUsaJnMc3wPHygGop4DV3jXRxpNDmhxre14jp1nSMJkj0DW11hvPOXpFGExuWe8V0kC+XzhaBrWHUzcwWcuShskMqixe9HL8Lp6V3sT0bUyHOUcvNntjExXWRZWzLGmYzHCsnJHiz8J4jp6UeOPTNDAdbU1wUEibOzxnWdIwmSFwAXytN57z9IwwmNy0Npio92+A+rPRBWMiY0nDdH+qTa5n2F1T0TjMAPapm9omVgdV66MNyETCkobp/qq3x8eaOksvm9Y1QmuDd1Htsy6qXBQqaYjIDBHZLiJlIjKvmc9FRJ5wn5eKyOR0bUVkoIi8JSI73fsAVz9IRJaLSK2I/CRpPyvctpJn9DPZbOeb8cVib4KNNRWhEm98/GaE5W8V+WeBJqekTRoikg88CcwEJgD3iUjyWNQzgUL3mgs8FaLtPGCpqhYCS10Z4BzwGPB3LYR0v6re4F5231+2izXArqXx4rLYDREGY7bqFdRqHwAGyUk4uiviiExXC3OmMQ0oU9VyVa0DXgFmJa0zC3hBfcVAfxEZnqbtLGCBW14A3AOgqqdVdRV+8jA57k++/QSlZXsprarhmPaj1IZBj5RHHus00D24zx70yzVhksYIoCJQrnR1YdZJ1XaYqh4AcO9hu5p+5rqmHhMRaW4FEZkrIiUiUlJdXR1ys6Y7uiN/Q3z5HW+SzdDXDTS53XlfcXSBmEiE+Q1s7j/m5I7MltYJ07Y17lfVicBt7vXl5lZS1adVdaqqTh0yZEg7dmciVXea6Xlb48XlnnVNdQfrvEK8xl/twx/BmWPRBmS6VJikUQmMCpRHAvtDrpOq7SHXhYV7T3t9QlWr3Psp4GX87i+TrXa/Sw8aANirw9jTOF+1iVQtF7HVG+0XVKFiTbQBmS4VJmmsBQpFZKyI9ATuBYqS1ikCZru7qKYDNa7LKVXbImCOW54DLEoVhIgUiMhgt9wDuBv4MET8JlMF7ppabhfAu5WmXVR2620uSXvvoqo2iMgjwBIgH3hOVbeIyIPu8/nAYuBzQBlwBvhKqrZu048DC0XkAWAf8MXGfYrIHqAf0FNE7gE+C+wFlriEkQ+8DTzTvq9vuq3aw/5Tx4AivONdH3FAJmiNjmcOS/xC5VqI1UN+j2iDMl0i1A3vqroYPzEE6+YHlhV4OGxbV38UuLOFNmNaCGVKmHhNFti+OP4MQKl3JUe5NOKATFCFDuWQDgA8qDsNB0thhP165gK7FcV0P54H2xJ/Z7zpTY0wGNM8adpFZQMY5gxLGqb7qVwLtYcAOKUXUewlP0tquoO1dl0jJ1nSMN3Ptlfji8u8G2zYkG5qs14JBW6O9ppKOFGRuoHJCpY0TPdy5liTCX6WeB+PMBiTSj0FPLt7QHwsKnvQLzdY0jDdy9Yi8GL+8mXXUak2JmV31mTU2z0rowvEdBlLGqb7aKiDLb+P/+U6e3XyaDWmu1ntXZN4OvxgKZw+Em1AptNZ0jDdx66lcPY4AEe1H+9510UckEnnBJewpXHucFUoXxFpPKbzWdIw3YMqbP51vPiaN50Y+REGZMJ615uYKARmWDTZyZKG6R72b4CjZQCcpwdvxOwCeKb4wLs20UV16EM4dTDagEynsqRhuofShfHFZbHJ1HJRhMGY1jhJXzZ6Vycq7Gwjq1nSMNE7spPSD96gtKqGTVUn+YN3c9QRmVZaGeyiKns7ukBMp7OkYaK3/oX44vvetVSqzYGSaYq9CZDf0y8cLYPqHdEGZDqNJQ0TrWO7Yfe78eLC2Keii8W02Wn6wJW3Jyq2vxZdMKZTWdIw0drwYnxxjXcNu22ipYx191uXJp4O3/k2NJyPOiTTCSxpmOicqGhy0fSV2KcjDMa014c6lgM6yC/U1UL5O9EGZDqFJQ0TnY0vgXoArPcKKdOREQdk2kd4ywvMqREYeNJkj1BJQ0RmiMh2ESkTkXnNfC4i8oT7vFREJqdrKyIDReQtEdnp3ge4+kEislxEakXkJ0n7mSIim922nhARaftXN5GqqYIdS+LFX9pZRlZYGrsx8czGgU1wdFe0AZkOlzZpiEg+8CQwE5gA3CciyRMczAQK3Wsu8FSItvOApapaCCx1ZYBzwGPA3zUTzlNu+437mhHqW5ruZ93z8bMMLr+RrTo60nBMxzhOP94PDv8SeMrfZIcwZxrTgDJVLVfVOuAVYFbSOrOAF9RXDPQXkeFp2s4CFrjlBcA9AKp6WlVX4SePOLe9fqr6gfeQB6MAABL4SURBVJte9oXGNibDHNsNZW/FL5rOeH98+jYmYxTFbkkUyt72h7s3WSNM0hgBBGdXqXR1YdZJ1XaYqh4AcO/pxsAe4dqnigMAEZkrIiUiUlJdXZ1ms6bLlTwXn/97nTeObXpFxAGZjrRNR/Gbykv8Pwr2VfvD3ZusESZpNHfdQEOuE6ZtWKG3papPq+pUVZ06ZIg9KNatVG9v8lzGz2OfiTAY0zmERcGzjS2/94e9N1khTNKoBEYFyiOB/SHXSdX2kOtyaux6OhwijuDtNc3FYbq7tf8VX3zfu5ZdanNmZKP3vOs4qv38wtnjNrRIFgmTNNYChSIyVkR6AvcCyeebRcBsdxfVdKDGdTmlalsEzHHLc4BFqYJw2zslItPdXVOz07Ux3czBzZQWvxkfY+ql2J1RR2Q6SYx8Xo0FxhDb+BJ4XnQBmQ6TNmmoagPwCLAE2AosVNUtIvKgiDzoVlsMlANlwDPAQ6naujaPA3eJyE7gLlcGQET2AP8B/LmIVAbuuPoa8Kzbzy7g9TZ+b9PVVGHNM/HiCu96KnRYhAGZzva6N40z9PYLNZWwe0Wk8ZiOURBmJVVdjJ8YgnXzA8sKPBy2ras/CjT7p6aqjmmhvgSw6dwyUdV6/759wCOPXzTcEXFAprOdoTevxqYznQ1+xYaX4MpPgz1eldHsiXDT+VRhbeIs483YFA4yKMKATFcpit0CBe5s42gZ7Psg2oBMu1nSMJ3uS//4fyndUExpVQ0N5LPQnv7OGSfpy7+UX5UYyHDDz+O3W5vMZEnDdC4vxuyCN+PFxbGbOMKlEQZkutrvYrfR0Djf+6EtsH99tAGZdrGkYTrX9te5Qvy7qc/Si1/Fbk/TwGSbY/Tj7djkREVg0i2TeSxpmM5Tf9Z/+tv5Tew2arg4woBMVH4dux2v8b+b/Ruhal20AZk2s6RhOs/mX8GZowAc10tYFPtExAGZqBxmQNOzjbXP2bWNDGVJw3SOM8dg4y/ixZ/HPsN5ekYYkInaL2OfYkNVrX9RfP17ULEm6pBMG1jSMJ2j5DmoPwPAPh3KUm9ymgYm21UzgDdiH09UlNjZRiaypGE63uFtTWZte77hjxL92San/Sp2O/WNzxRXb4O970cbkGk1+002HcvzYNUPE39BXnEzJWrzZRjfMfqxODYtUVHynI1JlWEsaZiOtf01/y9IgPyecMvXaX5Ue5Orfh27nfP08AtHy2DPu6kbmG7FkobpOGeOweqfAlBaVcM/7p7ImP+zMeKgTHdTw8VNR8Bd8wzE6qMLyLSKJQ3TMVRh1X/A+VMAHNIB/Dr2yYiDMt3Vb2O3QU/3zE5NJXxksxxkCksapmPsWga7V8aLP2m4h7rGLghjkpziIpg8O1Gx7nk4dzKyeEx4ljRM+505Bu/9KD4o3Q/2jWOTXh11VKabK/xFb96qzPcHMjx/yoYXyRCWNEz7qMLK/xv/K7Fa+/Oz2IyIgzKZoJ4Cngv+rGz5LZyoiC4gE0qopCEiM0Rku4iUici8Zj4XEXnCfV4qIpPTtRWRgSLylojsdO8DAp896tbfLiJ/FKhf4eo2utfQtn910yG2/Bb2rIoX/1/DFzjbOFubMWl84E3gI2+0X/BiTW/XNt1S2qQhIvnAk8BMYAJwX2D61UYzgUL3mgs8FaLtPGCpqhYCS10Z9/m9wLXADOA/3XYa3a+qN7jX4dZ/ZdNhqnewfuH34t1Sf4jdzEbrljKtIjwT+zwbq076P0drlsHON9M3M5EJc6YxDShT1XJVrQNeAWYlrTMLeEF9xUB/ERmepu0sYIFbXgDcE6h/RVXPq+pu/PnAA08DmW6h7gws/Q4FxADYpZfzvHVLmTbYpSP4Q+yWRMUHP4Gzx6MLyKQUJmmMAIIdjZWuLsw6qdoOU9UDAO69sasp3f5+5rqmHhNpfrJhEZkrIiUiUlJdXZ3u+5nW8jxY/q/+rZL482T8W/3/SAwPYUwrvRS7k2rt7xfOnYT3fxJtQKZFYZJGc/8xJ3c6trROmLat2d/9qjoRuM29vtzcBlT1aVWdqqpThwwZkmZ3ptXW/azJdYz/bJjFfgZHGJDJdOfoxX82/HGiouxt/2W6nTBJoxIYFSiPBPaHXCdV20OuCwv33nh9osU2qlrl3k8BL2PdVl1v17Imt0Yuit3CO971EQZkssU6Hc9y78ZExcr/gJMHogvINCtM0lgLFIrIWBHpiX+RuihpnSJgtruLajpQ47qcUrUtAua45TnAokD9vSLSS0TG4l9cXyMiBSIyGEBEegB3Ax+24TubtqpcB8u/lyiPmsbPYjOji8dknfkNd0O/y/1C3WlY9l1oqIs2KNNE2qShqg3AI8ASYCuwUFW3iMiDIvKgW20xUI5/0foZ4KFUbV2bx4G7RGQncJcr4z5fCHwEvAE8rKoxoBewRERKgY1AlduX6QoHP4Ql/wixekqranijsgfXvnuTDXluOtRZevPZ0tvjd1Nx6EO7DbebEc3yf4ypU6dqSUlJ1GFktqp1sOTb8UmVllXCN+u/SjUD0jQ0pm2+kLeSrxS8waQRl/oVNz8Ck74YbVA5RkTWqerU5Hq73cWkVrYUln+P0gp/ru8a7ctjljBMJ/uddytjvENMotyvKH4S+gyAws9EG5ixvgXTAi8GxfNh6f8GrwGAo9qPbzU8QBV2R5rpbMJPGu6BYdf5RVX/Nu/yd6INy1jSMM2oPQyvfoPSxfPjT3tX6hD+of4v2afDoo7O5Ih6CmDG92DgWL9CPVj6Hdj+RrSB5ThLGiZBFba+CgvnwIHSePU6bxz/UD+XavpHGJzJRWP+eRXXr5vJ65VumH0vBiv+D6xbYBfHI2LXNIzv1CF49wdQmbhpwEN4KfYZfh37JGp/X5iI1HAx36p/gPrKBYyRgwBM4jk4sgNu/yb07hdxhLnFkkauU4Wtf4Dip6D+jH+bI7BfB/Pjhq+yVUdHHKAxcIx+fLP+q/xjwctcn7fLr9yzyp9j/I5vw2UTow0wh1jSyGXH9/pzYRzYFK9ShEWxT/BS7E7O0zPC4Ixp6iy9+U7DbObkL2ES7nGvUweh6Osw4R6YNhd6XhRtkDnAkkYuajgPG16Ejb8AryF+dlGlg/lxw1y26RURB2hM8xoo4L9in2fL3rF8o+A3XMQ5/1mOLb+Dve/DzQ/B2Nuh+bFMTQewpJFrKtb6T9ierIpXeQi/j93Ky7E7bV5vkxGKvQl8ve5yHipYBFU7XG0Nk2r/Fwy/Hm75Kxhsc7t0BnsiPFecqIDV85uMTgvAsOu4ffVk9upl0cRlTLson8rbxFfzX+MSOZN4glwExs2EybOh3/BoQ8xQ9kR4rjpzDDb8HD76PaUVx+LVp7UPz8c+y5u7p9qdUSaDCSu8GyjxxnFv/nK8qmLy8ACYpIv9WQDHu+Rxsc0O3REsaWSrmkooXQjbX4dY01FCl3s38FzDTGq4OKLgjOlYtVzEs7HP84Y3jQfyFzMlz3VZeQ3+3YHbX4crPwXXfgGGXWvXPNrBuqeySU0VVK7x57xwD+c1XuQG+NAbw3OxmZTpyKgiNKZLTJRy7s9/mwl5ewES3VbgP2E+9nYY+0kYeKUlkBZY91S28GJwcj/UVPjXKWoqEstnjjbbZKeOYGHDp1mt19D8xIjGZJfNeiXzGr7K9bKL/5m/jEmcSHx4bLf/Wvc89B0Cl9+YeNn1j7TsTKO7UoXaQ/4P93H3Q36sHE7sg1hdkzMIaPqXVGlVDR7COm88RbGb2aRXYcnC5LKrpIrP563m9vxNTBnRN15/we/RuKth2AQYeq3/PngcFPTq6nC7hZbONCxpRC1W7w8QWFMJx/e4127/wTs3f0Wj5B/woEkjLoUefWD4DTyysoCV3kS7ZmFMkt6cZ0reDm7J28LUvB304XyTz5t0YwHk5cOgQhj6MT+BDLoKBoyFgux/8LVdSUNEZgA/BvKBZ1X18aTPxX3+OeAM8Oequj5VWxEZCPwSGAPsAb6kqsfdZ48CDwAx4K9UdYmrnwI8D/TBny3wrzXNF+iSpOF5/gW3xlf9WTh/CupO+e/na6GuFs6f9JfPn/QTxamDcOZIyoHXmksUx7QfFTqE/TqYKh1EpQ5hvw7iMANsJj1jQsrD40rZzyQpZ2Lebq7N20NvEjeNXJBAGkke9B8FA6+CSy7zu7j6DoYeff0/3Hr08c9OJN9POpIHeQXulR+o795n/21OGiKSD+zAn5K1En/e7/tU9aPAOp8Dvo6fNG4CfqyqN6VqKyI/AI6p6uMiMg8YoKrfFJEJwC+AacDlwNvAOFWNicga4K+BYvyk8YSqvp4q/jYljY0v++P2BxOBF2u5rF7rtu8EE0Jy9xJArfZhjw5jr17GXh3GPh3KPh1KLTZUgjEdLQ+P0XKIcVLBNVLB+LwKRkr1Beu1mExaS/KaJpH4cjC55AXqCgJJyL2n03803PqNtoXXjgvh04AyVS13G3oFmIU/h3ejWcAL7q/+YhHpLyLD8c8iWmo7C/iUa78AWAF809W/oqrngd0iUgZME5E9QD9V/cBt6wXgHiBl0miTUwegelu7NpGqKylIEY7pJbxS0Z99OpQKlxj26GXU0Be7FmFM1/DIY7cOZ7cOZwnTIAYXc4ZxUklhXhVj5SBj5CBadRShbd36TRKOehDzgPqO+QLNqT/X4ZsMkzRGABWBciX+2US6dUakaTtMVQ8AqOoBEWl88mYE/plE8rbq3XJy/QVEZC4w1xVrRWR7S1+unQYDRzpp29nEjlN6dozSi+QYbenqHbZP0jF6F5jf1m01O8R1mKTR3J+6yWm2pXXCtA27v9DbUtWngafT7KfdRKSkudM305Qdp/TsGKVnxyi9rjhGYa6aVgKjAuWRwP6Q66Rqe8h1YeHeD4fY1shm6o0xxnSRMEljLVAoImNFpCdwL1CUtE4RMFt804Ea1/WUqm0RMMctzwEWBervFZFeIjIWKATWuO2dEpHp7m6t2YE2xhhjukDa7ilVbRCRR4Al+LfNPqeqW0TkQff5fPw7mT4HlOHfcvuVVG3dph8HForIA8A+4IuuzRYRWYh/sbwBeFhVY67N10jccvs6nXERvHU6vQssS9hxSs+OUXp2jNLr/G75bH+4zxhjTMexJ8GMMcaEZknDGGNMaDmfNESkt4isEZFNIrJFRL7j6v9NRLaJSKmI/E5E+gfaPCoiZSKyXUT+KFA/RUQ2u8+ecBfscRf1f+nqV4vImK7+nu3R0jEKfP53IqIiMjhQl1PHCFIfJxH5ujsWW9xoCI31OXWcUvy+3SAixSKyUURKRGRaoE1OHaNGIpIvIhtE5FVXHigib4nITvc+ILBu1x0jVc3pF/7zHxe75R7AamA68FmgwNV/H/i+W54AbAJ6AWOBXUC++2wNcLPb5uvATFf/EDDfLd8L/DLq790Rx8iVR+Hf6LAXGJyrxyjNz9Kn8YfD6eU+G5qrxynFMXoz8B0/B6zI1WMUOFZ/C7wMvOrKPwDmueV5Uf2flPNnGuqrdcUe7qWq+qaqNrj6YhLPiMSHOVHV3fh3jE0T/1mTfqr6gfr/Eo3DnDS2WeCWfw3c2ZjxM0FLx8iVfwj8A00ftMy5YwQpj9PXgMfVHxoHVW18JinnjlOKY6RAP1d/KYlnsHLuGAGIyEjg88Czgerg91pA0+/bZcco55MGxE8DN+I/YPiWqq5OWuUvSNzem2rIlJaGOYm3cYmoBhjUkd+hszV3jETkj4EqVd2UtHpOHiNo8WdpHHCb6wZ4R0Q+7lbPyePUwjH6BvBvIlIB/DvwqFs9J48R8CP8P8aCo6E2GXoJCA691GXHyJIGoKoxVb0B/2ximohc1/iZiHwL/3mRlxqrmttEivpUbTJGM8doEvAt4J+aWT0njxG0+LNUAAzA74b5e/znk4QcPU4tHKOvAX+jqqOAvwH+y62ec8dIRO4GDqvqurBNmqnrtGNkSSNAVU/gj7Y7A0BE5gB3A/e70zto2zAn8TYiUoB/+n2sU75EJwsco1n4/aebxB+BeCSwXkQuI8ePEVzws1QJ/NZ1zazB/+txMDl+nJKO0Rzgt+6jX+GPrg25eYw+Afyx+716BbhDRH5Oxw691OZjlPNJQ0SGiLszSkT6AJ8Btok/edQ3gT9W1eAUem0Z5iQ4ZMqfAssCSajba+EYbVDVoao6RlXH4P8QTlbVg+TgMYKWf5aA3wN3uPpxQE/8kUhz7jilOEb7gdvdancAO91yzh0jVX1UVUe636t78eP/Mzp26KW2H6OwV8yz9QVMAjYApcCHwD+5+jL8Pr+N7jU/0OZb+HcobMfdjeDqp7pt7AJ+QuKJ+974fz2V4d/NcGXU37sjjlHSOntwd0/l4jFK87PUE/i5q1sP3JGrxynFMboVWId/F9BqYEquHqOk4/UpEndPDQKW4ifUpcDAKI6RDSNijDEmtJzvnjLGGBOeJQ1jjDGhWdIwxhgTmiUNY4wxoVnSMMYYE5olDWOMMaFZ0jDGGBPa/wdEixMYPsKLEwAAAABJRU5ErkJggg==\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydd3xb1fn/39q2JS9Z3nvG21keGSRlFQotBQJtCav9toUu6A/oHkBbuljfUqCD2UWgQPstowtoIZCdeDvx3lu2JQ9t6d77+0O2bHkknklo83m9FMWS7rlHV/dznnOe5/M8RyZJEudwDudw9kF+pjtwDudwDvPjHDnP4RzOUpwj5zmcw1mKc+Q8h3M4S6E80x04h5VDJpNpAB0QDIjABGCRJMl9Rjt2DiuC7GTeWsEjnnPlngGYzWaamppoaWmmqbmZrs5OxifGsVgsTExYsFotTExMMDExgcViQZIkgoOD0el0uN0e7HYbFosFpVLpez1YF4xWpyM4WIdOpyM0JJTUtDSyMjPJzMwiIyMDnU53pr/6fyUUSrlsvtfPkfMMwW6309LSQnNzE03NzTQ3TT43N+F0OsnKzCIzy0uc5KQkQkJC0AUHo9PqCA4OniZdcDAajcbX7nv73uW87TuQJAmHw+EjsI/MVguWiQlGR0dp7+jwnr+xiZbWFvR6PZmZWV7CZmWRmZlJVmYWaWlpqFSqM3i1/rNxjpxnGDabjQMHD7B371727t1LVVUlycnJZE5arilCZGVmER0djUw27+91SkyRc6kQRZHu7m6aW5ppbmqeHDSaaG5uxmg0UlZWxs6dH2Dnjp1s3rz5HFlXEefIeZoxm4zV1VUUFa1n586d7Ny5ky1lW9BqtSdtQ5IkRFHA5XbjcXtwu124Zz573Ljdbjxu9+Rn3IyOjRKhj0ChVKBQKFEqFCiUShQKhe//apUKrU6HNkiLQqE45XcxmUzs2/ce7+zdy7t799La1nqOrKuIc+RcY0iSRHV1NX/+vz/7LONiyeh2uyenneOMTz57PB4AlAolSpUK1cyHUolKpUalUvpeUyq9zwcPHaCkuBRBEPB4PAiCgCB48Hi8z4Ig4HS5sFosWG1WRFFEowlAp9Oi0+rQ6nTotDo0Gs2C1nshsp5//gVcs+sa0tLS1uw6/yfiHDnXCJ2dnTz/wvM8//weLBYL1177MS688EK2btk6LxntdjsjIyNMTIwzMTGBw+mYdNyEEBwcTEhwMLrgYNQq9bL6s9RprSRJOJ1OLBYLFqsF6+Sz0+lELpej1erQabXodMHo9Xq/9e0Upsj6xhtv8Kc//4mMjEyu372ba665FoPBsKzv8d+Ec+RcRZjNZl5++SX2PP88J04cZ9fVu9i9+3q2bt2KXO4fOhZFEbPZhHFoiOHhIdQqNREGAyHBIQSHBBOgCVj2+nI+LHfNOR8EQcBqtWKZ9A6bTCMIgoBeH4HBYCBCr0eh8I/Gud1u3njzDfbseY6///3v7Nixg93X7eYjH7mCwMDAVenXfxrOkXOFcDgc/O1vf2XP83v497//zcUXX8z1u6/nkksunWNNHA4HQ0NGjEYjFqsVfXg4UVFRREQYUCrXNrS8muScDx6PB5PJxPDwECMmE0qlEoPBgMFgICw0zG+gmZiY4C9/+T+e2/Mcx44d46NXfJTrdu/m/A+cv6i17n8LzpFzmRgYGOBnj/yMZ555msLCIq7fvZurr95FaGio7zOSJGE2mzEOGRkeHkKhUBIVGUVUVBQ6nW5VLeOpsNbknA2n08nw8BDDw8OMjo2iD9eTlJREaGiY3+f6+/v54x//yJ49zzFoHOS2227n1ltuJTg4+LT19WzFOXIuEW1tbTz00IO8+NKL7L5uN3fccScpKSm+9yVJYmjISE9vLxMT44SFhRMVGYXBYDijnsvTTc6ZkCSJoeEhuru7sNvtxMclEB8fj1rtv36urq7mgQfv58033+TWWz/HbV+6jcjIyDPS57MBC5HznHxvFmpra7n//p/yxptv8NnP3sLxuhNERUX53nc6nXT3dNPX10t4uJ70tDRCQkJPq3U8WyGTybwzhsgoXC4nPb29HD5yCK1WS2JiEoYIAzKZjKKiIv7w++doaWnhoYceJDcvhxuuv4E77riTpKSkM/01zhqcs5yTOHDgAD+9/yeUl5dz++1f5tZbbvVNXb3TVhMdnZ3YbFYSExKJj09Y8/XjcnAmLed8kCSJ0bFRuru7GB0dJSYmlsSERD/nUF9fHz975Gc8++wzXPGRK/jqV79Gdnb2Gez16cW5ae0CeOtfb/GjH/2Irq5OvnLXV7j55k/6bhy3201vby/dPd3odDpSkpMJCws/q63k2UbOmfB4PPQP9NPT3Y1cIScxMYnoqGifc8hkMvGLX/6Cxx9/jO3bt/Ptb32H9evXn+Ferz3OkXMWOjo6uPOuOzh+/Dh3f/duPv7xT/gsoSAItHe009fXS3xcAomJiXPWTWcrzmZyzoTFYqG7pxuj0UhKcjKJiUm+MJTVauWpp57i/gd+ypUfvZIf/OA+9Hr9Ge7x2mEhcv7X5XM6HA7u++F9lJaVsHlzMdVVNVx//Q0olUokSaK7u4v9B/Yhl8nYtnU76enp7xtivp+g0+nIyc5h65atOJ1O9u1/j+6ebiRJQqvV8uUvf5m62uMoFAoKCvP59RO/RhTFM93t04r/Ksv517/9lTvvvIOCggIeevBhkpOTAe+6yGg00tTchMFgICM9432rFX2/WM7ZcLlctLa1Mjw8RHpaBrGxsb7lQ2VlJRd/8CIMEQb+8Ifn2Lx58xnu7eriv3pa297ezp133UF9fT0/+9kjXHrJpb73zKNmGhoaCAwMZF3WurNexSJJEoIg4HZ7Re8ut8v7f5cbt8dNe3sb6ekZ00L3GQ+VSk1QUNAcFdPZBIfDQUtrC6OjZjIzMomK8mboHD58mEs/dAkajYarr7qaH/zgPiIiIs50d1cF/5WhFIfDwQMPPsCjj/6cO+64kxee/6NPzWO1WmlorEcQBPJy8wgJCTnDvfUSz263+wTwVqsFt8tLQEEQkCFDQkKpmBS8q1WTAngVapUKjdorVlcoFAiCgMvlmhS+ex9ulwub3YYoigRoAtDqtD7trFarIyBgdaWEy0FAQAD5efnY7Xaam5tobWslP7+A0tJSvnfv9/ntb3+DTC6noDCf73/v+/zP/3z6rB5sVoL/WMt55MgRbrr5RgoKCnjwgYd8U1in00lTcxPj42OsW5eNIeLMCLMlScJisWAyjfhlogQGBBIcEkJIcAg6nRaVWo1apUIuVyyKOIuZ1vrE7lYLVqvVm6FiteJwOEAGQUFBPtKGh+sJCgo6Y6QdHx+npraG6Kgo0tLSufGmG1ApVfy//3cHX/zSF5iYmOCVv7z6vs6E+a+Z1oqiyM8e+RkPPHA/jz/2OFdfvcv3elt7G319vWSkZ/qtaU4XBEFgeGQYo9GIyWRCp9USEWEgJMSbkbIa69yVrjlFUcRmt2G1eAXvJtMINpudsNBQIgwGDBEGAgICVtzPpfaptbUFo9FIZmYmH/7wh7n8wx/mzTffRKVU0tDYwGOPPsauXdec1n6tFv4ryDk8PMz/fPpTDA+PsOe5PT65ndVqpbqmCoMhkoz0jNM6DbLb7RiNgwwajTgcdgwRBqKiotHr9WvSj7VwCImiyNj4GCPDwwyPDON2uwkP1xMbG4s+XH/aBrnx8XFqa2vQaDRcdfVVXPahy/jlL3/FH/7wB27/8m1cf/0NPPTgQ6d98Fgp/uPJ+d6+97jxxhv4+Mc/wX0/uA+VSoUkSXR1ddLV3UVBfiFhYWGnbmiFEEWR0VEzg0Yjw8PDqNWqSRF89CkrH6wGToe3VhAETGYTvb29jI+PERcbR3xCAoEBa+9ME0WRtrZWevt6KSpaz5tvvMFtt9/Gww89zKuvvUpzUzN79jzPunXr1rwvq4X/WHIKgsBP7/8pjz/+GE8++RSXfegywOsMqqmpJkirJSc7Z81TlMbGxujo7GB01ExYWDjRUdEYDGufIjYbpzuU4vF46Ovvo6enG4VCSWJiop/qZ60wMTFBTW01Tz31FDfddDPnf+B8JEniySef4O577uahBx/i+utvWNM+rBb+I721AwMD3PzJm3A6XRw+dISEhAQA+vr7aGlpJicnl0jD2mU7CIJAf38fXV1dqDVqUpJTKCwoPOMez9MJpVJJUmISSYlJPtVPc3MTEREGkhITCQkJPXUjy0BwcDBbyrby9ttvMzIyzNjYGKGhodxyy62Ulpax+/rrePvtt3nkkZ+flhnLWuB964N+Z+87lJQWU1a2hbfefIuEhATcbjeVVRUMDAxQVrZlzYhps9mob6hn3/73sNpsbNiwkc2bijEYIv+riDkbU6qf87bvINIQSVNzE/v2v0d7Rzsul2vVzyeXy9m58wM8/vjj1B2vo7GxAVEUKSoq4vChIwiCQNmWUhoaGlb93KcD70vL+dJLL3L7l2/nd7/7PRdfdDEAdruN8opyUlNSiY9PWJPz2u12GpsasNlspCSnsi5r3X9sjG0lkMvlREdHEx0djdPppLe3h8OHD6HX60lPz1hVh832bdsxmUw4HA7kCgVHjh5m44ZN6HQ6nn32NzzzzNNceNEFvPzSn9iyZcuqnfd04H235nz88ce5/4Gf8tqrr1NYWAjA2Ngo1TXVFOQXEh4evurndLvdtLS2MDw8zLqsLCIjo85aC3m2yvckSaK/v4/W1lYiDAbS09LnLRa2HPzyl7/knb1v88cXXmRwcICmpiY2bNiATuetsvD3f/ydT33qkzz99DNcftnlq3LO1cT73iEkSRJ333M3L7/8En/7699JTU0FYHBwkKbmRjZu3IQ2aHXXFqIo0tnZQVd3N6mpqSQmJJ61pJzC2UrOKUiSRG9vD23tbURFRpGWnr7sSoNTGB8fJy09lZrqWuLi4hgfH6equpKc7BwiI72J8ocOHWLXNVfz4x/9mJtuunk1vsqqYSFyKu69994FD5JEaeE3TyM8Hg+f//znOHjgAG++8Rbx8fGAN+2ru6eL4s0lq+rGlySJ/oF+qqor0QZpKSwsIjz87M7jnEJXVyfJSclnuhsLQiaTERISSmJiEna7g7rjtbg9HsLCwpZ9fTUaDZ2dnTS3NLNzx040Gg0xMbEcP34cj8dNWFgYiYmJXPahy7j1c7fgcrnYumXrWfN7yuWy7833+llvOe12O7uvvw6Hw8FLL76MTqdDkiTqG+pxOBwUFRatqtveZBqhoaGB4JAQsjKzVm3qtVaQJAmPx+MTwFdWVrBp4ybkcgUKhRy5XIFysuL72QhBEGhta2VoyEhBfuGyNc51dXVc/MGLeO3V131ZK4IgUFdXi0wuIz+vALlcTnd3Nxs3beBjH/sYj/78sbPCZ/C+nNaazWauvOpKEhMTeObpZ1Gr1QiCQGVVJcE6HVlZ61Zt9LNYJqhvqEeGjOzsnDO+45bb7WZ8YhyrxeolnsuFy+3G5fL+3yN4fJ9VKlWoJ0Xwg4MDREVFIwoCgiggCiJujxtRFFGpVGi1WrRBOrRar342KCjorCDulPonMipq2SquV199hVs/dysvv/Qntm3bBngHr7b2NoaMRtLS0vnSbV+kpbkFrU5LSnIKTz/9zBnP133fkbOvr4/LLvsQF1x4IQ8+8CByuRyn00l5xTESEhJJSlydQlCiKNLU1MiIaYSc7Bz0+tOfhiQIAiaTCZPZxPj4ODabDZVSSXBICDqtDrVGjVql9hFQrVajUMwvhD/ZmtPlcnmF7lYrVpv32Ta5JYNKpfYSV6slLDSUsLDw007aKf3zwEA/+fkFhIUuXdH1zzf+yc0338SePc9zwfkX+F7fu/cdqmuq6ezs5If3/QhRFNl9/XU4nU7fjOxM4X1FTpPJxM4P7OC663bzzW98E5lMhsVioaKynJzs3FUrozgxMUF1TTWxMbGkpaWd1jWI0+nEaDQyaBzAZrOh10cQoY8gJDSEoMDlZ4EsxyEkSRJut3uSuBbMo6OMms2o1CoiJqu7h4aGnbYpoMUyQU1tDfpw75aESx0k9r67l49//GPs2LGDjRs3YbPZeOqpJ/nVL3+FIdJASXEpSqUSj8fDZz/7GQYGB3jlL6+eMQv6viGnzWbjg5d8kK1bt3L/T+8HvGStO17L+qINq5J3KUkSHR0d9PT2UFRYdFpyOSVJYmJigkHjIENGI8ggKiqa6KjoVS08vZreWofDwcjIMMMjI4yOmjEYIklMSDxt16u9o53e3h7y8wqWHCLr6enhvffepbyigpHhYb7//R+QmJhIb18v3V1dFBcXo1B4CXrNtbsICQ7hN7/57RlZg74vyOl2u9l1zdXow/U888yzyOVyRsdGqampprh4dTyyDoeD6poqgnXBrFuXvaZTN0EQGDGNMDg46E0R0+mIjoomKioStXptHE1rFUoRRZGhoSG6e7pwOBzEx8VPFoxeW4eZ1WqlpraG0JAQ1q1bN2dvluVgqu7w5k3FKBQKbDYbl37oUoqLi3nwgQdPuxf3rCenJEl85jOfZtA4yP/9+S+oVCoslgkqKivYvKmYoKCgFZ/DOGSkoaGe3JxcDGuoubXarHR2dDA0PITBYCA6KmbNUsRm43TEOV0uJ729ffT29RKg0ZCYmEhkZNSafb+p7KLOrk6KCtf7bYWxXHR2dmAcMrJp42bkcjkmk4nzL/gAN954E1+56ysr7/QScNYL37/9nW9zor6et958C5VKhd1uo6Kygg3rN64KMbu6uujp7aa0pGzZ4ZH5BrKpUda7PcMQHZ3tCIJISnIy2dk5Z4WrfrWhVmtITU0lNTWV8fFxenq6aWhsINIQScIaTHtlMhnJySlEGAxUVlaQkZ5BbGzcitpMTk5BEEUqq7z3mF6v56+v/40dO88jKjLyrBAqnBXkfOSRR3jllb+w95130Wq1OJ1OjpUfoyC/cMUb3UiSRFNTIxMWC6UlZcuexi40wxBFkaHhIZqaGgkLDSN7XTbBwSELTo0kSTprgt+rgZCQEHJz8xBFEeOQkcamRtxuN1mZmURMbr+wWtBpdZSVbqGysoKJiQkyM7NW1H5aahotgkB1TRXrizaQkJDAX1//GxddfCERBsMZl/qd8WntCy88zze++Q32vvMuycnJuN1ujhw5TNa6dSvOKhFFkeqaatRqNbk5uSv6Iee7TmPjYzQ0NKDRaBas3DfTss5+bS1wNsj3LBYLzS1NOBwOMjOzVr1OkyRJ1NefwO6wU1S4fsU5s41NjTgcDl+636FDh7jq6iv5vz//hbKyslXq9cI4K9ecb771JjfffBNv/PNN8vPzEQSBI0cPk5KcSmxs7IradrvdlFccIyoqmrTU5Rd/mu/6OBx2GhsbsTsc5GTnzNhTBRbi3ez31oqgZwM5p2CxTNDU3IzL6SQnJ2fOtoArRXdPN52dnWzauJHAwOUvfSRJoqGhHo8gkJ+Xj0wm4+//+Duf+cyn+ddb/17zfVvOOnJ2d3dTWlbCCy/8kR3n7UCSJI6VHyM6OnrFAgO73c6x8mNkZGQQG7N8ks++Nh6Pxys1MxrJyMwiajI7ZYpnkuT9RyaXTf89D9aSpGcTOacwPjFOXV0t4eF6spYRtzwZpsJsmzZuXlFStSRJHD/hrTCfk50DwDPPPM3PHvkZBw8cWtOE7bOKnB6PhwsuvIDLL7+cr3/t6wC0tLbgdrnIycldUdtjY2NU11StKH1s9jWRJImu7i46OztISvTf12PWB+dvcB4CrhVBz0ZywozYck83uXl5RKyiEmvqN1+pV1+SJMrLjxEXF09cnNfh9MlP3oxKpeLJJ59are7OwVm1V8q937sXrTaIr37lq4B39DMaB1m3bmXTB+OQkZqaajZt3LxqxBwaGmL/gX3Y7XbKSreSnJSMfAEyLWUom3makw2QS4EkSavW1mpDJpORmprKpk2baWlpprauFo/Hc+oDF4HQ0FAKC4o4Vn4Uu92+oj4WFa2ntbWFiYkJAB577HH2H9jPc8/9YVX6uqT+nG7L+cabb/CZz3yaY0fLiYqKwuVycejwQYo3F69o3dDd3UV3TzebN21edmB85rXweDycqD+By+UiNydv2tkz+RlRkkACuUKONHmZZl4suWyarFNUls0/QM67Tp2ypoIgYLVZcTlduFxOnC4XLpcLl9OJy+XC6XLh8biRIQOZ1xmj0+qQJnsjk8m9RaKDgggK0hKkDUIbpD2j1d0lSaKnp4f2jjay12UTFRW9Ku2azSbq6uooLi5ZUbWF8fFxqmuq2FK2FaVSSXV1NZdc+kHe3fseWVlZq9LXmTgrprX9/f2UlBbzu9/93lct7Vj5URITk4iJjll2uz29PfT29vgUH0vF7GswOjpKbV0tKcnJxMfPSrCWJC8xfX/P8MhOvSRKfkScfeXnI+nUKabq3JpHzVgsFmQyGVqtDo1Gg1rtFb1r1GrU6um/lUqlrw+zp7WCIGC327HarNisNmw2KzabDYfDgYREgCaAIK2WkOBgIiIMp7W6u8PhoO54HQqFgrzcvFXRto6MjHCi/jglxaUrSvfr6e1hcHCQjRs2IpPJ+PWvf8UTTzzB/v0HVr0u7hknpyAIXPqhS9m+fTv33H0PAG1trdgdDvJy85bdrtE4SEtLCyUlpct2qU9dA0mSfJXFCwvXe50AU9dHJkPwiMhmLwT8pqfT/z8ZOWe+L4kiE5YJBgcHGRo2olQqiY6KRq+PQKfTLXmwWcqaU5IkHE4HNqvNWzR6ZASbzUZISAjxcXEYDJGnRUTR399HU3MzmRmZvrXeSjA8PEx9Qz2lJSUrkhfW1tWi1WpJS01DkiSu2/0JIg2RPProYyvu40yccXLe98P7ePvtt3njn2+gUCgwj5o5ceI4W8q2LvsGMJlMHD9RR2lJ2YpHXY/HQ2VVJTqdjsyMLO+6Uibzm8ZKM7aH9JFUAmHqMknSnO8ik8vmkFMQvUWZh4yDmMwmdLpgoqNjiDREolbP3ZJhKZZspQ4hSZIYGxujp7eHkZHhSdVPwpqVuJyCy+Wi7ngdKqWSvLz8FQ8KxiEjTU2NlBSXLvveEASBw4cPkZ2djV4fwdjYGMUlm/nxj368qls/nFFy7n13L9dfv5sjh48SFxeH2+3m4KEDK/KuTdWJKd5csqxt+2Z+b5fLRXm5N080ISFx6gMAflNYadbereLsazfj75k3l1wuQxAEBo2DDBoHsFoshOv1REdFEx6uR6lQgEzm65N83mnv4gi6mt5ar9jdSE9PDza7nbi4OOLj4tdsuwNJkmhvb2PQaGTjho0rrkIxODhIS2szJcWly96Hxm63c+ToEUpLSgkICODo0aNc8dGPsH/fgVXbPOmMkXN8fJzCogJ+9atfc+kll3rd1RXHiI+LX7Y+0mazcfTYUTZt3LSiJFlJkrDZrJRXVLAua52vGNTkm8DJyQng8QjIFXK/YzyC94NqlRKb3U5vXzdG4yCRkZHERscRHBLsRza5TObnTJIxd116Jsg5Ey63i/4+r9hdoVCQEJ9ITEz0qmSJzMZUgsL6ovUrttgDA/20tbf5cjiX1R+jka7uLjZv8pY/eeSRR3jxpRd57933VmXaf8bIeddX7mJ8bMwXJ2rvaMdisVCQX7Cs9jweDwcPHaCgoHBZmfIzYTabqa2tobCwiNDQUP8w5eQfgjDp9ZT7k1OYJKDP2inkeDyC7323201nZxtjY2YSEpOJiY71bmIrl+F0edBovDeKKEggA4Vi+kc+G8k5E1arlZ7eHgYGBoiJjiY1NW3VE5UtFguVVRWsy8omKirq1AecBH19fXT3dFFSXLpsZ1dFZTnx8QlER0UjiiI7du7gU5/8JJ/+9GdW1Dc4Q+Ssra3lg5dcTE11LZGRkYyNjVJbV8eWsi3LVolUVVcRodeTuEIV0aBxkMbGRjZt3OSbWksS02ERyUua2etMSQS3W/BNPf2nxwJyBfT399Db10NiQjIx0XHeG2LW5VcpZ3z/yfcUCrnfx5ZD0NMpQhBFkZ6eHjo624mKiiYjPWNV94ZxuVwcOXp47qxmGWhsakQSRbIn1T9LhcPh4MjRw2zbug2FQkllZSWXf/gy6mqPo9frV9S30y5CkCSJ2798O3d/9x4iIyMRBIGamhrWr1+/bGL29PQgieL0unCZ6OrqorW1lZLiUgIDg3wWU5o9Fs360+MWcbu91lGc9VmXS8BkHqa8/DBOl4uNG0qIjYlHFL3WV/D4f949w8r6Tid5Y6NTj9n9OdtEBnK5nKSkJLZvO4/AwEAOHNzP4ODgqrWvVqsp3lxCQ2Mjw8NDK2orKzOL8fFxBgYGlnV8QEAAiYlJNLe0ALBhwwau2XUN3/3ud1bUr5Nhzcj5wgvPY5mY4JZbbgGgta2V+Ph4dNrlrREtFgvt7W0UrGCjIEmSaGxqxDg0SElxid9U7FT3/NQ0dgout4DD6cFqczM6Ok7d8UqMQwPk5hSRlJAGktw3JT4ZRElCFKcecxe1cwYMVk9RtFqQy+UkJyVTUlxKT28PFRXlOJ3OVWlbo9FQUlxCfUM9w8PDy25HJpOxfv0GmpobsVgsy2ojJTmFkeFhn3roe9/7Pq+8+grl5eXL7tfJsCbkHB8f5+vf+Do///mjvjIQg4MDpKR4q7Q3NTXx+l9fX/TmNoIgUFVdRWFh0YpimTU11bjdbjZu2DSv9ZbwN5aC4LWU8xETwO120d7eSFNLPcnJ6WRn5aPRTHsynS6P7wH4rKfHLeJxi9jtbjxu/7ZFUZxD0vkIejYiICCATRs3ER8fz+Ejh+ju7l6VgcRL0FLqG04wMjKy7HbUajVFheuprKpclnRQJpORl5dP3fE6JEkiPDyc+35wH7fdftu8A+tKsSbkvPd793Le9vN8uXD1DSfIXpdDY2Mju665mp0f2MGVV36UX/ziF4tqr76hnvi4+GWXp5Akidq6WgICAxfO65zplWXaUioUMp9lAy8xvRXhe6g7XkFISBi5ORsQJQ1jE15r4XB6cDg9flNfp8vjLQDtnvsjzvfafAQ9ywzmgoiOjmFL2VbGxkY5cuQwVqt1xW1qNBqKN5dw4sRxTKblEzQ0NJTUlBRqaquXNXCEhYWh0+no7bEWs4EAACAASURBVO0B4IYbbsTtdvHsb55ddp8WwqqTs66ujuf37OEzn/0sh48corOrE0mUiIyM5Hvfv5fU1DRamlupqz3OT376Y4xG40nbGxjox263+baQXw6ampuQyWRkzZc5L0mnntNOQhQl7DYnDQ012KwWCguKiYiIxu6YHoUdzrkj8hRZrXb3gm3PR1Bp1mOqu+8HqFQq8vMLyMrKoqKynNbWlhVbl4CAADYXF3P8+HHMZtOy20lISEQhV9DX37es47PXZdPa5q2ve/DQQb71zW/x3e9+B5Np+X2aD6vqrZUkiQsvupBrdl3DF77wBUZHzRw8dJCY6Bjy8vIoLinmt7/5HevXrwe8YRab1covf/mreduz273xzLLSLct21Xd0djAyMszGDZuQzQj0ezkp+XlHp76tOCtMMrV2HDGZaG6pJyE+laioaAaNVgID506zJUkiQDP9+tQ0eCpcog2cDojPDKGo1NNTbYVc5uetnaot63Q5vaJ3twun04nT6f3bOfnalPBdpVYRFBREUGAQgVPPgYFoNJrTLngXRZGWSVnk+qL1Ky7gbLfbOXrsyIrSAl0uFwcPHWTrlq1LFihYrVYqKstxu9yUlW0hKCiI2277EpIk8dhjjy+5L6cllPLii3/k/vvv5/DhIygUCtrb23C5XAQHB9Pc0sxjjz3Gb579rS9xdXR0lLz8XP76+t98hJ2CKIocOnyI7HXZy3ZV9/f30dnVSfHmEr81pjgzXDIJGZNe0hnry6nPOZ0eRkxDdHW3k5nuXVdabN718nxqnpntjo05CQ/3rkNnEnGKoAqF3DeFVqrkKCdDLDIkRkwjGI2DjE+MI0kiapXaK4DXaNBoNARoAtBM/n9KGL9v/3ts33YeLpcLm92G3WaffLZhs9t8jhqNJgCdTos+XE+EwbDinb4Wg/HxMaqqq1i3LpvoFWaiTA3chYVFy4539/T0YDabKCgoXNTnPR4PzS3NjIwMe61naytZWVmEh+sxm83kF+TNey+fCmtOTlEUSUxK5IMfvJhnnn4Wl8vF4SOHfHGhiYkJbrv9S3zqU58iJzvXF1h+8skn2PP88/z7X//2G9EbGupRKJVkZmQu+kvOhMVqobKigrKyLX4j4+wQyMzvL81YW878bHdPF/39feTmFAEKHzEBnM7pkMiUFW1pMWEwzJAlyiA8LMCPnE6ngEYzPWBoA1UgEzGZhhkxDWG1WNDrI4iKiiYsNAylytu2n7Jont90MXFOSZJwOp1YLBOMmEwMDRkJ0AQQHx9PVFT0mtbydblcVFSWY4gwkJ6esSIrbrPZOFZ+dNkVFSVJ4vCRQ6zLyj6pBZ5KcWtrbyMlOZmkpGTvLgRWC1VVlWzdsg25XM4jjzzC/gP7ePGPLy2pH2tOzldffYU77ryD0JBQPnLFFVyzaxcRBgNxMyR6Hzj/A3z9618nOjoaj8fty5PcvHkT9/3wh75qZybTCE3NTZSWlC3rxxNFkYMHD5CXnz/vqDqTgL5p7uRrgjg91XU63XR2tWOxTJCenoPdLmK1+a8bZxPEOGTF7ZomrMEQ5CdAsNu9a9KoSC0ajQKPx83Y2Aijo8MIght9eASRhmgiDOELEnHq9eWScz6MT4zT29uL0WgkLDSU+PgEIiIi1mQKLIoitbU1aAICWLfCzaiGhodoaWmmtKRsWVI6i8VCVfU0wWbDbDZxov4EYWHhZGVmzZkCNzY1olapSU1NxWq1kpGZztv/fmdJdYfWlJySJLFt+zbuuvNOduzYyRUf/Qh33HEn115zrd+F/+GPfojZbObBBx7EZDJRX3+C8PBwampq+PUTT/D2v99GkiT2H9jHhvUbl1235cSJ4wQEBJCWlj5vX2dXIZgZqpjKMBE8Ag2NJ5CQyMrIwTw6HbebSdBBozdmFhsTjHHI65WcSU4AQ2QQZrM3Qz8gQIUkSbhcY0jiKKIoEBoaQViYgeDJtVhggPcGUCrlKJTeG2Y2EadqF82+sVcjK8VkMtHb14PZbCYyMor4uHhCQhYu97nc89TU1vgqF66k7ZbWFpwOB3l5+cs6vqmpEYVCSXr69P1it9tpaKjH7XaTm5vr2yV7NjweD/sP7Gf7tu0oFAp+cN8P6Ozo4Kmnnl70+ddUIfTO3ncYHTVz5ZVXERkZyWOPPsarr77Cz3/+c7/PXXThRbz15psA6PV6tm7dRnBwMLGxscTFxXLw4EG6ujqJiDAsm5hG4yAWi4XUWRX3FlLXzBdDdLvc1NZVoVZryMrIYWjEjkcQfYL2KUwRE/ARE/wdO/09Zo5X99HXZQbA5Rpn1NyE1TKK06knODiDiIh4FIrpaZndMU1+wTPlnPJ/IE2vm1dTlCCTyYiIiKCwoIjztu9AH66npaWZgwcPMDg4uGrnkslkFBYU4nQ6vdK6FbSbnpaOw+mgp6dnecenZ9DX14vNZkMQBJqbmzh27Chx8fGUlJQuSEwApVJJbEyML7TyxS98kVdefYXu7u5l9WUmVsVyXvqhS/j4xz7Opz71PwiCh33795GclMJFF1/IV7/yVW699XOAd5QJDtExPDTiRz63282LL/4Rl9tFcnIyO87buSyxwZT+cb41iN/aUpp+bTY5rVYHtccriYmOIzoqzmclbbPCIINGK84ZYZPxUa9l7OkaBSA7L4b+Hi8h5XIZMrkHTeAYGo0alSoKmcxrHcPCvOlukZFBBAX6T5mCddPfQSYDxaSzaGqcnU97u1baWpvNRnNzEw6nk4L8glWpwg+rZ0E9Hg8HDu5fdgrh0PAQDQ31Xt9JQiIpKamLniZPldo5b/sOZDIZX/3aVxE8Hh5++H8XdfyaWc7y8nIaGhq4/vobAOju6SE+Lp6UlBT++Y83+PFPfszvfvdbABoaGkhMTJxjFVUqFVdeeRXNzc14PALVNdXYbLYl9UOSJKqrq8jNyT2lc2A+ayOIEhMTNqpqjpGUkEJcbMKc9eUUBo1eK2m1OLFapqe7U8QEqDo2NYqLqAPGUKqGMfYrGegJoLt9fE6bQ0Nzv6/F4kIUJV8JE2GWHnc+7e1aISgoiKKi9WRmZFBecYyWVYhbwupZUKVSSW5OLnXHa5fcxtjYGC0tLXg8HtJS00hLS1/S+lWtVmOIMNA/GTe94//dwe//8PsVyQ1hFcj50/t/wh3/707UarVvw5mkpGQA0tPT+cff/8ldX7mL/v5+3z6b800/RFEkJzuHd999l6TEJMorjtHY1LhomVVrWyshoaGL2qDIu16bHqwEUcLjcVN3vJKMtHUYDFG43IJX9zr5CAjwWvIpYs6E1eL0I6YXEna7mYAgI2OjAr1dOhx2JRMTXk9vV9u0ysVstmM226mqHsDpEnyP+bBAjbDTBr0+gm1btyOJIgcO7l+RGGAKq0VQgyESlVK1aHG70+mkpraaE/UnyM3JZeuWbbR3tC9L2peWlk5bexuSJBEXF8euq3fx6GOPLrmdmVgRORsaGnjvvff4zGe8OW2Dg4Po9RF+goHs7GxuuP4GHv7fh8nLy+P227/MjTfdiCD433zVNdX86te/4qNXfJTIyEi2bd2OWq3mwMH99Pb2nvQHM5vN3tKaWesW/MxC0yWZXIYoClTXVJKcnEp4uN4nGpiN+YjZ02GmusJfaaLWuElMsaIOEOju0DExpmF2zpjHI9LWNORzFAEEh0xb/AC1dzBQKKYLVE8W/PN7zA4NnQ7I5XIyM7PYsH4jTU1N1NbW4HYvrH5aDFaLoLm5uTQ1N520P6Io0tbWyuEjh4g0RFJWWkZoaCgajYaU5FQamxqXfF5vdQgZTz75JBdceAHv7H2Hn/3v/zI+PneWtFisiJwPPHA/X/zil3zT1Lb2NlJTU+d87vzzz+eZZ55GkiS+9tWvERgYwMUfvJijR48CYB4109bWSl5uHsXFxd6OyeWkpqRSVroFk9nEwUMHGB2dbZ2869XauhrWF2046VRkpjJoNuobThBpiCI6OnZBYsLczJSeDu+aMkIfyOioA5AIDbMRFm5joC+Q4cFA+hpH6Gvy14J6PNPttDf7p0L19Pr/mEtQF552aLVaSkpKCdfrOXDwAH19fSuaWvsI6nDQ2tqyrDbUag1paWk0NDbM+/6gcZD9+/chCALbtm4nNjbOb+BOTEzEZDJhdyy+/m1/fz+lpSV88YtfQK1Rcccdd/DM089SWFjIE08+sazvASsgp9Fo5JVXX+GLX/gi4CWYWq1GG+S/nvz973/HLbfewm+e/Q0ymQy5XM6rr7zGdZ/4BLuuuZrrdn+C2poaKisq2blz55zzqNVqCvILyM8voKGxnqrqKhwOh+/9E/UnSE/LOKmD4mQ3TF9fL6IgkJiY7POMTmGmg2ZszLu2FARxDkkBZDKJ0PAxHA4PxoEQuo+bGWieXnP0NY3Qeqib1kPd9PdNzDk+OESDw+HB4fBgm1zregTRr7oCnH0ZKjKZjIT4BLaUbWHQOEh1TfWcWdFS2ysoKGR4ZJhB4/JyQxPiE7BaLZjNZt9rFssER44c9m6aW1y84Hb2MpmM9LQ02tvaFnUup9PJtR+7lks/9CEqyispLCxi65atbN26lUcffYxHH/35sq/Hssn54osvcvlll/uUFW1tbaTPiCt6PB7u+spd/PBHP+Rfb/2bj3zkCt97SqWSz372FupPNLBhw0aOlR/jjy/+8aSWLyQ4hNKSMmKiYzhy9DAtrS3eEc5uW3Y5RZvNRlt7G9k5ed5yITAnZjE25vQRcyYEQfStC2UykdR0B+NjCgYH1DTu6/L7rN1kw27yd/j09034Hu3NQzgc/usc0+j0yD0fQc82kqrVatYXrSc8LIzDRw6vaJorl8vZuGEjjY0NWKxLz72UyWTk5xdw/EQdTqeT4yfqqK6pJjMzkw3rN55yh/TY2DiGR4ZxuU6dk/rlL99ObGwM99x9DwqFgrTUNDq7OgFYv3490VHRvLP3nSV/B1jB/pzP7XmOe+65B/De5E6ng/DwcDweD+Xl5dxz790AHNh/cEFtrFar5cMf/jByuZzh4RGKik6uSZTJZMTExBAZGUlbe5tXerXELRymql2KgkBtXTW5OXmolCqf1VSrlbhc00QJDdUwMeH1mkbogxiZJNlA7xgAHsFNbLyVvl4VzftnOHnavVPw8NQw+quNaILV6NO8aqWRxmEi1hl8oRGj0UaQ1kRCiv91Gp9wotN6rbcgiD75nygKOBxOnE4HDqcTp8v7f6fDgc1m84owAgMJDAiYfD49gvepTW5VajVHjx2heHPJsqveqdUaigrXU1VVxdYtSy+fqg3SolFr2PvuXnJzcsjNyVv095fJZKSmpNLe3n7S+6upqYnXXn+NhvpGX/9CQ0OxWCYQBA8KhZLrrtvNnuee48ILLlxS/2GZcc7m5mY+cP5OOju6UCgUvPveXpqbW3j99dd47733SEpMYtc11/CNr3/jpPHKqdjU9m3nLfni9/X1MWgcRC6T4XA4yMnNJSR44R2VZ8c5m5ubAIn09ExEQfJbB06R0zZpzaY8rAC11dPOH6PRjF4/TmcjtB8xoQmZdoTZRqYt32jXOJrg6fcEt0De1f4bNqWkeombkKInLtYb9FapvORVqwRGR4cwmYYRJRG5TEZAQACagAACNAHebJOAAAICNFRUlJOfl4/d4cBht2N32HHYHThdTiRJQqVSER0VTWxs3LLigYvFwEA/rW2tFBeXrEhUP59651QYHhmmoaGeCH0EQ8PDbNq4acmiFlEU2bf/PbaULZy18p3vfgen08kD9z/g93pjYwPBwcHExcXT19dHYVEB3V09C17vVd12fs/ze/jYtR9DqVRy1113kpmVSVVlFbuv280Tv35y0dXSurq7SIhPWDIxBUGgpbXZJzYYHR2lrq6WkOAQsrLWnTK9zGw2MTIyTElxqW86q1TKfQRVq5WMjjvmPVYul03GHl1ERIzTUge2ySWkc9zlR9D+qkEC9d4fxDnhQhOsRjiJwwmgp8NEsE6NVitjaGiYiYkRAjQaIgzRFBRsICjQv2asXCH3k/YpFIqThpOcTicDA/1UVVchSSKxsXHExsSuei3amJhYZDI5R44cmVMSZilIT8/gwMH9xMTEnJJgNpuNE/UnAIkNGzaiDdISbRqh7njtkivvyeVykhKT6OzsIGOe5AtBEPjDH37Pa6++Pue9uLh4GhrrfbuVbdq0iddff41rr/3Yos8Py1hzSpLE88/vYffu66mpqaGru4tt27bzq1/9mo997OOLJqZX6d9NQkIi7e3tSwpod/d0Exsb5xMbhIWFsaVsK+Hheg4dPkh7x8LtuVxujh+vo6hoPfJZDgGlcvpyKOQyFLMGNNNk2EOhdCBTGGmqmibmFJzjLpzjLvqrvM4Mu8mO3WTH3D7mR8zBurnODrncgyZwguHhBrq7m1Gp1CQn55OTs56oyFhEcdbPJfOvmDCfEH42NBoNyckpbCnbwsYNm0Dyln08dOggnZ0di1pnLRbR0dFkZWZx5OjhZdcU8u6jkk9tXc2Cjj2Px0NDYwPlFeWkJKeweVOxzzGp10egUWuWVSAsMTGJ3r7eee+lf7/9b6KjoikomFviNTg42JdfC3D97uvZ8/yeJZ9/yeQ8cuQIcrmcTZs2cddX7uLyyy/H5XRRVVVFZWUlFRUVlJeXU15eftIaQW3tbfT29FJWVkp+QR7f/8H3F3X+KaFD8qTQYQoymYz4+Hi2btmGy+Vi/4F9DA0N+b0vSRJ1x2tJS0uf3tFMNuuBf6hDIZdhMtt9xJTJ7UjSCCPGUKKSDMiVcy/hcLN/YN4xNvc6DDeNMFg3yGDdIAqFhC5kHG3IGAq5Cr0+A70+A1EM9hVtlssmBQiz+qqYMetYahgjICCA1NRUtm7ZRlHReu/2A0cO09BQjyCszvZ8UVFRrMvK5sjRI8smqF6vR6vV0dPrL16ZSuXaf2A/gQEBbNu6DYNh7hb3GRmZtLQsPTSjUCgwRBjmJfbvfvtbbrzxpgWPjYmJZaDfK4a48sqr2Lt375IVQ0ue1u7Z8xzXXbcbh8OB2+3Gbrfz2c96RQgzlTcTlglycnL408t/ntdl/de/vs7hw4e5/4EHyMvNY+u2LawvKuLKK6866fmHhoyEh4UvOE1SKpWsy1pHYkIi9Q31dHS2k5uTh1ar9VYrl8uJi4uf91hRlJArZMhFfwukmRSyW20TKJRmBntDfFYsIV1PT6uJiEw9I80mHzEDwgIY7fSPWY73WgiJn1EFQAYJmSqSU62IQiSCMxyQIZcrJ7+LHKVSdkpVkDTZ75UgMDCQtLR0UlPT6OzsYP+B/WRlriM6OnrFjqTIyEhkMhlHjh6mePPytufLXpfNwUMHiIqMQqPRYB41U19/gpCQULZs2XLSda1Op0MTEMDIyDAREXPJezLExcfT0dE+Z5vC4ZFhYmMX3hkvLjaOmtpqkpKSCAkJ4dJLL+Xll1/ic5/7/KLPvSTL6Xa7eenll9h93W6CgoL4w+9/z0UXXkRFRSUVFZWUl1dw7Fg5x46VU1tTx/j4BN/+zrfntGM2m5EkiSs/eiUXXXgRsbGx/PKXv+KHP/zhKfvQ3tHuq+J3MgQFBbFp4ybSUtOprKqgtq6WttZWv7QiaZ7E6ilog7w/tsPpISIiEFEUsFq60QUnz5lejjSPYG43M9Q04tPt2obtqLWqOdZsvNcbGohIUFN8cQDBEUoq33bT2yMyMOB9z0vK6XNMWKct73y1hvywwgiLTCYjJSWV0pIyBgb6OXrsKFbbygt0GQwGcrJzKa84tqy4n0qlYl3WOmpqq6mqqqSpqZGC/ELy8/IX5XDKzMikuaV5yecNCw1jYmJijqTv6quu5k9//tOCxwUFBSGKoi8mv3v39ex5/vklnXtJ5HzzrTdJS0v3ec76B/qJjY2lr68Ps9ns9wXUajUvPP8Cv/jF47z73rt+7bS3t6FWaXj77bd9r+m0OgICA5EkiUs/dAm3337bnF2KvfVCZQQHL5zCMxsREV4tqMUygdvjpr+/b1r4PrkNwkISuJnFusbHu1Eo9cjlGoK00967xoNzU4PsI/7OpJlC+4BgBZsvCyP/A6FU/ctC1T9HmXmvDgxML2L1+pN7U6emtLOt5mqVo1y/fgPpaelUVJTT1Ny0InEBeAmaEJ/A8RPHl3ysIAhMTEwwMjJCYGAQJcWlS7oPgoODUSpVfsKExUAmkxEdHTNHEHHllVfxj3/846QJGnGxcT4x/CUfvITm5iba29sXfe4lkfPll17iEx//BOBdhDscDnp7+0hKTiQjM50gbSCxcTGcf8H5fP7zn+PWz91CaWmZbwMY8KZ12R12tm3bxptvven7cn955S9kpKez9929k7set3P3PXf7nb+9o420eeSBp4LFYkEUJc7bvoMJi4VDhw74CbanQhazERYSgEqlYGxsGFEUUCqnqyoEaVXYJ9PIAkIDCAgNIHlLwhxizkRivo4rvpFKT72Vfz7Wg8Xsvdk1IRrCwgKIjQshNm7+cNCE1YVGo0Qm8+Z4yqeKla2xEH5qcFMoFOw/sI+xsbkSyqUgKSkZh8O+6PWXtwxpP/sP7EMml7NzxwcYNA4sa6BIS02la1IgsBTEx8XT29vr91pkZCTFxcX84x9/X/C42Ng4+vv7Aa/l33X1Ll56efElTBZNTkmSeOPNN9i4cQPgTWqOioomKyvLl33idLgoP1bBd7/7XQoKCtm0aTOv/OUVP2nd4OAAsbFx5Ofns33bdjZsXM8XvvB5Xn31Fb75zW/xm2ef5bOfuYVnn/kNL7zwPPv27QPA5XIyNja2qKyT2WhuaSI7Oxu1Wk1Odi4FhUW0tbVSVV3ps87KWY6dKYme2+3EaOwiNDSJkGANwTrvFGpkeO6IaRuxk7Rl/vVs5pYwNl4Ryd8e7qL9mNc6RmYbiMz2roEaKvqx21zYbS5qKnvnWE2Dfn554ulQCsnlctLT0tm0aTPVNdUrKgEpk8koyC/kRP3xU2Z/jI+Pc/jIIYxGI6UlZaSnpRMYGEhycgoty9DehofrGRsfWzKxdTodLpdrjoPzml3X8OBDD1FZWTnvcRqNBoVS4VsWXHzxxbzzztvzfnY+LJqc7e3t2Gw2dl2ziz17nqO7u5vYmFhkMhmJiUm0trYgk8mIi4vjgvMv4Atf+ALf+fZ35mhevTtTxSCXy3n22d/w4IMPMTo2yr739pOVlUVwcDA2u43e3l4CAwOpqakBoLOry1dYaSmwO+zY7Xb04XqfgFyn1bF5cwnxcfFUVJbTPDllm9m2QiGf9AY2ExOT6nPSgD8xk9Z7nQIzRQczCarRqcm7MIKiy6N47f52nNbJygYn2aohOd1Ae8cogiChUMhRKOSYxxa2yKcL2iAtxcUlHD9Rt6K9SwIDA0lOSl4w+8PlclJbW8Px43XkZOdSVFjkl6ObmJCI0Ti46B0DpiCTyYiOil7Wfi5xsbG+KeoUbrrpZi699FLKtpRy4403LHCcV4gAcN55Ozh48OCipY2LJufeve9w2WWX8cpfXuWpp5/m6LGjbNlaxtW7rqK1tYWrrrr6lG243C7cHo8fYT/y4Y+w57nnffHRa669locffohLP3QJ3/n2d/j85z+PKIr09fWREJ+w2O760NXZSXJSiq/mzkxERkZRVroFpUrJoUMH6B/wz6oYNPaiDdKi04UTETHd5/XFiYSGnjyhWxcdhC46iC03JJC5Tc/rD3Tgtk87cwzr9AzW+RfUTk43kJzu7010uQRckxpeQZDm3X9ltmB/LREYEEhJcSkNjQ0MDi5vUyDwTm8nxsf91oCiKNLW3sbBQ4eIiIigrGzLvFX+5XI5KckptLcvTpw+E/Hx8b6SIktBbFycj2TgrV370MMP8eSTT3DhhRdyy623zntcTEyML79Ur9eTlpq26L1VlkDOvezcsZOSkhJefullztt+Hn984UV2Xb2LP/z+uUWVJjQODhIdffJ6pdu3befu795DbU0dN910MzKZjP7+PqKjopZcslEQBAYGB4iNjZ3jJJneRVpOSnIqxcUlmM0mKiqOMjExjtVmYWCwj9TUTL8C0ROW6dE6NFRDaKiGjNL5B42NHzYQkaBh/wtDCO7p8xvWza81djk9uOapGA+gVip8XlyZTOaNxcomt7U/zRnYGo2GkpJSWltb/W7YpWAq+6TueC2CIGA0Gtm/fx8et5vt27YRFxd/0llSQkIiA4ODy5iiBk9OUZcWcw0MCEQmk2G3e2dNv/vd73j99df5+9/+wd/++nfO237evMepVCpUSqUvxrtz5072vrt3UedcVJxTkiT2vruXb33LGxYZHTWjn9wjMz9/8RXP+gf6yTnF/ogKhYLbbrvN77WOzg42bty06PNMoa+/j5gY76a1U2ScEr5PYeoGUKs1FOQXYDaP0dB4ApvNSm52AXa798cP0Cj9iDkFu8P7fmSO1+IN1XsdHZuviECjU3DoTyOIAiRvjUetW9jln5o7vZaeSVBd0PQxoij57QvqcjlxOBy+yu9OpwOHw0FXdxeBAYEEBAYQGBCIUqlcddG7WqWmpKSUo8eOIgieZe2XqtVqMUQY2PvuXsLCwti8uXjRel+5XE5MdDQDgwPELxC3XgheK9i/5C0+4uLi6e3rIyM9A5vdxtatWxd1/4eHh2Mym4iNiWXnzp38+olf8/Wvff2Uxy3Kcra3e0s3ZGZ6NYajY2OEhS2tDL7H48Fhd5y0ktl8sNqsqJSqU6b5zIYkSXR2dviURFMCCW9OqWzeUpOiKBEcHExMTCyhoWG0tDUxMNDlk29lZ01POVOzonzEnInIHAPrL9MTEh9MT5caSfJmpgC45iH3eN+EHzFnYqbFBm+cubevh8qqcg4fPUB1TTWdXZ2Mjo0ik3kzIhQKBaIoYjKN0NrayrHyY7z73l4qKssZGOhfcThkJpRKJSXFxfT1989R75wKbrebEyeOMzIygkKhICM9fclC/PiEhGVVY3/PoAAAIABJREFU3IuPi6evv/fUH5wFg8HAv/71Fvd+714cDseii9CF6/WYJ51oS1l3Lqr1vXvfYefOnb7R12KZWPJ+FyMjwxgil+5pHRgYICZmYSXGQjCbTWiDtCdVo0yJ2P0h0tvbzfrCYkCip7eLhoYK4uJSCNKGERI8PX1PSAqbUztIGywRZpBT/paXiDFF0Yz3zk2u9rWxMXbe1/PzvNN/p8vF6Oiwr+i0wRBFTnYOQUFaPw/z1NS2ubmZlOQUP0spSRLj4+P09ffR1NxEsC6Y2Lg4oiKjllWIeSYUCiWbNm7iwMH9ROj107LIBSBJEt3dXbR3tJOWmkZOTi4Wi4Xauhq2lG1dkoXXaXUIkyG9paiONBoNCoUCq9W6pGyVl19+CZvNRmNjA6+99hp33HHnoo4LDwufzILyX3dO7cK3EBb1y0ytNwFqampwudxLrrw2OjpK+BKtLXhDL9HRSydne0c7KYuIiU5Z0anQyeDgIIaISDQaFXK5gqTEVDIyCjCbh2hvP47b7fXKuj0i0bEhJCSFkZDktYwyuURaroyuDjWhidMx0fjiudOuhI2xPmKqVHKy8qa/Y1xcCCMmGyMjAzQ0VuJ0OYiPT6OsdCsZ6RkETYq65/X3zlNoWiaTERoaSk52Dudt30FqahqmkRH2H9i34h2jwWtB83LzqKldWJwOMDJ5TqvVytYt20hMTEIm84pKdFrdKXecmw/x8fH09S3dCsbHJ9C7hOMEQeAnP/0Jubl5/OTHP+Xhhx7mggsuWNSxKpUKUZR8s5bFrjtPaTlnrzd/8pMfE2GI4Nprr6GoaD2bNv1/6t47vM366v9/adqSLe+993Yc24lnBpBACRD2DC0QmoTRByijwMMo0JbR0gYo0LLKaqEFQkKfsiGMJHY8kjjDsR3vHW/L8pC1f3/Ili1bknWb9nv1974ucxFb99R97nM+57zP++SzKj+fdevWExnpPPZXj6mJjY1z+LePP/kYqVTKuT861+731kW0SPAcDK1Wy/T0tKCXgXSmkbm7u5OMdGunga2fUu5BdlYOA4NDdHa1MDEhR+kVZldeWb0+nuGhDoaHROh1c+881QypIChtLiT2CbF/W8enWDPVUVG+mM0W9PophoY68FGpSE1ZiVRqZSRZWMw5WPS7JcqeIpEIPz8//Pz80GqnOFlXR0dnJxnpGT+ovzMwMIjTfX3WpoQF3/PU1BT1DfWYzSanSv6JiUmcqD2+ZMJwIcLDI6iqrnSo7u8KYaFhlB8sIzkp2dYU8c2339BQ32BTMnji8Sdsoeunn36CTqdj9apVdHV1sWOH4+ysM/j5+qIeUxMYEOj2unNJz7lwvRkaGsqWa66lpbmVhx56iJCQEJ5/4QW2bf+p031YLBaXoccjv/wl11xzNd99/x1f7/2awsIC8vJyefbZ5amXdXS0E+fkReDw/GZC2zHNGHK5HJVqLmQ3mSx4e1mTMt7evqSn5SKVKVGrm9FODTJnDVpCwjzw9bGG7jExPmSXziVJPGbC4YgEf7y95YQEW8O/4jVz3t1sNjE83MHQUDuBgTHExCQjlcrwUlp/YPH0bRz92036nkKhZFX+KqKjo6k+VPWDZ2imp6XR0dlhY30ZjUZONZ7i8JFDxETHsHpVgdMw0svLC4vFYsuGugu5XI7CU4FGMyZoO6lUikqlQj3DeHrnnb+xffs26upOEhgYSHl5OS+9NDeasrR0DYmJSWy5dgvDyxjeu5x155Ke88SJ45hMJi6/4jLy81cREhqCUqkkICCAjRs2snHDRhrq611mU7VarVMBroaGBgYGB9i160MuvvgiQkNDeerJp4iIiKT8YBn/8z8/o6mx2e21iMVioX9ggBQXMpkwF/rNf5A75r31ZTIJhgWN0UqFjKGRKWJjojGZIujoaEUkHsRi9kMiVSOVRgHTrMgNRz1sZYXMKhwABIaqaG+Z+2JnDVPlLcdkMtLbW49KFUxKchoikQjlvNmf8+1tYZ+pO32crhASHEJgQCBNzU0cOXKY3Ny8ZU0ak0ikZGVmc+z4MaKjomltayEmOobSkjVurW2joqLp6u4mJTlF0HGjZhJDGRnCJp+HhoYxODDA6Mgo9/ziHj7/7Avb+L7NF2zmrA1ncuWVVxISEkJAQACffvIpD//yYWpqalAqvJZcM85HgH8Atb211v8PCCA4OJi2tjZSUpxf65J3rLGpifPPO5+rrryKsbExIiIiiI6OtvuMVCqls6vTyR6sIa2jaV8Wi4VXXn2Fyy+7nA1nbaC66hDHj53gkksuJT09HQ8PD5599llBSQKNZgyVSuV2omPWaxoMBjSaMQIDAm1/k8klyD0kNu3Y+ZBIJCQkJBMbm45YMjJjPBa8fazRgV+gYw8Rl2jdf2LqXPhmMhkYHm4iMjLetg5bCJFocQP4QnFskUjkcM3pDiQSCWmpaQQFBXPocPWy+znFYjGTkxN0dnZQVFgsaKxBeFg4fX2nBRP3g4NDGBwaFOz1AwMCGR4e5rrrr+O+e++3m6uZkZHBtdf+mIfmdVVJpVKefOJJkhKT+PVvfsUtt9zMnj27HUq2LoRCoUCr1dquLTk5mcaZJJEzLHnXmhobWb16NVdeeRWPPvII6ekZiwS7nnjiSd599x0qKysd7mNMrcZ3gXGePn2aSy69mO++/ZY77vi57YQ9PDzQarVcf8N1qLxVXHD+5qVO0Q59/f2ECUggzWY5u7q7iIqMsoaI85qZZyGRiJBIRATNcF7lMglymXWEn1LpjUQSgtHYi4/vJLD4IQkO97F5uLPOnntbms16BgYaiYtLxtfXel8X8nxlcikyudQ2AQ1mvKVlblQD/Hu6UeLi4oiMiKSqqkqQgt709DTHjh3l1KkGVuWvxmQ2CS7bSKVSfH18GREYNorFYgKdNEW7goeHB13dXfj4qLjjjjsW/f2XD/+Szz7/jKqqKrvfl5aW8tBDD5OSksprf3mNuPhY0jPSSEiMJzQshPiEuEXXLhKJ8PT0sFEOU5JTaGp03cLmludMnnG9o2o1fn6LPWBwcDDP7HyG7Tu2OeQ7qsfU+PnNhRzffvct+avyyMlZSUVFpV0xWK1Ws+m8TcTFxrFihXsTh+djaHCQYKElG5FVvzYywrGekWSBsQQFKPBRyRnT6NBo+vD3jyAuLpK0tDykUjlyjz4Cg0xkrYwgtyCG4PC5TpOYWH/btDKTScfUVAfJyemoVPb3VSoVozeYkcnd74f/dxENoqKiiYuLo6q6ckn+qlXPqZmq6kpCQ8MoKCjEz8+P+LiEZXWAREdH09UtfEJXVGQUPQLZSu3t7VRVVfHcs885/N59fX154vEnuP2O2+28sp+fP2aTiTvvvJNPPv6UvtP9fLhrN9/s/ZbaEycxm80O668KT4Wt0SI5OcVWXnGGpT1nU6NtDTCmdhyeAmzefCHt7e0zPZdzsFgsGAwG5HJrQkSv13PzzTfxyiuv8tijj9kpGvT19bFhw1mszMnhtttuw99f2Lh5s9mM2Wx2W45x1tMMDg7i5++PbOZcxGKxW6FYaIgHBsM0/v5+KJVWFo6XVzDBwSno9VMMDzeim54gKtqPqGj7+xYUpMBo7CEjPYuggAAUCwgHCk8pCk/pQlUSTPNYQvCf60oJD48gIT6R4yeOOR6daLHQ19dnbeUCSkvWEBYWZntBhIeH09ffJzjU9PcPQKPRCJ5X4uvri0ajERQ97N379aKGh4W49tofI5VKeeutN22/k0qlmMxzntHT05OMjAzi4uIIDg62zU1ZCGtoa014Jack09jk2nO6fC2PjY0xNjbGwYPlnHPOj1Cr1SQmJjn87CeffExxcTGBgYF2v19Y6H311VcBq27shx/uYmBgkLb2Ntrb26ioqGDHjpt48IEHOXbs6KIZm0thYmIcb5X75IjZL6W7p5uE+MWp+PkEBYlUbEerE4lF9J62kvH9fBU2jSFfX0/Gxqbx84vGYtEzOtqJeFxKTEwimZmhiMUiRtXTjIz04O0diJfXHGNK4SFFqzOi8Jz7Wubr1c5/hEQL1p7/CYSHhzMw0G+NKuY1HWjGNdTX1eHp6UnB6kKHWXiJREJgYBCDg4OCyiMikYjwsHB6T/cSI4ASKBKJUCoVTE1NuU0sOP/8C/jHe/9g957dHDlyhJDgEEJCQwkNCSEkJISQkFBCQ0N54okn2LLlGi655FJb5CiRSDAajQ5ZQslJSezZvZv169bbveRn150wE9Yu4TldGmdTUxNhYWG8/sbrbNu+jWeeeYaKikqSU5KJiY4hNjYWsVhMRWUFV19jbcKur68nPX2OP6ses19vzpZlXn/jdaRSKcFBwcTFx7MqP58HH3iInJwcwKpBJKTTHWBsTIOvj7CMndlsZmJiAh+fxU3Oiyh+8z2WxcLgYB+5OflIJSJ85s3SnP3/qWkDQUF+aDQjtLfX4eMTSGhoFHKZnulpDSuyF2e4/f0UtuOaHYx9mD+Y9/8FMjIyOVhRTmBgEGKxmMbGU2jGNWSkZzpc4sxHTHQMpxobBNcuo6KiqTlaI8g4wVprHR4edts4w8LC+Pu7/+BA2X7GNeO0trZiAbo6O+kf6Ke/f4DBwQH6+/vRarW88847/Oxn1vEjCoUSrVbr8Bl97LFfcfU1V3PhRZt56823bQ5LoVAwNmYt+czOZJmcnMTH1/Fz7tI4m5ubWL16NX9/9x9oxjV8//13HDl8hL3ffE1nZycdHR2YzWby8vL5xS/uJSM93VYPncWYWm3H8LnjjjscLr7nw2g0Iha5F1raHUszJpjqNz4+7vZI9dnyilQqYaB/wCZ9YbGAh4cEnc4+CTBbm4QAvL39GBk5TWNjDRaLmezMPDw9ZOjmqct7zsxmmSW4iyViGznC2emZzWYmJyeZntai01kbggeHBlF4euLl5f2DvapMJiMtNZ2KygpEIkhKTCIzM8ut/fr4+KDX6wXT66xkCAt6vV6Q3q2vj8+S4/+OHTvGc889S3NLCx0d7QwODrJz504qKiuRSiWMjI7w2aef221jsViYnJy0I2koFQqmtFMOjTM8PJyvv/qaBx96kIJCq/0UFBSgUCiZmvGcsw3szc3N5OXnOjxXl8ZZX19PW1sbJ06cICk5idDQMJ5//gW7kzabzbauj3/96/+oqalh1apVti9vTDNmV3M8duwYT/32SZoamxgcGmRsbIzQ0FCioqLJyszkuef+iEajcejJloJGM+ZyDKDDbcY1M8Y5e022i0MsEmGe+YVVSFqEXC7FbLbQN9BDbIx7kileShmTUwbiYuMxm/VMTo7T3NJAQnyKrRHAkUecNcyFsFgsDA0N0Xu6l7ExNT4+Pig8Fcg9PGY8+iBa7RQTExMEBQYRGRm1pJdzhsGhQU41NiCRiImJjrULb91BdFQMXV2dJAusXVq94BDh4e7PwVEqvZh0oumjVqt59NFHeP+D97n3F/ex9cYbiYuNIyIigsrKCrb9dDsmk4mU1GQOHz5Mfv5cVCMSiRZxyRVK5SKNq/mQyWT87re/o6S4mIsuvpDHf/M4N9yw1W4IV3JKMk1NjcszzsamRmJjY9mw8Sw+/+wLPBfQ6EQiERKJhN7eXnbctJ2e7h6mddMYjUauvPIqHvnlIxiNRluC5unfP82zzz7Dfffezz13/4Lg4GB8fHzo7+/niy++4M8v/QmwrnUdNdm6gsVisTuWuxjXaAieJ4S9sKVs/rXOQq/XMa3V4uPji0QsZnraWnKY7z1F4jlJTbCWXkwmExqNmvzcIiYnJ2hta0TppSQxPtmmTwtWVpLZbIH59jpzYuPj45ysO4lCoSAmOgb/FTl259bT001GunXUg9lsZmhokKamRiQzk5/d9WCTk5PU1Z9EIpaQn7cKqVRCVXUVMTGO67DOEBERQfnBMpJmaHLuwsfHh3GNhnDHfQEO4enpiW7aXjHCbDbzt7/9lQcefIALN1/IieO1i/IiCqWSqakpfHx8+Pkdd/KHnb/n3XdcK+UpFQq3dJDOP/8CcnJe4siRI9x440+xWOaWKsnJKS6TQi6Ns7Ozi98+9VvWrV3Pn/78J35xzz2LPnPy5EnOPmcjN910M3t2P4BUKqW6upq169Zwz933IELEyZMnaW1tob29jczMLG6//Xa7ffj5+fHCC89z4YUXATA2phacDJqamnI5BtAZNBqN4yTXjHKCBBEWswXTPA86MNBPeHgEErEYCxaHbCIPB55wcKifkGBrUsjf35c8v1UMDPZz9NghIiIiiY62ruFnSQ8iscguG9vW3kZf/2kyM7Lc8oRisZiQkFBCQkLpH+inqrqS5KQUwl088QaDgebmJoZHRkhPT7cjZSiVXqjHhDUwSKVS/Hz9GBkZFqQZ6+XlRd+MOJa7mM/6EolEHD16lNtuvw2DXs+e3R/ZZr8uhFKpZEprNU6lUuFWpnh2zbkUnnzqScRiMc888yxgn0hKiI93yg2AJUopGs0Yfn5+3HzzzUglEo6fOL7oM395/S/s2HETj/zyEWQy2UxfpJmszCy8vJQYTUbWrC3l+Rde4I033qC9vY3duz+0pdgbGhq4/IrL+PKrL7n3F/cCMD4xITwZpBkTnAyyWCzo9DrbBC4rAdrqpGxh7oxxSMQiJBIxMpmEyclxfOYdSywRIZNJkMkkeHvL8faWI5VKkM5MEZv1xH19PYSFRiCdIdRbNW3CKFhdjNFoorLqIIMzhXRbhlZsZf00Njag0YxRXFSyrBA1NCSU4qIS2tvb6HZQR7RYLHR2dVJ+sAyllxelJaV2hgkQExNDZ6dzJpgz+AcE2BIh7sJL6bUsvVwPT08mJia497572XTeuVz3k59QVlbu1DABlAqr57RYLLz66qts27Z9yeMoFArb+tEVxGIx2dkrbBGdQqFkeno2s+/LmAtOsEvjnJiYQCQW8eRTT2IymxhT25PQTSYTu3bNyWXOYv+B/axdu5ZpnY6hwSE2X7CZZ3Y+Q3x8PM8990d++cgvCQsPZcPGDZxx5noKC4s4WnPMFm5YLGZbMqipqYk//vGP3HXXnVxy6cWsWbuGx594nK4u+wdMMzZmZzDuQKfT4elhH+Y5jLzmW+vMffHx8XFLZX3WQKd1U8jkcry8vay7E88qMEiRyWUkJiSRuzKPvtO9HKk5xMTEhK1m19LSjF6vt853+QH9lzKZjIKCAnp6e+idJ1Y1MjJCeXkZExMTlBSXEutESC0wIJCxMbXg2Zve3t6MTwibsznbNC6Y9WSxcPMtN1NfX8eJ47Vs375jSZ6wUqlAOzVFdXU1mnENGzdsdPP8lmZAlZSUUF5ebvv3fKP2VqmYcHFfXIa1w8PDnHHGes7bdB7bt28nK8t+aMu+/fsICQ4hLW1uhuGuXR/w7LPP8NZbb6PT6WjvaGf9+vW0trUSFxfHpnM3sencTfT29lJdXcW6dettA3hhJiRBRG1tLU88+TjffPMNl15yKUnJyaxduxYfX1/27N5N/qo8CgoKuP32Ozh749mMacbcUoKfD824BtWCxNPC0Gj29TX/GTHMrG1FIpFN2FkiZlFoC1avKkVCT1sXkRHWZMqsV5xdU4qwTgtTKr3IyclFrR6ltvY4vr5+hIaGMjg0SHFRsa216YdkYCUSKXm5+RysKEepUNLa1orRaCQnZ+WSDfQikciq4drbI6jrx9vLm0mBxgnWjhODQW8jsCyFAwcO8MqrL1NSUswtN9/q9otsNqx95ZWX+elPt7m9nURszSO4Mv6C1QUcm5nG7unpaccS8vb2Znzc+X1xehYikUik1+uprKjirbfeRu7hYeur7Onp4e577ubKK6+wlUUMBgNbrr2GRx59hA937Wbjho3opnV0dXXz6Wef8rOf3cqF86ZbR0REcNFFF+Pv749er2dycnKGTaRHJpdx9jkbyclZSVNjM3/605+56867uOSSS9lw1gZeeOFF2ts6uOzSy7j77rsoLi5iXDMuuO9To9E4nem5uGHZ+qPT6fCQy10aiFgisv2AlZc7NjZKQECQnYC12Mk+AgICKCoqwVulovpQFcEzw3kWEt2XC7FYhErlQ0XlQaIioyhYXeC2skVYWDhDg8I4rDKZDIPRKNgLeim9mJhwL7Q1m81cfMlFrFu7ng1nbRAUYXh6Kujo6ODb777lxq03ur3dfFKBM3h7e5OWmkZFZcWibVTeKpcvLVdX4CESiXj5lZe55ZabaWio5/bbb+PKq65gZa6VKHC05pht0tKpU6eoqqriUPWc/IJOp+PGrVtReat48823uOmmm+0O0N/fzyOPPkJsXAyhYSGofLzZsHEDPd09DA8Pc8/d9zhdeyqVSrZuvZFjR49TVFREV3cXW7feIGjE9/gySjYTk+N4z5zTQurc7LpzIcxmEyKRGA8P67rDkZHN/muWgCAWi5DLZMTGxGIymSg/WMbwsLApVQthsVjo7e2hrLwMX19fAvwDBDdYz3oZofDwkAvWmfXycn/dKRaLeWbnM7z19ptudYnMwmKx8OhjjzI8PMy333wniJetcPNe/OQn13HxxReRl5fLybqTtoyySqVifMK5hI0r41R5enoSEBDAypXWkKekpISLLryIupP1/OH3f7BTPpicnCQoKNjuyzaajISHh/Pmm29x1pn2kg47n9lJZlYGQ4ODfP/dPibGJxnoH+Txx5+gpqYGT09Pt3oKxWIxd9zxc+RyOUPDQ3z33XdLbjOLiclJwVpIE+PjqLy9EYnsZ5TMtzWJWGxbawJotdOLWCviGUnLhdKW85e3XV2dxMbFkZ6eQW5uHm1tbRw+fMjlfA5nGBsbo6LiIMPDwxQVFpEQn0BsbCydXcLI6dYXi3hZkpQTLh5ER/ASGA7/5CfXUVBQyLHjx2hoaHD6OY1Gw2eff8YDDz5ASWkJn3z8Menp6YSFCWMyyeVyDG68cG677TYGB4a44sor+WjPHlsEYQ1rnd8TV2tO74CAAJuUQvnBMi7cfJFTxbEp7eJShm2exwLo9XqeeupJKg5WkpQ0V8ZQKpUz6fsAQVOAIyMjSUxM5LxN51MnYEjO/MST478vPvfxiQm7l5JdUmjBx6VSCVgsTE2Oo/JWuWyKXphcmp6exmQ24+1lfXl4Kb1YtWo1Q0NDHD5ymJDgYBITk5ZUgNPpdDScakCrnSIzM8suUggODqHhVIPgdaz3jEdztiRwBJW3NxMTE4LLKd09wjpUnv7d0+z5aDdnnnUGl192OQ8//EubYPmrr77Ca6+9RsOpBlatWs3atWv5zW9+Q2lJKbUnazEYjW6vb2Fu5qs7kMlkXHThRdx08w4sMw+K6gckhFQqOxnLuROxNiZrbNnVgwcPcv999y1KV4vmMWzm46uvvyItLd3OMGchnmkYnuXYLgWz2UxFZQXt7W08/sRvuO66621/s2q76p2uRUVLTAFypJYwMTGByom85yKW0QzGnXB352+zEOPjGrtE2SyCgoIoLSmls6uTsvIyEhMTiXQgwGw2m2lrb6Onp5vk5BTCQsMWfUYsFqNUKJmenhYU3nqrVEyMTwgyTm9vb04vQa1bCOv6TNgYCrFYTEhwCLUnTvL4E4+TvSKLu+66m/i4OB5/4nHefvuvFBYULnomxAtqyu5AiHECpKWlMTY2Zot8FAoFer0ekUgktVgsi4qrrsJab695Id/sF6vX69l49kbi4mMJjwijsLCAq6+5ittvv50//vH5RSfviG6z+8MPueLyKxweVCRyPpJvIfYf2E9CYjx33vlz5HIPvv3mOx584EEeevghoqIjUXop8PJW8pvHf+PW/pxh/hpxKW9r/by90U1OWcPnhYboylmNT0w4DblnxxEUFxWjVqs5WFHOqHp05vws9Pf3U1Z2AIvZTGnJGttMG0ewljmEhZveXl5MTArLvspkcsElGLFY7Fa5whECAwPZ+YedlB0o59CharbeuJX33/uAdWvXOXxZi0RizBZh7W1iJ87H6efFYlauXMnQ4NDMMW20QIdftGvPqZpvnNa38W23/Q9BQYGMqTX09fXR3NJMXm6ewwfJmecUi8XIZI4PLRaL7ShOzjA0NMSPf3wtz//xeTZvvpCv937Nno/28Kc/vcgZ68/g22++IyYmhrGxMdatX0tkRARbF2TiLAInzYpEIqRSKUaj0a31sEhk/c/8rOx8G1loMHOq9CKmp7VLUhjlcjlZmVmMj49TV38SkUjE5OQkp/t6Wb3avQnS7jJd7M4T55lmZzAajcjcFGGej2Vlp0VzpbCkpCQ+eH8XEy5edjCzBFuG5xRinACrVxfQ3tFuYwmpVCo0Go0KWJTFcuUCxPNvjMFgYMdN26mpqeGtN99GLBYTERGB2Wzmf//3fq798RbOO38Tl19xme3Ldub28/LyOXLkiNMLXuomWSwWtu/YxhVXXImHpyd5ebk0NTXS1dnJrg8+5K233rZJnoSEhPDnP7/EKzN9pI72JQRSibXRdr6C/FIljlmDXgrz9yORSDG7mXTx8PSwJk8mp6xvYy9vtznGJrPrOp0jOOtjdAWD0YBUIO/ZYDAgkwrbBnAoD7pU4k8sFu45hYa1AD++9seIRCKuu/4nGI3G2e/b4cPjyjjH59eYDh48SFRUNPv3H7C70Ddefx2tVssF51/AbbfdjsFgYOczO2dOXuzw5NesWcPuPbt57FePLRrH5s5Neumll+jp7iE1NZUbb9zKE08+SVZWFi+88KJDmlZRYRG9vT2LhHzFYongTn2pVIrR4NjQnBmoVCbDZDIKqlO6Y9Bms5n2jnYOHjyIn68fZ6w/A6VSiUgkpqz8AKfdEMtajkdbjnEaDcKPozfobeoUQmBBOFHDvAxyx3KM09PTk8KCQkZHR9m69YbZhJDDNYIr45yYTWOPj4+j1Wp56MEHF4VK/QP9XHHllVxzzRY2nbuJZ595jueee5aenh6MJiMd7e288cbr3Hnnz7nxxq1MTU2RnZ3NX9/+G0drasjITGf79m2XVyBiAAAgAElEQVTU1tbOXbADzzkxMcFnn3/Gffffx6OPPcK5mzbx298+xTd7v2XTuZtc3hCFQsELz7/Ajh3baWlpsf1eIpHYyU24A4lUitGFMt1CNTyRSIRMKsUwz9DcMVKFpyeTk85rfENDQ5QfLGN6eprSklKioqJs+01MTKSwoIiBgQEqqypcav9OTU7iKXAOjdU4hXk0o9Eg3NvqDcgFetvlipzpdTo8BGRqYUn9bscQgUgsZveHe+gf6J/lHDs0Tld3a3w2UdDT0zOziF78QPX32avdxcfHs2PHTeTl57J69WoiIyKZnp5mRU4OVdXV/PWvb7Njx008+ugjDA4N4uvrS2tbKz869xyys7N57NFf2QS1Jicn2fnMTr744guOHTtKbGws/v4BvPH6m/zl9df41WO/IjU1deaaXT/smzdfSF19PaVrSigpKeGuO+9G7iGzqg0I+P7d8WgLDc/dsHY+goKCOXLk8KI+yKmpKev6EhF5uflOO3E8PDzIWZHD2Jiakydr8VapSE1JtWtetlgstn5WITAalmFoRqNgwoPBoBfcAmiNBISHwrMNEEKg1+kEv9hE1gUxCoWCf/z9PSKjItA7YWe49JyzBdLmlmYUCsWit9IsdzY2Ntbu9489+hjfffs9z//xee6//37eeuttLtx8IS0tzWzceDb//OdHmM1mmpta2PXBh0RGRJKclMyll17KDVuvp62tjaamJkpKinnvvX9w/PgxQkJCyMzMZHBwAL1ex4oVOdTV1wu6Mffdex8tza0UFRVzxZWXWz3nMuQbTQINbVbuU+g2JrPJluE0Go00NNRzpOYwcXHx5OevcqtFztfXj6KiYgICAqioPEhbW6stlB/TjOGjck8FYj6001rBD7LRYBAe1uoNgpQQrNvokXsID4XN5qWz8Auh0y3DoPVzobrJZHLZfbXEmtPqbWtra/Hx8VlknPsP7CcjY7GWjEQiIT09HZlcjsViHeDy020/5cEHHuLIkcPc8fM7eOxXv0IsFpOXl8ebb76FykdFfV09e3Z/RHd3N/mr8tDpdMTExNBQf4qW5lbe+8f7pKamIpFIyFmxguPHjwm6MWD1xlVVlaxbt85GXBYCqUS4F5zt6heKqMgo2tpa6e7uoqy8DIVSSWnJGoIEFPJhjrBeWlKKwWCgrPwAA4MDtLa2CJ6rabFYlsx8OoJhGaGw1XMKNU6dICIBLD8Unl6GcU7rpm2iBRNLtEa6Mk6tXq+3vrHr662zLBZE2V999RXnnHOO853PLJhffPFFTpw4zt/+9leefvppXn/9Dc7bdN7c58Ri/vr23/jyqy+56OILsQrwevLjn/yEj//1CRERc1IVRqMRiVTKihU5HD++uL/UFUZGRshekcU///lPjEYjZeVl/P3vfxdkbB6eHminhRXG5XI5FrNFcJ1PpfKhpbWFkdERSopLnLZyuQuJREpKSiqr8lfT1tbK0NAQngr3tX1gVrDNV/B5WI1GmHHqDQbB2+gE6g6B1YMJDdNhxnN6CjNO3bQOj5m8zfj4uFNCC7hYc1osFouvry/j4+Ns3ryZ+oZ6Npy10UYnAytbpb2tzenORSIRw8PDPPjQA4SFhXHXXXdxxRVXOgwf/Pz8qDlylL6+Pk6dauDA/rJFcyQ0Gg0nTpwgMiISHx8fu0SHK6nCWfj7+7Nn90cMjwwzrtFgwUJZWTk/OvdHfLTnI7cavAMDAmlrawWEaRUFBVkVyd3RxNFOa2loqMegN5CSksLYmGZZD48zyOVyTEYTaalpHDt6FP+AAJKTkt1a3w0ODhIcFLLk5+Zj1tt6eQn0tvrleE7hxqlbRjIIZjy7wGhAp5u2jXAcHx/Hy1Xt1dWOZrl/l112OR4eHrz6mn2t8NJLLuWjf37kNDT08PCwllb+sJOTtXVcddXVLuN6uVxOTEwMHp4eDge8VFVXoVar2bdvH4ePHCY7e04RXqlUMrVEB4NIJKKkpITNF2xmy5ZrOevMs7j99tuJi43llltvdiu8Wa4XDA4OYXCJViuTyURTUyOHqquJiIikoKCQxIQkPORymppdCxC7C4vFwrHjx4iMiiImJpaSklK8vb0pP1hOZ2fnkvdgaGiQ4GBhYfVsa55wbyvc0PR6PR7LME65wPB0FkKvaXqet52YnLCbaLcQLo1zPmv+7I1nU15eRkfHXBdDYmLizJzDcofbKxRK5HI527fvEHSTxWLHXQ8bN2xk/74DfP7FZ1x66SWsyJ5r/lZ6eTE5Kaxbw98/gNHRUV544UUOHz7Mvv373NouMCiIIYFrSD8/P0bVo07V03tP91JWfgCJREJp6RpCQ+Y6JDIyMhlTq2lpaflB81BMJhO1tSfw9PAgNsaaxBOJRMREx1BSXMLk5ARl5WVOR9zp9XosZovgNd3g4ADBwcK8LVhpj0J1ofQ64WtO/TIytSaTaVmqFNPT0zb1jYlx1xPil/Ccc6FjSEgoxcXFfPrpp3afueTSS9n94YcOt/f09GRaJ2x9Bq6bWHNycvjk40/57LPPue22OaEwL6VSsOaMp6cnRoM1I3jDDVt5/7333NouOChYcMPxLCG7u8d+hoZGM0ZFZQVDg4MUFhSRkJC46EsXiUTk569ielrLocOHBHttsJZgKioP4uXtTfqMOt98yGQy0tMzWLlyJa2tLRw5cnjRrMyWlmZiFmTm3cHAMubXGAyGZWkXL2fNOdtALwTLMWjrsaZt22k0GpdrTpdXHhsbYxtGo1AoSEtL44sv7QV3L7n4Ev75f/90uP0sg0Lo216pWLqJdd3adXahr5+/P6OjI4KOA9Z16OjoKBduvpA9H+1x61z9/f1tRHMhSExMoq2tFZPJhE6n4/iJ45ysO0lGegYrVuS4/LLFYjGZmVlERUZSfrCcU42n0Ot1Sx5zcmqSE7UnOHzkEBnpGSTEJ7gMxby9vFm9qoDo6GgOHT7EqcZTmExGtNNahoaHiBKoW2swGDCbzYIf5MGhQYKChIXPAJOTE24rvs9C56JzyRmsmVphyTSwTz51dFr7dZ3BZZZhvq6mQqEgODiYffv22dV3XEl9wEw3gtGAXMDCXjmjIyoE8xXUhKwDgmZC1A92fUBUZJRb24rFYhQKxaI5MEtBLpcTFhbOkZrDaLVakpKSyc7KFnS+4eERhIaG0dPTTVV1FVKJFH9/f3x8fZGIJRgMBnp7e1Cr1ajH1IgQEZ+QQJabKu2zCA4OITAwiI6ODsrKy5DJ5CQlCtOeBRgaHlqWkQ0M9C8aX78UZgXBhHKFdTodcoHlqeXUOBe++JuaGjn/vPOdft6l50xOSqKp0TpsxWNGTTwoKIj29nbbZz7/4nN+dO65TvehVCjQTgkrwPv4+jImQGoCrF56YQbXHQQEBjI8PMTNN92MZlzDK6+87NZ2wUHBSyZ4FmJgoJ++vtOoR9UUFhQRER6xrNKIWCwmOjqGNaVryc3Nw9fPj3GNhpGRYUwmE5NTUwQHh7B6VQHFxSUOezndPU58fDzZWSuYmBinta1VkAQIWNebIQJDWovFwtjYmNOJds6wVGnCGZZa+zmCVblRmHEajPZE/qbGJpKSkp1+3rVxzvOcsyGqwWCw49d+8fnnnOvCON3VWZkPXx/Xep7OEBgQ5DSZ4QxymRyLxRqqfvLxpzz2q8c4ceLEktsFh4TQ1+ee6PHExARV1ZX09PSwenUBySkpM+WYHw4PDw/CQsNISUklLS0dT09PkpOSCQ4OFkx9cwSLxUJzSzO5K/PIzs6m4VS9TU1uKZjNZkZHR/ETIEIN1tEJfr5+gl8oIyMj+DloUHcFi8VijeyWs04VaJzza5wWi4XGpsZFs4Xmw6VxpqRYx5TNuuM5YVwrn3BgYIDGpkZKS0qd7mOmX03QRYjFYuRyuVsPwHwEBlrHiAtFYEAAIyMjJCYmUlJSQmPjqSW38VJ6IZPJXB7PYDBwsu4kx44fJTkpmdzcPBSe1jEKYxqNnXbsfyuam5tQKBQEBQXho/KhsKCI0NAwqqoraWlpdsmw6u7pJiQkVHBSZ2Cg3yYtIgT9A312WW53sJyMMFjHTSqFrm3nsYP6+vpQKBQO1S5m4fKuBQYGIhKJbDMh2lrbZsjn1h1WH6omJCRkkQ6K0Wjk7bff4vIrLmPXB7sEjxEHCAgIFLzdbK1TaAIqMGiOXhc0M1PSHaSmpnKqsWHR8SwWCx2dHZQfLMfHx4eS4lK7QcBisZi83DxaWloEh4n/L9Hb28uoetQ2ewWsEVRYWBilJWuwAGXlB+jr71t0D0wmE21tbSQ5mefqCtZkkPDsrtFoEkyut87lEa6g70quxhmm53lOq9d0PdzJpXGKRCJraNvYyM5ndrL/wH7++NzzNg+6ccNGztt0Hnn5uXz73be27RoaGrjl1ls4/7zz+WDXB9TV1dutU91BYEAAlVWVfPHlF25vIxKJ8FH5MD4ucN3pH2ALhyMiI+nu6XFrO29vFV5KLwYGBmy/Gxoeoqz8ANqpKUpLSomOinYYnslkMvLz8jl+4tiiksV/A0ZHR2ltayEvN9+h55NIJCQlJlGwupD+vj6qqivtlOTa29uJiowUHFprtVqkUqng7QYGBpblbZczNEtv0COVSgWH3ZOTEzYv3dTYSIqLkBbcGDsfGxvL/9z2M1555WXuuP3ndjKOHh4e7Nz5DC+//ArXX38d//vA/6LX68nMzCQwMJCSklK+/eZb5HIZF19yEf4BfkRFR1JUXERVVZXTY3Z1dTEyMsrIyAh33XWnIHJ6wDJCW6lUiqeHJ+Pj48TGxNAh4EWSkpJCU3Mjk5OTHD58iLa2NnJz80hLS1+ScqdUKsnOyubwkcPLql3+pzA1NcWJ2uPk5+UveQ2enp7k5KwkNSWNE7XHqa09weTkJN093YIV+GE2pBUWmoI1pA0TOKQXrOtUIYOZYPneVj2mtjWJNDY1keyABTcfSxpnVmYmCfEJVFcdIjk52WEY9qNzfsThQ0doaKhnzZpS9n6zl9jYOD7+5GMkEglnrD+TXR/soqO9k6rKau78+c+5+JKLrDNY5hmewWDgd0//jlWr81m3fi3Dw8N4e3vzz39+5PYNWO66Mzo6ho6OdlJSUjl1auk15yxkMjlYoKKygtjYWFavWo2X0v21iL9/AElJyVRUHBS8Nv9PYHBwgEOHq1mxIgeFwv21mJ+fH8VFJfj7+1NWfgCVanmDe/v6ha8bzWYz4+MTqASoAQJotVNIpRLBXnpsJmElBBaLxapyONP/2dT0b/CcKSkpiCViVCoVEokEsVjk8C0fHBzM7g/3sG3bNu6/7z7SUlNZv249YB0voFZbh7xGRERw1VVXU1lRxd69e4mLj2Xr1ht4+eWXKCws4Ntvv+FgeQXtbR0UrC7gwQce5KmnnnJ7HTk7nUroujM0NJRRtZqEhASblqsrWCwWunu6KSsvIyw8HJlMSsCCqVzuIiw0jNzcXI4dP0ZXl/ApXv8OWCwWTjWeorW1lcKCIsEPH1iXFYGBQXh6eqJUKCkrL2NoyP1y0+TkJBazRTCJYHhk2JYfEYK+vn7CQgUMAJ2BekyNr5+wUHihrnNTU9OSa84lWx0yMjLtWrP8/fxRq0cdciVFIhE7dtzEjh032f1eOUOtm08QiI6O5qsvv6KlpYWvv/6Kffv3cc8993DNNVtsn8nIyLSNxMvLy+XcTZv49a9+7TLUEolE+Pn5MzIysmhIqiuIRCKSEpPoH+hHq9XaJnY7wqh6lPr6Onx8fCkuKp7p8jBaFdoFFs5n4e2toriomNraEwyPjJCdlWU3UPc/CZ1OR83RI/j7B1BQUPiD2tKamhptM0Bjpqaob6inrb2djPSMJY2uublJ8FxWgP6+PkLDwpb+4AL09fexcuVKwdtNTk4Jio5gpjw0Ez5PTk7S2dnpULd5Ppb0nOnp6ajVarq7rZxQf/8ARkaFUddEItGMOtzkot8nJSVx88238M7f3mXLlmvtHgx/f3/Uo6OUlZXz8suvcPRoDbfcsnT3SFRk5CIOqzsICwtjbExNYGCgwyTI9PQ0R4/V0Nh4iuysFWRlZtnqYwkJiXR2dqJZRn12FlKplJyclQQEBFBWXkZ3d5dgATIhMBqNtLa2UFFZQWJCIqkpqT/IMAcGBpicmiRsxlCUSiX5efkkxCdQc/QI9Q31TntnNRoNU1qt4KSOxWJheGRk0SzRpaDX6zGbTbYw011MTk6iUHgKvk/qeaFw+cFycnPzlpQuXfLVLBaLWbduHfv2fc+WLdcSEBBAR0e7oBMDa2g7MjIiiIkhk8lsSmoFBQV88P4uNp69kVtuudmhyp5IJGLNmrUkJydzsu6k3Xi2K668nMnJSVasyOH+++53OIBWJBKh1xvYtm2b3c23lgVa6T3dS2pKGiEhIYu+HJlMRm5uHkdqDlOwutDt8e6OziEmOobQkFDa2tvYf2A/0TPtXf+unk6dTkdbWyv9AwNER0dTWlL6g/et0WhoOFVPUWHRonsTGBhIacmaGZX6AyTEJxC1IItd31BHelq64Ideo9GgUnkLrqX2D/QLXtsC9PWdJixMeCg8Nqa2EQ6+//571q9fv+Q2bl3R+nVn8P33VllJDw8P9Aa94Dd6QEAAIwuI6dPT07z55hvcdNMOVq7M4frrr1vkXQMDg2yhrbe3N//6v38hk8moqKiw/VRWVlJZWcm+ffs448z1bNi4gd7eXvbt30dfn7UGp/JW4eHhwbhGs6j0Mx+vv/4XVq5caePpnu47TVn5AUQiEaUlawgNDXWpnp6RnsmRI4cFy58shIeHB2mpaTaCR1l5GfUN9QwNDQmWSQFrF0V/fz/HTxyjqroSb29v1q5ZS0J8wg82zNmIIi8332m7lkgkIjYmlpLiUsbHxykvL2NkxPo8DA4OIJPJlzWxu/d077KMZblG1tcvPCtsNpsxGo02frm7xunWt7J+/XpeeHFu1MJsLVFIOtlH5bMo5Hv2uWf56KOP+MmPf8K2bdt5+eWXKC0tYdeuD23xeGREJI1NjYTP3Mjg4GCef/4Fp8fR6/X861//x8effExpaSnXXrvFptUil8vZs/sjPvv8M66//joiIiJZmZPDypUrSUlN5b1//INvvvmGxx77FXV1JzGajCgUSgoLitymagUFBTE1NcmxY0fJzc37QWEiWEPdhIRE4uLi6e/vo7/fqhRhtpjx9fXD388fpZfSpj1oNBoZGRnGYrEwPjHB6Ogo4+MaZFIZfv7+hIdFkJ214gef1yxMJqO14yUj062oSCaTkZGRycTEBHX1J5G0S5mcnGBVvvOx8M6PbWJgoJ/UFGGqFCaTkWnttGA+rXZai1gsEdwvOjExx/mdnJzk+PFjFBcVL7mdW8aZmZlpW3dGRUXhHxDAyMioIOOcT8nz9PTEaDTy8ssvsfvDPeTm5gKwatUqXn75JdauW8Nrr/2F8887Hx8fH7TaKav6t4OUd2NjI1988QVT2immpqbw8/VjzZo1vPrKazQ1N1FZWYWfrx9tbW02r7zp3E001J+ipqaGY8eOUVNTw1//9jfOPPNMao4cpb+/j6HhIXJyVtpeCkIQExPLxOQkjU2Ngh8cZxCLxYSHR9hkTkwmExrNGKOjo/SdnnvpGQwGTp8+DSLrCL2E+ARUKtUPGlfvDBaLhZqjNcRExwoWHfP29qZgdSEn62oZHh6iu6ebxIQEQUmwnt4ewsMjhNMDl9FfClbKXdgyEk9Dw8MEBFgZYmXlZaxcmesWZdCtq5pdd34/M5g2wD9gWb2TgYFBttT6xx//i6ioaJthgjX0ufnmW/hw125uvfUWfv2bX2M2mwkPj+C0Ex7qE088zhdffo5arUYqldLS0sz27dsICQ3mwQcfoKLiIEePHiU1NZXiYuvb6tChQ1x2+aXs27+PK664gpdffoX9+/Zz44030nCqgaCgIPJy8+kTOBVrPtLT0tFoxpaVmHIHEokEf/8AEhISyczMsv0oFArr/2dkERcbh6+v73/EMAHqG+rx9vImOjp6WdubTEaGhoZZv+4MpFIpB8rK6O3tcasMZrFY6Ohotyk6CMGys7vLNM6+vj5CZ7Sdv//+e8444wy3tnP7W5u/7rROphI2ZQogIjyCnl4rNe7FP/2Jn916q8PPlZSUUHGwkq+++orzLzgPi9li224+LBYL3+/7nj/8fidPPvEkDz/0MC+88CJHjx6jqbGZG66/gcHBQd5++y1i42LYvftDnnzqSTZfeAEXbr6I1pYW0jPSePjhh9i3/3uMBgNrSkuJiIgkODgYvV63LEIDWF80uSvz6Ghvd3ju/3/GbKeKVqslNTVt2ftpbW0lOjoaDw8PEuITKCosYnh4mIqKg4yNueYcd3V3ERgYJLgzxGw2M6YR3o6m1+uWld2dbd6YTRC6u94EIca5fr1t1ohIJEKh8HSbE6rT6bBYrMVlo9HI4OAgFRUHufTSy5xuEx4ezt6v93LGGWeyZm0pXV1di9hJ7e3tGAwGh2JggYGBXHjhRfz0xm1cd911fLTnn/zi3l/w9ddfU1VZzS233MIzzz7LRx/9k6joKB544AFMJrMtrBKJRGRnreBkXS16g7Bx6bOQSqUUFBTSd/o0tbUnfnCS6L8BBoOBw0cOM63Vkrsyd9lr16mpKU73nSZuXl3Yw8OD7OwVZGZmUldfz7Hjx9DpFqs9GI1G2traSHbRC+kMA4MDBAUGCScs9NtPNnB/uzlvOzExwYkTxykqLHJrW7eNMyMjg7GxMbq6rJOGAwOCGBxcWuRqfHycVavzufKqK9BqtUSER6IZ11g1YpYIt2QyGffdex/79x3gX//6mJdefsk2UwXg66+/Yv369S5vtIeHB4FBQURFR9HS3Mrer/cSFhZGXX0dx44eJSM9nZt23Mz2bds5a8OZfPnVl7ZtlUolKSmp1NQcWXa9USaTkZeXj9LLi8rKimWNjP9vwdiYdRZoREQEWVnZyw6XzWYzx44fIyszy+E+fHx8KSosIiQ4hMqqClpaW+zuf0tLM7GxscvqV21va3MpDeIMy87unj5ty1vM1jfdbVFz++7Orju/+WYvYPVsvaddh2sWi4WfbruRwoJC5HI55246F4VCQV/faQIDA91uzUpJSeHNN94kOTmJLddew+VXXEZaeioPPPgA11yzZcntExMSaWu1jiHonGnl8vb2pqSk1Ea527r1Rt57731uuOF6euZ1pYSFhhHgH0B9g7DRD/MhEolIiE8gPT2dQ4erGRjoX3qj/yJYLBbaO9o5UXuCvLx8ItzQ3nW1r5MnawkOCnJJdxSJRISHh1tb08xmysoO0N/fz9TUFAODA8QIVKoHq5iaRCKx0152BwaDAZ1OL5hWqNfrMJnn2tj27t3rdkgLAowT4LJLL+O9960KdUqlErPZ7LIheuczO+ns6OSFF17kr2//jdWrV3P2ORsxGIykp6Vbs4puQiKRsG7tOl566WXOP+98dn3wIf19A1xw/gVLbuvh4YHSy4vvvv+OiclJSopLiImOWeRx161dx623/owbtt5g96ZOSkpGNz39g3mv/v4BFBYU0dbexik3+Lv/DTAajdQcPYJGo6G4qETwg70QLa1Wec9EN/s8JRIJSUnJrF5dwOkZ+dCY6Nhlee22trZl0QMHBgcIXUY7Wt+8IV9ms5kPPnify1ws5RZC0BVu3nwhVVVVtpma1iyqYwP79rtveeaZnbz//gd4enoiFov5/dO/58atP+XPf/4TuXl59PULy4YGBQVjsZi56qqrycpyHBIthFY7xZEjhzEY9IjFItJS01yGQ/fdex9Hj9bY8YlFIhE5OTl0dHbaCufLhYeHBwWrrfzVisqDjAqkQv6/gsViYXBwkIMHywkNCWVF9grBwlkL0dPTzcjICFkCRc3AmlCJjY1DqVDS09PNyZPCcgE6nQ7N+LggvvUs+vpOExYuPKQ93ddr266svAwflQ8rVqxYYqs5CDJOpVLJ5gs2896MvmuEkxLHyMgI1133E9588y1iYuzDjzvvvBM/Xz9iYqLpF1iqmCWnN7c0L/lZk8nIqcZTHDp8iOjoaIoKiwkPj6C93fn4CIA9e3aTmJBIVlYWAwMD1NbW0tDQgEQiJT8vn9qTJwRPDHN0HSkpqWSkZ9Da2sLBioMMDPT/V3hSs9lMT083ZeUH6O3tIS8vn0iBcpiOMDQ0RHtHO3m5ecvyehaLhfqGOlasyKG4uARfX18OHjxIe0e7W/etuaWJ+Lh4wS8FnU7H1NSUYNUDvUGPwWCwEeTffecdtlx7raB9CL5L12zZwrvvvgNYvYBEIlmU5Pjt737L5gs2s3HDRof7ePjhX3L06FF0buiuLkRISCgazZhTA7FYrGWXsvIy5HI5pSVrbB00SYlJ9PT2uEzKPP/CC5yoPYG3yosVOdlcfMlFXHjRZsCqnZSVmT1DzzMyMDDA13u/5tlnn+XIkSOCr8XX14/8/FVkZ2XT19dHWfkBenq6/6Nkd2cwGo20trVyoGw/mvFxVuWvJidnpeB1liNoxjXU1Z8kP2/VsqmCp0/3olL5oFKpEIlEREVZOcHT09OUlR+wSek4wsTEOGq1msjISMHH7ejssBsgdcHm80lIjCcuPpb4hDg++OB9h9sN9Pfbaps6nY4Pd3/I1VddLejYgu/UWWeexY09Wzl16hSpqalERETQ29tjk/jr7u7mjTde52iN8/F8UVFRJCenEBAQKFhnViQSkZiQSEtrC1mZWXZ/U4+pqaurQ6VSUVRYvEhRTSKRkJmRyYna47bQciF2f7gbk8lEUFAQUqmUL7/6kqeffhqwepU33nyD7u4uvvr6K1588UWys1eQnJTE757+LV98/iXZ2dmcPHmSV199hQ0brTIuS4WD3t7erFiRw/T0NG3tbbS07ic6OoaY6Oj/eNuYTqejrb2N/v5+oqKiKC4q+beo9s1ienqao0dryHOjC8MZTCYjzS3NFBbYlyCkUmNZVEoAACAASURBVClpqWlMRcdQV19He4e1NW1hNrSuro6M9AzBXtNkMnH6dC9rStfa/r13715OHK/Fw8ODzs5OLr/iMoJDQjhj/Rl2257uO016mlV76fPPPyMrK3tRFLkUBH/zUqmUK6+8inf//i6PPfoYYWHhVFRW2Izzo48+Ys2aNXZj+xzhvnvvo66+jtOne4mIEPZGCwsLp7ml2UYFnJ6e5tSpBrTTWrIys1xOag4ICMRH1U9TcxMpDppdF9K6Gk81kpqSglar5fobrqO/f4B7772XkOBgLrroInJWrEQsFnPmP/7O5gsv4Jyzz+FfH/+LrVtv5PHHH+euu+7knrvvYceOm5Z8ODw9PUlPSycpMckmECaRSPDz88ffzx9/f/9lP+BgjSq0Wi2j6lHUo6OoZwr9sTGxpKxZ+29nEhkMBg4driYrMwvvZejJzqKuro7Y2DinhAOlUsmq/FUMDQ1xpOYwQYFBJCUlI5VKGRjoRyqT2QmsuYveXiupfvbl2t/fT0BAgI33HR0dzTvvvMvVV1/F6d4+2/drNBqZnp7j7r7z7jtsueYawccXuYrXTUazwz8ePnyYa7ZczamGRkQiEdWHqklLTUOlUjE0NERWdibf7P2WjIzFMznmY3p6mkOHqyktWSP4rdbT041arcbD05Pe3h5SklNddozMh8Vioaq6ktjYuCULy3fe+XMCg4L47LPPiI+P47VX/2IzkOaWZoaHhsjNzUUu9+DVV1+hpbWV++69z6ZQWFlZye2330ZySjKvvPyq8ME8Bj3qUTWj6lFGR0fQ6fSoVCr8/fzx8/O196wiOHLkMPl5q6z/FFnv8ejoKGq12tqNr1Di5++Pv58fvr5+/1YvaXfeeh3Vh6pJiE8kfBnJlFn09ffR3dVFfv4qt75bs9lMV1cn7R0dJMTH09bezupVqwWr8lksFg6UHaBgdYHtpVBZWckdd9xORUWl3WeDQ4Kor2uwKdv39PYwMTFBakoqarWahMR4WppbncpgSqRihxe2rJgpLy8PuVxORUUFxcXFxERH097RTnZWNkFBQfzv/Q9wzy/u5pOPP3V5Qz09PfH18WVgYC4+dwcWiwWxWEJHZwcJ8QmUlqwRlEmcpdZVVB7E28vL5Vs9MzOTO++6k5Urc/nr23+zu56kxCR8VCoqKitYsSKH7dt3LNq+sLCQ7777nh07tnPGGevZtetDQeGNXCYnJCTE1oRssVgYHx9nVD1KV3c3ZvMM68gCFqxhanNzk23MsYeHHD8/fyIjo1AoFP+2bhRXUKvVHD9xjPS09GVNF5uFVjtFY+MpigqL3T5vsVhMbGwc4RERHDpUjV6vs2r3CDTOoeEhfHx87Lx1cHAwnV2di+bARkdF09XVZTPOrq5OsjKtE/B27/6QDRs2uNSndXotgrfA+nBfc81cYigkJBS1Wm2red566610dHTw6WefutoNYB3uI2S03fj4OFXVlQwM9JOdlc2UVrusFL9cLmdlTi41R2tcKt/96EfnotVqOfvssx0+ICEhoeTnraK29oTTOqhCoeDtt//KVVdfTemaEg4cOCD4fGcxO3YiNiaWFdkrWJmTa/1ZmUvuTLdDbm4eeTM/mRlZREZEolQq/+OGabFYaG1t4WRdLfl5q36QYer1eg4dPkR21grBauwzJ4PRaGT1qtU0NTVSc/QI2mn3s+ytrS0kLKiJJiQkEBUZtWhUZGRUJN3dVubc2NgM0cEW0r7LFjeIMo6w7EXGlmu28MGuD9BqtYhEIuLj42mbKVPIZDJ+97unue++e5dsDFYqlSiVSpuoszPo9Xpqa09wovY4qSlp5OSsJDo6BrPZxMDggMttncHHx4ekxCSOHq1x+nKYzQq7YnZ4eXlRXFTM0NAQJ04cd8ihFYlE3H3X3bz22l+48qoreP31vyzrnP9bodPpqD5UxfT0NMVFJT8oy2symTh0uJrk5BR0Op1g5X+AU42NJCYm4evrR0FBIRERkRyqrqapqXFJjvPg4CAyqczhpPPLLr+cD3ftsvud1XNau49a2+aMur29nRMnjrNp03mCzx9+gHHGx8dTVFTEm2++CVhrnoODA7bC8HmbziM8PIK//OW1JfeVlJREc7Pj2qXZbKa9vY2DFQfx9/enuKjErmM+MzOLhoZ6TCbh6gBgJVKoVCqnIxhmW74KCwpd7kcikbJyZS7e3ioqqyqcvqV/dM6P+P67fdzzi3v+awkIQjE0NDQjDRpHRkbmD0osmc1mjtQc5sD+A6SmppCRmc4Fmy8QVAPWaDSMj4/b0QxDQ0IpLbUuf8rKD3D6dK/DfZrNZhpO1TucYQpWltyej/bY/S4qOpqu7i60Wi1Tk1M2WuLOnX/gpz/dtuwk3g9Kz9137/38YefvMRqN1lg/JtYmyCwSifjtU7/l17/5tZ0SuCN4e6uQyWSLekQHhwYpKz+ATq+ntKSUSAcj+mZnj8wOXFoOUlPTGNOMOWQ7zR7PnRs8G0GkpqZRXV3ltPaWnJzM2Wefze7djocO//8F1ge5gebmJgoLCpelyTMfFouF2pMn8PP1o+ZoDX987o8M9A9y7NhRm8CcO/uoqz/psHQiFotJSEiksKCIwaFBKisrFqlztHe0Exoa5nSNOjExsUhOJToqip7ubtrb24iLtxId+vv7+fs//s4dt98h4A7Y4wcZZ3FxMTExsbz33j+A/6+98w6Pskr78D0z6cnMpJNGCiQB0nvoYF1dwRWUbl2xfOvakaKuuuqqsGsBsay6q+tKgth1XQuCdEJCOgnpjfRk0maSTKa93x+TjLRAMjOhuNzXxTUhvJk5Q97fnHOe8zy/BwICxtPY1GRayiYkJHDF5Vfwt5f/dtbnCg0No3xw9uzt7SXrcCZ1dXUkJiQxKXzSGQ+vg4KC6erstLj2sqKinJ6TWjlcNvcydNrRlXp5uHuQkpxKZVUFBYWnL3tatnQZ6Vu3mjXeCwFjp+wMJGIxqalTLTriGaK8ohyRSEx7u4Ivv/ySH378AZ1OR19f34iLnBsa6nF0dDqjH5G9vT0x0bFMmRJBUXERBYUFDAwMMDAwQH39MSZOmDjsz+7Zs4fZs2af8L2A8eNpbmmhta3NVIGy6fVNLFm8xKzi7CEsPthas2YNG/66weTzGuAfQN1xgZHnnnueN97YfEqzo5ORy+UIgoHcvFxy83KYEDKRxITEER09iEQi4uLjKSo+MqpN//HY2toSFxdPXl7uKSZj5uDg4EBKciqeHl5kZh6iqLjohL3TNddcS35+3gkVMBcDGo2GktISsnOymTxpMmFh4VYJNNXW1dLT043URcqixTcRFxfPtm3b+PLLLwgLDRvRkY+qV0V1dTWREZEjek25XM7U1Gl4enpyKDODrMNZTJwQesYA4959e5g1a9YJ3/P18SU4OIjxAQGIxWK6u7t59913ePTRVSMax3BYLM6rr7oaO1s7/vOfbwBjb5Vjx37xWw0MDGTatGl8PUxrejAuRY4dO0Zvbx+dHR1Mmzp91AnKjg5Ge47c3Byzi5qlUimxMXFk5xy2ikBFIhF+fn7MnDkLNzc3DmdnUVhYQF9fHw4ODtzwuxtMVT4XOiqVkoLCAg4dysDJ0ZEZ02eYfHEspbmlmcaGBuLjEtj0+iYWLljIY6tW8eknn9Ld3Y2NjQ0azZmT3A0GA/l5eUTHxIwqRVAkEuHna6xP1WgGqKyqOKEx1fEIgsDevXuZPfvE4KCXlxexsbGMHyxje+vtt7jmN9cQEjL6XjHHY7E4RSIRa9as4aX16xEEARsbG3x8fE7wzlm+bDlp6Wmn/fmOjg4OHNyPUtnDrJmz8Pf3p7au1qyxeLh74OvjR1HxkbNfPAxyudwk0LPtlYfDaFj8S37s0A0wY/pMvL29yc3LIS8/j8VLlrA1Pd3ssY41Q5UpmZmHKCoqwmfcOGbOnEVgoHklW6ejvb2d8vJyEpOSkEgkpKen8eZbb3LLrbfg5u7OypV34eDoeMYPd4CSkqP4+vqa1UZCEARKS0pISkwmOSmF+oZ6MrMOnbLaG0oNPdkzqa+/j7y8PPR6Pf39/bz++iZWr14z6nGcjFX+hxcsWEhXVye7Bg3AQoJDqKmpNkXD5s+/noMHD55QXN3f309Obg4VleXExsQRERGJra0toaFh1A9GvswhODgYvd5YVG0ucrmc+LgE8vJyR1VzOsRtt99KWHgoz//l+ROWrSKRiHHjfJg+bQb+fv44Oztz7W+vJT9/+Dzk84Fer6euro59+/fS2NTI5ClTSE2dirf3yDKwRkpjUyMlpSUkJyVjZ2vHwMAAbW1tPP30M/zzH//koQcf5Ouvv0IikZxxqdnc3IxKpTKrVnNoHC5SKTKZDEdHRxLiEwidGEZ+QR5FxUWmc/DDh7NITT0xam90za9i7969tLW18f7775OSkkJUVNTpXmpUWCWrWiKR8Niqx1i//iUum3sZdnZ2eHh40tTchJ+vH87Oztjb26NWq9HrdVRWVdHS3MykSZNOafcmkUiMG/WiIyQljd7LVCQSERMdw8GMg0ilMrMyM8C4xJ06ddpgoXE34aNoVZCXl8fEiaHs2LGDjRtfY/r06ay8cyXXXvtbU19HLy8vZDIZb7zxBqWlJShVPYwb54OHuztyuavFtZOjRRAEenp6aGpqpKW1BV9fP1KSU0dtoDXS16qpqaG1tYXUlFTTftLe3h6VshcHB4fBqOtRtqRtQSwWEz6MxWhHh4LyinJSU83r8aLT6aisrGBq6ok+su7u7kyfNoP6+mMcOLif4OAQ3D08Tjn+Kq8oJzgoCKmLlIaGBl5+5W9s+ej0q8TRYlZu7enQaDSETwrjs08/JzExkf7+frJzDjNj+kzq6+tJnZpCVtZhKisrCAgYT3BQ8BmXRrm5Ofj6+ZllqgTGaOLh7KxRGUKfjqHuW8qeHuLi4kcUmNi/fz/btn3MN//5hoGBAaKjo1GpeqmqqmTJ4iUsX74CV1dXVqxYTmBQIO++8x7Ozs60tbXS0dFBV3cXEokN7m5uuLt74ObmNuIc2L379jBr5uyzXqdWq+nq7qK7q4vu7m761f1IXaT4+Pgybty4Mftw0Gq1FB4pQCKWEB0dY9HyuKenh7z8XJKTU0btijdEaWkJ9vb2Z+wlqtVqqaisoL6+nldeeZnv/vs9YAxA5eXlMmP6TObNv44pUyLIzc1lx087RjWG4XJrrSZOgM2bN/Ofb7/hu/9+j0gkorCwAA8PYySsra2NqVOnMSl80ojEMjAwwKHMDItKmNra2gZLjVIt3iM1NjVSWVlB/GCiwUgQBIFdu3dx880rqKyo4tixY6RvTSctbQuNjY2sf2kD995772k/8TUaDR2dHXR0dNDZ2YEgCMbqFDc3nJ2csXewRywSIxKLjI8iESKRiH379zJ71i8BC4PBYHQB6Omhu7uLrm5jLayDgz1yuSuucjlyuSsODqNvzjNaOjoUHCk6QujEsLNWLZ2NoQ/fhPjEUTu3D9Hb20tObg4zps8Y0f3R1dXFK6++zPz51xMZEcmRokJCJ4bi7u7BXXet5OtvviYtLZ0rLr9iVOM4J+LUarUkJSfyzNPPsGDBQmNPjIP7aW9vp6GhkdWPrR7VoJuaGqlvqCcpMdnsG6eisoKBATWREZbvAXqUPeTl5RIePmlUM/r86+dxw+9u4M47V6LT6aioqMDf3/+06WHDodPp6OoyVqf09fUxMKBGEAQMBgFBMJi+Vip7TFX7Q02g7O3tkUllyOWuyOXyc5YAP4TBYKC8opwOhYLY2LhRV+aczMDAAJlZh4iOijGrv8rQmA5mHCAyImpUzxEWHson2z6lo1OBRCxhxoyZ2NjY8NBDD1JUXMT2H38a9VisWpUyHLa2tmzcuInf//4OQsPC6OzswM/Pj6ysLOLi4s/+BCfh6+tHZ1cXFRVnbzQ6HBMnTCQnN9voDm5m78whZFIZU1OnmQyvwkLDRnSTP/LIo9x//x/p6+tn46bX6O7upq317Laix2NjY4Onp6ep8mE4RrqsPVf09fWRl5+Lp6cXqalTLV7B6HQ6DmcbSxTNFSYY3ep9fXxH/RyhoUY3DanUBV8fX5PjRvrWdA5lZJo9ntNhdZ/+uXPmkpqSyttvv8WM6TOJjIjC3t6eKVOmmPV8kydNRqFQ0GZmcrtIJCIuNp7mlhaLIrhD2NnZkZKcil6nIzPr0Ih8aOfOmYuPjy979+3h2T8/a9FNdTHR0NgwKKQphIeFWyxMg8FAds5hgoNCLKp4aWlpobe316zobujEUNrb2/Hz8ycsLJxpU6fxzNNPc/8f7yfYDD/cMzEmHhgvv/wK8QlxLLppEa9vfp0B9QDTp09nQsiEUS+nxGIx8fEJHMrMwNnZxawlkUQiISkxicPZWYhEItNhsbmIRCKmTImgo0NBds5h/P38CQ4OGfbmE4lE/LTduNzZv38/HqNs+nOxodPpOFJ0BIPBYDXbE0EQyMvPw8vL2ywvoCH6+/spLSs5bR/RkTBp0mRsbW0IGQwg7di5g9q6Wr5a9bXZYxqOMelw4+fnx5rVa7hpkTEN69NPP8Pd3cPs5AJ7e3uio2Isyv4ZEmhDY4Op9s5S3N09mD5tBlqdbtCc68xNeDQaDWvWrmHhwoWm72m1WgoKCvjww3+x6rFVbP9pu1XGdj4QBIGGhnoOHNyPp4cHCfEJVhNmUXERTo6Op9RYjgaDwUBefi6REVGjbuM3RMiEYH7+eRcSiYSBgQEefvghXn31NavkFp/M2LSfAu6//wF8fX2JjorCwcGBsNAw6upqT5sEPhLc3NwICAjgyJFCsy0kJRIbkpOSaWhooLKywipWlBKJhEnhk0hOSqG7u4t9+/ZSP4yD3iOPPIyHhzvhYWHcffddpKam4O7hxrLlS/nhxx+Qy+X8/vd3sH7D+gvCJnOkCIKAQqHg4MEDdHZ1MTV1GgEB5nUeO91zFx8twmDQW9Q0CaC0rBQPD0+zvGsBFB0KpFIp3377HwBefuVlIiIiuPaaay0a13BYNVp7Mj/v+pk77/w9RwqLcHJyoqWlmcbGRuLjE8x6PkEQyC/Ix83NzazWb0MYDAaKiovQ6bTERMda9UxvYGCAyqpK2tvbmBAyAT8/f8RiMR988D4r71qJXC4nIiKSpUuWkpCQQExMzAmFyfX19dy06CYmTAhhy0dpo156ncuA0JA7XW1tLY5OToSHhVlk5HW658/Lz8XF2WVUSSCno6q6iq6uLrObLwmCwP4D+4iLjSdkQjDff/cDv73uWg5lZFq81zwnRymnY8XNy5kwYSLPPfscANk52fj6+Jp9zqXX6zmYcZDIyEjcXM3L/hmitq6W+vpjJCQkmn2IPRwDAwNUVVfR1tZKSPAEFi9ZTEJ8PHfffQ+xsbFn/Fm1Wo2HpzvNTS2jOm6BcyNOtVpNbV2t0Ql9nA9BQcFWX9ZpNBqycw7j5+tncZS9vv4YjU2NJCUmmx2UKi8vwyAITAqfxGWXX4Zer+fqq6/mySeetGhscB7F2dDQQGJSAt/993vi4+PRarVkZBwkLi5+1DfeENbK/oFfDsajo6LNsk88GxrNAFVVVbS2tRIcFIK/v9+IvGgdHO3p6VaO2j9nLMXZ2dlJdU01fX19BAcF4evrNyaZRMbkgGyTo6IltLS0UFlVSWpKitkewAqFgrKyUtNR0O9u+B1ZWZlUVVZb5UPpvIkTYOvWdP787J/JPJSFVCqlR9lDfn4e06ZON9sBvLW1lcqqSlKSUyy+QYz9VHIIDAy0OJI7HBqNhpqaapqam3GVywkICMDd3eO0Syy9Xo+Doz2aAe15X9YaDAaamhqpqa3B0cGR4OAQ3NzcxiSJQRAEjtUfo7amhpiYGORmVJgcj6JDwdGjxaSkpGJna4ZJGL9kqiUnpeDo6Eh1dTXTZ0zj66++ITl59Lnfp+O8ihPgnnvupq+vjw8//DcikYj6hnpaWlpIiE8w+xddW1tDS2srSYlJFp+h6fU68gvysbezZ8qUiDFr1S4IAh2dHdTX19PV1YmTkzMeHh54uHsgk8moqanhrrvvQhCEUedognXEqdVqaWtrpbmlBaWyh3He4wgKCjZZdwiCwNPPPI2joyNBgYEEBgYRHR2NXC43+zU1Gg0FhQXY2dkSMSXS7A/tIXp6usnLzyMlOdXs2U0QBLKyMgkKCmbcuHFoNBrmzJ3DksVLeOihhywa3/Gcd3H29fUxbfpUHnrwIe644/cAFB4pxNnZ2aLweHV1NQpFOwkJiRYLShAEKisraFe0kxCfaJ4l4yhfr7evF4VCgULRTk1NDUVFRfj5+bNgwQLcXEc/Q41WnBqthp7uHrp7uunp7kalUiEWi/H08sJnnA8ymcw0BmOKoIH29nYmT5nEPffcS11dLbW1dRw9WsxvfvMbbr75Fq6+6upRHaG0tbdx9Ggx4WGTLLL1GKK3t5fsnMMkJiRZ5AJYUVmBRqMhYtDs67HVj1FWVsqXX3xl1ZXDeRcnGItVL7/iMpMbvF6v59ChDCZPnnzGRqpno6qqks7BSJw1ZryWlhZKy0qIjooxu+RsNFRXV3PX3XfR19fH22+9jaurK4qOdnp6lDg5OhqT3Z2dcXJyxtnZ+YyzynDiFAQB9YD6BCH29vVha2uLXCZDJpcjl8lxcXE57f9hd3c3V199Fe2Kdq695loyMzPJzMwy/btCoeCTT7ax+Y3N+Pj4mpIuzoRer6ektASVSkVsTKxV9m9qtZrMrEPExsRZNJN3dHRQUnqUqanTEIvFfPvfb/njH+/jcFa22Ucxw3FBiBPggw/e59XXXuXggQycnJzo7+8n63CmRcsPgIqKcpRKJXFmhspPpre3l8LCAqQy2VkNxizhH/94jyeefIJVqx7j4YcePmH/LAgCff19dHV20dvXS19fL319feh1wyRiiIzucDKpDLHYWLGi1WpNZ6b29vbIZXKTEEdqNK1Wq7lu3nVETJnCLbfcyqbXNzJunA8v/+3lU66tqKjgt9ddS1npmd0Qe5Q9FBTkm7KrrPE702q1HMrMYMrkKRZlYWk0GjIOHSQpMRknJydTyePHW7cxc+ZMi8d5MheMOAVB4Pbbb8PBwYG///0dANrb2yivsLy0q6yslL7+fmIsrBM8fqx1x+qoHexe5enpdfYfOgNarTHAI5FIEIlEfPLJNtauW8t/vvnW7Nzjk9mzdzczps/EYDCYbGMs+b/Q6/UsWboYOzs7/v3hR2cNvlVVVXH1b66iorzytP9uLLSupqGxgZiYWGTS4ZtOjQadTkdWViYhEyaYXQMMxgDY4ewsAscH4uPji06n48qrruSaa65h7Zq1VhnryZyTqpSRIBKJ2Lz5DVKnppCensayZcvx9PSis6uLktIS0/reHMLCwqmsrCDrcCbx8QlmR+iOH2tQYBDeXt4cOVJIfX09YeHhpoaoo6Gzs5PAoPEm+xWxWIyLiws7ftppNWEOjflsth6jIS1tC5mZmZSXVYzoOSUSybAplkqlkuLiIqRSKdOmTrfaGNVqNdnZhwkJCbFImIIgUFCYj7u7Bz6DFpfPPvcsDg72oy53tAZjlr53JqRSKWlb0nn4kYcpLTU6rYdODKW3V2WWZ88QIpGI0NAwggKDOZSRcVY7zpHi6OhIUlIyfn7+5OXlkl+QP6JqlONxc3OjIL+Qu++6G1dXV77++hvaWo0dyi5kbrzxJnx9fHnl1VdGdL1KpTpFnAMDAxQUFlB4pNDY0Tsi0mrCHOqdM3ny5FG3kjyeoTRBOzt7QicaW/xt/2k7H3zwPv/64MMxi96fiXM+cw4RFxfHiy+8yPzr57Fn9158fHyIi43nYMYBXFxczE5QAPDx8cHJyYmc3GwmT56CtwXlRUOIRCK8vb3x8vKita2V3LwcpFIpYaFhODqOrFImJCSEe+/9Pz7/4nMipkScc5+g0XLgwAFefuVvyGQy3nrrTRwdHFGqlPj5+nLnnStPuX7nzzu57bZbeWDQ5VyvN5pfNTU3ERYaTnRUtFWjnO2KdoqLi4iPS7DofgFjzEKn0xMTbfS8zc/P59Zbb2Hr1o8tToQwl/Mycw5xxx2/57bbbue6635Ld3c3tra2xMclkJuXQ2+fZb6xMpmM1JSpVFZWUFU18i5mZ0MkEjHOexzTp81gnLcPOTk5FBYWjMgtUKPRcMcdt7P+pfWj7nJ8LhAEgT179xgju2o1V119JVddeRWpU6cCIrZ+vJXOzk7WrltLc3MzYHxPTU1NPPf8c9x22628//4H3PeH+zhacpR9+/cjsbFh5oxZ+Pr6WlWY9Q31lJaUkJKcarEwa2pq6OnpISY6BpFIRHV1Ndf/bj6vb3qdOSd51J5LznlA6GQEQeDBBx+gqLiYb//zLQ4ODnR3d5NfkDfY9NQySwuDwcCRokIEg0BUVLTVZytBEGhpaaGishy53JXQiaHD9tl46umnKCjI54vPvxwzmxBLkhD++re/sm7dWrZ9vI0FCxbi6ORgSiFMS9vCI48+wmOPrTZ6IaWnodFo6O/vx8PDg8TERDas/ytKVQ86rc7YmNjHx+rLwSHLk+7ubhLiEyyOojc01JuscCQSCa2trcyZO5sH7n+Q//u//7PSqM/MBROtPe3r6PWsuHk5BoOB9LStSCQSOjs7KTxSYPERCwxGCGtraGpqIjEhcczsHptbmqmsrMBV7kpwcPAJFRqZmZksWHgD2YdzrHLQPhzmijM9PY0nnnyCp596mhdfepFXXnmV66+fT3NTi8kapby8nIcefpDa2lrWrl3HvOvm4eTkZKyRPXYMuVxOcHAwMpn554tnore3l/yCPFNbeUuF39LaQmVFBSkpqdjY2KBUKrnyyiu45tpr+fMzf7bSqM/OBS1OMAYN5s2fR3hYGJs3v4FIJKKjQ0FRcZHV/FOHMlEsPaA+E8aZtJm6Y3VotVoC/APw8PBk+oxpPPH4EyxZsnRMXncIc8S5Y+cObrnlwDWT6wAAFLdJREFUZrb/+BORkZEsWbqYxsYm7r3nHlasuPmEawVB4Jtvvmb9hvUsXLiQ2Ng4JoSEEDB+vMXR8eEwtuuoG+yebp3EEEWHguLiYlJTUrGzs0Oj0TD/+nmEBIfw1ltvn1MDtAtenGD0Ib388sv43Q038Kcn/wQY7fpLSo+SkpxqlXQ6Va+K3NwcJk4Itdie8Wyo1WoaGuo5nJ1Nc3MTNy68EU9PrzGpmh/CHHHOmDmDe++5h1tuufWM1w0MDNDa2kpDQz0CAhkZGaxfv55Vqx7j0UceHZMb2hjpzcfe3oGIKRFWSQYZ2jYNrcoMBgO33HozarWaj7duG7OEk+G4KMQJRmv92XNmserRVdx99z2AsQKlvKKMlORUq9heGI2NCwHBZEA2VgylLO7ZvRexWES7oh2tRourqysenp54eniYbZlxOswR53PPP0drSwuvv775hO8LgoBSqaSltYXW1pbBiPU4fH18TTmraWlbePyJx6mqrLb6/rKlpYWS0hImT5pstYipqldFTk62Ke9WEAQeeeRhcvPy+O6/3w0bLxhLLhpxgjEF7LLL5/LKy6+waNFiwCjayqoKkhKTrSam5uZmSstKCQsNxdfXz+qf/AaDgTlz57Bs6TL+8Ic/mL6v1+vp6uqiXdGOQqFAr9fh7uaOh6cnHu4eFn0AmSPOY8eOkZKaTH5+AZpBA+qenh40Wg3Ozs6M8x6Ht7f3KR8iWq2W5OQknnrqKRYuvNHsMZ+MTqfj6NFi1ANqYqJjrfb7VqlU5ORmExcbj0xmzEz6ywt/4dNPPuHnn3edN1fEi0qcYOw3Mv/6efzpyT+ZZtD29jaKjxYTFxtntaCDRqvhaHExGq2W6Khoqy45P/30E15+5RX279t/xllFr9fR0dmJor0dRYcCQTA2U3J2Nia6Ozs54+TkNKKZaSTi1Gg0gwLsNrZoVykpLy/H1dWV2Ng4ZDIZMqnsrKJ47bXX+P6H70wO/waDgRdefIHFixYTHm6ez3BnVydHCgsJCgpi/PhAq31gKhQKioqPmIRpMBhYs3YNP/zwPd9/98OYb3HOxEUnTvglifqWW27lySeeRCQSmfaMYaFhphQra9Da1kpJyVGCg0MYHzDeKjfFgoU3sHDBwrPu5U5Gq9WiVPbQ29tHb69qMOm9H4Z+VyJj+p9ELEEskSCRiBGLjY/Nzc2MHx+IRCxGEAQ0Wi1ajQaNVmNKgre1sUUmlyGTypHJZLi4uPDVV1+y+Y032Llj54jG2NzcTGxcDHt272XSpEkIgsCjjz7C5198jreXN/v27R9VjECv15u8l2Jj4iwq9TqZhoZ6qmtqSEpMwsHBAa1Wy8qVd1JVXc1XX35ltT6j5nJRihOMN8G8edcxdepUNm7chEQiQavVkpObjYe7BxMnhlrt01Wr1VJeXoaiQ8HEiaH4+ph/cN7R0UFo2ERqqmtNSyhrMVRXqdfr0Rv0GPQG02NuXg5RkVHo9QZEIrC1tcPOzg47O1tsbe2GfT8ajYbgkCB++P5HoqOjzzqGHTt38Mwzz7B3z14AXnjxBT7Zto2dO3/mjt/fToB/AKtXr2H8+DN/0AmCQGNjA5VVlfj7+RMSMsFqe1dBEKioKKezq8t0JqpSqVi8ZBG2trakp221uDWENbhoxQnGKO7CGxfi7u7Gh//6tynCVny0GK1WY3UHPbVaTUVlBV1dnYRODGPcuNH3pXzvvXf5acdPbE0/t52rLUlCePfdd3jzzTc5cODgWQMjzc3NxMRG09Lcil6vx8/fl1dfeZUVK26mra2NP95/H3v37qW1tZUjhUVMnnyqrWV7ezulpSW4uroSFhZm1cCYwWCgsLAAsVhMZGQUYrGY9vZ25l8/n8iICN5+++/nPCo7HMOJ87ym740UmUzGt//5FpFIxHXzrqO7uxuxWExUZBTubu4cysxArVZb7fUcHByIiowiMSGJ1rZWDhzYT2try6hSANO3bmXZ0mVWG9O5YOXKu4iMiuShhx4867Xjxo1DLBbT1NSEjY0Nn3/2BaseW0VWVhZeXl58vHUbx+rqsbW1paGx4YSfHUpWr62rIS4unshI802eT4dWqyUrKxMXqZSoqGjEYjG1tbXMmTubyy+/nHfffe+CEeaZuCjECcZC4bQt6URGRHD55ZeZqleCgoIJD5tEZtYhuru7rfqajo6OxETHEB+fQFNTEwczDtDW3nZWkdbX11NYWMA1Y2Q2PFaIRCLeevNt9u3fx5YtH5312qio6MEjKZg5cyZ/+MN9fPTRv03XfP75ZyQmJnH5ZZcbjac7FOTkZHPkSCFhoeEW24icDqVSScahDAIDA5k4YeJgK8pC5sydzT1338tfnv/LOU0wsISLRpxgrBXcuHETC2+8kdlzZpnKzTw9PUlMSKKgMJ+mZvNLzobDycmJ2Ng4YmJiqT92jIxDB1EoFMNev+2TbdzwuxvG9Px0rJBKpaSnbeXRVY9y9OjRM14bHRXFkSNHTH/vUCgIDDKafRsMBl544QXWrVtHbV0t+/bv5VhdHSEhE5g6dZrV7V8EQaCquoq8/DxiY2Lw9TVGX3fv2c1vrrma9S+t54EHHrDqa441F5U4wfiJ/cTjT/D4useZM3c26enGFt/Ozs5MTZ1GXV0d5RXlY9LOwMXZhfj4BKIio6mprSbj0EFaWppPab2wNT2dpcsuriXt8cTExPDCX15g6bIl9PYOXx3k6uZG/bFf+s44OTvz9NNPMSViMkuXLWX+9fORurig0+lITUklLi5+TGw1+/v7TFub6dOmI5PJMRgMvLT+JZYvX8aHH/57zNMmx4KLIiA0HHl5eSxfsYyZM2by2msbcXJywmAwUFJyFKVKSUx0jMVVLWdCqVRSd6yW9vZ2vLy8GR8wHhcXF+zsbenvU5+XfY21fGsFQeCPf7yPH7f/yHPPPsfixUtOiKL29/cTPimMr7/6xlQwbjAYaGhooLSslF6VCl9fP5KTzW98PJIxNjQ0UFVdSWRElMl4q6WlhdvvuI2+vn62fLSFgICAMXl9a3FRR2vPhFKp5L77/kB+Qb5xTxppLJZVKNopKi4iKCiYQCseZp8Og8FAS0sLx+rr0Gl1bPjrBt7/5wfnxLnvZKxtKr1r9y7WrlkDIhE7ftph2iNmZ2dz+x23UVhw5IT2DN5e3gQGBZll5TIaNBoNhUcKsJHYEBERacqq2vnzTm6//TZuu+12nn7q6Ysi8HPBeAhZG6lUyr/+9SH/+tcHXHHl5bz4wovcfvsdxkqQaTMoLS0xtSgfqzMtsViMr68vvr6+qNVq3N3dyc3LwcnRyZg/6+mJq9z1oglEHM/cOXP59tv/MjF0gkkAgiDg7+9PdHQ0Bw7sBxEEjg9k5oxZ58TdobW1laMlR5kUHm5KRNHpdDz3/HP885//4J//fJ+rrrxqzMcx1lz04gTjPvT22+8gJSWVZcuXsnPnTt588y2kUimRkVEoOhRkZx8mcNCdfCxF4uDgwMGDB1mzei1BQUG0t7dTW1tLQXcBTk6OeHp44enpiYuLy0Uj1j17djNr1iw6uzppbW2ls7MDZydnRCIRMpmcqKioczKO43Nup6b+0ienoaGBm2+5GVtbW7IyD49pvey55FchziEiIiI4eCCDhx9+iJTUZNK2pBMfH4+HuwfTp0+npLSUQ5kZxETHWnUWramp4eNtH1NRXk5SUhIGg4Genh4cHBwICAggICDgF3f39nbKykpR9Rr9ZT0HZ9ax3BuPFo1GQ4+yB2VPDz1KJf3qfhYtWkR7ezu+Pr5ERkQiFotRq9WUlBw9J+Ls7OzkyBFjzm3U+F+8iP773X+5666V3HffH1mzes0F78s0Gn5V4gTjscff//4OW7emc+1vr2HlyrtYt3Ydzs7OREZE0tHRQXbOYbw8vZgwYaJFNaLGouMNVFSUs3DBQuLjEzh06BBOjk6n+NqIRCJcnF1wcXYhKCgYQRDo6emhvb2dgsICBgYGcB50dHd2ccHF2RlnZ5cxbQlhMBhQqVRGISqV9PT0oFarsbOzRSaVIZXKCAoMwsXZhaVLlxAcEsy777xnCgxFR8dQUFjITTctGrMxqtVqysrL6O3tJSEx0bSXbW1t5fHH17Fj5w7S07cye9a56Ul6LrnoA0JnorGxkdVrHjO6yP3tZW64YYGpeqKhoYHqmip8fHwJCQ4xq0zrvvv+gEaj4c0337K4ztRgMNDf34+qV0WvSkVvby+q3l60Wg0SiQR7O3uTH+2Z/hwtOUpsTCxisQSdXod2MPFdq9Uak+C1WrRajcmQzMXFxSjE4ypRTrfcVqvVrF79GN//8D3paVtJTEzks88+ZUvaFj7/7AuL3vvp0Ol0VFZW0NraSmhYGD7jfBCJROj1ev7+97/z3PPPsmLFzTz1p6esnrt8rvnVRmtHwq7du3jggfvx9/dn42ubTOVMBoNh0NG9lgB/f4KDg0fVw3HnzztZt3Ythw5ljtXQAeONqtFoTInuep3e+PXQH8MvX9fUVOPn54dBb0BiY4OdrS22trbY2tn98rWtHQ4ODmYlmH/66Sfc/8D9LFu2nO+//464uDjStqRb7b0aDAZq62qpq6sjeLBsbGicBw4c4IEH7kcml7Np46Zzttcda/6nxQnGfMvNmzfz0voXWbnyLh5f97jpWECv11NbW8Ox+npj0Gh84Ij2LjqdjqDgQHbv2kNoaOhYv4URcS46W1dUVPD665tYsHAhc2bPsUpgSxAEmpqaqKgsx8fHlwkhE0zHIK2traxbt5btP21nw/oNLFmy9KIJpo2Eizrx3RrY2try8MMPk5uTR11tLdExUXz++WcIgoBEImHChInMmD4DnVbL/gP7qKurOyXz52RsbGxYuGAhn3726Tl6FxcGoaGhbNy4iblz5losEoPBQHNzMwcO7qejQ0FqylTCw8KxsbFBp9PxxhtvEBMbjbuHB0cKi1i6dNmvSphn4n9m5jyZXbt38eCDD+Dr68vzz/2FpKQk079ptBpjq/jWFiZOmIifn/+wN8TPu35m3dq1ZGQcOldDPyPnYua0Bmq1mrpjdTQ1NeHp6UFwUIhpJSMIAj9u/5HH163D1c2NTRs3mZJLfo38apMQzGXunLkczsrmnXfe4aZFNzJp0iTWrFnLZXMvw87WjsmTJhMSHGJ0jK+uYnzAePz8/U+xf+xQKPDwtG6/xl8rgiCgUCiora1BPaAeTFyYYdrn6/V6Pv/8MzZs2IBGq+HJJ57kppsW/c/MlCfzPztzHo9GoyEtbQsb/roBudyVdWvXMm/efFMgYmBggIaGehoaG3B2diFwfCAeHh6IRCLuvvsuoqNjuP/++8/zuzByIc6cGq2G+vp6GurrkcnlBAUF4Sr/xUxLo9Hw0Uf/5q9/+yvu7h6sXbOG666bd16aB50P/ucDQiNBr9fz5ZdfsH79etQDalY/tpolS5aekLbW1dVF3bE6uru78PHx5aabbuTTTz4z29DK2lxI4uzq6qK21tiHxJiMMf6EIyeVSsV7773HaxtfZcqUKaxZs9ZqAaaLiUviHAWCILD9p+2sX7+emppqHn1kFXfccccJ1h1DzVoPZR5i+rTp+AcE4OXpadWKfnM43+JUq9W0tLZQf+wYDo6OBAUF4eHucYLgOjo62PzGZt56601mz57N6sfWkJiYeN7GfL65tOccBSKRiKuvupqrr7ra6Gq+4SWefe7P3HTjTSxfvoJp06ZhY2PD/gMHqK2pYeWdK2lsaqK2thaDQY+Huweenl64u7v/qtLJTocgCHR2dtLa1kp7exs2Ehu8vL1JHHS6G0Kn07H9p+2kbdnCd99/xw2/u4FdP+9m0qRJ53H0FzaXZs4RUl1dTfrWdNLStjAwoMHR0YGuri7eeONN5s+bb7pOp9OhUChoV7TR0dGBrY3tYP6sF3K5fMyXbOdi5hwYGKCtrZXWtlaUSiWurm54e3vj6eF5wrJVEASysrJIS9vCtk+2ERwcworlK1i0aBHe3pb3TP21cGlZayUEQSA7O5tn/vw0GRkZTJgwkRXLV7BkyZLTVkOo1Wra29tpV7TR3d2Di7OzUaxeXmNS82hNcQqCwMDAAEqlEuVg/q1SpUQsFuPt5Y23tzdSqeyUD5yKigrS0tNMLhXLli1n+bLlF0yixoXGJXGOAXq9np0/7yQ9LY2vv/ma5ORklixewhVXXMn48eNPuV4QBFQqlUmsfX392NvbGRPiXVwGE95dcHBwMHuGNVecWq0WlUpJj1KJsscoRK1Oi729PVIXKVKZDKlUitRFespSXRAEysrK2L59O+lb06murmLJ4iUsW7Z8TJ0Qfi1cEucY09fXxzfffM3nX3zO7t27kcvlzJk9hzlz5jB79pxhxarRaFCpVKaEd1VvL2p1PyKReLAVgxO2traDie022NgYHyUSyUlfGx/37d9rEqcgCOj1erRaDVqt7rhHLVqddnBW7KG/X42tjY1RfFIpUqlRiMNVxAyJcffuXezes5vdu3dja2vLZXMvY9HixVx5xZVWaTj1v8IlcZ5DDAYDxcXF7N69m917drFnz54RifXk5+jt7aWvrxetTmdMbNfp0A0muOuGvqfXodMZH/V6PUqlEqmLFAEBESIkNpLBZPfBPza/fG1nZ49UKj3rTD2cGOfOmcucOXOYM2cuwcHBl2ZIM7kkzvPI6cQqk8mIiYkhLCyc8LAw42N4OF5eXhbd5JbsOQ0GA42NjZSVl1FeVkZZeTnl5WXk5uZiY2NzSYxjxCVxXkAYDAaOHj3K0aPFRgEMCqGsrBSDwfCLYMONj4GBQcgG93wuLi5IpdJhl43DiXMosKNSqVAqlXR3d1NdXWUSYHlZOeUV5chkslNePyoqmpCQkEtiHCMuifMiQaFQUFZWZhRMRQXl5WXU1NSiUv0iLJVKhY2NjUmoUhcpzi4uSKUudHR0IJHY0NurGoyyGq8XBMG0p3RxcUEmkzMhJMQkwLCwcEJDQy/6wuWLEbPEeYkLE5FxCnMApIDLSY8GQAmoBh+VgEoQhIHzM9pLmMslcV7iEhco/xtp/5e4xEXIJXFe4hIXKJfEeYlLXKBcEuclLnGBckmcl7jEBcr/AyyDGBPZrkeiAAAAAElFTkSuQmCC\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.3"
Maximilian Schanner's avatar
Maximilian Schanner committed
231
232
233
234
235
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}