{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
},
"colab": {
"name": "Foucault's Pendulum.ipynb",
"provenance": [],
"collapsed_sections": []
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "u4be2TojNEDO",
"colab_type": "text"
},
"source": [
"Importing the necessary libraries"
]
},
{
"cell_type": "code",
"metadata": {
"id": "97BfmBxuNEDU",
"colab_type": "code",
"colab": {}
},
"source": [
"from pandas import DataFrame as DF\n",
"from numpy import linspace, pi\n",
"from matplotlib import pyplot as plt"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "IfOVGEjQNEDi",
"colab_type": "text"
},
"source": [
"In 1851, the French physicist Jean Leon Foucault hung a 67-meter pendulum from the dome of the Pantheon to demonstrate the rotation of the earth for the first time. A pendulum swinging in a rotating frame exhibits deflection due to the coriolis force. \n",
"\n",
"Ignoring air resistance and losses due to imperfect pivot points etc., the equations of motion describing the tip of the pendulum on a rotating frame are given by\n",
"\n",
"$$\\frac{d^2x}{dt^2}=-\\frac{g}{l}x+2\\Omega\\frac{dy}{dt}+\\Omega^2x$$\n",
"$$\\frac{d^2y}{dt^2}=-\\frac{g}{l}y+2\\Omega\\frac{dx}{dt}+\\Omega^2y$$\n",
"\n",
"For casting the equation into first-order form, let's assume\n",
"$$\\frac{dx}{dt}=u(t)$$\n",
"$$\\frac{dy}{dt}=v(t)$$\n",
"\n",
"Therefore, we can write:\n",
"$$\\frac{d^2x}{dt^2}=-\\frac{g}{l}x+2\\Omega v+\\Omega^2x$$\n",
"$$\\frac{d^2y}{dt^2}=-\\frac{g}{l}y+2\\Omega u+\\Omega^2y$$\n",
"\n",
"Where, $\\Omega=\\frac{\\sqrt\\frac{g}{l}}{R_0}$ and $R_0$ is the dimensionless Rossby number.\n",
"\n",
"Let's write these equations in matrix form:\n",
"$$Y=\\begin{bmatrix}\n",
" x(t)\\\\\n",
" \\dot{x}(t)\\\\\n",
" y(t)\\\\\n",
" \\dot{y}(t)\\\\\n",
" \\end{bmatrix}$$\n",
" \n",
"Differentiating the above equation, we get:\n",
"$$\\dot{Y}=f(Y)=\\begin{bmatrix}\n",
" \\dot{x}(t) \\\\\n",
" \\ddot{x}(t) \\\\\n",
" \\dot{y}(t) \\\\\n",
" \\ddot{y}(t) \\\\\n",
" \\end{bmatrix}$$\n",
"\n",
"Substituting the equation for Foucault's Pendulum, we get:\n",
"$$\\dot{Y}=f(Y)=\\begin{bmatrix}\n",
" \\dot{x}(t) \\\\\n",
" \\dot{u}(t) \\\\\n",
" \\dot{y}(t) \\\\\n",
" \\dot{v}(t) \\\\\n",
" \\end{bmatrix} = \\begin{bmatrix}\n",
" u(t) \\\\\n",
" -\\frac{g}{l}x(t)+2\\Omega v(t)+\\Omega^2x(t) \\\\\n",
" v(t) \\\\\n",
" -\\frac{g}{l}y(t)+2\\Omega u(t)+\\Omega^2y(t) \\\\\n",
" \\end{bmatrix}$$\n",
" \n",
"It is given that\n",
"$$x(0)=1$$\n",
"$$y(0)=0$$\n",
"$$u(0)=0$$\n",
"$$v(0)=0$$\n",
"\n",
"# Runge-Kutta Fourth Order Method\n",
"\n",
"The RK4 method can be written as:\n",
"$$k_1=f\\left(Y_n\\right)$$\n",
"$$k_2=f\\left(Y_n+\\frac{dt}{2}k_1\\right)$$\n",
"$$k_3=f\\left(Y_n+\\frac{dt}{2}k_2\\right)$$\n",
"$$k_4=f\\left(Y_n+k_3\\right)$$\n",
"And,\n",
"$$Y_{n+1}=Y_n+\\frac{dt}{6}\\left(k_1+2k_2+2k_3+k_4\\right)$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "zXGcs073NEDl",
"colab_type": "text"
},
"source": [
"# Implementing RK4 Method"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "QyvvoXIONEDn",
"colab_type": "text"
},
"source": [
"Initializing the function `f` defined above in python"
]
},
{
"cell_type": "code",
"metadata": {
"id": "oCbe0TWyNEDp",
"colab_type": "code",
"colab": {}
},
"source": [
"def f(Y, l = 67, g = 9.81, R0 = 1):\n",
" \n",
" x = Y[0]\n",
" dx = Y[1]\n",
" y = Y[2]\n",
" dy = Y[3]\n",
" \n",
" omega = ((g/l)**0.5)/R0\n",
" \n",
" ddx = -(g/l*x)+(2*omega*dy)+(x*omega**2)\n",
" ddy = -(g/l*y)+(2*omega*dx)+(y*omega**2)\n",
" \n",
" return [\n",
" dx,\n",
" ddx,\n",
" dy,\n",
" ddy\n",
" ]"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "2gWYzSwPNED0",
"colab_type": "text"
},
"source": [
"Initializing a function to loop through the time and calculate the solution to the differential equation using RK4 method"
]
},
{
"cell_type": "code",
"metadata": {
"id": "YcPcWh0-NED2",
"colab_type": "code",
"colab": {}
},
"source": [
"def rk4(f, n, x0, dt, l = 67, g = 9.81, R0 = 1):\n",
" times = linspace(0, n*dt, n)[1:]\n",
" df = DF({\n",
" 'x': [x0[0]],\n",
" 'dx': [x0[1]],\n",
" 'y': [x0[2]],\n",
" 'dy':[x0[3]],\n",
" 't': [0]\n",
" })\n",
" for index, time in enumerate(times):\n",
" Y = df.iloc[index,:4]\n",
" k1 = f(Y, l, g, R0)\n",
" k2 = f([Y[i]+dt/2*k1[i] for i in range(len(Y))], l, g, R0)\n",
" k3 = f([Y[i]+dt/2*k2[i] for i in range(len(Y))], l, g, R0)\n",
" k4 = f([Y[i]+dt*k3[i] for i in range(len(Y))], l, g, R0)\n",
" data = [Y[i] + dt/6*(k1[i]+2*k2[i]+2*k3[i]+k4[i])\n",
" for i in range(len(Y))] + [time]\n",
" df.loc[index+1,:] = data\n",
" return df"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "wSDShb_5NEEG",
"colab_type": "text"
},
"source": [
"# Implementation\n",
"\n",
"$$\\textbf R_0 = 10$$\n",
"\n",
"## Self Convergence\n",
"\n",
"Testing the code for first second with different values of `dt`. The time is given as $t=dt\\times n$."
]
},
{
"cell_type": "code",
"metadata": {
"id": "N4iGjyr4NEEJ",
"colab_type": "code",
"outputId": "bebb6567-788e-457d-f340-1e8b9788cf78",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 335
}
},
"source": [
"l = 67\n",
"g = 9.81\n",
"R0 = 10\n",
"x0 = [1,0,0,0]\n",
"\n",
"ns = [10, 100, 1000, 10000]\n",
"dts = [0.1, 0.01, 0.001, 0.0001]\n",
"\n",
"for i in range(len(ns)):\n",
" n = ns[i]\n",
" dt = dts[i]\n",
" data = rk4(f, n, x0, dt, l, g, R0)\n",
" print('dt = {}'.format(dt))\n",
" display(data.tail(1))"
],
"execution_count": 0,
"outputs": [
{
"output_type": "stream",
"text": [
"dt = 0.1\n"
],
"name": "stdout"
},
{
"output_type": "display_data",
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
"\n",
"\n",
" | \n",
"x | \n",
"dx | \n",
"y | \n",
"dy | \n",
"t | \n",
"
\n",
"\n",
"\n",
"\n",
"9 | \n",
"0.941843 | \n",
"-0.128022 | \n",
"-0.001332 | \n",
"-0.004407 | \n",
"1.0 | \n",
"
\n",
"\n",
"
\n",
"
"
],
"text/plain": [
" x dx y dy t\n",
"9 0.941843 -0.128022 -0.001332 -0.004407 1.0"
]
},
"metadata": {
"tags": []
}
},
{
"output_type": "stream",
"text": [
"dt = 0.01\n"
],
"name": "stdout"
},
{
"output_type": "display_data",
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
"\n",
"\n",
" | \n",
"x | \n",
"dx | \n",
"y | \n",
"dy | \n",
"t | \n",
"
\n",
"\n",
"\n",
"\n",
"99 | \n",
"0.929769 | \n",
"-0.140265 | \n",
"-0.001769 | \n",
"-0.005311 | \n",
"1.0 | \n",
"
\n",
"\n",
"
\n",
"
"
],
"text/plain": [
" x dx y dy t\n",
"99 0.929769 -0.140265 -0.001769 -0.005311 1.0"
]
},
"metadata": {
"tags": []
}
},
{
"output_type": "stream",
"text": [
"dt = 0.001\n"
],
"name": "stdout"
},
{
"output_type": "display_data",
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
"\n",
"\n",
" | \n",
"x | \n",
"dx | \n",
"y | \n",
"dy | \n",
"t | \n",
"
\n",
"\n",
"\n",
"\n",
"999 | \n",
"0.928501 | \n",
"-0.141481 | \n",
"-0.001817 | \n",
"-0.005406 | \n",
"1.0 | \n",
"
\n",
"\n",
"
\n",
"
"
],
"text/plain": [
" x dx y dy t\n",
"999 0.928501 -0.141481 -0.001817 -0.005406 1.0"
]
},
"metadata": {
"tags": []
}
},
{
"output_type": "stream",
"text": [
"dt = 0.0001\n"
],
"name": "stdout"
},
{
"output_type": "display_data",
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
"\n",
"\n",
" | \n",
"x | \n",
"dx | \n",
"y | \n",
"dy | \n",
"t | \n",
"
\n",
"\n",
"\n",
"\n",
"9999 | \n",
"0.928374 | \n",
"-0.141602 | \n",
"-0.001822 | \n",
"-0.005415 | \n",
"1.0 | \n",
"
\n",
"\n",
"
\n",
"
"
],
"text/plain": [
" x dx y dy t\n",
"9999 0.928374 -0.141602 -0.001822 -0.005415 1.0"
]
},
"metadata": {
"tags": []
}
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "MdNgrn1_f736",
"colab_type": "text"
},
"source": [
"From the results above, we can see that the change is not very large after $dt=0.01$ which does matches the expectation, however, a little lower value of $dt$ was expected to ensure the convergence of RK4 solution to the second order differential equations for the Foulcault's Pendulum. Therefore, we will be using this value of $dt$ for the analysis ahead."
]
},
{
"cell_type": "code",
"metadata": {
"id": "8RYjAulZgMFB",
"colab_type": "code",
"colab": {}
},
"source": [
"dt = 0.01\n",
"n = int(round(((l/g)**0.5)/dt))"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "dxn2DakggRNH",
"colab_type": "code",
"outputId": "568058d8-8d00-4e29-ca6b-7d8c7d515911",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 851
}
},
"source": [
"R0 = 10\n",
"data = rk4(f, n, x0, dt, l, g, R0)\n",
"\n",
"data.plot(x='t', y='y')\n",
"plt.ylabel('Y')\n",
"plt.xlabel('t')\n",
"plt.legend('')\n",
"plt.title('y vs t')\n",
"plt.show()\n",
"\n",
"data.plot(x='t', y='x')\n",
"plt.ylabel('X')\n",
"plt.xlabel('t')\n",
"plt.legend('')\n",
"plt.title('x vs t')\n",
"plt.show()\n",
"\n",
"data.plot(x='x', y='y')\n",
"plt.ylabel('Y')\n",
"plt.xlabel('X')\n",
"plt.legend('')\n",
"plt.title('y vs x')\n",
"plt.show()"
],
"execution_count": 0,
"outputs": [
{
"output_type": "display_data",
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEWCAYAAACqitpwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXhV5bn+8e+TmYQwJIQxQEBARRGECIqzIGK14liHo0WrtbVaaz39ndraFqv2HGytWlurUofiWC0OoDghDlitSgCZkSCKBgMEwhCGkOn5/ZGFjbhDAuxk7ezcn+taV/Za+11rP687cme9azJ3R0REJJoSwi5ARETij8JFRESiTuEiIiJRp3AREZGoU7iIiEjUKVxERCTqFC4iIhJ1CheROGNmN5nZY2HXIa2bwkVERKJO4SLShMzs/5nZM7stu9vM/hSh7c/NbMpuy/5kZncHry81s5VmVmZmn5rZf0XYxljgl8D5ZrbVzOZHt0cijWO6/YtI0zGzbsAKoIe7bzKzJOBL4FR3n7Nb297AUqCLu5eZWSJQBJwFLASKgSPc/eNgu1nuvjjCZ94E9HP3i5uybyJ7oj0XkSbk7sXALOC8YNFYYP3uwRK0XQXMpTZMAE4Ctrv7+8F8DXCombVx9+JIwSISKxQuIk1vMrBrL+Ji4NE9tH0CuDB4fVEwj7tvA84HfggUm9l0MzuoacoV2X8aFhNpYmaWRu2Q1rHA+8BAd/+8nrY5wOdAf2ARcJS7L92tTRvgVmC4ux8bYRsTgP4aFpMwac9FpIm5ezkwhdq9kA/rC5agbQnwFvAw8OmuYDGzLmY2zswygJ3AVmqHySJZC+SZmf7/ltDol0+keUwGBrHnIbFdngBGBz93SQCup/ZkgFLgeOCqetb/Z/Bzg5nN3adqRfaThsVEmoGZ9QKWAV3dfUvY9Yg0Ne25iDSxYHjqeuAfChZpLZLCLkAkngXHSNYCq6g9DVmkVdCwmIiIRJ2GxUREJOo0LAZ06tTJ8/Lywi5DRKRFmTNnznp3z4n0nsIFyMvLo6CgIOwyRERaFDNbVd97GhYTEZGoU7iIiEjUKVxERCTqdMxFRESorKykqKiI8vLyb7yXlpZGbm4uycnJjd6ewkVERCgqKiIzM5O8vDzM7Kvl7s6GDRsoKiqiT58+jd5eqMNiZjbWzD42sxVmdkOE91PN7Kng/Q/MLK/Oe78Iln9sZqc0dpsiIvJN5eXlZGdnfy1YAMyM7OzsiHs0exJauASPcL0HOBUYCFxoZgN3a3Y5sNHd+wF3ArcF6w4ELgAOofaWGn81s8RGblNERCLYPVgaWr4nYQ6LDQdWuPtKADP7BzAOWFKnzTjgpuD1FOAvVtvLcdTeBHAn8KmZrQi2RyO2GTWzPyvlneUlYIYBZmAYZpBg//lC6i7fvV1QJ8Z/1tnVDjOSEozkxASSE42khNqfyYkJJCX+Z3lyYsLX3muTkkhaciJtkhNJTrR9+sUQEdkfYYZLD+CLOvNFwIj62rh7lZltBrKD5e/vtm6P4HVD2wTAzK4ErgTo1avXPnVg7qqN3P3Gin1at7kkJhhtkhNpk1IbNm2SE0lLSSQ9OZGM1ETapSWTmZZEuzbJtEtLpl2bJDLT/vO6XVoyWW1TyExNUkiJSKO12gP67j4JmASQn5+/T3fv/MHxB/CD4w/YtT3cwXe9BmqCZbXvg/PNNu5Anfdq6ix3nOoap6raqaiuoaraqayuobK6hqoap7KqhsoapypYVlntVFTVsLOqhu0VVZRXVrOjspodFTXsqKxiR0UwX1nDjooqVm+qZFl5GVt2VFK2s4o93cM0JTGB7LYptVNGKtltU+jUNpWsjBSyM1Lo0i6Nbu3T6NahDW1TW+2vlUiL5u4R/4jclxsch/mvwGqgZ5353GBZpDZFZpYEtAc2NLBuQ9tsEruGs4K55vjIqKqpcbZVVLGlvKo2bIKfm3dUUrqtgvXbdrJhawWl2yrYsHUnK9ZtZf3Wneys+uaTdjPTkmqDpn2br3726NiGvOx0emdn0KltivaCRGJMWloaGzZs+MZB/V1ni6Wlpe3V9sIMl9lAfzPrQ20AXABctFubacB44N/AucAb7u5mNg14wszuALoD/YEPqf1XvaFtSgQJCUZmWjKZacn06NCmUeu4O9srqlm/dSdrt+ykePMOijeXs2ZzOV9uqn29+MstrN+682vrZaQk0is7g7zsdHplp5OXnUFedgYDurQlu21qU3RPRBqQm5tLUVERJSUl33hv13UueyO0cAmOoVwDvAokAg+5+2IzuxkocPdpwIPAo8EB+1Jqw4Kg3dPUHqivAq5292qASNts7r61FmZGRmoSGalJ9M7OqLfdzqpqvtxUzmcbtrFq/TZWlW5n1YbtLF9bxsyl66io/s/eT3ZGCgO6ZDKgS1v6d8nkwK6ZDOicSfv0xl+8JSJ7Lzk5ea+uY2mIHhZG7TEX3RU5HNU1TvHmHaws2cbytWUUrt3Kx2vLKFxbxraK6q/adWmXyqHd2zMotz2DetROndvt3W66iESXmc1x9/xI7+nIq4QqMcHI7ZhObsd0jhvwn8dCuDtfbi5n+Zoylq8tY9maMhau3swbH6/76sSDLu1Sg6DpwNDeHTi8V0edTCASI/R/osQkM6NHhzb06NCGEw/q/NXybTurWFK8hYVFm1m4unaauaw2cBIMDu7WjvzeHRmWl0V+7450b+TxIxGJLg2LoWGxlq6svJJ5n2+iYNVG5qwqZd7nm9geDKl1b5/GUQd04pj+2Rx9QCcNpYlEkYbFJK5lpiVz3ICcr4bVqqprWFpcRsGqUmZ/VsrMZWt5Zm4RAP07t+Xofp04ul8nRvTNol2aThQQaQrac0F7LvGupsZZUryFd1es518r1jP7s1LKK2tITDDye3dk1MGdGXVwF/p2ytD1NyJ7YU97LgoXFC6tzc6qauau2sQ7hSW8sWwdy9aUAZCXnc5JB3Vh9MGdyc/LIiVJz9IT2ROFSwMULq3b6k07eGPpWmYuW8d7n2ygoqqGzNQkRg/swmmDunHsgE6kJiWGXaZIzFG4NEDhIrtsr6jiX4XrmbFkLa8tWcvmHZVkpiZx8iG1QXNMfwWNyC4KlwYoXCSSiqoa3vtkPdMXFPPq4jVsKa8iMy2JMQO7ctbhPRh5QDYJCTpGI62XwqUBChdpSEVVDe/WCZqy8iq6t0/j7KG5nDMslz6d6r/9jUi8Urg0QOEie6O8spoZS2pPb561vIQah2G9O3LusFxOO6ybTm+WVkPh0gCFi+yrtVvKeW7eaqbMKWLFuq2kJiVw+mHdufjIXgzp2UGnNktcU7g0QOEi+8vdWVC0macKvmDqvNVsq6jmkO7tuOTI3pwxpDvpKbpeWeKPwqUBCheJprLySp7/6Esef38Vy9aUkZmWxDlDc7l0ZB55OjYjcUTh0gCFizQFd2fOqo08+v4qXlpYTFWNc/LBXbj8mD4M75OlITNp8RQuDVC4SFNbV1bOo/9exWPvr2Lj9koG9WjPFcf24VuDupGcqDsBSMukcGmAwkWay46Kap6ZW8RD//qUleu30b19Glce15cLhvciLVkXZ0rLonBpgMJFmltNjfPGsnXcP+sTZn+2kU5tU7ji2L5cfGRvPfBMWgyFSwMULhKmD1Zu4C9vruCdwvW0b5PMZUfncenIPDqkp4RdmsgeKVwaoHCRWDD/i0385c0VzFiylszUJC4/tg+XH9OHTF2UKTFK4dIAhYvEkqXFW7jr9eW8ungtHdOT+eHxB/Ddo/Jok6JjMhJbFC4NULhILFpQtInbX1vOrOUldM5M5ZqT+nHBEb30nBmJGXsKF/2WisSow3I78Mj3hvP0D44iLzuD30xdzKg73mL6gmL0R6HEOoWLSIwb3ieLp35wJJO/N5yMlCSufmIu5973b+Z+vjHs0kTqpXARaQHMjOMH5DD92mO57ZxBfF66nbP/+h7XPDGXL0q3h12eyDeEEi5mlmVmM8ysMPjZsZ5244M2hWY2vs7yYWa20MxWmNndFtxHw8xuMrPVZvZRMH2rufok0hwSE4zzj+jFWz87gWtH9ef1pWsZ9ce3ue2VZWyvqAq7PJGvhLXncgMw0937AzOD+a8xsyxgAjACGA5MqBNC9wLfB/oH09g6q97p7kOC6aUm7INIaDJSk7j+5AG89bMTOX1wN+596xNG/fFtHY+RmBFWuIwDJgevJwNnRmhzCjDD3UvdfSMwAxhrZt2Adu7+vtf+X/RIPeuLxL2u7dO44ztDeOaqo+iYnsLVT8zlkgc/ZMW6rWGXJq1cWOHSxd2Lg9drgC4R2vQAvqgzXxQs6xG83n35LteY2QIze6i+4TYAM7vSzArMrKCkpGSfOiESK4b1zmLaNUfz2zMOYX7RJsbeNYv/e2kp23ZqqEzC0WThYmavm9miCNO4uu2CvY9o7cffCxwADAGKgT/W19DdJ7l7vrvn5+TkROnjRcKTlJjA+JF5vPmzEzjr8B7cP2slY+6cxZsfrwu7NGmFmixc3H20ux8aYZoKrA2Gtwh+RvrtXw30rDOfGyxbHbzefTnuvtbdq929BvgbtcdqRFqVTm1T+cN5g/nnD48iLTmByx6ezXX/mMeGrTvDLk1akbCGxaYBu87+Gg9MjdDmVWCMmXUMhrfGAK8Gw2lbzOzI4Cyx7+5af1dgBc4CFjVVB0Ri3RF5Wbz0k2O5dlR/pi8sZvQdb/Ps3CId8JdmEVa4TARONrNCYHQwj5nlm9kDAO5eCtwCzA6mm4NlAD8CHgBWAJ8ALwfLfx+corwAOBH4aTP1RyQmpSYlcv3JA5h+7bHkdcrg+qfn892HPqRoo66Nkaale4uhe4tJ61Bd4zz2/ip+/8oyzIzfnD6Q8/Jz9bhl2We6t5iIkJhgjB+ZxyvXHcch3dvxP88s4IrJBazbUh52aRKHFC4irUzPrHSe/P6R/Pr0gfxrxXrG3DWLFxd8GXZZEmcULiKtUEKCcfkxfZh+7bH0zs7gmifm8eMn57F5e2XYpUmcULiItGL9OrflmR8exX+fPICXFxbzrbvfoeCz0oZXFGmAwkWklUtKTODHo/oz5aqRJCUa37n/39z1+nKqqmvCLk1aMIWLiAAwpGcHXvzxMYwb0oO7Xi/kor99wOpNO8IuS1oohYuIfCUzLZk7zx/CHd8ZzOIvN3PqXbN4ZVFxwyuK7EbhIiLfcPbQ3K8uvPzhY3O5+YUlVGqYTPaCwkVEIsrrlMGUH47k0pF5PPTup1ww6X3WbNY1MdI4ChcRqVdKUgI3nXEIf77wcJYWb+G0u9/h3RXrwy5LWgCFi4g06NuDuzPtmqPpmJHCJQ9+wF/eKKSmRreOkvopXESkUfp1zmTq1Udz+mHduf215Xz/kQK2lOuiS4lM4SIijZaRmsSfLhjCb884hLeXl3DWPe+yskSPVJZvUriIyF4xq70B5mNXjGDj9krG3fMub+lpl7IbhYuI7JMj+2Yz9eqj6dGhDd/7+2wmzfpEDyKTryhcRGSf9cxK59kfjWTsoV3535eWcf3T8ymvrA67LIkBChcR2S/pKUncc9FQrj95AM/NW835k96npGxn2GVJyBQuIrLfzIxrR/XnvouH8fGaLZz113cpXFsWdlkSIoWLiETN2EO78tSVR1FeWcPZ976nCy5bMYWLiETV4J4deP7qkXRrn8b4hz7k6YIvwi5JQqBwEZGoy+2YzpSrRnJk32z+Z8oCbn/1Y51J1sooXESkSbRLS+bhy47g/Pye/OXNFVz/9HzdWbkVSQq7ABGJX8mJCUw8ZxA9s9pw+2vLKd1Wwb0XDyU9Rf/0xDvtuYhIkzIzrjmpPxPPHsQ7hSVc+LcPKN1WEXZZ0sRCCRczyzKzGWZWGPzsWE+78UGbQjMbX2f578zsCzPbulv7VDN7ysxWmNkHZpbXtD0Rkca6YHgv7rt4GMuKt3Dufe9RtHF72CVJEwprz+UGYKa79wdmBvNfY2ZZwARgBDAcmFAnhF4Ilu3ucmCju/cD7gRua4LaRWQfjTmkK49ePoKSsp2ce++/Wa5rYeJWWOEyDpgcvJ4MnBmhzSnADHcvdfeNwAxgLIC7v+/ukR7sXXe7U4BRZmZRrVxE9svwPlk8/YOjqHHn3HvfY86qjWGXJE0grHDpUicc1gBdIrTpAdQ9Qb4oWLYnX63j7lXAZiA7UkMzu9LMCsysoKSkZG9qF5H9dHC3djxz1UiygoePvfeJLraMN00WLmb2upktijCNq9vOa09+b/YT4N19krvnu3t+Tk5Oc3+8SKvXMyudp39wFD06tOGyh2fzpm7bH1eaLFzcfbS7HxphmgqsNbNuAMHPSL9Vq4GedeZzg2V78tU6ZpYEtAc27G9fRKRpdG6XxlM/OIp+ndty5SMFvLIo0mi3tERhDYtNA3ad/TUemBqhzavAGDPrGBzIHxMsa+x2zwXecF0WLBLTsjJSeOL7RzKoR3uufmIez89r6G9IaQnCCpeJwMlmVgiMDuYxs3wzewDA3UuBW4DZwXRzsAwz+72ZFQHpZlZkZjcF230QyDazFcD1RDgLTURiT/s2yTx6+QiG52Xx06c/4skPPw+7JNlPpj/sIT8/3wsKCsIuQ6TVK6+s5oePzeGtj0u4ZdwhXHJUXtglyR6Y2Rx3z4/0nq7QF5GYkZacyKRL8hl9cGd+PXUxj72/KuySZB8pXEQkpqQkJXDPfw1l1EGd+dXzi3j8AwVMS6RwEZGYk5qUyF8vHspJB3XmxucW8cQHOgbT0ihcRCQmpSYlcu/FQznxwBx++dxC/qGD/C2KwkVEYlZtwAzj+AE53PDsQp6arYBpKRQuIhLT0pITuf+SYRwXBMxz84rCLkkaQeEiIjGv9iyyYRzZJ5uf/XMBry5eE3ZJ0gCFi4i0CGnJifxtfD6DerTnx0/M451C3XA2lilcRKTFaJuaxOTLhtM3J4MrH5lDwWelYZck9VC4iEiL0j699lYx3dqncdnDs1m0enPYJUkEChcRaXFyMlN57IoRtGuTzCUPfkChnmgZcxQuItIide/QhsevGEFSYgIXP/gBRRu3h12S1KFwEZEWK69TBo9dPoIdFdV896EPKd1WEXZJElC4iEiLdmDXTB689AhWb9zB9/4+m+0VVWGXJChcRCQOHJGXxZ8vPJwFRZu4+vG5VFbXhF1Sq6dwEZG4MOaQrvzurEG8+XEJNzyzED2rKlxJYRcgIhItFw7vxbotO7nz9eXkZKZyw6kHhV1Sq6VwEZG4cu2ofqwrK+e+tz8hJzOVy4/pE3ZJrZLCRUTiiplx87hD2bC1glunL6FHhzaMPbRr2GW1OjrmIiJxJzHBuOuCIQzO7cB1T83joy82hV1Sq6NwEZG4lJacyAPj88nJTOWKybP5olQXWTYnhYuIxK1ObVN5+NIjqKiq4bK/z2bzjsqwS2o1FC4iEtf6dc7k/kvyWbVhG1c9NoeKKl0D0xwULiIS9446IJuJZx/Ge59s4JfP6RqY5lBvuJjZS2aW13yliIg0nXOG5fKTUf2ZMqeIe95cEXY5cW9Pey4PA6+Z2Y1mlhzNDzWzLDObYWaFwc+O9bQbH7QpNLPxdZb/zsy+MLOtu7W/1MxKzOyjYLoimnWLSMt23ej+nHV4D25/bTmvLNKjkptSveHi7v8EhgLtgAIz+5mZXb9r2s/PvQGY6e79gZnB/NeYWRYwARgBDAcm1AmhF4JlkTzl7kOC6YH9rFNE4oiZ8X9nD2JIzw5c//RHLPlyS9glxa2GjrlUANuAVCBzt2l/jAMmB68nA2dGaHMKMMPdS919IzADGAvg7u+7e/F+1iAirVBaciKTLhlGu7Rkvv9IAeu37gy7pLi0p2MuY4GPgHRgqLtPcPff7pr283O71AmHNUCXCG16AF/UmS8KljXkHDNbYGZTzKxnfY3M7EozKzCzgpKSkkYXLiItX+d2aUz67jDWb92pM8iayJ72XG4EznP3G9x9r68+MrPXzWxRhGlc3XZee9pGtE7deAHIc/fDqN3TmVxfQ3ef5O757p6fk5MTpY8XkZbisNwO3H7eYGZ/tpFfP79IZ5BFWb33FnP3Y/dnw+4+ur73zGytmXVz92Iz6wasi9BsNXBCnflc4K0GPnNDndkHgN83umARaXW+Pbg7y9eW8ec3VnBg10y+p5tcRk1Y17lMA3ad/TUemBqhzavAGDPrGBzIHxMsq1cQVLucASyNQq0iEsd+OnoApxzShVunL2HWcg2RR0tY4TIRONnMCoHRwTxmlm9mDwC4eylwCzA7mG4OlmFmvzezIiDdzIrM7KZgu9ea2WIzmw9cC1zajH0SkRYoIcG44ztDGNAlkx8/OU/3IIsS0zgj5Ofne0FBQdhliEiIVm3Yxrf//C96dEzn2atG0iYlMeySYp6ZzXH3/Ejv6fYvIiJA7+wM/nTh4Sxbs4VfPLtAB/j3k8JFRCRw4oGd+enoATz/0ZdMfu+zsMtp0RQuIiJ1XHNiP0Yf3IVbpy/lw09Lwy6nxVK4iIjUkZBg3HH+YHplpfOjx+eyZnN52CW1SAoXEZHdtEtL5v5LhrG9ooqrHp/DzqrqsEtqcRQuIiIR9O+Sye3nDWbe55v43XRdMre3FC4iIvX41qBuXHFMHx759ypemP9l2OW0KAoXEZE9+PmpBzG0VwdueGYBK0u2NryCAAoXEZE9Sk5M4C8XDSUlKYEfPT6X8kodf2kMhYuISAO6d2jDnecPYdmaMiZMXRx2OS2CwkVEpBFOOLAz15zYj6cKvmDKnKKwy4l5ChcRkUa6bnR/juybxa+eX8jHa8rCLiemKVxERBopKTGBuy84nLapyVz1+By27awKu6SYpXAREdkLndulcfeFQ/h0/TZ+o+Mv9VK4iIjspZEHdOLHJ/bjmblFPD9vddjlxCSFi4jIPrh2VH/ye3fkV88vYtWGbWGXE3MULiIi+yApMYG7LhhCgsG1T86joqom7JJiisJFRGQf5XZM57ZzDmN+0Wb+OOPjsMuJKQoXEZH9cOqgblw0ohf3v72SWctLwi4nZihcRET2029OH8iALm25/un5lJTtDLucmKBwERHZT2nJifz5wqGUlVfy3/+cT02Nh11S6BQuIiJRcGDXTH51+kBmLS9h8r8/C7uc0ClcRESi5OIRvRh1UGcmvryM5Wtb9+1hFC4iIlFiZkw85zDapiZx3T8+atWnJ4cSLmaWZWYzzKww+NmxnnbjgzaFZjY+WJZuZtPNbJmZLTaziXXap5rZU2a2wsw+MLO85umRiEitnMxUJp5zGEuKt3DHjOVhlxOasPZcbgBmunt/YGYw/zVmlgVMAEYAw4EJdULodnc/CDgcONrMTg2WXw5sdPd+wJ3AbU3bDRGRbzp5YBcuHN6T+2d9wgcrN4RdTijCCpdxwOTg9WTgzAhtTgFmuHupu28EZgBj3X27u78J4O4VwFwgN8J2pwCjzMyaqA8iIvX61WkD6Z2VzvVPz2dLeWXY5TS7sMKli7sXB6/XAF0itOkBfFFnvihY9hUz6wB8m9q9n6+t4+5VwGYgO1IBZnalmRWYWUFJiS58EpHoykhN4o7zh7BmSzk3tcK7JzdZuJjZ62a2KMI0rm47d3dgr08KN7Mk4Engbndfubfru/skd8939/ycnJy9XV1EpEFDe3XkmhP78ey81by44Muwy2lWSU21YXcfXd97ZrbWzLq5e7GZdQPWRWi2Gjihznwu8Fad+UlAobvftds6PYGiIHzaA61zwFNEYsI1J/XjreUl3PjcIobnZdG5XVrYJTWLsIbFpgHjg9fjgakR2rwKjDGzjsGB/DHBMszsVmqD47o9bPdc4I1gz0hEJBTJiQnc+Z3BlFdW88vnFtJa/kkKK1wmAiebWSEwOpjHzPLN7AEAdy8FbgFmB9PN7l5qZrnAjcBAYK6ZfWRmVwTbfRDINrMVwPVEOAtNRKS59c1py/875UBeX7qOZ+e2joeLWWtJ0T3Jz8/3goKCsMsQkThWXeNcMOnfLFtTxoyfHk/X9i1/eMzM5rh7fqT3dIW+iEgzSEww/nDuYCqra7jh2QVxPzymcBERaSZ5nTK4YexBvPVxCf8sKAq7nCalcBERaUbfPSqPEX2yuOXFJXy5aUfY5TQZhYuISDNKCIbHqt35+TPxOzymcBERaWa9stP5xakH8U7hev4x+4uGV2iBFC4iIiH4rxG9GXlANre+uISijdvDLifqFC4iIiFISDBuO+cwgLgcHlO4iIiEpGdWOr887WDeXbGBpwvia3hM4SIiEqILj+jF8D5Z/G76UtZtKQ+7nKhRuIiIhCghwZh49iDKq2qYMC1+bs2vcBERCVnfnLb8ZFR/Xl60hlcWrQm7nKhQuIiIxIArj+vLwd3a8Zupi9i8o+U/uVLhIiISA5ITE7jtnEGs37qTiS8vC7uc/aZwERGJEYflduCKY/vy5Ief8/7Klv2cQ4WLiEgM+enoAfTKSucXzy6kvLI67HL2mcJFRCSGtElJ5P/OHsSn67dx98zCsMvZZwoXEZEYc3S/Tpw3LJf7Z61k8Zebwy5nnyhcRERi0I2nHUzH9BRueGYhVdU1YZez1xQuIiIxqEN6CjedMZCFqzfz6Purwi5nrylcRERi1GmDunHcgBz++Npy1rawW8MoXEREYpSZccu4Q6ioruHmF5eEXc5eUbiIiMSw3tkZXHNiP6YvKObt5SVhl9NoChcRkRj3g+P70rdTBr+ZuqjFXPuicBERiXGpSYnceuahrNqwnb++uSLschpF4SIi0gKM7NeJM4d05763V/JJydawy2lQKOFiZllmNsPMCoOfHetpNz5oU2hm44Nl6WY23cyWmdliM5tYp/2lZlZiZh8F0xXN1ScRkaZ242kDSU1O4NfPL4r5xyKHtedyAzDT3fsDM4P5rzGzLGACMAIYDkyoE0K3u/tBwOHA0WZ2ap1Vn3L3IcH0QJP2QkSkGeVkpvI/Yw/ivU82MG3+l2GXs0dhhcs4YHLwejJwZoQ2pwAz3L3U3TcCM4Cx7r7d3d8EcPcKYC6Q2ww1i4iE7qLhvRjcswO3vLiEzdtj97kvYYVLF3cvDl6vAbpEaNMD+KLOfFGw7Ctm1gH4NrV7P7ucY2YLzGyKmQLgoVkAAAk0SURBVPWsrwAzu9LMCsysoKSk5ZzeJyKtW2KC8bszD6V0WwV/eC12n/vSZOFiZq+b2aII07i67bx24HCvBw/NLAl4Erjb3VcGi18A8tz9MGr3dCbXt767T3L3fHfPz8nJ2duPFxEJzaE92jN+ZB6Pf/A5C4o2hV1ORE0WLu4+2t0PjTBNBdaaWTeA4Oe6CJtYDdTd88gNlu0yCSh097vqfOYGd98ZzD4ADItmn0REYsVPTx5AdkYqE6YtpqYm9g7uhzUsNg0YH7weD0yN0OZVYIyZdQwO5I8JlmFmtwLtgevqrrArsAJnAEujXLeISExol5bMz8ceyLzPN/HsvNUNr9DMwgqXicDJZlYIjA7mMbN8M3sAwN1LgVuA2cF0s7uXmlkucCMwEJi72ynH1wanJ88HrgUubc5OiYg0p3OG5nJ4rw5MfHkZW8pj6+C+xfq50s0hPz/fCwoKwi5DRGSvLSjaxLh73uXyo/vwq9MHNutnm9kcd8+P9J6u0BcRacEOy+3ABUf05O/vfUbh2rKwy/mKwkVEpIX72ZgDSU9J5KYXFsfMlfsKFxGRFi67bSo/O+VA3l2xgVcWrQm7HEDhIiISFy4a3ouDumZy6/Sl7KgI/7b8ChcRkTiQlJjAb884hNWbdnDv25+EXY7CRUQkXozom80Zg7tz39uf8PmG7aHWonAREYkjv/zWwSQlGLdMXxJqHQoXEZE40rV9Gj8+qT8zlqzlX4XrQ6tD4SIiEme+d0wePbPacOv0JVSHdN8xhYuISJxJTUrkF6cezLI1ZTw1+4uGV2gCChcRkTh06qFdGZ6XxR9f+ziU+44pXERE4pCZ8evTB1K6vYJ73lzR7J+vcBERiVODcttz9uG5PPyvz5r91GSFi4hIHPufsQeSmGBMfKV5H2+lcBERiWNd2qXxw+MP4KWFa/jw09Jm+1yFi4hInLvyuL50a5/GLS8uabZHIitcRETiXJuURH4+9iAWrt7cbI9EVriIiLQCZwzuzuCeHfjDq8vYXlHV5J+ncBERaQUSEozfnH4wa7fs5L63Vzb95zX5J4iISEwY1juL0w/rxqRZn1C8eUeTfpbCRUSkFfn52IOoqYE/vra8ST9H4SIi0or0zErnsqPzeGZuEUu+3NJkn6NwERFpZX50Yj/at0nmf19ainvTnJqscBERaWXat0nm2pP6868V63l7eUmTfEZo4WJmWWY2w8wKg58d62k3PmhTaGbj6yx/xczmm9liM7vPzBL3ZrsiIq3ZxUf25sQDc0hJbJoYCHPP5QZgprv3B2YG819jZlnABGAEMByYUCcsvuPug4FDgRzgvMZuV0SktUtJSuDhy4Yzsl+nJtl+mOEyDpgcvJ4MnBmhzSnADHcvdfeNwAxgLIC77zoSlQSkALsGDhuzXRERaUJhhksXdy8OXq8BukRo0wOo+xi1omAZAGb2KrAOKAOm7MV2MbMrzazAzApKSppmzFFEpLVq0nAxs9fNbFGEaVzddl57usJen7Lg7qcA3YBU4KQI79e7XXef5O757p6fk5Oztx8tIiJ7kNSUG3f30fW9Z2ZrzaybuxebWTdq90B2txo4oc58LvDWbp9RbmZTqR0OmwE0ZrsiItKEwhwWmwbsOvtrPDA1QptXgTFm1jE4kD8GeNXM2gbBgZklAacBy/ZiuyIi0oTCDJeJwMlmVgiMDuYxs3wzewDA3UuBW4DZwXRzsCwDmGZmC4CPqN07uW9P2xURkeZjTXV1ZkuSn5/vBQUFYZchItKimNkcd8+P9J6u0BcRkajTngtgZiXAqn1cvROwPorlxLLW0lf1M/60lr42dz97u3vE020VLvvJzArq2y2MN62lr+pn/GktfY2lfmpYTEREok7hIiIiUadw2X+Twi6gGbWWvqqf8ae19DVm+qljLiIiEnXacxERkahTuIiISNQpXBrJzMaa2cdmtsLMIj3YLNXMngre/8DM8pq/yv3XiH5eamYlZvZRMF0RRp37y8weMrN1ZraonvfNzO4O/jssMLOhzV1jNDSinyeY2eY63+dvmrvGaDCznmb2ppktCZ5O+5MIbeLlO21MX8P/Xt1dUwMTkAh8AvSl9sFk84GBu7X5EXBf8PoC4Kmw626ifl4K/CXsWqPQ1+OAocCiet7/FvAyYMCRwAdh19xE/TwBeDHsOqPQz27A0OB1JrA8wu9uvHynjelr6N+r9lwaZziwwt1XunsF8A9qb/FfV90nYE4BRpmZNWON0dCYfsYFd58FlO6hyTjgEa/1PtBh1524W5JG9DMuuHuxu88NXpcBS6nzYMFAvHynjelr6BQujbPHJ2Lu3sbdq4DNQHazVBc9jeknwDnBsMIUM+vZPKU1u8b+t4gHR5nZfDN72cwOCbuY/RUMSR8OfLDbW3H3ne6hrxDy96pwkb31ApDn7odR+3C2yQ20l9g2l9r7Qw0G/gw8H3I9+8XM2gLPANe5+5aw62lKDfQ19O9V4dI4q4G6f6HnBssitgkeYNYe2NAs1UVPg/109w3uvjOYfQAY1ky1NbfGfOctnrtvcfetweuXgGQz6xRyWfvEzJKp/cf2cXd/NkKTuPlOG+prLHyvCpfGmQ30N7M+ZpZC7QH7abu1qfsEzHOBNzw4staCNNjP3caoz6B2vDceTQO+G5xhdCSw2d2Lwy4q2sys665jg2Y2nNp/E1raH0UEfXgQWOrud9TTLC6+08b0NRa+16Tm/LCWyt2rzOwaah+7nAg85O6LzexmoMDdp1H7ZT9qZiuoPYB6QXgV75tG9vNaMzsDqKK2n5eGVvB+MLMnqT2jppOZFQETgGQAd78PeInas4tWANuBy8KpdP80op/nAleZWRWwA7igBf5RBHA0cAmw0Mw+Cpb9EugF8fWd0ri+hv696vYvIiISdRoWExGRqFO4iIhI1ClcREQk6hQuIiISdQoXERGJOoWLSIwysw5m9qOw6xDZFwoXkdjVgdq7bYu0OAoXkdg1ETggeB7HH8IuRmRv6CJKkRgV3PH2RXc/NORSRPaa9lxERCTqFC4iIhJ1CheR2FVG7WNsRVochYtIjHL3DcC7ZrZIB/SlpdEBfRERiTrtuYiISNQpXEREJOoULiIiEnUKFxERiTqFi4iIRJ3CRUREok7hIiIiUff/ATKL4l4ByhrbAAAAAElFTkSuQmCC\n",
"text/plain": [
"