{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Distances from parallaxes\n", "\n", "In principle, if we have a very good measurement of the parallax, we can simply inverse it to get the distance (in kpc, i.e in 1000 pc, where 1 pc=1 parsec = 3.26 light-years). This would work typically if the relative error on the parallax is below 20% at most. If the relative error on the parallax is larger than, then the distribution of distances is no more Gaussian and we cannot make such an approximation. You can see this below, by playing with the relative error on the parallax (relative_error_plx) and see how this changes the distribution of distance." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "from scipy import stats\n", "import pylab as plt\n", "%matplotlib inline\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we define a function to compute and plot the distribution of the distance, assuming a normal (i.e. Gaussian) distribution of parallaxes, based on the relative error. The parallax is fixed to 1 mas (i.e. a nominal distance of 1,000 pc), as everything can be scaled to it (this is particularly true, because the error does not depend on the parallax, but on the brightness of the object, its colours, the number of observations and the position in the sky). \n", "The mean distance should thus be this nominal distance, and the standard deviation should be given by the error_plx (as the parallax is 1). Note, however, how quickly the distribution of distances becoms non-Gaussian and the standard deviation should be asymmetric, and how the mode becomes different from the mean or median." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parallax= 1.0, with a relative error of 1.0%\n", "Distances vary between 960.7 and 1045.9 pc\n", "Formal distance: 1000.0 +/- 10.0 pc\n", "Mean distance: 1000.0 +/- 10.0 pc\n", "Median of distance distribution: 1000.0 pc\n", "Mode of distance distribution: 999.0 pc\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd4FFX3wPHv2d3QpFfpoYnSERBERGwIqKC8KioiCIi9vvrqT5SiYO9dQKQIWFARBQu9t4ABaVJCgBAgJJCQhIQku/f3x0xigCQkkGSym/N5njzZvTM7e2Z2d87MnTv3ijEGpZRSxY/L6QCUUko5QxOAUkoVU5oAlFKqmNIEoJRSxZQmAKWUKqY0ASilVDEVkAlARD4XkZfyaVn1RCRBRNz288UiMjQ/lm0v7zcRGZhfy8vD+44RkWgROZTL+Y2INLYf59v29Uf296FhDtPDReS6fHqvfP2++QsRCba/cx6nYwlkfpcA7B9XkojEi0isiKwUkQdFJGNdjDEPGmNeyeWycvyhGmP2GWPKGmO8+RD7KBH5+rTl9zTGTD7fZecxjrrAf4FmxpgL8/r6/Ny+/sj+PoQBiMgkERnjdEy5JSIlRGSm/dkYEel2lvkri8hPIpIoIntF5O5CCtVv5ec2FpHWIrLFPlh7KlN5kIissX/L58zvEoDtZmNMOaA+8DrwHPBlfr9JAB991AdijDFRTgeiHLEcuAfIzdnfJ0AKUAPoD3wmIs0LMLZ8l9XvOK+/7XPYF+TXNn4NeAZoDbwoIukHbE8DPxhj9ucxrlMZY/zqDwgHrjut7DLAB7Swn08CxtiPqwK/ArHAUWAZVuKbar8mCUgA/gcEAwYYAuwDlmYq89jLW2x/KGuBOOBnoLI9rRsQkVW8QA/7Q061329jpuUNtR+7gBeBvUAUMAWoYE9Lj2OgHVs0MDyH7VTBfv0Re3kv2su/zl5nnx3HpGxe/yxwEIgEBtvv3fh8tq89//dYP4o4e/s2z/Sek7B+DHOAeGAN0CjT9ObAPPt9DgMvZNpuzwO7gRjgu0yfSSnga7s8FlgH1Mhife8Dfsn0fBfwXabn+4E29mMDNAaG2Z9nir2Ov2T6zJ8BNtnr+S1QKpvtPAhYAXxkz7sduDbT9Mzfj8+AmZmmvQEsAOQcf0sRQLccpl9gr9tFmcqmAq9nM/9lwCp7Ox8EPgZKZJpugAeBncAx+7MWe5obeBvrex0GPEKm310W71UL+AHr+70HeDzTtFHATPtzPw4MzaasJPA+1nc80n5cMvNvGevg8hAw1YltDGzLFNNqexvXw9r/BJ3PvtQYExgJwC7fBzxkP57Evzuo14DPgSD778pMX7pTlsW/O9kp9gdTmqwTwAGghT3PD8DXmb802cVrfwm/Pm36Yv79gQ/G2vE0BMoCP6Z/8TLFMd6OqzVwErgkm+00BSs5lbNfuwMYkl2cp722B9YONn0dp5N9Asj19s20juX498cXmmnaJKyd+2WAB5gGfGNPK4e1U/kv1k69HNDRnvYk1o+jjr3cL4AZ9rQHgF+AMlg7mXZA+SzWuSHWjssF1MRKmgcyTTsGuOznWW6L0z7ztVg7qcpYP+IHs9nWg4A04Cl7+/XDSgTpCSzz96OM/TkOsrdzNFDHnlbPjj+7v7uzeO+z7ZzaAkmnlT1DpkR52rR2QCf7swu21/vJTNMN1sFCRTveI0APe9qDWMmvrr3NFpFNArA/o/XACKCE/fmEATdk+p2lArfY85bOpuxl+3tTHagGrAReyfQbScNKsiXt+Qt9G2MdMN2M9d0+BFQBZuW0zLz8+WsVUFYisb44p0vF+kHXN8akGmOWGXvL5mCUMSbRGJOUzfSpxpjNxphE4CXgjvSLxOepP/CuMSbMGJMA/B9w52mnn6ONMUnGmI3ARqxEcAo7ln7A/xlj4o0x4cA7wIBcxnEH8FWmdRyVw7x52r7GmIl2TCft5bYWkQqZZvnRGLPWGJOGlQDa2OU3AYeMMe8YY5LtZayxpz2AdTYUkWm5t9nbLRXrR9PYGOM1xqw3xhzPIq4wrLOONsBVwB/AARG52H6+zBjjy2E7nO5DY0ykMeYoVgJqk8O8UcD79vb7FvgHuDGLGE9gVSu8i3Uk+5gxJsKets8YUzGHv+l5iD1dWaxklFkcVvI9g71tVxtj0uzv3BdY2y6z140xscaYfVg7+fTtcoe9Dfbb2+y1HOLqAFQzxrxsjEmxP7vxwJ2Z5llljJlljPFl+h2fXtYfeNkYE2WMOQKM5tTfiA8YaYw5af/mnNjGzwAPAbOxDhKuwPqehonIzyKyRERuP4f3BaxMHShqYx09nu4trB3CnyICMM4Y8/pZlnW2erXM0/diHblVzV2YOaplLy/zsj1YdYPpMtcpnsD6Ap2uKtaR0enLqp2HONaf9trs5Hr72olpLHA71hFX+g61Kv/+CLJbv7pYVTxZqQ/8JCKZd9BerO021X7tNyJSEWvHOdwYk5rFcpZgHfk1th/HYu3ALref58Xp61Erh3kPnJY092Y3vzFmrYiEYR21fpfHmPIqASh/Wll5rB3QGUTkIqzk1B7rbMXDqd8jyP7zrcWZv6vs1AdqiUhspjI3VvVjuqx+w6eXZfV7y7zdjxhjknOIIz/kuI2NMXuBXgAiUgbrLOUGrCrDb7GqSzeLyAI7ceZJQJwBiEgHrJ3b8tOn2UeK/zXGNMQ6lXpaRK5Nn5zNIs92hpD5yns9rKPMaCAR64ufHpcba0eX2+VGYn25My87Das6Ji+i7ZhOX9aBXL7+IGeuY5byuH3vBvpgXYeogFVNACC5iGk/0CiHaT1POxorZYw5YB9VjzbGNAM6Y51J3JvNctITwJX24yVYCeAqsk8A+dGdbm2xs6etHtZ34Qwi8ghWlUQk1nWr9PL05srZ/fU/h7h2AB4RaZKprDWwJZv5P8OqxmlijCkPvEDuPlvIw3cO6/Pec9rnXc4Y0yvTPFl9LqeXZfV7y7zdT5m/CGzjEcAEY8xhoCUQYoyJw6pmanwO7+3fCUBEyovITcA3WHXrf2cxz00i0tj+gR3HOjJMb9J5GKv+MK/uEZFmdkZ+GevCnBfrwywlIjeKSBDWhdeSmV53GAjO3GT1NDOAp0SkgYiUBV4FvrWrQ3LNjuU7YKyIlBOR+litBr7O+ZUZvgMGZVrHkdnNmMftWw7rukUMVqJ8NQ+r9StwoYg8KSIl7fXqaE/7HGtd69sxVRORPvbjq0WkpZ2Mj2Mlxuya9C4BrgZK21Ury7Cuh1QB/srmNef6HcqsOvC43bTvduASYO7pM9lH2GOwqoEGAP8TkTZwSnPl7P6mZVpOSREpZT8tISKlTktA2MtMxLoO9bKIXCAiV2Al8KnZrEc5rG2cYFedPZSHbfCdvQ3qiEglrIv62VkLHBeR50SktIi4RaSFfSCYFzOwWtZUE5GqWDvYbH8jTm5jEWmGdXDymV20B7hGRGoATbCugeaZvyaAX0QkHutIYDjWaed92czbBJiPdaq1CvjUGLPYnvYa1hcgVkSeycP7T8W6+HcI64Lk4wB2Nn4YmIB1tJ2IlZ3TfW//jxGRDVksd6K97KVYH3Ay8Fge4srsMfv9w7DOjKbbyz8rY8xvWBdoF2JdlF6Yw+x52b5TsC+uAluxLsDlijEmHrge6yzjEFZLkqvtyR9g1ZH+aX8vVgPpyeFCrNYfx7EuSi4hmx+5MWaHvR7L7OfHsbbfCpP9fSBfAs3sdZyV2/U5zRqs7RiNVUV2mzEmJvMM9vWMr4E3jDEbjTE7sY6wp4pIydMXeBb/YLXOqo11rSMJ+0hYRF4Qkd8yzfsw1gXQKKwd5kPGmOzOAJ7BOsuLx6qT/zYPMY23Y9kIbMDaKWbJ/ixuxrp+sAdru03AOqvMizFACFZrrb/t982vezryext/AjyR6Xv4f1j7nS3Aq8aYXN3Qebr01hpKKQeIyCCsVj5dnI5FFT/+egaglFLqPGkCUEqpYkqrgJRSqpjSMwCllCqmivSNYFWrVjXBwcFOh6GUUn5l/fr10caYamebr0gngODgYEJCQpwOQyml/IqI5HQndQatAlJKqWJKE4BSShVTmgCUUqqYKtLXAJRS+SM1NZWIiAiSkwu6c0tVmEqVKkWdOnUICgo6p9drAlCqGIiIiKBcuXIEBweTRZ9kyg8ZY4iJiSEiIoIGDRqc0zK0CkipYiA5OZkqVarozj+AiAhVqlQ5r7M6TQBKFRO68w885/uZagJQSqliSq8BKGU7mpjCqt0xrAs/Ssjeo+yNPpExLcjjomXtCnQIrkT74Mq0r18Jj1uPn/LC7XbTsmXLjOezZs2isO70nzRpEiEhIXz88cdnTPv9998ZMWIEx48fp1SpUjRt2pS33nqLevVyGpTs/I0YMYKuXbty3XXXFej75OSsCUBE6mIN5HEh1hiu44wxH4jIKOB+4Ig96wvGmLn2a/4PGII18tLjxpg/7PIeWIN3uLGGNjvb2LxKFbi9MYmMWxrG9+sjSEnzUdqTTFvfLv4jEbgirEGtEm9OI3R/LG//uQOAupVLc/+VDbm9XV1Kl3A7Gb7fKF26NKGhoXl+XVpaGh5PwRyrbt68mccee4zZs2dzySWXADB79mzCw8MLPAG8/PLLBbr83MjNVk0D/muM2SAi5YD1IjLPnvaeMebtzDPbQ5fdCTTHGmB5vj2UHVij2lyPNUrWOhGZbYzZmh8rolReHYk/yZg5W/llYyQel4u+l9amX4fXaJH6I0FT7dH89n5k/b/N+hc7shYrfC2YcKwXI35O4oP5O3nsmsbce3kwLpfWsedVcnIyDz30ECEhIXg8Ht59912uvvpqJk2axJw5c0hOTiYxMZERI0YwcuRIatSoQWhoKH379qVly5Z88MEHJCUlMWvWLBo1asQvv/zCmDFjSElJoUqVKkybNo0aNWpk+/5vvPEGL7zwQsbOH6B3794Zj8ePH8+4ceNISUmhcePGTJ06lTJlyjBo0CBuuukmbrvN+mKULVuWhIQEDh48SL9+/Th+/DhpaWl89tlndO7cmSFDhhASEoKIMHjwYJ566qlTlvHyyy/zyy+/kJSUROfOnfniiy8QEbp160bHjh1ZtGgRsbGxfPnll1x55ZX5tv3PmgCMMQexBmzGGBMvItuwhjnLTh/gG2PMSWCPiOwCLrOn7TLGhAGIyDf2vJoAVKH7ffMhXvjpbxJOpnF/14YMuaIB1cuXAg5CeHajP0JFSeRG9xp6udawdmAYHy3cxahftvLHlsO8fUdralcsXXgrcY5G/7KFrZHH83WZzWqVZ+TNzXOcJykpiTZt2gDQoEEDfvrpJz755BMA/v77b7Zv30737t3ZscM6y1q1ahWbNm2icuXKLF68mI0bN7Jt2zYqV65Mw4YNGTp0KGvXruWDDz7go48+4v3336dLly6sXr0aEWHChAm8+eabvPPOO9nGtGXLFp55JvvRYPv27cv9998PwIsvvsiXX37JY49lP0rr9OnTueGGGxg+fDher5cTJ04QGhrKgQMH2Lx5MwCxsbFnvO7RRx9lxIgRAAwYMIBff/2Vm2++GbDOgNauXcvcuXMZPXo08+fPz/b98ypP51UiEgy0xRrD9ArgURG5F2tczf8aY45hJYfMY71G8G/C2H9aeUdOIyLDgGFAgZ+CqWJmVAVOmJKMSBvETO9VtKhdnvfuaEOTGuXyvCgR6NiwCpc1qMx3Ift5+Zet9HhvKa/c0oJb2uZ0fFR8ZVUFtHz58owd6sUXX0z9+vUzEsD1119P5cqVM+bt0KEDNWvWBKBRo0Z0794dgJYtW7Jo0SLAut+hX79+HDx4kJSUlDy1j4+JieHaa6/lxIkTDBs2jGeeeYbNmzfz4osvEhsbS0JCAjfccEOOy+jQoQODBw8mNTWVW265hTZt2tCwYUPCwsJ47LHHuPHGGzPizmzRokW8+eabnDhxgqNHj9K8efOMBNC3b18A2rVrR3h4eK7XJzdynQBEpCzwA/CkMea4iHwGvAIY+/87wGAgq/NgQ9Ytjs4YjcYYMw4YB9C+fXsdrUblm2OmLPelPMsm04jHrmnMY9c0oYTn/C7kigj9OtTj8oZV+e/3oTz5bSj7j57g0WsaF9lml2c7Ui9MOQ1IdcEFF5zyvGTJkhmPXS5XxnOXy0VaWhoAjz32GE8//TS9e/dm8eLFjBo1Ksf3b968ORs2bKB169ZUqVKF0NBQ3n77bRISEgAYNGgQs2bNonXr1kyaNInFixcD4PF48Pl8GeuQkpICQNeuXVm6dClz5sxhwIABPPvss9x7771s3LiRP/74g08++YTvvvuOiRMnZsSQnJzMww8/TEhICHXr1mXUqFGntO1PX0+3252xnvklV99+EQnC2vlPM8b8CGCMOWyM8RpjfMB4/q3miQDqZnp5HSAyh3KlCtyB2CRuSxnJVlOfz4Le57/dm/678x9V4d+/8OXntPx6H9Vk+sGb6OtaxjvzdjBq9hZ8Pj1+OZuuXbsybdo0AHbs2MG+ffto2rTpOS8vLi6O2rWtM7DJkyefdf7//e9/jB07lm3btmWUnTjxb+uv+Ph4atasSWpqakacYHVVv379egB+/vlnUlNTAdi7dy/Vq1fn/vvvZ8iQIWzYsIHo6Gh8Ph//+c9/eOWVV9iwYcMpMaTv7KtWrUpCQgIzZ848x7XPu9y0AhLgS2CbMebdTOU17esDALcCm+3Hs4HpIvIu1kXgJsBarDODJiLSADiAdaH47vxaEaWysysqngFfriXBVGJKidfp5Npu7exzI3yZ9X/UTWedNUi8vB30OZXTjjNh1Y3EJKbwXr82BGlz0Ww9/PDDPPjgg7Rs2RKPx8OkSZNOOdLPq1GjRnH77bdTu3ZtOnXqxJ49e3KcP/1C8r333kt8fDxVqlShXr16jB49GoBXXnmFjh07Ur9+fVq2bEl8fDwA999/P3369OGyyy7j2muvzThbWbx4MW+99RZBQUGULVuWKVOmcODAAe67776MM4bXXnvtlBgqVqzI/fffT8uWLQkODqZDhw7nvP55ddYxgUWkC7AM+BurGSjAC8BdQBusapxw4IH0hCAiw7Gqg9Kwqox+s8t7Ae9jNQOdaIwZm9N7t2/f3uiAMOp8RMYm8Z/PVpLqNUxJeYpmrn05v2BQovV/kl39MOlXuzyLBDAqLtPjfxOKMfD51Rt44/ft3NauDm/d1srx6qBt27ad0tJFBY6sPlsRWW+MaX+21+amFdBysq7Xn5vDa8YCZ+zc7fsEsn2dUvkpLimVQV+tJT45jW8f6ESzcWfZ+edVNmcRIvBQt0YkpXr5cMFOLixfimduOPdqDaUKit4JrAJScqqXYVNC2BOdyOT7LqN5rVxW+eSjp65rQtTxZD5etIsaFUoxoFP9Qo9BqZxo5aQKOMYY/jdzE2v2HOWdO9rQuXFVR+IQEcbc0oJrXRsYOWsTC1+6ypE4lMqOJgAVcKas2svsjZE8e0NTereu5WgsHreLj4M+5BLZy1OpD7P/6Imzv0ipQqIJQAWU0P2xjJmzlesuqc5DVzVyLpBMTUtLSwqfBn2ADxePTt9ASprv7K9XqhDoNQAVMGJPpPDItA1UL1eKt29vjevlik6HlKG+K4q3gr7gwYineHXuNkb1Ljo3Y6niSxOA8n+jKmAM/Df1GaKkHTMf7EzFMiWcjuoMPdzrGNKxAV8u30OH4Mrc2Kqmc8Hk9j6IXC8vLsfJ6V1Bp6am4vF4GDhwIE8++SQul4uQkBCmTJnChx9+mOVrw8PDWblyJXffnfVtQ5GRkTz++OPMnDkzx26fszNp0iS6d+9OrVpWdeHQoUN5+umnadasWa6X4a+0CkgFhG+8V7PAdykv9LqE1nWLzpH/6Z7veTFt6lZk+Ky/iYovPgO0p/cDtGXLFubNm5fRsRlA+/bts935g5UApk+fnuW0tLQ0atWqdV53z06aNInIyH87JZgwYUKx2PmDJgAVAA6YKoxN68/lri0MvDzY6XByFOR28c4drUlK8TL8p8059oUTqKpXr864ceP4+OOPMcawePFibrrJutFuyZIltGnThjZt2tC2bVvi4+N5/vnnWbZsGW3atOG9995j0qRJ3H777dx88810796d8PBwWrRokbH8/fv306NHD5o2bZqRZE6f5+2332bUqFHMnDmTkJAQ+vfvT5s2bUhKSqJbt26k34A6Y8YMWrZsSYsWLXjuuecyXl+2bFmGDx9O69at6dSpE4cPHy6MTZfvNAEov2aM4fnU+/Hh4k3POKveP/0CbBHVqFpZnunelHlbD/NzaPHsDqthw4b4fD6ioqJOKX/77bf55JNPCA0NZdmyZZQuXZrXX3+dK6+8ktDQUJ566inA6ip68uTJLFy48Ixlr127lmnTphEaGsr3339PTr0J3HbbbbRv3z5j/tKl/+3OOzIykueee46FCxcSGhrKunXrmDVrFgCJiYl06tSJjRs30rVrV8aPH58fm6XQaQJQfu2bdftZ5mvF/3mmU9d15OwvKCIGd2nApfUqMnL2FqKOF5+qoMyyOvu54oorePrpp/nwww+JjY3NdiSw07uKPn1alSpVKF26NH379mX58nPr4G/dunV069aNatWq4fF46N+/P0uXLgWgRIkSGWctBdFNc2HRBKD8VmRsEmPnbONy1xb6uxc4HU6euF3CW7e3JinVy/BZm8/+ggATFhaG2+2mevXqp5Q///zzTJgwgaSkJDp16sT27duzfP3pXUVndnq/SyJySvfNwCndLWcnp+q5oKCgjPcpiG6aC4smAOW3xszZSprPZ1X9iP/VpTeqVpanrruIeVsPM3+rf9Yhn4sjR47w4IMP8uijj56xs969ezctW7bkueeeo3379mzfvp1y5cpl9MKZG/PmzePo0aMZQ0VeccUV1KhRg6ioKGJiYjh58iS//vprxvzZLb9jx44sWbKE6OhovF4vM2bM4KqrAutubm0GqvzSsp1HmPv3IZ7pfhF1l/pP1c/phnRpwA8bIhj96xa6NKlKqaBCGmD+LM0281v6cJDpzUAHDBjA008/fcZ877//PosWLcLtdtOsWTN69uyJy+XC4/HQunVrBg0aRKVKlXJ8ry5dujBgwAB27drF3XffTfv2VqeYI0aMoGPHjjRo0ICLL744Y/5Bgwbx4IMPUrp0aVatWpVRXrNmTV577TWuvvpqjDH06tWLPn365NMWKRrO2h20k7Q7aJWVk2leer6/DJ8x/PFUV0qOybou+JzkpTvo/DAqjpW7orl7whqevK4JT153UYG8jXYHHbjOpztorQJSfufL5XsIi05kVO/mlPQU0hFzAercuCo3tarJZ4t3sy9G+wpShUcTgPIrB2KT+GjBLm5oXoNuTauf/QV+4sUbm+F2CS//usXpUFQxoglA+Qe7bf8bb76Czxheuimw7tS8sEIpnri2CfO3RbF0h/9e01D+RROA8hubfA2Y7buC+69sSJ1KZZwOJ39k6jV00BXB1K1cmlfnbsOrA8qrQqCtgJRfMAbGpvanCnE8sPIqWJXkdEj5rqTHzf9uuJjHZvzFjxsiuL19XadDUgFOzwCUX1joa8sa04wnPD9STgJv55/uplY1aV23Iu/8uYOkFK/T4agApwlAFXlpXh+vpd1FQ4nkLveZfb8EEhFheK9LOHQ8mYkr9jgdTr4REQYMGJDxPC0tjWrVqmV0p5BbwcHBREdH53d4xZYmAFXkfRcSwS5Th/95viFIAv+o+LIGlbm+WQ0+W7yb6ISTToeTLy644AI2b95MUpJ19jZv3jxq167tcFRKE4Aq0pJTvbw/fwft5B9ucBWfmwKf63ExSalePl202+lQ8k3Pnj2ZM2cOYHWzfNddd2VMO3r0KLfccgutWrWiU6dObNq0CYCYmBi6d+9O27ZteeCBB07pn+frr7/msssuo02bNjzwwAN4vYF/cJDfNAGoIu3r1XuJij/Js0Hfclq3MQGtcfWy9G1bm6/X7OVQXEH0Ftotn//O7s477+Sbb74hOTmZTZs20bFjx4xpI0eOpG3btmzatIlXX32Ve++9F4DRo0fTpUsX/vrrL3r37s2+ffsA6+7Xb7/9lhUrVhAaGorb7WbatGl5WH8FmgBUEXYiJY3Pl+zmisZV6OTKulfIQPb4tU3w+QyfLNrldCj5olWrVoSHhzNjxgx69ep1yrTly5dnXCO45ppriImJIS4ujqVLl3LPPfcAcOONN2b0A7RgwQLWr19Phw4daNOmDQsWLCAsLKxwVygAaDNQVWRNXrmX6IQUvri+KXzldDSFr27lMtzRoS7frNvHA1fl970Pi/NxWbnXu3dvnnnmGRYvXkxMTExGeVZ9kqX3FHp6j6Hp8w8cOJDXXnut4IItBvQMQBVJ8cmpfLF0N1c3rUa7+jn3/hiQ7JvDHg3tjXhT+GhBYJwFDB48mBEjRtCyZctTyrt27ZpRhbN48WKqVq1K+fLlTyn/7bffOHbsGADXXnstM2fOzBhR7OjRo+zdu7cQ1yQw6BmAKpImLg8n9kQqT1/f1OlQCk8Ww1jWkqPc7V7A1A0leKhbI4KrZj8Qij+oU6cOTzzxxBnlo0aN4r777qNVq1aUKVOGyZMnA9a1gbvuuotLL72Uq666inr16gHQrFkzxowZQ/fu3fH5fAQFBfHJJ59Qv379Ql0ff6fdQasiJy4plS5vLOTylFWMK/Fe4b55YXcHnQtRpiJXej/nxlY1efeONue0DO0OOnBpd9AqoExZGU58chpPeH50OpQiobrEck+n+vwcGqndRat8pQlAFSmJJ9P4csUerr24Os1dWqebbljXhrhF+GxJ4NwXoJynCUAVKdPW7CX2RCqPXNPY6VCKlBrlS3FHhzrMXL+fg3GB2xeSKlxnTQAiUldEFonINhHZIiJP2OWVRWSeiOy0/1eyy0VEPhSRXSKySUQuzbSsgfb8O0VkYMGtlvJHyalexi3dQ5fGVbm0XjFs+XMWD3RthDHwxRJt767yR27OANKA/xpjLgE6AY+ISDPgeWCBMaYJsMB+DtATaGL/DQM+AythACOBjsBlwMj0pKEUwLfr9hOdcJJH9eg/S3Url+HWtrWZsXYfR+IDo48g5ayzJgBjzEFjzAYQ/RW+AAAgAElEQVT7cTywDagN9AEm27NNBm6xH/cBphjLaqCiiNQEbgDmGWOOGmOOAfOAHvm6NspvpaT5+HzJbjoEV6Jjg3wc5D3APNStEaleHxOW61mAOn95ugYgIsFAW2ANUMMYcxCsJAGkD9BaG9if6WURdll25ae/xzARCRGRkCNHdGi84mLWXwc4GJfMI1c3zvLOT2VpWK0sN7aqxder9hJ3ItXpcHItv7qD7tatG+lNw3v16kVsbGy+xlnc5DoBiEhZ4AfgSWPM8ZxmzaLM5FB+aoEx44wx7Y0x7atVq5bb8JQf8/kMXyzdTbOa5bnqIv3Mz+bBqxqSmOLl6zX+00qqILqDnjt3LhUrVsyP8IqtXCUAEQnC2vlPM8akN84+bFftYP+PsssjgMxj2dUBInMoV8Xc/G2H2X0kkQeuaqhH/7nQvFYFul5Uja9WhJOc6j9dIOfUHXRiYiKDBw+mQ4cOtG3blp9//hmApKQk7rzzTlq1akW/fv0yEgicOjjMLbfcQrt27WjevDnjxo3LmKds2bIMHz6c1q1b06lTJw4fPlwYq+o3ctMKSIAvgW3GmHczTZoNpLfkGQj8nKn8Xrs1UCcgzq4i+gPoLiKV7Iu/3e0yVcx9sTSMOpVKc2PLmk6H4jce7NqQ6IST/LjhwLktoFs+/+VCTt1Bjx07lmuuuYZ169axaNEinn32WRITE/nss88oU6YMmzZtYvjw4axfvz7LZU+cOJH169cTEhLChx9+mNHRXGJiIp06dWLjxo107dqV8ePH5y7YYiI3fQFdAQwA/haRULvsBeB14DsRGQLsA263p80FegG7gBPAfQDGmKMi8gqwzp7vZWPM0XxZC+W31oUfZf3eY4zu3RyPW29LyVbmfoJGxXF5oyq0qlOBcUt3069DXdyuon/mlFN30H/++SezZ8/m7bffBiA5OZl9+/axdOlSHn/88YzXt2rVKstlf/jhh/z0008A7N+/n507d1KlShVKlCiRcZ2hXbt2zJs3r6BWzy+dNQEYY5aTdf09wLVZzG+AR7JZ1kRgYl4CVIHtiyW7qXxBCe5oXzfLztBU1kSEB7o24pHpG/hzyyF65vXsaXGBhHVWOXUH/cMPP9C06Zmd/52tWnDx4sXMnz+fVatWUaZMGbp160ZysjWITlBQUMbr3W43aWlp+bg2/k8PuZRjdhyOZ/62KAZeHkzpEm6nw/E7PVpcSP0qZfh8ye4s+9MvirLrDvqGG27go48+yliPv/76Czi1m+jNmzdnDBWZWVxcHJUqVaJMmTJs376d1atXF/BaBA5NAMox45eGUTrIzb2Xaxe+58LtEu6/siEbI+JYs8c/alOz6w76pZdeIjU1lVatWtGiRQteeuklAB566CESEhJo1aoVb775JpdddtkZr+3RowdpaWm0atWKl156iU6dOhX4egQK7Q5aOSIqPpkury/iTn7n5aBJTofzryLYHfQZRsVlPExO9XL5awtoV78yEwZm3/uvdgcduM6nO2gdEEY5YuqqvaT6fAwO+s3pUPxPpmslpUbFMaBTfT5atIuwIwk0rFbWwcCUv9EqIFXoklK8fL16L9dfUoNgl7bLPi+jKjBgxXUEmRS+fP8lp6NRfkYTgCp0P2yI4NiJVIZe2dDpUAJCNTnOre4V/OC9kqOJKdnOV5Sre9W5Od/PVBOAKlQ+n2Hi8j20rlOBDsHaGWx+GeKeSzIlmbY66+4hSpUqRUxMjCaBAGKMISYmhlKlSp3zMvQagCpUC7dHERadyId3tdVuH/LRRa4DXOUKZfKqkgy7qiElPac2q61Tpw4RERFoB4uBpVSpUtSpU+ecX68JQBWqCcvDqEU0vX5sDj/5nA4noAx1z2VAQht+Do20bqzLJCgoiAYNGjgUmSqqtApIFZotkXGsDjvKIM8feER3/vmti2szTWuU46sV4VrVo3JFE4AqNF+tCKd0kJt+7sVOhxKQROC+K4LZdvA4q8P848Yw5SxNAKpQRCecZHZoJLe1q0MFSXQ6nIB1S9vaVCoTxFcr9jgdivIDmgBUoZi2eh8pXh+Drgh2OpSAVmpsZfqf/I55Ww+yL+aE0+GoIk4TgCpwJ9O8TF29l6ubVqOR3qla4AZ45uHGx6SV4U6Hooo4TQCqwP06ug/RCScZHPaUdvlcCGpILDe5VvNdyH7ik/1n3GBV+DQBqAJljGFiWg+aSARdXJudDqfYGOz5jYSTaXwfEuF0KKoI0wSgClTI3mNsMQ0Y5P4dve+r8LRy7aFd/UpMXhWOz6dNQlXWNAGoAjVpRTgVSOBW9wqnQyl2BnUOZm/MCRbviHI6FFVEaQJQBSYyNonftxziTvciyshJp8Mpdnq0uJALy5fiqxXhToeiiihNAKrAfL16L8YYBnh0IG4nBLldDLi8Pst2RrMrKt7pcFQRpAlAFYjkVC8z1u6je7MLqSPRTodTbN3ZoS4lPC5tEqqypAlAFYifQw9w7ESq3vjlsCplS9KndS1+WH+AuCRtEqpOpQlA5TtjDF+tCOfiC8vRsUFlp8Mp9gZ2DiYp1cv3IfudDkUVMZoAVL5bs+co2w/FM6hzsPb5XwS0qF2By4IrM3lVOF5tEqoy0QSg8t3kleFULBPELW1rOx2Ksg3sHMz+o0ks2q5NQtW/NAGofHUgNok/thyi38kfKDW2snb9UER0b16DC8uXYvKqcKdDUUWIJgCVr762x6Qd4JnvcCQqsyC3i3s61dMmoeoUmgBUvklO9fLN2n1c36yGNv0sCkZVOOXvziXXUIJUJn/wotORqSJCE4DKN7M3RnLsRCoDOwc7HYrKQlU5zs2ulfzg7cpx7SVUoQlA5RNjDJNWhNO0Rjkub1jF6XBUNgZ5/uQEpZipvYQqNAGofBKy9xhbDx5nYMx7yOiKToejstHStYdLZQdTVmkvoSoXCUBEJopIlIhszlQ2SkQOiEio/dcr07T/E5FdIvKPiNyQqbyHXbZLRJ7P/1VRTpq0MpzyJHKL9vpZ5A30/EF4zAmW7DzidCjKYbk5A5gE9Mii/D1jTBv7by6AiDQD7gSa26/5VETcIuIGPgF6As2Au+x5VQA4GJfE75sP0U97/fQLPV1rqVauJJO1f6Bi76wJwBizFDiay+X1Ab4xxpw0xuwBdgGX2X+7jDFhxpgU4Bt7XhUApq3eh88YBri16ac/KCFe7kn6msX/HCFsRFOnw1EOOp9rAI+KyCa7iqiSXVYbyNzhSIRdll258nPpvX5ee3EN6rn0LlN/cZd7IUGkMcXb3elQlIPONQF8BjQC2gAHgXfs8qw6fjE5lJ9BRIaJSIiIhBw5onWURd2cTQeJSUxhkDb99CvVJY4bXauZ6e1Kwsk0p8NRDjmnBGCMOWyM8RpjfMB4rCoesI7s62aatQ4QmUN5VsseZ4xpb4xpX61atXMJTxUSYwyTV4XTuHpZrmisTT/9zSDPHyRQhh/Wa5PQ4uqcEoCI1Mz09FYgvYXQbOBOESkpIg2AJsBaYB3QREQaiEgJrAvFs889bFUU/LU/lk0RcQy8vL72+umH2rh201p26cDxxVhumoHOAFYBTUUkQkSGAG+KyN8isgm4GngKwBizBfgO2Ar8DjxinymkAY8CfwDbgO/seZUfm7wynHIlPfS9tI7ToahzNMjzB2FHElm2S7vuKI48Z5vBGHNXFsVf5jD/WGBsFuVzgbl5ik4VWYePJzNn00EGXF6fC0qe9WukiqgbXasZW/ZpJq8M56qLtMq1uNFfrjon09bsw+vzMnDdrbD+sNPhqHNUQrz071iPDxbsZE90Ig2qXuB0SKoQaVcQKs9OpnmZvmYv17hCCXbpzt/f9e9YjyC3MGVVuNOhqEKmCUDl2ZxNB4lOSGGQ+3enQ1H5oHr5UtzYsibfh0Rok9BiRhOAypP0Ad8bVy9LF9fms79A+YVBVzQg4WSaNgktZjQBqDzZsC+Wvw/EMbBzMNryM3C0qVuRNnUrMnmlNgktTjQBqDyZtDKccqU89NUB3wOHPWLYfYfGEBadyFLtJbTY0ASgcu1QXDK/bdxHv9SfueA1vfM30PR0raF6uZJ8tSLc6VBUIdEEoHJt6upwfLgY6P7T6VBUASghXu7pVJ8lO46wKyrB6XBUIdAEoHIlOdXL9DX7uN4VQl2XVhEEqrs71qOEx8WklXucDkUVAk0AKldm/XWAYydSuc+jTT8DWdWyJenTuhY/rD9A3AkdOD7QaQJQZ2WMYeKKPTSrWZ6Ost3pcFQBu++KBiSlevlm3T6nQ1EFTBOAOqsVu2LYcTiB+67Qpp/FQbNa5enUsDKTV4aT5vU5HY4qQJoA1Fl9tWIPVcuW4ObWtZwORRU0u0no4P3DiYxL5s+t2tVHINMEoHK0JzqRhf9E0b9jfUoFuZ0ORxWSa10bqCeHmbhcLwYHMk0AKkdfrdhDkMtF/071nA5FFSK3GAa5/yBk7zFC98c6HY4qIJoAVLZiT6TwfUgEfVhI9XdqWNUDqti4w72YciU9fKlnAQFLE4DK1vS1+0hK9TLE/ZvToSgHlJVk7upYj7l/H+RAbJLT4agCoAlAZSklzcfkleFc2aQqF7v2Ox2OcsjAzsGANfynCjyaAFSW5vwdyeHjJxnSpYHToSgH1a5Yml4tazJjzT4dKyAAaQJQZzDGMGHZHhpXL6vjxCqGdGlA/Mk0vlunZ4KBRhOAOsPqsKNsiTzOkC4NEL3zq9hrU7ci7etXYuKKPXh1rICAoglAnWH8sjAqX1CCW7XPf2XfGDY0cgQRx5L4ffMhpyNS+UgTgDrFzsPxLNwexcDLg/XGL5XhelcIwVXKMG7pbozRs4BAoQlAnWLc0jBKcZIBy67OOPpTyi2GoVc2ZGNEHGv3HHU6HJVPNAGoDIePJzMr9AB3uJdQWeKdDkcVMbe1q0PlC0owbmmY06GofKIJQGX4akU4Xp9hqHuu06GoIqjU2Mrce3I6C7ZHsfOwHiAEAk0ACoCEkTWYtmQTPWU19VxRToejiqh73fMoxUnGL9OzgECgCUAB8I33GuK5gGGeX50ORRVhlSWe291LmPVXJFHHk50OR50nTQCKlDQfX6b1pKNspbVLj+xUzoa655Lm8/HlCu0kzt9pAlDMCj3AQarwkGe206EoP1DfFUWvljWZtnofcUk6brA/0wRQzHl9hs+X7Ka57OEq1yanw1F+4qFujUg4mcbUVeFOh6LOgyaAYu7PLYcIO5LIQ57ZOt6vyrXmtSrQrWk1vloRTlKK1+lw1Dk6awIQkYkiEiUimzOVVRaReSKy0/5fyS4XEflQRHaJyCYRuTTTawba8+8UkYEFszoqL4wxfLZkN8FVytDTtdbpcJSfebhbY2ISU/guRDuJ81e5OQOYBPQ4rex5YIExpgmwwH4O0BNoYv8NAz4DK2EAI4GOwGXAyPSkoRwyqgIrRnRhU0QcD8R9gFv09n6VN5c1qEz7+pUYtzSMVK/P6XDUOThrAjDGLAVOv/e7DzDZfjwZuCVT+RRjWQ1UFJGawA3APGPMUWPMMWAeZyYVVcg+9famOsfo617mdCjK39jdhDwc+X8ciE1idmik0xGpc3Cu1wBqGGMOAtj/q9vltYHM54MRdll25WcQkWEiEiIiIUeOHDnH8NTZrPc1YaWvBfd75lBSdKAPdW6udoVysezlk8W7tKtoP5TfF4Gzuoxocig/s9CYccaY9saY9tWq6WAkBeWjtFupzHH6uxc4HYryYyLwmGcWYUcSmfv3QafDUXl0rgngsF21g/0/ve+ACKBupvnqAJE5lCsHbNwfy2JfG4Z65lBGTjodjvJzPV1raVK9LB8t3IlPzwL8yrkmgNlAekuegcDPmcrvtVsDdQLi7CqiP4DuIlLJvvjb3S5TDvho4S4qkMC97nlOh6ICgEsMj17TmB2HE/hzqw4Y409y0wx0BrAKaCoiESIyBHgduF5EdgLX288B5gJhwC5gPPAwgDHmKPAKsM7+e9kuU4VsS2Qc87cdZojnN8qK9uWi8sdNPzWnoUTy4bQfdcAYP+I52wzGmLuymXRtFvMa4JFsljMRmJin6FS++3jhLsqV9DBQT8BUPnKL4WHPzzyT+hALtkVxXbMaToekckHvBC5G/jkUz2+bD3HfFcFUkBNOh6MCTB/XSurJYT5cuFPPAvyEJoBi5L15OyhX0sPgLg2cDkUFoCDx8qh7Fpsi4pi/TceU8AeaAIqJzQfi+H3LIYZ4v6Him9q8VhWMvu5lBFcpwzt//qMtgvyAJoBi4t15O6hAAoPdvzkdigpgHvHx5HUXsd2ublRFmyaAYmDDvmMs3B7FMM+vlJckp8NRAe7m1rVoUr0s783foXcHF3GaAALdqAq89/nnVCGOQW5t+aMKntslPHndReyKSuCXjXq/Z1GmCSDArfU1ZZmvFQ96fuECvetXFZKeLS7k4gvL8f78HdpTaBGmCSCAGWN4I/VOqnOMe9zznQ5HFSMul/Df7k0JjznB9yERToejsqEJIIDN23qY9aYpT3lmUlpSnA5HFTPXXVKd9vUr8f78HTpqWBGlCSBApXl9vPnHPzSUSG53L3E6HFWc2GMFyOiKPHfwSaLiTzJxxR6no1JZ0AQQoH7ccIBdUQn8z/MtHtE6WOWMDq5/uO6S6ny+eDfHEvUstKjRBBCAklO9vDtvB23qVuQG1zqnw1HF3LM3XExiShqfLt7ldCjqNJoAAtCkleEcOp7M8z0vRrIaikepQtT0wnL0vbQOk1fuJeKY9kFVlGgCCDDRI+vyyW8buNr1F52mNHQ6HKUAePr6ixCBN3//x+lQVCaaAALMe2m3cYKSDPdMczoUpTLUqliaYV0bMntjJOv3HnM6HGXTBBBA/jkUzwzvNQxwz6OxS+/AVEWE3SrowRVXUr1cSV75dat2FFdEaAIIEMYYxszZSjlO8ITnR6fDUeoMF8hJnr2hKaH7Y/llkx6gFAWaAALBqAosGtGNZTujecLzI5UkwemIlMrSfy6tQ4va5Xn9t+16c1gRoAkgAKQYN2PS7qGhRDJAB3pXRZjLJbx0YzMOxiXz+ZLdTodT7GkCCAATvL0IM7V4yfM1QaJHVapo69iwCje1qsnnS3azL0abhTpJE4CfOxCbxEdpt9LdtY6r3aFOh6NUzuwLwi/+0xe3S3j51y1OR1SsaQLwc2N+3YpBGBE01elQlMq1C+UYT17XhPnbopi/9bDT4RRbmgD82JIdR/ht8yEe88yijkQ7HY5SeXLfFQ1oUr0so3/dQnKqVl06QROAnzqZ5mXU7C00qHoBQ91znA5HqTwLcrsY3ac5+48m8eki7SfICZoA/NTHC3exJzqR0b2bU1LSnA5HqbwbVYHOUxvRx7WCzxZuZ+fheKcjKnY0Afih7YeO89ni3fRtW5uuF1VzOhylzstLQVO5gCSe+2GTDiJfyDQB+Bmvz/DcD39TvnQQL97UzOlwlDpvVeU4I4KmsmFfLF+v3ut0OMWKJgA/M2llOBv3xzIy5R0qv1XNalanlJ+71bWcK5tU5c3ft3MgNsnpcIoNTQB+ZP/RE7z9xz9c3bQavV2rnA5HqXwjAq/e2hIDDP/pb4zRqqDCoAnAT/h8hme+34jbJYy5taUO9KICTt0Pa/KsbyKL/znC9yERTodTLGgC8BMTV+xhzZ6jjLi5GbUrlnY6HKUKxED3n3RybWH0L1vYf1S7iShomgD8wI7D8bz5xz9cd0kNbm9Xx+lwlCowLjG8HfQFIsJ/v9+orYIK2HklABEJF5G/RSRURELsssoiMk9Edtr/K9nlIiIfisguEdkkIpfmxwoEupQ0H099G0q5tGO8FnYrMrqiXvhVAa2ORDPS+zFr9xxl4oi7nQ4noOXHGcDVxpg2xpj29vPngQXGmCbAAvs5QE+gif03DPgsH9474H2wYAdbIo8zNuhLqslxp8NRqlDc5l5Kd9c63krrx7aD+r0vKAVRBdQHmGw/ngzckql8irGsBiqKSM0CeP+AsXxnNJ8u3s3t7erQwx3idDhKFRoReDXoSyqQwKPTN5B4Uu92LwjnmwAM8KeIrBeRYXZZDWPMQQD7f3W7vDawP9NrI+yyU4jIMBEJEZGQI0eOnGd4/isqPpknvw2lUbWyjO7T3OlwlCp0VeU4HwR9Qlh0IiN+1m6jC8L5JoArjDGXYlXvPCIiXXOYN6uGi2dc4THGjDPGtDfGtK9WrXh2c+D1GZ76NpT45FQ+uftSypTwOB2SUo7o7N7KY9c04YcNEcxcr01D89t5JQBjTKT9Pwr4CbgMOJxetWP/j7JnjwDqZnp5HUBHhs7Cp4t2sWJXDKN7N6fpheWcDkcpRz1xbRM6NqjMS7M2sytKO4zLT+ecAETkAhEpl/4Y6A5sBmYDA+3ZBgI/249nA/farYE6AXHpVUXqX4v/ieLd+Tvo06YW/ea0zBhBSaniyv1yRT6MvJMyqUd54L0ZxCenOh1SwDifM4AawHIR2QisBeYYY34HXgeuF5GdwPX2c4C5QBiwCxgPPHwe7x2QwqMTeXzGXzStUY7X+urdvkqlqyGxfBz0IeHmQp76diM+vT8gX5xz5bIxJgxonUV5DHBtFuUGeORc3y/QJZ5MY9jUEFwuYdyA9lrvr9RpLndv4yUzlVHbBvHBS0N4KugHGBXndFh+TfcyRUB6Pz+7ohKY4nmVeh9tdjokpYqkge4/2Wwa8IH3PzRz7eUGpwPyc9oVRBHw9p//8NvmQ/xfz0vo4tadv1LZEYExnom0ll08mfowmyJinQ7Jr2kCcNi0NXv5dPFu7u5Yj6FXNnA6HKWKvFKSyvgS71CZeAZPCtFO486DJgAHLdoexUuzNnN102q83Ls5old9lcqV6hLH5BJvkJLmZdBXa4k7oS2DzoUmAIdsiojlkekbaMYePg6/Cc8rlbS5p1J50NgVyTjfCPYfieX+Vz4gOdXrdEh+RxOAA7YfOs69E9dS+YISTCzxFhfISadDUsovdXJt5+2gz1lnmvLQ1+tJSfM5HZJf0QRQyMKOJHDPhDWU8riZPrQT1UUvYil1Pnq7VzHWM5FF/xzhiW/+Is2rSSC3NAEUov1HT9B/whqMga+HdqRelTJOh6RUQLjbs5AXb7yE3zYf4n8zN+mNYrmk9wEUkj3RifQfv5rEk2l843uGxp/uczokpQLK0AVtSfLcwjt/3YHZ+A1vvfIqHrce4+ZEt04h2HE4nju+WEVymo/p93eimUt3/koVhMc8s3jG8y0/+a7k0el/6TWBs9AEUMA2H4ij3xerEODbYZ1oUVtb+ihVkB71/MxLnin8vuUQw0a8RvLIqk6HVGRpFVBBGVWBJd5WPJz6BBVJYPqzt1O/ygVOR6VUsTDE8ztlOMkLaUO4O2U44xNOUqVsSafDKnL0DKCATE+7hsGpz1JPDjOz5Gjd+StVyO7yLOLToA/YYoLpO/ZrwkY01XttTqMJIJ/5fIbXf9vOC2lD6eL6m+9LvExNOep0WEoVSz3d65hRYgzxpgx9U0azxnex0yEVKZoA8lHsiRTum7SOz5fs5m73fL4MepuykmxNTB/YRY9AlCpUl7p28VOJEVSWePqnvMBXK/Zg9U6vNAHkk80H4rjpo+Ws3B3N2FtbMNYzEY9oCwSlioL6rih+KjGCbq6NjP5lK09+G8qJlDSnw3KcXgQ+T8YYvl6zjzG/bqXyBSX47oHLaVuvEvzmdGRKqcwqyAnGBb3LZ97evB16O9s2ruGDoI+5xLX/35mK2QAzegZwHqLikxk8aR0vzdpMp4ZV+DV5IG0nBms1j1JFlEsMj3h+ZkrQ6xwzZemTMobxab3wmeLZE68mgHNgjGHu3wfp8f4yVu6O4eU+zZl0XweqSLzToSmlcuFK92Z+L/k83VyhjE27h/6pL7DPV93psAqdVgHl0YHYJEb+vJn526JoUbs87/drQ+Pq5ZwOSymVR1Ukni+C3uM7bzdeThtA95Q3eOLFexjqnkuQeItFdZAmgFxKSfMxeWU4783fgTEwvNcl3HdFsPY1opQfE4F+nsV0dW9iZOpA3ki7i5+9nXklaBIdTq/KDcCEoAngLIwx/L75EK//vp29MSe4xrWB0Z5J1F0YDQudjk4plR9qylHGlXiPP7ztGZU6kNtTRtLDtZbnPTMIdh12OrwCowkgG8YYVoXF8N68HawLP8ZFNcoy6b4OdJtxt9OhKaUKyA3uELq6NjHeeyOfp93MgpRL6e+ez0Oe2dRwOrgCoAngNMYYVu2O4f35O1kbfpTq5Uoy1jOBfrGL8czQdv1KBbrSksLjnp+4072Id9NuY6r3eqZ7r+GunzfzULfGXFih1Kkt/fy4akgTgC0lzcecvyOZuDycvw/EcWH5Uozu3Zx+HepSaux/nA5PKVXIqkssrwdN4GH3bD7x9mHampJMX7uPm1vVYrAvmBaucKdDPG/FPgHsP3qC79dH8M2CtURRiUZygLGe3/jPszMpFeR2OjyllMPquaJ4wzWeR32z+NLbk+//uoofeZXLZBt3eRbSM9Xrt/sKKcp9YrRv396EhITk+3ITTqYxf+thZq6PYMXuaACulI0Mdv9GV9ffuKTobhNVwAYlWv8n2b23TvrVLr/JmXhUkXPclOY779VM8V7PPlODciTS272KW93LuVR24hrt/DjfIrLeGNP+bPMVmzOAuKRUlu44wpxNB1n0TxQn03zUrliaJ65twm3t6lDnA724q5Q6u/KSxFDPXAa7f2O17xK+83Zjprcr07zXUYtoev26lZ4tL6RN3Uq4XUX7DuOATQA+n2HboeMs3xnNwu1RhOw9htdnqMYx7nKv4cYSa2iXtAPXcgPLnY5WKeVvXGLo7N5KZ/dWXjFfMd93KXO8nZi8vCITlu+hEvFc1aYp3ZpWp3OjKlQvX8rpkM8QMAngZJqXrZHH2bAvlpDwo6wKiyH2RCoAF8s+HnD9xTUl/qKt7MStVTxKqXxUTpK41b2CW6crue4AAAeBSURBVN0riDNlWOprxSJvWxaFwqzQSAAaVy/L5Q2r0K5+JS6tV4m6lUsj4uwZgl8mgISTaew8HM/Wg8fZEmn9bTt4PGMA6Noc4Tr3VjoHbeFy11YdkEUpVWgqyAludq/mZvdqvEbYYoJZ5WvOyuhm/BjVlKmrSwNQtWwJmteqQPNa5WleqwKX1CxHvcplCrV3gSKdAFK9PhZsO0x4zAn2xiSyJzqR3VEJRMYlZ8xTrpSH5rXKM6hzMJfWq0jbepWo8W4g3rKhlPI3bjG0kj20cu3hAX7Fa4R/TF3+8jXmr6QmbN4ZzIodtUmzd8UlSKWBHKSRHKS+HCJYDlN/8ERa16lI6RL539Ko0BOAiPQAPgDcwARjzOvZzRuTmMKQyVYroHIlPdRP2UlHiaCx5wCN73qLZjXLU6dSaWR0RYgE1hbKKiil1Dlxi6GZ7KOZax/97b5kkk0QO00d/jF12OmrzS5Th62mPn/62lmJYdxq5pd4hsauyHy/6axQE4CIuIFPgOuBCGCdiMw2xmzNav7KCbuYVGIEwXKYSsQjJTNN/H52IUSslFIFq5Sk0lL20JI91mGxLc24iDRV2WuqU0/s/ojyeayRwj4DuAzYZYwJAxCRb4A+QJYJoASpXOraVYjhKaVU0eARH/UkinpEFdh7FOqNYCJyG9DDGDPUfj4A6GiMeTTTPMOAYfbTpkAMEF1oQRaeqgTeegXiOkFgrlcgrhME5nqdyzrVN8ZUO9tMhX0GkFWbp1MykDFmHDAu4wUiIbm5o83fBOJ6BeI6QWCuVyCuEwTmehXkOhX2aCYRQN1Mz+tgXb5VSilVyAo7AawDmohIAxEpwf+3d3YhWpRRHP/9MxSsrNWwIgrsg0xCQyskTOkDvy7aTAojSrKbokAFwRWjGwmyj5uuIlSym5LAyIgyE8pAV1TwYyV1dytQEoWyIgJt4XTxnLcdlnd293VXZ2bf84NhZp85s8yPM7xnZp6ZZ2AJEL25QRAEBXBZbwGZWY+kV4HtpP7uTWZ2dIDNPhhgfVUZiV4j0QlGptdIdIKR6XXJnEo9GmgQBEFw6YgvmgdBEDQpUQCCIAialMILgKTlkjokHZW0wtu2SDro0y+SDmbi10jqknRc0rzi9rx/crzuldTuXvslPeDtkvSeex2WNL3Yva9PjtM0SXskHZH0haRxmfhS5krSJklnJXVk2sZL2iGp0+ct3p6bG0lLPb5T0tIiXLI06DXZ83Ze0qo+/2e+56xLUtvl9uizL404Pes5Oixpt6RpmW1K4+T704hXqzvVfjdmZbYZ2jFoZoVNwD1ABzCW1CH9LXBnn5h3gdd9eQpwCBgDTAK6gVFFOjTiBXwDLPCYhcB3meWvSO9JzAT2Fu3QgNM+YI7HLAPWlT1XwGxgOtCRaXsLaPPlNmB9f7kBxgM/+bzFl1sq5DURuB94A1iViR/luboNGO05nFIRpwdrOQAWZHJVKqeL8Lqa3v7aqcCx4ToGi74CuBtoN7N/zKwH+B5YVFspScDTwMfe1Ap8YmbnzexnoIs0vETZyPMyoHaGfC2970C0Ah9Zoh24TtJNl3unByDP6S5gl8fsABb7cmlzZWa7gL5jhLcCm315M/BEpr1ebuYBO8zsdzM7R3Kff+n3Pp9GvMzsrJntA/7tE///cC1mdgGoDddSCA067fZcALST3jOCkjlBw15/m//iA1fR+/LskI/BogtABzBb0gRJY0lnW9kXxR4CzphZp/99M3Ays/6Ut5WNPK8VwNuSTgLvAGs8vgpeeU4dwOMe8xS9+auCU5YbzOw0gM8nenueR1X88rzyqILXYJxeJF25QTWcoB8vSYskHQO+JF1pwzB4FVoAzOxHYD2pcn1NujTryYQ8Q+/ZPwxiKIky0I/Xy8BKM7sFWAls9E1K79WP0zLgFUkHgGuAC75J6Z0GSZ7HSPHrS+W9JD1MKgCra011wirlZGafmdlk0lXBOm8eslfRVwCY2UYzm25ms0mXRJ0Akq4EngS2ZMIrM5REjtdSYKuHfErvLZFKeNVzMrNjZjbXzGaQinW3h1fCKcOZ2m03n9eGYMzzqIpfnlceVfDKdZI0FdgAtJrZb95cBScYRK781tHtkq5nGLwKLwCSJvr8VtIPfu2M/zFSZ8epTPg2YImkMZImkTohS/kZmByvX4E5HvIIXuxIXs/7EyczgT9rl4Jlop5Tpu0K4DXgfQ+vTK6cbaQCjc8/z7TXy812YK6kFn9aY663lY08rzyqMFxLXSc/LrcCz5nZiUx8FZwg3+sO7w/Fn0IbTRoleejHYJE94d6v8QPpewCHgEcz7R8CL9WJX0s6yzyOP1FTxqmeFzALOOBte4EZ3i7Sh3K6gSPAfUXvfwNOy4ETPr2JP61Q5lyRivFpUgfoKdLtggnATlJR3gmMHyg3pNtfXT69UDGvGz3mL+APXx7n6xZ6PruBtRVy2gCcAw76tD/zf0rjdBFeq4Gj7rQHmDVcx2AMBREEQdCkFH4LKAiCICiGKABBEARNShSAIAiCJiUKQBAEQZMSBSAIgqBJiQIQBEHQpEQBCIIgaFL+AyQei39WsxKfAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def plotdis(rel_err):\n", " plx = 1. # in mas; this value doesn't really matter, so we can keep as such\n", " print(f'Parallax= {plx}, with a relative error of {rel_err}%')\n", " \n", " error_plx = rel_err / 100. * plx\n", " test_par= np.random.normal(plx,error_plx,100000) # Compute a distribution of parallaxes\n", " test_dis = 1000./test_par # Distances is 1000 / parallax in pc\n", " test_dis = test_dis [test_dis > 0.] # Remove the non-physical negative distances\n", "\n", " #print the miminum and maximum distances \n", " print (\"Distances vary between \",np.round(test_dis.min(),1),\"and\",np.round(test_dis.max(),1),\" pc\")\n", "\n", " # Compare to a Gaussian distribution of distances, centred on 1/plx \n", " # and with an error given by \n", " # error_dis = error_plx / plx**2. \n", " dis = 1./plx\n", " error_dis = error_plx *dis**2. \n", " \n", " # this is in kpc, so multiply by 1000.\n", " dis, error_dis = 1000.*dis, 1000.*error_dis\n", " \n", " xdis = np.linspace(-3.*error_dis,3.*error_dis,100)\n", " ydis= np.exp(-0.5*(xdis/error_dis)**2)\n", " \n", " \n", " #Compute the histogramme of distances\n", " bin0 = np.percentile(test_dis,.1)\n", " bin1 = np.percentile(test_dis,99.9)\n", " if rel_err > 20:\n", " bin0 = 0.\n", " bin1 = 3000./plx\n", " \n", " z, bin = np.histogram(test_dis,bins=100,range=(bin0,bin1)) #dis-3.*error_dis,dis+5.*error_dis))\n", " mo = bin[np.argmax(z)] # compute the maximum of the histogram, i.e. the mode\n", " plt.vlines(mo,0,z.max(),color='yellow',label='Mode') \n", " med = np.percentile(test_dis,50) # the median value\n", " plt.vlines(med,0,z.max(),color='magenta',label='Median')\n", " \n", " ydis = ydis * z.max()\n", " plt.plot(xdis+dis,ydis, label='Formal Gaussian')\n", " \n", " plt.hist(test_dis,bins=bin,label='Distribution')\n", " print(\"Formal distance: \",np.round(dis,0),\"+/-\",error_dis,\" pc\")\n", " print(\"Mean distance: \",np.round(np.mean(test_dis),0),\"+/-\",np.round(test_dis.std(),0),\" pc\")\n", " print (\"Median of distance distribution: \",np.round(med,0),\" pc\")\n", " print (\"Mode of distance distribution: \",np.round(mo,0),\" pc\")\n", " \n", " plt.xlim(bin0,bin1)\n", " plt.title(f'Distribution of distances with plx={plx} and error={rel_err}%')\n", " plt.legend();\n", " #print('------------------------------------------')\n", "\n", "\n", "relative_error_plx = 1. # in percents\n", "plotdis(relative_error_plx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With a small error, we can see that we have a Gaussian distribution still. Change progressively the value of the relative error on the parallax (from, say, 5% to 50%) below to see how the distribution becomes more and more skewed and it is no more possible to use the inverse of the parallax to define a distance. In particular, when the error is above 33%, it is not improbable that the distance becomes infinite!" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parallax= 1.0, with a relative error of 10.0%\n", "Distances vary between 701.3 and 1799.1 pc\n", "Formal distance: 1000.0 +/- 100.0 pc\n", "Mean distance: 1010.0 +/- 105.0 pc\n", "Median of distance distribution: 1000.0 pc\n", "Mode of distance distribution: 968.0 pc\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd4FNX6wPHvu7tp9BYQSEIXpCMBoiCCKCIq2AERwYJd79Wf91ovYrvYrt2LInIBRVBRERVEVEJRBALSlA4BAkivIXX3/P6YSVxCOkkmu3k/z7NPsmfaO7O7886cOXNGjDEopZSqeFxOB6CUUsoZmgCUUqqC0gSglFIVlCYApZSqoDQBKKVUBaUJQCmlKqigSAAi8q6I/KuE5hUjIidExG2/jxeR20ti3vb8ZovI8JKaXxGW+5yIHBCRPws5vhGR5vb/JbZ9A5H9fWiaz/BEEbm4hJZVot+3QCEije3vnMfpWCqScp8A7B9XiogcF5EjIvKLiNwlItmxG2PuMsY8W8h55ftDNcbsMMZUMcZ4SyD20SLyUY75X2aMmXSm8y5iHNHA/wGtjTFnFXX6kty+gcj+PmwFEJGJIvKc0zEVloiEish0+7MxItKrgPFriciXIpIsIttF5MYyCjVgFbSNxfKiiBy0Xy+JiOQzvxvtbZ8sIjNEpJbfsNdF5LCILBaRhn7lQ0XkjaLGXu4TgO1KY0xVoBHwAvAI8EFJLySIjz4aAQeNMfucDkQ5YhFwE1CYs793gHSgHjAUGCsibUoxthKX2++4qL/tYuwL8tvGdwBXAR2A9sAVwJ15LLcN8B4wDOszOAn81x7WFegMnGUv7zG7vDrwMDCqiDGDMaZcv4BE4OIcZV0BH9DWfj8ReM7+vw7wDXAEOAQsxEp0H9rTpAAngH8CjQED3AbsABb4lXns+cUDY4ClwFHgK6CWPawXkJRbvEA/rB9Shr28VX7zu93+3wU8CWwH9gGTger2sKw4htuxHQCeyGc7Vben32/P70l7/hfb6+yz45iYx/T/APYAu4Fb7WU3P5Pta4//GdaP4qi9fdv4LXMi1g7nW+A4sARo5je8DTDXXs5e4HG/7fYosAU4CHzq95mEAx/Z5UeAZUC9XNb3FuBrv/ebgU/93u8EOtr/G6A51g85w/5cT2RNb3/mDwOr7fX8BAjPYzuPAH4G3rLHXQ/08Rvu//0YC0z3G/Yi8CMgxfwtJQG98hle2V63s/3KPgReyGP8rsBiezvvAd4GQv2GG+AuYBNw2P6sxR7mBl7B+l5vBe7F73eXy7IaAJ9jfb+3AQ/4DRsNTLc/92PA7XmUhQGvY33Hd9v/h/n/lrEOLv8EPiypbQz8Atzh9/424Nc8pv838LHf+2b2Z1IVGASMscv7AbPs/98GbixWvMWZqCxf5JIA7PIdwN32/xP5awc1BngXCLFfF/h96U6ZF3/tZCfbX/4Ick8Au4C29jifAx/5f2nyitf+En6UY3g8f/3Ab8Xa8TQFqgBfZH3x/OJ4346rA5AGnJPHdpqMlZyq2tNuBG7LK84c0/bD2sFmrePH5J0ACr19/daxKn/9+Fb6DZuItXPvCniAKcA0e1hVrJ3K/2Ht1KsC3exhfwd+BaLs+b4HTLWH3Ql8DVTC2sl0Bqrlss5NsXZcLqA+VtLc5TfsMOCy3+e6LXJ85kuxdlK1gHXAXXls6xFAJvCgvf0GYSWCrATm//2oZH+OI+ztfACIsofF2PHn9Tpth0DBCaATkJKj7GH8EmWOYZ2BOPuza2yv99/9hhusg4Uadrz7gX72sLuwkl+0vc3mkUcCsD+j5VhHuKH257MVuNTvd5aBdZTtwvq95Fb2jP29qQtEYu2Yn/X7jWRiJdkwe/wS2cb259vN730scDyPbfoV8EiOshP2tm6LdeQfAbxsv2KBucXdvwZKFVBudmN9cXLKwPpBNzLGZBhjFhp7K+ZjtDEm2RiTksfwD40xa40xycC/gBuyLhKfoaHAq8aYrcaYE1indINznH4+bYxJMcasAlZhJYJT2LEMAh4zxhw3xiQC/8E6jSyMG4D/+a3j6HzGLdL2NcZMsGNKs+fbwT5lzfKFMWapMSYTKwF0tMuvAP40xvzHGJNqz2OJPexOrLOhJL/5XmdvtwygNtYO22uMWW6MOZZLXFuxzjo6AhcCc4BdItLKfr/QGOPLZzvk9KYxZrcx5hBWAuqYz7j7gNft7fcJsAG4PJcYT2JVK7yKdSR7vzEmyR62wxhTI5/Xx0WIPUsVrJ2Vv6NYyfc09rb91RiTaX/n3sPadv5eMMYcMcbswNrJZ22XG+xtsNPeZmPyiasLEGmMecYYk25/du8Dg/3GWWyMmWGM8fn9jnOWDQWeMcbsM8bsB57m1N+ID3jKGJNm/+ZKahvn3K5HgSp5XAfI8zMwxqzFOgD9FSs5vQi8ATwgIg+IyAIRmSIiNQoZV0AngIZYR485vYx1VP29iGwVkUcLMa+dRRi+HevIrU6hosxfA3t+/vP2YNX9ZfGvUzyJ9QXJqQ7WkVHOeTXMZdy84si5jnkp9PYVEbeIvCAiW0TkGNaRcla8WfJav2isKp7cNAK+tBsFHME68vRibbcPsXbm00Rkt33BLSSP+czHOvLraf8fj7UDu9B+XxSF+Zyy7MqRNLdjfQanMcYsxTraFayqrtJ0AqiWo6waVqI8jYicLSLfiMif9uf7b07/XeS1XYrynWsENMj6vO3P/HFO/Z3k9hvOWZbb781/u+83xqTmE0dx5dyu1YATeRw45fsZGGNeM8Z0MMYMwjroy6qCvQPog/VbKMw+DwjQBCAiXbB2botyDrOPFP/PGNMUuBJ4SET6ZA3OY5YFnSFE+/0fg3WUeQBIxjpNz4rLjXVqWdj57sb6cvvPOxOrOqYoDtgx5ZzXrkJOv4fT1zFXRdy+NwIDsa5DVMeqJgBrZ1aQnVj1n3kNuyzH0Vi4MWaXfVT9tDGmNXA+1pnEzXnMJysBXGD/P5+CE0BBn2lhNMxx9BeD9V04jYjci1UlsRvrulVWeVZz5bxeQ4sR10bAIyIt/Mo6AL/nMf5YrGqcFsaYalg75cJ8tlCE7xzW570tx+dd1RjT32+c3D6XnGW5/d78t/sp45fgNv6dU8/c89ump4xrNz8Ow/ps/GOrh3Um/AxW1dBqY0wG1jWv9oWMK7ASgIhUE5ErgGlYdetrchnnChFpbv/AjmEdGWY16dyLVX9YVDeJSGsRqYS1wacbq5noRiBcRC63jzKfxPqwsuwFGvs3Wc1hKvCgiDQRkSpYR1Cf2NUhhWbH8inwvIhUFZFGwENY1QaF8Skwwm8dn8prxCJu36pY1y0OYiXKfxdhtb4BzhKRv4tImL1e3exh72KtayM7pkgRGWj/31tE2tnJ+BhWYsyrSe98oDcQYVetLMS6HlIb+C2PaYr7HfJXF+u0PURErgfOAWblHElEzgaew6oGGgb8U0Q6winNlfN6TfGbT5iIhNtvQ0UkPLfqB7v67wvgGRGpLCLdsRL4h3msR1WsbXzCrjq7uwjb4FN7G0SJSE3yP2pdChwTkUdEJMI+s2xrHwgWxVTgSfv7UgfrmkKev5ES3MaTsQ6UGopIA6zrWhPzWOwU4EoRuUBEKmPtb74wxuQ8C3sVq7rqJNZF8S72PqQX1hljoQRKAvhaRI5jHQk8gbXyt+QxbgvgB6xTqcXAf40x8fawMVhfgCMi8nARlv8h1gf2J9YFyQcAjDFHgXuA8VhH28lYF4GyfGb/PSgiK3KZ7wR73guwPsRU4P4ixOXvfnv5W7HOjD62518gY8xsrAu0P2FV7/yUz+hF2b6TsS+uAn9g1V0Wiv2FvwTrLONPrJYkve3BbwAzsaqhjtvzzUoOZ2G1/jiGdTo8nzx+5MaYjfZ6LLTfH8Pafj+bvO8D+QBoba/jjMKuTw5LsLbjAeB54DpjzEH/EcS6nvER8KIxZpUxZhPWEfaHIhKWc4YF2IDVOqshVvVYCvaRsIg8LiKz/ca9B+si4z6sHebdxpi8jlYfxjrLO45VJ/9JEWJ6345lFbACK/Hkyv4srsS6frANa7uNxzqrLIrngASs1lpr7OWW1D0deW5jrGsjX9vLXIvV6u29rAnts4kLAOxtfRdWItiHlWTv8V+QiPQGahhjvrSnWWrPcyfWb+SFwgad1XpDKVUGRGQEViufHk7HolSgnAEopZQqYZoAlFKqgiowAdgXM5aKyCoR+V1EnrbLm4jIEhHZJCKfiEioXR5mv99sD2/sN6/H7PINInJpaa2UUuWVMWaiVv+o8qIwZwBpwEXGmA5YF2H6iUgc1k0IrxljWmDdNXmbPf5twGFjTHPgNXs8RKQ11o0bbbBaWvxXSuZmKqWUUsVQYIdH9s0KJ+y3Wbf/G+AirBYAAJOw7sYci9VsbLRdPh14224ONRDrNv80YJuIbOavvkRyVadOHdO4ceMirZBSSlV0y5cvP2CMiSxovEL1eGcfqS/H6hDrHaw7NI/4tVdP4q+7Thti34FnjMkUkaNY7aobcmozQP9p/Jd1B9ZdbcTExJCQkFCYEJVSStlEJL87q7MV6iKwsfpU6YjV+VZXrBtXThsta9l5DMurPOeyxhljYo0xsZGRBSYwpZRSxVSkVkDGmCNY/aXEATXkr07Lovjrluok7Fu87eHVsfrsyS7PZRqllFJlrDCtgCLF7l1ORCKw+nVZh9Wz33X2aMOxujEF6w7N4fb/1wE/2dcRZmL1dBkmIk2w7oRcWlIropRSqmgKcw2gPjDJvg7gwnpoxjci8gdWj4vPYfWbkvWErg+wblffjHXkPxisW5xF5FOsLgEygXvzud1eKVWCMjIySEpKIjW1NDq7VE4JDw8nKiqKkJC8OrzNX7nuCiI2NtboRWClzty2bduoWrUqtWvXJpd+4FQAMsZw8OBBjh8/TpMmTU4ZJiLLjTGxBc1D7wRWqgJITU3VnX+QERFq1659Rmd1mgCUqiB05x98zvQz1QSglFIVVKFuBFOqqE6kZbJo0wF+3XqQxVsOsvvIX49bjqwaRremtTmvWW16tqhDjUqhDkaqyorb7aZdu3bZ72fMmEFZ3ek/ceJEEhISePvtt08b9t133zFq1CiOHTtGeHg4LVu25OWXXyYmJr+HlJ25UaNG0bNnTy6++OJSXU5+NAGoIugFiYtgYmXr7eicz662dvyTfklk3IKtHE3JICLETWzjmpzfvDYyXTAYdvQ+yTerdjN16Q4qhboZcX5jRl7QlJqVNREEs4iICFauXFnk6TIzM/F4SmdXtXbtWu6//35mzpzJOedY97fOnDmTxMTEUk8AzzzzTKnOvzA0AagSYYzh46U7eGXOBg6fzKBPq7rcfkFTOjeqSajHrmn8jz3yCMj0+liz6yj/+zmRsfO3MHnxdu7p3Yw7ezbD7dK66ooiNTWVu+++m4SEBDweD6+++iq9e/dm4sSJfPvtt6SmppKcnMyoUaN46qmnqFevHitXruSaa66hXbt2vPHGG6SkpDBjxgyaNWvG119/zXPPPUd6ejq1a9dmypQp1KtXL8/lv/jiizz++OPZO3+AAQMGZP///vvvM27cONLT02nevDkffvghlSpVYsSIEVxxxRVcd511K1SVKlU4ceIEe/bsYdCgQRw7dozMzEzGjh3L+eefz2233UZCQgIiwq233sqDDz54yjyeeeYZvv76a1JSUjj//PN57733EBF69epFt27dmDdvHkeOHOGDDz7gggsuKLHtrwlAnbFjqRk89vkavl2zh/Oa1uaRy1rRMbpGvtN43C46xdSkU0xN7ruoOS/P2cBL323gl80HeX1wR+pUKepTD1VhPf317/yx+1iJzrN1g2o8dWWbfMdJSUmhY8eOADRp0oQvv/ySd955B4A1a9awfv16+vbty8aN1vPPFy9ezOrVq6lVqxbx8fGsWrWKdevWUatWLZo2bcrtt9/O0qVLeeONN3jrrbd4/fXX6dGjB7/++isiwvjx43nppZf4z3/+k2dMv//+Ow8/nPfTYa+55hpGjhwJwJNPPskHH3zA/ffn/dTWjz/+mEsvvZQnnngCr9fLyZMnWblyJbt27WLt2rUAHDly5LTp7rvvPkaNGgXAsGHD+Oabb7jyyisB6wxo6dKlzJo1i6effpoffvghz+UXlSYAVXyjq7POF83dGQ+y00TySL823NmzKa4iHsGfXa8q44Z1ZtqynTw183f6v7GQd4aeS5fGtUopcOWE3KqAFi1alL1DbdWqFY0aNcpOAJdccgm1av31HejSpQv169cHoFmzZvTt2xeAdu3aMW/ePACSkpIYNGgQe/bsIT09/bT28fk5ePAgffr04eTJk9xxxx08/PDDrF27lieffJIjR45w4sQJLr00/8eYdOnShVtvvZWMjAyuuuoqOnbsSNOmTdm6dSv3338/l19+eXbc/ubNm8dLL73EyZMnOXToEG3atMlOANdccw0AnTt3JjExsdDrUxiaAFSxrfY14ab0x4kgjWmhz9GlV/F79hARhnSNoUNUDe79eAU3jV/C+OGxXNBCOwQsaQUdqZel/G5ErVy58invw8L+Oit0uVzZ710uF5mZVsfE999/Pw899BADBgwgPj6e0aNH57v8Nm3asGLFCjp06EDt2rVZuXIlr7zyCidOWD3gjxgxghkzZtChQwcmTpxIfHw8AB6PB5/Pl70O6enpAPTs2ZMFCxbw7bffMmzYMP7xj39w8803s2rVKubMmcM777zDp59+yoQJE7JjSE1N5Z577iEhIYHo6GhGjx59Stv+rPV0u93Z61lStBmoKpaVvmYMTX+capLM52Gj6eLaUCLzbd2gGtPvOo8mdSpz26QE5m/cXyLzVeVTz549mTJlCgAbN25kx44dtGzZstjzO3r0KA0bWr3MT5o0qcDx//nPf/L888+zbt267LKTJ09m/3/8+HHq169PRkZGdpwAjRs3Zvny5QB89dVXZGRkALB9+3bq1q3LyJEjue2221ixYgUHDhzA5/Nx7bXX8uyzz7JixYpTYsja2depU4cTJ04wffr0Yq590ekZgCqylb5mDEt/jJpynKmhz9FQDpbo/GtXCePjkXHcNH4JIycnMG5YZ3q1rFuiy1Dlwz333MNdd91Fu3bt8Hg8TJw48ZQj/aIaPXo0119/PQ0bNiQuLo5t27blO37WheSbb76Z48ePU7t2bWJiYnj66acBePbZZ+nWrRuNGjWiXbt2HD9+HICRI0cycOBAunbtSp8+fbLPVuLj43n55ZcJCQmhSpUqTJ48mV27dnHLLbdknzGMGTPmlBhq1KjByJEjadeuHY0bN6ZLly7FXv+i0r6AVBH0Ytea9Qyc8iqVJJVpoc/SQA4VPFlWc9Fe9vv4wi3tcHI6Q8cvYfvBZD6/53xanVWtGDErgHXr1p3S0kUFj9w+W+0LSJW4k+khjPzuCdIIYULIy4Xb+edmdPW/XvmoWTmUCSO6UDnMw8jJCRxKTi/e8pRSudIEoArFGMM/PuvHuoONeTPkLZq7yuZZPmdVD+e9YZ3ZeyyNe6YsJ8PrK5PlKlURaAJQhTJ2/ha+XdOKx+Im0du9qkyX3SmmJi9e245ftx7i+W/XFTyBUqpQNAGoAq1OOsKr32/k8vbrGdlhhiMxXN0pilu6N2biL4nEb9jnSAxKBRtNACpfKele/v7JSiKrhvHvq77HyR6FH+nXipb1qvKP6av1eoBSJUATgMrXmNnr2Lo/mVeu70D1SmmOxhIe4ua1QR05cjKdx79Yk+9NREqpgul9ACpP8zfuZ/Li7dzavQndm9cp/oyyWvskfnPGMbVuUI3/69uSF2av5/MVu7iuc9QZz7NCKqAFVtHnd3rPsP6yuoLOyMjA4/EwfPhw/v73v+NyuUhISGDy5Mm8+eabuU6bmJjIL7/8wo033pjr8N27d/PAAw8wffr0fLt9zsvEiRPp27cvDRo0AOD222/noYceonXr1oWeR6DSBKBylZLu5Ykv19AssjL/7Ff8OzOLJOdOKY+dysgLmvLjur08+80f9G4ZSW3tOK7c8+8HaN++fdx4440cPXqUp59+mtjYWGJj826ynpiYyMcff5xrAsjMzKRBgwZndPfsxIkTadu2bXYCGD9+fLHnFWi0CkidbnR13nr6bpIOp/D80UcID3E7HdEp3C7h+avbkZyWyZjZ650ORxVR3bp1GTduHG+//TbGGOLj47niiisAmD9/Ph07dqRjx4506tSJ48eP8+ijj7Jw4UI6duzIa6+9xsSJE7n++uu58sor6du3L4mJibRt2zZ7/jt37qRfv360bNky+47enOO88sorjB49munTp5OQkMDQoUPp2LEjKSkp9OrVi6wbUKdOnUq7du1o27YtjzzySPb0VapU4YknnqBDhw7ExcWxd+/esth0JU4TgDrNZl8D3vdezjWuBcS5yucO9ux6Vbn9gqZMX57E0m3FvCFNOaZp06b4fD727Tu1Rdcrr7zCO++8w8qVK1m4cCERERG88MILXHDBBaxcuZIHH3wQsLqKnjRpEj/99NNp8166dClTpkxh5cqVfPbZZ+TXm8B1111HbGxs9vgRERHZw3bv3s0jjzzCTz/9xMqVK1m2bBkzZlit4JKTk4mLi2PVqlX07NmT999/vyQ2S5nTBKBOYYzhycxbqEQqj4d8bBVm3bWbuMjZ4HJ4oE9zGtaI4MkZa/QGsQCU20X87t2789BDD/Hmm29y5MiRPJ8ElrOr6JzDateuTUREBNdccw2LFhXve7ts2TJ69epFZGQkHo+HoUOHsmDBAgBCQ0Ozz1pKo5vmsqIJQJ3iq5W7+dXXhn96plFHSvahISWtUqiH0QPasHHvCSYsyr/TL1W+bN26FbfbTd26p3by9+ijjzJ+/HhSUlKIi4tj/frcz0BzdhXtT3K0VRaRU7pvBk7pbjkv+bUyCwkJyV5OaXTTXFY0AahsqRleXvxuPe1lC0Pc85wOp1AuaV2PPq3q8vZPmzl4wtlmqqpw9u/fz1133cV999132s56y5YttGvXjkceeYTY2FjWr19P1apVs3vhLIy5c+dy6NCh7EdFdu/enXr16rFv3z4OHjxIWloa33zzV4u0vObfrVs35s+fz4EDB/B6vUydOpULL7yw+CteDmkrIJXtg0Xb2HM0lddDp+CSctDG3r9VUD7NDB/r34pLX1/IWz9tZvSA8vOwk3KtgGabJS3rcZBZzUCHDRvGQw89dNp4r7/+OvPmzcPtdtO6dWsuu+wyXC4XHo+HDh06MGLECGrWrJnvsnr06MGwYcPYvHkzN954Y3YLo1GjRtGtWzeaNGlCq1atsscfMWIEd911FxERESxevDi7vH79+owZM4bevXtjjKF///4MHDiwhLZI+aDdQSsADp5I48KX44lrWpvx2/rkPtKIZOvvxLxPv/M10T7qGnFF0actYIf1+Jdr+HTZTuY+dCFN6hQzviCm3UEHr1LtDlpEokVknoisE5HfReRvdvloEdklIivtV3+/aR4Tkc0iskFELvUr72eXbRaRR4u0lqpUvfnjJlIyvDx6WauCRy6H/n5xC8I8Ll76rny2WlKqPCrMNYBM4P+MMecAccC9IpJ1i9xrxpiO9msWgD1sMNAG6Af8V0TcIuIG3gEuA1oDQ/zmoxy0df8JpizZwZCu0TSvW8XpcIqlbtVw7rywGbPX/klCojYLVaowCkwAxpg9xpgV9v/HgXVAw3wmGQhMM8akGWO2AZuBrvZrszFmqzEmHZhmj6sc9p+5GwnzneRvv/Uv+S4CytDtFzShbtUwXvpug/YTpFQhFKkVkIg0BjoBS+yi+0RktYhMEJGsKzMNgZ1+kyXZZXmV51zGHSKSICIJ+/frA8FLjd22f/2otsxavYtb3N8RWc6bfRakUqiH+y5qztLEQ/yypWSfU6xUMCp0AhCRKsDnwN+NMceAsUAzoCOwB/hP1qi5TG7yKT+1wJhxxphYY0xsZGRkYcNThZHLoxjfyLyGKqRyu2eWg4GVnEFdoqlfPZxX527UswClClCoBCAiIVg7/ynGmC8AjDF7jTFeY4wPeB+rigesI/tov8mjgN35lCuH/OGLYbavG7e4Z1NDkp0Op0SEedzc27s5y7cfZuGmA06Ho1S5VphWQAJ8AKwzxrzqV17fb7SrgbX2/zOBwSISJiJNgBbAUmAZ0EJEmohIKNaF4pklsxqqOF7PvJaqJHObZ7bToZSoG2KjaVgjQs8CyhERYdiwYdnvMzMziYyMzO5OobAaN27MgQOa2EtKYc4AugPDgItyNPl8SUTWiMhqoDfwIIAx5nfgU+AP4DvgXvtMIRO4D5iDdSH5U3tc5YC1vkZ87+vC7Z5ZVJeTTodTokI9Lu67qDkrdx4hfqNeRyoPKleuzNq1a0lJSQGsu3UbNsyvLYkqC4VpBbTIGCPGmPb+TT6NMcOMMe3s8gHGmD1+0zxvjGlmjGlpjJntVz7LGHO2Pez50lopVbC3M6+mGsnc4v7O6VBKxXWdo4iqGcGbP27Ss4By4rLLLuPbb78FrG6WhwwZkj3s0KFDXHXVVbRv3564uDhWr14NwMGDB+nbty+dOnXizjvvPOWz/Oijj+jatSsdO3bkzjvvxOv1lu0KBQHtC6gC2uxrwBxfLMPd31NNUpwOp3ByuYCdnxC3izt7NuW3HUdYot1F56JXCb8KNnjwYKZNm0ZqaiqrV6+mW7du2cOeeuopOnXqxOrVq/n3v//NzTffDMDTTz9Njx49+O233xgwYAA7duwArLtfP/nkE37++WdWrlyJ2+1mypQpRVh/BdoXUIX0rvdKwshghCc4j/6zXB8bzRs/bmJs/BbimtZ2OpwKr3379iQmJjJ16lT69+9/yrBFixbx+eefA3DRRRdx8OBBjh49yoIFC/jiiy8AuPzyy7P7Afrxxx9Zvnw5Xbp0Aay+hnL2LKoKpgmggtllajPD252b3D9QWwrfw2IgCg9xc0v3Jrw8ZwNrdx2lbcPAvcmt5MU7stQBAwbw8MMPEx8fz8GDf92rkVs1XVZPoTl7DM0af/jw4YwZM6b0gq0AtAqoghmfaR15jfR863AkZWPYeY2oGuZh7PwtToeigFtvvZVRo0bRrl1kHlhRAAAgAElEQVS7U8p79uyZXYUTHx9PnTp1qFat2inls2fP5vDhwwD06dOH6dOnZz9R7NChQ2zfvr0M1yQ46BlAsPOrMz9kqjLN25uBrp9pKBXjTtlq4SEMjWvEuAVb2HYgWXsKdVhUVBR/+9vfTisfPXo0t9xyC+3bt6dSpUpMmjQJsK4NDBkyhHPPPZcLL7yQmJgYAFq3bs1zzz1H37598fl8hISE8M4779CoUaMyXZ9ApwmgApmU2ZcUwrnb87XToZSpW3s0ZsLP2xi3YCtjrmlX8ASqxJ04ceK0sl69etGrVy8AatWqxVdffXXaOLVr1+b777/Pfv/aa69l/z9o0CAGDRpU8sFWIFoFVEGkmhA+8l7Mxa7lNHdVrBuw61YN59pzG/LFiiR9aphSfjQBVBAzvN05SHVucwdHnz9FdWv3JqRl+piyZIfToShVbmgCqACMgfHe/rSRbcS51jkdjiNa1KtKr5aRTF6cSGqG3jCkFOg1gOCTy41S833t2WyieC3kHXJpURd4cq5jIZ9ve3uPptz0wRJmrtrNDbHRBU+gVJDTM4AK4ANvf+pxiMtdvzodiqO6N69Nq7OqMmHRNu0eQik0AQS99b5oFvrac7Pne0KlYld9iAi39WjC+j+P8/PmitEMVqn8aAIIcv/z9iOcNIa6f3Q6lHJhQMcG1KkSxgeLtjodSoVSUt1B9+rVi4SEBAD69+/PkSNHSjTOikYTQBA7bKoww9udq92LguaBL2cqzONmaLcY4jfuJ/GAbpOyUhrdQc+aNYsaNWqURHgVliaAIPaJtxdphDLc/X3BI1cgQ7vF4BZh8mLtOqAs5dcddHJyMrfeeitdunShU6dO2TeFpaSkMHjwYNq3b8+gQYOyEwic+nCYq666is6dO9OmTRvGjRuXPU6VKlV44okn6NChA3Fxcezdu7csVjVgaAIIUl4jfJh5Cd3kD1q5djodTrlSt1o4/dvV57PlO0lOy3Q6HGf0KuFXIeTXHfTzzz/PRRddxLJly5g3bx7/+Mc/SE5OZuzYsVSqVInVq1fzxBNPsHz58lznPWHCBJYvX05CQgJvvvlmdkdzycnJxMXFsWrVKnr27Mn7779fuGArCG0GGqR+8HVmF5E86fnI6VBKX17PCMineejw8xsxc9VuvvxtFzfFaf8xZSG/7qC///57Zs6cySuvvAJAamoqO3bsYMGCBTzwwAPZ07dv3z7Xeb/55pt8+eWXAOzcuZNNmzZRu3ZtQkNDs68zdO7cmblz55bW6gUkTQBBapK3Lw04wCWu3I+YKrpzY2rStmE1Ji9OZGi3mFy7HA5q8c4sNr/uoD///HNatmx52jQFfTbx8fH88MMPLF68mEqVKtGrVy9SU1MBCAkJyZ7e7XaTmVlBz/jyoFVAQWijryG/+Noy1PMDHvE5HU65JCIMP68xG/eeYPEWbRJaVvLqDvrSSy/lrbfeyr4/47fffgNO7SZ67dq12Y+K9Hf06FFq1qxJpUqVWL9+Pb/+WrHvdykKTQBB6EPvJYSSzhD3PKdDKdeu7NCAWpVDmbQ40elQKoy8uoP+17/+RUZGBu3bt6dt27b861//AuDuu+/mxIkTtG/fnpdeeomuXbueNm2/fv3IzMykffv2/Otf/yIuLq7U1yNYaBVQMPCrAz9hwvnCewFXuJZQK8if+HWmwkPcXB8bxfiF2/jzaCpnVQ93OqSgVVB30BEREbz33nunjRMREcG0adNynWdiYmL2/7Nnzy5wuddddx3XXXddEaIOfnoGEGRmeLuTTAQ3efRiV2EM7doInzFMXaq9hKqKRxNAEDEGPvJeTGtJpJNsdjqcgBBTuxIXnh3JtGU7yPDq9RJVsWgCCCLLzdmsN40Y5p4bHL1+lpFhcY3YeyyNH/4I7puEtAO84HOmn6kmgCDyUebFVOUkA92/OB1KQOnVsi4Na0Tw0ZLgvTM4PDycgwcPahIIIsYYDh48SHh48a9d6UXgIHHQVGWWrxtD3D9RSfSxh0Xhdgk3dovh5Tkb2LL/BM0iqzgdUomLiooiKSmJ/fv3Ox2KKkHh4eFERUUVe3pNAEHiM++FpBPCTe4fnA6l/CjCg2NuiI3m9R828vGSHfzritalHFjZCwkJoUmTJk6HocqZAquARCRaROaJyDoR+V1E/maX1xKRuSKyyf5b0y4XEXlTRDaLyGoROddvXsPt8TeJyPDSW62KxWeEqd6L6CrraOHa5XQ4ASmyahiXtjmLz1ck6SMjVYVRmGsAmcD/GWPOAeKAe0WkNfAo8KMxpgXwo/0e4DKghf26AxgLVsIAngK6AV2Bp7KShjozv/jasN2cxY2en5wOJaDd2DWGIyczmL12j9OhKFUmCkwAxpg9xpgV9v/HgXVAQ2AgMMkebRJwlf3/QGCysfwK1BCR+sClwFxjzCFjzGFgLtCvRNemIhldPfs11XsRNTlOP9dSp6MKaOc1q02TOpX5eIneE6AqhiK1AhKRxkAnYAlQzxizB6wkAdS1R2sI+Pc/nGSX5VWecxl3iEiCiCToBauC7TfVmOOL5Vr3AsIlw+lwApqIMKRrNMsSD7Nxr95FrYJfoROAiFQBPgf+bow5lt+ouZSZfMpPLTBmnDEm1hgTGxkZWdjwKqzPvBeSiYchbq3+KZDfWVNerj03ilC3S88CVIVQqAQgIiFYO/8pxpgv7OK9dtUO9t99dnkSEO03eRSwO59yVUw+I0zzXkSc63eaubTeuiTUrhLGpW3P4gu9GKwqgMK0AhLgA2CdMeZVv0EzgayWPMOBr/zKb7ZbA8UBR+0qojlAXxGpaV/87WuXqWL62deGHaaeHv2XsBu7xnAsNZNvV2tSVcGtMGcA3YFhwEUistJ+9QdeAC4RkU3AJfZ7gFnAVmAz8D5wD4Ax5hDwLLDMfj1jl6limurtY1/8XeZ0KEElrmktmtaprB3EqaBX4I1gxphF5F5/D9Anl/ENcG8e85oATChKgCp3B0w15vo6M9w9hzDRpxyVJBFhUJdoxsxez6a9x2lRr6rTISlVKrQvoEDhfwFzdHU+9/YkAw+D9aEvpeLazlGEuIVpy3YWPLJSAUoTQAAyBqZ5e9NF1tPcpdfRS0OdKmFc0rqeXgxWQU0TQAD61XcO20x9Bnv06L80De4Sw+GTGcz5/U+nQ1GqVGgCCEDTvBdRjWQud+nDr0tTj+Z1iKoZwbSlWg2kgpMmgABzxFRmtq8LV7sX6Z2/Z6IQN4W5XMLgLtEs3nqQxAPJZRicUmVDE0CA+cJ7AemEMljb/peJ62Ojcbv0YrAKTpoAAogx8Im3Fx1kM+e4dIdUFupVC6d3y7pMX56kzwxWQUcTQABZaZqxwcRo088yNrhLNAdOpPHT+n0Fj6xUANEEEECmeS+iEqlc6V7sdCgVSq+WkdSrFsY0vTNYBRl9JGSAOGHC+dp7Hle6F1NFUp0OJ7gU8OhIj9vF9Z2j+W/8ZnYfSaFBjYgyDE6p0qNnAAHia+95nCScQVr944gbYqPxGfgsIcnpUJQqMZoAAsQ0b29ayg46yWanQ6mQYmpXokfzOnyasBOv77THWCgVkDQBBIB1e46xyjRnkDseyatbPlXqBnWJZteRFBZtPuB0KEqVCE0AAeCTZTsJJYOr3YucDqVC69umHjUrhfDJMr0YrIKDJoByLjXDyxcrkujnWkpNOeF0OBVamMfNNedGMfePvRw4keZ0OEqdMU0A5dno6sx55gqOpWZq2/9yYnCXaDK8hi9X7HI6FKXOmCaAcm6qtzcxspc41zqnQ1FAi3pV6dyoJlOX7cB69pFSgUsTQDm2zXcWv/raMMg9D5fozqa8GNQlmq37k0nYftjpUJQ6I5oAyrFPvL1w4+U69wKnQ1F+rmhfnyphHn1msAp4mgDKqQyvj+nenvR2/UY9OeJ0OBVLAV1FVwr1MKBjA2at2cPRFO2SWwUuTQDl1I/r9nGAGgzRi7/l0pAuMaRm+Ji5Ui8Gq8ClCaCcmrZsB/U4xIWuVU6HonLRtmE1WtevxtSlO/VisApYmgDKoV1HUpi/cT83uOfjEe2DvjwSEYZ0jeaPPcdYs+towRMoVQ5pAiiHPrWfPnWDVv+UawM7NSQ8xMVUfWawClCaAMoZ71M1+OzHxVwgq4h2aZ8z5Vm18BAub9eAmSt3kZyW6XQ4ShWZJoByZoGvA7upwxB95m9AGNI1muR0L9+s3u10KEoVmT4QppyZ6u1NHY7Sx7XC6VAUFPiwmM6NatK8bhWmLt3JoC4xZRiYUmeuwDMAEZkgIvtEZK1f2WgR2SUiK+1Xf79hj4nIZhHZICKX+pX3s8s2i8ijJb8qgW/fsVR+9J3Lte75hIrX6XBUIYgIg7tEs3LnEdb/eczpcJQqksJUAU0E+uVS/poxpqP9mgUgIq2BwUAbe5r/iohbRNzAO8BlQGtgiD2u8vPZ8iS8uBnsjnc6FFUE154bRajbxdQlemewCiwFJgBjzALgUCHnNxCYZoxJM8ZsAzYDXe3XZmPMVmNMOjDNHlfZfD7D1KU7ON+1liauP50ORxVBzcqhXNbuLL74bRcp6XrmpgLHmVwEvk9EVttVRDXtsoaAf5u4JLssr/LTiMgdIpIgIgn79+8/g/ACy8LNB0g6nKIXfwPUjV1jOJ6aqReDVUApbgIYCzQDOgJ7gP/Y5bk9sNDkU356oTHjjDGxxpjYyMjIYoYXeKYu2UHtyqFc6lrmdCiqGLo2qUWzyMraQZwKKMVKAMaYvcYYrzHGB7yPVcUD1pF9tN+oUcDufMoV1sXfuev2cl3nKL34G6CsO4NjWLFDLwarwFGsBCAi9f3eXg1ktRCaCQwWkTARaQK0AJYCy4AWItJEREKxLhTPLH7YweWz5Ul4fYbBXbUZYSC79twoQj16MVgFjgLvAxCRqUAvoI6IJAFPAb1EpCNWNU4icCeAMeZ3EfkU+APIBO41xnjt+dwHzAHcwARjzO8lvjYByPdUDaamv8b5rn00eftGp8NRBfG/LyDHPQE1K4fSv611MfjRy84hItRdxsEpVTQFJgBjzJBcij/IZ/zngedzKZ8FzCpSdBXAAl87kkxdHvFMczoUVQKGdI1hxsrdfL16NzfERhc8gVIO0q4gHDbFezG1OUpfV4LToagS0LVJLZrXrcIUrQZSAUATgIN2H0nhR9+53OCOJ0y0M7GAk8uTw0SEod1iWLXzCGu1m2hVzmkCcNC0ZTsxwI3uH50ORZWga86NIiLEzZQl250ORal8aQJwSIbXx7SlO+jl0m6fg031iBAGdGjAjN92cyxVnxmsyi9NAA754Y+97DuexlA9+g9KQ+NiSMnw8uUKfWawKr80AThkypIdNKwRQW/Xb06HokpB+6gatI+qzpQl2/WZwarc0gTggG0Hklm0+QBDukbjFt05BKubujVi494TLEs87HQoSuVKE4ADPvp1Ox6XcEMXbScezK7s0IBq4R4+/FUvBqvySRNAGTuZnsmnCTu5rF196lYNdzocVYoiQt1cHxvN7DV72Hcs1elwlDqNJoAy9tXK3RxPzeTmdXee/rhBFXRuimtEps8wdenOgkdWqoxpAihDxhgm/ZLIOZJIrGxwOhxVBprUqcyFZ0cyZcl2Mrw+p8NR6hSaAMpQwvbDrP/zODe75yK5PSFBBaWbz2vEvuNpfP/7XqdDUeoUmgDK0OTF26kW7mGg+xenQ1FlqFfLukTXimDy4kSnQ1HqFAX2BqpKxr5jqcxes4fh5zemUkKa0+Gokpbzeo5fV9Ful3BTt0aMmb2e9X8eo9VZ1co4OKVyp2cAZWTKkh14jeGmuEZOh6IccENsNGEeF5N+0SahqvzQBFAG0jK9TFmynd4t69KkTmWnw1EOqFk5lKs7NeTL35I4cjLd6XCUAjQBlIlvVu3hwIl0btnyd236WVHk0lX0iO6NSc3waZNQVW5oAihlxhj+98s2WkgSPVxrC55ABa1WZ1Xj/Ga1+XBxIpnaJFSVA5oASlnC9sOs3XWMEe7vtOmn4pbuTdh9NJU52iRUlQOaAErZ/37eRvWIEK5xL3I6FFUOXNSqLjG1KvG/n7c5HYpSmgBK064jKcz5fS+Du0YTIXrhT1lNQoef35iE7YdZnXTE6XBUBacJoBRNtI/ybj6vsbOBqHLlhtgoqoR5+GCRngUoZ2kCKCXHUjOYunQnl7erT8MaEU6Ho5yUo0VQ1fAQBneJ5pvVe9h9JMXh4FRFpgmglHyydCcn0jIZeUFTp0NR5dAtPZoA6LUA5ShNAKUgw+vjfz9vI65pLdpFabt/dbqGNSLo364+05bu5Lg+OF45RBNAKZi1Zg+7j6bq0b/K18gLmnA8LZNPlumNYcoZmgBKmDGG9xdupWlkZXq3rOt0OKocax9Vg65NajFh0TZ9VoByRIEJQEQmiMg+EVnrV1ZLROaKyCb7b027XETkTRHZLCKrReRcv2mG2+NvEpHhpbM6zlu85SBrdx1j5OHXcT1T47TuAJTyd8cFTdl9NJVZa/Y4HYqqgApzBjAR6Jej7FHgR2NMC+BH+z3AZUAL+3UHMBashAE8BXQDugJPZSWNYDN2/hbqVAnjar3xSxXCRa3q0rxuFcbGb8EY43Q4qoIpMAEYYxYAh3IUDwQm2f9PAq7yK59sLL8CNUSkPnApMNcYc8gYcxiYy+lJJeCtSTrKwk0HuK1HE8JFL+ypgrlcwl0XNmP9n8eJ37Df6XBUBVPcawD1jDF7AOy/WZXdDQH/K1pJdlle5acRkTtEJEFEEvbvD6wfxLvzt1A1zMPQuBinQ1Hllf89AXbV4MCODWhQPZz/xm92ODhV0ZT0ReDcujsz+ZSfXmjMOGNMrDEmNjIyskSDK03bDiQza+0ehp3XiGrhIU6HowJIiNvFyJ5NWZZ4mGWJOU+2lSo9xU0Ae+2qHey/++zyJCDab7woYHc+5UFj3IIthLhd3NK9idOhqAA0uEsMtSqH8m78FqdDURVIcRPATCCrJc9w4Cu/8pvt1kBxwFG7imgO0FdEatoXf/vaZUFh77FUPl++ixtio4isGuZ0OCoARYS6GXF+Y35cv491e445HY6qIArTDHQqsBhoKSJJInIb8AJwiYhsAi6x3wPMArYCm4H3gXsAjDGHgGeBZfbrGbssKLw7fwteY7izZzOnQ1EBbPh5jakS5uHteXotQJUNT0EjGGOG5DGoTy7jGuDePOYzAZhQpOgCwL5jqXy8ZAfXdGpIdK1KToejAlj1SiEMP78R/43fwqa9x2lRr6rTIakgp3cCn6FxC7aS6TPc27u506GoQJSjRdBtPZoSEeLmrZ/0LECVvgLPAFTe9h9P46Ml2xnYsQGN61TWO37VGatVOZSbz2vMewu28ECfFjSvW8XpkFQQ0zOAM/D+wq2kZ/q4T4/+VUmwzwRG/nox4R43b/+0yemIVJDTBFBMB0+k8eHi7Qzo0ICmkXqUpkpObTnOMN9XzFyZxNZRLfXMUpUaTQDFNDZ+C2mZXu67qIXToaggNNLzDWGk83rmtU6HooKYJoBi2HM0hcm/bueac6O0jlaVikg5xgj3HGb6urPOF13wBEoVgyaAYnjrp80YY/hbHz36V6XnTs83VCWZ/2Re73QoKkhpAiii7QeT+XTZToZ0jdF2/6pU1ZBk7vR8ww++WFbsOOx0OCoIaQIootd/2ITHLdryR5WJW9zfUZuj/Of7DU6HooKQ3gdQBOv/PMaMlbu4s2cz6lYLtwq1hYYqRZUljXs8X/Hs5ptZtOkAPVrUcTokFUT0DKAIxsxaT9UwD3ddqA97V2VnqPtHGtaI4N+z1uHz6VPDVMnRBFBICzbuZ/7G/TzQpwU1KoU6HY6qQMIlg3/2a8kfe47x5W+7nA5HBRFNAIXg9Rn+PWsd0bUiGHZeI6fDURXQle0b0D6qOq98v4HUDK/T4aggoQmgED5fkcT6P4/zz0tbEeZxOx2OqoBcz9Tg8X3/x56jqXywaJvT4aggoQmgACfTM/nP9xvoGF2DK9rXdzocVYHFudZziSuBsXN+48BT0doAQZ0xTQAFGBu/hb3H0nji8nMQye3RxkqVnUc9U0kllFcyb3A6FBUENAHkY/vBZN6bv5WBHRvQpXEtp8NRimauPdzi/o5PvL1Y5dPWaOrMaALIx7Pf/IHHLTx22TlOh6JUtgc8X1KHo4zKGKHNQtUZ0QSQh3kb9vHDun3cf1ELzqoe7nQ4SmWrKik8GjKNVaY505cnOR2OCmCaAHKRlunlma//oGmdytzao7HT4Sh1mqtdi+gsG3jxu/UcTclwOhwVoDQB5OLd+K1sO5DMqCtba7NPVS65xPB0yEQOn0znpe/WOx2OClDaF1AOm/cd5515m7myQwN6tax7+gja9E6VE21d27nFNYsPlvTnqt9uo4srR4dxo486E5gKGHoG4MfnMzz+xVoiQt2MuqK10+EoVaCHPJ/RkP08lnEbaUaP51TRaALw80nCTpYmHuKJ/ucQWTXM6XCUKlBlSeO5kP+x2UTxrvdKp8NRAUYTgG3vsVT+PWsdcU1rcX1slNPhKFVovd0rudL1C+9kXsUmX0Onw1EBRBMAYIzhkc9Xk+H18e+r2+kdvyrgjAqZTGVS+b+Mu8gw2nBBFY4mAOCTZTuJ37CfR/u1omlkLg95H139r5dS5VCkHOP5kA9YbZox1jvA6XBUgDijBCAiiSKyRkRWikiCXVZLROaKyCb7b027XETkTRHZLCKrReTckliBM7Xz0Eme/eYPzmtam5vPa+x0OEoVW3/3Uga4fubNzKtZ62vsdDgqAJTEGUBvY0xHY0ys/f5R4EdjTAvgR/s9wGVAC/t1BzC2BJZ9Rnw+w8OfrUJEePn69rhcWvWjAtszIROpxXEeyrib1Kfq6JmryldpVAENBCbZ/08CrvIrn2wsvwI1RMTR/pXHzt/Ckm2HGHVFa6JqVnIyFKVKRA1J5sWQcWw00YzJvNHpcFQ5d6YNhw3wvYgY4D1jzDignjFmD4AxZo+IZN1N1RDY6Tdtkl225wxjKJaExEO8OncjV7Svn3urHz1qUgGqt3sVt/pmMcHbn/Ncv9PPneB0SKqcOtME0N0Ys9veyc8VkfzuSc+tfuW0rgxF5A6sKiJiYmLOMLzcHTmZzgNTf6NhjQjGXKOtflTwecQzjWW+Vvwz407auhI57RDH/wBH7xiusM6oCsgYs9v+uw/4EugK7M2q2rH/7rNHTwKi/SaPAnbnMs9xxphYY0xsZGTkmYSXV8z8Y/pq9p9I4+0bO1E1PKTEl6GU08Ikk7dD3sSH8ED6fWR4fU6HpMqhYicAEaksIlWz/gf6AmuBmcBwe7ThwFf2/zOBm+3WQHHA0ayqorI0dv4W5v6xl0cvO4f2UTXKevFKlZlGrn2MCRnPCnM2z4/6mzZnVqc5kyqgesCXdvWJB/jYGPOdiCwDPhWR24AdwPX2+LOA/sBm4CRwyxksu1jiN+zj5TkbuLJDA27t3risF69UmbvS/Ssrfc35wNuftq5tXOde6HRIqhwpdgIwxmwFOuRSfhDok0u5Ae4t7vLOVOKBZB6Y+hutzqrGS9e213p/VWE85vmYdSaGxzNuo4XsooNrq9MhqXKiQtwJfCw1gzs+TMDtEsYN60xEqN4qryoOj/h4O+QtIjnKXekPss9o1aeyBH0CSM/0cc9HK9i6P5l3bjyX6Fra3l9VPLXkOONCX+Uolbkl/R8kG7/ebvXaQIUV1AnAGMOjX6xm0eYDjLmmHec3r+N0SEo5po1rO++EvMl6E8O9GX8j0wT1z18VQlB/A16bu5EvVuziwYvP5vrY6IInUCrI9Xav5DnPBOJ9HXky81bMaXfiqIokaB8h9MGibbz502ZuiI3igT7NnQ5HqXJjiGceu0wd3vZeTXVO8KhnGtomomIKygQwZcl2nv3mDy5re1bh+/fX+k9Vgfyf5zMOU5X3vAOoJGn8zfOlNSDn70DvEg5qQZcApi9P4okv13JRq7q8MbgTHndQ13IpVSwi8Kznf6QRwmuZ1xNOBnd6vsl/Iu0+IugEVQL4ZNkOHv1iDT2a1+G/Q88l1KM7f6Xy4hLDi55xpJpQxmTeSCYu7vXMPHUkPTMOakGTAMYv3Mpz367jwrMjefemzoSHFKKtv365VQXnFsNrIf/Fk+Hl5czBHDeVeESvCVQYAZ8AjDG8/sMm3vhxE/3bncXrgzrpkb9SRRAiXl4NGUvlzFTe9Q7gOJV4xvM/3JJPEyG9VhAUAjoBpGf6eOyLNXy+IonrO0cx5pp2Bdf561G/UqdxieE5zwSqcpJ3vQPYa2ryRsjbVJY0p0NTpShgD5WPnExn2AdL+HxFEn+/uAUvXddeL/gqdQZE4NGQaTzj+R8/+TpxffpT7DG1nA5LlaKAPANY/+cx7v5oBbsOp/D6oI5c1amh0yEpFTRu9swlRvZyX8YDDEx7lv+GvkGsa2P+E2kLoYAUcIfM05cncdU7P3MiLZOPR3bTnb9SpaCXezXTQ58mQtIYnP4k4zP7613DQShgzgCS0zJ55us/+CRhJ+c1rc0bQzpSt2q402EpFbRauXbydegT/CPjTp7LvIllvrN5IWQ8NeVE/hPqBeKAERBnAMsSD3HZGwv5dPlO7uvdnI9u76Y7f6XKQDVJ4d2Q13nS8yE/+c6lb9qL/OTt6HRYqoSU6zMAnzE8/+0fjF+0jaiaEXxyx3l0bVLEi1La6kepMyICt3tmc77rdx7KuIdbM/7JDb55POH5mOqS7HR46gyU6wSwce8J3l+4jSFdY3ji8nOoElauw1UqqLV27eCr0Cd5PfNa3vNeyY/ec3ks5GOudS3M/8axvA7CtGrIceW6Csgtwmd3nceYa9rpzl+pciBMMnkk5BO+Dn2CRrKXhzPu5ob0Uaz0NXM6NFUM5Xqv2rxeFbo01nbISpU3bVzbmR76NJ96L+TlzEFclf4s/V1LeNjzCU1dfxZuJnqx2DKcKtQAAAfjSURBVHHlOgEUuzsSrfdXqtS5xDDYE88V7l95P/Ny3vdeznfpXRjg+oW7PTNp6Uoq2gz1XoIyV64TgFKq/KsiqTwY8jk3eebyfuYVfOS9mBnpPbjYlcAt7jmc7/q96J3L5XcQp8mhxARHAtAjfqUcFynHeDzkY+72zOR/mZfyobcvP/hiaSa7uMn9AwPdv1BLjjsdpvITHAlAKVVu1JQTPBTyOfd4ZvKNL44PM/vydOZw/p05lN6u37javYherlVESHrxFqBVRSUmcBOAHvUrVa6FSwbXuRdynXshf/hi+NLbgxne7nzv60I4afRyreJS9zIucK2hjhwr3kIKW1WkSSNXgZsAlFIBo7VrB61dH/OIZxpLfa2Y7evKHG8XvvN1BaCtbKOHaw3dXOvo7NpINUk584XqQWKBAisB6AeqVEDziI/z3X9wvvsPnvZMYq1pzAJfexZ42zPe2593vQNw4aOV7KCDayvtZQvtXFtpLrsJl4ySCUKbn2YLrASglAoaLjG0l220d23jPs9XnDRh/OZrzlJfK1aYFnzr7cZULgLAjZcmsoeWkkRT2UNT126ayW7au7adeSDFObAMkqRR5glARPoBbwBuYLwx5oU8R979mx71K1VBVJI0urt/p7v7dwCMge2mHmtNYzb4YlhnollrGjPb1xWf10V9DrI4/H6Ho6bw+6hymDTKNAGIiBt4B7gESAKWichMY8wfZRmHUqr8E4HGspfG7OUK95Ls8jTjYYepxxGqOBdcsc4ainkwW4qJo6zPALoCm40xWwFEZBowENAEoJQqlDDJpIXscjqMslOKtSBiyvAxPyJyHdDPGHO7/X4Y0M0Yc5/fOHcAd9hvWwIbSmDRdYADJTAfJwXDOkBwrIeuQ/mg65C3RsaYyIJGKuszgNxuCD8lAxljxgHjSnShIgnGmNiSnGdZC4Z1gOBYD12H8kHX4cyVdXfQSUC03/soYHcZx6CUUoqyTwDLgBYi0kREQoHBwMwyjkEppRRlXAVkjMkUkfuAOVjNQCcYY34vg0WXaJWSQ4JhHSA41kPXoXzQdThDZXoRWCmlVPlRrh8JqZRSqvRoAlBKqQoqKBKAiDwoIr+LyNr/b+9cQuOqwjj++2MxbQRtkhK1JZJGumhxFboILrIwkj6QxPqAQMH4WOlCXYhGBgQRF1VBcJWNgorG+sS60DbW5yYVlKad0tZMTW0ttRWq3RRqg5+L801yM0xSppC592bODy733O+cge8/c+Z+9zyvpDFJK32g+YCkKUm7fdAZSU1+XfL8znS9n0PSU67hiKSn3dYqadx1jEtqcbskveE6DknqTsnntySdl1RM2Gr2WdKwl5+SNJwBDQ/67/CfpM0V5Z93DcclbUnYt7qtJGkkAxpelXTMv+vPJK3OsoZFdLzkGg5K2idprdtzU58Sec9IMklrMqHBzHJ9AOuAaWCVX38IPOznIbeNAo97+glg1NNDwO60NbgvdwBFoJkwOP81sAF4BRjxMiPALk9vB74krK3oAQ6k5Hcv0A0UE7aafAZagd/83OLplpQ1bCQsRPwO2JywbwImgSZgPXCCMKHhOk93Add7mU0pa+gHVnh6V+J3yKSGRXTcmEg/mfj/5qY+ub2DMAHmd2BNFjQsixYA4Ya5StIKwg30LHAX8LHnvw3c6+lBv8bz+6Sa31i6FGwEJszskpnNAN8DO5jvb6WOdywwAayWdGu9nTazH4ALFeZafd4CjJvZBTP7GxgHti6994FqGszsqJlVW4U+CHxgZpfNbBooEbY4md3mxMz+BcrbnNSFBTTs87oEMEFYdwMZ1eA+V9ORfFvMDcwtHs1NfXJeB55l/uLXVDXkPgCY2RngNeAU4cZ/EfgZ+CdR+f8gtBTw82n/7IyXb6unzwtQBHoltUlqJjwZdAA3m9lZAD+3e/lZHU5SY9rU6nOWtVSSVw2PEp40IYcaJL0s6TSwE3jBzbnRIWkAOGNmkxVZqWrIfQDw/uVBQlN2LeEJYVuVouWoe9XtKNLAzI4SmunjwFeE5vfMIh/JpI6rsJDPedKSOw2SCoS69F7ZVKVYpjWYWcHMOggaynuH5UKHP9AVmAtc87Kr2OqmIfcBALgbmDazv8zsCvApcCehKVVe6JbccmJ2OwrPv4nqzbW6Y2Zvmlm3mfUSfJoCzpW7dvx83otneVuNWn3OspZKcqXBBw/vAXaady6TMw0VvA/c7+m86Lid8IA6Kemk+/OLpFtIWcNyCACngB5Jzd6X30fYXvpb4AEvMwx87uk9fo3nf5P4Y6SKpHY/3wbcB4wx399KHQ/5LIIe4GK52yUD1OrzXqBfUou36PrdlkX2AEMKs8nWEwbqfyKD25wovHzpOWDAzC4lsnKjAUDShsTlAHDM07moT2Z22MzazazTzDoJN/duM/szdQ31GhlfygN4kVApisC7hNkNXYRKXQI+Apq87Eq/Lnl+V9r+J3T8SAhek0Cf29qA/YTWwH6g1e0ivFznBHCYxEyVOvs8Rhh7ueIV+7Fr8ZnQR13y45EMaNjh6cvAOWBvonzBNRwHtiXs24FfPa+QAQ0lQj/yQT9Gs6xhER2f+H/7EPAFsC5v9aki/yRzs4BS1RC3gohEIpEGZTl0AUUikUjkGogBIBKJRBqUGAAikUikQYkBIBKJRBqUGAAikUikQYkBIBKJRBqUGAAikUikQfkf+0Z2zSjnvPEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "relative_error_plx = 10. # in percents\n", "plotdis(relative_error_plx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the expert: \n", "\n", "Given the above issue, it is difficult to estimate the true distance of an object observed by Gaia if the relative error is larger than 20%. One can estimate some intervals, but the range will be large.\n", "To remedy to this, one can use some priors - for example the fact that the Milky Way has a given distribution of stars and thus a star is more likely to be somewhere (e.g., in the spiral arms) than somewhere else. This is the so-called geometric distance. \n", "One can also use the fact that a star occupy a given position in the colour-magnitude diagram (or Hertzsprung-Russell diagram) and thus needs to be at a given distance. This is the photo-geometric distance. It is, however, sometimes dangerous as without spectroscopy, we do not know if a star is a main sequence star or a red giant for example, and this could affect significantly its distance.\n", "A reprocessing of all Gaia EDR3 parallax to use such priors has been done by Coryn Bailer-Jones and colleagues and is explained at https://www2.mpia-hd.mpg.de/homes/calj/gedr3_distances/main.html. The corresponding distances are available via the Gaia science archive or using, for example, the incredible tool TopCat." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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", "version": "3.7.4" } }, "nbformat": 4, "nbformat_minor": 2 }