{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What's the terminal velocity of a raindrop? \n",
    "Nathan Moore, Winona State University\n",
    "\n",
    "This is a \"juypter notebook\" which is essentially a python program with \"markdown\" text output.  Its something like a Mathematica notebook or Matlab program, but free. For an intro, see https://swcarpentry.github.io/python-novice-gapminder/01-run-quit/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# These are some libraries I'm going to use\n",
    "from math import *\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are a few raindrop assumptions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "our raindrop has radius, r =  1.0 mm\n",
      "this raindrop has volume, v =  4.1887902047863905 mm^3\n",
      "this raindrop has mass, m =  0.0041887902047863905 grams\n"
     ]
    }
   ],
   "source": [
    "d=1.0 # water density, grams/cm^3\n",
    "r=1.0 # radius, in mm\n",
    "vol=(4.0/3.0)*pi*r**3 # volume, in mm^3, note 2**3 = 2^3 = 8\n",
    "print(\"our raindrop has radius, r = \",r,\"mm\")\n",
    "print(\"this raindrop has volume, v = \", vol, \"mm^3\")\n",
    "# note, to find mass, (mm^3) * (grams/cm^3) *(1cm/10mm)^3 to make the units work\n",
    "m = vol*d*(1.0/(10.0**3)) # mass (grams)\n",
    "print(\"this raindrop has mass, m = \",m,\"grams\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rain falls from clouds, and is subject to two main forces, weight, $W=mg$, (gravitational attraction from the earth), and a drag force from the air molecules it runs into as it moves, $F_D$.  For sufficiently large raindrops, we can model for drag force with the following term. \n",
    "\n",
    "$F_D=0.5 C d A v^2$,\n",
    "\n",
    "where $C=0.45$ for a sphere, $d=1.2kg/m^3$ is the density of air, $A=\\pi r^2$ is the cross sectional area of the raindrop, and $v$ is the raindrop's speed.  \n",
    "\n",
    "As the drop falls, it gains momentum because of the downward pull of the earth.  As the drop accelerates, the speed increases, and so does the drag force.  Eventually, the weight force from the earth and the drag force from the atmosphere equilibrate, and the drop reaches its *terminal velocity*, $v_T$, which can be computed by equating the two interactions.\n",
    "\n",
    "$mg = 0.5 C d A (v_T)^2$\n",
    "\n",
    "Or,\n",
    "\n",
    "$v_T = \\left[ \\frac{2 m g}{C d A} \\right]^{0.5}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, algebraically, $m=Volume*density_{water}=\\frac{4}{3} \\pi r^3 d_{w}$, and $A=\\pi r^2$, so, the terminal velocity can be expressed as,\n",
    "\n",
    "$v_T = \\left[ \\frac{2 \\left(\\frac{4}{3} \\pi r^3 d_{water}\\right) ~g}{C d_{air} \\pi r^2} \\right]^{0.5} \n",
    "= \\left[ \\frac{ 8 r g}{3C}\\frac{ d_{water} }{ d_{air} } \\right]^{0.5}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the drag coefficient, $C=0.45$ has no units.  The other units in the terminal velocity expression are mixed and we'll have to work them out to plot the terminal speed. Here's a version with mixed units: \n",
    "\n",
    "$v_T = \\left[ \\frac{ 2 (1mm) (9.8m/s^2)}{3C}\\frac{ 1 grams/cm^3 }{  1.2 kg/m^3 } \\right]^{0.5}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a version of the same expression with explicit unit conversions included.\n",
    "\n",
    "$v_T = \\left[ \\frac{ 8 (1mm) (9.8m/s^2)~ \\left(\\frac{1m}{1000mm}\\right) }{3C} \n",
    "\\frac{ 1 grams/cm^3 ~ \\left(\\frac{100cm}{1m}\\right)^3 \\frac{1kg}{1000g} }{  1.2 kg/m^3 } \\right]^{0.5}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated terminal velocity for a r =  1  mm raindrop is  6.956655929999346 m/s\n"
     ]
    }
   ],
   "source": [
    "g=9.8\n",
    "C=0.45\n",
    "r=1 # in mm\n",
    "dw=1 # grams/cm^3\n",
    "da=1.2 #  kg/m^3\n",
    "v_T_sqr=(8*r*g/(3*C*1000))*(dw*100**3/(1000*da))\n",
    "v_T=sqrt(v_T_sqr)\n",
    "print(\"estimated terminal velocity for a r = \",r,\" mm raindrop is \",v_T,\"m/s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so now that we have an expression that seems to work, we'll make a plot that shows the terminal speed as a function of the raindrop size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# plots require a list of points to plot\n",
    "list_of_r=[] # list of radaii, in mm\n",
    "list_of_v=[] # list of terminal velocities, in m/s\n",
    "\n",
    "# this is a for (or for-each) loop.  It repeats the lines below with \n",
    "# i=0, i=1, i=2 ... up to but not including i=30\n",
    "for i in range(30):\n",
    "    # this is the repeated commands tat the forloop runs\n",
    "    \n",
    "    # vary the radius (in mm, incremented by 0.1mm steps)\n",
    "    r=0.1+i*0.1\n",
    "    #print(i,r)\n",
    "    \n",
    "    # save this loop's radius in a list\n",
    "    list_of_r.append(r)\n",
    "    \n",
    "    # using the code above, compute the terminal velocity at this radius\n",
    "    # note that we're using C,g,dw,da defined above\n",
    "    v_T_sqr=(8*r*g/(3*C*1000))*(dw*100**3/(1000*da))\n",
    "    v_T=sqrt(v_T_sqr)\n",
    "    # and save this value as well\n",
    "    list_of_v.append(v_T)\n",
    "\n",
    "# end of tabbed in section means the loop is over"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEeCAYAAACOtbLLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVNX5x/HPwwJL770si3QQBEQEe8GISqJGTSxExURi\nYtToz5ZEY0xMoskvpvxiYhAViYoVjdFojMaC0gSkN+kdlg7LsvX5/TEXnV23zCwze3d2v+/Xa147\nc+eW59wzM8/ee+49x9wdERGRI+qEHYCIiFQvSgwiIlKMEoOIiBSjxCAiIsUoMYiISDFKDCIiUowS\ng4iIFKPEICLVkpkNN7MZZvahmU0xs3phx1RbKDGISHW1ETjL3U8D1gEXhhtO7aHEUIXMbImZnZGA\n9awzs1EJCCnh601UGWPYTh8zm29mB8zs5mRvr5TtV7qcZjbJzB5IcEhJVVX1Gs3dt7p7TvAyDyiK\nZbmq3r9hfxaToW7YAVQlMzsY9bIRkAsUBq+/6+7PJHP77j4gmeuvDqLLaGbrgO+4+ztJ2NSdwHvu\nPjgJ665QbajLaGGW18y6AV8BqmsyDfWzmAy16ojB3ZsceQAbgK9GTYs5KZhZrUqo1VQ3YEkyVhxW\n/epz9WVm1gz4O3Ctu+cnYH3J2MeV+ixW5/quVYmhImbWycxeNrMsM1sbfVgYnGa5y8wWAtlmVjeY\ndoeZLTSzbDN73Mzam9mbwWHlO2bWssQ6RkU9vz1Ydp+ZPW9mDaLmvdvMVgfrWWpmF8cQ/11m9lKJ\naX80sz9VVL5S1tXPzN43s73BaYSvRb3X1cymBuvZZWZ/LllGM/s7kAH808wOmtmdwb56ucR2/mRm\nf4wnBjP7L3Am8Odg3b1LWTau/VdK/d5T3vJx1uUQM5sXrOt5oEE5261b3r6PWuZHQVx7zOzJ6O1V\nUE43s55Rr4uddgli2RzEusLMzq5EeYea2afBOl4M3v/Sf/tm1iyIp03UtP5mtt3Mmlvkh/M54H53\nX1FOmcrcv+Xs4zI/H7HGH8z7pc9iBd+d0mIp7/sU83c2ody9Vj6INGaNinpdB5gL/BSoDxwDrAHO\njZp/PtAVaBg1bSbQHugM7ADmAUOIfDj/C9xX2jaD57OBTkArYBlwQ9S8lwXv1QG+CWQDHUuLPWqZ\nbsAhoGnwOg3YCoyoqHzR6wXqAauAHwfzngUcAPoE61wA/B5oHJTzlHLKGL2POwblaBG8rhvss+NL\nKUuZMQTvv0/kNFVZ9Vvm/ivn8/B5/Va0fKx1GcS+Hrg1KNOlQD7wQBnbLbfcUcssDpZpBXx8ZH0x\nfO4d6Bn1elJULH2INPh2Cl5nAj0qWd5bgrJ8nUj7QKnxETlyPy3q9SvAbcHzbwG7grp+H/hmKcuX\nu3/L+e6WWr/xxl/ys1hR/ZVS32V+n4jhO5u038dkb6C6Pvjyj9aJwIYS8/wIeDJq/utKWcdVUa9f\nBv4a9fom4NXSthk8Hxv13m+AR8uJdz5wYWmxl5jvI+Dq4Pk5wOpYyhe9XuBUYBtQJ+q9KcDPgJFA\nFlC3ov1aWpzAm8D1wfMxwNIy1lNmDMHzz7+MMdb35/uvnLivi3X5WOsSOA3YAljU+9Mpnhiui3qv\n3HJHLRP9T8T5R+o5hv1QXmLoSSRRjwLqVVCv5ZV3c4nyfkTZieFfRNr3jnxG1wPpcdRrufs3lrqN\nrt944y/5WYzhc1uyvsv8PhHDdzZZD51K+kI3oFNw+LfXzPYSyfrto+bZWMpy26Oe55Tyukk529wW\n9fxQ9LxmdrVFrnQ4EsuxQJuSKyjFs8AVwfMrg9cQW/mO6ARsdPfoq0DWEzkq6gqsd/eCGGIpzVPA\n2OD5WCLnj0tTXgwVquT++7x+K7F8WXXZCdjswbc6qhylbpfYy72xxPudyoktJu6+CvghkX8AdpjZ\nc2ZW1nrjKW9p35sjFgP9g+e/IvIDmhtH2LHs3y/FUE79xht/afFUVH/R6yvv+xTPdzahlBi+sBFY\n6+4toh5N3f38qHmqZFQji1yF8RjwA6C1u7cg8gWyGBZ/ETjDzLoAF/NFYoilfEdsAbqaWfTnI4PI\nf1IbgQyLreGstP31KjDIzI4lcsRQVqN/eTGU6yj2nx/l8qXZCnQ2s+hlM0rbbiDWcnct8f6WGOM5\nROSKvCM6FAvE/Vl3P4XIj5IDD8W43iNKK2/XsmYmSAxB+0UHYHICtldy/0LUPq6gfuONv6RY6q9k\n0inr+xTPdzahlBi+MBs4EDQMNTSzNDM71sxOCCGWxkQ+PFkAZjaOyH80FXL3LCKHtk8S+VAtC96K\np3yziPyA3Glm9Sxy/fpXiTQEziby5XnQzBqbWQMzO7mMcLYTOS8aHd9h4CUiCWu2u28oY9nyYqhI\npfdfgpaPNgMoAG4OyvF1YHg588da7hvNrIuZtQJ+AjwfxDrJzCaVs/75wJVB/Y8GTj/yhkWuxz/L\nzNKBw0SOeGO6dyDKDCKXgP8gaFi9kPLLe+SI4VfAT9y9sJx5y9pePPsXyq/feOMvKd7PbXnfp9B+\nk5QYAsEHcgwwGFgL7AQmAs1DiGUp8DsiH9LtwEAiDYyxepbIeeIjRwtxlc/d84h8mM8L5vsLkXaL\n5cF6vkrkfPQGYBORxrvS/Bq4JzgMvj1q+lNBmco6jVRuDBUV/mj3XwL2f/S68og0YF4L7Cayr6ZW\nMH8s5X4WeJtIY+RqvrjGv2sFsd4SrH8vcBWRI7gj0oEHg+1uA9oROacds6jyfjvYxljgdSL3DJVm\nKZEjhUJ3f7WMeWLZ3rXEsH+DZcqs30rEX1o8MX9uy/s+hfmbZMVPpYkkn5llAMuBDu6+P+x4Uo2V\nceOgmdUncoXLIE/ANf+JYmaziDROPxl2LJWR6vFXho4YpEoF515vA55TUkgsd89z935hJwUzO93M\nOgSnYq4BBgFvhRlTPFI9/kSotnfeSc1jZo2JHLqvB0aHHI4kTx/gBSLn8tcAl7r71nBDikuqx3/U\ndCpJRESK0akkEREpRolBRESKSck2hjZt2nhmZmbYYYiIpJS5c+fudPe2Fc2XkokhMzOTOXPmhB2G\niEhKMbPSugv5Ep1KEhGRYpQYRESkGCUGEREpRolBRESKUWIQEZFilBhERKSYKksMZvaEme0ws8VR\n035rZsstMqj4K2bWoqriERGR0lXlEcMkvtxx2n+AY919ELCSOPt+FxGpLT7bfoBfvrGU3dl5Sd9W\nld3g5u4fmllmiWlvR72cCVxaVfGIiFR3B3MLeH3BFp6fs5FPN+ylXppxYvfWjOqf3GGfq9Odz9cR\nDE9YGjMbD4wHyMgobUhXEZHU5+7MXb+H5z/ZyBuLtnIor5Ce7Zrwk/P7cfHQzrRpkp70GKpFYjCz\nnxAZt7WsgeFx9wnABIBhw4apr3ARqVGyDuQydd4mXpizkdVZ2TSqn8ZXB3XiGyd0ZWhGC8ysymIJ\nPTGY2bVExjU92zU4hIjUIoVFzocrs3jukw28u2wHBUXO8d1a8ptLenDBoI40Tg/nJzrUxGBmo4E7\ngdPd/VCYsYiIVJVNew7xwpxNvDhnI1v3HaZ14/qMOzmTb57QlZ7tmoYdXtUlBjObApwBtDGzTcB9\nRK5CSgf+ExwmzXT3G6oqJhGRqpJfWMS7y7YzZfZGPvwsC4BTe7Xlp2P6c3a/9tSvW31uK6vKq5Ku\nKGXy41W1fRGRMKzdmc3zn2zkpbmb2Hkwlw7NGnDTmT25bFhXurZqFHZ4pQq9jUFEpKY5nF/Iv5ds\nY8rsDcxcs5u0OsZZfdtx+QldOb13W+qmVZ+jg9IoMYiIJMjqrINMmbWBl+dtYs+hfLq2asgd5/bh\n0uO70L5Zg7DDi5kSg4jIUcgtKOStxV8cHdStY3xlQHuuGJ7ByT3aUKdO1V1mmihKDCIilbAm6yDP\nBW0Hu7Pz6NqqIXeOjhwdtGuaOkcHpVFiEBGJUW5BIW8v2c6zszYwY80u6tYxzukfOTo4pWdqHh2U\nRolBRKQCG3cf4tnZG3jhk43sys6jS8tI28Flx3ehXQq1HcRKiUFEpBSFRc5/l+/gmVnr+WBlFgaM\n6teeK0/M4LRebWvM0UFplBhERKLs2H+Y5z/ZyJTZG9iy7zDtmqZz01m9uGJ4Vzo2bxh2eFVCiUFE\naj13Z8bqXTw9az1vL9lOQZFzaq82/PSrkbuS61Xz+w4STYlBRGqtfTn5vDx3E0/PWs+arGxaNKrH\nuJMzufLEbnRv0zjs8EKjxCAitc6SLft4euZ6Xv10Czn5hQzJaMHD3ziO8wd2pEG9tLDDC50Sg4jU\nCrkFhby5aBt/n7meuev30KBeHS4a3JmxI7pxbOfmYYdXrSgxiEiNtnlvDs/MXM/zwaWmma0bcc8F\n/bjs+K40b1Qv7PCqJSUGEalx3J2PVu1k8oz1vLtsOwBn9W3P1SO71agb0ZJFiUFEaoyDuQW8PHcT\nT81Yx5qsbFo3rs8Np/fgyhMz6NKyenZxXR0pMYhIyluTdZDJM9bz0txNHMwt4LguzXn4G8dxwaCO\npNdVY3K8lBhEJCUVFTnvr9zBpOnr+XBlFvXSjDGDOnHNSZkM7toi7PBSmhKDiKSUfTn5vDhnI5Nn\nrGfD7kO0b5bObef05orhGbRtmh52eDWCEoOIpIRVOw4yafpaXp67mZz8Qk7IbMmdo/tw7oAOte7O\n5GRTYhCRaquoyPngsyye/HgdH67Mon7dOlx4XCeuPTmTAZ1070GyKDGISLWTnVvAy/M2MenjdazZ\nmU27pun8zzm9ufLEDFo30emiZFNiEJFqY+PuQzw1fR3Pz9nIgcMFHNe1BX+8fDDnHduR+nV1uqiq\nKDGISKjcnVlrd/PER2t5Z9l2zIzzB3Zk3MmZDM1oGXZ4tZISg4iEIq+giDcWbWHitLUs2bKflo3q\nccPpPfjWyG61ZtyD6kqJQUSq1N5DeTwzawOTZ6xj+/5cerRtzK8uHsjFQzrTsL5uRqsOlBhEpEqs\nyTrIEx+v5aW5mzicX8Spvdrw4CWDOL2GD5OZipQYRCRp3J0Za3bx+LS1vLt8B/XT6nDRkE5cd0p3\n+nZoFnZ4Uoa4E4OZNQYOu3thEuIRkRogv7CI1xdu4bEP17J0635aN67PLWf3YuyIbro7OQVUmBjM\nrA5wOXAVcAKQC6Sb2U7gDeBv7r4qqVGKSErYfzif52Zv4MmP17F132F6tmvCg18fyEVDOmtktBQS\nyxHDe8A7wI+Axe5eBGBmrYAzgYfM7BV3fzp5YYpIdbZlbw6Tpq/j2VkbOJhbwMhjWvOriwdyem+1\nH6SiWBLDKHfPLznR3XcDLwMvm5mGQRKphZZs2cfEaWv554ItOHDBwI5cf+oxDOyi7ipSWYWJ4UhS\nMLPLgLfc/YCZ3QsMAR5w93mlJY6SzOwJYAyww92PDaa1Ap4HMoF1wDfcfU8lyyIiVcDdmfbZTiZ8\nuIaPVu2kUf00rh6ZyXWnZGownBoinsbne939RTM7BTgb+C3wV+DEGJefBPwZmBw17W7gXXd/0Mzu\nDl7fFUdMIlJFCgqLeH3hVh79YDXLtx2gXdN07hrdlyuHZ2js5BomnsRw5CqkC4AJ7v6GmT0Q68Lu\n/qGZZZaYfCFwRvD8KeB9lBhEqpVDeQW88MlGHpu2ls17c+jVrgm/vXQQFw7urP6Laqh4EsNmM/sb\ncA6RBud04Gg/Fe3dfWvwfBvQvqwZzWw8MB4gIyPjKDcrIhXZk53HUzPW8dT0dew5lM+wbi25/2sD\nOKtvOzUo13CxXK46EpgJfAMYDfyvu+81s47AHYkKxN3dzLyc9ycAEwCGDRtW5nwicnQ27TnExGlr\nef6TjeTkFzKqXztuOL0HwzJbhR2aVJFYjhiuBh4BVgJvAQcAgv/0t5azXCy2m1lHd98aJJodR7k+\nEamkZVv387cPVvPPhVsx4KIhnRl/2jH0bt807NCkisVyVdL3AMysL3AeMMnMmhO5v+Et4OOjuAv6\nNeAa4MHg7z8quR4RqaS563fzyHur+e/yHTSun8a4kzL59qnd1cNpLRZzG4O7LweWA783s4ZEbm67\nDHgYGFbR8mY2hUhDcxsz2wTcRyQhvGBm3wbWEzldJSJJduSS00feW8Wstbtp2aget53Tm2tGZuoK\nI6lcJ3rungP8K3jEuswVZbx1dmViEJH4FRU5by/dxiPvrWbR5n10aNaAe8f054rhXWlUX31qSkTM\nnwQzGwb8BOgWvZy7D0pCXCKSQPmFRbw2fwt//WA1q3YcpFvrRjz49YFcPLQz6XXVh5EUF8+/CM8Q\nuQppEVCUnHBEJJEO5xfy4pyNPPrBGjbvzaFvh6b86YohnH9sB+qm6R4EKV08iSHL3V9LWiQikjCH\n8gp4ZuYGJkxbQ9aBXIZmtOAXFw3gzD7tMNM9CFK+eBLDfWY2EXiXSNfbALj71IRHJSKVcjC3gMkz\n1jFx2lp2Z+dxUo/W/OnyIYw4ppUSgsQsnsQwDugL1OOLU0kOKDGIhGxfTj6TPl7HEx+vZV9OPqf3\nbsvNZ/fk+G66KU3iF09iOMHd+yQtEhGJ257sPB7/aC1PTV/HgdwCRvVrz01n9eS4ri3CDk1SWDyJ\nYbqZ9Xf3pUmLRkRisvNgLo9NW8PTM9aTnVfIecd24Adn9WRAJ42DIEcvnsQwAphvZmuJtDEYkS6O\ndLmqSBXZeTCXv32wmr/PXE9eQRFjBnXiB2f1VLcVklDxJIbRSYtCRMq182AuEz5cw99nrCe3oJCL\nBnfmxrN60qNtk7BDkxoolt5VzSPWVzRPYkMTkV0Hc5kwbQ2Tp0cSwoWDO/MDJQRJsliOGN4zs5eB\nf7j7hiMTzaw+cAqRzu/eIzJCm4gkwO7sPCZ8uIbJM9aRk1/I147rxE1n9aJnOyUESb5YEsNo4Dpg\nipl1B/YCDYA04G3gD+7+afJCFKk99mTn8di0NTw1fR2H8gv56qBO3Hx2T3q2UxuCVJ1Yut0+DPwF\n+IuZ1QPaADnuvjfZwYnUFvty8pk4bQ1PfLSWQ/mFXDCwI7ec3YtealSWEMTVnaK753P0g/OISCA7\nt4BJ09fxtw9Ws/9wQSQhjOqlq4wkVOpnVyQEh/MLeXrmev76/mp2Zecxql87bjunD/07NQs7NBEl\nBpGqlFdQxAtzNvLn/65i2/7DnNKzDbd9pTdDM1qGHZrI5yqdGIIxmne7e26FM4vUcoVFziufbuaP\n765k4+4cju/Wkt9/czAje7QOOzSRLzmaI4a/Az3M7GV3vz1RAYnUJEVFzpuLt/Hwf1awOiubYzs3\n4+fjjuWM3m3V26lUW5VODO4+yiKf7P4JjEekxpj2WRa/eWsFizbvo1e7Jjw6dijnDuighCDVXtyJ\nwcyud/fHINJRErAk4VGJpLCFm/by0FvL+XjVLjq3aMjvLjuOi4Z0Jq2OEoKkhsocMQxJeBQiNcCa\nrIP87u2VvLFoK60a1+enY/pz1YgMjaksKacyiWG0mT0GzA0eC9UALbXZ9v2H+cM7n/HCnI2k163D\nzWf34vpTu9O0Qb2wQxOplMokhn8DPwaOB84CfghclcigRFLBvpx8Hv1gNU9+vJbCIudbI7px45k9\nads0PezQRI5KZRLDLHffA7wTPERqlcP5hUyesY5H3lvN/sP5XHhcJ247pw8ZrRuFHZpIQsScGMzs\nY+Aed5+UvHBEqq+iIuefC7fwm7dWsHlvDqf3bstdo/vqbmWpceI5YvgucL+Z3UMkQcxIUkwi1c7M\nNbv41b+WsXDTPvp3bMZvLh3EyT3bhB2WSFLEnBjcfTFwiZkNBX4eXIv9E3dfkKzgRMK2ascBHnxz\nOe8s20HH5g14+BvHcdHgztTRpadSg1WmjWEV8AsiYzTMreQ6RKq1rAO5/OGdlTz3yUYa1kvjztF9\nuO7k7jSop0tPpeaLp43hPaAXkAMsDR7XJicskXDk5BUycdoaHv1gNbkFRYw9MYObz+5F6ya60khq\nj3j+2/8fYJm75yQrGJGwFAWd3P3m38vZvj+Xcwe0567RfTlGYytLLRRPG8O8ZAYiEpa56/fw89eX\nsmDjXgZ1ac6frxzKCZmtwg5LJDTVon3AzG4FvgM4sAgYFwwpKpI0W/bm8OCby3ltwRbaNU3nd5cd\nx8VD1LAsEnpiMLPOwM1Af3fPMbMXgMuBSaEGJjXWobwCHv1gDRM+XI073HRWT244vQeN00P/OohU\nC9VloJ66QEMzywcaAVsSsE6RYoqKnNcWbOHBN5ezbf9hxgzqyN3n9aVLS92xLBIt9IF63H2zmf0v\nsIHIFU9vu/vbRxGXyJd8umEP9/9zKfM37mVg5+b835VD1I4gUobQB+oxs5bAhUB3YC/wopmNdfen\nS8w3HhgPkJGRcTSblFpkx/7DPPjmcqZ+upm2TdP57aWDuGRoF7UjiJQjnvsYbgKeDjrQAxI2UM8o\nYK27ZwXbmQqcBBRLDO4+AZgAMGzYMD/KbUoNl19YxFPT1/GHdz4jr6CI75/RgxvP7Kl2BJEYxPMt\naQ98YmbzgCeAfweJ4WhtAEaYWSMip5LOBuYkYL1SS01fvZP7/rGEz3Yc5Iw+bbnvqwPo3qZx2GGJ\npIx47mO4x8zuBb4CjAP+HFxB9Li7r65sAO4+y8xeAuYBBcCnBEcGIvHYui+HB95YxhsLt9K1VUMe\nu3oYo/q10xjLInGK67ja3d3MtgHbiPyItwReMrP/uPudlQ3C3e8D7qvs8lK75RYU8vhHa/m/d1dR\n5M4PR/XihtN7qF8jkUqKp43hFuBqYCcwEbjD3fPNrA7wGVDpxCBSWR+szOL+15awZmc2X+nfnnvH\n9KdrK11+KnI04jliaAV83d3XR0909yIzG5PYsETKt3lvDve/toS3l26ne5vGTBp3Amf0aRd2WCI1\nQjyJoUHJpGBmD7n7Xe6+LMFxiZSqoLCIJz9ex+/fWUmRO3ec24fvnNqd9Lo6bSSSKPEkhnOAu0pM\nO6+UaSJJ8emGPfz4lcUs27qfs/u242dfG6DTRiJJUGFiMLPvAd8HjjGzhVFvNQU+TlZgIkfsy8nn\nt/9ezjOzNtC+aQMeHTuUcwd00NVGIkkSyxHDs8CbwK+Bu6OmH3D33UmJSgRwj/Rt9IvXl7E7O5dr\nT8rkf77Shya6SU0kqSr8hrn7PmAfcEXywxGJWLczm3v/sZhpn+1kUJfmTBp3Asd2bh52WCK1Qiyn\nkj5y91PM7ACR8RI+f4vIrQ3Nkhad1Dq5BYVM+GAN//feKuqn1eH+rw1g7IhupKlvI5EqE8sRwynB\n36bJD0dqs3kb9nDXSwv5bMdBLhjYkZ9+tT/tmzUIOyyRWqdOrDOa2VNm1iLqdUszeyI5YUltciiv\ngJ//cymX/HU62bkFPHntCTxy1VAlBZGQxNOKN8jd9x554e57zGxIEmKSWuSjz3Zy99SFbNqTw9Uj\nu3Hn6L5qXBYJWTzfwDpm1vJIt9tm1irO5UU+t+9QPr/811JemLOJ7m0a88J3RzK8uwbOEakO4vlh\n/x0ww8xeDF5fBvwy8SFJTffW4m3c+4/F7M7O43tn9OCWs3upwzuRaiSebrcnm9kc4Kxg0tfdfWly\nwpKaKOtALj97bQlvLNpK/47NePJaXYIqUh3F2+32UkDJQOLi7kydt5mfv76UnLxC7ji3D+NPO4Z6\naTFf+yAiVSiuxGBmxwGnBi+nufuCxIckNcmOA4f58dTFvLNsO8d3a8lDlwyiZ7smYYclIuWIdzyG\n64GpwaSnzWyCu/9fUiKTlPevRVv5ySuLyM4r5J4L+jHu5O66UU0kBcRzxPBt4ER3z4ZIl9vADECJ\nQYrZeyiP+15bwj/mb2FQl+Y8/I3j6NlO90eKpIp4EoMBhVGvC4NpIp97b8UO7n55IbsO5nHrqN58\n/8weaksQSTHxJIYngVlm9gqRhHARoDufBYCDuQX88o1lTJm9gV7tmjDx6hMY2EVXHImkonguV33Y\nzN4HTg4mXePu85MSlaSUWWt2cftLC9i0J4fvnnYMt57TW/cliKSwWHpXje5V1Uq8p95Va7HD+YX8\n779X8PjHa+nashEvfHckJ2Tq7mWRVBdL76pqNZQvWbHtADdP+ZQV2w9w1YkZ/Pj8fjRWH0ciNUI8\nl6sacBXQ3d1/YWZdgY7uPjtp0Um14+78feZ6HnhjGc0a1OXJcSdwZp92YYclIgkUz794fwGKiHSJ\n8QvgIPAIcEIS4pJqaHd2Hne+tJB3lm3n9N5t+d/LjqNt0/SwwxKRBIsnMZzo7kPN7FP4vNvt+kmK\nS6qZ6at2cusL89mTnc+9Y/oz7qRM6uhmNZEaKZ7EkG9maQQN0WbWlsgRhNRg+YVFPPyflTz6wWq6\nt2nM49eo4zuRmi6exPAn4BWgnZn9ErgUuCcpUUm1sH5XNjdP+ZQFm/ZxxfCu3DumP43qq4FZpKaL\n5z6GZ8xsLnA2wQ1u7r4saZFJqKbO28S9ry4mrY7xl6uGcv7AjmGHJCJVJJ6rkm4Dnnf3R5IYj4Qs\nO7eAe15dzCufbmZ4Zit+f/lgOrdoGHZYIlKF4jkv0BR428x2A88DL7r79uSEJWH4bPsBbnh6Lmt3\nZnPrqN784Kye6g1VpBaKuXczd7/f3QcANwIdgQ/M7J1EBGFmLczsJTNbbmbLzGxkItYrsfvH/M18\n7c8fsy8nn6e/cyK3jOqlpCBSS1WmJXEHsA3YBSTqzqY/Am+5+6XBJbCNErReqUBuQSEPvL6Mv89c\nz/DMVvzflUNo36xB2GGJSIjiaWP4PvANoC3wInB9IsZ8NrPmwGnAtQDungfkHe16pWKb9hzixmfm\nsWDTPsafdgx3nNtHXWSLSFxHDF2BHyahR9XuQBbwZDB06FzgliMDAklyvLdiB7c+P5/CQufRsccz\n+tgOYYckItVEPG0MP0pSN9t1gaHAX919CJAN3F1yJjMbb2ZzzGxOVlZWEsKoHQqLnIffXsF1kz6h\nQ7MGvHbTKUoKIlJMdThvsAnY5O6zgtcvEUkUxbj7BHcf5u7D2rZtW6UB1hS7DuZyzROz+dN/V3Hp\n0C68euPJdG/TOOywRKSaCf02VnffZmYbzayPu68gcgPdUbddSHHzNuzh+0/PY/ehPB66ZCDfPCEj\n7JBEpJqPBopAAAARU0lEQVQKPTEEbgKeCa5IWgOMCzmeGuWluZv48dRFtG+eztTvnaS+jkSkXPGO\n4FbsLSAhI7gFbRfDjnY9UlxhkfPQW8uZ8OEaTurRmkeuHErLxuoQV0TKpxHcaqgDh/O5ecqnvLci\ni6tHduPeMf11KaqIxCSuU0lm1hLoBXx+B5S7f5jooOTorNuZzXcmz2HdzmweuOhYxo7oFnZIIpJC\n4rnB7TvALUAXYD4wAphBZEQ3qSamr9rJ95+dB8Dkbw/npB5tQo5IRFJNPOcWbiEyjOd6dz8TGALs\nTUpUUil/n7GObz0xm3ZN03ntxlOUFESkUuI5lXTY3Q+bGWaW7u7LzaxP0iKTmOUXFvGz15bwzKwN\nnN23HX+4fDBNG9QLOywRSVHxJIZNZtYCeBX4j5ntAdYnJyyJ1Z7sPL73zFxmrtnNDaf34I5z+6hX\nVBE5KvGM4HZx8PRnZvYe0Bx4KylRSUzWZB3k2ic/Ydv+w/z+m8dx8ZAuYYckIjVApW5wc/cPEh2I\nxGf+xr1cN+kTDHhu/AiGZrQMOyQRqSHiuSopHbgEyIxezt1/nviwpDzvLd/B95+ZR9um6Uy+bjiZ\n6u9IRBIoniOGfwD7iHSLnZuccKQiL87ZyN1TF9GvY1OevHY4bZumhx2SiNQw8SSGLu4+OmmRSLnc\nnb+8v5rf/nsFp/Zqw1/HHk+T9OrS1ZWI1CTx/LJMN7OB7r4oadFIqQqLnPv/uYTJM9Zz0eBO/ObS\n46hfV91biEhyxJMYTgGuNbO1RE4lHelEb1BSIhMADucXcuvz83lz8TbGn3YMd4/uSx1djioiSRRP\nYjgvaVFIqfbl5HP95DnMXrubey7ox3dOPSbskESkFojnPgbdzFaFtu07zDVPzGbNzoP88fLBXDi4\nc9ghiUgtEct4DB+5+ylR4zJEn8dIyHgMUtyqHQe4+vHZ7D9cwKRxwzm5p/o8EpGqE8t4DKcEfzUu\nQxVYse0AVz42EzPjufEjNNqaiFS5eG5wGwb8mC/f4KbG5wQ5khTS6hhTxo+gR9smYYckIrVQPI3P\nzwB3AIuAouSEU3st37afKx+bRb00Y8r1IzhGSUFEQhJPYshy99eSFkkttnTLfq6aOJP0umlMGT+C\n7uriQkRCFE9iuM/MJgLvEtUlhrtPTXhUtciRpNCgXhpTrh+hfo9EJHTxJIZxQF+gHl+cSnJAiaGS\nlmzZx1UTZ9GoXuRIoVtrJQURCV88ieEEd9eIbQmyeHMkKTRJr8uU60eQ0bpR2CGJiADxjfk83cz6\nJy2SWmTRpi+SwnPjlRREpHqJ54hhBDBffSUdnYWb9jJ24iyaNqjHc+NH0LWVkoKIVC8xJQYzM+C7\naIzno7Jg416+9fgsmjWMJIUuLZUURKT6iSkxuLub2SPuPjDZAdVUy7ftZ+zjs2jRqB5TrldSEJHq\nK542hnlmdkLSIqnBtuzN4donPqFx/bpKCiJS7cXTxnAicJWZrQeyURtDTPbl5HPtk7PJzi3gxe+N\nVFIQkWovnsRwbtKiqKFyCwoZP3kOa3dm89R1w+nbQR3Rikj1p/EYkqSoyLnthQXMWrubP14+mJN6\nqOtsEUkNMbcxWMRYM/tp8DrDzIYnL7TU9qt/LeONhVv50Xl9NciOiKSUeBqf/wKMBK4IXh8AHklU\nIGaWZmafmtnriVpnWB7/aC0TP1rLtSdlMv40DccpIqklnsRworvfCBwGcPc9QP0ExnILsCyB6wvF\n6wu38MAbSxk9oAP3julP5BYQEZHUEU9iyDezNCId52FmbUnQuAxm1gW4AJiYiPWFZeaaXdz2/AKO\nz2jJHy4fTFodJQURST3xJIY/Aa8A7czsl8BHwK8SFMcfgDspJ9GY2Xgzm2Nmc7KyshK02cRZuf0A\n4yfPoWurhky8ZhgN6qWFHZKISKXEkxgGEfnx/jWwFbgIGHa0AZjZGGCHu88tbz53n+Duw9x9WNu2\nbY92swm1bd9hrnliNun10pg0bjgtGiXyDJuISNWK5z6Gc9z9LmD5kQlmdh5w11HGcDLwNTM7H2gA\nNDOzp9197FGut0rsPxy5gW1/Tj4v3DBSneKJSMqr8IjBzL5nZouAPma2MOqxFlh4tAG4+4/cvYu7\nZwKXA/9NlaRQVOT88Ln5rNpxkEe/dTwDOjUPOyQRkaMWyxHDs8CbRE4h3R01/YC7705KVCnisWlr\n+O/yHdz/tQGc2qt6nd4SEamsChODu+8D9vHF/QtJ4+7vA+8nezuJMHf9bn7z7xWcd2wHrh7ZLexw\nREQSJp7GZwnsyc7jpmc/pVOLBjx06SDdqyAiNUo8jc8CuDu3v7iArIO5vPy9k2jWoF7YIYmIJJSO\nGOI0cdpa3l2+gx+f349BXVqEHY6ISMIpMcRh3oY9PPTWcs4d0J5rT8oMOxwRkaRQYojR3kORdoUO\nzRvwm0uPU7uCiNRYamOIQaRdYSE7DhzmxRtOonlDtSuISM2lI4YYPP7RWt5Ztp27z+vH4K5qVxCR\nmk2JoQLzN+7lobeWc07/9lx3cmbY4YiIJJ0SQzn2Hcrnxmfm0a5pA36r+xVEpJZQG0MZ3J07XlrA\n9v2HefGGkeoxVURqDR0xlOHJj9fx9tLt3H1eX4ZktAw7HBGRKqPEUIrVWQf59ZvLGNWvHd8+pXvY\n4YiIVCklhlI88PpS0uum8euvq11BRGofJYYS3l+xg/dWZHHTWT1p2zQ97HBERKqcEkOU/MIifvH6\nUrq1bsS1ujRVRGopJYYoT89cz+qsbH5yfj/S66aFHY6ISCiUGAJ7svP4wzufcXLP1pzTv33Y4YiI\nhEaJIfCHd1Zy4HA+947prwZnEanVlBiAldsP8PSsDVx5YgZ9OzQLOxwRkVDV+sTg7vzi9aU0rp/G\nbef0CTscEZHQ1frE8N/lO5j22U5uGdWbVo3V7YWISK1ODHkFRTzwxjKOaduYq0d2CzscEZFqoVYn\nhskz1rF2Zzb3XtCfemm1eleIiHyu1v4a7jqYyx/f/YzTe7flzL7twg5HRKTaqLWJ4eH/rORQXiH3\njukXdigiItVKrUwMy7buZ8rsDXxrRDd6tmsadjgiItVKrUsMRy5PbdawHj8c1SvscEREqp1alxje\nXrqd6at3ceuo3hqVTUSkFLUqMeQWFPKrfy2jV7smXHViRtjhiIhUS7UqMTz58TrW7zrEvWP6U1eX\np4qIlKpW/Tq2bZLOZcd34bTebcMORUSk2qobdgBm1hWYDLQHHJjg7n9MxrYuOb4LlxzfJRmrFhGp\nMUJPDEAB8D/uPs/MmgJzzew/7r407MBERGqj0E8luftWd58XPD8ALAM6hxuViEjtFXpiiGZmmcAQ\nYFYp7403szlmNicrK6uqQxMRqTWqTWIwsybAy8AP3X1/yffdfYK7D3P3YW3bqvFYRCRZqkViMLN6\nRJLCM+4+Nex4RERqs9ATg0UGWH4cWObuD4cdj4hIbRd6YgBOBr4FnGVm84PH+WEHJSJSW4V+uaq7\nfwRY2HGIiEiEuXvYMcTNzLKA9SUmtwF2hhBOstS08kDNK1NNKw/UvDLVtPLA0ZWpm7tXePVOSiaG\n0pjZHHcfFnYciVLTygM1r0w1rTxQ88pU08oDVVOm6tDGICIi1YgSg4iIFFOTEsOEsANIsJpWHqh5\nZapp5YGaV6aaVh6ogjLVmDYGERFJjJp0xCAiIgmgxCAiIsWkVGIws9FmtsLMVpnZ3aW8b2b2p+D9\nhWY2NIw44xFDmc4ws31Rd4X/NIw4Y2VmT5jZDjNbXMb7KVVHMZQn1eqnq5m9Z2ZLzWyJmd1Syjyp\nVkexlCll6snMGpjZbDNbEJTn/lLmSW4duXtKPIA0YDVwDFAfWAD0LzHP+cCbRO6kHgHMCjvuBJTp\nDOD1sGONo0ynAUOBxWW8n2p1VFF5Uq1+OgJDg+dNgZU14HsUS5lSpp6C/d4keF6PyDAEI6qyjlLp\niGE4sMrd17h7HvAccGGJeS4EJnvETKCFmXWs6kDjEEuZUoq7fwjsLmeWlKqjGMqTUjy2gbFSrY5q\n1GBfwX4/GLysFzxKXiWU1DpKpcTQGdgY9XoTX678WOapTmKN96TgcPFNMxtQNaElTarVUSxSsn7K\nGRgrZeuovMG+SKF6MrM0M5sP7AD+4+5VWkehd6InFZoHZLj7waDX2VeBXiHHJF9IyfqpaGCsVFRB\nmVKqnty9EBhsZi2AV8zsWHcvtZ0rGVLpiGEz0DXqdZdgWrzzVCcVxuvu+48cVrr7v4B6Ztam6kJM\nuFSro3KlYv3EMDBWytVRRWVKxXoCcPe9wHvA6BJvJbWOUikxfAL0MrPuZlYfuBx4rcQ8rwFXBy32\nI4B97r61qgONQ4VlMrMOZmbB8+FE6mxXlUeaOKlWR+VKtfoJYq1oYKyUqqNYypRK9WRmbYMjBcys\nIXAOsLzEbEmto5Q5leTuBWb2A+DfRK7mecLdl5jZDcH7jwL/ItJavwo4BIwLK95YxFimS4HvmVkB\nkANc7sFlCdWRmU0hcgVIGzPbBNxHpPEsJesohvKkVP3wxcBYi4Jz2AA/BjIgNeuI2MqUSvXUEXjK\nzNKIJLAX3P31qvytU5cYIiJSTCqdShIRkSqgxCAiIsUoMYiISDFKDCIiUowSg4iIFKPEICIixSgx\niIhIMUoMUq2Y2fQ45/+Zmd2erHjiVTKeeMuT4FgamtkHwY1SiVxvfTP70MxS5gZZiY8Sg1Sp4Bb+\nMj937n5SsrdRletLRHnKEsMP/nXA1KBDtoQJuoh/F/hmItcr1YcSgySdmWVaZJS6ycBioKuZvWpm\ncy0yQtX4qHkPRi2zzMweC+Z5O+g3BjP7iZmtNLOPgD7lbOM2M1scPH4YNd9yM3smWP9LZtboKGP+\nUjzR5QnWtThq2u1m9rPgeWMze8Mio3UtNrNyf2zN7EUz+5uZzQR+VMGuvwr4R4lyTwpifcbMRpnZ\nx2b2mZkNj2WeqHW/GqxfaqJEjvqjhx6lPYBMoIioUaiAVsHfhkR+eFsHrw9GLVMADA5evwCMBY4H\nFgGNgGZE+oq5veQ2ouZrDDQBlhDppz+TyKAnJwfzPQHcXtmYy4onapmDwboWR027HfhZ8PwS4LGo\n95pXsC+XAz+PYZ/XB7aVKE8BMJDIP4Rzg7IbkUFfXo1lnqj1pQFZYX+29EjOQ0cMUlXWe2SkqSNu\nNrMFwEwi3QeX1jf+Wnc/0inaXCI/XKcCr7j7IY/0uR/dG230Nk4J5sv2SHfLU4NlATa6+8fB86eD\neSsbc3nxxGIRcI6ZPWRmp7r7vrJmNLMGQCvg5zGstw2wt8S0te6+yN2LiCTKd93dgxgy45gHj5ye\nyjOzprEUUlKLEoNUlewjT8zsDGAUMNLdjwM+BRqUskxu1PNCKu4NOLuC948o2XNkWT1JVibm0hRQ\n/Lv2+XLuvpLImNKLgAes/EHqBxAZ27cghm3mlBJf9P4sinpdxBf7NpZ5jkgHDscQi6QYJQYJQ3Ng\nj7sfMrO+RAYzj9WHwEXBFTdNga+WMd+0YL5GZtYYuDiYBpBhZiOD51cCHx1FzLHEsx1oZ2atzSwd\nGHPkDTPrBBxy96eB3xJJEpjZu2ZWcqjGgcDC6AllzIe77wHSgqOMhDOz1sBOd89PxvolXLrcTMLw\nFnCDmS0DVhA5NRMTd59nZs8DC4iMh/tJOfNNAmYHkya6+6cWGRN4BXCjmT0BLAX+WtmYY4nH3fPN\n7OdBLJspPujKQOC3ZlYE5BMZM6AO0BPYXWJVA6PKQznzHfE2kdNk78RQvnidCbyRhPVKNaDxGKRW\nCRLD6+5+bMihlMnMjgWuc/fbjmY+MxsK3Oru30pCjFOBu4NTYVLDKDFIrZIKiSGRzOw64ClP4L0M\nFgxD6+6TE7VOqV6UGEREpBg1PouISDFKDCIiUowSg4iIFKPEICIixSgxiIhIMUoMIiJSjBKDiIgU\n8//jI9f+MV3ZuwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x4db8ee0438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# now, make a plot\n",
    "plt.plot(list_of_r,list_of_v)\n",
    "plt.ylabel(\"terminal velocity, $v_T$, (m/s)\")\n",
    "plt.xlabel(\"raindrop radius, $r$, (mm)\")\n",
    "plt.title(\"Terminal velocity of a raindrop, using $v^2$ drag force\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Note, the internet tells me that most raindrops are 4mm diameter (2mm radius) or smaller.  https://water.usgs.gov/edu/raindropsizes.html\n",
    "\n",
    "See also, figure 3 in `http://journals.ametsoc.org/doi/pdf/10.1175/1520-0477(2001)082<1169%3AWIARSD>2.3.CO%3B2`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# How much pressure does a raindrop exert on the ground?\n",
    "\n",
    "A moving raindrop comes to a stop when it hits the ground because there's a contact interaction between the rain and the ground.  I've heard that erosion is caused by wind and water, so let's investigate this with an estimate.  How big is the stopping force or pressure that exists between the raindrop and the ground?\n",
    "\n",
    "Before hitting the ground, the raindrop has a downward (terminal) velocity and a corresponding downward momentum.  The ground exerts an external force on the raindrop for a length of time, and the momentum disappears.  Mathematically, \n",
    "\n",
    "$P_1=mv_T$ and $P_2=0$ because $P_2-P_1=F_{external}~\\Delta t$.\n",
    "\n",
    "How long does it take the raindrop to splash?  Some engineers at Harvard recently investigated how Mosquitoes survive rainstorms (their conclusion, rain isn't a big deal for Mosquitoes).  In the course of that paper, they estimated that the raindrop-Mosquito interaction time is about 1ms. Here's the paper, http://www.pnas.org/content/109/25/9822.full\n",
    "\n",
    "So then, we can estimate the raindrop-ground force as,\n",
    "\n",
    "$F=\\frac{m~v_T}{\\Delta T}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, remember from above that:\n",
    "\n",
    "$m=\\frac{4}{3} \\pi r^3 d_{w}$ and $v_T = \\left[ \\frac{ 8 r g}{3C}\\frac{ d_{water} }{ d_{air} } \\right]^{0.5}$\n",
    "\n",
    "We could combine these into a single formula, but I don't think the resulting expression is particularly informative.  Instead, remember that the stopping force is spread over the ground, making a pressure, $P=F/A$.  The simplest estimate for the relevant area is that the force is spread over is the cross-sectional area of the raindrop, $A=\\pi r^2$.  This gives an expression for raindrop pressure,\n",
    "\n",
    "$P=\\frac{F}{A}=\\frac{m~v_T}{\\Delta T~\\left(\\pi r^2\\right)}=\n",
    "\\frac{\\left(\\frac{4}{3} \\pi r^3 d_{w}\\right)~v_T}{\\Delta T~\\pi r^2} =\n",
    "\\frac{4  ~r ~d_{w}~v_T}{3 ~\\Delta T}$\n",
    "\n",
    "$P=\\frac{4  ~r ~d_{w}}{3 ~\\Delta T} \\left[ \\frac{ 8 r g}{3C}\\frac{ d_{water} }{ d_{air} } \\right]^{0.5}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The units are again mixed up.  A different way to resolve this is to have two sets of variables in python code - one that you use in the measurement (eg $r=1mm$), and another in standard mks units that you use in the calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inputs:\n",
      "raindrop radius r =  0.001  m \n",
      "impact time dt =  0.001  s\n",
      "water density d_w =  1000.0  kg/m^3\n",
      "air density d_air =  1.2  kg/m^3\n",
      "Spherical drag coefficient, C =  0.45\n",
      "g =  9.8\n",
      "Computed pressure, P =  9275.541239999127 Newtons/m^2\n",
      "Computed pressure, P =  1.34530298951655 psi\n"
     ]
    }
   ],
   "source": [
    "delta_T_ms = 1 # ms\n",
    "delta_T = delta_T_ms*(1.0/1000.0) # seconds\n",
    "#\n",
    "r_mm = 1 # mm\n",
    "r = r_mm*(1.0/1000.0) # meters  \n",
    "#\n",
    "d_w_gpcm3 = 1 # grams per cm^3\n",
    "d_w = d_w_gpcm3*(1.0/1000.0)*(100.0/1.0)**3 # kg/m^3\n",
    "#\n",
    "d_air = 1.2 # kg/m^3$\n",
    "C = 0.45\n",
    "g = 9.8 # m/s^2\n",
    "P=(4.0*r*d_w/(3.0*delta_T))*(8.0*r*g*d_w/(3*C*d_air))**0.5\n",
    "# echoing out all of the calculation parameters as a check of the unit conversions \n",
    "print(\"Inputs:\")\n",
    "print(\"raindrop radius r = \",r,\" m \")\n",
    "print(\"impact time dt = \",delta_T,\" s\")\n",
    "print(\"water density d_w = \",d_w,\" kg/m^3\")\n",
    "print(\"air density d_air = \",d_air,\" kg/m^3\")\n",
    "print(\"Spherical drag coefficient, C = \",C)\n",
    "print(\"g = \",g)\n",
    "print(\"Computed pressure, P = \",P,\"Newtons/m^2\")\n",
    "# note, Newtons/m^2 is a \"Pascal (Pa)\", and 6894.76Pa = 1psi\n",
    "P_psi = P*(1.0/6894.76)\n",
    "print(\"Computed pressure, P = \",P_psi,\"psi\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "So, how does the Pressure depend on raindrop size? We'll use the looping structure from above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# plots require a list of points to plot\n",
    "list_of_r=[] # list of radaii, in mm\n",
    "list_of_p=[] # list of pressures, in psi\n",
    "\n",
    "# this is a for (or for-each) loop.  It repeats the lines below with \n",
    "# i=0, i=1, i=2 ... up to but not including i=30\n",
    "for i in range(30):\n",
    "    # this is the repeated commands tat the forloop runs\n",
    "    \n",
    "    # vary the radius (in mm, incremented by 0.1mm steps)\n",
    "    r_mm=0.1+i*0.1\n",
    "    #print(i,r)\n",
    "    \n",
    "    # save this loop's radius in a list\n",
    "    list_of_r.append(r_mm)\n",
    "    \n",
    "    # convert the radius to (standard) units of m\n",
    "    r = r_mm*(1.0/1000.0) # meters  \n",
    "\n",
    "    # using the code above, compute the terminal velocity at this radius\n",
    "    # note that we're using C,g,dw,da defined above\n",
    "    P=(4.0*r*d_w/(3.0*delta_T))*(8.0*r*g*d_w/(3*C*d_air))**0.5\n",
    "    P_psi = P_psi = P*(1.0/6894.76)\n",
    "    # and save this value as well\n",
    "    list_of_p.append(P_psi)\n",
    "\n",
    "# end of tabbed in section means the loop is over"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEaCAYAAAASSuyNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VFX6wPHvSxJI6AqI9NCkBgKEohTFriCIDRVcEQVF\nV3T5udgX1sWOrqIsiKIUEVRcRRFcUUEWFULovbdQQyehpby/P+5NdggpNyGTSXk/zzNPZuaee+97\nz0zeuXPumXNEVTHGGFP0lQh0AMYYY/KHJXxjjCkmLOEbY0wxYQnfGGOKCUv4xhhTTFjCN8aYYsIS\nvgkIEVkjIlflct0JIjIij0MymRCRcBFREQkOdCz+JiLbReTaQMfhL5bwc8l9Y5wSkXgR2e8mobKB\njquwUNVmqjov0HEUJAXlg6yoJ73izBL+hblFVcsCrYEo4IX0BcQR0HrO7zOzQJ4JFtazUBEJCnQM\nBZHVS96yhJ8HVHU3MBtoDiAi80TkZRH5DTgJ1BORCiIyXkT2ishuERmR+mYWkQYi8quIHBORgyLy\nufu8iMg/ReSAiBwXkVUi4ruPh1JjEJF+IrLA57GKyGMisgnY5D7XWETmiMhhEdkgIndldkwi8oCI\nrBOREyKyVUQezqJsPxH5zY31EDBcROqLyC8icsg9pikiUtFnnbSzSBEZLiJfiMgkd39rRCTKp2wr\nEVnqLvscCPVZdpWIxIrI0yKyD/jEfX6AiGx2j/VbEamerm4Gu8d1UETezOxDWURKiMgzIrLFPZYv\nRORid1lvEdkmIuXdxzeJyD4RqZJdfbtn82NEZJaIJAAPAn2Aoe63xu/cck+775cT7jauySTOCm79\nxYnIDhF5IfWYUt8bIjJSRI64Md+UyXYmA7WB79w4hvos7iMiO906e95LHWWyj6Hu/8EeEXnIfT0a\nZFIvXbM5tuEi8qnPts9pfhLn/+Qf7vvzhIj8KCKVfcrf527zkO8xFVmqardc3IDtwLXu/VrAGuAf\n7uN5wE6gGRAMhABfAx8AZYBLgGjgYbf8VOB5nA/gUKCT+/wNwBKgIiBAE6Cazz4e8omnH7DA57EC\nc4CLgTB3v7uAB9yYWgEHgaaZHF83oL673ytxPrhaZ1K2H5AEPO5uOwxoAFwHlAKqAPOBdzKpv+HA\naeBmIAh4FVjoLisJ7AD+4tbjHUAiMMJdfpW779fdfYUBV7vH1tp97j1gfrq6mevWTW1go29dpju2\nJ4CFQE13Wx8AU32WTwEmAJWAPUB39/ks69td5xjQ0ed1n5B6XG6ZRu42qruPw4H6mcQ5CZgBlHPL\nbQQe9Hl9EoEBbv0OcmOV7N7bPvtV4EO3flsCZ4AmXuoo3bZvBPbh/G+UBj51t90gi3rJ6tiGA59m\nEGuwz//JFuAyN/Z5wGvusqZAPNDFjfttnPfStRnFXhRuAQ+gsN7cf4p44ChOQvoXEOYumwe85FO2\nqvsPEubz3D3AXPf+JGAcUDPdPq5239wdgBLpls0j+4R/tc/j3sB/023jA2CYx+P9Bngik2X9gJ3Z\nrH8rsCxd/fkm/J98ljUFTrn3u5AuOQG/c27CPwuE+iwfD7zh87gsTsIL96mbG32WPwr8nEnc64Br\nfB5Xc7eVmlAq4ny4rwI+8FrfOIltUrrlEzg34TcADgDXAiFZ1G2QWwdNfZ57GJjn8/ps9llW2q2D\nS7N4b2eU8Gv6PBcN3O2ljtJt+2Pg1XTHmD7hT8rBsQ0n+4T/QrrX+gf3/t+AaT7Lyrj7KrIJ35p0\nLsytqlpRVeuo6qOqespn2S6f+3Vwzk73ishRETmK889/ibt8KM6ZdLTbnNEfQFV/Ad4HRgMHRGRc\navOBR+ljaJ+6fzeGPsClGa3oNk8sdJsjjuKcfVfOqGwG+0JEqorINLc54jjOmVxW6+/zuX8SCHW/\nllcHdqv7H+nakW7dOFU97fO4um8ZVY0HDgE1Mol3h7tORuoAX/vU2TogGedDHFU9CnyJ05z3Vrr1\nsqvvc+osPVXdDDyJk9QOuPWZUZyVcd5fvvWyg3OPN61+VfWkezennQzSv0ap62dZR+lU59zjzqgO\nfJ/zcmy5jfucWFQ1Aed9UmRZwvcf3wS1C+cMv7L7AVFRVcurajMAVd2nqgNUtTrO2cu/Uts0VXWU\nqrbBOeu9DPiru80EnDO1VBkl7vQx/Oqz/4qqWlZVB6VfSURKAV8BI4GqqloRmIXzoeTleAFecZ+L\nUNXyQN9s1s/MXqCGiPiuWzubfe/BSUIAiEgZnCaX3T5laqXb3p5M9r8LuCldvYWqc90GEYkE+uM0\ny41Kt1529Z0+7vOGrlXVz1S1k3s8itN0ld5BnDPqOj7P1U53vDmR0yF0s6yjdPbiNP2kqpVBGd/9\nZ3dsXv4PMrPXd/8iUhrnfVJkWcLPB6q6F/gReEtEyrsXueqLyJUAInKniKT+ExzBecOniEhbEWkv\nIiE4b+zTQIpbbjlwm4iUdj8cHswmjJnAZe5FqhD31lZEmmRQtiROm2YckORe4Ls+h4ddDqfJ65iI\n1OB/H1Q59QdOu+pgN+bbgHbZrDMVeEBEIt0Pr1eARaq63afMX0XkIhGphdMG/Xkm2xoLvCwidQBE\npIqI9HTvh+J8c3kOp62+hog86q6Xk/pOtR+ol/pARBqJyNXuMZwGTvG/1z+NqiYDX7hxlnNjHeLG\nlhvnxOFBpnWUgS9wXpsmboJ9MasNezi25UAXEaktIhWAZ3MQ93Sgu4h0EpGSwEsU8ZxYpA+ugPkT\nTiJdi5PUp+O0dQK0BRaJSDzwLU5b+VagPM6FsiM4X2MPAW+66/wTp71xPzAR5+JhplT1BE7Svhvn\nbHYf/7vQmVHZwTj/aEeAe924cuLvOBdNjwHfA//O4fqpsZwFbsNphz6M0zae5bZU9SecRPIVzllc\nfZzj9jUD54L4cje+8Zls7l2cY/9RRE7gXJxs7y57FdilqmNU9QzOt5gRItIwJ/XtYzzQ1G0a+cYt\n+xrOWe4+nCbAzBLa4zgnBVuBBcBnOO3lufEq8IIbx1MeymdVR+dQ1dk434TmApvdsuB8A85Mpsem\nqnNwPqxX4ryeMz3EmxrLGuAxd3t7cd7rsV7XL4zk3KZRY4o+EVGgodtGbgLI/cazGiilqkmBjqeo\nszN8Y0y+EpFeIlJKRC7C+dbznSX7/GEJ3xiT3x7G6W66Bac3z3kdB4x/WJOOMcYUE3aGb4wxxUSB\nGmiqcuXKGh4eHugwjDGm0FiyZMlBVa3ipWyBSvjh4eHExMQEOgxjjCk0RCT9L88zZU06xhhTTFjC\nN8aYYsISvjHGFBOW8I0xppiwhG+MMcWEJXxjjCkmLOEbY0wxUaD64efa7Gdg36pAR2GMMblzaQTc\n9Jrfd2Nn+MYYU0wUjTP8fPhkNMaYws7O8I0xppiwhG+MMcWEJXxjjCkmLOEbY0wx4beELyKNRGS5\nz+24iDzpr/0ZY4zJmt966ajqBiASQESCgN3A1/7anzHGmKzlV5PONcAWVfU8UL8xpnA6k5RMYnJK\noMMwGcivhH83MDWf9mWMCRBV5bl/r+a+8YtIsqRf4Pg94YtISaAH8GUmyweKSIyIxMTFxfk7HGOM\nH33y23a+WhpLh3qVCA6yPiEFTX68IjcBS1V1f0YLVXWcqkapalSVKp7m4TXGFEC/bT7Iy7PWcUOz\nqgy+umGgwzEZyI+Efw/WnGNMkbbz0Eke+2wp9auU4a27IilRQgIdksmAXxO+iJQBrgP+7c/9GGMC\nJ+FMEgMmxaAKH/4pirKlisYQXUWRX18ZVU0AKvlzH8aYwElJUYZ8sZxNB04wqX976lQqE+iQTBbs\nqooxJtfe+2Uz/1mzn+dubkKnhpUDHY7JhiV8Y0yu/LhmH//8aSO3tarBg53qBjoc44ElfGNMjm3c\nf4K/fL6cljUr8MptEYjYRdrCwBK+MSZHjp48y4BJMYSVDGbsfW0IDQkKdEjGI0v4xhjPEpNTeOyz\npew5eooP7mtNtQphgQ7J5ID1nzLGeDZi5lp+23yIN+5oQZs6Fwc6HJNDdoZvjPHks0U7mfjHDh7q\nVJe7omoFOhyTC5bwjTHZWrj1EH+bsZorL6vCszc3CXQ4Jpcs4RtjsrTr8EkGfbqE2pVKM+qeVgTZ\nsAmFliV8Y0ym4s8k8dDEGJJTlPH3t6VCWEigQzIXwC7aGmMylJKiPDltOZvj4pnwQFvqVrZhEwo7\nO8M3xmRo5I8b+Gndfl7s1oTODW3o8qLAEr4x5jwzlu/mX/O2cE+7Wtx/RXigwzF5xBK+MeYcS3ce\n4a/TV9Iu/GL+3qO5DZtQhFjCN8ak2X30FAMnLeHS8qGMva8NJYMtRRQldtHWGAM4PXIenLCYM4nJ\nTB3QnovLlAx0SCaPWcI3xpCcojw5bRkb95/gkwfa0bBquUCHZPzAvq8ZY3j9h/X8tO4Aw25pxpWX\nWY+cosrfc9pWFJHpIrJeRNaJyOX+3J8xJuc+X7yTcfO3cl+HOtYjp4jzd5POu8APqnqHiJQESvt5\nf8aYHFi49RAvfLOazg0rM+yWpoEOx/iZ3xK+iFQAugD9AFT1LHDWX/szxuTMjkMJPPLpEmpfXJr3\n721NcJC18BZ1/nyF6wJxwCciskxEPhIR+222MQXAsVOJ9J+wGMDGyClG/Jnwg4HWwBhVbQUkAM+k\nLyQiA0UkRkRi4uLi/BiOMQbcWaumLGXn4ZOM7duGcBsjp9jwZ8KPBWJVdZH7eDrOB8A5VHWcqkap\nalSVKtY7wBh/UlX+NmM1CzYf5JVeEXSoVynQIZl85LeEr6r7gF0i0sh96hpgrb/2Z4zJ3rj5W5ka\nvYvHutbnTpu1qtjxdy+dx4Epbg+drcADft6fMSYTs1ft5dXZ6+neohr/d12j7FcwRY5fE76qLgei\n/LkPY0z2lu86ypOfL6d17YqMvLMlJWzWqmLJ+mEZU8TFHjnJQxNjuKR8KT78UxShIUGBDskESLZn\n+CISBXQGqgOngNXAHFU94ufYjDEX6Phpp/vlmaRkpg1sT6WypQIdkgmgTM/wReQBEVkKPAuEARuA\nA0An4CcRmSgitfMnTGNMTqV2v9wal8AHfdvQ4BIbEK24y+oMvzTQUVVPZbRQRCKBhsBOfwRmjMk9\nVWXYt2v476aDvHF7C65oUDnQIZkCINOEr6qjs1rRvSBrjCmAxs3fymeLdjLoqvrc1da6XxpHpglf\nRIaq6hsi8h6g6Zer6mC/RmaMyZXvVuzh1dnruaVldf56vXW/NP+TVZPOOvdvTH4EYoy5cIu3H+b/\nvlhB2/CLePOOFtb90pwjqyad79y/E1OfE5ESQFlVPZ4PsRljcmBLXDwDJsVQ8+Iw635pMpRtP3wR\n+UxEyrsjXa4G1orIX/0fmjHGq4PxZ+j3STRBIkzo146KpW0+WnM+Lz+8auqe0d8KzMYZ9vg+v0Zl\njPHs1NlkHpwYQ9yJM4zv15balWyeIZMxLwk/RERCcBL+t6qaSAYXcY0x+S85RXli2jJWxh5l1N2t\niKxVMdAhmQLMS8L/ANgOlAHmi0gdwNrwjSkA/jFzLT+u3c+w7k25vtmlgQ7HFHDZDq2gqqOAUT5P\n7RCRrv4LyRjjxfgF25jw+3Ye6lSXfh3rBjocUwh4uWhbSURGichSEVkiIu8CFfIhNmNMJmau3MOI\n79dyU/NLee7mJoEOxxQSXpp0puHMTXs7cId7/3N/BmWMydzCrYcY8vkKoupcxD97R1pfe+OZl/Hw\nq6nqP3wejxCR3v4KyBiTuQ37TjBgUgy1K5W2vvYmx7yc4f8oIneLSAn3dhfwH38HZow5195jp+j3\nSTSlSwYxsb/1tTc55yXhDwA+A864t2nAwyJyQkSst44x+eDYqUT6fbyYE6eT+KRfO2pUDAt0SKYQ\n8tJLxwbRNiaAziQlM3BSDFsPxjPhgXY0rV4+0CGZQiqr0TLDVXV7FssFqKGqsVmU2Q6cAJKBJFW1\n+W2NyYGUFGXIFytYtO0w794dSUcb195cgKzO8N90B0ubASzB6Z0TCjQAugLXAMOATBO+q6uqHsyD\nWI0pdl6ZtY7vV+7luZsb0zOyRqDDMYVcVqNl3ikiTYE+QH+gGnASZ9jkWcDLqno6X6I0phj6cP5W\nPlqwjX5XhDOgc71Ah2OKgCzb8FV1LfD8BWxfcea/TQY+UNVx6QuIyEBgIEDt2jZFrjEAXy2J5eVZ\n6+gWUY0XuzfFaUE15sJ46aVzITqpaiRwE/CYiHRJX0BVx6lqlKpGValSxc/hGFPwzV1/gKFfraRj\ng0q83bslQfbDKpNH/JrwVXW3+/cA8DXQzp/7M6awW7rzCI9OWUqTauUY27cNpYLth1Um7/gt4YtI\nGREpl3ofuB5nAhVjTAY2HzhB/wmLuaR8KT7p145yoSGBDskUMV6GVjiHiFQDDqvqmWyKVgW+dtse\ng4HPVPWHnIdoTNG35+gp7hsfTUhQCSb3b0+VcqUCHZIpgnKc8IHJQH0R+UpVn8qskKpuBVrmOjJj\niomjJ8/yp4+jiT+dxLSHO9iMVcZvcpzwVfVa90dXTf0QjzHFyqmzyfSfsJidh08y8YF2NKtuI48b\n//HUhi8inUTkAfd+ZSBcVdf4NTJjirjE5BQenbKEZbuOMuruSC6vXynQIZkizssEKMOAp4Fn3adK\nAp/6MyhjirqUFOWpL1cwd0McI25tzo3NqwU6JFMMeDnD7wX0ABIAVHUPYAOqGZNLqsrw79YwY/ke\n/npDI/q0rxPokEwx4SXhn1VVxfnVbGoXS2NMLv3zp01M+mMHA7vU49Gr6gc6HFOMeEn4X4jIB0BF\nERkA/AR86N+wjCmaxi/YxqifN9E7qhbP3tTYhkww+crLePgjReQ64DjQCPibqs7xe2TGFDFfLYnl\nHzOdicdfuS3Ckr3Jd1kmfBEJAn5S1a6AJXljcunHNfsY+tVKOjWozDt3R9r4OCYgsmzSUdVkIEVE\nrHOwMbn0+5aD/HnqMiJqVOCD+2x8HBM4Xn54FQ+sEpE5uD11AFR1sN+iMqaIWBl7lAETYwivVJpP\n+rWlTKnc/LjdmLzh5d33b/dmjMmBDftO8KePo7moTEkmP9iei8qUDHRIppjzctF2Yn4EYkxRsu1g\nAn3HL6JUcAk+e6gDVcuHBjokY7JP+CKyDbcPvi9VtTnXjMnA7qOn6PvRIpJTlKk2GJopQLw06UT5\n3A8F7gQu9k84xhRuB06cps+HCzl+OpGpAzrQ4BL7UbopOLL94ZWqHvK57VbVd4Bu+RCbMYXKkYSz\n3PdRNAdOnGHCA+1oXsM6t5mCxUuTTmufhyVwzvitq4ExPk6cTuT+T6LZdiiBCf3a0qbORYEOyZjz\neEncb/ncTwK2A3f5JRpjCqFTZ5N5cEIMa/cc54P72nBFg8qBDsmYDHnppdM1PwIxpjA6k5TMw58u\nIWbHYd69uxXXNKka6JCMyZSX8fCfEJHy4vhIRJaKyPVedyAiQSKyTERmXlioxhQsZ5NSeGzKUuZv\njOO121pwS8vqgQ7JmCx5GS2zv6oeB64HKgH3Aa/lYB9PAOtyEZsxBVZScgpPTFvGT+sO8I+ezbir\nba1Ah2RMtrwk/NRRnm4GJrlTG3oa+UlEauL06Pkod+EZU/Akpyh/+WIFs1fv48XuTbnv8vBAh2SM\nJ14S/hIR+REn4f9HRMoBKR63/w4wNKvyIjJQRGJEJCYuLs7jZo0JjJQUZej0lXy3Yg9P39iYBzvV\nDXRIxnjmJeE/CDwDtFXVk0AI8EB2K4lId+CAqi7JqpyqjlPVKFWNqlKlipeYjQmIlBTl+W9W8dXS\nWP5y7WUMstmqTCHjJeFfDmxQ1aMi0hd4ATjmYb2OQA8R2Q5MA64WEZv83BRKqfPQTo3exWNd6zP4\nmgaBDsmYHPOS8McAJ0WkJfB/wBZgUnYrqeqzqlpTVcOBu4FfVLXvhQRrTCCoKi9/v45Jf+xgQOe6\nPHV9I5utyhRKXhJ+kjuJeU/gfVUdDdgAIaZYUFXe+M8GPlqwjX5XhPPczU0s2ZtCy8svbU+IyLM4\n3TE7i0gJnHZ8z1R1HjAvx9EZE0CqysgfNzBm3hbubV+bYbc0tWRvCjUvZ/i9gTM4/fH3ATWBN/0a\nlTEBlprsR8/dwj3tajOiZ3NL9qbQ8zJa5j7gK6CU+9RB4Gt/BmVMIKkqb/240U32tXj51uaUsEnH\nTRHgZbTMAcBAnDHw6wM1gLHANf4NzbvXo19n/eH1gQ7DFBG7jpxk95FT1GkWyr7SZXjwx0BHZIq6\nxhc35ul2T/t9P16adB7D6WJ5HEBVNwGX+DMoYwIlNdlfUi6UelXKBDocY/KUl4u2Z1T1bGr7pYgE\nk8GUh4GUH5+MpmhTVf45ZyNzV2zm7ra1eKVXhDXjmCLHyxn+ryLyHBAmItcBXwLf+TcsY/KPqvLP\nnzYx6pfN9I6yZG+KLi8J/xkgDlgFPAzMwvm1rTGFnqry9pyNjPp5E3dF1eTV2yzZm6IryyYdEQnC\nGSGzD/Bh/oRkTP5QVV77YT0f/LqV3lG1LNmbIi/LM3xVTQbqiEjJfIrHmHyhqvxj5jo++HUrfTvU\ntmRvigUvF223Ar+JyLdAQuqTqvq236Iyxo9SUpRh365h8sIdPNAxnL91t1/QmuLBS8Lf4t5KYGPo\nmEIudYjjqdG7eLhLPZ65qbEle1NseJnE/O8AIlLeeagn/B6VMX6Q7E5e8tXSWP7ctQH/d/1lluxN\nseJlEvMoEVkFrARWicgKEWnj/9CMyTtJySkM+WJ52uQlT91gQxyb4sdLk87HwKOq+l8AEekEfAK0\n8GdgxuSVs0kp/OXz5Xy/ai9/vaERj3W1yUtM8eQl4SenJnsAVV0gIkl+jMmYPHM6MZlHpyzll/UH\neP7mJgzoUi/QIRkTMF4S/q8i8gEwFWdIhd7APBFpDaCqS/0YnzG5lnAmiYcmxrBw2yFe7tWcPu3r\nBDokYwLKS8Jv6f4dlu75VjgfAFfnaUTG5IFjJxPpNyGalbHHePuulvRqVTPQIRkTcF566XTNj0CM\nySuH4s9w3/hoNh04weh7W3Nj80sDHZIxBYKXM/xcEZFQYD7OxCnBwHRVTf8twZg8te/Yafp8tJDd\nR0/x0f1tufKyKoEOyZgCw28JH2daxKtVNV5EQoAFIjJbVRf6cZ+mGNt56CR9xi/kSEIiEx9oR/t6\nlQIdkjEFit8SvqoqEO8+DHFvBWocfVN0bD5wgj4fLeJMUgpTHmpPy1oVAx2SMQWOlykOQ4FHgU44\nCXsBMEZVT3tYNwhYAjQARqvqogzKDMSZQpHatWvnKHhjAFbsOkq/T6IJKlGCaQM70PjS8oEOyZgC\nyct4+JOAZsB7wPtAU2Cyl42rarKqRgI1gXYi0jyDMuNUNUpVo6pUsfZWkzO/bT7IvR8upEypYKY/\ncrkle2Oy4KVJp7mqNvV5PFdE1uZkJ6p6VETmAjcCq3OyrjGZ+WH1XgZPXU7dymWY9GA7qpYPDXRI\nxhRoXs7wl4pIh9QHItIeiMluJRGpIiIV3fthwHXA+twGaoyvadE7eXTKUiJqVuCLhy+3ZG+MB17O\n8NsAv4vITvdxbWCDO6CaqmpmY+pUAya67fglgC9UdeYFR2yKvTHztvD6D+u58rIqjOnbmtIl/dnZ\nzJiiw8t/yo252bCqrsT5Na4xeUJVeXX2esbN30qPltUZeWdLSgZ7+ZJqjAFvv7TdISItgc7uU/9V\n1RX+DcuYcyUlp/Dsv1fx5ZJY7r+8DsNuaWZTEhqTQ17Gw38CmAJc4t4+FZHH/R2YMalOnk3i4clL\n+HJJLE9c05DhPSzZG5MbXpp0HgTaq2oCgIi8DvyB003TGL86FH+GByfGsDL2KCNubU7fDjbipTG5\n5SXhC5Ds8zjZfc4Yv9p1+CR/+jiaPUdPMaZvG25oZoOgGXMhvCT8T4BFIvK1+/hWYLz/QjIGVu8+\nRr9PFpOY7AyVEBV+caBDMqbQ83LR9m0RmYcztALAA6q6zK9RmWLtv5vieGTyEiqWLsm0ge1pcEm5\nQIdkTJGQZcJ3+9CvUdXGgM1sZfzu62Wx/PXLlTS4pCwTHmjHpRXsB1XG5JUse+moajLOj6xsVDPj\nV6rKmHlb+MvnK4gKv4gvHrnckr0xecxLG/5FwBoRiQYSUp9U1R5+i8oUK0nJKfzt2zV8tmgn3VtU\n4627WlIqOCjQYRlT5HhJ+C/6PQpTbJ04ncifP1vGrxvjGHRVff56fSPrY2+Mn3i5aPtrfgRiip+9\nx07xwCeL2XQgnldvi+CedtZyaIw/ZZrwReQEWcxQpao28LjJtTV7jtF/wmISziTzcT+be9aY/JBp\nwlfVcgAi8g9gL86kJwL0wRkJ05hcmbv+AH/+bCnlw0L48pHLaVLNzh2MyQ9e2vB7qGpLn8djRGQF\n8Dc/xWSKsMkLdzBsxmqaVCvPx/3a2jj2xuQjLwk/QUT6ANNwmnjuwae3jjFeJKcor8xax/gF27i6\n8SW8d08rypSyceyNyU9e/uPuBd51bwr85j5njCcnTicyeOoy5m6Io98V4bzQrQnBQTaOvTH5zUsv\nne1AT/+HYoqinYdO8tCkxWyJS7DRLo0JsGwTvohUAQYA4b7lVbW//8IyRUH0tsM88ukSklOUyf3b\ncUWDyoEOyZhizUuTzgzgv8BPnDtMcpZEpBYwCaiK0xQ0TlXfzU2QpvD5YvEunv9mFbUuLs34+9tS\nt3KZQIdkTLHnJeGXVtWnc7HtJOD/VHWpiJQDlojIHFVdm4ttmUIiOUV5bfY6PvzvNjo3rMz797am\nQlhIoMMyxuBhikNgpojcnNMNq+peVV3q3j8BrANq5HQ7pvA4fjqRAZNi+PC/27j/8jp80q+tJXtj\nChAvZ/hPAM+JyBkgEefHV5qTX9qKSDjQCliUwbKBwECA2rXtp/WF1eYD8QycHMPOQyf5R89m3Hd5\neKBDMsak46WXzgXNPiEiZYGvgCdV9XgG2x8HjAOIiorKdCgHU3D9tHY/T36+nFLBJZjyUHva16sU\n6JCMMRnIaiydxqq6XkRaZ7Q8tbkmKyISgpPsp6jqv3MfpimIUlKU937ZzD9/2kiLmhUY27cN1SuG\nBTosY0xorrI0AAAgAElEQVQmsjrDH4LT1PJWBssUuDqrDYuI4Mx9u05V3851hKZAOnE6kf/7YgU/\nrt3P7a1r8nKv5oSG2Bj2xhRkWQ2eNtD92zWX2+4I3AesEpHl7nPPqeqsXG7PFBBb4+IZOHkJ2w4m\nMOyWpvS7Ihzn890YU5B5GsxERJoDTYG0ka5UdVJW66jqApwLvKYI+WX9fp6YtpyQoBJ8+mB7Lq9v\n7fXGFBZefmk7DLgKJ+HPAm4CFuD8qMoUE8kpyj/nbOT9uZtpXqM8Y/u2oeZFpQMdljEmB7yc4d8B\ntASWqeoDIlIV+NS/YZmC5GD8GQZPXcbvWw5xd9taDO/RzNrrjSmEvCT8U6qaIiJJIlIeOADU8nNc\npoCI2X6Yxz5bytGTibx5RwvujLKX3pjCykvCjxGRisCHwBIgHvjDr1GZgFNVPv5tO6/OWkfNi8L4\n5NF2NK1uM1MZU5hlmfDdrpWvqupRYKyI/ACUV9WV+RKdCYgTpxN5+quVzFq1j+ubVmXkXS0pH2pD\nJARaYmIisbGxnD59OtChmAAIDQ2lZs2ahITk/n8xy4Svqiois4AI9/H2XO/JFArr9h7nsSlL2XH4\nJM/d3JgBnetZl8sCIjY2lnLlyhEebt1gixtV5dChQ8TGxlK3bt1cb8fL4GlLRaRtrvdgCgVVZfIf\n2+k5+jfizyTx2UPtGdilviWWAuT06dNUqlTJXpNiSESoVKnSBX+789KG3x7oIyI7cOayTR08rcUF\n7dkUGMdOOk04P6zZx1WNqjDyzpZULlsq0GGZDFiyL77y4rX3kvBvuOC9mAJryY7DDJ66nP3HT/P8\nzU14sFNdSpSwpGJMUZRtk46q7sjolh/BGf9JSVFGz93MXR8spEQJmD7oCgZ0qWfJ3mQpKCiIyMhI\nmjdvzi233MLRo0ezXeeKK67I0T7mzZtH9+7dcxtilsLDw4mIiKBFixZcf/317Nu3zy/7Kai8tOGb\nIubAidP86eNo3vzPBm5sfinfD+5MZK2KgQ7LFAJhYWEsX76c1atXc/HFFzN69Ohs1/n999/zZN9J\nSUl5sp25c+eycuVKoqKieOWVV85bnpzseSbXC5af+wKPY+mYouOntft5+quVJJxN4rXbIujdtpa1\nCxdCf/9uDWv3nDe9xAVpWr08w25p5rn85ZdfzsqVTg/t+Ph4evbsyZEjR0hMTGTEiBH07NkTgLJl\nyxIfH8+8efMYPnw4lStXZvXq1bRp04ZPP/0UEeGHH37gySefpHTp0nTq1CltH8OHD2fLli1s3bqV\n2rVr88knnzBo0CBiYmIIDg7m7bffpmvXrkyYMIGvv/6aY8eOsXv3bvr27cuwYcOyjL9Lly6MGjUq\nLcaHH36Yn376idGjRxMWFsaQIUOIj4+ncuXKTJgwgWrVqjFq1CjGjh1LcHAwTZs2Zdq0afz66688\n8cQTgNPOPn/+fJYsWcLIkSOZOXMmAH/+85+JioqiX79+hIeH07t3b+bMmcPQoUNp27Ytjz32GHFx\ncZQuXZoPP/yQxo0be3/hcsASfjGRcCaJEd+vZWr0LppUK8+7d0dyWdULmtvGFGPJycn8/PPPPPjg\ng4DTR/zrr7+mfPnyHDx4kA4dOtCjR4/zTiaWLVvGmjVrqF69Oh07duS3334jKiqKAQMG8Msvv9Cg\nQQN69+59zjpr165lwYIFhIWF8dZbbyEirFq1ivXr13P99dezceNGAKKjo1m9ejWlS5embdu2dOvW\njaioqEyPYebMmURERACQkJBA+/bteeutt0hMTOTKK69kxowZVKlShc8//5znn3+ejz/+mNdee41t\n27ZRqlSptOaskSNHMnr0aDp27Eh8fDyhoaGZ7jNVpUqVWLrUmVLkmmuuYezYsTRs2JBFixbx6KOP\n8ssvv3h8JXLGEn4xsHTnEYZ8vpwdh0/yyJX1+ct1DSkVbGPhFGY5ORPPS6dOnSIyMpLdu3fTpEkT\nrrvuOsDp1vvcc88xf/58SpQowe7du9m/fz+XXnrpOeu3a9eOmjVrAhAZGcn27dspW7YsdevWpWHD\nhgD07duXcePGpa3To0cPwsKciXUWLFjA448/DkDjxo2pU6dOWsK/7rrrqFTJGb31tttuY8GCBRkm\n/K5duxIUFESLFi0YMWIE4FybuP322wHYsGEDq1evTju25ORkqlWrBkCLFi3o06cPt956K7feeisA\nHTt2ZMiQIfTp04fbbrst7fiykvqhFh8fz++//86dd96ZtuzMmTPZrp9blvCLsMTkFN7/ZTPvz93M\npeVDmTagg00/aC5Iahv+yZMnueGGGxg9ejSDBw9mypQpxMXFsWTJEkJCQggPD8+wz3ipUv/r7hsU\nFOSpXb5MmTKeYkv/bSKzpsq5c+dSuXLlc54LDQ0lKMg5CVJVmjVrxh9/nD+CzPfff8/8+fP57rvv\nePnll1m1ahXPPPMM3bp1Y9asWXTs2JH//Oc/BAcHk5KSkrZe+rpIPaaUlBQqVqzI8uXLyQ920baI\n2nYwgTvG/sG7P2+iZ8vqzH6ysyV7k2dKly7NqFGjeOutt0hKSuLYsWNccsklhISEMHfuXHbs8N6R\nr3Hjxmzfvp0tW7YAMHXq1EzLdu7cmSlTpgCwceNGdu7cSaNGjQCYM2cOhw8f5tSpU3zzzTd07Ngx\nV8fWqFEj4uLi0hJ+YmIia9asISUlhV27dtG1a1def/11jh07Rnx8PFu2bCEiIoKnn36atm3bsn79\neurUqcPatWs5c+YMR48e5eeff85wX+XLl6du3bp8+eWXgPNhs2LFilzF7YWd4RcxKSnKp4t28Oqs\n9ZQMLsHoe1vTrUW1QIdliqBWrVrRokULpk6dSp8+fbjllluIiIggKioqRxcdQ0NDGTduHN26daN0\n6dJ07tyZEydOZFj20UcfZdCgQURERBAcHMyECRPSvjW0a9eO22+/ndjYWPr27Ztl+31WSpYsyfTp\n0xk8eDDHjh0jKSmJJ598kssuu4y+ffty7NgxVJXBgwdTsWJFXnzxRebOnUuJEiVo1qwZN910E6VK\nleKuu+6iefPm1K1bl1atWmW6vylTpjBo0CBGjBhBYmIid999Ny1btsxV7NkRVfXLhnMjKipKY2Ji\nAh1GobXjUAJDp69k0bbDdLmsCm/c3oJLK2R/AckUDuvWraNJkyaBDqNAmjBhAjExMbz//vuBDsWv\nMnoPiMgSVfX06ea3M3wR+RjoDhxQ1eb+2o9xzuon/rGdN37YQHAJ4Y3bW3BnVE3rbmmMOYc/m3Qm\nAO9jUyH61baDCQydvoLF24/QtVEVXrktgmoVwgIdljH5ql+/fvTr1y/QYRR4fkv4qjpfRML9tf3i\nLjlF+eS3bbz5nw2UCi7BW3e25LbWNeys3hiTqYBftBWRgcBAgNq1awc4msJhw74TPPvvlSzdeZRr\nm1Tl5V7NqVre2uqNMVkLeMJX1XHAOHAu2gY4nALt1Nlk3v15Ex/9dyvlQoN5p3ckPSOr21m9McaT\ngCd8483c9Qd4ccZqYo+c4q6omjx7UxMuKlMy0GEZYwoR++FVAbf/+GkenbKEByYsJjQkiM8HduCN\nO1pasjcB8cMPP9CoUSMaNGjAa6+9lmGZefPmUaFCBSIjI4mMjOSll17K0T7eeecdTp48mRfhpjlx\n4kRaPJGRkVSuXJknn3zyvHLR0dFpZVq2bMnXX3+dtuyqq66iUaNGacsPHDhwwXHdfPPNnoaYzjOq\n6pcbMBXYCyQCscCD2a3Tpk0bNY6k5BSd8Ns2bfa3H/Sy52fpez9v1DOJyYEOywTQ2rVrA7r/pKQk\nrVevnm7ZskXPnDmjLVq00DVr1pxXbu7cudqtW7dc76dOnToaFxd3IaFmq3Xr1vrrr7+e93xCQoIm\nJiaqquqePXu0SpUqaY+vvPJKXbx4sV/jyk5G7wEgRj3mZX/20rnHX9su6mK2H+bv361l1e5jdG5Y\nmRG3NqdOJW/jiZhiYvYzsG9V3m7z0gi4KeOzdnDOfhs0aEC9evUAuPvuu5kxYwZNmzbN1e4GDRrE\n4sWLOXXqFHfccQd///vfGTVqFHv27KFr165UrlyZuXPnUrZsWQYNGsSsWbOoVq0ar7zyCkOHDmXn\nzp2888479OjRI0f73bhxIwcOHKBz587nLStdunTa/dOnT+f4+ljqcM6bN2/m4MGDDB06lAEDBrB3\n71569+7N8ePHSUpKYsyYMXTu3Jnw8HBiYmLOG9vHX6wNvwDZc/QUr81ez7cr9lCtQiij7mnFLS2q\n2UVZUyDs3r2bWrVqpT2uWbMmixYtyrDs77//TosWLahRowYjR46kWbPzR/d8+eWXufjii0lOTuaa\na65h5cqVDB48mLfffvucAc4SEhK4+uqrefPNN+nVqxcvvPACc+bMYe3atdx///306NGDDRs2nDes\ncqp58+ZRseL/JviZNm0avXv3zvT/atGiRfTv358dO3YwefJkgoP/lybvv/9+QkJCuP3223nhhRcy\n3MbKlStZuHAhCQkJtGrVim7dujF16lRuuOEGnn/+eZKTk/O8ycorS/gFwKmzyYybv5Uxv25GFQZf\n05BHrqxH6ZL28phMZHEmHmitW7dm586dlC1bllmzZnHrrbeyadOm88p98cUXjBs3jqSkJPbu3cva\ntWtp0aLFeeVKlizJjTfeCEBERASlSpUiJCSEiIgItm/fDjgDnnkdcXLatGlMnjw50+Xt27dnzZo1\nrFu3jvvvv5+bbrqJ0NBQpkyZQo0aNThx4gS33347kydP5k9/+tN56/fs2ZOwsDDCwsLo2rUr0dHR\ntG3blv79+5OYmMitt95KZGSkp1jzml20DSBVZebKPVz79q/886eNXNO4Kj8NuZIh111myd4UODVq\n1GDXrl1pj2NjY6lRo8Z55cqXL0/ZsmUB56JkYmIiBw8ePKfMtm3bGDlyJD///DMrV66kW7duGQ6n\nDBASEpJ2Jl2iRIm0wdJKlCiRNrzyhg0bzrko63vzvSi6YsUKkpKSaNOmTbbH26RJE8qWLcvq1avT\njh+gXLly3HvvvURHR2e4XkbDNHfp0oX58+dTo0YN+vXrx6RJgRmAwBJ+gKzYdZTe4xby58+WUT4s\nhGkDOzC6T2tqXVw6+5WNCYC2bduyadMmtm3bxtmzZ5k2bVqG7ef79u1L7bhBdHQ0KSkpaROTpDp+\n/DhlypShQoUK7N+/n9mzZ6ctK1euXKajZWYm9Qw/o5tvc87UqVO5557MLy9u27Yt7UNkx44drF+/\nnvDwcJKSktI+tBITE5k5cybNm2c8RNiMGTM4ffo0hw4dYt68ebRt25YdO3ZQtWpVBgwYwEMPPZQ2\n21V+s9PIfLb5QDwj/7OBH9bs4+IyJXm5V3PubluboBLWTm8KtuDgYN5//31uuOEGkpOT6d+/f1rb\n/NixYwF45JFHmD59OmPGjCE4OJiwsDCmTZt23llvy5YtadWqFY0bN6ZWrVrnjF0/cOBAbrzxRqpX\nr87cuXPz9Bi++OILZs2adc5z3377LTExMbz00kssWLCA1157jZCQEEqUKMG//vUvKleuTEJCAjfc\ncAOJiYkkJydz7bXXMmDAgAz30aJFC7p27crBgwd58cUXqV69OhMnTuTNN98kJCSEsmXLBuwM34ZH\nzid7jp7i3Z828eWSXYSFBDGgSz0e6lyPsqXsM9d4Y8MjF3zDhw+nbNmyPPXUU37ZfoEdHtk4jiSc\n5V/zNjPxjx2g0O+KujzWtT6VypbKfmVjjMlDlvD95NipRCb9vp1x87cSfzaJ21rV5MlrG1obvTFF\n2PDhwwMdQpYs4eexQ/Fn+Pi3bUz6fQcnziRxbZOq/PWGRjS6tFygQzPGFHOW8PPI/uOnGTd/K58t\n2snppGRuan4pj17VgOY1KgQ6NGOMASzhX7Bdh08y5tctTI+JJVmVni2r82jX+jS4xM7ojTEFiyX8\nXFBVlu48woTfdzBr1V6CRLgjqiaPdKlP7UrWRm+MKZjsh1c5cDoxmS9idtH9vQXcPuYP5m04wANX\nhDN/aFde6RVhyd4UeV6GR54yZQotWrQgIiKCK664ghUrVqQtCw8PJyIigsjISKKiPPUkTDNv3jx+\n//33C4o/I/Pnz6d169YEBwczffr0TMtNnTqViIgIWrRowY033njer4e/+uorRIS86Fo+duxYv/TV\ntzN8D2KPnOTThTv5fPFOjpxM5LKqZXm5V3N6taphQyCYYiM5OZnHHnuMOXPmULNmTdq2bUuPHj3O\nGy2zbt26/Prrr1x00UXMnj2bgQMHnjPImu/AaDkxb948ypYtyxVXXHHBx+Krdu3aTJgwgZEjR2Za\nJikpiSeeeIK1a9dSuXJlhg4dyvvvv5/WK+fEiRO8++67tG/fPk9ieuSRR/JkO+lZtsrEmaRk5m2I\n46slsfy0bj8A1ze9lPuvCKdDvYttBEsTUK9Hv876w+vzdJuNL27M0+2eznS51+GRfRNyhw4diI2N\nzVEc3333HSNGjODs2bNUqlSJKVOmcOrUKcaOHUtQUBCffvop7733HuPHjycsLIxly5Zx4MABPv74\nYyZNmsQff/xB+/btmTBhgqf9hYeHA87YPJlJHU8+ISGBSpUqcfz4cRo0aJC2/MUXX+Tpp5/mzTff\nzHD97du3c+ONN9KmTRuWLl1Ks2bNmDRpEqVLl+aZZ57h22+/JTg4mOuvv56RI0f67QdclvB9pKQo\nMTuO8PWy3cxatZdjpxKpXLYkj1xZnz4d6lCjYligQzQmYHIyPHKq8ePHc9NNN6U9FhGuvfZagoKC\nePjhhxk4cOB563Tq1ImFCxciInz00Ue88cYbvPXWWzzyyCPnJMHx48dz5MgR/vjjD7799lt69OjB\nb7/9xkcffUTbtm1Zvnw5kZGR9O7dmw0bNpy3nyFDhmQ42mVGQkJCGDNmDBEREZQpU4aGDRsyevRo\nAJYuXcquXbvo1q1bpgkfnAHexo8fT8eOHenfvz//+te/eOCBB/j6669Zv349IuL32a8s4QOb9p/g\n62W7mbF8D7uPniIsJIgbm1/Kra1q0LF+JYKD7FKHKViyOhMvKObOncv48eNZsGBB2nMLFiygRo0a\nHDhwgOuuu47GjRvTpUuXc9aLjY2ld+/e7N27l7Nnz1K3bt1M93HLLbcgIkRERFC1alUiIiIAaNas\nGdu3bycyMpLPP//8go8lMTGRMWPGsGzZMurVq8fjjz/Oq6++ynPPPceQIUM8fZvwHTOob9++jBo1\niieffJLQ0FAefPBBunfvTvfu3S841qwUy4R/OjGZRdsOM39jHPM3xrHpQDwlBDo3rMJfb2jEdU2r\nUsbGuDHmHF6HRwZnEpCHHnqI2bNnnzNSZmr5Sy65hF69ehEdHX1ewn/88ccZMmQIPXr0YN68eVn+\netV3qOTU+6mPU0e9zIsz/NSx9uvXrw/AXXfdxWuvvcaJEydYvXo1V111FeCMFNqjRw++/fbb8y5K\nZzRscnBwMNHR0fz8889Mnz6d999/n19++cVTTLnh16wmIjcC7wJBwEeqGpBZG1SVTQfimb8xjl83\nxhG97TBnklIoGVyCduEXc2/72nRvUZ0q5Wx8G2My4zs8co0aNZg2bRqfffbZeeV27tzJbbfdxuTJ\nk7nsssvSnk9ISCAlJYVy5cqRkJDAjz/+yN/+9rfz1j927FjaB8PEiRPTni9XrhzHjx/Pcdx5cYZf\no0YN1q5dS1xcHFWqVGHOnDk0adKEChUqnNNb56qrrmLkyJEZ9kDauXMnf/zxB5dffjmfffYZnTp1\nIj4+npMnT3LzzTfTsWPHtOsj/uK3hC8iQcBo4DqcScwXi8i3qrrWX/vMyJy1+3nxm9XsO+5MrtDg\nkrLc2742XS6rQoe6lQgrGZSf4RhTaHkdHvmll17i0KFDPProo2nrxcTEsH//fnr16gU4vV7uvffe\ntJmsfA0fPpw777yTiy66iKuvvppt27YBTvPNHXfcwYwZM3jvvffy7LgWL15Mr169OHLkCN999x3D\nhg1jzZo1AERGRrJ8+XKqV6/OsGHD6NKlCyEhIdSpU8fzReFUjRo1YvTo0fTv35+mTZsyaNAgjh07\nRs+ePTl9+jSqyttvv51nx5URvw2PLCKXA8NV9Qb38bMAqvpqZuv4Y3jkVbHHGPPrZro0rELny6rY\nhVdTaNnwyIXX9u3b6d69e9rsWblVkIdHrgHs8nkcC5zXSVVEBgIDwekPm9cialbgX32yn87MGGOK\nuoB3P1HVcaoapapRVapUCXQ4xhiT58LDwy/47D4v+DPh7wZq+Tyu6T5njMmlgjRDnclfefHa+zPh\nLwYaikhdESkJ3A1868f9GVOkhYaGcujQIUv6xZCqcujQIUJDQy9oO35rw1fVJBH5M/AfnG6ZH6vq\nGn/tz5iirmbNmsTGxhIXFxfoUEwAhIaGUrNmzQvahl/74avqLGBWtgWNMdkKCQnJ8lenxmQn4Bdt\njTHG5A9L+MYYU0xYwjfGmGLCb7+0zQ0RiQN25HL1ysDBbEvlP4srZyyunLG4cqYoxlVHVT39iKlA\nJfwLISIxXn9enJ8srpyxuHLG4sqZ4h6XNekYY0wxYQnfGGOKiaKU8McFOoBMWFw5Y3HljMWVM8U6\nriLThm+MMSZrRekM3xhjTBYs4RtjTDFRqBK+iNwoIhtEZLOIPJPBchGRUe7ylSLSuoDEdZWIHBOR\n5e7t/Ik8/RPXxyJyQEQyHIg7gPWVXVyBqq9aIjJXRNaKyBoReSKDMvleZx7jyvc6E5FQEYkWkRVu\nXH/PoEwg6stLXAF5j7n7DhKRZSIyM4Nl/q0vVS0UN5wRN7cA9YCSwAqgaboyNwOzAQE6AIsKSFxX\nATMDUGddgNbA6kyW53t9eYwrUPVVDWjt3i8HbCwg7zEvceV7nbl1UNa9HwIsAjoUgPryEldA3mPu\nvocAn2W0f3/XV2E6w28HbFbVrap6FpgG9ExXpicwSR0LgYoiUq0AxBUQqjofOJxFkUDUl5e4AkJV\n96rqUvf+CWAdzlSdvvK9zjzGle/cOoh3H4a4t/S9QAJRX17iCggRqQl0Az7KpIhf66swJfyM5shN\n/6b3UiYQcQFc4X5Fmy0izfwck1eBqC+vAlpfIhIOtMI5O/QV0DrLIi4IQJ25zRPLgQPAHFUtEPXl\nIS4IzHvsHWAokJLJcr/WV2FK+IXZUqC2qrYA3gO+CXA8BV1A60tEygJfAU+q6vH83HdWsokrIHWm\nqsmqGokzhWk7EWmeH/vNjoe48r2+RKQ7cEBVl/h7X5kpTAnfyxy5gZhHN9t9qurx1K+Y6kwKEyIi\nlf0clxcFct7hQNaXiITgJNUpqvrvDIoEpM6yiyvQ7zFVPQrMBW5Mtyig77HM4gpQfXUEeojIdpym\n36tF5NN0ZfxaX4Up4XuZI/db4E/ule4OwDFV3RvouETkUhER9347nHo/5Oe4vAhEfWUrUPXl7nM8\nsE5V386kWL7XmZe4AlFnIlJFRCq698OA64D16YoFor6yjSsQ9aWqz6pqTVUNx8kTv6hq33TF/Fpf\nfp3iMC9pJnPkisgj7vKxONMp3gxsBk4CDxSQuO4ABolIEnAKuFvdS/L+JCJTcXojVBaRWGAYzgWs\ngNWXx7gCUl84Z2D3Aavc9l+A54DaPrEFos68xBWIOqsGTBSRIJyE+YWqzgz0/6THuAL1HjtPftaX\nDa1gjDHFRGFq0jHGGHMBLOEbY0wxYQnfGGOKCUv4xhhTTFjCN8aYYsISvjHGFBOW8I0xppiwhG/8\nTkR+z2H54SLylL/iyan08eT0ePI4ljAR+dX9UVFebrekiMwXkULzY0yTc5bwzQVzfwae6XtJVa/w\n9z7yc3t5cTyZ8ZDI+wP/VtXkvNyvO7T3z0DvvNyuKVgs4ZtcEZFwcWb5mgSsBmqJyDciskScWYYG\n+pSN91lnnYh86Jb50R3rBBF5XkQ2isgCoFEW+xgiIqvd25M+5daLyBR3+9NFpPQFxnxePL7H425r\ntc9zT4nIcPd+GRH5XpwZl1aLSJZJVES+FJEPRGQh8Gw2Vd8HmJHuuCe4sU4RkWtF5DcR2SQi7byU\n8dn2N+72TVGVm1lT7GY3IBxnTO8OPs9d7P4Nw0moldzH8T7rJAGR7uMvgL5AG2AVUBoojzOOyFPp\n9+FTrgxQFliDMzZ8OM4EFx3dch8DT+U25szi8Vkn3t3Wap/nngKGu/dvBz70WVYhm7pcD7zkoc5L\nAvvSHU8SEIFz8rbEPXbBmUjjGy9lfLYXBMQF+r1lN//d7AzfXIgd6szKk2qwiKwAFuIM8dowg3W2\nqWrqAGBLcBJSZ+BrVT2pzjjvvqON+u6jk1suQZ2hbf/trguwS1V/c+9/6pbNbcxZxePFKuA6EXld\nRDqr6rHMCopIKHAx8JKH7VYGjqZ7bpuqrlLVFJwPwJ9VVd0YwnNQBnWaic6KSDkvB2kKH0v45kIk\npN4RkauAa4HLVbUlsAwIzWCdMz73k8l+xNaEbJanSj8KYGajAuYm5owkce7/T9p6qroRZ87eVcAI\nyXqC7GY485YmedjnqQzi863PFJ/HKfyvbr2USVUKOO0hFlMIWcI3eaUCcERVT4pIY5wJmL2aD9zq\n9kApB9ySSbn/uuVKi0gZoJf7HEBtEbncvX8vsOACYvYSz37gEhGpJCKlgO6pC0SkOnBSVT8F3sRJ\n/ojIzyKSfrq6CGCl7xOZlENVjwBB7reCPCcilYCDqproj+2bwLMuWCav/AA8IiLrgA04TSSeqOpS\nEfkcWIEzB+niLMpNAKLdpz5S1WXizPO6AXhMRD4G1gJjchuzl3hUNVFEXnJj2c25E2xEAG+KSAqQ\niDPuegmgAedP3h7hczxkUS7VjzjNVT95OL6c6gp874ftmgLCxsM3hZ6b8GeqaoGYTzUj4syp2l9V\nh1xIORFpDfxFVe/zQ4z/Bp5xm6RMEWQJ3xR6hSHh5yUR6Q9M1Dzsiy/u9JyqOimvtmkKHkv4xhhT\nTIgaU6MAAAAzSURBVNhFW2OMKSYs4RtjTDFhCd8YY4oJS/jGGFNMWMI3xphiwhK+McYUE5bwjTGm\nmPh/jXv897pCCq4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x4db8f957b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# now, make a plot\n",
    "p_atm=14.69 #atmospheric pressure (psi)\n",
    "plt.plot(list_of_r,list_of_p,label=\"Raindrop Pressure\")\n",
    "plt.plot([0,4],[p_atm/2,p_atm/2],label=\"0.5 atm=7.35 psi\")\n",
    "plt.plot([0,4],[p_atm/4,p_atm/4],label=\"0.25 atm=1.84 psi\")\n",
    "plt.ylabel(\"raindrop pressure, (psi)\")\n",
    "plt.xlabel(\"raindrop radius, $r$, (mm)\")\n",
    "plt.title(\"Pressure a raindrop exerts on the ground\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Caveats**\n",
    "1. To find the terminal velocity of the raindrop, I assumed the drag force is quadratic in speed, $F\\approx b v^2$.\n",
    "2. I assumed constant air density of $1.2kg/m^3$.  Air pressure (and density) change with elevation and weather pattern.\n",
    "3. I assumed that the raindrop-ground collision takes $1ms$.\n",
    "4. I assumed that the raindrop collision was spread over the raindrop's cross-sectional area, $A=\\pi r^2$.\n",
    "5. Literature suggests that raindrops are roughly $r=0.2mm$ or smaller, so the horizontal scale of my graph is somewhat exaggerated.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What force does the raindrop exert on the ground?\n",
    " Recall that earlier, we estimated the raindrop-ground force as,\n",
    "\n",
    "$F=\\frac{m~v_T}{\\Delta T}$.\n",
    "\n",
    "Again, remember from above that:\n",
    "\n",
    "$m=\\frac{4}{3} \\pi r^3 d_{w}$ and $v_T = \\left[ \\frac{ 8 r g}{3C}\\frac{ d_{water} }{ d_{air} } \\right]^{0.5}$\n",
    "\n",
    "What does a plot of this look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#redefine constants so that this code can be stand-alone\n",
    "delta_T_ms = 1 # ms\n",
    "delta_T = delta_T_ms*(1.0/1000.0) # seconds\n",
    "#\n",
    "d_w_gpcm3 = 1 # grams per cm^3\n",
    "d_w = d_w_gpcm3*(1.0/1000.0)*(100.0/1.0)**3 # kg/m^3\n",
    "#\n",
    "d_air = 1.2 # kg/m^3$\n",
    "C = 0.45\n",
    "g = 9.8 # m/s^2\n",
    "\n",
    "# plots require a list of points to plot\n",
    "list_of_r=[] # list of radaii, in mm\n",
    "list_of_f=[] # list of forces, in N\n",
    "\n",
    "# this is a for (or for-each) loop.  It repeats the lines below with \n",
    "# i=0, i=1, i=2 ... up to but not including i=30\n",
    "for i in range(30):\n",
    "    # this is the repeated commands tat the forloop runs\n",
    "    \n",
    "    # vary the radius (in mm, incremented by 0.1mm steps)\n",
    "    r_mm=0.1+i*0.1\n",
    "    #print(i,r)\n",
    "    \n",
    "    # save this loop's radius in a list\n",
    "    list_of_r.append(r_mm)\n",
    "    \n",
    "    # convert the radius to (standard) units of m\n",
    "    r = r_mm*(1.0/1000.0) # meters  \n",
    "\n",
    "    # mass, in kg\n",
    "    m=(4.0/3.0)*pi*(r**3)*d_w\n",
    "   \n",
    "    # terminal velocity\n",
    "    v_T=(8.0*r*g*d_w/(3*C*d_air))**0.5\n",
    "    \n",
    "    # force, N\n",
    "    F=m*v_T/delta_T\n",
    "    list_of_f.append(F)\n",
    "    \n",
    "# end of tabbed in section means the loop is over"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEtCAYAAADnbHzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VGX2wPHvIQkkQCCGKgQMVVoKEJpiW0URCxZ+NhCx\noa66li2wrn2LrOvaKyK6VnRXV9EVewGkBqSGXhN6DSEhIeX8/rg3cUgyk0nIZDLJ+TxPnszc+957\nzy0zZ95b3ldUFWOMMQagQbADMMYYU3tYUjDGGFPCkoIxxpgSlhRqORGJFBEVkTgv49uLyBwRyRKR\nv9Z0fDVFRN4QkT9UcdpbReSb6o7JeCciO0VkaLDjCDQRmSYi9wc7jupkSaEKROSwx1+RiBzxeD+6\ngmmHi8j6agzn18BmVY1W1T9V43xrFVUdp6qPBzuO2qS2JLu6+MVYn4UHO4BQpKpNi1+LyGbgJlUN\n1ofzJCCtKhOKSLiqFlRzPCEVR23ZBpUlIvbZLUeo7s/axGoKASAiUSLygojsEJEMEfmHiESISAvg\nv0Bnj5pFCxE5VUTmi8hBEdkuIk/586EXkfeAK4EH3Hmd5m3ZbvnhIrJeRB4QkV3AS+7w/xORZe4p\nqHUicrY7PFZE3nRPBaSLyEMiUu4xU5l1EJEeIlIgIjeLSDrwuYiEi8iHIrLLncf3InKyxzQlv0Y9\n1uM+EdkjIts8a2gi0lpEPheRQyIyFydxFo8rPh13m4hsAFa4w88QkcUikiki80RkgMc080TkzyKy\nyB3/oYg097FfLnW350ERmSUivTzWe7+I9HHfnyQi+0RkSEXb260VfOfu2wPAH4GngTPdfb/TLTdS\nRFa7+zJdRH7jJcYwEXlERLa623yqiESX2j/Xu8fQHhH5vZf5/Aa4nF+OwX97jB4gIivcbfaOiDSs\naBt5WcYF7nF5UESedvfHGC/bZWIF61ampi4ep7pEZJIb63vuNlwmIskeZQeKyFJ33NtAQ+oaVbW/\n4/gDNgPnlBr2ODALaAm0ARYCf3LHDQfWlyo/EBgAhAFdgPXAre64SECBOC/LnwbcX4llFwCP4hzM\nUcDpwAHgLJwfCR2B7m75GcBzQGPgROBn4DovcXhdh3LK9nDXaYo77yicWutYoKm7zi8B88pbT3c9\n8oE/ARHApUAW0NQd/zHwtjvfZGAX8E2p7fk/IMYt0xo4BFzhxjEO2AM0d6eZB2xx424KfApM8bJu\ng4EdQH93W4wH1gLh7vg7gaXucn8A/uIxrdftDdzq7rub3flGucO+KbX8fcBA93ULoK+XOH8NrMJJ\nmM2Az4BXS+2fF9ztNQA4CnT25xh0h+0EfsI5Blu5x8M4f7ZRqfmcCBwGLnT39R/cfT/Gx3bxtW7l\nff52AkPd15OAHGCYO7+ngB/ccVFu3L92YxntLvv+8rZLqP4FPYBQ/6P8pLAN+JXH+5HAavd1mYOy\nnHlOBN5zX1c2KVS07GwgwmP8v4DHypnvSeWUvR6Y4ed2KVmHcsYVf+m08zF9W6AIiCy9nu56ZAIN\nPMofwkkAke508R7jnqRsUjjFY/zNwMxSy/8ZuMp9PQ942GNcPyDbS9yv4yZhj2FbgEHuawG+BJYD\nqcXbt6LtjfPlt7bUfMtLCrvd6aIr2D8/ATd4vE/C+TIUj/3T0mP8MuASf45Bd9hOYJTH+2eBp/3Z\nRqWGjwe+93jfwF1Hz6RQerv4Wjd/ksJnpfb1Qff1ucCmUtMuLr3uof5np4+qmYgIzhfaFo/BW4D2\nPqbpJSIz3KruIeBBnF/6gVj2TlXN93jfAdhQzuxOwvkC3eNW2w8Cz+D88quOdShS1e0e04eLyD9F\nZKM7/WqcD3ELL9PvUdUij/c5OL/i27rTpXuM89wexTzHtyunTOntVnp+jb2cQjoJuK94m7nbrVXx\nvNT5JpkC9AGe8dgX/mxvzxi8GYlzOmere1plgJdypdd5C84v4Vj3faGq7vUYX7x9K2Onl+l9bqNy\n4ixZb3efbytVpvR2qWjdqhp3OyCjVNnyjq2QZkmhmrkf+p14nMfGOSVTfCCX167Iqzi/OLqoajOc\n0zsSgGWXt/x0nNM9paXjVNtPUNUY96+ZqvbzsvjKrkPpOK7HqbKfBTTH+bVKBfMoz0533h08hnWs\nYPnbOXabFU/jud1Kzy9HVTPLmW868KDHNotR1caq+hGAiMQA/wReA/7ikVj82d6lt1mZY0lV56rq\nhTjJ5Cvg3XJiLG+dOwJHgP1eyvtS2bZyfG6jUnYAJbdju9dYSieP0sv3tW7ZOKfniucXgf/J4phY\nPOZdp1hSCIz3gIfEuYjcGufc99vuuF1AaxHx/NUVDWSq6mER6Y1zOiMQyy7PFOAWETldRBqISAcR\n6a6qm3BOmzwuItHuuG7i/d7z412HaCAX55x4E+AvlZweAFXNxTnn/4g4F90Tcc79+jId6Csio9wa\ny1icD/sMjzLjRKS7u98eBt73Mq/JwJ0ikiKOpiJysYgUfxG9gHM65CbgR5xrCFRhe4NzLHWQX24k\naCIiV4lIM5zz7lk4p9LK8x7wOxHp6F6E/QvwrvvDorJ2AZ0rUb6ibeRpOjBIREaIc+PCvcAJFczf\n17qtAmJF5Gx3uz2C/9+DM4FI9+J2uIhcDSQWj5RfbmIY7Of8aiVLCoHxIM5toiuBJTjnOIvvsV+K\nc6BvcavOscA9wE0ichjnS8PbF87xLrsMVZ2Fc172RZzz9N/yy6+hq3Euxq7G+ZX1Pl5OH1XDOryG\nc3F3J8759tmVnN7TLThx7gJewTmH7ZWq7gIuxkmg+4A7gAtL1QTewvmy2YbzRftbL/P6CfiNu9yD\nOBdQrwFURK4AhrrjwbnofIaIXO6+r8z2BvgC55rWbhEpPq1xA84pjUycC/djvUz7EvARMAfn9OF+\nnC/cqpiMc6fRQRGZVlFhX9uonLI7cLbLs8BenGNzOZDnYxFe1809JXYX8A7OqaCd7nwrpKpHcG5q\n+DXOzRkX4PwAKdbBHV6lW8RrC6naDwNj6g8RmQc8r6q+alymBri1hZ3ARao6N9jxeBKRm4D2qvpI\nsGM5HvYAjDGmVhOR83F+9efh1OZygEVBDaocqjol2DFUBzt9ZIyp7U4HNuHcino2cKmqHg1uSHWX\nnT4yxhhTwmoKxhhjSlhSMMYYU8KSgjHGmBKWFIwxAeO2KjpXRGa6LY9GBDsm45slBWNMIKXjNNB4\nOs6DdiODG46piCWFekhEVorImcGOI1SIyGYROaec135tx5rc3rVt36rqDvdJYHCa3/bW7EYZ4nTB\nWqXmTipLRE4WkSXi9JNQbh8U9YUlhUoQkR9E5ICINAp2LP7y/BIrpqq9VfWHIIWEiNwhIqkikici\nbwQrjuPl73YM1PaujfvWGxE5Cafp6U8rKhskf8BpkypaVZ8NdjDBZEnBTyISD5yG0z7LxUENJvRt\nx2mkbGqwAzGB5zbQ9xZOJzv5FZX3Y36BaInhJJz2wiotQPEEjSUF/43FacXyDeA6zxEiMkGcLiGz\nRGSN/NKdpbfhKiJdPaY/pprs/gL8vThdAWaLyGsi0kac/gqyROQbETmhVPk/ikiaW5N53W2x8S2c\n1j4/FaerxD94lC8+BdLTrQEddE89HJPw3LK/c2PJFJH3RSTS1/pVRFU/UtWPcRqfK0NEbhSRL0Xk\nJXd91orTX8Nd4nSxuFdELvMo31lE/ucOPyQiX3uZbwcR+Uic7iX3icjz/mwDb0r/Uvexv/3a3r62\ndTnL9mffVvY4aidOV6N7RGST+HkaRUSaucd0S49hvcTpW6M5Tic8j6jqmgrm01ecLlGzROR9nP4l\nPLfNBBFZBmSL00rpRBHZ4JZPE5FLS82vn4j87I7/t7s9y5yOEpHvcJpsf97dlt39/FyUjsfb8VWl\n7Ro0we7lJ1T+cLoT/DVOF4L5QBt3+Mk4F9Paue/jcfonKHe4+1qBrh7zfoNju2XcjJOA2uC0Hb8b\np6+CvjgflO+Ah0qVX4HTSmMsTsuof/EYV7pnuM3AOThdCq4H7sPpnvNXOM0tn1yq7AKcDkZicZoe\nvrWC9XsReNGPbfoX4I1yhj+LkzDOxukS8UOclj/vcmO+C1jsUX42cLtbNhI4tZx5huG0UPsUTtPc\nkTgtlvq7Dc6p4LWv7eHX9va2rX1sP6/7trLHEc4PxEU4rew2xGkKeyNwnp+fj63A6R7v/4vTMum1\n7r78wf270sv0Dd19fI+7nUbhfM48j+MlOMd4lDvs/9xt1QCnr/Js4MRS8ys+Zi7DuabxFy/L/wG4\nyX3t7zFREg/ej6/j2q7B+LOagh/EadP+JOADVV2E0xzvNe7oQqAR0EtEIlR1s6pu8DHcX8+p6i5V\n3YbT5/J8Vf1Znf4C/ovzwfb0vKqmq+p+4K84zQ1XZDBOr1KTVPWoqn6H059t6WmfVdXt7rw/xen2\n0uv6qeqvVfXXlVjX0pJwugj9VlULcZoiXq6qxT2VreDYxhy74Hwow1Q1V52mmUsbiPMF8ntVzXbL\nza7ENqiIP/vbn2WVt62Ph7/H0QCglao+6sa2EafjpKv8XM4KoCeAiAzC6cbyBVV9S1VbqOqZ7p+3\nJtUH43wZP62q+ar6H5z+xT096x7jRwBU9d/utipy57sOZz8Xzy/cnSZfnQ58Fvi5LpX5XBTH4+34\nOt7tWuMsKfjnOuAr/aV7wnfdYajqeuBunI5XdovINBFp5214JZa5y+P1kXLel+4asXR3kf4sqx2Q\nrsd2a1le16FluieshvXzJRHnQ1isVznvV3u8H41zq+N29xRJeT1pdQC2qGpBqeH+bgOf/Nwe/izL\nW1eQVeXvcXQS0E6O7SLzPnz35+BpBc5+AfgbTp/Wvvo8KK0dsE3dn+Gu0l1dHtPtpoiMFeeOoeJ4\n+/BLF7Dlzc+f7kyLp/XnmPCcn7fj63i3a42zpFABEYkCrsDpDGWniOzEqeImiUgSgKq+q6rFtQkF\n/u5rOM6H3bOXqbbVEGrp7iKL+z/21eLhdpyeuzyPg9LdUHrlY/2qTJy7VBridLxSLBmnql4s0fO9\nqn6nqmfjfCklAePKmXU60FHKXhQ8rm3gyY/tUW3LKl5kFacrTzpOp/SeXWRGq+oIP6dfgVNLOgfn\neH6zksvfAbQXEc/uV0t3dVmyvu5x8ipOh0gtVDXGjaF4+vLm5/kZ8cXf/VQ64ZR3fB3vdq1xlhQq\ndgnOqYFeOF9OyTjV5FnAWHHub/6VOLep5uL8+iryNtyd5xLgGhEJE5HhwBnVEOftIhLn/kr+E7/0\nfOarq8T5OAnqDyISIc797RfhXBj0qYL1q2jacPcCahgQJs5F8eIPUxLOqaIit2wznC/ZZR6zSMI5\nf4uIXCZOt5WC06XnCRybQIotwPmimCROt5WRInLq8WyDUuvkz/aolmV5qGw3mL4sALLci6dR7rHZ\nR0QGQMnNEG/4mL64pvA34E/uab/KmAsUAL9xt81l/HIqqDxNcL6U97jxXY9TU/CcXyFwh3u8jaxg\nfp6qsp+8HV8+t2ttZEmhYtcBr6vqVlXdWfwHPI9z2qIRMAmnS7+dQGvgjz6Gg3Px6yKcrghHAx9X\nQ5zv4nTUvhHnmkfxXRaPAfe7VdffeU6gTpv0FwHnu3G+CIxVVc9TM954XT8ReVlEXvYx7f04X5oT\ngTHu6/vdcUmUrRWsV9Ucd94NcD78xWWG4vR1nAV8jnMe+LvSC3S/pC4CuuJcFM3Aueh5PNvAk6/9\nXRxDdS2rmNd9W1nu9rkQ50fPJje+KUBzt0gHnBsYvEnDqSEUqnNnWWWXfxTnYvA4nO4zr8TpUtNb\n+TTgnzhf/ruABM/4POZ3I87nbAzOKcgKT2lVZT/5OL4q2q61jvWnUAeIyGacOye+CXYspu4RkYY4\nNbNErYbnDIJFROYDL6uqzz676zurKRhjfHLvmukZaglBRM4Qkbbu6aPrcGqdXwQ7rtquTj2JZ4wx\nHk4GPsC5/rARGKWqO4IbUu1np4+MMcaUsNNHxhhjSlhSMMYYU8KSgjHGmBKWFIwxxpQIubuPWrZs\nqfHx8cEOwxhjQsqiRYv2qmqrisqFXFKIj48nNTU12GEYY0xIEZHSDQyWK2Cnj0RkqojsFpEVFZQb\nICIFIjIqULEYY4zxTyCvKbwBDPdVQETCcFqS/CqAcRhjjPFTwJKCqs7EadjKlztxetXaHag4jDHG\n+C9o1xREpD1wKU7fqD6bkRWR8cB4gI4dSzexDvn5+WRkZJCbmxuASI2pXpGRkcTFxRERERHsUIwp\nI5gXmp8GJqhq0bH9YJSlqpOByQApKSll2uXIyMggOjqa+Ph4KpqXMcGkquzbt4+MjAw6deoU7HCM\nKSOYSSEFmOZ+ibcERohIQVXaYs/NzbWEYEKCiNCiRQv27NkT7FCMKVfQkoKqlvxMcnt0+qwqCcFj\nHtURljEBZ8eqqc0CeUvqezi9Ip0sIhkicqOI3CoitwZqmcYYU1c98806Fm6u6N6d4xfIu4+uVtUT\nVTVCVeNU9TVVfVlVy3TTqKrjVPU/gYol0Pbt20dycjLJycm0bduW9u3bl7w/evRomfL79+/n5Zd9\n9VbpKCgoICYmptxxTz75JD179mTs2LHHHb8/du3axYABA+jbty9z5sypkWX6MnToUE4++WSSkpIY\nOHAgy5Ytq3Ca66+/njVr1vi9DF/b35iatHZXFk99s5a5G/YFfFkh90RzbdSiRQuWLHG6DH744Ydp\n2rQpv/ud9y5zi5PCrbdWvdL04osvMnv2bNq2betX+YKCAsLDq767v/76a/r3719uMissLCQsLKzK\n866q999/n+TkZF599VUmTJjAjBkzfJZ//fXq6YXxeLelMZU1eeZGoiLCuHbwSQFfljWIF2CPP/44\nffr0oU+fPjz33HMATJw4kTVr1pCcnMzEiRM5dOgQv/rVr+jXrx+JiYl89tlnPud50003sXXrVoYN\nG8azzz7L3r17ufjii0lMTOSUU05hxQrnIfL777+fsWPHcuqppzJu3DgKCgq455576NOnD4mJibz4\n4osALFy4kDPOOIP+/ftz/vnns2vXrmOWl5qayn333ceHH35YUvuJiYnh7rvvJjExkQULFvDVV1+R\nnJxMQkICN998c0kNKS4ujvvuu4+kpCQGDBjA4sWLOffcc+nSpQuvvvpquet30UUX0b9/f3r37s2U\nKVMq3MZDhgxh27ZtJe/Hjx9PSkoKvXv35tFHHy0ZPnToUJYsWVJSA5g4cSJJSUkMGTKE3budR2U2\nbNjAoEGDSEhI4KGHHiqZ9ptvvuHMM8/kwgsvJCEhweu+Xb9+Pb179+aqq66iZ8+eXHHFFRw5cqTC\ndTDGm52ZuXyyZBtXDujACU0aBnx5de7nziOfriRt+6FqnWevds146KLelZ5u/vz5vPPOOyxcuJCC\nggIGDhzImWeeyaRJk1i/fn1J7SI/P5+PP/6YZs2asXv3bk499VQuvPBCr/OdMmUKX3zxBbNmzSIm\nJobbbruNQYMGMX36dL766ivGjRtX0j7U6tWrmTlzJpGRkTz33HNs376dpUuXEhYWxv79+8nLy+Ou\nu+5i+vTptGzZknfeeYcHHniAyZMnlywvJSWFBx98kBUrVvD0009TUFBAZmYmp59+Ok8//TQ5OTl0\n796dH3/8kS5dujB69GgmT57MHXfcAUCnTp1YunQpd955JzfeeCOzZ8/m8OHDJCUlcfPNN5dZv3/9\n61/ExsaSk5NDSkoKl19+OSeccILX7fHFF19wySWXlLyfNGkSsbGxFBQUcNZZZzFq1Ch69ep1zDSZ\nmZmcccYZTJo0iXvvvZepU6cyceJE7rzzTu666y6uueYannnmmWOmSU1NJS0tjY4dO3rdt1FRUaSl\npfHaa68xePBgxo4dyyuvvMLdd9/tNX5jfHn9p00UFik3Dq2ZW5itphBAs2fP5vLLLycqKoro6Ggu\nueQSZs2aVaacqjJx4kQSExM599xzSU9PZ+/evZVazrXXXgvAueeey/bt28nOzgZg5MiRREZGAs6v\n3VtvvbXkVE9sbCyrVq1i5cqVnHPOOSQnJzNp0iTS09MrXGbDhg259NJLAVi1ahXdu3enS5cuAIwd\nO5aZM2eWlL344osBSEhIYPDgwTRp0oQ2bdrQoEEDDh8+XGbeTz31VMkv+IyMDDZs2FBuDFdeeSXx\n8fH84x//4LbbbisZ/t5779GvXz/69evHqlWrSEtLKzNtVFQU559/PgD9+/dn8+bNAMydO5crr7wS\noGSbFhsyZEjJw5O+9m2nTp0YPHgwAGPGjGH27Nm+NqUxXh3Kzefd+VsZkXAiHWIb18gy61xNoSq/\n6IPtzTffJDMzk8WLFxMeHk5cXFy1PZ3dpEkTn+NVlcTExHKTlS9RUVF+31rZqFEjABo0aFDyuvh9\nQUHBMWW/+eYbZs6cybx584iKimLo0KFet8X7779PUlIS99xzD3fddRcffPAB69at45lnnmHBggXE\nxMQwZsyYcqdv2PCXanhYWNgxcXhbr4q2pbfp7RZUU1Xvzd9KVl4Bt5zepcaWaTWFADrttNP473//\ny5EjRzh8+DCffPIJp512GtHR0WRlZZWUy8zMpHXr1oSHh/P1118fc37c3+W88847gPOl2r59+3K/\nwIYNG8bLL79MYWEh4Fzw7tWrF9u2bWPBggUAHD16lJUrV1Zq+T179mTdunVs3LgRgLfffpszzjij\nUvMolpmZSWxsLFFRUaxcuZKFCxf6LC8i/O1vf2PmzJmsW7eOQ4cOER0dTbNmzdixYwdffvllpZY/\nZMgQPvjgA4CSbVoeb/sWYNOmTSVxv/vuuwwdOrRSMRgDcLSgiKk/beKULi1IiGteY8u1pBBAAwcO\n5Oqrr2bAgAEMHjyY2267jYSEBNq0aUP//v1JSEhg4sSJXHvttcyZM4eEhASmTZtGt27dKrWcRx99\nlLlz55KYmMiDDz7o9S6bW265hbZt25KYmEhSUhIffPABjRo14j//+Q/33nsviYmJ9O3bl/nz51dq\n+Y0bN+a1117jsssuIyEhgUaNGpV7rcAfF1xwATk5OfTq1Yv777+fQYMG+bX8e+65hyeeeIJ+/frR\nq1cvevToUXKRvTKeffZZnnrqKRITE8tccPfkbd+CkySLbxnOyclh/PjxlYrBGIBPlmxj16E8bjmj\n5moJAKJapimhWi0lJUVLd7KzatUqevbsGaSIjPnF+vXrGTVqVMlNBN7YMWt8KSpShj8zkwYizLjr\ntGo5BSkii1Q1paJyVlMwxpha5oe1u1m76zC3nNG5xq9JWVIwphp17dq1wlqCMRV55ceNtGseyYWJ\n7Wp82ZYUjDGmFlmSfpD5m/Zzw9BORITV/Fe0JQVjjKlFJs/cQHRkOFcNLNuhWE2wpGCMMbXE5r3Z\nfLFiJ9cOPommjYLzGJklBWOMqSWmzN5IeIMGjDslPmgxWFIwxphaYN/hPP6dmsGlfdvTullk0OKw\npGCMMbXAv+ZuIa+giJtP7xzUOCwpVKMvvviCk08+ma5duzJp0iSv5QoLC+nbt2+ZllC9Dfc27pVX\nXkFEWLVqVcmwnj17smnTpoDFf8MNN9C6dWv69OlTZty3337LmDFjAPjhhx/KNChX7Le//S1JSUnc\neeedVYrTmLom52gBb83dzDk929C1ddOgxmJJoZoUFhZy++23M2PGDNLS0njvvffKbZ0T4Jlnnin3\naVZvw72NW758OcnJyfzvf/8DIDc3l127dhEfHx+w+MeNG8cXX3xR7jyWLl1KcnJyyeu+ffuWKbNh\nwwZ++uknli5dWtIHgTH13b9TMziQk8+tZwS3lgCWFKrNggUL6Nq1K507d6Zhw4ZcddVVfPLJJ2XK\nZWRk8L///Y+bbrrJr+G+xi1btowJEyaUJIW0tDR69OhRpScg/Y3/9NNPJzY2ttx5LF26lKSkJACW\nLFnCtm3bGDRoEJ07d+aHH35gzZo1nHnmmWzZsoW+ffuWNO9tTH1WUFjElNkb6dcxhpT48j9bNanO\nNZ3NjImwc3n1zrNtApzv/XQQwLZt2+jQoUPJ+7i4uHIblrv77rt5/PHHj2kl1ddwX+PS0tIYOXIk\nf/7zn8nMzGT58uUkJiaWmf60004rd75PPPEE55xzTqXi92Xp0qU88cQTJa8vvvhi5s+fz1dffcUD\nDzzArFmzuO6664iPjy83+RlTH32xcifp+49w/wW9Ki5cAwKWFERkKnAhsFtVy5yAFpHRwARAgCzg\nNlVdGqh4aoPPPvuM1q1b079/f3744YcKh/sal56eTosWLYiKimLYsGF8+eWXLFu2rKSlTk+V7Suh\nKvLz88nMzKRVq1bk5+ezd+9e7rvvPgCSk5NLOg1avnw5I0eODHg8xoQCVeWVHzfSuWUThvVsE+xw\ngMDWFN4Angfe9DJ+E3CGqh4QkfOByUDF7SRXpIJf9IHSvn37Y3osy8jIoH379seU+emnn5g+fTqf\nf/45ubm5HDp0iDFjxtChQ4dyh7/99ttep7nmmmtKEsCIESN455132LFjxzHdUhbzp6bgT/y+eLb6\nuXr1arp27VrSkc3ixYtLTiutXLmy3IvUxtRHs9btZfm2TB67LIEGDWpJZ0yqGrA/IB5Y4Ue5E4Bt\n/syzf//+WlpaWlqZYTUtPz9fO3XqpBs3btS8vDxNTEzUFStWeC3//fff6wUXXOD38NLjHnvsMX34\n4YdVVTU3N1fj4+O1TZs2euDAgYDHv2nTJu3du/cxw958802dOHFiyev4+HjNzc3VrKwsHTx4sM6d\nO1cPHTqkPXr0qFJ8dU1tOGZNcBUVFemlL8zWIX/7RvPyCwO+PCBV/fiOrS0Xmm8EZngbKSLjRSRV\nRFL37NlTg2H5Lzw8nOeff57zzjuPnj17csUVV9C7t9M16IgRI9i+fXu1Lm/58uUlNYVGjRqRmJhI\nw4YNiYmJqdL8/I3/6quvZsiQIaxZs4a4uDhee+01oOydR5dddhmnnHIKAwcO5De/+Q2DBw9mxYoV\nVkswxjVr3V4Wbz3I7b/qSsPw2vJVHOBOdkQkHvhMy7mm4FHmLOBFYKiq7qtontbJjqkL7Jit31SV\ny1+aw86VWB1IAAAgAElEQVTMXH74/Vk1khRCopMdEUkEpgAj/UkIxhhTF9TWWgIEMSmISEfgI+Ba\nVV0brDiMMaYmqSpPf7OWds0j+b/+HSqeoIYF8pbU94AzgZYikgE8BEQAqOrLwINAC+BF92GrAn+q\nNsYYE8qKawl/vbRPraslQACTgqpeXcH4mwB7gskYU2/U9loCWDMXxhhTY2rztYRitTMqY4ypY0Kh\nlgCWFIwxpkaEQi0BLCkYY0zAhUotASwpVCt/Oql56qmn6N27N3369OHqq68mNzfX53Bf42q6k530\n9HTOOussevXqRe/evXnmmWeOGW+d7BhTvlCpJYAlhWrjTyc127Zt49lnnyU1NZUVK1ZQWFjItGnT\nvA73NQ3UfCc74eHh/POf/yQtLY158+bxwgsvHFPGOtkxpqxQqiWAJYVq428nNQUFBRw5coSCggJy\ncnJo166dz+G+xtV0Jzsnnngi/fr1AyA6OpqePXuybdu2kvHWyY4xZYVSLQHqYCc7f1/wd1bvX12t\n8+wR24MJAyf4LONPJzXt27fnd7/7HR07diQqKopzzz2Xc889F8DrcF/TBLOTnc2bN/Pzzz8zaNAv\nrZ1bJzvGHCvUaglgNYUadeDAAT755BM2bdrE9u3byc7O5u233/Y63Nc0le1kZ8mSJWX+ihNCZR0+\nfJjLL7+cp59+mmbNmgGV62SnuDZhTF0XarUEqIM1hYp+0QeKP53UfPPNN3Tq1IlWrVoBcNlllzFn\nzhwaNWpU7vAxY8Z4nSY2NjYonezk5+dz+eWXM3r0aC677LKS4dbJjjHHCsVaAtTBpBAsAwYMYN26\ndWzatIn27dszbdo03n333WPKdOzYkXnz5pGTk0NUVBTffvstKSkpXof7msazVnDGGWdwyy23cOTI\nkSp3x+lP/KrKjTfeSM+ePbn33nuPGVf6esKmTZvIy8sjPz+fRx55hKeeeoqsrCwiIiKIioryf8Ma\nE6JqextH3oROpLWcP53UDBo0iFGjRtGvXz8SEhIoKipi/PjxXocDXscFo5Odn376ibfeeovvvvuO\n5ORkkpOT+fzzzwHrZMcYT6FaS4AAd7ITCNbJjqkL7Jit275fs5vrX1/IXy7pw5jBJwU7HCBEOtkx\nxpi6prBI+fuM1XSMbcwVKaFVSwBLCsYYU60+WpzB6p1Z/GH4ySF1LaFY6EXsRaidBjP1lx2rdVdu\nfiH//GotSXHNuSDhxGCHUyV1IilERkayb98++7CZWk9V2bdvH5GRkcEOxQTA1J82sfNQLn8c0bNK\nLQvUBnXiltS4uDgyMjLYs2dPsEMxpkKRkZHExcUFOwxTzfZnH+Wl7zdwTs/WDO7cItjhVFmdSAoR\nERF06tQp2GEYY+qx575bR/bRAiYM7xHsUI5LnTh9ZIwxwbRlXzZvz9vClQM60K1NdLDDOS4BSwoi\nMlVEdovICi/jRUSeFZH1IrJMRPoFKhZjjAmkf3y5hvAGDbj7nO7BDuW4BbKm8AYw3Mf484Fu7t94\n4KUAxmKMMQGxJP0gny3bwc2ndaJNs9C/gSBgSUFVZwL7fRQZCbypjnlAjIiE5j1cxph6SVV57PNV\ntGjSkPFndAl2ONUimNcU2gPpHu8z3GFliMh4EUkVkVS7w8gYU1t8t3o38zft5+5zutG0UZ24byc0\nLjSr6mRVTVHVlOImpI0xJpgKCouYNGM1nVo24aqBHYMdTrWpMLWJSGvgVKAdcARYAaSqatFxLnsb\n4NkwSJw7zBhjar3/LMpg3e7DvDymHxFhIfH72i9e10REzhKRL4H/4VwUPhHoBdwPLBeRR0Sk2XEs\nezow1r0LaTCQqao7jmN+xhhTI3KOFvDk12vp1zGG83q3DXY41cpXTWEEcLOqbi09QkTCgQuBYcCH\n5U0sIu8BZwItRSQDeAiIAFDVl4HP3WWsB3KA66u8FsYYU4Nem7WJ3Vl5vDi6X8g2Z+GN16Sgqr/3\nMa4A+NjXjFX16grGK3B7RQEaY0xtsvdwHi//uIHzerchJT422OFUO69JQUTG+ppQVd+s/nCMMaZ2\ne/bbdeQWFPGHEG/Owhtfp48GeBl+Mc6to5YUjDH1yqa92bw7fytXD+xAl1ZNgx1OQPg6fXRn8Wtx\nTpqNBiYA84C/Bj40Y4ypPVSVRz9dSWREGHedHfrNWXjj85ZU94LyOOB3OMlglKquqYG4jDGmVvlm\n1W6+X7OH+y/oSavoRsEOJ2B8XVO4HbgL+BYYrqqbayooY4ypTXLzC3nk05V0b9OU606JD3Y4AeWr\npvAcsBsYCpzqcduV4Nw8lBjg2IwxplZ48YcNZBw4wrTxg+vUg2rl8ZUUrNcaY0y9t2VfNi//uIGR\nye1Cukc1f/lKClu1gk6PRUQqKmOMMaFKVXl4+koiGgj3jegZ7HBqhK960PcicqeIHNPSk4g0FJFf\nici/gOsCG54xxgRP8cXle4Z1rxN9JfjDV01hOHAD8J6IdAIOApFAGPAV8LSq/hz4EKvH3xf8ndX7\nVwc7DGNMiChSZWlGJi26wuzDMfz0RbAjgh6xPZgwcEJAl+HrOYVc4EXgRRGJAFoCR1T1YEAjMsaY\nWmD7wVzy8gvpdWIz6ljzRj751SuEquYDId2CaaCzqzGm7tiyL5thT81keO+2PHtp32CHU6Pq9r1V\nxhhTBY98mkZEA+FPF9SPi8ueLCkYY4yHb9J28d3q3dx9Tv25uOzJkoIxxrhy8wt5+NOVdGvdlHGn\nxgc7nKCodFIQkW9EZIaIXBiIgIwxJlhecp9cfnRknzr/5LI3fl1oLmUsTtecg6s5FmOMCZot+7J5\n6ccNXJzUjiFd6v6Ty974lRREJAroqKprVHU7sB1YFNDIjDGmBj1ajy8ue6qwfiQiFwFLgC/c98ki\nMj3QgRljTE35Om0X39bji8ue/Dlp9jAwEOeJZlR1CdZYnjGmjjiUm88DH6/g5DbR9fbisid/kkK+\nqmaWGuZXI3giMlxE1ojIehGZWM745iLyqYgsFZGVInK9P/M1xpjq8tfPVrE7K5fHRyXW24vLnvzZ\nAitF5BogTES6ichzwJyKJhKRMOAF4HygF3C1iPQqVex2IE1Vk4AzgX+KSMPKrIAxxlTVzLV7eD81\nnfGndyGpQ0yww6kV/EkKdwK9gTzgXSATuNuP6QYC61V1o6oeBaYBI0uVUSDa7QO6KbAfKPAzdmOM\nqbLDeQX88aPldGnVhLvP6RbscGqNCu8+UtUc4E/uX2W0B9I93mcAg0qVeR6YjnM3UzRwpaoWlZ6R\niIwHxgN07Nix9GhjjKm0xz5fxfbMI/zn1lOIjAgLdji1hj93H30tIjEe708QkS+rafnn4dzZ1A5I\nBp4XkWalC6nqZFVNUdWUVq1aVdOijTH11Zz1e3ln/lZuPLUT/U86Idjh1Cr+nD5q6dlctqoeAFr7\nMd02oIPH+zh3mKfrgY/UsR7YBPTwY97GGFMl2XkFTPhoGfEtGvPbc08Odji1jj9Jociz9zUROQn/\n7j5aCHQTkU7uxeOrcE4VedoKnO3Otw1wMrDRn8CNMaYq/vHlGjIOHOHxUUlENbTTRqX580Tzn4DZ\nIvIjIMBpuOf3fVHVAhG5A/gSp7e2qaq6UkRudce/DPwZeENElrvznqCqe6u2KsYY49uCTft5Y85m\nxp0Sz8BOscEOp1YSVe8/+t27guKAI/zS1tG8YH5xp6SkaGpqarAWb4wJUUeOFnL+MzMpVOXLu0+n\nccOqNP0WukRkkaqmVFTO51ZRVRWRz1U1Afis2qIzxpga9sRXa9i8L4d3bx5U7xJCZfhzTWGxiAwI\neCTGGBMgi7bsZ+pPmxg9qCOndGkZ7HBqNX/S5SBgtIhsAbJxzv2rqiYGNDJjjKkGufmF/P4/y2jX\nPIo/jqjfLaD6w5+kcF7AozDGmAB56pu1bNyTzZs3DKRpIzttVJEKTx+p6hYgBrjI/YtxhxljTK22\nJP0gr87cyJUpHTi9uz346g9/nmi+C3gH54G11sDbInJnoAMzxpjjkXO0gHs/WELr6Ej+dKGdNvKX\nP3WpG4FBqpoNICJ/B+YCzwUyMGOMOR4PT1/Jpr3ZvHPTIJpFRgQ7nJDhz91HAhR6vC90hxljTK00\nfel2PkjN4PYzu9rdRpXkT03hdWC+iPzXfX8J8FrgQjLGmKrbui+H+z5aTr+OMdxlTWJXmtekICKd\nVHWTqj4pIj8AQ91R16vqzzUSnTHGVEJ+YRF3TvsZEXjmqr7Wk1oV+Kop/AfoLyLfqurZwOIaiskY\nY6rkia/WsDT9IC+O7keH2MbBDick+UoKDUTkPqC7iNxbeqSqPhm4sIwxpnJmrt3DKz9u5OqBHRmR\ncGKwwwlZvupWV+FcVA7H6RWt9J8xxtQKe7LyuPeDpXRv05QHLyzdFbypDK81BVVdA/xdRJap6owa\njMkYY/xWVKT89t9LycrN552bBlkfCcfJnyeaLSEYY2qtKbM3MnPtHu6/sBcnt7WTGMfLLs0bY0LW\n0vSDPP7FGs7r3YYxgzpWPIGpkCUFY0xIysrN5873fqZ1dCP+fnkiTp9g5nhV+PCaiEQCv8Z5TkGB\n2cBLqpob4NiMMaZcqsr9H68g40AO798yhJjGDYMdUp3hzxPNbwJZ/NLW0TXAW8D/BSooY4zx5cPF\n2/hkyXbuHdadAfHW13J18icp9FFVz3u8vheRtEAFZIwxvqzcnsn9Hy9nUKdYbj+ra7DDqXP87Y5z\ncPEbERkEpPozcxEZLiJrRGS9iEz0UuZMEVkiIitF5Ef/wjbG1Ef7Ducx/s1FnNC4Ic9f04+wBnYd\nobr5U1PoD8wRka3u+47AGhFZjo9uOUUkDHgBGAZkAAtFZLqqpnmUiQFeBIar6lYRaX0c62KMqcPy\nC4v49TuL2Xs4j3/fOoRW0Y2CHVKd5E9SGF7FeQ8E1qvqRgARmQaMBDxPPV0DfKSqWwFUdXcVl2WM\nqeP+8lka8zft56krk0iMiwl2OHVWlbvjLP7zMWl7IN3jfYY7zFN34AQR+UFEFonI2PJmJCLjRSRV\nRFL37NlTUcjGmDrm/YVb+dfcLdx8Wicu7RsX7HDqtGB3xxmOc3rqAuA84AER6V66kKpOVtUUVU1p\n1cr6WTWmPlm05QD3f7yC07q1ZMLwHsEOp84LZHec24AOHu/j3GGeMoB97ryzRWQmkASs9SMuY0wd\ntzMzl1vfXsSJzaN47uq+hFv/CAEXyO44FwLdRKSTiDTEaXV1eqkynwBDRSRcRBoDg4BVfszbGFPH\n5eYXcsvbi8jJK2DKdSn2gFoNCVh3nKpaICJ3AF8CYcBUVV0pIre6419W1VUi8gWwDCgCpqjqiqqs\niDGm7lBV/vTfFSxNP8gr1/anextr6K6miKpWXEikH790xzkrmN1xpqSkaGqqX49JGGNC1NTZm3j0\nszTuOrsb9wwrc5nRVIGILFLVlIrK+awpuM8arFTVHlh3nMaYGvDT+r389fNVnNurDXed3S3Y4dQ7\nPq8pqGohzoNq1iatMSbg0vfncPu7i+nSqglPXplMA3tiucb5c03hBGCliCwAsosHqurFAYvKGFPv\nHMrN5+Y3UykqUiZfm0LTRv58PZnq5s9WfyDgURhj6rXc/ELGv5nK+t2HeeP6gcS3bBLskOqtCpOC\nqlojdcaYgCksUu79YAnzNu7n6SuTGdqtZbBDqte8JgURycLpVKdcqtosIBEZY+oNVeWRT1fy+fKd\n3H9BTy7pW7olHFPTvCYFVY0GEJE/AztwOtYRYDRwYo1EZ4yp0174fj1vzt3C+NM7c9NpnYMdjsG/\nJ5ovVtUXVTVLVQ+p6ks4rZ0aY0yVvb9wK098tZZL+7ZnorVpVGv4kxSyRWS0iISJSAMRGY3HXUjG\nGFNZX6ft4o8fLef07q14fFSi3Xpai/iTFK4BrgB2uX//5w4zxphKS928nzveXUxC++a8NLofEdbI\nXa3iz91Hm7HTRcaYarB2VxY3/iuVdjFRTB03gCb2LEKtU+EeEZFWwM1AvGd5Vb0hcGEZY+qa7QeP\ncN3UBTQMb8CbNwykRVPrTrM28idNfwLMAr7h2Ca0jTHGLwdzjnLd1AVk5Rbw/i2D6RDbONghGS/8\nSQqNVXVCwCMxxtRJOUcLuOlfqWzZl8MbNwygd7vmwQ7J+ODPFZ7PRGREwCMxxtQ52XkFjHt9IYu3\nHuCpK5M5pYs9rVzb+ZMU7sJJDEdE5JCIZInIoUAHZowJbYfzChj3+gJSN+/nqSuTuSDRnnkNBf7c\nfVQ3ujyaMRF2Lg92FMbUCwVFRWzemcVv8wro1r4pLX5uBEHrmqsOaZsA508K6CJ8tX3UQ1VXu72u\nlaGq1umOMaaMgqIiVu/MIjuvgG6tm9Kiid1lFEp81RTuBcYD/yxnnAK/CkhEgRLg7GqMgcwj+Yyd\nuoCV2Zk8f00/WvRpG+yQTCX5ahBvvPv/rJoLxxgTqjJz8rl26nxW7TjEi6P7cW5vSwihyK/ny0Wk\nj4hcISJji//8nG64iKwRkfUiMtFHuQEiUiAio/wN3BhTexzMOcro1+axekcWL43ubwkhhPnzRPND\nwJlAL+Bz4HxgNvBmBdOFAS8Aw4AMYKGITFfVtHLK/R34qgrxG2OC7ED2UUZPmc/63Yd55dr+nNWj\ndbBDMsfBn5rCKOBsYKeqXg8kAf48fTIQWK+qG1X1KDCN8ttQuhP4ENjtX8jGmNpif/ZRrpkyn/V7\nDjN5rCWEusCfpHBEVYuAAhFphvPl3cGP6doD6R7vM9xhJUSkPXAp8JKvGYnIeBFJFZHUPXv2+LFo\nY0yg7cnK45pX57Fxz2GmjE3hzJMtIdQF/iSFVBGJAV4FFgGLgbnVtPyngQlu0vFKVSeraoqqprRq\n1aqaFm2Mqar1uw9z6Ys/sWVfDq9dN4DTu9vnsq7weU1BRAR4TFUPAi+LyBdAM1Vd5se8t3FsjSLO\nHeYpBZjmLIaWwAgRKVDVj/1dAWNMzVq4eT83/SuViDBh2vjBJHWICXZIphr5TAqqqiLyOZDgvt9c\niXkvBLqJSCecZHAVpTrnUdVOxa9F5A3gM0sIxtRe/1u2g3s+WEJcTBRvXD+Qji2stdO6xp/TR4tF\nZEBlZ6yqBcAdwJfAKuADVV0pIreKyK2VnZ8xJnhUlVdnbuT2dxeT2L45H952iiWEOkpU1XcBkdVA\nV2ALTt/MglOJSAx8eGWlpKRoampqMBZtTL1UWKT8+bM03pizmREJbXnyimQiI8KCHZapJBFZpKop\nFZXzpz+F86ohHmNMCDpytJC73/+ZL1fu4qahnbhvRE8aNJBgh2UCyJ9WUrfURCDGmNpl3+E8bnoz\nlSXpB3nool5cf2qniicyIc96zTbGlLF5bzbjXl/AjsxcXhrdj+F9rC+E+sKSgjHmGAs37+eWtxah\nqrx782D6n3RCsEMyNciSgjEGcO4wem32Jh6bsZqOsY157boUOrdqGuywTA2zpGCM4XBeARP+s4z/\nLd/Bub3a8MQVSTSLjAh2WCYILCkYU8+t25XFLW8vYvPebP54fg/Gn94Zt5UBUw9ZUjCmHpu+dDsT\nP1xG44ZhvHPTYIZ0aRHskEyQWVIwph46WlDE3z5fxRtzNpNy0gm8MLofbZpFBjssUwtYUjCmntmR\neYTb31nM4q0HueHUTvxxRA8iwvzqhNHUA5YUjKlH5qzfy53v/UxufiHPX9OXCxPbBTskU8tYUjCm\nHsgvLOKF79fz7Lfr6NyqKS+P6UfX1tHBDsvUQpYUjKnj1u7K4rcfLGX5tkwu7duev1zShyaN7KNv\nymdHhjF1VGGRMmXWRv751VqaRobz0uh+nJ9gzVUY3ywpGFMHbd6bzW//vZRFWw5wXu82/PXSBFo2\nbRTssEwIsKRgTB1SVKS8PX8Lj32+mogw4ekrkxmZ3M4eRjN+s6RgTB2x7eAR/vCfpfy0fh9ndG/F\n3y9PpG1ze/bAVI4lBWNCnKry79QMHv0sDVXlscsSuGpAB6sdmCqxpGBMCNu45zAPf5rGzLV7GNw5\nln+MSqJDrPWdbKrOkoIxISg7r4DnvlvPa7M3EhkexsMX9WLskHjrKtMct4AmBREZDjwDhAFTVHVS\nqfGjgQmAAFnAbaq6NJAxGRPKVJXpS7fzt89XsetQHqP6xzFheA9aRdudRaZ6BCwpiEgY8AIwDMgA\nForIdFVN8yi2CThDVQ+IyPnAZGBQoGIyJpSt2nGIh6avZMGm/SS0b85LY/rTr6P1imaqVyBrCgOB\n9aq6EUBEpgEjgZKkoKpzPMrPA+ICGI8xISkzJ58nv17DW/O20DwqgscuS+CKlA6E2akiEwCBTArt\ngXSP9xn4rgXcCMwob4SIjAfGA3Ts2LG64jOmVissUj5ITecfX67hYM5Rxgw+iXuHdSemccNgh2bq\nsFpxoVlEzsJJCkPLG6+qk3FOLZGSkqI1GJoxNU5V+X7Nbv751VpWbj/EgPgTeOTiQfRq1yzYoZl6\nIJBJYRvQweN9nDvsGCKSCEwBzlfVfQGMx5haTVWZtW4vT369liXpB+kQG8UzVyVzcZI9kWxqTiCT\nwkKgm4h0wkkGVwHXeBYQkY7AR8C1qro2gLEYU6vN2bCXp75ey8LNB2jXPJLHLktgVP846/zG1LiA\nJQVVLRCRO4AvcW5JnaqqK0XkVnf8y8CDQAvgRfeXUIGqpgQqJmNqm4Wb9/PkV2uZu3EfbZo14s8j\ne3PFgA40Cg8LdmimnhLV0DpFn5KSoqmpqcEOw5jj8vPWAzz59VpmrdtLy6aN+PWZXbhmUEciIywZ\nmMAQkUX+/OiuFReajakPVJU5G/YxZdZGvl+zh9gmDblvRA+uHRxPVENLBqZ2sKRgTIDl5hfyyZJt\nTJ29mTW7smjRpCG/P+9krjslnqbWA5qpZeyINCZAdh3K5a25W3h3wVb2Zx+l54nN+MeoRC5Kamen\niUytZUnBmGq2NP0gr/+0ic+W7aBQlWE923D9qZ0Y3DnWbi01tZ4lBWOqQV5BIV+n7eL1nzazaMsB\nmjYKZ+yQeMadEk/HFtaUtQkdlhSMqSJVZWlGJh8uyuDTZds5mJPPSS0a89BFvRjVP47oyIhgh2hM\npVlSMKaSth88wn9/3sZHizPYsCebRuENOK93Wy7r157TurWyhupMSLOkYIwfco4W8MWKnXy4OIM5\nG/ahCgPjY7n5tM6MSDyRZlYrMHWEJQVjvMjNL2T2ur3MWLGTGSt2kHO0kI6xjbnr7G5c1jfOrhWY\nOsmSgjEe9mcf5bvVu/lq5U5mrdvLkfxCohuFc1FiOy7vH8eA+BPsDiJTp1lSMPXeln3ZfJ22i6/S\ndpG6eT9FCm2bRTKqfxzDerVhcOcWNAy3hulM/WBJwdQ7eQWF/Lz1ILPW7eHrtF2s3XUYgB5to7n9\nrK4M69WGhPbNrUZg6iVLCqbOyy8sYllGJnM37GXuxn2kbj5AXkERYQ2EAfEn8MCFvTi3Vxs6xNo1\nAmMsKZg6p7BIWbk9k7kb9jFnwz5SN+8n+2gh4NQGrhnUkSGdWzCoUwuaN7a7hozxZEnBhLzdWbks\nS89kacZBlmZk8vPWA2TlFgDQpVUTLusXx5AuLRjUKZYWTRsFOVpjajdLCiakZB7JZ3mGkwCWZRxk\nWUYmOzJzAWgg0L1NNBcknMiQLi0Y0rkFrZtFBjliY0KLJQVTKx0tKGLLvmzW7T7Mul2HWbs7i7Tt\nh9i0N7ukTHyLxgyIjyUxrjlJHWLo3a4ZjRvaIW3M8bBPkAmqvIJCNu3NZt2uw6zbfZj1u7NYu+sw\nm/dmU1Dk9AooAh1jG3Nym2hG9Y8jMa45ie1j7HqAMQFgScEEVGGRsiPzCOn7j5B+IIeM/TmkHzhC\n+v4cMg4cYVdWLsU9wjYQOKlFE7q2bsp5vdvQrXU0XVs3pUurptYzmTE1xJKCqRJVJSuvgN2H8tid\nlcuerLxjX2flkXHgCNsPHin5xQ/OF/+JzaOIOyGKU7u2pENsFJ1aNqF7m2g6tWxinc8YE2QBTQoi\nMhx4BggDpqjqpFLjxR0/AsgBxqnq4kDGZMoqKlKyjxZwMCefzCP5HMzJ50DOUQ4eyScz56j7Pp/M\nI0c5kJPvfunnkptfVGZejcIb0LpZI1pHR5LcIYaLkk4k7oTGdDihMR1iozixeZQ9HWxMLRawpCAi\nYcALwDAgA1goItNVNc2j2PlAN/dvEPCS+79eUlUKipSCQiW/qIiCQqWgsIj8Ivd/oZJf6Aw/WlhE\nXkEhefnu/4IicvM9/ucXkVdQxJH8QrLzCsg+6vw/nFdAdl4BOUcLj3ntS+OGYcRERRDTuCExjSNI\n7hBD6+hGJV/+xa9bRUfSLDLcngQ2JoQFsqYwEFivqhsBRGQaMBLwTAojgTdVVYF5IhIjIieq6o4A\nxuXT+t1Z3Pq2U1lR/eW0h5Z54bxUVfc/KEqR++O59PDCImdYoSqFRUpRkVKkUKjO60JVPBZXLRoI\nREaE0aRROE0bhdOkURhNGobTplmkO8x5Xzy+eeOIki//ExpH0LxxBM2jImgUbqd0jKkvApkU2gPp\nHu8zKFsLKK9Me+CYpCAi44HxAB07dqz2QD01Cg/j5DbRHgsv+9Lzl7Dg3B0jQAP3hSAlw8R936CB\nENbAKdNAhLAGzp8IhLnvi4dHhDUgIkwIbyCEl7xuQHiY0DCsAeFhzuvI8DAiIxrQqPh/RBiNwhsQ\n6f4PbyD2q90YUykhcaFZVScDkwFSUlKq+ff0sTrENuaF0f0CuQhjjKm1AnnFbxvQweN9nDussmWM\nMcbUkEAmhYVANxHpJCINgauA6aXKTAfGimMwkBnM6wnGGFPfBez0kaoWiMgdwJc4t6ROVdWVInKr\nO/5l4HOc21HX49ySen2g4jHGGFOxgF5TUNXPcb74PYe97PFagdsDGYMxxhj/2VNExhhjSlhSMMYY\nU8KSgjHGmBKWFIwxxpQQre62FQJMRPYAW6o4eUtgbzWGE0y2LrVTXVmXurIeYOtS7CRVbVVRoZBL\nCnZDJ04AAAc9SURBVMdDRFJVNSXYcVQHW5faqa6sS11ZD7B1qSw7fWSMMaaEJQVjjDEl6ltSmBzs\nAKqRrUvtVFfWpa6sB9i6VEq9uqZgjDHGt/pWUzDGGOODJQVjjDEl6mRSEJHhIrJGRNaLyMRyxouI\nPOuOXyYitbZXHT/W5UwRyRSRJe7fg8GIsyIiMlVEdovICi/jQ2mfVLQuobJPOojI9yKSJiIrReSu\ncsqExH7xc11CZb9EisgCEVnqrssj5ZQJ3H5R1Tr1h9NM9wagM9AQWAr0KlVmBDADp8fMwcD8YMd9\nHOtyJvBZsGP1Y11OB/oBK7yMD4l94ue6hMo+ORHo576OBtaG8GfFn3UJlf0iQFP3dQQwHxhcU/ul\nLtYUBgLrVXWjqh4FpgEjS5UZCbypjnlAjIicWNOB+sGfdQkJqjoT2O+jSKjsE3/WJSSo6g5VXey+\nzgJW4fSR7ikk9ouf6xIS3G192H0b4f6VviMoYPulLiaF9kC6x/sMyh4c/pSpDfyN8xS3CjlDRHrX\nTGjVLlT2ib9Cap+ISDzQF+dXqaeQ2y8+1gVCZL+ISJiILAF2A1+rao3tl4B2smNqxGKgo6oeFpER\nwMdAtyDHVN+F1D4RkabAh8Ddqnoo2PEcjwrWJWT2i6oWAskiEgP8V0T6qGq517CqW12sKWwDOni8\nj3OHVbZMbVBhnKp6qLiqqU5PdxEi0rLmQqw2obJPKhRK+0REInC+RN9R1Y/KKRIy+6WidQml/VJM\nVQ8C3wPDS40K2H6pi0lhIdBNRDqJSEPgKmB6qTLTgbHuFfzBQKaq7qjpQP1Q4bqISFsREff1QJx9\nuq/GIz1+obJPKhQq+8SN8TVglao+6aVYSOwXf9YlhPZLK7eGgIhEAcOA1aWKBWy/1LnTR6paICJ3\nAF/i3L0zVVVXisit7viXcfqNHgGsB3KA64MVry9+rsso4DYRKQCO8P/t3c1rXFUYx/Hvz0JfEmuh\nKYKCJQtBwQZFQai1oFB3CooLRe2mCxEE0ZKFLwtDcCP9A0QqpS0tYpUaoQUpdmFtpSix1ET74qKI\nCopiFZL6kpjHxTlzvR0zyU0yM82kvw9ccufOc899zp3kPsnNzDnweOS3Jywmkt4mvftjnaTvgVdJ\n/0DrqNcEKvWlI14TYBOwFRjJ968BXgbWQ8e9LlX60imvyw3AHknLSIXrQEQcatc1zMNcmJlZYSne\nPjIzs3lyUTAzs4KLgpmZFVwUzMys4KJgZmYFFwUzMyu4KJiZWcFFwRYNSZ/OMX5AUn+r8pmr+nzm\n2p8m57JK0sf5A1DNbHe5pGOSltwHXy1xUbC2yR/Jb/g9FxH3tPoY7WyvGf1ppMLFfhtwMA+s1jR5\nCPejwGPNbNcWDxcFaylJvUozx+0FRoGbJA1JGlaaVerpUuxYaZ8zknbmmCN5DBgkvSLpvKTjwC0z\nHGO7pNG8PF+KOytpf27/PUldC8z5f/mU+5PbGi1t65c0kNe7JR1WmmFrVNKMF1pJ70p6U9JJ4KVZ\nTv2TwAd1/d6dc90vaYukE5K+kXR3lZhS20O5fVuKmj1rkBcv5QXoBaYozRwFrM1fV5Euuj358Vhp\nn0ngjvz4APAUcBcwAnQB15HGfemvP0Yprhu4FviKNL5+L2mykk05bhfQP9+cG+VT2mcstzVa2tYP\nDOT1R4GdpefWzHIuzwKDFc75cuDHuv5MAn2kXwSHc99FmqxlqEpMqb1lwM9X+nvLS2sW/6Vg7fBt\npNmhap6TdBo4SRr+d7ox7S9ERG1gs2HSRWsz8H5EXIo0Vn55xNjyMe7NceORhko+mPcF+C4iTuT1\nfTl2vjnPlE8VI8ADkl6XtDkifm8UKGklsBYYrNDuOuC3um0XImIkIqZIRfJoRETOoXcOMUS6JfW3\npNVVOmmdxUXB2mG8tiLpPmALsDEibgdOASun2eev0vo/zD6i7/gsz9fUjwDZaETI+eQ8nUku/zkr\n9ouI86S5nkeA1zTzRPK3kebhnaxwzD+mya98PqdKj6f479xWialZAfxZIRfrMC4K1m5rgIsRcUnS\nraRJx6s6Bjyc31mzGnioQdwnOa5LUjfwSN4GsF7Sxrz+BHB8ATlXyecn4HpJPZJWAA/WnpB0I3Ap\nIvYBO0gFAklHJdVPrdgHfFne0CCOiLgILMt/XTSdpB7gl4iYaEX7dmX5bWXWbh8Cz0g6A5wj3Y6p\nJCK+kPQOcJo0d+3nM8TtBj7Lm96KiFNKc/eeA56VtAv4GnhjvjlXySciJiQN5lx+4PLJUvqAHZKm\ngAnSWP/XADcDv9Y11VfqDzPE1Rwh3Rr7qEL/5up+4HAL2rVFwPMp2FUjF4VDEbHhCqfSkKQNwLaI\n2L6QOEl3Ai9ExNYW5HgQeDHf/rIlxkXBrhqdUBSaSdI2YE808bMKytPCRsTeZrVpi4uLgpmZFfyP\nZjMzK7gomJlZwUXBzMwKLgpmZlZwUTAzs4KLgpmZFVwUzMys8C9bIJBLaLXEdgAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x4db7b6f5c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(list_of_r,list_of_f,label=\"Total force from a Raindrop\")\n",
    "plt.plot([0,3],[0.444,0.444],label=\"$0.444N = 0.1lbf$\")\n",
    "plt.plot([0,3],[0.888,0.888],label=\"$0.888N = 0.2lbf$\")\n",
    "plt.ylabel(\"raindrop force, (N)\")\n",
    "plt.xlabel(\"raindrop radius, $r$, (mm)\")\n",
    "plt.title(\"Total force a raindrop exerts on the ground, \\nAssumptions: $1ms$ collision time, $v^2$ drag force\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}