{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# One Variable Equations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Throughout this section and the next ones we shall cover the topic of solutions to one variable equations. Many different problems in physics and astronomy require the use of complex expressions, even with implicit dependence of variables. When it is necessary to solve for one of those variable, an analytical approach is not usually the best solution, because of its complexity or even because it does not exist at all. Different approaches for dealing with this comprehend series expansions and numerical solutions. Among the most widely used numerical approaches are the Bisection or Binary-search method, fixed-point iteration, Newton's methods.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- - -\n",
"- [Bisection Method](#Bisection-Method) \n",
" - [Steps](#Steps-BM)\n",
" - [Stop condition](#Stop-condition-BM)\n",
" - [Error analysis](#Error-analysis-BM)\n",
" - [Example 1](#Example-1)\n",
" - [Example 2](#Example-2)\n",
"- [Fixed-point Iteration](#Fixed-point-Iteration)\n",
" - [Steps](#Steps-FP)\n",
" - [Example 3](#Example-3)\n",
" - [Stop condition](#Stop-condition-FP)\n",
" - [Example 4](#Example-4)\n",
" - [Activity](#ACTIVITY-FP)\n",
"- [Newton-Raphson Method](#Newton-Raphson-Method)\n",
" - [Derivation](#Derivation-NM)\n",
" - [Steps](#Steps-NM)\n",
" - [Example 5](#Example-5)\n",
" - [Stop condition](#Stop-condition-NM)\n",
" - [Convergence](#Convergence-NM)\n",
"- [Secant Method](#Secant-Method)\n",
" - [Derivation](#Derivation-SM)\n",
" - [Steps](#Steps-SM)\n",
"\n",
"- - -"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"import numpy as np\n",
"%pylab inline\n",
"import matplotlib.pyplot as plt\n",
"# JSAnimation import available at https://github.com/jakevdp/JSAnimation\n",
"from JSAnimation import IPython_display\n",
"from matplotlib import animation\n",
"from scipy import integrate\n",
"from scipy import optimize"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- - - "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bisection Method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Bisection method exploits the [intermediate value theorem](http://en.wikipedia.org/wiki/Intermediate_value_theorem), where a continuous and differentiable function $f$ must have a zero between an interval $[a,b]$ such that $f(a)f(b)<0$, or equivalently, there must be a value $p\\in[a,b]$ such that $f(p)=0$. Below the algorithmm is stated explicitly."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Steps BM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
" \n",
"
\n",
"\n",
"\n",
"1. There must be selected two values $a$ and $b$ such that $f(a)f(b)<0$ and $p\\in[a,b]$ where $f(p)=0$. In other words, though we do not know the value of the root, we must know that there is at least one within the selected interval.\n",
"\n",
"2. To begin, it must be set $a_1=a$ and $b_1=b$.\n",
"\n",
"3. Calculate the mid-point $p_1$ as\n",
"\n",
" $$p_1 = a_1 + \\frac{b_1-a_1}{2} = \\frac{a_1+b_1}{2}$$\n",
"\n",
"4. Evaluate the function in $p_1$, if the [stop condition](#Stop-Condition) is true, go to step 6.\n",
"\n",
"5. If the [stop condition](#Stop-Condition) is not satisfied, then:\n",
"\n",
" 1. If $f(p_1)f(a_1) > 0$, $p\\in(p_1,b_1)$. Then set $a_2=p_1$ and $b_2=b_1$\n",
"\n",
" 2. If $f(p_1)f(a_1) < 0$, $p\\in(a_1,p_1)$. Then set $a_2=a_1$ and $b_2=p_1$\n",
"\n",
" 3. Go to step 3 using $p_2$, $a_2$ and $b_2$ instead of $p_1$, $a_1$ and $b_1$. For next iterations the index increases until the [stop condition](#Stop-Condition) is reached.\n",
"\n",
"6. The End!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Stop condition BM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are several different stop conditions for this algorithm. The most used are stated below:\n",
"\n",
"* A fixed distance between the last two steps (absolute convergence):\n",
"\n",
" $$|p_i - p_{i-1}|<\\epsilon$$\n",
"\n",
"* A fixed relative distance between the last two steps (relative convergence):\n",
"\n",
" $$\\frac{|p_i - p_{i-1}|}{|p_i|}<\\epsilon\\ \\ \\ \\ \\ p_i \\neq 0$$\n",
"\n",
"* Function tolerance:\n",
"\n",
" $$f(p_i)< \\epsilon$$\n",
"\n",
"All these conditions should lead to a desired convergence expressed by the $\\epsilon$ value. However, the first and the third conditions present some problems when the function has a derivative very large or close to $0$ as evaluated in the root value. When the function is very inclined, the first condition fails as a convergence in the $x$ axis does not guarantee a convergence in the $y$ axis, so the found root $p$ may be far from the real value. When the function is very flat ($dF/dx\\rightarrow 0$), the third condition fails due to an analogous reason.\n",
"\n",
"A final stop condition that does not have mathematical motivation yet computational is a maximum number of allowed iterations. This condition should be used not only for this algorithm but for all iteration-based numerical methods. This condition guarantees a finite computing time and prevents undesired infinite bucles.\n",
"\n",
"* If $N>N_{max}$, stop!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Error analysis BM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we suppose $f\\in C[a,b]$ and $f(a)f(b)<0$, the Bisection method generates a sequence of numbers $\\left\\{p_i\\right\\}_{i=1}^\\infty$ approximating a root $p$ of $f$ as:\n",
"\n",
"$$|p_i-p|\\leq \\frac{b-a}{2^n},\\ \\ \\ \\ \\ i\\geq 1$$\n",
"\n",
"From this, we can conclude the convergence rate of the method is\n",
"\n",
"$$p_i = p + \\mathcal{O}\\left( \\frac{1}{2^i} \\right)$$\n",
"\n",
"This expression allows us to estimate the maximum number of required iterations for achieving a desired precision. The next figure sketches the number of iterations required for some precision."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAGBCAYAAACq6fyLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VOXywPHvJKErIooFRGMBBS66NrCTUENHeg8qWPFa\nuD+xgWDFcoWLXRESeu+9JYoVUVAURFBRimADASVAyPz+2E3OGikpu3u2zOd58ty8h91zJnPXTN53\nThFVxRhjjAmkOLcDMMYYE32suBhjjAk4Ky7GGGMCzoqLMcaYgLPiYowxJuCsuBhjjAk4Ky7GGGMC\nzoqLMcaYgEtwO4DiEJFywKvAASBTVce7HJIxxhgif+bSFpisqrcCrdwOxhhjjFfYFRcRGSkiO0Vk\nbb7tKSLytYhsFJH+vs1VgC2+7w+HNFBjjDFHFXbFBRgFpPhvEJF44GXf9ppAFxGpAWwFqvpeFo4/\nizHGxKSw+4WsqiuAXfk21wE2qepmVT0ETARaA9OBdiLyKjA7tJEaY4w5mkhp6Psvf4F3xlJXVf8C\nbj7WG0XEbvtsjDFFoKpS1PeG3czlKIpVIFQVVSUnJyfv+1j8euyxx1yPIVy+LBeWC8vFsb+KK1KK\nyzac3gq+77cWZgdZ2Vm0nNCSkatHBjSwSLJ582a3QwgblguH5cJhuQicSFkWWwVUE5FEYDvQCehS\nmB3M+2Ye8zbOY/7G+QjCTZfeFPgojTHGAGE4cxGRCcAHQHUR2SIiN6lqNtAXWASsAyap6vrC7Ldd\nzXY81/A5FOWW2beQtiYt4LGHu169erkdQtiwXDgsFw7LReBIINbWwpmIqP/P+Ox7z/LgsgcRhJGt\nR9LL08u94IwxJkyJCBoDDf2A6X9df55p8AyKcvOsm0lfk+52SCGTmZnpdghhw3LhsFw4LBeBE3PF\nBeDB6x7k6fpPoyg3zboppgqMMcaEQswti/l7esXTPLL8EQQhrU0aPS/pGeLojDEmPNmyWDE8fP3D\nPJn8JIrSa2Yvxnw+xu2QjDEmKsR0cQF45IZHeCL5CRQldWYqY78Y63ZIQWPryQ7LhcNy4bBcBE7M\nFxeAR294lMeTHs8rMOO+GOd2SMYYE9FiuueS3xPvPMHAzIHESRyj24ym28XdghydMcaEJ+u5BNCA\negMYnDSYHM2h58yejF9rD7Y0xpiisOKSz8B6A/MKTI8ZPZiwdoLbIQWMrSc7LBcOy4XDchE4VlyO\nYGC9gQyqN4gczaH7jO5M/HKi2yEZY0xEsZ7LMQzKHMTgdwYTJ3GMazuOzv/qHODojDEmPFnPJYgG\nJQ1i4A0DydEcuk3vxqQvJ7kdkjHGRAQrLscxKGkQA24YkFdgJn812e2QiszWkx2WC4flwmG5CJxI\neZ6La0SEwUmDUVWeXPEkXad1BaBjrY4uR2aMMeHLei4FpKoMyBjAUyueIl7imdBuAh1qdQhAhMYY\nE36s5xIiIsITyU/w8HUPc1gP02VaF6atm+Z2WMYYE5asuBSCiPBk/SfzCkznaZ2Zvn6622EVmK0n\nOywXDsuFw3IROFZcCim3wDx03UNk52TTaWqniCowxhgTCtZzKSJV5eFlDzPk/SEkxCUwuf1kbqxx\nY8CPY4wxbrCei0tEhKcbPE3/a/uTnZNNx6kdmbF+htthGWNMWLDiUgwiwjMNnuGBax6IiAJj68kO\ny4XDcuGwXASOFZdiEhGGNBzC/13zf3kFZubXM90OyxhjXGU9lwBRVR5Y8gAvfPgCCXEJTO0wldYX\ntQ76cY0xJhhiuuciIheJyGsiMllEbnE5Fp5r9Bz9ru5Hdk42HaZ0YPaG2W6GZIwxrono4qKqX6vq\nHUBnoInb8YgIzzd6nn5X9+NQziHaT27PnA1z3A4rj60nOywXDsuFw3IROGFRXERkpIjsFJG1+ban\niMjXIrJRRPof5b0tgXlAWDx0JbfA3H/V/RzKOUS7ye3CqsAYY0wohEXPRUSuB/YBo1W1tm9bPLAB\naAhsAz4BugBXAJcBz6vqdr99zFLVfzQ5QtVzyU9V6be4H0M/GkqJuBJM7zSdFtVbhDwOY4wpiuL2\nXMLirsiqukJEEvNtrgNsUtXNACIyEWitqkOAMb5t9YC2QGkg42j793g8eDweEhMTqVChAh6Ph6Sk\nJMCZBgdj/N/G/+XHz39k2rpptJvcjmkdp3HC9hOCdjwb29jGNi7qODMzk7S0NHbs2EFWVhbFFRYz\nFwBfcZnjN3NpDzRR1T6+cXegrqreXcj9ujJzyaWq3L/ofoZ9PIyS8SWZ3nE6zas3dyWWzMzMvA9V\nrLNcOCwXDsuFI5rPFguPqldMIsKLTV7knrr3cPDwQdpObsv8jfPdDssYY4IqnGcuVwGDVDXFN34I\nyFHVZwu5X1dnLrlUlXsX3svwlcMpGV+SGZ1m0KxaM7fDMsaYI4rmmcsqoJqIJIpISaATELEXjogI\nw1KGcXeduzl4+CA3TrqRBRsXuB2WMcYERVgUFxGZAHwAVBeRLSJyk6pmA32BRcA6YJKqrnczzuIS\nEf6X8j/6Xtk3r8As3LQwZMfPbd4Zy4U/y4XDchE44XK2WJejbF8ARNWf9yLC8KbDUZRXPnmFNhPb\nMLPzTFIuSHE7NGOMCZiw6bkES7j0XPJTVfrO78urq16lVHwpZnWeRZMLXL/JgDHGANHdc4lqIsLL\nzV7mzivu5MDhA7Se2JrF3y52OyxjjAkIKy4uyi0wd1xxR0gKjK0nOywXDsuFw3IROFZcXJZbYG6/\n/HaysrNoPbE1S75d4nZYxhhTLNZzCRM5msOd8+7kjU/foHRCaeZ0mUPD8xq6HZYxJkZZzyVKxEkc\nrzZ/lVsvu5Ws7CxaTmjJ0u+Wuh2WMcYUiRWXMBIncbzW4jX6XNYnr8As+25ZwPZv68kOy4XDcuGw\nXASOFZcwEydxvN7i9b8VmOXfL3c7LGOMKRTruYSpHM3htjm3MWL1CMoklGFu17nUP7e+22EZY2KE\n9VyiVJzE8UbLN7jl0lvYn72fFuNbkPH9UR9ZY4wxYcWKSxiLkzjebPkmN3tuZn/2fpqPb07m5swi\n78/Wkx2WC4flwmG5CBwrLmEuTuJ4q9VbeQWm2bhmxSowxhgTCtZziRA5mkPv2b0ZtWYUZUuUZX7X\n+dRLrOd2WMaYKGU9lxgRJ3GMaDWCXp5e/HXoL5qNb8Y7m99xOyxjjDkiKy4RJE7iGNGy6AXG1pMd\nlguH5cJhuQgcKy4RJj4unhEtR5B6SWpegXn3h3fdDssYY/7Gei4R6nDOYW6efTOjPx9NuRLlWNBt\nAdefc73bYRljooT1XGJUfFw8I1uNpOclPfnz0J80HdeUFT+scDssY4wBrLhEtNwC0+PiHnkF5r0f\n3zvq62092WG5cFguHJaLwLHiEuHi4+IZ1XoU3S/unldg3v/xfbfDMsbEOOu5RInDOYfpNasXY78Y\nywklT2BR90VcU/Uat8MyxkQo67kYwDuDSWudRrfa3dh3cB8pY1P4YMsHbodljIlREV1cxOspERku\nIj3djsdt8XHxpLdJp2vtruw9uJeUsSl8uOXDvH+39WSH5cJhuXBYLgInoosL0AaoAhwEtrocS1jI\nLTBd/tWFvQf30mRsk78VGGOMCYWw6LmIyEigOfCzqtb2254CDAPigRGq+my+9/UHflfVt0Rkiqp2\nOMK+Y6Lnkl92TjY9Z/RkwpcTOLHkiSzusZirzrrK7bCMMREiWnouo4AU/w0iEg+87NteE+giIjVE\npIeIDBWRynhnK7t9b8kJZcDhLiEugdE3jqbzvzrnzWA+2vqR22EZY2JEWBQXVV0B7Mq3uQ6wSVU3\nq+ohYCLQWlXHqOp9qrodmA40EZHhQGZIg44ACXEJjLlxDJ1qdWLPgT00eLwBH2/92O2wwoKtrTss\nFw7LReAkuB3AMVQBtviNtwJ1/V+gqvuB3sfbkcfjwePxkJiYSIUKFfB4PCQlJQHOhylax++9+x59\nKvZBaymTv59M/cfr83yj57mzw51hEZ9b41zhEo+b4zVr1oRVPG6O16xZE1bxhHKcmZlJWloaO3bs\nICsri+IKi54LgIgkAnNyey4i0g5IUdU+vnF3oK6q3l3I/cZkzyW/7Jxsuk7rypR1UyhfqjxLeiyh\nTpU6bodljAlT0dJzOZJtQFW/cVXsjLAiS4hLYHy78XSo2YE9B/bQeExjVm5b6XZYxpgoFc7FZRVQ\nTUQSRaQk0AmY7XJMEe29d99jXNtxtK/Znj8O/EHjMY35ZNsnboflivzLY7HMcuGwXAROWBQXEZkA\nfABUF5EtInKTqmYDfYFFwDpgkqqudzPOaFAivgTj246nXY12/HHgDxqNacSq7avcDssYE2XCpucS\nLNZzObJDhw/RZVoXpq2fxkmlTmJpz6VcUfkKt8MyxoSJaO65mCAqEV+CCe0m0LZG27wZzKfbP3U7\nLGNMlLDiEkPyryeXiC/BxHYTufGiG9mdtZuGYxrGTIGxtXWH5cJhuQgcKy4xrkR8CSa2n0ibi9qw\nO2s3jcY04rOfPnM7LGNMhLOeiwHg4OGDdJzSkVkbZnFy6ZNZ2nMpl515mdthGWNcYj0XExAl40sy\nucNkWl/Yml1Zu2g4uiGrf1rtdljGmAhlxSWGHG89ObfAtLqwlbfAjIneAmNr6w7LhcNyEThWXMzf\nlIwvyZQOU2h1YSt+3/87Dcc0ZM2ONW6HZYyJMNZzMUd08PBB2k9uz5xv5lCxTEWW91zOJWdc4nZY\nxpgQsZ6LCYrcGUyL6i34ff/v1B9dn893fO52WMaYCGHFJYYUdj25VEIppnaYSvNqzfl9/+80GN0g\nagqMra07LBcOy0XgWHExx1QqoRTTOk6jebXm/Lb/NxqMbsAXO79wOyxjTJiznospkAPZB2g7uS3z\nN87n1LKnsrzncmqfXtvtsIwxQWI9FxMSuTOYphc05de/fqX+6Pqs3bnW7bCMMWHKiksMKe56cumE\n0kzvNJ2UC1LyCsyXP38ZmOBCzNbWHZYLh+UicKy4mEIpnVCaGZ1mOAUmPXILjDEmeKznYookKzuL\nNhPbsOjbRVQqW4mM1AxqnVbL7bCMMQFiPRfjitIJpZnZeSaNz2/ML3/9Qv3R9fnq56/cDssYEyas\nuMSQQK8nl04ozcxO3gLz858/U390fdb9si6gxwgWW1t3WC4clovAseJiiqVMiTLM7DSTRuc14uc/\nfyY5PTliCowxJnis52ICYv+h/bSe2Jol3y3h9HKnk5GaQY1KNdwOyxhTRNZzMWGhTIkyzOo8i4bn\nNWTnnztJTk9m/S/r3Q7LGOMSKy4xJNjrybkFpsG5DfIKzNe/fh3UYxaVra07LBcOy0XgRHRxEZGa\nIjJJRF4VkXZux2OgbImyzO4ym/rn1g/7AmOMCZ6I7rmIyP3ASlV9T0RmqWrrI7zGei4u+OvQX7Sc\n0JLl3y/njBPOIDM1kwtPvdDtsIwxBRQVPRcRGSkiO0Vkbb7tKSLytYhsFJH+R3jrGKCziDwHnBKS\nYE2BlC1Rljld5lD/3Prs2LeD5PRkNvy6we2wjDEhEhbFBRgFpPhvEJF44GXf9ppAFxGpISI9RGSo\niFRW1V9UtS/wEPBryKOOMKFeT84tMMmJyfy07yeS05P55rdvQhrD0djausNy4bBcBE5YFBdVXQHs\nyre5DrBJVTer6iFgItBaVceo6n2qul1EzhGRN4B04LkQh20KILfAJCUmhV2BMcYET4LbARxDFWCL\n33grUNf/Bar6A3Db8Xbk8XjweDwkJiZSoUIFPB4PSUlJgPOXSiyMk5KSXDv+3C5zaTGhBZkZmVwz\n4Bo+fOJDqp1SLazyE8vjXOESj1vj3G3hEk8ox5mZmaSlpbFjxw6ysrIorrBp6ItIIjBHVWv7xu2A\nFFXt4xt3B+qq6t2F3K819MPEnwf/pPn45rzzwztUPrEymamZVDulmtthGWOOICoa+kexDajqN66K\nd/Ziiij/X6mhVq5kOeZ1nccN59zA9r3bSU5PZtPvm1yJxe1chBPLhcNyETjhXFxWAdVEJFFESgKd\ngNkux2SKyb/AbNu7jaS0JNcKjDEmeMJiWUxEJgD18J5O/DMwUFVHiUhTYBgQD7ytqs8UYd+2LBaG\n9h3cR7NxzVjx4wqqnFiFzF6ZXFDxArfDMsb4FHdZLCyKSzBZcQlf/gXmrPJnkZmayfkVz3c7LGMM\n0d1zMQEWbuvJJ5Q8gfnd5nPd2dexdc9WktKT+Pb3b0Ny7HDLhZssFw7LReBYcTGuOqHkCczvOp9r\nq17L1j1bSU5P5rtd37kdljGmmGxZzISFvQf20nRcU97f8j5Vy1cls1cm5518ntthGROzbFnMRIUT\nS53Igm4LuKbqNWzZs4WktCSbwRgTway4xJBwX0/OLTBXn3U1W/ZsITk9me93fR+UY4V7LkLJcuGw\nXASOFRcTVsqXKs/C7gu5+qyr+fGPH0lKTwpagTHGBI/1XExY2nNgD03GNuGjrR9xzknnkNkrk8QK\niW6HZUzMsJ6LiUrlS5VnYbeF1K1Slx/++IGktCQ2797sdljGmAKy4hJDIm09+aTSJ7Go+6K8ApOc\nnswPu38IyL4jLRfBZLlwWC4Cx4qLCWu5BaZOlTps3r2ZpPSkgBUYY0zwWM/FRIQ/sv6g8djGrNy2\nknMrnEtmr0zOPulst8MyJmqF5N5iItIVKHGMlxxS1fFFDSKYrLhEj91Zu2k8pjGfbP/ECowxQWY3\nrjwOKy4O/yfsRSr/AnPeyeeRkZpRpAITDbkIFMuFw3LhCOnZYiLST0SWichXIvK0iBxrNmNMwFUo\nXYHFPRZzReUr+G7XdySnJ7Pljy3Hf6MxJqQKNXMRkRaqOldEBKgP1FPVgUGLLgBs5hKddu3fRaMx\njfj0p0857+TzyEzNpOpJVY//RmNMgYT6OpczRKQZUE5VlwGfFPXAxhTHyWVOZkmPJVx+5uV5M5it\ne+wp2MaEi8IWl6pATWCUiCwH/k9EuotI/8CHZgIt2s7hzy0wl515Gd/u+paktKQCF5hoy0VxWC4c\nlovAKWxxmQV8qKodVLU+cBMgQLOAR2ZMAeQWmEvPuJRvd31Lcnoy2/ZsczssY2JeQM4WE5EzVfWn\nAMQTcNZziQ2/7/+dhqMbsnrHai6oeAGZqZlUKV/F7bCMiViu3FtMRNr7j8O1sJjYUbFMRZb2XMql\nZ1zKpt832QzGGJcV9fYv5QIahQmJaF9Pzi0wnjM8bPx9I8npyWzfu/2Ir432XBSG5cJhuQgcu7eY\niSoVy1RkaY+CFRhjTPAUqeciIqmqmh6EeALOei6x6be/fqPhmIas2bGG6qdUJyM1g8onVnY7LGMi\nRkw9z0VEzhWRESIyxTduLSJvishEEWnkdnwmfJxS9hSW9ljKJadfwje/fUNyejI/7bXWoDGhUtTi\n4spUQFW/V9XefuNZqnorcDvQyY2YIkmsrSefUvYUlvZcysWnX/yPAhNruTgWy4XDchE4RS0uU4tz\nUBEZKSI7RWRtvu0pIvK1iGws5IWZjwIvFycmE51OLXsqy3ou4+LTL2bDbxuoP7o+O/btcDssY6Ke\nK3dFFpHrgX3AaFWt7dsWD2wAGgLb8N5apgtwBXAZ8Lyqbve9doqqdvDd42wIsNh3O5ojHct6LoZf\n//qV+un1WfvzWmqcWoOM1AxOP+F0t8MyJmxFZM9FVVcAu/JtrgNsUtXNqnoImAi0VtUxqnqfqm4X\nkYoi8jrgEZEHgb5AA6C9iNwW0h/CRJTcGUzt02qz/tf1JKcns3PfTrfDMiZqJRTnzSLSDfgL2K+q\nC4sZSxXA/97pW4G6/i9Q1d/x9lf8vXS8HXs8HjweD4mJiVSoUAGPx5P3zIbcNdZYGPuvJ4dDPKEe\nVypXiccTH+e+b+5j/ffrqS/1eeLcJ6hYpmJYxOfWeM2aNdx7771hE4+b42HDhsX074e0tDR27NhB\nVlYWxVWsZTERuRQ4BFxe2FOTRSQRmOO3LNYOSFHVPr5xd6Cuqt5d5ACxZTF/mfYgJAB++fMX6jxa\nh80VNlOzUk0yUjM4rdxpboflGvtcOCwXDreXxW7EO7vILOZ+wNtn8X8gR1W8sxcTIPYfjVelcpX4\n+MmPqVWpFut+WUf99Pr8/OfPboflGvtcOCwXgVPc4rIAWEq+5asiWgVUE5FEESmJ99Ti2QHYrzH/\ncFq501ieupxalWrx1S9fxXyBMSbQivWYY2CVqv6gqpMLuZ8JwAdAdRHZIiI3qWo23gb9ImAdMElV\n1xdmv+bYctdXjTcXuQWmZqWafPXLVzQY3YBf/vzF7dBCzj4XDstF4BR25rJBVRsA/wKWAQOKclBV\n7aKqlVW1lKpWVdVRvu0LVPVCVb1AVZ8pyr6NKYzTyp3G8p7eAvPlz19Sf3T9mCwwxgRaoRr6ItIb\n2A68q6r7RKSlqs4JWnQBYA19UxA79+0kOT2Z9b+u51+n/YvlPZdTqVwlt8MyxjXFbegXtrgMBvbi\n7bGcgvdU5jeBKqr6bFGDCCYrLqagduzbQf30+qz/dT21T6vNsp7LrMCYmBXqs8XsMccRzNaTHUfK\nxRknnMHy1OVcdOpFrP15LQ1GN+DXv34NfXAhZp8Lh+UicApVXFT1M1V932/8raqOwW4aaaLEGSec\nwfKesVdgjAm0Ai2LiUhXoMQxXnJIVccHLKoAsmUxUxQ/7f2J5PRkNvy2gUtOv4RlPZdxStlT3A7L\nmJAJac8lEllxMUVlBcbEMrev0DcRxNaTHQXJxZknnklGagbVT6nO5zs/p+GYhvz212/BDy7E7HPh\nsFwEjhUXY47Bv8Cs2bEmaguMMYFW0J6LAGep6pbjvjjM2LKYCYTte7eTlJbExt83cukZl7K051Iq\nlqnodljGBE0ol8UWFPUgxkS6yidWJiM1g2oVq7F6x2oajm7I7/t/dzssY8JWgYqL70//T0WkTpDj\nMUFk68mOouSiSvkqZKRmcEHFC1i9YzWNxjSKigJjnwuH5SJwCjNzuQr4UES+E5G1vq8vghWYMeHI\nv8B89tNnNBrTiF378z9U1RhT4FORfQ/3Ash9gwCo6uZABxVI1nMxwbB1z1aS0pL4dte3XH7m5Szp\nsYSTy5zsdljGBEyo7y3mAa7HW2BWqOrnRT1wqFhxMcHiX2CuqHwFi7svtgJjokbIGvoicg8wFqgE\nnA6MFZF/F/XAJvRsPdkRiFycVf4sMlIzOO/k81i1fRWNxzZmd9bu4gcXYva5cFguAqcwPZfeeJ9p\nP1BVB+DtwfQJTljGRIaqJ1UlMzXTKTBjIrPAGBNohem5rAXqqOp+37gMsFJVawcxvmKzZTETCj/+\n8SNJaUl8v/t7rqx8JYt7LKZC6Qpuh2VMkYXyOpdRwMciMsj3XJePgJFFPbAx0eTsk84ms1cm51Y4\nl0+2f0KTsU34I+sPt8MyxjUFKi6+K/Sn4n1+yy7gN6CXqg4NYmwmwGw92RGMXJx90tlkpGaQWCGR\nldtW0nhs44goMPa5cFguAqcwM5f5qvqpqv5PVYer6uqgRWVMhDqnwjlkpmbmFRibwZhYVZieSzrw\niqquDG5IgWU9F+OGzbs3k5SWxA9//MBVZ13Fou6LKF+qvNthGVNgIbvORUQ2ABcAPwB/+jarql5c\n1IOHghUX4xYrMCaShaSh7+u59AHOB+oDLX1frYp64KIQkXNFZISITPGNk0RkhYi8JiL1QhlLJLL1\nZEcocpFYIZHMXpmcc9I5fLT1I1LGprDnwJ6gH7ew7HPhsFwETmF6Lq+q6ub8X8EK7EhU9XtV7e23\nKQfYC5QCtoYyFmMKIrfAnH3S2Xy49UOajmsalgXGmEBzpeciIiOB5sDP/tfJiEgKMAyIB0ao6rNH\nef8UVe0gvjUvETkNeFFVux/htbYsZlz3/a7vSUpP4sc/fuSaqtewsNtCTix1otthGXNUobzOJZB3\nRR4FpPhvEJF44GXf9ppAFxGpISI9RGSoiFTOvxO/qrEb7+zFmLB07snnkpGaQdXyVflgywekjEth\n74G9bodlTNAUprg0IUA9F1Vdgfd6GX91gE2+5bZDwESgtaqOUdX7VHW7iFQUkdcBj4g8KCI3+saj\ngZeKEksssfVkhxu5OO/k88jslZlXYJqOaxoWBcY+Fw7LReAkFPSFIeivVAH8H6O8FaibL4bfgdvz\nvW/G8Xbs8XjweDwkJiZSoUIFPB4PSUlJgPNhsnFsjXOF+vg/fv4jQy4YwoObHuT9Le9zzcBreLbB\nszRr3My1fKxZs8b1/z/CZbxmzZqwiieU48zMTNLS0tixYwdZWVkU13F7LiLygKo+5/u+g6pO8fu3\np1X14SId2Pt8mDm5PRcRaQekqGof37g73htl3l2U/fsdx3ouJux8+/u3JKUnsXXPVq47+zoWdFvA\nCSVPcDssY/KEoufSxe/7/IWkaVEPfATbgKp+46rYGWAmSp1f8XwyUzM5q/xZvPfjezQb14x9B/e5\nHZYxAVOYnkuwrQKqiUiiiJQEOgGzXY4pquRfEopl4ZCL8yueT0ZqBlVOrMKKH1fQfHxzVwpMOOQi\nXFguAseV4iIiE4APgOoiskVEblLVbKAvsAhYB0xS1fVuxGdMqFxQ8QIye2VS5cQqvPvDuzQf35w/\nD/55/DcaE+YK0nM5DPzlG5YB9vv9cxlVLfBJAW6wnouJBJt+30RSWhLb9m6j3jn1mNd1HuVKlnM7\nLBPDQnZvsUhlxcVEio2/bSQpPYnte7eTlJjE3C5zrcAY14TyIkoT4Ww92RGOuah2SjUyUzOpfGJl\nMjdn0nJCy5AskYVjLtxiuQgcKy7GhJFqp1QjIzWDM084k4zNGbSc0JK/Dv11/DcaE2ZsWcyYMPTN\nb9+QlJbET/t+Ijkxmbld51K2RFm3wzIxJGTLYiIS57vP10Df+GwRqVPUAxtjjq76KdVtBmMiWqFu\nuQ9cDXT1jff5tpkIYevJjkjIxYWnXkhGagZnnHAGy79fTqsJrYJSYCIhF6FiuQicwhSXuqp6J75T\nkX33+SoRlKiMMcDfC8yy75cFrcAYE2iFeZ7Lx8A1wCpVvVREKgGLVfXSYAZYXNZzMdHg61+/Jikt\niZ1/7qRjFWRyAAAgAElEQVTheQ2Z3Xk2ZUqUcTssE8VCeSryS3jvQHyaiDwNvA88U9QDG2MK7qJT\nLyIjNYPTy53O0u+W0mpiK/Yf2n/8NxrjkgIXF1UdC/THW1C2433WyuRgBWYCz9aTHZGYixqVavyt\nwLSe2DogBSYScxEslovAKczZYs+q6npVfdn3tV5EjvgYYmNMcNSoVIPlqcs5rdxpLPluCW0mtbEZ\njAlLhem5rM7fXxGRtbnPYwlX1nMx0WjdL+tITk/m5z9/pvH5jZnZaab1YExABb3nIiJ3iMha4EIR\nWev3tRn4oqgHNsYUXc1KNVneczmVylZi8beLuXHSjWRlF//pgcYESkGWxcYDLYFZQAu/r8tUtVsQ\nYzMBZuvJjmjIRa3TapGRmkGlspVY9O0i2kxsU6QCEw25CBTLReAct7io6h+quhn4Gujl99U392p9\nY4w7ap1Wi+Wpy/MKjM1gTLgoTM/lP0Dui8vgnb2sU9WbgxRbQFjPxcSCL3/+kuT0ZH7961eaXtCU\n6Z2mUzqhtNthmQjm2vNcRKQU3oso6xX14KFgxcXEirU711J/dH1+/etXmlVrxvSO0ymVUMrtsEyE\ncvN5LuWAKsV4vwkxW092RGMuap9em2U9l3Fq2VOZv3E+bSe35UD2geO+LxpzUVSWi8ApzHUu/meK\nfQVsAP4XvNCMMYV18ekXs6znMk4pcwrzN86n3eR2BSowxgRaYXouiX7DbGCnqh4KQkwBZctiJhZ9\nsfML6qfX57f9v9GiegumdphqS2SmUFzruUQKKy4mVn2+43MajG5gBcYUSSguotwnInuP8rWnqAc2\noWfryY5YyMUlZ1zCsp7LqFimInO/mUuHKR2OuEQWC7koKMtF4BTkOpcTVPXEo3yVD0WQxpii8S8w\nc76ZQ4cpHTh4+KDbYZkYUKhlMRG5BLgB7/UuK1T182AFdpTjnws8Apykqh1E5CxgOLAL+EZV/3Ej\nTVsWMwZW/7SaBqMbsCtrF60ubMWUDlMoGV/S7bBMGAvZqcgicg8wDqgEnA6MFZF/F/XARaGq36tq\nb79NtYFpqnoLENYPLTPGTZeeeSnLei7j5NInM3vDbDpO6WgzGBNUhbnOpTfeRx0PVNUBwFVAn6Ic\nVERGishO3w0x/beniMjXIrJRRPoXYFcfALeKyDJgYVFiiSW2nuyIxVxceualLO25lJNLn8ysDbPo\nNLUTBw8fjMlcHI3lInAKexFlzlG+L6xRQIr/BhGJB172ba8JdBGRGiLSQ0SGikjlI+znJuBRVW0A\nNC9GPMbEhMvOvIwlPZZQoXQFZn49k85TO5N9ONvtsEwUKsx1LvfjvWHldECANkCaqg4t0oG9183M\nyX0ejIhcDTymqim+8YMAqjrE7z0VgaeBBsAIYAEwEPgF2KuqDxzhONZzMSafT7d/SsMxDdmdtZsb\nL7qRSe0nUSK+hNthmTBS3J5LQkFfqKovisg7wHV4G/q9VHV1UQ98BFWALX7jrUDdfDH8Dtye733t\nj7djj8eDx+MhMTGRChUq4PF4SEpKApxpsI1tHGvjJT2WkDQoiRnfz6CzdGZiu4m8v+L9sInPxqEd\nZ2ZmkpaWxo4dO8jKKv6dtQszc+kALFLVPSIyALgMeEJVPyvSgf85c2kHpKhqH9+4O94ez91F2b/f\ncWzm4pOZmZn3oYp1lguvVdtXkTQoiT+r/Em7Gu2Y0G5CTM9g7HPhCOWNKwf6Cst1eJel3gZeL+qB\nj2AbUNVvXBXv7MUYEyRXVL6CFxq9wEmlTmLa+ml0mdaFQ4fD/q5OJgIUZuayRlU9IjIEWKuq40Rk\ntaoW6RTgI8xcEvDeDLMBsB1YCXRR1fVF2b/fcWzmYsxxrNy2ksZjGvPHgT9oX7M949uOj+kZjAnt\nzGWbiLwJdALmiUjpQr4/j4hMwHsacXUR2SIiN6lqNtAXWASsAyYVt7AYYwqmTpU6LO6xmPKlyjN1\n3VS6Te9mMxhTLIUpDh3x/uJvrKq7gZOB/yvKQVW1i6pWVtVSqlpVVUf5ti9Q1QtV9QJVfaYo+zZH\nl9u8M5YLf7m5qFOlDkt6LKF8qfJMWTeFbtO7kZ0TW6cp2+cicApcXFT1T1WdpqobfeOfVHVx8EIz\nxoRanSp1WNx9cUwXGBMYhem5lAHuxDkVeQXwmqoW/5y1ILKeizGF99HWj2gytgl7DuyhU61OjG07\nloS4Al+5YKJAyJ7nIiJTgD3AWLwXUXbFdwPJoh48FKy4GFM0H239iMZjGrP34F46/6szY24cYwUm\nhoSyoV9LVW9R1QxVXe67gWStoh7YhJ6tJzssF46j5eKqs65icY/FnFjyRCZ+OZGeM3pG/RKZfS4C\npzDF5TPfLVoAEJGrgE8DH5IxJlxcddZVLOq+iBNLnsiELyeQOjM16guMCYzjLov53bk4AbgQ7y1a\nFDgb+FpVawY1wmKyZTFjiu+DLR/QZGwT9h3cR9faXUlvk25LZFEu6D0X38WO/hRvz+Vs4EFVbVbU\ng4eCFRdjAiN/gRndZjTxcfFuh2WCJOg9F1XdnPsFVMR7oWMm8DjeuxKbCGHryQ7LhaOgubim6jUs\n7LaQE0qewPi140mdmcrhnMPBDS7E7HMROMctLiJyoYgMEpH1wDDgRyBOVZNU9aWgR2iMCRvXnn0t\nC7otoFyJcoxbO45es3pFXYExgVGQZbEcYC7QV1V/9G37XlXPDUF8xWbLYsYE3ns/vkfK2BT+PPQn\nPS7uwajWo2yJLMqE4lTktsB+4F0ReV1EGuDtuRhjYtR1Z1+XN4MZ88UYbpp1k81gzN8UpOcyU1U7\nAf/Ce1X+fUAlEXlNRBoHO0ATOLae7LBcOIqai+vPuZ753ebnFZibZ98c8QXGPheBU5h7i+1T1XGq\n2gLvs1ZWAw8GLTJjTNi74Zwb8grM6M9Hc8vsWyK+wJjAKPDtXyKV9VyMCb53f3iXpuOa8tehv+jl\n6cWIliOsBxPhQnn7F2OMOaIbzrmB+V3nU7ZEWdLWpNF7Tm+bwcQ4Ky4xxNaTHZYLR6ByUS+xHvO6\nzssrMH3m9CFHcwKy71Cxz0XgWHExxgRMUmIS87rOo0xCGUatGUXv2b0jrsCYwLCeizEm4DK+z6D5\n+Obsz97PzZ6beavVW8SJ/S0bSaznYowJO8nnJjO361zKJJRh5JqR3DrnVpvBxBgrLjHE1pMdlgtH\nsHJR/9z6zOkyh9IJpXl79dvcNue2sC8w9rkIHCsuxpigaXBeA+Z2mUvphNKMWD2C2+feHvYFxgSG\n9VyMMUG37LtltJjQgqzsLG697FZea/Ga9WDCXNCf5xJORKQ10BwoD7wNbAIeAU5S1Q5HeY8VF2PC\nwNLvltJyQkuysrO47fLbeLX5q1ZgwlhMNfRVdZaq3grcDnRS1e9VtbfbcUUKW092WC4cocpFw/Ma\nMrvzbEonlOaNT9/grnl3hd0SmX0uAseV4iIiI0Vkp98jlHO3p4jI1yKyUUT6H2MXjwIvBzdKY0yg\nNTq/EbM6z6JUfCle//T1sCwwJjBcWRYTkeuBfcBoVa3t2xYPbAAaAtuAT4AuwBXAZcDzwE/AEGCx\nqi7z298UWxYzJnIs/nYxrSa04sDhA9xxxR280uwVROxJHuEkIpfFVHUFsCvf5jrAJt8jlQ8BE4HW\nqjpGVe9T1e3A3UADoL2I3CYiFUXkdcBznJmOMSaMND6/cd4M5rVVr3HX/LuwPwKjS4LbAfipAmzx\nG28F6vq/QFWHA8Pzve/24+3Y4/Hg8XhITEykQoUKeDwekpKSAGeNNRbG/uvJ4RCPm+PcbeESj5vj\nNWvWcO+994b8+E0uaMLj5z7Oo8sf5TVeQxDal22PiLiWj2HDhsX074e0tDR27NhBVlYWxeXa2WIi\nkgjM8VsWawekqGof37g7UFdV7y7mcWxZzCczMzPvQxXrLBcOt3OxcNNCWk9szcHDB7nryrt4qelL\nri2RuZ2LcBKxpyIfobhcBQxS1RTf+CEgR1WfLeZxrLgYE+YWbFxAm0ltOHj4IH2v7MvwpsOtB+Oy\niOy5HMUqoJqIJIpISaATMNvlmIwxIdC0WlNmdJpByfiSvPzJy9yz8B7rwUQ4t05FngB8AFQXkS0i\ncpOqZgN9gUXAOmCSqq53I75o5d9viHWWC0e45KJZtWZ5BeallS9x78J7Q15gwiUX0cCVhr6qdjnK\n9gXAghCHY4wJE82qNWN6x+m0ndyW4Su9S2NDmwy1JbIIFFG3fykK67kYE3nmfjOXtpPacijnEPfU\nvccKjAuiqedijDEAtKjegumdplMirgT/+/h/3L/ofuvBRBgrLjHE1pMdlgtHuObCv8AM+3gY/Rb3\nC3qBCddcRCIrLsaYsNWiegumdZxGibgSDP1oKP9Z/B+bwUQI67kYY8Le7A2zaT+5PYdyDtHv6n48\n3+h568EEmfVcjDFRr9WFrZjSYQol4krw3w//ywNLHrAZTJiz4hJDbD3ZYblwREouWl/UmikdppAQ\nl8ALH75A/6X9A15gIiUXkcCKizEmYvgXmOc/eJ4Hlz5oM5gwZT0XY0zEmbF+Bh2ndiQ7J5sHr32Q\npxs8bT2YALOeizEm5txY40YmtZ9EQlwCQ94fwsPLHrYZTJix4hJDbD3ZYblwRGou2tZoy8R2E/MK\nzCPLHyl2gYnUXIQjKy7GmIjVrmY7JrabSLzE88x7z/Do8kdtBhMmrOdijIl409ZNo9PUThzWwzxy\n/SM8kfyE9WCKyXouxpiY165mOya2985gnlrxFAMyBtgMxmVWXGKIrSc7LBeOaMlF+5rtmdBuQl6B\nGZgxsNAFJlpyEQ5ceZ6LMcYEQ4daHVCUrtO68uSKJxERBicNtiUyF1jPxRgTdSZ/NZmu07pyWA8z\n8IaBDE4e7HZIEae4PRebuRhjok7HWh1RVbpN78bj7z6OiDAoaZDbYcUU67nEEFtPdlguHNGai07/\n6sS4tuOIkzgGvzOYQZmDjvueaM2FG2zmYoyJWp3+1QnFO4MZ/M5gBOGxpMfcDismWM/FGBP1Jqyd\nQPcZ3cnRHAYnDWZgvYFuhxT2rOdijDHH0aV2FxSlx4wePJb5GIIwoN4At8OKatZziSG2nuywXDhi\nJRdda3dldJvRxEkcAzMH8uS7T/7jNbGSi1CIqOIiIq1F5E0RmSgijXzbyonIJyLS3O34jDHhrdvF\n3Uhvk06cxDEgYwBPvfuU2yFFrYjsuYhIBeAFVe0tIoOBvcB6VZ13hNdaz8UY8zdjvxhL6sxUcjSH\np+o/xcPXP+x2SGEnIu8tJiIjRWSniKzNtz1FRL4WkY0i0v8Yu3gUeNk3e1kH/BLMeI0x0aX7xd1J\nb5OOIDyy/BGeWfGM2yFFHbeWxUYBKf4bRCQeeNm3vSbQRURqiEgPERkqIpXF61lggaquAeoBVwFd\ngT5i93g4JltPdlguHLGaC/8C8/Dyhxny3pCYzUUwuHK2mKquEJHEfJvrAJtUdTOAiEwEWqvqEGCM\nb9u/gQZAeRG5QFUf9W1PBX6x9S9jTGH0uKQHitJrZi8eWvYQfSr2ISkpye2wokI4nYpcBdjiN94K\n1PV/gaoOB4bnf6Oqph9rxx6PB4/HQ2JiIhUqVMDj8eR9gHL/UomFcVJSUljFY+PwGecKl3hCOT6b\nsxnVehQ3zbqJtz59C16CN+9+M2ziC9U4MzOTtLQ0duzYQVZWFsXlWkPfN3OZo6q1feN2QIqq9vGN\nuwN1VfXuYh7HJjTGmONKW5PGzbNuRlGebfgsD1z7gNshuSoiG/pHsQ2o6jeuinf2YgIk/1+pscxy\n4bBcePXy9OL/Kv8fgtB/aX+ef/95t0OKaOFUXFYB1UQkUURKAp2A2S7HZIyJIU2rNWVEqxEIwgNL\nH+CFD15wO6SI5cqymIhMwHum1ynAz8BAVR0lIk2BYUA88LaqFvv8QFsWM8YU1sjVI7ll9i0AvNDo\nBfpd08/liEKvuMtiEXkRZWFYcTHGFMXbn71N7zm9Afhv4/9y/9X3uxxRaEVTz8UEma2tOywXDsuF\nwz8Xt1x2C2+1fAuAfov78eKHL7oUVWSy4mKMMUfR+7LefyswQz8c6nJEkcOWxYwx5jje+vQtbp17\nKwBDmwzl3qvudTmi4LNlMWOMCbI+l/fhjRZvAHDfovsY9tEwlyMKf1ZcYoitrTssFw7LheNYubj1\n8lt5vfnrgLfA/O+j/4UoqshkxcUYYwrotitu47XmrwFw76J7Gf7xP+5GZXys52KMMYX0+qrXuWPe\nHQAMTxnO3XWLdZeqsGQ9F2OMCbHbr7idV5u9CsC/F/6bl1e+7HJE4ceKSwyxtXWH5cJhuXAUJhd3\nXHkHrzR7BYC7F9zNKytfCVJUkcmKizHGFNGdV97Jy029s5a+C/ry6ievuhxR+LCeizHGFNMrK1+h\n74K+3u+bvcKdV97pckTFZz0XY4xx2V117sqbwdw1/y5eX/W6yxG5z4pLDLG1dYflwmG5cBQnF3fV\nuYuXmr4EwB3z7oj5AmPFxRhjAqRvnb4MT/Fe+3LHvDt4Y9UbLkfkHuu5GGNMgA3/eDj3LLwHgNeb\nv85tV9zmckSFZz0XY4wJM/+u+2+GNfHef+z2ebfz5qdvuhxR6FlxiSG2tu6wXDgsF45A5uKeq+5h\naBPvLfpvm3sbb336VsD2HQmsuBhjTJDce9W9vNjY+5CxW+feyojPRrgcUehYz8UYY4LsxQ9fpN/i\nfgCMaDmCWy67xeWIjs96LsYYE+buv/p+Xmj0AgB95vRh5OqRLkcUfFZcYoitrTssFw7LhSOYueh3\nTT+eb/Q8itJ7du+oLzBWXIwxJkT+c81/eK7hc3kFZtTqUW6HFDQR1XMRkdZAc6A88DawH+gGJAA1\nVfXaI7zHei7GmLDy/PvP88DSBxCEka1H0svTy+2Q/qG4PZeIKi65RKQC8IKq9vaNWwOnqeo/zvWz\n4mKMCUfPvf8c/Zf2D9sCE5ENfREZKSI7RWRtvu0pIvK1iGwUkf7H2MWjgP/TeboC44MRazSxtXWH\n5cJhuXCEMhcPXPsAQxoMQVFunnUzaWvSQnbsUHCr5zIKSPHfICLxeAtGClAT6CIiNUSkh4gMFZHK\n4vUssEBV1/jedzbwh6r+GeKfwRhjiqX/df15psEzeQUmfU262yEFjGvLYiKSCMxR1dq+8dXAY6qa\n4hs/CKCqQ/ze82+gJ/AJsEZV3xCRQcBCVf3oKMexZTFjTFh7ZsUzPLz8YQQhrU0aPS/p6XZIxV4W\nSwhkMMVUBdjiN94K1PV/gaoOB4bn2zboeDv2eDx4PB4SExOpUKECHo+HpKQkwJkG29jGNraxW+OH\nkh5CUR55+xFSh6Yi9wk9LukR0ngyMzNJS0tjx44dZGVlUVzhNHNpB6Soah/fuDtQV1XvLuZxbObi\nk5mZmfehinWWC4flwuF2Lp569ykezXgUQRh942i6X9zdtVgisqF/FNuAqn7jqnhnL8YYExMeueER\nnkx+EkVJnZnKuC/GuR1SkYXTzCUB2AA0ALYDK4Euqrq+mMexmYsxJqI8+e6TDMgYQJzEMbrNaLpd\n3C3kMUTkzEVEJgAfANVFZIuI3KSq2UBfYBGwDphU3MJijDGR6NEbHuXxpMfJ0Rx6zuzJ+LWRd6WF\nK8VFVbuoamVVLaWqVVV1lG/7AlW9UFUvUNVn3IgtmuU274zlwp/lwhFOuRhQbwCDkwaTozn0mNEj\n4gpMOPVcjDHG+BlYbyCD6g3KKzAT1k5wO6QCi8jbvxSG9VyMMZFuUOYgBr8zmDiJY1zbcXT+V+eg\nHzMiey7GGGMK7rF6jzHghgHkaA7dpndj0peT3A7puKy4xJBwWk92m+XCYblwhGsuRITBSYN59PpH\n8wrM5K8mux3WMYXTFfrGGGOOQkR4PPlxAJ5c8SRdp3UFoGOtjm6GdVTWczHGmAiiqgzIGMBTK54i\nXuKZ0G4CHWp1CPhxouneYsYYY45DRHgi+QlUlaffe5ou07oABKXAFIf1XGJIuK4nu8Fy4bBcOCIl\nFyLCk/Wf5KHrHuKwHqbLtC5MXTfV7bD+xoqLMcZEIBHhqfpP8eC1D3JYD9N5amemrZvmdlh5rOdi\njDERTFV5aNlDPPv+syTEJTCp/STa1mhb7P3adS7GGBPDRIRnGjxD/2v7k52TTaepnZixfobbYVlx\niSWRsp4cCpYLh+XCEam5yC0wD1zzANk52XSc2tH1AmPFxRhjooCIMKThkL8VmJlfz3QvnmjvR1jP\nxRgTS1SVB5Y8wAsfvkBCXAJTO0yl9UWtC70f67kYY4zJIyI81+g5/nP1f8jOyabDlA7M+npWyOOw\n4hJDInU9ORgsFw7LhSNacpFbYPpd3Y9DOYfoMKUDczbMCWkMVlyMMSYKiQjPN3qe+6+6n0M5h2g3\nuV1IC4z1XIwxJoqpKv0W92PoR0MpEVeCaR2n0fLClsd9n/VcjDHGHJWI8N/G/+XeuvfmzWDmfjM3\n6Me14hJDomU9ORAsFw7LhSNacyEivNjkxb8VmHnfzAvqMa24GGNMDMgtMPfUvYeDhw/SdnJb5m+c\nH7zjRXs/wnouxhjjUFXuXXgvw1cOp2R8SWZ0mkGzas3+8TrruRhjjCkwEWFYyjDurnM3Bw8f5MZJ\nN7Jg44KAHyeiiouIXCQir4nIZBG5RUTKiki6iLwpIl3dji/cRet6clFYLhyWC0es5EJE+F/K/+h7\nZd+8ArNw08KAHiOiiouqfq2qdwCdgSZAW2Cyqt4KtHI1uAiwZs0at0MIG5YLh+XCEUu5EBGGNx1O\n3yv7cuDwAdpMbMOiTYsCtn9XiouIjBSRnSKyNt/2FBH5WkQ2ikj/o7y3JTAPmAicBWz1/dPhoAYd\nBXbv3u12CGHDcuGwXDhiLRe5BeauK+/yFphJbdjyx5aA7DshIHspvFHAS8Do3A0iEg+8DDQEtgGf\niMhs4ArgMuB5Vd2uqnOAOSIyC5iMt8B8ToTNwowxJhyICC81fQmA6qdUp+pJVQOyX1eKi6quEJHE\nfJvrAJtUdTOAiEwEWqvqEGCMb1s9vEthpYEMYAbwsog0B2aHJPgItnnzZrdDCBuWC4flwhGrucgt\nMCJFPjnsn/t06zRdX3GZo6q1feP2QBNV7eMbdwfqqurdxTyOnYdsjDFFUJxTkd1aFjuSoBSB4iTH\nGGNM0YRTn2Ib4L/YVxWnWW+MMSaChFNxWQVUE5FEESkJdML6KMYYE5HcOhV5AvABUF1EtojITaqa\nDfQFFgHrgEmqur6Yxznuqc3RSkSqikiGiHwlIl+KyL992yuKyBIR+UZEFotIBbdjDRURiReR1SIy\nxzeOyVyISAURmSoi60VknYjUjeFcPOT7b2StiIwXkVKxkosjXRJyrJ/dl6uNvt+pjY+7/2i975bv\n1OYN+J3aDHQpbsGKFCJyBnCGqq4RkROAT4E2wE3Ar6r6nK/gnqyqD7oZa6iIyP3A5cCJqtpKRJ4j\nBnMhIunAO6o6UkQSgHLAI8RYLnwnFS0HaqjqARGZBMwHahEDuRCR64F9wGi/E6uO+N+EiNQExgNX\nAlWApUB1Vc052v7DaVks0PJObVbVQ3gvumztckwho6o7VHWN7/t9wHq8H4pWQLrvZel4C07UE5Gz\ngGbACCD3JI+Yy4WInARcr6ojAVQ1W1X/IAZzAewBDgFlfUW2LLCdGMmFqq4AduXbfLSfvTUwQVUP\n+S4X2YT3d+xRRXNxqQL4X2q61bct5vj+QrsU+Bg4XVV3+v5pJ3C6S2GF2lDg/wD/v7RiMRfnAr+I\nyCgR+UxE3hKRcsRgLlT1d+C/wI94i8puVV1CDObCz9F+9sr8/QSr4/4+jebiEp3rfYXkWxKbBtyj\nqnv9/833LIKoz5OItAB+VtXVOLOWv4mVXOC9/OAy4FVVvQz4E/jbkk+s5EJEzgfuBRLx/vI8wXd9\nXZ5YycWRFOBnP2Zeorm4xPypzSJSAm9hGaOqM32bd/r6MYjImcDPbsUXQtcArUTke2ACUF9ExhCb\nudgKbFXVT3zjqXiLzY4YzMUVwAeq+pvvhKLpwNXEZi5yHe2/ify/T8/ybTuqaC4uMX1qs3jv4/A2\nsE5Vh/n902wg1fd9KjAz/3ujjao+rKpVVfVcvHfUXq6qPYjNXOwAtohIdd+mhsBXwBxiLBfA18BV\nIlLG999LQ7xnqsZiLnId7b+J2UBnESkpIucC1YCVx9pR1J4tBiAiTYFhQDzwtqo+43JIISMi1wHv\nAl/gTF8fwvuBmAycDWwGOqpqzNwK1nd/un6+s8UqEoO5EJFL8J7YUBL4Fu8ZhPHEZi4ewPtLNAf4\nDOgNnEgM5MJ3SUg94FS8/ZWBQO4Ngf/xs4vIw8DNQDbeZfZj3p8/qouLMcYYd0TzspgxxhiXWHEx\nxhgTcFZcjDHGBJwVF2OMMQFnxcUYY0zAWXExxhgTcFZcTMQRkRwRecFv/B8ReSxA+04TkXaB2Ndx\njtPBd7v7Zfm2J+beAl1ELvFdqxWoY54kInf4jSuLyJRA7d8Yf1ZcTCQ6CNwoIqf4xoG8WKvI+/Ld\nWbegbgF6q2qDY7zmUrx3cg5UDCcDd+YOVHW7qnYozP6NKSgrLiYSHQLeBO7L/w/5Zx4iss/3v0ki\n8o6IzBSRb0VkiIj0EJGVIvKFiJznt5uGIvKJiGwQkea+98eLyPO+138uIrf67XeFiMzCexuV/PF0\n8e1/rYgM8W0bCFwLjPQ9P+MffPeFexzoJN4HnHUQkXK+Bzx97LujcSvfa3uJyGzfLGiJ73VLReRT\n37Fb+XY7BDjft79nReQcEfnSt4/Svjslf+Hbd5LfvqeLyALxPkDqWb98pPl+ri9E5N6C/V9nYkVh\n/tIyJpy8CnxxhF/O+Wce/uOLgYvwPsPie+AtVa0j3qd03o23WAlwjqpeKSIXABm+/03Fe0v2OiJS\nCnhPRBb79nspUEtVf/A/sIhUxvsL/TJgN7BYRFqr6uMikoz3NjSfHemHU9VDIjIAuFxVc58i+jSw\nTB0Qv1kAAALASURBVFVvFu8TAj8WkaV+MdRW1d3ifVDejaq6V0ROBT7Ee2+o/r44L/XtL9EvP3cB\nh1X1YhG50Bdr7v3HLgE8eGeMG0TkJby3Yq/s95Cpk470c5jYZTMXE5F8jw8YDfy7EG/7RFV3qupB\nvA87yr030pd4b7sO3l+2k33H2AR8h7cgNQZ6ishq4COgInCB7z0r8xcWnyuBDN9ddw8D44Ab/P79\niLf/z/fv/q9pDDzoiyEDKIX3HlAKLPG7/1Uc8IyIfA4sASqLyGnHOd61wFjfz70B+AGo7tv3MlXd\nq6oH8N7Y8Wy89yQ7T0SGi0gTvA/eMiaPzVxMJBuG92aDo/y2ZeP7o0lE4vDenDHXAb/vc/zGORz7\nv4Xcv+77+h4mlce3fPTnMd7n/wtd+PtMqij9nbaqujFfDHXzxdAN780IL1PVw+J91EDpAuz7aMXH\nP2+HgQTfDOkSoAlwO9ARbx/JGMBmLiaCqeouvLOMW3B+UW8GLvd93wooUcjdCtBBvM4HzsN7a/ZF\nwJ25DXMRqS4iZY+zr0+AeiJyim+pqjPwTiFi2YP3Dr25FuE3UxORS/1i9lce78PRDvuW387xbd+b\nb3/+VuAtSviWw87G+3MfqeCI72SKeFWdDgzAu/RnTB4rLiYS+f/F/1+8f6XnegvvL/Q1wFXAvqO8\nL//+1O/7H/E+mmA+cJtvGW0E3iWhz3ynCr+Gd7Zz1Kf1qepPeJ/ymAGsAVap6pxC/HwZwP+3d8cm\nCARBFED/CjZjbDtiLWJkAzZgYhdG14AIVmEFBmtwa3JoNuF7+cBs9Nkd2Nl8B/pJjknWY4D+SHL4\n0X8yP79tW2v3JLskz9HPK8k0hvCnRd05yWrUXJPse+/vP+frmVfc3sYT3SWLbZbgy30Ayrm5AFBO\nuABQTrgAUE64AFBOuABQTrgAUE64AFBOuABQ7gPc2rZUru5YjQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#Array of iterations\n",
"Niter = np.arange( 1, 100, 1 )\n",
"\n",
"plt.figure( figsize=(6,6) )\n",
"plt.semilogy( Niter, 2.0**-Niter, color=\"green\", lw = 2 )\n",
"plt.grid(True)\n",
"plt.xlabel(\"Number of Iterations\")\n",
"plt.ylabel(\"Absolute Error $|p_n-p|$\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Find the root of the function $f$\n",
"\n",
"$f(x) = x^3 - 2$\n",
"\n",
"for $20$ iterations, show the result and the relative error in each iteration."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Defining Bisection function\n",
"def Bisection( f, a, b, Nmax, printer=False ):\n",
" #verifying the STEP1, a and b with different signs\n",
" if f(a)*f(b)>0:\n",
" print \"Error, f(a) and f(b) should have opposite signs\"\n",
" return False\n",
" #Assigning the current extreme values, STEP2\n",
" ai = a\n",
" bi = b\n",
" #Iterations\n",
" n = 1\n",
" while n<=Nmax:\n",
" #Bisection, STEP3\n",
" pi = (ai+bi)/2.0\n",
" #Evaluating function in pi, STEP4 and STEP5\n",
" if printer:\n",
" print \"Value for %d iterations:\"%(n),pi\n",
" #Condition A\n",
" if f(pi)*f(ai)>0:\n",
" ai = pi\n",
" #Condition B\n",
" elif f(pi)*f(ai)<0:\n",
" bi = pi\n",
" #Condition C: repeat the cycle\n",
" n+=1\n",
" #Final result\n",
" return pi"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Value for 1 iterations: 1.0\n",
"Value for 2 iterations: 1.5\n",
"Value for 3 iterations: 1.25\n",
"Value for 4 iterations: 1.375\n",
"Value for 5 iterations: 1.3125\n",
"Value for 6 iterations: 1.28125\n",
"Value for 7 iterations: 1.265625\n",
"Value for 8 iterations: 1.2578125\n",
"Value for 9 iterations: 1.26171875\n",
"Value for 10 iterations: 1.259765625\n",
"Value for 11 iterations: 1.2607421875\n",
"Value for 12 iterations: 1.26025390625\n",
"Value for 13 iterations: 1.26000976562\n",
"Value for 14 iterations: 1.25988769531\n",
"Value for 15 iterations: 1.25994873047\n",
"Value for 16 iterations: 1.25991821289\n",
"Value for 17 iterations: 1.25993347168\n",
"Value for 18 iterations: 1.25992584229\n",
"Value for 19 iterations: 1.25992202759\n",
"Value for 20 iterations: 1.25992012024\n",
"Real value: 1.25992104989\n",
"Absolute error 9.29655615378e-07\n"
]
}
],
"source": [
"#Defining function\n",
"def function(x):\n",
" f = x**3.0 - 2.0\n",
" return f\n",
"\n",
"#Finding the root of the function. The real root is 2**(1/3), so a & b should enclose this value\n",
"a = 0.0\n",
"b = 2.0\n",
"Nmax = 20\n",
"\n",
"result = Bisection(function, a, b, Nmax, True)\n",
"print \"Real value:\", 2**(1/3.0)\n",
"print \"Absolute error\", abs((2**(1/3.0)-result))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the error analysis, we can predict the produced error at $20$ iterations by computing:\n",
"\n",
"$$ \\left( \\frac{1}{2^{20}}\\right) \\approx 9.53674316\\times 10^{-7} $$\n",
"\n",
"This value is very close to the obtained relative error.\n",
"\n",
"If we were interested in a double precision, i.e. $\\epsilon \\sim 10^{-17}$, the number of required iterations would be:\n",
"\n",
"$$ 10^{-17} = \\left( \\frac{1}{2^{N}}\\right) \\longrightarrow N = \\frac{17}{\\log_{10}(2)} \\approx 56 $$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### ACTIVITY "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"In an IPython notebook, copy the last function and find the first solution to the equation\n",
" \n",
"$ 7 = \\sqrt{x^2+1}+e^x\\sin x $\n",
" \n",
"CLUE: this solution is within the interval $[0,2]$.\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In orbital mechanics, when solving the central-force problem it becomes necessary to solve the Kepler's equation. This is a transcendental equation that relates the orbital parameters of the trajectory.\n",
"\n",
"*Kepler equation:* $M = E - \\epsilon \\sin E$\n",
"\n",
"where $M$ is the mean anomaly, $E$ the eccentric anomaly and $\\epsilon$ the eccentricity. The mean anomaly can be computed with the expression\n",
"\n",
"$$M = n\\ t = \\sqrt{ \\frac{GM}{a^3} } t$$\n",
"\n",
"where $n$ is the mean motion, $G$ the gravitational constant, $M$ the mass of the central body and $a$ the semi-major axis. $t$ is the time where the position in the trajectory will be computed.\n",
"\n",
"The coordinates $x$ and $y$ as time functions can be recovered by means of the next expressions\n",
"\n",
"$$x(t) = a(\\cos E - \\epsilon)$$\n",
"\n",
"$$y(t) = b\\sin E$$\n",
"\n",
"where $b = a \\sqrt{1-\\epsilon^2}$ is the semi-minor axis of the orbit and the implicit time-dependence of the eccentric anomaly $E$ is computed through the Kepler's equation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Problem:**\n",
"\n",
"For a stallite orbiting the earth in a equatorial trajectory with eccentricity $\\epsilon = 0.5$ at a geostationary distance for the semi-major axis, tabulate the positions $x$ and $y$ within the orbital plane in intervals of $15$ min during $5$ hours.\n",
"\n",
"**Parameters:**\n",
"\n",
"- $\\epsilon = 0.5$\n",
"\n",
"- $a = 35900$ km\n",
"\n",
"- $G = 6.67384 \\times 10^{-11}$ m$^3$ kg$^{-1}$ s$^{-2}$\n",
"\n",
"- $M_{\\oplus} = 5.972\\times 10^{24}$ kg"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#====================================================================\n",
"#Parameters\n",
"#====================================================================\n",
"#Eccentricity\n",
"eps = 0.5\n",
"#Semi-major axis [m]\n",
"a = 35900e3\n",
"#Gravitational constant [m3kg-1s-2]\n",
"GC = 6.67384e-11\n",
"#Earth mass [kg]\n",
"Me = 5.972e24\n",
"\n",
"#Semi-minor axis [m]\n",
"b = a*(1-eps**2.0)**0.5\n",
"#Mean motion\n",
"n = ( GC*Me/a**3.0 )**0.5\n",
"\n",
"#Hour to Second\n",
"HR2SC = 3600.\n",
"#Initial time [hr]\n",
"t0 = 0*HR2SC\n",
"#Final time [hr]\n",
"tf = 5*HR2SC\n",
"#Time step [hr]\n",
"tstep = 0.25*HR2SC\n",
"#Number of maxim iterations\n",
"Niter = 56\n",
"#Root interval\n",
"a0 = -10\n",
"b0 = 10\n",
"\n",
"#====================================================================\n",
"#Kepler Function\n",
"#====================================================================\n",
"def kepler( E ):\n",
" func = E - eps*np.sin(E) - n*t\n",
" return func\n",
"\n",
"#====================================================================\n",
"#Position function\n",
"#====================================================================\n",
"def r(E):\n",
" x = a*(np.cos(E)-eps)\n",
" y = b*np.sin(E)\n",
" return [x/1.e3, y/1.e3]"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"In 0.000000 hours, the satellite is located at (17950.000000,0.000000) km\n",
"In 0.250000 hours, the satellite is located at (17454.741542,5146.426647) km\n",
"In 0.500000 hours, the satellite is located at (16033.097675,10023.437750) km\n",
"In 0.750000 hours, the satellite is located at (13848.847528,14430.262364) km\n",
"In 1.000000 hours, the satellite is located at (11104.379909,18261.701894) km\n",
"In 1.250000 hours, the satellite is located at (7989.437466,21493.410338) km\n",
"In 1.500000 hours, the satellite is located at (4657.168469,24151.489610) km\n",
"In 1.750000 hours, the satellite is located at (1221.066166,26286.121177) km\n",
"In 2.000000 hours, the satellite is located at (-2238.747501,27954.872718) km\n",
"In 2.250000 hours, the satellite is located at (-5667.264143,29214.008624) km\n",
"In 2.500000 hours, the satellite is located at (-9027.373694,30114.739827) km\n",
"In 2.750000 hours, the satellite is located at (-12294.392344,30702.085866) km\n",
"In 3.000000 hours, the satellite is located at (-15452.164136,31014.965936) km\n",
"In 3.250000 hours, the satellite is located at (-18490.361033,31086.789919) km\n",
"In 3.500000 hours, the satellite is located at (-21402.633716,30946.195086) km\n",
"In 3.750000 hours, the satellite is located at (-24185.347784,30617.769816) km\n",
"In 4.000000 hours, the satellite is located at (-26836.717825,30122.702148) km\n",
"In 4.250000 hours, the satellite is located at (-29356.211228,29479.336137) km\n",
"In 4.500000 hours, the satellite is located at (-31744.135350,28703.638662) km\n",
"In 4.750000 hours, the satellite is located at (-34001.350021,27809.586870) km\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEQCAYAAAAZPssSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHX+P/DXyKBYXlBSSIZEBwQRArygXSQMYfCGpoVr\npqB0w/VSmV8rtaxEdHfd/aErbZnmLcWyxAuJmInQFpgImWHGEigMiAaEqFwE3r8/jpxERh1kZs7M\n4f18PHgw85lzzrzfjM278/l8zucoiIjAGGOMSaCD1AEwxhhrv7gIMcYYkwwXIcYYY5LhIsQYY0wy\nXIQYY4xJhosQY4wxyZisCNXU1GD48OHw8fGBh4cH3nzzTQDA8uXLoVKp4OvrC19fXxw8eFDcJyYm\nBq6urnB3d0dycrLYnpmZCS8vL7i6umLBggVie21tLaZOnQpXV1eMGDEC586dM1V6jDHG7oHJipCN\njQ2OHj2K7OxsnDp1CkePHsW3334LhUKB1157DVlZWcjKysKYMWMAADk5Odi1axdycnKQlJSEOXPm\noOmSpqioKGzcuBG5ubnIzc1FUlISAGDjxo2ws7NDbm4uXn31VSxevNhU6THGGLsHJu2Ou++++wAA\ndXV1aGhoQI8ePQAAuq6X3bt3L6ZNmwZra2s4OzvDxcUFGRkZKCkpQVVVFfz8/AAAM2fOREJCAgBg\n3759CA8PBwBMmTIFR44cMUVajDHG7pFJi1BjYyN8fHxgb2+PUaNGYdCgQQCAdevWwdvbG5GRkfjj\njz8AAMXFxVCpVOK+KpUKWq22RbujoyO0Wi0AQKvVwsnJCQCgVCrRvXt3lJeXmyo9xhhjraQ05Zt1\n6NAB2dnZqKyshEajQUpKCqKiovD2228DAJYtW4aFCxdi48aNRo1DoVAY9fiMMSZHxljlTZLZcd27\nd8e4ceNw4sQJ9O7dGwqFAgqFAs8//zyOHz8OQDjDKSwsFPcpKiqCSqWCo6MjioqKWrQ37XP+/HkA\nQH19PSorK9GzZ0+dMRCRLH/eeecdyWPg/Dg/zk9+P8ZisiL0+++/i11t1dXVOHz4MHx9fXHhwgVx\nmz179sDLywsAEBoaivj4eNTV1SE/Px+5ubnw8/ODg4MDunXrhoyMDBARtm3bhokTJ4r7bNmyBQCw\ne/duBAYGmio9s1FQUCB1CEbF+Vk2zo/dymTdcSUlJQgPD0djYyMaGxsxY8YMBAYGYubMmcjOzoZC\noUC/fv3w4YcfAgA8PDwQFhYGDw8PKJVKxMXFid1ocXFxiIiIQHV1NcaOHYuQkBAAQGRkJGbMmAFX\nV1fY2dkhPj7eVOkxxhi7Bwoy5nmWmVIoFEY9vZRSSkoKAgICpA7DaDg/y8b5WS5jfW9yEWKMMXZX\nxvre5GV7ZCYlJUXqEIyK87NsnB+7FRchxhhjkuHuOMYYY3fF3XGMMcZkh4uQzMi9T5rzs2ycH7sV\nFyHGGGOS4TEhxhhjd8VjQowxxmSHi5DMyL1PmvOzbJwfuxUXIcYYY5LhMSHGGGN3xWNCjDHGZIeL\nkMzIvU+a87NsnB+7FRchxhhjkuExIcYYY3fFY0KMMcZkh4uQzMi9T/pO+aUmJmKpRoPlAQFYqtEg\nNTGx1cc3xDHaoj1/fnIg9/yMQSl1AIwBwpd/8tq1UNbWor5TJwTPnw//cePuvBMRUFMDXL6M1L17\ncej99xFdVCS+vCQnB1iyBP4aDdCpU/MfpRJQKFrEcGjBAkTn5f15jBuP7xoLY+ye8JgQk1zq/v04\nNG8eos+dE9uW9OwJzahR8O/VC6isBC5f1v1z/ToAYCmAFTqOvQzA+7reVKEQilHHjmJhWlpWhhXV\n1S2P4eWF99euBZycAJVK2P5O+dxLQWXMzBnre5PPhJjB3PHL98oVID8fyMsDfvut2e/k3FxE33Ks\n6PJyLPviC/jf7U07dQK6dYOyqko4K7qFVefOQO/eQG1t85+GBmH7m/a53X8MVj/9BIwa9WdD795C\nQdLxk5qbi0MrVvDZFGN64iIkMykpKQgICDD5+6YmJuLQ/PmI/u03sW3Jd98BKhX8KyqA0tLb7nvb\nL//+/YHXXwe6dRN/Un79FQGBgcLzrl3Fs5J6jQZITm5xjAZ/fyApqeXBGxqAurpmhal++nTg229b\nburgAKjVQGEhoNUCFy8KP5mZLbZNBloW1Lw8LFu9Gv5jxgAd7jwMK9XnZyqcH7sVFyF2bxobgV9/\nBY4fBzIykLx9O6IvX262SfSVK1j2yy/C2UzHjkC/fsKXef/+zX7XL1gAHDnS4i0aXF2BqKjmjUTC\nfrcInj8fS/Lymp2BvKVWI2TePN3xW1kBnTsLP03HeOMNLLllTOgttRohsbFA01lMQwNw4YJQkHT8\nKLOzheJ269ulpQE9ewJDhgDDhgFDhwq/H3qoxdgUY+0KmUh1dTX5+fmRt7c3DRw4kN544w0iIior\nK6PRo0eTq6srBQUFUUVFhbjPypUrycXFhdzc3OjQoUNi+4kTJ8jT05NcXFxo/vz5YntNTQ2FhYWR\ni4sLDR8+nAoKCnTGYsK0Lc6xAwdoSXAwvfPEE7QkOJiOHTggvFBaSrRvH9HSpURBQUS2tkRCSSAC\n6J2bHjdr9/EhOn+eqKHhju/5llrdbL831eo/37sVsS/VaOidJ56gpRpNq/c3xDGWBAfr/Dss7dRJ\nZzv16kU0ZgzR228Lf9+Skmax6PwsGJOAsb43TfptfPXqVSIiun79Og0fPpzS0tJo0aJFtHr1aiIi\nWrVqFS1evJiIiH7++Wfy9vamuro6ys/PJ7VaTY2NjURENGzYMMrIyCAiojFjxtDBgweJiGj9+vUU\nFRVFRETx8fE0depUnXFwEdJNVzF4q0sXOtarl+4v0D59iJ56imjVKloydKjuL1+NRu/3bmsBMQd3\nLKhaLdHevUIh12iI7Ox0/11VKjo2YgS91bNn88/iHgozY4YiiyLU5OrVqzR06FA6ffo0ubm50YUL\nF4iIqKSkhNzc3IhIOAtatWqVuI9Go6Hvv/+eiouLyd3dXWzfuXMnvfTSS+I26enpRCQUugceeEDn\n+8u5CB09evTediwooCUDB+ouJABRly5EAQFEixcTffEFUVFRs90NdTZzN/ecnwnpXVAbG4l++43o\ns8+IFi0iCgigo507EwG05DZnlksff1zYz0JZwufXFnLOz1jfmyYdE2psbMTgwYORl5eHqKgoDBo0\nCKWlpbC3twcA2Nvbo/TGAHZxcTFGjBgh7qtSqaDVamFtbQ2VSiW2Ozo6QqvVAgC0Wi2cnJwAAEql\nEt27d0d5eTl69uxpqhQtx/XrwHffAV99BSQmAj//fPsJAkOHAunpwjjKbTTN/Fq2bh2samrQYGOD\nkHnz2uWMMP9x4/TLW6EQxsn69QOeeUZo++YboE8fKJ96Cvjllxa7WH37LTBgABAaCkycCDz6qHDN\nE2MWyqT/ejt06IDs7GxUVlZCo9Hg6NGjzV5XKBRQmGiQNiIiAs7OzgAAW1tb+Pj4iLNamq56tsTn\nAQEBSElJwY/p6bh09CiUtbXIu3oVQ6dMwYLZs4GDB5GyeTPwww8IuHpV2B8AOndGfdeuwMWLwnMA\nATd+/6ZQICUt7e7vf+PLt+m5vxHzM9bfT+rnAU8+iZSUFOR16YImKTd+BwBosLZGyv/+B/zznwj4\n5z8BOzukDBkCPPYYAl59Feja1azyaZGf3D8/GeXX9LigoABGZZTzKz2899579Pe//53c3Nyo5MZg\nbHFxsdgdFxMTQzExMeL2TV1tJSUlzbrjduzYQS+//LK4zffff09E7bc7jug2YzudOtGxW7t33N2J\nFi4kOnKEqLbWZF1q7O5u+1ns20eUlkb0+utErq7NP8+OHYVJDh98IIw/MWZAxvreNNm38aVLl8SZ\nb9euXaORI0fS119/TYsWLRLHfmJiYlpMTKitraXffvuN+vfvL05M8PPzo/T0dGpsbGwxMaGpIO3c\nubNdTkw4evQoLXnsMd3jCQoFUUgI0dq1RHl5Ovc39wkCcu5zJ2qen16fxZkzRKtXEz36KJFC0fwz\nHzqU6P33iX78kY7t328WM+3a0+cnNxZfhE6dOkW+vr7k7e1NXl5e9Le//Y2IhCnagYGBOqdoR0dH\nk1qtJjc3N0pKShLbm6Zoq9VqmjdvntheU1NDzzzzjDhFOz8/X2cssixCZWVE69fTUXf320+XHjlS\n6ijbTM7/kRO1Mb/SUqKNG4kmTiS6McGBADoG0FtKpVnMtOPPz3IZ63uT146zZNevC6sBbNkC7N8v\nXiS51MoKKxoaWmy+TKPB+7pWD2Dyc+2acAHw3r1YunUrVtxYY+9mywID8f7XX0sQHLNEvHZcO9Zi\nTbYJE+Cflwfs2CEsHwMIM62CgoDwcATb2GDJ4sX6rx7A5Oe++4AJE4AJE6D83/+AY8dabGKVkgLM\nmQP89a/AoEGmj5ExcBEyezpvL3BjjTR/AHB3B8LDgeeeA1QqpDStXWVjI8vp0mJ+MmWM/Opvs+p3\nQ0MD8MEHws+oUUIxmjjRqFO++fNjt+IiZOaSo6ObFSBAWCBzmZMT/HfvFtYf0zGtXe9rVZjs3XZd\nvVdfBU6fBrZtA44eFX4cHYGXXwZeeAG4cf0eY8bEY0LmiEhYEfof/8Dyr7/Gch2bLH/iCSznuzgy\nPaUmJuLwTWfGQTefGVdWCuOK69cLi9ICgLW1cAHtX/8KPPIIL7LKjPa9yUXInNTVATt3Av/4h/B/\nqACWduiAFY2NLTblSQbM4BobhckM69cLE12a/t35+grFaNo0YayJtUvG+t68881NmGlUVACrVgnL\nt0RECAXowQeBVasQvGMHlqjVzTZ/S61G0G0mGaTI/OyI8zOiDh2EyS0JCcINB994A3jgASArC3j+\neeGusq+/jtSPP8ZSjQbLAwKwVKNBamKi3m/Bnx+7FY8JmVizmW5ECLa1hf+RI8CNJXTg6SncyG3a\nNKBjR2HyQZcuspxkwMxY375ATAzwzjvAZ58JZ0fHjyN1zRocQvMb9/GdY1lbcHecCemc6QZAA8B/\n9Gih+AQHc/87M08//IClkyZhRXFxi5e4e1j+uDtOBpLXrtU50+3wo48Chw8DGg0XIGa+hg2D0tVV\n50tWJ08Ct/zbZkwfXIRMSFlbq7PdytraYO8h9z5pzk9at73m6NIlYOBA4JVXgLKy2+5v7vm1ldzz\nMwYuQiZ02/+AbWxMHAlj9yZ4/vyWE2X69kVQYCBQXw/ExgJqNfC3vwE1NRJFySwJjwmZkK4xobfU\naoTExvKgLrMYt73m6McfgUWLhK5lAHjoIWDlSmGSTQf+/11Lx9cJGZCU1wnd8aJBxuTg0CGhGP30\nk/B88GDh2rdRo6SNi7UJFyEDMtuLVQ1A7mtXcX4WoqEB2LoVWLoUaJpNN24cUp5+GgEREZKGZkyy\n+fx04NlxjDHLYWUFzJoF5OYCK1YAXboAiYnA7NnAiy8CJSVSR8jMBJ8JMcaMr7QUePdd4KOPhLOk\n++8Xrot7/XWkHjvW/FYl8+dzF7UZ4u44A+IixJhEfvlFWA5o714AQKqtLQ5ZWyP60iVxkyVqNTQ8\nWcfscHcc04vcr1Pg/CxbyoULwtp0x44Bw4Yh+Y8/mhUgAIjOy8PhdeskirBt5P75GQMXIcaY6fn7\nA+npUA4cqPNlK77GqN3gIiQzcp2Z04Tzs2zN8uvQAfVOTjq3a9Bx+xJLIPfPzxi4CDHGJKNzBQYA\nQRkZwIYNwg0emaxxEZIZufdJc36W7db8/MeNgyY2Fss0Gix/4gksCwxEyBNPwL+uTpjKPXky8Pvv\n0gR7D+T++RkD30+IMSYp/3HjWs6E+/RTYM4cYRJDRoZw+/GgIGkCZEZlsjOhwsJCjBo1CoMGDYKn\npyfWrl0LAFi+fDlUKhV8fX3h6+uLgwcPivvExMTA1dUV7u7uSE5OFtszMzPh5eUFV1dXLFiwQGyv\nra3F1KlT4erqihEjRuDcuXOmSs9syL1PmvOzbHrnN326sBbd448LF7YGBwMLFwK3WYneXMj98zMK\nMpGSkhLKysoiIqKqqioaMGAA5eTk0PLly2nNmjUttv/555/J29ub6urqKD8/n9RqNTU2NhIR0bBh\nwygjI4OIiMaMGUMHDx4kIqL169dTVFQUERHFx8fT1KlTdcZiwrQZY21RX0+0YgWRlRURQPTww0Sn\nT0sdVbtkrO9Nk50JOTg4wMfHBwDQpUsXDBw4EFqttqkQtth+7969mDZtGqytreHs7AwXFxdkZGSg\npKQEVVVV8PPzAwDMnDkTCQkJAIB9+/YhPDwcADBlyhQcOXLEFKmZFbn3SXN+lq3V+VlZAUuWAN99\nJ9wi4tQpYOhQYN06s5y0IPfPzxgkmZhQUFCArKwsjBgxAgCwbt06eHt7IzIyEn/88QcAoLi4GCqV\nStxHpVJBq9W2aHd0dBSLmVarhdONKZ9KpRLdu3dHeXm5qdJijBmLnx+QnQ1ERgr3KZo/Hxg3TlgO\niFk0k09MuHLlCp5++mnExsaiS5cuiIqKwttvvw0AWLZsGRYuXIiNGzcaPY6IiAg4OzsDAGxtbeHj\n4yP25zb934wlPg8ICDCreDg/zs+g+X38MVIeegj4+98RcPAg4OWFlFdfBR55RB75mdHzpscFBQUw\nKqN08t1GXV0dBQcH07/+9S+dr+fn55OnpycREcXExFBMTIz4mkajofT0dCopKSF3d3exfceOHfTy\nyy+L23z//fdERHT9+nV64IEHdL6PidNmjBlaURFRYKAwTgQQzZlDdPWq1FHJmrG+N03WHUdEiIyM\nhIeHB1555RWxveSmJd337NkDLy8vAEBoaCji4+NRV1eH/Px85Obmws/PDw4ODujWrRsyMjJARNi2\nbRsmTpwo7rNlyxYAwO7duxEYGGiq9MzGzf8XI0ecn2UzWH6OjkBysnCzPGtrIC4OGDoUqevWYalG\ng+UBAViq0SA1MdEw76cnuX9+xmCy7rj//ve/2L59Ox5++GH4+voCAFauXImdO3ciOzsbCoUC/fr1\nw4cffggA8PDwQFhYGDw8PKBUKhEXFweFQgEAiIuLQ0REBKqrqzF27FiEhIQAACIjIzFjxgy4urrC\nzs4O8fHxpkqPMWZqHToI07YDA4Fnn0XqmTM4NH8+om/aZEleHgDwitxm7K63ctBnYL9Dhw6wtbU1\nWFDGxrdyYExmrl3DUnd3rCgsbPHSMo0G7yclSRCUvBjre/OuZ0IPPvgg+vTpc8dt6uvrUajjw2eM\nMZO47z4o+/cHdHwP8Yrc5u2uY0IDBw5Efn7+HX/s7OxMESvTg9z7pDk/y2bM/Oo7ddLZ3mC0d2xJ\n7p+fMdy1CKWnp9/1IPpswxhjxnTbFbl//hk4c0aaoNhd8e29GWOykZqYiMPr1sGqpgYNVlYIKi6G\n/y+/AD17Al99BQwfLnWIFstY35t6F6EffvgBK1euREFBAerr68WgTp06ZfCgjI2LEGPtxLVrQFgY\nkJgI3Hcf8OWXgEYjdVQWyVjfm3pfJzR9+nTMmjULX3zxBfbv34/9+/dj3759Bg+ItY3c+6Q5P8tm\n8vzuuw/YswcIDxcK0vjxwI4dRns7uX9+xqD3dUK9evVCaGioMWNhjDHDs7YGPvkE6NVLuLh1+nTg\n0iXgptvAMOno3R2XnJyMXbt2YfTo0ejYsaOws0KByZMnGzVAY+DuOMbaqb//Hfi//xMev/kmEB0N\n3LgInt2Z5GNC06dPx9mzZzFo0CB06PBnL94nn3xi8KCMjYsQY+3Yli3CatwNDcLv//wHUPJNpu9G\n8iLk5uaGX375RVw6x5LJuQilpKSIq+HKEedn2cwmvwMHhAkL1dXAxInAzp1A585tPqzZ5GcEkk9M\nePTRR5GTk2PwABhjzOTGjwe+/hqwtQX27gVCQoAb9zJjpqX3mZC7uzvy8vLQr18/dLpxZTJP0WaM\nWbTTp4Up28XFwMMPA0lJwIMPSh2VWZK8O07XjY0UCgX69u1r6JiMjosQY0x07hwQHAz8+ivQrx9w\n6BDg6ip1VGZH8u64srIyODs7N/v56aefDB4Qaxu5X6fA+Vk2s8yvb1/g22+BYcOA/Hzg8ceBkyfv\n6VBmmZ+Z07sIvfDCC82Kzs6dO/Hee+8ZJSjGGDOpXr2Ab74BgoKAixeR+vjjWDpsmGQ3x2tP9O6O\n++233/D0009jx44dSEtLw9atW3HgwAF0797d2DEaHHfHMcZ0qqtD6ujROJSW1vzmeGo1NLGx7frm\neJKPCQHA2bNnMWnSJPTt2xdffvkl7rvvPoMHZApchBhjt7NUo8GK5OQW7e395niSjQl5eXmJP08/\n/TTKy8uRn5+P4cOH4+GHHzZ4QKxt5N4nzflZNkvIT1lbq7Ndn5vjWUJ+5uaulwnv37/fFHEwxphZ\nMIeb47UnfD8hxhi7SWpiIg4tWIDovDyx7S0AIf36wf/UKaBLF+mCk5BkY0KDBw/GybtMV9RnG3PC\nRYgxdictbo535gz8S0qA0FDhnkRWVlKHaHKSFaHOnTvDxcXljgeprKzE+fPnDRqYMcm5CMl57SqA\n87N0Fpvfr78CI0YAFRXAokXA3/6mczOLzU8PxvrevOuY0Bk97s2u5BVoGWNyNmAAsHu3sMTP3/8O\nuLkJK3CztiMTOX/+PAUEBJCHhwcNGjSIYmNjiYiorKyMRo8eTa6urhQUFEQVFRXiPitXriQXFxdy\nc3OjQ4cOie0nTpwgT09PcnFxofnz54vtNTU1FBYWRi4uLjR8+HAqKCjQGYsJ02aMycmGDUQAkVJJ\ndPSo1NGYlLG+N032bVxSUkJZWVlERFRVVUUDBgygnJwcWrRoEa1evZqIiFatWkWLFy8mIqKff/6Z\nvL29qa6ujvLz80mtVlNjYyMREQ0bNowyMjKIiGjMmDF08OBBIiJav349RUVFERFRfHw8TZ06VWcs\nXIQYY/ds4UKhEPXoQfTrr1JHYzLG+t7Ue9metnJwcICPjw8AoEuXLhg4cCC0Wi327duH8PBwAEB4\neDgSEhIAAHv37sW0adNgbW0NZ2dnuLi4ICMjAyUlJaiqqoKfnx8AYObMmeI+Nx9rypQpOHLkiKnS\nMxtyv06B87Nssshv9WrhVhAVFX/+vkEW+ZmY3oM5jY2N+PTTT5Gfn4+3334b58+fx4ULF8Ri0BoF\nBQXIysrC8OHDUVpaCnt7ewCAvb09SktLAQDFxcUYMWKEuI9KpYJWq4W1tTVUKpXY7ujoCK1WCwDQ\narVwcnISElMq0b17d5SXl6Nnz54tYoiIiICzszMAwNbWFj4+PuKAYtM/JH7Oz/k5P2/xPC0NiIpC\nwPnzwKlTSAkMBP72NwSMHm0e8RnoedNjXXdQMCh9T5leeuklioqKIjc3NyISxnKGDBnS6lOvqqoq\nGjx4MO3Zs4eIiGxtbZu93qNHDyIimjt3Lm3fvl1sj4yMpN27d9OJEydo9OjRYntqaiqNHz+eiIg8\nPT1Jq9WKr6nVaiorK2sRQyvSZowx3c6dI7K3F7rmXnyR6MZwgVwZ63tT7+64jIwMxMXFofONW+D2\n7NkT169fb1XBu379OqZMmYIZM2Zg0qRJAISznwsXLgAASkpK0Lt3bwDCGU5hYaG4b1FREVQqFRwd\nHVFUVNSivWmfpqni9fX1qKys1HkWxBhjbfbQQ8JdWW1sgI8+Av7f/5M6IoukdxHq2LEjGhr+XLji\n0qVL6NBB/yElIkJkZCQ8PDzwyiuviO2hoaHYsmULAGDLli1icQoNDUV8fDzq6uqQn5+P3Nxc+Pn5\nwcHBAd26dUNGRgaICNu2bcPEiRNbHGv37t0IDAzUOz65uPlUWo44P8smu/yGDwdufOdg4UKkrFwp\nbTyWSN9Tpm3bttGECROoT58+9Oabb5Krqyvt2rVL71OutLQ0UigU5O3tTT4+PuTj40MHDx6ksrIy\nCgwM1DlFOzo6mtRqNbm5uVFSUpLY3jRFW61W07x588T2mpoaeuaZZ8Qp2vn5+TpjaUXaFueozKeN\ncn6WTbb5vfceEUBHbWyIsrOljsYojPW92aq1486cOSPOOAsMDMTAgQONVBqNS84rJjDGJEAEPPcc\nsGMH4OQEHD8OODhIHZVBSX4/ocWLF2P16tV3bbMEXIQYYwZXUwM8+STw/feAnx+QkgLcGEOXA8nu\nJ9QkWcdNnr766iuDBsPaTnZ97rfg/CybrPOzsUHKokWAs7NwJjRrlnCGxO7orkXogw8+gJeXF86e\nPdvsBnfOzs58UzvGGLtZjx7A/v1A167Arl3Au+9KHZHZu2t3XGVlJSoqKvDGG29g9erV4ulY165d\nYWdnZ5IgDY274xhjRnXwIDB+PFIbG5Hs5QVlz56o79QJwfPnw3/cOKmjuyeSjwkBQEVFBXJzc1Fz\n021u/f39DR6UsXERYowZW+pLL+HQRx8h+qa2JWo1NLGxFlmIJB8T2rBhA/z9/aHRaPDOO+9Ao9Fg\n+fLlBg+ItY2s+9zB+Vm69pRfckFBswIEANF5eTi8bp1JYzJ3eheh2NhYHD9+HH379sXRo0eRlZWF\n7t27GzM2xhizWMraWp3tVjf1JLFWFCEbGxtxyZ6amhq4u7vj7NmzRguM3ZumRQjlivOzbO0pv/pO\nnXRu02BjY6JoLIPeRcjJyQkVFRWYNGkSgoKCEBoaKq5CzRhjrLng+fOxRK1u1vZWp04I+utfJYrI\nPLVqYkKTlJQUXL58GSEhIejYsaMx4jIqOU9MSJHxPe4Bzs/Stbf8UhMTcXjdOlhdvYqGH35AUG0t\n/NetA+bOlS7Ie2Ss70297ydUU1ODL774AgUFBaivrwcAZGdn4+233zZ4UIwxJgf+48b9ORNuzx5g\n8mTgrbeAp54CHB2lDc5M6H0mpNFoYGtriyFDhsDKykpsX7hwodGCMxY5nwkxxswUETBpErBvn1CM\nvvhC6ohaRfLrhDw9PXH69GmDByAFLkKMMUkUFgIDBwJXrwr3IgoNlToivUl+ndCjjz6KU6dOGTwA\nZljt6ToMOeL8LNtd83NyAlasEB7PnQtcuWL0mMyd3kUoLS0NQ4YMwYABA8T143jtOMYYa6V584Ah\nQ4SzIh71lAaZAAAftElEQVRT1787rqCgQGe7JU7T5u44xpikTp4Ehg0THv/wAzB4sLTx6EHyMSE5\n4SLEGJPca68B//qXcFaUng4o9Z6sLAnJxoQee+wxAECXLl3QtWvXZj/dunUzeECsbdp9n7uF4/ws\nW6vye+89YYwoMxNYv95oMZm7uxah//73vwCAK1euoKqqqtnP5cuXjR4gY4zJUpcufxafpUuFMaJ2\niLvjGGNMSlOmAF9+CUycCCQkSB3NbUk2JrRmzRqdQSgUCgDAa6+9ZvCgjI2LEGPMbGi1wrVDVVXC\nqgqTJkkdkU6SjQlVVVXhypUryMzMxAcffIDi4mJotVr85z//wcmTJw0eEGsb7nO3bJyfZbun/Bwd\ngZUrhcdz5wLtbJjjrtMxmm5cN3LkSJw8eRJdu3YFALz77rsYO3asUYNjjLF2ISoK2LYNOH4cWLYM\niI2VOiKT0fti1YsXL8La2lp8bm1tjYsXL7bqzWbPng17e3t4eXmJbcuXL4dKpYKvry98fX1x8OBB\n8bWYmBi4urrC3d0dycnJYntmZia8vLzg6uqKBQsWiO21tbWYOnUqXF1dMWLECJw7d65V8cmBnFco\nBjg/S8f53YaVFfDRR8LvdeuEa4faCb2L0MyZM+Hn54fly5fjnXfewfDhwxEeHt6qN5s1axaSkpKa\ntSkUCrz22mvIyspCVlYWxowZAwDIycnBrl27kJOTg6SkJMyZM0fsj4yKisLGjRuRm5uL3Nxc8Zgb\nN26EnZ0dcnNz8eqrr2Lx4sWtio8xxiTj7S1cO0QEvPgicONuBXKnVxEiIsyYMQOffPIJbG1t0bNn\nT2zevBlvvfVWq95s5MiR6NGjh87j32rv3r2YNm0arK2t4ezsDBcXF2RkZKCkpARVVVXw8/MDIBTH\nhBszSvbt2ycWxilTpuDIkSOtik8OuM/dsnF+lq3N+b3zDtC3L5Cd3W665PS+RHfs2LE4ffo0hgwZ\nYvAg1q1bh61bt2Lo0KFYs2YNbG1tUVxcjBEjRojbqFQqaLVaWFtbQ6VSie2Ojo7QarUAAK1WCycn\nJwCAUqlE9+7dUV5ejp49e7Z4z4iICHHJIVtbW/j4+Iin0k3/kPg5P+fn/Nykz++/HylRUcAbbyDg\n7beBp59GSn6+JPE0Pb7dkm0GQ3qaOXMmZWRk6Lv5beXn55Onp6f4vLS0lBobG6mxsZGWLFlCs2fP\nJiKiuXPn0vbt28XtIiMjaffu3XTixAkaPXq02J6amkrjx48nIiJPT0/SarXia2q1msrKylrE0Iq0\nGWPM9MLC6BhASx54gN554glaEhxMxw4ckDQkY31v6n0mlJ6eju3bt6Nv3764//77AQjjOW29vUPv\n3r3Fx88//zwmTJgAQDjDKbzpCuKioiKoVCo4OjqiqKioRXvTPufPn0efPn1QX1+PyspKnWdBjDFm\nzlLHj8ehzz9H9O+/A8eOAQCW5OUBwJ93apUJvScmHDp0CHl5eTh69CgOHDiA/fv3Y9++fW0OoKSk\nRHy8Z88eceZcaGgo4uPjUVdXh/z8fOTm5sLPzw8ODg7o1q0bMjIyQETYtm0bJk6cKO6zZcsWAMDu\n3bsRGBjY5vgszc2n0nLE+Vk2zk8/ydu3I/qWsfLovDwcXrfOIMc3J3qfCTk7OyM7OxtpaWlQKBQY\nOXIkvL29W/Vm06ZNw7Fjx/D777/DyckJ7777LlJSUpCdnQ2FQoF+/frhww8/BAB4eHggLCwMHh4e\nUCqViIuLE1dpiIuLQ0REBKqrqzF27FiEhIQAACIjIzFjxgy4urrCzs4O8fHxrYqPMcbMgbK2Vme7\nVU2NiSMxPr3XjouNjcWGDRswefJkEBESEhLwwgsvYP78+caO0eB42R7GmDlbqtFgxU3XRjZZptHg\n/VsuczEVye8n5OXlhfT0dHE86OrVqxgxYgR++ukngwdlbFyEGGPmLDUxEYcWLED0jXEgAHirTx+E\nfPSRZGNCkq0d12zjDh10Pmbmg/vcLRvnZ9kMlZ//uHHQxMZimUaD5c7OWAYgRKWS3aQEoBVjQrNm\nzcLw4cObdcfNnj3bmLExxli75T9unFB0SkuBhx4SlvLJywPUaqlDM6hW3U8oMzNTvMndyJEj4evr\na7TAjIm74xhjFiUiAtiyBXjlFeGW4BIwi+44KysrKBQKKBQK7o5jjDFTaVqoedMm4b5DMqJ3JYmN\njcVzzz2HS5cu4eLFi3juueewdu1aY8bG7gH3uVs2zs+yGS0/X19g5EjhXkObNxvnPSSidxH6+OOP\nkZGRgffeew/vv/8+0tPTsWHDBmPGxhhjrEnT2dC6dUBjo7SxGFCrpmgfP34cnTt3BgBUV1fDz8+P\np2gzxpgp1NcLkxLOnwcOHABMPFPOWN+bPDuOMcYsgVIp3P77//5PuM2DTKZr690d99prr+GTTz5B\njx49YGdnh82bN+PVV181ZmzsHnCfu2Xj/Cyb0fN7/nngvvuAw4eBnBzjvpeJ6F2EwsPD0b9/fyxY\nsADz589H3759+UyIMcZMqUcPYOZM4bFMJobpPSbk4+OD7Ozsu7ZZAh4TYoxZrDNnAA8PoHNnoKgI\nMNHtaiS/ToiIUF5eLj4vLy9HQ0ODwQNijDF2BwMHAsHBQHU18PHHUkfTZnoXoYULF+KRRx7BsmXL\nsHTpUjzyyCNYtGiRMWNj94D73C0b52fZTJZf03Ttf/9bmDVnwfQuQjNnzsSXX36J3r17w8HBAXv2\n7MHMpr5JxhhjphMSAgwYABQWAgkJUkfTJq1aO04ueEyIMWbx/v1vYN484PHHgbQ0o7+d5PcTkhMu\nQowxi1dVBahUwlI+mZnA4MFGfTvJJyYwy8B97paN87NsJs2va1cgMlJ4HBtruvc1sFYVoYiICCxa\ntAgJCQkoLS01VkyMMcb0MXcuoFAA8fHCfYcsUKu7486cOYP09HSkp6cjMzMTYWFheP311y3q1g7c\nHccYk41Jk4C9e4Hly4F33jHa25jFmFB6ejqICI888ggA4PPPP4e3tzdSU1Px/PPPGzw4Y+EixBiT\njaNHgSefBOztgXPngE6djPI2ZjEm9PXXXyM1NRVTp07FrFmzcPr0aWi1Wtjb2xs8MHZvuM/dsnF+\nlk2S/AICgIcfFrrjPvvM9O/fRnqvog0AkyZNwrVr17B48WKx7eOPP4aTk5PBA2OMMaYHhQKYP19Y\n3DQ2FnjuOaHNQph0ivbs2bORmJiI3r17i/chKi8vx9SpU3Hu3Dk4Ozvjs88+g62tLQAgJiYGmzZt\ngpWVFdauXYvg4GAAQGZmJiIiIlBTU4OxY8ci9sbMkNraWsycORMnT56EnZ0ddu3ahb59+7aIg7vj\nGGOyUl2N1N69kXzlCpS+vqjv1QvB8+fD34C3ezCL7ri2mjVrFpKSkpq1rVq1CkFBQfj1118RGBiI\nVatWAQBycnKwa9cu5OTkICkpCXPmzBH/AFFRUdi4cSNyc3ORm5srHnPjxo2ws7NDbm4uXn311WZn\nbIwxJlep33yDQ0olVgBYnpWFFcnJOLRgAVITE6UO7a5MWoRGjhyJHj16NGvbt28fwsPDAQi3i0i4\nsQTF3r17MW3aNFhbW8PZ2RkuLi7IyMhASUkJqqqq4OfnB0BYTqhpn5uPNWXKFBw5csRUqZkN7nO3\nbJyfZZMqv+S1axH9xx/N2qLz8nB43TpJ4mmNVo0JGUNpaak4scHe3l68/qi4uBgjRowQt1OpVNBq\ntbC2toZKpRLbHR0dodVqAQBarVYcn1IqlejevTvKy8vRU8dS5xEREXB2dgYA2NrawsfHBwEBAQD+\n/IfEz/k5P+fnlvC86KZrhFJu/A4AYFVTc8/Hb3pcUFAAoyITy8/PJ09PT/G5ra1ts9d79OhBRERz\n586l7du3i+2RkZG0e/duOnHiBI0ePVpsT01NpfHjxxMRkaenJ2m1WvE1tVpNZWVlLWKQIG3GGDOa\nJcHBRECLn6UajcHew1jfm5JfYWpvb48LFy4AAEpKStC7d28AwhlOYWGhuF1RURFUKhUcHR1RVFTU\nor1pn/PnzwMA6uvrUVlZqfMsiDHG5CR4/nwsUaubtb2lViNo3jyJItKf5EUoNDQUW7ZsAQBs2bIF\nkyZNEtvj4+NRV1eH/Px85Obmws/PDw4ODujWrRsyMjJARNi2bRsmTpzY4li7d+9GYGCgNElJ6OZT\naTni/Cwb52cc/uPGQRMbi2UaDZZ37IhlAEIWLjTo7DhjMemY0LRp03Ds2DH8/vvvcHJywnvvvYc3\n3ngDYWFh2LhxozhFGwA8PDwQFhYGDw8PKJVKxMXFQXFj7ntcXBwiIiJQXV2NsWPHIiQkBAAQGRmJ\nGTNmwNXVFXZ2doiPjzdleowxJhn/ceOEovPss8DOnUKHnAXgWzkwxpicbNgAvPgiMGUKsHu3wQ5r\nFmvHyQUXIcaYbP3vf4CrK2BnB1y8CHQwzKiLLC5WZcbHfe6WjfOzbGaRn1ot3OyurAw4fVrqaO6K\nixBjjMmJQiGsqg0A33wjbSx64O44xhiTm82bgVmzgNBQ4V5DBsBjQgbERYgxJmvnzgHOzkD37kK3\nnJVVmw/JY0JML2bRJ21EnJ9l4/xMpG9foH9/oLISyMqSOpo74iLEGGNyNGqU8PvoUWnjuAvujmOM\nMTn69FPhBndjxgBffdXmw/GYkAFxEWKMyV5xMeDoCHTpApSXA9bWbTocjwkxvZhNn7SRcH6WjfMz\noT59ADc34MoV4MQJqaO5LS5CjDEmVxYwLsTdcYwxJleffw6EhQGjRwOHD7fpUDwmZEBchBhj7cKl\nS0Dv3oCNDfDHH0CnTvd8KB4TYnoxqz5pI+D8LBvnZ2K9egGenkBNDZCRIXU0OnERYowxOTPzcSHu\njmOMMTlLSACeegrw9weOHbvnw/CYkAFxEWKMtRsVFcK9haythXGhzp3v6TA8JsT0YnZ90gbG+Vk2\nzk8CPXoAPj5AXR3w3XdSR9MCFyHGGJO7pvsLmeG4EHfHMcaY3CUmAuPHA488cs9nQzwmZEBchBhj\n7crly0DPnsJdVysqhPXkWonHhJhezLJP2oA4P8vG+UmkWzdg6FCgvh749lupo2mGixBjjLUDqU5O\nWApg+YsvYqlGg9TERKlDAmBG3XHOzs7o1q0brKysYG1tjePHj6O8vBxTp07FuXPn4OzsjM8++wy2\ntrYAgJiYGGzatAlWVlZYu3YtgoODAQCZmZmIiIhATU0Nxo4di9jY2Bbvxd1xjLH2JDUxEYdeeAHR\nJSVi2xK1GprYWPiPG6fXMWTfHadQKJCSkoKsrCwcP34cALBq1SoEBQXh119/RWBgIFatWgUAyMnJ\nwa5du5CTk4OkpCTMmTNH/ONERUVh48aNyM3NRW5uLpKSkiTLiTHGzEHy2rXNChAAROfl4fC6dRJF\n9CezKUIAWlTZffv2ITw8HAAQHh6OhIQEAMDevXsxbdo0WFtbw9nZGS4uLsjIyEBJSQmqqqrg5+cH\nAJg5c6a4T3thtn3SBsL5WTbOTxrK2lqd7VY1NSaOpCWl1AE0USgUGD16NKysrPDSSy/hhRdeQGlp\nKezt7QEA9vb2KC0tBQAUFxdjxIgR4r4qlQparRbW1tZQqVRiu6OjI7Rarc73i4iIgLOzMwDA1tYW\nPj4+CAgIAPDnPyR+zs/5OT+Xw/O8q1fRJOXG7wAADTY2t92/6XFBQQGMisxEcXExERFdvHiRvL29\nKTU1lWxtbZtt06NHDyIimjt3Lm3fvl1sj4yMpN27d9OJEydo9OjRYntqaiqNHz++xXuZUdqMMWZ0\nxw4coLfUaiJA/HlTraZjBw7ofQxjfW+azZnQgw8+CADo1asXnnrqKRw/fhz29va4cOECHBwcUFJS\ngt69ewMQznAKCwvFfYuKiqBSqeDo6IiioqJm7Y6OjqZNhDHGzEzT5INlc+fCqqAADQMGIOSf/9R7\nUoIxmcWY0LVr11BVVQUAuHr1KpKTk+Hl5YXQ0FBs2bIFALBlyxZMmjQJABAaGor4+HjU1dUhPz8f\nubm58PPzg4ODA7p164aMjAwQEbZt2ybu017cfCotR5yfZeP8pOM/bhze/+tfsRzA++PGmUUBAsxk\nTKi0tBRPPfUUAKC+vh7Tp09HcHAwhg4dirCwMGzcuFGcog0AHh4eCAsLg4eHB5RKJeLi4qBQKAAA\ncXFxiIiIQHV1NcaOHYuQkBDJ8mKMMbNyozcJFy9KG8dNzOY6IVPi64QYY+3SwYPA2LFAcDBw6FCr\ndpX9dUKMMcaMzAzPhLgIyYw590kbAudn2Tg/ifXqJfy+dEnaOG7CRYgxxtqLpiJ08aIwUdsM8JgQ\nY4y1J127AleuCLf67t5d7914TIgxxljbmdm4EBchmTH7Puk24vwsG+dnBm7ukjMDXIQYY6w9aToT\nMpPJCTwmxBhj7cnzzwMbNwIffgi8+KLeu/GYEGOMsbYzs2naXIRkxiL6pNuA87NsnJ8Z4IkJjDHG\nJGNmZ0I8JsQYY+1JcjKg0QCBgcDXX+u9G48JMcYYazueos2MySL6pNuA87NsnJ8ZMLMp2lyEGGOs\nPXngAeH3pUtAY6O0sYDHhBhjrP2xtQUqK4GyMqBnT7124TEhxhhjhmFG07S5CMmMRfRJtwHnZ9k4\nPzNhRtO0uQgxxlh7Y0ZnQjwmxBhj7c2LLwIbNgBxcUBUlF678JgQY4wxwzCjadpchGTGYvqk7xHn\nZ9k4PzNhRhesyrIIJSUlwd3dHa6urli9erXU4ZhUdna21CEYFedn2Tg/M8FnQsbT0NCAuXPnIikp\nCTk5Odi5cyfOnDkjdVgm88cff0gdglFxfpaN8zMPqXl5WApgeXIylmo0SE1MlCwWpWTvbCTHjx+H\ni4sLnJ2dAQB/+ctfsHfvXgwcOFDawBhjzAykJibi0EcfIRoA/vgDSE7Gkrw8AID/uHEmj0d2Z0Ja\nrRZOTk7ic5VKBa1WK2FEplVQUCB1CEbF+Vk2zk96yWvXIrqwsFlbdF4eDq9bJ0k8spui/cUXXyAp\nKQkbNmwAAGzfvh0ZGRlYd9MfWKFQSBUeY4xZLGOUC9l1xzk6OqLwpipfWFgIlUrVbBuZ1V3GGLNY\nsuuOGzp0KHJzc1FQUIC6ujrs2rULoaGhUofFGGNMB9mdCSmVSvz73/+GRqNBQ0MDIiMjeVICY4yZ\nKdmcCa1ZswYdOnRAeXk5xowZg7NnzyIyMhKbNm2Cu7s7kpOTxW0zMzPh5eUFV1dXLFiwQGyvra3F\n1KlT4erqihEjRuDcuXPia1u2bMGAAQMwYMAAbN261WR5LVu2DN7e3vDx8UFgYGCzrsaYmBi4urpa\ndH6LFi3CwIED4e3tjcmTJ6OyslJ8TQ75ff755xg0aBCsrKxw8uTJZq/JIT99WdK1e7Nnz4a9vT28\nvLzEtvLycgQFBWHAgAEIDg5uNhXbkJ+jsRUWFmLUqFEYNGgQPD09sXbtWunzIxk4f/48aTQacnZ2\nprKyMiIi+vnnn8nb25vq6uooPz+f1Go1NTY2EhHRsGHDKCMjg4iIxowZQwcPHiQiovXr11NUVBQR\nEcXHx9PUqVOJiKisrIz69+9PFRUVVFFRIT42hcuXL4uP165dS5GRkbLKLzk5mRoaGoiIaPHixbR4\n8WJZ5XfmzBk6e/YsBQQEUGZmptgul/z0UV9fT2q1mvLz86muro68vb0pJydH6rBuKzU1lU6ePEme\nnp5i26JFi2j16tVERLRq1Sqj/Ds1hZKSEsrKyiIioqqqKhowYADl5ORImp8sitDTTz9NP/74Y7Mi\ntHLlSlq1apW4jUajoe+//56Ki4vJ3d1dbN+5cye99NJL4jbp6elERHT9+nV64IEHiIhox44d9PLL\nL4v7vPTSS7Rz506j53WrlStXiv845Jjfl19+SdOnTyci+eV3axGSW3538t1335FGoxGfx8TEUExM\njIQR3V1+fn6zIuTm5kYXLlwgIuGL3M3NjYgM+zlKYeLEiXT48GFJ87P47ri9e/dCpVLh4YcfbtZe\nXFzcbFZc0/VCt7Y7OjqK1xHdfI2RUqlE9+7dUVZWdttjmcqSJUvw0EMPYfPmzXjzzTcByCu/Jps2\nbcLYsWMByDO/m8k9v5vJ4dq90tJS2NvbAwDs7e1RWloKwHCfY3l5ualSERUUFCArKwvDhw+XND+L\nmJgQFBSECxcutGiPjo5GTExMs35KssDp17fLb+XKlZgwYQKio6MRHR2NVatW4ZVXXsEnn3wiQZT3\n7m75AcJn2bFjRzz77LOmDq/N9MmvPZPbdXkKhcLic7py5QqmTJmC2NhYdO3atdlrps7PIorQ4cOH\ndbafPn0a+fn58Pb2BgAUFRVhyJAhyMjIaHG9UFFREVQqFRwdHVFUVNSiHRCq+fnz59GnTx/U19ej\nsrISdnZ2cHR0bLY6bmFhIZ588kmj53erZ599VjxTkFN+mzdvxldffYUjR46IbXLKTxdLyq+t9Ll2\nz9zZ29vjwoULcHBwQElJCXrfWADUUJ9jz549TZbL9evXMWXKFMyYMQOTJk2SPj9D9i9KTdfEhNra\nWvrtt9+of//+4oCan58fpaenU2NjY4sBtaa+9Z07dzYb+O3Xrx9VVFRQeXm5+NgUfv31V/Hx2rVr\n6bnnnpNVfgcPHiQPDw+6dOlSs3a55NckICCATpw4Idv87uT69evUv39/ys/Pp9raWrOfmEDUckxo\n0aJF4thITExMi4F7Q3yOptDY2EgzZsygV155pVm7lPnJqgj169dPLEJERNHR0aRWq8nNzY2SkpLE\n9hMnTpCnpyep1WqaN2+e2F5TU0PPPPMMubi40PDhwyk/P198bdOmTeTi4kIuLi60efNmk+RDRDRl\nyhTy9PQkb29vmjx5MpWWloqvySE/FxcXeuihh8jHx4d8fHzEWTVE8sjvyy+/JJVKRTY2NmRvb08h\nISHia3LIT19fffUVDRgwgNRqNa1cuVLqcO7oL3/5Cz344INkbW1NKpWKNm3aRGVlZRQYGEiurq4U\nFBTUrMgb8nM0trS0NFIoFOTt7S3+N3fw4EFJ85Pd2nGMMcYsh8XPjmOMMWa5uAgxxhiTDBchxhhj\nkuEixBhjTDJchBhjjEmGixBjjDHJcBFizASWL1+ONWvWiM8fe+yxezpOZWUlPvjgg1btU1BQgM6d\nO2Pw4MHi85tvU9BaNTU18PHxQadOnSRZ84zJCxchxgyMhIvAm7XduhbXf//733s6dkVFBeLi4lq9\nn4uLS4v7Gd0rGxsbZGdno0+fPgY5HmvfuAixdm/r1q3ijQNnzpwptv/zn/+El5cXvLy8EBsbe8f2\ngoICuLm5ITw8HF5eXigsLER0dDTc3NwwcuRInD17ttl7dunSBefOncPAgQPx4osvwtPTExqNBjU1\nNeI2Tz31FIYOHQpPT09s2LABAPDGG28gLy8Pvr6+WLx4MQBg+/btGD58OHx9ffHyyy+jsbFR79x/\n++03DB48GJmZmSgoKIC7uztmzZoFNzc3TJ8+HcnJyXjssccwYMAA/PDDD63/4zJ2NwZdE4IxM1Nf\nX0+ffvopvf/++7R582aaM2cO5eXlia+fPn2aBgwYIC73VF5eTkTCkiReXl507do1unLlCg0aNIiy\nsrJu256fn08dOnQQb/LVtF11dTVdvnyZXFxcaM2aNeL7dunShQoKCkipVNKPP/5IRERhYWG0fft2\ncZumWK5du0aenp5UXl5OBQUFzdY0y8nJoQkTJlB9fT0REUVFRdHWrVub/Q1uXQet6fkvv/xCvr6+\ndOrUKbFdqVTS6dOnqbGxkYYMGUKzZ88mIqK9e/fSpEmTmh335rUaGbtXFrGKNmP36scff8SUKVPw\nxRdfoLa2Fs888wwefPBB8fVvvvkGYWFh4iq/PXr0AAB8++23mDx5Mjp37gwAmDx5MtLS0kBEOttD\nQ0PRt29f+Pn5AQDS0tIwefJk2NjYwMbGBqGhoTpvM9KvXz/xXlhDhgxBQUGB+FpsbCwSEhIACKsU\n5+bmiqsbNzly5AgyMzMxdOhQAEB1dTUcHBzu+ne5ePEiJk2ahD179sDd3b1ZPIMGDQIADBo0CKNH\njwYAeHp6NouNMUPhIsRkrWkw/vvvv8drr72Gfv36NXtdoVDoLA63tt/pcdN4z/3336/X/jfr1KmT\n+NjKygrV1dUAgJSUFBw5cgTp6emwsbHBqFGjmnXV3Sw8PBwrV67U+drt2Nraom/fvkhLS2tWhG6O\np0OHDujYsaP4uL6+vlXvwZg+eEyIydoPP/yA33//HadPn0a/fv2QlpbW7PUnn3wSn3/+uTjLq6Ki\nAgAwcuRIJCQkoLq6GlevXkVCQgL8/f11to8cObJFkfH390dCQgJqampQVVWFAwcOtOpGYZcvX0aP\nHj1gY2ODX375Benp6QCArl27oqqqStwuMDAQu3fvxqVLlwAA5eXlOH/+/F2P37FjR3z55ZfYunUr\ndu7cqXdcjBkanwkxWUtKSoK9vT0ee+wx7NmzBw888ECz1z08PLBkyRI88cQTsLKywuDBg7Fp0yb4\n+voiIiJC7F574YUXxJsn6movKChoVmR8fX0xdepUeHt7o3fv3uL2TZq2vbUwNT0PCQnBf/7zH3h4\neMDNzQ2PPPIIAMDOzg6PPfYYvLy8MHbsWKxevRorVqxAcHAwGhsbYW1tjbi4ODz00EN3/LsoFArc\nd999OHDgAIKCgtC1a1d4eXndNh5dsTJmCHwrB8ZkrqCgABMmTMBPP/1k0OP269cPmZmZJr0rKJMf\n7o5jTOaUSiUqKyvF8bG2arpYtb6+Hh068FcIaxs+E2KMMSYZ/t8YxhhjkuEixBhjTDJchBhjjEmG\nixBjjDHJcBFijDEmGS5CjDHGJMNFiDHGmGS4CDHGGJPM/wekinOacVQ+iwAAAABJRU5ErkJggg==\n"
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#====================================================================\n",
"#Solving for different times\n",
"#====================================================================\n",
"#Time array\n",
"times = np.arange( t0, tf, tstep )\n",
"\n",
"rpos = []\n",
"for t in times:\n",
" #Finding the new eccentric anomaly\n",
" E = Bisection( kepler, a0, b0, Niter )\n",
" #Computing coordinates at this time\n",
" ri = r(E)\n",
" print \"In %f hours, the satellite is located at (%f,%f) km\"%(t/HR2SC, ri[0], ri[1])\n",
" rpos.append( ri )\n",
"rpos = np.array(rpos) \n",
"\n",
"#Plotting\n",
"plt.plot( rpos[:,0], rpos[:,1], \"o-\", color=\"red\", lw = 2 )\n",
"plt.grid(True)\n",
"plt.xlabel(\"$x$ coordinate [km]\")\n",
"plt.ylabel(\"$y$ coordinate [km]\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- - - "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fixed-point Iteration"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Although in many cases the use of Bisection is more than enough, there are some pathological situations where the use of more advanced methods is required.\n",
"One of the advantages featured by this method is one does not have to give an interval where the solution is within, instead, from a seed the algorithm will converge towards the required solution."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Steps FP"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. Take your function $f(x)$ and rewrite it like \n",
"\n",
" $f(x) = x - g(x)$\n",
"\n",
"2. Give a guest to the solution (root of $f(x)$). This value would be the seed $p_0$.\n",
"3. The next guest to the solution will be given by\n",
"\n",
" $p_2 = g(p_1)$\n",
"\n",
"4. If the stop condition is not satisfied, then repeat step 3.\n",
"5. The End!"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Defining Fixed-point iteration function\n",
"def FixedPoint_Animation( f, pini, Nmax, xmin, xmax ):\n",
" g = lambda x: x-f(x)\n",
" #Initial condition\n",
" pi = [pini,]\n",
" px = [pini,pini,]\n",
" py = [0,]\n",
" #Iterations\n",
" for n in xrange(Nmax+3):\n",
" pi.append( g(pi[n]) )\n",
" px.append( g(pi[n]) )\n",
" px.append( g(pi[n]) )\n",
" py.append( g(pi[n]) )\n",
" py.append( g(pi[n]) )\n",
" \n",
" py.append( g(pi[n+1]) )\n",
" pi = np.array( pi )\n",
" px = np.array( px )\n",
" py = np.array( py )\n",
" \n",
" print \"Result:\", pi[-1]\n",
" \n",
" #Array X-axis\n",
" X = np.linspace(xmin,xmax,100)\n",
" \n",
" #Initializing Figure\n",
" fig = plt.figure( figsize=(7,7) )\n",
" ax = fig.add_subplot(111)\n",
" #Graphic iterations\n",
" fixedpoint, = ax.plot( [], [], color=\"red\", linewidth = 3 )\n",
" #Function g\n",
" ax.plot( X, g(X), color=\"green\", linewidth = 2 )\n",
" #Identity funcion\n",
" ax.plot( X, X, color=\"blue\", linewidth = 2 )\n",
" ax.grid(True)\n",
" ax.set_xlim( (xmin, xmax) )\n",
" ax.set_ylim( (xmin, xmax) )\n",
" ax.set_xlabel( \"X axis\" )\n",
" ax.set_ylabel( \"Y axis\" )\n",
" ax.set_title( \"Fixed-Point iteration\" )\n",
" \n",
" def init():\n",
" fixedpoint.set_data([], [])\n",
" return fixedpoint,\n",
" \n",
" def animate(i):\n",
" #Setting new data\n",
" fixedpoint.set_data( px[:2*i], py[:2*i] )\n",
" ax.set_title( \"Fixed-Point. Iteration %d\"%i )\n",
" return fixedpoint,\n",
" \n",
" return animation.FuncAnimation(fig, animate, init_func=init,frames=Nmax, interval=500, blit=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Find the roots of the next functions:\n",
"\n",
"$$f_1(x) = \\frac{x^2-1}{3}$$\n",
"\n",
"and \n",
"\n",
"$$f_2(x) = x-\\cos x$$\n",
"\n",
"using Fixed-Point iteration."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Result: 0.999999989073\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"FixedPoint_Animation( f, pini = 0.4, Nmax = 50, xmin = -0.5, xmax = 0.5 )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ACTIVITY FP"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When a new planet is discovered, there are different methods to estimate its physical properties. Many times is only possible to estimate either the planet mass or the planet radius and the other property has to be predicted through computer modelling.\n",
"\n",
"If one has the planet mass, a very rough way to estimate its radius is to assume certain composition (mean density) and a homogeneous distribution (a very bad assumption!). For example, for the planet [Gliese 832c](http://es.wikipedia.org/wiki/Gliese_832_c) with a mass $M= 5.40 M_{\\oplus}$, if we assume an earth-like composition, i.e. $\\bar \\rho_{\\oplus} = 5520\\ kg/m^3$, we obtain:\n",
"\n",
"$$R_{g832c} = \\left( \\frac{3 M_{g832c}}{ 4 \\pi \\bar\\rho_{\\oplus} } \\right)^{1/3} \\approx 1.75 R_{\\oplus}$$\n",
"\n",
"That would be the planet radius if the composition where exactly equal to earth's.\n",
"\n",
"A more realistic approach is assuming an internal one-layer density profile like:\n",
"\n",
"$$\\rho(r) = \\rho_0 \\exp\\left( -\\frac{r}{L} \\right)$$\n",
"\n",
"where $\\rho_0$ is the density at planet centre and $L$ is a characteristic lenght depending on the composition. From numerical models of planet interiors, the estimated parameters for a planet of are $M= 5.40 M_{\\oplus}$ are approximately $\\rho_0 = 18000\\ kg/m^3$ and $L = 6500\\ km$.\n",
"\n",
"Integrating over the planet volume, we obtain the total mass as\n",
"\n",
"$$M = 4\\pi \\int_0^R \\rho(r)r^2dr$$\n",
"\n",
"This is a function of the mass in terms of the planet radius. \n",
"\n",
"Solving the equation $M(R) = M_{g832c}$ it would be possible to find a more realistic planet radius. However when using numerical models, it is not possible to approach the solution from the left side as a negative mass makes no sense."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" \n",
"In an IPython notebook, solve the previous problem and find the radius of **Gliese 832c** using your own version of the Fixed-point iteration algorithm.\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- - -"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Newton-Raphson Method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Although Fixed-point iteration is an efficient algorithm as compared with Bisection, the Newton-Raphson method is an acceletared convergent scheme where the roots of a function are easily found with just a few iterations."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Derivation NM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Although this method can be presented from an algorithmic point of view, the mathematical deduction is very useful as it allows us to understand the essence of the approximation as well as estimating easily the convergence errors.\n",
"\n",
"Let be $f(x)$ a continuous and differentiable function defined within an interval $[a,b]$ (i.e. $f\\in \\mathcal{C}^2[a,b]$), and $p$ is a root of the function such that $f(p) = 0$. If we give an initial an enough close guess $p_0$ to this root, such that $|p-p_0|<\\epsilon$, where $\\epsilon$ is adequately small, we can expand the function by using a second order taylor serie, yielding:\n",
"\n",
"$$f(p) = f(p_0) + (p-p_0)f'(p_0) + \\frac{(p-p_0)^2}{2}f''(p_0) + \\mathcal{O}^3(|p-p_0|)$$\n",
"\n",
"but as $f(p) = 0$ and $|p-p_0|^2<\\epsilon^2$ is an even smaller quantity, we can readily neglect from second order terms, obtaining\n",
"\n",
"$$p \\approx p_0 - \\frac{f(p_0)}{f'(p_0)} \\equiv p_1$$\n",
"\n",
"If we repeat this process but now using $p_1$ as our guess to the root instead of $p_0$ we shall obtain:\n",
"\n",
"$$p \\approx p_1 - \\frac{f(p_1)}{f'(p_1)} \\equiv p_2$$\n",
"\n",
"and so...\n",
"\n",
"$$p \\approx p_n - \\frac{f(p_n)}{f'(p_n)} \\equiv p_{n+1}$$\n",
"\n",
"where each new iteration is a better approximation to the real root."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Steps NM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. Take your function $f(x)$ and derive it, $f'(x)$.\n",
"2. Give a guest to the solution (root of $f(x)$). This value would be the seed $p_0$.\n",
"3. The next guest to the solution will be given by\n",
"\n",
" $$p_{n+1} = p_n - \\frac{f(p_n)}{f'(p_{n+1})}$$\n",
"\n",
"4. If the stop condition is not satisfied, then repeat step 3.\n",
"5. The End!"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Defining Newton Method\n",
"def NewtonRaphson_Animation( f, fp, pini, Nmax, xmin, xmax ):\n",
" #Initial condition\n",
" p = [pini,]\n",
" p_dash = []\n",
" p_der = []\n",
" #Iterations\n",
" for n in xrange(Nmax):\n",
" p.append( p[n] - f(p[n])/fp(p[n]) )\n",
" p_dash.append( p[n] )\n",
" p_dash.append( p[n] )\n",
" p_der.append( 0 )\n",
" p_der.append( f(p[n]) )\n",
" \n",
" p = np.array( p )\n",
" p_dash = np.array( p_dash )\n",
" p_der = np.array( p_der )\n",
" \n",
" print \"Result:\", p[-1]\n",
" \n",
" #Array X-axis\n",
" X = np.linspace(xmin,xmax,100)\n",
" \n",
" #Initializing Figure\n",
" fig = plt.figure( figsize=(7,7) )\n",
" ax = fig.add_subplot(111)\n",
" #Graphic iterations\n",
" dash, = ax.plot( [], [], \"--\", color=\"gray\", linewidth = 2 )\n",
" derivative, = ax.plot( [], [], color=\"red\", linewidth = 3 )\n",
" #Function f\n",
" ax.plot( X, f(X), color=\"green\", linewidth = 2 )\n",
" #Horizontal line\n",
" ax.hlines( 0, xmin,xmax, color=\"black\", lw = 2 )\n",
" ax.grid(True)\n",
" ax.set_xlim( (xmin, xmax) )\n",
" ax.set_xlabel( \"X axis\" )\n",
" ax.set_ylabel( \"Y axis\" )\n",
" ax.set_title( \"Fixed-Point iteration\" )\n",
" \n",
" def init():\n",
" dash.set_data([], [])\n",
" derivative.set_data([], [])\n",
" return dash, derivative\n",
" \n",
" def animate(i):\n",
" #Setting new data\n",
" dash.set_data( p_dash[:2*i+2], p_der[:2*i+2] )\n",
" derivative.set_data( p_dash[2*i+1:2*i+3], p_der[2*i+1:2*i+3] )\n",
" ax.set_title( \"Newthon-Raphson Method. Iteration %d\"%i )\n",
" return dash, derivative\n",
" \n",
" return animation.FuncAnimation(fig, animate, init_func=init,frames=Nmax, interval=500, blit=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Find one root of the function:\n",
"\n",
"$f(x) = x^2 - x$\n",
"\n",
"with derivative\n",
"\n",
"$f'(x) = 2x -1$\n",
"\n",
"using the Newton-Raphson method."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Result: 1.0\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#Defining the function\n",
"def f(x): \n",
" return x**2-x\n",
"#Defining the derivative\n",
"def df(x): \n",
" return 2*x-1\n",
"#Calculating root\n",
"NewtonRaphson_Animation( f, df, pini = 0.55, Nmax = 10, xmin = 0, xmax = 3.5 )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Stop conditions NM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The stop conditions are the same than Bisection and Fixed-point iteration:\n",
"\n",
"* A fixed distance between the last two steps (absolute convergence):\n",
"\n",
" $$|p_i - p_{i-1}|<\\epsilon$$\n",
"\n",
"* A fixed relative distance between the last two steps (relative convergence):\n",
"\n",
" $$\\frac{|p_i - p_{i-1}|}{|p_i|}<\\epsilon\\ \\ \\ \\ \\ p_i \\neq 0$$\n",
"\n",
"* Function tolerance:\n",
"\n",
" $$f(p_i)< \\epsilon$$\n",
"\n",
"* Computational stop:\n",
"\n",
" If $N>N_{max}$, stop!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Convergence NM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is possible to demonstrate by means of the previous derivation procedure, that the convergence of the Newton-Raphson method is quadratic, i.e., if $p$ is the exact root and $p_n$ is the $n$-th iteration, then\n",
"\n",
"$$|p_{n+1}-p|\\leq C |p_n-p|^2$$\n",
"\n",
"for a fixed ans positive constant $C$.\n",
"\n",
"This implies, if the initial guess is good enough such that |p_0-p| is small, the convergence is achieved very fast as each iteration improves the precision twice in the order of magnitude, e.g., if $|p_0-p|\\sim 10^{-1}$, $|p_1-p|\\sim 10^{-2}$, $|p_2-p|\\sim 10^{-4}$, $|p_2-p|\\sim 10^{-8}$ and so."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Defining Newton Method\n",
"def NewtonRaphson( f, fp, pini, Nmax ):\n",
" #Initial condition\n",
" p = pini\n",
" #Iterations\n",
" for n in xrange(Nmax):\n",
" p = p - f(p)/fp(p)\n",
" #Final result\n",
" return p"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"**ACTIVITY**\n",
" \n",
"In an IPython notebook, copy the latter routine NewtonRaphson and the Bisection routine provided in this notebook (*and if you have it already, the Fixed-point iteration routine as well*) and find the root of the next function using all the methods.\n",
" \n",
"$f(x) = x - \\cos(x)$\n",
" \n",
"Plot in the same figure the convergence of each method as a function of the number of iterations.\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- - -"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Secant Method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Newton-Raphson method is highly efficient as the convergence is accelerated, however there is a weakness with it: one needs to know the derivative of the function beforehand. This aspect may be complicated when dealing with numerical functions or even very complicated analytical functions. Numerical methods to derive the input function can be applied, but this extra procedure may involve an extra computing time that compensates the time spent by using other methods like Bisection."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Derivation SM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Retaking the iterative expression obtained from the Newton-Raphson method:\n",
"\n",
"$$p_{n+1} = p_n - \\frac{f(p_n)}{f'(p_{n+1})}$$\n",
"\n",
"the derivative can be approximated as\n",
"\n",
"$$f'(p_n) = \\lim_{x\\rightarrow p_n} \\frac{f(x)-f(p_n)}{x-p_n} $$\n",
"\n",
"As we know, the convergence of the NR method is quadratic, so $p_{n-1}$ should be close enough to $p_n$ such that one can assume $p_{n-1}\\rightarrow p_n$ and the previous term is:\n",
"\n",
"$$f'(p_n) \\approx \\frac{f(p_{n})-f(p_{n-1})}{p_{n}-p_{n-1}} $$\n",
"\n",
"The final expression for the $n$-th iteration of the root is then:\n",
"\n",
"$$p_n = p_{n-1} - \\frac{ f(p_{n-1})(p_{n-1}-p_{n-2}) }{f(p_{n-1})-f(p_{n-2})}$$\n",
"\n",
"In this consists the Secant method, what is just an approximation to the Newton-Raphson method, but without the derivative term."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Steps SM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. Give the input function $f(x)$.\n",
"2. Give two guests to the solution (root of $f(x)$). These values would be the seeds $p_0$, $p_1$.\n",
"3. The next guest to the solution will be given by\n",
"\n",
" $$p_n = p_{n-1} - \\frac{ f(p_{n-1})(p_{n-1}-p_{n-2}) }{f(p_{n-1})-f(p_{n-2})}$$\n",
"\n",
"4. If the stop condition is not satisfied, then repeat step 3.\n",
"5. The End!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"**ACTIVITY**\n",
" \n",
"In an IPython notebook and based on the routine NewtonRaphson, write your own routine SecantMethod that performs the previous steps for the Secant Method. Test your code with the function $f(x)$:\n",
" \n",
"$f(x) = x - \\cos(x)$\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**ACTIVITY**\n",
"\n",
"\n",
"It is known that light rays are deflected when they pass near by a gravitational field and that this deviation is proportional to the body mass which the light is interacting with and inversely proportional to the passing distance. \n",
"Since it is common finding very massive structures in the universe and the measures that are done to study it involve photons, it makes sense to study what happens to a light source image when the rays get close to a grumpy object like a dark matter halo. \n",
"\n",
"\n",
"In order to study the light deflection in these cases, it will be used the simplest model, gravitational lens theory, where the len is a very massive object. A sketch of a typical system is shown in the figure below. The source plane is the light source or image that is going to be affected, $\\eta$ is the distance from a image point to the line of sight and $\\beta$ the subtended angle by the point. \n",
"The lens plane corresponds to the mass that affects the light coming from the source, $\\xi$ is the new image point distance to the line of sight, $\\theta$ is the subtended angle by the new point position. Then, $\\alpha$ is the deflection angle. \n",
"\n",
"Since from observations $\\theta$ is known, the problem to be solved per pixel usually is \n",
"\n",
"\\begin{equation}\n",
"\\beta = \\theta - \\hat{\\alpha}(\\theta) \n",
"\\end{equation}\n",
"\n",
"but $\\alpha$ also depends on $\\theta$ besides the len halo properties. This would allow construct the real image\n",
"from the distorted and magnified one. \n",
"\n",
"\n",
"\n",
"This equation can also be written in terms of distances \n",
"\n",
"\\begin{equation}\n",
"\\vec{\\eta} = \\frac{D_s}{D_d} \\vec{\\xi} - D_{ds}\\alpha ( \\vec{\\xi }) \n",
"\\end{equation}\n",
"\n",
"The solution to the lens equation is easier to get if it is assumed that the len is axially symmetric. In this case, the deflection angle takes the next form\n",
"\n",
"$$ \\hat{\\alpha}(\\vec{\\xi}) = \\frac{\\vec{\\xi}}{|\\vec{\\xi}|^2} \\frac{8G\\pi}{c^2} \\int_0^\\xi d\\xi'\\xi'\\Sigma(\\xi')$$\n",
"\n",
"The quantity $\\Sigma$ is the surface mass density, i.e., the len's mass enclosed inside $\\xi$ circle per area unit. \n",
"It is important to notice that the direction of $\\alpha$ is the same as $\\xi$ and consequently $\\eta$. \n",
"\n",
"\n",
"The problem to be solved is the next: Given the positions of a square find the image distorsion due to gravitational lensing, i.e., find the root of \\xi in the trascendal equation it satisfies. Use the routines given below and all of \n",
"the data for the len and image that is going to be distorted. \n"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#Superficial density of the lens\n",
"def Sup_density(radio):\n",
" return radio*M*len(r[r"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"xc = epsilon*np.cos(theta)\n",
"yc = epsilon*np.sin(theta)\n",
"plt.figure(figsize=(12,6))\n",
"plt.subplot(1,2,1)\n",
"plt.plot(xc,yc,'b.');\n",
"plt.subplot(1,2,2)\n",
"plt.plot(X0,Y0,\"b.\");"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.10"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
![](\"http://upload.wikimedia.org/wikipedia/commons/thumb/8/8c/Bisection_method.svg/300px-Bisection_method.svg.png\")
\n",
"![](\"./figures/lente1.png\")
\n",
"\n",
"This equation can also be written in terms of distances \n",
"\n",
"\\begin{equation}\n",
"\\vec{\\eta} = \\frac{D_s}{D_d} \\vec{\\xi} - D_{ds}\\alpha ( \\vec{\\xi }) \n",
"\\end{equation}\n",
"\n",
"The solution to the lens equation is easier to get if it is assumed that the len is axially symmetric. In this case, the deflection angle takes the next form\n",
"\n",
"$$ \\hat{\\alpha}(\\vec{\\xi}) = \\frac{\\vec{\\xi}}{|\\vec{\\xi}|^2} \\frac{8G\\pi}{c^2} \\int_0^\\xi d\\xi'\\xi'\\Sigma(\\xi')$$\n",
"\n",
"The quantity $\\Sigma$ is the surface mass density, i.e., the len's mass enclosed inside $\\xi$ circle per area unit. \n",
"It is important to notice that the direction of $\\alpha$ is the same as $\\xi$ and consequently $\\eta$. \n",
"\n",
"\n",
"The problem to be solved is the next: Given the positions of a square find the image distorsion due to gravitational lensing, i.e., find the root of \\xi in the trascendal equation it satisfies. Use the routines given below and all of \n",
"the data for the len and image that is going to be distorted. \n"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#Superficial density of the lens\n",
"def Sup_density(radio):\n",
" return radio*M*len(r[r