{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integration of the Kepler problem with Taylor's method in `Julia`\n",
    "\n",
    "#### L. Benet, Instituto de Ciencias Físicas, UNAM\n",
    "benet < AT > fis.unam.mx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I. The Taylor Integrator\n",
    "\n",
    "-------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Taylor's integration method for ODEs is quite powerful, allowing to reach a precision comparable to round-off errors in one time-step. While the equations of motion can be in general easily implemented, an efficient implementation requires additional work. Below, we shall show two distinct implementations of the equations of motion for the Kepler problem, which shows how to optimize the running-time; see [TaylorIntegration.jl](https://github.com/PerezHz/TaylorIntegration.jl) for a yet non-optimized implementation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE** \n",
    "\n",
    "Below we use `julia` version 0.4.7 (0.4.8-pre+1) as the kernel, though it works using Julia v0.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4.8-pre+1\n"
     ]
    }
   ],
   "source": [
    "println(VERSION)\n",
    "\n",
    "using Compat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We load first the relevant package: [TaylorSeries.jl][1]\n",
    "\n",
    "[1]: http://github.com/lbenet/TaylorSeries.jl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: New definition \n",
      "    /(TaylorSeries.Taylor1{#T<:Real}, #T<:Real) at /Users/benet/.julia/v0.4/TaylorSeries/src/Taylor1.jl:261\n",
      "is ambiguous with: \n",
      "    /(TaylorSeries.Taylor1{Base.Rational{#T<:Integer}}, #S<:Union{Base.Complex, Real}) at /Users/benet/.julia/v0.4/TaylorSeries/src/Taylor1.jl:254.\n",
      "To fix, define \n",
      "    /(TaylorSeries.Taylor1{_<:Base.Rational{#T<:Integer}}, _<:Base.Rational{#T<:Integer})\n",
      "before the new definition.\n"
     ]
    }
   ],
   "source": [
    "using TaylorSeries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following *parameters* are set for the integrator:\n",
    "\n",
    "- `ordenTaylor`: Order of the Taylor polynomials considered\n",
    "- `epsAbs`: Absolute value set for the last (*and* the one-before-last) term of the Taylor expansion of the dynamical variables. This value is used to define an integration step and represents a measure of the error. Notice below that this value is smaller than `eps(1.0)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Taylor order = 28\n",
      " Eps = 1.0e-20\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Parámetros para el integrador de Taylor\n",
    "const _ordenTaylor = 28\n",
    "const _epsAbs = 1.0e-20\n",
    "\n",
    "println(\" Taylor order = $_ordenTaylor\\n Eps = $_epsAbs\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Taylor method works as follows: The equations of motion and the initial conditions $x(t_0)$, $y(t_0)$, $v_x(t_0)$, $v_y(t_0)$, are used to obtain *recursively* each term of the Taylor expansion of the dynamical variables, exploiting the relation\n",
    "$$\n",
    "x_{[n+1]} = \\frac{f_{[n]}(x,t)}{n+1}.\n",
    "$$\n",
    "Here, $f_{[n]}(x,t)$ is the $n$-th (normalized) Taylor coefficient of $f(x,t)$ in terms of the independent variable $t$, where $f(x,t)$ is the rhs of the equation of motion $\\dot{x} = f(x,t)$. Likewise, $x_{[n]}$ is the $n$-th Taylor coefficient for the dependent variable $x(t)$, expanded around $x(t_0)=x_{[0]}$. The latter corresponds to the initial condition.\n",
    "\n",
    "Once all Taylor coefficients are obtained up to a maximum order of the Taylor expansion, which is large enough to ensure convergece of the series, the last two terms of the expansion are used to determine the step size $h$ for the integration together with the value of `EpsAbs`. (Other methods exist that yield better optimized step-sizes, but these are usually more involved to compute.) Evaluating the Taylor expansion with the step-size yields the actual values of the dynamical variables at $t_1=t_0+h$. These values are then used as new *initial conditions* (at $t_0+h$) and everything is iterated.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Except for the actual implementation of the equations of motion (the *jet*) discussed below, the following functions do what we just described:\n",
    "\n",
    "- `taylorStepper`: carries one step of the integration from $t_0$ to $t_1=t_0+h$, returning h and the value of the dynamical variables evaluated at $t_1$. This routine depends on a user-supplied routine that performs the computation of the Taylor coefficients of all dynamical variables, i.e., that implements the calculation of the jet (`jetEqs`). This routine has as input value a vector with the dynamical variables.\n",
    "- `stepsize`: Returns the integration step ($h$) from the Taylor expansion coefficients of *one* dynamical variable and `epsAbs`, as given by \n",
    "$ h = \\min[\\, (\\epsilon/x^{[p-1]})^{1/(p-1)}, (\\epsilon/x^{[p]})^{1/p}\\, ], $\n",
    "where $p$ is the order of the Taylor polynomial, $x^{[r]}$ is the $r$-order Taylor coefficient (the $r$-th derivative divided by $r!$), and $\\epsilon={\\tt epsAbs}$. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "taylorStepper (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Returns stepsize of the integration and a vector with the updated values of the dependent\n",
    "# variables\n",
    "function taylorStepper{T<:Real}( jetEqs::Function, vec0::Array{T,1} )\n",
    "    \n",
    "    n = length( vec0 )\n",
    "\n",
    "    vec0T = Array(Taylor1{T},n)\n",
    "    @simd for i in eachindex(vec0)\n",
    "        @inbounds vec0T[i] = Taylor1([vec0[i]], _ordenTaylor)\n",
    "    end\n",
    "\n",
    "    # Jets\n",
    "    vec1T = jetEqs( vec0 )\n",
    "    \n",
    "    # Step-size\n",
    "    hh = Inf\n",
    "    for i in eachindex(vec1T)\n",
    "        @inbounds h1 = stepsize( vec1T[i], _epsAbs )\n",
    "        hh = min( hh, h1 )\n",
    "    end\n",
    "    \n",
    "    # Values at t0+h\n",
    "    @simd for i in eachindex(vec0)\n",
    "        @inbounds vec0[i] = evaluate( vec1T[i], hh )\n",
    "    end\n",
    "    \n",
    "    return hh, vec0\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "stepsize (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Returns the maximum step size from epsilon and the last two coefficients of the x-Taylor series \n",
    "function stepsize{T<:Real}(x::Taylor1{T}, epsilon::Float64)\n",
    "    ord = x.order\n",
    "    h = Inf\n",
    "    for k in [ord-1, ord]\n",
    "        kinv = 1.0/k\n",
    "        aux = abs( x[k+1] )\n",
    "        h = min(h, (epsilon/aux)^kinv)\n",
    "    end\n",
    "    return h\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II. The Kepler problem\n",
    "\n",
    "-------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a concrete example, we numerically integrate the cartesian equations of motion of the (planar) [Kepler problem](https://en.wikipedia.org/wiki/Kepler_problem):\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\dot{x} &=& v_x,\\\\\n",
    "\\dot{y} &=& v_y,\\\\\n",
    "\\dot{v}_x &=& - \\frac{G M x}{(x^2 + y^2)^{3/2}},\\\\\n",
    "\\dot{v}_y &=& - \\frac{G M y}{(x^2 + y^2)^{3/2}}.\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For concreteness, we fix $\\mu = GM = 1$, and Kepler's third law defines the units of time in terms of those of distance: $T= 2\\pi a^{3/2}$. The origin is the center of mass of the two bodies, so $x$ and $y$ above are actually relative coordinates to the center of mass. We choose the $x$-axis to be parallel to the major axis of the ellipse.\n",
    "\n",
    "The initial conditions for the particle are set at periapse, which we locate on the positive x-axis. Using the semimajor axis $a$ and the eccentricity $e$, we have\n",
    "$$ \n",
    "x_0 = a (1-e),\\\\\n",
    "y_0 = 0,\\\\\n",
    "v_{x_0} = 0,\\\\\n",
    "v_{y_0} = \\frac{l_z}{x_0} = m \\frac{\\sqrt{\\mu a (1-e^2)}}{x_0},\n",
    "$$\n",
    "where $l_z$ is the angular momentum. Below we set the mass $m$ equal to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " mass = 1.0\n",
      " a = 1.0\n",
      " e = 0.8\n",
      "\n"
     ]
    }
   ],
   "source": [
    "const mu = GM = 1.0\n",
    "\n",
    "const masa = 1.0\n",
    "const semieje = 1.0\n",
    "const excentricidad = 0.8\n",
    "\n",
    "println(\" mass = $masa\\n a = $semieje\\n e = $excentricidad\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following functions allow us to calculate the energy and angular momentum using cartesian coordinates or the semimajor axis and eccentricity of the orbit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "energy (generic function with 2 methods)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function energy{T<:Real}( x::T, y::T, vx::T, vy::T )\n",
    "    eneCin = 0.5*(vx*vx + vy*vy)\n",
    "    r = sqrt( x*x + y*y)\n",
    "    enePot = -GM*masa/r\n",
    "    return eneCin + enePot\n",
    "end\n",
    "energy{T<:Real}(a::T) = - 0.5*GM*masa/a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lz1 (generic function with 1 method)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lz{T<:Real}( a::T, e::T ) = masa * sqrt( GM*a*(1-e^2) )\n",
    "lz1{T<:Real}( x::T, y::T, vx::T, vy::T ) = masa*( x*vy - y*vx )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned above, we set the initial conditions in cartesian coordinates from the values of the semimajor axis and eccentricity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "iniCond (generic function with 1 method)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function iniCond{T<:Real}(a::T, e::T)\n",
    "    x0  = a*(1-e)\n",
    "    vy0 = lz(a, e) / x0\n",
    "    y0  = zero(vy0)\n",
    "    vx0 = zero(vy0)\n",
    "    return x0, y0, vx0, vy0\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III. Taylor integration of the Kepler problem\n",
    "\n",
    "-------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned above, Taylor's integration method is quite powerful. Yet, the equations of motion have to be implemented individually, and this involves a bit more than simply defining a function which contains the equations of motion, although the latter can actually be implemented just like that, which is paid back in performance.\n",
    "\n",
    "Below, the functions `jetKepler1` and `jetKepler2` are two implementations of the equations of motion for the Kepler problem. Both return a vector that contains the Taylor series for the dynamical variables, taking as input a vector of Taylor coefficients.\n",
    "\n",
    "The former function is an *almost* straight forward implementation of the equations of motion, except for the fact that we have to use `Taylor`-type variables. The latter yields the same results, but the implementation is done by splitting the equations of motion in unary and binary (elementary) functions and operations. This second method is quite *hand-crafted* to the actual problem, and it is not the most explicit way of doing this; yet, it turns out to be much more efficient than the former."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "jetKepler1 (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function jetKepler1{T<:Real}( vec::Array{T,1} )\n",
    "\n",
    "    xT = Taylor1(vec[1], _ordenTaylor)\n",
    "    yT = Taylor1(vec[2], _ordenTaylor)\n",
    "    vxT = Taylor1(vec[3], _ordenTaylor)\n",
    "    vyT = Taylor1(vec[4], _ordenTaylor)\n",
    "\n",
    "\n",
    "    for k = 0:_ordenTaylor-1\n",
    "        knext = k+1\n",
    "        # Taylor expansions up to order k\n",
    "        # This part makes it somewhat slower the implementations, since there\n",
    "        # are many operations which are completely superfluous\n",
    "        xTt = Taylor1( xT[1:k+1], k)\n",
    "        yTt = Taylor1( yT[1:k+1], k)\n",
    "        vxTt = Taylor1( vxT[1:k+1], k)\n",
    "        vyTt = Taylor1( vyT[1:k+1], k)\n",
    "        # Eqs of motion <--- This as straight forward as it can be\n",
    "        xDot = vxTt\n",
    "        yDot = vyTt\n",
    "        rrt = ( xTt^2 + yTt^2 )^(3/2)\n",
    "        vxDot = -GM * xTt / rrt\n",
    "        vyDot = -GM * yTt / rrt\n",
    "        # The equations of motion define the recurrencies\n",
    "        xT[knext+1]  = xDot[knext] / knext\n",
    "        yT[knext+1]  = yDot[knext] / knext\n",
    "        vxT[knext+1] = vxDot[knext] / knext\n",
    "        vyT[knext+1] = vyDot[knext] / knext\n",
    "    end\n",
    "    \n",
    "    return Taylor1[ xT, yT, vxT, vyT ]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "jetKepler2 (generic function with 1 method)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function jetKepler2{T<:Real}( vec::Array{T,1} )\n",
    "\n",
    "    xT = Taylor1(vec[1], _ordenTaylor)\n",
    "    yT = Taylor1(vec[2], _ordenTaylor)\n",
    "    vxT = Taylor1(vec[3], _ordenTaylor)\n",
    "    vyT = Taylor1(vec[4], _ordenTaylor)\n",
    "\n",
    "    # Auxiliary quantities\n",
    "    x2T = zeros( T, _ordenTaylor+1 )\n",
    "    y2T = zeros( T, _ordenTaylor+1 )\n",
    "    sT  = zeros( T, _ordenTaylor+1 )\n",
    "    rT3 = zeros( T, _ordenTaylor+1 )\n",
    "    Fx  = zeros( T, _ordenTaylor+1 )\n",
    "    Fy  = zeros( T, _ordenTaylor+1 )\n",
    "\n",
    "    # Now the implementation\n",
    "    for k = 0:_ordenTaylor-1\n",
    "        knext = k+1\n",
    "        # The right-hand size of the eqs of motion\n",
    "        # This is more adpated for this problem, and avoids many superfluous operations\n",
    "        x2T[knext] = TaylorSeries.squareHomogCoef(k, xT.coeffs)\n",
    "        y2T[knext] = TaylorSeries.squareHomogCoef(k, yT.coeffs)\n",
    "        sT[knext]  = x2T[knext] + y2T[knext]\n",
    "        rT3[knext] = TaylorSeries.powHomogCoef(k, sT, 1.5, rT3, 0)\n",
    "        Fx[knext]  = TaylorSeries.divHomogCoef(k, xT.coeffs, rT3, Fx, 0)\n",
    "        Fy[knext]  = TaylorSeries.divHomogCoef(k, yT.coeffs, rT3, Fy, 0)\n",
    "        # The equations of motion define the recurrencies\n",
    "        xT[knext+1]  = vxT[knext] / knext\n",
    "        yT[knext+1]  = vyT[knext] / knext\n",
    "        vxT[knext+1] = -GM * Fx[knext] / knext\n",
    "        vyT[knext+1] = -GM * Fy[knext] / knext\n",
    "    end\n",
    "    \n",
    "    return Taylor1[ xT, yT, vxT, vyT ]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following shows some benchmarks for 10 identic one-step integrations, using both implementations of the equations of motion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.017379273627668643,[0.19626418116550612,0.05181472066492753,-0.42543199148800787,2.944787769051677])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0, y0, vx0, vy0 = iniCond(semieje, excentricidad)\n",
    "\n",
    "taylorStepper( jetKepler1, [x0, y0, vx0, vy0] );\n",
    "\n",
    "timeJK1 = @elapsed begin\n",
    "    for i=1:10\n",
    "        taylorStepper( jetKepler1, [x0, y0, vx0, vy0] );\n",
    "    end\n",
    "end\n",
    "taylorStepper( jetKepler1, [x0, y0, vx0, vy0] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.017379273627668643,[0.19626418116550612,0.05181472066492753,-0.42543199148800787,2.944787769051677])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0, y0, vx0, vy0 = iniCond(semieje, excentricidad)\n",
    "\n",
    "taylorStepper( jetKepler2, [x0, y0, vx0, vy0] );\n",
    "\n",
    "timeJK2 = @elapsed begin\n",
    "    for i=1:10\n",
    "        taylorStepper( jetKepler2, [x0, y0, vx0, vy0] );\n",
    "    end\n",
    "end\n",
    "taylorStepper( jetKepler2, [x0, y0, vx0, vy0] )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "timeJK1 = 0.001156201   timeJK2 = 0.000211803 \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.458850913348725"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println( \"timeJK1 = $(timeJK1)   timeJK2 = $(timeJK2) \")\n",
    "tau = timeJK1/timeJK2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are identic, as it should be; yet, the elapsed time is somewhat shorter when using `jetKepler2`. The numbers clearly show that it is worth taking the time to construct the jet as in `jetKepler2`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we carry out a long integration of this rather eccentric keplerian orbit (eccentricity is 0.8). Everything needed is included in the function `keplerIntegration`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "keplerIntegration (generic function with 1 method)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function keplerIntegration( a::Float64, e::Float64, time_max::Float64, jetEqs::Function )\n",
    "    # Initial conditions, energy and angular momentum\n",
    "    t0 = 0.0\n",
    "    x0, y0, vx0, vy0 = iniCond(a, e)\n",
    "    ene0 = energy(x0, y0, vx0, vy0)\n",
    "    lz0 = lz1(x0, y0, vx0, vy0)\n",
    "    \n",
    "    # Change, measured in the local `eps` of the change of energy and angular momentum\n",
    "    eps_ene = eps(ene0); dEne = zero(Int)\n",
    "    eps_lz = eps(lz0); dLz = zero(Int)\n",
    "    \n",
    "    # Vectors to plot the orbit with PyPlot\n",
    "    tV, xV, yV, vxV, vyV = Float64[], Float64[], Float64[], Float64[], Float64[]\n",
    "    DeneV, DlzV = Int[], Int[]\n",
    "    push!(tV, t0)\n",
    "    push!(xV, x0)\n",
    "    push!(yV, y0)\n",
    "    push!(vxV, vx0)\n",
    "    push!(vyV, vy0)\n",
    "    push!(DeneV, zero(Int))\n",
    "    push!(DlzV, zero(Int))\n",
    "    \n",
    "    # This is the main loop; we include a minimum step size for security\n",
    "    dt = 1.0\n",
    "    while t0 < time_max && dt>1.0e-8\n",
    "        # Here we integrate\n",
    "        dt, (x1, y1, vx1, vy1) = taylorStepper( jetEqs, [x0, y0, vx0, vy0] );\n",
    "        t0 += dt\n",
    "        push!(tV,t0)\n",
    "        push!(xV,x1)\n",
    "        push!(yV,y1)\n",
    "        push!(vxV, vx1)\n",
    "        push!(vyV, vy1)\n",
    "        eneEnd = energy(x1, y1, vx1, vy1)\n",
    "        lzEnd = lz1(x1, y1, vx1, vy1)\n",
    "        dEne = trunc( Int, (eneEnd-ene0)/eps_ene )\n",
    "        dLz  = trunc( Int, (lzEnd-lz0)/eps_lz )\n",
    "        push!(DeneV, dEne)\n",
    "        push!(DlzV, dLz)\n",
    "        x0, y0, vx0, vy0 = x1, y1, vx1, vy1\n",
    "    end\n",
    "\n",
    "    return tV, xV, yV, DeneV, DlzV\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  4.623210 seconds (47.67 M allocations: 3.549 GB, 6.94% gc time)\n"
     ]
    }
   ],
   "source": [
    "#jetKepler1\n",
    "tV1, xV1, yV1, DeneV1, DlzV1 = keplerIntegration( semieje, excentricidad, 2pi, jetKepler1);\n",
    "@time tV1, xV1, yV1, DeneV1, DlzV1 = \n",
    "    keplerIntegration( semieje, excentricidad, 1000*2pi, jetKepler1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.882801 seconds (16.18 M allocations: 494.170 MB, 4.83% gc time)\n"
     ]
    }
   ],
   "source": [
    "#jetKepler2\n",
    "tV2, xV2, yV2, DeneV2, DlzV2 = keplerIntegration( semieje, excentricidad, 2pi, jetKepler2);\n",
    "@time tV2, xV2, yV2, DeneV2, DlzV2 = \n",
    "    keplerIntegration( semieje, excentricidad, 1000*2pi, jetKepler2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking the consistency of the results after 1000 periods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(true,true)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tV1[end] == tV2[end], yV1[end] == yV2[end]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The minimum value of the step-size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.017379273627668643"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minimum([tV1[i+1]-tV1[i] for i=1:length(tV1)-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which, in units of the orbital period, is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0027659973051900803"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ans/(2pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The average step-size is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.14405720621814588"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(tV1[end]-tV1[1])/(length(tV1)-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02292741645762644"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ans/(2pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us plot the trajectory using `PyPlot`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAroAAAILCAYAAAAHaz/JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X2MVfWd+PHPHaeZGR4cEXwYtClFQNuKQayztHZTSpN2JVmruzjuhIesS/CPTdvsbF3sg2ZbU5NtNkqz7U4amqZsy4bVBbvtrlGzxFAgBQsC0Z+xmIi6RkYrC0scGEDL+f3RMGWYpzsz9+l87+uVEOXce+6c4cz33vf93u+9U8iyLAsAAEhMQ7UPAAAAykHoAgCQJKELAECShC4AAEkSugAAJEnoAgCQJKELAECSGqt9AOV05MiRePrpp2PWrFnR0tJS7cMBAOACfX198dprr8XnP//5mDFjRklvO+nQffrpp2PFihXVPgwAAEaxcePGWL58eUlvM+nQnTVrVkT8/h/uIx/5SHUPhoro6uqKdevWVfswqBDnu7443/XF+a4fL730UqxYsaK/20op6dA9t1zhIx/5SCxcuLDKR0MltLa2Otd1xPmuL853fXG+6085lpl6MxoAAEkSugAAJEnoAgCQJKFLUjo7O6t9CFSQ811fnO/64nxTCkKXpLhjrC/Od31xvuuL800pCF0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSUIXAIAkCV0AAJIkdAEASJLQBQAgSY3VPgCAWlQoVPsIyiPLqn0EAJUjdIHkjTdaU4vCQmHs/xap/RsA9UXoArlWbLgJtrH/G4wljP37ArVI6AI1rZjQElnlUey/62hB7PwA1SJ0gaobLWaFUm0b6fyIYKCahC5QUcNFj+BJ03gj2M8DUApCFygLAcNohvtZEMBAqQhdYMJECaU01gD2cwYMR+gCYyJqqZahfsbELzASoQuMSERQy4qNXz+zUJ+ELtDPbC0puPDnVfhC/RK6UMfM1lIPhC/UL6ELdUTYgvCFeiJ0IXEewGFkwhfSJXQhMWZtYWJGC1/jCfJD6EICzD5B+Zw/nsz2Qr4IXcgpD7ZQeWZ7IV+ELuSIuIXaMtJsr/EJ1ddQqS905syZuO++++Lqq6+OSZMmxaJFi2Lr1q2j7vetb30rGhoaBv2ZNGlSBY4aquvcA+e5P1k28A9QOy4cm+ePXaA6Kjaju2rVqvjZz34WXV1dMWfOnNiwYUMsXbo0tm3bFp/85CdH3LdQKMQPfvCDmDx5cv+2iy66qNyHDFVjVgjyzUwv1IaKhO6vf/3reOyxx+Lhhx+Orq6uiIhYuXJlXH/99bF27drYuXPnqLfx53/+53HppZeW+1ChajwQQppEL1RPRZYubN68ORobG2PNmjX925qammL16tWxa9euePPNN0e9jbNnz8a7775bzsOEirvwpU1LEiBtljdAZVUkdA8cOBDz5s2LKVOmDNje3t7ef/lIsiyL2bNnR2tra0ydOjVWrlwZv/3tb8t2vFBOw8WtwIX6Mlz0AqVTkaULPT090dbWNmh7W1tbZFkWhw8fHnbfadOmxZe+9KX4xCc+EU1NTbFjx474/ve/H3v27Im9e/cOimeoVV6uBIYzXOy6r4CJqUjo9vX1RVNT06Dtzc3N/ZcP58tf/vKAv99xxx1x8803x/Lly6O7uzvWrl1b2oOFEvJxYMBYDLee130HjE9FQrelpSVOnz49aPupU6f6Lx+Lzs7O+MpXvhJbt24tKnS7urqitbV10G10dnaO6etCsczIABN1YfQOtR3yZtOmTbFp06YB244fP162r1eR0G1raxtyeUJPT09ERMycOXPMt/nBD34wjh49WtR1161bFwsXLhzz14Cx8EAElMtQSxvcz5BHQ0007tu3L2666aayfL2KvBltwYIF8fLLL0dvb++A7bt3745CoRALFiwY822+9tprcdlll5XqEGHcfGoCUCk+sQHGpiKhu2zZsnj//fdj/fr1/dvOnDkTGzZsiEWLFsVVV10VERFvvPFGHDx4cMC+R44cGXR73d3d8c4778Stt95a3gOHEQhcoFp8YgMUpyJLF9rb2+POO++Mr33ta/H222/3/2a0119/PX784x/3X2/lypWxffv2OHv2bP+2D33oQ3HXXXfF/Pnzo7m5OXbs2BGPPvpoLFy4MO65555KHD70szwBqDWWNcDwKvYrgH/605/GAw88EBs3boxjx47FDTfcEE888UTccsst/dcpFArR0DBwknnFihXxq1/9Kh5//PE4depUfOhDH4qvfvWr8fWvf73/Uxug3AQuUOsELwxWyLJ0h8G5xc3PPfecN6MxLgIXyCv3X+RFOXutYjO6kCceIIC88/FkIHRhAA8GQIosa6BeCV0IgQvUB8FLvRG61DWBC9QjwUu9ELrUJYELIHhJn9ClrghcgMEEL6kSutQFgQswuguD1/0leSd0SZrABRi7LDO7SxqELslyBw0wfpYzkAKhS3LM4gKUjuAlz4QuyRC4AOVj/S55JHTJPYELUDnW75InQpdcc0cLUHmWM5AXDdU+ABgPd64A1ZdlA6MXao3QJXfOD1yRC1B9Q83wQi0QuuSGWVyA2mV2l1okdMkFs7gA+WB2l1oidKlpZnEB8sfsLrXCpy5QswQuQL5dGLvuz6k0M7rUHLO4AGkxu0u1CF1qirW4AGmydpdqELrUBLO4AOmzdpdKE7pUnVlcgPpidpdKEbpUlVlcgPpkdpdKELpUhaUKAESIXcpL6FJxlioAcD5LGSgXoUtFmcUFYCiWMlAOQpeKsFQBgGKIXUpJ6FJ2lioAMBZil1IRupSVWVwAxsO6XUpB6FI2IheAibBul4kSupSFyAWgVMQu4yV0KTmRC0CpiV3GQ+hSMj5ZAYByEruMldClJHyyAgCVIHYZC6HLhJnFBaCSxC7FErpMiMgFoBrELsUQuoybyAWgmsQuo2ms9gGQTyIXgFpwYex6XOJ8ZnQZM3cmANQas7sMRegyJiIXgFoldrmQ0KVoIheAWid2OZ/QpSgiF4C88FjFOUKXUYlcAPImy8zqInQZhcgFIM/Ebn0TugxL5AKQZ9brInQZksgFIAVit74JXYYlcgFIgditX0KXQQoFkQtAWsRufRK6DOAOAIBUid36I3TpZ10uAKkTu/VF6BIRIheA+uGxrn4IXUQuAHXJrG76hC4RIXIBqC+WMNQHoVvnDHAA6pXYTZ/QrWOWLABQ7zwGpk3o1imRCwC/l2VmdVMldOuYyAWAPxC76RG6dchABoCBrNdNk9CtM5YsAMDQPDamR+jWEZELAKMzq5sOoVtnRC4ADM8ShrQI3TphwAJAcUwKpUPo1hEDFwCKZ5Io/4RuHSgURC4AjIUlDGkQuokzQAFgfEwS5Z/QrQMGKgCMn0mj/BK6CTMwAWBiTBblm9BNnAEKABNn8iifhG6ivAENAErDG9PyS+gmyEAEgNIyeZRPQjdRBiQAlJ7JpHwRuokxAAGgPEwi5Y/QTZCBCADlY1IpP4RuQgw8ACgvk0n5InQTYwACQPmZXMoHoZsIAw4AKsOkUn4I3YQYeABQOSaZap/QTYCBBgCVZXIpH4RuIgw4AICBhC4AwDhkmVdVa53QzTkDDABgaEI3AZYtAED1mHSqXUI3xwwsAKguk021TejmnAEGADA0oQsAMEFeZa1NQjenDCgAqA1eXa1dQjfHDCwAqB0moWqP0AUAmCCTT7VJ6OaQZ4wAAKMTujnlmSMAwMiELgBAiXjVtbYIXQCAEvBqa+0RujnjmSIAQHGEbg55xggAtcukVO0QugAAJWIyqrYIXQAAkiR0c8RLIQAAxRO6OeMlEQCA4ghdAIAS8ypsbRC6AAAl5NXX2iF0AQBIktDNCS+BAACMjdDNES+FAAAUT+gCAJAkoQsAUAaWHVaf0AUAKDHLDWuD0AUAIElCFwCAJAldAACSVLHQPXPmTNx3331x9dVXx6RJk2LRokWxdevWovY9fPhwdHR0xLRp06K1tTVuv/32ePXVV8t8xLXDYnYAgLGrWOiuWrUqvvvd78aKFSvin/7pn6KxsTGWLl0av/rVr0bc78SJE7F48eLYsWNH3H///fHggw/G/v37Y/HixXHs2LEKHX31WdQOADA2jZX4Ir/+9a/jsccei4cffji6uroiImLlypVx/fXXx9q1a2Pnzp3D7vvP//zP8corr8SePXti4cKFERHxJ3/yJ3H99dfHww8/HN/+9rcr8S0AAJAzFZnR3bx5czQ2NsaaNWv6tzU1NcXq1atj165d8eabbw6775YtW+Lmm2/uj9yIiGuvvTY++9nPxmOPPVbW4wYAqmPTC5uqfQgkoCKhe+DAgZg3b15MmTJlwPb29vb+y4eSZVk8//zz8fGPf3zQZe3t7fHKK6/EiRMnSn/AAEBVbfp/QpeJq0jo9vT0RFtb26DtbW1tkWVZHD58eMj9jh49GqdPnx5234gYdl8AAOpbRUK3r68vmpqaBm1vbm7uv3y4/SJiXPsCAFDfKvJmtJaWljh9+vSg7adOneq/fLj9ImJc+56vq6srWltbB2zr7OyMzs7OUfcFAMpv0wubBixX+M+X/zNu23Rb/987r++Mzvket/Nu06ZNsWnTwGUpx48fL9vXq0jotrW1DbnEoKenJyIiZs6cOeR+l156aTQ1NfVfb6h9h1rWcKF169YNeDMbAFBbOucPDNnbNt0Wv+j8RRWPiHIYaqJx3759cdNNN5Xl61Vk6cKCBQvi5Zdfjt7e3gHbd+/eHYVCIRYsWDDkfoVCIebPnx979+4ddNmzzz4bs2fPHvQGNwAAiKhQ6C5btizef//9WL9+ff+2M2fOxIYNG2LRokVx1VVXRUTEG2+8EQcPHhy07549e2Lfvn392w4ePBjPPPNMdHR0VOLwAQDIoYosXWhvb48777wzvva1r8Xbb78dc+bMiQ0bNsTrr78eP/7xj/uvt3Llyti+fXucPXu2f9tf//Vfxw9/+MNYunRp3HvvvdHY2Bjr1q2Ltra2+Nu//dtKHD4AUGGd11uPy8RVJHQjIn7605/GAw88EBs3boxjx47FDTfcEE888UTccsst/dcpFArR0DBwknnKlCnxy1/+Mrq6uuKhhx6Ks2fPxmc+85l45JFHYvr06ZU6fACggrzxjFIoZFmWVfsgyuXc4ubnnnsu129GKxR+/990zxQApKdQ8NhdjHL2WkXW6DIxBgkAwNgJXQAAkiR0AQBIktAFACBJQhcAoMTOvZGc6hK6AABl4M3k1Sd0AQBIktDNES+DAAAUT+jmhJc/AADGRugCAJAkoQsAUEKWGtYOoQsAUGKWHNYGoZszniUCABRH6OaIZ4cAAMUTugAAJEnoAgCUiCWGtUXo5pBBBAC1y1LD2iF0c8bgAQAojtAFACgBr7jWHqELAFAiXnmtLUI3pzxrBAAYmdDNIc8WAQBGJ3QBACbIK621SejmmEEFALXDK661R+jmlMEEALXBxFPtEroAABNkAqo2Cd2c8ywSAGBoQjfHPHsEgOoy4VTbhG4CDDIAqB4TT7VL6OacwQUA1WGiqfYJ3UQYbABQeSacapvQTYBBBgAwmNAFABgjr6Tmg9BNiEEHAJXjFdXaJ3QTYbABQGWYWMoPoZsYgw8Ays8EUz4I3YQYdABQXiaU8kXoJsggBIDyMbGUH0I3MQYfAJSHiaT8EbqJMhgBoPRMKOWL0E3QuUEodgGgNDym5pPQTZRnnABQGuci12Nr/gjdxHkGCgATJ3LzSegmzKAEgIkxYZRvQrcOGKQAMH4mjvJL6CbOG9MAYHw8duaf0K0DnokCwNh4A1oahG4d8cwUAIoncvNP6NYJgxUAimNiKB1Ct45kmcELACOxZCEtQrcOiV0AGJ7ITYfQrTM+hQEAhuaxMT1Ctw55pgoAA1mykCahW8c8cwUAkZsyoVunLGEAgD8QuWkSunVM7AJQ7zwGpk3o1jnPYAGoV5YspE/oEhGe0QJQX0RufRC6WMIAQF0SuekTukSE2AWgfhQKIrdeCF36iV0AUucxrr4IXQYQuwCkyrrc+iN0GcQdAACpEbn1SegypCwzqwtAGkRu/RK6jEjsApBnIre+CV2GZb0uAHkmchG6jEjsApBHIpcIoUsRxC4AeSJyOUfoUhSxC0AeiFzOJ3QpmtgFoJaJXC4kdBkTsQtALRK5DEXoMmZiF4BaInIZjtBlXMQuALVA5DISocu4iV0AqknkMhqhy4SIXQCqQeRSDKHLhJ0fu4IXgHITuRRL6FISWWZ2F4DyE7mMhdClpMQuAOUichkroUvJiV0ASk3kMh5Cl7IQuwCUwvnv/xC5jJXQpWzELgATcX7gilzGQ+hSVmIXgPEwi0spNFb7AEjfhbHrTguAkXi8oFTM6FIxZncBGIn1uJSa0KWixC4AQ7Eel3IQulSc36QGwPnM4lIuQpeq8JvUALBUgXITulSV2AWoT5YqUAlCl6qzlAGgvpjFpVKELjXBUgaA9FmqQKUJXWqK2V2ANFmqQDUIXWqO2V2AdJjFpZqELjXL7C5AvpnFpdqELjXN7C5A/pjFpVYIXXLB7C5APpjFpZYIXXLD7C5A7TKLSy0SuuSO2V2A2mIWl1rVWO0DgPEYKnbduQJU1vmTDe6DqUVCl1y7MHjd0QKUn8AlLyxdIAmWMwBUhmUK5IkZXZJhOQNA+bhfJY/M6JIcn84AUDomD8gzM7ok68LYdQcNUDzrcEmB0CV5WWZGAqBYApeUCF3qgvW7ACMTuKRI6FJXBC/AYO4PSZXQpS4JXgCzuKRP6FLXBC9QjwQu9ULoQgheIH0Xftyi+zjqgdCF8wheIDVmb6lnFfuFEcePH4977rknLr/88pgyZUosWbIk9u/fX9S+d999dzQ0NAz689GPfrTMR029uvCXTvjFE0DeXPiEXeRSjyoyo5tlWSxdujReeOGFWLt2bUyfPj26u7tj8eLFsW/fvrjmmmtGvY3m5ub40Y9+FNl5I7W1tbWchw1meIHcMYMLf1CR0P33f//32LVrV2zZsiXuuOOOiIi48847Y968efH3f//3sXHjxlFvo7GxMTo7O8t9qDCk8x8sPIgAtch9EwxWkaULW7ZsiSuvvLI/ciMiZsyYER0dHfHzn/883nvvvaJuJ8uy6O3tLddhQlEsawBqiSUKMLyKhO7+/ftj4cKFg7a3t7fHyZMn4+WXXx71Nk6ePBlTp06Niy++OKZPnx5f/OIX48SJE+U4XCiK4AWq5dx9TqHwh/sigQuDVWTpQk9PT3z6058etL2trS0iIg4fPhwf+9jHht1/5syZsXbt2li4cGGcPXs2nnrqqeju7o7nn38+tm3bFg0NFXtPHQwyXOx60AFKzX0MjM2YQzfLsjhz5kxR121qaoqIiL6+vv7/P19zc3NkWRZ9fX0j3s5DDz004O8dHR0xd+7cuP/++2Pz5s3R0dFR5NFD+Vy4jteb14BS8Pm3MH5jDt3t27fHZz7zmVGvVygU4qWXXop58+ZFS0tLnD59etB1Tp06FYVCIVpaWsZ6GNHV1RUPPPBAbN26ddTQ7erqGvQJDZ2dnd7cRtmY5QUmyn0HKdq0aVNs2rRpwLbjx4+X7euNOXSvu+662LBhQ1HXPbc0oa2tLXp6egZdfm7bzJkzx3oY0dzcHNOnT4+jR4+Oet1169YNuUYYys0sLzAWZm9J3VATjfv27YubbrqpLF9vzKF7xRVXxKpVq8a0z4IFC2Lnzp2Dtu/evTsmTZoU8+bNG+thRG9vbxw5ciQuu+yyMe8L1eAjyoChiFson4q8i2vZsmXx9ttvx+OPP96/7ciRI7F58+a47bbb4gMf+ED/9kOHDsWhQ4f6/3769OkhP1LswQcfjIiIW2+9tYxHDuUx1Cc2+NQGqC9DfSyYyIXSqsinLixbtiy++93vxt133x0vvvhizJgxI7q7u+N3v/tdfPOb3xxw3SVLlkRDQ0N/7L711ltx4403RmdnZ1x33XUREfHUU0/Fk08+GUuXLo3bbrutEt8ClMVwSxsuvAxIg9lbqKyKhG5DQ0M8+eST8Xd/93fxve99L/r6+qK9vT1+8pOfxNy5cwdct1AoROG8e4JLLrkk/vRP/zS2bt0aP/nJT+J3v/tdzJkzJ/7hH/4hvvKVr1Ti8KEiRC+kSdxC9VQkdCMiWltbY/369bF+/foRr/fqq68O2u9f/uVfynloUHNEL+SbuIXaULHQBcZH9EI+iFuoPUIXckT0Qu0Y6g2kxiHUFqELOTVS9F54OVAaxhnki9CFBFz4YGu2F0rDrC3km9CFBJnthfEzXiAdQhcSN9ps71DXgXpi1hbSJXShzghf6p2whfohdKHOCV9SJ2yhfgldYADhS56JWuB8QhcY0VCRICaoBUP9HEb4WQT+QOgCY1bMrO9Q14PxErXAeAhdYMKGiw1xwnj4uQFKRegCZTPcsgchw3A/A+f4WQBKQegCFTXS7O9I8SN88sk5BapJ6AI1YaToGS2CR9uf8nFegFomdIGaN1osFRPCxdwOf1DMv+c5/l2BWiV0gdwrJrSKjeHx3HatG8/3HZHG9w7UN6EL1IXxRtt4I7HWiFagHgldgBEIRID8aqj2AQAAQDkIXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQhcAgCQJXQAAkiR0AQBIktAFACBJQpekbNq0qdqHQAU53/XF+a4vzjelIHRJijvG+uJ81xfnu74435SC0AUAIElCFwCAJAldAACS1FjtAyinvr6+iIh46aWXqnwkVMrx48dj37591T4MKsT5ri/Od31xvuvHuU47122lVMiyLCv5rdaIf/3Xf40VK1ZU+zAAABjFxo0bY/ny5SW9zaRD98iRI/H000/HrFmzoqWlpdqHAwDABfr6+uK1116Lz3/+8zFjxoyS3nbSoQsAQP3yZjQAAJIkdAEASJLQBQAgSUIXAIAkCV0AAJKUTOg+88wzsXr16rj22mtj8uTJcc0118SaNWvirbfeKvo2Dh8+HB0dHTFt2rRobW2N22+/PV599dUyHjXj9dZbb8VXv/rVWLJkSVx88cXR0NAQ27dvL3r/b33rW9HQ0DDoz6RJk8p41IzXRM93hPGdN8ePH4977rknLr/88pgyZUosWbIk9u/fX9S+d99995Dj+6Mf/WiZj5rRnDlzJu677764+uqrY9KkSbFo0aLYunVrUfsaw/kz3vNdysfoZH4z2n333RfHjh2LO++8M+bOnRuHDh2K733ve/HEE0/EgQMH4vLLLx9x/xMnTsTixYvj3Xffjfvvvz8aGxvjkUceicWLF8eBAwdi2rRpFfpOKMbBgwfjH//xH2Pu3Llxww03xK5du8Z8G4VCIX7wgx/E5MmT+7dddNFFpTxMSmSi59v4zpcsy2Lp0qXxwgsvxNq1a2P69OnR3d0dixcvjn379sU111wz6m00NzfHj370ozj/EzRbW1vLedgUYdWqVfGzn/0surq6Ys6cObFhw4ZYunRpbNu2LT75yU8Ou58xnE/jPd8RJXyMzhKxY8eOQdu2b9+eFQqF7IEHHhh1/+985ztZQ0ND9txzz/Vv+81vfpM1NjZm3/jGN0p6rExcb29vduzYsSzLsmzz5s1ZQ0ND9stf/rLo/b/5zW9mDQ0N2f/+7/+W6xApoYmeb+M7Xx599NGsUChkjz/+eP+2d955J5s2bVq2fPnyUff/y7/8y2zq1KnlPETG4dlnn80KhUL2yCOP9G87depUNmfOnOyWW24ZcV9jOH8mcr5L+RidzNKFT33qU4O2/fEf/3Fceuml/b9DeSRbtmyJm2++ORYuXNi/7dprr43Pfvaz8dhjj5X0WJm4yZMnxyWXXDLh2zl79my8++67JTgiymmi59v4zpctW7bElVdeGXfccUf/thkzZkRHR0f8/Oc/j/fee6+o28myLHp7e8t1mIzR5s2bo7GxMdasWdO/rampKVavXh27du2KN998c9h9jeH8mcj5PqcUj9HJhO5QTpw4Eb29vaP+Orksy+L555+Pj3/844Mua29vj1deeSVOnDhRrsOkSrLrcS2PAAAGWElEQVQsi9mzZ0dra2tMnTo1Vq5cGb/97W+rfViUmPGdP/v37x8QNOe0t7fHyZMn4+WXXx71Nk6ePBlTp06Niy++OKZPnx5f/OIXnecqO3DgQMybNy+mTJkyYHt7e3v/5UMxhvNpvOf7nFI9RiezRnco69ati/feey/+4i/+YsTrHT16NE6fPh1tbW2DLju37fDhwzF37tyyHCeVN23atPjSl74Un/jEJ6KpqSl27NgR3//+92PPnj2xd+/eQQOT/DK+86enpyc+/elPD9p+/vn62Mc+Nuz+M2fOjLVr18bChQvj7Nmz8dRTT0V3d3c8//zzsW3btmhoSHqOp2b19PQMOw6zLIvDhw8PuZ8xnE/jPd8RpX2MrsnQzbIszpw5U9R1m5qahty+ffv2ePDBB+Ouu+4a8g7zfH19fcPeVnNz84DrUHqlON9j9eUvf3nA3++44464+eabY/ny5dHd3R1r164tyddhsEqfb+O7usZzvvv6+oY9X1mWjXq+HnrooQF/7+joiLlz58b9998fmzdvjo6OjiKPnlIa6byeu3y4/SKM4bwZ7/mOKO1jdE0+rd2+fXu0tLSM+mfSpElDvoT1m9/8Jv7sz/4sbrjhhvjhD3846tdraWmJiIjTp08PuuzUqVMDrkPpTfR8l0pnZ2dceeWVRX/UDeNT6fNtfFfXeM53S0vLsOerUCiM63x1dXVFoVAwvqtopPN67vLh9oswhvNmvOd7OON9jK7JGd3rrrsuNmzYUNR1L5wWf+ONN+Jzn/tcTJs2LZ544okBH0sxnEsvvTSampqip6dn0GXntg01/U5pTOR8l9oHP/jBOHr0aFm/Rr2r9Pk2vqtrPOe7ra1txPM1c+bMMR9Hc3NzTJ8+3fiuora2tiFfrh7tvBrD+TTe8z2S8TxG12ToXnHFFbFq1aox73f06NH43Oc+F++9915s27YtrrjiiqL2KxQKMX/+/Ni7d++gy5599tmYPXu2NZtlNN7zXQ6vvfbakG+CoXQqfb6N7+oaz/lesGBB7Ny5c9D23bt3x6RJk2LevHljPo7e3t44cuRIXHbZZWPel9JYsGBBbNu2LXp7eweMud27d0ehUIgFCxYMuZ8xnE/jPd8jGc9jdE0uXRiPkydPxq233ho9PT3x5JNPxuzZs4e97htvvBEHDx4csG3ZsmWxZ8+e2LdvX/+2gwcPxjPPPGM9V84Ndb6PHDky6Hrd3d3xzjvvxK233lqpQ6MMjO/8W7ZsWbz99tvx+OOP9287cuRIbN68OW677bb4wAc+0L/90KFDcejQof6/nz59esiPFHvwwQcjIozvKlq2bFm8//77sX79+v5tZ86ciQ0bNsSiRYviqquuighjOBUTOd+lfIwuZNl5vzYmx26//fb4xS9+EatXr47FixcPuGzKlCnxhS98of/vixcvju3bt8fZs2f7t/X29saNN94Y7777btx7773R2NgY69atiyzLYv/+/TF9+vRKfSsU6dvf/nYUCoV48cUX49/+7d/ir/7qr+LDH/5wRER84xvf6L/eUOd78uTJcdddd8X8+fOjubk5duzYEY8++mjceOONsXPnzv7F8tSOiZxv4ztfzp49G5/61KfixRdfjHvvvTdmzJgR3d3d8T//8z+xd+/eAe+unzVrVjQ0NPTH7uuvvx433nhjdHZ2xnXXXRcREU899VQ8+eSTsXTp0viv//qvqnxP/N5dd90V//Ef/xF/8zd/0/+bsvbu3RvPPPNM3HLLLRFhDKdkvOe7pI/RE/6VEzVi1qxZWUNDw5B/PvzhDw+47uLFi7OLLrpo0G28+eabWUdHR3bJJZdkF198cfaFL3whe+WVVyr1LTBGhUJhyPN94bkd6nzfc8892fXXX5+1trZmTU1N2bx587Kvf/3rWW9vbyW/BcZgIuc7y4zvvPm///u/bM2aNdlll12WTZkyJVuyZEm2b9++QdebNWtWNnv27AH7rVq1Kps3b142ZcqUrKWlJZs/f372ne98J3v//fcr+S0whNOnT2dr167NZs6cmbW0tGR/9Ed/lP33f//3gOsYw+kY7/ku5WN0MjO6AABwvmTW6AIAwPmELgAASRK6AAAkSegCAJAkoQsAQJKELgAASRK6AAAkSegCAJAkoQsAQJKELgAASRK6AAAkSegCAJCk/w8PM8PVWFgFuwAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x32d70b050>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "2-element Array{Any,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x32f8074d0>\n",
       " PyObject <matplotlib.lines.Line2D object at 0x32f807650>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "axis(\"equal\")\n",
    "plot(xV1, yV1, \",\", [0], [0], \"+\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we also plot the absolute change of energy and angular momentum as a function of time, in units of the \n",
    "local `eps` values for the initial quantities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAIUCAYAAAAHV9oiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xm8VXW9//H3RqbDgUwc8AA5kBD6QD0CEoFzGGqmlQziQHYRlR5dCctOWD5sEBPEMDNT1ELCMMIckuv1XuyiKRCXqWsGor/EBI4aDgh44ECu3x+4t2uvs+Zpr7X36/l4+JCz9ndNe1j7sz/r8/1+C4ZhGAIAAAAgSWpX6QMAAAAAsoQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADBJJUD+29/+pjFjxuiTn/yk6uvrdfDBB+vUU0/V448/3qbt+vXrddZZZ6lbt2468MADNX78eG3durVNO8MwNGPGDPXp00d1dXU6/vjj9eCDD6ZxOgAAAKhi7dPYyauvvqodO3bosssuU8+ePfX+++/roYce0nnnnafZs2fr8ssvlyRt3rxZJ598sg444ADdfPPN2r59u2655Rb99a9/1YoVK9S+/UeHO3XqVM2YMUNXXnmlBg8erEcffVQXXXSR2rVrpzFjxqRxWgAAAKhCBcMwjErs2DAMDRw4ULt379bf/vY3SdLXvvY1zZ07Vy+++KJ69eolSXrqqad05plnlgXSW7Zs0ZFHHqmrrrpKP/3pT0vbPPXUU7Vx40Zt3LhRhUIh/ZMCAABA7lWsBrlQKOgTn/iE3n333dKy3//+9zr33HNLwbEkffazn1W/fv20YMGC0rJHHnlEe/fu1aRJk8q2OWnSJG3atEnLli1L/gQAAABQlVINkN9//3299dZb+vvf/65Zs2bpiSee0IgRIyTtywq/+eabGjx4cJv1hgwZojVr1pT+Xrt2rerr69W/f/827QzDKGsLAAAABJFKDXLRN7/5Td19992SpHbt2umCCy7Qz372M0lSc3OzJKmhoaHNeg0NDXr77be1Z88edejQQc3NzerRo4dtO2lfsA0AAACEkWqAPGXKFI0ePVpbtmzRggUL9K9//Uu7d++WJLW0tEiSOnXq1Ga9zp07l9p06NBBLS0tnu2cbN26VU8++aSOOOII1dXVRT4nAAAAxKulpUUbN27UyJEjddBBB6W+/1QD5H79+qlfv36SpEsuuURnnXWWvvCFL+jPf/5zKVgtBsxmu3btkqRSm7q6Ol/t7Dz55JO65JJLop0IAAAAEjdv3jxdfPHFqe831QDZ6oILLtBVV12ll156qVQeUSy1MGtublb37t3VoUMHSftKKZYsWWLbTpJ69uzpuM8jjjhC0r4n/Oijj454BojblClTNGvWrEofBhzw+mQXr0128dpkG69PNq1bt06XXHJJKW5LW0UD5GIpxLZt29S3b18dfPDBWrlyZZt2K1asUGNjY+nvxsZG3XfffVq/fn1ZR73ly5erUCiUtbUqZpePPvpoDRw4MK5TQUz2339/XpcM4/XJLl6b7OK1yTZen2yrVDlsKqNY/POf/2yzbO/evbr//vtVV1enY445RtK+jPLjjz+uzZs3l9o99dRT2rBhQ9nkH+eff772228/3XnnnWXbvOuuu9SrVy8NGzYsoTMBAABAtUslg3zllVfqvffe0ymnnKJevXrp9ddf1wMPPKAXX3xRP/nJT9SlSxdJ0nXXXaeFCxfqtNNO0+TJk7V9+3bNnDlTxx9/vC677LLS9nr16qUpU6Zo5syZam1t1YknnqiHH35Yzz33nH7zm98wSQgAAABCSyVAvvDCC3Xffffprrvu0ltvvaVu3bpp0KBBuuWWW/T5z3++1K537956+umndc0112jq1Knq2LGjzj33XM2cObNUf1w0ffp0de/eXXfffbfuv/9+9e3bVw888IDGjh2bxikBAACgSqUSII8ZM6asRMLN0UcfrSeeeMJX26amJjU1NUU5NGTMuHHjKn0IcMHrk128NtnFa5NtvD6wUzAMw6j0QaRp9erVGjRokFatWkVRPgAAQAZVOl5LdappAAAAIOsIkAEAAAATAmQAAADAhAAZAAAAMCFABgAAAEwIkAEAAAATAmQAAADAhAAZAAAAMCFABgAAAEwIkAEAAAATAmQAAADAhAAZAAAAMCFABgAAAEwIkAEAAAATAmQAAADAhAAZAAAAMCFABgAAAEwIkAEAAAATAmQAAADAhAAZAAAAMCFABgAAAEwIkAEAAAATAmQAAADAhAAZAAAAMCFABgAAAEwIkAEAAAATAmQAAADAhAAZAAAAMCFABgAAAEwIkAEAAAATAmQAgRUKlT4CAACSQ4AMAAAAmBAgAwAAACYEyEBIeSszKB5v3o4bAIC0ESADIRlGpY8gmOLx5u24AQBIGwEyUCPMGeSksshkpwEA1YAAGagRcWaOnbZFdhoAUA0IkAEAAAATAmQgpLyVE9gdb9hzcCrTKC5PsowDAICkta/0AQBIl2F8FLwmURJBmQUAIO/IIAM1Iq7ANcngGgCALCBABmqEtfQhahmE17qUWAAA8ooSC6BGkQEGAMAeGWQgpLwHmGEzvG7nnffnBAAAiQAZCK0aSgjCnIPbOoxgAQCoBpRYACHlNVsadRQLr3WKjxMkAwDyigwyEFJeA0DzcUfJIDutSwYZAJB3BMhASHnOIEc59uK6dtswL8vr8wMAAAEyUAOSyOh6zcxHFhkAkFcEyEBIQQLAuILFODrVJTX+MQExAKBa0EkPCClICUFc5QZhO9UVCvbrhglqnY7B3PnPrR0AAFlHBhmoAUkEq4yHDACoVgTIgA92ZQpR1k+bXVlF2HINp1EszFnqSp8vAABRECADPlgzokEzpJXOqNqNXBH3MVX6HAEAiEsqAfLKlSv19a9/XQMGDFDXrl11+OGHa+zYsXrppZfatF2/fr3OOussdevWTQceeKDGjx+vrVu3tmlnGIZmzJihPn36qK6uTscff7wefPDBNE4HqBpJBbUEywCAPEulk9706dO1dOlSjR49Wscdd5xef/11/exnP9PAgQP15z//Wcccc4wkafPmzTr55JN1wAEH6Oabb9b27dt1yy236K9//atWrFih9u0/OtypU6dqxowZuvLKKzV48GA9+uijuuiii9SuXTuNGTMmjdNCDbF2cgs6E51TJ7mgxxBkn17rRZlNz7wNu7IKa4kFATMAIE8KhpH8V9fy5cs1ePDgsgD35Zdf1oABAzRmzBjNnTtXkvS1r31Nc+fO1YsvvqhevXpJkp566imdeeaZmj17ti6//HJJ0pYtW3TkkUfqqquu0k9/+tPSNk899VRt3LhRGzduVMGhCHL16tUaNGiQVq1apYEDByZ1yqgy1mAyjuAy6jFUSjEojiPoBwDATqXjtVRKLIYOHVoWHEvSUUcdpQEDBmjdunWlZb///e917rnnloJjSfrsZz+rfv36acGCBaVljzzyiPbu3atJkyaVbXPSpEnatGmTli1bltCZAOHE2WktTAfBpMZh9tou004DAPKoop303njjDR100EGS9mWF33zzTQ0ePLhNuyFDhmjNmjWlv9euXav6+nr179+/TTvDMMraAlkQ5zjIYToIJrF/P9uMc98AAKSlYgHyvHnztHnzZl144YWSpObmZklSQ0NDm7YNDQ16++23tWfPnlLbHj162LaT9gXbQBrSzo6GycgmdYx+tkv2GACQRxUJkNevX6+vf/3rGj58uMaPHy9JamlpkSR16tSpTfvOnTuXtWlpafHVDsiCYlAbddzhqMcQZLldG6fzcBsX2e8+AADIktSnmn7zzTf1+c9/XgcccIB+97vflTrT1dXVSZJ2797dZp1du3aVtamrq/PVzs2UKVO0//77ly0bN26cxo0bF+BsUCuilAlE6dCWdHmC3zIJu7Ze61qnngYAwM78+fM1f/78smXbtm2r0NHsk2qA/N5772nkyJF677339Oyzz+rQQw8tPVYsjyiWWpg1Nzere/fu6tChQ6ntkiVLbNtJUs+ePT2PZdasWYxigdDCBK5xjfqQ1ZrerB4XACDb7BKUxVEsKiW1Eovdu3frC1/4gl5++WUtWrRIn/rUp8oe79mzpw4++GCtXLmyzborVqxQY2Nj6e/Gxka9//77Wr9+fVm75cuXq1AolLUF4mAtLQhSNuE2PrCfdYI8lgS34/eaehoAgDxKJUD+4IMPNGbMGC1fvlwLFy7UkCFDbNtdcMEFevzxx7V58+bSsqeeekobNmwom/zj/PPP13777ac777yzbP277rpLvXr10rBhw5I5EdQs62gMxb/TnomukqNCBA14ySgDAPIqlRKLa665Rn/4wx903nnnaevWrXrggQfKHr/44oslSdddd50WLlyo0047TZMnT9b27ds1c+ZMHX/88brssstK7Xv16qUpU6Zo5syZam1t1YknnqiHH35Yzz33nH7zm984ThICVFKlanKj7tdPrTEAANUklQD5L3/5iwqFgv7whz/oD3/4Q5vHiwFy79699fTTT+uaa67R1KlT1bFjR5177rmaOXNmqf64aPr06erevbvuvvtu3X///erbt68eeOABjR07No1TQo2xm0nP2vnObqa7oOUVTm2izODnp1zDbXt202y7/e1nvwAAZFkqU01nSaWnLkQ+he1gZxckhg1Go0zxHKWDoHm/UtvtOB2X+dxr6yoDAIiq0vFaRWfSA/LEOh6wdVmc+3Hr/BanqNv2sz7BMQAgb1IfBxnIO7uAz61TnRQtgxtkamevbfhd7tTWGhBbOy6G3TYAAFlCBhkIIcxwbV7reM12l8WaXq+h7rJ4zAAAeCFABjzYjX9sfdxtXGSvIDlM570ggadbEOs3gPcTCAeZzjquKbQBAEgCJRaAT04lBm6d0cJOL+22n6Db9NverV3YoeLiKO0AACBtZJABH7yCx7i27fTvSmZb3eqMAQCoRgTIgA9u5RNBpmL22obb6BXmDn9BaqD91j67HZvXMQVFeQUAIMsIkAEPdtNK+x1Zwryenyx0HCNWhOVndA4yyACAWkCADJgE7VAWJiPs1lEv7vGUvXhlpYOMvOEnePbq8AgAQBYQIAMmTllUp8yuVyc0u3XtthFk+0E63Vkz325tw4jaYZCMNAAgiwiQgQDsaoy9gjyvLGsxkxp00g4/spqhJTAGAGQZATJgI8i4w0EeC7qdKGMIBx0rOSl2nRizPPkJAACMgwzYcCoF8DuubzEj7DVOcnGZXZuoYwj7HbvYXIccZl9+Oipa/108NjLJAIAsIoMM2AjSKc/r8TiCwKQzvGlNFx1lRkAAANJCBhkwsZYD2GVh7TKfbmUEkn2W1ilza92v07781j5HydRaz8VpO277CDu9NgAAlUKADJjYBabWQNZaMmANDs3trI+FHQs57Gx2cbUNUkbh9Lj1xwFBMgAgqwiQAQun4M1vTXDe6mq9AuMwo2xYOf24AAAgi6hBBj7kNLKC14gLbtNLWx+3Tt0c9Vi9Hvc72kXYiULCiut5AAAgCQTIgI0gmc4wj+cpkxrn8QWdeQ8AgEogQAZM3KaIdhsOLej00db9OE0g4rZ9P4IGoXHs0+8xkD0GAGQVNcjAh/yOzRtk+uco2eeoWeY4xkv2GpM56LGQQQYA5AEZZMBGXNnNMOMmxyXJ8YsBAKhmBMjAh+xKIqKWHPjJkiYVfMYxlFqcPxQIsgEAeUGJBfAhc4mFdSzj4r/9bifMY3EKWhqRdEfCvA+FBwCoLWSQAYukyyvykk1N4jgZ3g0AkAdkkIEPOY1/HHU7UduFPYYgWVqvWe7imijEaTkZZQBAlpBBBj4UV0mBYSQf8EUdm9muvfm4kyiJIAgGAOQFATJgEVdw6GfotzSC6awwZ6Fr5ZwBAPlEgAxYuE0V7Xd9t/XimIwjrSmg4w5kg0yBDQBApVCDDJiYA8KkZq9LYx9RuM0YGNe2nf4GACALCJCBD6WV1YxjP16d5qJ0qrNmuOMMYu1m0ot7HwAAREWADHwo7TGKowTKcXfSi2vdMNskOAYAZA01yADKJDVOMXXHAIC8IEBGTXML2qIGdElmRr2OO64yDgAAahElFqhplbz1n1QZRNTtMnkHAKDWkUEGPpR2xjTJodqytB0AAPKGABk1zW7M4zjGKQ6y/7DruU2NnbXxld32QSAOAMgaAmTAxDqzXZ7KLPIyHTRjIQMAso4aZKBCkgoMo2436YCVcY8BAFlHBhmokKilBU6lClFGsWAqaAAACJBR48wlFWlnNbOYQU7yuXDaJsE4ACBrCJBRs4qZ1kplTbM4ioX1OYlTkmNOAwAQJ2qQUbOKY/7mMYOc5DjISWG6aQBAXpBBRs2q1nrbrJ8PwTAAIOsIkAElV1aQtjycQx6OEQBQ2wiQAbUd/ziv8nAOjIMMAMg6AmRA+c8g5/nYAQDIGjrpIfOYWMKeXVCcx0C5eMy8xgCArCBARuYlPV5wHoNKqfx5yes5SATGAIDsIUBGzauWAK1azgMAgEojQEbmJVViUdyuOfuapyDTaarpPJ2DlN/jBgBULzrpIfOiBk5O5QfVNppCXo+/WkYQAQBUDwJkZJ7XCBN+62/N7fI+aoWdvJ5PNb4WAIB8o8QCmeeVXfSbfTS3c/p3XuW5w6G1zAUAgEpLLYO8c+dO3XDDDTr77LN14IEHql27dpo7d65t2/Xr1+uss85St27ddOCBB2r8+PHaunVrm3aGYWjGjBnq06eP6urqdPzxx+vBBx9M+lSQsjiCp2oPwPKchc3rcQMAqldqAfLWrVv1ox/9SOvXr1djY6MKDt+Kmzdv1sknn6y///3vuvnmm3Xttddq0aJF+tznPqe9e/eWtZ06daq+853vaOTIkbrjjjt0+OGH66KLLtKCBQvSOCWkKO4gyqmDGwAAQGolFj179tTrr7+uQw45RKtWrdKJJ55o227atGlqaWnR2rVr1atXL0nSiSeeqDPPPFNz5szR5ZdfLknasmWLZs2apX//93/XT3/6U0nShAkTdOqpp+raa6/V6NGjHYNw5E/UMohq65AnlZcm5LlMoRpeCwBAdUktg9yhQwcdcsghnu1+//vf69xzzy0Fx5L02c9+Vv369SvLDD/yyCPau3evJk2aVLb+pEmTtGnTJi1btiy+g0euVXsAVjy/aj9PAADSkqlRLLZs2aI333xTgwcPbvPYkCFDtGbNmtLfa9euVX19vfr379+mnWEYZW2RX8WsaLHGNsz0yuZ1rf/lmfW5ybO8Hz8AoLpkKkBubm6WJDU0NLR5rKGhQW+//bb27NlTatujRw/bdtK+YBvVwZwhtcuS+hnlwvpfNbA7n2o5NwAAKilTAXJLS4skqVOnTm0e69y5c1mblpYWX+1QHZIO/AgsAQBAUabGQa6rq5Mk7d69u81ju3btKmtTV1fnqx2yI8qU0X6mg3bavp/b90lNZ520aipNYMppAEBWZCpALpZHFEstzJqbm9W9e3d16NCh1HbJkiW27aR9o2a4mTJlivbff/+yZePGjdO4cePCHDp8SDI4Nrd1GrGi+JhdIJbXoCyvx22nms4FAODf/PnzNX/+/LJl27Ztq9DR7JOpALlnz546+OCDtXLlyjaPrVixQo2NjaW/Gxsbdd9992n9+vVlHfWWL1+uQqFQ1tbOrFmzNHDgwPgOHp7iyNLabcNtu15TVFdDUFYtmddqOQ8AQDB2CcrVq1dr0KBBFTqijNUgS9IFF1ygxx9/XJs3by4te+qpp7RhwwaNGTOmtOz888/XfvvtpzvvvLNs/bvuuku9evXSsGHDUjtm+BMm8PEzfrHbMGduHfuqJRCrhnOptg6UAIB8SzWD/POf/1zvvvtuKfh97LHH9Nprr0mSrr76anXr1k3XXXedFi5cqNNOO02TJ0/W9u3bNXPmTB1//PG67LLLStvq1auXpkyZopkzZ6q1tVUnnniiHn74YT333HP6zW9+wyQhGRQ2Y2sui4gr61st2eNqw+sCAMiCgmGk93V05JFH6h//+IftY6+88ooOO+wwSdK6det0zTXX6Nlnn1XHjh117rnnaubMmTr44IPbrDd9+nTdfffdam5uVt++fXXdddfpwgsvdDyGYsp+1apVlFjEzOsWuVfwE6STnbmu2LzMvA0/s8wRjGWDW505QTMA1J5Kx2upZpBfeeUVX+2OPvpoPfHEE77aNjU1qampKcphISZxTwdtXe4nULJ2vHMLrpEPvF4AgLRlrgYZCMJvwEx9az7wGgEAsiBTo1gg39xKLPxMB11kVz7h1M78t3X/TvsMMmxcXlRLGYK1RMZpWD4AAJJEgIxU+Mn0WgNau1pUp205LTMHWdb/V5NqCB7dRiyphvMDAOQHJRbIjSjDtbkNBYfK43UBAGQJATJi4Xc65yjbNA/1Zv23+W8/7ZE9Tq+33WMAACSJABmxcssEBs0SOt1yt8sa22WI3f5NxjL7eI0AAJVCDTJi5TQVdJHfKaGDdMZz2odbBtrtWFA5dp30AABIGwEyYuXWgc7tNrnXlNJO23Aa8cDPNpE9dll/XjcAQNoosUAmONUGBxkeLui6yA9eSwBAmgiQkQq7TnRe7fwGRX63bW2P7HHqWAkAQJoIkBELu45vYQKbMB35onb+QzZ4TTUOAEBaqEFGYoLUANu1tQuwzZN9+N02AABAEATIiIXTqBPmf/uditopOHYaF9etc541mGYEi+xzep0Z1QIAkBYCZMTCKzvsFdz4CWLtRrLwO/qFn+NENjj9wJIIkgEA6aAGGZnjNwAiUKo+fl5TOu0BAJJGgIxYBSmN8GrjZ8QLr8lA/B4jsoMppgEAlUaJBWLllQGMIzvsVrLhZ3IJMs/ZZi2xkAiUAQDpIoOMyKxj1kbJ3vrZjluW2CuDHGR8ZVSGNYPsdwIZXlcAQFzIICMypyHa7KYN9rMtr5EmwnTW8/sYKs/6+lKTDgBIGxlkZE7SgQ6ZxvzitQMApIEAGbGI0jnPqaOdn/3Z/Z8yinyzvn5Or7HdegAAxIEAGbEwd46z3hr3Mx201+Qhbvtz+7/detyKz7biaxS0TIbXFQAQFwJkJM6rVthu5jy/wU6Y7DMAAIAbAmREYi1pcJsm2LrMvL5XO7v9ui0jUM4vPyUylFgAAJJEgIxYOY1c4Wd84jDt/Owb+eP2WjqVX/CaAwDiwjBviMRpIoegw7qFCXgIkqoTryEAoNLIICOSoCUSXlND+y2x8NoHt9vzy89kM5TTAACSRICMWPidqMPpdnmQiT6C7gP55ed9wesNAIgbATJi4dT5zq1TnrlNmCyg3yw18sttlBLGvgaypfCDYB/CoO2TUPhBIRPHgeyhBhmRBekcF2dHOjKH1cuus6b53+a6dd4HQHYUg03jBj6YyDcyyIjMra44aC0xwQ6cBK1vB5AOuwys36xspTK45v0SzMMOGWREUhyFwmkkCivzLHtOY9kSJMMvgmIAQBLIICMS67TSZmECXYJjBO18yXsGqJxiFtYuGxsks5xEJtl8bNb/isgewwkBMiJxmrTB/HjavI6pVlVjZxReZyA7isGmNeiMUoKRJIJjuKHEAqHY3dq2C1bcMsxu240S+MSxDVSOU5mNdTklOkBlhAluCz8olALSqMGxeVt+2topHQvXDDggQEYo5hpiaz2x1zTBcZdloLp4lVK4Bcm8f4D0GTcYbQLR4jLzcrs2dsut3B536mjntK82x841Aw4IkBGZ3SgW5imonS5AjF6RHrsvqTzfXnQbyYL3EpAc23KJgmRE+OCZA2Xr9cmtZtn6d5gMNXcc4YQaZERiNxOe11jHSY9fSw2yPeMGI9dBsZn1PcZrDqTLfC3xHP/edO1xuw7FfX3ysz2uG3BCgIzI3MY7ruQwXNXYKa0azykqZtEDKsNcx+v2eBhhrnWh1uHaAQcEyAjMKyDxuuC4TfgQx8XKenx+6tdC7SflQDWuaVztbmfmDV9qgLOgn+0gk3pI9oGvnyDZ6d9erKNjhN0OEAQBMgKzK6kIW2Lhtm5cx5ekSgWYSXz55YnTa8ztUiAedmMG2ykFrQl/9uyCZKdh5azLvIJorhuwQ4CMwILc0k4r0+fVW9r8uFdbr/0EXScOYbPHftarpgCazDJqhZ/g1c/jbtfKQNuK6bNXDHyjZIbJMCMOBMiILOiFsVRiUcxQRLywOgW7cXZKy0o5QpBzcmvnZyD/MOx+fNi9PnHsz2kUFGqSgY84BcF+2tlJM+AMeg2vpo7IqDwCZEQW9PaU3djIbuyCKr+3/5y25fZvv9tIIsiMkt1243ULMs59Fbfl9KUc9w+NuEt0gCiS7NfgdC0s/j/IFMphjtP1R3eIz16cP2TDBsZcM+CEABmBBRmmLerFx/qFEKZkIGpGwetLJ+nMctCAPEgWPa1sSxJTzbrVmvOlh0pK4o6T3TXQbtzfKIGv03BsnnXIfN5QhQiQEZjbKBSu6xUv8AXLBd3niBNubXyNd+mRRY3ypWbN4oTZjp/boEkJGihXutykWEZh/q+43Px4qX0GymNQ3eIKit2yxGZ+r3nWUWucfvBHqdvNc1lTno8dySJARmCxZQu+/+GGvm8fVPr9svEz9I/bIPVOWWGnTLRX7+kwvM7d6Zjj+FKOOzg2d7KxPt9xl3WUssjfL5T9neZIJqhtcdbvx8Vz1IYA1zA/7fL8ecvzsSNZBMgIJOyv7bKsyPeD9ax2m3UpdN2ZR+BmPha3x+22a143zeA1zuxVmHWiZt/DKsveF5LtGAiYuZU9FJdb/3Pbjtt2i9v2+vHvp79BUJ71zDnMwubteJE+AmSEEtev7iCBp5+sR9SgOex6bscWJDCzC9rDTMsa5lziyCQHfT39bNN2fR/vv6yMPILa43YdcntPhh2hJmi5RdzCfh9UMnNL1hheCJARSpBf336ysEllPsKIu5NfGvuLdfi0ANtIcr9h69GBNCR9nfDbQS7PyOIiywiQUXHWrKnfmrgs83uMfjrMxLGfoJyGjkrqy9ptu463pU2LS8+Dj/IdIKogw6n5bZPGNuLcjlmUEotKBMnWjr0E6rBDgIwfyx6MAAAgAElEQVRAirelwt7eLgXA3zdct1GpAd/91POF4adjW1yibCvIrdog+4mrt3zZNovvRZ9lN2FueQNurO89u9KHOH7sZj0pEKVcoRKlDuZ90kkPTgiQETvPwe4zHo8kFawmlYmt5K3duAIAM7dJEHzVFZdGRwnWwZHOfNXHz/vF7+sd9f3h1IHXrZY/64FxHLL+fYDaRYBcAbX0BexWf1ztF0a7Lzm74C+u269JZN2jdroLu4/InxFTkOx2C7WSY08jPXGPUSyF7wgbtISsFoLkSqr27yGER4CcomrJUtndjgozQUat3tbK62uf9HEH6ajpWsJjuX1qt66XvL5GcOd3ODWv9ZI6njwq/gANO4pFJb4HCIrhR64D5NbWVjU1Nal3797q0qWLhg4dqsWLF1f6sCD3gehrseYrS6N0RJHYMFEOz03cHTbD7gcfyXtQV4nJdbzW5/2Xrlr7/kE4uQ6Qx48fr9tuu02XXHKJbr/9drVv317nnHOOli5dWpHjsbvoVkPG2CyOX970Gs5+UBZlLOdKsr63HP+d4hjJeXnugsjDOcVde5zmj8M8Kn72Ag0BGmKdsOxGraj17yG4y22AvGLFCi1YsEA333yzbr75Zl1++eV66qmndPjhh+vb3/52RY/Nq64xz+Nbev3ytgZ+thm7GswgF2U9MDarVEchu+fIXGPt2tvf5r1VXOb3vZvnz2eSKlWvnfSPGfPr7We0mbx8ftPmd4Qj6zppfR84XRcAJ7kNkBcuXKj27dtr4sSJpWWdOnXShAkTtGzZMm3evDnV47Fmil2HOLNZL6/MXypBvjj45Z4/WQoMgh6LOWtUqfMIU6dvXjcL7IJkp/OKO7ANsy2vH1VeHUSTzh5Xoyxe252y1Fk8VmRHbgPktWvXql+/furatWvZ8iFDhpQej1vQC35SbbMocC0o3ze5ksUAwU8mOez7LMyYyWEfc2ufZimIm6ijN0Tltc8gQxJGGcM7Tln8TIWRh0ys9RjN2e6sHzsqp32lDyCs5uZmNTQ0tFne0NAgwzC0ZcuWWPcXJavg9+Je+EGhai6aQKW5ZYeMGwxfnzdzm7Qyv053oSqRzcxKttvtGplE6UdaP0q43ldG2FE3UFtym0FuaWlRp06d2izv3Llz6XE3g+4e5Ds74/diaZeh8Jrdy25fWcgYOYnrlhS3tpAkP++vIBOXOE3y4vVYkP1kjVPfCeu/zX9X+toVaJhJy/kkNZZ4tbPr8Ob3+p5WJz27/aTZQRD5lNsAua6uTrt3726zfNeuXaXHXf2npN/s+++8887Teeedp/nz58d/oAFl/eJs/tUdJQDg1zuSZC2xCPN+c3tfBx7zO8ZhweKu63UK8L2Ow0+HXD/b9Do+6/adtucW0CM5dp8xv5+3tDrpUUqRffPnzy/FYsX/pkyZUtFjym2JRUNDg20ZRXNzsySpZ8+e7hs4S9KHTf6gP/i6uBdvyyahuO0s3FYNIuxxcYsLaTJnuaK87/xcA4KWbpiXmfdjXRY3v4GmF7e7ZGFqr4NMFpPV62Kty9r1PWvHg7bGjRuncePGlS1bvXq1Bg0aVKEjynEGubGxURs2bNCOHTvKli9fvlyFQkGNjY2Btuc3cxLkFpyfdk7lGEEzOwDKpXH71CmQ9TuKjZ/t21174gpuvdaNKwAN02HZz/E4JRDiOm4CcP9Ko8QEeMrsyjPs/g5zLHFvE7UntwHyqFGjtHfvXs2ePbu0rLW1VXPmzNHQoUPVq1cv1/VXXbkqcs/oKMM2hWE3pFKlAmcCdmRdHLdv7X7AutWqxhGohVnHek3wuj5YA0u3z3PQpICfPhhevEbNcDpeaojzxy5wjZrtDVriAdjJbYA8ZMgQjR49WlOnTlVTU5PuuecenX766Xr11Vc1Y8YM39sJe0HN0tA/BKuAt6S+DJ2CwbBj78YlbJY57s5q1iA87trtqME44hO05j+p4NStpCJKvwTUltwGyJL061//Wt/4xjc0b948TZ48Wf/617+0aNEiDR8+3HNdu7IWp5KGIBfppPgd0D5phQIBOfLDfKs1SslFEuPmhr0DFTTgdJvMI23WY87C3TA7ZKL9cxoZwuuz5tYmyCgYQR+L43qA2pDbTnqS1LFjR02fPl3Tp08PvO6qVeV/h+lQklSnPb8BeSW+UAxDKvwg9d0CoRhGuNrIwPuJWFoRRymCeVt+J86Isu8gvI4pSEIiqRE9EI41I5vG582676TXQW3KdYAcxaBB9h+UrGWOg+wrjREv4toHFylUQjX3Zk8iy50kgtvqYv1shfmsBV0n7Oe5mq8DiE/NBsjWDHJekU0B3Jm/CLPypWgukwjbB8Ltx6rddaG4TqWDY2s2OcpxVfpc8JGw9cfmMoegn0+vOmO/fwN2cl2DHIXX0Hp5vvD6DZSzUJMIJMlvPWSl5Pk6E5RbXW8tPQ+1IKufuSweE7KrZgNkqe2HJW8X6SgzfUUJjPP2PAFF1fQFGXZc5SzIynEgurDjDdutF+Tz6fbj19pZMOq+UJtqNkBetar6brPE2dEnSjsgK4pDu9Xy0E5ZHZEhq8eFYJzKF7w+a+Z2YUogrJ9tp8fsPv9JDfmI6lKzNchOnfRq4YJtF+hmfUprIAw64wDVy+3zzWcfUdV0BrkauWVlrGUXQSYsIHuMPLJ+QXJbFYiX00x4UT5rfsZQ9tPWPOYxEFTNBsi1wqsW2S6gtq5DcAwAMEtivGPrMHF+9h90u4BfBMhVyi7odXvcTz0gNYMAAMl/jXFc24vaHgiqpmuQi/ig2bMG1WnOkAQAyKY2I0DZTbrlUQNcfNxrW07bsZZZOGWf+d5CWDWbQS7WIFfzh8ZpKlm3LLDTpAOlbVXx84XaQD0iEE0co8LYjSjhNCKF3235/Rvwo2YD5GIG2c+XZZ6/UMPOTuXY0S/HzwUAILq4O7+Ztxd2XGVze76nEIeaLbFYtUoaONBf22r49Rm2drhN9rkKngvULt6/QHTm0gi/0z0H2bbkvv3i/r1KNBjqDVGQQeaXpiPbYeB4vgCgpvkdZi3qtp2275bBNmeQow43h9pWswGy5P82ER8woHrweQaS51bq4FUGYf1uNrf3GzB7rQN4qekA2a9av0Vj/jVe688F8o/3MBCNnymbg3zO7KaCtj7u9pjTv92mowa8ECCLDw8AAEGEGV3CT1u7QNhptr44RtMAnNRsJz0zv+M11iK3sSYBALXJq5NekPGH/Y4mZa4p9lMaae7wx3cXgiKD7EMtf7Ccbl8BAGpX3FNMW0sh3DLKfkon+O5CVATI4sPjhecHAGAVdrppv98pYb57qDlGXCixULy3iioh6dtH5ltZWX0OAL+43QpE5+d7wVwWYR3f2GubbsO9OY2BbG5nfZzPPIKq2QzyqlXVU+Sf9LHn+bkBrHg/A9FFyR47ZXm9vo/dRrewW5/POqKo2QBZCj7NdFbHUkxyak2m7QQAWIX9Xih+p8Q9TXWUYwLs1HSJRbX8ykzy+OkFDACwCjvMW9i6YqdSR6/SDWbTQ1g1nUGGP8xGBACIQ9TvkaDr872FsGo+QI77do/XvpLabpLbBqoJ7+l9eB4Qlp+kidt3q9tyK7uxj4O+d7n7iTBqPkD20ykgrg9XHj+k1dKRESjifbwPzwOS5DUOsdMyu2msnTrjee0XiKLmA+Qir1+ztT62Yi2fOwDgI5VMmgRNZvHdhbAIkD3EWYKRx1ua5vPP4/EDVryP9+F5QFhBSx3i/P502pbT9xTvc4RFgKzw4znGvZ8sMo9ZmcfjB6x4H+/D84A8CTv7Hu9zhFXTw7yZuX2I8jDEWdaPDwBQW+yGWIvyXUXtMdJEBln+ygeClhnYtU/i9k/SQ7AxxBuqDe/lfdxuVQNOokwQUvx/2NIMt3Wdvqt4PyMsMshyHkjcbXByv9uNo42dYlY7rWmm+XUOVBe/0/kCUVgn8gjy/gpSLuH0XcX7GWGRQbZwm7YyaEe1KL+Ug+6DDDJQ/eL8HNKpCWnwcyc16PaSbA8UkUH+kN2vz6hDxoT59ZvEtqPg1zeqTZ7f03EeOxlkhOF3PGI/60fdfxr7Q+0ig5ygIL+cs/qrmF/fqDZ5fk/HnUFOeh+oPgz5iVpBgPyhODvVxV2G4eeLLEgphN8LHBdB1Io8vNeTKHeyu+YBcTPXIUcJsMOsx3saYREgf8hvh7q4bw9F7cgXZApOp3X87JfbVKgWeZ1tiw6zyIKw34NxfJ9UoqQDtYsA2UPUMRupDQYQlzRqkIE4WQNj3nfICzrpfcitjMH8gQ4yaYjdsDbF9YOWODjt0+72qNfxRRm6DqhGeZgMKM7PLbedEVbYz4qf7zMgSwiQTfzceo2jJKL4f68Lht9A2hqA++HnPPIQNABxqNX3ud9rDADUGkosPpTGF2Te6qdqNWgAsirOgNbv0JYAUIvIIH/Ims31W3Lhta3i9szL7bZvV8rhd/te23Lbhp/SDb40UU2cPmtZfp8nNYJFkfUaBcSl0nciK71/5BcBskUcHySnKaqt/w5aPuG0bvHvuDvwcGFBNcrjezqJADavo3kgXyp9h4L3NMKixMIi7gyK07jKdsGn3b6DjMWcxADuZJRQbfL4nk7qs213bcnj84N84L2FPCGD/CFrxjSJX51e2WSnDLNdR0HraBhBp8D2m73mgoZqk9Y07XGLM4vMNNMIK80xjIFKIoNskUStn9t+ou7LWksZ55TVBMeoRnl9XydRfxxmplAAqAVkkG24lUEEWd/p77DtrBldP8u99umEL0pUM/NnOw+d9Iqsn/O4rk9567SIygn6votytxOoJDLIJn5uvfr9cJtLNdzKNqztrP9Zt2Fd12moJrdbqH7KSJIsNQGyJE/v8SQ6POXp/FF5Qd8vXt9jQFYRIFvEPRJEmvuM+wLExQzVKO/v6yQ/43l/bgAgLpRYmMRRVmBX2+d31jq3x5zKJdzGTvYqsfBbf8yXJqpJ3upunUoi4hiG0XrbO+vPBfKLYUORNwTIJkn0EnfrOOc1aoVbzZbXhcarfMLvxYoLGqqN3TjieRBX2ZNXPWheng/kC+8r5A0lFjailCq41QtbhZkoJM26Qy5oqFZud15qAcO8AYA7MsgWXpklr8yr18gSfvbr9G+7bTodi12GyG0IO7usMrfEUK3cypOy8J53mjQojrHJna5RRXGVb6Cykho1Isp7g/cV8iSVDPLrr7+u73znOzrjjDP0sY99TO3atdMzzzzj2H7p0qU66aSTVF9fr4aGBk2ePFk7d+5s0661tVVNTU3q3bu3unTpoqFDh2rx4sWRjtXrNmbQCTnsRqXw2q/dyBZ+9++nrdMxJdFDHsgipzszWXjfW4P1uI/N6RqTxL5QOUmNGhFlm7yvkCepBMgvvviibrnlFm3ZskXHHXecCi4pkLVr12rEiBHatWuXZs2apYkTJ2r27NkaM2ZMm7bjx4/XbbfdpksuuUS333672rdvr3POOUdLly4NfaxeE4WEmYjDzzrm/Vr/7XfyErsMs1Mmym5f1jZczFCt8lJW4bcPQ5Rtmq8zaU2UhGRYXz9eRyC8VEosBg8erLfeeksf//jH9dBDD2nZsmWOba+77jp1795dTz/9tOrr6yVJhx9+uK644gotXrxYI0aMkCStWLFCCxYs0K233qopU6ZIki699FINGDBA3/72t/Xss89GOuawGeQgNcNOGVuv4NQtyxOkA07QETCAapTl97pXOVeYY3e6Rpm3l+XnBM6srx+vIxBeKhnk+vp6ffzjH/dst337di1evFiXXnppKTiW9mWK6+vrtWDBgtKyhQsXqn379po4cWJpWadOnTRhwgQtW7ZMmzdvjvckPhSklthpmd2tL6df/nbZXi9+L4pkG1DL8vCjMO3PJNeAfLO7OwAgnEyNYvH8889r7969GjRoUNnyDh06qLGxUWvWrCktW7t2rfr166euXbuWtR0yZEjp8Si8xgj2Kl3wWma3Dae//XTgs9uu3X78XkC5sKIW2H1esibO4/JzTSguQ34lWS7DewO1IlMBcnNzswqFghoaGto81tDQoC1btpS1dWpnGEZZ2yD83poKWjrh1CnGbRthyjmc/guyXaAWZeVzYXcNiqsDrdP1zc81CfmRZIdr3iOoFYFrkA3DUGtrq6+2nTp1CrTtlpYWx/U6d+5cerzY1qmdeVthhA2Oo458YTfskldNcVB2dYtO+wWqmbn2Nmv8/FDnMwo75vdOFt/bQF4EziA/88wzqqur8/yvS5cu2rBhQ6Bt19XVSZJ2797d5rFdu3aVHi+2dWpn3lYS3Eos/Cx3G3vU2sbPqBRO+wtayhF0P0CeeZVK+d1G3Px+fsPs262kis98dQgyIknQEjveI6glgTPI/fv315w5c3y1tSuB8GpvGIaam5vbPNbc3KyePXuWtbUroyiua25rZ8qUKdp///3Llo0bN07jxo1zXa8YyCY9FqRdBilIJz23tmQYgHgkOc6s2/BuUUaw4DNf3YK8xnbvoajfb0AY8+fP1/z588uWbdu2rUJHs0/gALlHjx4aP358EseiAQMGqH379lq5cqVGjRpVWr5nzx6tXbtWY8eOLS1rbGzUkiVLtGPHjrKOesuXL1ehUFBjY6PrvmbNmqWBAwe2WV68OHhdJNwej7JukDZu67otsw7nRK9nwPkz4tY2iUDC6fNrDZz9Xqvctov8sr4Hw95ZsHsvEyQjbXYJytWrV7cZtCFNmeqk97GPfUwjRozQvHnzymbOmzt3rnbu3Fk2WcioUaO0d+9ezZ49u7SstbVVc+bM0dChQ9WrV69QxxClk55TRxe7jnJO24zjomTtEOjn+KyPc3FENfP6PPr5/Cf1GXHqXOv0eJA+DnbtvToFIpvcru9O3wF+tgNgn1QmCpGkG2+8UYVCQS+88IIMw9DcuXP1pz/9SZL03e9+t9Ru2rRpGj58uE455RRdccUV2rRpk2699VaNHDlSZ555ZqndkCFDNHr0aE2dOlVvvPGGjjrqKM2ZM0evvvqqfvWrXyV+Pna/sJ1+dQf9NR7HbVCnbSSZ+QLyIg/Z1KSOMa7rFLIp6DWe1x2wVzCMdD4a7dq1k90U04VCQXv37i1btnTpUjU1NWn16tXq1q2bxo4dq5tuuqls8hBpX8b4+uuv17x58/TOO+/ouOOO04033liabc9OMWW/atUq2xILP5wuQEGXx7V9t204cdq20+NANTG/351qfv2UQflpF5bdMdrtO8x1Ieo1CpXnVGJRZF7uVg5o1978t7U97xGkJY54LYrUMsgffPCB77bDhg0rZZfddOzYUdOnT9f06dOjHFpgQYdeCnpBcSvBCLINc42ieZlbe6AW2L3f4/qcxsXvdSbodcFpOZ//fHMKiL3eR06jKgG1LlM1yHAX1xcYmQHAWa0OecbnP5/sfuhFff/W4vsfsEotg1xNvG59+q1Ndtt+2BEy7Nr7XZeLImpF2PGPk862ud3h8WrrtyTErq3d7XoC5uzy8/718/7wGiOZ9wNqGRnkENx6gyddYhHmIuVnHXoyo5Y4jebit8d/0BEkghyX9f9xlFyZt+V0jWIUi/zwU34T5+vJ+wG1iAA5JLsMTNYzsH4DZaBW+AkW3dZN4jMfJPhOKlDP+rUM/vh9XwSpWQZqBQFySEl+gcQ1zae1V3LYW8hAtXK6Le32GTT/l+QxeW3ffAxxfLa9btFXUpaOJevM7wm394fTBFF274Ok3/NAFlGDnEFx3VK1ruP34kaWALXCrszCT/u0AgW32k9rBjtqnWiWP/dZPrZqYJdBZug/1DoyyBHY/aJO+osz7DSiQdYjS4BaYP38+smS+fnMB13fa9247ih5rZPVzz3Zy7bi/t6xZp3t/uY1QK0hgxxS0MxTUvv12z5MLRpQzcJ8hq1ZWz/jKZszu3GMWWzXJo5RBrLeSS+Lx1Qpfp4Lvxlgt9fdbeQToNqRQQ6JX9JA7bEbDi3sNqK2dQtmUBu87hAWf8B5PW63Xd5PqHUEyBEE7SSTB9VyHkAUUT4HYdcNG2wn0fk2i9eBLB5TJRXvGvjtiOe0Dev/eZ6BfSixCKl4YUpqmKVKqZbzAKLw6hjndus57GQ81muK17b8PB60Xdj2ScrSsWSJ3/ddmBILa/DNa4BaRAY5omrLIANwzqSl/TmPa395vj5xjY3Ob+fMrLzvgSwgQA7JbdzILMrDMQJZ5mc0iSi9/aOOOBNkrNug2/azPC4EaO6cRpXw82+nZVFHZwGqESUWIVlvP2X94sEtMiAauxEsJPvANq4p4aOOepHEuMhJX0vsSgeyfn2tBLcx7v104PQzYgXPO2oZGWQA8MFpKLW46oQBP+J4P/mdTpr3LmoZAXJIeSuxABCd33KKMNeHqLW2cV+HslRikda+88Bu5Ikw2whbngHUCkosQspbiQUA/9w+00Em8Ai77zBt/WYFox5HJbOKZDTbshv9JMq27P7m+w21iAxyDLh4ANXLLtvmNQZ62Oxe0I5RTiMQxDG+rVt2PKlrnt3zZvdvrrnuNfBFTu9dJ3F1OAWqARnkCBgnEqheYadpd8vmBc3yhp0m2Hrseckg+x0Lupavt0HG3bZ25HPq2MdzDLRFBhkAYhZlmLdKH0cltmndttMwZGQzw7HLDNsFyzy/wEcIkCOgMwNQvZzKHYLcoo66z6DrugWYUbZn9xgqw8+YxW7rhnksyD6AakGJRQR0YACql9OtbLdb3FFvS8e1fpDb8EGPJY1OyXQWc2b3/AfpOGpXAuTnNaXkArWGDDIA2IiSAQ5bEhA1AIwz4+tnCLukEAg747kB0kGAHBEXK6C6eY1a4TUurXXdKLfIgxxvXNt2GgnDz3bjOC+756xW65Gj/ADyGnkl7L6BakWAHFGUXuIA8sGtjMLPiBFebQzjo//illTZRhL7divrsL4GXHM/4me0Ez8jg8Q9ljaQZwTIMeDCAVS/oBOEOAUcaV8v4sr8FWtXgwTzYfbN9dRZ0tNMA/gIAXIMuPUEVCc/E2043br2s9y6LIvXEqfJI8KUcfgZKSFouUCtsBueLej6UR4Hag2jWMSAX+RA9QsycYbThAxhth+EdbSHqCNOuJ2bdaIJP8Fb1MdrXdRyF7f3As89UI4AOSKv2YgAVAenANDvSA9+MnhRriHWACiuDnJuy8J29nI6T7ug288x1Yo4MshR32dArSBAjogLDVCd/GZf3aZHDhKMxHktibtW1ZyRdvvbz3HElUWuxWtv1B9QAPyjBhkAIvAze52fGuMsZkadjtlvsJXkOWXx+UpaUuNkA2iLABkAbEQJJoJ0iEoiaLEbezmu7TrtJ+hxhUUmNByCYyAYSiwAwIa1xCLI1MtOndbsOlklXVqR1NBg1mV+Ms1Rx0U2P38EfMHwwwIIhgwyAESUleAj6cyxUxlJXMO8RW0PZzyXQDAEyABgI8h00H4DR7dxkeMMYOKYltnr/L2CY6fzijLOcZbHi05SUiUyAJwRIAOAjSBTG/sdI9mpxCKuqZPdth9mW+Z1/ZyL2/p+9+nn71qbarrWzhfIAmqQAcBBXMFdpYKbuPbrJxiO83jcaoxrOVCMc5xsAO7IIAOAgyi39N3KHJIaxSLqKBNO23MqDTHvw8+kIn6HunMr36jFMos4zrmWni8gDgTIAOAgSvmD3wxykqNYRN22NXMc9DyCzoYXZCrlWsokx1GGQ5kGEAwBMgA48DtKg1dnNLcgMalxkOPah9/MpTnLac14+h0XOulRMdKQ5LjWabyOAPahBhkAHPjNaPrpXGY3VrKfbQdlnebabjzmMNvzu7x4nnbHUVzu9Vz4ncI6ixnRpI4prrsBAPwhgwwAFRR3Vi9ISUNS+/NaHrVTH5lQAEkjQAaAiPyOGWxdlkSgF/doBX4DYLcyk7BjIcfRplrU0rkCWUCADAAR+SmxiGNs4kqxC86CjItsd/5e/HQIzNNzCCBfqEEGgAqw1ubGqVJjMJsDVz9BdZR91JJaPGeg0giQASAia8c0pzbWzmdJ3Ta3dohLKsDyO56zXac8px8IXsfq57mupCSe76yfM1CNKLEAgBgEGekijUAnrjIEv2Mf+y2bsBuxIuix+m1XLUOb5a0kB6gGBMgAEFHUiUTiltTQcX4eswt8/RxT3MExAERBgAwAEQXNUiad2UxqKmuvfZn/9ju9tt3EIkH257SskpIYui9r5whUO2qQAaBC0phUIumMq3VfdjW4QeuM/e7Pz/YqkXFOM4MPIBlkkAGgQtLICqaVQbZOLe2W9YwydbJdJjpIB8E0pDEONYBkkUEGgIiCZviSHl3Cuq+0t++1zPzvoM+D3XaCToedJLLHQHUggwwAiFXQOuGoWVKyrADiRoAMABVSS4GdXV1y1M6EXmUalejcFqV8xGl7ANJHgAwAFVBtY9t6TTVtZg76ojwPWX7+4jw2gmQgfakEyH/84x81YcIEfepTn1J9fb0++clPauLEiXr99ddt2y9dulQnnXSS6uvr1dDQoMmTJ2vnzp1t2rW2tqqpqUm9e/dWly5dNHToUC1evDjp0wEA2AgylnGWg9s4ENQC+ZZKJ72mpia98847Gj16tPr27au///3v+tnPfqZFixZp7dq1OuSQQ0pt165dqxEjRuiYY47RrFmztGnTJt1yyy16+eWXtWjRorLtjh8/Xg8//LCmTJmio446SnPmzNE555yjJUuWaNiwYWmcGgCEUm3TB3sFhHZTbftZz+92gqwHAF5SCZBnzZqlk046qWzZyJEjdeqpp+qOO+7QD3/4w9Ly6667Tt27d9fTTz+t+vp6SdLhhx+uK664QosXL9aIESMkSStWrNCCBQt06623asqUKZKkSy+9VAMGDNC3v/1tPfvss2mcGgCEUhzJolpYR6Zwe9y6PEjwajfuctD1khT361pt7xMgL1IpsbAGx5J08sknq3v37lq3bl1p2fbt27V48WJdeumlpeBY2pcprtGcd4oAACAASURBVK+v14IFC0rLFi5cqPbt22vixImlZZ06ddKECRO0bNkybd68OaGzAQC4STNT61auUalSjrj3S+YbSF/FOunt3LlTO3bs0EEHHVRa9vzzz2vv3r0aNGhQWdsOHTqosbFRa9asKS1bu3at+vXrp65du5a1HTJkSOlxAEC2RcmO+p2QJOp+wh4TI1kA+VWxAHnWrFnas2ePLrzwwtKy5uZmFQoFNTQ0tGnf0NCgLVu2lLV1amcYRllbAEA2xZEd9ZpqOs364ySm+SaDDKQvcA2yYRhqbW311bZTp062y5955hn98Ic/1NixY3XqqaeWlre0tDiu17lz59LjxbZO7czbAgBUnleQGiWIdVvXnEFOI9C0dkC0G+LO6ziC1HADSEbgAPmZZ57R6aef7tmuUCho3bp16tevX9ny9evX68tf/rKOO+443XPPPWWP1dXVSZJ2797dZnu7du0qPV5s69TOvC0nU6ZM0f7771+2bNy4cRo3bpzregAQl1oKetwyq1GfB6+pptMsUfDap59zraX3BSBJ8+fP1/z588uWbdu2rUJHs0/gALl///6aM2eOr7bWEojXXntNn/vc53TAAQdo0aJFZR3xiu0Nw1Bzc3ObbTU3N6tnz55lbe3KKIrrmtvamTVrlgYOHOjrPAAA8UpjmLu0h3artqH7gLTYJShXr17dpk9amgIHyD169ND48eMD7+jtt9/W5z73Oe3Zs0dLlixRjx492rQZMGCA2rdvr5UrV2rUqFGl5Xv27NHatWs1duzY0rLGxkYtWbJEO3bsKOuot3z5chUKBTU2NgY+RgBIU62OzWsNJKMGlk6lDOb/p8lr+ms/51lsR9ANVEYqnfTef/99nX322WpubtYTTzyhPn362Lb72Mc+phEjRmjevHllM+fNnTtXO3fu1JgxY0rLRo0apb1792r27NmlZa2trZozZ46GDh2qXr16JXdCABCDWg16rMOgJTHVdHEflRjqzW2/foNj8/9r9X0CVFIqE4VcdNFF+t///V9NmDBBL7zwgl544YXSY127dtX5559f+nvatGkaPny4TjnlFF1xxRXatGmTbr31Vo0cOVJnnnlmqd2QIUM0evRoTZ06VW+88UZpJr1XX31Vv/rVr9I4LQBADKJkeZPq+Bf3sYTdHkO8AZWRSoD8l7/8RYVCQb/85S/1y1/+suyxww8/vCxAPuGEE7R48WI1NTXpmmuuUbdu3TRx4kTddNNNbbb761//Wtdff73mzZund955R8cdd5wWLVqk4cOHJ35OAIBwzOUDUYJKa/BoN5uf16gScbCWQxSPJewoFua2ACojlQD5lVdeCdR+2LBh+tOf/uTZrmPHjpo+fbqmT58e9tAAACkrBohxjV7hZ8SIpLPJSW2f8gqgMio2UQgAAFIydcJxDyUXZF9R9xnXDwgA4REgAwAqKmopgV0gaS6xiHvqZ+t+rNv1GkHD6zisxw4gfQTIAICKSrI0oRLZWK/RNeLaFoDkECADAHLLK0ub5HjITpnjIOvE3R5APFLppAcAQBKcMsTWcZaTCDTDZHf9HEscI3wAiIYMMgCg6liHd8sCssdAfhAgAwCqltMU1HGwdv6z23bQMgxrxzyCZKAyKLEAAFS1NEoW7ALZMMO/2U12AiB9BMgAgKqT5jjIbtu3BuZBjoMaZKByCJABALlXzLaaZ9dzm2kvzimu7fZvbWctw3ALqIuP01EPqBwCZABA1XHK2sZVthC0ttlPwGvdBsExUDl00gMAIAC3zLRdO6e/vdoDqBwCZABAVUl6imbryBV+ppT2Gu3CaTkd9YDKoMQCAFA10q7b9TsJCR31gHwhgwwAqCpOWd00ppt2O54w016TQQYqgwwyACD3rB3x7DKvcWVjrR3ugu4vyHIyyEBlECADAHLPnGn1M1JE2oEnQ7YB+UKADACoKkkHwXalEm4BsNN0025DxRFMA5VFgAwAyD1z2UNadbtRgljKKYBsI0AGAFSdoPW/AGDGKBYAgNyzlj04jWJh/S/sfvwst+7D6d9u6wCoDDLIAIDcs5ZWeJUwRO0052eWvDBTXJPhBrKBABkAUDPiyM5ah3mLa3+VGmEDQFuUWAAAqkIxsHQLMA0jnsyx3204ZZapkQayjQAZAICMoP4YyAYCZABAVTB30PPqpGduH3QfQTrS+RnzmKAYyB5qkAEAVcGrPCGO8oWg2zC39/NvANlABhkAUBMqMYSan32SQQayhwAZAFBzwo6FHCWYdVuXIBnIFgJkAEBNCxKchimH8Br5wvw45RZANhAgAwBqUhrBqN+xkgmMgWwhQAYAVJWgI0wksW1ze6cRNfy0A1AZjGIBAKgqWZqEI8ixkEUGsoMMMgCgplUiMPWTVQZQOWSQAQA1xzCyE5CSOQayhwwyAAAAYEKADACoOWmOf+xn23TSA7KFEgsAQM0JO2V0XOUQTDUNZBsZZAAAAMCEABkAULP8ljYU28RVBmHeDqUVQPZQYgEAqFl+yxuSLLEAkD1kkAEAqBCmmQayiQwyAKBmFcsb3IJUawlE3AGtefsEy0A2ECADAGqWn4A06aCVoBjIHkosAAA1rdKd9Ip/01kPyA4CZABATQuSRY67kx7jIQPZRIAMAICHuLO7dhlpMshAdhAgAwBqmluQai59iLsMwhokk0EGsoMAGQBQ09zKHAyjvBwiiSC2uE0yyEB2ECADAGqa3zIHOtIBtYMAGQAAOQe/BMVA7SFABgBAla8BrvT+AXyEiUIAAPAhzgDWWndMcAxkCxlkAADUdlQJa81x3CNYMMQbkF2pBMh/+tOfdP755+uwww5TXV2dGhoadPbZZ2vp0qW27ZcuXaqTTjpJ9fX1amho0OTJk7Vz58427VpbW9XU1KTevXurS5cuGjp0qBYvXpz06QAAqojdyBXF/yc1goV1e2SQgWxJJUDesGGD9ttvP02aNEl33nmnrr32Wr3xxhs65ZRT9F//9V9lbdeuXasRI0Zo165dmjVrliZOnKjZs2drzJgxbbY7fvx43Xbbbbrkkkt0++23q3379jrnnHMcA28AAKzcOuclNQay234BVF7BMCrzu7WlpUV9+vTRCSecoP/4j/8oLT/nnHP0f//3f3rxxRdVX18vSbrvvvt0xRVX6Mknn9SIESMkSStWrNDQoUN16623asqUKZKk3bt3a8CAAerRo4eeffZZ2/2uXr1agwYN0qpVqzRw4MCEzxIAkGdJ1ggXCmSOASeVjtcqVoNcV1engw8+WO+++25p2fbt27V48WJdeumlpeBY2pcprq+v14IFC0rLFi5cqPbt22vixImlZZ06ddKECRO0bNkybd68OZ0TAQBUJWqEgdqVaoC8fft2vfXWW3rxxRd13XXX6YUXXihlhCXp+eef1969ezVo0KCy9Tp06KDGxkatWbOmtGzt2rXq16+funbtWtZ2yJAhpccBAPBSLJ+w/mfXLs59AsiuVId5GzNmjJ588klJUseOHXXllVfqe9/7Xunx5uZmFQoFNTQ0tFm3oaGhrGyiubnZsZ1hGNqyZUsCZwAAqDbmMgdzSYXdv5PYJ4DsCRwgG4ah1tZWX207depU9vf06dP1rW99S6+99pruv/9+tba2as+ePerYsaOkfXXJdutJUufOnUuPF9s6tTNvCwCAsJIMZKlBBrIrcID8zDPP6PTTT/dsVygUtG7dOvXr16+07Ljjjiv9++KLL9bAgQP11a9+tVRbXFdXJ2lfZzurXbt2lR4vtnVqZ94WAABu3LLDxSA2iSwygOwKHCD3799fc+bM8dXWrgSiqEOHDjrvvPM0ffp07d69W506dSqVRzQ3N7dp39zcrJ49e5Zt266Moriuua2dKVOmaP/99y9bNm7cOI0bN851PQBAbSC7C6Rj/vz5mj9/ftmybdu2Veho9gkcIPfo0UPjx4+PZefvv/++DMPQ9u3b1alTJw0YMEDt27fXypUrNWrUqFK7PXv2aO3atRo7dmxpWWNjo5YsWaIdO3aUddRbvny5CoWCGhsbXfc9a9YshnkDACReb+y2XwD2CcriMG+VksooFv/85z/bLHv33Xf10EMP6bDDDtNBBx0kSfrYxz6mESNGaN68eWUz582dO1c7d+4smyxk1KhR2rt3r2bPnl1a1traqjlz5mjo0KHq1atXgmcEAKgWdsO5uY1mkcR+AWRLKqNYnH322erdu7c+/elP65BDDtGrr76qOXPmqLm5uWxsY0maNm2ahg8frlNOOUVXXHGFNm3apFtvvVUjR47UmWeeWWo3ZMgQjR49WlOnTtUbb7yho446SnPmzNGrr76qX/3qV2mcFgCgylgzyElmeckgA9mVSoA8YcIEPfjgg7rtttv07rvv6oADDtBnPvMZXXvttRo2bFhZ2xNOOEGLFy9WU1OTrrnmGnXr1k0TJ07UTTfd1Ga7v/71r3X99ddr3rx5euedd3Tcccdp0aJFGj58eBqnBQCoIQS0QO2o2FTTlVLpqQsBANljHXKNIdiAyqp0vFaxqaYBAMgCaoEBWBEgAwCg8g565v8DqD2pTjUNAEDWWEspih31KLEAahcZZAAALAiOgdpGgAwAgEnS4x8DyD5KLAAAMCF7DIAMMgAAAGBCgAwAgAklFgAIkAEAAAATAmQAAADAhE56AACY0EkPABlkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAEwJkAAAAwIQAGQAAADAhQAYAAABMCJABAAAAk4oEyJdffrnatWun8847z/bxxx57TIMGDVJdXZ0OP/xwff/739e//vWvNu22bdumK664Qocccoi6du2qM844Q2vWrEn68AEAAFDFUg+QV61apblz56qurs728SeeeEJf+tKX1L17d91xxx360pe+pBtvvFFXX311WTvDMHTOOefowQcf1NVXX61bbrlF//znP3Xaaafp//2//5fGqSAB8+fPr/QhwAWvT3bx2mQXr0228frATuoB8tVXX62vfOUrOuSQQ2wf/+Y3v6nGxkY9+eSTmjBhgm677TZNnTpVd999tzZs2FBq97vf/U7Lli3T/fffr+9973uaNGmS/ud//kf77befbrjhhrROBzHjQpVtvD7ZxWuTXbw22cbrAzupBshz587VCy+8oGnTptk+vm7dOq1fv15XXHGF2rX76NC+9rWv6YMPPtDChQtLyx566CEdeuih+tKXvlRadtBBB2nMmDF69NFHtWfPnuROBAAAAFUrtQB5x44dmjp1qr773e86Zo/XrFmjQqGgQYMGlS1vaGhQ7969y+qL16xZo4EDB7bZxpAhQ/T++++XZZsBAAAAv1ILkH/wgx+orq5O3/jGNxzbNDc3S9oXEFs1NDRoy5YtZW2d2kkqawsAAAD41T7oCoZhqLW11VfbTp06SZI2bNig22+/Xb/97W/VoUMHx/YtLS1l65l17txZ27dvL2vr1M4wjNK2nPaxbt06X+eAdG3btk2rV6+u9GHAAa9PdvHaZBevTbbx+mRTMU5ziueSFjhAfuaZZ3T66ad7tisUClq3bp369eunb3zjGxo+fLi++MUvuq5THNli9+7dbR7btWtX2cgXdXV1ju0KhYLjKBkbN26UJF1yySWe54DKsJbYIFt4fbKL1ya7eG2yjdcnuzZu3Kjhw4envt/AAXL//v01Z84cX20bGhr0xz/+Uf/5n/+phx9+WK+++qqkfVnovXv3qqWlRa+++qq6d++ubt26lcojmpub1atXr7JtNTc369Of/nTZtoslGdZ2ktSzZ0/bYxo5cqTmzZunI444wjGIBgAAQOW0tLRo48aNGjlyZEX2HzhA7tGjh8aPH++7/WuvvaZCoVA22oS0L8O8efNm9enTR7NmzdLVV1+txsZGGYahlStXavDgwaW2zc3N2rRpk6688srSssbGRj377LNt9rd8+XJ16dJF/fr1sz2egw46SBdffLHv4wcAAED6KpE5LioYhmEkuYNNmzbZ1vZMnDhRRxxxhL73ve9pwIABOvLIIyVJxxxzjDp37qxVq1apUChIkq6//nr9+Mc/1l//+lf1799fkrRgwQKNGzdOv/vd7/TlL39ZkrR161b169dPZ599th544IEkTwsAAABVKvEA2cmRRx6pY489Vo899ljZ8kWLFun888/XaaedpgsvvFDPP/+8fv7zn2vixIn6xS9+UWr3wQcf6KSTTtILL7ygb33rWzrooIN055136h//+IdWrlypvn37pn1KAAAAqAIVC5D79OmjY489Vo8++mibxx577DH94Ac/0Lp163TwwQfrq1/9qq6//nrtt99+Ze22bduma6+9Vo888ohaWlo0ZMgQzZw5UyeccEJapwEAAIAqU7EAGQAAAMiiVKeaBgAAALKuZgLk1tZWNTU1qXfv3urSpYuGDh2qxYsXV/qwqtbKlSv19a9/XQMGDFDXrl11+OGHa+zYsXrppZfatF2/fr3OOussdevWTQceeKDGjx+vrVu3tmlnGIZmzJihPn36qK6uTscff7wefPDBNE6n6t14441q166djjvuuDaPLV26VCeddJLq6+vV0NCgyZMna+fOnW3a8RmLz+rVq3XeeefpwAMPVNeuXXXsscfqjjvuKGvD61IZL7/8si688EJ94hOfUH19vY4++mj96Ec/ajOZAa9Psnbu3KkbbrhBZ599tg488EC1a9dOc+fOtW2bxHeM323WIj+vjWEYmjNnjs4//3wddthhpevctGnTbOe4kKT77rtPxxxzjOrq6tSvX78218SiLVu2aMyYMTrggAO0//7764tf/KJeeeWV4Cdi1IixY8caHTt2NJqamox77rnHGD58uNGhQwfjueeeq/ShVaVRo0YZPXv2NCZPnmzcd999xrRp04xDDz3U6Nq1q/HCCy+U2m3atMk46KCDjL59+xp33HGH8eMf/9jo3r27ccIJJxh79uwp22ZTU5NRKBSMq666yrj33nuNL3zhC0ahUDB++9vfpn16VWXTpk1G165djW7duhnHHnts2WNr1qwx6urqjEGDBhl33323cf311xudO3c2zjnnnDbb4TMWjyeffNLo1KmT8ZnPfMa47bbbjHvvvdeYOnWq0dTUVGrD61IZr732mvHxj3/cOPLII43p06cb99xzj/Fv//ZvRqFQML74xS+W2vH6JG/jxo1GoVAwjjjiCOOMM84w2rVrZ9x///1t2iXxHRNkm7XIz2uzY8cOo1AoGMOGDTNuuukm49577zUmTJhg7LfffsYZZ5zRZpu/+MUvjEKhYIwZM8a49957ja985StGoVAwZsyY0Wa7ffv2NQ499FBj5syZxm233WYcdthhxmGHHWa8/fbbgc6jJgLkP//5z0ahUDB+8pOflJbt2rXLOOqoo4zhw4dX8Miq17Jly9pcKF566SWjU6dOxqWXXlpaNmnSJKO+vt7YtGlTadnixYuNQqFg3HPPPaVlmzdvNjp27GhcffXVZds85ZRTjMMOO8z44IMPEjqT6jd27FhjxIgRxmmnndYmQD777LONXr16GTt27Cgtu/fee4127doZ//3f/11axmcsHu+9955x6KGHGqNGjXJtx+tSGdOmTTPatWtnrFu3rmz5V77yFaNdu3bGu+++axgGr08aWltbjTfeeMMwDMNYuXKlUSgUbAPkJL5j/G6zVvl5bVpbW41ly5b9//buNqapKw4D+HPv2kI3ZFtFa0kwQMSXBE0nBhOnkS8GCeElOkGJwajTRaP105hRpwmg4ltMNqWM8dFFiYnJlmjMjEs0GqduwsQYnJjIlFWJYgLOWkp99sH02ustSB21zv5/ST/03HNPbu/D4fyht7eGfaurq6mqKk+fPq21eb1epqSksLi4WNd36dKlHDVqlDbvSHLXrl1UVZW///671tbe3k6TycTNmzdH9DriokD+8ssvaTab2dfXp2vfuXMnVVXV/ZCL6MrJyeGMGTO053a7neXl5YZ+kyZN4rx587TnBw8eDLswHT58mKqqyn9bXtOZM2doNpvZ1tZmKJB7e3tpNpu5ceNG3T79/f0cNWoUV61apbXJHBsZbrebqqryxo0bJMl//vnH8Mef5BI7GzdupKqqfPjwoa79q6++oslk4pMnTySfGBiqQI7GGjPcMcXQ2YTT1tZGRVF44MABre3EiRNUVZUnT57U9b1w4QIVReEPP/ygteXm5nLmzJmGcfPz85mVlRXRscfFNcitra2YOHEikpKSdO25ubnadvFm3L9/HykpKQCeXyfU3d2t+9bEoNzcXLS0tGjPW1tb8cEHH2hfFBPaj6SurxieZ8+eweVyYdWqVcjOzjZsb2trw8DAAHJycnTtZrMZTqfTkI/Msf/u9OnTSE5Oxp07dzB58mQkJSUhOTkZa9eu1a7Lk1xiJy8vDySxYsUK/PHHH7h79y6am5vR0NCADRs2wGq1Sj5vkWisMZGMKSLn8XgAQKsTAGjn9OU5lZOTA1VVte0kcfXq1UGzuXXrVtjPAQwmLgpkj8cDh8NhaHc4HCCJv//+OwZHFX8OHTqErq4uLF68GMCLiTBYNj09PfD7/Vpfu90eth8AyfA1uN1u/PXXX6ipqQm73ePxQFGUQfMJPecyx0bGzZs34ff7UVJSgoKCAhw7dgwrV65EQ0MDVqxYAUByiaX8/HzU1NTg1KlT+OSTTzB+/HhUVFTA5XJh7969ACSft0k01phIxhSR2717Nz788EMUFBRobR6PB++9956uaAae/9E5evRoLZuenh74fL5BswEiqxVMr/MC/m+8Xi8SEhIM7YmJidp2EV3t7e1Yt24dPv30U1RWVgJ4cd5flY3ZbJYMR1hPTw+2bduGrVu3wmazhe3zqnxCz7nkMzIeP34Mr9eLNWvWYP/+/QCA0tJS+Hw+NDY2orq6WnKJsfT0dMydOxefffYZbDYbjh8/ju3bt2PcuHFYu3at5PMWicYaE8mYIjI7duzAL7/8ArfbjeTkZK3d6/XCYrGE3Sd0Tg03m+GKiwLZarWGvW3I06dPte0ierq7u1FYWIiPP/4YR48ehaIoAF6c9+FkIxmOrM2bN2P06NFYt27doH1elU/oOZd8RkbwPAXfZQmqqKjAd999hwsXLkguMXTkyBGsXr0aHR0d2n+kSktLEQgEUFVVhSVLlkg+b5ForDGRjCmGr7m5GV9//TU+//xzrF69WrfNarWiv78/7H6hc2qks4mLSywcDof2tkioYFtqauqbPqS40dvbi/z8fPT29uLkyZMYN26cti24wAyWjc1m0/4KdzgcuHfvXth+gGQYiY6ODnz//fdwuVzo6upCZ2cnbt++jadPn8Lv96OzsxOPHj3S3uYdLJ/Qcy5zbGQEz9PLb/WOHTsWACSXGHO73Zg+fbrhLdzi4mJ4vV60tLRIPm+RaKwxkYwphufUqVNYtmwZioqK4Ha7DdsdDgcCgYDhPtN+vx8PHz7UsrHZbEhISBhyToW7/GIwcVEgO51O/Pnnn3j8+LGu/ddff4WiKHA6nTE6snebz+dDUVEROjo6cPz4cUyaNEm3PTU1FWPGjMFvv/1m2PfSpUu6XJxOJ548eYL29nZdP8kwcl1dXSAJl8uFjIwMZGRkIDMzExcvXsSNGzeQmZmJmpoaZGdnw2QyGfLx+/1obW015CNz7L8Lfgilq6tL1x68bm7s2LGSSwzdv38fgUDA0O73+0ESAwMDks9bJBprTCRjile7dOkSFixYgNzcXDQ3N0NVjWWp0+kEScM5v3z5Mp49e6adc0VRMHXq1LDZXLx4EZmZmYYPxA4ponte/E8F7zW5b98+rc3n8zErK4uzZs2K4ZG9uwKBAIuLi2mxWAy3Zgk11P0kGxsbtba7d+/SbDZz/fr1uv3nzJnDtLQ0uQ9yBB48eMAff/zR8MjOzmZ6ejp/+uknXrt2jeTQ93P9+eeftTaZYyOjpaWFiqJw6dKluvaKigpaLBZ6PB6SkkusFBUVMTExkTdv3tS1l5aW0mQyST4x8rr3QX7dNWa4Y4qhs7l+/TpTUlI4bdo03b2MX+b1emmz2cLeBzkpKYmPHj3S2oa6D/KmTZsiOva4KJBJsqysjBaLhVVVVWxsbOSsWbNosVh47ty5WB/aO2nDhg1UFIUlJSU8dOiQ4RF0584djhkzhhMmTOC3337LHTt20Gaz0el0sr+/XzdmVVUVVVXlF198waamJhYWFlJVVR45cuRNv7x3UrgvCrly5QqtViunT5/OhoYGbtmyhVarlQUFBYb9ZY6NjJUrV1JVVZaXl7O+vp6LFi2iqqrcsmWL1kdyiY2zZ8/SbDbTbrezpqaG9fX1LCgo0H4vBUk+b8aBAwdYW1vLNWvWUFEULly4kLW1taytrWVvby/J6KwxkYwZr16VTV9fH9PS0mgymbh7925DjfDyl4jU19dTVVUuWrSITU1NrKyspKqqrKur0/Xr6+vjhAkTaLfbuWfPHu7fv5/jx49nWloaHzx4ENFriJsC2efzsaqqiqmpqbRarZw5c6buG43EyMrLy6OqqoM+Ql2/fp3z589nUlISbTYbKysr2d3dHXbcuro6ZmRkMDExkVOnTuXhw4ffxMuJC3l5eZw2bZqh/fz585w9ezbff/992u12ulwu3X/GgmSOjYyBgQFWV1cz//9WqgAAAQNJREFUIyODCQkJnDhxIr/55htDP8klNi5fvszCwkKmpqYyISGBkydPZl1dHQOBgK6f5BN96enpg64xnZ2dWr9orDGRjBmPXpXN7du3h6wRli9fbhizqamJU6ZMYWJiIrOyssL+XiSffytiWVkZP/roIyYnJ7OkpIS3bt2K+DUoJBnZFSNCCCGEEEK8u+LiQ3pCCCGEEEIMlxTIQgghhBBChJACWQghhBBCiBBSIAshhBBCCBFCCmQhhBBCCCFCSIEshBBCCCFECCmQhRBCCCGECCEFshBCCCGEECGkQBZCCCGEECKEFMhCCCGEEEKEkAJZCCGEEEKIEFIgCyGEEEIIEeJfD6RB6gSsDBQAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x32f7c4150>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "2-element Array{Any,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x32fe99d90>\n",
       " PyObject <matplotlib.lines.Line2D object at 0x32fe99e10>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(tV2/(2pi), DeneV2, \",\", tV2/(2pi), DlzV2, \",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tV2, xV2, yV2, DeneV2, DlzV2 = \n",
    "    keplerIntegration( semieje, excentricidad, 2pi*10000.0, jetKepler2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAIUCAYAAADhfhFmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XucVWWh//HvxuEyzAgJBAyQmEcJPagjcIgjpGYYSYlaCIwKljiYvgzESn54NNHEQuVgZpqgQUihiGIaL/M0FpoBcbh1PB5M6QQJTJTpIS4DA7h+f+De7Mvae6+1132tz/v18iWz9tprr+uzvvvZz3qelGEYhgAAAAAU1SboFQAAAADCjtAMAAAAlEFoBgAAAMogNAMAAABlEJoBAACAMgjNAAAAQBmEZgAAAKAMQjMAAABQBqEZAAAAKIPQDAAAAJThW2jesmWLxo8fr4997GOqqanRaaedpm9/+9tqaWnJmW/VqlUaPny4ampqVFdXp6lTp2rfvn0Fy2ttbdX06dPVp08fdezYUUOHDlVTU5NfmwMAAIAESRmGYXj9Idu3b9cZZ5yhE044QV/96lfVpUsXrV69WgsWLNAll1yi5cuXS5I2bdqkc845R6effromT56s7du367777tMFF1ygFStW5Cxz/PjxWr58uaZNm6ZTTjlFCxcu1Nq1a7Vy5Uqdc845Xm8SAAAAEsSX0HzPPffo9ttv1xtvvKH+/ftnpn/5y1/WE088offee0+dO3fWqFGj9F//9V/6wx/+oJqaGknS448/rsmTJ+ull17SiBEjJElr167V0KFDNWfOHE2bNk2SdPDgQQ0YMEA9evTQa6+95vUmAQAAIEF8aZ6xZ88eSVL37t1zpvfs2VNt2rRRu3bttGfPHjU1NWnChAmZwCxJEydOVE1NjZYuXZqZtmzZMlVVVamxsTEzrX379po0aZJWr16tHTt2eLxFAAAASBJfQvP5558vwzB0zTXX6Pe//722b9+up556Sj/84Q81depUVVdX6/XXX9fhw4c1aNCgnPe2bdtW9fX12rhxY2bapk2b1K9fP9XW1ubMO2TIkMzrAAAAgFt8Cc0jR47Ut7/9bf3yl7/U2WefrRNPPFFXXHGFpkyZovvvv1+S1NzcrFQqpbq6uoL319XVaefOnZm/m5ubi85nGEbOvAAAAIBTVX590EknnaTzzjtPY8aMUZcuXbRixQrNmjVLPXv21A033JDpRaN9+/YF7+3QoUNOLxstLS1F50u/Xsy7776rl156SSeddJKqq6udbhYAAABc1tLSoq1bt2rkyJHq1q1b0KsjyafQ/OSTT2ry5MnasmVLpob40ksv1ZEjR3TLLbeooaEhE2APHjxY8P4DBw7kBNzq6uqi86VfL+all17SVVdd5Wh7AAAA4L3FixfryiuvDHo1JPkUmh955BENHDiwoEnF6NGj9eMf/1gbN27MNK1obm4ueH9zc7N69eqV+Tu/uUb2fJJy5s130kknSTp6EE477bRKNgcemjZtmubOnRv0aqAIjk94cWzCi2MTbhyfcNq8ebOuuuqqTG4LA19C865du9SlS5eC6YcOHZJhGDp8+LAGDBigqqoqrVu3TmPGjMmZZ9OmTRo3blxmWn19vVauXKm9e/fmPAy4Zs0apVIp1dfXF12XdC30aaedpoEDB7qxeXBR586dOS4hxvEJL45NeHFswo3jE25hakrry4OA/fr108aNG7Vly5ac6T/96U913HHH6cwzz1SnTp00YsQILV68OGcEwEWLFmnfvn0aO3ZsZtqYMWN0+PBhzZs3LzOttbVVCxcu1NChQ9W7d2/vNwoAAACJ4UtN8ze/+U394he/0PDhw3XjjTeqa9eueuGFF/TSSy+psbFRPXv2lCTNmjVLw4YN07nnnpsZEXDOnDkaOXKkLrzwwszyhgwZossvv1wzZszQrl27MiMCbtu2TQsWLPBjkwAAAJAgvtQ0f+pTn9KqVas0ePBgPfLII5o2bZr+9Kc/6Z577tHDDz+cme/ss89WU1OTOnbsqJtvvlnz589XY2Ojnn766YJlPvHEE7rpppu0ePFiTZ06VUeOHNGKFSs0bNgwPzYJAAAACeJbl3ODBw/Wz3/+87LznXPOOfrNb35Tdr527dpp9uzZmj17thurh5BoaGgIehVQAscnvDg24cWxCTeOD6xKGYZhBL0SftqwYYMGDRqk9evX0/AfAAAghMKY13xpngEAAABEGaEZAAAAKIPQDAAAAJRBaAYAAADKIDQDAAAAZRCaAQAAgDIIzQAAAEAZhGYAAACgDEIzAAAAUAahGQAAACiD0AwAAACUQWgGAAAAyiA0AwAAAGUQmgEAAIAyCM0AAABAGYRmAAAAoAxCMwAAAFAGoRkAAAAog9AMAAAAlEFoBgAAAMogNAMAAABlEJoBAACAMgjNAAAAQBmEZgAAAKAMQjOAxEqljv4HAEA5VUGvAAAExTCCXgMAQFRQ0wwgsahpBgBYRWgGAAAAyiA0AwAAAGUQmgEAAIAyCM0AEod2zAAAuwjNAAAAQBmEZgAAAKAMQjMAAABQBqEZAAAAKIPQDCDxeDAQAFAOoRlA4jGcNgCgHEIzgEShVhkAUAlCM4DES6UI0wCA0gjNAAAAQBmEZgCJQvtlAEAlCM0AEoVmGACAShCaASQKNc0AgEoQmgEkGiEaAGAFoRlAoqSbZ+T/v9S8AAAQmgEkitWaZQIzACAboRlAolCzDACoBKEZQKJk1zTn1zrTvhkAUAyhGUCiZNcm59csU9MMACiG0AwAWQjOAAAzhGYAiWOlGQZNNQAA2QjNAAAAQBmEZgCJYLfZBc00AADZCM0AAABAGb6G5g0bNmj06NHq2rWramtrdcYZZ+ihhx7KmWfVqlUaPny4ampqVFdXp6lTp2rfvn0Fy2ptbdX06dPVp08fdezYUUOHDlVTU5NfmwIgBmi3DACwqsqvD/qP//gPjR49WgMHDtS3vvUt1dbW6o9//KO2b9+emWfTpk0aMWKETj/9dM2dO1fbt2/Xfffdpy1btmjFihU5y5s4caKWL1+uadOm6ZRTTtHChQs1atQorVy5Uuecc45fmwUgwmiCAQCwypfQvGfPHl199dW6+OKL9fTTTxed79Zbb1WXLl30yiuvqKamRpLUt29fTZ48WU1NTRoxYoQkae3atVq6dKnmzJmjadOmSZImTJigAQMG6JZbbtFrr73m/UYBiBxCMgCgUr40z/jJT36iv/71r5o1a5Ykaf/+/TLyfhfds2ePmpqaNGHChExglo7WKNfU1Gjp0qWZacuWLVNVVZUaGxsz09q3b69JkyZp9erV2rFjh8dbBCAuyjXRIGgDACSfQvPLL7+sTp066Z133lH//v1VW1urTp066YYbbtDBgwclSa+//roOHz6sQYMG5by3bdu2qq+v18aNGzPTNm3apH79+qm2tjZn3iFDhmReB4ByygVmw6DdMwDgKF9C89tvv61Dhw7pkksu0UUXXaRnn31WkyZN0g9/+ENdc801kqTm5malUinV1dUVvL+urk47d+7M/N3c3Fx0PsMwcuYFgGLK1SKnUtQ0AwCO8qVN8969e9XS0qLrr79ec+fOlSRdeumlOnjwoObNm6e77rpLLS0tko42s8jXoUOHzOuS1NLSUnS+9OsAAACAW3wJzdXV1ZKk8ePH50y/4oor9Oijj2r16tWZedLNNbIdOHAg83p6ecXmy/68UqZNm6bOnTvnTGtoaFBDQ0PZ9wKIH8OgVhkAgrBkyRItWbIkZ9ru3bsDWpvifAnNvXr10v/8z/+oR48eOdO7d+8uSXr//fd18sknyzAMNTc3F7y/ublZvXr1yvyd31wje77055Uzd+5cDRw40NZ2AAAAwF1mlZYbNmwoeM4taL60aU5vdH6vFung2717dw0YMEBVVVVat25dzjyHDh3Spk2bVF9fn5lWX1+vt956S3v37s2Zd82aNUqlUjnzAgAAAE75EprHjh0rwzD0+OOP50x/7LHH1LZtW5133nnq1KmTRowYocWLF+eMALho0SLt27dPY8eOzUwbM2aMDh8+rHnz5mWmtba2auHChRo6dKh69+7t/UYBiLx0zxj0kAEAKMeX5hn19fW65pprtGDBAh06dEjnnXeefv3rX+uZZ57Rrbfeqp49e0qSZs2apWHDhuncc8/V5MmTtX37ds2ZM0cjR47UhRdemFnekCFDdPnll2vGjBnatWtXZkTAbdu2acGCBX5sEoAYMgvPBGoAgOTjMNqPPvqo+vbtqwULFui5555T37599cADD+hrX/taZp6zzz5bTU1Nmj59um6++WYdf/zxamxs1D333FOwvCeeeEK33367Fi9erPfff19nnnmmVqxYoWHDhvm1SQBiKP9hwFSK4AwAkFJG/tB8MZduWL5+/XoeBAQSJDsMp3vKMCv90tPz5wcA+CeMec2XNs0AECV0PQcAyEdoBoAiqGEGAKQRmgHEHjXHAACnCM0AUAahGwBAaAYAAADKIDQDSAzaKAMAKkVoBpAYUW5mEeV1B4A4IDQDAAAAZRCaASRSFJtqUNsMAMEhNANAlnSYDmOoDuM6AUBSEJoBIEu52lw/anuLfQY1zQAQHEIzgFhzM2gSWgEguQjNAGItv0lDuSYONIEAAJghNAOINbu1w2GpTQ7LegAAjiI0A0AWqzXNfoZaAjQABI/QDAAhR5MRAAgeoRkAKuB3kCU4A0CwCM0AUAEvmkyUWiZNNAAgWIRmALGVHTSpqQUAOEFoBhBLbtXMmoXt9LL9qv2llhkAgkdoBhBLZmHXrdpmv4faDvPQ3gCQFIRmAInhVo1tUDXN1DgDQHAIzQBiKSoBMyrrCQBJR2gGEEteN2WI+vIBAPYQmgHEktc1uE6Wb9bcIn951EADQLgQmgHEXhhrba2EYoIzAIQHoRkAfEQQBoBoIjQDgI+s1npnzxfGmnIASBpCMwBY4HZ3dXbmK9X2GQDgD0IzAEQAtc0AECxCMwDkcSugulkrTG0zAASL0Awg1twIwJUuw6vaYWqdAcB/hGYAsRaFWlm76xiFbQKAuKkKegUAIG7yQ62dmuFU6uj8BGMACBdqmgEkhl9NNdLzGEbh/OXeX+x1s2UBAPxDaAaQGG7U3toZyS+Vsj88ttl7rH4uAMA7hGYAiUFNLQCgUoRmAHCZ1+Gc8A8A/iM0A4ANYQisNNUAAP8RmgHABjttmit5LwAgnAjNAGKJgUUAAG4iNAOIJS9qdc0Cs5OaZwBAdBCaAcSaFzXD2X0m21m+3XUp1WczAMBfhGYAcMBOLbJbw2VTcw0A/iM0A4BNhFYASB5CM4DYCsvDgPlNOaw27aAZBgCEB6EZQGyFpUbY7lDaducDAHiP0AwADtkNwXZqkAnOABAOhGYA8IGd8GsWqmmqAQDBIjQDiC2/gmapz7Habjm7G7tKPwsA4B1CM4DY8btJQ6nPK7cu6dfdmg8A4A1CM4DIiGNgtFq7XMlgKgAA9xCaAaBC+QHWabtlAEB4EZoBwIJKQm7+e0q1WyZEA0C4BRaa7777brVp00ZnnnlmwWurVq3S8OHDVVNTo7q6Ok2dOlX79u0rmK+1tVXTp09Xnz591LFjRw0dOlRNTU1+rD6AgES5iUaU1x0Aki6Q0Lxjxw7Nnj1btbW1Ba9t2rRJI0aM0IEDBzR37lw1NjZq3rx5Gjt2bMG8EydO1AMPPKCrrrpKDz74oKqqqjRq1CitWrXKj80AgAJuBeNyNc8EcADwV1UQH/r1r39dQ4cO1eHDh/X3v/8957Vbb71VXbp00SuvvKKamhpJUt++fTV58mQ1NTVpxIgRkqS1a9dq6dKlmjNnjqZNmyZJmjBhggYMGKBbbrlFr732mr8bBQAuKheKac4BAP7yvab51Vdf1bPPPqu5c+cWvLZnzx41NTVpwoQJmcAsHa1Rrqmp0dKlSzPTli1bpqqqKjU2NmamtW/fXpMmTdLq1au1Y8cObzcEAEw47eXCyvsJzADgP19D8wcffKApU6aosbFRAwYMKHj99ddf1+HDhzVo0KCc6W3btlV9fb02btyYmbZp0yb169evoInHkCFDMq8DiLcwNVHI70fZy3VLpcK17QCQBL42z3jkkUf05z//Wb/61a9MX29ublYqlVJdXV3Ba3V1dTlNLpqbm4vOZxiGdu7c6d6KA0AZhpEbZO3UBlNzDADh51tN83vvvac77rhD3/rWt9SlSxfTeVpaWiQdbWaRr0OHDpnX0/MWmy97WQDiw6wLt2LzhSWIerUeYdk+AEgK30Lzv/3bv6lr16668cYbi85TXV0tSTp48GDBawcOHMi8np632HzZywIQH/lNEoo1UQii+UKpdan0vW6/BwBQOV+aZ2zZskXz58/X9773vcwDeoZh6MCBAzp06JC2bdumTp06ZZpWNDc3FyyjublZvXr1yvxdV1dn2gQj/d7sec1MmzZNnTt3zpnW0NCghoYG29sHwB/5TSDKzRsVUVpXAHDbkiVLtGTJkpxpu3fvDmhtivMlNO/YsUOGYWjKlCn62te+VvD6ySefrKlTp2rmzJmqqqrSunXrNGbMmMzrhw4d0qZNmzRu3LjMtPr6eq1cuVJ79+7NeRhwzZo1SqVSqq+vL7lOc+fO1cCBA13YOgB+KReYU6lkBtCkbjeAeDCrtNywYUNBxxBB86V5xoABA7R8+XItX75czz33XOa/f/7nf1bfvn313HPPadKkSerUqZNGjBihxYsX54wAuGjRIu3bty9ngJMxY8bo8OHDmjdvXmZaa2urFi5cqKFDh6p3795+bBoAn0WhKza3hsq2On9YthsA4syXmuauXbtq9OjRBdPnzp2rVCqliy++ODNt1qxZGjZsmM4991xNnjxZ27dv15w5czRy5EhdeOGFmfmGDBmiyy+/XDNmzNCuXbt0yimnaOHChdq2bZsWLFjgx2YBCEC5WtUwt/WlRhgAoiuQYbSzpfLucGeffbaamprUsWNH3XzzzZo/f74aGxv19NNPF7z3iSee0E033aTFixdr6tSpOnLkiFasWKFhw4b5tfoAfOLVA3VxkNTtBgA/pQwjWfUe6TYy69evp00zEEHpgJh+KDC7BKu0n2Q316vYZ2evt93l5r8nf1nUYAOImzDmtcBrmgHArqgGRLfW22p/1QAA9xCaAURCdt/LUW2O4NV6R3V/AECUEJoBxEqYa13DvG5AVPGlEX4hNAOIhHTgjHLwrOTmHuXtBfzANQK/EJoBRILVphlxq3WK2/YAbuMagV8IzQBiJehap6A/HwDgDUIz4CFqQPwX9D738/OzH44EYI5rBG7xZURAIKmodfRfuv/moD7byetefCaQdFwjcAs1zYCHqOHwX5D7PIj21tQ2A4A/CM0AQoUACMAqK+UFZQrcQmgGgBDjp2UACAdCM4BIiXKIjPK6A0DSEZoBj/HToH35+yz772L7Mwr7OQrrCCQN1yWsIjQDiAVqcQEAXiI0Ax4jzDmXvQ8No/g+DXKo7SC6mwNQHtce3EJoBjzGT3/O5TfPMOtmLXtaEPs8DMN7c64Bhaw07wKsIDQDCJVStUJmtcxmN0FqmgEAbiM0AwiVuNYExXW7gDDg+oIfGEYbQKiUa68cZlFYRwBAZahpBhAqxYaFpiYJABAkQjMAxABfKoBCPAQINxGaAQBALNgNxgRp2EGbZgChY7VtMG2IAUjB9tGO5KCmGUDoWK39KdZnc5KxL5BEQfbRjuQgNAMAAABlEJoBhI7dn1ij8JNsFNYRiLv869BswCSgGEIzgNCJ48M8UVhHIA5KXWv5r9G8C3YQmgFEVpRqiNxe1yhtO+AlJ9cCNc2wg9AMeITaCwDwnpOylppm2EFoBjxGLYY9VvaX2TxR6HLKq9rm7OWGefsBt1mpKXb6OpBGaAY8Ri2GdXa6moN5N1vsGySJlfOdawJuITQDPqDQtq7SfZXEfUwNGXBMEssA+IvQDCBySjXPSBoeZAIAfzCMNoDQKBb+yrVhTtcwER4BFJNdVgCVoKYZQGgUa6Nb6gn37PdwQwSix4vr1qw8oHyAU4RmAJFQaa8aAJKJ8gBuIzQDHqLQ9lfSa5I43xBVXtQKl+pVhl+mUAnaNAMeIcD4j30ORFP+tevWtZxux1xs+QRn2EFNM+ARCmMAsMar8tLs+QgvPw/xRmgGfEIhDQDWuF1emi2PMhl2EZoBn9B0AH7gPEOSmA0ln/93sb7MuVZgF6EZ8ElcazXiul1RxfFAkpTrirLUfFwrsIvQDKBipZ5OB4BKWC1LrJY/lE1wC6EZgCu8/KmTn1GBeKP8QBQQmgFUjJsRgCCVK4Moo+AmQjOAink9hDU/qwLJ4OW1TjkCtxCaXcJFCXAdAHCuktphapThB0KzTV4M9Yl44ZxAEAgNiBqrPV8AYUFotsmroT4RH0k8R9wc8rbS9yVhP1vBfkDUOT2Hk1gGwx+EZofocgv5+DWicuwrAE6VK0coZ1ApQrNNdJCOYoqdB0mo5eAaCI9yD2V69dAmUAnuqYgSQrMLkhCKUF72eZCUQp9zP5xKHReasiBMip2LVs7RYkNjl1om5z6cIDS7ICkBCaUl8TxI4jZHAaOjISrcfhiQX1ngJUIz4CIK5MrRFtw/7FvEQSXnMTXNcILQDMBXbrf95iYIRE+pwFtpMwp6zYDXfAnN69at04033qgBAwaotrZWffv21bhx4/T2228XzPvmm2/qc5/7nI4//nh17dpVEydO1Lvvvlswn2EYuvfee3XyySerurpaZ511lp588knPt4WLEMVwbniD/QrAjBtlA+UL7Kjy40Nmz56tVatW6fLLL9eZZ56pv/zlL/r+97+vgQMH6ne/+51OP/10SdKOHTv0qU99SieccIK++93vas+ePbrvvvv03//931q7dq2qqo6t7owZM3Tvvffquuuu0+DBg/Wzn/1MV1xxhdq0aaOxY8d6ti3pb8eGwU+cyJW086FUW8RKbkSVLq/Sz0sK9g+iJvs+W24+u+e2WTMwrg9YlTIM70+XNWvWaPDgwTmhd8uWLRowYIDGjh2rRYsWSZJuuOEGLVq0SH/4wx/Uu3dvSdLLL7+sCy+8UPPmzdO1114rSdq5c6c+/vGP66tf/aq+973vZZZ53nnnaevWrdq6datSRe7AGzZs0KBBg7R+/XoNHDjQ1nZkLzK91/I/hosvWbIL3FLngtWbQNRU0tSi2E2q2BfSuO0zr+Wfa/n7O67nIqLFrOzIL0udhOZi570ZroVwcpLXvOJL84yhQ4fmBGZJOuWUUzRgwABt3rw5M+3ZZ5/VF77whUxglqTPfOYz6tevn5YuXZqZ9txzz+nw4cO6/vrrc5Z5/fXXa/v27Vq9erVHW3IMT+GinLj3PxqnbYkrjhHCyOp56cf5S2CGHYE+CLhr1y5169ZN0tHa47/+9a8aPHhwwXxDhgzRxo0bM39v2rRJNTU16t+/f8F8hmHkzOuF7Iss+4EFLj5kS8L54PY2JmGfAShUrM9lrz6DsgaVCCw0L168WDt27ND48eMlSc3NzZKkurq6gnnr6ur03nvv6dChQ5l5e/ToYTqfdDSAe6nY0NnU6iBp3D7nuYaAZPLjlznKFzgVSGh+8803deONN2rYsGGaOHGiJKmlpUWS1L59+4L5O3TokDNPS0uLpfn8woUINyThPErCNgIA4smX3jOy/fWvf9XnP/95nXDCCXr66aczD+xVV1dLkg4ePFjwngMHDuTMU11dbWm+UqZNm6bOnTvnTGtoaFBDQ0PR9xTrMYOeNJLL6k98Vs4Rfi5kH7iBcw1hF8TDvunP5H4dTkuWLNGSJUtypu3evTugtSnO19D8j3/8QyNHjtQ//vEPvfbaa+rZs2fmtXTTinQzjWzNzc3q0qWL2rZtm5l35cqVpvNJUq9evcquy9y5cx31nmFlOuLPandFVs4Ruj46tg+4sVUuvd/KDSec9HMNwQni2rZyXSA4ZpWW6d4zwsS35hkHDx7UxRdfrC1btmjFihX6xCc+kfN6r1699NGPflTr1q0reO/atWtVX1+f+bu+vl779+/Xm2++mTPfmjVrlEqlcuYFooCC/CiCHJBMlIGIAl9C8wcffKCxY8dqzZo1WrZsmYYMGWI635e+9CX9/Oc/144dOzLTXn75Zb311ls5A5ZccsklOu644/Twww/nvP+HP/yhevfurXPOOcebDQFQwE7vMYRiAOUw0h/CypfmGTfffLNeeOEFjR49Wu+++65+8pOf5Lx+5ZVXSpJuvfVWLVu2TOeff76mTp2qPXv26P7779dZZ52lL3/5y5n5e/furWnTpun+++9Xa2ur/uVf/kXLly/Xb3/7W/30pz8tOrCJF/gZGVJhW7li3SdxrnAzA1B8UKhKlmH3PQzwg0r5Epp///vfK5VK6YUXXtALL7xQ8Ho6NPfp00evvPKKbr75Zs2YMUPt2rXTF77wBd1///2Z9sxps2fPVpcuXfToo4/qxz/+sU499VT95Cc/0bhx4/zYpAyG5IRZWzmz8yBJ54rVUf/iuv0ASivX5j7731ZH/bPymtVlA2Z8Cc2//vWvLc972mmn6cUXX7Q07/Tp0zV9+vRKV8t11CTCjjgX1nHeNgD+KVWWVPIa92k4EeiIgHHDhZgspY63Wz1qREW5WvV8cdr2KGG/I2h2ywarNdKVfD5f7mEXodlFXIDJ4nRI1jieL3b6ra50GirH/kRYZT9Q7KQccWNeoBhCs4uoxUkWp0Oox/F8cbJNfgyjm3TsT4RVKpX7n9X3AH4iNAM+s3NTiKukb38YcAwAwB5CswP83INKxfXcsfOzav685f6GdU6bDgFesdIMw04TjUo/G6gEodkifjqGm6J+7jhdfyu17VHfR2GRvx8JDAhSdhedxa5xr36No0yBU4RmwGVOQknUC3Unbbujvu1ILs7dytgtK9nPCBqhuULU1sAOL54GjwNuggDMeFUWJq2MhbsIzRaVantV6nUkA+eBdeWuJTjDfkSUmLW/5xxGWBGaLSrV9qrU64gOt7tLC3I5YZDfLrHctRKnbQ9SUvdjENud1H3tFTv3U/Y9gkBotoiaZsCeYr1h0FTFW+w3wBzXBpwiNFtE7RivfYVeAAAgAElEQVScKlZgx+kcyt5GO8PhAlHH+V1aJfuHkIuwITRbVKxfWfqWhVNRPGes1ho73bYo7psw8bsf3CRjf1pX6b6i/3EEjdBsUbp9Zn5NGrULsCMuBb3VX17c6M8Z3mDfImrMnpEA/ERodhkXMuyI2xevYjWYDGTintSdqZz/l5zXp/0apuPn1aAYnMPhwH5GkAjNNlipJYxLTSLscePnxjRuCsdwPVUuiU0wvLp2zPYl12nl2HeIKkKzRaXaUpV6+Anx5mYbu6gFnKitb5xYqWU2w4ARlTEr1+O+zW5zo6xknyNohGaLSrWlIignV3a/ok6GkM7/dxRYXd+obVeYVRqWJW+bAiXhGFPuO+NGWcc+R9AIzRXI/8bMt994shIyknzsk7ztUUPYcIfZfmTfAslBaHYRATrZ7Bz7JHZVmIRtDCv2vTvYj+6grERUEZrLsFuLQK0DKhGl86bSoeOjtI1h5LRpBpxjPzrjtAlbpcsA3EJoBkKCGhRYZdxx9GRxEqThjN1h4ZGL/YYoIjRbwDfb5AjqWHMDQSnZ4TgdmO3g/HJfpb+44Cj2G6KI0AzAtnLt9+nT3D9Wa5vT+5vaaXdQ01yZYvut1JDv5ZYF+IXQDMA2N2qJqGmyL7uW2U6Nc35PMARn+MlJd3Ol5qcMgd8IzQAAyxhO+iiaZ7jLzn5knyMohOYiuCgBBCl1Z8pyjXB6XtP5ZxZZRrHpZXg9nHTYyl4r60MzAeecdEMH+IXQXARtMmGV3fOAtpDH0Le5dVabY5QM2hUGZT+F7Xwotz5hW9+wq7S8TP/bMML3xQrJQWguwspFyYUbX14eW37WPcbLoZ2ToFyQzgnQJoG5krbNXl8bYTsfyvUTHLb1DTsnbZrDeH4gWQjNJVQ6hDIXdbSlUs5qQ6zOG+UaKqfrHod94KXsMFsu2Nrtgq6SLutQHueyNfyKiygjNJdQ6YXLBR8Pfvwsy7kCLxS0b86rZXY7OCfxPKZpkXVW95PVfcq+R1AIzSVU+pQ4Nc3RV+5nQDd+JozyT41ubHtSlHugL/v1YvOVC7np163O57ZKmixEvdvCgm78EnRO22V131gtE6NcdiLaCM0O8E0XTiS5tiQp222nzXC54FxOsUBcNCjP9O4gJOX4ZkviNlvldk0zEBRCMxIvyFqL/M/lhhFfZmG4VEB2vVa4SEj2YqCTJP4SEcV19ksSzwfEE6HZIsJMfKVrN+w+zEeb5tLMhsnN7z6q2LxRUa7ZRZDKNtUI8T4P2/MkVoZ4DvP+jLpKhtgGvEBotohvuvGVrmlmRCrrKrlZhbUdop/h1k6vGG4oFZwZVts6uiANFvsWYUFojin6AvZeEvet022OWs2QG2HSuMMoCK/5yzWbJz290s/MYdI0I2zHgofqUArnBMKA0BxzYbsxxgn7trww7iOz3irMerhw8mBe+j1modcsMOf/22yaVbbm9/BhQL/4EabCeB4D8B+h2aKoFZpRW18kV1DnarEH8+yGZCe10W7WLgNh5eQZEO5lCBNCswmzmgsnQ38GIbt5Bu0V3eHnMQ36/HGLne0I082xWK2zlflL1TIDSeTkeYa4lIWIB0KziTDdvBEefj+ZHwfZw2XHZTuLNduwW0sdllDt1nrEMdyU2qa4nM9+cVLTzL5GWBCaLbLTObud+f1CbXN5ZoVz2I4j3FFpm2GnA5DY/Tw/OS0j3LpWnHb96PU1S5kQHLr4Q9AIzSac1JiEpteKmR/Wes0MekWiw+wnxOzu6AI/pnDMyWAiZu+1Mwpfdg8ZYQrObnW1WMlynDaFM7s2vb5OQ1PGh5jb+8Zsn7P/EQRCs0eC/BZcrDChtrk8u78ooLww7KvsphOleqjInxamcOuHiofw/nA3xT3IhOFcDiu/ezHhWCAIhGaPeFWAlOoaC/6JezhwU5j3lZXeK6wG5ygHbKfr7kaNsZs8K3/LLDfM53rcsK8RBEJzGWH9NluybSVNMjzn5XkR1nOuUnHYnlKhMr9GOqzhOczDvrvxLIgXz5Pk12yaLTsO57dfvDo2gF8IzaCmugJe1nJQgxK8UjXQYWyXbIXV88pKeeBHm1Xby0hXJLi4blyLALIRmh0KMnCW/Ozskb6y/p1u3kETD9hFgDjGzgOAYWGnZq5YecA5ACc4fxB1hGaXefmTkdvBNj88E5xzBfXzX9gfqoraz6LZ57XTUBvmUOxUXLYtqC7nonZdBCHMTYQAK6qCXoFImJmSVHilhilkmo5CNpPSxQ4v2kSm7kzFJoxEGcfAPcXa9aa/5GX/u5JlW25Gkiqc/1iZ7N7xJqQ5R1hGXFDTbFH+g3d+1tDm1AjPrKwbLKvvCdMXATvcGOY8vw/Qol332f2siO7TYsJaAw4ffNj/u1l/5mb/9ppb/Ubb+QzOf/vc2Gfsd4QBobmUmfaaL9gdRrfcskp9ZrGHkawM3kCtmzkrT8rbEbewnEaNT7zknPcu9E0d1lFRKxWX7QDgHKG5iDAUlOUexnGzxjnqQbrS45X+iReA/5x04ZaZbyYjnyYNZTaCQmguwm4tYcH8Dgpx2zWUFtsuF/t5K+qBWXJn6HM3udU7CT9Jwmt+DwRSrimHnTbNJRGkA0W3nIijSIfm1tZWTZ8+XX369FHHjh01dOhQNTU1ef65ZQOQR4W1WbittNbZ7GfYKIXnMBeabjbLoEYFfrByPZmd19kP/znh1vDImTKMwBwIv8plykUEJdKheeLEiXrggQd01VVX6cEHH1RVVZVGjRqlVatWufo5pdr5lQqaboUn28P4Ouw1w8222chl9wHSMH85QLz5WQZ4WitJeeY5yikkRWRD89q1a7V06VJ997vf1Xe/+11de+21evnll9W3b1/dcsst7nxIiZ4qvKiVLVW42/k8Ow+x+RH6w8h2u8lKPydr/1byICk1KpWJ87lbiVJ9C1spL9xsflZsndy81tK86OWIa7JQ3B7+BIqJbGhetmyZqqqq1NjYmJnWvn17TZo0SatXr9aOHTssLcdOQZpf42yleUOlBXUlPV2kUkfnd/KtP0pNNNLCVMthNpCGlfOCEdjcRw3jMeW6Ziva/rggLOedy0E2gyjy2V6XYVyThYoNg86+QtxENjRv2rRJ/fr1U21tbc70IUOGZF4vx9Mb6szyNYxOPr9Y4VSs8EKunP0zs3jf25aaUJjME/d+scPCypePKH4RdJvVGsB0gC71C0n6tWN9x3/45cTlAG2pDEt/ptWHobnePMe9B3EW2dDc3Nysurq6gul1dXUyDEM7d+60vCy7Pz2WqmEudoMuNmR1erqtGm8X+hCGjdo1K8uycQzLBT1u7PBTsWYSTr5slC2fKgzYlr/MMhqqr8odb+5XiIvIhuaWlha1b9++YHqHDh0yr5cy6NFBJV837jBsXehWgrOfaGNmk8ObeDbbD26WmidExy9M62KGLxvOlTzGbgVRC9daqYqBSo5zGMrkJHC7nbrbywGcimxorq6u1sGDBwumHzhwIPO6U45/Zipyk7Fa6BettS7RBtFs3uz/2/3cSmtPveT1k/Zm/04PoJBKybQ5RzlWu/fL7zIrTD91+rkudprF0HbZmfx+k/06zqWOma1jWizMl3iQG+4q16bZ6TkVpnIQyVYV9ApUqq6uzrQJRnNzsySpV69epRfwC0kdpIs/cbFe+MMLSv00pZ/O/KkaGho8WNtjrNwI3CzoDaOyAse4wwhlYM753JQLNRBWaphnpgr/XUHNtNk+tbOfkyI7DBe0n7XxfokaRs/NNKSZqWPn8cyUJOddXnqhkmuNGk53sB9RzpIlS7RkyZKcabt37w5obYqLbE1zfX293nrrLe3duzdn+po1a5RKpVRfX196AZ+TdIX0/PPPS1cc/bfXgTmf08FFvP72HbXAEYbaiEr3Wam+wJOsXMgptc/Yn/aUun4K9mV27W5eTW/2A4Jl5c1j1r6/kuc+Mu9PHV2/kr2HlKrxLvbcQwjKmihhf6GchoYGPf/88zn/zZ07N+jVKhDZ0DxmzBgdPnxY8+bNy0xrbW3VwoULNXToUPXu3TvAtTvKMOy3dbbbH3P2/70WthpRt9rPVdK9X/57K31/lESptijux8ILxdoR55zj2a9baONc6S8Elsw0jjbBsPAQWsVlQ4n+re1MhzM8/I6wiGxoHjJkiC6//HLNmDFD06dP1/z58/XpT39a27Zt07333mtpGflhye2BJdLfrs1u4G7c1P3sXs6LQQKccmu7nS7HTi0xYc4b7FfnHLdnLtHNZvneLrwrU/K3y1Z3ktQ0h4Kfbe2BUiIbmiXpiSee0E033aTFixdr6tSpOnLkiFasWKFhw4ZVvEwv29M5eT1IZvvEzo3HL7YK1pnH+pXN/kUg+/8lf5JW5cfM8rmQ9TBgGPazn1/OKnkY1ckvBigiv+lFBedAzgObdt5vte9lh+vkpiSHO2qCkQSRDs3t2rXT7NmztWPHDu3fv19r1qzRiBEjKlqWFzfaYj915v+kb/ezgyycwlTjnN08xc2f74r2bJGe7HIfsFEJefn72PXAYfchLY9+wUFxTptEZc9vqymayZfa9PLs9BHs9fmR1GYExc4LN/dFEvcrwifSodmJ9det93T5lfZaEWVB9qLhVQ1PsZu1b4MnZNWIh0FBl1I+HfNStcgEZf8U61Iso0w3b6V+tTJj9sxA2B6Yzv+spJX7Es1YkByJDc1R5mdBVLK5Qt48knn/ql4Fq0puUGY3+XJyts+nfW+2r4Ou2TdrF5r/by/R/CLc7NQEFvvyafVLqdvXYdDXFoBoSGxoHlR6QMBQc/KTaEWfV0GvH8WGDXeb7W2byUAYThiGN8fSjaYZ8E5+SM1p8lCmx41K+NmkIs3sCz9NAgBkS2xoXm/SOsPtn94pcL0Ny473r4Pj7MuxtdIXbgyYDfiS/j/9VwfLURvmD7uEszoyn6OuHyup5S7C6bVFuX8U+wFxlNjQXKym2c0btFc/5QfZTsxqMw2vedmO2ernJxE/iydL2TbMNpbj5fMAVtYrex6r/atX1DMHp7Qk9gPiKbGh2Q9J+KYd9lrASkcSCwPDkGnACLS7vyIPJbq1LmE/n+BQscDsMEhXPHiJg1rtJJTvTpRrwgNEEaH5Q158K477N22rNc5BBCGzXh7cDJm+3gRKDBrhh1LnsVvHlrAcXnbKsVLXRdhqbb1upoFj4n4vRHIQmj/kRQjyK1gF8S0+6CYadvpmNeVXl3ExUK5mzXYXYEW+wBCcw6niWlyr7aFLXItm77ESwGw/LJ3XnprAbI/X9yBqqhEWhOYP8U3YfUE+xFVyFD0Csy35w7UXO6blgkaUm8okmZN2zFH6bNPzOkR9pIdZxcfJYnnA/RlhkfjQ7OXF6FeH72EsUOyEZTtBKozbGoRKv4zYDa1ehlyzIbMRDklrrpazbjHuqSaM2L+IksSHZi9/9in6c7YLn2n2kIXTIW795mnN48zigcxxv9UB7Veng3tY7Tc7fVwy88zM+z9QgVLNfOyUlVauv0rfV2w+q+WUK11hxlBUy1wgX+JDc1RrmnNGZ4tIlnEzJBdvelFi+R/WIEVlf5kxG0TGKSu9cZj1pVxuuGRqkJAtv5mPnfeUm+amUssvd14ndRjtctgniItEh+ZKCnG3P9vr90SJ7aBVrPszs8npn1xdaM8cl+Pg5pDn+UOpm31W9rw0yQg/P89zPz7L1mfw3EPF8muFre53J1+y+YIOvyQ6NPvNzRtDqZ+rSv00GdbA50qhZ6MZQVT6WrWyfn62U3a7my7Cc3S50RShkuswqGvWThMNHONFGQYEJdGhOeqFWyU3raDb40qFtY75KipASwTlgn6jI37cvVBsn7sZavNrmAnM4RC1ZyHsslJOmr5uYwjwksuBKYIyoijRodlsiNiw1sSaqWR9I7V96QfSXD4u6eUl8QZn9WHCTH+1dga2KPIFiJtjuJmVg5W8zy9eNG1zui2Z5wIidg/xg2m7dHrOQUQlOjRHUX6vGVEKfuXavWYr+Xp2rbKNGubM9Ijtt2LMHgp0ddTD9PIraN/JzTA68q8HJ71M+KGSz3V7pMLMe/MfkOW0t81pWcGXcviJ0Bxx+QV+FAvt7NrPUgMMpFIyD8mm00r07BCjMtbJCGa2blYf/lQNOOW0jLJ6/dodTbBULXGpX2is9D7jB7/GBQCSjNBswo/g6dZn2Kk1DUOg9nL4bSvLcLuPbL84Dhpl2iw77QO61LIRDZXUNttZplvvrXSZFbdtTr8W0fM5DOV+Nq++XAT9pQXJQGg2EbrujyL8mV4pegMrUxvqdh/ZfnG155X8ByM9CMx+LBvu8qLv9zCNLGilvbGtPqRDFNKKdZ8a1jbWbnc7SfkCvxCafRbGAsxvVh5CyxSqmf6VU4XNMLICcs4yE96MwGnbZrdvQNzQ4iOIWsuiTSYsrEvYaln9FMS9puJfASgjEBGEZhNeFrRu/FyZpBuBpRsjBW5JYaoRQ7Q4bbLhZl/OTtbD6vui2v1eGNbTzXIm03OSkwfGAQ8Qmk14+Q09yOYBYShYs1kJu2W31eda5TD2vGG3GznAKjeabDjuzq3CLvHc+KyiZuZ2mRb0tRXkL5iudgdKF5UIOUKzibCFIreG+45s05ASwdgw4jfcbyWK3bT9vgExRHa8ePFwYFC8unZDc75/2IQt3be9758dkPyuTAnd8BKhOUKS2DzDVF6I9nN/BNVzhuOfuMNyY0doeRWQnfRYEUalrqXAAlvB8x6F6+HLA+4EVsQcoTkAXtwk4nTjKcpBX8FR2z+VsjOADOA3u9dhpLr/TOoXU7OAnh7xz8YuSez+Q6QQmgMQVNvAJHKraUuQ7HYbxc0HTnl1vdhdrtfXbVi7ZHNFVpj1YxtNB37JbhNv4Ut89jLMBm4q9m/AL4kNzevX+/dZUe0bOFQS3I1cRb0WmNx8gHL8bJYRBrHvIs2sq06PmYZZH9bBya9sBHBYldjQPGiQP5+TSjm7ebg6qEVEynmEQ+K/qMGRKJw/Ttax+ABLPofUCOznNLvh1IsvJ/nrQGCGHYkNzX7VNDvt95Oge5TThyDD2FWcHUEO3R3l/YbgmD1YGKVzyVFlhw9BLNNTRDqkl3rmI4DeLVJ3pjyryS+2fyPzCwAiK7Gh2S/5bebC1qYvrIr1zerlELphFvWhuwEpWudUlNY1W2YkVR+atGXCa5HPMgu3TvuTL/eFpFRbaGqV4RShGbESpZosu6J6EwcApyoJ4MXelz+dGmpYlfjQnP7ZPu59/UZNqX3kxlDkUVPJOcpAI3DC7WslKteek3tCqK85F5to5ATRmUZuU5wQ7AOrn0/NM+xKfGiW/P/p3uyzqEXM5XR/FLvhsZ8Ba9y+VqJy7SV25FUL7ITM/ODqdpAOOpgjmQjN8qcGpFRBGudCNihRb8NsJm7bA4SR3zXintZ2etC22Wx902WT0/bK5T7T7vvp4xluIzR/yIsmGlH4OTKMQczqOlndv5kCPeI9aABBiPM1E+S2ZXq/sDg9fx5JxwKxxw/9FVvPtJL7MYCeOwCvEJrzeF2I2u16ya/1CbNK2xhGYduihP0Jp1wbsjri56JZX8HlgqkZs/1gWulgEqrt9ldctJs3r+9RJs08aJqBoBCas/hR65rddZqVz/NjGNmwS+8ruusDILl3bQfSlaNLTQOMO4yj62+jljkdzkt1wWZ7ND0Ls4cx5IZxnRB+hOYsUa+9SAKOEeA/N667pA3RHaV2s5WGaK/aMLuxrGLzZ0+n1hp2EZoTLCo1sVG4QQJwX/6173SQIzuf5Yb8ml0r85v9OywqCpgutLd2EmzTwTg/LAOVIDR7LH8o2UraNHsVGuMcRuO4bUFuUxz3J6LLzbbRnpSxFoJiJX0JFwvg+etvuj02Ru0rWH6JdS237wIttwjHcBmh2WP5Q2jbrSmJY9dpfojjPovjNgGVcLNNs19lbHbgtduHcbna6vz1r6SdsdWmFvnzUS4hSaqCXoGkC7LACXthF/b1S5JUitpmxJNhBFfW5DcZsNokw2kNatDvD0pU1xvhQU1zBNht0mFlWXGXlO30C/sTQYrSF+iCphIVjIznRrizct8o9jm2Rv6jbECCJDo0+32xU7jYl98m3Mn7nSwnDIJY76juKzgT9uNe6foVe58X21tJzwxO+iS2sg1+9nYBxFGiQ3OUai/cEuZtNmtb6HR9zd4f5n1QSiB9ykZ0X8GZuB73Ytvl5vY6XVaxtsN2P7fSnkacfi4QZ4lu0xz22pSk8aJtYZDtFQFEg9tN4PLLHD/6K87/3PTfBdOpTQYqluia5qQL+5cGp00z4Bz7HWHjdtOMKLDaTMPN8A+gEKE5wbwcKKASpZpmhGUdk4L9nWxhDl1unJtRPb/Ldk3nU5OM/M8DkoLQjEgI800cgH/cKguSFPgq3WdhHJUQCBKhOaSSVKADgFWUjcUVC8el9pkb3c4BSUFo9lixIVrtDj3qSZdIEai99aKbuShsNxC0MF4nbg55Hcbtc0t+22Yn28qDg8AxhGYfhK3bs7C1ZS4lCusIwB9h6houzgjKgDlCc0h5XaDH+YYRti8pAMIjCmVBFNYRSKJE99McNuk+Nf2QtP6L4/xTLJAESSqvwoDaZqCQLzXNv/rVrzRp0iR94hOfUE1Njf7pn/5JjY2N+stf/mI6/6pVqzR8+HDV1NSorq5OU6dO1b59+wrma21t1fTp09WnTx917NhRQ4cOVVNTk9eb40i6TZ7VNstehb0k3YAIzAD85LTMSeozLEDY+RKap0+frldeeUVf/OIX9f3vf18NDQ1aunSpBg4cqL/+9a85827atEkjRozQgQMHNHfuXDU2NmrevHkaO3ZswXInTpyoBx54QFdddZUefPBBVVVVadSoUVq1apUfm2VZfmFlp/lAUsKtl9tpNjx3ErEPAH84vda8uFYpBwHnfGmeMXfuXA0fPjxn2siRI3XeeefpoYce0l133ZWZfuutt6pLly565ZVXVFNTI0nq27evJk+erKamJo0YMUKStHbtWi1dulRz5szRtGnTJEkTJkzQgAEDdMstt+i1117zY9Mqwjd+/7HPAQCAE77UNOcHZkn61Kc+pS5dumjz5s2ZaXv27FFTU5MmTJiQCczS0RrlmpoaLV26NDNt2bJlqqqqUmNjY2Za+/btNWnSJK1evVo7duzwaGsqE9bQFpb1KtYlXFjWL4nY9wiD7G7T3G72EPVz3KxLuVLbFPXtBYIWWO8Z+/bt0969e9WtW7fMtNdff12HDx/WoEGDcuZt27at6uvrtXHjxsy0TZs2qV+/fqqtrc2Zd8iQIZnXwySsP4uFZb2KDZldyfpxYyjOzv4My7mBZMvuItPtZg9RP8fNug8ttU1R314gaIGF5rlz5+rQoUMaP358Zlpzc7NSqZTq6uoK5q+rq9POnTtz5i02n2EYOfOGVRDhLgqBkppmAGlJKwPsBFu7A2clbV8CbrPdptkwDLW2tlqat3379qbTX331Vd11110aN26czjvvvMz0lpaWou/r0KFD5vX0vMXmy15WlCW5ViBpXeIBgF12y0g/uzUF4sh2aH711Vf16U9/uux8qVRKmzdvVr9+/XKmv/nmm/riF7+oM888U/Pnz895rbq6WpJ08ODBguUdOHAg83p63mLzZS+rmGnTpqlz58450xoaGtTQ0FDyfZVyazhoN6SXmy5A4x5O43KT8Hs74rLfgGx+lHlhvFaTUNYjupYsWaIlS5bkTNu9e3dAa1Oc7dDcv39/LVy40NK8+c0n3nnnHX32s5/VCSecoBUrVuQ87Jee3zAMNTc3FyyrublZvXr1ypnXrAlG+r3Z85qZO3euBg4caGk73FDJN/x0Aed2Aey0zXDUxKV2JS7bAQTJjzKPaxWwx6zScsOGDQXPuAXNdmju0aOHJk6caPuD3nvvPX32s5/VoUOHtHLlSvXo0aNgngEDBqiqqkrr1q3TmDFjMtMPHTqkTZs2ady4cZlp9fX1Wrlypfbu3ZvzMOCaNWuUSqVUX19vex29RAEKN3AzBiAlo8IDCBtfHgTcv3+/LrroIjU3N+vFF1/UySefbDpfp06dNGLECC1evDhnBMBFixZp3759OQOcjBkzRocPH9a8efMy01pbW7Vw4UINHTpUvXv39m6DXGQWgCgM3RHHcBnHbQL8wvXjPvYpksSXwU2uuOIK/ed//qcmTZqkN954Q2+88UbmtdraWl1yySWZv2fNmqVhw4bp3HPP1eTJk7V9+3bNmTNHI0eO1IUXXpiZb8iQIbr88ss1Y8YM7dq1S6eccooWLlyobdu2acGCBX5slmuCaGsWtvZtYVufMOLmhKSKWvng5Fq1u62UC4B/fAnNv//975VKpfSjH/1IP/rRj3Je69u3b05oPvvss9XU1KTp06fr5ptv1vHHH6/Gxkbdc889Bct94okndPvtt2vx4sV6//33deaZZ2rFihUaNmyY59sUdWF7CDAs6xFmXrVxB8LO7fLB67LPybVKWQiEly+h+U9/+pOt+c855xz95je/KTtfu3btNHv2bM2ePbvSVQuNMAVYAEDlKM+BeApscBOEQ5hrLd1atzBvY6W4ISNJongNR3GdAZRGaA4Js+FQ/fzcqH9GmD8fQLK4MeR3ejlh49a2AVFEaE44akMAAADK86VNM2CG2orw4ZgA7qJvdSA+CM0h5ldBG2SBnj+kN4rLfrjIjyHWgbjzo8yJc7kW520DzNA8I8SS0HYsfxu92F4KduuScM4BfuJ6AuKD0AxfWL1xeBFw43LTyn5Y1Mk2xWV/AGEX92st7tsH5CM0h0SxsEgtqXPsw9ybG/sDiAeuZcBfhOYQ8rogDKKgNftMCnx70vvL6X4rVzvEcUFYuXluenmecw0B8URoDiGvf/LKX74f7VjNls9Pe/b41Zc3xwVhFZVzMyrrCcAeQjOQINSAIcrcOH+j0mMGwRsIH0JzCN87iKgAAB8cSURBVPkdbIIKUtmfS5izx8n+Yl8jiqJ43layzm41w/JKWNcL8AOhOYTiOJR2uc+lVsWeSvcXXcoB3nPSlMrKe4PsPYfyA0lGaA4Ral5hVaXnB+cV4L0wX2dhXjcg7AjNKIoahXDJvtl5UVtEt3SAO9wsO90uhynXgcoRmkPGMMLTps3P7p2C3taocdJWstL3A7DG67Iz+z7hN8oRJBmhOWSy25wGXSPgZvvXoLclbpy0laz0/QCs8aNbyDA8iwIkDaEZAACfWQ2+XoRUKkOAyhCaETh+4rOOfQVER9xHHQzDOgB+IjSHjF+FkJV2aVEZshbWZLeVZ1hzwDtJuJaoZUYSEZpDxq+CyEr71jA/AQ77+zS7rTzDmgPe8/KaCsP1GoZ1APxEaEamVsSPApBCFkDSeFHuBVGbTfmNpCM0I5DabQCAuSQ07wCiiNCcQPltWr0uoOnX013sQyBa4jSCZ7lnI4A4IzQjh9fdG1Hb7Bz7EIiW/Gs2ytdwuWcjgDgjNMOUG4UhD5yFA/sbCK9SlQphGegKwFGEZphyYyjvsA3/Ghd2j42d+Tg2gDvcbJLh9nXptFx34/4ARBGhGaa8quFI1z5Tc+KcV8cGgP/8vPacfFZ2OUF5gaQhNMMz6QKV2ohgleqHm6AMeMPOdVXuuY+wXqOU7UgaQnOIRb1Aivr6J4VhhPemDCQV5ScQPoTmEAqisPTjM83awXFj8E+xfZ0dmDkegPvKXVd+PU/g1mdQTiCpCM0hFEStnx/tY83awVHD6R8r+5rjAbivXDtgv5pJOW3L7MZygCgjNIdQ0DXNfhaI1FiEC8cDcEf+tRSGaysM6wBEGaEZkqg5wDHcWAH3UcYC0UdohiSCUtzRBzPgP6vPbxR7zYtr1o1lUpYgqaqCXgEAAOLIajvgUt1Cus2t0V6BJKKmOWT4Bg8vcJMDkI0yAbCP0IwCXvx8x5cBd9ndn1bm5xgB7nOj6YUXw2hXWoZQTiDJCM0o4MXPd9RqhB/HCHCfk6YX6YDqRZegfrwHiBtCMxBB3MAAOEEZAthHaAYiyOlPpPzECvjD7Z4xCLtAcAjNyOFnW2SCWzDM9jvHAkA52e2aKTOQRIRmlEStRvi50T6R4wyEWxiu0VJDgQNJQGiGL6iV8IfVMMzxAKKBaxUID0IzyqJWAQD84XX3kGbvtToIixufD0QZoRlF0Z45/jgOgP+cPBzo9dDalAlAcYRmBILa63DgOADe8SKAplLe9NtMe2WgPEIzEBPc7ADYYbXMoGwBjiI0oyx+rosOL4bXBuBMqessCtdgFNYR8AOhGWVRyxBedB0HAIA/CM1wXangRo2FNwjLQHJ4VY4WK0coX4CjCM0w5VaXRvnLSaUIzn5jfwPhFNS16fbQ3kBSEJqR4ebwysWWRaHsr6i3pQQAICwIzSip0iGazd5Hl0b+YV8DKIZmGEBlAgnN1157rdq0aaPRo0ebvv78889r0KBBqq6uVt++fTVz5kwdOXKkYL7du3dr8uTJ6t69u2pra3XBBRdo48aNXq9+omTXRjotUNPLoobTe5XsY44LEG5uhVo3f0EEksT30Lx+/XotWrRI1dXVpq+/+OKLuuyyy9SlSxc99NBDuuyyy3T33XdrypQpOfMZhqFRo0bpySef1JQpU3Tffffpb3/7m84//3z98Y9/9GNTEoGah+jxYvADAKBcQdJV+f2BU6ZM0dVXX62mpibT17/+9a+rvr5eL730ktq0OZrpjz/+eH3nO9/R1KlT1a9fP0nS008/rdWrV+uZZ57RZZddJkm6/PLL1a9fP91xxx1avHixPxsUc17VLFBj4S36awbCxTCi37NQFNYR8JKvNc2LFi3SG2+8oVmzZpm+vnnzZr355puaPHlyJjBL0g033KAPPvhAy5Yty0x75pln1LNnz0xglqRu3bpp7Nix+tnPfqZDhw55tyEJQq1lNHHMAABwl2+hee/evZoxY4b+7d/+Td27dzedZ+PGjUqlUho0aFDO9Lq6OvXp0yenvfLGjRs1cODAgmUMGTJE+/fv11tvveXuBgAxRe0R4B+uNyC6fAvNd955p6qrq3XTTTcVnae5uVnS0ZCcr66uTjt37syZt9h8knLmhTM8wBddHDMgPKLcP3IU1hHwmu02zYZhqLW11dK87du3lyS99dZbevDBB/XUU0+pbdu2RedvaWnJeV+2Dh06aM+ePTnzFpvPMIzMsmBfscKRgUmih2MGhEex6zEK12kU1hHwmu2a5ldffVXV1dVl/+vYsWOmicRNN92kYcOG6dJLLy257HSPGgcPHix47cCBAzk9blRXVxedL5VKFe2dA/bRPhYAnIty6IzyugNusV3T3L9/fy1cuNDSvHV1dfrVr36lX/ziF1q+fLm2bdsm6Wht9eHDh9XS0qJt27apS5cuOv744zNNK5qbm9W7d++cZTU3N+uTn/xkzrLTzTny55OkXr16lVy3adOmqXPnzjnTGhoa1NDQYGnbkqLcE98ID44VEE0EUiTdkiVLtGTJkpxpu3fvDmhtirMdmnv06KGJEydanv+dd95RKpXK6eVCklKplHbs2KGTTz5Zc+fO1ZQpU1RfXy/DMLRu3ToNHjw4M29zc7O2b9+u6667LjOtvr5er732WsHnrVmzRh07dsx0TVfM3LlzTR8khHu4EXiLkAygUnbKZ8pyeM2s0nLDhg0FHUMEzfN+mj/zmc9o+fLlBdMbGxt10kkn6bbbbtOAAQMkSaeffrr69++vefPm6brrrlPqw0Tw8MMPq02bNvrSl76Uef+YMWP0zDPP6Nlnn9UXv/hFSdK7776rZcuWafTo0SXbTsMftIHzVn5gLvc3AGSzWkZTlgNHeR6a+/Tpoz59+hRMnzp1qnr06KGLL744Z/p9992nSy65RBdeeKHGjx+v119/XT/4wQ/U2Nio/v37Z+YbM2aMHnjgAX3lK1/RG2+8oW7duunhhx/WkSNHNHPmTK83CxZQyPorv+aZmmggfvwuVylDgGN8H0Y7LZVKZWqSs33+85/Xs88+q/fff19TpkzRc889p9tuu00PPfRQznxt2rTRiy++qHHjxun73/++brnlFnXv3l0rV67Uqaee6tdmoAgCs/fy9zE1zQDcQNkBmPN9GO20//3f/y362ujRozV69Oiyy+jcubPmzZunefPmublqQOSkAzQ3OyB+/PjVqNhn8IsVcExgNc2IPrOCND2NQtZf7G8AALxFaIYldkMZzTMAIPzKle18IQeOITTDEkIwAAQnyDKY8h84itAM11AjEbzsmxs3OgAA3ENohiXFAjFBOVyyjwfHBoAbKEuAowjNcITazHDheABwG+UKcBShGRWjIA0exwBIDq+vd8oToDRCMyqWSvGzXdDY/0ByeH29U54ApRGa4RgFbXDMaoaoLQLiiWsbCBahGa4gOIcHxwJAJQjlQGmEZpRlpSClq7Pg5O9vjgWASlFmAMURmlGWlW7MqN0EgHihXAdyEZoBAIAkgjJQCqEZiAF+UgVgV7rcoEkXYA2hGWWVK0QpZAEgfijbgVyEZjjGz3nB4xgAsCtdbvCsCmANoRmWUesQXhwbAHaZNc8AUFxV0CsAwH3cBAEAcBc1zUBM8dMqAADuITQDAAAAZdA8A7ZRgwkAAJKGmmZYkt+Pp1mb2WLTAQDhVKwsL/YakGSEZtiWSlHbDABxVa4rOiCpCM1wDQUsAACIK0IzXMNPeQAAIK4IzbCFYAwAyUB5D+QiNMMWmmAAQDJQ3gO5CM0AAABAGYRmOMZPeAAAIO4IzXCMn/AAAEDcEZrhGDXNAAAg7gjNcMQwqGkGgLihXAcKEZrhGDXNABAvlOtAIUIzAADIQU0zUIjQDMQUNUUAnKIcAY4hNAMAAABlEJoBAACAMgjNAAAAQBmEZgAAkEE7ZsAcoRm2UaACAICkITTDNroiAgAASUNoRsWya5ypfQYAAHFGaAYAABn8mgiYIzTDNmqVASC+KOMBc4RmAACQQU0zYI7QDAAAAJRBaIYt/GwHAMlAeQ/kIjSjIhSmAAAgSQjNsIW2bgCQDJT3QC5CM2yhhhkA4i1dzlPeA7kIzQAAAEAZhGYAAACgDEIzAAAAUAahGQAAACiD0AwAAACU4Wtobmpq0mc+8xl95CMfUadOnTR48GA9/fTTBfM9//zzGjRokKqrq9W3b1/NnDlTR44cKZhv9+7dmjx5srp3767a2lpdcMEF2rhxox+bAgAAgATxLTQvWLBAI0eOVLt27fSd73xH999/v8477zy98847OfO9+OKLuuyyy9SlSxc99NBDuuyyy3T33XdrypQpOfMZhqFRo0bpySef1JQpU3Tffffpb3/7m84//3z98Y9/9Guz4LIlS5YEvQoogeMTXhyb8OLYhBvHB5YZPti6davRsWNHY9q0aWXnPe2004yBAwcaR44cyUy77bbbjOOOO874wx/+kJn21FNPGalUynj22Wcz0/72t78ZJ5xwgnHllVcWXf769esNScb69esr3Bp46eKLLw56FVACxye8ODbhFcVj4086CIcoHp8kCGNe86Wm+ZFHHtEHH3ygO++8U5K0b98+0/k2b96sN998U5MnT1abNsdW7YYbbtAHH3ygZcuWZaY988wz6tmzpy677LLMtG7dumns2LH62c9+pkOHDnm0NQAAxBsDmwCFfAnNL7/8svr3768VK1boYx/7mI4//nh17dpV3/rWt2RkXZkbN25UKpXSoEGDct5fV1enPn365LRX3rhxowYOHFjwWUOGDNH+/fv11ltvebdBAAAASBRfQvPbb7+tP//5z7rmmmt07bXX6plnntGoUaN0991367bbbsvM19zcLOloSM5XV1ennTt35sxbbD5JOfMCAAAATlTZfYNhGGptbbU0b/v27SVJe/fulWEYmj17tr7xjW9Iki677DL9/e9/1/e+9z3deuutqqmpUUtLS877snXo0EF79uzJ/N3S0lJ0PsMwMsvKl56+efNmS9sAf+3evVsbNmwIejVQBMcnvDg24cWxCTeOTzilc1qxPBcE26H51Vdf1ac//emy86VSKW3evFn9+vVTdXW19u/fr/Hjx+fM09DQoJdeekkbN27U8OHDVV1dLUk6ePBgwfIOHDiQeV2Sqquri86XSqVy5s22detWSdJVV11VdhsQjPzmOQgXjk94cWzCi2MTbhyf8Nq6dauGDRsW9GpIqiA09+/fXwsXLrQ0b7qpRK9evbRlyxb16NEj5/Xu3bvLMAy9//77OfM3Nzerd+/eOfM2Nzfrk5/8ZM6y08058udLf6aZkSNHavHixTrppJOKBmsAAAAEp6WlRVu3btXIkSODXpUM26G5R48emjhxoq33DBo0SFu2bNGOHTt00kknZabv2LFDqVRKH/3oRyVJ9fX1MgxD69at0+DBgzPzNTc3a/v27bruuusy0+rr6/Xaa68VfNaaNWvUsWNH9evXz3RdunXrpiuvvNLW+gMAAMBfYalhTvPlQcBx48bJMAw9/vjjmWmGYWjBggXq0qVL5meR008/Xf3799e8efNyetV4+OGH1aZNG33pS1/KTBszZox27dqlZ599NjPt3Xff1bJlyzR69Gi1bdvWhy0DAABAEqQMw5/eGC+88EL9+te/1rXXXquzzjpLy5cv18svv6x58+Zp0qRJmflWrFihSy65ROeff77Gjx+v119/XT/4wQ/U2NioRx55JDPfBx98oOHDh+uNN97QN77xDXXr1k0PP/yw/vznP2vdunU69dRT/dgsAAAAJIBvoXn//v267bbb9NRTT+m9997TJz7xCf2///f/Ch4OlKTnn39ed955pzZv3qyPfvSj+spXvqLbb79dxx13XM58u3fv1je/+U0999xzamlp0ZAhQ3T//ffr7LPP9mOTAAAAkBC+hWYAAAAgqnxp0wwAAABEWWJCc2trq6ZPn64+ffqoY8eOGjp0qJqamoJerVhYt26dbrzxRg0YMEC1tbXq27evxo0bp7fffrtg3jfffFOf+9znMkOpT5w4Ue+++27BfIZh6N5779XJJ5+s6upqnXXWWXryySdNP9/qMiHdfffdatOmjc4888yC11atWqXhw4erpqZGdXV1mjp1qvbt21cwn51ryeoyk2zDhg0aPXq0unbtqtraWp1xxhl66KGHcubh2ARjy5YtGj9+vD72sY+ppqZGp512mr797W8XDLbA8fHOvn37dMcdd+iiiy5S165d1aZNGy1atMh03iDvL3aWGSdWjo9hGFq4cKEuueQSnXjiiZlybtasWabjbUjS448/rtNPP13V1dXq169fQZmYtnPnTo0dO1YnnHCCOnfurEsvvVR/+tOfHC2zJCMhxo0bZ7Rr186YPn26MX/+fGPYsGFG27Ztjd/+9rdBr1rkjRkzxujVq5cxdepU4/HHHzdmzZpl9OzZ06itrTXeeOONzHzbt283unXrZpx66qnGQw89ZHznO98xunTpYpx99tnGoUOHcpY5ffp0I5VKGV/96leNxx57zLj44ouNVCplPPXUUznz2Vlm0m3fvt2ora01jj/+eOOMM87IeW3jxo1GdXW1MWjQIOPRRx81br/9dqNDhw7GqFGjCpZj9Vqys8ykeumll4z27dsb//qv/2o88MADxmOPPWbMmDHDmD59emYejk0w3nnnHeMjH/mI8fGPf9yYPXu2MX/+fOOaa64xUqmUcemll2bm4/h4a+vWrUYqlTJOOukk44ILLjDatGlj/PjHPy6YL+j7i9Vlxo2V47N3714jlUoZ55xzjnHPPfcYjz32mDFp0iTjuOOOMy644IKCZT7yyCNGKpUyxo4dazz22GPG1VdfbaRSKePee+8tWO6pp55q9OzZ07j//vuNBx54wDjxxBONE0880XjvvfcqWmY5iQjNv/vd74xUKmX8+7//e2bagQMHjFNOOcUYNmxYgGsWD6tXry4oQN5++22jffv2xoQJEzLTrr/+eqOmpsbYvn17ZlpTU5ORSqWM+fPnZ6bt2LHDaNeunTFlypScZZ577rnGiSeeaHzwwQe2l4mjN+wRI0YY559/fkFovuiii4zevXsbe/fuzUx77LHHjDZt2hi//OUvM9PsXEtWl5lU//jHP4yePXsaY8aMKTkfxyYYs2bNMtq0aWNs3rw5Z/rVV19ttGnTxvi///s/wzA4Pl5rbW01du3aZRiGYaxbt85IpVKmoTnI+4udZcaNlePT2tpqrF69uuC9d911l9GmTRvj5ZdfzkxraWkxunXrZoz+/+3db0xT1xsH8O+9awudyGb9U0qCAQKKCZpOF0ycZr5ZGCH8yTY1I4Zl/lskrrySmelmIkzRzZg4RxnjpYuSJSZb4mJGtmTLFqdOwbEYnJjIhHUSwQTUWkr9/l6YXrneQm/94eq4zyfpi5577pPmPD2ch/b23PJyXd/169dz5syZ2rwjyf3791NVVZ4/f15r6+7ups1m486dOx8rZjyWKJq3b99Ou93OkZERXfu+ffuoqqpuQoips2zZMr744ovac7fbzXXr1hn6LVy4kK+88or2/LPPPou5WB07doyqquo+lTEb0+p+/PFH2u12dnV1GYrm4eFh2u127tixQ3fO6OgoZ86cyc2bN2ttZudSIjGtyu/3U1VVXr58mSR5584dw+IquUmeHTt2UFVVDg4O6trfe+892mw23r17V/LzL5usaE7m+pJIzOlssvzE0tXVRUVReOTIEa3t22+/paqqPHXqlK7v6dOnqSgKv/zyS62tqKiIy5cvN8QtLi5mfn7+Y8WMxxLXNHd2dmLBggVIS0vTtRcVFWnHxdS7ceMG5syZA+DBdUcDAwO6Oz1GFRUVoaOjQ3ve2dmJGTNmoKCgwNCPpNY3kZhWdv/+ffh8PmzevBmFhYWG411dXRgbG9NuMhRlt9vh9XoNuZloLpHU5lIiMa3q+++/R3p6Oq5fv46CggKkpaUhPT0dNTU12nV+kpvkWb16NUhiw4YNuHjxIvr6+tDW1obm5mbU1tbC6XRKfp4SyV5fzMYUeoFAAAC0OgGANlaPvv+XLVsGVVW14yTx+++/T5ifq1evar8BMBvTDEsUzYFAAB6Px9Du8XhAEn///XcSXtX0dvToUfT392v7cEcnx0R5GBoaQjgc1vq63e6Y/QBo+UokppX5/X789ddfqK+vj3k8EAhAUZQJx3H8/JhsLgH63JiNaVVXrlxBOBxGRUUFSkpKcOLECWzcuBHNzc3YsGEDAMlNMhUXF6O+vh7t7e144YUXMH/+fFRVVcHn8+GTTz4BIPl5WiR7fTEbU+gdOHAAzz33HEpKSrS2QCCAZ555RldIAw/+aZw9e7Y2lkNDQwiFQqbnlJmYZthM9/wPCwaDSElJMbSnpqZqx8XU6e7uxrZt2/DSSy+huroawMMxjpcHu91uOl+JxLSqoaEh7N69Gx9++CFcLlfMPvHGcfz8mKrcyJwDbt++jWAwiK1bt+LQoUMAgMrKSoRCIbS0tGDPnj2SmyTLzs7Gyy+/jDfeeAMulwsnT57ERx99hIyMDNTU1Eh+nhLJXl+kxkjc3r178cMPP8Dv9yM9PV1rDwaDcDgcMc8Z//43m59EYpphiaLZ6XTG3Nbk3r172nExNQYGBlBaWopZs2bhq6++gqIoAB6OsZk8mM1XIjGtaufOnZg9eza2bds2YZ944zh+DKcqN1bPC/BwjB69K2pVVRU+//xznD59WnKTRMePH8eWLVvQ09OjfXJVWVmJSCSCuro6vPnmm5Kfp0Sy1xepMRLT1taGDz74AJs2bcKWLVt0x5xOJ0ZHR2OeN/79n2h+zMQ0wxKXZ3g8Hu2rlvGibZmZmf/2S5qWhoeHUVxcjOHhYZw6dQoZGRnaseiiM1EeXC6X9omwx+PBP//8E7Mf8DBficS0op6eHnzxxRfw+Xzo7+9Hb28vrl27hnv37iEcDqO3txe3bt3SLlOaaBzHzw+zcymRmFYVHYNHv9adN28eAEhukszv92Pp0qWGr3/Ly8sRDAbR0dEh+XlKJHt9MRtTAO3t7XjrrbdQVlYGv99vOO7xeBCJRAx7YYfDYQwODmpj6XK5kJKSMumciubQbEwzLFE0e71e/Pnnn7h9+7au/ddff4WiKPB6vUl6ZdNHKBRCWVkZenp6cPLkSSxcuFB3PDMzE3PnzsVvv/1mOPfs2bO6HHi9Xty9exfd3d26fo/mK5GYVtTf3w+S8Pl8yMnJQU5ODnJzc3HmzBlcvnwZubm5qK+vR2FhIWw2m2Ecw+EwOjs7DbkxM5cSiWlV0R+l9Pf369qj19fNmzdPcpNEN27cQCQSMbSHw2GQxNjYmOTnKZHs9cVsTKs7e/YsXnvtNRQVFaGtrQ2qaixBvV4vSBrG/dy5c7h//742loqiYPHixTHzc+bMGeTm5mo/ujUb0xTT+2z8h0X3xzx48KDWFgqFmJ+fzxUrViTxlU0PkUiE5eXldDgchi1dxptsz8uWlhatra+vj3a7ne+++67u/FWrVjErK8v0PprjY1rRzZs3+fXXXxsehYWFzM7O5jfffMM//viD5OT7wn733XdaWyJzyWxMq+ro6KCiKFy/fr2uvaqqig6Hg4FAgKTkJlnKysqYmprKK1eu6NorKytps9kkP0nwuPs0P+n1JZGY09lk+bl06RLnzJnDJUuWTLovcjAYpMvlirmnclpaGm/duqW1TbZP8/vvv/9YMeOxRNFMkmvXrqXD4WBdXR1bWlq4YsUKOhwO/vzzz8l+af95tbW1VBSFFRUVPHr0qOERdf36dc6dO5d5eXn89NNPuXfvXrpcLnq9Xo6Ojupi1tXVUVVVvvPOO2xtbWVpaSlVVeXx48d1/RKJKR6IdXOTCxcu0Ol0cunSpWxubuauXbvodDpZUlJiON/sXEokplVt3LiRqqpy3bp1bGpq4po1a6iqKnft2qX1kdwkx08//US73U632836+no2NTWxpKRE+7sUJfl58o4cOcKGhgZu3bqViqLw9ddfZ0NDAxsaGjg8PEwy+euL2ZjTUbz8jIyMMCsrizabjQcOHDDUCI/e+KSpqYmqqnLNmjVsbW1ldXU1VVVlY2Ojrt/IyAjz8vLodrv58ccf89ChQ5w/fz6zsrJ48+bNx4oZj2WK5lAoxLq6OmZmZtLpdHL58uWWu7PSk7J69WqqqjrhY7xLly7x1VdfZVpaGl0uF6urqzkwMBAzbmNjI3NycpiamsrFixfz2LFjMfslElM8yNeSJUsM7b/88gtXrlzJZ599lm63mz6fT/dJV1Qic8lsTKsaGxvjnj17mJOTw5SUFC5YsICHDx829JPcJMe5c+dYWlrKzMxMpqSksKCggI2NjYxEIrp+kp8nKzs7e8L1pbe3V+uX7PXFbMzpJl5+rl27NmmN8Pbbbxtitra2ctGiRUxNTWV+fn7Mv4vkg7sxrl27ls8//zzT09NZUVHBq1evxuxrNuZkFJI0fzGHEEIIIYQQ1mOJHwIKIYQQQgjx/5CiWQghhBBCiDikaBZCCCGEECIOKZqFEEIIIYSIQ4pmIYQQQggh4pCiWQghhBBCiDikaBZCCCGEECIOKZqFEEIIIYSIQ4pmIYQQQggh4pCiWQghhBBCiDikaBZCCCGEECIOKZqFEEIIIYSI439+xwE3x3iopgAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x32fed1150>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "2-element Array{Any,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x330da16d0>\n",
       " PyObject <matplotlib.lines.Line2D object at 0x330da1750>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(tV2/(2pi), DeneV2, \",\", tV2/(2pi), DlzV2, \",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(704,176)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "maximum(abs(DeneV2)), maximum(abs(DlzV2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.5631940186722204e-13,3.907985046680551e-14)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "maximum(abs(DeneV2))*eps(1.0), maximum(abs(DlzV2))*eps(1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The absolute change of this quantities, after 10000 periods, is for the energy $704\\ {\\rm eps}\\sim 1.56\\times 10^{-13}$, and for the angular momentum $176\\ {\\rm eps}\\sim 3.9\\times10^{-14}$. The changes in time of these queantities have a random walk like behaviour, pointing out that they are due to the round-off errors in the computations.\n",
    "\n",
    "The numerical precision attained is **really good**."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.8-pre",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.8"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}