{ "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 planet. See [Bond et al. 2004](http://adsabs.harvard.edu/abs/2004ApJ...606L.155B). The data were downloaded from the [NASA Exoplanet Archive](https://exoplanetarchive.ipac.caltech.edu/cgi-bin/DisplayOverview/nph-DisplayOverview?objname=OGLE-2003-BLG-235L+b&type=CONFIRMED_PLANET)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import basic packages\n", "import MulensModel\n", "import matplotlib.pyplot as plt # MulensModel uses matplotlib for plotting.\n", "import os.path" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Define a point lens model:\n", "my_pspl_model = MulensModel.Model(\n", " {'t_0': 2452848.06, 'u_0': 0.133, 't_E': 61.5})\n", "\n", "# Or a model with 2-bodies:\n", "my_1S2L_model = MulensModel.Model(\n", " {'t_0': 2452848.06, 'u_0': 0.133, 't_E': 61.5, 'rho': 0.00096, \n", " 'q': 0.0039, 's': 1.120, 'alpha': 223.8})\n", "# Since rho is set, define a time range and method to apply \n", "# finite source effects:\n", "my_1S2L_model.set_magnification_methods(\n", " [2452833., 'VBBL', 2452845.])" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1kUlEQVR4nO3dd3wU1frH8c9DEogIhBYNLYRiaKEoIBBaABEEFEUpFgQE8UfzCgIigoqIgnK9XCwUAZUiKFZQkCLVK0jvEYlIVYTQIZSU8/tjJzGEQHaTbGY3+7xfr32xMzs782SzfHP27JkzYoxBKaWU78hjdwFKKaVylga/Ukr5GA1+pZTyMRr8SinlYzT4lVLKx/jbXYAzihcvbsLCwuwuQymlvMrmzZtjjTHBadd7RfCHhYWxadMmu8tQSimvIiIH01uvXT1KKeVjNPiVUsrHaPArpZSP8Yo+fqWUckV8fDxHjhzh8uXLdpeSIwIDAyldujQBAQFOba/Br5TKdY4cOULBggUJCwtDROwux62MMZw8eZIjR45Qrlw5p56jXT1KqVzn8uXLFCtWLNeHPoCIUKxYMZc+3WjwK6VyJV8I/WSu/qwa/Mrrbdu2jV9++cXuMpTyGhr8yuvdeeed1K9f3+4ylLqGn58ftWrVIiIigo4dOxIXFwfAmDFjqFatGjVq1KBWrVopjZaoqCgqVapEzZo1adiwIXv37k1Zn90nsGrwK6WUG9xyyy1s27aNXbt2kTdvXiZPnsy6dev47rvv2LJlCzt27GD58uWUKVMm5Tlz5sxh+/btdOvWjSFDhritNg1+pZRys8aNGxMTE8Nff/1F8eLFyZcvHwDFixenZMmS123fpEkTYmJi3FaPBr9SKveLioKPP3bcj493LM+e7ViOi3Msf/aZY/nsWcfyV185lmNjHcsLFzqWjx1z6dAJCQksXryY6tWrc++993L48GHCw8Pp27cvq1evTvc5CxcupHr16i4dxxVuC34RmSEix0VkV6p1b4vIryKyQ0S+FpHC7jq+UkrZ6dKlS9SqVYs6deoQGhpKz549KVCgAJs3b2bq1KkEBwfTuXNnPk7+gwQ8/vjj1KpVi//973+MHz/ebbW58wSuj4H3gJmp1i0DXjTGJIjIOOBF4AU31qCUUrBq1T/3AwKuXc6f/9rloKBrl4sXv3Y5JMSpQyb38afl5+dHVFQUUVFRVK9enU8++YTu3bsDjj7+OnXqOLX/rHBbi98YswY4lWbdUmNMgrW4HijtruMrpZSn2bt3L/v27UtZ3rZtG2XLls3xOuycsuEp4DMbj6+UUjnqwoULDBgwgDNnzuDv70/FihWZOnVqhs9r27Ztyjw8DRo0YP78+Vmqw5bgF5GXgARgzk226Q30BggNDc2hypRSKntcuHDhunW1a9fm559/Tnf7Vam7k5xYnxU5PqpHRLoD7YDHjTHmRtsZY6YaY+oYY+oEB1935TCllFKZlKMtfhFpDQwFmhpj4nLy2EoppRzcOZxzLrAOqCQiR0SkJ45RPgWBZSKyTUQmu+v4Siml0ue2Fr8x5tF0Vk931/GUUko5R8/cVUopH6PBr5RSPkaDXyml3OCpp57itttuIyIiImXd+vXrqVevHrVq1aJKlSq8+uqrgOOM3Ro1alC9enUiIyPZvn17ynMKFCiQ7bVp8CullBt0796dH3744Zp13bp1Y+rUqSnTNXfq1AmAcuXKsXr1anbu3MnIkSPp3bu3W2vTi60rpZQbNGnShAMHDlyz7vjx45QoUQJwzNlTtWpVACIjI1O2qV+/PkeOHHFrbRr8Sqlc7bnnnkt3srSsqFWrFhMmTHD5eQMHDqRSpUpERUXRunVrunXrRmBg4DXbTJ8+nfvuuy+bKk2fdvUopVQOefnll9m0aRP33nsvn376Ka1bt77m8ZUrVzJ9+nTGjRvn1jq0xa+UytUy0zJ3pwoVKtCnTx+efvppgoODOXnyJMWKFWPHjh306tWLxYsXU6xYMbfWoC1+pZTKId9//z3JU5Tt27cPPz8/ChcuzKFDh+jQoQOzZs0iPDzc7XVoi18ppdzg0UcfZdWqVcTGxlK6dGlGjRrFsmXLGDhwIPnz58ff3585c+bg5+fHa6+9xsmTJ+nbty8A/v7+bNq0CYC4uDhKl/7n0iWDBg1i0KBBWapNbjJBpseoU6eOSX4RlEpLRADwhveyyhnR0dFUqVLF7jJyVHo/s4hsNsZcd0kv7epRSikfo8GvlFI+RoNfKZUr+VLXn6s/qwa/UirXCQwM5OTJkz4R/sYYTp48ed2JYDejo3qUUrlO6dKlOXLkCCdOnLC7lBwRGBh4zcifjGjwK6VynYCAAMqVK2d3GR5Lu3qUUsrHaPArpZSP0eBXSikfo8GvlFI+RoNfKaV8jAa/Ukr5GA1+pZTyMRr8SinlYzT4lVfzhVPylcpuGvzKqyUlJdldglJeR4NfeTUNfqVc57bgF5EZInJcRHalWldURJaJyD7r3yLuOr7yDYmJiXaXoJTXcWeL/2OgdZp1w4AfjTF3AD9ay0plmga/Uq5zW/AbY9YAp9Ksbg98Yt3/BHjQXcdXviEhIcHuEpTyOjndx3+7MeYv6/4x4PYcPr7KZbTFr5TrbPty1zjG4d1wLJ6I9BaRTSKyyVcupqBcp8GvlOtyOvj/FpESANa/x2+0oTFmqjGmjjGmTnBwcI4VqLyLBr9Srsvp4F8AdLPudwO+zeHjq1xGg18p17lzOOdcYB1QSUSOiEhPYCzQUkT2AfdYy0plmga/Uq5z2zV3jTGP3uChFu46pvI9GvxKuU7P3FVeTYNfKddp8CuvpsGvlOs0+JVX0xO4lHKdBr/yatriV8p1GvzKq2nwK+U6DX7l1TT4lXKdBr/yahr8SrlOg195NQ1+pVynwa+8mga/Uq7T4FdeTYNfKddp8CuvpsGvlOs0+JVX0+BXynUa/Mqr6Zm7SrlOg195NW3xK+U6DX7l1TT4lXKdBr/yahr8SrlOg195tfj4eLtLUMrraPArt4mNjeW3335z6zEuX77s1v0rlRu57dKLSoWHh3P69GmMMW47xpUrV9y2b6VyK6eCX0TCgSFA2dTPMcY0d1NdKhc4ffq024+hwa+U65xt8c8HJgMfAvptmvIYGvxKuc7Z4E8wxkxyayVKZUJy8OfJo19XKeUsZ/+3LBSRviJSQkSKJt/cWplSTkgOfn9//bpKKWc5+7+lm/XvkFTrDFA+e8tRyjXJwe/OL5CVym2cCn5jTDl3F6JUZmgfv1Kuc3ZUTwDQB2hirVoFTDHG6Nkzylba4lfKdc529UwCAoAPrOWu1rpe7ihKKWdp8CvlOmeDv64xpmaq5RUist0dBSnliri4OLtLUMrrODuqJ1FEKiQviEh5sjCeX0QGishuEdklInNFJDCz+1K+7cKFC4C2+JVyhbMt/iHAShHZDwiOM3h7ZOaAIlIKeBaoaoy5JCKfA12AjzOzP+X5kpKS3DbO/uLFi4AGv1KucHZUz48icgdQyVq11xiTleEU/sAtIhIP5Af+zMK+lIdLTEx0W/Ant/iVUs67afCLSHNjzAoR6ZDmoYoigjHmK1cPaIw5KiLjgUPAJWCpMWZpOsfuDfQGCA0NdfUwyoMkJiYSEBDgln1ri18p12XUDGtq/Xt/Ord2mTmgiBQB2gPlgJLArSLyRNrtjDFTjTF1jDF1goODM3Mo5SHcebEUbfEr5bqbtviNMa9Yd18zxvyR+jERyexJXfcAfxhjTlj7+QqIBGZncn/KwyUlJblt3xr8SrnO2Y7XL9NZ90Umj3kIqC8i+UVEgBZAdCb3pbyAu1r8SUlJnD17NmVZu3uUck5GffyVgWpAUJp+/kJApoZgGmN+EZEvgC1AArAVmJqZfSnv4K7gP3/+PMYYChUqxLlz5zDG4GhLKKVuJqNRPZVw9OUXxtGvn+w88HRmD2p1Ib2S4YYqV3BX8J85cwaAIkWKcO7cObccQ6ncKKM+/m+Bb0WkgTFmXQ7VpHIZdwV/8hW+ihQpwsGDB7WrRyknOXsC11YR6Yej2yeli8cY85RbqlK5iruC/9SpUwAULeq4NES2B//OnfD333DPPdm7X6Vs5uyXu7OAEKAVsBoojaO7R6kMuWtUz7FjxwAoUaJE9u308uV/7o8bB717/7M8eDCMHJl9x1LKJs4Gf0VjzEjgojHmE6AtUM99ZancxF0t/r/++gv4J/iz3OL/8ksoUwb273csjxoFCxf+8/jFi2D9sVHKmznb1ZM87/4ZEYkAjgG3uackldu4K/iPHTtGYGAghQsXBrIh+OvXh9atwc/PsVyhwrWPT5oEyZ9e9u2DiRNh/HjIly9rx1Uqhznb4p9qnXE7AlgA7AHecltVKldxZ/CHhIRkbQjnuXPw7rtgDJQqBbNmQdmyN94+ec6hZctg/nw4eDDzx1bKJs5O0jbNursGvc6ucpE7u3pCQkJSljPV4v/kExg0CJo0gZo1r3vYGMPRo0c5ePAgfn5+VKhQgeDgYOjbF7p2hYIFHRsmJf3zR0EpD+fUO1VE3hCRwqmWi4jI626rSuUq7vpyNzn4k1v8mQr+/v1h06brQv/8+fOMGTOGihUrUqZMGRo1akSDBg247bbbuOuuu5g0aRJX8uZ1bDxxInTp8k83kFIeztkmyn3GmDPJC8aY00Abt1Skch13tPiNMRw4cICwsLDMdfW8/TYcPQoi14X+woULqVixIiNGjKBcuXJMnDiRxYsX89133zF27Fj8/Pzo27cvlSpVYvny5XDliiP0r17Npp9OKfdyNvj9RCTlGywRuQXQb7SUUxISErJ9n8eOHSMuLo6KFSumrHO6xX/wILz2GsyYcc1qYwwvvvgiDzzwAKVKlWLDhg0sX76cAQMG0Lp1a9q2bcsLL7zAhg0bWLp0KYGBgbRs2ZJhsbEkzZsHgXohOeUdnA3+OcCPItJTRHoCy4BP3FeWyg2SW+LuCP6YmBgAKlSo4HpXT9mysH07vPhiyqqkpCT69OnD2LFjefrpp/n555+pW7duuk8XEVq2bMnWrVvp3bs34956i86PPsrVQ4egY0f9wld5PGe/3B0nIjtwzKQJMNoYs8R9ZancwN/fn/j4eOLj4zPe2EW///47ABUrVmTHjh3OPenSJVizBlq1gvLXjlEYOnQoU6ZMYdiwYbzxxhtOdR/dcsstTJ48mfDwcAYPHkzC2bN8vn07ATt23HxkkFI2c3YcP8aYxcBiN9aichl3Bv+ePXvImzcvZVMFbIYt/gkT4KWXYM8eqFw5ZfWHH37Iv//9b/r16+d06CcTEZ5//nny5cvHgAEDeLJjR+a0bev0R2ml7JDRtMw/GWMaich5IPX/KgGMMaaQW6tTXs3POhHKHV09W7ZsoXr16gQEBDjf1TNoEEREXBP669ato2/fvrRu3ZoJEyZk+pyA/v37ExcXxwsvvEDFSpUY3agRlCgBNWpkan9KuVNGLf4nAYwxBXOgFpXL+Ps73l7Z3eI3xrB161YefvhhgIzD+vJlxxj7fPng/n9mFz937hyPP/44pUuXZt68eSn1ZtaQIUPYt28fr7/+OlWCgnisdWuYNy9L+1TKHTL6RDofQER+zIFaVC7jruA/dOgQp06d4s4777xm/Q1b/CNHQt26EBd3zeoBAwZw8OBBZs+eTVBQUJbrEhHef/99mjZtylOXLrFt4MAs71Mpd8ioiZNHRIYD4SIyKO2Dxph33FOWyg2Sgz+7u3o2b94MwF133QWQcVdP06aOoZb586es+uGHH5g5cyYjR46kYcOG2VZb3rx5mT9/PrVq1aLzk0+yef16CsTHw206tZXyHBm1+LsAiTj+QBRM56bUDbmrxb927VoCAwOpaZ14lWFXT7t2MHp0yuKlS5fo168flSpV4qWXXsrW2gCCg4P59NNPiYmJoV+VKtCpk2MuIKU8REZX4NoLjBORHdaoHqWclvzlbnYH/4oVK2jYsCGBaU6Yuq7F/9FHjrNqe/e+Zh6dN954g/3797NixQryuWlmzaZNm/Lyyy/z6quv0rJLF55wy1GUypybtvhFJPn9WlVEBqW95UB9you5o6vnxIkT7NixgxYtWqSsu2FXz8KFjjn2U30iiImJYdy4cXTt2pVmzZplW13pGTFiBI0aNaL/xx9z5OhRtx5LKVdk1NVzq/VvAbSrR7nIHV09S5cuBaB58+Yp627Y1fPll9cF/4gRIwgICGDcuHHZVtON+Pn58fHHHxMfH8/T992HSXWmsFJ2yqirZ4r176icKUflJu7o6vn6668pUaJEutMppLT4f/0VQkKgcGEo9M+pJhs3buSzzz5j5MiR2Xu5xpuoUKECb731Fv3792d6QgK9Ro+GLA4bVSqrnHoHikgw8DQQlvo5erF1dTPZ3dVz6dIlFi9eTPfu3cmTqs/+mq4eYxxTJOfLB+vXp7T2jTG88MILBAcHM3jw4Gypx1l9+vThqy++YNDmzbQ8evSas42VsoOzZ5Z/CwQBy4HvU92UuqHs7upZsGABcXFxKSduJbumq0fEMevmuHHXdPEsWbKElStXMnLkSAoVytkTzvPkycP0jz7CGEPvp57CrF2bo8dXKi1nP3PmN8a84NZKVK6T3cE/Y8YMQkNDiYqKSvdxk3whFGt8f8p6Y3j55ZcJCwvjmWeeyZZaXBUWFsabb77JgAED+HTDBh4/fhxuucWWWpRytsX/nYjohVeUS7Iz+A8dOsSyZcvo0aPHNd08kKqrp39/eOWV6567bNkyNm7cyPDhw8mbfNUsG/Tp04f6NWvynL8/sRcv2laHUs4G/79whP8lETknIudF5Jw7C1PeLzmQL1++nOV9ffjhhwB069bthsfBz++a7p1kr7/+OqVLl+bJJ5/Mch1Z4efnx9RZszhz4YLjewY3XYtYqYw4Ox+/Dt1ULkseZXPhwoUs7efixYt88MEHtG/fnnLlyt34eBMmQPHi16xbvXo1a9eu5d1333XbyVquqF69OkOHDuWNN96g64YNtNi50/EHS6kc5OzF1u9K51ZBRDI1Lk1ECovIFyLyq4hEi0iDzOxHebbk4L+YxW6NGTNmcOrUKYYMGZLu47J79z/HS9Pif/3117n99tvp2bNnlmrITiNHjuSOkBCeOXKES6dO2V2O8kHOdvV8AKwHPrRu63HM3LlXRO7NxHH/C/xgjKkM1ASiM7EP5eGyI/ivXr3KO++8Q2RkJJGRkddvEB+P3GDq4y1btrB8+XKef/55bvGgL1IDAwOZ8umn/H7+PK+9o/McqpznbPD/CdxpjKltjKkN1AL2Ay2Bt1w5oIgEAU2A6QDGmKvGmDOu7EN5h+zo6pk2bRoHDhxgxIgR6W8QEIAMHQo4rpub2n/+8x8KFChA7969M318d2nWrBk9evTg7bffZrue0atymLPBH26M2Z28YIzZA1Q2xuzPxDHLASeAj0Rkq4hME5Fb024kIr1FZJOIbDpx4kQmDqM8RWZb/BcuXOC1116jadOmtG7dOr0dA5DXmvI49eiho0ePMm/ePHr27Jktc+27w/jx4ymWLx+93nqLxL/+srsc5UOcDf7dIjJJRJpatw+APSKSD3B1rJ4/cBcwyRhzJ3ARGJZ2I2PMVGNMHWNMneDgYBcPoTxBcov//PnzmXr++PHj+fvvv3nzzTevn48nMRGiomDAAAICAgBHt1Cy999/n8TERJ599tlMHTsnFC1alP++/z6bkpJ497PP7C5H+RBng787EAM8Z932W+viAVenODwCHDHG/GItf4HjD4HKZZKDPzOf2GJiYhg7diydO3emQYN0vvtPSoL27aFx45Sx+cnBf/HiRSZPnsxDDz1E+fLlM/8D5IDO3brRpk0bXnrpJQ5s3Gh3OcpHOBX8xphLxph/G2Mesm7jjTFxxpgkY4xLHbjGmGPAYRGpZK1qAexxsW7lBZKD/9ixYxlfCD3N8/r370/evHl550ZffgYEwIgR0KnTdcE/c+ZMTp8+zUAvuPShiDBp0iQkIYH/i4zE/P233SUpH+DscM47rOGXe0Rkf/ItC8cdAMwRkR04vih+Iwv7Uh4qOeyvXr3K6dOnnX7e3LlzWbJkCa+//jolS5a8foPx4+Gnn1IWUwe/MYaJEydSp06dbL2kojuFhoby5pAhLElI4NNFi+wuR/kAZ7t6PgImAQk4unZmArMze1BjzDar/76GMeZBY4zzqaC8RupW/rFjx5x6zqFDh+jbty+RkZH07dv3+g0uXoT334f581NWpQ7+1atX8+uvv9KvX7+ML8noQfqOGkW9evV4buhQYmNj7S5H5XLOBv8txpgfATHGHDTGvAq0dV9ZKjcwxlCgQAEA/nJi1EpiYiLdu3cnMTGRWbNmpcz1c41bb4WdO2HMmJRVqYN/0qRJFClShM6dO2fPD5FD/Pz8mDZtGmdOn2bQ3XdDNl+gXqnUnA3+KyKSB9gnIv1F5CEcV+VS6oaMMSkXPHEm+F955RVWrlzJxIkT0/9Sdvt2x5e6BQo4bpbk4D98+DBfffUV3bt396gTtpwVERHBiw8+yKw//mDJrFl2l6NyMVcmacsPPAvUBroC18+WpVQaZcuWJU+ePOzdu/em23355ZeMGTOGXr160b179+s3+OsviIyEdE52Sg7+yZMnk5CQYNvUy9lh+KxZVAoP5/9eey3LU10odSPOjurZaIy5YIw5YozpYYzpYIxZ7+7ilHczxpA/f36qVavGxpsMVfzf//7Hk08+Sf369XnvvffS75sPCYEpUyCdfv/k4F+3bh3NmzenUqVK123jLQJvuYUPrbOVX+7Rw+5yVC5100nWRGTBzR43xjyQveWo3MQYg4hQt25dvvnmG5KSkq6bS3/btm20bduW0qVL880336Q/g2by5GtPPJHucVLPsd+nT59s/Rns0LhxY/6vUSMmzJ9Pl5kzqWvzdNIq98moxd8AKA2sBcYD/05zU+qGkoO/RYsWnDp1ilWrVl3z+MqVK2nWrBmFChVi2bJl3H777dfv5MwZqFMHFi++4XGSgz8kJIT27dtn409gn7Hz51OiaFG6jR3LpUuX7C5H5TIZBX8IMByIwDGjZksg1hiz2hiz2t3FKe+WHPwPPfQQRYoUYfTo0cTHx3PlyhXGjBlDq1atKFmyJGvWrCE0NDT9nZw+7ThZK808+6kljxzq1atXyvQN3i4oJIQZc+cSHR3NSJ3ETWU3Y4xTNyAfjmkaTgD9nX1edtxq165tlPepVq2a6dChgzHGmGnTphnAlClTxhQuXNgApmPHjub06dMZ7ygpKcNNVq5caS5fvpzFij1Pn0ceMQJm9Xvv2V2K8kLAJpNOpmb45a6I5BORDjhO2OoHTAS+dtPfIZWLGKvFD9CzZ0+++OIL6tatS4cOHVixYgWff/45hQsXTv/JZ8/CqFFw6VK6l1NMKyoqyiOusJXd3po4kXL589P9zTczPdmdUmll9OXuTBzdPIuAUcaYXTlSlcoVTJr5eR5++GEefvhh5578/fcwejS0bevo4/dRBUqU4JMlS2jSpAlDhgxh8uTJdpekcoGMWvxPAHfgGMf/s3Whdb3YunJapqdNeOwx2LfPp0M/WaNGjXj+ueeYMmUKP+gVu1Q2uGnwG2PyGGMKWrdCqW4FjTGFcqpI5Z1Sd/U4LTYWfv3Vcf8mF1b3NaOHDaNaQAA9XnmFv3UGT5VFzp65q5TLMhX8gwdDgwZwTj9QphZ4223MXbWKMwkJdOvW7brLTCrlCg1+5TaZCv433oDp06GQfqBMq3pkJP/5z39YsmQJ7+SCE9WUfTT4ldu4FPxnzzrO0C1ZEjp0cG9hXuyZZ57h4eBgXpw6lQ0//2x3OcpLafArt3E6+K9ehWbNoH9/9xfl5USED5cupWSpUnR54gnOnj1rd0nKC2nwK7dJO5zzhvz94ZFHoGVL9xaUSxSpVYu5n3/OoUOH6Pnwwy5d1lIp0OBXbpZhiz8pCfLkgeHD4cEHc6Sm3CAyMpI3H3mEL3/8kfH9+tldjvIyGvzKbTLs6tmzB2rWdFxRS7ls8IwZdKxVi2FTprB8+XK7y1FeRINfuU2GwX/pEtxyCxQrlnNF5SKSPz8z1q6lcuXKdOncmYO//WZ3ScpLaPArt8kw+GvXhl9+cYzkUZlSoEABvp4zh/gzZ+jQsKFO4aycosGv3OaGwT96NIwb988FVlSWhNeqxeyuXdl68iQ9evTQk7tUhm46SZtSWZFu8CclOfr2c+FMmna6/+OPGVetGkOHDqVCmTKMefttu0tSHkyDX7lNusGfJw98+ikkJGhrP5sNHjyYfUuW8Mb48VQsUoQew4fbXZLyUNrVo9zmmvHlhw9Dx45w4oQj8HPJlbI8iYjw/pQptCxRgt6vvMLKlSvtLkl5KA1+5VYpLf7oaFizBo4ft7egXC6gQgXmR0cTHh7OQw89xPb16+0uSXkgDX7lNtd09dx7L+zfD9Wq2VuUDwgKCmLxokUUvHKFVk2bEqPDPFUaGvzKbYwxyMqVsHixY8Wtt9pbkA8JLVuWZSNHkpg3L/fcey9Hjx61uyTlQWwLfhHxE5GtIvKdXTUo9zKJici5c47+fZXjKg8fzg+rVnHq1ClaNmtG7IkTdpekPISdLf5/AdE2Hl+5S1ISJCVhRJD27aF3b7sr8lm1a9dmwbRp7N+3j3tr1uTkyZN2l6Q8gC3BLyKlgbbANDuOr9zIGOjbF3r1cnT1+OuIYbtFdezI148/zp6TJ2nRogWxsbF2l6RsZleLfwIwFLjhKYYi0ltENonIphP6EdV7iEBICISE6HTBnkKE+2bP5tsFC9i7dy/N69XjhI6u8mk5Hvwi0g44bozZfLPtjDFTjTF1jDF1goODc6g6lWmJifDnn477r74KY8Zk7tKLym1atWrFwhEj2Ld/P83r1uXYsWN2l6RsYkeLvyHwgIgcAOYBzUVktg11qOz07LOOi6QnXxHKCnwNfs9yz4sv8v2QIeyPjaVhw4bExMTYXZKyQY4HvzHmRWNMaWNMGNAFWGGMeSKn61DZrGdPGDgQgoIwxjBx4kROnjxJUFCQ3ZWp1PLkoflbb7FixQrOnjlDw5o12aInefkcHcevMi82Fj77zHH/rrvgueeIjY3lgQce4F//+hdt2rRh6NCh9tao0lWvXj1+euklAuPiaNqiBT/++KPdJakcZGvwG2NWGWPa2VmDyoI33oAePeCvvwBYuXIlNWrUYOnSpUycOJEFCxZQtGhRm4tUN1J50CB+3rCBsPLlad26NdM++MDuklQO0Ra/cl3yfO+vvw5r1xJfvDgjRoygRYsWFCpUiF9++YUBAwZo/74XKFW3LmvXrqV5zZo83a8fA594gsTERLvLUm6mwa9cM2UKtGgBV65A/vwcKFaMpk2bMmbMGHr06MHmzZupVauW3VUqFxQuXJjvZ87k2fLlmTBnDu3ateNs8pf0KlfS4FeuKVwYChWC+Hjmzp1LzZo12b17N3PnzmX69OncqvPxeCX/qlX57++/M8W6cHu98HB2bd1qd1nKTTT4Vcb274flyx33O3fm/KxZdOvXj8cee4yIiAi2b99Oly5d7K1RZYvevXuz/LXXOHP8OHfXr8/MmTPtLkm5gQa/ylifPvD00xAfz4YNG7jzrruYPXs2L7/8MqtXryYsLMzuClU2avrii2xbtIh6kZF069aNXk8+qRdxz2U0+FX6Ll2Cixcd9z/8kCuLF/Py6NE0bNiQq1evsmrVKkaNGoW/zsWTK4Xcdx/Lli3jpX79mD5rFndXrMj27dvtLktlEw1+db3Ll6FuXXj+eQA2/v03tR95hNGjR9OlSxe2b99O48aNbS5SuZu/vz+vjx3LD/ffT2x8PHXr1mXcuHE66icX0OBX1wsMhK5duXz//QwbNoz69etz+vRpFi5cyKxZsyhSpIjdFaqcUqAArRYsYOeePbRv355hw4YRVakS+/fvt7sylQUa/Mrh99+haVPYuROApXfeSY2BAxk3bhw9evRg9+7dtGun59r5quLFi/P5nDnMqlGDnUeOUKNGDSZMmKCtfy+lwa8cgoLg5EkObd3KI488QqtWrTDGsHTpUqZNm0bhwoXtrlDZTPLm5Ylt29i5Zw9NmzZl4MCB1AsNZcsvv9hdmnKRBr8v+/FHx4gdY7hSsCBjH3+cKn36sGjRIsaMGcOuXbto2bKl3VUqTyJCmfLl+e677/jskUc4euwYdSMjef755zl//rzd1SknafD7sl27SFq+nE+nTKFy5cq8OHw4rVq1Ijo6muHDh5MvXz67K1QeSkToNH8+0dHRPP3007zzzjuElyzJR2+/TVLSDa+vpDyEBr8viYuDIUNg8WIAVlSpQt2CBXm8Tx8KFy7M0qVL+eqrryhbtqzNhSpvUTg8nMmTJ/PLZ59RLi6Op4YOpa41/4/yXBr8vsTfH77/nk1ff02bNm1o0aoVJ2JjmTlzJps3b9ZuHZVpd3fqxP/++INPZ8/m+PHjNGnShEdq1ybaGiygPIsGf263YAG0bOk463bbNtqGhlL3ww9Zv349b731Fr/99htdu3YlTx59K6iskdBQHn38cfbu3cuomjVZsm0bEbVq0a1bNx3+6WH0f3tulJQE8fEpi+uPHuW+e++lXr16rN+4kTFjxnDgwAGGDBlCYGCgjYWq3Ch//vy8vHUr+3fuZNCgQXz++edUqliR/2vfnsOHD9tdngLEGGN3DRmqU6eO2bRpk91leIdTp6BBA5L69OG78uV55513WL16NcWKFWPw4MH069ePggUL2l2l8iF/rl3LG+3aMTUuDkR44tFHGTJ4MFWqV7e7tFxPRDYbY+qkXa8t/tzg3Dn46ScA4gIDmXT77VR++23at2/P/v37GT9+PAcOHGDYsGEa+irHlWzcmPdOnmTfvn0888wzzJs7l6o1avBQu3as1+v92kKDPzfo04cDbdvy0tChhIaG0nftWoJKlmTu3Ln8/vvvPP/88xQoUMDuKpUv8/enbFgY7777Lgdnz+bl+vVZ/fPPNGjQgCYREcyfNIn4VN2Tyr00+L3Rnj3w4IMkHD7Mt99+S5vDhyl//jxj//1vGjVqxJo1a9iwYQNdunQhICDA7mqVukZwp06MWreOQ4cO8Z+33+ZIdDSd+vYlLCyMUaNG8Zd1DWflPtrH7y2OH4erV6F0aQ6tWcOMNm2YdsstHI2NpWTJkvTs2ZNevXoRGhpqd6VKuSRx/34Wr1jB+19+yQ8//IA/8HCLFvQePpyoqCgdcZYFN+rj1+D3BleucOa22/jyzjuZBaxevRoRoVWrVjzzzDO0a9dO58VXuULMl18yqX9/ZsTFcebcOUJLlKBbs2Z0e+UVKoSH212e19Hg9zbjxnF161Z+eOwxZs+ezYJvvuFKfDzh4eF07dqVJ554Qq98pXKtS5cu8e233/LR0KEsO3wYAzRp0oTunTrR4fHHCdJJA52iwe/p9u2DL77g8nPPsXTZMr589VUW7NrFmfh4goOD6dKlC127dqVOnTqIiN3VKpUzrlzhyPLlzNqxg48++oh9+/aRN08eWrVtS6dOnXjggQcoVKiQ3VV6LA1+T2OMY+778uW5ACx+4QW+/OADvs+fnwtxcRQuXJgHHniAjh070qpVK/2SVvk8YwwbXnqJz7du5fNduzhy5Aj58uThvjp16PTcc7Rp04agoCC7y/QoGvye4MoVx2UNg4LY/803LHroIRbVqMHK337j8uXLBBcrxoMdOvDwww/TrFkz8ubNa3fFSnmkpKQk1i9axOdPP838K1f48/Rp/P39aVKiBPd37cr9Tz1FhQoV7C7Tdhr8domPh4AArpw+zdoyZVgUEcGiM2fYu3cvAOEVKnBfu3Y8+OCDNG7cGD8/P5sLVsq7JCUlsW7dOha+9x4LP/uMPVamValQgXbh4bR79lkatGjhk5+aNfhzSmIi+PmRmJjIlshIVl6+zIoSJfjpp5+4ePEi+fLmpVnz5rRp04b77ruPihUr2l2xUrnHhQvsP3aM7xYtYuG777I6JoZ4oECBAjStWZN7IiK4p08fqtWo4RPflXlM8ItIGWAmcDtggKnGmP/e7DkeHfwJCeDvT1JSEru6dmXF2rWsqFWLNWvWcPbsWQCqVatGs2bNaN26Nc2aNSN//vw2F62UD0hK4tyWLfx4+DDLly9n+axZ/GZdJSwkJIQWVavS4u67adKrF+XLl8+Vfwg8KfhLACWMMVtEpCCwGXjQGLPnRs/xqOA/eRKKFePChQtsGDyYn2fOZF1UFOvWr+f06dMAVKxYkebNm9OsWTOioqIICQmxuWilFBcvcmjtWpb/+afjD8H8+ZxISACgRIkSNCpThkYNGtC4Wzdq1KiRK7pdPSb4rytA5FvgPWPMshttY1vwJyRAdDSmXDn+OH6cdW+9xc9TprCuWjW2R0enXGKuWpUqNGjYkEaNGtG8eXPKlCmT87UqpVySdPIke9at46cjR1i7di0/zZvHIev/dMGCBWlQrBiRTZtSt1Mn6tatS3BwsM0Vu84jg19EwoA1QIQx5lyax3oDvQFCQ0NrHzx40P0FnTtH0rff8nupUmw+fpwtX3/N5s8/Z0uBApy5cAGAAnnzUr9+fSKjomjQoAH16tWjSJEi7q9NKeVeFy5w6Ndf+em33/hpxQrWfvwxu5OSSM7IsoGB1K1bl7rt2nF33brUrl2bgh5+DoHHBb+IFABWA2OMMV/dbNtsb/EbAyJcjo3l1+eeY2epUmxLSGDzzz+zdf16kv8C5c2blxqlSlG7cWPuatiQevXqERERkSs+AiqlMmAM50+eZMvu3WxcsoSNkyezMW9e/vj7bwAECC9dmpqRkdSsVIkaxYtTs107Spcr5zHfF3hU8ItIAPAdsMQY805G22cp+A8cIOHSJWJE2LVjB7v692dX8eLsSkpi3759Kd01gYGB1KxZk7vCwqjdvDm1776bqlWr6lh6pdQ1YmNj2fjNN2ycOpUtQUHs+P13/vjjj5THixQpQo3y5amZLx81Onemev36VK5c2ZYzjD0m+MXxp/AT4JQx5jlnnpPZ4P/vf//Lx0OHEp2QwBUr4AWoGBJCRIMGREREEFG1KhE1anDHHXf45DhfpVTWnYuJYeecOWwPCmJ7dDQ7li1j5x9/cDHVNiUDAqgSGUnliAiq3H47VUJDqXzPPZQoWdJtnxBuFPx2TOnYEOgK7BSRbda64caYRdl9oKSkJEJq1OCeqlWJaN6ciIgIqlSposMplVLZqlDFijR85RUaplqXdPo0v8fGsnv3bn794gui16whOi6OmTNnct4aVgoQFBRE5WLFCA8MpGLnzlSsWJE7QkKoWL06Rdz0hbLto3qc4VHDOZVSKguMMfy5ejW/rllDdNGiREdHE71wITF//83hq1ev2bZo0aLMmzePli1bZupYntTiV0opnyUilIqKolRUFC2SV77/PuCYjnr//v3smz2bmP37iSlShLJly2Z/DdriV0qp3OlGLX69pplSSvkYDX6llPIxGvxKKeVjNPiVUsrHaPArpZSP0eBXSikfo8GvlFI+RoNfKaV8jFecwCUiJ4DMTshfHIjNxnKyi9blGq3LNVqXazy1LshabWWNMddN+OMVwZ8VIrIpvTPX7KZ1uUbrco3W5RpPrQvcU5t29SillI/R4FdKKR/jC8E/1e4CbkDrco3W5RqtyzWeWhe4obZc38evlFLqWr7Q4ldKKZWKBr9SSvkYrwt+ESkjIitFZI+I7BaRf1nra4nIehHZJiKbRORua31lEVknIldEZHCafbUWkb0iEiMiw3K4rsdFZIeI7BSRn0WkpofU1d6qK3l9o1T76iYi+6xbt5ysK9Xz6opIgog84gl1iUiUiJy11m8TkZdT7Svbfo+ZqS1Vfdus7Ve7o7ZMvGZDUr1eu0QkUUSKekBdQSKyUES2W9v3SLUvO99jRUTka+v/5QYRiUi1r8y9XsYYr7oBJYC7rPsFgd+AqsBS4D5rfRtglXX/NqAuMAYYnGo/fsDvQHkgL7AdqJqDdUUCRaz79wG/eEhdBfjnu58awK/W/aLAfuvfItb9IjlVV6rXZgWwCHjEE+oCooDv0tlPtv4eM1lbYWAPEJr8f8ET3mNpnns/sMIT6gKGA+Os+8HAKasOu99jbwOvWPcrAz9m9fXyuha/MeYvY8wW6/55IBooBRigkLVZEPCntc1xY8xGID7Nru4GYowx+40xV4F5QPscrOtnY8xpa/16oLSH1HXBWO8q4FZrO4BWwDJjzCmr7mVA65yqyzIA+BI4nmqdJ9SVnmz9PWaytseAr4wxh6znJL9utr7H0ngUmOshdRmgoIgIjgbQKSAB+99jVXE0eDDG/AqEicjtZOH18uqLrYtIGHAn8AvwHLBERMbj6MKKzODppYDDqZaPAPVsqqsnsNhT6hKRh4A3cXxaanuTukrlVF0iUgp4CGiG4xNcMlvrsjQQke04/qMONsbsvkFd2fJ7dKG2cCBARFbhaFn+1xgz0521ufLeF5H8OAK0v7XK7rreAxbg+D0WBDobY5Ks956d77HtQAdgrdX9UxZHQzHTr5fXtfiTiUgBHK2/54wx54A+wEBjTBlgIDDdG+oSkWY4gv8FT6nLGPO1MaYy8CAw2kPqmgC8YIxJcmc9mahrC475UGoC7wLfeFBt/kBtHH+8WwEjRSTcA+pKdj/wP2PMKXfV5GJdrYBtQEmgFvCeiBS6boc5X9dYoLCIbMPxqXcrkJilg2e2n8rOGxAALAEGpVp3ln/6pgU4l+Y5r3JtH38DYEmq5ReBF3OyLhx96L8D4Z5UV5rn7scxSdSjwJRU66cAj+ZUXcAfwAHrdgFHd8+DdteVznMPWK9Xtv8eM/GaDQNGpdpuOtDRU95jwNfAY57y3ge+Bxqn2m4Fju4Uj3mPWesP4OgSyvTrlaU3oR036wefCUxIsz4aiLLutwA2p3n8Va4Nfn8coVaOf74YqZZTdQGhQAwQmWZ7u+uqmOrNdxdw1NpHURzhW8S6/QEUzenfo7X+Y679cte2uoCQVK/X3cAhax/Z+nvMZG1VgB+tWvIDu4AIu99j1nIQjj70Wz3ovT8JeNW6f7v13i/uAe+xwkBe6/7TwMysvl6ZfhPadQMa4fgSZAeOj2XbcHwD3gjYbP3wvwC1re1DcPR9nQPOWPcLWY+1wfGN+u/ASzlc1zTgdKptN6Xal511vQDstrZbBzRKta+ncPyxigF65GRdaZ77MVbw210Xjv7p3db69aT6Q56dv8fMvmbAEBwje3bh6FKw/T1mPac7MC+dfdn53i+JY2TNTuv1esJD3mMNrNdkL/AVqUYUZfb10ikblFLKx3jtl7tKKaUyR4NfKaV8jAa/Ukr5GA1+pZTyMRr8SinlYzT4lccTkWKpZnM8JiJHrfsXROSDHK5ljjUb4i4RmSEiAWkeT2/20MRU9S9Itb6ciPxizaz4mYjktdbns5ZjrMfDUj3nRWv9XhFplWp9ts4EqnI3Hc6pvIqIvApcMMaMt+n4bfhnXqVPgTXGmEnWY344JvC6DMwwxnxhrb9gjCmQzr4+xzGJ2jwRmQxsN8ZMEpG+QA1jzP+JSBfgIWNMZxGpimNCs7txjDlfjmM+HnCM5W6J4zyVjTjOLN3jjtdAeT9t8SuvJY655r+z7r8qIp+IyFoROSgiHUTkLXFc7+CH5Ja5iNQWkdUisllElohICVeOaYxZZCzABv6ZVRXSnz30RrUL0Bz4wlr1CY4pKMAxw+In1v0vgBbW9u1xnPR0xRjzB46Tie7GDTOBqtxNg1/lJhVwhOkDwGxgpTGmOnAJaGuF/7s4zvqtDczAcZ0Gl1n76gr8YC0nzx46KZ3NA60La6wXkQetdcWAM8aYBGs59YyPKbMuWo+ftba/0SyRbps9UuVOXj0ts1JpLDbGxIvIThwXqfjBWr8TCAMq4ZirZpmjAY0f8Fcmj/UBjm6etdbyBKzZQ619p1bWGHNURMoDK6z6zmbyuEplmQa/yk2uAFjhG2/++QIrCcd7XYDdxpgGN9qB1U+/2VpcYIx5OZ1tXsFxhaZnUq2uA8yzQr840EZEEowx3xhjjlp17bfmxr8TR5dQYRHxt1r1pXFMCob1bxngiIj445jQ7GSq9clSP+dG65W6jnb1KF+yFwgWkQbg6K4RkWqpNzDGJBpjalm39EK/F4552x81qa4NYIwpZ4wJM8aE4eiX72uM+UYc10vNZz23ONAQ2GP9UVoJJI/+6QZ8a91fYC1jPb7C2n4B0MUa9VMOuAPH9wwbgTusUUJ5gS7WtkqlS1v8ymcYY65awywnikgQjvf/BByzazprMnAQWGe17r8yxrx2k+2rAFNEJAlHQ2tsqtE2L+D4lPA6jotrJF94YzowS0RicExd3MWqf7c1EmgPjksC9jPGJAKISH8c87v74RhR5MrPpHyMDudUSikfo109SinlYzT4lVLKx2jwK6WUj9HgV0opH6PBr5RSPkaDXymlfIwGv1JK+Zj/B8+g4YFkZIulAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "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", "plt.legend(loc='best')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Suppose you also had some data you want to import:\n", "\n", "OGLE_data = MulensModel.MulensData(\n", " file_name=os.path.join(\n", " MulensModel.DATA_PATH, \"photometry_files\", \"OB03235\", \"OB03235_OGLE.tbl.txt\"),\n", " comments=['\\\\','|'], plot_properties={'label': 'OGLE', 'color': 'orange'})\n", "\n", "MOA_data = MulensModel.MulensData(\n", " file_name=os.path.join(\n", " MulensModel.DATA_PATH, \"photometry_files\", \"OB03235\", \"OB03235_MOA.tbl.txt\"),\n", " phot_fmt='flux', comments=['\\\\','|'], plot_properties={'label': 'MOA', 'color': 'cyan'})" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/jyee/MulensModel/source/MulensModel/utils.py:149: UserWarning: Flux to magnitude conversion approached negative flux\n", " UserWarning)\n", "/Users/jyee/MulensModel/source/MulensModel/utils.py:150: RuntimeWarning: invalid value encountered in log10\n", " mag = zeropoint - 2.5 * np.log10(flux)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOy9eVxU5734/35mANlxQWAIuMAwIAqaVZsmCopJs5k2baK5XRJrXJKYWJv87m2b2+/NvbfNvW3SW7MaE26StrdVm6RJzB4XRmPSpMYVERgGXEAGcGdAEJzz/P44ZyYDDItsKj7v1+u8ZuY8Z/k855x5Puf5fD7P5xFSShQKhUKh6A9M51sAhUKhUAwdlFJRKBQKRb+hlIpCoVAo+g2lVBQKhULRbyilolAoFIp+QykVhUKhUPQbSqkoFAqFot9QSkWhUCgU/YZSKhcoQoh7hRCFQojTQogaIcRKIcRwv/I0IcQaIcQRIUS9EKJMCPGsECLJKM8RQlR1cuzXhBAtQogGv2V3F7I8ZRzfLYQoEUL8qF35FCHEdkPW7UKIKe3Klxt1OCWEeEUIMcyv7P+EEC6jDg4hxH1dyJEjhND8ZD4shPj3dttIIYS1k/0tQoiXhRDVxv4VxrXI6E3dhRCxQojPhBDHhBAnhRB/F0J806/8XiGEp911zml3/ETvfRJCHDDuS2y7bXYZ9Rrnt+5aIcQmQ65TQoh3hRCZAeS/19j3rs7qaGxnE0K8YzxPx4UQHwsh0ttt09V9HCmEeEsI0SiEOCiE+Ce/skwhxFdCiBPGsiGQrH7b+z+fbuOZmtGuTlu72H+2EKLA2PeYcf3+RQgR2sU+0pDde5/y/crmCSFKjXrXCSH+IISI7up6XtJIKdVygS3AI0At8C0gGBgHfABsA0IAK3Ac+B8gydgnDvgJMM/4nQNUdXL814BfnYM8/w5koL+ETAVOANcaZSHAQWA5MAx42PgdYpTfaNRlIjACsAP/7XfsicAw43sGUANc2YkcbeoEjAeqgG/7rZOANcC+o4D9wJ+BVEAAw4H5wEO9rHsokG6UCeDbxn0JMsrvBbZ2c23vA/KN7weAUn95gCxjnQTGGeu+ATQAy4AoYCTwK0O2lHbHLwCOAe93I8c1wALjWMHAfwIlfuXd3cfVwFogErgOOAVMNMqGoz/DAjAbz8ieLmR5DeP5NK7tfUYdzN1dV+BO49wLgZHGunTgWSCti3MGfG6MsmQg1vgeaTxDz5zvduJCXc67AGppd0Mg2mgw7mq3PhKoA34M/B/wbjfHyaGflEqA/dcBjxjfbwAOA8Kv/BDwLeP7X4An/MpmATWdHDcdcLWve1d1Av4K/MLvd2dK5VfAbsDUx/vjq3u79SbgNuP8cca6Ths/v/3+BtxhfD8A/Cuwza/8KeAx2iqVT4EXAhzrQ+CPfr/HAhrwXeAsEH8O9RxpnHNUd/cRiABaAJtf+Z/wUzp+64OAB4HTXZy7zfMJhBuyJHZ1XdGVVmWg+9OD+naqVNptFwn8EfigL8/RUF6U+evC41r0N+C/+a+UUjagNxqzgTzgzcEXDYQQYcDVQJGxaiL6W6d/Erk9xnpvub9pbTcQL4QY5XfMF4QQp4ESdKXyQQ9lSQO+CXzRg83zgLeklFpPjt3J+drX3bt+D9CMrnDypZR1fsWXCyGOGqa9Xwohgvz2CwamA+v9tv8CiBZCTBBCmIG56C8R3n3C0Z+R1wOI+Ff058PLj4CvpJRvAsXA98+hutPRlcYx43dX99EGeKSUjnblE/1+I4Q4iX6dngWe6IkQxjX4EXovs7abzdOBJHr/39himPf+5m9qNOS4TghxCnCjK+kVvTzHkEcplQuPWOColPJsgDKXUR6LbiYCQAix1LDpNwghXu7heR419vEuf+jhfi+iNxgfG78j0c0N/pxCN8sEKvd+95YjpXzA+H09ujI908X5Ew156wEH8CXQqX3dj/bXbI5xHLcQ4pMe7A8d6+6VPxu9h/lP7WTZAkxCN01+F7gb+P/8yqcDu6WU7nbn+RN6QzobXdEe9isbif6/dQWQz/t8ePkReg8D4/OeLmtnIHS/3PPAT/1Wd3Ufu3sGAJBSDgdigKXAzm7EeNRQQo3oDfgvpZSebvbx1t3/Pq8x7vNpIcQPu9h3BrqJLgOoBt7zfwGQUm6VUsagK60n0XuVigAopXLhcRSI9X+g/bAY5ceM7wBIKZ8z/rAr0O3hPeEpKeVwv+UeACHEi37Oyl/47yCEeBK9kbzLr2fSgN6g+hON/kYXqNz7vU1DKqX0SCm3ov9p7zfO96GfLN637GpD3mh0W30T0BOF2P6arTOu2XJ0v1Bv6u4vf7OUcjXwMyHEZGNdhZRyv5RSk1IWAv8BfM9vt5sJ3Cv7E7qCuhfd1OLPCXSTloWOeJ8PhB4wMB5YY5T9BcgSRhCFaBs8MMavnqOBT9DNa6v9jt3VfezuGfAhpWxEV85/FELECSHG+Mvit+lTxv0JA64CnhRC3BSgzv54e1X+93mecZwd6P6cgHWXUm6RUrZIKU+i+6rGAxMCyH8Y+Iivr6uiHUqpXHj8Hf1N/Q7/lUKICOAmYKOx3NFx174jpVwipYw0Fp+JQuhRVjcBN0gp6/12KQKyhRDCb102X5uIioDJfmWTgVo/s0p7gtAd6Ugpb/KT5c8BZD2F3lje1oOqbQS+LYTo9JnvRd0DEQykdHYKdLu/l5uB9wPIcRDd3HMzHc2gjejPyJ0Bjn8Xej1B75UIYJcQoga9Rwd67wW/ekZKKQ8Z9RyBrlDWSSl/3e7YXd1HBxBkmCP9y9uYCf0woftJLpNSHvKXJcC1kFLKvcBnwC2dHM+Lt1fX5X8jUN0DbUbbe+WP7xlVBOB8O3XU0nEB/pnA0V870COsbMBJ9Oivy4x9YtFtya8Zv3PQI6NC2y2Cc4/++jlQBlgClHmjv5YZsi2lbfTXt9DNEZnoUUObMBy46GaheejmEzN6hFEjcHsncuTQNvorEv2t/ku/ddI4l3+dzcb1OUTb6K8o4DnA3su6T0OPdApBf6P+F/S3c69D+SYM5zi6WWUv8G/G7/FARbvjHQDyjO+pwFXG9yDaOuqvM67Tw0YdRqAHIpwE0ow6n0SP5krwWx40nqugAHWJBv4BPNfJdej0Phrla9AjwCLQ/Vz+0V+zgcuN+xANPINuYgrt5Fyv0dZRn4Fu2rvf+H0vupJp82wbZXOBevTorxHGfU4z7uG9nZxvIjDFkC8SvcdfCgQb5d8HxhjHGgtsBv52vtuJC3U57wKopZMbozcIe9HNO7XAKmCEX3kGumP2qNGQlaI7QJON8hyjIWq/WI0/bQu62cK7HO1CFonee/Lf3j/i6nJguyHrDuDydvv/1KhDPfAqX4cQjzb+oCeNskJgYRdy5KCbfrwyHEN/07e2k7X9cp9Rlgj8r9FANQDl6KazCb2pO7odfrdx/Y8bdZnut+9TRr0bgQp085e3oVpKuwYcP6XSbn0bpWKsuw49rLfBuHbvA5OMsnlGHYPbHSfUeF5uDXCOe4xzNLar65ju7qNRNhJ429j/EPBPfmV3ovciGoAj6C9I2V1c89f4+vn0Hu8JjMg9dKUS6D57Q7m/ZdwL7zOyE92XFdHJ+Wai/38a0SMs38Yv/Bj4NfoLWqPx+RJGVJxaOi7CuGgKhWIQEUJ8gK5UehTpplBcLAyYT8UYcVsnhNjbbv1DxujUIiHEbzvZ94DQR5PvEkJ8NVAyKhTnETv6wESFYkgxYD0VIcR09O7nH6WUk4x1uegDuW6RUp4RQsTJtjH93n0PoNuTjw6IcAqFQqEYEAaspyKl3IJuZ/bnfnTn3hljmw4KRaFQKBQXL4HGQgwkNuB6IcSv0UfWPiql3BZgOwl8IoSQwCop5UudHVAIsQhYBBAREXFlRkanuQEVCoVC0Y7t27cflVKO7q/jDbZSCUIP85uGnu7ir0KIFNnRBvdNKWW1ECIOWC+EKDF6Ph0wFM5LAFdddZX86ivlglEoFIqeIoQ42J/HG+zBj1Xo8d1SSvkP9PDQ2PYbSSmrjc864C30DKoKhUKhuMAZbKXyNnpMOEIIG/qgsTbOeCFEhBAiyvsdPQvuXhQKhUJxwTOQIcWr0dNJpAshqoQQC4BXgBQjzHgNcI+UUgp9oiJvvH48sFXok0b9A30eiI8GSk6FQqFQ9B8D5lORUt7dSdEPAmxbjZ7nCCllBW1zDCkUCkW/0draSlVVFc3NzedblEElNDSUpKQkgoN7mnO2dwy2o16hUCjOK1VVVURFRTFu3Dja5kEdukgpOXbsGFVVVYwfP35Az6WyFCsUikuK5uZmRo0adckoFAAhBKNGjRqU3plSKgqF4pLjUlIoXgarzkqpKBQKRTfkGIuie5RSUSgUikFGCMEPf/j17MZnz55l9OjR3Hrrrb51b7/9NtnZ2WRkZJCVlcXbb7/d5hhHjhwhODiYVatWDZrcPUEpFYVCoegCTdMY5XAwdvNmHA4Hmqb1+ZgRERHs3buXpqYmANavX89ll13mK9+9ezePPvoo77zzDiUlJaxbt45HH32UPXv2+LZ5/fXXmTZtGqtXr+5w/POJUioKhULRCZqm8eratcTb7VS3trLKbufVtWv7RbHcdNNNvP++Ppv06tWrufvur0dhPPXUU/ziF7/wRWqNHz+en//85zz55JO+bVavXs3vfvc7qqqqOHz4cJ/l6S+UUlEoFIpOcDqd7HO7WbVgARvy8lixYAHF9fU4nc4+H3vevHmsWbOG5uZm9uzZw9SpU31lRUVFXHnllW22v+qqqygqKgKgsrKSmpoarrnmGu666y7Wrl3bZ3n6C6VUFAqFohNcLheFKSloZjMAmtnMntRUampq+nzs7OxsDhw4wOrVq7n55pvblEkpO0Rr+a9bs2YNd911F6ArpwvJBKaUikKhUHSCxWIhq6ICk8cDgMnjIbu8nISEhH45/pw5c3j00UfbmL4AJk6cSPuM6zt27CAzMxPQTV+vvfYa48aNY86cOezevZuysrJ+kamvqBH1CoVC0QlWq5XM7dtZnJ9PWWoqk8vLmRAdjdVq7Zfj//jHPyYmJoasrCzsdrtv/aOPPsqdd97JzJkzGTduHAcOHOCJJ57gjTfeoLS0lMbGxjZ+lH/7t39jzZo1/PKXv+wXufqC6qkoFApFJ5hMJubPnUttbi6JISEsys1l/ty5mEz903QmJSWxbNmyDuunTJnCb37zG2677TYyMjK47bbb+O1vf8uUKVNYvXo13/nOd9ps/93vfveCMYEN2Bz15wM1SZdCoeiO4uJiJkyYcE775Bif9v4WZpAJVHchxHYp5VX9dQ5l/lIoFIpusJ9vAS4ilPlLoVAoFP2GUioKhUKh6DeUUlEoFApFv6GUikKhUCj6DaVUFEOaHFTKckU/sCFHXxTdopSKYsgyENllFYr+oqqqittvv520tDRSU1NZtmwZLS0tAPzjH/8gJyeHtLQ0rrjiCm655RYKCwsBePzxx3nqqac6HM9sNjNlyhTf8t///d+DWh8vKqRYMSTxZZd1uylLSWGV3U7m9u39OnBNcWmgaRrO46NwNUZicTiwWq19foaklNxxxx3cf//9vPPOO3g8HhYtWsRjjz3Go48+yl133cVf/vIXrr32WgC2bt1KeXk5WVlZnR4zLCyMXbt29Umu/kApFcWQxD+7rGY2s8njYXl+Pk6nE5vNdr7FU1wkaJrG2j+/irs6npSwMuzrVrE9PpO535/fJ8WyadMmQkNDmT9/PqD3Mn7/+9/7Ut3fc889PoUCcN111/WtIoOIemVTDEkGMrus4tLB6XTirt3HgqRV5MVuYIFlBfW1xX1OfR8otX10dDRjxoyhvLycK6644pyP2dTU1Mb8db7S4aueimJIYrFYyLLb2ejxoJnNX2eXzc0936IpLiJcLhcpIYWYhe6PMwuN1JA91NTU9KnHGyi1vXd9+9RZU6dOpb6+nhtuuIGnn36602NeKOYv1VNRDEmsVitjhGDhihXMWr+e5fn5/ZpdVnFpYLFYqGjJwiP1ptIjTZS3ZPc59X2g1Pb19fVUVlZitVrZsWOHb/2XX37Jf/7nf3Lq1Kk+nXOwUEpFMSQxmUxMHD+egpdeIunPf+737LKKSwOr1UpUfCb5VYvZcDSPfNdyouMn9PnlZNasWZw+fZo//vGPAHg8Hh555BHuvfdeHnnkEV577TU+//xz3/anT5/u0/kGE/UPUwxZNE3D4XBQXV2NzWZTCkVxzphMJuZ+fz6542sJiUokd86iPjvpAYQQvPXWW7z++uukpaVhs9kIDQ3liSeeICEhgbVr1/Lzn/8cq9XKtddeyxtvvMHSpUt9+//qV78iKSnJt0BHn8rPfvazPsnYW5RPRTFkOXv2LABBQeoxV/Qek8mEbeQxbCOPQT9GDiYnJ/Puu+8GLJs2bRqbN28OWPb444/z+OOPd1jvMWanPN+of5tiyOJVKmYjAkyh6DV59vMtwUWDsgcohizeNzelVBSKwUMpFcWQRZm/FJ0xlGa87SmDVWelVBRDFtVTUQQiNDSUY8eOXVKKRUrJsWPHCA0NHfBzqVc4xZDFq1RUT0XhT1JSElVVVRw5cuR8izKohIaG+iLFBhL1b1MMWZT5SxGI4OBgX44tRf+jzF+KIYtSKgrF4DNgSkUI8YoQok4Isddv3VohxC5jOSCECJioRgjxLSFEqRDCKYQ4PyN4FBc9yqeiUAw+A9lTeQ34lv8KKeVcKeUUKeUU4E3gb+13EkKYgeeBm4BM4G4hROYAyqkYoiilolAMPgOmVKSUW4DjgcqEnp7zLmB1gOJrAKeUskJK2QKsAW4fKDkVQxdl/lIoBp/z5VO5HqiVUpYFKLsMqPT7XWWsC4gQYpEQ4ishxFeXWjSHomvUiHqFYvA5X0rlbgL3UgA6TjIAnQaUSylfklJeJaW8avTo0f0inGJooEKKFYrBZ9D/bUKIIOAO4MpONqkCkv1+JwHVAy2XYuiheioKxeBzPnoqeUCJlLKqk/JtQJoQYrwQIgSYB6wbNOkUQwblqFcoBp+BDCleDfwdSBdCVAkhFhhF82hn+hJCJAohPgCQUp4FlgIfA8XAX6WURQMlp2Loohz1CsXgM2D/Ninl3Z2svzfAumrgZr/fHwAfDJRsikuD1tZWQB9BrVAoBgc1ol4xZPEqFTXjo0IxeKh/m2LI4lUqmqadZ0kUiksHpVQUQ5aWlhbg0pw7Q6E4Xyilohiy/M3oqSilolAMHiosRjEk0TSN5KAg0qZPR9M0NE1TvhWFYhBQ/zLFkEPTNF5du5bL09MJfuwxjsfG8uratcq3olAMAkqpKIYcTqeTfW43Lz/4IBtvuIEXHniA4vp6nE7n+RZNoRjyKKWiGHK4XC4KU1LQjJH0mtnMntRUampqzrNkCsXQRykVxZDDYrGQVVGByUjTYvJ4yC4vJyEh4TxLplAMfZSjXjHksFqtZG7fzsJnn8WZmcmE0lImxMZitVrPt2gKxZBH9VQUQw6TycT8uXP5oqCA1l//mh0uF/PnzlXRXwrFIKD+ZYohiclkwlRZSU1NDSOFwOl0qugvhWIQUEpFMWj88z//M2lpaYNyLk3TuPyb3yRn0SKacnJYZbcPWlhxjrEoFJciyqeiGDSefPLJQTuX0+kkePx4Xl62DM1spsDjYXl+Pk6nE5vNNmhyKBSXGqqnohiSuFwunJMmqbBihWKQUUpFMSSJi4sjtbBw0MOKczSNWoeDYZs343A4lB9HccmhzF+KIUlcXBzu4mIWrliBMzubKRUVTIiOHtCwYk3TyFy7lnS3m/KUFFbZ7WRu364izxSXFEqpKAYdKSVCiAE9R0NDA2teeYW0tDRs3/42ixYswGq1Dmjj7nQ6we0mf8EC5cdRXLKo1yfFoOOdO34gOXbsGFJKHA4H406fxmazDXhvweVyUdbD9DA5qAgxxdBEKRXFoOOdPGsgOXLkiO/7QM+nkmMsFouFNJUeRnGJo8xfikGnpaWFiIiIAT2Hv1IZLGe51WqF7dtZmJ+PMzWVKeXlnI6OZpHVit1vuxxgFzBlUKRSKAYX1VNRDDqD2lMxm3l7kGZ+NJlMrJ47l4LcXEwhISzKzWXf3LnQ3uymaYx3OBirIsQUQxDVU1EMOoOhVKqrqwkNDcUTHQ3dKJUc49Pei/NomsYop5NIl4uS+HjigITaWiotFlJSUnxlDovFF3mWuXYtmooQUwxRlFJRDDoDqVRyjM+EykqSk5MJHjGC8REROByOfo/+8s4wGe92UzZ+PK++9x55oaEciYvjqm3b+BUQHx1NmZ/y+Obll4PbzctGhNgmj4fF+fl84HTypooQUwwBlFJRDDqtra0Dfo7Kykpm3347Zy0WyrOyBqRH4J1hctWCBVjLy0kpL8cdFUVCTQ1N4eGc0TRW+SmP5fn5jNi3D0e7CLGy1FQSa2pAKRXFEED1txWDhrcxH6ieiqZpjDJ8FcOHD0ckJ5O/bBkbZ89mxYIF/T6lsP8Mk5bDh2kNCmLcwYMANIWGUpqe3iG8GCC1XYRYWnk5DSpCTDFEUD0VxaARFBRES0vLgCiVNqaolBSyr7uOPRkZWMvLSXC5qLFY2JOSQk1NTZuBiP4+Ea/fY6ah/OzdnNNisZBlt7OptZWrC79iTEMlGRHF7GnJ5mDrOFLKy0lPTia+tpba+HiynU4ypk/n048/5v6VKylNTyejtJTQ1laOpaT0+zVRKM4HSqkoBo3g4OABUyr+pijNbKYqMZHb3n+flAMHKE9JIcduZ0R9PXG33urbp70i8prICBSxFQDvDJNLX3iB5NNVPDD2ecxCY6bcxO+rf8LRs7HMWbeOA2PHklVYSExrK1JKmkNDKcjNJb6ujvV5eeTY7ZyqqFDmL8WQQCkVxaARHBwMDIz5y98UBaCZTJwOC2uTMmXpiy+22ae9Itrk8bD85ZcZt3kz0mTiuxYLm6xWJptMbXotmqbhfONOyo6E0ei5hiNVVdyStA+z0EODBZKw1jNchov00BKKDk7icGgixzDz3scf48zKojQjg9KMDACSqqqI6QefSo7xae9iG4VioFE+FcWg0RelkkPXaU0sFgtZfr4KS01NB59GcXo6dXV1vn3aKyIpBE0eD+ElJVS3tJD80Uf86IUXGFVS4htLomkaq174PR+WjCa4tYbW40WMbG6gyG3DI/W/k6PRhgQeGPs8N8Su5+HLniaxsYrw06dpbGjg8p07CTpzBtB9KrayMmqVT0UxRFA9FcWgMZA9Fa8panF+Po7x48navRstPJxNM2eimc1fp0zJzfXt4/WJbPR40Mxm0ktKEFKy6r77uHfNHxjbvJ9Jw/ZS8lYZr4xMZcbMm3j88ce5Ii2Sn6TlYxYas0YX8EzEck7JGJ6rXEpm2D62N1zFlZFf+XouUoLnbAiJZhfpMaWUNqbz+JP/QcHV08kuKiK4tRU0DU3T1FgVxUWPUiqKQSMoSH/cBiKk2GQyMX/uXN4rLKR12TKipk9nzGWXsfCZZ3BOmhQw9b1XES3Mz8eZksK0bdv48uqrSS0vJ/twIQ8mP+fzkTx/aClPPfssR44cIfuaOp/CMAuNScP28HHQt2k+aaauIY6zZjPFpzOZKTdhFhpbjs8gzHya+8eu1I8Xu4kXDjzAN3d8wUhxnMywfVz2VjUvRiRx35JlhISEdFvfHOPT3u9XUqHoG+q1SDFoeJVKoJ5KDn3P2jvTZOKTgwfZsnkz06dP58fz5vHZO+8Qsm4di3JzO4xR8SqiL3NzGeV20xIRQcr+/UwqKmJCWFsfSXxwDaPi4/nBD37A/tZsn6nLI01UtGQzdfRogs6e5WRMDGdCQzlijuWFgw+y/kgen5+6lvSI0jaKKG5YLbHyKEuTn2P26A08OOY5zjZW8l8rVgxKFmeFYqBQSkUxaPTW/OU//iRQrixN03AY5ZPWrycmJobp06djMpkor67mq6NHO019bzKZ2G+zcXLUKHZlZuKOisLqdFJ8OhOPNKFJwZrqeRxpiSM3ajMN1TtolMPJr1rMhqN55LuWEx0/gdjYWI6PGEFEUxNN4eH8/fJpVJkS+aTlBtzBkZQ2prdRRAebxpHRTtFMCC+m5exZtmzZ0pvLq1BcEAyY+UsI8QpwK1AnpZxkrFsLpBubDAdOSik7JGsVQhwA3IAHOCulvGqg5FQMHr1RKp2F/Xp7HZqmsfbPr+Ku3s4PwsrYE5vKsUX3+npFCNFl7q8c4BRwKD6eGzZu5PlFi7j+00/J+WIzzx9aSnxwDcdbR3H/2BcMU9hGXq5eju2aOzGbzeQmJGC1WnE6nQz/8ktOREf7Is6qkpPJW7+eTTNnMm/dWlYefID0iBJKGzNoCI6k+HQGM+VGzELDI03sac3mdFgYDoeDmTNndnlN2o+tUb4YxYVCj55EofMDIcT/M36PEUJc081urwHf8l8hpZwrpZxiKJI3gb91sX+usa1SKEOE3igV/7DfDXl5HUbGO51O3LX7WJC0itmxG/iJ9WUuGyl95cJk8imV4cYSCAGENjVx3yuvYPZ4aDCHU6PFsedMNumRbXsU1mF7MJvNTJ8+3dcDslqtxAQHU5aW5osmc9hsSCG48eOPqYmKo641FntDDn+b9h2e+OnPOGIezbNVD/HJ0dmscC2nOD6TvRMnEhoa2un18ClZu53q1lZW2e28unZtnzMd56AmDVP0Dz19vXkB+AZwt/HbDTzf1Q5Syi3A8UBlQp9L9i5gdQ/PrxgCdOVTaU+OsXQYf9JuNkWXy0VKSGGbRj89vPjr2RaFQHbT4JqB5Npadk6eTEFuLq3DhrH96qs5FR1Nqc1GceOENqar8pbsgBNvhUVFkVFa2iYFS8Tp0wR5PGAy8fdp0/CYzUzct4/czZtpFKEcCh7DE9m/4KU5i3ht3o/IKC3l6quv7lTW7pSsQnG+6alSmSqlfBBoBpBSngC6D1HpnOuBWillWSflEvhECLFdCLGoqwMJIRYJIb4SQnzlPzGT4sLjXJSKl/bjT9rPpmixWKhoyWrT6DuaJjJ69Gg0TSM1MZGrRo/G4XBQr2mc8jt2DvpkWQB1FgtZ+/fjTE1ly/TpSACTiVHHj1MVkcTzh5ay/uhy09gAACAASURBVEgeq6oeJjp+QpsoMtAb+wNScmzkSBa99BJ569fz8LPP4gkOpj4mhlWLF7Phxhv53fLlhDY1kVpezjfGjSPY42FCcTFjDh5k4csvE+F2s/uTZzqdZ6U7JatQnG966lNpFUKY0Rt7hBCjgb70t++m617KN6WU1UKIOGC9EKLE6Pl0QEr5EvASwFVXXTU4szEpeoXX/HUuIcX+40/KUlOZXF7eJjTYarXyVdwEXjj4IBkRxRQ1T+SkNpy/fvABl+3cyTduuYXySZNYVVDAgs2bKUtL47uJiaxNSWFURQVTDx8m+tgxot1uhh075svJddVXX+GOiiL/vvuQQmBzOJizbh0T4kYz9/v3dPBhuFwuClNT2TBzJmlOJwk1NZwZNoymsDAOjRnjUwJnhw1j15QptAwbxt1SsnP4cA6OHUuCy0Ww2UOCuY7RpgoK3lnFV/ETuOKa66itrcVi+E7aj63xKtm4GTMY5XAoP4vivNNTpfIM8BYQJ4T4NfA94F97c0IhRBBwB3BlZ9tIKauNzzohxFvANYAKibnI6U1PxRv2+4HTSWJNDYtyc9s0mCaTiZFxSXx1tJ43J9+BKzGR8vHjefjZZzlx6hT5Dz+MFIK4I0cYeeIEQWfPEl9QwH+9/z7xYWGkuN00hYVRN3o0YZGRbMjLI76ujpMxMTj9/COlGRnsrKxkYnV1wMba19jn5uKw2XCmpjJt505C3G5SKyoo8FMCaU4nf581iziPh/SNG9loOOXnbFvHQ0nPYhYarVoBKw48Qn31TqzhDuzbstken8mdd9/z9dgaY8rijKgotu7YQXxDQ8BgBoViMOmRUpFS/lkIsR2Yhe7T/LaUsriX58wDSqSUVYEKhRARgElK6Ta+3wD8Ry/PpbgAOdeQYpPJxDGbjY02GwfpOOBv27ZtFE6bxuacHN+6k8OHU5WcjGY2Y3M4iGpo4MVFi0jdv5+YU6cwu1yUWq1MKC1l5ZIlXLd1K8dHjfLl5DJ5PEwsKmozIj+jtJQxEycGlLHNQMrUVGzl5VhjY9l7+jTRJ0+yZNUqHGlpZDoctI4YQYzVistuJ6ypiYX5+ZwJCSE7eI/PN7T/dCoRZjcLk1cZ0WGbyHctp6Kigvlz5/Jboze0KDcXTdN4ecuWDnO3OJ3ONhmZFYrBoMvXGCHESO8C1KGbrP4C1Brrutp3NfB3IF0IUSWEWGAUzaOd6UsIkSiE+MD4GQ9sFULsBv4BvC+l/OhcK6a48JBGFFZ/pmmpq6vjs88+I72dg3z4yZNkOp2YPB4s1dVUjB3LPX/8I7M2biS6vp4Sm40rd+yg1GZDM5upsVjazHOy9dpriWpo4IGVK8lbv54HVq5kdGsr119/fUA5vD2qI8b89LUzZiCEoFXTODhmDJFuN1cWFrJk5kz2zZsHJhOVlZXszM6m1GYj9MwZipon+nxD1WcspIWXtQlASA3Z09Z3YlzPmpqaDn6W3amp/Fr5WRTnge56KtvR/SgCGAOcML4PBw4B4zvbUUp5dyfr7w2wrhq42fheAUzuXnTFxYZXqTQ2Nna5XaBxGHaTKWDI6y9/+UsKCgqYNnt2mzlKhgOJMTEsfOYZWqKiiD9yBHd0NCuXLMFaXs6NH32Ex2Qi3eFg08yZlFmtXLltG8uefpqTw4czur6e0OBgzjY3M2nvXkI9HhISEro0J3l7VLtsNmY5HBQ3NLDqvvt8PZ3F+fmYzWZfWv3k5GQu376dUzExlKWmEre9zucb2td0BWbtNDNjN/nGsZS3ZDMjLo5X1qzhxqNHcUdF8b+7d3MmNJQUdGXqM7GVl1Prl+dMoRgsulQqUsrxAEKIF4F1UsoPjN83oZuxFIoe441mOnXqVJfbdDbYsf0cJ19++SUvv/wyy5Yt41+XLePWTz8lubKSOyZO5Prrr8dkMvHkrbcyc9IkPLGxlNhs+qRd1dWYPR72ZGUxad8+XRnZbCQdPkxTeDiVycnEFhXRqGnUx8RQnpJCakUFsqYGh8NBhpGyvjOmAA+7XPw6wLTBNX4p7i0WC01+6fntOTkse/ZZEpNmcOPEiezYtpUX9i8iPbKM0tOZxCZPQNM0yqqrORMdTWVyMiEVFUSfOkXDiBEszM+n3AhmqIuOZlO7CDWFYjDoqaP+ainlEu8PKeWHQoj/HCCZFEMUb0+lvVLx75ls1rSOc5wY/gH/+UZOnz7Nj370I5KSknj88ccJCgriC+PN/BW/YzuOHOEbHg/Hxoxhyu7dpO7fT3lKCh6TiYyyMp574AGu+/xzMkpKaAoP54X770czm4murye+trbNfCxLVq1i37593SoV6JgB2dt78GZJ3gX8d10drowMn+LxBAdTmJ3N7SEhZGRkYLPZKCu7npUrV/Lhhx/yz/88meLiYtzh4R3kOj5iBEVZWUwz/CxTrVbkOTjpd3W/iULRI3r61B0VQvyrEGKcEGKsEOIx4NhACqYYeniVysmTJ33r2o8Q/2THDgrHj+92HMZjjz2Gw+HglVdeISYmBtB7CP45fzRNIyM+ntGhoSS6XL5ewca8PF544AHCTp/mvpdewmwMUvSffwUh2oyO18xmHGlpPa6r1WolMyqKxfn55G3YwOL8fGiXJbmhmzE4JpOJ9PR0nnzyScaNG8eSJUuoq6sLLJcQOGw2Dhqj/HuiUHKMRdM0xjscTO0kt5pCcS70VKncDYxGDyt+G4jj69H1CkWPCNRTaT9C/K1bbiHd4ei0oQX46KOPWLFiBQ8++CB5eYGtsJqm8cqaNcyeMYPgsDAimppwGE550BvjA2PGgNErKbr8crL27/edd19GRofR8ZkOB5mZmT2qq9dxX5ubS2JICLW5uewzQnzt6MrvmKF4FubnM2vDBpbn53dIzw/6+J61a9dis9l49913uWLXLn3uF03D5PGQ7nBQ5CdXDtDQIykBQ6lPtdvx9GPaF8WlS09Dio8DywZYFsUQx6tUjh//OntP+xHiDpuNvI0buX/VKkpsNmzGOAxt58+5/7NI/i91Ad///vfJzs7mt7/9bafncjgcuu8hJgZzfT2hZ86QUVJCZXIy8TU1pDmdhDY3U5KRQWpFBU3NzUyIjPx6kKXTiWhtZdFLL+G0WslyOskcMeKcQnS9jvtjne3jNwbnsgBjcPyJjo7mkWX3c9JVSFbIJyS9VUmdeTSN5nDMgJCS6Zs3M8pi4ZjV2sH/1BmjDKX+sp85TYUjK/pCj5SKEKIAYzS9P1LKzlOpKhTt8CqVuro6PB4PZrO5g+9BSEmo2czR9HRMZjO1M2YQvm0rWyrjGR9WRq59NZd95xYe+ZdfEh4e3um5ioqKfL6HpU8/TfOwYUQ0NpK3cSN1sbGENjfz4uLFpO7fz9mgIKbs3s21s2fzYVCQPshy5kwWpaRQW1FBck0NS2bN6tEodfs5rvdXPF014U6nk7P1+1lmzfdNHPZs1UO4zVGENDXxrU8+4VR0NDG7dxM/ahSuK68ksba229H1kV2kfVFKRdEbeuqof9TveyjwXUDNJKQ4J7wmFY/HQ11dnS/1iG8aYGPQ4ISYGP50/fWcqqggtbCQetc+FibpgwBnjtrEqpCHfArKH7vfd+HnEwlraUFomi+k+LqtWzk+YgTzX3uN6Pp6KpOTaQkK4oPNmzm2eDHHbDasRvCAcLk40oe0J/YA69qHTNODnkWgxJmZYfvYmngdKfv30xQWpkeDlZczorKSWadOUZaWxoubNpH51Vf8eN68gPI3dJL2JUGFIyt6SY/+JVLK7X7LZ1LKnwJTB1g2xRDDXxFUVekJFfx9DybD93DPnXeS+frrTLXbOVNfT+qwto2pLXwfLpery3NlZmaSafhmzprNNIeF+RzxtfHxXP/V52QdLeTG4I+ZeuhLWqKHUS0lLzmdbPILHjC3tjK1H/0MgVLXZ65dC90cO1DizD1nskkrK0MCm6dPZ9PMmWy5/noaIyP1BJazZ/P0woUUnjiBw+EIeNye+nUUip7SU/OX/+h5E3rero65vxWKLpBSMnbsWA4ePEhxcbEvxbv/oMEpQIXDAYad31pezsR1e33zvXukiX2NEzA5HL6xKIGw2WxM2r6dJc8/j8dkIvT0aTJKS9k0cyYSGCGO8+CY5zELjVlyIytcy/nHmGt8UWb+Yc0FHg+R/eRn8A9M8A+ZfqldyHR7rFYr2+MzebHyIWxhRexuyabRE0midJExrISE92oojs/kRPhwSvyi2DSzmX02W4dQ6Bz0MOIphlJ/zulkdDd+HYWiJ/TU/OU/sv4ssB9Y0OUeCkU7pJSkp6cTGRFK8Rev45g27esGTNMY73Qy1uVi7941OG136AMGrVYORo3j+UNLyYgoprAlm+KkDIY3nu2ykdc0jeraWoYBZ0JCiGxoaJNna0J4cZveT1bIHo5VWkiYMqXL9PJ9VSq9PbbJZGLu9+fzxhtvsNkRS0tICEnyMA+Mfd43I+WzVQ9TGpVGfG1tm5xltrIyaDf/y25Nw2Jcb6fFwgGrlf3d+HUUip7QU6UyQUrZ7L9CCDFsAORRDGGklEzJsnHTteOZONyJfd2qrzPvvv46mttN+fjxbGscy+U7d3IoORmHzcYe2yRq6uJ40/JdahISKLNaydu0qcuG2G63czQ4mPWzZxN/+DAzPvuMnZMnc2jcOCbu3cveI5Pa9H5KGieQGPu12Weg/Aydpa7vybFNJhONjY0cHzGC2oQEZh/e0EYxTgwrovZUHB6zmaXPPUd9dDTR9fWEtbS0CYXWNI3b1q4lzO2mwshacNv27bwbIGuBQnGu9FSpfA5c0W7d3wOsUyg6JS4ujsQYjaUpLxuNeQH5ruV8+umn4HaTP38+9675I7GmI2QH7+Gyt6qoCk9GyGCag4N583vf63FDXFhYiBYSwowtWwhqaaF52DDSy8rYNGsWZVYrI/98nOcPPciEiBL2NU0kcmQK98yf75sauKs5XPpCX4/d3NxMmc3GobFj2bM/q+0c9y3ZvD1nDre/+y4es5mqpCQySksJDwpqc3yn00m4XxjxJo+Hhfn5jOvGBKdQ9IQulYoQIgG4DAgTQlyObv4CiAY6j+dUKALQGh3NhGhHh8y7hw5lUZaSQur+/Uyo3cdPLCt8YbPPVT1M1tW3UVFTc04N8ZkzZxDBwXpI8bPP0hgeTpTbzQMrV1KSnk7YySZqtXjez76DyMREXm83R0tXc7j0hb4eOz09nYx9+yjIyaE4PtOnGPe0ZFMcPwEpBO6oKF5cskRXGDNnsjw/n4qKCl+vrrq6Gmc7E5wzNZUYv7xkCkVv6a6nciNwL5AE/I/fejfwiwGSSTFEqT15kr31Nma1ybybhW3MGNLKymgNCiKrQ9hsEUFB3+HH8+bx4Tk0xKGhob5or/qYGCIaG/nd8uVc9/nnjDl0iLCmJlISEvgsJ4czdAyD7On4kd7Ql2NPnz6drbt3+waHJhyopqo5iXVz5lBmtepK0y+fWHufjaZp7C0rI6OlhYLcXF/Pz1pezmYVRqzoB7rLUvwH4A9CiO9KKd8cJJkUQ5T9Lhdl9TbyXctJCdnDnuMpVB1vYEfZGhpDQph84gT7yGSWn0mnqCWbm42U8+fSEGdlZeEsLmbTzJkcGDOGIE3j7LBh2I2GM2/9er7lzfN1EREUFMTPHnyQf1q3jsziYhJGjKCqoYG8DRvI3rWLsKYmsvfsodLwRwkpmVRWRsKsWYBu+trv8dA8bBiLV63CYbMxobgYKQQJ1dU4wKewc4xz2s9TXRUXJ92Zv34gpfw/YJwQ4qfty6WU/xNgN4UiIFJK3nBWsnzhb6ipqSGqpIT1bzyDy+VCTplC3OTJ3HrazEuVS0gLd7CnZTLV8ROwWq1fh8D28FwzZszgsz17eGDlShrDwghvamoTEZVeWoroZBbHCxlN0/jTm28yyu2maOJEwnftwh0eTqTbTfaJGqJCGsgM3sfYtw5SFZFMozmcRMNPBLrp64zHg1kIjsbGMnnXLoI8HtzR0ZjPnu1yqgGFoid0Z/6KMD4jB1oQxSWAlCAENpsNm83G9OnTWbRoEaCPm9gMlGsaRU4nv66poSEhoY2v41wwmUxY4uI45XIR4Xbjjo72TfVrLS9HCnFex2LYe7mfd5yLN/X94cREbnv/fc6Eh5PYVMPSpOeMzAMbeabyYWqIZ3xaGk6nE6vV6huAmm9MHpZeUkLehg2sWrSoQ+4v5V9R9IbuzF+rjM9/HxxxFEOVHKBRSs4IEbDcO5+HNN6qG9DzUjnRzTHn+tbsdDopaWxkxfLlWMvLySkoYPOMGcTX1WGfMYOZdjsWi6UPNTo/dBjnYjJxOiyMvZMmcfOe933+KIEkRGshicPEHN6IfX822+MzSUhKaZOtOb62tk3K/zZTDVwgSiXH+LSfRxkUPadH/1QhxGghxC+EEC8JIV7xLgMtnGKIYfRUukIESGPyag/SmLTHv/Ets1pxR0cza9MmQs6cIWfzZoiJ0acppuskkJ2VnS8s7eZgsdTUUJqejisxkUK/NC6ORhsSwQNjXyAvdgMLLCuory0GaJPivzY+noxuphpQKM6Fno5TeQf4FNgAeAZOHMWQRtM6VSpTNI1ap5PJhYUUnTrFqoUL26QxGXWO5hiLxUJWQQGHkpKIr61l5+TJ3LB+PVH19Ryw2YgXwmcSuphSkrQf53JFURHNQUG+EONnDi9j0rBCdriv4IqonR3Ct02m23y5vpypqVidTuJaWny/pxjh2ousVrai7N6Kc6enSiVcSvkvAyqJYugjJSJAA65pGplr15LudtMSHExhamoHc4zlHM0xKSkpxLz3Hjd/+CEnY2K4cscOpKYRduYMCWVllBkjyb1O6YtFsbQf5zL/xhvZumMH973yCs6UFBIOH6GgMZcaSzwhx1raDI4sbZzAbIuF66+/nuecTqJratg8cyaFKSlcXVHRJvfXny6S66G48OipUnlPCHGzlPKDAZVGMbTpxPzlcDgQp06Rv3Ch7v+w2zH5pzFxOjmans7UzZuJtFjQetC7cDqdNGoajdHRHBozhuCKCkYcP050fX2HXtDFNiGVf3h1BnryzCedThJqaphitfL5oUOMPH6ck9oInq16iImhRRQ3TmC/q5nZxv77bTZO2WyY0RuB/Tabyv2l6Bd6+jqyDF2xNAkh6oUQbiFE/UAKphhiaBq2xESui41tMw+6pml8uGEDDqN3Uma14o6KYsmqVeStX8/yl18mpqmJ8LIyPK2txPcwDb3/JF0b8/L43/nzAXAG6AV5MxNfrJhMJhw2G1umTyc2Npaqyy6jKTycXz/yc1Z+5wF+PfkxDoaO49O/b+f2229vM51zf5u3cvjasa64NOnpfCpRUkqTlDJMShlt/I4eaOEUQwOveStnzhxMt9/eZh50p9NJjaaRYjifpcnEX7/3PUJaWrC63dySns7JiAhWGcph1YIFFNfX6yGvXeA/SZfQNO564w2ahw0jraxsSDulExMTGVtVRWl6Op7gYBw2G5tzcijMzuahhx6ivLycJUuWMHXTJnIKCphqt+s9RU2jgc4VQk4XZQqFPz2dTyVQ4shTwEEppZoBUtElTqdTTxj58MMdxkK4XC52ZmQwsaiI+1eupDQ9nfTSUkJaW9k5Zw6mzz7rVar4zMxMMjduZNPMmVjLy4lyu3lu6VLueuMNFubnU56SQpbTyYThw4fEhFQxxqfVaqU1OpqM0tI2aViyy8uZOmMGT/zHvyKaD5Nd+F8Un55AVUQyq0qjuW37dlarAY+KfqCnT9ALwBfAy8byBbAGcAghbhgg2RRDBJfLRVknisFisXB5SQmnw8PZkJdHy7BhbMjLoyEiglEVFR1CaHvau7DZbGSNGMGil17im1u3UjF+PJ7gYNbMnUtBbi7JVVVcFx9/UTnp/bHTNtx5irGYTCZ23XMPmhDc/+KL5K1fz0+M2RwBooNO8ZO0fGaP3sCDY54nsamaXampjKir4/a332aU1zSpaaQ5HFy+eTPfdTjOOaRbcenS03/TAeByKeWVUsor0Z/fvUAe8NsBkk1xgZLDuZlCLBYLaZ0oBqvVSoLJREVKCqUZGWyZPp3SjAwcNhuRNTV6CG1UFIvz88k7h+luTSYTP543jy9mzaJx+HCyyst95jVnaiqhLS1MmjTpolQoXaFpGj988008ZjN1o0czec8exptM3HPnndTW1nac5z60iG988QXDzpzB1NpKfEEBr6xZQ+aaNUw3plOO7+GUxwoF9Dz6K0NKWeT9IaXcJ4S4XEpZIboZzKZQpKSkEPH++z7zVkZpKTGtraSkpGAymbgpLw+H3d4mN1daebk+b30fUsV7o5wOWK08tXatbyyGrbwchtg87Hbj0+FN42KkYTEZpsa5FRVgsXDjtiw8Rphxq2Zmp/ty4s11pIeWUFiZRXH8BMrq6hChoeQbUXIFHg+LezFWSHFp0lOlUiqEWIlu8gKYi276Gga0DohkCmBopKioqKigMTycTUaalPV5edxqt/vm+LDZbDTu3NlmAN7p6GiOGY1+X9PQS0MxPeV0El9TQ21uLscuskGPPaWz6YotNTUcvO46ouIzWXXoftIiHOw4fTnh5iYeGKNPSTxLbmSFazkV4Sm+aDzvMcpSUxE1NeTYbBf1s6gYeHr6r7oXcAI/AZYDFca6VkBNwqDoEpfLRVlqahvz1h6r1RfKazKZeNfwdZhCQliUm8u+fnYam0wmXDYbu6ZP55jNNmQd0p35oBoSEsCY5z5j6h183HID/0i7hvSIkjbmsKyQPQx3u7EZ5kLvMVLLy6kcQlFyioGjRz0VKWUT8DtjaU9Dv0qkGHJYLBbS7HY2dTEvuzTGWtQavRH7AMjR07T5FzPeNC7+GZknREfzlZ+pr6KmBgGMOXSI4pa289eUNE4gYdQITkRGtjEXNkRH47Raie/i3JqmccrpJM7lwmGxnHMKnBzj096LeisuHHoaUpwG/BeQCYR610spUwZILsUQwmq10rB9OwuffhpnVhZTKiraONs1TWOW08lxl4uTPRwxrwiM1wflTcPyZW4uL/ilXXE6nRQ3NPD0smVYy8u59ZP3eOHQg2SEF1NUn0aVKYYn7r2XuysqkEVFRNfXUztjBptsNmQX90QzEoFOdbtxnkMKnBzj095/l0BxnumpT+VV4N+A36Obu+bz9Xz1CkWXeM1bcddeS2JVFYuWLPG9xXobo3i3m1MpKXzDbufVfszHFah3Yu/zUS9svAEKDTYb19HWxu1yuSgcP57U/fuJr63lvRtuZfLu3exsvpLPdmzC+eeXGBkxkviwMI6OGEFyZSWeo0dJkZLIujpOdaL0vfO8vGzM81JwkabAUfSdnv5rw6SUGwEhpTwopXwcmDlwYimGGtJkwnHwIJ81NGCz2XyNkrcx6m7EvJ3eKYPe7jdUiY+P5/Ldu8kpKCC+tpZb3n+f1P37OXX11dS8+CJZt93GwbNnORkVxYiTJ9mdnU3IqVN8Y9MmgltbmdpJmpzOAgQu9hQ4inOnp0qlWQhhAsqEEEuFEN8B4gZQLsUQQtM0xjscTL/iCtJCQ9s0SIPZGNlRCgagOTQUt1dpTJ5MfWQkSZ9/jmnYMLKWLOFEcjJRDQ3kL1hA5dixnIqJYeWSJWzMy+PlTpR+bwepKoYePVUqPwHCgYeBK4EfAvd0tYMxkVedEGKv37opQogvhBC7hBBfCSGu6WTfbwkhSoUQTiHEz3oooyIAOZx7zqbe7NMZPlu73U7wsmVcZ7W2edNVjdHgUltbS01cHFFuN5unT6c1OJgNeXl4zpxhnNPJ6fHjSa2spNxQ9Akul+87dK70vYNUF+bnM+scBql2h6ZpjHI4GLt5c5tEpIoLl54mlNwmpWyQUlZJKedLKe+QUn7RzW6vAd9qt+63wL9LKacA/48Ao/GFEGbgeeAm9MCAu4UQmT2RU3Hu5NB7BdKTff1t7RtvuIH8Bx9s86bb2xHziq4RRpoVb2OMnxJPOXSIoLNnmbFlC8GtrczYsgUJjHa5OGa1cllMDBmlpZg8HmosFlJ7oPS9AQJf+oWF99Uv5vO3tZsFVCmWC5suHfVCiHVdlUsp53RRtkUIMa79asCb3TgGqA6w6zWAU0pZYciwBrgd2NeVLIoLk67MW17fSm9HzCsCo2kat61dS5jbTYVfJNa+uXNJSUkhtLUVd2go+X5O9ftXrsSkaWAysWT+fB75wx9Y/OKLOKxWok+c4P4XX6QkPd0XohxI6c80mSi02Yjsp3lZ/P1tF/P8N5ca3UV/fQOoBFYDX9L3iK+fAB8LIZ5C7yVdG2Cby4xzeqkCpnZ2QCHEImARwJgxY/oo3oWFpmmMcjqJ7GXc/0DL1GCx6KPeu5DJYrGQZbezsYsxKn0dMa9oi9PpJNKvMS5obWXpCy8w4e23+XTECJpDQihNT2+j6EvT033PVlBQELvnz8fldJJQUMC7f/gDIjOTuIwMX4jyYDyH3b2QKC5MunsyEoBfAJOAp4HZwFEp5WYp5eZenO9+YLmUMhl9ZP7/BtgmkOKSnR1QSvmSlPIqKeVVo0eP7oVIFyYXYtffX6bDra2MNhINdiVTG1v7J5+w6PnnlXlrgPHPCu2dS6YlOBhnVBSf7NhBRVISKfv3tzFppZeVUWexfH0Q78Rfixdj+vnP0QoLYc0a5FljposNOfoygCh/28VJlz0VKaUH+Aj4yMjzdTdgF0L8h5Ty2V6c7x70WSQBXgfyA2xTBST7/U4isJlsSHMhdv3by+RNNNiVTG0G4z30ELsnTeL5pUvPe49rKOOfwcA7l8yqxYvRzGYqk5OZvWEDx0aO9I2Yn1BSgiksjANWK5O9B9E0bE4nCdXVTGlpofnhhylPS2Pq+vW8smsXPx4NpgEeqebNDrA4P5+y1FQmd2F6uxDJMT7t51GG80G3gx8NZXILukIZBzwD/K2X56sGZqBf55lAWYBttgFpQojxwGFgHvBPvTzfVHwuSAAAIABJREFURcuF2PUPJFNZD2TyzYn+5Zc8nJGhFMoAY7VaudJojJtDQqgYP953zxw2G3kbNzL8xAmOjh7NlN27GRcdzRv33OMbMe+dqTPd7ebE8OF4mpvJX7JEf5GYNYvg557D3jQWswksDgfWg4t0BZNn79d6KH9b38kxPu2DeM7uHPV/QDd9fYgetbW3q+3b7bsavU6xQogq9BH5C4GnhRBBQDOGL0QIkQjkSylvllKeFUIsBT4GzMAr/mn3LxV64osYSHah3zx7NzKlnYtMUqKmShh4/BvjuL17ySov900rIKQk1Gxmu+FTKZsyhe1WK2/4NdS+mToXLOC6rVs5PnIkSElOQQFjDh5kRHM9RQ2RTIhyUPDOKjaLKaSNPMEohwPRjY+tN3VR/raLi+56Kj8EGgEb8LBfgyAA2dU89VLKuzspujLAttXAzX6/PwA+6Ea2Ic2F2PX3l8lhJCvsybwkHo+H5sceg/p64uO7Skmo6C98jbHVyg/Xrm3zHGVER/O+xUJcbW3Aff19MjUWCzmbNpFVWEhrSAh1sbGM1o7wwLgXEUjqqkdzomUkZxuquX7dKkbGZ/LG9+cP2SzQiu7p8s5LKU1SyihjifZborpSKIq+433brM3NJbGXcf/9PXDMXyZTSAhfGinqu5LJ7XYzZ84czvzmN4QsXsxPf/rTPsmgOEfaPUf3zZgBwHRjjMo3jACQTZrmy5PmP1NnmdVKSEsLrcHBrFyyhNqEBDKMdPnOxjQaPFEsGbuSvNgNPGRZQUZtMePajbZXXFr0NKGk4jzQl66/f6LGsnPIGusNGZ7qchEZIHmgV6Zdhg+lq3Tyhw8f5tZbb6WwsJDQF19k2OLFDOtkW/s51k/Rc/yfI5PDQbGRgsU/8aPD4WCUyUSky4UWHw9+qe+jGhrYdvXVvp5L4bYsZsmNuM4kkBJW3mE+luyaGjVD5CWM6qMOUfwjtTbk5bGik5xN/viHDHuMucl7G8a8Z88epk2bhtPp5P3332fY4sV9qY6inwgYAJKSwnsFBb7w9Ze3bAEpKU1L47LDhzn9/7d35uFRlVn+/7w3G1kRAkkqECBJpbKwu9vIEgguiG27ICjjKDuoKPRuM/PrnsV2mZkeRBSjEXRoARVb29ZulwAVwF3ZE5JKJWwhlbAIZIGQpd7fH3WrqCoSyFIJRXg/z5MnVTf33jpVqXvP+55z3u8xGEi1WNCamrAmJXEgciAvH3iMirOxFNam0SQdt5EmqbGnfpijIVgbUXIs3Qc1U+mmtKV6zJmUf7WZkuH2lDF/+umnTJkyhaioKLZu3crw4cMvfpCiS2iu2GJkQQHlISEeiyUXv/ACaVVVHI+OJubAAUKbmljw8suE2evoV3uImOAKDtUNhMBwXiqdS1pkMXtq07EMSHO1gW417ZxVK/wT9R/rprRn4ZgvFINzcnK44447SExM5Ouvv1YOxc9oTmstLiCAPUaj6/+evG8ftaGhnOrZk14nT7J95Eiqw8MJO3OG/qfLeHzgS0yNf5fFiX8iIqCGIaOm8PnhkWT/32dszbcycMuWNs02otsxq1b4L8qpdFPaI9TYkRXMdrudJUuWMGfOHLKystiyZQv9+/f32ftR+IbmCkBuz8ry+L8byss51qePS/5+Q1YWKx59FCkEaeF7PXIoxpBdBAYG8uyzz3LTAw8wqm9fDtfXt0kBIkL1YulWqPBXN6U1C8e8k/JJSUkeJcOmVpYxnz17lkceeYR169Yxd+5cli9fTlBQUGe/RUU78S4AsdvtZGzf7io7vjo/H+rr+WHkSI8b/YGBA9m7z7OnfUn9MDLj4igpKSEsPt4jdLqolaHTmjasyRqn/zZ38DPw1Xlawh91+7oK5VS6MS1Vj40DsNt5yK2Nb6zZzJs//MDDU6bw99JStIoKKjMzmXGRi6Hhxx+55e672bx5M8888wy/+c1v1AJHP8Hcyv28ByAzbr2VzzdvpraoiE2Zma4bfe8ff2TfVcnklM0jObSYkvrhRMWmYzQa2bJli2ttC5wrACgrK8NoNGJdPwVbbQSGUUvOu8EeNxqZ5GdrsjpCeysvuwvKqXQzxum/zRfZL7qFpHxpaamrZHgEF46P2vftY/vtt9O0bx9r1qzhgQeaX+9qt9tJtFqJuQJHbZcL7gOQNBzh0/lvvsm87GyKTSZSSko42bMn70yZwtd/mUpFbTyZo+a6/pfuemNOJ5S8axevfPMN5QeLOV0ZS1JoMeYPs/khNoOp02ec+w5cYFY9DkchyYVK1/0Nf9Tt60qUU7nCGIfjIr3rQnFso5FEq5WBF3ACjd99x+nJkwloaGDD558zZsyYZl/P1fmxuhrrFThqu1zRNI0vR40iLT8fY1UVlWPH8qHJRLimYep9HFPv4x5rUYxGI/zwA3NycijRZxtXCUHe7t38ZEgfFiZn6yGzjeTYFmO1WjEemMtvT0SzNvYJrHobhebWZNXg+M76inFc3FGN03+b23F+f9Tt60rUVe3nmGn7F9tZItwSNcCuFpLyMTExZOjtf1uS3D/24YfUjh0LYWGM/PLLFh0KeHV+VJU9lwwzrf8eOQcCozdvpqpnT8SxY0zavp2IixxXOXIkJ/v0wVhVxeyxY/ndokU8+eSTZERZPJL7ycG7KN/yNC/vHI75wHWcqawke9MmMt5+29Wh8nLmSpfsVzOVK5T9XjpezqQ8ALoTaG7qvnz5cgoWLSLttttIePRRIjUNu93e4qzjSh+1XY44BwLeq+4HWa3sM5nOUyN2zyFYkpLoXVrKF9u3YzKZGD9+PBveL6RJbnIl94vrhlFfK9BkHeMjNtL34FH2xqYTfOoU0VYrmEzY7Xas66ewoDaCylFLqPCxUGVn4o+6fV3J5fFfUvgc6aXjValri1VWVp6fcE1O5vDhwzz22GM88cQTPP7LX5J5553Yg4Mvuur+Sh+1XY60NBDo20KJr3sOYUNWFtmzZrH3RCXW9VMwGo1cFT+UVw/O47OjE1haPIdT9WEEyjM8OvBlbunzOYsMS0mvLORY795EVFRgt9t5+61VmPfFEl1dzuwPs5n61irsdjs78G0oDByz+nE+PJ8vdPsuZ9RMxY8Yp/82d9Hruet4OZPyzSVch1qtvJSby3vvvcdvf/tb6hMTeamVq+6v9FGbv2K+wN9aaruwoYUWB806oZTBVJRtx6RpTJ0+A6t1NBaLhe/+/GeGp4WT2qvwPM2wikMGSkeMwGq1Ul1ZwKz+jjxMptzIC3oeBpOJGs5vy+BrOuq4rmTJfuVUFB54J1yHFhdzoqCAjz76iDfffJOBAwfydEMDUghMFgtxNhuHevXCZrM161RUo6XLD+dAwCko6QyNfm80Npvc9nZCAQ0NjNyzg2PBYVgsFoxGIyaTCZPJRHBwMCd25VB6JtkjJFZUk0ZD3yiOG43YtmwhKXi3h9MZErzLUUTSzHdsnP7b3FkfiKJNKKei8EDTNAqmTuWU1YoxL4+N2dlUVlZiNpu58cYbsVgsDN20iT5HjhBZU0NJUhIxR4+yq7qa0aNHN+ssruRR2+WIewvonm7rlVa3MBDw6LOTlMSoH74kRh4jKsB6XgnxoEGDKN2ZQFSAjZxDc0gOtVJYm05ARAI7Zjj6sBgMBszfDaXJbZFlQf0wJsXFUQM04fsQmMJ3KKfSzTE3s60JRwVYS0ghKPzgA3b89rdcPXIk3377LQkJCYDjBpKYl8fZEydY4Wwxm5l5RdXhXwk4W0DXmEzczIWTr+6z0eg9e+itnWJ+v5fPKyE2mUwYjUZ+iMvg1OEz9Ao8xncnRmCraeBXD00j6S9TOVUbgf2mp4iITXctstxVP5xyfZGlL9iB5/dfOSjfopyKwoPq6mpiZs6kbv16gu6/ny2rVhEWFub6u6ZpDElJ4a+NjaqiS+HCORsNsdnIKNtzXgmx87uhOXMs6/9ORW0Er0UPJm/hQuw1VaQPjCE13Mrmj14jIiadsXc9wTNHjrArLo6eukNJsViIsdmoaqbXT1tQs53OQzmVFhin/zZfQhu6moKCAu655x6sVis9/uu/CP7FLwhrRnIlPj7+vGR+S1pNiiuLIwYD+d8NZXwz+mBONE3DdP97mIBGYGJcHMnf/pX5A151zW5eK1+M7YutTNIku+KWOOTx332XMdXVlCQlccOmTbz88esM6VlB9M1LHHL7Kk/nFyin4seM03+b23hce1YgH1m3jutnzyYiIoINGzZwl952tjmaWz2tKrq6P+ZW7LPfaKQsNqNZfbCW6BsWxtBepa7ZjUDS1FBH4ZEIUsKLeOjDbI5FDqKAYHJmz0YKwSDbfuSPZ2gMKGf0h9mUxWZgd5d+uQDeod8LhYIVbUc5FT/DuRre3I5jncqoo202TrYyPGCvr+fMr37F3mXLGDVqFO+88w7x8fEXPMY9mT9MVXQp3JCaxvfTZzB3/d/P0wdrjh3AKYOBUrfEvKXWhATmDljhKileXvYEXybehD0gAJPFQnrlXh4d8JLr7y+65W0UlxblVLoJrlXNVVUE9OrF1d9/z8tRUcyfMYPAwOb/zQcPHmTnAw9Q/+WX9Fu0iE3PP996yXo9kdtTVXR1W0bQzrxDC/pgLbHfaCQyNoNXy+aRElrMtqqrGRm5zSMvkxGaT8Wh/qQWFnLDN98w1KvkeLBb3kZxaVFOpZtgtVopqKriZGQkvU6eZOfw4dQXFbHizTd5bIYjLODsA353fj7i0CHuXr2a08XFhK5bh3HqVFQHFIXP8JJyuRBST95PsY7m64oKHmlqovDbs4yXG115mb216QTKs2Rt2MCRPn0oLE9jQrRbyfGZwQxrarqgZBA4Bl8mi4W0/HyChaAgLY0UTSOuslIpaPsI5VS6CTabjbJevbjq5Mlzmk2ZmTyRne1QhDUaWbluHUnl5dSGhVE8eDA3PfII14WF8faUKZfafMUVjvtaprF2O6WWPbx08DHSwgvZXT+MA70HYG/QWDF/PlIIot6q5qWDj5Eetpe9pzM4HB7PQYuF0oqKFiVR7HY7K9etY1J5OafDwig2Gnnj448ZHxZGUWqqh4K2L5L+5g6f4fJEOZUL0JH8RlcxTv/9qsFA0vff893w4R6lvvkmk6sta9HRo1RFRbmczsbx45n/6qsM0uUvFAp/QNM0klOHsqIyhvcM91ERF4fh8GGC3MrYV06fyYxVqyj/sR8f3H0XFpMJIeUF10tZrVaKjx2jWr8GjCUlJO3ff55wplVdDx1CzfO6CUajkX5RUaQVFZ0n3hgbG8tHH33EschISrw0moqNRoapXuAKP8DMuQFcfHw8/U+dYuuoUVhMJjQpSXX7bgspiT5+nG0jrqYoLQ2paa5uk67e9rnjHD86NpuNo27XQJzNdt714OoppGg3yql0EzRNY/6MGQT06MG87Gwm5OayOCeH5B49+PWvf012djY9jx/3uDC1pibSioqoiYk573wjuLy67Sm6F0ajkYzISObl5DAhN5drfviBoIYGFi9dypR33mHOa6/RGBREhtXqOYgqzicuLg673Y7lx2jyDg3EYrFgt9sxGAz0ra4mWVfNrjAYSC4pIbWwkDF5eaQWFjLMalUK2h1Ehb+6EYGBgeycMQOb1YqhooIhISH8etFcTtU2MODZ/yYpPJyyqiqXUKCxpISwM2eILC8nQiUqFX6Eh/6YzUaPEIip3UdqeCEF+zIovSqJM01NJDY1Ob7PSUlcs2cHcfWVNG57inXf3E5NuWcL4ykPPExKnz4UlZczPzsbS3IyvU6cICs3l6LUVCZu2EDP+nqSkpIu9du/rFF3jw4yDt/2YugwmoYlJobNOTnMnDmTa0eksfS5aQybPJmhJhPbR4xgU2YmDcHBmMeOpSYsjIiioha7PCoUbWEEvstBOvXHKg0GousreHTAS9zS53MW9n8R43ErIZWV/NtvfoP51ClGFO5kUN1++gUc5G/FKZwo28Ws/tlk9clllmEpVZV7KS0tZea0aXw0aRI1YWEMKSykLiqKFQsWkDtxIi/Pn8/JsDBKS0t99A6uTJRT6Wac2LABhg2DtWv5j98/xV23jSIl4CijP8ym1LKHlNJSrMnJbNZbAAtgxdy55KpWvwo/Jc5mY7DXupT0iEJC4+KYuXAh8vXXiTtRSs/AE5xs7EVMcCWmsL3N6o/Z7XZ+8tVXhJ85w8nISPJTUjxyKruTkrDZbJfsvXYHlFNpAW8lU3+nqqqKxx57jF1ZWRAWhun99+kbXsdsfbS20LCUxuoDBGsac/Q49T0ff0yRyXReotJmsxFtsTAwL88Vj1Z0f8z4Z6VjhcFAfv1QmqTjdtUkNXbVD+P9O+4gJDGRIbNnkxB2jJqmSGYl5DCq15fsO5PksX9J/TBiYmJ4ZdUqQuvqyJ47l/2JiaRaLB45mdSiIux2OzX43/U/jrZFRZzr0rr6OlY5lW7A8Y8/ZvD8+Rw+fJiQRYsYOHMmN27bdl6jo+TgXXxpupOvDQYMFRVMvPpqzhYXe3b4s1rZFRREbFERxUlJHrX7KteiuBRYdT2x7IMLSAm3sKt+GHtj07GYTOwqKyOyqorSwh+5Lvw7AoQdY3gxP1RdQ/bB+aSEWShtcOiPARyuqqJQL7u3C4EEjxyjFKJbfM9dChvV1V1+HV/+n14n4WwG5M8cPXqUvdOns2fyZHr27MkXX3zB/ddfz7ivvuJ0aCh7a9POG63VxMYCjpJMg8FAekQE83JyyNKrxQYFBLC/qYnsWbNUSEzRphyJuQ37tgWnnlhqTA25Zyfw2uTZrJw+AyElw0pK2JuRQXnvfhTp33dNSO6Le4ez9lCqExaQ+dO5TJ0+g8rKSkoTEs5Vf8XH0xgYiHnsWFeOsUdAAMtiY0m2WBh1Gc/UrVYrBdXVl+Q6Vk6lA4zDdz0ZnFPVkXl5nLrIF1lKyZo1a8jIyODou+8y8A9/YNu2bURHRxNWXU3OrFmsv+8+dvYfxksHH+PzY1m8aFtMZGwasdu2Mc5sJqChgdc2bwagcuxY4oODmZuZyZCUlPP7jTdTu29GlRx3d8z4UTjMbueLU8OQaEzYuJEJGzbw2IoVpEVGUmwykT17Pj1CAsg+OJ/Pj2axvHQu+cXlvPHGG4SFhaHpHSX7nzhBdUQEc3JySDhwgJ5VVWTl5hJUX8/kvDzSoqJc10jQZVy8YrPZWnUddwad5lSEECuFEEeEEHvcto0QQnwthNghhPheCHF9C8fuF0Lsdu7XWTZeiK6cpbimqvrN/gb9i0wzX+QbSkvpc8cdTJ8+neTkZK7Zvp1Bv/89wcHB2Gw2rImJGEtKGL1lC1uvu5kvkkZhsWew5adzufq6m6GmhpxZs9igj14Kq6tB0zgwZgwmk4n4+HiG6iM5OLeAUtXuKy4lg7Zs4UhICH9atIjcrCwagoIIrK8nKS7OsfAxMJAZw3aSlVxByMh/5fb7nyB96DXk5uaSkZHBK6+8QlJSEumRkfQ8eZIfr7qKETt30j8qitwJE2jQB1U3X331edfI5ThTNxgMl+w67sycyhvAcuD/3LY9D/yblPIfQohJ+vNxLRyfKaU85itjnC9i9tUJfYj7VNUpFxGRk0O01cpxXS6irq6O5557ju+eeQYtKIilS5fy+OOPM0EfiQDExsYy8qOPSN63j5KkJMZu3kyvUydIGlDGuyYTlXl5FDczejFUVLhex73feLEPe6WYO3S0ojswTv9tbsMxwm4nxWoldtcutmVk0BQUhMVkwmIyEdTQwIiyMte+msChjqxXNqalpXH33Xczb948FixYwJo1a7j3n/+ZoMZGwk+fpj44mLDISCypqRSlp2MC8lq4Ri43BeTOuo5bQ6fNVKSUm4EfvTcDUfrjnkB5Z71+RxG6mmlXVE60NFWN0Keq//jHPxgyZAh/+MMf6POzn3FdYSFPPvkkQgiP6g4pJadDQ12jrJxZs6gJi2DJ8KcBx+glpZnRS41XV74ZU6dSmZnpCompJL3iUmC327nz7bcZYzZzNCKiWTWIhIQEIoAIcCgje6kjJyUl8dlnn7Fq1Spqa2spOHWKlx57jDdmzOClxx6jsLqaFLdZiKFyGUMsey77mfqlvI67uvprEfCpEOK/cTi0n7SwnwQ+E0JIIFtK+WpLJxRCzAXmAgwYMMAnRtrtdqa9/TaR1dWUtrdywqk51AoJcIPBwFCz2bMKq6SEoiFDyL/nHia9/z6pqank5ubyHxMmuGz0ru7oX19PUVqah3MqSk2lz5Ej1KeleXRstCYnM0IfvXzvNXpxV4y9fMZmisuRC+XlrFYrYdXVvDZrFkjJEy++yKMrVlCYmkpqURExDQ2MHj36oq8hhOCRRx6hb9++LA0OPm/wFldRgcU5U+91nIFHy867Ri7HrqaX6jru6uHnAmCxlDIBWAy83sJ+o6SUVwO3A48JIca0dEIp5atSymullNf27dvXJ0ZarVYi9YR3eyonmtMduhDeOkdzcnLo8eOP/OPWW/nx00955pln2LVrFxN0h+K00bu6o6KpiVQvLaSUkhKO6qMsZ8fGTZmZNLqNXlRvb4U/YrPZsOozeHtgIMsWLqS2Rw+G79zJmYgI7rvttjaNvFNSUhi2b5/n9ZGfT4UQgOO6tZ6IJi70KEUmEw1BQWqm3g66eqbyMPCk/vhdIKe5naSU5frvI0KI94Hrgc1dYiGOL3NL6qUXi6va7XbefmsV1V66Q1Mv0D/bOVX9uLiYhs8/Z9O6dRR/+SXRd99N8p/+xG8HDnTt66w2ay5ktj09neF79zIvOxuLycSIkhJOR0Wx32hk+LkXw2IyUaJmIQo/x2AwYDSb2aTP4IWURJ88SU1EBAf69ydnyxYyduxATJ2KbMVN35lncK1LKSzkVEEBxU89hXjqKdYOiKO2wnHd/vS7v3E4NB577GTPk7QhAnEh9LN0y1xjVzuVcmAsjs9yPFDsvYMQIhzQpJTV+uNbgH/vSiMNBgPJbl9mV1w1M/Oix1qtVqorC5jVP1vvSreRnFb0z/7222/JXbyY6q+/Rhs5kg0bNrC8Xz9q9u/n3rNneddL6NE7ZBbQ0MDVu3ZRHR7OsehoRuzcyaCoKNZPmeJxwZlRK14VlwdGoxFNdwIlyclcnZ9PTWgor86Z49H/ZJDVyr5WJNGdg7fnrVbiKirIu+UW/n7HHXzX0MDg7ds5ShwLk7LP9b0vW8iqzz7DtH27Y0YPWH+MxlYbgcFiUeKrLdBp9xchxFocDrmPEKIM+D0wB3hBCBEI1KHnQoQQ8UCOlHISEAu8LxxT0kBgjZTyk86yszmMRiPVbl/mtlRO2Gy2ZleytzTLOXjwIE899RRr1qyhT+9gHv1FKq888w2l69cTW1zMqaQkYs1mVnl1pPOu7rg6P58zoaG8MneuyxEuzskhurS00xoOmTvlrAqFA2e49pTVyrCKCm5ISODlqKjzIgh9Kypa5VSc53RWjwXa7Wx8+23G3Xgj9YGBDK5e63HdZoQW8MGwn9FUVITFYmH7d1+0KQJxpdJpTkVK+UALf7qmmX3LgUn641I4F625FGiaxrqpU0mxWrmxooK5mZnnjUrsdjvRVis32GxEGAzY9b8bDAbM3w2lSZ7rn11SP4zMuDiPqfOJEyd47rnneOGFFwBYsmQJ16Xlc7Ihihu3bqWgqors2bNb7EjnHHX93WolXr/gXmnmgjNUVJznVCI6/yNUKHyDrlTc02RiiMXiEQ5zRhA2tCKC4E4AjnVoRj0v6ewCOeTDfMbLja7rdnf9MGzx8exqbKSgoKBdEYgrEeVinbh1ibPrtfFxNhs1cXHNOpRVb79NX7OZpoYGx0xCX3VrNBqJjM0gp2weuceyyLEtJio23TXLqT3TxDPPPENiYiLPP/889913H3v37mVwaiLlhw3EV5dz+7d/JeBsFdKZQGxhNayzuuPAmDEMGTKk2cVONTExJCpxSEU3wGg0cjoy0iWIujgnh3Q9Z9geBrjlJYuNRg5EDnIoUBzNYqltEXtj0ylJTCS1YA8NDQ3NRyC+eNqju2RnsgPfKXh0Jiq87oXdbmflunXceuwY1ZGR9Ny5k5V9+vDnadNA0zBzrvLqtWZ6W5tMJqZOn4F1/d+pqI0nc9RcjEYj9fX1rPgH5OYHYd33BmPGjOHpp59m6NChWCwWqisLmKuPgsbLjbx08HFMFgtFaWkuB/FCZiY1ND/TcA+HWXRxvLTISE5v20bfmhpK3EqjW5vYVCj8CU3T+NvUqcRYrQx0iyD8op3f5YMGA5Pc8pK7TEOoOBKD5aiJAJrYH5vAwuXLCTx6lNWf5XLrT9Lo3+MglWdjiQ2ppLhuGGnhx8g7NPDCOZYWkvs7cOQHPLe2jL+pJreEciropYR6Aq7JbKb48GHO9uzJoYQEgktLsZaXE22xcDwtDbiwro7JZELTNEz3v4cJaGpqYu3atXy9dQODYlNZdHsJ1rpJXBU/hMGDB7vO5z0KSg/fy50ff0z/sjJMLawnccc9HNZYUcE3mZksstvZtnmzy/ltbGNiU6HwN6SeE6nsQPWiU4LJ6lUNdk1+PnUBGi/Nm4+xpITJn33EwLqDpMXuxXBdClX1oeQeyyI1vIjco1nUyh4UHwtTORYvrux3z7kSYPO+WBqqyyn89q/0aDzN6zNmuFalV4eFEVNQ4DqmNbo6TmcydOhQ/v3f/50kQwiLUnKY2HcDc/svo6pyr2vdi8FgoNSrX0Rp/TCOXn01WnAwla1cT+IMh20dM8bRMa+yslnJib5dICqnUHQUp8jqDXl5RHdC6DZcH4h9k5mJFhzMjFtvZWjvaGavXMmwXbsYUHuQRwe8xMS+uUwybKRncC0LBq5gYt9csvrmEq6dPK+7ZLNr2U7scPxcIVzxTsW9BDirTy7zEpbTSztB8r59gONGbElJ8TjGuVjRGdudp8d2jUYjjY2N/PnPf2bw4ME8+OCDaJrGE088wdCrSpqtCHOeLzI2g1fL5vH5sSwbdKdEAAAeD0lEQVRe0PMw+8eO5dsxYziuz37aSkuyLJFxcUphWOEXtJQncBdZdc9bCi/HYrfbSdQdT2tzhgFuj50ti78dM4a0tDRmTpvGN5mZhNbVkRZe6LpmK8/Gkhpe1OJz72va3T5LdSJ5R2/wbU7TLQfsb1zxTqW50NPgHvkYyh2yZFpTE+lFRRxJT3cdo7mNcBr1mcRD997L6tWrycjI4KGHHiI4OJh3332XXbt2MXHixPNmIiX1w1wzG03TmDp9BomJlZRHxrNa7//Q0ZXuRqMRdOeX5ZbYPH4ZSk4orizcFSM2ZGWRratajLda6YlD3sXpeG7QHU9rZeqdun5OR+TuqJxO5usbb2S32zUbG1KJxa0/UWxIJUW1qR7XdNHpDGJiYlznstvtrHtrJRsrb6S+wc7GD15h3Vsru32xzBWfU2muBLiwNp3hO3cS2NiIsaSE0Lo60CuxnDi/fNUDB2JcvZr022+ntLSUESNG8Je//IW77rrLNbswGo38oFeEJYcWU1I/3KMizHm+uN7HCe19nH0mU4e8fQTnNJUqR44kKD8fY1UVs8eOxWQysboVzsrcgddXKDpKS3lL9xL5ixXMNIfdbmeKruvnLF6584cf+JtXeLnYaGSv2zVrPTuc+qCeHs8bgnqSc2geSaEW9pxKIb/0CC/fey9Llixh6tSpFBcXc+LwbuYPeNVVgPPKocexWCyg52e7I1fmTMVt6uhdApxd9gSHwvvz6a230hAczKbMTLYPG0ZNZaXHNP3kyZPUPfss9kGDsMyZw1VXXcUHH3zAtm3buPvuuz3CVc6ZSGZiJcGR8a5OdN4hrUVZZiZ0UP7BhTN8sHkzx3r2RBw7xhfbt/vm3ApFJ9NS3tJdUbs9jajcdf2c/VLCqqoY5JULkZrGSrdrdvxdc1m4+KnznydVEhLVj9vvf5Ks236KpmlMnz6d9PR0/va3v2EKLfCIgphCCyhwy8+2FrvdTorFwui8PAoLCyk83nptwa7mypip5I5jKY6btjfOG76zBDgmfhQFFRUUp6S4ynlv+ewzigcPpga4qayMUUuXkp2dzdmaGpg4kWGrV/P9hAkIr9mM9+s4K8K6guhmerR4L6BUKPwV7xJ5ZwXk80Yj4/V9WlL3vpCcUnO6flZn8Yp+XdTa7Zj0dWqMeIab3UqFTb2PY+p93LWv+zWdlpbG/fffzwcffMB//Md/8NVXXxF+XYrHgsri0ybaKqLvDPONqa6mJDGRdz54h2gZT2p4oV9WnV0ZTuUiaJrm+rKYA/tAZeU50bmSEqQQNB07hv2RR/jmrbf4Tkruv/9+/varX3Fm5Eh6AS27kzagL7ocYLNhMRjAaGx3XiXiQqM45VQUfo57ibxWUeGogPRaB+ItENkamfrmdP2MJSV8ozui5sJjHm0vLhJJ0DSNe+65h7vvvps1a9ZQVrKDnENzSA6zUnLaSG1TBBkZGW36LKxeK//Dd37CfMPLfruy3z9cWxdyMVn6+Ph4QgIDMY8dS0NwMOaxY5H19Rz8l39Bvvsu8QsWUFxczJo1awgYOdKndmWsW8fETz4hvrSU1z/5hIx1686rdvHGjGf+wwycBJZcwnaiCkVruFjllrNEvqUKSM2rJLg1MvVGo5Fqr1X5p91W5TcXHmtPO2EhBA888AAJg1I50xjCwTMJnKwLwVJSxjPPPEPN7t0XPoFbiN49zBdnszG0BW1Bf6HbOxWnEyk6NJDowkLWvbXStSbF/GE2b7+1yvFl1rvGGY1GTGFhjMvNJaiujnEff8xpq5VxNzfy19UjMC5bRmJiIqC3OvWRBIrFYiG8vJyzPXpwKCGBkz16EF5eziCLpV3n8+7R4l72bEYl4hWXlvZWbnnjXhJsamXp/baRIznSpw/RevHK39wUJi7U9qLV6A5B0zSm/dNMJsVvxhh5mJ9MuI+wqD689957VN93HyNnziRo9WoaGho8Dvce+MbGxroGiBUGg0dVmnclqT/QrcNf3r1NbvxgFSeoY35C86JwBQUFvPjii6xevZp+/foRd8MN3JCZyfSFC/l2c75rJa7z3He+/TahHekO6UZ+fj61YWHkuOVA5mdnYyoooNJZKeIVHruQ9HZrwgcKxaWiPZVbHcXpyMbqoS1jaamjeMXt9TrS9qI5NE3DFLkPU+Q+GD2aUaNGkX7NNeyoqsJqMpG0ezfzFi8mbcAAZs2aRa9evc7rxxQRm056RIQjzJeURG1jFCsOPkpaeGGzlaSXmm7tVLx7mwQeb6DBHnTe1PHzzz9n/vz5bNq0iZCQEKZPn87KhQuxjBhBHo4v47EfozlVG+FY2Ws0erQ6dZdAae9FIYSgOCXFY4RkSUmhd3W147ndTsbbb5PaUqy3GZzhgx0mEyO4AqalisuGi0kddQbuuQmnI4vwki0yGo2caWOexh13ySeXHtiUkx42FNfV8dqCBQ4bJkxg3osv8vqKFfz+979n5syZDBkY7NIBdA58xybso1cv2BvyJPf/bCrsXMKR0+e0Bf1psOg/lnQC3gsb40NsFJ9O8Zg67jqRxLJly9i3bx9//OMfKSsr4/XXX4cRjpUeztmObV8s8dXljNJDZocPH+74NNmNjIwMUi0WjxxIalERVVddxcC8PPLy8sAHsd7mMKPCYYqupTVSR77G6cikEJgsFm7eupVDvXrR12Zz7eMUrdyUmUlDK/M0TrwlnzzC6142eFSfDRnCf/7nf/Lwww9TVFREcjM5kyOnI4jufZzN+sr/tOjjjEk40OqQX1fSrWcq3gsbE8NKqD1yF8tL55EeaWH3SSPHT4ewbNkyJk6c2Ow/xznbcVcQzrEtxmazkhrWj42ZmT6ZJptMJuzbtjH31VexGo0MLS4moKGBgVYrJUlJfLZzJ8VDh2IsKSHOZqPCYGBXUlKnjuwUis6iPZVbHcVgMDB00yb6HDlCpK7cHXP0KD2rq/lu9GhXpaVTtNLSRtHK1nR9bakMenhmJlOmTGHbtm18+ekbTJCbXGXI+VWpxNQfp56BmPRIiX+5EU+6tVNxrGRP59WD8zCGWdh10kjRvjKKi3Zzy03JTJw2k2uvvfaCnr6lTo4bzo7nx969XRdFemEhA0ND231RaJpGwbRpVFqtJFRUcGd6Ou8WF7vCa2Xx8dz58cck7d9PSVIS48xmelVVETN58sVPrlD4Gc6c33Krlb4tNMJrDXa7nUSrlZhW5BmNRiOJeXmcPXGCFfPnO8JPmZnM0UNgHS21b03XV6czfSL7RfJTBmMsLfVwpiNGjKAof7hrpX5+VRo/ngmksTGSoUElzH1vOW/FDWb6wxv9bobixD+t8gHl5eX8z//8D398/gWWvrGBP37QwLoD9Uz7p5ls+rfeLJl0kuuvv/6i/5jmFIQtZwZTOiCRddOmuabJR2JiGNbB2KZZ07DplSyapnkqDGsap0NDXeGvnFmzOBMa2u7XUiguNe2p3HKnrRVkmqYxJCWFwtTU5hc/dpDm7hXelVmapvHwlClcH7KDyfkfcdpk4uHobLSN411/nzp9BplJlfwY1Y8dN99Lv16CRSk53NJ3A4sGvMSRA9tZuHAheXl5freaHrqZU7Hb7axdu5bbbruNhIQEfv3rXxMREcHi2yR/WNhI+MqVZGVloWmtX6rolHHJdlMQDu6ZyKDyA4zZvBmAL266iYQTJzAYDD57L94Kw4aKCoq8Loa9qakcOXLEZ6+pUFxOuFeQtTbPGBcXR4rV6rquAhoaGJyfT+9jx85bFhDQ0kla4GJdX8Fxj3p37ZucPBZGulbAiB8+4t3CDOzy3Hmci7FTEw4gNY3kkD0es59hvUqwWCyMGzeOpKQklixZQmFhYRut7Ty6lVPZuXMnDz74IHv37uV3v/sdRUVFfPXVV8y/M56eUUHndtTXpLQG58gh3qkgPHk2YWFhDKg7xJKdT7Pg/ZdZ8j9/JC0y0qfxYG+F4Rt27jwvka8WMyquZNqq/WW329m6bRvhNTXMz84m67PP+PnSpTQFBnI8Kso107nYguOWaI3Gn3uOdmKfXBYallJ1Bqwnoj1PlmVmUZaZmmZmP9azQ/nv//5v3nrrLdLT03n22WdJT0/nuuuuY9myZZd8oNmtciq9evVi/fr1jNHDRy6yzCzqwHk1TaNP7+Nc1fs4aBo1R/Yyf8DLrsT98sNPcvPVV/s0xqlpGgVTp3LKamVYRQUzbr2VZdu2dWliU6HwZ9qq/WW1WtlbU8MLTz5J8r59DNm9mzOhoa78intn1B3tzK+4Sz41l6NpNu8SWkxFbbxnUYBTr3D8RiJi01lx8FFSwwrJrxtMlbyK7/LzmTltGg8++CAVFRWsXbuW1atX8+STT/Lzn/+cW2+9lenTp9N4550ERka26720l27lVAYNGsS4ceM6fB67m6CcM/k3OcvMKeCWvLzzvhQZPfZw5MgR0nwtZ63HnHuaTKQBBSaTK5Hf3sSmQtFdaGsFmXNm0xQUhMVkIs5moyoqqvnOqB1J2l8gCuJdkdpgDyC/ZjAJAWexNNfnXtO4+rqbWXXsBOuH34stPp6SxEQWrVzpqiqLi4tj8eLFLF68mPz8fP785z/z1ltvMX36dERICL1vv52199/P5MmTiewCB9Ot7kh1dXUdTlzZ7XZWrlvHrZ98QpKuwfX62rUkFRYyJi8PYbefNx3N7yqZBE2juB2JzRGodSiK7kdbtb+818ZUxsaS1kxI+WgnXsvuXV4/O5rFnw78kkDRSFST1bWupbGx8Zy0lMVCRUUF24eMIG/cOCwmE01BQS2G+QYPHswzzzzD/v372bp1K/Hz5lH9zTc8+OCDxMTEcM8997Bu3Tpqamo67T12q5lK9aljvP3Wqg7JQFssForLyzkbFcWhhASCS0ooPXSIG6qrKTYaGVRSQjW9ySmbR1JoMbvqh2PzM5kEheJKwVlBtq8Va0q8ZzZGq5WY+nrm5eRQnJzM8JIS0iIj2Wi3MyYvjyMGg2NNiI/D2lOnz+Cz9X+n9JiRHkH1zO+3wrWu5bXyxbz5+gqaTjpkWkZ/mE1J5CCGNjayYfz4Vq+J0zSNUaNGYRw1iuT//V+e/vJL3nnnHdavX8/7779Pjx49mDRpEvfdd5/P3pvrtX1+xktI36BKqir3dmiVeX5+PtW6BteGrCzyxo6lOjycTePG0RAUhHnsWGzhvTH1OU1FZDw5P53L953Uy8AMqpe8QuEjnDObPH0ZwDfjx/PUwoVUZmYSHxzM7LFjAbh+82aCGhoY006Ry9bY8ex972KOm0hGD8/Krt7iEPUnS5jVP5ssPZHfWH2AQUHhPL58OY+sWsUvly9vU2GQ0DRuvvlmli1bRllZGZs3b2b27Nl8+eWXPPjggz59b9DNnArgExlodw2u2IoKpKYxVv+ijd28mTq7HU1IEhMOUGwydbiXfGdiRoW+FJc33uFbM+3/Tmt6CHnzmDHsM5kIDAzkuMnEAb24Z29NjatEOceHUkjuOHX8YsvKKDiT4RFKP1SfTIpXt8jk4F001NdTHxJCWf/+NISEtPq1zHh+VpqmMXr0aF588UUOHz7MF1984bP35XoNn5/xEtNRGehevXqRVlR0ro7dbkdIyeszZnBowACKTCaCz56lye6TtlwKheIS4R0FaE974vZgtVqhupqXFixgd/xQltoW8fnRLF49/CRRvePPy9kW1w2h0m7n1TlzyJ04kRfmzKGwurrDzk7TNH7yk5/44i150K1yKkcbYjssA61pGlIIlwbX8B072Dl8OPevX+/qBlfXowfb6q/n6fErfWi9QqG4lLSnPXF7sNlsFOtVaCunzyDFamXU1q1MTriKn/70p7y79k1yyuaRHFrMnvrhBIfGsT0tvUsVnTtCt5qpRPbs0+Fezc7Oj5syM6kPCWHbNdcwND/foxtc9rx5lAvh0AvqIux2O9EWC6Pz8kjsYEMwheJywUzXhW+dje2cXSHnuDW28yXuahlS07AmJ9O7roohQ4YQGBjoWkB5PDKeLT+dy/iJkxm6b99ls/C5W81UevTo0eGEubNCpDYvz1H7brUSdfYsuwcP9lAI3mM0dryevbXoGkex1dWcSkriBrOZVR1oCKZQKM7HXeQyqqKCbzIzebkT1oIZjUbQq9BKkpMZbtlDuixzOS9N0zDd/x5z9f1NdjsZ27dfNgufu5VT8QWapvHQvffyjw8/ZHB+Pjf260fqqFFUf/LJeQrBn3aRQnC0rnGU7dVcqDO75CkUlwu+rJB0liifMpnoSeeEcrzVMuZmTWpVF9fnrVbiLoOFz8qpeNHY2MhzL71EXHAwBRkZNFks5JeWcjo83KNj3OOvvEJXpeojLkGXPIVC0Ym4qWW05grW2tnj5VLgn67uErJlyxYqg4NZMX8+uRMnsmL+fGoBi8l0nkJwny4Sbqu5BF3yFApF52Gm+65BU07Fi4MHD1Lk5UD2DxzYrEJwZ8o5uHNcTyDO6+QEokJxpRBA8wUAdrudRIuFMX5YEBNA2+X4LwXKqXgxYMCA8xxI7JEjIISrKmSxflPf3wk39RE0M4LRY6qVmZk0BgfzTRv6ZisUitbh3vQrqKGBmzppRX13p9NyKkKIlcBk4IiUcoi+bTjwChAB7AemSymrmjn2NuAFHI45R0r5bGfZ6c3o0aPZumMHC1asoCg1lbSiImIaG1k6bx599u/H4JYoW61p7OgiuzRN47jJxA6TiRGo0YBC0R7sdjspXgrkTkUM96ZfztzpYlUQ02Y68970BnCb17Yc4LdSyqHA+8CvvA8SQgQALwG3AxnAA0KIjE6004PAwECeWriQ/MGDibfZuGfwYJ5auJCm4GAsJhNb29n6VKFQXFqcM5ExZjMBbu2H0WciXbWivrvTaTMVKeVmIcQgr82pwGb98efAp8C/eu1zPWCVUpYCCCHWAXcBBZ1lqzeBgYGY9VW0b3bViyoUik7FORPJ8ZqJRFutHDeZumxFfXdHSCkvvld7T+5wKh+5hb++BJ6TUv5VCPFz4N+klJFex9wH3CalnK0/fwi4QUr5eAuvMRdc64SGAHs64734kD7AsUttRCtQdvoWZadvabOdUVFRBhkbG18dFeXaFllVhaisLK+qqrIB9O7b1xQQEhJeFxqq9Thzxt509mztj0ePWrrSzktAqvd9uCN09TqVmcAyIcT/Az4E6pvZp7nlHy16Pinlq8CrAEKI76WU1/rC0M7icrARlJ2+RtnpW5SdvkMI8b0vz9elTkVKWQjcAiCEMAF3NLNbGZDg9rw/UN751ikUCoWio3RptlkIEaP/1oB/wVEJ5s13QIoQIlEIEQxMwzGrUSgUCoWf02lORQixFvgKSBVClAkhZuGo5LIAhThmH6v0feOFEH8HkFI2Ao/jSOLvBd6RUua38mVf9fHb6AwuBxtB2elrlJ2+RdnpO3xqY6cm6hUKhUJxZaEWWygUCoXCZyinolAoFAqf4ddORQiRIITYJITYK4TIF0I8qW8fIYT4WgixQwjxvRDien17tL5/jRBiude5rhFC7BZCWIUQy4QQPlGub4eNE4UQP+i2/CCEGN/ZNrbTzuv1bTuEEDuFEHf7o51uxw3Q/++/9Ec7hRCDhBBn3D7TV9zO5Td26n8bJoT4St9/txCih7/ZKYSY7vZZ7hBC2IUQI/zQziAhxJu6PXuFEE+5nctf7knBQohVui07hRDjOmSjlNJvfwADcLX+OBKw4JBu+Qy4Xd8+CTDrj8OBm4H5wHKvc30L3IRjHcw/nMdfAhtHAvH64yHA4c62sZ12hgGBbscecXvuN3a6Hfce8C7wSz/9PAcBe1o4lz/ZGQjsAobrz6OBAH+z0+vYoUCpn36eDwLr3K6p/cCgzrSzHTY+BqzSH8cAPwBae23065mKlNImpdymP67GUQ3WD8diSOey2J7o61iklLVSyq1Anft5hBAGIEpK+ZV0fFL/B/zsEtm4XUrpXHeTD/QQQoR0po3ttPO0dFTiAfTQ9+vUz7I9duo2/QwoxfF5Orf5nZ3N4Yd23gLsklLu1I85LqVs8kM73XkAWAt++XlKIFwIEQiE4ljwXeVP9yQcDmeDvv8R4CRwbbtt9JUH7+wfHCO9g/qHkq4/PgQcBgZ67fsIbjMV4Fog1+35aBzyMZfMRn3/+5x2dZWNbbETuAHHjboGuNsf7cQxO/0Kh/L1H9BnKn5o5yCgFtgO5AGj/dTORcBqHCX924Bf+6OdXvuXAEP80U4gCFgHHNX//3O70s5W2jgXxyw/EEjE4VTuba+Nfj1TcSKEiMAR3lgkHVL5C4DFUsoEYDHw+sVO0cw2n9ZSt9VGIcRg4DlgXlfZ2FY7pZTfSCkHA9cBT+mxdX+z89+A/5VS1nifws/stAEDpJQjgZ8Da4QQUX5oZyCOEPJ0/ffdQogJfminc/8bgNNSSqfmn7/ZeT3QBMTjuGH/QgiR1BV2tsHGlTiUTL4HlgJfAo3ttrEzPLiPPW0QjlHTz922neLcGhsBVHkd8wieMxUDUOj2/AEg+1LZiEN6xgKM6iob2/tZuu23CcfIxa/sBLbgiFPvxzHC+hHH4lm/srOZY81++nlOA95w2+9fcbSo8Cs73f7+v8Dv3J77lZ042ng85LbfSuD+zrazg9/NL3GExNplo1/PVPRKg9eBvVLKP7n9qRwYqz8eDxRf6DxSShtQLYS4UT/nPwN/vRQ2CiGuAj4GnpJSftEVNrbTzkQ9DowQYiCOtgX7/c1OKeVoKeUgKeUgHKOsP0opl/ubnUKIvsLRKwh9pJqCI7nsV3biuBENE0KE6f//sUCBH9rplHuagiO8BPjfdYQj3DReOAgHbsRxo/ane1KYbhtCiIlAo5Sy/f9zX3pwX//gmH5LHNUoO/SfSfr2H4CdwDfANW7H7McxWq3BMaXL0Ldfi0MWvwRYju6xu9pGHJpntW777gBiOtPGdtr5EI58yg4csfWfuZ3Lb+z0OvYPeFZ/+Y2dOGLU+fr2bcCd/minfsw/6bbuAZ73YzvHAV83cy6/sRNHru9d/fMsAH7V2Xa2w8ZBQBGOhH4unnnVNtuoZFoUCoVC4TP8OvylUCgUissL5VQUCoVC4TOUU1EoFAqFz1BORaFQKBQ+QzkVhUKhUPgM5VQU3QLhUKh2KtZWCCEO649rhBAvd7EtbwkhioQQe4QQK4UQQV5/v04I0SSEuM9t235dDXaHEOJ7t+29hRCfCyGK9d+93P72lK4eWySEuNVte7PKsrrG3Nv69m+EEIM683NQXJkop6LoFkiH8OEIKeUI4BUc0i0jpJQRUspHu9ict4A0HOq5ocBs5x/0BZDP4Vhk6E2mbvO1btt+C2yQUqbgEP37rX6eDByr3wcDtwEvOxdXAitw6Dml6D+36dtnASeklEYcK9Gf6/hbVSg8UU5F0a0RQowTQnykP/6DcPS2+EyfGdwjhHheH9V/4pxR6CP9POHod/OpcKi1thop5d+lDg7p8P5uf16IQ4/pSCtPdxfwpv74Tc6pxN6FQ1L9rJRyH2AFrhcXVpZ1P9d6YIJzFqNQ+ArlVBRXGsnAHThusH8GNkkphwJngDt0x/IicJ+U8hocWk1Pt+eF9HM9BHyiP+8H3I1jJuWNBD7THdlct+2x0iGXgf47Rt/eD4farJMyfVs//bH3do9jpKOtwSkc/VIUCp8ReKkNUCi6mH9IKRuEELuBAPQbPrAbh1xFKo7maZ/rg/gAHArD7eFlYLOUcov+fCnwG+noT+K97ygpZbkQIkZ/7UIp5eYLnLslBdkLKct2iYKv4spGORXFlcZZACmlXQjRIM/pFNlxXA8CyJdS3tTSCfTcxQ/60w+llP+vmX1+D/TlXGsDcOgordMdSh9gkhCiUUr5gdQbt0kpjwgh3schmb4ZqBRCGKSUNj205QyblQEJbufuj0MwsAzPcJtzu/sxZbpYZE8cOnkKhc9Q4S+FwpMioK8Q4iZw9Rgf7L6DlLLJWRTQgkOZDdwKPCCltLsdlyjPKSqvBx6VUn4ghAgXQkTqx4bj6L7o7A/yIfCw/vhhzqnEfghM0yu6EnEk5L+VF1aWdT/XfcBGN6eqUPgENVNRKNyQUtbrpb7LhBA9cVwjS3FrVdwKXgEOAF/ps5K/SCn//QL7xwLv6/sGAmuklM6w3LPAO0KIWThk1KfoduYLId7BoXzbCDwmpWzSj1kAvIGj8uwf+g845NBXCyGsOGYo09rwnhSKVqFUihUKhULhM1T4S6FQKBQ+QzkVhUKhUPgM5VQUCoVC4TOUU1EoFAqFz1BORaFQKBQ+QzkVhUKhUPgM5VQUCoVC4TP+Pzy+NuDT0qgCAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Now suppose you wanted to combine the two together:\n", "my_event = MulensModel.Event(\n", " datasets=[MOA_data, OGLE_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", " data_ref=1)\n", "my_event.plot_data(\n", " subtract_2450000=True, data_ref=1, s=5, markeredgecolor='gray')\n", "# MulensModel automatically fits for the source and blend flux for the \n", "# given model.\n", "\n", "# Customize the output\n", "plt.legend(loc='best')\n", "plt.title('OGLE-2003-BLG-235/MOA-2003-BLG-53')\n", "plt.ylim(19., 16.5)\n", "plt.xlim(2810,2890)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Chi2 of the fit: 1774.42\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": {}, "outputs": [], "source": [] } ], "metadata": { "interpreter": { "hash": "4893225011d183c1b6f4cb363cb51655b562391822e07fee7770eae97b6f9d53" }, "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.9.0" } }, "nbformat": 4, "nbformat_minor": 2 }