1{ 2 "cells": [ 3 { 4 "cell_type": "markdown", 5 "id": "e6e99643-fa2f-40c5-b748-19ae6baae4f1", 6 "metadata": {}, 7 "source": [ 8 "# Inverting Kepler's equation in ODEs\n", 9 "\n", 10 "In the [previous tutorial](<./Comparing coordinate systems.ipynb>) we considered the numerical solution of the Stark problem, whose Hamiltonian, in Cartesian coordinates, reads\n", 11 "\n", 12 "$$\n", 13 "\\mathcal{H}_\\mathrm{cart}\\left(v_x, v_y, v_z, x, y, z \\right) = \\frac{1}{2}\\left( v_x^2+v_y^2+v_z^2 \\right) - \\frac{1}{\\sqrt{x^2+y^2+z^2}} - \\varepsilon z.\n", 14 "$$\n", 15 "\n", 16 "We saw how a reformulation of the Stark problem using Delaunay elements, instead of Cartesian coordinates, can help reducing the total number of timesteps per unit of integration time. The Delaunay elements are the Hamiltonian version of the [Keplerian orbital elements](https://en.wikipedia.org/wiki/Orbital_elements), to which they are related via the following equations (valid in adimensional units):\n", 17 "\n", 18 "$$\n", 19 "\\begin{aligned}\n", 20 "L & = \\sqrt{a}, & l & = M, \\\\\n", 21 "G & = \\sqrt{a\\left( 1 - e^2\\right)}, & g & = \\omega, \\\\\n", 22 "H & = \\sqrt{a\\left( 1 - e^2\\right)} \\cos i, & h & = \\Omega.\n", 23 "\\end{aligned}\n", 24 "$$\n", 25 "\n", 26 "The Hamiltonian of the Stark problem in Delaunay elements reads:\n", 27 "\n", 28 "$$\n", 29 "\\mathcal{H}_\\mathrm{Del} \\left( L, G, H, l, g, h \\right) = -\\frac{1}{2L^2}-\\varepsilon L\\sqrt{1-\\frac{H^2}{G^2}}\\left[ L\\left( \\cos E - \\sqrt{1-\\frac{G^2}{L^2}} \\right)\\sin g + G\\sin E \\cos g \\right].\n", 30 "$$\n", 31 "\n", 32 "In this Hamiltonian, $E$ is the eccentric anomaly, which is a function of $L$, $G$ and $l$ implicitly defined by [Kepler's equation](https://en.wikipedia.org/wiki/Kepler%27s_equation):\n", 33 "\n", 34 "$$\n", 35 "l = E - \\sqrt{1-\\frac{G^2}{L^2}} \\sin E.\n", 36 "$$\n", 37 "\n", 38 "Because Kepler's equation is trascendental, it cannot be inverted to yield an explicit expression for $E$ in terms of elementary functions. In the [previous tutorial](<./Comparing coordinate systems.ipynb>) we worked around this issue by leaving $E$ as an unspecified function in the Hamiltonian; we then augmented Hamilton's equations via the derivatives of $E$ (which can be formulated explicitly) and numerically solved the augmented system of differential equations.\n", 39 "\n", 40 "Here we will employ a different approach, and we will use instead the function ``kepE()`` provided by heyoka.py's expression system. As the name suggests, ``kepE()`` represents symbolically the inversion of Kepler's elliptic equation. That is, ``kepE(e, M)`` is the bivariate function $E\\left(e,M\\right)$ implicitly defined by Kepler's elliptic equation:\n", 41 "\n", 42 "$$\n", 43 "M = E - e \\sin E.\n", 44 "$$\n", 45 "\n", 46 "For the evaluation of ``kepE(e, M)`` heyoka.py uses an iterative numerical scheme. The high-order derivatives necessary to implement Taylor's integration method are built using automatic-differentiation techniques starting from the evaluation of ``kepE(e, M)``.\n", 47 "\n", 48 "Let us begin by writing down the Hamiltonian:" 49 ] 50 }, 51 { 52 "cell_type": "code", 53 "execution_count": 1, 54 "id": "c82ace35-08ea-4c28-b04c-81b468c5193b", 55 "metadata": {}, 56 "outputs": [], 57 "source": [ 58 "import heyoka as hy\n", 59 "\n", 60 "# Numerical value for the eps constant.\n", 61 "eps = 1e-3\n", 62 "\n", 63 "# Symbolic variables for the Delaunay arguments.\n", 64 "L, G, H, l, g, h = hy.make_vars(\"L\", \"G\", \"H\", \"l\", \"g\", \"h\")\n", 65 "\n", 66 "# Define E as the solution to Kepler's equation.\n", 67 "E = hy.kepE(hy.sqrt(1-G**2*L**-2), l)\n", 68 "\n", 69 "# The Hamiltonian.\n", 70 "Ham_del = -0.5*L**-2 - eps*L*hy.sqrt(1.-H**2*G**-2)*(L*(hy.cos(E)-hy.sqrt(1.-G**2*L**-2))*hy.sin(g)+G*hy.sin(E)*hy.cos(g))" 71 ] 72 }, 73 { 74 "cell_type": "markdown", 75 "id": "a6a91824-ca5b-4140-8485-2434733ad532", 76 "metadata": {}, 77 "source": [ 78 "Note how, in terms of Delaunay arguments, $e=\\sqrt{1-\\frac{G^2}{L^2}}$ and $M=l$.\n", 79 "\n", 80 "Unlike in the [previous tutorial](<./Comparing coordinate systems.ipynb>), and thanks to the fact that ``kepE()`` supports symbolic differentiation, we can now use the standard formulation of Hamilton's equations:" 81 ] 82 }, 83 { 84 "cell_type": "code", 85 "execution_count": 2, 86 "id": "a1b17f85-4efc-4ec0-8b20-fb1aa8a3bddb", 87 "metadata": {}, 88 "outputs": [], 89 "source": [ 90 "# Define the system of ODEs.\n", 91 "sys = [(L, -hy.diff(Ham_del, l)),\n", 92 " (G, -hy.diff(Ham_del, g)),\n", 93 " (H, -hy.diff(Ham_del, h)),\n", 94 " (l, hy.diff(Ham_del, L)),\n", 95 " (g, hy.diff(Ham_del, G)),\n", 96 " (h, hy.diff(Ham_del, H))]" 97 ] 98 }, 99 { 100 "cell_type": "markdown", 101 "id": "72fe3a74-e035-440c-b420-3646aabb4516", 102 "metadata": {}, 103 "source": [ 104 "We can now proceed to the creation of the integrator object. The initial conditions are taken from the [previous tutorial](<./Comparing coordinate systems.ipynb>):" 105 ] 106 }, 107 { 108 "cell_type": "code", 109 "execution_count": 3, 110 "id": "56aedf2a-767a-4333-a202-778a99c8430f", 111 "metadata": {}, 112 "outputs": [], 113 "source": [ 114 "ta = hy.taylor_adaptive(sys,\n", 115 " [1.0045488165591647, 0.9731906288081488, -0.9683287292736491,\n", 116 " 2.776991035843252, 4.314274521695855, 3.3415926535897924])" 117 ] 118 }, 119 { 120 "cell_type": "markdown", 121 "id": "5cfa5b10-aec7-41d6-bb2f-4a31ac51c6a3", 122 "metadata": {}, 123 "source": [ 124 "Let us now integrate up to $t=250$ on a time grid. As usual, we will provide a callback to periodically reduce the $l$ angle modulo $2\\pi$:" 125 ] 126 }, 127 { 128 "cell_type": "code", 129 "execution_count": 4, 130 "id": "b0c52e16-f32c-49d4-9bb3-0cd07e7ae1af", 131 "metadata": {}, 132 "outputs": [], 133 "source": [ 134 "import numpy as np\n", 135 "\n", 136 "# Callback to reduce l to the\n", 137 "# [-pi, pi] range.\n", 138 "def mod_cb_del(ta):\n", 139 " \n", 140 " l = ta.state[3]\n", 141 " if l < -np.pi or l > np.pi:\n", 142 " ta.state[3] = (l + np.pi) % (2 * np.pi) - np.pi\n", 143 " \n", 144 " return True\n", 145 "\n", 146 "# Define a time grid for the integration.\n", 147 "t_grid = np.linspace(0, 250, 1000)\n", 148 "\n", 149 "# Integrate.\n", 150 "_, _, _, _, out_del = ta.propagate_grid(t_grid, callback = mod_cb_del)" 151 ] 152 }, 153 { 154 "cell_type": "markdown", 155 "id": "26208487-5c97-4412-9f6a-9efad835a9c2", 156 "metadata": {}, 157 "source": [ 158 "Let us take a look at the time evolution of the Delaunay elements:" 159 ] 160 }, 161 { 162 "cell_type": "code", 163 "execution_count": 5, 164 "id": "77b9c512-ad7d-43ab-b11b-9a15b00d6d28", 165 "metadata": {}, 166 "outputs": [ 167 { 168 "data": { 169 "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAAGsCAYAAADqs/chAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABi3ElEQVR4nO3debwkWVnn/++Ty11rr65e6IWmWQVkLVsQBQUUBLRHxnFwRsVten4uiDouIM7AzOj8hEHE0dGxRRxGUGfYGVRWBfSnNFQ3TS80zdZN0xtdXbe2u+Z2fn9EROaJiJMZ93blzcgb9/N+vepVdePG8mTkvafOk885J8w5JwAAAADAcLWyAwAAAACAaUfiBAAAAAAFSJwAAAAAoACJEwAAAAAUIHECAAAAgAIkTgAAAABQYCoSJzP7BTO7xcxuNrO/MLO5smMCUF20OQAmhfYGqI7SEyczu1jSz0k66px7vKS6pJeUGxWAqqLNATAptDdAtZSeOMUakubNrCFpQdI9JccDoNpocwBMCu0NUBGNsgNwzt1tZq+XdKekNUkfcs59KLufmV0t6WpJWlxcfOpjHvOYyQYKIOW66657wDl3pOw4tmozbQ7tDTBdqtzeSLQ5wDQZ1d6Yc27S8aQDMDso6Z2S/qWkU5LeLukdzrm3Djvm6NGj7tixY5MJEECQmV3nnDtadhxbtdU2h/YGKN9uaW8k2hygbKPam2kYqvdcSbc7544759qS3iXpW0qOCUB10eYAmBTaG6BCpiFxulPS08xswcxM0nMk3VpyTACqizYHwKTQ3gAVUnri5Jy7VtI7JF0v6SZFMV1TalAAKos2B8Ck0N4A1VL64hCS5Jx7taRXlx0HgN2BNgfApNDeANVResUJAAAAAKYdiRMAAAAAFCBxAgAAAIACJE4AAAAAUIDECQAAAAAKkDgBAAAAQAESJwAAAAAoQOIEAAAAAAVInAAAAACgAIkTAAAAABQgcQIAAACAAiROAAAAAFBgKhInMztgZu8ws8+b2a1m9vSyYwJQXbQ5ACaF9gaojkbZAcR+V9IHnHPfb2YzkhbKDghApdHmAJgU2hugIkpPnMxsn6RnSvpRSXLOtSS1yowJQHXR5gCYFNoboFqmYajeFZKOS/pTM/uMmb3JzBbLDgpAZdHmAJgU2hugQqYhcWpIeoqkP3TOPVnSiqRXZHcys6vN7JiZHTt+/PikYwRQHYVtDu0NgDGhjwNUyDQkTndJuss5d2389TsUNTIpzrlrnHNHnXNHjxw5MtEAAVRKYZtDewNgTOjjABVSeuLknLtP0tfM7NHxpudI+lyJIQGoMNocAJNCewNUS+mLQ8ReJult8WozX5H0YyXHA6DaaHMATArtDVARU5E4OedukHS07DgA7A60OQAmhfYGqI7Sh+oBAAAAwLQjcQIAAACAAiROAAAAAFCAxAkAAAAACpA4AQAAAEABEicAAAAAKEDiBAAAAAAFSJwAAAAAoACJEwAAAAAUIHECAAAAgAIkTgAAAABQgMQJAAAAAApMTeJkZnUz+4yZvb/sWABUG+0NgEmizQGqYWoSJ0kvl3Rr2UEA2BVobwBMEm0OUAFTkTiZ2SWSXijpTWXHAqDaaG8ATBJtDlAdU5E4SXqjpF+R1Cs5DgDV90bR3gCYnDeKNgeohNITJzN7kaT7nXPXFex3tZkdM7Njx48fn1B0AKqE9gbAJNHmANVSeuIk6RmSvtfM7pD0l5KebWZvze7knLvGOXfUOXf0yJEjk44RQDXQ3gCYJNocoEJKT5ycc690zl3inLtc0ksk/a1z7odKDgtABdHeAJgk2hygWkpPnAAAAABg2jXKDsDnnPuYpI+VHAaAXYD2BsAk0eYAOx8VJwAAAAAoQOIEAAAAAAVInAAAAACgAIkTAAAAABQgcQIAAACAAiROAAAAAFCAxAkAAAAACpA4AQAAAEABEicAAAAAKEDiBAAAAAAFSJwAAAAAoACJEwAAAAAUKD1xMrNLzezvzOxWM7vFzF5edkwAqos2B8Ck0N4A1dIoOwBJHUn/zjl3vZntlXSdmX3YOfe5sgMDUEm0OQAmhfYGqJDSK07OuXudc9fH/z4r6VZJF5cbFYCqos0BMCm0N0C1lJ44+czscklPlnRtyaEA2AVocwBMCu0NsPNNTeJkZnskvVPSzzvnzgS+f7WZHTOzY8ePH598gAAqZVSbQ3sDYJzo4wDVMBWJk5k1FTUob3POvSu0j3PuGufcUefc0SNHjkw2QACVUtTm0N4AGBf6OEB1lJ44mZlJ+hNJtzrn3lB2PACqjTYHwKTQ3gDVUnriJOkZkn5Y0rPN7Ib4zwvKDgpAZdHmAJgU2hugQkpfjtw59w+SrOw4AOwOtDkAJoX2BqiWaag4AQAAAMBUI3ECAAAAgAIkTgAAAABQgMQJAAAAAAqUvjjEdvnI576ua28/UXYYQCmce/DH/tLzHq25Zn18wewC6+2ubrkn90zLITb/5mzlfdzKW76l825h583uubXXtT03YdvuV8nv73a8X1vdeSv34OFH9uihhxe3EgkAlKayidNnvnZSb7v2zrLDAErzYJdx+rnnPpLEaYvuPb2uf/6H/1h2GMCO88rvfoz+7bMeXnYYALAplU2cfvl5j9EvP+8xZYcBYBe4YN+s3vLjV256/60ktbaFnW0LZ97aebdgkztPQ6y2hRNvXwxb2HmTZ56K+7XJ/S7aP7eFCACgXJVNnABgUhZmGnrWo46UHQYAANhGLA4BAAAAAAVInAAAAACgAIkTAAAAABQgcQIAAACAAiROAAAAAFBgKhInM3u+md1mZl8ys1eUHQ+A8pjZITN7yDZfgzYHgMzsUbaVddYf3DVob4CKKH05cjOrS/rvkr5T0l2SPm1m73POfa7cyACU5PWSvijp/5UkM/tHRW3D9ZL+zDl397mcnDYHgOddki41sy9IuknSjcnfzrnj53py2hugWs6p4mRm/98YYrhS0pecc19xzrUk/aWkq8ZwXgA701Ml/Zb39V5JfyLpPEmvHMP5aXMASJKcc4+XdL6kn5b0PZIeLunXJN1kZveN4RK0N0CFnGvFaRzDaS6W9DXv67skfXN2JzO7WtLVknTZZZeN4bIAptSGc855X/+tc+6DZvYhSf80hvMXtjm0N8Du4ZzbUFQJWnbOvSzZbmYHx3B6+jhAhRRWnMzs98zsajN7upntzXzbBQ/amtDY4tx5nXPXOOeOOueOHjlyZAyXBTCl1s3sockXzrmXx387Sc0xnL+wzaG9AXalbDtwcgznpI8DVMhmKk43SXqCpH8t6fFmdibedpOiITTn6i5Jl3pfXyLpnjGcF8DO9JuS3mNmP+ic+3yy0cwu0njmZdLmAJAkmdnvS/qMojmU27FIBO0NUCGFnRDn3DX+12Z2iaJE6hslfXAMMXxa0iPN7GGS7pb0Ekn/agznBbADxcPy9kn6OzO7QdLN8bdeLOnXx3AJ2hwAiRslPVnSj0jaa2afk3SLpM9J+pxz7n+f4/lpb4AK2fKnt865uxR9gvLX4wjAOdcxs59VlITVJb3ZOXfLOM4NYGdyzr3dzP5K0gskPU7SmqQXO+c+O4Zz0+YAkFT44fCLJJ1T4kR7A1RL6cuRS5Jz7q81pkQMQDU451YlvSP+M+5z0+YAyBn3h8PxOWlvgIqYigfgAgAAAMA0I3ECAAAAgAIkTgAAAABQYCrmOAEAAGBz/uP/vUWfu+dM2WEAO8pjH7JPr/6ex53TOag4AQAAAEABKk4AAAA7yLl+ag7gwaHiBAAAAAAFSJwAAAAAoACJEwAAAAAUIHECAAAAgAIkTgAAAABQgMQJAAAAAAqUmjiZ2X81s8+b2Y1m9m4zO1BmPACqjTYHwKTQ3gDVU3bF6cOSHu+ce4KkL0h6ZcnxAKg22hwAk0J7A1RMqYmTc+5DzrlO/OUnJV1SZjwAqo02B8Ck0N4A1VN2xcn345L+Ztg3zexqMztmZseOHz8+wbAAVNTQNof2BsCY0ccBKqCx3Rcws49IujDwrVc5594b7/MqSR1Jbxt2HufcNZKukaSjR4+6bQgVQAWMo82hvQGwGfRxgN1l2xMn59xzR33fzF4q6UWSnuOco7EAcE5ocwBMCu0NsLtse+I0ipk9X9KvSnqWc261zFgAVB9tDoBJob0BqqfsOU6/L2mvpA+b2Q1m9j9KjgdAtdHmAJgU2hugYkqtODnnHlHm9QHsLrQ5ACaF9gaonrIrTgAAAAAw9UicAAAAAKAAiRMAAAAAFCBxAgAAAIACJE4AAAAAUIDECQAAAAAKkDgBAAAAQAESJwAAAAAoQOIEAAAAAAVInAAAAACgAIkTAAAAABQgcQIAAACAAlOROJnZL5mZM7Pzyo4FQPXR5gCYFNoboDpKT5zM7FJJ3ynpzrJjAVB9tDkAJoX2BqiW0hMnSb8j6VckubIDAbAr0OYAmBTaG6BCSk2czOx7Jd3tnPvsJva92syOmdmx48ePTyA6AFWz2TaH9gbAuaKPA1RPY7svYGYfkXRh4FuvkvRrkr5rM+dxzl0j6RpJOnr0KJ/cAAgaR5tDewNgM+jjALvLtidOzrnnhrab2TdKepikz5qZJF0i6Xozu9I5d992xwWgmmhzAEwK7Q2wu2x74jSMc+4mSecnX5vZHZKOOuceKCsmANVFmwNgUmhvgGqahsUhAAAAAGCqlVZxynLOXV52DAB2D9ocAJNCewNUAxUnAAAAAChA4gQAAAAABUicAAAAAKAAiRMAAAAAFCBxAgAAAIACJE4AAAAAUIDECQAAAAAKkDgBAAAAQAESJwAAAAAoQOIEAAAAAAVInAAAAACgAIkTAAAAABQoPXEys5eZ2W1mdouZva7seABUG20OgEmhvQGqpVHmxc3sOyRdJekJzrkNMzu/zHgAVBttDoBJob0BqqfsitNPSfot59yGJDnn7i85HgDVRpsDYFJob4CKKTtxepSkbzOza83s42b2TcN2NLOrzeyYmR07fvz4BEMEUCGbanNobwCMAX0coGK2faiemX1E0oWBb70qvv5BSU+T9E2S/o+ZXeGcc9mdnXPXSLpGko4ePZr7PgBI42lzaG8AbAZ9HGB32fbEyTn33GHfM7OfkvSuuBH5lJn1JJ0niY9bADwotDkAJoX2Bthdyh6q9x5Jz5YkM3uUpBlJD5QZEIBKe49ocwBMxntEewNUSqmr6kl6s6Q3m9nNklqSXhoqYQPAmNDmAJgU2hugYkpNnJxzLUk/VGYMAHYP2hwAk0J7A1RP2UP1AAAAAGDqkTgBAAAAQAESJwAAAAAoQOIEAAAAAAVInAAAAACgAIkTAAAAABQgcQIAAACAAiROAAAAAFCAxAkAAAAACpA4AQAAAEABEicAAAAAKEDiBAAAAAAFSk2czOxJZvZJM7vBzI6Z2ZVlxgOg2mhzAEwK7Q1QPWVXnF4n6T86554k6T/EXwPAdqHNATAptDdAxZSdODlJ++J/75d0T4mxAKg+2hwAk0J7A1SMOefKu7jZN0j6oCRTlMR9i3Puq0P2vVrS1fGXj5Z02yYucZ6kB8YQ6nYixvEgxvHYSowPdc4d2c5gxm2zbc6DbG+k6r3HZSHG8ahSjJVtb+J96eOUhxjHo0oxDm1vtj1xMrOPSLow8K1XSXqOpI87595pZj8g6Wrn3HPHeO1jzrmj4zrfdiDG8SDG8dgJMRahzRmNGMeDGMdjJ8Q4Cu3NaMQ4HsQ4HuOIsTGuYIYZ1UiY2f+S9PL4y7dLetN2xwOg2mhzAEwK7Q2wu5Q9x+keSc+K//1sSV8sMRYA1UebA2BSaG+Aitn2ilOBfyPpd82sIWldg/G943LNmM+3HYhxPIhxPHZCjOeCNocYx4UYx2MnxPhg0d4Q47gQ43icc4ylLg4BAAAAADtB2UP1AAAAAGDqkTgBAAAAQIHKJk5m9nwzu83MvmRmryg7noSZ3WFmN5nZDWZ2LN52yMw+bGZfjP8+OOGY3mxm95vZzd62oTGZ2Svj+3qbmT2vxBhfY2Z3x/fyBjN7QVkxmtmlZvZ3Znarmd1iZi+Pt0/NfRwR49Tcx52K9mZLMdHejCdG2pxdivZmy3HR5px7fLQ3Cedc5f5Iqkv6sqQrJM1I+qykx5YdVxzbHZLOy2x7naRXxP9+haTXTjimZ0p6iqSbi2KS9Nj4fs5Kelh8n+slxfgaSb8U2HfiMUq6SNJT4n/vlfSFOI6puY8jYpya+7gT/9DebDkm2pvxxEibswv/0N48qLhoc849Ptqb+E9VK05XSvqSc+4rzrmWpL+UdFXJMY1ylaS3xP9+i6R/NsmLO+c+IWlpkzFdJekvnXMbzrnbJX1J0f0uI8ZhJh6jc+5e59z18b/PSrpV0sWaovs4IsZhSnmvdyDamy2gvRkP2pxdi/Zmi2hzzh3tzUBVE6eLJX3N+/oujb55k+QkfcjMrjOzZGnSC5xz90rRGy/p/NKiGxgW07Td2581sxvjMndSIi41RjO7XNKTJV2rKb2PmRilKbyPO8g03yfam/Gayt8T2pxdZZrv0U5pb6Qp/T0JmLrfk93e3lQ1cbLAtmlZd/0ZzrmnSPpuST9jZs8sO6AtmqZ7+4eSHi7pSZLulfTb8fbSYjSzPZLeKennnXNnRu0a2FZWjFN3H3eYab5PtDfjM5W/J7Q5u84036Od3t5I03V/p+73hPamuonTXZIu9b6+RNETvEvnnLsn/vt+Se9WVBb8upldJEnx3/eXF2HfsJim5t46577unOs653qS/liDEmspMZpZU9Ev69ucc++KN0/VfQzFOG33cQea2vtEezM+0/h7QpuzK03tPdpB7Y00Zb8nIdP2e0J7E6lq4vRpSY80s4eZ2Yykl0h6X8kxycwWzWxv8m9J3yXpZkWxvTTe7aWS3ltOhCnDYnqfpJeY2ayZPUzSIyV9qoT4kl/SxPcpupdSCTGamUn6E0m3Oufe4H1rau7jsBin6T7uULQ3525qfk+GmbbfE9qcXYv2Zjym5vdkmGn6PaG98RStHrFT/0h6gaIVNb4s6VVlxxPHdIWiFTw+K+mWJC5JhyV9VNIX478PTTiuv1BUvmwrysB/YlRMkl4V39fbJH13iTH+maSbJN0Y/wJcVFaMkr5VUYn3Rkk3xH9eME33cUSMU3Mfd+of2pstxUV7M54YaXN26R/amy3HRptz7vHR3sR/LD4QAAAAADBEVYfqAQAAAMDYkDgBAAAAQAESJwAAAAAoQOIEAAAAAAVInAAAAACgAIkTzpmZHTazG+I/95nZ3fG/l83sD8qOD0C10OYAmBTaG/hYjhxjZWavkbTsnHt92bEAqD7aHACTQnsDKk7YNmb27Wb2/vjfrzGzt5jZh8zsDjN7sZm9zsxuMrMPmFkz3u+pZvZxM7vOzD6YeeIzAAxFmwNgUmhvdicSJ0zSwyW9UNJVkt4q6e+cc98oaU3SC+OG5fckfb9z7qmS3izpN8sKFsCOR5sDYFJob3aBRtkBYFf5G+dc28xuklSX9IF4+02SLpf0aEmPl/RhM1O8z70lxAmgGmhzAEwK7c0uQOKESdqQJOdcz8zabjDBrqfoZ9Ek3eKce3pZAQKoFNocAJNCe7MLMFQP0+Q2SUfM7OmSZGZNM3tcyTEBqC7aHACTQntTASROmBrOuZak75f0WjP7rKQbJH1LqUEBqCzaHACTQntTDSxHDgAAAAAFqDgBAAAAQAESJwAAAAAoQOIEAAAAAAVInAAAAACgAIkTAAAAABQgcQIAAACAAiROAAAAAFCAxAkAAAAACkxF4mRmv2Bmt5jZzWb2F2Y2V3ZMAKqLNgfApNDeANVReuJkZhdL+jlJR51zj5dUl/SScqMCUFW0OQAmhfYGqJbSE6dYQ9K8mTUkLUi6p+R4AFQbbQ6ASaG9ASqiUXYAzrm7zez1ku6UtCbpQ865D2X3M7OrJV0tSYuLi099zGMeM9lAAaRcd911DzjnjpQdx1Ztps2hvQGmS5XbG4k2B5gmo9obc85NOp50AGYHJb1T0r+UdErS2yW9wzn31mHHHD161B07dmwyAQIIMrPrnHNHy45jq7ba5tDeAOXbLe2NRJsDlG1UezMNQ/WeK+l259xx51xb0rskfUvJMQGoLtocAJNCewNUyDQkTndKepqZLZiZSXqOpFtLjglAddHmAJgU2hugQkpPnJxz10p6h6TrJd2kKKZrSg0KQGXR5gCYFNoboFpKXxxCkpxzr5b06rLjALA70OYAmBTaG6A6Sq84AQAAAMC0I3ECAAAAgAIkTgAAAABQgMQJAAAAAAqQOAEAAABAARInAAAAAChA4gQAAAAABUicAAAAAKAAiRMAAAAAFCBxAgAAAIACJE4AAAAAUIDECQAAAAAKTEXiZGYHzOwdZvZ5M7vVzJ5edkwAqos2B8Ck0N4A1dEoO4DY70r6gHPu+81sRtJC2QEBqDTaHACTQnsDVETpiZOZ7ZP0TEk/KknOuZakVpkxAagu2hwAk0J7A1TLNAzVu0LScUl/amafMbM3mdlidiczu9rMjpnZsePHj08+SgBVUdjm0N4AGBP6OECFTEPi1JD0FEl/6Jx7sqQVSa/I7uScu8Y5d9Q5d/TIkSOTjhFAdRS2ObQ3AMaEPg5QIdOQON0l6S7n3LXx1+9Q1MgAwHagzQEwKbQ3QIWUnjg55+6T9DUze3S86TmSPldiSAAqjDYHwKTQ3gDVUvriELGXSXpbvNrMVyT9WMnxAKg22hwAk0J7A1TEVCROzrkbJB0tOw4AuwNtDoBJob0BqqP0oXoAAAAAMO1InAAAAACgAIkTAAAAABQgcQIAAACAAiROAAAAAFCAxAkAAAAACpA4AQAAAEABEicAAAAAKEDiBAAAAAAFSJwAAAAAoACJEwAAAAAUaJQdAADsdPeeXtM//4N/1JG9szqyd1b75pra6PS00elp33xDhxZmdHBxRr2eU7vb0wX75/SQ/fPaO9eQk9Tu9iRJhxej4xt103q7K+ek/fNNzTXr5b5AAAAwPYmTmdUlHZN0t3PuRWXHA6C6tqO9edrDD+uB5ZbuPrWuz6+f1VyzrkbNdPbejk6utrTa6sbXlpzb2rkXZ+raO9fUWrurTren+ZmGFmfrWphpaL4ZDRxYmGlo71z0Z65ZV6vTU7Ne08JsXYszDS3M1NWs11SrmfbNNbRvvql9cw1tdHrq9pzmm9H5FmbqmmvWZSbNz9S1Z6ahWs3GcYuAXYs+DlANU5M4SXq5pFsl7Ss7EACVN9b25qL983rDDzxp5D7r7a5qZqrXTPefXde9p9e1vN7pb3NyOrHc0vGzG+o5p9lmXTWTTq22tbTS0tn1tuabdTXqNa22ulptdbSy0dV6O0rIVlod3XdmXWfX21pv9zTTqKnd7Wm11VWr03vQr81M2jvbUL1mWm111azXNNesaa4ZJVgz9ZpmGjXtm2/qwHxTCzN1bXR6anV72jPT0J65hvbMDiprs42aFmbqmm/WNdusq9dzUYI3U9d8vL3Tc1Hi1kySw3o/4ZybqUUJXryt03Nq1ExmJHeYavRxgAqYisTJzC6R9EJJvynpF0sOB0CFldXe+MPtLto/r4v2z0/q0v0Eqttz6nR7OrPe0em1tpY3Opqp19Som9biZGy11dV6uycnp7VWV2fW2jq91lbPRRWoTtdpvdPVequrtXaUlLW6PZ1ebemrJ1a02upGCV7NtNLqaHm9o5W42tasm9rdLZbbNsFMmm3U+olczUyzzZr2zEZJW71mWosTzNlGTbONumYb0es2M+2ZaWhxNqq2tXs9tTtOM16CV6tZXMGzQcLYqKnb66lm0bb5ZnTOTs+p51y0X6Ouubgi2Ok5zTZqmp+JttfMtNHtqlmL4p5tRNVAF5cjSQSrgz4OUB1TkThJeqOkX5G0d9gOZna1pKsl6bLLLptMVACq6I3aZe1Ns17T/vnBWkDnT/gz727PySTVaqZuz2mt3dVaK6qW1WvWT+xW422NmqnnoipdUl2TomQiOjaqtjlJzZqp1Y3mk220u9ro9NRzThudnpbXOzq70VGr09Oe2ei/u41OT6fW2tpod6Mkp+e02upqeaOjlVZHzXpNM/VaPyGcpEbN+tW22UZN83GC1u46OeeihK9Z02yjJuei+zrTqGm2WddcI3p/W92eGjWLXkejpkat1j/3wkxdC7N11c200en1k7kkEex0e2p3XZxcRueVosS7Wa/1t8uiRLJmiq7fqPeTYhdXS2cb8X3sRkNBk+Nn4tg7veg9uezQgg7vmZ3ofS7BG7XL2hygqkpPnMzsRZLud85dZ2bfPmw/59w1kq6RpKNHj47/I0sAlUd7U466N0eqXrN+JWjaOOdSlZ5Ot6e1dle9ntRsREnNelyRa3W7qtdq6jmn9XY3/hMlLbWa9b9ea3dlil53qxN9vd7uqufUT1bW2z2tt7tqd3tq1Gty/XP21OpEwy4lqdXpab3T1Ua7p1pNqsUJzHonOt4k7ZltqNN1anV6WtnoqBVX+Lq9QXLa7TnNNWv9BGq93Y0TrqgK145jmoRfff5j9FPf/vCJXKsMtDlAtUzD/1zPkPS9ZvYCSXOS9pnZW51zP1RyXNiBHlje0Kvfd4t+68XfqL1zzbLDkRR9WmuSGvXxrf5//Z0nNduo6XEP2T+2c+4StDcYKjs8rlGvaW/m93bflLQr2805p1a3J5P1K2FJ0pZUv5yiyl5SnWvWajJTvKJkV+1uNP+sUTe1O06tbnS8xfP6VjY6uvy8xbJf6najzQEqpPTEyTn3SkmvlKT405hfokGZPksrLb3wv/29/vhHjurxF09vZ/1Tty/pr268Vz/xrQ/TUy47WHY4kqRffvtn1e45/fd/9ZSxnfPfv+dmXbR/Xm966dH+tlOrLa23e7pw/9zYrlM1tDfA5piZZhuDeXkzNYsqXzQvW0KbA1QLD8CVdMcDK/rDj315U/v+4v++QX//xeOF+/2fT39Nz3/jJ841tEKdbk/3nV7f9uvc/sCy7j29rtsfWHlQx//5tXfqB6/55JijyltaaUna+nLPo9zxwIr+6sZ7H/TxX/j68tjfoweWN9TppYfS/MZf3ap/+9brxnodAAAARKYqcXLOfayM5xu854a79doPfL4/AXmYtVZX7/rM3frU7UuF57zx7lP6wtfPprbddt9Z/eRbPq2NTrfw+B/4o3/Sm//h9sL93nn9XfqO13+svyTxVp1ebRe+bkk6sdx6UOdPXPfVk7rp7tMP+vh/+Uf/pHded1fhfidXkjjHlzn9z3+8Q7/yjs+mtm10urrxrlObOv7k6rnduyznnE6utHPJ4V0nV3V2vZ3a9tZPflXP+e2PjfX6VVFWewNgd6LNAXa+qUqcynJyk1WKpdXNVzNOrrRz2/7hSw/oI7ferwcKkhDnnK7/6kl9+fhy4XXuOLGqtXglqcTn7jmjo7/xER0/u1F4/I/+z0/pv/z1rYX7JZ3/7Eu//8zmKinnkjy0Oj1de/uSPn/fmdT2P/r4l/XeG+5ObQu9RzfddVrP+e2P5ZKKzTqx0lI386a/+/q79X1/8I86vTb6nM45La20+ksMJ+4+tbap9ydkpRVN5M6+FydX2rk36JZ7zujOpdXMttP6kTd/6pye7QMAKM81n/iy/tP//Zz+x8e/rHded5fedf1d+otP3al3f+YuffwLx3Xz3af1lePL+tw9Z3T3qTV1e6w3AYxD6XOctsvyRkf/+k3X6rde/I36hotGr727tBp1foualX6CtYlqxtJKK9CxTTr1o48/u9FRZ5ONXL/C4u1+8z2n9cDyhr5+Zl1H9o5e5vWrJ1b1kMzzZP7LX9+qkyst/dd/8cT+tqU4EfRj//x9Z/T8N/69/urnvrVwkYITgeQhWcmpaLL1qSFJ15998qt6/EP266onXdzfNniPBm742kl9+fiK7j+70V8w4ux6W7/+npv16u95nA4tzoy8/smVVi5Zvuf0uro9F1cPh8efJLXZd/Nn//x6XXZoQb/7kif3tyX/sfkrkIUsDUm8l1Zb2ptZqWxQgRv4py+f0Ce+cFxLK63+fKjVVke/+9Ev6hee+6jU84YAANPnM3ee0ie+cLz/jLQizbrp4MKMZho19XpOPSftnWvowv1z2j/fVKfr1O72tGeuocOLszq8Z0Zm0no7emj13vhB1ouzDTXrppl6Xfvnm9o339C+uaacJJO0ONvorwIJVFFlE6d/+OJxffZrp/SGD39Bf/wjgwn0b/r7r+jbH31Ejzh/8DiFUELzsdvu19/cdJ9e+/1P6G8LzZ9Zb3f1nb/zcb36RY/Tcx97weCcq/nO9mYrVqHOvyTdf3Zd++ebqQm7S4GOcaizHNLtuSjOzJWu++rJXDUiVDG644HVXAwrGx298L/9vX7rnz9BT7vicCqm7Ot57Qc+r+vvPKX3/swzRsY57L5F50xvPBF4j04E7sdnv3Za773hHn3fky/Wtz/6/P72P7/2Tn3X4y7Qed5zRU6MSIKLcuhkiGM29rtPrulwJmH7z+//nG5/YEVv+fErR55zcD8GJ42G77VySzwvBX4OTwQ+APinL5/QH338K3rOYy7QlQ87NPpFAQBK9Yc/9FRJ0YfEx89uyBQtb7/R6enE8oYeWG5ptdXRfLOuU2tt3bm0qqXlltrxQ5trJp1ea+vrZzZ018k1NevRs7+Wj3f0wNmNTSdkITON6OHT88266jXT3rmGDiw0dWB+pr8s/1yzpsXZQTJWr5lmGzXtn2/qwMKMFmfr/X7Iwkz0cOroT6N/fhI0lKGyiVNSITm0MOicrrW6+o2/ulVn1zv6he8cJE6hjvVHb71f777h7nDi5O13/OyGvra0pjtOpBdNCJ0zlNB84Ob79NoPfF4f/oVn9perDiVovZ7T837nE/qZ73iEfvLbrhicczXfCQ4lGs45ff3MRmrFtdNr+XkyyfVzHfBATCcD17n71JruOLGqLx9fzidOmWvd/sCKTiynh6u9/8Z79Ka/v13v8ZKp0H1fb3e10urmk6lAghcairkUGHp4z6k1/dq7b1LPOf3Q0x6aO963FLjOfafX9ak7lvS9T3xILh7/Os65YGL9+fvO5BLhv7npXv3VTffq970V+ULxnFmPqpTZql4oYQ3ej01WQwEA0yP0TLSHjWGJ92Te9GycjC1vdLSy0dHZ9U482qKn02ttnV5r68xaWzWL/p9b2YgeOr2y0ek/M2x5vaOTqy3de+qMei56aPN6O3rO2PJGJzXVYCv2zEYJmXPRg5/nmjXtnY2qYPPNev/ZZIuzdS3ORAlazUz1WpSMDZK26P/D+Zm69sw2tDDT0Gyz1l92P9lvrlnrP9NtvlnPPb4Au0NlE6ekw3pgcTCMKtRZlsIVnqXVVm7HUPIQqvgkn/4n/05+uULVkBvvOqXbH1jRRqeXS5z8AM6st3VytZ1LDIIJTaCi8MFb7tPP/cUN+tSrnqMDcTI5bAW6pZWWFmfTw7VC5wy99tC2VqensxsdzTXTnw6FkqnrvnpSN3ztVOq+heI8NWR45cnAkMJkKKa/91KSsAXey16mkhP8WVjO/8z8+afu1O/97Rf1gsdfmH8vvXMub3TU7uYTlNC8uE988QF99Nb709ce+Z5nzrmaHyIZSkSHzWEDAOw+/pDtuWZdc816aiTGOHW6PXWd03q7pzNrbZ1abevsRluzcUUpeXDzejv6e6Pd1Zk4GTu92paZaaZh2mj3dGa9rTNrHT2w3IofMN3RXSc7WtnoamWjI6dotM3ag1xQK2EmLTTrmp9pyEyqmfrJ2cJMXWbR/9Fzzbrmm3XNNWuaa9ZVq5kWmnXtm29q71yUyG10oodp+/s169HDqPfONrR3rqmFmSgR7PacZhvRPnONumabNdXM1KwbidyEVDZxSubFHPQqTicDndh+x1j5jmh2GFiomhE6dtgcpVClIHjOkRWBIfumtrVz+97+wKpa3ehToyRxClWMOt3oU6RLXHreU6iCFr52KMkZPicnK3iPRiRtm7kf4YpTkniNPudqq6tWp6dmPd0ghe7dieUNORdOSNKvJ5z0nVhp6eBCM7Pv8J/Doipjr+d0crUdTKaifQffCSX1AABst0a9poak2UY0b+rSCYwW7/acVlpRZaznovlZUYIWJVnrna5M6lfGVlodrbejBKfd7UX7bnTiIY1O3Z7TaitKzlY24tEwJp1aa+u+0+taa0eJX8+5fiI4Tn7iVq+Z2t2eGjXTXLOumTjRatSiB08vxhW0+WZN3V7Ul5hr1rU4EyeC8f2ZbdQ0P1PX/EyUyPV6Ts36YNtso6Zuz6lmlkoQE3PNaHhlkoT3nNNMvdb/YHmnqmzidDLuGO+f9ypOgU510jHOfmMpUA1ZCnTgQ516v/PvXPTJRBTTiAQgFfvwYXG+bs/p1Fq+whKqHoTOGXo9g/OlrxWMKTC8K5hMDZmjtLTc6idxg30H10/uWygRDFXl1lrd/qdI4dcZij2wX+j1bOZnIXiP80nS0DlKqy0dyCROo+Yo+cLD99rq9pyyH0KNmhe3mYVPAADYyeo10765ZuHiVNul3e3p7HonGg7YrKvbdVrvdPv9mG7PRR92r3d0Zr2t1VZXs42a6rWosrbe6Wq93dNGp6teL6rWJYlgpxclKJ2e03q8QNV6fM5O1+n+s+tafSBK3pLFqDY6Xa1sDPpQNZO2ayHGes00U69pplGTmTRTr2khTtqa9WgOXFRFrGk2/lOvmepJgjYTJWmdnlO311OzXtNsI0oQZ+LKW5IsLszUVTdTq9vTfLOuC/bN6Vsfed45xV/ZxCnJeuter7FoaFk2oRn2SX1oeNew4XvJ5uTT/+y+waFlgSFboQ74qDlKW9m2mSFfwdcZStACVb3QsLZkjtL++fycnNw5R1T6NGS/UEJTdHwwoRgxRyn6d/749M/RRmC/jdw5z8TjxkPD6jYzRym04ENh0hdI8AAAwPZq1mu5FX33j1ihd1J68YetZlHVaq3d1XorSr6SStZaPGSy1YmqWt2e03qnp7V4OGUyTHG93e3vK0UJU7vT00YnSvhanZ56ToMKXqurTq+nmXpNTtJGp6dWp6uz6x31XFTVW28Pkst6raZGHFMrPm+rO3q+3OWHF/SxX/6Oc7pHlU2cfurbH6G3fvLOwo5kKslxyd8uODfkRCABGDVcLTmXZP0Jldl9TwQ60aOqWEUJWvp1hqtgo2I/ETg2+WRkWJypRDKQaCwFEs5hc5TSsafnOPmK7oe89zI0DC10fChJClV3/DlK6fsRqoyFhgSGEughc5RWAvOrtvh6fH6VMnWdwHsEAAB2j5r3OJRmvaZmvVZaVe7Bci6q1vnVu5l6LR4uee7Pr6xs4pS89eGOcX5uiC/pGGcfpxOaGzJqjlJ0rfx10scHOtuh+TMjhob5+yZzlPxrDzvnqPlEqW1DYp/UHKWRQ+gC+/n8oZhFQyRHJyT5Y7NGDaErSpJCyWUvXi4++1ynzc5x8qtIyWIbfpWyaDgjAADATmJmmm3UNduo68A2nH9nz9AaIRmht9mOsb/vsMn7o+baDPukP9meDM/yz9vtucHCCUVzcgL7nVj2Y4++4VcTCoeRhTrbm0wEJX8Inn/O4VUXX6jzv6U5SgVJbD9hDVQUU+csmBu2uYpievieL3g/A7GHhjOeWW8HxxgX/RwP9vN+5vo/h8OHpmbPCQAAgIHqJk5xzaloKFaoEzzsOUijV0jLD8XyhbadWQt3jAuHoSX7+fG49LGpjSrqbG9+SGAimaOU3zeZ0xNISAIVtKGvJ7A9XIELx5nsG0piU3OUvOMHD6sdHXvonGf8oZhFQ+hCSc4mV1dMbU9tGz1Xrv+hQGGVkswJAAAgpPTEycwuNbO/M7NbzewWM3v5eM4b/e13BAcd48F+4SpF+qGsUqZj7G0vGgKXXD+70p5/bDbOpUA1Jji8K9ABPxG4zrA5SqMWCUhf20v64s2nVsOVrVC1LpTkFM/ZGlRyRq5mWHDfQ/do6HOUgsnL5pLg0LX9c4buR+j4ouGIw+YohX5mQ6/Hr1ImQuersu1qcwAgi/YGqJZpmOPUkfTvnHPXm9leSdeZ2Yedc587l5OGHgNW1DEeDGfKzw0JJT7RvqOrB6EkSYFOfejT/9A5hy344ALb+vsVPEcpPARuyOsJxe6fc4vzhIoqRmvxUpq52DeZNIaSsZOB9zzad3RCMupnIfSep1dSLKjqjRom6O2YWkkx9QFAYF5bYJGS0AcFw362K2xb2hwACKC9ASqk9IqTc+5e59z18b/PSrpV0sXjO//g3+EVzgKf1Ac6kicCHXB/jlLR4gHBznYg+RjaMR7R2U6dM9gxTieCymwPPgi1IPZQMubPUfKdy7Oqhs1RGvWsKn976PgTqfc8P0cpVEFLxRSKczl/7eQ5SsOOL6ygFQ0tLRiKGVykJDA09cSQJLiqtrvNAYAE7Q1QLaUnTj4zu1zSkyVdG/je1WZ2zMyOHT9+fBMni/7q5yD+HKUhHeN+NaVgtbj+vJaCOUr+vqHO6clAJ3ZYJacwoRk1JDBwzo1OV8sbnZGxJ0LHnwjEHkpc/O3hTn044RwVu3NuyPOvWmoMWYFu2LbgHKUhqy72YwosyrHZ50JJBQ/fDcyVCyVTyXMSpOFzlIJzvjb5c7hbDGtzttzeAECBsfZxAJRiahInM9sj6Z2Sft45dyb7fefcNc65o865o0eOHCk+Xz9zijqCRXOU/G+EOpejVkfzLjN036IFAUZVBPw5Sqnjg4lGvroU2i81RylwzmxCMtuopY4P3qMhq81t9eG70TmHx77S6gYfcnZypd1/oNzgtRcMxQzGPiymeN/A+x5KPpIkZ2Gm3t/Pn6O0ldUZE8kcpQPzzX48w1ZSDD0QeLNVyt1gVJuz1fYGAEYZdx8HQDmmInEys6aiBuVtzrl3jeec0d9Fn6qfXGn1n9c0spoSb2vWLXVsf7/47+TT/7lmnGiM6LAWPjMptG3YUL1kW2DZ883Or/Kvla0YHc484bqoOpRI5igliVfo+NDrycaTin15+P3oJ07J8edQlUueo9T/WXKB4wOxD96LKKE5tDjTv/ap1VZwHtHJ1dD8quH3w3/ieLLf3rn0lMVRCWvqnEMqY1W2HW0OAITQ3gDVUXriZGYm6U8k3eqce8PYzpv5OrSSWdIxHlQpou3DOsazjZoWZhojKznJHKXDi7O5cyZJ16CaMrwaMtuoecfmhxMm2/fMZjrLgQ74ZhOF9XZXq/ES437mtLTS0qE9SUKSH5anwLb+EMU4yTmcqQQFE7RN3o/QYhXR8L2WDsdxKrBv//VssiJ4Oh6KeXAhnYyFHmYcvsdREnt4cSY4N2xLc5QyP3OHFmdy2w4vzvTPmFQp55v11HlDy9UvBVZNrLLtanMAIIv2BqiW0hMnSc+Q9MOSnm1mN8R/XjCuk4eHUkV/Jx3jQZVixHCm5SjBMhsyvEv5jm20fbBvNkFLV4fSCcXhQEXBPzaZo5QdmnYyMKxuVOVhz+yQRNDb9+RKK5cIFs9HSr+efuLlhs9ROhGolp1caaleM+2bb6bikZSaz3RmvaNOz+lQEqf3fsyMHGZYPITtUKbadnKlpX2jqjuZpO+Ql9CEkuWhKykOmaO0MFPXXLOeq6od9BO0lXDsS6teAj/i57DitrXNAQAP7Q1QIaUvR+6c+weFVw8/J2ZJ5zDQMY7/HtUxHuw7SAAOLsxo/fRasPKRyCVOzqnd7enMekcXH1zQ189seMfnP+lfWmlpz2xDs37HOHCdZI7SocUZ3bm0mkpeDi/O6J7T68HqUPacBxcHCUkSe1Td8aoh3lA9//g9s43UAhOjFpYYJDTSqjdHaXjFKfrOiZXovtcsn6AdCiUKC834+MH2w4szuvf0emooZqNm6vTcyGQqm3w4N5ijdPnhRZ1Z74xMok+utjTXDFcp9841hsxR8iqKQ35m+xWwzM/h4cUZfW1pNb1tz4zuPrXmvc5oHtjXz2x452yrZlLPpSteVbVdbQ4AZNHeANUyDRWnbZG0UtlP5f0O+KiOsT+XSRpUjMws1bHuz2XKDMXyK0bZbcG5Nt51Di42ZUonD/39km3L+esk1zq8ZzZ10lHVoUMLw+OUpNVWR+vt3iC59K6fJF1+QpJbRCI+53l+BS2OfdFbNEEatpR6S4cWmzJZLhH0E6dBojDbP0cyFPPwnnxV7mC2IpiayxS+R06uP0fpYDZBGzJH6dDCTOq/zNQ9zlxnYaaeSiSHDdU7vCf5ORzco9z9yFU+nVelzA8jzVZDAQAAkFbdxCnbCV6NhmwtzjRyw+oOe9WQpGOcm/e02tbBxZk4oYnPudJOHZtsk9JD9U5mtiWGDatLOtvZjvHiTH1w7Go+6UvmKIWGCSb8pO3AQlP1muU723tmcscO5jgNrp/rgAcqU7n74dyg0rdnME8nmaOUux9xhcVffvtEPF9s31xzZHXozHoyFDPzHvlxJgnJ8qCSk73Hh71hhoP7niRog5+l7M9ckqCZ8u/FQW/BiP62hUHik8xRmqkni4x451yYSX18eSKuUs40at5+2fs+2HZeZr5aahipAAAAEFLdxEnJQgyRpWUvIclUQ/xOfb5jHB+/Eh6u1h/qlln0wJ/Tkx++J7U6PZ3d6Oi8PenrROeMO8ZeQrN3tqGm1zHOJjT+tlwVKu5Ye2HqRJygZSto0fGzuYUp/MpUMkfpvMVAB3xPOuFcWtlIzVFy8pOcwWs/uxHNUQrFfiiTfCTJgwLD9/z3aNT9OLwnnyTl5gPFrz093yyTfLjBHKXsIhInvITEfy8XZ+qaa3hJ8Eo6QfO3heYoJXPt+j9z/Xtkqf3ysQ8Zmhp47QAAAEirbOKUHVGcSkhi/Y5xKslJd5adiz79TzrG5p0gmX+TXTBicaau2aRj7FUpQsP3cgtThBKa1Sgh8atdoeF/oQRNiqs2i9mhZa18BS2e63JgYVDJyQ6Bc27wHKX8aoTtwOIMbR1caA6WfHfphMZPbCUvEczEmXo9gQpJaG5ZqCqXDMUMJcb5Z0BtaK5Z0/zMIOkMLf5xcjVfVUtiHwzvTO5RfN8t/f4M7kf6vofucVJx8ucoDc45qKBJSg1J7M9r8+ZIJVXKweIf1JwAAABCqps4xfy5IYe8hRCibXHHuFn3tuUTmsFCDE1J5nViB9UQ/3g/QUslNAVVKP+c/ralfoJmyg7vOuB1ggcLPgxiX2t1td7u9TvGgzjbuSFwSysbOrAwo7qlX2M6TjfYtidfDTkvW83wYk9iSiU0meTBT7yiOUptHV4cHC8NFupI3ffVlmbqNS0mlTUN5oENqimDOUq5RSTi66TveztOYv37lq/a5JPgwb79JMdbpCS7OmNoVbylTBUqPUcpvwBG//V4Me2da6hZS4Y+unxVbkgiCAAAgLzKJk6Wqzi1B8lHP5lq96s7g/3yCY2fkJg3hs5PCrLVEAtUWPxFBk5mO7GKPv1faXUDlTFvuJp3nf3zTe/ZUIPr+MPITnjPEkquncSeW3Qhrg7Ju86JQJyhYXFJQhWudvnD/AZzlPyV5XL3I56j1O05bwhcekigeWWXk8miGt59z1XlAnOURs3zOZmNPZv0acj76wZzlPyfhX6cC9GwOn/xj72zjWg1w/5++crnYN7SbC6Jzs+/a8UJZxJ7aL7a8ColAAAA0qqbOMV/hxKaoR3jIZ/A+x3jpHOazFHKJSShakjcMU6eJySFhmKlO+V+JWiQoKXP6c9rSVWCvMUQQosEJAs0ZBd36CcPmY5+vWbaPz9IqPzFDJJrJ89R8u9ncny2apO8nloqiR1x3xczSdJqnPDKTx6i4XepeT6BhTr62/akk5zTa6GKU7qimLz25DlK0TnDFbhUsj0IvT/vKVtxiu6bpRLOKPZBpXBwP5q5JOlQ5kOBkT/bgSrlYH4VmRMAAEBIdRMnb2hY0jEOVW1SHeP4E/jFmXp/mXFJqY5xktCcynaMA/NvonOGO7HZznY2IUklBavJMENLJ30LzVxFwcxfKju9gl1itdVVq9OLzplZKju5dmIpvk7NQpWtQac+XzFKz8/yDavuSOFE8qCXaHTjJcaz84ROxvdokKC51HOUsvf4PC9JOpWdo+TPLctUbbKr2jmXX4jBr+QkwwxHDsWM5ygNruLNUfIqgP4cpSSmtVZXa+1uLmH1k6nEydXBSor+fv5rBwAAQFh1Eyfv337HONuRzA6rGywhPajk9JMPb8jYUrY6lJxzpZW6TjK35GCmA34iU7WJjm2nz+nUn6OUrw61gyuuHZhvquZ1jEND4NIJWnpOzuF+cpgejpgYtuhCaDGDZI7SoQX/fg6ryrU106j1V/9LtiX3PRnadnqt3Z+jlF2Uw08okteZq3YFHno8WAkxMw9sOZ1IDip1metk51KFqnJy/aGYoWXpDy00c3Hum2uoUfeqlKkkx3LDEXP3wx/y6cLJFEP1AAAANqeyiVMi90m9vHkty5nOtpTvGGvQMT7gdaIHFYGkU69UxzhVpYifG5QdVrd/vqmG/+l/v1M/mOjvbwvNr0q9zpV2Zk6O3zEezOkJDYFLEqJQMpUdvre00lKjFs1RSs45WFJ7MEywP0cpdT9cethkPxHcSCU5SXVH8ip93usZrLSXnqPkL0OfJA+DexSY0+Od008uk6GYoTlKSeLSvx+r6aGY2SGSSUUySeAHC33493g2XA317luoAuffD4sf/pW8ztDxocqnmXRgnsUhAAAARqls4hScFD+sY+x1gkMJid8xTjr7qQpLfHg2QesfnzxcNZWgtTOVqfzzjfxtg9XZvMrHHr964AYVhX7w0Wuv10z75rzhaqvpzraT09mNjtpdl59flUkkk20HF2f6w/ekfPLhb0vm5PSv31/IIT1HKbUaoRtU5fxnSKXmgcV7Z5+j5F8nPRQzMEcpNTdskDz4QzG925maoxRtc4EqpX+Pm/3nTeXuR2bFx9D8u9T9WE6SHG8oZiD2tXZXG53BcvH+8amfzbhyenBhpj98j4oTAABAWHUTJ+UTEv+ZS8GOsXNex3iQ5Zz0KhdJNSbf2XapjrGfaSRzlPyhaVGC1kx1Yk/En/7vn2/2r39iJV0dci5+jlKnl0qSkspaKhnrd4zT81pOeglJcs50ghZesS25zonkYcLetnwlKDwnp+vi5yh5iaA0mKMkL/ZkjtL8TL2/b7Za5iSdWhs8jytbYfET2/5QzIV8MuS/lwoMzxx6vMvOUUpX4PyhfkuBn0N/jpK/0l52mGFSudw/31SjXuu/H9l7nEqmFtJJY3YeWL9KmZorR+YEAAAQMhWJk5k938xuM7MvmdkrxnPO6G9/PlK0cpilK0aB4WHpjrVLdYyTif7J/JsDSfLjBsOz/E/1kzlKBzPVpdCy5SdXojlK9dpgwJn/jJ9ka2pbZt5UdhW47Ip8qQrcnvxwtaFzlDIJTXbp7+Q5SntnB8P3Bg8THqx2d3o1nqMUJzT+3LCDgWpZkrgkQ9tCcWarcv49Ts1RUn6OkrzjD3jzptIJWrTzRicZiplPgrNzlJZWojlKzXqtvz01Ly6JPfNzmJsr58+vSt2PwXX8c/rb0j9zrv+A4vTPzEbuZ2Y32I42BwBCaG+A6ig9cTKzuqT/Lum7JT1W0g+a2WPHeY3BHKXBPKHQamIb7Z5WWt2ocuEdf3KllRqCliQPfsfYn4902FtFLnmOkl+5SK4fml81SNDSFZbQim25RRcyi1BIfsUo3QFv1Ex7Zxv9ffPPqpLOrncGc5SSyF3yjKBZ+fOJQs9RWopf+0HvwcMnsolgMqdnNf/coVylz6um+POEcnO2JHV6I56jlFkkZGm1pT2zDc026oNKTmaOUnIvB/conQSH7tGhTOyjqpR+pc6vfCaSitHgflj/Pa/1q5TZpG/wfqxsdNXq9vrP7hrcj3bq53A3mESbAwAS7Q1QNaUnTpKulPQl59xXnHMtSX8p6apzPWm2GpLtGIcqCtmhVP3jM1Uo5TrG0RCrpeWN3PHpyfvpmPxluvsPkPUqCklnuWbSvrlmbn5VOslpq9NzwRXsDmXmVw2W87Z+7CeW86uzZRerGBzfziRJUXXpYKoylX6O0uB+bPS3Rfs5dbo9ncoOd0vuUbLNBslUMkcpmSeUXsEueS/b/XuUet9WA0lwnPRF991ySV8/dv8eeccPG1YXWqgjOxTz5Gq+YpTMUconwYN5XNF7FFXQDsRzlJLYQ3PtkgQ+u0Ji6h5rUPEqk5l90czeZWavNrOrzOzyMV9iW9ocADuXmT3bzP7EzH7bzH7MzJ5qZrPFRxaivQEqZBoSp4slfc37+q54W4qZXW1mx8zs2PHjx7d0Ab9jLGWGd3nDrlIr2CX7ajB5P4oj3nclXR2Kjm8POsZKz1HyqynLGx1vjpJX+QglY3HSVouH7/nzq/wE4ESmcuHHHqyGBBK0wfGWqnb5w7u6PadTq61MBW0wvyq7bbZR03y8EIOkfoLmJ5fJHKXsQh3Zqk2S+BzMxJ5dCj06dlDpS72X2QUS4jlKqblM3jDDdPIx/DlKh/aMHlaXJDT+HKX8nC1Tdo5SP6Y4mTu8mD7nUjxXLrlH2SGSynwo4M/56rn8KoHJHKuS/ZGk+ySdUPQp7c1mdpOZ/Scza44+dFMK25xzaW8A7EhvlfR+SZ+UdIWk/yDpljGcd9v7OAAmpzBxMrP/ZWa/GH8ac3gbYggNEsr13pxz1zjnjjrnjh45cmQTJ/U6sV7HeNDhHHSME+mKUXT8Wjv9gFF/ONRhr1OfJGP9jnGw4mT5bZm5RyOHZ1k6ofETotDwrJ6L5ygt+n1NFy9bnk4kl1bamqnXtDhTzw9n9BK8M2tt9Vw+QfMf7GrxDRnM47LwYhfKJLGpqs0gaRxcJ6oOHd6TTihSSXB8/InldEIiZeYoZYYEHsq+l6st7Y2HYiZCwzv9OUq5RSi8JDh57fmfw3SSFMUzqJalFykJLzefXT3v5Eq0kmKyXLx/P/yfmX6VMlPZmgI/5Jz7aefc7zvn/h9J3yrpbyWdkfSGMZy/sM3ZansDYMf7knPu3c65tzvn/r1z7irn3CPGcN5t6eMAKMdmKk5vif9+qaSPmtmXzez9ZvYbZvYvxhDDXZIu9b6+RNI953rS7KILqUpOPEdpb2byfnYFO8l7uOmCd7z3zKNEUvkY7Bcfn1nBLrutf7zLJklxTMv56o6/xHh+kYDByILTa9FzlKKKUXpOz+F4v2Se0GCOkuUTGj8hCdwjOYXn5GQrRkqvRhhciCHer91x/TlKSZzJ6zyYu8ft3FDMUJzp5CE/Rym579mEJJv0+QlvqNKXbE9Vh7xlyxVHn9yjmkn74jlKSVUtez/ObrT7c5Sio62/4EP/59AGi4QcXGj2q5Sp++F9KJB9EHJyP6bAaTN7QvKFc+4GSU9zzr1e0jPGcP5taXMA7Dzxh8M/L+mfzOzfbcMlaG+ACilMnJxzH3XOvcE591Ln3JMkPVrSr0m6TdHY3XP1aUmPNLOHmdmMpJdIet+5njQ71C606MLhQEVASlcUTmTm5AzrWGerDP3jlwcd40RokYGz6+3oOUpeUuAvMT64tvegWrP+Kz2xPCJBG/IcpT7nUkPg+vctMATuxPLgfiSxt3vp5yglsWeTqeQeD+YoDZLQwTnj2DPLxfeHtq1m7ru3yl+yLX2PB689u4Kd/36kE+P8fY/2i177AX8o5nI+yVltDeYoJcdn58r1lw735ygpWezCv8eDGCWlkqTkPiUVOH/VRf86qde+J1SlTC+WMQX+raQ/jecbvMzMfl9SL/7ezIjjNmtb2hwAO9JbFPWFLpT0w2b2VTN7n5n95zF9OEx7A1RIo3iXNOdcR9KN8Z9z5pzrmNnPSvqgpLqkNzvnznlccarCspKvBIWqO0vL3uR9bz9JqcrJykamY6zB0LJLDy2kz7my4XWMI35n+55Ta/0YJa9y0V8xrq2nPjS0MEV2ztVgMYPljU4Ue2AIXHaOkl8J6nfAbbAtmaPU6vbScXpVilOrgzlKyb1PEprLcvcjP0cpPXwvSXzSqxH25wmtDEnQAq8nifNrS6vRfQ9UtrJDMZNvnFhu6aL9c6lzLq20dGAhGoqpQELzVUuuk11J0frDLp94yQH5kucoJZKFNtLH55N6vyqYSpKSeWDZDwWWByspZs+ZfQZV2ZxznzezKyW9WNITJH1J0qvNbFHRxOpzPf+2tDkAdh7n3EclfTT52swaila/e6KiD4fffo7np70BKmTLidN2cM79taS/3o5z5zrGijvroY6xN3k/6en7w7sS2U69vIpC0jH25zMdynbKg5Wt9IIPUcc4eYjrIElK5sAczAzfO7HSUrNu2hPoGB9anNGZtSiZOpWZo+RX0B77kH392JNFCg5nqhHB4YwrQ6pDgaF6uSF9meco9a+znE0UTK1uT8sbnfT9cMrNe0quk5ujtJy/79lV/qRBVS+5H/7PwrChmP5Ke7n7kcTpPw/M23cwbDLa5s9R6u+XvR9mOrPe6c9RSs6ZDAl8+JE9qXOeWNkYrKSYGYZ6eHFG3V6UME3JUD0557qKOizZTstvjOn829bmANi5xv3hcHxO2hugIqZhVb1tMarDOWw+Uagy5c8NSY7PLhIQqmL5SZI/DCy6zsbQOUp+TGfiOUoHU8PIXGZ41uD4ZPjeYE7P8OXRs5WLpdyQwPxzg6Tsw2bTw+L8OT2dXjRH6aB335KY8snDYI5StpLjx5SvkFhg2GT0ndNrbW/VwfwcpaHXkbyqXn5uWXYo5lJ2jlLqnINEcDk7R8kfvrfor4qXrBwYz1EKrYoX7+svO+6f86S31L3/85H9OQzN2ZqWxAkAAGDaVDdxyk7ezy03HZh/43fqvU/l/TlKpvAQp7PrnWiOktcJljId4/58k8Ezj7LJx8jhWV51KJuQhJKHfkLjVWOyc2XMTJ2eK5yjlJxzpdXVfLOu+RlvifHAPc4mD3415VBSWfKeZXQwcN+ie+zPDYsSwdRzmJKqXCYpiI5Nx57EOWyOUrLvSqsTDcUMVJdyieCK/xwlZc452z9+vd1L3aNkkZJcwppdzCT7c+y99qRK5P8c91dSzFYkl/M/h0srG4OVFJPFMgQAAICQCidO6bky/vyd1WSOkjcHRYqGsR0MfCqfdIyTXc+uR8Pe/CrWUuYBo8n1V1vd/HC15Y1cQjEYEjhIIEJVhv4cpRFD4PzYkzlKoQUBkuNPrbbkXLYq59LzZ7x7m+vUL2cW0JBflUs/PzBaDnyQUPRjDyQk/v2UBtVDv8Ky0elptdXNJUmSAsPqNnLLxWev48d+OHPO1LDJJPbljdQcpeic+flZ/Zi8c/Zyc5SiZOpEcOXAjdQcJf+cflWw56Kfkez9OOmvpNg/p7eSYr/iROoEAAAQUtnEKZGvsAxWbMt2jJ1LVzOi4zdSn+inEgj/nIHqTiI7H+mMP4TNu44/R8nMtNrqRtfxrn86maOU6ZS3Mqu4SYPKhZkNHZomhedcDVuS298vm4z15yilhtWlE43o+HzlI3ROf46SmWklcz9CyaXvYCb2kyvt1Byy5B6lzmn5eW1JIpqdGyZlVt/T4Dr+c5T8Kph/Hf85SklM2YpTcqH+g4zjr/37mf05jq4zWLZcUpxcDt4fSakEDQAAAKNVOnHy58WEqiHZuTapbYF5T9H26BvpjnGocjEwMvnwruM/eDcV04j5Vb7s8KyVVjdXQXtgOVsZi5Ku7Dk7XRdM8CQF53EtxkuM9+9Hpgrl9+qzCcnScuD5VSutVBI77D06vZZd0c+7H5l7tLzRCT5vKVlJMTFsmKGk3BL27a7LJbEnVjb6c5SGxZSuyqUT1tRcOS/O1P0I3s+B7LLlUv7nw39Oln99AAAA5FU6cZKizmG2Yzzo1Kef/ZPaluy7knnmUezgwozXMQ4sGDFkeFZ2Wz+mzJyrcGc7vGhCP6ZABzqUoKXnKHkdcC8BWGvH1Z1ABS2Zo5QkD61OL7WfmVfJCSWSmQ78A5nkNIpzI5U8+K/zwHz6+v7rTMU5IqFI3vNkblcyFNO/x8GkMTD8LltlXMpUcoJJp6KkK7XNkjlb7eDqfaFzJnOUskFlh+X51xlWAZNE5gQAADBEpROnpC94YL456BhbNJdJyg+h87clHcnsp/LJOf0V00zqD6sLzbUZldAkW7OLVSSd25lGTQv9jnG4gpYIDdkKLX9dnKB5HetAopCthvj7JZJhdQdGHD9qmOHScrbSp/h88XLx2XNmEqJhsWcT4+hY7700yz27K/U6Q0MCRyzU4Us9RylQCTKZ2l2nrj98z78f/jnj7ckcpdxrCn4oMCqpjzZOw3OcAAAAplG1E6e4M5hNfJL576GKRCh5OBjowIeSrtTk/aLkI7PNX+DA5z9HyTQYVpcdLpfaFkimErnKlve97Kpr6W2bS9Cic0bf2DvX0Eyjltombe6+rwy5H9kEbdQ5QxWW4CISmfvRc4qXiw/dj1GxR/LJqfWPDSU5IxOaeM+zG51U5TPZPmyO0pYrcPHXrA0BAAAQVu3EKf471Ikt6hhraMc6kIz1P/1PPyw2dHx229AhbCMSNH/f0PGjKjF+dcfft/8cpew5Ry26MCT25JyhY/3tWxnO2E8UArH7QzGD1aUh1Z3+6wned2+OUiAJVvAeR1/7i0j4QnPl/NcUTPoCsfsO78nfj2Qlxdw5Q1W5wHUAAACQV+3EKZB8JP3DYZP3sxUaaXiS5F0pd2y4GhOqfOS3+XuGkr65Zq0/Ryl9fPEKdtJgjpIvNI8r2jeQEAU64KHYQ5U6yZujFOjAp2IqqPT5VZd6IMkJJZ2HChPjwLH+vnvyiVfRQh2hRDI52p+jNGop9dzxI+7HIS+BHzWvTsoneBScAAAAwiqdOCUKO+WhhMbfd5Md41BHfaax+Y5xKPEKJSSHM89G6p9zk9Wh0LC60OsZdv3Rc7aGD5GU1H+OUi7OQDIWWlXPn1vmV4dCQvdpM8lDdlty7c08Ryl8/PAqZSrJSR0fSm6LEuvh26QhC1tkYmeoHgAAQFilE6fQ8C4FEpJkW7pjHBie5QkvGBFIHoYsMR5arS6c0PgLFyT7pRem6B+/2QUSQklfIHb/OUq+4mF16Xj8fYuG7/nS71GSfAySoVHJQ2q5+FCSM6TK6L9voXhCQzFDc77CSfDoZMhXeI9HVFNDQxz9lRTDPzPR1ywOAQAAEFZq4mRm/9XMPm9mN5rZu83swHgvEP0VrFwUdIy3Uk0ZHJ9PaEIJiX/9ooQmnSgMhqZl98s+R8l/TcpsDc0TOlTQAffjL6rg9V9PcDGDfCKYmqM09PjANgvcj348o+cohZIHP6bQ+xa6R/6+W6pSBl9P9Pdcs6aFmfwiI8OG5eViLxomGEjGkk1Vrzhte5sDADHaG6B6yq44fVjS451zT5D0BUmvHOfJw5WPfBUqlExtdriav2+4MpWvMqTnKPnXye8b6lhvdtEEKTxn63DB8cMWpkhOW/gcpUDSF9ovOEcpMLTMN3Le0pDYwxWncOyD4ZCh9zyftBXNUUofH/g5DCVTw+ZXFQ7LG57YprcNTnogSWT7FafK29Y2BwA8tDdAxZSaODnnPuSc68RfflLSJeM8/6hJ+cHhWYGOqT9HKXXOB/Ppv4YfO/T4TS6aEEpIpNHPUfKDKlo5MDlHao5SqlM/ujoUvE7/HuUrU8P2LRwiOSKRbNRM+5Lhe97rCiY0geP9OVP+8LtBEuZX9YoS3vw2jUjqc/uO2BZMOAOvJ12lzA9BrKLtbnMAIEF7A1RP2RUn349L+pvtOHGo01g00T7hz1FKHx+ae5TvwBct+OD3V4NDsYoqLIFjE3tn/ecojT5+M+c0WeHCA9H2EUPoNrnAgb9c/NDjRyRjoUrb0OXii4YZht7z5JybeI5SsBK02UQwPnbBS3JSMRXGnvwcBaqZgXl2//49N+v0ajv0kqpo29ocAMigvQEqoLHdFzCzj0i6MPCtVznn3hvv8ypJHUlvG3GeqyVdLUmXXXbZ5q49Mvko6hgHhrDJT4g2t0hB4epq8d+LwzrGm1ytrvDa/tLhwTlBgeQj27G2bHVosP+Bwnk1gXlCI6py/hyloccHYu+/xi2sHOgnEMH7MeL4oucoFZ5zk1XK/GscfnxoAY2tDE2d9yqsO9E42pwH094A2H3K7OMAmLxtT5ycc88d9X0ze6mkF0l6jnPDp6Y7566RdI0kHT16dFNTMUYmRKFkKtj5Hz1HyVe4wlkgoegfO2R1tfCS3KMrF6Fz+nWWAwuBikTBqniJ0LylA/PN/hylousHF03YVLUr3jeVqIxIJIMLPuR/Dvw5Sv6+m162PPCeH85Wtkadc5MVxdCwSSnzXvaT/cAwwU3O85PUr1LuVONocx5MewNg9ymzjwNg8rY9cRrFzJ4v6VclPcs5tzr28yvqGO+ZHbzMzS6wkMh14C1fwRo1R6pogYNhE/qT7aFKzmYrF6Fhgvvnm6klxkPzkYYljZbbb1hVbvj1Q/O7ioZSJtv95eL96/j7zjRqqpl0wb7NVgSb6SRnZCVnpmC/4T9H0ubnpm1mrptpeJUyeM5AVS2UTFXddrc5AJCgvQGqp9TESdLvS5qV9OG44/ZJ59z/M66Tm5kOLjbC81o2uUJZdq6MqfjZO1vZNqq6489RKoq9aI7RsHg2u+y5JF2wb04PP7InF3sokayZUnOU+tcPDG07HEpIHsRzlCRp71xTf/FvnqYnXHIgd51NVbY2ue/IxT8yxz78vD36hov2paqUo4ZIhldsTN/LFz/lEj36wr3p2JPrhxK0outo19jWNgcAPLQ3QMWUmjg55x6xnec3hTrG+U/1k2FmoXlL2SRpcbbRf8bO4Jz5yfsX7p9Xs266/LzF1H7Zcw7rqIcStFGLFBQ+86ifDOWrSMlqeUUxffDnn6lm3XL7haohBxdmUnOUzt87p5l6TZccmFdWUfKQxJ9Lpobs+81XHM5dQxo2FDMfe3aOUm2rc8sycf7AN12qH/imS9PXSZLGokpQvwqVHiL51Ice1FMfejC17UmXHtCzH3O+Du/JD6fc7BDJqtvuNgcAErQ3QPWUXXHadqHOYbZjfPShB/Wf/9nj9bQrDhUe//p/8UQ1MvN5zPKd5SddekA3veZ5qWQqGTL4kP1zg2OHLP194b65/HUClaDLz1vUoy7Yo8c9ZH8u9s0kaNLwOUrZfbPzukIrBybXOpC5H894xGF9+tefm0nQhick2XMe2TubW6zCZLk5SmHDh9ANqwj6la1nPvKIfvE7H6VvuGhfLvYHm5Ak9+5A4IHA4aXM89W7rG++4nAuabz4wLxmGzVdfthL4EcMCQQAAEBYtRMnC1dDsh3jRr2mH37aQ9P7DakePCRQMYmW6Q4tIpHu0D/08KLe/dPfoiduYhjZb37f49XLTA81k/bNNVJzlC7YN6cP/cKzUvvtnWtovlnXI7xhdYlQglezbIIWjilr+LOm8tUhM0slTf51ip6jJEUJa9b8TE3n75stnJ+zf76pRi1T/RuSsIaS4P0LTf3ccx6Z2nbFkUV9x6OP6JsuP5g6NnTOELMokZ5t+HOU8pXPxLDhoUWeeOkB3fqfnp+q/u2Jn2Xl/ywn9yObrAMAACBS6cRpz2xDF2cSnVDHOGRYx3qYzXZsn3zZweD2bEyNen5uU5L0Fdk719R1//656SWxk059pmrzo894WO748/fN6bw9s8E5Sr5mvaanXHZARx+af01bGQZWtFy8FF4i+2XPfqT+9Tc/NLc968jeWf3TK5+j87zXXjOpWTedvzedpGSfVTXMvrmm/vTHrswdK23uZ+EZjzgvl1gnDgaS8NA92qxaJhl62HmLev/LvlWPe8igguYUZel+VQ0AAAADlU6c3vaT35yb79GomfbNFw97WpyNOrUXesPqhrlw/5weeX6+urMZoflIw3e2TSdouXlYIxZdyPrBKy/T9z354lyHO6teM73rp5+R237+3lldcWQxcETaoy/cq2+4aJ8u2B+YW7aJOC/YN6cL9hW/P1KUPPka9Zre9pNPyy2wUAtUKTctsBDDMN/zxIfoe574kNS2hx5e0CUH58MLfYx5PtLjL04P7TywMKPffcmT9K2POG+s1wEAAKiKSidOVwSGqv3aC75hUw/4fOYjj+i9P/OM1Cpyw7z1J75ZD3aE0/n75lSvWTDWrCsvP6j5mXN7yzbTAa/XTIuzD/467/ipp6tRK34W0DMecZ7+5uXfltp2ZO+sahYNa9xuVz4sP6ftV7/7MQ/62rV+krOJJDjgRU94iF70hHQytX++KTPlKqfb4aonXbzt1wAAANipKp04hQxbcS2rVjM98dIDm9r3XB4Y+ojz9+jm1zxvU8ncq1742Ad9nYsPzOtRF+zRkzb5ms5Fttq1FU+89ICu+/XvHHuFZbPOJXk4+tCD+pnveHhutbtz8W2POE8f+cVn6dJDC2M7JwAAALZu1yVO02gzSdO52r/QzC0iMa3KSprO1eJsQ7/8vMeM9Zy1mm2q6gkAAIDt9eBLJQAAAACwS5A4AQAAAEABEicAAAAAKEDiBAAAAAAFSJwAAAAAoMBUJE5m9ktm5syMp28C2Ha0OQAmhfYGqI7SEyczu1TSd0q6s+xYAFQfbQ6ASaG9Aaql9MRJ0u9I+hVJruxAAOwKtDkAJoX2BqiQUhMnM/teSXc75z5bZhwAdgfaHACTQnsDVE9juy9gZh+RdGHgW6+S9GuSvmuT57la0tWSdNlll40tPgDVMo42h/YGwGbQxwF2F3OunOqxmX2jpI9KWo03XSLpHklXOufuG3Xs0aNH3bFjx7Y5QgCjmNl1zrmjZcexWQ+2zaG9Acq3W9obiTYHKNuo9mbbK07DOOduknR+8rWZ3SHpqHPugbJiAlBdtDkAJoX2BqimaVgcAgAAAACmWmkVpyzn3OVlxwBg96DNATAptDdANVBxAgAAAIACJE4AAAAAUIDECQAAAAAKkDgBAAAAQAESJwAAAAAoQOIEAAAAAAVInAAAAACgAIkTAAAAABQgcQIAAACAAiROAAAAAFCAxAkAAAAACpA4AQAAAECB0hMnM3uZmd1mZreY2evKjgdAtdHmAJgU2hugWhplXtzMvkPSVZKe4JzbMLPzy4wHQLXR5gCYFNoboHrKrjj9lKTfcs5tSJJz7v6S4wFQbbQ5ACaF9gaomLITp0dJ+jYzu9bMPm5m31RyPACqjTYHwKTQ3gAVs+1D9czsI5IuDHzrVfH1D0p6mqRvkvR/zOwK55wLnOdqSVdL0mWXXbZ9AQPY0cbR5tDeANgM+jjA7rLtiZNz7rnDvmdmPyXpXXEj8ikz60k6T9LxwHmukXSNJB09ejTX6ACANJ42h/YGwGbQxwF2l7KH6r1H0rMlycweJWlG0gNlBgSg0t4j2hwAk/Ee0d4AlVLqqnqS3izpzWZ2s6SWpJeGStgAMCa0OQAmhfYGqJhSEyfnXEvSD5UZA4DdgzYHwKTQ3gDVU/ZQPQAAAACYeiROAAAAAFCAxAkAAAAACpA4AQAAAEABEicAAAAAKEDiBAAAAAAFSJwAAAAAoACJEwAAAAAUIHECAAAAgAIkTgAAAABQgMQJAAAAAAqQOAEAAABAARInAAAAAChQauJkZk8ys0+a2Q1mdszMriwzHgDVRpsDYFJob4DqKbvi9DpJ/9E59yRJ/yH+GgC2C20OgEmhvQEqpuzEyUnaF/97v6R7SowFQPXR5gCYFNoboGLMOVfexc2+QdIHJZmiJO5bnHNfHbLv1ZKujr98tKTbNnGJ8yQ9MIZQtxMxjgcxjsdWYnyoc+7IdgYzbpttcx5keyNV7z0uCzGOR5VirGx7E+9LH6c8xDgeVYpxaHuz7YmTmX1E0oWBb71K0nMkfdw5904z+wFJVzvnnjvGax9zzh0d1/m2AzGOBzGOx06IsQhtzmjEOB7EOB47IcZRaG9GI8bxIMbxGEeMjXEFM8yoRsLM/pekl8dfvl3Sm7Y7HgDVRpsDYFJob4Ddpew5TvdIelb872dL+mKJsQCoPtocAJNCewNUzLZXnAr8G0m/a2YNSesajO8dl2vGfL7tQIzjQYzjsRNiPBe0OcQ4LsQ4HjshxgeL9oYYx4UYx+OcYyx1cQgAAAAA2AnKHqoHAAAAAFOPxAkAAAAAClQ2cTKz55vZbWb2JTN7RdnxJMzsDjO7ycxuMLNj8bZDZvZhM/ti/PfBCcf0ZjO738xu9rYNjcnMXhnf19vM7HklxvgaM7s7vpc3mNkLyorRzC41s78zs1vN7BYze3m8fWru44gYp+Y+7lS0N1uKifZmPDHS5uxStDdbjos259zjo71JOOcq90dSXdKXJV0haUbSZyU9tuy44tjukHReZtvrJL0i/vcrJL12wjE9U9JTJN1cFJOkx8b3c1bSw+L7XC8pxtdI+qXAvhOPUdJFkp4S/3uvpC/EcUzNfRwR49Tcx534h/ZmyzHR3ownRtqcXfiH9uZBxUWbc+7x0d7Ef6pacbpS0pecc19xzrUk/aWkq0qOaZSrJL0l/vdbJP2zSV7cOfcJSUubjOkqSX/pnNtwzt0u6UuK7ncZMQ4z8Ridc/c6566P/31W0q2SLtYU3ccRMQ5Tynu9A9HebAHtzXjQ5uxatDdbRJtz7mhvBqqaOF0s6Wve13dp9M2bJCfpQ2Z2nZklS5Ne4Jy7V4reeEnnlxbdwLCYpu3e/qyZ3RiXuZMScakxmtnlkp4s6VpN6X3MxChN4X3cQab5PtHejNdU/p7Q5uwq03yPdkp7I03p70nA1P2e7Pb2pqqJkwW2Tcu6689wzj1F0ndL+hkze2bZAW3RNN3bP5T0cElPknSvpN+Ot5cWo5ntkfROST/vnDszatfAtrJinLr7uMNM832ivRmfqfw9oc3Zdab5Hu309kaarvs7db8ntDfVTZzuknSp9/Ulip7gXTrn3D3x3/dLereisuDXzewiSYr/vr+8CPuGxTQ199Y593XnXNc515P0xxqUWEuJ0cyain5Z3+ace1e8earuYyjGabuPO9DU3ifam/GZxt8T2pxdaWrv0Q5qb6Qp+z0JmbbfE9qbSFUTp09LeqSZPczMZiS9RNL7So5JZrZoZnuTf0v6Lkk3K4rtpfFuL5X03nIiTBkW0/skvcTMZs3sYZIeKelTJcSX/JImvk/RvZRKiNHMTNKfSLrVOfcG71tTcx+HxThN93GHor05d1PzezLMtP2e0ObsWrQ34zE1vyfDTNPvCe2Np2j1iJ36R9ILFK2o8WVJryo7njimKxSt4PFZSbckcUk6LOmjkr4Y/31ownH9haLyZVtRBv4To2KS9Kr4vt4m6btLjPHPJN0k6cb4F+CismKU9K2KSrw3Sroh/vOCabqPI2Kcmvu4U//Q3mwpLtqb8cRIm7NL/9DebDk22pxzj4/2Jv5j8YEAAAAAgCGqOlQPAAAAAMaGxAkAAAAACpA4AQAAAEABEicAAAAAKEDiBAAAAAAFSJxwzszssJndEP+5z8zujv+9bGZ/UHZ8AKqFNgfApNDewMdy5BgrM3uNpGXn3OvLjgVA9dHmAJgU2htQccK2MbNvN7P3x/9+jZm9xcw+ZGZ3mNmLzex1ZnaTmX3AzJrxfk81s4+b2XVm9sHME58BYCjaHACTQnuzO5E4YZIeLumFkq6S9FZJf+ec+0ZJa5JeGDcsvyfp+51zT5X0Zkm/WVawAHY82hwAk0J7sws0yg4Au8rfOOfaZnaTpLqkD8Tbb5J0uaRHS3q8pA+bmeJ97i0hTgDVQJsDYFJob3YBEidM0oYkOed6ZtZ2gwl2PUU/iybpFufc08sKEECl0OYAmBTam12AoXqYJrdJOmJmT5ckM2ua2eNKjglAddHmAJgU2psKIHHC1HDOtSR9v6TXmtlnJd0g6VtKDQpAZdHmAJgU2ptqYDlyAAAAAChAxQkAAAAACpA4AQAAAEABEicAAAAAKEDiBAAAAAAFSJwAAAAAoACJEwAAAAAUIHECAAAAgAL/P9XxzagX9DyIAAAAAElFTkSuQmCC\n", 170 "text/plain": [ 171 "<Figure size 864x432 with 6 Axes>" 172 ] 173 }, 174 "metadata": { 175 "needs_background": "light" 176 }, 177 "output_type": "display_data" 178 } 179 ], 180 "source": [ 181 "%matplotlib inline\n", 182 "from matplotlib.pylab import plt\n", 183 "\n", 184 "def plot_t_evol(out, labels):\n", 185 " fig = plt.figure(figsize = (12, 6))\n", 186 " \n", 187 " max_abs = 8\n", 188 " \n", 189 " ncoord = out.shape[1]\n", 190 " ncols = 3\n", 191 " nrows = ncoord // ncols + (ncoord % ncols)\n", 192 "\n", 193 " for i in range(0, ncoord):\n", 194 " ax = fig.add_subplot(nrows, ncols, i + 1)\n", 195 " ax.plot(t_grid, out[:, i])\n", 196 " ax.set_xlabel(\"Time\")\n", 197 " ax.set_ylabel(labels[i])\n", 198 " ax.set_ylim(-max_abs, max_abs)\n", 199 " \n", 200 " plt.tight_layout()\n", 201 "\n", 202 "plot_t_evol(out_del, [\"$L$\", \"$G$\", \"$H$\", \"$l$\", \"$g$\", \"$h$\"])" 203 ] 204 }, 205 { 206 "cell_type": "markdown", 207 "id": "6f5585c7-2429-4f26-9cc8-8e8269e254d5", 208 "metadata": {}, 209 "source": [ 210 "It can be confirmed by direct numerical comparison that the integration of the Stark problem using ``kepE()`` produced results matching those from the [previous tutorial](<./Comparing coordinate systems.ipynb>) to machine precision." 211 ] 212 } 213 ], 214 "metadata": { 215 "kernelspec": { 216 "display_name": "Python 3", 217 "language": "python", 218 "name": "python3" 219 }, 220 "language_info": { 221 "codemirror_mode": { 222 "name": "ipython", 223 "version": 3 224 }, 225 "file_extension": ".py", 226 "mimetype": "text/x-python", 227 "name": "python", 228 "nbconvert_exporter": "python", 229 "pygments_lexer": "ipython3", 230 "version": "3.8.6" 231 } 232 }, 233 "nbformat": 4, 234 "nbformat_minor": 5 235} 236