IntroNP.ipynb

intronp.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "mimetype": "text/x-python",
      "nbconvert_exporter": "python",
      "name": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.6-final",
      "file_extension": ".py",
      "codemirror_mode": {
        "version": 3,
        "name": "ipython"
      }
    },
    "colab": {
      "name": "IntroNP.ipynb",
      "provenance": [],
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/RichardSear/7cc9d153bfd8c91b03513452d41bb525/intronp.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vfz8XItJjKXt",
        "colab_type": "text"
      },
      "source": [
        "# <center>Introductory numerical physics & Numerically solving the equations for radioactive decay and for oscillatory motion\n",
        "\n",
        "\n",
        " ## <center>Part of Energy, Entropy and Numerical Physics (PHY2063) module\n",
        " \n",
        "## <center> 2020 notes by Richard Sear, with contributions from Jake Frazer, and evolved from earlier notes by Jim Al-Khalili and Jeff Tostevin"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X7qsuO7zjKXu",
        "colab_type": "text"
      },
      "source": [
        "#### Learning Objectives for this notebook:\n",
        "\n",
        "> Learn Python Jupyter notebook basics\n",
        "\n",
        "> Learn basics of plotting in Jupyter notebook (using pyplot)\n",
        "\n",
        "\n",
        "## Introduction to the Numerical Physics part of Energy, Entropy and Numerical Physics\n",
        "\n",
        "This numerical physics course is part of the second-year\n",
        "Energy, Entropy and Numerical Physics module (PHY2063),\n",
        "and is online at the EENP module on SurreyLearn.\n",
        "See there for assignments, deadlines etc.\n",
        "\n",
        "Very few models in physics can be solved analytically with just pen and paper.\n",
        "Therefore most predictions\n",
        "and modelling done in physics requires numerically solving the equations of the model.\n",
        "This is true for electromagnetism, quantum physics (including nuclear and particle physics),\n",
        "astrophysics, soft matter physics, etc.\n",
        "This course introduces numerical physics through\n",
        "some simple common examples, and discusses general\n",
        "points about numerical physics, in particular the key problem of estimating\n",
        "how accurate the numerical solution is.\n",
        "\n",
        "It is taught using Python, which is a good language for numerical physics;\n",
        "One of the main strengths of Python is its modularity and so it often makes it trivial to solve \n",
        "otherwise complex problems by simply importing a module that has a premade solution. Additionally\n",
        "it is widely used throughout scientific communities meaning there are many very useful and efficient\n",
        "scientific modules that you will become familiar with during this course. Additionally Python is\n",
        "widely used outside of academia and so it is a very useful skill to have regardless of what career\n",
        "path you intend to follow. However almost all high performance computing (e.g., solving large PDEs,\n",
        "quantum matrix problems, molecular dynamics, fluid mechanics, etc) is done using Fortran or C. So Python\n",
        "is by no means the best at every aspect of computing, however for the size of the problems you will be \n",
        "handling python will be more than sufficient. The numerical physics skills you will learn using Python \n",
        "will be transferable to number crunching applications in other computer languages. And the coding \n",
        "skills are transferable to a huge variety of other tasks people use computers to do.\n",
        "\n",
        "The course looks at numerically solving ODEs (ordinary differential equations) and\n",
        "PDEs (partial differential equations), and introduces the numerical modelling\n",
        "of random processes, using both \n",
        "a very widely used numerical\n",
        "technique called Monte Carlo (MC), and the use of Bayesian statistics.\n",
        "MC is used a lot in fields from statistical physics,\n",
        "to nuclear and particle physics. \n",
        "Bayesian statistics is used from particle physics to predictive text on mobile phones.\n",
        "A lot of modelling in physics is based on PDEs,\n",
        "from heat and mass flows in statistical physics, a lot of electromagnetism and\n",
        "quantum physics, etc. And the numerical techniques we will learn\n",
        "are also used extensively outside\n",
        "physics, for example to model share prices in stock markets.\n",
        "Increasing numbers of graduates in physics (and related subjects) are going\n",
        "into data science, which uses many models of random processes such as Monte Carlo modelling\n",
        "and Bayesian statistics.\n",
        "\n",
        "Examples of ODEs are the equations for radioactive decay and for a mass attached\n",
        "to a spring.\n",
        "PDEs are even more common in physics than ODEs, for example, Schrödinger's\n",
        "time-independent equation in three dimensions is a PDE, Maxwell's equations\n",
        "which govern electromagnetism (including light), and the wave equation,\n",
        "are all PDEs. In many cases (all except very simple cases) we cannot\n",
        "solve them analytically, and so can only solve them numerically.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T3dSM1Q7jKXu",
        "colab_type": "text"
      },
      "source": [
        "### Running a Jupyter notebook\n",
        "\n",
        "A Jupyter notebook (https://jupyter.org/) looks a bit like a webpage and consists of a set of cells, one after another. Each cells is (usually - I just cover the basics here) one of two basic types: 1) Python code, 2) Markdown.\n",
        "\n",
        "Some basic commands are:\n",
        "\n",
        "> Double-click on a markdown cell to see the actual markdown\n",
        "\n",
        "> Clicking the \"Run\" button above, or Ctrl-Enter, runs the current cell - cell with box around it\n",
        "\n",
        "> Clicking the fast-forward double arrow button above, starts the notebook running afresh and then runs every cell one after another, starting with the first cell. Starting afresh means that the imported modules are forgotten, as are the values of variables that have been set. Once you are happy that the code is OK, it is good to do this, to check that the code all runs, and that it does so without some variable having a value you set a while ago and have forgotten about.\n",
        "\n",
        "> 'Insert' above is how you add new cells. A new cell is code by default, use the pull down with the little arrow on the right that says what the current cell is ('Code', 'Markdown', etc) to change the current cell from a code cell to a Markdown cell, etc.\n",
        "\n",
        "More details [here](https://jupyter-notebook.readthedocs.io/en/stable/examples/Notebook/Notebook%20Basics.html).\n",
        "\n",
        "## Markdown and code cells\n",
        "\n",
        "\n",
        "1) Python code cells just contain the usual Python code, just as you would have in a .py file. Each cell can be run individually: just select the cell, then click the play symbol above. The notebook server will then run the cell, putting the output below, and move on, ready to run the next cell.\n",
        "\n",
        "2) Markdown code cells are for text, and you can just write in them. Formatting etc is done by using markdown ...\n",
        "\n",
        "See for example https://jupyter.org/ for more details on the Jupyter project, and [some markdown basics are here](https://jupyter-notebook.readthedocs.io/en/latest/examples/Notebook/Working%20With%20Markdown%20Cells.html), but you can just Google stuff, eg Google \"markdown headings\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UFZGCYjkjKXv",
        "colab_type": "text"
      },
      "source": [
        "This is a very simple markdown cell, which illustrates what markdown is, and describes the Python code below. The Python below prints out the sum of two numbers, $a$ and $b$."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "tags": [],
        "id": "1ilVUp7vjKXw",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "0058214f-2ab8-4ad7-84ad-9d3537dd846c"
      },
      "source": [
        "# this is a very simple Python cell\n",
        "a=5.0\n",
        "b=4.0\n",
        "c=a+b\n",
        "print(' the sum of a and b is c = ',c)"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            " the sum of a and b is c =  9.0\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "j7AVqOLVjKX0",
        "colab_type": "text"
      },
      "source": [
        "Above we used very simple variable names, $a$, $b$ and $c$. This is OK to demonstrate simple arithmetic, but single character names like this are mostly not suitable for  code.\n",
        "\n",
        "\n",
        "A good variable name should make it obvious and unambiguous what it is, e.g., the variable name $mass$ clearly tells you this is the mass, not the $x$-coordinate, or value of the spring constant, etc. With a good variable name,\n",
        "if you have started a program and added a variable, and then come back to it after a week, you quickly know what the variable is. If you just call a variable $a$, then a week later you will not know what it is and will have to spend time working out what it is, or you risk making a mistake and using the variable incorrectly. Good variable names (and comments) avoid this problem. Some single character names, like $x$, are OK, because of the convention that in three dimensions we refer to the axes as $x$, $y$ and $z$.\n",
        "\n",
        "Examples of good variable names:\n",
        "~~~~\n",
        "mass   time     x       energy      pi\n",
        "~~~~\n",
        "Examples of bad variable names:\n",
        "~~~~\n",
        "a    thingy   d2\n",
        "~~~~\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pY-ewTM_jKX1",
        "colab_type": "text"
      },
      "source": [
        "Below is a simple example of plotting using matplotlib, in particular the pyplot version of matplotlib. It is a just a simple example in which we plot a quadratic function of $x$. Note that essentially whenever you want to use matplotlib you will also want the line \"%matplotlib inline\" as well, as that says to show plots in the notebook."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "WrgyNnz9jKX1",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 307
        },
        "outputId": "643e85cb-1a51-4c56-c90b-76a2728c3058"
      },
      "source": [
        "# useful imports for arrays and plotting capabilities\n",
        "import numpy as np\n",
        "### this is typically how matplotlib - the most common plotting package in Python is imported\n",
        "# this line tells the notebook to show the plot inside the notebook\n",
        "%matplotlib inline\n",
        "# this imports pyplot\n",
        "from matplotlib import pyplot as plt\n",
        "\n",
        "# constants and initial boundary conditions\n",
        "n_xpts=100\n",
        "x=np.linspace(0,3.0,n_xpts)\n",
        "y=1.0+2.0*x-x**2\n",
        "# the line below is the basic command to plot two arrays (that must have the same lengths)\n",
        "# the first array is the values along the x axis, the second is the values along the y axis\n",
        "# i.e., the first point plotted is (x[0],y[0]), the second is (x[1],y[1]), and so on\n",
        "plt.plot(x,y,linewidth=6,label='$y=1+2x-x^2$')\n",
        "# these lines label the axes (here x and y are unitless so no units are needed)\n",
        "plt.xlabel(\"x\",fontsize=24)\n",
        "plt.ylabel(\"y\",fontsize=24)\n",
        "plt.xticks(fontsize=24)\n",
        "plt.yticks(fontsize=24)\n",
        "# show label in a legend\n",
        "plt.legend(fontsize=24)\n",
        "plt.show()"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEiCAYAAACGKbfRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXwU5f3A8c83CeEICVeCnCGA3CBXVFBEEW8UseJdq7Y/r2q9a7FWQatW1NZ6VtG2Hi2tVz1QPBEV5dBwKDeC3Gc4QkIucnx/f8xEVpjd7Ca7m032+3695jXZfZ555hlm2e/OzHOIqmKMMcbEi4S6roAxxhgTTRb4jDHGxBULfMYYY+KKBT5jjDFxxQKfMcaYuJJU1xUwh0pPT9esrKy6roYxxtQr8+fP36mqGdXls8AXg7KyssjJyanrahhjTL0iIuuDyWe3Oo0xxsSVuAl8IpIpIjeJyDQR2SAipSJSICLfisiDItK+luW3E5HHRGSNiJSIyHZ3X6PDdQzGGGNqLy5udYpIZ2AdID5v5wMpwBHucpWInKuqM2tQ/hHAp0Abn7LTgTOBMSLye1V9sOZHYIwxJlzi5Yov0V2/B5wHtFbVFkAz4AxgLdAKeEtE2oVSsIg0Bd7BCXoLgf5u2a2AP+ME2wdE5JRwHIgxxpjaiZfAtwcYrKpnqurrqroHQFX3q+r7OMGvBEgDrg6x7KuBLsA+4CxVXeqWna+qtwFv4QS/P4XnUIwxxtRGXNzqVNW9wLcB0leIyFzgBGBoiMVf4q6nqupmj/SHgXHAEBHppaorQyzfxBBVZXfhfrbnl7K9oISdBaXkFZWxp2g/+SVlFJVWULS/guKyCioqlfLKSiorISEBkhISSEwQmjZKpFlyIs0aJ5LWpBGtmiXTslkj0lMbc1hqEw5La0zrlGREpPoKGWNCFheBL0i73HViwFw+RCSVA4HyQz/Z5gJ7gRbAaMACXz2wp3A/3+/Yx/c7Cli9Yx8bdhWxfncRG3cXUVpeGfH9N05KILN1M2dp04webVM5vG1zerRtTquU5Ijv35iGzAIfICJJwLHuyyUhbNqHAw1mlnplUNVKEVkJHAX0rXElTcTsLSpj4cY9LNqYx5LN+Szdspete0vqtE6l5ZVu4N13SFqHFk3o26EF/TumMahzSwZ3bkWLZo3qoJbG1E8W+BzXAe2ASuDFELbz7QKxJUC+qjS/XSZE5CrgKoDMzMwQqmBCtSO/hLlrdzNnzS6+XruLNbmFdV2lkGzZW8KWvSV8snz7j+91z0jh6G5tGNatDcO6tqZtWpM6rKExsS3uA5/bFaGq4cmTqroshM1TfP4uDpCvyF0395dBVacAUwCys7NtduAwKi2v4Ju1e/h81Q4+X5XLqu2HXkXVd2tyC1mTW8jUeRsA6HVYKiN7pnN8z7Yc2bUVjZOCvoNvTIMX14HP7bT+FtAUmA/8rm5rZMJlb3EZn63cwUdLt/PZyh0U7q+o6ypF1crtBazcXsBzs9bSvHESx/fK4JS+hzGqd1vSmthtURPf4jbwiUhr4COgK/A9MEZVQ32w43uPrClQ4CdfM3fd8C41YkhBSRmfLN/Oe99t5YtVO9lfEZlGKE0bJdKuhdP6sm1qE1qnOK0yWzZtRLPGSaQkJ9E0OYGkhASSEoSEBKGyUil3W3kW76+kcH85RaXl5BWXkVdUxu7C/ewoKGF7filb9xZTUha+uu8rLee977by3ndbSU5MYGTPdM48ogMn9T2M5o3j9ivAxLG4/NSLSAucVpj9gQ3ASaq6PfBWnnyf63XAf4vNDu56aw32YQKoqFS+XL2TN+Zv4sOl28LW4lIEstqk0POw5hzetjndM5rTpU0KXdo0o02EuxqoKjv37WfD7iLW7ypkTe4+vt/uNHRZt6sQrcWN8P0VlXyyfAefLN9Bk0YJnNqvHecO6cSxh6eTmGDdJ0x8iLvAJyIpwHQgG9iGE/Q21LC4FYDitOzsh0fgE5EEoJf7MpTnhyaATXuKePWbjbyas4lt+bVvgdm5dVOGZLZiUOeWDOjYgj7t00ipo6shESEjtTEZqY0Z2qXVT9L2lZazYms+izfvZdHGPBZs2MPG3YEeL/tXUlbJ24u28PaiLbRLa8L52Z04/8jOdGrVrPqNjanHRGvz87GecYcXew8YhdNv7/iqkVZqUebXwJHAM6p6rUf6cGC2+7J3MB3Ys7Oz1aYlOlRlpfLZqh28OHs9X3yfW6srn67pKQzv3obh3dpwdLfWtE2tv60gd+SXMG/tbub8sIs5a3axdmfNW6mKwPE9M7hseBbH98wgwa4CTT0iIvNVNbvafPES+EQkGXgbOA3IA0ar6oIwlHsT8CjO871eqrr1oPQ3gJ8BQZ0QsMB3sH2l5bz6zUZemrOOdbuKqs3vpXFSAsf1SOeEXm05vmcGnVs33KuaDbuK+Pz7XD5fuYNZ3++s8e3frDbNuOyYLM7L7mzPAk29YIHPh4gkAq8A5+IEqFNUdW6Q22bhDGINcIWqvnBQelNgOc54nQuAS1V1mTuqy13Ab92sp6rqR8Hs0wKfY0d+Cf/4ah3/nreegpLykLdPbZLEyX0P49R+7TiuRzrNkuPvy7tofzmzvt/Jh0u28fGy7RSUhv7vmNYkiUuGdeGKY7Ksf6CJaRb4fIjISOBz92UJzhBi/mxU1SN9ts0iQOBz8wwEZvDTaYma4wwCrkBI0xLFe+DbsKuIv32+hjfmbwq5ZWaTRgmc0rcdYwd24Lie6dZ/zUdpeQVfrNrJO99u4aMaNARKTkzg3KGd+PUJ3Rv0FbOpv4INfPHyE9h3Foom7uJPyC0lVPVbEekP3IEzB19HnGeIXwOPquqMUMuMR2t3FvLkp6t5a9FmKipD+0F2VFZrxg/txOkD2pFq/dQ8NU5K5OS+h3Fy38MoKClj+uKtvD5/E9+s2xPU9vsrKvnP1xt4NWcj4wZ15DcnHk5Wekr1GxoTY+Liiq++ibcrvk17inhixmpeX7AppIDXslkjxg/pxIVHZXJ4W7+D4phqrN5RwH++3sgbCzaRV1QW9HaJCcJ5Qzvxm9E96NiyaQRraExw7FZnPRYvgW/XvlKe+HQ1/563nrKK4D+HfdqnccWxWYwd2IEmjexWZriUlFXwzqIt/HP2OpZvzQ96u+TEBC4+OpMbRvegtc0cYeqQBb56rKEHvuL9Ffzjq7X87bM17AuhscXo3m25cmQ3ju7a2uaqiyBVZe4Pu3lu1g98umJH0NulNk7imhO688tju9I02X6QmOizwFePNdTAp6q88+0WHnx/RdDT/iQlCOMGd+Tqkd3ocVhqhGtoDrZqewFTvviBtxZupjzI29DtWzRhwum9GTuwg/1AMVFlga8ea4iB79uNedwzbSkLNuQFlT85MYHzsjtxzfHWgjAWbNzttLR9PSf4lrZDu7Ri4ll9OaJTywjXzhiHBb56rCEFvj2F+3now5X895sNQY20kpQgnH9kZ64fdTgdrMFEzNmSV8wTn67mtZyNQV0BisBFR2Vy+6m9aNnMnv+ZyLLAdxC3Q/konOHFst11Vb+7Pqq6ooblngDMDCJrhqruDKbMhhD4VJXXcjbxp/eXsyeIloIiMG5QR246qQdd2lgT+Vi3flchj368ircWBZp/+YDWKclMOL035w3tZLc/TcRY4DuIiIwD3vSTHI7AVwnkBsjaV1V3B1NmfQ98a3L38fv/LWbe2qAOl+N6pDPh9N7069AiwjUz4bZk814mf7CCWd8H9ZuOYd1a88A5A+iWYd1PTPhZB3ZvO4Ac4BtgM+6M52GyUVWzwlhevVNWUckzn63hiZmr2R/EqCA92jbnD2f25fieGVGonYmE/h1b8PKvjubzVbn88d1lrN4ReMrJuT/s5rTHZnHDiYdz9fHdaZSYEDC/MZEQT1d8iapa4fM6iwNDkYXjim99uAJffbziW7plL7e99l1Q/b9aNG3ELSf35JKjM0myL74Go6yikn/PXc+jn3zP3uLqb2/3bZ/GI+cNpG+HtCjUzsSDYK/44uZbxzfomfDZX17Jox+v4uwnvwoq6F10VGdm3nYClx2TZUGvgWmUmMDlx3Zl5m0ncEF252rzL9uaz9gnv+TRj1dRFuKYrMbURrzd6jRh9P32Am5+dRFLNlcf8Pp1SOO+cf0ZnNmq2rymfmudkszk8UdwwVGd+cObS1gW4AdReaXy2Izv+XTFDh69YJANPWeiwn5yh0+GiCwQkUJ3WSUiU0RkQF1XLNwqK5V/fLmWMU98WW3Qa9ookTvP6MPb1x1rQS/ODMlsxTvXH8vvz+hNk0aBv2oWb97LmMdn8cJXa6kMcYByY0JlgS98mgGDgVKcK+kewJXAQhG5rS4rFk47Ckq47J9fc++7y6ptwHJcj3Q+unkkV47sZrc141RSYgJXjezOxzcfz3E90gPmLS2vZNK0ZVzxwjfkFpRGqYYmHtm3Ue3lAQ/j9A1sqqqtcYLg8cBsIBF4WEQuDlSIiFwlIjkikpObG6hXRN2ZuXIHZzw2q9qm66mNk3jo3CN46ZdH2agrBoDOrZvx0i+PYvK5A0itZjb3z1flcvpjs/h8VWz+PzD1X9y06jxYuFp1VrOPZOBT4FhgE9BFVat9ih9rrTrLKip5+MOVTPnih2rzHt8zgwfPHUD7FjbqivG2Ja+YCf9bzBdBBLarR3bjtlN7WbcHExRr1RkDVHU/cJf7shPOrdB6ZUteMRc8O6faoNe0USL3n9OfF6440oKeCahDy6a8eMWR/HFc/2qf/T37xQ9cOGUuW/KKo1Q7Ew8s8EXePJ+/u9VZLWpg5sodnPH4rGoHlh7UuSXTbzyOS47uYsNRmaCICJcO68L0G45jYOfAg1jPX7+HMY/P4rOVwU+RZEwgFvjMISorlb9+sopfvvBNwBm5ReD6UYfz2jXD6Zpu42ua0HXLaM7r1wzn1yd0J9Bvpj1FZVzxwjc8PuN7a/Vpas368UXe0T5/r/WbK0bsLSrjplcWMnNl4Ocv7Vs04dELBjGsW5uA+YypTqPEBG4/rTcjDk/n5lcXsT3fu0WnKvzl41Us2pjHo+cPokWzRlGuqWko7IqvliTAvT0RaQTc677cCiyISqVqaNX2AsY+9WW1QW9Urwym33CcBT0TVsccns77N46sduzWT1fs4OynvuT77QVRqplpaOIq8IlIetUC+PambumbJiIJB22n7jLJo9glIvIbEelRFQRFJFFERgAzgBFuvjuCadFZVz5auo1znvqK9buK/OZJTBAmnN6bv192JK1SbG41E36tU5L55+VHcvtpvUhM8H/vc92uIs55ejYfL9sexdqZhiLebnX6u5SZc9DrrsC6IMvsCzzu/l0qIgVAGlAVGcqBP6jqiyHUM2pUlSc+Xc1fPl4VMF9GamOeungIR3VtHaWamXiVkCD8+oTDGZrZiuumLmTnPu9bn/tKy7nypRxuPbkn1594uDWsMkGLqyu+CLkaeAlYCuQDLXFGb1kMPAkMVNXJdVc9/0rKKrjxv4uqDXpHZrXivd+MsKBnourobm1474YRZHcJPNTdnz9exc2vLKKkzMahN8GJ2w7ssSwaHdhzC0q56uUcFlbTVeGKY7P4/Rl9rAOxqTNlFZXc/95yXpi9LmC+IZktmfKLbNKbN45OxUzMsQ7sxq9V2wsY99RXAYNeclICj5w3kIln9bOgZ+pUo8QEJo3tx8PjjyA5wGdxwYY8zn7yK2v0Yqpl32hxZvaanZz7t9lsDjASxmFpjXn16uGMH9opijUzJrDzsjvz36uH0TbV/xXd5rxizv3bbOas2RXFmpn6xgJfHHlz4SYu+8fXFJSU+81zRKcWvHP9CAZVM5qGMXVhSGYrpv1mBAM6tvCbJ7+knF/8Yx5vLdwcxZqZ+sQCXxxQVf722RpufuVbyir8P9M9Y0A7XrlqOIelNYli7YwJzWFpTXj16uGc3r+d3zxlFcpNryzi2c/XYO0YzMEs8DVwlZXKve8uY/IHgSefuH7U4Tx50RCaJidGqWbG1FzT5ESeungI143qHjDfn95fwX3vLbdhzsxPWOBrwPaXV3LjK4v451fr/OZJShAeGn8Et53ai4QAHYaNiTUJCcJvT+3N5HMHBOzs/vcv13LTK4uqnTjZxI+4CXwikioiY0XkjyLyvojs9BmRpXcYyk8TkftEZLmIFInILhGZISLjw1H/UBXtL+f/Xsph2rdb/OZJbZzEC1ccxfnZnaNYM2PC64IjM/nn5UfSPMAEt+98u4UrX8qheL/19TNxFPiA0cDbwB+A04CwDTQpIp2ARcCdQG+gAmf0lhOB10Tk6XDtKxh7i8r4+fPzAk70eVhaY169ZjgjeqRHsWbGRMbInhm8evXwgC0+P1+Vy6V/n8feYv8zjpj4EE+BD2AHMB24B7gqHAW643O+zoFhzo5V1VQgFbgdqASuFZErw7G/6uQWlHLBlDkB59DrnpHCG9ceQ5/2adGokjFR0bdDGv/79TF0y/A/RVbO+j1cOGUuuQXew6CZ+BBPgW+aqh6mqmNUdRLwcZjKPRtn6qFK4BxVnQ2gqiWq+jAHxvG8V0QiOrLz1r3ObOkrtvnvwDsksyWvX3MMnVo1i2RVjKkTnVo14/VrjgnYHWf51nwumDKHrXttVvd4FTeBT1UjdXP/Enf9iaou8kh/BFCgHc6tz4jYuLuI856Zww87C/3mOaFXBv/+v2E2s4Jp0FqnJDP1yqMZGWB6ox9yCzn/2Tls3O1/NhLTcMVN4IugUe76Q69EVd2MM4A1RCjwrcndx3nPzGHTHv+/YM88oj1TLs227gomLjRLTuL5X2QzZkB7v3k27i7m/GfnsDbAj0XTMFngqwURacuBRjJLA2Rd5q77hrsOK7blc8Gzc9iWX+I3z8VHZ/LYhYNJTrLTbeJHclICj180mIuO8t9qeeveEi54dg5rcvdFsWamrtk3Ye34/pz032/gQJr/n581ULS/nEv//jU79+33m+eqkd24f1z/gP2cjGmoEhOEB84ZwP+N6Oo3z46CUi6cMpfVOyz4xQsLfLXj23ws0JPyqgcJzf1lEJGrRCRHRHJyc/13Q/DVLDmJe8b2w19Mu3F0D+44vbdN0Gnimohw55g+3HDi4X7z5P4Y/Gxmh3hggS9GqOoUVc1W1eyMDP8P5Q92xoD2/Pn8gRwc2yac3pubT+5pQc8YnOB3yym9uP20Xn7z7NxXykXPzbPbnnHAAl/t+D4VbxogX1XfgYj8jzpncCfuHzfgx9f3jO3HNccHHsPQmHj06xMO584z+vhNzy0o5eLn5rLOGrw0aP7H+DHB8H2u1wFY7CdfB3e9NVIVufjoTErLK2iWnMgFR2ZGajfG1HtXjuyGCNz33nLP9O35pVz03FxeuWo4mW2sv2tDZFd8taCqucBO92W/AFmrWnMuC5Cn1q44tqsFPWOC8H/HdWPiWf4bWW/dW8Ilf59rndwbKAt8tTfTXZ/slSgiHTkQFGdEpUbGmGpdcWxX7j7Tf/DbuLuYS56fx859NrxZQ2OBr/amuutTRGSgR/otgODc5pzpkW6MqSO/HNGV35/hf3KWH3IL+fnz88gr8t9lyNQ/cRX4RCS9agFa+SS19E0TkYSDtquavmiSR7FvA/Nw/i3fFJFh7jaNReRW4CY330RVtf89xsSYq0Z257en+m/tuWJbAZf/8xsKS8ujWCsTSXEV+IBcn2WBz/tzDkoL+kGZqiowHliLM0PDHBEpwGnB+QjOv/EzqvpcOA7AGBN+1406nBtG9/CbvmhjHtf8az6l5TafX0MQb4EvIlR1EzAIeABYgdNatgDn1ub5qnptHVbPGBOEm0/qEXCEl1nf7+TmVxZRUalRrJWJhLjqzqCqNerNHcx2qpqPMxHtnTXZhzGmblWN8FJUVsHUeRs880xfvI0WTRfzwDkDbHCIesyu+IwxxiUi3Hd2f84e1MFvnv98vZFHP/k+irUy4WaBzxhjfCQkCI+cN5ATe7f1m+fxGd/z8tz1UayVCScLfMYYc5BGiQk8dfEQjspq7TfP3W8vYfriiA3GZCLIAp8xxnhompzI85dn06d9mme6Ktz0yiK+Wbc7yjUztWWBzxhj/Ehr0ogXrziSTq28x6DfX17J/72YY3P51TMW+IwxJoC2aU146ZdH0Tol2TN9b3EZl//za3YUlES5Zqam4i7wiUg7EXlMRNaISImIbBeRaSIyuoblneAzskugJT3cx2KMiY5uGc355+VH0iw50TN9055ifvVCDkX7bXSX+iCuAp+IHAEsAW4AugGlQDpwJvCxiEyoRfGVwPYAS2UtyjbG1LGBnVvy1MVDSEzw7r+3ePNebvyvdXCvD+Im8IlIU+AdoA2wEOivqi1wxuz8M85A0g+IyCk13MVGVW0XYLEn4MbUc6N6t+W+cf39pn+8bDt/mu49z5+JHXET+ICrgS44Y2iepapLwRlxRVVvA97CCX5/qrsqGmNi3UVHZXL9qMP9pj//5VpenrMuavUxoYunwHeJu56qqps90h9210NExP9Q7caYuHfrKT0ZF2B0l0nTlvHFqtwo1siEIi4Cn4ikAkPdlx/6yTYX2Ov+XaOGLsaY+CAiTB5/hN8O7hWVynVTF1g3hxgVF4EP6INzGxNgqVcGVa0EVrov/U/L7F+GiCwQkUJ3WSUiU0RkQA3KMsbEuMZJiTx76VC6pqd4pheUlPN/L35jk9jGoHgJfO19/t4SIF9VWvsAefxpBgzGaSmaBPQArgQWisht1W0sIleJSI6I5OTm2i0SY+qDVinJ/OPyI2nZrJFn+rpdRVz7rwWUVVij7lgSL4HP9ydZcYB8Re66eQhl5+E8H8wGmqpqa5wgeDwwG0gEHhaRiwMVoqpTVDVbVbMzMjJC2L0xpi51TU/hmZ8PJclPN4c5P+zi3mnLolwrE0i8BL6IUdVFqnq7qs5X1RL3vQpV/QIYBXzlZp0sIvbvbUwDNKxbG+4/x383h5fnrvc7x5+Jvnj5Ii70+dt70D1HM3cdlifSqrofuMt92QnnVqgxpgG64MjMgDO43/32Er5ea915Y0G8BD7f53r+2yAfSAvnXCPzfP7uFsZyjTEx5o4z+jCql/ejivJK5dp/zWdzXqCnLSYa4iXwrQCqxhHq55XBvQ1Z1X/PbsgbY0KWmCA8dtFgumV4t/TcVbifa16eT0lZRZRrZnzFReBT1QIgx315sp9sRwMt3L9nhHH3R/v8vTaM5RpjYlBak0Y8/4tsUpskeaYv3ryXO99cgqqN6VlX4iLwuaa660tExKu7QlWXg/mqutIj3ZOIeDflctIaAfe6L7cCC4It1xhTf3XLaM4TFw3GT0NP3liwiZfnro9upcyP4inwPQusB1KBd0WkLzijuojIQ8DP3Hy/P3hDn6mFJnmUu0REfiMiPaqCoIgkisgInCvHEW6+O9xO8saYOHBCr7bcflpvv+n3Tltms7fXkbgJfKpaDJwN7AKGAEtFZC9OP7zf4jwDvENVPwqx6L7A48AqoFhEcnH6A84CjgPKgQmq+mJYDsQYU29cPbIbY47wHg+jvFK57t8LbALbOhA3gQ9AVb8F+uMEqh+AxjiB8D3gZFV9sAbFXg28hDMUWj7QEmf0lsXAk8BAVZ1c+9obY+obEeHh8UfQu12qZ/qOglKun7rQRnaJMrEHrLEnOztbc3Jyqs9ojKkX1u8q5KwnviS/xHuG9iuP68qdY2oyRLDxJSLzVTW7unxxdcVnjDF1oUubFB67aDD+msI9N2st0xeHs/uwCcQCnzHGRMGoXm254cQeftNvf/071u4s9JtuwscCnzHGRMmNo3twfE/vkV32lZZz7b+sc3s0WOAzxpgoSUgQHrtwEJ1aeQ8ZvGJbARPf9pwy1ISRBT5jjImils2SefqSISQnen/9vpKzkTfmb4pyreJL3AU+EWknIo+JyBoRKRGR7SIyTURG17LcNBG5T0SWi0iRiOwSkRkiMj5cdTfGNAxHdGrJXWf28Zv+h7eWsHpHWCaJMR7iKvCJyBHAEuAGnJkSSoF04EzgYxGZUMNyOwGLgDuB3kAFkAacCLwmIk/XvvbGmIbk58O6cNZA78liissquH7qAnveFyFxE/hEpCnwDtAGWAj0V9UWQCvgz4AAD4jIKSGWK8DrQFdgHXCsqqbiDI12O1AJXCsiV4bpUIwxDYCI8KefDfA7k8OKbQXcYzO3R0TIgU9ETopERaLgaqALziSzZ6nqUgBVzVfV24C3cILfn0Is92ycGRgqgXNUdbZbbomqPowzSgzAvSKSXPvDMMY0FM0bJ/H0JUNonOT9VfyfrzfwzrdbPNNMzdXkiu8jEflBRCaKSJew1yhyLnHXU1V1s0f6w+56iIj08kivrtxPVHWRR/ojOOOAtsO59WmMMT/q3S6NSWM9pwkF4M7/LWbDrqIo1qjhq0ngKwKygLuBNSLysYhcKCKNw1qzMBKRVGCo+/JDP9nmAnvdv0Np6DIqULlukK1qn2yBzxhziAuP7MxYP8/7CkrLueG/Np5nONUk8B0GXAnMcbcfDfwb2CoiT4rI0EAb15E+OLcx4UAQ+gl3yqCqefiCGjRPRNriPDP0W66r6ka9DcZnjDmEiHD/Of3JatPMM33Rxjz+8vGqKNeq4Qo58Klqoar+XVVHAL2AyTiTrLYErgW+FpFv3TnqWoe3ujXmOy9IoBvmVWne84hEr1xjTJxJbdKIJy8eQqNE7wE9n/l8DV9+vzPKtWqYatWqU1W/V9U7gEycLgFvAmXAAOCvwGYReUVETgs0U3kU+DabKg6Qr+pGevNolysiV4lIjojk5ObmBrl7Y0xD0r9jC37nZ/JaVbjl1UXsLtwf5Vo1PGHpzqCqlao6XVXHAx2Bm3Bu7zUGxuPMd7febRDTLhz7bGhUdYqqZqtqdkaG91h+xpiG71cjujKql/d3wI6CUn73xnfYdHK1E4l+fFk4t0A74LRmFHfphNMg5gcRmRiB/QbiO+S59yB5jqob7MEOmRCpco0xcUpEePi8gWSkercX/HjZdqZ+vSHKtWpYwhL4RCRDRG4Wke+Ar4FrcDqGfwtcjxMELwVmA02Au0Xk9nDsO0i+z9+8m079NC3YibEiVa4xJo6lNzwM02EAACAASURBVG/Mo+cP8pv+x3eX2ZBmtVDjwCciCSJyloi8CWzC6a/WHygApgDZqjpEVZ9W1W2q+m9VPQ6nRagAV4Wh/sFagXP1CeDZYUZEEnCuVOFAK8yAVDUXqHra7L8jzoHWnDYMgzEmKCN6pHPlcV0900rKKrnplYXsL7cuDjVRk5Fb+orIw8BmnNFOzgYa4XRvuAJor6rXquoCr+1V9e/AbpxRVKJCVQuAHPflyX6yHQ20cP+eEULxMwOVKyIdORAUQynXGBPnbju1F33bp3mmLdmcz2MzrItDTdTkim8JcAtOf75dwF+Avqo6QlVfVNVArRur7Kvhvmtjqru+RES8uhXc5q7nq+pKj/Tqyj1FRAZ6pN+Cc4W7lQNB0hhjqtU4KZHHLxpEk0beX5d/+2wNOet2R7lW9V9Ng88nwAVAR1W9TVVXhLj9sTizI0TTs8B6nMGj3xWRvuCM6iIiDwE/c/P9/uANRUTdZZJHuW8D83D+Ld8UkWHuNo1F5FacFq4AE1XV2iEbY0JyeNtU/jDGe+yLSoWbX11EQUlZlGtVv9Uk8HVV1VNV9TVVrdG/tqpuVtX1Ndm2ptwr0bNxrlKHAEtFZC+QB/wW5xngHar6UYjlKk6XjbU4MzTMEZECnKvaR3D+jZ9R1efCdSzGmPhyydGZnNi7rWfaxt3F/PFdaz4QipqM3BLVgBVOqvotTgOcx4EfcPoZ7sLpZ3iyqj5Yw3I3AYOAB3Aa0iThNPKZCZyvqtfWvvbGmHglIjx47gBap3hP8PJqziY+WbY9yrWqv8Q6Qsae7OxszcnJqT6jMSaufLh0G1e/PN8zLb15Yz66eaTf4BgPRGS+qmZXly9uJqI1xpj67tR+7Tg/u5Nn2s59pfzhrcU2qksQLPAZY0w9cvdZ/ejUynugqOmLt9nEtUGwwGeMMfVI88ZJPHLeQPwN+3/XW0vYnl8S3UrVMxb4jDGmnhnWrQ2/OtZ7VJf8knJ+/z+75RmIBT5jjKmHbju1F4e39Z7pbMaKHfxvweYo16j+sMBnjDH1UJNGifzl/IEkJnjf85w0bSnb9totTy9xFfhEJE1E7hOR5SJSJCK7RGSGiIyvRZlZPiO7BFqqbWJrjDGhOKJTS359QnfPtIKScib8z+bu8xI3gU9EOgGLgDuB3kAFkAacCLwmIk+HYTfbAyw2ppAxJux+c2IPerdL9Uz7bGWu3fL0EBeBT0QEeB1nSLF1wLGqmoozbuftQCVwrYhcWZv9qGq7AMu3tTwMY4w5RHJSAo+cN5AkP7c875m2lB3WyvMn4iLw4YzReTROgDtHVWcDqGqJqj6MM4QZwL0iEr/DHhhj6qX+HVtw3ajDPdPyS8q5860ldsvTR7wEvkvc9Sequsgj/RGcQarb4dz6NMaYeuW6UYf7veX58bLtTPtua5RrFLviJfCNctcfeiWq6mZgqfvSAp8xpt5JTkrg4fH+W3lOfHsJu/aVRrlWsanBBz4RaQu0cV8uDZC1al4P74mvgtvXHBHJF5FiEVkrIv8SkRE1Lc8YY0IxoFMLrh7pPdXpnqIy7rXpi4A4CHyA72zrgQaxq0rzmp09WMNwniMCZOHcYp0lIn91G9gYY0xE3TC6h9+O7W8v2sKnK2z6ongIfCk+fxcHyFfkrr0/Mf6VAE8DI4FUVW0JNAOGAtPcPDcCdwQqRESuEpEcEcnJzc0NsQrGGONo0iiRh8Yf4XcszzvfXBL3M7bHbOATkbtFpLyGy/3RqqeqblPV61R1lqruc99TVV2gqmOB19ysvxeRlgHKmaKq2aqanZGREY2qG2MaqCGZrbj8mCzPtK17S3jog5XRrVCMidnAh1O3xFosVQp9/vaey8PRzF3vC0Pdff3OXacAo8NctjHGeLrtlF50bOn9lffy3PV8s253lGsUO2I28KnqJFWVGi4TfIryfa7XIcAuq9LC2uZXVdcCVfcuvZ86G2NMmKU0TuJPPxvgN/2O/y2mtLwiijWKHTEb+MJFVXOBne7LfgGyVrXmtGZPxpgGYWTPDH42pKNn2uod+3j28x+iXKPY0OADn2umuz7ZK1FEOnIgKM4I545FpCtQ9dBubTjLNsaY6tw1pi9tUrwHpHry09WsyQ33053YFy+Bb6q7PkVEBnqk3wIIzm3OmR7pfgXRTeEBd10MfBpK2cYYU1utUpK5+yzv7sn7KyrjctLaeAl8bwPzcI73TREZBiAijUXkVuAmN99EVd1/8MYiss6dWugFj7I/E5E7RKS/iCS6+UVEBovIm8CFbr7Jqhq/T5ONMXVm7MAOjOzp3Vp83trdvDZ/U5RrVLfiIvCp83NmPM6txq7AHBEpwGnB+QjOv8MzqvpcDYrvgnNVtxgoFpGdOC1JFwDj3DxPAPfW6iCMMaaGRIT7x/WnaaNEz/QHpi+Pq+HM4iLwAajqJmAQTpBaASQBBTi3Ns9X1WtrWPRvgeeAb4HdOHP8VQIrgX8Aw1T1Bo23ewnGmJjSuXUzbjm5p2daXlEZD0xfEeUa1R2x7+PYk52drTk5OXVdDWNMA1NeUcnYJ79i2dZ8z/T/XDmM4d3beKbVByIyX1Wzq8sXN1d8xhgT75ISE3jgZwMCDGcWH337LPAZY0wcGdS5JZcO6+KZ9sPOQqbEQd8+C3zGGBNnbju1F21TG3umPTlzNRt2FXmmNRQW+IwxJs6kNWnkt29faXklE99Z0qD79lngM8aYODRmQHuO99O3b+bKXD5a1nDn7bPAZ4wxcUhEuGdsP5KTvMPAPe8spbC0PMq1io64CHzuCC2nisgfRORtEdnijsSiInJamPaRLCK3i8giEdknInkiMsedYNZmXzfGxJys9BR+fUJ3z7Qte0t44tPVUa5RdMRF4AP6AB8AfwTGAu3DWbiIpAGzgcnAQJxxP5sCw4BngXdEJCmc+zTGmHC45vjuZLVp5pn29y9/aJCDWMdL4APIw5l54UHg3DCX/RwwFGfklrOA5jgT214OlABnAveEeZ/GGFNrTRolcu/Z/T3TyiqUSe8sbXANXeLlKuQ7oLXvsGHhuvsoIoOB892XV6jqu+7fFcCLItIS+Ctws4g8pqo7wrLjg6gqBQUF5OfnU1RUREVFw++EaoxxJCUl0aJFC1q3bk1SUuhf6yN7ZnDGgHZMX7ztkLRZ3+/k/SXbOGNAWG+U1am4CHyqWhnB4i921ytV9R2P9Ck4V3stgJ8Bz4S7AqrKjh07KCwspHXr1rRr147ExMSwBXdjTOxSVfbv38+uXbvYuHEjXbp0ISEh9Jt5fxjTl5krcikuO/RH8x/fXcYJvTJoltwwQkY83eqMlFHu+iOvRFUtBma5L0+MRAUKCgooLCykS5cutGzZkqSkJAt6xsQJEaFx48a0b9+epKQk9uzZU6NyOrRsym9GH+6ZtnVvCU82oIYuFvhqwW2t2dt9uTRA1mXu2rvHaC3l5+fTunVrEhO9pxwxxjR8IkLLli0pLCyscRn/N6Ib3dJTPNOen7WWdTtrXnYsscBXO2lA1adkS4B8VWl+b5K73R5yRCQnNzc3pEoUFRXRvHnzkLYxxjQ8zZo1o7i4uMbbJyclMGlsP8+0/RWV3PvuMs+0+sYCX+34/jQK9GmrGvjOb3RS1Smqmq2q2RkZ3qMp+FNRUWFXe8YYEhISqKysXZOGkT0zOK1fO8+0T1fs4NMV9X9El5gNfCJyt4iU13C5v67rH232TM8YE67vgTvH9KGxnxFd7p22rN5PXRSzgQ+nbom1WKLB94Z30wD5qnqHNryeoMaYBqdz62Zc62dEl3W7ivj7l2ujXKPwitnAp6qTVFVquEyIUjXzORD8OgTIV5W2NbLVMcaY8Ljm+O50bOn9e/6pT1ezI78kyjUKn5gNfPWB2yF+ufvS+4mwo6o1Z8N4MmyMafCaNErkrjP7eKYV7q9g8gcro1yj8LHAV3sz3fXJXoki0gQ4zn05Iyo1MsaYMDi1XztGHJ7umfbGgk0s2pgX5RqFhwW+2vuPu+4tImd6pF+JM2pLMfBm1GpljKlXvvnmG66//nr69etHSkoKmZmZnH/++axatarO6iQi3H1WXxITvBvNTHpnKZWV9W8cz7gJfCLSSkTSqxafpDTf90Wkkce269wpjF44OE1VFwKvui9fEJEz3G0SReQXODM2ADwaqXE6jTH13+TJk3njjTcYPXo0jz32GFdddRVffPEFQ4YMYfHixXVWr56HpXLpsC6eaYs25vHWos1RrlHtNYyB14KzEPA6e68c9HoU8FmIZV8JdMeZoeE9ESnCaVna2E1/F5gYYpnGmDhyyy23MHXqVJKTk39874ILLmDAgAE88MAD/Oc//wmwdWTddFIP3lq0mbyiskPSJn+wgtP6t6tX43jGzRVfJKlqPnAMMAH4FlCgFJgLXA2MVdWGOZWxMSYsjjnmmJ8EPYAePXrQr18/li2r23ZxLZslc+vJPT3TtueX8sznP0S5RrUTN4FPVbOC7ArxWYBtLw9Q/n5Vnayqg1S1uaq2UNXh7ogs9e8muDGmzqkq27dvJz3du4FJNF10VCa926V6pj37+Ro259V8qLRoi5vAZ0xDUVBQwDvvvMNdd93F6aefTnp6OiKCiLBixYq6rt4hNmzYwF//+lfOOussMjMzady4MampqQwcOJAJEyawdat1b/Xn3//+N5s3b+bCCy+s66qQlJjAH8Z4j7NfWl7JQx/E3mfPn/pzU9YYA8CMGTM455xz6roaQdm4cSNZWVk/mcE7LS2NwsJCvvvuO7777jumTJnCG2+8wahRowKUFH9WrFjBddddx7Bhw/jlL39Z19UBYESPdE7q05ZPlh/aTu/tRVu47JgshmS2qoOahcau+Iyph9q2bcsZZ5zBxIkTmTJlSsT2U3UlWVMVFc6YjmPGjOG1115j9+7d7N27l6KiIqZPn07Xrl3Zs2cP48aNY9u2Q2f/jlfbtm1jzJgxtGjRgjfeeCOmBqH//Rl9SPLTveGP7y6jPjzZsSs+Y+qZs846i3Hjxv34et26dXVXmWq0atWKhQsXMnDgwJ+8n5yczOmnn8706dMZPHgw+fn5PPvss0ycaI2f9+7dy+mnn05eXh6zZs2iQ4dAoyFGX7eM5vxieBb/+OrQ8ToXbsjj3e+2ctbA2KrzweyKz8S1JUuW/HhVs2HDBr/5Nm3aREpKComJiSxcuDCKNTxULP36r06LFi0OCXq+evfuzbBhwwCYP3++Z55XX30VEaFRo0asWbPGM88vfvELRISuXbuyfXvdTZtT27qWlJRw1llnsWrVKt5991369o3I3NW1duPoHrRsdkiXZ8Dp3lBSFtuzN1jgM3Gtb9++pKWlAZCTk+M334QJEygqKuLyyy9n8ODB0apeXGjTpg1w4Lbowc477zwGDhxIeXk5999/6Ixjd999Ny+//DKtW7fm/fff57DDDotofQOpTV0rKiq44IILmDNnDq+99hrDhw+PZtVD0qJZI24+ybt7w6Y9xbw4e110KxSiuAh8ItJYRE4VkT+IyNsissUdiUVF5LRalp3lU1agJTtcx2PCJyEhgSOPPBJwhozyMm/ePKZOnUrz5s257777olm9Bq+8vJyvvvoKgP79+3vmEZEf/91ffvll1q49cIvtn//8J3/84x9p3Lgxb7/9Nr179458pQOoTV1vvfVW3nnnHU4//XR2797Nv/71r58ssebiozPplpHimfbkzNXsLtwf5RoFL16e8fUBPojCfgLdYzl0yIMIy5rwXrR3GRXrHhwT1vKGDRvGjBkz/F7x3XzzzagqEyZMoH379gHLmjRpEvfcc0+N6jFx4kQmTZpUo23rq6eeeopt27aRkJDAZZdd5jffmWeeybBhw5g7dy73338/zz//PB9//DFXXXUVIsLLL7/MiBEjolhz/2pa10WLFgEwbdo0pk2bdkj6z3/+84jXPRSNEhO44/Q+XPnSof9vCkrKeeyTVdxztvePmboWL4EPIA+YD3zjLm+Eeweq2i7cZZrIq3rGlJOTg6r+pBXj1KlTmTNnDpmZmdxyyy3VltW8efMa32pr3rx5jbarr7777jvuuOMOAK6//vpqn2fdf//9jB49mpdeeomxY8dy6aWXUl5ezp///GfOO++8aFQ5aDWp62effRbdSobBSX3aMrxbG+b8sOuQtH/P28Blx2TRLSP2PtdxcasT+A5oraonqeodqvq/uq6QiR1VgS8vL4/Vq1f/+H5xcTETJjhzGj/44IM0beo9Kaev2267jW3bttVoue222yJzgNV45JFHaNeunedSxV/6I488UqN9bt26lXHjxlFcXMzQoUOZPHlytduceOKJjBo1irKyMs4++2zy8/O54YYbgvpBUuXee+8lKSmpRsudd94Z9H7CUdf6QES4c0wfvHq8lFcqk2O0U3tcBD5VrbRhw4w/6enpdO/eHfhpA5eHHnqIjRs3MmzYsJgYOSNS9u3bx/bt2z2XKv7S9+3bF/L+du/ezSmnnMLatWvp0aMH7733Hk2aNAlq2+uvv/7Hv8eOHcujjz4a0r4rKyupqKio8RKK2ta1vujfsQXnDunkmfbh0u18s253lGtUvbgIfMZUp+qqr6qBy+bNm3nooYcAePTRR2vViTvWTZo0CVX1XKr4Sw/1meTevXs59dRTWbJkCZmZmXzyySdB3xretWvXj7dGwWkFmZAQ2ldYoGOtbnnwwQeD3k846lqf3HpKT5o08j6+B6Yvj7lO7Q33TNQBEZkjIvkiUiwia0XkXyISG0/cTUAHB76q7gsXXXTRj2nBCHTbsLqlprcN64vCwkLOOOMMcnJyaNeuHZ988gmZmZlBbVtSUsLZZ5/NqlWrGDx4MAkJCbz33nvMmTMnwrUOXX2qa7i0b9GUX43o6pm2cEMe0xfH2Kg8Nf31U98XnKmDFDitluVk+ZSlOI1oig9676+AVFPOVUAOkJOZmamhWLZsWUj5zaG++eYbBTQlJUVnz56tIqJNmjTR9evXh1TOxIkT9aBzH/QyceLEGtV97dq1P5axfPnyGpXhT1W5tVVUVKSjRo1SQNu0aaNLliwJetvKykodP368AtqrVy/dtWuXXnjhhQroqFGjal23cIqFutbV90F+8X4dcu9H2uV37x6yHDf5Uy0tq4h4HYAcDeZ7O5hMDXEJY+BrBzwFHAc0d98TYAjwjs9+fh9smUOHDq32BPuywFd7+/fv16ZNmyqgXbp0UUDvvPPOuq5WUGI98JWWluppp52mgLZs2VLnz58f0vY333yzAtq2bVtds2aNqjqf+YSEBAV0xowZtapfOMVCXevy++DF2Ws9A1+X372rL3y1NuL7DzbwxeytThG5W0TKa7gcOmRChKjqNlW9TlVnqeo+9z1V1QWqOhZ4zc36exFpGa16mdA0atSIIUOGALB+/XratWv3Y4vOWLRz584flz179vz4fl5e3k/SKisr67CWzrOtiy++mA8++IDU1FTef//9H/+dg/H444/z6KOP0rRpU6ZNm0a3bt0A6NOnD+effz5ASK0tI6k+1TVSLjoqk67p3p3aH5/xPQUlUe/O7ClmAx9O3RJrscSK37nrFGB0XVbEBOb7LO/++++P6X51GRkZPy6+gWT48OE/SQs0/mg0fPXVV7zxhtNltqysjHHjxvl9xlk1gk6VN998k5tvvpmEhASmTp3KUUcd9ZP0u+66i4SEBObOncu7774btWPyUp/qGkmNEhP43Wm9PNN2Fe7nuVmHDmxdF2I28KnqJA1uxnSvJWZ+qqvqWiDXfdmtLutiAqvqpzdo0CAuv/zyuq1MA+F7xVlSUuK3W8T27dvJzc39Me/cuXO55JJLqKys5C9/+ctPZqOo0rdv3x87g991111Vjx6irj7VNRpO7deOwZneN7een/UDOwpKolyjQ0lDPgGBiEjVgZ+uqhEdzkxEdgAZwO2q+nB1+bOzszXQgMkHW758OX369KlFDU1eXh7du3dn9+7dzJw5kxNOOKGuq2RMjcTC98HXa3dz/rPerVgvOTqT+88ZEJH9ish8Va12XOSYveJrKESkK07QA4iN63xziDvvvJPdu3czfvx4C3rG1NJRXVtzUp+2nmn//WYja3cWRrlGP2WBr5ak+p7ND7jrYuDTCFfHhEhVefzxx3n66adp0aJFgx1dw5ho+91pvfGaqL2iUnnko5XRr5CPuAl8ItJKRNKrFp+kNN/3ReSQ2RVFZJ07tdALHkV/JiJ3iEh/EUl084uIDBaRN4Gqsa4mq2rsjd0Tpz788EOysrJo0aIFN954IyLC888/T6dO3kMvGWNC0+OwVM4b2tkz7b3vtrJ4094o1+iAuAl8wEKcRiZVS5VXDnr/2BDL7YJzVbcYKBaRnUAhsACoesr9BHBvjWtuwm7u3LmsX78eVWX48OFMmzaN8ePH13W1jGlQbjq5B42TvMPMQx/W3QDW8RT4IuW3wHPAt8BuIA2oBFYC/wCGqeoNGq+tiGLUxIkTUVUKCgqYPXs2Y8aEd44/Y4wzlNnlx2R5ps36fidfrd4Z3Qq54mY+PlXNisS2qvoaBzqpG2OM8XHtCd2Z+vUGCkrKD0mb/MEK3r7u2KgPAm9XfMYYYyKmZbNkrjm+u2fad5v28uHS6A9gbYHPGGNMRF1xbBYZqY090x75aBUVldF9EmSBzxhjTEQ1S07ihtE9PNNW79jHmws3R7U+FviMMcZE3AXZncls3cwz7dGPV1FaHtoM97Vhgc8YY0zEJSclcMvJPT3TNucV89+vN0atLnER+EQkQ0SuFpHXRGSNiJSISKGILBeRJ0Xk8DDsI01E7nPLLBKRXSIyQ0Ssc5gxxgBnDexAr8NSPdOe+HQ1RfsPbfkZCXER+IAtwDPAeJwZEspwunL0Bq4DFovIRTUtXEQ6AYuAO90yK3D6850IvCYiT9eq9kGwboLGmFj/HkhMEG49xfuqb+e+Ul6cvT4q9YiXwJcEfAFcBrRX1VSgGTACJ2A1AV4SkSNCLdgdq/N1oCuwDjjWLT8VuB2nM/u1InJlGI7DU2JiIhUV0bs/boyJTZWVlSQkxPbX+sl9D2NQZ+9pi575fA35UZisNrb/hcLneFU9XlVfUtVtAKpaoapfAacAO3CC4801KPts4GicAHeOqs52yy9xpyB63M13r4gk1/ZAvDRr1ox9+/ZFomhjTD1SVFT047ySsUpEuP1U78lq9xaX8fcoTFYbF4FPVb8IkJYLTHdfDq1B8Ze4609UdZFH+iOAAu1wbn2GXVpaGrt377arPmPimKqSl5dHSkpKXVelWsccns7wbm080/7+5Vr2FO6P6P7jIvAFYZe7TqzBtqPc9Ydeiaq6GVjqvoxI4EtNTSUlJYX169eTl5dHeXl5zN/rN8aEh6pSWlrK1q1bKS8vp1WrVnVdpaDcdqr3s759peU8+8UPEd133IzVWY3j3fWSUDYSkbZA1c+WpQGyLgP6A31Dr1pQ9aBt27YUFBSQn5/Pjh077OrPmDiSlJREixYtaNu2bcw/46sytEtrRvXKYObK3EPSXpi9ll+OyKJtapOI7DvuA5+InA1UTVX/zxA3b+/z95YA+arS2gfIUysiQlpaGmlpaZHahTHGhNWtp/TyDHwlZZU8PXMNk8b2i8h+68dPgwgRkY7AFPflO6r6QYhF+N5MLw6Qr8hdNw9Ql6tEJEdEcnJzD/0gGGNMQ9O/YwtO79/OM23qvA1szgv0tVpzMRv4RORuESmv4XJ/EOU3B94C2gLrgV9F+pgCUdUpqpqtqtkZGRl1WRVjjImam0/uidesRCN7ZlBREZm2CrF8qzOBmjU2obrtRKQJ8DbOLc5c4FRVrcmMiIU+fwdqQ1w1QJ31OTDGGB89D0tl3KCOPw5UfVyPdG45uSeDMyPXSCdmA5+qTgImhbtcty/d6zgtLPOAU1R1ZQ2L832u1wFY7CdfB3e9tYb7McaYBuvG0T3YureYm07qyTA/3RzCKWYDXySISBLwH2AMztXXGX763gVFVXNFZCeQDvTDT5cGDrTmXFbTfRljTEOVlZ7Cf68aHrX9xewzvnATkQTgReBnOA1RxqrqnDAUPdNdn+xnvx1xgiLAjDDszxhjTC3EReBzx9OcAlwM7Ad+pqozA28VtKnu+hQRGeiRfgsgOLc5w7VPY4wxNRQXgQ94FKfVZjlwfqjdFkRE3WWSR/LbwDycf8s3RWSYu01jEbkVuMnNN1FVIzsOjzHGmGo1+Gd8IpIJ3Oi+VOBZEXnWX35V9e5U4j+/unPufYEzQ8McEdmHM+ND1b/vM6r6XMiVN8YYE3YNPvDx06vaRsBh4d6Bqm4SkUHA73CeIWYBBThTHv1NVV8L9z6NMcbUjNhgxrFHRHJxOtXXRDpQkz6JsaihHEtDOQ6wY4lVDeVYanscXVS12hFALPA1MCKSo6rZ1eeMfQ3lWBrKcYAdS6xqKMcSreOIl8YtxhhjDGCBzxhjTJyxwNfwTKk+S73RUI6loRwH2LHEqoZyLFE5DnvGZ4wxJq7YFZ8xxpi4YoHPGGNMXLHAZ4wxJq5Y4ItBItJORB4TkTUiUiIi20VkmoiMrmW5aSJyn4gsF5EiEdklIjPcIdciItzHIiIn+IydGmhJD+MxpIrIWBH5o4i8LyI7ffbTOwzlR+28ROpY6ui8ZIrITe7naYOIlIpIgYh8KyIPikj7WpYfkf+HHvuJyHGIyOVBnI+wTo4tItnuZ+sDEVktInvd49ksIm+LyLhalh+ec6KqtsTQAhyBM3KBusteoML9uxKYUMNyOwE/+JRbAJT5vH66PhwLcIK7fQWwLcDSOozHMc7nGA5eetey7Kiel0gdS7TPC9DZ/Qz51n8vzkD0Va93A6Ni5bMb7eMALne33x/gfKwJ8+frmYOOpQBnGjjf914HGtXlOQnbAdsSlg9NU2CdeyIXAP3c99OAR3xO8CkhlivAXHf7tcAx7vtNgN/6fHiurAfHUvUFuqvkZwAACOVJREFUuy6K52UcsB14D5gEXOnzn682waIuzkukjiWq5wVnPNxK4F1gPNDKfT8ZOJ0DPyb2Au1i4bNbB8dxubvtZ9E4J+4+L8OZkWYI0Nzn/c7AQz6ftbvq8pxE5R/DlqBP7k0c+JXU0SP9TTd9fojlVv3KrwAGeaQ/6qZvBZJj/Fii+gXr7jPxoNdZYQoWdXFeInUsUT0vQAtgYID03hy40pgYYtkR+ezWwXFEPfAFUaeX3TqFdKUZ7nNiz/hiyyXueqqqbvZIf9hdDxGRXjUo9xNVXeSRXvWLqR1wYgjlBrPPcB9L1KlqRYSKjvp5ieCxRJWq7lXVbwOkr8C5mgYYGmLxUfvsRvg4YtE37rpDiNuF9ZxY4IsRIpLKgQ/2h36yzcW55QEQysPcUYHKdT9IS92Xtf6CjfCxNCRRPS9xaJe7Tgx2gxj97IZ8HDHsGHe9NtgNInFOLPDFjj44z3zgwJfdT6hqJbDSfdk3mEJFpC3QJlC5rmWhlFuNiBzLQTJEZIGIFLrLKhGZIiIDalBW1NXReYmGmDgvIpIEHOu+XBLCptH47AatFsfhq5+ILBWRYre16BIReVREuoapmgGJSHMROUJEngIucN9+MoQiwn5OLPDFDt8my1sC5KtKC7aJc6TKret9NgMGA6U4Eyr3wGmosVBEbqtBedFWF+clGmLlvFyHc4u4EngxhO1i7bzU9Dh8peMEjyKchlP9cJ6ZLRWRi8NRyYOJSKeqLhM4z+W+BX4NlOA0bHk6hOLCfk4s8MWOFJ+/iwPkK3LXzeu43LraZx7O/fxsoKmqtsb5sj0emI1zO+jhSP2HDqO6OC+RFDPnRUSOAP7kvnxSVZcFyn+QmDkvtTwOcALBRKA/0ERV2+DUdwzOnYSmwIsiMjJMVfZVgdN6eDtOdwpwumj8CXgqxLLCfk4s8Jl6RVUXqertqjpfVUvc9ypU9QucZ2ZfuVkni4h9vqMkVs6L29n7LZwv9fnA7yK1r0gKx3Go6keqeq+qLlXV/e57pao6HedZ22qcq/IHw1fzH/e9VVXbqWo7nGPoBbwE3AMsEpF+4d5nKOyLIXYU+vzdNEC+Zu462BEXIlVurO0T9z/3Xe7LTji33GJVnfwb1YVonRcRaQ18BHQFvgfGVAXhENT5eQnTcQSkqnuBB9yXw8I5oo7HvipVdZWq/gr4C5AJvBzCD6CwnxMLfLHD9951oKa+VWlb67jcWNtnlXk+f3cLY7nhVpf/RnUhoudFRFrgtPjrD2wATlLV7TUoqk7PSxiPIxhV50Rwgmw0POGuBxP8D6CwnxMLfLFjBU6fLXAePh/C/YVU1UclqPv9qpqLM8yP33JdVS2hQn2O4CUix9KQ1NF5aZBEJAWYjvN8cRtOsNhQw+Lq7LMb5uOIVb598LoHuU3Yz4kFvhihqgVAjvvyZD/ZjsYZ6QFgRgjFzwxUroh05MAHKpRyPUX4WKpztM/fQfcVqiNRPS91LCLnRUSaAtNwnlntwgkW39e0vLr67Ib7OILke07WRXhfVXyvLIO6TRyRc1LXQ9jY4jksTz7Q3iP9DTc9J8RyfYfGOmR4JODPbvoWwj9kWbiPRQKkNQJm+RxLQoTOUxbhH7IsKuclgscS9fOCM57l+265e4AhYSo3Ip/daB5HoPPhpqfh9HtTYF6YjiMxiP0+5+6zDHdc0ro4J7U+WFvCt/DTgVjnA33d91P56QCvhwzE6pM2ySPNdzDkH4Bh7vuNgVuJ/CDV4TyWpcBvcPqHifteIjAC+MJn28vCfG7SfZbBPvsZdlBaQgjHEvXzEsFjiep5cct+3efLcFgI22b51OfycH52Y+U43LS5wK+ATJ/3k4HTgMUc+NF1YpjOSZb77/VLoJPP+wnAIODfPvX9S12ek7D9Z7IlPAswkBpMvRHoS8lNr276m7/Vh2PxSVOczrC5OJ2lq94rA34XgWPRIJesenBewn4s0T4vwEifsosJPBXSNwdtm+Wz7eXh/OzGynEclFZV9k6cPnVV7xUCl4bxnPx/e3cTakUZx3H8+yOFopWUYbRxIxEUiUWLtkIv6MpNi6SijWDl1nTTLigI0YXQou5CV9laRZFwIYUYGRHhokUkSiRCREG4+Ld4zqWDnXuu5+i94+X5fuAs5oUz/+GZmR8zz7xMWubvo+1hfPwCsG7INlmH7itV9X2Sp4EDwE7gCdo1/4vAoaqaq0+hqq4m2Up7HmgXbUP7E7hMO7ieuAfl377MlViXPbRXOD0HPAZsoO1gV4DztHVZMzeCDNEuK2S122X8/oQHR7+lzPwowErthxOs1Hr8BuyjnXE/C2yk9YH9RXtE4hytTX6ZqdrprtFeSbYdeIH2BpVHaHX/DHwNLFTVhSX/YYp72SaLlyQkSeqCd3VKkrpi8EmSumLwSZK6YvBJkrpi8EmSumLwSZK6YvBJkrpi8EmSumLwSZK6YvBJkrpi8EmSumLwSZK6YvBJkrpi8EmaKMmHSSrJjSSbJkxPktOjeb5Nsn6IOqVZGXySlvIB8B3tm2qfT5j+DvAy7bt7u6vq1irWJs3N4JM00SjIXqcF26tJ9i5OS/Ik8PFocH9V/TRAidJc/BCtpKmSvAccAf4GtvHf17SfB84Ar5QHEq0hBp+kqZIEOEW7rHkJOAscAG4Cz1TVtQHLk2Zm8ElaVpLHgR9o/X2LXquqLwYqSZqbfXySllVV14GDY6NOGHpaqww+SctK8gDw5tiorUkeHqoe6W4YfJLuxPvAi8AfwK/AFuCTQSuS5mQfn6SpkmwDvgHWA28AV4FzQIAdVXVywPKkmXnGJ2lJSR4CjtNC78uqOlZVXwGHRrN8luTRwQqU5mDwSZrmI+Ap4DqwZ2z8QeBHYBPw6QB1SXMz+CRNlOQl4N3R4NtVdXNxWlX9A+wGbgG7kry1+hVK8zH4JP1Pkg3AAq0f72hVnb59nqq6THufJ8DhJJtXrUDpLnhziySpK57xSZK6YvBJkrpi8EmSumLwSZK6YvBJkrpi8EmSumLwSZK6YvBJkrpi8EmSumLwSZK68i8I5S4/SFgCIgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true,
        "trusted": false,
        "id": "ZGCeNUS4jKX5",
        "colab_type": "text"
      },
      "source": [
        "### Writing to a file\n",
        "\n",
        "Often we need to save data to a file, for later/further analysis, and with functions typically this is done as one column each per function (array), as that is the format most plotting/analysis software expects data to be in. Just below is the Python for saving $t$ and $x$ arrays as two columns, which each number written out to 4 decimal places."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "trusted": true,
        "id": "-BocF1LijKX6",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# we need to specify the filename to be written to, this writes the data to this filename in tbe current directory\n",
        "filename='xy.txt'\n",
        "# the numpy savetxt commands saves data in arrays\n",
        "#note we need to use numpy's column_stack to put the arrays into columns\n",
        "# also the fmt argument allows us to specify the format of the numbers\n",
        "np.savetxt(filename,np.column_stack([x,y]),fmt=\"%9.4f\")"
      ],
      "execution_count": 3,
      "outputs": []
    }
  ]
}