{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Modelling and fitting two unresolved emission lines with a bayesian approach\n", "\n", "We will show how to model a spectrum with two superimposed lines and then try to retrieve the modelling parameters. this example is based on the preliminary examples :\n", "\n", "1. [Modelling and fitting one emission line](./script_example_model+fit_1_line.ipynb)\n", "2. [Modelling and fitting two resolved emission lines](./script_example_model+fit_2_lines.ipynb)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import orb.fit\n", "import pylab as pl\n", "import numpy as np\n", "from orb.core import Lines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Third step: modelling and fitting a spectrum with two unresolved lines (classic fit and bayesian fit)\n", "\n", "Now the two lines are set to nearly the same velocity but the other parameters are unchanged." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Model" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "incident angle theta (in degrees): 15.445939567249903\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD4CAYAAAD//dEpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1d348c93liyTlSwkQJAEBFnDYgRR3IuitW51wbrUXR+r1F9brd3UWnlc+tT6tPpoUSm27lpbcaWCdd8AQWTfEiCBTDbIvs3M+f0xMyGEyTZLMoHvm1dembn33HO/CTP5zrn3LGKMQSmllOqOpb8DUEopNTBowlBKKdUjmjCUUkr1iCYMpZRSPaIJQymlVI/Y+juAzmRkZJjc3Nz+DkMppQaUlStXVhhjMiNRd9QmjNzcXFasWNHfYSil1IAiIjsiVbdeklJKKdUjmjCUUkr1iCYMpZRSPRK19zCUOty0trZSXFxMU1NTf4eiBoC4uDhycnKw2+19dk5NGEpFieLiYpKSksjNzUVE+jscFcWMMVRWVlJcXExeXl6fnVcvSSkVJZqamkhPT9dkobolIqSnp/d5a1QThlJRRJOF6qn+eK2EJWGIyEIRKRORtZ3sv0xE1vi+PhORyeE4r1LRYnXZatZXru/vMJSKqHC1MBYBc7rYXwicZIzJB34HLAjTeZWKCvd+cS9/WPGH/g4jZImJiQDs3r2bCy+8MCx1vvLKK0yYMAGLxaKDcQe4sCQMY8xHQFUX+z8zxuz1Pf0CyAnHeZWKBsYYimuL2VO/p79DCZuhQ4fy6quvhqWuiRMn8tprr3HiiSeGpT7Vf/rjHsa1wDuBdojIDSKyQkRWlJeX93FYSgWnqqmKRlcjznonh8oKlkVFRUycOBGARYsWccEFFzBnzhxGjx7NHXfc0Vbu3//+NzNnzmTatGlcdNFF1NXVHVTXuHHjOOqoo/osdhU5fdqtVkROwZswZgXab4xZgO9yVUFBwaHxzlOHvJK6EgBaPC3sa97HoLhBIdf52zfWsX53Tcj1tDd+aDJ3f29CUMeuXr2aVatWERsby1FHHcWtt95KfHw89913H0uXLiUhIYEHH3yQhx9+mLvuuiuscavo0WcJQ0TygaeAM40xlX11XqUizZ8wAJwNzrAkjGhz2mmnkZKSAsD48ePZsWMH+/btY/369Rx//PEAtLS0MHPmzP4MU0VYnyQMETkCeA24whizuS/OqVRfaZ8wSutLGZs2NuQ6g20JREpsbGzbY6vVisvlwhjD7NmzeeGFF/oxMtWXwtWt9gXgc+AoESkWkWtF5CYRuclX5C4gHfg/EVktItpVQh0yimuLsVm8n72c9c5+jqbvHHvssXz66ads3boVgIaGBjZv1s+Dh7KwtDCMMZd2s/864LpwnEupaFNSV8KYQWPYVLUJZ8PhkzAyMzNZtGgRl156Kc3NzQDcd999jBkz5oBy//znP7n11lspLy/nu9/9LlOmTGHJkiX9EbIKkURrr46CggKjfbbVQHDWa2cxIX0Cq8tXMz17OvNnzQ+qng0bNjBu3LgwR6cOZYFeMyKy0hhTEInz6dQgSoXA7XGzp34PwxKHke3IprS+tL9DUipiNGEoFYKyhjJcHhfDkoaRlZB1WF2SUocfTRhKhaC4rhiAYYnDyHJkHVKD95TqSBOGUiHwd6nNScwhOyGbJncT1c3V/RyVUpGhCUOpEJTUlSAIQxKGkOXIAtDLUuqQpQlDqRCU1JaQnZCN3WonK0EThjq0acJQKgQldSUMSxwG0NbCGMg9pSIxvfntt9/O2LFjyc/P5/zzz2ffvn1hqVf1PU0YSoWguK64LWFkxmdiFeuAThh+4ZzefPbs2axdu5Y1a9YwZswY7r///rDUq/qeJgylgtTibqG8oZxhSd6EYbVYyYjPOCQuSYVzevPTTz8dm807qcSxxx5LcXFx3/wQKuz6dHpzpQ4lu+t2YzDkJO5fDyxsYzHeuRNKvw29nvayJ8GZDwR1aLimN1+4cCGXXHJJsD+B6meaMJQKkr9Lrf+SFEC2I5vNew+9CfjCMb35/PnzsdlsXHbZZX0Sswo/TRhKBSlQwshKyOLjko8xxiAiwVceZEsgUkKd3vyZZ57hzTffZNmyZaH9XlS/0nsYSgWpuK6YGEsMmY7Mtm1ZjiwaXY3UtIR3tbxo1NPpzd99910efPBBFi9ejMPh6OswVRhpwlAqSMW1xQxNHIpF9r+NDqexGO2nN8/Pz+fYY49l48aNB5W75ZZbqK2tZfbs2UyZMoWbbropQG1qINBLUkoFqf0YDL9sRzbgHYsxZtCYQIdFNX8vp9zcXNauXQvAVVddxVVXXdVW5s0332x7fOqpp7J8+fIu6/S3QNTApy0MpYIUMGEkeBPG4dDCUIcfTRhKBaGupY7q5uq2MRh+GfEZWMRyWC3Vqg4fmjCUCkKgHlIANouNjPiMQ2K0t1IdhSVhiMhCESkTkbWd7BcR+ZOIbBWRNSIyLRznVaq/+NfBaD9ozy/bka2XpNQhKVwtjEXAnC72nwmM9n3dADwepvMq1S9KagO3MCCMo72VijJhSRjGmI+Aqi6KnAv8zXh9AaSKyJBwnFup/lBSV0KCPYGU2JSD9mU5siitL9WV99Qhp6/uYQwDdrV7XuzbptSA5O8hFWjUcnZCNo2uRmpba/shstCEa3rzRx99lCOPPBIRoaKiIiyxnXzyyaxYsQKAs846K+hp0h9++GHGjx9Pfn4+p512Gjt27AhLfL0Rzunj+1JfJYxAcwEc9PFLRG4QkRUisqK8vLwPwlIqOIG61Pq1rbw3gHtKhTq9+fHHH8/SpUsZMWJEj8rfc889LFq0qMf1v/3226SmpgYV29SpU1mxYgVr1qzhwgsvPGD23UA++OCDA8ahhEM4p4/vS32VMIqB4e2e5wC7OxYyxiwwxhQYYwoyMzM77lYqKhhjuk4Yh8Bo71CnN586dSq5ubkRiy83N5eKigqKiooYN24c119/PRMmTOD000+nsbERgG3btjFnzhyOPvpoTjjhhLZR6KecckrbFCXhmm79qquuYt68eRx33HGMHDmyLRkYY7j99tuZOHEikyZN4qWXXgIO/P2uW7eO6dOnM2XKFPLz89myZQsAzz77bNv2G2+8EbfbHXKcoeqrkd6LgVtE5EVgBlBtjNnTR+dWKqyqmqpodDWSk3RwDyk4cLR3sB786kE2Vh08zUYoxqaN5efTfx7UseGa3jwStmzZwgsvvMCTTz7JxRdfzD/+8Q8uv/xybrjhBp544glGjx7Nl19+yc0338z7779/wLFPP/00Z555Zlji2LNnD5988gkbN27knHPO4cILL+S1115j9erVfPPNN1RUVHDMMcdw4oknHnDcE088wY9//GMuu+wyWlpacLvdbNiwgZdeeolPP/0Uu93OzTffzHPPPceVV14ZlliDFZaEISIvACcDGSJSDNwN2AGMMU8AbwNnAVuBBuDqcJxXqf7Q2RgMvwxHBoIM6BZGR+GY3ryjb7/9liuuuAKA0tJSYmJieOSRRwBYtmwZ6enpPaonLy+PKVOmAHD00UdTVFREXV0dn332GRdddFFbuebm5gOOe/bZZ1mxYgUffvhhwHpnzJhBc3MzdXV1VFVVtZ3jwQcf5Iwzzjio/HnnnYfFYmH8+PE4nd7/+08++YRLL70Uq9VKVlYWJ510EsuXLyc/P7/tuJkzZzJ//nyKi4u54IILGD16NMuWLWPlypUcc8wxADQ2NjJ48OAe/T4iKSwJwxhzaTf7DfCjcJxLqf7WXcKwW+zelfdCuIcRbEsgUkKd3jyQSZMmsXr1asB7DyM3NzeoewUdY2tsbMTj8ZCamtpWf0dLly5l/vz5fPjhhwcc396XX34JeO9hLFq0qNt7LO3r8feQ60lPuR/84AfMmDGDt956izPOOIOnnnoKYww//OEPo245Wx3prVQvdZcwwHvj+1BqYQTS0+nN+0NycjJ5eXm88sorgPcP9zfffAPAqlWruPHGG1m8eHHEP7WfeOKJvPTSS7jdbsrLy/noo4+YPn36AWW2b9/OyJEjmTdvHueccw5r1qzhtNNO49VXX6WsrAyAqqqqfunN1ZEmDKV6qbi2mLS4NBz2ztd2yE7IPuSnB+lqevM//elP5OTkUFxcTH5+Ptddd12fx/fcc8/x9NNPM3nyZCZMmMDrr78OwO23305dXR0XXXQRU6ZM4ZxzzolYDOeffz75+flMnjyZU089lYceeojs7OwDyrz00ktMnDiRKVOmsHHjRq688krGjx/Pfffdx+mnn05+fj6zZ89mz57+v+0r0Tq4qKCgwPj7XCsVTa7/9/XUt9bz/Hef77TMA189wL+2/osvfvBFj+vdsGED48aNC0eI6jAR6DUjIiuNMQWROJ+2MJTqpa661PplObKob62nrqWuj6JSKvI0YSjVC26Pmz31e7pNGLouhjoUacJQqhfKGspweVwHrYPRkX+0d2/vY0TrJWIVffrjtaIJQ6le8E9r3u0lqSBGe8fFxVFZWalJQ3XLGENlZSVxcXF9el5d01upXvB3qQ20DkZ7g+MHewfv9WIshr9Xkc6jpnoiLi6OnJyuX4fhpglDqV4oqStBEIYkdD07v91qJz0+ndKGnl+Sstvt5OXlhRqiUhGjl6SU6oWS2hKyErKwW+3dls1yZA3oGWuV6kgThlK90JMutX6Hw2hvdXjRhKFULxTXFgdOGNs/gB2fH7ApOyFbWxjqkKIJQ6keanY3U9ZYFviG95JfwT9vgHY9nLISsqhtraW+tb4Po1QqcjRhKNVDu+u8a34FHINRuwf27YSd+6cCORRW3lOqPU0YSvVQp7PUulqgodL7eM2LbZvbBu/1oqeUUtFME4ZSPVRS20nCqPdOQY0tHtb9E1qbgHbTg2gLQx0iNGEo1UMldSXYLXYGOzqsoVDrSwjTroCmatiyBKCtnLYw1KFCE4ZSPVRcV8zQxKFYpMPbpta3TkH+JZCYBWteBiDGGkNaXJq2MNQhQxOGUj3U6RiMOl8LInkoTLoINi+BhirA17VWx2KoQ4QmDKV6qNOEUesEBBIGe1sZnlZY9xrgvfF9qK+8pw4fYUkYIjJHRDaJyFYRuTPA/iNE5D8iskpE1ojIWeE4r1J9pa6ljurm6s5bGAmZYLVB9iQYPB6+eQnQ0d7q0BJywhARK/AYcCYwHrhURMZ3KPZr4GVjzFRgLvB/oZ5Xqb7U1qU24BgMJyR5u9AiAvkXQ/FXULnNO3ivpZaG1oY+jFapyAhHC2M6sNUYs90Y0wK8CJzboYwBkn2PU4DdYTivUn3Gvw5GwFHedaWQmL3/+aSLAYE1L7d1rdWeUupQEI6EMQzY1e55sW9be/cAl4tIMfA2cGsYzqtUn+l0DAYc2MIASBkGeSfAmpfIitfR3urQEY6EIQG2dVwy7FJgkTEmBzgL+LtIx76JICI3iMgKEVmhi8ioaFJSV0KCPYHU2NQDd3jc3oF77VsYAPlzYW8h2TXeRKH3MdShIBwJoxgY3u55DgdfcroWeBnAGPM5EAdkdKzIGLPAGFNgjCnIzMwMQ2hKhYe/h5RIh89H9RVgPJDUIWGM+x7Y4hm8+T1AWxjq0BCOhLEcGC0ieSISg/em9uIOZXYCpwGIyDi8CUObEGrA6HYMRmLWgdvjkmHsWcSuf5202EF6D0MdEkJOGMYYF3ALsATYgLc31DoRuVdEzvEV+ylwvYh8A7wAXGV0pXs1QBhjuhmDwcEtDPBelmrcS5Y1XlsY6pAQljW9jTFv472Z3X7bXe0erweOD8e5lOprVU1VNLoayUnqpIcUHNzCABh1KiRkktVUx269h6EOATrSW6ludDqtOexvYQRKGFYbTLyQrOo9OHW0tzoEaMJQqhtdJoy6UohLBXtc4IPzLya7tZXqlhoaXY0RjFKpyNOEoVQ3um5hlAa+f+E3dCpZ8d4ef3ofQw10mjCU6kZxbTFpcWk47I6Dd9Y5A1+O8hMhK+8UAJzObyIUoVJ9QxOGUt3otIcU+EZ5d9HCALInXARA6ZZ3wh2aUn1KE4ZS3eg0YRjjm0eqixYGMDh7CgDO4i+8xyg1QGnCUKoLbo+bPXV7AieMxr3gbum2hRFniyPVGo+zZR+UfB2hSJWKPE0YSnXB2eDEZVyBpzWv66JLbQfZSTk47TGw5sUwR6hU39GEoVQXuu0hBd22MACyEodS6kiFtf8Ad2s4Q1Sqz2jCUKoLxbVdrIPRljCGdFtPliMLp1WgoRK2LgtniEr1GU0YSnWhpK4EQRiSECApdDUtSAdZCVnsczXQFJ8G374S5iiV6huaMJTqQkldCVkJWdit9oN31johJhFiE7utx7/ynnP4NChbH+4wleoTmjCU6kKXYzB60KXWL8vhW3nPkQrVxeEKT6k+pQlDqS7sqd/D0IShgXf2YNCeX1vCiI2H5hpoqglXiEr1GU0YSnXC7XFT3lDedjnpIL1pYST4EobN6t1QUxKOEJXqU5owlOpERWMFbuNuax0cpBctjHhbPCmxKZTi8W7Qy1JqANKEoVQnnL5FjwK2MJprobW+xy0M8HWt9fimONeEoQYgTRhKdcKfMPyXkw7Q1dKsnchyZOFsrgax6CUpNSBpwlCqE6W+VfICXpLqxRgMv+yEbEobnN6BftrCUAOQJgylOuGsdxJrjSU1NvXgnb2YFsQvy5HF3ua9NKcM04ShBqSwJAwRmSMim0Rkq4jc2UmZi0VkvYisE5Hnw3FepSKptKGULEcWInLwzl5MPOjnv7RVlpihl6TUgGQLtQIRsQKPAbOBYmC5iCw2xqxvV2Y08AvgeGPMXhEZHOp5lYo0Z72z8y61taVgjYX4QT2uz19XqSOF4dUl3rUxAiUjpaJUOFoY04GtxpjtxpgW4EXg3A5lrgceM8bsBTDGlIXhvEpFlLPB2XmXWv/SrL34g++vqzTWAe5mqK8IR5hK9ZlwJIxhwK52z4t929obA4wRkU9F5AsRmROoIhG5QURWiMiK8vLyMISmVHDcHjdlDWWBe0iBt4WR1PPLUdButLfV97ar0fsYamAJR8II9BGr4zqUNmA0cDJwKfCUiBx0J9EYs8AYU2CMKcjMzAxDaEoFp7KpErdxk+3obJS3s1f3LwAcdgdJMUk4xe3doDe+1QATjoRRDAxv9zwH2B2gzOvGmFZjTCGwCW8CUSoqtXWp7bKF0fMeUn7ZCdmUuhq8T6r1xrcaWMKRMJYDo0UkT0RigLnA4g5l/gWcAiAiGXgvUW0Pw7mVioguR3m3NkHTPkjsfcIYmjCU4sYysMXpJSk14IScMIwxLuAWYAmwAXjZGLNORO4VkXN8xZYAlSKyHvgPcLsxpjLUcysVKc563yjvgIP2/KO8e3dJCiAvJY+dNTtxJ+vgPTXwhNytFsAY8zbwdodtd7V7bICf+L6Uinql9aWdD9prSxjdL83aUV5KHi2eFnanDPN2rVVqANGR3koF4O9SG3DQXm3vpwXxy03OBaAoPkUH76kBRxOGUgGU1peSlZCFMQZvA7mdIKYF8ctLyQOgMDYGaveA2xVqqEr1GU0YSgXgbHCS7cjmqr8u5zevrz1wZ10piBUcGb2ud1DcIFJjUykUDxiPN2koNUBowlCqA/9Ke1kJWXy9cy//Xuc8sJVR64TEwWAJ7u2Tm5xLkbve+0QvS6kBRBOGUh1UNlXiMi6SbRnUNrkoq21mR2XD/gK9WJo1kLyUPAqbfNOCaE8pNYBowlCqA3+XWotnfw+pLwvb9QLvxdKsgeSl5FHZUk2NRTRhqAFFE4ZSHZQ2eG9qu1qS27Z9WVi1v0CILYy2nlKOQXpJSg0omjCU6sDfwqivTwLguFHpfOVPGG6Xd5bZEFsYAIXJ6drCUAOKJgylOnA2eFfaK6+2khRnY/b4LIr3NlKyrxHqywATUgtjWNIwbBYbhXEJmjDUgKIJQ6kOSuu9K+0V721k+CAH0/PSAFheWBXSGAw/u8XO8KThFNmseklKDSiaMJTqwNngJCvBlzDS4hmbnUxSnM17H6NtadbgEwZAXnIehbRAQyW0NHR/gFJRQBOGUh04651tLYycQQ6sFuGY3DRvT6m2Fkbwl6TANwmhqxYXQE3H1QCUik6aMJRqx7/SXrItk8ZWN8MHxQMwPS+N7eX11Ff6LiElhLYsfW5KLi7jocRm02nO1YChCUOpdqqaqnAZFza8YzCGpzkA2u5jVJbuBEc62GJCOk9bTym7XW98qwFDE4ZS7fhX2vO0HpgwJg1LId5upaGqJOT7F7B/LEZhjE1X3lMDhiYMpdrxr7TX1JgIwLBU7yUpu9XC0SMGIXXOkO9fAKTEppAWl0ZRfJJeklIDhiYMpdrxJ4zqWgfpCTEkxO5fY2x6XhpJrkpa4kO7f+GXl5JHYWycXpJSA4YmDKXaKa0vJcYSg3OfjRzf5Si/6bmpZFJNiSslLOfKS8mj0GL0kpQaMDRhKNWOs947BqNkb2NbDym/Keke7OJmS0NCWM6Vm5zLPtzsrdsNHRdpUioKacJQqp3ShlKyHNmU7POOwWgvrqkcgFX74sJyLn9PqSLTDE37wlKnUpEUloQhInNEZJOIbBWRO7sod6GIGBEpCMd5lQo3Z72TFHsGrW7D8LQDWxjUeXtQrayMob459KVVD+xaq5elVPQLOWGIiBV4DDgTGA9cKiLjA5RLAuYBX4Z6TqUiwWM8lDWUEcMgAIZ3aGFQ670hvseTwtc794Z8vqEJQ4mx2Ciy23VOKTUghKOFMR3YaozZboxpAV4Ezg1Q7nfAQ0BTGM6pVNhVNnpX2hO3L2GkdUwY3vW3KyWNL7dXdTy816wWK0ck5lBot0H1rpDrUyrSwpEwhgHtX+3Fvm1tRGQqMNwY82ZXFYnIDSKyQkRWlJeXhyE0pXrO36W2pTkJERia2uFeRZ0TYlMYPTRj//oYIcpLPZLCmBi9JKUGhHAkDAmwra3Lh4hYgD8CP+2uImPMAmNMgTGmIDMzMwyhKdVz/lHedfWJZCXFEWuzHligthSSspiel8bqXftoanWHfM681JEU22y0agtDDQDhSBjFwPB2z3OA9tNvJgETgQ9EpAg4FlisN75VtPG3MCqr4w++4Q3eFkZiFjPy0mlxe/hmV+g9m3KTc3EL7KrZEXJdSkVaOBLGcmC0iOSJSAwwF1js32mMqTbGZBhjco0xucAXwDnGmBVhOLdSYeOsdxJjiaF0r/XgG97ga2Fkc0xuGiKE5bLUyJSRABT61hFXKpqFnDCMMS7gFmAJsAF42RizTkTuFZFzQq1fqb5SWl/KYEcWpdVN5HQYtIcxbS2MFIedo7KS+Koo9ISRm5ILQGFrDXg8IdenVCTZui/SPWPM28DbHbbd1UnZk8NxTqXCzdngJDUmA4/hoGlBaKoGV1Pb0qwz8tJ4ZWUxrW4Pdmvwn7sS7AkMtiVQaKuD+vKwTGyoVKToSG+lfErrS0mwZgABxmB0WJp1el46DS1u1pZUh3zePEe2dyyGTkKoopwmDKXYP2jP6vGvg9HhklSHpVn9CyqF4z5GbkoehXY7RntKqSinCUMp9q+052pJxmoRspMDjMGAthZGZlIsIzMTwpIw8tLHUWu1UFm1NeS6lIokTRhKsX8MRmNjEkNT47B1vC/RoYUB3vsYXxVV4faENtNsXsYEAAr3bg6pHqUiTROGUni71ALsq3UE7lJb5wRbPMQmt22anpdGbZOLTaW1IZ07z9e1tqhO72Go6KYJQym805oDlO2N72IMRhbI/okNpuelA/BlYWVI585KyCIOobCpIqR6lIo0TRhK4W1h2C12KqptB4/BAN8YjOwDNg1LjWdYanzI9zEsYiHX4qDQVR9SPUpFmiYMpfC2MNJiMwE5eJZaaBvl3dGMkWl8VViFCXHFvLy4DIosHnC1hFSPUpGkCUMpvC2MJJt3wstO55EKlDDy0qisb2FbeWitg9zEYZTYrDTv0zmlVPTShKEU3lHenS6c1NIAzTWQePAobP99jFAvS+WlHokRYUfpqpDqUSqSNGGow57HeHA2ODGuVGJsFjISYw8sUOfvUntwCyM33UFmUixfhXjjO29wPgBFletDqkepSNKEoQ57VU1VuDwumpsSyRkUj8XSYYmXWv+gvYNbGCLC9Lw0vgzxPsYRQ44GoLB6e9B1KBVpmjDUYc8/BqOuPrGTMRidtzAAjs1LY091E8V7G4OOweHIYIjLQ2H9nqDrUCrSNGGow55/lHdFZwsn1R44LUhH+8djhHgfQ2Ioag19USalIkUThjrs+Qft1dQmkNNZC8NiB0dawONHD04k1WHny+2h3cfItSdT6GkKuYuuUpGiCUMd9pwNTmxix7gTOh/lnXjgKO/2LBZhRl4an22rDOmPfV78YBoEyhrKgq5DqUjShKEOe6X1paTEZOAdtBfoklRptwsbzRqdScm+Rgorgh+PkZc8AoCiyo1B16FUJGnCUIc9Z72TOPHehwh8SergaUE6OuFI78JLn2wNfj6o3LSxABSWfRN0HUpFkiYMddhzNjixeFJJiLEyyGE/uEAPWhgj0h3kDIrn4y3BJ4zB6WNxeDwUVmkLQ0UnTRjqsOYftOdqSWJ4mgPpeJ/C1QKNVd22MESEE0Zn8MW2SlxuT1CxSOpw8lpbKarVlfdUdApLwhCROSKySUS2isidAfb/RETWi8gaEVkmIiPCcV6lQuUftNfQkNT55SjotoUBMOvITGqbXXxTHGTX2OSh5La6KGzUm94qOoWcMETECjwGnAmMBy4VkfEdiq0CCowx+cCrwEOhnlepcPAP2ttbHd/5tObQbQsD4LhR6YgQ/GUpq508iWWPu4GG1obg6lAqgsLRwpgObDXGbDfGtAAvAue2L2CM+Y8xxv8O+ALICcN5lQqZfwxGY1Ny59OaQ49aGIMSYsgflsInIdzHyIvxjvXYWbsz6DqUipRwJIxhQPuLrsW+bZ25Fngn0A4RuUFEVojIivLy8jCEplTX/KO8TWsKwwO2MHwJowctDIBZozNYtWsftU2tQcWTmzgUgMLqwqCOVyqSwpEwAo1mCjh6SUQuBwqA3wfab4xZYIwpMMYUZGZmhiE0pbrmbHBiFVMqTL0AAB5/SURBVDvG7eikheEEBBJ69nqcdWQmbo/hi+3BTRMyImUkYgyF+3QSQhV9wpEwioHh7Z7nALs7FhKR7wC/As4xxjSH4bxKhcxZ7yTBmgZYOrmHUQqJg8Fq61F900akEm+38smW4FrIsakjGOZyUbR3c1DHKxVJ4UgYy4HRIpInIjHAXGBx+wIiMhX4C95koV1AVNQorS/FbgaR6rCTFBdoDIYz4LTmnYm1WZkxMo2Pgx3AlzLM21Nq37bgjlcqgkJOGMYYF3ALsATYALxsjFknIveKyDm+Yr8HEoFXRGS1iCzupDql+pSzwYmnNTXwHFIA+3Z2Oq15Z2YdmcH28np27wtiuvPkHO9YjPrdeExw4zmUipSetbO7YYx5G3i7w7a72j3+TjjOo1Q4+Qft2ZvGMTHQHFIlX0P5Bjj6h72q94TRmcAGPtlSwcXHDO+2/AFSvAmjydOKs97JkMQhvTteqQjSkd7qsOUftFdT18m05sufAnsCTPlBr+odk5XI4KRYPgrmPkZCJrku70PtKaWijSYMddhyNngH5bmakw/uUltfCd++CpPnQlxKr+oVEWYdmcFn2yrxeHo53bnFQl6cdyLEwhpNGCq6aMJQhy3/GAyPK4Wcjl1qv34G3M0w/Yag6p41OoOq+hbW76np9bHpyTkkGdEWhoo6mjDUYcs/LYjpeNPb7YIVCyHvRBg8Nqi6Z/mmOw9mmhBJGU6e27Bl75agzq1UpGjCUIet0oZSLNgwbseBYzA2vwvVu4JuXQAMTo7jqKwkPtkaxH2M5GHMqq3h67Kv2a4D+FQU0YShDlvOeiexMojMpHji7Nb9O776CyTnwJgzQ6r/hNEZLC/aS1Oru3cHpgzjkpoa4qyxLFq3KKQYlAonTRjqsOVscGJxDzrwhnfZRij8CI65tsejuzsza3QGLS4PXxX2cpqQlOGkeTycl308b2x/Q9f4VlFDE4Y6bJXWl9LSnHTgHFLLnwRrLEy7MuT6Z+SlE2O19H7Z1mTv3J1XDpqEx3h4dsOzIceiVDhowlCHJY/xUNZQ5ls4ydfCaKqG1S/AxO9DQkbI54iPsXL0iEG9v/Gd4k0Yw5sbmT1iNq9seoW6lrqQ41GHtl21uyI+O4AmjAHIYzx8XPyxLrITgr1Ne2n1tOJuTd7fQ2r1C9BaDzOCv9nd0azRGWzYU0N5bS/m24xLhZhEqC7h6glXU9dax6ubXw1bTOrQ83Hxx5z12lnc9eld3RcOgSaMAcblcfGbT3/Dzctu5vr3rqe6ubq/QxqQ/AsnmdYU7yUpj8d7OSrnGBg6NWznOWG0t6Xy2bZetDJEvJelqncxIWMCM7Jn8Pf1f6fVHdwaG+rQ9+S3TxJvi+f1ba9H9DyaMAaQFncLd3x0B4u3Lea7I7/LhsoNXLPkGioag1/h7XDlH4PhcaV4Wxjb/wOVW0PqShvIhKEppDrsvb8slTYSdn4Oe3dw9cSrKWss463Ct8Iamzo0rHSuZFXZKm6bdhu3TbstoufShDFANLoamff+PN7b8R53HHMHD5zwAI+e9ii7andx1btXsaduT3+HGBS3x40xvZw+Iwz8o7xxpTAkNQ6+WgAJg2H8eWE9j9UiHD8qg0+2VPTu5/zOPd4BhM9dxHGDxjFm0BgWrV2kM9iqgzz97dOkxaVx/ujzuXbStRE9lyaMAaC2pZab3ruJz3Z/xj0z7+GK8VcAcNzQ4/jL7L9Q2VjJle9eyY6aHf0cac8ZY3h96+uc8NIJLFizoM/P72xwItjITszAXr0DNi+Bo68CW0zYzzVrdAalNU1sLevFjevBY2Huc1C1HXnpCq4adznbqrfxScknYY9PDVybqjbxccnHXDbuMuJtAWZcDjNNGFFuX9M+rvv3dawpX8NDJz7E98d8/4D9UwdPZeEZC2l2NfPDd37IpqpN/RRpz1U0VjDv/Xn8+tNf4/K4WLRuUZ/3AnI2OLF5Uhk+KAFWPA0WKxRcHZFzBT1NSN4JcN7jsOMT5qxeTLYjm4VrF0YgQjVQPf3t0yTYE5g7dm6fnE8TRhQrbyjn6iVXs3XvVv731P9lTt6cgOXGpY9j0ZxFWC1WrllyDWvK1/RxpD33buG7nPf6eXy2+zN+VvAzFp6xkLrWOl7Z/EqfxlFaX4q7NZmRqRb4+u8w7nuQPDQi5xqe5iA33dH78RgA+RfBaXdhX/caV9qzWOlcyTfl34Q/SDXg7KrZxZIdS7h4zMUkxyT3yTnDsoDSQGCM4bUtr/H8xue5efLNnDbitP4OqUsldSVc/+/rqWys5PHvPM70IdO7LD8ydSTPzHmG6/99Pdf/+3oePe1Rjsk+po+i7d7epr3M/3I+S4qWMDF9IvNnzWdk6kgAZgzx9gK6bNxlxFjDf0movd11u3n8m8dZVbaK5oYCTnV8CE37wn6zu6MTRmfyj6+LaXF5iLH18nParJ/A3h18f9XfeHzkkSxau4g/nvLHsMe4qmwVj656lHWV6xjsGEyWI4vshOy27+0fJ9oTEZGwx6B67q/r/opNbG2XqPvCYZEwGlobmP/lfBZvW0ySPYnbPriNC0ZfwM+P+TkOeydLc/aj7dXbuf7f19PkauLJ058kPzO/R8flJOXwzJnPcMO/b+C/lv4XD5/8MCfmnBjhaLv3/s73+e3nv6WmpYZ5U+dx9cSrsVn2v/SumXgNN753I29tf4vzR58fkRgqGit4cs2TvLz5ZSxYOHvExTy3cQzHeH4PWRPhiJkROa/frNEZ/P2LHazauZcZI9N7d7AIfPdhHDUlzK1cyVOeZRRVF5GbkhuW2NZXrufPq/7MJyWfkB6Xztkjz6aqqQpnvZPPdn9GRWPFQTfbHTYHkzIncc/Me8hJyglLHOHU5GpifeV64mxxJMUkkRyTTKI9EavF2v3BA0B5Qzn/2vovzjvyPDIdmX123kM+YRRWF/KTD37Ctn3buHnyzVwz6RoeX/04C9cuZHnpch444YEe/0HuCxsqN3DT0psQhIVnLOSotKN6dfxgx2D+Ouev3LT0Jn78/o+5/4T7O72UFWnVzdU8+NWDvLH9DcamjWXB7AUBf56ZQ2YyLm0cC9cu5Nwjz8Ui4btSWt1czTPrnuHZDc/S4m7hvCPP46bJN7G5xMom8zdSazbBSX/y/lGOoJmj0rFahE+2VvQ+YYB3XquLFvGDv87hGVPN35b/kbu+878hxbRt3zYeW/0Y7+14j+SYZG6bdhuXjr30oA9RrZ5WKhsrKa0vpbShFGe9kz31e1i8dTEXv3kx84+fzylHnBJSLOHi9rhZvG0xj61+rG2BrPYS7AkkxSR5v+zeRJKVkMW1E68dUMvh/n3D33EbN1dPiMx9t85If3Rp7ImCggKzYsWKkOp4t/Bd7v7sbmKtsTxwwgMcN+y4tn0rSlfwy09+SVlDGTdOvpHrJ11/wKfenmpobaDV04qI0PbP98dH8H0XwRjDvuZ9VDVVUdVURWVjpfd7k/d7VaN3+46aHaTGpfLk7CdD+gRZ21LLLctuYVXZKu6cfic/GNe7ZUaDZYyhsqmSlc6VPLT8ISobK7lu0nXcmH8jdqu90+PeKXyHOz66g0dOeYTTjgj9cmFDawPPb3yehWsXUttSy5m5Z/KjqT9iRPIIAJ77cgfJb97Adx0bsPx0I8REvqV5wf99isfAv350fPCV1Ozhty/MZnEsLJnzHBnZk3tdxa6aXTz+zeO8uf1N4m3xXDnhSq4cfyVJMUm9q6d2Fz/78Gesr1zPVROuYt60edgtnf8fd6XR1Yiz3smI5BFBXeoyxvBxycf8ceUf2bpvKxPTJ3L1xKuxWqzUtdRR21JLbUstNS01bY9rW73fi6qLEBFunnwzl42/LOifoa/UtNRw+qunc+KwE3nopIcO2i8iK40xBZE4d1gShojMAf4XsAJPGWMe6LA/FvgbcDRQCVxijCnqqs5QEkaLu4X/WfE/vLDxBSZnTuZ/TvofshOyDypX21LLf3/537y5/U3yM/N5YNYDDE8e3m39O2p2sHTHUpbuWMrayrVBxejnsDlIi0sjLT6NtLg0shzh+7TT6Grkjo/u4INdH3DJUZdw5/Q7g0qKgRhjcDY42bZvG9v2bWN79fa27zUt3lXmRqWMYv6s+UzImNBtfS6Pi7P/eTbp8ek8e+azQV8fb3G38MrmV3hyzZNUNlVyUs5J3Dr11oNaNo++/jE3fn0u1pn/hWXO/KDO1VsPv7eZR9/fwqrfnE6KI/g/SkWFyzjnwx9zXaudeZf/B+JTe3RcaX0pf1nzF/615V9YLVYuHXsp10y8hkFxg4KOpcXdwu+X/54XN73IlMwp/P6k3wd8r3WmobWBFze9yDPrnqGqqYojko5gTt4czsw9kyMHHdmjOtZVrOPhlQ/zVelXDE8azrxp8zhjxBk9fg3trtvN/V/ezwfFHzB60GjuOvYupgye0uOfoa8tWLOAP6/6M69+79WALfaoThgiYgU2A7OBYmA5cKkxZn27MjcD+caYm0RkLnC+MeaSruoNNmHsrtvNTz/4KWsr13LF+Cv4f0f/v24/MbxT+A6/+/x3uI2bO6ffyXlHnnfAi80Yw5Z9W1i2Yxnv7XyvbSW0CekTOCnnJBJjEjHGYDAHHNP2zxhEhNTYVNLj0g9IEJHuO+32uHnk60dYtG4Rxw09jt+f9Puge1RUNFawYM0C1lasZdu+bTS49s9llRaXRl5KHqNSRjEydSSjUkcxbfC0Xt3EfnHji8z/cj6L5izi6Kyjex1fQ2sDV7xzBZv3bqYgq4AfT/vxgW98VzM07oXGvXzw/IOcuO9fWOatgrS8Xp8rGMuLqrjoic954vJpzJkY2geC2968nOVlq3jPMgLH5f8EW2zAcsYYVpWtYvG2xbyx7Q08ePj+6O9zQ/4NDHYMDimG9vyt+RhrDPefcD+zhs3qsnx9az0vbHyBv637G3ub93Lc0OM4MedEPtj1AV+VfoXHeDgy9Ujm5M5hTt6ctpZhe7tqdvGnVX/i3aJ3GRQ7iJsm38RFYy7qsiXbGWMM7+96n/u/vB9ng5MLx1zIbdNuIyW2d+u5R1qjq5E5/5jD+PTxPP6dxwOWifaEMRO4xxhzhu/5LwCMMfe3K7PEV+ZzEbEBpUCm6eLkY4/INAt/3rsboN9a97Iwbgse4IdNo5jq7vm14ippZlHsVjbbapjqSuOyppFUWJpZZatkla2KMksTYmCUO4kprnSmutJIM4HfpKHp/v8jmP+yT+1lvBBXSKYnlpsbxpJp4npwJi8Pho/tThbH7cKFh1HuJLLdDoZ44sn2xJPliSfJdPMm7cHJmnFxd9JqRrgTuamhd/duAP4ev43l9gpuqD+CWU1CgrsWh7sWh7sGh7uWWNN0QPmVcTM5+s53e32eYLW6PUy99z1GZSYwbUTwn+oBKl1b+E/Nb/h55V5Ob8mg1p6BR2y4xYpb7Oyxevgopo5PYqopt7QQYyxMd2VwVmsOGSaul2fr2Sd1pzTyRNxGSiwNnNWSw/daj8Da4dhGXLxv38NS+27qxcUEVypntw5nlGf/h5gaWlhpq2SFrYItVm9r9Qh3Ase4MyhwZRBjrLwVs4sPbaVYEWa3DuX01mHE+2/JhnA7qgk3b9h3ssy2mwTsXNiSy7HuzLbLy10dV2SpZbullt2WBmxYsBsLMfi+jNX33b/NSqKxMdqT3G3d7b1v28MLMdu5vWkiYzyBk9nMm/8S1QnjQmCOMeY63/MrgBnGmFvalVnrK1Pse77NV6aiQ103ADcATB1iO/qjG7r/g18vQoXNwjuJsTyTmsCRzS7+u6yG4a5ernIGeIDnU+L5y6AEDOAWwWoMRze1cnJ9Myc2NJPujo57PsG8J1bG2fnF4GQEeKCshqlN3U9mtzHGxoMZiWyItXNMYwu3V9RxRBC/255amOpgwaAEniuuYlQvVqp7OzGWezOTuXZvPZdXe6ghiVpJpEYSD3pc53s8Y9ZsLjlhYsR+lkDuWbyOf3xdHJ7KhjyGw7aHR3ZZicNDq7SyPMHFR0mG9XGCGMOUJjdn1DVzWn0jScFMK9LLl3uTwB/Sk3gjOZ5pjS3cW1ZDhttDnQgvp8TzYoqDGquF4+ubuWZfPROaXV3W57RaWJYYy9KEONbHeT+UxHoMrQLfq23iur31ZLrDP13KlhgbD2YksTbOzrTGFu6oqCXX93o0wC67lW9j7ayNs7E21s62GBse31WJbF+5ZovQJNAkgunk8thxDc38pqyGQZ7uf9Eu4MLh6Qx2u1mwe1+n5RLvLYvqhHERcEaHhDHdGHNruzLrfGXaJ4zpxpjKzuotKCgw73z0DkU1RZQ3lFPeWE55QzlljWVtz8saymh0NbYdc8HoC/jF9F8QZ+vtJ6gDbazayKubX2VSxiROHn5y1DVLQ7GzZic/WvYjiuuKuXvm3Zx3ZOC5k2pbavnzqj/z0qaXSItL445j7mBO7pyI972vbq5m9quzmT1iNvNn9ezewvbq7cx9cy4T0ifw1OlPHTJdJ7vz4a4PueX9W7h6wtWU1pfy/q73aXY3k5ucyzmjzuHskWf3W8+fxdsWc98X9+GwOfjuyO/yz63/pLallpNzTuamyTf16L5WR7tqd7GkaAm763Zz+bjL28bxRIrHeHh186s88vUjNLoaOXfUuZQ1lLGmYk3bLNGJ9kQmZUxi8uDJTM6czKSMSQf9vTDG0OJpocnV5P1ye78vL13OH1f+kZTYFO4/4X5mDJnRZTyLty3mV5/8isdOe6zL7vKH5SWpxJGJJu/uA68tx1njyHRkkhmf2fZ9sGMwmY5MRiSNYFLmpJB+lsNFdXM1P/vwZ3yx5wuunng1t027ra0rqzGGd4vebevhNHfsXG6demuve9CE4sGvHuTFjS/y9gVvd/sHr8nVxGVvX0Z5QzmvfO8VshKy+ijK/ucxHi54/QK2VW8jOSaZM/PO5JxR5zApY1JUDKrbuncrP/3wp2yv3s6pw0/lpsk3MS59XH+H1WsVjRX8YcUfeLvwbfKS88jPzGdypjdBjEwdGVI38E1Vm7j9o9spqi7i2knXcvOUmwPec/X/X1ssFv7xvX90+f8byYThvTkbwhfesRzbgTwgBvgGmNChzI+AJ3yP5wIvd1fv0LFDzXPrnzOflnxqtu7daqqbq43H4zEqPFrcLeZ3n//OTFw00dy67FZT31JviqqLzHVLrjMTF000F79xsVlbvrZfYttdu9tMeWaKeeDLB7ot6/8ZPtz1YR9EFn2279tu/rPzP6bZ1dzfoQTU2Npodtbs7O8wwsLldkWk3vqWenP3p3ebiYsmmh+89QOzq2bXQWWW7VhmJi6aaN7c9ma39QErTIh/1zv7Cle32rOAR/B2q11ojJkvIvf6Al8sInHA34GpQBUw1xizvas6wzEOQ3XNGMPzG5/noeUPkZOYQ2l9KTHWGOZNm8fFYy7u10s7v/z4lyzduZT3Lnyv00uCS4qW8LMPf8ZVE67ipwU/7eMIlQqvdwvf5bef/xaAu2fe3Tbg1hjD5W9fTmVTJW+e/2a3XeOj+pJUpGjC6DsfF3/MLz/5JTOHzOT2Y27v06kGOrN572a+v/j73DLlFm6cfONB+3fV7uLiNy5mZMpIFs1ZFFRXSqWiTXFtMT//+OesKV/TNn3Rusp1XLPkGn4949dcMrbL0QiAJgzVB4xvrEg0uXnpzayrXMeS7y85oCNDq7uVK9/xrv/xyjmvMCxxWD9GqVR4tXpaeXz14zz17VPkpuSSaE9kd91ully4hFhr9135I5kwdHpzBRB1yQK8kxJWNVXxr63/OmD7I18/wtrKtdx7/L2aLNQhx26xM2/aPBacvoC6ljq+rfiWK8Zf0aNkEWmaMFTUOjrraPIz83lm3TO4PN7++h/u+pC/rf8bc4+ay3dGfKefI1Qqco4dciyvnvMqd06/k8vGXdbf4QCaMFQUExGumXANxXXFLN2xlNL6Un716a8YmzaWnx3zs/4OT6mIS4tL47Jxl4U8tixcDvnpzdXAdsoRp5CbnMvCtQuJt8V7J7s78fdR0TxX6nCjLQwV1Sxi4eqJV7OhagNfl33NXTPvCtvCQUqp3tGEoaLe2SPPZlTKKOYeNZezR57d3+EoddjSS1Iq6sVYY3jt3NfCuhKfUqr39B2oBgRNFkr1P30XKqWU6hFNGEoppXpEE4ZSSqke0YShlFKqRzRhKKWU6hFNGEoppXpEE4ZSSqke0YShlFKqRzRhKKWU6hFNGEoppXpEE4ZSSqkeCSlhiEiaiLwnIlt83wcFKDNFRD4XkXUiskZEul/FXCmlVNQJtYVxJ7DMGDMaWOZ73lEDcKUxZgIwB3hERFJDPK9SSqk+FmrCOBd4xvf4GeC8jgWMMZuNMVt8j3cDZUBmiOdVSinVx0JNGFnGmD0Avu+DuyosItOBGGBbJ/tvEJEVIrKivLw8xNCUUkqFU7cLKInIUiA7wK5f9eZEIjIE+DvwQ2OMJ1AZY8wCYAFAQUGB6U39SimlIqvbhGGM+U5n+0TEKSJDjDF7fAmhrJNyycBbwK+NMV8EHa1SSql+E+olqcXAD32Pfwi83rGAiMQA/wT+Zox5JcTzKaWU6idiTPBXfkQkHXgZOALYCVxkjKkSkQLgJmPMdSJyOfBXYF27Q68yxqzupu5aYFPQwfW/DKCiv4MIgcbfvzT+/jOQYwc4yhiTFImKQ0oYkSQiK4wxBf0dR7A0/v6l8fevgRz/QI4dIhu/jvRWSinVI5owlFJK9Ug0J4wF/R1AiDT+/qXx96+BHP9Ajh0iGH/U3sNQSikVXaK5haGUUiqKaMJQSinVIxFPGCKyUETKRGRtu233iEiJiKz2fZ3l2z5bRFaKyLe+76e2O+Zo3/atIvInERHf9m6nWI+S+OeLyC4RqetQf6yIvOT7ub4Ukdxoil1EHCLylohs9E1R/0BfxB6u+H373hWRb3zxPyEiVt/2AfHaaXfs4g51DYj4ReQDEdnU7pjBvu0D5fUTIyILRGSz733w/UjHH6b3blK7sqtFpEJEHgkpdmNMRL+AE4FpwNp22+4Bfhag7FRgqO/xRKCk3b6vgJmAAO8AZ/q2PwTc6Xt8J/BglMZ/LDAEqOtwzM3AE77Hc4GXoil2wAGc4nscA3zc7ncfsdjD/LtP9n0X4B/A3IH02vFtuwB4vkNdAyJ+4AOgIMAxA+X181vgPt9jC5AR6fjD+dppV24lcGIosUe8hWGM+Qio6mHZVcY7BTp4R4bH+TLhELxv+s+N9yf8G/unUu92ivVQhCN+374vjG9m3w7ax/8qcJqIt/UUqnDEboxpMMb8x1emBfgayIl07OGK37evxrfdhjfp+Xt6DIjXjogkAj8B7utw2ICIvwsD4vUDXAPc7yvnMcb4R4FH9Xu3fRkRGY13NvGPfZuCir0/72HcIt4V+BZ20pT+PrDKGNMMDAOK2+0r9m2DXk6xHka9ib8rw4BdAMYYF1ANpIc31IMEFbt4F776Ht7FsqB/Yocg4heRJXgnx6zF+waBgfPa+R3wB7yLkbU3UOIH+Kvvsshv2v1hivrXj+xf7O13IvK1iLwiIlm+bQPmvQtcircV4f+wFFTs/ZUwHgdGAVOAPXjfDG1EZALwIHCjf1OAOvqzP3Bv4+9KX/9sQcUuIjbgBeBPxpjt/s0B6o/0/0tQ8RtjzsB7STAWOOj+QB/qVfwiMgU40hjzzz6OszPB/P4vM8ZMAk7wfV3hLx6g/mh7/djwtqg/NcZMAz4H/sdfPED9Uffe9ZmL9/3bVjxAme5jD9c1t26ux+XS7lpcV/vw/udsBo5vt20IsLHd80uBv/gebwKGtCu3Kdri71C+4z2MJcBM32Mb3knPJNpiBxbiTRZ9Fnu4f/e+Mj8EHh0orx3gv4DdQBHelnUL8MFAiT/AMVe1+/1H/esH7x/WesDiez4cWNcX8YfxvTsZ2NxhW1Cx90sLw3dPwu98YK1veyredTN+YYz51F/AeJvbtSJyrK85eyX7p1Lvdor1cOtt/N1oH/+FwPvG978YCcHELiL3ASnAbR2q69PYfbH0Kn4RSfQf42slnQVsDBB/VL52jDGPG2OGGmNygVl43/gn+3ZHffwiYhORDN9jO3C2/xgGwOvHF88bwMm+TacB632Po/6963MpB7YuINjYw5nNO8mEL+BtPrXi/YR0Ld6V974F1vgC939K+jXebL663ddg374C3y9oG/Ao+0epp+O9pr7F9z0tSuN/yHe8x/f9Ht/2OOAVYCvenmAjoyl2vJ9cDLCh3fbrIh17GOPPApb7yq8D/gzYBtJrp119uRz4qTLq4wcS8PbO8f/+/xewDpTXj2/fCOAj3zHLgCMGwnu3XV3bgbEd6g8qdp0aRCmlVI/oSG+llFI9oglDKaVUj2jCUEop1SOaMJRSSvWIJgyllFI9oglDKaVUj2jCUEop1SP/H0Y2741PtVZVAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "halpha_cm1 = Lines().get_line_cm1('Halpha')\n", "\n", "step = 2943\n", "order = 8\n", "step_nb = 840\n", "axis_corr = 1.0374712062298759\n", "theta = orb.utils.spectrum.corr2theta(axis_corr)\n", "print('incident angle theta (in degrees):', theta)\n", "zpd_index = 168\n", "\n", "# model spectrum\n", "velocity1 = 50\n", "broadening1 = 15\n", "spectrum_axis, spectrum1 = orb.fit.create_cm1_lines_model_raw([halpha_cm1], [1], step, order, step_nb, axis_corr, zpd_index=zpd_index, fmodel='sincgauss',\n", " sigma=broadening1, vel=velocity1)\n", "\n", "velocity2 = 10\n", "broadening2 = 30\n", "spectrum_axis, spectrum2 = orb.fit.create_cm1_lines_model_raw([halpha_cm1], [1], step, order, step_nb, axis_corr, zpd_index=zpd_index, fmodel='sincgauss',\n", " sigma=broadening2, vel=velocity2)\n", "\n", "spectrum = spectrum1 + spectrum2\n", "\n", "# add noise\n", "SNR = 22\n", "spectrum += np.random.standard_normal(spectrum.shape) * 1. / SNR\n", "\n", "spectrum_axis = orb.utils.spectrum.create_cm1_axis(np.size(spectrum), step, order, corr=axis_corr)\n", "\n", "pl.plot(spectrum_axis, spectrum1, label='line 1')\n", "pl.plot(spectrum_axis, spectrum2, label='line 2')\n", "pl.plot(spectrum_axis, spectrum, label='line1 + line2 + noise')\n", "pl.xlim((15200, 15270))\n", "pl.legend()\n", "pl.savefig('gvar_model.svg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Classical fit\n", "\n", "The classical fit will be be unable to make any difference between an infinity of different possibilities which all gives approximatly the same chi2. the best fit will be very badly constrained and can give random sets of parameters depending on the noise." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING:root:nan in passed parameters: {'amp0': 1.18644 +- nan, 'amp1': 0.0516463 +- nan, 'pos_def1': 26.9033 +- nan, 'pos_def2': -57.1967 +- nan, 'sigma0': 34.6287 +- nan, 'sigma1': 0.5303 +- nan, 'cont_p0': 0.00261949 +- nan}\n", "WARNING:root:nan in passed parameters: {'amp0': 1.18644 +- nan, 'amp1': 0.0516463 +- nan, 'pos_def1': 26.9033 +- nan, 'pos_def2': -57.1967 +- nan, 'sigma0': 34.6287 +- nan, 'sigma1': 0.5303 +- nan, 'cont_p0': 0.00261949 +- nan}\n", "WARNING:root:nan in passed parameters: {'amp0': 1.18644 +- nan, 'amp1': 0.0516463 +- nan, 'pos_def1': 26.9033 +- nan, 'pos_def2': -57.1967 +- nan, 'sigma0': 34.6287 +- nan, 'sigma1': 0.5303 +- nan, 'cont_p0': 0.00261949 +- nan}\n", "WARNING:root:nan in passed parameters: {'amp0': 1.18644 +- nan, 'amp1': 0.0516463 +- nan, 'pos_def1': 26.9033 +- nan, 'pos_def2': -57.1967 +- nan, 'sigma0': 34.6287 +- nan, 'sigma1': 0.5303 +- nan, 'cont_p0': 0.00261949 +- nan}\n", "WARNING:root:nan in passed parameters: {'amp0': 1.18644 +- nan, 'amp1': 0.0516463 +- nan, 'pos_def1': 26.9033 +- nan, 'pos_def2': -57.1967 +- nan, 'sigma0': 34.6287 +- nan, 'sigma1': 0.5303 +- nan, 'cont_p0': 0.00261949 +- nan}\n", "WARNING:root:nan in passed parameters: {'amp0': 1.18644 +- nan, 'amp1': 0.0516463 +- nan, 'pos_def1': 26.9033 +- nan, 'pos_def2': -57.1967 +- nan, 'sigma0': 34.6287 +- nan, 'sigma1': 0.5303 +- nan, 'cont_p0': 0.00261949 +- nan}\n", "WARNING:root:nan in passed parameters: {'amp0': 1.18644 +- nan, 'amp1': 0.0516463 +- nan, 'pos_def1': 26.9033 +- nan, 'pos_def2': -57.1967 +- nan, 'sigma0': 34.6287 +- nan, 'sigma1': 0.5303 +- nan, 'cont_p0': 0.00261949 +- nan}\n", "WARNING:root:Nan in model\n", "WARNING:root:Nan in model\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "velocity (in km/s): [26.9033 +- nan -57.1967 +- nan]\n", "broadening (in km/s): [34.6287 +- nan 0.5303 +- nan]\n", "flux (in the unit of the spectrum amplitude / unit of the axis fwhm): [2.33652 +- nan 0.0583332 +- nan]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEICAYAAABMGMOEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXiU1dn48e/JnpCVLAQIJCBhJ4CAgNaVgoiWultqq/wuFZdSbfvWamtrsXVvX+3rUtEqxbpURC2yKoUWEcElKFtYQ0xIIJB1Jvus9++PGWL2dbLJ/bmuXGSe88w5JyEz9zzPOec+RkRQSimlWuPX0x1QSinVN2jAUEop1SYaMJRSSrWJBgyllFJtogFDKaVUm2jAUEop1SYaMJTqBGPMRcaYPB/Wl22M+W4zZaHGmDXGGKsxZqUx5kZjzEZfta1UawJ6ugNKtcQYswWYCCSKiK2Hu9PTrgUGALEi4vQee+N0oTFGgFQRyeyJzqlvP73CUL2WMSYFOB8QYH43tdmbP0QlA4frBAulupUGDNWb3QR8CiwHbm7uJGPMD4wx6Q2O/dwYs9r7fbAx5s/GmGPGmFPGmKXGmFBv2UXGmDxjzH3GmJPA340x+4wx36tTV6AxpsgYM6mFPvzGe062MebGOscvN8Z8ZYwpM8bkGmOWNHjej40xOcaYYmPMAy3U/xDwIHCDMabCGHOLMWahMWabt3yr99Td3vIbmqtLqY7SgKF6s5vw3HJ5A7jUGDOgmfNWA6OMMal1jv0QeNP7/RPASGASMAIYjOfN97REoD+eT/CLgH8AP6pTPg/IF5FdzbSfCMR5670ZeMkYM8pbVun9OaKBy4E7jTFXAhhjxgIvAD8GBgGxQFJTDYjI74FHgRUiEi4irzQov8D77URv+Ypm+qpUh2nAUL2SMeY7eN7A3xaRncBRPEGgERGpAt4HFnifmwqMBlYbYwxwG/BzESkRkXI8b7w/qFOFG/i9iNhEpBp4HZhnjIn0lv8YeK2VLv/O+/yPgHXA9d6+bRGRvSLiFpE9wD+BC73PuRZYKyJbveMzv/P2RaleSQOG6q1uBjaKSJH38Zu0cFvKW77A+/0PgVXeQBIPhAE7jTEWY4wF+MB7/LRCEak5/UBETgCfANcYY6KBy6gzuNyEUhGprPM4B88VA8aY6caY/xpjCo0xVuAOPFcjeM/JrdNuJVDcQjtK9ajePMCnzlDe8YXrAX/vuAJAMBBtjJkoIrubeNpGIM47zrAA+Ln3eBFQDYwTkePNNNlUyuZXgVvxvEZ2tPBcgBhjTL86QWMosM/7/ZvAc8BlIlJjjPkL3wSMfGBMnZ87DM9tKaV6Jb3CUL3RlYALGItn3GESnjfWj/GMBzTinTn0DvAnPOMR//YedwN/A542xiQAGGMGG2MubaUPq4CzgXvwjGm05iFjTJAx5nzgCmCl93gEUOINFudQ/7baO8AVxpjvGGOCgD/QudfkKWB4J56vVIs0YKje6Gbg7yJyTEROnv7C80n9xhamvr4JfBdY2WDq6X1AJvCpMaYM2ASMauL5tbxjGe8Cw4D3WunvSaAUOIHn1tUdInLQW3YX8AdjTDmegfa367SRAfzE2+98bx2dWQS4BHjVe+vt+k7Uo1STjG6gpFTTjDEPAiNF5EetnqzUGUDHMJRqgjGmP3ALnhlSSin0lpRSjRhjbsMze2mDiGxt7XylzhR6S0oppVSb6BWGUkqpNum1YxhxcXGSkpLS091QSqk+ZefOnUUiEt/6me3XawNGSkoK6enprZ+olFKqljEmp6vq1ltSSiml2kQDhlJKqTbRgKGUUqpNeu0YhlJnGofDQV5eHjU1Na2frM54ISEhJCUlERgY2G1tasBQqpfIy8sjIiKClJQUPNt4KNU0EaG4uJi8vDyGDRvWbe3qLSmleomamhpiY2M1WKhWGWOIjY3t9qtRDRhK9SIaLFRb9cTfik8ChjFmmTGmwBizr5nyG40xe7xf240xE33RrlJKqe7jqyuM5cDcFsq/Bi4UkTTgj8BLPmpXKeVD4eHhAJw4cYJrr722w/U899xzjBgxAmMMRUVFrT+hDS666KLaxbzz5s3DYrF0qJ6nnnqKsWPHkpaWxqxZs8jJ6bJ1bs3q7O+3p/gkYHgzepa0UL5dREq9Dz8FknzRrlK9gYjwyiuv8MILL/R0V3xm0KBBvPPOOx1+/nnnncemTZtITk5u0/lLlixh+fLlba5//fr1REdHd6hvkydPJj09nT179nDttdfyq1/9qsXzt2zZwsKFCzvUVnM6+/vtKT0xhnELsKGpAmPMImNMujEmvbCwsJu7pVTHVFVVkZeXR0FBQU93xWeys7MZP348AMuXL+fqq69m7ty5pKam1nuD3bhxIzNnzuTss8/muuuuo6KiAvC8KXdlLriUlBSKiorIzs5mzJgx3HbbbYwbN445c+ZQXV0NwNGjR5k7dy5Tpkzh/PPP5+BBzyaIF198MWFhYQDMmDGDvLzObHLosXDhQu6++27OPfdchg8fXhsMRIR7772X8ePHM2HCBFasWAHU//1mZGRwzjnnMGnSJNLS0jhy5AgAr7/+eu3x22+/HZfL1el+dla3BgxjzMV4AsZ9TZWLyEsiMlVEpsbHd0nuLKV8LiQkhF/k5HB/VZVP612+fDm7du0CwOVysXz5cvbs2QN41mwsX76cffs8w4Y1NTUsX76cAwcOAJ4gtnz5cg4dOgRQ+0beUbt27WLFihXs3buXFStWkJubS1FREQ8//DCbNm3iyy+/ZOrUqTz11FOdaqcjjhw5wk9+8hMyMjKIjo7m3XffBWDRokU8++yz7Ny5kz//+c/cddddjZ77yiuvcNlll/mkH/n5+Wzbto21a9dy//33A/Dee++xa9cudu/ezaZNm7j33nvJz8+v97ylS5dyzz33sGvXLtLT00lKSuLAgQOsWLGCTz75hF27duHv788bb7zhk352RretwzDGpAEvA5eJSHF3tatUV/P39ydwzRrsMTEEP/FET3enS8yaNYuoqCgAxo4dS05ODhaLhf3793PeeecBYLfbmTlzZpvr3Lt3Lz/+sWdDw5MnTxIUFMRf/vIXADZv3kxsbGyb6hk2bBiTJk0CYMqUKWRnZ1NRUcH27du57rrras+z2Wz1nvf666+Tnp7ORx991GS906dPx2azUVFRQUlJSW0bTzzxBJdeemmj86+88kr8/PwYO3Ysp06dAmDbtm0sWLAAf39/BgwYwIUXXsgXX3xBWlpa7fNmzpzJI488Ql5eHldffTWpqals3ryZnTt3Mm3aNACqq6tJSEho0++jK3VLwDDGDAXeA34sIoe7o02lusuJo0cZVFSEu7KSzMxMRowY4ZN669439/f3r/c4MDCw3uOQkJB6j8PCwuo9Pj2Y3VHBwcH1+uJ0OhERZs+ezT//+c8O1TlhwoTaK6glS5aQkpLSobGChn2rrq7G7XYTHR1dW39DmzZt4pFHHuGjjz6q9/y6PvvsM8AzhrF8+fJWx1jq1nN6Y7q2bFD3wx/+kOnTp7Nu3TouvfRSXn75ZUSEm2++mccee6zV53cnX02r/SewAxhljMkzxtxijLnDGHOH95QHgVjgr8aYXcYYzVuuvjUObvAMyYVVVxPscPRwb7rPjBkz+OSTT8jMzAQ8t8EOH+4dnwcjIyMZNmwYK1euBDxv3Lt37wbgq6++4vbbb2f16tVd/qn9ggsuYMWKFbhcLgoLC9m6dSvnnHNOvXOysrIYPnw4d999N/Pnz2fPnj3MmjWLd955p3ZcrKSkpEdmczXkq1lSC0RkoIgEikiSiLwiIktFZKm3/FYRiRGRSd6vqb5oV6neYIJ3ABUgwTvgeiaIj49n+fLlLFiwgLS0NGbMmFE7sPzMM8+QlJREXl4eaWlp3Hrrrd3evzfeeINXXnmFiRMnMm7cON5//30A7r33XioqKrjuuuuYNGkS8+fP77I+XHXVVaSlpTFx4kQuueQSnnzySRITE+uds2LFCsaPH8+kSZM4ePAgN910E2PHjuXhhx9mzpw5pKWlMXv27EZjHz2h1+7pPXXqVNENlFRfcOrXv2bA448DUPr668TceGOH6jlw4ABjxozxZdfUt1xTfzPGmJ1d9aFcU4Mo1Qlutxvrrl24vWkaMtav7+EeKdV1NGAo1Qnl5eVw9CiFCQm4/PzoX17e011SqstoenOlOiEyMhIpL6dk8GCC7HaiNWCobzG9wlCqE4zbTXhhIY6hQ6mKiyPwxIk2TaVUqi/SgKFUJ2R9/DEBLhcMH05lXBzB+fm1i7aU+rbRgKFUJ2Rv3gxA4OjRuJKSiKioIFj3tFDfUhowlOqEsd7VvaHjxxNw1lkYEcJKmk3c3Ov5Kr15XStXrmTcuHH4+fmhU+X7Ng0YSnWC6/BhXH5+RI0fT5A3JUj5vib3EetTfJl+e/z48bz33ntccMEFPqlP9RwNGEp1kNPppCYjA2tMDKHh4YSMGgXAwQ8/7OGedV5n05vXNWbMGEZ5fzeqb9OAoVQHlZeXE5qfT1l8PMYYIsaOBSDWV2nOL7oITie8czg8j19/3fO4qsrz2Lu/Alar5/F773keFxV5Hq9Z43l88mSnutKb05ur7qPrMJTqoJiYGEItFvImTAAgNCaGivBwIoq/fdn7uyK9uep7NGAo1UGuwkJCqquR4cMBMMZQEReH//HjiAims7Oltmz55vvAwPqPw8LqP46Kqv84Lq7+4wYJ79qrK9Kbq75Hb0kp1UEZq1cD4D9yZO2xiv79CcrPp6ampqe61W16c3pz1TU0YCjVQeXezXlCvYPDAK7Bg4myWs+IF1ZL6c3r+te//kVSUhI7duzg8ssvb3K3OtU3aHpzpToo9847GbJ0KRWnThHu3Ygn6xe/YPjTT1OTlUXIsGHtqk/Tm6v20vTmSvUR5uhRyiMi6BcfX3ss4KyzAChtZmtQpfoyDRhKdYDNZoOjR2un1J4W7B3PyN66tae6plSX0YChVAdUVlYSWVRE5YAB9Y5HeMczdF8M9W3kk4BhjFlmjCkwxjSZE8F4PGOMyTTG7DHGnO2LdpXqKdEhIUSWleGfmlrveGhiIjUhIYRoxlr1LeSrK4zlwNwWyi8DUr1fi4AXfNSuUj2iYu9eAPwaBAxjDBWxsfjl5vZEt5TqUj4JGCKyFWgpRef3gX+Ix6dAtDFmoC/aVqonHPLu3R0yblyjsrLoaILy83UjJfWt011jGIOBuh+58rzHlOqbjh4FIHzixEZF7qQkoqzWPhkwuiK9+b333svo0aNJS0vjqquuwmKx+KRe1f26K2A0lSOh0avJGLPIGJNujEkvLCzshm4p1TEJFRXYgoMJT0lpVBZw1lmE2GzY+vA4hi/Tm8+ePZt9+/axZ88eRo4cyWOPPeaTelX3666AkQcMqfM4CTjR8CQReUlEporI1Pg6c9uV6m0Cc3I8U2r9Gr+E/L25pYq/+qq7u+UzvkxvPmfOHAICPGnrZsyYQV5eXvf8EMrnuitgrAZu8s6WmgFYRSS/m9pWyqeqqqoIzsujosGU2tNk6FAA8j/9tFPt9KLs5j5Lb75s2TIuu+yyznVG9RifZKs1xvwTuAiIM8bkAb8HAgFEZCmwHpgHZAJVwP/zRbtK9YSaykqiSkspPPfcJstjz/bMGo/6Ft2r90V680ceeYSAgABuvPHGbumz8j2fBAwRWdBKuQA/8UVbSvW0gFOnCHC5mpwhBRCWkoIzIICAE43uurZLL8pu3un05q+++ipr165l8+bNnU/7rnqMrvRWqp0qvWswgptJFGj8/amIicHk5HRnt7pdW9Obf/DBBzzxxBOsXr2asLCw7u6m8iENGEq1U85//gNAv7S0Zs+xRkcT1NmBg16urenNFy9eTHl5ObNnz2bSpEnccccdPdBb5Qua3lypdtpzxRWM27ABU12NX1BQk+d8PWsW8Z9/Tr+ysjbfgtH05qq9NL25Ur1cdEkJ5bGxzQYLgIDhwwmvqKDmWzTwrZQGDKXaKeT4caoHtpzZ5vTU2qI+vBZDqYY0YCjVDhaLhYiCgmbXYJxm805LKtq5szu6pVS30IChVDtU5eURWlODu5XtVweccw4A/YqKuqNbSnULDRhKtYPJygIgcvLkFs+LGD0aMUbTnKtvFQ0YSrVDzf79AIR58yw1xwQHUxkZiWRnd0OvlOoeGjCUaofCHTsAiGgirXlDlj64FuOZZ55hzJgxxMTE8PjjjwOwatUq9nsDpTqz+SQ1iFJnisjiYioiIgiPiGj1XBk6lMh9+xCRPpMO469//SsbNmxgWJ0xmlWrVnHFFVcwduzYHuyZ6g30CkOpdogsLKSqlSm1p/mlpBBptVLdRMrv3uiOO+4gKyuL+fPn8/TTT7N48WK2b9/O6tWruffee5k0aRJHvRtHqTOTXmEo1UZut5vwU6com9q2RbSOQYPwd7sp3b+fsOnT29XWBx98wEkf385KTExk7ty5zZYvXbqUDz74gP/+97+sXbsWgHPPPZf58+dzxRVX+GwHPtV36RWGUm10IiuLCKuV8oSERmXXXANBQbB9+zfHKmJjASjdtau7uqhUl9IrDKXayHbwIAYIapC7x+WCTZs8mxzV3Uwu+YILAAjsQJrzlq4ElOopGjCUaiPz9dcAxHkX5Z3m7+/Z8U4E6o5tR5zeL6OPpzmPiIigvLy8p7uhegG9JaVUGzkOHAAgvJkptaeDRXW190B4ONVhYbj7+FqMH/zgB/zpT39i8uTJOuh9htMrDKXaqGrvXmxBQQQ3yCO1Ywe8/DL88Y9w660QHAz/+penzBIV1afWYmR7g9vChQtZuHAhAOedd56uw1CAXmEo1WYxJSWUx8fXv+8EHD8OH3zg+X7ePLj00jqFyclElpbSW/edUao99ApDqTaKKCzEnpra6Pi113q+ABYvblA4dChRO3dSXVVFWL9+Xd9JpbqQT64wjDFzjTGHjDGZxpj7mygfaoz5rzHmK2PMHmPMPF+0q1R3qSgrI6qkBHdKSqvn2u3w0Uee76sSEghyOLB6B8yV6ss6HTCMMf7A88BlwFhggTGmYQ6B3wJvi8hk4AfAXzvbrlLdKffTTwlwubDGxTUq+/3v4aGHvnn87LNw0UWQmQnWqCgAKjIyuqmnSnUdX1xhnANkikiWiNiBt4DvNzhHgEjv91FA+yemK9WDTqc1j54ypVFZdjbUvYD44Q9hzRpITobRc+YA4PI+X6m+zBdjGIOBukn/84CGeRCWABuNMT8F+gHf9UG7SnWb0wEjYebMRmWvvlr/8cCBcMUVnu8DvYv83HpLSn0L+OIKo6k0nA2nhCwAlotIEjAPeM0Y06htY8wiY0y6MSa9sLDQB11TyjfcR47g8vMjoJWd9k4rL4enn4bPsyJxBAb2+bUYHZGSkkJRKzsONnfOAw88wJAhQwgPD++q7qkO8EXAyAOG1HmcRONbTrcAbwOIyA4gBGh0M1hEXhKRqSIyNT4+3gddU8pHsrKwxsRAQP2L8sJCmDMHNm9u/JQHHoCN/w7wrMXIz++mjn47fO973+Pzzz/v6W6oBnwRML4AUo0xw4wxQXgGtVc3OOcYMAvAGDMGT8DQSwjVZ0QXF1PTRFrzigpPWhCns/7xiAg4ehR++1t/SE4moqSk16/FyM7OZvTo0dx6662MHz+eG2+8kU2bNnHeeeeRmppa+wZeUlLClVdeSVpaGjNmzGDPnj0AFBcXM2fOHCZPnsztt99e7+d9/fXXOeecc5g0aRK33347Lperxb7MmDGDgW1MI6+6T6fHMETEaYxZDHwI+APLRCTDGPMHIF1EVgP/A/zNGPNzPLerFkpvf/Uo5VVTU0N0cTGWadMalQ0bBp991vTzTr/fuZOSiDxwgKqqKvq1dS3Gz34Gvs5yO2kS/OUvLZ6SmZnJypUreemll5g2bRpvvvkm27ZtY/Xq1Tz66KOsWrWK3//+90yePJlVq1bxn//8h5tuuoldu3bx0EMP8Z3vfIcHH3yQdevW8dJLLwFw4MABVqxYwSeffEJgYCB33XUXb7zxBjfddJNvfz7V5XyycE9E1gPrGxx7sM73+4HzfNGWUt2t6MgRkmpqsJx1Vruf+9BDMD53MNdUVXHixAn6NbHwrzcZNmwYEyZMAGDcuHHMmjULYwwTJkyoTRuybds23n33XQAuueQSiouLsVqtbN26lffeew+Ayy+/nJiYGAA2b97Mzp07meYNuNXV1SQ0kSJe9X660lupVhTs2EESUJ6QQMObJH/7G7z9tic1iL9/4+dWV0OODAag8sABaGvAaOVKoKsEBwfXfu/n51f72M/PD6f3vltTNwdOb0Hb1Fa0IsLNN9/MY4891hVdVt1Ic0kp1YqgXM+s8QHnNb5IPp3SvKlgAfD44/CTJz3p0O1HjnRZH7vTBRdcwBtvvAHAli1biIuLIzIyst7xDRs2UFpaCsCsWbN45513KCgoADxjIDl9POX7mUoDhlKt8POuoYiaNKlR2aJFsHFjy88PHjkSANe3ZGrtkiVLSE9PJy0tjfvvv59XvQtRfv/737N161bOPvtsNm7cyNChQwEYO3YsDz/8MHPmzCEtLY3Zs2eT38qssV/96lckJSVRVVVFUlISS5Ys6eofS7WB6a1jz1OnTpX09PSe7oZSHDzvPIbu20eY1dqh57/yooWFd8aya85lTPlgbbPnHThwgDENdvNTqiVN/c0YY3aKSNs2nm8nvcJQqhWhJ05Q0r9/k2Vz5nhuO7Vk2kywhEcTWtB39sVQqikaMJRqgd1uJ7qkxDN/tglxcRAZ2WRRrbS0aEzKYMJLinr9WgylWqKzpJRqQWl+PgllZdiauVX05pttq8c1OImonE/Yv7+KceN0XwzVN+kVhlItKP3yS0+ytA6swajrZFAYEeXl/P0lTXCg+i4NGEq1oMy72rqqiTQVH30EY8aANzNGi2oG9MdPhB9dfNjXXVSq22jAUKoFoSc8eTQHfuc7jctCYfx48C5obtHE+fMB8D/+7ViLoc5MGjCUaoF/Tg6O4GCCk5IalZ1zDqxcCUOGNPHEBoK8K7xPfp7FU0/5upe9U0fTm1dVVXH55ZczevRoxo0bx/33N9r1WfUQDRhKtSA4L4/yhATPcu5OyPU+33Ekp80D5WeyX/7ylxw8eJCvvvqKTz75hA0bNvR0lxQaMJRqlogQUVREcXR0k+W33AJz57axstBQKsPDGSon+OIL3/XRl3pLevOwsDAuvvhiAIKCgjj77LPJy8vrwp9ctZVOq1WqGZWVlURZrdgvuKDJ8kmToIk7VU0aMmQI1kGDCCsswJPhv+Urlh7Kbt7r0ptbLBbWrFnDPffc44tfgeokDRhKNaM8L4+BdjtBw4c3Wf7Tn7avPtfgwUTu2cMvf1nJ3XeHk5zsg076WG9Kb+50OlmwYAF33303w5v5P1DdSwOGUs0o27+fgYDTRzu/HQ8IYLTVyvurSpk/v+WA0UPZzXtVevNFixaRmprKz372s3Y9T3UdHcNQqhkVBw4AUNnEvNnycs82rC++2Pb67AMHEuh0smZZFhde6Ktedr/uSG/+29/+FqvVyl96KnKqJmnAUKoZocXFACQ3sQbD5YJbb4WxY9teX9r3vgeA7XDfXrzX1enN8/LyeOSRR9i/fz9nn302kyZN4uWXX+6Wn021TNObK9WMAzfcwOiVKzE2GwQGdr7CXbtg8mRWXn8/p85/jMWLG7Sn6c1VO2l6c6V6i9xcKiMjmwwWLpdnt7322FdeDkDgiZxeO7VWqZb4JGAYY+YaYw4ZYzKNMU0uyzTGXG+M2W+MyTDG6NIl1euFFBZSFhHRZNkzz0B4OLRnT6XgAQOwBQeT4ncc710cpfqUTs+SMsb4A88Ds4E84AtjzGoR2V/nnFTg18B5IlJqjGl9Tp1SPUhECLdacXm3V23o7LPhzjtb3wujrtSRI6lITCS0oAARaXZGUVPHlWqoJ4YTfHGFcQ6QKSJZImIH3gK+3+Cc24DnRaQUQEQKfNCuUl2muqqKSKsVGTy4yfILL4Q//7n9GUOcSUlElZbywx+W0zDNUkhICMXFxbrJkmqViFBcXExISEi3tuuLdRiDgdw6j/OA6Q3OGQlgjPkE8AeWiMgHDSsyxiwCFgG1MyyU6gnW3FwG2u1Ux8U1WV5V5clW256A4XQ6yXa7GWW1kptroaQkkrrVJyUlkZeXR2Gh7pmhWhcSEkJSW1MN+IgvAkZTL5mGH5ECgFTgIiAJ+NgYM15ELPWeJPIS8BJ4Zkn5oG9KdUhZRgYDgYpm8khNmwYTJsBbb7W9zoCAAExKCqE7dvDCk3mMHFn/Q1FgYCDDmtkKVqnewBe3pPKAugmek4ATTZzzvog4RORr4BCeAKJUr2SOHwcg5fzzmyxfvBhuuKH99Y6//HIAag4d6nDflOopvggYXwCpxphhxpgg4AfA6gbnrAIuBjDGxOG5RZXlg7aV6hJu70rk8NGjmyy/80646qr21+vvzYm0e3UWzz3X4e4p1SM6HTBExAksBj4EDgBvi0iGMeYPxpj53tM+BIqNMfuB/wL3ikhxZ9tWqqvYsrJwG4MkJjYqs9uhtLT96zAAvvCmx4guy8FiaeVkpXoZnyQfFJH1wPoGxx6s870Av/B+KdXrBeTnUxEeTkSdZHyn7dwJ554L69fDZZe1r97+Y8fi8vcnWfK49rc+6qxS3URXeivVhAirFUdiYpNrIoYMgaee8gx6t9dZqanYExIILSxscRMhpXojDRhKNSAihBUX4xgwwHPfqcG9p6Qk+PnP2755UkPOwYOJsliYPduKN2O4Un2CBgylGrDV1BBpsVAeHQ2jRsGcOfXKCwuhpKTj9ecAUVYr/ftb8KaXUqpP0IChVANlubkEORxURkdDWho0mCl1330dux11Wr+xY4koL+c3vyygia02lOq1NGAo1UBNZiYASTNnwjvvwLPP1itfuBCefLLj9SdfcAEGqO7j+2KoM49u0apUA/ajRwEISUyEzz+HgQPB7eb0nqoXXNC5+k1KCgCfvX2UvVVwxx2dq0+p7qJXGEo1UO7dmtXk5sL06TB0KDxYO0uc/fuhrKzj9We73QAMdmYTHt6prirVrfQKQ6kG/E6cwMA4iSEAACAASURBVG0M5vLLYdgwT3Tw7mrmdML48fC738FDD3Ws/rjJkwFIcp3gvB/5qtdKdT0NGEo10K+0lKqoKMJHjIARI+qVicCKFY3GwdslIi4OW2wsYUVF2O12goKCOtljpbqH3pJSqoHgggJscXGwbRtkZHjGL7Zvh337CAyE667r3CwpAMegQURZrVxyieYHUX2HBgylGggtKsIaFQW33w6//a3nsmLePHjmGQoKYPduTz6pzigMDSXaYmHixFK8QxpK9Xp6S0qpOmw1NURarZyKjYX//V/PQX9/T+Ko1FRW/csTR44d86QI6ajICROI3LmTm35Ugp9+bFN9hP6pKlVHeV4eQQ4H/dPSPKPb48d7Cs49F+LjufRSWLkSmkhi2y79J08mwOWiKkuz/Ku+QwOGUnWcXkwXEBsL775LbQ7ymhr4+99JPvU5114LgYGda8ftTUT1yZuZ/OMfnatLqe6iAUOpOkr37gXAUVgI114LeXmeAn9/+OlP2fnsJ/hiszzHoEEADPPLRrevV32FBgyl6srNBSDgyivhyy9h5EjP8cBAyMjgzkM/4+67O99MyKhRAAx05HPhhbp9veobdNBbqTpCiopwG0PCzJkQ0ODlkZzM0qUd22mvkchInBERRJSUUFlZTXh4mA8qVapr6RWGUnX45+dTFRUFq1bBli31C0tLOfvt+5lSvc0nbdkHDiTKYmHevFKf1KdUV9OAoVQd/vn5lEVFwa9/DS+8UK+sSkJZ93+ZnNpywCdtVfTvT5TVyqxZRT6pT6mupreklKojvLQU65AhsGmTZ2ZUHZl5IVxR8w4rR8O1PmgrYvx4/L78kksutvqgNqW6nk+uMIwxc40xh4wxmcaY+1s471pjjBhjpvqiXaV8yWG3E2m1EpyaCv37g3cm02kjRsCOHXDxxb5pL3TUKILtdiryTuBw+KZOpbpSpwOGMcYfeB64DBgLLDDGjG3ivAjgbuCzzrapVFco9+60Z6Kj4S9/gYKCeuVhYTBj0DFib70Ktm7tfIPe/TU+eTOTjRs7X51SXc0XVxjnAJkikiUiduAt4PtNnPdH4EmgpokypXrc6TUY1WVl8POff7Noz2vXLljzaTySsR+KfDDu4A0YI0OzGDas89Up1dV8ETAGA7l1Hud5j9UyxkwGhojI2pYqMsYsMsakG2PSCwsLfdA1pdrO4U3TYS6/3BMQzjqrXvmyZfDjRaGYw4fg6qs736A3YCRUFzBmjK7FUL2fLwKGaeJY7V+/McYPeBr4n9YqEpGXRGSqiEyNj4/3QdeUajv//HwAhp53HsTGelZ31/G73/nmTlSthATcwcFElpSQm1vuw4qV6hq+CBh5QN28nUnAiTqPI4DxwBZjTDYwA1itA9+q18nNxW0MgatXe/JINRAfD2lpwN69MHZs56OHMTgSE4myWrnrLp1aq3o/XwSML4BUY8wwY0wQ8ANg9elCEbGKSJyIpIhICvApMF9E0n3QtlI+48jOpiIiAlm61LNwr4E334SdO4GkJM/WrQ2uQDrU5uDBRFssXHDBsU7XpVRX6/Q6DBFxGmMWAx8C/sAyEckwxvwBSBeR1S3XoFTvEF5aSmVMDJH79zdagyECt9wCP/kJTPlzDKxb55M2Q0eNgt27Oeccn1SnVJfyycI9EVkPrG9w7MFmzr3IF20q5Wv9SkuxjRwJxkBoaKPyI0capDW32z1XGZ240vAfNozwykrys05hmwnBwR2uSqkup6lBlAKcDgcRFgvu6Gj4n/+B48frlRvjuRM1YID3wJYtEBUF6Z28s+qdKbVnXRa7d3euKqW6mgYMpfhm0V6Z0+nJIdVg6XVWludw7fKLsWPhzjshJqZzDZ9eixF8RPfFUL2eBgylgMqDBwHwO+88qKysfSM/bccOuOuuOgEjIQGeeuqb/TI6yttObEURcXHOztWlVBfTgKEU4MzOBiBx2jTP/SdTf3nRDTd47lKNGFHnoAhkZoLL1fGGBw9G/PyItlrZv1+TEKreTQOGUnyzyjt03Tp46aVG5QEBnlyE9fZUevttSE31rMvoqMBAnAMGEGWx8Oij+R2vR6luoOnNlQIqDx3CbQwBX34Jbnej8pUrPcMaP/xhnYMXXghLl8LgwY3Obw+TnExUfj7Tp2fiWeOqVO+kAUMpPFNqKyIjidi2rdHtKPhmHLxewEhMhNtv73Tb/sOGEX3wIBMmhHe6LqW6kt6SUgoIKy7GFh+PaSJYAPz737C2qdSZVits2NCpjb5NSgqRZWUcPVyMU8e9VS+mAUMpILS4GHtUFFx/PeTkNCr39/csu2jknXdg3jzwzrLqkORk/NxusrZlk52tEUP1Xhow1BnP5XQSbrFgDQjwbHoRFFSvvLoaHnjAm0eqoe99D/77Xzq1oUXtWoxMwsIqO16PUl1MA4Y641Wc3mlv4kQ4fBgGDqxXfuoUPPEE7NvXxJMTEuCiiyAkpOMd8AaMmLISYmM7UY9SXUwDhjrjVR0+DED0+KZnKKWkgM3WYMC7riNH4K9/7fg4hneJd5TFwpdfWlo5WameowFDnfGqvQEjeNMmePzxJs/x92+QeLCuf//bk8b2WAdTlPfrh6t/f6KtVv7618yO1aFUN9CAoc54Zfv3A+BfU9MorTnApk3wm994rjKatGAB5OU1SifSLkOHEmWxMG3a1x2vQ6kupgFDnfFCi4tx+/kR/N57sGRJo/LPP4f/+79GY+HfiInp9OI9v2HDiC4rIzU1tlP1KNWVNGCoM15wQQGVUVGE9OvXZPlvfgMVFU2u5/vGJ5/AL3/Z4XEMk5xMlNXKgf2lHXq+Ut1BA4Y64wWePEl1ZCTMnOnJY96EFoMFePJJvfIKFBR0rBPJyQTa7eTvO86BA7pdq+qdNGCoM15IURFlYWEQFtbk6rwHHoDXXmulkoULobCwzg5L7eQd/0gNyiI4uOOrxpXqShow1BnN7XIRWVaGjBgBmzdDbOMxhPXr27CxXkhIg1S27eQNGFGWEgYMiOt4PUp1IQ0Y6oxWkZtLoMNBcL2NLur76ivPoHer1q2DSy+lQwmhvAEj2mrliy90LYbqnXwSMIwxc40xh4wxmcaY+5so/4UxZr8xZo8xZrMxphPzD5XynfLTU2p37ID77utcZTYbFBfDyZPtf27//rhDQ4myWFi2rKkcJEr1vE4HDGOMP/A8cBkwFlhgjBnb4LSvgKkikga8AzzZ2XaV8gWLN9+HKy6OpjbVPnAAfvxjz7+tuvpqz72rpKT2d8SY2plSM2cWIZ3IfqtUV/HFFcY5QKaIZImIHXgL+H7dE0TkvyJS5X34KdCBV5RSvhfkndUU/uijntXaDRQUwMcfN7mer3kdnVqbkkJMeTlJSc2nWVeqJ/kiYAwGcus8zvMea84twIamCowxi4wx6caY9MLCQh90TamW+Z84gdvPj5gxY5osv/BCyM6GyZPbWOHbb3uuMKwd2J87OZkoi4WMDIteYaheyRcBo6mPQk3+tRtjfgRMBf7UVLmIvCQiU0Vkanx8vA+6plQrcnOpjIjApKRApg/yOA0ZAhdfDGVl7X9ucjKhlZUUHzvJli0fdb4vSvmYL7ZozQOG1HmcBJxoeJIx5rvAA8CFItJcVh6lupX/yZOURUQQ8d3velKVN/D00/D11/DMM22scOZMz1dHeGdKDfPLYeDA73WsDqW6kC+uML4AUo0xw4wxQcAPgNV1TzDGTAZeBOaLSAeXwirle5EWC47Bg+Hvf4fIyEblJ054bkm1W2kHUnzUrsWwMLDBnhxK9QadDhgi4gQWAx8CB4C3RSTDGPMHY8x872l/AsKBlcaYXcaY1c1Up1S3EbebcKsVM2RIs+f86U+wur1/ra+95lkA2N5053X2xfjss0Idx1C9ji9uSSEi64H1DY49WOf77/qiHaV8qeLYMSIcDti/H26/HV58sV653d5ChtqWzJzpyXrb3icPGoQEBBBltfLyG9uYPj2eqCY3EleqZ+hKb3XGOr0Go3rYMLjkkkblixbBvHkdmCU7YgQ8+CAkJrbvef7+kJRElNXKJZcY/P3929mwUl1LA4Y6Y4n3llH4//t/cMMNjcqnToVzz21DptqmOJ3w6afgdrfraSY5mdiKCiIjIwkPD+9Aw0p1HQ0Y6ozlyskBoP/YhokJPBYvht/+toOVr1jhuTW1a1f7npecTESJhYyMUioqKjrYuFJdQwOGOmPZjx7F5edHyPjxcORI7XGnEz78sN0XB/VdeqknaLSQ1LBJycmEl1kpLylk5cqVneiAUr6nAUOdsSQ3l4qICMyvf10v/9O6dTB3rieteYfFxcH11zc5VbdFycn4iZtB7nzOP//8TnRAKd/TgKHOWJFlZVTHxWEefhhCQ2uPX3YZrFzpCRqdUlQEr74KVVWtn3va6am1VivR0dGd7IBSvqUBQ52xQouLccTHN9q/IigIrr22c/shAZ7MtQsXevb7bqs6i/e++qqQrVu3Ul5e3smOKOUbGjDUGUncbsJLS3EVF8PNN9ce/9OfPBcFPnHhhZ5B71mz2v4c7xVGtNVKVlYxW7Zs4Uid8RWlepJPFu4p1ddU5uYS7nRSPnRobcAQ8azqPuusejGk40JDYeLE9j0nJAQZMIDosjKGD69mwYKfEdnecRCluoheYagzUtWhQwCEz50Lc+YAnvUWW7fCX//qw4a+/hp+/Ws4darNTzm9FsNisdQGC3enpmypb7uqqqpuSSWjAaOPcbvdtfe07XY71dXVPdyjvsmelQVAv/h4cLtxuz2pQIyBsDAfNlRcDH/+s2dj8LZKTqZfkZX9+4+wb98+0tPTef7553G5XD7smPo2eeutt3jrrbeoas8Eiw7QgNHHrFu3jmXLlmGz2Vi5ciWvvfaafvrsgMqDBwEIu+su+PprNm3yDB/s3evjhs4+G0pK2jflKjmZiNJSXI4E3n33XTIzM0lKSsJm010BVGMiwpQpU5gwYQKHDx/u0rZ0DKOPmTJlCgkJCQQHBzNjxgwqKyvx8+u7cd/tdmOM6fYtSe1ff43Lzw+/hx+GoUOJKfGMTY8c6eOG/PwgIqJ9zxk6lCC3jQfvnMfmfXv59NNPSUpKwtlgNpdSAMYYJrZ3rKyD+u47zRnEZrNx0PuJeNCgQUyfPh2As846i7S0NAByc3MpLi7usT521I4dO/j444+7vd0oq5XKqChCf/YzCAxk2jR44w0IDu6Cxg4cgGuuAe//Yau8U2sDjucxfvylWCzXUlBQwAsvvMDOnTu7oIOqryooKGD37t3ddrtSA0YfsHXrVt555x2szewT7Xa7WbVqFWvWrOkTeyhUVFRwyjsIbLVaOXjwYLf3O7CggOqYGExBAVu2eO4adZnQUM+ajLy8tp3vDRjk5LBjB7zyyjjOPfdW7HY7a9euZevWrX3i/1l1vd27d7N+/XocDke3tGd66x/e1KlTJT09vae70Ss4HA5OnDhB8uk3kiYUFhYSFBTU6/dPEBH+9re/4Xa7uf3226mpqSE4OLjbbqtV5ORQ9KtfkfTuu+QNHcrA8y5h4LqXmTsX3nyzCxsWaXPa283vV3DflQexJSbz6EvxzJgB8fGQmZnJzp07OXjwIKmpqVx11VWE1lmhrs48IkJxcTFxcXG1x4wxO0Vkale0d8ZcYYgIu3fvprKysqe70iYWi4X3338fp9NJYGBgi8ECID7es9mOiLB582ayO7SvaNeprq5GRDDGMG/ePK655hqMMYSGhuLn54fb7e7SwXtrbi5HfvQjAkeNIvnttzk0aRLHr76a4HvuYMsW+N3vuqxpjzYEi8OHYf58+O6V4RSbeBw1bubPhxtvhH374MsvR1BWdj3z5s3j6NGjvPjii5w4ccJnXTz9+6+oqOD555+vvQ1aU1PDzp07sVgsgOe11Fs/aJ5JTr+e6gaLrnbGBAyr1cqaNWv47LPPACguLsZut/dwr5p3/PhxDh06REk775XY7XYOHTrU5bMl2sNqtfL888/X/u6TkpKIj4+vLS8vL+e5555jV3tTgbdBcW4uGQsXEjh6NKlvvEHhpElYt20jafNmzv7jH2HqVNLSYMwYnzdd3969MG6cZ6FHA6Wl8ItfeIq3bIHHH4cDY69h75jr+b8HCkhPFyZOhHvvhWXLnHz9dQ7XXHMNAMuWLeOLL77o1Bu4iPDPf/6TjRs3Ap7AkZCQUHv1UlxczNq1a2tvI+bn5/P444/T2+8A2O128vPzv5Wzy8rLy3n22WfJ8k4P7zanPy30tq8pU6aIL1RVVdV+f/z4cXG5XOJyueS5556TV1991Sdt+JLT6az9vrq6ukN1VFdXi9vtblRfdzvdB7fbLRs3bpSTJ082e96qVavk6NGjPms7PydHdt56q5SFh4uAnJw8WV6+4w7Zs2eP54Tjx2X70t1yz08cUlTks2abZ7GIzJkj8vHHtYccDpHnnhOJjRUxRuS220Rqf0ULF4p4bmRJUXiyLB70rvgbp0SF1sj1F6+WXV9lSGVlpbz22muyZMkSef3118VqtbarS+Xl5bXff/jhh7J9+/Ymz3O5XGKxWMRms4mISHFxsfz73/8Wh8Ph/Tkc7Wq3K7lcrtp+5ubmypIlS+Tw4cMiInLs2DF59NFH5euvvxYRkRMnTsibb74pBQUFPdXdDisoKJBXX31ViouLG5UB6dJF78s9Hhia+/JFwDh27Jg89thjkpmZ2agsOztbcnJyRMTzhuVyuTrUhs1mq33BuFwuqampqX2TdrvdYrPZauu22WySl5dXGwiKiorkgw8+kJKSEhEROXr0qDzxxBNy4sSJDvWlocrKSnn++edl586dPqmvPbKysuSFF16QioqKbmvTbrdL5qFDsv3226UkOloEJH/ECKn64AMR8fw+ai1dKs/yE0mIc0o3drHWhg0iY8Z4XoEXXyyya1eDE5xOz8FXXhG56y6Rc86RjMCJcikbBERSOSzvjrxPXPf8TA49+aQ8umSJPP7447Jnz57aQN2S3bt3yx//+Mcm33Daw+FwyIsvvihbtmzpVD12u10sFkun6nA6nfLiiy/K2rVrRcTzYfHAgQO1f4NFRUWyYcOG2tdbTk6OLF26VE6dOiUingDakx+wfKXXBwxgLnAIyATub6I8GFjhLf8MSGmtTl8EDJvNJqtWrWr1k9f27dvl5ZdfbvMn+tMB4OTJk7JkyRL58ssvRcTzB7lkyZLaT7GnTp2SJUuWSEZGhoiI5OXlyZIlS+TQoUMi4glojzzyiGRlZYmI51PD6tWrO/3COc1ut8vbb78t2dnZLZ5ns9nk1KlTYrfb21y32+0Wi8VS+yZssVjk5Zdflu3bt8uqFSvk9ZtvlhcWL5Z1a9e2+U3JZrO1+Q1PxPMCz8jIkP+8+aZsvuUW+fLss8USGSkCUjZypHz20EPy1j//2fSHgfx8kfXrpaq8+94g3G6Rz7fZ5LJLXQIiI0aIrFrlOd4mdrvIrl3y+k0fSqLJFxCZ7fdv+YIp4ozpL/suuUReuu02WfHWW00GaqvVKqWlpSLi+d1t3LixfhBtz8/icIi7ulocDoesW7eu9m+6I44fPy6PPPKILFu2TPbt2ydlZWW1b+JtUff1snXr1trXW3u43W5ZtmyZLF++vM1/fz3l2LFjLb5XdWXA6PQsKWOMP3AYmA3kAV8AC0Rkf51z7gLSROQOY8wPgKtEpPEmynW0NkuqqqqKwsJCioqKKM3Lw5GXR0B8PH4xMRzP92w+Ex0dTWRkJP7+/i3+DBkZGRw5coTvf//7TS4gczqdlJeX1w5EJ4SFMUQEyc6m+tgxAseOJfbii4kePJi8vDxGjhxJfHw8lZWV7Nq1i1GjRhEXF4fNZiMnJ4fBgwfTr18/3G43FouFY8eOcTIjA/tnnxGYmYkrJga/s84iZNQookaOJCExkfj4eEJCQlr8OVoiIuzZs4fQ0FDKrFYqMjKQvXsJPHiQyNxcYouLKUpMpHLWLGJuuIGzJkwgKCioXh0VFRU4nU6io6Oprq7mySefZPbs2UybNo1dX33FgWXLGJuezrj9+wn1piyxREWRmZpK6cyZxF53HWOmTm12Zk96ejrr1q3jtttuY9CgQY36X1hYSG5uLvkHDyJbtpCwdy/DsrJIKCwEwBEeTs7w4UTfdRdxixbhFqk3+0q8E5XKyz0pnrxLWBqpqoKsLDh61PNvdDR85zuezfPau77Q5YJt2+Bf/4JVK+3knAgiKszOg38MYvFiTyr1jnjvvU9YvryabdtmUVpqSOlXwPdr3uIq17uMjDvMkWmTGXzffZx14YXefrh4+umnSU5O5rrrrmtUn4jnd7JxYxU5OVXEx1eQkFBGTEwJIhXIiRP027eP6EOHiM3KIjE3F3+Xi5OTJuGYP5/YhQuJHDKEffv2ERISwohWdhq0WCykp6dTtnMnkZs2MfrAAWJLSshPTOREUhJ5gwYRM3cuE+bMITExsdlZdLt27WLNmjXccccdteNiIkJVVRVlZWU4nU5cLhdOp/Ob7x0O3FYrfkVF+BUVYex2QiZOpCQ4mNCwMCZMmICIYLfbCW5lYY7T6eTUqVOcOH6ckj17cH3xBQE5OcjgwfiPHk2/iROJGzaMxMREwsPD271A1W63U1lZicvlQkRwOBy89tprJCYmcskll9ROFnG7XIjDQXhMDAMHDuyyWVK+CBgzgSUicqn38a8BROSxOud86D1nhzEmADgJxEsLjcfFjJdx/Vcw8Tv9MdFBFOw/TsnREvoF2pFKG1Q5MTYBG7hdfjgJIBgbIVQT6O/EP9CFf6ALvyA3AaHgH2RwuA3RQyJwBwZjLXVTYvWnX0I4Dr9gyq2CtdyP4JgwKm2C0x6Aqa4mzFQTbK8m2FFDqL2GUHs1wS4b/rjqfRkEe2gQ1RFhuOJiqAqIoKAmlpk/GEdkZBgZa7LJznIy6dpEio4VUrznGJRX4V/lwlS7cTkCqCYUO0GEUEM4FfSjkjBTSUCwC/9QFyYE+iWG0m9QFMWn3FQ5ghk0dSBV9gCO7iqn2hlAZEp/qmv8OZljx+4KJCAyDFNZQUCFhTB7JeG2CvrZKghx2wjASSAOCAvCTiBRtmJCXNX4+bkoi4miLHYA8bOnk5iazKZnD+Mmn6RZYUyefC7vPnYIvwg7YfFu7Afz8Su2U+MIxWqiKY8YxKnKCCKCbQwMKia27Gvi3aeIMaX4xbgJHRVLdkksE6anMO+XE6gqd/Lkgq8YNS2C0IlxlJcG8N/38gmLCcAhAViKoaIsgACng3BHORGOCkKoIdjYCezfj+C4KLKPB5I4IZqKkCKS+ifw5SYHs+YEMH72QAqOWPjD//bjd7ed4pKbkjiyvYCF9ybw0xuLmXl5LEePUu+ruYlHCQmewHH6a9IkCAxsfF5NDWza5AkSq1d79lEKDobZs9xc5VjBlb+dQP9BIfDcc56Nw0eM8Czqe/FFuPtuGDYMMjLg5Zfh5z/35CzZvRuWL/eMfA8aROGHH+JYtoz91z7HovviSU2w8lF6GDZXINHGwpXyL77PKkYkHyN9xN3ccM8Qqt99na2X/omP0+N5dM4WMv6xk6XBd7J1u4PCkwGUWPs1/XObU6TKEc7iKMPM18RFWRmaUk2k6yQhWaegUijziyRnwHD2xE9mT+BU5o81VOeW8ln4LEpKYELoEaqKK9nvl4Kr9BQh1hICyp2Iw2AjmOrACExAAP1cZSTY84ilmBhKCQuuxEQbAqJdDBkezf5hNxEY2Y8LS9+lwB7B+4UTqSgqoL/dSmllCAUVEVRWB1NpC6PGGUI/U0mEKSOKMqLcVmJcpcS4SogWK5GUEUkZYVThwp+qwDBKIvpTFhlPTXAolqAwcvvNJHZQHCNOfYU7NpH/Hh9KTHQFYScP4qhxIeXVBFXZCax21L4PufAnlGrCqCKMKgICnZhQwYQbnGH9CBoQTrzzJP0mTmJrZiKDYioIyDmMNSCWk6VB+Ilgr/S8vm2OIBzOQIKwE0o1oVJNmFQRJlX0c1cS5q6in6uCCFcFoVRjCwziPsefe3XAuBaYKyK3eh//GJguIovrnLPPe06e9/FR7zlFDepaBCzyPJoyBb65wvDHSRB2grATjI1A4wSEsFAhNCIAW5WL4vJgwvsJdncA1TY/atxB2AjGTctXGIHeOut+BWHHjZ83HATgMIHYCMIvOBC3XyA1NkONK5DgUD9cTsHpENz4IW2YeBaEzfOfTzXB/g5CQw32GqHSFcLYc6OoqXSRk1FBjSOAwLAAKmyBVLpan29vcHuCjPcPNYwqQqhBMDgJwEYwruBQnIFhWCoDcZpAgqNCcLr8qCjzTKl0t2PinB+u2hddgHHiHxzI8O8MIjImgLXvVJMSV8HgifEUF7k5uKsGMX7USNuuksKopD8ltV8xlAJgi0qgKjyB3cfjiI4LJHxgBBUWJ9m5fkSEuiAwgOoqcDjb90lu0CBPWvOGX8OHQ2Gh5yph2zb4+GPPJ3HwJCmcMeObAFJU5AkSGzZARYVnd9bLL4errvKkkqqXIWTbNrjiClizBs4/H/7zH8+JGzfC9OnwwQdwww2eaVOTJ8P778NNN8H27Z7pVO+8A7fcwrL7D/OPDwfw3oKVBN57D3de9Alrt6aAzYa1JoRQqriUD5kfsI6h7mOsj7ye7WXj2SfjqZBwAIb4HWOm/w7O4xNm+KczrOYwuQwlkxHsDZ7Gp84pOCLjOGrpT54kNfzVtcjPz42fnxAeUEGwu4pwZyWh7mpCqCHQ30lpaBKl/nGcG7kPqajkY7+LsNe4iXYWUmIPp1KaDmINhVJFLMXEUUQsxcRSTD8qKTORlPrFUkUo5X7RFLhiqTb9qHK1/r5QVwAO/I2LQBwEioNAHN984PI3hLirCQgLJqcqgYg4Q0zh11SFxVNUFYbbz48ad9vaC6Oy9sNiOBWEU0Hw/2/v/mOjvu87jj/ftrEd+wD7TCAHgQaCMwtSlaUIQcM6QoLKa3U2iAAADTlJREFUqnV0/bGELdRVg5Zsa6SpmqpE6ySWZtqKNqnJLIUgLVUBLQKaVTVitVNYQxsaCgaMi1NDDAm/bGoMzD/OPv+4e++P7+fMcdj47Lvz3S3vh3S6r7/fz/frl8/f773v++M+XxlgUAvpzy+lJ1LKYP49DIQLCEkJfZFihojfTU3f9zBSUTC+CnwurmCsUNXnY9o0uzaxBWOFqo7Zl8WywDytW72WSOd1CFQQnDOLxevWIYsXc6yjgxPNzTz77LMA1NXV0djYyPPPP09p6e0rmKp3Q7VQvxLq6Gaw438plCGK8oa8AiRD5IWHvEZD3vNAb6+3Kzp3LrpgAeL3J3wsQhXCHdcZajhJz69+TejYCfLeP0Nx+++IBMop+FQVJatXUbxqlfdGkOAX7SIR6O+HYEeQ7pbLXH//I4oYwFc4jK9wmOmFwxTnh8kXvOMg0Yeq94n14YchEAARIpEIe/bsoaqqimXLlt2RPxKBoUHlw9r/pvXV16i62Ebp5asMU8D16X5m93RQRheRB+dS+Mw3KKqu9t5xExAKQcex85zb+VM6f97E8Pku7on04ysJ4ZsZZua9+RT4hqC8hIUrV1IQCMCcOd4315YsGbVfps7OThobG3nsscdGDj8ODkIwePdHWZlXFBYunFgPtVeueDfRixaRU6e81wy8qBs2eO/9a9dO/pBTIlR15OZKDz30EC0tLezevZtNmzYxf/4iDh1Sdu3sZd9/hbnR693uVYhQVdjCp+85wcryZlbfd47K+4a4Z8YMpLDQC1xRAStWeI/Zs2/7naGQVzBbW72OeH0+KC11jxKl8OwJqH0T/8Fa7u38iPy8MKHiYkr7+gjn5XFz+XIqNm9GNmy4Y9mjGRz0LjvuPHONi28fo/29FnouXKPCD4HANGbPK2beJ3yUzZ9Fnt/vbU9lZbceYxz+VPUOPXZ3e+tCQYG3tzhtGhTkK9PaLyLNp4icaiBysoHhkyeY0d3NYGUl/VVV/CoU4pFnnmH++vVc7+mhqamJ5cuXM32U9VNVGR4OEw7n09/RS+eJFq4ePU2ko4vSGdNo7/yIqmWVPPjwIvJmTvde1OnTvYfPN7Ib29zcTFtbG2vWrGFazK5tJOL9X0Ihr+uy8vL0FYxUnPBeBdTH/Pwi8GJcm3pglRsuADpxxWqsx7x58zQUCqmq6uHDh3XLli0jJ3qampp07969I1c0dHR06OXLl8c8CZQVsuhEWigU0h07dmhDQ8Md0yKRyMgJ0+HhYT106JB3me6HH+qNl17Sa5/5jF7bvFnDp06lJMtQX5+G4y7LrKur0z179iR88vHdd9/VrVu3and3d0oyTVRXl+rbb6sePuxd3DRVIpGI1tTUaH19vap6/69Lly7d8boNDAxqTc0R/fa3D+j+/e/qxYsX0381UCSinfX1emHjRu3ZsEHDu3bpjXEuvshW/f39+tZbb41c9h0KhfT48eMpWd+6urq0pqZGT58+PW7bgwcP6rZt28bdLsjyk94FeCe9Hweu4J30/nNVbY5p8zfAJ/XWSe8vqeqf3W25lZWV2tjYSGlpKQMDA6gqRUVFU96r6f9XqjryWl69ehW/309hYSH79u3jwoULPPfccxQkfVPr5LMlKhgM3rF3+XHQ2dlJX18fC9ytXU3uCYfDI3vG58+fp7y8nPLy8nHbjiWruwZR1WHgm3h7Eb8F9qhqs4i8JCJ/4pr9B1AhIq3At4AXxlvuzJkzR94AioqKKC4utmKRQtHXcnBwkF27dlFbWwvAkiVLWLFiRUa7TI9mCwaDhEKhMdudO3du5JvwH8diATBr1iwrFjkuWgDC4TC1tbXs37//tumqSm9v721tM8U6HzS0trbi9/vx+/2ZjjKit7eXV155hdWrV/OH7tLQWOFwmJqaGvx+P5s2bcpAQmNS7+bNm4gIZWVlDA0Noaq0tbWxc+dOnn76aRYuXDjuMtK5h2E3UDLjXjOfCT6fjyeeeGLMbPn5+VRXV09xKmPSK/ZQVF1dHZcuXeLJJ59k5cqV3H//xK5QSwcrGCZrRW8UFa+9vZ1AIEBZWdkUJzJm6ixduhS/309FRQXr1q3LdBzgY9RbrclNN2/e5MCBAyNdb1+4cIHt27fT1NSU4WTGpNeiRYt49NFHMx3jNlYwTFa7evUqR44cob29HYD58+ezfv16lixZkuFkxnz82Elvk9VUlWAwSElJCcPDw3f0b2WMuV1WX1ZrTDqJCD6fj6NHj/L666/f9TJbY0x6WcEwOSEYDLJ48eKkeuw1xiTHrpIyOWHt2rX2xU1jMsz2MExOsGJhTOZZwTDGGJMQKxjGGGMSYgXDGGNMQqxgGGOMSYgVDGOMMQmxgmGMMSYhVjCMMcYkxAqGMcaYhFjBMMYYkxArGMYYYxJiBcMYY0xCkioYIuIXkZ+JyAfuuXyUNstE5D0RaRaRJhF5MpnfaYwxJjOS3cN4ATioqpXAQfdzvD7ga6q6FFgPfF9E7GbMxhiTY5ItGBuAH7rhHwJfjG+gqmdV9QM33AZ0APcm+XuNMcZMsWQLxhxVbQdwz7Pv1lhEVgCFwLkxpv+liDSISMO1a9eSjGaMMSaVxr2BkogcAO4bZdLfT+QXiUgA2AlUq2pktDaquh3YDt49vSeyfGOMMek1bsFQ1SfGmiYivxORgKq2u4LQMUa7GcB+4DuqemTSaY0xxmRMsoekaoFqN1wN/CS+gYgUAj8Gdqjq3iR/nzHGmAwR1ckf+RGRCmAPsAC4CHxVVW+IyHLgOVXdLCJPAz8AmmNm/bqqNo6z7B7gzKTDZd4soDPTIZJg+TPL8mdOLmcH+D1VnZ6OBSdVMNJJRBpUdXmmc0yW5c8sy59ZuZw/l7NDevPbN72NMcYkxAqGMcaYhGRzwdie6QBJsvyZZfkzK5fz53J2SGP+rD2HYYwxJrtk8x6GMcaYLGIFwxhjTELSXjBE5A0R6RCR0zHjtojIFRFpdI/Pu/HrROS4iPzGPa+NmefTbnyriLwqIuLGj9vFepbk/ycRuSQivXHLLxKR3e7v+rWIPJBN2UWkRET2i0iL66L+X6Yie6ryu2l1InLK5d8mIvlufE6sOzHz1sYtKyfyi8g7InImZp7ZbnyurD+FIrJdRM667eDL6c6fom13ekzbRhHpFJHvJ5VdVdP6AD4LPAKcjhm3Bfi7Udr+PjDXDT8MXImZdhRYBQjwU+CP3PitwAtu+AXge1mafyUQAHrj5vlrYJsbfgrYnU3ZgRLgMTdcCPwy5rVPW/YUv/Yz3LMAbwFP5dK648Z9CfjPuGXlRH7gHWD5KPPkyvrzj8DLbjgPmJXu/Klcd2LaHQc+m0z2tO9hqOovgBsJtj2pXhfo4H0zvNhVwgDeRv+een/hDm51pT5uF+vJSEV+N+2Iup5948Tm/xHwuIi395SsVGRX1T5V/blrMwicAO5Pd/ZU5XfTut34AryiF73SIyfWHRHxAd8CXo6bLSfy30VOrD/AN4B/du0iqhr9FnhWb7uxbUSkEq838V+6UZPKnslzGN8U7w58b4yxK/1l4KSqDgDzgMsx0y67cTDBLtZTaCL572YecAlAVYeBLqAitVHvMKns4t346gt4N8uCzGSHSeQXkXq8zjF78DYQyJ1157vAv+HdjCxWruQH+IE7LPIPMW9MWb/+yK2bvX1XRE6IyF4RmePG5cy2C2zE24uIfliaVPZMFYzXgAeBZUA73sYwQkSWAt8Dno2OGmUZmbweeKL572aq/7ZJZReRAuBN4FVVPR8dPcry0/1/mVR+Vf0c3iHBIuCO8wNTaEL5RWQZsFhVfzzFOccymdf/L1T1k8AfuMemaPNRlp9t608B3h71YVV9BHgP+Ndo81GWn3XbrvMU3vY70nyUNuNnT9Uxt3GOxz1AzLG4u03D++ecBR6NGRcAWmJ+3gi87obPAIGYdmeyLX9c+/hzGPXAKjdcgNfpmWRbduANvGIxZdlT/dq7NtVATa6sO8BfAW3AR3h71oPAO7mSf5R5vh7z+mf9+oP3xhoE8tzP84Hmqcifwm33U8DZuHGTyp6RPQx3TiLqT4HTbnwZ3n0zXlTVw9EG6u1u94jISrc7+zVudaU+bhfrqTbR/OOIzf8V4H/U/RfTYTLZReRlYCbwt3GLm9LsLsuE8ouILzqP20v6PNAySv6sXHdU9TVVnauqDwCr8Tb8NW5y1ucXkQIRmeWGpwF/HJ2HHFh/XJ59wBo36nHgfTec9duus5Hb9y5gstlTWc3HqIRv4u0+DeF9QnoG7857vwGaXPDop6Tv4FXzxpjHbDdtuXuBzgE13PqWegXeMfUP3LM/S/NvdfNH3PMWN74Y2Au04l0JtiibsuN9clHgtzHjN6c7ewrzzwGOufbNwL8DBbm07sQs7wFu/1SZ9fmBUryrc6Kv/ytAfq6sP27aJ4BfuHkOAgtyYduNWdZ5oCpu+ZPKbl2DGGOMSYh909sYY0xCrGAYY4xJiBUMY4wxCbGCYYwxJiFWMIwxxiTECoYxxpiEWMEwxhiTkP8DGXk7r09XhdQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nm_laser = 543.5 # wavelength of the calibration laser, in fact it can be any real positive number (e.g. 1 is ok)\n", "fit = orb.fit.fit_lines_in_spectrum(spectrum, [halpha_cm1, halpha_cm1], step, order, nm_laser, theta, zpd_index, \n", " wavenumber=True, apodization=1, fmodel='sincgauss',\n", " pos_def=['1', '2'],\n", " pos_cov=[velocity1, velocity2], \n", " sigma_guess=[broadening1, broadening2])\n", "print('velocity (in km/s): ', fit['velocity_gvar'])\n", "print('broadening (in km/s): ', fit['broadening_gvar'])\n", "print('flux (in the unit of the spectrum amplitude / unit of the axis fwhm): ', fit['flux_gvar'])\n", "\n", "# independant plot of the two lines models and the real lines\n", "pl.plot(spectrum_axis, spectrum, label='line1 + line2 + noise', ls=':', c='0.5')\n", "pl.plot(spectrum_axis, spectrum1, label='line 1', ls=':', c='red')\n", "pl.plot(spectrum_axis, spectrum2, label='line 2', ls=':', c='blue')\n", "models = fit['fitted_models']['Cm1LinesModel']\n", "pl.plot(spectrum_axis, fit['fitted_vector'], label='fit', ls='-', c='0.5')\n", "pl.plot(spectrum_axis, models[0], label='model 1', ls='-', c='red')\n", "pl.plot(spectrum_axis, models[1], label='model 2', ls='-', c='blue')\n", "pl.xlim((15200, 15270))\n", "pl.legend()\n", "# In fact this \"very bad fit\" may be not so bad, and will be, in general, not so bad... \n", "# but its outputs are not constrained as in the bayesian fit and can sometimes be very far from anything realistic\n", "# if you obtain a good fit, redo the model and the fit multiple times\n", "pl.title('A very bad fit')\n", "pl.savefig('gvar_bad_fit.svg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bayesian fit\n", "\n", "Now let's say you have some informations on the broadening and the velocity of one or both of the unresolved lines e.g. there is some diffused ionized gas in the foreground which is everywhere in the field of view and you are interested into the point-like source emitting in H-alpha at a slightly different velocity. \n", "\n", "LSQFIT, a fitting module which integrates gaussian random variable as priors (initial guess) has been developed by G. Peter Lepage (Cornell University) (see https://github.com/gplepage/lsqfit and http://pythonhosted.org/lsqfit/index.html). This module gives the perfect answer to this problem. We can now inject some more information and help the fitting algorithm to find a unique and better constrained best fit.\n", "\n", "This algorithm has been implemented into ORCS. To use it you have to :\n", "\n", "- guess the SNR of the lines (yes, this is not so easy, but you can try with one rough SNR, do the fitting, compute the real SNR from the residual and then fit again, the only thing that will change is the uncertainty on the parameters)\n", "- define the initial guesses as random variables (we will use the package gvar which is intimatly linked to lsqfit - same author)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== velocity ===\n", "input velocity (km/s): 53(10) 7(10)\n", "fitted velocity (km/s): [54.8(5.1) 8.5(4.8)]\n", "real velocity (km/s) 50 10\n", "=== broadening ===\n", "input broadening (km/s): 18(10) 27(10)\n", "fitted broadening (km/s): [18.000(25) 27.000(37)]\n", "real broadening (km/s) 15 30\n", "=== flux ===\n", "flux (in the unit of the spectrum amplitude / unit of the axis fwhm): [0.81(18) 1.50(18)]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEICAYAAABMGMOEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXxU1fn48c+Z7PtGgECAsO8hICC4ohREVNywSq1Lq6JVv1Zbbf1+bRWt1qqt+nOpK4p1BVwRES20gCCgQQgCAQl7SMi+7zPz/P64mZCdLJNJAs/79corc++5954zkMyTe885zzEiglJKKXUits5ugFJKqe5BA4ZSSqkW0YChlFKqRTRgKKWUahENGEoppVpEA4ZSSqkW0YChFGCMOWiM+VkLj73RGLO+o9vUHsaYAGPM58aYAmPMUmPMtcaYrzu7Xap704ChOoUxZo0xJs8Y49fZbfGkxoKNMWaRMeZRN1c1F+gFRInIVSLyrojMrFWnGGOGuLlOdZLTgKE8zhgTB5wNCDCnUxtzEjDGeDeyewDwk4jYPd0edfLSgKE6w/XAJmARcENzB1bfiTxqjPnWGFNc/ZglyhjzrjGm0BjzfXUAwhgTV/2Xs3e982+utX2LMSbZGFNkjNlljJlQq7oEY8z26sc4i40x/s03zTxffexuY8z0WgVhxpiFxph0Y8zR6vZ7GWNGAi8DU6vfS74xZj5wLfAH1/urvkYfY8xHxpgsY8wBY8xdta6/wBjzoTHmHWNMIXBjvYY9DDwIXF19zZtq39kYY9ZVH5pUXX51c/8HSrlowFCd4Xrg3eqvC4wxvU5w/DXAdUBfYDCwEXgTiASSgYdaUqkx5ipgQXX9oVh3Nzm1Dvk5MAsYCMRT74O4ntOB/UCP6vo/NsZEVpe9BdiBIcB4YCZws4gkA7cBG0UkWETCReRVrH+HJ6v3XWKMsQGfA0nV73k6cLcx5oJa9V8KfAiEV59fQ0QeAv4KLK6+5sJ65edUvxxXXb64mfepVA0NGMqjjDFnYT0uWSIiW4B9wC9OcNqbIrJPRAqAL4F9IrKq+nHLUqwP5Za4GeuD+XuxpIjIoVrlz4lImojkYn1gJzRzrUzgWRGpqv7A3QNcVB38LgTuFpESEckEnsEKei01CYgWkUdEpFJE9gOv1bvGRhH5VEScIlLWimsr1WaNPftUqiPdAHwtItnV2+9V73ummXMyar0ua2Q7uIV198MKUE05Vut1KdCnmWOPSt3MnYeqjx8A+ADpxhhXmQ040sI2Un2NPsaY/Fr7vIBvam235npKuYUGDOUxxpgArMc+XsYY14ezHxBujBknIkntrKKk+nsgUFj9unet8iNYj7Tcoa8xxtQKGv2BZdV1VAA9muhwbiw9dP19R4ADIjK0mfo1zbTyOH0kpTzpMsABjMJ63JMAjMT6y/n69l5cRLKAo8AvqzuZf03dAPE6cK8x5jRjGWKMGdDG6noCdxljfKr7RkYCK0QkHfga+IcxJtQYYzPGDDbGnFt9XgYQa4zxrXWtDGBQre3vgEJjzB+r51N4GWPGGGMmtbGtjalfp1InpAFDedINWP0Rh0XkmOsLeAG4tonhoa11C3AfVmf2aOBbV4GILAUew3oMVgR8itVx3habgaFAdvU154qIqwP9esAX2AXkYXVOx1SX/QfYCRwzxrgeyy0ERlWPmvpURBzAJVgB9UB1Ha8DYW1sa2MWAG9V1/lzN15XncSMLqCklFKqJfQOQymlVItowFBKKdUiGjCUUkq1iAYMpZRSLdJl52H06NFD4uLiOrsZSinVrWzZsiVbRKI74tpdNmDExcWRmJjY2c1QSqluxRhz6MRHtY0+klJKKdUiGjCUUkq1iAYMpZRSLdJl+zCUOtVUVVWRmppKeXl5ZzdFdQP+/v7Exsbi4+PjsTo1YCjVRaSmphISEkJcXBy1UqMr1YCIkJOTQ2pqKgMHDvRYvfpISqkuory8nKioKA0W6oSMMURFRXn8blQDhlJdiAYL1VKd8bPiloBhjHnDGJNpjNnRRPm1xpjt1V/fGmPGuaNepZRSnuOuO4xFwKxmyg8A54pIPPAX4FU31auUcqPgYGu127S0NObOndvm67zwwgsMGTIEYwzZ2dknPqEFpk2bVjOZd/bs2eTn55/gjMY9/fTTjBo1ivj4eKZPn86hQx02z61J7f337SxuCRgisg7Ibab8WxHJq97cBMS6o16lugIRYeHChbz00kud3RS36dOnDx9++GGbzz/zzDNZtWoVAwa0bEHDBQsWsGjRohZff8WKFYSHh7epbePHjycxMZHt27czd+5c/vCHPzR7/Jo1a7jxxhvbVFdT2vvv21k6ow/jJuDLxgqMMfONMYnGmMSsrCwPN0uptiktLSU1NZXMzMzOborbHDx4kDFjxgCwaNEirrjiCmbNmsXQoUPrfMB+/fXXTJ06lQkTJnDVVVdRXFwMWB/KHZkLLi4ujuzsbA4ePMjIkSO55ZZbGD16NDNnzqSsrAyAffv2MWvWLE477TTOPvtsdu/eDcB5551HYGAgAFOmTCE1NbXd7bnxxhu56667OOOMMxg0aFBNMBAR7rvvPsaMGcPYsWNZvHgxUPffd+fOnUyePJmEhATi4+PZu3cvAO+8807N/ltvvRWHw9HudraXRwOGMeY8rIDxx8bKReRVEZkoIhOjozskd5ZSbufv789NN93EPffc49brLlq0iG3btgHgcDhYtGgR27dvB6w5G4sWLWLHDqvbsLy8nEWLFpGcnAxYQWzRokXs2bMHoOaDvK22bdvG4sWL+fHHH1m8eDFHjhwhOzubRx99lFWrVvHDDz8wceJEnn766XbV0xZ79+7ljjvuYOfOnYSHh/PRRx8BMH/+fJ5//nm2bNnC3//+d26//fYG5y5cuJALL7zQLe1IT09n/fr1LF++nPvvvx+Ajz/+mG3btpGUlMSqVau47777SE9Pr3Peyy+/zG9/+1u2bdtGYmIisbGxJCcns3jxYjZs2MC2bdvw8vLi3XffdUs728Nj8zCMMfFY6xJfWGvtY6W6PS8vr5N+OOz06dMJC7OWFB81ahSHDh0iPz+fXbt2ceaZZwJQWVnJ1KlTW3zNH3/8keuuuw6AY8eO4evry7PPPgvA6tWriYqKatF1Bg4cSEJCAgCnnXYaBw8epLi4mG+//Zarrrqq5riKioo6573zzjskJiaydu3aRq97+umnU1FRQXFxMbm5uTV1PPHEE1xwwQUNjr/sssuw2WyMGjWKjIwMANavX8+8efPw8vKiV69enHvuuXz//ffEx8fXnDd16lQee+wxUlNTueKKKxg6dCirV69my5YtTJo0CYCysjJ69uzZon+PjuSRgGGM6Q98DFwnIj95ok6lPCUtLY0PP/yQkpISrrrqKoYMGeKW69Z+bu7l5VVn28fHp862v79/ne3AwMA6267O7Lby8/Or0xa73Y6IMGPGDN5///02XXPs2LE1d1ALFiwgLi6uTX0F9dtWVlaG0+kkPDy85vr1rVq1iscee4y1a9fWOb+2zZs3A1YfxqJFi07Yx1L7OiJS53tzfvGLX3D66afzxRdfcMEFF/D6668jItxwww08/vjjJzzfk9w1rPZ9YCMw3BiTaoy5yRhzmzHmtupDHgSigH8aY7YZYzRvuTppJCYmkpdnjenw9fXt5NZ4zpQpU9iwYQMpKSmA9Rjsp5+6xt+DoaGhDBw4kKVLlwLWB3dSUhIAW7du5dZbb2XZsmUd/lf7Oeecw+LFi3E4HGRlZbFu3TomT55c55j9+/czaNAg7rrrLubMmcP27duZPn06H374YU2/WG5ubqeM5qrPXaOk5olIjIj4iEisiCwUkZdF5OXq8ptFJEJEEqq/JrqjXqW6gjPOOAOw+hn69evXya3xnOjoaBYtWsS8efOIj49nypQpNR3Lzz33HLGxsaSmphIfH8/NN9/s8fa9++67LFy4kHHjxjF69Gg+++wzAO677z6Ki4u56qqrSEhIYM6cOR3Whssvv5z4+HjGjRvH+eefz5NPPknv3r3rHLN48WLGjBlDQkICu3fv5vrrr2fUqFE8+uijzJw5k/j4eGbMmNGg76MzmJbcMnWGiRMnii6gpLqDI0eO8MYbbwBwzz33EBoa2qbrJCcnM3LkSHc2TZ3kGvuZMcZs6ag/yjU1iFLt4HQ666wM2R3H1ivVUhowlGqHoqKimqGuACEhIZ3YGqU6lgYMpdohNDSU4cOH14xC6gpDH5XqKLoehlLtYIyhsLCQXr16AZCXl4eInNRzMtSpS+8wlGqHlJQUsrKyiIiIwN/fn+3bt9dM2lLqZKMBQ6l22LZtG3a7naioKCIjI/Hx8WlyIphS3Z0GDKXawTUJKzIykl69elFVVdXmYbVdgbvSm9e2dOlSRo8ejc1mQ4fKd28aMJRqB9cM78jISCIiIhCRmn3dmTvTb48ZM4aPP/6Yc845xy3XU51HA4ZSbWS320lMTMQYQ0RERM36DF9+2Wj2/m6lvenNaxs5ciTDhw/3WNtVx9GAoVQbFRUVkZqaSkBAAF5eXjUBIyIiwj0VTJsGroR3VVXW9jvvWNulpdZ29foKFBRY2x9/bG1nZ1vbn39ubR871q6mdOX05spzdFitUm0UERFBTEwMAQEBAISFhWGMqdk+mXREenPV/WjAUKqNRITc3NyatQ1sNhthYWHk5OS4Zy7GmjXHX/v41N0ODKy7HRZWd7tHj7rb9RLetVZHpDdX3Y8+klKqjbZs2UJFRQWRkZE1+2w2G8nJyZSXl3diyzyjK6c3Vx1DA4ZSbbRv3z6AOgEjOjoaX19fbLaT/1erufTmtX3yySfExsayceNGLrrookZXq1Pdgz6SUqqNRowYwe7du+ssJdqnTx/27NnTbQOGa5RTXFxczXrhN954Y52V8JYvX17z+vzzz+f7779v9pqXX345l19+ufsbqzyue/5UK9UF5OTkYIypGR0F1LzW9CDqZKQBQ6k2qKioYPv27QQHB+Pl5VWz3zWSyLUetFInEw0YSrVBSUkJxcXFBAYG1tnv6s/ozulBlGqKWwKGMeYNY0ymMWZHE+XGGPOcMSbFGLPdGDPBHfUq1VkiIiLw9vamf//+dfa77ji66tLHSrWHu+4wFgGzmim/EBha/TUfeMlN9SrVKUpLSxsMqQVrfYywsDAyMzM7qWVKdRy3BAwRWQfkNnPIpcC/xLIJCDfGxLijbqU6wzfffAPQIGAAOBwO9u3bp3cZ6qTjqT6MvsCRWtup1fuU6pZyc62/j2oPqXXp06cPPj4+nm6SW3REevP77ruPESNGEB8fz+WXX05+fr5brqs8z1MBo7EcCQ3+/DLGzDfGJBpjErOysjzQLKXapnfv3g2G1Lr07duXqqoqKioqOqFl7uHO9OYzZsxgx44dbN++nWHDhvH444+75brK8zwVMFKBfrW2Y4G0+geJyKsiMlFEJkZHR3uoaUq1Xm5uLuHh4XWG1Lq4Rkilp6d7ullu48705jNnzsTb25ojPGXKFFJTUz3zJpTbeSpgLAOurx4tNQUoEJHu+9ukTmmlpaXs3bsXf3//Rst9fX0B2LVrV7vq6ULZzd2W3vyNN97gwgsvbF9jVKdxS2oQY8z7wDSghzEmFXgI8AEQkZeBFcBsIAUoBX7ljnqV6gzl5eXNLsXap08fgAZzNLozd6Q3f+yxx/D29ubaa6/1SJuV+7klYIjIvBOUC3CHO+pSqrP5+fkhIsTFxTVaHhwcjK+vb7sz1nah7ObtTm/+1ltvsXz5clavXt3+tO+q0+hMb6VaqbkRUmDNxQgJCenWfRgt0dL05itXruSJJ55g2bJlJ9Vd16lIA4ZSrbR+/Xqg8TkYLhUVFaSlNRjXcVJpaXrzO++8k6KiImbMmEFCQgK33XZbJ7RWuYOmN1eqlSorKwEaHVLrEhcXx+7du92z8p4HdUR6c9cdiOr+9A5DqVYKCgoiIiKi0SG1Ln379sVut1NaWurBlinVsTRgKNVKubm5zT6OguMjpI61dzyrUl2IBgylWiE/P5/09PQTrqjndDoBOHDggCeapZRHaMBQqhVcj5iaGiHlMmTIEICaGc5KnQw0YCjVCna7HYBBgwY1e1xwcDCBgYEUFRV5ollKeYQGDKVawTUH40R9GAABAQEn/VwMdWrRgKFUK2zatAlofkitS0lJCdnZ2R3dJLd67rnnGDlyJBEREfztb38D4NNPP213Xix1ctCAoVQr2Gw2/P39mx1S6zJq1CicTme3Wkjpn//8JytWrCAvL4/7778f0IChjtOAoVQrxcbGtui4mJgYHA5Ht+nHuO2229i/fz9z5szhmWee4c477+Tbb79l2bJl3HfffSQkJLBv377ObqbqRDqEQ6kWcjqd5OTktDhguO5CsrKymsxs25SVK1e6fQ5H7969mTVrVpPlL7/8MitXruS///1vzWzuM844gzlz5nDxxRe7bQU+1X3pHYZSLXT06FEqKyupqqpqUHblleDrC99+e3xfWVkZwEmfU0qdOvQOQ6kWcuVZ6tu37nL0DgesWmUtclR7Mblx48bx73//u2Yobms0dyegVGfRgKFUC7nW6B44cGCd/V5e1op3IlA7z2BQUBAhISEUFBR4spluFxIS0m36YVTH0kdSSrVQVlYWxpgmh9S6gkX1kyjAmund3fNJXXPNNTz11FOMHz9eO71PcXqHoVQL7dixA5vN1mBI7caN8Prr8Je/wM03g58ffPKJVVZYWNiiIbhdxcGDB4G6Kc3PPPNMHVarAL3DUKrFbDZbozO8jx6FlSut17NnwwUXHC+bNGkSVVVVOBwOD7VSqY6jdxhKtYCIUFpaSkJCQoOyuXOtL4A776xb1rNnT0SEgoKCFqUTUaorc8sdhjFmljFmjzEmxRhzfyPl/Y0x/zXGbDXGbDfGzHZHvUp5SmFhIZWVlURERJzw2MpKWLvW9dpanS8/P78jm6eUR7Q7YBhjvIAXgQuBUcA8Y8yoeof9CVgiIuOBa4B/trdepTzJta6Fa2htbQ89BPPnQ1KStf388zBtGqSkWIEG6HY5pZRqjDvuMCYDKSKyX0QqgQ+AS+sdI4BrqmsYoDOZVLdSXl4OHF/noraDB+Gjj2DGDMjJgV/8Aj7/HAYMgHPPPRdjjA5LVScFdwSMvsCRWtup1ftqWwD80hiTCqwA/scN9SrlMSUlJRhj6NevX4Oyt96y5mJkZcG990JMDFx8Mfj4gK+vL2FhYfpISp0U3BEwTCP76qfnnAcsEpFYYDbwtjGmQd3GmPnGmERjTGJWVpYbmqaUe2RkZBAWFtboENnKSitYREfDokXWrO+iInjmGfjuu0pEhIyMDM83upPFxcWd8FFcU8c88MAD9OvXj+Dg4I5qnmoDdwSMVKD2n12xNHzkdBOwBEBENgL+QI/6FxKRV0VkoohMjI6OdkPTlHKPw4cP1yzPWltWFkyfbr1+8EEYOhRuvRVKS+GBB+Drr70oKCjo9rO9Pe2SSy7hu+++6+xmqHrcETC+B4YaYwYaY3yxOrWX1TvmMDAdwBgzEitg6C2E6hZEBIfDwYABAxqUFRdb/RYAgwbBa6/B/v3wj3/Avn3wpz95MW3atCaTFnYlBw8eZMSIEdx8882MGTOGa6+9llWrVnHmmWcydOjQmg/w3NxcLrvsMuLj45kyZQrbt28HICcnh5kzZzJ+/HhuvfXWOuuAvPPOO0yePJmEhARuvfXWE85LmTJlCjExMR33ZlWbtHsehojYjTF3Al8BXsAbIrLTGPMIkCgiy4DfA68ZY+7Belx1o3SnVWXUKa2kpAS73c7gwYMblA0caM3wnjsX+vaFceOs2d7/+Adcc43Vn+Gaf5Gfn0+L75zvvhu2bXPn24CEBHj22WYPSUlJYenSpbz66qtMmjSJ9957j/Xr17Ns2TL++te/8umnn/LQQw8xfvx4Pv30U/7zn/9w/fXXs23bNh5++GHOOussHnzwQb744gteffVVAJKTk1m8eDEbNmzAx8eH22+/nXfffZfrr7/eve9PdTi3TNwTkRVYndm19z1Y6/Uu4Ex31KWUp2VmZgI0uaaFK3t5nz7W9yefhOXLrcBx8cVw8KBh8GDIy8trecDoJAMHDmTs2LEAjB49munTp2OMYezYsTVpQ9avX89HH30EwPnnn09OTg4FBQWsW7eOjz/+GICLLrqoZs7K6tWr2bJlC5MmTQKstO89e/b08DtT7qAzvZU6AdcHZUlJSYOy116Dp5+2RkRFRVn7IiLghResu46wMPD392Pw4FZO3jvBnUBH8fPzq3lts9lqtm02W02a9sYeDpjqzIvGNBwDIyLccMMNPP744x3RZOVBmktKqROoqqrCGMPw4cMblIlARYX16MlW67fpiivg0kth0yZ4/vkheHt7nzRDa8855xzeffddANasWUOPHj0IDQ2ts//LL78kLy8PgOnTp/Phhx/W3Knl5uZy6NChzmm8ahcNGEqdQGFhIREREYSEhDQomz8fRo06/jjKxRh48UVrFb7bbjMEB4efNAFjwYIFJCYmEh8fz/33389bb70FwEMPPcS6deuYMGECX3/9Nf379wdg1KhRPProo8ycOZP4+HhmzJhBenp6s3X84Q9/IDY2ltLSUmJjY1mwYEFHvy3VAqar9j1PnDhREhMTO7sZSvHCCy8QGBjIr3/960bLR42CkSOt2d71vfQS3H47TJmykZ///AfuueeOJutJTk5m5MiR7mq2OgU09jNjjNkiIhM7oj69w1CqGSJCbm4uubm5jZbPnAkHDlgjpBpz662QkGBn69YJZGa2fqlWpboSDRhKNaOkpAQRIT4+vtHy8HAoL2/4SMrFZoP33vPG6fRh+fLza3JSKdUdacBQqhmuO4tBgwY1Wv7Xv1rfmwoYYD2umj8/mx07xvLSS2VNH6hUF6cBQ6lmuNbj9vf3b7S8/hyMpkydmkh0dCYPPRSMJq5V3ZUGDKWakZqaCliTzepbuxbmzbNenyhgBASUM2fO5xQXe/Pww+5upVKeoQFDqWZUVVURFhZWM0S0toAA6FGdQrOpTm+Xyy+/nCFDshg/Ppv//KcDGqqUB2jAUKoZrvxPtWdAu0yeDD/7GQQGQhNZQ2oYYwgPD8fXN5+9ezuosV1MW9Obl5aWctFFFzFixAhGjx7N/fc3WPVZdRINGEo1Iy8vr9F0Fy5padbjqGYOAaz06EVFRXh55VFcbK2hoZp27733snv3brZu3cqGDRv48ssvO7tJCg0YSjXJbrdTUVFRk9Kivptugq++OnH/BVh3GDabjehoK6v/CSY6d4qukt48MDCQ8847D7BWLJwwYUJNX5LqXJp8UKkmFBcXAzBxYuOTZhMS4JNPWhYw+vXrx1lnncXWrSkAHD1qrfndlE7Kbt7l0pvn5+fz+eef89vf/tYd/wSqnTRgKNWEourxr7179260/M474Y9/PHGHt0t4eDghIYWAFQzOOMMtzXSrrpTe3G63M2/ePO66664m58Eoz9KAoVQTXJ2xTeVbKyiAsrKW3WEAbN68mdBQKwjt2tX8sZ2U3bxLpTefP38+Q4cO5e67727VearjaB+GUk3IyMgAaDSdR1HR8TuLlgaMoKAgAgLK8PV1EhDgrlZ6nifSm//pT3+ioKCAZzsrcqpGacBQqgk2mw0vL69G18FwOGDWLOt1SwPGlVdeSVBQIJGRZRw96saGelhHpzdPTU3lscceY9euXUyYMIGEhARef/11j7w31TxNb65UEz7++GOOHDnSZIfrW2/BjTfC3r0wZEjLrvn666/zl79cgq9vrwbzMTS9uWotTW+uVBeRkZHR5BwMh4Oau4SW3mEkJSWRlZVFYGAB1U+7lOpW3BIwjDGzjDF7jDEpxphGp2UaY35ujNlljNlpjHnPHfUq1ZHy8/MbXccb4Lnn4KGHrDW7AwNbdr3AwEBCQkKIisqhqkroojf3SjWp3aOkjDFewIvADCAV+N4Ys0xEdtU6Zijwv8CZIpJnjDnxmDqlOpGIICIkJCQ0Wj5hAgwcCD4+Lb/m0KFDycvLY/nyPMrLDXl5EBnZsN7mZpYr5dIZ3QnuuMOYDKSIyH4RqQQ+AC6td8wtwIsikgcgIo1PnVWqi6ioqKCqqorw8PBGy889F6KiWv44yqX2XIwdO+qW+fv7k5OT0ykfBKp7ERFycnKaTLvfUdwxD6MvcKTWdipwer1jhgEYYzYAXsACEVlZ/0LGmPnAfKDR7KBKeUp+fj5gZattTGmplUdq2rSWX9Nut/PFF18QGhoGwE8/wTnnHC+PjY0lNTWVrKystjZbnUL8/f2JjY31aJ3uCBiN3T/X/xPJGxgKTANigW+MMWNEJL/OSSKvAq+CNUrKDW1Tqk1ck/aaWlJ14kQ4cqR1dxje3t4MGDCAQ4cOAw0TFvr4+DBw4MA2tVcpT3DHI6lUoF+t7VggrZFjPhORKhE5AOzBCiBKdUmuWc0TJkxotPzGG0Gk5WlBXK644gqio627lu48F0OdmtwRML4HhhpjBhpjfIFrgGX1jvkUOA/AGNMD6xHVfjfUrVSHKCy0+hma6sOYOdP63to+DIAePUIJDCzn66/b2jqlOke7A4aI2IE7ga+AZGCJiOw0xjxijJlTfdhXQI4xZhfwX+A+Eclpb91KdZTU1FS8vb2x2Rr+ilRWWv0P0PqAsWHDBjIzMwkOLkK7KlR345bkgyKyAlhRb9+DtV4L8LvqL6W6vLy8PBwOR6MBY8sWuPpq63VrA0bv3r3p1asXEREFBAT0oPEuQKW6Jp3prVQjfH19GTRoEObee+Gpp+qU9esHF1xgvW4i83mTBg8ezNixYwkOLtA+DNXtaMBQqhFFRUWEhITAv/4Fjz8OFRU1ZbGx1uJHPXuCr2/rrx0WFkZISCFZWYYmJpIr1SVpwFCqHqfTSVFREY7cXMjOhrw8qLWmdFYWHDrUtg5vsFJ/u9bFSElxR4uV8gwNGErV41ppLyCt1ujwd96pefnHP8J//tP2gDFp0qSa2d7Vq8Aq1S1owFCqHtda3vHBwdaOc8+Fzz+H6tnfN95oJRxsa8A4++yziYwsA3QuhupeNGAoVY/rDiPEtSrcb35jjaX98EMApsuqEDcAACAASURBVE6FwsK2Bww4fu7Spe1pqVKepQFDqXoOHjwIgM+WLdaOa66B0NCax1LffNO2Wd4ue/fupbw8DW9vu3Z6q25FA4ZS9RQUFADgW1kJI0fC22/DL34Ba9di33+Y6dOt49p6hxETE0P//v0ICysiLMxNjVbKAzRgKFWPn58foaGheKWmwqRJ8Mtfwn33AWDef49777WOa2vACAkJYdiwYQQHF3DkiNNNrVaq42nAUKqeoqIiwv38IDUVQkLA6YRjxyAhAa/33mZgnJVIuT19GNbqe0Vs2aJJmVX3oQFDqXrS0tLwcw2p3bTJ6rCYPRuCgmDXLn5cnYGXF0RHt72O7du3ExpaRFWVweFwT7uV6mgaMJSqp7Kyksi8PGvjd78DLy9YsQIWLsRp8ybxkyNER1u72+r0008nJKQQh8PmGq2rVJenAUOpWiorK3E6nQxzJR10JY064wwYPpzy82cT7l1Cv9j2PUoaMWIE4eHWECmdi6G6Cw0YStXimoMRmppqzc5zBY7ycnjzTQLPn8IxexR9/dqXm9zpdBIVZS3S9K9/tetSSnmMBgylajlUPVnPd/dua+Fu15//Xl7wP//DD9ttpNGXPtk/tqueyspKvL3TgYZLtSrVVWnAUKoWV1oQ/5ISmDEDhg2zCnx8YOdO5u/9A7lE0ufAeiugtFFAQAAJCVaveUhIu5utlEdowFCqFi8vL4zDgU96ujUHo3b+8gEDePBB63agT+VBWFZ/JeLWGTQolsDAEg4f1rkYqnvQgKFULYWFhUSXlWHsdqiqqluYl0fkkpcA6Nujsk4G27bw9fUlNLSQjz6yt+s6SnmKBgylaklJSSE8p3q5+e+/r1NWKgEsX2I9hupzyWmwciXtWZg7NzeXkJAi/P31DkN1DxowlKrFbrfTs9Baq4K//71OWUqqP09U/R6APr+eBQ4HLF7c5roSEhIIDS2ioqIdEzqU8iC3BAxjzCxjzB5jTIox5v5mjptrjBFjzER31KuUuxlj6F9VBf7+MH58nbIhQ+Daa8HPDyLOHAXx8e16LBUTE0NoaBF5eT66kJLqFtodMIwxXsCLwIXAKGCeMWZUI8eFAHcBm9tbp1IdQUQoKioiLC0NwsOt5VlrCQwEKS6hj9cxzDfrrKSEmzfD3r1tqs9msxEeXg7AkiXtbr5SHc4ddxiTgRQR2S8ilcAHwKWNHPcX4Emg3A11KuV2WVlZOJ1OAo4csZIN1svZsW0b/LjXn77OVCuYzJtnTaJ499021xkcXJ1K3fcEByrVBbgjYPQFjtTaTq3eV8MYMx7oJyLLm7uQMWa+MSbRGJOY1Y7ORKXaorCwEEQIzM+3VtkbPLhO+RtvwM7dXvSZMxGuuAJiY+G886zHUtK2VCFjxoQCGjBU9+COgNHYPNWa3x5jjA14Bvj9iS4kIq+KyEQRmRjdnlSgSrWB0+kkuKgIr4oKGDOmQXbBP//Z6tqok9b8l7+EffusR1NtMGxYEAA//aQjpVTX546AkQr0q7UdC6TV2g4BxgBrjDEHgSnAMu34Vl1NUVHR8Sy16ekNyv39rcndfUw6jBoF69ZZdxr+/m3u/O7Rwwtv7ypee62iPU1XyiPcETC+B4YaYwYaY3yBa4CaKbAiUiAiPUQkTkTigE3AHBFJdEPdSrlNUlISkbm51saPDXNFvfKK9b3PiFAYONC6AwkLgzlz4IMPGk70awFfXx9CQoqIjCxrT9OV8oh2BwwRsQN3Al8BycASEdlpjHnEGDOnvddXylO8vLzoUVBgBYJ6KWRF4P/+z3rdZ2gQfPEFnHmmteOXv4ScHPjqq1bXGRcXR2hoITabdmKors8t8zBEZIWIDBORwSLyWPW+B0WkQbIdEZmmdxeqK/L29qZ3cTHExUFoaINy1zy+vq4hHZWV1uS9Cy6AqCh4//1W1xkZGUlISBFpaV5U6FMp1cXpTG+lqhUVFRGelQV2e4NVjYyxlsQAiIkB1qyxHkclJlpDnM46C374odV1ent7ExZWRmamL9u2tf89KNWRNGAohTVpLyMjg+CMDEhNbdAfsX+/lToqOLg6HfmoUdbQ24gI64CRIyElpU39GOHhJTidXprmXHV5GjCUAioqKvAvLcWvvByeeAIGDKhTvnEj/Pe/UDPau2dPePrp4+tljBxp3ZmkpLS67qFD/QDrdKW6Mg0YSgFlZWXHh9QOHdpgGbyrr4aJE+vFERErQDgc1h0HQHJyq+seMMAbgC1b2rdOuFIdTQOGUkB+fv7xIbWN9EV4e1vZQGJja+1cssQKLj/+CCNGWPvaETCefbaw1ecq5UkaMJTCmoMR4QoYjUzaW7LE6tqoM8v73HPh5ZetYVPBwdCvX5sCxpAh1mzvvn1z29J0pTzGu7MboFRX4OfnR2RuLs4+fbC5ZujV8vzzVh9DnYDRuzfceuvx7ZEj2xQwevWKICiomMDAsDa0XCnP0TsMpbAm7UXm52OGDGm0/Pnnre91AgZAQQF8+aXVnzFyJOzeDc7W5YUKCwsjJKSIlBTt+FZdmwYMpYC8vDyicnMxhw/DoUMNyjMyrO8NAsaHH8Ls2VagGDnSSjZ15EiD85vj5+dHaGgJBw54cfCgRgzVdWnAUAo4uHMnQUVFUFTUINd4WdnxO4y+feudeMkl1njbgQOtgAFteiwVGVlGVZUPgYElbWi9Up6hAUMpoJdrjdR//rN6KvdxGRmwYoX1ul6RNR9j2jQrY60rYOza1er6+/WDsrJAoqK0H0N1XRow1ClPRKwZ3tBg0SSwUkvdequVLsrPr5EL7N1rBZoePayD2nCH4bpzWbtW52KorksDhjrlFRYWEupa4XFZg3yZgDXStkH/hcu//w133AGHD1sT+NoxF+PRR4+e4EilOo8GDHXKS05OJjI3F7ufX4MZ3gCrVsF331mjaBs1b541SWPAgONDa1u5ZOugQdaty/DhGa1tvlIeowFDnfL8/f2JyMvDPnw4LFjQoPy77+DYsUY6vF0iIo4XjhwJubnQyjXphw4NBCAyMvYERyrVeTRgqFOeiBCZm4txJRKs549/tG48Ypv7LN+wAe69t80pQgYMCMPbu4otWxytOk8pT9KAoU55WUePElZQgM+mTVYe83oyM625eE32YYCVT2rhwuPpbFsZMAIC/AkJKWLnTkNy8uFWnauUp2jAUKe81PXrsYlgCwmxFkWq54EHrO/NBowbb7QeQ02YAEFBrQ4YxhgiI8vw83Pg79+qU5XyGA0Y6pTX17WU3muvWcNi61mzxvrebMDw97dS2hpjPZZqw0ipnj2rKCkJYeDA/q0+VylP0IChTnnBx45ZL5rII/WHP1jfmw0YAF98Ya3vPXx4mwJG794O8vODSEpq9alKeYRbAoYxZpYxZo8xJsUYc38j5b8zxuwyxmw3xqw2xgxo7DpKeVpVVRV+qanYvb3hH/9o9Ji0NLDZoFevE1ysogJycqze8dRUKGzd+hZ9+xrsdm/++MfEVp2nlKe0O2AYY7yAF4ELgVHAPGPMqHqHbQUmikg88CHwZHvrVcodMjIyCMnKojQsrMGyrGDdKLz3nvWkyvtEiwFccQUkJsLUqdb27t2taktcnFXBOefkI62cx6GUJ7jjDmMykCIi+0WkEvgAuLT2ASLyXxEprd7cBOhgc9Ul+Pr6EpmbS9WYMdZs7XoyM+HoUSvrR4u1cWjt4MFWb3d09AhMIxMIleps7ggYfYHa+ZxTq/c15Sbgy8YKjDHzjTGJxpjErFZOfFKqLSrLyojIy2tyHYxzz7W6JJoobmjJEjj/fPDxaXXAGDbMWnlv9eoqvcNQXZI7AkZjfwo1+tNujPklMBF4qrFyEXlVRCaKyMRo13h2pTpQ5tateDsc2FauhJSURo9JS2tBh7dLv35WwIiLa8MdhjXbOzHRwZo1a1t1rlKe4I6AkQr0q7UdC6TVP8gY8zPgAWCOiFS4oV6l2u3oN98A4Dd2rJWqvJ6nnrKmVzSZFqS+qVPhnXcgPr7VAcPPzxAcXEp0dAV9+7Y0QinlOe4IGN8DQ40xA40xvsA1QJ2Un8aY8cArWMEi0w11KuUW/Susv138X3oJQkMblO/da31v8R2GS1wc7NtnjZxqhaioMkpKQhjWRJoSpTpTuwOGiNiBO4GvgGRgiYjsNMY8YoyZU33YU0AwsNQYs80Y03gOaaU8LCAtDYeXF6Zfv0bLf/Ur63urAsbbb1tDdJ3O4xGnhXr0qCIry58fftB+DNX1uGUehoisEJFhIjJYRB6r3vegiCyrfv0zEeklIgnVX3Oav6JSHU9E8Dp4kKKQELj99gbllZVW/wW0MmBMnQq33Wa9buVjqehoJ4WFQTz00JcUtnIeh1IdTWd6q1NWaWkpgenpFEdFWR3V9cyfD3/+s/W6VQFjyBDrDsOYVgeMQYMMpaVBzJ/vj5eXV6vOVaqjacBQpyx/Pz8i8/KoHDECrr66QfnEidb8Cx+fRlNMNc/X11oAvJXre8fF+QAQEjKc4ODgVlaqVMfSgKFOWfb0dPwqKpBBgxotv/NOq++6Tx8rNUirLF5sPc/64YdWnTZokDV574MPSikuLm5lpUp1LA0Y6pSVvn49AGb58jqd03Y7fPWV1WfdqjkYtV1wAVx8sbXOt6PliyK5Vt5LTCxn6dKlbahYqY6jAUOdsjK//RYAn1mz6iyn98UXMGsWfPYZbN8OAwe24eI9esBll1nDag8ebPFp/fpZv5JDhlQxbdq0NlSsVMfRgKFOWbEVFQgQ+Kc/QUBAzf4LL4SlS6G42Jq0d/31bawgJsb6vm1bi08JDwdfXzt5eYEMbFOkUqrjaMBQpyzvgwcpCAsjJDy8zn5fX5g7F156CYYOhRkz2lhBSYn1feXKFp9iDERElJOa6sP27RWsW7eOoqKiNjZAKffSgKFOWWb/fvLDw/G95ZaafU89BW+9BVu2wMaNVgLbVnd4u1x8sTW8ym5v1WmRkVUUFQWwe3cZa9asYW8rJ/8p1VFOlOFfqZOSiBCYnk76iBFwww3V+2DZMhg82FqWNSjIWqq7zQICYOzYVq+LMWCAsGlTCDNnwqxZdxPaSMoSpTqD3mGoU1NhIUGlpVQMHgwzZwLW46B16+CRR+D9962+i7CwdtbTty9s3QquZWBboF8/G0VFIeTl5dcEC6fT2c6GqJNZaWmpR1LJaMDoZpxOZ80z7crKSsrKyjq5Rd2T2b8fAImJAacTp9NKBWKMFSwqKhpdT6n1+vSxLvaf/7T4lLg4HxwOb15/3fp/TkxM5MUXX8TRiuG56tTywQcf8MEHH1BaWnrig9tBA0Y388UXX/DGG29QUVHB0qVLefvtt/WvzzbIT7TWzZZNm+DAAVatgv79rZuBf/7TyhQyerQbKnL1mPfu3eJTXCvvpadbmW4jIiLo378/Fa3MfKtODSLCaaedxtixY/npp586tC7tw+hmTjvtNHr27Imfnx9TpkyhpKQEW5t7ZTuf0+nEGOPxJUnzt2whHJCLL4b+/YnIhenTrTWUDh+GZ591U0WuqJOc3Gi+qsb072/lkOrbtwqAwYMHM3jwYDc1SJ1sjDGMGzfOI3V130+aU0hFRQW7qztO+/Tpw+mnnw5YHyTx8fEAHDlyhJycnE5rY1tt3LiRb6oXMfKk6MJCioOCCJs7F3x8mDQJ3n0XXnnFWjTvkkvcVFFMDAQHw9NPt7jz27VY04EDlRw+7OD2262nWnl5eWRkZLipYepkkJmZSVJSksceV2rA6AbWrVvHhx9+SEFBQaPlTqeTTz/9lM8//7xbrKFQXFxc88FXUFDA7t27Pd5us28fuZGRhFZVsWYN5OZaeQJXr4bf/Aa83XXvbYyVvfboUUhNbdEpMTFgjJCaKixZ8h6LF1ewY4fw9ttv89VXX7mpYepkkJSUxIoVK6iqqvJIfaarfsBMnDhREqufM5/qqqqqSEtLY8CAAU0ek5WVha+vL2HtHtbTsUSE1157DafTya233kp5eTl+fn4nfqy2b5/VDxAU1P5G/PQTjjFj2DpuHANGT2D4sleYNQsiImDhQjhyBNy6pPyvfmVN3ktPb/EpvXvDGWfkkJDwIlFRvbjhhmvJyckhIiJCh9mqGiJCTk4OPXr0qNlnjNkiIhM7or5T5g5DREhKSqLENfu2i8vPz+ezzz7Dbrfj4+PTbLAAiI6OJiwsDBFh9erVHGxF/iJPKCsrQ0QwxjB79myuvPJKjDEEBARgs9lwOp1Nd96vXw+jRlkJ/Vo5Ca4BhwNuuIEqHx/WnnMOUbffxJo1cPfd1oS9a65xc7AAGDnSGlabn9/iU/r2hdLSKH7xi19QUJDDwoULWbcuiCVLOi5YuP79i4uLefHFF2seg5aXl7Nlyxbyq9svIt3iTvZk5/p9qh0sOtopEzAKCgr4/PPP2bx5MwA5OTlUVlZ2cquadvToUfbs2UNubm6rzqusrGTPnj0dPlqiNQoKCnjxxRdr/u1jY2OJrvWpXFRUxAsvvMC2xnIupaRYSfzCwmDDBnj44WbqgSuugEcftSbh1VZSUsLKlSupeOwx2LSJpFtuQWJisE2eTHw8bN5sZfK48063vOW6Ro60vo8fb030aIHTTrMy5j733BCuvvqG6v/XN/j004MsWfKhW/urRIT333+fr7/+GrACR8+ePQmozq+Vk5PD8uXLax4jpqen87e//Y2u/gSgsrKS9PT0k3J0WVFREc8//zz7q4eHe8pJHzBc8xTCw8P59a9/zbRp03A6nTXjlrsaV+fV6NGjueuuu+jZs2erzvfz8+PXv/41M6qHc3bm2H3XX6GhoaGMGzeuyWR6wcHBDBgwgPB6OZ3Iy4OLLrI+/TdssB7tPPZYo3MacnPhZz+DTz6xVsl76KG65RUVFRz+4gt8Hn0UrrySvWPG4FXhw9132snKghdegClTrEWT3M4VMIKCWpxn5Pnn4be/tb5feWVfzjnnFqKj/Zk69V32708hMzOzVU0QaZhl3bXehjGGqKgo63Hmli2ErljBVXPn1tzVxsTEcPfdd9f8//n7+zNp0iQSEhIAsLf3rs+NnE5nzR+CmZmZvPrqqxw+fBiwBoY8/vjjNXff6enpvP/++2RlZXVWc9usvLyc8PDwhr8zHc11e9nVvk477TRpr8OHD8vjjz8uKSkpDcoOHjwohw4dEhERp9MpDoejTXVUVFRIVVWViIg4HA4pLy8Xu91ec92Kioqaa1dUVEhqaqqUlZWJiEh2drasXLlScnNzRURk37598sQTT0haWlqb2lJfSUmJvPqPf8iWxES3XK81Du7bJy+//LIUFxc3e9zhwyLl5Y0UVFSITJsm4uMjsm6dZGZmStGxYyLDh4vExIhkZtYcmpkpMm6ciK+vyGefidx0kwiI3HLLIfnyyy+tgyorxTFunDijo2Xfpk3y1COPyGPX/K/07GGXzz6zjn/nHTf+A9Rmt4v4+Ynce2/Lz0lMFHnmGfl8QaL0iLRLYKBTXnqpXF555VV5+OGH5Xe/S5TqH99mVVVZ72vECJH+/UW+/dban5SUJH/5y18kJyfH2uFwiDz5pIi3t/WPMWuWyNGjjV6zqEjkV78SGT1a5O237fLSS6/ImjVrWv7eGlFZWSn5+fntuobdbpdXXnlFli9fLiIipaWlkpycXPMzmJ2dLV9++WXN79uhQ4fk5ZdfloyMDBERKSoqqvnd7c6AROmgz2X3XARmAXuAFOD+Rsr9gMXV5ZuBuBNd0x0Bo6KiQj799FMpKCho9rhvv/1WXn/99ZoPcnE6RdLSrO+NcAWAY8eOyYIFC+SHH34QEesHcsGCBbJ9+3YREcnIyJAFCxbIzp07RUQkNTVVFixYIHv27BERK6A99thjsn//fhERyczMlOXLl7foF6ey8gQHJCWJ48orxWmMVAwbJvLKKyIlJY0eevCgyD/+IfLJJ9aHQUs4nU7Jz8+Xkupr5qelyYp77pG0G2+UwlGjxOHlJT+NHi0HtmxpcO7OnSILFoiMGmX9BA4fbn0+VlRUyPbt28XpcFifSCDpf/+7vPnmm7JgwQJ55JFHZOUTT4jD11ecs2eLOBySlmZdx99f5KuvrP8bh0Pk+uuta8+bt83a9+c/i4CsmD9fFixYIP946inZ/8EHUlpkl0suEenZs4nA5S7x8SIXXWRVUlHR/LG7dokEBVlvAOQoMXKebY2AyDUDNsjqi26Q96++Wl6/b6E4m/hBqKwUefNNkaFDrcuMGSMyYIBDvL2d8vTTIoWFRfL1119b/39ZWSKzZ1sHzp0r8uyzIgEBIpGRIkuW1LluUpIVfIwRGTLEOmXgwAJ59dXUpn5dTqigoECefPJJ+e6770REpKqqquZDvCVq/76sW7dOkjdvFlm4UOTuu0XWrm3y91hEZP9+kaVLRdLTnfLGG2/IokWLxNnWN+Ihhw8fPv5Z1YiODBjtHiVljPECfgJmAKnA98A8EdlV65jbgXgRuc0Ycw1wuYg0XES5luZGSdnLyyk5cIDy/fupOnQI++HDSFERVePHUz5uHD8dOsSsWbMICgqiuNiQkmI9Ct+711rfYMIEOOcccPUj79y5k5TkZOaEhyOffALLlmFLTcXRsydVZ59N1VlnYT/3XOyxsSxZsoTevQdhs53FunXebNhQwvDhAUyeHMiwYZU4HNuZNCmO6OhoSkpK2LZtG8OHD6dHjx5UVFRw6NAh+vbtS5DTSdn27ZQmJVG5cydm7158Dh3CLzub7LFjKbzsMnxmzKBn797YbBFs22YjMdHKopqYaA0amjQJfvlLq6O25snV1q3II49gPv0UZ1AQ2RddRPC2bQT+9BOVISHkX3UVJTfeiCMmlm++Cea990JYvdoPEWvinK8vTJtmJVq96CKovXppcXExdrud8PBwyvLyWHzPPZwHRO/YgX9SEja7HYfNRlqfPhT068fI776jKCSE737/e5wTr2fTpn4sXWrYtcsabXr22VYap5degowMmD8/lR493uB3FaWEPfEEibNn88XkyYSEhDBlyhSKi4tJSkpi9Jo1zF6xgm9u+BO/Wv8wx47ZWL4czjijkjfffJNLL72U6OjeXHed8P77ht/9fDNPfXgGP44Zw8bf/Ibx48/Az280CQleHDhgJRt84AH4y19a//PfYtdcYz1Wy82Ft9+2OlsaU1wMkydDdrY1xjcvD5KTcexI5m9fjuOhfdcxgEN8wDVMIhGHtze2s87C3H8/zJxJZZXhX/+Cv/4VDhyAhATrEd0llzj429/+yeefX8H33/flssvgjTcgYsc3MG+eVd8zz8Btt1n/OXv2wHXXwfff45w3j9Inn+K1peH83/8FEhLi4O9/P8bkyaWsXh3N3/8exsGDNs4+G265ZT9TpzoZMniw9X7/3/+D776DXr2s8cLVX5VRUWR4eeHdvz95fn5sPnCAnn36EB0dTV5eEf/+91bOPPMKAgLiyMw01V/Wz4mINYH+wgshNXUbn3/+ObfdcgvRSUnw1lvIJ59gyssRb2+M3Y5j0CCqrr2WqnnzMP36UVxs45NPvHn3XW82bLAeERojjB1bzrnnFnLLLb0YPVqoqqrEz8/Pepb3ww/WuibBwdZwuuovZ2goZf7+lNntlJaWUlpaSmVlJf7+/gQFBREYGEhQUBA2mw9lZYaSEmq+7HYYNsy6ZB1lZcjGjTj/+19k7VpsO3fiGDsWx7nnYp82jZc2b6ZfXBxXXXUVWVk21q0zrF1rNW/oUHjrrY4bJeWOgDEVWCAiF1Rv/y+AiDxe65ivqo/ZaIzxBo4B0dJM5T2DBsuv+9zJ6F65BFUVYMvKwbcgj+DyYgJLS/DGgUGw4cQgVODHfgaxxwxne2A8u71GcKhyAIXldYeZ+vo6qKy0ZtKGBRUyKSaZc+z/YU7aYuIrk7B7e7Nv8GAO9+9PTHo6Aw8cwLekks2czhf+s/m390y2l8Zjd/pgszno0SObwsJQysuPL8ATGlpGTGgmvX3SODuhgLFhqUTs2UzAsYP088sj8OhRAosK67SrICSUvLAoKqN7k7XDxo9VY/nWZyrfeZ9Oalm/muN69y4nPr4Kjmaw42gP0vLDsdmEMb32cINjIbdlvojNz8mmKVPYfPrplAcGggj9Dx3i9M2bCUvO5w1u4iWf35Be1Yfg4CImTNhKQsI2CgpCOXBgNHv2DCUjw3o2GhdXxvTzSghMXMdpVV8yfmQOI7KyMBs34+OowmkM6TExHBgwiMz+8Qy84xoGxcfzt/M2MaRvNgfWpPFJ8Wx2MRqDMC6hiOuv9yf5pQ1cOc+PCx4+g8w0O2ePzOKnwhjG9zvIx0emUTQmhJVzr2fGpTMY1HMgC362kctvjmLqXeP54ZtE8q9fwK2pL3HMtzf3/W4t5ZvzmXVxLD8F7iC+92Sef6CQ869L4/m3x/HjrtE87vNHBj9wNnMfvIiNH6Vz5lV9eOreDDLoxTPPwKFDxyfMdYiHH7a+7rgDbr7Z6s944QWrl33IEGtS38svw/79sHw5vPYa7NgB99xj5SxJSoJFi/h22v8x73+iSDsq3DFoIVdG/4uxO3YQUFTGs6H38Hf7fWSXRhHbN52fnbma/sMPI/gS4l+Js7ycKp9gvlk/iVWrzqGXbyafls9hUOQhVl/3C3JCQykOiqGkVIjwL6CqpIRJG74lYc0WbvZeyEf2uQwatI8rrviE4ODjow3tdhtbt05g3dpzKCoOYULPrTxn/sCZGasoDwxh3+Cp9Ko6jHdeHr7FJQSWWH0nFfiylfFsYgqbOJ0fbWM5ZmLIc0QgjXSv+njbCQquwm73orjYF5txMrB3GjN8VnJT9iucVppIuX8AO8aOIWncODJ79mRkcjIJW7fS/+BhvmYG/wy+g69KL6DK6UtUVDbjxiUxYMAhDh2M46efhpN61Poh6BGczZnh67lYvuLq7PcJqWg4B0qAY/QmhSEke41gj9cIfjJDOSL9KTQhlBJEqQRSPCzJDwAAFDJJREFU7vCnyuHbxA+G0LNHIcMDdjPGZxeTCtdzXu4qBjgPI8ZwrHdvMnr1IiY9nd4ZGaTSl3/7zODLoFlsrDqD1BLrs8HXp5LeMRnk5ERSUhLcpQPGXGCWiNxcvX0dcLqI3FnrmB3Vx6RWb++rPia73rXmA/OtrdNOg7aNwujtnckwezJDSGGwbR+9ehbi1buC8qE2SvsE4Z3kIHtnD3Zlj2C9nE0GVp6fIO9ixozNJipmD70rvfnqx2GERdhISelJpd0Hg5Pxtm1Md67iPP5LH7+j9J09BEdJKfuSKthf1Ic0R3/2VAxhJ6PZyWiKCWnTewDobdKZKpuYyPf0DTuKjINjU2OoqLU6nElycPS/ffgyfzZH6E+gdxlnn52NI/sI4bYSXl87GRHDrJGHSS2JJqM0GrvDxvms5g5e5GdjD1H4y6v5/Ucj6Nm7iEsur8IrOZmUVfmkHBvC5uzT2VQ5hSp8CSePmXxNQGApOUE9yPHrgXdsf7z8I9i01kZkcCW9h4ZSVgZ7kh048cIY4ZzInVyV8xJjB+xi4zWnYw8J4YfvRjDBN5AHll5IYXYhl1/+PQPNYZZ+cwVO40XPXse469wMfvvBmeTuz///7Z17dFXVncc/vyQESCAFgoTwEGQCrUJLFMujFosRR6FaUBQUqDIOfagji1pdrXV01DJrBB1FdCqyVpWH1WpwLEQdsYAi5U0gQniG98MggSB5kJtwc3/zxz4Jl5CQS573tr/PWmfdc/fZZ9/vOXfv8zt779/em54p0cz4aTa/nDeEvy3OZ9ioBBIoZHHH8Wx4aAgFosSUt+K6tEEc2nmIvbkHiBYYs3kLTy4ezwfcwRN372HaOymcXr2N9WOm02/qjfSdfh/Dh8N779X5bwqN996DcePca2D//s5V+NZbISPDVbWWL3dVOp/PVXWuvdbF//xz5121aJGbPnf1ak516cvPRhzh/XXd+E5KPh07HmBXVg/yfIkMZg0PtplNyq0BjvuKWHn5CA5+1YWJA3YSs24VWf3vZe+GQu7cksEDR14gV5K5+ntZ/PrOv8LOHWzv9WP2HYjijt67ic3M5IPLHif9re9Q6GvD73mS6675iv+IeYy3/vnPkJHBi2l/ZufOIp7ovYC4pWv4S86tvHx2CgUkMKzL51wzehNtOxUjAaXsaEuyjv2Ag/u7UH48ipyTl+PXFgB0a3GMq8kiQb/hCn8OnTlGEl8T19KHv0MLvtX5Gwp6tGdXj96UxyjxfzvD/uwr+Kz4BjZzDQBdW+eS0N1HnyuL+E3sdOITWjFlz6MUFrbn8LYY8koTaccp7uEdJsS8Tc+ripjebiqtB/fiJ9ufJ+HsWfZlt+ZwXnc+LRvOUoZTQhzxUcVc3Wkz3fsc43BeEiUl7Tmd14qjvq6UlJ8rizGcJUm+5orYwySVHqUNRcRTTDzFxHGGOM5ArBLdBqJan6UsrjXHczuyv+gKvqQ/+zg3/Uu76G/o0f00eWUt+dGgYs78bTvrzw4l9xv3IpfAaYaykuv5gu/HZ3FFuxyODhjKkR1l3J3zblgbjLuAm6sYjIGq+nBQnG1enGCDMVBVa/QNTL3yuzr3f1bw9Zm2nDhdyLFjx+nT59uoCrt27Wbfvv3cdNPNqMKmTZs5dGgfDz00kn79Wrsq3unTrlCuWOG2zMzz3UR69oTRo9FRo8npdB1frI5h5Urn9Rg8hKFfPyUtTUhLc81Y7dv63Qx1y5e7JoONG53LZ+fOrtrduXPlviZ15kCgG1l5ndhxLIGo1glER0df9H5GR7shBwMGeOMBvvrKzVkxbx5s24a2bEnJ8OGcGDiQxCVLiF+9Gu3YEX3kUVamPsxb/xtHerq7/ORk90z67DPXJJeYCPffD5MnB9iyej6Ddu+m+/vvw+7d7sfKy12TiYf/sss42LYt/tTr+SJ2JB8d7s+qnd3xl8cQHx9N69ZRtGpFjdvgwa7lJbmzwqxZ8Nhj+JOSWD1lCivLyvD7/XTp0oUTJ07Q+uuv+cWbb3KkTT/GtltM1vZExo51zVUdOpy7P9u3uzmf/H5YOmMT/Sd/n5Lbb2fRmDGUlpVx4MABWrZsyYABA/gBEH/zzZT9yy+4M+81MjLci/vkyS6tP/7R7a9Y4f7bRmXrVvje9+Dtt10TUFXWr4cf/tC1tWRk1OpNpQqvv65MnaqUlrrmoEk/3c+3Pv41t2VvJXbPHrRPH+R3v4Px46GFezCzfDlMmADffEP+c3O4b+lEPvxQuOsud28qxn2qukkYH3nEZY133vQxdPFjrlZ05ZXw1luubffLL12z09tvu/lLRowg//5HmbHhBl6aGaC8HAYPDrBjRzT5+e6a4uOdPRw0yOWRQYPcpL6VFBa6dDMzz207d0LVcToDBsC993L0R+P5cG1HMjJckfT5oG1bN8XL9u1u1P6IEW7Zk1tHlNNy5VLXHveXv7gpii+7zLVVgyvDN94IaWmc+cFwlu3pwaJFARYuLOH06XhiYqB79zJiYg5w3XWdGDCgHUlJBZSWbuPHP+5H+/ZBL4hlZXDkCLp/P4F9+4g6dAg5eJDyvXvh4EGi8/PRq69mR1ISHUaPJm74HWzZ14asLCq37Gx3Wzt0gNTUAnr3zuX++3sx4JoYog/shaVL3UVXNF9GRSGBQKMZjIbo8B4CLAn6/jjweJU4S4Ah3n4McALPWNW0de3aVX1eL+SqVav06aefruzo2bJli6anp1d6NBw/flyPHDlSYyeQqqoWFKh+8onqyy+7nruLdGwdPOiiXkK/W+MTCLie4YcfVk1MdL2NSUmqL7ygWsUTqaTEdeSNGuWcjIYOVf3Tn1y4qqrP59P58+frxo0bnXfMxx+rjhun+vOfa2DWLD3z8ceqeXnq9/t1xYoVF+1guyTWrVPt2VO1RQstnTFDV69apbNnz9YP5s7Vsj59VNu1U92xQz/66BMdP/5LjYkJaNeuqkuXutOzslQ7dlTt3Nl1nKuq6rPPqoJ+NHasFhQUaEFBgcs3RUWuV7ZnT9WCAvX5nOOPiOrcue52pqa6vugm6eP0+VSjolSffPLCYydOOBemHj1UK7yWQiAQCOhTTy3QF19cq6rOS+jw4cMa8PtVFy50FwjuHrz2mupTT7kbcOWVqlu3qqr7+6dPV42Odrdr82bVU6dUx4xxp44c6frEK1myRLVLF+dNNXCgixQXp/rAA6o7dpynLysrT8eMOa6pqc6H4cUXizQry3luXTJFRaqrVqnOmqX6zDOV+qtSXKy6eLHqz37mnOxmzrxIOT5xQvWVV1QnTHDPhezsajNDSUmJpqe/r6tXH9CzZ135yczM1IKCgjpcyPmcPn1aX331Vc3Ozq72+Nmz7nlUXq66bNkynT17dvWd8n6/ak6O6pkz4e0l5RmAfcAVQCzwJdC3SpyHgNne/t3Ae7Wlm5KSUukO5/P5tKSkJOy9F5qM0lLnH1mD11MwNd2y4HuZm5urpZ7nzuLFi/WVV16pdBVucPLzVUePdllv1CjnE3vTTe4BtGzZedo2bnQeVOBcZdu3V+3WTXX37qD0/H7VYcM0EB+vunPnufCHH3YnfvZZZdCZM6rDh7tn5tSp7vCcOY1zmdWSkuK8kIIpL3eWLDZWdcOGS04yLy+v0j38AgIB1Q8/VB00qNLjSidNuuAFQ1V15UrVrl2d92+3bu7veP55J+8CTp5UnTjRuafNmOH+U6POBLvy7t27t9Ltt7a4NRHWBsPpYyTOU2ov8IQX9izwE2+/FZCOc6tdD/SqLc2GcKs1aqe0tFSff/55TU9PV1XVPXv26Lp16+o8LiUkAgH36teihXPfBNU33rggWlFRkZ48WaIPPnjuRdnzQK7UevLkSdUjR1ytKzXVvckvW+ZOmDLlgjSLi1VvuMEdbteu2mdn43HbbW7wQjBeDUn/8IfG+91AwN2TRYsuGu34cVej6NVLdc2axpNjVI/f79eXXnpJFyxYcF54IBDQwlD93TUCDEZjbGYwmo6cnJxzA7iaknXr3ACBZ5654FBhYaFOmzatckDYhg3njdVTv9+vM2fO1Pnz57uAjAyXnSdPdk07vXvXWAMrKlIdO9a1bjQpv/mNM5IVtbdPP3XVnQkTmqhdLDTCSMo/HPn5+Xrq1ClVdYMZS0tLdf/+/frss89WjteqjcY0GDZbrRG2rFu3jpSUFBITE6s9XjEZXuX0CL/6lVv5KCrKOTwMGdJUUkNj3jyYNMmNcWjd2nUaJyW5iawaYhZe4++KjIwMDh8+zLhx49i0aRPDhg2jRYXzwkVozNlqbcU9I2ypWCiqKrm5uSQnJ184j85zz7kRjUOHhp+xgHNzSm3Z4hZU8vlg4UIzFka19O3blw4dOpCYmFg5N1xzYzUMI6w5deoUmZmZpKWlERUVxcGDB5k7dy6333575WqDEUNBgfNbTU52a2O8+y6MHdvcqoy/M2w9DOMflmPHjrF27VpyvcWHunfvzi233MJVV13VzMrqQEKCG06emwtTppixMCIOq2EYYY2qUlxcTFxcHH6/n9jYmqZYiBAmTnRLtX76qZu4yzAaGKthGP+wiAht2rRh/fr1vP766/h8vuaWVD8WLHCjrc1YGBGIdXobEUFxcTEpKSm0atWquaXUDxG3GUYEYgbDiAjS0tIQe9AaRrNiTVJGRGDGwjCaHzMYhmEYRkiYwTAMwzBCwgyGYRiGERJmMAzDMIyQMINhGIZhhIQZDMMwDCMkzGAYhmEYIWEGwzAMwwgJMxiGYRhGSJjBMAzDMELCDIZhGIYREvUyGCLSQUT+KiI53mf7auKkisgaEdkmIltEZFx9ftMwDMNoHupbw/gtsExVewPLvO9VOQPcq6p9gVuAmSLSrpp4hmEYRhhTX4MxCpjn7c8DRleNoKq7VTXH2/8KOA5cVs/fNQzDMJqY+hqMJFXNBfA+O10ssogMBGKBvTUc/7mIbBSRjXl5efWUZhiGYTQktS6gJCJLgc7VHHriUn5IRJKBBcB9qhqoLo6qzgHmgFvT+1LSNwzDMBqXWg2Gqg6v6ZiIfC0iyaqa6xmE4zXESwA+Av5dVdfWWa1hGIbRbNS3SWoxcJ+3fx+wqGoEEYkFPgDmq2p6PX/PMAzDaCZEte4tPyKSCLwHXA4cAu5S1XwRuRb4papOFpGJwJvAtqBTJ6lqVi1pFwK76iyu+ekInGhuEfXA9Dcvpr/5iGTtAN9W1baNkXC9DEZjIiIbVfXa5tZRV0x/82L6m5dI1h/J2qFx9dtIb8MwDCMkzGAYhmEYIRHOBmNOcwuoJ6a/eTH9zUsk649k7dCI+sO2D8MwDMMIL8K5hmEYhmGEEWYwDMMwjJBodIMhIm+IyHERyQ4Ke1pEjopIlreN9MJvEpFMEdnqfaYFnTPAC98jIrNERLzwWqdYDxP9/ykih0WkqEr6LUXkXe+61olIz3DSLiJxIvKRiOz0pqh/rim0N5R+79gnIvKlp3+2iER74RGRd4LOXVwlrYjQLyKfi8iuoHM6eeGRkn9iRWSOiOz2ysGYxtbfQGW3bVDcLBE5ISIz66VdVRt1A64HrgGyg8KeBh6tJu7VQBdvvx9wNOjYemAIIMD/ASO88BnAb7393wLTw1T/YCAZKKpyzoPAbG//buDdcNIOxAE3ePuxwMqge99o2hv43id4nwK8D9wdSXnHC7sDeLtKWhGhH/gcuLaacyIl/zwDTPP2o4COja2/IfNOULxM4Pr6aG/0GoaqfgHkhxh3s7op0MGNDG/lWcJkXKFfo+4K53NuKvVap1ivDw2h3zu2Vr2ZfasQrH8hcKOIqz3Vl4bQrqpnVPUzL04ZsAno1tjaG0q/d6zAC4/BGb0KT4+IyDsi0gZ4BJhW5bSI0H8RIiL/APcD/+XFC6hqxSjwsC67wXFEpDduNvGVXlCdtDdnH8a/iVuB740aqtJjgM2qWgp0BY4EHTvihcElTrHegFyK/ovRFTgMoKp+4DSQ2LBSL6BO2sUtfHUbbrEsaB7tUAf9IrIENzlmIa6AQOTknd8D/41bjCyYSNEP8KbXLPJk0IMp7POPnFvs7fcisklE0kUkyQuLmLIL3IOrRVS8LNVJe3MZjNeAfwJSgVxcYahERPoC04FfVARVk0Zz+gNfqv6L0dTXViftIhIDvAPMUtV9FcHVpN/Y/0ud9KvqzbgmwZbABf0DTcgl6ReRVCBFVT9oYp01UZf7P0FVvwsM9bafVkSvJv1wyz8xuBr1KlW9BlgDvFARvZr0w67setyNK7+V0auJU7v2hmpzq6U9ridBbXEXO4b7c3YD1wWFJQM7g77fA7zu7e8CkoPi7Qo3/VXiV+3DWAIM8fZjcJOeSbhpB97AGYsm097Q996Lcx/waqTkHeAB4CvgAK5mXQZ8Hin6qzlnUtD9D/v8g3uwFgNR3vfuwLam0N+AZbc/sLtKWJ20N0sNw+uTqOB2INsLb4dbN+NxVV1VEUFddbtQRAZ71dl7OTeVeq1TrDc0l6q/FoL13wksV+9fbAzqol1EpgHfAqZWSa5JtXtaLkm/iLSpOMerJY0EdlajPyzzjqq+pqpdVLUn8ENcwR/mHQ57/SISIyIdvf0WwK0V5xAB+cfTkwEM84JuBLZ7+2Ffdj3u4fzaBdRVe0Na8xos4Tu46tNZ3BvSv+JW3tsKbPGEV7wl/TvOmmcFbZ28Y9d6N2gv8CrnRqkn4trUc7zPDmGqf4Z3fsD7fNoLbwWkA3twnmC9wkk77s1FgR1B4ZMbW3sD6k8CNnjxtwGvADGRlHeC0uvJ+W+VYa8fiMd551Tc/5eB6EjJP96xHsAX3jnLgMsjoewGpbUP+E6V9Ouk3aYGMQzDMELCRnobhmEYIWEGwzAMwwgJMxiGYRhGSJjBMAzDMELCDIZhGIYREmYwDMMwjJAwg2EYhmGExP8DpsXxtOqc+ywAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import gvar # library used to define gaussian random variables\n", "# now we can define our random variables, we are purposely biasing the inital guess \n", "# and giving a large error of +/- 10 km/s on both the velocity and the broadening\n", "\n", "velocity1_gvar = gvar.gvar(velocity1+3, 10) # velocity1 is known at +/- 10 km/s\n", "velocity2_gvar = gvar.gvar(velocity2-3, 10) # velocity2 is known at +/- 10 km/s\n", "broadening1_gvar = gvar.gvar(broadening1+3, 10) # broadening1 is known at +/- 10 km/s\n", "broadening2_gvar = gvar.gvar(broadening2-3, 10) # broadening2 is known at +/- 10 km/s\n", "\n", "fit = orb.fit.fit_lines_in_spectrum(spectrum, [halpha_cm1, halpha_cm1], step, order, nm_laser, theta, zpd_index, \n", " wavenumber=True, apodization=1, fmodel='sincgauss',\n", " pos_def=['1', '2'],\n", " pos_cov=[velocity1_gvar, velocity2_gvar], \n", " sigma_def='free',\n", " sigma_guess=[broadening1_gvar, broadening2_gvar],\n", " snr_guess=SNR)\n", "\n", "print('=== velocity ===')\n", "print('input velocity (km/s): ', velocity1_gvar, velocity2_gvar)\n", "print('fitted velocity (km/s): ', fit['velocity_gvar'])\n", "print('real velocity (km/s)', velocity1, velocity2)\n", "print('=== broadening ===')\n", "print('input broadening (km/s): ', broadening1_gvar, broadening2_gvar)\n", "print('fitted broadening (km/s): ', fit['broadening_gvar'])\n", "print('real broadening (km/s)', broadening1, broadening2)\n", "\n", "print('=== flux ===')\n", "print('flux (in the unit of the spectrum amplitude / unit of the axis fwhm): ', fit['flux_gvar'])\n", "# independant plot of the two lines models and the real lines\n", "pl.plot(spectrum_axis, spectrum, label='line1 + line2 + noise', ls=':', c='0.5')\n", "pl.plot(spectrum_axis, spectrum1, label='line 1', ls=':', c='red')\n", "pl.plot(spectrum_axis, spectrum2, label='line 2', ls=':', c='blue')\n", "models = fit['fitted_models']['Cm1LinesModel']\n", "pl.plot(spectrum_axis, fit['fitted_vector'], label='fit', ls='-', c='0.5')\n", "pl.plot(spectrum_axis, models[0], label='model 1', ls='-', c='red')\n", "pl.plot(spectrum_axis, models[1], label='model 2', ls='-', c='blue')\n", "pl.xlim((15200, 15270))\n", "pl.legend()\n", "pl.title('A much better fit')\n", "pl.savefig('gvar_good_fit.svg')" ] } ], "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.5" }, "toc": { "colors": { "hover_highlight": "#DAA520", "navigate_num": "#000000", "navigate_text": "#333333", "running_highlight": "#FF0000", "selected_highlight": "#FFD700", "sidebar_border": "#EEEEEE", "wrapper_background": "#FFFFFF" }, "moveMenuLeft": true, "nav_menu": { "height": "117px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": false, "toc_section_display": "block", "toc_window_display": false, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 4 }