{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "The colour-magnitude diagramme\n", "\n", "One important tool in astrophysics is the Hertzsprung–Russell (HR) diagramme, in which we plot the luminosity of a star as a function of its surface (or \"effective\") temperature. In practice, we use a logarithm-logarithm plot, as the ranges covered are rather large (10 orders of magnitudes in luminosity, 2 in temperature). Such diagrammes have proved very useful, as it was seen that the majority of stars (including the Sun) follow a line called the \"main sequence\", while others will occupy different regions in the diagramme, corresponding to giant stars, horizontal-branch stars or white dwarfs, for example.\n", "\n", "In principle, we do not know necessarily the temperature (which is obtained mostly with spectroscopy) nor the total luminosity (as we would need to observe it at all wavelength) of a star and so a different diagramme, the colour-magnitude diagramme (CMD), is used as a proxy. \n", "\n", "In astronomy, for historical reasons, we use magnitude in a different way as in other part of physics. Indeed, the Greek astronomer Hipparchus of Nicaea devised a numerical scale to describe the brightness of stars as seen with the unaided eye: the brightest star in the sky was assigned an apparent magnitude of 1, and the dimmest star visible assigned a value of 6. The difference between them corresponds roughly to a factor of 100 in brightness. Thus, as 5 magnitudes correspond to a factor 100, the magnitude is given as -2.5 times the logarithm of the brightness in a given wavelength range. Because of the minus sign, a star is brighter the smaller (or more negative) its magnitude is, as devised by Hipparchus.\n", "\n", "We also use two kinds of magnitude - the apparent (or observed) magnitude as described above and the absolute one. The apparent one tells us how we perceive the star as seen from Earth. However, a star can be close or far, and thus its intrinsic luminosity may be very different. This is taken into account into the absolute magnitude, which takes the distance into account. The later is thus a good proxy for the total luminosity of the star. Thus, the brightest star in the night sky, Sirius, has an observed magnitude of -1.46 because it is so close to us, but an absolute magnitude of 1.4. By comparison, one of the intrinsically brightest stars in our Galaxy, Eta Carinae, has an absolute magnitude of -8.6 (it is thus intrinsically 10,000 times brighter than Sirius), but an apparent magnitude of 4.3 (so it seems to us , as it is much farther away.\n", "\n", "The apparent magnitude of a given band, for example in the visible, is given by \n", "V = -2.5 log FV + CV,\n", "where CV is a constant and FV is the observed flux of the star in this band.\n", "To determine the absolute magnitude, we need to take into account the distance, as well as the fact that the flux of the star is dimmed by all the matter and dust that lies between the star and us (the \"extinction\"). By definition, an object's absolute magnitude is defined to be equal to the apparent magnitude that the object would have if it were viewed from a distance of exactly 10 parsecs (32.6 light-years), without extinction. \n", "MV = -2.5 log (FV * d102) + CV + AV, \n", "where d10 is the distance of the star, in units of 10 parsecs, d10= d / 10 pc, with d expressed in parsec, and AV is the extinction.\n", "\n", "Thus,\n", "MV = -2.5 log FV - 5 log d10 + CV + AV = V - 5 (log d - 1.) + AV = V + 5 - 5 log d + AV. \n", "As the distance is often expressed in kilo-parsecs (kpc), i.e. 1000 pc, this is then rewritten\n", "MV = V - 10 - 5 log d* + AV.\n", "Using the parallax, varpi, in milli-arcseconds, instead, we have\n", "MV = V - 10 + 5 log varpi + AV, as varpi = 1./d*.\n", "Instead of looking into the V band, we can also look in any other band, and the same relation exists.\n", "\n", "One can use such absolute magnitude as a proxy for the total luminosity of a star. What about the temperature? We know that stars of different temperatures have different colours: the hotter stars are bluer, while the Sun with a surface temperature of about 5,800 K is yellow, and cooler stars are red. Such a colour can also be seen as the difference between the magnitude obtained in two different bands, a blue and a red one. Thus, if a star is bluer, its magnitude in the blue band will be smaller than in the red band, and the difference between the two magnitudes will be negative. The opposite is true for a red star. \n", "\n", "The Gaia satellite provides photometric observations in three bands: the G band, which is very wide, and in the blue (Bp) and in the red (Rp). One can thus create a colour-magnitude diagramme using Bp-Rp as the colour and an absolute magnitude based on the magnitude in G.\n", "One needs to also correct these magnitudes for the extinction. This is done with:\n", "A_G = 0.86117 * Av\n", "A_Bp = 1.06126 * Av\n", "A_Rp = 0.64753 * Av\n", "\n", "Let us for example look at the CMD for the open cluster Messier 7. As this is a cluster, it is defined by a single age and chemical composition and is thus representing a single population." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAAD8CAYAAADkM2ZpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXl8VNXd/99n9snMZN8JSdgF1IAEZFNBQEFaQBFFfVrQitblZ9Xqg9rF6mOtPvrUWlttcSmLtgpKFUQBsdadTQRl35JAEiB7Mslk9vP7I7m3k5BAyEw2uO/XK6/M3LnLueee87nfs32/QkqJhoaGhkZTdF2dAA0NDY3uiCaOGhoaGi2giaOGhoZGC2jiqKGhodECmjhqaGhotIAmjhoaGhot0OniKISYKoTYJ4Q4KIR4qLOvr6GhodEWRGfOcxRC6IH9wBSgENgC3CCl3N1pidDQ0NBoA51tOY4CDkopD0spvcCbwMxOToOGhobGaTF08vV6AUdDvhcCF4fuIIS4DbgNwGazjTjvvPM6L3Ua3RYpJUII9XswGESn6/ou85ZaXqHp1Oh+fPPNN2VSyqTT7dfZ4thSqWlSuqSUi4BFALm5uXLr1q2dkS6NboqUEq/Xy+9//3suvvhiJk6cyPHjx3niiSe499576d+//ynFqLl4CSGabNu3bx/r16/n9ttvx2w2n1Ha9u7dyzPPPEN9fT1GoxEAm83GwoULyczM1ESymyKEKGjLfp396i0Eeod8zwCKOzkNGj2M4uJiXnvtNdasWYOUkl27dvHmm2+yfft2oEEAW/sLBoOUl5fjcrnUfX0+H6Wlpfh8PmJiYhgyZAh6vf6U5wn9U0hPT+fWW29lwYIFXHvttRQVFeHxeIiNje2SfNKILJ1tOW4BBggh+gBFwFzgxk5Og0YP48iRI1RVVXHgwAHcbjdbt24lEAiQn58PgN/v58svv+SLL74gOTmZmTNnkpCQwObNm3n33XfZsmULN910EzfddBPr1q1j7dq17Nmzh6effpqUlBSio6PR6/WUl5ezatUqiouLGTVqFBMmTGDHjh14PB727NlDRUUFt956KwkJCQA4HA5Gjx4NwAcffIDRaOSXv/wl0dHRXZVVGhGkUy1HKaUfuBtYB+wBlkspd3VmGjR6Hrt372bIkCHU1dVRWFjIN998w8yZM8nPzycYDPL111/zhz/8gdjYWD788EPeeustDh06xMKFC4mOjmbu3LlYLBbWr1/PM888w4gRI7jyyisRQvD111/z5ptv4vF4eO655/j2228xGo387ne/o6CggPfee4+7776bgwcPsm7dOg4dOqSmS2k2O51OFi9ezA033EBWVlaT3zR6Lp3eoy2l/EBKOVBK2U9K+dvOvr5Gz0FpAm/atInp06djsVj46KOPqK2tZdasWRQXF+N0Olm9ejU5OTmcOHECg8HAiBEj8Pv9AOTl5ZGcnMysWbPw+/0EAgH27t3LZZddxkUXXYTT6cRkMlFUVMSOHTu44IIL2LNnDyNGjCAxMRG/38/MmTN57LHHeOCBB8jMzDwpnd988w11dXWq4GrCeHbQ2c1qDY0zoqysjLy8PG677TaKiop48cUXmTJlCsOHD+fFF1/k6NGjFBcXc+LECa677jpuv/127HY7FouFv/zlL2zatIkXXngBm83G1KlT6dWrF59++imPPfYYf/3rX/F6vRiNRpxOJ3l5eWzbto158+YxfPhwdDodQgj69++P2Wxm2rRpJ6XP7/fz7rvvcsUVV5CUdNoBUI0ehCaOGt2a4uJigsEgWVlZjBkzhmXLlnH55ZeTnJyMw+Fg9+7djBw5ks2bN5ORkcEnn3zC5s2bycnJwePxkJmZSVRUFIWFhSxevJjY2Fh69+6N2+2mpqYGk8lEMBgkLS2NgQMH0r9/f3w+H88//zzp6elN0tKSRVhQUMD333/PggULNKvxLEMTR41ujcPhYMaMGSQkJHDJJZfwy1/+kvHjx2O1WvnBD34AwPz584mKiuL9998nJiaG66+/HpPJxNKlS/niiy8YM2YMV111FZs2beLdd9/F7/dz8803M3DgQEwmE06nk6SkJB577DGWL1/Ou+++y+DBg5k+fTr//ve/ycjIaDV9brebH/zgB/Tv37+zskSjk+jU5YNnSlfMc2xvfigWQ3uOj6S1cSbXbz7nLzQtrU1uPt32SN+LMh1Hr9er25RrKP2KBkPDOz4QCKhNYeV7IBDAaDQ2OUZKqc5LbOmaodcLBAIIIZqcN5RgMNhk/3DKQSSIRP6fKu2tlYHT/daW35vv215aK9Mh37+RUuae7jya5dgCgUBA7bz3+Xz4/X68Xi8ej4f6+no8Hg8+nw+fz4fX68Xv96srNvR6PUajEbPZjNlsxmQyNflvsVgwGo3o9Xr0en2HNcOUSq7ci8/nw+PxqH9utxu3263en3IPwWAQoMm9GI1GoqKiiIqKUu/LaDRiMBgwGAzqfSgrViJ1T8p5lDmIgUAAr9er/rndburr69VnotyHci/KOUwmU5PnYDKZMBqNWCwWrFarej+KACpiCzT53BI6ne6UK3WUPA0tL6FlyeVy4fV61eekCL6S/yaTCYvF0qT8KOlVnoFSniKNUoaU8uP1etVyFHoPyj0pLx7lBda8DClpt1qtTZ6DXq9Xy5HyDJqviGpv2sMpk5o4tsCKFSv4+OOP8fl8BINB9U2k0+lUQVD+lIerPNBQMfL5fKr1EggE1HMYjUZMJhM2m43Y2FhsNht2ux2Hw6GKpyKwSqWFpoVVEe66ujpcLhd1dXXU1taqBVb5Uwp0IBBQC5xSGEMFTqnkodcKrbDKC0G5D6XwK+dRBNRutxMVFaWKjiJAofcV+mJQ8k15ASl/9fX11NbWUldXh9PppKqqirq6Ourr61UBUYRMyU/lvAaDoYkFpzwL5UUQCARU0RJCqGm32WxER0djtVqxWq3qCyFUQENfBsrLJ1Toamtrm/wpAq6Id2hZUtIdmieKyIW+2BThUdLe/BxWq5WYmBji4uKIiYlRB6SUl0HzPFfKUOizVe5BSXd1dTUulwu3260+k9AXTmgZCn2uzZd4KmW/+X2EtjSUcynlRSkzcXFx6j0pghqaP0o9U9Kt1IOamhqqqqqQUnL77bczYMCAdumAJo4tMGzYMHQ6HTk5OURFRTUpxKGFoLm1pDxwRcRCC4fy5lUeYH19PTU1NVRWVqrCVlZWplYipUAqFVm5hiJgyhtZqcQ2m42kpCT1zWyz2bDZbFgsFqKiotRCZzKZTrL4QgcSmjcLlZeDcg+h1o9ivbndbjweD7W1tTidTlwuF/X19dTV1VFeXt5EpEPFKdTKCLUwlPtSXhjp6enExcURHR1NdHQ0NputiRWoiHvoPYTeR2jzPPSZeDweXC4XTqeTmpoa9U95sZSXl1NYWKimP9QqbZ7u0OfgcDhIS0vD4XCo2xTBCn0GitiGPoOW8j803c0tOJfLRU1NDRUVFVRWVlJVVUV5eblq3Sn7hrYKQl9ser1ebQ1YLBbsdjtpaWmcd955ahlSypOSdkWkTtdqaK0+hJYjJZ1ut1stN8r/yspKjh49yq5du9RWTnOBVixSm82mvpzT0tKIi4sDoLq6ut1WqCaOLeD1evn888+ZNWvWGa23DX0AZ9rMCe0nCS1UoQISep3mYtb8c7go51IKvVIIz4TW7qn5b8r1OuKeItF3daq0t5TuSDyH5vnfVtqa56dKc6S7ReDM6wOcXA9Ctynnb63cQMP802XLlqnTss4UTRxbwOPxqJ37nUVHiVxX0tPvKZJi11n09DwPJdz8d7vdatdJe+h6n0/dEGVicE8vXBoa5zLh1mNNHFtAydTu4C9QQ0OjfXg8HkwmkyaOkcTpdGKz2dqcqc1dWWloaHQtUkq1HrcXrc+xBaqrq0/rdkoRw7o6SUWFxG6H2FgdQgACRIt+fTU0NDqL6urqsHxrauLYAhUVFarPvlNRVxdk2zYPZosXr9tE3/5G0tP0IEXLPs81NDQ6jfLycnr37n36HVtBa1a3QHl5eZvEsaoqSDDo5NN/LyIhsZbjx4KdkDoNDY3TIaWksrKyTfW4NTRxbEYwGKSmpqZN5nhMjA4poxg9+jq8XhthdG+chNaPqaHRfgKBAE6nk5iYmHafQxPHZgQCAdxuN3a7/bT72u2C4cOtBALJBAJ6srIMRKI9LaVkz5495OXlhX0uDY1zEWUFVDgDMu0WRyFEbyHEJ0KIPUKIXUKInzVujxdCfCSEOND4P65xuxBC/FEIcVAI8Z0Q4qJ2p7oDUZY3WSyWNoxWC2Ji9FxyqZXhw81YLDplc9gcO3aMsrKy8E+koXEOoiz1PNNVXaGEYzn6gZ9LKQcDo4G7hBBDgIeAj6WUA4CPG78DTAMGNP7dBrwUxrU7BMUBQiAQOO2ywdClSwa9QKcLWR8bpjoKIbj88ssZOXJkWOfR0DhXUeqxyWRq9znaLY5SymNSym2Nn500BMzqBcwEljTutgSY1fh5JrBUNrARiBVCpLU75R2E4nkmnDfOmaD0LSpLFhVC14xqaGicGV6vFynlGcciDyUifY5CiGxgOLAJSJFSHoMGAQWSG3frBRwNOaywcVvzc90mhNgqhNhaWloaieSdER6PByCsN86ZEggE+Nvf/sb+/fs77ZoaGmczHo9H9eXZXsIWRyGEHXgHuFdKWXOqXVvYdtJwrJRykZQyV0qZ2xUBi/x+v+rjr7PQ6/VcddVVYc3J0tDQ+A8+n+8kx8VnSliTwIUQRhqE8Q0p5crGzSeEEGlSymONzeaSxu2FQGjtzwCKw7l+R+D1elV/iR1FqPslpemsCaOGRuRQ6nE44hjOaLUAXgX2SCl/H/LTKmBe4+d5wHsh23/cOGo9GqhWmt/dCa/Xqzry7EgOHDjAihUrVJdKWh+jhkbk8Hg8YTuPCUcBxgE/Ar4XQmxv3PYI8BSwXAjxE+AIMKfxtw+Aq4CDgAu4OYxrdxh+v/+0cUEiQWxsLIMGDWoxIJMmkhoa4eH3+9WYNO2l3eIopfyC1mf0TWphfwnc1d7rdRbNI9h1FElJSWoQeGUKUW1tbVgL5TU0NBqIRD3WVsg0IxzPwW2lpZghBQUFLF68GK/X26HX1tA4FwiNNdNeNK88zQgNlNURtNZ8zsrKYv78+ZhMJq2JraERJi11V50pmuXYDL1e3yQca0fQUuAgp9PZZITc7XarIUg1NDTOjEjUY00cm6EEkFfM8kgTCAT4xz/+QV5eHocPH6a8vByAjz76iC1btgANovn3v/+d7du3n+pUGhoarWAwGMIWx3OqWR0av1ix0po3WzvactTpdJx33nk4HA42bNjAoEGDSEhIYMaMGU1i//7whz8kKiqqQ9KgoXG2E4l6fM6Io5SSmpoa/vGPf1BQUMCECROYNGmSGpRcwWQyqUHfOwIhBMnJybzzzjtceeWVZGZmAmC1Wpvsl5yc3NLhGhoabcBkMqnOJ9rLOdOsllKyYsUKKisrmTFjBv/85z85dOiQ+pvyhjGbzQSDwQ4dNbbb7ZhMJtauXQtogy4aGpHGYrGo7gfbyzljOXo8Hr777jt+9rOf0bdvXzZt2sTevXsZOHAgGzZsYMOGDbjdbkaNGqXu3xEIIYiNjeXHP/6xOuFcQ0MjslgsFqSUYRk554w4BgIBgsGg2nw1GAyqyf3tt9+yePFiDAYDCQkJCCE6fL6hXq9X+0WKiopISUkJy72ShobGf1Acx4RTj88Zs8VsNhMVFUVxcTFer5f8/HzS09MB+PGPf8z69evZsGEDN998MwaDAbfb3WFpCZ387fF4ePvttzWv3xoaEcRoNKLX68NqAZ4zlqPBYGDKlCm8/PLLJCYmIqVkyJAhAKSkpJCSkgI0uDoymUy4XK5OSZfVauWOO+5QrUal/1NK2SnLGBWklNTV1TUZMTcYDKoDjkikQ0pJSUkJJ06cYNCgQRiNRoqKinA6nQwcOPCkwTENjfZiNBoxGo1h1eNzRhwBJk6cSFpaGpWVlQwZMgSHwwE0rfh6vR6z2YzL5WriUqwjae51fMuWLRQWFnLNNdd0+LUV6uvreeGFFzhw4AA1NTUkJiaSnJzM3XffTWt+Nc90JU9VVRWPPfYYffv2JSsrixMnTvDoo48yduxY+vbt26Fu4jTOLfR6fdhGzjkjjkII9Ho9Q4cOPWl7KDqdDofDQXV1dYekQxGUYDCoWoahK2agYSlh83i7p5qvFQkBt1gsLFiwgAMHDvDCCy9w5513YrVacTqdOBwOvF6v6kDU7/dTXl5OWloax48fV63Mvn37otPp1LS63W6KioqIj48nOjqanTt3EggEmD17Nmazmc2bN5OSksL06dO1/laNiKLX67HZbNTUnMr/9qk5Z8QR2i4idrsdp9PZYemorq7mH//4BzfeeKMaV7eoqIgdO3YwderUJs38UI4fP05BQQEmkwmLxYLFYsFsNjf5MxgMJ91nW+5bCKEKssPhICMjA7vdzqOPPsqtt97Kzp07KSsrY8+ePdTU1FBdXc3111/PX//6VzIyMpBS8uyzzxIXF6fGxHn++efZv38/MTEx3HLLLbz88svk5eWxatUqxo4dy+LFiwkEAqxdu5Yf/ehHHe5DU+PcQQiB3W6ntra23efQSmMLhJupp8NisTB69OgmzWmPx6NaXK2JWVJSEjabDY/Hg9vtxu12U1lZicfjwePxqJad2WzGarUSFRWFzWbDbrdjtVqbNFtbE1CTyaSO1gcCAWprazEYDDidTvx+P3V1dcyfPx+z2UxtbS1CCO6//35efvllDhw4wPbt2/H7/QwdOpSCggKefvpp/vCHP+D3+5k7dy67d+/mrrsaPNddffXVmM1mbrjhhk7tX9U4N7DZbNTV1bX7eE0cW8But3dYsxoaRs6HDRvWZNu2bdtwu914PJ4Wlw0q3QIOh0PtK21OMBjE5/PhdrtxuVy4XC5KSkrIy8tTY3FHR0cTExNDdHQ0UVFRJ82zVAZFfD4fBoMBo9FIWVkZdXV1CCGw2Wz06tWLzMxMvvvuOwYPHkx2djZSSiwWC5dffjk6nY7q6moMBgM+nw+fz4fVasVqtar3Af9ZFaQJo0ZHoFmOHYDdbqeoqKhDzt3cJZryefLkybz22mscOHCAnJycVo89FYrnY7PZrDbXFRRnujU1NZSXl1NcXMyQIUNOEmKDwUBOTg5WqxWDwcCIESP44x//SF1dHVdffTXJycnqHDKbzUafPn0wmUxkZmZiMBjo378/0DDAk5aWxiOPPMKFF15IdnY2hYWFTQKXdWYIXI1zD4fDQXFx+8NUiXAdLAgh9MBWoEhK+QMhRB/gTSAe2Ab8SErpFUKYgaXACKAcuF5KmX+qc+fm5sqtW7eGlb4zRUrJ2rVr+eyzz/jtb3/bIStYlGkzq1evZvr06aolWF9fj9FobLHfMBLXbInm11G8kiv9f36/n8rKSqChOyAYDGK329UJ7G63m6ioKGpqarBarU38UXo8HiorK4mPj8dkMnHixAnq6+vJzs4GID8/H7PZTFpammY5akQUKSWrVq3i22+/5dFHH21SvoQQ30gpc093jkjU/J8Be0K+Pw08J6UcAFQCP2nc/hOgUkrZH3iucb9uSVRUFPX19R3q01Gv15OamtqkHzAYDPLxxx+H1U/SGs29j7cW0EsJS6v8bjQa1ZAODoeDmJgYtemtjAgCREdHq1ahcl6z2Uxqaqraj5mSkqIKI0B2djZpaWkRv1cNDaULSJmS1x7CEkchRAYwHXil8bsALgfebtxlCTCr8fPMxu80/j5JtNNcCJ0oHfoXKaxWKx6Pp8N8OkKDFTZhwgSioqIoKSnB6XTi8/n46KOPOqxJ315OJ6otbW9p/9YEWrMaNToCq9WK2+3uGnEE/gD8N6CoSAJQJaVUXFgXAr0aP/cCjgI0/l7duH8ThBC3CSG2CiG2lpaWnnRBRQiPHz9OZWUlfr+flStXcvjw4YgJpMViwefzdZg4NheF9evXs2vXLmJjY3n88cfp169fxAVfQ+Ncw2Kx4PV62+3XMZy41T8ASqSU34RubmFX2Ybf/rNBykVSylwpZW5rKzP279/PbbfdxsMPP0xVVRUlJSW8+uqrEfPBqHgR7kjLMZTrrruO3NyGLpCoqCj0ej3r16/XPIFraISB0WjE7/e3ux6HYzmOA2YIIfJpGIC5nAZLMlYIoYyCZwDKcFEh0Bug8fcYoKI9F96xYwfZ2dmYTCY2bNjA+PHjOXz4cMScRSiDDZ0hjkIITCYThYWFLF++XBX41NRULUyrhkYYKCFPOr1ZLaV8WEqZIaXMBuYC/5JS3gR8AlzbuNs84L3Gz6sav9P4+79kO1MdFxdHZWUlEydOZM2aNezatQuz2Ryxtbld0QfmcDjo27cvQgjKysqorq6md+/enZ4ODY2zBWXaXHvFsSPmOS4E3hRCPAF8C7zauP1VYJkQ4iANFuPctp6w+c2NGTOGTZs2sXTpUvbv309xcTH33XdfxObMdUVfX0JCAvX19WzatAmv18vbb7/N6NGjT3mMFsJVQ6N1FMcx7a0bERFHKeW/gX83fj4MjGphHzcwJ4xrUFNTg5QSk8nEfffdx5EjRygtLSUhIYG+ffu299QnEQwGO3UUVXnDKcsBzz//fG666aY2rTUOBAKaNxsNjRYItx73mBUyNTU1PP744+Tn56vz6aKjoxFCUFlZyahRo7jvvvsi4rxACV/Q2SEM+vfvr64wUQJvnco69Hq9LF68mKuuuoqMjAzNetTQCMHv94flI7THiKPD4eCBBx6gqqqKzZs3s2TJEubOnYvdbufDDz+koKAgYs3hcDO1PSjXUqYpKZOma2tr2bZtGxdffHETt16Km7Dx48cTFxfXaenU0Ogp+Hw+dUlte+gR4qisxkhPTyctLY3vv/+e4cOHM3PmTIQQxMXF8dRTT+F2u5us3W0vXq8Xg8HQZcGvduzYAcCVV16Jy+XijTfeID09ncLCQs4//3wSExMBWvRPqaGh0YDX622y2utM6RHi2JxBgwbx1ltv8fnnn5OWlsbq1avJyMiI2ICM2+3GZDJ1mThOnjwZaHgpxMfHk5OTg9vtpqSkJKxQkxoa5xJutxuz2Xx2W47NOf/887nllltYunQpdXV1pKamcs8990TMWWp9fT0Wi6VL+vCU5jKgxpHJzc1l+/btxMbGamuRNTTaiMvlCstg6pHiWF5eTmxsLPPmzVN9BdbV1eH3+yPSrK6vr1d9D3YVUko+//xzUlJSGDVqlLqkEFqettPSwI2Uskm/izZgo3GuIKXE5XK16Bu1rfRIcdy7dy9//vOf8Xg86HQ6KioqcDgcLFu2LCKDE4rl2NU4nU4sFgvffvst5513nirYwWBQfRF4PB727t1LdnY2//73v5k0aRJ2u13db/HixYwbN06NtKihca6gGDntpUeK4/jx48nNzVWX2m3cuJFFixZhMplOeZyUkmAwiNfrxWw2t9pRG26mRgqn08mJEycoKCggOjqafv36AXDkyBHWrFnDggULcLlcbNy4kfT0dLXzWUGn03HllVcSHx/fVbegodFlhFuPu2bEIUwqKio4cOAAx48fp6ioiIKCAvx+/ymn8kgpKSsr47nnnuPBBx9k6dKl+Hy+Fo+pqqrqFuuaL7/8cgoLC4mLi1PdmEkpiY2NJTs7m3Xr1lFVVcWCBQtISkriqquuwmazNRHIzMxM7Ha7eqzm6UfjXKGqquokj/hnQo+0HLdv386f/vQn1WN1VFQUt9xyy2n7FwoKCujduzeTJk3i2Wef5YorriAtLQ2/39/Ee0dVVRUDBgzojFs5JYmJidxzzz0YjUb++c9/snHjRpKSkvjmm284duwYFouF4uJi5s2bp1rCoYT2PdbV1bFmzRquuuqqVmPQaGicLSgr6sIxcnqkOA4aNIhf//rX9O7dG4PBgNVqbVMf4YgRIxgyZAhvv/02ycnJ6ltlxYoVvPnmm7jdbi6++GIqKiq6fGK1ImwxMTFIKYmLi6O6upr09HTS09O5+uqrKSsr4+9//zuvv/46t9xyyykHXHQ6HQkJCV02PUlDozMJBAI4nc5zw3JUmoPBYJC1a9eSn5/Pr3/9a9W1mNfrPW2fY3V1NX/605/w+Xzcf//9qqVZUVHB4cOHMZlM6PV66uvrw8rUSBEqdldddZX6WVlamJaWxoIFCygtLVXXkbZ2vNVqZdKkSR2cYo2zmdBume7uwd3n8+Fyuc4NcYSG5u7//d//8cknn1BeXs7hw4fVkduhQ4fy0EMPqTFNWuKTTz6htraWu+66C4fDoXrtmDdvHnPmzEGv1+P3+3nwwQdxOBzd5uG3lA6/309ZWRkJCQl89dVXFBUVMWzYMPbs2YOUkpEjR6rzJdtyH6EDVVJKAoHAGR2vcW6wZcsWvv/+e+bOnXtS/3Z3wuPx4Pf71f729tCjxNFsNnPppZdSU1PD0aNHWbBgAWazGYPBgM1ma7L2uCWCwSBHjhzh6aefVvvz4uLisNlsqqiWl5cTDAa7xVSeU1FTU8PSpUu58cYbOXLkCJ999hn5+fnqtJ4RI0accl5kcyoqKliyZAm33nor5eXlrFu3jp/85CcRmTeqcXYQDAbZsmULHo+H9957j+uvv77TfRC0Fa/Xi5TytJpwKnqUOFqtVqZMmYLBYKC0tJTLL7/8jI6fPn06l156qVrhlYGJ0IELv78h/E2kVtt0FLGxsdx5553odDp69eqF1WolLy+PefPmkZ2djcvloq6ujvT09DadLzo6mmuuuUYN0zB58uRunwcanYcywFFXV8esWbNYu3YtmzZtYuzYsV2dtBbx+XxNVpu1hx5V+vfv34/f72fkyJHk5eWxd+9ejEYjwWCQ+Ph4EhJOitelIoTAbDaf9CZp/taLRKZ2BkroSYD58+cTDAapqKjg0KFDLFu2DK/Xy+7du3nuueeaWH+tveWNRqMaNtVgMKiu0zQ0FI4dO0afPn3YuXMnw4cPP20ff1dyzonjvn37cLlcADz11FM4nU51IGLWrFmnHbFti/nv9/t7jDiGotPpSExMxGw2k56ezssvv8z48ePR6/WUl5fj8/lITU1t8/k0NJqTnZ2s23HRAAAgAElEQVRNZmYmX331FS6Xi0suuaSrk9QqPp8PnU4XliPosBRACBFLQ8zq82mIJHgLsA94C8gG8oHrpJSVjTGqnweuAlzAfCnltjO53vTp09XPL774IvX19fh8PnV9dSTwer3o9foe511b8SbucDhwOBzExcXx9ddfc+WVV7J8+XKSk5OZPXv2addnBwIBhBDaWmyNk1DqWGpqKtu2NVTd7lpGlHocztS1cM2j54G1UsprhRAmIAp4BPhYSvmUEOIh4CEa4spMAwY0/l0MvNT4v00oPh2llOTl5bFo0SKKi4vx+/3U1dWRnZ3N448/HvYUHK/X2yVewCNBaN9pv3796Nu3LzU1NQBMmTKFvLw8bDYbycnJ6kh9KFJK3n77bfr27cvIkSPVbS1dQ+PcIrRspaSkUFFR0WRGQ3dD8cnaJZajECIauBSYDyCl9AJeIcRMYELjbktoiC2zEJgJLG2MOLhRCBErhEiTUh4702uvXbuWgwcPcscddyCEICoqirS0tLCG7RWUmCw9URxD+eEPfwg0FOrbb7+dnTt3cs8993DDDTeQnZ3NmDFjiI6ObnKMEILRo0eflI8ejwcpZbcfwdfoHBwOBz6fD7fbHZE61xFEwpt/OArQFygF/iaE+FYI8YoQwgakKILX+D+5cf9ewNGQ4wsbtzVBCHGbEGKrEGJraWlpixfOycmhrq6OYDDIyJEjGTFiBL169YqIoLU0mbqnoTSLlaaxTqcjOzub//3f/6W6uppnnnmGzz77DI/HQ0VFBYWFhaqFmJWVddLA1meffcaHH37YFbei0c1QBjZNJhNOp7PbrtWPRJC8cGxiA3AR8P+klJuEEM/T0IRujZZSeVLOSikXAYsAcnNzW8z5o0ePkp+fz09/+lOGDh1KQkICQ4YM4Z577glrXhP8p+/ubCM6OppRo0aRkpLCnDlzWLlyJS6XS5239stf/pLt27dz2WWXnZSHo0ePVtedw3+8GzV/GfX0l4pG29DpdMTExFBRUXHKQb6uJBL1OBxxLAQKpZSbGr+/TYM4nlCay0KINKAkZP/QKPUZQHF7Ljx27Fj++te/4vf78Xg8+Hw+MjIyMBgMYcdyPhsreGh/UVZWFlJK5syZg06n4/Dhw+rAVmlpKSdOnFCXJyrHNm9+19XV8frrr3Pttdeyf/9+LBYLw4cPb7Efs7vSUsVpLe2tlanT9ceezf218fHxVFRUdHUyWkURx3AEst3tUCnlceCoEGJQ46ZJwG5gFTCvcds84L3Gz6uAH4sGRgPV7elvhIYJ0OvWrWPZsmV8/vnn7Nu3j23btrF+/XpKSkpOf4JTYDAYCAQCTSylswWlmaHT6cjKyiIxMZEhQ4bws5/9jISEBC688EJWr16t+slsDZPJxOjRo7HZbFitVoxGIwcOHMDtdnfSnUSGgoICdbJ8MBhsUpGUiqVs83q9uN3ukypbfn4+1dXVrV7jyJEjVFVVdcwNdCHx8fFUVlZ2dTJaxWAwnPRMz5RwO+n+H/CGEOI7YBjwJPAUMEUIcQCY0vgd4APgMHAQeBm4s70X3bFjB8XFxcycOZPY2FjeeecdvvzyS9555x0eeOABTpw40e4bUkbEz8amdXMsFgszZswgOjqaV199FZ1Ox+TJk1H6ekMFIjQ/jEYjOTk5WCwWhg0bRv/+/Vm3bl23riwKyr0EAgH+9re/UVBQwPLlyzl+/HiT3wF27drFvn37ANi2bRuffPLJSfu8++67bN++vcn20N/XrFmjnuNsIjY2tluLviKO4Rg5YY3DSym3A7kt/HSS+5fGUeq7wrmeQiAQICkpienTpzNjxgxycnJYu3Ytjz76KA888ADfffcdU6ZMade59Xp92JnaU1CaeXq9nlmzZhEfH8/HH38MwBVXXIHRaKS2thabzYZOp2uxmaiMYv/0pz/tttM6FJS+0kOHDlFSUkJpaSkWi4WxY8eSkJDA8ePH2bhxI6WlpVxyySUsXboUgJtvvpmsrCz69OlDMBhk7969lJaWMmzYMAKBAFVVVezYsYOsrCwcDgcHDhzgyy+/JCMjA4/Hc1aWpZiYGGpqarptV4pSj8Mxcrp3aW6F8847j1deeYVFixaRk5PDtm3bMJvNREVF4XA4wmremUwmAoGAusb6XECn05GRkYGUkiuvvJJFixaxb98++vTpw69+9SsmTJjAlClTsNlsJ608UCpGT3FQ8dVXX/HOO+9gt9vJz88nGAyyZMkSFixYwFNPPUVKSgomk4ni4mKOHj2Kz+fj+PHjbNq0ibS0NFJTU1m0aBG9e/dGCIHH4+Hll18mKSmJiy66iHHjxvH0008zbNgwdDrdabsoeio2mw23200gEOiW097MZjN+vz+s/O9+d9UGUlJSePjhh8nLy+OFF16gtLSUm2++GSEE/fr1C2tdsOKy61yMDy2EwGg0cuONNzJw4ECsVivTpk1j27Zt7Ny5E5/Px4oVK9i8eXNXJ7VdKBEd/+u//osHHniA/v37EwwGqa6uRqfTYbFYqKioYPTo0UyYMIHZs2czZcoUJkyYQEVFBcFgkPXr13PNNdfw4IMPMm7cOAKBAJdeein33HMPx441dKEbjUZcLheXXnppj1tp1RaU6Tzd2Ygwm82qG7720iMtx9raWvbs2cO4ceMIBoMUFRWxYsUK7rvvPm6//fawCqQyjaWnDS5ECiGEGpBLsSTT0tL41a9+xbXXXqvOCKipqSE6OrrJSLhCdXU1er1ejYLYVbQ0yqzX6ykoKMBisVBVVaVWbikl8+fP59ChQ7z55puMHTsWk8mkviSVc7ndboLBIMePH0dKicFgIDc3F7vdTl1dHXa7nfvvv5/Vq1ezdu1aRo0aRe/evTnbMJlMquXcHYLRNcdisSClxOPxtPscPVIcN23axAsvvMCAAQOwWCyUlJSQnp6OyWQKu3kX6g38XEfpU8zMzCQ5OZnHHnuMqVOnsnjxYv73f/+XnJycJvvX1NSwZs0a/H4/KSkpXHHFFUgp8Xq9GI3GTmt+SSmprq4mPz+fzMxMNeSFEIKrr76aV155hU8//ZSMjAwcDgfjxo3D6XSyePFi6urquPjiizGbzcTFxaniOWzYMFJSUkhLS2Px4sWsXr2aq6++moEDB5KZmUl8fDy9e/dm+/btbNiwASEEU6ZMoaqqil69Tlrr0OMxGAyYTCbq6+u7RTC65phMJnQ6XVhGTo8Ux+rqanJycnj66acxmUwcOHCAJ598kvr6+rDEUfHGYzQaz1nLsTlCCGJjY5kyZQrBYJB9+/YRFxdHfX09+/fvJz09XV1CZjAYSEtLY9iwYao14XK5WLRoEXPnziUtLa1D06oEEjt48CBOp5O+ffs2WWsvhGDAgAE8/vjjBAIBTCYTBoOBG264AYDf/e53+Hw+7HY7er2+yeT3iRMnqtbnb3/7WwKBAHa7XR2QEEJw1113YTKZmDp1KjqdDpvNpjpAONtQuiEUL1ndDaPRiF6vP/fE8aKLLmLx4sW88sorjBo1ik2bGuahR2K0VIlmWF1d3W1H4rqCOXPmcOmll3LvvfeSmZnJ/fffj91u52c/+xnTp09X/UtOnDhRPUZpdlZXV1NdXd0h4qg0dz0eD4cPH6a0tJSsrCzOP//8FtfWCiFOWiOu7NM8xEZrfjBbC8WhxCQKFeSzeT262WxW5352t3piNBqxWq04nc52n6NHimNWVhaPPvoor7zyCv/617+Ij4/nrrvuikjfh16vJzo6ukfM2etMhBAkJyfzi1/8grVr13LrrbfSu3dvjh49qoa11el0GAyGJhXFYDAwY8YMkpOTT3H2M0cRRb/fz5EjRygsLCQ1NZUxY8aootbdKuzZhtVq7bbdT3q9HofDEdZczB4pjlJKkpOTufHGG8nLy6OsrIwtW7YwePDgk5a6tYfY2FhNHENQREan03HBBRdQUlLCW2+9xbhx4xg/fjxffvklq1atYsyYMVx77bVNLAmdTsfw4cMjmh5lknVpaSn79+8nOjqakSNHqi9HTRQ7B4vF0m27n4QQxMTEnHviuHnzZh566CGSk5OJj49Xl8FFahJyXFycJo4toAzQ5Obm4vf7Wb58OYMGDWLChAnodDq2bdvGvn37GDRokLq/giJo4QiXYi3W19ezd+9ePB4PQ4cObRJjXBPGzsNoNHbrKW/h1uMeKY4nTpwgNjaWP//5z8TFxakrHyIV08LhcKjLyTROJiYmhilTprB+/XoKCwtJSEhgx44dxMfHN+lvUwQxGAxSUlKCx+NR49ScqYgp5yosLCQvL4/MzEwyMzPVEXBNFDuf0KlO3ZFw63GPFMd+/fpRUFDA448/TlRUFOXl5eh0Oh5++GH69u17ymOVOXqBQIC4uLgWK5XNZuu2o3BdTWhz+fbbb6empobKykoqKyv54Q9/iNPpJCkpSbXid+zYwd69e+nduzeVlZWqOLYVxVqsq6tj9+7dCCEYMWKEOiiiiWLXYbVaKSsr6+pktIrNZqOurq7dx/dIcdy5c6caIS8pKYn4+Hji4uJO2+kvpaS2tpZf/epXZGRk8OCDD7a4nyKOwWDwrJyGESkGDhxIbW0tf/zjHzly5AiHDh3i3Xff5cknn1Sn98TGxuJ0OsnNzVUnDreVUGvx8OHD9OvXj169eqnn0ISxa+nOAzLK7IlwjJweKY4XXHABKSkpWK1WJkyYgNlsRq/XtzrFAv7jdODdd9/l2LFjJCYmqttdLhe1tbX4/X68Xq+67EgTx9ZR+h/1ej2pqans37+flJQU7rvvPoxGI0uWLOGSSy4hLi5O7RNs6yTw0L7FPXv2EAgEGDlypDpVRhPF7oHRaOy2ywehQbwVxx/tqcc9ThyllHz11VeUl5fz4osv8t5776HT6Rg2bBi/+MUvTjmd5/vvv+eLL75g9uzZHDx4UN2+ZMkS/vKXv6hOFe644w4CgcA54bYsXKxWKzfffDMzZsxgy5YtrF+/njFjxvDmm29iMBi46aabuPbaa086TkpJRUUFZrMZm83WZBmilJLi4mIOHjxIVlaW1rfYTTGZTGEtz+toFPFubz3uceIohGD27NlcdtllmEwmgsEgHo+HuLi4U4ZI8Hq9vP766xQXF7N69WpOnDjBNddcw+DBg4mNjaVPnz6YzWYSExNJSUlh+/btmjiehlChSkxMZNCgQWzZsgW/389tt91Genp6q8ceOnSI3/3ud8yfP5/x48cDqA4/du/ejdvtZsSIEWrzXBPF7odimXXHSeDQVBzbk8YeJ47Q4IVYcY7QVoxGI/feey9Op5MtW7awdetWtfJed911zJ49Ww1KtX//fj744IOz1t1UpFGa2Pv27SMnJ4cZM2ZQWVnJypUrKS4uZuDAgQwZMkTdFxpWV9x0002MGjWqybkUX52pqamatdjNsVgsarO1O7otM5lM6gKF9tDjxLGlwDltqTxCCNUBQCAQwGg0qtNODAaDOrqq9KOdK97AI8nkyZPVPPvggw/YsGED+/btw2Kx8Itf/ILp06er+2ZkZJCRkQE0fX4Wi0V9aWmi2L1R3JZ11/jVioPm9tbjsOReCHGfEGKXEGKnEOIfQgiLEKKPEGKTEOKAEOItIYSpcV9z4/eDjb9nn8m1mrufdzqdfPfdd9TW1rYpA0LDNA4ePFjtB9MqYGQQQmAymdQR6fj4eNWpw913343ZbGbt2rXU1NQ0OSY0/5Xv4YbU1Oh4lOfdnX2fhhtkq93iKIToBdwD5Eopzwf0wFzgaeA5KeUAoBL4SeMhPwEqpZT9geca9zsjvvvuO9auXYuUkl27dnH99dczb9481qxZg8/na1MmCCHQ6/UYjUatAnYASp5OmzaNF154gV69evHBBx9w6NAh9u/fz4cffkh5eTnQcswVZVZBW5+nRtehuKHrroMyoR6T2kO4HQUGwCqEMABRwDHgchrCtAIsAWY1fp7Z+J3G3yeJNqZaeTu99NJLHD58GGhYWJ6YmMill17Ks88+ywcffBDmrfyHSAQEP5dR8k5xJ7dv3z6WLVvGyJEj+fvf/87ChQv59NNP8Xg8fPXVV+zcubPJ8bt27eL1118/K2OvnE3o9XrVp2N3JNyXazihWYuAZ4EjNIhiNfANUCWlVCY/FQKKp89ewNHGY/2N+yc0P68Q4jYhxFYhxFYlCh40zHk7evQogwcPRghBYmIiRqORadOmcccdd7Bs2bKwZsOHooljZLDZbCxcuJCXXnqJK664goEDBzJjxgzy8vL4/e9/z8qVKwGa9PcCZGZmMnHixG7Zya/xH4QQREVFqV1b3Q1loKjTLUchRBwN1mAfIB2wAdNa2FXJtZZSeFKOSikXSSlzpZS5SUlJ6nabzca4ceN48cUX2blzJ8nJyZx//vkUFRVx/vnnU1NTQ21tbXtvpwmBQEATxzBR8i8qKoqRI0cycOBAPvzwQxwOB3PmzOGKK67AZrORm5tLfHx8k6lT0dHRZGVldfEdaLSFmJiYU8bt7krCrcfhvJonA3lSylIppQ9YCYwFYhub2QAZQHHj50KgN0Dj7zFARVsupPQT3nnnneTk5HDffffx3HPPcdlllyGE4N133yUhIUGdExcufr+/RUepGu2nX79+FBUVUVZWRllZGXPnzuXIkSOUl5dTWVnJgQMHgAbrccWKFU0m6Wt0X7qzez+fz4fBYGh3CySc8fcjwGghRBRQT0Os6q3AJ8C1wJvAPOC9xv1XNX7/uvH3f8kzsMUV/2wLFy5k6tSpvP/++7z11lt4PB7S0tJYuHDhKZcPngmKOGrNusgghOCiiy5S3cpt2rSJe++9F71ezw033MCAAQMYOHAgQghqamrIyMg4KS5Je6ZvaXQ8CQkJ7Nmzp6uT0SKhRk57yku7xVFKuUkI8TawDfAD3wKLgDXAm0KIJxq3vdp4yKvAMiHEQRosxrlnek3RGDo0NzeXESNG4Ha78fv9WCyWiI4+KwGhtAoYOYxGo+qhOzMzk+LiYu655x5WrlyJzWZj7tyG4vDmm29y7NgxLrroopPOUVBQQFlZGSNGjOjUtGu0TmJiIhUVFd1ylYzH4zljZyehhDVzU0r5KPBos82HgVEt7OsG5oRzPfiPxaD0Z3VER7Db7Vajl2lEnszMTJ555hlee+01Vq5cybRp0xg/fjy9e/dm1qxZ5OXltTipuK6ujoqKNvXEaHQSsbGx1NfXd8sQrW63G7PZ3GVTeboURRgjPeVDyVSNjsFgMJCTk8PFF19MVFQUGzdu5Mc//jFbtmwhPj6e0tLSFkVwyJAhTJkypQtSrNESQgisVis6nS5ig6GRJNx63OPEUUqJ3+9nz549FBYW4vV6efXVV9mxY0fErMhwzXGN06PX65k0aRLZ2dlYrVZKS0v55JNPOHToEMXFxSdNy9JWz3RPDAYDMTExlJeXd6vpPFJKPB7PuSWO0OB67J577mHhwoWUlZURFRXFkiVLIuZbrr6+vts1Ec42hBCkpaXxf//3f1x77bX8+c9/JjExkXvvvZeSkhLNE3sPIjU1lWPHjnV1Mk7C5XKpPkDbQ48Ux71793LhhReSlZXFBx98wLBhwzh27FjEIqFp4tg56HQ6cnJyWLhwISNHjuSzzz5DCEEwGMRsNuP3+3G73fzrX/+irq5O9fNYU1OjrZ7pRqSlpVFcXHz6HTuZcOtxjxTH9PR0jh49Sk5ODuvWrePzzz/HbrdHLMBWdXV1k0BRGh2HEEJdxeBwODh27Bj79u3jtttu46WXXsLj8VBSUqK6nnrvvfd48skn2bp1K1999dVJ63qbr9XW6HhSU1MpLS3tdi+smpqasEI190hxHDVqFJMmTWL58uUcPHiQ999/n3nz5kVMHGtra3E4HBE5l0bbiIqK4rHHHmP27Nn069eP/Px8/vKXv7B//36uu+46oqOj0el0zJ8/n3nz5pGYmMj+/ftP8ghTUFDARx99pAlkJxIbG4vH4+lWMayVeFHh1OPu54StDdTU1HDFFVcwa9YsqqqqiI2NJXSpYbi43W4sFkvEzqfRNhISErjjjjt45JFHiI2N5fLLL+epp55i4sSJjBgxgtGjRyOlxOFwkJ6eTnZ29kmDMzqdrlv6FjxbUUaszWYzVVVVEVuIEQk8Hk9Y9bjHWI6hzaXVq1fz2muvERcXR3Z2NtHR0RHz2h0MBqmvrw+rI1fjzFFELj4+nmeffZa///3vPPDAA5jNZl599VVqa2spLi7m6NGj/POf/8Tj8bQ4cp2amkpVVRWhTks0OhadTkdSUhLHjx/vNtZ6MBjE7XaH1efYo16x1dXVvPDCC2zYsIHy8nKOHDlCIBDA7/czePBgHnjggbDfXEpMGk0cOx9F6BwOB4MGDSIQCHD11VcTFxdHYmIib731FnfccQfXXHMNr776KvPmzTup2aT0XWrWY+dywQUXqCugugOBQACv13vuiKPRaOTCCy+kqKiIwsJC5syZg9lsxmAw4HA4ItLnqDhaDWdmvUb7Cc1zvV7P+PHj2bt3L3/5y19wOBzs3LmTYDDIZZdd1mKTyWAwMHnyZPVcoZaM9jzbRnvWsV9wwQVNHBZ3dV4rRlM48xx7lDhGRUUxY8YMrFYrJSUlTWKSRAJlgnm4maoROVJTU7n66qs5ePAghw8fZs+ePXz99dfcddddaowQaJi2Ebp+W0HxQD516tQur7A9ASklR48eZc+ePaSlpTF48GAMBkOT0LlerxeXy0VsbKzatSGlpKSkhKqqKgYNGqTuG0pn5n8k6nGPEcfQh5OSksKOHTt46aWX1GBYMTExzJw5M+zmsM/nIxAIaPMcuwFKpRs6dChDhw5FSslLL73EG2+8wbfffsvSpUtJT0+nurqaVatWMXbsWC688MImx5vN5rCmc5yLLFu2jOLiYjweD1OnTmX27NlNhK6goIAtW7Zwww03NDlu3759bNq0iQcffFDdv7i4mO3btzNtWkuuXk8mUgLq8/kIBoNhDcj0GHEMRa/X43a71Vgk+fn5OJ3ONlmSp3ubKbFLIjUtSCM8mj+fkpISRo8ejc1m44033iApKYn8/HzuvvtuTCYTeXl5xMXFqS7PQqMcarQNt9vNj370I/x+PytXrmTy5MkUFRVRUlKCzWYjOzubMWPGIKUkLy+PwsJC+vbtSzAYxOVy8fXXXxMbG0tiYiLLly/n448/RqfTcdlll7Fv3z6Kioq48MILcTgcFBcXU1xcTFJSEjk5ORG7B6/XCxBWP2iPFEfFklDYuHEjTz/99GnfOkqzOT8/n6ysrBYzrr6+HiGENpWnC2lpxFN5tv/1X//FxIkT+eSTT8jLy+P7779n3LhxvP766+Tl5VFaWsqjjz7KRRdddNJ52tIXdq73UQaDQSorK1mzZg2HDh1i4sSJ5Ofn8+CDD3LxxRdTW1vLxIkTOX78OBdccAErV66kX79+vPPOO0yePJn169dTUVFBfX09CxcuxOv14na78Xg8fPPNN6xYsYKoqCg+//xzpk6dyu9+9ztyc3MBOP/88yM2kFZfX49Opzs3mtWh7N69m40bNxIdHY3P52PLli3Y7fZTZqwS1e69995j586dPPjggxiNxpMqkMfjUaMTanQeoc/B7XZTUlJCSUkJffv2JT4+Xv2tf//+9OnThxEjRvDGG29gt9tJSUnh5z//OS6Xi6uvvhqn08nXX3/NmDFjcDqdrFq1ipkzZ57WU7xSRmpqatDpdERHR59zAhkMBvF6vRgMBpKSkvj+++/p27cvw4YN46GHHuLYsWN8++23AHz00UfMmDGDrKwsdu7ciU6nY9CgQfziF7/giSeeICoqitmzZ1NcXMy0adP44x//iMvlwmw2M336dILBIGPHjuXnP/85J06cQK/XR+w+3G43BoPh7LccmwuY2+1m9+7d1NbWYrVaSUtL49Zbbz2ttbdr1y7WrVvH7bff3uRBlJSUcOLECfx+Px6PB51Op/ly7ARCn6vX66WkpITjx4/jcrlISEigb9++TTyCK0JlMBiw2+385Cc/YcWKFWzbtk2NgvfVV18xePBg1SGuwWAgIyOjzRbJRx99xPr165k4cSI/+MEPIni3PYNgMIjBYCA3N5eqqirWrVuH1+slOjqaqKgo+vfvz/bt21XH08ePH2fPnj34fD70ej2DBg0iOjoaKSUulwu9Xk8wGCQYDGIymYiKiuKGG26gqKiI+vp6YmNjsdvtEQtxoqCkJ5yXW48QR2i42YKCArxeLzabjdtuu019QEajsU2Zu3btWvLz83nxxRfp27cv999/P1arlSVLlvDss89iMpn47//+b00cO5DmglhWVsaJEydwOp3qpP64uLgmL6+WCrgQAoPBwDXXXENxcTGHDx/G4XDw3nvvsWzZMux2O0eOHOGSSy7hsssua3P6hg4dysCBA0lLSwvvRnsoisC9//77ZGZmsnDhQkwmE1VVVepzGDp0KC6Xi7i4ON5++21qa2tJTk4mIyMDm82G2Wxm1KhRqkPqCy+8EIPBwJw5c1i6dClvvfUWAwYMYMiQIWrffqQtdL/fH3Y9Fqeb0S6EeA34AVAipTy/cVs88BaQDeQD10kpKxvjUD8PXAW4gPlSym2Nx8wDftl42ieklEs4Dbm5uXLr1q3qNIFf//rXFBcXqytllNHMYDDIqFGjeOSRR1rtYwgGgzzyyCNMmjSJnJwcfvWrX/Hggw/Sr18/NmzYwMcff0x0dDRDhgzhgw8+4E9/+pM2KBMmLZWt+vp6ysrKKCkpoa6ujpiYGFJSUkhMTGxi3bW1soT69/zXv/6Fw+Fg8eLFeL1eRo4cidFo5NFHH1UdiZzqvK0N1p1L8WuklAQCAdWCVOpYIBA4yfo+evQotbW1fPzxx9hsNm6++Wb1t0AgoDoU8fv96rHBYJBAIIDRaCQYDCKl7JBgdps3b2bx4sU8//zzTaYiAQghvpFS5p7uHG2xHBcDfwKWhmx7CPhYSvmUEOKhxu8Lad7OA78AACAASURBVAjNOqDx72LgJeDiRjF9FMilIRzrN0KIVVLKNoctS0pK4plnnlG9syiZrAzZW63WU4qZEIKUlBSOHTumrslVputMnjxZnTj87bffntWFvyNpSQy9Xi/V1dVUVFRQUVGB3+8nLi6OzMxM4uLi2iWIoSithwsuuIDBgwdz5MgR8vPz2bx5M5999hmxsbH87W9/Y+7cuaSmpjYZlDkT0auvr+fIkSP079//rF59o1jkzbc1t8CklOzevZsvv/ySpKQktQsitOtDIbTfT6/Xq62CSPYxNkd5zh3arJZSfiaEyG62eSYwofHzEuDfNIjjTGBpY1TBjUKIWCFEWuO+H0kpKwCEEB8BU4F/tCWRyk0qfRk+n49PP/2Ujz76CLfbzbBhw5g1a9ZpzzNr1ixefvllvvvuO6ZNm0Zqaqp6/tBrdZf1od2Z1vKovr6e6upqqqqqqKqqwuPx4HA4iI+PZ+jQodjt9iYVLVIvIkUk+/Tpw8MPP8ySJUv47LPPMJvN/PnPf2bLli089dRT6rQe5bpKy2TkyJGnPH9VVRUbNmwgKyvrrBbHM+HKK69kypQp3dZDe7j1uL1POUVKeawxAceEEMmN23sBR0P2K2zc1tr2kxBC3AbcBg2BmFpi27ZtPPnkk/zoRz8iMTGR1atXc/DgQX7zm9+0aj0KIcjOzuY3v/kNgUBAXR4Y+kC7w7Kn7kpLBc3j8eB0OqmursbpdKpxRBwOB3FxcaSnpzcRw87IW51Oh9Vq5ac//SnXX389TzzxBDfffDMbN25kwYIF3HvvvUydOlW9n9Dycqp45Wlpafz0pz/tUGunJ6HkUXfND2X1VDgCGelXYEulX55i+8kbpVxEQ4hXcnNzW9yntLSU+Ph45syZg91uJzExkf/5n/+hrq7utE3r0/UjKqNr57r12Pz+A4EAtbW1qlXodDoJBALYbDZiYmJIT0/H4XBgsViaRIjsCpTrxsbGsnDhQmJjY/nnP//J8uXLiYqKoqKigj179jBixAh1Rc2xY8dYtmwZd955J3a7/aTWBKBZjD2IrhTHE0KItEarMQ0oadxeCPQO2S8DKG7cPqHZ9n+389oMHz6c1157jYULF3LBBRfwxRdfcOGFF0bEl5ySqd3Nq3FrtPbwm1vErf3W2j6BQACn00llZSWVlZXU1tZiNpuJiYkhOTmZ/v37Y7Vau1wIT4UQguTkhkbNtddey8SJE3niiSf461//SlVVFU899RRDhw5Fp9NRVVVFXV1dj3nuGqcmEkZOe8VxFTAPeKrx/3sh2+8WQrxJw4BMdaOArgOeFELENe53BfDwmV5UudH09HR+//vf8/7773P48GEuv/xyZs2aFZGJ2waDocdYjoqIl5SUkJiYCDRYQL169VJHChX8fj9VVVUkJCS0OiDhcrkoLy+nrKwMp9OJ2WxWB09iYmKa5G93FMOW2LJlC6mpqfTu3ZuEhATuuusuysrK+Pzzz3n22Wd55pln0Ol09OrVi+nTp/P2228zf/78bttc1GgbBoNBHXVv9zlOt4MQ4h80WH2JQojC/9/eecdHVaWN/3uml/SEHlIg9N4EFBAQAStgW7srrLJid9f28119d1131bWsbe2IviquItJEEFlYECH0UJIACb0lIT2TTD+/P2budRISCOmB++UzZObMnXuec+45z+nPQ2DV+UXgayHEDOAwcGPw8qUEtvFkEdjKczeAlLJACPE8sCl43V+UxZlzxev18tVXX5Gbm8uMGTPYvHkzGRkZ5OfnVzpJUVcU5dhQxnMbmxMnTvDoo4/yt7/9DSklTzzxBG+99Rbx8fGVlF9BQQHff/89t912W6WphbKyMk6ePEleXh4+n4/o6GjatWtHnz59zjgF0ViNR0MqXSklBQUFquEJg8FAjx496NGjB0OHDiUtLQ273c5rr71Geno6Dz74IImJiZW2imm0TppEOUopb6nhq8uquVYC99dwn9nA7HOSrhr8fj/Lly8nOTmZHTt28MILL9CpUye2b9/OO++8U+8z0SaTCb/f32BuXhub7OxsNmzYQHZ2NuXl5ezYsYPjx4/TsWNH9u7dy549e+jZsycdOnRg5MiR5OXlsXv3bkpLS1WjDEIIevXqxbFjx1i1ahVlZWU8+OCDZGZmkpWVRd++fenWrRtZWVls2bKFqKgo+vbtS35+Pn369GHr1q3069eP/fv3q/G1b9+e1NRU1cLR2LFjOXHiBGvXriUlJYVBgwaxb98+MjIy6N69O3369GlwZSSEYNKkSaeFAVgsFoYPH46UkosvvpiysjKys7P5/vvvGTBgABaL5bS5R43Wg9lsrncnp1UdA1EWVMaNG8eJEydYtmwZo0eP5rHHHqvWEXxdMJlM6nah1kBmZiY9e/YkMzOTrVu3MmTIEI4dO8b27dv55z//yZ49e3jvvffYs2cPX3/9NYcOHeKZZ55h69atzJkzh+joaFJSUnC73Tz//PN4vV5iY2PZtWsXH330ETt37uSNN95g9+7dvPLKK1RUVPD111+zfv16fvzxR8rKypgzZw5bt27l/fffJz09nddff53Dhw/zl7/8hTVr1vD5559z8OBB3n77bUpLS8nMzCQzM5N3332XjIwMXnvtNUpLSxs8b0K3mJxJyY0bN44pU6bw0UcfsXHjRp577jkefPBBTp48SUlJSaUK5vP5cDqdrWLa5UJG6eTUpx63KuWoMH78eHJycli5ciXjx4/n6NGjhIWFNYglHZPJhBDiNJefLRGPx8OuXbuYNm0av/zyC0ePHuXqq69m3759rFq1ikmTJnH11Vdjs9nU9Ljdbvr168d9992HXq9XTTuZzWZSUlLYvXs3HTp0YMeOHRQWFuJwOJg6dSppaWn079+fq666ivDwcKSUausMgc3z+fn5FBcXM23aNAwGA507d+bhhx8mMjKStLQ0rFYr06dP5/bbbycjI4O8vDyKioqYNm1as7mlUDY4Dx8+nKeffpoZM2Zw8OBBMjIyeOedd7jzzjvZvXs3K1aswOVysWfPHj755JNWM+1yoaKclFPKd11olcoxMTGR9957j/fee4+BAwcSFxfHgw8+2CAVzGg0otPpWoVyLCwspKCggDFjxpCTk0P79u0ZOnQo2dnZ+Hw+8vPzycnJobCwUC0kTqeTnj17Eh4ejs/nU9NZXl7OpEmTGD9+PIsXL0av12M2m7nxxhupqKigsLCQ3Nxctm/fjsfjoV27dhw+fJj09HQcDgc2mw2DwcANN9yA3+8nKyuLxMREdQHI6XRSXl7OwYMH2bRpEzqdDr1ez3XXXYdOp+PgwYPNlo9CCPR6PRMnTuTxxx/n2Wef5Z133mHkyJHk5+ezdOlSjhw5QmZmJrGxsVx99dXqrgaXy6W5gW2BGI3GendyWuXGLa/Xy6ZNm1iwYAFut5uBAwcyY8aMet9XOWVhMpkaZIje2AghGD16NF26dOGee+6hT58+JCQk0L17d8aNG8fixYs5efIkiYmJtGvXjssuu0xdtLLZbAwaNEjdoO10Olm6dCkOh4Nrr72WESNGUFFRwYIFC0hOTmbKlCmsWLGCefPmYbfb6devH+vXr+e7774jPj6eiRMn4vF4WLhwIQkJCYwaNQqPx4PRaGTIkCEMGDAAj8fDu+++S+/evbn22mspKChg8eLFxMfHM2LEiObMSlVB6vV61QWs3++nQ4cOzJw5k/j4eFauXMm9995LSUkJmzdvZuDAgXz33XfMmDGjTtbGpZRkZWURERGBw+EgPj5eO8/fQCgGacrLy+t8j7ManmhOFMMToUgpSU1N5a9//SsPPPAAcXFxzJkzh44dO/Lkk0/WewuG1+vl4Ycf5vbbb2fEiBEtekJeMQiguIqAQCV3u93qwX4FZWtP6HXKiRClF6ScW1cqaOj9lXxYvXo1mzdv5o9//KN6xl35TdXr/X4/Op1ONUIAqAoz1KCBcn1Ly2spJRUVFaxbt45169YRHh5OZmYmBw8eJDc3l7lz52I2m0lISDjNuEFt8Pv9fPnllyQlJbFx40buuusuYmNjGyk1FxZut5uHHnqIe+65h8GDBzea4YkWR3Z2Nt27d2fixIkIIfD5fLz66qtUVFTU2y6cTqfDZrOpR+FaOkpjEPrwFeVWtaGoWnmrnviozuBA1bDQ3k1Vk1DK9YriU74LlSPURFVLVIhVsVqtTJgwgfHjx+P3+1m1ahV//vOfueGGGwgPD6/XQpIQgptvvhkhBEOHDtV6jQ2IXq/HYrHUawTYKpSjMsTZv3+/aqV4165dpKamEhUVxY4dOzCZTA2yCVyx1lNRUdEAkjcudbViU9d7SSlJSUkhJSXlrL+prWwtWTmGbpZXGoIJEybQo0cPoqOj+emnn3juuef4xz/+wcCBA4mKijrNZN7Z0qc0HMrvQk3xeTwetbFpyfnUEhEi4OqkPsPqVqEcIdBN/vbbb9m5c6c6XHvzzTexWq04nU769+/fYKca6tsjOF+5UCto1XQnJCQgpeTEiRMMGTKEuXPn8uKLL9KrVy/+/Oc/I6UkJyeH7t27q79XdkFUd1+l8f/hhx/o378/CQkJeL1ePvnkE8aPH09KSsoZF3wu1OdyJoQQhIWF1WsE2GqUo8Vi4bHHHjttaV6ZL1Mm1BsCu93eaobVGk1LqCIaM2YMU6dOxW6389lnn3H48GFmz55NamoqXbp04c477+T9999n7Nix3HjjjWe4a4Dw8HB19KOsnrdp0wYpJatWrSI5OZnk5GTy8/NZt24dV1xxhebr6AzUtx63CuWoDDOUFSjFQsyJEyfIysri559/pmPHjtx///0NYjnFarVSXFzcAJJrnM/07dsXQLWSPXz4cFauXElSUhKlpaV8++23TJw4keTkZNUadk29PJ1Ox5gxY9TPiok9QPXHkpeXR3JyMi6Xi+XLlzNq1KgGOTJ7vqKMKutKq1COCocPH2bBggVkZGRw9OhRsrKySElJYfjw4YwdO7bB/L6YzeZWsc9Ro/kIVXJ6vZ67774bnU7H5MmTMZlM/PTTT2RmZrJgwQLS0tJ4/PHHueWWW6p1/VqdwqwaFhsby969e7nooouIjo7miiuu0BZwzoLFYqGoqKjOv29VyjEvL485c+ZgMBh47LHH+PTTT/nd737HlClTGjQeu93eKvY5arQclNNZisK68sormTRpEt9++y0rV65k9+7d5Ofn4/f7Wb16NZdffrnqWfFsc4ZCCIYNG8bgwYOBQI/oQvSMeK7Y7XaOHTtW59+3qhMygwcP5osvvmDYsGF8/vnnFBYWkpmZqX7fUBPTdrudiooK7dSDRq2oeoZbeen1eqZNm8bs2bM5cOAAv//971m1ahXffPMNCxcuJCcnB+CsJ2yUk0zKinZtz4yH3vtCO8UjhMBut18Yq9VKIejVqxevvPIKP/30E6+//jrz5s0jMTGRyy+/vJKtwvqgzFUoq+IaGnVBOXE1YcIEBg4cyIkTJ+jatSsGg4HU1FT++9//csstt9ChQwd69+5do2+dkydP8uWXXzJr1qw6HZHNyMjA6XQyaNCgBklXa8Fms6mdnLrohVajHOHXhRmr1co111zDgAEDWLJkCfPmzePo0aM89thjDbJibbVacbvdmlVojXqjlNm4uDjVIPHUqVOZPHky33zzDYcPH+bIkSP06NGDtLQ0SktLGThwoOpKFgKeN++44w6sVqtqaUbxgVQbnE5nq9i329BYrVZcLledOzmtSjlC5b1hCQkJzJo1izvuuAO3233WDKhuWFFdAVMMZV5IwxCNxqNqGVNOYd11111A4Bihy+Xiww8/pKioiD//+c+VzmobjUbatWsHQFZWFitWrODee++ttVvb0B6j3+/H6XRWcnFxvhJaj+vSe2y1Y8bQOZeIiAji4uLOmHjlHO+OHTtYtmwZx48fr/Fa5ayyphw1GoOq84U+n4/XX38di8XCzJkzOXjwIJ9++in79u3j1KlT6s4JIQTx8fFMmTIFg8GA2+0+65xa1XnJnJwc/vWvf10QC4719SNzVuUohJgthMgVQuwKCfuHECJTCLFDCPGdECIq5LunhRBZQog9QohJIeGTg2FZQoin6iRtPdm8eTOzZ88mPT2dN998E7fbXe2EtZKp2rBaoykwGAxcf/31PPLIIxQUFPDll1/yyiuv8Oijj/KHP/yBvXv3AqhTSp06dUJKyerVq/nmm29qHU9ZWRnbtm1j6tSpWK3WxkpOi0FxeVLXelybYfUc4G3gs5CwFcDTUkqvEOIlAs6ynhRC9AZuBvoAHYGfhBDdg795B7icgCfCTUKIRVLK9DpJXUcOHTpEXFwcAwYM4LvvvsPr9WIymdizZw979uzB6XRSVlZGx44dtZ6jRpOg9Op69OgBQOfOnbnqqqs4fPgwGzZsYP369axZs4b09HTGjh2LzWYjLCyM0tJS1q1bxx133AHUzue64lGyY8eOZ5ybr423ytZAqBWqulAbHzJrhBBJVcJ+DPm4Abgh+H4K8JWU0gUcEEJkARcFv8uSUu4PCv1V8NomVY5JSUn861//YsOGDfTt21ctIPPmzePvf/87ENgA/qc//akpxdK4wKliTku1yh4bG0tycjKbNm3ijTfeoHPnznz33Xfcc889xMfHM2vWLKxWK7t27aq02l2TMouMjOSmm26qlUwVFRV4vV7Cw8Prn8BmQlGOdVWQDbEgMx34d/B9JwLKUuFoMAzgSJXw4dXdTAhxL3AvBA74NyQ///wz9957L2PGjOFPf/oTBw4coGfPnlx99dUkJCRgsViIiIggIiKCtLS0Bo1bQ+NcEEIQHR3NJZdcQnh4OGvXriUjI4OdO3fy3HPPMWvWLNq2bcuJEyfYsGEDKSkp+P1+LBbLaW55lftVpbrTOgobN27k6NGj3HbbbY2TwFZAvZSjEOIZwAt8oQRVc5mk+rnNatW5lPID4AMIGLutj3xV0ev15ObmcuDAAeDX0wwDBgxgwIAB6nXKHE9rHU5onB8o5S8+Pp4nnniCiIgIevfuzeLFi2nfvj1/+MMfsFqtPPvssxw5coSFCxdy22230aFDh1rd3+v1kpeXR9u2bdWVb0VhDh8+nCFDhjROwpoIv99fL3NvdV6tFkLcBVwN3CZ/bYKOAp1DLosHjp8hvEm55ZZbCA8PZ/PmzUyfPp3ExETg9NVDr9dbbeurodEcxMTEMHz4cHr16kWvXr3o3r07bdu25aGHHmLcuHE89dRTbNy4kVtvvVXdS1mb4WRRURGff/75aSvXXq9XdU2r3KsqreHUTail+7pQp56jEGIy8CRwqZQydC/BIuBLIcRrBBZkugEbCfQouwkhkoFjBBZtbq2TxHVECEHbtm1P8zVTnQJUDOpqp2M0WgKhZVTZmmI0Ghk7diwZGRnceuutJCUlMXfuXDp37szAgQPJy8tj5MiRleYhqw6jY2Njuf/++yutXBcVFfF///d/3HnnnURGRla70KPcJ/RIY0ukvvW4Nlt55gLrgR5CiKNCiBkEVq/DgRVCiO1CiPcApJS7ga8JLLQsA+6XUvqklF7gAWA5kAF8Hby2Sanu/Gt1KBaYNeWo0dKIiYnht7/9LTabDZ/PR5cuXbj44otZuHAhhw4dYvbs2XzwwQccOHCAJUuWsGfPHuBXfzih+3t1Oh12u71SObfb7UyYMAGbzUZpaSmzZ8+mpKTktB5iUVER77zzTr2s3jQ29bWkftbaL6W8RUrZQUpplFLGSyk/llKmSCk7SykHBl+/D7n+BSllVyllDynlDyHhS6WU3YPfvXDOkjYhikMobVit0dIQQqhl0+FwkJWVRVhYGDfeeCN2ux2fz8fgwYOZPHky27dvV5VjQUEBn376Ka+//jper7fSPd1uNz/99BOFhYUYjUZ69+6NyWTCbDbTv3//anuHYWFhTJs2rUafTYp18+YcdivD6rrS6o4PNgUul0v1kKeh0VKJiYlh1qxZGI1GPB4PDz30EEVFRSxYsACdTsejjz6q9jCzsrJITU0lKSkJv99PaWkpYWFh6nC7pKQEn88H/Lpn0mw2M2zYMDU+JVwxqNGlS5dq5VLmIr/++mv69+9P7969myQ/quJyuTCZTI03rL4QKS0txW63a8pRo0Ujgr5pANauXUtWVhY9evRg5syZFBcX8+GHH5KVlcWcOXMYMGAA06ZNY+TIkeTm5vLyyy+r/phMJhOXXXYZ//nPf6o9VpiZmckPP/xwWvjZ6Nu3L23atKl3OuuClFJtAOqK1nOsBofDoSlHjVbF1KlTVQ+JMTExhIWFcdttt2GxWBg5cqTqXuTQoUP07NmTiRMnsnHjRlwuFxUVFQwePJjo6Ohqh6F2u11dBa8NSu9ScSPRXDgcDk05NjRlZWXY7fbmFkNDo9ZUdbRlMplUSz6K+bNrrrmGzMxMFi1axL333svIkSPJzs7GaDRitVoZP3488OtqtNI56Ny5M507d6a1oYwA64qmHKuhsLCw1htpNTQaC8WSlNKbq2kkczZrVAoGg4FevXrRqVMn9uzZo+6ZBFixYgUQ6G21bduWIUOGVLI8XhfZaytjY1FYWKieWa8L2pxjNRQVFan+PTQ0mosTJ07wwQcf1MvZm8/nY/78+Rw7dkx13SCl5Oeff8blcqlKKykpSd3XWFxczKuvvkp5eTmnTp2isLCw2u08CsoCTNXVabfbzaZNmyp5AAy9rjFXspVFpvrUY005VkHZD6YNqzWam6ioKMaOHVsv39RCCNq0aVNpO05MTAy/+93v8Hq9qmXxZcuWkZeXh8/nY9++ffTr148ff/yRt956i0cffZT3338fj8dTYzxut5vPP/+cEydOVArbsWNHJeVYXFzMl19+qS78OBwOPB5PgytKxahvXdxKKGjD6ip4vV7Ky8tbtTUSjfMDq9VKr169gLoPS3U6HaNHj1Y/K1t3SktL+fDDD5k5cyaRkZFMmzaN2NhYLBaL2qNcs2YNAwcOpHPnzlxzzTUYDAa1R1ZSUkJ8fLx6X71eT48ePSp1KsLCwpg+fXol2Y1GIwkJCaqV7k8//ZRRo0bRr1+/OqWvJpR6XJ8FGa3nGIKUEo/Hg9PpJDw8XFut1mhWauthsC73iIyMZMaMGURERKgWxhXXCRaLBbPZzOWXX05ZWRl6vZ6MjAw2bdoEwK5du3jxxRdxuVz4fD7S09PJz88nNjZW7VQocSp7DJWeoc1mY/To0ZjNZvR6PTfccAMpKSkADTrcdrvduN3uenVyNOVYBY/Hg8/nU/0Qa2icbyhzj3Fxceppm+peUkoKCwspLy+nbdu2xMTEANC/f39mzpyJwWDg6NGj/O1vfyMjI4NVq1ZVsrqtKLrjx4+r7hyqKus2bdqojsPmzZvHvn37GiSNHo8Hv99fr7PfmnKsgjL/0ZIP1Gto1Jfa9kp79+7NpEmTaN++PTt37sTtdhMWFka/fv3Q6/XExsbyxBNPcPHFF3PDDTec1uvz+XwsWrSIQ4cOnVEGIQQpKSmVvC6Gcq7+t5X5UWWTfF3QlGMVlEwN9eymoXEh4ff7KSsrA2DChAmkpKSg1+vVgxGhSjUsLIz+/fsjpeSzzz5TbaUq6PV67r777rNuqRFCMGjQINq2bVvjNbt27WLt2rW1SoPb7UYIUa96rCnHKng8HvXsqIbGhUhZWRn/+te/yM/PBwKKKzIykokTJ1ZbLw4ePMjcuXP5zW9+o9pIVX6nzGGezQBEdUP6qj1Fo9HItm3bOHny5Fl7j0o91pRjA+LxeNDpdPWy5qGh0Zqx2+3ccccd6h5Br9erbsepalvy2LFjHDt2jEOHDhETE6MOY6WUOJ3Oaj181naInJuby48//qgaxOjevTvJycmVtgbVhNvtRq/X16sea8qxCg2Rqa2RGguv8k8pyOrrDL85hwrQmjmX9DfcK5D9jYlOp6N9+/Zqr2vv3r3Mnj1bVVIKFRUV/OUvfyE3N5f8/Hzeeecdjh49qj7vxYsX89///pc9e/awdetWdbFGSkleXh5ut/uMcii7R+DXnuU111xDUlLSWdOgGLrVTJY1IIqLhAvT0K3E5wedALWDICHUNZCs9E6q7/x+8HoClddkEpyefefrtiglRwLpk9IPPj/S4wG9DmE0BfOyJvdKotJf6ZdIrxfp86IzmhAG/a9fq9fWdL+GIbR3KKUkOTmZmJiY0xSNxWJh0qRJREZGcskll9ClSxdsNhtHjhyhoKCAUaNGUVpayssvv4zD4eC+++5j2LBh6HQ6vvjiC6699lq6dOly2lluhXbt2nHVVVdVK9fZUGw51mcblKYcq+D3+y9YQ7fr17v57jsP7dsLbrrJTHiEYNNGF+PG2dDrwe+X7N7toVs3ExYLOCsk6RluIiJ0fPyxh9JSN36/ZPRoEwmJelI3uOnR3cToS81EhFdfAVozzt27KVq2DHNSEuGXXoo+NoaSJUsoW/kfvG4XOouV2OnTcWXtQ0gIGzMaQ2wsIHAdPkLJ90sIGzMGS/duYDTir6gg/933cB08iM/jwRTfmejrplL84wqMMdGEX3opps4JIS1X4yOEwGq1VnKloKDX65k2bRpz584lOTmZsLAwDhw4wMKFC4mJiWH69OnExsbyP//zP/h8PpYvX47b7eayyy5j+vTpmM1m1q1bR3JyMh07djwt3vrQEAara+MmYbYQIlcIsaua7/4ohJBCiLjgZyGEeFMIkSWE2CGEGBxy7V1CiH3B1111lriRUTyWXYhs3OgmPh5sVsG77zrJOSlZvNiLYjja5YKPP/JxKjBPT16+n09mV3D8uCQnx8+sWWE88kg4wy8y8/NaN3m5ftZv8PDSi2WUlPhrjriV4tq6Fc/+A7gyMjnx/F/x5OTi2rUba9cU2j32B+KmT8fYpg0lPyyjeMFCjv/xcUqWLUd6vbgOZFO8aBF5r7xCzt9fxJuTg/R4cKftIHrqVNr/4TFipk7Bm5OLc8MGfPn5nHz+bF7iYwAAIABJREFUr1SkN6mr9zOiDHVvvPFGhg4dyubNmzGZTDzxxBPMmjWLn376iePHj5OUlKQOvZOTkwHYvHmzan28PmfHa6Km3ui5UJux4xxgctVAIURn4HLgcEjwFQScanUj4Hv63eC1McBzBHxVXwQ8J4SIrrPUtaSmeZsz7cTX6XSNeyAeKs3hVZ7La/iYfp0jrGYuMWTeEMDt8tGli47evfWUlRFUihIRFMzvl1Q4BTKo53xeidfrx++X2Gw+/H4oLPBjsYHJKBg1ysT/+392Skt1pO/2Vitla0TNM50O+5DBtH3icfRREZSlpiIliHA7/rIy/BXlCKsVwsKIe+Qh4h64n6Jvv6Vs3c9IocM8fDgdXn4JfXQMeW+8iSwrRxr0YDLhLSoK+DQ2mbB07UrcAw8QNmoUJcuWVZLh13nIyp/VZx76r7q5y+rmjvm1TKrXniE/jEYj+/fvJzs7m5KSEsLCwjAYDHTo0AGbzYaUkp49e9K5c2fSg8q9d+/eDBw4kKuvvlpVmA2JohTrU5fPOqyWUq4RQiRV89XrwBPAwpCwKcBnMiDRBiFElBCiAzAWWCGlLAgKvoKAwp1bZ8lryZEjR1i5ciU33ngjAAsWLODAgQOMHTuWiy+++LSut16vb1zfF+p9JRIoLPBhseqwWXUEFFFD9Frl2Sftpfqfen1BgZ4P3oeCAjdTpwmk9KPTeRHBJlQI0OsD82IADofEZNLhcvrZkSb44P1y7GF+rrveQkGB5JdffOza7aS83EeHjudfb9xXVERFRga+4iI8h44Qdd11ONasxbllEyJtB8ZOnYhJSEQiwGjE1rcvMbfeSsnyH7FPmggCDHFxxM64m+NPPoUzIwPfsWMUzvkUGRWJbcgQjFFRuLKyyf/gQ8o3byHqxhuqSKHM/Ypfy5asofwKgkNyoU5ZVu88vup36qRntezevZtNmzbRrVs3OnXqFBjSAiOGD0foAr+Liori9ttvV+VqbJOAer0en8/XuMqxOoQQ1wLHpJRpVbqtnYAjIZ+PBsNqCq/u3vcS6HWSkJBQF/FU3G43c+bMYd26dVxzzTUsXbqU7OxsxowZwxdffEGHDh3o2rVrpQxUDsRXXZlrWCQVTsm6tU527HRx041hmNoLhE4gRKAgCl3V4ljFPt4ZQwP1xO+TeH1+3G6oqACHw09pqcThgLJSDyWlEBenY9QoC0KARMf0GTpsNh1zPqmgSxeJyQg6nUBKMBjAbvdw4qSOzokm9u7x0SnegEQw+Qo9f/yDHUNwG9zX//aRmenn5AkvM39voEOH8231X+JzuvBkZVGxcSPh112PObkL0mQk9rE/YL/oItDrkV4vOr8fnT5Y1QyGwMPxeNErYTodQqdD+jzokpLo8Je/oI+JAb2e4qVL8eTmUDL3K8z9+mIfOhQAv6Oc4h+X4y8qwufz4ff58fl9+H0+9D4fumADFrqEg07g1enQGYyYzCZ0VivYbAirFZPdjrDbwWrDEGZHH2ZHmC3ojEbQ64O9UIXKpc9oNNKnTx+GDRtGSUkJc+fO5cpu3fDt2UP4JRdjTkpC6A1NerBCsXweepzxXDlnaYUQNuAZYGJ1X1cTVlOzU32jJeUHwAcAQ4cOrVf3zWg0cvvtt3PkyBE8Hg9btmzh/vvvp1u3buzYsYP09HS6du3KunXr+Pnnn3G73fTr1w+/31+vTD0bpaV+PpntIDpacMkoMwsWunE5vQidDp0And6HweDEYAgYOjUZTZjMOvR6gRABZRXg1+GR3y/xegUejw+Xy43LpcftNuH3eZC4MRjMWK1+wsLAbrcQFqYjPl7Qrp0hqPwk0g8bNnjweiXxnfWEhRnIzfMxb14Fe/f6GDPGwLhxRj77zElGhpfUDR5+f58NR5mPY0e9/LzOhc8n6dXLQEyMl8cftxIeLvjwgwo6d3YyaZKV82fVWqCPjCTqnnuw9OnDqQ8+oHDRIvQ2O+XbtuF3OgO9xQED0LndOH7+GXdWFsXfLyXy+uvQ6XR4Dxyg9KeVOHfvArMFc89e4PmGsvW/IGx29G3aojObCbvqKmJu/g0Fn37GqY9n0/axxxAmI9a+fcHrRRgMCJ0e9Dqk3oAw6BGVVpaVlXA/0ucDtxtcbnwV5XjLy5Hl5fjy8nAdOoi3wonOUY7X7UKPRG8w4LeHYYyIQISFoQ8LwxAejs5iAbMZncVCl5gYRNu2+MtKMUvJsAEDsHXsgNvroXDxEvQmIzG33ooxOqbJnk6zKEegK5AMKL3GeGCrEOIiAj3CUHvq8cDxYPjYKuGr6xD3OaHT6TAajZhMJrWbrZgwMplM6h6qtWvX8vLLL2M2m3nyySebRDmePOnDHqZHJwQ33GDEbNLh84PfL/D5DLjdBrxeic8ncLsFbrfE5wPpB3/ISEoEXzpdoFNiNBowm01YLAKLRYfJZMZkAoNBoNcrC50BBVV5ywbccquRfXt9xMfr6NPXgNMpSU93kZvrZchgI336GImMFHTooOPwYR8PPmSlZ08jDoeei7J8rFvnxGDwEx9v46bfWImJ1hMVrefPz+twu888NGtdBJRN2ORJCIMBQ/t2dHj2WTzFxYQNHUrh999Tum4d+vbtsA0YQMyM6ZT89794MvcQfesthI0cic/hwH3yJGWbNmPq2J62N9+CoU0bom65mbKNG5E+H9ZBgwgfPhxT1y4YO3Wi7aOP4jp2LLBFSK/H0q17gyxcKyMniUTIwHYk6fYgXS78zgp8JSX4ikvwOhz4SkvxnTiBx+nC63Khd7vA40H6/fgAnwQrkB8VSeToUZg7d6Z8+3akr2kX5AwGQ73r8TkrRynlTkA9ACmEOAgMlVKeEkIsAh4QQnxFYPGlWEp5QgixHPhbyCLMRODpOkt9Dij7nZTtCDk5OcTFxXHkyBEGDhwIwO23387YsWNVX71btmw5zbdvQ9Kxo5Fn/ieSzEw3W7e68Xr9jBxpo7LyaPph6MCBJoJZAkBEBDzyyOmWlAcNMjNo0K+fw8P13HWXHSlF6HSWSpu483HHmMAUas8wKgp98ERJB+UccbAhMgwejG3woGCDFsgdQ1QUsbfeGmiVQjRcxIQJRFx2WfD3gXBj0IOfsFqxBs17NU6Kgs9P6BEWPVgs6CMjMbZrX/0PlEUbvx98vkDP1O8Hvx/XwYMULV6CtVs32j/2GIboRl9/rYTJZMLv99erHp+11Aoh5hLo9cUJIY4Cz0kpP67h8qXAlUAWUA7cDSClLBBCPA9sCl73F2VxpjFxu92sXLmS9PR00tLSuPTSS/noo4/o2LEjDodD9aerOBCSUqrnSRtjewH82luz23UMGWJh0CAzfr887fumpv7xiqbcftes/JrOGhJcbUZU02rUdG0TZ2Sdn33wd0Kvh5Az11JKTPHxtPntbzHGxSGa4UCF2WxGSnnWUzhnojar1bec5fukkPcSuL+G62YDs89Rvnqh1+sZOXIkvXr1IikpiTZt2tC2bVtOnTrFoEGDiIqKOq1gKE7AG0s5VkWnEyFziBoa5wdN3VOsislkQghRr3p8Po53VPR6PX369KkUdtFFF9V4vWLFQ6/X1+pwe324UDeaa5z/tISybTQa0el09eo5XogHiM+IXq/HYDDUK1M1NDSaF8V4jNZzbED0ej1RUVE1ekQ702mbUHt0yubymgxYNFXrWtMm2KqyV7WgUzUtZzqn2txpUVYla0pH6LOoyfp1Y6ehprJUVW5ldbW6vG+u/K9NGaru1FlzliGDwUBERETjLshcaOj1eq644gqOHDnCkSNHcDgclJeX43A4KCsro6ysDKfTic/nw+v1qu4tQ5WjXq/HZDJhNBoxm81YrVbsdjthYWFER0cTFRWlfg4LC8NisWAymdQtR0pBqq5SVK1IPp8Pj8eDx+PB7Xarf51OJ+Xl5VRUVFBaWkppaSnFxcVqOsrLy3G73aeloapSUXrSFouF8PBw7HY7ERER2Gw2rFYrYWFh2Gw2zGazGqak3WAwqK+q6akuXYqSU1YZlfQ4nU7KysooKSmhsLCQwsJC9XNJSQkVFRW43W51837oCSfleRgMBjWPzWazmp7Y2FhiY2OJiopS06LIH5qGqgpWQZFXKQ+h+e90OqmoqMDhcFBcXKzK7HA41OdQUVGh+i2q2tAq+W82mzEajarMYWFh2O129VlERkYSGRmpliVF7jOVo9C4lHzz+Xy43W5cLhfl5eVqGVLkLy0txeFwUFpaSkVFBRUVFbhcLjweT6UypFC1DJnNZsxms1qOlDpgtVqx2WzY7XasVqtaH5QdJqH1IvRYYNU6oJQZpcxffPHFxMbG1lkXaMqxGoqKitiwYYP6wGw2G1FRUcTHxxMeHq4WwFAloFgvViqJYuhTKUQOh4P8/HwOHDigKrPQSg1UamGVQlXV7FJoL0lRIqG/D1UGihc5RQlHRkYSHR2tVi4lHU6nk+3bt5OYmEi3bt0A1LQosioVpLy8nJKSEk6ePKmmS0mry+U6rZKEyhSatuqUjPLbqgpOp9OpyjciIkJNQ1JSEhEREWpaFEWmxBV6X5fLpcqoKK2SkhKOHz/Orl271AqvOGZSUJ5D6N9Qj3qhijH0uJpOp1MVcqhSi4iIIC4ujsTERFVJmEymSrIrciv573Q6cblclRq6oqIijh07pn4uLy+vZPswtAyF5nmo3KGKvapiNhgMlRr3iIgIVV7FKZbSGCrKW6kHoc9UyZucnBz27t1L27Zt0ev1lJWV4XA4yMvLo6KiQm3IlQZRkUmRpepzPVveK0p4UOies3NENKaRhfoydOhQuXnz5iaNs2rvqSFZs2YNzz//PI8//jiXXXaZ2tIpf51Op1owQr8LHWqFKs3Q1jhUUSvfVW1tayI9PZ2rr76aW265hb/+9a/nlO7QFlzphSo9qKq9WqUQK6/QyqjX6yv11pStGAsXLkSv13PnnXdis9kaxQhxaKOm9EBCe66hvSPluuzsbIqLixkwYADh4eHqcwh9FkajUX1WjWUfNLQRc7lcbNu2jS1btnDJJZcQHR1dqVen5HnVMqTIqihDpeE/lzJ0Nn766SfuueceXnjhBW699dZq06HIqOSx0iAovdhQJVi1Dii9TUX2UPmhcl0WQmyRUg49m8xaz7EKZ5rbqQ9SSrxeL6mpqWRnZzNx4kR1yNSY1CYtFosFv99PXl5epemBc7m/0lNqSMrKypg/fz5+v5877rij3sZLa0JJr9LbPhtSSv7973+zfPlyFi1aRNeuXc96/8YiVHar1Up2djYvvPAC33zzDSNHjmywOOqL2WymsLCQnJycau8Z2vAr1MfndHVxnCuacmxCzGYzOp2uWn8czYnSa3A6nXW2YtIYaVEaj6KgcYXGVjK1RcmjoqIidVdDcz7L0Hk4i8WiDsWbW65QzGZzrbbItRR5QVOOTUqnTp146KGHGDJkSHOLUomIiAjuueceEhISWlbhNBi48cYbqaioqFWPrqkQQjB+/HhiYmJUR/cthT59+vDII4/QpUuX5halEh06dODBBx9k+PDhzS1KrdHmHBuRc8nbplRK5/rMW6JsTa3Ez4dn2VLlgqaVTZtzbGakDHhO27x5Mzk5OQwbNoxOnX41YXny5ElSU1OJi4tj2LBh6nGnppDL7/eTnp7Onj176Nu3L927d0en0+Hz+diyZQunTp3C7/fTqVMn1ThHU6DItmPHDtq1a6caRJVSkp2dzY4dO+jatSt9+/ZttPnHmuQCyMrKwufz0aNHD4QQOJ1OUlNTcTgc+P1+evbsedb5x4aWy+12s3PnTkpLSxk4cKDqTlVKyf79+0lLSyM5OZn+/fs3WZ4pz3Hv3r0cOXKEvn370qFDBzXPNm7cSFlZGX6/nx49epDSiMY06oN2QqYRWbJkCV999RUHDx7k1VdfpbCwEICCggJeeuklsrOz+fbbb1myZEmTyKNU8u3bt/PWW29x8uRJ/vnPf5KVlQUELBh9+OGHrF+/nqysLEpKShrPInoNrF+/ngceeIDdu3erYdnZ2bzyyiucOHGC999/n61btzapTAD79+/nkUceYfny5WpYSUkJb775Jmlpaezfv7/Rj5xWRUrJggUL+Oqrr1i6dClvv/22uq3rwIED/OMf/+D48eN8+OGHNPUIbOPGjbz77rusX7+e559/npKSEiCwyPbWW2+peVZRUdGkcp0LmnJsJCoqKli9ejWzZs3iwQcfxGazkZGRAcC2bduIi4vj4YcfZvr06axdu7bJjitKKVm+fDnXX389999/PxdddBHr168HAh7bnE4nxcXFHD9+XO3NNqWCTElJYcSIEZXCVq1axahRo5g1axZXXXUVq1evbjJ5FGJjY9UdBgrKFqX8/Hxyc3MbfedBVYQQXHzxxTz77LNMmDCBwsJC9VmtXr2akSNHcv/993PttdeyatWqJpWta9euPPvss9x0001UVFSoJ1WULUenTp1qljw7FzTl2Egop2iUTa8RERGUlpYCkJubS3x8vBoeuhG8sfH5fBQWFtK5c8AmcUxMDMXFxep3NpuN0aNHM3HiRL788kvKy8ubRC6Ftm3bYrVaK23jyMvLIzExESEE0dHRlJaWNnmPNjIyUj2NogxNfT4f0dHRXHnllfTu3ZuvvvqqUY0kV0fbtm1ZtWoV33zzDTfffLO6nSo0z2JiYpo8z2JiYti3bx9vvPEG1157rTrc9/l8REVFccUVV9C3b1/mzp3b5HlWW7Q5x0ZC6XU5HA4iIiIoKCggOmjGKSIign379gGBIbay8bYp0Ol02Gw2VSHm5OQQGxuLlBK73c4LL7xAdHQ0+fn5zJs3D5fLhc1maxLZ4Nd5NLvdroaFh4erUxK5ubnVmpprCpxOJ+3bt6/kJOqf//wn0dHR7Nq1iy1btuD3+xtlo3pNLFq0iNTUVJ555hkSExPV8PDwcAoKAiZTc3JymjzPtm3bxpw5c5g5cyYDBgxQvXq2a9eO119/nejoaDIzM0lNTcXn8zVpntUWTTk2EjabjW7dujF37lwSExMpKCggMTGR/fv306NHDxYtWsT333/PmjVrGDJkSJM5H9LpdAwbNoz58+dTUFDAhg0b+OMf/8jBgwcJCwtj4cKF9OnTh61bt9KmTRvVrURTcejQIQ4dOqQeZ8zPz2fIkCF8+eWXmM1mfvjhB+68884mlQkgPz+fffv2kZuby7hx48jPz8dsNrN06VIGDRrE4sWL6dWrV5M6kfL5fKxbt47Y2FhWrVpF165dGTJkCLm5uQwZMoTPPvsMu93O999/X+2plMZk06ZNmM1m0tLSOHbsGOPGjSMnJweLxcKSJUsYPHgw33//PT169GiyjsG5og2rGwkhBHfccQexsbEcOXKEhx9+GKvVyrJly4iLi2PGjBns3r2bQYMGMXXq1CaTCeCyyy5jxIgR7Ny5k7vvvpuUlBRWrlxJUVERCQkJLFmyBI/Hw3333YfRaGzSHkdFRQXjx48nMjKSvLw8VqxYwcCBA7nmmmvYtm0b06ZNO6NNzsbC5XLRt29fUlJScDgcLF++HIPBQHh4OPPnzycpKUl1/9tU6HQ6rrvuOhISEggLC8Pn85Gfn8+PP/5Iv379mDJlCtu2bePaa69tsNMytWXs2LEMGjRIHXU4HA6WLVuGTqcjKiqK+fPnEx8fz80339ykcp0LZ93nKISYDVwN5Eop+4aEPwg8AHiB76WUTwTDnwZmAD7gISnl8mD4ZOANAs5RPpJSvng24c7HfY5er7fa3kVzm/1SfO0IIfD7/ZXOAjeXbMq54ebML0WOqoQ+x5aSX0pYc+fZmcpYc+dZMK4G2+c4B3gb+Czk5uOAKUB/KaVLCNE2GN4buBnoA3QEfhJCdA/+7B3gcgKeCDcJIRZJKdNrn6TWR3UPvLmHEDUVwlC5mmv+p6bzts3N2Z5jS8kvJay586wll7FzoTY+ZNYIIZKqBN8HvCildAWvyQ2GTwG+CoYfEEJkAcoYKEtKuR9ABLwTTgHOa+Wo0Xg01Emamu5T0++q2ivUOH+paxPTHRgthEgVQvxXCDEsGN4JOBJy3dFgWE3hpyGEuFcIsVkIsTkvL6+O4mlcCPh8PjIzM1WjsVlZWezbtw+Hw3FO91H2Kjocjl99ONegNIuLi9m3b99phl01zj/qqhwNQDQwAngc+FoEmtHqmtKavLlXW7KklB9IKYdKKYe2Cfrr1dAIRVFKP/74I8899xxFRUXMnz+fO+64g1mzZvHwww9z8OBB9dqzKbH58+dz3XXXcfvtt/PGG29w6tSpSr8N/X1OTg5PPfWUunFeU5DnL3VVjkeB+TLARsAPxAXDO4dcFw8cP0O4hkadOHXqFB9//DEzZsygU6dOHDlyhMsvv5yPPvqIyMhI3n//ffx+P263+6wK8vjx43Tr1o2ZM2eSlpbG008/TWFhoWpENpSUlBR+85vf8N5776mb+jXOT+q6KWsBMB5YHVxwMQGngEXAl0KI1wgsyHQDNhLoOXYTQiQDxwgs2jTtxiuN84pNmzZhtVoZPXq0atAgLi6OhIQE+vXrxy+//EJ5eTmvvfYaiYmJ5OTkMGTIEMaOHXuaAQaz2UybNm2YPHkyQ4cOZfr06aSmpiKlJC0tjYiICIxGI9ddd516jPDf//43u3btavItMhpNx1l7jkKIucB6oIcQ4qgQYgYwG+gihNgFfAXcFexF7ga+JrDQsgy4X0rpk1J6CWz7WQ5kAF8Hr9XQqBObN2+mb9++WK1WIGAEYtu2bbzyyiu8//77TJo0Cb/fz6pVq/j8888pLi7mb3/7G5mZmUBlz3mFhYVYLBZOnTpFeno6paWlWCwWDh48yIcffsj+/ftZvnw5n3zyCVJKIiMjSU5OJi0trTmzQKORqc1q9S01fHV7Dde/ALxQTfhSYOk5SaehUQ1SSvLy8lSjwX6/n7KyMo4dO6Z6Vdy1axdjxowhJiaGmTNnMn78eB544AF27drFoUOHSE1Nxe/3c/3115Ofn8+SJUtYvHgxADNmzGDEiBEcPXqU8ePH8/e//501a9bw8ccf43a7MZvNtG/fXjX5r3F+oh0f1Gh1CCEwmUzqfKDin2fWrFlcccUVbNq0iaeeeoorr7xSdXil4Pf7SUlJUb3+KQY4Lr30UgB27NihnupQzryHegRU5i4VJalx/qIpR41WSWJiItnZ2aqyUobQ6enpbNiwgcGDB5OUlITH42Hu3Ln88ssvZGdn8/DDD9O9e3e6d++u3ks5Bz99+nRWr17Nq6++itPpJCYmhq1bt/Lmm2+ydu1aRo8ejdlsxuv1cujQIVWhapyfaMpRo1UyfPhw/vd//5ejR48SHx/PlVdeyerVq9Uz4ZdccgkWiwUhhOri87nnnqNHjx5AZadUkyZNwmazYbFYmDx5MklJSeTm5pKbm4ter+fkyZNcddVV3HDDDQghyMrK4tChQwwePLg5s0CjkdF8yGi0OqSUuFwunnnmGUpLS3nppZeIiopS3coquN1ufvvb3zJz5szTenmhyrEm5s2bx3/+8x/efvtt9Uhebm4ujzzyCH369OHJJ588zZG9RsuntmerW7RyFELkAQ4C24Sai7gLPP6WIENzx98SZGju+FuCDA0Vf6KU8qwnTFq0cgQQQmyujZbX4j9/ZWju+FuCDM0df0uQoanjb36TJxoaGhotEE05amhoaFRDa1COH2jxNzvNLUNzxw/NL0Nzxw/NL0OTxt/i5xw1NDQ0moPW0HPU0NDQaHI05aihoaFRDS1WOQohJgsh9gghsoQQTzVSHJ2FEKuEEBlCiN1CiIeD4f8rhDgmhNgefF0Z8pungzLtEUJMaiA5Dgohdgbj2hwMixFCrBBC7Av+jQ6GCyHEm0EZdggh6nVMQwjRIySd24UQJUKIRxo7D4QQs4UQuUHLTkrYOadZCHFX8Pp9Qoi76hn/P4QQmcE4vhNCRAXDk4QQFSF58V7Ib4YEn11WUMZa7wivQYZzzve61pUa4v93SNwHhRDbGysPzlD/mqwcnJFQa8ct5UXAQ2E20IWArcg0oHcjxNMBGBx8Hw7sBXoD/wv8sZrrewdlMQPJQRn1DSDHQSCuStjLwFPB908BLwXfXwn8QMBG5gggtYHz/SSQ2Nh5AIwBBgO76ppmIAbYH/wbHXwfXY/4JwKG4PuXQuJPCr2uyn02AiODsv0AXFHPPDinfK9PXaku/irfvwo821h5cIb612Tl4EyvltpzvIigQy4ppZuAzcgpDR2JlPKElHJr8H0pAVuT1fq2CaI6EJNSHgBCHYg1NFOAT4PvPwWmhoR/JgNsAKKEEB0aKM7LgGwp5aGzyFXvPJBSrgEKqrn3uaR5ErBCSlkgpSwEVgCT6xq/lPJHGbA9CrCBgMX6GgnKECGlXC8DtfSzEJnrJMMZqCnf61xXzhR/sPd3EzD3TPeoTx6cof41WTk4Ey1VOdbaIVdDIQIeFgcBqcGgB4Jd99lKt74R5ZLAj0KILUKIe4Nh7aSUJyBQiIC2jSwDBCy0h1aGpswDOPc0N6Ys0wn0UhSShRDbRMCh3OgQuY42Qvznku+NlQejgRwp5b6QsEbLgyr1r0WUg5aqHGvtkKtBIhMiDPgWeERKWQK8C3QFBgInCAwvGlOuS6SUg4ErgPuFEGPOJG5jyCCEMAHXAt8Eg5o6D84oXg1xNlZePAN4gS+CQSeABCnlIOAxAq5AIhop/nPN98Z6HrdQuaFstDyopv7VeGkNcTVKHrRU5dhkDrmEEEYCD+YLKeV8AClljgy4d/ADH/LrsLFR5JJSHg/+zQW+C8aXowyXg38V3+CNlTdXAFullDlBWZo0D4Kca5obXJbgZP7VwG3BYSLBoWx+8P0WAnN83YPxhw696x1/HfK9MfLAAFwH/DtErkbJg+rqHy2gHEDLVY6bCDrkCvZobibgvKtBCc6rfAxkSClfCwkPncObBiireYuAm4UQZhFwFqY4EKuPDHYhRLjynsBbxqFNAAABYElEQVSiwK5gXMqq213AwhAZ7gyu3I0AipUhSD2p1FNoyjwI4VzTvByYKISIDg4/JwbD6oQQYjLwJHCtlLI8JLyNEEIffN+FQJr3B2UoFUKMCJalO0NkrqsM55rvjVFXJgCZUkp1uNwYeVBT/aOZy4FKfVd0GutFYGVqL4EW6plGimMUge73DmB78HUl8H/AzmD4IqBDyG+eCcq0h3NYmTyDDF0IrDCmAbuVtAKxwEpgX/BvTDBcAO8EZdgJDG0AGWxAPhAZEtaoeUBAEZ8APARa/hl1STOBucGs4OvuesafRWDuSikL7wWvvT74bNKArcA1IfcZSkCBZQNvEzx1Vg8Zzjnf61pXqos/GD4H+H2Vaxs8D6i5/jVZOTjTSzs+qKGhoVENLXVYraGhodGsaMpRQ0NDoxo05aihoaFRDZpy1NDQ0KgGTTlqaGhoVIOmHDU0NDSqQVOOGhoaGtXw/wGjRvrzYpPWIAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import matplotlib.image as mpimg\n", "%matplotlib inline\n", "\n", "img=mpimg.imread('./cmd_M7.png')\n", "imgplot = plt.imshow(img,interpolation ='bilinear')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The main sequence is clearly visible, as well as the \"turn-off\", which is the highest stars that are still on the main sequence for the age of the cluster. The life of a star is shorter the more massive it is, so as a cluster ages, his more massive stars will evolve and leave the main sequence. Thus, there is a maximum mass that one can see for a given age of a cluster. There is one star, however, which seems to be above this - it is surrounded by a blue circle - and is called a Blue Straggler Star. You can read more about it in this paper. \n", "In the CMD of Messier 7, we also see some stars that are above the main sequence. These are districhs (or binaries), i.e. they are composed of to stars - the companion is adding some flux and makes the whole system redder (do you see why?), hence, moves it above the main sequence." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The CMD of the question mark stars\n", "\n", "So, what about the position in the CMD for the planetary nebula and for the stars that form the conundrum? \n", "In the paper, we showed that they seem to be all at the same distance. This was based on the Gaia DR2 dataset, which is less precise that the newest Gaia EDR3. Let's see what we have now.\n", "\n", "As we saw in this notebook, the parallax are not precise enough to be able to directly obtain the distances from them. So, we will use the photo-geometric distances instead.\n", "\n", "We need to get the magnitudes in the three bands and the distance. \n", "We also need to estimate the errors on these. In the Gaia archive, we have the inverse of the relative errors on the flux in the various bands, not on the magnitude, so we need to convert this. This is done with\n", "m = -2.5 log F + c\n", "dm = |2.5 dF/F / ln(10.)|,\n", "where dF/F is exactly the relative error.\n", "\n", "One needs also to take the extinction into account, which is about Av=2.35 (i.e. almost a factor 10 in extinction!) for this region. One can check this using dust maps for example.\n", "\n", "I also compare the resulting data points with PARSEC isochrones (http://stev.oapd.inaf.it/cmd) for two different metallicities and two ages. " ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAETCAYAAAAs4pGmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XlcVPX++PHXZ1gEFGRVwBVBEAXc0HA3l1xyq2w3WyzLtlu3/dd6763rvXrL6ltW3ptd7ZqWpWmZW5nZ4i4qrrkrqCgobqggvH9/zEAiow7KMAO8n4/HeTDnzOec8x4Y5j3n8/mcz8eICEoppdSFLK4OQCmllHvSBKGUUsouTRBKKaXs0gShlFLKLk0QSiml7NIEoZRSyi63SBDGmLHGmC3GmPXGmJnGmEBXx6SUUtWdWyQIYCGQICJJwO/ACy6ORymlqj23SBAiskBEztlWlwH1XRmPUkopN0kQF7gPmOvqIJRSqrrzrKgTGWO+B8LtPPWiiMyylXkROAdMucRxRgIjAWrWrNm2WbNmTohWVTa5mZmcycrCLzwcn5AQV4ejlFtbvXp1loiEXa6ccZexmIwxdwMPAT1FJNeRfZKTk2XVqlXODUy5vf0//8zihx4i5uabaf/aa64ORym3Z4xZLSLJlytXYVcQl2KM6Qs8B3RzNDkoBZB/8iTLXnqJ2k2b0ub5510djlJVilskCOA9oAaw0BgDsExEHnJtSKoy2PDhh5zJyqLbe+/h6ePj6nCUqlLcIkGISIyrY1CVz9mcHLb+7380GTKEkMREV4ejVJXjjr2YlHLInu++ozA/n7i77nJ1KEpVSZogVKW1Z+5cAmNjCdKebEo5hSYIVSnlnzxJ1rp1RHbr5upQlKqyNEGoSilzxQqkoICIDh1cHYpSVZYmCFUpHVy6FA9fX0Jbt3Z1KEpVWZogVKV0ODWVsFat8PD2dnUoSlVZmiBUpVSYl4dXrVquDkOpKk0ThKqUvPz9OZuT4+owlKrSNEGoSskvIoLczExXh6FUlaYJQlVKterV41R6ul5FKOVEmiBUpVSve3eksJDM5ctdHYpSVZYmCFUphSQk4Onnx0FNEEo5jSYIVSlZvLyoe801ZCxeTGFBgavDUW7i1o+WcutHS10dRpWhCUJVWlEDB3I6M5PMZctcHYpSVVL1TBCfXG9dVKVWr3t3vAMC2DlzpqtDUapKqp4JQlUJHjVqEDVoEHsXLuTYzp2uDkdVIloV5RhNEKpSazFyJJ4+PqwZM8bVoShV5WiCUJWaT0gICaNGceDnn8lYssTV4Sgn0G/7rqMJQlV6sXfcQUBUFGv++U8K8vJcHQ6gH2qqatAEoSo9D29v2jz3HCd27yZt/HhXh+N01TX5VNfX7UpulSCMMX2NMVuNMduNMc+7Oh5VeUR26UL00KFs+ve/OfDbb64OR6kqwW0ShDHGA3gf6Ac0B243xjR3bVSqMmn7wgsENGnC8pdfJu/4cVeHo5RbKjxX6HBZt0kQQHtgu4jsFJE8YBow2MUxqUrE08eHlDfe4HRWFr8995zeYa3UeUSEtKlpjG/heDWsOyWIesC+89bTbduUclhoUhLJL7zA/iVL2DxxoqvDUcptbP5qMzPumIFHDQ+H93GnBGHsbJNShYwZaYxZZYxZdfjw4QoIS1U2MbfeSsM+fUgbP57stDRXh6OUW1g3eR3+9fx5MPVBh/dxpwSRDjQ4b70+sP/CQiIyQUSSRSQ5LCyswoJTlYcxhuQXX8S3Th1+HDmSo1u3ujokpVyqIL+A7fO2E39TPBYPxz/23SlBrASaGmOijDHewG3AbBfHpCopn5AQen78MZ5+fvz4wAM6FIeq1iyeFhDwruldtv2cFE+Zicg54FFgPrAZ+EJENro2KlWZ1apfnx7/+Q8Ai+6/n5Pp6S6OSCnXMMbgF+rHqcOnyrSf2yQIABH5TkRiRSRaRN5wdTyq8guIiqLHf/5DwenT/HDffeQePOjqkJRyCf9If05knCjTPm6VIJRyhsDYWK6dMIGzOTn8MGIEp7OyXB2SUhUuqEkQWVvK9t7XBKGqhZDERLp/8AG5mZksuv9+Tu0v1f9BqSqtUfdG5OzKIXtbtsP7aIJQ1Uadtm3p9t575B44wPzbbtMusKpaadq/KQAbpm1weB9PZwWjlDsKT0nhus8+48eRI/npkUfo9v77hCQmujost1CQX0B+bj4FZwsoPFdYYinIL72t8FwhhfmltxWVtceYC253unD1gueNxeC/+iBYDOJhwBh2LNiBxdOCxdOC8TDFjy2eFmpknEA8LBzZfgTjYfDw9sDTx7N4KUsXz6omKCqImH4xrHxvpcP7aIJQ1U7t6Gi6jR/PT488wsLhw2n38stE33ijq8O6KufOnOPEgROc2H+C3KxczuScKbGczTlb/DjvVB75ufnk5+Zz7vS54sdlGaOnIjW4YP1/76detGy07ef/vfar3ectnhY8angQC4iXhXfHrvwjgdT4I5F41fTCu5Z36cXfG99gX/xC/PANsf0M9rV2I60EOj3XiUndJzlcXhOEqpaC4uLoO306vz79NMtffpnsDRto+/zzeHiXrZ94Rck7lUf279lkbcki7JvteGWdZvL0bZzYf4KTB05yJufMRff19vfGJ9DHutT2wS/ED68GXnj5eeHp54mXnxdevtZ1Lz8vPLw9sHhZ8PDyKPHt3OJpweJlKbXNbjlPS+mxES4YF0Hkwg2lYy8sKOSpqakgYAoEU1jI3wYllLxqKfjj8Ztzt2AKhce6x1ivZPILOXf2HOfOlFy+W52ByS+kTVRwqefO5JwhPzefvJN55J3KI+9E3mWTZ43aNahVtxb+kf7417MuAfUCin8G1A/AP9IfY7E3YETFadytMZ4+nnDxt0sJmiBUteUTFMS1H33EurffZvMnn5CVmkqnsWOpHRPjspjyc/M5tPEQh9IOkZmWSdamLLK2ZHFs77HiMqEG8oN8yI8LJSw+jKgeUdSKsH04RfjjF+aHb5AvPoE+1AioUWm+3V7M2Xr+JdYbdLzwmuIPJw7lAJA0LOmSx5xom1fixgc7OBRDQV4BeSfzOHPsDKePnOZ09mlys3M5nX2a00esj08dPMXxjOPs+3UfJ/afoCCv5GCRnr6eBMcEE9I0hOBY68+QuBDCW4WX+Qa2qzFs/jBe6vaSQ2U1QahqzeLpSeunn6ZOu3Yse+kl5t1yC62feYamt91Wur68jGpknOAv5i/2nwuoQf7pfArzL/3N1FgMkcmRNL2+KRFtI4hsG8kTv2xHvDx4w8EPN3X1PLw98A32xTfYl6CooMuWFxFys3I5kXGC4xnHObb3GEe2H+HItiMc2niIrd9sLf7bG4shND6UyOTI4qVuy7p4+Xo55bWEtw53uKwmCKWAet260X/GDJa99BKrXn+d/T//TMrf/oZPSMgVHzM/0AffEF9OZ58u9dzZ42cdOoYUChkrMshYkVG8Lc7TQmEND8a9sQzvmta6ca+aXtSsU5OgJkElltqNalfrhllXMcZQM6wmNcNqEt6q9Ady4blCcvbkkLU5i/2r9rN/1X62z93OuknrrPt7GOok1CmRNMJahF1x0ig8V8jun3azafomNn7u+AAVmiCUsvENC6P7Bx+wdcoU1r75Jt8MGEDSo4/S9NZbsXiW/V+lsKYXz2Y9+8d6QSF5J6zVFAVnC/AL86OGfw3OnTlH/ukLGo1t6/mn8q2Nyqf+qBP/8tfdWM4UkBwVYq0jP2l9PnNdJltmbilRX95yeEuGTBpSLr8fVX4snhaCo4MJjg4mdkAsYL3qOJFxojhh7F+1ny1fbyH1Y2ujvLEYgmOCqZNYhzoJfywBDQLw8Pb4o+orK7fEcmTHEX6f/Tu5Wbl41fQibmCcdbYdB2iCUOo8xmKh2V13EdGhA6v/8Q9W//3vbP/iC9q+8ALhKSlXdWyLh6W4sfh8RT1kHDX+I2vZwQ92QEQ4uvMou37Yxa4fdhXXixdJHKZdeCsLYwwB9a0N2s2GNAOsSePYnmPsX7WfzLRMa9vU+kw2z9hst1HfHp9AH2L6xtD85ubE9I3By89LE4RSV6N2TAzX/vvfpC9axJoxY1g0YgQNevWi9bPPUquea+ex8sw5Q82tR5i1fBa7fthV3IDtH+lP0/5NieoZRVSPKGo3qO3SONXVM8YQ2DiQwMaBNB/6xwzM+afzydqcxaENhziZeZJzZ85Zu9+G+pVcQvysvZaukCYIpS7CGEODnj2J7NyZzf/9Lxv//W8yBgwg5pZbaDZ8eIUmiuzfs9n01SY2f7mZ2DUHANgS5EPUtVF0eq4TUT2jCIkNueqGdVU5ePl6EdEmgog2EU49jyYIpS7Do0YNEh58kCaDB7P+//6PbdOmsW3qVBr27Uvze+8lKD7eKec9degUa/6zho2fbyRzfSYA9a6pR+YNTTkVH8Inf+2tDdDKqTRBKOUgv/BwUt54g6THHmPLp5+yffp09syZQ3jHjjS/7z7qpqRgjOHr1AxS9+aQV1BIp38s4pk+cQxp7fjVRtaWLH4d8ytpn6VRcLaABp0a0OftPsTfGE/tBrW51daHX5ODcjZNEEqVkV94OG2eeYaEBx9k+xdfsOXTT1l0//0ExceTfv29vLXTQl6BtSdRRs5pXphhHRTQkSSR+kkq3z38HcZiaD2iNdc8dg2hzUKd+nqUuhhNEEpdIe+AAJrffz9xw4ez+5tv2PzJJ7yfepTTfiVvpDqdX8DY+VsvmSAKCwqZ+/hcVo1fRVTPKG6cciO16tZy9ktQ6pL0GlWpq+Th7U30TTdx/ezZHPOzf5ft/pzSN8sVERFm3DmDVeNX0fGZjgybP0yTg3ILmiCUKifGYiEy0Nfuc7VzjzL/ttvYMmlSqWlPt87eysbPN9L9r93pPUYbnpX70HeiUuXomT5x+Hp5lNjm42l4oImFwoIC1owZw9c9ezLvlltIGz+eI5s38/PflhDaLJQuL3RxUdRK2ec2bRDGmAbAZCAcKAQmiMg7ro1KqbIpamd49sv15BUUUi/Q949eTE/cwfE9e9i3YAEZixeTNn48ae+/j09hTWoFx3NwaRJ127fHo0YNF78KpazcJkEA54CnRGSNMcYfWG2MWSgim1wdmFJlMaR1Paau2AvA5xeMuBrQqBEtHniAFg88wJnsbDKWLOG31/5LYc5aFj/0EJ6+voQlJ1O3fXvqtmtHUHx8iXGgrrYLbVWmv5vy5zYJQkQOAAdsj08YYzYD9QBNEKpK8gkJIfqGG9i3Poglf13EbdNSOLZlNZkrVrD2zTcB8KpVi7C2banbvj2rQ+N5Y/mRK+5CW5V9nZrBCzPS9HdTztyyDcIY0xhoDSx3bSRKOV/SsCS8a/mx4KWdxI98kgHffMMNixfTcexYGvXvz4k9e0gdO5Z/LdjG6QvmjyjqQlvdjZ2/ldP5JSfo0d/N1XO7BGGMqQV8BTwhIsftPD/SGLPKGLPq8OHDFR+gUuUsODqYO767g+Ppx5ncczJZW7PwDQujcf/+tH/1VQbOmcOQRYs45hdod/+Mo7msfftt9i1cSG5mZgVH7x4u1o34Ut2L1eW5TRUTgDHGC2tymCIiM+yVEZEJwASA5ORkBwe8Vcq9NezUkNu/uZ3pN0/no9Yf0eufvWj/SPviOYz96tYlMtCPDDsfeMHnTrH5k0+Qc+esZcPDCW3ZkpCkJEJbtiS4efMq3/AdGehr93dzsW7HyjFukyCMdRjKj4HNIvKWq+NRqqJF9Yhi1IZRzB4xm3mPz2PtJ2vp9mo34gbFYYzhmT5xvDAjrURViq+XB6/c2pmBf1/B0a1byV6/nqx168hav5698+cD1mlVA5s1K04aIYmJ+DdsWKVGfr3Y7+aZPnEujKryc5sEAXQC7gLSjDFrbdv+n4h858KYlKpQ/hH+3DHnDtKmpPHTX37i8yGfE946nO6vdWfwQOvMY3a70AKhSUmEJiURN2wYAKcPHyY7Lc2aMNatY8eMGfw+ZQpgHSYkJDHRutiShk9wsGtedDm4ZPdidcXcJkGIyC9A1flKo9QVMsaQNCyJhNsSWD9lPUv+toRpg6cR0SaCbq91o3X92mAxpbrQXsg3LIz6PXpQv0cPAArPneP4zp1krV9Pdloa2evXs3HCBKTQ2vBds359a8JISCA0KYmg+Hg8fStPFc2luherK+M2CUIpVZLF00Kru1uRdGcS6/9nSxSDphET4sOxdhEc7hpDWHxYGY7nSWBsLIGxscQMHQrAudxcjmzaZE0YaWlkrV3L3rlzATAeHgTGxpa4ygiIisLi4XGp06gqRBOEUm7O4mmh1T2tSLwzkU3TNzH5jZ8Inb+L8c3HE946nBa3tqD5Tc0Jjil7FZGnnx91kpOpk5xcvK2oaio7LY3sDRvYM3cu27/4wlq+Zk1CWrQokTT86tYtt9eq3IsmCKUqCQ8vDxLvSGTviZN4HDvL4zV8SJuSxg/P/8APz/9A3aS6xN8UT9PrmxLROqK4B1RZXVg1JYWFnNizp0TV1JZJkygs6jUVEUFYq1aEtm5NWKtWBMbGYvHyKrfXrVxHE4RSlVBB7RqkPJhCyp9SyNmTw+YZm9n81WYWv7aYxa8upmadmsT0jSGmXwzR10XjG3zlbQnGYiEgKoqAqCiaDB5sPf/ZsxzdupWstWvJWreOQ2vWsMdWNeXh60tIQgJhrVsT2qoVoS1bUiPQ/j0cyr1pglDKCSqykTSwUSAdnuxAhyc7cOrQKbbP3872udv5/dvfWTd5HcZiqHdNPaL7RBPVI4r619THw/vq2hE8atQo7jVV5NSBA2StXcvhtWvJSk1l08cfIwXWbqcBTZoQ2qqV9UqjVSsCoqIwFre7T1ddQBOEUlVIzTo1aXlXS1re1ZLCgkL2r9zP9nnWhPHTX37ip9d+wsvPi4ZdGhLVM4qoHlGEtwovlzkoakZEUDMigkb9+gHWBvDsjRutSSM1lfQffmDnDOv9r94BAdari1atCGvdmpCEBDz9/C567KwTZ9l39LQOxFfBNEEoVUVZPCzUT6lP/ZT6dH+tO6ePnmbPT3vYtWgXu37YxffPfg+AT5APjbs3JqqHNWGExoeWy010nn5+1G3Xjrrt2gHWmfNO7N5dfIVxODWV/UuWANYeU6FJSYR37EhEp04Et2hRPIpt1omz7Mo+RaFt3AQdiK/iaIJQqprwDfKl2ZBmNBvSDICTB0+y60drsti1aBdbZm4BoFZ4LaJ6RNG4R2Oa9GxCYOPyaT8wxhS3ZUTfcAMAZ3NyyFq/nsNr1nBw6dLiOTK8AgIIb9+e8E6dSD/iR6GUTFiOzPOtrp4mCKWqqVrhtUi8PZHE2xMBOLrrKLsW7WL3ot3sWrSLtM+s39IDowKtVxc9o4i6Nopa4eU3X3aNwEDqde1Kva5d4YknOJuTw8Flyzj4228c+PVX9n3/PWcH/NPuLbRXOhCfzhvhOE0QSikAgqKCCBoRRJsRbRARsjZnFVdHbf5qM6kfpwIQ1iKsuDqqcffG+AT6lFsMNQIDadS3L4369i2uknr9g7XkGu9SZev6Ggrz88vUpVbnjSgbTRBKqVKMMYQ1DyOseRjtH21PYUEhB1MPWhPGol2kfpzKiv9bgbEYItpEFFdHNejUAO+apT/MrzSGgKgoIursKdEGAeBVkEeXX77kq65/o3737jTo3ZuIzp0ve8yisZrOp9VVF6cJQil1WRYPC5HJkUQmR9Lp2U4U5BWQvjy9uP1i2bhl/DbmNyxe1obxJr2bENMnhoi2EVfdQyrU3zpUeVEvpnqBvjzVoznt+9Vm38KFpC9ezK7Zs/Hy96dVgyQymran8Fy7ElO1FrkwORTReSPs0wShlCozD28PGnVpRKMujej+WnfyTuWx95e9xVVSi19dzOJXFuMb7EuTXk2I7hNN9HXRBNQPuKLzhfrXKE4Uf9xjEkX9Hj0ozM/n4LJl7PnuO3LnLqDRpl+Z+fMkGvbuTaPrryesdeviey7q6bwRZaIJQil11bxrehPTJ4aYPjEA5GblsmPhDnYu2Mn2+dvZ+MVGAMKahxHdJ5qYfjE06toIzxpX/xFk8fIisksXIrt04d36/ai7J42hBTvZOWsW2z7/HL/wcKIGDyb6xht13ogy0gShlCp3fqF+xT2kRIRDGw6xY/4OdizYwcrxK1k2bhleNb2I7h1N7MBYmt3QDN+g0t/iz+9x5O1hoYGdMueb+khXoCsA+adOkWGrfto4YQIbP/qIuikp/LnzjYzdbsgrEJ034jKMSOWdtTM5OVlWrVpV9h0/ud7689455RuQUhXg1o+WApV3zoP83Hx2/biL37/9nW1ztnF833EsXhZi+saQcHsCcQPj8K7lXdzj6Pxv+xYDb93Sqswf6KcOHGDn11+zc8YMTu3fT14NP9KbpfDIq48SFFf9rh6MMatFJPmy5TRBKFW5VPYEcT4R4cDqA2yYtoEN0zZwIuMEXn5exA6MZUysP4fyzpXap16gL78+3+PKzldYSOby5Uz+50dE7EjFo/AcwS1aED10KI369cPb3/9qX1Kl4GiC0CompZTLGGOKe0f1HtObvb/sJW1qGpumb+JQo+ZgZ8iPq+lxZCwWwjt0YHU/8Dp9kpdDM9nx1Ves/MtfWPPPf9KwTx+ib7qJsDZtqtSc3VdKE4RSlUxVuHKwx1gMjbo2olHXRvR7tx9fvv49h86WvoKICCifG/PyfWsRN6w3sXfeyZENG9gxYwa758xh16xZ+DduTPSNNxI1aBC+YY7P2lfV6Hi7Sim34+Hlwf8bkoCvV8lhyT3yC4j7Yiuz7pvF3l/3Uh5V5MYYQhITaf/qq9y4eDEpb7yBT3Awa996i6979mTJ44+TsXhx8QRJ1YnbXUEYYzyAVUCGiAxwdTxKKdcoaoguuvvZ28NCtMAN7eqzcdpG1n6yltBmobS6rxUth7ekVt2rHyPK08+PJkOG0GTIEI7t3MnOGTPYNXs26T/8gG+dOtbnbrgB/4YNr/pclYHbNVIbY/4MJAMBl0sQ2kitVNVX1Chf5PMHO5B3Mo+NX2wk9eNU9v22D4unhdiBsbQd2ZboPtGXbT8oS0N/YX4+GT/9xI4ZMzjw889IYSF127cnatAgIrt1wye47HOBu1qlbKQ2xtQHrgfeAP7s4nCUUm7Ku5Y3re9rTev7WnN482FSJ6ayfvJ6tszcQljzMFKeTCFpWBKePuVzI16DXr1o0KsXuZmZ1u6yM2ey7KWXMBYLoS1bUu/aa6nXvTsBTZpUqcZtd2uDeBt4FrA/YIpSSl0gLD6M68Zex5P7nuSGT2/Aw9uDbx74hrcbvc1Pf/2JU4dPldu5/OrWJeHBBxk4dy59v/iCFg89xLmzZ1n71lvMGTSIb/r1Y82YMeRs315u53Qlt7mCMMYMAA6JyGpjTPdLlBsJjARoWE3qAZVSl+fh7UHSsCQS70xk94+7WfrmUha/uphfRv9Cy7tbkvJkCqFxoeVyLmMMwS1aENyiBUmPPELuwYNk/PQTGYsX8/tnn7Fl0iTCWrcm+uabadinD54+5TckekW6qjYIY0wE4CUie686EGNGA3cB5wAfIACYISLDLraPtkEoVfXZa4Nw1OHNh1k2bhnrJq+j4GwBsQNi6fhsR57dmAHGOKXL8JmjR9k1axbbp0/nxO7deAcEEDVoEDG33ELt6OhyP9+VqJA7qY0xm4FYEfG4bOGyHbc78LQ2UiulriZBFDl16BQrx69k5fsryc3KJbdJbbKvi+Lf7w/CWJzTZiAiHFq5ku1ffMG+hQspPHeOsLZtibn5Zhpedx0eNWo45byOcDRBXG0bxAvAfVd5DKWUcqqadWrS/bXuPLHnCfq91w/PY2dp8OFaxrcYT+rEVM7ZuSHvahljqNu+PZ3+9S+G/PgjrZ56itOHD7P0+eeZee21rBk7llMHDpT7ecuT23VzLQu9glCq6iuPK4hSxxz/KwFrMmm76hCZ6zLxj/TnmieuIfnBZGoEOO+bvRQWkrliBdunT2ffwoVgDA379CFm6FDqtGtXYT2gyrWbqzHGAiAihbb1cGAAsElEfruaQJVSqsJ5WDjeLoIH/30DOxfu5Nd//sr3z37Pz6//TOv7W3PN49cQ2Ciw3E9rLBbCU1IIT0nh1P79bJk0iZ2zZrFnzhz8GzUi+qabaDxwIH516pT7ua+Eo1VMc4DHAIwxtbDe6TwW+MkYM9xJsSmllFMZY4i+LprhPwzn/hX307R/U5a/s5x3m7zL9Fums2/pPqedu2ZkJG1feIEbFi+mw+jR+ISGWof36NGDH0aMYMfMmeSdOOG08zvC0QTRFlhke3wjcByoAzwAPO2EuJRSqkLVa1ePm6bexJ92/okOT3Vgx4IdTOw4kf+k/IcN0zZQcN68FOXJ08eHqEGD6D15MgPmzCHhoYc4lZHB8pdeYkbXrvzy5z+TvmgRBXl5Tjn/JWNzsJw/kGN7fB0wU0TyjTGLgPedEplSSrlA7Ya16T2mN91e6cbaSWtZ/vZyvrr9KwIaBND+0fa0eaCN3dnvykNA48YkPfooiY88QnZaGru//ZY9c+eyd/58vAMCaNinD40HDLAOR25x/n3OjiaIvUAnY8w3QB/gZtv2YCDXGYEppZQredfypv0j7Wk3qh2/z/mdZeOW8f1z3/PTX3+i1T2tuOZP1xDSNMQp5zbGEJqURGhSEm2eeYaDy5ax+9tv2fXtt2yfPh2/iAga9+9P44EDCWza1CkxgOMJ4i3gU+AksAdYYtveFUhzQlxKKeUWjMUQNzCOuIFxHFx7kGVvL2P1hNWsHL+S2AGxpDyZQuPujZ3WA8ni5UVkly5EdunCudxc0hctYte337L5v/9l08cfExgbS+MBA2jUvz81IyLK9dwOJQgR+cgYsxpoACws6s0E7ABeLteIlFLKTYW3CmfIf4fQ6x+9WDl+Jas+WMXkbyYT3iqca564hoTbEvCs4bwRjDz9/Gg8YACNBwzgTHY2e+bNY/ecOax96y3WjhtHneRkGl53HQ169y6XiY70PgillFtzyn0Q5TSvd/7pfNKmpLHs7WUc3niYWuG1SH44meSHkqkZVvOq43TUiT172D1nDnvmzuX4zp1gzCWTRbnfSW2MGWKOS51rAAAeVElEQVSMWWKMybItPxtjbriC16KUUlWCl68Xbe5vw6i0UQybP4zwVuEsfmUxbzd8m9kPzObg2oMVEod/o0YkPvwwA775hv6zZpEwahRnjxxh1RtvMPPaa/n+7rvZOmUKpw8fLtNxHb1R7ing78Bk4L+2zR2Az4wxL4vIv8p0VqWUqkKK7qeIvi6aw5sPs/yd5aybvI7U/6QSmRxJm5FtSLgtgRr+zh9/KTAmhsCYGJIeeYRj27ezd8EC9s6fz+q//53Vo0cT1qaNw8dytLLsaeBREfn3edsmGmNWAH8FNEEopRTW+SkGfDiAnqN7sv5/61kzYQ3fjvyWBX9eQMIdCbQd2ZbItpEVEkvtmBgSY2JIfPjhEsnCUY4miFrAj3a2/2h7Timl1Hl8g3y55rFraP9oezKWZ7B6wmrWf2pNGOGtw2k7si2JdyQ6deyn852fLHCwx5WjbRBfA0PtbL8JmO3gMZRSqtoxxlA/pT6DJw7mqQNP0f/9/kiBMGfUHN6MeJPZ988mfXk67thh6KJXEMaY8+eE3g48b4y5FijqUpBiW95yXnhKKVV1+NT2od3D7Ugelcz+lftZ/e/VbJi6gdSPUwltFkrC7Qkk3J7gtBvwyupSVUyPXbB+FIi1LedvuwdrO4RSSikHGGOo174e9drXo8+bfdgwbQNpn6Wx+LXFLH51MRFtI6zJ4tYEAuoHuCzOiyYIEYmqyECUUupSnDE9qDuoEVCDtiPb0nZkW46nH2fD5xvYMHUDC59eyMJnFtKoSyMSbk+g+dDm+IX6VWhszrvlTymlVJkE1A+g41Md6fhUR7J/z2bDNGuymDNqDnMfm0v0ddEk3J5A3OC4Cuky63CCMMbEYm2obgh4n/+ciOi0o0opVY5CYkPo9ko3ur7clcx1maRNTWPjtI3MvGsmnj6exA6MJeG2BGL6xuDl5+WUGBy9Ue564CsgFevcECuBaKAG8LNTIlNKKSepTNVVxhjCW4UT3iqcXqN7sW/pPtI+S2PT9E1smr4JT19PYvrEEDckjtgBsfiFlF81lKNXEH8F/iIio40xJ4C7gP1YR3hdesk9lVJKlQtjMTTs1JCGnRrS751+7F68my1fbylejIehUddGNLuhGc0GN6N2w9pXdT5HE0Qc8LntcT7gJyJnjDF/xTodabl0dTXGBAL/ARIAAe4TEU1ASil1AYunhSa9mtCkVxP6/V8/9q/ab00UM7cw7/F5zHt8HhFtI2g2pBnNbmhGWPOwMg9J7miCOAH42B4fAGKADbb9g8p0xkt7B5gnIkONMd5AxTbZK6VUJWSMoV67etRrV4+eb/Qka2sWW77ewtavt/Ljyz/y48s/EhwTbL2yGNLM4eM6miCWA52BTVivGN40xrQEbqCcqpiMMQFYJyC6B0BE8oCKn4RVKaUqudC4UDo/15nOz3XmxP4TbJ29lS0zt7Bs3DJ+G/ubw8dxNEH8mT/GXHoN6xzVNwG/254rD02Aw8AntuSzGviTiJwqp+MrpVS14x/pT/JD1jkqzuScYdt323jtztcc2tehsZhEZKeIrLc9zhWRUSKSJCJDRWTvlYdegifQBvhARFoDp4DnLyxkjBlpjFlljFl1uIxjmyulVHXmE+hD4h2JDpd3eMKgCpAOpIvIctv6l1gTRgkiMkFEkkUkOawcptRTSilln9skCBE5COwzxsTZNvXE2uahlFLKBdxtqI3HgCm2Hkw7gXtdHI9SSlVbbpUgRGQtcNmJtJVSSjmf21QxKaWUci+XvIIwxtQEBorINNv6eP64YQ6gAHhCu6IqpVTVc7kqpnuBHsA02/pdwAog17beEngEGOOU6JRS1V5lGlivqrlcFdNtwOQLtj0gIgNFZCDwHNYb5pRSSlUxl0sQTbHeLV0kB2u1UpFVQHx5B6WUUsr1LlfFFAAUFq2ISAM7+ztnpgqllFIudbkriH3Ape7Lbmkro5RSqoq5XIKYA7xmjPG58AlbD6dXbWWUUkpVMZerYhoN3AJsNca8xx/tEc2AR7EmmNHOC08ppZSrXDJBiMghY0xH4EPgH0DRdEQCLAAeFpFDzg1RKaWUK1x2qA0R2QP0M8YEYe3VBLBdRI44NTKllFIu5fBYTCJyFOtNckoppaoBHYtJKaWUXZoglFJK2aUJQimllF2aIJRSStmlCUIppZRdmiCUUkrZpQlCKaWUXZoglFJK2eVWCcIY86QxZqMxZoMxZqq9QQKVUkpVDLdJEMaYesDjQLKIJAAeWGe0U0op5QJukyBsPAFfY4wn4Afsd3E8SilVbblNghCRDOBfwF7gAHBMRBa4NiqllKq+3CZB2EaLHQxEAZFATWPMMDvlRhpjVhljVh0+fLiiw1RKqWrDbRIE0AvYJSKHRSQfmAF0vLCQiEwQkWQRSQ4LC6vwIJVSqrpwpwSxF0gxxvgZYwzQE9js4piUUqracpsEISLLgS+BNUAa1tgmuDQopZSqxhyeMKgiiMirwKuujkMppZQbXUEopZRyL5oglFJK2aUJQimllF2aIJRSStmlCUIppZRdmiCUUkrZpQlCKaWUXZoglFJK2aUJQimllF2aIJRSStmlCUIppZRdmiCUUkrZpQlCKaWUXZoglFJK2aUJQimllF3VL0Gs/wLSV8KeX2BcgnVdKaVUKdUrQaz/Ar55HArOWteP7bOua5JQSqlS3GpGuXL1yfWlt6Wv/CM5FMk/DbMehdWTSm6/d47zYlNKqUqgel1BXJgcLrddKaWqsap7BWHvCmBcgrVa6UK1G+gVg1JKXaDCryCMMRONMYeMMRvO2xZsjFlojNlm+xnklJP3fAW8fEtu8/K1bldKKVWCK6qY/gv0vWDb88APItIU+MG2Xv6SboGB74JHDet67QbW9aRbnHI6pZSqzCq8iklElhhjGl+weTDQ3fZ4ErAYeM4pASTd8keDtFYrKaXURblLI3VdETkAYPtZx8XxKKVUtecuCcJhxpiRxphVxphVhw8fdnU4SilVZblLgsg0xkQA2H4eulhBEZkgIskikhwWFlZhASqlVHXjLgliNnC37fHdwCwXxqKUUgrXdHOdCiwF4owx6caYEcA/gN7GmG1Ab9u6UkopF3JFL6bbL/JUzwoNRCml1CW5SxWTUkopN1N1h9pQqhLJz88nPT2dM2fOuDoUVYX4+PhQv359vLy8rmh/TRBKuYH09HT8/f1p3LgxxhhXh6OqABEhOzub9PR0oqKirugYWsWklBs4c+YMISEhmhxUuTHGEBISclVXpZoglHITmhxUebva95QmCKUUALVq1XLKcWfOnIkxhi1btjjl+EXy8vK49957SUxMpGXLlixevNhuuenTp9OiRQssFgurVq267HGXL19Oq1atSiw+Pj588MEHDsV15MgRevfuTdOmTenduzdHjx61W27SpEk0bdqUpk2bMmnSHxOYrV69msTERGJiYnj88ccREQBefvllkpKSaNWqFddddx379+93KJ4yEZFKu7Rt21auyMT+1kUpN7Fp0yZXhyA1a9Z0ynFvvvlm6dy5s7z66qtOOX6R9957T+655x4REcnMzJQ2bdpIQUFBqXKbNm2SLVu2SLdu3WTlypVlPs/8+fMlLi5OTp486VD5Z555RkaPHi0iIqNHj5Znn322VJns7GyJioqS7OxsOXLkiERFRcmRI0dERKRdu3by22+/SWFhofTt21e+++47ERE5duxY8f7vvPOOPPjgg3bPb++9BawSBz5j9QpCKXVRe/bsoWfPniQlJdGzZ0/27t0LwI4dO0hJSaFdu3a88sorF736OHnyJL/++isff/wx06ZNK95eWFjIww8/TIsWLRgwYAD9+/fnyy+/BKzfmLt160bbtm3p06cPBw4ccCjWTZs20bOn9XaqOnXqEBgYaPcKIT4+nri4uFLb9+/fT//+/S95jqysLB544AGmTJlCzZo1HYpr1qxZ3H23daCIu+++m6+//rpUmfnz59O7d2+Cg4MJCgqid+/ezJs3jwMHDnD8+HE6dOiAMYbhw4cX7x8QEFC8/6lTp5xSRam9mJRyM/OemMfBtQfL9ZjhrcLp+/aF07Bc3qOPPsrw4cO5++67mThxIo8//jhff/01f/rTn/jTn/7E7bffzocffnjR/b/++mv69u1LbGwswcHBrFmzhjZt2jBjxgx2795NWloahw4dIj4+nvvuu4/8/Hwee+wxZs2aRVhYGJ9//jkvvvgiEydOZOzYsUyZMqXUObp27cq7775Ly5YtmTVrFrfddhv79u1j9erV7Nu3j/bt2zv0WiMjI/nuu+8uWWbEiBE8/PDDtG3btnhbly5dOHHiRKmy//rXv+jVqxeZmZlEREQAEBERwaFDpYeay8jIoEGDBsXr9evXJyMjg4yMDOrXr19qe5EXX3yRyZMnU7t2bX788UeHXmdZaIJQSl3U0qVLmTFjBgB33XUXzz77bPH2om+yd9xxB08//bTd/adOncoTTzwBwG233cbUqVNp06YNv/zyCzfffDMWi4Xw8HCuvfZaALZu3cqGDRvo3bs3AAUFBcUfrs888wzPPPPMRWO977772Lx5M8nJyTRq1IiOHTvi6Vl+H3Effvghx48fLxXDzz//fNXHFlu7wvmMMRfdXuSNN97gjTfeYPTo0bz33nv85S9/uepYzqcJQik3cyXf9CtKWaoxsrOzWbRoERs2bMAYQ0FBAcYYxowZY/eDD6wflC1atGDp0qWlnrvcFYSnpyfjxo0r3t6xY0eaNm3qcLyXsmXLFl5//XWWLVuGxVKyZv5yVxB169blwIEDREREcODAAerUKT3dTf369Us0qqenp9O9e3fq169Penp6ie2RkZGl9r/jjju4/vrryz1BaBuEUuqiOnbsWNx2MGXKFDp37gxASkoKX331FUCJtoXzffnllwwfPpw9e/awe/du9u3bR1RUFL/88gudO3fmq6++orCwkMzMzOIPx7i4OA4fPlycIPLz89m4cSNgvYJYu3ZtqeXdd98FIDc3l1OnTgGwcOFCPD09ad68ucOvNSMjo7gN43x5eXnccccdjBs3rkR1T5Gff/7Zbly9evUCYNCgQcW9kiZNmsTgwYNLHaNPnz4sWLCAo0ePcvToURYsWECfPn2IiIjA39+fZcuWISJMnjy5eP9t27YV7z979myaNWvm8Gt1mCMt2e66aC8mVVW4Qy8mY4zUq1eveHnzzTdl165dcu2110piYqL06NFD9uzZIyIiv//+u7Rv317atWsnr732mkRGRpY6Xrdu3WTu3Lkltr3zzjvy0EMPSUFBgTz44IMSHx8vgwcPlr59+8qCBQtERCQ1NVW6dOkiSUlJ0rx5c5kwYYJD8e/atUtiY2OlWbNm0rNnT9m9e3fxcyNGjCjusTRjxgypV6+eeHt7S506deS6664TEZGVK1cWPz7fZ599Jp6entKyZcsSy1tvveVQXFlZWdKjRw+JiYmRHj16SHZ2dvH5RowYUVzu448/lujoaImOjpaJEycWb1+5cqW0aNFCmjRpIo888ogUFhaKiMiNN94oLVq0kMTERBkwYICkp6fbPf/V9GIycpFLvcogOTlZHOnHXMon11t/6pzUyk1s3ryZ+Ph4V4fhsNzcXHx9fTHGMG3aNKZOncqsWWWbxuXkyZPUqlWL7Oxs2rdvz6+//kp4eLiTIr689957j4YNGzJo0CCXxeAM9t5bxpjVIpJ8uX21DUIpVWarV6/m0UcfRUQIDAxk4sSJZT7GgAEDyMnJIS8vj5dfftmlyQGsPbZUSZoglFJl1qVLF9atW3dVx7jYnc7KfWgjtVJKKbs0QSillLJLE4RSSim7NEEopZSyq8IThDFmojHmkDFmw3nbxhpjthhj1htjZhpjAis6LqWqu8o+3Hd+fj533303iYmJxMfHM3r0aLvldu3axTXXXEPTpk259dZbycvLu+Rx33///RJDfSckJGCMYfPmzQ7FdbHhus8nIjz++OPExMSQlJTEmjVrAFi7di0dOnSgRYsWJCUl8fnnnzt0zvLiiiuI/wIXjiWwEEgQkSTgd+CFig5KKeUcU6dOpXPnzhe947q8TJ8+nbNnz5KWlsbq1av56KOP2L17d6lyzz33HE8++STbtm0jKCiIjz/++JLHfeSRR0rcIT1o0CDuvPNOh+9bGTVqFBMmTGDbtm1s27aNefPmlSozd+7c4ucnTJjAqFGjAPDz82Py5Mls3LiRefPm8cQTT5CTk+PQectDhScIEVkCHLlg2wIROWdbXQaUvp9dKVXhKtNw38YYTp06xblz5zh9+jTe3t4lhsQG6zf1RYsWMXToUKDk8NuzZ8/mlVdeueQ5lixZwhdffMH48eMdiulSw3Wfb9asWQwfPhxjDCkpKeTk5HDgwAFiY2OLx5OKjIykTp06HD582KFzlwd3vA/iPqBir6OUciOrR4/m6Nat5XrMoLg42r5Q9gvzyjTc99ChQ5k1axYRERHk5uYybtw4goODS5TNzs4mMDCweJTX84fPHjRo0CXvos7JyeHee+9l8uTJxYln69at3HrrrXbLL168+LLDdRe52HDfRSPZAqxYsYK8vDyio6MvGmN5c6sEYYx5ETgHlH4X/FFmJDASoGHDhhUUmVLVU2Ua7nvFihV4eHiwf/9+jh49SpcuXejVqxdNmjQpLmOv/t/REWpHjRrFsGHD6NSpU/G2uLg41q5de9F9HD3f5codOHCAu+66i0mTJpUaTdaZ3CZBGGPuBgYAPeUSA0SJyARgAljHYqqg8JSqMFfyTb+iuPNw35999hl9+/bFy8uLOnXq0KlTJ1atWlUiQYSGhpKTk8O5c+fw9PS86PDZF5o0aRK7d+/m008/LbH9clcQjg7XXb9+ffbt22e33PHjx7n++ut5/fXXSUlJuWys5cqREf3KewEaAxvOW+8LbALCynIcHc1VVRXuMJqrvTmpBw4cKJMnTxYRkU8++USGDBkiIiL9+/eXadOmiYjIRx99ZHffDz/8UEaOHFliW9euXWXJkiXyxRdfyPXXXy8FBQVy8OBBCQoKkunTp8vZs2clOjpafvvtNxERycvLkw0bNjgU/z/+8Q+55557pLCwUE6ePCnx8fGybt26UuWGDh0qU6dOFRGRBx98UN5//30RsY7y+vzzz5cqv2PHDqlXr57s2LHDoTgulJycLEuXLi2eU3rOnDmlynz77bfSt29fKSwslKVLl0q7du1EROTs2bPSo0cPGTdu3BWdW+TqRnN1RXKYChwA8oF0YASwHdgHrLUtHzpyLE0QqqpwhwRR2Yf7PnHihAwdOlSaN28u8fHxMmbMmOLn+vXrJxkZGSJi/cBv166dREdHy9ChQ+XMmTMiIjJ27Fj5+9//Xuq4I0eOlMDAwFLDfS9ZssShuC42XPcHH3wgH3zwgYiIFBYWysMPPyxNmjSRhISE4qHJP/3001JDjaempjp03iI63HdZ6XDfys3ocN+uH+572LBhjBs3jrCwMJfF4Aw63LdSqkJVxeG+//e//7n0/O5IE4RSqsx0uO/qQcdiUkopZZcmCKXcRGVuD1Tu6WrfU5oglHIDPj4+ZGdna5JQ5UZEyM7OxsfH54qPoW0QSrmBohuqKnKcHVX1+fj4lBjqo6w0QSjlBry8vIiKinJ1GEqVoFVMSiml7NIEoZRSyi5NEEoppeyq1ENtGGMOA3vK6XChQFY5Hau8aEyO0Zgc545xaUyOKc+YGonIZccUqdQJojwZY1Y5MjZJRdKYHKMxOc4d49KYHOOKmLSKSSmllF2aIJRSStmlCeIPE1wdgB0ak2M0Jse5Y1wak2MqPCZtg1BKKWWXXkEopZSyq9omCGNMsDFmoTFmm+1n0EXKjTHGbDTGbDbGvGvKMmu782JqaIxZYItpkzGmsatjspUNMMZkGGPec1Y8jsZkjGlljFlq+9utN8bYn1n+6mPpa4zZaozZbox53s7zNYwxn9ueX+7Mv1UZYvqz7X2z3hjzgzGmkbNjciSu88oNNcaIMcbpPXYcickYc4vt97XRGPOZq2Oy/f//aIxJtf0N+zstGEfmJa2KCzAGeN72+Hngn3bKdAR+BTxsy1Kguytjsj23GOhte1wL8HN1TLbn3wE+A95zg79dLNDU9jgS6zzogeUchwewA2gCeAPrgOYXlHkY2xzrwG3A507+3TgS07VF7xlglLNjcjQuWzl/YAmwDEh2dUxAUyAVCLKt13GDmCYAo2yPmwO7nRVPtb2CAAYDk2yPJwFD7JQRwAfrH6oG4AVkujImY0xzwFNEFgKIyEkRyXVlTLa42gJ1gQVOjMXhmETkdxHZZnu8HzgElPdkw+2B7SKyU0TygGm22C4W65dAT2dehToSk4j8eN57Zhlw5cN9lmNcNn/D+gXgjJvE9ADwvogcBRCRQ24QkwABtse1gf3OCqY6J4i6InIAwPazzoUFRGQp8CPWb58HgPkistmVMWH9ZpxjjJlhu8Qca4zxcGVMxhgL8CbwjBPjKFNM5zPGtMea5HeUcxz1gH3nrafbttktIyLngGNASDnHUdaYzjcCmOvEeIpcNi5jTGuggYh8WwHxOBQT1v+3WGPMr8aYZcaYvm4Q02vAMGNMOvAd8JizgqnSw30bY74H7M2E/qKD+8cA8fzxDWuhMaariCxxVUxY/2ZdgNbAXuBz4B7gYxfG9DDwnYjsK68vx+UQU9FxIoBPgbtFpLA8Yjv/8Ha2Xdgt0JEy5cnh8xljhgHJQDcnxlN8OjvbiuOyfckYh/W9XFEc+V15Yq1m6o71c+BnY0yCiOS4MKbbgf+KyJvGmA7Ap7aYyvv9XbUThIj0uthzxphMY0yEiBywfYjYu3S8AVgmIidt+8wFUrDWkboqpnQgVUR22vb52hbTFSeIcoipA9DFGPMw1jYRb2PMSRG5aENkBcSEMSYAmAO8JCLLrjSWS0gHGpy3Xp/Sl/tFZdKNMZ5YqwSOOCGWssSEMaYX1mTbTUTOOjEeR+PyBxKAxbYvGeHAbGPMIBFZ5aKYisosE5F8YJcxZivWhLHShTGNAPqCtZbDGOODdZymcq/+qs5VTLOBu22P7wZm2SmzF+hmjPE0xnhh/ablzComR2JaCQQZY4rq03sAm1wZk4jcKSINRaQx8DQw+WqSQ3nEZIzxBmbaYpnupDhWAk2NMVG2891mi+1isQ4FFomtddFVMdmqcj4CBlVAnbpDcYnIMREJFZHGtvfRMlt8zkoOl43J5musjfoYY0KxVjntdHFMe4GetpjisbaTOmcqQme2yLvzgrUe+Adgm+1nsG17MvAf+aNHwUdYk8Im4C1Xx2Rb7w2sB9KA/wLero7pvPL34PxeTI787YYB+cDa85ZWToilP/A71vaNF23b/or1ww2s/7zTge3ACqBJBby3LxfT91g7WxT9XmY7OyZH4rqg7GKc3IvJwd+VAd6y/f+nAbe5QUzNsfauXGf7+13nrFj0TmqllFJ2VecqJqWUUpegCUIppZRdmiCUUkrZpQlCKaWUXZoglFJK2aUJQimllF2aIFS1ZoxZZBtaumg5aoz5zhjT4PJ7OyWGfGMdxvyBijq/UhejCUJVd62B/wdEYB0U7UagHfC6i2KIwTqC50e2O56VchlNEKraMsZEA4HAYhE5KCL7ReRHYANQs4JjmGeLYQ/Wu/cN0KIiYlDqYqr0YH1KXUZb4BzW4Qqwjbd1M9YriIEVGMNxrEOnYIwJB8YChcCaCopBKbs0QajqrC3W8bYO20YQ9QVygBtsVxIVFUMt4JhtyGtfIA94SkScOQijUpelVUyqOmsLzABa2ZaOWEfTnGAbQtkhxpjXL2jotrd0v0QM/7GdvzMwH/i3iLx9Fa9LqXKhg/WpassYkw28LiLjztvWHessgm1FZI1tvo2zWOcIDsE6mueKC44TinU8/kvZK3amhrXF8GcRmWRbbwTsAlqKSNp55S4bh1LlTa8gVLVkjIkCgrFOSH++JrafRePrJwHrRKQd1gl2Sk2rKiJZIrLlMou95FAUQ9p5x9pji+muC4pfNg6lypsmCFVdtbX9PGiMCTfGNDHG3An8A5gi1ulTawE1gDG2shuBoHKOoZDSk1AtxDqbIQAVEIdSdmkjtaquihJE0YfzMayT+rzMH9O3JgEbReScbb01tt5G5RjDDhE5fcH2hcBzxpgWIrKxAuJQyi5tg1DqIowxDwF/xno/Qi2sDcg3ikh6dYxDVT96BaHUxSUBXwG/YL1x7jkXfSi7SxyqmtErCKUuwhjzC3CXiOzSOFR1pI3USl1cQ2C3q4PAfeJQ1YxeQSillLJLryCUUkrZpQlCKaWUXZoglFJK2aUJQimllF2aIJRSStmlCUIppZRdmiCUUkrZpQlCKaWUXZoglFJK2fX/AdnBMBil9SErAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import pandas as pd\n", "import numpy as np\n", "from scipy import stats\n", "\n", "def read_and_error(file):\n", " stars = pd.read_csv(file)\n", " stars_G = stars['phot_g_mean_mag']\n", " stars_Bp = stars['phot_bp_mean_mag']\n", " stars_Rp = stars['phot_rp_mean_mag']\n", " \n", " Av = 2.35 # This is the estimated extinction in this region\n", " Ag = 0.86117 * Av\n", " Abp = 1.06126 * Av\n", " Arp = 0.64753 * Av\n", " stars_Bp = stars_Bp - Abp \n", " stars_Rp = stars_Rp - Arp\n", " stars_G = stars_G - Ag\n", "\n", " stars_error_G = 2.5 / stars['phot_g_mean_flux_over_error'] / np.log(10.)\n", " stars_error_Bp = 2.5 / stars['phot_bp_mean_flux_over_error'] / np.log(10.)\n", " stars_error_Rp = 2.5 / stars['phot_rp_mean_flux_over_error'] / np.log(10.)\n", " \n", " stars_rpgeo = stars['rpgeo']\n", " stars_rpgeo_low = stars['b_rpgeo_x']\n", " stars_rpgeo_up = stars['B_rpgeo_xa']\n", " \n", " stars_colour = stars_Bp - stars_Rp\n", " stars_error_colour = np.sqrt( stars_error_Bp**2. + stars_error_Rp**2. )\n", " stars_absG = stars_G + 5 - 5 * np.log10(stars_rpgeo) \n", " stars_absG_low = stars_G - stars_error_G + 5 - 5 * np.log10(stars_rpgeo_low)\n", " stars_absG_up = stars_G + stars_error_G + 5 - 5 * np.log10(stars_rpgeo_up)\n", " \n", " plt.errorbar(stars_colour,stars_absG,xerr=stars_error_colour,yerr=[stars_absG_low,stars_absG_up],marker='o',linestyle='None')\n", " \n", " \n", "file1 = './stars_eDR3_dis.csv'\n", "read_and_error(file1)\n", "\n", "file2 = './hen237_dis.csv'\n", "read_and_error(file2)\n", "\n", "# Reading the PARSEC isochrone file\n", "file ='./iso.dat'\n", "Z, la, m, l, Gm, Bp, Rp = np.loadtxt(file,usecols=(0,2,5,9,25,26,28),unpack=True)\n", "\n", "colors = ['purple','brown','blue','blue','g','r']\n", "\n", "metalkeep = 0.0027383 # a metal-poor isochrone of 1.26 Gyr\n", "agekeep = 9.10000\n", "\n", "cond = np.logical_and(la == agekeep, Z == metalkeep)\n", "Gmm = Gm[cond]\n", "BRp = Bp[cond] - Rp[cond]\n", "\n", "plt.plot(BRp, Gmm, c=colors[0],label='Log Age=9.1; Z=0.003')\n", "\n", "metalkeep = 0.020322 # a metal-rich isochrone of 100 Myr\n", "agekeep = 8.00000\n", "\n", "cond = np.logical_and(la == agekeep, Z == metalkeep)\n", "Gmm = Gm[cond]\n", "BRp = Bp[cond] - Rp[cond]\n", "\n", "plt.plot(BRp, Gmm, c=colors[1],label='Log Age=8.0; Z=0.02')\n", "\n", "plt.xlim(-0.9,0.9)\n", "plt.ylim(13.,-2.)\n", "plt.xlabel(r'$B_p-R_p$',fontsize=14)\n", "plt.ylabel('G abs.',fontsize=14)\n", "plt.legend(loc='lower right');\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, we can see that because of the large errors on the distance, it is very difficult to say anything about the positions of the stars from the eroteme in the CMD. The error on the distance of the planetary nebula is also very large, but it appears that the central star is much bluer than the other (which makes sense as it is a hot core of a star) and it occupies in the CMD the position of stars that become white dwarfs." ] }, { "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 }