[latexpage]
In particle tracking simulations you often need to interpolate particles onto a grid in one or more dimensions. Recently I decided to write a linear particle-to-grid interpolation in one dimension in python.This is an educational introduction into interpolating particles onto a one-dimensional grid.
Interpolation of particles onto a grid
Particles are usually described by a vector of coordinates in an $ n$-dimensional phase space. Often people want to compute a density of particles along one of the coordinate axes. Let us first start with the example of particles in a two-dimensional space, they have coordinates $ (x,x’) $. A common question will be: What is the density of the particles projected on the $ x $ axis?
note: code blocks from here on are executed one after another in an ipython notebook. You can download the notebook here.
To answer the question you first have to think about: How do you compute the density? A builtin way in python/pylab is the histogram function.
%pylab inline
x=10*tan(random.random(1000000))/tan(1)
# Generate 1 Million particles
figure()
hist(x,bins=10) # Calculate the histogram
xlabel("x")
ylabel("number of matches")
show()
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZgAAAELCAYAAADkyZC4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9QFGeeP/D3GCYmGwUU4xBnyKEwiMgIJgFZb71gyOCPb4msqAQ3guJuZU2dMTnL87bq9qJ7p+DmzGoSvdym9IJcasctrhZIThFjjtWtxIlAfriSWjkdhBmQS2BAXAUc/Hz/IHRglYwQenqE96uqq2Yeprvf3SV8fLqffkYnIgIiIqIRNk7rAERENDqxwBARkSpYYIiISBUsMEREpAoWGCIiUgULDBERqWLEC0xubi4MBgMsFsuA9tdffx2zZs1CbGwstm3bprTn5eXBbDYjOjoa5eXlSntVVRUsFgvMZjM2b96stHd1dSEzMxNmsxlJSUm4fPmy8rOCggJERUUhKioKhw8fHulDIyKioZARdurUKamurpbY2Fil7YMPPpCnn35auru7RUTk//7v/0RE5Pz58xIXFyfd3d3icDgkIiJCbt26JSIiCQkJYrfbRURkyZIlcuzYMRER2b9/v2zcuFFERGw2m2RmZoqISEtLi8yYMUPcbre43W7lNRERaWPEezALFizApEmTBrT927/9G372s59Br9cDAB5++GEAQElJCbKysqDX6xEeHo7IyEjY7XY0NTWho6MDiYmJAIDs7GwUFxcDAEpLS5GTkwMAyMjIwMmTJwEAx48fR2pqKoKDgxEcHAyr1YqysrKRPjwiIrpLPrkHU1tbi1OnTiEpKQnJycmorKwEADQ2NsJkMimfM5lMcLlct7UbjUa4XC4AgMvlQlhYGAAgICAAQUFBaGlpGXRbRESkjQBf7MTj8cDtduPMmTM4e/YsVq9ejUuXLvli13ek0+k02zcR0b1KhjizmE96MCaTCStWrAAAJCQkYNy4cfjqq69gNBrR0NCgfM7pdMJkMsFoNMLpdN7WDvT2Zurr6wH0Fq729naEhITctq2GhoYBPZq/JCJcRPDyyy9rnsEfFp4Hnguei29fhsMnBSY9PR0ffPABAODChQvo7u7GlClTkJaWBpvNhu7ubjgcDtTW1iIxMRGhoaEIDAyE3W6HiKCwsBDLly8HAKSlpaGgoAAAUFRUhJSUFABAamoqysvL0dbWBrfbjRMnTmDRokW+ODwiIrqDEb9ElpWVhd///vdoaWlBWFgYfvGLXyA3Nxe5ubmwWCy4//77lSHEMTExWL16NWJiYhAQEIADBw4ol68OHDiAdevW4caNG1i6dCkWL14MANiwYQPWrl0Ls9mMkJAQ2Gw2AMDkyZPx85//HAkJCQCAl19+GcHBwSN9eEREdJd0Mty+zz1Mp9MNu8s32lRUVCA5OVnrGJrjefgGz8U3eC6+MZy/mywwRETk1XD+bvpkFJk/Skv7kab7f+65Z/H//t8STTMQEalpzBaYd99dquHej+LRR0+wwBDRqDZmCwygZQ+mGYDT66eIiO5lnE2ZiIhUwQJDRESqYIEhIiJVsMAQEZEqWGCIiEgVLDBERKQKFhgiIlIFCwwREalizM5FBmh52K9Cr/85bt68rmGGXhMnTsLVq61axyAiP8e5yO4hvcVF+9re0cFv9yQidfASGRERqYIFhoiIVMECQ0REqmCBISIiVbDAEBGRKka8wOTm5sJgMMBisdz2sz179mDcuHFobf1mWGxeXh7MZjOio6NRXl6utFdVVcFiscBsNmPz5s1Ke1dXFzIzM2E2m5GUlITLly8rPysoKEBUVBSioqJw+PDhkT40IiIaChlhp06dkurqaomNjR3QXl9fL4sWLZLw8HBpaWkREZHz589LXFycdHd3i8PhkIiICLl165aIiCQkJIjdbhcRkSVLlsixY8dERGT//v2yceNGERGx2WySmZkpIiItLS0yY8YMcbvd4na7ldd3AkAA0XDZ4wcZ+pYR/ydARKPQcP5WjHgPZsGCBZg0adJt7X/3d3+HX/7ylwPaSkpKkJWVBb1ej/DwcERGRsJut6OpqQkdHR1ITEwEAGRnZ6O4uBgAUFpaipycHABARkYGTp48CQA4fvw4UlNTERwcjODgYFitVpSVlY304RER0V3yyT2YkpISmEwmzJkzZ0B7Y2MjTCaT8t5kMsHlct3WbjQa4XK5AAAulwthYWEAgICAAAQFBaGlpWXQbRERkTZUf5L/+vXr2LVrF06cOKG0iV/MTrO93+vkrxciIgKAiooKVFRUfKdtqF5gLl68iLq6OsTFxQEAnE4nHn/8cdjtdhiNRjQ0NCifdTqdMJlMMBqNcDqdt7UDvb2Z+vp6TJs2DR6PB+3t7QgJCYHRaBxwMhoaGvDUU099S7LtI3mYRESjSnJyMpKTk5X3O3bsGPI2VL9EZrFY0NzcDIfDAYfDAZPJhOrqahgMBqSlpcFms6G7uxsOhwO1tbVITExEaGgoAgMDYbfbISIoLCzE8uXLAQBpaWkoKCgAABQVFSElJQUAkJqaivLycrS1tcHtduPEiRNYtGiR2odHRESDGPEeTFZWFn7/+9+jpaUFYWFh+MUvfoH169crP++dybhXTEwMVq9ejZiYGAQEBODAgQPKzw8cOIB169bhxo0bWLp0KRYvXgwA2LBhA9auXQuz2YyQkBDYbDYAwOTJk/Hzn/8cCQkJAICXX34ZwcHBI314RER0lzhdvyZeBbBF4wx9hj4FNxGNPcOZrp9P8hMRkSr4fTBjXsCAy5Za4JeeEY1OLDBjngdaX6rjl54RjU68REZERKpggSEiIlWwwBARkSpYYIiISBUsMEREpAoWGCIiUgULDBERqYIFhoiIVMECQ0REqmCBISIiVbDAEBGRKlhgiIhIFSwwRESkChYYIiJSBQsMERGpggWGiIhUMeIFJjc3FwaDARaLRWnbunUrZs2ahbi4OKxYsQLt7e3Kz/Ly8mA2mxEdHY3y8nKlvaqqChaLBWazGZs3b1bau7q6kJmZCbPZjKSkJFy+fFn5WUFBAaKiohAVFYXDhw+P9KEREdFQyAg7deqUVFdXS2xsrNJWXl4uPT09IiKybds22bZtm4iInD9/XuLi4qS7u1scDodERETIrVu3REQkISFB7Ha7iIgsWbJEjh07JiIi+/fvl40bN4qIiM1mk8zMTBERaWlpkRkzZojb7Ra32628vhMAAoiGyx4/yNC3+EOOEf9nSEQjbDi/pyPeg1mwYAEmTZo0oM1qtWLcuN5dzZs3D06nEwBQUlKCrKws6PV6hIeHIzIyEna7HU1NTejo6EBiYiIAIDs7G8XFxQCA0tJS5OTkAAAyMjJw8uRJAMDx48eRmpqK4OBgBAcHw2q1oqysbKQPj4iI7lKAr3d46NAhZGVlAQAaGxuRlJSk/MxkMsHlckGv18NkMintRqMRLpcLAOByuRAWFgYACAgIQFBQEFpaWtDY2Dhgnb5tDW57v9fJXy+kjQDodDpNE0ycOAlXr7ZqmoHIn1RUVKCiouI7bcOnBWbnzp24//77sWbNGl/udhDbtQ5ACg8A0TRBR4e2BY7I3yQnJyM5OVl5v2PHjiFvw2ejyN5++20cPXoU77zzjtJmNBrR0NCgvHc6nTCZTDAajcpltP7tfevU19cDADweD9rb2xESEnLbthoaGgb0aIiIyLd8UmDKysrwyiuvoKSkBA888IDSnpaWBpvNhu7ubjgcDtTW1iIxMRGhoaEIDAyE3W6HiKCwsBDLly9X1ikoKAAAFBUVISUlBQCQmpqK8vJytLW1we1248SJE1i0aJEvDo+IiO5kpEcaPPPMM/LII4+IXq8Xk8kkBw8elMjISHn00UclPj5e4uPjlVFgIiI7d+6UiIgImTlzppSVlSntlZWVEhsbKxEREbJp0yalvbOzU1atWiWRkZEyb948cTgcys8OHTokkZGREhkZKW+//fagGaH5yCmOIvPHDEQ0uOH8jui+XnFM6b2hrOVhvwpgi8YZ+mh9Lvwnwxj8VSC6azrd0H9H+CQ/ERGpggWGiIhUMaQC09PTg6tXr6qVhYiIRhGvBSYrKwtXr17Fn//8Z1gsFsyaNQu//OUvfZGNiIjuYV4LTE1NDQIDA1FcXIwlS5agrq4OhYWFvshGRET3MK8FxuPx4ObNmyguLsayZcug1+s1n9aDiIj8n9cC89xzzyE8PBzXrl3D3/zN36Curg5BQUG+yEZERPewIT8HIyLo6elBQIDP58kcMXwOpj+tz4X/ZOBzMESDU+U5mCtXrmDDhg1YvHgxAOCLL75QpmohIiIajNcCs27dOqSmpqKxsREAYDab8atf/Ur1YEREdG/zWmC++uorZGZm4r777gMA6PX6e/ryGBER+YbXAjNhwgS0tLQo78+cOcOb/ERE5JXXrsiePXuwbNkyXLp0CfPnz8eXX36JoqIiX2QjIqJ72F2NIrt58yb+9Kc/AQBmzpwJvV6vejA1cRRZf1qfC//JwFFkRIMbziiyu7qZ8vHHH6Ourg4ejwfV1dUAgOzs7KEnJCKiMcNrgXn22Wdx6dIlxMfHKzf6ARYYIiL6dl4LTFVVFWpqajg9DI1yAX7xb3zixEm4erVV6xhEI8LrKLLY2Fg0NTX5IguRhjzovQ+k7dLR4Vb9SIl8ZdAezLJlywAA165dQ0xMDBITEzF+/HgAvTd7SktLfZOQiIjuSYMWmC1btgC488iBb7uUkJubi//+7//G1KlTce7cOQBAa2srMjMzcfnyZYSHh+O3v/0tgoODAQB5eXk4dOgQ7rvvPrz22mtITU0F0Htpbt26dejs7MTSpUuxb98+AEBXVxeys7NRXV2NkJAQHDlyBH/1V38FACgoKMDOnTsBAP/4j//I+0RERFoSLy5evCjXr19X3l+/fl0uXbo06OdPnTol1dXVEhsbq7Rt3bpVdu/eLSIi+fn5sm3bNhEROX/+vMTFxUl3d7c4HA6JiIiQW7duiYhIQkKC2O12ERFZsmSJHDt2TERE9u/fLxs3bhQREZvNJpmZmSIi0tLSIjNmzBC32y1ut1t5fScABBANlz1+kKFv8YcczNA/B5E/Gs6/Ta/3YFatWjVg9Ni4ceOwevXqQT+/YMECTJo0aUBbaWkpcnJyAAA5OTkoLi4GAJSUlCArKwt6vR7h4eGIjIyE3W5HU1MTOjo6kJiYCKB3xFrfOv23lZGRgZMnTwIAjh8/jtTUVAQHByM4OBhWqxVlZWXeKywREanC6yiynp4e3H///cr78ePHo7u7e0g7aW5uhsFgAAAYDAY0NzcDABobG5GUlKR8zmQyweVyQa/Xw2QyKe1GoxEulwsA4HK5EBYW1hs+IABBQUFoaWlBY2PjgHX6tjW47f1eJ3+9EBERAFRUVKCiouI7bcNrgZkyZQpKSkqwfPlyAL29jilTpgx7hzqdzi+Ggw4sMERE1F9ycjKSk5OV9zt27BjyNrxeInvzzTexa9cuhIWFISwsDPn5+fj3f//3Ie3EYDDgypUrAICmpiZMnToVQG/PpKGhQfmc0+mEyWSC0WiE0+m8rb1vnfr6egC9X+fc3t6OkJCQ27bV0NAwoEdDRES+5bXA3HfffbDb7aipqUFNTQ0++ugjjBvndbUB0tLSlC8pKygoQHp6utJus9nQ3d0Nh8OB2tpaJCYmIjQ0FIGBgbDb7RARFBYWKj2o/tsqKipCSkoKACA1NRXl5eVoa2uD2+3GiRMnsGjRoiHlJCKiEeRtFEB8fPxtbY899tign3/mmWfkkUceEb1eLyaTSQ4dOiQtLS2SkpIiZrNZrFbrgNFdO3fulIiICJk5c6aUlZUp7ZWVlRIbGysRERGyadMmpb2zs1NWrVolkZGRMm/ePHE4HMrPDh06JJGRkRIZGSlvv/32oBmh+YghjiJjhsFzEPmj4fzbHHQ25S+++AI1NTXYunUr/vVf/xUiAp1Oh6tXr+KVV17B+fPnfVYERxpnU+5P63PBDANxVmfyTyM6m/KFCxfw7rvvor29He+++67SPnHiRLz11lvDT0lERGOC1++D+fDDDzF//nxf5fEJ9mD60/pcMMNA7MGQf1Ll+2Dmzp2LN954AzU1Nbhx44YyxPjQoUPDS0lERGOC1+Fga9euRXNzM8rKypCcnIyGhgZMmDDBF9mIiOge5vUSWXx8PD799FPMmTMHn3/+OW7evIkf/OAHsNvtvso44niJrD+tzwUzDMRLZOSfhnOJzGsPpm+amKCgIJw7dw5tbW348ssvh5eQiIjGDK/3YH7yk5+gtbUV//Iv/4K0tDRcu3YN//zP/+yLbEREdA/zeolsNOIlsv60PhfMMBAvkZF/UmUUmdvtxuHDh1FXVwePx6Ps6LXXXhteSiL6FgGaTwY7ceIkXL3aqmkGGh28FpilS5fi+9//PubMmaP8w9f6F4Bo9PJA655URwd/v2lkeC0wXV1dePXVV32RhYiIRhGvo8jWrFmDX//612hqakJra6uyEBERfRuvPZgHHngAW7duxc6dO5Vp+nU6HS5duqR6OCIiund5LTB79uzBxYsXv9O3WBIR0djj9RKZ2WzGgw8+6IssREQ0injtwXzve99DfHw8Fi5ciPHjxwPgMGUiIvLOa4FJT09Henq6MjS574vHiIiIvg2f5NcEn+RnhsH4Qw7OJkC3U2WySyIiouHwaYHJy8vD7NmzYbFYsGbNGnR1daG1tRVWqxVRUVFITU1FW1vbgM+bzWZER0ejvLxcaa+qqoLFYoHZbMbmzZuV9q6uLmRmZsJsNiMpKQmXL1/25eEREVE/gxaYtWvXAgD27t07Ijuqq6vDW2+9herqapw7dw49PT2w2WzIz8+H1WrFhQsXkJKSgvz8fABATU0Njhw5gpqaGpSVleH5559XumcbN27EwYMHUVtbi9raWpSVlQEADh48iJCQENTW1uKll17Ctm3bRiQ7EREN3aAFpqqqCo2NjTh06NCAJ/iH+yR/YGAg9Ho9rl+/Do/Hg+vXr2PatGkoLS1FTk4OACAnJwfFxcUAgJKSEmRlZUGv1yM8PByRkZGw2+1oampCR0cHEhMTAQDZ2dnKOv23lZGRgZMnTw45JxERjYxBR5H99Kc/RUpKCi5duoTHH398wM+G8yT/5MmTsWXLFjz66KN48MEHsWjRIlitVjQ3N8NgMAAADAYDmpubAQCNjY1ISkpS1jeZTHC5XNDr9TCZTEq70WiEy+UCALhcLoSFhfUeWEAAgoKC0NraismTJ98h0fZ+r5O/XoiICAAqKipQUVHxnbYxaIF54YUX8MILL+CnP/0p3nzzze+0EwC4ePEi9u7di7q6OgQFBWHVqlX4z//8zwGf0el0PhwCvd1H+yEiuvckJycjOTlZeb9jx44hb8PrczBvvvkmPvvsM5w6dQo6nQ4LFixAXFzckHdUWVmJ+fPnIyQkBACwYsUKfPTRRwgNDcWVK1cQGhqKpqYmTJ06FUBvz6ShoUFZ3+l0wmQywWg0wul03tbet059fT2mTZsGj8eD9vb2QXovRESkNq+jyPbt24cf/ehH+PLLL9Hc3Ixnn312WE/xR0dH48yZM7hx4wZEBO+//z5iYmKwbNkyFBQUAAAKCgqQnp4OAEhLS4PNZkN3dzccDgdqa2uRmJiI0NBQBAYGwm63Q0RQWFiI5cuXK+v0bauoqAgpKSlDzklERCNEvIiNjZVr164p769duyaxsbHeVruj3bt3S0xMjMTGxkp2drZ0d3dLS0uLpKSkiNlsFqvVKm63W/n8zp07JSIiQmbOnCllZWVKe2VlpcTGxkpERIRs2rRJae/s7JRVq1ZJZGSkzJs3TxwOxx1zABBANFz2+EGGvsUfcjCDf+UI+DqHtsvEiZOG9XeG1HEX5eI2Xp/kt1gs+Pjjj5UJL2/cuIHExEScO3fuOxU2LfFJ/v60PhfMMJA/5PCHDABnFPAvw3mS3+s9mPXr12PevHlYsWIFRATFxcXIzc0ddkgiIhob7mousqqqKvzhD39QbvLPnTvXF9lUwx5Mf1qfC2YYyB9y+EMGgD0Y/zKcHgwnu9QECwwzDMYfcvhDBoAFxr9wsksiIvIbLDBERKSKby0wHo8HCxcu9FUWIiIaRb61wAQEBGDcuHEDptAnIiK6G16HKT/00EOwWCywWq146KGHAPTe7BnO0/xERDR2eC0wK1aswIoVK5RJKEXEhxNSEhHRvequhilfv34d9fX1iI6O9kUm1XGYcn9anwtmGMgfcvhDBoDDlP2LKsOUS0tLMXfuXCxevBgA8MknnyAtLW14CYmIaMzwWmC2b98Ou92OSZMmAQDmzp075C8bIyKiscdrgdHr9QgODh640jg+PkNERN/Oa6WYPXs23nnnHXg8HtTW1mLTpk2YP3++L7IREdE9zGuBef3113H+/HmMHz8eWVlZCAwMxN69e32RjYjGtADla9S1WgID+Y2438VdT3bZ3t7+9QkPVDuT6jiKrD+tzwUzDOQPOfwhA+AfOTiSrY8qo8jOnj0Li8WCOXPmwGKxIC4uDpWVlcMOSUREY4PXBy1zc3Nx4MABLFiwAADwhz/8Abm5ufj8889VD0dERPcurz2YgIAApbgAwA9+8AMEBHitS0RENMYNWmCqqqpQVVWFJ598Es899xwqKipQUVGBjRs34sknnxzWztra2rBy5UrMmjULMTExsNvtaG1thdVqRVRUFFJTUwdMrJmXlwez2Yzo6GiUl5cPyGaxWGA2m7F582alvaurC5mZmTCbzUhKSsLly5eHlZOIiL67QW/yJycn33H+sb7X//M//zPkneXk5ODJJ59Ebm4uPB4P/vznP2Pnzp2YMmUK/v7v/x67d++G2+1Gfn4+ampqsGbNGpw9exYulwtPP/00amtrodPpkJiYiDfeeAOJiYlYunQpXnjhBSxevBgHDhzAH//4Rxw4cABHjhzB7373O9hsttsPmjf5+9H6XDDDQP6Qwx8yAP6Rgzf5+wznJj/ER9ra2mT69Om3tc+cOVOuXLkiIiJNTU0yc+ZMERHZtWuX5OfnK59btGiRfPTRR9LY2CjR0dFK+29+8xt57rnnlM+cOXNGRERu3rwpU6ZMuWMWAAKIhsseP8jQt/hDDmbwrxz+kMFfcmBk/gCOAsM5F15vprjdbhw+fBh1dXXweDxKJRvqdP0OhwMPP/ww1q9fj88++wyPP/449u7di+bmZhgMBgCAwWBAc3MzAKCxsRFJSUnK+iaTCS6XC3q9HiaTSWk3Go1wuVwAAJfLhbCwMAC9946CgoLQ2tqKyZPvNJZ9e7/XyV8vREQEQLkt8l14LTBLly7F97//fcyZMwfjxo2DyPCm6/d4PKiursYbb7yBhIQEvPjii8jPzx/wmb6Hm3xju4/2Q0R070lOTkZycrLyfseOHUPehtcC09XVhVdffXXIG/5LJpMJJpMJCQkJAICVK1ciLy8PoaGhuHLlCkJDQ9HU1ISpU6cC6O2ZNDQ0KOs7nU6YTCYYjUY4nc7b2vvWqa+vx7Rp0+DxeNDe3j5I74WIiNTmdZjymjVr8Otf/xpNTU1obW1VlqEKDQ1FWFgYLly4AAB4//33MXv2bCxbtgwFBQUAgIKCAqSnpwMA0tLSYLPZ0N3dDYfDgdraWiQmJiI0NBSBgYGw2+0QERQWFmL58uXKOn3bKioqQkpKypBzEhHRCPF2k+b111+XwMBAefTRRyU8PFzCw8PveLP+bnz66afyxBNPyJw5c+SHP/yhtLW1SUtLi6SkpIjZbBar1Sput1v5/M6dOyUiIkJmzpwpZWVlSntlZaXExsZKRESEbNq0SWnv7OyUVatWSWRkpMybN08cDscdc0Dzm4e8yc8M/pzDHzL4Sw4M62/daDScc+F1LrLp06fj7NmzmDJlirqVzoc4TLk/rc8FMwzkDzn8IQPgHzn0ADwaZwAmTpyEq1eHfuVoJA1nmLLXezBmsxkPPvjgsEMREd27PNC+yAEdHb4a/DSyvBaY733ve4iPj8fChQsxfvx4AMMbpkxERGOL1wKTnp6u3Hjv47uhxEREdK+66++DGU14D6Y/rc8FMwzkDzn8IQPgHzn8IQPgD1PWqHIPZvr06Xfc0aVLl4a0IyIiGlu8FpizZ88qrzs7O1FUVISWlhZVQxER0b1vWJfIHnvsMVRXV6uRxyd4iaw/rc8FMwzkDzn8IQPgHzn8IQMwai+RVVVVKTf1b926hcrKSvT09AwvIRERjRleC8yWLVuUAhMQEIDw8HD89re/VT0YERHd27wWmO86XTMREY1NXgtMZ2cn/uu//gt1dXXo6emBSO90/f/0T//ki3xERISAe/L5Q68FZvny5QgODsbjjz+OBx54wBeZiIhoAH+YsmboBc5rgXG5XDh+/Piw4hAR0djl9ftg5s+fj88//9wXWYiIaBTx2oM5ffo0/uM//gPTp08fMNkliw4REX0brwXm2LFjvshBRESjjNcCEx4e7oMYREQ02ni9B0NERDQcLDBERKQKnxeYnp4ezJ07F8uWLQMAtLa2wmq1IioqCqmpqWhra1M+m5eXB7PZjOjoaJSXlyvtVVVVsFgsMJvN2Lx5s9Le1dWFzMxMmM1mJCUl4fLly747MCIiGsDnBWbfvn2IiYlRnkrNz8+H1WrFhQsXkJKSgvz8fABATU0Njhw5gpqaGpSVleH5559XZvLcuHEjDh48iNraWtTW1qKsrAwAcPDgQYSEhKC2thYvvfQStm3b5uvDIyKir/m0wDidThw9ehQ//vGPlWJRWlqKnJwcAEBOTg6Ki4sBACUlJcjKyoJer0d4eDgiIyNht9vR1NSEjo4OJCYmAgCys7OVdfpvKyMjAydPnvTl4RERUT9eR5GNpJdeegmvvPIKrl69qrQ1NzfDYDAAAAwGA5qbmwEAjY2NSEpKUj5nMpngcrmg1+thMpmUdqPRCJfLBaB31oGwsDAAvTM/BwUFobW1FZMnT75Dmu39Xid/vRARUa+Kr5fh81mBee+99zB16lTMnTt30BmadTqdDyd02+6j/RAR3YuSMfA/3juGvAWfFZgPP/wQpaWlOHr0KDo7O3H16lWsXbt+YmxRAAAKwUlEQVQWBoMBV65cQWhoKJqamjB16lQAvT2ThoYGZX2n0wmTyQSj0Qin03lbe9869fX1mDZtGjweD9rb2wfpvRARkdp8dg9m165daGhogMPhgM1mw1NPPYXCwkKkpaWhoKAAAFBQUID09HQAQFpaGmw2G7q7u+FwOFBbW4vExESEhoYiMDAQdrsdIoLCwkIsX75cWadvW0VFRUhJSfHV4RER0V/w6T2Y/vouhf3DP/wDVq9ejYMHDw74tsyYmBisXr0aMTExCAgIwIEDB5R1Dhw4gHXr1uHGjRtYunQpFi9eDADYsGED1q5dC7PZjJCQENhsNm0OjoiIoJO+4VxjSG+h0vKwXwWwReMMfbQ+F8wwkD/k8IcMgH/k8IcMgH/k0GGo5YJP8hMRkSpYYIiISBUsMEREpAoWGCIiUgULDBERqYIFhoiIVMECQ0REqmCBISIiVbDAEBGRKlhgiIhIFSwwRESkChYYIiJSBQsMERGpggWGiIhUwQJDRESqYIEhIiJVsMAQEZEqWGCIiEgVPiswDQ0NWLhwIWbPno3Y2Fi89tprAIDW1lZYrVZERUUhNTUVbW1tyjp5eXkwm82Ijo5GeXm50l5VVQWLxQKz2YzNmzcr7V1dXcjMzITZbEZSUhIuX77sq8MjIqK/JD7S1NQkn3zyiYiIdHR0SFRUlNTU1MjWrVtl9+7dIiKSn58v27ZtExGR8+fPS1xcnHR3d4vD4ZCIiAi5deuWiIgkJCSI3W4XEZElS5bIsWPHRERk//79snHjRhERsdlskpmZeccsAAQQDZc9fpChb/GHHMzgXzn8IYO/5PCHDP6SA0P+u++zHkxoaCji4+MBABMmTMCsWbPgcrlQWlqKnJwcAEBOTg6Ki4sBACUlJcjKyoJer0d4eDgiIyNht9vR1NSEjo4OJCYmAgCys7OVdfpvKyMjAydPnvTV4RER0V/Q5B5MXV0dPvnkE8ybNw/Nzc0wGAwAAIPBgObmZgBAY2MjTCaTso7JZILL5bqt3Wg0wuVyAQBcLhfCwsIAAAEBAQgKCkJra6uvDouIiPoJ8PUOr127hoyMDOzbtw8TJ04c8DOdTgedTuejJNv7vU7+eiEiol4VXy/D59MCc/PmTWRkZGDt2rVIT08H0NtruXLlCkJDQ9HU1ISpU6cC6O2ZNDQ0KOs6nU6YTCYYjUY4nc7b2vvWqa+vx7Rp0+DxeNDe3o7JkycPkma7KsdIRDQ6JGPgf7x3DHkLPrtEJiLYsGEDYmJi8OKLLyrtaWlpKCgoAAAUFBQohSctLQ02mw3d3d1wOByora1FYmIiQkNDERgYCLvdDhFBYWEhli9fftu2ioqKkJKS4qvDIyKivzTkYQHDdPr0adHpdBIXFyfx8fESHx8vx44dk5aWFklJSRGz2SxWq1Xcbreyzs6dOyUiIkJmzpwpZWVlSntlZaXExsZKRESEbNq0SWnv7OyUVatWSWRkpMybN08cDscds0DzERkcRcYM/pzDHzL4Sw5/yOAvOTDkv/u6r//gjim993m0POxXAWzROEMfrc8FMwzkDzn8IQPgHzn8IQPgHzl0GGq54JP8RESkChYYIiJSBQsMERGpggWGiIhUwQJDRESqYIEhIiJVsMAQEZEqWGCIiEgVLDBERKQKFhgiIlIFCwwREamCBYaIiFTBAkNERKpggSEiIlWwwBARkSpYYIiISBUsMEREpAoWGCIiUgULDBERqWJUFpiysjJER0fDbDZj9+7dWschIhqTRl2B6enpwd/+7d+irKwMNTU1+M1vfoMvvvhC61hERGPOqCswH3/8MSIjIxEeHg69Xo9nnnkGJSUlWsciIhpzArQOMNJcLhfCwsKU9yaTCXa7/bbPBQUt82WsAbq6LqKzU7PdExH5xKgrMDqd7q4+197+nspJ7sbdZVWfP+Rghm/4Qw5/yAD4Rw5/yAD4T467N+oKjNFoRENDg/K+oaEBJpNpwGdExNexiIjGnFF3D+aJJ55AbW0t6urq0N3djSNHjiAtLU3rWEREY86o68EEBATgjTfewKJFi9DT04MNGzZg1qxZWsciIhpzRl0PBgCWLFmCP/3pT/jf//1f/OxnP1Pa+XxMr4aGBixcuBCzZ89GbGwsXnvtNa0jaa6npwdz587FsmXaDf7wB21tbVi5ciVmzZqFmJgYnDlzRutImsnLy8Ps2bNhsViwZs0adHV1aR3JZ3Jzc2EwGGCxWJS21tZWWK1WREVFITU1FW1tbV63MyoLzJ3w+Zhv6PV6/OpXv8L58+dx5swZ7N+/f8yeiz779u1DTEzMXQ8SGa02b96MpUuX4osvvsDnn38+Znv/dXV1eOutt1BdXY1z586hp6cHNptN61g+s379epSVlQ1oy8/Ph9VqxYULF5CSkoL8/Hyv2xkzBYbPx3wjNDQU8fHxAIAJEyZg1qxZaGxs1DiVdpxOJ44ePYof//jHY3oASHt7O06fPo3c3FwAvZebg4KCNE6ljcDAQOj1ely/fh0ejwfXr1+H0WjUOpbPLFiwAJMmTRrQVlpaipycHABATk4OiouLvW5nzBSYOz0f43K5NEzkH+rq6vDJJ59g3rx5WkfRzEsvvYRXXnkF48aNmV+HO3I4HHj44Yexfv16PPbYY/jJT36C69evax1LE5MnT8aWLVvw6KOPYtq0aQgODsbTTz+tdSxNNTc3w2AwAAAMBgOam5u9rjNmfqPG+qWPO7l27RpWrlyJffv2YcKECVrH0cR7772HqVOnYu7cuWO69wIAHo8H1dXVeP7551FdXY2HHnrori6DjEYXL17E3r17UVdXh8bGRly7dg3vvPOO1rH8hk6nu6u/qWOmwNzN8zFjyc2bN5GRkYFnn30W6enpWsfRzIcffojS0lJMnz4dWVlZ+OCDD5Cdna11LE2YTCaYTCYkJCQAAFauXInq6mqNU2mjsrIS8+fPR0hICAICArBixQp8+OGHWsfSlMFgwJUrVwAATU1NmDp1qtd1xkyB4fMx3xARbNiwATExMXjxxRe1jqOpXbt2oaGhAQ6HAzabDU899RQOHz6sdSxNhIaGIiwsDBcuXAAAvP/++5g9e7bGqbQRHR2NM2fO4MaNGxARvP/++4iJidE6lqbS0tJQUFAAACgoKLi7/5jKGHL06FGJioqSiIgI2bVrl9ZxNHP69GnR6XQSFxcn8fHxEh8fL8eOHdM6luYqKipk2bJlWsfQ1KeffipPPPGEzJkzR374wx9KW1ub1pE0s3v3bomJiZHY2FjJzs6W7u5urSP5zDPPPCOPPPKI6PV6MZlMcujQIWlpaZGUlBQxm81itVrF7XZ73Y5OZIxfeCYiIlWMmUtkRETkWywwRESkChYYIiJSBQsMERGpggWGiIhUwQJDRESqGHXfB0N0r+np6cGRI0dw6dIlhIWF4eOPP8aWLVswY8YMraMRfSfswRBp7LPPPkNGRgZmzJiBW7duYdWqVXjkkUe0jkX0nbHAEGnssccew/jx4/HRRx8hOTkZycnJePDBB7WORfSdscAQaezs2bP46quv8Mc//hHTp0/H6dOntY5ENCJ4D4ZIY2VlZTAYDPjrv/5r/O53v8OUKVO0jkQ0IjgXGRERqYKXyIiISBUsMEREpAoWGCIiUgULDBERqYIFhoiIVMECQ0REqmCBISIiVbDAEBGRKv4/dRbWj41SWk4AAAAASUVORK5CYII=)
The histogram function computes the number of particles between the left and right edges of the bars respectively. We can also look at the data:
numbers,bins=histogram(x,bins=10)
print "there are %d particles with x between\
%0.2f and %0.2f" % (numbers[0], bins[0],bins[1])
there are 154841 particles with x between 0.00 and 0.1
The histogram function does exactly this, count the number of particles between the edges of a bin. Visually the edges are represented by filling a bar between the left and the right edge up to a height proportional to the resulting number of particles.
Another way to look at this would be to say: the histogram function computes the density at the centers of the bins by attributing each particle to its nearest neighbouring bin. We could then compute a value for the density as follows:
numbers,bins=histogram(x,bins=10)
# compute the width of a bin
binwidth=bins[1]-bins[0]
# compute the positions of the centers of the bins
bincenters=(bins[1:]+bins[0:-1])/2
# plot the density (particles per length unit)
plot(bincenters,numbers/binwidth)
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAYsAAAD9CAYAAABN7FvjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4lPW5//F3NIMbQiBCAjNBIJkQYiLogRC1/EwNCaIlLAFCQBOJ24FKkdqW0+0IbVmsp1WOQHtqgwRcQsWW0AohgEZUZBQilYKn5NCBJJNA2yyIZYmB7++PR0Z2yEKeLJ/XdXEZvpnnyf3MJXPn/q4BxhiDiIjIRVxldwAiItLyKVmIiMglKVmIiMglKVmIiMglKVmIiMglKVmIiMglXTRZZGVlERISQmxs7BntL7zwAv379ycmJoZZs2b52+fPn4/b7SYqKoqCggJ/+/bt24mNjcXtdjNjxgx/+/Hjx0lLS8PtdhMfH8/+/fv938vJySEyMpLIyEiWL1/e6AcVEZFGMBexefNmU1RUZGJiYvxtb731lhk2bJipra01xhjz97//3RhjzK5du8yAAQNMbW2t8Xq9Jjw83Jw8edIYY8zgwYONx+MxxhgzYsQIs27dOmOMMYsXLzZTp041xhiTm5tr0tLSjDHGVFZWmr59+5rq6mpTXV3t/1pEROxx0cpi6NChdOnS5Yy2X/3qV3z/+9/H4XAA0K1bNwDy8vJIT0/H4XDQu3dvIiIi8Hg8VFRUcPjwYeLi4gDIyMhg9erVAKxZs4bMzEwAUlNT2bRpEwDr168nOTmZoKAggoKCSEpKIj8/vwlTpIiI1Ee9xyyKi4vZvHkz8fHxJCQksG3bNgDKy8txuVz+17lcLnw+3zntTqcTn88HgM/nIywsDIDAwEA6d+5MZWXlBe8lIiL2CKzvBXV1dVRXV7N161Y++ugjJkyYwN/+9rcrEdtlCQgIsO1ni4i0VqaeOz3Vu7JwuVyMHTsWgMGDB3PVVVfxz3/+E6fTSWlpqf91ZWVluFwunE4nZWVl57SDVWWUlJQAVhI6dOgQwcHB59yrtLT0jErjbMYY/TGGp59+2vYYWsIfvQ96L/ReXPxPQ9Q7WYwePZq33noLgD179lBbW8tNN91ESkoKubm51NbW4vV6KS4uJi4ujtDQUDp16oTH48EYw4oVKxg1ahQAKSkp5OTkALBq1SoSExMBSE5OpqCggJqaGqqrq9mwYQPDhw9v0AOKiEjjXbQbKj09nXfeeYfKykrCwsL4yU9+QlZWFllZWcTGxtKhQwf/tNbo6GgmTJhAdHQ0gYGBLFmyxN9FtGTJEh566CGOHj3Kfffdx7333gvAww8/zIMPPojb7SY4OJjc3FwAunbtyo9//GMGDx4MwNNPP01QUNAVexNEROTiAkxDa5IWIiAgoMFlVVtTWFhIQkKC3WHYTu/DV/RefEXvxVca8rmpZCEi0s405HNT232IiMglKVmIiMgltetkUVcHO3bYHYWISMvXrpPF3/4GKSlw553wyitw/LjdEYmItEztOllERloJ43vfg2XL4Oab4Yc/hC/XCYqIyJfadbIACAyE0aNhwwYoLITPP4fbboMxY2DjRtBEKxERTZ09r88/t7qlFi+G2lqYNg0yM6Fz5yb9MSIittA6iyZmDLz3HixaBAUFkJYG3/wmnHUWlIhIq6J1Fk0sIACGDoWVK2H3bujRA+69F+6+G373O/jiC7sjFBFpHqos6umLL2D1aquLas8eeOwx60/Pns0WgohIo6iyaAYOB4wfbw2GFxTA3/8OMTEwYQK8844GxEWkbVJl0QQ++wyWL7eqjcBAa0D8wQehY0dbwxIROS8NcNvMGHj7bStpvP02TJ5sDYhHRdkdmYjIV9QNZbOAALjnHnjjDfjkEwgKgoQEGDYMfv97a3sREZHWSJXFFVZbayWPRYusleGPPw6PPgohIXZHJiLtlSqLFqhDB0hPh/ffhz/+Efbvt7qlJk+GLVs0IC4irYMqCxtUV1t7US1ZYg2Cf/Ob1oD4NdfYHZmItAca4G5lTp609qR67jkoLYWXXoK4OLujEpG2TsmilTLGWiX+5JOQkQFz5sB119kdlYi0VRqzaKUCAmDiRGsG1b591q63W7bYHZWIyFdUWbRAb7wBTzxhDYz/7Gdw/fV2RyQibYkqizYiNRV27oSDB2HAANi82e6IRKS9U2XRwuXlWduHpKbC/Plwww12RyQirZ0qizZo1Ciryjh0yDpH4+237Y5IRNojVRatyJtvwr//O3zjG/Dzn8ONN9odkYi0Rqos2rj777eqjNpaq8rYsMHuiESkvVBl0UqtX28dupScDP/1XzofXEQuX5NXFllZWYSEhBB7nkOnf/GLX3DVVVdRVVXlb5s/fz5ut5uoqCgKCgr87du3byc2Nha3282MGTP87cePHyctLQ232018fDz79+/3fy8nJ4fIyEgiIyNZvnx5vR6qPRg+3Koyrr7aqjLWrbM7IhFp08xFbN682RQVFZmYmJgz2ktKSszw4cNN7969TWVlpTHGmF27dpkBAwaY2tpa4/V6TXh4uDl58qQxxpjBgwcbj8djjDFmxIgRZt26dcYYYxYvXmymTp1qjDEmNzfXpKWlGWOMqaysNH379jXV1dWmurra//X5XOIR2oWNG43p08eYzExjqqrsjkZEWrqGfG5etLIYOnQoXbp0Oaf929/+Nj//+c/PaMvLyyM9PR2Hw0Hv3r2JiIjA4/FQUVHB4cOHifty06OMjAxWr14NwJo1a8jMzAQgNTWVTZs2AbB+/XqSk5MJCgoiKCiIpKQk8vPzG5sX26zERGv19403Wke8rlljd0Qi0tbUe4A7Ly8Pl8vFrbfeekZ7eXk5LpfL/3eXy4XP5zun3el04vP5APD5fISFhQEQGBhI586dqaysvOC95MI6doQXXoBXX4VvfxseeAAqK+2OSkTaisD6vPjIkSPMmzePDadNwzEtYHB59uzZ/q8TEhJISEiwLRa73X03/PnP8KMfWWMZixbB2LF2RyUidiosLKSwsLBR96hXsti7dy/79u1jwIABAJSVlfFv//ZveDwenE4npaWl/teWlZXhcrlwOp2UlZWd0w5WlVFSUkLPnj2pq6vj0KFDBAcH43Q6z3iw0tJS7rnnngvGdXqyEGuV93PPwbhxkJVl7Wi7aBF062Z3ZCJih7N/iZ4zZ06971GvbqjY2FgOHjyI1+vF6/XicrkoKioiJCSElJQUcnNzqa2txev1UlxcTFxcHKGhoXTq1AmPx4MxhhUrVjBq1CgAUlJSyMnJAWDVqlUkJiYCkJycTEFBATU1NVRXV7NhwwaGDx9e74dr7+66C3bsgF69rCpj5UqdzCciDXSx0e+JEyeaHj16mA4dOhiXy2WWLl16xvf79Onjnw1ljDFz58414eHhpl+/fiY/P9/fvm3bNhMTE2PCw8PN9OnT/e3Hjh0z48ePNxEREWbIkCHG6/X6v7d06VITERFhIiIizLJlyy4Y4yUeQb60dasx/fsbM3asMQcO2B2NiNipIZ+bWpTXjhw7Bj/5CWRnwy9/CZMmWWdpiEj7opPy5LJs2wZTpkCfPvDrX0PPnnZHJCLNSXtDyWUZNAi2b7dO5Bs4EHJyNJYhIhenyqKd27EDHnrIqi7+53/gy2UvItKGqbKQehs4ED76CO64AwYP1ql8InJ+qizEb8MGmDzZ2sU2I8PuaETkStEAtzTa7t3W4UqTJlkzp65S7SnS5ihZSJP4+99h9Ghr/GLZMrjuOrsjEpGmpDELaRLdu8Nbb1lnZXz963DwoN0RiYjdlCzkvK69Fl55BUaMgCFD4C9/sTsiEbGTuqHkkl55BWbOhOXL4d577Y5GRBpL3VByRUyeDH/4g7Xqe/Fiu6MRETuospDL9re/wf33Q3KytbfU1VfbHZGINIRmQ8kVV1MD48fDNdfAa69ZR7mKSOuibii54oKCYO1acDrha1+DkhK7IxKR5qBkIfXmcFi71WZmWtuEfPSR3RGJyJWmbihplLw8eOQRK3mkptodjYhcjoZ8btbrDG6Rs40aZa30HjUKioth1iwdqCTSFqmykCbh88HIkdYutr/+NXToYHdEInIhGuAW2zid1vbmVVXW1NqqKrsjEpGmpGQhTaZjR3jjDetcjPh4q1tKRNoGJQtpUldfDc8+C9/5DgwdCu+8Y3dEItIUlCzkinjsMXj5ZWsBX06O3dGISGNpgFuuqE8/tQ5TmjgRfvpTHaYk0hJouw9pkf7xDxgzBnr0sHau1WFKIvbSbChpkbp1g40brem0CQlw4IDdEYlIfSlZSLO49lprDOP++62ZUjt32h2RiNSHuqGk2b32GsyYYQ18jxhhdzQi7Y+6oaRVSE+H1ashKwsWLbI7GhG5HKosxDZer9UtlZgIzz0HgdqpTKRZNHllkZWVRUhICLGxsf627373u/Tv358BAwYwduxYDh065P/e/PnzcbvdREVFUVBQ4G/fvn07sbGxuN1uZsyY4W8/fvw4aWlpuN1u4uPj2b9/v/97OTk5REZGEhkZyfLly+v1UNI69OkDW7bAX/8KKSnw2Wd2RyQiF2QuYvPmzaaoqMjExMT42woKCsyJEyeMMcbMmjXLzJo1yxhjzK5du8yAAQNMbW2t8Xq9Jjw83Jw8edIYY8zgwYONx+MxxhgzYsQIs27dOmOMMYsXLzZTp041xhiTm5tr0tLSjDHGVFZWmr59+5rq6mpTXV3t//p8LvEI0grU1hrz+OPGxMYas3+/3dGItH0N+dy8aGUxdOhQunTpckZbUlISV325smrIkCGUlZUBkJeXR3p6Og6Hg969exMREYHH46GiooLDhw8TFxcHQEZGBqtXrwZgzZo1ZGZmApCamsqmTZsAWL9+PcnJyQQFBREUFERSUhL5+flNliClZXE44Fe/gilTdJiSSEvVqF7ipUuXkp6eDkB5eTnx8fH+77lcLnw+Hw6HA5fL5W93Op34fD4AfD4fYWFhViCBgXTu3JnKykrKy8vPuObUvS5k9uzZ/q8TEhJISEhozGOJDQICYOZM6NsX7rsPXnrJWvktIo1XWFhIYWFho+7R4GQxd+5cOnTowKRJkxoVQFM4PVlI6zZqFISGwujRMHs2PP643RGJtH5n/xI9Z86cet+jQVNnly1bxtq1a3nllVf8bU6nk9LSUv/fy8rKcLlcOJ1Of1fV6e2nrikpKQGgrq6OQ4cOERwcfM69SktLz6g0pG0bMgTeew9+8Qv4wQ9Ak91E7FfvZJGfn8+zzz5LXl4e1157rb89JSWF3Nxcamtr8Xq9FBcXExcXR2hoKJ06dcLj8WCMYcWKFYwaNcp/Tc6XW5KuWrWKxMREAJKTkykoKKCmpobq6mo2bNjA8OHDm+J5pZUID4f334e334YHH4TaWrsjEmnnLjb6PXHiRNOjRw/jcDiMy+Uy2dnZJiIiwvTq1csMHDjQDBw40D+byRhj5s6da8LDw02/fv1Mfn6+v33btm0mJibGhIeHm+nTp/vbjx07ZsaPH28iIiLMkCFDjNfr9X9v6dKlJiIiwkRERJhly5Y16ai+tB7/+pcxo0cbc889xtTU2B2NSNvQkM9NLcqTFu/ECWvw++23Ye1a+HJOhIg0kLb7kDbp6qth4UJ46CG480745BO7IxJpf1RZSKuyciVMnw6vvgrDhtkdjUjrpMpC2ry0NFi1CiZPtg5SEpHmocpCWqVPP7UW7z38MPzwh9aiPhG5PDpWVdqVigprlfftt1vbhWjXWpHLo2Qh7c7nn8P48VZl8bvfQceOdkck0vJpzELanY4dYc0a6NkT7r5b53uLXClKFtLqORzw4ovWflJ33AH/+792RyTS9qiXV9qEgAD48Y+hVy+rwli1CoYOtTsqkbZDlYW0KZmZ8PLLkJoKr79udzQibYcqC2lzkpJgwwZrplRpqbVViKbWijSOZkNJm1VaCiNGQGIi/PKX1rYhIqKpsyLnqKmBsWMhKAheeQWuu87uiETsp6mzImcJCoJ16+D6660K45//tDsikdZJyULavGuusfaRSkiwdq3du9fuiERaHw1wS7tw1VUwb541tXboUFi9GuLi7I5KpPXQmIW0O3/8o7UBYXY2jBxpdzQizU9jFiKXYeRIePNNePxxawNCEbk0VRbSbu3da02tTU2FuXOtriqR9kBTZ0Xq6Z//hJQU6NMHli61BsNF2jp1Q4nU0003waZNcOwY3HuvtS5DRM6lZCHt3nXXWWdhDBgAX/salJTYHZFIy6NkIYK1Fcjzz1uzpO66C3bssDsikZZFyULkNDNnwnPPWZsRrlxpdzQiLYcGuEXOY8cOa0+psWNhwQKd7y1ti2ZDiTShqiqYNAmOH7eqjO7d7Y5IpGloNpRIE+ra1Vq8d9ddMGgQfPih3RGJ2EeVhchlWL0aHnvM2l/qkUfsjkakcZq8ssjKyiIkJITY2Fh/W1VVFUlJSURGRpKcnEzNaRPT58+fj9vtJioqioKCAn/79u3biY2Nxe12M2PGDH/78ePHSUtLw+12Ex8fz/79+/3fy8nJITIyksjISJYvX16vhxJpaqNHw7vvWocoPfaY1TUl0q6Yi9i8ebMpKioyMTEx/rbvfve75plnnjHGGLNgwQIza9YsY4wxu3btMgMGDDC1tbXG6/Wa8PBwc/LkSWOMMYMHDzYej8cYY8yIESPMunXrjDHGLF682EydOtUYY0xubq5JS0szxhhTWVlp+vbta6qrq011dbX/6/O5xCOINKnPPjMmNdWYuDhjSkrsjkakYRryuXnRymLo0KF06dLljLY1a9aQmZkJQGZmJqtXrwYgLy+P9PR0HA4HvXv3JiIiAo/HQ0VFBYcPHybuy/2gMzIy/Necfq/U1FQ2bdoEwPr160lOTiYoKIigoCCSkpLIz89vsgQp0lA33givv27tJxUXB4WFdkck0jzqPSHw4MGDhISEABASEsLBgwcBKC8vJz4+3v86l8uFz+fD4XDgcrn87U6nE5/PB4DP5yMsLMwKJDCQzp07U1lZSXl5+RnXnLrXhcyePdv/dUJCAgkJCfV9LJHLFhAA3/se3H47TJxofT1zptUu0hIVFhZS2MjfbBo1ezwgIICAFvAv5PRkIdJchg0Dj8dai/Hhh9b5GDfcYHdUIuc6+5foOXPm1Pse9Z46GxISwoEDBwCoqKig+5eTz51OJ6Wlpf7XlZWV4XK5cDqdlJWVndN+6pqSLzfiqaur49ChQwQHB59zr9LS0jMqDZGW4uab4b33rP2l4uPh//7P7ohErox6J4uUlBRycnIAa8bS6NGj/e25ubnU1tbi9XopLi4mLi6O0NBQOnXqhMfjwRjDihUrGDVq1Dn3WrVqFYmJiQAkJydTUFBATU0N1dXVbNiwgeHDhzfJA4s0teuus7Y3nzbNOuP7zTftjkjkCrjY6PfEiRNNjx49jMPhMC6XyyxdutRUVlaaxMRE43a7TVJS0hmzlObOnWvCw8NNv379TH5+vr9927ZtJiYmxoSHh5vp06f7248dO2bGjx9vIiIizJAhQ4zX6/V/b+nSpSYiIsJERESYZcuWNemovsiV8v77xjidxsyebcyJE3ZHI3J+Dfnc1KI8kSZ24ABMmACdOsHLL0NQkN0RiZxJ232ItAChodaBSuHhMHgw7Nxpd0QijadkIXIFOBywcCE8/TTccw/k5todkUjjqBtK5Ao7td35mDHwzDPa7lzspy3KRVoobXcuLYnGLERaqLO3O/d47I5IpH5UWYg0s7w8ePRRmDvX+q9Ic1M3lEgr8de/WmMYd90FixbBNdfYHZG0J+qGEmkl+vWzuqKqq2HoUDhtdxuRFknJQsQmp7Y7HzfO2u787bftjkjkwtQNJdICbNwIDzyg7c6leWjMQqQV27/fWo/hdsNvfwsdO9odkbRVGrMQacVO3+78jjuguNjuiES+omQh0oKc2u78m9+0Zkq9/rrdEYlY1A0l0kJ9+CFMnmzNllq40BoQF2kK6oYSaUPi4uDjj+Gqq+C222DrVrsjkvZMlYVIK/DGG9ZJfE88Ad//vjYjlMbRbCiRNszng8xMOHYMVqyAPn3sjkhaK3VDibRhTicUFFjbhMTFWafwiTQXVRYirdCOHdaW5wMHwpIlOrpV6keVhUg7MXAgbNsGXbpYX7/7rt0RSVunykKklfvTn6ytzrOyYPZs60hXkYtRZSHSDn3jG1a31McfWwv5tPJbrgQlC5E2ICTEOokvIwPuvBOys0EFtzQldUOJtDG7dlmD3xER8JvfQHCw3RFJS6NuKBHhllusrUJ694YBA6ztz0UaS5WFSBu2YQNMmQITJ1pnfuv4VgFVFiJylqQka/B7714YMgR277Y7ImmtlCxE2ribboLf/97aV+ruu61FfCrGpb4anCzmz5/PLbfcQmxsLJMmTeL48eNUVVWRlJREZGQkycnJ1NTUnPF6t9tNVFQUBQUF/vbt27cTGxuL2+1mxowZ/vbjx4+TlpaG2+0mPj6e/fv3NzRUkXYvIAAeeQTef986L2PkSDh40O6opDVpULLYt28fL774IkVFRezcuZMTJ06Qm5vLggULSEpKYs+ePSQmJrJgwQIAdu/ezcqVK9m9ezf5+flMmzbN3182depUsrOzKS4upri4mPz8fACys7MJDg6muLiYmTNnMmvWrCZ6ZJH2KzIStmyBW2+1tj1fu9buiKS1aFCy6NSpEw6HgyNHjlBXV8eRI0fo2bMna9asITMzE4DMzExWr14NQF5eHunp6TgcDnr37k1ERAQej4eKigoOHz5MXFwcABkZGf5rTr9XamoqmzZtavTDigh06ADz5sFrr8HUqVb31NGjdkclLV2DdsXv2rUrTz31FL169eK6665j+PDhJCUlcfDgQUJCQgAICQnh4Jd1bnl5OfHx8f7rXS4XPp8Ph8OBy+XytzudTnw+HwA+n4+wsDAryMBAOnfuTFVVFV27dj0nntmzZ/u/TkhIICEhoSGPJdKu3H23Nfg9dSoMGgSvvmpNtZW2p7CwkMLCwkbdo0HJYu/evTz//PPs27ePzp07M378eF4+a7/kgIAAAgICGhXc5To9WYjI5evSxaowXn4Zhg2zDlZ68knrdD5pO87+JXrOnDn1vkeD/pfYtm0bd955J8HBwQQGBjJ27Fg++OADQkNDOXDgAAAVFRV0794dsCqG0tJS//VlZWW4XC6cTidlZWXntJ+6pqSkBIC6ujoOHTp03qpCRBonIAAefNBayPfGGzB8OJSX2x2VtDQNShZRUVFs3bqVo0ePYoxh48aNREdHM3LkSHJycgDIyclh9OjRAKSkpJCbm0ttbS1er5fi4mLi4uIIDQ2lU6dOeDwejDGsWLGCUaNG+a85da9Vq1aRmJjYFM8rIhfQpw+88w4MHQq33w5/+IPdEUmLYhromWeeMdHR0SYmJsZkZGSY2tpaU1lZaRITE43b7TZJSUmmurra//q5c+ea8PBw069fP5Ofn+9v37Ztm4mJiTHh4eFm+vTp/vZjx46Z8ePHm4iICDNkyBDj9XrPG0cjHkFELmDLFmP69jXmkUeMOXzY7mikqTXkc1PbfYjIeR0+DNOnw3vvwW9/C5o30nY05HNTyUJELmrNGpg2DVJS4Jln4MYb7Y5IGkt7Q4lIk0tJgb/8BY4fh5gYWL/e7ojEDqosROSyFRTAY4/B178Ov/ylNfVWWh9VFiJyRSUnw86dcMMNVpWRl2d3RNJcVFmISINs3mxtTnj77fDCC9Ctm90RyeVSZSEizeb//T9ru5CwMIiNhdxcbX3elqmyEJFG+/BDyMqyzv1esgR69rQ7IrkYVRYiYou4ONi+3dr6fOBAeOklVRltjSoLEWlSf/6zde53t27wm9/AzTfbHZGcTZWFiNhuwADweKwV34MGwa9+BSdP2h2VNJYqCxG5Yj791BrL6NABsrOtMQ2xnyoLEWlR+ve39pYaMwbi4+EXv4ATJ+yOShpClYWINIu9e611GUePwtKlEB1td0TtlyoLEWmxwsNh0yZr8Pvuu2HuXPjiC7ujksulykJEml1JCTz+OBw4YFUZt91md0TtiyoLEWkVevWCtWut876HD4cf/cja1VZaLiULEbFFQABkZlrrMnbvtvaY8njsjkouRN1QImI7Y2DVKvjWt2DSJPjpT+H66+2Oqu1SN5SItEoBATB+vLX9+YED1sK+d96xOyo5nSoLEWlx/vhHmDpVR7leKaosRKRNGDnSOsq1ttba/lxHudpPlYWItGgbNljTbPv1gx//GO680+6IWj9VFiLS5iQlWXtMjRkDkydDYqLGM+ygykJEWo0vvoCXX4Z586BHD6vSGDbMGiCXy9eQz00lCxFpderqYOVK+NnPoHNn+M//hBEjlDQul5KFiLQrJ07AG29YSaNDB2sleEoKXKUO9otSshCRdunkScjLsxbznThhJY3UVCWNC1GyEJF2zRhrz6mf/AQOH7aSRloaXH213ZG1LM06G6qmpoZx48bRv39/oqOj8Xg8VFVVkZSURGRkJMnJydTU1PhfP3/+fNxuN1FRURQUFPjbt2/fTmxsLG63mxkzZvjbjx8/TlpaGm63m/j4ePbv39/QUEWknQgIgPvvh61b4fnnrSNd+/eHZcu0HXpjNThZzJgxg/vuu49PP/2UTz75hKioKBYsWEBSUhJ79uwhMTGRBQsWALB7925WrlzJ7t27yc/PZ9q0af6sNnXqVLKzsykuLqa4uJj8/HwAsrOzCQ4Opri4mJkzZzJr1qwmeFwRaQ8CAiA5GTZvht/8BpYvt9ZpvPiitdBPGsA0QE1NjenTp8857f369TMHDhwwxhhTUVFh+vXrZ4wxZt68eWbBggX+1w0fPtx88MEHpry83ERFRfnbX3vtNfP444/7X7N161ZjjDFffPGFuemmm84bSwMfQUTamffeM2b4cGPCwoxZtMiYo0ftjsg+DfncbFBl4fV66datG1OmTOH222/n0Ucf5V//+hcHDx4kJCQEgJCQEA4ePAhAeXk5LpfLf73L5cLn853T7nQ68fl8APh8PsLCwgAIDAykc+fOVFVVNSRcERHuugvy863dbfPzrZP7nn8nY0eKAAAKVUlEQVQejhyxO7LWIbAhF9XV1VFUVMSiRYsYPHgwTz75pL/L6ZSAgAACmmnS8+zZs/1fJyQkkJCQ0Cw/V0Ran7g4a6PCjz+2ptwuWABPPWVtXNixo93RXRmFhYUUFhY26h4NShYulwuXy8XgwYMBGDduHPPnzyc0NJQDBw4QGhpKRUUF3bt3B6yKobS01H99WVkZLpcLp9NJWVnZOe2nrikpKaFnz57U1dVx6NAhunbtet54Tk8WIiKX47bbrDUaO3da54H37QszZsATT1gL/dqSs3+JnjNnTr3v0aBuqNDQUMLCwtizZw8AGzdu5JZbbmHkyJHk5OQAkJOTw+jRowFISUkhNzeX2tpavF4vxcXFxMXFERoaSqdOnfB4PBhjWLFiBaNGjfJfc+peq1atIjExsSGhiohcVGws5OZag+F//StERMDs2VBdbXdkLUxDB0h27NhhBg0aZG699VYzZswYU1NTYyorK01iYqJxu90mKSnJVFdX+18/d+5cEx4ebvr162fy8/P97du2bTMxMTEmPDzcTJ8+3d9+7NgxM378eBMREWGGDBlivF7veeNoxCOIiJyjuNiYrCxjunY15gc/MOYf/7A7oqbXkM9NLcoTETmPffus8YzXX4esLPjOd+DL+TutnrYoFxFpIr17w69/DX/+Mxw/bi3umz4diovtjsweShYiIhfhcsF//zfs2gWdOllTcEeOhI0bre1F2gt1Q4mI1MPRo/Dqq9YaDWOsGVSTJ8P119sd2eXTRoIiIs3EGHj7bStpbN0KjzwC06ZZlUhLpzELEZFmEhAA99wDa9bAli3wr3/BrbfCxIlW8mhrVFmIiDSRQ4fgpZfghRfgppvgySdh3DhwOOyO7EzqhhIRaQFOnIA//QkWLoQ9e6zuqccesxJIS6BuKBGRFuDqq2HUKHjrLeswpr17we2GRx+Fv/zF7ugaRslCROQKuvVWyM62Koybb7bO2Rg2zNrM8ORJu6O7fOqGEhFpRrW11qrwhQuhqgq+9S2YMgVuvLH5YtCYhYhIK2HMV8e/btwIGRnWCvG+fa/8z9aYhYhIKxEQAHfcAStXwo4dcM01MGQIjB5trd9oab8Dq7IQEWkhjhyBFSusLiqHw1odPmkSXHtt0/4cdUOJiLQBxlhdU88/D9u2WdNup06Fnj2b5v7qhhIRaQMCAiApCd580zqUqboaYmLggQfgo49sikmVhYhIy1dTY03BXbwYCguhV6+G30vdUCIibdzJk3BVI/uE1A0lItLGNTZRNPjn2vNjRUSkNVGyEBGRS1KyEBGRS1KyEBGRS1KyEBGRS1KyEBGRS1KyEBGRS1KyEBGRS1KyEBGRS1KyEBGRS2pUsjhx4gS33XYbI0eOBKCqqoqkpCQiIyNJTk6mpqbG/9r58+fjdruJioqioKDA3759+3ZiY2Nxu93MmDHD3378+HHS0tJwu93Ex8ezf//+xoTaLhQWFtodQoug9+Erei++oveicRqVLBYuXEh0dDQBAQEALFiwgKSkJPbs2UNiYiILFiwAYPfu3axcuZLdu3eTn5/PtGnT/JtYTZ06lezsbIqLiykuLiY/Px+A7OxsgoODKS4uZubMmcyaNasxobYL+sdg0fvwFb0XX9F70TgNThZlZWWsXbuWRx55xP/Bv2bNGjIzMwHIzMxk9erVAOTl5ZGeno7D4aB3795ERETg8XioqKjg8OHDxMXFAZCRkeG/5vR7paamsmnTpoY/pYiINEqDk8XMmTN59tlnueq0LRAPHjxISEgIACEhIRw8eBCA8vJyXC6X/3Uulwufz3dOu9PpxOfzAeDz+QgLCwMgMDCQzp07U1VV1dBwRUSkEQIbctGf/vQnunfvzm233XbB0i4gIMDfPXWlNdfPaQ3mzJljdwgtgt6Hr+i9+Irei4ZrULLYsmULa9asYe3atRw7dozPPvuMBx98kJCQEA4cOEBoaCgVFRV0794dsCqG0tJS//VlZWW4XC6cTidlZWXntJ+6pqSkhJ49e1JXV8ehQ4fo2rXrObHo4CMRkSuvQd1Q8+bNo7S0FK/XS25uLvfccw8rVqwgJSWFnJwcAHJychg9ejQAKSkp5ObmUltbi9frpbi4mLi4OEJDQ+nUqRMejwdjDCtWrGDUqFH+a07da9WqVSQmJjbF84qISAM0qLI426luoP/4j/9gwoQJZGdn07t3b373u98BEB0dzYQJE4iOjiYwMJAlS5b4r1myZAkPPfQQR48e5b777uPee+8F4OGHH+bBBx/E7XYTHBxMbm5uU4QqIiINYVqpdevWmX79+pmIiAizYMECu8OxTUlJiUlISDDR0dHmlltuMQsXLrQ7JNvV1dWZgQMHmm984xt2h2Kr6upqk5qaaqKiokz//v3NBx98YHdItpk3b56Jjo42MTExJj093Rw7dszukJrNlClTTPfu3U1MTIy/rbKy0gwbNsy43W6TlJRkqqurL3mfVrmC+8SJEzzxxBPk5+eze/duXnvtNT799FO7w7KFw+HgueeeY9euXWzdupXFixe32/filLPX/7RXM2bM4L777uPTTz/lk08+oX///naHZIt9+/bx4osvUlRUxM6dOzlx4kS76qmYMmWKf/3aKRdaE3cxrTJZfPjhh0RERNC7d28cDgcTJ04kLy/P7rBsERoaysCBAwHo2LEj/fv3p7y83Oao7HO+9T/t0aFDh3j33XfJysoCvpp+3h516tQJh8PBkSNHqKur48iRIzidTrvDajZDhw6lS5cuZ7RdaE3cxbTKZHH6Ggz4at1Ge7dv3z4+/vhjhgwZYncotjnf+p/2yOv10q1bN6ZMmcLtt9/Oo48+ypEjR+wOyxZdu3blqaeeolevXvTs2ZOgoCCGDRtmd1i2utCauItplf+i2nv3wvl8/vnnjBs3joULF9KxY0e7w7HF6et/2nNVAVBXV0dRURHTpk2jqKiIG2644bK6GtqivXv38vzzz7Nv3z7Ky8v5/PPPeeWVV+wOq8W43DVxrTJZnL1uo7S09IyV4O3NF198QWpqKg888IB/unJ7dGr9T58+fUhPT+ett94iIyPD7rBs4XK5cLlcDB48GIBx48ZRVFRkc1T22LZtG3feeSfBwcEEBgYyduxYtmzZYndYtjq1Jg44Y03cxbTKZDFo0CCKi4vZt28ftbW1rFy5kpSUFLvDsoUxhocffpjo6GiefPJJu8Ox1fnW/yxfvtzusGwRGhpKWFgYe/bsAWDjxo3ccsstNkdlj6ioKLZu3crRo0cxxrBx40aio6PtDstWF1oTd1FXarrWlbZ27VoTGRlpwsPDzbx58+wOxzbvvvuuCQgIMAMGDDADBw40AwcONOvWrbM7LNsVFhaakSNH2h2GrXbs2GEGDRpkbr31VjNmzBhTU1Njd0i2eeaZZ/xTZzMyMkxtba3dITWbiRMnmh49ehiHw2FcLpdZunSpqaysNImJifWaOhtgTDvv3BURkUtqld1QIiLSvJQsRETkkpQsRETkkpQsRETkkpQsRETkkpQsRETkkv4/6SmlQxz22aEAAAAASUVORK5CYII=)
This looks nice and smooth, but you might say. But if you think about it, for a particle in the middle of two bins, say at the right edge of the first bin at 0.95, is it really justified to assume that it only contributes to the grid point at 0.5? It virtually has the same distance to the grid point at 1.5!
[latex]
\begin{figure}\begin{center}
\begin{tikzpicture}
\fill [red] (0.95,.15) circle (3pt);
\draw [line width=2pt] (0,0) — (11,0);
\foreach \x in {0,…,11}
{\draw [line width=2pt] (\x,0)–(\x,0.3);
}
\foreach \x in {0,…,10} {
\node () at ($(\x,-.3) + (0.5,0)$) {\x};
}
%\fill [red] (0.95,.15) circle (2pt);
\end{tikzpicture}\end{center}
\caption{A grid with one particle at 0.95}
\end{figure}
[/latex]
This problem can become even worse. Let us define point positions $$x_i=10(1-1/i)$$. Here see a figure with particles with $i$ between 10 and 25:
[latex]
\begin{figure}\begin{center}
\begin{tikzpicture}
\foreach \x in {10,…,25}{
\pgfmathparse{(10*(1-1.0/\x)-9)/.05833}
\fill [red] (\pgfmathresult,.15) circle (3pt);
}
\foreach \x in {0,1,…,10} {
\node () at ($(\x,-.3) + (0.5,0)$) {\x};
}
\draw [line width=2pt] (0,0) — (11,0);
\foreach \x in {0,…,11}
{\draw [line width=2pt] (\x,0)–(\x,0.3);
}
\end{tikzpicture}\end{center}
\caption{A grid with particles of decreasing spacing}
\end{figure}
[/latex]
You can see how the particle distance (which is inversely proportional to the intuitive density) changes smoothly. The density per bin can be calculated analytically to be $\rho(x)=\frac{x \delta}{(-10 + x)^2}$ with $\delta$ the bin width. Let us plot the point positions and the resulting histogram for $i$ between 10 and 25. The histogram however has a step-shape as you can see below and does not compare too well to the analytic result (grey vertical lines indicate particle positions, the green curve is the analytic particle density):
x=10*(1-1/numpy.arange(10,25.)) # the points
density=lambda x: (10/(10/(-10 + x) + (10/(-10 + x))**2))**(-1)*.05833
for i in x:
axvline(i,color="grey",linewidth="1") # grey lines for the point positions
hist(x,bins=10) # histogram of point positions
plot(numpy.arange(9,9.6,.01),density(numpy.arange(9,9.6,.01))) # analytic density function
xlim(8.99,9.6)
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXIAAAEACAYAAACuzv3DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVGX/x/H3sLuhhoqmoimkCImWW+ijWJap5ZKPmrmWCmakZf20njAhLbXSUNPM0jST0kjKXdNE3MV9SXNBxQVCQEBGQWDO74+jo4TsMxwGvq/rOhdn5mzfI/jh5p5z7qNTFEVBCCGExbLSugAhhBDFI0EuhBAWToJcCCEsnAS5EEJYOAlyIYSwcBLkQghh4STIhRDCwuUb5Ldv32bo0KG0a9eOjRs35ljep08f6tSpg5+fn1kKFEIIkTddfjcEnTp1ikcffZT4+HimTZvGd999Z1wWGRmJXq/Hx8fH3HUKIYTIRb4tcnd3d6pWrUpCQgITJ07MtmzTpk34+voyZcoU0tPTzVakEEKI3OXbIgeIiopi5MiRuLm58c0332RbdufOHcaMGUP79u157bXXsu9cpzNttUIIUQ4UeuQUpYAMBoPi5eWl6PX6HMtu3rypDBo0KNt7gEwyySSTTEWcCqPAV63odDqeeeYZKlasSHJyMqhHAiAxMRFvb++HbqcoSo6pU6dOD33fkqbAwMAc702ePFnzugpSZ1EmU5ybqWox9XFN+X0z1zkWdb95nZspatVyH0XJkcIeqzDrF3Td/NYrCpv8Vpg3bx7nz59nwIAB+Pr6cuzYMaZPn05ISAjdunXDw8ODBg0a4O/vX6QChBBCFE++Qf7mm2/meC8kJATgoZcjCiGEKFma3BDUsGFDLQ5rdmX5Mkw5N8tUls+trOZIUeTbIje1Hj16sH79epYuXVrShza5wMBArUsokJKss3v37qxbt67EjmcKZTnsyvK5PfbYY1qXUGqUeIt8/fr1JX1IUYLk+ytEySvxFvk9Rf10VpRect+AENqQQbOEEMLCSZALIYQGjv9znCnbp5B6J7XY+9Ksa0UIIcqzD7Z+wLqz69Bn6JneZXqx9iUtchP75ZdfCn1Z1M6dO+nbt2+hj3XkyBHatWvH9u3bC73tw7zxxhuEhYWZZF9CiNxtv7iddWfXUdmuMuOfHl/s/UmQm1j37t2Jjo7Od70//viDS5cuAdCmTRtmzpxZ6GO1aNECg8Fgsg8ZAwIC6Nq1KwBZWVksWbLEJPsVQtynoDBxizqS7P95/x+1KtUq9j4lyE2sYsWK+a6TkpLC+PH3fwvb2dkV+eaGghyvoOrWrWvcX2BgoPEXjRDCdE5xin1X91GrUi2TtMahFPaR64JMdwmbMrlwlzhmZGTw4Ycf4uLiwrp161i+fDnJycm8/fbbdOzYkX379nH48GF27tyJs7Mz33//Penp6WzdupVXX32VPn36ZGsdf/vtt/j5+bFp0yaee+45QkNDWbVqFYMHD+batWt899139OvXj9DQUI4ePcrvv/8OwMqVK0lKSiIiIoLAwEBcXV2z1bllyxa2bNlC3bp1OX/+vPH9Xbt2ce7cOSIiIhgyZAju7u7873//o3bt2sTHx7Nx40bWrVuHp6cnGzZsICEhgYULF/Lmm2/SsmVL3n77bV555RU6d+7M/v37jbWsXbuWLVu2sHfvXlxcXPD396dZs2aMGTOmqN8aIcqlTEMmW9kKwOROk6lsV9kk+5UW+QN27txJZmYm/v7+WFtbs23bNh577DEcHBw4c+YMS5cupW3btoSGhgIQHBzM6NGj6d+/P4sXL86xv1GjRtG5c2cMBgMAN2/eZOrUqXTv3h1HR0dGjRpF8+bNadu2LSkpKQDs3r2bQ4cO4evrS5cuXVi7dm22fd6rb/r06bz11ltUrnz/B2HevHkMGzaMt99+m7Fjx+Ls7IyTkxMnTpxgxowZDBkyxHhH7dKlS3n55Zf5+eefcXBw4PHHH6dSpUoA1K9fH29vb9q3b0///v2ZP38+BoOB6tWrA1C9enUJcSGK4PvD35NAAo2rN2bUk6NMtt9S1yIvbCvalDp37oynpychISGkpKSg1+sBqFChAu3bt6dChQq4urqSkJAAwOHDh9m4cSNHjhwhNfXhlxC9+eabfP3113Tt2pXz58/nePgGZO8e+emnn+jSpQsAw4cPz7Hu9u3b8fT0NL6uXbs2AIcOHeLGjRssXbqUzMxM3NzcuHnzJhUrVuSpp56iWrVqNG7cmB07dgDQtGlTvLy8+Oyzz+jTpw+AMcjvuXfTVuXKlenfvz+LFi3Cz8/PpN05QpQXtzJuMTl8MgCfPPMJtta2Jtu3tMgfcOTIET744AMGDBhA48aNje8/2F1ybz4zM5NXX32VJk2a8Pzzz+e6z549e3LkyBF27dpF/fr1H7rOg/vPyMggMjLS+DouLi7buhcuXODKlSs59nHr1i2Sk5MZNmwYI0aMIDQ0lAoVKqDT6YyB/OD8pEmT+PTTT/Hz82P+/Pk56vi3MWPGMH/+fMLCwujdu3eu6wkhHm723tnEpMZQhzr08+hn0n1LkD9g2bJl1K9fnzt37nDt2jXu3LmDXq8nKyvLGICKomAwGDh16hSRkZE0bNiQc+fOkZGRwa1bt3IMPWBjY8PIkSMZNGgQ/fv3N75vZ2dHamoqycnJ2bbp2bMnS5YsYcuWLRw6dIjjx49n29+LL77I4cOHOXv2LAaDAb1ej16vp3379ly4cIEvvviCy5cvGx/Jl5WVlW37e8eaOXMm/fr1Y+XKlRw6dAgAg8FgXG5nZ4derzc+RKRp06a4uLjw008/4e7uXux/ayHKk4RbCUzfpV4r/hzPYaUzbfRKkD+gb9++LF68mLfeegsvLy9OnjxJQkICBw8eZPv27URFRbF//37279+Ps7MzLi4utG3bFhsbG9LS0rh27Rpr165Fp9MZP7gEta+8VatWxj5mgK5du/Lee++h1+vZsGED58+f5+zZszz33HOMHj2aESNGsGTJEjp37pytxtq1a7NgwQJ69erF2LFjsbOz4+DBg2RkZLBs2TLCwsIYMGAAHh4eJCcnExERwf79+4mKiiIiIoLjx49z+fJlFi9ezPjx4zl48CDvvvsuV65c4dixY2zfvp1bt27RsWNHfvvtN2NXDKjXmd/r9hFCFNynOz4lJT2F5xs/TyMamXz/pa6PXEve3t4PveTu5MmTxvkHR/fbtm2bcX7YsGEAuLq65mgFJycnM3DgwGzvzZkzxzj/2Wef8dlnnxlfBwQEEBAQkGudr7322kP72rt06ZIjaMPDw43zD34ge/r06RzbHz582Djfvn17zpw5k215SkpKjvMQQuTtwo0LfBX5FQDTn53O6vOrTX4MaZGbUXx8PIcOHeKnn36iZ8+eWpdTZH///TcnTpwgNjaWWrWKf/OCEOXJe3+8x52sOwxuPpiWdVqa5RjSIjejDRs2EBQUxJIlS7C1Nd0n1CVt1qxZnDx5Um7fF6KQ/rzwJ6tOraKibUWmP1u88VTyIkFuRkOGDGHIkCFal1Fs9z44FUIUXKYhk7c3vg3A/zr8j7qOdc12LOlaEUIIM/j24LccjztOw2oNTXYrfm4kyIUQwsQSbycyadskAGY+P5MKthXMejwJciGEMLGg8CASbifQuWFn+jTtY/bjSZALIYQJ/XX9L+ZFzsNKZ0XwC8El8izbPIP89u3bDB06lHbt2rFx48Zsy2JiYpg6dSorVqwo0PjbeenRowc6nc6sU48ePYpVoxBC5EdRFN7e+DZZShZ+T/nR3Ll5iRw3zyC/ePEic+fOZfny5cYR/+7x9/fH19eXXr16Zbu5pSgevMnGXApyjJ9//hkXF5cc77/44oscPHjQHGUJIcqQNWfW8EfUH1RzqMbHnT8usePmGeTu7u5UrVqVhIQEJk6cmG1ZVFQUtWrVwsHBIcd4IEWnmGkqmJdeeumhA1LNnTsXLy+vXLf79ttvC3wMIcqzh/31DRTpr+yibFfYbQqz/h3u0GtBLwCSwpKoWalmrvs0tXyvI4+KiuL999/Hzc3NeD1xQkICVlb3fwfExcWRlJREtWrVcmwfGBhonPfx8Sl+xWaU2/Csjz32WK7bXL16lcmTJzNqlOnGFhairCqJv761sp3tUA2IASLzWzu78PDwbMNpFFa+Qd6oUSO2bt1Ky5YtuXXrFhUrVsTJycn4sAQAvV7/0BCH7EFe2t37bbllyxYCAgJo06YNw4YNMz7IoVOnTvz6668kJSUxZ84cFi5cyIULF4iNjWXWrFn0798fg8HAjz/+iMFgwMnJiddffx17e3uuXLnC77//zqVLl/jiiy8YNWoUgwYNYsKECQwcOJD58+ezZMkSduzYwSOPPMLq1av54IMPcHd356OPPiItLY369euzePFiAgMDsbe3Z8aMGbRv356vvvpK4385IQrrwb+UgyjMX87F2y7oIcc3wTFqHWcPLUDRwdq9YGiTy4oPb5H7+PgYG7pBQUEPXScvBbpqRafT8cwzz1CxYkXjk2xcXV1JTk5Gr9eXuWFN69Spw/bt21m0aBHNmjUjMzPTGPLLly9n+PDhLFq0CCsrK1555RUAxo8fT7169ejfvz/vvPMOAQEBREZGMmPGDEAdNrZ27dp89tlnODo6MnToUDp27GgcYGvz5s20aNGC3377jZEjR9KxY0dCQkKoVq0aTZo04dy5c0yYMIGFCxfy4Ycf4u3tzd69e1m0aBG3bt3S5h9KCAE6A7zkhwEDRL4BV3MLcfPJM8jnzZvH+PHj2bdvH76+vhw7dozRo0cDMG3aND766CMWLFjAvHnzSqTYkuLh4YG9vT1OTk4kJiZme3JO1apVad26NXq9ntatW2fb7uLFi1y7do0KFdSL/5999ll++eUXACIiIrC3twfA09PT+BdNpUqVaNmyJQ0aNKBChQps27aN3377jTNnzhifOlSxYkUaNGiAvb09rq6u2Nra4uLigr29PbVr1zY+sUgIoYEnv4P6e6hMZdj6qSYl5Nm18uabb+Z4LyQkBFBb5LNnzzZPVaWEtbV1ttY4qB9s3hsPfNmyZbz00kvGZTY2Nly7do2bN29SpUoVPDw8jMs6d+7M7t27efHFF7l16xYdOnQAsj+VJzk5mZEjR/L111+TnJxs7DN72BOKHvTvh1kIIUpIpTjool4I8gIvEJpeVZMyStkNQTozTQXzsEBUFMU4gToSoL+/P8HBwRw9ehQAW1tbbt26ReXKlXnyySeN19wfOXKEfv3URzr5+fkRHx/Pzz//zIYNG7IF8r3ule3btxMfH0+NGjU4f/48GRkZ6PX6fIP6wc8rhBAl6Pl3oUISnOuKBx75r28mpSLIu3fvXiqOsXr1anQ6HWvWrGHfvn3Ex8ezbt06zp49y6ZNmzAYDEybNo1JkyZx7do1fH19AXjhhRcYNWoUNjY2LF68mFmzZjFp0iSio6OZNEkdb2HRokWsW7cOPz8/mjRpwqeffsrx48c5d+4coaGh6PV6OnbsyI0bN+jSpQt16tQhLi6O5ORktm3bxvHjx7l06RKrVq3i+vXr7Nmzh927d3P9+nXWrFlj1n87IcRDPLYVvH6EDAdYNw9dIRqNplYqhrFdt26d1iUA0KtXr2xP97l58yaQvYvpxo0bObZ78LFunp6e7NmzJ8c6Dg4OXL16FYDU1FTmzJnDE088weXLl7Otd+TIEeP8G2+8AcAPP/xgfG/8+PGMH39/JLV7/ehCiBJkcxteVP9/EhEANxrnvb6ZlYoWeVl38uRJNmzYYPwlkJGRQaNGpn9unxCihHSeDE5n4bo77P4/raspHS3yss7Dw4O+ffvSpUsXHB0d8fDw4PPPP9e6LCFEUdTbC0/PBIMV/PY9ZNlpXZEEeUmZMGECEyZM0LoMIURx2KRBr9fAygA7J8DVtlpXBEjXihBCFJxPINQ8DdebQnjh78A0FwlyIYQoiLr7wftztUvl9+8h00Hriow061opicHWhRDCJGzSoPdwtUtl13twpZ3WFWVT4i3ykrhmXGhHvr+iTOoUBDVPQfzjsK3kxhkvqBIP8nXr1hEYGJjtjklLnCzlHEq6ztJyT4AQJvNoJLT/TB3Z8PfvIdO8D1IuCukjF0KIXNzhDrw8WO1S2fMOXPbWuqSHkiAXQohcbGIT1DgDcR7w51Sty8mVBLkQQjxM0984yEHItINfQ0pll8o9EuRCCPFvVa5Bz5Hq/JYZ8E9zbevJhwS5EEI8SGdQLzWsmEBjGsO+sVpXlC8JciGEeFDb2dD4D9DXoDe9QSn9MVn6KxRCiJLifBS6vK/Or15EFapoW08BSZALIQSArR76vgo2dyByNPzdU+uKCkyCXAghUKDHGKj1lzog1uaZWhdUKBLkQgjR8nto8QNkVIBffoGMilpXVCgS5EKI8s35KHS/+zjHtV9DnKe29RSBBLkQovyyT4H+/cA2DQ6NgKPDtK6oSCTIhRDllKLe9ON0FmKbw/q5WhdUZHkGeUZGBq+//jru7u5MmTIlx/I+ffpQp04d/Pz8zFagEEKYRZt54PELpFdR+8VL8S34+cnzwRKRkZEEBwejKArNmjVj0KBBxqe/R0ZGMm7cOMLCwkqkUCGEMJm6+6HreHV+9XeQ8Li29RRTni1yb29vHB0dqVq1Kp6entjb2xuXbdq0CV9fX6ZMmUJ6errZCxVCCJOoDAx4GawzYJ8/nOyvdUXFVqBHvcXGxuLp6UndunWN7wUEBDBhwgTGjBlDSEgIr7322kO3DQwMNM77+Pjg4+NTrIKFEKKoMsmE/oDjVYhuX2quFw8PDyc8PLzI2xcoyENDQ/n0009zvG9nZ0dwcDCjR48uUJALIYR2FNazHlyA5Hqw4lfIstO6KCB7IzcoKKjQ2+d71cqqVavo27cv9vb2xMXFkZKSAoCiKAAkJibi7V06n5ohhBBGrRZwiEOQAawIA72z1hWZTJ4t8q+++oqZM2dSo0YN0tLSGDVqFHv37iUkJIRu3brh4eFBgwYN8Pf3L6l6hRCi8BpEQLe7w9GuAa610rQcU8szyP39/XOE9Nix6j/Gxo0bzVeVEEKYStVo6P9fsM7kaZ5mz7E9WldkcnJDkBCi7LLVw4A+UOk6nH+OLnTRuiKzkCAXQpRNuix1WNpHD0FiYwj9GWusta7KLCTIhRBlU9d3oelquF0dlq+D249oXZHZSJALIcqeNnOh3WzIsoWfwyChidYVmZUEuRCibHl8Dbzwtjr/+2K41EnbekqABLkQouyocxD++wpYGWBbEBwbrHVFJUKCXAhRNlSNhldfBLtbcGQobJ+kdUUlRoJcCGH5KsbD4K5QJRYu+MCabwGd1lWVGAlyIYRls0uFV3tAzdPwjyesWFVqxlApKRLkQgjLZZ2uDklbbz/caAg/boK06lpXVeIkyIUQFsmAAfoMhcZ/QGotWPYH3HxU67I0UaBhbIUQonS5OySt5wFIc4QfN0Kiq9ZFaUZa5EIIy+MTyAEOQKY9/LQaYltqXZGmJMiFEJal/Wfg8zE6dBD6c7m44Sc/EuRCCMvRLhiemwiKjl70gtO9ta6oVJAgF0JYhtbz4YV31Pk1C2lBC23rKUUkyIUQpd+T30GPN9X5dfPg0Eht6yllJMiFEKWb11J4yVed3/glRI7Rtp5SSIJcCFF6PbEcer0OOgX+mAF739a6olJJglwIUTq1XAQvD1FHMvzzY9g1QeuKSi25IUgIUfq0nn+/T3zLp7DzA23rKeUkyIUQpcvTs9THtIHaJy7dKfmSIBdClB7/+QSeDVDn186HA29oW4+FkCAXQpQCCnT+CDpNBUUHq7+Dw69rXZTFyPPDzoyMDF5//XXc3d2ZMmVKtmUxMTFMnTqVFStWEB0dbdYihRBllwED9BijhrjBGlb9KCFeSHkGeWRkJMHBwezdu5cFCxYQFRVlXObv74+vry+9evVizpw5Zi9UCFEGWacTSii0XqAOgLUyFI6/qnVVFifPIPf29sbR0ZGqVavi6emJvb29cVlUVBS1atXCwcGB48ePm71QIUTx9ejRA51Oh06nPgbt3nxRp8LuIxv7FBjcjb/4Sx2KdtkmGTuliArURx4bG4unpyd169YFICEhASur+78D4uLiSEpKolq1ajm2DQwMNM77+Pjg4+NTvIqFEEW2fv16rUtQVfoHBneDOoepTGVSv4+Af7y0rkoz4eHhhIeHF3n7AgV5aGgon376qfG1k5MTBoPB+Fqv1z80xCF7kAshSgsFCLr7tTgKuw8dVAcGdwCnc5DYmNcf6cGcchzikL2RGxQUVOjt872zc9WqVfTt2xd7e3vi4uJISUkBwNXVleTkZPR6Pe7u7oU+sBCiHKoHjEQN8ZiWsGgXj/CI1lVZvDxb5F999RUzZ86kRo0apKWlMWrUKPbu3UtISAjTpk3jo48+wsXFhXnz5pVUvUIIS9XsF3gZNXXOPQ+//ALpjlpXVSbkGeT+/v74+/tne2/s2LGA2iKfPXu2+SoTQpQRivpUn+feV18eANavBYOtplWVJXJDkBDCfKwy1DFTnvpWfb0Z2A0gIW5KMvqhEMI8HG7AoB5qiGc4wMpf7oa4MDVpkQshTK/GKRjYU/1QU19TfdL9lXZaV1VmSZALIUzr8bXQ91WwvwmxXvDT75DcQOuqyjQJciGEiSjwn2nwTID6RJ8T/eH3xZBRSevCyjwJciFE8dnqodcI8Fyhjl649RPY8QGgy3dTUXwS5EKI4nE6A/37gvMJSK8Cvy6HMy9pXVW5IkEuhCiyv/gLfFup/eHxTWDFKrjeTOuyyh0JciFE4VllQJf3WclKsEftD1/9HdyponVl5ZIEuRCicKpcg/8OgAY7scIKw4ZZsG8s0h+uHbkhSAhRcK4bYXQLaLATUh5lOMNh3zgkxLUlQS6EyJ/1HXj+PXUM8UrXIepZ+OYwLrhoXZlAulaEEPl55Bz0HQh1D6jP1PxzCuyaCIq0A0sLCXIhRO6eWA4vjgb7VEhqAKE/wZWnta5K/IsEuRAiJ4cb0N0fmoeor0/2gzULIe3hTwIT2pIgF0Jk12gL9B4OjlfhTkXYGAyHRiIfaJZeEuRCCJXNbejyPrSbo76+0hZWLYNEN23rEvmSIBdCQJ2D8PJgqHkasmxg+2TY+T4YJCIsgXyXhCjPbNKg08fqo9issuC6u9oKj3lK68pEIUiQC1FOXeYy+D0JNU+pIxbueUcdtTCzgtaliUKSIBeivLEFOr/LIhZBTeB6U3XccLms0GJJkAtRnjwGvAg4zUKHDmXHRLU/PNNB68pEMcitWUKUA9f11xkaNhSGAU7AP08wkpGwdZqEeBkgQS5EGaYoCosPL6bpvKYsO7YMMoGtwMID1KWu1uUJE8k3yCMiIujSpctDl/Xp04c6derg5+dn8sKEEMVz6vopfJb6MGL1CBJvJ/Jco+dgPrADyLLTujxhQvn2kXfs2BG9Xp/j/cjISMaNG0dYWJhZChNCFM3N9Jt8vP1jgvcFk2nIpGbFmgS/EMxAz4FYDZU/wsuiAn1XHRxy9qFt2rQJX19fpkyZQnp6uskLE0IUjqIoLD+2nCZfNeGLPV+QZchi1JOjOO1/mlefeBWdTm6xL6uKfNVKQEAAEyZMYMyYMYSEhPDaa689dL3AwEDjvI+PDz4+PkU9pBAiF0djj/LWhrfYEb0DgLZ12/JV969o9WgrjSsTBREeHk54eHiRty/W5Yd2dnYEBwczevToAgW5EMK04vRxTN42mYWHFmJQDNSsWJMZXWYwrMUwrHTSjWIpHmzkBgUFFXr7Qgd5cnIyVatWRVEUdDodiYmJeHt7F/rAQoiiS8tMY86+OXyy4xNS0lOw1lnzVpu3+Ljzx1RzkKFmy5t8g/zo0aOcOXOGEydOkJWVxYwZMwgJCaFbt254eHjQoEED/P39S6JWIco9RVEI/SuUiVsmciHpAgDd3brz+XOf06xmM42rE1rJN8i9vLy4evWq8XVIiDrQ/MaNG81XlRAih53RO5m4ZSK7L+8GwKOmB7O6zuL5xs9rXJnQmtyiL0QpdzLuJCGEEPh9IAA1K9ZkSucpjHhyBDZW8l9YSJALUWpdTr7M5PDJLD26FAMGKtlW4j3v93j36XepYl9F6/JEKSJBLkQp80/qP0zfOZ2vD3xNelY6NlY2PKU8xZqxa3Cu7Kx1eaIUkuuThCgl4m/FM/GPiTSa04jgfcGkZ6UzwGMAp948RQ96SIiLXEmLXAiNJaUlMXP3TIL3BZN6JxWAXk16EeQThFdtL42rE5ZAglwIjSTcSuDLvV8yd/9cUtJTAHjB9QU+9vmY1nVba1ydsCQS5EKUsDh9HDN3z2T+gfnGFvgzjz3Dxz4f096lvcbVCUskQS5ECbmScoWZu2fyzcFvuJ15G1Bb4JM6TsK7vtwdLYpOglwIMzt1/RSf7/6cH4/9SIYhA4CXHn+JSR0nSReKMAkJciHMZO+VvczYNYPfTv8GgJXOiv4e/fmgwwe0qN1C4+pEWSJBLoQJGTDw++nfmblnpnFIWXtre4a3GM573u/h+oirxhWKskiCXAgT0N/Rs+TIEr7iKxJXJALgaO/ImFZjGNduHLUr19a4QlGWSZALUQyXky8zP3I+3xz8hhtpNwBoWK0h49qO4/WWr+No76hxhaI8kCAXopAURSHiUgRz988l7HQYBsUAgHd9b+pdrsfyt5bLYFaiRMkt+kIU0K2MW3x78Fu8Fnjhs9SHX0/9ipXOilc8X2HPiD3sen0XzWgmIS5KnPzECZGPv67/xTcHvmHp0aUkpycD4FzJGb+n/PBr5cejVR7VuEJR3kmQC/EQ6ZnphJ0OY8GBBWy/tN34ftu6bfFv40+/Zv2wt7HXsEIh7pMgF+IBp+NPs+jQIpYeXcr1W9cBqGRbicHNB+P3lB8t67TUuEIhcpIgF+We/o6eX/76he8Ofceuy7uM7zd3bs7op0YzqPkgufpElGoS5KJcUhSF3Zd3s/ToUn4+8TM379wEoLJdZQZ6DmTkkyNp/WhrdDqdxpUKkT8JclGuRCdH88PRH/jh6A+cTTxrfN+7vjcjW46kn0c/KttV1rBCIQpPglyUeWmk8f3h7/nx+I+1Hfq4AAAPgklEQVRsu7ANBQWAR6s8ypDmQxjmNQz3mu4aVylE0UmQizLpTtYdNp7byI/HfiSMMDJXZwLquCd93Psw3Gs4XRp1wdrKWuNKhSg+CXJRZmQZsoi4FMHPJ34m9FQoibcTjcs6NejE4OaD+W+z/1LNoZqGVQphevkGeUREBB9//DFbtmzJ9n5MTAyLFi3Czc2Np59+GhcXF7MVKURuDIqBvVf2suLEClb+tZLY1FjjMs9angx+YjCxW2P5cviXGlYphHnlG+QdO3ZEr9fneN/f35+vv/4aR0dHAgIC+OKLL8xSoBD/ZlAM7Luyj9C/Qgk9FUp0crRxWaPqjRjgMYBXPF+huXNzAIK2BmlVqhAlokBjrTg4OOR4Lyoqilq1auHg4MDx48cLdLADBw7g7OwMgL29vSaTra0tOp2u2BNgkv2YeypNdRarFisdugY6dN10WL9njfdib2btnaWGeDKwG1gIUeOimNZlGl61vTT5NzDXsUy1X1E2FamPPCEhASur+78D4uLiSEpKolq1nH2PgYGBxvn4+Hji4uIAuHPnTlEOLcoTa+AxoOnd6cGrApOBv+5OV+DuhShCWKTw8HDCw8OLvH2RgtzJyQmDwWB8rdfrHxrikD3Ig4ODH1iSVpRDF9N8YPzd+eL+zw8q4j7utYpKKnn+XWdJH//ftdyTy/HtU8B1IzQNA7f14JByf9mNhvDXf9Xpahvun0tBjnvveOY+//x+Lop6/KL+vOV2fFGa+Pj44OPjA0BQUOG7Agsd5CkpKTg6OuLq6kpycjI2Nja4uxflGlwtBhyy1eCYIl/VLsLja6DJGmgYDtYZ95f98wSc6gOn+0CsFxJEQuSUb5AfPXqUM2fOcOLECbKyspgxYwYhISFMmzaNjz76CBcXF+bNm1cStYqyQgfU26W2uB9fC7WP3V9msIJL/4HTveB0b7jRWLMyhbAU+Qa5l5cXV69eNb4OCQkBwNXVldmzZ5uvMlG2VIznGMegL9AYqNjh/rL0ynDuBfi7J5ztDredtKpSCIskNwQJ87DKhLr71f7uxpugbiSrUOCJu8sTG8OZHnC2B1zsBFkytrcQRSVBLkyn2kVotAUab4ZGf0CFpPvLMu1oZFOPqI1RcBZIOKdVlUKUORLkougckqDhNmj8hxrcTv8K5wQ3ONdV7Ta56MPQD78gcG+gJqUKUZZJkIuCs70F9XdBo63w2J9Q5yBY3b8MlTRHuPAMRD2nBrh8UClEiZAgF7mzSYO6+9RLAh/bBvX2gM0DN3Jl2cCl9nD+OTW8r7UCg/xICVHS5H+duM/mNtTbBw22q+Fdfw/YpN9frujg2lN3W93PQnQHyKikWblCCJUEeXlmnwwuu6BBBLjsgLqR2W/GAfUmnIud4KIPXOoEtx/RpFQhRO4kyMsNBapGgwvqVN8LnI+D7oHbvhUdxLSASx3hYmf1xhy5pluIUk+CvKyyvgO1D7OHPfDfAWrL2/HqAyscgyxbuNpaDe7o/8Blb0iThy4IYWkkyMsEBRyvqDfg1Nur9m0/egBs0tkE4Hl3tdvV4fINiAaiI9QPJzMraFe2EMIkJMgtkX2KeulfvX3qVSX19kGVmJzrXW9Ki5oOHFk9Rm1tx7uDcu8Zlf8p0ZKFEOYjQV7a2aSB81H1g8hHI9VWd42/s/dtA9yupg7terUtXH5a/Xr7EXoHBnHk0ChtahdClAgJ8tLEJg2cj6mt7UcPqt0jNU+CdWb29bJs1atJ7gX3lbaQ6AZKgR74JIQoYyTItWJ3E2ofhTqH1Kn2Yaj5V87QNlhBXDO41lr9YPJqG/inuQwyJYQwkiA3OwWqXIPaR+5OqNMjVXN2j9wL7Zin1BtvrrWC2BZy040QIk8S5KZkc1ttVTsfU6faR9WvFRNyrptlA3GeEPPk3aml2tKW0BZCFJIEeVHoskggAdxXQa3jUOuEenPNI2ezDyJ1z+3qass6tgXEfgmxQHwqZNmVeOlCiLJHgjwvOoN6N2TNk2pY1zqpztc8xVxuw4C52dc3WMF1d7Vl/U9z9QPJf5pDSj3uP2vyy7tfJcSFEKYhQQ7q02yqR0GNU2rXSM27X2ucArtbD92kClW4ebY9xD2hPiA4zlO9TjvToYSLF0KUd+UryO3SoAbgBNQMgBqn1bB2OptzsKh7Up3VkI7zgOse6vz1Zrz7/mwCl08uyeqFEOKhyl6QW2WqjxxzOnN3+lu9gabG6X/d/fhJ9u2SXNQW9fVmdyd39bWM9ieEKOUsM8h1WVD1MjxyTv2A0ens/a/Vz+e8FvueTBtIyIQE4HoAxDdVwzq+iVwtIoSwWKU3yK0yoNolNayrn78b2nen6lHZn1Tzb0kukPC4OiW63Q3sJpC8BpRxd1eaUiKnIYQQ5qZhkCtQIRGqX1CD2TidV79WjX74pXz33KyjPtw30Q0SXdX5hMfV+VxH9JNb2IUQZU+BgnzhwoU4OakPGOjbt2+2ZX369GHv3r307NmTb775Js/9nMg4AQPga76G92eBQ0ruKys6tWWd6Ko+xDfR9f50oxHcqVyQ0oUQoszLN8gPHDjAyZMnmT17NpMnT8bHx8cY6pGRkYwbN46wsLACHeyG4Qa4wz/8Aw5AehU1lG88dvdr47tfG0FSQ7lhRgghCiDfIA8PD6dNmzYANG7cmM2bNzNw4EAANm3axA8//MCQIUOYMGEC9vZ5D+TU1KYp/AKj+o3i2xnT7l4RostzGyGEEHnLN8ijo6Np0qQJAI6Ojmzbts0Y5AEBAUyYMIExY8YQEhLCa6+9lmP7wMBA43xyfDKchLr96sqzIIUQ4q7w8HDCw8OLvH2+Qd6sWTNiYtTrr2NjY2nQoEG25XZ2dgQHBzN69Oh8gzw4OLjIhQohRFnl4+ODj48PAEFBQYXePt/LODp06MDp06cBtXXeoUMHUlLUDykVRR2GNTExEW9v70IfXAghRPHl2yL39PTEycmJ7777jkaNGuHg4MDo0aMJCQmhW7dueHh40KBBA/z9/UuiXiGEEP9SoMsPP/zww2yvQ0JCANi4caPpKxJCCFEocoeMEEJYOAlyIYSwcBLkQghh4STIhRDCwkmQCyGEhZMgF0IICydBLoQQFk6CXAghLJwEuRBCWDgJciGEsHAS5EIIYeEkyIUQwsJJkAshhIWTIBdCCAsnQS6EEBZOglwIISycBLkQQlg4CXIhhLBwEuRCCGHhJMiFEMLCSZALIYSFkyAXQggLJ0EuhBAWLt8gX7hwIb/++iu//vprtvdjYmKYOnUqK1asIDo6ulAHvXDhQuGqtBjhWhdgRuFaF2BG4VoXYEbhWhdgNmU3RwovzyA/cOAAJ0+epG/fvhw7doyEhATjMn9/f3x9fenVqxdz5swp1EEvXrxYpGJLv3CtCzCjcK0LMKNwrQswo3CtCzCbspsjhWeT18Lw8HDatGkDQOPGjdm8eTMDBw4EICoqilq1agFw/PjxIhz6YBG2Ka7C/eUghBAWQcnDW2+9paxevVpRFEUJCwtTxo4dqyiKosTHxytPPvmkcb0WLVooN27cyLYtIJNMMskkUxGnwsiza6VZs2bExMQAEBsbS4MGDQBwcnLCYDAY19Pr9VSrVi3btmqWCyGEKKwvv/yyUOvn2bXSoUMHFi9eDEB0dDS9e/cmJSUFR0dHXF1dSU5OxsbGBnd394duL2EuhBDmp1PySdtPPvkEZ2dnANq0acP06dMJCQnh3LlzzJ07FxcXFwYMGEC9evVKpGAhhBDZ5RvkQgghSje5IUgIISxcnn3kxZWZmcnatWtJS0sjKSmJ0aNHG5fFxMSwaNEi3NzcePrpp3FxcTFnKSaX17klJSUxdepUWrZsyaBBgzSssvDyOq81a9bw+eefY21tzfz583P9bKQ0y+v8fvrpJ77//ntq1KjBsmXLsLa21rDSwsnrvO4ZMWIEQ4cOpVOnThpUWHT5nVufPn3Yu3cvPXv25JtvvtGoyqLJ79yio6PZvn07bm5utGvXLtf9mLVFvn79emJjY3nllVc4cuSI8QoYKN4NRaVBXudWrVo13NzcyMzM1LDCosnrvGxsbIiIiKBbt25MmzZNwyqLLq/zc3FxYfPmzdja2nLt2jUNqyy8vM4L1F/Cer1eo+qKJ69zi4yMZNy4ccTExFhciEPe53bt2jXmzZvHkCFD8gxxMHOQ79ixA1tbWwC8vLwICwszLrt3Q5GDg0MRbyjSVl7nBuDg4KBFWcWW13l169YNAG9v7zJ5fu3bt0dRFLp27Ur9+vW1KrFI8jqvixcvkpWVZZF/QUHe57Zp0yZ8fX2ZMmUK6enpWpVYZHmd2xtvvIGXlxfvvPMOJ0+ezHM/Zg3yNm3asHnzZgBSU1OpWbMmAAkJCVhZ3T90XFwcSUlJ5izF5HI7N0tXkPPasGEDAQEBJV2aSeR3fvPnz+fDDz/k2LFjWpRXZLmdV2ZmJhs2bKB3794WezlwXt+zgIAATpw4waVLlwgJCdGqxCLL7dyysrK4cOECr776KiNHjqRfv35576hQtw8VUnp6uvL5558rCxYsUFq1aqVcv37duKxFixbGeTc3N3OWYRZ5nZuiKMqSJUuUJUuWaFRd0eV3Xn///bfy+++/a1Rd8eV3foqiKCtXrlQmTpyoQXVFl9t5/fnnn4q3t7fi4+OjNGzYUGnZsqVy9epVjastnIJ8z27evKkMGjRIg+qKJ7dzS01NVdq3b29cz83NTUlMTMx1P9aBgYGB5vptY21tjbe3NwaDgfr169OhQwdSUlKwt7dn69atdO3alaysLCIiIoxjuFiKh51bcnKyscvhyJEj6HQ6WrRooXGlhZPXecXExHDo0CF69+5NamoqmZmZxj8LLUV+3zeA+Ph4nJ2dadKkiYaVFk5u59W0aVNGjBjB8OHDSU5O5o033sDT01Prcgslr++ZoijodDpiY2O5ffs2rVu31rrcQsktI6tUqcKJEyfo0KEDdnZ2bNq0iREjRuS+I3P+trly5Yqyb98+ZeXKlYqiKMqRI0eUgQMHKoqiKGfPnlXGjh2rfPHFF8rly5fNWYZZ5HVuqampyujRoxV/f39Fr9drWWah5XZe8fHxSvPmzZVWrVopTz31lNKqVSuNKy2a3M7v5s2bSrdu3ZSwsDDlzz//VDIzMzWutHDy+nm8JzAwUAkPD9eivGLJ69y6du2qjB8/Xpk9e7aSlZWlZZlFkte5nTp1Spk+fboSHh6unDx5Ms/9yA1BQghh4eSGICGEsHAS5EIIYeEkyIUQwsJJkAshhIWTIBdCCAsnQS6EEBZOglwIISycBLkQQli4/wcdwFag+FO9DwAAAABJRU5ErkJggg==)
You see how steppy the histogram looks when you compare it to the analytically calculated density? Maybe we can do better. If we think about the single particle from above again
[latex]
\begin{figure}\begin{center}
\begin{tikzpicture}
\fill [red] (0.95,.15) circle (3pt);
\draw [line width=2pt] (0,0) — (11,0);
\foreach \x in {0,…,11}
{\draw [line width=2pt] (\x,0)–(\x,0.3);
}
\foreach \x in {0,…,10} {
\node () at ($(\x,-.3) + (0.5,0)$) {\x};
}
%\fill [red] (0.95,.15) circle (2pt);
\end{tikzpicture}\end{center}
\caption{A grid with one particle at 0.95}
\end{figure}
[/latex]
how about we say: It should contribute to its neighbouring bins according to its distance to the bin center. There are the two bin centers at $x_1=0.5$ and $x_2=1.5$ neighbouring the particle. We should pay attention that the total amount of density created by the particle stays the same. Let us for the time being focus on the assumption that a particle only contributes to two bins. So let us zoom in on the particle and its two neighbouring bins:
[latex]
\begin{figure}\begin{center}
\begin{tikzpicture}
\fill [red] (2.85,.15) circle (7pt);
\draw [line width=2pt] (0,0) — (6,0);
\foreach \x in {0,3,6}
{\draw [line width=2pt] (\x,0)–(\x,0.3);
}
\foreach \x in {0,3} {
\draw [line width=2pt] ($(\x,0) + (1.5,0)$) –($(\x,-.3) + (1.5,0)$);
}
\node () at ($(0.5,-.7) + (1.5,0)$) {0.5$\delta$};
\node () at ($(3,-.7) + (1.5,0)$) {1.5$\delta$};
%\fill [red] (0.95,.15) circle (2pt);
\draw[line width=2pt,||] (2.85,.5) — (4.5,.5) node[above,pos=0.5]{0.55$\delta$};
\draw[line width=2pt,||] (1.5,.5) — (2.85,.5) node[above,pos=0.5]{0.45$\delta$};
\end{tikzpicture}\end{center}
\caption{A grid with one particle at 0.95}
\end{figure}
[/latex]
The easiest way to distribute the particle density between its neighbouring bins is now to just measure its distance $d_i$ to its neighbouring bins, in the example $d_1=0.45\delta,d_2=1-\frac{d_1}{\delta}=0.55\delta$ and add $1-d_i$ to the grid point $i$. This way we make sure that the particle number is conserved (the total number added to the grid is 1, as for the histogram). In our example, the particle contributes $\rho_1=.55, \rho_2=.45$ to the grid bins 1 and 2.
Now expressing this in python and adding the bin width to calculate a proper density we write a function pics2gridpy(particles, left, right, bins). The function is called with four arguments:
particles: the Array of coordinates
left: The left edge of the grid (also the zeroth grid point)
right: The right edge of the grid (no grid point here)
bins: the number of bins (including the one at the left edge)
def pics2gridpy(particles, left, right, bins):
grid=numpy.zeros(bins,dtype=numpy.float64)
binwidth=1.*(right-left)/bins
invbinwidth=1/binwidth
leftIndexF = 0.
binPosition= 0.
for i in range(particles.shape[0]):
leftIndex=particles[i]
(leftIndexF,binPosition)=divmod((particles[i]-left),binwidth)
binPosition*=invbinwidth
leftIndex=int(leftIndexF)
grid[leftIndex % bins]+=1-binPosition
grid[(leftIndex + 1) % bins]+=binPosition
return numpy.array([numpy.arange(left,right,binwidth),grid])
Let’s try to run it on our example particles and again compare to the analytic density!
# because actually the method makes a circular grid, we need to add one bin
# to the left and right of the distribution
bins,density=pics2gridpy(x,min(x)-.05833333,max(x)+0.05833333,12)
plot(bins,density,label="density (grid)")
plot(numpy.arange(9,9.6,.01),density(numpy.arange(9,9.6,.01)),label="density (analytic)")
xlim(min(x),max(x)-0.05833)
pylab.legend(loc="upper left")
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAW0AAAEACAYAAAB4ayemAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVdX+//HXYRAcAgXEcABRUUEcw3lCc7Yc0tTyXk1TnBusX97bYHatr3ZNM6erllMqaTmXOSs55IDzrAkChiCKMiMIZ//+WHnwhMwH4cDn+XicB7D3PnuvJfJms/YadJqmaQghhDALFkVdACGEELknoS2EEGZEQlsIIcyIhLYQQpgRCW0hhDAjEtpCCGFGJLSFEMKM5BjaycnJDBs2jFatWrFz585M+/v374+LiwtjxowplAIKIYTIoMtpcM2VK1eoWrUq9+7dY8aMGXz33XeGfYGBgSQmJuLr61vY5RRCCEEu7rQ9PT2xt7cnOjqaKVOmGO3btWsXfn5+TJ8+nZSUlEIrpBBCCCXHO22A4OBgRo0ahYeHB0uWLDHal5qayvjx42nbti0jRowwPrlOZ9rSCiFEKZBtLGu5pNfrtcaNG2uJiYmZ9sXHx2tDhw412gbIS17ykpe88vnKSq57j+h0Ojp37ky5cuWIjY0FdVYA7t+/T5s2bZ76Pk3TSu3r008/LfIySP2l7lJ/86p/TqxyOmDhwoUEBQUxePBg/Pz8OH/+PDNnzsTf35+ePXvSoEED3NzcmDhxYm7zXwghRD7lGNoTJkzItM3f3x/gqV0AhRBCFB4ZXFOISntXyNJc/9Jcd5D6F2b9c9V7JF8n/qvnyN9P37t3b3799dfCuKQwA7169WL79u1FXQwhiq2sstOw/1mHtnQDFIX0X06IEiGn0M6xTbuwyA9u6SO/sIUoOGnTFkIIMyKhLYQQZkRCWwghzIiEdhYCAgLw8vIiLCzMJOd76aWXOHXqVIHOcevWLf78889cHx8SEkLLli0zbX9ct9DQUA4ePEhwcHCByiWEeHYktLPg6+vLvXv3THa++fPn07hxYwBiYmLYsGFDnt5/4cIFHjx4QPXq1XP9Hjc3N9asWZNp++O66XQ6OnTowNq1a4mLi8tTeYQQRUNCOxvlypUz2bnc3d2xsrJC0zTefvttEhMTc/1eTdN45513aNSoUZ6uqdPp8PDweOq+J+s2cuRIpk2blqdzCyGKRpF1+cuKKXuF5adX4dq1awkODqZs2bLExMQYtu/cuZPw8HCOHj3KpEmTsLe35+2336ZDhw4cP36cM2fOcPjwYapUqcLGjRuJiYlh3rx5LF26FEtLSyZNmsTMmTNxcnLiypUrpKam4uTkxLx587hz5w779u2jYsWKDBw4kNGjR9OrVy/DtS9dukSTJk2Myvm///0PW1tb/vWvf1G3bl0mT57M/v37SU5OJiYmhqpVq1KrVi1++OEHjh8/nm3dqlWrxu7du0lISKBChQp5/0cTQjw7WiEhi+kFs9qesd90r7y6e/eu1rx5c03TNC01NVV77rnntNDQUC0hIUEbNWqUpmmatm3bNu3ll1/WNE3TBg0apPn5+WlJSUna0KFDtQULFmiapmn9+/fX0tLStMDAQO3EiROapmmaj4+P9ttvv2mapmlvvPGGtmrVKk3TNO3mzZuaq6urptfrNU3TtA8//DBTuRYsWKB99dVXhq/Pnj2r9e7dW9M0TZs0aZKhbPPnz9d8fHy0+Ph4LTg4WLt48aJWs2bNbOv2WOfOnbXff/897/9oeZDT914IkfPPSbFrHjFlbOfV1q1b6dChAwDW1tY4OjoCaoWehw8fsmrVKkJCQrCzswOgbNmytG3blrJly1KnTh2io6MBsLe3p3nz5iQmJtK8eXMAypcv/7d6qgLWrFmTxo0bs3XrVsLCwnB3d89Ursd38I/99ttv2NjYANCwYUPS09MN1/D29qZChQq4u7sbNYFkVbfHHB0duXjxYt7/0YQQz1SxC+2idPPmzaf2zkhKSiI9PZ3hw4czadIkVqxYARiP8NPpdIYg/vbbbxk5ciR9+/bl559/znTs37+eMGEC8+fPZ+PGjQwcODDT9a2srEhOTjZ83b59e06fPo2maYSEhPDKK68AYGFh/O188hrBwcE59jx5/vnns90vhCh6EtpP6Nu3L/v37ycqKoqUlBRSU1NJTEzk5ZdfZseOHXz33XfcunWL5cuXA5Cenm4Iau2JCcznzJnDxIkTmTt3LufOnQNAr9cb9pcpU4aEhATDYhLdu3fn1q1bBAcHU7FixUzl8vb2NurJ0qhRI4YOHcq8efN4+eWXeemllwxl0Ov1huO0J/7c6NevX6a6JSQkGPbfuXPHcCcuhCi+it2DyKLUvHlz3n//fTp27Ejfvn1xcnJi//79TJgwgSVLlvD555+zceNGZs2aRWhoKKdOncLKyor27dtz4sQJAMOq9fHx8ZQtWxY/Pz+uXr3KjRs32LVrF+3bt6dz587MmDGDxo0b07ZtWwDGjBlDnTp1nlqutm3bsnLlSsPXFy5cYPny5aSmppKSkkLTpk1Zu3Yt+/fv58SJE1y7do169erx888/c/fuXU6dOvXUuh04cAAvLy80TaNChQrY29sX+r+xEKJgimyWv0K6rNmaPXs2b731FtbW1k/dP378eBYsWICFhQUbNmygfv36eHt7o9frOXr0KAkJCXTv3j1f1z58+DDHjh3j/fffL0gVciTfeyFyltPPiTSPFLFz585x7tw5rK2tswxsgPfff5/Vq1cDMGvWLEOXPU3TCAoKyvek6+np6aSlpRV6YAshTEPutItYjx49sLKyYuPGjYYeIVm5deuW4eHjnDlzuHHjBh4eHnz66aeZ+nHn1vXr16lbt26+3ptX8r0XImfFdhEE+cEtfeR7L0TOpHlECCFKEAltIYQwIxLaQghhRiS0hRDCjGQb2snJyQwbNoxWrVqxc+dOo30RERF8/vnnrF+/3mQLBQghhMhetqEdEhLC/PnzWbt2baZJ+ydOnIifnx99+/Zl3rx5hVrIolASVq7JyqNHj/jPf/5D37598/S+mTNnMnfu3Cz3b9++naioqIIWTwiRjWxD29PTE3t7e6Kjo5kyZYrRvuDgYJydnbG1teXChQuFWsiiUBJWrsmKtbU1Pj4+uVqt5rvvvjN8PmrUKIYNG5blsb1792bOnDmGWQeFEKaX49wjwcHB/Otf/8LDw4MlS5YAEB0dbTSjXFRUFDExMU+d7OjJFVHyO2qvqJh65RrAsHJN586dc/1e7a+Va/bt22ey8pQtWzbHYy5cuMDcuXMZNWoUAE5OTjm+p0+fPixcuJC33nqrwGUUorQICAggICAgV8fmGNq1atVi3759NG3alKSkJMqVK4ejo6PRbHKJiYlPDWwgz8tY6T4z3dI12qd5H8RhDivXBAUF8f3336NpGuHh4SxbtozDhw/j5+fHZ599xty5c3FwcGDbtm3odDo+/vhjXF1d2bRpEwsXLqR27dpGHfiHDx/O3r17OXbsGK6urkycOBEvLy/i4+OJiIhgzpw5DBw4kA8//BBPT08++ugj4uLi2LJlC7GxsZw4cYLly5djbW1NixYt+Mc//iGhLUQe+Pr6Gm5qP/vss2yPzVXvEZ1OR+fOnSlXrpzhT+o6deoQGxtLYmIinp6eBStxMXHv3j2++eYbPvnkE95++23DL6bExEQ2btzIm2++Sd++ffnkk0+oWbMmtra2XL9+nVWrVtGyZUtDk8fatWt54403WLZsGRYWFvj4+JCWloZOp6NBgwY0aNCAnj170rt3b5YsWcKDBw9wcHDA0tISLy8vo8AGtehB1apVDV9///33vPDCC0ybNo0NGzZw//592rVrR3R0NPHx8Rw6dIjz589z5coVbty4wdWrV/Hz86NOnTps2bLF6Nw6nY5Fixah1+upVKkSAJUqVWL8+PEMHjwYOzs7Jk+ejKurKy4uLoamjwULFlClShUmTZpE+fLlDe31VlZW6HQ6IiIiCuebJEQpl+2d9sKFCwkKCmLw4MH4+flx/vx5Zs6cib+/PzNmzGDq1Km4urqycOFCkxUoP3fHppKblWvi4uJyvXLN119/TceOHYHcrVzTrFmzLFeu6d27t+Hrzz77jCtXrrB27VrD3NwODg6ULVuWLl26YGFhgbu7O9HR0bRv354VK1bw008/cfv2bZydnTOdv0KFCgwaNIhly5YxZsyYLJuFnqzDunXrmDBhAgCLFy82Os7R0ZFLly7h4uKS1T+1ECKfsg3txz+UT/L39wfUnfY333xTOKUqIrlZuQZg7NixQPYr1yxevJi+ffuyevVqXn755RxXrvnvf//LSy+9xIgRIzJd/+8r16xatYoHDx7wzjvvMHXq1Kee83F5wsLCmDJlCsuXL+fy5ctZzmcwfvx4+vTpg7OzM/369XvqMU+ePzU1lcDAQLp06QKo5xpP/kJ4cnk0IYTpyOCaJ5jLyjXz5s2jadOmREZGkpSURFxcHOnp6UbPGR5fc+vWrdja2mJtbU1ISAiPHj0iMTExU3jXr18fV1dX/P39Dc1dZcqUMRybkJBgVIc+ffowbdo0rl27xubNm41+qaSkpNCwYcN8fheEENmR0H7Ck6u7fPrpp4aVa+zt7VmyZAnz5s3Dz8+Ptm3bGlau+e233wgODubEiROcOHHCsHLNJ598wu3btzOtXKPX6+ncuTNLly41Wkh3zJgxhrvWv2vbti1BQUGGr0eMGMHw4cNZtmwZPj4+nDp1iuPHjxMREcG2bdu4evUqwcHB7Nq1i549e3LkyBGGDBlC/fr1uXnzJikpKezYsYOgoCAuX75sOO+4cePo2rWr4WsXFxeqV6/O5MmTSUpK4tixY5w4cYLo6GgmT55M/fr16dWrF1FRUbi5uQGq/b9evXqm/tYIIf4iU7MWE3lZuaawrFixgt69ez+13Tu3/P39sba25tVXX820T773QuQsp58TWSOyiD1uPsntyjWP29VN6dq1azx69IjIyMgCBXZSUhLu7u60bt3ahKUTQjxJ7rSLWH5WrnF1dTVpGcaMGcOlS5fYvHkzlStXzvd5cloFR773QuRMVq4RxYZ874XImaxcI4QQJYiEthBCmJEiexD598EmQgghcvbM77T/Pq+GKF3k+y9EwTzzB5FCCCGyJg8ihRCiBJHQFkIIMyKhLYQQZkRCWwghzIiEthBCmBEJbSGEMCMS2kIIYUYktIUQopgIDs75GAltIYQoBm7fhiwWrzIioS2EEEUsOhq6dYObN3M+VkJbCCGKUHw89OoFly6Bl1fOx0toCyFEEXn4EPr1gxMnoGZN2L075/dIaAshRBFIS4PXXoP9+6FKFdizB6pVy/l92Yb2o0ePGDlyJJ6enkyfPj3T/v79++Pi4sKYMWPyXXAhhCht9Hp4803YsgUqVlR32HXq5O692U7N+vvvv+Pt7Y2maXh5eXHo0CFq1aoFQGBgIImJifj6+j79xDI1qxBCZKJp8M47MG8elCsHe/dC69YZ+ws0NWubNm2ws7PD3t4eb29vo9XCd+3ahZ+fH9OnTyclJcUEVRFCiJLvs89UYJcpo+60nwzs3MjVcmORkZF4e3tT7YkGl48//pgPPviA8ePH4+/vz4gRI5763mnTphk+9/X1zfLOXAghSrpvvlGhbWEBP/wAXbuq7QEBAQQEBOTqHLlauWbBggWMHj3a6E77sYSEBMaOHcuaNWuMTyzNI0IIYbByJTy+t12xAt544+nHFXjlmk2bNjFgwABsbGyIiooiLi7O6IT379+nTZs2eSu9EEKUIps3qwePAF9/nXVgb7qyKcdzZds8smDBAmbPno2TkxMPHz5k9OjRHDt2DH9/f3r27EmDBg1wc3Nj4sSJea2DEEKUCnv3wpAhqsfI1KnqIeTT7A7azZANQ3I8nyzsK4QQheTYMTWfSGIiTJqk2rT/ikYjh8MO0211N5LTkmGa2iYL+wohxDN04QL07KkCe9gwmDv36YF98vZJevv3JjktmRFNnt6h40lypy2EECZ24wa0bw+RkdC3L2zYAFZPaYw+efskXb7vQmxKLAO9BrJuwDqsLNWBWWWnhLYQQphQeDi0awchIdC5M2zfDra2mY97MrD71+/P+oHrsba0LnjvESGEELnzeIrVkBBo0UINnnlaYAeGBxoC+xXPVwyBnRsS2kIIYQLx8aoN+/JlaNAAfv0Vnnsu83GB4YF0Xd3VENjrBqzLdWCDhLYQQhTYw4eq7TowENzd1QRQjo6ZjytoYIOEthBCFMijRzB4MBw4AM8/r6ZYrVo183Enwk8YAnuA54B8BTZIaAshRL7p9TByJGzbBpUqqcCuXTvzcQEhAbz4/YuGwP5hwA/5CmyQ0BZCiHzRNHj7bVizBsqXhx07wNs783G/XP+FHmt6kJCawGverxUosEFCWwgh8uXTT2HBAjXF6tat0LJl5mN+uPAD/df3JyU9hbEvjGV1/9UFCmyQ0BZCiDz7+muYPl1NsbpuHbz4YuZjFp9czNBNQ0nTpzGl7RQW9V6EpYVlga8toS2EEHmwfDlMnpzxef/+mY/58vCXjNs+Dg2NGS/OYGaXmYZBMwWVq0UQhBBCwMaNMHq0+nzuXBg+3Hi/pml8uO9DZh6ZiQ4dC3stZFzzcSYtg4S2EELkwp498PrrqsfItGnqIeSTHqU/YvTPo1l1bhWWOktW9VvF0EZDTV4OmXtECCFycPSommI1KUmF9ddfG8/YF5cSx8AfB7IneA/lrMvx48Af6V23d76ulVN2SmgLIUQ2zp+Hjh0hJkatOLNsmXoA+djt+Nv09u/N2cizOJd3Zvvr2/Gp6pPv60loCyFEPt24oWbsu3NHPXD88UfjKVYv371Mz7U9CYsNw8PBg53/2EmtSrUKdE0JbSGEyIfwcGjbFkJDVZe+X34xnrHvYOhB+q7rS8zDGFpXb82217bhVM6pwNeVqVmFECKP7t2Drl1VYLdsmXmKVf8L/nRd3ZWYhzH0q9+PvcP2miSwc0NCWwghnhAXp6ZYvXJFDUv/9VeoUEHt02t6Ptr3EUM3DSU1PZUJzSew4dUNlLMu98zKJ13+hBDiL8nJaorVkyehVi01xaqDg9qXkJrAsM3D2Hx1MxY6C77p8Q0Tmk8w2aCZ3JLQFkIIMqZYDQgAFxfYu1d9BLgVe4s+6/pwNvIs9jb2/Pjqj3Sr3a1IyimhLYQo9fR6GDECfv5Z3Vnv2aMWMwA49ucx+q3rx53EO3g4ePDzaz9Tz6lekZVV2rSFEKWapsGkSbB2rWq73rFDLRcGsPrcanxX+nIn8Q4vur/IsVHHijSwIYfQfvToESNHjsTT05Pp06cb7YuIiODzzz9n/fr1hIWFFWohhRCisHzyCSxaBDY2aorVFi0gNT2VSb9OYtiWYaSkpzDOZxw7hu7AoaxDURcXtGwcOXJEi42N1WJiYrSqVatqQUFBhn2vvPKKdufOHS05OVl77733Mr0X0HI4vRBCFKmvvtI00DRLS03bskVtC48L19osa6MxDc36P9ba4sDFz7RMOWVntm3abdq0MXzu7e2NjY2N4evg4GCcnZ0BuHDhQiH8OhFCCNPT62H7djVL3/79atvy5arXyKHQQwzaMIjIhEiqPVeNjYM20rL6U1Y3KEK5ehAZGRmJt7c31apVAyA6OhqLJwbfR0VFERMTQ8WKFTO9d9q0aYbPfX198fX1LViJhRAiHxISYMUKmDdPDU8H1YY9Zw78858a3xybx/t73idNn4ZvTV/WD1yPc3nnZ1K2gIAAAgICcnVsroaxL1iwgNGjRxvdaTdt2pQzZ84AULduXa5fv258YhnGLoQoBkJC1LJg330HsbFqW82a8NZbalFey7IJ+P3sxw8XfwDgvdbvMbPLTKwsiqZzXU7ZmWOpNm3axIABA7CxsSEqKgpbW1vs7OyoU6cOsbGxWFlZ4enpadpSCyFEAWga/P67mkJ182bVJALQvj28845qCrG0hHOR5xi0ehDXo69T3ro8y/suZ1CDQUVb+Bxke6e9YMECZs+ejZOTEw8fPmT06NEcO3YMf39/bty4wfz583F1dWXw4MFUr17d+MRypy2EeMZSU+Gnn1R79cmTapu1tRo088478MILapumaSw+uZh3d71LSnoK3s7erB+4Hq/KXkVX+L/ILH9CiBLv3j1YuhQWLoTbt9U2R0cYN069qlbNODbmYQyjfx7NhssbABjdbDRze8x9pvOHZKfAzSNCCFFcXboE33wDq1fDw4dqW4MG6q566FAoW9b4+MDwQAZvGMzNmJs8V+Y5lr68lCHeQ559wQtAQlsIYVb0eti1SzWB7N6dsb1XLxXWXboYLwUGana+OUfn8OG+D3mkf0Qzl2asH7ieOg51nm3hTUBCWwhhFhIT1R31N9/A1atqW7lyakX0t9+GelmMLg+LDWP4luEEhAQAMKnFJGZ1nYWNlc3T31DMSWgLIYq1W7dUW/XSpfDggdpWvbqaL2TUqIypU5/G/4I/47ePJzYlFufyzizvszzfC+4WFxLaQohi6fhx1WVvwwZIT1fbWrVSTSCvvKJ6hWTlQfIDxv86nnUX1wHQt15fvn35WyqXr/wMSl64JLSFEMVGWhps3Kjaq48dU9ssLTO67LVqlfM59t/cz/Atw/kz7k/KW5dnbo+5vNn0zWe+WEFhkdAWQhS5Bw/g22/VyMVbt9S2SpXAzw8mTIAaNXI+R0JqAlP2TGHRyUUAtKzWkjWvrDHLh43ZkdAWQhSZa9fUg8VVqyApSW2rV089WBw2DMqXz9159gXvY9TPowiJCcHKwopPOnzCh+0/LLKh6IWp5NVICFGsaZpaymvuXLVo7mNdu8K770L37mCRy+VZ4lPi+X97/h9LTi0BoOnzTVnZbyWNqjQqhJIXDxLaQohnIjlZrQ4zd64aFANgawv//KeavMnbO2/n2xu8lze3vUlYbBjWFtZM7TiVKW2nYG2ZzRPKEkBCWwhR6H74QQXzvXvqaxcX1Vbt5weV89ihIzopmvf3vM/KsysBeMHlBVb0XUHDKg1NW+hiSkJbCFFo0tJgyhQ1ZzWoCZvefRdefRXKlMnbuTRNY/X51by3+z3uJd3DxtKGqR2n8kHbD0pk23VWSk9NhRDP1N27MGSIWh3Gyko1i4wfn3mIeW78Ef0H47aPY9/NfQB0qtmJxS8tpq5jXROXuviT0BZCmNzp09C/P4SFQZUqaoBMu3Z5P09qeiqzjsxi+sHppKSn4FjWkdndZjOs8bAS0+86ryS0hRAmtXq1aqt++FANhtmwAf5aqTBP9gbvZdKOSVy9pyYaGd54OF91+wqnck4mLrF5kdAWQpjEo0fw/vtqDUaA0aNh/nywyeO8TGGxYby3+z3DfNd1Hevyv97/o7N7ZxOX2DxJaAshCiwqSj1cPHhQzQmyYIG6286LlLQUZh+dzReHviDpURLlrMvxSYdPeLfVu2Y7I19hkNAWQhRIYKCawOnPP1VXvo0boXXrvJ1jxx87eGvnW9y4r5ZJH9RgEF91/Yoa9rkYv17KSGgLIfJtxQq1nFdKCrRpo9qvXVxy//7Ldy/z3u732HljJwBelb2Y33O+NIVkQ0JbCJFnqamqv/UiNTcT48apLn257Xt9L+kenx74lCWnlpCupWNvY8/UjlOZ1GJSiR/RWFAS2kKIPImMhIED4cgRFdKLFsGbb+buvSlpKcw/MZ/PD35ObEosljpLxvuMZ5rvtBIx1/WzIKEthMi1Y8dgwAC14nm1arBpE7RokfP7NE3jp8s/8eG+Dwl6EARAjzo9mN1tNl6VvQq51CWLhLYQIle++07NF5KaCu3bw08/qYEzOdkXvI9/7fsXJ2+fBFS79exus+lRp0chl7hkktAWQmQrJUXNb71EzX7KpEkwe3b2y30BnIk4w7/2/YvdQWrJdJcKLkzzncbIpiNL1VwhppbjrLUHDx6kS5cuT93Xv39/XFxcGDNmjMkLJoQoerdvQ6dOKrBtbGDlSjV4JrvADrofxOsbX6fZ0mbsDtqNnY0d/9f5/7jx1g38XvCTwC4gnaZpWk4HtW7dmqNHjxptCwwMJDExEV9f36ef+K95AXJxeiFEMXTkiHrgGBmplvvatAl8fLI+PjQmlOkHp7Py7ErStXTKWJZhYvOJfNj+QxzLOT67gpu5nLIzV+tD2NraZtq2a9cu/Pz8mD59OikpKQUoohCiONE0WLxY3WFHRoKvL5w6lXVgh8eFM377eDzme7DszDI0NN5o8gbXJ15ndvfZEtgmlqs77U6dOnHgwIFM21NTUxk/fjxt27ZlxIgRxif+67fFp59+atjm6+ub5Z25EKLoPXwIEyfCsmXq63ffhf/+V02t+neRCZHMPDyTxScXk5Kegg4drzd8nakdp5bKKVMLIiAggICAAAA+++wzIOs77QKFNkBCQgJjx45lzZo1xieW5hEhzMqff6rufCdOqGXAvvsOhg7NfFx4XDizfp/F0lNLSU5LBuBVr1eZ5jtNuu+ZQE7ZmecnArGxsdjb26NpGjqdjvv379OmTZuClVIIUaQOHVLt11FR4OYGmzdD06bGx4TGhDLz8EyWn11OanoqAP3q9+Mz389K9EK6xU2OoX3u3DmuX7/OxYsXSU9P58svv8Tf35+ePXvSoEED3NzcmDhx4rMoqxDCxDQNFi5UzSBpafDii7BuHTg9MWX1H9F/MOPwDFafX02aPg0dOgY1GMSH7T6k8fONi67wpVSumkfydWJpHhGiWEtOVnOGrFqlvn7/fZgxI6P9+nTEaWb9PosfL/2IXtNjqbPk9Yav8+92/8azsmfRFbyEM3nziBDC/IWFqelUT52CsmVh+XK1nqOmaewN3seXR75kb/BeAKwsrBjRZAT/bvdvajvULuKSCwltIUqZgAC1YMG9e+DurtqvGzRMY93FDfz3yH85E3kGgAplKjDmhTG80+odqttVL9pCCwNpHhGilNA0NZrxvfcgPR26dYMlq+LYfHMZ807MIyQmBADn8s683fJtxvmMo1LZSkVb6FJImkeEECQlwZgx8Lhn7th/36RM+3k0WrGM+NR4AOo41OH91u8zrPEwylqXLcJ4DYsRAAAZEUlEQVTSiuzInbYQJVxIiGq/PnNGw7buERqN+ZqTCVvQa3oAOrp15N1W7/JS3ZewtLAs2sKKHLNTQluIEmzvXhj8jyTuV/2BMu0XkOpwFgBrC2uGeA/hnVbv0MylWRGXUjxJQluIUkjT4KOvgpmxbxE0WQ5lHwDgVM6JMS+MYULzCbg8l4fFHMUzI6EtRCmSrk/nx9O7eH/9Im6X/xV06uevedUWTGoxkVcbvIqtVeYJ4ETxIQ8ihSihkpLgyhW4cAGOXrzNvvvLCHH4jvTnwqACkGZDp8pD+HLABJpXa17UxRUmIqEtRDGXng43bqhwvnhRfbxwAf4ISodau+GFpVDvZ3guHQBdTC3c7/ux6p2RtGsqi+WWNNI8IkQxoWkQEZERyo9D+vJlNWWqgX0oNFkJTVdAxVAALLCipX1fxvqM4fXWL2Jlmaup8kUxJM0jQhRDsbFw6VLmgL5//+nHV6/5EKd2m3ngvpwwy31oqB/omhVrMrrZaEY2HcnzFZ5/hjUQRUVCW4hClJoKV69mbtoIC3v68RUrQsOG6uXtrWHtdpLjKavYcH0tfz6MAcDG0oYBXgMY2WQkndw7YaGTu+rSRJpHhDABvR5CQ43vmi9cgGvX1JSnf2djA15e4O2dEdING0LVqnArLoy159fy/fnvuXrvquE9PlV9GNlkJEO8h8jw8hJMuvwJYWJ6vQrj48fV68wZ1dSRkJD5WJ0Oatd+8u5ZfaxTx3gJr/iUeDZe2cj3574nICTA0PzhXN6Z171f540mb8jc1aWEtGkLUUB372YE9PHjajmu2NjMx1WpYnzX7O2t7qbLl3/6eVPSUthxYwf+F/z5+frPPExTTxttLG3oV78fwxoPo1vtblhZyI+pyCD/G4R4wsOHcPYsHDuWEdI3b2Y+rlo1aNlSvXx8VEhXzkXvunR9OgEhAfhf8GfjlY3EpmSkf3vX9gxrPIyBXgOpaFvRhLUSJYmEtii1NE31f37yLvrsWXj0yPi4cuVUMLdsCa1aqY/VquX+OnpNz5GwI/x46Uc2XNlAZEKkYV+T55vwuvfrDPYejKu9q4lqJkoyCW1Raty/r5o2Ht9FnziRuYudTgcNGmTcRbdsqb62yuNPil7Tc/TWUUNQ346/bdhXq1ItXvd+ndcaviarl4s8k9AWJVJqKpw7Z3wX/ccfmY+rUiUjnFu1UnfUdnb5u2a6Pp3DYYfZdGUTG69sJDw+3LDPzd6NQQ0GMajBIF5wecHwsEmIvJLQFmZP09Sc0Y/D+dgx1aMjJcX4OFtbeOEF47toV1d1d51fKWkp7L+5n01XNrH12lbuJt017HO1d2WQlwpqn6o+EtTCJCS0hdmJjVVNG0/eRd+9m/m4evWMA7pRI7C2Lvj141Li2HljJ1uvbeWX678QlxJn2FfHoQ4DPAfQv35/WlRrIUEtTE5CWxR7d+/Cvn1qQv/ff1cz2/2dk5NxQDdvDpVMOP7kz7g/2XZtG1uvbeXAzQM80mc8rWzo3JABngN4xfMVvJ29JahFoZLQFsVOcjIcPgx79qjX2bPG+8uUgaZNjUO6Vq2CNXP8nV7TczriNL9c/4Vfrv/CqYhThn0WOgvau7anT70+9K3XFw9HD9NdWIgc5BjaBw8e5D//+Q979+412h4REcGyZcvw8PCgdevWuLpKdyWRP3q9emj4OKQPHzae1c7GBtq3h65doWNHaNJEbTO1+JR49gTvYfv17Wz/Yzt3Eu8Y9pWzLkf32t3pU68PvT16U7m8THkqikaOod2hQwcSExMzbZ84cSL/+9//sLOz4+OPP+arr74qlAKKkiksLCOk9+2De/eM9zdtqkK6a1do2xbKFsLi4JqmcfnuZXbc2MHOGzs5GHrQqNmjhl0NXqr7Er09etPZvbOsUC6KhVw1j9jaZl6eKDg4GGdnZwAuXLiQ5XvDwtQTelG6xcXBgQMZQX39uvH+GjUyQvrFF3M3ujA/Yh/Gsu/mPnbe2MnOGzu5FXfLsM9CZ0GbGm14yeMlXqr7krRPi2IpX23a0dHRWFhkTAcZFRVFTEwMFStmHnrbsOE0hg8HBwfw9fXF19c334UV5uPRI9XD43FIHz+uVmB57LnnoFOnjKCuW9e0bdKPpenTCAwPZHfQbvYE7+HYn8dI1zIK4lzemR51etCjdg+61e6GYzlH0xdCiBwEBAQQEBCQq2NzNctfp06dOHDggNG2pk2bcubMGQDq1q3L9b/dOmXcoWi4uKgn/14y+KvE0jQ1893jkA4IgPj4jP2WlmrwyuOQbtEi76MMc1cOjRv3b7A3eC97gvew/+Z+o/k9LHWWtKreip51etLToydNnm8i81GLYsXks/zFxcVhZ2dHnTp1iI2NxcrKCk9PzyyP9/VVP8AdO6of5iZN8npFUVxFRan26D171C/lW7eM99erlxHSvr75H2mYk9vxt9kXvI99N/ex/+Z+oyYPAA8HD7rW6kq32t3wremLva194RREiGcgxzvtc+fO0atXL3bt2kV6ejpffvkl/v7+3Lhxg/nz5+Pq6srgwYOpXr268Yn/+m2RlKQxYADs2KFW5di5U3XREuYnORkOHcoI6b93xatcGbp0USHdpYtqpy4MUYlR/BbyGwEhAewP2W+0UACAY1lHOrl3olutbnSt3ZWaFWsWTkGEKATFYhGElBR47TXYvBkqVIBfflF33qJ40+tVMD/ZFe/JoeG2thld8bp2VSMOLQqhpeFu4l1+C/2NAzcPEBAawOW7l432VyhTgQ5uHXjR/UU6u3emUZVG0uQhzFaxCG1QSy4NHw7+/qr71ubN0L17YVxZFERysvqraMMGFdR/74rXrFnG3XS7diq4Te1W7C0Ohh7kYOhBDoUd4so94yGQZa3K0ta1Lb5uvvjW9KVFtRZYW5pgfLoQxUCxCW1QvQfGjoXvvlOj2n78Efr2LYyri7xITVUBvW4dbNlivGyWq6txVzwnJ9NeW6/puXL3CkduHeFw2GEOhh4kNDbU6JiyVmVpU6MNvjV96VSzE82rNaeMZRnTFkSIYqJYhbb6Gt55B+bNUz0KVq9WTSfi2UpLUw+I162DTZvgwYOMfT4+MHgw9OkDHh6m7YqX/CiZwNuBHAk7wpFbR/j91u88ePjA6Bh7G3vaubajg1sHOrh1oJlLMwlpUWoUuzUidTqYO1etmzdjBgwdqv4kHznyWZek9NHr4cgRFdQbNqjeH481bAhDhsCgQWrRWVPQNI2w2DCO/nmUo7eOcvTPo5yJPEOa3nh58mrPVaOdazva1mhLB7cOeDt7Y2lhaZpCCFHCFOlq7F98AR9/rD6fPx8mTiyMkpRumgaBgSqof/wRwjPm5aduXRXUgwebpg99YmoiJ2+f5ET4CY6FH+PoraNEJEQYHaNDR8MqDWlXox1tXdvStkZbXO1dZeShEH8pds0jfzd3Lrz7rvp85kyYMqUwSlO6aBqcP6+Cev1644Vp3dwygrpJk/w3faTr07l89zInwk9wPPw4x8OPczHqInpNb3RcJdtKtKreitbVW9O6RmtaVGuBnU0hddgWogQo9qENsHSpekCpaTB1KkybVjhDmku6q1dVSK9bpz5/zMVFhfTgwaqPfF7/bTVNI+hBEIHhgQTeVq/TEadJepRkdJylzpJGVRrRslpLWlZvSavqrajrWFe63wmRB2YR2gBr1qgugXo9vPcezJolwZ0bN29mBPW5cxnbnZzg1VdVULdrpx765sbjgD51+xSnItTrdMRpYh7GZDq2ZsWaNK/a3BDSzVyaUc66nIlqJkTpZDahDbBxo+pJ8uiRuvNeuLBwBmuYuz//hJ9+UkF94kTG9ooV4ZVXVFB37pzz3B7p+nSuRV/jTMQZzkSe4XTEaU5HnDaaq+OxKuWr0Lxac5pXVS+fqj4yp7QQhcCsQhvg119V8KSkwLBhsGxZ4UwsZG6iolSPj3Xr1FDyx8qXh379VFB365b14gBJj5K4GHWRs5FnORt5ltMRpzl/5zzJacmZjn2+wvO84PKCelVVH6s+V1UeFgrxDJhdaIOahKhPH0hKUn/ir1mjBuOUNvfvq5Gj69bB/v2q6QjUKMTevdUDxV69oNwTLRKapnE7/jbn7pzjXOQ5zt05x9nIs/xx/49MDwkB3OzdaOrSlKbPq9cLVVVACyGKhlmGNqj+xL16qcnzX3pJNQcUxpDp4iYuDrZtU0G9e7dqKgK1inj37iqo+/RR81EnpCZw+e5lzt85b/T6+2AVUA8JPSt70uT5JjSu0phmLs1o8nwTHMo6POMaCiGyY7ahDXDqlPqT//59Nd/Fli2qOaCkSUqC7dtVUP/6a8b6iJaWqm164OBUGnS4RtjDi1yMusjFuxe5cOcCN2NuPvV8DmUdaOjckMZVGquQfr4xXpW9sLUqBb/1hDBzZh3aABcuqHkv7txRvSC2by+8eZmftevXVT/177+HxETAMgWcrlO//SXcW1zGosolguIv80f0H0arrTxmbWFNfaf6NKrSiEZVGtHQuSGNqjSS9mchzJjZhzaocHvxRdVrwscHdu1Sy5eZI02DnfsTmPHtVQ5duQKVL0PlK9hWv0JqhSD0ZA5nHTpqVapFwyoNaejcEG9nb7ydvfFw8JDZ7YQoYUpEaAOEhKimgps31TwZe/ZAlSomOXWh0DSN8Phwrt27xtV7V7kUdZVDl69xLfoqj8rdeup7dOio7VAbr8peNKjcgAaVG+BV2Yv6TvVlJXAhSokSE9qg5s3o0kWN9qtXT62e8rcFc565+JR4rkdf51r0NcPHa/fU54mPEp/+pnRrnC3r0qq2F02qeeJZ2RNPJ0/qOtaVcBailCtRoQ2qv3LXrmpuDXd31T3Q3d2kl8gk6VESQfeDuB59nT/u/8Ef0X+oj/f/IDIhMsv32eodSb1dH31UfbhXH9dy9Rk7sB6T/ulOhXLS+VwIkVmJC21QvUl69FCz11WrpoK7Xr2CnTMuJY6g+0HcuH/D8Ap6oL4Ojw/P8n02ljZ4OHpQz7EeHg51IbouR7bW49DWupDkCKiuepMnq1828nxQCJGdEhnakNF/+9AhcHZWbdyNGmV9vF7Tczv+NsEPggl+EEzQ/SCCHvz1uh9EdHJ0lu+1trDGvZI7Hg4e1HWsi4eDBx6O6vPqdtVJT7Pgp59gzhzVTRHUyMR//EMt+ODtbeLKCyFKrBIb2qC6yfXvrwK7YiWNDb88oGLNm4TEhHAz5iY3H9wkOEaFdEhMCKnpqVmeq6xVWWpVqkUdhzrUcahD7Uq1DZ/XsK+BlUXm5owHD+Dbb9UqPI/nqXZyggkTYNy44v2gVAhRPJWY0NY0jZiHMYTEhBheobGhBN0P4eC5m8RZhIBtXLbnqFK+CrUq1cK9kju1K9VWL4fa1KpUC5cKLrnu23zjBnzzDaxY8Vf/asDTUzWBDB2qFi4WQoj8KHbLjWUlXZ9OZEIkobGhhMaEEhYbRmhsxsfQmFDiU+Of/ubHc2+kVEAX606r+jVpXscd90ru1KpUSwV1RXfKl8n/cEpNg8OHVRPI1q3qa1Dt1JMnq5GbMiOhEKKwPZM7bU3TuJ98n1txt7gVe8vwMSwuTH2MDSM8PjzT2oF/V6FMBWpWrEnNijVxs3czfHSv5I7rc+588JYDq1bqsLFR07z27l3wejx6pGbXmzMHTp5U28qUUXfU776r+owLIYSpmKR5ZOnSpTg6qp4QAwYMMNrXv39/jh07Rp8+fViyZEmmC9edX5dbsbeeOgXo3zmXd8bN3g23im642rmqj/auuNmrjw5lHbJtwtDrYdIkWLRITefq769mCcyPmJiM9uo//1TbHB1h/Hj1ev75/J1XCCGyU+DQPnnyJKtXr+abb77h008/5a233jIEeGBgIImJifj6+mZ5YaapD3Y2dtSwq0F1u+rUsKuBq72r0au6XXVsrLKYDDoPNA0++AC++ko1V6xYoeblzq3gYNVevWxZRnt1/frqrvqf/5T2aiFE4Spwm3ZAQAAtWrQAoHbt2uzevZvXXnsNgF27dvH999/zz3/+kw8++ACbp8zAf3HcRWrY13hmi7nqdPDf/6rZAD/7TC1hlpSkVsLJiqbB77+rJpAtWzLmrX7xRdVe3aOHtFcLIYqHHEM7LCyMen+NXLGzs+PAgQOG0P7444/54IMPGD9+PP7+/owYMSLT+39a9JPhc19f36felZuaTqcWBy5fXt11jxungnvyZOPj0tJU2/ecORnLdllbqzvqd9+Fxo0LvahCCEFAQAABAQG5OjbH5pHFixdjYWGBn58fixcvJikpicl/S7+EhATGjh3LmjVrMk78DPpp58bChTBxovp8+nT46CM1MOe771R7dViY2ufgoMJ9wgS1erkQQhSFAjePtGvXjuXLlwPqrrtfv37ExcVhZ2eHpmnodDru379PmzZtTFhs05kwQS3HNWoUfPKJ6rZ35AgkJKj9deuqu+phw4yX7RJCiOIoV71HvvjiC6r8NbyvRYsWzJw5E39/f3r06EGDBg1wc3Nj4sSJWDzR8Ftc7rQfW7dODStP/2u66k6dVHNJr17SXi2EKD5KzIhIU9i7V618M2wYNG1a1KURQojMJLSFEMKM5JSd0jAghBBmREJbCCHMiIS2EEKYEQltIYQwIxLaQghhRiS0hRDCjEhoCyGEGZHQFkIIMyKhLYQQZkRCWwghzIiEthBCmBEJbSGEMCMS2kIIYUYktIUQwoxIaAshhBmR0BZCCDMioS2EEGZEQlsIIcyIhLYQQpgRCW0hhDAjEtpCCGFGJLSFEMKM5BjaS5cuZePGjWzcuNFoe0REBJ9//jnr168nLCys0ApozgICAoq6CEWqNNe/NNcdpP6FWf9sQ/vkyZNcunSJAQMGcP78eaKjow37Jk6ciJ+fH3379mXevHmFVkBzJv9xA4q6CEWmNNcdpP5FFtoBAQG0aNECgNq1a7N7927DvuDgYJydnbG1teXChQuFVkAhhBBP0LIxadIkbdu2bZqmadrmzZu1t956S9M0Tbt3757WrFkzw3FNmjTRHjx4YPReQF7ykpe85JXPV1ayvdP28vIiIiICgMjISNzc3ABwdHREr9cbjktMTKRixYpG71W5LYQQIq8OHDiQ5T6r7N7Yrl07li9fDkBYWBj9+vUjLi4OOzs76tSpQ2xsLFZWVnh6ej71/RLcQghhWjoth2T94osvqFKlCgAtWrRg5syZ+Pv7c+PGDebPn4+rqyuDBw+mevXqz6TAQghRmuUY2kIIIYoPGVwjhBBmJNs27dxKS0vjl19+4eHDh8TExDB27FjDvoiICJYtW4aHhwetW7fG1dXVFJcsVrKrf0xMDJ9//jlNmzZl6NChRVjKwpNd/X/++WdmzZqFpaUlixYtyvL5hznLrv4//PADK1aswMnJidWrV2NpaVmEJTW97Or+2JtvvsmwYcPo2LFjEZSwcOVU//79+3Ps2DH69OnDkiVLTHJNk9xp//rrr0RGRjJkyBDOnj1r6HECpWMQTnb1r1ixIh4eHqSlpRVhCQtXdvW3srLi4MGD9OzZkxkzZhRhKQtPdvV3dXVl9+7dWFtbc/v27SIsZeHIru6gfmknJiYWUekKX3b1DwwM5O233yYiIsJkgQ0mCu1Dhw5hbW0NQOPGjdm8ebNhX2kYhJNd/QFsbW2LoljPTHb179mzJwBt2rQpsf8O2dW/bdu2aJpG9+7dqVGjRlEVsdBkV/eQkBDS09NL5F9Xj2VX/127duHn58f06dNJSUkx2TVNEtotWrQwjJZMSEigcuXKAERHR2NhkXGJqKgoYmJiTHHJYiWr+pcWuan/jh07+Pjjj5910Z6JnOq/aNEiPvroI86fP18UxStUWdU9LS2NHTt20K9fvxLd9Te77/3HH3/MxYsXCQ0Nxd/f33QXzW5EZG6lpKRos2bN0hYvXqz5+Phod+/eNRot+ZiHh4cpLlfsZFd/TdO0lStXaitXriyi0hW+nOp/7do1bevWrUVUusKXU/01TdN+/PFHbcqUKUVQusKVVd3379+vtWnTRvP19dVq1qypNW3aVAsPDy/i0ppebr738fHx2tChQ012Tctp06ZNK2jwW1pa0qZNG/R6PTVq1KBdu3bExcVhY2PDvn376N69O+np6Rw8eJDXXnvNBL9qipen1T82NtbQHHD27Fl0Oh1NmjQp4pIWjuzqHxERwenTp+nXrx8JCQmkpaUZ/pwsKXL6/gPcu3ePKlWqUK9evSIsqellVff69evz5ptv8sYbbxAbG8u4cePw9vYu6uKaXHbfe03T0Ol0REZGkpycTPPmzU1yTZM0j4SHh3PixAlCQ0MZNWoU586dMzxFnTFjBlOnTmXx4sUsXLjQFJcrdp5W/3HjxgFqiP/x48c5efIkSUlJRVzSwpFV/aOjo+nRowczZszAx8eHTp06Ua5cuaIursllVf+EhAR69erFli1b0Ov19OrVq6iLanLZ/d9/klZCm0iyq3/Pnj1577332LJly1N71eSXDK4RQggzIoNrhBDCjEhoCyGEGZHQFkIIMyKhLYQQZkRCWwghzIiEthBCmBEJbSGEMCMS2kIIYUb+P/o1IKwihBIDAAAAAElFTkSuQmCC)
We see, this is much better. But one thing is notable: The spike for the bin at 9.0, where does it originate? In this region, the density of perticles per bin is smaller than one. Because of that you have to expect aliasing effects that can produce these spikes. In fact, normally you try to have at least a few tens of particles in the bins that matter to you because you want to have a good approximation of the real density. In our case we needed very few because we had a smooth distribution. In most cases however you sum random particles and therefore you expect some noisiness due to the random fluctuations of your distribution.
Now that we have the method let’s get an estimate of its performance for very many particles:
x=random.random(100000)
%timeit bins,density=pics2gridpy(x,0,1,100)
gives a time of 844 msec on my machine. The histogram however only takes 8.5msec. There has to be a way of optimizing this. Next post we will look at how fast we can make it using cython!