position_distribution.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "from numpy import linalg as lin\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "from tqdm import tqdm\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#####  Function to compute time for wall collision ######\n",
    "def wall_time(coord, velcomp, rad):  \n",
    "    if velcomp > 0.0:\n",
    "        del_t = (1.0 - rad - coord) / velcomp\n",
    "    elif velcomp < 0.0:\n",
    "        del_t = (coord - rad) / abs(velcomp)\n",
    "    else:\n",
    "        del_t = float('inf')\n",
    "    return del_t\n",
    "\n",
    "\n",
    "# Function to calculate time it takes for a pair of particles to collide. \n",
    "# pos1 and pos2 are positions of particles 1 and 2, vel1 and vel2 are their \n",
    "# velocities. rad is the radius of the particles. \n",
    "def pair_time(pos1, vel1, pos2, vel2, rad):  \n",
    "    rel_pos = pos2 - pos1\n",
    "    rel_vel = vel2 - vel1\n",
    "    rel_dist_squar = np.dot(rel_pos,rel_pos)\n",
    "    rel_speed_squar = np.dot(rel_vel,rel_vel)\n",
    "    scal_prod = np.dot(rel_pos,rel_vel)\n",
    "    a = scal_prod ** 2 - rel_speed_squar * ( rel_dist_squar - 4.0 * rad **2)\n",
    "    if a > 0.0 and scal_prod < 0.0: ## Conditions for collision.\n",
    "        del_t = - (scal_prod + np.sqrt(a)) / rel_speed_squar ## Collision time.\n",
    "    else:\n",
    "        del_t = float('inf')\n",
    "    return del_t\n",
    "\n",
    "\n",
    "# Function to generate an image of the disk configuration in the box.\n",
    "def create_picture(positions):\n",
    "    plt.cla()\n",
    "    plt.gca().set_aspect('equal')\n",
    "    plt.axis([0, L, 0, L])\n",
    "    plt.setp(plt.gca(), xticks=[0, L], yticks=[0, L])\n",
    "    for x,y in positions:\n",
    "        atom = plt.Circle((x, y), Ratom, fc=\"r\")\n",
    "        plt.gca().add_patch(atom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def x_component(positions):\n",
    "    A = []\n",
    "    for i in range(len(positions)):\n",
    "        a = positions[i][0] \n",
    "        A.append(a)\n",
    "    return np.array(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This function calculates the coordinates of intersection points of a regular grid pattern. \n",
    "#These coordinates are used to fix the initial position-configuration of the disks. \n",
    "def positions_f(N_atoms):\n",
    "    global Ratom\n",
    "    \n",
    "    density = 0.2 # Fraction of the area of the box covered by the disks\n",
    "    Ratom = np.sqrt(density/(N_atoms*np.pi))   ## Radius of disk.\n",
    "    \n",
    "    spacing = 1/(np.sqrt(N_atoms)+1)\n",
    "    \n",
    "    x_spacings = np.linspace(0+spacing,1-spacing,int(np.sqrt(N_atoms)))\n",
    "    y_spacings = x_spacings\n",
    "    \n",
    "    position_list = np.array([(i,j) for i in x_spacings for j in y_spacings])\n",
    "    \n",
    "    return position_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#By using random numbers between 0 and 1 as the x and y velocity components,  \n",
    "#the function generates the initial velocity list for the particles.\n",
    "def velocities_f(N_atoms):\n",
    "    velocity_list = np.array([(random.random(),random.random()) for i in range(N_atoms)])\n",
    "    return velocity_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evolve(t,next_event,positions,velocities):\n",
    "    global x_components_list\n",
    "    \n",
    "    pos = np.copy(positions)\n",
    "    vel = np.copy(velocities)\n",
    "    if t+next_event >= int(t)+1:\n",
    "        t_previous = t\n",
    "        for snapshot_instant in range(int(t)+1,int(t)+2):\n",
    "            delta_t = snapshot_instant - t_previous\n",
    "            for i in range(Natoms):\n",
    "                pos[i] += vel[i]*delta_t\n",
    "            t_previous = snapshot_instant\n",
    "            x_components_list = np.append(x_components_list,x_component(pos))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "######################### Initialization #################################\n",
    "L = 1. ## Box edge length\n",
    "Natoms = 4 # Number of disks\n",
    "\n",
    "## List indexing all pairs of atoms.\n",
    "pairs = [[i,j] for i in range(Natoms) for j in range(i+1, Natoms)] \n",
    "\n",
    "t = 0\n",
    "n_events = 1000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "positions = positions_f(Natoms)\n",
    "velocities = velocities_f(Natoms)\n",
    "\n",
    "x_components_list = np.array([]) \n",
    "time_list=[]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████| 1000000/1000000 [01:58<00:00, 8453.38it/s]\n"
     ]
    }
   ],
   "source": [
    " ############### Event Loop ##############################################\n",
    "\n",
    "for event in tqdm(range(n_events)):\n",
    "    # Wall collision times for all disks and their velocity components.\n",
    "    wall_times = [wall_time(positions[i][j], velocities[i][j], Ratom) for i in range(Natoms) for j in range(2)] \n",
    "    # Pair collision times\n",
    "    pair_times = [pair_time(positions[i], velocities[i], positions[j], velocities[j], Ratom) for i,j in pairs] \n",
    "    # The next collision event is the minimum of wall and pair collision times.\n",
    "    next_event = min(wall_times + pair_times) \n",
    "    \n",
    "    if event> 500:\n",
    "        evolve(t,next_event,positions,velocities)\n",
    "        \n",
    "    t += next_event \n",
    "    for i in range(Natoms):\n",
    "        positions[i] += velocities[i]*next_event # Evolve positions to collision event\n",
    "        \n",
    "    if min(wall_times) < min(pair_times): # Check if next event is a collision with a wall\n",
    "        wall_index = wall_times.index(next_event)\n",
    "        particle, component = divmod(wall_index, 2)\n",
    "        velocities[particle][component] *= -1.0 ## Velocity component normal to wall changes sign\n",
    "    else:\n",
    "        pair_index = pair_times.index(next_event)\n",
    "        particle_1, particle_2 = pairs[pair_index] # Indices of particles participating in collision.\n",
    "        rel_pos = positions[particle_2] - positions[particle_1]\n",
    "        rel_vel = velocities[particle_2] - velocities[particle_1]\n",
    "        distance = np.sqrt(np.dot(rel_pos,rel_pos))\n",
    "        unit_perp = rel_pos/distance\n",
    "        scal_prod = np.dot(rel_vel,unit_perp)\n",
    "        velocities[particle_1] += scal_prod*unit_perp # Change in velocities of atoms colliding with each other\n",
    "        velocities[particle_2] -= scal_prod*unit_perp\n",
    "##################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD0CAYAAACPUQ0CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVWElEQVR4nO3df4xdZ33n8fdMTDCpbVoCtUNSNf0Rvh12qYHJCjvrgAkNVdNEuEFVhZugxskGQ5VCU5VN0JZGWyqVKl7IUlLIDysSapr+SC2VAmkQSbtxTKRlll0CvXyDU1CTCLvBBQdIJ4nt2T/umeTOnbn3nhnfuffOM++XZGnOc86555szmc997nOfc87YzMwMkqRyjA+7AElSfxnsklQYg12SCmOwS1JhDHZJKozBLkmFWTOsA09NTTnPUpKWYHJycqzb+qEFO8Dk5OSc5UajwcTExJCqqWfUaxz1+sAa+2HU6wNr7Jf2Gqempnru41CMJBXGYJekwhjsklQYg12SCmOwS1Jhus6KiYgXAXuBs4EXAx/KzL9tWX8J8EHgGLA3M2+NiHHgZmAz8AxwVWYeXJ7yJUntevXYLwOOZOb5wC8BfzK7ogr9jwBvBd4EXB0Rm4AdwNrM3ApcB+xZhrolSR30Cva/An6vZflYy88TwMHM/G5mPgvsB84HtgH3AGTmQ8C5/StXktRL16GYzPwBQESsB/4a+G8tqzcAR1uWvw+8dIH24xGxJjNb3xSA5sT7VtPT0/PaRs2o1zjq9YE19sOo1wfWCPDm89/M4SOH57VvPH0j9z9wf63XWEqNPa88jYifAPYBN2fmnS2rngLWtyyvB763QPv4QqEOzLviayVeBTZqRr0+sMZ+GPX6wBoBDh85zA3cMK/9hiM31D5u3688jYiNwL3Af83Mve3HA86JiJdFxKnAG4EvAg8CF1X7bwEerlX9Ip216SzGxsbm/Ttr01nLcThJWjF69dg/APwY8HsRMTvWfivwI5l5S0RcC/w9zTeIvZn5RETsAy6MiAPAGHDFchT+xOEnFn4nPDy/TZJWk15j7O8F3ttl/aeBT7e1nQB296U6SdKieYGSJBXGYJekwhjsklQYg12SCmOwS1JhDHZJKkxxwX4Kp3jhkqRVbagPs14OxznuhUuSVrXieuyStNoZ7JJ0kjrdu2pYihuKkaRB63jvqgXaBsEeuyQVxmCXpMIY7JJUGINdkgpjsEtSTaM2+6UTZ8VIUk2jNvulk1XTY/dWA5JWi1o99oh4A/DhzNze0rYJuKtls9cC12XmJyLiy8DRqv2bmbkszz1dDG81IGm16BnsEfF+4HLgh63tmXkI2F5tsxX4Q+DWiFhbrd/e51olSTXUGYp5FLi008qIGAM+Brw7M48Dm4HTIuLeiLgvIrb0p1RJUh09e+yZeXdEnN1lk0uAr2VmVstPAzcCtwHnAJ+LiMjMY+07NhqNOcvT09Pz2gZhMcccVo11jXp9YI39MOr1gTX2Uve4S6mxH7NiLgNuall+BDiYmTPAIxFxBDgDeKx9x4mJiTnLjUZjXtsgLOaYw6qxrlGvD6yxH0a9PrDGXuoet73Gqampnvv0Y1bMJHCgZXkXsAcgIl4JbAC+3YfjSJJqWHSwR8TOiLi6+vkVwPer3vms24EfjYj9wF8AuxYahpGkUbVSLkTqpNZQTGZ+C9hS/XxnS/uTNKc5tm77LLCzbxUus9n57e3O3Hgmjx96fAgVSRq2lXIhUier/spT57dLKs2qufJUklYLg12SCmOwS1JhDHZJKsyq//JU0up01qazeOLwE8MuY1kY7JJWpU5TGmHlTGvsxKEYSSqMwS5JhTHYJakwjrF30OlWAxtP38ih7xwaQkWSVI/B3kHHWw0cmd8mSaPEoRhJRVvpd2pcCnvskoq20u/UuBT22CWpMAb7Is1+qdr+76xNZw27NEkCHIpZNO/fLmnU2WOXpMLU6rFHxBuAD2fm9rb2a4ErgSerpncB3wBuBjYDzwBXZebBfhUsSequZ7BHxPuBy4EfLrD69cA7M3OqZftLgbWZuTUitgB7gLf1qV5JUg91hmIeBS7tsG4SuD4i9kfE9VXbNuAegMx8CDj3pKuUJNXWs8eemXdHxNkdVt8FfBx4CtgXERcDG4CjLdscj4g1mXmsfedGozFneXp6el7bSjIKta+Ec2iNJ2/U64OVUeMw1T03SzmPS54VExFjwEcz82i1/BngdTRDfn3LpuMLhTrAxMTEnOVGozGvbaU4hVN49atfPa/9zI1n8vihxwdWx0o4h9Z48ka9PlgZNQ5T3XPTfh6npqa6bN10MtMdNwBfjYgJmuPvFwB7gZcAlwB/WY2xP3wSx1gxnAYpDVfJT0RarEUHe0TsBNZl5i0R8QHgfpqzX76QmZ+NiHHgwog4AIwBV/S1YklawGq8dUAntYI9M78FbKl+vrOl/VPAp9q2PQHs7l+JkqTF8AIlSSqMwS5JhTHYl5k3DZP6azXeX32xvAnYMnO2jNRffknamz12SSqMwS5JhTHYJakwBrskFcZgl6TCGOxD4jRIScvF6Y5D4jRIScvFHrskFcZgHzEO0UhNXmG6dA7FjBiHaKQmrzBdOnvsklQYg12SCmOwS1Jhao2xR8QbgA9n5va29ncA7wOOA18B3pOZJyLiy8DRarNvZqaPx5OkAekZ7BHxfuBymg+sbm1/CfAh4DWZ+XRE/DlwcUTcC9D+JiBJC/Eh1P1Xp8f+KHApbc82pfkA6/My8+mW15oGNgOnVQG/BvhAZj7Up3olFcbZL/3Xc4w9M+8Gnlug/URmHgaIiGuAdcDngaeBG4FfpPlQ6z+LCKdVStKAnFTgRsQ48MfAq4C3Z+ZMRDwCHMzMGeCRiDgCnAE81r5/o9GYszw9PT2vTU2zFy4tZOPpG7n/gfuBlXEOrfHkjXp9sDJqHKa652Yp5/Fke9KfpDkksyMzT1Rtu4DXAO+JiFcCG4BvL7TzxMTEnOVGozGvTU2dLlwCuOHIDc+ft5VwDq3x5I16fbAyahymuuem/TxOTU313GfRwR4RO2kOu3wJuBJ4ALgvIgBuAm4H7oiI/cAMsCszjy32OJKkpakV7Jn5LWBL9fOdLas6jdHvPLmypMHoNCNj4+kbOfSdQ7W3P3PjmTx+6PFlqbEUzn4ZHL/U1KrWcUbGkflt3bb/g8N/sOB3IAb+C5z9MjgGu1aF5e4tevO2F9gzHz6DXSvOUoZDFttb7DYLSd3ZMx8+g10rzmKHQ5aiYw/ccNIKYLAXoFPvcrWN73adEjqkQO70uzl1/FSePfHsvPaV9DtzyGV0GewFWOnju50ColP4rSQdfzcnbljxX8I65DK6DHYNXceA6BB+JQfHSn+T1mjwfuwF8/mp0upkj71go9b7c0xWGgyDfRVa7i9buwX4ahtaKYFvyCuPwb4KLXdP3i/V+m8Qs2t8Qy6Hwa6evD/K8C12ds1S3qR9Qy6Hwa7ndbvacpTG6iV1Z7DreV5tuTo4Zl4+g11L5v1UViaHXMpnsGvJ7OGPLt90VzeDXSrQKN43R4PjlaeSVJhaPfaIeAPw4czc3tZ+CfBB4BiwNzNvjYhx4GZgM80HXV+VmQf7WrUkqaOePfaIeD9wG7C2rf1FwEeAtwJvAq6OiE3ADmBtZm4FrgP29LlmSVIXdYZiHgUuXaB9AjiYmd/NzGeB/cD5wDbgHoDMfAg4t0+1SpJq6BnsmXk38NwCqzYAR1uWvw+8dIH24xHhl7SSNCAnE7hPAetbltcD31ugfTwzjy30Ao1GY87y9PT0vDZJKlHdrFtKLp5MsDeAcyLiZcAPgDcCNwIzwCXAX0bEFuDhTi8wMTEx9wUbjXltklSiulnXnotTU1M991n0dMeI2BkRV2fmc8C1wN8DX6Q5K+YJYB8wHREHaH65+tuLPYYkaelq9dgz81vAlurnO1vaPw18um3bE8Du/pUoSVoML1CSpMIY7JJUGINdkgpjsEtSYQx2SSqMwS5JhTHYJakwBrskFcZgl6TCGOySVBiDXZIKY7BLUmEMdkkqjMEuSYUx2CWpMAa7JBXGYJekwhjsklSYno/Gi4hx4GZgM/AMcFVmHqzWbQLuatn8tcB1mfmJiPgycLRq/2ZmXtHPwiVJC6vzzNMdwNrM3BoRW4A9wNsAMvMQsB0gIrYCfwjcGhFrq/Xb+1+yJKmbOkMx24B7ADLzIeDc9g0iYgz4GPDuzDxOs3d/WkTcGxH3VW8IkqQBqNNj38ALQyoAxyNiTWYea2m7BPhaZma1/DRwI3AbcA7wuYiItn1oNBpzDjQ9PT2vTZJKVDfrlpKLdYL9KWB9y/J4e0ADlwE3tSw/AhzMzBngkYg4ApwBPNa608TExJwXaTQa89okqUR1s649F6empnruU2co5kHgIoBqSOXhBbaZBA60LO+iORZPRLySZq//2zWOJUk6SXV67PuACyPiADAGXBERO4F1mXlLRLwC+H7VO591O3BHROwHZoBdC/TyJUnLoGewZ+YJYHdb89db1j9Jc5pj6z7PAjv7UJ8kaZG8QEmSCmOwS1JhDHZJKozBLkmFMdglqTAGuyQVxmCXpMIY7JJUGINdkgpjsEtSYQx2SSqMwS5JhTHYJakwBrskFcZgl6TCGOySVBiDXZIKY7BLUmF6PhovIsaBm4HNwDPAVZl5sGX9tcCVwJNV07uAb3TbR5K0fOr02HcAazNzK3AdsKdt/euBd2bm9upf1thHkrRM6gT7NuAegMx8CDi3bf0kcH1E7I+I62vuI0laJj2HYoANwNGW5eMRsSYzj1XLdwEfB54C9kXExTX2AaDRaMw50PT09Lw2SSpR3axbSi7WCfangPUty+OzAR0RY8BHM/NotfwZ4HXd9mk1MTExZ7nRaMxrk6QS1c269lycmprquU+doZgHgYsAImIL8HDLug3AVyNiXRXyFwBTPfaRJC2jOj32fcCFEXEAGAOuiIidwLrMvCUiPgDcT3P2yxcy87PVTJo5+yxT/ZKkNj2DPTNPALvbmr/esv5TwKdq7CNJGgAvUJKkwhjsklQYg12SCmOwS1JhDHZJKozBLkmFMdglqTAGuyQVxmCXpMIY7JJUGINdkgpjsEtSYQx2SSqMwS5JhTHYJakwBrskFcZgl6TCGOySVJiej8arnl96M7CZ5nNNr8rMgy3r3wG8DzgOfAV4T2aeiIgvA0erzb6ZmT73VJIGoM7DrHcAazNza0RsAfYAbwOIiJcAHwJek5lPR8SfAxdHxL0Ambl9WaqWJHVUZyhmG3APQGY+BJzbsu4Z4LzMfLpaXgNM0+zdnxYR90bEfdUbgiRpAOr02DfwwpAKwPGIWJOZxzLzBHAYICKuAdYBnwf+I3AjcBtwDvC5iIjMPNb6wo1GY86Bpqen57VJUonqZt1ScrFOsD8FrG9ZHm8N6GoM/o+BVwFvz8yZiHgEOJiZM8AjEXEEOAN4rPWFJyYm5hyo0WjMa5OkEtXNuvZcnJqa6rlPnaGYB4GLAKohlYfb1n8SWAvsaBmS2UVzLJ6IeCXNXv+3axxLknSS6vTY9wEXRsQBYAy4IiJ20hx2+RJwJfAAcF9EANwE3A7cERH7gRlgV/swjCRpefQM9mocfXdb89dbfu7U69+51KIkSUvnBUqSVBiDXZIKY7BLUmEMdkkqjMEuSYUx2CWpMAa7JBXGYJekwhjsklQYg12SCmOwS1JhDHZJKozBLkmFMdglqTAGuyQVxmCXpMIY7JJUGINdkgrT89F4ETEO3AxsBp4BrsrMgy3rLwE+CBwD9mbmrb32kSQtnzo99h3A2szcClwH7JldEREvAj4CvBV4E3B1RGzqto8kaXnVCfZtwD0AmfkQcG7LugngYGZ+NzOfBfYD5/fYR5K0jMZmZma6bhARtwF3Z+bnquV/AX46M49FxDbgmsz8tWrdfwf+BdjSaZ/Z152amup+YEnSgiYnJ8e6re85xg48BaxvWR5vCej2deuB7/XYp1ZhkqSlqTMU8yBwEUBEbAEeblnXAM6JiJdFxKnAG4Ev9thHkrSM6gzFzM5w+XlgDLgCeD2wLjNvaZkVM05zVszHF9onM7++fP8ZkqRZPYO93+pMhYyI04DPA1cO4w2hxhTPdwDvA44DXwHek5knRqzGt9OckTQD3JKZt41SfS3b3QL8W2ZeN8j6qmP3OofXAlcCT1ZN78rMHLEa/xPwP2h2oA4Bl2Xm9KjUWM2Su6tl89cC12XmJ0ahvmr9rwO/Q/PveW9m/umgaltEjZcDvwscBe7IzNu7vd4wLlDaQZepkBFxLvC/gJ8ZfGnP20HnKZ4vAT4EvDkzzwNeClw8YjWeAvwR8AvAVuB3I+Llo1LfrIh4F/CaAdfVagfda3w98M7M3F79G2ioV3bQ+fc8BtxK8xPx7Ey0nxylGjPz0Oz5A64H/g/NmkeivsqNNP9W/jPwOxHxY4MtD+j+e345zczZTnNa+a9HxNndXmwYwd5rKuSLgV8Bhjl0063GZ4DzMvPpankNMNAeUqVjjZl5HJjIzKPA6TR7cz8YlfoAImIrzdlTnxxwXa16/b84CVwfEfsj4vpBF1fpVuOrgCPA+yLiH4GXDenNp+f05upN6GPAu6v/PwepV31fodlBW0vzb2UYM/a61fjTwP/NzH+rRgb+N82/nY6GEewbaH6cmHU8Ip6fnZOZD2bmY4Mva46ONWbmicw8DBAR1wDraA4bDVqv83gsIi4F/h/NT0DPjUp9EXEGcAPwmwOuqV3Xc0hzCGE3cAGwLSKG8cmsW40vB86j+RH+F4C3RMRbBlwf9D6PAJcAXxvSG0+v+r4KTAFfA/4uM783wNpmdavxG8B/iIiN1TD1W4Af6fZiwwj2nlMhR0DXGiNiPCJuBC4E3p6Zw3iH73keM/NvgDOBU4F3DrA26F7fr9IMpc/S/Ni5MyJ+Y7DlAV1qrHqYH83M71QX330GeN0o1Uizt34wM/8pM5+j2eObHHSB1Pubvgy4ZXAlzdHt9/zzwC8DPwWcDfx4RPzqwCvsUmNmfhf4beBuYC/N4azvdHuxYQT7SpgK2avGT9L82LajZUhm0DrWGBEbIuIfI+LF1Ue3HwID/XK3W32Z+T8zc7Iad/0j4M7MvGPA9XWtkWYP6qsRsa4K+Qto9uoGrVuN/wysi4ifrZbPp9nrHLQ6f9OTwIFBFtWiW31HgX8H/r0aIvpXYBhj7N3+ntfQHHp5I80O2s9V23c0zFkxC06fbNnuH4DdQ54VM69G4EvVvwd4YSzupszcNyo1VtNQr6Y5o+M5mmOI1wxybHMRv+ffAH5uyLNiOp3Dy4Hfovm9yhcy8/dHsMYLaL45jgEHMvO9I1jjK4DPZ+ZrB11bzfp2A7uAZ4FHgf9SfUobpRp/n+YXrNPAnsz8626vN/BglyQtL+/HLkmFMdglqTAGuyQVxmCXpMIY7JJUGINdkgpjsEtSYQx2SSrM/wfXkfVsDnFdcwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_bins = 50\n",
    "plt.hist(x_components_list, bins = n_bins, density=True, color='purple', edgecolor='black');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPMAAAD3CAYAAADIZ2IpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUYklEQVR4nO3de3CU1f0G8GdJIAm3AKEUpCBycZFCAmwC/Cq2DBUcIMhlisOIiCByt9KxYMDriG1FpEorhLQUUxUoKv6gOiNSbJC7NgsoiK7cKQ5SkMotmwBh+8e3KRAI2WTf95zznn0+MzsR2bzn+4Y8e97Lec/xRSIREJH31dBdABE5g2EmsgTDTGQJhpnIEgwzkSUYZiJLRBVmv9/f3e/3r3O5FiKKQWJlb/D7/dMBjARwzv1yiKi6oumZ9wEY6nYhRBSbSnvmUCi0wu/3t6ro74PBIIeQESkQCAR8N/r7SsMcZSNObIaIKhAMBit9D69mE1mCYSayRFSH2aFQ6CCAHu6WQkSxcOScmQwWDgM7dgCFhUAwCOzdCxQVASUlgM8HJCcD9eoBP/whEAjIq0MHIJG/Gl7DfzEb7d8PLFwIvP8+8MUXQGlp5d+zbt3l/05JATIygGHDgNGjgYYNXSuVnMNzZltcugS89x7Qrx/Qti0wZw6wa1d0QS4vHAa2bgUefRRo3hwYM0Z6dTIaw2yDpUuBNm2AgQOB1asBJ2ePCYeBV18FMjOBHj2ATz5xbtvkKIbZy775Bhg8GBgxAjh40P32Pv4Y+NGPgMcek3NuMgrD7FVLlshFq1Wr1LZbWgq88ALQpQt7acMwzF5TUgIMHw7cdx9w8qS+Or74Qnrp55/XVwNdhWH2knPngOxsYPly3ZWI0lJgxgxg2jTdlRAYZu8Ih+VK9dq1uiu51osvAj//ue4q4p5d95mPHZNbKIcOyS9/SQmQkCD3TVNTgfR0Oc+sWVN3pVVz4QIwdCiwYYPuSir2+9/Lz3jWLN2VxC1vh3nzZuDDDyXAwSBw5Ejl35OUJKEOBIBu3YAhQ4AGDVwvNSZTp8otJ9M99xzQvr1cXSf1IpFITK/CwsKIUmfORCK5uZFIenokIndUY3vVrh2JjBkTiQSDavcjWh9+GIn4fM7sq4pXw4aRyNGjun9q1vlvzm6YRe+cMx88CEyZIiOSJk4EPvvMme0WFQGLF0tP3b273PIxZcmes2eBBx80p55o/PvfwPjxuquIS+aHORIB5s8HOnaUr6dPu9fWJ5/ILZ+f/hQ4cMC9dqI1fbqawSBO++tfgTfe0F1F3DE7zAcOAL17S498TuF8ggUFQKdOwIIF+nrF9evlYQmveuQR4MQJ3VXEFXPD/Ic/SKCufJpHpXPngMmTpZf++mv17c+e7a3D6/JOnpR/Q1LGzDDPmCHnXSp744oUFAC33w7s2aOuzQMHvHH1ujJ5efI0FylhXpgffti8IYKHDgF33AF8/rma9hYutCMEhw/LY5mkhFlhfuwx4JVXdFdxfceOAX36yIP/biopkavrtliwQHcFccOcML/0kjyNY7KjR4E77wS+/da9NlautOvC0Zo1cmRDrjMjzLt2ATk5uquIzoEDcirgFpOHbFZHJAJs3Ki7irigP8wXLwIPPACcP6+7kugtWyY9qBsKC93Zrk6cckgJ/WGePdub/9gTJjj/PPHFi8Cnnzq7TRN48d/Xg/SGedcu4NlntZZQbceOOX+4vXs3UFzs7DZNsH27t++Ze4TeMD/xhLcOr8tbulTmpHaKrT3YmTPAV1/prsJ6+sJsyz1IJ2+9HD3q3LZMY/O+GUJfmPPyqjens2mWLAFOnXJmWzYeYpexed8MoSfM588Df/qTlqYdV1QE5Oc7sy2bp6+1ed8MoSfM77wjF5BskZvrzHa8Np1RVdi8b4bQE+Z339XSrGtCIWcexEhJiX0bprJ53wyhJ8w2XrV1Yp/S0mLfhqls3jdDqA+zrbcpnAhzly6xb8NESUnAbbfprsJ66sO8bZudAwicGIaZkWHnusjp6TxnVkB9mG08xAac+ZBKTpZ5vW2Tmam7grigPsw2HmIDMtGgE1foA4HYt2EaG/fJQOrDXFSkvEllnNi3Hj1i34ZpbNwnA6kPs80jgcLh2Ldxzz1AnTqxb8cUWVl2njoYSP8jkDbx+WLfRmoqcO+9sW/HFJMm6a4gbqgPs82DB5zaN1sC0KiRrCVNSqgPs02HkOU5tW+dO9txnjl6tFyhJyXUh7lDB+VNKpGWBjRp4tz2Hn3UuW3pkJQkiwiQMurDbOttCqf362c/AwYOdHabKj3zDHDLLbqriCvqw5yRIQug28aND6m8PKBhQ+e367asLGDaNN1VxB31Ya5d285xum6EuVkzYN4857frpqQkeb7bxg9sw+m5NWXj8D63Th9GjvTW4fbTT9t7XcRwesI8ZIiWZl3TtSvQqpV728/Pl/WpTXfPPbLEEGmhJ8zZ2cDNN2tp2hUTJ7q7/UaNgL/9DWjTxt12YtGvnyywXoPjkHTR85OvUQMYN05L045r0EDNiK2mTWWt6vbt3W+rqu6+W6aC4mOOWun7GB07FqhVS1vzjnngAbmop8IPfgCsX2/W7b2RI4EVKzg4xAD6wtykidxL9TKfz/1D7PK+9z1g82Zg5ky9V4wbNAD+/GfgtdfsnFDBg/Se4PzqV0DdulpLiMnkycCtt6pvt1Yt+dlt3arniaTsbFl4/v771bdNFdIb5latzF+TuSKtWwPPP6+3hsxMmbllxgy5v+u2xo3lyvq77wI33eR+e1Ql+i89TpgA9O6tu4qq8fmAxYvNeGgkKQn49a+BI0fkw8WNIZTdu0uI//lPYNQo57dPjtAfZp9PVrfw0uH2pEnAT36iu4qrNW4s93j37pU1vAYMiK23Tk0FxoyRnn/rVgkxL3IZzYwrF61ayaoQ999v/sydnTvLmtKmqlFDgjxgAHDhgpzbFhZKKINBmaw/HJblYnw+CWjdujIoJRC4/Grb1pnJFkgZM8IMAPfdB3z7LTB1qu5KKnbrrcAHH5hxeB2NmjXlw6dzZ7kVeKVIhGG1jP7D7Cs98gjwm9/oruL62rYF1q519pllnRhk65gVZgDIyQHmzzfrly09Hdi4EWjRQnclRBUyL8yAXGBatUqGMOo2YgTw0UfA97+vuxKiGzIzzIA89vf55/pmqmzaFFi5Uh4eaNBATw1EVWBumAF5WmjJEgmVyl56xAj5IBk0SF2bRDEyO8xlBg0Cdu+Wi2NuPTecmAgMHQoUFEhv3KiRO+0QucQbYQZkLqycHGDfPhlO2K+fM8/ONmsGPPUUcPCgPP3Tq1fs2yTSwJz7zNGqUUMG+mdnSwALCmQwRGEh8OmnlS9/07Ll5YERWVkylJRP/ZAFvP1b3KqVjBrLzJTXzp3A4cPAyZPAxYtyeyslRZ4DTk+Xw/WWLXVXTeQKb4Z5/Xrg7belN96xI/oF237xC5lsLhCQw/QhQzg7BlnDO2E+c0YehM/NlSvN1VFaKr33zp3yFFDTpjLMcfx46b2JPMz8C2DHjwNTpgDNm8vX6gb5er75BnjuOTlcHzJEQk7kUWaH+a23ZCaN+fOlZ3ZLaancy87MBGbNkvNtIo8xM8zHjwPDhsk8zMePq2v3/Hm5TdW9O3tp8hzzwrx6tfTGb7+tr4Zt26SXnjNHXw1EVWRWmJculTHZKnvjipw/D0yfLlfATZ8wgQgmhXnxYpmD2bTz1Zdflgn7GWgynBlhfust4KGHgEuXdFdyfYsWeX/xc7Ke/jBv2SJTBpka5DIvvSS9NJGh9IY5HJZZH8+f11pG1GbMAL78UncVRNelN8yPPy6zRXpFcTEwerT5RxEUl/SFeeNGYN48bc1X29atwG9/q7sKomvoCXNxsUyw7tUe7skngVBIdxVEV9ET5tdf99bhdXnFxeZOCUxxS0+Yc3O1NOuo5cvluWkiQ6gP85YtwPbtypt1XHGxDHQhMoT6MNvQK5dZuJAjw8gYasN88iTw5ptKm3TVvn3AmjW6qyACoDrMGzbI6oM2WbtWdwVEAFSHORhU2pwSNu4TeZLaMBcWKm1OiW3beN5MRmDPHKtTp+TcmUgzdWE+cgT417+UNaeUjR9S5DnqwnzwoLKmlDtwQHcFRArDHO1E9V5k876RZ6gLc2VrQHkZw0wGUBfmhARlTSnHhefIAOrCnJysrCnlUlJ0V0CkMMz16ytrSjmb9408Q12YO3Sw91A7PV13BUQKw1y7NnDbbcqaU8bnA7p21V0FkeIRYIGA0uaUaNMGSE3VXQWR4jBnZiptTgkbP6DIk9SGOStLaXNK2PgBRZ6kNszdugGtWytt0lUJCbLsLJEB1IbZ5wMmTFDapKsGDABattRdBREAHXOAjRljzwCSyZN1V0D0P+rDnJZmx6Fpu3ZAnz66qyD6Hz3zZk+ZoqVZR02eLKcNRIbQE+asLFmP2asyMoBJk3RXQXQVfQvHzZ3rzYtHNWsC+fnylcgg+sJcrx6waJG25qtt5kygc2fdVRBdQ+/6zH36eOtwOyND1pQmMpDeMAOy1rEXhkSmpQHLlvHwmoylP8x16wKrV5v9RFW9eubXSHFPf5gBoHFjWebFxLCkpgLvv88x2GQ8M8IMADfdBKxfb1ZomjQBCgqA22/XXQlRpcwJMyA99EcfAQ8/rH9ARt++spxOly566yCKkllhBmRGkt/9Dli3Th78V61+feCPfwQ++ABo0UJ9+0TVZF6Yy/z4x8Bnn6ntpfv2BXbtAsaOVdMekYPMDTNwuZfevVtC7cb0PImJwLBhcm7M3pg8zOwwl2nfXkL99ddAXp4M3ohV8+bAM88Ahw8Db74J9OoV+zaJNPLWUgx16gDjxslrzx5ZfbHstW2bLK96PUlJMh1uIHD51akTV6Igq3j3t7ldO3kNHy5/jkRk2dizZ2Xtp0uXZKWJlBQ5dObILbKcd8Ncns/H8123RSLAoUPAmTOXPzCTk+UDs2VLLtOjmT1hJuft2SP32oNB+bp9O3D69PXfm5Agq5aUncZkZsriALVqqa05jjHMdLWiImDpUiA3V65DRKu0FNi5U175+fL/GjcGRo8GJk4EbrnFlXLpMm9czSb3ffUVMHWqXOV/6KGqBbkiJ04Ac+YAbdsC/fsD770nh+bkCoY53hUXA9OmyUMu8+YB333nfBuXLsnDKgMHAnfcIR8c5DiGOZ5t2SKzprz4oroec/NmaXPuXPbSDmOY41FZb9yzJxAKqW8/HAZ++Uv20g5jmOPNqVPAnXeq7Y0rsnmzLFm0YYPeOizBMMeT48dl2OqmTboruezUKeCuu+ScmmLCMMeL776Tp8J27NBdybXCYWDIEODvf9ddiacxzPGguBjIzjYzyGVKSoBBg4B//EN3JZ7FMMeDxx8369C6ImfPylj7oiLdlXgSw2y7TZuAl1/WXUX09u8HcnJ0V+FJDLPNwmFZQlf3VeuqeuUVmQuOqoRhttkTT3jzPm4kIh9C587prsRTGGZbffmltw6vy9u/H3jhBd1VeArDbKv58713eF1eXh5w4YLuKjyDYbbRuXPAa6/priJ2x44BK1borsIzGGYbvfFGxZMIeM2CBbor8AyG2UY2BWDDBpnLnCrFMNtm925ZPMAmf/mL7go8gWG2TWGh7gqcZ+M+uYBhtk0wqLsC59m4Ty5gmG1jYy924oRM8Us3xDDbpLTU7CejYsHeuVIMs03277f3iSPbLuq5gGG2iS33lq/nzBndFRiPYbZJcbHuCtxj8745hGG2SWmp7grcY/O+OYRhtklysu4K3GPzvjmEYbaJzaswMsyVYphtcvPNsrStjVq31l2B8Rhmm9SvLwvQ2ygQ0F2B8Rhm29j4S1+rFtCpk+4qjMcw28bGMHfsyEXbo8Aw2yYzU3cFzrPxA8oFDLNtsrKA1FTdVTirTx/dFXgCw2yb2rWBUaN0V+GcZs2AwYN1V+EJDLONJk3SXYFzxo4FatbUXYUnMMw28vuB3r11VxG7xERg/HjdVXgGw2wrG3rngQOB5s11V+EZDLOtBg/29lXgxETgqad0V+EpDLOtEhKAV1/17v3ZmTOBzp11V+EpDLPNOnUCnnxSdxVVl5Ehi95RlTDMtsvJ8dbhdmIikJ/PK9jVwDDbriwcderoriQ6Tz/Nw+tqYpjjQceOwDvvmH/+PH48D69jwDDHi759gWXLzD18vfdeu9bI0oBhjidDhwIrV5o3I8m4ccDrrwM1+OsYC/704k3//sCaNUCLFrorkfP5Z5+VRdUZ5JjxJxiPevaUZVIffFBfDZ06AR9/7M1bZ4ZimONV/frAokXA6tVqe+nERLnIVVgIdO2qrt04wDDHu7vukl56+nQgLc29dhISgLvvlt541izzr6x7EMNM0kvPng0cOSL3pLt1c27bTZrI0Mz9+4FVq9gbuyhRdwFkkORkmdhg1Cg5DF6+XL5u2xb9OlaJiUCHDjLqrG9fuYLOXlgJhpmuLzPz8nxikQiwd68sq7pzpyziFg7LkjEpKfJq3VoCnJ5u3q2vOMEwU+V8PpmPu107YPhw3dVQBXjOTGQJhpnIEgwzkSUYZiJLMMxElmCYiSzBMBNZgmEmsgTDTGQJhpnIEgwzkSUYZiJLMMxElmCYiSzBMBNZgmEmsgTDTGQJhpnIEgwzkSUYZiJLMMxElmCYiSzBMBNZgmEmsgTDTGQJhpnIEgwzkSUYZiJLMMxElmCYiSzBMBNZgmEmsgTDTGSJxMre4Pf7awBYACADQAmAsaFQaK/bhRFR1UTTMw8GkBwKhf4PQA6Aua5WRETVUmnPDKAngNUAEAqFtvr9/szybwgGg07XRURVFE2Y6wM4dcWfS/1+f2IoFLoIAIFAwOdKZURUJdEcZp8GUO/K7ykLMhGZI5owbwLQHwD8fn8PADtdrYiIqiWaw+z/B9DH7/dvBuADMNrdkoioOnyRSKRa38hbVkRq+P3+7gBmh0KhXjd6XyyDRgaDt6yIXOX3+6cDWAQgubL3xhLmq25ZAbjmlhURxWwfgKHRvDGWMF/3llUM2yOickKh0AoAF6J5byxh5i0rIoPEEmbesiIySCyHxbxlRWSQat+aIiKz8HlmIkswzESWYJiJLMEwE1mCYSayBMNMZAmGmcgS/wFw/wiog7rbsAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "create_picture(positions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "92775.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(x_components_list)/4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}