A jupyter notebook which demonstrates the use of f2py for computing radial distribution functions

f2pyExample.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Radial Distribution Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [radial distribution function](https://en.wikipedia.org/wiki/Radial_distribution_function) (RDF) or the pair correlation function $g(r)$ is a statistical mechanical property which quantifies density variation as a function of distance from a reference particle as a function in a system of particles.\n",
    "\n",
    "It is a fairly fundamental property that can easily be computed from molecular simulations. The Fourier transform of the RDF is experimentally accesible via light scattering.\n",
    "\n",
    "For the purposes of this demonstration, we are given a bunch of points in a periodic simulation box. Our goal is to find the pair correlation function $g(r)$.\n",
    "\n",
    "Let's generate some random points in a cubic box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(12311)\n",
    "npart = 100  # number of particles\n",
    "Lbox  = 10. # box size\n",
    "r     = np.random.uniform(0., Lbox, (npart,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basic idea is quite simple:\n",
    "\n",
    "* consider every particle pair\n",
    "* compute the distance $d_{ij}$\n",
    "* histogram the set of distances.\n",
    "\n",
    "The basic algorithm is goes as the square of the number of particles, and for the most general case, there is no way around it. There are a couple of nuances to do with (i) [periodic boundary conditions](https://en.wikipedia.org/wiki/Periodic_boundary_conditions) (nearest image convention) and (ii) \"edge\" effects, which have to do with trying to fit a round hole in a square peg (simulation box).\n",
    "\n",
    "Let us write some python code to compute the RDF. Luckily we there is code on github by [Craig Finch](https://github.com/cfinch/Shocksolution_Examples/tree/master/PairCorrelation) that we can use without much change:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def py_rdf(r, S, dr):\n",
    "\t\"\"\"5/12/2015: Based on code by Dr. Craig Finch @\t\n",
    "\thttps://github.com/cfinch/Shocksolution_Examples/tree/master/PairCorrelation\n",
    "\t\n",
    "\tModifying it to use (i) all particles, and \n",
    "\t(ii) assume periodic boundaries in all 3 directions.\n",
    "\t\n",
    "\tCompute the three-dimensional pair correlation function for a set of\n",
    "\tspherical particles contained in a cube with side length S.  This simple\n",
    "\tfunction finds reference particles such that a sphere of radius rMax drawn\n",
    "\taround the particle will fit entirely within the cube, eliminating the need\n",
    "\tto compensate for edge effects.  If no such particles exist, an error is\n",
    "\treturned.  Try a smaller rMax...or write some code to handle edge effects! ;)\n",
    "\n",
    "\tArguments:\n",
    "\t\tr               an array of r = [x, y, z] positions of centers of particles\n",
    "\t\t\t\t\t\tr is of size r[nparticles, 3]\t\t\t\t\t\t\n",
    "\t\tS               length of each side of the cube in space\n",
    "\t\tdr              increment for increasing radius of spherical shell\n",
    "\n",
    "\tReturns a tuple: (g, radii, interior_indices)\n",
    "\t\tg(r)            a numpy array containing the correlation function g(r)\n",
    "\t\tradii           a numpy array containing the radii of the\n",
    "\t\t\t\t\t\tspherical shells used to compute g(r)\n",
    "\t\"\"\"\n",
    "    \n",
    "\tfrom numpy import zeros, sqrt, where, pi, mean, arange, histogram, absolute\n",
    "\t\n",
    "\tnum_particles  = len(r)\n",
    "\trMax\t\t   = S/2.0;\n",
    "\tedges          = arange(0., rMax + dr, dr)\n",
    "\tnum_increments = len(edges) - 1\n",
    "\tg              = zeros(num_increments)\n",
    "\tradii          = zeros(num_increments)\n",
    "\tnumberDensity  = len(r) / S**3\n",
    "\n",
    "    # Compute pairwise correlation for each particle\n",
    "\tfor index in range(num_particles):\n",
    "\n",
    "\t\td = 0.0\n",
    "\t\t\n",
    "\t\tfor i in range(3):\n",
    "\t\t\tdp = absolute(r[index,i] - r[:,i])\n",
    "\t\t\tmask = dp>S/2.0\n",
    "\t\t\tdp[mask] = S - dp[mask]\n",
    "\t\t\td += dp*dp\n",
    "\n",
    "\t\td = sqrt(d)\n",
    "\t\td[index] = 2 * rMax\n",
    "\n",
    "\t\t(result, bins) = histogram(d, bins=edges, normed=False)\n",
    "\t\tg += result\n",
    "\n",
    "\tg = g/(num_particles * numberDensity)\n",
    "\n",
    "    # Average g(r) divide by apropriate radii\n",
    "\tfor i in range(num_increments):\n",
    "\t\tradii[i] = (edges[i] + edges[i+1]) / 2.\n",
    "\t\trOuter = edges[i + 1]\n",
    "\t\trInner = edges[i]\n",
    "\t\tg[i] = g[i] / (4.0 / 3.0 * pi * (rOuter**3 - rInner**3))\n",
    "\n",
    "\treturn (radii, g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let me now try this algorithm on the points generated above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0xae02f34c>"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEZCAYAAACJjGL9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XGd97/HPLNp3yVq8K5L9OHYsO4mzEUI2CEtJCqQJ\nTYDQ5LKES3sLpfdyU1poL9DlVQovlhbahi1sgToQCEkaAlnIQvZNju3ksSx5kTW2pdE+o3Xm3D9m\nRhKOVlszZ5bv+/XyK9Yc6cwvx6P5zrOc5wEREREREREREREREREREREREREREZGs4XHriY0xecC3\ngfVAAfB5a+0vZxy/Cvg0MAl821r7TVcKFRHJcV4Xn/u9QLe19mLgrcC/Jg7EQ+RLwBXAJcCHjTF1\nrlQpIpLj3AyKncBnZtQxOePYZqDNWjtgrZ0AHgMuTnF9IiIC+N16YmttCMAYU0YsNP56xuFyYGDG\n10NAReqqExGRBDdbFBhj1gIPAt+z1v54xqEBoGzG12VAXyprExGRGNdaFMaYeuB+4KPW2odOOPwK\nsNEYUwWEiHU7fWG+80WjUcfjcW1sXkQkI3kW8cbpWlAAnyLWnfQZY0xirOJWoMRae6sx5hPAr4i1\ner5lrQ3MdzKPx0N391BSC84UtbVluhZxuhbTdC2m6VosTdZ8BHccx9E/fIx+CabpWkzTtZimazGt\nrq58wRxwdYxCRETSn4JCRETmpaAQEZF5KShERGReCgoREZmXgkJEROaloBARkXkpKEREZF4KChER\nmZeCQmQej7zUxQPPdbpdhoir3FzrSSStOY7DzofaGB2PcOHWBooK9OsiuUktCpE5DIYnCI1OEok6\n7D2oVe4ldykoROYQ6AlN/b11f9DFSkTcpaAQmUOgNzz1913tQRzHcbEaEfcoKETmkGhRNFQX0zc0\nRmd3aIGfEMlOCgqROQSCsWB487lrAWjd3+NmOSKuUVCIzKErGKaqrIAdm2rxALs0TiE5SkEhMouR\nsUn6hsZYVVNMWXE+TavKaTsySHh0wu3SRFJOQSEyi6PxgeyGmhIAWppriDoOuw9omqzkHgWFyCy6\n4gPZq2qKAdjWXANonEJyk4JCZBaBYKxFsTLeolhXX0Z5ST672nuJapqs5BgFhcgsEjOeVq6IBYXX\n46GlqZrB0DiHjg25WZpIyikoRGYRCIYpKfRTXpw39di25hWA7tKW3KOgEDnBZCTK8b4RVtaU4PF4\nph4/o7EKr8ejabKScxQUIic41jdC1HFYGR/ITiguzGPDmgrauwYZCo+7VJ1I6ikoRE6QWLojMZA9\n07bmGhzg5Y7eFFcl4h4FhcgJEgPZq1YUv+bYtqbYNFl1P0kuUVCInCAxNbZhlhbF6toSqsoK2NUe\nJBrVNFnJDQoKkRN0BUPk+b2sKC98zTGPx0NLUw2h0UnaA4MuVCeSegoKkRmijsPRYJiG6mK8Xs+s\n35O4S1vdT5IrFBQiM/QOjjI+GX3NjKeZNq+vwuf10NquoJDcoKAQmSExPrFqlvGJhKICP2ZtJQeP\nDjEwPJaq0kRco6AQmWFqauyKuYMCZnQ/tWuarGQ/BYXIDF1TiwHO3fUEM1aTVfeT5AAFhcgMgWAI\njwfqq+YPiobqYlZUFLK7o5fJSDRF1Ym4Q0EhMkMgGKausog8//y/Gh6Ph23NNYyMTbL/yECKqhNx\nh4JCJG4wPM7wyMSsS3fMRt1PkisUFCJx02s8zd/tlLBpXRV5fq/up5Csp6AQiQv0/v6udgspyPNx\n+roqOrtD9A6OJrM0EVcpKETiAj3xoJhlMcC5tDRVA+p+kuymoBCJm9r+tHpxLQqAFi3nITlAQSES\nFwiGqCzNp7jQv+ifqa8qpr66mD0H+5iY1DRZyU4KChFgdHyS4ODYoscnZtrWVMPYeIR9nf1JqEzE\nfQoKEeBo78JrPM1lapqsup8kSykoRJgeyG5Y5NTYmczaSvLzvOzSgLZkKQWFCLHNigBWnURQ5Pm9\nbFlfTSAY5nj/yHKXJuI6BYUI08uLL7Rq7Fy0mZFkMwWFCLEZT0UFfipK8k/q51uaEsuOKygk+yx+\nHmCSGGPOB/7JWnvZCY//BfABoDv+0M3WWpvq+iT7TUaiHO8bobGhDI9n9u1PF1JTUcjq2hL2Huxj\nfCJCfp5vmasUcY+rQWGM+STwPmB4lsNnAzdYa19IbVWSa7r7R4hEnZOaGjvTtqYa/vupQ7xyqH+q\nK0okG7jd9dQGXA3M9jFuB/ApY8yjxphbUluW5JKuk1i6YzYap5Bs5WpQWGt/BkzOcfh24GbgcuAi\nY8zbU1aY5JSppTtOsUXRvLqCogIfre09OI6zHKWJpAXXxyjm8RVr7SCAMeYe4Czgnvl+oLa2LBV1\nZQRdi2kLXYve0DgAWzfWUXuSs54SztpUx+9aA0x4vKyuLT2lcyWDXhfTdC0WLy2DwhhTAbQaY7YA\nYWKtim8t9HPd3UPJLi0j1NaW6VrELeZadBwZwO/z4o1ETvm6bVpdwe9aAzz8zCHefO7aUzrXctPr\nYpquxdK4PUaR4AAYY643xnzIWjsA3AI8BDwCvGytvc/NAiU7RR2Ho8EwDdVFeL0nN+Nppq2JabL7\ne075XCLpwvUWhbX2AHBh/O+3z3j8dmLjFCJJ0zc4xthE5JTHJxKqygpYV1/Kq4f7GR2fpDDf9V8x\nkVOWLi0KEVcEepe2/elibGuuYTLisPdg37KdU8RNCgrJaYnFAFed4iD2TNuaVgCaJivZQ0EhOW25\npsbO1LSqnJJCP63tQU2TlaygoJCc1hUM4/FAQ3XRsp3T6/WwtamG3sExjvSElu28Im5RUEhOCwRD\n1FYUkedf3rWZtjXpLm3JHgoKyVnDIxMMhSdOarOihZzRVI0H7Xon2UFBITmrqyexWdHyjU8klBfn\nc9qqcvZ1DhAenWuVGpHMoKCQnDU9kL38LQqI7VERdRz2HOhNyvlFUkVBITnrVHe1W0hiNdlWbWYk\nGU5BITkrERQns0/2YqxvKKOsOI9d+zVNVjKbgkJyViAYoqIkn+LCvKSc3+vx0NJUw0BonEPHZtub\nSyQzKCgkJ41NRAgOjCZtfCJB3U+SDRQUkpOOBsM4JG98IuGM06rxeHQ/hWQ2BYXkpMSMp2RMjZ2p\npDCPDasr2N81wPDIRFKfSyRZFBSSk7riA9nJuNnuRNuaa3AceLlDrQrJTAoKyUmpalFA7H4KUPeT\nZC4FheSkQDBMYb6PytL8pD/X2rpSKkvz2dXeSzSqabKSeRQUknMi0SjHesOsrCnB4zn17U8X4vF4\n2NZcw/DIBB1HB5P+fCLLTUEhOae7f5RI1EnajXazadFmRpLBFBSScwLxxQCTPTV2pi2NVfi8Hnbp\nfgrJQAoKyTldSV4McDZFBX42rqmgIzDEQGg8Zc8rshwUFJJzptd4Sl2LAmBbc6z76WW1KiTDKCgk\n5wSCIfw+DysqC1P6vC3x5TzU/SSZRkEhOcVxHALBMPVVxfi8qX35r6oppqa8kJfbe4lEoyl9bpFT\noaCQnNI3NMboeCSl4xMJiWmy4bFJ9h/RNFnJHAoKySmB3vhmRSken0hQ95NkIgWF5JTpqbGpb1EA\nbF5Xhd/npVX3U0gGUVBITnFrxlNCQb6P09dVcvj4MH1DY67UILJUCgrJKYFgCA/QUO1OiwLU/SSZ\nR0EhOaUrGKamopD8PJ9rNUzteqfuJ8kQCgrJGaHRCQZD46xK4dIds6mvKqa+qog9B3qZjGiarKQ/\nBYXkjEBPfLMiF7udElqaahgdj7Cvc8DtUkQWpKCQnJFY48ntFgVMdz9pNVnJBAoKyRkBFxYDnMum\ndZXk+720akBbMoCCQnJGYmqsWzfbzZTn97F5fRVdPSF6+kfcLkdkXgoKyRmBYIjy4jxKi/LcLgWY\n0f2kVoWkOQWF5ITxiQg9/aNp0ZpIaGnSNNmlchyH4/0jWlQxxfxuFyCSCkd7wzikdle7hayoLGLV\nihL2HuxjYjJCnt+9ezsyxW+e6+T23+yjIN/HxtUVbFpXyaa1VTSuLMPv0+feZFFQSE6YHp9wfyB7\npm1NNdz39CFePdTP1ngLQ2Y3GB7n5492UFTgo7K0gJc7enm5oxeA/DwvG1ZXsGltJZvWVXHaynLy\n/AqO5aKgkJyQmPHk1hpPc2lpjgVF6/6ggmIBv3i0g5GxSa5/00auOGctA8NjvHq4n1cP92MP9bPn\nQB97DvQBHeT5vTSvKsfEg6N5Vbmrd+NnuiUFhTGmAHgd0AisACLAMeAQ8KS1dnK5CxRZDl1p2qLY\nuKaCwnwfre1B3uN2MWnsSPcwD794hIbqYi47azUAFaUFnLe5nvM21wOxFse+w/28eigWHq8e6ueV\nQ/3w+AH8Pg+nrSyPdVWtq+KC8iI3/3cyzoJBYYzxAFcCfwZcCsw1ZWTYGPMg8E1r7d3LVqHIMggE\nQxTk+6gqK3C7lN/j93k5o7Ga52w3x3rD1KfBXePpxnEcfvxgG44D7758w5xjEeXF+ezYVMeOTXUA\nDI9MsK9zOjjajgywr3OAu393EN9/vUTjyjI2ra1i07pKNqyuoKhAHSxzmffKGGPeBnwFaAKeAL4A\n7ALagUFis6ZqgDXABcBFwF3GmFeAT1lrf5680kUWJxKNcqw3zJraUjwej9vlvEZLcw3P2W5a24Nc\noaB4jV3tQXZ39HJGYxXbmxffPVdalMdZG2s5a2MtAOHRyVhwHO6nPTBI2+EB9h8Z5N4nD+L1eFjf\nUDoVHBvXVFJcqOBImPNKGGN+CryBWFB821obWOBcP47/3Ebg/cA3jTHvt9ZevVzFipyMnoFRJiNO\nWk2NnSkxTXbX/iBXnLPW5WrSy2Qkyk8ebMPjgT9+48ZTCvriQj/bN6xg+4YV1NaWcaizj/1HBqa6\nqToCg3QEhrjv6UN4PLCuriw+q6qSjWsr0+b+GzfMF5kvAO+31oaWckJr7T7g08aYLwAfP5XiRJZD\nYjHAVS7tareQqrIC1taV8sqhfsbGIxTka9A14bcvdhEIhrn0zFWsqS1d1nMXFfjZ2lQzNYlgbCIS\nC454V1V71wAHjw1x/zOH8QCra0ungsOsq6S8OH9Z60lncwaFtfbzM782xpwB7LXWLupOF2vtIPDZ\nUytP5NRNr/GUni0KiN2lffj4MHsP9XHmhhVul5MWQqMT/PzRdooKfLzzDU1Jf76CPB9bGqvZ0lgN\nxG7SbO8ajLc4+tjfNUhn9zAPPNcJwPlb6vnglZvxebN/Gu5SOuEeBL4D3JKkWkSSoiuNFgOcS0tT\nDfc8cZBd+4MKiri7HjtAaHSSay9rprwk9Z/e8/N8nL6+itPXVwGnMTEZpSMQC45nXznOU3uOUVaU\nx3uuMCmvLdWWEoUlQEeyChFJlkAwjM/roa4qfadENq8up7jAT+v+II7juF2O6472hnnw+U5qKwt5\n0470GLfJ83sxayu56sJGbnnv2ayuLeE3z3VOtTCy2VKC4svAJ4wx5y5nAcaY840xD83y+FXGmKeN\nMb8zxnxwOZ9TcofjOASCIeqri9O6i8Dn9bK1qZrg4OjUPR+57L8ebCMSdbj20g1peYd1UYGfj12z\njfKSfH70G8tLbT1ul5RUS/kX2AGsAp4yxgwbYw4aY9pn/OkwxrQv5cmNMZ8EbgUKTng8D/gScAVw\nCfBhY0zdUs4tAtA7OMrIWISVGTDtdObsp1y250AvL7b1YNZWsmNTrdvlzGlFRRF//kfb8Pu8/Ptd\nuzl0bMjtkpJmKUFRBDwLPBL/bwexO7ITfw7G/yxFG3A1cOKct81Am7V2wFo7ATwGXLzEc4vQeWwY\ngJVpOuNppunVZLP70+l8olGHHz/Qhge4/hSnw6ZC06pyPnTlFsbGI3zljlb6hsbcLikpFj2Yba29\ndLmf3Fr7M2NM4yyHyoGZmwkPARXL/fyS/TqPxz7lpfOMp4TyknxOW1nGvs4BRsYmc/JO4Udbu+js\nHub1LQ2sbyhzu5xFOef0Oq65tJk7Ht7PV3/ayi3vOTvrpjjPd8Pd66y1T5zKyY0xF1lrHzuJHx0A\nZr5KyoC+hX6otjYzXlipoGsRc/jR2PyLLRtqM+KaXNCyio7Aq3T2jnDhtlXLfv50vgbh0Ql+8dgB\nCvN9fOhd26ipSO7kg+W8Fu+/8gwGwhP8+ulD3Hb/q9zyJ+fh86Z3a2gp5vvI8hNjzIvA31trn1rK\nSY0xlxObRns6sO4k6noF2GiMqQJCxLqdvrDQD3V3Z28f4VLU1pbpWsQdjvcbF3oy4/XRHP8U/fiL\nnWxcubxv6un+utj5cBv9w2O86w2nER2fTGqtybgW117SROexIZ58+Sj/vvNF3n35hmU9v5vmC4rN\nwN8BjxpjDgK/AO4BWq21vzfaFh9ovoDYkh/vBlYC3yA2/rAYTvw81wOl1tpbjTGfAH5FbBzlW4tY\nQkTkNTqPD1FTXpgxXQGNK8soLcqbmiab7n30y6W7f4RfP3OY6vIC3nLeyXy2dJ/f5+Wj79rKP3z/\nOe57+hB11UVceuZqt8taFvPdmR0C/o8x5hvA/wI+CHwCwBgzTGxRQB9QzfSKsv3Ad4EvW2sPLaYA\na+0B4ML432+f8fjdgFahlZMWHp2kd3CMrU3VbpeyaF6Ph5amap7YfYzDx4dZV5++XUXLaedDbUxG\nHK65pDmj940oKczjY9du5/O3PcsPfmWprSjijNMy5/U3lwVnPVlr2621f0GslfAW4HPEWhetwHPA\nD4i1PC4Gaq21n1hsSIgkU7puVrSQlvgKqbvac2OarD3cz7OvdtO8qpzzt9S7Xc4pq6uMTZv1ej18\n/ee7ONI97HZJp2wps55GgF/H/4ikvcTSHQ1pvHTHbLaeVoPHA637g7z9dY1ul5NUUcfh9gf2AXBd\nBkyHXawNayr4H28/nf+8aw9f3tnK3/zJOVS4sAzJcll0UBhjOoiPJczBAcaA48DTwBettcdOrTyR\nk5fYJzvTWhSlRXk0r6qg7cgAodEJSgqzd3nrJ14+ysGjQ5y/pZ7m1dk1A/6CLQ0c7xvh54928LWf\ntvLJ68/K2G61pdxw9wCx+xsaiQXCS8BTQG/8sZVAD7GNjP4SeMkYk5mjUpIVAj3pvxjgXFqaa3Ac\n2N3R63YpSTM2HuGnv91Pnt/LNZc0u11OUlx1YSOvO6OB9q5BvnnPXqIZuo7XUoLieWILA77TWrvZ\nWvsua+311tpziA1GTwK3WWtbgDOJtTA+P/fpRJIr0BumvCSfsgzcN2Db1F3a2TtO8d9PHaR/eJy3\nnLeOmopCt8tJCo/Hw41vOx2zpoJnXznOnY8saZWjtLGUoPhL4KvW2rtOPGCtfRL4KvBX8a93AV8H\n3rQcRYos1cRkhO7+EdZm6KyhdfWlVJTms6s9mLGfQufTOzjKfU8doqI0nz+4ILs7HvL8Xv7sj7ZR\nV1XEPU8c5NHWLrdLWrKlBEU9MN96ut3E9s5OCACVJ1OUyKk61juC48CauuXdFS1VPB4PLU01DIUn\nOHg0fW+SO1l3/HY/45NR/ujiZgrzs3+pktKiPD5+7XZKCv18775X2XtwwYUm0spSgmI3cKMxpuDE\nA8aYfOAGYndUJ5zN0hcJFFkWiRlPmdqigOztfmrvGuTJ3cdYX1/GhS0NbpeTMg3VxfzZ1S0AfP3O\nXVPTtzPBUoLi74DtxAapP26Mebsx5k3GmD8FngTOIj4mEb9J70PA95e5XpFFScx4WluXuUGxpbEa\nn9eTVfdTOI7D7Q9YAK574wa8WTIddrE2ravixredTmh0kq/sbGUoPO52SYuy6KCw1v438E5iU2q/\nBPwSuB/4GrGZTu+21t5hjKkFPgD8EPiXZa9YZBESn9bW1Gdm1xNAcaGfjWsq6OgaZDBD3lAW8vTe\n4+w/MsgOU8umdVVul+OK17es5MoLGzneP8LXfraLicmo2yUtaElbR1lr77bWbiA2q+la4D3E1nhq\ntNb+NP5tQWLrNd1orc2OV7dknK6eMPl5XlYkeQXSZGtpqsEBdrdn/jTZ8YkIdzzcht/n4drLsnM6\n7GK98w2ncd7mOto6B/jOvXvTfvvbkxpFsta2ElvCY7ZjUUABIa6JRh2O9oZZvaIEb4Yv9dzSXMPO\nh/fT2h7kdVszuz///mcOExwc463nr6OuKvPubVlOXo+HD7x9M8HBUZ7cc4z66mLecdFpbpc1p/Tb\njFbkFPUMjjIZiWbErnYLWb2ihOryAl5uDxKNpvenzvn0D49xzxMHKSvO48osX5ZksfL8Pv7X1dtY\nUVHILx7r4IndR90uaU4KCsk603dkZ9bSHbPxeDxsa6ohNDrJSxm8RerPHmlnbCLCu97QRHFh9k+H\nXazyknw+du12igr8fOfevdjD/W6XNCsFhWSd6TWeMr9FAfC6rQ14PR6+fufLPPh8Z9r3Z5/o4NEh\nHm8NsLq2hDdsX+l2OWln9YoSPvqurTgO/OvPdnGsL+x2Sa+hoJCsk7iHIhtaFAAb11Tyv687k+JC\nPz+43/Lte/cyPhFxu6xFcRyHHz+wDwe47vKN+Lx6y5nNGY3VvO/NhuGRCb68s5XhkQm3S/o9+leT\nrBMIhvB5PdRVZfaMp5lOX1/F3954Lo0NZTy+6yj/+MPnCQ6Mul3Wgp63Pbx6uJ/tzTVZsYFPMl1y\n5mreev46jvWG+fqdu5iMpM+0WQWFZBXHcQj0hKmrKsLvy66Xd3V5IX/1vrO5aNtKDh4d4v999xn2\nHkjfabMTk1F2PtSGz+vJqv2jk+maS5vZYWp55VA/t933Stp0M2bXb5LkvMHQOOGxSRqqs2N84kR5\nfh83ve10bnjLJkbGJvmXn7zIfU8dSps3lJkeeK6T4/0jXHbW6qzpBkw2r8fDB6/aMtVyvPfJ9FgF\nSUEhWaUrMZC9InvfmDweD5edtZr/+96zKS/J578eauM/7trN2Hj6jFsMhsf55e86KCn084dpfH9A\nOirI8/Hn12yjuryAn/62naf3ur//m4JCskogmLmbFS3VhtUV/O2N57JhTQVP7z3O33//WY6nyYyZ\nXzzawchYhD+86DRKi7J3h75kqSwt4OPXbKcw38c3797L/iMDrtajoJCskpgamytdHZWlBXzy+rO4\n/OzVdHaH+Ox3n3V9tdnO7mEefvEIDdXFXHbWaldryWRr6kr5yDu2EolG+dpPW+npH3GtFgWFZJVc\nalEk+H1e3vfmTXzg7ZsZn4zylZ0v8cvHO1zZ8MhxHH7yYBuOA+++fEPWTShItW3NNbz3CsNgeIIv\n39FKeHTSlTr0ryhZJRAMU11ekBOb4Zzo9S0r+dQNZ1NdXsCdj3bwbz/bxchYat9YdrUH2d3RyxmN\nVWxvrknpc2ery89ew5vOWUNXT4hv/NydabMKCskaI2OT9A2N5Uy302waG8r59I3nsnl9FS/s6+Fz\ntz1LV09qNsiZjET5yYNteDzwx2/ciCfH9ppIpusu38j25hp2H+jjh7+2KZ/lpqCQrDE9PpE73U6z\nKS/O5xN/vJ23nreOo71hPve9Z3nu1e6kP+/DLxwhEAxzyfZVrKnN3H1A0pHX6+Hmd5zBurpSfvti\nF796+nBqnz+lzyaSRInxiVU53KJI8Hm9vPvyDXzkHWfgOA7/ducuvnfvnqStQBsaneAXj3VQVODj\nnW9oSspz5LrCfD9/fs02Kkvz2flQG8/b5Id/goJCskZXDg5kL+S8zfX8zQ3nUFtZyM4H9vHlO15K\nyjpCdz12gNDoJFde2Eh5Sf6yn19iqssL+dg128nL8/Kfd+3mwNHBlDyvgkKyRqAnt6bGLtaaulI+\nc+O57Di9jpfbe/ncbc9w6NjQsp0/EAzx4POd1FYW8qYda5ftvDK79Q1l3PyHZzAxGeUrO1sJjSZ/\nAUEFhWSNQG+YkkI/ZcW6wetEJYV5fPoDF3DlhY1094/yD99/jif3LM9GOTsf2k8k6nDtpRvI8+st\nJRXO2ljLe64wRKJOSlaazb05hJKVJiNRuvtGaFpdrtk2c/B5PVx9cRONDWV88+49/OddezgQGOLa\ny5pPevnv3Qd6ebGtB7O2kh2bape5YpnPG3es4fKzV6fk9a74l6xwrDdM1HGyZrOiZDrb1PLpPzmH\nlTXF3P/MYb744xcZDC99m/to1OEnD+zDA1yv6bCuSNU1V1BIVsi1pTtO1cqaEv7m/edwdnxJ689+\n9xk6AksbGH2ktYvO7hAXtjSwvqEsSZVKOlBQSFbItl3tUqGowM9H37WVqy9uom9wjH/8wfM82tq1\nqJ8dGZvk54+0U5Dn4+qLm5NcqbhNQSFZIdv2yU4Vr8fDlRc28vF3byff7+U7977C93/16oLLRNz9\nxAEGwxP8wQXrqCorSE2x4hoFhWSFQE+IfL+X6opCt0vJSC1NNXzmxnNYU1vCQy8c4Z9/9AJ9Q2Oz\nfu/x/hF+/cxhqssLeMt561JcqbhBQSEZL+o4HO0N01BdjFcDqietrqqYv77hHM7fUk/bkQE++91n\n2NfZ/5rvu+OhNiYjDtdc2kx+ns+FSiXVFBSS8YIDo4xPRlmZxbvapUpBvo8PX7WF6y7fwFB4gn/+\n0Qs8+Hzn1CJ09nA/z77aTfOqcs7fXO9ytZIqCgrJeFoMcHl5PB7efN46/vK6Mykq8POD+y3fvncv\nYxMRbn9gHwDXaTpsTtENd5LxtBhgcmxeX8Xf3ngu/3bnLh7fdZQ9B/roGxrj/C31NK+ucLs8SSG1\nKCTj5eKudqlSU1HIX73vbC5qWUnf0Bh5fi/XXKLpsLlGLQrJeF3BMF6Ph/pqBUUy5Pl93PQHp7N9\nQw2FBX5qNLMs5ygoJKM5jkOgJ0RtVZH2Z04ij8fDjk11bpchLtFvlmS0ofAEodFJ3WgnkkQKCslo\nifGJBgWFSNIoKCSjdU0t3aEZTyLJoqCQjBbo0WKAIsmmoJCMFujVzXYiyebarCdjjBf4OrANGAM+\naK3dP+P4XwAfALrjD91srbUpL1TSWiAYoqqsgKICTeATSRY3f7veCeRbay80xpwPfDH+WMLZwA3W\n2hdcqU7S3uj4JL2DY2xprHK7FJGs5mbX0+uB+wCstU8B55xwfAfwKWPMo8aYW1JdnKQ/7Wonkhpu\nBkU5MHNI1GGjAAAHd0lEQVTvxUi8OyrhduBm4HLgImPM21NZnKS/6TWeND4hkkxudj0NAjM32vVa\na2duq/UVa+0ggDHmHuAs4J75Tlhbq317E3LhWgyMHAZgc3PtvP+/uXAtFkvXYpquxeK5GRSPA1cB\nO40xFwCtiQPGmAqg1RizBQgTa1V8a6ETdncPJanUzFJbW5YT16LtUB8ARb65/+1z5Voshq7FNF2L\npXEzKO4ErjDGPB7/+iZjzPVAqbX21vi4xEPEZkT9xlp7n1uFSnoKBMMUF/gpL8l3uxSRrOZaUFhr\nHeB/nvjwjOO3ExunEHmNyUiU7v4RGleWaQMdkSTTDXeSkY73jRCJOprxJJICCgrJSNrVTiR1FBSS\nkbq0T7ZIyigoJCNNbX+6Qi0KkWRTUEhGCvSEyfN7WVGubTlFkk1BIRkn6jgEekM0VBfj9WrGk0iy\nKSgk4/QOjjI+EdX4hEiKKCgk42gxQJHUUlBIxgloxpNISikoJOPoHgqR1FJQSMYJ9ITweKC+Wi0K\nkVRQUEjG6QqGqa0sIs+vl69IKug3TTLKUHic4ZEJdTuJpJCCQjKKBrJFUk9BIRmlKz6Q3aCgEEkZ\nBYVklEBPrEWhrieR1FFQSEYJ9MYXA1RQiKSMgkIySqAnTEVpPsWFbu7iK5JbFBSSMcbGIwQHR9Xt\nJJJiCgrJGEd7NeNJxA0KCskYiRlPGp8QSS0FhWSM6TWe1KIQSSUFhWSMxNRYbX8qkloKCskYXcEQ\nRQU+Kkry3S5FJKcoKCQjTEaiHO8bYWVNCR6Ptj8VSSUFhWSE7v4RIlFHM55EXKCgkIyQWAxQ91CI\npJ6CQjJCQFNjRVyjoJCM0DU140ldTyKppqCQjBAIhvD7vNRWFLldikjOUVBI2nMch0BvmIbqIrxe\nzXgSSTUFhaS9vqExxsYjNGh8QsQVCgpJe11aukPEVQoKSXtTS3eoRSHiCgWFpL2AlhcXcZWCQtJe\noCeEB2ioVlCIuEFBIWkvEAyxorKQ/Dyf26WI5CQFhaS14ZEJBsMTGp8QcZGCQtLa9GZFCgoRtygo\nJK0lFgPUQLaIexQUkta6erQYoIjbFBSS1qZaFFoMUMQ1CgpJa4FgiPKSfEoK89wuRSRnKSgkbY1P\nRAgOjGrpDhGXKSgkbR3tDeOg8QkRtykoJG11Te1qpxaFiJsUFJK2phYDXKEWhYib/G49sTHGC3wd\n2AaMAR+01u6fcfwq4NPAJPBta+03XSlUXKOb7UTSg5stincC+dbaC4FbgC8mDhhj8oAvAVcAlwAf\nNsbUuVKluCYQDFOY76OyNN/tUkRymptB8XrgPgBr7VPAOTOObQbarLUD1toJ4DHg4tSXKG6JRKMc\n7Q2zsqYYj0fbn4q4ybWuJ6AcGJzxdcQY47XWRuPHBmYcGwIqklmM4ziERieT+RQpUxAaZ3hkwu0y\nTkl3/wiRqKMZTyJpwM2gGATKZnydCAmIhcTMY2VAXzKLue2+V3jkpUAyn0JOgmY8ibjPzaB4HLgK\n2GmMuQBonXHsFWCjMaYKCBHrdvrCfCfzqH8iK939JbcrEBHX3lyNMR6mZz0B3ATsAEqttbcaY64E\nPkNsHOVb1tpvuFOpiIiIiIiIiIiIiIiIiIiIiIiklYyfUrrQmlG5yBhzPvBP1trL3K7FLfFlYL4N\nrAcKgM9ba3/pblXuMMb4gFsBAzjAR6y1u92tyl3xJYGeA95orbVu1+MWY8zzTN/c3G6t/cBs3+fm\nfRTLZWrNqPgb5Bfjj+UkY8wngfcBw27X4rL3At3W2hvi9+O8CORkUABXAlFr7UXGmEuAvye3f0fy\ngP8gdo9WzjLGFAIs5gNlNiwzPt+aUbmoDbiaLGgtnqKdxO7DgdjrPDvWZzkJ1tpfADfHv2wkyasc\nZIAvAN8Acn0phu1AsTHmV8aYB+IftGeVDUEx65pRbhXjNmvtz8jhN8UEa23IWjtsjCkjFhp/7XZN\nbrLWRowx3wW+CvzI5XJcY4y5kVhL8/74Q7n8gSoEfMFa+xbgI8AP53rvzIY31PnWjJIcZoxZCzwI\nfM9a+2O363GbtfZGYuMUtxpjilwuxy03AVcYYx4CzgRuM8bUu1yTWyzwQwBr7T4gCKyc7RuzYYxi\nvjWjJEfFf/nvBz5qrX3I7XrcZIy5AVhjrf1HYASIxv/kHGvtJYm/x8PiZmvtMRdLctNNxCYB/akx\nZhWx3plZu+OyISjuJPYJ4fH41ze5WUwacdwuwGWfIrY0/WeMMYmxirdZa0ddrMktdwDfNcb8FsgD\nPmatHXO5JnHft4DvGGMeiX99k3pjREREREREREREREREREREREREREREREREREREREREREREkiKX\nl9gVSSpjzAFiCxP6gPcAPcB2a22vm3WJLFU2LAookq4c4HpgD/DnQINCQjKRgkIkeTxAIfAOa+1R\nt4sROVnZsHGRSDprU0hIplNQiCTXcbcLEDlVCgqR5Iq4XYDIqVJQiIjIvBQUIiIyLwWFiIjMS0Eh\nkjyO2wWIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiOv+P0OP0tN5lyEcAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xae0290cc>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute self-pair correlation\n",
    "dr = 0.5\n",
    "rad_py, g_r_py = py_rdf(r, Lbox, dr)\n",
    "\n",
    "plt.plot(rad_py, g_r_py)\n",
    "plt.xlabel('r', fontsize=18)\n",
    "plt.ylabel('g(r)', fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $g(r)$ is only computed to $L_{box}/2$ = 5, due to nearest image convention."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## F2PY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I am going to write a Fortran 90 subroutine called gr.f90 which does the expensive calculation in Fortran. I use the same algorithm, so that I expect an asymptotic scaling that goes as the square of the number of particles.\n",
    "\n",
    "The \"%%file FILENAME\" magic command, write a file in the same directory from where this jupyter instance is fired."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write Fortran Subroutine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting gr.f90\n"
     ]
    }
   ],
   "source": [
    "%%file gr.f90\n",
    "! Fortran File gr.f90\n",
    "    subroutine grSelf(npart, num_increments, r, S, radii, gr)\n",
    "    \n",
    "        integer, intent(in)     :: npart\n",
    "        integer, intent(in)     :: num_increments\n",
    "        real, intent(in) \t\t:: r(npart, 3)\n",
    "        real, intent(in)        :: S\n",
    "\n",
    "        real, intent(out)      :: radii(num_increments)\n",
    "        real, intent(out)\t   :: gr(num_increments)\n",
    "\n",
    "        ! local variables\n",
    "        integer :: i, j, k, ind\n",
    "        real    :: dr, rMax, numberDensity, d, dp\n",
    "        \n",
    "        rMax   = S/2.0\n",
    "        dr     = rMax/(num_increments)\n",
    "\t\t\n",
    "        numberDensity  = npart / S**3                \n",
    "\n",
    "\t\t! Compute pairwise correlation for each particle\n",
    "\t\tdo i = 1, npart\n",
    "\t\t\n",
    "\t\t\tdo j = i+1, npart\n",
    "\t\t\n",
    "\t\t\t\td = 0.\n",
    "\t\t\t\tdo k = 1, 3\n",
    "\t\t\t\t\tdp = abs(r(i,k) - r(j,k))\n",
    "\t\t\t\t\tif(dp > S/2.0) then\n",
    "\t\t\t\t\t\tdp = S - dp\n",
    "\t\t\t\t\tendif\n",
    "\t\t\t\t\td = d + dp*dp\n",
    "\t\t\t\tenddo\n",
    "\t\t\t\t\n",
    "\t\t\t\td   = sqrt(d)   ! distance between particles i and j\n",
    "\t\t\t\tind = int(d/dr) + 1\n",
    "\t\t\t\t\t\t\t\t\n",
    "\t\t\t\tif(ind <= num_increments) then\n",
    "\t\t\t\t\tgr(ind) = gr(ind) + 2.\n",
    "\t\t\t\tendif\n",
    "\n",
    "\t\t\tenddo\n",
    "\n",
    "\t\tenddo\n",
    "\n",
    "\t\tgr = gr/(npart * numberDensity)\n",
    "\t\t\n",
    "\t\tdo i = 1, num_increments\n",
    "\t\t\tradii(i) = (i - 0.5) * dr\n",
    "\t\t\trOuter =  i * dr\n",
    "\t\t\trInner = (i - 1) * dr\t\t\n",
    "\t\t\tgr(i) = gr(i) / (4.0 / 3.0 * 3.14159 * (rOuter**3 - rInner**3))\n",
    "\t\tenddo\n",
    "\t\n",
    "    end subroutine grSelf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general you can also add directives to f2py by using \"!f2py DIRECTIVE\" tags. These are ignored by the Fortran compiler, and interpreted by f2py. This is similar to OpenMP directives."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Signature and Module"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate python interface file `grFor.pyf` and module `grFor` for later importation into python. I use the `rm grFor.pyf` command to remove any existing version of grFor.pyf in the current directory that might mess up my process. The command is:\n",
    "\n",
    "`f2py -m module_name -h sig_file.pyf list_of_fortran_files`.\n",
    "\n",
    "Note, I am not compiling the source code yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading fortran codes...\r\n",
      "\tReading file 'gr.f90' (format:free)\r\n",
      "Post-processing...\r\n",
      "\tBlock: grFor\r\n",
      "\t\t\tBlock: grself\r\n",
      "Post-processing (stage 2)...\r\n",
      "Saving signatures to file \"./grFor.pyf\"\r\n"
     ]
    }
   ],
   "source": [
    "!rm grFor.pyf\n",
    "!f2py -m grFor -h grFor.pyf gr.f90"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I can take a sneak peek at the signature file. The file provides an interface so that the python program and Fortran subroutine can talk to each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "!    -*- f90 -*-\r\n",
      "! Note: the context of this file is case sensitive.\r\n",
      "\r\n",
      "python module grFor ! in \r\n",
      "    interface  ! in :grFor\r\n",
      "        subroutine grself(npart,num_increments,r,s,radii,gr) ! in :grFor:gr.f90\r\n",
      "            integer, optional,intent(in),check(shape(r,0)==npart),depend(r) :: npart=shape(r,0)\r\n",
      "            integer intent(in) :: num_increments\r\n",
      "            real dimension(npart,3),intent(in) :: r\r\n",
      "            real intent(in) :: s\r\n",
      "            real dimension(num_increments),intent(out),depend(num_increments) :: radii\r\n",
      "            real dimension(num_increments),intent(out),depend(num_increments) :: gr\r\n",
      "        end subroutine grself\r\n",
      "    end interface \r\n",
      "end python module grFor\r\n",
      "\r\n",
      "! This file was auto-generated with f2py (version:2).\r\n",
      "! See http://cens.ioc.ee/projects/f2py2e/\r\n"
     ]
    }
   ],
   "source": [
    "!cat grFor.pyf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Compile the File"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compile the file with the \"-c\" flag. You may see a bunch of warnings thrown up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[39mrunning build\u001b[0m\n",
      "\u001b[39mrunning config_cc\u001b[0m\n",
      "\u001b[39munifing config_cc, config, build_clib, build_ext, build commands --compiler options\u001b[0m\n",
      "\u001b[39mrunning config_fc\u001b[0m\n",
      "\u001b[39munifing config_fc, config, build_clib, build_ext, build commands --fcompiler options\u001b[0m\n",
      "\u001b[39mrunning build_src\u001b[0m\n",
      "\u001b[39mbuild_src\u001b[0m\n",
      "\u001b[39mbuilding extension \"grFor\" sources\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpU1hSmB/src.linux-i686-2.7\u001b[0m\n",
      "\u001b[39mf2py options: []\u001b[0m\n",
      "\u001b[39mf2py: grFor.pyf\u001b[0m\n",
      "Reading fortran codes...\n",
      "\tReading file 'grFor.pyf' (format:free)\n",
      "Post-processing...\n",
      "\tBlock: grFor\n",
      "\t\t\tBlock: grself\n",
      "Post-processing (stage 2)...\n",
      "Building modules...\n",
      "\tBuilding module \"grFor\"...\n",
      "\t\tConstructing wrapper function \"grself\"...\n",
      "\t\t  radii,gr = grself(num_increments,r,s,[npart])\n",
      "\tWrote C/API module \"grFor\" to file \"/tmp/tmpU1hSmB/src.linux-i686-2.7/grFormodule.c\"\n",
      "\u001b[39m  adding '/tmp/tmpU1hSmB/src.linux-i686-2.7/fortranobject.c' to sources.\u001b[0m\n",
      "\u001b[39m  adding '/tmp/tmpU1hSmB/src.linux-i686-2.7' to include_dirs.\u001b[0m\n",
      "\u001b[39mcopying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpU1hSmB/src.linux-i686-2.7\u001b[0m\n",
      "\u001b[39mcopying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpU1hSmB/src.linux-i686-2.7\u001b[0m\n",
      "\u001b[39mbuild_src: building npy-pkg config files\u001b[0m\n",
      "\u001b[39mrunning build_ext\u001b[0m\n",
      "\u001b[39mcustomize UnixCCompiler\u001b[0m\n",
      "\u001b[39mcustomize UnixCCompiler using build_ext\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n",
      "\u001b[39mFound executable /usr/bin/gfortran\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler using build_ext\u001b[0m\n",
      "\u001b[39mbuilding 'grFor' extension\u001b[0m\n",
      "\u001b[39mcompiling C sources\u001b[0m\n",
      "\u001b[39mC compiler: i686-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC\n",
      "\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpU1hSmB/tmp\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpU1hSmB/tmp/tmpU1hSmB\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpU1hSmB/tmp/tmpU1hSmB/src.linux-i686-2.7\u001b[0m\n",
      "\u001b[39mcompile options: '-I/tmp/tmpU1hSmB/src.linux-i686-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c'\u001b[0m\n",
      "\u001b[39mi686-linux-gnu-gcc: /tmp/tmpU1hSmB/src.linux-i686-2.7/grFormodule.c\u001b[0m\n",
      "In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4,\n",
      "                 from /tmp/tmpU1hSmB/src.linux-i686-2.7/fortranobject.h:13,\n",
      "                 from /tmp/tmpU1hSmB/src.linux-i686-2.7/grFormodule.c:18:\n",
      "/usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning \"Using deprecated NumPy API, disable it by \" \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-Wcpp]\n",
      " #warning \"Using deprecated NumPy API, disable it by \" \\\n",
      "  ^\n",
      "/tmp/tmpU1hSmB/src.linux-i686-2.7/grFormodule.c:150:12: warning: ‘f2py_size’ defined but not used [-Wunused-function]\n",
      " static int f2py_size(PyArrayObject* var, ...)\n",
      "            ^\n",
      "\u001b[39mi686-linux-gnu-gcc: /tmp/tmpU1hSmB/src.linux-i686-2.7/fortranobject.c\u001b[0m\n",
      "In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4,\n",
      "                 from /tmp/tmpU1hSmB/src.linux-i686-2.7/fortranobject.h:13,\n",
      "                 from /tmp/tmpU1hSmB/src.linux-i686-2.7/fortranobject.c:2:\n",
      "/usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning \"Using deprecated NumPy API, disable it by \" \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-Wcpp]\n",
      " #warning \"Using deprecated NumPy API, disable it by \" \\\n",
      "  ^\n",
      "\u001b[39mcompiling Fortran sources\u001b[0m\n",
      "\u001b[39mFortran f77 compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran f90 compiler: /usr/bin/gfortran -Wall -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran fix compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -Wall -fno-second-underscore -fPIC -O3 -funroll-loops\u001b[0m\n",
      "\u001b[39mcompile options: '-I/tmp/tmpU1hSmB/src.linux-i686-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c'\u001b[0m\n",
      "\u001b[39mgfortran:f90: gr.f90\u001b[0m\n",
      "gr.f90:6.26:\n",
      "\n",
      "        real, intent(in)   :: r(npart, 3)\n",
      "                          1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:10.26:\n",
      "\n",
      "        real, intent(out)    :: gr(num_increments)\n",
      "                          1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:22.1:\n",
      "\n",
      "  do i = 1, npart\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:24.1:\n",
      "\n",
      "   do j = i+1, npart\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:26.1:\n",
      "\n",
      "    d = 0.\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:27.1:\n",
      "\n",
      "    do k = 1, 3\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:28.1:\n",
      "\n",
      "     dp = abs(r(i,k) - r(j,k))\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:29.1:\n",
      "\n",
      "     if(dp > S/2.0) then\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:30.1:\n",
      "\n",
      "      dp = S - dp\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:31.1:\n",
      "\n",
      "     endif\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:32.1:\n",
      "\n",
      "     d = d + dp*dp\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:33.1:\n",
      "\n",
      "    enddo\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:35.1:\n",
      "\n",
      "    d   = sqrt(d)   ! distance between particles i and j\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:36.1:\n",
      "\n",
      "    ind = int(d/dr) + 1\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:38.1:\n",
      "\n",
      "    if(ind <= num_increments) then\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:39.1:\n",
      "\n",
      "     gr(ind) = gr(ind) + 2.\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:40.1:\n",
      "\n",
      "    endif\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:42.1:\n",
      "\n",
      "   enddo\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:44.1:\n",
      "\n",
      "  enddo\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:46.1:\n",
      "\n",
      "  gr = gr/(npart * numberDensity)\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:48.1:\n",
      "\n",
      "  do i = 1, num_increments\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:49.1:\n",
      "\n",
      "   radii(i) = (i - 0.5) * dr\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:50.1:\n",
      "\n",
      "   rOuter =  i * dr\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:51.1:\n",
      "\n",
      "   rInner = (i - 1) * dr  \n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:52.1:\n",
      "\n",
      "   gr(i) = gr(i) / (4.0 / 3.0 * 3.14159 * (rOuter**3 - rInner**3))\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "gr.f90:53.1:\n",
      "\n",
      "  enddo\n",
      " 1\n",
      "Warning: Nonconforming tab character at (1)\n",
      "\u001b[39m/usr/bin/gfortran -Wall -Wall -shared /tmp/tmpU1hSmB/tmp/tmpU1hSmB/src.linux-i686-2.7/grFormodule.o /tmp/tmpU1hSmB/tmp/tmpU1hSmB/src.linux-i686-2.7/fortranobject.o /tmp/tmpU1hSmB/gr.o -lgfortran -o ./grFor.so\u001b[0m\n",
      "Removing build directory /tmp/tmpU1hSmB\n"
     ]
    }
   ],
   "source": [
    "!f2py -c grFor.pyf gr.f90"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are lucky and your program is compiled you have a usable Python module called grFor.\n",
    "\n",
    "Let us make sure we understand the python module by calling the `help(grFor)` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on module grFor:\n",
      "\n",
      "NAME\n",
      "    grFor\n",
      "\n",
      "FILE\n",
      "    /home/sachins/Desktop/Summer Learning 2016/scientific-python-lectures-master/grFor.so\n",
      "\n",
      "DESCRIPTION\n",
      "    This module 'grFor' is auto-generated with f2py (version:2).\n",
      "    Functions:\n",
      "      radii,gr = grself(num_increments,r,s,npart=shape(r,0))\n",
      "    .\n",
      "\n",
      "DATA\n",
      "    __version__ = '$Revision: $'\n",
      "    grself = <fortran object>\n",
      "\n",
      "VERSION\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import grFor\n",
    "help(grFor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes, Note that the order in which variables (input or output) are expected may differ from the Fortran subroutine. In our case, the Fortran subroutine had:\n",
    "\n",
    "`subroutine grSelf(npart, num_increments, r, S, radii, gr)`\n",
    "\n",
    "while the python interface expects:\n",
    "\n",
    "`radii,gr = grself(num_increments,r,s,npart=shape(r,0))`\n",
    "\n",
    "Note that npart is the last argument expected. Also note that the function is all lowercase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Run the Fortran code from python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "num_increments = 10\n",
    "rad_for, g_r_for = grFor.grself(num_increments, r, Lbox, npart)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us plot the results, and compare it with the plain python result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0xae312f4c>"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb8AAAEvCAYAAADVdc8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XGd97/HPzGi0j/aRZVu25MR+HHl3ErLYzlZICzRA\nSLjgFAIJSW6aBEIvpLwKFLqkBXpJ4CYhodmAlNLmNoWUQCktFxLSOIlDiOTdfrzEi7xp39eZOfeP\nkWR50TLSaI5m5vt+vfySZs7RmZ+P5fnO85znOQ+IiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIjMGh63\nCzDGXAp83Vp7zRnP3wR8BggB24C7rbWOCyWKiEiK8br54saYzwNPAllnPJ8D3A9cba3dABQC1yW+\nQhERSUWuhh+wD7iBs1ugfcDl1tq+occZQG8iCxMREZkxxphqY8xr42z/tDHm3xNZk4iIpLYMtwsY\nizHGC/xvYDFw40T7RyIRx+Nx/RKmiIi4zDOJMJi14Qc8TrT784OTGeji8XhobOyc+apSTDAY0Hmb\nAp232OmcTY3O28yYLeHnwMgIz3zgTeCTwMvAr40xAA9Za//NtQpFRCRluB5+1tqDwLqh7/951Caf\nKwWJiEjKc3u0p4iISMIp/EREJO0o/EREJO0o/EREJO0o/EREJO24PtpTREQS54Fna9l1sBWAmupi\n7tu41uWK3KGWn4hImnjg2Vp2HmzFITq5eufBVj736CYOnZj5SfR1dW+xf/8+AN7//j+Y8debiMJP\nRCRNDLf4Rmvt7OfhH22d8df+939/gaamRgBmw50o1e0p06ZuFJH09fOf/5TNm1+lra2d9vY23vWu\n3+c3v3mRJ598BoCvfOULbNz4UTZvfo29ey3V1YsYGBjkr/7qzzl58gSFhYXcf//f0dvby/33f5me\nnh7C4RB33HE3F154MZ/4xEbWrr2Iffv24vF4+PrXHyQvL3/adavlJ9NyqhvFwckYSGg3iojEpqa6\n+KznigNZ3Hvjqikf0+PxEIk4PPTQYzz44MM8//y/4vf7OXjwbTo62jl+/BjLlq3gssvWcffd9zJn\nTgW9vT3ceeeneOyxp+jq6mLv3j0888zTXHLJZXz7209w//1/x9e+dj8APT09vOtd7+bb336CYLCc\n119/dcq1jqbwk2kZbvFlzN9H9tpf4yurT1g3iojE5r6NaykOnFo7vDiQxYP3rKeqIjCt41500TsA\nKC0tIxAo4PrrP8TPf/4Cv/zlf/Lud7/3rP0LCgqpqKgAoKSklL6+Pg4fPsjq1dFeo7KyIHl5ebS2\ntgBgzFIAysvnMDAwMK1ahyn8ZPq8ITLmHMLjAf+i7fjKjrpdkYiM4d4bV1EcyJp2i2+03bt3AtDS\n0kxfXy9XXnk1b7yxmZdffonf//1o+Hk8HsLh8ND3Zx+jqmoRW7bUAtDY2EBXVycFBYUjPxtvuuYn\n01JTXcye3jo8GSFCTfPwFTWSuWgb76w83+3SROQcqioCPHjP+rges77+CJ/5zN309HRx331fIDs7\nmzVrLqS9vY1AINqqXLZsBY8//ijz5s0HTg8zj8fDzTffyte+9te89NKv6e/v4/Of/xI+n++sfeNl\nFoy5iQ/HcRyteRW76a4VFnEi3PuL+4n4e+mru5rC4jC+JW/QF+rj1uU3cdGcNXGsdvbQGmux0zmb\nmtl+3v7jP35GW1sbN930sdOe/9a3/jdXXfV7XHjhxQmvqby8YMJsU7enTMvO5j04Wd1kdFRSnBPg\nT/7wSj695nayfFl8f+ezvNWga38iqe7MXsnPfvZTdHZ2uhJ8k6VuT5mWF4+8AsDnr72BysC8oWcD\nfGrNbXy77im+t+Of8Hm8rA6ucK9IEZkx73nPdWc9981vftuFSmKjlp9M2fHuk+xu3cuSovNGBV/U\nosIq7l59GxneDJ7e/kO2Ne10qUoRkbMp/GTKXhpq9V29YMM5t59fVM3dqz6Jz+PlqW0/YEfz7kSW\nJyIyJoWfTEn3YA+bT7xFaXYxq8qWjbnfkuLz+ONVt+LxeHhi2z+wq9kmsEoRkXNT+MmUvHrsDQYj\ng1xZuQ6vZ/xfo6Uli7lz1S0APL7t++xp2ZeACkVExqbwk5iFI2F+U/8qmb5M1s29ZFI/U1Ni+J8r\nP4HjOHxn6/fY27p/hqsUkUQIh8N8+tN3ctddt9HV1TXuvh0dHfzyl79IUGXj02hPidmWph209rdx\n5fzLyfXnTPrnlpcu5faVN/Pkth/w2Nbvcc/q21hctGgGKxVJLz/e9zNqG7bF9Zhry1dyw+KzR3QO\na2xspKenh6ef/sGEx9q3z/LKKy9z7bXvjmeJU6Lwk5iNDHSpjP0uESvLlnHbio/x1PYf8NiWp/nU\nmjs4r7Aq3iWKSII88MBXqa8/zDe+8dWhIOw+bVWGm2/+MAsXVpGR4aejo519+/bywgvPs23bFjo6\n2uno6ODv/u6bPPbYwzQ0NNDc3MSGDVdyxx138bd/+5dkZmZy/Phxmpub+NKX/gJjLohL3Qo/icnh\nznr2tx9kWclS5uSVT+kYq4PL+eTyj/LdHT/k0bqn+fTa26kuWBjnSkXSzw2Lrxu3lTYT7rvvC/zF\nX3yR3Nw8LrnkPD70oY00NTVy112389xzP6Gvr49bbrmDJUsMtbW/4yc/+THvf/8H2b59KxdddAkf\n/vBNnDhxnBUrVnLdddfT39/PjTf+IXfccRcej4eKinn86Z9+kZ/+9N944YXnue++L8Slbl3zk5i8\ndGQTMPb0hslaW76SW5ZtpD/cz7frnuJwR308yhORBHMcB4BDh94ec1WGhQurRvYd3n/084FAgF27\ndvLXf/1lHnnkWwwMDI7sM7yiQzBYHrcVHUDhJzFo7+/kzZN1zMkNUlOyZNrHu2jOGj6xbCN9oX4e\nqXuSI53H4lCliLihunrsVRm83mjU+Hy+08JveLWGn//8Z+TnB/jKV+5n48aP0t/fN+P1qttTJu2V\no68RdsJcXbl+wukNk/WOirVEnAg/2PUvPFL3BJ9Zeyfz8+fG5dgikhgej4ePfWziVRnmz6/kwIF9\n/Mu//PPIzwFcfPEl/NVf/Tl79uyiomIuS5fW0NTUeNo+8V7WSKs6pLnJ3jF+MBLiy5u+SsgJ8Tfr\nvkR2RtaEPxOLV4/9lh/ufo58fx6fWXsn8/Ir4nr8eJvtd9qfjXTOpkbnLXZa1UHi5q2TW+gc7GLd\n3EviHnwA6+a9g5uW3kDXYDcP1z3Bie6GuL+GiMgwhZ9MyHEcXqx/BQ8erqpcN2Ovs2H+ZXzEXE/n\nQBcP1z7OyZ7GGXstEUlvCj+Z0P72gxzpPMqq4HJKc0pm9LWurFzHh5a8n/aBTh6ufYLGnuYZfT0R\nSU8KP5nQ8KT2a6YwqX0qrlmwgRsWX0dbfzsP1T5OU29LQl5XRNKHwk/G1dLXypamHczPn8viovMS\n9rrvXHgl15//Xlr723io9nGae1sT9toikvoUfjKul+tfI+JEuKZyQ9yHGk/k2qqred95f0BLXysP\n1z5Oa19bQl9fRFKXwk/G1B8eYNOxzeT787h4zhpXanh39Tt576Jraepr4aHax2nrb3elDhFJLQo/\nGdMbJ96iJ9TLhvmX4ff5XavjvdXv4t1Vv0djbzMP1T5Oe3+Ha7WISGpQ+Mk5OY7DS/Wb8Hq8XDH/\nMldr8Xg8XHfeH3Dtwqtp6Gni4don6BjQpF8RmTqFn5zT7ta9nOg+yYXlqyjKKnS7HDweDx84/z38\n3oIrONHTwCO1T9I5MP7CmSIiY1H4yTmNTG+Y5uoN8eTxeLhh8XVcXbmeY90neKTuSboGu90uS0SS\nkMJPztLQ08j25t0sKqiadevseTwePrTk/Vw5/3KOdh3n27VP0jPY43ZZIpJktKqDnOWl+lcBuGZB\nYia1x8rj8fA/zAcIOxE2HdvMI3VP8ek1d5Drz3G7NBGJ0QPP1rLrYHQeb011MfdtXJuQ11XLT07T\nG+rl9eO/pSirkDXBlW6XMyavx8vGpR/k8rnv4HBnPY9ueZre0MyvASYi8fPAs7XsPNiKAzjAzoOt\nfO7RTRw6MfMD2hR+cprXjr9Jf3iAK+dfjs/rc7uccXk9Xv7oghu5tOIiDnYc5rEtT9OnABSZ9SJO\nhOPdJ9nTtR1/1U4ya14no+IAAK2d/Tz8o60zXoO6PWVExInwmyOb8HszWD//UrfLmRSvx8vHav4H\nYSfMmyfreGzL97hnzW1k+TLdLk1EiE6baulr5VBnPYc6jnCo4wiHO+vpDw+QOXTHRCfiIdJeltC6\nFH4yYnvTLpr6Wlg39xLy/XlulzNpXo+Xj9d8BMdx+F3DFr6z5bvcvfqTZCoARRKuY6BzKOTqOdR5\nhMMd9aeNyvbgYU5ukKqCBezbC0cP+3F6AuBEe5qKA1nce+OqGa9T4ScjXqzfBMDVs3Sgy3h8Xh+f\nWLaRsBOhrnEbj299hjtX3UKmi3emEUl1vaFeDncc5VDnUNh1HKG1//R78JZkF7O2+HyqApVUFSxg\nQWA+ORnZ0Y3L4HOPbqLV6QeiwffgPYl5/1H4CQBHu45jW/dhihczP3+u2+VMic/r45PL/4intv8j\nW5t28MS2Z7hz5SdcvTWbSKoYDA9S33WMQx31HOw4wuHOI2ctOB3w57Oi9AIWFiwYCbtAZv64x733\nxlUj1/gS0eIb5nr4GWMuBb5urb3mjOffB3wZCAHftdY+5UZ96eKlI9FWX6LW7JspPq+P21Z8lCe3\n/YDtzbv4wn89SuvWFeB4EzqMWiSZhSNhjnefHGnRHe44wtHuE0ScyMg+2b5sTNH5VBUsGPpTSXFW\nUcyrv1RVBBLW2hvN1fAzxnwe+BjQdcbzfuCbwMVAD7DJGPOCtbYh8VWmvq6Bbn578i3KsktYUVbj\ndjnTluHN4PaVN/OF/3yE3qxj+BeHGNi3ZmQY9b03rqKqIuB2mSIJM95cOsdxaOxtGrlGd6ijniOd\nRxmMDI7sk+HNoCpQeVqLrjy3DK8neScMuN3y2wfcAPzgjOdrgH3W2nYAY8wrwJXAvya2vPSw6dhm\nBiMhrlqwPql/mUfzezNo3boCvxnAV9xA5tLfEm6dQ1tXEQ/9uJZv3n2l2yWKJMTwXLphO48e50+e\n2c9Fa/y0RRo41FlPb6h3ZLvX42Vu3pyRkKsqWMC8vIpZP/UpVq6Gn7X2x8aY6nNsKgBGL9zWCbh/\nd+UUFIqEefnoa2T5Mrl87sVulxNfjo+BvReSubgWX1ETvoLoG0B/xMs33txKdcECFhUspLqwitLs\n4oQv1iuSCLuOHsdX1oS3sAlfoBVPZj+DwOvN0e3lOWUsL10aDbrAAhYE5qXFSGm3W35jaQdG90sF\ngNYx9h0RDKorK1avHn6Ttv523r3kahbOLXe7nLhavSRI3d5GBuxFeLJ68Oa3k1vSyZwFAxzprOdg\nx2FeInqtsyArn8WlizCli1hcUs3ikmpyM8e/XZp+32KnczY1sZy3wfAgu5v2U3d8B1tO7CJ77dGR\nbc5AFuGWciLdheQT5Dv33kB+ZvJMa4qn2Rp+u4ElxphioJtol+c3Jvqhxkat8Rarn9sX8eDh0tJ3\npNz5u/fGldFh1J39OP15FGSW8OCN0QvrA+FB6ruOcrD9MG93HOZgxxHeOraNt45tA4bmIuWVR1uG\nBQtYVFjF3Lw5I93CwWAg5c7XTNM5m5rJnLeGniZ2tuxhV/MebOt+Boau1/m9GeT0z6X9RCGR9jKc\nvjzAMzKXrrc9Qi/p+W8yW8LPATDG3ATkW2ufNMZ8FvhPordge9pae9zNAlPRoY4j2OYDrCi9gPLc\noNvlzIixhlFn+vycV1jNeYXVI8+193dwcCgI324/xKHOek50n+S1478d+plMqgKVLCqsYlW/oYQg\nhVkFCf37iAD0hfqwrfvZ1WLZ2byHpr6WkW0VueUsK13KspKlnF+0iEyfP/ohsC/xc+lms5S5yOE4\njqNPlbH5/o5/5rcna/nUmtupKTFulzPrDN9/8FTr8DAnuhtwop/VACjOKmJR4UKqCxayqHAhlfnz\nNbF+DGr5TU0wGKChoYOjXceHWneW/e0HCTthIDrl4IKSxSwrWUpNqaEku/isYxw60Xnah8BUH+1c\nXl4wYbYp/NJUe38HX371a8wNlPNnF/2JBntMUm+ol0Md9TSGT7Lj+F7ebj982q2bvB4vlfnzRsKw\numABwZwynV8UfrHqGuxmd8teDnQfoPbYDjoGTp27hYFKlpUupabEsKhgYcqNxJyuyYTfbOn2lAT7\n76OvEXbCvGfJNXpjjkFORg4XlCzhiuCFXBHcgOM4NPe1nHbtsL7zKIc763n5aHRdxDx/LlXDI0uH\nriHm+nNdW8dMZqdwJMyhziPsbN7DzhbL4Y76kV6GgD+fSyouZFnJUi4oWTLhXVNkYinzrqeW3+QN\nhgf581e/SsSJ8PgHvk5Ha7/bJSWd8Voxg5EQ9Z3Hhq4fHubt9sM0j7omA+APBehtDRDpLsLpy8MJ\nZVCQnccf/+EalsxL7snDY1HL72ytfW0j1+12t+4bmW/n9Xg5v7CamhLD+sUXkjtYkJK/EzNFLT85\npzdP1tE12M21C68mKyMTUPjFk9+bwaLCaLfnsM6BrpEgPNhxmN1Nb5MR7ITgsZF9BoCH9/w/PHs8\nZGdkk5uRQ64/J/p16PucUd9Hn8899bw/hxxfdly6wNQqnRmD4UH2tb/NrmbLzpY9HO8+ObKtNLuY\ni+asZlmJwRQvHrn5c7BYHxpmgsIvzTiOw4v1r+D1eLmy8nK3y0kbgcx8VpYtY2XZMgBu+/qvIKcL\nb147nsw+PBmD4BskMzvCosocekO99IR6OdndMDJsfbKyfVkjYRgNytwJAvTUPj6v7+w7gui2cJNy\nrg8MjuPQ0NvEzuY97Gqx2Nb9I7cN83v9I6Myl5UYynODugSRQAq/NLOv7QBHu46zNrjynKPCJDFq\nqkvYedBDuPdUmAzPvTozYEKRED2hXnoGe0dCsWdw9NceekK99A4/N/R8c28LR8OxteozfZn0F3nI\nWuHHCfkhlIkzkEXnYBbf/H+H+J/vuYjCzACFWQXk+XPVFTfktA8M3hC723Zz749eI7+8jfbBU0v8\nzM2bQ02JYVnpUhYXLtKKIy5S+KWZU2v2bXC5kvR238a1IxPwYfy5VxneDAoyAxRkxt7qCkfC9Ib6\nRgJyODR7zwrQoecHezjU04onsw9v7mn3mycEPLalbuSx1+MdqaswK0BBZgGFmQEKsgpGnivMjH6f\nzKMRB8IDdA120zXYTfdgD90D3XQN9tA9eOrrvqx6spYPQMYgHn8/Hq9DGGjvy2BtxUpqSg3LSpZS\nnF3k9l9Hhij80khzbwtbG3ewIDCf80dN7hZ3JGIdM5/XR35mXky3sDrVinHAP4DH30+gIMI7Ly0j\nI3uQjoEO2gc66ejvpGOgg2PdJzjcWT/uMfP9eZTkFpHnyxsKxoLTvg5/nzWJe0pO53rkaUE2cCrA\nRoJtVMhFv3YzGAlNeFxfIThhH04ok0h3IZGOUiLtZRR4yrn92ismXZ8kTsp0MGu058R+vPdn/OrI\ny3y85iNcOvciQCPwpirVz9tkW6UQvY7cG+qlfaCT9v4OOs7xtWOgk47BTnoH+8Z93WxfFgWjWoxn\nhuSPf3WUfQf7IOwHbxhPxiAFhfCBq+aTn+/EJciG68jz55LnzyPfn0eeP3fo69D3mXnkD23P8+fy\n5PN72XWw47RjjNWNHatU/12bCRrtKSP6Qv28evwNApn5XDhntdvlyCwXS6vU4/FEB8z4c5mbN2fM\n/YLBAPUnmuno76R9YOyQbO/voKGn6dwHKYWcUnAcGB4bMgA8d3js+qJBlsfcvIqRABsOrvzM3NMe\nDwee3xvbW+Ofbrw4pg8M4j6FX5p448Tv6A318d7qd8X8H1vSz0ytrp3lyySYW0owt3Tc/cKR8Glh\nGO1m7eCFN3bj8fdDxiCEM3BCmRDyk+3N4UNX1JwebP48cv25Cft9T0Q3tsSP3gXTQMSJ8FL9Jnwe\nHxvma3qDzH4+r4/i7KKzBojseiPIzr2nr24Wr+7F6ZqpDwwyMzROOQ3satnLyZ5GLpqzmsIszdOS\n5HXfxrUUB7JGHg93L7odfJJ8FH5p4KUjrwBwTaWmN0jyu/fGVRQHskZafCJToW7PFHeiu4GdLXs4\nv7CahQWVbpcjMm3qXpR4UMsvxf1Gk9pFRM6i8EthPYO9vH7idxRnFbG6bLnb5YiIzBoKvxT26vE3\nGAgPcFXluqS+vZSISLwp/FJUxInwcv2r+L1+1s27xO1yRERmFYVfitratJPmvlYurbiQPH+u2+WI\niMwqCr8UNTy9QQNdRETOpvBLQUc6j7G37QAXFC8Z916LIiLpSuGXgl6qH5rUrlafiMg5KfxSTOdA\nF2+erKM8p4xlpUvdLkdEZFZS+KWYV45uJhQJcVXlerwe/fOKiJyL3h1TSCgS4r+Pvkq2L5vLhhar\nFRGRsyn8UkhtwzbaBzq5fN7FZGdku12OiMispfBLIS/Wv4IHD1dX6qa/IiLjUfiliLfbD3Go4wgr\nymooyxl/lWwRkXSn8EsRL2rNPhGRSVP4pYC2/nZqG7cxL68CU3y+2+WIiMx6Cr8U8HL9a0ScCFcv\nWI/H43G7HBGRWU/hl+QGwoO8cux18vy5vGPOhW6XIyKSFBR+Se7Nk7V0D/awft6lZPr8bpcjIpIU\nFH5JzHEcXjzyCl6PlyvnX+52OSIiSUPhl8T2tu3nWPcJ1gZXUpxd5HY5IiJJQ+GXxH6tNftERKZE\n4ZekGnua2d60i6rAAhYVLHS7HBGRpKLwS1K/OboJB0fTG0REpkDhl4T6Qn28duxNCjMDXFi+yu1y\nRESSjsIvCb1+/Hf0hfu4Yv7lZHgz3C5HRCTpKPySTMSJ8Jv6TWR4fGyYf5nb5YiIJCU1G5LIA8/W\nsqd1D5lLm8jvXUQgM9/tkkREkpJafknigWdr2XmwFV/FIQCa9s/lc49u4tCJTpcrExFJPgq/JLHr\nYCue7C58hc2EO4pxegpo7ezn4R9tdbs0EZGko/BLIr6yYwCEGzSvT0RkOly75meM8QKPAauAfuB2\na+3+Uds/CHwRcIDvWmv/3pVCZ4kLqos4EDiBE/YSbgsCUBzI4t4bNdVBRCRWbrb8rgcyrbXrgD8D\nHjxj+zeBa4H1wOeMMYUJrm9W+eh18/Dm9BBpD0Ikg+JAFg/es56qioDbpYmIJB03w2898AsAa+1m\n4OIztg8CRUAO4CHaAkxbtY3bAMjurVSLT0Rkmtyc6lAAdIx6HDbGeK21kaHHDwK/A7qBH1lrO848\nQDqpa9hGhsfH1276ADkZ2W6XIyKS1NwMvw5gdJ/dSPAZYxYCnwKqgB7gH40xH7LW/ut4BwwGU7ML\n8FjnSY51n+DCeStZODcY9+On6nmbaTpvsdM5mxqdt/hzM/w2Ae8DnjPGXAaMHrOfDYSBfmttxBjT\nQLQLdFyNjak55+3XB18HYHlhTdz/jsFgIGXP20zSeYudztnU6LzNjJjCzxiTBVwOVANlRAPqJHAY\neN1aG4rhcM8D1xpjNg09vtUYcxOQb6190hjzDPCqMaYP2Ad8P5ZaU0ld4za8Hi+rypa5XYqISEqY\nMPyMMR7gOqLdkFcD/jF27TLG/Bp4ylr7s4mOa611gLvOfHrU9m8B35roOKmuubeFw51HqSkx5Ppz\n3S5HRCQljBt+xpj3AA8B5wGvAd8AtgEHiF6z8wKlQCVwGbABeMEYsxv4orX232au9PQwPMpzTXCF\ny5WIiKSOMcPPGPMj4Aqi4fdda+3xCY717NDPLQE+DjxljPm4tfaGeBWbjuoatuPBw2qFn4hI3IzX\n8qsFPm6t7Y7lgNbavcCXjTHfAP5kOsWlu7b+dt7uOMSSovO0goOISByNGX7W2r8Z/dgYsxzYNWoe\n3riG5uX99fTKS291jdsBWFO+0uVKRERSSyx3ePk18NWZKkTOVteg630iIjMhlvDLA96eqULkdJ0D\nXexre5vzCqsoykrr25qKiMRdLOH3f4DPGmPeMVPFyClbGrfj4LAmqC5PEZF4i2WS+0XAPGCzMaYH\naCY6yX2YB3CstefFsb60VasuTxGRGRNL+OUAbxINubGk9coL8dI92INt28/CQCWlOSVulyMiknIm\nHX7W2qtnsA4ZZWvTTiJOhLXq8hQRmRFjXvMzxlw+3YMbYzZM9xjpaGSUZ7m6PEVEZsJ4A17+rzHm\nBWPMpbEe1Bjze8aY/wL+aeqlpafeUB+7Wyzz8iooz43/8kUiIjJ+t2cN8JfAfxtjDgE/Af4d2Gqt\nbR69ozGmnOi9Pa8APgzMBb4D6NZmMdretIuQE2atJraLiMyY8e7w0g38qTHmO8CngduBzwIYY7qI\n3tjaB5RwaqWHNqJLD/0fa+3hmSs7ddWN3Mha4SciMlMmHPBirT0A/C9jzBeJrtqwgegqD6VAhOh6\nfoeI3gHmNWtteKxjyfj6wwPsaN7DnNwgc/PmuF2OiEjKimW0Zy/wy6E/MgN2NO9mMDLI2uBKPJ7x\nZpSIiMh0TDr8jDFvM/48PgfoBxqAN4AHrbUnp1deejk1ylNdniIiMymW25v9CigAqomG3BZgM9Ay\n9NxcoIlod+jngC3GmIVxrDWlDYYH2d68i9LsEirz57ldjohISosl/N4ienPr6621NdbaD1prb7LW\nXgysA0LAM9balcAaoi3Bvxn7cDLarhZLf3iANeUr1OUpIjLDYgm/zwEPW2tfOHODtfZ14GHgC0OP\ntwGPAe+KR5HpYHjtvrXBVS5XIiKS+mIJvzlA/TjbG4HKUY+PA0VTKSrdhCIhtjbtpCirkKqCyol/\nQEREpiWW8NsB3GKMyTpzgzEmE7gZ2D3q6QuJToGQCdjW/fSGelkTXIHXE8s/iYiITEUsqzr8JfBT\nogNZ/h7YS3Tgy1LgNmA18BGAoYnxtwN/Ec9iU9Wp5Ys0ylNEJBFimef3H8aY64kuavvNMzYfAT5s\nrf2RMSZINAx/CDwQt0pTVDgSZmvTDgKZ+ZxfVO12OSIiaSGWlh/W2p8BPzPGrAKWEL2t2QHgt9ba\n4TmAzUAGgpNuAAAMM0lEQVS+tXYgrpWmqP3tb9M12M2G+Zepy1NEJEFiCr9h1tqtwNYxtkUABd8k\n1TYMj/JUl6eISKKoqeGiiBNhS+M28vy5LCk6z+1yRETShsLPRW+3H6Z9oJNVZcvxeX1ulyMikjYU\nfi46tXyRVmwXEUkkhZ9LHMehtmEb2b5slpYscbscEZG0ovBzyeHOelr721hZtgy/d0rjjkREZIoU\nfi4Znti+tlxdniIiiabwc4HjONQ1biPTl0lNyVK3yxERSTsKPxcc7TpOY28zy0svINPnd7scEZG0\no/BzwfAoT01sFxFxh8LPBbWN2/F7M1heeoHbpYiIpCWFX4Kd6D7Jie6T1JQsJTvjrNWhREQkARR+\nCTZyL89ydXmKiLhF4ZdgdY3b8Hl8rCitcbsUEZG0pfBLoKbeZuq7jrG0ZDG5/hy3yxERSVsKvwQa\nmdiuUZ4iIq5S+CVQbeM2vB4vq8qWu12KiEhaU/glSGtfG4c6jrCk6DzyM/PcLkdEJK0p/BKkrjE6\nynONujxFRFyn8EuQ2oatePCwWmv3iYi4TuGXAO39nRxoP8R5hdUUZgXcLkdEJO0p/BJgS+N2HBxN\nbBcRmSVcW0XVGOMFHgNWAf3A7dba/aO2vwN4EPAAR4GPW2sH3Kh1uoZvZL1GXZ4iIrOCmy2/64FM\na+064M+IBh0AxhgP8ARwi7X2CuBXwCJXqpymroFu9rYdoLpgIcXZRW6XIyIiuBt+64FfAFhrNwMX\nj9pmgGbgs8aYl4Aia+2eRBcYD1ubdhBxImr1iYjMIm6GXwHQMepxeKgrFKAMWAc8ArwLeKcx5poE\n1xcXtcNr9+l6n4jIrOHaNT+iwTd66KPXWhsZ+r4Z2Dfc2jPG/IJoy/DF8Q4YDM6ukZTdAz3sad1H\ndVElNQur3S5nTLPtvCULnbfY6ZxNjc5b/LkZfpuA9wHPGWMuA7aO2nYAyDfGnD80COYK4KmJDtjY\n2DkjhU7V5uO/IxwJs7Jk+ayrbVgwGJi1tc1mOm+x0zmbGp23meFm+D0PXGuM2TT0+FZjzE1AvrX2\nSWPMbcA/DQ1+2WSt/Q/XKp0i3dVFRGR2ci38rLUOcNeZT4/a/iJwaUKLiqO+UB87W/ZQkTeHirxy\nt8sREZFRNMl9huxo3k0oEmKtRnmKiMw6Cr8ZUjvU5bm2fJXLlYiIyJkUfjNgIDzIjubdBHNKmZdX\n4XY5IiJyBoXfDNjVsoeB8ABrgivxeDxulyMiImdQ+M2A2gZNbBcRmc0UfnE2GAmxrWkXJdnFLAxU\nul2OiIicg8Ivzva07KUv3Mea4Ap1eYqIzFIKvzjTxHYRkdlP4RdH4UiYrY07KMwMsKhwodvliIjI\nGBR+cbS37QDdoR5WB1fi9ejUiojMVnqHjqNTyxfpri4iIrOZwi9OIk6ELQ3byffncX5hUi46LyKS\nNhR+cbK/7SCdg12sDi7H5/W5XY6IiIxD4RcndUNdnhrlKSIy+yn84iDiRKhr3E5ORg6m+Hy3yxER\nkQko/OLgUMcR2vrbWVW2jAyvm+sDi4jIZCj84uDUKE91eYqIJAOF3zQ5jkNdw3ayfJlcULzE7XJE\nRGQSFH7TVN91jOa+FlaU1uD3+d0uR0REJkHhN02nli/Siu0iIslC4TcNjuNQ17gNv9fPstKlbpcj\nIiKTpPCbhuPdJznZ08jy0qVk+TLdLkdERCZJ4TcNtZrYLiKSlBR+01DXsI0Mj48VZTVulyIiIjFQ\n+E1RQ08jx7pPcEGJIScj2+1yREQkBgq/KaprGFqxXRPbRUSSjsJvimobt+L1eFlVtsztUkREJEYK\nvylo7m3hcOdRlhYvJs+f63Y5IiISI4XfFNQ1DnV5BrViu4hIMlL4TUFd4zY8eFit8BMRSUoKvxi1\n9bdzoP0Qi4sWEcjMd7scERGZAoVfjEa6PDXKU0QkaSn8YlTXMHxXF3V5iogkK4VfDDoHutjX9jaL\nCqooyip0uxwREZkihV8MtjRux8HRiu0iIklO4RcDTXEQEUkNCr9J6h7sYU/rPhYG5lOaU+J2OSIi\nMg0Kv0na2rSTiBNhbVArtouIJDuF3ySNjPIsV5eniEiyU/hNQm+oj90tlnl5FZTnBt0uR0REpknh\nNwk7mnYRcsKa2C4ikiIUfpNQ2xjt8lwbVPiJiKQChd8E+sMD7Gjew5zcIHPz5rhdjoiIxIHCbwI7\nm/cwGBlkTXAlHo/H7XJERCQOFH4TqGvUKE8RkVSj8BvHYHiQbU07Kc0uYUH+fLfLERGROFH4jWN3\n6176wwOsKV+hLk8RkRSS4dYLG2O8wGPAKqAfuN1au/8c+z0BNFtrv5DgEqlt0ChPEZFU5GbL73og\n01q7Dvgz4MEzdzDG3AmsAJwE10YoEmJr006KsgqpKliQ6JcXEZEZ5Gb4rQd+AWCt3QxcPHqjMWYd\ncAnwOJDwPkfbup/eUC9rgivwetQ7LCKSSlzr9gQKgI5Rj8PGGK+1NmKMmQt8Bfgg8JHJHjAYDMSt\nuN0H9wBw9ZJL43rc2SjV/34zRectdjpnU6PzFn9uhl8HMPpf1GutjQx9/yGgDPg5UAHkGmN2WWv/\nYbwDNjZ2xqWwcCTM5iO1BPz5lFIet+PORsFgIKX/fjNF5y12OmdTo/M2M9wMv03A+4DnjDGXAVuH\nN1hrHwEeATDGfAK4YKLgi5cHnq1lT/M+Mmu6CfQsVpeniEgKcvOd/Xmgzxiziehgl/9ljLnJGHPH\nOfZNyICXB56tZefBVrwlJwFoOlzE5x7dxKET+tQlIpJKXGv5WWsd4K4znz7Hfs8kpiLYdbAVcPAV\nn8QJ+Yl0ltDq9PPwj7by4D3rE1WGiIjMMPXpnckXAn8/4eYKcHR6RERSkd7dR6mpLoawn/5tGxg8\ncgEAxYEs7r1xlcuViYhIPCn8Rrlv41qKA1k4ffkQ8VEcyOLBe9ZTVaFhxiIiqUThd4Z7b1xFcSBL\nLT4RkRTm5lSHWamqIqDBLSIiKU4tPxERSTsKPxERSTsKPxERSTsKPxERSTsKPxERSTsKPxERSTsK\nPxERSTsKPxERSTsKPxERSTsKPxERSTsKPxERSTsKPxERSTsKPxERSTsKPxERSTsKPxERSTsKPxER\nSTsKPxERSTsKPxEREREREREREREREREREREREREREREREREREREREUkYj9sFTJcxxgs8BqwC+oHb\nrbX73a0qeRhjLgW+bq29xu1aZjtjjB/4LlAFZAF/Y639qbtVzX7GGB/wJGAAB/hja+0Od6tKDsaY\ncuB3wDuttdbtepKBMeYtoH3o4QFr7W3n2i8jcSXNmOuBTGvtuqE38geHnpMJGGM+D3wM6HK7liTx\nUaDRWnuzMaYYqAMUfhO7DohYazcYY64C/hb9H53Q0Ietx4Fut2tJFsaYbIDJfJhPhdubrQd+AWCt\n3Qxc7G45SWUfcAMp0AOQIM8BXxn63guEXKwlaVhrfwLcOfSwGmh1r5qk8g3gO8BxtwtJIquBXGPM\nfxpjfjXUIDqnVAi/AqBj1OPwUFeoTMBa+2P0Bj5p1tpua22XMSZANAi/5HZNycJaGzbGfB94GPgn\nl8uZ9YwxtxDtZfivoaf0AXVyuoFvWGv/APhj4Idj5UEqhEQHEBj12GutjbhVjKQ2Y8wC4NfAP1hr\nn3W7nmRirb2F6HW/J40xOS6XM9vdClxrjHkRWAM8Y4yZ43JNycACPwSw1u4FmoG559oxFa75bQLe\nBzxnjLkM2OpyPZKiht58/gu421r7otv1JAtjzM1ApbX2a0AvEBn6I2Ow1l41/P1QAN5prT3pYknJ\n4laigx/vMcbMI9ozeM5u41QIv+eJfkLaNPT4VjeLSVKO2wUkiS8ChcBXjDHD1/7eY63tc7GmZPCv\nwPeNMb8B/MBnrLX9Ltckqelp4HvGmJeHHt+qnkARERERERERERERERERERERERERERERERERERER\nERERERERERERmWFaJkMkBRhjDhK96bYP+COgCVhtrW1xsy6R2SoVbmwtItGbk98E7ATuBSoUfCJj\nU/iJpAYPkA18wFp7wu1iRGa7VFjMVkSi9in4RCZH4SeSOhrcLkAkWSj8RFJH2O0CRJKFwk9ERNKO\nwk9ERNKOwk9ERNKOwk8kNThuFyAiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIgL/\nHxwpL5w8i/LUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xae11d0ec>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(rad_py, g_r_py, 'o', label='python')\n",
    "plt.plot(rad_for, g_r_for, 'x-', label='fortran')\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.xlabel('r', fontsize=18)\n",
    "plt.ylabel('g(r)', fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voila!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Time Tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us perform some timing tests. Let's start with `npart = 10` particles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(12311)\n",
    "npart = 10\n",
    "Lbox  = 10.\n",
    "r     = np.random.uniform(0., Lbox, (npart,3))\n",
    "dr    = 0.1\n",
    "num_increments = 50 # Lbox/(2*dr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 1.97 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit rad_py, g_r_py = py_rdf(r, Lbox, dr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 5.25 times longer than the fastest. This could mean that an intermediate result is being cached \n",
      "100000 loops, best of 3: 4.54 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit rad_for, g_r_for = grFor.grself(num_increments, r, Lbox, npart)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Fortran enabled code ran nearly **500x** faster! One can repeat this with different number of particles. Note that the basic algorithm is $O(N^2)$, so expect to wait a longer time for large $N$.\n",
    "\n",
    "I did a few tests with up to 10000 particles. Here were my timings in seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = np.array([10, 100, 1000, 10000])\n",
    "pyt = np.array([1.92e-3, 19e-3, 366e-3, 22.2])\n",
    "fort = np.array([3.97e-6, 163e-6, 20.5e-3, 2.01])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0xadf26d6c>"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEgCAYAAACXa1X+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd0m9d5+PEvwD3AIW5S3NSl9qIkalCy9rCtYTt1bNdp\nRuM0jlPHTR13nWb+GruJnTRpM932pHVrO3btxFOyZEm29l7U4pVEcVMkSHGTINb7+wMgRUmkRIBY\nJO/nHB0TwDse8jXw4L3PHaAoiqIoiqIoiqIoiqIoiqIoiqIoiqIoiqIoiqIoiqIoijIonb8DcIcQ\nYhXwWSAS+JGU8oyfQ1IURVGc9P4OwE0RUsqvAC8Ca/0djKIoijIGCCGihBD/KYRI9HcsiqIoyg3B\n/g7gVkKIYuAFKeUKIYQe+CUwE+gFviylvOJMJj8Cvi2lbPJjuIqiKMotAqopTAjxHPAyEOZ8agsQ\nKqVcDPwt8JLz+ZeAFOB5IcRDPg9UURRFGVKg3bFcBh4EXnE+LgG2AUgpDwsh5jl//rx/wlMURVHu\nJqDuWKSUbwPWAU8ZgPYBj23O5jFFURQlQAXaHcut2nEklz56KaXdlQPY7XZNpxuVvaoVRVH8RjeC\nD85ATyz7gY3Am0KIhYDL41V0Oh1GY4fHA3NVUpLBo3G4czxX9hnOtnfbZqjXXXne0383dwTCtXN1\nv5FeP1dfC9Rr5404vH39AuW9NxIB91VeCJEDvCqlXCyE0HGjVxjAF6WU0pXjaZqmeThERVGUMW8k\ndywBl1g8TdM0TX1rcn2fQPnWFAjfegPh2rm6n7pj8V4c4+WOJTk5xu38oArhiqIoikepxKIoiqJ4\n1LhoCvN3DIqiKKPNWO4V5hGqndf1fQKlnTcQ2ukD4dq5up+qsXgvjvFSYxkJ1RSmKIqieNS4uGMZ\nqRdfP8mFihYApuTE8+wjc3xy3lOnTmAwxJCfX8CmTet4992PfHJeRVGUkVB3LHfx4usnOV/RggZo\nwPmKFv76F/upvOb9W/wPPniXpiYjAGryAEVRRosx/3F1t+L9f753jv2na4d8vbGlZ9Dn9TpIjIsY\n9LUlszL40sZpQx7z7bffZu/evbS0tNDS0sJ9993Hjh07ePPNNwF45pln+NKXvsSTTz5JYmIiv/71\nr9m8eTNLly6lvr6euLg4fv7zn9Pd3c23vvUturq6sFqtPPPMMyxcuJCNGzdSXFxMWVkZAL/61a+I\njo6+059BURTlJqp4fxd3Kkr1dJux2dzrODbUfj3d5jsWMjs7e+npMfPjH/8rzc1NfOUrXyAlJZWj\nR88wYcIEKiqqSEvLZcGCRaxevY7g4Gi6urr4whe+SmpqKn/5l3/B+fPneeutd5g9ex6f+cwjNDUZ\nefLJL/Pmm+/Q0dHJkiUr+Yu/+Abf//4/8v77H7Fq1VpVvHfTaCv+DndbVbz37fHGU/F+XCSWO3l4\nZQEPrywY8vW+prCB4g1hPP3QTLJTDUPsdXdFRfMBSEhIxGCIYcuWz/Dhh++SkpLG+vX33rZ9TEws\nqampAEyYkIDJZKKqqoJ16zYAkJiYRFRUFC0t1wEQohCA5OQUzGaz23EqiqK4StVY7uLZR+YQbwjr\nfxxvCOOlp5aMKKkAXLx4HoDr15sxmXpYtmw5R44cZs+eT1i71pFYdDodNpvN+fPtx8jOzuX06ZMA\nGI2NdHZ2EBMT27+voiiKP4z7O5bhePqhmfz8rTP9P3tCTU013/jG1+ju7uTZZ/+O8PBwZs+eS1tb\nKwaDI2lNnTqd3/zmF6SnZ3BrOUyn0/G5z32R55//Pp98soveXhPPPfcPBAUF3batoijKcPX1gr3/\nm3+0v/+TLW7dfKjEMgzZqQZeemqJR4+5ZMkyHn308Zue0zQ7Gzdu6X+8efODbN78IADvvLOt//nv\nfe+H/W2izz//4m3HfvPNd/p//upXv+7RuBVFGbtuafpXk1CONre2VH3zm1+no6ODuXPn+ScgRVHG\nvQu31JPdpe5Y/GDDhvtve+4nP/k3P0SiKIoC1Q0d/H77RTw1seKYb4xXk1AqiqLcTtM0zlxq4o97\nrnDsQgMAoSF6zJYbq7+//5MtbuWIcXHHovrSu75PoPSlD4SxEIFw7VzdT41j8V4co30cS119G0cu\nNPDRkWpqjJ0ATM2dwIrZGcyZlMi3fnWAlo7eu8Z5J6M6sQghVgKPSimf8HcsiqIogayzx8Kuj8t4\nb285bZ1m9DodC6Yks3Z+FsWzMvoTUV8v2JaO3qGnJLmLUZtYhBD5wGwg3N+xKIqiBKpr17vZcbSa\n/aX1mK12IsKCWL8gi1VFE0mIvf3js68XbHJyzER3zzlqE4uU8grwEyHEK948z1O7nhvytV+s/JFb\nx7TZbDzzzNewWq38+Mc/u+M8Xu3t7Rw+fIA1a9a7dS5FUcYfTdO4WNnC9qPVnLrcBEBibDhblhcw\nJ28CEWHe/egPyMQihCgGXpBSrhBC6IFfAjOBXuDLzqQyahmNRrq7u/mP/7h7Trx8WbJv3x6VWBRF\nuSurzc7RC43sPFlLeW0bAPkZMaybn8UckUhqSqxP6l4Bl1iEEM8BjwOdzqe2AKFSysXOhPOS8zmP\nePvy+5xsLHVr33888Pygz89JnsGDBbd3Ke7z4os/pKamih//+IfOJNOFzWbliSe+xty58/jc5x4m\nKyub4OAQ2tvbuHz5Eu+++wdKS0/T3t5Ge3s7//7vv+WFF35IY2Mjzc1NlJQs44knnuSf/um7hIaG\nUl9fT3NzE//wD99BiMlu/X6KoowOXSYLn5ysZefxGlo7zeh1MG9yMmvnZ1KQEevzeAIusQCXgQeB\nvq/zJcA2ACnlYSHETSMIpZSf8214I/fss3/Hd77z90RGRrFgQd5tsxObTCa+8IUnmDRJcPLkcd55\n5202bXqAs2fPUFS0gIcffpSurnamT5/B/fdvobe3l4ceuo8nnngSnU5Hamo63/rW3/Pee3/k3Xf/\nwLPP/p2/f2VFUbygoaWbj4/WsLe0DrPFTnhoEGvnZ/Lw2snonfMM+kPAJRYp5dtCiJwBTxmA9gGP\nbUIIvZTSzjAlJQ09YeRfJD0KPDrk6w///skhX/v15h8ON4Sb4ujtbSMkJIj6+mo++9mHSEoykJRk\nIDbWgF5vRq/XUVQ0nbCwMGJjIwgLCyYpyUB4eAgzZkwmKclAZ6eOq1cv8c///D2io6OxWCz928yf\nP5ukJAMFBdlcunT+pt//Tn+LoeIdyTZDve7K867E7C2ejsHd4/ny+rn6WqBeOxh91+9O22maxrny\nZv746WUOn7uGpkFSfASbluaxZkE2UREhbh3bk3+jgEssg2jHkVz6uJRUwHvjWFw57sB+5devd2Gx\n2EhPz+STT/aRmDgRo7GR1tY2LJYg7HaN5uYuQkLMdHQ41m4xGjswmSy0t5swGjvYuvUPBAeH8zd/\n8y1qaqp54403+rdpa+vBaOygra0Hk8nSf141jsU9o20cxHC3VeNYfHu8kY5jsdrsHCtrZPuRaiqc\nK9jmphlYtyCLosIkgvR6ujtNdHeaXHqPaZpGaH0F146fIWHjZpd/r8GMhsSyH9gIvCmEWAic8XM8\nHqHT6Xj88bvPTpyRMZHy8su88cZr/fsBLF68mN///q8oK7tAamoahYVTBixjrLvpv4qijF7dJguf\nnq7j42M1tHT0ogMWzUhj+aw0CjJi3X6fa5pGz8ULNL/7R3ouSdDpiFm0mJDEpBHHHJCfPM6msFed\nBXsdN3qFAXxRSimHeyw1pYuiKKPRteYu3ttbzo4jlfT02ggPDWJNcTYbS/JIS4wa0bFbz5RS/drv\naT9/AYD4+fPIeuRhogvy+7cZydLEAZlYPEnTNE3djru+j2oK814MqinMt0bT9dM0jeYuK7/fcZET\n0oimORYXXD1vIvfMSicyPGRYx7rbe8z45u9p+WgrUbNmk7BxC5nzZ9y2fXJyjFrzXlEUZbSy2e0c\nLzOy/Wg15XWOvkrZKQbWLchk3uRkgoM8u8JJ/PoNGOYXE56T49Hj9lGJRVEUxU+6TVb2nnHUT5rb\nTeiA4mmpLJ+VhsiMG1GdVNM0TFevEp6be9txgg0xBBtiRhj90MZFU5i/Y1AURRmo8Xo37+4tZ/vh\nSnp6rYSFBrF6fhabluaRnjT0FE/DoWkabWdKqX79DdrPX2DaD75L3MwZLh9nJDWWcXHHotp5Xd9H\n1Vi8F4OqsfhWIF2/Q6dr2H6kmmNljWgaxEWHcu/CPO6ZnUF0RAg4l9py59rd1ssLiJo1m057MJa7\nDDnw9N9oXCQWRVEUf7HbNU5II7tfP8WFiusAZCVHs3ZBJgumpHisftJx6ADX/uNlwNHLy7Dufq/V\nUO5GJRZFURQv6Om1su9MPTuOVdPUZgJgVn4CaxdkMTlrZPWTwUTPnYfh7Fni16wbtJeXL6nEoiiK\n4kHNbSZ2Hq/h09N19PRaCQnWs3xOBp9dW0iYB3KJpmkMVjrWh4WR9sRfjPwEHqCK94qiKB5wqbqF\nP356hX2n67DbNeINYdxXksv6hTnERoeN+PgDi/Jp920gsWSJB6Iemire34UqILq+jyreey8GVbz3\nLW9eP7td49TlJrYfqULWONY/mZgUxboFWSyYkkJIsB5zjxljj9ntucIGK8pHF+SjFc4c6hAeee+N\nxLhILIqiKJ5kMlvZX3qNHUeraWztAWBGXgJrF2QyNTveY/UTa2sL9b/51U29vIYaKR9IVGJRFEUZ\npuvtJj44XMXWAxV091oJDtKzbFY6a+ZnkjHC+bsGE2SIwdra2p9Q/NXLy1UqsSiKotxF5bUOth+t\n4siFRmx2jZjIELaU5LJ8TgYxUaFeO68uKIisb3+PoIgIr53DG1RiURRFGYRd0zhzuZntR6u4WNUK\nQHpiFA+tnMT0rFhCgoM8cp6+Gord1EP0nKLbXvdlUnlq13P9P//J61/V3nzk12616anEoiiKMkCv\nxcaB0nq2H62mocVRP5mWO4F18zOZljuB5OQYj9Q3bi3KB8fHEzVj1oiPGwhUYlEURQFaO3vZebyG\nT07W0mWyEhyko2RmGmvnZTIxeWTzdw006NQrM2eRsGkLuuCx8ZGsxrEoijKuXa1r44+fXmHPyRqs\nNo2YqFDuXZzLvUtyiDeEe/x8mqZx5lt/S+ely8TPLyLzsw9jmFTg8fO4wq7ZOVxzkp8e+PebnldN\nYXcQCN3yAmEshBrH4p5AuHau7qfGsdw5Drumcba8mY+OVHOhsgWAtIRI1szPZPG0VEJDgrCaLBhN\nlmEdz9U44v/kMSbodYTn5GICTAOO58v3nl2zc7KxlB3Vu6hur3fxNxraqEwsQojFwFecD78hpWzz\nZzyKoowOZouNA+cc40/qm7sBmJIdz7oFmUzPS0Dvwfm7NE3Der2ZkITE216LyMvz2Hnc4UgoZ9ha\nsZP6rgZ0Oh3FqUUcvnbcI8cfdmIRQkwE1gElQC6QCNiABqAK2Alsl1I2eySyO3sCR2IpBj4L/NYH\n51QUZZRq6+xl14ladp+spbPHQpBex5LpqayZn0lWisGj5xpYQ+mtqSb3hRcJivL8GBd32O12jjWc\nYmvFTq51NaDX6SlOLeKxuZsINkX4LrEIIUqAvwI2AUFAN1ABdAB6IA9YAnwJsAsh/gi8KKU85JEI\nBxckpTQLIeqBlV48j6Ioo1iNsZNXd17mkxPVWG0aUeHB3L84m5VzJxLngfm7BhqqKG839fg9sdg1\nOycaTrP92G5q26+h1+lZmDqPdTkrSY5MJMlgwGjq4Bcrf9S/j1fWvBdCpAH/AjwIfAx8FfgUuCKl\n1G7ZNgiYgeNu5lHggBDiXeApKWWtKwEJIYqBF6SUK4QQeuCXwEygF/iylPIK0C2ECAXSgWuuHF9R\nlLFN0zTOXb3OR0erOXfVsf5JSnwEa+dnsnh6GmGhnhl/cqumt96kZduHwI1eXuE5uV4513DZNTvH\nG06ztWInDd2N6HV6FqXNZ132SpIiE7x23jvdsZQCbwD5UsqqOx1ESmkDTjn//ZsQogD4W+AMMOzo\nhRDPAY8Dnc6ntgChUsrFzoTzkvO53wK/ccYfGPNEK4riVxarjYPnGthxtJrapi4AJmfF8SerC8lO\nivRo/WQw0XPnYa6rDYyEYrdz5NoJtlXspKHbiF6nZ3HafB6duxF9j+d7ut3qTollgZSy3J2DSikv\nA18WQjzv4q6XcdwhveJ8XAJscx7zsBBinvPnE8AX3YlNUZSxpb3LzO6Ttew6UUNHt6N+smhaCmvn\nZ5GdavBZ77SIvDwynv4rr5/nTmx2G8cbT7P96C7qOxqdCWUB63JWkhgxgaRoA8Ye7/8tAm4cixAi\nB3hNSrlICPEy8JaUcpvztUogV0ppH+7x1DgWRRmbqq618+7ecnYdq8ZitRMVEcKGRTncX5JLQqzn\np0HRNI220rPUvPkWk55+irCkJI+fw102u439Vcd469yH1Hc2EqTTszx3MQ9MXU9ylHtNXj5bj0UI\nEQNMk1IedD5eBjwNWIFfSin3uBvIENqBgV029K4klT6B2pfe18dT41jcEwjXztX9xuo4Fk3TOF/Z\nwvYj1ZSWOzqgJsdFsGZ+JktmpBIeGozdbL3pvCONY7CifNXu/cStWOXScdxdj+VO29jsNmcvr48x\n9jSj1+lZkl7Mo3M3ousOhW4wdt/9b+G39ViEEFOBT3B0L54hhMgDtuO467EADwohNkgpd3osOtgP\nbATeFEIsxFGzURRlnLFY7Rw+38D2o1XUGB31EzExlrULsphdkIhe753Gl96aahr/95Wb1kMp+LPH\n6IlN9sr5hstmt3Go/hjbKnZi7GkmSBdESXoxa7NXkhART1KU4aaE4mvDvhpCiLdxdCv+gpRyqxDi\nhzgK9CU4iva7gS4p5Yi6/zqbwl51Fux13OgVBvBFKaV05XiqKUxRRq+2zl62Hazg/f1Xae3oRa/X\nUTIrnc3L8hFZ8V4/f099PSef+gZxc+eQ9cjDRBfke/2cd2K129hbcZi3L2yjodNIkD6IlbmLeWDK\nehKjJnj0XCNpCnMlsTQBL0kpn3c+PguESCkLnY+/DjwvpfTsaKMR0jRN83dTCgRGc4pqCnNPIFw7\nV/cb7U1h9c1d7DhWw4HSesxWOxFhwdwzO53VRROZEONar6aRXj9LczMhCTfqFN6+foNtZ7PbOHzt\nBB9V7KTJdJ1gfTCL0uazNns5E8JvT7CeeO95ZRzLICKARgAhRDYwFfjXAa9rOEbiK4qiuEzTNC5W\ntbL9SBWnrzjqJ4mx4ayZl0nJzDQiwrwzA1VfDSU4Lo7QtPTbXh+YVHzNkVCOs61iF82m6wTrgliW\nsYhH5t6P1hXit7juxpUrdQVHs9d/AJ93PvcOgHMg42eASx6NTlGUMc9qs3PkQgPbj1RT1egYwlaQ\nEcva+ZnMFUleq5/cWpSPLppH+pNf98q5XGWz2zh07RgfVeyi2dTiTCiLWZu9nPjwOBIjDRi7/N8S\nMxRXmsK+iqPecRaYAlwAZuG4c/kf589/JqX8Hy/E6TZVY1GUwNTRbXbUT/aVc729F70OFs9MZ/M9\n+UzO9my9YCBN02g7U0r162/Qfv4CAPHz5wVGDcVm5ZOKQ/zhwjaMXc2E6INZlVfC5ilrSYj0fk1p\nIJ/UWACEEH8K/ClQDXxfSlkrhJgCvAn8SEr53+4G4i2qxuLePqrG4r0YxnuNpeF6N9uPVbO/tB6z\nxU54aBDLZjnqJ4lxnh9/cmscts5Oyp/7JprZTNSs2SRs3EJ4To7bx/PEfla7lUP1x/iocjfXTS0E\n64MpSS9mTfZy4sJiXY5hNNVYkFL+L/C/tzx3AZjubgCKoox9mqYhq1vZfrSaU5ea0ICEmHDWLJ3I\n0lnpXqufDCYoOpqUxz9PaHqGSwnFG6x2KwfrHU1eLb2thOiDuXfSCpYkLx40oYwWd5qE8u+Bn0op\ne9w5sBDCgGOtlP/nbnCKoow+L75+kgsVLaCDyVnxLJ2ZxkdHq6m85vhGnJsWw7oFmRQVJhGk13st\nDk3TsPf0cPMYa4eYxUu8dt7hsNitHKo/ykcVu/sTyorMEtZkLadgYobf79JHashbHee4lcXAT4H/\nlVLWDOeAQoh8HMX9rwH7pJRbPBGou1SNRVF85x9/fYBTl4y3Pa8DFs1MY8uyAibnxDOC5vu7GlhD\nQadj+j9936vnc4XFZmH31QP84fxHNPe0EBIUwtr8ZWyavIb4iMC6Q/HKlC5SygeFEPfhmDr/eSHE\nCeADHKPfr+KYbkWPY/biTGAhsBQoAiTwhJTyD+4G5kmBkP0DoZ1e1VjcEwjXztX9/FFjOScbBk0q\nAIaoUL587xQAmpo6B91mpAZdD2XWbOy9vVzvuH2JYXe5c/0sdiul7Wd469xWWnvbCNGHsDJzKauz\nlhMbZsDaCcbOjmEf3xfvvZG4Y8OmlPIDIcRWHNOqPIVjpP1Qnae7cIy+3yKlfNdjESqKEpA0TaO2\nqYvjZUZOX2mmor59yG2DvNRleKD6X/0bnSccKyAOLMoHhYeDBxOLKyw2Cwfqj7K9cnd/QlmVuYzV\n2fcQExpQY8k96q4VM+ekj+8A7wghInEsB5yH407FjmPusErgsJSy14uxKoriZ5qmIata+PhwBSfK\njDS0OEqwwUF6ZuUnYGztoc65lnyfeEMYTz80c7DDeVR4fgGazeZyLy9vsNgs7K87wvbK3bSZ2wnV\nh7CxcDWLkxaN6YTSx9VeYd047kp2eyccRVECjd2ucammleNlRo5LIy0dju+PYSFBzJuczLzCJFYs\nyKarwwTAX/9if/828YYwXnrKN4Xy+LXrmbBug0/ONRSzzcL+usPsqPykP6GszrqH1Vn3kJeR5vdm\nXV/xXR8/RVFGDYvVTml5M8fLjJy8ZKSj29GUFBUezMp5mUzLjmNazgRCQxzL/EaGh/QnlqcfmsnP\n3zqDXq/j6w/M8FhMfTWU9oP7SfnCn6O7pUeZPwv0NxLKbtrMHYQGhbImazmrspZhCI32W1z+ohKL\noigA9FpsnC1v5rg0UnqlmS6TFYCYqFCWz8mgSCRRmBVHWmrsHb95Z6caeOmpJR4rCA9WlI9ZspTI\nwskjPvZImW0W9tUdYkflJ7SrhNJPJRZFGce6TVZOX2niRJmR0vJmzFbHOnrJ8REsmZFGUWES+emx\nXpuv6256Ll2i6e03b+rlFQg1FLPNzL7aQ+yo+pR2cwdhQaGszV7BqsxlRIdG+TW2QKASi6KMM+3d\nZk5daqL06nVOSSM2u2OoV1pCJEWFSRSJZIqmp3mtW7ArzI0N9FySAZVQ3i/7mD+c/4gOcydhQaGs\ny17JyqylRIeohNInMEYNeZEaIKkoYGzp4eDZOg6W1nO+vBlnLiF/YiyLZqSxeEY6mSmB11tJs9no\nqqgkOj/Pr3GYrL3suLyXdy9up623g4jgcDaI5dwnVmEIG5tNXj6bhBJACDETuB/HoMifAd3ANCnl\nVneD8CY1CaV7+wTKIC01QNK9/ZKSDJwta+C4NHK8zMhV5xgTHY5kUiSSWLMoF71t8CWUfLnQV18N\nJbygAH1I6DB+O9eM5Pr12szsrT3IjspP6LR0ER4Uxr2FK1mYUExUSKRX4giU957PJqEUQvwbjqla\nwLGw15tALPCWEOJD4DNSSpO7wbgYy0rgUSnlE744n6IEOk3TqDF2cbyskTPl1/sHLOp1OqbmxFNU\nmMycSYnERYcBkDQh0q9J+9aifNJjjxO/crXf4hmo12ZmT80BPq761JlQwtmQs4oVmUvJSU/x+5ed\nQDfsxCKEeBpHUnkex4DJQ86XdgE/Br6FY2T+dz0b4qCx5AOzAdfWKFWUMcauaVytb+eEc4xJo3PA\nYkiwntkFiRQVJjGrIJHoiMBZbXDQqVdmziIiv8DPkTmavPbWHrwloaxmZWYJkS7eoYxnrtyxfBX4\nPynlPwghEvuelFK2AX8jhEgGHsEHiUVKeQX4iRDiFW+fS1ECjc1uR1a3caLMyIlLAwYshgaxYEoy\nc8XNAxYDjenyJWpe+hHgSCgJm7YQnpPr35isveypPcDOqj10WrqICA7n3pzVrFAJxS2uJJY8HDWV\noRzEkVjcIoQoBl6QUq5wLnX8S2Am0At8WUp5RQjxA6AAeFJK2eruuRRltLFY7Ry70MCuI5WcvNRE\nZ8+NAYtLZqRSVJjMtJx4QoJvH7AYaMILJhG/dj2GBcUBkFBMfFpzgJ3Ve+iydDsSSu4aVkwsITLE\n84uOjReuJJZmHAX7oUwFrrsThBDiOeBxoK9/4xYgVEq52JlwXsIxueU/unN8RRmNes02x+h3aeT0\n5SZMZkehPTY6lBVzHQMWRWYcwUHeW9NkJDRNQ7Na0QXf/DGj0+lIetjt76Ae0eNMKLuq9tBl7SYi\nOIL7ctewXCUUj3AlsbwNfE0I8RpwbeALQoh7cTSV/c7NOC4DDwJ9TVslwDYAKeVhIcS8wXaSUn7O\nzfMpSkDqMlk4fbmJ42VGzl69jsU5YDExNpz1i3KYmhlHXkYM+gBZX2QwfTWUs1vfIzhfkLj5AX+H\n1M+RUPazq2pvf0K5P3ctyzOXEBGsEoqnDPv/TiFEPLAHmASU4lh3ZS8QA8zCsUbLIillozuBCCFy\ngNeklIuEEC8Db0kptzlfqwRynTMtu0SNY1ECXUuHicNnr3GwtJ7Tl24MWMxMMbB4RhqLZ6aTmx4T\nMItVDWXgAlvt5y8AkHTPMsQ3v+HnyKDb0sNWuZv35U66zN1EhUZyv1jFhkkriAxVCWUwXlno61ZS\nyhYhxEIcvb8eAkw4Fve6CrwIPC+lbHE3kFu0c/N6onp3kkqfQOgaGAhjIdQ4Fvd449pdvGzkhDRy\nvKyRSzVt9H37yU41MK8wibkiibSEGyO5m5o6A3qhrwkxoZz+h+/e1Msr/88ewxTn3665PdYejlw/\nxvsXP6bb2kNUcCQb89Zxz8QlRASH09VmpQvvvY/c2S9Q3nsj4eq0+V04en1912MRDG4/jsXF3nQm\nszNePp+ieF19cxcnpJHT5de5XO3oe6IDJk2MZW5hMnNFIomxo/Pbc1BYGPqIiJt6eRmSDJj8lFR6\nrD3srt51h9VOAAAgAElEQVTHrup99DgTyqa89dwzcTHhwWqUgre5dKsjhAgCcoA0HMsS30ZKuced\nQJxNYa86C/Y6bvQKA/iilFK6c1zVFKb4i6ZpXK1r50CpYyqVqmuOD9kgvY6ZBYksmpnOwumpxBvG\nxged3WJBH+Lf8TJd5m4+lLv4UO6iy9KDITSKjZPXsK7gHiJCxsbf2Vd8MqWLcyqX/8PR3XcompQy\nyN1gvEFN6eLePoFyOz7amsLsmkZ5Xd+AxUaMrY4uvyHBeqbnTqCoMIlVxTn0dLm+2Kq/m8L6ivIh\nbUZCFt5z1+19ee26LT3srt7L7pp99FhNRIVEsjrrHpZlLCIzLcnv7z1X9guU956vpnT5JZAK/BNQ\nAQw+yZCijDM2ux1Z1cox6aibtHWaAQgPDaJ4agpFIonpeRMID3W83aIjQ91KLP5y60h5XXAwedPm\nEmTw/6SV3ZZudlXvY3f1Pkw2E9EhUWzJv5elGYsIDw7zd3jjlit3LJ3AD6SU/+zFeDxONYUp3mC2\n2Dh1ycjBM/UcPlffv8KiITKUhdNTWTQjjVmTkvpXWBytWs+UUv3a7/t7ecXPn0fWIw8TXZDv17g6\ne7v4QO7iw0u76LGYiAmLZtPkNazNX0a4avLyCJ/0CgNqGaV3Kf5uSgHVFObu84HUFGYyWyktv+6Y\n5PFKc/+AxbjoUFbOzaCoMBmRGUuQc8ncttbuOx7P3Tg8te1wmsKubf2Y9vMXbloPJdrDsxu7osvS\nza7qvXxSvQ+TrRdDSDQPFNzH0oxFhAWF0tFqoQOLV+MYL01hI+FKYvkh8EMhxAdSygsei0BRAliX\nyULpsSp2H63m7NXrWG2OXu9JceH9y/Xmpgf2gMWRSNi8hbiVq/2+wFanpYvdVXv5pGZ/f0LZkLu6\nP6EogcWVxPIK8ChwRgghgUbgtmYmKeVKD8WmKH7R1tnLiUtNnChr5GJVa/+AxYykKIqEY4xJZnJ0\nwA9YHC5N0zDX1hA28fYZm0ISEglJSBxkL9/oNHexs3oPn9bsp9dmxhAazb25a1iasZBQlVACliuJ\n5QVgLWDFsQbLYOtwqnqGMio1tfY4BixKI5cHDFjMTTOwdM5EJk+MJXXC2JrlVtM0ui+c7y/KZ3/n\nB4Rl3mk6QN8ZLKHcn7uWEpVQRgVXiveNwHHgs1LKdu+F5FmqeK8Mpbqhg4Ol9RworeNKTRsAOh1M\nzU1g8cw0Fk5PIzl+bCUTcE69Unr2tqJ8zhc+R+TEiX6Nrd3UwXtlH7Pt8qf0WnuJC49h8+S1rM5f\nSliwSii+5KvifTiO+btGTVLp4+/iLwRGAXG8F+81TaOqoZPjspHjZUbqmx3F9SC9jul5EygSScye\nlERslPMDzGrDaOwIiGvn6n532rZl18cYX/0fgJuK8l1Al3MfXy5NDNBh7mRn1R4+rT2A2WYmNtTA\nxknrWJJeTGhQCO0tvThW0HDdaLt+gfLeGwlXEst2YCXw7x47u6J4mV3TuFLbxvEyxxiTpjbHgMXQ\nYD1zRRJFIolZBQlEhgfOCoveZpi3gJ6yi+Q/9jA9scl+jaXD3MnHVZ+yp+YAZruF2FADm/M2sDh9\nAaFB4+eajDWuJJbvAx8IId4A/gA04Ki33MTdKV0UxVOsNjtl1a0cLzNyUhpp63IMWIwIC2Lh1BSK\nCpOYnptAWOjoHmNyN32twLe2aATHxJD+5NeJTjLQ46e7+Q5zJzuqPmFvzUFnQolhc869LElbQIhK\nKKOeK4nllPO/n3H+G4wGjO13qxKQLFYbZ69e50SZkVOXm+gyOb7zREeEsGxWGnNFMlOy4wkJDsxF\nsTxJ0zRaT5+h5pXXSNj8AJGTp/g7pH7t5g4+rvyUPbUHsdgtxIXFsiV7BYvT5quEMoa4kli+5LUo\nFMUNPb1WxwqLZUbOXGmm1+IYsBhvCGPhtFSKRBKTBgxYHOtunXoFoPvieZ8mlqd2PTfkayszl7K3\n9lB/QlmXvYJFKqGMSa6sx/I7L8ahKDd58fWTXKhoAR1MyY7n2UfmANDZY+HUpSaOlzVyrqKlf8Bi\nclwERYVJzC1MIjdt7A5YHIrFaOTaf77cn1Di58/DsO5+vw9sHGhX9V5nQlnJovT5hOhdWrVDGUWG\nvLJCiGXAxb4VIZ2P70rVWJSRevH1k5yvcK4Zp8H5iha+/tM9JMdHUNXQid1ZO5iYFEVRYTJFIomM\npKgxM2DRHUExMZivXevv5ZU5f0ZA9IYc6JHCB1iYphLKeHCnK/wJ8Djw6oDHd6NqLMqIXai4fSHS\n7l4rFdc6yEuP6R/9njLGBiyOhD4sjJwf/JCg6Gh/hzKkpRmL/B2C4iNDfsUTQnwB+FRKeXXA47sK\ntCYzNUAy8Fltdq7UtFJ6pZnSK02cuNg46HbxhjD++7vrfRxd4OhbUx4gbtbMu2ztWyZrL1vlbl4r\nfWfIbd747K98GJEyUl4ZIHlrgrhbwnCuLpntbiDeFAhNAoEwSCtQBkjGT4jiWGkdF6taKKtq5VJN\nW3/hHSA0RI/ZYr95H0MYTz800y/X0t/Xrq8o3771PdrPXyA0PZ3s7/3THZv+fLXmvdVuZV/dYbZX\n7abNdOex0/56H/r7+rm637gaICmEsAOPSylfHWKTzwM/xTGPmNcIIVYBnwUigR9JKc9483zKyFlt\ndirqOyirbuFiVStXatv6p5wHSEuIZHJWPPOnpZEWF0ZsdBh//Yv9tHQ4RlrHG8J46akl/grfbwbr\n5dVXQ/F3Pcmu2Tly7QQfXt1Bs6mFsOAwNuSsYmvFTr/GpQSGOxXv04E1OOomff8XLxVCDLaPHvhT\nz4c3qAgp5VeEELNxTIqpEkuAsdrsXK1vZ9fpOk5eaOBSbdtNdyCZKQYKMmIozIyjMCu+fwqVgd+a\nnn5oJj9/6wx6vY6vPzDDL7+H39ntNLzyX1gaG4iaNZuCP3vM7yPlNU3jcM1J/vfkH7nW3UiwLogV\nmSU8Nncj5g4d9+et6982ENbSUfzjTncsTcDfA5MGPPcXzn9D+VdPBHUnUsr3hRBRwNPA0J3mFZ+x\nWB2JpKyqhbLqVi7XtGG23kgkGYlRFGbFMTkrHpEZR35Owl0/cLJTDbz01JJx/eGkCwoi5XOfRx8R\n2b/Alr9GymuaxsWWS7x7ZRtVHTXodXoWpy1gQ+4qJoTHExtuwNgxPq+Tcrs71VjMQoi1QK7zqV04\nFvv6eJDNbYBRSnnRnSCEEMXAC1LKFUIIPfBLYCaOWee+LKW8IoT4AVAAfAPHFP7fllI2uXM+ZWT6\nEsnOU3WcuNDAldpbEklSFJMz45k/I43U2DBiItWstHeiaRq2tjaC4+Juey1yylQ/RHSz8rZK3r2y\nlUut5QAszixidfoKUqL8e/ekBK471liklJVAJYAQ4ksM6CXmKUKI53B0a+50PrUFCJVSLnYmnJeA\nLVLKf3Ru/19AIvC8EOKPUsq3PBmPcjuL1UZ5XTtlVa1crGrhSl07lgGJZGJSFIVZ8UzOikNkxmGI\nvL1pS7ndwBqK2dhI7vM/Qh8SOEm4trOe98q3UdrkmFp/WsJkNuatY27eZHVdlTsKhJH3l4EHcaxQ\nCVACbHOe87AQYt4tcXzeS3EoTharjdLLTRwuraOsqoXLte39I9wBMpOjKcyMY8GMNFJiwvoTiTI8\nty6wBY6ivL2rG32c//+Wjd1NvHrwDQ5UHUdDIz82l0356ymIy737zoqCCwt9eZMQIgd4TUq5SAjx\nMo51X7Y5X6sEcqWU9jsdYyhqHMvd9VpslFVe56xzHElZZUv/HYlOB7lpsUzPT2B6fiLT8hKIifL/\nh99oduXXL3Nt6zbAMfVK1iMPE12Q7+eo4Hp3K/937gN2XT2AXbOTG5fJozM3Myt1qt97oSm+56uF\nvnylHTAMeKx3N6n0CYTb9kDoS9+3j9li40ptGxerWimrbqW8rg2rzTnFOpCZEs2cwhSyEiOZlBlH\ndMSNSQJ7u3sxdvcOK4ZAW+jLXZ6OYcL8IjrqrvUvsNUDwyrKe2scUqe5i+2Vu/m09gBWu5WUyCQe\nm72ZvLAC9Do9TU2dg+7nznP+EAjvPVf2G1fjWHxoP7AReFMIsRDVnXjEep2JpPpYDScvNlBe335T\nIslKMVCYFUehs0YSFR4SMB8KY1F80VysWZPuvqGXdVt6+ODqDnZV7cFk6yU+LI57c9dQnDqX1JQ4\ndf0VtwXE/a2zKexVZ8Fex41eYQBflFJKd489HpvCTL1WLlZed0yRcrmJS9Ut/YlEr4O8jFim5ycy\noyCRqbkJN92RKCPXN/VK3TvvIr75DMEBNn+X2Wpm+5U9/OH8NjrMXcSERfPg1A2syV+qprBX+o2k\nKSwgEos3aZqmBcI3L2/ejveabVyubeufIuVqfTs2e9/qgZCdYmByVjwLZqSRbAgjMvzuN6qBcjse\nCHdOw41hsJHyqX/+BDGLbp41wF9r3tvsNg7WH2VrxU5ae9uIDIlgVeYylk8sITw4bNjnUk1h46Mp\nLDk5xvs1FueULgO//feddOBzZhxLFh8BvielPOduYMrQTGYrJy42cuSsY76tivqO/kSi1+nITjU4\nByTGUZAR159IAuWNPhaZKiswvv7qjV5eM2eRsGkL4Tn+70ll1+ycaDjN+1e3Y+xpJkQfwpqs5Tw6\n93562kdUvlSUQQ07Iwkhvo1jtHscsAO4CJhwjMzf4DzWW87XVzh3WySlLPVkwK4aC01h3SYLFypu\n9Nq6XN16I5HodUyaGNffa2tq7gQiw1Vzhq91XLrMmWf/hvj5RWR+9mEMkwr8HRKapnGi/iyvn3mH\nyrZagnR6VuWX8NDUe4mP8OqUfsoY4KteYTocc4LNk1KeGviCECIXOAScl1I+L4RIAfYB3wcecDc4\nTwmEb+mu3C309Fq5VNPWP0VKRX1H/+JWep2O3DQDcyankJkYSUFGLBFhNy5jV4eJrg7TiGMIlNvx\nQLjLGlYMcSnk/L8XCE1NxQSY3Pi7eCQOp0atnv8+/jZX2yvRoaM4tYh7c1eTGJGAtROMnR0emd3Y\n1ef8QTWFuf+8u1xJLE8AP7s1qQBIKa8KIX4OfA14XkrZIIT4LWour2FxJJJWR/ffqlYqr91IJEF6\nHbnpjhpJYWYcBRNjCQ8NDpg37XjSV0MJSUkhZELCba+Hpqb6IaqbVbZX8175R1y47miSm5U0nftz\n15Ie7f/YlPHDlcQSw41pVwZjwjHVSp9WIMKdoMa6bpMjkfRNkVLZ0EFfg12QXkdeekx/99+CDEci\nUfzn1pHysfesIOVzgTUBRH1XA++Xf8Qp41kAZqRMZkPmGrJjMv0cmTIeufKJdQx4UgjxOyll88AX\nhBBxwFeBkwOeXo5jupZxr9tk4cj5axwpreNiVStVtySS/IxYJmc5ppAvSI8lLFSt7hwI+u5QzjoX\n2ALH1CuxS5f5ObIbmnuu88HVHRy5dgINjZyYLDblraekcI66o1X8xpXifTGwG+gG/ge4hGP24ULg\nURx3K+ullJ8IIbYC64CnpJR+XY/UH8X7zm4z58qbOVvuKLaX17b1J5LgIB0iK54Z+YnMyE+kMCde\n3ZEEqF5jE8e+8iTY7QE19QpAa08bb5/fxo7yvdjsNjJj03l0xiaK0meq6VcUj/DZOBYhRBHwI+Ae\nHIX8PgeAZ6WUh5yF+8PA73B0OfZrryxfjGPp7LFwqdoxPcrFqhaqGzr7+2AHB+nIS49lzuRkshIi\nycuIJSxk5Hcko3lp4tFUvG/Z9THpRTM9usDWSIq/lXUN7Kj6lE+q92G2W0iMSOD+3LUUpcxCr9Pf\ntK0vliZ29Tl/UMV79573yTgWACnlcWCVECIeyANCgHIpZeOAbRqAHHcDGg06eyxIZxKRVa1UN96c\nSERmnLNGEk9+egyhIUEB8yZTbqdpGprZjD4s7LbX4leu9usCW316bWbePr+Vdy5sp8dqIjY0hgdz\nN7I4bT5BetV0qgQWt9pgpJQtwHEPxxKwOnsslFW1UlblWLO91jgwkej759ianBVPnjORKIFv4Ej5\noJgY0p/8ur9Duo3FbmV/7WG2Ve6kw9xJVHAkDxTcx7KMxYSq6VeUAOXKyHsdjgL9nwIpwK2fnjpA\nk1LmeS48/+joNjvvSBzJpMbY1f9aSLC+f5ndwqw48tJjCAlWiWQ0GWzqlajZc9CsVnTBgVHvstlt\nHGk4yYdXd3Dd1EJYUCifmXYvCxOKiQhWnS2VwObKu+jbwHeAFkDiKNzfalSOcm/vNiOdY0guVrdQ\ne0simZId359MctNiCAnW3+FoSiDTNI3an/2E7rOOCSGiZs3un74+EGiaxkljKe+Xb6ehu5FgfTAr\nM5eyNnsFeRlpqjlVGRVcSSx/DnyKo+fXYEklIG169h2mZMfz7CNz+p9r7zJTVu24GymraqW26UYi\nCXUmkr7uvyqRjC06nY6wiZnogoICLqFcuC55r3wbVR216HV6lqQvYEPOauLD4/wdnqK4xJXEkgR8\nfzQlFQBNg/MVLfzlv+xhSnY8dc3d1A1MJCF6pubE96/ZnpsWQ3CQSiRjWeKDn0GnD5xrXN5WwbtX\ntnGptRyAouRZ3Je3lpTIJD9HpijucSWxnAOEtwLxti6TlWNlRkJD9EzLnUChs9iek2ZQiWSM6auh\ndJ46SdIjj902riNQkkpNRx3vlW/jbPNFAKYnTOb+vPVkGtL9HJmijIwrAyTXA78HHpNSfuC9kDzr\n/m/+sb/uExsVyu++s04lkjGqb4Gt6tff6B8pP+unLxKd5/+p6weq72jk92ff40DVMQCmJE3i0Rmb\nmZwUGIMvFQV8N7vxN4AO4D0hRBfQDAxczCGge4XFG8J4+qGZtFzvuvvGXhAIg7TG8gDJbllG8x/e\nutHLy1mU7zEkjngMiqeuXYupla0VH3Ow/hh2zU6mIYNNeeuZMkGgQzfiv7er26oBkr493mgbIDkS\nriSWCBzTuNxp/i+v9wpzjv7/Oo5E9tzAwZlDiTeE8dJTS+62mTKKma6W03NJBlwvL4AOcyfbK3ez\np/YgVruVDEMq67NXMydphpp+RRmThp1YpJTLvRiHK8KAZ4C1wCLgnTttnBAbztcfmOGLuBQ/ilu+\nksjCyQGxYmOfHquJnVV72FW9h16bmfiwOO7LW8t905dxvbnb3+EpitcExmgwF0gpDwghFgHPAg/f\nbfvffXtdQNyOKyOnaRo9lyQRBZNuK8Drw8ICJqmYbRb21B5ge8VuuqzdGEKi2ZS3gSUZxYTog9UU\nLMqYN2RiEUJcBb4hpXx3wOM7NXW5XWNxzpz8gpRyhRBCD/wSmIljEOaXpZRXhBDfx7EM8k9wTOG/\nAceAzW+4ej5ldLl1pHzaV57EsKDY32Hdxma3caD+KFuvfkybuZ2I4HA25q1n+cQlhAffPg+ZooxV\nd7pjqcQxRf7Ax3fjco1FCPEc8Dg3FhHbAoRKKRc7E85LwBYp5bed268A/hMwA79x9XzK6KFpGq2n\nz1Dzyms3FeVD09L8HNnN7Jqd4w2nef/qdpp6mgnVh7A2ewVrsu4hMiTS3+Epis8NmVhural4scZy\nGXgQeMX5uATY5jznYSHEvFvi2I1jXRhljOs6c5pL//ovQOBNvQKOxFfadJ73yj+irusaQbog7pm4\nmHXZq4gNM/g7PEXxm4DokiKEyAFek1IuEkK8DLwlpdzmfK0SyJVS2u90jKH4Y6EvxTM0m40rv3mZ\n1LVrAmaBrT5nG8p4rfQdLjVfRafTsSy7mD+Zdh/J0Yl331lRRgFfjWNBCLGBoWc3BkBKudLdYJza\ngYFf9/TuJpU+gVC8D4S+9KNxHEvB176K0djh1/VQBsZV2V7Nu1e2cbHlEgCzk6Zzf9460qJSoAeM\nPZ4ZpzDS/dQ4Fu/Focax3J0r0+Z/Dfg3HHWUBhw1jlt54u5gP7AReFMIsRA444FjKsqI1Hc18F75\nR5w2ngVgygTBxrx1ZMdk+jkyRQk8rkzpIoEuHLMbN3gyCGdT2KvOgr2OG73CAL4opZTuHls1hSkj\n0djZxBvn3mdvxRE0NCYl5PLYzC1MSx610+YpyrD4ZM17IUQP8IyUclT1xPLFmvfDEQi346OxKcxf\nzSltve1sq9jF/rrD2DQb6VGpbMpfz/SEKSMeLa+awnwrEN57ruwXKO89X615Xw6kunsiRRkNuizd\nfFz1Kbur92GxW0iMSOCxWZuZFCHQ69TkpYoyHK7csTwO/Ay4R0p51nsheZZqClOGw2Qx8eGl3bx7\ncQfdlh4mRMTxmWn3sjx3McFqpLwyDvmqV9gSHLMbnxJClAFGbp7dGPBIrzCPU7fjru8TKLfj3m5O\nsdit7Ks9xEcVu+iwdBIVEskDBfexLGMxoUEhtDR3B8S1c3U/1RTmvTjGS1PYSLiSWDbg6PVVA0Q5\n/91K3R0oo4LNbuPItRN8cHUHLb2thAWFcm/OalZmLSMiONzf4SnKqObK7MY5XoxDUXzCrtk5ZTzL\n++Uf0dBtJFgfzKrMZazNXkF06GDflRRFcdWom91YUdyhaRrnr0veK99GdUctep2eJenFbMhZRXx4\nnL/DU5QxxaXijI9G3nuUKt4rF41XeK30HS4YHaPll2TN4+HpG0kzJPs5MkUJXD4p3vtw5L3HqQKi\n6/sESgFxJH+36o463ivfxrnmiwBMT5jCxrx1TDSkgwmMJs/9LVyhive+NdquX6C890bClaawZ4DT\neGHkvaJ4UkO3kQ/Kt3O88TQAk+Ly2JS/nrzYHP8GpijjhCuJJRPHyHuVVJSA1GJq5cOrH3Po2jHs\nmp0sQwab8jYwecIktba8oviQGnmvjHod5k62V+5mT+1BrHYrKZHJbMxbx+yk6SqhKIofuJJYngd+\nJoR4azSNvFfGrh5rDzur9rCrei+9NjMTwuO5L3cNC1LnqulXFMWPxsXIe2VsMdssfFqznx2Vn9Bl\n7cYQGs2m/A0sSS8mRK960CuKv6mR98qoYbPbOFB/hK1Xd9JmbiciOIJNeetZnllCWFCov8NTFMVJ\njbxXAs5Tu5674+uh+hDWZa9kddYyIkMifRSVoijDNeYrm2qA5Ojz8O+fHPK19ZOW8+CU9cRFxPow\nIkUZf3w1QHI3d27q0gFaINZY1CAt1/cJ1EFaGzPvxdIJxk7fXdNAuHau7qcGSHovDjVA8u5cqbHk\n4kgsA7NYEJAIhAGVgE96iwkhUoD3pZTzfXE+xTeqWmt5T+7ydxiKoozQiGssQohgYBPw78CPPRPW\nXX0LqPDRuRQvMtssnGw8w766Q5S3Vfo7HEVRPGDEfTOllFbgbSFEMfDPwKIRR3UHQogngf8B/tqb\n51G861pXI/vqDnG4/jjd1h506JidOpUFSfP4bel/+zs8RVFGwJOd/i8DT7uzozMpvSClXCGE0AO/\nBGYCvcCXpZRXhBDfByYByc7XFgghHpJSvuWZ8BVvs9itnGos5XDpsf6Zhg2h0azNXsGS9GKmZGUH\nRJu8oigj45HEIoQIxzGdfqMb+z4HPA50Op/aAoRKKRc7E85LwBYp5bdv2e+/VVIZHRq7jeyrO8zh\n+uN0WroAKIwvoCRjITMTpxJ8y6DGX6z8Uf/PgVIAVhRl+DzRKywMmAzEA99xI4bLwIPAK87HJcA2\nACnlYSHEvMF2klL+mRvnUnzEardysPo4H174hLKWywBEh0SxOuseNk5fSbApws8RKoriLcPupyyE\nqBjiJRtwDXgV+KWU0uVxI0KIHOA1KeUiIcTLwFtSym3O1yqBXCnlbdPHDIcax+JbDZ1GdpbvZ3f5\nAdp6HXcaU5MmsTp/KcUTZxMSFOLnCBVFGQ6fjGMBiqSUze6eyAXtgGHAY727SaVPIDSlBEJfem+N\ng7DZbZQ2X2Bf7SEuXr+EhkZkcAT3iVXMjZ9DapRjpcbW6ybA5LOFvjwlEK6dq/upcSzei0ONY7k7\nVxLLKSHEb6WUP/DY2Qe3H9gIvCmEWAic8fL5FDc1dV3n/fJdHKg7QpvZ8T9lXmwOJenFzEmeSUbq\nhID4YFEUxbdcaQrrAZ6WUr7s6SCcTWGvOgv2Om70CgP4opRSunts1RTmWXa7nRP1Z/n4yl5OXjuH\npmlEhkSwLLuY1fklZMVl+DtERVE8YCRNYa4klv8ApgObpZTX3D2hr2mapgXCt+ZAuB0fSVNKa28b\nB+qOcKDuKC29rQBMmpBDcfJ8ilJmETrI7ML+XvPeUwLh2rm6n2oK814c46UpLDk5xic1FhswFagW\nQlzG0bXYdutGgThXmOIeu2bnwvVL7K89RGnzBeyanbCgUEoyFlKSvpC5eYUB8cGhKEpgcbVX2K1z\nhd1Kk1LmjjQoT1JNYa5r7Wlj19UD7Czfj7HL0V8jNz6TNfnLKMmaR3hIuJ8jVBTF23zSFDZaqaaw\n4e1j1+zIlivsqz3E6aZz2DU7ofoQ5qXMoSSjmOyYTJdjUE1hnj2eagpzz2i7fuOtKUwZgzrMnRyq\nP8b+usMYexx3JxnRaawX9zAlegoRwWogo6IorlGJZRzSNI1LreWOuxPjWayajRB9MMWpRSzNWEhO\nTBbJyTEB8W1TUZTRRyWWcaSzt4tdVXvYV3eYhm4jAKlRKZSkF1OcOlct86soikeMixqLv2PwJ03T\nKGsqZ8eVPRyqPoHFbnXcnWTOZU1+CZMTCxhBjU5RlDHKV1O6jFqB0KTj6wJit6WHI9dOsK/uEPVd\nDQCkGZJZlLKA4rQiokOiAGhq6hx0/0ApIAZCAXi0FX+Hu60q3vv2eKOteD8S4yKxjBeaplHRXs2+\n2kMcbzyNxW4hSBdEUfIsSjKKWTxp9pCJRFEUxVNUYhkDeqwmjl47yb66Q9R21gOQGD6BJRnFLEqb\njyE0GkA1eSmK4hMqsYxiVe017Ks7xNGGU5htZvQ6PbOTZlCSUUxhfAF6nd7fISqKMg6pxDLKmKy9\nHG88xaGTRylvqQJgQng8S7JXsChtPrFhMX6OUFGU8U4lllGipqOOfXWHOXrtBCZbLzqdjhmJUylJ\nLwnx7dcAAA+ySURBVGZqQqG6O1EUJWCoxBLAzDYzxxvPsL/2EFfbHXcncWGxrMxaxsbpK7B3qcun\nKErgGfPV3NE4jqW6rY4dV/ayt+IwXZYedOiYnTaNNfklzEmbTpA+yN8hKooyxqlJKO9gtExCabFZ\nOGksZV/tIa60VQAQE2pgcdp8FqcvICFigkvHcycGV7dV41h8ezw1jsU9o+36Bcp7T01COYo1dDWy\nr+4wh+uP02XtBmDKBEFJejEzEqequxNFUUadUZdYhBCzgH8FrgD/JaX8xL8Ruc5it3LaeJZ9tYe4\n1FoOQHRIFGuylrMkvZikyAQ/R6goiuK+UZdYgAVAPWAFzvk5FpcYu5vZX3eYg/VH6bR0ASDi8inJ\nKGZW0nSC9aPxciiKotxsNH6S7QNeB1KBZ4G/8W84d2az2zjTdJ4j545xpuECAFEhkazKXMaSjGJS\nIpP8HKGiKIpnBURiEUIUAy9IKVcIIfTAL4GZQC/wZSnlFSHE94FJwLs47lhaCZD4B9Pcc539dUc4\nWH+UdrOjKJYfm0tJRjFzkmYQEhTi5wgVRVG8w+8fzEKI54DHgb7ZEbcAoVLKxc6E8xKwRUr5bef2\ni3DUWCzA9/wQ8pBsdhtnmy+yr+4QF5olGhoRwREsn7iETdNXEWaO9neIiqIoXuf3xAJcBh4EXnE+\nLgG2AUgpDwsh5g3cWEp5EDjo0wjvosXUyoG6IxyoP0prbxsAuTHZlGQUMzd5JqFBoSTFBkbXS0VR\nFG8LiHEsQogc4DUp5SIhxMvAW1LKbc7XKoFcKaXdnWN7a4Ck3W7n1LVz7LiylxP1Z9E0jYjgcJbm\nLGBN/lKy4yZ647SKoig+MdYW+moHDAMe691NKn08eafQ2tvGwbpj7K87TEtvKwBZhokszVhI0f9v\n796DrqrKOI5/AYGEaFRkQDBN0Sd0Kt9EfEGgi5dMTWPSvOTdkSayvKUmZpZmFzUZL6Rjltdg1Lyn\nxmiJCSivYmaGlydJLTWE10uoeOG99MdaG7abc973nMN52eec9/eZcV7OPmut/Zy9nL32WmuvvYc3\nMbDfAFi19j5rYZGWFkhWphbqrtx8WiDZc3H0lgWS66IWG5YFwL7A781sPPD3nOOho7ODZ19/jnmv\nLOTJ1qfo6OxgQL8BTBzZzKRRzWwxRL0TEZFELQ2FzY4T9n1Yc1cYwNHu7pWWvS5DYW++t4IHnn+Y\nPy+Zz6vvtALwiY02Z/fRk5m05TgG9d+w0qJFRGqanhXWhXKfFdbZ2Ym/sYT5ryzkieWLae9sp3/f\n/owdvgOTR41nyyEfr+hNjLXQHddQWGVqoe7KzaehsJ6Lo7cMhelZYVXw9gfvsHDpIha83MKyd0Pv\nZOTgEUwc1czOw3dU70REpES9osfSxXc8vfw5/rRkHgtfepy2jjb6992ACR8fyx7bTMaGbq33xItI\nr9Rod4VVXbaLt3LVSlqW/pX5Ly9k6cplAAwfNIxJo8bTPGIsg/sPgk5obX27UHEVqYXuuIbCKlML\ndVduPg2F9VwcvWUobF00fMNy4I3TPvS5ecRY/rrsCVZ1tLFBn37sNLyJSSOb2WYj9U5ERKqh4RuW\nrJaljzFsw6FMHNnM+M12YsgAPWZFRKSael3D8t2mqdjGo+nbp2/eoYiINKSGH/v5+g3f+tDk/U0H\nXZ5XKCIidUOT92XIazKxFiYQNXlfmVqou3LzafK+5+LQ5H33NB4kIiJVpYZFRESqquGHwm466PKa\n6I6LiPQW6rGIiEhVqWEREZGqUsMiIiJV1fDrWHrq1cQiIo1M61i6UQuT97VwL73WsVSmFuqu3Hxa\nx9JzcWgdS/fqrmExs+2BE4ABwC/dfXHOIYmISEo9zrEcC7wEvAe8kG8oIiKSVY8Ny2jgUuBm4Iic\nYxERkYyaGAozs2bgF+7+RTPrC1wGfAZ4HzjW3ZeY2TnAtsByYCXwBvXZMIqINLTcGxYzOw04DEhe\n1zgFGODuu8QG50JgirufFdOPBa4k3NF2Qg4hi4hIF3JvWIDngK8B18fPk4A5AO7eYmY7pRO7+2PA\nkes1QhERKVnuQ0nufivQlto0BFiR+tweh8dERKQO1EKPJWsFoXFJ9HX3jkoLW5dFPiIiUr5a7Aks\nAPYGMLPxwN/zDUdERMpRSz2W5NErtwF7mNmC+PnonOIRERERERERERERERERERERERGpP71yjYeZ\n7Qoc4u5T845FymNmuwEHAYOA891dt6PXifg4pu8QzjunufuynEOSMpnZcOAudx/XVbpaXMfSo8xs\nNNAEfCTvWKQiG7r7N4FfAl/KOxgpy0DgROBuYELOsUiZzKwPcColvK6k1zUs7r7E3WfkHYdUxt3v\nMrPBwPHANTmHI2Vw94eA7YFTgL/lHI6U71vA7wjvwupSLS2QXGelPH4/1wClSyW+PmFT4HzgLHdv\nzTFcSSmx7sYBi4C9gB+hp5PXjBLPnbvHbTub2f7ufkux8hqmxxIfv38lobsNqcfvA6cTHr8vNaqM\n+rsQGA783Mz2X++BylrKqLuPAlcBFwCz1necUlip9efu+7v7NKClq0YFGqvHUu7j9w9fv+FJN0qq\nP3fXKxNqT6l1NxeYm0uE0pVyz53dvrm3YXosevx+fVP91S/VXX3rifpr5Mqu6uP3Zb1T/dUv1V19\nW+f6a+SGRY/fr2+qv/qluqtv61x/jTTHktDj9+ub6q9+qe7qm+pPRERERERERERERERERERERERE\nRERERERERERERKThmFmHmV2ddxxdMbOBZnaVma0ws/+Z2T41ENPWmc8PmNnzZZZRdh5pHI34SBeR\ntM7uk+RqKnAUcB3wIPBYnsGY2dHAr4BBqc3nZj6XqtaPvfQQNSwi+fpM/Hucu7+TayTB51nzwicA\n3P1POcUidaqRn24sUg8GANRIo5Lok3cAUt/UY5GqMrMXgD8SHr09Hdga+A9wkbtflkn3vLt/sUD+\n1dvj5z8ATwCnAZsD/wCOi+VeAnyZ8A6Ja4Ez3T09BNPHzM6I6TcCFgLfd/dFmf1+BTgD2IHwnu/7\ngenu/s9Umg7CsFAT8CXCm/d2cPf2Isfiq8D3Y/r3CUNdZ7r7k6ny0mX/JXs8Mt9PA7YCjiE0SPPj\nb1mcStcfOAU4GNiG0Eg4cLG7X50pL/tb3gAmpr6/xt2PMbMHgC3dfatU/jHAOcCuhPPI48AP3X1+\nofhjnu2BnwJfiPE/Dpzj7vem0gwEzgP2A0YCy4A743F7s1jZUlvUY5Fq6wT2Ai4GbgJOBN4BZprZ\nXpl0hcbgs9s7Ce/gPhv4dfw7BrgFuI/w5ruTCY3NdCD7yukDgJOAywgnwu2AB+JJDgAzOwq4A3gL\nOBWYAUwAWsxs20x5JwH9ge8CV3bRqBxHePx4vxjXDKAZeCj1qtfDgXnx34cRTvRdmQ4cEcs6H9gZ\nmG9m26TSXE04RnNjjGcT3jX/28zxX+u3xLTpeK5IpV1dJ/GYtBAaiEsIDfImwH3Z19im8nwaeJhQ\ndz8FfhD3fY+ZHZhKOhM4FphNaEhvBr4J3Fj0qEjNUY9Fqq0PoVfR5O7/ADCz24FXgEMJvZkkXbH8\n2c+bEXoGi2N5mxAagPnu/o24bTbwOrAHYSI8MRAYn8p7M/A0oZE5wMw+RmgEb3D3Q5NMZnYl8BTh\n6vlrqfI+AKa4+/vFDoCZDSWc+FuAye7eFrdfBywmTI43u/ssM9sjppldrLyUYcD27v5iLO92Qk/u\nHOAbZjYCOAT4hbv/IBXPbcAzwJ6sOf4Ff4uZHVZCPOcSGszx7v6vmO8GYAlrektZlwKvAju6+7sx\nz6WEnuHFZnZrPE6HAr9x9zNTMb0N7Glmg9x9ZXcHSfKnHov0hGeTRgXA3V8lnFSGV1jekvRwD5AM\nT92W2sdKYDmhEUqbk87r7kuAewgnqr6EhmgIcIeZbZr8B7QTrvqTdIlHumpUot2ADYELk0Yl7vtF\n4HpgnJlVcixmJY1KLO8pYA6wT/y8NP6W1T0fM+tDnMch9FzSSvktHxKPxd7APUmjEvf9OjAJOL5A\nnqHA5wjHfXDqGG8M3E74/2JcTP4f4GAzO9LMNopln+XuzWpU6od6LNITlhfY9gHhKrcSr2Y+Jyfr\nZZnt7ax9sfRMgfL+RRjDHwaMjttuKLLvzpguiSG7z0KSuYhnC3yXxLMla/+u7iwusO2fwD5mNtTd\nXwNWAYeZ2Z6AEX5f8v7y7LEp5bdkDQUGs6ZxXy3T+Kclx/h4CjQ8hGO8BWGobBphCPVqoM3MHiZc\nQFzl7isqiFdyoIZFekJH90mKKtT4tBXYBqWtkyiUJhlua0vtbypQbEFfetK44JxKkfILSU7uH5RQ\nTtaqAtuS+NvM7COEOZImwhDTvYT3lT8I/LtA3lJ+S7H9lbNGJckzk9BDKeQpAHe/38y2APYFvkK4\nsWAGcJKZjXX31vJDlvVNDYvkpZ3Megkz2wDYlAJXw+vgEwW2GfCmu78W7zoDaHX3+zPxTAb6lTtc\nxJoGajvgycx3n4x/XyqzTAh3eWVtCyx39/+Z2RHAWOAYd78mSWBmIyvYVzGtwLuFYjGzU4AR7n5K\n5qsX4t/2Asd4DKGHtzLe0dYEvOzuNwI3xqG8k4ELCHM3M6v4W6SHaI5F8rIUGBOvshP7kWlsqmBv\nMxuVfDCzTxEmse+Mm+4D3gNOjQ1bkm4zwm3O51Wwz6TMk+PJMilzc8LdVi2ZK+9Sr/4PN7ONU+Xt\nEH/LzXHT0Pj36Uy+E+LfUi4k22PZBXtdcc7oXsJx3TwVy8aEGyq2KpDnv8Ai4Kh4XJM8GwC/Jdzh\n149wZ9lCwt1vSd7OmBeK91ylxqjHInmZTbhTaI6ZzSJcAU8FXqS6C/TeA+aZ2cWEuYYTgdeAMwHc\nvTWuc5kBPBzvLusHfJvQyGWvvrvl7q+nylwQyxwSy4S15xlK/b2DgUfM7PLUb3kF+HH8/l7Cyfd6\nM5sZ/70v8FnCvNfHSthHMu9ytpnNdfe5BWKcTrjj7ZG4n7cIdTeIeFwL5DmeMDz3WIy/ldADmQCc\n7u5vAJjZtcC3zWwwYc5lKPAdwoXITSXELzVAPRaptmJX39ntlwE/IlzhXkK4a2gKYT1Kdh1Lqfsp\ntO0KYBZhrcXphIWbu7j76qEod78IOJBwIj43pnNgV3eft1aJJYhlHhRj+hnhxDqfcJvxo5mYS+2x\nzCQuFiScbO8GJrj78rjPxcD+wNvAz4GzCENuTcCfgYlm1t0NFJcDjxIWo55aKEZ3f4bQIDwS050N\nvAxMcveni+RZSFh8uYg1Q1uDgSPd/fzU/qcBPwF2IdwG/j3CvNGkeOeZiIhUQ3xS88/yjkOkFOqx\niIhIValhERGRqlLDIiIiIiIiIiIiIiIiIiIiIiIiIiIiIkX8H7Mx5mK+qN22AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xadf2498c>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(N, pyt, 'o-', label='python')\n",
    "plt.plot(N, fort, 's-', label='fortran')\n",
    "plt.plot(N[-3:], 1e-7*N[-3:]**2,'--')\n",
    "plt.xlabel('number of particles', fontsize=18)\n",
    "plt.ylabel('running time (s)', fontsize=18)\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The dotted line has a slope of 2, and eventually both python and Fortran gravitate to that. Note that the Python code exploits vectorization, so its performance is not as bad as it would have otherwise been. Indeed, its relative performance in this case seems to improve as the number of particles increase from 500x worse to only 10x worse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEtCAYAAAAvAJYdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VFX6wPHvTHojCSG0IJ0jTaRIRxELiKigYO/u2ndZ\nV7Ht7s9dV1fRVddVV3eXXcG1oagUxYJKaKFDKFI8EHooCRAIkEYy8/vj3gmTyYRkkknuzOT9PA/P\nMLec+87czLxzT7sghBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQIcTarAxCiPiilFNBLa/2Z2zIH\nsF5r3ce6yM5QSm0DZmmtH1dKxQFvAmOBGOBHrfXVfjxWfyBJa/29v8oUwpPd6gCE8Del1PnARmCI\nl9XOBg7HK6VUR6AT8K256A/AXcB24G/ADD8eawywHOjmrzKF8Cbc6gCEqAfJQITVQVRjFFAILDaf\n9zUfb9Za7/DzsVKR2gjRAOQKRYSyQP4SHQUs1FqXmM+jMK6ejtTjMQP5/RAhQP7AhOWUUtOAO4AW\nwGvA1UApxq/332utNyulbgfeA/6itf4/j/1jgUPAeuB74I8eh7hYa73IbENZBzwEPA8MAoqB+cAk\nrfVuj3IV8AxwOZAE7AU+N2PI9xJ/U+BF4FogEdgEvKC1/sKj3AjgsFn2evP4ntprrfcopaKAx4Db\ngI7ASWAJ8JzWeo1bmXcB7wI3APcCw4GDwELgds+yzbLmm+/FcGAckAeM11ovU0q1A54CRgKtMc7H\nz8AUrfW/vBz3MoyrrPuBc4B95vLJWmuHl9cnQpBcoYhA8i1wMfAf4EfgKiBDKdUL44v8JHCzl/3G\nAXEYCWeB+QhGu8GfgF1u23YA0gEHRiP4OmA8sMhMTAAopQYCa4EbgQxz2xzgcWC5UirZSxzfY1x5\nTAc+BHoAM5RSl3tsNxhIAOYBO4FnAVcym2zGfFwpFQ38gJH8TgNvux1jqVLqGi8xvAmkAK8DKzHe\nt9nmum9dZbtt/0egH/B3YA2wVinVHliNkYgyMJL8FxhtMO8opR72ctyXzLIWAm8BsWbcf/ayrQhR\n0oYiAkkScL7W+iiAUuo64DPg71rrEUqpz4E7lVIDtNYr3fa7FeNK41Ot9XGllA24E1iutfb8QksE\nfqe1nuxaoJSaidG7aiQwSykVBryP0Q4zRms9z23bF4Engb8Cv/QouxToobUuNLf9ESOx3IORCFxG\nAfu01lvM588qpUYAbTF+0eeb+/8fMBSYCtzr+qWvlOqDcZUyTSnVTmt9wq3sEmCY1rrILeam5uv7\nVmv9hrnMtTrefM9z3LZ/CuNq63Kt9Xy35W8BKzCS+j88Xnsns5wd5rZvAtp87X9ANApyhSICyXOu\nZAJgVhVlAMOVUq2B/5mrbnFto5RKxUgEX2mt3X95V6UAeMVj2VfmYwfzcQjQGfjYPZmY/gjsB241\nq67cveVKJqZvzMd2HtuNAr6rQax3AaeAie7VRlrrTIwv9CTgOo99vnFPJjWQ4Z5MTO8D97gnE/O4\nq4AioLmXcj5370xgVh9uAVoopSJ9iEcEMUkoIpCke1nmuhLppbVOB/YANyilXH+7NwKuK4qa2KO1\nLvVY5moIjzcfe5uPizx3NhvRV2E0onf1XO3x3JXgolwLzATYB6O6q0pKqQSMBJeptT7lZZMM87GX\nx/KdZyvXi0rba60ztNbvKaWaKqUuUUrdp5R6VSmVgfFawryU4/nawcvrF6FNEooIJNlelh00HxPN\nxw+AlsAI8/ltGA3cX9fwGGf79e7qpNLEfKzqime/+RjrsbzY/YnW2jXmxb3zy+UYvbl+OHuYtY6h\n0HPDalTaXimVbHY0OIgR5zsYVWabMF6jt848xV6WeXv9IoRJQhGBJMbLsiTz8bD56Gpwv8HsiTQA\n+MTLVUdduNok0qpY72qQr00X31HAKq11noUxVOcDjF5rUzCq/5porTtrre9DkoM4C0koIpAM8LJs\nMEYPpzUAWuttGA3DV5n/oHJ1V11Hw2eajxd6rjCr2oZhfOHv9lx/NmZngZFUU90FYDbM7wTOVUo1\n87LJRebjphocusbvh1IqCRiNkfQe1lovd1W5mb2/opCkIqogCUUEkj8rpVztGCilJmB0I56ttT7m\ntt17QCtgErDNo8cXGAkIoLaNwUswpkC5Tik12mPds0AbjB5lpyvteXa9MMba1KRBHmAaxlXb38ye\nZwAopfoCv8YYN/JlDcpxxVmTtowSjC7Vye6dDpRSMRjdgSHwZyEQFgmqbsNKqRYYvXn6Wx2LqBfd\ngXVKqbkYX9rjMAbIPeqx3XSMcRbtMAYHetpnPt6olCoApmmtN9c0CK21Uyl1J8YX/5dKqS+BHRjV\nPwOBzRjjUXw1CsjHGB/jjecv/5fNfW4Feiml0jES0jiMq47btNYna3Bc1/vxoFLKNUbFK611gVLq\nC2ACsFIp9T1GZ4WrMdqfNHCOUsrm1kZ0NnI104gEzRWKWV3wOBUHqYnQcgPGQMN7ML68pwIDtdb7\n3Dcyr1bmY3ypfuBZiNZ6D8bYByfGSPALqjlupS9GrfUyoD/wiRnLQxjtOc8BAzyumJzeyvBiJMYs\nwt5GjlcqQ2tdjDEC/RmMq60HMK/YgMFa6y/Ptr9bOYswuhk3NV9Ht2ri/QVG0kkCJpoxvIcxEv5b\nINqMw/3Y3tT0fREhImh+PSilHgSWAY9prT2nkhBBzG3qkp41uZIw2zH2AFla6+H1HJ4QooYCosrL\nnOZisjka2o4xxUQvjK6Iv9RaZ2H8SuoFDFBKjddaf25dxMJi92LML/WU1YEIIc6wPKEopZ7AGEvg\nqgseB0RqrYeYieZVYJzWery5/f8kmTROSqlPAIXxw2IL8LG1EQkh3AVCG8p2jOkjXNVvwzBvOqS1\nXoFH/bfW+o4GjU40hJrWtR/CSCgrgLFa67J6jUoIEXyUUu2VUsvM/09RSl3htm632zQbQgghApTl\nVV5e5GNM7e1ir839FBwOh9NmC5o+B0IIERBsdfjiDMSEkoHR532GUmoQsKE2hdhsNnJzT1S/YT1L\nTU3wWxy1LcuX/Wqy7dm28XVdVdv7832rLX/HEAjnr7brfVkeCOfO33EEwrmrbpvarPP3uQqkhOKq\nQ58JXG7ObApwt0XxCCGE8EHI1gk5nU4ZUCWEED4KtSovv5HLbt/3kyqv+oshEM6fVHk1bFmB/Nnz\nNb6akN5TQggh/EISihBCCL+QNhQhhBDlpA2lClKP6/t+0oZSfzEEwvmTNpSGLSuQP3u+xlcTUuUl\nhBDCLyShCCGE8AtpQxFCCFFO2lCqIPW4vu8nbSj1F0MgnD9pQ2nYsgL5s+drfDUhVV5CCCH8QhKK\nEEIIv5CEIoQQwi+kUV4IIUQ5aZSvgjQM+r6fNMrXXwyBcP6kUb5hywrkz56v8dWEVHkJIYTwC0ko\nQggh/EISihBCCL+QhCKEEMIvpJeXEEKIctLLqwrS08T3/aSXV/3FEAjnT3p5NWxZgfzZ8zW+mpAq\nLyGEEH4hCUUIIYRfSEIRQgjhF5JQhBBC+IUkFCGEEH4h3YaFEEKUk27DVZCui77vJ92G6y+GQDh/\n0m24YcsK5M+er/HVhFR5CSGE8AtJKEIIIfxCEooQQgi/kIQihBDCLyShCCGE8AtJKEIIIfxCEooQ\nQgi/kIGNQgghysnAxirI4Crf95OBjfUXQyCcPxnY2LBlBfJnz9f4akKqvIQQQviFJBQhhBB+IQlF\nCCGEX0hCEUII4ReSUIQQQviFJBQhhBB+IQlFCCGEX0hCEUII4ReSUIQQQviFJBQhhBB+IQlFCCGE\nX0hCEUII4Rcy27AQQohyMttwFWTGU9/3k9mG6y+GQDh/Mttww5YVyJ89X+OrCanyEkII4ReSUIQQ\nQviFJBQhhBB+IQlFCCGEX4RuQnnySSJ/+A5KSqyORAghGoXQ7eX18sskAo7EJEpGj6H4mnGUXDQC\nIiOtjkwIIUJS6F6hZGRQcP9DOOPiiJ7+IYm3XE/ihGusjkoIIUJW6F6hDBnCqS7ncerZFwhfvYqo\nL2dS1qmL1VEJIUTICt2E4mK3UzpgIKUDBla5SezrrxCWtZ3isddKtZgQQtRS6CeUGoiY/wORy5cS\n/clH0uYihBC1FLptKD44Putr8r76vlKbi/3AfqtDE0KIoCFXKFChWszV5hKxagWOdu0rb+t0wunT\ncuUihBAe5ArFk5lcCh+e6HV1+JpVpPToTMLEB2WcixBCuJGE4iN7bm6FarGUHp1J+PUDhK9cYXVo\nQgSc8DWrrA5BNKCgSShKqX5KqalKqWlKqeZWxVEyegxH124ib65bm8snHxGxVj44Qriz78gifMtm\nq8MQDShoEgoQBTwCzAUGWxqJ3U5p/4Gcem5yeXIpGn+j9223bJFqMRGS7Dt3wD33VLk+6pu5FI+5\nGoDIr+YQ88+3SLjvLmL/+mJDhSgaWNA0ymutlyqlBgOTgBusjqecmVy8cjjgsstIOVUgXZFFSIn+\n77+IWL8ODmZXuY39yGGcyU2x79yBPf84hQ/8CoqKaDqkH2UdO1E8PnA+xsI/AuIKRSk1UCmVbv7f\nrpT6p1JqqVIqXSnVyVzeH1gNjAYetTDcmisogBtuqNTmEv/Iw1BWZnV0QtRa0S/up+imW6tcH5a1\njbLOxswU4Vu3EPvyC8aK6GhK+/QjYuXyhghTNDDLE4pS6glgCkaVFsA4IFJrPQR4CnjVXB4PvAv8\nFfiwoeOslfh4+NvfjGoxt3EuYTuyICzM6uiEqBuns8pVkd98XV7dVXLZSI5//Hn5Ovv+bErVufUe\nnmh4gVDltR24DnjffD4M+BZAa71CKXWB+f90IN2SCOvKY5yL7cgRr5uFbd1C2L49Ui0mgp497yjO\nxCTjSUQEZd26AxC2cQO2Y3kU3XKHhdGJ+mKzOgAApVR74GOt9WCl1BTgc631t+a63UAHrbXDlzKd\nzrP8fApUDz0E77wDSUkwbhxcfz1cdpkkF1F7q1bBq6/CggWQk2M8XnTRmfXLlsETT0BGBigFO3bA\n73/vvawePWDChDPPFyyAZ5+FdI/feVu3wooVcOedFZcXFsJNN8Hf/w7t29f9tYl6YbPZap0XAuEK\nxVM+kOD23O5rMnHJzT3hn4jqIDU1ocZxhI+7gSiHjag5swibNg2mTcORmET+u+9z+sLhPpVV2xhq\nsu3ZtvF1XVXb1/a1+pO/Y7Dk/LXvCm9OIeGBe4j88Qfsr79Obrc+Z9Z37glffE2TW68n/4NPSW3e\n5Oznw21dxLECkqj8OUudPZvD19+G02N57At/pujZyTjiUrCvWIejY6cavaba8uf5C8XPnq/x1YTl\nbSheZABXAiilBgEbrA2n4ZT26Wd0Rc7cfKbNpUkTSrt2tzo0EcTs+7Mp66wouuV2mDMH+949Fdfv\n3sXpoRdB7X+YVpSXhzOhSYVF0dP+S8nlo3CGR2A/sJ/IRQv8cywRUAKpyusjrfUQpZQNeBvoZa6+\nW2utfS0zKKu8vHE6vX/QS0th4kS46iqpFhNn98EHRhVTmzbQuTP89rfw17+eWf/uu9Cnj/GvpqZM\ngdmzYflyuP9+o3pLKdi8GdasgdtvP7PtkiVw8cVGN3qXzz6D666r6ysT9aAuVV4BkVDqg9PpdFpd\nZQL1d9kdsXghSeONXjTVTbkfyJfdUuXl3/28bRv/9CROPjcZwsNJvf9OHOnpHFm3FWJjz6x/8ZWz\nHqumy2PeeI34Jx4lt6jGL6/eSJVX9eu8LW/evEmt80IgVnmJGjg99EKvU+43ue9uq0MTgeZ0KYSb\nzaW/+Q2248eJ/vTjM+v9OCbKnp8PCQnVbyhCkiSUYGV2RS6f/sVMLkXXTah+X9Fo2Hfvosz9NgzD\nh1PWvScx//0XAGH6Z8q6KP8ca9dOSnv09EtZIjiFdJWX1TEElEcfhbw86Yrc2EydCr16Qb9+FZf9\n4hfw3XegNQwfDj0lEQiDtKF4EeptKD7t1yye0h7nEb5lE2C2uVxxJcVjr6Vk+CUQEeHTMaQNpWHL\nq0s9fNzvHufU8y+B3X5m/b7DpPTpxuk+/XC0bM3JV/9e7bF8WR4I587fcUgbSs1IlVdjYLORl55R\nsc3lk49ocuct2ApOWR2dqEc2h6M8mZSLiqLwjruJ/PF7KD1tTWAiJElCaSw821zmfs/JF185Mz2G\nu9JSmXI/BISvXEFZq9Ze1xXdfS+Eh1Pau28DRyVCWUhXeVkdQ9CaM8cYVyDTvwSnrVvh8ceNKVGi\no42xStOmVd7uzjuNaVaUfxrlRWiQNhQvpA2ldvulpiaQ/85/iHv+T4TtN+514RrnUnjPvZT27itt\nKA1cnj/r4Wu7XtpQrD931W0jbSgiIBVPuLHSlPvR0z8kbPs2q0MTQgSwQJwcUgQCjyn3w1evoqy7\n9znFwjeup/TcblItJkQjJwlFVM9MLl4VFJB09SicEZGUXHEl3HErnD9QkosQjVCN68qUUj2ALbWd\nSr6hSaN8AzlyBP7yF5gxA/btM5YlJcHNN8Pbb1sbmxDCZw3SKK+UOgRM1Vo/VduDNSRplK/dfrVu\nGHQ4CF+ziuTv51L26QxK+15A/rvvV7ufNMr7dz9plK+fOALh3FW3TSA0yvtS5RUH7KztgUSIs9sp\n7T8QrryMo0/9CVv+ca+bha9ZhT3vqNzmWIgQ5Esvr9eBR5VS/esrGBEi7HacScleV8X+4w0Sb7me\nlB6dSZj4IHz9tQyiFCJE+HKF0g9oDaxQShUARwD3ea9tgFNr3dGP8dXaNZNm061dMpNu8uGmQaLe\nFTw8kbK0NKK+nE309A9h+oekJCZx/LPZlJ4v50qIYOZLQokBVnP2dpeAaQh3OmHzrjwe+0cGE8f3\nol1LuUdDICjt15/Sfv3LuyIn/zAXx9yvKVVdrQ5NCFFHITtS/qpHZ5Unt5TEaKY9M8rKcERt5OfD\nI4/AhAky/YsQDaQuvbx8HoeilGoKXA60A0qAvcD3Wuv82gZR30pLHZb1OpGeJmdfd7beJ/mffEGT\nqVNh6tRqb3NcH6SXV+2WSy+v2u0XKL286sKnqVeUUg9hJJCPgcnAa8AM4KBS6mG/ReVn+adKeOuL\njWzadRSHDE8JGsXXXEve3Mq3OY5/8lGrQxNCeFHjKxSl1FjgLWAt8FdgK0ZCOhd4FHhDKbVHa/1l\nfQRaG0nxkYwd1oH0zGzW6lzW6lxaJMcwok8aQ3u1Ii46ovpChHXMrsil/c3pX9asImrOTEouv8Lq\nyIQQXvhS5fUUkAkM1VoXuy3PVEp9ASwFHgcCIqGkJEbzq2vPo13LBC46vzU7DuSTvjablVtymD5/\nO58v2sHAbi0Y0TeNDq2aWB2uqI5bcqlKwoO/hIiIBq0WE0Kc4UtCOR942iOZAKC1LlFKfQA857fI\n6mjaM6PK6wZtNhudWifSqXUiN13ahSUbDrAgM5slGw+wZOMB2rVM4JI+aQzo3oKoiDCLIxe1cvo0\nEatXErZ7F9HTP6zY5nLxpRAu09YJUd98+ZQVA/FnWR9PxXEpASk+JoIrBrZl5IBz2LzzKOmZ2azb\nfpip32zlk/nbGXpeKy7u05pWKXFWhyp8ERHB0RXrCF+9iqgvZxI1ZxbR0z8k8pu5HNm03erohGgU\nfEkoC4GHlFJTtdb73VcopdKAh4DF/gyuPtltNnp2TKFnxxSOHC9i4fr9LFq/n+9X7+X71Xvp1i6Z\nEX3S6N2lGeFhctuYoOBlyv2wfXu8V30VFRn3WpdqMSH8xpfJIc8DlgOlwPvAz+aqbsBtGMlpqNY6\n099B1kZtZhs+Xepg+U8H+HrpTn7KOgJA0ybRjBrUjlGD2pGSGOP3OIVFpkyBJ56Q2xwL4aHBbgFs\nzuP1JjDAY9VqYKLWenltA/G3us42nH34FAsys1n60wEKi8uw22z06dKMEX3T6NYumZq+59IX/uzr\nrJptOPp/U4l97eVKtzkuuP9hynr0rJcYAuH8yTiUhi0rkD97VS1vqNmG0VqvAgYppVoA7TES0i6t\n9cHaBhCo0prFcevlivHDO7Ji8yHS12azRueyRufSomms0fX4vJbS9ThIFd1xN0W33XmmzcWcW6xo\n/A2B3xAoRIDyZRxKOvC81vpHrfUh4JDH+quByVrrHn6O0VLRkeEM751mdD3en8/8tdms2prD9B+3\n8cXCLAZI1+Pg5aXNpbRvP6+bhq9aYUxeKdViQlSpyoSilIoFmplPbcBwYKZSapuXzcOA0UBAzDRc\nH2w2G53SEumUlshNl3ZmycaKXY/bt0xgRN80BnSTrsdB6Sy3Obbl5Bi3OU5o0uDTvwgRTM52hRIP\nrAOS3Ja9bv6ryg/+CCrQJcRGMnpgO0YNaMvmnUeZvzab9VmHmfr1Vj75cTvDerXi4j5ptGwaa3Wo\nwg9sjjIK732gvFrMNc6l6NY7OPWn560OT4iAUWVC0VrnKKVuBVw/254BZgIbvWxeBuQA0/0eYQCr\n3PU4m0Xr9jNv1V7mrdpL9/bJjB3emY4t4gizS9fjYOVo2YpTz02uMP1L1JezsRUWWB2aEAHlrG0o\nWutvgG8AlFLtgX8GUk+uQJKSGM11F3XimqEdWKtzSV+bzeZdeWzetYqk+MjydpjkhCirQxW15TG3\nGAXeE0rEwnRsp0ukWkw0OjVulNda36WUaqeUmgy8pLXOA1BKPQk0N5fl1FOcQSM8zM6Abi0Y0K0F\n2bknWbE1lx9W7WH2kp18mbGLPqoZI/r41vVYBCC7HeK9TxwR98pkIlYsqzzlvhAhzpeBjT2BBRht\nKv201uvN5S8BvwGOYgxs3FkPcfqsNgMb60thcSkL1+7j66U72bnfuG1MWmo8o4e059ILziE+Vn7F\nhpRly+DTT+Gzz2DfPmNZUhKsWgWdO1sbmxDVaJCBjUqpr4AewEit9TaPdR2B+cAqrfX1tQ3Gn+o6\nsNFf3AcOOZ1Osvbnk752H6u25lBa5iQy3M6A7i24pG8a7VuevetxKA6usmpgY03UOQaHo3ycS8Tq\nlUSsXEHukVP1GocMbKyfOELxs1fV8oYa2DgIeM4zmQBorXcopd4EnqxtII2BzWajc1oindMSufHS\nLmRsOEB6ZjZLNhxgyYYDdGiVwIg+bRjQrTmR0vU4+LmNcwFI9dIxw75vL3EvvyBdkUVI8CWhhAFn\nm8zKVs164aZJbCSjB7Vj1MC2bNp5lHSz6/G7X2/hk/nbzFmPpetxqIv87psKXZFlnIsIZr70ZV0G\n3KeUSvZcoZRKAH4JrPBXYI2F3WbjvI4pTJzQi5ceGMyYwe0Is9uYt2ovv/v3cl6Znsman3Mpczis\nDlXUg6K7f0neV5Vvcxz30l+sDk0In/lyhfIssAjYqJT6CNgGOIHOwM1AS+Aev0fYiDRLjGH88E6M\nHdaBNT/nkp7p6nqcR3JCFKOHdKBf5xTpehxKPKd/Mce5FI+7zurIhPCZL92GVyilLgNeASZ5rF4P\n3Km1XurP4Bqr8DA7A7u3YGD3FuzLPWnOenyQj77byvR5NvqaXY+7Stfj0FKD2xwnThiLo3VrqRYT\nAcnX2YYXAwOVUs2BdhjtKns8b7gl/KdNajy3jTyX8cM7sWnvceYszGL1z7ms/jmXVimxXNzbmPU4\nVmY9Dnm2E/mEbddELkqv3OZyyeXG2BghLFSrG22bAxgrDWJUSqVqrXPrHJWoJCYqnNGD29OvU1Oy\nsvNJzzS6Hn/84zY+X5jFwO7GrMfVdT0WwcuZ0ISjazdVmnI/YlkGR1eutzo8IXxLKEqpB4FRGBNH\nuv8cCgeaAN0BuQavRzabjc5tEuncpmLX48UbDrB4wwE6tGrCiD5p0vU4VHmZct9+7Ch4q/o8edKo\nEpNqMdFAfLkfyhPAZKAYOIExtf1e8zEG2I3RviIaiHvX4592HCV97T42ZB3h3QP55V2PR/RJIzU1\nwepQRX04y5T7ALFT3iHm7Tdl+hfRYHy5Qrkbo/H9IiAV2A5cgpFI7gVeAqb6O0BRPbvNRq9OKfTq\nlMLhY4UsXL+fxevPzHrcW6UytEdLendJkVmPGxFnTAzO2NgK41wYN5aw+35NWRdldXgiBPny7dIe\n+J/W+oTWegeQB1yotS7VWr8DzAWeq4cYhQ+aJRldj195eCj3X9MD1SaRdTqXf8zcyBPvLGPOkp0c\nO1lsdZiiARQ+8CuOZm6uMM6F996DkhKrQxMhype5vPKB32qt/2s+XwEs0Vo/Zj6/D3hRa51SL5H6\nKJAmh7Ta7gP5fL10J+lr9lFYXEqY3cagnq24cmh7zuvUTLoeNxYOB6xZAxdc4L3NZf58GDZM2lwa\nubpMDulLlddWYCjwX/P5z8AFbuuTgIDquyoT1BnatWrChIs6MmZgW5ZvPkT62mwyNuwnY8N+o+tx\nnzSG9jS6HsvkkPUXQyBMMJjav7/X9WFbt9D00kshKYmiKypP/yKTQ1p/7qrbxl+TQ9aFLwnlXeBt\npVQUcD8wG5ihlPoTsAV4BNjgt8iE38VEhTOiTxoX927N9uzjpGdms3prDh//YHQ9HtS9BddeokiM\nkt5hjY0zNpaC+x8idu6cCm0uhfc9SMHjT1sdnggSVSYUs4vwD67ZhbXW/1RKtQF+BZRg3A74K4xb\nAwPkI7MNBwWbzUaXNkl0aZPETZd2YcmGAyzIzGbR+gMsWn+Ajq2Nrsf9u0rX48bC0bYdp56bTOzb\nb5L3zXxjnMucWVaHJYLM2a5Q/gpMxJizC6XUTowbaaVqrU+by8YCFwIpQIbcsTH4NImN5MpB7bhi\nQFt+2nmEjE2HWL35EDv25zP9x20M62XMetwiWWY9bhQ8xrlQ7L0DR+RXcyA6CsZf08ABikB2toRS\nDIwzG99PYky10h5opVSFLoe7zH/RSqm2Wus99RKpqFd2u41enZpx6aAObNmWw8L1+1m0fj/frdzL\ndyv30qNDU0b0SeP8zgHR50I0BLsdYrzckcLpJO75PxK+IwuSkkjw0uYiGqezJZT/AI8DV7kte938\nVxUnxvxeIoi5uh5fM7QDa3QO6Wuz2bTzKJt2HiU5IYorhxqzHifFy6zHjdWJN/9J1JxZxM6dXaHN\n5ejyTJCOyKoCAAAerklEQVSBtI1WlQlFa/2kUmox0AtjOpVnMNpNNp6lPOmqG0Iiwu0M6t6SQd1b\nsi/nJOmZ2SzddJAPv93KdLuNPirVmPW4bZJ0PW5MbLbyWZFj336DvG/nEzVnJmE7snCmyBVsY3bW\nXl5a668wGt5RSt2FMbBxdgPEJQJMm+bx3D7qXCZc3IlNe44xZ1EWq7fmsHprToWux6KRqcGU+2Eb\nNxA75R2pFmsEfLkfSvt6jEMEiZio8PIbfW3bd5wFmdnGrMdm1+OL+57D4G7NaddSqj2EIeq7r+U2\nx41EraavF8Jms6HOSUKdY3Q9XrxhPwvX7Wfeit3MW7G7vOvxgG7NiQiXZrXGrODRJyi5aER5V2RX\ncjn5x+cpfHii1eEJP5KEIuqsSVwkYwa3Z/TAduw5UsCsBdvZmHWkUtdjmfW4kfIy5X7UlzMpvqqK\nLsdOp/epYUTAk4Qi/MZut9G/e0vap8aRe6yQhesqdj3uo1IZ1rMlvTrLrMeNllty8crhIPmSYZT2\nOl+qxYKQJBRRL1KTYphwcSfGDuvAmp9zmJ+ZTabOJVPnkpwQxcW9W3PR+a3lqkVUYD+wH9uxvMpt\nLmOvpeTSkVaHJ6ohCUXUq4hwO4N6tGRQj5acPO1g5vxtLN10kJmLdzInYxeDz2vFkO4tONfsevzK\n9Ey27MoDG3Rrl8ykm/pY/RJEA3KktfF6m+Ow7dskoQQBSSiiwXRonVje9Xj5poPMz8xmyfr9LFlv\nzHrsdDo5eLTQ2NgJm3fl8dg/Mpg4vpf0GmtMvLS52EpPe93UdiwPZ2ycVIsFiKBJKEqpS4EbgVjg\nZa21zGwcpGKiwhnRtw0X90kj9+RpZqZvY/XWHMoclcfF5p0o5o3PN/Dqw0MtiFRYrrrbHL/8AtEz\nPpGuyAEimFpGY7TW92Hct16ufUOAzWajR8cU7r+mhyQMUSvO5i1wxsURPf1DEm+5npQenUn49QPY\nd++yOrRGKWgSitb6K6VUHMYMyNMsDkf4WZO4SLq3T/a6rnfnZhSXlDVwRCIYFDwyiaNrN525zXFs\nLFGffgxRMs+cFQIioSilBiql0s3/25VS/1RKLVVKpSulOpnLmwFvAs9orQ9bGa+oH5Nu6kNywpkv\ngpiocOJjIkjPzOapfy9j4bpsyhwOCyMUAcmsFjv13GSOZm4mb34GjpatKm9XVkbkj/OgpKThY2wk\nLE8oSqkngCmA65tkHBCptR4CPAW8ai5/FWgBvKiUGt/ggYoGMXF8L5ITokhJjOaJm/vw0gODuWpI\newqLS3nv25955r8rydS5OJ0yD6nwwm6nrEdPr6silmWQePMEo1ps4oNE/vCdJBc/C4RG+e3AdcD7\n5vNhwLcAWusVSqkLzP/faU14oiG1a5nAqw8PrXCv6+su6sglfdOYvWQni9cf4M0vNtKlTSLXj+hM\n57REiyMWwcLRujUF9z9U3hXZNc6l4JFJMgWMnwTE/AZKqfbAx1rrwUqpKcDnWutvzXW7gQ5aa5/q\nOpzyEzYk7T10gv99vZnlPx0EYPB5rbjjym60aS7dikUNORywfDnMmAGffQa//z088IDVUQUMWx3u\nRREIVyie8gH3bwe7r8nExfUL10ruv7StKsuX/Wqy7dm28XVdVdtXtTzaDvdd1Z0RvVszIz2LZRsP\nsOKng1zUuzVjh7Yn0Y83/fLnuatLef48f7Vd78tyf79vtXXWOLqcB787D576E5SWgpftoqZ/iLNZ\nM0ouGkFqWorl5666bWqzzt/nKhATSgZwNTBDKTUIkPEmopIubZJ4+ra+ZG47zGcLsliQmc2ynw4y\nasA5jBrQlpioQPzTFgHHbvc+buX0aeL/9HvsR4/iSEyCa8cROXKMjHOphuWN8m5cVVQzgSKlVAZG\nQ/xvrQtJBDKbzUZflcpzvxzAHaPOJToyjDkZu3j6X8v4cc0+SsukR5iopbAwjv/vE6MrclwcTJtm\njHPp3RWKiqyOLmAFRBtKfZA2lManqLiUWYuy+CJ9G4XFZbRqFscdV3ZjaK/WcotiUXsOB6xYAZ9+\nCseOwdSpVkdUr+rShhKynzKn0+kM+HrcBiorkOtx66NuN/9UCV9m7GLBumzKHE46tGrCDSM6cW5b\n7wMnqyJtKLVbHhRtKH4uK2LRAqI/+6TS9C+B/Nmrannz5k1qnRcCqcpLCL9oEhfJrSMVz987kP5d\nm7PzQD4vfZTJ6zPWsy/3pNXhiRAUOe+bitO/NNJxLpJQRMhqkRzLg+N68n93XkDXtklsyDrCH99d\nybtzt3A0X+rBhf+c+vOLZ6Z/cZtbjI8/tjq0BhXSVV5WxyACh9PpZM3WHKZ9tYndB08QGW7n6gs7\nMuFSRXxMhNXhiVDiPs7lmWcg2UtVawDf5ljaULyQNpTa7RcKbShn43A4WfrTQWYu3kHeiWLiosO5\nekh7RvRtQ0R4xQt2aUOp3fLG2Ibi035FRTQdegGnh15Yoc1F2lCECDJ2u41hvVrx4n2DuP7iTjic\nMH3+dn4/ZTnLNh3EIRe2op6F7doJpaWVptxn3jyrQ6szSSiiUYqMCGP0oHa89MBgRvY/h2Mni5ny\n5Wb+PHUVm3YetTo8EcLKunYzZkV2b3P55CN4+22rQ6uzkK7ysjoGETwOHS3gg2+3sHDtPpxO6K1S\nuWtMdzq1SbI6NBHqXONcIiLgggsqrz90yGiHaaAR+tKG4oW0odRuv1BvQ6nOnkMnmLEgq/wqZVCP\nFlx3YUeaJcXUuexAOH/ShtKwZfnj3CXcdxeR6fMpGT2G6NtvIff8gV6Ti7ShCBFg2rZI4LEbe/PY\nTb3pmJbI8k2H+N2U5Uz/cRsnC09bHZ5ohMo6dCzvisyYMeXjXOyHDlodWiWSUITwokf7pvztkeHc\nd3V3kuKjmLdqL0/+cxlzl+2i5LTcjlg0nIKnnym/zTGPPIIzLo6o2V/giA+8WzZIQhGiCna7jUE9\nWvKXewdx06VdsNvg84U7ePrfy1m8fj8OhzTTiQZi3uaYv/3NSC4/LoG4uMrbnTxp6Qj9kG5DsToG\nEVpOFp7mi/RtzF6YRUmpg7YtE7hzTHf6d2shk0+KwDB9Otx8MyQlwdixcMMNcNllPjXoS6O8F9Io\nX7v9GnujfE1iyDtRzKzFO1iy8QBOJ6hzkrh+RCc6tT777YgD4fxJo3zDltXQn72jS1YR/cF7RH05\ni7D92QA4EpOwT36R3PG31qg8aZQXogElJ0Rx95Xd+PM9Azi/Uwp67zH+8r81vD1zI4eOFlgdnmjE\nys7tyqnnXjSqxeaeGedC8+YNcny5rZ0QtZSWGs9vrj+fn/fk8Wl6Fqt/ziVz22GG927NNUM70CRO\n7uwnLGK3U9p/IKX9B3Lq2RdIbRYP3n7svPEGkamt/XYnSkkoQtTRuW2T+cMd/Vjzcy6fL8xi/tps\nMn46yOgBbRk54ByiI+VjJixkt0NYWKXFthP58MQTJBYX40hMouSKKykee22dDiV/6UL4gc1m44Ku\nzendpRmL1u9nzpKdzFqyk/mZ2Ywd1oELe7WyOkQhKnDGxUN6OgXvfUDUl7OJ/uQjwjesr1OZId0o\nb3UMovEqKDrNrIVZzFywnaKSMtJS47jjyu4MPq+V9AgTgcc1/cuJE9hGjZJeXp6kl1ft9pNeXv6N\n4fjJYuZk7GLhuv04nE46pTXh+os7o86p+Rxh0surduSzV/066eUlRBBJjI/i9lHn8vy9AxnSqxVZ\n2flM/nAtb3y2gf2HT1kdnhB+JW0oQjSAlk1jefrOASxbt48Z6dtZt/0w67MOc2GvVowd1pHkhCir\nQxSiziShCNGAOqcl8tStfVm3/TCfLchi0foDLN90iJEDzmH0wHbERMlHUgQv+esVooHZbDb6dEml\nV6cUMjYeZNbiHXy1dDcLMvdz9dD2jOiTRniY1EaL4CN/tUJYJMxu56LzW/Pi/YO57qKOlDkcfPzD\nNn4/ZTkrNh+S2xGLoBPSvbysjkEIXxw/WcynP2i+XrqT0jInndskctdVPTi/S6rVoYlGRCaH9EK6\nDdduP+k2XH8x1LS8nGOFzFy0gxWbDwHQt2tzxg5pzznN4+t8DOk23LBlBfJnr6rldek2LG0oQgSY\n5kkx3H9ND0YNOIcZ6Vms3ZpD5tYcBvdsybUXdiQlMdrqEIXwShKKEAGqfcsmTLqpN/uOFvGf2RtZ\n+tNBVm7J4bJ+bRgzpB1x0RFWhyhEBZJQhAhgNpuNvl2b88em/Vm+6SAzF+3g25V7WLR+P2OGtOOy\nfm2ICK888Z8QVpCEIkQQsNtsDOnZiv5dm/PjmmzmLtvFjPQsflyzj2sv7MjgHi2x20O2SVQECek2\nLEQQiQgP44qBbZn8wGCuGNiW/FOn+e/cLfxp6io2ZB1BOjcKK0lCESIIxUVHcMOIzrx43yCG9mxJ\ndu5JXp+xnj/8cyk7D+RbHZ5opKTKS4gglpIYzS+u6s7IAW35bEEWG7YfZsP2wwzo1pzrhneieVKM\n1SGKRiRkK11lYKNojDZsz2XqV5vZvvcY4WE2Rg/pwI2XKRLjZfJJUTMysNELGdhYu/1kYGP9xdBQ\n5+9QTj6rt+bw+cIsco8VER0ZxuhB7Rh5wTm0SUuSgY0NWFYgf/aqWi4DG4UQ5ew2GwO6taCvSmVB\nZjZzMnYxc9EO5q/dx+2ju3F+h2TC7NJ8KvxPEooQISo8zM5lF5zD0PNa8c2KPcxbuYe3ZqynVUos\nE4Z3oneXZnI7YuFXklCECHExUeFcd1FHRvRJY96afcxbsZs3v9hIlzaJXD+iM53TEq0OUYQIue4V\nopFITojiV9f35rlfDKRPl2Zs23ecF95fw1tfbOTAEbkdsag7uUIRopFp3SyOX4/vxbZ9x/g0fTtr\ndS7rth3mot6tueeanlaHJ4KYJBQhGqkubZL43W39WKsP8/nCLBZkZrN800FG9j+HUQPayu2Ihc/k\nL0aIRsxms9Hv3FR6d0lh8foDfLl0F3MydrEgM5urh3ZgeO/WcjtiUWPylyKEIMxu5+I+afzr6csY\nd2EHiksdfPi95g//WcGqrTkyR5ioEUkoQohyMVHhXDO0Ay/dP5hL+7bhyPEi3pn1E5PeWMTPe/Ks\nDk8EOEkoQohKmsRFcutIxfP3DuSCrs3Re47x0keZvD5jPftyT1odnghQ0oYihKhSi+RYHhrXk7zC\nUv79xQY2ZB1h444jDO3ZinEXdiA1NcHqEEUACdlhsjI5pBD+5XQ6Wb3lEO/N3czugyeIDLdz9YUd\nmXCpIj5GbkccKmRySC9kcsja7SeTQ9ZfDIFw/mq73n25w+Ek46cDzFq8k7wTxcRFh3P1kPaszzrC\n1t15YINu7ZKZdFOfGsVcX+SzV/06f08OKW0oQgif2O02LuzVmhfvG8SdY7rjcML0+dvZsjsPJ+B0\nwuZdeTz2jwx2H7T+R51oOJJQhBC1EhkRxoRLuvDSA4O9rs87UczLH2eyIeswB48WUFrmaOAIRUOT\nRnkhRJ3Ex0RgA7w1WhYWl/L6jA0A2GyQ0iSa5skxNE+KoXlybPn/U5NjiIoIa9C4hf9JQhFC1Fm3\n9sls3lVxnEpCbAQj+5+Dwwk5eQXk5hVy6Fghm3flsZnKY1oS4yNpYSaa1OQYWiTHkJpkPMZGS6N/\nMJCEIoSos0k39eGxf2SQd6IYMGY2fvXhoV63LS4pI/dYIYfyCsk9VkhOXgE5xwrJyStkW/Zx9L7j\nlfaJj4koTy6pSTE0T46hhZl4msRGyH1dAoQkFCGEX0wc34s3Pt+A3W7jV9eeV+V2UZFhtGkeT5vm\n8ZXWlZY5OHy8yEgyeUaScSWbPYdOsPNAvtfymptJxr06zRkehtPpxC7JpsFIQhFC+EW7lgm8+vDQ\nOnXXDQ+z07JpLC2bxlZa53A4OXqiqFKiyckr5FBeAXtzKo/gDw+zk5oUXbHNxvyX0iRaJr70M0ko\nQoigYLfbaJYYQ7PEGLq3r7jO6XRy/FRJhWSTX3CaPQfzyckr5MCRAuBIxfJsNlISo8qTjatKrblZ\nrRYpnQR8JglFCBH0bDYbSfFRJMVHoc5JAioO2jtZeNpstzGq0lwdBHLzCtm0K49Nuyp3EkhOiCpv\nr+nYJonYCLvRbpMUQ2y0fHV6I++KECLkxcdEEB8TQYdWTSqtKyopNZLMMc+qtAK27T2G3nuMJRsO\nVCqvRbLR3flM+41RpZYQ03g7CUhCEUI0atGR4bRtkUDbFpUnujxd6uDw8UKKHbBt11G3hFPAroMn\nyNpfuZNAdGRYeeeA9mlJxEed6TSQlBAV0p0EJKEIIUQVIsLttEqJIzU1gfapcRXWlTkcHM0vrpBk\nXP8/eKSAPYdOsvrn3ErlpSbFeOmVFkPTphXLD0aSUIQQohbC7EZySE2KoYfHOqfTybGTJZQ4Qe86\nUqln2v7Dp7yUZyufSSA1OabCIM/mSdFEhAd+JwFJKEII4Wc2m81o1E9NoEWTqArrnE4np4pKK3YQ\nyCvk2KkSsnNO8NPOo7DTozwgKSGKNs0TSI53DfKMLe80EBMVGF/lgRGFD5RSlwA3a63vtToWIYTw\nlc1mMzsJJNKpdWL5clevtMLiUo8OAmeq0jZmHfZaZpPYCFqnxpMcH1VhrE3zpJgGvVdNUCUUpVQn\noDcQbXUsQghRH2Kiqu4kkJgUy5btueW90HLcEo/eewyHo/IUnTFR4bROjaOpK9m49UpLSTkzW8Er\n0zPZsiuPqx6d5fjqtXG1GvEZVAlFa50FvKaUet/qWIQQoqFFRoTRulkcrZtVbsBPbhrHz1m5FWYQ\ncCWbvQdPkFVaeY60yIgwUhOjyT9VwonC067Fte6GFjAJRSk1EJistR6hlLIDbwO9gGLgl2YyEUII\n4UV4mN0cCxMLHSquS0mJZ9vOw+WTcroSzdETxezPPUlRSZl/YvBLKXWklHoCuA1wTcYzDojUWg8x\nE82r5jIhhBA+stttNG0STdMm0ZzbNrl8eWpqAjk5+fzypXSv97Px+Th+KMMftgPXceZSaxjwLYDW\negVwgfvGWuvbGzQ6IYQIUTabjW7tk6vfsAYCIqForb8ASt0WJQDuQ1DLzGowIYQQfjbppj4kJ0RV\nv2E1AqLKy4t8jKTiYtda+3RDaltjnUxHCCFq4apHZ/UF5tSljEBNKBnA1cAMpdQgYIPF8QghREj7\n6rVxa4E2dSkj0BKKq11oJnC5UirDfH63RfEIIYQQQgghhBBCCCGEEEIIIYQQtdaoutbKTMXBSSl1\nKXAjEAu8rLWWXn9BRCnVD/gVxvfNE1rrHItDEj5QSrUAvtJa969u20YzWFBmKg5qMVrr+4BXgJFW\nByN8FgU8AswFBlsci/CBUsoGPA7sqsn2jSahaK2ztNavWR2H8J3W+iulVBwwEZhmcTjCR1rrpUB3\nYBKwzuJwhG8eAD4AimqycaCNQ6kVmak4eNXk3CmlmgEvA89orb3fYUhYoobnrz+wGhgN/BH4jWUB\ni3I1/N68zFw2QCk1Xmv9+dnKDPorFHOm4ikYl9XgNlMx8BTGTMUiAPlw7l4FWgAvKqXGN3igwisf\nzl888C7wV+DDho5TVFbTc6e1Hq+1fhBYUV0ygdC4QnHNVOy66VaFmYqVUjJTceCq0bnTWt9pTXii\nGjU9f+lAuiURiqr4+r15R00KDforFJmpOHjJuQtucv6CV32du1A82XWeqVhYRs5dcJPzF7z8cu5C\nMaFkAFcCyEzFQUfOXXCT8xe8/HLuQqENxUVmKg5ecu6Cm5y/4CXnTgghhBBCCCGEEEIIIYQQQggh\nhBBCCCGEEEIIIYQQQgghhBBCiICglHIopaZaHcfZKKWilFLvKqXylVLHlVJjAiCmjh7PFyildvpY\nhs/7iNARSlOvCOHOWf0mlroXuAv4H7AIWGNlMEqpu4F/ALFui5/3eF5Tgf7ei3oiCUUIa/QyHx/W\nWp+yNBLDcM7cbAkArfUPFsUiglQozjYsRDCIBAiQZOJiszoAEdzkCkX4hVJqF/ANxjTYTwMdgb3A\n61rrtz2226m1HuFl//Ll5vMvgfXAE0Ab4CfgYbPcN4ArMO7j8B7wB621e1WLTSn1O3P7JGA58KTW\nerXHca8Cfgecj3Ev7fnA01rrbW7bODCqf3oDIzHudne+1rqsivdiLPCkuX0xRpXWH7TWG93Kcy97\noef74bH+QaADcA9GIlpivpZNbttFAJOAm4DOGMlBA3/XWk/1KM/zteQBQ93WT9Na36OUWgC001p3\ncNu/K/Bn4BKM749M4P+01ku8xW/u0x34C3CxGX8m8Get9Ty3baKAl4BrgNZADjDHfN+OVVW2CCxy\nhSL8xQmMBv4OfAo8ApwC3lJKjfbYzlsdu+dyJ8Z9rp8F/m0+dgU+B77HuNvcoxhJ5mnA89bOE4Df\nAm9jfAF2AxaYX24AKKXuAmYDJ4DHgdeAwcAKpVQXj/J+C0QAvwamnCWZPIwxFXiYGddrwEBgqdtt\nVW8HFpv/vw3jC/5sngbuMMt6GRgALFFKdXbbZirGe5Ruxvgsxr3c/+vx/ld6Lea27vH8y23b8nNi\nvicrMBLDGxiJuCnwvectY932OQ9YhnHu/gL83jz210qpG9w2fQv4JfARRgL9DLgP+KTKd0UEHLlC\nEf5iw7iK6K21/glAKTUL2A/cinH14tquqv09n7fCuBLYZJbXFOOLf4nW+hZz2UfAUeByjAZulyhg\nkNu+nwFbMJLLBKVUE4zkN11rfatrJ6XUFGAzxq/l69zKKwHGaa2Lq3oDlFIpGF/4K4ALtdal5vL/\nAZswGr0Haq0/VEpdbm7zUVXluUkFumutd5vlzcK4cvszcItSqiVwMzBZa/17t3hmAluBUZx5/72+\nFqXUbTWI53mMRDlIa73D3G86kMWZqyNPbwKHgL5a60JznzcxrgT/rpT6wnyfbgX+o7X+g1tMJ4FR\nSqlYrXVBdW+SsJ5coQh/+tmVTAC01ocwvkxa1LK8LPdqHcBVDTXT7RgFQC5G8nH3rfu+Wuss4GuM\nLyg7RgJKAGYrpZq5/gFlGL/yXdu5rDxbMjFdCsQAr7qSiXns3cD7QH+lVG3eiw9dycQsbzPwLTDG\nfH7QfC3lVzpKKRtmOw3GlYq7mryWCsz34krga1cyMY99FBgGTPSyTwpwEcb7Huf2HicDszD+Lvqb\nm+8FblJK3amUSjLLfkZrPVCSSfCQKxThT7lelpVg/KqtjUMez11f0jkey8uo/ONoq5fydmDU0acC\nncxl06s4ttPczhWD5zG9cbU1/OxlnSuedlR+XdXZ5GXZNmCMUipFa30EOA3cppQaBSiM1+e6R7jn\ne1OT1+IpBYjjTFIv55H03bne44l4STgY73FbjCqxBzGqSqcCpUqpZRg/HN7VWufXIl5hAUkowp8c\n1W9SJW9Jp9TLMqjZOAdv27iq1UrdjncvUNVAPPfGYK9tJlWU743rS72kBuV4Ou1lmSv+UqVUNEYb\nSG+MqqR5GPcEXwTs8bJvTV5LVcfzZYyJa5+3MK5IvNkMoLWer5RqC1wNXIXRYeA14LdKqX5a68O+\nhywamiQU0dDK8BjvoJQKB5rh5ddvHbT3skwBx7TWR8xeZACHtdbzPeK5EAjztVqIM4mpG7DRY925\n5uM+H8sEo9eWpy5Artb6uFLqDqAfcI/WepprA6VU61ocqyqHgUJvsSilJgEttdaTPFbtMh/LvLzH\nXTGu6ArMHmq9gWyt9SfAJ2aV3aPAXzHaZt7y42sR9UTaUERDOwh0NX9Vu1yDR5LxgyuVUmmuJ0qp\nnhiN03PMRd8DRcDjZkJzbdcKo7vyS7U4pqvMR80vSVeZbTB6T63w+KVd01/7tyulkt3KO998LZ+Z\ni1LMxy0e+/3GfKzJD8cys2yvV1lmm9A8jPe1jVssyRgdJTp42ecAsBq4y3xfXfuEA//F6LEXhtFT\nbDlGbzbXvk5zX6j6SlUEGLlCEQ3tI4yeP98qpT7E+MV7L7Ab/w6sKwIWK6X+jtGW8AhwBPgDgNb6\nsDlO5TVgmdlbLAx4CCO5ef7arpbW+qhbmRlmmQlmmVC5HaGmrzcOWKmUesfttewH/mSun4fxpfu+\nUuot8/9XA30w2rWa1OAYrnaVZ5VS6VrrdC8xPo3Rg22leZwTGOcuFvN99bLPRIxquDVm/IcxrjgG\nA09prfMAlFLvAQ8ppeIw2lRSgF9h/AD5tAbxiwAgVyjCX6r6te25/G3gjxi/aN/A6AU0DmM8iec4\nlJoex9uyfwEfYoyVeApjwOUQrXV5lZPW+nXgBowv4OfN7TRwidZ6caUSa8As80YzphcwvlCXYHQX\nXuURc02vUN7CHOSH8SU7Fxistc41j7kJGA+cBF4EnsGoWusN/AgMVUpV1zHiHWAVxiDSx73FqLXe\nipEIVprbPQtkA8O01luq2Gc5xqDJ1ZypwooD7tRav+x2/AeB54AhGN25H8NoFxpm9iQTQghRF+bM\nyS9YHYcQNSFXKEIIIfxCEooQQgi/kIQihBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghauL/\nASPWsLjy2HRSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xadeda66c>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(N, pyt/fort, 'o-')\n",
    "plt.plot(N, 5000./np.sqrt(N), 'r--')\n",
    "plt.title('python/fortran', fontsize=20)\n",
    "plt.xlabel('number of particles', fontsize=18)\n",
    "plt.text(500,700, r'$N^{-1/2}$', color='r', fontsize=18)\n",
    "plt.ylabel('factor', fontsize=18)\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}