Section_09_Plotting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Plotting\n",
    "\n",
    "- [Overview](#Overview)\n",
    "- [Example - Cannonball Trajectory](#Example---Cannonball-Trajectory)\n",
    "- [Example - Numerical Integration](#Example---Numerical-Integration)\n",
    "- [Example - Approximation pi](#Example---Approximation-pi)\n",
    "- [Introduction to Matrices](#Introduction-to-Matrices)\n",
    "- [Example - Potential Flow Around a Cylinder](#Example---Potential-Flow-Around-a-Cylinder)\n",
    "- [Meshgrid Explanation](#Meshgrid-Explanation)\n",
    "- [Example - 3D Visualisation of a function](#Example---3D-Visualisation-of-a-function)\n",
    "- [Example - 3D Representation of a geometry object](#Example---3D-Representation-of-a-geometry-object)\n",
    "- [Publication quality plots](#Publication-quality-plots)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overview\n",
    "\n",
    "- Done by using the `matplotlib` module\n",
    "    - `from matplotlib import pyplot`\n",
    "\n",
    "\n",
    "- Types of 2D plots:\n",
    "    - `pyplot.plot` $\\to$ line or scatter plots\n",
    "    - `pyplot.semilogx` $\\to$ line or scatter plots $\\to$ log scale on the x axis\n",
    "    - `pyplot.semilogy` $\\to$ line or scatter plots $\\to$ log scale on the y axis\n",
    "    - `pyplot.loglog` $\\to$ line or scatter plots $\\to$ log scale on both the x and y axis\n",
    "    - `pyplot.bar` $\\to$ vertical bar chart\n",
    "    - `pyplot.barh` $\\to$ horizontal bar chart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overview\n",
    "\n",
    "- See [matplotlib](https://matplotlib.org/gallery.html) for many example plots\n",
    "\n",
    "\n",
    "- General usage of the `matplotlib` module:\n",
    "    1. `from matplotlib import pyplot`\n",
    "    2. `pyplot.figure(num)` $\\to$ create a blank figure canvas, given a figure number\n",
    "    4. Collect / generate data\n",
    "    3. `pyplot.plot(xdata, ydata, plotstyle)` $\\to$ plot data to the figure canvas\n",
    "    4. `pyplot.title(Title)` $\\to$ create a figure title\n",
    "    5. `pyplot.xlabel(X Label)` $\\to$ create a x axis label\n",
    "    6. `pyplot.ylabel(Y Label)` $\\to$ create a y axis label\n",
    "    7. `pyplot.legend(loc=num)` $\\to$ create a legend at a given location on the figure\n",
    "    8. `pyplot.show()` $\\to$ show the figure canvas with plotted data\n",
    "\n",
    "\n",
    "- See the `help()` information for more on each of the above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "xdata = [0, 1, 3]\n",
    "ydata = [0, 1, -2]\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(xdata, ydata, '-b')\n",
    "\n",
    "plt.figure(2)\n",
    "plt.plot(xdata, ydata, 'or')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "help(plt.plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overview\n",
    "\n",
    "- `pyplot.plot(xdata, ydata, plotstyle)`:\n",
    "    - `xdata` is a `list` or `array` of x axis data\n",
    "    - `ydata` is a `list` or `array` of corresponding y axis data\n",
    "    - **both lists or arrays need to be the same length**\n",
    "    - `plotstyle` $\\to$ string that defines the style of the line and the markers $\\to$ see `help(pyplot.plot)`\n",
    "\n",
    "\n",
    "- For example\n",
    "    - `pyplot.plot(xdata, ydata, b-)` $\\to$ blue (`b`), solid (`-`) line\n",
    "    - `pyplot.plot(xdata, ydata, g--)` $\\to$ green (`g`), dashed (`--`) line\n",
    "    - `pyplot.plot(xdata, ydata, ro)` $\\to$ red (`r`), circle (`o`) markers\n",
    "    \n",
    "Note the effect of the different `%matplotlib` magics in the next cell.  These magics do not execute any plotting, but\n",
    "they determine how Jupyter Notebook handles figures.\n",
    "- The `%matplotlib inline` magic causes the figure to appear seamlessly inside your Jupyter Notebook and is ideal for report writing in Jupyter Notebook.\n",
    "- The `%matplotlib notebook` magic puts the figure inside a frame inside the notebook, with a upper bar that indicates the figure number and which also includes a switch on the far right of the bar, with which the specific figure can be switched to the `inline` mode.  The `notebook` mode also includes a bar with tool buttons in the figure frame, which includes a zoom and a pan tool.  It also shows the horixontal and vertical coordinates of the mouse pointer position on the graph, if the mouse is hovered over the graph.\n",
    "- The `%matplotlib qt5` creates a mode that is essentially the same as the `notebook` mode, except that each figure is in its own window that is opened by JUpyter Notebook, outside of the notebook and browser.  The `qt5` magic does sometimes not run properly on some machines.\n",
    "\n",
    "Only one of these `%matplotlib` magics can be activated at a time.  Also, you need to restart the kernel if you activate a different magic and expect Jupyter Notebook to behave properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%matplotlib notebook\n",
    "#%matplotlib qt5\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "xdata = np.linspace(-10, 10, 100)\n",
    "ydata = xdata**2\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(xdata, ydata)\n",
    "#plt.plot(ydata)\n",
    "#plt.plot(ydata,'.')\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: with the %matplotlib inline magic, a plt.plot command in a new cell creates a new figure.\n",
    "#       On the other hand, with the %matplotlib notebook magic, if the figure number is not\n",
    "#       incremented, the new graph is plotted on the last figure that was plotted in a previous cell.\n",
    "xdata1 = np.linspace(0, 10, 100)\n",
    "ydata1 = xdata1**0.5\n",
    "#plt.figure()\n",
    "plt.plot(xdata1, ydata1,'g')\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example - Cannonball Trajectory\n",
    "\n",
    "- A cannonball is fired from a cannon with an initial $x$ and $y$ velocity ($V_x^0$, $V_y^0$).\n",
    "\n",
    "\n",
    "- The cannonball exits the cannon at an initial $x$ and $y$ coordinate of $x^0 = 0$ and $y^0 = 0$.\n",
    "\n",
    "\n",
    "- The $x$ and $y$ coordinates of a cannonball trajectory are updated by:\n",
    "$$\n",
    "\\begin{align}\n",
    "    x^{n+1} &= x^n + \\Delta t \\left( V_x^n \\right) \\\\\n",
    "    y^{n+1} &= y^n + \\Delta t \\left( V_y^n \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "- The velocities of the cannonball are updated by:\n",
    "$$\n",
    "\\begin{align}\n",
    "    V_x^{n+1} &= V_x^n - \\Delta t \\left( \\frac{2 V_x^n}{m} \\right) \\\\\n",
    "    V_y^{n+1} &= V_y^n - \\Delta t \\left( \\frac{2 V_y^n}{m} + 9.81 \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "- where $\\Delta t$ is a constant time step and $m$ is the mass of the cannon ball.\n",
    "\n",
    "\n",
    "- Write a python program that visualises the trajectory of the cannonball for various different masses ($m$) and initial conditions ($x^0$, $y^0$, $V^0_x$, $V^0_y$).\n",
    "\n",
    "\n",
    "- Take $\\Delta t$ as 0.01 s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Outcomes:\n",
    "\n",
    "-   Basic 2D line plotting\n",
    "\n",
    "\n",
    "-   Figure annotations and ``fontsize``\n",
    "\n",
    "\n",
    "-   Multiple figures\n",
    "\n",
    "\n",
    "-   Sub-plots\n",
    "\n",
    "\n",
    "-   Multiple plots on one figure\n",
    "\n",
    "\n",
    "-   Legend and legend location\n",
    "\n",
    "\n",
    "-   Pyplot `axis('equal')` command to ensure a circle is plotted as a circle and not an ellipse (i.e., forcing the horizontal and vertical axes to be on the same scale).\n",
    "\n",
    "\n",
    "-   Saving figures in different formats in different ways, including the use of the Pyplot `savefig` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "def next_coordinates(x, y, Vx, Vy, dt):\n",
    "    x = x + dt*Vx\n",
    "    y = y + dt*Vy\n",
    "    return x, y\n",
    "\n",
    "\n",
    "def next_velocities(Vx, Vy, dt, mass, g):\n",
    "    Vx = Vx - dt * (2*Vx / mass)\n",
    "    Vy = Vy - dt * (2*Vy / mass + g)\n",
    "    return Vx, Vy\n",
    "\n",
    "\n",
    "def trajectory(x0, y0, Vx0, Vy0, dt, mass, g):\n",
    "    x = x0\n",
    "    y = y0\n",
    "    Vx = Vx0\n",
    "    Vy = Vy0\n",
    "\n",
    "    xcoords = []\n",
    "    ycoords = []\n",
    "    while y >= 0:\n",
    "        x, y = next_coordinates(x, y, Vx, Vy, dt)\n",
    "        Vx, Vy = next_velocities(Vx, Vy, dt, mass, g)\n",
    "        xcoords.append(x)\n",
    "        ycoords.append(y)\n",
    "    return xcoords, ycoords\n",
    "\n",
    "\n",
    "def plot_trajectory(x0, y0, vel, theta):\n",
    "    Vx0 = vel * np.cos(numpy.radians(theta))\n",
    "    Vy0 = vel * np.sin(numpy.radians(theta))\n",
    "    xcoords, ycoords = trajectory(x0, y0, Vx0, Vy0, dt=0.01, mass=10, g=9.81)\n",
    "    \n",
    "    plt.plot(xcoords, ycoords)\n",
    "    plt.title(\"Cannonball Trajectory\", fontsize=16)\n",
    "    plt.xlabel(\"Distance [m]\", fontsize=12)\n",
    "    plt.ylabel(\"Height [m]\", fontsize=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%matplotlib qt5\n",
    "#%matplotlib notebook\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Figure 1\n",
    "plt.figure(1)\n",
    "plot_trajectory(x0=0, y0=0, vel=40, theta=25)\n",
    "#plt.axis('equal')\n",
    "#plt.savefig('figure_Python-generated.jpg')\n",
    "#plt.savefig('figure_Python-generated.pdf')\n",
    "\n",
    "# Figure 2\n",
    "plt.figure(2)\n",
    "plot_trajectory(x0=0, y0=0, vel=40, theta=45)\n",
    "#plt.axis('equal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Figure 1 - Subplot 1\n",
    "plt.subplot(1, 2, 1)\n",
    "#plt.subplot(2, 1, 1)\n",
    "plot_trajectory(x0=0, y0=0, vel=40, theta=25)\n",
    "\n",
    "# Figure 1 - Subplot 2\n",
    "plt.subplot(1, 2, 2)\n",
    "#plt.subplot(2, 1, 2)\n",
    "plot_trajectory(x0=0, y0=0, vel=40, theta=45)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Figure 1 - Plot 1\n",
    "plot_trajectory(x0=0, y0=0, vel=40, theta=25)\n",
    "\n",
    "# Figure 1 - Plot 2\n",
    "plot_trajectory(x0=0, y0=0, vel=40, theta=45)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "x0 = 0\n",
    "y0 = 0\n",
    "vel = 40\n",
    "for theta in [25, 35, 45, 55, 65]:\n",
    "    Vx0 = vel * np.cos(numpy.radians(theta))\n",
    "    Vy0 = vel * np.sin(numpy.radians(theta))\n",
    "    xcoords, ycoords = trajectory(x0, y0, Vx0, Vy0, dt=0.1, mass=10, g=9.81)\n",
    "    \n",
    "    label = \"Theta: {}\".format(theta)\n",
    "    plt.plot(xcoords, ycoords, label=label)\n",
    "    \n",
    "plt.title(\"Cannonball Trajectory\", fontsize=16)\n",
    "plt.xlabel(\"Distance [m]\", fontsize=12)\n",
    "plt.ylabel(\"Height [m]\", fontsize=12)\n",
    "plt.legend(loc=\"upper right\")\n",
    "plt.grid(True)\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example - Numerical Integration\n",
    "\n",
    "- In numerical analysis, the trapezoidal rule is a technique for approximating the definite integral:\n",
    "\n",
    "$$ \\int_{a}^{b} f(x) \\;dx \\approx \\sum_{n=0}^{N-1} \\frac{\\Delta x}{2} \\left[ f(x_{n+1}) + f(x_n) \\right]$$\n",
    "\n",
    "where \n",
    "$$ \\Delta x = \\frac{b - a}{N}$$\n",
    "\n",
    "<img src=\"./figures/trapezoid_rule.svg\" alt=\"Numerical Integration\" style=\"height: 250px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Write a python program that calculates the numerical integral of:\n",
    "    $$ f(x) = x^2$$\n",
    "    between\n",
    "    - $a = x_l = -5$ and\n",
    "    - $b = x_u = 5$, \n",
    "    - using the trapezoidal rule.\n",
    "\n",
    "\n",
    "- The program must also create a similar graph to one shown in the previous slide, to visualise the numerical integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Outcomes:\n",
    "\n",
    "-   Using a function as an input to another function\n",
    "\n",
    "\n",
    "-   Multiple plots on one figure\n",
    "\n",
    "\n",
    "-   Using a function to plot data / lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "def my_function(xval):\n",
    "    fx = xval**2\n",
    "    return fx\n",
    "\n",
    "\n",
    "def integrate(func, lower, upper, num):\n",
    "    integral = 0\n",
    "    dx = (upper - lower) / num\n",
    "    xvals = np.linspace(lower, upper, num+1)\n",
    "    yvals = func(xvals)\n",
    "    for i in np.arange(0, num, 1):\n",
    "        area = 0.5 * dx * (yvals[i] + yvals[i+1])\n",
    "        integral = integral + area\n",
    "        plot_square(xvals[i], xvals[i+1], yvals[i], yvals[i+1])\n",
    "    return integral\n",
    "\n",
    "\n",
    "def plot_function(func, lower, upper, num):\n",
    "    xvals = np.linspace(lower, upper, num)\n",
    "    yvals = func(xvals)\n",
    "    plt.plot(xvals, yvals,'b-')\n",
    "\n",
    "\n",
    "def plot_square(x0, x1, y0, y1):\n",
    "    plt.plot(\n",
    "        [x0, x1, x1, x0, x0], \n",
    "        [ 0,  0, y1, y0,  0],\n",
    "        'g-'\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "num = 4\n",
    "#num = 40\n",
    "low = -5.0\n",
    "upp =  5.0\n",
    "# low = -4.0\n",
    "# upp =  8.0\n",
    "\n",
    "plot_function(my_function, low, upp, 100)\n",
    "integral = integrate(my_function, low, upp, num)\n",
    "\n",
    "plt.title(\"Trapezoidal Integration\", fontsize=16)\n",
    "plt.xlabel(\"x [radians]\", fontsize=12)\n",
    "plt.ylabel(\"f(x)\", fontsize=12)\n",
    "#plt.show()\n",
    "\n",
    "print(\"Integral:\", integral)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check result with a simple analytical calculation:\n",
    "\n",
    "$$ \\int_{-5}^{5} x^2 \\; dx = \\left. \\frac{x^3}{3} \\right|_{-5}^{5} = \\frac{5^3}{3} - \\frac{(-5)^3}{3} $$\n",
    "\n",
    "and\n",
    "\n",
    "$$ \\int_{-4}^{8} x^2 \\; dx = \\left. \\frac{x^3}{3} \\right|_{-4}^{8} = \\frac{8^3}{3} - \\frac{(-4)^3}{3} \\;.$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(5**3/3 - ((-5)**3/3))\n",
    "print(8**3/3 - ((-4)**3/3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example - Approximation pi\n",
    "\n",
    "<img src=\"./figures/darts_pi_approx.svg\" alt=\"Pi Approximation Darts\" style=\"height: 300px;\"/>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Areas: \n",
    "\n",
    "    $$ \\frac{Ac}{As} = \\frac{\\frac{\\pi R^2}{4}}{\\frac{4 R^2}{4}} = \\frac{\\pi}{4} $$\n",
    "\n",
    "\n",
    "- Approximation (Throwing darts):\n",
    "\n",
    "    $$\\frac{\\text{darts in circle}}{\\text{darts in square}} \\approx \\frac{\\pi}{4}$$\n",
    "    $$\\therefore \\; \\pi \\approx 4 \\left( \\frac{\\text{darts in circle}}{\\text{darts in square}} \\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Outcomes:\n",
    "\n",
    "- `numpy` module and `linspace` function\n",
    "\n",
    "\n",
    "- Idea of `arrays` and working with arrays\n",
    "\n",
    "\n",
    "- Line plot with a Scatter plot\n",
    "\n",
    "\n",
    "- Histogram plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "def random_dart_position(radius):\n",
    "    x = np.random.uniform(-radius, radius)\n",
    "    y = np.random.uniform(-radius, radius)\n",
    "    return x, y\n",
    "\n",
    "\n",
    "def in_the_circle(x, y, radius):\n",
    "    dist = (x**2 + y**2) ** 0.5\n",
    "    return (dist < radius)\n",
    "\n",
    "\n",
    "def pi_approx(radius, tol):\n",
    "    error = 1\n",
    "    cnt_circle = 0\n",
    "    cnt_square = 0\n",
    "    xcoords = []\n",
    "    ycoords = []\n",
    "    while (error >= tol):\n",
    "        x, y = random_dart_position(radius)\n",
    "        xcoords.append(x)\n",
    "        ycoords.append(y)\n",
    "        cnt_square = cnt_square + 1\n",
    "        if in_the_circle(x, y, radius):\n",
    "            cnt_circle = cnt_circle + 1\n",
    "        approx = 4.0 * cnt_circle / cnt_square\n",
    "        error = abs(np.pi - approx)\n",
    "    return approx, xcoords, ycoords\n",
    "\n",
    "\n",
    "def plot_square(radius):\n",
    "    plt.plot(\n",
    "        [-radius,  radius, radius, -radius, -radius],\n",
    "        [-radius, -radius, radius,  radius, -radius],\n",
    "        'b-'\n",
    "    )\n",
    "\n",
    "\n",
    "def plot_circle(radius):\n",
    "    xvals = np.linspace(-radius, radius, 100)\n",
    "    yvals = (radius**2 - xvals**2) ** 0.5\n",
    "    plt.plot(xvals,  yvals, 'r-')\n",
    "    plt.plot(xvals, -yvals, 'r-')\n",
    "\n",
    "\n",
    "def plot_darts(radius):\n",
    "    distances = []\n",
    "    approx, xcoords, ycoords = pi_approx(radius, tol=1e-3)\n",
    "    for i in range(0, len(xcoords), 1):\n",
    "        dist = np.hypot(xcoords[i], ycoords[i])\n",
    "        distances.append(dist)\n",
    "        if dist <= radius:\n",
    "            plt.plot(xcoords[i], ycoords[i], 'r*')\n",
    "        else:\n",
    "            plt.plot(xcoords[i], ycoords[i], 'b*')\n",
    "\n",
    "    offset = radius + 0.1\n",
    "    plt.title(\"Pi Approximation\", fontsize=16)\n",
    "    plt.xlabel(\"Width [m]\", fontsize=12)\n",
    "    plt.ylabel(\"Height [m]\", fontsize=12)\n",
    "    plt.axis([-offset, offset, -offset, offset])\n",
    "    plt.axis('equal')\n",
    "    plt.grid()\n",
    "    return distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "y = \\pm \\sqrt{R^2 - x^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%matplotlib qt5\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "radius = 1\n",
    "plot_square(radius)\n",
    "plot_circle(radius)\n",
    "distances = plot_darts(radius)\n",
    "\n",
    "# create histogram\n",
    "plt.figure(2)\n",
    "plt.hist(distances, bins=20)\n",
    "\n",
    "plt.title(\"Frequency of Point From the Origin\", fontsize=16)\n",
    "plt.xlabel(\"Distance from origin [m]\", fontsize=12)\n",
    "plt.ylabel(\"Frequency\", fontsize=12)\n",
    "plt.grid(True)\n",
    "\n",
    "binedges = np.linspace(min(distances), max(distances), 21)\n",
    "print(binedges)\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(plt.hist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to Matrices\n",
    "\n",
    "- 2D data structure $\\to$ 2D indexing\n",
    "$$\n",
    "\\text{col} \\\\\n",
    "\\text{row}\n",
    "\\begin{array}{l|cccc}\n",
    "      & 0 & 1 & 2 & 3 \\\\\n",
    "    \\hline\n",
    "    0 & 10 & 12 & 14 & 16 \\\\\n",
    "    1 & 20 & 22 & 24 & 26 \\\\\n",
    "    2 & 30 & 32 & 34 & 36 \\\\\n",
    "    3 & 50 & 52 & 54 & 56 \\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "\n",
    "- Indexing $\\to$ `matrix[row, col]` $\\to$ `row` first, then `col`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "mat = np.array([\n",
    "    [10, 12, 14, 16],  # 0\n",
    "    [20, 22, 24, 26],  # 1\n",
    "    [30, 32, 34, 36],  # 2\n",
    "    [50, 52, 54, 56],  # 3\n",
    "])\n",
    "#     0   1   2   3\n",
    "\n",
    "print('Matrix mat:')\n",
    "print(mat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Numpy array **slicing**; first 1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aaa = np.linspace(0,10,11)\n",
    "print(aaa)\n",
    "print(aaa[2:4])\n",
    "print(aaa[-1])\n",
    "print(aaa[-2])\n",
    "print(aaa[5:-1])\n",
    "print(aaa[5:-2])\n",
    "print(aaa[5:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Numpy array **slicing**; 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Matrix mat:')\n",
    "print(mat)\n",
    "print('\\nElement 2nd row, 3rd column = ',mat[1,2])\n",
    "print('\\nColumn 3:')\n",
    "print(mat[:, 2])\n",
    "print('\\nRow 2:')\n",
    "print(mat[1, :])\n",
    "print('\\nsub-matrix, rows 2 & 3, columns 2 to the last:')\n",
    "print(mat[1:3,1:4])\n",
    "print('\\nsub-matrix, rows 2 & 3, columns 2 to the last:')\n",
    "print(mat[1:3,1:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `array` operators:\n",
    "    - `+ − ∗ / ∗∗` etc $\\to$ same as normal number mathematics → done on an element-by-element bases\n",
    "        - Can do operations with `arrays` and numbers\n",
    "            - E.g. `arr3 = arr1 + 12.5`\n",
    "        - Can do operations with `arrays` and `arrays`\n",
    "            - E.g. `arr3 = arr1 + arr2` $\\to$ **arrays must be the same shape !!**\n",
    "    - Operator priority $\\to$ same as mathematical priority"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.info(np.zeros)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.zeros( (2, 5) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "mat = np.array([\n",
    "    [10, 12, 14, 16],\n",
    "    [20, 22, 24, 26],\n",
    "    [30, 32, 34, 36],\n",
    "    [50, 52, 54, 56],\n",
    "])\n",
    "\n",
    "print(mat * 2, \"\\n\")\n",
    "print(mat * mat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional Types of Plots\n",
    "\n",
    "- Contour plots\n",
    "    - 3D data represented on a 2D plane\n",
    "    - 3 matrices required: `X`, `Y`, `Z`\n",
    "\n",
    "\n",
    "- Quiver plot\n",
    "    - Vector fields\n",
    "    - 3D data represented on a 2D plane\n",
    "    - 4 matrices required: `X`, `Y`, `Zx`, `Zy`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "x = np.linspace(-2, 2, 5)\n",
    "y = np.linspace(-2, 2, 5)\n",
    "\n",
    "Xmat, Ymat = np.meshgrid(x, y)\n",
    "print(Xmat, '\\n')\n",
    "print(Ymat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Meshgrid Explanation\n",
    "\n",
    "- Physical \"mesh\" points are structured (we are trying to create a structured \"grid\" of points):\n",
    "\n",
    "$$\n",
    "\\newcommand{\\point}[2]{\\begin{array}{cc} \\bullet \\\\[-0.5em] \\left( x_#1, y_#2 \\right)\\end{array}}\n",
    "\\begin{align}\n",
    "    \\point{0}{2}\n",
    "    \\point{1}{2}\n",
    "    \\point{2}{2}\n",
    "    \\point{3}{2}\\\\\n",
    "    \\point{0}{1}\n",
    "    \\point{1}{1}\n",
    "    \\point{2}{1}\n",
    "    \\point{3}{1}\\\\\n",
    "    \\point{0}{0}\n",
    "    \\point{1}{0}\n",
    "    \\point{2}{0}\n",
    "    \\point{3}{0}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Thus, along the $x$ axis $\\to$ `x = np.linspace(`$x_0$, $x_3$`, 4)`\n",
    "\n",
    "\n",
    "- and, along the $y$ axis $\\to$ `y = np.linspace(`$y_0$, $y_2$`, 3)`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The \"mesh\" points can also be represented by 2 matrices\n",
    "    - One matrix of $x$ coordinates, and\n",
    "    - Another matrix of $y$ coordinates:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\point{0}{2}\n",
    "    \\point{1}{2}\n",
    "    \\point{2}{2}\n",
    "    \\point{3}{2}\\\\\n",
    "    \\point{0}{1}\n",
    "    \\point{1}{1}\n",
    "    \\point{2}{1}\n",
    "    \\point{3}{1}\\\\\n",
    "    \\point{0}{0}\n",
    "    \\point{1}{0}\n",
    "    \\point{2}{0}\n",
    "    \\point{3}{0}\n",
    "\\end{align}\n",
    "\\Rightarrow\n",
    "\\begin{bmatrix}\n",
    "    x_0 & x_1 & x_2 & x_3 \\\\\n",
    "    x_0 & x_1 & x_2 & x_3 \\\\\n",
    "    x_0 & x_1 & x_2 & x_3 \\\\\n",
    "\\end{bmatrix}\n",
    "\\; \\text{and} \\;\n",
    "\\begin{bmatrix}\n",
    "    y_0 & y_0 & y_0 & y_0 \\\\\n",
    "    y_1 & y_1 & y_1 & y_1 \\\\\n",
    "    y_2 & y_2 & y_2 & y_2 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "\n",
    "- **Note:** The first row in a matrix is at the top (`mat[0, c]`), thus the $y_0$ coordinates appear in the first matrix row\n",
    "\n",
    "\n",
    "- The `meshgrid` function thus does this for us\n",
    "\n",
    "\n",
    "- Given the $x$ and $y$ coordinates along the $x$ and $y$ axis, `meshgrid` \"fills in the blanks\" and returns 2 matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "x = np.linspace(-2, 2, 5)\n",
    "y = np.linspace(-2, 2, 5)\n",
    "\n",
    "Xm, Ym = np.meshgrid(x, y)\n",
    "Zm = np.sin(Xm) + np.cos(Ym)\n",
    "\n",
    "plt.contourf(Xm, Ym, Zm, 10)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "x = np.linspace(-2, 2, 5)\n",
    "y = np.linspace(-2, 2, 5)\n",
    "\n",
    "Xm, Ym = np.meshgrid(x, y)\n",
    "Vx = 2 + np.zeros((5, 5))\n",
    "Vy = 2 + np.zeros((5, 5))\n",
    "\n",
    "print(Vx)\n",
    "plt.quiver(Xm, Ym, Vx, Vy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Example - Potential Flow Around a Cylinder\n",
    "\n",
    "- Suppose you have a fluid flow field $\\to$ E.g. fluid moving around on object $\\to$ and you want to visualise this flow field at discrete points\n",
    "\n",
    "\n",
    "- For example $\\to$ the potential flow around a cylinder is given in polar coordinates by:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    V_r &= U \\left(1-\\frac{R^2}{r^2}\\right)\\cos\\theta \\\\\n",
    "    V_\\theta &= -U \\left(1+\\frac{R^2}{r^2}\\right)\\sin\\theta\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $R$ is the cylinder radius and $r$ is the distance from the cylinder centre"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Outcomes:\n",
    "\n",
    "- `numpy` module functions:\n",
    "    - `linspace`\n",
    "    - `arange`\n",
    "    - `meshgrid`\n",
    "\n",
    "\n",
    "- Quiver plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "def polar_coordinates(x, y):\n",
    "    r = (x**2 + y**2)**0.5\n",
    "    #theta = np.arctan(y / x)\n",
    "    theta = np.arctan2(y , x)\n",
    "    return r, theta\n",
    "\n",
    "\n",
    "def polar_potential_flow_velocity(r, theta, U, radius):\n",
    "    Vr =  U * (1 - (radius**2 / r**2)) * np.cos(theta)\n",
    "    Vt = -U * (1 + (radius**2 / r**2)) * np.sin(theta)\n",
    "    return Vr, Vt\n",
    "\n",
    "\n",
    "#def cartesian_velocities(Vr, Vt, theta):\n",
    "def cartesian_velocities(Vr, Vt, Theta):\n",
    "    Vx = Vr * np.cos(Theta) - Vt * np.sin(Theta)\n",
    "    Vy = Vr * np.sin(Theta) + Vt * np.cos(Theta)\n",
    "    return Vx, Vy\n",
    "\n",
    "\n",
    "def remove_in_cylinder(Vx, Vy, R, radius):\n",
    "    for i in np.arange(0, num, 1):\n",
    "        for j in np.arange(0, num, 1):\n",
    "            if R[i, j] < radius:\n",
    "                Vx[i, j] = 0.0 #np.nan\n",
    "                Vy[i, j] = 0.0 #np.nan\n",
    "\n",
    "\n",
    "def plot_cylinder(radius):\n",
    "    xvals = np.linspace(-radius, radius, 100)\n",
    "    yvals = (radius**2 - xvals**2)**0.5\n",
    "    plt.plot(xvals,  yvals, 'b-')\n",
    "    plt.plot(xvals, -yvals, 'b-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#%matplotlib inline\n",
    "#%matplotlib notebook\n",
    "%matplotlib qt5\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "U = 20\n",
    "num = 30\n",
    "radius = 1.0\n",
    "\n",
    "# create x, y points\n",
    "x = 1.0*np.linspace(-3, 3, num)\n",
    "y = 1.0*np.linspace(-3, 3, num)\n",
    "\n",
    "# create X, Y meshgrid\n",
    "X, Y = np.meshgrid(x, y)\n",
    "\n",
    "# convert X, Y meshgrid to polar R, Theta meshgrid\n",
    "R, Theta = polar_coordinates(X, Y)\n",
    "\n",
    "# Calculate Vr, Vt polar velocities\n",
    "Vr, Vt = polar_potential_flow_velocity(R, Theta, U, radius)\n",
    "\n",
    "# Convert Vr, Vt polar valocities to Vx, Vy Cartesian velocities\n",
    "Vx, Vy = cartesian_velocities(Vr, Vt, Theta)\n",
    "\n",
    "# Remove vector components inside the cylinder\n",
    "remove_in_cylinder(Vx, Vy, R, radius)\n",
    "            \n",
    "# Plot the cylinder\n",
    "plot_cylinder(radius)\n",
    "\n",
    "# Plot the vectors using a quiver plot\n",
    "quiver = plt.quiver(X, Y, Vx, Vy)\n",
    "key = plt.quiverkey(quiver, 0.9, 1.02, U, 'Key: 10 [m/s]')\n",
    "\n",
    "plt.title(\"Velocity Field\")\n",
    "plt.xlabel(\"Distance [m]\")\n",
    "plt.ylabel(\"Height [m]\")\n",
    "plt.axis('equal')\n",
    "\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "``meshgrid`` graphical coordinate picture of storing values in a matrix - x & y stored separately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example - 3D Visualisation of a function\n",
    "\n",
    "- Given the following function: $$z = \\sin(2x) + \\cos(y)$$ where $x$ is between 0 and $2\\pi$; and $y$ is between $-\\pi$ and $\\pi$.\n",
    "\n",
    "\n",
    "- Write a python program that visualises this function as a 3D surface plot\n",
    "\n",
    "\n",
    "- The visualisation must also include a 2D contour projection onto the x-y plane\n",
    "\n",
    "\n",
    "### Outcomes:\n",
    "\n",
    "- 3D surface plot + optional inputs (``rstride``, ``cstride``, ``cmap``)\n",
    "\n",
    "\n",
    "- projected 2D contour plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import cm\n",
    "from mpl_toolkits.mplot3d import axes3d\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "def plot_3D(X, Y, Z, **kwargs):\n",
    "    ax = plt.gca(projection='3d')\n",
    "    ax.plot_surface(X, Y, Z, **kwargs)\n",
    "    return ax\n",
    "\n",
    "def project_3D_contour(X, Y, Z, **kwargs):\n",
    "    ax = plt.gca(projection='3d')\n",
    "    cset = ax.contourf(X, Y, Z, **kwargs)\n",
    "    return cset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "x = np.linspace(0, 2*np.pi, 30)\n",
    "y = np.linspace(-np.pi, np.pi, 30)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "\n",
    "Z = np.sin(2*X) + np.cos(Y)\n",
    "\n",
    "ax = plot_3D(X, Y, Z, rstride=10, cstride=10)\n",
    "project_3D_contour(X, Y, Z, zdir='z', offset=-2.5, cmap=cm.coolwarm)\n",
    "ax.set_xlabel(r'$x [rad]$')\n",
    "ax.set_ylabel(r'$y [rad]$')\n",
    "ax.set_zlabel(r'$f(x, y)$')\n",
    "ax.set_zlim(-2.5, 1.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example - 3D Representation of a geometry object\n",
    "\n",
    "- Simple geometry objects can be expressed as 1 or more 2 variable mathematical functions, thus visualising these simple geometry objects is similar to the previous example.\n",
    "\n",
    "\n",
    "- Write a python program the visualises a simple sphere (`R = 2`).\n",
    "\n",
    "\n",
    "### Outcomes:\n",
    "\n",
    "- 3D surface plot + optional inputs (``rstride``, ``cstride``, ``color``)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# create supporting points in sherical coordinates\n",
    "u = np.linspace(0, 2*np.pi, 100)\n",
    "v = np.linspace(0, np.pi, 100)\n",
    "\n",
    "# transform them to cartesian system\n",
    "x = 10*np.outer(np.cos(u), np.sin(v))\n",
    "y = 10*np.outer(np.sin(u), np.sin(v))\n",
    "z = 10*np.outer(np.ones(np.size(u)), np.cos(v))\n",
    "\n",
    "ax = plot_3D(x, y, z, rstride=5, cstride=5, color='g')\n",
    "ax.set_xlabel(r'$x [m]$')\n",
    "ax.set_ylabel(r'$y [m]$')\n",
    "ax.set_zlabel(r'$f(x, y)$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Publication quality plots\n",
    "\n",
    "- For creating publication quality plots please see the textbook and the following links\n",
    "    - https://www.bastibl.net/publication-quality-plots/\n",
    "    - https://python4mpia.github.io/plotting/publication.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}