From 0c7cdf5db1e544de3912ba44fb3d8bc8cbeb40a9 Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Sat, 28 Sep 2024 00:41:15 +0200 Subject: [PATCH] trajopt mga --- doc/notebooks/udp_mga.ipynb | 154 +++++++++++ kep3_devel.yml | 4 +- pykep/docstrings.cpp | 2 +- pykep/trajopt/CMakeLists.txt | 2 +- pykep/trajopt/__init__.py | 7 +- pykep/trajopt/_mga.py | 484 +++++++++++++++++++++++++++++++++++ src/core_astro/ic2eq2ic.cpp | 2 +- 7 files changed, 648 insertions(+), 7 deletions(-) create mode 100644 doc/notebooks/udp_mga.ipynb create mode 100644 pykep/trajopt/_mga.py diff --git a/doc/notebooks/udp_mga.ipynb b/doc/notebooks/udp_mga.ipynb new file mode 100644 index 00000000..1168e115 --- /dev/null +++ b/doc/notebooks/udp_mga.ipynb @@ -0,0 +1,154 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Multiple Gravity Assist (MGA)\n", + "\n", + "In this tutorial, we show the use of the {class}`pykep.trajopt.mga` to design an interplanetary trajectory through multiple planetary encounters, \n", + "allowing for propulsion manovures in the form of instantaneous $\\Delta V$ immediately after each fly-by.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "import pykep as pk\n", + "import pygmo as pg\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0\n" + ] + } + ], + "source": [ + "udp = pk.trajopt.mga()\n", + "prob = pg.problem(udp)\n", + "uda = pg.cmaes(1500, force_bounds=True, sigma0=0.5, ftol=1e-4)\n", + "#uda = pg.sade(4500)\n", + "algo = pg.algorithm(uda)\n", + "res = list()\n", + "for i in range(1):\n", + " pop = pg.population(prob, 100)\n", + " pop = algo.evolve(pop)\n", + " res.append(pop.champion_f)\n", + " print(i, flush=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAGVCAYAAADZmQcFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOt0lEQVR4nO3dd3yT1eLH8c+TpC2l7D3K3iBDAUVAVEQRFdwLt1733nvvdV0/13VdB9eJqIgMFZCtoGzZu0BZpYXuJnl+f5yWgqw2OUnT8n2/Xn1l9jwHhXxztuO6rouIiEiYPGVdARERqRgUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFihQBERESsUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFihQBERESsUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFihQBERESsUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFihQBERESsUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFihQBERESsUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFihQBERESsUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFihQBERESsUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFihQBERESsUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFihQBERESsUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFihQBERESsUKCIiYoUCRURErFCgiIiIFQoUERGxQoEiIiJWKFBERMQKBYqIiFjhK+sKiESc60J+PmRlQWbm/m+9XqhZc++fypXBcVi0OYXf1y1mc1YaW3PSScvdzo68Hews2EGWfwd5wZ3kBzPxO1m4TjY4Lk4wHod4PG4CXicBL5WI8yQQ51QiwZtIgrcSib7KNKmaTNf6bendtCOtajco6/9iIiFxXNd1y7oSImFzXdi4ERYv3vNnyRLYsAH8/pCLLvB62JHkZUdlhx1JXjbXiGNpk0q7ftbXicP1OPb+LIEkEmlE3UpNaF6tBV3qteWoJh3pUr8pHo86FSR2KVCk/CkogD/+gEmT4O+/i4Nj586D/25CAiQlQZUqkJSEW6UK2zwum8kjJ3cH8TszqJKVQ7VsP9WyA8QFDl5kbqV4NrZMZnvrFuR07IjbuStJ3XvgrV6d9JwsMnKzSM/NZEdeFjvzs8nMzyazIJvsgmxyCnLYWbCTLbnr2RlYj+tL2+913GACiW4y7at35/S2/RnS4SjifepkkNihQJHY57qwaBH8/DP88gtMnGi6qP7J64WWLaF9+z1/mjbdFSDExZFdkMewORP5cfk4VmX/getL3/uSwTgqBRvT3GlC17hkelVuTK8q9amakwurVsHcuTBvHixcaLrT9qV1azj5ZDj9dOjXD+LjD/pH3Za9k2lrFvHXxiUsSVvBhqzVpPtT8Hu34DjBPd8cqEy9uC4c3bA3F3c9kfZ1kw/+31IkghQoEpvWr4dffzUB8ssvpjtrd7VrQ//+cPjhxcHRqtV+P7RTd27nw7/G8Ova8WwumAPe3F2vucE4qtGWJlVa07luR/o27ULvpu1L9u2/oACWLi0OmKLbDRv2fF+1ajBoEAwZYm5r1izVf46svDx+T1nCuBW/83vqNLb6F+zxZwCICzShXdUenNL6OM4+rA+V4xJKdQ2RcClQJHZkZsLXX8OHH8KUKXu+VqkSHHMMnHgiDBgAXbvCQcYTFm5ax4ezf2R66iR2sAjH2a3/KlCF5IQeDGoxgEsPH0CNxCS7f5atW2HaNBg50vxs2lT8mtdrWixDhsDgwSYISym3IJ/hf09j1LKJLNkxk3zv2j3fEKxE04Te3ND9Ek5t1yPMP4xIyShQpGy5LkydakLkq6/MjCsAx4Hu3YsDpHdvEyoHEQwG+d+833h/7sdsdf/CcYr/env99WhbrRdntRvIWR17R2/8IRiEmTPhhx/Mz4IFe77esSNccglcfbVpeYVg2daNDJv3M1PXTyW1YC54s3a9lhBowaCmZ3FX3/OoXqlyOH8SkQNSoEjZWL8ePvkEPvoIli0rfr5NG7jySvMB27hxiYvLysvjxalfMnL1F+R71+16PiHQgsNrH8NFnQdxXMvDbP4JQrdypWm1/PCDmVhQNAMtMREuvRRuvRU6dAi5eH8gwGdzJvDJwi/YHJxV3DILVKZdlf7ccdRl9G7W3sIfRGRPChSJrunT4Zln4KefzDd3MIPl558PV1wBffqY1kkJrdiWylOTPmTW9h/Ba2Z5uUEfzRL6cedRV9G/VZdI/CnsSU+H776D11+H2bOLnz/5ZLjtNjjppFL99/inJVs28MLUT5i57Sdc3/Zdz1cNduTM1udyc68hVIo7+GQBkZJQoEh0/P47PPoojB1b/FzfvqY1cu65ZhZWKfy8bA6v/PEBa/On4HjMN3wnUJ2etQbzYL/LaVmrvs3aR57rwuTJ8Mor8P335jGY7rBbbzUttsTEkIvP9/t5549RfLX0S9JZsKsr0AlU5/gGF/DMCVeTlKBBfAmPAkUi648/4LHHYPRo89jrhcsvh7vvhnbtSl3c2GWzeWLqC+xwisch4gPNOL3FBdzZ+9yK8aG4cqVpsXzwQfH06Nq14dpr4cYboVGjsIqfmbKcl6d/zMKdv4DXlO/x1+WiNtdxV99ztHhSQqZAkciYNcsEyahR5rHXa8YHHnrIrBUppdVpm7ll7HOszPsFx3FxXYe6nu5c0/Vyzu98TMX8EMzIMJMVXn8dVq82zyUmwoMPwp13lmiSwoHszMvhsQn/ZdyGT3d1F1YKtOTOHndyQZd+YVZeDkUKFLFr6VK46y4z6Axmau8ll5ggad261MVlF+Rx37j/MGHTp+DNAaCO04PnjnuAo5q2sVnz2BUImG6wF1+EGTPMcy1bwmuvwWmnhV38pswM7hr7OrN3jMDxFABQm+482e8ejmnRMezy5dChQBE78vPNB96TT0JengmSiy82QdImtA/+92aO5a15L+H3pQLg8zfm1m53cXn3ATZrXn64LnzxhQnsooWTp5wCr74a8n/j3S3anMJdv7zAmvyJha1AD60SBvDSiXfRpk7DsMuXik+BIuGbMcOsoShaXzFwoPn2HMIYCcCMtUu4Z8LTbKdw1lOgMgMbXcEzA/6lvavA7Fn21FNmAL+gwOwOcNdd8MADZsZcmMavmMejk18g3ZkLmD3Ejqp5Lm+eertmhMkBKVAkdDt3mg+xN980357r1DHflocODWmq6/bsTK4b9RwLM3/E8QRwXQ9tK53MayffQ5MaoS34q9CWLIFbboFx48zj5GR4+WUzay6MqcZF/vvnL7wx95Vdq/ATAi14/YQXtIZF9kuBIqH54Qcz4yglxTy+7DLzYRbiSu+fl83hnkl37+reqhrsxNPHPsjxLTvbqnHF5LpmfOX224sH7o8/Ht54Azp1Crt4fyDAkxM/49s1b4I3BzcYz6mNruPZE6+qmBMhJCwKFCmdjAwzffXLL83jli3h3XfN9ighCAaD3Pvze4ze8I5ZTxKoxlXt7+aWXkP0gVUaOTnwwgvw3HOQm2u6wd54w3RFWmitzNm4mmvH3Em2ZykANTmC/w5+ofyt95GIUqBIyS1cCGedZWZyeb2m3/6RR8yJhiFIyUjjku/vYqs7E4Bqbmf+e9rLGgAOx+rVpuX400/m8WWXwVtvhfz/aHf+QIAbR73C1LTPzHYugarccNhDXH/kKWGXLRWDAkVK5ssvzar27Gxo0gS++QaOPDLk4r6ZP5Un/ngA15eG63roW/tS/u+U2/B5vRYrfYgKBs2MuwceMPe7dIHhw0Oatr0vo5bM4sEp9xMo7J5slXAy/z39Mfs7Nku5o0CRAysogHvuMYPtACecAJ9/DnXrhlRc8bfcT3GcIB5/bR7p9Qxnd+ptr85iTJgAF1wAmzeb81g++cQc9mVBek4Wl3//GCvyxgDg8zfgqb7Paqv8Q5wCRfYvNRXOO8/sMQVw//1mnUmIrYhlWzdy+Y937No2pY7Tk2FnvESjarVs1Vj+af168/9w2jTz+N57zZRjS9Ov35wxknf+fha8O3GDPi5udR/39TvfStlS/ihQZN+mTjXTTzduNN9uP/4Yzjgj5OJ+WvIn9069Bbw7cIM+Tml0Pc+d+C8NvEfDP1uZxx1nFkjWtzOgvmJbKpeNvIMMZz6u63BsnSt587TbrJQt5YsCRfY2fDhceKH5IOrUCb79Ftq2Dbm4T2eP54U594InF6+/AS/2e5ET23SzV18pma++gquuMhtONmxoTsfs08dK0fl+P2d/fTer838BoG2lU/nynKc1JnaIUaDInj75xJxLEgyaGV0ff1zqreV398rUEXyw9Akcj5/EQBu+Pft9kquri6vMLF4MZ58Nf/9tur3eestMLbYgGAxy7ciXmZH+CQB1nSMZcd4bOiXyEKJAkWJvvWWmnIKZ0fWf/4Q8XgLwwM8f8MP613GcIDXoxsjz3tFMoFiQmQnXXGMmV4D5/3799daKf3z8p3y99mUcJ0DlYBu+OfM97XRwiFCgiPHCC2bAFsx2Hq+8YjZ4DNG1P7zEtO0fA9DIewzfn//qobMPVMAPOzdCbjpUTwbHAzhmgeGu+x7wJVhZdBgS1zXjKi+9ZB6/+SbccIO14j+YNZZX5z+0q5vzo0H/4fBGLayVL7FJgXKoc12zOPGpp8zjBx80M7lC/KALBoOc/81DLM4x29e3TxzMl+c8VfEG33emwrYVsGM9ZKyDjPWF99fDpvklKyOuMtRoWvjTrPh+zWbmcWLNyAZOhENlzNK/uGfKLbjeDAhU47k+r2lacQWnQDmUua7ZA+q118zjZ5+F++4LubjcgnxO//I2NgTMNOPeNS/j3SF32ahp2XJd2LoM1k4v/tm+OvLXrdoImvct/qnV0n7AuK5pmb74onn8f/9X3O1pwZyNq7nip2vw+zbiBhN48IiXubDrsdbKl9iiQDmU3X578VTSMD9IsvLyOPmLf5HOHFzXw5DGt/DMiVfZqWdZ2LwIlv9aGCAzIHvrnq87HtOKqJ5c/FOtcfFttUaQUBXcoPnQdoOAW3zfDULWFkhfA+lrzc/23e5npu5dp2qN9wyYmi3sBEyEQyUlI42zv73a7AMWSOSlY95lYJvDrZUvsUOBcqh69VUTKGCOmb3iipCLCgaDDBx2PanBabhBH1e1fYTb+5xpp57RtHMTLPgG5n4Oqf/otvImQHIPaHo0NDsako+EStUiV5f8LEiZBaunmJ+UmRAs2PM99TrC4RdDl/MhqU5413Nd0zp94QXz2HKobM/OZOAXF5PjXYETqM5HAz+me+NW1sqX2KBAORR9+y2cc475EHnhBbj77rCKu/DrR1iQPQLX9XBDh6e54ajwj6WNmoIcWDwK5n4BK8aDGzDPe+Kg5XHQvA807Q2NuplB9LKSn21CZV8B4/FBu0HQ7WJoPQC8Ia6C/2eovPEG3HSTnfoD69K3MWT4UPy+DXj99Rh+xjBa1W5grXwpewqUQ8306dC/v9ni/PrrzUBsGN0md4x+m583vwXA4Ea3lZ9urk0LYcZbsPB7yN9Z/HzjHtD1AjjsbKgcw+tlcrbDguEwexhs+Kv4+Sr1Tf0PvwTqhHAssOuaLXaef948thwqC1LXMnTUxbi+7SQEmjH6vP9Rt0oEW3oSVQqUQ8ny5XD00bB1K5x2GowYEdaeTq9MHcEHyx7FcVyOqHoBH5/1oMXKRsimv+G35+Hv74qfq9EUulxguo7q2NmRN6o2/Q1zhplW1u5jPZ3OguMfLP2f6Z+h8u67Zt2KJRNXLuDmiVeBN5uqwU78PPRjkhLKsPUn1ihQDhVbt5owWb4cuneH334L6/zxb+ZP5bFZN+N4Ckj2HceoC1+L7anBmxfDb8/Bwu+Awr/yHU+HI6814yKxXPeS8ufDsnEw+1NYanYBxvFCt6Fw7L1Qo0nJy9o9VOLiYNIk6NXLWlWHL5zGo3/cjOPJp56nF2OHvqNtWioABcqhICfHbDs/fTo0awYzZkCD0Puup61ZzLW/Xg7eLKq5nfn5oo+oHBej3zC3LDEtkgXfsitIOgyB4+6D+uEfkRuzUufD+Kdh6Wjz2BsPPa6CY+6AKvVKVobrmg1Chw8359X/9VfIxxbsy7t//MQbf9+P4wRpnTCI4ec9F9tfSuSgFCiHguuuM90WNWqYbcw7dAi5qBXbUjnru6EEfVuIDzRhzHlfxGYfeG4G/PwI/PkxxUEy2HxTb3AInVO/7g/49QlYXXgEQVxl6HU99L3dTGs+mB07oGdPc0rniSfC6NFhbcfzT09M+Iyv15qutSOrX8wHZ9xrrWyJPgVKRfftt2YzQMeBsWPNh0KItmdncuIXQ8nzrsLx1+KrIcNoXzfZYmUtWToWRt4GOzeYx+1PM0HSsEuZVqvMuC6snGiCpWgAv0YzOONtM4vtYBYsgKOOMqd1PvwwPPGE1erd9ONr/LbtfQBu7vgC1/QcZLV8iR4FSkW2bh107Qrbt5stNooGWUM0aNjNpPgnQqAybx7/If1axFiXUXYajLkf5n1hHtdqCUPeMIsAxQTL4h9hzAOQsRZwoNcNcMLDEJd44N8dNgwuvtjcHzUKTrF7jvxp/7uNNQW/QqAq3wz+hnZ1G1ktX6JDgVJRBQJm3OS336BHD3NgVnzomzO+MOkrPl31JK7rcGfnl7mie+gtnYj4+3sYdadZfe54zAfl8Q9CvLZO30vuDhj7gBm8B6jTDs58Gxp3P/Dv3Xij2Zm4Zk0zntK8ubUqpedkcfz/zsLv20DVYCcmXTpMg/TlkAKlonr6aXjoIXOWyezZ0Dr06bDzUldz0U/ngTeHwyqfyefn2u3yCEvWNhh1uwkUMB+Op78JTXqWbb3Kg6Vj4YebIXOTmQ12zJ3Q727w7eeLR14e9OsHf/xhZgpOmQKVKlmrzvgV87hl0uU4noKKsw/cIUaBUhFNnw7HHGNaKR9/DJdeGnJR+X4//T69gCzPEuIDzZh8yfDYmdG16W/4/AKzH5bjNTOY+t1dtivay5vsNPjpbrPlDECDLnD+Z2bH431ZuxaOOAK2bYNrr4V33rFanXvG/ofRqW/guh4e7/k2Z3fqbbV8iSwFSkWTkQHdusHq1eYY32HDwloJf/X3LzIj/RPcYDzv9v+MPs1CnyFm1ZIxMPwqyM+Ems3h3I/N9igSmoUj4Mc7ICcNkurB0C/23wU2diwMGmTGZML8wvJPwWCQAZ9dzRb3Dzz+2ow+91saVYvhHQtkD5r0XdHceacJk+bN4e23wwqT7//+nenbPwNgSOMbYiNMXBemvmZaJvmZ0PwY+Nd4hUm4Op0J1081U6qzNsNHp8Kikft+78CB8Nhj5v5118G8edaq4fF4+OyMF3H8tQj6tnHxd/cQDAatlS+RpRZKRfLHH2Z6J8DkydA39NlN27J3csLnZxDwbaY2PRh/yQdlv+isIBd+vM3sBgzQ/Qo45UXwxpVptSqUvJ3wzZVmxT0OnPQUHH3j3l9MgkGzfc/o0dClC/z5Z1jb+PzTV/On8MSfN+I4QU5teAvPnWTn3HuJLLVQKopgsHgTv0svDStMAC797kECvs04gep8cvoLZR8mmZvh48EmTBwvnPISnPaKwsS2hKpwwedmVT0ujHsQfrrLHGu8O48HPvnEzPiaNw/ee89qNc7r3Jeja5ppyj+uf5sJK0t4CqaUKQVKRfHRRzBzJlStGvZ6kxcnf8Paggm4rsPtXR+laQ17222EZGcqfDgQUv6AStXh4uFw5NVldx57Ref1wakvw0lPAw7MfB++uBDyMvd8X506xYscH34Y0tKsVuPt0+6garAjjqeAuyc+qK6vckCBUhFs3158dO/jj4e1T1fqzu18suxlADomDS779SaZW+DjIZC20uwKfPUEaHV82dbpUOA40PsmOO8T8CWaLrAvLjTdjru77jro1MnM+ioaV7HE5/Xy3ikv4QbjyfOu4omJn1ktX+xToFQEjzxidhPu2DHssytuHfMyeDPx+uvx0ZBHLFUwRNlp8MnpsHWJOf72spFQW6f8RVXHIea/e3wVWDXJzKzbvfvL54PXXjP333oLFi60evlO9ZvQq9Z5AHy7+j9sydxhtXyxS4FS3s2bZ/4hA7z+utlqPEQzU5azMPNHAC7vcHPZnlGRkw6fngmbF5pDoy4baaYHS/Q16QkXfm6OQV78I/x4q5ltV+SEE+DMM826p9tu2/M1C14eeDMef21cbwa3jH7JatlilwKlPHNduPlmMyB/zjnmH3YY7hn/NI4nQJVgB27pNcRSJUOQtxOGnQMb50DlOnDpD2qZlLUW/eCcD822NrM/g58f3jM4XnrJbO3zyy/www9WL129UmUua3cLAPMzf2DG2iVWyxd7FCjl2dix5uCjxER4+eWwivp09ni2urNwXYfH+t5fdrO68rNg2HnmzPTEmnDp91CvfdnURfbU4TQY8n/m/rQ3YMorxa+1bGnWQAHccYc5Ytqi23qfUThAH+CeCc9YLVvsUaCUZ88U/sO6/npo2jTkYvyBAK/ONl0JzeP7M7DN4TZqV3quC9/fBGunQUI1uGQENDisbOoi+3b4RYWzv4BfH4dZHxW/9sAD0LAhrFwJr75q9bIej4en+j2I63rYzl+8N3Os1fLFDgVKeTV5svmJjzffCMPw6PiPyfeug2AlXjmpDA84+v0dWPgteHww9EtoVEbBJgfW+yazkSTAqDvMYD2YjUiLpqw/9RRs2GD1sv1bdaF1pZMAeGveS2QX5FktX8KnQCmvilonl18OjRuHXMyWzB38sO4DAPrUuYA2dRpaqFwI1kyDcQ+Z+wOfgWbaFDCm9X8Yug4FNwjDrzbTuwEuusicPZ+VZc6kt+yNk++DQGX8vlTuHfeu9fIlPAqU8uivv2DMGLNa+Z57wirqtrH/Bu8OPP46vHjSDZYqWEo7U+HryyHoh87nwpHXlE09pOQcB059yRwXkJkK311nJod4PMXTiD/5BH7/3eplm9SozcBGVwAwcdNnrEzbZLV8CY8CpTx69llze+GF0Cr02U+zN6xi7g5zjsjFbW+kasJBTu2LhECBCZPMTVCvIwx+TSvgy4v4JDj3v+CrBMt/gelvmOePPBIuu8zcf8T+WqZnBvyLOH8yeHN48Nc3rZcvoVOglDeLFsHw4eZ+0er4ED0x6W0cj5/KwTbc2ecsC5ULwc+PwtrpZhD+vE/Nh5SUH/U7wqDCcZNfn4B1M839Rx4xXwzGjYPly61eMt7n46J2/wJgwc4xWuwYQxQo5c3zz5vZUKefDoeFPgNqw440luWMB+DyTv8qm2nCi0bCjMJvmGe8DXVCP1VSytARl0Gns0yX5TdXQs52M4144EDz+rv2xzpu7nU6Hn9d8Obw5G8fWy9fQqNAKU+2bYP//c/cD3PA86nfPsbx5OHzN+DaHidbqFwp5WaYM+AB+txq1jhI+eQ4pquyZgvIWGuOFXZdM50d4MMPISfH6iXjfT5OaHQ2AL+lDiff7z/Ib0g0KFDKk88/h4ICOPzw4nNPQpBbkM+Uzd8BcFLyeWXTOvn1CTNuUrs1HP9g9K8vdlWqZlbSe+JMy3Pxj3DqqWZ9VFoafP219Us+fOxlEKhM0LeNV6d9a718KT0FSnnycWHTvmjAM0SvTh+B60uDQBL39xtqoWKllDILZpqpypz2is6ArygaHwF9bzP3xz4AwXy4pnDG3ttvW79czcpV6FJ9EABfLx9mvXwpPQVKebFwIcyaZXZ3HRpeCAxfbk487Fp9EDUSozwIHvDDyNsAF7peaPaIkoqj7+1mZ+j0tTDt/+Cqq8yGpTNmwOzZ1i/3cL9/4Qa95HpX8sW8SdbLl9JRoJQXRa2TU0+FuqEfePXN/KnkelfgBr081O8qS5Urhd/fhk3zzT5dJz0V/etLZMUnwYmFh25NfhkS/XBW4QzCCLRS2tdNpkm8OZ30ndkfWi9fSkeBUh74/fDpp+b+5ZeHVdRbs83eS43jetO+bnKYFSul9LVQtLHfiU9CUp3oXl+i47CzoWlv8OfAuIfhhsIFs8OGQUaG9cvdeZQ5b36r+5d2Ii5jCpTy4OefITUVateGU04JuZjZG1axOWjWCdzaswxaJ2Puh4Js82Fz+MXRv75Eh+OYtSmOx+zN1sRjTnXMzjar5y0b0Lor1dzOOI7Ls1P/Y718KTkFSnnw3/+a24suMptBhujZKe/jOEGSgu05pV13O3UrqfV/mZk/jgdO+7dWw1d0DbtA98vN/TH3wXXXmvtvv239AC6AKw8zE1VW5E4gJcPu2fZScgqUWJeZCd+b7VHCmd2VXZDHosxfALig3UU2alY6kwvPaznsHKjXIfrXl+g7/iGoVMOMmXWrDElJZqeH336zfqkrjjgRr78BjqeAd2faPeBLSk6BEusmToS8PGjRwqw/CdGwORPBmw2BJK7rGXq3WUhSF5jWCQ70uyu615ayk1S7eBrxvA+KZycWHVltkcfjoXMNM2NwYsov1suXklGgxLqxhQcJDRwYVjfR98tGA9C00pFUigu92ywkkwvPAe94OtRtF91rS9nqfjnEJcHmhTC4h3luxAjYvt36pS7rao6t3u4uYMMOdXuVBQVKrNs9UEKU7/ezJtdsIz6kdZS3WdmyBBZ+Z+73uzu615ayl1gTjrjE3E8fDe3bm1mLEej2GtC6a2G3V4B3Z460Xr4cnAIllq1aBcuWgdcL/fuHXMz/5k4EbyYEErnk8BPs1a8kJr8MuNDuVB3ne6jqdb2ZjLFiPPTqZp4bPz4ilzqshlmTMl7dXmVCgRLLxo0zt0cfDdWqhVzMiKWmu6txfE8qx0Vxm5NtK2B+4R5Ox6p1csiq2Rw6mO4o6hV2RUUoUC7pUtjtFZxP6k773WpyYAqUWGahu8sfCLAyZzoAp0W7u+uP98wRsa1P1Pnwh7reN5tbp/C8lIULYZP90xZPbNUVr79eYbfXj9bLlwNToMQqvx9+/dXcDyNQvlwwGbw7IVCJyw4fYKlyJRAoKG6dHHl19K4rsSm5BzTpBZUC0KqBeW7CBOuX8Xg8dKpxDAC/rvvZevlyYAqUWPX777BjB9SqBUccEXIx3yz6CYCGcd2je8Tv8l8heysk1YVWoY//SAXS+yZz2yjX3Eaq26uz6fZKU7dX1ClQYtXkyea2f38zKB8CfyDAimzT3TWoZeitnJDMNTsa0/lc8MZF99oSm9oMhErVITnPPI5QoJzUultht5ef/8waFZFryL4pUGLV3LnmtnvoW6R8+/d0XG86bjCBK7tHMVBytsMS0zKi64XRu67ENl+8me3XzAceB1asgDVrrF/G4/HQsXpht9dadXtFkwIlVs2ZY267dQu5iDErTCunrrcz1StVDr9OJbVwBATyoV4naNA5eteV2NfpDEhwoEnhbMMIjKMAnNHOjBemBZYQDAYjcg3ZmwIlFuXkwNKl5n7XriEXszR9AQAda4VeRkjmfmFuu16gTSBlTy2Pg4Rq0LTwQz5C3V6D2nbHDXrBm8X0ddrSPloUKLFowQIIBs1BWg0ahFREMBgkI7gcgOOb9bRZuwPbtgLW/W4WsnU+N3rXlfLBlwDtBkELn3k8fnxEdh+umpBIotsMgDHLfrdevuybAiUWFY2fdO0a8jf8qWsWgzcbN+hjUNseFit3EEsL18606AfVGkbvulJ+dDwDkr3gc2D9erMbRAQ0TTK7Ws/ePCci5cveFCixqGj8JIzurrErzLeyRLcZSQlRXB2/qnCPplZR3uJFyo9W/SGpKiQXfvxEqNurZ0OzmHZDrrq8okWBEouKWihhDMjP2WzKaFYlimePBPyweqq536Jf9K4r5UtcJWg7cM9urwg4re3RAOR71rMp0/7Rw7I3BUqscV2YN8/cD6OFsj5nMQBHNojilicb50D+TnOokmZ3yYE06wONCtdXLV4ckUsc1qApjr8GjuMyarHGUaJBgRJrtmwxK+QB2oV2dsimzAwKvBsAOLXd0bZqdnBF3V3N+4IntMWYcohI7gnVC8cH162L2GXqxLUFYErKnxG7hhRToMSa1FRzW7duyOfH/7h4Bo7j4vhr0ql+E4uVO4hVk8xti2Ojd00pn+p1hFpJ5n56ujnqOgI61DIt5aIp9BJZCpRYs3GjuQ1xujDAlHXm21jRt7Oo8OfB2hnmvsZP5GC8PmjZA4rmi6SkROQyRVPmM4LLtcAxChQosaaohdIw9Cm3yzIWAtCxVhTHMVJmgj8XkurpmF8pmeQeUK3wIyhCgbJryrw3m2cnfRGRa0gxBUqssdBC2REw/zh7JXexUaOSWf+XuW12tFbHS8kk94RqkR1H2X3K/JjVP0XkGlJMgRJrwmyh5Pv9BL3pAHRp0MJSpUpga+FWMXWjOE1ZyrfknsUtlNUrInedgDm2oXnVKHYBH6IUKLGmKFBCbKH8vSUFxwniuh7a10m2WLGD2Fq42rlOm+hdU8q3qvWhbg1zf+m8iF2mQ9UTAdhZsDNi1xBDgRJrirq8QmyhLNy8GgBvoCbxPp+lSpXANgWKhKBJY3O7bm3ELtGoSiMA0vI2R+waYihQYk2YLZTl20xfdCVPbVs1OrisbZC9zdyv3Tp615XyL7lwWvsG++fLF2lZ3bTUMwNbInYNMRQosSY93dzWrBnSr6/ZYQbka8TVt1ShEihqnVRvAvFJ0buulH/NW5nbzZE7qrdtHRNaBWyL2DXEUKDEmqK58iEe+7sp23SZ1UuM4k6/ReMnap1IabXqaG6z8iK2uLFLg+bmjjebbdkaR4kkBUqsCTNQthf2Ezep1shWjQ6uaIZXHc2ikVJq1DriixsbVasFwUoAzEtdHZFriKFAiTVFgeIJ7X9NdtD0E7euGcUtV9ILzwWv1TJ615SKoXpy8dThCO7p5QvWAmDx1shdQxQosScQMLchBEowGMTvSQOgY71mNmt1YLmFm1km1ojeNaViqNa4eHHjyqURu0yStw4Aq7ZHphUkhgIl1oTR5bVsWyqOx4/rOnSu39xuvQ4kr7BfOqFq9K4pFUN8ZUgq7PPaErnWQ82EegCk7NwQsWuIAiX2hNHltWq7mXLsBJOie0rjrkCpFr1rSsURZ8Y3yM2K2CUaVDZjiltzUiN2DVGgxJ4wuryyC/IAcNw4mzU6uLzCLi+1UCQU3sIur4A/YpeollAFgLxgTsSuIQqU2BNGl1eePx8AhyiukAd1eUl4ir48FX2ZioA4j/k3EXQjdw1RoMSeou1SCgpK/as5hS0UTzQDJRiA/ML1A+ryklAUBYo/ci0UX1GgoECJJAVKrKlduGXKttKv6s0NFLZQnCgGSt5uC8XUQtmD35/J6tVPMm1aEyZO9DJtWhNWr34Svz8yC/jKrV0tlMgFSrzXdAOrhRJZCpRYU8dMb2Tr1lL/ap7ftGqi2kIpyDa3jrd4cFXw+zOZM+dYVq9+jPz8FCBIfn4Kq1c/xpw5xypUduct/BgKRu7D3ucxXchBIhdaokCJPUWBsqX0G9nlFbZQvNFsoXgLZ5O5gYh+IJQ3KSmvkJk5B/jnsbNBMjPnkJLyShnUKkblFf692Rq5/bziPEUtFB0DHEkKlFgTRgslP1AGLZT4yrtVIHLTPsubDRv+w95hUiRY+PohLi0NTj4ZZhRO5f1itHm83X6wxHnNvwlXLZSIUqDEmjACJddvBuWj2kLxVQIKp30WaEpmkfz8Ay+gO9jrh4ShQ+GXX/Z87pdf4MILrV8q3ls0KK8WSiQpUGKNjRaKE8V1KI4DcYWtlAK1UIrExx94c86DvV7hLV0KY8fuPVU4EDDPL1tm9XIJhYPyrgblI0qBEmvq1jW3YQSKL5qBAsXdXvnZ0b1uDGvU6Br2/8/LU/j6IWzFQc6QX77c6uWKpg2ryyuyFCixJowWij9o/rF4nCj/b41LNLfq8tolOfl2qlTpxt7/xDxUqdKN5OTby6BWMaRVqwO/3tru2ToJvsIWirq8IkqBEmvqmU3sWL++1L9aNb5we4lAlFsKcYWnNKrLaxefrwrduv1G8+aPER+fDHiIj0+mefPH6NbtN3y+KmVdxbLVti0MHLj3jhBer3m+TRurl3OKxvlwrZYre4ryHh1yUO3bm9sVKyAvD0qxyWOdyjUAyA1G+VS6omN/87S2Ync+XxWaN3+Y5s0fLuuqxKbPP4ezBsPEqcXPDRhgnrdsU6Y51sHn6IjqSFILJdY0aAA1apg9vZYsKdWv1ksyhwgVuFH+YC86HTJDhxdJKdSsCR88BT0Lx/yGDoUxY8zzlqVmmS7kSh5tDxRJCpRY4zjQqZO5v3BhqX61UVWzbUvAiXLXU9FJjWkro3tdKf+ytkJCYXdUUXdvBGzNMS2UJG/1iF1DFCixKdRAqWYCxfVk44/gzq172RUoq6J3TakYMjexa7/G+PiIXWZ7rlksWTW+RsSuIQqU2BRioDStYWaIOY7Lxp3plit1AGqhSKi2LI5KoOzITwegZiX73WlSTIESi0IMlKoJibhB849y/Y7STzsOWa0W5jZ9TUR3jJUKaOM8CBTOvIqL3PqpLH86AHUSa0XsGqJAiU1FgbJiBeTmlupXPUGzyHD9ztJvfx+yqo3MJpFBvwbmpeQCBbB5UVRaKLlBc6po/aQ6EbuGKFBiU/36UKuWmem1eHGpftWHWd+QGs1A8XiKWynq9pKS2roUAnlAYcskgoGS75qp9I2qKlAiSYESixwHDjvM3J89u1S/muAxh1xtzU63XKmDKBpH2WZ3ywypwFLnm9u4wqm8EQyUoMdMpW9SvW7EriEKlNjVp4+5nTixVL9WyWtaKFuy0yxX6CAadjW36/6I7nWl/CoKlKJdAyI0hpKek4XjMWcFtaxVPyLXEEOBEqv69ze3v/4Kbsm3i6ib2BCAdTvXRqJW+9est7ldM61U9ZVD2Ma55tZTuHo9Qi2UVds3A+AGvdRL0jqUSFKgxKo+fcw/sPXrS7WVd4dabQFIzV0doYrtR+Me4ImDnRtge5SvLeWP6xa3UDyFR0dHKFBmbzD/frzBGng8+siLJP3XjVWJidC78Fv/r7+W+Nd6NOoIQLabQjAYxZ1V4ytDo8PN/bXTo3ddKZ82zoXcdLOxqFMYJBHq8pq10Uy/r+FrFpHypZgCJZadcIK5LUWgHNO8I67rgDebpduifCrgrm6vqQd+n8jSsea25XFQULh2KUItlBUZpoXSrKrdHYxlbwqUWFYUKBMmmCnEJVAjMQlvwEyNnLy6dAsjw7b7OIrIgSwrDJS2AyHfDJhHKlC25JktgbrW6xCR8qWYAiWW9egBVapAWhrMnVviX6vhawLAnE2LIlWzfWtyFOCYtSg7U6N7bSk/MrfA+r/M/TYnFS/ejUCgZOXlke/dCEC/5l2tly97UqDEsrg4OPZYc3/8+BL/WpMqZk3IyoyDHLNqW2INaFC4fmbFhOheW8qP5T8DLjToAlXqw8rCxbBNm1q/1OQ1C3GcAAQr0b1RS+vly54UKLGuaPrw6NEl/pWOdcxMry15ayJRowNrd6q5XTA8+teW8mHpGHPb9mRYtQpycsxBcgc7FjgE09eZmWSJbhPN8IoC/ReOdWecYW4nTIANJRtk75Vs9gLLdTZEdxt7gMPONrcrJ0BWFLd/kfIhUFDcem07EBYsMPc7dNj7OGALFm4zWxc1TGxhvWzZmwIl1rVsaaYPB4PwxRcl+pVeTdvhul4cTx5zNkb5jJK6baFBZ7NR5KLvo3ttiX1rpkLeDqhcBxodUbyjdtFWQ5atzzLdvh1qt49I+bInBUp5cPHF5vazz0r09spxCcQFzBYTk9fMj1St9u+wc8ztfHV7yT/89am57XCa2VS0qIUSgUAJBoNkumbHiF7Jna2XL3tToJQH550HPp/ZKPLvv0v0Kw0qmTn3U1J+j2TN9q2o22vNVNgR5bUwErsyt8Dfha3WHlea26JAKTqywaJFW9aDNwvXdTiuhQIlGhQo5UHt2jBokLk/bFiJfqVvY7MmZGVm6XYrtqJGE2jSC3BhwbfRv77EpjmfQbAAGnc3m4kWFBQfzxCBFspPS82XqbhAfWokJlkvX/amQCkvirq9hg0r0SLHi7oMwHUd/L4NLNxUBodedS7q9vo6+teW2BMMwqyPzP2i1sny5SZUqlSJyJThiesmAdAsqZv1smXfFCjlxeDBULUqrFkDUw++tUnzWvWoFDT/SD+f/0uka7e3TmeCNx42zoG1M6J/fYktK8ebI6IrVYdOZ5nnirq7OnY04ykWBYNB1uWaxZMDWx5vtWzZPwVKeZGYCOcUfusv4eB822o9AJixsQw2a0yqA10vMPenvhb960tsmfmhue061GwkChEdkB+zbDauNwM3GMfQLsdZL1/2TYFSnhR1e33xBWRkHPTtg1r1A2BTwfzor0cB6H0L4MCSn2DLkuhfX2JDxnpYWrgwt8cVxc9HMFCGL/oZgJqeTlSvVNl6+bJvCpTy5LjjzAKwHTvg3XcP+vazOvbBDSaAN5PRS/+MfP3+qU4baF+4cn7a69G/vsSG6W+CG4RmfaBuu+Lni9agRGCG17w00816VP0+1suW/VOglCceD9x9t7n/6quQl3fAtyclJFDTY3ZY/X5pGe2t1ec2czv3S9ixsWzqIGUnfR3MfN/c73tH8fO5ucUHx1luoaxO20yOx+wPdlGXgVbLlgNToJQ3F10EjRvDxo3w6acHfXv3er0AWJA2M9I127cmPaFpbzNd9Pe3y6YOUnZ+ew4CedCsL7Q+ofj5hQvNzK+aNaFhQ6uX/GTuOBzHxedvxOGNtOVKNClQypv4eLij8Jveiy/CQcZGzu1oNpfMdJazLXtnpGu3b31uNbezPoKc9LKpg0TflqUw53/m/oBHwXGKX/vxR3Pbp8+ez1swOWUyAG2r9bRarhycAqU8uvpqqFEDli6F7w+8X9bRTdrh+GvhOAH++9e46NTvn9qcBHU7mD2cJr1YNnWQ6Bv/pBk7aXcqNDlyz9eGF27Lc/bZVi+Z7/eTWjAHgMFtTjjwm8U6BUp5VLUq3Hijuf/cc+C6+32rx+OhbVWzav77FT9Eo3b7qgSc9JS5//s7mvF1KFj/Jyz6AXCg/0N7vrZsGcyfb7YTGjLE6mVH/D0dvNkQSOScThqQjzYFSnl1yy1QqRLMnAkTJx7wrdcdYdaDpLlzWLKljPbWajMA2p1idiEefe8BQ1AqgF+fMLddL4D6Hfd8rah1cvzxUKuW1ct+vfgnAOr6ulApLjJHCsv+KVDKq3r14MrCLSyee+6Abx3QuisJgRY4TpBXZ5RsC/yIGPgMeBPMWSmLR5VdPSSyFv0IKyeCJw6Ou3/v178t3N/trLOsXnZnXg5LMs3JpoNbnWq1bCkZBUp5dued5lCicePg558P+NbjG5t/YDM2jyZYgr3AIqJWC+h9s7k/9n4oyCmbekjkZKfBj7eb+71vgprN9nx97VrTqnac4sPjLHll2nDwZuP4a3DjUYOtli0lo0Apz1q2hJsLP6BvvdVstLcftx99Lm4wDr8vlW8XlsFWLEWOuQOqJUP6WpiqxY4Vzk93Q9ZmqNv+wK2Tvn2hQQOrl/5xlSn7iFonE+/zWS1bSkaBUt49+ijUrQuLFsGbb+73bY2q1aKRz8y0+Wj+V9Gq3d7ik+CkJ839Kf+GbSvKri5i16KRsOAbcDxwxlvgS9j7PRGa3fXL8rnkeJfhuh7u7n2p1bKl5BQo5V2NGvDMM+b+o4/C5s37fevQjmZzyTV5U9menRmFyu1HpzOhxbHgz4VvrgR/ftnVRezIToMfC9dH9bnVnHnyT6mpxTtlWx4/+b9ZZsPU2k43OtVvYrVsKTkFSkVwxRVwxBFmj68HHtjv2y7udjwef20cTx7/nvZNFCv4D44DZ7wNiTXN9va/Pl52dRE7DtbVBTBihJndd+SR0MTeh/727ExWZE8EYGiH862VK6WnQKkIvF54vXA84sMPYdasfb7N5/XSvbbZ22jcupHRqt2+VW8Mp79l7k//P1g6tmzrI6H7+/vCri7v/ru6IGKzu16c+iV4c/H4a3NVd+3dVZYUKBVFnz5me3vXNWtU9jOT67ZeF+K6DtmepUxdsyjKlfyH9qfAUdeZ+99dr80jy6PNi+C7wkW2++vqAti2DSYUblBqefxk3LrvADiq7qn4vF6rZUvpKFAqkuefh6QkmD59v2fPd2nQnBp0BuDpKW9Fs3b7duIT0KAzZG+Db6+GYBmc2yKhyU6Dzy+A/J1m88f9dXUBfPed2XeuSxdo3dpaFUYu+oM872pc18vdvS+2Vq6ERoFSkTRqBA8VbnNx551mEHQfbjriWgDWFkxiZsryaNVu33wJcM5HEJcEqyfDby+UbX2kZAIF8NWlsH011GgG530Cvv2sTA8G4eWXzf2hQ61W453ZZjC+vrcnberY3bVYSk+BUtHcfjt07gxbtsDll++z6+uCLv2oGuyE4wR55LdXo17FvdRpA6cWfuD89lzxDrUSm1zXDMKvngzxVWDol5BUe//vHz7cTGuvUQOuu85aNWZvWMWavEkAXH7YhdbKldApUCqahARzRHClSjB2rDmIax9u63ETAOsKJvP72mVRrOB+dLuw8Mhg4PubYGkZ7YwsBzfzffjzI8CBs9+Heh32/95gEJ4q3Bj0llugenVr1bh//Es4ngBJwXZc1PU4a+VK6BQoFVHHjvDKK+b+fffB7Nl7veW8zn2p5h6G4wR5dNKr0a3f/gx4HLpcAG4Avr4MUvY9W03K0IrxZnNPgAGPQbtBB37/jz/CvHlQpYrZzcGS8SvmkeI3557c1eN2PB59lMUC/V+oqK691uyVVFAAF14IWVl7veX27qaVkuKfwrQ1i6NcwX3weOD0/4PWA6AgG4adC1tjoPUkxooJ8PlQE/hdLig+OG1/XBeeLNwV4aabrO4s/PiUl3Ecl1ocwTmdtU19rFCgVFSOA++/b44LXrJkn98Oz+nch2puZxwnyOOTY2RfLW8cnPsxNDoCctLg07M0nTgWLP/FzOjy55gD0wa/dvCTFseONWuiKlcuPmXUgm/mTyWNv3Bdh0f73mmtXAmfAqUiq13bnDvvOPDBB/D113u95Y6eppWy3j81NlopAAlV4KKvoVYryFgLn54BGevLulaHrqXj4PMLzVY57U6B8z+DuEoH/p3dWyfXXWf2m7MgGAzy4qx/A5DsO4b+rbpYKVfsUKBUdMcfD/cXrg+4+mpYs2aPl8/u1JvqbhccJ8hjk18rgwruR1IduORbqNoQtiyGD07SSY9lYfFP8MVQCORDh8Gm9bi/lfC7mzABpk0zk0Tuustadd6bNYZsz1LcoJdn+9srV+xQoBwKHnsMjjoKMjLMtheZe24MWdRK2eCfyoSV88uggvtRszlcORZqt4YdKfDhQFj3R1nX6tDx9w/w1SUQLICOZ5j1Qvtba/JPRa2Tf/0LGtpZH+IPBPjPgv8DoF3lkzm8UQsr5Yo9CpRDQVwcfP451KkDf/0F5567x9kpZ3U6mpocjuO43Pfbo/gDMbRavWYzuHKc2dIjZzt8PASWjCnrWlVsrgu/vwtfX26ObO58Lpz9gRnfKokpU8yx1HFxcO+91qr14pRvyPeuww0m8MKA26yVK/YoUA4VLVqYKZyJiTBmjOnX3u1c93+f8AhuMI5szzIe+OWDMqzoPiTVhstGQusTzaDwF0Phr0/LulYVU0EOjLgORt9jZnN1uwjOfBe8pTiwqqh1cvnl1nYVzi7I48vl/wGge/UzaVXb7uFcYocC5VBy1FHw5Zdmeu6HH8ITT+x6qUdya/rVuQSA0evfY9HmlLKq5b7FJ8GFn0PXwmmrP9wEPz9itgARO7avMWNV874wOwcPfAZOfxM8pdhwcdo0cyS112vWQFly3cgXCfg2Q6AKL5x0k7VyxS4FyqFm8GB4q3BTyMceM8FS6NVBNxMfaAbeXG4Y83DZ1O9AvHFme/S+hWeWT33NjKukrSrbelUEKyfCf46D1HlQuTZc+h0cfePBpwbvLicHrrrK3L/8cnNEtQUjF83krx1mhuL5LW+mfhV7q+3FLgXKoejaa4sP4rrmGhg9GoB4n48n+zyO63rY6s7ilakjyrCS++E4ZoX2eZ9Apeqw/k94tx/ML8MDw8qzYNAE86dnmnU/DbvBNb9Bi36lL+vhh2HxYjMI/4KdTT6z8vJ4ZNpDOE6Q2vTggX4XWClXIsNx3d060uXQ4bpw2WVmnUpSEvz2G3Q3Z1mc//VD/J39PU6gOuPOHUmDqjXLuLL7kb4Whl8N62aYx90uhlNeMN1jcnBbl8HIW2FN4bG83S4ym3TGJZa+rKlT4ZhjzN+rkSPhtNOsVPGi4Y8xL3M4BKrw3ekjNHYS49RCOVQVraQ/4QSzLcupp5oV9cC7p92Px18H15vBtT8+VcYVPYAaTeHyUXDsveB4YM5n8O6x2gPsYPz55piAt3ubMImrDKf+24yXhBImWVmmi8t1za2lMPnu7xnM3WlayZe0vl1hUg6ohXKoy8iAfv3MBn516pgZYN2788Gssby68C5c1+HhI97k/C7HlHVND2z1FNNa2bnBPO46FE54BKrpjIw9rP0dRt5iFouCmTl36stmenaobrkF3ngDkpNhwQIrOwpn5GZz3LDT8ftSqefpxa+XvBd2mRJ5ChQxZ6cMGgR//glVq8IPP8BxxzHwsxvYEJiMz9+QSRd/T9WEEL69RlN2Gox9AOZ+bh7HJUG/O6HXjQffKqSiy82AX5+AmR8ALiTVhZOfg8POLt3A+z9NmAD9+5v7Y8fCSSdZqW5RtyuBqow84zua16pnpVyJLAWKGDt2mN2JJ0zYdabK6n69Gfzd6eDNpGlcf0YNjaGtWQ4k5U8Ycy+kzDSPazSDgU9D+9PC+/Asj3LSzSLFGW+aUAE4/GI48UmoHObuvzt3miN9V682kzvefTfc2gLw1fwpPPHnDTiOyxWtHueOvmdZKVciT4EixXJzzVb3331n1qq8/z5vdazLW4sewHFcTmlwC88PvLqsa1kywSAs+MasVdlZuFtxsz7Q5zazPX5FPz8jOw1mvA2/vwN5O8xzddrBqS+FNoNrX667zoRI8+amy7Rq1bCLTM/J4rj/DSHg20xDTx/GXfJO+PWUqFGgyJ78fvNt86OPzOOXXuLKVvnMzPgfbtDLy/XuY6CnOrRuDW3alG1dSyIvE6a+ClNfh0Ceea5OW+h1vTnTI75ymVbPuqxtpjXy+38gf6d5rl5H6Hc3dDy9dIsUD2TcOBg40NwfP95sQhqmYDDIKf+7mfWBSRCoxk9n/UCTGgc4WlhijgJF9ua6cM898NJLAATvuYfBTVZx/7uj6btgt40lBw40e4TVjNFpxbtLX2e+rf/1SfE39sSa0ONK6Hl1+R68DwZh7TSY+wUs+BYKCg9Tq98Zjr3HdPXZbJGlp0PnzpCSAjffDK/bOUvn2h9eYtr2j3Fdh1s6Pc81PQ9yGqTEHAWK7N/zz+/aPiPQsAFsSsUb3O11rxcGDDAzw8qL3B0wZ5jpDkov3MrfEwftT4EOQ8zhUZWqlW0dS2rrMhMi874y58YUadjNTKVuN8j+mFEgAOecY7pFW7eGOXPMOqYwvTJ1BB8sexTHcelf91peO0Xbq5RHChQ5sPfeMyvrD/TXZOnS8tH9tbtgABaPghlvwdrpxc9746HlceZbfbtToIqdg6GscF3YvhqW/Wz221r/Z/FrCdVMl1bXC8xYUSQmH7gu3HADvPMOxMebHYWPPjrsYn9a8if3TLsGx5NP8/gBfH/+yzojvpxSoMjBPfts8VYt+/LTT2bacXm1cS4sHAGLRsK25cXPOx5oejS0PsEcSdyom+kmixbXNa2QNVPNz+qpxetswGzg2HoAdD3fhF8oixJL46mnzPYqjgNffWVaKmFavCWF8364ENeXTpVgB369+FMqx5XgAC+JSQoUObilS6Fdu/2/Pnas6QopLwP1++O65lTIxSNNuGycu/d7arYwwdLocPNTr5OZfhtOi8B1IXub2e13+yrTCkmdD2umQdbmPd/riTNnw3Q6Aw47J3otqPffNyd+glnEeFP4XVLpOVkM+PxC8ryr8PrrMfLsrzQIX84pUKRkTj4ZfvnFBMeBlKeB+oPZvgaW/ATrfocNc8yH/b54fGahYFJdqFKv+H7l2oALAb859TCQb7bbD/rN/czNJjy2r4b8zH2X7U2A5J7QvI/pykruGf2ZaSNHmjVKwaBpqT79dNhFBoNBThx2LZuDMyCQyDsnfEyfZh3Cr6uUKQWKlMz27WaNytixB35feRyoL6nsNNNq2TgHNsw2P+lrD/prJeNA1Ybm2OOazaF2K9Pd1rh72a7ynz7d7PeWk2P26frwQyvjM5ePeJo/d3yB63q4s/NLXNH9xPDrKmVOgSKls2wZLF9uzqU/77z9v688DtSHwp8PWVtM11Rm0e1m81x2mhmH8cYV/sSb1ow33jyuXLs4QKo3ib3tYRYvhj59IC0NTjnFzOyKK+ExwAfwxITP+Hrt8wAManAzLwy8JuwyJTaU4lxPEUxItGmz6wyV/Vq+/NAIFF88VG9sfiqSDRtM92VaGhx5pBmEtxAmD//yX0ak/BvHgbaVTlWYVDAKFAlNq1YHfr116+jUQ+xLTzdjZmvXQtu2MGqUlbUmd415h7Gb3sRxINl3HJ+fHcNHI0hINNlbQtO2rfkG6917K4+UOolsqBfmxoNSNtLSzDHR8+dDgwZmLKxOnbCLvfHHVxm76U0AWiWczMgLXiXep++zFY0CRUL3+edmAP4fkrfm8Ncp3di0fWsZVEpCtnQp9OoFU6ZAtWqmW7NFi7CLveq755m07QMAOlQewrfnPY9vH19EpPzToLyEr2igvnVr5n4zjE4PPY4vCFM616fDr3OoXVcn7cW8CRPg7LPNbL5mzcxU4c6dwyoyGAxyyYgnzBG+wOFVz+e/ZzygVfAVmAJFrJvx9r/petvdJOYHmd+yBr6RE+nQsWtZV0v25/334frrzU7TvXqZ2Vz164dVZDAY5NyvH2Bp7igAete8jHeH3GWhshLL9FVBrOt1/R3M+PBtMip76bwynYQBR/PTL9+WdbXknwIBuPtuswLe74cLLjBb0YcZJv5AgNO/vHNXmPSve63C5BChFopEzPzxo6l7zpk02J5Has04Rrz0ONdfeX9ZV0vArCO66CJz3DPAY4/BI4+EvWhxU2YGF3x7G1vdWbiuw2mNbua5k8rJoWwSNgWKRNSWJX+T2b83LTZksKOyh7dvO587n/hUg7JlKSXFzOSaM8cc9/zRR2YXhDD9snwud/12BwHfZlzXy7lN7+TR/peEX18pNxQoEnH5mzezql932i1JAWD4Ce3o+/lE6muwPvpmzYIhQ2DjRqhXD77/3oybhOnx8Z/y9ZpXcDwFOIEaPNTzWc7r3NdChaU8UaBIdOTmMuOi0+n17TgAVjRKIvfD/9Fp4JAyrtghIj8fXnwRnnwS8vLgsMPgxx/NjK4w7MzL4aJvH2RV/s8AVAl24LMhr9Oqtr4sHIoUKBJVk9/5N+3vu4+6GQUUeB3m3XwN3V9+y+4RtbKn6dPNwPvChebx4MHw2WdmrUkYZm9YxdVjbiHPuxqArlXO4cPTH9SCxUOYAkWibuGCP9lyySkcN8ec9TGvUzOajhhLjTYHOHNFSi8jw2w3//bb5syVunXh1VfNeEmYg+9v//ETby14HLzZEEjk2o4Pc1OvwXbqLeWWAkXKxPasnXx8y2Cu+WwSlfNdMir7mPPoAxx7z+NlXbWKYcQIcwjWhsITHq+4wnR51Q7vAKvcgnyuGfkCf+34CsdxiQ805d2TXqNHsvZuEwWKlLHPvnqHbvfezWGrzQFTfxzRktbv/I9aPY8q45qVU+vXmyD57jvzuHVrePdd6N8/7KK/mj+F5/54mgKfmVzRLO4Ehp31DNUrRfnAL4lZChQpcxu2bWLSVadxzshZ+IIQcGDFqSfS9s33oWnTsq5e+RAMwjvvwH33wc6d4PPBvffCgw9CYnhnza9O28wNo59knX+ieSKQyNnNb+Kx/peGX2+pUBQoEjO+/OY9aj13Pyf+uQ2A/DgvBdddR9Kjj4fdVVNh+f3wzTfwwgswe7Z5rlcveO89M5MrnKIDAR769UNGpbxvxkqAxt5+vDnoYc3ikn1SoEhM2ZSZwQuvXMsFw0bSc4n5EMutUpn4+x/Ec9ttUFndK4BZ6f7BB/DKK7BmjXmualV47jm47rqwZ82NWjKLx6Y+Sa53JQA+fyPu6nE/F3U9LsyKS0WmQJGY9NGsccz87E5u/XY57dblApBVtzaJjz+J57LLDt1gSU2FN94wM7e2bzfP1a0LN99sNngM8+yS1J3buX7UsyzLHYPjuLjBBPrVuYiXTr6BynEJFv4AUpEpUCRmbcveyY0jn6P5r8O4acQGkrcWAJCfVJn4Sy6FK6+EHj3CngJbLixaBC+/DJ9+ahYpgjli+c474dJLwx4nSc/J4snfPmHchs/AuwOAus6RvH7SoxzWQONYUjIKFIl5i7ek8PCYf3PE6K+45OfNu4IFMGd2XHklXHyxlZMFY4rfD5MmmbUjI0cWP9+7t9klePDgfZ6YWRrr0rfx2MT3+CPte/CamXYef11u6HwX1x55Slhly6FHgSLlxrzU1Tz066vUmTuSsydvY8CsHST4C//6xsXB6afDVVfBiSeG/UFbZjIyzLG7I0eaExPT0szzjmP+fHffbQIlTAtS1/LYpLdZnPUzjifPXMJfk2MbnMNTJ/xLU4ElJAoUKXdmpSzn4d9eJSN9Aqf8vp2zJm2n45rc4jckJ8Npp0G/fnDssdCoUdlVtiRWrDABMnKkaZH4/cWv1awJ550Hd9wBbduGfakJK+fz/LR3SPFPxXECgBlwP73FRdxzzPkaJ5GwKFCk3Jq6ZhGPTXqNjYFptF+bzZmTt3PatB1Uz/bv+cbWrYvDpV8/aN68TOq7S3Y2/PVXcYgsWrTn6+3bm+6s004zrZEw98YKBoN8tWAKb81+n+3M3vV85WAbLulwBTcceaqO5RUrFChS7k1YOZ9Xfv+IVTlTiA9k0WdBJj2WZNFzcYB26zLwBP/xV7xpUxMsfftCq1amRZOcDFWq2KuU68LWrSYsFi/e87Zomm8Rr9fUZ/Bg89M6/G1MgsEgo5b8yf8WjOTvHVMI+rYUVsuhlnM4Nx1xtbaXF+sUKFJhZORm838zvmPU6pHsYCGO41IlO0C3pfkcv6wqx6/No86iZTh+/74LqF69OFx2/2nY0Ixh5Ofv+6egwNzm5ZnDq4qCo2j8Y1/q1jVjPYMHw8knQ40aYf/59xciAG7QR5P4Ptx79HUc1zK8BY8i+6NAkQppXupqXpvxBbO2jdvjg7VyVnUGpNTh5BSHI1LSSdq8xYRARob9SjiO6V5r3x46dDC3RfctzUjbM0QmE/Rt3fWaG/RRx9OVE5udxDU9TqNulfC2qxc5GAWKVGjBYJAv50/mkwXfsC5/+q4ZTUV8/gY0rdyF42t24dwarWmcmWMCZt06c5uSYhYTejwQH29+4uKK7+/+ExdnTkHs0MH8tGljfQFmdkEeE1bOZ+LqP1mwdQEbchcoRCRmKFDkkLEteyfD5oznt3XTWZU5l3zvehxnz7/+cYFkmlfuQt/koziqSQcOb9SyzGY+BYNBZqas4OeVM5m9aR7rshaT7azF8RTs8T6FiMQKBYocstamb+HL+b8xJWU6a7Pn4/dt3Os9ruvgCdQkyVOfWgkNaVwlmdY1mtGpfkt6NG5D/SrVQ76+PxBgTfoWlqdtZE16Kut3bCY1awtbc7aSmrOWHcEVuzZl3EOgElWcFrSo2p4eDbpxSbcBChGJCQoUkUJLtmzgm4UTmb7hD9bnLKLAs2Wv1sBeAlXwuVVx8OJxfHh23frwOF68jq/wJw4HyArsICeYjp8dBD07cZzgAYt3g14quU1onNiWbvW6cELLHvRu2h5feV24KRWaAkVkP4LBIEu2buCvjctZtGUVqzLWkJq9noyCVPLYvO/WQygCSfjc6lTyVKeKrxY1EmrRuGpj+jU9ghNbd6NqQnj7dIlEiwJFJEQbdqTx5/rlbMlKJy/oJ8+fT0GggLxAQeFtPgVBP/mBAvxBP/5ggNqJNWhctR5NqzegVa2GtKxdX6vTpcJQoIiIiBXab0FERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKxQoIiIiBUKFBERsUKBIiIiVihQRETECgWKiIhYoUARERErFCgiImKFAkVERKz4f3wEeM1yn6NrAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "ax = udp.plot(pop.champion_x, figsize=(5,5))\n", + "ax.view_init(90,0)\n", + "ax.axis('off');\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(-1.150591947301098,\n", + " 1.1494267321191918,\n", + " -1.3449433682343668,\n", + " 1.1413261458890887,\n", + " -0.049310762081661616,\n", + " 0.04878496706352893)" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcoAAAGFCAYAAAB9krNlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB0u0lEQVR4nO3ddXhUR9vH8e9ZiYfg7u4OxaFAkVIqQCnWB1qgAnV526fu9tQo9RYq0EKRUty1BYq7OwQIhIS4rJzz/jEJJCTZJLCShPtzXbkie3Z3Atn9nZkzc49mGIaBEEIIIbJl8nUDhBBCiIJMglIIIYRwQYJSCCGEcEGCUgghhHBBglIIIYRwQYJSCCGEcEGCUgghhHBBglIIIYRwQYJSCCGEcEGCUgghhHBBgjKfDMNA13Wk8p8QQtwcLL5uQGHidDpJSEjAYrFgsVgwm82YzWY0TfN104QQQniIJkXR80bXdRISEnA4HAQEBGQKx/TANJlMEppCCFHESI8yD3RdJzExEYfDgcViuRKIhmFgGAYOhwOn04mmaZl6mRKaQghR+ElQ5sIwDJKSkrDb7VgslkzhlzEM00PTbrfjcDgwmUwyNCuEEEWADL26kB6SqampmULSarXmGH7p/5wZ/1llaFYIIQovCcocGIZBcnIyKSkpV0IunZ+fX54fI/0DyNLLlNAUQoiCT4IyG4ZhkJKSQkpKCiaTKVNIgusepavHTP/QNE2GZoUQopCQa5TXMAyD1NTUHEPyeqX3INPPS3RdR9f1TNczZWhWCCEKHgnKa9hsNpKTk6/0+twtPQivnTWbHpgZZ9VKaAohhO9JUGaQMSTNZrPHny+7WbM2m02GZoUQogCRoExjt9tJTk7GMAyvhOS1ZGhWCCEKJglKwOFwkJSUhK7rPu/B5TQ0m17QQIZmhRDCu276Wa/p9VudTmeWggI5SQ8rb5JZs0II4Rs3dVA6nc5MpenyGjjXszzEXa4taJAxNGVoVggh3O+mHXrVdZ2kpKR8h6SvydCsEEJ4100ZlK7qtxYmMmtWCCE876Ybes2pfmt++OIaZV7J0KwQQrjXTdWjTK/fmpqaWmR7WjI0K4QQ7nXTBGV6/db0kCyoPUJ3cjU0KzuaCCFE3twUQZlbkfObwbUFDZxO55VepgzNCiFEzm6KoEwvcu6p+q2FiQzNCiFE/hT5oLTZbFdC0hel6QoyGZoVQojcFemg9HX91sJEhmaFECJ7RTYoC1L91sIkr0OzN/sQthDi5lEkgzK9NF1+6reKrGRoVgghimBQpoekruseC8mbrEYDkHVoNmMvU4ZmhRBFWZGqzKPrOomJiR4vTZfem7rZpfcy08nQrBCiKCoyPcqiUr+1MJGhWSHEzaBIBGV6SNpsNglJH8ltaDb9/0X+b4QQhU2hH3pNr9+akpLitdJ0MvSaN9cOzWa8lilDs0KIwqJQ9yhvxvqthcm1Q7O6rl/pZcrQrBCisCi0QSn1WwsXGZoVIhfx8XDoEFy8CImJkJCgPtvtULw40f4WtiRdYp8jhqMlLFzSUoi3x5LoiMdAx88UgL8pAH9zIIGWIKoWq0Tz8vXoWLURVYqX8vVvV6gV2qHXlJQUkpOTfVKaToZe3UOGZsVNKzYW1qyBVatg7144eBDOncvz3XUNTpf143CVAA5XCWBfjUC21A8m1S+H140zlFCtKk1Lt2Fg/R70qNVUXmP5UCiD0mazkZSUBOCTwCqQQZn+31gIe2QZN5s2DEOGZkXRdPw4TJ0KS5bA5s3gdGY5JLlUCc4V8+OyNYWkACdJASZ0k0ZokpNiiU6KJULJeCdhibYs903192Nf89psblmL9U0rcsz/MgnGOQxzTJZjNWdxKvs35646fRndqjeWgvZ+VsAUuqC02+0kJiZeqd/qizfR9KFCrzMMtOijaCf/Ros+hinqKNrlE6DbIfESmjMVAEfrhzDKNUav3gWKVfR+O29A+p+jrutXhmLTQ1OGZkWhk5QEs2fD5MmqB5lR3brQsyexTRrxfeIR5vptIbZY7JWbDd1CCa0RTUq1ok3FxvSo1ZyqxcuoGy9cgN271ceuXeqxz5y5+tiaBn37wlNPEXFLK9afOciak1vZHvkvscYBNJPjyqEmR2nal72D/+s4gpoly3nu36IQK1RB6XA4rlTd8WX9Vm8HpXZ2G6YDf2E6ugzT5RN5vp+hmdBr9kBvPhy9dm8wFa6zxuyGZv38/CQsRcGXkABffgkffwxRUepnmga9esHgwdCzJ9vMdt775wcOJ60EU4o6xhlIBb9W3FatO6Na9KZMSLG8PZ9hqMCcN099bNt29bZGjeCpp2D4cAgMJCY5kWl71rL42ApOpPx95bkN3UIla3ve7fYMrSvXdt+/RRFQaILS6XSSkJBQIOq3eiUonXZMhxZg3vI9pnNX/+gNsx9GlXbo5RpjlKyFUbIWWALBmYoWeQBT5CEMkwXTuW2Z7qdXaInj9k8wyjbybLs9IOOQbEBAgK+bI0TOkpNh4kT43//g0iX1sxo14MEHYeRIqFKFCwmxPLX4E/YkzEMzqeFXs6MsPSoN4qUuIygVFHrj7Th6VLVj8mQV2gBly8L778OoUZB2fTIqKZ6P/p7OirNzsJlVj9TQzdQP7stnvZ6TSUBpCkVQZixyXhB2AvFoUBo6pn2zsax9Hy0uXP3I7Ide/070erejV+8G/iF5eigt6him3b9h3vELWmo8hsmCs/0TODv/H2iF60K+ruuYTCb8/f193RQhsrdiBTzyCBw7pr6vXRteew2GDgWLBYfTyZurp/DX6R/AHAdAkF6X/zQYxcNtbvfMdcKYGBWWX3wBp06pn7Vtq3q7bdpcOUzXdWbt28BnW78gwXRA/dAZSOcyQ/m0z3gCrH7ub1shUuCD0lv1W/PDU0GpnduOZfnLV3qCRlBpnC0fwNlyJASXvf4Hjo/Asvy/mA8tBMDZaBCOO74AU+FZHZQ+3O7nd3O/YEUBFBUFzz4Lv/yivq9UCd59Vw11pr1PbAk/yvhlz5FsViFqcpRmVP0neLL9Xd6ZfWqzqR7mG2+oHqamwejR8OGHULLklcN0XeebzQv5cd+XOCxqFq6/swZf9PiIDtXqe76dBVSBDkpd1wtkaTq3B6UtEcuqNzDvUC80wxqEs8NTONs8DNZAtz2Nac8fWBY+hWY4cdbpg+OeH8FcOIInfchdglIUKBs2wL33qqUdmgbjx6uQLHb12uIH6/7gt6MfgzkFQ/enQ6n7+KT3Y4T6u++1nWfnz8P//Z+afQtQrZqabNSqVabDbA4Hr636iYXh36e12487Kj3Cez1H35TLSgpsUKbXb01NTS1QIQnuDUotYjeWuY9gij4KgLPJfTi6vgyh5d3y+NcyHVmG5a8xaI4UHK3H4rztXY88j7s5nU6sVitWq9XXTRFCTZ756it4+mlwOKB+fTXE2b79lUNikhMZMedlTtlXAhDgrMm3vT+lVaVavmr1Vf/8o65VHjsGfn5qKHbMmCzLy3acO8EjS58jyXQYgJK05I8Bn1M+tIQPGu07BTIofVG/NT80TbvxN2zDwLz5G8xr3kXT7Rgh5bH3/xKjehe3tNEV05GlWGfdD4D9nkno9ft7/DlvlK7rWK1W3yzLESIjux3Gjr061Dp4MEyaBCFX5w4cuXSeIfPGYDOfxjA0mgTfzaS7XybIWoCuscfEqAlG8+ap70eNUuEfFJTpMJvDwWOLPmND9G9omhOrozJT7vieRuWqeL3JvlLggrKghyS4ISgdqVgWP4N570wAnHVvx9H3Uwgqmcsd3ce8+i0s/36J4R+Kbcw6KFbJa899PZxOJ35+fhKUwrdSUlQwzp8PZrOa3frUU5l6YlvDjzJ66UPolkhwBvN4k9d5qE1fz7bLkQqW6whhXVe/w0svqa/btVMFEcLCshz6576NvL7pWTDHozlK8Hm3L+leq6kbGl/wFaigLCz1W28oKJOisP75AKYz/2JoZhy3vYve8gHvV9Rx2rFOvRPTuW04mw7F0W+Cd58/H9KXh/j5+RW8ikji5pGQAHffDStXQkCAurZ3++2ZDll+ZCfP/D0ezHFojpJ81f0bOtdo6J7nNwy4fBJOb4RzOyDmDIRvgaS0ZSjWIEADswWKVYbiVaFENSjbEKp3gpI1c36fWbMGBg6E6Gi45RZYujTbsNwSfpSxSx/GabkIzkBeaPkhI5rf6p7frwArUEHpy/qt+XHdQXn5BH7T70OLOYnhH4r9nkkYNbq5u3l5pp3dit+vt2NoJuxj1mKUrueztriSvvOIv79/gf67EEVYair07g1r16oh1vnzoVu3TIcsPLSVF9ePA3MyFkcFpvT7kcblq97Y8+pOOPk37PoDjq+G+PPX/1ihFaBGF2gyGGrdmrUAyc6d0KNHrmF5Mvoi9/71ECnmYxi6H2+2ncjARh2uv12FQIEJSpvNRmJiYoEPSbjOoIw+jt/v96DFn8cIq4p98G8FIpgss0dhPrxIDf8O/NnXzclWeo/S39+/wI4yiCJM12HECJg2DUJDYdkyNUSZwbazxxi19D9gjiPAWZPZA368Wm7uesRfgE3fwq7pEJ+hWLrJChVbQJW2UKqWCj/NrEpV+oeCoath2NhwiDmleqBnt6mepzNDfdhilaDZUGgzOnOZy4xh2bat+l2zCcuY5ERun/4g8ab94Aziq1sn06VG4StmklcFIigz1m8tDNeg8huUWvRxrL/djZYQgV66HvahsyCkYNRU1CIP4vdjFwzNhO2xXQWmXRlJUAqfeuklVdHGYoFFi+C22zLdfCYmiv6zh+C0RGBxVGThvdOoWOw65xskXIT1E2DLJHAkq58FhEHjgdDoHqjc5vqWjNmSVFgeXAC7Z0BKjPq52R/ajoVOT0NwafWzjGHZvbsKy2w6LxcSYrl9xnBs5lNojhL83m/qjfegCyifB2VBqd+aH/kJyiwhOWx2luIBNoeDs3HRnI2PItlhz3SbWdOoGFqSKmGlCfbzzIw566+3Yzq7FUf3N3DeMs4jz3Ej0oMyICCgUPx9iCJk2jQYNkx9/dNPamZoBrEpSdz2+wiSzUfQnGFMvX0qTctXz//z6E7Y+BWseR/samckKrWGjk9A3T7XN1EnJ45UOLQINn2nrncC+IVA52ehwxPqGufOndCpk9oP86WX1NrQbByLimDgX8NxWi5icVRgyX0zKReStQda2Pk0KAtS/db8ytPC98SL+P1yO1rsac6VqsPCZs+wJ+4Sp+JPcTElnCT9ErqWgGFKRtNy/28wdH9MehABWilK+1WmSkg16peuSduK9WlZocZ197ZMO37BuuR59DINsI9eU+C26tJ1HUCCUnjXiRPQvDnExeUYFnf8/pRaJ6kH8HGn7+ldp0X+n+fiAZg7Xg2RAlRqBd1egto9PPtaNAw4uhJWvQXnd6U9d2u451soXQemT1fl90AtIemf/TKyreFHeWDZSDDHUcHUkWX3f+u5NvuIz4KyoNVvzS9XQelwOll9ZDMR/zzHEaLZGhDERUvuv5+hB6AZ1/ZUnXkLUmcxSpvr0bhUc3rVuIWeNfOxMWtKLH5fNEJz2rA9tBGjVAFYEJ1B+pZbUhBdeI3dDl26wL//QseOalboNZeFPvlnNj8fewPD0Hi84Qc83Pb27B/LlW0/w8Ln1FZ5/mHQ+11oMcK7J6uGATt/hyX/hdRYsARAnw+g9QPw5JOqTmxYGGzfDjVrZvsQv+1aw/s7nkDTDO6u9Cxv9xzlvfZ7gU+CsiDWb82va4MyITWF3/asYenJVZxJ3YqRVvQ4nWFomJ1lCDVXoFxAFWqGVaNuyWpULlaGKmGlqVysFIE5FB52OJ1EJMQQHhdFeNwlDkef5ujlE5xLOk2M4xx2cwSalnkTWM1RglrB7RhUty93178l14LL1ql3YTqzEXvfT9Gbj8j/P4gHSUF04XVvvQWvv64CYtcuVeotg30XzjBk4SAwJ9Eg8E5mDM5nhSunA5a9rCbsgBpeveMz3+4fGxsOcx9Ts2sB2j8GXV9V1yk3blS963//hRxehw/MeY+tcdMwdD++6TbFfctiCgCvB2VBrd+aX35+fui6zsz9G5iyfybnHVvRTKlXbg/SdZqk2vErdhtNa/SlX922lPXQ2H1sSiLzD29h3ektHI7bQzxH0ExXr3VqjuI0LHYrL7QbScOylbN9DPPa97Fs+Axn48E4+n/pkXZeLymILrzq+HFo2FAtCfn996vDj2kcTiedfx1KgukAVmcV/h4+h+D8nMTZkuCPEXBMlbaj+yvQ+bmCccnDMGDdx7D6HfV9g/7Q5nW4pYPaNuz99+HFF7O9q83hoPOU+0gyHcbqrMI/988pWJWIboBXg7Ig12/Nj/C4aL7cNocNkYtwWDKsa3IWo4VWgYcu/k3b5BQ0H5WHi01J5Oedy1l6agUXnDuuBLhhmChnasPYpvdzT4PM09u146vw+2MIRvGq2B7d6vU2uyIF0YVX3X03zJ2rZn4uX54lwF5aNon55z/H0K182fVXutVsnPfHtqfAtPvg+BpVIGDA9yqMCpo9s+CvR9WSkto9wdkfRj0AwcFw+DBUzL7nm7Gn3bXUGL6840kvN9wzvBaUhaE0XW5OxkTy9oZJ7E9aejV8dCuVre0Z0WAgA8qWJnBKHzRbIo72T+Ls9rKPWwxxKUl8t30R807MItl89MrPQ/SGvNT2GXrVbq5+kBqP/6fq2mTq00fUlPQCQgqiC69ZsgT69lXXI3ftUj3LDC4kxNJzRl8wx9O55IN83f/pvD+2wwZ/DIcjy8AaDPfPgaq3uPkXcKOT6+G3QWoWbsMB8NkB+HcT/Oc/V+vcZuP5pd+xJOJLcAbw113zqVXKMxs8eJNXgrKwh+TJmEje2vAjB5KWXQlIi6MCXcr15+l2g6kQWhzsSVh/6Ysp8gB61Q5qrWQB2+9x0eFtfLtzCmedG9E0J4ahUcnckXe6PE3T8tXwm9AQLekStgdWYpRv4uvmXiEF0YVXGIbazHjbNrUryKefZjlk6MzX2Js0B5OjDBtGLM77kKthqB7armlgCYQRs1RZuYLuyAqYNkRNNireH57+Tf38339V9Z5s2BwO2v3SH7slnGrWHiwY9rn32ushHg/KwlK/NTs2p4O3/pnKqotTwKwW/1odlRlc4z880e6uTBWELIufxbxzCkZwGWwPrsq6cN9uh3Pn0CIj0SIj4eJFtNhYSE6GpCT1WdfVwl6zGUwmCAnBKF5cTSgoUQKjQgWMqlWhePEbup6xJfwor67/lEvGZkD1im8rO5qPw2dgPr8T+4DJ6PXuuO7HdzcJSuEVy5apMnWBgXDqFJTJXFlna/hRRi2/F83k4IFab/JMpwF5f+wtP8LCZ1UVneEz1HBmYbF3NswaDRiwoyXMW6Oq9mzcqN6nsjF152o+3PUEhqHxYfvJ9KvX2qtNdjePB2Vhqd96rUVHtvPRjv9hs5wCwOKoyJCaD/JIyz5X9qNMD/0r1/fQsPeZBDHBmPbuRdu/H+34cbQTJyA8HC1tPeCNMkJDMapVw2jQAKNpU4wmTdCbNYMKFfL1OHMO/Msn2z8jJW3X9capgUy4cJSSXV8tMIUHpM6r8Jpu3VQt1yefhM8/z3LzbVMeIUJfT7Benw0j/8j7Sf+ZLfBTX9Uru+1tVUSgsFn9Pqz9AFIDYWKCKkTgYm0lQPdfxxBpbKI4zfl75BQvNtb9PBqUhal+a7oUh53Hln3KvuR5aJqBoQfQueRQ3uoyigCLFV3XcTgcWCwWLIaBZdu/+H0xEu1oNMaFALSohBwf2/DzgzJlMMqWVZ/DwiA4GCMoSJ3FapoaonE6weFAS0iA2Fi0mBi4fBnt7FnVG82BXqcORpcu6F27ovfoAaVK5fr7OpxOnlvxLX9HTUEzOQhzOnne2oI+gydfzz+f20n5OuEVW7aoXpLVqma9Vs48O3x3xEmGLbkLTdN5u80P3N2wXQ4PdA1bEnzdTtVdbXAnDP61YMxuzS/dCb/cCaf+gfWhsOKsupa7aFGOd1l3Yh/j1g5F0wy+6jK9UNeC9VhQFrb6rQC7L5ziqbWvkGJRPawy3MJHXZ6nXmk1w0vXdfSTJ/FbtQq/VauwrFuHFheX5XGMmjXRmzTBaNQIo3ZtjBo1MGrUgHLlbvxFkpyMduaM6qnu3Yu2ezfanj1oBw+iZfivNMxmjO7dcQ4ciH7nnVDSde3JNSf28r9/xnPBmohmQO9y43mn+wM31lY3kKAUXjFuHHzzDQwfDlOnZrn5vpmvsD9pLkF6XTY9MDvvj7viDfjnM7Xt1biNEFDMfW32trhz8E1HOHtJ9So1DY4ezbEIAUDnX+4nhp1UtXZn4bCCu5VfbjwSlIWxfuv325bw8/H/gTkZQw9gUOXHebbdPQBoJ09imTsXy59/YtmxI/MdA8CoasF++39w9r4XvXlzzMWLq/t58/eOicG0fj3a2rWYVq/GtGfPlZsMPz/0IUNwPv44RpOcJ+mkrv+c9w98x4KQYAAaBNzFT/1fyrVYgSdJnVfhcSkp6rJFTIxaDtIz8/XDqKR4uk3vAebk/F2bvLAfvusMugOG/A71+7m/7d627ReY/wT8ngpHUuH//g8+/DDHw6fsWMVHu5/E0C3MvXNxoZ0B6/agLIz1W19Z+zMrL/2Aphn4Oarxcad3aF28Apa5c7H++iuW9euvHGuYTOi33IKzR3csSTMxBZ/H3no0th7vkPGfUtM0TCYTmqb55N9AO3oU06xZ6mPv3is/17t1w/HGGxjtsg4dmbZOwrL8v7xdoRkzAy4DUNuvD7/f/ZbPenMSlMLjZs6EwYPVcOvJk1l2ynhm8Vcsv/gtJkdptoxcjl9eR8h+vVtVual/Bwz5ze3N9gndCT90h9Vb4I9kdXknPFxtZJ3d4bpOm5/vxGY+RcvQIfwywPdL5q6HW9/90uu36rpeKEJS13UeWvw/VkV9j6YZVDR1YUnrt+k48UdC6tUj8OGHsaxfj2EyYe/UiaRPPiHpyBFSV6yAniUxhUZghJTB3vmFLMFoGAZOp/PKR/obvrcYtWvjfPFF7Fu3Ylu7FuegQRhmM6Y1a/Dr1g3LyJHqDzwjawAa8HJQWfqUfQyAo7YljFv8mdfafS3DMAr835Eo5GbMUJ9HjMh2O6m15xcD0KX8PXkPyTNbVEiaLND7PXe11PdMZrj9f1DXAsU0iIqCWbNyPtxk4s7q9wGw8/LKKxscFDZuC8r00nQOh6NQDLfqus6wBa+zN3kOAH0udWL+/EhKt2qD31dfocXGolevTsorrxCzaxcJ8+ahP/QQWvnykBiJdb1aY2Xr8l/wv3rdIT0o00MzPTAdDgdOpxNd170amADGLbfgmDoV24EDOEeNwtA0zH/8gV/Tpph++unqgekbu5qtvNN9FJ1KqGuUW+On8fLKn7J5ZO8p6H9PopByOGDFCvX1nXdmuXlr+FFs5lMYhsYTt9yb98dd95H63GwIlKjm+tjCpkpbqNMTWqZVyko/0cjBuLZ3YegWdEski49s90ID3c9tQTlixAgWLVpUKHqSAA8v+ZjTjpVUiLQz6dcgPnr+O/xmzkRzOnF07UrSrFnEb99O4lNPYVSsiNVqvTL86LfuAzRbPM5yTXE2GZLjc6QH5rWh6ateJlWr4vj2W+wbNqB36ICWlIT10UexjB0LSUlotrQZu34hAHzeezyNgtR12iUXv2XhId+UtpNJPMJjtmxR1yaLF1fFBq4xacc8AEKMutQpncflVxF7VPUdzQSdnnFfWwuSDo9DnbTe9bq1aqZ+DsqEFKOUqSkAv+2Z543WuZ3b3oGio6MJDw8vFCH5f6t+4Gjcnzw++wKL/nuUtqs2oxkG9jvvJHHdOpLnz8fesyeOtO2d/Pz8rq6ZvHQI855pANh7vqNeDLlw1cvUdd3rvUyjRQvsK1bgePNNDJMJ85QpWHv2hOhL6gC/4CvH/nTHfylJSzTNyVtbXuNyUqLX2imExy1bpj737JllGy2ALZFrAOhQ/ta8P+aOtFmzDe6EArZlndvU7AbNmoI/EBsH105yvEaPqrcBsC/2n0I5/Oq2oCxTpgyRLtb4FRRfbVmA/vcX/PXyER6aH4nF4cTRrRuJq1aRMnUqevPmV9ZKapqWqScJYN3wGRoGjjp90StlPQPNzbW9TF3XfdPLNJlwvvAC9kWLMEqXxrR9O+aXp4DNwEjrUarDTPzQ511whuG0XOTBRW96vm0ZyDVK4VEbNqjP3btnuWl3xElSzScwDI2xre7K2+M57bBnpvq6RcHars6tNA3aPgjV004uVq1yefjDbfpfGX5dfmyXFxroXm4NykuXLrnr4Txiz+nDVH7vKb775BSVI+3olSqR/NtvJM+di976aoklp9OJpmlYLJZMhRK0yIOYD6qhA3vH526oLQWll2l066bCsnhxtAPnYEYSWEMyHVOteBkea/gSAGccq5h3cIvH25XOV7OGxU1iV9qbdsuWWW7668A/AAToVWmQw/Z0WRxdAUlREFwWauajF1oYNbgTaqQF5bKcCw8AlAsJoxj1AFh6dIOnW+Z2N01Q2o4dpXTf7oxYeRGApLFjSNyyBUf//pmKADgcDgzDwGw2ZymUYN3wiepN1u2HUdZ9m5L6updpNG2Kfe5cDD8zHHOirTyc5ZhRLXpQTmsPwKfbv/BYW7K0zfv7ioubRUQEXLigXv+Ns26VtS1CDSdWCcrHa/3AAvW58UAwF45CK9ctpCx0SBtVW78RbDaXh9cqpirz7Lm0x+VxBZHbgrJs2bJERUW56+HcyrRpE9aunal7Jo6oUAu7vvsC5yefQkjmnlP6rFSz2Zyl5J4WeQDLIfUisHd81iPt9GUv07jlFrhdFYE2fTkTzp/PcswbHZ/EMMwkmA7w607XQy3uIj1K4THpvck6ddQ+i9cITz4EQKvyzfP2eIYBJ9amPWYhKnp+I3oOhSANUmywe7fLQztUaQXARdshb7TMrdx+jbKg9QDMS5cS2K8fxeIS2V8tgB+/eZ+aQ0dlOS69F3dtwfN01q0/AKjeZJkGHm/3tb3MjD1Mj/QydQc0TYGKJrT4BMwff5zlkDaVa1PDT13Lmbzf88tFvD4rWNxcTpxQn+vXz3JTTHIiqaYzAPStk8c9Iy+fhNgzau1k1fZuamQBV609lEx7r0z/98xB/3qqyIluucSqY65DtaBxW1CWK1euwPUozStXEjhiBCabjdXNQ3nkhdt45I6xWY5zNXkHgKRLmPf/CYCj9UPeaPoV6T2q9LWpHutlxp1F03SMW1Uv2/znn9lO+X6x3RgA4rWD7Lt45safNxfSoxQeExGhPmez687iw1vRNB2cxWhRoUbeHu+kuqZJ5TaZZo4XaeUaQVjaEPPhnS4PrRx2td70hxu/9WCj3M+tPcqoqCicLtbTeJNp1y4Chw1DS01lZctQnh5fjcdueTZLCGYMyYzLQDKy7JqK5kzFWa7pdc10dRdP9jJNl44AYLSqjREWhnb+PNqGrBfdW1eqRaCzNppm8M22P2/o9xHCp9IvL5TPWn90X+RxAEK0Snlfx3txv/pcMevEoCLLbIVKaScah3bm+W7RtrOeaY+HuPUapWEYREdHu+shr190NIH334+WnMzmRhV5blwVyvh1pG+dzH/AufYkAZx2LDt+AcDRemyB2CLHE71MLUJdrzEqN0Pv3RsAU4Yatxl1KqfWRG2L8t51SiHczkWP8kzcOQCK+5XLcluO0k42KV3nRltWuFRJqzx05nSuh9YNUIXhqwRlHe4uyNwWlEFBQYSEhBSItZQBTz+N6eRJUqpU4alHi+OwmHi8+cgsx+W0DCQj89GlmBIiMILL4qyftcSVr7mrl3klKMs3x6iothUjmy3EAMa0UH/sdstZzsfHuOX3yI5cnxQelZhWPCM0NMtNF5NVb7NcUD52u4i6SYOyWlX1OSL3VQ+VQioBEJ16wZMtcju3BaWmaZQqVcrnS0TMK1dinTMHw2zmrXF9iA8xEeisw601Mk//drUMJNPj7VWLhx2NB4PZz6NtvxGuepl5qTFrilAX1/UKzaCYql2rxcZme2ytUuXRHCUAWH7Ms7Ub5Rql8Jj0CjHZnCTH2NQysiqhlfL+WDFpPaoSebymWVRUr6s+R2b/fpHp0DD17xnv9H2HKj/cGpQ+r85jGPi/rLZxsY0dy4IyBwHoUzlzT9DVMpBMki5hPqGGFx2N81EQ2cdcrcvMNjATI9Hiz2GgYZRrrPbnAzXdPQclLao016bznp29JiEpPCY9KLO55JKsq4mJtUrksdCAPRGMtMcLLOGO1uWJw5HAyZNvs2FDFdasMbNhQxVOnnwbhyPBa22gUtrGzYmpuR5ar4zqfdooAJfo8sGt1aZLly7t0x6l+e+/Me/fjxEczPqR92FYLmPoZh5s1vvKMbktA8nIcmAumu7AWb4ZRqm63vgV3Cq7XqbD4bjykR6apjP/AmCUrgd+IZi2quLneosWOT529ZDaAJyOP+nx30EIj0ifeJjNe4CuqcXzJQKzDstmKzU+7bEsYA10R+ty5XAksHNnV06efAObLRzQsdnCOXnyDXbu7Oq9sPQLSmuQAUeOuDy0XHBxAAwt91AtSNwalGXLlvVpj9KatmWUfcgQ/ry4E4BgoxalgtSShzxN3snAvE/ts+ZsNMhzjQaczgTOnv2AHTvqsHlzKDt21OHs2Q9wOt33h+6ql6kdXw2AXqMr6DpaWlAa2eymkK5EQHEAUpyeK5Iu1yiFRwWlvcEnZv0bNnAAEGDN4+WW9KD0C/HahL/w8M9ISNgJXFtkXCchYSfh4V7YRzY6Gp5/S33tMKBuXejTBy5fzvbwK/t5agVjdUReFZ0epWFgXquqYjgGD2Z3tLp21iBMzXTNyzKQjLTYM5gjdmJoZhz17/ZYs53OBA4c6M3Zs+9it58DdOz2c5w9+y4HDvR2a1hCNr1MXcd0Yg0AjqqdYNEitLg4jNBQjGzKeqULC1Bn2jY9ya3ty669QnhE6dLqc7bvWSoogyz++XtML/69njv3PVlDMp2edruHDRsGm68pcr5iBQwdmu3h/mnzPDRNL1S7iLi9R+mroNSOHcN06RKGvz/Oli2Jdap1OrdUaJrvniSA+dhyALVuMri0x9odEfElSUm7ye6sMClpNxERX3rsuTVNwxx7ClNcOIbJir1iWyzvvw+AY+xYDBfXb8P8VS/dbiR7tH0SlMJjXASlkdbjCcxrjzJ9uNXm2RPHjGy2czd0+w07fBiWLr16rTed06l+ns0wbKDVerV9Todn2+dGbg3K9KIDvmBKWxOlV62Kw2rBaVLtaFS6+pUiCK6WgVzLfFTtU+es3csDrb3q4sVJuDorVLd7jvmk6oXrldpg/ftfzNu2YQQFYXvssSzXMjNKsqsJPybNmuUx3UWGXoVHpQflxYtZb9PUa9KS18Lm1rRKPM5U0L0zrOjnV/GGbr9hx465vv3o0Sw/8svw75nisLu7RR7jkWuUPnmDS0gbogwJ4WhUBJrJiWGYqFuiHIZhYLFYXC4DySQ1AdOZjQA4a3m2uLHdHnFDt98o82G1PY6z9C34Pf00AI4xY6BsWeDqkPW1gRmRqM7Cg8xhHmlX+j6U0qMUHlMtbaF8Nm/4mq6GXC8l5r7kAchcsi4lj/e5QRUrPkTOb+GmtNs9qFYum1LXrp3lRzEp6nqwYWgEWfM5rO1Dbg9Kn9V7Tbswr8XGcjJWTSjS9BCsJnPuy0CuYT61Ds1pQy9eA6Nk1v9sd7JaXS9ozu32G5IYienMBtANzBNWYjp5Er16dezPP39l8o/ZbL4yVJ1xxuylJBWUoRbPTYWXHqXwqIZp22ft25flJrOhLi2cj8/j+5k1AELSqvjEnHJH63JVufLThIQ0J+vbuImQkOZUrvy0ZxtQty707g2ma05mzWb18zpZCy+cikl/bw66OrGnEHB7UEZHR+NweH/sWU87u9FOncKSPvZtkKdlINcynVal25w1b/X4xfmyZUfj6qxQ3e5+2pEjWH/5EO2SHX1rCcyr/8EIDCR12jQoWTLzsRlmzILqZV5MUZVLSgaU8ligSY9SeFSDtF2AIiPVRwYWTfUQIxLysd6vZNp6wujj7mhdriyWEJo3X0v16m/g51cZMOHnV5nq1d+gefO1WCwhuT7GDZs2DepfU+avZ0/182ycTevEmI08LrspINwa6WXThuuioqIon02hYU8yKlZEL1MGU2Qk5TZtp1pCKlUiovBrdQqtbv7WQJrDNwOgV87j9jo3oHz5x7h8eX42E3pMBAU1pXz5x9z7hNHR+D/wAOYVKzI8kxoOsU2ciNG0aY53zRhckfZDYIE25ZpemSiVvvREwk0UCsHBUKOG2h5q71649dYrNwWYQkkFopKzX+aQrZI14fRGiPJOUIIKy+rVX6V69Ve99pyZlCgB95aDNyMgwB9278m2J5nufIIaifI3Fa6gdGuPMiAggNDQUN+spdQ0HAMGANDy9fdY8OIRvvn8EMEtWuB/1105ruvJIjUBLVLtAuCNnULM5hAaNFhKpUovY7VWBExYrRWpVOllGjRYitns3rNC/wcewLR6dZaf69Wq4cxhSve1dl84rYo5GCburNdW3T/tWmZeyuXlh4Su8KjWrdXnf/7J9ONgi7r2fjEpH7P4y6QV+j7n2bKOBYo9GSLT9qEMCnIZkgAXk1QPPdDkmbkNnuLWoNQ0zadrKVPT3uj9ozIPl5hWrcJ/1Kg8PYbp/DY0Q0cPq4IRmnVXAU8wm0OoVOlFWrQ4Qtu28bRocYRKlV50e0hqR45gXrECLZut0EynTqFlM0stO/MPq+23/J2VKRMclulaZnr1n4y7mVxPaMr1SeEV3dVG5KxcmenHVUNVvdbT8a43I86kWgf1+dSGrEsmiqqz269WOPIPyPXwU3Hq+m2pgHzsylIAuD0ofVnv1Rmo1jJd2wfRdF0FxPbcz/RMZ7cB3ulNept23PWQkJbbdO80q8+ppTN1QjNvW5bxWqZhGDfcy5RhXOFxPXqozxs3QtLVNZBNyqjLNVG23LeOuqJCM7AGQUoMRB5wYyMLsCNLIf282y/3Nafnk9V7UINS9TzYKPdza1CC76rzOBwOTCdcn/35Pf54ro9juqT+wPVyOV+rK6yMSq53QjBym+4N7Io4RZymhqYfaZ59ofhrZ8xeby9TQlJ4XO3aUKUK2Gywbt2VH3eo2ggAuymCJHse65KarVAlbV7DkeXubmnBdHgZxKf1nq+ZBHgtXddJMtSJR4cqhev91e1B6Yt6r+k9ltzW9Zh37sx1eNF06TAAeunCVwTdpchI/J56CoBrI8owm3H27ImRzbqnaz28+iE0zcDiqEjbyrkfn10v0+VOJtfcVwiP0jTop/ZX5Y8/rvy4eYUaGLo/msnJpjOH8/54Dfqrz3tnu7GRBVT0cdVzjkx7DbsoeQmw4/wJMKdgGGa6Vnd9bEHjkR6lN9dSZtwNxFS/Ps7mzV0e73J40WlHu6yGBgrjbiHZMgzMs2YR2Lo15o0bMUKD0Spm/m/Xb72V1J9/zvWhzsfHYFhiALBo+dshIWMvM7udTK7dYFquUQqvGTFCfZ49+8rwq8VsJtBQW2wtP7Y574/V8G61g0jEbrjkeieNQm/7r+pzUin1uVEjl4evPaFqwlqd5Qn2LzzFBqCQ9yizq+Fq++ILl/dxNbyoxZxA0x0Y1mCMUA+Xf/IC7fx5/IYMwX/kSLRLl9AbNMD5Zh8YG4Ljf31I+fNPknftInXuXDXNOxcvrf3mytez+n/j4shc2pXNuszsrmVKj1J4RYcOUL06xMfD/PlXflw3rDkAmyL+zftjBZeCmmnLTHb+7r42FjQOG2yfor6OTltlmEuPcl24qnZWMbDwdULcHpTeqvea024gRqtWahjxmko8BuBs187l8KIpWl3j1EvW8uouAG6XmIjls88IaNUKy4IFGFYr9pdeInXRNMxxaRtR3/MKeu/eeRpuBTgQeZZ9iQsBGFrlJcoE3/j07px6mU6nk7lz53L27Nkbfg4hcqVpMHy4+jptqz6AvrW6AnDRvgdHNjPFc9RqlPq8ddLV7beKmgPzIOkSBJWHE2nF110Epa7rnEhU2/d1r9rVGy10K48Epafrvea2G0jqzz+jZ1g8DGomrOngQUxLluT4uFriBQCvLQtxu4QELJ9+SmDDhvi98gpabCzOVq1I+ecf7C+/jGXH92iGE2f1rugVct6U+Vq6rvP8mg/QTHYCnLV4om1/tzc9Yy9z586djBkzhhO5TM4Swm0efFBt4Lx0KezeDcCARh0wdH8wJ7Lw0Na8P1a926FUbVXzNX14sijRnbDuf+rrMn3BboeQEKhaNce7rD91EN1yCcMwM6K5Z+tne4JHhl493aPMdTeQEiVInTuX5F27eOqpDtz3Wk0O1yyPFhNDwMCBWF95JdNU8CsS1ZCxEVzGk813O+38eSwffURgo0b4vfqqGmatVYvU774jdfVqjMaN0eLPYdmjJivYO+SvBuR766cTyWYMw8SzzZ/OVznA/EpMTOTBBx/k6aefplu3bh57HiEyqVkT7k2bxf3RRwAEWf0paVL1YOceXpX3xzKZoMMT6usNE8Hmuc3NfWL3HxB5EAKKQ0zaTPoWLVyOwk3buxSAUKMu5UIKV7EB8FBQxsTEYLd7ZguV9IkfedkNxKhdmzqDXmB/zSDue6Ec5+5XBQmsn31GQPPmmKdNy7QwWEsLSoLLeqTtbmW3Y54/H/9BgwioVw+/N99UAVm7Nqnff0/K9u04R4xQBYoBy4bP0XQ7zirt81Wab+OZw8w79y0ArULv5c76nltfahgGL7zwAmXKlOH111+Xa5TCu154QX2ePh1OngSgXflOAOyIWpW/jYab3gdhVSH+PPzzmZsb6kP2FFit9qyl09Mwf7H6+q67XN5t2yVV+ahF6faebJ3HeCQoAY+spUyf7JGf3UBGNutOkLMuDn8nA7rZSf7tN/QqVTCdPYv/mDEEtG2L5ccfISEBLVHtS1dge5Q2G6ZVq7A+/zyBdergP2QI5sWL0ZxOnO3bk/rjj6Rs24Zz+HDIcBKhXTqEZfdvANg7PZ/np4tKSuD5f15CM6US5KzLhF5PuP1Xymj27NnMnTuXKVOmYLV6bp9LIbLVogXcdpuqNPPGGwA81X4Qhm7FYYlg1r4NeX8sawD0fld9vf4LiC4ilxHWfgCxpyG0AtS4B9aq/WxJKx+anQ2nDpJkOoxhaIxueaeXGupebg/KgIAAihUr5vagzLgMJD+7gZhMJj7o9BqGbiXRfJBXi8eSsmMHtjffxAgNxXTgAH5PPklgnTqYJv0LxxwYWu6lmLxFO3cO8y+/4Dd0KIFVqhDQvz/Wr79Gi4zEKFsW+9NPk7xjB6krVqhardn0sv1Wv4lm6Djq3o5eJW9ndIn2VO6b/wR2yxlwBjGx+7uZNl11txMnTvDEE0/w7bffUqNGDY89jxAuvfOO+vzLL7BpExWLlaSiRdUz/mXPjPw9VoP+ULOb2sx5wVOFv6xd+DZYP0F9ffvHsHi5+p1atlTF5XPw+SY1O7Y4jWlVKfeiJgWRRy42ubuMXW6Td3JzS+U6dCyhhl2XXJjEmgvHcTz3HMkHD2L76CP02rXR4uIwrT4FU5Pw6zUev//8B8uPP6Lt3Ane2jYsKQnT+vVYPv8cvxEjCKhXT/Ucx43DMm8eWkICRrlyOP7zH1JnziT58GHs77yD4WJ3FNPxVZhPrMYwWbF3fSVPzXDoTu6b+zzxpn0YupUXm71P47I5X6i/UXa7nQcffJDBgwdz773ZV/sRwivatoX0utCPPw66zrCGgwA4lbqBqKR8zGLVNBUo1iA4vgY2THB7c73Glgh/PQqGDk3uhQZ3qHWnAAMH5ni3+NRkDiSonYoG1R3sjZZ6hGZ4YHpqu3btGDNmDPfdd98NP1ZOy0Dyy+Z00GvmKJLNR9AcJfmt98/UKlku/UkwrV6N38fjMG0/CwmZ/0mMwED05s3RGzbEqF4do0YN9OrVMSpXhrCwPNU4xDAgMREtOhqiotAiIjAdPYp2/PiVz9rp02jXnHUaJhN6ixY4+/bF2acPRrNmarJAXjjtBPzcE1PUYextHsF+6+t5+ncaNu8VzjhWYxgmxtZ+g7Ete+ft+a6DYRi88cYbLFy4kC1bthAcHJz7nYTwpIgItSlxfDz8+COOUaNo9XMPdEsUt5Udx6d9H83f423/FeY9DpoZHlwCVdp6pt2eouswc6RaEhJcFsZvgiQdKlZUM14PHID69bO961urpzLz9IdojuJsHbm6UG3WnJFHWu2uHuWN9iQz8jNb+LnPpwxd/AC65RIPLH2S+fdMJiwgCEwm9B49MCLqQbdYbHVeQtsfh2nTJkzbtqHFxmLeuBHzxo3ZPrYRGIgRFqZC02pVf1iGoT7rurr+GR2Nlpp7zUi9QgX0Nm3QW7dWn1u0gNDr27vNsuVbTFGHMQJLYG//VK7Hx6Ykcd+8Z7is7cAwNAZWetqjIQmwdu1avv76azZs2CAhKQqG8uXh9dfhuefg6aex3Hor3SoMYlXkd6w4/zvxqaMI9c9HZaoW96se5d7ZMH0YjF52dZPnwmDthyokTVa4bwoElYT3XlEh2apVjiHpcDr56/hvYIEWJfoU2pAEDwalO5aI5LoMJJ+qFy/Lxx0/4ZkNj5BqOcE9f43jj/4TMi+e1zSMRjVxDLhDfa/raEeOYNq2DdOxY2gnTqCdPInpxAm0i2ryj5acjJacrM5Ec2H4+WGULAlly6LXqoVRqxZ6zZoYtWuj16qlXqRuoEUdxbr+EwBs3d+EANdTsk/GXOQ/i58g1XICQzdzb5Vneb5DzkMq7hAZGcmYMWP48MMPadasmUefS4h8efJJ+OsvtU/lsGG8tXwpq2dOxzDH8Pqqn/PXq9Q0uONzuHQYIvbAlAEqLEMKwez6rT+pCTwA/SdA1XYQEwMTJ6qfvfRSjnf94O8/sFvCQQ/gzW5jPd9WD/JYUEbkITRcyc8ykPzoWLU+T8W+zWf7XybRfIAB88Yype+XVC9eVq0LArVQOJ3JhFGvHs569chSm8PhgLg4tLg4iIlBi41VPzOZMn0YQUFQqpQKyOBgz1f9MXT8ljyD5kzFWeNWnGnXWHKy4NAW3t3+BrrlEjgDeabx2wxp3NmjTdR1nUcffZS2bdsybtw4jz6XEPlmscBvv0HTprBpE2EffET3TkNZefHb6+tVBhSD4bNg0m1w+QRMuQdG/AmhBXhfxs0/wKLn1Ncdn4QWadWLvvgC4uJUJZ677872rkn2VGYe+wEs0KbEAKqXLAQnBS54LCj37dt33fe/nmUg+TG0SWeCrZ/w7s4XsFlOMWzxaD7r/Cmd03pdWsagdMVigZIlVQCSdVcOX7Hs+AXz2S0Y1mBsvT7MMZhtTgfPLP+KzXHT0Sw6mqMkH7X/hC7VG3q8jd988w179uxhx44dHi1gIMR1q1oVfvgBBg+G99/n/Wm/cYszDMMcwwvLvuXr/vkr3EFoebj/L5jcGy7sVaF5/xwoVcBmghoGrP8cVryhvu/wOPR8U30dFweff66+fvnlHOdLvLpyctqJdygf9hzv6RZ7nEfeoW6kMPr1LgPJrzvrt+Hj9l+hOYrjtFzk8fWj+SnhMgBaaoxHntMbtKjDWNe+DYC968sYYVWyPW5XxCl6z3yALfG/o2k6JWnFjNuneiUkd+3axRtvvMGUKVMoXbq0x59PiOt2770wfjwYBoGjHuSRWFWAYN2lqew4dx1rI0vVUsOuJapDzCmY1AtO/O3eNt8IWyLMHn01JDs9Dbe9ffVk+8sv4fJlqFfvaiWja5yJiWLZObUkpEeF4ZQJKeaFhnuWx5aHREVF5bveqzsn7+RF52oNmNxjEoHOOmimVL4wHeG/ZUoRExPu0ef1GHsS/nMfRrMn46zWGUeLkVkOuZycyMOLPmbs2uEkmQ9j6H7cXu5xFg78kqphpTzexISEBEaNGsUzzzxD166FrziyuAl9/jn07w8pKTz6zg/UO1sOzWTj8WWv5K9aT7qSNWH0cqjQTBUW/6U/rHoHnF5ahpaTC/tVcO+drbYK6/s/6PH61ZA8cuTqOtNXXrlS9etaoxe8AuZ4zI6yvNt9jJca71keWR6ybds2unfvzpkzZ/Jchsxdy0Cuh83p4ImlE9iRMAM0KOGAzpWf4oWOg7GY3D/06yl+S57Fsvt3jOAyJI9aCRkqDOm6zif//snsUz9c2VMy0FmbDzu9wS2V63ilfYZhMH78eI4dO8bq1avdeu1ZCI9KTITu3WHzZpIrVWTA48UIL29hSLX/8nLF1nDsGNSuDXXy8VqyJcLiF2BH2nZVlVqpdZeVWnrmd8iJPUUVOV//OegOtQRk8C9QrcPVY5xO6NIFNmyAHj1g2bJsh10/Wz+HyUdfwzA0Xmv1NYObdPLe7+FBHgnK8PBwqlSpQmRkJP552KDz2p6kJ65L5sW8TdP55dhHhKeVT/NzVOOZ5s9wT4N2PmlPfpj3/4n/gvEYaKTe9wd6NTUZx+Z0MHHzPP46OQ2b5TQAmqMEA6uN5dl293j1hGTGjBk888wzbN++nerVq3vteYVwi4sXoVMnOHKEmOIhPPloGcbOi6bTvgxzGnr3hmnT8rS/6xV7ZsGCpyE1Tn3ffDh0fxWKeXgXI6dD9R7XfgDRasN66vWDfp9kfe6PP4bnn1dL1fbsgWrVsjzcyeiL9P/rbjDH0yDoTmbc+65n2+9FHgnK1NRUAgIC2L9/P5UrV871eLvdjmEYWK1W3/YynDZMn9VkamgwXxQvi2GyARCqN2J4veGMbHprgZx4okXsJmDa3Wj2ZOwdnsHe6XliU5L4YMN01lyYrS6qA4ZupUXo3Xxw66OUCPTumsXjx4/TqVMnfvjhB6m+IwqviAgVhrt34zCBZoA54zuo2Qw9e4KL7fyyFXcOVrwJu6er701WaDII2o2DCk3d1nxA7ZG590/Vg0wPyJDy0O9jVXbvWvv3qzJ1qalqctOYrMOpDqeTHlMfJJrtWBzlWTN8rlqjXkR4JCgNw6BkyZLMmzeP5s2buzzW4XCg6zoWi6VAFMIOmHwrpksHOdj9I544tYtzzn/QNPVPZHFUpHelQTze5m6vB01OtLiz+E/phynxAvbqXfi9wXhmHF3M6dSNYE5WBzmDaBral/+2H0nNkt6fjm6z2ejVqxctWrTg22+/lV1BROEWE6OGH7dvz/mYw4fzNwyb7swWWP4qnM5Q3KRSKxVgDe68/hmyqQlwagPsmQkH5oMj7b0hsKSa1dp2LPhnU9gkPh46d4Zdu6BPH1i0KNtZ9MNnv8HuhNkYupl3bvmWuxsW/FG4/PBYUNarV48PPviA2267LcfjnE4nTqcTs9ns0Rmu+WFd+RrWbT/gaDocW5+P2Rx+lE+3/szx1LVoaT1MQ7dSxtSCXlV7MqrZbb47c0pNwPL7nRyJPcbCsAr8HlAcp+XylZs1Rwk6lxnAix2GUSooxCdNNAyD119/ncWLF7N582apviOKhjlzXO6YQcuWsGJF/oZgMwrfChu/gv1zwciwgrt4VajYQn2UrqeKFgSXUes0dSc47WBPUjNqL5+EqGMqdM/tzPw4pWpDy5HQ+kHwz+G9wWaDfv3U71GmjDoxyGaE8M1VU5h1Ru3heVelZ3in5wPX9zsXYB4JSoAOHTowatQohg4dmu3t6dclTSaTV2a45pXp2EoCZo9AL1aZlIc3Xzl7Co+N4oONU9l6eQmGJfrK8YbuT0mtIY1KNKN7tbb0qNmEAEsear9eJ13X2XLuGPMPryP27C8csiYRm+GarqH7U9HSmrtq9WN4k64e3fEjL1avXs19993Hhg0baNrUzUNIQvjK4cNqiUROrncI9lrxEXBwIRxcACfWqck21yusKtTtDc2GqglDrkZ2dB1GjFDXW4ODYc0aaN06y2F/7P6bt7c9gWZy0DDoLv64953rb18B5rGgvOuuu2jbti1PPJF1D0NfznDNlS2JwC/qo+l2kh9YjVEmcx1DXddZcHgbMw4t4mjShiszSNMZuj+h1KRsQFWqh1ajfsnqtKpYh4alK+fr94xNSeJI1HmOxZxnz8WjHI45zIXUkyRrZ9FMKZmONel+FNMa0r1SDx5q2a/ADAtHRkbSrl07Xn75ZR577DFfN0cI9+rTB5Yvd7191vUOwWYnJVb1DM/vhHM7VI8xIRISL4JTjXahmcESAMWrQPFqar1mpZZqBmvxfOwA9Oyz8OmnqqjKggXquuw1/tr/L6/++ySYkyhBS1aNmIzFRxMxPc1jQTlmzBiKFSvGm2++mennBTok0/jNeQDLkSXYW43G3iPnMySH7mTh4W0sP7mRQ7G7iTOOgjkl22MNw4SmB2IyArEQjJ8pGDMWdJzohhMdJwZO7EYiTi02x8cBNXmgWWoKbZPtVKr/BF06POjRXuz10HWde++9l4CAAGbPnl0g/5+FuCGXL6teo6trlYsWQd++nm2HYaghV5Ml7zsL5UTX4b//hY/UUCq//gr335/lsGm71vLe9ufAlIK/swaLB08tEoUFcuKxcbnSpUtz9uzZTD/zdkGB6+Vo/h8sR5Zg2TsLe5eXwZp9TUeLycxd9dtyV321bY7N6WDViT38c2YXJ+NOcTElnHg9Aqf5EpqmgzkRnURsXMKWh3YYuhWzXoxgU3kqBlanQYna3HN+AU1OrsZi9iP17snotXq48Td3n6+//pq9e/eyc+fOAvv/LMQNKVFCDU26GoKtXdvz7dA0cMeJcmKiCsU5c9T3H3+cbUhO3rqMT/e8iGayE6TXYd7gn4p0SIIHg7Js2bLs3r0708/SdwPxVA1Xd9Grd0UPq4op9jTmfbNwNs/6x5IdP7OFPrVb0Kd2i0w/T7bbOHH5IhGJl4lMiuVSUiyXU2Kx606sJjNmzYLFZMZqslAysBjVw8pRp1RFygUXuxoyqQn4z38Y80m1CXPq3ZMKbEju3LmTN998k4ULF1KqlOer/QjhM3XrqmHJFSvUovwMdOCnH97mgQ9+Lvgni+fOwZ13wrZtan/dSZPUNcprvLV6KjNOfoJmchCqN2ThfZMo4aOJgt7ksaHXKVOm8MUXX7Bq1Sqg4C0DyY1l6/f4rXodI7gsyWPW5zwzzAu0+HP4z7ofU+R+DEsAqXd+j14759nEvpSQkECnTp0YMmQIb775piwFEUXf5cswdCgsXXrlR9HFAikZp5ZgLOzRiC4z1xJaooCeNK5bB8OGwdmzULq02l6sY8dMhySmpjLsz5c4blsGQHGas+C+74rUWklXPHaaU65cuSv1Xj29G4gnOFqMQi9eAy3xItZ/J/isHVrEbrVOMnI/RnAZUofOKbAhaRgGzz//POXKleO1116TkBQ3hxIl1OzWw4fVNcnDhyl+KZblA7oD0G/lPqIaVmP7vBk+bug1Ll+Ghx6Crl1VSDZsCJs2ZQnJnedP0vW3e6+EZJPggawc/tNNE5LgwaAsU6YMkZGRXtsNxO3MfthvfR1QvUvT+R3efX7DwLLtRwJ+648pIQK9dD1SRixEr9Dcu+3Ih5kzZzJ//nymTp0qdVzFzadOHTVxp04dTFYrt81eybyJ73ChuJXqEYk0GXAfywbciu06d1ZyG8OAGTOgQQNVaQfg4YdVHdeaNa8c5nA6eXHZD9y/+D5SzSfAGchDdd/h90Fv4HeTvb49NvR67tw5KlWqRHh4OIGBgQV2hqtLhoHfX6OxHFmMHlKBlJFLMxUa95jES/gvfgrz8ZUAOGvdRuodX4J/wb1gfvz4cTp27MikSZMYNMj1RtFC3Ex2791O1Ki7uHWb2pUoLsjK+fGPUu/NDyAwH5s/3yhdh8WL1YzWdevUz+rXh++/V9V3Mlh0aBuvr3+LFLMqcefnrMq3vSbQprIXJicVQB4LSpvNhr+/P0uWLKFjx46FLyTTpcYTMKUvpuhjOCu1IXXQVM8FlqFj3j8bvzVvoyVGYpj9sd/6Oo4Wo1wvDvax9BJ1LVu25JtvvpEhVyGuoTudTHnvGdp//wN1w9W1y8gSwSSPeoCqTz2vNon2lNRU+O03+OQTVbcV1ISdl16CF1+EDBtXHLgYzn9XTeBoylI0zcDQ/elSejgf9xlHkDX3DS6KKo8FZWpqKnXr1uX8+fP069ePxx57jHbt2hXKN1Et6ggBU/qh2eLRyzQkddBvGKHl3focpvM7sK58FfO5bQDopeuR2v/bLAUPChrDMHjttddYunQpmzZtkhJ1Qriw9+wJ5r35H0b9uYmKUXYAdA3iOnei+LjHoFev6y97l1FCgiqGsGABzJ8P6cO9oaFqmPXJJzOVo1t3Yh8fbPiW07a/0Uxq9m4ZrS1f9HqdxuU9GOKFhMeCEtSb6PHjx/niiy/4+eefqVu3LuPHj+eee+4pFDNfM9Iu7CFg1gi0xIvooRWx9foQvVbPG35c07ltWLZ8i+XQAgAMazD2Dk/haDUWLAX/DG7VqlUMGTKEjRs30qRJE183R4hCYeGuf9j69f/R++/dtDuQeOXnhqbhbNwYS7duav/HBg2gShUo5mIUy2ZTmyofPAgHDsD69bB6tepJpqtUCZ56CsaOhbAwAOJTk5m0bSlzjswhythxZfOHQGdtxjZ5hLFtslbjuVl5NCgziomJYdKkSXz55Zc4HA4efvhhHnjgAUq44+zJS7SY0/jPGoYp+hgAjjp9sXf6v/z3+lLjMR9bjmX75Cs9SABHo3uxd3nJ7b1VT7l48SLt27fn1VdfZdy4cb5ujhCFzow9/zBz8Yf0WL+J3ptjqRGRQymS0FDVAyxRAux2FY42GyQnw5kzWdZwAmpiTv/+6qNLF7BaiU1J4qfty1h0bAnnHdsgQznM4jRnfIuxDGnaxUO/beHltaBM53A4+PPPP5kwYQK7du1i+PDhjBs3jtq1axeOYVlbItYNn2LZ+j1aWoFiZ7mmOBsNRK/QEr10/cxrLg0DkqIwRR/BdHEf5mPLMZ3eiKarYRfD7Iez4QDsrR/CKNPAF7/RddF1nUGDBhEUFMSsWbMK7zVoIQqAreFH+WLzdM6cXEjLY2dpdSiRFkeSqHTJRrEkF7Vk0xUrpibmNGgATZqg9+7NucoV2Bd5mnUnd7D70h7OJR8m1RSuqoSl0Zxh1AxqzxNt7qd7Ldm0ICdeD8p0hmGwadMmPv/8c+bMmUPPnj157LHH6Ny5c6F409UiD2Jd/z/MR5dfCb10RlBpDM0MmoZmT0ZLjc1yf71kLZz178LeYpR3ZtK62cSJE/nqq6/YsWOHVN8Rwk0cTidTdq5i7pFlnErcjcNyjsAUJ+UuOyh32U5wso7drGE3W9FNYTgtViJKhRJdPAiTyYrDSMVmxKKb4tBM2e80ojmLUyuoPYMb9OPexp2KbCFzd/JZUKYzDIMzZ84wceJEJk2aROXKlRk/fjyDBw/G37/gX6MjKQrL/j8xn1iNFnkQU8L5LIcYaBhhVTBK1cZZtTPO2rdhlLzODVgLgB07dtCrVy8WLVpE165dfd0cIYqsI5fOM3PfGv49t4WLKadJ0i+im2OuXE/MlTOAIK0KVYPr07JcM3rVbk2LCjUKRWekIPF5UGYUHx/Pzz//zMSJE4mLi2Ps2LGMGTOGMmUKUY8rORot/rwacsUAsx9GWNUcC6sXNvHx8XTq1Ilhw4bxxhtvFI7hciGKkPjUZHacO87RqLOkOG3YnA5sDhspThvB1kCqhZWneoly1C5VgeIFZMu9wq5ABWU6p9PJggUL+Pzzz9m0aRODBw/mscceo0GDBvLG7EOGYfDoo49y8uRJVq1aJdV3hBA3hQLZ/zabzdx1112sWrWKf/75B8Mw6NSpE3fddRfLli1Dd7VRqvCYGTNmsGDBAqZMmSIhKYS4aRTIHmV2zp07x9dff813331H6dKlGTduHEOHDiUo6OYpzOtLx44do1OnTkyePJmBAwf6ujlCCOE1hSYo0yUlJTF16lQmTJjAhQsXGD16NA8//DDlyxeOtYeFkc1m47bbbqNVq1ZSok4IcdMpkEOvrgQFBfHQQw+xZ88epk6dyq5du2jYsCFjx45l165dFLLcL/AMw+Ctt94iOTmZzz77TEJSCHHTKXRBmc5kMtGnTx8WL17Mtm3bCA4OpkePHvTt25cFCxbgzK5Shci3VatW8f333zNt2jQCvbnTgRBCFBCFbujVlcjISL799lu++eYbgoKCePTRRxkxYgShoaG+blqhdPHiRdq1a8frr7/Oo48+6uvmCCGETxSpoEyXmprK9OnTmTBhAidOnGDkyJE88sgjVKlSRYYO80jXdQYOHEhISAgzZ86UBcpCiJtWkXz38/f3Z+TIkWzdupU5c+Zw/PhxmjVrxsiRI9m8ebNcx8yDL7/8koMHD/LDDz9ISAohbmpF+h3QZDLRrVs35s6dy549eyhXrhz9+/enR48ezJ49G7vdnvuD3IS2b9/O22+/zZQpUyhZsqSvmyOEED5VJIdeXbl8+TLff/89X331FZqm8cgjjzBy5EiKFy/u66YVCOkl6oYPH87rr78uQ9VCiJveTReU6ex2O7NmzWLChAns37+fESNG8Oijj1KzZs2bNhwMw+CRRx7h9OnTrFy5UqrvCCEERXzo1RWr1crQoUPZuHEjixcvJjIyktatWzNkyBD+/vvvm7JM3h9//MGiRYukRJ0QQmRw0/Yor2UYBidPnmTixIlMnjyZmjVrMm7cOAYNGoSfn5+vm+dx6SXqfvrpJwYMGODr5gghRIEhQZmNuLg4Jk+ezJdffklSUhIPPfQQo0ePLrIbFNtsNnr27EmbNm34+uuvb9qhZyGEyI4EpQsOh4N58+YxYcIEtm7dypAhQxg/fjz16tUrMmFiGAavvvoqy5cvZ/PmzVJ9RwghrnHTXqPMC4vFwoABA1izZg1r1qwhNTWVDh06MGDAAFauXFkkrmNKiTohhHBNgjIPNE2jTZs2TJ06laNHj9KyZUtGjRpFu3bt+OWXX0hOTvZ1E6/LhQsXGDt2LB9//DGNGzf2dXOEEKJAkqHX65SYmMgvv/zCF198QXR0NGPGjGHs2LGUK1fO103LE6fTycCBAylWrBgzZsyQ6jtCCJEDeXe8TsHBwYwbN479+/czefJktmzZQsOGDXnkkUfYs2dPgS+T9+WXX3Lo0CG+//57CUkhhHBB3iFvkMlk4o477mDZsmVs2rQJq9VKt27duOOOO1i0aFGB3O5r+/btvPPOO1KiTggh8kCGXj3gwoULfPPNN3z77bcUK1aMcePGMXz4cIKDg33dNOLj4+nYsSP3338/r732WpGZvSuEEJ4iQelBycnJ/P7770yYMIHw8HAeeOABHn74YSpVquSTgDIMg4cffpgzZ85IiTohhMgjGXr1oMDAQEaPHs3OnTuZMWMGBw8epEmTJjz44INs27bN69cxp0+fzuLFi6VEnRBC5IMEpReYTCZ69uzJggUL2LVrFyVKlKBv37706tWLv/76C4fD4fE2HD16lGeeeYYff/yRqlWrevz5hBCiqJChVx+Jioriu+++4+uvv8ZqtV7Z7qtYsWJuf67U1FR69uzJLbfccmV7MSGEEHkjQeljNpuNGTNmMGHCBA4fPsx//vMfHn30UapVq+aWQDMMg1deeYWVK1eyadMmqb4jhBD5JEOvPubn58eIESPYtGkT8+fP5+zZs7Ro0YIRI0awcePGG76OuXLlSn788UcpUSeEENdJgrKAMJlMdOnShdmzZ7N//36qVq3KgAED6NatG3/88Qc2my3fj5mxRF2jRo080GohhCj6ZOi1AIuJiWHSpEl8+eWXOBwOHnroIR544IE8FQlwOp0MGDCA4sWL88cff0j1HSGEuE4SlIWAw+Fgzpw5TJgwgZ07dzJ8+HDGjRtH7dq1c7yOOWHCBL777ju2b98u1XeEEOIGSDejELBYLNx77738/fffrFy5kri4ONq2bcvgwYNZu3Ztlu2+tm3bxttvvy0l6oQQwg2kR1kIGYbBmTNnmDhxIpMmTaJy5cqMHz+ewYMHk5qaSseOHRk5ciSvvvqqLAURQogbJEFZyMXHx/Pzzz8zceJE4uLiKFeuHKGhoaxZs0aq7wghhBtIUBYRTqeTWbNm8fjjj7N8+XKaNWvm6yYJIUSRIEFZxBiGIcOtQgjhRjKZp4iRkBRCCPeSoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCCFckKAUQgghXJCgFEIIIVyQoBRCCJE/uhMuHYGkaF+3xCssvm6AEEKIQiI+Av79Grb/CsmXwRIA9/4C9fr4umUepRmGYfi6ESIb9hSIOQ2XT0LsaXCkgu5QZ3IBxaBYJQitAKVqg3+Ir1srhCjKHDYVkGs/Anti1tvfiPV+m7xIepQFxaWjsPpdCC4DpzbAhb1AHs5hNDNUaApVO0CdnlCjK5jMHm+uEOImcX4XzHkELu5X31dqDV2eg+WvwaXDvm2bl0iP0pfsKXBkGez8DQ4vyXq7XwiUqAHFq4JfMJgsoJnUkEfcWfWRGJn5PsUqQ/Oh0PpBKFbRO7+HEKLo0Z3wz2ew5n01mhVUGnq9A82GqNvfLK4+V2gGD6/zWTO9QYLSF5Ivw+YfYdM3kBSV+bZyjaHzs1C1PYSWB01z/VgxZ+D0v3ByHeyfCylpQyCWAGg3Djo9BQFhHvk1hBBFVEIk/DkWjq9W3zfoD3d8DsGl1feHl8Hv94I1CJ7aC8GlfNZUb5Cg9KbkGHWGtmUS2OLVz0IrQtN7odE9ULHFjT2+PQUOLYRN38OZf9XPAktC/8+h4V039thCiJvD6X9hxkhIiABLIPT7GJoPv3rS7nTA913V5aEOT0Cvt33bXi+QoPQGXYfd09WYfvpQadlG0PkZaHg3mN18qdgw4NBiWPH61WsIrUdD73fBGuje5xJCFA2GAVsnweIX1FBrmfpw789QtkHm4zb/AIueg4Di8Pj2It+bBAlKz7t0BP4aB+Gb1fel68Jtb0HdPrkPq94op11NEPrnM/V9xZYwYjYElfTs8wohChdHqgq/7b+q7xvdA3d9peZGZBR3Dr5qB6mxcPvH0Has99vqAxKUnmIYsGOKOjuzJ6mJOV1fgFseAYufd9tydCXMHgPJ0eoa6P1/QUgZ77ZBCFEwJUbBHyPg9AY1WbDH69Dxyawn8oYB04aoiYcVW8Lo5e4fDSugJCg9ITUe5j0O++ao72t0hbu/gbBKvmvTxQPwy52QeBHKNIDRy9R6TCHEzSvyMPw+GC6fAP9iMOgntcwsOzunwV+PgMkKj/yddUi2CJMSdu4WexYm91UhabJAzzdUD86XIQnqj/qBxapIQeQB+OtRde1UCHFzOvE3TOqpQrJ4VXXynFNIXj4Ji55XX3d74aYKSZCgdK/zu+HHHnBhDwSXhQeWQKenwVRA/plL14b7fgOzHxxcAH9/4usWCSF8Yc8smDpALSer3BbGrMo5/JwO+PMhNVO/anvo9Ix321oAFJB38CLg9L/wU1+IP6+GNseuhCptfN2qrCq3gn6fqq/XvAcRe3zbHiGE9xgGrJ8As0eD0wYN7oSR813PWVjzPpzZpIZm7/nupqz8JUHpDme2wNRBYEuA6p3hwSVqKKOganm/WpZi6LDwWRmCFeJmoOuw5EW1TA2g3XhV0NwakPN9jq64OvJ0x2dQoprn21kASVDeqLPbYepANSxRvTMMmwGBxX3dqtz1fg+swepMcdfvvm6NEMKTHDZVaWfTt+r73u9Bn/dcXxaKO6eGXDFUScwmg7zS1IJIgvJGXD4Fvw1Sa4qqdoBhf4BfkK9blTdhldRFeYB1H6u6jkKIoseWqJZ17J2lJhgOnATtx7u+jyMVZvxHldgs3xR6v++dthZQEpTXKzUepg29+oc0fEbWxbkFXZsxqrrG5RNwaJGvWyOEcLekaLUs7NhKVZd12B956xku+S+Eb1F1ogfnMjx7E5CgvB66roYkLu6DkHIwdDr4h/q6VfnnFwxtRquvN0z0bVuEEO4VfwF+vgPOboXAEmrSTu0cln9ktGOqKmWHpnqfJWt6vKkFnQTl9djwheqBmf1hyO++XyN5I9o+pPa0PLMJok/4ujVCCHeIOaNm4V/cByHl1Rrqyq1zv9+ZLbDgafX1rS9Bnds8285CQoIyv87vglXvqK9v/1/e/vgKstDyUL2j+vrgAt+2RQhx4y4dhcl9IPoYhFWFBxfnrUBA3Dn4Y7haNlKvH3R+zvNtLSQkKPPDnqyGXHU71L8DWv7H1y1yj/r91ecDEpRCFGoXD6ieZFw4lKqjlqrlZejUngzTh0PCBSjbEAZ8V3AKpRQA8i+RH2s+gMiDqupO/wme3/3DW+r3U5/PbIKUON+2RQhxfSL2wM/9VD3nco3VcGteLgsZBsx9DM5tV9cyh/xeOOdceNDNUfrdHS4dgY1fqa/7T7i607e3GAZERcHp03DmDISHq4+ICLh8+epHbCzY7WrCkdOp7hcSAsWLq4+yZaF+fWjQAJo1U5/DKkGxyuosNGI3VO/k3d9NiILCMGDvXvVx+DCcOgUJCeojKQnMZvD3Vx8BAVCyJJQpc/WjfHmoXh0qVVLHesu5nTDlbki+DBWaw/1z8r6d3toPry4dGfwrlKzhwYYWThKUeWEYqqKFboc6vaD+7Z59rpMnYedO2L9fvVgPHVIfMTHX95hRUeoFn50aNWDAAKhcCwhXLzgJSnGzSUiA99+H337L+bWSHxYLVK2qXl81aqgT0kaN1EelSu4djQrfBlPuUeu5K7eB4bPyXvRkzyxVog5UacsaXdzXriJEgjIvDi9VpZzMftDnA/c9rmHAiROwcSNs3Qo7dqiAjI3N+T7ly0OVKlC5svpcoQKUKHH1IywM/PzU9YX0M9r4ePWYMTFw9iwcOKA+tm1Tz//JJ+BngQ5maLjdfb+fEIVBeDh06aJeCwCBgdCqFdSrp0IuLEyNygQGqlGa1FSw2SA5WZ2ERkZe/Th3To362O1w/Lj6uFaxYtCwITRurJ6nTRto0kS9bvPd9gwhWbU9DJ+Z92HT05vUpvIAHR6HViPz//w3CdmPMjeGAd91UUOSHZ+E2966/sdyOFQ4rV0LGzaogLx4MetxVqt6ETVurF6s6R+1a6sXq7skJsLSpfDNN7BihfpZm/Lw2zqoU8d9zyNEQTZyJPz6q+rpffop3HEHBN1AhS2nE86fV8F74gQcO6ZGh/btgyNH1PvAtfz8oGlTaN1afXTsqF7z2fU8Dx9WjxmSAv88qUKyWkdVPtM/JG9tjDoGP/ZUm7nX6wf3Tbkpi53nlQRlbg4vg9/vVXVRn96b93F/UCG7ezesXAmrVsG6dap3l5HVCi1bwi23qM/Nm6thmus5u7xeUVHQrjUcPXn1Z717w7RpqpcqRFHWvDns2gVTp8Lw4Z59LptNheW+feo5t25VH9HRWY8tUwY6d1YfXbqoEaT771cnt+lqmeHJbvDQX3kPycQotQ9l9HGo2AJGLSx8VcW8TILSFcOASbepUk4dHode7+R+n4QE1TtbtEh9nD2b+fYSJaBbN3XG2L69CscAH5eH6tNHtdmZod6r2Qw9e8KSJb5rlxDeMHgwzJwJ99wDs2d7fzZ7+ryELVtUaG7aBJs3Q0pK5uPMZjVJL+NbtkmDHj1g2fK8PZc9WZW0C9+sdjgasxJCyrrtVymqJChdObURfuqjKvA8tQdCy2V/XHQ0zJ0Ls2apwLHZrt4WGKiCsUcP6N5dDa94czZcbg4fVkM8rm6XYVhRlO3apU5YdR1eeEFN6vH10q/UVHWZZt06+Ptv9TkhIefjP/wQBg6EmjVzbrvuhJkj4cB8VcN19HIo4+K1L66QoHRlzqNqC6rmI+DurzLfFhengnH6dFi9OvN1h5o1oV8/9dG1q+97jK4sXgy3u5jFu2gR9O3rvfYI4Qvffw8PP6y+HjoUfvzxxq5TutuCBdC/f+7HVax4dai2Sxc1achkUr3QRc/Blh/VpMT758js9nyQoMxJShx8Ug/sSfDgMqh6iwrD5cvVhf+//so8NNKsmTqjGzhQXWP09RlpXkmPUgjl22/h8cfV67xuXZg0CToVkDD5dym075Pz7c2bq+uednvmn5cooYKzqgH6KihtVruBNLrbk60tcmR5SE72zlYhWboemCvDW2+ps86M1xzr11cX1++9t/CGSd26athp+zXLQtKvURbW30uI/HrkEdUDGzJEnSB26aJ6mW++qQp1+Er0cfjnCTVx54QT9Ay3ZZxLkJysrm2uW6c+NmxQRUjmzbt6fLkScHYB9Lary0El8zE58SYmPcqM0qdd164NG56BdavgeB34Z9/VodVSpdTQzMiRag1UYek55mTWLBg2LOuZqMx6FTery5fh+edVjxIgNBReekn1NoO9PDs09qwqcB57GoLrwLIwWLHq6u2uXqd2O8yZCD+8DMftcAawZ5iwZzJB27Zw111w993qxF9kS4IS1FnYo49m7lUFABknnXXsqI4ZNEiVr3IzhyOB8PDPOHfue2y2c/j5VaRixYeoXPlpLJY8TvvOj+RkePZZtYYSoF8PqPYvOErDc8ulJynEunXwzDNqUg2oXuVzz6keXESEOqE2jKsn1+5+zSREqgLnUUdUYfMHluAIDCbi71eI2z6V+HLROGtWyvl94uQ/MGUAOFOh+XC47X9qYtCyZWqJyf79mY+vV0/N/L37blUEQYqiX3FzB2V0tOpNZVyXdK0HH4Qnn1SzVT3E4Uhg586uJCTsJPO4iomQkOY0b77WvWG5Zg2MG6eq84B68d/fHv58AKrcAqOXue+5hCjMdF2VtXvjjeyr7GTkzlGY5Mvwc3+4sEfVYX5wMY6Qknl/nzi/S23anBoH9W6HwVPAfM2VtvBwWLhQzbdYuTLzqFLFimoIevhwaNGi8I+c3aAiHZQHI8MZMX8cqebsNyT+5uOTtNuXgMXFv0C/D+pwurz7e5AZ9Q67yO1hFzFl87eoG7AotixLY2/8Gkn5KBvPTY+g9xa1Q0hkmIWXx1ZmY+PrC2FDtzD7joXUK1PxhtsmhDdtO3uM0UsewWmJyNPxFofBH28cpU54KjlFhsME/zYM4dHnqrutnRl58n0iJMlJp93xdN8eT+fd8YSkXA3iE+X9WNi+OAvbhxFe1g3vhc4A3mn3JXc1vOXGH8tLinTf+qP1U3IMyWoRqXTa6zokAapetLk+wA06hkRn+8cPoKXdfiNKx9j5v9/PM//FI/TeEodTg+ndS3LPu7WvOyQBNJODD9f/ckNtE8IX3v772zyHJEClSzbqughJAIsOnfYmUDUi9cYbmA1Pvk8kBJlZ0q44/zeuCl0m1ueJJ6qypG0xUqwaNSJsPDbnIov/7wiT3z/O7Rtj8LPpuT9oTswpvLfpo+u/vw8U6Vmvz3e4n+EL/sZuPpPltk674rO5R1any3q+lFyYOZvaj2k0zfXtrpS9bGfU4kvcuzqaALs6I9hSL4gPRlTkcBX3rO18tv0ItzyOEN703w4PMXb5ZgxL3sKlSj5OmKtetHlkFMpT7xPXsltNrG5ZjNUtixGc7KT79jj6bYyl3b4E2hxKos2hJGKCzczrVJyZ3UpyskL+flfD0Pi/ts+6pa3eUqSHXrPQdZg/Hz76SE2ddkUDevX2Sgm3DRuqYLOF53i7n19lOnTIGvY52roVPvsMZsy4Olu3Qwd1naVnz6zXG5x2eK8iOG3wxE7Zj06Ia+W23jijhQtdF/HIjmHAwmdg62QwWWHYdKjdM9Mhub5P2P3osD4AwqrAg0sgrHL+2pCb8HA1E/jHH9XX6e64Q80S7ty5yF7LLNJDr1c4nSo0mjdXM7o2bFDFyCtWzHlmV9MK6sK8F1Ss+BA5/1eY0m7PRUqKam/nzmrG2u+/q5Ds0kXNcvvnH7jttuz/kC/sVSHpHwYlqt/AbyJEEVW3rpqsk5fyk/36qbKV06dnrdeakzUfqJBEgwHfZQlJyOV9woCK4RoEl4X/zHV/SILa2u/111Vd2vnzVUBqmqoa1LUrtGsHf/6ZuRZtEVH0g3LZMrXe8b77YM8etSbqv/9Vm7Pu3avCI6NaZWFMMHz8sNfWEFau/DQhIc3J+t+hZrNVrvx0znfevx+eflptETRsmApEqxVGjFA9y7Vrcw7IdEdXqs/V2hfZM0Ihbti0aWpEJietW6vwMJvV627oUHUyPn68ei3mFCBbJ8PatH1u+30MjQdme1iO7xMGhCSYqHyprCpNV6pWvn+1fDGb1e85f76aOf/QQ2rJ3ObNqjJZmzZXt+0rIoru0OuRI/DYYyooQW2++vTT8MQTWQPwyBE4elSthdr2BhxcALd/DG3Heq25+VpHeeEC/PGHmra+efPVn1epopazPPSQeoHm1aTecOZftcN5m9Hu+YWEKKoyvl/A1a/T11GGh8MPP8BPP8GZDJdM6tVT4Tl0qOqhAhxYADPuB0OHri/ArS+5fOos7xMOPyqeMah8sRSW+xeobbN84cIFmDgRJky4Wry9Rw/44gtV7aiQK3pB6XDABx/AO++oCvxWqwrMl19WVXVy8303OLcDhkyD+vm8zuBJly+rM7hp01S92fQtscxmVSx57Ni8Dw1llBQN/6ulXqhP7VFb7wghbpzTqfah/eknmDMn8zBsixbQow2kzoSSTmg5EvpPyPuIjsOmAvbwErVX7v1zVD1qX4uMhPfeg6+/VrsoWa3w4ouqslFB3hwiF0UrKE+fVmdr6RN1evVS/2G18jEU8WlDiDsLY1dDpZaeaWdeRUaq7btmz1ZDGRl3KLnlFrUY+L77bqwO5c7f4a9HoUwDGP/vjbdZCJFVXJxa2H/tiS5ApWIw4lFVSq5t29xPdh02mDkKDi0ESwAMnwU1Onuy9fl38qQavZs/X33fpIm6fpneCy9kik5QbtmiZppdugTFiqmAHDYs/9fcPqoFSZfg0Y1QzstDBoahdgBYsEB9bNyoZuqma9xYldAbNsx95bJ+uh1OrYfur0CX593zmEKInJ3YD/+9DbZdVEXOnRnegsuUUe9jffqouQXXjoI57SokDy5Q++QOnQa1e3i1+XlmGCocx49XQ7NhYWqD7GvnhRQCRSMot21Ts64SE9WQxqxZak/I6/FBNUiJgfFboExdtzYzWzExanhm+XK1FOXkycy3t2ihwnHgwLxPT8+rqGMwsSVoJnh6HxSTCjtCeJQtEX7upy7vlKwJA2fB2k2q57V4McTGXj1W09QEodtuU9f72raGhePhwDy1p+SQaVDHxeSiguLcObXD0oYNavh1+fKCs31ZHhX+oIyNhUaN1PZXt96qhipDQ6//8d6rBLYEeCLtD9ndkpJUT3HNGlVfcdOmzL1Gf3/1oujfX00zr1LF/W1It+IN+OczqNMLhs/03PMIIUB3wvThcHgxBJaEMSsyz1C122H9erUOc+lSNUs/I6sJKmlQwx9GvQ6DHitYm0u7YrPBgAHqdytdWk2IKl7c163Ks8IflG+/Da+9psa+t25V3fsbkd6jHLcJyrph25moKPj3X3U2tW6dCsZrt7SqX1+dNfbqpcLeG1v5pMbD501U8eXBU6DhnZ5/TiFuZotfhE3fqOuKI+dDlbaujz93Ts3aX7YMlvwFl5Mz326xqKVvHTuqOQu33AJVqxbcJV5JSaq9Bw+q/X1ffdXXLcqzwl/CLn29znPP3XhIAgSVUkGZFJX/+6akwO7dKrC3bFEBefBg1uMqV1YLkrt1UwFZ1QczTbdMUiFZqjbU7+f95xfiZrLlRxWSAPd8m3tIglriNWIY+C+DOhaIKQElhsOhS2q7rLNn1Yn3pk1X71O+vBqubdVKbcjeqpV6nIIQnkFBanb+s8+q98hCpPAHZWCg+rx3r7p4fKN/EMGlIfqYmtCTE8OAixfV0Mju3epj1y7VhowzU9PVq6dKyHXsqHqMNWr49g/XlgQbv1Rfd34WTPlcUiKEyLujK2HR/6mvu78Kje7J2/3sKTBzpFoCYgmAx3+/WrHHMFTRlL//Vifkmzap96CIiKuTAdOVKqW2CWzaVM0+bdBAjWKVLOne3zM38fFqmQyoQC9ECv/Q6+zZarILqBAaP1710ooVu77HmzZMTbvu+z+ocTecOKE+jh9X4+oHD6qPmJjs71+6tKpM0aqVKunUrl3e1m9604aJsOwVtWby8e1gtvq6RUIUTZGH4Meeal/IZkPh7m/ydpJsS4Lpw+D4ahWSQ6dBre6u75OcrDaf37bt6seBA5nnQGRUpowqfFCr1tWPqlXVvIhKldQayBtlt8POnWruyPffqyVvQUEq2Bs3vvHH95LCH5SgqkE8//zVa39ms/oDaNhQDXOWKwchIeo/yGRSZ2NOpxozT0xUa5yiotTSksNb4dxZSDCBw5nzc2qa+sNq1uzqmVqrVuqPrCAMc+Qk/gJ82Vq9cO+cCC3/4+sWCVE0JV+GH7pD9HGo2l7VYLXkYaeN1Hj4/T61bMsarAqk1+hynW1IVmUu00e/9u5VJ/oZKwZlR9NUkJYrp9ZplymjJt8UL646IQEBauKh1aqCWNfVc8XGquIo4eFqBv/+/arwS7rateHXX6F9++v7fXykaAQlqP+Ur75SZy5HjrjvcStVUktNatZUwVi/vvqoU6dwVpqY8wjsmqZKXY1ZKcOuQniC0wG/3wvHVqndPMauhpAyud8v+TJMHQRnt4J/MVVMwBMVdxIS4NAhVX7v2DH1cfy4CtAzZ9QsVXcJC4Pu3dUSkUGD3NNT9bKiE5QZnT2rzqAOHYLz59X1xIQE1YPUdXW2ZDarHmZwsFpOUrq0GiK1psLGl6FkCLx9qnCGYU5O/wuTe6uvx6yCyq182x4hiqqlL6t5ANYgeHApVGia+30SLsKUe9RuPoElVFk6X9RuTZ+DERGhCgVcvKiGTGNj1SWnuDjVS0xNVaN4JpN6P/X3V6EYFqY6GFWqqKV7NWvmvEtTIVE0g/JGOB3wYXWwxcPYVVCpiISJLRG+6wJRR6HF/XDXl75ukRBF0+6Z8OcY9fW9P+dt8k7Mafj1LjVMG1JOhWS5Rh5tpsi7wh3znmC2QO20i+aHFvu2Le605EUVkqEV4ba3fN0aIYqmiL0w73H1dedn8xaSkYdhch8VkmFV4YHFEpIFjARlduql7RpSVIJy/1zY/itXNoUN8vK0cCFuBskx8McIcCSrGaq3vpz7fcK3qsshcWehdD0YvdTz+0mKfJOgzE6dXqCZ1bWCiL2+bs2NuXTk6hlup6euf/acECJnhgF/jYPLJ9Syq4GTcp8od2Q5/NIfkqOhYkvVk5R6ywWSBGV2gkpCg/7q6/RqGoVRYhT8di+kxELlttDN9aawQojr9O83av212Q8G/5r7qM2u6TBtCNiToFYPVdIuuICttxZXSFDmpN049Xn3TEh0UaWnoHKkqgXLl09A8WpqwbLFz9etEqLoObsNlr+mvu79nuuZqoYB6/4Hcx4G3QFNBsPQ6eAf4p22iusiQZmTKm3VcIgzFf7+xNetyR+HDWaPhjP/gn+Y2hkkuLSvWyVE0ZMaD7MeBN0ODe6ENmNyPtbpgPlPwqp31PcdHod7vpMT2EJAgjInmqY2MwbY9B1cPODb9uSVI1Vt7HpgvhoGuu9XKOPmfSyFEMriF+HySVVU4M6JOVflSolTQ63bf1H7v97+MfR6p9CvL7xZyP+SK7V7QP07wHDCoudzrplYUNhT4I/7066V+KuNXWt283WrhCia9s+DnVMBTfUMA4tnf9zlkzCpFxxdDpZAuG8qtB3rxYaKGyVBmZve76mixCf/hn8K8BBswkX49U44slS1d9gfhWP3cyEKo6RoWPC0+rrTU1C9Y/bHnf5X1XuNPAChFeCBRbKtXSEkQZmbEtXUMAnAqnfh6Arftic753bC993gzKar1yRr3errVglRdC19WW3FV6Y+dPtv9sdsnQw/36H2tq3QLK3SV0vvtlO4hQRlXrS8H1qNAgx14f7sNl+3SDEM2DFVVfWIOwul6qgXo6yVFMJzTvwNu34HNHVd8todQRypMO8J1ePU7dDwLlkjWchJrde8cqSqWoynN3q2qn9exZ2HBU+pTV0Bat8GgyZBQJjv2iREUac74buucGEPtB4Nd3ya+fbo42oy3fldgAY9X4eOTxXsrfdEriQo8yM1IW2fuH/UPnF3TYTGA73bBt0JO39XGy+nxKiZrbe+BB2ekC2zhPC0XdPVGsiAMHh8R+YiAXtmqV5kahwEloSBP0BtmSdQFEhQ5pctCf4YrvaZA7UTR98PwS/Ys89rGHBkGax4Ay7uVz+r0Bzu+RbKNvDscwsh1Gvwmw7q9dfjNVX0HNTEnkXPw95Z6vsq7WDQZAir5Lu2CreSoLweTgeseT+tEIGhKv7f+hI0Hez+Xp3TAYcWqRJZpzeonwWEQefnoN2jYC58m6AKUSid3aZmsFqD4Zn96nW4azosfxUSI1V96M7PQtcX1C5EosiQoLwRx9fCX4+qiTQAZRtB+/HQ8E7wD72xx446Bntnw9afIP6c+pnZH9o9Ap2eVhu7CiG8Z/0EVaquUmt1krrhi7RrkaidP+7+RjZDL6IkKG+ULQk2fwf/fKaKj4NaVNzgDlXsuHJrKFnLdQUOw4DYcLiwD06shcNLIfrY1duDSkOrkWrygAznCOEb/34LS17I/DNrMHR9HtqNl1J0RZgEpbskX4Ytk9RQTNSRzLcFhEGJGhBUSn1Y/MCWqCYHJV+GyENgi898H5MFqnVU10Ab3pl1CroQwruSL8O0oWrme8ma0HgQ3PKI7PpxE5CgdDfDgLPbYf8cOLMFzu8ER0ru9zNZoHRdtSC5Tm9Vei6gmKdbK4TIL90pM8xvMhKUnua0q4Lq8edVhY6kKLUm0z8U/EJUGJaqrYZnZehGCCEKHAlKIYQQwgUpYSeEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrggQSmEEEK4IEEphBBCuCBBKYQQQrjw/4rTAK57caYYAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure()\n", + "\n", + "ax = fig.add_subplot(2,2,1,projection=\"3d\")\n", + "ax = udp.plot(pop.champion_x, ax=ax)\n", + "ax.view_init(45,30)\n", + "ax.get_xaxis().set_ticks([])\n", + "ax.get_yaxis().set_ticks([])\n", + "ax.get_zaxis().set_ticks([])\n", + "\n", + "ax = fig.add_subplot(2,2,2,projection=\"3d\")\n", + "ax = udp.plot(pop.champion_x, ax=ax)\n", + "ax.view_init(90,0)\n", + "ax.axis('off')\n", + "\n", + "ax = fig.add_subplot(2,2,3,projection=\"3d\")\n", + "ax = udp.plot(pop.champion_x, ax=ax)\n", + "ax.view_init(0,0)\n", + "ax.axis('off')\n", + "\n", + "ax = fig.add_subplot(2,2,4,projection=\"3d\")\n", + "ax = udp.plot(pop.champion_x, ax=ax)\n", + "ax.view_init(0,90)\n", + "ax.axis('off')\n", + "\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "kep3_devel", + "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.11.10" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/kep3_devel.yml b/kep3_devel.yml index 97b27605..380121f9 100644 --- a/kep3_devel.yml +++ b/kep3_devel.yml @@ -9,8 +9,8 @@ dependencies: - cmake >=3.18 - boost-cpp >=1.73 - fmt - - heyoka = 5* - - heyoka.py = 5* + - heyoka >= 5* + - heyoka.py >= 5* - spdlog - xtensor - xtensor-blas diff --git a/pykep/docstrings.cpp b/pykep/docstrings.cpp index 9f1cf996..223f0382 100644 --- a/pykep/docstrings.cpp +++ b/pykep/docstrings.cpp @@ -1012,7 +1012,7 @@ std::string planet_elements_docstring() { return R"(elements(when = 0., el_type = KEP_F) -The period of the planet in seconds. +The elements of the planet at epoch. If the UDPLA provides a ``elements(float, pk.el_type)`` method, then ``planet.elements`` will call it. Otherwise, if the UDPLA provides a ``get_mu_self()`` method ``planet.elements`` will return the elements as computed by the diff --git a/pykep/trajopt/CMakeLists.txt b/pykep/trajopt/CMakeLists.txt index 4c52977b..abf279a6 100644 --- a/pykep/trajopt/CMakeLists.txt +++ b/pykep/trajopt/CMakeLists.txt @@ -1,2 +1,2 @@ -set(PYKEP_TRAJOPT_PYTHON_FILES __init__.py _direct_point2point.py _direct_pl2pl.py) +set(PYKEP_TRAJOPT_PYTHON_FILES __init__.py _direct_point2point.py _direct_pl2pl.py _mga.py) install(FILES ${PYKEP_TRAJOPT_PYTHON_FILES} DESTINATION ${_PYKEP_INSTALL_DIR}/trajopt) diff --git a/pykep/trajopt/__init__.py b/pykep/trajopt/__init__.py index 09cb7c62..b72864d8 100644 --- a/pykep/trajopt/__init__.py +++ b/pykep/trajopt/__init__.py @@ -9,6 +9,9 @@ User defined problems (compatible to pagmo) that represent interplanetary optimization problems """ -# Direct methods +# Direct methods for low-thrust problems from ._direct_point2point import direct_point2point -from ._direct_pl2pl import direct_pl2pl \ No newline at end of file +from ._direct_pl2pl import direct_pl2pl + +# Evolutionary encodings for high energy transfers (chemical propulsion) +from ._mga import mga \ No newline at end of file diff --git a/pykep/trajopt/_mga.py b/pykep/trajopt/_mga.py new file mode 100644 index 00000000..0bf49a7f --- /dev/null +++ b/pykep/trajopt/_mga.py @@ -0,0 +1,484 @@ +import pykep as _pk +import numpy as _np + +from typing import Any, Dict, List, Tuple +from bisect import bisect_left + + +class mga: + r""" + The Multiple Gravity Assist (MGA) encoding of an interplanetary trajectory. + + This class may be used as a User Defined Problem (UDP) for the pygmo (http://esa.github.io/pygmo/) optimisation suite. + + - Izzo, Dario. "Global optimization and space pruning for spacecraft trajectory design." Spacecraft Trajectory Optimization 1 (2010): 178-200. + + The decision vector (chromosome) is:: + + direct encoding: [t0, T1, T2 ... ] in [mjd2000, days, days ... ] + alpha encoding: [t0, T, a1, a2 ...] in [mjd2000, days, nd, nd ... ] + eta encoding: [t0, n1, n2, n3 ...] in [mjd2000, nd, nd ...] + + .. note:: + + The time of flights of a MGA trajectory (and in general) can be encoded in different ways. + When they are directly present in the decision vector, we have the *direct* encoding. This is the most 'evolvable' encoding + but also the one that requires the most problem knowledge (e.g. to define the bounds on each leg) and is not + very flexible in dealing with constraints on the total time of flight. The *alpha* and *eta* encodings, instead, allow + to only specify bounds on the time of flight of the entire trajectory, and not on the single legs: a property that is attractive + for multi-objective optimization, for example. + + In the *alpha* encoding each leg time-of-flight is decoded as follows, T_i = T log(alpha_i) / \sum_n(log(alpha_n)). + In the *eta* encoding each leg time-of-flight is decoded as follows, T_i = (tof_max - \sum_0^(i-1)(T_j)) * eta_i + + The chromosome dimension for the direct and eta encoding is the same, while the alpha encoding requires one more gene. + + .. note:: + + The resulting problem is box-bounded (unconstrained). + """ + + def __init__( + self, + seq=[ + _pk.planet(_pk.udpla.jpl_lp("earth")), + _pk.planet(_pk.udpla.jpl_lp("venus")), + _pk.planet(_pk.udpla.jpl_lp("earth")), + ], + t0=[0, 1000], + tof=[[30, 200], [200, 300]], + vinf=2.5, + multi_objective=False, + tof_encoding="direct", + orbit_insertion=False, + e_target=None, + rp_target=None, + ): + r"""mga(seq, t0, tof, vinf, multi_objective=False, alpha_encoding=False, orbit_insertion=False, e_target=None, rp_target=None) + + Args: + - seq (``list of pk.planet``): sequence of body encounters including the departure planet. + - t0 (:class:`list` of :class:`float` or :class:`~pk.epoch`): lower and upper bounds for the launch epoch. When floats are used MJD2000 is assumed. + - tof (``list`` or ``float``): defines the bounds on the time of flight. If *tof_encoding* is 'direct', this contains a list + of 2D lists defining the upper and lower bounds on each leg. If *tof_encoding* is 'alpha', + this contains a 2D list with the lower and upper bounds on the total time-of-flight. If *tof_encoding* + is 'eta' tof is a float defining an upper bound for the time-of-flight. + - vinf (:class:`float`): the vinf provided at launch for free + - multi_objective (:class:`bool`): when True constructs a multiobjective problem (dv, T). In this case, 'alpha' or `eta` encodings are recommended + - tof_encoding (:class:`str`): one of 'direct', 'alpha' or 'eta'. Selects the encoding for the time of flights + - orbit_insertion (:class:`bool`): when True the arrival dv is computed as that required to acquire a target orbit defined by e_target and rp_target + - e_target (:class:`float`): if orbit_insertion is True this defines the target orbit eccentricity around the final planet + - rp_target (:class:`float`): if orbit_insertion is True this defines the target orbit pericenter around the final planet (in m) + - max_revs (:class:`int`): maximal number of revolutions for lambert transfer + + Raises: + - ValueError: if *planets* do not share the same central body (checked on the mu_central_body attribute) + - ValueError: if *t0* does not contain objects able to construct a epoch (e.g. pk. epoch or floats) + - ValueError: if *tof* is badly defined + - ValueError: it the target orbit is not defined and *orbit_insertion* is True + """ + + # Sanity checks + # 1 - All planets need to have the same mu_central_body + if [r.get_mu_central_body() for r in seq].count( + seq[0].get_mu_central_body() + ) != len(seq): + raise ValueError( + "All planets in the sequence need to have exactly the same mu_central_body" + ) + + # 2 - We try to build epochs out of the t0 list (mjd2000 by default) + for i in range(len(t0)): + if type(t0[i]) != type(_pk.epoch(0.0)): + t0[i] = _pk.epoch(t0[i], _pk.epoch.julian_type.MJD2000) + + # 3 - Check the tof bounds + if tof_encoding == "alpha": + if len(tof) != 2: + raise TypeError( + r"When the tof_encoding is 'alpha', the tof must be in the form [lb, ub]." + ) + elif tof_encoding == "direct": + if type(tof) is not type([]): + raise TypeError( + r"The selected encoding is 'direct', the tof must be of type list." + ) + if len(tof) != (len(seq) - 1): + raise TypeError( + r"The selected encoding is 'direct', the length of the tof list must be exactly " + + str(len(seq) - 1) + + r" while it seems to be " + + str(len(tof)) + ) + if not all(isinstance(elem, list) and len(elem) == 2 for elem in tof): + raise TypeError( + r"The selected encoding is 'direct', the content of the tof list must be 2-D lists." + ) + if not all(elem[1] >= elem[0] for elem in tof): + raise TypeError( + r"The selected encoding is 'direct', the content of the tof list must contain lower bounds < upper bounds " + ) + elif tof_encoding == "eta": + try: + float(tof) + except TypeError: + raise TypeError( + r"When tof_encoding is \'eta\', the tof must be a float." + ) + if not tof_encoding in ["alpha", "eta", "direct"]: + raise ValueError("tof_encoding must be one of 'alpha', 'eta', 'direct'") + + # 4 - Check that if orbit insertion is selected e_target and r_p are + # defined + if orbit_insertion: + if rp_target is None: + raise ValueError( + "The rp_target needs to be specified when orbit insertion is selected" + ) + if e_target is None: + raise ValueError( + "The e_target needs to be specified when orbit insertion is selected" + ) + + # Public data members + self.seq = seq + self.t0 = t0 + self.tof = tof + self.vinf = vinf * 1000 + self.multi_objective = multi_objective + self.tof_encoding = tof_encoding + self.orbit_insertion = orbit_insertion + self.e_target = e_target + self.rp_target = rp_target + + # Private data members + self._n_legs = len(seq) - 1 + self._common_mu = seq[0].get_mu_central_body() + + def get_nobj(self): + return self.multi_objective + 1 + + def get_bounds(self): + t0 = self.t0 + tof = self.tof + n_legs = self._n_legs + + if self.tof_encoding == "alpha": + # decision vector is [t0, T, a1, a2, ....] + lb = [t0[0].mjd2000, tof[0]] + [1e-3] * (n_legs) + ub = [t0[1].mjd2000, tof[1]] + [1.0 - 1e-3] * (n_legs) + elif self.tof_encoding == "direct": + # decision vector is [t0, T1, T2, T3, ... ] + lb = [t0[0].mjd2000] + [it[0] for it in self.tof] + ub = [t0[1].mjd2000] + [it[1] for it in self.tof] + elif self.tof_encoding == "eta": + # decision vector is [t0, n1, n2, ....] + lb = [t0[0].mjd2000] + [1e-3] * (n_legs) + ub = [t0[1].mjd2000] + [1.0 - 1e-3] * (n_legs) + return (lb, ub) + + def _decode_tofs(self, x: List[float]) -> List[float]: + if self.tof_encoding == "alpha": + # decision vector is [t0, T, a1, a2, ....] + T = _np.log(x[2:]) + return T / sum(T) * x[1] + elif self.tof_encoding == "direct": + # decision vector is [t0, T1, T2, T3, ... ] + return x[1:] + elif self.tof_encoding == "eta": + # decision vector is [t0, n1, n2, n3, ... ] + dt = self.tof + T = [0] * self._n_legs + T[0] = dt * x[1] + for i in range(1, len(T)): + T[i] = (dt - sum(T[:i])) * x[i + 1] + return T + + def alpha2direct(self, x): + """alpha2direct(x) + + Args: + - x (``array-like``): a chromosome encoding an MGA trajectory in the alpha encoding + + Returns: + ``numpy.array``: a chromosome encoding the MGA trajectory using the direct encoding + """ + T = _np.log(x[2:]) + retval = T / sum(T) * x[1] + retval = _np.insert(retval, 0, x[0]) + return retval + + def direct2alpha(self, x): + """direct2alpha(x) + + Args: + - x (``array-like``): a chromosome encoding an MGA trajectory in the direct encoding + + Returns: + ``numpy.array``: a chromosome encoding the MGA trajectory using the alpha encoding + """ + T = _np.sum(x[1:]) + alphas = _np.exp(x[1:] / (-T)) + retval = _np.insert(alphas, 0, [x[0], T]) + return retval + + def eta2direct(self, x): + """eta2direct(x) + + Args: + - x (``array-like``): a chromosome encoding an MGA trajectory in the eta encoding + + Returns: + ``numpy.array``: a chromosome encoding the MGA trajectory using the direct encoding + + Raises: + - ValueError: when the tof_encoding is not 'eta' + """ + if self.tof_encoding != "eta": + raise ValueError("cannot call this method if the tof_encoding is not 'eta'") + + # decision vector is [t0, n1, n2, n3, ... ] + n = len(x) - 1 + dt = self.tof + T = [0] * n + T[0] = dt * x[1] + for i in range(1, len(T)): + T[i] = (dt - sum(T[:i])) * x[i + 1] + T = _np.insert(T, 0, x[0]) + return T + + def direct2eta(self, x): + """direct2eta(x) + + Args: + - x (``array-like``): a chromosome encoding an MGA trajectory in the direct encoding + + Returns: + ``numpy.array``: a chromosome encoding the MGA trajectory using the eta encoding + + Raises: + - ValueError: when the tof_encoding is not 'eta' + """ + if self.tof_encoding != "eta": + raise ValueError("cannot call this method if the tof_encoding is not 'eta'") + from copy import deepcopy + + retval = deepcopy(x) + retval[1] = x[1] / self.tof + for i in range(2, len(x)): + retval[i] = x[i] / (self.tof - sum(x[1:i])) + return retval + + def _compute_dvs(self, x: List[float]) -> Tuple[ + float, # DVlaunch + List[float], # DVs + float, # DVarrival, + List[Any], # Lambert legs + float, # DVlaunch_tot + List[float], # T + List[Tuple[List[float], List[float]]], # ballistic legs + List[float], # epochs of ballistic legs + ]: + # 1 - we 'decode' the times of flights and compute all epochs (mjd2000) + # of the various planetary encounters + T: List[float] = self._decode_tofs(x) # T is [T1, T2 ...] + ep = _np.insert(T, 0, x[0]) # T is [t0, T1, T2 ...] + ep = _np.cumsum(ep) # [t0, t1, t2, ...] + + # 2 - we compute the ephemerides + r = [0] * len(self.seq) + v = [0] * len(self.seq) + for i in range(len(self.seq)): + r[i], v[i] = self.seq[i].eph(ep[i]) + + # 3 - we solve the lambert problems (and store trajectory r,v) + l = list() + for i in range(self._n_legs): + lp = _pk.lambert_problem( + r0=r[i], + r1=r[i + 1], + tof=T[i] * _pk.DAY2SEC, + mu=self._common_mu, + cw=False, + multi_revs=0, + ) + l.append(lp) + + # 4 - we compute the various dVs needed at fly-bys to match incoming + # and outcoming + DVfb = list() + for i in range(len(l) - 1): + v_rel_in = [a - b for a, b in zip(l[i].v1[0], v[i + 1])] + v_rel_out = [a - b for a, b in zip(l[i + 1].v0[0], v[i + 1])] + DVfb.append(_pk.fb_dv(v_rel_in, v_rel_out, self.seq[i + 1])) + + # 5 - we add the departure and arrival dVs + DVlaunch_tot = _np.linalg.norm([a - b for a, b in zip(v[0], l[0].v0[0])]) + DVlaunch = max(0, DVlaunch_tot - self.vinf) + DVarrival = _np.linalg.norm([a - b for a, b in zip(v[-1], l[-1].v1[0])]) + if self.orbit_insertion: + # In this case we compute the insertion DV as a single pericenter + # burn + DVper = _np.sqrt( + DVarrival * DVarrival + 2 * self.seq[-1].mu_self / self.rp_target + ) + DVper2 = _np.sqrt( + 2 * self.seq[-1].mu_self / self.rp_target + - self.seq[-1].mu_self / self.rp_target * (1.0 - self.e_target) + ) + DVarrival = _np.abs(DVper - DVper2) + return (DVlaunch, DVfb, DVarrival, l, DVlaunch_tot, ep, T) + + # Objective function + def fitness(self, x): + DVlaunch, DVfb, DVarrival, _, _, _, _ = self._compute_dvs(x) + if self.tof_encoding == "direct": + T = sum(x[1:]) + elif self.tof_encoding == "alpha": + T = x[1] + elif self.tof_encoding == "eta": + T = sum(self.eta2direct(x)[1:]) + if self.multi_objective: + return [DVlaunch + _np.sum(DVfb) + DVarrival, T] + else: + return [DVlaunch + _np.sum(DVfb) + DVarrival] + + def to_planet(self, x: List[float]): + """ + For a decision vector *x*, constructs a virtual planet with the spacecraft ephemerides. + + Args: + - x (``list``, ``tuple``, ``numpy.ndarray``): Decision chromosome, e.g. (``pygmo.population.champion_x``). + + Example: + udpla = mga_udp.to_planet(population.champion_x) + r, v = udpla.eph(7000.) + """ + + class mga_udpla: + def __init__(self, seq, lambert_legs, mjd2000s): + self.seq = seq # p0, p1, .., pn + self.lambert_legs = lambert_legs # l0, l1, .., l(n-1) + self.mjd2000s = mjd2000s # ep1, ep2, .., epn + + def eph(self, mjd2000: float): + if mjd2000 < self.mjd2000s[0]: + raise ValueError( + "Ephemerides out-of-bounds, requested at mjd2000 " + + str(mjd2000) + + " which is before the launch date " + + str(self.mjd2000s[0]) + ) + if mjd2000 > self.mjd2000s[-1]: + raise ValueError( + "Ephemerides out-of-bounds, requested at mjd2000 " + + str(mjd2000) + + " which is after the arrival date " + + str(self.mjd2000s[0]) + ) + + if mjd2000 == self.mjd2000s[0]: + # exactly at launch + return self.lambert_legs[0].r0, self.lambert_legs[0].v0[0] + + i = bisect_left( + self.mjd2000s, mjd2000 + ) # finds the corresponding leg from seq[i-1] to seq[i] + + # the starting conditions of the leg + r0, v0 = self.lambert_legs[i - 1].r0, self.lambert_legs[i - 1].v0[0] + elapsed_seconds = (mjd2000 - mjd2000s[i - 1]) * _pk.DAY2SEC + + # propagate ballistically the starting conditions + r1, v1 = _pk.propagate_lagrangian( + rv=[r0, v0], + tof=elapsed_seconds, + mu=self.seq[0].get_mu_central_body(), + stm=False, + ) + return r1, v1 + + def get_name(self): + return "MGA spacecraft" + + if len(x) != len(self.get_bounds()[0]): + raise ValueError( + "Expected chromosome of length " + + str(len(self.get_bounds()[0])) + + " but got length " + + str(len(x)) + ) + + _, _, _, lambert_legs, _, mjd2000s, _ = self._compute_dvs(x) + return _pk.planet(mga_udpla(self.seq, lambert_legs, mjd2000s)) + + def pretty(self, x): + """pretty(x) + + Args: + - x (``list``, ``tuple``, ``numpy.ndarray``): Decision chromosome, e.g. (``pygmo.population.champion_x``). + + Prints human readable information on the trajectory represented by the decision vector x + """ + DVlaunch, DVfb, DVarrival, lambert_legs, DVlaunch_tot, mjd2000s, T = ( + self._compute_dvs(x) + ) + print("Multiple Gravity Assist (MGA) problem: ") + print("\tPlanet sequence: ", [pl.get_name() for pl in self.seq]) + print("\tEncoding for tofs: ", self.tof_encoding) + print("\tOrbit Insertion: ", self.orbit_insertion) + + print("Departure: ", self.seq[0].get_name()) + print("\tEpoch: ", mjd2000s[0], " [mjd2000]") + print("\tSpacecraft velocity: ", lambert_legs[0].v0[0], "[m/s]") + print("\tHyperbolic velocity: ", DVlaunch_tot, "[m/s]") + print("\tInitial DV: ", DVlaunch, "[m/s]") + + for pl, e, dv in zip(self.seq[1:-1], mjd2000s[1:-1], DVfb): + print("Fly-by: ", pl.get_name()) + print("\tEpoch: ", e, " [mjd2000]") + print("\tDV: ", dv, "[m/s]") + + print("Arrival: ", self.seq[-1].get_name()) + print("\tEpoch: ", mjd2000s[-1], " [mjd2000]") + print("\tSpacecraft velocity: ", lambert_legs[-1].v1[0], "[m/s]") + print("\tArrival DV: ", DVarrival, "[m/s]") + + print("Time of flights: ", T, "[days]") + + print("\nTotal DV: ", DVlaunch + _np.sum(DVfb) + DVarrival) + + def plot(self, x, ax=None, units=_pk.AU, N=60, figsize=(5, 5)): + """plot(self, x, axes=None, units=pk.AU, N=60) + + Plots the spacecraft trajectory. + + Args: + - x (``tuple``, ``list``, ``numpy.ndarray``): Decision chromosome. + - axes (``matplotlib.axes._subplots.Axes3DSubplot``): 3D axes to use for the plot + - units (``float``, ``int``): Length unit by which to normalise data. + - N (``float``): Number of points to plot per leg + """ + import matplotlib.pyplot as plt + + if ax is None: + ax = _pk.plot.make_3Daxis(figsize=figsize) + _, _, _, _, _, mjd2000s, _ = self._compute_dvs(x) + for i, item in enumerate(self.seq): + _pk.plot.add_planet(pla=item, ax=ax, when=mjd2000s[i], c="r", units=units) + _pk.plot.add_planet_orbit(pla=item, ax=ax, units=units, N=N) + + _pk.plot.add_sun(ax=ax) + + _pk.plot.add_planet_orbit( + pla=self.to_planet(x), + ax=ax, + plot_range=[mjd2000s[0], mjd2000s[-1]], + c="r", + units=units, + N=N, + ) + return ax + + +del Any, Dict, List, Tuple diff --git a/src/core_astro/ic2eq2ic.cpp b/src/core_astro/ic2eq2ic.cpp index 333e4d4a..df6a6fa2 100644 --- a/src/core_astro/ic2eq2ic.cpp +++ b/src/core_astro/ic2eq2ic.cpp @@ -29,7 +29,7 @@ namespace kep3 // Cefola: Equinoctial orbit elements - Application to artificial satellite // orbitsCefola, P., 1972, September. Equinoctial orbit elements-Application to // artificial satellite orbits. In Astrodynamics Conference (p. 937). - + std::array ic2eq(const std::array, 2> &pos_vel, double mu, bool retrogade) { {