{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction to MulensModel\n", "\n", "How to create and plot a model and then add some data and fit for the source and blend fluxes.\n", "\n", "This example shows OGLE-2003-BLG-235/MOA-2003-BLG-53, the first microlensing planets. See Bond et al. 2004 " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Import basic packages\n", "import MulensModel\n", "import matplotlib.pyplot as pl # MulensModel uses matplotlib for plotting." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Define a point lens model:\n", "\n", "my_pspl_model = MulensModel.Model(t_0=2452848.06, u_0=0.133, t_E=61.5)\n", "\n", "# Or a model with 2-bodies:\n", "\n", "my_1S2L_model = MulensModel.Model(\n", " t_0=2452848.06, u_0=0.133, t_E=61.5, rho=0.00096, q=0.0039, s=1.120, \n", " alpha=43.8)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8FNX6x/HPkxAIoTelBAQFQgsBE+lCqCKIgICAiCAo\n10ax/pR7RcCOqKCIghSDgiKIoiKCUr00TehBKSpCKJfQS4C08/tjNxgwkM0ms7Obfd6v176yM5mZ\n8yXAPjlnyhFjDEoppfxXgN0BlFJK2UsLgVJK+TktBEop5ee0ECillJ/TQqCUUn5OC4FSSvk5LQRK\nKeXntBAopZSf00KglFJ+roDdAVxRtmxZU7VqVbtjKKWUT4mLiztqjCmX3XY+UQiqVq1KbGys3TGU\nUsqniMhfrmynQ0NKKeXntBAopZSf00KglFJ+zifOESilVE6kpKSQkJDAhQsX7I7iEcHBwYSGhhIU\nFOTW/loIlFL5TkJCAsWKFaNq1aqIiN1xLGWM4dixYyQkJFCtWjW3jqFDQ0qpfOfChQuUKVMm3xcB\nABGhTJkyuer9aCFQSuVL/lAEMuT2z6qFQPm01NRUZsyYQVpamt1RlPJZWgiUT5s4cSKDBw9m+vTp\ndkdR6jKBgYE0aNCAevXq0atXL5KSkgB4+eWXqVu3LvXr16dBgwZs2LABgOjoaMLCwoiIiKB58+bs\n3Lnz0nqrb6jVQqB8WmJiIgDHjx+3OYlSlytcuDCbN29m+/btFCxYkA8++IB169bx7bffsnHjRrZu\n3cqPP/5I5cqVL+0ze/ZstmzZwoABA3j66ac9ltWyQiAiM0TkiIhsz7TuDRH5TUS2isiXIlLSqvaV\nfzDGAP41Hqx8z6233sqePXs4dOgQZcuWpVChQgCULVuWihUr/mP7li1bsmfPHo/ls7JH8BHQ8Yp1\nPwD1jDH1gV3Acxa2r/yAFgLlkuho+Ogjx/uUFMfyJ584lpOSHMtz5zqWT51yLC9Y4Fg+etSx/M03\njuXDh3PUdGpqKosXLyY8PJwOHTqwf/9+atasySOPPMKqVauy3Oebb74hPDw8R+3khmWFwBizGjh+\nxbqlxphU5+J6INSq9pV/0EKgvNX58+dp0KABUVFRVKlShcGDB1O0aFHi4uKYOnUq5cqVo3fv3nyU\nUaCAfv360aBBA9asWcP48eM9ltXOG8oGAXNtbF/lA+np6YAWApWNlSv/fh8UdPlySMjlyyVKXL5c\ntuzly+XLu9RkxjmCKwUGBhIdHU10dDTh4eHExMQwcOBAwHGOICoqyqXj5yVbThaLyL+BVGD2NbYZ\nIiKxIhKbcUJQqavRQqB8wc6dO9m9e/el5c2bN3PDDTfYmMjB4z0CERkI3AG0NRn9+iwYY6YCUwGi\noqKuup3ybzo0pHzJ2bNnGTp0KCdPnqRAgQJUr16dqVOnZrtf586dLz1HqGnTpsybNy9Pc3m0EIhI\nR+AZoJUxJsmTbav8SQuB8lZnz579x7rIyEjWrl2b5fYrMw8/ubA+L1l5+einwDogTEQSRGQwMAko\nBvwgIptF5AOr2lf+QQuBUrlnWY/AGNM3i9V6+6fKU1oIlMo9vbNY+bRrnGZSSrlIC4HyaRmFICBA\n/ykr5S7936N8mg4NKZV7WgiUT9NCoFTuaSFQPk0LgfJWgwYN4rrrrqNevXqX1q1fv57GjRvToEED\nateuzejRowHHHcX169cnPDycZs2asWXLlkv7FC1a1PKsWgiUT9OTxcpbDRw4kO+///6ydQMGDGDq\n1KmXHk999913A1CtWjVWrVrFtm3beP755xkyZIhHs+rk9cqn6cli5a1atmzJ3r17L1t35MgRKlSo\nADieOVSnTh0AmjVrdmmbJk2akJCQ4LGcoIVA+TgdGlLZGTFiRJYPf8uNBg0aMGHChBzv9/jjjxMW\nFkZ0dDQdO3ZkwIABBAcHX7bN9OnTuf322/Mqqkv01yjl07QQKF8yatQoYmNj6dChA3PmzKFjx8un\nbFmxYgXTp0/n9ddf92gu7REon6aFQGXHnd/crXTTTTfx8MMP8+CDD1KuXDmOHTtGmTJl2Lp1Kw88\n8ACLFy+mTJkyHs2kPQLl0/RksfIlixYtuvRvdvfu3QQGBlKyZEn27dvHXXfdxccff0zNmjU9nksL\ngbLEsmXLqFGjBufPn7e0He0RKG/Vt29fmjZtys6dOwkNDWX69Ol8/PHHhIWF0aBBA/r378/s2bMJ\nDAxk7NixHDt2jEceeeTSrGYZkpKSCA0NvfR666238jyrDg0pS4wYMYI9e/awe/du6tevb1k7WgiU\nt/r000//sW7w4MFZbjtt2jSmTZuW5fcyZuGzkvYIlCU89cGsl48qlXv6v0dZyuoxfO0RKJV7WgiU\nJTzdI1DqSv70byO3f1YtBMpS2iNQdggODubYsWN+UQyMMRw7duwfN6blhJ4sVpbwdI9AC4HKLDQ0\nlISEBBITE+2O4hHBwcGEhoa6vb8WAuXTtBCorAQFBVGtWjW7Y/gMHRpSltKhIaW8nxYCZQkdGlLK\nd2ghUJbSHoFS3k8LgbKE9giU8h1aCFS+oIVAKfdpIVCW0qEhpbyfFgJlCU99MGc8kEsLgVLu00Kg\nLKU9AqW8nxYCZQk9WayU77CsEIjIDBE5IiLbM60rLSI/iMhu59dSVrWvvIM/POtFKV9nZY/gI6Dj\nFeueBZYZY2oAy5zLKh/ydI9AC45S7rOsEBhjVgPHr1jdFYhxvo8BulnVvvIPWgiUyj1PnyO43hhz\nyPn+MHC9h9tXHuapk8VKKffZdrLYOP4HX/V/sYgMEZFYEYn1l0fJ5ieevnxUC4JS7vN0IfifiFQA\ncH49crUNjTFTjTFRxpiocuXKeSygylv6Aa2U9/N0IfgaGOB8PwBY6OH2lYfoyWKlfIeVl49+CqwD\nwkQkQUQGA68B7UVkN9DOuayU27QQKJV7ls1QZozpe5VvtbWqTeV99GSxUt5P7yxWltChIaV8hxYC\nZSntESjl/bQQKEt4qkeQmpoKaEFQKje0EChLWf0BrYVAqdzTQqAs4ekegVLKfVoIlE/THoFSuaeF\nQFnKU0NDSin3aSFQltCTxUr5Di0EylJ6slgp76eFQFnCUz2ClJQUj7SjVH6mhUD5NO0RKJV7WgiU\npfRksVLeTwuBsoSeLFbKd2ghUJbSk8VKeT8tBMoSGT0CHRpSyvtpIVCW0h6BUt5PC4GyhPYIlPId\nWgiUpbRHoJT300KgLJHRI0hPT7e0neTkZEALgVK5oYVAWcJTQ0NaAJTKPZcmrxeRmsDTwA2Z9zHG\ntLEol8onPPVBrQVBKfe5VAiAecAHwIdAmnVxVH7hqR6BUir3XC0EqcaY9y1NovIlKwtB5iuGtOAo\n5T5XzxF8IyKPiEgFESmd8bI0mcoXrPyAPn/+vEfaUSq/c7VHMMD59elM6wxwY97GUfmFJ4aGLly4\nYNmxlfInLhUCY0w1q4Oo/El7BEp5P1evGgoCHgZaOletBKYYY3RWEJUl7REo5TtcHRp6HwgCJjuX\n+zvXPWBFKJV/WFkIkpKSPNKOUvmdq4XgFmNMRKbl5SKyxYpAKn/wRI/gzJkzl95rIVDKfa5eNZQm\nIjdlLIjIjeTifgIReVxE4kVku4h8KiLB7h5LeTdPFQKllPtc7RE8DawQkT8AwXGH8f3uNCgilYBh\nQB1jzHkR+RzoA3zkzvGUd9MegVLez9WrhpaJSA0gzLlqpzHmYi7bLSwiKUAIcDAXx1JeyNNDQ0op\n912zEIhIG2PMchG564pvVRcRjDELctqgMeaAiIwH9gHngaXGmKVZtD0EGAJQpUqVnDajvIT2CJTy\nftmdI2jl/Noli9cd7jQoIqWArkA1oCJQRETuvXI7Y8xUY0yUMSaqXLly7jSlbOSJHsHp06cvvddC\noJT7rtkjMMa84Hw71hjzZ+bviYi7N5m1A/40xiQ6j7MAaAZ84ubxlBez8gP6+PHjlh1bKX/i6lVD\nX2Sxbr6bbe4DmohIiDh+bWwL/OrmsZSXs7oQlCxZ0vJ2lMrvsjtHUAuoC5S44jxBccCtSz6NMRtE\nZD6wEUgFNgFT3TmW8l6eGBo6fvw4pUuX5uTJk5a1oZQ/yO6qoTAc5wJK4jgvkOEM8KC7jTqHnF7I\ndkPl8zxRCP744w/tESiVC9mdI1gILBSRpsaYdR7KpPIBT/QIDh8+TL169SxvR6n8ztUbyjaJyKM4\nhokuDQkZYwZZkkrlG1Z9QBtjOHjwIB06dLDk+CQnw6ZNUKcOFCsG6ekg4ngplc+4erL4Y6A8cBuw\nCgjFMTykVJas7hGcOHGC5ORkKlasaE07ixdDkyawfbtjeckSqF4d4uPzth2lvICrhaC6MeZ54Jwx\nJgboDDS2LpbKL6wqBIcOHQK4VAhyLTUVxoyBuXMdy9HRMH8+F2+8kd27d7Pz5Eku1K0LNzkfuZWi\nT2BX+YerhSDjX/1JEakHlACusyaSyk+sKgQHDzqeSlKhQoW8aSctzdELWLcOYwyL/vtfOn74IcWr\nVKFmzZrUuuceiixaRJtOnZgbE0N6kyYwaVJu/xhKeQVXzxFMdd4R/B/ga6AoMMqyVMrnWT00dGWP\nwO12jh2D4sWhUCFYvpzDp09zf6dOfP/994SGhjJ06FDq16+PiBAfH88XX3xBn4EDeaN0aWYFBVEn\nr/5AStnI1YfOTXO+XY3OU6xywKpC8OeffyIiuXsO1dmz0KoVREZCTAwbf/uNrl27cvz4cSZMmMCj\njz5KgQKX/xd5+eWX+eyzzxgxYgQ3Dx/OjGLFuKdxY7jxRj2RrHyWS0NDIvKKiJTMtFxKRF6yLpby\ndVb3CPbs2UPlypUpXLiw++0ULQqPPAIDB7Jx40batGlDQEAAa9asYfjw4f8oAgCBgYH069ePbdu2\n0aRJE/r168ebYWEwbVoWDSjlG1w9R3C7MebS7ZvGmBNAJ2siqfzEykJQvXp193Y2Bg4fdrx/5BF2\nh4bSvn17SpQowU8//USDBg2yPUT58uVZsmQJd/fqxVNpaUzS5x4pH+ZqIQgUkUIZCyJSGCh0je2V\nAqwvBG71PN55B+rXh/37OXPmDN26dUNEWLZsWY6GmgoVKsTsOXPo2rUrQ599ltmffAJaEJQPcvVk\n8WxgmYjMdC7fD8RYE0nlB1YODR09epSjR49Ss2bNS+ty1E779pCQgKlUiQfvuYfffvuNH374wa0e\nRoECBfjss8+4/fbbGTxgADWrV+eWrVsdJ5+V8hEu9QiMMa8DLwO1na8XjTHjrAym8gcrCsGWLVsA\niIiIuFRwXJKe7vhapw688QZzP/+cuXPn8uKLL9KmTRu38wQHBzNv3jzKly3LXUeOcEQfgqd8jKtD\nQxhjFhtjnnK+llgZSvk+K3sEWRUCl9oZPBiefRZwXH766KOP0rhxY5555plcZypbtixffv89Ry9c\noE/fvqRnFB2lfMA1C4GI/Nf59YyInM70OiMip6+1r/JvGR/QVnwgbty4kYoVK5KjmevS0yEkBIKC\nABg+fDhJSUnExMRkeXWQOxo2bMh7773HihUreOumm+C0/hdRviG7HsF9AMaYYsaY4plexYwxxT2Q\nT/moIOcHbooFj2JYs2YNzZo1A3LQ8wgIgPfeg7FjWbFiBfPmzWPkyJGEhYXlabb777+f7q1aMfKv\nv9iycmWeHlspq2RXCOYBiMgyD2RR+UhGIUhOTs7T4yYkJLB3715atGhx2fprFoJ334WdOwFITUtj\n+PDhVK1alaeeeipPs4GjME2dP58y119Pv5EjuXDhQp63oVRey64QBIjISKCmiDxx5csTAZVvKliw\nIJD3hWDNmjUAlwpBtieLjx2D0aPhww8BmDp1Ktu2bWP8+PGXbkbLa2XLlmXmzJnEx8fzYteucPGi\nJe0olVeyKwR9gDQcl5kWy+KlVJas6hEsX76cokWLEhERAbgwNFSmjOPR0WPGcO7cOcaMGUN0dDR3\n3XVX1tvnkY4dOzLwttsYt3QpWydOtLQtpXIruxnKdgKvi8hWY8xiD2VS+UBAgON3jLw8R2CM4dtv\nv+W2225z7QTv3r1QtSqULw/Ae+PGceTIEb788sucXXbqpvGzZ7OoRg0emD+fdU8+SWBgoOVtKuWO\n7K4autf5to4ODSl35GWPYOPGjRw8eJAuXf6ePvuqPYKDB6FePXjlFQDOnDnDuHHj6Nix46UTzVYr\nU6YM70yezC+//MK7r77qeLSFUl4ou6GhIs6vRdGhIZUDGR/MeVkIvvnmG0SETp3++ZirfxSCcuUc\nE8306QPAxIkTOXbsGGPHjs2zPK7o3bs3nW+9lX8//zx/TZ7s0baVclV2Q0NTnF/HeCaOym/yamjI\nGMMXX3xBs2bNLrt/4KpDPEFB8OSTAJw6dYo333yTO++8k1tuuSVP8rhKRJgcE0PtmjUZ8dVXfPno\nox5tXylXuPoY6nIiMlJEporIjIyX1eGU78rrHsGmTZvYvn07/fr1u2z9P4aGUlKgRw9YvfrSNlOm\nTOHkyZOMGmXPXEpVqlVj1Esv8dWPP7J4sZ5qU97H1UdMLMQxPeWPwKJML6WuKa8KQUxMDAULFqR3\n797X3nDfPtiyBZzP+0lOTmbixIm0bduWyMjIPMnijscff5ywGjUY2qcPFzZtsi2HUllx9d76EGPM\n/1maROUrGb+h58XQ0IULF5gzZw533nknpUuXvux7/+gR3HQT7Nhx6VESc+bM4eDBg8yYYW8HtmDB\ngkx65RXa9+rF+NGj+c/ChbbmUSozV3sE34qITkSjciwvegRz5szh6NGjPPzww1fdxhgDS5ZAaioU\nLAgiGGMYP3489evXp0OHDrnOkVvtevakV/fuvLx0KXv37rU7jlKXuFoIhuMoBuf1oXPKFRm/oV/M\n5V21xhgmTJhA/fr1ad269T++f+lk8eHD0LEjTJp06XuLFy8mPj6ep556yiP3DbjirXfeITAwkBGD\nBsGZM3bHUQpwfT6CYsaYAGNMYX3onHJFRiFISkrK1XGWL1/Otm3bGDFixDU/zM1118GXX8KDD15a\nN2HCBCpVqkQf5yWk3iA0NJRRw4ezcMUKvsuUVSk7uXrV0M1ZvG4SEbee3ysiJUVkvoj8JiK/ikhT\nd46jvN/Zs2fd3tcYw4svvkj58uXp27fvtTcWgW7doIjj1pddu3bxww8/8NBDD1163IW3GPHCC9Sq\nWJFhP/+sD6VTXsHVoaHJwHrgQ+drPY4nk+4UEXcGXycC3xtjagERwK9uHEN5sYwewZlcDH8sX76c\nVatWMXLkSIKDg7Pe6NgxBDC7dl22evLkyQQFBfGgF/7WXbBgQd6NieH3P/9k3Ouv6x3HynauFoKD\nQENjTKQxJhJoAPwBtAdyNGWliJQAWgLTAYwxycYYndsvn3K3EBhjeP755wkNDb32h7lz8hdT/O+R\nynPnzvHRRx/Rs2dPrr/+erfat1q7du24u1s3Xh07lj/0oXTKZq4WgprGmPiMBWPMDqCWMeYPN9qs\nBiQCM0Vkk4hME5Ei2e2kfEtGj8DdoaG5c+eybt06Xnjhhav3BgCqVUMCAuC66y6tmj17NqdOneJR\nL7+L962JEykgwrCPPrJkSk+lXOVqIYgXkfdFpJXzNRnYISKFgJxeKF4AuBl43xjTEDgHPHvlRiIy\nRERiRSQ2MTExh00ou+VmaOjcuXM89dRT3Hzzzdx///1Xa8Ax4czp04jzUtGMdt977z0iIiI89nA5\nd1WqUoXRr73Goi1b+Prrr+2Oo/yYq4VgILAHGOF8/eFclwL885q+a0sAEowxG5zL83EUhssYY6Ya\nY6KMMVE5mptWeZWUlJQcF4NXXnmFAwcO8I7zUsss/fILDB8O8+dftnrt2rVs3bqVRx991GsuGb2W\nYcOHU7duXYYPGUJSXJzdcZSfcvXy0fPGmDeNMd2dr/HGmCRjTLoxJkd9f2PMYWC/iGRMFtsW2JHD\n3MrLZR7qOHDggMv7bd68mXHjxtG/f3+aN29+9Q0bNYJNm2DAgMt6BNOmTaNo0aLcc889bmf3pKCg\nICa/8QZ/HTnCy4MG2R1H+SlXLx+t4bzcc4eI/JHxykW7Q4HZIrIVx4nnV3JxLOXlDh486NJ2ycnJ\nDBgwgLJly/L2229ffcNTpxxfIyLA2WMwxnDmzBk+//xz+vTpQ5EivnPaqeXtt9P/jjt449df2emc\nW1kpT3J1aGgm8D6QimMoaBbwibuNGmM2O4d96htjuhljTrh7LOWdMvcIXC0EY8aMYevWrUyZMoUy\nZcpkvdFvv0Fo6GVDQhlDQHPnziUpKYnBgwe7H9wmb0ybRkhICI898gjm3Dm74yg/42ohKGyMWQaI\nMeYvY8xooLN1sVR+cJ3zSh5XCsHixYt55ZVXGDRoEHfeeefVNyxdGu65B2699dKqjKGh6dOnU6dO\nHRo3bpzr7J52/fXX89KoUfy4fDnzsnvCqlJ5zNVCcFFEAoDdIvKYiHTHMWuZUlkyxlCsWDGKFStG\nQkLCNbf966+/uPfee6lfvz6TMj0rKEvXXQdTpsAV9wfs2LGD9evXM2jQIJ84SZyVh4cPp2GFCjy+\nbl2ubsRTKqdy8tC5EGAYEAn0BwZYFUr5PmMMIkJYWBjbtm276nYnTpygc+fOpKamMn/+fAoXLny1\nDWHAAMd8A1cQERYtWkSBAgXo379/Xv0RPC4wMJDJCxZw8PhxRo8ebXcc5UdcvWroF2PMWWNMgjHm\nfmPMXcaY9VaHU76vSZMm/PLLL6Smpv7jexcuXKBbt27s2rWLr776iho1alz9QLGxsHChoyBkIS0t\njS5dulwajvJVTZo04cEHHmDC228T9/rrdsdRfuKahUBEvr7Wy1Mhle/J6BE0bdqUc+fOEXfFNfJn\nzpyhU6dOrF69mpiYmCwfMX2Z9u1h/37HlUJXyBgK8sWTxFkZ9/rrXF+gAINffz3P5nxW6lqy6xE0\nBUKBn4DxwJtXvJS6KhGhU6dOhISEMHny5Evrd+/eTatWrVi9ejUff/zxtZ8smpICa9Y43hcrdtV2\nKlSowG233ZaX8W1TsnRp3ps+nS0nTjB+/Hi74yg/kN1jpMvjeLBcX+AeHPMUf5r5uUNKZSXj8tGS\nJUvy0EMP8dZbb1GqVCkuXrxITEwMhQsX5uuvv6ZTp2wmvpsyBYYOddxJHBWV5SY33ngj/fv3p0AB\nt56K7pW69+9Pj4ULGTNmDHc1akRY27Z2R1L5mTHGpRdQCMdjJRKBx1zdLy9ekZGRRvmWu+++24SF\nhRljjLlw4YLp06ePAUzhwoVN//79zb59+1w70PnzxnzyyTU3SUtLM+np6bmN7HUOHThgSgYEmBYl\nSpi0tDS74ygfBMQaFz5jsz1ZLCKFROQuHDeQPQq8A3xpVWFS+U+hQoX49NNPOXv2LKdOnWLWrFlU\nrlz52julpTmGhYKDoV+/a24aEBDgs5eMXkv5ihV567HH+O+pU0yZMsXuOCofy+5k8SxgHY6Hwo0x\nxtxijHnRGOP6w2OUXzLOk8WZFSlSxPXZwt56y/E8oZP+PVXFwAkTaNeuHf/3f//HX7//bncclU9l\n1yO4F6iB4z6Ctc6J63XyeuWSXP2WXrMmNG4MJUrkXSAfJCJ8+OGHkJzMgMhI0vUqImWBaxYC45iw\nvpjzVTzTSyevV9dkcjvRSteu8MEHjrmI/VzVqlWZeN99rDp1irff1Iv1VN5z9c5ipXLMrR7B+PEw\nc6bO43uFgVOm0K1bN0a+8ALbt2+3O47KZ7QQKEu41SNIT4elS2HZMu0JXEFEmDp1KiWLFePeVq24\neFpHZlXe0UKgLJHVyeJsBQTA4sUwdao1oXxcuXLlmPbQQ2w5fpxRw4bZHUflI1oIlHeYMQNOn3ZM\nNBMSYncar9XlpZcYMmAA42JiWLJkid1xVD6hhUBZIkc9gl27YMgQeP99a0PlE29Pnky9evW49+67\nObB5s91xVD6ghUBZxuVCULMmrF0LTzxhbaB8IiQkhM/fe4+k06e5p0uXLJ/sqlROaCFQlnDpZHFS\nEmQ8lbRRI3D1ZjNF7ZYteX/UKFYnJDBmzBi74ygfp4VAWcKloaHRo6F5c8hmBjOVtfvGjOH+++/n\n5Zdf5vvsZnZT6hq0ECj7PPccfPSRYzJ65ZZJkyZRr3hx+g4bxp4tW+yOo3yUFgJliWv2CPbudTxU\nrlQp6NPHo7nym5CQEBZ+9RUBRYvS9Z57dK5j5RYtBMqzTp6Epk1Br4PPM9Wio/n8yy/ZuXMn/bt3\nJz0tze5IysdoIVCWuGqPoGRJGDUKHn7Y86HysbZt2/Lmk0+ycNkyxnbpYncc5WO0ECjLXFYI0tPh\n8GHH+4cfhnr17AmVjw179VUG3nwzYxYvZs6cOXbHUT5EC4GyxD8uH33tNcfE8wd0KgurSEAAH6xd\nS3R0NAMHDmTFp5/aHUn5iPwzyavyKv8YGurRw3HfQMWK9oXyA4UKFeLLL7+kRfXqdOvXjzWVKlGv\nZUu7Yykvpz0CZa29ex1fw8LgpZf0qaIeULJkSb6bN4+ixYtze79+HNBemMqGbYVARAJFZJOIfGtX\nBmUdYwxy/jzUquWYX0B5VJXWrflu1SpOnTpFh9atOarTXKprsLNHMBz41cb2lcUkJASeecYx25jy\nuIiICL7+4gv+2LOHDjffzIkTJ+yOpLyULYVAREKBzsA0O9pXFtu+HZOS4hgGGjsWSpe2O5Hfim7f\nni+ff57tSUncfvvtesOZypJdPYIJwDNAuk3tK6skJUH79pitW3M3eb3KMx3HjOHzefOIjY3ljqZN\nOXfsmN2RlJfxeCEQkTuAI8aYuGy2GyIisSISm5iY6KF0KtdCQmD6dMfJYeU1unXrxuw33uC/8fHc\nHhnJaZ3qUmViR4+gOXCniOwFPgPaiMgnV25kjJlqjIkyxkSVK1fO0xlVTsXFwQ8/AJDesSOnUlK0\nR+Blej/+OHOee451Bw7Qvn17jh8/bnck5SU8XgiMMc8ZY0KNMVWBPsByY8y9ns6h8pAx8PjjMGIE\np44fp3v37qxbt4527drZnUxdofcrr7BgwQI2b95M61q1OLJ7t92RlBfQ+whU7onA558T/+673NKk\nCd999x2VKXbuAAASw0lEQVQTJ07kpZdesjuZykKXLl1Y9NJL7ElM5NbWrfnzzz/tjqRsZmshMMas\nNMbcYWcGlQszZsBDD4ExLFi7liZdu3L69GmWL1/OsGHDdGjIi7V7+mmWfvUViUlJNG3alNh16+yO\npGykPQLlvr17Sf/9d54fOZIePXpQp04d4uLiuPXWW+1OplzQvGtX1qxZQ+GAAFo1b863EybYHUnZ\nRAuBypmkJPjjDwBOPv44dxYsyEuvvcagQYNYtWoVlSpVsjmgyonatWuzbsECahctStcnn+Tdd991\nbb5pla9oIVA507s33HYbv27ZQqPGjVmydCmTJk1i2rRpBAcH251OuaF8kyasOnSIzp07M2zYMAa3\nasWFs2ftjqU8SAuBypn//IfZXbvSqEULTp06xbJly3j00Uf1fICPK1KkCF999RXP9+vHzJ9+olWD\nBiQkJNgdS3mIFgJ1bcbAuHEwcSLnzp1j8NSp3Pvmm0RERBAbG0tLfcRxvhEQEMDYTz7hyzFj2PG/\n/xEZGckq570hKn/TQqCyt2ED8UuX0qhRI2bOnMm///1vVq5cSeXKle1OpizQbdQofv75Z0oWK0ab\nDh0Y0707aToPcr6mhUBlbfVqOHyYdGN4p1kzblmxgqNHj7JkyRJeeuklChTQOY3ys9q1axO7fDn9\nbriB0V99RZs2bXSoKB/TQqD+6dgx6NSJP594grZt2zL8qaeIjo5m8+bNtG/f3u50ykOKVanCrL17\niYmJIS4ujoiwML7497/tjqUsoIVA/e3IEQBM6dJMefBB6n/zDXFxcUybNo1FixZRoUIFmwMqO9x3\n331sXL+easbQ85VX6NOnD0ePHrU7lspDWgiUw7JlUKUKO2bMoHXr1jw0YQKNGzdm27ZtDB48WK8K\n8nM169Vj3ZEjvPj88yxYsIA6tWrxxZgxdsdSeUQLgb+7eBGAcxERPBceTsS//sXWrVuZMmUKS5cu\n5YYbbrA5oPIWQUWL8p+xY4mLi6OyMfQcPZpe3brpuYN8QAuBPxszBtOiBQsXLKBuVBSvxcZy7733\nsnPnToYMGUJAgP7zUP8UHh7O+j//5OV//YtvlyyhVq1ajB8yhBTnLxXK9+j/dH+TmgrpjonhNhYq\nRNvDh+nWowdFixZl9erVzJw5E53/QWUnqHhxRn7wAfHx8bSOiODpDz+k4U03sWrVKrujKTdoIfAn\nhw9Dw4bsf/tt7rvvPiKfe46t58/z7rvvsmnTJn1YnMqxG2+8kW9++omvn3yScwUKEB0dTa/27dn9\n0092R1M5oIXAHzgnLD8CPJOaSs2RI/n888959tln+f3333nssccICgqyN6PyXQEBdBk/nvgdOxgz\nejSLly+nTsuWDH3sMY44r0RT3k0LQX736qsk1qjB/z3xBNVuuok3d+2iV+/e7Ny5k1dffZUSJUrY\nnVDlEyEhIYx64QX2rF/PA9278/4HH1C9enVe7NKF04cO2R1PXYMWgvzor7/g9GkSExN5Nj6easeP\n88aECXTv3p0dO3Ywa9YsvRpIWab8Lbfw/oIFbN++nTb16jHq22+pWrMmL774IqdOnbI7nsqCFoL8\n5uBB9lSvziNt23LDDTcwbs4c7uzZkx07dvDJJ58QFhZmd0LlJ2rVqsVXa9fyy6xZ3NqmDaNGjeKG\nihUZ3aEDJ48ftzueykQLQX6wfj3MmMG6devoMXQoNdPSmL51K3379iU+Pp45c+ZQq1Ytu1MqPxXV\nvz8LFy5k48aNtC5dmjE//EDlG25gxIgR/Pnbb3bHU4D4wmxEUVFRJjY21u4Y3iU9HQICuHDhAvM7\ndGDyhg2sS06mVKlSPPzwwzz22GP6SAjllbasX8+bkyfz6aefkp6aSo/ISJ6YNIkmTZrYHS3fEZE4\nY0xUdttpj8AXrVzJnkqVeHrIEEJDQ+n/008crVyZiRMnsm/fPl5++WUtAsprRTRpwqxZs9i7cSPP\nNGzID7t20bRpU5pGRhLz8MOcP3nS7oh+R3sEviA9HVas4GLJkny7dy9T332XpatWERgYSLdu3Xjo\noYdo06aN3gmsfNLZs2eZOXMm7734IjsTEylZvDj3DRzIvx54gDrh4XbH82mu9gi0EHizixcxBQuy\nftUqZt12G3OBE8nJVKpUiSFDhvDAAw9QsWJFu1MqlSdMWhqrpk1jysqVfPHFF6SkpNCiTBkGvfEG\nPXr0oHjx4nZH9DlaCHzc73ffzZyff2ZWUBB79uyhcKFC3NW9O/fdfz9t27YlMDDQ7ohKWSYxMZGP\n7r+fD9evZ/exYwQHB9O1Rg3uHTCA24YN0xsgXaSFwNf89BO733yTeVFRzF+wgE2bNgHQunVr7rvv\nPnr06EGxYsVsDqmUZxlj2LBhA59Mm8ZnM2ZwzBjKli1L71696FW3Li2GDCFQi8JVaSHwdmfPYhYt\n4rcqVZj/44/MnzaNrfv2AdCkSRN69uxJr169qFKlis1BlfIOySdPsmTpUj6eP5+vv/qKiykpXFei\nBN1696bHnXfSul07ggoVsjumV9FC4I327uVCWhqr9uxh0Ucfseizz/jD+a3mzZrRs2dPevTsqZPC\nK5WNMwcPsvjVV/ni0CEWff89586do5QId959N11796Zt27Z6TgEtBN4hJQVOnCAhOZnvFixg0YgR\n/FigAEkpKQQHB9OmYUM633MPXbt3p1KlSnanVconnT9/nh/efJMvPvmEr//3P06ePEmBgABuLV2a\nTs88w+2dOlGnTh2/nGXPawuBiFQGZgHXAwaYaoyZeK19fKYQGAPHjnE8IIAVK1awfMgQliUns/Ps\nWQCqlC1L59tuo3PfvrRu3ZqQkBCbAyuVv6SkpLB27VoWjxnDd1u2sM35KIsqRYtye+3atHvmGaKj\noylbtqzNST3DmwtBBaCCMWajiBQD4oBuxpgdV9vHawuBMXDkCKcLF2bNmjUsf+IJlu3dy+aLFzHG\nUKRQIVrWq0ebPn3o2LEjdevW9cvfSpSyy/79+1m8eDGLx47lx8REziYnA1C/WDFaN2tGm0ceoWXL\nlpQsWdLmpNbw2kLwjwAiC4FJxpgfrraN1xSCtDRMfDz7ihdnzdq1rHn7bdbExbFNhPT0dAoWKEDT\nG2+kbb9+tGnblkaNGullbkp5iZSUFGJjY1nx/fcsf+cd1pw7x4WUFAICAmgYHEzLDh1o1q8fTZs2\nzTdDtT5RCESkKrAaqGeMOX217WwrBGfPkrR8OVtCQvhlxw7WzJrFmrg4Dji/XTQkhKZVqtCse3da\ntGlDs2bNdLhHKR9x8eJFNmzYwPL581kxaxY/nz/PBWePoUpgIM1at6bZnXfSNDKSiIgIgooUsTlx\nznl9IRCRosAq4GVjzIIsvj8EGAJQpUqVyL/++svyTBcOHGDLuHHEFi9O3IEDxK5axY4//iDN+f0q\nlSrRvHJlmvfoQfN27QgPD9cbu5TKJ5KTk9myZQtr581j7Zw5rDOG/QcPAlAYaNiwIZEtWnBz1apE\nVqtG7c6dKVCwoL2hs+HVhUBEgoBvgSXGmLey2z5PewTGQGoqiSdPsm39eraNHcvWkiWJS0wkPj6e\n1NRUAMqVK0dUw4ZElitH1B13EHXrrfmmu6iUcs3+/ftZN2cO6+bNIzY4mE2bN3Pu3DkAgoODiYiI\nIPL664ksVYqGQ4dSu25dgoODbU79N68tBOI4WxoDHDfGjHBln9wUgvPr1rFj7162XbzIti1b2Dp5\nMtuCgvif8y8ToEyRIkQ2b05UVBRRNWoQ2aYNlStX1hO7SqnLpKWlsXvlSuK++Ya4gADi4uLYtG4d\nZ1JSAAgICKBGqVKEFytGvYEDCQ8Pp16VKtwUEWHLHdDeXAhaAD8B24B05+qRxpjvrraPu4XgiSee\nYOLbb19qJDg4mLqlShFepw7hnToRHh5OeHg4119/vX7oK6Xckp6ayp4NG9iUkEB8fDzb5s5l+6FD\n/H72LBmfr8EBAdRp0IC6desSdv48YbVqUbNXL2rUqEHhwoUty+a1hcAd7haCBQsWsOXbbwmPjCS8\nXTuqV6+uY/pKKY9ISkpix44dbJ8+ne3797M9NZX4+HgSEhIu265KcDBhlSpRs2NHwsLCqJmcTFiL\nFlSOisr155UWAqWU8kJnjx5l99at7EpMZOdvv7ErJoadaWnsPHGCM2fOXNouKCiIatWqMWXKFKKj\no91qy9VCUMCtoyullHJL0bJladimDQ0zVrzwAuB40ur/Dh5k14IF7Dp1it/PnWPPnj0euQtaC4FS\nSnkBEaF8pUqUHzqUlh5uW+c2VEopP6eFQCml/JwWAqWU8nNaCJRSys9pIVBKKT+nhUAppfycFgKl\nlPJzWgiUUsrP+cQjJkQkEXB3QoKywNE8jJNXNFfOaK6c0Vw5k19z3WCMKZfdRj5RCHJDRGJdedaG\np2munNFcOaO5csbfc+nQkFJK+TktBEop5ef8oRBMtTvAVWiunNFcOaO5csavc+X7cwRKKaWuzR96\nBEoppa7B5wqBiFQWkRUiskNE4kVkuHN9AxFZLyKbRSRWRBo515dxbn9WRCZdcaxIEdkmIntE5B3J\nxcTFbuRqLyJxzvbjRKSNl+Rq5Fy3WUS2iEh3b8iVab8qzr/Lp7whl4hUFZHzmX5mH3hDLuf36ovI\nOuf220Qk2O5cItIv089qs4iki0gDL8gVJCIxzvZ/FZHnMh3LzlwFRWSms/0tIhJtRS6MMT71AioA\nNzvfFwN2AXWApcDtzvWdgJXO90WAFsBDwKQrjvUz0AQQYHHG/h7K1RCo6HxfDzjgJblCgAKZ9j2S\nadm2XJn2mw/MA57ykp9XVWD7VY5lZ64CwFYgwrlcBgi0O9cV+4YDv3vJz+se4LNM/wf2AlW9INej\nwEzn++uAOCAgr3P5XI/AGHPIGLPR+f4M8CtQCTBAcedmJYCDzm3OGWP+C1zIfBwRqQAUN8asN46f\n6iygmwdzbTLGHHSujwcKi0ghL8iVZIxJda4Pdm5n+8/LmaEb8CeOn1fGOttzZcULcnUAthpjtjj3\nOWaMSfOCXJn1BT4Dr/h5GaCIiBQACgPJwGkvyFUHWO7c/ghwEojK61xuVQ9veeH4bWyf8wdY2/l+\nP3AAxx11mbcdSKYeARAF/Jhp+VbgW0/ncm7fMyOLN+QCGuP4sD0LdPeGXEBRYJ3z62icPQIvyFUV\nOAdsBlYBt3pJrhHAx8ASYCPwjDfkumL734F63pALCMJRlBKdf59DvCTXEBw94AJANRyFoEde5/K5\nHkEGESkKfAGMMMacBh4GHjfGVAYeB6b7Qi4RqQu8DvzLW3IZYzYYY+oCtwDPZYwt25xrNPC2Meas\nVVnczHUIqGKMaQA8AcwRkeJZHdPDuQrgGBLt5/zaXUTaekGujO0bA0nGmO1WZcphrkZAGlARxwfu\nkyJyoxfkmgEkALHABGCtM2feyovK5ukXjuq9BHgi07pT/H05rACnr9hnIJf3CCoAv2Va7gtM8WQu\nIBTHGGFzb8p1xb7Lcfz2YWsu4Ccc47Z7cfxWdBx4zO5cWey70kt+Xn2AmEzbPQ88bXeuTN9/Gxjp\nLf/ugfeA/pm2mwHcbXeuLPZdi2O4KE9z+VyPwHlmfDrwqzHmrUzfOgi0cr5vA+y+1nGMMYdwjAE2\ncR7zPmChp3KJSElgEfCsMWaNF+Wq5hwnRURuAGoBe+3OZYy51RhT1RhTFcdvRq8YYybZnUtEyolI\noPP9jUAN4A+7c+H4oAkXkRDn32crYIcX5EJEAnB8yH6Wsc4Lcu1zLiMiRXCchP3N7lzOv78izvft\ngVRjTJ7/Pbpd1ex64ejmGhxXRGx2vjo518cBW4ANQGSmffbi+A3yLI5uVh3n+ihgO46xykk4K7In\ncgH/4e+x5YzXdV6Qqz+O8wObcYwtd8t0LNtyXbHvaC6/asjOn1ePK35eXbwhl3Ofe53ZtgPjvChX\nNLA+i2PZ+fdYFMdYfDywA3jaS3JVBXbiOKn8I5efy8uzXHpnsVJK+TmfGxpSSimVt7QQKKWUn9NC\noJRSfk4LgVJK+TktBEop5ecK2B1AKW8iImWAZc7F8jju4kx0LicZY5rZEkwpC+nlo0pdhYiMBs4a\nY8bbnUUpK+nQkFIuEpGzzq/RIrJKRBaKyB8i8po4nrP/s/P58Dc5tysnIl+IyC/OV3N7/wRKZU0L\ngVLuicAxx0VtHHdj1zTGNAKmAUOd20zE8aC8W3DcgTzNjqBKZUfPESjlnl+M43kviMjvOCYWAdgG\ntHa+bwfUyTRxVHERKWo89ARVpVylhUAp91zM9D4903I6f/+/CgCaGGMumxRJKW+jQ0NKWWcpfw8T\nIc65eZXyNloIlLLOMBzTCm4VkR04ziko5XX08lGllPJz2iNQSik/p4VAKaX8nBYCpZTyc1oIlFLK\nz2khUEopP6eFQCml/JwWAqWU8nNaCJRSys/9Pw78XGtU8U8pAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot those models:\n", "my_pspl_model.plot_magnification(\n", " t_range=[2452810, 2452890], subtract_2450000=True, color='red', \n", " linestyle=':', label='PSPL')\n", "my_1S2L_model.plot_magnification(\n", " t_range=[2452810, 2452890], subtract_2450000=True, color='black', \n", " label='1S2L')\n", "pl.legend(loc='best')\n", "pl.show()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/jyee/MulensModel/source/MulensModel/utils.py:52: RuntimeWarning: invalid value encountered in log10\n", " mag = MAG_ZEROPOINT - 2.5 * np.log10(flux)\n" ] } ], "source": [ "# Suppose you also had some data you want to import:\n", "\n", "OGLE_data = MulensModel.MulensData(\n", " file_name='../data/OB03235/OB03235_OGLE.tbl.txt', comments=['\\\\','|'])\n", "\n", "MOA_data = MulensModel.MulensData(\n", " file_name='../data/OB03235/OB03235_MOA.tbl.txt', phot_fmt='flux', \n", " comments=['\\\\','|'])" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/jyee/MulensModel/source/MulensModel/utils.py:52: RuntimeWarning: invalid value encountered in log10\n", " mag = MAG_ZEROPOINT - 2.5 * np.log10(flux)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXecFOX5wL/P7lXKCVKVAw4UUKTJIYKKBdRQrFixBNT8\n7F00URNLNGgQTUxMLFFRQhErFqxoFFAQOQWkWhA8hAM8yh1c3d3n98fM7s7u7d7tld1r7/fzmc/u\nzDvlnd2Zeeapr6gqBoPBYDDUBa767oDBYDAYmg5GqBgMBoOhzjBCxWAwGAx1hhEqBoPBYKgzjFAx\nGAwGQ51hhIrBYDAY6gwjVAwGg8FQZxihYjAYDIY6wwiVBoiITBKRb0WkSETyRORJEWnjaO8lIi+J\nyE4RKRCR70XknyKSabefKCJbouz7BREpE5F9jmllJX2ZZu+/UETWi8hvw9oHiUiO3dccERkU1n6L\nfQ4FIvK8iKQ62mY62r4Tkd9V8btsEpFiu8+7RWS+iHQNO7cHo2wrInK9iKxy/K6fisiFlRxvnIgs\nFpE99vrPikhrR/tUEcm1+79ZRO4K215FZL/jd342wjE+EJFTReQ+e/2bwtpvspff51jWxr4m8uxz\n+VZELotyDp/av1VqpHbHesNE5CMR2WVfV6+IyEFhv99fRSTfnv4qIuJozxKR/9n9WS8iJzvaTrL7\nuMfe9g0R6VJJXyaJiNfxu20UkWvCjqUikhRl+0rvj0p+pxLHMTc42vqKyHL7d9wtIgtEpG9lv2ez\nRlXN1IAm4DZgOzAaSAaygHeBr4AU4FBgF/AYkGlv0xG4GbjQnj8R2BJl/y8AD1ajP/cDh2G9gBwN\n7AaOsdtSgM3ALUAqcKM9n2K3/8Y+lyOAtsCnwMOOffcDWtjfDwPygOxK+rIJONn+ngY8D8yL5dyA\nfwI/AKcA6YAbOA54oZLjXWT/Dy3s/r8HPOVoPwzIsL93AdYA4x3tChxayf5bAvn2b3cfsAHICVvn\na3v5fY7ffLl9TfSwr5HR9u98a9i2WYDXvl7Oq+J/HgOcB2TY5/s88L6j/Sq7H5n2ua4Frna0L7Gv\nyXTgHGAP0MFu6wR0BcQ+16nAW5X0ZRKw2DF/JFAIHOk4LwWSImxb5f0R5ZifAr+L0tYGOMS+ZtxY\n1/mq+nxONOSp3jtgJsefYd3Q+4Dzw5a3AnYClwMzgber2M+J1JFQibD9W8Bt9vdTgV8AcbT/DIy2\nv88GpjjaRgJ5UfbbB9gWfu5h62zCFir2/Fjgu6rODehtP1yH1PL/GQ98G6WtC/AtcIdjWVVC5Qz/\nwxVLqMwE1gFH2MuOsB/eMwkKlSuAHUDLsH1dYF87GY5l9wCf2w/Yd6p5roOBQsf8F8CVjvnLgaWO\n37cUaO1oX4hD6DiWpwIPAWsrOfYkHELFXrYMuMj+nkV0oVLl/RHlmJ8SRaiErZcEXAcU1eZaasqT\nMX81LI7BegN/3blQVfdhvZmeApwMvJb4roGIpANHYb2Rg/XQW6X23Waz0l7ub18Z1tZJRNo59vlv\nESkC1mMJlXdj7EsLrAfp0hhWHwnkquryWPZdCccTPHd/P/4gIvuALViax+ywbRbaZqrXRSQrrG0s\nMD9s2X8Bv4lxoj3v5BTgPVXdH7b8NaxrZ7hj2W+BWfb0GxHpFP3UKhB+rpH+S+f/vFFVC6O0IyLd\nRGQPUAxMxtJWYkJEjsISXLH8f7W5Px4SkV9F5HMROTFCP/YAJVha75QaHqPJY4RKw6I98KuqeiK0\nbbPb22OZiQCw/QR7bDvwf2I8zmR7G//0YozbPYX1sPjAnm8F7A1bpwBoHaW9wP4M+CVU9Vp7fgSW\nMC2tog/z7Jt7L9YD9pEY+h3ymwGIyBb73EtEpHtVOxCRU7Ae8vc4l6vqw3b/B2MJAOf5noD1Vn0Y\nsBV4J8wPMJaKQnQmMEFEkoEL7fnwc9kW3j/7mvnVbkdEjgO6Ay+rag7wI5Y5r0pEZIB9nrc7Fkf6\nL1vZfpWqrgNU9WdVbWP3749YLxGVMcz+fwqxtJT/At/H0P2a3h+/B3piaZzPAG+LyCHOFez+HwBc\nD3wTQ1+aJUaoNCx+BdpHcUAeZLfn298BUNUn7Iv971j29ViYpqptHNNEABF5yuGoDHc6P4LlAznf\noZnswzLZOTkAy/4dqf0A+9P5RouqelV1MZa9/hr7eO85+nKxY/Wz7PNNw7q5PxORzlWcb8hvZh8z\nE+sBlGodTkY4jheujQzD0kDOVdXvwneuFt9gvYXf71i+UFXLVHUPcBOWgDnc3md/YK+q5obt62cs\n388U4Pvwdqxr4KCwZdjXTHu7HSwB+KGq+udn28v8WkMgUCNsP4di+Y5uUtVFjqZI/+U++1qo6jpw\nnt8u4EXgTRFJquR3X2pfm62BzlhaTyzaQZX3R6RrS1W/VNVCVS1V1RexzIZjI/R/P9bL1QwR6RhD\nf5odRqg0LJZgvamPdy4UkVZYjtSP7Wl8xU1rj6peraqt7ClwA4vI/fbxT1XVAscma4ABziggYABB\ns8kaYKCjbSCwXVXzo3QhCcshiqqOcfRlVoS+elX1dSxfyXFVnNonQKaIDIm2gqouchzPabY5EsuP\ndLmqflzFcQL9rwT/bxVJS/EzAytgY0aEtgXAGBFpGbb8HKxrZ6ltpjwfOME2veVhBVMMFJGBttbg\nP9dWgY5ZGtsC4AFVDTe7Rfovnf9zT3FExoW1h5OE5TzPiPa7O1HV7VgmrdOj7M9JlfdHVdeWfzWC\n/1U4LqxghqgRbM0ZI1QaEKq6F+tN958iMlpEkm07/MtYNvv/Yjl0R4jIY/6wTBFpj/0G7ERE0sKm\naDdJVETkTiyzyckRhMGnWA/1G0UkVURuxLoZP7HbZwBX2CGZbYE/YTnTEZGOInKhiLQSEbeI/AaY\ngPVQiKVfIiJnYkVlrXM0ucPOOUVVNwBPAy+JyCkiki4ibiwfVmXH6Ae8D9ygqm+HtblE5CoRaWv3\nZSiWA/dju/0IscKt3fZLwWNYQQ3+vkbyp/iZixUE8XKEtv9iXQuviBVam2z/dv/AcubvBc7C+l/6\nAoPs6XBgEUF/Tfi5dsH6355Q1acirDIDuFVEutjr3ob9X9ra2wrgXvs3Hw/0x/ZtiMh4Eelj/2Yd\n7N/iG1trqRKxfHBnU1FIpYb91y6qcX849t9GRH5j7yPJ1l6Ox/rvsa+ZI+3/MsPu/25CrzuDn/qO\nFDBTxQkrwmc1ljllO9YDsa2jvQ/WA+dXLBPDBiznYVe7/USsh3v4dCjWg6AMy2Thn36tpC+K9Qbs\nXP8uR/uRQI7d16+xwz4d7bfa51AATAdS7eUdgM+wQk8LsCKn/q+K32WTfZx99nmvBi52tL8Q4ZwX\n222CFQr6rb2PbfbxzwdcUY43HfCFnfsau82F9dDZZS//DrgLOxIOKzhgA7AfK1prHtDLbmuDFc2X\n5DjWfcDMKP0IRH/Z8wfa18R2+1zW4Ihcsvv1aIT9nI/lb4gUNXWv/Xs5z3Wfo12wnOu77Gmq/1zt\n9iysl4xi+7ydUXo3AD/Zv0Ue8BLQvZL/eRKWUPT3YwcwB+joOFak69sfbl7p/RHheB2wQvYLsa7H\npcApjvbzsHxA++z/bT4woL6fEw118t8ABoMhQYjI+Vj+mfPruy8GQ10TN/OXWNnTO0RkddjyG8TK\nuF0jIhHDCsXKnP5WRFaISG3DQA2GhsYe4G/13QmDIR7ETVMRkeOx1MUZqtrPXnYScDcwTlVLRaSj\nqu6IsO0mrES1X8PbDAaDwdBwiZumoqoLsWyvTq7BKtNRaq9TQaAYDAaDofESsSBbHOmNFZnxF6zM\n1Mmq+lWE9RRYICJe4GlVfSbaDkXkSuBKgJYtW2Yfdthhcei2wWAwNE1ycnJ+VdUOdbW/RAuVJKzI\nlWFY5T5eFpGeWtEGd5yq/mInF30kIuttzacCtsB5BmDIkCG6fLlxwRgMBkOsiMjmutxfovNUtgCv\nq8UyrHDN9uErqeov9ucO4A1gaEJ7aTAYDIYakWihMg84CUBEemOV8Q5xxotIS39mrp01fCpWPoLB\nYDAYGjjxDCmeg1V2pI9YxfuuwBqjoacdZvwSMFFVVUQOFhF/yYpOwGKxBo5aBsxX1ffj1U+DwWAw\n1B1x86mo6oQoTZdEWHcrdvE2Vd1IaI0hg8FgqDPKy8vZsmULJSUl9d2VhJKWlkZmZibJybHWna0Z\niXbUGwwGQ72yZcsWWrduTVZWFjUoh9coUVXy8/PZsmULPXr0iOuxTEFJg8HQrCgpKaFdu3bNRqAA\niAjt2rVLiHZmhIrBYGh2NCeB4idR52yEisFgMFTBBU8v4YKnl9R3NxoFRqgYDAZDPbBlyxbOPPNM\nevXqxSGHHMJNN91EWVkZAMuWLePEE0+kV69eDB48mHHjxvHtt98CcN999zFt2rQK+3O73QwaNCgw\nPfzwwwk9Hz/GUW8wGAwJRlUZP34811xzDW+++SZer5crr7ySu+++m8mTJ3P++ecze/ZsjjnGGkdu\n8eLF/Pjjj/Tv3z/qPtPT01mxYkWiTiEqRqgYDAZDFRSWlFNQ4iFn826yu7et9f4++eQT0tLSuOyy\nywBLy/jb3/4WiMyaOHFiQKAAHHdcVSNmNxyM+ctgMBgqIWfzbtbnFbJldzEXP7uUnM27a73PNWvW\nkJ2dHbIsIyODbt268eOPPzJ48OBq77O4uDjE/DV37txa97MmGE3FYDAYKmHpxnx8dsnbco+PpRvz\n60RbiZWjjz6agoICTj31VB5//PGo6zUU85fRVAwGg6EShvVsh8uOxk1OcjGsZ7ta77Nv377k5OSE\nLCsoKODnn3/mkEMO4euvvw4s//LLL3nggQfYu3dvrY+bCIxQMRgMhkrI7t6Wwzq3JrNtOrN+N6xO\ntJRRo0ZRVFTEjBkzAPB6vdx2221MmjSJ2267jRdeeIEvvvgisH5RUVGtj5kojFAxGAyGKmidlkyX\nNul1ZvYSEd544w1eeeUVevXqRe/evUlLS2PKlCl07tyZuXPncuedd3LooYdyzDHH8Oqrr3L99dcH\ntn/wwQfJzMwMTFDRp/KHP/yhTvpa7XOL1xj19YEZpMtgMFTFunXrOPzww6u1jT/xce5Vw+PRpYQR\n6dxFJEdVh9TVMYyj3mAwGKqgsQuTRGLMXwaDwWCoM4xQMRgMBkOdYYSKwWAwGOoMI1QMTZvp46zJ\nYDAkBCNUDE2b0gLYmwu5y+q7J4bGjHk5iRkjVAxNl9xlsH017NkML55hBIuhwSAiXHLJJYF5j8dD\nhw4dOO200wLL5s2bx4ABAzj88MPp378/8+bNC9mHf5v6ykeJhhEqhqbLpkWgPuu7t8yaNxgaAC1b\ntmT16tUUFxcD8NFHH9GlS5dA+8qVK5k8eTJvvvkm69at46233mLy5MmsWrUqsM5HH31E7969eeWV\nV2hI+YZGqBiaLlkjQOxL3J1izRsMNSEOZtSxY8cyf/58AObMmcOECRMCbdOmTeOuu+4KlMLv0aMH\nd955J4888khgnTlz5nDTTTfRrVs3lixpOKNSGqFiaLp0HQqd+kGb7jDxLWveYKgucTKjXnjhhbz0\n0kuUlJSwatUqjj766EBbpNL4Q4YMYc2aNQCUlJSwYMECTj/9dCZMmMCcOXPqpE91gREqhqbN1Yvg\n5lVGoBhqTpzMqAMGDGDTpk3MmTOHsWPHVmvbd955h5NOOon09HTOOecc5s2bh9frrZN+1RYjVAwG\ng6Ey4mhGPeOMM5g8eXKI6Qsil8bPycnhiCOOACzT14IFC8jKyiI7O5v8/Hw++eSTOutXbTC1vwxN\nlu+++44ff/yRMWPG1HdXDI0Zvxm1ZC+c82ydar2XX345bdq0oX///nz66aeB5ZMnT+a8885j5MiR\nZGVlsWnTJqZMmcKrr75KQUEBixYtIjc3l9TUVACmT5/OnDlzOOWUU+qsbzXFCBVDk6VPnz4ADSoy\nxtBISc2wpjo2o2ZmZnLjjTdWWD5o0CD++te/cvrpp1NeXk5ycjJTp05l0KBBvPjii4wcOTIgUADO\nPPNM7rjjDkpLS0OW1wem9L2hySJiDdfXlK5xQ+2pSen7QOLjZfPrvkMJxJS+NxgMhoZAIxcmicQ4\n6g0Gg8FQZxihYjAYmh3N0SSaqHM2QsVgMDQr0tLSyM/Pb1aCRVXJz88nLS0t7scyPhWDwdCsyMzM\nZMuWLezcubO+u5JQ0tLSyMzMjPtxjFAxGAzNiuTk5EBNLUPdY8xfBoPBYKgz4iZUROR5EdkhIqsd\ny+aKyAp72iQiK6JsO1pENojIDyLSsAYLMNSIN998ExFh165d9d0Vg8EQR+KpqbwAjHYuUNULVHWQ\nqg4CXgNeD99IRNzAv4AxQF9ggoj0jWM/DQlg6tSpAKxdu7aee2IwGOJJ3ISKqi4EIr6WipXqfD4Q\nqV7zUOAHVd2oqmXAS8CZ8eqnITGY7HaDoXlQXz6VEcB2Vf0+QlsXINcxv8VeFhERuVJElovI8uYW\nzdGYcLmsS80IFYOhaVNfQmUCkbWUaqOqz6jqEFUd0qFDh7rYpSEO+DUVn89Xzz0xGAzxJOEhxSKS\nBIwHsqOs8gvQ1TGfaS8zNGKM+ctgaB7Uh6ZyMrBeVbdEaf8K6CUiPUQkBbgQeCthvTPEBSNUDIbm\nQTxDiucAS4A+IrJFRK6wmy4kzPQlIgeLyLsAquoBrgc+ANYBL6vqmnj105AYjFAxGJoHcTN/qeqE\nKMsnRVi2FRjrmH8XeDdefTMkHiNUDIbmgcmoNyQEI1QMhuaBESqGhGCEisHQPDBCxZAQjFAxGJoH\npkqxISHUh1AZlunmxCw3+vOXSLejE3Zcg6E5Y4SKISH4M+oTlvyYu4yPf9uCFDcw40yY+BZ0HZqY\nYxsMzRhj/jIkhISXadm0iBQ3JLkEvGWwaVFijmswNHOMUDEkhESbv3zdjqXMC+VeBXcyZI1IyHEN\nhuaOESqGhJDo2l/F7fszakYR93xaiuei143py2BIEEaoGBJCojWVoqIilm7x8vDiMnxdhiTkmAaD\nwQgVQ4JItFApLi4OfNcXzXA8BkOiMELFkBDqQ1Pxk/DUmOnjrMlgaIYYoWJICPUlVIZlunEX/gK5\nyxJyXIOhuWOEiiEh1If5a1imm49/24LkfT/Di2cYwWIwJAAjVAwJIdFCZd++fZyY5SbFDQKJy1WZ\nPg5+XgJ5q+J/LIOhAWKEiiEhJFqo7Nmzh083eSnzggK4UxKXq6IKnlKjGRmaJUaoGBJCojPqd+/e\nzdItXkbNKKK0ZWbiyrSUFgA+8JYak5uhWWKEiiEh1IemArB0i5eSFl0Sl/xYsjf43ZSHMTRDjFAx\nJIREC5Xdu3cHvutFcxNyTADSDgh+r8zkZsKODU0UU6XYkBASXaZl165dge9xF2R+4XDZfEjNAFxW\nvTFTGdnQDDGaiiEhJFpT2b59e+B7wgcGE4GkVEugRNJIpo8z0WGGJosRKoaEUJ9CxTfzgspXrktT\nVN4qUG/V6/k8sDfXOPINTQ4jVAwJIdFCJS8vL/A9oXpKaYH9uS84Hy48SgugvAj2bDYRYoYmhxEq\nhoSQMKEyfRyeZ8ewdetWDjjgAPuY8T1kQHAsf8Gx0Acf3QvbV1cUHjs3BFfzlsFrvzNOe0OTwQgV\nQ0JIpKbyc34JHo+HXr16MSzTTYtfV8RPG8hdFhQc794W2vb546B2YIIzvNhbHlzHnRIaMWYwNHKM\nUDEkhEQKlR937gdgTL92fPzbFrRyl1fPzFQdH8umRUHBUSGyTUO/+8OL7d8Cd6oVIZaaEduxDIZG\ngBEqhoTgFypebwxO7NpQWkBGWR7DMt2M6CpW7S8hfomIWSNA7NvIHSFCX9yAQKd+wfDiWBz5BkMj\nxQgVQ0Lwl2mJa56KbYoa0nYvH/+2BSltD8Ljs30qLnf0RMRIznSITWPpOtQSGG26M989Cp9TOREX\nuCpJBfOXctn6tQkxNjQZTPKjISEkRKjYpii3QIpb6FC+1dEokbfx+0TUZz3g2/eC/B+gZQfL17F7\nkyVYLpsfspnX62XRokUsXLiQlW99z8/ubow6qRNjU0MDA8Rban3JWwVPjYBxj4Ye31tqbSB18H7n\nTMI0GOoJI1QMCaHWGfWxPDCzRqDiwuPxopJk5SC6bPOXz2MJnfAMd6dPxFsG+3YEw33FBUlpIavv\n2bOHv//970x9/N8U79mJiNCrYwtaDD2ctq4iFHDZ8sunvsC8ArJvR0UTnLiCxzcYmgDG/GVICAkp\n09J1KHmph3LPp6XMb38VP7QaapW+V6LX4QrxiaRAUX6wTX3g8wb6/cQTT9CjRw/uv/9+2vY7gX43\nPMOn325mw0Mn8OxoN/vcGQjW8VTBh+DBbR1fwVe4HUoKQo8//Pqg495gaAIYoWJICImq/fXpTx4e\nXlxGt2PP5de0HoyaUcSelIOj1+HqOtTSRvyRWA4zlA8XP5a15ecf13D9mUO54YYbOOqoo5jz0VJS\nj5tEYYuDufKl1eSMnEn2tc9z9dC2lkbiD+5C+drXK7DMp8rXn7wamjdTGiZkYsUUpDQ0UIxQMSSE\nKoVKXTwkp49jQc53tGmRxKBBg3C5XCzd4mVPcufKCzu6kqxaXWD5OAAf4FXo6cqjq3s3jw36gbef\n/BMffPAB+e52AYd8ucfH0o2WdtPloMyA6csvWNpSENBe3AId3UUhh97/4xKWl/ekxOuykiUNhkaO\nESqGhFBrTSVahJYDVWXBJh8n/eZ03G53MDfm3OdiO0aYv8ONDxFLQKS44bQjDkBEGNazXUB4JCe5\nGNaznTVTnB+6P4GWyfZXAS9CO1dQM1HAlf8d2e7vSZVy9PO/G8FiaPQYoWJICLXKU3FmrVeSxJiT\nk8PP+SWMGzcu5Jj6yuXWCtEqBpcWQNn+EJ+LTwn4Q1TBQxJXLUonZ/Nusru35bDOrclsm86s3w0j\nu3tba6Mwn40Amb2PtE1fUE4ye2kVXEEhyU6QDLhV1r1V/d/HYGhAxC36S0SeB04DdqhqP3vZXKCP\nvUobYI+qDoqw7SagEPACHlUdEq9+GhJDcrL1yl5WVlb9jcMjtCJEcY19fCH8egotMnM4++yzAYd2\n5PdhxJgL4lPF44Vp69pz/uHCga79POo5nw9KuvPZs0uZ9bthvHvT8RU3jGRia9UJcOESH4tajiY3\npQeX73rcriwgLPIezklJa7Bn4fAzYuojpQXWKJO5y+pmzBYTjmyoI2LSVMTiEhG5x57vJiJVXckv\nAKOdC1T1AlUdZAuS14DXK9n+JHtdI1CaACkpKQCUlJTEtoFTqwiP0ArTCHI272ZdXiFrup9Pxwv+\nwk+FljCJuTSMuCGlJWxahKriEiHJJdxxlNI7KY8DpZB7kv7LYPkuxIdSJe5UGDgBkq2w5FNlGVeM\nP802qQl62t94seRU5vzSmcJywdvvPDjl/qr3G6PmZjDUB7Gav/4NDAcm2POFwL8q20BVFwK7IrWJ\ndbefD8yJ8fiGJkJxcXH1N3JkrUeK4lq6MR/1KdnuH7km9V1++uZ/gPXgHpbppu3e1daDt7TAmvzC\nKmywrLz0XpR6wetT3EkppCRZt4dLIJkyhrnWhfpQwgl/uI+Zan2W2875fXnwwmmBZtdRl/PCTWM4\n76BttEjy4V31Gt7NS6r+PSJpbgZDAyFWoXK0ql4HlACo6m4gpRbHHQFsV9Xvo7QrsEBEckTkysp2\nJCJXishyEVm+c+fOWnTJEE/82kKNhApYRRcP6BrR1DOsZzuOZD1zUh5kctIrjF91FeQuo33JT3z8\n2xa0T/VYb/TRO4enrITnH7odAJeI5ePoETRxuQUOSisN9aGEE/5wL863yto78Yaa/2TzYpLdlmbk\nxstjd13N+U9+Hr2vUKXmZjDUJ7EKlXIRcWOXXRWRDlhRlzVlApVrKcfZJrIxwHUiEsGAbaGqz6jq\nEFUd0qFDh1p0yRBP/EIlZvNXNXDv3szonx4lBQ8uAZevDFbOoXPx96EFJf0EBtIqgPJiwIfLW8rt\nfbeR6l/f54Vflgc2EeDSrILoAgUqPtyzRlQsa+9OgeQWwcrE6UGtxy3Cxm27Wb74Y3I2RVTyLboO\n9d+JoZpbDBFyBkO8iVWo/AN4A+goIn8BFgNTanJAEUkCxgNzo62jqr/Ynzvs49aBJ9JQn9RaUwHL\nVBUhl+W+++4jyVMUtlTZ2bJPMKPe5Q7dz/IXLL+Er9zyo7gsbUECmYspFZ3mh59Zef+6DrUEBgKd\nB1jzqRnWMlcytOoMk96Bgwdb7QB5K4PbCxw7fDDeLoM4/6nF5GzezdjHF3LcXz8hZ/PusIP5AA0K\nFONnMTQQYhIqqjoLuAN4CNgGnKWqr9TwmCcD61V1S6RGEWkpIq3934FTgdU1PJahgVAnQiUCy5cv\n580330QGTgiG5bpTYOBF7Gp5KKNmFLGjJAmOvCR0w29mVKy5JY4vE98KdZqf9jgMmVR1hw4eDKmt\ng5pIaYGl9biTod2hEcx3YUEE4kJcbjxe5cGZH7E+r5Atu4u5+NmlEQSLg3A/ixlN0lBPVCpURORA\n/wTswDJZzQa228sq23YOsAToIyJbROQKu+lCwkxfInKwiLxrz3YCFovISmAZMF9V36/uiRkaFrU2\nf102P/h2b+Pz+bjhhhvo2LEj59z8V0sjcKfCpPnQdSgiwtItXnYUeqDzwND9bV2BOvplFQpOtrZ3\nR3AXduobe187D7D669cevKWWsz5SSZbOoRH1azQLAJcoixYtxGfHQ1eMOnMBYh1j+jhYOC3YZEaT\nNNQjVWkqOcBy+3Mn8B3wvf09p7INVXWCqh6kqsmqmqmqz9nLJ6nqU2HrblXVsfb3jao60J6OUNW/\n1OzUDA2JamsqMfgHZsyYwdKlS5k6dao1Hr2/3IqtDYSEFOetCOuQj1+LrBwW8Tvme59qCQD/OCfO\nMedfGFd9k5JTewArryQcRxa+iIs+rcvJbJvOjMuO4qDSXHyeMkBDo85ylxEwf71gJ2/6I8z8EXK7\nN5kxWgz1QqVCRVV7qGpPYAFwuqq2V9V2WEmNHyaig4amgV+oFBWF+z4iEIN/YOfOndxxxx0MHz6c\nSy+91FqQDr5nAAAgAElEQVTYeUCINhMQKgD7QiMDFfjLwmLKfY48lladgit4y+DDu0PnV86uuu9O\nnFFa4grVHvz+Iadz353KpIsvZfHvRzLi8C58OPsp9r77KIVLX+Gx07OCQQLOfvjL9fuJEiFXJXmr\njBAy1AmxOuqHqarfPIWqvgccE58uGZoi/gf33r0R3tYhVDNxvuF7iiuE6qoqV155JXv37uWZZ54J\nDAAWTkCoJKVDcag/QlXJ2epjW6HPWi+5hZWo6AzVrWAGq2aJemd+Tad+kceidwqAsBycgw46iCWv\n/QfPN/O4+6qL2L9/f836EQ1T6dgQB2IVKltF5I8ikmVPdwNbq9zKYLDxC5XduyM4m8M1k/R2oW/4\nYaG6L774IvPmzWPKlCn069cv6jFdLhfDMt10TCmG/b+GtAnKnad0oMeBtuBwJVVMsjwgM7iBPzu+\nuvjza5wCJdw/lJphTRE0jL59+/LSSy+xatUqJk2aZBXkdPbDlQytOoZuNH2cVcusOvg84Ck1UWOG\nWhOrUJkAdMAK730D6Egwu95gqBK/UNm1a1fFsinhkUvF+dbD3Z1qfToetuu27uOGG27ghBNO4JZb\nbqn0mAfu/4H/TWzBQa2AXzcE+qGq+FQ4qW8n+7iu4EPemWSZmmG1uVOtUOC6qLEViTCzXThjxoxh\n6tSpvPrqqzzwwAOhjXUxwFfuMssn4/clGcFiqAUxFZRU1V3ATXHui6EZ4PF42L9/P61aOar1+n0P\najufs0bAD59YTnf7Df+Cp5dw6y/5TPzPalq0aMnMmTOjmr38dN212E5mDD54RcTOS4H0wk3W8dDo\ng2WJhDj/qyRaQcZaFmq89dZbWb16NVOeeYkD9q/jpnR/UqcnNACgtMD6zdQbu7biNC9GKdhpMMRK\nTEJFRP5HhYB6UNWRdd4jQ5PEqZ3k5+eHChW/2WnnBujQJ+IDraC4jGtW9OHnpL28N/MpMjMzK6xT\n4cEd9havqnakl9gXs+OSdjq765I6qvorIlx198N80u4L8tyfAXYYtPgsweIn71vo3N/6rl4rgq2q\n/JqQYAFT9sVQO2ItfT/Z8T0NOAfwRFnXYKiAU6hs2bKF7t27h66QmhGimQCWnX9vLuu/WsC6rfvR\nw87h4MPOolVWdFORk1/aj6DztgWkujVEW4EIru7i3ZbZxykESgvscerjcKk7S9fHKHhycguRpGT6\nuX8GHDKzyBnZprDrp+DsO7aBoTLB4q8E4PNGH3bZYIiRWDPqcxzT56p6K3BifLtmaEo4hUpubm7V\nG/hzL/Zspsf8CznS9QPZ7h+5KjlYhbgq9rbuw0kvFrGlwEf5wRUflPtIC44X7/OEmoH8wQOo1Y+6\n9DPUsKSKNeKk0J49gWUKaFiRSv+QyAHWvVlxZw91hZ8dFZHDcnwMhpoS63gqBzqm9iLyG8Ck7Bpi\nRlVp2bIlAD/99FMVaxPiJ3D5yhnvXsSslCnckvQq47+9JqYHsYjQr6OL1qnCFz+X43XkpCiwxtON\nElLwqAufKznU7BOeuFiX5eVrWLo+u3tbXrn6GLp0zQoutF1CIbQIK6xaVc0yg6EOidX8lYN16QqW\n2esn4IpKtzAYHKgqGRkZdOjQgX+/uoAVB57I3KuGh65km7v++M/nuHyfm262EPDiom16Cinl5bhF\nwVcekzO5685PeOY0a4Cs4/maPSU+2qYHC0t+TyZ/LbuI4a619DhyDOc69+cMHogQ1lwrnPuupg8j\nu3tb6NkVthHQsrwIbvUPS+yClu2g8BersbKaZaqmqrGhzok1pPhwVe1pZ9j3UtVTga/i2TFD08Kv\nIWRnZ7OvpJxf9hSHFkh0mLv+lH8nCzfs5vgX9vPCGje+S+dx2sTJuF2OCsIxPIh7bnkDCEZ/tUlz\nXO7i4i3vsazUQ3jON44eR54UurEzZ+XyD+rWLFTFoGNVYme++30qhZrumPfBXket1mg1y3wea12/\nCS4efiNDsyRWofJFhGUxDFFnMFioKlK8iy5pBbQefXPFyrsOc5dby9jsbY/7oL5cfOrRtOh9Yo0e\nxEleZ0ithIUWu/j91Vdwa/uvmNX1rcjjpFQyMFitqc2+bXOWX1MpJTm03RkeHc1n4y13fC+zxpWp\nbsKkwRCBqqoUdxaRbCBdRI4UkcH2dCLQIiE9NDQJ/OG8rboPRNyW1TWk8m7aAaiCxwflJLOy2wR2\nnnQ/63wHB3dSzQfxvpZZwZnwcK+kFmR3b8t17b4mO317zU+sPrDNWSLwj5bXkecNujcVQrWOCGVu\nyF1mmRD9uJKos9IvhmZPVZrKb4BpQCbwGPCoPd0K3BXfrhmaEqqKAOMP/hV8XtTnJdltVd5VVd5e\nV8zZc/cz5cdeXFx2J19zGOW4WVrUxdpB2HjysfBjj0vx+PwCzR3aeGBWnZxXfXPTHVMY1CkZsCoF\noBESysJNheFC5siLImfmm9pghhpQVZXiF1X1JGCSqp7kmM5Q1dcT1EdDE8B6sEN2ix3cNbwVexfN\nJHnxk7z46L0MGjSIM/6Rw/pflZ4Zygr64MJLsvgYdt6tNT7mAUeczJhZReQW+GDcY3V4Ng2I3GXg\n8wSSOr0+xeeUKu4IYcIhQkZg4EWWLU19xmlvqDWVRn+JyCWqOhPIEpEKd7eqNtE71RAP/C/DV51z\nCh2kkD/+8Y88//wnDBw4kGcn9ee3mZtITvqJ2W1ac9SutzmjUz7Z3c+q8fF69+7Nh5e2wCUSOmwv\nRB7bJJHUJtPe+eB/8QzIsEyE/ki529cczoVdtnFM291W7klVbF+LNT6LvT+TAGmoBVWZv1ran62A\n1hEmgyEm9PuPEE9ZwIQ1fvx41q5dy759+/j888+5ovdukt2W1Hn3puO5P2tt7X0ducssgQLwzczQ\ntsY8MmJ4ra4yR0l8l5sth57LhKTHeH1Xz6q3R0OTI6uRN5MwjBmuUVGppqKqT9uf91e2nsFQFVad\nqmps4CxjUtO3ZufD0ecNfk9uESwHU5nGUEd1u+qc8DwXGxFwq49h7g0oSnFhGT/vd/OXKf/myosn\nBCPcwn0sh58JP35ifTe1vwy1JNaCkh2A/wOynNuo6uXx6ZahqaFaRXxR5wFW2ZCUlsEyJuoLmmNq\nQnihRI89lPHBg2u2v4aCP7y6ZC+c86y17LlTUIVyktitLZibOoXkgyxB+o/SPzHx2XJu/91vLcHi\nr/VVXgS4rGiyD+4Mrf01fVzw/zAYqkGseSpvYpVlWQDMd0wGQ0woWrWmomoNFLVyTo3KmFTAqeGM\nfhjEbU3OUSYbK87wavs8ReDujL8wPnkpyXgRsZYl4yFb1wTDt8EOIyaoPpraX4Y6Ilah0kJVf6+q\nL6vqa/4prj0zNCmq1FRKCwCfVQzxm5kE1q6NOcYpNN673Y5w0hoVc2wsPHbblRzVendAVqiCD2GZ\n73CG9WxXv50zNAtiFSrviMjYuPbE0KSpkDsRjnM8E58XWnUKzZ4vLbC0mGiDaUVi5ezgd28ZVoST\n1o0W1JBxDi8soCrsXvAUkh+hkGddm7eMU73ZE6tQuQlLsBSLSIGIFIpINe5uQ3NHVTnyIFfkcdBz\nl8G+vOC8uKwHo9+84/exeEvh1++roV1E0o2UOtGCGjKOMWkESHILp3TczVlnncWOHZUMRla2P7pA\nMMLCECOxjqfSWlVdqpquqhn2fEbVWxoMFoce4GXGWemRx0EP1xacJef97TXRLgZOCH53hdXHatW5\naeVjdD/OmiIg7hTOuvlRdu7cybnnnktZmSMMOZKQNxhqQazjqQyOMB0iIrGWzjc0JWrw1nr4gV6S\n/ZVSwgVDBW1BQ5MT/SG0UD3toutQrEtcYPCloftv7E7py+bHFvLsToWJb9Hn5Et5/vnnWbRoEYef\nfSNzy4/HpxIU8tGqFFemvRgMEYjV/PVvYCnwH3taCrwCbBCRU+PUN0MTYu0uN+X+VJFwwdB1qPXw\n8+NODU1OrE2p+NRWkNo6VGsRV+NOfqyKy+YTuLVd7sDvNWHCBK74/YOUH3Yqm7U9Pr+ny1NqVS1W\nX6jPyl+6pTp+LEOzJ1ahshU4UlWzVTUbGARsBE4Bpsarc4amw3e73Ty0qNQK6R39cEXB4C8nYr9Z\nh4xVD7UvQ+/UWjr1q7j/xo5Tc8ldRqDsin8oZFu7HHTqeYg7mV2agRtbbuBDfeWAHRmXuywYjedf\nVlpgVUMwWouhCmIVKr1VdY1/RlXXAoep6sb4dMsQoIk4SA89wMOdI1JBvfD+H6Lb8eNplhKxtJSm\nJlDCCfc5OeaHH9IBl8vFgRTgQxCxI/P84Xnqs9Z3mh/VV/+10gyNhliFyhoReVJETrCnfwNrRSQV\nKK9qY4OhbztfdJ9Koug23JqaOuE+J8d8dve2HH5QazLKdlDqVWv8Gk2ijCQ86sKH21rfaR4UF+zf\nGdsgXnmrqj1EgaFpEaujfRJwLXCzPf85MBlLoJwUZRuDIcCaX12UHwJJLmJztjfUuluNAWcZluQW\nFTS/d286Hv5yAZt2ePm4sDvvdLwCEWG4aw09MuDcrkNDtblO/WD3JqA0vv32a+Tmv2/UxBpSXKyq\nj6rq2fY0TVWLVNWnqvvi3UlD42dD3j7+7+3ioM8k3MTVeYDlbzHUDf4yLK4I7425y6C8iKy2bi7N\n/AX1evjaewjPecfSo7Wv4vrVMRf6PNULU24i5l1DkFgLSvYCHgL6Amn+5aoapba2wRCKqrJquy+6\nz6S0wHoAts2q2wN3HlBxWXN5ExZ35PN3mB5TXMpxe+fzv2+Xc1XmRrL7dLIaIkV8qbdy05YtrAAz\nLkszJlafynTgScCDZe6aAcysdAuDwYFSSel7Z8a8P/oonFjzMupqu6ZMmOnx1rv/ysnlC5n2Wg6v\nLd8W/D/8FG4L5rFEy2eBiuO8NMUSOIYqiVWopKvqx4Co6mZVvQ8wOqshZrSy4l/OjHl/9FE8MAIm\nIq6d65n5f4M4OjOFS/6zks0LZ4VWNdj1Y1AD8YcoRyJ8qIGmWALHUCWxCpVSEXEB34vI9SJyNtZo\nkIaGTk1s1nGwc/duJ1xxZHLkN11nxry4zMMo3oQL7XdvI13389al7ejSNo2rH5pJBM9K9O39+AME\novnNDM2C6hSUbAHcCGQDlwITK9tARJ4XkR0istqxbJCILBWRFSKyXEQiXnUiMlpENojIDyLyhxj7\naKgJtREgsW6bu4wHTkji2qNSIr/pOjPmL//APIziTbjQ9nmhZC8dWrl575aj+CpP+NOnHjSgXkrl\n2zuJx7gsTWH8m2ZErNFfX6nqPlXdoqqXqep4VV1axWYvAKPDlk0F7lfVQcA9RMjGFxE38C9gDFZg\nwAQR6RtLPw0NlE2LSHJBkksC8xWobca8oSKqkR/GFX5jDUSJ9erUkrfffpu/LynlsjeLKSlX6Nzf\n0kAgYohygOnjYstlqQ5+/04THf+mKVJp9JeIVDqOq6qeUUnbQhHJCl8M+OMTD8Aq/xLOUOAHf7a+\niLwEnAmsrawvhgZM1gg8Pmv0xyS325i3EkHbLCtSy/8wbt8rGBpc4cHsCjFLDh8+nNeuG8zpj3/F\npj1FvH9XS9IqC1GOJ5EqVJsXjwZNVVfIcCAXmAN8SRWD98XAzcAHIjINS0s6JsI6Xexj+tkCHB1t\nhyJyJXAlQLdu3WrZvQZIaYFVIiN3WcO5mUoLYOcGq1hhLOO9dx3K3f/zkHWAcv0ZAxvOeTRlnGVV\nvGXWIGj+6yhcU3TZBTaLdwcWje7fgRlnpXPx68Vc+NQKXh2fFHOmdJ3i97epzzj/GwlVmb86A3cB\n/YDHsQpI/qqqn6nqZzU43jXALaraFbgFeK4G+whBVZ9R1SGqOqRDhw613V3DoiGq/s7w3/KimCvY\nrstXXlxZ3vTrbjUUnGVWXG7Ytz14HaWHDSs89tGI/8uE/sn8Y0wab36znYc/3W35WGyN5oKnl3DB\n00vieQYWtalQXd8008TOSoWKqnpV9X1VnQgMA34APhWR62t4vInA6/b3V7BMXeH8AnR1zGfay5of\nNR2cKp44+wQxFxpUBXG5TEhvorh6kZX42KY7HHkJgYqR3jIozg+ul9wChkyKupvrh6bw5qQu3Hq0\nG59SeUhxvDD+tppTD4KtSke9iKSKyHisZMfrgH8Ab9TweFuBE+zvI4HvI6zzFdBLRHqISApwIVCp\nb6fJUtPBqeqCaGXOnX2CmMclUbTWtlNDNfE/jAdOiH4dRfORODTQ07vvIy0J3C7Bq1R4uUmY1mJo\nFFTlqJ+BZfp6Fytqa3Vl64dtOwc4EWgvIluAe4H/Ax63R4wswfaFiMjBwLOqOlZVPbYm9AHgBp53\nlt1vVvhV/5K9cM6zDeNNzd8nv08lBnNWeXk5eYU+UpOMWKkXIl1H4rbKrkTDoYEKCiJ4fYoCXy78\nmMJWR1FQ4iFn8+7o+zA0S6ryvV0C7MfKU7lRgnU2BNDKxqlX1QlRmrIjrLsVGOuYfxdLkBlSM6yp\nIQgUP6kZwUG1qmD//v2cc845rNhWzr/OaBPnjhmiEn4d+XNQ/LXWLpsfqpmGlb4XBJdY/pShZZ8z\nIO91XtJTuPjZpXTOSMPjU3Jad6p4cxuaHVX5VFyq2tqeMhxT68oEiqEBUY+JY/n5+YwaNYqPPvqI\n/5zdlmuHVVKEwZRQSRzOkSGj1VoLL30PiAj+F8sx7q8AKCv3sTm/iC27i7k49wxyXP0iF7E0NBti\nzag3NEZqGj1WWmCVL6/F2OQ///wzxx13HCtWrOC1117jd0e1rPG+DHWM0yfirLUW7QUkNSPU9yLw\nnuco1OvB8pZZlKubpZ7e8ey5oRFghEpTJtboMadTPpaKwVWwdu1ajj32WLZt28aHH37IWWedVcMT\nMMQFZ7BFUro1H/4C4nyhKC2AzKMAAXHzTMaNzHOPwrPhM3a+/0TAN5MsXoYlfVfzfplyLE0CI1Sa\nMjWJHqtlxeAlS5Zw3HHH4fF4WLhwIccff3w1O22IO5FyP8JfQPbtCK6f961V/l5ckNKSN92n4JMk\n+h0/liFtS8mb9QdaFeYyq+tbZCf9VLM+NcScLEONMEKlKVOTxLFaVAyeP38+o0aNol27dnzxxRcM\nGGBs6w2W8NyP8BeQw8Y6VlbYtRFUyfH0YH1eIWUeH9/nl/LgU7M4ts9BrH3yWr5e9pWVHFkTbaMh\n5mQZaoQRKk2d6iaO+QWRO9X6jHG7GTNmcOaZZ9K3b18+//xzevToUYtOGxJO+AvIwPDgTQWUpZ7e\nVhIkVgjoiq1FvPPOO4wePZrnP15HeWlRzbSNWLXqZpql3piol3I+hnrGf1NeNj+yU94fMhxjSZVp\n06Zx++23M2rUKN544w1at24dh04bqk11o+kqhK8LgUx8cYHCsKTvcJWBTyE5ycWwnu1IT0/njTfe\n4OVLD0Yot7arbvHHROdkOe+BeNEQ6/YlAKOpNDWq8yZXS6e8qnLHy+u5/fbbOf/885k/f74RKE0J\nv+bQprv1wBchO+knDuvcmsy26cz63TCyu7cFIDU1lQnHH47XB+VepVwlRNuIKes+klY9fZwVSNLY\naMY+IiNUmhvOm7QWTvkyjzLx1d088v5GrrvuOmbPnk1qaiUJkT6PpRE1o5uryXBA1xCttXVaMl3a\npAcEip+klgfgTm3JGxtTOf65vdzxz1cDA30VlpTzy57immXgl+2vW8ESi6CqrZmtGfuIjFBp6NQk\nKTBa3S4//pu0Mqd8JQ7XguJyxs34lf9+U8QDZ/fmn//8J263O/rxcpdZhQi9pc3ura3BUJPryJ91\nX1pgbZtSda5RUlIS5448isFnXMUjjzzC5ZdfzrKNv7I+r9BKkHx2aUCwXPD0Ejb+srVphhHXZ92+\nesb4VJozzjpeHfoEzQ6lBZYQ8KvujsixrVu3Mvbhpaz5pZTp57Rl0rhDQaqo6eV8SzMDLTUOYsm6\nD8f2z7nKC3niiSfo2LEj9913H2vpiq+TNSRSucfH0o35TH1/PcnblpPFT7BHKw4k1thpiHX7EoTR\nVJo7kZzy4QM82UJh7dq1DB8+nB93K4tu7cekgUmxZd0347e2Rku0rHtg7lXDmXvV8ND1w/xzsuUr\n7r33Xp588klWfvQK6ikDgs59gGzf6oAyhLeswjAKFzy9hDXbwoZWyFtV9z6WeEWUNdOS/UaoNCQa\nSriks5igLQQWLVrEscceS1lZGcvf+DfDWm2J3cHfmAdaaqpUZQ5zCv5Y8pXC/HNzXp7FBU8v4eqr\nr2bm3x9gx8v3UPzNWzx+5iFkd2/L9oISFpUfRhnJeNSFz5VcYRiFwpJyXi4dTo730NBjle2P/30S\nD+HVTDBCxVCR1Axr8CZbCLz4yTpOPvlkOnXqxJIlS+iTuqP6Dv5m+tbWaOk61LoGkGC+UucB0YtF\nhvnn1qQE1zv33HOZ/8LjFC+Zw+/OHsXLH3/F5vwivtbeXFx2F3/znMvr/Z8M0ZZzNu9mfV4h/y09\nnovL/kCOp5K8p4byMmYAjFAxRMOVhLd1Jnf881UmTZrEiBEj+OKLL8jKygobqEsqDk9raBq4kqz/\nORY/R1jS7PcpfUOaR44cycKFC1FVrvvz44GosK+1F097T6fHkSeFrL90Yz4+BR8uykliabo9tl/Z\nfqvWmNEiGixGqDR1Ipk51GvdnJVQWOrj7CdyeOSRR7j22mt57733OPDAA63GrkOhrf/NUeH9PzS9\n6B1D9XH45yKFEA8aNIgvv/ySzlKAz1OGqg83Pv7caWGF8ORhPdvhEnDhJRkvw1rU4YjieatCr39j\n6qpTTPSXoQKbTvoXZwzvzdod5TzxxBNcd911FVfyeYLfTUSXwUFOcSfW5xXiU7j42aUhSZJdunRh\nyTuz6TP2cvaVw8g227jgXCth1u+UPwLI7t6Wwzq35qhdb3NG0pdkp9dxUq3RduKG0VSi0UzttJ9+\n+ilDhw4ld4+H9ya2jyxQIKIz32AAWFrUJVAfzB9C7KRVq1YcM+FGOqYr895fzLn/+pp9+/axztPF\ncszb2s27Nx3P/VlrQysf+8PFnC81hgaFESoNmZoKthpkIKsq06ZN4+STT6Zdu3YsvaYjp/RKi75B\nmDPfaCnNgBgTKIe1+AWXnbrkDCF28so1x/H9J6/w+IS+vLViO0PHXMCd+y+0HPN2gmSFkGJn7kx5\nEX/853NVl36JRLjpt2x/leZgQ+wYoWKgsNTH+U9+w+23385ZZ53FsmXL6NMhueoNXUkmostQgez0\n7RHrg0XixlOy+ODWoexKbk+Zui3HvK3dBEKK/ZFfYVGGh+5fYfltijvF83QM1cQIlYZGVSVWqqK0\nwLIXl+6LafX12/Yx9N87eD0nj0ceeYRXXnnFFIU0WHQeEFNplkhEqw9WgbxVnNxuK//5y+2Iz4N6\nPYj6aJueHAwp3n+LJTgcJlYfwtsFPa3SL7lnGMHSgDBCpSnhz2oGwAfLX6h09ddee42j/vw5+UU+\nFkw+msmTJyNVlVwxGOLA6cP78WD6bDp89wa5/72Dp1+chU9hkPzAFe75/FToCtGIP9XB5Ph6A1Cu\nbpYWdamvrhvCMNFfTQlnVjPAu7dBp76h5qnSAjxe5YstpUx79gKO6NKaV89LI7OPyTUx1B0VyrjE\nwCU9i7m4RxFThl7CA0/+lzEXJTMrfSrJlOMqdsHy4wLrnihfk+36jhxfb45ybeAs/RhyT6+eKdZf\nNVu91nxdj3sSz7FaGjBGU6mM2pqiEoHTmZ81AmtgJZvwbPfcZfjyVuF2wYhubhZe3ppFtxxB5gHm\n3cLQMBAR7r77bl57ehoDNvybZC0jSRQXXlj3ZmA9F8rpGRv5TcZmZqdMocv+NbFXwA6UpC+FFxz3\nt6mgXScYodKU6DoUDuwZnHenBuzQqkrOa4+jPkVEEBGSRUkeeHb00hsGQyJxRJeNGzeO/UddQ7lP\nKfcqZV7wdDgiZPXSlDacmLrBEjgQ27glucsIjGYJ4C13fG9e457ECyNUmhqtD7K/SCDUt7CwkEsv\nvZTrH52LxyeBEhmV5pdUVufJYEgA0/54B0mdB/DBT8oJL+znoX+9GLx2gauOasOE8y8OzJdqEn9c\n0SYwH3G0yXChIY5HoMm3qhOMUGmqiOXYXLJkCUce2pk5s2cx9sp7STq4f9AZP2aqCQc2NGiSCzdz\nWq9kbj9zIB+u3U2xB3x+wZI1Ilj40p3KA+0eqlBzrALhQmPcY8Hvox8290MdYIRKE6XMo9x9990c\ne+xxFJd5+fT3w/jTn/6E2+XwufhrdpUWNM3R9wyJp/OAuDiox2d3Jum3L3LFumHc81UGD33uZV9b\nW4B4SsHn4fuUvqzdVhDQTiIOYewUGsktrEAWP6aGXZ1ghEptaYDlXNbs8DLs2UKmTJlC1vCxvH7D\nIEb0totBhg/AtXK2FYbsH+XR3FSGBkhOcSc2l6TxZe/rWd//RopS2nD5qf346quvAusUlpRT6vFR\nWFIeKJ0fPoRxCD4vrJwTnDc+lTrBhP1EI7ySaSPA5/Px91c/564P9pORKsybN4+nN7Xl3V1vk1Sc\nTzZUrNmFOKJhTGHIZk0DDoH11xMbLD/y39RHSBlSRpkXTj3vOMYc4uKUY/oFiliu21bI05/9WKH+\nWIVETG8pfDMzOO9KCprHGuK97395bcD/ExhNpcnwww8/MHLkSG57dy+jD01i9bUZZA46PpiVnHuG\n9bYWGBvDduQPnBB5qF9jEjM0IPz1xIa71pBMOW6BtGQ3144dwN2flPC75YcEfC0KfLxueyC4PqT+\nWPj17Iz+ckaFNRXq4T42QqWR4/H6mPb+Rvr378+KFSt4/py2vHFBOh1bwk/f/C840JG6Q6vF2o78\niEP9+jPzjUnMUB0fSYwFJ2uCv55YD9kWeGiJO4UL//AEM85uQe66lfi83kB0mALtW6dUrD9WIfrL\n4WP0ea326eOshEh/UmRjpZ7uYyNUouEfYa4ypo+r+zEZqlFh+Ntvv+WYKUu4/eX1/OY3v2Ht2rVc\ndi7RApgAAB5FSURBVGSqHd2ljP/2GrJd31kDHYk3tFqs+oIXWfhQv87M/Gh25svmm5Djpk4chURN\naJ2WzDlJi3GJgrhh4ltIt6O5dGAqK87ehufb98DnRX0+kl3CU5cMYfHvR4aavcKjv465Mfg9Ukhx\nY36hiuU+jgNGqESjob6lTB9H6TOjuffeexk8eDCbfi1m7tVH8sYbb3DwwQdbb1s2Ll85p2ds5NLU\nRczq+pZ1c5UW2K0a/e3FOVywid03NBDmZm8I1otQL2xfG2jr1sbFcWf9Fu/PX1Ow5CW2z72bpe/M\nxufzhe4kxF/oglPuD876NfXAPULj1tTr6T6Om1ARkedFZIeIrHYsGyQiS0VkhYgsF5GIHmER2SQi\n3/rXq5MONcAorQqUFlhvFpUMQLTou10Mvv9z/vznPzNhwgTW/eV4zh96UDD3xOUOruxy80PLQZyf\nuoTs9O3WsvDor0hvL5FMYrWlgb31GuqB2t6DjjItkeZFhGNHjmbV7IcZdmgnrr/+ekaOHMkPP/wQ\n2/7913ks90hjIB73cQzEU1N5ARgdtmwqcL+qDgLuseejcZKqDlLVIXHqX8MiUGFYobyowtvRjh07\nmDRpEsc/vJR9pR7mz5/PjBkzaNcqpZKdCpfvfYKs8o3BRbGO2BhuEjMY6pvDz4w8n9IypER/165d\nee+993juuef45ptv6NevHw888AClpaVh95UvshbSlEY1rYf7OG5CRVUXArvCFwP+8KMDgK3xOn6d\nkajIifAKw/bbkdenPPnkk/Tp04fZs2dz57hDWPvg8YwdO9ZaLzy6w2H+wuehZysvLVMc2kt49JcR\nGoaGjqp1jXfqS6BgqjsVhkwKWW3uVcMD1ZFFhMsvv5x169Zx5plncs899zBw4EA2fjIjdN9hWsgF\nTy9hzS5HFJi5R6pNon0qNwOPiEguMA24M8p6CiwQkRwRubKyHYrIlbYpbfnOnTvruLvUOHIiYt2h\nynDaP+355T/tYdiDX3DttdcyePBgVq1axZRz+tAy1U4vihTd0aFPcB/ulNC3LifiMjeLof6pqsac\nzwP4gtd4crq1PCk1pt0ffPDBzJ07l/fee4+ysjIu/uNTlHo0WEOsKi3E3CPVJtFC5RrgFlXtCtwC\nPBdlveNsE9kY4DoROT7aDlX1GVUdoqpDOnToUPc9hsTYVf32T4Q9njSumTKdoQ9+wZbdJcyePZsF\nCxZw2GGHWZFh/uiwSNEdqRmAy3qTm/iWQzMxGBohTs3bWxY6Xw1Gjx7N6tWrOenSOzh5RhEPLizj\nrQ0efF2ah3U9kSRaqEwEXre/vwJEfA1Q1V/szx3AG9HWSxg1sKtGrDtUBWXuVjy2pIysR3fyn//8\nhxtGZbH+pnZMKJkZeUTGSNEdpQWAWg77SG9Zl823wjENhsZAuObtD0TxhPtHqqZFixZMmTKFp05L\n46ONHs58qYihQ4fyRW4wMKawpJyXS4ez3pdpLQg7RrUtENFoDIFDNSTRQmUrcIL9fSTwffgKItJS\nRFr7vwOnAqvD10sYNYiciKnukANV5c033+SIPy7ktg+LOaZbKqtWreLxczI5IKk8NMTRSXh0B1Tq\n7K9zTESXId6kZliFH53XOFglVmoY7ntERzefTWrBrPEtycvL47YPS5m7upylH77C+rxCVpcdRHfZ\ngSqNO6S4nohnSPEcYAnQR0S2iMgVwP8Bj4rISmAKcKW97sEi8q69aSdgsb3OMmC+qr4fr35WSQ0i\nJ5ZuzK9QdygaK1euZNSoUZx11ln4VHn7opa8O6k9fVvvswSEt9T6jHZhO6M7ojj7KxAWLWMwNGhc\nScFrPNwcVhOztLgRES4akMJ3n8zks0ktOKdvEgMWXckgNnC2ezGplFnJ9mHHqIkForkRt4KSqjoh\nSlN2hHW3AmPt7xuBgfHqV52y9WsoL66gSQzr2Q6XgE9D6w751ea5Vw0nNzeX++67j+nTp3PggQfy\nxBNPcOyuWfT3rrN24hQQ/mGBqxJufnOYf7usEfDDJ3V2ugZDveNy4x/oMWDyXfz3mu0rpSUttueA\n2zYtq3K2fMZ57sWINYu4g0Um/RYIn8LFzy4NLf9iCGAy6m2i2kr9AiPcBpq7zDIxoRU0iezubTms\nc+uKdYeA0n17uO222+jVqxczZ87k5ptv5vvvv+e6667jB+3Kk57TyfH0CPWXiMuanz6u8uqpDmc/\nyS2CQsjnMcUhDU0Dl/0e7A9EqU10Vtn+EF9pUnIaLSjBjQcRazCweUX9ueBdS4pFskDUyseSt6p6\nfhVnkE4DxpS+j4Tz4Zv3beSHsVPtjqBJtE5LpnVackCg7Nu3jzXzp7Pho1n4ykqYOHEi9913H926\ndQOst6C7iibgwUWKx8ssXy+yO/WDnRssZ2WsN09qhiWE/DdfaYEl/PwhmU67tMHQWElKrblAcZZg\n8o8c6fPCxLf4Zt63jP91KT5Vyrzw15mf8VPGzaw4+mmG9ewe2MxvgVj4XdVpDE4LRa1oiOX4I2CE\nSiRWznbMaNi8jTMazK9JOPBfQKWlpTzzzDM8+OCD7Ek6kM5jrufvN13AOScMDll/6cZ8PLitisL2\nfHZqhnXz1CYsuKmUnDA0TerzzdvvV3QlWVPXoXyf4qVE0kinBJcohSfcSVnuJo4eN4HzTsrG3fE0\nSMvgntOOILt7WwpLyiko8ZCzebcxhdkY85dNqAMuPHw3Qjiv/w0HsUxOYW9NpaWlPPnkk/Tq1Ysb\nb7yRQ44+mYMvnYq79/HcvWB7BUffsJ7tSMKLGy/JhFUUrg1NqeSEoflRWmCFD0eLgKwJlYTUz71q\nOOmpVumj1a7eFHUdRsbwC8m89K+8u34P5ZKCx+Pl/rdXM/vLn2OK8ly7rYC12+qw/w0cI1SoGAK8\nvpPTzinWQFaRcCVZWopDkygpKeFf//oXhx56KNdeey1du3blww8/5OJb7kPtnztSRFh297ZMaTGH\nW5JeZVbLv9XdW8//t3fu4VHU5x7/vLu5CAEUDHIzF8MDCF5QAwHlwAl4qViLeJd6qtQLoFVRxPN4\nwWq91LanVlF7RCoochTBejSt0qoFDaAGmnggCkTkFgFpEEUhgCSb/Z0/ZjY72WxIspnNTsL7eZ48\nmZnMzH53srPv/N7fe4kMyczIazyLWVG8QKhihDMCsuep1mfa+fkNBmLKW2mMouDA2jkUxM85P78L\n8ScjPh8/VFXzh9eW1jYGayzK021cy5eJA2pUqD8Bt6QyO/zHbjlN8t0ePHiQp59+mr59+3LLLbeQ\nnZ3Ne++9x4oVKzj33HMZnpOOzx7w1OlE52Bg0g5uSvoruUlbXHhXDpwhmYrSVogWARlJKGAmlryV\nRgzRcN/6Ovfs+CEn4Pf7AINfhC+WLiJYXQUmSJJfGvQu1AQNhwLBIyYMWedUqB8CfHanreE/frvZ\n+uA14Pvdd8jwp3c28/v7+rJz505GjRrF/PnzGT16dJ0s+FBE2N4fAsy86vT4+l9T0sJPcrG6DzSp\nUUk0zhD50LxlZIi809CE5gwbe3gK1f0KGSKRcGCLg1zfF/Xu2f8p2lq7fuxNpzL6pkfZ8/1eanau\nZ97efLpPn14bfAOWF+RAlRUYMOFPRSy4sf2HIR+RI5XIoWNkCPCJFc4vVANrFliLVftrQwArKiq4\n793vyXzie+5cWMaJJ57I+++/T2FhIWPGjIlaVqXzUcn0OaZDgx+qk3odjT9aOZZYieY+UJS2QihE\n3p8anreM7DjqnCNs8pyhowpxtHpijiiryHvWuZ6dnU3e+OsYccHlXD56CM8++yx9+/Zl4sSJrF9v\n5Zv97yfbw6cNBOusNxtj6nZs9ShHhFG58rmPWbvz+8PuU/fDYyL+Gl7f8K9KJk+eTFZWFo8V7uPs\nE5JYOeMsli5dSn5+/mFfw1mau1VoivtAUbxMYxGQoYCZZuWtOO5vZz2xKDR2zy6cfCZvz7iMuXPn\nsnHjRm6++WYWLVrEoEGDuPjii6moqGjolZvHtlVAkFDH1j6Vn3o2s/+IMCrNZvBPHSsCg3/Kym2H\nuHRhJSfet4x58+YxceJEPr+jB3++shN5OcckTGpUep5qPdFFS6BUlPaGLym2vJWQIYri+qolYlRw\nOCOTlZXFzJkzKS8v5/7776ewsJD5D9+KCVSDMaT4hUvPOL55GkM4HghNTRW99hQ3ubZga6NGJRr2\nhzMQNLy7uxdnXjGV4c9+zdItAe79cV/Ky8uZNWsW/dKT3X9tY5oXyRJZ1NG5Hs19EO0YRUkUcYrc\nahINGSLjGE/EUFCye/fuPPTQQ5SXl/PYnTeyu+Ax9ix7CV/h05S88xoHDhxo9Bz1orscD4QBSeLj\nmoFA60edNYUj0qhEKwrnfAKpqKjg4cJDZD5RyY/+WMY333zD787vwpbbu/DIJQPo0aNH+GShrnRu\n3BSH9gJBaw7kxR+7E5vvRgKlosSDlkRutQj7ay8YaODvEXMuzXQbhwxC586dmTZtGtlXzqBzzyw6\nHqhg8uTJZGRkcO/rn7P9+4ZePwoZebZuYdPYV1hNf6DhSNJEckQYlVCPhJKDPQ5blr6kpIRrr72W\nzMxMfvnBIQb39LF48WLKysq4eXganVMjJtEju9K19Kao3BVerqmqu64o7Y1okVvxpnZugsO0h3Dc\n5y4kDPuSUug6+BxKSkpYtmwZ+fn5FJZu5dmPvuWXrxazcuXKesdEr4ZszamcOPScBmsLeoF2H1Ic\nMiJlZhQLt9Vw6Sfb6+SkLC/bSenSN5k9ezZFRUV06tSJyZMn84vgXAak+2Hs2IZPHq0Md7xzQZzu\nMc07UdoyMUVutZBIwxXNkKV2CnsJXOhR/+mDP6pdHjlyJCOzUwnOWYIJBjlUc5CzLxtBdY/TmDJl\nChMmTKBsd1W9asi/+3sZrxI2d5G1Bb1Eux+phBIbg/ioNn4M2AlNBhOs4aFbr+G6667ju+++48kn\nn2T79u089dRTlkFxsLomh1mBCynZvi+8MbIrXUtvik7HOc6XWncd6rrH3HQX6ByLkghiitxqIZH3\naIdGXEcxaGq058rW5fgw+H1Ch2ThD7ddzKFDh7jxxhvp3bs39z35QsIy9d2g3RuVUGKjjxqSpYaU\n7Z+w791n2FP4El8vmsGFwwexfPly1q1bx9SpUzn66KPrnaOkfA/XV93J44HLuXr/HZSU7+HK5z5m\n8+5K8CVDp57u3BTO/vIT37LWnWXrncUhA4c0RFhp+8QauRUrGXnUcW/9/e7DzK00nyZ1fXVEZYr4\nOPOquygtLWXFihWMGzeOD998kWB1FSZYgw/DGcd3dk1fa9Du3V9nZB5Dr8Au2LKcdSUr+dXWdXTp\ndQI5Iy9iyduz6NatW6PnKNr8DVX4CeKvrSDcr2od2YEtgIH9Ls19/PxteCzDWs7Iq1+2voNTa7Dx\npyxFORIIBixXdCwu4ZqqcJViF4jWc6Weiyojj83+bI4PbCOlx0DIyEOAESNGMGLECHbv3k3+7TOp\n+NdO/rX+Qy6cvYtjzrqS+074CZcevZ5cXCijH0farVHZsmUL8+fP56WXXmLTpk10SPZx+dBeTH75\nQ54sBRFpkkEBa7STQg3VUFtBeMMnpUgoSqSpnRmbS2TZ+mpHPwXxwcG2NSxWFNepbZZHuF9Qo/eh\nUBvh5U+xyrRE3a3hpMiGaKjrayQHfWlUSwopqV3q9VtJT0/npJEXsODr6VSc1oUnvz6H0vTRvOxP\n4pWDAS7/7X8zfeKldaNQPUS7cn/V1NQwZ84cRo0aRU5ODg888ACZmZnkXTuDD+4ZyrwbBnPWWWdF\nLaFyOHKzujIn5XGmJb1WW0F4bYqjVES8Egsjy9Z3dHxA/amazKgosUSQhe7/piQ/NpPDdX11sj7Q\nhz9Vj6XkYMOGwecTzjkpnbETp+JLTkV8fvD5mf3mB/Tp04dx48bxxhtvUFVV5Zp+N2hXRmXNmjXc\ncMMNVFRU8Oijj7J161aWLl1K9pkXkJYa/uDEUi7lNP9mbk76S20F4U+C/VgUGM2OYHrUfiquEFm2\nvnMvaudcWmtiU1G8TCwRZKHkRp/fuoe6ZkOHrq4FvjRW4y/U5fWpwEVcvW0c+36orrfPwslnclIv\n66EyNPoZImXcllxAwT3jmD59OsXFxVxyySX07t2bKVOmUFhYSDAYdOU9tIR2ZVTS09MpKiqirKyM\ne++9l6wsq/2n8x/kBqHJuHsC13F21X9REuzn2rnrEVm2XqR1JzYVxcs0N4IsMk+l+EWr2KpbuWZN\nINzl1U+18bP3h/qBAs56hblZXRmfvoMFKY9wR9Lr5JdO5ze/uIQvv/ySxYsXc9555zF//nzy8/PJ\nzMzkzjvvpLi4GGNirjTWItqVUcnMzGTYsGHNdm8dFvtDlpbir40ZCYcp+6kmiaIDfdx7vXgQqgWm\nKO2R5kSQRbrH1heEi666lIDZmCekTpdXqaHLUY2734ayliQJWp47W2dSUhJjx47llVdeYdeuXSxY\nsIDc3Fyefvpphg4dyoABA/jlGxtY/1Vli99Tc2hXRsU1nE8rUZ5eQsNRq/VvgOEdd7SyQEVRYiLS\nPTbwonDR1VZKwAx1eb01qYCXM/4C0GjF4TpzuFF0pqWlcdVVV1FQUEBFRQXPP/88GRkZPPLWRgbN\nWMYpp5zCAw88QGlpadxHMO3KqISa4bQY59NKKB8kFMseDNROxt2aVMDLKb8ht0NF9PO0lMjeEYqi\ntIzaGlpYbrMhE605UWe77dC9Hscv34FJO7gx+W8AUfNanKWlAL5IGRQ+uBE3X9euXbn++utZsmQJ\nOx4fw8wJg+jWrRsPP/wwgwcPpn///tx9992sWrUqLgamXRmVzbsr3SkDXSf/Iwg/7A2HLdr1gjof\nlcyNyX8j17+x5a8XTzRbXmnruO2+DbnHQ1FfqV3C85bOEGXi1xDrpF5Hk5bip+hAn3p5LaE52/mH\nRnH1tnGUlO9h4eQzicWp3+uYo7jt3GwKCwvZuXMnzz33HDk5OTz++OMMGzasdt7ZTdqVUTEGd0oa\n1Mn/8NX3s655peWvoShKbMTzQakptcHcIhhgvFlCrm8DEM5riSwtVbT5m0Zd8g3iuFY9evRg0qRJ\nvPPOO+zatYt58+Zx+umnu/2u2pdREcGdMtCRI5Wk1MhXavlrKIqSWKK5liPnVOI1x2JXy+izfy0v\np/yaH3Upr81riSwtNTznWNcrOnft2pVrrrmGgoKCFr6R+rSrjPqc9E7uVO2MHKl0HwDlH4U39RzM\nwgvPhMf84K28I0VpX7S26zYUolx9wPodr9B9R7WMZALkp35e+90VmrMd+u1fGdfjG3KzxoMvARWd\nY6RdjVQ6pjS/rEJUnP+wpFSrvbDfMVr5+92J6VR3aK8V/uhiATxFUSIIzbW4mGlfD0e1jADJdaO7\nsBIor0j9OBwE5DRuHk98bldGxTUa+we2VkMhJ9tWWUlamMM0F1IUpU3gqJbx8LGP1Y3uopGEbQ8b\nFFCj0jgV66zfPscoyOdv/eHn1uXhJK3QuqIo1txIWwy9t6tlPHLr9Z6uOtxc2tWcimsUvxhefmtq\nlB0SMFEf6sEQMiwe9qkqitIIP38bXvhxolXEBR2pRGN9Qf11Z+vgYKD1RwoZeVaSFhLfCURFUbyJ\n+GMqx9/aqFGJxsCL6q873V+Jir5I7WKNVuI5gagoRwLibziyLBgAxKperDQbNSrRGDIxvHzhTGu9\n9otcwpP3PU+FlDR3X7ut+ocVpT1Qm1FvrMAYDYhpNnEzKiIyV0R2ichnjm2DReRjEflURP4qIl0a\nOPZ8EflcRDaKyN3x0tgknAYGrJFCIl1PKWlqdBSlJRhjzU1GMxhOt3aoo6vSLOI5UnkROD9i2/PA\n3caYU4A3gLsiDxIRP/BHYCwwCJggIoMi91MURWk2tf1UTPRyJ063drw6urZz4mZUjDHLgG8jNvcH\nltnL7wGXRjk0D9hojNlsjKkCXgUuirKfoihK82is3Ekoox6JX0fXdo7Es7a+iGQDbxljTrbXPwJ+\nZ4x5U0SmAb8yxnSOOOYy4HxjzA32+s+AYcaYWxp4jUnAJHv1ZOCzaPt5iHRgd6JFNAHV6S6q011i\n0tk5hbR+3Xz9EQSD+eLb4IZ9VeyPg74QbeF6Doj8Hm4JrR1GdB3wlIjcD/wFFypnGWNmA7MBRKTY\nGDOkpeeMJ21BI6hOt1Gd7qI63UNEit08X6saFWNMGXAegIj0B6Jl/+wAMhzrx9vbFEVRFI/TqiHF\nInKc/dsHzABmRdntn0A/ETlBRFKAq7BGNYqiKIrHiWdI8QLgY2CAiGwXkeuxIrk2AGXAV8AL9r69\nRWQxgDEmANwCvAOsBxYZY9Y28WVnu/w24kFb0Aiq021Up7uoTvdwVWNcJ+oVRVGUIwvNqFcURVFc\nQ42KoiiK4hqeNioikiEi74vIOhFZKyJT7e2niUiRiKwWkWIRybO3H2vvXykiz0ScK9cuD7NRRJ4S\nEVfq18eg8VwRKbG1lIjImHhrjFFnnr1ttYisEZGLvajTcVym/X+f7kWdIpItIgcd13SW41ye0Wn/\n7VSxyimttXUd5TWdInK141quFpGgiJzmQZ3JIjLP1rNeRO5xnMsr30kpIvKCrWWNiOS3SKMxxrM/\nQC/gDHu5M7ABq3TLu8BYe/sFwAf2chrwb8AU4JmIc60ChmM1Q/lb6PgEaDwd6G0vnwzsiLfGGHV2\nBJIcx+5yrHtGp+O4PwOvAdM9ej2zgc8aOJeXdCYBpcBge/1YwO81nRHHngJs8uj1/CnwquOe2gpk\nx1NnDBp/AbxgLx8HlAC+WDV6eqRijNlpjPnEXt6HFQ3WBzBAqBjl0ViRZBhj9htjVgA/OM8jIr2A\nLsaYImNdqZeA8QnS+H/GmK/s7WuBDiKSGk+NMeo8YKxIPICj7P3iei1j0WlrGg9swbqeoW2e0xkN\nD+o8Dyg1xqyxj/nGGFPjQZ1OJmCVc/Li9TRAmogkAR2wEr73euk7CcvgLLX33wV8BwyJWaNbFjze\nP1hPel/aF2WgvbwNKzEyK2LfiThGKsAQ4B+O9ZFY5WMSptHe/7KQrtbS2BydwDCsL+pK4GIv6gQ6\nYYWudwIexB6peFBnNrAfWA0UAiM9qvN2YD5WSP8nwH96UWfE/puAk72oE0jGMnhf2///Sa2ps4ka\nJ2GN8pOAE7CMyqWxavT0SCWEiHQCXgduN8bsBW4C7jDGZAB3AHMSqQ+ar1FETgJ+C0z2qk5jzEpj\nzEnAUOCekG/dYzofBJ4wxlS2lrYYde4EMo0xpwHTgFekgdYPCdaZhOVCvtr+fbGInO1BnaH9hwEH\njDGtWvOvGTrzgBqgN9YX9p0ikuMxjXOB7UAx8CTwka05NuJhwV22tMlYT03THNu+J5xjI8DeiGMm\nUnek0gsoc6xPAJ5LlEas0jMbgBGtpTHWa+nYbynWk4undALLsfzUW7GesL7FSp71lM4ox37g0et5\nFTDPsd/9WC0qPKXT8fcngHsd657SidXG42eO/eYCV8RbZws/mx9hucRi0ujpkYodaTAHWG+M+YPj\nT18B/24vjwG+ONx5jDE7sfyYw+1zXgMUHO6YeGkUkWOAt7H6ynzYGhpj1HmC7QdGRLKAE4GtXtNp\njBlpjMk2xmRjPWX92hjzjNd0ikh3sXoFYT+p9gM2e00n1hfRKSLS0f7//zuwzoM6Q+WersCeTwHv\n3UdY7qYx9rFpWJPeZR77Tupoa0NEzgUCxpjY/+duWnC3f7CG3wYrGmW1/XOBvb0EWAOsBHIdx2zF\nelqtxBrSDbK3D8Eqi78JeAbbYre2RqyaZ/sd+64Gjounxhh1/gxrPmU1lm99vONcntEZceyD1I3+\n8oxOLB+183r+xIs67WP+w9b6GVarCq/qzAeKopzLMzqx5vpes6/nOuCueOuMQWM28DnWhP4/qDuv\n2myNWqZFURRFcQ1Pu78URVGUtoUaFUVRFMU11KgoiqIorqFGRVEURXENNSqKoiiKa7Rqj3pFaS+I\nyLHAEnu1J1YG8tf2+gFjzFkJEaYoCUZDihWlhYjIg0ClMeb3idaiKIlG3V+K4jIiUmn/zheRQhEp\nEJHNIvIbsfqArLJ7VPS19+suIq+LyD/tnxGJfQeKEjtqVBQlvgzG6u8zEKtKQX9jTB7wPHCrvc9M\nrKKYQ7Ey759PhFBFcQOdU1GU+PJPY9VQQkQ2YTVKAvgUGG0vnwMMcjTV6yIinUyCKi8rSktQo6Io\n8eWQYznoWA8Svv98wHBjTJ3mcorSFlH3l6IknncJu8IQu9e6orRF1KgoSuK5Dat9a6mIrMOag1GU\nNomGFCuKoiiuoSMVRVEUxTXUqCiKoiiuoUZFURRFcQ01KoqiKIprqFFRFEVRXEONiqIoiuIaalQU\nRVEU1/h/VCrW6z+gUAUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Now suppose you wanted to combine the two together:\n", "my_event = MulensModel.Event(\n", " datasets=[OGLE_data, MOA_data], model=my_1S2L_model)\n", "\n", "# And you wanted to plot the result:\n", "my_event.plot_model(\n", " t_range=[2452810,2452890], subtract_2450000=True, color='black')\n", "my_event.plot_data(subtract_2450000=True, label_list=['OGLE','MOA'])\n", "# MulensModel automatically fits for the source and blend flux for the \n", "# given model.\n", "\n", "# Customize the output\n", "pl.legend(loc='best')\n", "pl.title('OGLE-2003-BLG-235/MOA-2003-BLG-53')\n", "pl.ylim(19., 16.5)\n", "pl.xlim(2810,2890)\n", "pl.show()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Chi2 of the fit: 2194.16\n" ] } ], "source": [ "# If you want to see how good the fit is, output the chi2:\n", "print('Chi2 of the fit: {0:8.2f}'.format(my_event.get_chi2()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you want to optimize the chi2, we leave it up to you to determine the best method for doing this." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "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.6.2" } }, "nbformat": 4, "nbformat_minor": 2 }