Skip to content

Instantly share code, notes, and snippets.

@hannorein

hannorein/masses.txt

Created Nov 12, 2018
Embed
What would you like to do?
1.000000000000000000e+00
1.660114153054348795e-07
2.447838287784771540e-06
3.040432648022642030e-06
3.227156037554997664e-07
9.547919152112404250e-04
2.858856727222416731e-04
4.366243735831270249e-05
5.151383772628674303e-05
7.361781606089468627e-09
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"y0 = np.loadtxt(\"nbody.txt\")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"m = np.loadtxt(\"masses.txt\")"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"Nd = len(y0)\n",
"N = Nd//6"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"#. 0. 1. 2. 3. 4. 5. 6. 7. \n",
"# y = (x_0, y_0, z_0, vx_0, vy_0, vz_0, x_1, y_1, ...)\n",
"def G(y):\n",
" r = np.zeros(Nd)\n",
" for i in range(N):\n",
" # Right hand side of positions\n",
" # Time derivative of position = velocity\n",
" r[i*6+0] = y[i*6+3]\n",
" r[i*6+1] = y[i*6+4]\n",
" r[i*6+2] = y[i*6+5]\n",
" \n",
" # Right hand side of velocities\n",
" # Time derivative of velocity = acceleration\n",
" for j in range(N):\n",
" if i!=j:\n",
" dx = y[i*6+0]-y[j*6+0]\n",
" dy = y[i*6+1]-y[j*6+1]\n",
" dz = y[i*6+2]-y[j*6+2]\n",
" \n",
" d = np.sqrt(dx**2 + dy**2 + dz**2)\n",
" \n",
" r[i*6+3] += -m[j]/(d**3)*dx\n",
" r[i*6+4] += -m[j]/(d**3)*dy\n",
" r[i*6+5] += -m[j]/(d**3)*dz\n",
" \n",
" return r"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"def euler_step(y,dt):\n",
" return y + dt* G(y)"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"steps = 20000\n",
"dt = 0.0005\n",
"y = y0\n",
"orbits = np.zeros((steps,Nd))\n",
"for i in range(steps):\n",
" y = euler_step(y,dt)\n",
" orbits[i] = y"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x10bcecb38>]"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAD8CAYAAACPd+p5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXeYW9d1r/3ugw4MML33xt6rSJHqvVuWbMuxJTuW5bg7cYpj59qJk3sdlyRucveXSLJiq1ldVBclURQp9j6d0xumAAMMOs7+/jiYESVLGsocEjPD/T7PPATOHJ69Dgb4Ye21115LSClRKBSK90JLtwEKhWLmo4RCoVBMiRIKhUIxJUooFArFlCihUCgUU6KEQqFQTIkSCoVCMSVKKBQKxZQooVAoFFNiTrcB70ZeXp6sqqpKtxkKxZxmz549Q1LK/KnOm7FCUVVVxe7du9NthkIxpxFCdJzMeWrqoVAopkQJhUKhmBIlFAqFYkqUUCgUiilRQqFQKKZECYVCoZgSJRQKhWJKlFAoFIopUUKhUCimRAmFQqGYEiUUCoViSpRQKBSKKVFCoVAopkQJhUKhmBIlFAqFYkqUUCgUiimZFqEQQlwhhGgUQrQIIb72Hud9UAghhRBrpmNchUJxZjhloRBCmIA7gSuBRcAtQohF73CeG/gysPNUx1QoFGeW6fAo1gEtUso2KWUM+ANw/Tuc96/Ad4HINIypUCjOINMhFKVA1wnPu1PHJhFCrALKpZRPTsN4CoXiDHPag5lCCA34T+CrJ3HuHUKI3UKI3V6v93SbplAoTpLpEIoeoPyE52WpYxO4gSXAViFEO3AO8Ng7BTSllL+SUq6RUq7Jz5+ygrhCoThDTIdQ7ALqhRDVQggr8BHgsYlfSin9Uso8KWWVlLIK2AFcJ6VUtfgVilnCKQuFlDIBfAF4BjgG3C+lPCKE+LYQ4rpTvb5CoUg/09IASEr5FPDU2459813OvWA6xlQoFGcOlZmpUCimRAmFQqGYEiUUCoViSpRQKBSKKZmx3cwV6SeR1AlEEoxF4oyFE4TjSRK6TlKXJHSJlBKryYTNomE3G/+67WZyXFZsZlO6zVdMI0oozmJCsQStg+O0eAO0D4Xo90fo9Yfp80cY8EcIRBN/9rUzbIZg5GVYKct2UpHjpDzHQXmOk7r8DPLdNoQQ03g3itOJEoqzhOFglAPdPvZ3+jjY46d5IEiPL3zaxgtGEwSjCTpHQuzt9P3J73NdVhYWe1hQ5GZxqYc1lTmUZTuUeMxQlFDMUfr9EV5rGeK1liF2dYzQNfKmKJg1QUKXJ3UdsyaQQPIkzz9ZhsdjbGsZYlvL0OSxAreNtVU5rKnKZlNdHnUFGUo4ZghKKOYI8aTOG8dHeO7oANtahmgZDAJgMQk0IdAETHzWJ0TCYhLYzSYiiSS6BJMmMAlBQteJJ+Vbzj0TDAaiPHmojycP9QFQmuXggvn5XDi/gI11uTit6u2aLoSUZ+6N8H5Ys2aN3L1bbQd5LyLxJK82D/H04X5eaBjAF4pjNWnYzBrRhE5c15n485o1gcdhIZHUiSZ0NCFISkksoU9ez20zIwRICWaT4UkkkpKErmMShuAkdEkydVFT6hp66phZM779dWl4IEJgjDMNYuOwmLhoYQHXLivmgvkF2C0qWDodCCH2SCmnrDinJHqWIaVkb6ePh/Z288SBXsYiCWxmDU0InFYTupQEoglMmqDAbSOW0Enoxod5ZDwGgNtuxm0zk0itXghA0wS+UGzSk5jAatbIdVnRhCCW1EkkDZExmwwBiCdTqx9mE1JKYkkdPSUM47HkpKicKuF4kicP9vHkwT5cVhOXLy7iw2vLWVedo6YnZwDlUcwSRsZj3Leriwd2d9E2NI7FJHDbLWhCEIjEiSZ03HZzallSktQlo6E4AEUeO1azllrW1BmPJgmmVjTMmqAix4ku5aRImE0CXUoSSUkwkvizVj/MmsBlM6fERMesCUzam+JyordzKtTku/jI2nI+uKqM3AzbqV/wLONkPQolFDOcI71+7trezqP7e4kmdEoy7cSSkmA0TiSu47SayHJY0DSBPxwnEElgt2iUZDnQUx6DNxAlmppi1OS5sJo14kkjDqFLSZ8/8pbpQV6GjWynZXIaEk/qCCGY+N5O6BKJRCAQwphm6Lo0zkn9n4RueB4T/8cfjjMeS06OMREfSUpDSKSEpJR/tnhYTRo3rCzh05trqC90/3kXOQtRQjGLkVLyavMQd77Uws7jI1jNGpkOC2ZN0D8WQQDVeS50CUPBKIFIApfVRHGWA4BoIjm5ylGd5yLDZiaW0IkldbpHQ5OeQ2mWgyynhaQujQ+1MGIS/WMR/OH4O9pmNWm4bCZMmpaKQRjHYwmd8VjyLTGPE7GZNTwOC6bUOBPiEIolCccNAbFbjClULKGfkmhctKCAT2+u4ZwaNS2ZCiUUsxApJS8cG+QnL7VwoMtHjstKlsOCLxxnZDxGjsuKx25Gl9A5EsKkCeYXujFpAl84RtdIGLMmWFDsRtdBl5Lu0TDBaAJNwKISDyZNY+Jv3usLMxSMTY5f5LFTV5BBdZ6L4iw7xZl2ijwOCjw2PHYLbrt5yiBiPKkTiiXxhWIMBqIMjkUZGIvQ6wvTNjROmzdI50hocgXGZTXhtluIJ3V0aUyZxiLGVMdq1jAJQTyp/1mrL+uqc/jqpfNYX5P7vv/v2YISilmElJKtTV6+93Qjx/rGyHFZybCZ8Yfj+MNxKnOduO1m+v0RhoIxSrMc5GVYCceTNA0Yy6ALitzYUh/iIz1+ErqkNMtBcaZ9Mg+ioX+MSNz4xq/Nd7GiPJuVFVksLvFQW5CBx245I/cbTSQ5PjTOwW4/B7t9HOjy09A/RjxprJzkZlhJJOWk6EyIxIRwTHggJ8vm+jz+5tJ5rKzIPh23M6tRQjFLONTt5ztbjrG9dZhCj41sp5UeX5hAJMH8QjdCGN5DKJZkQZGbDJuZrtEQA2NRSjLt1ORnEEvoHOj2EU3olGU7qMhxoglB92iI9uEQYAjD5vp8NtfnsaYqh0zHmRGFkyUST7KnY5RtqSSxQz1+pDRWaBwWEwld4g/HU9MkYwn2xNyQk+G65SV87coFlKSmaAolFDOePn+Yf9/SwKP7e/HYzRRnOhgKRhkej7Gw2IPLauJgjx9dl6yqzEZKycFuP9GEzqqKLKxmjR5fmK6RMG6bmdVV2QigxRukaySMSRNsqMnliiVFXLSgYNZ9OEbHY7zQMMjTh/t4pXmIWELHYTFht2iYNI2hYBRNGF5GLKGftGDYLRqfPb+OO86rwWFVuRhKKGYo8aTOXdvb+a/nmkjokkUlHgb8EXr9EeoKMnDbzTT0BYglDUEQQrC3YxQh4JyaXGxmjZ1tIwSiCZaXZ5HjtBCIJNjTOQrAubV5XLeihEsXFpLtsqb5bqeHYDTBiw2DPLy3m5ebvOgSsp0WEklJNBWkBSPQOvF4KkqzHPzbB5Zw4fyC02n6jEcJxQxkd/sI//TIYRr6AywocmMxaRzq8VOSaae2IIPDPX584TjnVBvBtzfaR7CYBJvq8gjHk+xoGwHg3Lo8XFYT+zp99I9FKM1ycPOaMm5aXUZZtjOdt3ja6fdHeGhvNw/s7qJ9OITLasJlMzMeTTAeS04u/Z7s2/qGFSV889rF5MwRUX2/KKGYQYRiCb73dCP/s72dQo+N8mwnh3v9AKyqyKbPH+H40DhLSj0UuO1sbx0iqUvOn2d8273UOIhZE1yysJBoQue1liHC8SSb6vL41KZqzp+Xj6adXcuAui55ucnLb7a18VrL8GRehi7lZL7GxLRkKnJcVv7lusVcu7zkdJs941BCMUPY0zHCV+8/QPtwiI21uXgDUZoHg6ytykYIwRvHRyjNcrC0NJM9naN4A1HOn5ePxaTxctMgAsEliwqQEl44NohEcv2KUj61qZqFxZ50396M4GjvGL95tY1H9vdgNmlkOSwEowlCseTkvpeT4cZVpXz7+iVk2M6enQ1KKNJMLKHzH8818qtX2ijy2KnKdbG7YwSP3cLqymx2d4wSiMS5dFEh/f4Iezt9LC7xsKjYw/PHBvCH41yxpAghBM8fHSCpS25eU8bnL6yb89OLP5c2b5Afv9DMowd6MaX2vkQTE5vgQJzEBrXKXCc/+shKVpRnnSGr04sSijTSPRriC/+7j/1dPjbV5TEYiNA0EGRjbS7heJJ9KVGoyHHy3NEBXDYzFy8s4EjPGI0DAdZWZVNXkMGWw/2MhePcuKqML11UT0WuEoiToXkgwA+ebeSZIwOTSWJDwSgWTUPTmMwleTdMmuAfrpjPpzfXzPnMTiUUaeK5owP87QMHSOqSzfV5bGseAgGXLirk5UYvwWiC65aXcKjHT0N/gMsWFWKzmHjiYC/FHjsXLyxk5/FhmgaCnFOTwzevWcyiEjXF+HPY3jLEvzx+lMaBwOR0YmIz3MnkYFyzrJjv3bRsTtfBUEJxhknqku893cAvX2mjviCD3AwrO9pGWFziochj54WGQeYXullU4uGJg71kOixcsaSIlxq89PrD3LCilFhS58mDfZTnOPjGVYu4fHHhnP9GO90kkjq/f6OT7z3TSCSeJNdlY2Q8RlJKnBbTlDtjFxS5+dXH18xZb04JxRlkLBLnS7/fx9ZGLxcvKKDHF6ahP8B1y0to6B+jaSDIdctL6PGF2dMxyiULC3FYTTx+oJe6ggwumJfPI/t7GQ3F+PTmGr5ySb0qzDLNDIxF+KdHDvPcUWM6ktQloVjypHIvMh0Wfn3rGtZV55wha88cSijOEO1D49x+927ah8a5Zlkx21qGiMR1blhZwuMH+pBS8qE15Tyyv5dgNM7Hz6nklaYhGgcCfHhNOf5wnKeP9LO4xMN3P7iMJaWZ6b6lOYuUkicP9fGtR4/gD8fJcVkZThXzsZq099xDYjVp/OeHl3PNsrm1hKqE4gzwxvER7rjHsPHC+QU8fqCX8hwnK8uzeHh/D/ML3SwtzeTBvd3U5mewuT6Pe3d24rGb+ej6Sv64t5s+f4SvXFzPZy+oxWxS/ZjOBN5AlK8+cIBXmrwUuG0EIgniSSNFfKqpyDeuWsjtm6vnzJRQCcVp5tkj/Xzx9/sozXawvCyLh/f1sLE2F4fFxAsNg1y2qJCkLnmhYZCrlxZjNgke3d/L5vo8avMzuGdHByVZdn70kZWsUrsazzi6LvnttuN875kGTJrAZjbhD8cxa2LKWhif3lzN169aOCfE4mSFQn2F/Rnct6uTv/rdHuoLM6jJc/Hwvh6uWFxEPKnzQsMgHz+nkj5/hBcbB/nLc6vp9Yd5dH8vt26oxKwJ/md7O9csK+bJL21WIpEmNE3w6fNqeOizG8l12QjFEpTnOEjo0uh6Zn73j8avXz3Otx47Mlkb9GxACcX75Bcvt/IPDx1iTVUOTouZ548NcvPqMlq8QfZ3+fjM+TU8c6SfNm+QL19cz+MHe2nsD/C3l81jW8sQrzYP8S/XLeaHH15xxuo/KN6dZWVZPPqFc1lZkU3XSJhCj41YUkcIIy7xbtz9egdff/jQWSMWSijeB3e+1MK/b2ngkoUFJHXJ7o4RPnNeDa82DzEwFuGLF9Xzvzs7MWmCL11cz69eacNq0vjKJfX84uU2/KE4996+nts2Vs0Jt3WukJdh497b1/PxcyoZGIuS5bAgMPqbvJdY/GFXF3//0MGzQiyUUJwkd77UwvefaeSShYWMjMfY3+XjsxfU8tDeHhK6zl+dX8udL7WQ77bx4bXlfO+ZRqrzXNyyrpzvPt1IZa6Tx7+4SZVlm6FYTBr/esMS/s81ixgejxktDewWdCnfs8jPg3u6+fYTR5mpsb7pQgnFSfCmSBQwGIhwsNvPZ86rSXkP8NF1Ffznc03ML3JzycJCfvh8Mxtrc9lYm8sPnm3i3Lo87vvMhllXPOZs5FObqvnxLSsZDcWQUmK3mBiPGsWL4c1iwifyP9vb+dELzWfY0jPLtAiFEOIKIUSjEKJFCPG1d/j93wghjgohDgohXhBCVE7HuGeCu7a38/1nGrlicREj4zGO9Y3xhYvq+N2ODpxWMx8/p5KfvtTC2qpsNtTk8qtX2rhmWTGlWQ5+/epxPriqjN/etuas2pE427lueQn/88l1k93OPA4L4XgSu0WbLMH3dn74fDP//drxM2/sGeKUhUIIYQLuBK4EFgG3CCEWve20fcAaKeUy4EHge6c67pngiYO9/PPjR7hwfj7jsQT7unx86aJ67nm9A5fNzK0bKvnP55pYU5XDkpJMfvlKGx9YWYrdYuIPu7r4zPk1/ODmZVhUfsSs49y6PO76y3WT8Yccl41EUpLjsr6rWHz7iaM8e6T/DFt6ZpiOd/A6oEVK2SaljAF/AK4/8QQp5UtSylDq6Q6gbBrGPa1sax7ir+/bz+qKbCwmjVebh/jChXX8/o1OAD6+oZLvPdPI6spslpZm8pttx7lptXFbD+7p5ksX1/O1KxaooOUsZk1VDnf95Tqi8STxpNFsKRhJ4LYZLRMsprf+baWEL/9hP4d7/Gmy+PQxHUJRCnSd8Lw7dezd+BSwZRrGPW0c7R3jM/fspjY/g6o8F88eHeCO82p48mAfgWiCz19Yxw+fa2ZZWSbn1uXx223H+dCaMhJJnYf39fA3l87jby6dp0RiDjAhFomkjtmk4bQZ3d9tZo2ELv9ELMLxJLfftZuBsUiaLD49nFGfWAjxMWAN8P13+f0dQojdQojdXq/3TJo2iTcQ5fa7duFxWLhsUSEP7unmo+sr2NMxSo8vzD9csYD/er6JilwnVy8t5ofPN3P10mJsZhOP7O/lq5fO40sX16fFdsXpYU1VDr/4+GrGwnFcVqN9gCYEGVaj0bPpbfOQ/rEId9y9m2ji/fUfmclMh1D0AOUnPC9LHXsLQohLgG8A10kpo+90ISnlr6SUa6SUa/Lz86fBtPdHNJHkr363h5FQjFs3VHHn1lYunJ+PPxxnT8cof3f5fH76YgtOq4mPn1PJ/3vqGJvr86jJd3HPjg5u31TNFy6qO+N2K04/m+vz+cHNy+nxhXHbjTaMZpPAYtLQxJ/GLA50+/m/Tx5Lj7GngekQil1AvRCiWghhBT4CPHbiCUKIlcAvMURicBrGnHaklHzj4cPs6Rjls+fX8bOXWqgvyKA2P4MnD/bx+QtreXBPN+PRBH9/+QL+fUsDS0szOa8+n5+82MKNq0rnTP6/4p25YWUp/3T1Qnp8YUqzHYyG4jitRnMiu8XE2//0d7/ewaP7/+Q7c1ZyykIhpUwAXwCeAY4B90spjwghvi2EuC512veBDOABIcR+IcRj73K5tPG7HR08uKeb2zdV89iBHmwWjZtWl/GbbcYSZ5t3nKaBAN+6bjE/eLaRbKeFWzdU8Z0tx7hoQQHf/eCys64S9tnI7ZtruGVdBceHxllc4sEXipOfYZss5Pt2/vGPh2gZDKTB0ulF7R4FDvf4ufFn29lYl0uGzcxTh/r41xuW8J2nGqgryOC8efn8+IVmvnrpPJ5vGKSpP8D3blrG1x8+RGmWg4c+uxGXypM4a4gldD766x0c7vVT5LHTNRqm0G2j1x/BbtGIxPW3lNpbWOzh0c+fi/U9NpqlC7V79CQZi8T5/P/uJTfDysrybJ442MfnLqjj7u0dWEyCm9eU8ZMXm7lhRQntwyEOdPn4l+sX85/PNWE1afzmtjVKJM4yrGaNn31sFVkOK75wnEyHBX84TkZq2dRq1jjx6/dY3xg/eqEpbfZOB2e1UEgp+ceHDtE9GuaO82r46UvNXLSggO7REE2DAf7xqoX84JlGFhR5WF2ZzUN7u/niRXU8ebCP7tEQv/j4alU6/yylwG3nFx9fTTCSIMtpmVwqlVJifwfP4edbW9nTMZoGS6eHs1oo/ri3hycP9fH5C42U7FyXjfPq83hkfy9fvLCOh/Z0E03ofOmiOv7vU8fYWGv0/ny5ycs3r1nE2qq5V0NRcfKsKM/iby+fT5t3nPlFbkZDccqynYxFEpMp++ZU3EqX8NX79xOKvXcFrZnKWSsUvb4w//z4EdZWZTM6HqPVO87fXDaP/3iuidWVRhevncdH+NqVC/jBs4247RZu3VDJf6XyJj52zqzZrqI4jdyxuYZNdXk0DQQoy3bQ4wtT5LETiSdxWU1ImMyzaB8OzdrNY2elUEgp+YeHDpJIGu357tnRwSc2VvHw3h50XXLrhsrJuERjf4C2oXG+fd1i/uXxo5RlO/jOB5eqZVAFYFTK+s8PLcdpNaPrEluqSbIQAotZQ0rJie+U3756nIb+sbTZ++dyVgrFvTs7ebV5iC9dXM9PX2yhriCDAo+N19uG+dqVC/jJiy0UeexcvriIe3d2cvumap4+0o83EOWnt6xSlakUb6HAY+db1y6i1x+hriCD4fEYRR47vlCcbKeVhC4nC+AkdMnX/zj7KmOddUIxOBbhu1saOLcul35/mIFAhL++ZB4/eaGFixYU0OeP0DIY5OtXL+TbTxylriCDJaWZPLq/ly9eVM/SMlVOX/GnXLe8hIsXFHCsb4z6ggy8gShl2Q6CqVoWQrw5Bdnb6eO+3V1TXHFmcdYJxb8+eYxoUudDa8q5e0cHt55TyUN7uxECPrSmnF++0sbNq8t4udHLYCDKN65eyLcfP8qSUg+fu7A23eYrZihCCP7tA0uwaMbUQ2LsAYkndVw2M9GEsfsUDMH4j2cbCUTiabb65DmrhOLVZi+PH+jljs01/OLlNgrcNhaXZPJiwyBfuaSen7zYTK7LymWLi3hgTzef3lzDH/f2MBaJ8x83r1B1JRTvSXGmg7+9fD7twyEWFHnoGA5Rlu1kKBilwG1kb07sNh0Kxvjly21ptvjkOWve+dFEkm8+eoSqXCcZdjPH+sb4u8uNFY0lpR5sZhNHesf4x6sW8P1nGijNcrCmMpvHD/Ty+QvrmF/kTvctKGYBf7G+gnmFGfT7I5RmOQjHk6nam8bv7WYTSV2iCfj1q230+cPpNfgkOWuE4u7tHRwfGuevL53Hz7e2srk+jzZvkMFAlK9eNp//er6Jc2pyGByL0jQQ5J+uXsh3thyjIsfJX52vphyKk8Ns0vjWtYvpH4uQl2HFG4hSkuVgKBil0G0jFE9OFr6JJnR+8MzsyNg8K4TCF4rxkxebOW9ePge6/AQicT55bhW/2XacD6ws5fmjAwQiCb5wYT0/eqGZSxYW0DkSotU7zreuXaQaBiveF+fW5XHZokKaB4PU5Lvo84fJy7ASS0rMmpjc82G3aDy8r5vjQ+Nptnhqzgqh+OmLLQSiCT66rpy7X2/nw2vLeXR/LwJj6/Dv3+jk4+dUsuVwH7GEUXp/QjAuXliYbvMVs5C/v2I+kXiSXJcVXyiOx25hKBilONNutAOwmdF1Iwj60xdb0m3ulMx5oegaCXH36x3ctKqMR/f3YjNrXLygkEf393LHeTX8706jmvY1y4r5w64u/mJ9BY8d6CWa0Pn6VQvTbb5illJX4Ob6FaUc6vGzvDwLb9BYLg1EEjgsJixmjVjSaDD0yP4eOoZntlcx54Xixy80IwRctbSYLYf7+dTmGu7Z0UG208KGmlyeOTLA7Zur+dUrbTgsJq5fWcr/7uzkw2vLqcnPSLf5ilnMly+uJ56abgQiCbKdVobHYxR4bPhCMXJdVpK6RErJnS/NbK9iTgtF10iIh/f1cMu6Cu7b1YXbbmZVRRYvN3m547xafra1lRyXlbVVOTx7dIDPnFfDXdvbMZsEX1Z1LxWnSFWeixtXGl7FklIPXaMhCj024gmjUK/FZHgVLpuZR/b1MhiYuQV557RQ/GxrK5oQnD8/n6eP9POX51bz36+1k+20sKDYzbaWIT53QS3//Vo7WU4Lm+rzeHR/L5/YWE2hx55u8xVzgM+cX0ssoZPtNGIVWQ4rvf4IFTlGfsVEu8JYUufeHZ1ptvbdmbNC0eML8+CeLj68tpz7J7yJyuxJb+Ke1zvIcVlZVZnN88cG+OTGau55vQOHxcQd59Wk23zFHKGuIIOLUqndC4s9BKMJcl1W4kmdhC5xWk0EIgk8djP37uyYsZW756xQ/OZVI+vtyqVFPH2kn1s3VPKHNzrJdFhYX5PDiw2D3Lahiv9v23EybGYuWVTAowd6+ci6cnJc1jRbr5hL3L6pmqFgjEyHmR5fmMpcJ50jIarzXIRiycnaFUPBGI8f6Euzte/MnBSKQCTOA7u7uWZZCc8eGcCsCS5aUMgzR/q5ZV0Fv0t5DhfMz+fJQ3187JxK7t/VhSbg05uVN6GYXjbU5rKw2MNQ0AhgRhM6AqPTmD8cJzfDSjCaINNh4Z4dHek29x2Zk0LxwO5ugtEEH1xVxv27u7h2eQlPHepDE4LLFxfyWMpzePxALyYh+OCqUu7b3cX1K0pVx3HFtCOE4KPrymkZDLK0LJOjfWPMLzKEYyK4OVFr80CXj+aBmVe1e84JRVKX3PV6O6srsznS6ycUS/LhNeXct6uLq5YW83KTl6SUfGRtBffv7uKKJUW80jxEJK7zyXOr0m2+Yo5y3YpSbGYNgdGj1GU1MTIeI9tppW/M2BcST+qYNcGDe7rTbe6fMOeEYmvjIB3DIW7bWMW9OztZV51D02CQYDTBbRsr+cMbXWyuz2dv5yhjkQS3bqji3p0drKzIYnGJqjWhOD1kOixcuaSI3R2jLCvLZGQ8RoHbBhjCYTNrxqqI08of9/WQSOpptvitzDmh+P0bXeS7bWQ7LXSOhLhlXTkP7O5iQZGboWCM/rEIf7G+gnte72BhsYd4UqfNO87HVQ1MxWnmQ2vKUyscFtqGxqnKddHmHae+IINoQsdu0TBpRv/bV5uH0m3uW5hTQuENRHmpcZAbV5by0J5u3HYzVbkuDnb7+fDacn7/RieFHhulWQ6O9o1xyzrjWJbTwlVLi9NtvmKOs74ml7wMK0KA1aSR0HViSaOgTY8vTHGmg0hcn2xCNZOYU0LxyL4ekrrk8iVFbDncz/UrSnjsQC9Wk8amujxeafJy0+oyHjvQi1kTXDi/gOeODnDd8hK1Q1Rx2jFpgssXF7GnY5T1NTn0+yNU5jpJpOpTCMAfjpPltPDcsQHiM2j6MWfIUwc5AAAgAElEQVSEQkrJ/bu7WFmRRUNfgGhC58ZVZTyyr4dLFxWyvXUYXcI1y0p4eF8PFy4o4I3jI0QTOtevKE23+YqzhKuWFhOKJbGaNHr9ESpzXTT2B5hf5DGK8Ka2oPtCcXa0DafZ2jeZM0JxtG+M5sEgN60u48lDvdTkuQhGEoyG4ly3ooRH9/ewsNiDNxDFG4hy48pSHtnfQ3mOg1UVWek2X3GWsL46Z7KzmBCg6zLVDV2jcyRETSoJy2k1seVwf7rNnWTOCMXTh/vRBKyryuH11mGuWlrMU4f6cFlN1Oa72Nvp47rlJWw5bBxbXp7Fay1DXL+8VPXoUJwxzKlp8NG+MVaWZ03u9zCl3oNCCEbGYxRl2tnaMMhMaSI+Z4TiqUN9rK/O5Y32EXQJly8u4pkj/Vy8sJCtjV4Arl5azPPHBrlgfgHbWobQJVyxpCjNlivONs6rz0+V83fS0B9gaWkm7cPjVOe53tLvo9cfodU7M+pUzAmhaB4I0Ood56qlRTx1qI+aPBf+cJzRUJyrlhbxwrFBY3l03Jh2XLa4kBePDVLksbO4xJNu8xVnGZvq8wAjORDAbjExFIyR7bTQPRqiMvfNxtfbmr1psfHtzAmh2HK4HyFgY10eO9tGuGxxES82DGI1a6yqyGZX+wgXLTBWOMyaYGNtHq82e7loYYGadijOOCVZDuoKMvCH4zitpsnphVnTGI8lyXJa6fWFKcm0z5h8ijkhFC81DrKsLIuWwSAJXXLB/HxeafayvjqHN9pHSOiSixcW8FLDIGurcmjoH2M8luSShQXpNl1xlrK2KptDPX5WV2bTORKiJNPOxHeWWRNE4jpuu4XX24ZnRJbmrBcKfyjOgS4f59fn8XKTlwybmZJMBy2DQc6rz2dro5dsp4WKHBcN/QE21efxWsswFpNgQ01eus1XnKWsLM/GH45T4LbTPBikvtBN29A48woziCUMYTBpglAsSUN/+jeJTYtQCCGuEEI0CiFahBBfe4ff24QQ96V+v1MIUTUd4wK81moEJTfPy+flRi8ba3PZ3mq4a+fPz2dH2zDn1OSyq30EMLb87jw+zLKyLBxWlWSlSA8rU0vyE9MOl82ENxAl12Wj621xin1dvrTYeCKnLBRCCBNwJ3AlsAi4RQix6G2nfQoYlVLWAf8FfPdUx53glSYvbpuZbKeFHl+Y8+bls61liEKPDZfNTPdomHXVOWxvHcJlNVFXkMGhbj/rqnOmywSF4n1Tm5+B22Ym9rZphUkT+EJxCtw2BgNR8jJs7OsYTZOVbzIdHsU6oEVK2SaljAF/AK5/2znXA3elHj8IXCymKYr4WusQG2pz2dtpqO45NTns7RhlXXUubxw3MtvWV+eyo22EtdU5HOzyk9Al65VQKNKIpgkWlnjo90coz3EQS+iYNIGW6nhu1rTJnqV7O+eGUJQCJ/Zw704de8dzpJQJwA/knurAg2MRukbCrK3KYU/7KFlOCw6rmV5/hNUVWexsG8FjN1Oa7aDVG2RVRTZ7OkYRAlZXZp/q8ArFKVFfkEHzYJAlJZk0DwapzXchpcRiEuipKYnNotE+HEp75/MZFcwUQtwhhNgthNjt9U69fjyhtKsqs9ndMcLqimz2nXBsf5ePlRXZHOsbQ0pYWpbJ4V4/1bku3HbLab0XhWIq6lNLpFlOK10jIeoKMmgZDFJX4J4Uigm3u3kwmD5DmR6h6AHKT3heljr2jucIIcxAJvAnO16klL+SUq6RUq7Jz8+fcuC9nT6sJo2ybAet3nFWVxkeg92iUZtvvOhLSj0c7DamJctKMznaO8YilWSlmAHUF7oBSCSNUng2s4k+f4TiTLvR3DjTPpnn05TmlY/pEIpdQL0QoloIYQU+Ajz2tnMeA25LPb4JeFFOQxL7no5RlpZlTi4frSjP4kjPGItLMmn1GjkVi0syOdjtpzTLgUkT9PjCqpKVYkZQned6y3N9MvFK0DUapjrfxXg0gdNqojHNdTRPWShSMYcvAM8Ax4D7pZRHhBDfFkJclzrtt0CuEKIF+BvgT5ZQ3y9JXXK0d4xlZZk09I0BsKjYQ0P/GAuK3BzuMY4tKcmksT/AwmIPR1PnqbRtxUyg0GPHpInJWigT2zxMmiCpS6wmjeND49QVZNA8kN6ph3k6LiKlfAp46m3HvnnC4whw83SMNUHnSIhwPMnCIg+vtw1TnGknFEsyFkmwoMhNQ/+YkXyVZadjOMRFCwsm28vXFaieoor0Y9IEhW4b4zGjAZApFZCYWPmIJyXRhFHxqms0lEZLZ1gw8/3Q2G94B/OL3BzrM7yIxtQUZH6Rh+NDxm68Hl+YWFKnNi+DjuEQVrNGkWoXqJghFGc56PNFKPDYiSV1PHbz5A5Sk/ZmBkGvL/yWnaVnmlkrFA39AYSA6nwXrd4g84s8NKXmcfML3bQPj1OVZxQvBajJd9ExPE5FjnNSsRWKdFOUaWcgECE/w4Y3ECXfbUNK0ITRIAhACMO7GAxE02bnrBWKpoEAVbkuRsdjxJOS2nwXXaOhVC6FiZ7RMNW5zsnpRnWei47hEFUnpMYqFOkm02FhLJwgP5WJmZthYzQUI99tw6wZH8+JsH93Gqcfs1Yo2oeMD33niPHiVeQ46RoJU5btoGs0hC6NtvP9YxFsZo0cl7F1t1R1AlPMINx2M2OROHkZNoYCUXJdVobHYxR67ITjSWxmbbKOZq8/kjY7Z61QdI2GKM95UyjKc5x0j4Yoz3bS5zNe0NIsBwNjEQpT87+xSIK8DFs6zVYo3oLHbiGW0LGYBaF4kmyXlaFglByXlZHxGFlOC1oql8IXiqXNzlkpFP5QnEAkQUVKKKwmjUKPne5Rw6PwBg2hyHfbUkJhYzhovMi5SigUMwiP3Vh4TCQlUhr9PoKRBC6bmfFogmyndfJcXyh9adyzUigmvIiybCe9vgjFWXbGwnGiCZ3iTAfeVNAn321jcCxKgcfOUNA4lpdhfdfrKhRnGpvZyKGwmCbiEUZVbofFNNnhPBhJ4LaZGVUexfujxxcGoCzbwXDQ2Io7EprwGKx4A1HsFo2M1Iub4zTcuInfKxQzhYk91HaL8VFMpiKXFpNgPCUU/nCcLJdFeRTvl+HxCe/AmFLkuqyT87cspzV1zIYQgmA0QYbdTCSeBMBpnZYcM4ViWpjIlRCp7V8TnoUmBOOxJBazRiypk2GzEIwm0mbnrBSKkVS8IdtlYXg8Rm6GlZFxQ21znFaC0QRuu5loIkk8KcmwmQmnhEK1DlTMJCYClUndKGAz8f6cqNBt0QSxhI7VrE2WyEsHs1MoQjEybGYsmmZMLVxWRscnPAoLoVgSh9VEMGIocIbNTDhmvMgOJRSKGYREpv5NPT9hvwcYHkY8qWM1CSUU75eRcUMcQvEkSV2S6bAQir0pCqFYApfVTChmeBEOq+kEj2JW3rJijjLxBTbhScgTdpCC0Vksnkx5FGmsxj0rPzVj4Tgex5txB7vFRDSltlazNulRTKizJsTkH0ClbytmEhNfcBPp2hOYUlmZQhi7Sic8i3QxK4UimtCxm01vCoXZNOmWWc0aCV2+5YUXvDkXTOfGGoXi7Uy8hyc0YMKzmFgNSaS8CV2S1mZVs1IoIvEkdouJSNx4dW0WjWhCRwjDZRMYcz3Jm6IwMedLKqFQzCDGY0nMmiCaMARjQgukZFIgrCaNeMKIU6SLWSoUOnaL9papR1zXsWgaQgiEeDMoBMaLPzHlSM6Q7tAKBcBwMEpuhtXIxrSaGI8mcFhMxJM6LqtpMj6R0PXJTWLpYFYKRTSRnMxoA2NqYRJiUgQ0IZDIyTXpiagxkNbIsULxdoaCMfIybAQiRr7PWDiBx2FmPJbAaTUb+0BMgnhSYjEroXhfSGl4CRPTCV0aopDUJVJKhBDoElyp5KpgNIknVXU7EElf0opC8Xa8qSY/3lSG8VgkbqziRZM4rSbGInHcdguReBKrSQnF+0LTjL4HE0JxYvAyoUucVhPh1MoHQCiawOMwhGIsnN7+CArFiQwFjWI1/f4IRR47o6EYmQ4LI+Mxsl1WRsfjZDstjIUNAUkXs1IoTMIoPvpmVpvEnFLbRCoTMxBNYDVrWE1GK/mJF3lMeRSKGcSVS4rZVJdH/1iEwkw7/X6jLII31SXMF4qR6bCm+n8ooXhfGFWK31x7jicNLwJgPGbM9YKpzkoZdjOBSHxy6pHOPf0Kxdv55rWLuGJJESPjMQrddvr8EUqyHJNl8XzhOBk201u+7NLBrBQKs0mQ0I3qxADByJtumT8cx20zT26gyUvtJs13G3Uo0ll3UKF4J3pTu6GdViNxMNNhbABz24zs4onkK+VRvE+cqWUk9wkByokYhC9kiMZoKI6UkgK3ncFAFIfVNFkOT6GYSbSmCkBPxNQmptSJiY1hZuN5rit9RZdmpVB47BYCESMGYTNrBFL79sEIVhZ47MQS+mT7+IlCNsWZhmunUMwkWr1Gc5/JZKtUoqD5hC3nACVZ6WszMSuFwm23TK5euO0WApE4WSmhGBmPUegxlHcgECE/JRS6LinOdCiPQjHjaB0MUuC2MeCPoIk3a1NM5P5M5GOmszD0LBUK82Q+RK7LylAwRlGmobb9Y5HJBj/9/ghlOU5iSZ2BQISKHCcdwyG130Mxo2gaDBptAweDVOa66PWFcdsn4hPGCp/FJNJaGHpWCkWW00IgmiCW0ClKLSk5rWaynBZ6fWEKTxCK6lyjEexx7zj1hRmE48nJUnoKRbqJJXSO9Y6xtDST5pRgtA8bXe5avUEqc530+SMUZdrTuvN5VgrFhMcwGIi8Je5QkppaFGfasZgE7cMhqvNTQjE8zrxCo+doy2B6G74qFBM09geIJXXmFbppHzLeo8f6xiYbE88rcNPqDVKTl95+ubNTKDLf9BgKPXaGx6PEEjolWQ76/BHMJo2qXEORiz12bGaN495x6vLdAJOtBxWKdLO/2wdMLPlLijIdDAVjhnAMG53MWwaD1Ke5sfasFIriTCOo0+c3PAopYWAsQnmOYzIGUZPvos0bRNMEtfkZNA4EyHRaKPTYaOhXQqGYGRzo8pHrsk6uzFlS0wuHxYQuwWkzcivqC5VQvG9ODFZW5aWmFkPj1Be4J2MQtflG9/J4UmdJqYcjvWNIKVlelsX+Ll86zVcoJtl5fJjVldns6/KlOttFjbIIqXDExNJoXYE7jVbOUqHwOMy47WY6R0LU5htK2+oNTsYgmgcDzCt0k9Alrd4gS0szGRmP0euPsLIim+ND45PFeBWKdNE1EqJrJMy5dXns7/SxojyLg90+avJcNA0EybCZJ8siKI/iz0AIYzrR6g2Sl2HFYzfT6g1SXzARgwiytCwTgIPdfpaUGo8PdftZWZEFoLwKRdp5rWUIMBpZ9fjCrKnK5o32EdZV57Cva5Tl5Zkc7PZRm++a3KuULmalUACTQiGEoLYgg9bBcTKdFgrcNhr7A1TnusiwmTnY7WNhsQezJjjY7WNZWSYmTbCnYzTdt6A4y3mtdZgCt43+MWPVbqKAzdLSLI71BVhRnsX+Lj/Ly7PSbOlsFooCFwNjUQKROPMK3DQOBJBSsqwskwPdPjRNsKTUw6FuP3aLiaVlmew8PoLTamZ5WSbbUmquUKSDRFLn1WYvm+ryeK1liJJM++SGRZfNRFKXFKZ65q6c7UIhhMgRQjwnhGhO/Zv9DuesEEK8LoQ4IoQ4KIT48KmMOcFEbKJlMMiSMiMG0T0aZmVFNm3ecfyhOMvLDGWOxJNsqMnlQJeP8WiCTfX5HOz24U9jL0fF2c0b7SP4QnEuWVTI9tZhzq3LY2fbMOU5Do4PjacCmkYgcy54FF8DXpBS1gMvpJ6/nRBwq5RyMXAF8EMhxCnf+aJiDwCHe8dYfkI8YkXqRd3f7WNddQ6xpM7ezlHOqckloUv2dIyyuT4PXcL2VuVVKNLDs0cGsJmNRtq+UJxzanJ5rWWIzfX5bGseYmlpJsf6xnBZTSxMvdfTyakKxfXAXanHdwE3vP0EKWWTlLI59bgXGATyT3FcyrId5LisHOr2Mb/IjcX0ZgxCCNjfaQiFSRO83mosQZk1wfbWYVaUZ5FhM/NKs/dUzVAo3jdSSp490s958/J5pcmL1aQZpRNiSdZX57Cvy8emujy2tw6zviZ3skh0OjlVCwqllH2px/1A4XudLIRYB1iB1lMcFyEES0szOdjtx2Y2VHd/lw+33cL8QjdvtA/jtltYWprJ9tZhXDYzqyuz2do4iMWkcd68PJ47Oqj6fCjOOHs7ffT6I1y2qJAth/vZVJ/HzuMj2C0aWqrMY21+BseHxtlYm5tuc4GTEAohxPNCiMPv8HP9iedJo2ffu37qhBDFwD3AJ6WU71gzXwhxhxBitxBit9c79bf9srJMmgYChGNJ1lYZShyJJ9lUl8eu9lEi8SQba43YRDCa4NJFhTT0B+gaCXHlkmKGglF2t49MOY5CMZ08tLcbu0WjLNtJjy/MFYuLeKFhgHNr83i9bdgoDp3qWbOpPi/N1hpMKRRSykuklEve4edRYCAlABNCMPhO1xBCeIAngW9IKXe8x1i/klKukVKuyc+fenayvCwLXRo5ERtrc4kldPZ1+thUn0csofPG8RE21+eT0CWvNnm5dJHh8Dx/bIALFxRgM2tsOdw/5TgKxXQRiSd54kAvVywuYmvTICZNUJxlp2skzIULCnj2yAAXzi9gR9sweRk25hemNyNzglOdejwG3JZ6fBvw6NtPEEJYgYeBu6WUD57ieG9hbXUOmoDX24bffNw6xPrqXKwmjW0tQ6ytyibLaeG5owNU5rqoL8jguaMDZNjMnD8vny2H+1R9CsUZ4/ljA4xFEtywspRH9vVw4fx8Xm0ewmISFLhtDAWjXLSggK2NXi5dVJDWfqMncqpC8e/ApUKIZuCS1HOEEGuEEL9JnfMh4DzgE0KI/amfFac4LgCZDgtLSjPZ0TqMx25haVkW21uHcVhNrK7M5uVGL2aTxkULCnihYZB4UufyxUXsaBvGG4hy7fISBsaibG8dng5zFIopuX93N0UeO7qUDIxFuXFVGY/t7+X8efnsPD6C1aRhs2gEowkuW1SUbnMnOSWhkFIOSykvllLWp6YoI6nju6WUt6ce/05KaZFSrjjhZ/90GA+woSaXfV2jhGNJNtXlsr/LyI+4ZFEhjQMB2ofGuWxRIf5wnF3tI9ywsgRdwmMHerl0USGZDgv37e6aLnMUinel1RvklSYvH11fwR/39pDttOC2m+kfi3Dt8hKePtzP5vo8XmsZxmU1sWGGBDJhFmdmTnBObS7xpGRX+wgXLywkoUu2Ng1y+WIjHvF0ahnKbtF44mAfdQVulpZm8si+HuwWEx9YWcozR/pVvw/Faefu7e1YTRpXLini2aMDXL+ilMcP9OK0mvDYLfT4wly1tJjnjg5wwfwC7BbT1Bc9Q8x+oajOxWbWeLFhkBVlWeS7bTx7ZICybCdLSzN5+nA/TquZKxYX8cSBXqKJJDesLOVQj5+WwQA3rykjltB5dH9vum9FMYcZi8R5cE831ywv5tmjA8QSOtcuL+HxA31cv6KEJw724baZybCbGQpGuXpZcbpNfguzXigcVhOb6/N47ugAQsAlCwvY2jhINJHkiiVF7O/y0esL84FVZYxFErx4bJDrlpdg1gR/eKOLxSWZLCvL5O7X21VQU3HauH9XF+OxJB87p5J7Xu9gc30eh3v8hONJrl9RypbDfVyzvJgth/rw2M1cvLAg3Sa/hVkvFACXLCykxxemoT/AZYuKGI8l2d4yzDUpVX54Xw/n1uZS4Lbxx3095LttXLGkiPt3dxGOJfnUpmpaveO83KQyNRXTTySe5JevtLGhJpfu0TD9YxE+sbGKe3d2sLwsk66REKFYkiuWFPP0kX6uWV6CzTxzph0wR4Ti4oWFCAHPHx1gY10umQ4LD+/roTLXxbqqHB7a041JE9ywspSXGgYZHItw64YqxiIJHt3fw1VLiyny2PnNtrZ034piDnLfri68gShfvLiO3247TlWuE4fVRNNAkI+ur+B3OzupyXcxOBYhEtf54KrSdJv8J8wJoch321hVkc2Th/qwmU1cs6yYZ4/2E4wmuGl1GW1D4+ztHOWj6ypI6JJ7d3aytiqbBUVu7n69A7Mm+MS5VbzWMsyRXn+6b0cxh4gmkvx8aytrq7JJJCUHunzcvrmGX73SRq7LSkWOiwNdPj6xsYr7dnVRletkVcWfbMJOO3NCKABuWFFCQ3+AY31j3LiqlEhcZ8uhPq5aVozDYuKB3d1U5bm4YH4+//tGJ/Gk5LaNVRztG+P11mFuWVuB22bmpy+2pPtWFHOI+3Z10T8W4YsX1fPjF5opzrSzrCyTrY1ePrGxij/s6sRtM7OgyMPujlE+dk7ljEmyOpE5IxRXLzMClI/s72FVRTaVuU4e3tdDhs3MNcuKeexAL/5wnNs2VuENRHn6SD8fWFlKgdvGnVtbyHRa+OSmarYc7ldehWJaCETi/Oj5ZtZX52DWBLs7RvnsBbX892vtOK0mLl9SxJMH+7h5TTkP7enGYTFx8+rydJv9jswZochxWblgfj6P7utFSrhpVRnbW4dp8wa5bWMVoViS+3d1cX59PlW5Tn677Tg2s8anN9fwWssw+zpH+dSmatx2Mz98vjndt6OYA/zi5VaGx2N84+qF/PD5Zgo9NjbW5vHYgV5uWWckXSWl5PoVJTyyv4cbVpaS6Uxvbcx3Y84IBcANK0vpH4vwWusQH15XjsUk+N2OTpaUZrKuKoe7Xm9HArdvruFAl4/trcN8dH0FWU4Ld77UQqbDwqc31/Dc0QEOqOK7ilOg1xfmN68e5/oVxjaBN9pH+MJF9fx8aysmTXDT6jLufr2da5eV8FrrENGEzq0bKtNt9rsyp4TikoWF5Lis3LujkwK3nSuWFPPAni5CsQSfPLeK7tEwzx8b4KbVZRR6bPzkxWZcNjN/eW41zx8bZH+Xj0+eW0VehpV/e/Ioxs55heL9850tDUjgry+Zx3eeOkZtvos1ldk8vK+b2zZU8tShPkKxJLdvrua3rx5nc33ejKhk9W7MKaGwW0x8aE05zx0boM8f5tYNlQQiCR7ZZ+zrKMt28POtrdjMGnecV8uOthF2tY/wl5uqycuw8p2njpFhM/PVy+azq32UJw/1TT2oQvE2Xmny8viBXj53QS2vNHtpGxrn61ct5McvNOO0mvmL9ZX8z2vtXLmkiH2dPobHY3zhwrp0m/2ezCmhAPiL9RXoUvL7N7pYU5nN4hIPv9nWhhCCz11Qx/4uH682D3HLunJyXVZ+9HwzGTYzX764np3HR3ixYZAPrSlnYbGH7zzVQCRVQEShOBki8ST/9MhhavJcfHRdBT98vpkNNbnkZdjYcrif2zdX88e93QSiCf7q/Fp++XIrayqzWVedk27T35M5JxTlOU4unF/A71NLoJ+9oJY27zjPHunng6tLKcm086MXmnFYTHzuwjq2tQzxSpOXj6yroDrPxXefbkBKyTevWUSPL8zPt55y1T7FWcSdL7XQORLi3z6whP94tgl/OM7/uWYR//z4EfIybFy7vIRfv3qcq5cVc6xvjF5/hM9fVDcjl0RPZM4JBcAnUkugj+zr4colxVTnubhzawtWk8ZnL6hlT8co21uH+dg5FZTnOPjOlgY0IfiHK+bTNBDknh0dbKjN5foVJfxsa4vqfq44KQ51+/n51lZuXFmKWdO4b3cXt2+q5mjfGPs6fXztygX86uU2ErrOly+u54fPN7O8PIsL5p1yrenTzpwUis31eSwu8fCLVwxv4LPn13K4Z4xXmoe4eU05RR4733+mEatJ4+8uX8CxvjEe2dfD5YuLOG9ePv/xbBMDYxG+ec0iMmxm/v7Bg6oIr+I9icSTfOW+feRl2PjHqxbyjYcPUZrl4FObqvn3LQ2srMhicYmHB/Z0ceuGKl5qGKR/LMI/XrlgxnsTMEeFQggxOeV47mg/N6w0phz/9VwTNrPGVy+bx/4uH08c7OOapcUsL8vk+880Mh5L8u3rFhNL6vzbk8fIzbDxzWsXsb/Lx13b29N9W4oZzHefbqDVO873b17GvTs7aB4M8u3rF/OLl9sYHo/yL9ct5v+lguW3bqjkzpdauGB+PufUzJziNO/FnBQKgCuXFFOV6+RnW1uxmARfucQQh2eO9HPjqjIWFnv47tMNxHWdb123mIFAhB8+10RVnovPnl/L4wd6ebnJyw0rSrlgfj7ffbqBZjUFUbwDrzZ7+e/X2vnExircdsv/3955hzdZtX/886S76d57L9rS0sEGFWSqyBJUkCE4EReigry++LpeBDeKiIrsLaKCbFR2oS3dlO69VzrTNMn5/ZHKq/zUOoAy8rmuXH2SnCbf+0nyfc65z2L5kRzG9XLDVm7MlyfzmdrXi8LaVo5l1zBveBCb4opoalfzwsiQ7pb+p7lhjcJApuvlSClRcDCjkgnR7gQ6WbB0/wWEECy6owcl9W2sPVlAtJct9/X24suTBZwvb+Tx2/zxd5Sz4KsUGpVqlt4TgYWJIU9uPqfvBdHzK8oVbTy9JYlAJwueuj2QeVuTcLY0YdGdobywIwU3azOeGBLAa7sz6OluzYAAB1afyGditAehbtfuuIlLuWGNAmBCtDt+jnKW7b+AJEk8PzKYvOoWtieUMCjQgSHBjnx4OIcKhZIXRwVjY2bEoq9TMTaQ8e7kXlQ1tfOf79JxsjRl2aQIMiuaWLrvQneHpecaQaXW8sTGRNo7NKycFsP7h7LIq2nh7cmRrDmZT05VM29O6MmnP+VR3dzOG+PDefW7DEyNDFgw+vqpTcANbhSGBjLmjwgmu6qZr8+VMjzUmVhvW945kEWjsoNX7g6jQ6Pl1d3p2Jgb89IdPUgsamDNyQIiPW144jZ/diaWsj+9gqEhzswc4MPqE/kczKjs7vPKS4UAACAASURBVND0XAP8d+95EosaWHpPJFkVTaw7VcjsQb6YGhmw8qc8JsV4YGFiyLpTBUzr501xXRvHc2qYPyIYBwuT7pb/l7ihjQJgdLgLPd2tee9gFu1qLYvHhFHb0s67B7Lwtpfz1O2BfJ9awZFMXfNkaIgTb+3LJKeqmblDAwlzs2LBVymUNbSxYHQIPd2tmbc1idzq5u4OTU83si2+mC9PFPDgQB96uFry/I4UennaMOc2f57afA5Xa1OeHxnM/O3JuFqbMec2XfMj1NWKqX29ulv+X+aGNwpJklgwOoTShja+OJ5PTw9rpvb1Yt2pAtLLFDw82I9AJwte3pVOq0rDkgk9MTM24LltScgk+PD+KFRqLU9uPoeBTGLltBiMDWU8vC6eRmVHd4enpxs4mVvDSztTGRzowLPDg5izMREjA4mPp0az+Nt0yhVKPrw/io9/yCG/RtcT8sHhbCqblLw2LhzDa2DT4b/K9af4bzAwwIFRYS4sP5JNaUMbz48IwcbcmJd3pWEok3hzQk/KFG288f15nKxMeX1cOMklCpYfycHf0YIlEyNIKKxn6b5M3G3M+GhKNIW1rczbmqQfX3GTkVvdzOMbEvF1kPPRlGj+vSuNC5VNfHBfFMezq9mdUs684UG0qTSsPVXIgwN90GgFm88U8dAgX2K8r73Vq/4MN4VRAPzrrh4IAW/uOY+1uRELR4eQWNTAxrhCevvY8fBgPzbFFXEks5K7ItyYEOXOh0eyOZFTw5hIN6b39+azY/nsTS2nv789/74rlEPnq3j1u3T9LNObhKpGJQ9+eRZDmcTqmb1Ze7KAXUllzB8RjK25Mf/+Jp0B/vZMivHg2a1J+DnKmXNbAC/uSMHPUc5zI4K7O4S/zU1jFB625jwxJIA9qeWcyKnhnhgPBgc68Ob3mRTVtvLciCBCXCx5YUcqdS0qXhsXjr+jBU9vOUdlo5JFd/agl6cN87Ylk1aqYMYAHx4a5MvaU4V8elS/KO+NTn2Lige+iKOmuZ0vZvYmsaiedw9mMSHancmxnjyyPh4HCxPev68X87Ylo2jr4OMp0by9/wIVjUqW3RN5TW3o81e5aYwC4JFb/PCxN2fhzlTaOjS8NTECQ5nE/B3JGMl0XaKKNhUv7EjB3NiAT6ZG09Ku4clN5zCQJFZNj8FObszstWepUCh56Y4e3BXhypK9mXx9rqS7w9NzhWhSdjDjyzMU1Lby+YxY1Botz29PoY+vHa+NDeeJjYnUt6r4dFoMm+OKOZ5Tw6tjw8iuamZrfDGP3up/3TY5fuamMgpTIwOWTIygqK6VZfsv4GZjxstjQjmTX8eakwWEulnx4qgQDp2v5LNjeQQ6W/LmhHDOFNTx+p7zOFma8vmMWJqVamavPYtSreGdyZH097Nn/vYU9qXp16+40WhpVzN7bTwZZY2smBKNjZkxs9acxd3WjE8fiOH1Pec5U1DHWxMjaGjt4P3DWUyIcqevrz0v7Uwl2suGecODujuMf8xNZRQA/fzsmd7fmzUnC4gvqGNSjAe3hzixZG8maaUKZg/yZXS4C2/tu8CZ/DrGR3kwe5Ava04WsO5UAT1crfhoSjTnyxt5dH0CAJ/NiCXSw5q5m86xP72iewPUc9lQtHUw7Ys4EgrreffeXvg5ypm+Og65iSHrZ/dhY1whm88UMec2f3q6W/PEpkQCnSxYPCaMJzefu9hrZnQd9nJcyvUfwd/gxVEhuNuY8fyOFNo6NCybFImd3Ji5mxJpbtcN2fayM2fupkSqmnRNjNtDnPjPdxn8lFXNkBAnlkyM4Fh2Dc9sScLUUMbaWX0Id7dm7qZE/YCsG4C6FhVTPjtNaqmCj6dEEeNtywOfxyEEbHioL6dya3n7gK728Mgtfjy0Nh4DmcQXM3qz7EAmqaUKlk2KxMPWvLtDuSzclEYhNzFk2T2RFNS28J9vM7CTG/Ph/VEU1bXy0tdpWJgYsmJqNE1KNY+sS6BDo+WD+6MIdLJg7sZEzpc3MjnWk5fvCmVvWgULd6YiNzZk3ew+hLpZ8/iGBL5JKu3uMPX8TcoVbdy36hQ5Vc2smhZLmJs19606RVO7mrWz+lBc18rCnakMCnDgjfE9eWJTIsX1rXw6LYaj2dVsOF3Eo7f4MTLMpbtDuWzclEYB0N/fniduC2BrfDHfJJXSx9eOecOD+C65jPWnC+nhasX79/UiuaSB57YnY25kwBczeyM3MWTaF2fIr2lh9iBfnro9kO0JJSzalYaFsa5KGutjy9Nbklh9PL+7w9TzF0kvUzDu4xOUNSj58sHeeNubM/nTUzS2qdkwuy+Ktg4eXZ9AsIslKx6IZvG3aZzIqeXN8T0RAhZ/k86tQY68MOr6msvRFTetUQA8MyyQGG9bFn2dRmFtC3NuC2BoZxPjZE4NI8NcWDAqhD0p5bx3KAt3GzM2PNQHrRA88Hkc5Yo2nh0WyJzb/Nl8poj5nYay5sE+jAxz5tXdGSzdl6nfJf064YcLVUxeeQqZJLH9sf7Yy02Y/Olp2tVaNj/cj7YODbPXnsXHXs762X35+IcctsWX8NTtgQwIcODxDQl42Znz4f1RGMiu/cVo/go3tVEYGsj44L5eyCR4YlMiKo2WD+7rhZ+DnDmbEimsbeGRW/y4r7cny4/ksCmuiAAnS9bN6kNjWwcPfB5HTbOKF0aF8NzwIHaeK+XpLUkYyCRWTI3h/j6erPgxl7mbE2lVqbs7XD2/gxCCNSfyeWhtPD4OcnY9MZAmpZp7V51CJsHWR/rRolIza81ZPGzN2fhwX7bFF/PpT3lM6+fNrIE+zFh9BpVGy6rpsVibXZub+PwT/pFRSJJkJ0nSQUmSsjv//m5nsSRJVpIklUiS9NE/ec/LjYetOe9O7kV6WSMLvkrBwsSQz2fEAjB7bTxN7WpeGxfOkGBHFu1K5dvkMsLdrfliZm/KGpTct+oUFQolT94eyKI7erAntZxZa87SqlLz5vievHRHCHvTKrjnk1OUNrR1c7R6LqVVpeaZrUm88l0GQ4Kd2PZof87k1/HA53HYmRuz/bH+lCmUTP/iDC5Wpmx6qC8HMypZsjeTMZFuLLwjhNlr4ymqbeWz6bEEOFl0d0hXhH9ao1gAHBZCBAKHO+//Hq8BR//h+10RhoU6M29YELuSyvj8WD7e9nJWTImmoKaFR9cloNEKVkyNobePHfO2JnEks5I+vnasndWHCoWSe1edoqS+lYdv8WPpxAhO5dYyaeUpKhqVPHKLP6tn9Ka4rpW7lx/nZG5Nd4erp5O86mbGf3yS75LLeH5kMKumxbD+dCFPbj5HpKc1Xz0+gNRSBQ+tPYuPg5ytj/bn4PlKFu5M5bZgR5bdE8HTW5JILKrn/ft6XTfL2v0d/qlRjAXWdh6vBcb9ViFJkmIAZ+DAP3y/K8bcoQGMDnfhv3vP81NWNQMCHHh7UiSn8mp5dmsSxoYyvpgRSw9XKx7fkMjx7Br6+Nqx4aG+1LeouPfT0xTUtDC5tyerZ/ampL6NcR+fIL1MwZAQJ75+YiDW5kZM/TyOdw9m6SeTdSNCCHYklHD3RyeoalKydlYfZg305bntySzZm8ldEa6sn92X79PKeXLzOXp52rDlkX7sTStn0ddpDA1xYuUDMfxrVxoHMypZfFcod/R07e6wrij/1CichRA/D0esQGcGv0KSJBnwDjD/H77XFUWSJN6eFEmQsyVzNiSQVqpgXJT7xS7Qf3+j6zZdO6sPvg5yZq09yw+ZVUR52bKpM9E14ZOTJBTWc0uQIzse749Mkpj4yUl2nSslwMmC7+YOYkKUBx8ezub+z05ToVB2d9g3HfUtKuZsTGT+9mTC3KzY/dRgvOzMGb/iBLuSSpk3PIj37+3F+4eyWfR1GkOCnVg3qy9fJZTw72/SGdbDmRVTo3nl23R2JJTw9O2BzBzo291hXXmEEH94Aw4Bab9xGws0XFK2/jf+fy7wQufxTOCjP3ivR4B4IN7Ly0t0B+UNbaL/m4dE7OsHRVFtixBCiCV7zwvvF3eLN/dkCK1WK+qa28VdHx4TAS/tEXtTy4QQQuRVN4tblx4RQYu+v/hYZWObmLTypPB+cbdY/E2aaO/QCCGE+CqhWPR4ea/ouXif2B5fLLRabbfEerNx+HyF6P36QRHw0h7xyY85Qq3RikMZFaLn4n0i4pX94khmpWhWdoiH154V3i/uFgt3poj2Ds3Fz//RdfFC2aEWC3emCO8Xd4u392de958dEC+68AAhRNdG8Yf/DBcA185jV+DCb5TZCBQBBUAN0Ags6eq1Y2JirugJ+iOyKhpFz8X7xJBlP4ja5nah1WrFy7tShfeLu8UbnWbR0KoS4z8+LvwW7hFfJRQLIYSoaVKKcR8fFz4LdotVP+UKrVYrVGqNeO27dOH94m4xYcUJUVynM5+86mYxccUJ4f3ibjFzdZwoa2jttnhvdCoVbWLOhgTh/eJuMfzdH0VqSYNoU6kvfqaj3z8qCmtaREl9qxj1/lHhu2C3+PJ4nlCpNeLZred+ZRo/m8SSveeve5MQ4s8bhST+wVoKkiQtA2qFEEskSVoA2AkhXviD8jOBWCHE3K5eOzY2VsTHx/9tbf+UswV1TP08jh4ulqx/qC+WJoa88m06a08V8tAgXxbd2YNWlYaH18VzMreW54YHMXdoAO1qLfO2JfF9agV3Rbjy1sQI5CaGfJdcxsKdqUjAa+PCGRfljkYrWHuygKX7MzGSyZg3Iohp/byvyxWQrkU0WsGmM0Us3ZtJu0bLU0MDeOQWf7KrmnhmSxLZVc3MHuTL8yODOZVXy7ytSag1guWdQ7af2HSOo1nVzBsexGO3+jNvWxK7U8p5/DZ/XhgZfF1s3NMVkiQlCCFiuyz3D43CHtgGeAGFwGQhRJ0kSbHAY0KIhy4pP5PrxCgADmVU8vjGBMLdrVk3qw8WJob857sM1pwsYHp/bxaPCUOjFSz4KoWd50qZHOvBG+N7YiiTWPlTHsv2ZxLgZMGn02LxdZBTXNfKvG1JnC2o5+5IN14bG461uRGFtS0s+jqN4zk1BDtbsvjuUAb4O3Rr7Nc7x7NreH1PBpkVTfT3s+eN8eF42Jqz8qdcPjqSg425EW9PimSAvz3vHsxixY+5hLhYsmJqNJIk8fC6ePJrWnhjXDh393LjsQ2JHM2qZuHoEB691b+7w7tsXBWjuJJcC0YBsC+tgrmbEon0tGHtrD7IjQ34795MVh3N486errwzORITQxnvHcziwyM5DApwYPn9UdjKjTmWXc1Tm8+h1gheHx/O2F66WsTKn3J572AWtnJjXhkTxh09dXMC9qdX8vqeDErq27ijpwvzRwTj53hj9stfKbIrm3jz+/P8cKEaD1szXhgVwpgIVxKLGli4M4WsymbuinDl1bHhtHVoeHZrEmfy67i/jyeLx4RxOq+WpzrXR/14ajTBzpY8tC6e5OIGlkyIYHJvz+4O8bKiN4rLyPepum6yaC8bvpjZGytTIz47mscb35+nn58dq6bHYmVqxPb4YhZ9nYaTlQkrH4gh3N2akvpWntp8jsSiBsZHufOfsWFYmRqRVqpg4c5UUksV3B7ixKvjwnG3MUPZoWHlT7msOppHu1rLxGh3nro98IaZhXilyKlqYvmRHL5LLkNuYsjcIQHMGOBDm0rDuwez2BBXiKuVKa+PD2dIsBM7Ekp49bsMNELw+rhwxvVyZ+XRXN7ef4EgZ0s+mx5Lq0o3ZLu6qZ0P7otiVPiNM8nrZ/RGcZnZk1LO01vOEeRsydpZfXC0NGHXuVKe35GMv6MFn02PxdPOnKTiBh7fkEBdi4o3xvfknhgP1BotH/2Qw/IjObhYmfLu5Ej6+tmj1mhZc7KAdw5kIUnw2K3+PDzYDzNjA2qa21nxQy4b4gpBwOTeHjw82A9ve3l3n4prigsVTSw/ks2e1HJMDQ2Y3t+bR2/1x9LUkA2nC3n/UDZNyg6m9/dh/shgWlVqXtqZxqHzukFz70zSLVH33PZkjmZVc2eEK8vuieB0Xi1PbjqH3MSQz6bHEulp092hXhH0RnEF+PFCFY9vSMTZyoT1s/viaWfO8ewa5mxMwNBAxoqp0fTzs6emuZ0nN53jVF4tk2I8WHx3GBYmhiQU1vPs1iSK6lqZ0teLF0eFYG1mREl9K6/vPs++9ApcrEx5YVQw43q5I5NJlDW0sfxIDl8llKDWahkV7sLDg/2I8rq+l1b7J2i1gh8uVLHmZAHHsmuQGxswvXMNUzu5MYfOV/HfvefJq25hUIAD/7qrB0FOlmw+W8TSfRdo69DwwshgZg305WRuLc9sTaJJ2cG/x4Ryf28vPjuWx1v7MunhasXnM2JxtTbr7pCvGHqjuEIkFNYza81ZTAxlrJ7Zm3B3a/Kqm3l4XTyFta0svjuMaf28UWu0vH8om49/zMHLzpz37+1FlJctrSo17x7IYvWJfBwsTPjP3WGMCndBkiTO5Nfx+p4MUkoUhLlZ8cywIIb1cEKSJKoalaw5WcCG04U0KtVEedlwf28v7oxwRW5i2N2n5apQ29zO1+dKWX+6kMLaVlysTJnW35spfbywMTfihwtVvH8om5QSBX4Ochbd2YOhIU6klzXyr11pJBU30NfXjjfGh+NuY87bBy6w+kQ+/o4WfDQlClcrM57bnsyh85WMDnfhncmRmBvf2OdWbxRXkKzKJmasPkNDawfv3RvJqHBXGpUdPLMliSOZVUyIdue1seHITQyJy6tl3rZkKhqVPDU0kDlD/DEykJFaomDBzhTSyxoZHOjAv+4MJdjFEq1W8E1yKe8dzKaorpVQVyueuj2QEaHOyGQSLe1qtsUXs+F0IbnVLViYGDIm0o17YtyJ8rRFdoNNb1aptRzJrGRHQik/XqhCrRVEe9nw4EBfRoW7YCBJ/HChig+P5JBc3ICHrRlPDg1gQrQHirYOPjiUzca4QmzNjVl0Zw/GR7lztqCeF3YkU1DbygP9vFh0RyjZVU3M2Zh4cdHkBwf63BDdn12hN4orTFWTkkfWJZBU3MD8EUE8MSQArYAPD2ez/Eg2PvZylk+JIszNGkVbBy/vSuPb5DJCXCxZMjGCXp42qDVa1p4q5INDWTS3q7m/jxfzhgdhb2FCh0bLN0llfHQkm4LaVoKdLXlwoA9je7ljZmyAEIKEwno2nylmT2oZyg4trtamjA535c4IV6I8ba5b02hVqTmaVc2B9EoOZ1ahaOvA0dKECVHuTIzxIMjZklaVmq8SS/nyRD551S242+gMYmKMBx0aLZ8fy+fTn3JRqrVM6ePF/BHByGTwzoEs1p4qwMPWjLcmRtDX157Vx/NZtv8CjpYmLJ8SRfRN1KzTG8VVQNmh4cWvUvgmqYwxkW4smdATuYkhp3JreWbrOepbO1h0Rw+m9fNGJpM4kF7Bv79Jp7JJycwBPswfEYzcxJD6FhUfHM5m/elCzI0MeHCQL7MH+mJtboRao+Xb5DJWHc0js6IJG3Mj7uvtxbT+3rjb6NrOjcoODmVU8n1qOUezalBptDhamjA40IFbAh0ZFOhwTW+KK4Qgu6qZEzk1nMip4Vh2De1qLdZmRtzew4kxkW4MDnDA0EBGbnUz2+KL2XKmGEVbBxEe1swe5MsdPV3RaAXb44tZfiSHqqZ2RoW58PyoYPwc5OxKKuXN7zOpbmpn5gAfnh8ZTF2Liue2JXOmoI4Roc68NTECW7lxd5+Oq4reKK4SQghW/JjLOwcu4OsgZ8XUGIJdLKltbmf+9mR+uFDNAH973poYgaedOU3KDpbuu8CGuEKcLE14YWQI46N0icucqmbe3n+BfekVWJoaMmugL7MG+WJtZoQQgrj8OtacKOBARgUCGOjvwMQYd0aGuVxsSzcqOzh8vpIjmdUcz66mvlW3P2qIiyXR3rZEedoQ7W2Lr72822ocbSoN6WUKkksUJBU3cCq3lprmdgA87cwYGuzEyDAXevvaYWQgo7ldzZ6UMrbFl5BQWI+BTGJ4D2dmD/Yl1tuWVpWGTXFFrDqWR3VTO7HetiwYHUKsjx3pZQpe+TadswX1RHpY85+x4UR6WLPlbDGv785AJkm8cncYE6Ldb4qmxqXojeIqczK3hqe36LLnr44NZ3KsJ0IItpwt5o0959EKwYLRITzQV1e7SCis59XdGSQXNxDhYc3Ld4XS28cOgIyyRj48nK0zDBND7u3tyYwBPnja6cZSFNe1sj2hhJ2JJZTUtyE3NmBUuCujwl0YHOhwcUcqrVaQXtbI0exqTufVklTcQJNSt9KWpakhQc6WBDpZEOBkQaCzJR62Zrham16WBJ4QAkVbByX1beRWN5NX3UJeTQvZlU1kVzVfnGbvYmVKXz87Bvo70N/f/mKMjcoOjpyv4vvUcn7KqqZdrcXfUc7kWE/GR7vjZGlKZaOSjXFFrD9VQH1rBwP87Zk7NID+fvaU1Lfx3sEsvk4qxdbcmBdHBTMpxpO8mmYWfZ1GXH4dA/ztWTYp8mLN7GZEbxTdQFWTkqc3J3Eqr5a7Ilx5bWw4tnJjShvaWLgzlaNZ1cR62/Lq2HBC3awuJi7f2qvbdm5EqDPPDAsi1M0K0C30+smPuexNq0AIwcgwFx4c6EtvH1skSUKrFZwpqOOrhBL2pVfQpFRjZmTALUEODA91YVCAAy7Wphf1abWC3OpmzhU1kFLaQFZlM9mVTRdrHT9jbWaEq7UpdnJjLE0NsTQ1wtLUEHNjAyQkJEk3LR8haOvQ0KLS0NqupkWloba5naom3U2l1l58TUkCD1sz/B0t6OluTYSHDZEe1jhZ6fQJodN2NKuGo9nVnMypRaXR4mxlwuhwV8ZEuhHtpRvLkFBYz5qTBexLq0AjBEODnZgzJIAYb1uqm9r5+IccNsYVIpMkZg7wYc5tARgbyvjoh2xWHc3D3NiQhaNDmBzred3mcS4XeqPoJjRawSc/5vDB4WxszI15a2JPhoY4I4Rge0IJS/Zm0tCqYlo/b+YND8ba3Ig2lYbPj+Wx6lgeTUo1o8JceHpYID1cdYZR1tDG+tOFbIorQtHWgb+jnHtiPJkQ7Y5z5w9NpdYSl1/LgfRKDmRUUNmoq8r7Ocjp729Pf397orxscbM2/X9V7NrmdnKqmilTtFGuUFLeoKRcoaShVUWTUk2TsoMmpZrWDo1uNiHw89fGzMgAc2MDzE0MkBsbYic3xsnSBGcrUxwtTXC3McPP0QJve/Nf7b2p0Qqyq5o4V9RAQmE9x7NrqGjUrc/hY2/OsB7OjO75v6RshULJN0ml7Ews5UJlE1amhkyO9eSBft74OMgpbWjj82N5bDlTjEqjZXKsJ0/fHoizlQnfpZSzdF8mJfVt3BPjwcLRIdhfwzmbq4neKLqZ9DIFz21LJrOiicmxHiy6IxRrcyMUrR28e/AC608XYmNuzLPDArm3txfGhjIUbR2sPp7P6uP5NLWrGdbDiVmDfOnvZ48kSbSq1OxOLmd7QjFnC+qRSXBLkCN3RbgxrIcTNua6RJxWK8gob+R0Xi0nc2s5k19Hc7uuyWEnN+68olsT5GyJr4McP0f5X25u/Py9+TPtekVrB9lVuiZHdmUz58sbSSlpoEWlAcDG3IgB/vYMDnRkUIDDxeZHfYuKw5lVfJNUyvGcGoSAaC8b7onxZFyUG+bGhmRVNrHyp1y+TSoD4O5ebswdEoCfowVxebW8+f15kksU9HC1YvGY0Bt6ubq/g94orgHa1Ro+OJTNp0fzsDEz4qU7elxMmmWUNfLKd+mcya/Dy86c50YEMSbCDZlMQtHaweoT+aw/XUhdi4pQVytmDfJlTKQrJoa6q3J+TQs7Eor5OrGUMoUSA5lEfz97Roa7MCTY8VdzQ9QaLelluh9naqmClBLFr/IEoMsVeNmZ42hlgqOFCU5WJjhYmGBhomtymBvr/sokCa0QCAFaIVBptBdrHI1tHTQq1VQolJQr2ihrUFLRqKSuRXXxfcyMDAh0tqCXpw29PG2I8rLFx978ouEU1rZwMKOSgxmVnC2oQyt0Cc7xUR6Mj3LH10FOh0bLgfRKNsYVcjK3FjMjA+7r48lDg/1wtzEjrVTB+4eyOXS+EhcrU+aPDGZ8lPsNt4T+5UBvFNcQGWWN/GtXKolFDfTxteP1ceEEOVsihODHrGqW7rvA+fJGQlwseWZYICNCXZDJJJQdGnadK2X1iXyyKpuxNTdiXJQ7k2I8L+YxhBCklCjYn17BvrQK8mpaAF31fUCAA4MCHOjra/f/qtrKDg35NS0Xb3nVLRTXt1LTmV/4uQbyd/g5x+FqbYqrjRleduYEOVsQ6GSJu43Zr/ICFQolp/JqOJmjq/38vFJ5sLMlI8KcGR7qTE93ayRJoqCmha8SS9hytpjqpnbcbcyY0teLKX28sJUbk1TcwPLD2RzOrMLSxJDHbvNn1kBfzIwNfk/qTY/eKK4xtFrBtvhiluzL1O0Z0duTZ4YF4mRpilYr2J1aznsHs8ivacHPUc5jt/gzNsoNE0Pd4KrjOTVsOVvMwfRKVBot4e5WTIjyYHRPl4tzEX5OCB7L1o1HOJ33vyaHh60ZkZ429PKwIdzdmkBnC+zlxr/bdGhVqaltVtGiUtPSrqFNpaFVpUYrQCaBTJKQycBQJsPKzKgz6WmIlanRr3IRPyOEoKqpncyKJlJLGkgpUZBaqqC8c91QG3Mj+vnqcilDgp3wstfViKqalOxOLueb5DKSixuQJBga7MTUfl7cGuSETIITObWsOpbH0axqbMyNmDXQlxkDfG7I/TUuN3qjuEapa1HxwaEsNsYVYWwo45Fb/Hh4sB9yE0M0WsHetHI++TGX9LJGnK1MmN7fh8mxnjha6moE9S0qvkkqZVt8CRnljQD08rRhdLgLo8JdfjW7tEOjJaVElyxMLtaNWfjl3iI25kb4O1rg7yjH3cYcF2tdEtLF2hR7uQmWpoaYGMq6zEMIIVB2aFG0dVDVYa5tjQAABPlJREFUpKSqsZ3KJiWVje0U1upqK/k1Lb+qpfg5yOnpoev96OdnRw8XK2QyqdPsWjiSWcnh81UXmx9hblaM7eXGXRFuuNmY0apSszOxlDUnC8ipasbBwpjZg/yY1t8bi5tk7svlQG8U1zj5NS0s25/J96kV2MuNmT3Yl2n9vLE0NbpYg1j5Uy4ncmoxMpAYEebC1L5eFxObALnVzexLq2BvWjlppTrT8LY3Z1CAA4MDHejv7/D/rqrVTe2klynIrW4ht7qZ3KpmcqtbLg54uhQDmYSFiSFyYwMkSdc1CrruzvYOLa0qDS0qNb/1NZIkcLM2w89Rjr+jBX6OcgKcLAh3t8bK9H+6apvbOZNfR1x+HT9eqKKgthXQDRIbEerM3b3cCHDSNdVSSxXsSCjh63OlNCnV9HS35sGBPtwZ8b/8jZ4/z3VvFJIkVaNbXu9y44Bukd/rhetJ7/WkFa4vvVdKq7cQwrGrQtesUVwpJEmK/zMOeq1wPem9nrTC9aW3u7Xql3vWo0dPl+iNQo8ePV1yMxrFqu4W8Be5nvReT1rh+tLbrVpvuhyFHj16/jo3Y41Cjx49f5Eb3igkSbKTJOmgJEnZnX9/d50zSZKsJEkqkSTpo6up8RINXeqVJKmXJEmnJElKlyQpRZKke6+yxlGSJF2QJCmncyvJS583kSRpa+fzcZIk+VxNfZdo6UrrPEmSMjrP42FJkry7Q+cv9Pyh3l+UmyhJkujcle+Kc8MbBbAAOCyECAQOd97/PV4Djl4VVb/Pn9HbCkwXQoQBo4D3JUm6KhtPSJJkAHwMjAZCgfslSQq9pNhsdDvbBwDvAW9dDW2X8ie1nkO3zWUEsANYenVV/o8/qRdJkiyBp4G4q6XtZjCKscDazuO1wLjfKiRJUgzgDBy4Srp+jy71CiGyhBDZncdlQBXQ5aCZy0QfIEcIkSeEUAFb0Gn+Jb+MYQdwu9Q968x1qVUI8YMQorXz7mnA4ypr/CV/5tyC7oL2FqC8WsJuBqNwFkKUdx5XoDODXyFJkgx4B5h/NYX9Dl3q/SWSJPUBjIHcKy2sE3eg+Bf3Szof+80yQgg1oAC6YyGIP6P1l8wG9l5RRX9Ml3olSYoGPIUQe66msBti9owkSYeA39oYctEv7wghhCRJv9XNMwf4XghRcjUufJdB78+v4wqsB2YIIbS/V05P10iS9AAQC9za3Vp+j84L2rvAzKv93jeEUQghhv3ec5IkVUqS5CqEKO/8YVX9RrH+wGBJkuYAFoCxJEnNQog/ymd0p14kSbIC9gCLhBCnr4TO36EU+OWW3h6dj/1WmRJJkgwBa6D26sj7TR0/81takSRpGDqTvlUI8duz464OXem1BMKBHzsvaC7At5Ik3S2EuLIzKIUQN/QNWAYs6DxeACztovxM4KNrWS+6psZh4Jlu0GcI5AG+nTqSgbBLyjwBrOw8vg/Y1k3n8s9ojULXbAvsrs/8r+i9pPyP6BKxV15bd5+cq3Dy7Tt/VNnAIcCu8/FY4PPfKN/dRtGlXuABoANI+sWt11XUeAeQ1fkDW9T52KvA3Z3HpsB2IAc4A/h14/nsSushoPIX5/Hbbv6+/qHeS8peNaPQj8zUo0dPl9wMvR569Oj5h+iNQo8ePV2iNwo9evR0id4o9OjR0yV6o9CjR0+X6I1Cjx49XaI3Cj169HSJ3ij06NHTJf8HsAyK7jN/3mIAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fix, ax = plt.subplots(1,1)\n",
"ax.set_aspect(\"equal\")\n",
"ax.set_xlim([-.5,.5])\n",
"ax.set_ylim([-.5,.5])\n",
"ax.plot(orbits[:,6::60],orbits[:,7::60])"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [],
"source": [
"def midpoint_step(y,dt):\n",
" k1 = dt*G(y)\n",
" k2 = dt*G(y+0.5*k1) \n",
" return y + k2"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [],
"source": [
"steps = 500\n",
"dt = 0.02\n",
"y = y0\n",
"orbits = np.zeros((steps,Nd))\n",
"for i in range(steps):\n",
" y = midpoint_step(y,dt)\n",
" orbits[i] = y"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x10cbe3588>]"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAD8CAYAAACPd+p5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VVXaNvD7SSekkUZJIRBC6KEEEASlSfFVYUQRkFEZHT51bOOMiOM4zifiWAbLB3ZHZRxeARUVFUWIhCYtdAKShBIIJQ1IICH1PN8fHMaQ7JLknOxz9j7P77q4ONlrkfMQwp291157LWJmCCGEFi9XFyCEcH8SFEIIXRIUQghdEhRCCF0SFEIIXRIUQghdEhRCCF0SFEIIXRIUQghdPq4uQE1kZCQnJCS4ugwhLG3Hjh1FzByl189tgyIhIQEZGRmuLkMISyOi3Mb0k0sPIYQuCQohhC4JCiGELgkKIYQuCQohhC4JCiGELgkKIYQuCQohhC4JCiGELgkKIYQuCQohhC4JCiGELgkKIYQuCQohhC4JCiGELgkKIYQupwQFEY0nokNElENEczT6TSYiJqJUZ7yvEMIYDgcFEXkDeBPABAA9AEwjoh4K/YIBPApgq6PvKYQwljPOKAYByGHmI8xcBWAJgIkK/eYCeAlAhRPeUwhhIGcERQyAE3U+zrMf+y8i6g8gjpm/c8L7CSEM1uKDmUTkBeBVAH9qRN9ZRJRBRBmFhYUtXZoQopGcERQnAcTV+TjWfuyKYAC9AKQT0TEA1wBYoTSgyczvMXMqM6dGRemuIC6EMIgzgmI7gCQi6kREfgCmAlhxpZGZS5g5kpkTmDkBwBYAtzCzrMUvhEk4HBTMXAPgIQCrABwEsIyZM4noOSK6xdHPL4RwPadsAMTMKwGsrHfsbyp9RzjjPYUQxpGZmUIIXRIUQghdEhRCCF0SFEIIXRIUQghdEhRCCF0SFEIIXRIUQghdEhRCCF0SFEIIXU6Zwi2srbrWhvKqWlRU18KLCMEBPgjw9XZ1WcJAEhQeqqS8GseKy3C48CIycs9h29GzyCm46JTPPaxLJG4bEIshiRGIDvYHETnl8wrXkaCwOJuNkV1wEduOFuOHzDPYlFPc4u+5MacIG3OKFNtenZKCsT3bIchfvvXMRP61LKa61obMU6VYn1WIhT/loKrW5uqSrvL4sj0A9vz346kD4/D3W3rKpYybk6CwgJLyaqzKPINFm48h81SpUz7ndV2jMKRzBFJiQxEXHoiwQF/4+XjBz9tL81KisqYWZ8uqkHfuEn45XYqNOUVYlZmv2n/J9hNYsv3XJVefuakH7h3WySl/B+E8xMyurkFRamoqZ2TIIlhqyiprsPpAPl5bk4Xc4vJmfY7IIH/cN7wThidFIjEqyLCf6ufLq/DyqkP4363HNftNGxSPf9za25CaPBUR7WBm3X12JChMpLKmFmt/KcTirbnYkK08BqClX3wYZgzuiMGdwxET1sqtBhlPnb+EoS/+pNo+sW8H/PP2FPh6yx19Z5KgsJBzZVX4z5ZczF+d1eQ/+9f/6Y5bUjogOiSgBSprGbU2RuJfVqq2f/WHa9E3LszAiqxLgsICcovLsPCnHHy2I69Jf+7PY7tiUr8YxLYJbKHKjFN0sRKpz69RbHtgRCKeGJsMLy/3OTMyGwkKE9uRew5vpGVjfVbj9zZ5cEQiJg+IRWJUUAtW5lqvrc7CG2nZDY6HtvJF+p9HoE1rPxdUZW4SFCa0/2QJnl2RiR255xr9ZxZM64cJvdrBx4Ou3TOOncVt72xWbNsweyTiws1/JmUUCQoTyS+twEvf/4Llu07qdwaQGNUacyf2wpDECLcakDTanhPnMfHNTQ2OB/v74OenRiE4wNcFVZmLBIUJXKqqxfsbjuDVRg5SjuoWjT+PTUaPDiEtXJm5rMsqxN0fbmtw/NZ+MXjptj5yp0SDBIUbs9kYK/acwmNLdzeq/9DECLx4ax/ER8gptZZ/bz6Gv32d2eD43Ik9MeOajh599qVGgsJNnTp/CX9atgebj+g/c+HrTfh45iBc2yXSgMqsoabWhutfScfJ85catH3z0DD0jg11QVXuS4LCzTAzvtx10v6sg76/39wDM67p6FGDlM60/2QJblqwscHxKamx+MetfeAtt1QBSFC4lbNlVZj9+R6sOVig23dKaiyeHN8NEUH+BlRmbTW1Nty7KAPrFG4zb3xypCXmmThKgsJNpB3Mx72L9P8eXgR88cBQ9ItvY0BVnkXt7shf/6c77h3WyaPHLiQoXKy61obnvjmAT7bk6vad1LcDnv9Nb1mjoQXV1Now4Pk1KLlUfdVxImDXMzcgLNAzJ2s1NijkArgFnC2rwp0fbG1USLwxtS9en9pPQqKF+Xh7Yc+zY/HYmKSrjjMDfZ9bjS2NGFz2ZBIUTpaVfwEj/5mObUfPavbr1i4YG2aPxMS+MQZVJgDgsTFdsWTWNQ2OT31vC5Zu137s3ZNJUDjRmgP5GPva+gant/U9MjoJ3z48TKYau8g1nSOw8cmRDY4/+cU+/H1FJtz1ctyVJCicgJnxdvph3Pdv/TGVt+7sj8dv6Cq3PV0stk0gMv/vuAbHP/75GG57ZzOq3WwJQVeT71YH2WyMv361Hy/98Itu3y8eGIIbe7c3oCrRGK39fXDkhRvRo/3VU+J35J5Dt2d+wIUK7TNDTyJB4YBaG2PO8r1YrLOkW5vAy49BD+gYblBlorG8vAjfPTIM13WNuup4rY3R++8/4nRJwxmensgpQUFE44noEBHlENEchfbHiegAEe0lojQi6uiM93WlWhvjic/2YFmG9qIyKbGhSPvTCCREtjaoMtFURIRFMwdiTPfoBm1D/vETzpRUuKAq9+JwUBCRN4A3AUwA0APANCLqUa/bLgCpzNwHwOcAXnb0fV2pptaGPy7drftY+I2922Hp/xmCcFlQxe0REd6/KxXje7Zr0HbNP9JQfLHSBVW5D2ecUQwCkMPMR5i5CsASABPrdmDmtcx8ZanoLQBinfC+LlFda8OjS3ZjxZ5Tmv1u6tMeC6b1l/0qTISI8Nad/TGkc0SDtgHPr0GpB49ZOCMoYgCcqPNxnv2YmnsBfO+E9zWczcZ4bOlufLfvtGa/Md2j8dodfeXBIxPy8iIsvm8wEqMaXioOmLsal6pqXVCV6xk6mElEMwCkAnhFpX0WEWUQUUZhYePXizTKSz/8gu/2aofEsC6RWDi9vyyWYmJeXoQf/3g9WtU7G6yuZYyan46qGs+7deqM7+aTAOLqfBxrP3YVIhoD4GkAtzCz4gUfM7/HzKnMnBoVFaXUxWUWb83Fu+uPaPYZmNAG79+VKpcbFuDtRdj5zA0Njp8uqcCtb29Crc2zJmU5Iyi2A0giok5E5AdgKoAVdTsQUT8A7+JySOg/a+1m0g8V4Okv92v26RMbio9mDkIrPwkJq2jl5431TzScwbn/ZCn++pX294PVOBwUzFwD4CEAqwAcBLCMmTOJ6DkiusXe7RUAQQA+I6LdRLRC5dO5nQOnSnHPR9s1+8S2aYVPfjdYHuyyoPiIQHw0c2CD459uO44f9p9xQUWuIY+ZazhTUoFbFm5EwQXtW2NrHr8eXaKtu5+GuDw+9Xb64QbH0/9s7jky8pi5g2ptjEeW7NINiffvSpWQ8ABPjE1GTFirBsdH/DMdFdXWvxMiQaFiwU/Zuo+KPzYmCTf0aGtQRcKVvLwIKx8drtj26JJdBldjPAkKBVuPFOP1NQ23rqvrhh5t8cioJM0+wlpCW/ni24eHNTi+KjMfy7afUPgT1iFBUc+5sir84X+1f0J0imyNV6ekyOa4HqhXTCgm9284sXj2F3tx6MwFF1RkDAmKOpgZs7/YiyKdef3v/XaAbFfnweZO6ql4fPr7W2Cz6PwKCYo6Fm89jtUH8jX7zJnQDUltgw2qSLijQD8ffKxwy7S4rAqfWnQ5PQkKu/zSCsz77qBmn5S4MNw3rJNBFQl3NiI5usEaFgDw9Jf7Uahzp8yMJCjs5n57AJd0bnPNv72PLGEn/mv+7SmKx59avtfgSlqefNcD2JBdiG91HvZ6cnw3dImWSw7xq6hgf8z7Ta8Gx9ccLMCGbPd7qNERHh8UFdW1uvP2U+LC8PvhcskhGpo+KF7x+G//tc1SE7E8PijeWXcYucXlmn1eniyXHEIZEeHDe5RnQL+zruGUb7Py6O/+48XleCNNe2LV9MHxSG4nlxxC3cjkaAQrPBD4+ppslJRbY1Usjw6K19OyoPdMXP0t6ISoj4jw2h19Fds+3HTU4GpahscGRU7BRSzfqb047mNjkhAdHGBQRcLMRnePRlhgw0l4b6RlW2KtTY8NCr1LjsggP/x+eGeDqhFmR0R45Tbl26WLNh0ztpgW4JFBcejMBXyjs4r2n8Ymo7UsRCOaQGlfEACYvzrL9LuOeWRQvLY6S7O9S3QQbh9g2h0FhIsQEd6+s79i27835xpcjXN5XFBknirBD5naS5jdf32i3A4VzTJGZX2SV1YdQnlVjcHVOI/H/W/4WOd6MTLIHzenyEbConl8vb3wu2uVJ+et3GfeNTY9KijOlVXh693aYxP3DO0Ifx9ZSVs03x0D4xSPf7jRvLdKPSooPttxAlW16pu3eHsRpg82/f7JwsWS2wUjPjywwfEDp0txrKjMBRU5zmOCwmZj/GeL9loBdwyMkw2FhVOoPRv0+Y48gytxDo8JinXZhTh+VvuZDrVrSyGa6uaUDorHF67NMeUuYx4TFJ/o3J4a1iVSlt0XThMW6IfhSZGKbWZ8BN0jgqLoYiXWHtLeyXBSP60N2IVoujsHKz+CvtSEK3Z7RFD8mJmv+fCXtxdhXE/Zn0M415BE5TOK7/efMd1aFR4RFN/v11696obubWVVbeF0oa180TlKebvBrTqbS7kbywfF2bIqbMgu0uwzsa/ywJMQjrqpt/LkvXWHzDVOYfmgWH1AezZckL8PRnZTfphHCEcNS2q4UjcAfLNXe+Kfu7F8UOhNmx3bsy0CfGUmpmgZ/eLDFI8XXqhE3jnt2/XuxNJBUVZZg0052pcdY7rLIKZoOb7eXhiaGKHYtj5L+3vTnVg6KHbknkONzuSWIZ2V/xGFcBa1He/XZWnfsncnlg6KrUeLNdt7dghBG5myLVpYn1jly49Vmfmm2avU0kGx5Yj2LSi1U0IhnKlrW/UZv7k6jxW4C6cEBRGNJ6JDRJRDRHMU2v2JaKm9fSsRJTjjfbWUV9Vgz4nzmn2GdlGeECOEMwUH+MJPZSGkg6dLDa6meRwOCiLyBvAmgAkAegCYRkQ96nW7F8A5Zu4C4DUALzn6vnp25p7XHJ/w8SIMSghv6TKEAAAMUTl79ZigADAIQA4zH2HmKgBLAEys12cigEX2158DGE1E5IT3VrVNZ3wiJS5MFs8VhlHbROrAKc8JihgAdZ9yybMfU+zDzDUASgC06ADBgdMXNNt7dghpybcX4ipJKk8mbzPJVG63GswkollElEFEGYWFjk1xzcrXDopu7SQohHHUziguVNbgfHmVwdU0nTOC4iSAuosExtqPKfYhIh8AoQAaXBsw83vMnMrMqVFRylNfG6O8qgYndGa9yX6iwkhaa50c1Dn7dQfOCIrtAJKIqBMR+QGYCmBFvT4rANxtf30bgJ+Y9Xb9bL6cgou6e4pKUAgjBfqpj4cdNcE6mg6P5jFzDRE9BGAVAG8AHzJzJhE9ByCDmVcA+BeAT4goB8BZXA6TFpOVf1GzPS68FYJkIFMYzN/HC5U1DRd3LrpY6YJqmsYp/1uYeSWAlfWO/a3O6woAtzvjvRojW2d8IrmtjE8I43WJDkKmwl2OwgvuHxRuNZjpLHqL6CaqLCYiREtSW+FdgsJF9L7w0SEBBlUixK/CAlWCwgSXHtYMCp0vfHSwv0GVCPGr8EDl5RbPlFQYXEnTWS4omBkFpRIUwv2onVGcPH8JLXgT0CksFxRlVbW4pLPCsVx6CFcIDlC/d6D3PetqlguKxgwMyRmFcAWtHcKqa+WMwlBny7Snw7b285aHwYRLaD3NLJceBqvW2K0cAAIlJISLaJ1RuPtCVx4XFGoLiAjR0mo0vjdtckZhrBqdaz0f7xZdBkMIVVqXHhIUBqvSOaPw8ZKgEK6hdenh5jlhvaDQO6PwlUsP4SJadzbkjMJgNTadMwq59BAuUlpRrXi8R/sQBPi49251lgsKref+AcCrZZfqFEJVfqnyVO0P7k51+/1lLBcUrf21k/liZY1BlQhxNbVnOtSeKnUnlgsKvQVpLlRIUAjXOFasvJKVGTbJtlxQ6M26vChBIVygrLJGcTCzY0SgC6ppOssFhd4ZxaXqWt1JWUI42xmV8YnIIHM8d+RxQQFcTnchjJSvMj4RYYLxCcCCQRHo541gnbAovSRBIYx18vwlxeORJnmS2XJBQUSIadNKs0/+BfdfUUhYywGVPUY7R5pj/VbLBQUAxLbRHiA6XmyOreaFdWSeVA4Ks+wvY9Gg0D6jyNVZpVsIZ7LZGDuPn1NsS24rQeEyekFxQoJCGOhYcZnik6OBft6IkjEK19G99JCgEAbar7DpDwD0igkFmeSRAksGRTed6z4JCmGkzJMlisfNctkBWDQoOkYEaq54XHih0hRbzQtr2KcSFF1NMpAJWDQoiAi9OoRq9tmbp/yPJ4QzXaqqxc+HixXbBia0Mbia5rNkUABArxjtjYj3nDhvUCXCk205ohwSrXy90TVazihcrleM9hnFnjwJCtHy0g8VKB4f1S0aXiZaltGyQdFbJyh2nyhx+70UhPmtPVSoeHxIYoTBlTjGskHRKbI12mlsHVh0sRKnTLA5rDCvo0VlqnfYhkpQuAciwvCkSM0+GcfOGlSN8ERqlx1tQ/zRySTPeFxh2aAAgOFdozTb1/6i/A8phDOoXXZcmxhpmolWV1g6KK7VOb1bl1WoudeCEM1VcKEC67OUg+LG3u0NrsZxDgUFEYUT0Woiyrb/3uDGMBH1JaLNRJRJRHuJ6A5H3rMpIoL8NW+TniuvlrsfokV8ufOk4nEvAoZ31b4kdkeOnlHMAZDGzEkA0uwf11cO4C5m7glgPIDXiSjMwfdttBFdozXb0+XyQzgZM+OzHXmKbbcNiIW/m+/hocTRoJgIYJH99SIAk+p3YOYsZs62vz4FoACA9uCBE93St4Nmu9p1pBDNtSevBDkFFxXbbk7R/n50V44GRVtmPm1/fQZAW63ORDQIgB+Aww6+b6N1bRus+ZDYvpMl8ti5cKrPMk4oHg8L9MWQzua6LXqFblAQ0Roi2q/wa2Ldfnx59pLqyCARtQfwCYCZzKy4DDYRzSKiDCLKKCx03k/6iX1jNNu/3KV8PSlEU1VU12Lx1uOKbTf36QAfk+59q1s1M49h5l4Kv74GkG8PgCtBoHjBT0QhAL4D8DQzb9F4r/eYOZWZU6OinHd1cnOK9ijz8p15MktTOMV3e0+rtk0eEGtgJc7laLytAHC3/fXdAL6u34GI/AB8CeDfzPy5g+/XLLFtAjWf1DtWXI4ducpLlQnRWDYb4+11ylfV/ePD0DfOsDF8p3M0KF4EcAMRZQMYY/8YRJRKRB/Y+0wBcB2Ae4hot/1XXwfft8nuGBiv2f7FTuVRaiEaa/XBfNVBzJnXdjK4GudyKCiYuZiZRzNzkv0S5az9eAYz32d//R9m9mXmvnV+7XZG8U1xc0p7hLbyVW3/ctdJXKqqNbAiYSXMjLfSlc8mIoP8Mb5XO4Mrci5zjqw0g7+PN+4emqDaXlFtw/JdclYhmufnw8Wqa5z8blgCfE06iHmFuatvohmDtS8/PthwVKZ0i2Z5c22Oats0ncteM/CooIgOCcAkjQlYR4vKsPpAvoEVCSvYefyc6nJ3M66JRxuT7C+qxaOCArg8qOTrrf7k3nvrDZsLJiyAmTHvu4Oq7Q+PSjKwmpbjcUGREheG7Hk3YlxP5UmkO4+fl3UqRKOt2HNK9db6/dcnoq3G4klm4nFBccVjY7qqti3UuN4U4opLVbV4XuVsggh44PpEgytqOR4bFN3bh+DG3sq3rNIPFeLnnCKDKxJm8866wyi8UKnY9sS4ZIQGqt+ONxuPDQoAeHS0+lnFC98fhE3ugAgVJ89fwhtp2Ypt4a39MHOouSdY1efRQZHcLhh3pMYptu0/WYpv9p4yuCJhFi+sVB/AnD0uGa38zLfmhBaPDgoAmD0+WbXtxe9/QUW1zNYUV1u577Tqw1+DEsIxReWHj5l5fFBEBPlj7sSeim2nSyrw0aZjxhYk3FpBaQUeX6b+BMILt/Yy1cY+jeXxQQEA0wd3VN1Z+qUffsGxojKDKxLuiJkxZ/k+VFQrLqeCx8YkoYuJtglsCgkKAN5ehBdu7a3a/tTyfbJehcDS7Sfwk8oaq50jW+OBEda5HVqfBIXdgI5tcI/KQ2ObjxRj6Xbl5c2EZzheXI45y/eptr84uY8pF81tLAmKOuZM6IbYNq2U25bvwxnZgtAjVdbU4uElu1TbZ13XGYM6hRtYkfEkKOoI8PXGOzMGqLY/tXyvXIJ4GGbG377KVH2EPCU2FE+MU79zZhUSFPX0ilH/h197qBAfbDhqcEXClT7ZkoulKqtqA8DC6f1Nv9ZEY1j/b9gM91+fiH7xyusbzlt5UNbX9BBbjhTjb19nqra/+9sBiAsPNLAi15GgUODtRXj3t+qXIPct2o5zZVUGViSMlneuHNPeV10wHjOvTcC4nuZe3q4pJChURAcH4Os/XKvYdq68Gn9ctlueBbGo0opqzPxoO9SGowYlhGPOhG7GFuViEhQaUuLCMP/2FMW29EOFmL/6kMEViZZ2qaoWv/toO7JVVtPuEBqAd387wNK3QpVIUOiYPCAWM69NUGx7c+1hfLpNeVcoYT6VNbWY9UkGMjTGoJbMGmKJpe2aSoKiEZ6+sTtGJivvXPbU8n1Ye0h2RDe7mlobHv10NzZkq69DsvzBoYiP8IzBy/okKBrBx9sLb88YgJTYUMX2mR9tx/6TJQZXJZzFZmM8+cU+/JB5RrXP23f2R/949d3mrE6CopECfL2x+PfXICZMeebmTQs24qg8PGY6NbU2zFm+V3OnuH/c2hsTemvvX2t1EhRNEOTvg28fHgY/H+Uv28h/puNIofIgmHA/FdW1uP8/O7EsQz0kXprcG9MGmX9fDkdJUDRRm9Z+2DB7JPxUZuONmr8OhyUs3F7JpWrM+GAr1hxU38fl5cl9dPes9RQSFM3QNiQAW/4yGhEqo9+j569T3axWuF5+aQUmvblJ8+7GK7f1wZSB1lupqrkkKJopvLUf0p8Ygc6RrRXbx7y6TgY43VBW/gWMnr9Oczzp1SkpuN2Cy9k5QoLCAcEBvlj56HD0jVN+LuSmBRuRpnFqK4z1zZ5TGPvaelysrFHts/i+wbi1f6yBVZmDBIWDAny98dn9QzAlVfmb695FGVj08zFjixJXqa61Ye63B/Dwp+prSgT4emH1H6/DtV0iDazMPCQonMDX2wsv35aCF1WW03t2RSae+Wq/7JTuAgUXKnDHu5vxr43qywP0aB+CDbNHIUll3VQhQeFUUwfFY/mDQxXbPtmSiwlvrFfdWUo436acIgyal4adx5UXnQGAm/q0x/IHhyIq2N/AysxHgsLJ+se3wda/jFa8fZqVfxED562R7Qpb2MXKGjz95T7c+cFWzX7zftMLC6b1Q4CvZz3g1RwOBQURhRPRaiLKtv+uOseViEKIKI+IFjrynmbQNiQAB54bh1nXdVZsn/7BVrz8wy9yKdICfj5chIHPr8HireoP64W28sX3jw7HnYM7gsh6e3C0BEfPKOYASGPmJABp9o/VzAWw3sH3Mw0fby/85cbu+EplTYu30g9jwPOrcby43ODKrKmssgbPfLUf09/fiksau7tN7h+LzU+NQvf2IQZWZ36OBsVEAIvsrxcBmKTUiYgGAGgL4EcH3890+saF4cBz4zA8qeFo+vnyalz3ylq8uTYHNbXKm8oIbTYb44sdeej57Cp8siVXs+//m9YP86ekINDPx6DqrIMcWVWaiM4zc5j9NQE4d+XjOn28APwEYAaAMQBSmfkhvc+dmprKGRkZza7NHW05Uoyp76kvr/btw8PQK0b5CVXRUMaxs5j9+V4c0XkY7zf9YvDszT0QFuh560joIaIdzJyq1083WoloDQClxQGfrvsBMzMRKaXOgwBWMnOe3vUgEc0CMAsA4uOtN8f+ms4RyJk3AW+kZWPBTzkN2m9asBE39WmPeZN6IzTQ1wUVmsOJs+WY++0B/HhAezKbv48XPpo5EEMTZW6Eoxw9ozgEYAQznyai9gDSmTm5Xp/FAIYDsAEIAuAH4C1m1hrPsOQZRV0FFyowaeEmnFLZVOixMUm4//pEGZGv41hRGd5Kz9F82vOKR0Yn4cER8vXT09gzCkeD4hUAxcz8IhHNARDOzLM1+t8DD770UHLwdCkmvLFBtX3eb3ph6sB4eFtwh+zGOni6FK+vycKqTP3p8LcNiMXsccmIDgkwoDLzMyooIgAsAxAPIBfAFGY+S0SpAO5n5vvq9b8HEhSK9pw4j4lvblJtf/bmHpg+ON5jFnW12RhbjhTjtTVZ2H5Mfx+VkclReOamHugcFWRAddZhSFC0JE8Liiu2HzuL29/ZrNp+z9AE/HFMV8uOYRSUVuCzHXl4ZVXjVjgfmhiBJ8d3Q4rKg3lCmwSFyeWXVuCuf23DofwLiu3tQgLw/l2p6BUTYvpJQ5U1tdiQVYSFa3OwW2WPz/pu7R+DR0YlIUHlMX/ROBIUFlFVY8Pra7LwVvph1T5TB8bhT2OTTfW8Qkl5NdYeKsCn245j69Gzjf5zT4xLxp2D4+VWp5NIUFhQfmkFfv/vDOzNU18QZ3L/WDwxLhntQt1rMK+m1oaDpy9g69FiLPgpByWXqhv9Zwd1CsfDo7pgaGKkRw/qtgQJCovLKbiImxds1JyuDAALp/fDhF7tDf8PVnKpGpknS7A+uwgfbjqKqpqmzTyNCvbHUxO6YVzPdmjtLzMpW4oEhQepqrHhhZUH8XEjF8iZM6EbJvePRWSQX7PHN5gZFyprUHShEmdKK3DgVCnSDhZg85Euf1BYAAAEVElEQVTiZn0+AOjePgQPjeyCYUmRCG1lzcFadyNB4cEqqmvx16/24/Md+hOTXCnQzxv3De+MsT3aokf7EHjJZYXhJCjEVc6VVeGBxTuw5UjjBw6dbUz3tpjUrwP6x7dB+9AA09+tsQKnPeshrKFNaz8smTWkwfGqGhs2HS7Cx5uOYV1WocPvMyghHP3iwzCoUzg6RwUhJqyV6oZJwjwkKDycn48XRiZHY2RytKtLEW5Mol4IoUuCQgihS4JCCKFLgkIIoUuCQgihS4JCCKFLgkIIoUuCQgihy22ncBNRIS4vr+dskQDMtKefmeo1U62AueptqVo7MnOUXie3DYqWQkQZjZnb7i7MVK+ZagXMVa+ra5VLDyGELgkKIYQuTwyK91xdQBOZqV4z1QqYq16X1upxYxRCiKbzxDMKIUQTWT4oiCiciFYTUbb99zYafUOIKI+IFhpZY70adOslor5EtJmIMoloLxHdYXCN44noEBHl2LeSrN/uT0RL7e1biSjByPrq1aJX6+NEdMD+dUwjoo6uqLNOPZr11uk3mYjYvitfi7N8UACYAyCNmZMApNk/VjMXwHpDqlLXmHrLAdzFzD0BjAfwOhEZslUWEXkDeBPABAA9AEwjoh71ut0L4BwzdwHwGoCXjKitvkbWuguXt7nsA+BzAC8bW+WvGlkviCgYwKMAthpVmycExUQAi+yvFwGYpNSJiAYAaAvgR4PqUqNbLzNnMXO2/fUpAAUAdCfNOMkgADnMfISZqwAsweWa66r7d/gcwGhyzQKZurUy81pmLrd/uAVArME11tWYry1w+QfaSwAqjCrME4KiLTOftr8+g8thcBUi8gIwH8CfjSxMhW69dRHRIAB+ANS3EnOuGAAn6nycZz+m2IeZawCUAIgwpDqVOuyUaq3rXgDft2hF2nTrJaL+AOKY+TsjC7PEmplEtAZAO4Wmp+t+wMxMREq3eR4EsJKZ84z4weeEeq98nvYAPgFwNzM3bYcdcRUimgEgFcD1rq5Fjf0H2qsA7jH6vS0RFMw8Rq2NiPKJqD0zn7b/xypQ6DYEwHAiehBAEAA/IrrIzFrjGa6sF0QUAuA7AE8z85aWqFPFSQBxdT6OtR9T6pNHRD4AQgE0f2eg5mtMrSCiMbgc0tczc6VBtSnRqzcYQC8A6fYfaO0ArCCiW5i5Zfe2YGZL/wLwCoA59tdzALys0/8eAAvduV5cvtRIA/CYC+rzAXAEQCd7HXsA9KzX5w8A3rG/ngpgmYu+lo2ptR8uX7YluerfvCn11uufjssDsS1fm6u/OAZ88SPs/6myAawBEG4/ngrgA4X+rg4K3XoBzABQDWB3nV99DazxRgBZ9v9gT9uPPQfgFvvrAACfAcgBsA1AZxd+PfVqXQMgv87XcYWLv181663X17CgkJmZQghdnnDXQwjhIAkKIYQuCQohhC4JCiGELgkKIYQuCQohhC4JCiGELgkKIYSu/w/whaynNevJMwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fix, ax = plt.subplots(1,1)\n",
"ax.set_aspect(\"equal\")\n",
"ax.set_xlim([-.5,.5])\n",
"ax.set_ylim([-.5,.5])\n",
"ax.plot(orbits[:,6::60],orbits[:,7::60])"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [],
"source": [
"def rk4_step(y,dt):\n",
" k1 = dt*G(y)\n",
" k2 = dt*G(y+0.5*k1) \n",
" k3 = dt*G(y+0.5*k2) \n",
" k4 = dt*G(y+k3) \n",
" return y + 1./6.*(k1+2.*k2+2.*k3+k4)"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [],
"source": [
"steps = 500\n",
"dt = 0.02\n",
"y = y0\n",
"orbits = np.zeros((steps,Nd))\n",
"for i in range(steps):\n",
" y = rk4_step(y,dt)\n",
" orbits[i] = y"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x10cc62e48>]"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAD8CAYAAACPd+p5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHoRJREFUeJzt3Xd8FVXaB/DfEyD0Ir2poUQQEBEjWOhFQBZwd30tq4soyOqurg3dAFYQBQu6a9tl1V3Wjrr6oiAtygsWSpCmtCAgRCmhSDWBJM/7Ry5z597cZJLccu7c+X0/Hz6cmXuSeQjJL3NnzpwjqgoiotIkmS6AiOIfg4KIHDEoiMgRg4KIHDEoiMgRg4KIHDEoiMgRg4KIHDEoiMhRZdMFlKRhw4aakpJiugyihLZq1ar9qtrIqV/cBkVKSgoyMzNNl0GU0ETkh7L041sPInLEoCAiRwwKInLEoCAiRwwKInLEoCAiRwwKInLEoCAiRwwKInLEoCAiRwwKInLEoCAiRwwKInLEoCAiRwwKInLEoCAiRxEJChEZLCKbRWSriKSX0u+3IqIikhaJ4xJRbIQdFCJSCcCLAIYA6ADgOhHpEKJfbQB3Alge7jGJKLYicUbRDcBWVd2mqicBvANgRIh+kwFMA5AbgWMSUQxFIihaANhl28727bOISFcAZ6rqnAgcj4hiLOoXM0UkCcB0APeWoe9YEckUkcycnJxol0ZEZRSJoPgRwJm27Za+fafVBtAJwGIR2QHgYgCzQ13QVNUZqpqmqmmNGjnOIE5EMRKJoFgJIFVEWolIMoBrAcw+/aKqHlbVhqqaoqopAJYBGK6qnIufyCXCDgpVzQdwO4D5ADYCmKWq34nIJBEZHu7nJyLzIrIAkKrOBTA3aN9DJfTtE4ljElHscGQmETliUBCRIwYFETliUBCRIwYFETliUBCRIwYFETliUBCRIwYFETliUBCRIwYFETliUBCRo4g8FEaJTVVxLC8fX39/AF9s3Y8V2w9i056jAX2SKyXhsrYN0CO1EYZ0aopmdatBRAxVTJEmqmq6hpDS0tI0M5NTVsTK0qwc/P7VFTE51therXHfoHaoUokntKaJyCpVdZwVn2cUHrTguz0Y+/oqY8efsWQbZizZFrDvyi7NMe2qzqhauZKhqqg0DAoPWLPrZ1z54pdl7v/IsA64tttZqFYlMj+0hYWKr7cdwPWvlLxSw0drfsJHa36ytidecS5u6dU6Isen8PGtR4KaNm8TXl78fal9rrqwJab+5jxUNvQW4JeTBfj3Vzswbd6mUvu1blgTC+7uZazORFbWtx4MigTy/qpsjHtvbYmvj7v8HNzWpy0qJcXnRcbCQsUn63fjz2+vLrHPby5ogWeuPp8XSiOEQeERh385hfMfXVDi68vG90fTutViWFHk/HDgOHo/tbjE15fe3xdn1q8Ru4ISEIMiwc37dg9ufSP0BclF9/RG28a1YlxRdGUfOoEe0z4P+dp9g9rhT33bxriixMCgSFCTP9mAV7/YXmy/V35YCgsVf1/yPZ6ct7nYaxelnIFZf7iEb0vKgUGRYP789mrMXvtTsf1fpfdD83rVDVRk3tZ9RzFg+pJi+3umNsTro7sbqMh9GBQJ4qn5m/Di58XvXnz76CDUqsq72wBw4mQ+Ojw0v9h+r5xlhYNB4XIljX3YNHlwxMY3JJojuafQ+ZHiF3Y/vr0HzmtZ10BF8Y9B4VKnCgqROvHTYvsZEGW36+AJ9Hyy+IVPfg2LY1C40JiZK7Fo476AfYvH9UFKw5qGKnK3JVtyMPK1wOdX7hqQirsGnGOoovjDoHCRvUdy0f3xjIB9d/ZPxd0D+Q0dCaP+tQKLN+cE7Nv82GA+VwIGhWvc/O+V+GxT4FkEv4kjL+doHi6asihg36QRHTHykhQzBcUJBkWcC3Ut4q0x3XFp24aGKkp8qoo73l6NT9btDti//YkrPDv2oqxBwadsDFialVMsJLKmDGFIRJmI4IXfdcVn9/YO2N9q/FzsPZJrqCp3YFDE2MDp/xcwQcwd/dpix9ShnMQlhlo3qoVtj18RsK/74xn4dP3uEj6C+N0ZI/kFhUhJn4OsfcesfUvu64t7L29nsCrvSkoS7Jg6FBOuaG/tu+3NbzBmZuK+3Q0HgyIGjuflo23QW43vH78CZzXgk4+mje3VBvPu6mltL9q4FynpcwxWFJ8YFFG2Y/9xdHzYP7y4U4s62DF1aNzOCeFF7ZvWwZqHBgbsS0mfg3i90G8CgyKKvsjajz5PL7a2H/xVB3xyR8+SP4CMqVcjGVseGxKwr9X4ucjLLzBUUXyJSFCIyGAR2SwiW0UkPcTr94jIBhFZJyIZInJ2JI4bz95avhM3vOqfI/KtMd0xukcrgxWRk+TKSdj+ROBFznYPzGNYIAJBISKVALwIYAiADgCuE5EOQd1WA0hT1c4A3gfwZLjHjWd/y8jChA/XW9sL7+7FW58uIVJ0kbN+zWRrX7sH5uFUQaHBqsyLxBlFNwBbVXWbqp4E8A6AEfYOqvq5qp7wbS4D0DICx41L0xduwfSFW6ztZeP7I7VJbYMVUUV88+BANK3jn0IwdeKnKCj07jWLSARFCwC7bNvZvn0lGQ2g+OORCeBfX27H3zKyrO21D1/u2vkqCVg2oX/AmUWbCXNR6NGwiOnFTBG5AUAagKdKeH2siGSKSGZOTk6oLnFr3rd78OjHG6ztbx4ciLrVqxisiCLhmwcD74a0njDXk3dDIhEUPwI407bd0rcvgIgMADARwHBVzQv1iVR1hqqmqWpao0aNIlBabKzeeShgotuvx/cL+E1E7hZ8gbPV+LmGKjEnEkGxEkCqiLQSkWQA1wKYbe8gIhcA+AeKQmJfiM/hWnsO5+LXL31lbX96Z080q+vNOSwTlYgUG/J91ctfldA7MYUdFKqaD+B2APMBbAQwS1W/E5FJIjLc1+0pALUAvCcia0RkdgmfzlV+OVmAi5/wzyPx75suwrnN6hisiKIlKUmwdYp/nEXmD4fw0epiJ84Ji4+ZV5CqBpyCpg9pj1t7tzFYEcVC8CRDKyb0R+M67r1gzcfMo8weEm0b12JIeESTOtXw6o3+n6tuj2d44k4Ig6ICnrfdAgWKVuYi7+h/bhN0b1Xf2m49IfEvbjIoyumHA8fxjG1A1abJgw1WQ6a8M/bigO0XP99qqJLYYFCUQ2GhBiyaO/+uXpz+3aNEBN8+Osjafmr+5oSeJYtBUQ72U8ybL2uFdk05NNvLalWtjNdHd7O2g2dSTyQMijL63zWBt8IeGhb83Bt5Uc/UwIGBE20PAyYSBkUZnCooxJ3vrLG2Nz/G6xLkZ5/H4s3lO3EsL99gNdHBoCgD+4zZf722C9fcoADJlZPwykj/LdNODxdfMNntGBQOvsjaH7A9oktpD8aSVw3o0CRgO9Fm9GZQlEJVA2apWvfI5QaroXhnn3fztje/SaiBWAyKUox48Uur/ae+bVCnGh8bp5LVq5GMX1/gP+O8eeZKg9VEFoOiBIeOn8S67MPW9n2D2pfSm6jI9KvPt9qLN+ckzHybDIoSXDB5odW2r/tAVBoRwZO/7WxtXzb1c4PVRA6DIoQfDhwP2G7flI+OU9ldfZF/Hqf9x/IS4qyCQRGCfZj2yokDzBVCrvX3G7pa7bTJiwxWEhkMiiDrbdclAKBR7aqGKiE3G9ypmdU+mpfv+rMKBkWQYS98YbWDl5kjKg/7hc2rXv7aYCXhY1DYbNx9JGC7Xg1OkEsV95uu/uVr1v942NWLCDEobIb8danVXs/BVRQBf+7X1mrPWLLNYCXhYVD4/HziZMB2bQ6uogi4c8A5Vvup+ZsNVhIeBoVPl0n+cROLx/UxVwgllEpJgqqV/T9m2YdOlNI7fjEogGJj8lMa1jRUCSWihXf751TtMc2dA7AYFACeW+SfA/OJ35xnsBJKRGc1qBGwHa9LZJSGQQHgb5/5J0a9Ju3MUnoSVcyYHq2s9sfr3PcIuueDYtdB/3vGejWqIClJDFZDiWrcoHZW+89vrzZYScV4Pijsj5IvuLuXwUookQXP1p7vsjEVng+Kg8f9t0Ub13bv0nAU/+4f7D+rWLhhr8FKys/TQWF/SnRg0FRmRJE26tIUq33bm9+YK6QCPB0Uw573P9cx5dedDFZCXlAjubLpEirM00FxJNc/rTrfdlAsXHj2GVbbfiE93nk2KA4cy7PabRpxgBXFxqQRHa325E82GKykfDwbFFPmbrTaM2xrMhBFU4dm/tnSFrjogqZng+K/3/iXCGzTqJbBSshLRALH6bhllKZng4LIlJ6pDa327sPuWAE9IkEhIoNFZLOIbBWR9BCvVxWRd32vLxeRlEgct6LsT/D9/uKzDVZCXnSX7dHz9zKzDVZSdmEHhYhUAvAigCEAOgC4TkSCl/oeDeCQqrYF8CyAaeEeNxzTF/gfArvZNgafKBY6tfBfp3jW9kBiPIvEGUU3AFtVdZuqngTwDoARQX1GAJjpa78PoL8Ev1mLof+u9l+fSAl6so8o2ty4yHUkgqIFgF227WzfvpB9VDUfwGEADSJw7LAZzCsi14iri5kiMlZEMkUkMycnJyrHOHHSP8iqdjX3jpQjd2vXpLbVPp6XX0rP+BCJoPgRgH0Sh5a+fSH7iEhlAHUBHAj+RKo6Q1XTVDWtUaNGESituE17jlrt9CFcT5TMuOmyFKu9ZtfP5gopo0gExUoAqSLSSkSSAVwLYHZQn9kAbvS1rwLwmRq6gfzuCv+7pD7tGpsogQi92/l/ES7aGP8Dr8IOCt81h9sBzAewEcAsVf1ORCaJyHBft1cBNBCRrQDuAVDsFmqsvJvpD4pmdfh8B5lhf7bIDbdII/ImXVXnApgbtO8hWzsXwP9E4liRxNmsyJRKtu+9Yx65RkFECc5TQeGWcfVE8cZTQXHoxCmrfXHr+gYrIXIXTwXFweP+OSj6tecdD6Ky8lRQ7Dr4i9U+r0U9g5UQAQ1rVbXauacKDFbizFNBsWH3EavduE7VUnoSRV/LM6pb7fzC+L5+5qmgWGsbAVevOlcrJ7Na2IKiMM4vtHsqKLbuO2a16zIoyLA61fzfgxrn6wF5Kii27fev41G5kqf+6RSHqlb2fw/yjIKIHBUwKIgoFPudDp5REFFIR23PeMR5TjAoiEz56Wf/uB779Yp4FN/VESWwHbaL6/VqJBusxBmDgsgQ+7NH8c5TQVEj2X2zHxPFA08FRaptQtO8/PgeW08UTzwVFOc29QfFoePuOe0jMs1TQXFey7pWe88Rd6z5SImvY/M6zp0M81RQ2Fct32R7kpQo1gptT4sOO7+5wUrKxlNBkdKgptVe9cMhg5WQ1/1oG0PRo23DUnrGB08FRYNa/nvVn367x2Al5HUrth+02m0b1yqlZ3zwVFBUsT0x6oYp0ilx/XPpNqtdrUr837b3VFAQxQv70pZuwKAgIkeeC4qz6tew2kdzOZaCYs++vszvup9lsJKy81xQ3GxbRZp3PsiELXv9UzKO7tHKYCVl57mgGNixqdWe+dUOc4WQZz2zYLPVbt2wZik944fngqJ5Xf8q0p9vzjFYCXnVgg17rbaIOxbK9lxQuOU/hhJf5ST3fC96LiiCHed4CoqhXQdPWO3nru1isJLy8WRQ/NX2HzSPIzQphqbM2Wi1B5zbxGAl5ePJoBhku6B573trDVZCXjPvO/8vJjeMyDzNk0Hhpv8gShz5Bf7lwHqmxv+DYHaeDIpgv5zkbFcUfR+v+8lqP3lVZ4OVlF9YQSEi9UVkoYhk+f4+I0SfLiLytYh8JyLrROSacI4ZKa/emGa1n1u0xWAl5BV3v+t/m9usbvVSesafcM8o0gFkqGoqgAzfdrATAEaqakcAgwE8JyL1wjxu2Pq0a2y1/7FkWyk9icJXYJuo5owa7lsgO9ygGAFgpq89E8CVwR1UdYuqZvnaPwHYB6BRmMcNW6Wge9gn8+N8OWlytbdX7PS3x15ssJKKCTcomqjqbl97D4BS7/eISDcAyQC+D/O4EfHQrzpY7TeX/2CwEkp0D3z0rdVu3zT+58gM5hgUIrJIRL4N8WeEvZ8WPRJX4gqKItIMwOsAblLVkL++RWSsiGSKSGZOTvSHV99w8dlW+9GPN0T9eORN9ovlvc8xfjJdIZWdOqjqgJJeE5G9ItJMVXf7gmBfCf3qAJgDYKKqLivlWDMAzACAtLS0qC/bmhy03uPJ/MJi+4jCdct/Mq32S9d3NVhJxYX7UzEbwI2+9o0A/je4g4gkA/gQwH9U9f0wjxdx9lGaD3y03mAllIhUFV9s3W9t16zq+Ls5LoUbFFMBDBSRLAADfNsQkTQRecXX52oAvQCMEpE1vj9xM8h9WGf/VOmzMrMNVkKJaL5tJOZro9JK6Rnfwoo3VT0AoH+I/ZkAxvjabwB4I5zjRFNSkqB+zWQcPH4SALAu+2d0bmn87i0liFvf+MZq97XdkncbviEHMO+unlZ7+AtfGqyEEon9SdGbLktx9RQHDAoAjWtXC9jmkG6KhJ5Pfm61HxjaoZSe8Y9B4fPKSP/7R/t/MFFF7DvqX9u2W6v6xQb4uQ2Dwqf/uf73j/uP5XGkJoWl25QMq/3mmO4GK4kMBoWPiGDSiI7W9u/+WeJwD6JS7TzgvzbRsXmdgBXq3Mr9/4II+r1tpGbmD4dwqoBnFVR+vZ7yv3X94LZLDVYSOQwKGxHBXwa3t7aHPf+FwWrIjdbu+tlqDz+/ecJMksSgCHJr79ZWe9Oeo8g9xTsgVHYjXvTfXn/umrgZVxg2BkUQEcHT/3O+td3+wXkGqyE3sS/s89w1XZDk8jsddgyKEK66sGXA9vb9xw1VQm6Re6oAz3+21dq+8oIWBquJPAZFCRbe3ctq9316sblCyBXsZ55L7+9rsJLoYFCUILVJ7YDtlxfHxVw7FIeWbPHPndK6YU2cWb+GwWqig0FRio2TBlvtafM2IS+fFzYpUEGhYuRrK6zthff0NlhN9DAoSlE9uRLSh/hvl7Z7gBc2KVCbCXOt9id39HD9UO2SMCgc3Nq7TcD2uyt3ltCTvGbWyl1Wu3XDmujUoq7BaqKLQVEGax++3Gr/5YP1OJJ7ymA1FA9+PnES93+wztrOuDcx33KcxqAog7rVqwQ8B9L5kQUGq6F40GXSQqv9ZXo/V881URYMijIaeUlK4LbtAhZ5S0r6HKv98LAOaFHPXat+VQSDohw2P+a/C7JkSw6+tE2aSt7wwmdZAds3XdbKUCWxxaAoh6qVK2GR7fbX9a8st+bapMS3ZtfPeHqBf53a7U9cYbCa2GJQlFPbxrXwoG2Fsa6TFwasK0mJ6cCxPFxpe+Dru0cHJfx1CTsGRQWM7tEKTepUtbbt99Ip8ZwqKMSFjy2ytjPu7e3a9TkqikFRQcvGB65SYL/ARYlDVZE68VNr+7VRaWjTqJbBisxgUFSQiGDT5MEB+66bwenzEomqotV4/9ni+CHt0a99qetwJywGRRiqVamE9Y/4B2N9ve0Axv+XyxImguCQuCbtTPwhaJSulzAowlS7WhWsnOhfx/ntFTvx5LxNBiuicAWHxKVtGmDaVZ0NVmQegyICGtWuisXj+ljbLy3+HtMYFq4UHBKdW9bFW7dcbLCi+MCgiJCUhjWxwDbZzcuLv8e499YarIjKK7+gMCAk+rVvjNm39zBYUfxgUETQOU1q4//u62Ntv78qG4OfW2KuICqzo7mn0NZ2d2PkJWfjtVEXGawovjAoIuzsBjWxYoL/1ummPUd56zTObd13FOfZHvR7ZFgHTBrRyWBF8YdBEQWN61QLeDQdKBpnwRGc8efdlTsxYLr/rO+tW7pjlEee3ygPBkWU1K1epdg4izYT5mLfkdwSPoJirfMj8/GXD/y3s5eN749L2zQ0WFH8YlBEUbUqlbBj6lCM6eH/DdXt8Qx8sCrbYFWUl1+AlPQ5OJKbb+3LmjIETetWM1hVfGNQxMADv+qAWX+4xNq+9721aD2e1y1MWJqVEzD3ad3qVbBj6tCEWEg4msL66ohIfRFZKCJZvr/PKKVvHRHJFpEXwjmmW3VrVT/gukWhFl23OHAsz2BV3pKSPge/f9U/4dBL13ctdi2JQgs3RtMBZKhqKoAM33ZJJgPw9L3CutWrYPsTV6C57RT3wscW4f73Od4imrL2Fr/ztGHSIFxxXjNDFblPuEExAsBMX3smgCtDdRKRCwE0AeD5ySZFBF+N7493x/pH+83KzEZK+hwcPsFJeyNJVZGSPgcDn/X/fhp2fnPsmDoUNZK99Zh4uMINiiaqutvX3oOiMAggIkkAngEwLsxjJZTurRtgy2NDAvadP2kBBj27BKq8jRquD1dnB4yyBIDMBwbg+esuMFSRuznGqogsAtA0xEsT7RuqqiIS6jv8jwDmqmq204xAIjIWwFgAOOuss5xKc73kyknYMXUovvp+P373z+UAgM17j6LV+Ll4+fquGMJT43LbczgXFz+REbBv1KUpeGR4xxI+gspCwvntJSKbAfRR1d0i0gzAYlVtF9TnTQA9ARQCqAUgGcBLqlra9QykpaVpZmZmhWtzG1VFt8czkHM08OLm0vv7JuRalpGWe6ogYKHg07KmDOEdjVKIyCpVTXPsF2ZQPAXggKpOFZF0APVV9f5S+o8CkKaqtzt9bq8FxWkHjuUFTLt22sqJA9CodtUQH+FtpwoKA2agOm3xuD5IaVjTQEXuUtagCPeKzlQAs0RkNIAfAFztO3gagFtVdUyYn99zGtSqih1Th2Lj7iMY8tel1v6LphSFB88wihzLy0enh+cX2//aqDTPzkIVTWGdUUSTV88ogi3Nygm493/aP0emYWAH7/1ArM8+jGEvfFFs/9+uuwDDz29uoCJ3i8lbj2hiUARasf0grv7H18X2165aGasfGojKCfw+XFUx6Lkl2LL3WLHXXh/dDT1TGxmoKjEwKBLU3iO56P54RsjXHhh6Lsb0bB3jiqLnreU7MeHD0HOQrpjQH43r8NmMcDEoEtypgkLc9K+V+KKEZQ3v6NcW9ww8x3WL1DyzYDOe/2xryNcuad0Ab93S3XX/pnjGoPCQfUdy0a2Es4zTFt7dC6lNaseoorLbvv84+j69uNQ+ax4aiHo1kmNTkMcwKDyqLD94ADC2V2vcN6hdTMcY5BcUYuKH3+LdzF2OfVc/OBBn1GQ4RBuDgnCqoBAPfvQt3lnp/INpN2lER/y2a8sKLZt34mQ+3lq+E4/N2Viuj7tvUDv8sU8bvq2IMQYFFVNQqJiVuSsuFil6bVQa+rZrzGAwLFYDrshFKiUJrut2Fq7rFvgczcHjJzF94Wa8sWxnxI9536B2GN2jFapVqRTxz02xwzMKIg8r6xlF4o7SIaKIYVAQkSMGBRE5YlAQkSMGBRE5YlAQkSMGBRE5YlAQkaO4HXAlIjkoml4v0hoCCP1sdnxyU71uqhVwV73RqvVsVXWc+SdugyJaRCSzLCPR4oWb6nVTrYC76jVdK996EJEjBgUROfJiUMwwXUA5ualeN9UKuKteo7V67hoFEZWfF88oiKicEj4oRKS+iCwUkSzf32eU0reOiGSLyAuxrDGoBsd6RaSLiHwtIt+JyDoRuSbGNQ4Wkc0istW3lGTw61VF5F3f68tFJCWW9QXV4lTrPSKywfd1zBCRs03Uaaun1Hpt/X4rIupblS/qEj4oAKQDyFDVVAAZvu2STAawJCZVlaws9Z4AMFJVOwIYDOA5EakXi+JEpBKAFwEMAdABwHUi0iGo22gAh1S1LYBnAUyLRW3ByljrahSth9sZwPsAnoxtlX5lrBciUhvAnQCWx6o2LwTFCAAzfe2ZAK4M1UlELgTQBMCCGNVVEsd6VXWLqmb52j8B2AcgVstldQOwVVW3qepJAO+gqGY7+7/hfQD9xczkmI61qurnqnrCt7kMQMsY12hXlq8tUPQLbRqA3FgV5oWgaKKqu33tPSgKgwAikgTgGQDjYllYCRzrtRORbgCSAXwf7cJ8WgCwT+ud7dsXso+q5gM4DKBBTKoroQ6fULXajQZQfGn02HGsV0S6AjhTVefEsrCEmFxXRBYBaBripYn2DVVVEQl1m+ePAOaqanYsfvFFoN7Tn6cZgNcB3KiqhZGt0ltE5AYAaQB6m66lJL5faNMBjIr1sRMiKFR1QEmvicheEWmmqrt9P1j7QnS7BEBPEfkjgFoAkkXkmKqWdj3DZL0QkToA5gCYqKrLolFnCX4EcKZtu6VvX6g+2SJSGUBdAAdiU17IOk4LVStEZACKQrq3qubFqLZQnOqtDaATgMW+X2hNAcwWkeGqGt2ZqFU1of8AeApAuq+dDuBJh/6jALwQz/Wi6K1GBoC7DNRXGcA2AK18dawF0DGoz58A/N3XvhbALENfy7LUegGK3ralmvo/L0+9Qf0Xo+hCbPRrM/3FicEXv4HvhyoLwCIA9X370wC8EqK/6aBwrBfADQBOAVhj+9MlhjVeAWCL7wdsom/fJADDfe1qAN4DsBXACgCtDX49nWpdBGCv7es42/D3a6n1BvWNWVBwZCYROfLCXQ8iChODgogcMSiIyBGDgogcMSiIyBGDgogcMSiIyBGDgogc/T/4x80IQOC5ewAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fix, ax = plt.subplots(1,1)\n",
"ax.set_aspect(\"equal\")\n",
"ax.set_xlim([-.5,.5])\n",
"ax.set_ylim([-.5,.5])\n",
"ax.plot(orbits[:,6::60],orbits[:,7::60])"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [],
"source": [
"def leapfrog_step(y,dt):\n",
" yc = y.copy()\n",
" # Advance positions for half a timestep\n",
" yc[0::6] += yc[3::6]*dt/2.\n",
" yc[1::6] += yc[4::6]*dt/2.\n",
" yc[2::6] += yc[5::6]*dt/2.\n",
" # Advance the velocities for a full timestep\n",
" rhs = G(yc) # updated values!\n",
" yc[3::6] += rhs[3::6]*dt\n",
" yc[4::6] += rhs[4::6]*dt\n",
" yc[5::6] += rhs[5::6]*dt\n",
" # Advance positions for half a timestep\n",
" yc[0::6] += yc[3::6]*dt/2.\n",
" yc[1::6] += yc[4::6]*dt/2.\n",
" yc[2::6] += yc[5::6]*dt/2.\n",
" return yc"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [],
"source": [
"steps = 100\n",
"dt = 0.1\n",
"y = y0\n",
"orbits = np.zeros((steps,Nd))\n",
"for i in range(steps):\n",
" y = leapfrog_step(y,dt)\n",
" orbits[i] = y"
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x10ce5ec50>]"
]
},
"execution_count": 73,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAD8CAYAAACPd+p5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4VHX2h9+bTHrvvUFIQhJCSYCABRWkVymK3bXturv6s5e1ra5r3dV11bVhL1QVFFBBAaWTBFII6b33mbTJtPv7404mMyEUJQQm3Pd5eHJbZu6Qmc+cc76nCKIoIiMjI3MybM71DcjIyJz/yEIhIyNzSmShkJGROSWyUMjIyJwSWShkZGROiSwUMjIyp0QWChkZmVMiC4WMjMwpkYVCRkbmlCjO9Q2cCF9fXzEyMvJc34aMzLAmPT29SRRFv1Ndd94KRWRkJGlpaef6NmRkhjWCIJSfznWy6yEjI3NKZKGQkZE5JbJQyMjInBJZKGRkZE6JLBQyMjKnRBYKGRmZUyILhYyMzCmRhUJGRuaUyEIhIyNzSmShkJGROSWyUMjIyJwSWSgucORxDTKnw3lbFCYzNGzLreeJjTkkBnuQEOJBQrA7iSEeBHs4IgjCub49mfMEWSgucHxcHZg60pejNUp25DdgMBoYns52RvFwJyHYg8RgdyJ9XLCxkcXjQkQWiguc5AgvkiO8AOjW6DlWp+JojYqj1UpyapR8uLsMjd4AgIu9LfHBknD0Wh7R/q7Y2coe7HBHFgoZE072tkwI92JCuJfpmEZnoLCh3Uw8VKxNq6RLowfAXmFDXKCbhXjEBbrhaGd7rl6GzFlAGIxgliAIs4H/ALbA+6IovnCC65YC64GJoiietCtNSkqKKDeuGVwMBhFB4IxjD3qDSGlTJ0drlBytUZFTLf1UdmsBsLURiPZzJSHEXXJfgt2JD3bHzdFuMF6GzCAiCEK6KIopp7zuTIVCEARboAC4EqgCDgErRVHM7XedG7AZsAf+IgvF2UfZpSWjspWM8lbSylrJrGrDIIpE+rgQ5etCpK8LUT4uRPm5EOnjgq+r/e8WEVEUqWrtthCPnBoVje09pmsifZxJCPEwiUdCsDs+rg6D9XJlfgenKxSD4XpMAopEUSwxPvFqYBGQ2++6Z4EXgQcH4Tll+iGK0rd8WrkkDOnlrRQ2dADSN3x8kDvLk0NR2NpQ1tRJfl0723Lr0Rn6vihcHRRE+joT5etKlI8zkUYxGeHrgqez/UmfXxAEwrydCfN2ZnZikOl4g0ptYXVkVraxOavWdH6EnwtPzI/n8lj/Qf4fkRlMBkMoQoBKs/0qYLL5BYIgTADCRFHcLAiCLBSDgFqrJ7OyjfSKPmFo7ZJMfw8nOyaEe7JoXDDJEd6MDfPA2f74P7VOb6C6rZvSpk5Kmzopa+qktLnL+GGuwUxD8HS2M1ki5tZIpK/zSV0Kf3dH/N0duTyuTwjaujTk1qjIqVGyNq2KWz48xOJxwTwxP162MM5TznowUxAEG+DfwM2nce0dwB0A4eHhZ/fGrIw6pZp0oyCkV7RytFppsgZG+LkwY3QAKZHSCsYIX9fTWsZU2NoQ4eNChI8Ll8VantPoDFS0dFHW1ElZcyclRiE5UNLM14erLa71dXUgyteZSJ8+CyTSV3JnnOyPD2p6OtszNdqXqdG+3DQ1krd2FPPWziJ+KWziyfnxLBoXLOdwnGcMRoxiCvC0KIqzjPuPAoii+Lxx3wMoBjqMvxIItAALTxanuJBjFDq9gby69j5hKG+luq0bAAeFDWPDPEmO8CIlwovx4V54u5zcLRhsujV6yluMFkhTl9ESkawS85gEQJCHo0lAonydGeXvxqUxftj2E7L8unYe3pDFkco2Lov147klYwjxdBrKl3VBMpTBTAVSMHM6UI0UzLxWFMWjJ7h+J/CAHMzswzzomF7eypHKNtPyY4C7AykR3qZ8h9FB7tgrTp63UKdU80tBIzsLGsisVOLhZEeghyMB7o4EeTgS6O5IgIe0HeDuiLujYtC+wTt6dCYrpLRREhBpv4uWTg0ACcHuPLMo0ZS/0YveIPLx3jJe/iEfGwEemh3HDakRcpLXWWTIhML4ZHOB15CWRz8QRfE5QRCeAdJEUdzU79qdXOBC0djew878BjIqJGEoqO8LOo4OciMlwpsJRmE4nVRqjc5AWnkLuwoa2ZXfSF5dOyCJzMRIb7o0euqUaupUatOH1RwnO1sCjQJiLigBxv1Ad0f83ByOswJ+K8ouLTsLGnh+Sx51KjUrUkJ5eHbccXGJypYuHvs6m18Lm5gQ7smLS5MYFeB2Rs8tMzBDKhRng+EoFAaDyBcHK3hhax4dPTrcHRVMMLoQEyK8GBvqiYvD6YWNqlq72JnfyK6CRvYWNdGp0WNnK5AS4c20WD8ui/UjNsDtOJHp0elpUPVQaxSOeuPPOrOfDe1qtHrL94WtjYCfq4NkiZgJSqCHA4HuTiZBGSgm0Z/OHh2v/1zIql9Lcba35cHZcVw7KdxCiERR5OvD1TzzXS5dPXr+fHk0f7ps5CmtKZnfhiwU5xkljR088lU2B0tbuCjah8fmjmZ0oPtpm9VqrZ6DpZLVsDO/geLGTgBCPJ24LNaPaTF+TI32xfU0heZkGAwizZ0a6o3CUTuAoNQr1bT36I77XQ8nO5NrE+juQKCHE4FGC2VilLfF/RU1tPPEN0fZV9JMYog7zy5KZHy4pTvS1NHDM9/msimzhpgAV15YmmSROSpzZshCcZ6g1Rt479cSXtteiKPChsfnx7M8OfS0YgKlTZ3sym9gV0Ej+0qaUWsN2CtsmBzlzWWx/kyL8WOkn8s5WyHo6NFRp1SbBKW/ZVKnUtPU0UPvW8zX1YGHZseybEKoSSBFUeS7rFr+sTmXelUPV6eE8dDs2OPckZ+O1fP4NznUqdTcPDWSB2bGnrb1JXNiZKE4D8ipVvLQ+ixya1XMSQzk74sS8HdzPOH1XRod+0uaTS5FeXMXAFG+LkyLkayG1BE+JzXvSxo72JJdy+bsOoobO/BytsPL2R5vl75/J9r3crHDQTG4NRpavYHG9h6KGzt4dVsBGRVtjAnx4KkF8aREepuu6+jR8fpPhXywuxQXBwUPzoplZT93pF2t5aXv8/l0fzkhnk7886oxTIs55XxdmZMgC8U5RK3V89r2Qt77tQRvF3ueXZTI7MTA464TRZGihg6jO9HIwdIWNHoDTna2TB3pwzSjSxHh43LS5ytqkMRhS3atKZA5IdyTCeFetKt1NHdqaO3S0NqpoaVLQ5sxMWsgXB0UeLnY4e1sj1evmJhv9xMXTye703afRFFkU2aNKZi5cGwwj8yJI9hsGbSgvp0nN+awv6SFpFAPnlmUyLgwT4vHOVTWwsMbsihp7OSq8SE8MT8eryFeIh4uyEJxjthf0syjX2VT2tTJ1SlhPDZ3NB7OfZmL7Wote4qa2VXQyC8Fjab8iFH+rkyL8eOyWH9SIr1OWX1Z1NDO5qw6tmTXkl8viUNKhBdzxwQxZ0wgQR4nzkHQ6Q20dWsl4TCKSHOnUUg6tf32pfO9y7X9sRGkBCovZzuTgPi4mlkpzvZ4u0piE+7tjJeLPV0aHW/vLOadX0oQBPjjtJHceelIk6XUKyjPbT5GY0cP10wM48FZcRb5Imqtnjd3FPG/ncV4ONnx1MIEFiQFyYlavxFZKIYYlVrL81vy+PJgBeHezrxw1RimRvsC0jf+ttx6duY3kF7eis4g4uqg4KJoH6bF+DMt1u+0kosK69vZbLQcCuo7EAQzcUgMItDjxG7NmaLW6mkxE44WMyFp6dLQ2qk1nW8xWi/mdSQglaTfenEUd102EjdHO6pau3h+ax6bs2oJ9nDkkbmjLT7s7Wot/9leyId7y3BzVPDQrDiumRhmYcEcq1XxyIYsMquUTI/z59nFiRYWiszJkYViCNmWW8/j32TT2N7DbZeM4N4ZMTjZ26I3iLy5o4jXthdgEGF0kLtphSI5wuu0Gr4U1LezOUsSh8IGSRwmRngzd0wgc8YEEeB+vDiIosjhyja+OVxNvUqNh5Md7o520k8nO9ydFKZj7k52pm1HO5tB+0YWRZH2Hh2tnZJ10tKhYUtOLV9lVOPr6sADM2NYnhKGrY3AgZJmnvkul6M1KlIivHhqQQJjQj1Mj5Vf184TG3M4WNrCWKM7MtbMHdEbRD7cU8q/fizA1kbg4dmxXDdZTtQ6HWShGAIa23t4+tujbM6qJS7QjZeWJZEUKr2B65Rq/m/NYfaXtLB4XDCPzh094Ie6P6IoUlDfYbIcinrFIdKbeWOCmJ0YeMLHqWnr5uvD1WxIr6KkqRNHOxvCvJxpV+tQqbUndB96sbe1wd1JgbujHW4mAVH0CYxJbBSWwuOowN3J7rSE70hlG89+l0t6eSujg9x5cn48U0b6oDeIrEur5OUf8mnp0rAiOYwHZsXi5+Zg+n/ZeKSG57Yco6mjh5WTwnlwZqxFbMI8USslwosXliYR7e96ynu6kJGF4iwiiiJfZVTz7GYpGeju6dHcOW2k6YPyc149D6zLoluj59nFiSydEHLSb2pRFMmvb2dLVi2bs2spbuxEEGBSpDfzkoKYnRCI/wnEoUuj4/ucOjZkVLG3uBlRhElR3iybEMqcMYEWlZ0anYF2tRaVWoeyW4uqWyv9VGtRdevMtnuP62g3biu7tce5Ev1xtrc1s1QsrRZzQfF2tqe1S8Nr2wupbutmVkIAj80dTYSPCyq1ljd+LuLDPaU4KGz56xXR3HxRpGk1RqXW8tq2Qj7eV4a7o4KHZ8exIiXMYrl1Q0Y1z36XS7dGz1+vkP42cqLWwMhCcZY42beWRmfgpe/zeH93KaOD3Hnj2vGM9Bv4G00URfLq2o1LmbWUNHZiI8DkKB/mJgUxKyHghEupBoPIwbIWNqRXsSW7lk6NnjBvJ5ZOCOWq8aGE+zhbXN/Zo+PXwkZAMH1Yez/Mrg4KFKdhCYiiiFprGEBMzETGuC9tWwpPe48O87daiKcT98+Moaatm7d2FqPVG7jloij+ckU07o52lDZ18tzmXLYfayDCx5nH58UzY7S/SXCP1ap4cmMOh8paGRfmybOLEi3clf7W3otLkyzcFRkJWSgGmd6CpVd+zEcAHpkTZ+EHlzV18tcvD5NdreTmqZE8MifuuJULURQ5VttuWsosaZLEIXWED3PHBDErIdBkag9EeXMnGzKq+SqjiqrWblwdFMwdE8jSCaFMjPQ+zic/WqPky4MVfHO4ho4Bsih7cbG3xd3JDjdHRZ8F4KjAzbHPzTjZ+dPJvTAYpJiFqltLUUMHL36fR15dO5OjvPnTZSPZnFXL+owqvJ3tuW9mDFenhKGwtWFXQSPPfpdLUUMHl4zy5Yn58cQY6z5607z/uSWP5s4erpsczgMzYy2a7JjHj265KIr7Z8YM2JvjQkUWikGkoF4qgT5c0cblsX78o18J9MYj1Tz2VTYKWxteXpbEzIS+nAlRFDlao2JLdi1bc+ooNYrDlJF94uB7kmYtKrWWLVm1bMio4lBZK4IAF0f7snRCKLMSAo9Lvurs0fFtZg1fHqwgs0qJvcKG+WOCWJYSioeTnRSv6NZS3txFvUqNjY2ATi/SrtaaYhm9VkKvm6I/hcvhoLA5iagYf5pZMm6OClwcFGzNruWdX0rQ6A1cMzGcWQkBvLWjmINlLcQFuvH4vHguHuWLVm/g8/3l/HtbAZ0aPddPDufeK2NMgqBSa3l1WwEf7y3D09meh2fHsjy5zx1RqbW8uDWPzw9UEObtxGtXjyM5wvtkL+mCQRaKQUCjM/DWziLe3FGEq4OCpxcmsHBsX1OVLo2OpzYeZV16FRMjvXjtmvEmAVGptbyzq5jNWbWUNXdhayMwxWQ5BJy0k5PeILKnqIkNGVV8n1NHj87ACD8Xlk4IZcn4kAGX/3KqJeth4xHJehjl78q1k8NZMj7E9IEyGET2Fjfz6f4yth9rMAmAwkbA01kKTno6S0lUHsZ9Dyc77BU22AoCNoKAIICNINBrvHRq9CZhMbkZvYJjFJseneG0/88fnzcaPzcHXvkxn8qWbmaM9uexuaMZ4edKS6eGV7cV8PmBctwc7bjvyhiumxxucp1yayR3JK28lfHhkjuSGNLnjhwoaeaB9Zl0a/T8eO+0Ie/jcT4iC8UZcriilYc3ZFFQ38GiccE82a9NW26Nir9+mUFJUyd/vTyau6ePMr1hS5s6ufXjQ5Q1dXJRtK/JcjjVG7OooZ316dV8c7iaOuOy5oKxQSydEMq4MM/jAqL9rQcHhQ3zkoK4dlI4yRFepuuVXVrWpVfyxYEKSpo68XaxZ3lKKCN9XWnrljI127q1KLu0tHVrUHZraeuS9gcq/DLH3VEhiYuZsHg62+Hp1HfMyd4WW0FAMAqMjSCg1ulN1k1JYydr0vq6KUb6OPPInNGUNnXy5o4i1Fo9N02N5O4rRuHhbEdenYpnvs1lb3EzMQGuPDE/nktGSancBoPIV4ereWHrMVo6NVyfGsH9V8aakt6O1apY8N/dzB0TxOsrx5/mu2H4IgvF70QURV7Ymse7v5YQ5O7IP5YkckVcgMX5z/aX8+zmY3g62fHaNeOYOtLXdP7Xwkb+/HkGtjYC/7s+mdQRPid9vrYuDd9m1rA+vYrMKiW2NgKXxfixNDmU6aP9B/T/c6qVfHGwgo2Hq+nU6IkJcOXaSeEsGR9qkQWaXaXk0/1lbMqsQa01kBzhxQ2pEcwZE3jaNR1avQFVtyQkbV1a47ZRXLq0RlExiotJbKRjJ/NYnOxszawYO+pVPZQ2dVpcM8rflWcXJ7LxSDWrD1Xi6WTHvVfGmErSf8yt57nNx6ho6WLG6AAenzeaSF8p3V3ZLbkjn+wrw8vZnofnxJmK0f6zvZBXtxfw9vXJA6bWX0jIQvE7eXtXMS9szWPlJCn92nx5Udml5aENmfxwtJ7LYv341/KxJitDFEU+2lvGPzYfI9rPlfdvSiHM23nA59DqDezKb2RDRhU/HWtAozcQF+jGsuRQFo0LGTCg2dmjY5PResgysx6umxzOhPA+60Gt1fNdVi2f7i8ns7INJztbFo8P4frUcBKCPY573LOFwSDSodGhNIlJn8D0Lre2dVlaM61dGhr6tdIDuHFKBHPHBPGf7YXsK2lmlL8rj8+PZ1qMHz06PR/sLuONnwvR6A384eIo/nJ5tOnvdrRGyRPf5JBR0cbM+ADeuSEZnUFk0Rt7aGjvYdu9l17QdSKyUPwO9hU3c937+5kzJog3Vo63MPXTylq4Z/URGtrVPDw7jj9cFGUKlml0Bp7cmMPqQ5VcGR/Aq1ePG7AvxNEaJRvSq9mUWU1ThwYfF3sWjQthaXLICT/Ep2s9lDV18vmBctalV9HWpSXa35XrJ4dzVXIo7lY2eKepo4dHNmSx/ViDxfHnliTi42LP81vzKG/u4rJYPx6fN5pofzcaVGpe/iGfdelVUjn7rFiWJUsWhMEg8tbOIl75sYDnliRy3eQIjtYoWfTGHuYnBfHaNReuCyILxW+kXqVm3uu7cXdSsOkvF5s+6HqDyNu7ivn3tgJCPJ3478rxFuvxzR09/OmzDA6WtfDny0dy/5WxFsuUje09bDxSzfr0KvLq2rG3tWH6aH+WTghlWqzfgNmMHcbYwxcHKsiulqyH+UnBXDs5zMJ60BtEfs5r4NP95fxS0IjCRmBWQiDXp0aQOsLb6gukihramfHvXyyOBXk48uLSJPLr2nn9p0K6tHpuSI3gnumj8HKxJ6uqjb9/K2V+Joa489SCBCZGemMwiFy/6gBHKtv44f8uJczbmVe3FfCfnwp594Zki5WqCwlZKH4DWr2Ba9/bz9EaFRv/fJGpP2ODSs29a4+wp6iZBWOD+eeSRAtX5Fitits+TqOpo4eXliWxaFyI6dyRyjb++1MhOwsa0RtExoZ6sDQ5lAVJwSc0dXOqlXx+oIJNRyTrITbAjWsnh7N4XIiF9dDY3sOaQxV8caCCGqWaQHdHVk4K55pJYaeVJm5N9GZaPrAu0+L4nMRA/jhtJGvTKvnyYAVujnbcM30UN0yJQGEjsCmzhhe25lGrVLPAWM4uiiKzX/uVhGB3vrw9FZ1BZOEbu2nu1LDt3ktPOeRoOCILxW/gH9/l8v7uUv5zzTjTh31nfgP3r82kU6PjmYWJLE+x7Er1fU4d9609gpujgndvSLGwMrbl1vOXLzJwd7Jj6YRQliWHEO0/cHPYjh4dm45IsYfsaiWOdpL1sHJSOBPC+1Y6RFHkUFkrn+4v5/ucWrR6kYujfbk+NYIZo/1PK7vSmunS6Lj7y8PHuSN/uTyaK0b78+8fC9hd1MQIPxcenzeay2P96dbqeXtXCe/sKkYQ4NE5o3G0s+HhDdk8tSCeWy6KIqdayeI397BwbDD/vnrcOXp15w5ZKE6TLdm13PV5BjdNieDvixLR6Ay88mM+7/5SQlygG29cO97iQy6KIm/8XMS/thUwNsyTd29ItvgWX5tWyaNfZZMY7M6Ht0w64ZJodpUUezjOehgfgoeTZf+Kbw5X8+n+cgrqO3BzVLA8OYzrUsNPmB4+nClu7GD6v3ZZHPNwsuOpBfG4Odrx/JZjlDR1cskoXx6fF09soBtVrV08/k0OuwoaWXvnFN7aUcS+kma23H0JI/xc+fe2Al7/qZD3b0xhRnzACZ55eCILxWlQ3NjBojf2EO3vyto7p1CnVPPX1YfJrGzjhtQI/jZvtEUadrdGz4PrM/kuq5bF44J5YWmS6bwoirzzSwkvbM3jklG+vH198nE9HU9kPVw7OZzx/fIk8upUfLa/nK8zJCFJDHHnxtRIFowNPq1O18MZURRZm1bJwxuyLY6PD/fkb3NHk1Wl5LXtBXT06Lh2cjj3zojBwc6W2a/9gsJG4KNbJrHwjd1E+7uy7o9T0RtdkJZODdvunWbh5g13ZKE4BV0aHYvf3ENTh4bv/nox6eWtPPZVNoIALy1Lshi0C1Cr7OaOT9LJqVHy8Ow47rx0hOmDbTCI/HPLMd7fXcqCscH8a/lYi2pFyXooZ+ORGro0euICJeth0ThL60GjM7A1p5bP9pdzqKwVe4UNC5KCuWFKBGNDPaw+ODnYqLV6bv7wIPtLWiyOXzUhhDsuHcGXByr47EAFzva23DN9FKMC3Ljpg4PcclEkSaEe3Lsmk0fmxPHHaSPJqVay6M09LB4Xwr9WjD1Hr2jokYXiJIiiyH1rM/nmSDUf3TKJH47W8cWBCiaEe/L6yvGEelnmP2RUtHLnp+l09eh4feV4po/uM0+1egMPrc/i68PV3Dw1kifnx5tWPSqau/jrlxlkVknWw4KkYFYOYD1UtXbx5cEK1hyqpKlDQ4SPM9dPjmBZcugFvcZ/uhQ1dDDj35buiIPChrunj2JajB8v/5DProJGIn2c0egM1CjVrL4jlQ/3lLIjr5Hv7r6YmAA3/vVjPv/9uYgPbk6xSLIbzshCcRI+3V/OE9/kcP+VMTg7KHj2u1zuvHQED8yKPW658quMKh75KptAd0fevynFVLkIklVy1+cZ7Mxv5IGZMfz58miTAJQ3d7Ly3f10afXcd2UMi8eHWOQzGAwivxQ28tn+cn7OkwJ0V8QFcMOUCC6J9pW7M/0OPt5bxlObLCdZhnk78be58TjY2fDc5mMUNUhT2cK9nfns1sksfmsPIZ5OfHXXVAyiyML/7qGtW8OP906zsPaGK7JQnIAjlW0sf3svF0f78rd58cx7/Vcujvbl/ZtSLL7l9QaRl77P451fSpgywoe3rptg8e3e1qXhlo8OkVnZxnNLxrByUt/09bKmTla+tx+1Vs/nt6USH+xuOtfaqWFdeiWfH6igvLkLX1d7rpkYzsrJ4fJQ3kFArdWz5K29HKtVWRyfOlIauvTpvnJTXcmNUyJIHeHDXZ9ncO+MGO6ZMYqsqjaWvLWXq8aH8PLy4e+CnK5QXFCF+S2dGu76LB1/N0deWT6W2z5Jw9HOluevGmMhEu1qLXd/eZgd+Y3ckBrBkwviLSyNWmU3N646SHlLF29dN8EinlHaJFkSGr2BL25PZXSQJBLKLi3PbpYmXml0BiZFenP/zFhmJwTK3ZcGEUc7W7becwn5de3Meq0vWWtvcTML39jN4nEhuNjb0qnR88m+cmYnBLJgbDD//bmQ6aP9SQr15I/TRvDmjmLmJgVxeaz/OXw15w8XjEWhN4jc/OFBDpS0sOFPU9lT3MQLW/MscidAsgZu+ySNsqZOnlqYwA2pERaPU9TQwY2rDtCu1vHujSlMGdlX9FXS2MHK9/aj1Yt8cftk4gIlkejs0XH9qgMcrVZx9cQwrk+NIDZQHro7FLy1s4iXvs8/4fkQTydW35HKVf/bi4+LPRv/chEA81/fTbtax4/3XWp1KfC/hdO1KC6Yr7LXfyrk18Imnl6YgIOdDf/+sYDZCYEsHBtsumZPUROL3txDU0cPn9w66TiROFzRyvK396LRi6y+M9VCJIobO7jm3f3o9CJf3p5qEgm1Vs/tn6SRVaXkv9eO59nFibJIDCF3XRZN3rOzCT7BKIPqtm7+t6uY55eMIc+YFu6gsOWV5WNpaFfz3HfHhviOz08uCKHYmd/A6z8XsnRCKMtTQrl/bSaujgr+sSQRQRAQRZFP9pVx4wcHCXB3YNOfL7YoHQfYVdDIte8dwM3Rjg1/mmJRxFXU0MHKd/djEEW+vCPVJARavYG/fHGYvcXNvLI8iVkXaD3BucbRzpa9j05n458vGvD8FwcqsFfYsCw5lP/tLOZwRStjwzy5c9pI1qRVsqugcYjv+Pxj2AtFVWsX/7fmCLEBbvxjcSJv7ywmu1rJPxYn4uvqgEZn4G/f5PDkxqNcHuvHhj9NPa457cYj1dz60SGifF1Y/6cpFiP+ihraWfnefgwifHl7qmlVxGAQeWBdJtuP1fPsogSWjA8d0tctczxjwzwpe2Eef5w28rhzj2zI4t4rYwhwd+T+dZmotXrumT6KaH9XHtmQhUp94jGMFwLDWih6dHru+jwDvV7k7euTKW3q5PWfC1kwNpi5Y4Jo6dRww6oDfHGggj9dNpJ3bkixKPoC+HBPKfesPkJyhBd1v7t6AAAgAElEQVSr70y16IxdWN/ONe8eQBRh9R2TTcVkoijy+MYcNh6p4aHZsdwwJXIoX7bMKXhkThzZT8/E2SzDtUap5p9bjvHi0iRKGjt55Yd8HO1seXlZEvUqNc9vubBdkGEtFM98m0tWlZJXVowl2NOJ+9dl4uFkzzMLE2hoV7Pozd0crmzjtavH8fDsOIvJ2aIo8vIPefz921xmJQTw8R8mWQS1CuolS0IQYPUdqaZ6kN4OWV8cqOCuy0Zy12XRQ/66ZU6Nm6Mduc/MtmiHtzmrln9vK+C6yeGs2lPKwdIWxod7cfulI/jyYCW/XMAuyLAViq8yqvj8QAV3ThvBrIRA3vi5kGO1Kp6/agxeLva8sDWPemUPq+9IZfH4EIvf1ekNPPpVNm/uKGblpDDeui7ZouYjv66dle/ux0YQjCLRV5z15o4i3vmlhBunRPDgrNghe70yv4+FY4O5/ZIo0/6Ryja25dbj42LPg+sz6dLouHdGDCP9XHj0q2zaL1AXZFgKRV6dise+zmZylDcPzowlu0rJmzuLuWpCCFfGB5Be3sJXGdXcdkkUE8K9LH5XrZXcldWHKvnrFdH8c8kYC0sjr07Fte/tR2EriYR5BedHe0p55ccCrhofwtMLEuTaDCvhsbmjLfYb2nto6tBQ3tzFC1vzJBdk+Vhqld08vzXvHN3luWVQhEIQhNmCIOQLglAkCMIjA5y/TxCEXEEQsgRB+EkQhIiBHmcwUKm1/OmzDNwd7fjvtePRiyL3rzuCr6s9T81PQG8QeXpTLoHujvz58ujjfvemDw7yY249Ty+I5/6ZsRYf9mO1Kq597wB2tjasvmMKI8xEYl1aJU8b3ZSXliXJKdhWhCAI7HjgsgHPfbKvnLVplUwI9+K2S0bwxYEKdhc2De0NngecsVAIgmALvAnMAeKBlYIgxPe77DCQIopiErAeeOlMn/dEPLw+i4qWLt64dgL+bo68uq2QgvoOXlyahIezHevSKsmuVvLo3DiLMvAGlZqr39lPRkUr/7lmHDdfFGXxuLk1kiVhb2vD6jtSifLtW/nYml3LwxuyuGSUL6+vHD/sm8gMR6J8XZjdb/n6SmNviofWZ1GnVHPflTGM8HXh4Q1ZJ528NhwZjHf0JKBIFMUSURQ1wGpgkfkFoijuEEWxy7i7Hzgra4VFDe1szanj7itGMSnKm4yKVt79pZhrJoZxWaw/yi4tL/2Qz8RIL4tEq/LmTpa9vY/y5k5W3TTRIlMTpKa4172/H0c7W1bfkWpqCQ9Sjsbdqw8zPtyLd25IPu02+DLnH/+7foLF/rbcep6YL33npT7/k9EFSaJG2c0LWy+sVZDBEIoQoNJsv8p47ETcCmwdhOc9js1ZdQgCXDMpDLVWzwPrMgnycOJv8yQf9NXtBbR1aXh6YV/8IKdaydL/7aVdreWL21O5NMbP4jFzqpVc9/4BnAYQiYOlLfzxs3RG+bvxwc0T5ZmWVo4gCHx111SLY/l1fcVlb+8qJjnCm1sviuKz/RXsLbpwXJAhtZEFQbgeSAFePsH5OwRBSBMEIa2x8bcvRW3JrmVihDcB7o68/EM+JY2dvLQsCTdHO/Lr2vl0fznXTu6bb7G3uIlr3t2Pg8KWdX+cyrh+0657RcLFXsHqOywTrbKrlNz60SGCPZ345NZJF0RJ8oXAhHAv4sxS7NemVZnE44WtebR1abh/ZixRvi48tCGLzgvEBRkMoagGwsz2Q43HLBAEYQbwN2ChKIrHT3kBRFF8VxTFFFEUU/z8/Aa65IQUNXSQX9/O3DGBHCxt4YM9pdyQGsFF0b6IosjTm47i6qDg/iulJcvdhU3c/MEhgjwcWf+nKRZLnCAJwbXv7cfVQcHqO1ItsjUL69u58YMDuDvZ8fltk086ZFjG+th89yUW+1e9tdfURmDcM9twsrflpWVJVLd18+L3F8YqyGAIxSFglCAIUYIg2APXAJvMLxAEYTzwDpJINAzwGGfMluxaBAGmxfrz4PpMwryceWROHABbc+rYV9LMAzNj8HKxl1Y+vj1KqJcT6/44hSAPyz4QWVVtXPf+ftyd7Fh9R6rFxK+K5i6uX3UAha0Nn982+bjflbF+bG0E3r0h2eJYoFkD5XVplUyM9OaWqVF8sq+cfcXNQ32LQ84ZC4UoijrgL8APwDFgrSiKRwVBeEYQhIXGy14GXIF1giAcEQRh0wke7nezJbuWlAgvPtxTSkVLFy8vS8LFQUG3Rs9zm48xOsidaydLq7LfZtZQ1NDB/TNjj5vlkFnZxnXvHxhQJOqUaq5btZ8enYHPbp1sEa+QGV7MTAjE0a7v4/Hq9gKTeDy4Povmjh4enBVLhI8zD2/IQqc//Ynt1sigxChEUdwiimKMKIojRVF8znjsSVEUNxm3Z4iiGCCK4jjjv4Unf8TfRlFDB3l17Xg62/PJvnJumRrFZONw4P/tKqa6rZunF8RjayOg0xv4z0+FxAW6MaffgNrDFa1c//4BPJ3tWHPnFIvemS2dGq5fdYDWTi0f3zJJLhW/AMh44kqL/Ts+TSclQkrQW/72PhztbPjzZdFUtHRR0dI10EMMG4bFgv+W7FoAduQ1EOXrYkqdrmzp4u1dxSwYG2wSjq8PV1Pa1Mm9V8ZYJEVlVLRy46qDeLnYs+aOKRZt6XoTsSpbunj/JsthPzLDF2d7Bc8sSrA4dnmc1PGqpKmTTZk1RPlJVmW5LBTnP71CYRBFXlk+1jT34h+bc7EVBB6bK8UqtHoDr/9cSGKIOzPNBr2kl0si4e1qz5o7Uwk2E4lujZ5bPzrEsVoVb1+fTOoIH2QuHG7sV/n78g/53DsjBoB7Vh/ByVgDVN7UOdS3NqRYvVAUN0puB8Dtl4wg2Wga/lrYyA9H6/nLFdGmgOP69CoqW7q578oYUx5FenkLN31wEF9Xe1bfkWoRnOzR6bnzs3TSy1t57Zpxpm8TmQuLg3+bbrH/xcFy3B2lnJm/f3sUJztb2aI439mcJVkTbo4K7r1SUnqt3sDfv80l3NuZWy+WUrF7dHr++1Mh48I8TQ1TK1u6uHHVQfzcHFh9h+Xqh05v4P9WH+GXgkaev2oM85OCkbkw8Xdz5LaL+1L661U9LE2WkosPlbXSrdVT0SwLxXnN6z8VAvC42fi/j/eWUdTQwZPz403H1hyqpEap5v6ZfdbEqt2laPQGPr11EoFmPRUNBpFHvspma04dT8yP5+qJ4chc2PRm9/by4Z4ykkL72iHKFsV5TEljBzqD1EW8t8dlY3sP/9leyLQYP6aPliwHtVbPmzuKmBjpxcXR0nXKbi3r0ipZkBRssbohiiLPfJfL+vQq/m/GKJNFInNhIwgC2++71OJYW1dfb4qK5i4MhvOzo/1gYNVC0RvEtFfYmPIdXv4hD7VOz5ML4k2Ww+cHKqhX9XDflX1l42sOVdCp0fOHfkLw720FfLS3jNsujuKe6aOG8NXInO9E+7uZKkoBiyVRjd5Afbv6XNzWkGDVQrEpswaApROkGrQjlW2sTaviDxdFmRrKdGl0/G9nEVNH+pja6+v0Bj7aU0bqCG8SQ/rMx3d/Kea/PxdxzcQw/jZvtNx4RuY43rpuwgnPZVUph/BOhharFYqSxg4K6qU5klNG+mIwiDy16Sh+bg785Yq+hjSf7iunqUPDfcZAJ0gp3TVKNbdePMJ0bOORav65JY/5SUE8t2SMLBIyA2Jna8OXt6cOeO6e1YeH+G6GDqsVil63A6S5khsyqsisbOPROXGmTtodPTre3lXMpTF+pER6A1IM4v3dpUT6ODPduNwpiiJv/FxEQrA7/14xzqL1nYxMf6aM9CE24PjMXLXWQE718LQqrFYoNmfXARAX6Ia9woYXv89jQrgni82azny0p5TWLq2FNZFR0UpmZRt/uDjKlJl5uLKNwoYOrk+NkOeAypwWG/r1rejlyY05wzKoaZWfitKmTtO06ikjfXh9eyHNnRr+vjDR9OFXqbW8+0sJ0+P8LfpMrNpdioeTHcuS+5psrT1UiZOdLfOTgpCROR1cHRS8evXx084zKtpYn1F1Du7o7GKVQmHudgS6O/LR3jKumRjGGLN17VW/lqJS60xJWCAlWH2fU8fKSeGmblSdPTq+zaxhXlLQccN/ZGROxkDT31IivHhhax7KruHV1t8qhaI3G9PWRuCHo3U42dvywMy+GRptXRo+2F3KrIQAi1WNj/aWYSMI3DS1rwn45uxaOjV6rp5o3ntHRub0WHvnFIv9P102krYuDa/8eOIJ6taI1QlFWVMnuUa3I9LHmYyKNm65KAofsy5T7/1aQofG0ppoV2tZc6iSeUlBFqna69IqGeHrYioflpH5LYwN87DYX7W7lBunRPLZgXKyh9FyqdUJxWYzt8PL2HRmQnhfDKK5o4cP95Qxb0wQcYHupuNrDlXS0aOzyLQsbuzgUFkry1PC5OVQmd9F/67re4ubWTA2GB8XBx4fRoFNqxMK8/iEp7MUU+htlgvwzi8lqLV6/m9GX1alTm/go71lTIr0Jim0T1TWplViayOwNPlkTcNlZH4bL32fx6Nz4sisbGNtWuWpf8EKsCqhKGvq5GiN5HbYK2xwUNji7+aAn5vkdjS0q/lkXxmLxoWYhgYD/JhbT1Vrt0W6tlZvYEN6NZfH+ltMKJeROVMOlLYAMDHSixe/z6O1U3OO7+jMsSqh2JIjWRPujgpSIrwoauggIbjPvfjfzmK0epG7+9VorNpdSri3s0We/o68Bpo6euQgpswZ4+tqf9yxpzYd5YGZsajUOl4eBoFNqxKK4oZOHBQ2qNQ6SSgaO0xuR62ym88PVHDV+BCLcX+HK1pJL2/llosiLTIu16ZV4ufmwOWxv20sgIxMfyYas37N6ejR8W1WDTdPjeTLgxVkVradgzsbPKxKKBQ2Aj06qduxj6sDeoNosije2lGMwTCwNeHmqGB5Sp/l0KBSsyO/kaUTQuU5oTJnTMoAQgFS1fIVcf74ujrwxMYc9FYc2LSqT4mtrWQRuDko6DUOEoI9qGrtYvWhCpanhFm0169u62arMcHK1Wwg8fqMKvQGkRUpZ2UEqswFhrn728uy5FD8XB14fusxHp0TR1aVkjWHrDewaVVCYWdUh8kjvMmra8fNUUGYtxNv/FyEgMBfzapGQep0BXDT1EjTMVEUWZdWxcRIL0b4WU4Hk5H5PQQPMASqR2fg8fnx5FSrUHVrmRzlzUs/5NFipYFNqxKKepU0iXDqSF+O1qiID3KnXtXDuvQqVk4Ks+ie3dGj48sDFcxJDLRovX+orJXSpk5WpMhBTJnBIcDj+JGSzR09LEgKYnKUN+/9WsqzixNpV+t4+QfrHEFoVUJxoFQa3SZZFCoSQzw4VqtCbxBZOM6y+e26tErae3TcdskIi+NrDlXi6qBgnlwAJjNI9E+6Amju0CAIAskRXtSr1ET7ufKHiyJZfaiSwxWt5+AuzwyrEopWY6GNwsYGtdZAQrA7la1SO7Iws76XeoPIB3tKSY7wsqgcbVdr2ZJdy4KxQaaiMBmZs0Fzp2T9+rg6oDOIqNRa7pkRg7+bA09uPGp1gU2rEQpR7PuPza2VcuilQGY39gobi4ni23LrqWzpPq4x7reZtXRr9bLbIXPWaenUoDeIphyLpo4eXB0UPDE/nikjfdBa2axSqxGK4sYO03ZOtQoHhQ0j/Vyoau0i1MvJYjzgqt0lhHo5WUwDA1iTVklMgKuFlSEjMxgk9ysqNIjQ2qUxfYE1dUhBzPlJwTw2t2+0hLVgNUKxp6hvtHxWVRtxgW4obG2oau22aLefWdnGobJWbp4aaZEjkV/XTmZlGyvkAjCZs8BASVfNHRp8jBZFc4d1rnb0YjVCsbe4ybSdWaUk3piRWdkiWRS9rNpdiquD4rjU7LVpldjZCiwZLxeAyQw+cQNMt2/u6DGzKHqG+pYGFasRiqKGPtdDo5MCmR09Olq7tCahqGnrZkt2LVdPDLPoVqXRGfj6cDUzRgdY9K2QkRkszCfN9dLUqcHL2R5BkETDmrEaofBwsmxTlxDsTnVrN9C34vHxvjIMosjNZglWANuP1dPSqWGFXAAmc5YYKOmquaMHWxsBb2d7GmXXY2jwdrG0BOIC3akyLo32WhSbs2q5Is7fIo0bpNyJIA9HLh0lF4DJnB383QdKupLEwdfVQbYohgofl75SXndHBU72tlS29AqFJAytnRoifFwsfq+mrZtfChtZlhwqz+uQOWsMtIrRl0thT7OVpm73MihCIQjCbEEQ8gVBKBIE4ZEBzjsIgrDGeP6AIAiRv/U5fMxq/ntb3FW1duNoZ4Ovqz1avYFOjf44F2V9ehWiCMuTZbdDZmhpbO+zKC74YKYgCLbAm8AcIB5YKQhCfL/LbgVaRVGMBl4FXvytz+NtZlHEGiPMvUujgiDQrtYBkrXRi8EgsjatkqkjfQj3sXRHZGQGmzBvyzhFfr3Ujc3H1V5eHgUmAUWiKJaIoqgBVgOL+l2zCPjYuL0emC78xmQGc4vCJBRtfUujqm4pvdvdzKLYX9JMVWu33MVKZkjon0tR2SIF231dHejo0aHW6s/FbQ0KgyEUIYB5oX2V8diA14iiqAOUgM9veRLzYGavUFS2dJuEQmkUCnPXY01aJe6OCmYlBP6Wp5KR+V1E+w/ctsA8jdtaOa+CmYIg3CEIQpogCGmNjY0W58yDmS72ClRqLcpurWlpVKW2tCiUXVq25tSxaFyI1aXLylgnAy2Rdmv0+Bi/5KzZ/RgMoagGzG37UOOxAa8RBEEBeADN/a5BFMV3RVFMEUUxxc/PcinT3PXQGQymHIreFQ9Vd2+MQhKKPcVNaHQGFo+3LD+XkTlbDJR01dzZg6+b9WdnDoZQHAJGCYIQJQiCPXANsKnfNZuAm4zby4CfRfNy0NPAvOZfZxDNlkaNMQqTRSEFMzXG3pq9Q4JkZM42QQMJRYfGZA1f0BaFMebwF+AH4BiwVhTFo4IgPCMIwkLjZasAH0EQioD7gOOWUE9FnVJt2tYbRKpMFsXAMQoRSYds5AIwmSEiwL1PKMYaB2Y3d/bVezRasUUxKN1bRFHcAmzpd+xJs201sPxMnsNcrXV6SSic7W1Ny6aqbi0KGwEnYzzCYCz3l3VCZqgwj4X5uzsCSpo6NDjZ2+Jib3thWxRDhZdZMFOyKKSl0d5VVpVai7uTnWm/168RkJVCZmgwb5zbG4/oFQcfVwdTpqY1YjVCYY7WYKCyXx8KVbfOItmqNwQiWxQyQ4X5XNzRQe4429uaBMPX1f6CD2YOOXq9ZFGEmfWhUHZrLXIoTBaFLBQyQ8SGjCrTtlZnMGZk9vXOlF2PIeZojYp2tc7SojC6Hr30WRSyUsicfapauzhc0Tc2sLylS6oa7Rwe9R5WJRQTwqVel69uLwCw6Gyl6taacigAehdfZZmQGQq+zay12M+qasPHxcHUK9PX1d7UcNcasSqhmBRlmfVtaVHoLC0K40/ZoJAZCj7ZVwZAfJBU2azWGvA1dz1c7DGI0NZlne6HVQnF+HDL7tnm1XrKbq0p2Qr6LAo5j0LmbJNf106tMc9nhlnnd09nyYowGESz7ExZKM46/XPpe4OXaq0ejc5g4XoYemMUQ3d7MhcoG49IFQu+rvbk1ihNxwUB0/CfvnoP64xTWJVQeLtapmNnVkl/lP4FYdDneshKIXM2MRhE3tpZDMCS8SHsyG8kytfFdA4kK8LPTXrvWmt2plUJhXkFKcD6dKm6vbcgzKK7lcmikJVC5uyRYTZH1CBKyYB3XTbSuC+9B5s7eqy+gtSqhKJ/ufhn+ytQa/WmOg+LhCvjT7lNpszZ5OvDktsxMdKLnfkNpER4kTpCCrr3Tg1s6tDg4WSHrY1gtdmZViUUA/HTsYYBXY9es0/Oo5A5W2j1Bj4/UAFAYogHxY2dFiMhnOylj1dzZw82NgI+LvY0tcsWxZDg72bZFn19emVfGzzHAZZHh+rGZC44dhf2Ta9r7dTgbG/LvDFBpuDmsuQwBKFvpcOa6z2sTigSgt0t9ncVNFJsnCJmkcItL4/KnGV6rYl5Y4LYllvP/KQgnOxsWZtWxZQRPkT5uuDt3JdL4etqvYOArE4ogjwtl0gNIny8rxwAtwFiFLJJIXM26NLo2H6sHpCaJXVq9KxICWN/aTMVLV2smBgKWHbgtuZBQFYnFEsnhFrsuzsqUHZrsVfYWAQ75epRmbPJttx603ZRQwcjfF1IjvBiXVoVbo4K5iQGAVJntt5gu68Vt+23OqFIjvCy2B9p7Hzc2/quF7nWQ+Zs8j9j7sTiccEcKmtleUoYKrWOLdm1LBwbjKOdLTVt3eTUKJkUJbXx93F1oFurp7NHdy5v/XdhdUIBcKVZmmyw5/Gdj0FuhSdz9mjp1JBX1w5Irq+tjcDSCSF8m1lDj85gmiPz9eFqRLHPCrbm3plWKRR/nDbCtN0byATo0fUNWOkt0pN1Qmaw2WxsUOPqoOBAaTOXxfjh7+7I2rRK4gLdGBPigSiKrE+vYlKUt2lKnco4zc4a35NWKRTRfm6m7by6diKNf4iNR2pMx/tcDyv8q8ic1zy96SgAl8b4Uq/qYcXEMI7VqsiqUrIiJQxBEMioaKO0qZNlyX0xte+yahgd5E6Yt/WNt7RKofBwtjM10QVYOE4aTPbQ+izTsV7XwxrVW+b8paq1y9RTol2tw9fVnivi/FmbVom9rQ1LxkvvxfXpVTjZ2TJ3jBTUrGyRGtssHGudc2asUigAEkP68il6hwEBNLZLy0+i7HrInAXWpUnt7kb4ubCvuJkl40MwiCJfH67myoQAvFzsUWv1fJdZw5zEQFwdpCX7TZmStbtgbNA5u/czwWqFordCD6RehXHGeaTfGHPvRbkoTOYs8J+fCgGI9HFBZxBZnhLG9twG2rq0rEiRgpg/5tbT3qOzcDu+zawhOcLLotmSNWG1QjHCz3IgbO8f5bktxxBFUbYoZAadvDqVabuypYtxYZ7EBLixJq2SYA9HLo72BSS3I8TTyVQcVlDfTl5du9W6HWDFQtFrUSQZJzLtL2kxnTtao5JrPWQGnfd/LQWkVvyFDR2sSAmjuq2bXwsbWZYciq2NQJ1Sze7CRq6aEIKNsXR505EabARM8QprxGqFotfV6O16tf1YPcHGaWLr06vkWg+ZQcVgkJY7ATyd7HC0s2H+2CA2GN9ry1P6cicMZrkToiiyKbOGi6J98etX0GhNWK1QRPi4MGN0ALsKGk3HekvKP9lXhtbYDECjNwz06zIyv4l0swY1OdVK5iYG4WqvYG1aJRdF+xDm7WzMnagkJcKLSKPFm1mlpKKliwVW7HaAFQsFwEOzY1GbJVlVt0mrHwYR0zwFcyGRkfm9PL/lGCBZsu09OlZMDGN/STNVrd2mIGZmlZLiRsvciU1HarC3tWFWQuA5ue/BwqqFIibAjavGhw54rl6lxtfVnq8zqof4rmSGG1q9gQzjcB9BEIjwcWZylDdr0ipxd1SYRGB9eiUOChvmJkmxCL1B5LusGi6L9bNs02iFWLVQAPzfjFEDHt+Z30ByhBc/5zWg7NIO8V3JDCd25DWYto/VqlieHIqqW8fWnDoWjQvB0c4WtVbPt5m1zE4MNDVQOlDaTEN7DwvHWbfbAcNAKMK8nZmfdHw02SDCD0fr0egNbMmpHeA3ZWROj/vWZgIQE+CKIMDS5FA2ZVajMSsA++lYA8pu7XG5Ey72tkyPCxjwca0JqxcKgKcXJhx3bGyoh6mx7ob0quPOy8icDl0aHR3GsvC2Li2XjvIjyMOJNWmVxAe5kxgiLc+vT68k0N2RqSOlXAqNzsCW7DqujA/Ayd72hI9vLQwLofB17Vt2GhsmTRNraO/hqQWSgKSVt1LV2nVO7k3GuulN2QbpPbUiJYyjNUpyqlWsSJGshwaVml8Km7hqQgi2xm+nXwsbUXZrh4XbAcNEKADuuzIGgGM1UvZcrVLN3DFBzBjtD8CKt/eds3uTsU5yqpU8ZawUHeHrgqezHTPi/VmXVoW9wobFxgKwb45UozeILDVf7ciswcPJjouj/c7JvQ82w0YoLhllNPnM8ia+P1rHOzekAFCjVLMzv2HA35WR6c+hshZWvrvftF/V2s3icSGIopRUNSshEE9ne0RRZEN6NePDPRlpLCvo1ujZllvP3DGB2CuGx0fsjF6FIAjegiBsEwSh0PjTa4BrxgmCsE8QhKOCIGQJgnD1mTzniYjv150b4IlvcrC1EXhsbhwAN394yCJfX0ZmIHYVNHLDqgO0G2MToV5OaPQGVqSEsS23HmW31uR25FSryK9vtwhibj9WT5dGb/VJVuacqdw9AvwkiuIo4Cfjfn+6gBtFUUwAZgOvCYLgOcB1Z4SDwnbAP0yPTm9KiAG49aM0GtrVg/30MsOErdm13PbxIUb4urIiRarf0OgMjAnxID7YnbVplYR4OnGRMWi5IUNyQ+Yn9b33NmXW4O/mwOQon3P1MgadMxWKRcDHxu2PgcX9LxBFsUAUxULjdg3QAJwVx+3xeaNx7hdhfun7fDyd7Zlp7LPZ2NHD7R+n0a3RD/QQMhcw69Iq+fMXGSSFevLF7ZP5paAJX1d7YxAzlKrWLnYXNbEsORQbG4EenZ5vjlQzMz7AlFCl7NayK7+R+UnBpsDmcOBMhSJAFMXeJIU64KQLxoIgTALsgeIzfN6Bb8bdkYdnx1kcW7Vbqvjr7Tx0zcQwsqqV3LvmiGnsoIzMB7tLeXB9FhdF+/LxHybxrx8LqFOpEUWwV9iwcGyIqShsudHt2JEn9aEwD2L+kFOHRm8YNqsdvZxSKARB2C4IQs4A/xaZXydKnWJO+MkTBCEI+BS4RRTFASu1BEG4QxCENEEQ0hobf1+NxvWpEabS816+PlzF5XH+uDkq6FDreHxePN8frePF7/N+13PIDB9EUeQ/2wt55rtcZiUE8Pb1yfzt690gTU4AACAASURBVGw+3V/OzVMjUWv1zE4IxM1Rwbq0Ki6O9jU1n1mfXo2/mwOXGPtQgOR2RPg4M7bfe9DaOaVQiKI4QxTFxAH+bQTqjQLQKwQDLisIguAObAb+Jori/oGuMT7Xu6IopoiimOLn9/u8E1sbgX8uGWNx7N41mezIa2DemCC+P1rHyklh3JAawTu/lPDlwYrf9Twy1o8oijy3+Rivbi/gqgkh/GvFOP765WE2HqnhodmxTIjwQqXWsSIljL3FzVS3dZvKyZs6etiZ38CSCSEobKWPUUO7mr3FTSxICh52w7HP1PXYBNxk3L4J2Nj/AkEQ7IGvgU9EUVx/hs93WiSGePCHi6Isjt29+jAeTnZ0GZeunloQz7QYPx7/Jsdi2KzMhYHeIPLoV9m8v7uUm6ZE8NSCBP7w0SF25Dfw3JJE7rosmnXGwOXUkT6sSavEw8nOFOvaeKQGnUFkmdnkui1ZtRhEhp3bAWcuFC8AVwqCUAjMMO4jCEKKIAjvG69ZAVwK3CwIwhHjv3Fn+Lyn5L6ZMRb7o4Pc+WCPFK/4+nA1Clsb3rh2PKP8XfnT5+kU1ref7VuSOU/Q6Azcs/owqw9V8tcrorl7+iiue38/GeWt/Oea8Vw3OYLKlr7ApUqt5YejdSwZH2IaW7k+vYqxoR6MCugbHbEps4a4QDdizI4NF85IKERRbBZFcbooiqOMLkqL8XiaKIq3Gbc/E0XRThTFcWb/jgzGzZ8MVwcFb1+fbNp//qoxRPtLf8Cd+Y00tvfg5mjHqpsn4mhnyy0fHTJ18JYZvqi1eu78NI3vsmp5bG4cKyeFs/ydfRTWd/DejSksHBtMvUrNLR8dwt7WhuUpoWw8UoNGZzAFMY/WKDlWq7IIYla2dJFR0TascifMGR5pYydgVkLfIsxVb+3ls1snmfaf+S4XgBBPJ1bdlEJTRw+3f5KGWisvmw5X2tVabvzgIDsLGvnnkjHMGB3A8rf30ajq4dNbJ3N5nD+VLV0sf3sftW3dfPyHSYR6ObPmUCWJIe4kBEsByg3p1djb2rDALHfi2yypHb81N9A9GcNaKARBYNu9lwLQozOgUus4+Nh0QCoBzqlWApAU6slrV48ns6qN+9dmysumw5CWTg3XvX+AjPJWXrt6HGPDPFjxzj7UWj1f3pHKpChvihraWfb2XpTdWj6/PZXUET7kVCvJrVWZkva0egMbj1QzI94fL+MsUZA6WY0P97TKKWCnw7AWCoBRAW7Y2UoR6Mtf2Ym/u6Mp0Dn/v7tNKd2zEwN5dE4cm7Nr+de2/HN2vzKDT71KzdXv7COvrp13bkgm+P/bO+/4qKq8/79vek9Ib6QSSkgIoSUICLiANAVRUAQUQYFVH9vaVn/uPrvrs4+Lj7qu4qIUUYoUFRVBpFdJ6EkgEBJCep/0Saaf3x8zGTJkQthdSYDc9+s1rzn33DOT79zM/cwp3/P9ejnzyGcpONjasHnJcOJCPDlXXMesT1PQG2DT4mQGmnYhbzpRiIOdDdMSjH44B7IqUSg15uC5ANl3QDj+jrjjhQIg888TzeW1x/JY3CrJ8ZwVqeSYEh0/NSqK2cN6smz/ZTafLOxsM2VuAgWKJh5a/gsltc2seWIotjYS81al4ufmyJbf3kW0nxun8quZvSIFJzsbtiwZTt9A476hr08VsS41n2kJwXi62JvqCvF1c+Tu3leX77elGcPxT7ESQOlOoVsIhb2tDe/PSgDgre/P4+lsb95tKkkwZ2UK+QolkiTx52lxjIrx5Y1vM/jlsrxsejuTXd7AzE9/oUGlY/1TySgaNTz15Umi/dzYvGQ4IV7OHMmuYu7K4/iahKMlX8zmE4W88nUaI6J9+fO0OMA4fNl3sYLpA4OxN/lOtITjHx7tg7+7U5d91ptNtxAKgBmtuop939rJdFNi4xfH90ajM/DoilSKapqwt7Vh2ZxBRPq6smTtKXNvQ+b2IqOojlmfHsMgYNOi4Vworee5jWcY2NOLrxYl4+vmyO7MchasOUG4jwubFicT4mXMEbMhtYBXv0lnVIwfKx8fYo5Q9e3pIrR6y7gTGcV15Cma7uhhB3QjoQA4/dZ4c/mPP5zHwc6GzJJ61i5MokGlZc7KVMrqVHg42bN6/lAc7GxYsOYEikZ52fR2IjVXwewVKbg62vH1kuHsz6rg999mMLq3H18uSMLDyZ7vzxazZN0p+gV7sHFRsrk3sDYlnze2ZjC2jx+fzRts9ps4kl3F0p1ZJEd50y/oakiDH86WYG8rMbH/nTvsgG4mFN6uDix9aAAAjWodGp2Bb08X0zvAnS8WDEPRqOHRlSlUNqjp6e3CiseGUF6vYtHaU/Ky6W3C/qwKHlt9nAAPR7YsGc7GE4W889NF7ksI5rN5xt7BhtQCXth0lqERPVj/ZBJeLsbVizVHr/DWd+cY18+f5a1E4kReNU99eZIoP1cL35ycigY2nSxkTB9/8xzGnUq3EgqAmYNDzakHAZq1en5MLyExrAer5w+ltFbFvFWp1Cg1JIb14P1ZAzmVX8P0ZUc5JCcTuqX5Mb2Ep744SS9/NzYuGs5H+3L454HLzEkK4+8PD8TBzoYVh3J5Y6uxd7HmiWG4OdoBsPJwLv+9LZMJsQF8MmcwjnZGkUgvqmXB5ycI8nJi7cKrolJRr+Lx1SdwtLPlD1Nju+wzdxbdTigkSWLDU8kWdS9tTqOktplhkd6sfHwIuVVK5q1Opa5Zy5QBQSyfOxilRsdjq48zb1Wq2f9C5tZh04kCnvvqDIlhXqxdmMSftp1nQ2oBT4+J5u3pcdhI8MHuS/zPjgtMiQ/is3lDzD2Gzw5d5u3tF5gUF8iyOYPM4euyyhp4bPVxPF3sWf9kkjl3qFKtY8EXJ6hWalg9f8gd6zvRmm4nFAARvq78brzlXpC73tlHgaKJEb18+XTeYLLKGpj/+XEa1TomxgWy56XRvDU1loziOu77+AgvbTorR/a+RVh5OJfXvslgZIwfn84bwkubz/JjeimvT+rLq6b4JG9vv8CHe7OZOTiUf8xONIvBJwdy+OuOi0wZEMQ/ZieaVzOuVCmZszIVRzsbNjyZTJApGbZOb+DZDafJLKln2ZxEBoT+6sHabkm6pVAALBodha+bg0Xd3e/uJ7OknrF9/Plo9iDSi+pYsOYEzRo9jna2LBwZycFXxrJkdDTbM0q5572D/HXHBTkTWRchhOCD3Zd4e/sFJscH8sGsBBavPcnBS5X874x4loyONu8SXXXkCvPviuBvDw4wR576aG82S3dmcX9CMB8+PNAsEkU1TcxZkYJBCNY/mUSYj4v57731/Tn2Z1Xyl+lx3HMHJPa5UbqtUDja2fLxo4Pa1E/+x2FScxVMjAvkg4cHctI0kdUymenpbM9rE/uy/+Ux3J8QzIrDudz97n5WHMqVJzw7ESEEf/nxai/hv+/rz7xVxzlbWMtHsxOZPSwMrd7AC5vOsvFEIc+O7cUf74vFxkYyC8x7uy8xIzGEDx4eeDWmRL2KOStTaVTrWLtwmHkjIcAnBy7z1fFCnh4TzZyk8K766F1CtxUKgOQoH4voyS08/FkKu86XcX9CMEsfSuBIThVPrz+NRnc1MFewlzP/NzOBHc+NYmBPL/5nxwV+895BvjtTLO8VucnoDYLXvkln9dErPDEigufHxfDIZynkVhl3gE4dEIxKq+e3606xLa2E1yf15eV7+yBJRpF4f/clPtybzUODQ3l3ZoK5h9GyH6SyQc2aBcPMm8DAGCXt3Z+zmDYwmJcn9Omqj95ldGuhAHhjcj/cTTPfrVm09hSbThTw0OBQ3p4ex76LFTy/8Qw6vWUUv35BHnyxYBjrFibh5WLPC5vOct/HR+RgODcJjc7Af311ms0ni3juNzHMTQ5n1vJjVDYad4CO6eNvnGxcc4I9Fyr4y/Q4loyOBoy9kHd/zuKjfTk8MrQnS1sNQ+qatcxblUpBdROrHh/KoLCrmSd+yani1a/TSY7yZulDA7C5g4Lm3ijdXii8XR14d2aC1XOvfZPBsv05zEkK462psfx0royXNqdZHWKMjPFl27Mj+fCRgdQ2aZm7KpXHVh8ns0TOI/Jr0azR89SXJ9mRUcb/m9KPCbHGbeIavYGNi5IZGuFNnenap16p5v1ZCcxLNg4RhBC889NFPjlwmUeTwvjrA/HmG75FWC6VN7B87mCGR18Ns59V1sDidaeI8HHl03lDzMum3Q3JGBP31mPIkCHi5MmTnfb3ThfUMOOTX6yee2JEBG9NiWX5ocss3ZmFr5sDC0dGMTc5DHento42Kq2edSn5fLQvh3qVlhmJofxuQm+CTS7CMjeOWqfnl8sKdmeWszuznKpGNf/7QDxRfm4sXHMCD2d71i4cRpSfG1WNauatOk5ORQMfzR7ExLhAwCgSb2+/wKojV3hseDh/ur+/OaalSqtnwZoTpOQqWPboICbFX/WwLK9X8cCyo+gMgq3PjDC7eN9JSJJ0SggxpMN2slBcpaJexbC/7rWoc3O0o1GtY9rAYN59KME4WbYvm8PZVXg42TF/RCRP3BVhEZughbomLZ8cyOHzX/IAWDAikt+OiTbngJCxTr1Ky/6LFezKLOdgViWNah2uDraM6eNvnqRcsu4UoT2cWbswiWAvZ0rrmpmzMpWS2mY+nTeE0abdnUII/rQtkzW/5DH/rgj+eF+sWSQ0OuP77LtYwfuzEiz2AzWqdcxafox8hZJNi4ebs5bfachC8W9S1ahmyNt7LOru6evPvosV3N3bj+VzB+HiYEdaYS3L9uewK7McFwdb5iaH8+SoSKs7CItqmnh/1yW2ni3G09meZ8f2Yt7w8G7bjbVGWZ2K3RfK2XW+jJRcBVq9wNfNkfGxAUyIDWB4tA9O9rb8kFbCS5vO0jfInS+eGIaPmyP5CiWPrkilvlnL6ieGMjTCGwCDQfCHH86xLqWAhSMj+X9T+plFQm8QPPfVGbZnlPL29DjmJl9dxdDqDSz84iRHc6pY9fgQxvTx75Jr0hnIQvEfcCKvmpnXZD9/7p5efLw/hwGhXnw+f6i5B5FV1sAnB3LYllaCna0NDw/pyeLRUebcD605X1LHOz9d5HB2FaE9nHnl3j7cNyC4W06OCSHIqWhkV2Y5uzLLSSusBSDS15UJ/QOYEBtIYk8vi2uzLiWft74/x9AIb1Y9PgR3J3sulTcwd2UqWr2BLxckEW/Kp2EwCN787hxfHS9g8d1RvD6pr1kkDAbBq9+k8/WpIt6c3I+n7o6ysOu1b9LZfLKId2bE88iwsE68Kp2PLBT/IR+YltBa8/eHB/LqN+n07OHMlwuTLMaseVVKlh+8zDenixACpieG8PSYaKJMGa5bc+hSJf/700UulNYTH+LJ7yf35a5o3zbt7jQMBsGZwhp2nTeKw5UqJQAJPb2YEBvAvf0DiPZzs8iJ0azRc/BSBdvSSo1Obn39+WTOIJzsbckoquOx1anY29qw7skkc/Rrg0Hw+rfGm/3pMdG8YloaBaMQ/PGH83x5LJ8XxsXwwjhLD91/7M3m/d2X+K97evG7brAMKgvFf4hOb+Dhz1I4lV9jUb/+ySSWrD2Fm5MdXy4YZhGuHaCktpnPTImFNHoDU+KDeGZsL4utyWD8Mn93tpj/+zmLkjoVY/v48fqkfvQJvLNCvau0eo5dVrArs4zdmRVUNaqxs5EYHu3DhP6BjO8XQKCn5XCtSaNj/8VKdpwrZf/FCpo0erxdHXhocCiv3NsHe1sbjl+pZsGaE3iZ9mGE+xgDzugNgle/Tueb00U8d08vXhzf20Ik3tl5kU8P5rLo7ih+36qXAcaIVi9vSWNGYgjvzUq445L4WEMWil+BwuomJn94mEg/V9KLrm4E2/BkEs9vOotGZ2D1/KEMDu/R5rVVjWpWHbnC2mP5NKp1jOvnzzNje5EYZtlWpdXzxS95fLw/B6Vax0ODQ3lxfG/z3oLbkbpmLQeyKth1vpwDWRUoNXrcHO0Y08eP8bEBxm3Z10zoKtU69mdVsCOjlP0XK2nW6vF1c+De/oFMiQ9iWKS32Xvy4KVKFq89SbCXM+ufTLLYh/HyljS+O1vCi+N68/y4GIu/8dHebN7bfYk5SWG8PT3OQgiOZFcx//PjJEV58/n8Yea9IHc6slD8Snx/tpjnN54lOcqblNxqc/3yuYN456eLlNWr+PvDA7m3f6DVX6C6Ji1rfsnj81+uUNukZUQvH54Z24vhUT4W7WubNCzbn8MXv+QjSTAkogfBns4EeTkT4uVEkKczwV7OBHs54eLQ1kGsqymta2Z3Zjm7zpeTkqtAZxD4uVtORl47eduo1rHvYgU70ks5cKkCldaAr5sjk+ICmWwSh2szgu88V8p/fXWGGH93vlw4DF83445Ond7Ai5vT2JZWwssTevPsPZYisfJwLm9vv8CMxBD+b2aCxdzHhdJ6Zi4/RoiXM1t+OxwPK0vedyqyUPyKbD5ZyJ9+OI/OIFC3cuN+a2os354u4nxJPb383ZibFMaMwaFWv2hKtY71qfmsOHyFygY1g8K8ePaeXozt428hGIXVTXxy4DIXy+opqW2mokHNtf8iLxd7o3B4OhHs5UyQlxMhXs4mMXEiwMPJvMHpZiGEILuikV3ny9iVWW7ucUX5ujKhfyAT+gcwMNSrzURtg0rLvosVbE8v5eClStQ6A/7uRnGYFB/E0Ii24mAwCE7kVbP1TDFbTBm6Pn9imLlXotUbeGHjWbZnlPLaxL78dky0xes3pBbwxtYMJsUF8tHsRHPPBIwC98Ayo//Mt0/f1e18XWSh+JUpUDTx0uaznLxmzmJecjj9gz346kQhaYW1ONvbMj0xmDlJ4VbX3lVaPVtOFrL8YC7Ftc3EBnnwzNheTIwLbHODgPEmKKtTUVqnoqS2mZK6ZkpqmymtVVFiqqtrtty9KkkQ4O5EkJcTwSbxaN0jCfZyxsfV4V8eg+sNgtMFNaaeQxl5CuM2+4E9vcwrFb38207e1qu07L1Qzvb0Mg5lV6LRGQjwcGRSXBCT44MYEt7D6srP5cpGtp4uZuuZYoprm3FxsOX+hGDemhqLq8ntvsWl++fz5W1WMAC+O1PMi5vPMrq3H5/NG2IxpKhXaZm1/BhFNc1sXjyc2GDLeaTugCwUNwG9QfCpyTuzNXdF+/Dxo4MormlmXUo+36cVo9IaSAzzYm5SOFMGBJmDpLSg1Rv47kwx/zxwmdwqJdF+rvx2TC+mtYrwfKMo1TpK65opqW0RExWlJlEprVVRXNts0RMCcLCzIcjTyTS8seyRBHs5E+TphLuTPSqtnqM5Vew6X86eC+UolBrsbSXuivZlfGwA42MDCPBo6ztS16xlT2Y5OzJKOZxdhUZvIMjTySQOgQwKsy4OikY129JK2HqmmLSiOmwkGBnjx4zEECb0D7AYdql1ep5Zf4Y9F8p5a2osC0daJqbeea6MZzacZmhED9Y8Mczif6DRGcwemZ8/MZRRMX50R2ShuImcL6ljyj+OWNSFeDnz6bzBxIV4Utek5ZvTxpwQuZVKerjYM3NIT+YkhZln51vQGwQ/nSvl4305XCxrILSHM4tHRzNzcGgbcfl3EUJQ06Q1ikhtc6veicrUO2mmrF7FtZte3Z3s0OkFzVrjZOTYvv5MiA1gdB8/q8OruiYtuzLL2JFRypGcKrR6QYiXs3lYca1fRAsqrZ49F8rZerqYg5cq0RkEsUEezBgUwv0JwfhbEaKyOhVvbM1g38UK/nR/fx6/K8Li/IGsCp768iRxIZ6sXZhkDnnXcj1+tyWNb08X8+5DA5hpygLWHZGF4iaj0urp+9ZOizpHOxv+9uAApicaUwEIITh2WcG61Hx+Pl+O3iC4u7cfc5PCuKevv8VYWQjBvosVfLQvh7OFtfi7O/LUqCgeHBxKDxf7m75Up9MbqGhQW/ZIapuRJImxff1JjvK26klao9SwO7Oc7RmlHM2pQmcwisOUAUFMigtkYE8vq7YbDILjedVsPV3MjoxSGtQ6Aj2cmJYYzIzEUKvLxAaD4OjlKtal5LPnQgUGIfjztDjzxq8WUnIVPL76ONF+bnz1VHKbwLfv777EP/ZmW/Wj6G7IQtEJGAyCQW/vpvaaCFcLRkTyxuS+FkJQXq9i4/FCNhzPp7xeTbCnE7OHhfHwsJ4Wbt9CCH65rODjfTkcy1UAYGcj4eXigLervfHZxYEersbjHi4O9HBxwNvVVOfigJerPe6OdjdNXKqVGnadL2N7RinHLhtXOEJ7ODMl3jjnMCDUs92/nVPRyNYzRXx3poTi2mZcHWyZGBfEjEEhJEf5WJ2nqVFq2HKqkA2pBeQpmvB2dWDmkFAeHda2h3a2sJY5K1II8nJm06JkfEyrIi1sPlHIq9+kM3NwKEsfGtAtfCWuhywUnUSTRseEDw5RVNNsUZ8c5c2yRwe1+aLq9Ab2XKhgXUo+R3KqsLORuDcukHnJ4SRFelt8cc8U1HAqv4aaJg3VSi01Sg3VTRpqW46bNOjbCZJjZyOZhaNHi6C4WopMi+i0iIyrg227N46iUc3P541zDsdyFegNgjBvFybHBzElPoi4EI92X1vVat4h3TTvMCrGjxmDQhgfG2B1uVcI48TpupQCtmeUotEZGBrRgzlJ4UyKD2zTu1E0qo3L0Efz8HZ1YPPi4W0cuQ5eqmTBmhPcFe3D6vlDb/rK0O2ALBSdSGldM89/dZbjedUW9cGeTvxuQh/G9w+wOqbPrWxkQ2oBW04VUdesJcbfjbnJ4TwwKOSG1vKFENSrdNQoNdQ0adoRFA01Si3VTRpzu/YCcDnY2uDlYm8UDpOAeLnYc6VKSUquAoOACB+jOEyOD6J/cPvioNLq2Z1ZztYzxnkHvUHQP9iDBxJDuH9gcLvp9xrVOraeKWZ9Sj4Xyxpwc7RjxqAQHk0KM+cEbU1hdRMrDuey6UQhGr2B8f0C+MN9sRZ7bXR6A9+eKeZPP5wnzMeVzYuTrYYH6I7IQtHJCGF0yX5xU1qbcw62Nozp48fUhGDG9fNv8wuq0urZllbCupR80orqcHGwZdrAYOYmh1uEY/s1MBgEDSod1WYRub7IVCs1eLs6mJcy+wW5tysOBoMg9Uo1W88U8VNGGQ1qHUGeTkwbGMKMQSHmvRjWyCypZ11qPt+fKUap0RMb5MHc5HCmDQw2L4Ve2375wctszyjFRoIHEkNYdHe0xfKsEIJdmeW8+3MWORWNJPT04tO5g9v0NLozslB0EfUqLR/svsTnR/PMdYEeTggE5fVqnO1tuaefP/cNCGZMH782KxvpRbWsS8nnh7QS8xLrvORwJse3XWK9VcipaODb08V8f/bqvMOk+CBmJIaQ1M68AxgFcnt6KetS8zlTUIujnQ33JQQzJynM6iSoEIKU3Gr+efAyhy5V4upgy5zkcBaMiGxz86fkKvjbzoucKaglys+VV+/t0673bHemU4RCkiRvYBMQAeQBs4QQNe209QAyge+EEM929N63q1C0cKG0nrkrU1EoNQD08nfj95P6ciCrkh0ZpSiUGtwc7RgfG8B9CUGM7OVn4QzU3hLrzMGh+Lk74mhni4OdTbs34c2mqlHND2eN8w4ZxXXY2kiMivHlgcQQJsQGmhP7WuNKlZL1Kfl8fbqI2iYtUX6uzEkK58FBIeZMXK0xGAS7Msv458Fc0gpr8XVz4IkRkcxNCm+zonG+pI6lO7M4eKmSQA8nXhgXw0ODQy0mlmWu0llCsRSoFkK8I0nS60APIcRr7bT9EPAztb/jhQKMv4BbThXx6tfp5rrlcwcxrl8AKbnVbEsrYef5MuqatXg62zOxfyBTE4IYHuVj/mK3LLGuTclnV2Z5m8lLe1sJB1sbHO1tcbSzMT1scbRvVbazMR3bmtpaP3dt2cFKm4ziOraeLuJQdhV6gyAuxIMHEkO5PyHYnEnLGlq9gT2Z5axPLbg6ids/kDlJYQyP9rH6S6/W6dl6upjPDuWSW6UkzNuFRXdH8ZAVH5N8hZL3dl3ih7QSPJ3teXpMNI/fFXHL9sJuFTpLKLKAMUKIUkmSgoADQog2m/glSRoMvALsBIZ0F6FooV6l5fHVxzlTYAzO4u5kx4k3x+Fkb4tGZ+BITiU/ppWyK7OcRrUOH1cHJsUHMnVAMMMivM1OSuX1Kg5kGbddq3UG1FoDap2prNObjo1lja6l3PZc69f+O5kFgj2dmJYYwozEkDbb7K+lpLaZjccL2HiikIoGNSFezswe1pNZQ3padaQC436QDakFrDpyhYoGNf2DPVgyOppJcYFtegYVDSo+3pfDhtQC7GwlFoyIZPFoOdzgjdJZQlErhPAylSWgpuW4VRsbYB8wFxhHNxSKFk7l1/DgP68G8P3bg/E8PPRqBCWVVs+BrEq2pZew90I5Kq1xT8Tk+CDuSwgmsR3npX8XIQQ6g2glKu0LilpnQKMzEOjpZCFe1jAYBIeyK1mfWsDeC+UIYExvP+YkhTO2r3+7w6WKBhWfH81jXUo+DSodI3r5sGR0NCN7+bb53PUqLSsO5bLy8BU0egOPDO3J87+JaVd8ZKzzqwmFJEl7gEArp94EvmgtDJIk1QghLAIuSJL0LOAihFgqSdJ8riMUkiQtAhYBhIWFDc7Pz+/I/tsOIQTPbzzLD2kl5roPHxnYZlyvVOvYe7GCH9NKOJBViUZvIMTLmakJQdw3IPi6S5NdhaJRzeaTRWw4nk9hdTO+bg7MGtKT2cPCrpvIN69KyWeHc/n6VBFavYFJcYEsGR1tNa9nS4TzZftzqGnSMnVAEL+b0IdIX1cr7yzTEbfM0EOSpPXAKMAAuAEOwCdCiNev9953Yo+iNQWKJu5+d79FXXKUN6N7+zMqxpfYIA/zr3a9Ssvu8+VsSy/hSLbRTTrS15WpA4w9jestO95MWnokZwpqWZ+az08ZZWj0BpIivZmTHM7E/oHXDQCTUVTH8oOX+elcKXY2xXmrxgAACMpJREFUNjw4OJRFd0dZvelbfCH+vvsSJXUqRsX48uq9fc0xMmX+PTpLKN4FFK0mM72FEK9ep/18uvHQwxpnC2uZvuxom3pvVwdG9vJlZIwvo2J8zVGcapQadp4v48f0Eo5dNjpB9Q5wY+qAYIK9nNHqDaaHMD7rDGgNrcp607G1cstr2i1frdPpBZpWWdPcnex4cFAoc5LCrjtvIYTgaI6C5QcvcySnCndHO9MSZ4TVYUMbX4hQT16b2Je7et35MUY7g84SCh9gMxAG5GNcHq2WJGkIsEQI8eQ17ecjC0UbWpy1/mf7Raoa1bg42JIU6U1GcT1VjWrAuLw6yiQaSZE+uDraUdGgYue5MrallXAiz+qqtBkHWxvsbCXsbW1MD+ma56tlBzsb7GxM9XY22LcqO9iaztkZX+NgKxHk6cyk+MDrRt5q2SX76cFcMorr8HN3ZOHISB5NCmvXC/VaX4hXJvRhYpzsC/FrIjtc3YbUq7T8fXc2XxzLw8PJjlfu7cvAnl4czaniUHYlx69Uo9YZsLeVGBzeg1ExfoyK8aV/sCc1TRqaNXqrYmBnI3XZzaXS6vnmdBErDuWSp2gi0teVRXdH8UBiSLtLl+dL6nj35ywOZMm+EDcbWShuYy6W1fOH785zPK+ahJ5e/GVafwaEeqHS6jmZV8PhnEoOX6ois9SY19TLxZ4RvXwZ1cuXUb39Oi31nU5voKZJS7VSg0KpNu4pUapRmFzDFUoNKbnVVDWqSQj1ZMnoaCb0tx7JC4zzNu/tzuL7s7IvRGchC8VtTuvhiEKppre/O/GhngwI9SQ+xJN+QR40qnUczanicHYVh7MrKa83DlOi/FyNohHjR3K0j0XQluvRpNGZ93e03Oytywrl1f0hCqWmTQi+1ng42eHj5ki0nysLRkS261QFUNmg5qN92bIvRBcgC8UdQr1Ky9pj+ZzMqya9qM7sEm5nI9E7wJ34EE/iTeJhayOReqWaw9mVpOZW06zVY2cjMSisByNjfPFzdzTe+I3GTWDGG9/YE1Ao1ai0Bqs22NlIeLs63PCjh4tDu1u4DQZBRYOaK1VK8hRKLpbWs+VUEWqd0Rfiud/EWA2tJ3NzkIXiDkQIQWmdivSiOjKKa03PdebAOfa2En0DPYgP9aRvoDsanTFq1bHLCs6V1Jmjebs62OLt5oC3qyPeLvbGZ1fjs09LABxXB3PZw+lfC4IjhHEDXJ5CSV6VkisKJflVTcZjhdJCkBxsbZjQP0D2hegiZKHoJgghKKppJqO4zkJAGlQ6wBhEt1+QB2HeLsT4uzG6tx/9gz3+44lBIa72DPIVSq5UNZFn6iXkK5po1urNbR1sbejp7UykryvhPq5E+LoS6eNKhK8LQZ7OXbaxTUYWim6NEIJ8RRMZxXUmAanlXHE9jWqjeDjZ2xAb5MGAUC/iQ4zzHlF+bm1uWCEElWYxaOJKSw+hqq0Y2NtKhHm7EGESgghfVyJ8jMfBXrIY3KrIQiFjgcEguKJQklFk7HmcK67jXEkdTRrjze7iYEtcsCe9A92oVmq4UtVEvkJpPg9GMejp7UKkj7FnEOnrYhIEWQxuV2ShkOkQvUGQW9lonutIL6olu6IRPzdHInxdCfdxIdIkBEYxcJJ9Ge4wblQobr0kljKdhq2NREyAOzEB7jw4OLSrzZG5hZF/HmRkZDpEFgoZGZkOkYVCRkamQ2ShkJGR6RBZKGRkZDpEFgoZGZkOkYVCRkamQ2ShkJGR6ZBb1jNTkqRKjOH1fm18gaqb8L43i9vJ3tvJVri97L1ZtoYLIfw6anTLCsXNQpKkkzfisnqrcDvZezvZCreXvV1tqzz0kJGR6RBZKGRkZDqkOwrFZ11twL/I7WTv7WQr3F72dqmt3W6OQkZG5l+nO/YoZGRk/kXueKGQJMlbkqTdkiRlm557XKethyRJRZIkfdyZNl5jQ4f2SpI0UJKkY5IknZckKV2SpIc72caJkiRlSZKUY0olee15R0mSNpnOp0qSFNGZ9l1jS0e2viRJUqbpOu6VJCm8K+xsZc917W3V7kFJkoQpK99N544XCuB1YK8QIgbYazpuj78AhzrFqva5EXubgMeEEP2BicDfJUlqm/r7JiBJki2wDJgExAKzJUmKvabZQqBGCNEL+AD4W2fYdi03aOsZjGkuBwBfA0s718qr3KC9SJLkDjwPpHaWbd1BKKYBX5jKXwDTrTWSJGkwEADs6iS72qNDe4UQl4QQ2aZyCVABdOg08ysxDMgRQuQKITTARow2t6b1Z/ga+I3UNTkNO7RVCLFfCNFkOkwBujLU141cWzD+oP0NUHWWYd1BKAKEEKWmchlGMbBAkiQb4D3g5c40rB06tLc1kiQNAxyAyzfbMBMhQGGr4yJTndU2QggdUAf4dIp17dhhwpqtrVkI/HRTLbo+HdorSdIgoKcQYntnGnZHxMyUJGkPEGjl1JutD4QQQpIka8s8TwM7hBBFnfHD9yvY2/I+QcBa4HEhhPU0XzI3hCRJc4EhwOiutqU9TD9o7wPzO/tv3xFCIYQY1945SZLKJUkKEkKUmm6sCivNhgOjJEl6GnADHCRJahRCXG8+oyvtRZIkD2A78KYQIuVm2NkOxUDPVsehpjprbYokSbIDPAFF55hn1Y4WrNmKJEnjMIr0aCGEupNss0ZH9roDccAB0w9aIPCDJEn3CyFubsh6IcQd/QDeBV43lV8HlnbQfj7w8a1sL8ahxl7ghS6wzw7IBSJNdqQB/a9p8wyw3FR+BNjcRdfyRmxNxDhsi+mq//m/Yu817Q9gnIi9+bZ19cXphIvvY7qpsoE9gLepfgiw0kr7rhaKDu0F5gJa4Gyrx8BOtHEycMl0g71pqvszcL+p7ARsAXKA40BUF17PjmzdA5S3uo4/dPH39br2XtO204RC9syUkZHpkO6w6iEjI/MfIguFjIxMh8hCISMj0yGyUMjIyHSILBQyMjIdIguFjIxMh8hCISMj0yGyUMjIyHTI/weFAicD9WVEmwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fix, ax = plt.subplots(1,1)\n",
"ax.set_aspect(\"equal\")\n",
"ax.set_xlim([-.5,.5])\n",
"ax.set_ylim([-.5,.5])\n",
"ax.plot(orbits[:,6::60],orbits[:,7::60])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
-3.664205888551948080e-04
7.319825676938241862e-03
-6.720660212985565319e-05
-4.451940140172100206e-04
1.366938569783342845e-04
1.113841358960815781e-05
3.570859982638058816e-01
-1.101735683190560172e-01
-4.246003048011558201e-02
1.921983878259755363e-01
1.627470754700197642e+00
1.153135053892231920e-01
3.745013392806120223e-01
6.240723868150592768e-01
-1.323735288432680880e-02
-1.008948002233915897e+00
6.056540210341178065e-01
6.651754105772547665e-02
6.489778509777003324e-01
7.545969441461461980e-01
-1.030728902624701533e-04
-7.716917064269775972e-01
6.523785058502404866e-01
-1.609427106253405731e-05
1.374707626818506778e+00
2.867095949934717614e-01
-2.795414324610703827e-02
-1.313395414670954497e-01
8.668450713601270063e-01
2.138413120653607710e-02
-2.472789988893362434e+00
-4.752512391238646217e+00
7.502383660041166136e-02
3.839039322450402691e-01
-1.815664756258509593e-01
-7.834119987025187776e-03
1.698487500660619487e+00
-9.909826200814681130e+00
1.046962230745202699e-01
3.017917659737404668e-01
5.375610208183433103e-02
-1.294773460018492327e-02
1.711798200594520125e+01
1.008988708557628478e+01
-1.842916696809539512e-01
-1.177781346127944467e-01
1.863128486686037866e-01
2.217843197649187539e-03
2.894283286410058054e+01
-7.634846177534076794e+00
-5.097915391345879677e-01
4.533290414756331599e-02
1.775343222299327339e-01
-4.700696384524625550e-03
1.173043539184388351e+01
-3.156283476648605912e+01
-1.572132791376496322e-02
1.750592885386928199e-01
2.508528944010257264e-02
-5.332210597456527035e-02
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.