Practical_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Practical 1) Creating binned size distributions and calculating PM2.5\n",
    "\n",
    "In our first practical, we are going to practice creating discrete size distributions before then calculaing the total mass loading for a range of ambient environments. Creating a size distribution can often feel like an ardous task and is a useful exercise in translating equations to code. This will also serve as a useful warm-up exercise for deploying a range of Python functions. You will have looked at Python as part of your training already; However learning to code wholst focusing on a domain relevant problem is often the best way to learn. With that in mind, today's first practical is comprised of 4 exercises across 3 main objectives.\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Objectives:</b> The objectives of todays practical are:\n",
    "\n",
    " - 1) [Objective 1 - Practice changing the shape of a normal and log-normal distribution ](#Part1)\n",
    "      * [Exercise 1: Create 3 normal distributions and plot the results](#Exercise1)\n",
    "      * [Exercise 2: Create 3 log-normal distributions and plot the results](#Exercise2)\n",
    " - 2) [Objective 2 - Populate a discrete 'binned' size distribution from a continuous function](#Part2)\n",
    "      * [Exercise 3: Create a discrete size distribution with 30 bins](#Exercise3)\n",
    " - 3) [Objective 3 - Calculate PM2.5 mass](#Part3)\n",
    "      * [Exercise 4: Sum to the total particulate mass for a range of environments](#Exercise4)\n",
    "\n",
    "Please note that you should not feel pressured to complete every exercise in class. These practicals are designed for you to take outside of class and continue working on them if useful. Proposed solutions to all exercises can be found in the 'Solutions' folder. There are often many ways to arrive at the solution of any particular problem, where choice of final method can rest on computational efficiency. We will discuss this later.\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>Please note:</b> After reading the instructions and aims of any exercise, if some boilerplate is already provided, search the code snippets for a note that reads ------'INSERT CODE HERE'------ to identify where you need to write your code \n",
    "</div>\n",
    "\n",
    "## Objective 1 - Practice changing the shape of a normal and log-normal distribution  <a name=\"Part1\">\n",
    "\n",
    "The probability density of many phenomena in the natural world often follows a normal, or Gaussain, distribution. Algebraically, the probability density for the normal (Gaussian) distribution is given by:\n",
    "\n",
    "\\begin{equation}\n",
    "p(x) = \\frac{1}{\\sqrt{ 2 \\pi \\sigma^2 }} e^{ - \\frac{ (x - \\mu)^2 } {2 \\sigma^2} },   \\tag{1}\n",
    "\\end{equation}\n",
    "\n",
    "where $p(x)$ is the probability density, $\\sigma$ is the standard deviation and $\\mu$ is the mean. In some cases you may hear reference to $\\mu$ as the location parameter and $\\sigma$ as the scale parameter. Indeed, if you were to use the Python package Scipy's ```python scipy.stats.lognorm ``` function you would find such reference in the documentation. We are going to create our own normal and log-normal distributions.\n",
    "\n",
    "A lognormal distribution is a normal distribution on a logarithmic scale. For aerosol particles that can span large ranges in magnitude, this is the most appropriate and perhaps commonly used distribution. Distinct from a normal distribution, a lognormal distribution has \"geometric\" mean diameter and a \"geometric\" standard deviation. As noted by [Jacobson 2005]() fifty percent of the area under a lognormal distribution lies below the geometric-mean diameter.  and 68 percent of the area under the lognormal distribution lies within one geometric standard deviation of the geometric-mean diameter.\n",
    "\n",
    "The probability density for the log-normal distribution is\n",
    "\\begin{equation}\n",
    "p(x) = \\frac{1}{\\sqrt{ 2 \\pi log(\\sigma_{g})^2 }} e^{ - \\frac{ (log(x) - log(\\mu))^2 } {2 log(\\sigma_{g})^2} },\n",
    "\\end{equation}\n",
    "\n",
    "A standard normal distribution can be visualised as per the following figure:\n",
    "\n",
    "<img src=\"images/Normal_Distribution_PDF.png\" alt=\"By Inductiveload - self-made, Mathematica, Inkscape, Public Domain, https://commons.wikimedia.org/w/index.php?curid=3817954\" style=\"width: 500px;\"/>\n",
    "\n",
    "\n",
    "\n",
    "This brings us to exercise 1. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<b> Exercise 1. Create 3 normal distributions and plot the results <a name=\"Exercise1\"> </b> \n",
    "    \n",
    "Using the defition of a normal distribution, you are tasked with now creating three different distributions and plotting the results. We are going to assume that each distribution has the same mean of 0 but 3 different standard deviations of 0.3, 1.0 and 1.3.\n",
    "\n",
    "We are going to use the Numpy and Matplotlib libraries to create our arrays and plot results accordingly. As already noted above, this warm-up exercise will help you get a renewed feel for Python programming. You may find the following functions useful, assuming we have imported ```python numpy as np```:\n",
    "\n",
    "```python\n",
    "   np.linspace(lower, upper, num=xxx) # Creates an array with a lower and upper limit with xxx entries\n",
    "   np.exp(xxx) # Calculates the exponential of a constant, variable or Numpy array\n",
    "   np.sqrt(xxx) # Calculates the square root of a constant, variable or Numpy array\n",
    "   np.pi # Gives us the constant 'pi'\n",
    "\n",
    "```\n",
    "\n",
    "In the code snippet below, as part of this exercise, I have imported the relevant libraries for you and also provided a plotting function. Notice I have assumed our x variable is called 'x' and our 3 different distributions are called 'distribution_1','distribution_2' and 'distribution_3'.  Once complete you should arrive at a plot that looks like the following:\n",
    "\n",
    "![title <](images/Prac1_ex1.png)\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU9b3/8dcnC6vsBGUTULEqXheMiIotiqBSK9a6gK1L3bVY71Vs4UeLreLaXq1t3eJ6cUO9YkWvgBa1tCpIcEHEIhFEIggogRC2JMzn98eZgSHMJJNksszM+/l4zGNmzvnOOZ85DO/55jtnMXdHRERSX1ZTFyAiIsmhQBcRSRMKdBGRNKFAFxFJEwp0EZE0kdNUK+7atav37du3qVYvIpKSFixY8K2758Wa12SB3rdvXwoLC5tq9SIiKcnMVsSbpyEXEZE0oUAXEUkTCnQRkTShQBcRSRMKdBGRNKFAFxFJEwp0EZE0oUAXifbCCzB0KNx2W1NXIlJr1lTnQ8/Pz3cdWCTNzoABsHhx8HjLFmjdumnrEanCzBa4e36seeqhi0SEQvDFF7ueL1vWdLWI1IECXSRi1SrYvn3X8+hwF0kBCnSRiKKi6p+LNHMKdJEIBbqkOAW6SERkzHzQoOBeQy6SYhToIhFr1wb3RxwR3K9b13S1iNSBAl0k4rvvgvsDD9z9uUiKUKCLRKxfH9z377/7c5EUoUAXiYj0yPv0gZwcKCuD8vKmrUmkFhToIhGRQO/aFTp33n2aSApQoIsAuO8K786doUuX4LECXVKIAl0EYPNmqKgIzt3SuvWuHrrG0SWFKNBFYFdPPNIzVw9dUpACXQQU6JIWFOgisGtoJTLUoiEXSUE1BrqZPWZma81sUZz5PzWzheHbu2Z2ePLLFGlgJSXBfadOu99HpoukgER66E8Ap1YzfznwA3c/DLgFKEhCXSKNa9Om4L59+93vI9NFUkBOTQ3cfY6Z9a1m/rtRT+cCvepflkgjiwR3u3a73yvQJYUkewz9UmBGvJlmdoWZFZpZ4Tqd+Eiak9LS4L5qoEemi6SApAW6mZ1IEOi/jtfG3QvcPd/d8/Py8pK1apH605CLpIEah1wSYWaHAY8Ap7m79vOS1KMeuqSBevfQzWxfYBpwgbt/Xv+SRJqAeuiSBmrsoZvZs8BQoKuZFQM3AbkA7v4gMAnoAtxvZgCV7p7fUAWLNAj9KCppIJG9XMbUMP8y4LKkVSTSFDTkImlAR4qKwJ5DLpFALyuDUKhpahKpJQW6COzZQ8/OhjZtgtPqbt7cdHWJ1IICXQT2HEOPfqxxdEkRCnQR2NVDjwy5RD/WOLqkCAW6SGUlbN0KZtC27a7p6qFLilGgi5SVBfd77RWEeoT2dJEUo0AXifzouddeu0+PPN+ypXHrEakjBbpIJLDbtNl9euS59nKRFKFAF4kX6JHxdPXQJUUo0EXUQ5c0oUAXiQS2euiS4hToIpHAjt5lEdRDl5SjQBfRGLqkCQW6SE2Brh66pAgFukhNP4qqhy4pQoEuUtOPouqhS4pQoIvU9KOoeuiSIhToIhpDlzShQBfRGLqkCQW6iHrokiYU6CLqoUuaqDHQzewxM1trZovizDcz+7OZFZnZQjMbmPwyRRqQ9nKRNJFID/0J4NRq5p8G9A/frgAeqH9ZIo1Ie7lImqgx0N19DrC+miajgCkemAt0NLPuySpQpMElMuQSCjVuTSJ1kIwx9J7AyqjnxeFpezCzK8ys0MwK161bl4RViyRBvEDPyoJWrYLH27Y1bk0idZCMQLcY0zxWQ3cvcPd8d8/Py8tLwqpFkiBeoIPG0SWlJCPQi4HeUc97AauSsFyRxlFdoGscXVJIMgJ9OnBheG+XwcBGd1+dhOWKNI5I77vqj6LR09RDlxSQU1MDM3sWGAp0NbNi4CYgF8DdHwReA0YCRcAW4OcNVaxIg1APXdJEjYHu7mNqmO/AL5JWkUhjqqiAykrIzobc3D3nq4cuKURHikpmi+6dW4zf93XVIkkhCnTJbNUNt0RPVw9dUoACXTJbvMP+I9RDlxSiQJfMFu+w/wj10CWFKNAls9U05KIeuqQQBbpkNo2hSxpRoEtmUw9d0ogCXTKbeuiSRhToktmqO+w/erp66JICFOiS2dRDlzSiQJfMpjF0SSMKdMls6qFLGlGgS2ZLtIeuQJcUoECXzFbTof+R6Vu3Nk49IvWgQJfMVtOh/xpDlxSiQJfMpjF0SSMKdMlsiQa6euiSAhToktkU6JJGFOiS2WoK9Nzc4PJ0FRXBTaQZU6BLZqvp0H8z7ekiKSOhQDezU81siZkVmdn4GPP3NbO3zOxDM1toZiOTX6pIA6iphx49T8Mu0szVGOhmlg3cB5wGHAKMMbNDqjT7DfC8ux8JjAbuT3ahIg0ikUDXwUWSIhLpoQ8Citx9mbuXA1OBUVXaONA+/LgDsCp5JYo0IPXQJY3kJNCmJ7Ay6nkxcEyVNr8DXjeza4G2wMlJqU6koSnQJY0k0kO3GNO8yvMxwBPu3gsYCTxpZnss28yuMLNCMytct25d7asVSabKSigvD374bNkyfjsFuqSIRAK9GOgd9bwXew6pXAo8D+Du7wGtgK5VF+TuBe6e7+75eXl5datYJFmiD/u3WP2WMAW6pIhEAn0+0N/M+plZC4IfPadXafMVMAzAzA4mCHR1waV5S2S4JXq+Al2auRoD3d0rgbHALOAzgr1ZPjWzm83sjHCzG4DLzexj4FngYnevOiwj0rwkGujay0VSRCI/iuLurwGvVZk2KerxYuD45JYm0sDUQ5c0oyNFJXMp0CXNKNAlc9V0cYsIBbqkCAW6ZK6aLm4RoUCXFKFAl8ylIRdJMwp0yVzay0XSjAJdMpd66JJmFOiSuRTokmYU6JK5arq4RYQCXVKEAl0yl3rokmYU6JK5FOiSZhTokrkU6JJmFOiSubTboqQZBbpkLh36L2lGgS6ZS4f+S5pRoEvmSnTIpWXL4IpG27fDjh0NX5dIHSnQJXMlGuhmu9ps3dqwNYnUgwJdMleigR7dRsMu0owp0CVz1SbQtaeLpAAFumSuRA/9B/XQJSUo0CVzachF0owCXTLTjh3BXisArVrV3F6BLikgoUA3s1PNbImZFZnZ+DhtzjWzxWb2qZk9k9wyRZIssrdKmzbBXiw1UaBLCsipqYGZZQP3AcOBYmC+mU1398VRbfoDE4Dj3b3EzLo1VMEiSVGb4Zbodgp0acYS6aEPAorcfZm7lwNTgVFV2lwO3OfuJQDuvja5ZYokWaKH/UdoLxdJAYkEek9gZdTz4vC0aAcCB5rZO2Y218xOjbUgM7vCzArNrHDdunV1q1gkGSLBvNdeibVXD11SQCKBHmuA0as8zwH6A0OBMcAjZtZxjxe5F7h7vrvn5+Xl1bZWkeQpKwvuE9llERTokhISCfRioHfU817AqhhtXnb3CndfDiwhCHiR5kk9dElDiQT6fKC/mfUzsxbAaGB6lTZ/A04EMLOuBEMwy5JZqEhS1eagIlCgS0qoMdDdvRIYC8wCPgOed/dPzexmMzsj3GwW8J2ZLQbeAm509+8aqmiRetOQi6ShGndbBHD314DXqkybFPXYgevDN5Hmr65DLtrLRZoxHSkqmam2PfRIO/XQpRlToEtm0o+ikoYU6JKZNIYuaUiBLplJPXRJQwp0yUzabVHSkAJdMpOGXCQNKdAlM9V2yEUn55IUoECXzKQeuqQhBbpkJv0oKmlIgS6ZqbY/ikYuU7d1K4RCDVOTSD0p0CUzRYZcEu2hZ2VB69bB423bGqYmkXpSoEtmqm0PHXQ+F2n2FOiSmWrbQ49uG3mtSDOjQJfMU1ER3LKzoUWLxF+nQJdmToEumSd6uMViXWExjnbtgnsFujRTCnTJPHUZbolur0CXZkqBLpmnLj+IggJdmj0FumQe9dAlTSnQJfOohy5pSoEumUc9dElTCnTJPPXtoW/alNx6RJIkoUA3s1PNbImZFZnZ+GranW1mbmb5yStRJMk05CJpqsZAN7Ns4D7gNOAQYIyZHRKjXTvgl8C8ZBcpklR1HXLRfujSzOUk0GYQUOTuywDMbCowClhcpd0twF3AuKRWKJJscXroy5bBSy/Bl19CXh6MHAn50X9rqocuzVwigd4TWBn1vBg4JrqBmR0J9Hb3V80sbqCb2RXAFQD77rtv7asVSYYqPfTt22HcOLj//t3PjHvTTfDDH8LjjwcBr0CX5i6RMfRYx0b7zplmWcA9wA01LcjdC9w9393z8/LyEq9SJJmieujr18OIEfDXvwZnAfjpT+Gee2DsWGjfHv7v/2DwYPj3v1GgS7OXSA+9GOgd9bwXsCrqeTvgUOBtC86LsQ8w3czOcPfCZBUqkjThQN/esj0/+hG8+y706AGvvAIDB+5qNmECjBoFhYUwbBjML+hED1CgS7OVSA99PtDfzPqZWQtgNDA9MtPdN7p7V3fv6+59gbmAwlyar7IyHLjq+ZN4913o1Qvmzt09zCEI+bffhiFDYNUqOHP8QWyjpXZblGarxkB390pgLDAL+Ax43t0/NbObzeyMhi5QJOlKS3mO83jin/vTpg1Mnw69e8du2rYtTJsGffrA/EWtmcTNUFrauPWKJMjcveZWDSA/P98LC9WJl8a39thRHDL3Ub6jKwUFcPnlNb/m/ffh2GMdQiHebX0yx2x5q+ELFYnBzBa4e8xjfXSkqGSccUsu5zu6MmxQKZddlthrBg2CG66HENlcvvVedpTvaNgiRepAgS4ZpbAQniw5nRZsp+COklpd3+L3Nxt9bQWfcBiPP7i94YoUqSMFumQMd7ghvHPtddzLfoe2qdXrW7eG2zvdBcBvJ7fUzi7S7CjQJWO8/jrMmQNd+Jb/x23Bjua1dN7eb3MMc/lmXTZ//nMDFClSDwp0yQjucMstweNfcRcdW2yFli1rvRzr2IHJ/AYIDkCKHKMk0hwo0CUjvP02vPMOdO4U4moeqFPvHID27RnGbAYfVMK338KDDya1TJF6UaBLRpg8Obj/zwtLaEcZdOhQtwW1b48BvznjEwD+8AfYujU5NYrUlwJd0t6778Kbbwad8mvPDJ9nLqqHvm3bNqZOncqVV17Jsccey4ABAzj66KMZM2YMjzzyCOvWrdu1sPAXwcj9lzBwIKxZA4880pjvRiQ+BbqkvbvvDu7HjoWOXhI8ad+eLVu2MHnyZHr06MGYMWMoKChg7ty5LF68mMLCQqZOncrll19O7969ueaaa1i9evXOLwIr3cjEicGi7rkHdmi3dGkGEjk5l0jKWrEiOMd5Tk4Q6LwfHLY/p6KCCw85hBUrVgAwcOBARo8ezaBBg+jSpQulpaUsXLiQl19+mZkzZ/LAAw/wzDPPcPcPfsDPASstZdQo6NcPli8PTux15plN9z5FAHD3JrkdddRRLtLQfvUrd3A///zgeeh//sdvBc8yc8APP/xwf/PNN6tdxqeffuojR450gtNG+8XgW37xC3d3v+eeYPlDhzb0OxEJAIUeJ1cV6JK2Nm9279Qp+JTPneteUVHhlwwe7IAb+MSJE728vDyhZYVCIZ8yZYq3btHCAf/+3nt7aWmpb9zo3q5dsI4PP2zgNyTi1Qe6xtAlbT39NJSUBOdhyc/fwc9//nMemzuX1sC0s89m8uTJ5ObmJrQsM+OCCy7gvdtvpycwZ80ahg8fzo4dJVxySdDm3nsb7K2IJESBLmnJnZ1Hco4dG+Kyyy7jqaeeYq/cXN4Azjz22Dot9/CBA5kD9G3Zknnz5jF8+HAuuaQMM3jmmWCvF5GmokCXtPTWW7BoEeyzDyxYcCNPPPEEbdq04bWTTuJ4gE6d6rbgzp3ZD5jTpw/77bcfCxYs4Ne/Pocf/rCC8nIdaCRNS4EuaSnSO8/Pf4h7772b3Nxcpk+fzgktWgQz6hro4df1Litj1qxZdO3alZkzZ7Jjx1WA88ADUF5e//pF6kKBLmln+fLgKkTZ2X9nxoxfAFBQUMCwYcOCQXWAzp3rtvDIF0FJCQcccACvvvoqrVu3ZsaMx9hnn3tZswZeeCEJb0KkDhToknbuuw/cvyYnZzQ7duxg/PjxXHzxxcHM9euD+7r20Nu2hdzc4Hj/bds45phjmDJlCgBr144D/qmzMEqTUaBLWikrg4cf3gH8jO3bv2PEiBHceuutuxrUt4dutlsvHeDss89m3LhxhEI7MDuX999fzbx5dX8PInWlQJe08uSTUFp6G/A2e++9N1OmTCErK+pjHgn0uvbQo18bWRZw++23M3ToUNy/Ac7hT3+qqPvyRepIgS5pIxSCO+74J/A7AJ588kn23nvvXQ3CwyS0aBFcfqiuIoEeGb4BcnJyeO6559hnn57AOzz33CRWrar7KkTqIqFAN7NTzWyJmRWZ2fgY8683s8VmttDMZptZn+SXKlK9adPW89VX5wMhbrxxPMOHD9+9QfRwS20uJlpVZLgmqocO0K1bN154YSqQhfudjB//Zt3XIVIHNQa6mWUD9wGnAYcAY8zskCrNPgTy3f0w4H+Bu5JdqEh13J1f/OISoJjevQdz660379movj+IRkQCPaqHHjFkyBAuvvi3gPP00z/j66+/rd+6RGohkR76IKDI3Ze5ezkwFRgV3cDd33L3LeGnc4FeyS1TpHqTJt3H2rUvAx14+eVnYx/SHzmveV5e/VbWtevuy6uioOA3tGlzPKHQak4//ZLgpEkijSCRQO8JrIx6XhyeFs+lwIxYM8zsCjMrNLPCdXH+M4jU1kcffcRtt90AwMknP8yRR/aN3XDt2uC+W7f6rTDy+sjyqsjNzeGmm54GOvDRR69w33331299IglKJNBjDTbG7HKY2c+AfOAPsea7e4G757t7fl59e0kiQFlZGeeeO5pQqBy4kj/96Zz4jZMd6NV0Sq69tg977fUwADfccAOffPJJ/dYpkoBEAr0Y6B31vBewx+/3ZnYyMBE4w923J6c8kepde+21LF26BBjASSfdw4AB1TRupB46BDvRXHvtOcCllJdvZ/To0WzZsiVue5FkSCTQ5wP9zayfmbUARgPToxuY2ZHAQwRhHv9TLpJETz31FE888QRmrYHnuP76GnZFbMRAB7jmGsjKuhf4HosXL+b666+v33pFalBjoLt7JTAWmAV8Bjzv7p+a2c1mdka42R+AvYAXzOwjM5seZ3EiSbF06VKuvvpqANzv5eCDB3DaaTW8qJEDvVcvGD26LTCVrKwWPPTQQ7z44ov1W7dINRLaD93dX3P3A919f3e/NTxtkrtPDz8+2d33dvcjwrczql+iSN1t3x4MYZSVldGu3bnAZYwfD1k1fZojY97JDPQa9mAZPx7gCLKzg5+VLrvsMr766qv6rV8kDh0pKiln3LhxfPDBB3Tr1pdNmx6iTx9jzJgEXpisHnrbttCmTXDUaVlZtU3/4z/g9NOhouJaDjjgdDZs2MD5559PZWVl/WoQiUGBLinlhRde4K9//Su5ubl06fI80JFx44ITINbom2+C+/oGevQyErhE0YQJAMaaNY/TvXsP3nnnHW655Zb61yBShQJdUkZRURGXXnopAJdc8t989tnR5OWx85qe1SothU2bgp51x471L6ZX+Ni54uIamx53HHz/+7BpU1dOOeUpzIzJkyfz1ltv1b8OkSgKdEkJ27Zt49xzz2XTpk2cddZPeO+9sUAwRt2mTQILWBk+Nq5Xr/qdxyUiEugrV1bfLuz3vw/up007keuvn0goFOK8885jZYKvF0mEAl2aPXfnuuuu48MPP2S//fZj5MhHWbjQ6NUr2DUwIZGedO/e1bdLVGQ5CfTQAYYOhREjgj8U3H/H8OHDWbduHT/5yU/Ytm1bcmqSjKdAl2bv/vvvp6CggJYtW/LMM89z++0dALjpJmjVKsGFRPfQk6GWPXSA224L7u+/P5u7736WPn36MH/+fH75y18mpybJeAp0adZmz57NddddB8Cjjz7K/PlH8cUXcOCBELmqXEIiwZvsHnotAv2oo+Dss4OdY+6+uwvTpk2jVatWPPzww9x/v873IvWnQJdma+nSpZxzzjns2LGDCRMmcPLJP+U3vwnm3XEH5OTUYmHNINABbr012CPn8cdh27aBFBQUAMEpDF555ZXk1CYZS4EuzdKaNWsYOXIkJSUljBo1ismTJ3PjjbBxI5x2Gpx5Zi0XuGJFcJ+sQN933+D+yy9rPLgo2oEHwq9+FTy++moYM+YCbrrpJkKhEKNHj2b+/PnJqU8yk7s3ye2oo45ykVg2bNjgRx55pAN+5JFHemlpqf/jH+7g3rKle1FRHRbas2ewgDq9OIZQyL19+2CZa9bU6qWbN7v36RO89E9/cg+FQn7xxRc74N26dfOiZNUoaQko9Di5qh66NCubN2/mjDPO4MMPP6R///7MnDmTrKx2XH55MH/CBNh//1outKwMvv46uJZo377JKdQMvve94PGSJbV6aZs28Je/BI8nToRly4yCggKGDx/O2rVrOfHEE1m2bFly6pSMokCXZqO0tJRTTz2VOXPm0KNHD15//XW6devGuHHw+ecwYAD8+td1WPDSpcH9AQdAdnbyCo4E+uef1/qlP/oRnHsubN4MF1wAZrm8+OKLHHfccaxcuZITTzyRL7/8Mnm1SkZQoEuzsH79ekaMGMG//vUvevXqxdtvv03fvn158UV48MGgc/3007XYTTFapAd94IFJrXnn8mrZQ4944AHo2RPeey/YBbNdu3bMmDGDwYMH89VXXzFkyBAWLVqUxIIl3SnQpcktXbqUwYMHM2/ePPr27cucOXPo378/n34KF10UtLnzTjj88DquIBKKBx2UlHp3Ovjg3ZdfS507w5NPBmeJvO02mDYN2rdvz8yZMxkyZAhff/01Q4YM4R//+EcSi5Z0pkCXJjV79mwGDx7M0qVLOfzww/nnP/9Jv379WLUqOEvh5s1w/vkQ3hW9biJ7jhx9dFJq3ik/P7h///1a7ekS7cQT4a67gscXXADz5kGHDh144403OOuss9i4cSPDhw/nvvvu08WmpWbxfi1t6Jv2cslsFRUVPnHiRDczB/z000/3TZs2ubv7t9+6H3posBfIoEHBXiF1Fgq5d+wYLGzlyuQUH73svLxg2V98Ua/FXHRRsJjOnd0/+SSYXllZ6f/1X//lBNfw9fPPP983btyYnNolZVHNXi4KdGl0H3/8sR9zzDEOeFZWlk+aNMkrKyvd3b242P2QQ4JP5sEHB+FeL0uWBAvr3j1IzmT74Q+D5T/9dL0WU17ufvrpu0J93rxd86ZOnept27Z1wHv16uUzZsyoZ9GSyqoLdA25SKMpKSnhxhtvZODAgcybN4+ePXsye/Zsfv/735Odnc28eXDMMbB4cbBHy9//Dl261HOlM2cG9yeckJyzLFZ1wgnB/axZ9VpMbi48/3yw98v69cFQzNSpwbzzzjuP+fPnc/TRR1NcXMxpp53G6NGjtWuj7Cle0jf0TT30zFFSUuKTJ0/2jh07OuBm5mPHjvUNGza4u3tFhfsdd7i3aBH0UI8/Pgk984hhw4KFPvVUkhZYxeLFwfK7dAneSD2Vl7tffHGwSHC//HL3yChLRUWF33nnnd6qVSsHPDc316+99lpfvnx5vdcrqQMNuUhjC4VCvmDBAr/qqqt2DhcAftJJJ/n7778fbuM+fbr7YYftCrCxY923b09SEStWuGdnB7fvvkvSQqsIhdz79w+Kf+WVpC3yL39xz80NFtujh/uDDwZh7+6+YsUKv/DCC3f+/pCVleXnnHOOz5w5c+fQlaQvBbo0iq1bt/rs2bN9woQJ3r9//50hDviwYcN89uzZ7u5eWur+yCPuRx+9K8j33dd95swkF/TLXwYLP//8JC+4ij/+MVjPkCFJHaf/5JPgR+HINurXz/3OO91XrQrmL1y40C+44ALPycnZuZ27d+/uV1xxhb/00kv6ATVNVRfoFsyvnpmdCtwLZAOPuPsdVea3BKYARwHfAee5+5fVLTM/P98LCwsTHxuSZsPd2bBhA8uXL2fRokV89NFHfPDBB8ybN2+3izV069aNc889l4suugr3Abz3Hrz+Orz1FmzZErTJywsOf7/yyjoeNBTP3/8Op5wSZOEHH8ARRyRx4VWUlkK/fsHg9x//CDfckLRFh0Lw4ovw29/uOn4pOxuOPz64YMYJJ0CnTiv529+eYMqUKRQVFe18bU5ODocddhj5+fnk5+dz0EEHsf/++9O9e3esIX5PkEZhZgvcPT/mvJoC3cyygc+B4UAxMB8Y4+6Lo9pcAxzm7leZ2Wjgx+5+XnXLrWugr/p8FSXfbAAI/4kRTPfwg0h/JiK0I7TzsXukVdAw+Fbb9Zzw86hWhELRj3ct13EIRdYZve5dBXjUvOD1uy1g51o8VOX1MdvvqtU9FF6E7Wq/c7VVat9ZS/TrI9MNx9levp3y8u1sK9/O9u3bdj7fsm0LG8s2sLF0Axs3lbBh0wbWrV/DqrVfU7ZlE7H07DaAnnt/n47tR+AM5cviNhR9lYv77gFywlFbuOSsDZxzSiltW4diLot4n8140ysrg7Mqzp4NDz0E5eXBuQLuuCN2+2R64YXgWH4ITno+alRwJGnHjsHJW8x23SD+4zgqK2HG7BY8/mwrXnm9BZWVu7fv1WMH/ferpGWLBWwqm8XXq2exonj+zs9KtFat2tBzn9507ZxHl05d6NypC107d6Fjhw60ad2a1q1a06pVcN+6VWtatGhJdnYWWVnZZGdnk5OVRVZ2NtlZ2WRnZ5GdlY1lZZFlWTu/KKLfTqwvj3hv2bCqE+J++cRebmp8UbVq24p+R/St02urC/Qah0aAY4FZUc8nABOqtJkFHBt+nAN8S/jLIt6trkMuh3Y+c7c/5XVrqtteDgMcznaY7PCqw1rf+c0Ydcuh3Afwif+MKf44F3kxPfZs1BC3665zb8wx5QcfDE4H2cDvq4QOPo0z/Sru96OZ563ZHKdpqcM/HP7b4acOgxy6NIPPjm5tsw6r88eMaoZcErlEQE8g+iz+xcAx8dq4e6WZbQS6EAT7TmZ2BXAFwL6R80nXUpuWrTC6xpy3+7d7vMc1zYtMSWRZez42vBbriykwDaYAAAesSURBVK448qjq62tqH0t0q6rL8j3aZFlLjOCWRUuyrBVGC4xW5NCBbOtINh3IsY7k0IUW1pNsOpKbFaKtbaVt1lbaZm2hrb3DXllb2CfnW7pnr6VH7jp656zmey2W0zKrIqqM9kD7GGXH6rLFeZ/x2vbqBYceGhxeeuSR8TZQw7jyyuBk7c8/D3PmwKpVsGEDbN26K2Mh/uMEdQR+zL/4Mf8CYIdn8UWoL1/u2JdVoX1YHdqb1aG92ejt2extKPMj2OzHscVbU0k25b6FbXxDhW+ggg1UUEKll7CDzTjbcLbudoPtODuAELADD99HbsHzyK2qWO8t9vv1ONNjS3y5zVGr7ESubF57iQR6rP9RVbdcIm1w9wKgAIIhlwTWvYd5q56ty8tEGse++8K4ccGtkWQDB4ZvktkSObCoGIi+zEsvYFW8NmaWA3QA1iejQBERSUwigT4f6G9m/cysBTAamF6lzXTgovDjs4E3w2M9IiLSSGoccgmPiY8l+OEzG3jM3T81s5sJBuenA48CT5pZEUHPfHRDFi0iIntK6Lrp7v4a8FqVaZOiHm8DzkluaSIiUhs6OZeISJpQoIuIpAkFuohImlCgi4ikiYROztUgKzZbB6yo48u7UuUo1GaiudYFzbc21VU7qqt20rGuPu6eF2tGkwV6fZhZocc7OU0Taq51QfOtTXXVjuqqnUyrS0MuIiJpQoEuIpImUjXQC5q6gDiaa13QfGtTXbWjumono+pKyTF0ERHZU6r20EVEpAoFuohImmi2gW5m55jZp2YWMrP8KvMmmFmRmS0xs1PivL6fmc0zs6Vm9lz41L/JrvE5M/sofPvSzD6K0+5LM/sk3K7Br4xtZr8zs6+jahsZp92p4W1YZGbjG6GuP5jZv81soZm9ZGYd47RrlO1V0/s3s5bhf+Oi8Gepb0PVErXO3mb2lpl9Fv78XxejzVAz2xj17zsp1rIaqL5q/20s8OfwNltoZgMboabvRW2Lj8ys1Mz+s0qbRtlmZvaYma01s0VR0zqb2RvhLHrDzDrFee1F4TZLzeyiWG1qFO/adE19Aw4Gvge8DeRHTT8E+BhoCfQDvgCyY7z+eWB0+PGDwNUNXO9/A5PizPsS6NqI2+53wLga2mSHt91+QIvwNj2kgesaAeSEH98J3NlU2yuR9w9cAzwYfjwaeK4R/u26AwPDj9sRXKC9al1DgVcb6/NUm38bYCQwg+AqZoOBeY1cXzbwDcHBN42+zYDvAwOBRVHT7gLGhx+Pj/W5BzoDy8L3ncKPO9V2/c22h+7un7n7khizRgFT3X27uy8HioBB0Q0suPT3ScD/hif9D3BmQ9UaXt+5QCpdH28QUOTuy9y9HJhKsG0bjLu/7u6V4adzCa5+1VQSef+jCD47EHyWhlkDX1be3Ve7+wfhx5uAzwiu2ZsqRgFTPDAX6Ghm3Rtx/cOAL9y9rkeh14u7z2HPq7VFf47iZdEpwBvuvt7dS4A3gFNru/5mG+jViHXR6qof+C7AhqjwiNUmmU4A1rj70jjzHXjdzBaEL5TdGMaG/+R9LM6feIlsx4Z0CUFPLpbG2F6JvP/dLn4ORC5+3ijCQzxHAvNizD7WzD42sxlmNqCxaqLmf5um/lyNJn7Hqqm22d7uvhqCL2ygW4w2SdluCV3goqGY2d+BfWLMmujuL8d7WYxpdbpodSISrHEM1ffOj3f3VWbWDXjDzP4d/iavs+rqAh4AbiF4z7cQDAddUnURMV5b731YE9leZjYRqASejrOYpG+vWKXGmNZgn6PaMrO9gBeB/3T30iqzPyAYUigL/z7yN6B/Y9RFzf82TbnNWgBnABNizG7KbZaIpGy3Jg10dz+5Di9L5KLV3xL8qZcT7lnFapOUGi24KPZZwFHVLGNV+H6tmb1E8Od+vQIq0W1nZg8Dr8aYlch2THpd4R97TgeGeXjwMMYykr69YqjNxc+LrREvfm5muQRh/rS7T6s6Pzrg3f01M7vfzLq6e4OfhCqBf5sG+Vwl6DTgA3dfU3VGU24zYI2ZdXf31eHhp7Ux2hQTjPNH9CL4/bBWUnHIZTowOrwHQj+Cb9n3oxuEg+ItggtWQ3AB63g9/vo6Gfi3uxfHmmlmbc2sXeQxwQ+Di2K1TZYqY5Y/jrO+RC7+ney6TgV+DZzh7lvitGms7dUsL34eHqN/FPjM3e+O02afyFi+mQ0i+H/8XUPWFV5XIv8204ELw3u7DAY2RoYbGkHcv5SbapuFRX+O4mXRLGCEmXUKD5GOCE+rnYb+1bcevxb/mOBbazuwBpgVNW8iwR4KS4DToqa/BvQIP96PIOiLgBeAlg1U5xPAVVWm9QBei6rj4/DtU4Khh4bedk8CnwALwx+m7lXrCj8fSbAXxReNVFcRwTjhR+Hbg1XrasztFev9AzcTfOEAtAp/dorCn6X9GmEbDSH4U3th1HYaCVwV+ZwBY8Pb5mOCH5ePa+i6qvu3qVKbAfeFt+knRO2h1sC1tSEI6A5R0xp9mxF8oawGKsL5dSnB7y6zgaXh+87htvnAI1GvvST8WSsCfl6X9evQfxGRNJGKQy4iIhKDAl1EJE0o0EVE0oQCXUQkTSjQRUTShAJdRCRNKNBFRNLE/wcZ535Bt5CjVwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot 3 different normal distributions with the same mean [0] but 3 different standard deviations of 0.3,1.0 and 1.3\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#------'INSERT CODE HERE'------\n",
    "\n",
    "# Define lower and upper limit of the size distribution [microns]\n",
    "lower_size =  \n",
    "upper_size =  \n",
    "\n",
    "# Use the above information and the fact we want to create 400 values to create a numpy array\n",
    "x = np.linspace(lower_size, upper_size, num=400) # values for x-axis\n",
    "\n",
    "# Define the mean and standard deviation for each normal distribution\n",
    "sigma1 =  \n",
    "mean1 =  \n",
    "\n",
    "sigma2 =  \n",
    "mean2 =  \n",
    "\n",
    "sigma3 =  \n",
    "mean3 =  \n",
    "\n",
    "# Now create each distribution using the mean, standard deviation and x array\n",
    "distribution_1 =  \n",
    "distribution_2 =  \n",
    "distribution_3 =  \n",
    "\n",
    "#------------------------------\n",
    "\n",
    "plt.plot(x, distribution_1 , linewidth=2, color='r')\n",
    "plt.plot(x, distribution_2 , linewidth=2, color='b')\n",
    "plt.plot(x, distribution_3 , linewidth=2, color='k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In exercise 2 below, we now create 3 lognormal distributions with the same mean but different geometric standard deviations. This is very similar to exercise 1. However, this time you will need to create a linearly spaced array in log(x) which is then used in your expression for defining the lognormal. Please read the instructions carefully.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<b> Exercise 2. Create 3 lognormal distributions and plot the results <a name=\"Exercise2\"> </b> \n",
    "    \n",
    "Using the defition of a log-normal distribution, you are tasked with now creating three and plotting the results. We are going to assume that each distribution has the same mean [in log space] but 3 different geometric standard deviations of 0.3, 1.0 and 1.3. As before, we are going to use the Numpy and Matplotlib libraries to create our arrays and plot results accordingly. I have imported the relevant libraries for you and also provided a plotting function. Notice I have assumed our x variable is now named 'x_log' to go into our definition of the lognormal. Also note that this time we are assuming we are using microns as our x axis. Again I assume our 3 different distributions are called 'distribution_1','distribution_2' and 'distribution_3'.  \n",
    "\n",
    "You should define an x array that has a lower limit of 0.01 microns [10nm] and 20 microns and has 400 values. Your geomtric standard deviations are 0.3, 1.0 and 1.3. Your mean diameter in this problem is 800nm.\n",
    "\n",
    "Once complete you should arrive at a plot that looks like the following:\n",
    "\n",
    "![title <](images/Prac1_ex2.png)\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1fn48c+TDQKCLAl7ZJGAsiQgCCp1xSrg7ldUXFp3FNfWarVotVKtL+tXq1ZUpC71p6VWqQIuICLq1wVkl6CRCAJhMYFEiEkgCfP8/jiZbGSZydzJTMLzfr3mdWfucu6TZPLMmXPPPUdUFWOMMc1fTKQDMMYY4w1L6MYY00JYQjfGmBbCEroxxrQQltCNMaaFiIvUiZOSkrRPnz6ROr0xxjRLy5cv36mqybVti1hC79OnD8uWLYvU6Y0xplkSkU11bbMmF2OMaSEsoRtjTAthCd0YY1oIS+jGGNNCWEI3xpgWwhK6Mca0EJbQjTGmhbCEbkxTWb8eLrsMfvc7yM+PdDSmBYrYjUXGHFTy8uD44+HHH93rFSvgww9BJLJxmRbFaujGNIU//ckl8169oG1b+OgjeP31SEdlWhhL6MaEW1ERvPSSez53Lvz1r+75s89GLCTTMllCNybc5syBPXtg1CgYNgwuvRTatIHFi+GHHyIdnWlBLKEbE27z5rnlxRe7Zfv2MGGCe/7ee5GJybRIltCNCSefDxYscM/Hjatc73/+/vtNH5NpsSyhGxNOX38NubmQkgJHHFG5/vTT3fKjj2D//sjEZlocS+jGhNPnn7vlCSdU76LYqxf07g0FBZCREZnYTItjCd2YcPryS7c85pgDtx17rFt+8UXTxWNaNEvoxoSTP1n7k3dVxx1XfR9jQmQJ3Zhw2bPH3e6fkABpaQduP/pot1y5smnjMi2WJXRjwmXNGrccPBji4w/cPnSoa1dftw727Wva2EyLZAndmHBZvdot09Nr3962LaSmQlkZfPNN08VlWixL6MaES0MJHdydowCrVoU/HtPiWUI3Jlz83RGHDq17nyFD3HLduvDHY1o8S+jGhEtmpltWvaGopoEDq+9rTAgaTOgi8oKI5IjI2jq2Xyoia8ofn4tIPd8vjTlI7NrlHm3bQo8ede83YIBbfvdd08RlWrRAaugvAePq2b4ROFFV04BpwAwP4jKmefMn6AED6p/EIjXVLb//3l0cNSYEDSZ0Vf0EyKtn++eq6p9P60ugl0exGdN8+RO6v0mlLm3bunFeSkth48bwx2VaNK/b0K8GbDxQY/xt4v4mlfpYO7rxiGcJXUROxiX039ezz3UiskxEluXm5np1amOiT6A19Kr7WEI3IfIkoYtIGjATOEdVd9W1n6rOUNWRqjoyOTnZi1MbE52CqaH797GEbkIUckIXkcOA2cDlqmqX6o3x+dwYLmBNLqZJxTW0g4j8CzgJSBKRbOA+IB5AVZ8F/gh0BqaLu5pfpqojwxWwMVFv82Y3Nkv37m66uYZYQjceaTChq+qkBrZfA1zjWUTGNHdZWW7p75LYkMMOc4N3/fgjFBa6ni/GNILdKWqM1374wS379g1s/5gYN3sRwKZNYQnJHBwsoRvjNX9C79Mn8GP8+/qPNaYRLKEb47Vga+hgCd14whK6MV5rTA3dn/wtoZsQWEI3xmuhNLnY7f8mBJbQjfHSvn2wbRvExkLPnoEfZ00uxgOW0I3x0pYtoOoG3IprsFdwJUvoxgOW0I3xUmOaWwC6dYOEBNi5E37+2euozEHCEroxXmpsQre+6MYDltCN8VJjEzpYTxcTMkvoxngplIRu7egmRJbQjfGSPxn7m0+C4T/GErppJEvoxnhp61a37NWImRj9x/jLMCZIltCN8Yqq64MO0KNH8Mf7E3p2tncxmYOKJXRjvLJrF5SUQIcO0KZN8MdbDd2EyBK6MV7x186DuUO0Kv9x2dmutm9MkCyhG+MVf826Mc0t4Ca26NjR1fJ37vQuLnPQsIRujFf8Cb2xNXSwdnQTEkvoxngllAuifpbQTQgsoRvjFauhmwizhG6MV0K9KAqW0E1IGkzoIvKCiOSIyNo6touIPCkiWSKyRkSO8j5MY5qBUC+KgiV0E5JAaugvAePq2T4eSC1/XAc8E3pYxjRD1uRiIqzBhK6qnwB59exyDvBPdb4EOohId68CNKZZKC2FnBw3DG6XLo0vxxK6CYEXbeg9gS1VXmeXrzuAiFwnIstEZFlubq4HpzYmSmzf7pbdugU3U1FNVRO63VxkguRFQpda1tX6TlTVGao6UlVHJicne3BqY6KEF10WAdq3h3btoKgIfvop9LjMQcWLhJ4NpFR53QvY5kG5xjQfXrSf+1mzi2kkLxL6HOBX5b1djgF2q+p2D8o1pvnwosuinyV000gNNvaJyL+Ak4AkEckG7gPiAVT1WeBdYAKQBRQBV4YrWGOilhddFv0soZtGajChq+qkBrYrcKNnERnTHFmTi4kCdqeoMV7w6qIoWEI3jWYJ3RgvWA3dRAFL6MZ4wcsauv9DwWYuMkGyhG5MqAoK3CMx0U0/Fyr/h4IldBMkS+jGhKpql0Wp7T67ICUlQUKCu7GoqCj08sxBwxK6MaHysssiuA8Ff1nb7B49EzhL6MaEyssLon6W0E0jWEI3JlReXhD1swujphEsoRsTqnDW0C2hmyBYQjcmVOGsoVuTiwmCJXRjQmU1dBMlLKEbE6pwJHSroZtGsIRuTCh8vsrZirp7OPOiXRQ1jWAJ3ZhQ5OZCWRl07gytW3tXbtVuizYVnQmQJXRjQhGOC6IAbdvCoYfCvn2QV98c7cZUsoRuTCjC0X7uZxdGTZAsoRsTinDV0MEujJqgWUI3JhRWQzdRxBK6MaEIZ0K3GroJkiV0Y0LRFE0uVkM3AbKEbkwomqLJxWroJkABJXQRGScimSKSJSJ31bL9MBH5SERWisgaEZngfajGRCGroZso0mBCF5FY4GlgPDAImCQig2rsdg/wuqoOBy4GpnsdqDFRZ98+2LkTYmOhSxfvy7eLoiZIgdTQRwFZqrpBVUuAWcA5NfZRoH3580MB+45oWr6qt/zHhKH1sls3N3tRTg6UlnpfvmlxAnkX9gS2VHmdXb6uqvuBy0QkG3gXuLm2gkTkOhFZJiLLcnNzGxGuMVEknO3nAHFx0LWru/V/x47wnMO0KIEk9Npmva05uMQk4CVV7QVMAF4RkQPKVtUZqjpSVUcmJycHH60x0STcCb1q2XZh1AQgkISeDaRUed2LA5tUrgZeB1DVL4DWQJIXARoTtcJ5QdTPLoyaIASS0L8CUkWkr4gk4C56zqmxz2ZgLICIHIlL6NamYlq2pqih24VRE4QGE7qqlgE3AfOBb3C9WTJE5AERObt8t9uBa0VkNfAv4ApVG/PTtHBNWUO3JhcTgLhAdlLVd3EXO6uu+2OV5+uAMd6GZkyUsxq6iTJ2p6gxjeVPslZDN1HCEroxjaFamWSthm6ihCV0Yxpj924oKoJDDoH27Rvev7Gshm6CYAndmMZoiguiAB07urlK9+yBn38O77lMs2cJ3ZjGaIoLouBu/bdmFxMgS+jGNEZT1dDBml1MwCyhG9MYTVVDB6uhm4BZQjemMZqiy6Kf1dBNgCyhG9MYTdFl0c9q6CZAltCNaYymbHKxGroJkCV0YxojEhdFrYZuGmAJ3Zhg7d9fOeFE9+7hP581uZgAWUI3Jlg5OS6pJydDQkL4z+dP6Nu3g88X/vOZZssSujHBCrL9PCcH/vIXOPlk6NsX+vWD886D55+HvXsDKCAx0d0xWlrqJqU2pg6W0I0JVoBdFvfvh4cecgn8D3+AxYvhhx9g40Z46y247jpITYU33gjgnHZh1ATAEroxwQqgy2JeHowfD1OnQmEhnHGGS9zr18O6da52npYG2dkwcSLcdpv7AKiTtaObAAQ0wYUxpooGmlzy8mDsWFi1yjWz/7//B6edVn2fI4+Eq66Cv/8dfvc7eOIJd9yLL0JsbC2FWg3dBMBq6MYEq54ml+JiVzNftQr694flyw9M5n4xMXDLLbBwIbRtC6+8AlOmuKHWD2BdF00ALKEbEyx/Uu3Vq9pqVdcuvnSpu/j50UeQktJwcSecAO+950bJnTEDnn22lp2sycUEwBK6McHKznbLGgl9+nTXvNK2LcyZc8Dmeh1/PMyc6Z7fcgssWVJjB2tyMQGwhG5MsGpJ6JmZri0c4IUXYMiQ4Iu99FK49VYoK4Nf/cpNiFTBaugmAAEldBEZJyKZIpIlInfVsc+FIrJORDJE5DVvwzQmSuzZAwUF0KYNdOgAuN4pv/6161P+q1/BhRc2vviHH4ZBg+C771xXxwpWQzcBaDChi0gs8DQwHhgETBKRQTX2SQXuBsao6mDgtjDEakzkVa2diwCuRr5kicu5TzwRWvGtW7uLo7Gx8NRTsGJF+YYuXdzK3FzYty+0k5gWK5Aa+iggS1U3qGoJMAs4p8Y+1wJPq2o+gKrmeBumMVGiRnNLfn5lTfqxxyoq7SE56ijXju7zuV4vPh8umXfr5nbYvj30k5gWKZCE3hPYUuV1dvm6qgYAA0TkMxH5UkTG1VaQiFwnIstEZFlubm7jIjYmkmok9Pvuc3fjn3iiu0HIK/ff78b9WrIEXn21fKU1u5gGBJLQpZZ1NXvKxgGpwEnAJGCmiBxQV1HVGao6UlVHJicnBxurMZFXJaFnZLieLTEx8OSTFS0wnmjf3o3/AnDvveWtLHZh1DQgkISeDVTtTdsLqFlFyAbeVtVSVd0IZOISvDEtS5WE/sc/uguikye72/i9dtllrrfMpk3lfdOthm4aEEhC/wpIFZG+IpIAXAzMqbHPW8DJACKShGuC2eBloMZEhfKEvnLfIGbPdhcx7703PKeKjXWDewH8+c9Q0LmPe2E1dFOHBhO6qpYBNwHzgW+A11U1Q0QeEJGzy3ebD+wSkXXAR8AdqrorXEEbEzHlCf2+/w4D3EXLcM5xceaZMGaMa6d//OuxbqXV0E0dRGsdOCL8Ro4cqcuWLYvIuY1ptE6dWJrfn9EspU0bNxRuly7hPeXHH8NJJ0HHdqVsKuhEu5NGunEFzEFJRJar6sjattmdosYEqrAQ8vOZJvcBcPPN4U/m4MZ6GTMG8gvieY7J1uRi6mQJ3ZhAbd3KNxzBPD2D1q3ht79tmtOKVPZ1/19uZ+/WXXUMyWgOdpbQjQlUdjb/y+0AXHFF09TO/caPh2HDlB1056WiiW4IAmNqsIRuTIC2r93FK1yO4Guy2rmfq6W7ju6PcCf7t9iFUXMgS+jGBOip2T0ooRXnpmaQGoG7LM4/H/q13spG+jH3zZKmD8BEPUvoxgTg55/hmS+GA3DHGd9EJIbYWLh50CIAnpjVNSIxmOhmCd2YAPzjH/BTSRvG8H8ce3zkpuK9akwm7djD4m+7sWpVxMIwUcoSujEN8PncZM4At/O/wU1F5LH2/ZK4khcBN36MMVVZQjemAR9+CFlZkBKTzdkEObec13r14maeQvDx6quQYwNVmyosoRvTAP+kzdf5niU2vsq45JHQuzf9+Z4z239CSQk891zkQjHRxxK6MfXYuhXefhvi4pSr+QekpLjxciOlTx8Abtv/GOCG7y0tjVw4JrpYQjemHv/4hxsi99xjfqQ7O6B378gGlJQEbdpwcuFcBh2xnx073AeOMWAJ3Zg6lZXBjBnu+fUjyweSi3RCF4E+fRDg+nN2AJVNQsZYQjemDvPmuSaXAQPglMQv3MpIJ/QqMVw+dBWJie6i7XffRTgmExUsoRtTB3/Nd/JkkM2b3ItoSOjl7egddn3PpElulf+bhDm4WUI3phbffw/z57sZia64AjcPHERVQueHH7j+evf0xRehuDhiEZkoYQndmFr4a7wXXQSdOhG1Cf3oo2HECMjLgzfeiGhUJgpYQjemhn374IUX3PPrr8f1C9y61V2QTEmp99gm4f9Q+eEHgIpaul0cNZbQjanhjTfcHJ7DhsHo0bhk7vO5yUMTEiIdXmUNvfxbw8UXQ/v28PnnsGZN5MIykWcJ3Zga/DXd6693lfKoam4BN7NG69aunWXPHg45BC6/3G2yO0cPbgEldBEZJyKZIpIlInfVs98FIqIiUusEpsZEu6+/hv/7P2jXDi65pHxltCX08r7oQEVs/maXV15xQ/2ag1ODCV1EYoGngfHAIGCSiAyqZb92wC3AEq+DNKap+Gu4l13mkjoAGze6pT+JRoMaCX3IEPjFL6CgAF59NXJhmcgKpIY+CshS1Q2qWgLMAs6pZb9pwCPAXg/jM6bJ/Pwz/POf7vkNN1TZsGGDWx5+eJPHVKcaF0ahspb+3HM2h/TBKpCE3hPYUuV1dvm6CiIyHEhR1Xn1FSQi14nIMhFZlpubG3SwxoTTa6+5Gu6YMTB0aJUN33/vlv36RSSuWlXpuuh3wQVuqJeVK2Hp0ohEZSIskIQutayr+PwXkRjgcSifDr0eqjpDVUeq6sjk5OTAozQmzFSrXwytJhpr6LUk9Fat4Mor3XPrwnhwCiShZwNVO9/2AqpOOd4OGAIsFpEfgGOAOXZh1DQnS5e6mm3nzq6mW6GoCLZvh/j4yE5sUVONNnS/665zy1mzID+/aUMykRdIQv8KSBWRviKSAFwMzPFvVNXdqpqkqn1UtQ/wJXC2qi4LS8TGhIG/Rnvlla5HYAV/7bxPHzdLc7Twt6H7L9iW698fTjsN9u6Fl1+OQFwmohpM6KpaBtwEzAe+AV5X1QwReUBEzg53gMaEW16eq9GCG4irmmhsbgE3a1KbNrBrF/z0U7VNVe8ctYujB5eApi9X1XeBd2us+2Md+54UeljGNJ2XX3Y12l/+0tVwqwnyguiePXtYsWIFGRkZZGRkkJWVRU5ODrm5ueTn56PlGTY2NpakpCSSk5Pp1q0bAwcOZPDgwQwZMoT09HQSGrojVcQFu2YNrF8PRx9dsemss6BHD8jMhMWL4eSTA/xFmGYvoIRuTEvl81U2t1TrqujXQA29rKyMTz75hPnz57N48WKWL1/O/v37Azp3YWEhm8rbwOfNq+wglpiYyJgxYzjxxBM566yzSEtLQ6SWvgmpqS6hZ2VVS+hxcXDttfCnP7mfzRL6wcMSujmoLVrkJofo2dPVbA/gr6FXSeiqyuLFi3nttdd466232LlzZ8W22NhYRo4cSVpaGoMHD2bgwIF069aNLl260LFjR2LL2+FLS0vZtWsXOTk5bN26lXXr1pGRkcHKlSvJzMxk4cKFLFy4kHvvvZfDDz+cCy64gMsvv5zBgwdXxub/OrF+/QFhX3MNTJsGs2fDjz9C164h/6pMc6CqEXmMGDFCjYm0885TBdUHHqhjhwED3A5r1mhOTo4+8sgj2r9/f8V13VVAU1NT9c4779T58+drQUFByDHt2LFDX3/9db322ms1OTm52rnGjBmjL7/8shYXF6vOnOliu+yyWss55xy3+cEHQw7JRBFgmdaRVy2hm4PW5s2qMTGqcXGq27bVskNpqWpCgm4Cvfn66zUxMbEisfbs2VPvueceXbNmjfp8vrDFWFZWph999JFOnjxZ27VrV3H+rl276kPXXqv5oDp6dK3Hvv+++w/v3Vu1rCxsIZomZgndmFrcc4/7D7jootq3r1+4UK8AjatSQx4/frzOmTNHS0tLmzZYVS0oKNCZM2fqsGHDKuJpBzo1MVF/+umnA/bfv1+1b1/3M77zTpOHa8LEEroxNezbp9q1q/sP+Pjj6ttycnL0pptu0rjYWAU0BnTSpEm6evXqyARbg8/n0wULFujYsWMrEnunjh310UcfdU0xVTz8sPsZzzwzQsEaz1lCN6aGf/3LvfuHDFH1t5gUFRXpgw8+WNG0ISJ6Bej6Sy+NbLD1+Lx/fz2hyjeIlJQUfe211yqagXJyVOPjVUVUf/ghwsEaT9SX0G2CC3NQmj7dLadMcV2633//fYYMGcLUqVMpKChg/PjxrD7/fF4E+h9zTERjrc+xw4axGHj3jjtIT09ny5YtXHLJJZxyyilkZGSQnOyGMlCtnCfVtFyW0M1B5+uv4dNP3XjnY8du48ILL2T8+PFs2LCBIUOGsHDhQt59912G5uS4AwYOjGzA9UlNRYDxrVuzYsUKZs6cSVJSEosXL2bYsGHcfvvtXHGFm/Fixgx3A5VpuSyhm4PO008D+Bg27ClGjjyC//znP7Rp04ZHHnmEFStWMHbsWLdjZqZbHnFEpEJt2KDyuWbWrSMmJoarr76azMxMpkyZwv79+3nssceYPHko/ft/yM6dbohg04LV1RYT7oe1oZtIyM1VbdXqe4UTKtqdzznnHP2hZgPzrl2ukb1tW9ddJFqtWOHiPPLIAzYtX75chw8fXqUf+3U6aNBPGsZelqYJYG3oxrjKy9VXP8e+fWnAJ3Tt2pX//ve/vPXWW/SuOV+ov3Y+YADERPG/ycCB7iLA+vVQUlJt01FHHcWSJUv485//XD42zAzWrRvCI4+8H5lYTdhF8TvVGO9s2bKF004bx5w51wOFnHTSRWRkZHDuuefWfkBzaG4BN+Jiv35QVlbrEADx8fFMnTqVFStW0LPnKCCbu+4az0033URxcXHTx2vCyhK6afFmzZrF0KFDWbhwAdCJlJRZLFo0i86dO9d90Lp1bhntCR0q29EzMurcZfDgwSxd+hmxsQ8D8Tz99NOMHDmSVatWNU2MpklYQjctVmFhIVdddRWTJk1i9+7dtG9/FpDBn/50EbUNXljNmjVuWW1y0ShV5cJofXr0iOPyy38PfEnHjgNZt24do0eP5tFHH8Xn84U/ThN2ltBNi7Ry5UpGjBjBiy++SOvWrbnttmfZs+dtunbtxiWXBFCAP6Gnp4c1Tk/4R2BsIKED/OY3AEexd+8KrrzyBkpKSrjjjjs47bTT2Lp1a1jDNOFnCd20KKrKE088wTHHHENmZiZDhgxh2bJlZGZOBoQpU9xkyvXKzXXziB5ySOXcndEsgCYXv7Q0mDABiovbkJIynblz55KcnMyHH35IWloab731VpiDNWFVV/eXcD+s26LxWk5Ojp5xxhkV3fRuuOEGLSoq0uXLXc++Nm1ct8UGffihO+DYY8MesycKC929/XFxqnv3Nrj7p5+6H69jR9U9e9xwvePHj6/4vU2ZMkWLioqaIHDTGFi3RdPSLVq0iPT0dN555x06duzI7NmzmT59OomJiTz4oNvnhhsgKSmAwvzNLWlpYYvXU23auO6LZWUB1dJ/8QsYMwby893do127duWdd97hb3/7GwkJCUyfPp1Ro0axdu3aJgjeeMkSumnWysrKmDp1Kqeeeirbt2/n+OOPZ9WqVZx33nmAy2+zZ7tmlttvD7DQ5pbQAYYPd8sVKwLa/e673fKxx2DfPhARbr31Vr788ksGDhzI2rVrOfroo3n22WfdKH6mWbCEbpqtzZs3c+KJJ/LQQw8hItx///0sWrSIww47rGKfhx5yy2uuge7dAyy4OSb0o45yywAT+oQJ7sfbtg3++c/K9cOHD2f58uVcffXV7N27lxtuuIHzzz+fvLy8MARtPFdXW0zVBzAOyASygLtq2f5bYB2wBvgQ6N1QmdaGbkIxe/Zs7dChQ8XsQYsXLz5gn+++q5yRaNOmAAsuLVVt3do1MtcyaUTU8rf71zF7UW1ee80d0revGx++plmzZmn79u0V0F69etX6OzZNj1DGQwdige+BfkACsBoYVGOfk4E25c9vAP7dULmW0E1jFBcX64033lhxAe/MM8/U3DqudF55pXuHX311ECfwj41y+OHeBNxU8vJc3K1buw+lAJSVqR5xhDvsmWdq32fjxo167LHHuok+YmL03nvvjchsTaZSqAn9WGB+ldd3A3fXs/9w4LOGyrWEboL17bffanp6ugIaHx+vf/vb3+qcz/Obb1ztPDZWdf36IE7yzDPu3+KSS7wJuin555v7+uuAD3n9dXdIjx6qdXVsKS0t1XvuuUdFRAE97rjjDhzMzDSZ+hJ6IG3oPYEtVV5nl6+ry9XAe7VtEJHrRGSZiCzLzc0N4NTGOC+//DIjRoxg9erVHH744XzxxRfceuutSB23fE6dCj4fXHst9O8fxImWLnXL0aNDD7qpBdmODvA//+Pundq2DZ57rvZ94uLimDZtGosWLaJnz558/vnnpKen85///MeDoI2n6sr0WlnjngjMrPL6cuCpOva9DPgSaNVQuVZDN4HYs2ePXn755RVNLJdcconu3r273mOWLHG1zsRE1a1bgzzhoEHu4C+/bHzQkfLQQy72KVOCOmzuXHdYcrJqA79a3blzp55zzjkVf49rrrlGf/755xCCNsGiKZpcgFOBb4AuDZWpltBNAJYsWaKpqakKaJs2bfSFF16os4nFz+dTPflk986+664gT7h7t7tBJyEhoBt0os7HH7sfPD09qMN8PtXjjgv8d+bz+fTpp5/WVq1aKaBHHHGErlq1qpFBm2CFmtDjgA1AXyovig6usc9w3IXT1IbKU0vopgGlpaV6//33a2xsrAI6dOhQXbduXUDHvvWWe1d36OCuEwZl4UJ38KhRwQcdDYqKKmeEDrKHzpdfuh+9VSvVjRsDO2bNmjU6ePBgBTQhIUGffPLJBj9wTehCSujueCYA35Un7anl6x4Azi5/vhD4EVhV/pjTUJmW0E1t1q9fr6NHj674Sv/b3/5Wi4uLAzq2uLjyuuCTTzbi5A8+6A6++eZGHBwljj3W/QzvvRf0oZdc4g698MLAjyksLNTrr7++Wq+jnJycoM9tAhdyQg/HwxK6qcrn8+nzzz+vbdu2rej3/OGHHwZVhj8fDx4ccM+96k491RUwa1YjDo4Sd9zhfoY//CHoQzdtquyC//HHwR07e/Zs7dixowLavXt3XbhwYdDnN4GxhG6i2rZt2/Tss8+uqOVNmjRJ84JsL9m82Q2+BaqLFjUiiOLiymzWnGuYb7/tfoYTTmjU4ffd5w4/4ojgLyNs3rxZTzjBzdUqIvq73/3OBvkKA0voJir5fD594YUXKu74PPTQQ/XVV19tRDmqEya4d/PEiY0MZtEiV0BaWiMLiBI7d2pFY3iATVVVFRerDhjginjggeBPX1ZWpqGwwmQAABIvSURBVA888IDGxMQooAMGDNBPPvkk+IJMnSyhm6izceNG/eUvf1lRK58wYYJu3ry5UWW98opWXAjdtq2RAU2d6gq57bZGFhBF0tPdzzJ/fqMO93+2tWqlmpnZuBCWLFlSccEU0BtvvFH37NnTuMJMNfUldBucyzQpn8/H3//+d4YMGcIHH3xAp06deOWVV5g3bx4pKSlBl/fjj3Drre75Y48FMQBXTR9+6JZjxzaygCgyYYJbvlfr/X0NOvlkuOIKNwrjVVe5UXmDNWrUKJYvX869995LXFwcTz/9NEOGDGH+/PmNiskEqK5MH+6H1dAPPkuXLtWjjz66otY2ceJE3bFjR6PL8/lUzz7b1SZ/+Uv3ulHy890YAbGxDd9Z0xz4Z7AYMKDRReza5YYDANVp00ILZ9WqVXrUUUdV/N0vvPBCzc7ODq3QgxjW5GIiaefOnTp58uSKsUB69Oihb775ZsjlPvmkewcfeqhqSEOL+NtsTjop5JiiQmmpa3+CIAeyqW7BAldEbKy7+za0kEr14Ycf1sTERAW0bdu2+sgjj+i+2oZ5NPWyhG4iorS0VJ977jnt3LmzAhoXF6d33HGHJ22py5e7GzpB9T//CbGwc891BT31VMhxRY2LLnI/0xNPhFTMbbdpxeCTQd+oVYtNmzbp+eefX1FbP/LII3XBggWhF3wQsYRumpTP59O5c+fqoEGDKv5xTz75ZM3IyPCk/J07XYIB1euvD7GwgoLK7ootqRnA/63jF78IqZjiYtVhw1xR48a5IXe98N5772n//v0r3h+nn366rl692pvCWzhL6KbJfPXVV3rSSSdV/KP27dtX//3vf3t2S/i+fa5lBFyiCbmbs3/82OYyIXSgCgoqO+Zv2BBSURs3qnbu7Iq6+25vwlNV3bt3r/7lL3+pmERDRPSKK67QLVu2eHeSFsgSugm7lStXVvsq3alTJ3388cd1r4eDXPl8qtde69613bqpevJ/f845rsBHH/WgsCjjv5c/1Kua6iZEio11xc2c6UFsVeTm5uqtt96q8fHxCmjr1q31lltusQundbCEbsLmq6++qnaXZ6tWrfTOO+/U/Px8z8/l7yreurXq0qUeFLhlS+UcdSH0tola772nFb1dPPiGNH26Ky4mRnX2bA/iqyErK0svvPDCivdSQkKCTpkyRTcFPH/gwcESuvGUz+fTBQsW6Pjx4yv++RITE/U3v/mNbmv0nT31mzZNK3pceJZM/vhHDXo0quaktNR9lQFXxfaAf2iAVq0afd9Sg1avXq0TJ06s6BUVHx+vV111la5ZsyY8J2xmLKEbTxQWFupzzz1X7Q7Atm3b6p133qk//vhjWM7p81UmkZgY1X/9y6OCS0pUu3d3BX/0kUeFRiH/J+G4cZ4U5/O5+TPA9TLyoPdpndauXauTJk2qSOyAjh07VufMmaP79+8P34mjnCV0E5K1a9fq7bffrp06dar4x+rRo4c++OCDdU7Q7IWSEjfBsz+Zv/SSh4X72w8GDfKkOSJq7dpVeXHUo14k+/er3nJL5d9lxgxPiq1TZmam3njjjRUjcQLav39/ffjhh8P2jTCaWUI3QcvLy9Pp06dXu7MT0FGjRulrr72mJSUlYT3/jh2qY8e6d2hiohtE0DOFhZW185A7sTcDt97qftZzz/WsyKrfnED1hhtcD6Rwys/P10cffVR79+5d8X6MjY3VM844Q994442D5iYlS+gmIHl5efryyy/rWWedVTG9GOWjIE6ePFmXhHq7YIAWLaps+k1OVv3iC49P4B84fcSIll0799u6VfWQQ9zP3IiJL+rzwguuPR1UjzlG9bvvPC2+VmVlZTp37lw9//zzNS4uruJ92rFjR73iiit03rx5nvauijaW0E2dNmzYoM8884yOGzeu2j+HiOipp56qr776apONaZ2frzp5cmWt74QTGjHJc0PWrKm8xfSDDzwuPIr99a/uZ+7f331D8dDSpaq9elV+m3r88UZOMNIIOTk5+vjjj2taWlq1b5Lt27fXSy+9VGfNmqU7d+5smmCaiCV0UyE/P1/nzZunN998c8UEzP5HTEyMnnLKKTp9+nTdvn17k8VUUqL67LOVtfL4eNcBxfOkUFhYObTsddd5XHiUKylxUzmB65/u8TeTvDzVyy+v/DAeMkT1/fc9PUWD1q1bp9OmTdP09PRq72sR0VGjRum9996rn376abOvvVtCP0j5fD7NysrSl19+WSdPnqxDhgyp1mPAX5M5//zz9fnnn2/yuSALC90FtX79KhPBcceprl0bhpOVlKiecYZWDExSUBCGk0S5tWtV27Z1v4MHHwzLKd5+W7VPn8q/5y9+oTpnjruQ2pTWr1+vjzzyiJ5yyimakJBQ7T3fqlUrHTNmjP7+97/XuXPnhq2HVrjUl9DFbW96I0eO1GXLlkXk3C3R3r17ycjIYM2aNaxevZrVq1ezZs0a8vLyqu2XkJDAiBEjOPXUUzn99NMZPXo0cXFxTRanzwdLlsC//w3//Cfk57v1AwfCtGnwP/8DMV6P0r97N1x2GcybB507w6efwpFHenySZuKNN2DiRPf87rvhz3/2/Be+dy889RQ89BD89JNb16+f+xNcdhmkpnp6ugYVFhayePFi3n//fRYvXszatWsP2Kd79+4MHz6cYcOGMWzYMNLS0ujXrx/x8fFNG2wARGS5qo6sdZsl9OZBVcnPzyc7O5sNGzaQlZVFVlYW69evJysriy1btlDb3zIpKYnjjjuOMWPGMGbMGEaMGEHr1q2bNPadO+GTT+Cjj+Dtt2HLlspto0fDLbfAhRdCWD5XPvgAbrwR1q+HDh1gwQI4+ugwnKgZeeUVuPJK2L8fTjwRZsyAAQM8P01BAcycCY8/Xv1vPnw4nHaaexx3HDTx25G8vDy++OILPvvsMz777DNWrlxJQUHBAfvFxcXRr18/BgwYwIABA0hNTaV3796kpKRw2GGH0b59+6YNvFzICV1ExgFPALHATFV9uMb2VsA/gRHALuAiVf2hvjIP9oSuqhQXF5OXl3fAY9euXezYsYNt27axbds2tm7dyrZt29i3b1+d5cXGxpKamkp6ejppaWmkp6eTnp5Oz549EZEm+Zl+/hk2b3a58+uv3WPNGvj22+r79erlEvikSTCy1rdliHJzYf58eP5590kCMHQo/Pe/cPjhYThhM/T++25aoh9/dDX0iRPh17+GU06BVq08PdX+/e7D/JVX4M03obCwcltcnPvTjBgBaWnQv7/7E/XpAwkJnoZRJ5/Px4YNG1i1ahUrV65k5cqVrFu3js2bN9daSfJr3749KSkppKSk0KVLF5KTkyseSUlJFcsOHTpw6KGHkuDRDxRSQheRWOA74JdANvAVMElV11XZZwqQpqrXi8jFwHmqelF95TY2oWd/m82ubXn4fD58qvh8vvL2Ix8+n+LzlT/f79b7KtaX7+c/zr/dV30/9flQXDn7fT7KSkspKyulpLSU0rJSykrd87KyUkrLyigtK6G0tJSysjJKy0opKS2heG8RRcXF5csiivYWUVxcRNFet654bzGFRYWUlNadoGvTrm17kjt3oVfXFA7r0ZvDuvfhsO69OaxHH7p36UV8XOXXQ/VV/l0r/sSqKFWSuypV//yqUFICe0uEfSUx7N0n7CuRiuW+EqGgMIa83bHk74khf08sebtjyc2LY8uOeH4qiK017lYJPo5LL+KkkT9z6jE/c0xaUeW3/NrefzXX1bZPaan7Pv/TT67dZvv26p8kFb+0dnDPPe5rQFNXBaPdzp2u2eWllyrnmYuPh2HDID0dUlKgZ0/o0gXatKl8JCZCbKz7IBCpfFR9XfV5FcXF8OmSBBYsTuCDj+P5+ps4VA+scMTEKF2TfXRJ8tElSemS5COpk492hyht2yiHtK2+jI9X4mJd+HGxSlyc+7Co+jxG3PvIH1LV8Gpbt3dfMT9s2cjGzd+zYdP3bNy8gW3bt7Ltx61s27GN4r3FQf26WyW0ot0h7Wjfrj3Dhg7jzblvBnW8X6gJ/VjgflU9vfz13QCq+pcq+8wv3+cLEYkDdgDJWk/hjU3oQzqfS0be20EfF51aAZ2BTjUeHYFuQA+gZ/myO3BIZMIMUGuKSWELfdnIUL6ueAxiHa0J7sMrtEBau+/yF1wAl14KEfpq3Gxs2uQuaLzxhvswbMJm2AIOYRXDWM4IvuFIvudwsujPFlLwUXsFITookAdswdVzc+t47AR2lz/2Vxx9SMxwCvavaNSZ60vogbRa9iyP2i8bGF3XPqpaJiK7cZlqZ41ArgOuAzjssMMCCr6mdq3aEkM3QICYKkuQaq8FEOSA/Wpbf+BriClfGw/EIcTjPqvi3fPydVQ8j0NIAOKJIZEY2iC0IUYSiSEREbcuhkSERGI4hFip/tW2aj1FRKusV+DH8of/df3HHLhf7cfUPK6VlJQ/Smkt+9zzmBJal68/JKaIjrF76BS7u2LZOfYnUuJ3kBT3U80KWXl1J7XG69r2aWBdzdexsa5NvEMH6NjR1SL793dXV486ymrjwejdG+691z327IGvvoLMTNi6FbKzXU2+uBiKiiqX+/e7xO/z+Tu0HPja52vw1O2A48nkeDKrrS/ReHJ8Se6hSeT4ktnp68TP2pZCbUOhtql8ThtKNZ4yYinTOEqJo4w4yjSOMmIpJZ4yjUOr/Be451LxvL51Fdu05jpwqa9ngz+nq9sW46MAZQ+D229r8JjGCCSh19YAW/MjPJB9UNUZwAxwNfQAzn2AL7a92pjDjDGBaN8exo51jwhKAHqVP0zgAumvlA2kVHndC6j58VKxT3mTy6G47yPGGGOaSCAJ/SsgVUT6ikgCcDEwp8Y+c4Bflz+/AFhUX/u5McYY7zXY5FLeJn4TMB/XbfEFVc0QkQdwdyzNAf4BvCIiWbia+cXhDNoYY8yBArqVQ1XfBd6tse6PVZ7vBSZ6G5oxxphgeH2TtTHGmAixhG6MMS2EJXRjjGkhLKEbY0wLEbHRFkUkF9gUkZMHJokad7pGGYsvNBZfaCy+0IQSX29VTa5tQ8QSerQTkWV1jZcQDSy+0Fh8obH4QhOu+KzJxRhjWghL6MYY00JYQq/bjEgH0ACLLzQWX2gsvtCEJT5rQzfGmBbCaujGGNNCWEI3xpgWwhJ6A0TkdyKiIpIU6ViqEpFpIrJGRFaJyAIR6RHpmKoSkb+KyLflMf5XRDpEOqaaRGSiiGSIiE9EoqKLm4iME5FMEckSkbsiHU9NIvKCiOSIyNpIx1KTiKSIyEci8k353/XWSMdUk4i0FpGlIrK6PMY/eVm+JfR6iEgKbnLszZGOpRZ/VdU0VR0GzAP+2NABTewDYIiqpuEmGb87wvHUZi1wPvBJpAOBignZnwbGA4OASSIyKLJRHeAlYFykg6hDGXC7qh4JHAPcGIW/v33AKaqaDgwDxonIMV4Vbgm9fo8Dd1LLdHqRpqp7qrxsS5TFqKoLVLV8Knm+JApnE1PVb1Q1s+E9m8woIEtVN6hqCTALOCfCMVWjqp8QpbORqep2VV1R/rwA+IZAJvxsQur8XP4yvvzh2f+uJfQ6iMjZwFZVXR3pWOoiIg+KyBbgUqKvhl7VVcB7kQ6iGahtQvaoSkjNhYj0AYYDSyIbyYFEJFZEVgE5wAeq6lmMAU1w0VKJyEKgWy2bpgJ/AE5r2oiqqy8+VX1bVacCU0XkbuAm4L5oiq98n6m4r8IRmd07kBijSECTrZv6icghwJvAbTW+yUYFVd0PDCu/rvRfERmiqp5ckzioE7qqnlrbehEZCvQFVosIuOaCFSIySlV3RDq+WrwGvEMTJ/SG4hORXwNnAmMjNcdsEL/DaBDIhOymHiISj0vmr6rq7EjHUx9V/UlEFuOuSXiS0K3JpRaq+rWqdlHVPqraB/ePdlRTJvOGiEhqlZdnA99GKpbaiMg44PfA2apaFOl4molAJmQ3dRBX+/oH8I2qPhbpeGojIsn+Hl8ikgiciof/u5bQm6+HRWStiKzBNQ1FWxetvwPtgA/Ku1Y+G+mAahKR80QkGzgWeEdE5kcynvKLyP4J2b8BXlfVjEjGVJOI/Av4AhgoItkicnWkY6piDHA5cEr5e26ViEyIdFA1dAc+Kv+//QrXhj7Pq8Lt1n9jjGkhrIZujDEthCV0Y4xpISyhG2NMC2EJ3RhjWghL6MYY00JYQjfGmBbCEroxxrQQ/x/63wUmDoKWFQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot 3 different log-normal distributions with a mean diameter of 800nm and geometric standard deviation of 1.1, 1.3 and 2.0\n",
    "# This time assume our smallest starying point is 10nm and largest is 20 microns\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Define lower and upper limit of the size distribution [microns]\n",
    "lower_size =  \n",
    "upper_size =  \n",
    "\n",
    "#------'INSERT CODE HERE'------\n",
    "x_log = np.linspace(np.log(lower_size), np.log(upper_size), num=400) # values for x-axis\n",
    "sigmag1 =  \n",
    "mean1 =  \n",
    "\n",
    "sigmag2 =  \n",
    "mean2 =  \n",
    "\n",
    "sigmag3 =  \n",
    "mean3 =  \n",
    "\n",
    "distribution_1 =  \n",
    "distribution_2 =  \n",
    "distribution_3 =  \n",
    "#------------------------------\n",
    "\n",
    "plt.plot(x_log, distribution_1 , linewidth=2, color='r')\n",
    "plt.plot(x_log, distribution_2 , linewidth=2, color='b')\n",
    "plt.plot(x_log, distribution_3 , linewidth=2, color='k')\n",
    "#plt.savefig('sample_curve.png', dpi=800)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objective 2 -  Populate a discrete 'binned' size distribution from a continuous function <a name=\"Part2\">\n",
    "\n",
    "When modelling aerosol particles we represent the variation in particle size from the following:\n",
    "\n",
    " - Track each particle\n",
    " - Place particles into distinct size ranges\n",
    " - Assume a continuous distribution\n",
    " \n",
    "The first, whilst [possible](https://www.atmos-chem-phys.net/14/5327/2014/), would require too much computational resource for complex systems. The third, whilst still used in large scale schemes and often called the modal approach, has limitations when we come to tracking invidual components. The second is the most popular and is known as the setional approach. This approach requires us to split the continuous size distribution into discrete size bins. We have used the log-normal approach above and will continue to do so. Limiting the number of bins in a model is often the easiest method of reducing computer-time and memory limits. We will discuss this in class. For now we need to cover all particle sizes with a fixed number of bins; The bins should be spread geometrically and we can thus look at a geomteric distribution known as the volume-ratio distribution. In this distribution, a fixed volume ratio between bins defines the centre of the bin and bin width. We will see in this exercise that we can also simply choose linearly seperated values in log diameter space. The figure below, taken from the [Fundamentals of Atmospheric Modelling]() illustrates the concept:\n",
    "\n",
    "\n",
    "<img src=\"images/jacobson_fig131.png\" alt=\"Drawing\" style=\"width: 300px;\"/>\n",
    "\n",
    "\n",
    "There are a number of formaulas we need to use to create this discrete distribution.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "V_{rat}=\\left ( \\frac{d_{Nb}}{d_{1}} \\right )^{3/(N_{b}-1)}  \\tag{1}\n",
    "\\\\\n",
    "D_{rat} = V_{rat}^{1/3}\n",
    "\\\\\n",
    "v_{i}=V_{rat}v_{i-1} \\text{ for i=1 to} N_{b}   \\tag{1}\n",
    "\\\\\n",
    "v_{i}= \\frac{1}{2}\\left ( v_{i,hi} + v_{i,lo} \\right )   \\tag{1}\n",
    "\\\\\n",
    "v_{i,hi}= V_{rat}v_{i,lo}   \\tag{1}\n",
    "\\\\\n",
    "v_{i,lo}= \\frac{2v_{i}}{1+V_{rat}}   \\tag{1}\n",
    "\\\\\n",
    "\\Delta d_{i} = d_{i,hi}-d_{i,lo}=d_{i}2^{1/3}\\frac{V_{rat}^{1/3}-1}{(1+V_{rat})^{1/3}}    \\tag{1}\n",
    "\\\\\n",
    "n_{i} = \\frac{N_{L}\\Delta d_{i}}{d_{i}\\sqrt{2\\pi }ln\\sigma _{g}}exp\\left [ -\\frac{ (log(x) - log(\\mu))^2 } {2 log(\\sigma_{g})^2} \\right ]   \\tag{1}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where each variable is defined as follows:\n",
    "\n",
    "| Variable | definition |\n",
    "|------|------|\n",
    "| $V_{rat}$ | Volume ratio between bins | \n",
    "| $D_{rat}$ | Diameter ratio between bins | \n",
    "| $d_{Nb}$ |  Diameter of largest bin | \n",
    "| $d_{1}$ | Diameter of smallest bin  | \n",
    "| $d_{i}$ | Diameter of particle in bin i  | \n",
    "| $N_{b}$ | Number of size bins | \n",
    "| $v_{i}$ | Volume of a particle in size bin i | \n",
    "| $v_{i,hi}$ | Upper volumne of a particle in size bin i | \n",
    "| $v_{i,lo}$ | Lower volumne of a particle in size bin i | \n",
    "| $\\Delta d_{i}$ | Diameter width of size bin i | \n",
    "| $\\sigma_{g}$ | Standard deviation of lognormal distribution |\n",
    "| $\\mu$ | Mean diameter of lognormal distribution |\n",
    "| $x$ | Diameter array of lognormal distribution |\n",
    "| $N_{L}$ | Total number of aerosol particles |\n",
    "\n",
    "The following schematic illustrates how each variable maps onto a hypothetical log-normal discrete distribution:\n",
    "\n",
    "<img src=\"images/Binned_size_distribution.png\" alt=\"Drawing\" style=\"width: 600px;\"/>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<b> Exercise 3. Create a discrete \"volume ratio\" log-normal distribution and populate number concentrations to each bin <a name=\"Exercise3\"> </b> \n",
    "    \n",
    "This is quite an involved exercise. We need to go through a series of steps to create our discrete lognormal distribution that has a mean of 150nm, standard deviation of 1.7 and total number of particles of 600 per cubic centimetre:\n",
    "\n",
    " - Define a lower and upper limit to our size distribution\n",
    " - Calculate a volume ratio for a specific number of bins\n",
    " - Calculate a diameter ratio\n",
    " - Create an array of diameters\n",
    " - Calculate the diameter width of each size bin\n",
    " - Calculate the log-normal probability density for our diameter array\n",
    " - Combine the last two elements with the total number to populate a discrete number distribution\n",
    " - Plot the results\n",
    " \n",
    "In the code snippet below I have started the exercise with the lower and upper bin boundary and desired number of size bins [30]. I have also defined the name of each required variable and array before leaving the plotting routine. Please note we are still working in microns and the number concentration is given in number per cubic centimetre.\n",
    "\n",
    "Once completed, your plot should look like the following:\n",
    "\n",
    "![title <](images/Prac1_ex3.png)\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.6102620275609394\n",
      "1.1721022975334803\n",
      "[0.01       0.01172102 0.01373824 0.01610262 0.01887392 0.02212216\n",
      " 0.02592944 0.03039195 0.03562248 0.04175319 0.04893901 0.05736153\n",
      " 0.06723358 0.07880463 0.09236709 0.10826367 0.1268961  0.14873521\n",
      " 0.17433288 0.20433597 0.23950266 0.28072162 0.32903446 0.38566204\n",
      " 0.45203537 0.52983169 0.62101694 0.72789538 0.85316785 1.        ]\n",
      "[0.01       0.01172102 0.01373824 0.01610262 0.01887392 0.02212216\n",
      " 0.02592944 0.03039195 0.03562248 0.04175319 0.04893901 0.05736153\n",
      " 0.06723358 0.07880463 0.09236709 0.10826367 0.1268961  0.14873521\n",
      " 0.17433288 0.20433597 0.23950266 0.28072162 0.32903446 0.38566204\n",
      " 0.45203537 0.52983169 0.62101694 0.72789538 0.85316785 1.        ]\n",
      "[0.00915058 0.01072542 0.01257128 0.01473483 0.01727073 0.02024306\n",
      " 0.02372694 0.0278104  0.03259663 0.03820659 0.04478203 0.05248912\n",
      " 0.06152262 0.07211081 0.08452124 0.09906754 0.11611729 0.13610135\n",
      " 0.1595247  0.18697927 0.21915883 0.25687657 0.30108562 0.35290314\n",
      " 0.41363859 0.48482674 0.56826653 0.66606651 0.78069808 0.91505802]\n",
      "[0.00157484 0.00184587 0.00216355 0.0025359  0.00297233 0.00348388\n",
      " 0.00408346 0.00478623 0.00560996 0.00657544 0.00770709 0.0090335\n",
      " 0.01058818 0.01241044 0.0145463  0.01704975 0.01998405 0.02342335\n",
      " 0.02745457 0.03217956 0.03771774 0.04420905 0.05181753 0.06073544\n",
      " 0.07118815 0.0834398  0.09779998 0.11463158 0.13435993 0.15748359]\n",
      "0.9916248890828715\n",
      "594.974933449723\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAZX0lEQVR4nO3debxddX3u8c/DJAooU6BM4TAE0FoIcJjqhCJeLaNUplJvoNRAW1tbsRqlWriWe6NWKcqgsaCpMqmFV5C0FRpApCoYKCgQJjFMCRAQFFDBwHP/WL9Tdk72yd4nOWvvnLOe9+u1X3uv+bv2SZ699m+t9duyTURENMca/S4gIiJ6K8EfEdEwCf6IiIZJ8EdENEyCPyKiYRL8ERENk+CPlSbpi5I+3u86xkIv9kXSaZK+Xl5PlvSspDXHaN3/U7+k/SU9PBbrLet7k6S7x2p90X/KdfzRjqSFwObAUuBF4E7gX4BZtl/qY10Gpti+r181rCxJpwE72v7jUSxzPPCntt84imX2B75ue+vR1liWH7fvcXQnR/yxIofY3gDYFpgJfAQ4v66NSVqrrnVPZGP1rSGaI8EfHdn+he0rgKOBaZJeDyDpq5L+obzeVNKVkp6W9HNJ35O0Rpm2jaTLJC2R9KSks8v44yX9l6QzJf0cOK2M/xNJCyQ9Jek7krYt468vJd1WmkmOLuMPlnRr2fb3Je3abj9UOVPS45J+IenHI+zLt8v6hx4vlSNvJO0i6eqyj3dLOmqk903SdpK+K+kZSVcDm7ZMG5DkoQ+78l7cX+b9maTjJL0W+CKwX6nj6ZZaz5P0b5KeA97aWn/LNj4m6QlJCyUd1zL+Okl/2jJ8vKQbRnqPhzcdSXptWcfTku6QdGjLtK9KOkfS3LIvN0raYaT3KPojwR9ds30T8DDwpjaTTynTJlE1EX0McDkavRJ4ABgAtgIuaVluH+B+YDPgDEmHl2WPKOv6HnBx2f6byzK72V7f9qWS9gAuAE4CNgG+BFwh6RVtanwH8GZgJ2BDqg+yJ9vs5yFl/esD7wEeBeZJWg+4Grio1HsscK6k3x3hLbsIuJkq8D8JTGs3U1nv54F3lW9Yvw/cansBcDLwg1LPhi2L/RFwBrABcEOb1f5O2e5WZbuzJO08Qp2t+77cezys1rWBbwNXUb0HfwlcOGzdxwKnAxsB95U6YzWS4I/RWgRs3Gb8b4EtgG1t/9b291ydQNob2BL4W9vP2f6N7dagWmT7C7aX2v41VYD/P9sLbC8F/i8wdeiov433AV+yfaPtF23PBp4H9h2hxg2AXajOby2wvXikHZW0E9V5jaNtPwQcDCy0/ZVS7y3Av1J9OAxfdjKwF/Bx28/bvp4qMEfyEvB6Sa+0vdj2HSuYF2CO7f+y/ZLt34wwz9C2vwvMBUb8djIK+wLrAzNtv2D7GqoP9mNb5rnM9k3l73chMHUMthtjKMEfo7UV8PM24z9DdXR3VWmymFHGbwM8UEKgnYeGDW8LnFWaEZ4u21LZbjvbAqcMzV+W2Ybqw2YZJaTOBs4BHpM0S9Kr261U0muAOVTh+b2Wbe0zbFvHUR1dD7cl8JTt51rGPdBuW2Weo6mO7heXZpJdRtjfIcPft+HabXu592QlbAk8NOwE/wMs+/d5tOX1r6g+KGI1kuCPrknai+o/+HJNC7afsX2K7e2BQ4APSjqAKqAma+QTt8MvK3sIOMn2hi2PV9r+/gjLPwScMWz+V9m+uO3G7M/b3hP4Xaomn79ts59rUDXTXGv7S8O29d1h21rf9p+12dRiYKPSjDNk8gj7gO3v2D6Q6lvTXcCXhyaNtMhI6yrabXtRef0c8KqWae0+uEayCNhm6PxNy7ofGcU6os8S/NGRpFdLOpiqbf7rtn/SZp6DJe0oScAvqS4BfRG4iSoEZ0paT9K6kt6wgs19EfjoULu5pNdIOrJl+mPA9i3DXwZOlrRPOXm7nqSDJG3Qpsa9ynxrU4Xfb0qNw50BrAd8YNj4K4GdJL1X0trlsVc5CbsM2w8A84HTJa0j6Y1UH4jLkbS5pENLUD8PPNtS12PA1pLWabdsB0PbfhNVM9U3y/hbgSMkvUrSjsCJw5Yb/h63upHqvftw2f/9y35dMsL8sRpK8MeKfFvSM1RHuqcCnwNOGGHeKcB/UoXWD4BzbV9n+0WqYNgReJDqBPDRI23Q9uXAp4BLJP0SuB14V8sspwGzS1PLUbbnU7Xznw08RdXcdPwIq3811QfFU1TNE08C/9hmvmOp2rKf0stX9hxn+xmqE8THUB35PlpqbXciGaoTsPtQNVf9PdX5gnbWoDo5vqjM+xbgz8u0a4A7gEclPTHC8u08WvZzEVU7+8m27yrTzgReoAr42WV6q9NoeY9bJ9h+ATiU6m/yBHAu8L9b1h3jQG7giohomBzxR0Q0TII/IqJhEvwREQ2T4I+IaJhx0SnWpptu6oGBgX6XERExrtx8881P2J40fPy4CP6BgQHmz5/f7zIiIsYVSW3vFk9TT0REwyT4IyIaJsEfEdEwCf6IiIZJ8EdENEyCPyKiYRL8ERENk+CPiGiYBH9ERMPUdueupJ2BS1tGbQ98gurHKC4FBoCFwFG2n6qrjog6DcyY2+8SRm3hzIP6XUL0WW1H/Lbvtj3V9lRgT6ofXb4cmAHMsz0FmFeGIyKiR3rV1HMA8NPyO6SHUf3cG+X58B7VEBER9C74jwEuLq83t70YoDxv1m4BSdMlzZc0f8mSJT0qMyJi4qs9+CWtQ/XjzN8czXK2Z9ketD04adJyvYpGRMRK6sUR/7uAW2w/VoYfk7QFQHl+vAc1RERE0YvgP5aXm3kArgCmldfTgDk9qCEiIopag1/Sq4ADgctaRs8EDpR0b5k2s84aIiJiWbX+ApftXwGbDBv3JNVVPhHRB93ce5Br/Se23LkbEdEw4+I3dyP6YTzelRvRjRzxR0Q0TII/IqJhEvwREQ2T4I+IaJgEf0REwyT4IyIaJsEfEdEwCf6IiIZJ8EdENEyCPyKiYRL8ERENk+CPiGiYBH9ERMMk+CMiGibBHxHRMAn+iIiGSfBHRDRMgj8iomFqDX5JG0r6lqS7JC2QtJ+kjSVdLene8rxRnTVERMSy6j7iPwv4D9u7ALsBC4AZwDzbU4B5ZTgiInqktuCX9GrgzcD5ALZfsP00cBgwu8w2Gzi8rhoiImJ5a9W47u2BJcBXJO0G3Ax8ANjc9mIA24slbdZuYUnTgekAkydPrrHMiBhuYMbcjvMsnHlQDyqJOtTZ1LMWsAdwnu3dgecYRbOO7Vm2B20PTpo0qa4aIyIap87gfxh42PaNZfhbVB8Ej0naAqA8P15jDRERMUxtTT22H5X0kKSdbd8NHADcWR7TgJnleU5dNUSMpJumjIiJqs42foC/BC6UtA5wP3AC1beMb0g6EXgQOLLmGiIiokWtwW/7VmCwzaQD6txuRESMLHfuRkQ0TII/IqJhEvwREQ2T4I+IaJgEf0REwyT4IyIaJsEfEdEwCf6IiIZJ8EdENEyCPyKiYRL8ERENk+CPiGiYBH9ERMMk+CMiGibBHxHRMAn+iIiGSfBHRDRMgj8iomES/BERDZPgj4homFp/bF3SQuAZ4EVgqe1BSRsDlwIDwELgKNtP1VlHRES8rBdH/G+1PdX2YBmeAcyzPQWYV4YjIqJH+tHUcxgwu7yeDRzehxoiIhqr7uA3cJWkmyVNL+M2t70YoDxv1m5BSdMlzZc0f8mSJTWXGRHRHLW28QNvsL1I0mbA1ZLu6nZB27OAWQCDg4Ouq8CIiKap9Yjf9qLy/DhwObA38JikLQDK8+N11hAREcuqLfglrSdpg6HXwDuA24ErgGlltmnAnLpqiIiI5dXZ1LM5cLmkoe1cZPs/JP0I+IakE4EHgSNrrCEiIoapLfht3w/s1mb8k8ABdW03InpjYMbcjvMsnHlQDyqJ0cqduxERDZPgj4homAR/RETD1H0df0TPddP2HNFkOeKPiGiYBH9ERMMk+CMiGibBHxHRMF2f3JV0BPBGqh43b7B9eW1VRUREbbo64pd0LnAy8BOq/nZOknROnYVFREQ9uj3ifwvwetsGkDSb6kMgIiLGmW7b+O8GJrcMbwP8eOzLiYiIunV7xL8JsEDSTWV4L+AHkq4AsH1oHcVFRMTY6zb4P1FrFRER0TPdBv9twJTy+h7bv6ipnoiIqNkKg1/SOlS/e3sY8DOqcwLbSrocONn2C/WXGBERY6nTyd2/A9YGJtvew/ZUqpO8awEfr7u4iIgYe52C/wjgfbafGRpRXv858O46C4uIiHp0Cv6XbP9q+Ejbz1LdwRsREeNMp5O7lrQRoDbTXqqhnoiIqFmn4H8NcDPtg7+rI35JawLzgUdsHyxpO+ASYGPgFuC9OUkcEdE7K2zqsT1ge3vb27V5bN/lNj4ALGgZ/hRwpu0pwFPAiStXekRErIxuO2nbo81jB0mdLgfdGjgI+OcyLOBtwLfKLLOBw1e+/IiIGK1ub+A6F9iDqn8eAb9HdVPXJpJOtn3VCMv9E/BhYIMyvAnwtO2lZfhhYKt2C0qaDkwHmDx5crtZIiJiJXTbSdtCYHfbg7b3BKZSdc/8duDT7RaQdDDwuO2bW0e3mbXtuQLbs8r2BidNmtRlmRER0Um3R/y72L5jaMD2nZJ2t31/1XrT1huAQyX9AbAu8GqqbwAbSlqrHPVvDSxa+fIjImK0uu6WWdJ5kt5SHucC90h6BfDbdgvY/qjtrW0PAMcA19g+DrgWeE+ZbRowZ9V2ISIiRqPb4D8euA/4a+BvgPvLuN8Cbx3lNj8CfFDSfVRt/uePcvmIiFgFXTX12P61pC8AV1G1yd9te+hI/9kulr8OuK68vh/Ye2WKjYiIVddV8Evan+rSy4VUJ2i3kTTN9vX1lRYREXXo9uTuZ4F32L4bQNJOwMXAnnUVFhER9ei2jX/todAHsH0PVXfNERExznR7xD9f0vnA18rwcVR9+ERExDjTbfD/GfAXwF9RtfFfD5xTV1EREVGfbq/qeR74XHkAIOlS4Oia6oqIiJp028bfzn5jVkVERPTMqgR/RESMQ526Vd5jpEnkqp6I6GBgxtyO8yyceVAPKolWndr4P7uCaXeNZSEREdEbKwx+26PthyciIlZznZp6jljRdNuXjW05ERFRt05NPYeU582A3weuKcNvpep0LcEfPdVNm3FErFinpp4TACRdCbzO9uIyvAW5gSsiYlzq9nLOgaHQLx4DdqqhnoiIqFm3XTZcJ+k7VD1ymuoXta6traqIiKhNt102vL+c6H1TGTXL9uX1lRUREXXp9oh/6AqenMyNiBjnOl3O+TOqpp12bHuHsS8pIiLq1OmIf3DY8BrAUcCHgP+upaKIiKjVCq/qsf2k7SeBp4CDqU7o7gccZPsPV7SspHUl3STpNkl3SDq9jN9O0o2S7pV0qaR1xmhfIiKiCysMfklrSzoJuJPqxO5htv/Y9p1drPt54G22dwOmAu+UtC/wKeBM21OoPlBOXKU9iIiIUenU1PMzYCnwT8CDwG6SdhuauKIuG2wbeLYMrl0eBt4G/FEZPxs4DThvJWqPiIiV0Cn4ry7Pu5bHEFGF+Aqv8pG0JtVv8+5IdafvT4GnbS8tszwMbDXCstOB6QCTJ0/uUGZERHSrU/DfThXwanmGka/0WYbtF4GpkjYELgde2262EZadBcwCGBwc7Gp7ERHRWafgX7887wzsBcyhCv9DqH5wvSu2n5Z0HbAvsKGktcpR/9bAotEWHRERK6/TVT2n2z4d2BTYw/aHbJ8C7EkV2iOSNKkc6SPplcDbgQVUVwa9p8w2jerDJCIieqTbO3cnAy+0DL8ADHRYZgtgdmnnXwP4hu0rJd0JXCLpH6juBTh/dCVHRMSq6Db4vwbcJOlyqjb5d1NdkTMi2z8Gdm8z/n5g71HWGRERY6TbTtrOkPTvvNxJ2wm2c+duRMQ4NJpO2m4BbqmxloiI6IFuf4glIiImiAR/RETDJPgjIhomwR8R0TAJ/oiIhknwR0Q0TII/IqJhEvwREQ2T4I+IaJgEf0REwyT4IyIapuu+eiIi6jAwY27HeRbOPKgHlTRHjvgjIhomwR8R0TBp6onVRjdf+SNi1eWIPyKiYRL8ERENk+CPiGiY2oJf0jaSrpW0QNIdkj5Qxm8s6WpJ95bnjeqqISIillfnEf9S4BTbrwX2Bf5C0uuAGcA821OAeWU4IiJ6pLbgt724/EA7tp8BFgBbAYcBs8tss4HD66ohIiKW15M2fkkDwO7AjcDmthdD9eEAbDbCMtMlzZc0f8mSJb0oMyKiEWoPfknrA/8K/LXtX3a7nO1ZtgdtD06aNKm+AiMiGqbW4Je0NlXoX2j7sjL6MUlblOlbAI/XWUNERCyrzqt6BJwPLLD9uZZJVwDTyutpwJy6aoiIiOXV2WXDG4D3Aj+RdGsZ9zFgJvANSScCDwJH1lhDREQMU1vw274B0AiTD6hruxERsWK5czciomES/BERDZPgj4homAR/RETDJPgjIhomwR8R0TAJ/oiIhknwR0Q0TII/IqJhEvwREQ1TZ189ERFjYmDG3I7zLJx5UA8qmRhyxB8R0TAJ/oiIhklTT/REN1/VI6I3csQfEdEwCf6IiIZJ8EdENEyCPyKiYRL8ERENk+CPiGiY2oJf0gWSHpd0e8u4jSVdLene8rxRXduPiIj26jzi/yrwzmHjZgDzbE8B5pXhiIjoodqC3/b1wM+HjT4MmF1ezwYOr2v7ERHRXq/b+De3vRigPG/W4+1HRDTeantyV9J0SfMlzV+yZEm/y4mImDB6HfyPSdoCoDw/PtKMtmfZHrQ9OGnSpJ4VGBEx0fU6+K8AppXX04A5Pd5+RETj1Xk558XAD4CdJT0s6URgJnCgpHuBA8twRET0UG3dMts+doRJB9S1zYiI6Cz98UfEhJCfZ+xegj9WWX5kJWJ8WW0v54yIiHok+CMiGibBHxHRMAn+iIiGSfBHRDRMgj8iomES/BERDZPgj4homAR/RETDJPgjIhomXTZER+mSISaK9OdTyRF/RETDJPgjIhomwR8R0TAJ/oiIhknwR0Q0TK7qiYho0YQrfxL8DZdLNSOaJ009EREN05cjfknvBM4C1gT+2fbMftQREbEyxntzUM+DX9KawDnAgcDDwI8kXWH7zl7XMtGlGSci2unHEf/ewH227weQdAlwGJDgLxLYEePfWPw/rutbQz+CfyvgoZbhh4F9hs8kaTowvQw+K+nuldzepsATK7nseDHR9zH7N/5N9H2sZf/0qVVexbbtRvYj+NVmnJcbYc8CZq3yxqT5tgdXdT2rs4m+j9m/8W+i7+N4279+XNXzMLBNy/DWwKI+1BER0Uj9CP4fAVMkbSdpHeAY4Io+1BER0Ug9b+qxvVTS+4HvUF3OeYHtO2rc5Co3F40DE30fs3/j30Tfx3G1f7KXa16PiIgJLHfuRkQ0TII/IqJhGhH8kk6T9IikW8vjD/pdUx0kfUiSJW3a71rGmqRPSvpx+ftdJWnLftc0liR9RtJdZR8vl7Rhv2saS5KOlHSHpJckjZvLHjuR9E5Jd0u6T9KMftfTrUYEf3Gm7anl8W/9LmasSdqGqhuMB/tdS00+Y3tX21OBK4FP9LugMXY18HrbuwL3AB/tcz1j7XbgCOD6fhcyVlq6n3kX8DrgWEmv629V3WlS8E90ZwIfps3NcBOB7V+2DK7HBNtP21fZXloGf0h1f8uEYXuB7ZW9+3519T/dz9h+ARjqfma116Tgf3/5Gn2BpI36XcxYknQo8Ijt2/pdS50knSHpIeA4Jt4Rf6s/Af6930VER+26n9mqT7WMyoT5IRZJ/wn8TptJpwLnAZ+kOkr8JPBZqv9c40aH/fsY8I7eVjT2VrSPtufYPhU4VdJHgfcDf9/TAldRp/0r85wKLAUu7GVtY6Gb/Ztguup+ZnU0YYLf9tu7mU/Sl6naiMeVkfZP0u8B2wG3SYKqieAWSXvbfrSHJa6ybv+GwEXAXMZZ8HfaP0nTgIOBAzwOb7AZxd9vohi33c80oqlH0hYtg++mOtE0Idj+ie3NbA/YHqD6x7jHeAv9TiRNaRk8FLirX7XUofw40UeAQ23/qt/1RFfGbfczE+aIv4NPS5pK9TVsIXBSf8uJlTBT0s7AS8ADwMl9rmesnQ28Ari6fHP7oe0Js4+S3g18AZgEzJV0q+3/1eeyVkkfup8ZM+myISKiYRrR1BMRES9L8EdENEyCPyKiYRL8ERENk+CPiGiYBH80nqRnx2AdA5J+Lem/JS2QdFO5IStitdOU6/gjeuGntncHkLQ9cJmkNWx/pc91RSwjR/wRbUjaVtK80rHfPEmTy/gdJP1Q0o8k/Z+Rvi3Yvh/4IPBXZbnTJH1N0jWS7pX0vt7tTcSyEvwR7Z0N/EvpH/9C4PNl/FnAWbb3onO/LLcAu7QM7wocBOwHfGKi/ZhMjB8J/oj29qPqDA7ga8AbW8Z/s7y+aPhCwwzvvXGO7V/bfgK4lqo/94ieS/BHdGdl+jbZHViwgnWkv5ToiwR/RHvfp+ptEaoffrmhvP4h8Ifl9THDFxoiaQD4R6qOyYYcJmldSZsA+1P17hjRc+mkLRpP0kss217/OeAy4AJgU2AJcILtB0v30F+nasaZC0y3vVUJ+gVU3UWvCzwDnDd0RY+k04AtgR2AycCnbX+59p2LaCOXc0bj2R7pm+/b2ox7BNjXtiUdA8wv61gIvLLDpu6xPX2lC40YIwn+iNHZEzhbVaf5TzPOfsIzAtLUExHRODm5GxHRMAn+iIiGSfBHRDRMgj8iomES/BERDfP/AeYJkzCMeJ7WAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Volume ratio discrete distribution\n",
    "\n",
    "d1 = 0.01 # Lowest size diameter\n",
    "d_Nb = 1.0 # Diameter of largest bin [microns]\n",
    "Nb = 30 # Number of size bins\n",
    "\n",
    "#------'INSERT CODE HERE'------\n",
    "di = np.exp(np.linspace(np.log(d1), np.log(d_Nb), num=30)) # values for x-axis\n",
    "#log_di = np.linspace(np.log(d1), np.log(d_Nb), num=30)\n",
    "#print(di)\n",
    "\n",
    "V_rat =  # Volume ratio between bins\n",
    "d_rat =  # Diameter ratio between bins\n",
    "\n",
    "# Use the volume ratio to create an array of diameters as follows\n",
    "d_i=np.zeros((Nb), dtype=float) # Diameter array \n",
    "d_i[0]=d1\n",
    "for step in range(Nb):\n",
    "    if step > 0:\n",
    "       d_i[step]= \n",
    "\n",
    "log_di = np.log(d_i) # Log of Diameter array \n",
    "\n",
    "# Create an array of lower bin boundaries\n",
    "vi = (4.0/3.0)*np.pi*np.power(d_i/2.0,3.0)\n",
    "vi_low = 2.0*vi/(1.0+V_rat)\n",
    "di_low = 2.0*(np.power(vi_low/((4.0/3.0)*np.pi),(1.0/3.0)))\n",
    "\n",
    "d_width = \n",
    "\n",
    "# Parameters of log-normal distribution\n",
    "sigmag1 =  # Geometric standard deviation\n",
    "mean1 =  # Mean particle size [150nm]\n",
    "distribution_1 =  # Probability density distribution\n",
    "\n",
    "Ntot = 600.0 # Total number of particles [per cm-3]\n",
    "N_dist =  # Discrete number distribution\n",
    "#------------------------------\n",
    "\n",
    "plt.bar(log_di, N_dist)\n",
    "plt.ylabel(r'N')\n",
    "plt.xlabel(r'LogDp')\n",
    "plt.title(r'Discrete size distribution')\n",
    "plt.show()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objective 3 - Calculate PM2.5 mass  <a name=\"Part3\">\n",
    "    \n",
    "You will have come across the metric PM2.5 when discussing air quality. As the total mass of sampled particulates below 2.5 microns, it has been used historically as a measure of ambient concentrations. Whilst a relatively crude metric, certainly as we learn more about [differential toxicity and size effects on human health](http://www.euro.who.int/__data/assets/pdf_file/0006/189051/Health-effects-of-particulate-matter-final-Eng.pdf), we can play with multiple size distributions and relate these to expected levels of PM2.5. For example, in an urban environment dominated by a large number of fine more aerosol, how does the PM2.5 mass compare with an environmnt where a lower number of accumulation mode aerosol dominate?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<b> Exercise 4. Calculate the PM2.5 mass associated with the size distribution  <a name=\"Exercise3\"> </b> \n",
    "    \n",
    "Using the same distribution as above, calculate PM2.5 mass assuming the particles have a density of 1400kg.m-3 which is a good approximation for secondary organics. Please note your final answer should be given in $\\mu g.m^{-3}$, so you should perhaps check unit coversions on paper. Do you arrive at a mass of ~5.14 $\\mu g.m^{-3}$?\n",
    "\n",
    "Following this, try to find some different size distributions and calculate the total PM2.5 mass. As an added challenge, rather than doing this within the jupyter notebook environment, create a seperate Python script and run this from the command line.\n",
    "\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.142161209310858\n"
     ]
    }
   ],
   "source": [
    "# Total mass in micrograms per cubic metre\n",
    "#------'INSERT CODE HERE'------\n",
    "density = 1400.0 #kg.m-3\n",
    "total_mass = \n",
    "#------------------------------\n",
    "print(total_mass)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}