examples/WendyScaling.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import numpy\n",
    "import time\n",
    "import wendy\n",
    "%pylab inline\n",
    "import matplotlib.animation as animation\n",
    "from matplotlib import cm\n",
    "from IPython.display import HTML\n",
    "rcParams.update({'axes.labelsize': 17.,\n",
    "              'font.size': 12.,\n",
    "              'legend.fontsize': 17.,\n",
    "              'xtick.labelsize':15.,\n",
    "              'ytick.labelsize':15.,\n",
    "              'text.usetex': False,\n",
    "              'figure.figsize': [5,5],\n",
    "              'xtick.major.size' : 4,\n",
    "              'ytick.major.size' : 4,\n",
    "              'xtick.minor.size' : 2,\n",
    "              'ytick.minor.size' : 2,\n",
    "              'legend.numpoints':1})\n",
    "_SAVE_GIFS= True\n",
    "numpy.random.seed(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scaling of ``wendy`` with particle number *N*\n",
    "\n",
    "We will investigate how ``wendy`` scales with the number $N$ of particles using a self-gravitating disk. The following function initializes a self-gravitating disk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def initialize_selfgravitating_disk(N):\n",
    "    totmass= 1.\n",
    "    sigma= 1.\n",
    "    zh= sigma**2./totmass # twopiG = 1. in our units\n",
    "    tdyn= zh/sigma\n",
    "    x= numpy.arctanh(2.*numpy.random.uniform(size=N)-1)*zh*2.\n",
    "    v= numpy.random.normal(size=N)*sigma\n",
    "    v-= numpy.mean(v)\n",
    "    m= numpy.ones_like(x)/N\n",
    "    return (x,v,m,tdyn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time/collision for exact solver\n",
    "\n",
    "We believe that this should go as $\\log(N)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ns= numpy.round(10.**numpy.linspace(1.,5.,11)).astype('int')\n",
    "ntrials= 3\n",
    "T= numpy.empty((len(Ns),ntrials))\n",
    "ncoll= numpy.empty((len(Ns),ntrials))\n",
    "tdyn_fac_norm= 3000\n",
    "for ii,N in enumerate(Ns):\n",
    "    for jj in range(ntrials):\n",
    "        x,v,m,tdyn= initialize_selfgravitating_disk(N)\n",
    "        g= wendy.nbody(x,v,m,tdyn*(tdyn_fac_norm/N)**2.,maxcoll=10000000,full_output=True)\n",
    "        tx,tv,tncoll, time_elapsed= next(g)\n",
    "        if tncoll > 0:\n",
    "            T[ii,jj]= time_elapsed / tncoll\n",
    "            ncoll[ii,jj]= tncoll*(N/tdyn_fac_norm)**2.\n",
    "        else:\n",
    "            T[ii,jj]= numpy.nan\n",
    "            ncoll[ii,jj]= numpy.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x106ed02e8>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAFLCAYAAABSqd7ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8XHWd//HXp7RAKbTTlCK2UGAKFBAEJumCyEVgQsEb\nrCatuPtTf0ITXHbX9bfbxupe1HW3pHVd113UBFbXy4o04SJ4q5kiXipCm1RAkW5JCkILizSZQiGU\nNvn+/jhnQm4zmZPMzJnL+/l45PHIzJw58/026TvnfK/mnENERLIzLewCiIiUEoWmiEgACk0RkQAU\nmiIiASg0RUQCUGiKiAQQSmiaWYP/FTGzqJk1h1EOEZGgwrrSjAAtQB/Q4X8vIlL0pof0uUlgLoBz\nLhlSGUREAgsrNBWWIlKSQgtNM2sAeoGlwO3Oua6wyiIiki0LY+65mUWdcz3DHncD1aOvPv1gbUg9\nnjVrVvVpp51WuIKKSEXo7Ox8wTk3P5tjQwnNMYUw6wRanHOtmY6rqalxW7duLVCpRKRSmFmnc64m\nm2ML3nvuDzHqG/V0D7C40GUREQkqrCFHTaMeR4DuMAoiIhJEwUPTb8uMpB6bWQSITnRrLiJSDMLq\nPW81s9X+94uB2pDKISISSCih6feSrwvjs0VEpkILdoiIBKDQFBEJQKEpIhKAQlNEJACFpohIAApN\nEZEAFJoiIgGEtjScSD61t7cTiUTo6ekhHo8TjUazPqa9vX3omGg0SiwWK1i5pfgpNKXs9PT0sGXL\nFpqbva2n6uvraWtry+qYnh5vxcK6ujoAmpqaFJoygm7Ppey0t7ezePHri2Z1dY1d3zrdMVVVVaxd\nu5aenh6SyeSIY0RAoTlCMpmkurqaZLJ4duJYt24d7e3ttLa20tqa3Zom7e3tNDWNXkhq8ucrBolE\ngvr6esyMRCIx4rVUAFZXV5NIJNizZw9VVVUjjhn9M013TCQSobGxkcWLF9PU1ERDQ8OIY7q6umhq\namLduvzOAg5SXyksheYwqV/AtWvXTnhsIa5AmpqaiEaj1NXV0dDQQHd394j2ttESiQTr1q2jpaVl\n3OAPer5iEo/Hqa2tpaGhYeiWOqWuro7GxkY6OzuJx+MA9Pb2TnjOdMd0dnbS1tbGhg0bxvzxicVi\nLF68mI6OjknWJDtB6yuFo9D0pa4ympubWbduXcarzUQiMW7HQq61trYOta0B1NbW0tKSfrfjeDzO\n6tWr07bBBT1fMWpsbCSRSIz5+UQiQ6sNMm/evBGv9fb2jng90zGtra00NjZSV1fHzp07x/2jUoif\nfUo29ZXCUmj6EokE8XiceDxOLBbLeLXZ1tZGbe3Eq9mtW7eO2tpampqahjoYRhuvvS3d81VVVZO+\nHcv1+YJob29P+0coaBNBLBYb8/NJJpMjbrXr6uro7Owcem341ViqHJmOSQVS6lY9TNnUVwpLocnr\nV5kp6a42U7e/ra2t7NmzJ2O7VlNTE3v27KG5uZmlS5dSX18/5lavq6uLdHse9fb2jvmPkSrjZNpc\ng56vp6dnzB+G4e2969atI5FIDLWfpgt/8K6AV65cOeZzGhsbs/7P39XVRU1NzdD7hoft1q1bR4Re\nNBodau9rbW0dcXtbXV2d8ZiGhoYRbb7Z3P4mk8mhtuL29vYxvxc9PT00NTUN/Vulzp3pD0aQ+kqB\nOedK5qu6utrlQ0tLy5jnotGoW7169Zjn+/r6XCQSmfCc4723paXFRaNRV1dX5+LxuIvFYmnf39bW\nNuZz+vr6HOC6u7sn/OyGhoYpna+5uXnEObq7u4fe39LS4tra2oZe6+jocJ2dnRnL1NfX5+LxuOvr\n63POOdfQ0DDiHBMZ/TMCXEdHx7iv5VtHR4eLx+NDj2Ox2FC9Uq8P/7eLxWJD/z6j35tOMdW3EgBb\nXZY5VPFXmslkctw2qkxXm6krgEzWrFkz5rlU50tjYyNNTU1Dt4fjGa/NKtVxMZlbs6Dn6+joGHGl\nmWq+AO8qLXXFlBoYPtFYxkgkQltbG/X19TQ2NlJbWzuifTWo8TpIwpBq3hj+7xuPx2ltbR363enq\n6hr6HYtGo2nvLjIplvqKbs/ZsGHDuLc6dXV1RKPRMW2bo8Mkk/r6eubOnUt1dfWIDoVU2ymQtve6\nqqpqTGCnHk+mEyDo+YaHJIysdzwep6WlhY6OjqGhL9k0GUQikaHQmOrtZaqDpKenJ9T2va6urnE/\nPxKJDIVjXV3dULgmEgmWL18e+HOKpb5S4aE5UYP6eFebw8MkXecOeG2ajY2N9PX10dbWxpYtW6iu\nrh7R9pepTSsWi40Js97e3kmHTZDzpa6Mhh+fqndPT8/Q921tbTjnhq6sJtLY2Eh9fT2bNm2ivr4+\n67bZ4e17w+sTjUapr68PtX0vGo2OO3Rp+B3MihUr6O3tpb29naqqqglHLBRzfaXCQzPV0J9MJsf9\nSs1HHt6B09vbSywWI5lMZgzN6urqEbezzc3NtLW10dTUxOLFi5k7dy7d3d0Zb1FTnRIpHR0dI3pz\ne3p6Ao2znOh8KYlEYkxgpkKgq6uLrq6uEb3uK1asmPCzU4EZj8dH3KpnE5yJRGLc2/+mpqYxnXiF\nVldXN+Z3ob29fehOBWDLli00NDRQV1eXVZNEMddXqNyOoFQnSLZfqYb+1atXu7a2tkCdGFPR3Nzs\n2traXHNz85gOgJaWlhGdCp2dna65udlFo1EXiURcc3PzmA6aTOdLicViLh6Pu5aWlqFOnlTHTWdn\np2tpaRnqDEqdK5PUeUbr6+tzdXV1ad/X19fnGhoaHDCmYyv1+ngdbvnU2dnp4vG4A4bqnSpHW1ub\na2lpGfPv0dLS4iKRiItGoy4ajbp4PD7u708x1rdSEKAjyLzjS0NNTY2bTCO6BGNmlNLvRTHr6uri\n9ttvZ82aNUQikaGr0pUrV3LLLbdoMZAiYWadzrmJe3ip8NtzGSvdraFMTiKRYOnSpSMGzMdisaxG\nYEhx0pWmjJDq0Bm9UIVMXmr40fDJBKk1AKQ4BLnSVGiKSMXT7bmISJ4oNEVEAlBoiogEoNAUEQlA\noSkiEoBCU0QkAIWmiEgACk0RkQAUmiIiASg0RUQCUGiKiASg0BQRCUChKSISgEJTRCQAhaaISAAK\nTRGRAEIPTTOLm5mWsBaRkhB6aALNQPrNx0VEikiooWlmcSD95uEiIkUm7CvNCNAbchlERLI2PawP\nNrM651y7mdWGVQYRKT93b9vF+o3b2Z3sZ0FkJquWLeGacxfm7PyhhKaZRYBkGJ8tIuXr7m27WHPn\no/QfGABgV7KfNXc+CpCz4Azr9ny5cy4R0meLSJlav3H7UGCm9B8YYP3G7Tn7jIJfaZpZFMhq83Iz\nawAaUo8XLVqUr2KJSBnYnewP9PxkhHF7HgOifs85QA1QZWY451qHH+g/HnqupqbGFa6YIlJqFkRm\nsmucgFwQmZmzzyh4aDrn2oc/NrOlQMfowBQRCWrVsiUj2jQBZs44hFXLluTsM0LrPQcws9VAHO/K\ns3d0oIqIBJHq7Mln77k5Vzp3vDU1NW7r1qyaQ0VEsmZmnc65mmyODXtwu4hISVFoiogEoNAUEQlA\noSkiEoBCU0QkAIWmiEgACk0RkQAUmiIiASg0RUQCUGiKiASg0BQRCUChKSISgEJTRCQAhaaISAAK\nTRGRABSaIiIBKDRFRAJQaIqIBKDQFBEJQKEpIhKAQlNEJIDAW/ia2TlAFRDxn0oCPc65J3NYLhGR\nopRVaJrZe4FawAEGdOOFJcBiYLmZ4b/e4Zy7M/dFFREJX8bQNLNzgTjQ6Zy7IZsTmtnlZvY3QMI5\n9+sclFFEpGikbdP0ry6dc269c+6+bE/onNvknPucdwq7LBeFFBEpFmmvNJ1zd0zlxM65bVN5v4hI\nMQrUEWRmJ/rf9jrnXvQ7hd4HPOGcuzXHZRMRKTpBhxzdgNfGWeW3d7YB3wE6zez6XBdORKTYBB1y\n1OGc2wRgZs1Aa6qzx8yqcl04EZFiEzQ09wz7Pg6sHPbYTb04IiLFLejt+WIzm21mq/CGIb0IoF5y\nEakUgULT71Fv9B/WA5jZTcByIJrboomIFJ/A0yidc+tHPdWBN41yZ26KJCJSvAJdaZrZTWZ2vZmd\nnXrO7xiKqvdcRCpB0DbNLcAVwE/MbMDMNvpTJucAJ+e8dCIiRSbQ7bnfpnkHgJmdBMTwFvJoxG/j\nFBEpZ4HbNFP8NsydwB1+gJ6Us1KJiBSpnCxC7AdoZMIDRURKXNCOoCfMbIeZfdnM3jNsLjrA0pyW\nTESkCAUdp3ky3pjMHryFOrr8DqEBIDa8V11EpBxNZpzmNmBo2Tczi+JNqYzxevtmF3C7v66miEjZ\nmHRHUIpzrgdoTT02szlADV6IioiUlSmH5mjOub3AJv9LRKSsZNzuYlRHTyBmdqKZvSfNaxEzW21m\nDWbWZmbxyX6OiEghZdzuwt8kLY63SdqT2ZzQD9p6oDvDrpRrnHNN/vEJoNvM5jrnkmmOFxEpChlv\nz4ctOLzK7/BJ4q2p2cPrW/hG8FY4Otr//olxFvUYrcHMOpxzCedcj7/9bxSvA0lEpGhl1aaZCsFh\nnTxRvP3OwQvPbamAzVK134GU6n0HL4hFRIpa0LnnOenkSQWmrxFo0q25iJSCnPeeZ8u/wqzDu2pd\nG1Y5RESCCC00/avNdX54dppZ9eirTTNrABpSjxctWlTgUoqIjJSTBTuCMrOhxT388EwCa0Yf55xr\ndc7VpL7mz59fyGKKiIxR8ND0hzD1jfOSVkkSkaKXs9A0s2zbJXuAplHPRYG2XJVFRCRfArVpmtnl\nQDNjFxw2vC0vxtxij+aPy+wys9V4t+XVwErnXCJIWUREwhC0I6gWuNwfejSCv5VvVvyAVEiKSMkJ\nenveMV5gAjjnPp6D8oiIFLWgoenMbPZ4L6RbnENEpJwEvT2/AmgysyTQO+x5Ay4H0i3QISJSFoKG\nZpz0s3eiaZ4XESkbQUNzpb/dxRhmpgU3RKTsBd1YbfjeQLOHL1KcLkxFRMpJ4MHtZrbSzJ4AduLt\nRrnHzK7PfdFERIpP0H3P3wveVr7OuXnOuSq8tkxTcIpIJQh8pemcu2XU473+c5azUomIFKnA4zQz\nvLZnKgURESkFQUNz3niD2/3nFo9zvIhIWQm63cUtZrbBzByv7+lT7b3kluW8dCIiRSbwyu3OueX+\nakcx/6nmgJuqiYjk1+AgDB6E6Yfm/NST2u7CD8kRQWlmlznn7stJqUREJmv/PrirEQ6fA1ffDJbb\nPuqMoekPI0o45570H69l/BXWU3PPT8lp6UREgujdCd95P/zhcbjis3n5iImuNGuArcMe1+LNPR9v\nu13NPReR8PT8FNo+CG4Q/qQdTr48Lx+TMTSdczeMeirT3PPe8Z4XEckr5+ChW+BHH4d5J8O1t8G8\n/A3mCdqm2efPN+91zr1oZucAK4Bu59ytuS6ciEhGB/fD9/8atn0TTr0K3tMKh4+75G/OBA3NG4An\ngISZLQY2APX40ygVnCJSMPueh9v/FJ5+EC76G7j0kzAt/xvsBg3NjtTwIjNrBlqdc7/2H1flunAi\nIuPa1eUF5iu9UPc1OLNwG0cEDc3hUyXjwMphjzNNsRQRyY1H2uCeP4dZ8+G6H8Mb31zQjw8amov9\nxYYbgU7n3IvgjdHMeclERIYbHIBNn4bN/waLLoDl34Aj5xe8GEGnUd5hZqv8h/UwtHVvBG9okga3\ni0ju9SfhjuvhiQ6o+TBc2ZyX2T7ZmMw0yvWjHn8cdLUpInnywg647Vro2wnv+DwsvS7U4uRqRhB4\nbZyaESQiubOjA9qvg0OmwwfugRPfGnaJcjojSEvDiUhuOOe1XSY+BW84E679NkQWhV0qQDOCRKTY\nHOiHe/4CHm2DM66Ba74Eh84Ku1RDgnYEZdpx8iRAO1KKSFp3b9vF+o3b2Z3sZ0FkJquWLeGacxe+\nfsDeXd6CG88+DJf9HVz01zlfpWiq0oammb0H73Y8G6lVju7MRaFEpPzcvW0Xa+58lP4DAwDsSvaz\n5s5HAbzg/P2D3oD1A/3e/PElV4VZ3LQyXWleAXQDXVmeS6sciUha6zduHwrMlP4DA6zfuJ1r3Cb4\n3v+DOcfBB++FY04LqZQTyxSaLRPcjo+gNk0RyWR3sn/Mc9M5yMp9X4F7NkL0Uqj7KhxR3DOy085u\nDxKYvpOmWBYRKWMLIjNHPI7wEl+f0cyHpm+E82/01sAs8sAEtWmKSIGsWrZkqE1zif2eW2b8C8da\nH53n/jPVV94YdvGypjZNESmIVC/5lh98jU+89kX67QgeuOibXHL520MuWTBq0xSRwjjQzzW7Psc1\nB74Kx9Uwa8W3uGT2G8MuVWBpQ3OiwDSz2UBVaorlJNpARaRSPP87aP8wPP8YXPAXcNnfh7bgxlQF\nXubYzFaa2RPATqDLzPb4c9RFREZyDjr/C1ov9VZa/5M7vF0iSzQwIWBomtl7AZxzJzvn5jnnqvDa\nMk3BKSIj9Ceh7UNw70dh0Xnwkc1wSjzsUk1Z4CtN59wtox7v9Z8rrrlOIhKep7dAy0Xw+Pcg/in4\n07vgqGPDLlVOBF1PM9OWFnsyvCYilWBwEDZ/Ae77LMxZCP/3R3D80rBLlVNBQ3Oemc1ObXOR4ncK\naWk4kUr20nNwZwPs/Cm86Y/hnV+AmemW3y1dQVc5usXMNpiZA3r8p6u9l9yynJdORErDjgTc1Qiv\nvQzv+iLEPlB0qxPlymS2u1huZpcDMf+p5tS2viJSYQ6+5m129sB/wDFneNvpFvFiG7kQODR93cOD\n0sxOTI3XzIaZRYAG/+FSYK1zLtuZRyJSDPZ0wx3Xwe5tUHMdLPsnmDFz4veVuEChaWbn4LVd3mRm\ntcOC0szsMudctrtRNjvnGv03RoFOM6t2zvVM8D4RKQaPbPCWcps2DZZ/E854d9glKpigQ45qnHN3\nOOdOGX5l6ZzbSfoN10bwQ7J72Ht78NpH6wKWRUQKbf8+uOsjcOdKeMOb4IbNFRWYEDw0M80vz3ZN\npwjQPM7z8wKWRUQK6dlHoPUSePg2uHg1fOj7EDk+7FIVXNDQXOzfoo/g73l+cjYn8Nsuq0c9HQM6\nApZFRArBOfjVV+DWy73e8Q/eA5d90ttWtwIFHXK03h9ydC7eknERvGmUXc65FQHOM9TpY2YNeHur\nJ0Yf57+W6jBi0aLi2MJTpByNu+nZqYfDd2+E//khnHolXP0lmFXZN4XmXKZJPmneZHYS3tVhFV7g\n7ZzUh3u96G3OuawWO66pqXFbt26d+EARCWT0pmcAF8/YTsusrzDzQB/UfgbOu6Fsx16aWadzriab\nYzOt3J62N9wPyUkF5SjNQH0OziMiUzB807NDGOAvp9/Fn0+7i137j2XRyg5YMKZVrmJlatOc69+K\nf9nf+iKnzGw13tCjpP84NsFbRCRPUpueLeAFbjv0s3x0+p3cNXgRb+//rAJzlEyLEN8B3AHeknBm\ntgFvUY4O59yU9gIyszq8NtFe/xY9CtSQ/dYaIpJDC+ccxsX7fsCa6d/GcPzVa3/G3YMXsjBS/oPV\ng8qqIyhDgLYFGNCO//4o0DbOS9lu4iYiudTbw11HrWP+/gfZPPAmPn7wep52b2DmjENYtWxJ2KUr\nOpOZez6lAPUHs5dna7JIKRkcgAdbYNNnmD9tOtvO/jSrHz+L3XtfZWGq99zfDE1eN6WBVqkANbM5\nQHwqV6AiUkB/2A7f/XN45iE45Qp45xc4d85CNoddrhKQk9Gpzrm9eFefqQBdbmY34C3ssSYXnyEi\nOTBwAH75Rbj/Jjh0FvxxK7x5edkOJcqHnA/p9wP0Fv9LRIrFc496A9WffRjOuBre/jk48piwS1Vy\ngq5yFGgJOBEpAgf3w88+B7/4PMycC8u/4YWmTErGuedmNtvfyiKlcZxjLtdOlCJF6plOaLkEfrYO\nzqyDGx9SYE7RRFeaS4EOM+sDEkDEzE5wzj2VOiC1GLGZXe+cuzV/RRWRrB3oh5/8EzxwMxx5LLx/\nA5yqHWlyIWNo+oE4zV+gowZoArb5nT0JvJWJEnjrYZbfDkoipeipX3o9473dEPsgXPGPcPicsEtV\nNrJaGs45t83f27zdOVcFnAK04y0H105u5qGLyFTsfwm+/zfwtatg8CB84Lvw7i8qMHMsaO95B4xY\nbV095CLFoPs+uOejsPdpOO8jcPnfeUOKJOeCrqeZdtfJ8fZDF5E860/Cjz8J274F806BD/8IFp0f\ndqnK2pTHafq9658AVgGHTLlEIpKdx38A3/sYvPw8XPgxuOTjMOPwsEtV9oJudzHEH450E157pjqB\nRArl5T3Qfh1851o4Yh5cvwnin1JgFsikQtPMVgFPAnOAqHPuhlwWSkTG4Rz85k64+Y/gse/C29ZA\nw/2wUEvRFlLQGUGr8PbsaQdO8qdMiki+7emGH62BHRthwblw9T3eFrpScFmFppmtxBuj2Ya397nC\nUqQQ9r/kTYF84GaYfhhc8Vmvd7xCd4IsBhP+y/vLvT3knMtqi14RyQHn4JEN0PH3sO85OPv9EP8H\nOOrYsEtW8SYMTefc8kIURER8u7fBD1Z7a10uiMH7/huOy2qjRCmAoG2acwCn8ZgiebDvD3DfZ6Dr\nmzDraLj6Zu8Kc9qkB7lIHgQd3L7XzOaY2YlAr8JTJAcGDsCWW+Ena+HAy/CWG+GS1Zr+WKQms0fQ\nXiAVnpfjrc7+ZM5LJlIJeu6HHzbBHx6H6KVwVTPM12ZmxWzSXXB+eG4ys5PM7L1oszSR7PU95U1/\n/N29EDkB3vdtWPJ2bTtRAqY8bsE5txPYaWaLc1AekfL22iuw+Quw+d/ApsFlfwtv+QvN5ikhORvs\n5YeniIzHOXjsbtj4t/DiM94q6rWfgTnaIrfUaISsSL7972+9dssnfw5vOAve0wonvjXsUskkKTRF\n8uWVXrh/rdczfvgceMfnofpDME2LgZUyhaZIrg0OQNfXYdM/wqtJqLkOLv0EHFE16VPevW0X6zdu\nZ3eynwWRmaxatoRrztWtfRgUmiK59NQD8MNV3h7jJ1zoDSE69swpnfLubbtYc+ej9B8YAGBXsp81\ndz4KoOAMgaYaiORC8mm443r42pXwSh/UfQ0+9L0pBybA+o3bhwIzpf/AAOs3bp/yuSU4XWmKTMW+\nP8DP/wW2/idgcPFquPCvcro/z+5kf6DnJb8UmiKT0Z+EX/47/OrLcPBVOOf9cEkTRI7P+UctiMxk\n1zgBuSAyM+efJRNTaIoE8drL8OBXvMHpr+6FM98Lb/sEHJ2/lRNXLVsyok0TYOaMQ1i1TNMtw6DQ\nFMnGwf3Q+V/egsAvPw+nXgmXfhLe+Oa8f3Sqs0e958VBoSmSycBBeOQ7cP9N3p7iJ1wIK74Fi84r\naDGuOXehQrJIKDRFxjM4CL/7Ltz3T7Bnh7cY8Lu/6K1EpEU1KppCU2Q452BHh7cY8HOPwvzTYcV/\nw2nvUFgKoNAUed2Tm2HTZ+DpX8HcE+GPW+GsOk17lBEUmiK7t3lTHrs3wZHHenPEz/0/MP3QsEsm\nRUihKZXr+cfhJ5/1FgKeWQW1/wh/tBJmaPyjpKfQlMrT96TXG/7I7TBjFrxtDZz/Z3D47LBLJiVA\noSmV46Xn4GfrofPrXjvlW26Et34MZs0Lu2RSQhSaUv5efgF++UV4sBUGD0DsA3DxKpi9IOySSQkK\nLTTNrBnocM4lwiqDlLm+J+GX/wHbvuXND3/zcnjbx6EqGnbJpIQVPDTNLA7EgDqgo9CfLxXguUe9\nueG/udPbvOzsFXDBR2H+qXn5OC0QXFkKHpr+lWXCzGoL/dlSxpyDpzbDL/4VnkjAoUfC+R/x2i3z\neBuuBYIrj9o0pbQNDsL278MvvgC7tsKs+XDZ38HS62Dm3Lx/fKYFghWa5UmhKaXp4H54ZIN3G75n\nhzeD5x3/Auf8SUHHWWqB4Mqj0JTS8uqL3hJtv/oSvPQsHHsW1H0VTr8aDin8r7MWCK48Ck0pDfue\n9xb/fehW2L8XTroYrr4ZFl8W6kIaWiC48hR1aJpZA9CQerxo0aIQSyOh6O3xtpXY9t8w8Bqc/i5v\nD56F1WGXDNACwZWoqEPTOdcKtKYe19TUuBCLI4X07MNe585jd8O06XD2tXDBX+Z1W4nJ0gLBlaWo\nQ1MqjHOw86deWPb8BA6b7QXl+R+Bo44Nu3QiQDiD22PACiAOVJnZ7c65dYUuhxSRwQFvpaHNX/CW\naTvyDRD/FNR8GA6fE3bpREYIY3B7F9AFNBX6s6XIvLoXHv4OPNgCvd3e9MZ3/Ru8+X0w4/CcfIRm\n60iu6fZcCu/ZR2Drf3rjLA+8AgtroP7rXidPDldJ12wdyQeFphTGgVfhse96Yfn0gzB9preVxNLr\nYMG5eflIzdaRfFBoSn71PQlbvwbbvgmv7IGqxbBsLZxzLXc//grrv7Gd3cnv5+XWWbN1JB8UmpJ7\ngwPwxCbYcivs+LE3+HzJ22Hp9XDSJTBtWkFunTVbR/JBoSm58/Ie74py61ch+ZTXC37xKqj+IMw5\nbsShhbh11mwdyQeFpkyNc/DMVu+q8rd3wcB+OOFCb8jQae9Mu6NjIW6dNVtH8kGhKZPz2svwaLsX\nls89Aoce5W0jsfQ6OOb0Cd9eqFtnzdaRXFNoSjAv7IAt/wm//ra3cMYxb4J3/iucVQ+HHZX1aXTr\nLKVKoSkTGzgI23/gXVXu/ClMmwFnXO117Cw6f1KrDOnWWUqVQlPSSz7tXVF2/he8tBvmHO+tih77\nABx5zJRPr1tnKUUKTRlp/0vw2D3w8G3w5M+9506Owzs/D6dckdMZOyKlSKEp3rjKnvu9eeC/uxcO\n9rNv1iJum34t33j5PAafOYFVb1rCNQpMEYVmRfvfx7wrykfbvK0jDp8D51zLT2fGueH+afQfGPSO\n05xtkSEKzUqz73lvqNDDt3lDhaZN9267z26GU5bBjMP5xE330X9g5HAgzdkW8Sg0K8GBftj+Q+/2\n+4kEuAFvkYyr1sGZ74VZR484XHO2RdJTaJaoCdeJdA5+/yvvivK3d3tjKo9aAG/9S2+9ymNOS3tu\nzdkWSU8vI8/jAAAIbklEQVShWYIyLnZxwn54+HYvLJNPwYxZcMa74ez3wYkXZdX7rYHnIukpNEvQ\n6MUuZrOPdw4+yEn3/AO4xwGD6CVw6Se8+d+HHRno/Bp4LpKeQrME7U72cygHuGjaI7znkJ8Tn9bF\nYXaQHQML4YpPwVnLYc7UAk4Dz0XGp9AsJa+9At2b+MoRrbxlYAuzrZ897ii+PXA5dwxcRN/sM9h8\n4eVhl1KkrCk0i93+l+B/NsLv7oEdHXDgFS6dMYd7D57PvQdq+MXgWRxkOjNnHMLaK9N37ohIbig0\ni1F/H2z/kbenTvd93hqVs47xOnPOuJpDT7iQQx75X3Zs3M5Asp+FanMUKRiFZh5MatvYl1+Ax7/n\nzfve+VMYPAizj/P2/j7j3XD8eSN6vtXmKBIOhWaOBdr75sVnvbnev7sHntoMbhDmngRvuRFOvxoW\nxia17JqI5I9CM8cm3Pum76nXg/LpB70Djl4CF/21t0blG85UUIoUMYVmjo031fAke5arXnoIWj4N\nz/7ae/LYs+DSv/Vuvedr0LhIqVBo5pg3BfEVTrVnuGraQ1x5yEOcPu1p78VpNVD7GTj9XVAVDbeg\nIjIpFReak+qkycb+fbDzZ3zr2Ls4rP8+FtgLDDpji1vCPw9+kOorP8CyC2qm/jkiEqqKCs1AnTQT\ncQ5e+B/Y8WNv/OTvH4CB1zjp0CPZ/cbzaH7hVNr3ncWhkTeyatkSlqmnW6QsVFRoTthJMxH/apIn\nOmBHAvb+3nt+/ulwXqO3LuXx57Ng+qE0AU25r4KIhKyiQjPwOpHOeVvW7vixF5RP/RIGXvNWDoq+\nDS76GJxcC5Hj81ZmESkuFRWaWa0T+drLsPPnrwdlMnU1eRr8UQOcUguL3gLTDytQqUWkmFRUaI6/\nTuQ0Pn3BDHjgZq9t8qnN/tXkEd7V5Fv/ygvKyKLQyi0ixaOiQjPVbvnvP3qYE1/q5O0zf8OVh/2G\nWfc94x1w9Kne1eTJcTjhAl1NisgYFRWaANecMoNrvvchOHQ/TDsCjrsYTvHbJueeEHbxRKTIVVxo\ncuQxcPEqOK4aFl0AMw4Pu0QiUkIqLzQBLlkVdglEpERNC7sAIiKlRKEpIhKAQlNEJACFpohIAApN\nEZEAFJoiIgGENuTIzFYDPUAUSDjnusIqi4hItkIJTTNrA9amgtLMOoDaMMoiIhJEWLfn8VFXlj1m\nFg+pLCIiWSt4aPrh2DPq6SS60hSREhDGlWZknOf24LVtiogUtTBCsyqEzxQRyYkwOoJ6x3lu3ngH\nmlkD0DDsqVfN7LejDpsD7M3y8fDvjwZeyLLMmYz+vKkcm+718Z4Pu97pyjXZY3NV90qt9+jH+l0P\nVu/s14V0zhX0C4gD3aOeawaas3hv60TPZXo86vutOarPmDJN9th0rxdjvYu17pVa70LUvVLrPfqr\n4LfnzrkEY2/Ro0BHFm+/N4vnMj0e7/1TFeScEx2b7vVirHfQ8xaq7pVa79GP9bueJ+YnckGNM06z\n0zlXXeAybHXO1RTyM4uB6l15KrXu+ap3WOM0VwIrzKzOzJr9x4XWOvoJM4uFUI5CG1FvM2vwv1rC\nKlCBjK533P9qNrPxRnSUkzG/6wD+/71yNvpn3mJmEf/3fdI/81CuNIuRP360xTm3OOyyFEpqzKxz\nrsef1opzbl3Ixco7/4/jGudcvf/HotM5N26wlCv/3+CWQt/hhcnM+vxvm6by8y7bBTv8K4gxs4zM\nbLV/hbt6+JWl39Y6etB9yQlY7yhQ53/fA5TsH4wg9XbOdTnn6v1DokCikGXNtaC/674qxh/JUjIm\nUe+Vzrm5U/4DmY/epTC/8HrnVwPdeNM1h7/WBsSGPe4Y9XpHvstXjPX2n2sB6sKuR4F/3g1AQ9h1\nKHTdgRjeJJOS/H2fQr1Xp947lc8vuytN51zCebeY4101lu2c96nU28yi/jna81zMnJtKvZ13xVFd\nqr8DU6h7lXMumf8S5sdk6+2cW+e8O8pkqjlqMsouNNOp1DnvWda70TnXWLhS5V+meptZbNhtWydQ\nUXUHtha+VPk3Qb3r/MkyMMWmqIoJTSp3znvGeptZHbDW/74kr7jSyFTvOK//3COUQVv2KJnqHgXi\n/s89WkE/8x5gg/9cFO+P5aRUUmhmnPM+7JdodZkNQUlbb/8/zC1Ap9+zWE5DrtLW27+1q/J/5ovx\n/2iUkUx1bx/WDFNu60BkqncXsDz1M3dT6AwKbeX2EGSc8+7/IpVcm14W0tbbb9+ZW9jiFMxEP+/U\nf5qK+pmnlOnve0F+5pV0pZlk/Mv3crs1G031Hqnc6w2VW/eC1LtiQtNNbc57yVK9Ryj7ekPl1r1Q\n9a6Y0PQlRg12jfr/0OVO9fZUSr2hcuue93qXXZum/w+2Aq+HtMrMbnevTw1cCazxxyUuJZw573mh\neldWvaFy6x52vTX3XEQkgEq7PRcRmRKFpohIAApNEZEAFJoiIgEoNEVEAlBoiogEoNAUEQlAoSll\nw8yi/uZZzt/xdPTrbanXKmQTPckDDW6XsuLPBGnC28pi7ugVys2s2TnXFErhpCzoSlPKTQwvNJN4\nwTnEXyd1SxiFkvKh0JRyk9r/ppWx21jUAF1j3yKSPYWmlKsWvJX4R694U+5rSkqeKTSlbPi3370A\nfjh2AWtCLZSUHYWmlJM4MHztxLVAXUhlkTJVdutpSkUbsZ+3c67dzPC3bu1hZKCKTIquNKXcpTqE\n1J4pOaHQlLIwvD1zlBa8YUiLC1siKVcKTSkXyxlnOJG/33UP0F3wEklZUmhKSTOziD9lsgVo8WcE\njdaM2jMlRzSNUkQkAF1piogEoNAUEQlAoSkiEoBCU0QkAIWmiEgACk0RkQAUmiIiASg0RUQCUGiK\niATw/wFHGC2B/93BugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111a9d278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(Ns,numpy.nanmean(T,axis=1)*10.**6,'o')\n",
    "p= numpy.polyfit(numpy.log(Ns),numpy.log(numpy.nanmean(T,axis=1)*10.**6/numpy.log(Ns)),deg=1)\n",
    "plot(Ns,numpy.exp(numpy.polyval(p,numpy.log(Ns)))*numpy.log(Ns))\n",
    "pyplot.text(10,4.5,r'$\\Delta t  \\approx %.2f\\,\\mu\\mathrm{s} \\times N^{%.2f}\\log N$' % (numpy.exp(p[1]),p[0]),size=16.)\n",
    "gca().set_xscale('log')\n",
    "xlim(5,150000)\n",
    "ylim(0.,5.)\n",
    "xlabel(r'$N$')\n",
    "ylabel(r'$\\Delta t / \\mathrm{collision}\\,(\\mu\\mathrm{s})$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The behavior seems to be more like $(\\log N)^2$ or $N^{1/4}$, probably because the implementation of the binary search tree to determine the next collision is not optimal.\n",
    "\n",
    "## Number of collisions per dynamical time\n",
    "\n",
    "In a dynamical time, every particle will cross every other particle, so the total number of collisions per dynamical time should scale as $N^2$, which is indeed what we observe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x106f3eba8>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAFKCAYAAABch5ftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW5P/DPQ0ggKDAkBoUAQkBwY0kIKGqvVoO0eqkU\nAqho22tlE62tLSD2+uutXTBE673WqgTbXhVEFjEWa00T7W21dcuCgEpAwhqUfcKSBLI8vz/mDE6G\nSTIzOTPnzJnP+/XKiznnTGaew0ye+c5zvouoKoiIKHI6WR0AEZHTMdESEUUYEy0RUYQx0RIRRRgT\nLRFRhDHREhFFGBMtEVGEMdESEUUYEy0RUYQx0RIRRVhnqwOIFBGZCGBi9+7dZw4dOtTqcIjIYcrK\nyg6palow9xWnz3WQnZ2tpaWlVodBRA4jImWqmh3MfVk6ICKKMCZaIqIIY6IlIoowJloioghjoiUi\nijDHJloRmSgiBTU1NVaHQkRxzrGJVlXXq+qsnj17Wh0KEcU5xyZaIiK7iJmRYSKSAWA7gCpjl1tV\nR1sYElHcK6yoRn5RJfa569DXlYz5E4ZhUma61WHZTswkWgAZqirAmaSbYXE8RHGtsKIai9ZtQl1D\nEwCg2l2HRes2AQCTrR9LSwcikiciOQH2LxCRXOPfLABQ1RKfu+T4bRNRlOUXVZ5Jsl51DU3IL6q0\nKCL7sqRFayTXLAC5AIr9jq0BsFhVy43tYgDjfY5nAODkBUQW2+euC2l/PLOkRauqJaq6BF/VW33l\neJOsocqv1Tvb7zgRRdmRk6fROUECHktOSkBjU3OUI7I3W/U6MBKqf/J1w6dFC+CsUgMRRc+B4/W4\nreB9NDcrkhJappDOnQS1p5swZ3k56v3KCvHMVokWgCvAvsMwLnyJSKDjRBQlX9TU4dal72P3kVq8\n8P0rsCR3BNJdyRAA6a5kPDZ1JH7+rcvw1pb9mPHcB3DXnrY6ZFuwW6+DlCDuE6jcQEQRtvtwLW5/\n7n3U1Dbgxe+PRfZAz59roB4Gad274Icvb0Dus+/hhbvGoq8rOdrh2ordWrRHAuxL9d5QVbeqTo1i\nPEQE4PMDJzBt6Xs4caoRK2ZecSbJtuam4X3w/F1jsb+mHpOf/hcqvzwepUjtyW6J1o3A5YOgW7Ei\nMktESr0/Bw8eNC86ojj02RfHcGvBe2hsbsbLs67EiH7BVfDGDU7F6jnj0KyKqc/+Cx9UHY5wpPZl\nq0Rr9I31/6jMgF8XsHYeo0BVs70/aWlBLelDRAFs3OvGrQXvo3OnTlg1exwuvqBHSL9/SZ8eWHfP\nVTivexfc+YcP8ebmLyMUqb3ZKtEaSryDFAwZ4QxO4OxdRB1TuvMIZiz7AN27dsaaOeMwOO3csB6n\nX69ueGXOVbisbw/cs6IMy9/fZXKk9mdJohWRLBHJg6erVp6ILPA5PBPAdGNkWJ6xHTLO3kUUvn9+\nfgh3/v5DpHXvgjVzxqF/SrcOPV6vc5Kw4u4rcN2w3vjPws34zV8r4fSFYX05dhVc73LjQ4YMmblt\n2zarwyGKGW9v2Y85y8sxKPUcLL/7CqR172LaYzc2NeOhVzdhdele3DqmP3456XJ0TrDjF+v2cRVc\nsEVLFI6/bPoCs18sw7Dzu+PlWVeammQBoHNCJ+RNGYH7rh+Clz/agznLy1B32vkDGxybaFmjJQrN\nqxV7Me+lcozo58KKmVeg1zlJEXkeEcGPbxyGX9xyGd7acgAznnsfR086e2CDYxMtW7REwXvpg914\nYPXHuDIjFS/cNRY9uiZG/DnvHDcQT9+ehc37jmHq0vdQ7eDJaBybaIkoOH94dwceenUTrhuahj98\nbwzO6RK9AaPfHN4HL9w1FvuP1WPy0//Eli+PRe25fRVWVOPqR9/GoAf/jKsffRuFFdWmPr5jEy1L\nB0Tt+93fPscjr3+Kb1x2AZbemY2uiQlRj+HKjFSsmTMOADD12feiPrDBO4F5tbsOiq8mMDcz2To2\n0bJ0QNQ6VcXjf61EflElJo3qi6duz0RSZ+vSwcUX9MArc69Cb2Ngw182fRHR52tqVnz2xTGs/HA3\nHnp1U8QnMLfbpDJEFGGqil/++TP8/t0duHVMf/zq28OR0Cnw3LLR1K9XN6ydcxW+//xHuOelcjxy\ny+W488oLTXnsL2rqsGG3Gxv2urFhtxubqmtQ205vBzMnMHdsovXpR2t1KES20dysePi1zVjxwW58\n76qB+NnESyFifZL18gxsuBL3rSzHw4WbceBYPR4YPzSkGE+easTGvTXYsMeNDXuOYsMeN/YfOwUA\nSErohEv69sC07P4Y1d+FUf1dmPHc+6h215/1OGbOOObYAQte2dnZWlrKlW+IGpuaseCVjVhXXo17\nrhuM+ROG2SrJ+mpsasZPX92MVaV7cOWgFOw+Wosv3PVnrbTb1KzYuv+4J6nudmPDHje2HTiOZiOt\nDUzthlH9XRhpJNVL+/ZAl84t69D+i0wCQHJiAhZPHt7mIpOhDFhwbIuWKN75LgXep2dX9O7eBRv2\n1uAnNw7FvddfZHV4beqc0AmPThkOd+0pFH164Mz+ancdFqz9GH/6uBonTzW1KAG4uiViZD8XvnH5\nBRg1wIVR/VxB9QX2JtNILpvOREvkQP6ttH019dhXU49Jo/raPsl6iQg27zt7HtvTTYq3txzEyP6u\nFiWAC1O7hd1Cn5SZHtEl0h2baFmjpXgWaClwAPho51ELoglfaxekBMBr866ObjAdwO5dRA7klKXA\nW7sgFWtL4zg20RLFM6ckqPkThiHZbxBFcmIC5k8YZlFE4WGiJXKgyVln1xtjMUFNykzH4snDW6y0\n215vADtybI2WKF7tP1aPlR/uQdq5Seic0Alf1pzdLSqWRPpCVTQ4NtHyYhjFo9ONzbhnRTlqTzei\ncN7VGHp+d6tDIji4dMCLYRSPfv3GZyjbdRR5U0YwydqIYxMtUbx5bUM1/vdfO/H9awZh4si+VodD\nPphoiRxgy5fH8OArmzB2YAoe/ObFVodDfphoiWLcsfoGzHmxDN27dsZTMzKRGKOLHTqZYy+GEcWD\n5mbFj1d/jL1H67By1pXo3b2r1SFRAPzoI4phz/x9O4o/3Y+f3nwJxgxMsTocaoVjEy2XsiGne2fb\nQTz+10p8a2RffO+qgVaHQ21wbKJl9y5ysr1Ha/GDlRW4qHd3PDpluG3nlSUPxyZaIqeqb2jCPSvK\n0dikePbO0eiWxEstdsdXiCjG/Hz9J9i4twYFd47GoPPOsTocCkJMJVoRyQCQBcANoEpVqywOiSiq\nVn20Gys/3IN5Xx+MGy+7wOpwKEixVjrIU9W1AFwAcqwOhiiaNu2twcOvfYJrhpyHB8bH1ixc8c7S\nFq2I5AEoVtUSv/0LAFQByABQoqrlIpJr7IORbInixtGTpzFneRnSzu2CJ2/LtMXy4BQ8SxKtiOTA\nUwLIBVDsd2wNgMWqWm5sFwMYD0/SdYlIFjyt2QJVdUc1cCILNDUr7l+1AQePn8KaOeOQEsSCg2Qv\nlpQOVLVEVZfAaKH6yfEmWUOVkZhTAbiNY+UAFkUhVCLL/XfJVvxj60H8/JbLMLK/y+pwKAy2uhhm\nJFT/5OuGp0W73WffEXhaxESOVvLpfvz27c8xLbsfbh3T3+pwKEx2uxgW6OP6MIxaLYDBxr4MeFq1\nRI6189BJ/Gj1Blye3gOP3HI5ByXEMLsl2lYHaxtdubYbF8XGqOrC6IVFFF11p5swZ3kZEjoJnpkx\nGl39Fiik2GKr0gE8JQF/qd4bqlpg3Gy114GIzAIwy7s9YMAA04IjigZVxUOvbkLl/uP43/8Yi/4p\n3awOiTrIbi1aNwKXD4IemKCqBaqa7f1JS0szLzqiKHjx/V14taIaD+QMxbVD+f51AlslWqM/rX/5\nIAN+XcCCwdm7KBaV7TqCR9Z/ihsu7o15X+fCok5hq0RrKDH6ynpl+A9oCAZn76JYc/D4Kdyzohzp\nvZLxm+mj0ImDEhzDqgELWQCmwzPwIEVEVhn9agFgJoBFxrwGY4ztcJ6Dy41TzGhsasa9L5Wjpq4B\nf/zeWPRMTrQ6JDKRqKrVMURUdna2lpaWWh0GUZt+9edPseydHXhi+kh8O7Of1eFQEESkTFWzg7mv\n3XodmIYtWrKzwopq5BdVYp+7Dq5uiTha24DvjLuQSdah7FijNQVrtGRXhRXVWLRuE6rddVAAR2sb\nIAKMSOd71akcm2iJ7Cq/qBJ1DU0t9qkCT5RssygiijTHJlp27yK72ueuC2k/xT7HJlqWDsiu+rqS\nQ9pPsc+xiZbIriZnpZ+1LzkxAfMncNUEpzI90YpID7MfMxwsHZAdHThWj5c/2oPzzk1Cn55dIQDS\nXclYPHk4JmWenYDJGTrcvctIrDn4aujseHgGI1hKVdcDWJ+dnR3WgAciszU0NePelypwor4RhfOu\nxrALulsdEkWJGf1ol6DlpNy9THhMIsdZ8uYWfLjzCP57+igm2ThjRqJdo6pveTdEhAsnEvl5Y9MX\nWPbODnxn3IUsEcQhM2q0KiKjRKSHUUaYYsJjdhhrtGQX2w+ewPw1H2NUfxd+evMlVodDFujwXAci\ncgTARwC8Uw0NUtWLOhqYWTjXAVnp5KlGTPrdP3H45Gm8ft817MLlINGe62CqX+kg04THJIp5qopF\n6zZh+8ETeOGuK5hk45gZpYMWiVVVK0x4TKKY9/y/duJPH+/Dj28chmsuOs/qcMhCZiTasSIyWURG\nmfBYRI5Qtusofvnnz5BzSW/MvXZw+79Ajtbh0oGqTgM8/WlFZDGAj1R1XYcj6yBOk0hWOXTiFOat\nKEdfVzIen8qVEsiEFq2IXC8ikwG8Dc/CijtEZIrVLVzOdUBWaGxqxn0vVeBo7Wk8c0cWenbjSglk\nzsWwtQB+DeAGVfX2paoQketNeGyimPJ48Va8V3UY+bkjcFlffsiTh+m9DoAzPQ9Gw9PKJYoLf/3k\nSzzzf9tx29j+mJrd3+pwyEbMqNH6J9keRs8D9j6guLHz0En8ePXHGJ7eEz+beJnV4ZDNhJVo2ykL\nzIYNJpUhipa6002Ys7wMCQmCp2dkoWtigtUhkc2E26J9EEAZPKPBMgBUGftdxg9RXFBV/LRwEyr3\nH8cfvzcG/VO6WR0S2VC4iXa2qu4AABG5wW9k2A2mRNZB7N5F0bDig91YV16N+2+4CNcN6211OGRT\nYXXv8iZZg/+lVVtcamX3Loq0DXvceGT9p7h2aBruv8E203uQDZnR6yDVGKiwHcBgAIdNeEwiWzty\n8jTuWV6GtO5d8N/TOSiB2tbhAQuqugxAKYAhAEpV9bEOR0VkY03NivtfrsChE55BCb3OSbI6JLK5\nkBKtiAwMtF9VX1HVB1X1FTOCauP5l4qIS0RmiQgvupEl/qdkK97Zdgj/9a3LMKIf34bUvjYTrc9k\n3l6zA9znBhG52/TIApsGYAcAqKo7Ss9JdMbbW/bjybc/R+7ofrhtLAclUHDaq9GOAVAsIkcBlABw\niciFqrrLewdvjwMRuVtVnwvlyUUkD0Cxqpb47V8AT5exDAAlqlpuHJqpqlwqhyyx50gtfrTqY1zS\npwd+ccvlEGFdloLTZotWVd9S1U7wrHJbAs/FrgoRaRKRIhH5iXcZG4TQf1ZEcoxkmhvg2Bp4kuta\nVV0CIM/ncIbP7xJFTX2DZ1BCsyqevSMLyUkclEDBC6pGq6oVxkWvtaqaAuAieCaTGWL8u6Ot3w/w\neCVGEq0KcDjHpwULAFUikmP83hKj9etmsqVo+tlrn+CTfcfwxLRRuDD1HKvDoRgTaveuYgBQ1Sp4\nkuQyM4MxEqp/8nUDGG9c/EpR1QLjPlPNfG6i1qz6aDdWle7BvK8PRs6l51sdDsWgVhOtiFyvqi1m\n3/KfQCYCApUfDsNTK14FT/kC8NRuyyIcC8Wxwopq5BdVYp+7Dgrgot7n4IHxw6wOi2JUW6WDXiKy\nWkSeMSb2joaU1g4Y5YRpIpILYLDRsiUyXWFFNRat24RqI8kCwJ6jdVj/8T5L46LY1WqL1ugT+woA\nGCsmrIandVkcwaVqjgTYl+oTkze5ttrzQERmAZjl3R4wYIBpwVF8yC+qRF1DU4t99Q3NyC+qxKTM\ndIuiolgWVI22jaS7xr+80EFuBC4fBLpoFpCRjM+0drOzs7WNuxOdZZ+7LqT9RO0JeQiuMQpsmqrO\nRcvyQoeXrjF6FPiXDzJgXIQLhYhMFJGCmpqa9u9M5COllSG1fV3JUY6EnKJDcx14ky4889OalXRL\nRCTLZzvDf0BDkLFx9i4K2d6jtag93Qj/oQjJiQmYP4EXwyg8HZ5UBgBUtcYv6Q42ku7iQPcXkSxj\nVFgOgDy/PrEzAUwXkVzjPjPDiYktWgpVfUMT7llRjs6dOuGhmy5BuisZAiDdlYzFk4ezPkthE1Vn\nlzCzs7O1tLTU6jAoBixatwkrP9yNpXeOxoTLLrA6HLI5ESlT1exg7hvSgAWf2buOqOoxERkF4FYA\nn4c6zwGRnawu3YOVH+7GnGsHM8mS6UItHcyB5+t+irGk+BoALwMoi+IMXkFh6YCCtbm6Bg8Xbsa4\njFT85MahVodDDhRqoi1W1edUdSc8tdgCVd1gLC8e0nwHkcaLYRSMmtoGzF1Rhl7dkvDb2zPROcGU\nyxZELYT6rvJdpiYHwFKfbVsVe9mipfY0Nyt+uKoCX9bU4+k7snDeuV2sDokcKtREO9iYDHw+gDJV\nPQZ45kUwP7SOYYuW2vPU3z7H3yoP4uF/vxRZA3pZHQ45WEiJ1hgh5l1lYSoAiMij8Kx8kGFuaESR\n8/etB/FEyVZMGtUXd155odXhkMOFvAququb7bT9oXjhEkbf3aC3uf7kCQ3t3x68nD+dKCRRxplX+\nWxucYBXWaCkQ76CEpibFs3eORrekkNsaRCELdRXcG0SkVEQO+/0cAWCrFQ9Yo6VAfr7+U2zcW4PH\npo3EoPO4UgJFR6gf5+MB3KCqZzUTjVotkW1xUAJZJZx+tAG/i9utVsvSAfnyDkq4ajAHJVD0hZpo\n1Vjx9ixRXIUhKCwdkJfvoIQnb+OgBIq+UEsHNwJYKCJutFwNQQDcACBSKy8QhcV3UMKq2eM4KIEs\nEWqizQHQWu8C9qMl2/EOSnjklss4KIEsE2qinWnMa3AWEQl6uRmiaOCgBLKLUEeGnUmyxlDcgYGO\n2QEvhsU3DkogOwn5qoCIzBSRz+GZravc6EdrqykSAV4Mi2cclEB2E+qAhSkAoKpDVDVVVVPgqc2K\nHZMtxScOSiC7CWcV3GV+2zXGPn43I8t5ByXMvY6DEsg+Qu5H28axw20cI4o430EJPx7PQQlkH6Em\n2tRAAxaMfYPNCYkodByUQHYW0lUCVV1mLCOuALzduUZ7DukE06MjCgIHJZDdhTMf7TQRuQFAlrEr\nT1XfMjcsouBxUALZXVj9XozE2iK5isj1qvq2KVGZQEQmApg4ZMgQq0OhCOKgBIoFotr69S2jy1aJ\nseqtd3JvV6C7wjN94kWRCLIjsrOztbS01OowKAL2Hq3Fv//2XZzfvStenXcV+8tSVIlImapmB3Pf\n9t6Z2QB8s9R4eOY6cAe4L+c6oKgorKjGkje3YF9NPQTAvOuGMMmSrbX57lTVOX672prr4Eig/ZEg\nInmqujBaz0f2UVhRjUXrNqGuoQmAp7/hb4q3Iq17F0zKTLc2OKJWhNoH5qiIDPR28RKRUSKyWETu\njtZcByKSBc8sYhSH8osqzyRZr7qGJuQXVVoUEVH7Qk20c+BJcikikglgDYBVAMrCGYIrInkiclbS\nFJEFIpJr/JvldzgFLefCpTiyz10X0n4iOwi1sFXs7colInkAClR1g7GdEuyDGMk1C0AugGK/Y2sA\nLFbVcmO7GJ7asLc1yytbcaqmtgGdOgmams++gNvXlWxBRETBCbVF6zvMNgfAUp/ttobntqCqJaq6\nBF8NevCV402yhiqfVm+Kqga6EEcO19yseGD1BjQ3K7p0bvm2TU5MwPwJwyyKjKh9obZoBxsTfM8G\nUKaqxwBPH1ozgjESqn/ydQMYb1xsY2s2Tj3z9+14a8sB/NfES+HqloT8okrsc9ehrysZ8ycM44Uw\nsrVQh+C+IiLzjc2pwJllxl3wJMGODlgI1Ef3MIAx8HQfyzAmcM4QkRxVLeng81EMeHfbITz+10p8\na2RffPeqgRARJlaKKeEMwc33234QMK1V22qdV1XXGs+T29b9yFn2uevwg5crMDjtXCzmSgkUo9pM\ntCGMDAM8NduOjgwL1Jsg1XfDSLhrO/g8FANONXpWSjjd2Ixn7xyNc7pwUALFJjNHhpkxTaIbgRN5\n0As/isgsALO82wMGDDAhLLLCr/78GTbsceOZGVkYnHau1eEQhc1WI8NUtSRAN7EMtOzd0N5jFAAo\n8E4qk5SUNLqjcVH0FVZU44X3dmHm1wbhm8P7WB0OUYeEvQpuAIM6GItXid8ghYxwLnpxccbYteXL\nY3hw3UaMHZSChd+42OpwiDqs1RatiEyGMVAgCALgBgDrgrqzJ5FOx1ejzFYZ/WoBYCaARSKSAU9v\ng5lBxuD/HJwmMQYdq2/A3OXl6N41EU9xpQRyiFanSRSRZwFsB1Ae8A5nW6iqN5oVmFk4TWLsUFXM\nWV6Gks8OYOXMKzF2EDuXkH2ZNU3i0lAmionm7F3kTAX/qELRJ/vxnzdfwiRLjtLq97IwZuMyq0Zr\nChGZKCIFNTU1VodCQXhv+2HkvbkFNw2/AN+/xlZvJaIOa6t0EHKNlissUDj2H6vHzU++ix7JnfGn\ne6/BuewvSzHArNLBjQitRssVFihkDU3NmLeiHLWnG7Fy5hVMsuRIjq3RstdBbFj8xhaU7jqKJ2/L\nxEXnd7c6HKKICLtGKyI9RGRgsPePNvajtb/XN+7DH/65A9+7aiC+NbKv1eEQRUzInRRFZKaIfA5g\nB4ByETkczuoKFN8+P3AcC9ZuxOgLe+Ghmy6xOhyiiAop0YrIFABQ1SGqmqqqKfDUZsVuyZa9Duzr\nxKlGzH6xDN2SEvC727OQ1JmDEsjZQn6Hq+oyv+0aY5+t5q9j6cCeVBULX9mIHYdO4snbMnFBz65W\nh0QUcaEm2raWqzncxjEiAMAf/rkTf974BeZPuBhXDT7P6nCIoiLURJvqXWrcl7HPjGkSycE+2nkE\ni9/4DDdeej7mXMvegBQ/Ql3KZpmIrBYRxVdzxI72HNIJpkfXAezeZS8Hjtdj3opy9OuVjMemjeRK\nCRRXwqnRTgNQAM9qCEcA5NktyQKs0dpJY1Mz7nupAsfqG/DMHaPRo2ui1SERRVW4w3C2q+pb3g0R\nGehd7obIX35RJT7YcQRPTB+JS/qcVXkicrxQu3eNMrp4FfsOVvAcMmfJcXKWNzd/gaX/qMIdVw7A\ntzP7WR0OkSVCLR1kq+orqnqRbwtWVXeg9UUbKU5VHTyBn6zZiJH9XXj43y+1Ohwiy4SaaNuaz8BW\nE4hywIK1ak83Yu7yciQmCJ6ekYUunROsDonIMqHWaAeLyChV3eC70ygb2OryvqquB7A+Ozs7rKVw\nKDyFFdXIL9qCanc9AGDOtRlIdyVbHBWRtULt3pVvdO/KhGf6RBc8Q3DLVXV6JAKk2FFYUY1F6zah\nrqHpzL7n/7ULF1/QA5My0y2MjMhaIfc6UNVpIjIIQBY85YISo0ZLcS6/qLJFkgWAuoYm5BdVMtFS\nXAure5eRWJlcqYV97rqQ9hPFC06bRKZobtZWZ+HqyxotxTkmWjLF7/72OU41NiMxoeXQ2uTEBMyf\nMMyiqIjswbGJlt27oufdbYfwm5Kt+HZmOpZMGYF0VzIEQLorGYsnD2d9luJeq6vgOgVXwY2sL2vq\ncfOT7yD13CQUzrsa3ZK4uCLFh1BWwW2zRRtoSkQir4amZtz7UjnqGprw9IzRTLJErWivdPCc7xwG\nIjIqwvG0SURyjJ88EeGQX4stedOzgu2jU0ZgSO9zrQ6HyLbaS7SrVPVtn+2c1u4oIj8xJ6RWHz8L\nwGxVLYFnoMS0SD4fte3NzV9g2Ts78N1xF3IFW6J2tPddL0VEVsEzybcAuEFEAq2kIACmAngslCcX\nkTwAxUby9N2/wHjODHgGRJSrarnxHDD254XyXGSeHYdOYr4xWcxDN3MFW6L2tJlojRUVSuBJbDD+\nXdvK3YNem0REcuAZWZYLoNjv2BoAi43EChEpBjDe5/gsAGtUtQoUdfUNTZi7vAwJCYLf3Z7JyWKI\ngtDu1QvfUWAiAt8Jv32JSFsze/k/ZgmAEhEZH+BwjqpO9dmuEpEcb6tXVQtEZKnvPoqe//faZmz5\n8jj++B9j0K9XN6vDIYoJIfWj9VtVoYfv5N+qWtHRYIyWrn9L1Q1gvIhkGXVaACgDMLujz0ehWf3R\nHqwu3Yv7rh+Crw/rbXU4RDEj5AELIjJTRD6Hp5VbLiKHReRuk+IJ1JPgMDxliRx8VZ5w4eyETBH0\n6b5jePi1zbhqcCp+mDPU6nCIYkqoS9lMAQBVHaKqqaqaAk/yE5OSbauTh6vqEnguzuXCs7T5YhOe\nj4JwrL4B96wog6tbIp68LRMJnbiCLVEowpkmcZnfdg2AZSJixgTbgeq8qT7PVWDcbO2CHJlMVbFg\nzUbsOVqHl2ddifPO7WJ1SEQxJ9RE29Z43cMdCcTgRuDyQdBlAqNXwizv9oABA0wIK379/t0dePOT\nL/HTmy7BmIG2Wq2IKGaEWqNNDTQs19gXqH9tSIxeBP5/zRnw6wLWzmMUqGq29yctLa2jYcWt0p1H\n8OhftuDGS8/H3V8bZHU4RDEr1KVslhlL2Si+amWO9hzSCSbFVCIiWd5+tAAywunGJSITAUwcMsRW\nS5nFjEMnTuHelyqQ3isZ+VNHQoR1WaJwhdzrQFWnASiAp556BEBeqEnW6KqVB09PgjxjJJjXTADT\nRSTXuE9YtV9VXa+qs3r27BnOr8e1pmbFD1/egCO1p/H0jCz0TE60OiSimObYaRJ9WrQzt23bZnU4\nMeU3xVvx5FvbkDdlOKaPYY2bKBDTpkmMZWzRhuf/Kg/gt29vQ+7ofpiW3d/qcIgcwbGJlisshK7a\nXYcfrdqMp/IXAAAPuklEQVSAYed3xy9uuZx1WSKTODbRskUbmtONzZi3ohwNTYqnZ2QhOYmTxRCZ\nhVPiEwDg1298hg173Hh6RhYy0jiJN5GZHNuiZekgeK9v3If//ddO3HX1INw0vI/V4RA5TliJVkSe\nCXTbTlg6CM7nB05g4dqNGH1hLyy66WKrwyFypPYWZ7y+lQUaU1u5TTGk9nQj7llRhi6JCXjq9kwk\nJjj2Cw6Rpdqr0QqAh4yFEBWeeWBL0HLOA1t2xOXIsLapKn766mZsO3ACL9w1Fn16JlsdEpFjtbeU\nzVsAfCf7zoRn3a7RRsmgCoEngbGcqq4HsD47O9uMWcUcZ+WHe/BqRTV+lDMUX7uI80EQRVKocx1U\nAKgQkTGqOldEBsEzXNZbp90OoNxv5VyyicKKauQXVWKfuw4KYNgF3XHf9WzxE0Vam4lWRB4FcAie\nlWg3+BxSwLOemIhUqepc4/6DAGSJyLMAjqrqogjFTSEqrKjGonWbUNfQdGbfrsMn8aeP92FSZrqF\nkRE5X3ulgwd9kucceBJsDQCXiPRQ1WPwqdH6LOT4SgRjDgprtC3lF1W2SLIAUN/QjPyiSiZaoggL\nZRXcVwBARHoCWANgkYj0gicJT4an1XssksGGgjXalva560LaT0TmCWcpmxoRcXvLAiKyGp5W7lm9\nE1R1p5nBUvhSz03CoROnz9rf18XeBkSRFu4Q3Bbdu1rpnZAFYGf4oZFZ9h+rR31DEwQtX7jkxATM\nnzDMqrCI4ka4PdSPtnIbgKd3gqquC/OxyUQNTZ7JYpoVWPiNYUh3JUMApLuSsXjycNZniaIg3Bbt\nwlZu2wYvhnnk/WULSncdxf/cOgq3jErHnOvi+/+DyAphtWiNJcbPum0nnOsAeGPTF3ju3R347rgL\nccsotlyJrMLB7Q61/eAJLFi7EaP6u/DTmy+1OhyiuMZE60C1pxsxd3kZEhMET8/IQlJnvsxEVuLE\n3w7jP1kMu28RWY9NHYdZ/sFuThZDZDNMtA6yYY8bv1j/Ka4bloZ7v87eBUR24dhEG29L2Rw9eRrz\nVpQjrXsXPDFtFDp14gq2RHbh2EQbT927mpoV96/agIPHT+GZO7LQ65wkq0MiIh+OTbTx5Ldvb8M/\nth7Ez751KUb0s+U87ERxLaZ6HYjILOPmaFWdbWkwNvH3rQfxP29tw+SsdNw+doDV4RBRADHTohWR\nHHhmBCsAsF1EFlgdk9X2Hq3F/S9XYNj53fGrScMhwroskR1ZmmhFJM9IoP77F4hIrvFvlrE7A0Cu\ncbsKwOBoxWlHpxqbMG9FOZqaFM/cMRrJSQlWh0RErbCkdGAk1yx4Emex37E1ABararmxXQxgvNGS\n9Rrv/3vx5pevf4aP99bg2TuyMOi8c6wOh4jaYEmLVlVLVHUJPC1TfzneJGuo8m31ikiG8RhrIxym\nbRVWVOPF93dh1r9l4BuX97E6HCJqh61qtEZC9U++bnhasF6z4/lCWOWXx7Fo3SaMHZiCBZy0mygm\n2CrRAgjUN+kwPPVZiEgugMXG7bNqu053vL4Bc5eX4ZwunfHU7ZnonGC3l4+IArHbX2pKaweMxLoM\nQJmIHIWnxhs3VBULX9mIXUdq8dTtmejdo6vVIRFRkOzWj/ZIgH2pgKeuC6BXdMOxj9+/uwNvbPoS\ni755Ma7MSLU6HCIKgd0SrRuByweBLpoFZAxq8A5swIABsd+J/6OdR/DoX7bgxkvPx6x/y7A6HCIK\nka1KB0ar1b98kIEQunKpaoGqZnt/0tJie6rAg8dPYd6KcvTrlYzHpo3koASiGGSrRGso8RmkAAAZ\nRgIOiRNm72psasYPVlbgWH0DnrljNHp0TbQ6JCIKgyWJVkSyRCQPQA6APL/htDMBTDdGhuUZ2yFz\nwuxdjxdvxXtVh/HLScNxSZ8eVodDRGGypEZrDEgoR4ClylXV7bM/7EEJsb7cePGn+/HM/23HbWMH\nIHd0P6vDIaIOsGPpwBSx3KLddfgkHli9AZen98DPJnIFW6JY59hEG6s12vqGJsxdXo5OInhmxmh0\nTeRkMUSxzm7du0yjqusBrM/Ozg6rxhtthRXVyC+qRLW7DgAw82uD0D+lm8VREZEZHNuijSWFFdVY\ntG7TmSQLAMvf343CimoLoyIiszg20cZS6SC/qBJ1DU0t9tU1NCG/qNKiiIjITI5NtLF0MWyfT0s2\nmP1EFFscm2hjRd3pJnROCDzaq68rOcrREFEkODbRxkLpQFUxf+3HaGhSJPlNeZicmID5nG+WyBEc\nm2hjoXTw1Nuf4/WNX2DhNy7GktwRSHclQwCku5KxePJwTMpMtzpEIjKBY7t32d2bm7/A48Vb8e3M\ndMy5NgMiwsRK5FCObdHa2Sf7avCjVR9jVH8XFk/mMuFETufYRGvXGu3B46cw8/lSuLolouA7HPlF\nFA8cm2jtWKM91diE2S+W4kjtaSz7TjZ6d+dyNETxgDXaKFFVLFq3CeW73fjd7Vm4PN0+HwBEFFmO\nbdHazbJ3qrCuvBo/zLkIN4/oY3U4RBRFTLRR8PaW/Vj8ly24eXgf/OD6i6wOh4iizLGJ1i4Xw7bu\nP44frNyAy/r2wGNTR6JTJ/YwIIo3jk20drgYduTkadz9fCmSkxKw7DvZSE5iDwOieMSLYRFyurEZ\nc5eX4ctj9Vg160r06cl5C4jilWNbtFZSVfzsT5/ggx1HsGTKCGQO6GV1SERkISbaCHjhvV1Y+eFu\n3HPdYA6rJSImWrO9s+0gHnn9U4y/9Hz85EbOvkVETLSmqjp4AvNWlOOi3ufiiemj2MOAiAA4ONFG\nu3tXTW0D7n6+FJ0TOmHZd7JxbhdeZyQiD8cm2mh272psasa9K8ux52gtnr1jNFevJaIW2Owywa/e\n+AzvbDuEvCnDMXZQitXhEJHNxFyLVkSyrI7B18oPd+OP/9yJ718zCNPHDLA6HCKyoZhKtCKSA2CN\n1XF4vV91GA8Xbsa1Q9Ow6JsXWx0OEdmUpYlWRPKM5Om/f4GI5Br/nmnBqmoJgKqoBtmK3YdrMXd5\nGS5M7Ybf3p6Jzgkx9ZlFRFFkSY3WSK5ZAHIBFPsdWwNgsaqWG9vFAMZHPcg2HK9vwN0vfIRmBZ77\n7hj06JpodUhEZGOWNMNUtURVlyBw6zTHm2QNVYFavVZpalb88OUN2H7wJJ6ekYVB551jdUhEZHO2\n+r5rJFT/5OuGjVq0S4q24K0tB/BfEy/F1UPOszocIooBduve5Qqw7zCAMdEOxFdhRTXyiypR7a4D\nAFw9OAV3jhtoZUhEFENs1aIF0GYnVBHJBZBhXCQLlJRNV1hRjUXrNp1JsgBQttuNworqaDw9ETmA\n3RLtkQD7Ur03VHWtqg5W1SWq6o5GQPlFlahraGqxr76hGflFldF4eiJyALuVDtwIXD4IukuXiMwC\nMMu7PWBAxwYR7PNpyQazn4jIn61atEY/Wf/yQQb8uoC18xgFqprt/UlLS+tQTH1dgVdGaG0/EZE/\nWyVaQ4nfMNsMIwGHxKzZu+ZPGIbkxJZrfSUnJmD+BM41S0TBsSTRikiWiOQByAGQJyILfA7PBDDd\nGBmWZ2yHzKzZuyZlpmPx5OFIdyVDAKS7krF48nCunEBEQRNVtTqGiBCRiQAmDhkyZOa2bdusDoeI\nHEZEylQ1O5j72rF0YAo7LDdORAQ4ONFGe4UFIqLWODbRskVLRHbh2ERLRGQXjk20LB0QkV04NtGy\ndEBEduHYREtEZBd2m+vANN5+tABqReQzv8M9AdQEue17+zwAh0wK0f85w71va8cC7W9vXyydd1vH\ngzn31s7Vf9usc7fjeftv2/01t9t5Xxjk/QBVdfQPgIL29rW17Xe7NJJxhXPf1o4Fc97tnKutz7uj\n597auUbq3O143rH2msfaefv+xEPpYH0Q+9raDvT7Zgjlcdu6b2vHgjlv/32xdN5tHQ/1NW/v/WAG\nO563/7bdX/NYO+8zHDsENxJEpFSDHHLnJPF63kD8njvP21zx0KI1U4H/Dr+ZxpyqxXmLyCzjZ6lV\nAUWR/7nnGD950VrlwyJnvdcBwJjoycn8X++lIuIy3u9hv95s0XaAsZjkUlUdbHUs0eJdQFNVq7yz\nrqlnRWPHMz5UF6nqVONDpkxVAyYkJzLOf5mqjrY6lmgRkaPGzYUdea3ZovVhtFLOWtrcWKMs1/j3\nTAtWPfPkBr36g12FeN4ZAHKN21UAYvpDJpRzV9VyVZ1q3CUDQMjzJNtFqO91QwoCLzcVM8I475mq\n2qujH6iO7d4VCuM/PgueBFLsd2wNgMWqWm5sF8NGy593RDjn7feGG+//e7GiI6+5sVzSGlWNuQ/Z\ncM/bSD6l0Y3WPB14vTO8v9uRb25s0cLTMjX+EwP94eR4XwBDVaBPxFjUkfMWkQzjMdZGOMyI6Mi5\nGx82o2PxfdCB807RKC2IGgnhnrd6FoItAeD2W6AgJEy0bfDWI/12u+GQFm1rgjzv2ao6O3pRRUdb\n526sDOL9WlkGwDHn3955I4Zbs21p57xzjW8vQAfLZEy0bQt0lfEwPPU5J2vzvEUkF8Bi43bMtera\n0da55+Cr194FB9TnfbR13hkAcozXPcNhr3lb510FYLWxLwOeD9ewMNG2zX9F3hZ83ngLHNbVp9Xz\nNv7IlgEoM67IOq17W6vnbnz1TDFe98EwPmwcoq3zXutTImrzbyIGtXXe5QCmeV/vjlwQ48WwtgW6\nwprqvWG8+WKyRtmOVs/bqFf1im44UdXea+79Y3Pa697meQOOfb9H5fVmi7ZtbgT+auGkr4yBxOt5\nA/F77jzvlkw9bybaNhitN/+vFhmI0S5NwYrX8wbi99x53i2Yft5MtO0r8evAnGG8OE4Xr+cNxO+5\n87w9TD9v1mhxpjP2dHiuKqeIyCqfzskzASwy+o2OMbYdIV7PG4jfc+d5W3PenOuAiCjCWDogIoow\nJloioghjoiUiijAmWiKiCGOiJSKKMCZaIqIIY6IlIoowJlqKayKSYSzAp8ZM+/7H13iPxclCnBQB\nHLBAcc8YEbQQwCwAvfxXEhCRPFVdaElw5Ahs0RJ55tRdCM9MTrN8DxjzDH9kRVDkHEy0RF+th1WA\ns5enyQZQfvavEAWPiZboK0vhWTHDfyYnp8/JShHGREtxzSgNHAEAI6GWA1hkaVDkOEy0FO9yAPjO\nPboYQK5FsZBDcT5aincpvr0MVHWtiMBYZroKLZMwUVjYoiU6m/eiGOuzZAomWopbvvVZP0vh6fI1\nOLoRkVMx0VI8m4YAXbdUtRyessH2qEdEjsRES3FHRFzGcNulAJYaI8P85YH1WTIJh+ASEUUYW7RE\nRBHGREtEFGFMtEREEcZES0QUYUy0REQRxkRLRBRhTLRERBHGREtEFGFMtEREEfb/AQa785EhNi1V\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111dcd1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(Ns,numpy.nanmean(ncoll,axis=1),'o-')\n",
    "gca().set_xscale('log')\n",
    "gca().set_yscale('log')\n",
    "xlim(5,150000)\n",
    "#ylim(0.,2.)\n",
    "xlabel(r'$N$')\n",
    "ylabel(r'$\\#\\ \\mathrm{of\\ collisions} / t_{\\mathrm{dyn}}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The full, exact algorithm to run a system for a dynamical time therefore runs in order $N^2\\,\\log N$ time, which becomes prohibitively slow for large $N$.\n",
    "\n",
    "## Approximate solver\n",
    "\n",
    "``wendy`` also contains an approximate solver, which computes the gravitational force exactly at each time step, but uses leapfrog integration to increment the dynamics. The gravitational force can be computed exactly for $N$ particles using a sort and should therefore scale as $N \\log N$. The number of time steps necessary to conserve energy to something like 1 part in $10^6$ should be relatively insensitive to $N$, so the overall algorithm should scale as $N\\log N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ns= numpy.round(10.**numpy.linspace(1.,6.,11)).astype('int')\n",
    "ntrials= 3\n",
    "T= numpy.empty((len(Ns),ntrials))\n",
    "dE= numpy.empty((len(Ns),ntrials))\n",
    "E= numpy.empty((len(Ns),ntrials))\n",
    "tdyn_fac_norm= 3000\n",
    "for ii,N in enumerate(Ns):\n",
    "    tnleap= int(numpy.ceil((1000*tdyn_fac_norm)/N))\n",
    "    for jj in range(ntrials):\n",
    "        x,v,m,tdyn= initialize_selfgravitating_disk(N)\n",
    "        E[ii,jj]= wendy.energy(x,v,m)\n",
    "        g= wendy.nbody(x,v,m,tdyn*(tdyn_fac_norm/N),approx=True,\n",
    "                       nleap=tnleap,full_output=True)\n",
    "        tx,tv, time_elapsed= next(g)\n",
    "        T[ii,jj]= time_elapsed*(N/tdyn_fac_norm)\n",
    "        dE[ii,jj]= (E[ii,jj]-wendy.energy(tx,tv,m))/E[ii,jj]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm indeed scales close to $N\\log N$ with a fixed time step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x111d98390>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWIAAAFKCAYAAAAno1UhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8VPWd//HXFwSMN4YAXqBFnKj0ZqtD6Nra33Yrg7Re\nKGoC3vBakna73Xa7JZF2218vu8WE2ovr/vpLqPcLQiKiaBUz+rN1ta0mwS5WS23iFWq1kFFroyB8\nf3+cM3EymSQzmcs5M+f9fDzmQebM7fNlks9853O+F2OtRUREvDPO6wBERIJOiVhExGNKxCIiHlMi\nFhHxmBKxiIjHlIhFRDymRCwi4jElYhERjykRi4h4TIlYRMRj+3kdgNemTZtmZ8+e7XUYIlJmurq6\n/mKtnZ7JfQOfiGfPnk1nZ6fXYYhImTHGPJ/pfVWaEBHxmBKxiIjHApuIjTFnGGNaX3vtNa9DEZGA\nC2wittZustbWTZ482etQRCTgApuIRUT8QolYRMRjSsQiIh5TIhYR8VjgJ3SIiIxk45btrN68jR3x\nfmaEKlixcA6LT5iZ19cIbCI2xpwBnHH00Ud7HYqI+NTGLdtZuWEr/Xv2ArA93s/KDVsB8pqMA1ua\n0PA1ERnN6s3bBpJwQv+evazevC2vrxPYRCwiMpod8f6sjo9VWSViY0yLMSZkjKkzxoS8jkdEStuM\nUIX7kx3meH6UVSIGlgDPAlhr4x7HIiIlbsXCORw9YScbJ36LY82LAFRMGM+KhXPy+jq+PllnjGkC\nOqy1sZTjDUAvEAZi1tpu96bl1tr2IocpZai3t5fKykpCoeG/WKW7T3t7O6FQiI6ODpYuXUokEilG\nuFIgiw/fyacP+A67336Lg+hnZoFGTfiyR2yMibrJtibNbW04ybfdWtsMNCXdHE56rMiYxGIx6uvr\nR1ynOt19ent7aWlpIRqNMm/ePFatWlWMcKVQen8B153K/pP255AvPsiGVf/CI5efnPckDD5NxNba\nmJtke9PcHE3qAQP0GmOi7uOa3d5zXMlYxioajRIOh7O+Tzgcpq2tDYDHH3+cBQsWFCxGKbAnN8At\nNTD5PXDZ/TA9v6WIVL5MxMNxE25qco4DC4wxNcaYOvdYL1BVzNji8Thz584lHvdXabq9vZ3GxsZR\n71dfX09vb7rPvXfF43EaGxtpb2+ntbWV7u7uEe/vF7FYjNraWowxxGKDqly0t7dTVVXF3Llzh9w2\nFqFQiPb2duLxOHV1dYNu6+7uprGxkebm5pxfZyTFbG9Z+vX/hfZLYeZcuPRemJz/HnCqkkrEQLqC\n3U6cWnEvsN49Fga6ihUUMPBLncnX0aqqwn9GxGIxmpubaWlpyejDIRaLUVVVhTFm0KW1tRVwkvD8\n+fNpamqipsapGJXKV+9oNMqCBQuoq6ujqalp0G01NTXU19fT1dVFNBrNy+vV1NSwYMGCIR+AkUiE\nqqoqOjo68vI6wyl2e8uGtRD7NtzXCO87DZbdARVTivLSpZaIK4e7wS1XLDHG1ABV1trWYgUVj8cJ\nhUI0NTXR3Nw8YuKLxWKjfu3Nh2g0SkNDQ8Yni6LRKF1dXfT09AxcGhoaBnp1jY2N1NfXD9y/rq6O\nNWvWFCT2QqmvrycWiw15f0Y6IZeN7u7ugQ/kSCSStudbjPc+odDtLSt798DGf4T//hHMvQSW3AgT\n8jtEbSSlloh3pTk2NfGDtbbVPYk37Hdxd4xxZ+Ly6quv5hxULBYjGo0SjUaJRCIj9hTb2toyqh02\nNzcP9KqGKxnkqzSQKDlEIhHC4TDhcJhYLMbKlSsH7tPa2jqkB1WMP+jE1/x0Er31TEUikSHvTzwe\np7Jy2M/3QUb7ZtHZ2TnwnsTj8aIm3XRybW9g7H4T1p4Lv70V/uHrcPqPYNz4ooZQaok4TvryxMjF\nzSRusq5OXKZPz2i36+EDcnvDCcP1ihOlgtbWVnbu3DlinbCxsZGdO3fS1NTEvHnzqK2tHfI1t7u7\nO2+7T4dCoUFJo7u7m3A4PNCuxAdBb2/vQH14pPh7e3uHfNgk18+bm5uJxWID9euRPlCi0SjLly8f\n8v9ZX1+fcULp7u6murp64HHJCbyzs3PIB0x7ezudnZ20tbUNim3u3Lkj3qeurm6gRtzS0jJw4m4k\n8Xic5uZm2tvbaW9vH/L/2tvbO1CXb2xspLW1deCSr/YG2ps74YYzoOcBOP3H8A+NYEzx47DW+vYC\ndOCMkkg+1pdyvS31Phk+9xlA69FHH21z0dLSMuRYOBy2DQ0NQ4739fXZUCg06nOme2xLS4sNh8O2\npqbGRqNRG4lEMoqvoaHB1tXVZXTfhNT7d3R0WMB2dHQMHGtqakobZ+K25Ofo6ekZaHdLS4tta2sb\n9NxdXV0jxtPX12ej0ajt6+sbiC/5OUaT+h4ltyXd+1dIHR0dNhqNDlyPRCID7Urcnvx/F4lEBv5/\nUh87HD+119d2PWftVRFrv3eotU/fnfenBzpthvmo1HrEADFjTHLhM2xTJnxkwuZh0Z/hvn6O1CtO\n9FRGklwSSKirq6Onp4f6+noaGxvp6irMucjESbtkiZ5ncuzRaHTYXnFHR8egHnGidANOjTTRs+vt\n7R0o54wkFArR1tZGbW0t9fX1LFiwYOCE4VikO4nlhUQ9OfkbVTQapbW1deB3J/HtBJz/u7F8C/JL\ne33l5a1wzQJ48y9w4Z3OyTkP+TIRG2Mi7qy6KNCUMiZ4ObDUHa7W5F4fy2vkvIvz+vXr037Nq6mp\nIRwOD6kVpyaokdTW1jJlyhTmzp1Le/u7kwUTtWhg0PF8aWlpGfLhkkgUyQkj8XO6umly4oXB7Y5G\no7S0tNDR0TEwjCqTUR2J8kk+vlonTmIlZsZ5pbu7O+3rh0KhgYRbU1MzkLBjsRhLlizJ+nX80l7f\nePaXcN2pMG4/uPQ+mHWi1xH5uzRRjMvcuXPH9LWjr69vxK/HbW1tFhj0tTMcDg98zezp6Rn2sXV1\ndQNfJXt6emxDQ8Ogr6jWDv2KP5xsSxNA2lJBKBQaFHNXV5d1fn0G6+rqsuFwOO1je3p6BpU3EvE1\nNTWNGlfi/yS1TDGarq6utO0Jh8NDygLFkFxeaGtrS1tiAgb+r9va2gbe60zeb7+113ee3GDtd6dZ\ne/VHrY2/WNCXosxLE3mRa484MYogHo+nvSRmXiWfZNu1axeRSIR4PD7i5Im5c+cO+irf1NREW1sb\njY2NVFVVMWXKFHp6enL6ep448ZYs0TNNNxpi5cqVgyYArFu3Lu3X3VgsNujxieFT4XCY7u7uQUO8\nAJYuXTpqrPX19dTW1hKNRgeVKTIdH52u9NHY2DjkRGux1dTUDPldaG9vH/hGBc4Mvbq6OmpqajJ6\nv/3cXs/9phXaLnEnatznzJrzi0wzdrlextIj7uvrszjr4mV0SfRCGhoaMu7Z5Kqrq8s2NTXZcDhs\nQ6GQbWpqGtKjTj3x09fXZ8Ph8LC9pqampkGXdCKRiI1Go7alpWXgRFzi5FpXV5dtaWkZ1MMbrTec\neJ5UfX19tqamZtjH9fX12bq6Oguk/UbQ19c37MnGQunq6rLRaNQCA+1OxJHo+ab+f7S0tNhQKGTD\n4bANh8M2Go2m/f3xY3t9Y98+a2PfsfZ/H2Ltredau/tvRXlZsugRG+f+wVVdXW3zNQxMwBhD0H+n\n8qW7u5t169axcuVKQqHQQO95+fLlrFmzRiu7ZWLvHtj0ZXjiFph7MZx6JYwvzqKTxpgua+3oZ+fx\n+TKYhaQ96/JvuK/FMjaxWIx58+YNOlkaiUQyGnkjOBM12i6GZ+6Hf1gJn/RojHAG1CNWjzhvEhMH\nUhe7kbFLDGVLHqUSDodzOj9QLkbcXfnNnXDrEtjRDaddCdWXFj2+bHrESsRKxCIlJ3V3ZXB2zlh1\n1nEsnv0O3Hw2xF+Ammvh/ad7EmM2iVijJnIYRywi3hhud+UN926Ga06BN19xJmp4lISzFdhEbPMw\ns05EvJFuF+UTxz3F1W9/Hcw4uOQ+OPJjHkQ2NoFNxCJSulJ3Uf7MuN9ww4Qr2DluKnyuAw77gEeR\njY0SsYiUnBUL51AxwVmqctn4+/mvCVfxO8I89en1/pqokSENX9PwNZGSs/iEmWAt8Xv+Nxfvbefh\ncfPoO7WFRdVF3SEtbzRqQqMmRErP3nfg7q/AlpsgciGc9qOiTdTIlCZ0iEj52tPvbO657efw9w3w\nqa/7dqJGppSIRaR09Pc52xq98Gs49Qfw0TGtgus7gU3EqhGLlJjXd8BNZ8GuHqi9Dj54ptcR5U1g\nR01oHLFICXn1D85EjddegvPbyyoJQ4B7xCJSIl7qhFtqYNwEuOQeOOIjXkeUd4HtEYtICXimw9ll\nef8QXLa5LJMwKBGLiF89sRZuXQrTjoHL7ofKoRv1lgslYhHxn0d+Ahs/D7M/ARfdDQcd6nVEBaUa\nsYj4x7590PFN+NXVzgm5M1tgv0leR1Vwge0RaxlMEZ95ZzfcUe8k4Y/Ww9nXBiIJQ4ATsYavifjI\n23+FtefA1vVw8jfhM00wLjjpSaUJEfHWm3+BW2rhT0/Aov901o4ImLL8yDHGNHkdg4hkoO95uHYh\nvPIULL0lkEkYyrBHbIyJAFGv4xCRUbz8pLO33Dv9zrZGs070OiLP+DoRuz3bDmttLOV4A9ALhIGY\ntbY76eZKYFfxohSRrD33iLN4z8QD4dLNcOj7vY7IU75MxMaYKBABaoCOlNvagFWJ5GuM6QAWuD9H\nAC0uLOJnT2+C9stgypFwwQYIvdfriDznyxqxtTZmrW3G6fWmiqb0gHvdxA1Qaa2NFz5CERmTzutg\n/YVw+HFOT1hJGPBpj3g4bsJNTc5xYIExZhfqDYv4k7Xwi2Z46PtwzClQe71TlhCgxBIxEEpzbCcw\nD6deHDbOSv1hY0w0tbYsIsWxcct2Vm/exo54P++ZPJEbjmgn/Nxt8JFznSFq4yd4HaKvlFoirhzu\nBmttO4Axpmak+4lIYW3csp2VG7bSv2cvk9jNyr/9mPBzj/GHoy/j2MVXlvy2RoVQaok43WiIqclX\n3ITcPtwTGGPqgLrE9VmzZuUtOBGB1Zu30b9nLwfzN1on/JCPjX+K7+05n/teOp1HlITTKrVEHCd9\neSLdSb20rLWtQGvienV1dbC3sRbJsx3xfqYT5/qJTRxrXuLLu/+RO/d9AhPv9zo03/LlqInhuDXf\n1LJDmJQhbpnQoj8ihTHvkDjtE7/NUeZlPrfna9y57xMAzAhVeByZf5VUInbF3PHCCeGxnJTToj8i\nBfCn33Kj+SaHmH7O2/0NfrHP2VGjYsJ4Viyc43Fw/uXL0oSbaJfiTFWuNMasc8cVAywHVhpjwjij\nJca0n7Z2cRbJs2cfhrXnsv/+k/nvk67n1UffwcT7mRGqYMXCOSw+YabXEfqWsTbYJdLq6mrb2anh\nxyI5eeouuP0ymHIULNsAk9/jdUSeM8Z0WWurM7lvKZYm8kI1YpE86boe2i5yNva89D4l4TEIbCJW\njVgkR9bCL1fDpi9D1XxnBbUDNIR/LHxZIxYRn9u3D+67HB5rgQ8vhc/+l2bL5SCwPWKVJkTG6J3d\nsGG5k4RP/CIs/r9KwjkKbCJWaUJkDN7+K6xdCk+2Q/TbsPA/ArW3XKGoNCEimXlzJ9xaCzu2wKKr\nIbLM64jKRmATscYRi2Qh/iLcdCa89iIsvRned5rXEZWVwH6nUGlCJEOv/B6uOQX++oqzo4aScN4F\ntkcsIhl48TFnq/v9JsElP4fDP+R1RGUpsD1ijZoQGcUf7ocbFjljgy+7X0m4gAKbiFWaEBnBb2+D\ntefAtGOcveWmzPY6orIW2EQsIsN49Gq4ox6O/DhcfA8cdKjXEZU91YhFxGEtxL4Nj/wY3r8IzloD\nE/b3OqpAUCIWEdj7Dtz9ZdhyM8y9BE67EsaN9zqqwAhsaUIn60Rce/ph/TInCX+yEU7/kZJwkQW2\nR2yt3QRsqq6uHtPC8iKlKHmb+xmhCr5+8hGc9uRX4YVfwak/gI/qz8ELgU3EIkGTvM09wO74Do6+\n58vsG7eDcTXXwIfO9jjC4ApsaUIkaBLb3APMNn9iw8Rv8x7+zFf3+4aSsMeUiEUCYoe7nf0HzbO0\nT/wOB5i3OHf3v3HnG9rU02sqTYgExIxQBe99vYs1E67kNQ7kwt2X02tnMFPb3HtOiVgkIH744Rc5\n/jdNPG8P5cLdl/MyU7XNvU8EtjSh4WsSKFtu5u8e/wpvVr6fL1d8nz8zlZmhCladdZy2ufcBY631\nOgZPVVdX287OTq/DECmcR34CHd+CqpNhyU0w6SCvIwoEY0yXtbY6k/uWVWnCGBN1f1wArLLWxr2M\nR8RT1joJ+NGr4INnwZktsN9Er6OSNMqmNGGMiQD11toYEAKWeBySiHf2vgN3/ZOThKsvg7N/piTs\nY75OxMaYpqRebvLxBmNMjftvBMBa222trXXvEgZixYxVxDf2vAVtF707ZVnrRvieL0sTbvKNADVA\nR8ptbThlh273egdOKSJxex3QZq3tLV7EIj7x1utw23nw3MPw6SY48fNeRyQZ8GWP2Fobs9Y2A+mS\naTSRhF29yb1ma20rMDddT1qkrP31VbjhdGfdiLPWKAmXkKx7xMaY44FKnDosQBzotdY+l8e4hnvt\nKEOTcxxYYIzZBU6JAugC6lF5QoIi/gLcuBhe3wHnrIVjT/E6IslCRonYGHM2ztd/CxigBycBAlQB\nS4wxuLd3WGs35D9U4N3kn2wnMA9IJOlu934qTUgwvPI03HQW7HkTLtwIs070OiLJ0oiJ2BhzAk6C\n67LWZvQ9xxgz3xjzNSBmrX0iDzEmqxzuBmttszGmzhhTg/Ph0Jjn1xbxnxcfh1tqYL/94ZJ74bAP\neh2RjMGwidjtBfdYa1dn84TW2geAB4wxJxhjTrbWPphrkEl2pTk2Nem1W90f2/P4miL+9McYrFsG\nBx0Gy+6AyqO8jkjGaNhEbK29PZcnttZuyeXxw4iTvjyRcRnCHVVRl7g+a9asPIQlUmRP3g4b6mH6\n++CC2+Hgw7yOSHKQ1agJY8xs93KIe/14Y8wVxpjPFSa8wdzJGqnliTApQ9xGeY5Wa2114jJ9+vS8\nxihScI//DNovg/fMg4vvVhIuA9kOX/s8Ts240q0ftwG3AV3FSsZALDGJwxV2E3RWtOiPlBxr4aEm\nuOdf4dhPw7INUJHuC6KUmmyHr3W4NWCMMU1Aa+KEnDFm2BNp2XIT7VLeTfrr3HHFAMuBlcaYMM5o\niTFtsqU966Sk7NsH910Oj7XAR86FRf8J4yd4HZXkSbaJeGfSz1EGJ8G8LePmjgXuJs3IB3chn8Tx\nMZ+UM8acAZxx9NFHj/UpRIpj7x7Y+AXY2gYnfhFO+XcY58u5WDJG2b6bVcaYQ4wxK3CGtL0OYIw5\nOf+hFZa1dpO1tm7y5MlehyIyvN1/g7XnOkl4/rdg4X8oCZehrHrE1trb3SQMUAtgjLkCZyRDJ5DP\noWoFpR6x+M2Qre4/dTinPfkVeOlxOP3HUH2J1yFKgeRtYXhjzPEFmMBRcFoYXvwgdav7Q+nj5klX\nUDX+ZcbXXAMf+KzHEUq2CrowvDFmNs6QsVSNwMJsn09EBm91f6R5mZsmrKKSN/iX8f/GVUrCZS+r\nROyWISKkn0CRLjn7lkoT4ieJre4/YJ7jholNjGcv5+3+BlvfruIqj2OTwsu2R/y4tfbydDe4U6JL\nhoaviZ/MCFVwxGtbuHbiav5KBefs/jd67ExtdR8Q2SbiYfeAy3VKtEiQ/eD4lzn+V1eww05l2e6V\n7GCatroPkGzHwfS66xEPYYxZlYd4ikYz68Q3trbzsce+xNuhKv654vv8iWna6j5gsho1YYyZD7QA\nRzG4d2yAydbaktsYS6MmxFOPX+NMWT7y43DuWthf49rLRSFHTUSAudbaId1I90SeiGTCWvjvH8ID\n34VjFsKSG2CC6sFBlW1pojtdEgYY7iSeX6k0IZ6xFjq+5STh42rhnFuUhAMu20RsE0tgpjLGnJWH\neIpGU5zFE/v2wl1fgkevgnmfgzNbtXiPZF2aOAVoNMbEGbxbhgHmA4Xaq06k9L3zNmxYDk/dCX+/\nAj71DXD2epSAyzYRR4HhRkeU1IQOkaLa/SasuwB6HoRT/gM+/k9eRyQ+km0iXj7cFkjGGO2aLJJO\nfx/csgS2d8KiqyGyzOuIxGeyXX1t2H3oCrRHXcFoirMUxRsvO1vd73wGam+ADyzyOiLxoWFP1hlj\nznYX+BkTY8xRfj6Bp5N1UnB9z8G1n3b+PW+9krAMa8RdnI0x840xUSBmrX0ukyd0k3ct0GOt1ck7\nCaZXnoYbF8M7b8FFd8F7MhrXLwE1YmkiaX+6Fe4ecXGc7ZJ6eXdmXQjnRN009+c/WmtXFyxiEb97\nqRNuqYHxk+CSe+GwD3gdkfhcRjXiRGI1xkwGqnESb5V7cxzYkkjaIoHW+xCsPQ8Omg7LNkLlUV5H\nJCUg25N1rwEPuBcRSfb0Jmi/FKYe42x1f/DhXkckJUK7EIrkw5abYf2FcMRH4OK7lYQlK4FNxFpr\nQvLmV/8Fd34RjvokXHgnHFDpdURSYgKbiDV8TXJmLTz477D5687mnuetg4kHeh2VlKCsNw8djTHm\nEGvt6/l+3gxfu879ca61tt6LGCQg9u2Dexvg8TVwwjI44ycwruSW4xafyDkRu6uxRYHE97EFwNJc\nn3cMcSTGO/caYxqMMQ3W2uZixyEBsHcPbPwCbG2Dj38JFnxPi/dITvJRmmjGGco2JemSF8aYJjfB\nph5vMMbUuP9G3MNhoMb9uZd3h9eJ5M+efrjtfCcJz/+WkrDkRT5KE23JY4iNMe25PqGbfCM4ibUj\n5bY2YJW1ttu93gEssNa2Jt1tQerjRHL21muw9lx4/lE47Ycw7zKvI5IykY8esTXGHG+MOcQtU5yd\n8xNaG3PLCulWdIsmkrCrN7nX7M4AxFqb8weCyIA3/wI3nAEv/gbO/pmSsORVPnrE7cDjOIvDg7Ox\n6A/y8LxDuAk3NTnHcXrAMfd6vU7USV69vgNu/CzEX4Bz1sKxp3gdkZSZfCTi2pTSxAl5eM7hhNIc\n2wnMc1+7BnfhemNM1FobS3N/kcztetZJwn/bBRfcDrM/4XVEUobyUZoYlHgLvC7xsCPl3d7yGqDL\nGNOHU2MWGbtXfu8sY/n26zx04s846bbdHHX5PZx0xYNs3LLd6+ikjOSjR/xRd93hXmvtE3l4vpHs\nSnNsKjh1ZfI4YkMCbscWZ0H38RN44MTr+afYW/Tv6Qdge7yflRu2ArD4hJleRillIucesbV2ibvu\ncK8xZlWBF4OPk748kfE2TcaYOmNMZ+Ly6quv5i86KQ/PPwo3LIKJB8El9/KtX+2jf8/eQXfp37OX\n1Zu3eRSglJucE7Ex5mQ3+T6IkySfdXf3OD7n6FK4vd7U8kSYLIaqWWtbrbXVicv06dPzGqOUuD/G\nnJ7wQYfBpffB1Cp2xPvT3nW44yLZykeNuB0nGc631n7BWrvFWns7I9RzcxRLmsQBEB7LSTkt+iND\nPHUX3HoOTDvaWdB9slN2mBGqSHv34Y6LZCsfibjWWvsDd61iYGDkxNyxPqExJmKMacKZOt1kjGlI\nunk5sNSdWdfkXs+aFv2RQZ64FdoughknwEV3Owu7u1YsnEPFhMHrSFRMGM+KhXOKHaWUKWOtze8T\nerjoTzaSdnFe/swzz3gdjnjpN61w7wpnGctzboVJBw25y8Yt21m9eRs74v3MCFWwYuEcnaiTERlj\nuqy1GW1WOKZEbIw5eYSb6621RV/0Z6yqq6ttZ2en12GIVx6+Eh74Lsw5DWquhQn7ex2RlIlsEvFY\nh69dDnThzKYL8+6ohRDpRzX4TlKP2OtQxAvWQuzb8MiP4bhaWPxTGD/B66gkoMaaiOuttc8CGGPm\np8ysm5+XyArMWrsJ2FRdXT2mGrOUsH374Odfg85roPpSOPVKGBfYPRLEB8b025dIwq7Us10lcfZL\noyYCau87zlrCndfASV92VlFTEhaP5WNm3VRjzCqgB2cN4J15eM6CU484gN5529ll+fd3w8nfhP/1\nr1pLWHwh50RsrV1jjDkbZ+Gdx90xxCL+svtNZ0H33v8Hn2mGv9MCfeIfWSViY8xsa+1zqcfd5FtS\nCVgn6wKkPw63LoWXHoPP/h844XyvIxIZZMTiWNJi7wlDuhHGmPnGmM/lPbIC04SOgEgs6L69C2qu\nUxIWXxrtLMU8IG6M2WmMWQdEjDFHJt/BWvuAtfZnpZiMpcy9th2u+wz85Q9w7lr44GKvIxJJa8RE\n7CbZcThTjWM4J+O2GGP2GmM2G2O+ltgmiRIZPywBsasXrvs0vP4nuGADHLPA64hEhpXRuB13IZ81\nQLu1thI4Bmexn6Pdf58d6fF+pOFrZeyVp+Haz8Dbb8BFd8Hsk7yOSGREWU1xTp28UQ40xbnMDCzo\nPhGW3QGHfcDriCSgspniPGyPON16EuWWhKXMPP8oXH+Gs2jPpfcqCUvJGGn42hRjzHqcCRod7i4c\nIv70TAzWXQCh98KyjQNrCYuUgmETcfLYYHfHDSVl8aen7nJmzB36PrjgjkFrCYuUgkxP1t3u7k33\nBcAYY9YbY346ynKYvqaTdWXif9qg7eK0C7qLlIqcFoZ3pzYvxekpt1lrH8xXYMWik3UlrPsmuOtL\nMPsTcO5taRd0F/FKMdYjBt4tXxhjJgPRpPJFSSZlKSGPrXGWsqyaD0tvhokHeB2RyJjlY/U13P3q\nkpPyEmPM54Eea+3KfLyGyIBHroKObzq7atReB/tN8joikZzkJREnc5PyGvcikj/Wwi+a4aHvwwfP\nhLPWaFcNKQtZrYhtjJldmDBERmEtPPAdJwl/5Dw4+xolYSkbgV19TUqItXDf5fDfP3K2Nvrsf8G4\n8aM/TqREjFaamAd0GGP6cBb9CRljjrTWPp+4Q2K2nTHmc9banxUu1PzSesQlYt8+uPsr0H0DnPiP\nsPD7g3bV0Db3Ug4Cu/qa1iMuAYn95bpvcLY1SpOEV27YyvZ4PxbYHu9n5YatbNyy3buYRcag7FZf\nM8ZEvI4Q3IahAAATF0lEQVRB8uCd3XD7pfA/t8HJ/wbzvzVkf7nVm7fRv2fvoGP9e/ayevO2YkYq\nkrNsR010AFhre4FefDYywhgTBVpweu5Sqva85cyW+8O9Ti/4Y19Me7cd8f6sjov4VVajJkZafS3l\npF5eGGOa3OSaerzBGFPj/jvQA7bWxnA+IKRU7f4brD3HScKn/XDYJAwwI1SR1XERv8oqEafjjqy4\nAujLQzyJ54waYxqAmjS3tQExa227tbYZaMrX64rH3n4DbqmBZ38Bi38K8y4b8e4rFs6hYsLg0RMV\nE8azYuGcQkYpkndjTsRJCfhZ8nyizlobc5Nsut5t1FrbnXS9N12vWUpMfxxuXAwv/BrO/hkcf96o\nD1l8wkxWnXUcM0MVGGBmqIJVZx2nURNScsY0s84YswJYCawDwtba14wxy/MaWfrXjTI0OceBBTij\nOqQUvbkTbloMr/4elt4E7zst44cuPmGmEq+UvGxn1q0wxjwDVAJHWWu/4E5pLpZ0Pe+dQLiIMUg+\nvfEyXH+as9PyOWuzSsIi5SKjRGyMWW6M+SNOAq621q4scgJOqBzpRmNMDRB2T+KV1LjmQHrtJbju\nVIi/AOe3wTGqMEkwjVqacJe2fMxa64cpaLvSHJua+MFa244zrln8btezcOMipza87A6Y9XdeRyTi\nmVETsbV2STECyVCc9OWJjIesGWPqgLrE9VmzZuUhLMnKX56BGxbBO/3OdvczTvA6IhFPZVsjnlyI\n8cKZcscJp5YnwrgTTTJ8jlZrbXXiMn26ttYpqj//Dq77DOzb42xtpCQskvWEjtdw9qyb7WFCjqVM\nYw67CTor2rPOAzuecE7MjdsPLv45HP4hryMS8YWsxxFba1+z1j6Hk5DnF2KNYmNMxBjThLPYUJM7\nuSNhObDUnVnX5F7Pmhb9KbIXH3PKERMPhkt+DtOP9ToiEd/IafNQAGPMUUAEZ5+6nGfqFUvSMpjL\nn3nmGa/DKW/PPgy3LoWDD4ML74LQe72OSKTgstk8NOfEaa191t1EtKQW2lGPuEj++IAzbTn0Xrjk\nXiVhkTTy1oO11vpmKcxMqEZcBM90wNpzYeoxcPE9cPDhXkck4kslU0rIN/WIC2zbfXDbeXDo+5wh\nagdO8zoiEd8KbCJWj7iAfv9zWHcBHPZBuPBOOGDECZEigRfYRKwecYE8fTesvxCO+DAs2wgVU7yO\nSMT3ApuIpQCeugvaLoIjPuJMW67Qch8imQhsIlZpIs9+t9HZ3mhGxEnC++ubhkimApuIVZrIoydv\nh/ZL4T3zYNkG2N+zWfAiJWlMC8OLDNjaDhuWw3tPhPPXw6SDB27auGU7qzdvY0e8nxmhClYsnKNF\n3EXSUCKWsfvtOtj4eZj1cThvHUw6aOCmjVu2s3LD1oHt7rfH+1m5YSuAkrFIisCWJlQjztETa+GO\nejjyJLcnfNCgm1dv3jaQhBP69+xl9eZtxYxSpCQENhGrRpyDLbfAxi9A+JNw3nqYeOCQu+yI96d9\n6HDHRYIssIlYxqj7Jrjzi1D1KTj3Nph4QNq7zQhVZHVcJMiUiCVzXdfDXf8ER893NvqcMHxSXbFw\nDhUTxg86VjFhPCsWzilwkCKlRyfrJDOPXwP3fBWOOQWW3AQT9h/x7okTcho1ITK6wCbipPWIvQ7F\n/x5bAz//Ghz7aVhyI+w3KaOHLT5hphKvSAYCW5rQyboM/abFScJzTs0qCYtI5gKbiCUDv/4p3NsA\n7zsdam9QEhYpECViSe/Rq+G+y+H9i6D2ethvotcRiZQtJWIZ6pGfwP3fgA8shpprYfwEryMSKWuB\nPVknw3j4h/DAd+BDZ8OZrTBevyIihRbYHrGmOKfxy9VOEj6uVklYpIgCm4g1aiLFL5rhwX+HDy+F\nM1uUhEWKSH9tAg9dAQ+tgo+cB5+9GsaNH/0xIpI3ZZWIjTE1QByIWGubvY7H96yF//d9+GUzHH8B\nLLpKSVjEA2VTmjDGRIBKa20M6DbG1Hkdk69ZCw9+z0nCkQth0X8qCYt4xNeJ2BjTZIyJpjneYIyp\ncf+NuIejwC73515gQbHiLEkPrYKHr4S5F8PpP4Fxvv5VEClrvvzrM8ZEjTENQE2a29qAmLW23S0/\nNLk3VeGUJcBJyNpCeDi//AH8oolN46OEH4lyUvNDbNyy3euoRALLl4nYWhtzk2xvmpuj1trupOu9\naXrNlYWLrsQ9ejU8+D3u3Pe/+PKbF7OPcQPbGCkZi3jDl4l4OG7CTU3OcZwyRA+De8FxZLDH1sD9\n3+CBcSfx1d117Et6+7WNkYh3SioRk77csBMIAzH3X9x/O4oVVEnout5dRe00Pv+3evYy9MSctjES\n8UapJeJhSw5uuSLu9poj1trW4oXlc0+shU1fcRZ1r72OQ0MHp72btjES8UapJeJdaY5NTfxgrW1N\nqi8LwJO3w53/6Gz0ueQm2G+StjES8ZlSm9ARJ315It1JvbTc8cUDY4xnzZqVh7B86ulNcPtymPUx\nd485Z3sjbWMk4i8llYittTFjTGp5Igy0ZPEcrcBA2aK6utrmKTx/2XYftF0CM+fCeeuG7LasbYxE\n/KPUShMAsaRJHABhdzZdVsp69bU/PgDrl8HhH4IL2mFS+pqwiPiDLxOxMSZijGnCmS3X5E7uSFgO\nLHVn1jW517NWtquvPfsw3HYeTJsDF2yA/cusfSJlyFhbnt/MR5O0i/PyZ555xutw8uOFX8NNZ0Ho\nvXDxPXDgNK8jEgksY0yXtbY6k/v6skdcDGXXI36pC26ugUOOgAvvUhIWKSGBTcRlVSP+02/h5jPh\nwKlw0SY4+DCvIxKRLAQ2EZdNj/jPT8GNi2HSIU4SPmSG1xGJSJYCm4jLokf86h/gxkWw3yS46C4I\nlfGYaJEyFthEXPI94p09cMMZgHF6wpXhUR8iIv5UUhM6xNX3PNywCPbtcUZHTDvG64hEJAeB7RGX\nbGnite1OT3j3G7BsIxz6fq8jEpEcBTYRl2Rp4o2XnSTc3wfL7oAjPux1RCKSBypNlIq/vuqUI954\n2UnCM+d6HZGI5IkSsQ9t3LJ90MpoX//UYZzWXQfxF5y1I2b9ndchikgeBTYRJ01x9jqUQTZu2c7K\nDVvp37MXgDfif+HIn/8re8e/xPjz18PsT3gcoYjkm2rEPqsRr968bSAJH0g/109s4liep3H8Cqj6\nlMfRiUghBDYR+1Vi37gK3uLaiav5sOnlS3v+mdvf+KDHkYlIoSgR+8yMUAWT2M2aCVdSbbbxlT1f\nZPO+edpPTqSMBbZG7FcrTjmaA+68jE+Y3/HV3Z/n7n0f035yImUusInYlyfrrGXx9h+CeYwf73cp\nd7z198zUfnIiZS+wC8MnVFdX287OTq/DcDx0BTy0Ck76Ciz4jtfRiEgOtDB8Keq81knCHzkPot/2\nOhoRKSIlYj94ehPc869wzCmw6CowxuuIRKSIlIi99twj0H4ZzIhA7fUwfoLXEYlIkSkRe+nPv4O1\n58KUI+H8Nph4oNcRiYgHApuIPV8GM/4C3Hw2TDzA2fb+gEpv4hARzwU2EXs6xfnNnc6297v/Bhfc\nDqH3Fj8GEfGNwI4j9szuN+HWJU6P+MKNcJimLosEXdn1iI0xEa9jGNbePdB2Mezohppr4ciPex2R\niPhAWfWIjTFRoAWoKsTzp64TnNWMN2vhrn+GZ+6H038M7z+9ECGKSAnyNBEbY5qADmttLOV4A9AL\nhIGYtbY7k+ez1saMMb35j3ToOsHb4/2s3LAVILNkHPs2/PZW+IevQ/UlhQhRREqUJ6UJY0zUTbY1\naW5rw0m+7dbaZqCp6AGmkbxOcEL/nr2s3rxt9Af/+qfwyI+h+lL4ZEOBIhSRUuVJIrbWxtwkm673\nGk3pAfe6JQdPJdYJzvT4gK3tcN/l8P4z4NQfaNaciAzhqxqxm3BTk3McWADEjDF16R5nrW0tdGwz\nQhVsT5N0R1wnuOdBuOPzcORJcNbPYNz4AkYoIqXKV4kYCKU5thOYB8VJuMNZsXDOoBoxMPI6wTu2\nwLplMO1YOOdWmLB/kSIVkVLjt+FrOU0vM8bUAGFjTIMxJl1SH7PFJ8xk1VnHMTNUgQFmhipYddZx\n6U/U7eyBW2qhotKZsFGR11BEpMz4rUe8K82xqZk+2FrbDrTnL5zBFp8wc/QREm/8GW4+C/bthWUb\n4JAjChWOiJQJvyXiOOnLE3kbkubWmQdqzbNmzcrXU8Nbr8MtNfDXV+CiTTDtmPw9t4iULV+VJtzx\nxKnliTDQkcfXaHVXzf8O0D1x4sT8PPE7b8O6C+CVp2DJjfCejBbmFxHxVyJ2xVKmKYdTJ3zkQ14X\n/dm3D+6oh2d/AYuuhmMW5P6cIhIYnpQm3ES7FIgClcaYde64YoDlwEpjTBhntMTyAsWQn81DrYXN\nK+F3d8CC78Lx5+YlPhEJDm0emuvmoQ//EB74Dpz4RVj4H5qwISKANg8tni03O0n4uFo45d+VhEVk\nTAKbiHPeoWPbfc5qalUnw2f/D4wL7H+liOQosNkjp5N1Lz7mrCt8+HHOCIn98jTyQkQCKbCJeMw9\n4le3OTtsHHIEnN8Okw4uTIAiEhiBTcRj6hHv3QNrz4FxE5wNPw+aXrgARSQw/Dazzt/GT3B216gI\nQeVRXkcjImUisIl4zOOIw58sSDwiElwqTeRjZp2ISA4Cm4hFRPxCiVhExGOBTcQ5T+gQEcmTwCZi\n1YhFxC8Cm4hFRPxCiVhExGNKxCIiHgtsItbJOhHxi8AmYp2sExG/CPwOHcaYV4HngclAons82s/T\ngL+M8SWTny/b+6Q7nnpspOuJn5OPlWJb8v2ejBRnJvfJti1+/f0a7rZSbIsf/laOtNZmtjKYtVYX\n58OoNdOfgc58vE6290l3PPXYSNeT4k8+VnJtyfd7Uuy2+PX3q5za4re/ldEugS1NpLEpy5/z8TrZ\n3ifd8dRjI13fNMx9xsqrtuT7Pcn0efLVFr/+fg13Wym2xW9/KyMKfGliLIwxnTbDTQH9rlzaUi7t\nALXFrwrZFvWIx6Y19YAxJuJFIHkwqC3GmDr30uJVQGOU2o6oe2kyxoS8CmqMhvx+ARhjmoodSB6k\nvi8txpiQ+ztW0u+LMSZsjKlxf8/CuTyxesR5YIyJAi3W2iqvY8mF245ea22vMaYBwFrb7HFYWXM/\nFFdaa2vdD5Qua23a5FYq3DatsdbO9TqWXBhj+twfG8vgPWlzf8dqgMpc2qMecRpuLyqa5niD+wnY\nkNwDttbGgN6iBpmhLNsSBmrcn3sB33ywZNMOa223tbbWvUsYiBUz1tFk+/vlqgR2FSfCzI2hLcut\ntVP8mISzaYubfHsBrLXtObenUGcBS/ECRIEGoAeIptzWBkSSrnek3N5R6PiK1Rb3WAtQU8rtAOqA\nOq/bkGtbgAgQ8tPvWA5taUg81us25NIW9/4t7nvTAIRyiUE94iTW2ph1voqn691GrbXdSdd70316\n+kUubUnUu6y17QUOc1S5tMM6vZS5fnmfcmhLpbU2XvgIMzfWtlhrm63zDTKeKH95bYxtmQrE3du6\ngZW5xKBEnIFE7TTlcBxY4EE4OcmwLfXW2vriRZW9kdphjIkkfR3uAkq6LUBn8aMam1HaUmOMqXOP\n+ar0lc4ofys97gWcklFOJ+uViDOT7uzuTpz6Y6kZsS1u7WuV+7MvepLDGKkdUd59b0L4tH6fZKS2\nhIGo+76Eff6ewMht6QXWu8fCOB+SfjZSW2K8+0ESxukVj5kScWYqR7ox6Y+koQSG5AzbFvePfA3Q\n5Z7d9vOQvGHb4X7NrHTflyrcDxYfG6kt7UklohF/D31ipLZ0A0sS74v14Qm7FCO1pRfocdsyz1rb\nmMsL7ZfLgwMk3dnqqYkf3D8Uz+upGRq2LW7tbkpxwxmz0d6TxB95KbwvI7YFSup3LDDvSz7boh5x\nZuKk/5ri96+86ZRLW8qlHaC2+FXR2qJEnAG3p5j6NSUMdHgQTk7KpS3l0g5QW/yqmG1RIs5cLGVg\neth9o0pRubSlXNoBaotfFaUtqhEncf/Dl+Kcda80xqyz707xXQ6sdMfYznOv+1a5tKVc2gFqi1/5\noS1aa0JExGMqTYiIeEyJWETEY0rEIiIeUyIWEfGYErGIiMeUiEVEPKZELCLiMSVikRG4G0S2GGOs\nMaYtze1tidtKeANZ8ZgmdIiMwp1V1Yiz9dKU1N0yjDFNuS6DKMGmHrHI6CI4iTiOk4wHuOtPP+5F\nUFI+lIhFRpfYM66VodsuVZPj7gwiSsQimWvB2YkldTWuUlxrV3xEiVhkBG7pYRcMbI+T8469IqmU\niEVGFsXZKDJhFVDjUSxSprQescjIKpNHSVhr240xuNvC9zI4SYuMiXrEItlLnLRTfVjyQolYZBjJ\n9eEULThD2qqKG5GUKyVikeEtIc3QNGttN05ZoqfoEUlZUiIWSWGMCbnTmVuAFndmXaomVB+WPNEU\nZxERj6lHLCLiMSViERGPKRGLiHhMiVhExGNKxCIiHlMiFhHxmBKxiIjHlIhFRDymRCwi4rH/D0Vl\nO1fMHVU/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111d5b128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(Ns,numpy.nanmean(T,axis=1)*10.**3,'o')\n",
    "p= numpy.polyfit(numpy.log(Ns),numpy.log(numpy.nanmean(T,axis=1)*10.**3/numpy.log(Ns)),deg=1)\n",
    "plot(Ns,numpy.exp(numpy.polyval(p,numpy.log(Ns)))*numpy.log(Ns))\n",
    "pyplot.text(10,10.**4.8,r'$\\Delta t  \\approx %.2f \\,\\mu\\mathrm{s} \\times N^{%.2f}\\log N$' % \n",
    "            (numpy.exp(p[1])*10.**3.,p[0]),size=16.)\n",
    "gca().set_xscale('log')\n",
    "gca().set_yscale('log')\n",
    "xlabel(r'$N$')\n",
    "ylabel(r'$\\Delta t / t_{\\mathrm{dyn}}\\,(\\mathrm{ms})$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, with the same time step, energy is much better conserved for large $N$ than for small $N$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x111ab56a0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAFLCAYAAABx8o8MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG61JREFUeJzt3T9vG1ma7/Hfs742QAywzZavEwvw7rAFaJLGwmo62nBp\nKDJg9FW3X4Hk2fxKvQo70sh6A2vNK+i2YQho3EAQZ/YFrCwFjgRfq/cGcuK1l96EgWE8N+ChWi6X\n+LdYdcj6fgDBZJFVPAeUfz5+6lQdc3cBAPL1N0U3AADKiPAFgAIQvgBQgP9RdAPGZWYrklqSltz9\nUdHtAYBBTPXI18yWJM25e1PSkZmtFd0mABhEdOFrZttm1kjZvmFmK+HPpbC5IeldeHwq6W5e7QSA\ncURTdgiBuyRpRdJB4rUnkrbc/Sg8P1AnaL+SdBTe9k5SNbcGA8AYohn5unsz1GxPU15udIM3OE0Z\nHc9NrnUAkK1owvcyIWSTgdxSZ+T7Sp+Odlt5tQsAxhF9+Cq9lPBWUk1SM/yp8OdBynsBIDrTEL6X\nlhNCKaLVrRe7+25+zQKA0UVzwq2HdynbrncfXAjcZq+DhGlo51PRfve7333zhz/8IZMGAkDX8+fP\n/9Pdb/R73zSEb0vppYe0E3OXCiF9PjKu1+t+eHg4ZtMA4FNm9v8GeV/0ZYdwAUWy9DByfdfM7pnZ\n7vv378duGwCMKvrwDZoXLqyQpFoIZQCYStGEr5ktmdm2OletbZvZxoWXVyU9CFe4bYfnI3H3X9x9\n7YsvvhizxQAwOivbzdTN7J6kewsLC6svX74sujkAZoyZPXf3er/3RTPyzQsjXwAxKF34csINQAxK\nF76MfAHEoHThCwAxmIaLLDJ14YTbQO/fOz7Tzv6JXrfaulmtaH15Ufdvz0+2kQBmXulGvsOUHfaO\nz7T57IXOWm25pLNWW5vPXmjv+GzyDQUw00oXvsPY2T9R+8PHT7a1P3zUzv5JQS0CMCtKF77DzHZ4\n3WoPtR0ABlW68B2m7HCzWhlqOwAMqnThO4z15UVVrl75ZFvl6hWtLy8W1CIAs6J0sx2G0Z3VwGwH\nAFkrXfgOO9Xs/u15whZA5kpXduAKNwAxKF34AkAMCF8AKADhCwAFIHwBoAClC1/u5wsgBjMRvonF\nNXtitgOAGEx9+JpZQ9KTotsBAMPIPXzNbDsEZnL7RlideGPIkWxT0mmmjQSACcvtCrcQuEuSViQd\nJF57ImnL3Y/C8wNJd/NqGwDkLbeRr7s33f2R0kepjW7wBqdpo2MAmBWF39shhGwykFvqjHybZraW\ntp+77066bQAwKYWHr6Rqyra3ku5IhCyA2RTDbIe5cXY2sxVJtXCiLi3IASA6MYx836Vsuz7ozu7+\nVNLT7JpTDFZJBsolhvBtKb30kOn0sVA7Pq8f37p1K8vDj6W7SnJ3sc7uKsmSCGBgRhVedgjzdJOl\nh5oS09Ey+Jxdd69L+lHS0bVr17I8/FhYJRkon8LDN2gmLqyohVDOXIyXF7NKMlA+eV5ksSTpgaSG\npDkz+ynM+5WkVUmbZlZTZ5bD6gTbMdQyQnm4Wa3oLCVoWSUZmF3m7kW3oRD1et0PDw+Lboakz2u+\nUmeV5K1vv6bmC0wZM3seSpw9xXDCLVcxjnxZJRkoH0a+AJChQUe+sZxwyw03UwcQg9KFb4yzHQCU\nT+nCFwBiULrwpewAIAalC1/KDgBiULrwBYAYlC58KTsAiEHpwpeyA4AYlC58ASAGhC8AFKB04UvN\nF0AMShe+1HwBxKB04QsAMSB8AaAAhC8AFIDwBYAClC58me0AIAZTH75mthZ+Hg/yfmY7AIjBVIev\nmTUkNd19V9IrM9souk0AMIjcw9fMtkNoJrdvmNlK+HNpwMPVJK2Ex6eSvsqqnQAwSbmtXhwCd0md\nsDxIvPZE0pa7H4XnB5Lu9jtmGPF23U0eFwBildvI192b7v5InRFqUqMbvMFp2uj4MmZWC5/xdMxm\nAkAuchv5XiaEbDKQW+qMZJtmtpa2X2LU+9DdH06oiQCQucLDV1I1ZdtbSXekz0L2M2a2ImkrPG64\nezPzFgJAxmII37lRdwyj5j9Lemdmc+qEMOF7ib3jM+3sn+h1q62b1YrWlxd1//Z80c0CSimG8H2X\nsu36IDuGUe6Xg7w3lC/OSxi3bt0aqHGzYu/4TJvPXqj94aMk6azV1uazF5JEAAMFiGGeb0vppYe0\nE3Mjc/ddd693f27cuJHl4aO3s39yHrxd7Q8ftbN/UlCLgHIrPHzD6DVZeqhpQtPGynp58etWe6jt\nACar8PANmokLK2qTOnFW1suLb1YrQ20HMFm5ha+ZLZnZtqSGpO3EpcCrkh6EK9y2w/NJtaOUI9/1\n5UVVrl75ZFvl6hWtLy8W1CKg3Mzdi25DIer1uh8eHhbdjFwx2wGYPDN77u71fu+LYbZDrszsnqR7\nCwsLRTcld/dvzxO2QCRiqfnmpqw1XwBxKV34lrXmCyAupQtfRr4AYlC68GXkCyAGpQtfRr4AYlC6\n8AWAGJQufCk7AIhB6cKXsgOAGJQufAEgBoQvABSgdJcXY/K4hwTQX+nCt8z3dsgDK2YAgyld2YET\nbpPFihnAYEoXvpgsVswABkP4IlOsmAEMhvBFplgxAxgMJ9yQqe5JNWY7AL2xjBAAZKg0ywiZWSM8\nvCtpy91bRbYHAAYx1TXfsNz8w7DMfFXS9wU3CQAGknv4mtn2hdHqxe0bYen4jRCqfbn7kbt/F57W\nJDWzbCsATEpuZYcQuEuSViQdJF57ok7J4Cg8P1CnjDDosdckPXH30+xaDACTk9vI192b7v5IUlpA\nNrrBG5ymjY57HHtX0jfD7AMARSr8hFsIzGQgt9QZ+TbDqPYz7r7bLU+E4H4u6aEoPQCYAoWHrzon\nypLeSrojnY9qL9MN7qNwHMoOAKZCDLMd5kbdMZQx5sxsRdJXkrYyaxUATFAMI993KduuD7rzhZHx\n017vC+WL8xLGrVu3Bv0IAMhcDCPfltJLD5mWENx9193r3Z8bN25keXgAGErh4RsukEiWHmpKTEfL\nCqsXA4hB4eEbNBMXVtRCKAPATMotfM1sycy21ZmhsG1mGxdeXpX0IFzhth2eTwQrWQCIQenuanbh\nlpKrL1++LLo5AGbMoHc1i6XskBtGvgBiULrw5YQbgBiULnwZ+QKIQenCFwBiULrwpewAIAalC1/K\nDgBiULrwBYAYlC58KTsAiEHpwpeyA4AYxHBLSWBoe8dn2tk/0etWWzerFa0vL+r+7fmimwUMjPDF\n1Nk7PtPmsxdqf/goSTprtbX57IUkEcCYGqUrO1DznX47+yfnwdvV/vBRO/snBbUIGF7pwpea7/R7\n3WoPtR2IUenCF9PvZrUy1HYgRoQvps768qIqV698sq1y9YrWlxcLahEwPE64Yep0T6ox2wHTjPDF\nVLp/e56wxVQrXdmB2Q4AYjAz4RvWfuuL2Q4AYtCz7GBm+5JOJb2SdOruz3Jp1ZDCyseNotsBAIPq\nV/O9q84y7v+RfMHM/knS7yX97O7/PegHhhHqQXJp+LCa8amkmqSmux8NekxJc5LeDfF+AChUv/B9\nmha8kuTuf5EkM1s3s7qkf3X3f7vsQGbWkLQkaUXSQeK1J5K2uoFrZgfqBH9fYdR7OMh7ASAW/Wq+\nb7sPzOy2mf198g3uviPpz5KaydcS72u6+yN1RrdJjcRI9zSE9SDm3L014HsBIAr9Rr7WfeDux2GU\n+y+SflandPAsvNY0s+NRGhBCNhnILXVGvk0zW0vbz913GfUCmFb9wtc/eeK+Y2Z33f2fU947aghW\nU7a9lXQnfOZuj31rkmpmpvBnI1lLBoAY9Ss71M3s7xLbLjsR9l8jtmFuxP3k7k/d/em4xwGAvPUb\n+X6jTv21pU5N90DSZRNk7ZLt/aTNUrg+zAFCAD/t+0YAiES/ke+upAVJm+qE6yNJD83srZn9ZGb/\n28z+IbzXLztIHy2llx7STsyNzMzWzOyw+/PmzZssDw8AQ+kXvo/d/Vd333X37919TtKXkh6qU2b4\no6RjM/soaWOUBoQabbJkUFNiOtq4Qh/qkn6UdHTt2rUsDw8AQ+kZvu7+2QwGd38faq1/dPcFdcL4\ngcYbqTbDzIWu2qROnHF5MYAYjH1XM3d/L+lpv3m5IVwfqHMZ8JyZ/RTm/UrSqqRNM6upM8thddx2\n9WjHPUn3FhYWJvURANCXuY9aqk0cyOyLEMRToV6v++EhU4QBZMvMnocSZ09Z3tUsmxSfMG4pCSAG\nY4evmf2tmf1Jo8/zzRU1XwAxGDl8L4Tur0qfKhYlRr4AYjBS+JrZuqT/UOeCi5q7/zHLRk0SI18A\nMRhqtkMI3TV1rib7/TSdYAOGtXd8xiKdmJiBRr5mtmpm/1ediyHq7r45rcFL2QGD2Ds+0+azFzpr\nteWSzlptbT57ob3js6KbhhnRN3zN7GdJX7j7wjSHbhdlBwxiZ/9E7Q8fP9nW/vBRO/snBbUIs6Zv\n2cHdv8+jIUBMXrfaQ20HhjXUCTcz+8LM/nZSjckDZQcM4ma1MtR2YFhDhW8oOZiZ/f20hjBlBwxi\nfXlRlatXPtlWuXpF68uLBbUIs2boezuEAH4fRsH/JOnVZYtsAtOqO6uB2Q6YlJFvrBNC+C9m9nsz\n+18a/WbqQJTu354nbDExWdzV7FdJv5rZVxm0Z+K4qxmAGGR2Y50QwtGj5gsgBlne1QwAMCDCFwAK\nQPgCQAHGPuEGYDzcwKecShe+zHZATLo38OneR6J7Ax9JBPCMm/qyg5k9NrOqma2ZWd+bujPbATHh\nBj7lNfXhK+l7dVbTkLu3Cm4LMBRu4FNeuYevmW2nLTNvZhtmthL+XBrikKvu/qW772bYTCAX3MCn\nvHILXzNrmNmGpJWU155Iarr7U3d/JGl7iEPXLhwbmCrcwKe8cgtfd2+GYD1Nebnh7kcXnp+mjY4v\nOe4jd29KahHAmDb3b89r69uvNV+tyCTNVyva+vZrTraVQOGzHULIJgO5JemupKaZraXt5+67ZrYi\naS6UHE4lfTfRxgITwA18yqnw8FX6svNvJd2ROiHbY99TSc3wuCbpebZNA4DJiCF850bd0d2PwhSz\nd5K+cvcfMmwXAExMDOH7LmXb9UF3vjAyfppNcwBg8mII35bSSw9pJ+ZGFmrH5/XjW7duZXl4ABhK\n4RdZhJkKydJDTdJBxp+z6+51ST9KOrp27VqWhweAoRQevkEzcWFFLYRy5ri8GEAMcis7hHB9IKkh\nac7MfgrzfiVpVdKmmdXUmeWwOsF2cGMdAIUzdy+6DYWo1+t+eHhYdDMAzBgzex5KnD3FUnbIjZnd\nM7Pd9+/fF90UACVWuvCl5gsgBqULX0a+AGJQuvBl5AsgBqULX0a+AGJQuvBl5AsgBqULXwCIAeEL\nAAUoXfhS8wUQg9KFLzVfADEoXfgCQAxiuJ8vgAnbOz7Tzv6JXrfaulmtaH15kXXjCla68OWuZiib\nveMzbT57ofaHj5Kks1Zbm89eSBIBXKDSlR2o+aJsdvZPzoO3q/3ho3b2TwpqEaQShi9QNq9b7aG2\nIx+ELzDjblYrQ21HPghfYMatLy+qcvXKJ9sqV69ofXmxoBZB4oQbMPO6J9WY7RAXlhECgAwNuozQ\n1I98w6KbS5Jakk7d/bTgJgFAX7NQ891296eSquqsjAwA0ct95Gtm25IO3L2Z2L4h6VRSTVLT3Y8G\nONZK2EchgAFgKuQWvmbWUKc8sCLpIPHaE0lb3cA1swNJdwc4bE1S1cyW1Bn17rp7K9OGA8AE5FZ2\ncPemuz9SGKkmNBIj3dMQ1v1cl9QK+x5J2sygqQAwcYWfcAshmwzkljoj36aZraXt5+67kl5d2PRO\nnZE1AESv8PBV50RZ0ltJd6TzkL1MU9LD8LimzugXAKIXw2yHuVF3DNPKXoUTb3fc/YfsmgUAkxPD\nyPddyrbrg+58YWTcc7ZDKF+clzBu3bo16EcAQOZiGPm2lF56yPRiCXffdfd69+fGjRtZHh4AhlJ4\n+Ib5vsnSQ02J6WhZYQFNADEoPHyDZpir21VLXoQBALMkt/A1s6VwdVtD0na4oq1rVdIDM1sJ71md\nVDtYyQJADEp3V7MLt5RcffnyZdHNATBjBr2rWSxlh9ww8gUQg9KFLyfcAMSgdOHLyBeYjL3jM/3j\nn/6q3//L/9E//umv2js+K7pJUYvhIgsAU27v+Eybz16cL1F/1mpr89kLSWK5okuUbuRL2QHI3s7+\nyXnwdrU/fNTO/klBLYpf6cKXsgOQvdet9lDbUcLwBZC9m9XKUNtRwvCl7ABkb315UZWrVz7ZVrl6\nRevLiwW1KH6lC1/KDkD27t+e19a3X2u+WpFJmq9WtPXt15xs64HZDgAycf/2PGE7hNKNfAEgBqUL\nX2q+AGJQuvCl5gsgBqULXwCIAeELAAUgfAGgAIQvABSgdOHLbAcAMZjq8DWzmpm5mb0KP8/77cNs\nBwAxmPYr3GrublIniNVZch4Aopf7yNfMts2skbJ9I6xevJFYRv5SieXlGyw3D2Ba5DbyDYG7JGlF\n0kHitSeSttz9KDw/kHR3iGPXJB1m11oAmKzcRr7u3nT3R5JOU15udIM3OE0bHffwMLE/AESt8Jpv\nCNlkILfUGfk2zWwtbT93373wdJigBoDCFR6+kqop295KuiN9FrKfMbO0/QEgajGE71wGx0grZQCY\nQXvHZ9rZP9HrVls3qxWtLy9O5X2EYwjfdynbrg+6s7u3JH2XXXMAxGqWlqiP4SKLltJLD5mOZs1s\nzcwOuz9v3rzJ8vAAcjBLS9QXHr5hbm6y9FBTYjpaBp+z6+51ST9KOrp27VqWhweQg1laor7w8A2a\niQsrapO6YILLi4HpNUtL1OcWvma2ZGbb6kwL2zazjQsvr0p6EK5w2w7PJ9UObqwDTKlZWqLe3L3o\nNhSiXq/74SEXxQHTJvbZDmb2PJQ4e4phtkOuzOyepHsLCwtFNwXACGZlifpYar65oeYLIAalC19q\nvgBiULrwZeQLIAalC18AiEHpwpeyA4AYlC58KTsAiEHpwhcAYlC68KXsACAGpQtfyg4AYlC68AWA\nGBC+AFCA0oUvNV8AMShd+FLzBRCD0oUvAMSA8AWAAhC+AFAAwhcAClC68GW2A4AYTH34hkU3G2a2\nNsj7me0AIAZTHb5huflWWGb+50EDGACKlnv4mtm2mTVStm+EUexGCNVBtNRZhr4qqS6J5YgBTIXc\nwjeUBjYkraS89kRS092fuvsjSduDHNPdTyU1Jf0qacndj7JsMwBMSm5Lx4fSQNPM7qa83HD37y48\nPzWzRtjnUmHEK0k/SHpsZkf99gGAfvaOz7Szf6LXrbZuVitaX17MfLn63ML3MqEEcZrY3JJ0V52w\nTq3juvuupDVJW+7eMrOmOiNmwhfAyPaOz7T57IXaHz5Kks5abW0+eyFJmQZw4eErqZqy7a2kO9J5\nyPbl7qdmdpBlwwCUz87+yXnwdrU/fNTO/snMhe/cqDu6+6Nwgu40HOfn7JoFoIxet9pDbR9VDOH7\nLmXb9UF3DifoACATN6sVnaUE7c1qJdPPiWGeb0vppYdkHXgsZrZmZofdnzdv3mR5eAAzYn15UZWr\nVz7ZVrl6RevLi5l+TuHhG2YnJEsPNUmZ1m/dfdfd65J+lHR07dq1LA8PYEbcvz2vrW+/1ny1IpM0\nX61o69uvZ2+2Q9A0s4vzdGuTmjLm7r9I+qVer69O4vgApt/92/OZh21SbuEbrlp7IKkhac7MfrpQ\nr12VtGlmNXVmOUwsGM3snqR7CwsLk/oIAOjL3L3oNhSiXq/74SFXIwPIlpk9DyXOngqv+eaNW0oC\niEHpwpdbSgKIQenCl5EvgBiULnwZ+QKIQenCl5EvgBiUdraDmb1R5+q6iyn8xYXn3cfdP/+npP8c\n4yMvHnvY96RtT25La/tlj8fpyzj9uOy1aezLsP1IPk/+fknT05dJfie92jnIe2Loy9+5+42+73L3\n0v5I2r3seffxhT8Ps/ysYd6Ttn2Qtvfo08h9Gacfs9SXYfvR7/drmvoyye9k1vrS66d0ZYeEX3o8\n/+WS92T1WcO8J237IG3v9XhU4/TjstemsS/D9iP5nN+vy81SXy5V2rLDsMzs0AeYOD0N6EucZqUv\ns9IPabJ9KfvIdxipN3UfYrHPmHzSl3DHtzUze1xUg8aQ7Esj/HQXVp0ml/2ODbSmYUSS38ljM6uG\n37Gp/k7MrBYW+m2E2yGMjJHvGMISSI/d/aui2zKq7jJO3lkJZEOa3nskh38IN939u/APyXMfcCWU\nWIU+/dndvym6LaMys/8KD3+Yge/jSfj9WpE0N05/GPleMOyy9t6581qm9x3OwpD9qOm3FaVPJUX1\nD8kwfXH3I/9tIdaaIlvPb9jfr2BO6QsOFGaEfqy6+5cxBu8wfQmBeypJ3llpfbz+TOpM3jT9qHOn\ntQ1Jr9RZSfnia0/UWZa++/wg8frBpNuXRz/CtseSVoruRwbfyZqktaL7MG5fJC2ps9BAFL9jY/Rj\no7tv0X0Ypy/h/Y/D97IhqTpOGxj5qjOC9c5/tdNGsQ3/7T7DUljWPqemDWWcfnTrV+7+dMLNHMg4\nffHOiOSbWL6nMfoy5+6tybdwMKP2w90feed/ia1uaatoI/bluqRWeO1I0uY4bSB8e+izrP3UGLAf\nD939YX6tGk2vvpjZ0oX/7j6XFHV/+vVF0lTc87RPP1bMbC1si66sldTn78qr8CN1SkFjnWwnfHu7\nbFn7sc5yFqBnP0Itays8jmK02EOvvjT023dTVYT1+IRefalJaoTvphb599KrH6f6bVXxmjr/KMas\nV1+a+u0fj5o6o9+REb699VzW/sJfjI3Ip9Bc2o/wl/rPkp6Hs9KxT527tC/hv5Fz4Xv5SuEflIj1\n6svTCyWgnr+HEejVjyNJ33e/E4/wpFtCr76cSnoV+nLH3X8Y54NiWcMtVj2XtQ9/OaKokfZxaT9C\nLe7LfJszln7fSfcv91R/L11T8jtWmu8ky74w8u0tl2XtczAr/ZDoS4xmpR9Sjn0hfHvwnJa1n7RZ\n6YdEX2I0K/2Q8u0L4dtfMzFhfGLL2k/YrPRDoi8xmpV+SDn1hZqv4lnWflyz0g+JvsRoVvohxdEX\n7u0AAAWg7AAABSB8AaAAhC8AFIDwBYACEL4AUADCFwAKQPgCQAEIXyBFWCjxsZm5mT1Jef1J97Up\nXUQVBeMiC+AS4QqnH9RZlujL5KoSZrY97m0FUV6MfIHLLakTvi11AvhcuH/zvxfRKMwGwhe4XHcN\ntV19viRRXWOuZIByI3yB/h6rs2JJ8k5X03i/WkSC8AVShLLCO+l8+ZixV6sFLiJ8gXQNdRZM7NqS\ntFJQWzCDuJ8vkG7u4uwGd39qZgrLoJ/q02AGhsbIFxhc98Qb9V6MjfAFEi7WexMeqzP97Kt8W4RZ\nRPgCn/teKdPI3P1InZLDq9xbhJlD+AKBmVXDpcSPJT0OV7glbYt6LzLA5cUAUABGvgBQAMIXAApA\n+AJAAQhfACgA4QsABSB8AaAAhC8AFIDwBYACEL4AUID/D1j1tUsRbZ27AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111d7f1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(Ns,numpy.nanmean(numpy.fabs(dE),axis=1)*10.**3,'o')\n",
    "gca().set_xscale('log')\n",
    "gca().set_yscale('log')\n",
    "ylabel(r'$\\Delta E$')\n",
    "xlabel(r'$N$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can design the time step such that energy is conserved to about the same degree for different $N$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ns= numpy.round(10.**numpy.linspace(1.,6.,11)).astype('int')\n",
    "ntrials= 3\n",
    "T= numpy.empty((len(Ns),ntrials))\n",
    "E= numpy.empty((len(Ns),ntrials))\n",
    "dE= numpy.empty((len(Ns),ntrials))\n",
    "tdyn_fac_norm= 10**4\n",
    "for ii,N in enumerate(Ns):\n",
    "    tnleap= int(numpy.ceil(((3000.*tdyn_fac_norm)/N)**1.))\n",
    "    tdt= tdyn*(tdyn_fac_norm/N)**.15/10.\n",
    "    for jj in range(ntrials):\n",
    "        x,v,m,tdyn= initialize_selfgravitating_disk(N)\n",
    "        E[ii,jj]= wendy.energy(x,v,m)\n",
    "        g= wendy.nbody(x,v,m,tdt,approx=True,\n",
    "                       nleap=tnleap,full_output=True)\n",
    "        tx,tv, time_elapsed= next(g)\n",
    "        T[ii,jj]= time_elapsed/tdt\n",
    "        dE[ii,jj]= (E[ii,jj]-wendy.energy(tx,tv,m))/E[ii,jj]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x1120ee2e8>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAFKCAYAAADooaOnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHAFJREFUeJzt3U9vFNme5vHnN1QtvOmbmOsNlqrrpltiVRphJ6tethFb\nVEPBK8C+vR9g/Aq4hhcwg/sVFKASUmsWyHl79tMJLGqFZux7N9TGDeSdjRcl9JtFnHSFgziZEZkZ\nkZHB9yOl5IxzIvIcbOLJc+KfubsAAMjznxbdAABAcxESAIAoQgIAEEVIAACiCAkAQBQhAQCIIiQA\nAFGEBAAgipAAAEQREgCAqK/q/kAzuy/pWFJXUt/dX09Tt4qyguX7kg7dvZ9a9kTSvrsfR/qxE358\nKmlV0q67P4j1O+33v/+9f/vtt0WqAkBhr169+g93X5tY0d1re0l6Jmkz9f5wmrpVlBVYd1vSfUlH\nkrYz632U5JnXx1T5/dTyI0ndov9mW1tbDgDzJmngBfZBdU83bfv5b+bHZrY9Rd0qysaWu3vf3R8p\nGWVkHUjaSL2uS7qbKh9KuijportveGTEAQBNU9t0U9jZZneOQyU71H7RumameZdJ6pdpX6atHUlP\n0jt+M9t294N0PXcfxrYBAE1V50iik7PsvZK5/zJ1qygr274z7j7MBMRONiBGy83slpntm9nmuG0C\nQFPUeeB6dU51qygrUj5RGFXkhU0/FSTPzezIzLYYXQBoujpHEh9yll2aom4VZUXKi9hTztRUzjGI\noaTbeRsII47B6HVyclKyCQAwP3WGxFD537LzDuKOq1tFWdn2xez456fMds3sY842N/I24O4H7t4b\nvdbWJp+hBgBVqS0kPLmuIDul05V0WKZuFWVl25fHzGLHPSQpe01ER8mpsADQaHWfAtvPHLTthp2z\nzGwzUxatW1FZkfJxukpGI+eEqaaz8AjHLbp5B7cBoGnqvuL6rqS98K37ms5fS3BHyc50t0DdKsrG\nlofwuKPkorpVM/sxXDeRNoj0+yBcyS39dh3FXL14806PX77VL8NTXe6s6N6NK7p5dX3eHwPgC2PJ\nhXdoql6v54NBLHsSL968095PP+v0109ny1a+vqCH339HUADIZWav3L03qR43+GuBxy/fngsISTr9\n9ZMev3y7oBYBaAtCogV+GZ6WWg4ARRESLXC5s1JqOQAURUi0wL0bV7Ty9YVzy1a+vqB7N64sqEUA\n2qL250lg/kYHpzm7CcC8ERItcfPqOqEAYO6YbgIARBESAIAoQgIAEEVIAACiCAkAQBQhAQCIIiQA\nAFGEBAAgipAAAEQREgCAKEICABBFSAAAoggJAEAUIQEAiCIkAABRhAQAIIqQAABEERIAgChCAgAQ\nVfszrs3svqRjSV1JfXd/PU3dKsoKlu9LOnT3fmrZTvjxqaRVSbvu/mDafgNAU9QaEmb2TNLD0Q7S\nzA4lXS9bt4qyAutuS9qUdEvSYaa5HUn7kp4oCYJzfSrTbwBokrqnm7Yz36CPw863bN0qysaWu3vf\n3R8pCYGsoaSLki66+4a7Z+uU6TcANEZtIRF2itmd51A536jH1a2irGz78rj70N2HZfpSZLsAsEh1\njiQ6OcveK5mjL1O3irKy7fuMme2Y2S0z2zezzVTRTNsFgEWq85jE6pzqVlFWpHycfmqK6bmZHZnZ\nVhhZzLJdAFioOkcSH3KWXZqibhVlRcqjco5BDCXdnma7YUQyGL1OTk6KNAEAKlFnSAyVP/USOxAc\nq1tFWdn2nTGzrpl9zFlnY5rtuvuBu/dGr7W1tXEfDwCVqi0kwnUF2amXrj4/nXRs3SrKyrYvx4PM\n+46kozlsFwAWqu5TYPuZg7rd0UVpZraZKYvWraisSPlnwlTT2UjBzDphvYNZtgsATVD3Fdd3Je2Z\nWVfStfB+5I6Sne1ugbpVlI0tDzv5O5K2Ja2a2Y/huglJOghXVEvJNFP29NZJnwsAjWTuvug2YIxe\nr+eDwWDRzQDQMmb2yt17k+pxgz8AQBQhAQCIIiQAAFGEBAAgipAAAEQREgCAKEICABBFSAAAoggJ\nAEAUIQEAiCIkAABRhAQAIIqQAABEERIAgChCAgAQRUgAAKIICQBAFCEBAIgiJAAAUYQEACCKkAAA\nRBESAIAoQgIAEEVIAACiCAkAQNRXdX+gmd2XdCypK6nv7q+nqVtFWcHyfUmH7t5PLetI2glvr0l6\nmPnMUdlTSauSdt39QfQfCQAaotaQMLNnSu1AzexQ0vWydasoK7DutqRNSbckHWaau+/uu6FeV9Ir\nM9ty9+NQ3pG0L+mJkgDK7TMANE3d003bmW/mx2HnW7ZuFWVjy9297+6PlOzkz4RQOBq9D8FwrCRM\nRoaSLkq66O4bqfAAgEarLSTCzja7cxwq51v1uLpVlJVtX8ZolJB1Kf3G3YfuPpywLQBolDpHEp2c\nZe+VzP2XqVtFWdn2nQkjj63M4k1lpqTMbMfMbpnZvpltjtsmADRFncckVudUt4qyIuVROQep++kD\n2+H9aJTy3MyOwjELRhYAGq3OkcSHnGWXcpZNqltFWZHyicJZTj+4+7kpqpxjEENJt8tsGwAWoc6Q\nGCp/SifvIO64ulWUlW1fzL6kH9ILzKxrZh9ztrmRt4EwLTUYvU5OTkp8PADMV20hEaZfslM6XX1+\nOunYulWUlW1fnnB9xf5oCilz3CF7TURHqTOi0tz9wN17o9fa2lqRjweAStR9Cmw/s/PsjubuzWwz\nUxatW1FZkfJcZnZL0mtJH8ysE7bRk86mmjqpup2w3YNJ2wWARav7iuu7kvbCtQXXwvuRO0p2prsF\n6lZRNrY87PjvSNqWtGpmP7r7o1D3WU5f08clDsJIQ0qmmbiYDsBSMHdfdBswRq/X88FgsOhmAGgZ\nM3vl7r1J9bjBHwAgipAAAEQREgCAKEICABBFSAAAoggJAEAUIQEAiCIkAABRhAQAIIqQAABEERIA\ngChCAgAQRUgAAKIICQBAFCEBAIgiJAAAUYQEACCKkAAARBESAIAoQgIAEEVIAACiCAkAQBQhAQCI\nIiQAAFGEBAAg6qu6P9DM7ks6ltSV1Hf319PUraKsYPm+pEN375dcr3C/AaAx3L22l6RnkjZT7w+n\nqVtFWYF1tyXdl3QkabtMv8r0O/va2tpyAJg3SQMvsA+qe7pp289/gz42s+0p6lZRNrbc3fvu/kjJ\naKBsv8r0GwAao7aQCDvF7A52KOl6mbpVlJVtX9G2zrJdAGiCOkcSnZxl75XM0ZepW0VZ2fYVbess\n2wWAhaszJFbnVLeKsiLldW8XABauzpD4kLPs0hR1qygrUh5T1XYBYOHqPAV2qPypl7wDwePqVlFW\ntn1F21p6u2a2I2ln9P6bb76Z8PEAUJ3aRhKeXFeQnXrpSjosU7eKsrLtK9rWabbr7gfu3hu91tbW\nxn08AFSq7lNg+2a2mXrfDTtRmdlmpixat6KyIuWl+zXjdgFgoeq+4vqupD0z60q6Ft6P3FEyLbNb\noG4VZWPLw07+jpKL6lbN7Mdw3cRM2wWAJrPkwjs0Va/X88FgsOhmAGgZM3vl7r1J9caOJMzspZID\nrEeSjt39pzm1DwCwBCZNN11XMn/+12yBmf2TpD9Ieuru/6+CtgEAFmxSSDzPCwhJcvc/S5KZ3TOz\nnqT/4e7/a87tAwAs0KSQeD/6wcyuSvqYDQ13fxzuT9SXdGHuLQQALMykU2Bt9IO7v5H0g5m9N7P/\nbmbfp8r6kt5U1EYAwIJMGkmcO/UpjBquu/s/59TlFBwAaJlJI4memf19ZlnsiWof59AeAECDTBpJ\nbCl5QM5QyTGHQ0m/i9S1yHIAwJKaNJI4kPQPkvaUhMAjSbvhuMSPZvZfzew/h7pclQcALTMpJJ64\n+1/CTeduu/uqpItKbp3xUdIfJb0xs09Knv8MAGiRsSERzmjKLvubuz939z+6+z8oCY07mnxLbQDA\nkpn5LrCj0JD05zm0BwDQIPO8VfiDOW4LANAA8wwJDlwDQMvMHBJm9ndm9idxnQQAtM7UIZEKh78o\n/xnOAIAlN1VImNk9SX9VcmFd193/OM9GAQCaoVRIhNuC/x9Jq5L+4O7/7O5/q6ZpAIBFKxQSZnbX\nzP6vknDoufse4QAA7Tfp3k0ys6eS/ne4cA4A8AWZGBLufruOhgAAmqfsMYnfmdnfVdUYAECzlAqJ\ncBzCzOxbwgIA2q/0KbDhXk1/VRIW/2Rm3867UQCAZpj6YroQFn9WEhb/RTx0CABaZ+KB60nc/S+S\n/mJmG3NoDwCgQWYOiZEQFhOZ2X0lz57oSuq7e+yZ2WPrVlE27bpm9kTSvrvnPlPDzHbCj0+VXGuy\n6+7cNRdA87l7bS9JzyRtpt4fTlO3irIZt/tRyV1w06+PqfL7qeVHSm5lUujfbGtrywFg3iQNvMA+\naJ63Ci9i289/cz82s+0p6lZRNsu6B5I2Uq/rku6m6g6VPMHvortveGTEsQxevHmnf/zTv+kP/+1/\n6h//9G968ebdopsEoEJzm26aJOxQszvHoZIdar9oXTPTvMsk9Wf4zIGSZ4GflZvZtrsfpCu7+1BL\n7sWbd9r76Wed/vpJkvRueKq9n36WJN28ur7IpgGoSJ0jibzbib9XMr9fpm4VZVN/prsPMwGxkw2I\n0XIzu2Vm+2a2mbO9xnv88u1ZQIyc/vpJj1++XVCLAFSttpGEkgO286hbRdms60qSzKyj/EDpp4Lk\nuZkdmdnWso0ufhmelloOYPnVOZL4kLPs0hR1qyibdd2RPWWmziQp5xjEUFLuPbHCiGMwep2cnORV\nW4jLnZVSywEsvzpDYqj8b9l5B3HH1a2ibJbPTNvxz0+p7ZpZ9tGux0oOcH/G3Q/cvTd6ra2t5VVb\niHs3rmjl6wvnlq18fUH3blxZUIsAVK22kHD3vj6ftulKOixTt4qyWT5z9MbMYscuJCl7TURHyamw\nS+Xm1XU9/P47rXdWZJLWOyt6+P13HLQGWqzOYxJSchbRZurbdjfsgDU6mJsqi9atqGzWdbtKRhzn\nuPtxOFah0M9OWPezg9vL4ObVdUIB+ILUHRJ3Je2Fb93XdP5agjtKvmHvFqhbRdms60rSINLvg3C1\ntvTbdRQA0HiWXHiHpur1ej4YxLIHAKZjZq/cvTepXt1XXAMAlgghAQCIIiQAAFGEBAAgipAAAEQR\nEgCAKEICABBFSAAAoggJAEAUIQEAiCIkAABRhAQAIIqQAABEERIAgChCAgAQRUgAAKIICQBAFCEB\nAIgiJAAAUYQEACCKkAAARH216AYAi/DizTs9fvlWvwxPdbmzons3rujm1fVFNwtoHEICX5wXb95p\n76efdfrrJ0nSu+Gp9n76WZIICiCD6SZ8cR6/fHsWECOnv37S45dvF9QioLkICXxxfhmelloOfMlq\nn24ys/uSjiV1JfXd/fU0dasom3ZdM9sJVZ5KWpW06+4Ppu03qnW5s6J3OYFwubOygNYAzVZrSJjZ\nM0kPUzvXQ0nXy9atomzGdTuS9iU9URIE5/pUpt+o3r0bV84dk5Ckla8v6N6NKwtsFdBMdU83bWe+\nQR+b2fYUdasom2XdoaSLki66+4a7H5foC2p28+q6Hn7/ndY7KzJJ650VPfz+Ow5aAzlqG0mEnWJ2\n5zlU8o26X7SumWneZZL6037mqO3uPlSOMv1GfW5eXScUgALqHEl0cpa9VzJHX6ZuFWWzfKak5LiE\nmd0ys30z20zVK9NvAGiUOo9JrM6pbhVls67bT00xPTezIzPbCqOLMv0GgEapcyTxIWfZpSnqVlE2\n07o5xyCGkm4X/NxzwohkMHqdnJzEqgJA5eoMiaHyp16yO9hJdasom/ozzaxrZh9z1tkosN3PuPuB\nu/dGr7W1tbxqAFCL2kLC3fv6fOqlK+mwTN0qymb5zPDzg0xZR9LRpO0KABqu7lNg+5mDut2wE5WZ\nbWbKonUrKptq3TDVdDZSMLNOKDsouF0AaKy6r7i+K2nPzLqSroX3I3eU7Gx3C9StomyWdQ/CFdVS\nMs2UvVBu0ucCQCOZuy+6DRij1+v5YDBYdDMAtIyZvXL33qR63OAPABBFSAAAoggJAEAUIQEAiCIk\nAABRhAQAIIqQAABEERIAgChCAgAQRUgAAKIICQBAFCEBAIgiJAAAUYQEACCKkAAARBESAIAoQgIA\nEEVIAACiCAkAQBQhAQCIIiQAAFGEBAAgipAAAEQREgCAKEICABD1Vd0faGb3JR1L6krqu/vraepW\nUTbtumbWkbQTql2T9DCz3qjsqaRVSbvu/mDCPxVa4MWbd3r88q1+GZ7qcmdF925c0c2r64tuFlBY\nrSFhZs+U2oGa2aGk62XrVlE247r77r4blnclvTKzLXc/DuUdSfuSnigJmdw+o11evHmnvZ9+1umv\nnyRJ74an2vvpZ0kiKLA06p5u2s58cz82s+0p6lZRNtW6IRSORgtDMBxLupWqO5R0UdJFd99IhQda\n7PHLt2cBMXL66yc9fvl2QS0CyqttJBF2ttmd41DJt+p+0bpmpnmXSepP+5mSflQySniUKb+UfuPu\nQ+GL8svwtNTypmPq7MtU50iik7PsvZL5/TJ1qyib+jPD6GIrU7Yp6TC9wMx2zOyWme2b2WbO9tAy\nlzsrpZY32Wjq7N3wVK7fps5evHm36KahYnWGxOqc6lZRNtO6OQep++6eHh313f3A3Z+HA9bPwsFu\ntNi9G1e08vWFc8tWvr6gezeuLKhF02Pq7MtVZ0h8yFl2KWfZpLpVlM26rqSzs5x+cPdzB6ZzjkEM\nJd3O2eZoxDEYvU5OTvKqYQncvLquh99/p/XOikzSemdFD7//bimnaNo2dYbi6jy7aaj8aZu8g7jj\n6lZRNstnpu1L+iG9YHS2k7tfzKy3kbM9ufuBpIPR+16v53n1sBxuXl1fylDIutxZ0bucQFjGqTOU\nU9tIIky/ZKdtusrM3U+qW0XZLJ85ehOuodgfHaDOHHfIXhPRUeqMKKDp2jR1hnLqPgW2n9l5dkdz\n92a2mSmL1q2obOp1zeyWpNeSPphZJ9TrSWdTTWejkDAl1Q0jBmAptGnqDOWYe32zGWEHuSfp35Vc\nmfxj6uK0fUmd1EVp4+rOvWzadbPXSaRcT4VI+orsDSUjjkLXSvR6PR8MBkWqAkBhZvbK3XsT69UZ\nEiiPkABQhaIhwQ3+AABRhAQAIIqQAABEERIAgChCAgAQRUgAAKIICQBAVO2PLwUwXzznAVUiJIAl\nxiNSyyNUy2G6CVhiPOehHB6eVB4jCTQO3/SK4zkP5YwLVf7G8jGSQKPwTa+cNj0itQ6EanmEBBqF\n6ZNyeM5DOYRqeYQEGoVveuXwnIdyCNXyOCaBRuExmeW15RGpdRj9O3HMqzhCAo1y78aVc6d0SnzT\nw3wRquUQEmgUvul92dp0Zltb+kJIoHH4pvdlatOFgW3qCweuATRCm85sa1NfCAkAjdCmM9va1BdC\nAkAjtOkahjb1hZAA0AhtuoahTX3hwDWARmjTmW119aWOM6jM3ee6QcxXr9fzwWCw6GYAaJjsGVRS\nMlopesW9mb1y996kekw3AcASqusMqtqnm8zsvqRjSV1JfXd/PU3dKsoWtV0AKKuuM6hqDQkzeybp\n4WgHaWaHkq6XrVtF2aK2CwDTqOs+Z3VPN21nvkEfm9n2FHWrKFvUdgGgtLrOoKptJBF2iseZxUMl\n36j7ReuameZdJqlfxWdO2q4y/QaAouo6g6rO6aZOzrL3kq6VrFtFWVWfOWm7ADC1Ou5zVud00+qc\n6lZRtqjtAkCj1TmS+JCz7NIUdasoW9R2P2NmO5J2Ru+/+eabWFUAqFydI4mh8qdesvP1k+pWUVbV\nZ07a7mfc/cDde6PX2tpaXjUAqEVtIeHufX0+9dKVdFimbhVlVX3mpO0KABqu7lNg+2a2mXrfDTtR\nmdlmpixat6KyRW0XABqr7iuu70raM7OukrN77qbK7iiZltktULeKskVtFwAaixv8NRw3+ANQhaI3\n+CMkGs7MTpQc/P5bWPS7CT//XtJ/zPCR6W2WrZO3PLtsUvvTP8/Sl1n6ESsr25f0skX1pWw/su+X\nuS+L/r9S9d9X+udp+vL37j75zBh359Xwl6SDoj9LGszrs8rWyVueXVZXX2bpx7z6klm2kL6U7Ueb\n+rLo/ytV/33Nsy/jXtwqfDn8a8mf5/VZZevkLc8uq6svs/QjVla2L8v4O8m+X+a+LPr/StV/X0Xa\nMDOmm1rGzAZeYJ5xGdCXZmpLX9rSD6navjCSaJ+D7ILM6bfL5FxfzGwnvJ4sqkEzyPZlO7z2zSzv\nYssm++xvTJLMbL/uhswo+zt5Ymad8De21L8TM+ua2a3wN9adZcOMJFou3IX2ibtvLLotsxjdTdfd\nj8MDnOTujxbcrKmE0N5z9x9C4L1y99wd77IIffoXd99adFumZWYfw48PWvD7eBb+vm5JWp2lP4wk\nlkz45vnZsyjM7H745nA/PXLw5KK93FuALFrJvnQl3Qo/H0tqVOiV6Yu7v3b3H0KVrhp2y/iyf2PB\nqvLvU7YwU/TjrrtfbGJAlOlLCIZjSXL357P2p/bHl2I64Q9kU8mO8jBTtlRPvpumL5k/9OvZ9RZl\nlt9LuJnjM3dvRIhP25ewc2rMxTwz/E66o3WbMkqdsi9dSZ3we9lWcgbUcNo2MJJYEu7eD3+4eTuU\npXry3Sx9Gc2vuvvziptZyCx9CcG31ZTf1Qx9WZ1lJzRv0/bD3R+FkfdwNKW5aFP25ZKkYSh7LWlv\nljYQEktuwpPvlkrBvuy6+64ablxfMvcpe6XfbkXTSJP6ogaNIsaZ0I9bYWQnNXA6M2vC/5Wj8JKS\nKcCZTlwhJJZf7Ml3M53RsCBj+xLmWh+Gnxvx7XuMcX3Z1m+/n44aeswoZVxfupK2w++m2/Dfy7h+\nHEt6GpZ1lYR3k43rS1+/hVxXyWhiaoTE8hv75LvUf977S3BaX7QvYefzL5JehbNQmn5ab7QvYfpg\nNfxuNhSCr8HG9eV5auqv6U9hHNeP15Juj34nTTx4nTGuL8eSjkJfrrn7g1k+iAPXy2/sk+/Cf+BG\nzN8XEO1LmCu+WG9zZjLp9zLaCS3D72bi0xWX5O/si/mdzLMvjCSWX6kn3zUcfWmmtvSlLf2QauwL\nIbHkvEVPvqMvzdSWvrSlH1K9fSEk2qFNT76jL83Ulr60pR9STX3hmMSSCH8Md5ScGbNqZj+mLvhZ\nqiff0Zdmaktf2tIPqRl94d5NAIAoppsAAFGEBAAgipAAAEQREgCAKEICABBFSAAAoggJAEAUIQEs\nsfDA+ydm5uFJZdnyZ6OynEeOAhNxMR2w5MIVtw8k7Ui6mH1KnJntz3q7aHy5GEkAy29TSUgMlQTF\nmfAMkX9fRKPQDoQEsPxGz5g+0OePQu1pxieT4ctGSADt8UTJUwizdwZdxucloCEICWCJhemkD9LZ\nYytfS9pbaKPQKoQEsNy2lTz4fuShpFsLagtaiOdJAMttNX02k7s/NzOZ2Y6SR1ku6wN10BCMJID2\nGR3A5ngEZkZIAEsqfTwi44mS02I36m0R2oiQAJbXbeWc3urur5VMNR3V3iK0DiEBLBkz64RbcDyR\n9CRccZ21L45HYA64LQcAIIqRBAAgipAAAEQREgCAKEICABBFSAAAoggJAEAUIQEAiCIkAABRhAQA\nIOr/A6gjdTlFrUg/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x114d97da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(Ns,numpy.nanmean(numpy.fabs(dE),axis=1)*10.**3,'o')\n",
    "gca().set_xscale('log')\n",
    "#gca().set_yscale('log')\n",
    "ylabel(r'$\\Delta E$')\n",
    "xlabel(r'$N$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the algorithm for integrating the self-gravitating $\\mathrm{sech}^2$ disk scales much better, approximately as $N^{1/4}\\log N$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x112394a58>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAFKCAYAAABch5ftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XHW9//HXl1qgLO00ZS0IZYq/gsrSJEUWvYpMrbJZ\nIGmRHUsSVLxer7ahXter3pJUuFy9V0zKdWG1TbpAESiZckHEhTapiKClNGVpUZAmU6wEWprv749z\nJmSZTGYyZ+acOfN+Ph55NHNmcubzzUw/+c7nfBdjrUVERPJnL78DEBEJOyVaEZE8U6IVEckzJVoR\nkTxTohURyTMlWhGRPFOiFRHJMyVaEZE8U6IVEckzJVoRkTx7l98B5NtBBx1kp0yZ4ncYIhIy7e3t\nr1lrD87ksaFPtFOmTGH9+vV+hyEiIWOMeSHTx6p0ICKSZ0q0IiJ5VlSlA2NMrftthbW2ztdgREQy\nVDQ9WmNMDIhba5uBzcaYBX7HJCKSCV8TrTGmwU2gg48vMMZUuf+Wu4ejQJX7fScwtVBxiojkwpfS\ngZtcy3ESZ9ug+1qARdbaDvd2GzDT7ckmzRz8cyIiQeVLj9ZaG7fWNuL0TAeLJZOsq7N/r9cYE3XP\n0ZrnMEVEPBGoGq2bUAcn3wRODzapThfCRKSYBG3UQSTFse3ADABjTBWwyP0+Zq2NFzA2EQmZVRu2\nsXjNRl5O9DA5Mo75s6Yxe/oRnj9PoHq0QNlwd7i93SVAuzGmG6fGKyIyKqs2bGPhiqfYlujBAtsS\nPSxc8RSrNmzz/LmC1qPtSnFsEjh1XWBiYcMRkbBavGYjPbv3DDjWs3sPi9ds9LxXG7REmyB1+SDV\nRbOU3EkNyYkNHHXUUR6EJSJh83KiJ6vjuQhU6cDttQ4uH0TJYiiXtbbZWluZ/Dr44IwW1xGREjM5\nMi6r47kIVKJ1xftNUgCI6qKXiHht/qxpjBs7ZsCxcWPHMH/WNM+fy68JC+XAXCAGlBljlrrjagFq\ngIXueNkZ7m0REU8l67CFGHVgrLWenzRIKisrrdajFRGvGWParbWVmTw2iKUDTxhjzjPGNO/YscPv\nUESkxIU20VprV1traydMmOB3KCJS4oI2vEska62trUQiETo7O4nFYkSj0SGPaW521iRqb2+nvr6+\n7zGJRIJly5b13Y7FhiwmJ5IzJVopap2dnaxbt46GhgYAqquraWlpGfCYjo4OKisrKS8vJx6PU11d\nTXt7e9/j29raSCQS1NTUKNFKXoS2dCClobW1lalT31mauKOjY8hjOjs7aWpqAqCyspLOTmf+Szwe\n7+vJRiKRIQlaxCuhTbReXwxLJBJUVFSQSCQ8OZ8XGhsbaW1tpbm5ue+j8UhaW1upr68f9f2FkOx1\nGmOIxwcOoU4m1oqKCuLxONu3b6esbOAcl8GvUVVVVV+PNx6P9/Vak0m5o6ODxsbGIUm6o6OD+vp6\nGhsbyads2ivFKbSJ1uuLYck3+aJFi0Z8bP8eVr4k64xVVVXU1tayefNmWluHX6I3Ho/T2NhIU1NT\nyj8WI91fSLFYjJkzZ1JbW9uXIJOqqqqoq6ujvb29L2F2daVaImOgSMSZ2b106VKWLFnSd7yrq4vy\n8nJqa2uprq4e8DPl5eVMnTqVtrb8rjGfbXul+IQ20XopkUgQiURoaGigsbExbSLq/3E0n5qbm6mq\nquq7PXPmzL6Px6nEYjEWLFhAeXnqRc9Gut8PdXV1xOPxIb/vZNIEmDRp0oD7urq6BtzfX2NjI0uW\nLOm7PxKJDCgdJEsK/RXitUzKpL1SnJRoM5D8uBmLxSgvL0/bq21paWHmzJnD3p/U2NjIzJkzqa+v\nT/kfHFLXG4c7XlZWFoiPlq2trcP+Icq0vJFUXl4+5PedSCQGlAqqqqr6LmwlEokBvb7+cbS2tlJb\nW0skEun7Pc2ZM6fvd59IJAqaVFPJpL1SnJRoR5DszSYN16tNfvRubm5m+/btaet69fX1bN++nYaG\nBmbMmEF1dfWQumhHRwfDzWjr6uoa8p8vGWO+PvY3NjYSj8f7arjD/RGIxWLU1NQMiaOuri7jhJEc\nJZD8uf4Jev369QOSaTQa7atfNjc3D/joXVFR0Xe+mpoaKioqmDhxYt9jIpEIc+fO7atxZ3IxLJFI\n9NXGW1tbh7zOnZ2d1NfX9/2ekudO90cmm/ZKkbLWhvqroqLC5qKpqWnIsWg0ahcsWDDkeHd3t41E\nIiOeM9XPNjU12Wg0aquqqmwsFrPl5eXD/nxLS8uQ5+nu7raA3bx584jPXVtbm9X9TU1NtqWlpe92\nW1ubbW9vH/Yc3d3dNhaL2e7ubmuttbW1tQN+fiSDf+eAbWtrS3lfvrW1tdlYLNZ3u7y8vK9dyfv7\n/77Ky8v7fjeDf3Y4QWpvkKzs2GpPX7TWTqm/z56+aK1d2bHV75AGANbbDPNQaHu0Xow6GO7jZLpe\nbbJnks7ChQuHHEte0Kqrq6O+vr7v43AqqWp2yQtC+fiYGY1G+3pnyUkB6Wq5yaFS1dXV1NXVMXPm\nzAH15Gylukjkh2TJof/vPxaL0dzc3Pde6Ojo6HvPRKPRYT+VpBOU9vqpkLsfFEJoE631YNTBsmXL\nUn5sq6qqIhqNDqnVtrW1ZVSfBWeg/MSJE6moqBgwWiBZCwaGHUVQVlY2JMknb+fjwkksFqOpqYm2\ntra+oUYjlSiSF5q8+OibvEjU2dnpa72yo6Mj5fNHIpG+hFpVVdWXkOPxOHPmzMn6eYLSXj+l2/2g\nGIU20eZqpIsQqXq1/cdoDneBC5wabV1dHd3d3bS0tLBu3ToqKioG1D3T1fTKy8uHJNSurq681fKS\n7WppacFa29eLS6euro7q6mrWrl1LdXV1xrXj/vXKpPLycqLRKNXV1b7WK6PRaMqhZP0/+cydO5eu\nri5aW1spKytLOxIEgt1ePxVy94NCUKIdRnNzM7FYjEQikfIrOae+/0Ws5JjMRCKRNtFWVFT0/QeK\nRqM0NDTQ0tJCfX09U6dOZeLEiWzevDntx+3a2toBPd62tjbq6t7Zhb2zszPtuNpsdHR0DBjRMHfu\n3LSPTybZWCw2oIyQSbKNx+MpyxL19fVDLkwWWlVV1ZDXtrW1te8TDsC6deuora2lqqoqo3JJkNvr\np0LuflAQmRZzi/VrNBfDkheWMv1KXhxZsGCBbWlpyerCTy4aGhpsS0uLbWhoGHLRpKmpacCFmPb2\ndtvQ0GCj0aiNRCK2oaFhwAWtdPc3NTX1XRBLPt9wmpqa+i7k9Nfd3W2rqqqG/bnu7m5bW1trgZQX\n67q7u1NeRMyn9vZ2G4vFLNDX5mQcLS0ttqmpacjvoqmpyUYiERuNRm00GrWxWCzl+yGI7Q2SlR1b\n7XFffcAeXX9f39dxX30gUBfEyOJimBb+FvFIR0cHS5cuZeHChUQikb7eb01NDUuWLAnUZJBisGrD\ntoLsfjBa2Sz8rdW7RDwSj8eZMWPGgJln5eXlGY1EkaFmTz8iUIk1F6Ht0RpjzgPOO/bYY2s2bdrk\ndzhSIpJDvfpPIEmuSSHhkk2PNrSJNkmlAxHJB+0ZJiISIEq0IiJ5pkQrIrJ9c15Pr0QrIqXrrZ2w\n+gvw35WwNX/XcjS8S0RK04u/g5V10P08nP55OOyEvD2VEq2IlJa3d8GjN8Cv/hPGHwlX/QKmnJHX\np1SiFZHS8eqfYEUt/PUPcPJl8PFFsO/4vD9taBNtvwkLfociIn7r7YXf3QLxb8E+B8DcO+H4cwv2\n9KG9GGY93gVXRIpU4iW47XxY8xWY+lH47G8LmmQhxD1aESlx1sIflsL988H2wnnfh/IrwJiCh6JE\nKyLh80YX3Pcv8Mw98O5T4YIfQdkxvoWjRCsi4bKpDe75nJNsz/oGnPEF2GuMryEp0YpIOOz6Bzz0\nVVj/YzjkvXBpKxx+ot9RAUq0IhIGL62DlbXQtQVOuw4++jUYu6/fUfVRohWR4rVnNzzaAI/dCOOP\ngKvugykf9DuqIZRoRaQ4vfpnpxf7lyfhpEvgEzfAvsEczqlEKyLFpbcXnmiCtm+4kw/ugOPP8zuq\ntEKbaDUzTCSEdmyFVZ+FLY/Ce2bB+T+AAw/1O6oRaWaYiASftfCHZfDD053lDM/7L7hkaVEkWQhx\nj1ZEQuKNLrjvi/DMKnj3B9zJB1G/o8qKEq2IZG3Vhm0sXrORlxM9TI6MY/6safnZGvy5OKz6HLzx\nGpz1dTjjX3yffDAaSrQikpVVG7axcMVT9OzeA8C2RA8LVzwF4F2y3d0D8W/C734EBx8Hly6Dw0/y\n5tw+CG2NVkTyY/GajX1JNqln9x4Wr9nozRP89Y/QfKaTZD9wLdQ+UtRJFtSjFZEsvZzoyep4xvrW\njP0mjJsIly2HY2O5nTMglGhFJCuTI+PYliKpTo6MG/1JX/8LrLoWOh+BaWc7w7b2P2j05wsYlQ5E\nJCvzZ01j3NiBF6TGjR3D/FnTRnfCZ+6FW06Dl56Ac2+Gi+8KVZIF9WhFJEvJC145jzp4ayc8WA8b\n7oDJ0+HCW+GgcE4wUqIVkazNnn5EbiMMtq6H5dc4W31/6EvwkYUwZqxn8QWNEq2IFM6et+FXN8Ej\nN8D4yQXZ6jsIlGhFpDC6n3e2+n7pd3BCNZz9PRgX8TuqgghtotWiMiIBYS08+XNnk0RjnFrsidV+\nR1VQoR11oEVlRAKgpxtar3aGbh12Alz7q5JLshDiHq2I+GzLL2HltbDzFWdrmQ9+sSjXKfCCEq2I\neOvtXfB/34HHv++ssjWvDY4o9zsqXynRioh3/rbRGbb11z9AxVUw6z9g7/39jsp3SrQikjtrYd2t\n8NDXYOw4Z3bXcef4HVVgKNGKSG52/g3u+RxsWgNTz4LZP4QDD/M7qkBRohWR0Xv2Ibjns/Dm6/Dx\nBjilFvYK7WCmUVOiFZHs7e5xygTrlsCh74cr7oVD3+t3VIGlRCsi2XnlGVg+D159Bk79nLPFzNh9\n/Y4q0JRoRSQz1sL6H8Oar8A+B4ZqYe58U6IVkZG90QX3fh7+fJ9zweuCH8EBh/gdVdFQohWR9J5/\nHFbUwM5X4WPfccoFuuCVFSVaEUltz9vwy0b45WKYOAWuaXMW6JasKdGKyFCJF2F5Dbz0WzjpEji7\n0anLyqgo0YrIQE+vgtX/7OxKe+ESOHGO3xEVPSVaEXHsegMevB46fgZHVMBFtzqLwkjOlGhFBP76\nFLTOg9eedZYzPPPfQr2HV6Ep0YqUMmvhiSXw0FedbWUuXwlTz/Q7qtAJbaLVVjYiI/jHdmcxmGcf\ngPd8DGbfAvsf5HdUoRTawXDaykYkjc5H4UdnwOa18PEb4JJlSrJ5FNoerYiksGc3PLIIHrsJJh3r\nJNjDT/Q7qtBTohUpFd3PO7sfbF0H0y+HTzRo94MCUaIVKQVPtcJ9X3S+r/oxvP8if+MpMUq0ImH2\n1k54oB5+fwcceYozNnbi0X5HVXKUaEXC6i9PQuunYftm+Kf58OHrYYz+y/tBv3WRsLEWfnsLxL8B\n+02CK1fDMR/yO6qSpkQrEiZvdMGqz8CzD8K0s+H8/4b9J/kdVclTohUJixd/55QK/vEqfKLR2SjR\nGL+jEpRoRYpfby/85gcQ/xZE3g3zHtK6sQGjRCtSzP6xHVZdC5segvfOhvO/D/tqNmTQKNGKFKsX\nfuOUCt54Dc7+Hsy4RqWCgMo60RpjTgbKgIh7KAF0Wmuf9zAuERlOby88fjM8/B2IHAXz2mDyyX5H\nJWlklGiNMRcBMwELGGAzToIFmArMMc5fUgu0WWtXeB+qiPCP12BlHTwXh/ddAOd9H/Yd73dUMoK0\nidYYMx2IAe3W2mszOaEx5ixjzJeBuLX29x7EKCIAL/zaWZz7je1wzk1Q+WmVCorEsInW7cVuttYu\nzuaE1tq1wFpjzHRjzEettQ/nGqRISevthcf/Ex7+rjN99po2OPwkv6OSLAybaK21y3M5sbV2Qy4/\nLyI4pYIVtc66se+7EM77L5UKilBWF8OMMVPcb7usta+7F8YuBp6z1t7qcWwipe35x2H5PGe217n/\nCRVXq1RQpLLdYeFanJptmVu/bQF+DrQbY67xOjiRktTbC79cDD8711kvtmat6rFFLtvhXW1uDRZj\nTAPQnLzgZYwp8zo4kZKz82+wogY6/w/eXwXn3Qz7HOh3VJKjbBPt9n7fx4Cafrdt7uGIlLAtjzk7\nILyZcGqx5VeqFxsS2SbaqcaYTqAOZ8jX6wDGmI96HplIqejdA4/d6OzlVRaFy5bDYe/3OyrxUFaJ\n1lq73Bgz371ZDWCMuQFnlth6IO9DuYwx5dbajnw/j0hB7HzVLRU8AifMgXNvyqlUsGrDNhav2cjL\niR4mR8Yxf9Y0Zk8/wrt4ZVSynoI7eFyttfZ678JJzxgTA5pwZqOJFLctv3RLBTvg/B84GybmUCpY\ntWEbC1c8Rc/uPQBsS/SwcMVTAEq2Pht21IEx5qJ+w7myZow5xhhz4QiPaXCT5+DjC4wxVe6/5cnj\n1to40DnamEQCoXcPPNIAt33SWWmr5mEovyLneuziNRv7kmxSz+49LF6zMafzSu7STlhwp9PGcKbT\nPp/JCd3kXI0zqyzlmgfuOcuBKqBt0H0twKJkecAY04azzoJI8fv7K7DiGqc3e+LFcM6NsM8Bnpz6\n5URPVselcNKWDvoN5ZpvjIniLCSzHadXmVxUJgJEgYPc758badqu2zONG2NSJdCYtba63+1OY0zM\n/RmR4vXCr6HlKnjzdfjk/8DJl3o6qmByZBzbUiTVyZFxnj2HjE5GNdpk4jTGTAAqcRJrsk6aADYk\nk3Iu3J7u4NJAAqdHq0Qrxcla+M3/QNvXYeIUuHwlHPo+z59m/qxpA2q0AOPGjmH+rGmeP5dkJ9tR\nBzuAte5XPkRSHNsOzMjT84nk11t/h3uug2dWwXHnwuxb8rZWQfKCl0YdBE/QdlhIO7vMGFMFRI0x\nC3BmpSXSPV7EV6/+GZZdDtufg5n/Dqf/c94nIMyefoQSawAFLdF2pTjWt1eytbYVaE13AmNMLVCb\nvH3UUUd5FpxIxv64HO75POy9H1xxLxzzIb8jEh8FLdEmSF0+yHhIl7W2GWhO3q6srNTUYCmcPbvh\noa/B726Bd38Aqn8K4yf7HZX4zPNEa4wZn5yamy1rbTzF4jRRnEkKIsH2+l+cUQUv/RY+cC3M/Da8\na2+/o5IAyDnRGmPG4y6d6B6aCczN4ZTxQdNsoxraJYH3/K+g5WrY9Q+46H/hhCq/I5IA8aJH24iz\nWWPSxJF+wJ3tNZd31rZdaq1tdO+uARa643ZnMHCFMJFgsRZ+/QOIf9NZEObKe+GQ4/2OSgLGWJtb\nCdMYc1b/MbTGmGOstVtyjixHxpjzgPOOPfbYmk2bNvkdjoTRm6/DPZ+FP62G4893JiFom5mSYYxp\nt9ZWZvLYbHdYSMUaY042xox3ywgXeXDOnFlrV1traydMmOB3KBJGr/4JlpwJf74fPvYdmHObkqwM\ny4vSQSuwDkgOEDwG+J4H5xUJpqda4d7Pw94HwJWrYcoZfkckAedFoq0eVDqY7sE5RYLn7V3w0Ffh\niSY46jSo+gmMP9zvqKQIeJFop9NvSq62GZdQ2rHNGbq19Qk49XMw81swZqzfUUmR8CLRnuKuO9uZ\n3KhRJFS2/NIZurW7x+nFvj/tMssiQ+ScaK21c8AZT2uMWQSsG24d2kLqN+rA71CkWFkLj98Ma/8d\nJh0Lc++Ag7USlmQv51EHxpiPuj3ah3Gmz25xd2c4OefocqBRB5KTN3fA0suc8bHHn+/sgqAkK6Pk\n1aiD/wDOcpdRBNignXGlaL3ytJNkEy/CrEVw6me07bfkxPNRB9A38qCCAuyKK+KpJ5fC6i84Y2Kv\nvA+OPs3viCQEvKjRDk6y492RBxp9IMXj7V2w5iuwbgkcfYZz0evAQ/2OSkJiVIl2hLJAHbktKiNS\nWDv/BsuugBd/DaddB7FvauiWeGq0PdrrgXac2WBR3lkvNkLq9WQLTqMOJCN/eRLuvgTeeE2rbkne\njDbR1iUXjkmxqMxZnkSWI2vtamB1ZWWlVv+S1P64HFZ9DvYrg08/CJM1qVHyY1SJdtDqXIPHT2k8\nlQRb7x54+Dvwq5vg3afC3NvhgEP8jkpCzItRB5PciQqbcbYg3+7BOUXy480dsLwGNq2B8ivh7O9p\nFwTJOy9GHSwxxlyEs0j3Omvt8tzDEsmD156Dn38KujqdBDvjGo2PlYLIKtEaY6ZYa58ffNxNrkqw\nElyb4tD6adhrDFy+SrvSSkGlnYLbbzHvpLoUjznLGHON55GJeMFaePz7cFc1RN4NtY8oyUrBjdSj\nnQG0GWO6gTgQMcYcba19IfmA5IgDY8w11tpb8xdqdjS8S9jd48zy+sNSeO8nYfYtsPf+fkclJSht\nj9Zau9ZauxfOJopxnItdG4wxe4wxa4wxX05uY0NAxs8maVGZErdjG/zkE06SPfOrUP0zJVnxTUar\nd1lrN1hrlwCt1toy4D04i8kc6/7r+2aMIn1e/B00fwRe2wQX3wUfnq+LXuKrbEcdtAFYaztxZoMt\n8TwikVx03A6/+FcYf4S2/pbAGLZHm2o9g8ELyIgExp7dcP8CuPc6OPp0Z/1YJVkJiHQ92onGmGU4\nExDagrBrgkhKb3RBy5XOljOnfhZmfhvGeDEXR8Qbw74b+4+NdXdMUNKV4Hnlabj7U/D3v8AnfwjT\nL/U7IpEhMvqznybptlhrtbi3+OOZe2HltbDPgXD1A3Bkpd8RiaSU9eerYkm6GkcbYr298GgDPHoD\nHFEBc++E8Yf7HdWIVm3YxuI1G3k50cPkyDjmz5rG7OlH+B2WFICx1uZ+EmMm4Iy1nUvAkm5lZaVd\nv36932GIV976u9OL/fN9cNKn4NybYey+fkc1olUbtrFwxVP07N7Td2zc2DEsuvAEJdsiZYxpt9Zm\n9DEq511wAay1O6y1y92tx68Hphpjlrmreol4o2sL/O/HYOP9zqaJs28piiQLsHjNxgFJFqBn9x4W\nr9noU0RSSJ5fmnV3wl2CxtiKlzofdUYWWAuXLYepxbXJ8suJnqyOS7hk1aM1xkzJTxgiaTyxBG6/\nAA441BkfW2RJFmByZFxWxyVctHqXBFdvL6z5N7j/y3BsDOa1waSpfkc1KvNnTWPc2DEDjo0bO4b5\ns6b5FJEUUmhX75Iit+sNWFkLf1oNp9TCx29w1pItUskLXhp1UJrSJlo3ie5ljJkOVAL1OKt3TcBJ\nvG3uv50EbPUuKWI7X4W7L4ZtHc5Fr1M/E4pFYWZPP0KJtURp9S4Jlr9thFvPgleecTZNPO2zoUiy\nUtq0epcEx5bHYOmlMGZvuOoXcGSF3xGJeCKrUQfpVu8adNHMd8aY84wxzTt27PA7FMnEkz93RxYc\nBtfElWQlVHKesOCOTLgB6PYgHs9oh4UiYS08cgOsrIOjToV5a2DiFL+jEvHUqBNtvwS7BV0Ik9F4\nexes+gw8ssiZTnvZChg30e+oRDw3qkRrjJkPPA9MAKLW2mu9DEpKQE833HEhPHk3fOQrznTad+3t\nd1QieZHVxTA3wdbijDQ4xp1uK5Kd7ufhzmpn7YILmuCki/2OSCSvMkq0xpganDG0LUClEqyM2tZ2\nuHsu7NkFl6+EYz7kd0QieTdionXXm33CWquFXSU3f1oNy2vggEOc4VsHB2f6qdaKlXwaMdG6Sx+K\njJ618NsfOusWHFEBn/o5HHCw31H1GbxW7LZEDwtXPAWgZCueyHb1rglBGy8rAbfnbbh/Pqz5Chx/\nLly5OlBJFrRWrORfVhfDrLU73GQ7Beiy1r6el6gkHN7aCcvnwbMPwmnXObvT7uXJWvOe0lqxkm9Z\nv+vd3RSeB4y7ROIUr4OSEHj9L/CTT8Cmh+Ds78Gs7wYyyYLWipX8G/U73024a3ES7kWAVv4QxytP\nw60x2L7ZqceeUuN3RGlprVjJt5y3srHWbgG2GGOKc0Vm8dbmh2HpFbD3/vDpB+Dwk/yOaERaK1by\nzbM9w9yEK6Ws4za474tw0DS4dBlMONLviDKmtWIln4JZNPOAVu8qoN5eWPvvcO/n4ZgPw6cfLKok\nK5JvoU20Wr2rQHa/CSuugcduhPIr4ZKlsK9GAIr05/l241JC3twBd10ML/4aYt+EM/5FuyGIpKBE\nK6Pzj9echbpf/RNc9L9wQpXfEYkElhKtZG/HNrh9NiRegk/dDe+Z6XdEIoGmRCvZ2b4ZbpvtrCd7\n+Qo4+nS/IxIJPCVaydxf/+iUC+weuGo1TJ7ud0QiRUGJVjLz0jq48yIYuz9ccV/BljjU8oUSBkq0\nMrLOR+DuS5x1ZK+4ByYeXZCn1fKFEhahHUcrHvnzL5xtZyYe7UxEKFCSBS1fKOGhRCvDe3IpLL0c\nDjvB2RHhwMMK+vRavlDCQolWUntiCayshSlnOOWC/coKHoKWL5SwUKKVgax1ptPe/2WYdjZc0gL7\nHOhLKFq+UMJCF8PkHdZC/Bvw+H/BCXNg9g9hzFjfwtHyhRIWSrTi6N0Dv/gStP8EKuc5uyIEYEcE\nLV8oYaBEK7BnN6y8Fv7YCh/8Vzjr61ocRsRDSrSlbncPLLsSNq1xVuD64Bcz/lFNJhDJjBJtKXvz\ndbj7U/DC43DOTTBjXsY/qskEIpnzvwiXJ9phYQRvdMFtn4QXfwMXLskqyYImE4hkI7SJVjsspJHc\nCvyVp+HiO+HE6qxPockEIpkLbaKVYXRtgR/Pgh1b4bJWmPaJUZ1GkwlEMqdEW0pe/bPTk33rdbji\nXjjmn0Z9Kk0mEMmcLoaVim0dcMdFMGZvuOp+OPS9OZ1OkwlEMqdEG1L9h16dPX4zN/fewNgDDoIr\nVkFZ1JPn0GQCkcwo0YZQ/6FXZ+61gRvfupkXOITnTvkpH/coyYpI5lSjDaHk0Kvz9vo1zWNv4ll7\nJNVvfY1v/zLhd2giJUk92hB6OdHD+Xs9zs1jf8g6O415u77MTvYjoaFXIr5QjzaE5h74B24aewu/\n7T2eK3fVs5P9AA29EvGLEm3YbH6Y7+65kaeJUrP7S7zJPoCGXon4SYk2TF74Ddx9CWMOPo6tZ99O\nJFKGAY7LI7X/AAALgUlEQVSIjGPRhSdohICIT1SjDYuXN8Bdc2DCkXD5Ss454GDO+UBuY2VFxBvq\n0YbBq3+C2y+EfSPOONkDDvY7IhHpR4m22HV1wm2znRlfV97j9GhFJFBUOigwTxfL3rEVfvZJ2LML\nrn7AsxlfIuItJdoC8nSx7J2vOuvJvpmAK1fDIcd5Ha6IeESlgwLybLHsnm64/QJ4/WW4tAUmn+xh\nlCLiNfVoC8iTxbLf+jvcUQWvPQuXLIWjTvUoOhHJF/VoCyjnxbJ39zh7fL28Aap/ClM/6l1wIpI3\nSrT9rNqwjTNueJhjrv8FZ9zwMKs2bPP0/Dktlv32Llh2BTz/K7igCY47x9PYRCR/VDpwFWJX11Ev\nlt27B1bUwKaH4NybR7XHl4j4R4nWle5ClZdTV7NeLLu3F+79Z3hmFXzsu1B5tWexiEhhFFWiNcZU\nAQmg3Frb6OW5A7mrq7Xw4PXw+zvgw9fD6df5F4uIjFrR1GiNMeVAmbU2DnQYY2q9PH8gd3V9+Nvw\nRBOcdh185Hr/4hCRnPiaaI0xDcaYWIrjC4wxVe6/5e7hGNDlft8JzPQylsDt6vrYTfDYjVBxFXzs\nO2CMP3GISM58KR24ybUcqALaBt3XAiyy1na4t9twkupUoMN9WBcQ8TKmQO3q+sQSWPstOKEazrlJ\nSVakyPmSaN2P/3FjTKpeacxa2/+yemeKXm9ZPuIKxK6uv78L7v8yTDsHZt8Ce40Z+WdEJNACVaN1\nE2rnoMMJnB7tZgb2YsO30+DTq+Cez0H0I1D1Yxgz1u+IRMQDgUq0pC4HbAeiQNz9F/ffthSPLV7P\nPgTLr4EjT4GL74Kx+/odkYh4JGiJdtiSgFuzTSTru9ba5sKFlWdbHoNll8Oh74VLl8He+/sdkYh4\nKGjjaLtSHJuU/KZfco0XJpwC2NoOd18ME6fAZSth3wl+RyQiHgtaok2QunwwuG47LHd8bd8Y26OO\nOsqDsPLkr3+EOy6E/Q+Cy1fB/pNG/hkRKTqBKh24oxEGlw+yqsdaa5uttZXJr4MPDuj+Wa89B7fP\nhrH7wRX3wvjD/Y5IRPIkUInWFe83SQEg6ibg8Nix1dkdwVq44h6YeLTfEYlIHvk1YaEcmIsz26vM\nGLO039oFNcBCY0wUmOHeDo/db8LSy+DNHfDpB+Dg/+d3RCKSZ35NWOjAmeVVn+K+RL/jrYWMK++s\nhfu/5CzcPfdOOOwEvyMSkQIIYunAE8aY84wxzTt27PA7lHe0/xQ23AEf+jIcf67f0YhIgYQ20Vpr\nV1traydMCMhwqZfWwf3zYepZcOZX/I5GRAootIk2UHa+6mxDM34yXHSr1i8QKTFBG0cbPnt2Q8vV\nzhbh8x6C/fKyHo6IBJgSbb61fQNe+BVc0AyHn+h3NCLig9CWDgJxMeypVvjt/8AHroWT5voXh4j4\nKrSJ1veLYX/9I9xzHRx1urNDgoiUrNAmWl/1dMPSS50FYqp/qnVlRUqcarRe6+2F5TWwYxtcfT8c\neKjfEYmIz5RovfboDfBcG5xzI7z7FL+jEZEAUOnASxsfgEcb4ORLoXKe39GISECENtEWfNTB9s2w\nohYOP8npzWrnWhFxhTbRFnTUwVs74eeXwl7vgrl3wNhx+X9OESkaqtHmylq49zp4bSNctgIiAd7R\nQUR8EdoebcH85r/h6ZVw1jdg6pl+RyMiAaREm4vOR6Ht63D8+XDGF/yORkQCSol2tHZshdarYdJ7\nYPYPdfFLRIalRDsaye1o3t4FF98J+xzod0QiEmChTbR5G97VfzuaC34EB73H2/OLSOiENtHmbXiX\ntqMRkSyFNtHmhbajEZFRUKLNlLajEZFR0oSFTPRtR9MF89q0HY2IZEWJNhPajkZEcqDSwUiS29Gc\nUqftaERkVJRo0+m/Hc2s7/odjYgUKSXa4Wg7GhHxSGgTbU4TFvpvRzP3dm1HIyI5CW2izWnCQnI7\nmk/coO1oRCRnoU20o6btaETEY0q0/b3RBSvqtB2NiHhK42j7268Mzv8+TJ6u7WhExDNKtIO9b7bf\nEYhIyKh0ICKSZ0q0IiJ5pkQrIpJnSrQiInkW2kSbt61sRESyFNpEm7etbEREshTaRCsiEhRKtCIi\neaZEKyKSZ0q0IiJ5Zqy1fseQV8aYvwEvABOA5BCEVN/3P3YQ8Noon7L/ebJ9TKrjg4+lu13MbRnp\n+1zakS7OTO4PUltyeU1S3Vcq76/Btwe3ZTTtONpae3BGj7TWlsQX0Jzu+0HH1nvxPNk+JtXxwcfS\n3S7mtmTw+oy6HZm0Jd39QWpLLq9Jtu+nML2/RmpLru+vkb5KqXSweoTv+x/z6nmyfUyq44OPpbtd\nzG3J5PtcjHSedPcHqS25vCap7iuV99fg2163Ja3Qlw5Gwxiz3lpb6XccXghLW8LSDlBbgijf7Sil\nHm02mgcfMMaU+xGIBwa0xRhT6341+RXQKA1uR8z9ajDGRPwKapSGvL8AjDENhQ7EA4NflyZjTMR9\njxXT6zK4HVFjTJX7HovmenL1aDNgjIkBTdbaqX7Hkgu3HZ3W2k5jzAIAa22jz2Flzf2jt9BaW+3+\nwWi31qZMXsXCbdMSa22F37HkwhjT7X5bX8yviTGmxX1/VQFlubalJHu0bi8oluL4Avev2IL+PVhr\nbRzoLGiQGcqyLVGgyv2+EwjMH45s2mGt7bDWVrsPiQLxQsY6kmzfX64yoKswEWZuFG2psdZODFqS\nzaYdbnLtBLDWtnrSlnxeaQvaFxADFgCbgdig+1qA8n632wbd35bv+ArVFvdYE1BVzO0AaoFav9uQ\na1uAciASpPdYDm1ZkPxZv9sw2na4j29yX5cFQCTXOEqqR2utjVvno3Kq3mnMWtvR73Znqr+AQZFL\nW5I1J2tta57DHFEu7bBOT6MiKK9TDm0ps9Ym8h9h5kbbFmtto3U+ASaS5Sk/jbIdk4CEe18HsDDX\nOEoq0Q4nWbscdDgBzPQhnJxk2JY6a21d4aLKXrp2GGPK+31cbQeKui3A+sJHNTojtKXKGFPrHgtU\naWqwEf6fbHa/wCnn5HwhXInWkerq6Hac+l+xSdsWt/60yP0+ED3BYaRrR4x3XpsIAa2f95OuLVEg\n5r4u0YC/JpC+LZ3AMvdYFOePYFCla0ecd/5IRHF6tTlRonWUpbuz33+CBUUwZGXYtrj/iZcA7e7V\n4SAPWRu2He5HwTL3dZmK+4cjwNK1pbVfCSft+zAg0rWlA5iTfF1swC6IDZKuHZ3AZrcdM6y19bk+\nmbYbd6S62jsp+Y37H8H3emaGhm2LWzubWNhwRm2k1yT5n7gYXpe0bYGieo+F5XUpaDvUo3UkSP1R\nIugfSVMJS1vC0g5QW4KooO1QoqWvpzf4o0QUaPMhnJyEpS1haQeoLUFU6HYo0b4jPmjgddR9MYpR\nWNoSlnaA2hJEBWtHSdVo3V/qXJyr1mXGmKX2nSmoNcBCd4zpDPd2YIWlLWFpB6gtQRSUdmitAxGR\nPFPpQEQkz5RoRUTyTIlWRCTPlGhFRPJMiVZEJM+UaEVE8kyJVkQkz5RopaS5m/A1GWOsMaYlxf0t\nyfuKeINO8ZkmLEjJc2cG1eNsjTNx8G4HxpgGL5bKk9KlHq2Isy5vPc6KTrX973DXH17nR1ASHkq0\nIu/s2dXM0G1xKvFghX0pbUq0Iu9owtlJY/CKTsW21qoEjBKtlDS3NNAFfVuYeLLrqUh/SrRS6mI4\nm/ElLQKqfIpFQqqk1qMVSaGs/ygDa22rMQZ32+xOBiZhkVFRj1ZkqORFMdVnxRNKtFKy+tdnB2nC\nGfI1tbARSVgp0Uopm0OKoVvW2g6cssHmgkckoaREKyXHGBNxp9s2AU3uzLDBGlB9VjyiKbgiInmm\nHq2ISJ4p0YqI5JkSrYhIninRiojkmRKtiEieKdGKiOSZEq2ISJ4p0YqI5JkSrYhInv1/7I9nZMvS\nosAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1120f6780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(Ns,numpy.nanmean(T,axis=1),'o')\n",
    "p= numpy.polyfit(numpy.log(Ns),numpy.log(numpy.nanmean(T,axis=1)/numpy.log(Ns)),deg=1)\n",
    "plot(Ns,numpy.exp(numpy.polyval(p,numpy.log(Ns)))*numpy.log(Ns))\n",
    "pyplot.text(10,10.**1.8,r'$\\Delta t  \\approx %.2f \\mathrm{s} \\times N^{%.2f}\\log N\\,$' % \n",
    "            (numpy.exp(p[1]),p[0]),size=16.)\n",
    "gca().set_xscale('log')\n",
    "gca().set_yscale('log')\n",
    "xlabel(r'$N$')\n",
    "ylabel(r'$\\Delta t / t_{\\mathrm{dyn}}\\,(\\mathrm{s})$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}