diffsimviz.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global/Universal Initializations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "import numpy\n",
    "import numpy as np\n",
    "import os\n",
    "import imageio.v2 as imageio\n",
    "from math import sqrt\n",
    "from scipy.stats import norm\n",
    "from IPython.display import display, clear_output\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "numofmol = 8\n",
    "N = 30\n",
    "radius = 0.2\n",
    "vradius=0.2\n",
    "theta = 45\n",
    "\n",
    "delta = 0.4\n",
    "# Total time.\n",
    "T = 7\n",
    "# Time step size\n",
    "dt = T/N\n",
    "\n",
    "traj = np.zeros((2,N,numofmol))\n",
    "initialrot = np.zeros((numofmol,2))\n",
    "where = np.zeros((numofmol))\n",
    "\n",
    "ratiointra = 0.5\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function is also used in the stick cross generation\n",
    "def make_stick_line(theta,radius):\n",
    "\n",
    "    y = np.linspace(-2,2,100)\n",
    "\n",
    "    xlow = np.zeros((100,))\n",
    "    xlow[:] = -radius\n",
    "    xhigh = np.zeros((100,))\n",
    "    xhigh[:] = radius\n",
    "\n",
    "    line1 = np.c_[xlow,y]\n",
    "    line2 = np.c_[xhigh,y]\n",
    "\n",
    "    c, s = np.cos(theta), np.sin(theta)\n",
    "    rot = np.array(((c, -s), (s, c)))\n",
    "\n",
    "    line1rot = line1 @ rot\n",
    "    line2rot = line2 @ rot\n",
    "    \n",
    "\n",
    "    return line1rot, line2rot\n",
    "line1, line2 = make_stick_line(theta,radius)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vertical Stick"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "brownian() implements one dimensional Brownian motion (i.e. the Wiener process).\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "def brownian_singlestep(x0, n, dt, delta, radius, theta, out=None):\n",
    "    \"\"\"\n",
    "    Generate an instance of Brownian motion (i.e. the Wiener process):\n",
    "\n",
    "        X(t) = X(0) + N(0, delta**2 * t; 0, t)\n",
    "\n",
    "    where N(a,b; t0, t1) is a normally distributed random variable with mean a and\n",
    "    variance b.  The parameters t0 and t1 make explicit the statistical\n",
    "    independence of N on different time intervals; that is, if [t0, t1) and\n",
    "    [t2, t3) are disjoint intervals, then N(a, b; t0, t1) and N(a, b; t2, t3)\n",
    "    are independent.\n",
    "    \n",
    "    Written as an iteration scheme,\n",
    "\n",
    "        X(t + dt) = X(t) + N(0, delta**2 * dt; t, t+dt)\n",
    "\n",
    "\n",
    "    If `x0` is an array (or array-like), each value in `x0` is treated as\n",
    "    an initial condition, and the value returned is a numpy array with one\n",
    "    more dimension than `x0`.\n",
    "\n",
    "    Arguments\n",
    "    ---------\n",
    "    x0 : float or numpy array (or something that can be converted to a numpy array\n",
    "         using numpy.asarray(x0)).\n",
    "        The initial condition(s) (i.e. position(s)) of the Brownian motion.\n",
    "    n : int\n",
    "        The number of steps to take.\n",
    "    dt : float\n",
    "        The time step.\n",
    "    delta : float\n",
    "        delta determines the \"speed\" of the Brownian motion.  The random variable\n",
    "        of the position at time t, X(t), has a normal distribution whose mean is\n",
    "        the position at time t=0 and whose variance is delta**2*t.\n",
    "    out : numpy array or None\n",
    "        If `out` is not None, it specifies the array in which to put the\n",
    "        result.  If `out` is None, a new numpy array is created and returned.\n",
    "    sticks : int\n",
    "        The number of sticks (pairs of wall) in the simulation.\n",
    "        Used for conditionals as more collisions will need to be checked.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    A numpy array of floats with shape `x0.shape + (n,)`.\n",
    "    \n",
    "    Note that the initial value `x0` is not included in the returned array.\n",
    "    \"\"\"\n",
    "\n",
    "    x0 = np.asarray(x0)\n",
    "    # For each element of x0, generate a sample of n numbers from a\n",
    "    # normal distribution.\n",
    "    coordhold = np.copy(x0)\n",
    "    coordfinal = np.zeros((2,n))\n",
    "    #coordfinal = []\n",
    "\n",
    "    c, s = np.cos(theta), np.sin(theta)\n",
    "    rot = np.array(((c, -s), (s, c)))\n",
    "\n",
    "    if x0[0] < -radius:\n",
    "        where=0\n",
    "        for i in range(n):\n",
    "            if i == 0:\n",
    "                coordfinal[:,0] = x0\n",
    "            #print(i)\n",
    "            else:\n",
    "                r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                if (r[0] + coordhold[0]) < -radius:\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                else:\n",
    "                    while (r[0] + coordhold[0]) > -radius:\n",
    "                        r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "    elif x0[0] > radius:\n",
    "        where=0\n",
    "        for i in range(n):\n",
    "            if i == 0:\n",
    "                coordfinal[:,0] = x0\n",
    "            #print(i)\n",
    "            else:\n",
    "                r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                if (r[0] + coordhold[0]) > radius:\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                else:\n",
    "                    while (r[0] + coordhold[0]) < radius:\n",
    "                        r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "    else:\n",
    "        where = 1\n",
    "        for i in range(n):\n",
    "            if i == 0:\n",
    "                coordfinal[:,0] = x0\n",
    "            #print(i)\n",
    "            else:\n",
    "                r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                if (r[0] + coordhold[0]) > -radius and (r[0] + coordhold[0]) < radius:\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                else:\n",
    "                    while (r[0] + coordhold[0]) < -radius or (r[0] + coordhold[0]) > radius:\n",
    "                        r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "\n",
    "    coordfinal = np.squeeze(np.array(coordfinal)) \n",
    "    # Add the initial condition.\n",
    "    #coordfinal += np.expand_dims(x0, axis=-1)\n",
    "    coordfinal = np.transpose(rot) @ coordfinal\n",
    "    initialrot = np.transpose(rot) @ x0\n",
    "    \n",
    "    return coordfinal,where,initialrot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "initialinx = numpy.random.uniform(low=-radius, high=radius, size=(round(numofmol*ratiointra),1))\n",
    "initialiny = numpy.random.uniform(low=-1.5, high=1.5, size=(round(numofmol*ratiointra),1))\n",
    "initialincomb = np.concatenate((initialinx,initialiny),axis=1)\n",
    "\n",
    "initialleftx = numpy.random.uniform(low=-1.5, high=-radius, size=(round(numofmol*(1-ratiointra)*0.5),1))\n",
    "initiallefty = numpy.random.uniform(low=-1.5, high=1.5, size=(round(numofmol*(1-ratiointra)*0.5),1))\n",
    "initialleftcomb = np.concatenate((initialleftx,initiallefty),axis=1)\n",
    "\n",
    "initialrightx = numpy.random.uniform(low=radius, high=1.5, size=(round(numofmol*(1-ratiointra)*0.5),1))\n",
    "initialrighty = numpy.random.uniform(low=-1.5, high=1.5, size=(round(numofmol*(1-ratiointra)*0.5),1))\n",
    "initialrightcomb = np.concatenate((initialrightx,initialrighty),axis=1)\n",
    "\n",
    "initial = np.concatenate((initialincomb,initialleftcomb,initialrightcomb),axis=0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAASI0lEQVR4nO3dX6hlZ3nH8e/POHMTZEab/3+UCENtvGhJDkNSoURqJRmEUVFIvDCIMETMZS8Cgt62vRRT41CCyUUN3kSHOhqjUGIv0uaMJJoYU6fBNuOEZIxkJChOgk8vzho5mexzzj5nv2fvtdb+fuBw9p81e73rWcn6vc9ee6+TqkKSpLctegCSpH4wECRJgIEgSeoYCJIkwECQJHUMBEkS0CgQktyf5OUkT2/w/C1JziZ5svv5Yov1SpLaeXuj1/k68BXgwU2W+VFVfaTR+iRJjTXpEKrqMeA3LV5LkrQYrTqEadyc5CngNPD3VfXMpIWSHAGOAFx88cU3vu9975vjECVp2E6cOPHrqrp0J/92XoHwY+A9VfVakkPAt4ADkxasqqPAUYCVlZVaXV2d0xAlafiS/O9O/+1cPmVUVb+tqte628eBPUkumce6JUnTmUsgJLkiSbrbB7v1vjKPdUuSptPkLaMk3wBuAS5Jcgr4ErAHoKruAz4BfC7JG8DvgdvLy6xKUq80CYSqumOL57/C2sdSJUk95TeVJUmAgSBJ6hgIkiTAQJAkdQwESRJgIEiSOgaCJAkwECRJHQNBkgQYCJKkjoEgSQIMBElSx0CQJAEGgiSpYyBIkgADQZLUMRAkSYCBIEnqNPkTmpKkxUoy82sYCJI0UC1CYD0DQZIGpnUQnGcgSNIAbBUCVTXVcpsxECSpxzY7wJ8PgVYMBEnqmWm7gdYMBEnqiXl2A5MYCJK0QIvqBiYxECRpARbdDUzS5JvKSe5P8nKSpzd4Pkm+nORkkp8kuaHFeiVpSJL86WeSqlpYGEC7S1d8Hbh1k+dvAw50P0eArzZaryT13jQhsMggOK9JIFTVY8BvNlnkMPBgrXkc2J/kyhbrlqQ+6ns3MMm8Lm53NfDCuvunusckaVSG0g1MMq9AmFSdiRVJciTJapLVM2fO7PKwJGl2Q+wGJplXIJwCrl13/xrg9KQFq+poVa1U1cqll146l8FJ0k4MuRuYZF6BcAz4dPdpo5uAs1X14pzWLUnNjKUbmKTJ9xCSfAO4BbgkySngS8AegKq6DzgOHAJOAr8DPtNivZI0L3383kBrTQKhqu7Y4vkCPt9iXZI0L336FvE8+E1lSbrAMnQDkxgIksTydQOTGAiSltqydgOTGAiSlo7dwGQGgqSlYTewOQNB0qjZDUzPQJA0SnYD22cgSBoNu4HZGAiSBs9uoA0DQdIg2Q20ZyBIGhS7gd1jIEjqPbuB+TAQJPWW3cB8GQiSesVuYHEMBEm9YDewePP6i2mS9Bb79+8f7V8fGyI7BElzZzfQ3lZvtU3DDkHSXNgNtLfV33feLjsESbvKbqC9VgFwITsESc3ZDbQ3j5raIUhqxm6gvXnW1A5B0kzsBtpbVE3tECTtiN1Ae4uuqR2CpKnZDbTXp5raIUja0qJnrmPUx5raIUiaqE8z17Hoe03tECS9SR9nrkM3lJraIUjq/cx1iIZY0yaBkOTWJM8lOZnkngnP35LkbJInu58vtlivpNmcP2CdPXv2Lc+dP2D17aDVd0Ou6cxvGSW5CLgX+DvgFPBEkmNV9bMLFv1RVX1k1vVJms3+/fsnHqzO6+vBqs/GUtMW5xAOAier6nmAJA8Bh4ELA0HSAg3lfewhGVtNW7xldDXwwrr7p7rHLnRzkqeSfDfJ+zd6sSRHkqwmWT1z5kyD4UnLa4jvY/fdmGvaokOYVJULq/Fj4D1V9VqSQ8C3gAOTXqyqjgJHAVZWVoZZVWnBxjZz7YNlqGmLDuEUcO26+9cAp9cvUFW/rarXutvHgT1JLmmwbkmdMc9cF2XZatqiQ3gCOJDkOuBXwO3Ap9YvkOQK4KWqqiQHWQuiVxqsW1p6yzBznbdlrenMgVBVbyS5G3gEuAi4v6qeSXJX9/x9wCeAzyV5A/g9cHuNuarSLhvLp1r6xJpC+ryRKysrtbq6uuhhSL2x2cx13759vPrqq/MbzEiMrRtIcqKqVnbyb/2mstRz076PbRhMb9nODUzLaxlJPWU30N7YuoHWDASpZzxoteW5gen5lpHUA+ffvpgUBvv27VvatzBmsdk1hazpZHYI0gLZDbRlNzAbOwRpzuwG2rMbaMMOQZoTu4H2rGlbdgjSLrIbaM+a7h47BGkXOHNtz5ruPjsEqRFnru1Z0/myQ5Bm5My1PWu6GHYI0g44c23Pmi6eHYK0Dc5c27Om/WGHIG3BmWt71rSf7BCkDThzbc+a9psdgrSOM9f2rOlw2CFIOHPdDdZ0eOwQtLScubZnTYfNDkFLx5lre9Z0HOwQtBScubZnTcfHDkGj5sy1PWs6XnYIGp3NZq6wvH9AfRZ2A8vBDkGj4cy1PWu6XOwQNGh2A+3ZDSwvOwQNkjPX9qyp7BA0GHYD7dkNaD07BPWeM9f2rKkmadIhJLk1yXNJTia5Z8LzSfLl7vmfJLmhxXo1Xn3tBvbu3fumse3du3fuY9ipvtZU/TFzICS5CLgXuA24HrgjyfUXLHYbcKD7OQJ8ddb1apymOWAt6qC1d+9eXn/99Tc99vrrr/c+FPpcU/VLiw7hIHCyqp6vqnPAQ8DhC5Y5DDxYax4H9ie5ssG6NQJDmbleGAZbPb5IQ6mp+qVFIFwNvLDu/qnuse0uA0CSI0lWk6yeOXOmwfDUV85c27OmmkWLQJj0X9+F/8VNs8zag1VHq2qlqlYuvfTSmQenfnHm2p41VSstPmV0Crh23f1rgNM7WEYjNoZPtezZs2fi20N79uxZwGjGUVP1S4sO4QngQJLrkuwFbgeOXbDMMeDT3aeNbgLOVtWLDdatHhvbzPXcuXNvOfjv2bOHc+fOzW0MY6up+mXmDqGq3khyN/AIcBFwf1U9k+Su7vn7gOPAIeAk8DvgM7OuV/015pnrPA/+6425puqPJl9Mq6rjrB301z9237rbBXy+xbrUT5sdsMCD1k5YU82b31TWTJy5tmdNtSgGgrbNmWt71lR9YCBoas5c27Om6hMDQZty5tqeNVVfGQiayJlre9ZUfWcg6E+cubZnTTUkBoKcue4Ca6ohMhCWlDPX9qyphs5AWDLOXNuzphoLA2EJOHNtz5pqjAyEEXPm2p411ZgZCCPjzLU9a6plYSCMhDPX9qyplo2BMGDOXNuzplpmBsIAOXNtz5pKBsJgOHNtz5pKb2Yg9Jwz1/asqTSZgdBDzlzbs6bS1gyEHnHm2p41laZnICyYM9f2rKm0MwbCgjhzbc+aSrMxEObImWt71lRqx0CYA2eu7VlTqT0DYZc4c23Pmkq7y0BozJlre9ZUmg8DoQFnru1ZU2n+DIQZOHNtz5pKi/O2Wf5xkncleTTJL7rf79xguV8m+WmSJ5OszrLORdu/fz9JNjxwVZUHrm06X09rKi3WTIEA3AP8sKoOAD/s7m/kg1X1V1W1MuM6F+L8Aevs2bNvee78AcuD1vZMEwLWVJqfWQPhMPBAd/sB4KMzvl6v2A20Zzcg9desgXB5Vb0I0P2+bIPlCvh+khNJjmz2gkmOJFlNsnrmzJkZh7czdgPt2Q1I/bflSeUkPwCumPDUF7axng9U1ekklwGPJvl5VT02acGqOgocBVhZWZnbEWL//v0TA2DduOY1lNGwptKwbBkIVfWhjZ5L8lKSK6vqxSRXAi9v8Bqnu98vJ3kYOAhMDIR581Mt7VlTaZhmfcvoGHBnd/tO4NsXLpDk4iTvOH8b+DDw9IzrnYnnBtqzptLwzfo9hH8Avpnks8D/AZ8ESHIV8C9VdQi4HHi4O1C8HfjXqvrejOvdEWeu7VlTaTxmCoSqegX42wmPnwYOdbefB/5ylvXMwvex27Om0jiN9pvKzlzbs6bSuM16DqFXfB+7PWsqLY9RdAjOXNuzptLyGWyH4My1PWsqLbfBdQjOXNuzppJgIB2CM9f2rKmkC6XP/9Mn2XBwfR53n9kNSOOW5MROryo9iA5hPWeu22c3IGkagziH4MFqZ+wGJG1HrzuEG2+80QPXNtkNSNqpQXQI2prdgKRZ9bpD0ObsBiS1ZIcwQHYDknaDHcJA2A1I2m12CD23WTewb98+Xn311fkNRtKoGQg95N8bkLQIBkKP2A1IWiQDoQc8SSypDzypvCDnTxBPCoN9+/Z5kljS3NkhzJndgKS+skOYA7sBSUNgh7CL7AYkDYkdQmN2A5KGyg6hEbsBSUNnhzADuwFJY2KHsAN2A5LGyA5hSnYDksZupkBI8skkzyT5Y5IN/6hzkluTPJfkZJJ7ZlnnvE1zhVEvKSFpDGbtEJ4GPg48ttECSS4C7gVuA64H7khy/Yzr3VV2A5KW0UznEKrqWdj8PXXgIHCyqp7vln0IOAz8bJZ17wbPDUhaZvM4h3A18MK6+6e6x3rBbkCS1mzZIST5AXDFhKe+UFXfnmIdk6bdGx5hkxwBjgC8+93vnuLld8ZuQJLebMtAqKoPzbiOU8C16+5fA5zeZH1HgaMAKysrTY/M/r0BSdrYPL6H8ARwIMl1wK+A24FPzWG9f2I3IElbm/Vjpx9Lcgq4GfhOkke6x69Kchygqt4A7gYeAZ4FvllVz8w27KnG5rkBSdqGWT9l9DDw8ITHTwOH1t0/DhyfZV3TshuQpJ0ZxaUrtvjYq0EgSVMYdCDYDUhSO4MLBLsBSdodgwkEuwFJ2l3p88E0yaaD6/PYJWkRkpyoqg0vNrqZwXQI5xkCkrQ7BhMIBoEk7a5e/4GcG2+80S+QSdKc9DoQJEnzYyBIkgADQZLUMRAkSYCBIEnqGAiSJMBAkCR1DARJEmAgSJI6BoIkCTAQJEkdA0GSBBgIkqSOgSBJAgwESVLHQJAkAQaCJKljIEiSAANBktQxECRJwIyBkOSTSZ5J8sckK5ss98skP03yZJLVWdYpSdodb5/x3z8NfBz42hTLfrCqfj3j+iRJu2SmQKiqZwGStBmNJGlhZu0QplXA95MU8LWqOrrRgkmOAEe6u39I8vQ8BrgAlwBj7pjcvmFz+4brz3f6D7cMhCQ/AK6Y8NQXqurbU67nA1V1OsllwKNJfl5Vj01asAuLo926V6tqw3MTQzbmbQO3b+jcvuGa5TztloFQVR/a6Yuve43T3e+XkzwMHAQmBoIkaTF2/WOnSS5O8o7zt4EPs3YyWpLUI7N+7PRjSU4BNwPfSfJI9/hVSY53i10O/EeSp4D/Ar5TVd+bchUbnmsYgTFvG7h9Q+f2DdeOty1V1XIgkqSB8pvKkiTAQJAkdXoTCGO/DMY2tu/WJM8lOZnknnmOcRZJ3pXk0SS/6H6/c4PlBrX/ttofWfPl7vmfJLlhEePciSm27ZYkZ7t99WSSLy5inDuV5P4kL2/0XaaB77uttm1n+66qevED/AVrX6j4d2Blk+V+CVyy6PHuxvYBFwH/A7wX2As8BVy/6LFPuX3/BNzT3b4H+Meh779p9gdwCPguEOAm4D8XPe6G23YL8G+LHusM2/g3wA3A0xs8P8h9N+W27Wjf9aZDqKpnq+q5RY9jt0y5fQeBk1X1fFWdAx4CDu/+6Jo4DDzQ3X4A+OjihtLMNPvjMPBgrXkc2J/kynkPdAeG/N/aVGrty6+/2WSRoe67abZtR3oTCNtw/jIYJ7rLXIzJ1cAL6+6f6h4bgsur6kWA7vdlGyw3pP03zf4Y6j6bdtw3J3kqyXeTvH8+Q5uboe67aW17383rWkbA/C+DMW8Ntm/SVQJ787ngzbZvGy/T2/03wTT7o9f7bBPTjPvHwHuq6rUkh4BvAQd2e2BzNNR9N40d7bu5BkKN/DIYDbbvFHDtuvvXAKdnfM1mNtu+JC8lubKqXuza7pc3eI3e7r8Jptkfvd5nm9hy3FX123W3jyf55ySX1HguYz/Ufbelne67Qb1ltASXwXgCOJDkuiR7gduBYwse07SOAXd2t+8E3tIRDXD/TbM/jgGf7j6xchNw9vxbZz235bYluSJZu7Z9koOsHS9emftId89Q992WdrzvFn22fN1Z8Y+xlth/AF4CHukevwo43t1+L2ufhngKeIa1t2IWPvZW29fdPwT8N2ufABnS9v0Z8EPgF93vd41h/03aH8BdwF3d7QD3ds//lE0+Ide3nym27e5uPz0FPA789aLHvM3t+wbwIvB69//eZ0e077bath3tOy9dIUkCBvaWkSRp9xgIkiTAQJAkdQwESRJgIEiSOgaCJAkwECRJnf8H7KASxGaawDgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "filenames = []\n",
    "\n",
    "for i in range(N):\n",
    "    print(i)\n",
    "    for j in range(numofmol):\n",
    "\n",
    "        new_y = traj[1,:i,j]\n",
    "        new_x = traj[0,:i,j]\n",
    "\n",
    "        if i == 0:\n",
    "            ax.scatter(traj[0,0,j],traj[1,0,j],color='black')\n",
    "        #elif i==N-1:\n",
    "            #ax.scatter(traj[0,i,j],traj[1,i,j])\n",
    "        if where[j] == 1:\n",
    "            ax.plot(new_x,new_y,'deepskyblue')\n",
    "        else:\n",
    "            ax.plot(new_x,new_y,'lightskyblue')\n",
    "        ax.plot(line1[:,0],line1[:,1],color='black')\n",
    "        ax.plot(line2[:,0],line2[:,1],color='black')\n",
    "        plt.xlim((-1.5, 1.5))   # set the xlim to left, right\n",
    "        plt.ylim((-1.5, 1.5)) \n",
    "\n",
    "\n",
    "    plt.savefig(f'test_{i}.png')\n",
    "    filenames.append(f'test_{i}.png')\n",
    "    #ax.clear()\n",
    "    #display(fig)\n",
    "    clear_output(wait=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stick Cross"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "def brownian_singlestep_cross(x0, n, dt, delta, radius, vradius, theta, out=None):\n",
    "    x0 = np.asarray(x0)\n",
    "    # For each element of x0, generate a sample of n numbers from a\n",
    "    # normal distribution.\n",
    "    coordhold = np.copy(x0)\n",
    "    coordfinal = np.zeros((2,n))\n",
    "    #coordfinal = []\n",
    "\n",
    "    c, s = np.cos(theta), np.sin(theta)\n",
    "    rot = np.array(((c, -s), (s, c)))\n",
    "    # Quadrant 1 check\n",
    "    if x0[0] >  radius and x0[1] > vradius: \n",
    "        where=0\n",
    "        for i in range(n):\n",
    "            print(i)\n",
    "            print(\"Q1 For Loop\")\n",
    "            if i == 0:\n",
    "                coordfinal[:,0] = x0\n",
    "                print(\"Q1 If Statement\")\n",
    "            else:\n",
    "                print(\"Q1 Else Statement\")\n",
    "                r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                if (r[0] + coordhold[0]) > radius and (r[1] + coordhold[1]) > vradius:\n",
    "                    print(\"Q1 Else If Statement\")\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                else:\n",
    "                    print(\"Q1 Else Else Statement\")\n",
    "                    while (r[0] + coordhold[0]) < radius or (r[1] + coordhold[1]) < vradius:\n",
    "                        print(\"Q1 Else Else While Statement\")\n",
    "                        r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                    print(\"Q1 Else Else Post-While Statement\")\n",
    "    # Quadrant 2 check\n",
    "    elif x0[0] < -radius and x0[1] > vradius:\n",
    "        where = 0\n",
    "        for i in range(n):\n",
    "            print(i)\n",
    "            print(\"Quadrant 2, Iteration:\", i)\n",
    "            if i == 0:\n",
    "                coordfinal[:,0] = x0\n",
    "            else:\n",
    "                r = norm.rvs(size=(2,), scale=delta*sqrt(dt))\n",
    "                if (r[0] + coordhold[0]) < -radius and (r[1] + coordhold[1]) > vradius:\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                else:\n",
    "                    print(\"Entering Quadrant 2 While Loop\")\n",
    "                    while (r[0] + coordhold[0]) > -radius or (r[1] + coordhold[1]) < vradius:\n",
    "                        r = norm.rvs(size=(2,), scale=delta*sqrt(dt))\n",
    "                        print(\"Inside Quadrant 2 While Loop\")\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "    # Quadrant 3 check\n",
    "    elif x0[0] < -radius and x0[1] < -vradius:\n",
    "        where=0\n",
    "        for i in range(n):\n",
    "            print(i)\n",
    "            print(\"Q3 For Loop\")\n",
    "            if i == 0:\n",
    "                coordfinal[:,0] = x0\n",
    "                print(\"Q3 If Statement\")\n",
    "            else:\n",
    "                print(\"Q3 Else Statement\")\n",
    "                r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                if (r[0] + coordhold[0]) < -radius and (r[1] + coordhold[1]) < -vradius:\n",
    "                    print(\"Q3 Else If Statement\")\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                else:\n",
    "                    print(\"Q3 Else Else Statement\")\n",
    "                    while (r[0] + coordhold[0]) > -radius or (r[1] + coordhold[1]) > -vradius:\n",
    "                        print(\"Q3 Else Else While Statement\")\n",
    "                        r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                    print(\"Q3 Else Else Post-While Statement\")\n",
    "    # Quadrant 4 check\n",
    "    elif x0[0] > radius and x0[1] < -vradius:\n",
    "        where=0\n",
    "        for i in range(n):\n",
    "            print(\"Q4 For Loop\")\n",
    "            print(i)\n",
    "            if i == 0:\n",
    "                coordfinal[:,0] = x0\n",
    "                print(\"Q4 If Statement\")\n",
    "            else:\n",
    "                print(\"Q4 Else Statement\")\n",
    "                r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                if (r[0] + coordhold[0]) > radius and (r[1] + coordhold[1]) < -vradius:\n",
    "                    print(\"Q4 Else If Statement\")\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                else:\n",
    "                    print(\"Q4 Else Else Statement\")\n",
    "                    while (r[0] + coordhold[0]) < radius or (r[1] + coordhold[1]) > -vradius:\n",
    "                        print(\"Q4 Else Else While Statement\")\n",
    "                        r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                    print(\"Q4 Else Else Post-While Statement\")\n",
    "\n",
    "    # horizontal stick check\n",
    "    # Does not need to check is x0 is above or below the stick as it is already checked in the quadrant checks\n",
    "    elif x0[1] < vradius and x0[1] > -vradius: \n",
    "        where=1\n",
    "        for i in range(n):\n",
    "            print(i)\n",
    "            print(\"Horizontal Stick For Loop\")\n",
    "            if i == 0:\n",
    "                coordfinal[:,0] = x0\n",
    "            else:\n",
    "                r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                if (r[1] + coordhold[1]) < vradius and (r[1] + coordhold[1]) > -vradius:\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                else:\n",
    "                    while (r[1] + coordhold[1]) > vradius or (r[1] + coordhold[1]) < -vradius:\n",
    "                        r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "    # else case should be the vertical stick check\n",
    "    elif x0[0] > -radius and x0[0] < radius:\n",
    "        where=0\n",
    "        where = 1\n",
    "        for i in range(n):\n",
    "            print(i)\n",
    "            print(\"vertical stick case\")\n",
    "            if i == 0:\n",
    "                print(\"vstick if\")\n",
    "                coordfinal[:,0] = x0\n",
    "            else:\n",
    "                r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                print(\"vstick else\")\n",
    "                if (r[0] + coordhold[0]) > -radius and (r[0] + coordhold[0]) < radius:\n",
    "                    print(\"vstick else > if\")\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "                else:\n",
    "                    print(\"vstick else > else\")\n",
    "                    while (r[0] + coordhold[0]) < -radius or (r[0] + coordhold[0]) > radius:\n",
    "                        print(\"vstick else > else > while\")\n",
    "                        r = norm.rvs(size=(2,),scale=delta*sqrt(dt))\n",
    "                    coordhold += r\n",
    "                    coordfinal[:,i] = coordhold\n",
    "    else:\n",
    "        print(\"conditionals failed\")\n",
    "        print(x0)\n",
    "\n",
    "    coordfinal = np.squeeze(np.array(coordfinal)) \n",
    "    # Add the initial condition.\n",
    "    #coordfinal += np.expand_dims(x0, axis=-1)\n",
    "    coordfinal = np.transpose(rot) @ coordfinal\n",
    "    initialrot = np.transpose(rot) @ x0\n",
    "    \n",
    "    return coordfinal,where,initialrot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "finished generating\n"
     ]
    }
   ],
   "source": [
    "# What's probably more efficient is to have a counter of how many particles need to be in each area, and then reduce a counter as particles are placed in each area.\n",
    "# Then, if the counter is 0, then the particles can't be placed in that area.\n",
    "# Two Pairs (CROSS)\n",
    "# In vertical stick, generate initial position/coordinate in top area only.\n",
    "initialinx = numpy.random.uniform(low=-radius, high=radius, size=(round(numofmol*ratiointra*0.25)*2,1))\n",
    "initialiny = numpy.random.uniform(low=vradius, high=1.5, size=(round(numofmol*ratiointra*0.25),1))\n",
    "initialiny2= numpy.random.uniform(low=-1.5, high=-vradius, size=(round(numofmol*ratiointra*0.25),1))\n",
    "initialiny = np.concatenate((initialiny,initialiny2),axis=0)\n",
    "initialinVStick = np.concatenate((initialinx,initialiny),axis=1)\n",
    "\n",
    "initialinx = numpy.random.uniform(low=-1.5,high=-radius,size=(round(numofmol*ratiointra*0.25),1))\n",
    "initialinx2= numpy.random.uniform(low=radius,high=1.5,size=(round(numofmol*ratiointra*0.25),1))\n",
    "initialinx = np.concatenate((initialinx,initialinx2),axis=0)\n",
    "initialiny = numpy.random.uniform(low=-vradius,high=vradius,size=(round(numofmol*ratiointra*0.25)*2,1))\n",
    "initialinHStick = np.concatenate((initialinx,initialiny),axis=1)\n",
    "\n",
    "initialQ1x = numpy.random.uniform(low=radius,high=1.5,size=(round(numofmol*(1-ratiointra)*0.25),1))\n",
    "initialQ1y = numpy.random.uniform(low=vradius,high=1.5,size=(round(numofmol*(1-ratiointra)*0.25),1))\n",
    "initialQ1 = np.concatenate((initialQ1x,initialQ1y),axis=1)\n",
    "\n",
    "initialQ2x = numpy.random.uniform(low=-1.5,high=-radius,size=(round(numofmol*(1-ratiointra)*0.25),1))\n",
    "initialQ2y = numpy.random.uniform(low=vradius,high=1.5,size=(round(numofmol*(1-ratiointra)*0.25),1))\n",
    "initialQ2 = np.concatenate((initialQ2x,initialQ2y),axis=1)\n",
    "\n",
    "initialQ3x = numpy.random.uniform(low=-1.5,high=-radius,size=(round(numofmol*(1-ratiointra)*0.25),1))\n",
    "initialQ3y = numpy.random.uniform(low=-1.5,high=-vradius,size=(round(numofmol*(1-ratiointra)*0.25),1))\n",
    "initialQ3 = np.concatenate((initialQ3x,initialQ3y),axis=1)\n",
    "\n",
    "initialQ4x = numpy.random.uniform(low=radius,high=1.5,size=(round(numofmol*(1-ratiointra)*0.25),1))\n",
    "initialQ4y = numpy.random.uniform(low=-1.5,high=-vradius,size=(round(numofmol*(1-ratiointra)*0.25),1))\n",
    "initialQ4 = np.concatenate((initialQ4x,initialQ4y),axis=1)\n",
    "\n",
    "initialQcomb = np.concatenate((initialQ1,initialQ2,initialQ3,initialQ4),axis=0)\n",
    "\n",
    "initial = np.concatenate((initialinVStick,initialinHStick,initialQcomb),axis=0)\n",
    "print(\"finished generating\")\n",
    "# There appear to be no issues with the initial positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mol # 0\n",
      "[0.18337839 0.60482748]\n",
      "0\n",
      "vertical stick case\n",
      "vstick if\n",
      "1\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "2\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "3\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "4\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "5\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "6\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "7\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "8\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "9\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "10\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "11\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "12\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "13\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "14\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "15\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "16\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "17\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "18\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "19\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "20\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "21\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "22\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "23\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "24\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "25\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "26\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "27\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "28\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "29\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "mol # 1\n",
      "[ 0.1595574  -1.17884967]\n",
      "0\n",
      "vertical stick case\n",
      "vstick if\n",
      "1\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "2\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "3\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "4\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "5\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "6\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "7\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "8\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "9\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "10\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "11\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "12\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "13\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "14\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "15\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "16\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "vstick else > else > while\n",
      "17\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "18\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "19\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "20\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "21\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "22\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "23\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "24\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "25\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "26\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "27\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "28\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > if\n",
      "29\n",
      "vertical stick case\n",
      "vstick else\n",
      "vstick else > else\n",
      "vstick else > else > while\n",
      "mol # 2\n",
      "[-1.03814726  0.18423003]\n",
      "0\n",
      "Horizontal Stick For Loop\n",
      "1\n",
      "Horizontal Stick For Loop\n",
      "2\n",
      "Horizontal Stick For Loop\n",
      "3\n",
      "Horizontal Stick For Loop\n",
      "4\n",
      "Horizontal Stick For Loop\n",
      "5\n",
      "Horizontal Stick For Loop\n",
      "6\n",
      "Horizontal Stick For Loop\n",
      "7\n",
      "Horizontal Stick For Loop\n",
      "8\n",
      "Horizontal Stick For Loop\n",
      "9\n",
      "Horizontal Stick For Loop\n",
      "10\n",
      "Horizontal Stick For Loop\n",
      "11\n",
      "Horizontal Stick For Loop\n",
      "12\n",
      "Horizontal Stick For Loop\n",
      "13\n",
      "Horizontal Stick For Loop\n",
      "14\n",
      "Horizontal Stick For Loop\n",
      "15\n",
      "Horizontal Stick For Loop\n",
      "16\n",
      "Horizontal Stick For Loop\n",
      "17\n",
      "Horizontal Stick For Loop\n",
      "18\n",
      "Horizontal Stick For Loop\n",
      "19\n",
      "Horizontal Stick For Loop\n",
      "20\n",
      "Horizontal Stick For Loop\n",
      "21\n",
      "Horizontal Stick For Loop\n",
      "22\n",
      "Horizontal Stick For Loop\n",
      "23\n",
      "Horizontal Stick For Loop\n",
      "24\n",
      "Horizontal Stick For Loop\n",
      "25\n",
      "Horizontal Stick For Loop\n",
      "26\n",
      "Horizontal Stick For Loop\n",
      "27\n",
      "Horizontal Stick For Loop\n",
      "28\n",
      "Horizontal Stick For Loop\n",
      "29\n",
      "Horizontal Stick For Loop\n",
      "mol # 3\n",
      "[0.3454049  0.08070741]\n",
      "0\n",
      "Horizontal Stick For Loop\n",
      "1\n",
      "Horizontal Stick For Loop\n",
      "2\n",
      "Horizontal Stick For Loop\n",
      "3\n",
      "Horizontal Stick For Loop\n",
      "4\n",
      "Horizontal Stick For Loop\n",
      "5\n",
      "Horizontal Stick For Loop\n",
      "6\n",
      "Horizontal Stick For Loop\n",
      "7\n",
      "Horizontal Stick For Loop\n",
      "8\n",
      "Horizontal Stick For Loop\n",
      "9\n",
      "Horizontal Stick For Loop\n",
      "10\n",
      "Horizontal Stick For Loop\n",
      "11\n",
      "Horizontal Stick For Loop\n",
      "12\n",
      "Horizontal Stick For Loop\n",
      "13\n",
      "Horizontal Stick For Loop\n",
      "14\n",
      "Horizontal Stick For Loop\n",
      "15\n",
      "Horizontal Stick For Loop\n",
      "16\n",
      "Horizontal Stick For Loop\n",
      "17\n",
      "Horizontal Stick For Loop\n",
      "18\n",
      "Horizontal Stick For Loop\n",
      "19\n",
      "Horizontal Stick For Loop\n",
      "20\n",
      "Horizontal Stick For Loop\n",
      "21\n",
      "Horizontal Stick For Loop\n",
      "22\n",
      "Horizontal Stick For Loop\n",
      "23\n",
      "Horizontal Stick For Loop\n",
      "24\n",
      "Horizontal Stick For Loop\n",
      "25\n",
      "Horizontal Stick For Loop\n",
      "26\n",
      "Horizontal Stick For Loop\n",
      "27\n",
      "Horizontal Stick For Loop\n",
      "28\n",
      "Horizontal Stick For Loop\n",
      "29\n",
      "Horizontal Stick For Loop\n",
      "mol # 4\n",
      "[0.25131477 0.35250245]\n",
      "0\n",
      "Q1 For Loop\n",
      "Q1 If Statement\n",
      "1\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "2\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "3\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "4\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "5\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "6\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "7\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "8\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "9\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "10\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "11\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else Else Statement\n",
      "Q1 Else Else While Statement\n",
      "Q1 Else Else Post-While Statement\n",
      "12\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else Else Statement\n",
      "Q1 Else Else While Statement\n",
      "Q1 Else Else Post-While Statement\n",
      "13\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else Else Statement\n",
      "Q1 Else Else While Statement\n",
      "Q1 Else Else While Statement\n",
      "Q1 Else Else While Statement\n",
      "Q1 Else Else Post-While Statement\n",
      "14\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "15\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "16\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "17\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "18\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "19\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "20\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "21\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "22\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "23\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "24\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "25\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "26\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "27\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "28\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "29\n",
      "Q1 For Loop\n",
      "Q1 Else Statement\n",
      "Q1 Else If Statement\n",
      "mol # 5\n",
      "[-0.34898986  0.45274132]\n",
      "0\n",
      "Quadrant 2, Iteration: 0\n",
      "1\n",
      "Quadrant 2, Iteration: 1\n",
      "2\n",
      "Quadrant 2, Iteration: 2\n",
      "Entering Quadrant 2 While Loop\n",
      "Inside Quadrant 2 While Loop\n",
      "Inside Quadrant 2 While Loop\n",
      "Inside Quadrant 2 While Loop\n",
      "Inside Quadrant 2 While Loop\n",
      "3\n",
      "Quadrant 2, Iteration: 3\n",
      "4\n",
      "Quadrant 2, Iteration: 4\n",
      "5\n",
      "Quadrant 2, Iteration: 5\n",
      "6\n",
      "Quadrant 2, Iteration: 6\n",
      "7\n",
      "Quadrant 2, Iteration: 7\n",
      "8\n",
      "Quadrant 2, Iteration: 8\n",
      "9\n",
      "Quadrant 2, Iteration: 9\n",
      "10\n",
      "Quadrant 2, Iteration: 10\n",
      "11\n",
      "Quadrant 2, Iteration: 11\n",
      "12\n",
      "Quadrant 2, Iteration: 12\n",
      "13\n",
      "Quadrant 2, Iteration: 13\n",
      "14\n",
      "Quadrant 2, Iteration: 14\n",
      "Entering Quadrant 2 While Loop\n",
      "Inside Quadrant 2 While Loop\n",
      "Inside Quadrant 2 While Loop\n",
      "15\n",
      "Quadrant 2, Iteration: 15\n",
      "16\n",
      "Quadrant 2, Iteration: 16\n",
      "17\n",
      "Quadrant 2, Iteration: 17\n",
      "18\n",
      "Quadrant 2, Iteration: 18\n",
      "19\n",
      "Quadrant 2, Iteration: 19\n",
      "20\n",
      "Quadrant 2, Iteration: 20\n",
      "21\n",
      "Quadrant 2, Iteration: 21\n",
      "22\n",
      "Quadrant 2, Iteration: 22\n",
      "23\n",
      "Quadrant 2, Iteration: 23\n",
      "24\n",
      "Quadrant 2, Iteration: 24\n",
      "Entering Quadrant 2 While Loop\n",
      "Inside Quadrant 2 While Loop\n",
      "25\n",
      "Quadrant 2, Iteration: 25\n",
      "26\n",
      "Quadrant 2, Iteration: 26\n",
      "27\n",
      "Quadrant 2, Iteration: 27\n",
      "28\n",
      "Quadrant 2, Iteration: 28\n",
      "29\n",
      "Quadrant 2, Iteration: 29\n",
      "mol # 6\n",
      "[-0.22471978 -0.77210639]\n",
      "0\n",
      "Q3 For Loop\n",
      "Q3 If Statement\n",
      "1\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "2\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "3\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else Else Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else Post-While Statement\n",
      "4\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "5\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else Else Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else Post-While Statement\n",
      "6\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "7\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "8\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "9\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "10\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "11\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else Else Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else Post-While Statement\n",
      "12\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else Else Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else Post-While Statement\n",
      "13\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "14\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else Else Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else Post-While Statement\n",
      "15\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "16\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "17\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "18\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "19\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "20\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "21\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "22\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else Else Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else Post-While Statement\n",
      "23\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "24\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else Else Statement\n",
      "Q3 Else Else While Statement\n",
      "Q3 Else Else Post-While Statement\n",
      "25\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "26\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "27\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "28\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "29\n",
      "Q3 For Loop\n",
      "Q3 Else Statement\n",
      "Q3 Else If Statement\n",
      "mol # 7\n",
      "[ 1.06923158 -0.28633121]\n",
      "Q4 For Loop\n",
      "0\n",
      "Q4 If Statement\n",
      "Q4 For Loop\n",
      "1\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "2\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "3\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "4\n",
      "Q4 Else Statement\n",
      "Q4 Else Else Statement\n",
      "Q4 Else Else While Statement\n",
      "Q4 Else Else Post-While Statement\n",
      "Q4 For Loop\n",
      "5\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "6\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "7\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "8\n",
      "Q4 Else Statement\n",
      "Q4 Else Else Statement\n",
      "Q4 Else Else While Statement\n",
      "Q4 Else Else Post-While Statement\n",
      "Q4 For Loop\n",
      "9\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "10\n",
      "Q4 Else Statement\n",
      "Q4 Else Else Statement\n",
      "Q4 Else Else While Statement\n",
      "Q4 Else Else While Statement\n",
      "Q4 Else Else Post-While Statement\n",
      "Q4 For Loop\n",
      "11\n",
      "Q4 Else Statement\n",
      "Q4 Else Else Statement\n",
      "Q4 Else Else While Statement\n",
      "Q4 Else Else Post-While Statement\n",
      "Q4 For Loop\n",
      "12\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "13\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "14\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "15\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "16\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "17\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "18\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "19\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "20\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "21\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "22\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "23\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "24\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "25\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "26\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "27\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "28\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n",
      "Q4 For Loop\n",
      "29\n",
      "Q4 Else Statement\n",
      "Q4 Else If Statement\n"
     ]
    }
   ],
   "source": [
    "where = np.zeros(numofmol) # initialize where to a default value\n",
    "for i in range(numofmol):\n",
    "    print(\"mol #\",i)\n",
    "    print(initial[i,:])\n",
    "    [traj[:,:,i],where[i],initialrot[i,:]] = brownian_singlestep_cross(initial[i,:], N, dt, delta,radius,vradius,theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_horiz_line(theta,vradius):\n",
    "    \n",
    "    x = np.linspace(-2,2,100)\n",
    "    \n",
    "    ylow = np.zeros((100,))\n",
    "    yhigh = np.zeros((100,))\n",
    "    ylow[:] = -vradius\n",
    "    yhigh[:] = vradius\n",
    "    \n",
    "    line3=np.c_[x,ylow]\n",
    "    line4=np.c_[x,yhigh]\n",
    "    \n",
    "    c, s = np.cos(theta), np.sin(theta)\n",
    "    rot = np.array(((c, -s), (s, c)))\n",
    "    \n",
    "    line3rot = line3 @ rot\n",
    "    line4rot = line4 @ rot\n",
    "    \n",
    "    return line3rot, line4rot\n",
    "line3, line4 = make_horiz_line(theta,vradius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABf30lEQVR4nO2dd3wcxfn/33OnU7Ut2XKRmyz3bmxsDKYEETohECCEYkogxCEEEgghceJf+peElhAIEDAEAsGEQEIvoQtsMLgg915luUi2qtVXuvn9sXu6vbu9XnR3mvfrpZe2zO7O3N7NZ56ZZ54RUkoUCoVCobD1dAYUCoVCkRwoQVAoFAoFoARBoVAoFAZKEBQKhUIBKEFQKBQKhYESBIVCoVAAMRIEIcSTQohqIcQGP+dLhRANQog1xt+vYvFchUKhUMSOjBjd5x/AQ8AzAdIslVKeH6PnKRQKhSLGxMRCkFJ+AtTG4l4KhUKh6BliZSGEwjwhxFrgAPATKeVGq0RCiAXAAoC8vLzZkyZNSmAWA7N69erubbvdzsyZM3suMzHAXB6A2bNn91BOFApFrFi9evURKeWgSK4VsQpdIYQoAd6QUk6zONcPcEopm4QQ5wEPSCnHB7vnnDlz5KpVq2KSv1hQUFBAQ0ND9346hP0QQnjsp0OZFIrejBBitZRyTiTXJsTLSErZKKVsMrbfAhxCiIGJeHYsqa+vJz8/v3vfuzJNRbwFIB3KpFAoIiMhgiCEKBJGTSOEmGs8tyYRz4419fX1HvvpUIEqUVAoFBCjMQQhxL+AUmCgEKIS+DXgAJBSPgp8E/i+EKITaAUulyncNyGl9Kg0hRAp39WSjmVSKBThERNBkFJeEeT8Q+huqWlDOlag6VgmhUIROmqmchSkY1dLOpZJoVCEhhKEKEnHCjQdy6RQKIKjBCEGpGMFmo5lUiQnTin5fL/W09lQoAQhZqRjBepdpoKCgp7JiCKtuWdNJ2XVcFe5EoWeRglCDElHUTDPu2hoaFCioIgrShR6FiUIMSbdRMF7Mp55prZCEQuOK/Tcv6tc48gRJQw9gRKEOJCOomAm1cujSC5OL3b4HHtiH7y9R4lColGCECfSTRTSrTyK5CTLtL22Dh5Zp0QhkShBiCPpVommW3kUyYPrm9QOLJzlthgau9S4QiJRghBn0q0STbfyKJKDG4o9982iAEoUEoUShASQbpVoupVH0fMUFroFYEONXvkvnOUg05RGiUL8UYKQINKtEk238iiShzcq3Ns/nuWgONe9rzyQ4osShASSbpVoupVH0bOMzLE+fuVEB+eNdO8/sQ/K9ilRiAdKEBJMulWi6VYeRc8xf5Kv+ynAR/s03tnneezzI/D4BiUKsSaRayorDLzDTBcUFPj4+qcSKmy2IhY0N7sreH/jBYOzoLpd367R4J5yjZ/OshYSRfgoQegh8vPzu2f9psPsXyUKinCobtb4z25oDKORPzATrp+iV/4uwXAa295eSYrIUILQQ9TX11NQUNAtBulQgSpRUFixtlrjgwPQEcFXIS8DbpnuW9kvnOXg3nKNLmNfiUJsUILQg9TX16ddBapEoffS0qLx1n7Y1aS33KPlu8We7qje3DHLwVObNKqMLqS7yjVuGAkDByphiBQlCD1MOlag6VgmhScVjRqv74WjnaFfMygTjnSA1TfBBlwzHor6OLq7gx6vgIWFFolNXDfFwWf7NT6p1vef2AdndGrMKep9ohALhw4lCElAOlag6Vim3srbezU21NLdPROMTBscNxBOGe6gqUnjse2gAYc7fNN+czSMK/CsvPva4WioDwNOHO5gQpbGE4Yn0vsHYXOdxtWT018UYu3VJ5L5Rzpnzhy5atWqns5GwvB+ucn8bkIlHcuUrjil5KnNnRxuD/2a/plwyWgYmOtZ+S6t1Pj0sPU1Zw2FY4O04F1WQrjjAt7eSek6rhBECFZLKedEcl9lISQR6diqTscypSv3ren02/dvB6YNgHNHhVbBeovB7AFwZojXmnl1p8aFY0O/buEsh4co3FWucf1wGDw49YUhmDXg+l1FYzWoiWlJRjpO9ErHMqUj5470f64L2NcEra2h+YleXOLe/sG4yMQAYHNj+Nd4WwVP7ocXt6XuJDYhhN/fjJSy+y8WKEFIQtKxAlXrMyc/0wc6WDhL/7t4FDi8vna1HfDAFr3VfVe5xrsV/ivZCf3dlfLDO8LPS+mQ8K9xYTWpbWdzagXHc4lAMCGINUoQkpR0FAW1PnPqMGGAg9tnugVicj/fNF/WuMXh/rUa+496Vrg/GOfe3lYXXmV8wrDou3iKsuCCEZ7Hkl0UEmkNWKEEIYlJN1FQ6zOnLheO1YXhnBHW59ud8M8dboH4xxaNDNMI5Ut7In/2fWFU4g+Y0n57ioMpgxxJv7ZCT1kDVihBSHLSURTMpHp5ehszB7kr16vGwmA/EUoPterdS2aWhhmh1KUnoUx1OFSjUVen0ernvJUofFnds8LQ09aAZZ5i8UAhxJPA+UC1lHKaxXkBPACcB7QA35ZSfhnsvr3N7TQQsXLfdJS5f2CHJsOQKPpqo0G5o6YuD67TaOnS5xv8+Bi9oq3t0PjX1tAmquXa4arxMCAneLdQKO6n9fUaj+4O/txsoS/V2Wr6qmWir7mQKEL1FIryGRG7ncbKQvgHcE6A8+cC442/BcDfYvTcXkOsLAXz77VoM4gysJVFnK2ISTfLpzfx7Yn6/w6Tj+qATAc/mK53z1wxDrIC1CwtXbDYGJyubQqtlf5llf90oYgBQJv0FAOADhLThZSM1oAVMREEKeUnQG2AJBcCz0idz4ECIcTQWDy7NxGLSlSWWhxDFwZRBtllEWQsQqzKU1JSwpIlSxKXCUXY9Mtyt6g/rPStTEf1dXDbMbo4nD0CMgJ8TQvyQpsK9e6B0PN33TDdovheCYzNC62Si4coJNPYQKgkagxhOGBe4qLSOKYIk1iJwr8mWp9rxy0OQ5eGfeuwefbZZz329+7dy4IFC5QoJIj1+zXuKdc4dCi8CnF8X/3/SmMCWrOmsblW46P9Gi/v0nhmq8ZjGzWWHQSbALufr+k9azoDVoqu5/jDqiIfMkQXrP79HVw6wcFPZzlCmoF7V7lGdQzGFVLFGrAiUTOVrT4dy09ECLEAvVuJ4uLieOYpZYnF7N/Lh+p/2WW6CFhxqEsXBoB5efDZcRFlNyCLFi2C/xyCD/4Ff7sNgJaWFhYtWsT8+fNj/0CFB28aQeH+cRA46FkZ2tAr8k7jq2X1DZNE17q2i8CNmkvGOcK6v3kZzh21Gv/Za51uwSgYMEAXjkfWaTQasZOe3A+F1RrfnRbeuEIixgYSQaIEoRIwz4McAVgagVLKxcBi0AeV45+11CRWISHaSuHQIRi6JXC65c1ucbi9CO6bFPajLKmoqIB+A+Bbt8IXb8CXHwC6paDoWZyAM8SvlN2wAhwCMu2QkwH9MmFAJgzK0VuEr5heaY4dfjQjvEr3kXUaN5muucdCKOZPcrB4g0athYYUZMCNFmsr3DTDQXW1xpP79f2aMPQtkBCkigiYSZQgvAbcLIR4HjgeaJBSHkzQs9OWWIlCURHIIpi+HDZ4mQuHJuuDz2b+dEj/A3hhElxaFPYjuykuLmavqwz3vANnuL+SKu5RfOjs7OTjQ7CpLnafbZfU/zqAkj7w9THuire+TeNRr+9QaxjRTF00el1jFXfJypr45igYNyCw+Awe7GDhYH3Ow+lBRjfTxRqwIlZup/8CSoGBQBXwa8ABIKV81HA7fQjdE6kFuE5KGdSfVLmdhob5C5qfnx/1+swuS8DFIBtUfwVePATfCmBJROLGumTJEq4quhzsdtA0OCvTJ00q/8CSgc7OTlbVQPkRSaOf9QgALhktGF9g3UZ8fafGxgjiCnljDj5ng5DWQ66p0Xi8wn09wF/WaLQF+FrY0RfQiSWpYg1E43aqwl+nAealOCE2orBoG/zBq1OvfBLMNKyBX2+D3wXw/HCeCqGOd4sPnWCzQVsLo26cwp133slVV13lkSaZv6fJiFNK/r29k33Nvi1pARRkwtzBglmDMjxa1aGGi/78gEZZVWzyOiwbLhwG+fn+n+3Ko2sVNX/jCjP7wzklsROCVLQGlCAofEQhZtEPyzxblAJwlnqmOXs1vHvU+vpMoL3U+lz3Pcv0/9OzYN0803E1eS1ipJTcvcY966RPBhxTKJg3GDIyPK2ASATBRXu7xv2bosurPwod8K2RulAEG1i+bpjbuygWpIo1YEUyTExT9DDxCgnhLNUtAxeuOQuLtrmPvTNbd2WVpVDs9ZvswO3GOviTwM96yyv8spq8FjlCCC4YJbrnADR1QoZN+IgBRLeITJZpTsK0/vq9zh4W8e08qNHgb7tC82KKhRik4ryBWKMEIY2IVwU6s0iv7AeZvi1/OOA71gCw9yS3OOR5nTvsdIvDnC98rx1hEThNiULkTBmQwU9mOhjXT//MPj7o5OENGk0d/kd0/xSBC+mMAfr/DXX6/9U17nNnjYB+cXZdGd8nuutTed5ArFErpqUZ8VyhrPorxj3L3MdEGUzLgvXzfNM3lbq37WWefdmrW60FxQq16lp0fHNsBgebu3huh5OjGjy00cm8wZJTh/v+/COZUXDeKAfrDD/P/+7QOGJ4qs0ogGMHOTh2kDvtllrNw/3UhbeV8uI2jZ3NoT3/kvHhWwepODaQCNQYQpoS7/73S9bAS/Xu/QmZsPXE4NdVVfm6sZp5dAx8z898RDWmED2v7ensdjfNzYBrJ9rIz7Rz6JCmT04jsi6k+9Zo3RPYAIblwjUTg9/nyY0aeQ64bELwtIcOaTx70DMeVxZwWxj5TeWxgVBRg8oKSxJRgZpb+VZxkkK5zh9rJ8MMLzdWJQrRc6S1i2e3O2kzeo5mFQrOLnZ7G/W1wQ+OCU8UdjVovLBL3+6TATdbTADrKXqbNaAGlRWWJKL/vTCCW3qLgSyFX1lEtjpms97V5JFWjSlEzcAcO7fOcDDLeHnlNZK/rHN3Fh21mvEVhDH5DjJtkGdPHjFQYwPho8YQ0px497+vnwDDtoaW9pfb4P8s5i4Eshas1l/xLlNBQUHU8y56I2cXZ3BCURdPb3XSEsoqNEH4cZhWRTzobdZArFEWQi8gnq3qoSEEMd9QrVf6VmJgRSawfopuOZgHps2Yy6TWZ46c/Ew7P5zu4ITBnt+JJzcm1zKTwVDWQGxQFkIvIVGeOocPw6iN+F3K0Bs78O4U+Org8J+Vn5/fPRlPrc8cHaXDM5gzqIuHNur9RdUd8J9dnVxQYiPTlpztRmUNxB4lCL2IeIjCZWvd26EMFP93Mlwco2U76+vrlTtqDOmTaefrRU5eNwIX7miQPLCui1OKJCcU2Xs2cyZ6g6dQT6EEoZcRjShsOgzHbAxt0XNvLsiHV2dFcGEQ1ByF2DJ1qIPXjcVyBmfBkXYoO+hk1WEnF5TYKe7bM9aCsgYSgxKEXkg4leh5X8LbYUS5lKW+lkLleBgex/XxlCjEh+p2uHW6nVf2ONlzVPLcji5G5nVx0Rg7uRmJEQZlDSSW5OwcVMSdUAearcQgA9g41R2i4sES033KvJ5TGl8x6H6OckeNGScOdG9nZ9i4fFwGV0+w088B+5rhr+u7+GB/F05nBP6pIaBiCvUcShB6MaFUoju9VkY7PA20UphiCkdwS4nvvSc4wpuoFguUKMSGr4x0u49+aSyrOTzPxk3THJw+3IZNwMpqJw9u6GJ7fexEQXkK9TxKEHo5wSrRMUWweZp7f9AG33tYdRFtPSlGGQwTJQqx5d1DnvvHDbZz2ww7EwsEbV3w391dPLlFo6EjMmEoKChQ1kASoQRBEbQSnTTQUxS8g9uZeXZCYrqIAqFEIXpG5/o/l2GzcdHoDL472U5hFlS3wt82dvHm3s6Qu5FcImDlLqysgcgIJKyhogRBAYQmCku9RMHKzfTqbb7HegIlCtFxmSkw3cGD1pPUCrNtfHeKg6+PspNpg/W1kvvXd7HuiHV47XhYAx9WQ58y9/dRlMEpK8O6RUoTbLwlXJSXkaKbYCEhTjZE4RSvbqMnx8H1O4x7JCCfoaK8j2LD04dgYYAZ6VMH2JhcAO9WOllbI3lrn5PlVU6+OcZOYY4tKk+h6mqYvgWqw+iRWhZi2OxUJl4NHCUICg/MlagrJIRZFLzFAOC6EXBvBWzuSFAmw0CJQni8vF1ja1P41wkhmF5oo1M62VArOdLaxTFzTmT/Jutoxa53UF0NpTthc3s0uXbjADpKY3OvZMN7mVxvXJ9pNGKhBEHhQ7ghIUSZ9fyDZEGJgn+e2qhRFaGQSyk51CLZXC/ZUuekUYNdqz/hhf/3bRqqKj0TZ+XA/1q6dyP5rgjg4gJ4ox689SMLaCsN/56pQCLnYihBUPgQKCSEv4o/WcXAhRIFnS/2a3xUHTxdHwE3z/SNXiqlpKoVttQ52VzvpKEDOtrbeO6nl7P9s3dwdpnGD+x2+O1LcNIFMcm7BP5b75VP4GhpTG6fVIRqDcQaJQgKSwJVoDOyYV2bfvysvvDu0Z7IYfj0ZlF4ZI1GY4CiDs+Cq6dYh6+WUnK4zS0Cde16a33pkgf5YPHvaTtqUXF9FN/PNR+oL43rI3qEnp6ZrQRB4Rd/FejaE9wWwbtH4Y4iuNfkrz5iKVSekti8hkpvEwXXKmjeTOsL548LvH7B4VbJlnonW+qc1BgisHPlR7x6z21U79xkeY2UMmJrMQc4vg/8fSyMX+u5BreLQTb32t7pQk9ZA1YoQVAExF8FenAaDDUGmO89BDcPhoeMroj91l6HSUO6i8L+/Rr/9NMtFGy95Jo2yeY6J1vqnRwxrMDiPoK/3PxNNn/8BtLp+3IPtUuGZOrbmw57nssCZubBc+NgTP/gee90guMT3+PFDtjbQ5Md40VPWwNWKEFQBMWqAh36iecX9lGvCsg10JysJKMoLFmyhEWLFlFRUUFxcTF33nkn8+fPD/n6Zfs1llkIgQ34aQAhqGvXRWBznZPDhgiMyBO8c/8dLHvhMTrbrVe3uGS95KUjUPQZjM+BiwbCPfs808zrAx+FuLqvlWUxPhO2nRja9alAMlkDVoiezkAg5syZI1etsnZbUyQenxZNCP3EySwK4Fumnvo9LFmyhAULFtDS4vbEyc3NZfHixUFFwd/4wOhczwlmZurb9e6gzXVOqoz6fnieYFKB4LQJRTQ3HLG8TkrJbdvhwf3uLp3CDKgJEBM92HfA0kkBcAa5LpVIpDUghFgtpQxRhr2uVYKgCAfvL/awpZIDQRZIUKIQnJKSEkafdQ1nfm8RFRtX89i1ev/IqFGj2LNnj+U1/sYHrh4Mw4f7CkFDh+4euqVecrBFL+PQXMHk/oJLTp5B5XbrcQHQP5O/7IOf7YIO4+M5oS+8dwz0yYCWTshbZn3tvH7w7SL9L9MUG8HfWEO+gPpT/WYlZegpa6DHBUEIcQ7wAPqKiE9IKe/yOl8KvArsNg69JKX8XbD7KkFITrwr0H4fSYItmZBKopCfn+8xGS8R2Gw2/m9lKzabja6uLv7fcTnd+TLHBwp3fKCxQ7K13snmOskBQwSKcgST+gsmFdjon+0/eo2rbnjpMHx7Mxw1sjE+B96fAcU57rSTl8MWY3KALIV3auCc9Z73E8DobNjV5v9zSAcx6OmxgWgEIeoxBCGEHXgYOBOoBFYKIV6TUno3N5ZKKc+P9nmKnse7/73xNBG0+0iUwXG5sGJunDMXIT29PvMvy6q7xzFcYgBQXFwM+B8fsAN3eAlBk+byDpJUNuvvZXAOnDrUxqT+NkYP6R9Sy3VFA3x9A1QbhsgQB7w5HWb3873GJQauCuXsQve5346AfS3wRG36ikGyjw2ESiwGlecCO6SUuwCEEM8DFwL+7U9FyuMtCoQgCitbdGEosEFdkrkO1tfXe/yoEz3InNOnHxLo6HDPwc3NzeXGf3xm2TXkPT7QrBmWQL1kX5Oe70HZcMpQG5MKbBRmBw6AZi7r7lY4cy3sNCrvfnb45yS4YJD1tZll7m2t1Pf8ryt9jw20Q02XO/ZV3xQVg562BmJNLARhOGD2LagEjrdIN08IsRY4APxESrnR6mZCiAXAAnC3jhTJiY8onN8f3qjzSecAzFVavVMXBhvQVRrfPIZDoBna8eSucg2EQEjJU5dPRwjBrz9rJCszE7wqHO/xgfU1TjbUOqlokkigMAtOKrIxucDGwBwRVsu1oQPOWAerjFhGWQLuHQu3jAicf9e7HW6qTaqqrNPeXgT3TdJFxCwGjSkkBuliDVgRC0GwkkjvT+RLYJSUskkIcR7wCjDe6mZSysXAYtDHEGKQP0Uc8RCF5npYdDHc+ZJHGuuhT91LxTWwuH0GjBsQp0yGQaLdUddXuj+dIZmCG1/ebpnO3/yBHY1OGjXJvCF6d9Cg7OBx8b3L09UFF2yAtwwttwO3jYB7xwXPv3lguPJkWFcFx2z2TXdvMfxkjL6dWeb+TqSSGKSbNWBFLAShEhhp2h+BbgV0I6VsNG2/JYR4RAgxUEpp7dumSCk8KtHPXoZHfgI33YcsBXuZ9YxTb8av0///tQRuLolDJsMgkaLwpmkiV5WXclqND3jztWI7Dpuex0hart/bCk8cNMQZ+NYgWDJJD0MUjMOmvF+QHzie1XcMsa/qcItBHskvBulsDVgRiwVyVgLjhRCjhRCZwOXAa+YEQogiYfzChBBzjefWxODZvY5k/QJ65OvFP+ljCuhdQifmWV8jgP9M9Dx2yx69YrloTezzGA6JWGDnYT9uo6NzdYsgmBgAZNoFNpst7NXH/rgXMj+GxYYYfCUfmk+B56eGJgYAg02dvq95PfqFSZ6eZSVr9f+uGc0AvwjSFdWT9NYV3aK2EKSUnUKIm4F30Bs1T0opNwohbjTOPwp8E/i+EKITaAUul+n4acaZ7XUa/90Dg7LgO34CkfUk/lrWnx4H22pgopcbogS+tVWvOGproXCd+9wr9bowjMiAfScnIPMWxNtS8I4J6G/+gBWRtlyfPQTf2wYthtk2NRfePwaKskJ6bDc3Wo4AwruT4cwh7v18GzQ4sZw4t6gSfhFCt1Si6G3WgBVqYloKse+oxhJjZbK8DLhlevKJAgSe6GXVrWAHOkvd+xll4B0xJwPdg6WuDn64D96qhzpn8BXa6mZCQUFI2fZLvCauVVZq/Psw3B6CJeAvL2YCzZ8oq4NLNkKtMYlwRCa8PQOm9QknxzrnfQlve008WTsZZgzxTVtfD/3X6Nsui6GmBgYajYPnJ8Blw8LPQyxJt7GBHp+YFi+UIPjy3FaNCiO6gUPA7RYx65OBQJXo8GX4zG62Wulq4CdQE8bSif6IxaS4npzNHE3LdXOTPkGswvBm7Z8B/54CZ0Y4gO8tBv6EwIyrEWB+D+aGQU9MWkxnayAaQYjFGIIigVw50UFfo6NPk3CPn37onsb7B1VgaqbvP1lvGZrRgKwyz2NHvqJXFqcGaMXagaIMWDRUT+v6K5/pTlPkJ6RCOCRiTMGbQP3Y+fn5AfuxD3fAMSthyipdDHJs8MQEqD05dmIgCC4GZk5Y4d4+Mt29vS2Bo4m9dWwgVJSFkKLcu0ajy/TqbpsCWVnJZy0ECwnh3YWUA7SUxubZJZ/CXkMvPzoGSkMIvxyMRFgK0XRhdHTBuevgQ6O+yxDwi5Hw2zHR5cmqmwigv4DaAJ5C39kIT5q8kXrCSkhna8AKZSH0Qu7w6iq6fxO0tyeftZCfn9+9bdkqK4Vs034r0KcsNs/eY4qff9ra2NwzXpaCq+VqdT9/1kCdBnkfw2926XMJrt4E2Ut1MRDAdUOg7eT4iQFAndQrdtefrcxz3ywG3nxmshJq4mAlRGNh9VaUIKQwN0/13L9/ExxuTi5RqK+vp28/tyhYVXitpfB9U1iEZqBfWWyeb255ZsTqnjEUhUCTyFwVlr+B4h3N0CLhtxXgWArPVuuD7GcVQNsp8OTk0F1I/eEtBq4uufv8uIz6q16/2tfXCphninc00MsDLVIKCgqi+kx7O0oQUpg+mQ4u82r9/X0b7GxIHlG4q1zjF2WHGXXJZd3HrH6sj0z17Fc+ChwNElY7VC41uoq6gB9vic09oxGFSKwBK8xrEEhgfDbUnAjvzITMCIRgcx18bGqpW1kGrpb/TyziEwXiw6OwxWIa6uNRWi8ulDUQG9QYQhqwdL/Gp16RMM8cBrOH9PyYgjkw279/cDFrlr/Rve/vu+dyO41lv7KtzN16jeV9wxlTiJV7Y2cnHLMaNlksZFZzIgzI9D0eDPPnkwu0BEhrZnIWLJ8Cpp5BthyBKRsCuwSXmlZSi2YsId1cRmOBGkPo5Zwy3MFoL0+c9w7A+xXJYykAXPbwS0w44xvd+/5+zJ2lsR9kNK++Feki8FYEsxRiZQ24uGM7ZC5zi8ENXl4+hZ9BR0do91rfoIuvKPOsvP2JQbED6md5enNtmucpBgCTBuqftyvNtGzvO0FZk9vaCJdYf6YKN0oQ0oTLxrvdUV2sqoGXdkYmCtXNGv/YrLGjLnJROXDQ99rr7nmBcSed172fCPdNF/8yhVOc9UXs7mslCvHox+50wn373ZX3qCxYYDGpK+cz32MVDVC41HPAd0a57wRAn7yWuv/2nuRb+YfC+hPc93h9auC0tjL/59TYQPxRXUZphrc7KkBRDnx7knX30eoqjWVV0BqgZvAXaTMYfynX8LceyhM/vJCdy97u3k/U97Dfx3DUeFQsZjGbiXQWcah0SbhgHfyvzjpgoA338WwBrafCyavg06bQ7n9uP98B5Hgz+BM47FUYD9fUOH+m6YiaqazwwGpBFZcpGO7E3/F94ZJx4QvC39ZpNARpfj5x8wXs/Ox/3fuJ+i7G2v+9J/qxj3TAhRvg80b/79R7HYpAfD0fXjeNx/bE7OF2J2hOfY1mNTYQOUoQFN18XKmxPIDvd6hk2+DWY0ITgi3VGq/sD+2+5lYswJM3X8j2zxJrKZjj60Qajz/Urq54tGIbO6Gv3b12ToMGBZ/6Ty+AGwbD437WYj45D5Y1u/d7av1rZQ3Ehh5dU1mRXIQiBrdMgLw8vbK/p1zzaWHa0FuXmw5pTCnyFYUnN2pUhzhwCTArH8qN1uc3RsKEgQ7uK9foBK5/6FWeuuUbbPv0LSAxq5QVFMDMHFjTqncfPb8fLh8e2rWhtFzNM2NjvT6zlJBvEYrDjns8YGy2e/lLF2YxKMqAQ4bL6ihHz4uBsgaSB2UhpBmu7qLjB+ktyM+NimB6f1hvrIglgBw7tAQbUYyAc4bCTENEXHlZOMvh0Y3lGpP4a7mGqy566kffYNvSt7rTJOJ7GWrXUSQtV3HVb2HJb7r3sz6S2ARk2CATyLXpf3kZesC5Sdlw7wT9vQSiS0LGx4HT+KOPDb5V6J49nCeg2fQxJ1IMlDUQP1SXkaIbV8U7Kg9KB0JZDew1BhW/PhJe3+d7jSB4GOloWDjLwZEjGk8Yzz5vJMwYqIvCUxs1qgxr4+lbL2bLJ8HnKcQSlygIPF1TIbqWqyhD75u6qL/rZvBh8BGcluMhJydoMr66Cj4KcbDYCu+uO7MYdEn9fDwcwJQ1EH/UPARFN64G5sEWeHqvLgZ9jI5BKzGA0MXg+uF65b5wloOJfcP78gwc6O56esuUj+umOphozKG49i8vMenU87vPJcIl9eeG26YEvrU2dj7un89A75t6bpfxAAmnCXIFZKH31WagC5H5SblfQKvFhDNvPpyjV+IlEXb6WonBxE91Icv4GGwfe7qojl5qfZ8lS5ZQUlKCzWajpKSEJUuW+KRR8wZSB2UhpBlPbdKoag+cpjATaoKMAYzIgav8uKpGitlKuG4EDBnkvn9ZhcbnRtiERFsKGWXQ9cuLYdnLluejeX5OGbStWwY/OiXg/VpbdTFwEaqlALCpEaZ+6Xt8sA2qvqJv+5sAlgc0W58KyqX7P+bNBefR0uKeypabm8vixYuZP3++sgZ6CNVl1EtxSomUYLe5f3i76zT+vSfwdQ4BVxbB0wfdx64fDoMHxzfUxbYajZcq9O1vDoFxwzyft+GAxhtV+vbTP7qILUvf7D6X0DDTWbnkZzti1o99wyb4+zPPwN3Xdh+zKk9tKxSGKQqVjTDSQgz8cXQO5K8K7H78xXiYawyyz1gO6/01MMxlkBKuHQWV/oMcqbGBxKC8jHop/9jSSXUbOGwwIg/mDbExur8D9vh6n988HsiEhzbqC+u8WqPHyu+U+v94iwHQLQbgKwYA04Y5GIDGM1Vw7QMv89QtX2fbp+8AsfU+8isE446Fx1cD0AAs3gsLRkX/vCemgPNb1/DUgCHws3P0PNgdyC7P9zQgB2qOd4tC7heBB3obLcTgrUlwXoAAfn39tK8qx8NwC0+rdfO80lVCyQ6LGc7fm+1XDMJ9b38u1+gA5hTCGcU9H4+rN6HGEFKUhjZdDECfzLP7KDy3w+kzKe3b4/V+/z59HPTJdHBJiX68vgNGGn33nRLa2uIb9+g+U76+W+w/3bBhDq4x4vNc99fXGX/Cmd3nYhFmOmA/9vbV/NoU1vl7uyG3LOJHevBUIzD3bPj1i/oBZycir4Avaj3TuUShO98Bnp9vEoM+Nl08zi3SZymDezzJCjueYSmsxMCKESNMsaa+anP/7VjjmbcoxgZcvZm7/KzBoIgfqssoCqSUdDpld6gIJ7rlLKWxbdqXxr7TlVaaznsdd90LwGlaSN5pus8be7vQnDAgC44fYmNznZMDzdDh1RcggL4OmDpAcHKRHbtN8NF+jS+8Jinl2OFHM+LTGqup0XjcsA5yBfwwxHWgXeL2+IIz2bXK7WsZznc20n5s74p43SyYHkEcH4D+n0C9+b2c5gCMiQBFJRz32m5WHOd5zf5WGGHqPrKyFKZ+Dpva4K2JcO5Q9/F3jsA5G9z74x2wXYN8oN7iPuES6DM1jyFEgpV7siI8VJdRD7G9wclLu2OwCnwU1LbD2xX+8yCBRg2WV0lG5jkZk2/ntOEODrZoVJjcFgPFMoqWx01dRaGKgXmt6O8ufo9fzutHZ7tuEgXrPgpmSYQiKLIUvl4ObxjzymaUw3AHVJ4U8DJL6r1ej5SaO4+H9rDy6rnYHl3h4fY6PAcqj3eLgijzFYWNJ1g/76lD7u1MYFsEefYm2GcqhKC4uJg777wzYjGoqUmu6Ly9ESUIUTAoRzC2n6C1U9JhxGHp6NJb6d4B5vxhEzA+X+AwOu+E+U8YLonC89jqwxIJ9HPA5P627uOu6z+r0mugTAHzimx0OiWdEkb3c/cQXjnewV/XazSbFln5+yaN70yJbavskXWhdRW5sIrDBKC1tXpUSgUFBb7rM8fYq+X1WdDQAAXl+v5+Ta+Y62eFHvXTX5ePlNKd360rkT8/H/HHN3hwDNxifE6hiIKZu3fDwr2ex9YfG1o+/ZFITyFzw0HRMyhBiIL+WTYuHWs9DNPplDR3QrPm+g/NndLjf4vx//xRdhy2wC0wFxVHO1llzDS9aZp15e0SBAnMK/Lfk3zLdIdH6IrDQdxVI6HRsDxsQGGhf7HxJwTmboP8/HyfkBCxsAYCkZ+vV8ITl8M24/MpKIcT8+DT4wJeypOmCu5f4+GK7b55687/52/CH6/lhz9/mjsr4NDJ+uHhObB9Loxfoe9bicIj++AHO32fH+nM456YRbzdyzo4Y6ifhIq4osYQUow/rdXQnDB7oODMkb56vqxSY5kpnlGwftiODo0/b3TvRxrd1IpQ+oNDEQIz5jhB/ojXd7qiAUaVez2r1H96f9ZBFtBmus6jAr7oh/DDBwBoOg7y8vTDO1rcouB67j/2w/XbrScWHp0DffpYnAhAT84b8P4eqPGDyFFjCL2ErXWdaE69W8hKDAAPMQiFzEwHP5iq8bAhCtuPRpdHF89vc//Az7foKrISAjtwR4CKIN7WQDCKDWuh4BNoMMwqUQbz8uCz4+CJCvjeruAhxr0NMQ9L4eUHIcMBN91Hn5UwfyA8Ow3G5cK6mTBjjfu5/nhvcuhikAwxhT6pVGMHyYJyO00hXturV3inDIltSIe+mQ4uLnHvL94Y/Q90j2n66zRTV9Fd5ZqPGGShtwj9iUGglbJcJNLSrf8KPDTSvb+8Wa+gvxtADMyrSP5nou95j/y/+Ce4/xYAlhyBrDL98PQCXRTMTMyGsVnu/QVD4AyvZTWtSKbVxz7zasSoVmrPERNBEEKcI4TYKoTYIYRYaHFeCCEeNM6vE0JEOdTV+yg/3NkddOzEYbH/yUzo7yDXGG6o7YCVVZGLglVXkZUQ5Bnnb7MQglCWofROHy9qavT5CObYPjf7iQvljcvPv8D0yi7x0z/uUabXHoI7zgV0v3xRBl/U6qKwdy6c2Fe/70l9YadhckzKhscm+89LqEt7JpIlW3y/Zz9R3UU9RtQ1ixDCDjwMnAlUAiuFEK9JKTeZkp0LjDf+jgf+ZvxXhMh7lfoP9cwR/iu+cotK3FUJ93fA+aNgeF//P7YfznCHqf7gAAzP0xjWJ7wfp9n8n15g3TVUYIcb/cx5CKcf26OrhehnM0sJtbUwaSMcCdGbOAt4axScvtf33G1F7u1Dnb7nrfNgKtOq/8FpfbrDmp6wDo7JhTVz4dPZ8NohdyjrHAGb/bihJnNMoX0hBPJTJI5YNDXnAjuklLsAhBDPAxcCZkG4EHhG6t++z4UQBUKIoVLKg763U3iz7EAnTsAuYNYg/6/s3QPubfOCKQB1GvxzB5gXVXQImDEAzjSFB+ibAUeNyuuZ7fDjqRqZmaGLgtn8X1/veW5oJlw71doaCESgSitWouCUYA+yzoAdWDIBLrNY2F6OhgnLYbtpgOCfh+HPk/TtObmwqkW3ioIx7lPJjlMywNkFNMNpAttHEiewtkVfiL5xDlxoClHR4rXqW6jzBpYsWRLxvIFoeXCNGjtINmIhCMMBswFdiW/r3yrNcMBHEIQQC4AFAMXFITiu9wI+rdIruAtGBak4jf9T+sEFY/WKt7ZV4+U9cKTN1xtFk7C6BlabXP68n/DnjbBwVmj59Ocx5MJbDGLVco2FKHjP8Aa4cwT8Ylzo99g2Dz7dBycbLqBHuvSunvenwsq5rrwGvse45bCzA/igE05zl8l5muDEVZLPjurv0RyTKNRF6Z999lkWLFhAS0sLUkr27t3LggULAHpEFFqMz8K89vPUCGeDK2JDLATB6hvo/bUPJY1+UMrFwGLQ3U6jy1rq896+TiSQaYOJ/UN7XS4xABiQ4+A7Fv3KH1ZqrKnxrQitPvA/rdG43WKGcUODxnvVsKspuGdNf+PyeA0ORysK2XY4PA0GmUI+/KEyPEEAuK/W99gZG/UunZZTAy86M3ulezwAgI+khyhsPL2A9ZX1TDeJwdE5oX+mJSUltLS08H8rW7HZbNQcqOBPXx/PokWLEi4Id5saD7ebVtT7+pjoxw9WHNL40NTUVC6soRMLQagETD4XjAAORJBGYcHqI/qP+ZtjAv/ot9SGZ35/dYSDrxqB3GpbNA61wa4G2NKgB7sD90pqmgze+g/EwlkOhBDc6Od8rPqxoxWFgQP11rbLpbOZ4LODvfnC5LZrvler1LcfHg03WURQPaccvjQ8swrtUGP093lPxjt5RAGyvp5TV8MncwR9/eTDqtwVFfpMOZtN9yUpHFbMH1a3c7SuJvQCxghX7gZlwv1RfLfM3Fuu+UZhVYRF1BPThBAZwDbgdGA/sBK4Ukq50ZTma8DNwHno3UkPSinnBrt3b5+Y9spujS311oHnNE3jYCtUtcCRDlhr+k3n2vVKvcsUaK8n+MVxucgu/6Op8RrQ9G4xR/KcvDJoMe0fnqYLRjCyytzROl1Ccu0GeOaIO40dPWKoi6s26O6loAf/az7VLSSvTIZrJwafjOciUFlLSkrYu1cf/f7t5004HO7vlOszS0Rr2tsLzbux8Z0RMGhQ+PmwarT0RuugxxfIEUKcB/wF/bv+pJTyTiHEjQBSykeF/m17CDgH/Xd2nZQyaE3fmwVBSsnda0J0TUkAmQIGZMKhEMJb/OGcMRytjk1s/EiJhShctQ6WmLqAbh4Mf50S5Lll+n8b0FXqec5W5tkl96MiPZbV/Ub3hgPoMK5xlOnxUF3HYjHesmTJku4xBBe/WX6UrMxMn76seFWk5si30wvga6PdbsnehJOHP5VreN+hN4oBJIEgxIveLAjgv5vGFejOLiDDpv8djdLqHpgJZ42A4nzfH9G/tmjstXAPvGoQPGt4FTU31PLQlcdTf9DC/5KecW+MhSgcOeI5rtAXaCwN8Mwy/X8fAUdP9T3/h12wyCKIm7eAPH8QrnhvG9w4G1qbfC8gsvIsWbKERYsWUVFR4RGdNNwQIpESKJzJ4cMaf/dqRxTnwpUTPdOtOqDxflXg5/RWMQAlCAoC9/FnAFMK4LzR4f1IHl+vUWNhpFw1CEaM0O8lhCC7Tz5tTb5dGsnw3YqFKIBnqIgZ2bDWn8+/ke6YbFjjJw0Y6zib9i/uD/89xrhHGJPsYv0Zx1MYzEuoXlwME/wEO/zHJi0kS9QfvVkMQMUyUgA/m5kRs9m6D5VrWLVJXT807+d4i0EyCIGLWM1TMA8Qb2oLnv7bAcJHvF/luwTlS+t3IY6bDFqH54kho+D0m5FLftJ9KJaT8bxxveP7yjXMbQGXUERT2ZqXUDWLwWcHND4J0uL3x5BMqDJ9ZL1dDKJFCUKaEAsxcK1l640/ITCTTCLgTaxnNI/LCp7m1tHWx3c0wZmbTQdO8/+ZfmWV5BOLYIOxLo8VrvAR95drHsH4IhWGd/Z6Wh7ReK2B3sX27eH6WuCue51hMWFQER5KEBTcXa5Zzj9wuYv+PMC1ySwEZmJZif43zLkJLg41wfhVwKEKuG4ytLV4JsjpC//cDYWFALw+3nPdZDOJEAWgO87Uw+UaZm0KJAzv7tb4sj6y543KgSsm+d7zgTUarabiOYEn96OvWmQwZ4iyDqJFCUIvJlB/cSAhSBUR8CZWleiUQeE/u6kJhp5xHqx4B6SvI7ArH0OXueMemcXgzl2waIzvNYkQBYAfGBW/t69/pC39iX3govGBK/Aup8RuLBz1I2Ni5Pt7NVZZTP6LX2jD3oUShF6Ivx/xz4/N1P/7uS5VhcBMIitR1/2D5QdgVSMcZ2ENZKLPa/jDPl9BcF0faXme3aJxWiEM9+Pzv3+/xj+rQ7qVX6b1gw2N7v1gXU0f7NVYaarwvzUWxvRzX3PGKAdneE3sa++SZNmVJMQCJQi9iGBCYEU6iIA33pWo1frM0WItBPrcb/Nn+qsd8HuLKRuPj4UbRsL3t8Cjh9xxf6yIRBSaNUllK/yzEm50aBQUOHhivcaRCKe+HFsAZ1l4sd0fwM3URUODxt92Wd/3hZ1w01iNfv38C4kSg9ihBCEN8fY1v/GFTeDw/EH99+4bWPXvZ/zeIx2FwIx3SIhYiIIQAopGw6HdnicGjoB73ubgpdMoMlYym7sCVnoNIQigcjYMM8Wj+NskXRCCEa4o5DncaR/dDT+bKUMSgxMKobQ49L5614B0lkWd/Zc1Gm1+svj9MXSLxCM7Qw+wqIgOJQhphnk26u9XtGC32/VZqFKCEPymdCjtjdaxa9JdBMzU19d7rM8camgIKzysAbMYTD0RHvoUgPLZUNQH+n8C9V5DCHkCmiwmsYVLuKKQLeiukL1nxYdb8VthtkhvM8YAXt2hsdnPMq2lQ+CEYe5n3jRW45Gd7nspl9L4owQhzVi0aFF3REu73Y6UkqpdG3nwiuNxevu4G/QmITBTX18fcgX67fWe+37HBjIy4cV9UDC4+9CrE2H2at+YUpOyYPO8wHlcsmQJDL0CbDbyr/slj5wxKWBk0nBE4daZvnGE4lHpDsz03105MBNusFgjA6BfPwcnFGp8brRf7l+rcdsxShTiiRKENMMV0XLpc3+lfv9uvnjxUcsg/L1VBLwJtQJ90WVU3XY6Ys2Hfu8FnrOaAS7c6rn/3cGwOEhMJHBbe/z3IsjOpfHiH7Pgcj1EbaxEIT8DGkzGwf3lmuWSpuFiFoAjXu0QO7BgDORbhEnxprTYwdZGjToN2p3wwjaNb01QohAvVOiKNCOZFqNPJQKFuBBCgLBZuot6p/UOSWHmrUnwvzpYXAUhTHbWcRrPtNn07dPtjBo1ij179gS9NNSwHd6t99F5cFmEle4/t2js97Ms5gXFMMVPuIpg3FOudVtY54+AaRFEQ+0tRBO6whbrzCh6hkCLp+fm5vLss88qMQiA92fjsyC9SwzsDrjmV/CRpKrK02Moq8y/GACctwUeDEcMQBcCm617DAjcVmAwrMpkxZBsz/3dzbDhSOjzC9ZUa9xVrv95i8GoXL0bauEsR8RiAPBTk9XyRiU0NqrlN+OBEoQUxqfSsjg/atQoFi9e3GPr5qYSoQjmnE874LrfAjBkE8xfpx9v7sIj7EdBLH5ZHR3QfBQqd3ocDmdp2VBE4brJvhX1G/ugvt5/pdvQoHG3IQL/22+d5pJiuGJi7Fry5vGNR3YGSKiIGDWGkIKkakyhZCacz3RDNUzfpG8/VwsvlIFWCoPtMCYblh+nn2ts9B96whvvRXMAlix50b1+wYeGhXLW1dx5zdmh3dSU/2BjCqPzdMvAzKO7fd09H9+gUeNHJ64cA8+Z5hOMj8Ii8MdNY1GeR3FEjSGkCGpsIPZEG2bae/DYe6lNe1ng1epm58Cq4wM/1zWnZO9jWyArm7zOVprOyAmaXyvM5c3Pz/eZd+HPE+iUwbDUz4zlqfnudZBf3amx2ZiV/N1iKIyDIAB8fkCjzIiOmoE7EJ9CR62HkMYoayD2hPqZhjIo268Mj6BvffHc9+bVCXBBBFE5L1oLr9QZ+SgN/3oXgUTh9V0aG43pGIFEINcGP/Ry/1xxUOND0wS6eLfcn9qkUWXMehuWDddYdHv1VtSgcpoRbGxASqnEIEwi+UxD6X9vLIXfjXDve4uBHc8K/KJtoefZzJLxkV3nTX5+fve292Q8V0sfdDEY6BXR5Ptj9IreLAZPbdLHERIpBgDXTXF093cfaIPyKjXIHAvUGEISoayB2BPtZxqs/9219rHltaXu7TuL9aUzncBNm+GRyUEf7UFubnjp/RFsMt7xA+GLI/r2JcOgT34GDpvvZ+gd9dTF1cNjk89Q+Mks98S6dw7A2JzAMY8UwVEWQg+jrIHYE+vP1MpSkFIfQwgU/keUQbXR7fKLMeBaV+dvEa4O5uLlKK8PZPmcNtJdoT62Bw8xqK11u5eaxUDgdi0dPjixFbLyPIotShB6iFAqLCUE4RHPz9T7OptFq9mKIZt0YRBleKw85jo2oAz+F0LwOnCb8zduDy19IAKJwhlD3cfr6jTWHNZFYPFez3sMcOgV8s96eFD3prHu7WhXYuvtqEHlBKI8hWJPoj9TfwPNhw7B8TugIsLw0VZ4Dx6fUQ4fNFifixR/5QlUsR4/0NOSSAbWH9Z40wgjbsNzIltvQw0qJznKGog9PfWZ+mtZFxXB3pP1itr8582/Jkb+7NeiuNYf/srzrRLftAtG6RZBsokBwPRBDoqNcRYn8ORGZSlEghKEOKHGBmJPsnymoYaEAF0U+ptOX7EVnhjrmeaQxQDz25N8j8VqYNkbq/KM6e8gC72CcI0PDBiQfEJg5sqJDlzLPFR36PMVFOGhvIxijPIUij3J+JmGE1G09lT4y264zeiDv2EnFAqoMZIXbXantZqxbMWHR+CrAyPLuxWJXlo0XtxuCuldVgVT+mh0AouNAeeSPnB5kLWcezPKQogBydJyTSdS4TMNx1K4dbRnF1KNRdbXTg4uBq4W3FVbAyaLiHDKk8x4ex512t3n9jTp4yOPqS4lS5QgRIEaG4g9qfaZeueloKAgcPpS/z+6GUOCP2+2sQTnwTjVZ+kiCmbPoyctJgPWdejC8Oe1ShjMRCUIQogBQoj3hBDbjf/9/aTbI4RYL4RYI4RIabehgoKCpG+5phqpYA0Ewpw31/rMgegqhbOMdZPN4wc5ZcGf9WEIC+tESzqIQr9+Dr4xyvf4wlkOZhW69zucujDcU67R0qLEIVoLYSHwgZRyPPCBse+P06SUMyN1h+ppXBWW1dq7PdVyfbEaBn0K6xqjv1dbGxy3UveNz/04+vuFQqpZA4EIFBLCindm69bCkCFwjTEW0AY8sjfQVfEbWPYmHURh0gAH4/t4Hjt6VOPsYn2Q/FxTyBEn8OBWXRwO9GJhiGoeghBiK1AqpTwohBgKlEkpfZzjhBB7gDlSyiPh3L+n5yGYF2G3IpTPrroahhmhkjtOhRDnMwWlqRP6LtO3X5wC3xwcOL0VbW1wynpY1ex7LlZ+7t6k81yMUFcos8JWBq7UwT57V5TV1cfAsZY2eeyIpkzJwgd7NVbWuvcn9IWLx7nHGQ60aDxjMSZz1jA4dkjqDUD3WLRTIUS9lLLAtF8npfT5igohdgN16N/5x6SUiwPccwGwAKC4uHj23r1BmkxxwG+lNeUEeHg5ByfpfudmjvsCVvlZOtDFfSVwe0kscuiuQEL1SnHR1gYnrYcvLUQAdJOxbh70y7I+HynJ6CkUD6KpQF0VvQCcpf7T2cv0Fu2oTNhzYljZi4h0EAXviXZWk9daWjQe2uobsvyYAXDuqNQRhrgKghDifaDI4tQi4OkQBWGYlPKAEGIw8B5wi5Tyk2CZS6SF4NcasDvgj2/CcWdG/YzmUyDXHjxdME76Ej4zuolaT4HsIPcMJgIuihxw8KTo8+ciFhZWKhJKBbqtBSZ4df88shd+sFvfvmYgPD3N+v7HrIB1LcGFI5akuii4BMGO5zKnPxgHffv6Vvb3r9Vo91KG+WNhZAoEz4tGEILOQ5BSnhHgwVVCiKGmLiPLCOpSygPG/2ohxMvAXCCoICQCvy3Xkqnw1IaI73t2P/jfsRFf7peHK91i8Oyk4GIw/FOwmp9zx3D40353a2j+IHh2amzy2FusAX8E8+l/dj9cbcQjmpQDm41Fcm4aBbfv1scSnjkCT/u5/wfTYNAKdxdTIkj2eQrVrZJDLZLpAzzHpZxSYv42ju8HWTZYW6/vP7wDphdofG20Z0V/mxHi+/GNGjUdrmfAyH5xLEQSEG2X0b1AjZTyLiHEQmCAlPKnXmnyAJuU8qix/R7wOynl/4LdP14WQsCW64RZ8FiI6x4afG8QPGqqTJ0SWp36OrtNXfr/5i5odnruN5mOB0vb3AU1pjg55ookEEOWQbVXfJ0RDqg0icT/psHZUU5y6q3WQCD8tarbuiBnqWfabwyAl2cY15WZrin1ve/1W+CpQ/7Px5NEWQpSSg63dlLVAnXtUK9BkwatnbpnUEcXdEn9t+aUgVemM9MnA26e7uDoUY2Hd7iP24E7LOIfVdRrPGdYbamyXGdPjiEUAi8AxUAFcKmUslYIMQx4Qkp5nhBiDPCycUkG8JyU8s5Q7h9rQQil5frlIZi9RT82LhM+HQeDB7vSwPUb4B81vtdnA62l+pfU8XF4rTeHgDw79LFDnk3f7t43juXa4G8H9fT9M+DAvODWgZmsMs9F4EH/ETTPg6woxgt6uzUQjEAV6JTPYXObZ/rfFsNZhTCvXN+fng3rTnCfd0q4chP8+7C+P7cv3DESLhwIjgTNKkqEKPx5rUZHqLV8CLi6irzHDrzXdfDuQrq7XEOi/0Zvn6kEoUeJhSBE03Kt0+DePXD9UBhncl9bXgMnrre4V6n+/7d7IMfmrsy7K3a7dcWfGcIPefIXsKUVJufAphAsAzO3bIGHvEIsl2TB7nnh3ceFsgbCI1AFuqIejl/je01/dC8MKzLwXYdhiAO+MxS+OxRKIltyOSziLQqv79bYWK9v24G+mZBlh5wM6OuAfg4YmA0FDnh+NwzJEVwxzu6Rr531Gi8arfu+GXDU+NC8W/qv7dTYZHLdnjUAzjYGkV1jDz+aBDk5ShB6lGgEIRYtV+9F1EPBAeTbYUw2XDoIbhwIffoEvSwg398Gjx7QReboSWAPwzJoaYG8FZ7HfjgUHoggcqayBiInWAU6bjnsbCcmzM2DRaPhvAGQEUerwVwm7/WZY8FH+zv5olr/nE4aIjhlmO+Q52t7OtlaL7l+UgaF2Z6f8V/XazR3QmEWDM6GzUYbxqrrx7sLySGgKAf2tfi/JlmJ66ByKhHrluuV/eE5f800P2jAkS440gwrmuGOPaFdZwOyBRRmwLF94JYiOH0IvH5YFwOANbPDEwOAPl5isGIWHJdvndYKZQ3EhmCDsjvmwQc1cIaF5ZkjoGYu1HXCeZtgfWvgPvMVzXCh4Q/Rz6ZblMNj7EYMugi4vhuhTMYLl9OGZ2Cjk+XVkk+rJE7ZyanD3VXW7kYnm+okJxXZfMQAoNmwCC4ZB52aWxCs6NvXwcJZcE+5hhPQpFsMJqf5QLKZtBCEWLZca2pg8PrQB6kEcP1g+Lulf1XoOIEWCS0a7KuDV+uAzZ5pJq703D84C4oCVO4DlrrHMhxAR2no+VHWQOwJJgqnF+rdjkM/g0OmAZ9WCblfwHAHVJrcgl0WrCyF27bC4kP6d8hMoxNGLYfzC2H+ELiwEDJj4PoM+vrM5gZDPDyPTh2egU108mmVZHm1xEknpw3PQHNK3tnXRf8smDfE1wzq6HB7TQzIdEAm6M21wPx0loMXt2vsbHIfSwFP05iRsoIQy5ZrTQ0UrQ+8Pq7f5+BfDIJ5gKw5Cn+tgI8boVqDFieWC5f7Y2i5/2f0+RiajY8gVDFQ1kD8CcV98+CJ8K8DcKVXULb9mi4C8/rAZ6YOgSMtcP9E/c+MSzB+MhKePASv1uhdId8dCvePsx672nEUxvUNvTz19fVxd0c9ZVgGgk6WVUm+qNYtBYdNUN8Bl4+1k2Ex/f8lYz5rZgSRAS4d7+i2FAC+qIE1dVq3K2o6k3LRTmMZUyirTP/RDLQQg7m57lWvIo02EWwMYmZf+PtUvbug8Sv6jGPXM08xzNShmb6rcDUd51kGMzua9Oc2mz6CYGKQjHGa0plQ4gRdMUx/1wMsWvPLmzy/W6dZdDOZWdMIlfPgxyPALuCRA9B3Kdy6HTpNpvBvd8H41fq9J34OW5v83jLs8kTLycMyOLVIv+/Kw5LPqpxM7S8o6Wddhe018j7PakptCHj3ELQbQfCamtI7zlFKCEIsI4zmlrkXOPd2w5yS5a50v5jrPu4s9b8kYjBEGRR8orushsof9sLSRsgQsG2u7/m8PHhsnL7dAcw1upJ+tRPGe43BH7a4HlTU1p4m1Eq05hR4eKzlqW42+QmZUmiIyTsNujXwp3Fw9GS4bYTeyHlgP/RZCj/dqQvDbyrc125rg0mrIOdjuGOH5e0jKk80zBuaQelQkzVikeZws8ZdhqsowPFDwu8EeWmnu9JfOMtBiWlG+UPbYfn+9BWFpPYyEkL4zVw4+S74GBr8JC9xwO4IwjWc/yW86SfKqA3fFsZJfWBZCOP+qxrhOGNe3MfHwFcCBC+buQLWGgNfQzKgysvMOTsf/jfL85gaG0guwnHf7PsJNPkZ3LJqrBxqgqGrrM93OOH2HfDYQX0A1UWuDRpOhMu3wCtHPLswJ2TDa9NgYgCvuXi7o66rcfJWhTtXMwsF5xTrlf7z2zX2mKyaC4phSqG7m8flQhrMY8iVriQPLp+gpz1wVOMZkzDm2OFHM5KzCylt3U6tBCHU/JZ8Cnv9CHmRHQ6eElXWPAjVPdUqKJ6Z1i7IX6b/QP9fMfx+TPB7Osqsxz7McW7U2EByE04l+rtd8GtTSz5LwD/Gw+XD/Ny7TP+/fY7nXBoXbV1w7RZ44bD7fj8vhl+OApsN/lcD394CVabfUraAF6bA1wdFX55waOmUPL6pk8JswcR8+OCAft+pBbCp3u1AkWXDsr/fVdFfPwEG5/meb2vTeO8g3fMfrITDO0heMrqjpr0ghJtHq9AAuUBzaYwyZuLML+H9MNYjCCQKo5ZDRTsc1xdWzA5+r38ehGv8LKXYPBfy8pQ1kCqEW4nmfAxtpiR/HwfXj/BN5xKEvjZ9nMqK3E/0UCu5NuiQ0Cn1yZOLiuGnI3Vh6OyEy7bAq4bVMD4btp1gfb9IyhMKb+7tZGOt5LpJGQzKEXx5uJN3Kz3ve/IQOHmYdSXtXZkHwioaqotnNmscMM0wP2soHFuUPMIQjSAk9RjC7NmzI/oiZdt9+xdbcI8d7Kj1vcabd6rd6e/Z6Xnu1i3uc95iMNYYBJ6aaX3foVvgrNW+xy/fqItBPzssnxk8f6XlbjHwcMFua4Ov2v2KgRobSE7C7YNvPRWeGO/e/84O93dy2DKoNiLbXmbEqDrqp6vpnr26GAA0fwXqToL5g3Vr9ee7of+n8Od9kJEB/52me0BBcG+4WI8pVDQ5WV8rmTvYxqAcQXu7xof7PZ8xKd+/GGyt08JyDnECD6/XqLAYRL5msoNrxrn33z0ID61Lj3GFpLYQYhG64qE9cMse63N2oHoGDBjgedwq7k8g8m1Q76f1BYG7lPZPgPcFXLtVF7F9x0PZS0tYtGgRFRUVFBcXc+eddzJ//nwAOjqg/+e6iyrAyf1g6bEgcvtBezM4fX/5yfyOFZ6EO/u3pQvylgZM0o3VOINrXY3fFsOvTF2UjZ2wYCu8eFivHAsy4PclcPMIwzMvAw6fHPyZsbAUupySJ7d20umEGyZnsO5wJ+8dsE57w0RBp1PyegXUtscnIqxDwNyBcMoIR1J2IaVtl1Gsg9st2ACP+1mzLRe4ZRjc7eeLZsVbU+DcMFYqW3UAjrNY8NvFC5Oh4/0lLFiwgJaWFnjoC7j5eHJzc1m8eDFjzpvPiWtNF/zhanjvWb/3S+Z3q/BPuKKwpRkmrwyYBID/TIBLTGMNc1fByibItkGrnwZNvQbf2aoPMDuBARlQ26l3LzUHaASZiVYUPjvUxScHnVwy2sare5x0Rvi1zhQwsxC+OjJwpb20UmPFEc/B9nA4bzjMGNxzwqAEIQJKV8HHIfpZh0L5TJhZEHr6k1bAZy0WJ7q64MEfwg8fdHfenpVJ7m9foOUrl+ppHrkdXvyz5X2T+X0qQsPbCSCUd/r7nfCrfRb3ElBvXH5gHgw1QlhUtMAoI6zJluNgYl7g+x/p0IXh9Rq91Z0hQDs1hMIYRCoK+5q6WLI9/LCndmBsP71yzs6OTeVc0aTx+m53kLxA/GxmRo+tQ60EIQqOXwErrCpmEw+P1hcv6d7fDTfvtU5rA3bOhJIC97EGDZbVw8pG3b+7ol3/gW1ts76HB67309EB35sNezf6SZa871ERPpGIwuyV7lXxBO7ukjwBTV6Vd59P9DU3puTAxjCi51Z3wGlr9OiqH84KmtyDcEVBSsnda4LXvgIYmgsXjIKCGFX+4dDUpPHKPtjfpn/mgQakE4EShAhoaICC8vCuuXoAPDPD89gvtsEfw+hmCpt/3glP/QqkbyspHhEmFclDJK3q3I/12EcA1xfCk8baHb8uht8YYwQP74ObDUeJZF5g55XdnWyp9zwvgIJMOH8EDM/v+f76ZEQJQpgc7YR+y3yPz8uDz4ywELW1MGSd//hG94yEO7xmkH5nAzxpMUaRJfT+yxwb9MsAm4BtxuzSe0fDVYOgyJgNKd7v1LuNfnAC7Fxj+exkfmeK2BKJKJidGG4ogie8VldzDSQvHAl/DDILOhr2t8KtO+DDBmjoNHkmXTkWDu7qTmdVpoZ2J3/bpF/x/al28kNZNEQBKEEIGynB9rG+7R1B0oodtTB+nf/z70zXV7kyY/5R5tigYq6+oEdDBxR+pv847hsDtxebrhECMhx6rF4vlDXQewlXFMxrYAiM/n4JwxwwLhc+adAbKO1hjAFY0doKP9kLL9fAYS3M4JBnZkGn25fPu0wu751p/QXnl6RsDM4eQa2HECZChGcqjxvgTv9uDZztFUzMvF9juLGWOGCPUa+3OmHQ5/DEBPjFbl0MTi9wi4HHD95LDJJZsBWJIdwF7nNz4dNpcNIG3RLoC9QCBzQ4YAxLfDEztGfftRv+fggqOvRJa+EigD42mJarB9e71FieVmrtfsu09KB7loMSg8SiPu0wOavQLQ737oSfenl2FBqWhPmDnZUH5c1wg+FyWpgBH8wSfifK5PXLp6mhPnaZVqQ84YrCiQPh5iJ96dRaCQU2qDeGocZlw0xjHY1/H4I/7tWdHdpk+H77An1i5NhsfVbz1X5CaIA+Ixp0l1V/ZWrv7OLTQ3pGr5sYo4UbFCHTK7uM4sFla+GFOutzu6bAv1v1mZ9cPRkqt1imS+Z3oUgOwu0+Gr8cdkS5NGemgOGZsKAIFo6O7B5tbZDzub598AQoynafM5fp9Jt/xxnXL2R4nuDqCaq9GgmqyygJ+Pcx8G9je/YK+NLkyjpmiH9/ZDU2oAiHcC2F7fOCz7zPAAY54KJCuG8U5OTELLvdTDAi+NrwFAPwLNMHD/0K0d7Gu4/9X+wzoQiKshDiiAo1rYgX4VoKp38JBXZ4cAIMj0OFHwyXk8UrU+BCi9n9kcy7UFijLIQkIpAIKGtAESvCtRQ+ODYRubLmPFO4FSsxgMSsz6wIjhKEGKGsAUWiCVcUeoq3jbG1W4YGTpeI9ZkVgVGzPaLAtQSllRjk5+erUNOKuJOIpSuj4T5TiJcHJwZPn+zlSXeUhRAByhpQJBPJbCn8dLf+f06Q4HkAWifcvEMPd223JWd50h0lCCGixgYUyUwyisKKeve8hpXHBU8/7HM40gnLGpOzPL2BqLqMhBCXCiE2CiGcQgi/o9pCiHOEEFuFEDuEEAujeWai8dclBO7Vx5QYKJKBZOtuOWmN/n9ACPPLLlqviwHAphY9DlKylac3EO0YwgbgYuATfwmEEHbgYeBcYApwhRBiSpTPjStqbECRqiRLJdrW5o5ttD+IdfB8FbxiRGW92IgJdoYx4z9ZytNbiEoQpJSbpZR+lnnvZi6wQ0q5S0rZATwPXBjNc+OFsgYU6UAyVKLDjRXcMoDsbP/p6jW4YrO+/fOR8N/p+szoLa2w2livPBnK01tIhJfRcMAc8afSOJYUKGtAkY70dCVaa8Sn+3Rm4HQjluv/p+bCH4xQ3H8wwmN8zRQ0sqfL01sIKghCiPeFEBss/kJt5Vu9Ob81rBBigRBilRBi1eHDh0N8RPgoa0CR7nhXogUFBQl57nGGdSCAuQEeeeJqfdW2TAEb5rqP314MfW1QpcHLpipAiUL8CSoIUsozpJTTLP5eDfEZlcBI0/4IwO8aY1LKxVLKOVLKOYMGDQrxEaGhrAFFOrFkyRJKSkqw2WyUlJSwZMkSnzT5+fnd2+bQEPHiSAesMpbxvCdAILw/VcDyo/r2giLf809P1v9fs9nzuBKF+JKILqOVwHghxGghRCZwOfBaAp7bjbIG0p+mNo0/rdEor/ZdXCgdWbJkCQsWLGDv3r1IKdm7dy8LFizwEYX6+noPUYh3BXrlJvf2D0dap9nRDD9xL5jGk1W+aS4aBEUOaHLCPRWe55QoxI9o3U4vEkJUAvOAN4UQ7xjHhwkh3gKQUnYCNwPvAJuBF6SU1ivFxxBlDfQuNjXoq4K9sx/erUh/UVi0aBEtLS387otm/rC6nd+vaKGlpYVFixb5pPVu7MSzAn3P9KjLNoHmuxQ4k1Z67rdYpAF421i//P/t9j2nRCE+ROtl9LKUcoSUMktKOURKebZx/ICU8jxTureklBOklGOllHdGm+lAKGugdzJ3iHvB9S9r4Kkt6S0KFRV6s9lutyOEoPVog8dxbxJRgb5a7d5+aDy8cgSu3Aydpgp/7HL32sqjA3gfAczsC5NydKG/0cKXUYlC7EmLWEaBrAFAWQO9hCvGuberWuH+tekrCsXF+vqrrhm8fzh9qMdxK+JdgV5sdBeNzIQfDIf7x8J/DsPVW6BLwq3bYJexWE++HXadoK/3DPCGH/+Rj2bq/x8/CB1dvueVKMSWlBaEUERACUHvYVRfB31NwVjanXB3uUaLln7CcOedd3LsRd/1OJabm8uddwY2wONVgR5qA5chsM0ItX3rSLhnDDxfDWeuhQcMVxI7UH+Kvj0qU///011YUpQFp+br975gg3UaJQqxI+UEQVkDikB8b6pneC4JPLgBdtSnlyjMnz+fS3/+AEIInE4no0aNYvHixcyfPz/otfGoQMeu0P/n2Dwnot1RDL8bBR/Vu4+1nuTe/nWJ/n9Hm/97vz1Dd2F9pw5q/Sz9pkQhNqSMIChrQBEKGTbBwlkOsry+2f/ZDR/sSy9RIEMXvxy7nT179vgVA6vfRawrUNfA8CaLMBW/MoXAPjUf2k2PutpYI0EL8NPNscN8Y2Gdr671n06JQvQktSCsXr1aWQOKiLjtGAdzCj2PrTwCz2xNM1EAfjzL4ffcPzZr3L2mk4fX+5Y7VhXoOGO2sQBKvJbnfLTSvX16PnzcAP2Wwbwv4ZO60J/xj4n6eMPaZtjc5D9dT03GSxeSWhCsUNaAIlTOKHZw2xTPqfIHWuCBdakvCm/tDF6GZZUah4yumKZO6zSxEIWdxkDxsxN8z10+CAZlwDMT4Z1j4P9KYIgDPm+EU9fCoE/d7+c3fsYRAOx2+H/GePlZ6wLnx3synhKF0BHJXLEKIbozl8z5VCQ/T27WqPbqp/7hNMh1+G9dJzN3lbsFYaGFhdDYqPHITvf+xSUwob//snoLQai/t1/vhN/t8zxmA7IFDHTAsX3gxmFw9kDPNCsb4I5dsLTBPRgtgC9mwXH5+CX3E2h1wjsz4KwB/tOZ12cOpzzpgBBitZTS73IEgUhqC2H27NnKGlDEhOsnO7i4xPPYgxtgV0NqWwt9/Bw3iwEEFgOI3FKYbZEBJ9AioaIDXqmFczaAKPP8O75ctxIGZcAAww9AAnPLYein8Ls90GExYW2xYYVctsn3nJlETsZLJ5JaEBSKWDKhv4OFsxw4THXDC7vgo/2pKwo3W1gHZushHCIRhQuGgCz1/PtwOlxSCCVZkCusKxkJtEuo6oRaU3fWvH56ULtf74G8pXD2WthgGjO4qggKM6C+Ex7dH/vy9HaSustozpw5ctWqVT2dDUUa8vZejbW17v2CTLhxamp0H724VWNni77t3V3kTwxGZsP8yaGVL9Luo3DY2Qh/OgAfNcD+dt1LqTgLds3Txzt+uRueroI6QyxGZMGtw+G2EbC8EU5eA1kC2k5NjvIkE2nbZaRQxItzRzm4ybRuX31H5C3rROMSA2/MM7NnFHie29cGBw+GVr5EtKzH9oNHJsHm46HxK9BZqosBQJ8MuH881J4Mb03Xu6X2t+sB8bKXwl37YGimbmH8ZEdylCddUIKg6LX0y9K7kPqY5rLdVa5R05IawtDftFbx89s02o0+934ZcN5oX2vg6UOh3ztZKtFzC2HVHKg9CW4cCtk2eKMGDhoT1O6vhC6LkBbeJEt5kh0lCIpez83THZxtWsPv8a3w4vbkF4XvzdAr/fIqjT3GGgQCuGm6g8c3Wuc/HCsomSrRAgf8bSI0ngL/mQLTc/XjTuC3ewNe2k0ylSdZUYKgUACzBjs8+uN3NsE9SdiF9NQGzzw1Nmq8Y1pu6mdGGWqMFrT3jG1IXVFwcclgWDcX9p8AVwyCn/uP5+dDMpYnmVCCoFCYWDjLQabxq3CiV57b6pJHGKq8smJ2L7Waj3DDGOv7PL0ptUUBYFg2PDcVcjKCpzWTrOVJBpQgKBRe/PgYB9P6u/df2qNPbEsmhmd5tvRvGus+d/So+3jfvtaeRQfbQx9khvQLCaFEwRolCAqFQWdnJ5tqO3lhp8buo54hL6rb9FDa7e3JIQwH2t3bXxsB/fq5K/5nvFYYm2ua0Wu2IsIZZIb0CwmhRMGXMI0thSL16erqYkejZF2tpKoVWjv1BVyCIYH7N8GcQo0zihM/Z+FvJrdSV3bH94HpgzzzctTwuulreCF9dZSDFbXua68tcovBXeWaZVeTFfX19R4hIcyhIVIVKaWHELgWHOqtKEFQpDXb6ztZW6NX/C1BKn4ben/0kFyY1l8woR9kGCGm29s1HtykL/+4qgbW12ncdkxiRaHBK5RDNnDJeP95uGa077EXt2lcOsHB5GaNzUf1Y+GKQrpVoEoU3ChBUKQFFUc7WXVYcqgFmkOs+AflwPQBgommit8fWVkO7pgFL+7Q2HlUX43trnKNy8dCSb/4CENtq8aS7TBnIMwb5iADMActbcM9juCq0IONH+w03FMvHOdgS7nWbWk8uVHj+hBnaqdjBZqOZYoEJQiKlOf9fZ2sOuL747UB2Xa94p/SXzC1IHjFH4xLxzk43Kzx9236/vM7YUI/jYvHRi8Kta0a/9oBR71CVX9cBVP7aPiJYA0EdyUdkOEZMwh0F1XXddUd+iDz0KFKFFykQ5nCRcUyUqQ0VS1dPLVV70sZ2Qem9hdMK4i+4g+FB9dptBj99aF2uZhpbdV4eicEW91zYBYcMQ0im591X7l/oRiQAQumu9N6WxPex63OBSMd4wSlepmiiWWkBEGRsnR1dfHn9U66JJwwWFA6PPEGr79K1orWVo1ndkGdn3WBQfdsmlIAlc1gFZn7prGeHkVWebHi5MGwrFrfvmAoTCly3+PgQc3D42hBCQwIEi7bI88pXoFakcplikYQVJeRImV5cqsuBgOz6RExCEZrq8azu6Gm3X8aAYzvh0+Xk6tyH90Hdhvhn/s7/IsBwDXj4Bk/wd5cYgDw2kGYUuTeHzrUwYxWjXWG09DiPbCwPyGTjl0t6VimUEi+X5FCEQIfVHZS0w52AddNSI7pNK2tGs/tgcNtQZN2I4HGdqiq0hgyxF3ZC+NcRRM4AA0INmF6iUkMXBbLhxUaK2p803pbNuMG0C0IAE9s0LhhWuhWQjpWoN5lKigo8Fl4J91QgqBIOfYd7WLlYb2yuXaiDbvdHuSK+BNN6OxD7fDUAeCA7z26jD/zc64bhod4mNMClOS6j3212MFXjVg/i9drPgPL/vJ9RAtvkBnSXxRck/HSWRSUIChSio6uLv61Qx9EPnWojcE5PS8G3ozMhYuLISdHr0xdlW6eDc4ogvcO6ktMRoo/8XCxp0V/Zv8M+MZgt3gsmO72KrLhXsvYH08fgoVDw8tbOopCfn5+Wk3GC4QSBEVK8fctTpzA0FyYV9TzYnDSYL1bxywA/mh2wuQhDiYPCe3e0S7YU9fpXzyc6N1FwZ4RzqQ1F+kmCuk4Gc8fUXkZCSEuBX4DTAbmSiktXYKEEHuAo+hWbWeoI+DKy0hh5p2KTsprJA4b3DotObqKQsXcMv9piBXsveUaXejjCT+b5eCFbRq7jIll3pX03aZJZvEgkCi8uVtjfb1v2lT21LEiVcrTk15GG4CLgcdCSHualPJIlM9T9FJq2pysqZEI4DuTUksMALKAdoJ305i5aBT8Z687btG3Jrhb9N4td1ea6QXwNYvV0rrTSUmnhEfWd9JqkRnXPV/aprGt2X380XUaN85w8PA6rTtWUjDSzVJIt/JYEZV7hpRys5Rya6wyo1BY4XQ6eW57FxI4f5SdgqzUEgOA22b5ThALxrgBvhX7KYPd2w0WExUCiQHoldirO91iYE49Isf0HC+30/ouPd+hioGLdIsomm7l8SYmE9OEEGXATwJ0Ge0G6tAbMo9JKRcHuNcCYIGxOw3dCklHBgLpbDGp8qU2qnypy0QpZd9ILgzaZSSEeB8osji1SEr5aojPOUlKeUAIMRh4TwixRUr5iVVCQywWG89eFWlfWLKTzmUDVb5UR5UvdRFCRDzwGlQQpJRnRHpz0z0OGP+rhRAvA3MBS0FQKBQKRc8Q9ymeQog8IURf1zZwFunbDaRQKBQpS1SCIIS4SAhRCcwD3hRCvGMcHyaEeMtINgRYJoRYC6wA3pRS/i/ER/gda0gD0rlsoMqX6qjypS4Rly2po50qFAqFInEkR1QwhUKhUPQ4ShAUCoVCASSRIAghLhVCbBRCOIUQft3BhBB7hBDrhRBronGvSjRhlO8cIcRWIcQOIcTCROYxGoQQA4QQ7wkhthv/LSPqp9r7C/Y+hM6Dxvl1QohjeyKfkRBC2UqFEA3Gu1ojhPhVT+QzUoQQTwohqoUQlk4sKf7ugpUtsncnpUyKP/R4SBOBMmBOgHR7gIE9nd94lA+wAzuBMUAmsBaY0tN5D7F89wALje2FwN2p/v5CeR/AecDb6CGHTgC+6Ol8x7BspcAbPZ3XKMr4FeBYYIOf8yn57kIsW0TvLmksBJnmYTBCLN9cYIeUcpeUsgN4Hrgw/rmLCRcCTxvbTwPf6LmsxIxQ3seFwDNS53OgQAgRZtDoHiGVv2shIfXJr7UBkqTquwulbBGRNIIQBhJ4Vwix2ghzkU4MB/aZ9iuNY6nAECnlQQDj/2A/6VLp/YXyPlL1nYWa73lCiLVCiLeFEFMTk7WEkarvLlTCfncJXQ8h0WEwEk0MymcVKStp/IIDlS+M2yTt+7MglPeR1O8sAKHk+0tglJSySQhxHvAKMD7eGUsgqfruQiGid5dQQZBpHgYjBuWrBEaa9kcAB6K8Z8wIVD4hRJUQYqiU8qBhdldbpUvm92dBKO8jqd9ZAILmW0rZaNp+SwjxiBBioEyfMPap+u6CEum7S6kuo14QBmMlMF4IMVoIkQlcDrzWw3kKldeAa43tawEfiygF318o7+M14BrDY+UEoMHVdZbkBC2bEKJICD2+sxBiLnp9UZPwnMaPVH13QYn43fX0aLlpVPwidMVuB6qAd4zjw4C3jO0x6N4Qa4GN6F0xPZ73WJXP2D8P2IbuAZJK5SsEPgC2G/8HpMP7s3ofwI3Ajca2AB42zq8ngIdcsv2FULabjfe0FvgcOLGn8xxm+f4FHAQ047f3nTR6d8HKFtG7U6ErFAqFQgGkWJeRQqFQKOKHEgSFQqFQAEoQFAqFQmGgBEGhUCgUgBIEhUKhUBgoQVAoFAoFoARBoVAoFAb/H104wEZ9XcC0AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, clear_output\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "filenames = []\n",
    "\n",
    "for i in range(N):\n",
    "    print(i)\n",
    "    for j in range(numofmol):\n",
    "\n",
    "        new_y = traj[1,:i,j]\n",
    "        new_x = traj[0,:i,j]\n",
    "\n",
    "        if i == 0:\n",
    "            ax.scatter(traj[0,0,j],traj[1,0,j],color='black')\n",
    "        #elif i==N-1:\n",
    "            #ax.scatter(traj[0,i,j],traj[1,i,j])\n",
    "        if where[j] == 1:\n",
    "            ax.plot(new_x,new_y,'deepskyblue')\n",
    "        else:\n",
    "            ax.plot(new_x,new_y,'lightskyblue')\n",
    "        ax.plot(line1[:,0],line1[:,1],color='black')\n",
    "        ax.plot(line2[:,0],line2[:,1],color='black')\n",
    "        ax.plot(line3[:,0],line3[:,1],color='black')\n",
    "        ax.plot(line4[:,0],line4[:,1],color='black')\n",
    "        plt.xlim((-1.5, 1.5))   # set the xlim to left, right\n",
    "        plt.ylim((-1.5, 1.5)) \n",
    "\n",
    "\n",
    "    plt.savefig(f'test_{i}.png')\n",
    "    filenames.append(f'test_{i}.png')\n",
    "    #ax.clear()\n",
    "    #display(fig)\n",
    "    clear_output(wait=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "with imageio.get_writer(f'{numofmol}mol_{theta}deg.gif', mode='I') as writer:\n",
    "    for filename in filenames:\n",
    "        image = imageio.imread(filename)\n",
    "        writer.append_data(image)\n",
    "        os.remove(filename)  # delete the file after appending its data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now plot all points from the origin (universal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "origtraj = np.zeros((N,2,numofmol))\n",
    "\n",
    "for i in range(numofmol):\n",
    "    origtraj[:,0,i] = traj[0,:,i] - initialrot[i,0]\n",
    "    origtraj[:,1,i] = traj[1,:,i] - initialrot[i,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABGJElEQVR4nO29eZwcVbn//z7VXbNvyUwm+2Tf9xBC2CTIegMKCCqIIqhEQVBR1Nzr76vXe6/eiAqKoIiIiiiIl0VkRzCsAglM9n2ZTGayTpLZt+qu8/ujqtPVPd09PTM9S0+e9+vV6VpOVZ2aSp9Pnec5z3OU1hpBEARBMPq7AoIgCMLAQARBEARBAEQQBEEQBBcRBEEQBAEQQRAEQRBcRBAEQRAEIEWCoJR6UCl1WCm1Mc7+pUqpOqXUWvfz3VRcVxAEQUgd/hSd5/fAPcBDCcq8obW+NEXXEwRBEFJMSnoIWuvXgWOpOJcgCILQP6Sqh5AMpyul1gH7gdu11ptiFVJKLQeWA+Tm5p4yffr0PqyiIAhCevP+++/XaK2HdefYvhKED4BxWutGpdQy4ClgSqyCWuv7gfsBFi1apNesWdNHVRQEQUh/lFJ7u3tsn4wy0lrXa60b3eXnAFMpVdIX1xYEQRCSo08EQSk1Qiml3OXF7nWP9sW1BUEQhORIiclIKfUIsBQoUUpVAd8DTACt9X3AVcBNSqkA0AJcrSXNqiAIwoAiJYKgtb6mk/334AxLFQRBEAYoEqksCIIgACIIgiAIgosIgiAIggCIIAiCIAguIgiCIAgCIIIgCIIguIggCIIgCIAIgiAIguAigiAIgiAAIgiCIAiCiwiCIAiCAIggCIIgCC4iCIIgCAIggiAIgiC4iCAIgiAIgAiCIAiC4CKCIAiCIAAiCIIgCIKLCIIgCIIAiCAIgiAILiIIgiAIAiCCIAiCILiIIAiCIAiACIIgCILgIoIgCIIgACIIgiAIgosIgiAIggCkSBCUUg8qpQ4rpTbG2a+UUncrpXYqpdYrpRam4rqCIAhC6khVD+H3wMUJ9v8bMMX9LAd+laLrCoIgCCkiJYKgtX4dOJagyGXAQ9rhHaBIKTUyFdcWBEEQUkNf+RBGA/s861XuNkEQBGGA0FeCoGJs0zELKrVcKbVGKbXmyJEjvVwtQRAEIURfCUIVMNazPgbYH6ug1vp+rfUirfWiYcOG9UnlBEEQhL4ThKeB69zRRkuAOq31gT66tiAIgpAE/lScRCn1CLAUKFFKVQHfA0wArfV9wHPAMmAn0AzckIrrCoIgCKkjJYKgtb6mk/0a+HIqriUIgiD0DhKpLAiCIAAiCIIgCIKLCIIgCIIAiCAIgiAILiIIgiAIAiCCIAhCkqwst3h0u9Xf1RB6EREEQRA6pcFyMs1UNMGPy0UUBisiCIIgdEq+GU5HFsTpLQiDDxEEQRC6hYjC4EMEQRCEbrOy3KKmRoRhsCCCIAhCUtxYFl5eUhJefmAfrBdRGBSIIAiCkBTFxeaJ5dE58AVPQvvn9sGft4kopDsiCIIgdJnHK6GkxGTFgrBIVDbDneJXSGtEEARBSJphGR23eUWhHXE2pzMiCIIgJM3nZ5kxt3tFAUQU0hURBEEQusXDWyMb/RULTAp84XURhfQjJRPkCIJw8lHV4nxvPGzx8gFoszuWWVlucflYmF4Su2chDCykhyAIQtK8tNdCedZXlls8U91RDD7mGYH01D54cqf0FtIB6SEIgtCBpiaLx/fB/pauHTe3CJZNcHoDK0rCZqNtDfDztRZfnS89hYGMCIIgAKsPWRxvgwvLTr4Ga3+zxRO7oDHQvePLcuBT0+I7m0Oi0KIdgYh2QAsDBxEEQQBe2e98b6+zuGXO4G2w/lFpsfYYBHTqzhlPDEJ4RQFEFAYyIgiCAMwbAuuOO2/Jd66z+Pq89G6wapotnqqAmrbkjynww4JieO1Q7P3ZPvjq3PDfpSujiFYsMPlZuUWr51gRhYGH0jqFrwopZtGiRXrNmjX9XQ3hJOEPWy0OuDZzA/hWmjRY7xywePswtMcY5RMLHzCjCC51bf2VjRZ/3hG/7E1TIC+v498iJAinFsN5SZra/rrDYldjeP0LY52IZyF1KKXe11ov6taxIgiCEObn6y1ags6yAr42EzIzB06DpbXm7g2BE3XsjDw/XDIGJgzpeA/tQc2d6zs6DhTw2SkwIoYIePH2ELrytr+1xuKpfZHbpLeQOnoiCGIyEgQPX51rcke5hQ1o4K7NcNtMa0CJQjwxmJAHHxkNOTnJ1dUfNej8ijKYVpz8fS4YCuXHki5+guklJl/A4gGPKIgJaWAgcQiCEMXXZkWu37UZatoGxjh6pRTzhhARCxBiTyPcvQ1+utZi+7HO62uoyLN0RQwALhrX/QY8lploZbnF4cMD4+98siKCIAhRZGSYLI8ShQc2Q0X9wGis/m28ybfmO537PGBIVMI5S8MTe50GdmW5xQObLVpaYtf9y5PDy/+q7v793dXFNBWHDsUu/2A1/H33wPg7n4yIIAhCDIZmmHy0LHLbo7tg7ZGB0Vj9dK1j+28EmtodG/yKBSYfKQN/VPehpg1+vtURiB+VW7xaFb6H/HzzRCPw2uHu16cLg5kA+N3+8PKKBSYlHlHbVAc/kjxI/UJKBEEpdbFSaptSaqdSakWM/UuVUnVKqbXu57upuK4g9CYzi00WFkdue6GKiAa1v/C6gr0pp2cVm9w+34xpUgLHL/LekXDv4WfrLa6bHt7/VBdTTNxY1nmZRCxx/75fmGVyxZjIekpyvL6nx4KglPIB9wL/BswErlFKzYxR9A2t9Xz38189va4g9AUXlpmMyInc9t4ReLKPzRrXrgdjFSj3s8lymvwaj4PZ24DO8wjZ12fBxPzY520Nwu+3hte3NnStXt5Z1LYfTfw32XrEOiFEIZZ6hqtOG2ZKGu1+JhWjjBYDO7XWuwGUUo8ClwGbU3BuQeh3rp9mRgxHBdhW58QtfHZ6547VpiaL/6uEQ63gDRW4aTIU5jvH//dO+K+qyDf/RITOU+KL3B4arXNxmclat4HWGj4x2blOe7vFH3fBkVbiEmqEDeCcUXDa8OScx09Uwori+Pufqop/rUSsLLe4ZDTMKZVRSL1NKgRhNOAdVVwFnBaj3OlKqXXAfuB2rfWmWCdTSi0HlgOUlfWwPyoIKcI7HDXEgRb41SaLm9xJY9YdtnjtIDQnESNwNAjF7/tIMpygA9nEjx8KiYKpHAfzb7bDLXOcfRkZJp+f4Swfa7b4c4IcRjbwz/3w+n6LbyYYEjomO5wKOx4fxHEiJ8uz1fBytcXXZWhqr5IKQYhlroz+3/oBME5r3aiUWgY8BUyJdTKt9f3A/eAEpqWgfoKQEr42C+6Meo2pa4/9lqs1NGmoscFGMd6nMRQct+GlVh/rLNAJLLY/Hge/OQzZwLoWyMJmrB9m+G3mmprcToy9K8stPj0JHk7Q4A/NMU8Ixe46i79VxJ7TwNfJtT493ez0Tf+l/R23RZuHjh+3+Gs1xBsxG/KVSLxC75EKQagCPNnPGYPTCziB1rres/ycUuqXSqkSrXVNCq4vCL1O+SGLF2M0al5ebjU4HFQcsRVHbchVmkuybGabmnbgvVbFy20GVieuO4Xmjn1Q5tOM89ucm6cpNUDF8xTj/JBvX2Dyk3LrhNnp4V3h/X/cZuE3oCkAbUGwgk6CO1s7b2+J3rzabXh0u8XVUztviH+7yeowzWZFbccWPlaj/uuK2Oc8qxTe9IyAElHoPVIhCKuBKUqpCUA1cDXwKW8BpdQI4JDWWiulFuOYJ4+m4NqC0Cfsa4q9fUw27G6CDMMxz2wMGGRic1GWzWkZGgW82W7wbrvBKabNldk2zVrTrDVtWtGmnS72CJ/Gp6DYgJE+Tbbb+Le6jXYiMQDH95DoLb26uTt37fxQDQVnDk+u/JH2jtse3RO//HsHLF49GHuft9E/azQdMqZ+e74f1dkfRugSPRYErXVAKXUL8CJOLqwHtdablFJfcvffB1wF3KSUCgAtwNV6ICdREk4Kai345Hp4egFkdmIWOX8EbK5zllcsMHmp0uKDo47t/O5GHxdk2SzLtjk/yyaoIduAXZbi1VbFaZmar+UF8AFBwK9C//UT/wS0hixFbKNsNyjNglwTCjKc5VG5UJIBphk/g+lNk51Yhc64sgwer+y4vTZG7+BDpc5kOS0xbn9sDlybYG6FX663qHcdL21ByJLkOylFktsJJy1qlfM9Kws2Lum8/MpyiyYb6gwYBRRnOEFfe4PwbruPT+aEXcShxr8z2jXsDShqbMUxGwoNzSy/Zoiv82N7igIWDYXz3BQU8XoYyZpnQsffWBYejprssNHl42Do0OSuo7WmKQB5pvQOYiHJ7QShi1z0fnjZTuKd6J4q2BtwooCHa0DBL2oNJvo0SzI14/1BDgdhmGvr70wM1rTBE61+tPv670ezLDPAOdnJp7H+1nz/iXxE3obXT3LDVzWw+hisjuHFLcuBStfM1NBgJdVLCPEbd/hprN6BlywFX+vGlJpKKTpJxCp0ExEE4aTkJU8A1ldHxS/39BH4zFaoD8LlWYpTMvSJ1BBzTc1sM6wmGyyD19t8fL+wY0NYmgGfm2WStwq87ogSQ/Px7ABjfY6QhMTANKAsD84brRia5efAAYs/uLb2WG/sirAB6nbP/vp6i1/u6lC8Uyo9Pod7dwI49zQuB66JY9LJMaDZI2b3xfEdfHgELB4pLfpARARBOOnIXhW5/sVxHcv8ZKPF0w0G77QpLAxKTZjt1xF5gqb6NS+1Gqxvh3OzNOdl2dTYip81+rg1N0iLhv85xWn4slbB51d1vM4Ynz4hBj4F10xWjMnr+LMcOdKEg06j/PgOiyunRDaoITGIdoUUFJisWAA/Lre6HfPgZW9zbDPQ8Ez47Hj41W5nvafmJ6F/kOR2wklFbS0kCNKlsgXO+FcQqx2WZdl8ryDIfxVY/GeRRa7Hrh/QcEeDn3+2+ag5x8cPJ2h2BxRXZgcpUvBeu0GWgsLXbNSq+Mnf/jFdcfFYAwUENfx5p2Z7beKme0dj5Lq38Y03y9v1UztuK0h4la5xqC0sBvEQMRj4iCAIJxVD1sbe3hiA09+Hce/CJJ8TN9DohwITTOUEoHlZ3a4oUpr1053Ge0qpyaPNBnU2XJ8bRGkbAzgzM76D4uhcKC5WLBjm47a5BsOyHH/GE3ts/rDNIhiMFIYLRnQ8x6Pbw2JwaYLA/mG5HRvjUHBQaZbTWM8rin88OCOAhMGNjDISBiSvVAXY2+j839TuPycGa7oL2rPP8xX57dlXa8F2S2EAczI1m9sUCijxa44EFAFgvF87Qz2BLF/42IDtvMHHYlQOXDfNRK2CoYbmlrwAjTbU2orxfs3/NcCyUpOVnsnrj86FoUM7nmvzsSDP7LWxcfwCy8YZzBka7pqEegNnlMCHxkZGCHf2Bv7wNouqGPEIsY6rr7d4YLczCipEkQlfmt2x7K/WW9R1Yo86fxQsSjInktAzZJSRMOjYWa853tUk+1EoIgO6mmynwc4ADgYg19BooNGGPENTaHBCDHwqcvSRTzlDSUMNZFUAJpjOaCLTCA9hPWYr/tDk48bcICWGRim4Og/+IwkxAJg51Me0QvjTTpv9zfDsXpt3D9lcN9UgwxcWhrdr4O2ars1pPDWPmIKgte4Q4FVQYPL1+eH19qDGjGNPuGlux2sfPGjx+wPhdRGD9EAEQRiQfHGmidaaloBjzmm2nLHnTQFNk9XxuznQMczrrTaDv7cmN6C/+Wz4e4XF7nrI9cOtczo2YHeuDTfA9zWZBJa6QlAbWW5v0CDHFzyRF6hNOcNKA6iEYhDC5/Nx3TQfO2qDPFlhU9MKd623OX80zC2E9XWR5ecUJT7fq3st3osx9/FnSmH06OQa6gxf18b8jxhhsiKGiUsY2IggCAMWpRQ5JuSYcCKXQwyer4HLNmoyFOQpyFOafAOOJjme/46xYOoAu12j+memxL5WZQBG+JwUFUHCvYIQ87Ngreuxvm2eyY/KLTROr+NDGTZ/ne7rVAy8TCny8Y058Nhum72N8HK1TVFGx3KXTIjdqD+902JzjPkNcoEbZ0BWlry1C5GIIAhpxfajMH+Dk/8kEoXlZhg9RDgyrEhBbVTXYVQGLMiCZ10B+OYk+OM2p1BhBhTFyYfQZCvwaTcFthudBpyeC2+f6tZilfN94ftwaQk0HXXMVhdk27zebHN5ktG4IXw+H9dM8fHTHUEO1tlOyk8PsUxFf9lqsSdGOuqu9AiEkxMRBGHAsv0gnL4djiX5pp8FrJsDU91JWnY1w/nroLYtvP/x2bCsBPJXhY9rtwInkr/dMC1270CtCvIf+Y5oZCoowmZxvo8XT4ldl5cb4OUGk+qFFg+5gWEv7IfLx8QuH4u7K+Ebu0NRxz7Ax1VZFqdkhsvUt9sUZDjG/d9utIg15bMIgZAsIgjCgCSoYdrW+Pt9wN9nO/6FT7jlyjIcMTjWDuetg7VRGUoDwCUbO57ro2tgaZYz/DLL3/EnoVZBiWGQb9gcCsJwH/x3sc0tczr6Jz4zBP543Fk+swBGFZjMzHdMN+NN+NMWi2tnxG+cI0Wg4z1Hdwh+uSmIn2DM8jLuX+gqIgjCgCSWD1MvDS+fsRqWRTXu29s72vW9xGo0DTRnuLEC18XwHYTON8HnlBnuakBjAKobLEZH5fh5rDa8/OZC5/ujk01eWG1R5od9rbD5qMVMz1zEiUQgxKQsuCnTIlb8b/RxIgRCdxFBEAYseqkzsXzIBaBWQSbxo367wrn58Oop8KtNAeraYUwu+KN6B15xmey30RrOGQ6vu5O1/NGT4wfg7Tbw46fN9S2sOg5Lhzj7ziqFSnekz9/2wohsi5HvmdEuAQDGZ0KF5yZvGwXfKbb4TYz00tEMzYSaFpuSbIk5FbqOCIIwoLGXwrS3nbd/6JkYeHsYAMdaAycikD8+QdEQgMYgNAVhynsRRzLOr9keUPzHDj+nmQEm+CFHafIMKFROYrczMuGMTCcorTqo+Psuzbs+KMmCK8bD1w8qpmc4sQn3b4N2bKKTBfxoPHy7IrxevhDmF8DK8vj39a35fhoszRO7gxxqgQe2BplaZPPRcQZ+Q4RBSB6JVBbSgj9Vw6d3xN8/RMHxTv4rPz0b7qoKN/pjgkHOzbJ5sMnHjkD8hnOI0nyrIMDfWgzeaY8d15CrNCMMzUifZrRPM8anKY6a9lLrjutFJSY3x8gBlKmg9RxnOZk5BUJmoh21Ns9WBmkNOma3c0YaLB7eB5MrCAMGiVQWBj3XjobnjsOfo2bhLsuEF+bAjDxnvaYGhsVwHAPYODEEQ/0wymezwLY5HITLRhrk+SDXgG96GudxJvx6PDxW7VjpDwfix0I0acWuoKIiCOULYE6h42N4phKOu70QpZxI5/1BRaFyJsF5+4AF+MhRzhDZYkNzVqZNmT9xryAeU4oMWvc4Y25tDa/ut1l9xObyCT5G50pvQUiMCIKQFrQG4J0YQVaVbTAzyU7kZSXOB+DOtUHaFWgf/HiSsy3aIb3Xgot3wJXZiiZbs8cd/qqwyUNTojQBpUApMtBkKUWW0nx9Mww3NKN9dJj5LEM5+ZJCzMqA/83oemLqDODrC8wOvYdXKsPrt87x8dQem8pGzR+3BxmfZ3P5BIMsvwiDEBsRBGHA8/098J97e36eUIOfo2y+le9EMr/SanDaW0FQcHaGY6rJUpCpNHkK8pVmtPsr+UFBEKWiG++em1zfbDPwoVmcoWOOroom1iiiH5VbfHuByeqj4W05foNPTTGobLB5em+QikbN3RuCnDZcc/YIhSH+BSEKEQRhwKI15L4BLe6b+eQsaLegMskXagOIFdM23e80/KN88JncyBK2doKBWzW0aadciACQ7IDOoIY6Gw7aioNBOG4b1NqaZq04ZmtaCXcdzjAtPhKVWjqDDkHJXDUOJseJdNbAw1tj+xrK8g1umW3wzsEgbxy0+dchm7U1cOk4mFQooiCEEUEQBiyvHg+LAcDORDPbxCBegPN6S9Hc5KMQTZNWHNKKNq0p0EEsP9QGTdrddNvnZAS4ONvpBYTmdG/RcNxNb+18oFa737aiUWt+N8XgO6Od8rdshnsPd6zHeJ/N1TlBvG1yHnDLAhNba+5Y6/guhpqwPEbaaXB6CyGzUZUnXcW0/I5ll4zwsbBU8fcKmx11mr/uDjI6J8gnJvnIFDOSgAiCMED52Fp4srZ3zh1A0RxQeAOhZ2cqNpzesVF8pVKfMMOYCmYXwaerTGoCcPEQ+N0kuKMaXqmFhnZo0ACK63fAVUWQmwv3zHQ+6xtg3vtBpvgV52baTPDrEym2i/1woyfDqqEU357vxwZ8nqFJ+/dbPORJpR2PKyZ3FJCArdnbABmGk7LbsmF/c/x5HoSTDxEEYcBwwwb4/dHY+xTwp7Fw+0HY3/kozE7xxnjNzoQNp8cud16ZyXmemcguWAs1ASg14fl5zrY7p0XVdZXznbc6HPtw/fsWa1sMbs7VjPXbJ3o+U3Lh43EmrVdKcShJAYhH0NZUNGi21Dq9gragM/HPjCGKGUUGpdmQE2+iA+GkQwRB6Fe+vgXuStDgbZsDczc4eXyW5kD1mc72RCkquko8MYjmp5Xwj1qnLttPjV/u0KkwfLWzXLDK4qJsGOODq3NtjtnwbIvBx/Jtbosx5wLAC7st1tbF3BXBdcNh1KiOI40A9tTbbDlus71O0xqETB9MKXREYHy+wmd0bX4D4eRABEHoc364E75TFX//Nk/G0vKKcHTyqG3AtsTnjuWM7YxlH8BzCxOXWd8At7sxCi/NddJkx6M0F5YXOD2CszJhmA+OBeFvLQbvtSvKMg1ui5EY72flFoncJCEB6IxMA/6yK0iG4YjA9CEGE/IVfhEBoRNEEISU4FsV34mbDG/PgdOLI7dVH4CFFV07j1cM/jQerk3i+OfrYc6/4vcU2oNw6gfO8jfHwIcTTHLz3C6L9fUwDpiQY3MgCH9q8rE5YGO7I4v2RF0nXiRyZ2mrvY5nLxMLFDOGGEwsEBEQuoYIgtBjtO6ZGACcswHal0ZuuysJs0kikhGDEBvbnER69tKO+2asdiKM5+fCHZNjH//IVou9nlE+B9wYhy0BhTMA1rHT713s7K+utvhjjJFHkDhbqdaaqibN1lrNljgTRVw2QX7WQvdIyf8cpdTFwM9xzKsPaK1XRu1X7v5lQDNwvdb6g1RcW+h/lILGUx0naggfEFjqLFc2wwcNsLUZ9rZBdRv8/XjkOSzCfoEiA45/CH4yHZ48CrtT4EQOEXIgD38NDkeNrtFE+ibygY8Nh92tTlqLNTHMSr9eb3HcExdxzIYnmg12BRWXDjXYPDd8zlnZsPmgxZ9jCF0mcFscIdBas79Zs+W4ZlutTYMFfgWBGKODhiUwZQlCZ/Q4uZ1SygdsBy4AqoDVwDVa682eMsuAW3EE4TTg51rr0zo7tyS3G9i8d9DiQAtc5pnTd9o7sN1jCH9qBlw2PPbxyTiGJ2fAjjO670Qu8UFN1zNDxOWaIfDneXBPuUWjZ3tAwyPNPja7SfKenAmXl7r7bPjmuxbDszqeb34hXDyxoxBorTnYrNlSq9l63KbecpLVTShwHMMEg/zd9cN8rAyecIdNpXIuBK01SonJKd3o7+R2i4GdWuvdbmUeBS4DNnvKXAY8pB31eUcpVaSUGqm1PpCC6wv9xKvu09ux1uL2+U5DtG0JvHAY/s19+pdvgQm7YPcZHY9vPg1y3g2vHz0Npq+BI54GfGcnk950RirFAOCR4zZTPwiSAYTM82ta4fE25/4zFDScDhnum3rIPxAtBrH8A1prDrXA1uM2W2pt6tqda0zIV5w90mBKkSLLzW2xsjx8Y+/EMT11l5+WhyfiKTDg5nky4c7JQioEYTSwz7NehdML6KzMaKCDICillgPLAcrKyqJ3CwOIEdlwsMV5O76j3OK22WCaJheXgi4Fc5WT7mGP26g3nuoEaoXIzo4UheJ3nfXsbGfdOznOQGGYocnEMZNZGsrbFS+1O72Cubmw7tSu+Qe01hxphS3HbbbW2hxvc7wN4/MVZ44wmFqoyPJHvqU/uydsQ7uxjKQmzkmGWM5te6A9AKFXSYUgxOpTRv83SqaMs1Hr+4H7wTEZ9axqQm9y/fTwGHgb+OlGuG6qxahcp9GzlsLl6+Bvrr8gbzV8cxTcMTV8jmhRyHkXDkyHkQnmU+4/NB/LslEKXmtVlPrg1AzNwowg7TpIjhE7ZXUW8LUoITjSotla68QKHGtzfiDj8hVLSg2mFimy/fFNNRtqw8vFxSZ4MpyuLLf4/BgYNqxrb/VvV8d21NwyX3oHJxOpEIQqYKxnfQywvxtlhDTk+inwe8/ENQ9th3NGWpw+wmlInpoHh5vCgVo/3g+/OggNHwofEy0KscQgE/hkHmxshQ8STT7ssmkyzNrZvXsKMTHLSSkxOhMmZcNbh23Gu+3jhnaDqnMMVpZb+LQzY1o00f6Bo636RE+gxvWzlOUpTi1VTC00yDU7t9f/dlO44Q71NlZEpcH+bRVQZXXJn/B6VI8mlF5bOLlIhSCsBqYopSYA1cDVwKeiyjwN3OL6F04D6sR/MDgYkWdSYFrUe14wXzsAVY0WH3fz6ZTmOikchr4Bx4PQaDsmpN2nwgTXhBQtCiH+tQCWFMa+9u/3wQ27Yu+btROuz4bfnQaHDsGILYnvwwe8Og/OWeesr1kIpxSE949fZfGJfCcLaouGaq1PNMLx/K5r6+BioLrJ5oXKIEdcERiTq7hgjGJakUFeEiLg5Uhosp2o7SsWmGw7YvGkJ+BvZbnFgiFw0fhww75/v8Ujh7wzQXekxIQvxEmmJwxuUjKFpjuK6Gc4v6sHtdY/UEp9CUBrfZ877PQenN9HM3CD1rrT4UMyyih9eGKXxfb6yG35Jnw5qmH54W74jsfmfX4BvOwZztnSEikKXp+Cl611MKfc8VGkmh9OgH8fF7ntU+8GmJuhCWjYHVAMNTQlUcHG0W/qXkbnKqYXKaYXGeRndG/kjtfZm+jt/551Fo3dDAyZlBM/t5KQHvRklJHMqSykjCNNFr/dHrnNgBPO5hBag/FaZJng0vB6IlGoaIGp73Z8w52WBVuXOMuJRiXlAk1J3IsiHKS2stxCu+mwQyOL2rUzoiiWfyDRHMg9GRYaOm+2gq+6tv39+y0ePtS9wEA/8KnhnEieZwLfEDNR2iOCIAwo7t1o0RDVJl4/DUbkRDY2s96FzZ7o3kemwdUjneVoUdh/Gkx4N5zXyEu2guZz4MABN99Rivl/+VaEj+CZFoPydk3T0o4W1x+VWydGS3xznp+71gVi9mK6KgyJRKa7XD4ajllh/0EqYxiE/kMEQRhwrD5o8UqUl+i8kXDqiMhG57XDsNQTsTLehD1uRtNjLc5Q1GhyDGj2vBKHUkz3JF5hVha8dgoUmzDyDTjoiV/4VE6AXKUpNhxn6381mBHXDXHvBosGT+vvbWDvKrdiill0I5zsfAexSBTt/ONyi85CMkQQBgciCMKApK3N4uebI80ZUwvhYzEiczNWRZqBQjELXlHIVHDwNBjyTrict1FOZUpsE/hsCfyxBtrQKDT/ryDIBkvxZIvTMzg9F95202A/s8tiY5QPJbqBTRSf0BVimamisbXGsiEzapLmjYctnqnuWP6yUTBjuAjCYEAEQRjQPLrDosKT5yHHB1+Z27HxuXoD/MUzQc5XRsLPpznBUQ1BKPRHNvoNiyAvr/PrL/sAXqzvWQK+CT6b5XlBHm32sc5y7EdDDWd+gwxsvleg8SYWvXAkbDgKB7qai9tDLvApT+DZGcPgQ2MSN9p3rrVo9/ykb54EBQWJjwlqHTErm5DeiCAIA56qeouHo4aIfiPK2QyRMQvgmIea3JiFrNecie8BXp4B58fJkZSI6F7E50rglHz48p7Ex52dEWRZts2+ZvjuZJMLd1hckwuvtSmyFSzO0Gxth+ndSC6XKD2E13cQz6TzfzssdjbG3JXwOGFwIoIgpA0/W2fR6nlVv2EKDM/r2GCVvAFHPUbvcSbsddvG5cPh1zO6d/1YZqUMnMRxLe5P4bpieChqKs9v5FnkG/Cf9aG62tyUa1Pmdw56q01xZmbi31KhATd1IS/QxqMWz7i9gyvLYEpx+NiddRb/tzv2ccMz4cox8EuPAIsonDz0d3I7QUiar80zeaPa4i3Xlv67HTClwOLKSZENVs3Z8NMKuL3CWQ+JwfSs7ouBl3ygwV1uh4hEKtFiYKDJN+CDdiNia30wCH5nGG08MZiSB1dO6V5j/IwnXiMkBvGcwz7gm1GN/lnDLN484iz/dK3FNyQNhdAJMru20OecPdrktpnh9R31TnK8aL4xPtJpnK1gy5LU1KGBjqOE4jHFr8lUkIXN5AznuB8XWszOdPYnMr/vaIQDB7o+ZPTxneFjRmQ6pqOVMcTgqonO23+0GACcNcak1DVhWRr+vC31Q1eFwYUIgtAvZGaarFhgEgratXEavSNNHRstvRTqz3JiDXrCkSOR652NSipwr319rtMML8l05mb4zSarSxHSfzjYhcIuOxrCywejxqtOznNEYMUCk8mFid/6PzfLPPEjr2yG8kMiCkJ8RBCEfuXr800WlYTXf7sd/rqzY6OV3w3j5mXlTqMf+pRu6trxUaNIuW2ByapKi6PuyCFTxbbNl8X4VSUTWFZXZ53oCUQTutaKBSZXddEE9S1PHV/cD/X1IgpCbEQQhH7n/LEmn5sWXt/V4JiQ2tq61nA1BCIF4OlO5mTWS+Hs3MhtIw2o8+RWGrMqXIeqeot3PP6FeDb5So/TfKxnYpx4ovCrjY4I/CqGk/hTrkmop/Z/r3D9Mk5CQEEQQRD6ndrWAOuOKoozw9ts4K7NsPZI8qKw/kjs7QuyncY/9AmhVsEbnsRG706B/R+CggI4y41vyPWFM9g9HGfUTqIRPBlRCfD+4trxX9sX7g3UxbnFFQtMyjoxCXWFmyeFl3sjFYaQ/sgoI6HPqG0NsOYIVDRo6iywkogUe6EK/nXY4qZZnTeMZ44EPLmMLiyAFxfGLR5BtIP5jUWOYIz1a7SOdBx3ZQjnrqbILKh7mmM3xrk+uHVuuFy+r0ORHlNQYHLRKIsX3ZlIflJucbsMRxU8SByCkHLqWwOsrnEb/nZo76ThzzBgSCZMzFecMhzy/H601vxua4DDreFyX5gJJZmdN2BeZ/GlhfD3BeH1r26Buz25grKAlqXxz/XVNRYjPY3zpyfBmKjI3wc2WdR4IpLzfNDoGQ4UnXvJy00TodDtBSQThJYKHtpisd/9uw7PhBtmiigMJiQOQegXGlsDrDkKu+s1tUk2/EUZMKFAcWox5GXF/++nlOJzM0y2H7d4osLZ9sBmmFVk8ZEJiRswvTQsCs/UwafXw8Nzw3M8e4kWg8Ymi3V1cOYoE601JR6jqtYhs5HFZ0ph9GinHjVR6Skag05De8gdHRRPDLyN/s/XhsVgYhLpOHrCdTPME3MrHGqDd/ZbLBkloiBID0HoJneus+IKQIYBBRkwIV+xqAQKEzT8yfLTtRaW+1/Vr+D2JJysySS7az7N4uE9cDRqaOe/jYWNNbCvBYIanm81mGfajE1wK1+eDPe603YW+KDe00soNuHG2SZ/2WaxpzlxnfoqqtjbI0km55GQHkjqCqFPeXx3gB11GgUMzYTx+YpFw2BIChr+RDy312L9MWc52UYzWhQUFr8ZDr8/bjDdr5lkxv7/f94IeMWNH/AD36xzrleiLK7LhRKjY0BavnIytrbGOKW3vvEcutkKbpnn79NEc31lphL6DjEZCX3G5mNBdtQ5Ld6XZxnkRQ+j6UWWjTNZf8xpwPbWWowrStyAvVABYFNiKCb6bWb6NVP9cKQVLsm2qff0cCbkwh53xFG+LywGALcvMCneB5/bBTXa5E43kVyxsvlIts00V1TqbWd0lC9Oe/7MbouNMYbCjgKu66fG+OZJ4WGoK8stEYWTHBEEIWka24M8vddpRT86rm/FIITCSTu0rQ7GFXXc39Ji8XglVDU7Nv//yId81w9QZ8P2gGK8T/NAk48qG96fZrBwJByus9jjxgEEPW/4oQbyhrGwqxV+VB32QxzVBtsDmmmm5pFmH2N8mnmmTYFyrt2onWs32PDd9yNnXfOyH3h4q8Wnp/d9Y1xQYHL5OIun9jrrd5RbEYFswsmFCIKQNA9sdcRgWqFi5tC+FwNw518GDrhTb2qt+cuOABUxJkpWyjHjbLMU1QGoCCgqbAOL8Cv8km3QPhJKC01CU/R8ZV54+Offd1t8xJ3Q53+mOB8vK8tttIaKQJD1lsnzrTDPD+dl2RS7ApBvOAKx04LdAcWrZzg/O2+iuqoW+o3pQ00m1ljsbnJ6OPdvtFg+W0ThZER8CEJSPLLDYm8j5Prh1jn911jctc6irScz3eDMqdCkod5WFBmaok7CM+OZURoarBNO5OgyD2yyOOI6qr0ugaNBKE6gpf1psvH+bZcUw9IyEYV0RJzKQq/y/pEgL1fZKOArcw2yff3TOwD49UaL4zF8siWZcO0EyM42IxylJrDbgn+2GxQqWJZlU2lDg63IM2CEoclLQbz+qCy4aCgM90xDGT2P8uttigk+nXCkUmmGk5Cuv5CRR+mPCILQa9S2Bblvs/PaeNUEg8lF/ScGAI3NFr/c5sQzfGaiIwCxCDVsoQbWO9rIG6cA8MZsyApavLq/Z9NsdsZ/1JmA5rsFAbISDCT67AgYOXJgiMKKBSYPb7VOmLSmFsDHJolIDGREEIReQWvNnesDWDbMHapYNi59XE7eRu22eSZZr0fuPzoXiteH17fMhukl4eNCppvQeo7h+BYA1h9yxKOVrtGcb7KmJsCHOplZzXv9eGw7YvFklbP8jXl+TCN1Q1Xr662IBHiZBh3MdKWZ8DmJcB6QiCAIvcKL+wKU12gKTLg5DZ2MocY8A/hFs0lFJ/nctsyGp/Y5hc4ZDqePMvndZutExLG3kd7fYPFQHP9BxDkPObmDFg+H192UGe+3K07JcH530wvg8kkm71Zb/PNw5LGh8z67x2JDbeK6p9r3sLve4rEksqJmG/DVLkwLKvQ+PREEyXYqxKSywaa8RpNhwPKZ/Wsm6i6h/9ztwJ4z4bnZicvP2AjaHfYTary9eX5+7Ol1/GlncnWYMdzkawvME+c7GIRnWnyE3sMud80vp402OzTqoWyonYlBb4SxTSwwmVnQcfuKBSZzhoTXW2ynnj8ut2hulgyq6U762ACEPqPdtvnrbqdl/PhEH34jPd8bvuXJMhoKuor2H0TzvUYf3y8IRrSyc4tgfS0R01eGlsdmdzxHwevOhDZrFjoBb97pQZ9qgrkZ9omRR11JQ+0Drh8Dw4ZFmrMykj5D1/joJJO96yyaPOaihgaLS8abXDIeVh+0eOWAsz0I3L0NwOIL06AkR3oN6UiPfulKqaFKqZeVUjvc7yFxylUopTYopdYqpcQGNMB5bGcQy4ZFwxRj89NTDELESiNdmuCWLAx2BRXtHkvqMk8yvegG/NqoYLKmoBOIdiwIE1dD8Wv2iVQWQRuW58MV2V1zXX9meHje5JAYeCnOinFQirh1nsm4nPD6vTvhDdesduoIp1fzqUmRxzywzfk7bayRHkO60SMfglLqDuCY1nqlUmoFMERr/e0Y5SqARVrrmq6cX3wIfc/qQ0Fe2W8zJAO+2I/DH1NJrHw9sXoJoSjoczKDXJxls7LeR502ODoXKgMWL1Q75UoywhlOo808v9lo8WQ9NGAw3q+ZFSdXkq2dj9+AxcXw4Rhj/kP1zsSZvjPe/vNHwqIRvfusooUwy4CvRfkOmpstfrHN+Rt6WVwCHx47OP4vpQP9mcvoMmCpu/wHYBXQQRCE9OB4q82r+20MBddOSU+/QSym5cM2d9L6w4ctSktNqqfC6O2R5UIN2YGAY88p88MGKzQayeT7+RYZRmS661gmnw9lg3cAa0DDO+1Q3m5wZqZmtqn5Yb2fU/IUby3ocPgJViwweWKHxcc6mUO5t8UgFq12x9xHOTkm317gCMN9O8Lp0N+rgfdqLG6f58efwtFQQurpqT1guNb6AID7XRqnnAZeUkq9r5RanuiESqnlSqk1Sqk1R47EmRNRSDm2bfOnHUE0sKzMR15GepuKvFwxOdxoPei+5Y8aBfmeMt673RHUtGuY7HdaNB82Iw3NC63J/01CHe+nWgy+X2/yzBKT+8bYLDCdcwdRPDq58/N0JgapYuMxm621NtEWg7agjti2YoEZ4TdZWW6x+kCkKObkmHx9nmNO8sa1xUv6JwwcOu0hKKX+AYyIses7XbjOmVrr/UqpUuBlpdRWrfXrsQpqre8H7gfHZNSFawg94NlKm8YATCpQzB46eMQgxIWj4SVXDLYetphealK/NGw6ajobProW3m6wKfUpMhQsztCM8VmMMCDWi21tEIp8MKcILonyMyjl9AzebXd6WmVvwc05ThqLPAXW0l68WaC21WLjMaizoMGClgC0BcEKOvUKuiYrTUcTTyKamiyunW5GDLt95SC8c8Ti1rkdxcv7P0n1YVpvoXt0Kgha6/Pj7VNKHVJKjdRaH1BKjQQOxyqntd7vfh9WSj0JLAZiCoLQ9+yotdl0XJPlgysnDD4xAFhYavJStfMm+1Q1fF5ZPFQF/1vo7P/ZeviwAR8ujDwuC3i9zQA0ZxmawwaMcn81dzb6sDDIr4f6Cc42rwnpMxNhezX88Qjss5yG2O8ZXaSAc0fA4l6ISr5vS8pPCUB5LZyVC6PyTVYsCN9vUzB2+uxa989xZjzbgTCg6KkP4Wngs8BK9/tv0QWUUrmAobVucJcvBP6rh9cVUkRrwOZvFc4gymsm+zDSdIhpMnxudNhk9NuqxGXfalGcma1Z1Waw2nLe8l8kyP8Whn0Dl2bZPNlq0KCdnsaFGRbnuuaUM4bBuCKTh4pg1Vuwz9KstRSnmPrEkFMNvHoQXj3Y0Q8xLBOuGAFDh3ZPLEZkwUE3lDrHgAM2rG2DSktxRCvqbEU7MNUPN+QGWdVm8MKSSL/RbzZbHG1zEho2uTm/t9bCWaPDZVYsMHlws8VhN3hvZbnFstEwt9Rk67HwfZ09WpzK6UBPBWEl8JhS6vNAJfBxAKXUKOABrfUyYDjwpNtd9AN/1lq/0MPrCini6b02AQ3ZPthyXNMStBmTqwal86+01ITqjo3vjAK4LCo/T1OTxc+3wUxTs9o9ZHSUEXxxpqbWsvhn0KRQ2ZyS6Qw5faRZcZWb5uP4cZibDfssxeMtfh5vgRyCTPBDqU8zxIBCBcWGZqhnFrYjbXD/XmBvcjOa/bYKflINFa2hGdtMZvuCfCLXpsKC3zT5adSR9TfRXJYd4HAQduiOLwKhaUUvLYPHdjsCVhdjJOnnZppUNFg86pqQnquGNw9Z1LsikjV4xicMenokCFrro8B5MbbvB5a5y7uBeT25jtB7nDfKQGHTFoT3Dtu8c9gxa4zJU4zPV4zPNyjNBmOQ2H9XLDDRWndqz87NNTGUxTS/RmGjMZhvagIa7qo3+I8iG0vDhXlwW6ETTWwD9zb6OGgbTFkf/9zN+NgUgE2ByO0Kmx/mBJ0UrTH4+hqL+5tNmm2ncb4tD5p1/EC7jUEfdpPimtwgX8sL8KdmRYbp56aR8NVx8MQem+21cH+jj91nx/97TCg0yTCc1NiBOA6H8VEmpHrPvX1xSuxjhIGHRCqf5BRnG3x8kvN22B7U7GvUVDRoKhpsVu3XgE2WD8blKybkG4zPVxRlprc4dMW5qRT8fGiQ247ZzDU12wOKY/i4vdbH/xY6jd+mOsfpbAD3TDG4alvn5zWA4FLPdVaBxuDfmw0MIrOuKmw+mW0zL0MzzWfzgW242yOjp01gZAZcNQz+ZwxkZztXenJ3kG11sDxPc8nYAHNK/Bxp0Ww5brPWUtTE6B28VOncW6hTVJIJ1S2dO6BXLDD5zSaLo56huZWtMC1GRLcw8BBBEE6Q4VNMKlRMKgTw0Whp9rriUNGg2VbrND+FGZzoPYzLV+T401sg4jEmx5mKsykIPyzU2MCW1vC9/rrB4Iv54aY7ZNLRI+Of89EDcM02p8GP92YfHcesMXiixabMD1dlBznDDPL9OSb3b4VhOOm8E3HFRJOn9lhsrYVn92lsAqyrcSYKeq7VR/VpHY9Zf8z5nlnkfC8YCtXVia8T4sZZkXNSPFkBQ6otvpiGCRJPNgavB1HoMXmmYtZQg0vG+bl5lp8bZ/i5YIxBabZi63HN3yqC3L0hwO+2WvyzOsieehvLHjwjhT89LdyAhRrpJZ60ExW2wXttjkAkm2306pFwSVHn5VaOdRr60KdtqZ8zipwey2gT7t/atbQQl08wmekmlnl+n2Z/i+a5Vh+tWlEUlfqirc06YRq6ZLxzX7NLk2/MK2s71u241bW8TUL/ID0EISmUUhRnQXGWj1OGga01B5s1exocE9PqIzbvHnZMDGNyQ/4HxfAcNWj8DwC/WGjyn0ehZAOA4slWP0+2alZ04RzPzIdD7TA8Rla6UK9hxT74dlSOoMsmmRzzpOPuKh8db2LbFlvrnPXdAXhpTmSZFyot1h7t3vlDPLLH+TYVfGO+yX0brRPDT1eWW1w92fE5CAMPEQShWxhKMSpXMSoXzhzh+B+qmsL+h9cOaF47wAn/Q8jEVJSR/gFKxcXOW7tvVajnoMhdBU1Lkz9HLDGAyNncjFVgR53zhpkmd621aOtmR8zvC2VsglvyghQEbMBPW5vFL7aEncY+4Cszu3eNUNVumOZ8f2m2yfrDFs+5JqdHd8KwTIvPywQ7Aw4RBCElZPgUEwsUEwsAfDR18D84DuqCDDh9uMGCkvQYi7hwKHzg2tOvmhC5L7gUzl4DbzbCd8ak7pq3jYC7DjoN6yfWwWNRY/Rum2/yo3KrSxHGAHsbbDYe0/yz1WCjpbglL8gr+zV7Gy121ofLTSt0/A5dZWetxdN7w+tDPdObzi01mVsaNhsdaYsdyCb0LzJjmtDraK053sYJcZhaZKRNegytNT9a64yh7MvGy78qPIIontO4psZiaLE/KZNcwNY8uDVAjQXfO+4ngOLFSQE+qAn//hVw00woyEw8T3UyTCmAK+PMvXzPeotGz/Co6yY7kc9CapApNAWhF9FaY9lOL6gv8Y5C6mwkUWe8eSDImwdtftfkY3vAYGmOxUVRbfA5IxSnjwwbDWpaYd778PvJFofrYF9T165pKjh9GJwRI0r5g4MWLx0Ir4/KgutmiCikAhEEQRiErDoO565zlstM2Htm985zrFXz260BhmQrbjzgZ5LP5gt5wZhlb53uJKvbXgflbZBtKMp8mowUaqEBjMuFS0bBPTsi94kJqeeIIAjCIGXEm3DIjfotnw/zi7p2vNaaR3cGOdisuafe4NJsm5FJum8CGg7bsC+gGO/XzM2DK8ZDUVbiRntnrcXz+8L5j7qKmJB6hgiCIAxiemI6er4ywLqjXfuNHw7C/7X4OBiEP88w+PgWyFLQck7Xrh1NY6PF4/ucpHvJRDwL3aM/Z0wTBKGXOT4fhqx1ljNXQdvS5I5rtOykxCDHB+eMhHnDTLY1wnT3HezXU+Cq4TBuN+xtg+/sgh9MSnyuROTlmXx2Rux9L++1WHvMcaSPljQX/Yb0EAQhDbjofXjJnQb0+6Phu0kkjAuNCjKA7RZUBRXnZGqWlMK5Yzq+gQcCkP0mBIDzi+Dl+c727U0wbbUTm9B2NvjSY8TwSUtPegjpMfZPEE5yXjwlvPy9JHIK/etg2Gn8J8vkt80mL7X5+feFZkwxAKfRDwAFRlgMAKbmwsI85+392q3dqr6QJoggCEKa4PUfxEuMB2AFbV474MRQn13mY32zs/2eBOaem7fBbjclxsgYUdSvznW+HzsCLbEHKAmDABEEQUgj7hkXXs5bFbvMvZucFnt4Nvy0OvwTv7A4dvnXj8OvPDEB21s7linMgEuHOs7gC9Z1rc5C+iCCIAhpxJcnQMjn2oTTU9hQF96/5bhNq/sG/9mpPp48Ft537jrY3RJ5vtYAnBPVwMfzKj41y2kw3qqHypY4hYS0RgRBENKM5qVwbn54fW45jHvLWX7L9R1cNs7HowfDP++1i6A5CB9eC3s9PYDSt8PL5xQ66SsADsbIqOrzwc3ufMrnJ5gRTkhfRBAEIQ159RSoXRBer7Sc3sJHhvu4YrzBjKEGn97u7JuUCfPy4OV5UOeKQlUrXLzOmQMaYFoWrFoARe4IohvjzPr2iymQoWBHC7xXF7uMkL6IIAhCmlJY6DiayzyDhkasN7h5j48Kz3SXO093vhfmw4tzocaCKe/Bi8ed7QU+2LrEWb5ymPO9qjb+dX800fn+yMYU3YgwYBBBEIQ0Z++ZUOHpLbxaDxPedZZzo37hiwvgsZnQ6vYMFFB3dnj/vW5j3xg9j6eHr411ROSwBX891OPqCwMIEQRBGASMc3sLuVFJ6D5Z2LHsxRvCy0P88PDB8HpGnIl7ovnTdOf7hjimJSE9EUEQhEFE4zkwxCMKDx535lYI8U5tePnuSXAsAJ/ZCrmvw/Vb4HB7uFF43TNCKZpLhznxCk02rNwbv5yQXoggCMIg43jUuNEgjsP5/kpYUgS3joLaM+DWsVBxGlxVAkENfzgEw98OTQsKX9mV+DovuPMxf7cipdUX+hERBEEYRDyfwKb/xd2OMPz1CHxsE/x+P4ww4K+zoflsuHcyTMgKl1/XBDdtg9o4E6XNzYdZOWBpuFFSWgwKJLmdIAwi6pqh7H1oDIbf9LuKXzlzIYRQOLmMfjABLoqKdj7S7sQyKKD1bMiQxHf9jiS3EwQBgMIcZ9RQcKnjZPZ+vjgsuXN4xeCOiTA6E95vdJzRQ9+E23ZAozv5zbAM+HChM8T1kg0xTyekET0SBKXUx5VSm5RStlIqriIppS5WSm1TSu1USq3oyTUFQege981yhCE61+kni8OiUXs6/L+xsDAXPj8cvlkG+06HdafABUOgPgA/q4aCN+GMD5w8SM/PdXoI/6h1egxC+tIjk5FSagZOz/TXwO1a6w72HaWUD9gOXABUAauBa7TWmzs7v5iMBKF3+N0++FyU0ziZ2djabVhZCb+qhoMe38JQHxwLwuwc2LA4pVUVuki/mYy01lu01p2NRF4M7NRa79ZatwOPApf15LqCIPSMG8Y6ApDnGaKqVsEZqxMfl2HAd8fDgTPhmWnh7cfchHobm2FjY4orK/QZfeFDGA3s86xXudsEQehnGs6Bf3minP/V5AhDfX3i4xob4dKoV8FC16F8576O5YX0oNM5lZVS/wBGxNj1Ha3135K4hoqxLa6dSim1HFgOUFZWlsTpBUHoCUvcKOcRb8Ih11lc+AFMzIBdZ3Qs39gI+TEsucuHO3MwL44RHS2kB532ELTW52utZ8f4JCMG4PQIxnrWxwD7E1zvfq31Iq31omHDkhwWIQhCjzl4FuxbGF7f3e70FjZH9Ra8YtCwCPLdVuQn+0UM0p2+MBmtBqYopSYopTKAq4Gn++C6giB0kTEFTm9hRmZ426wPYNgbzrJ36s4DiyAvzxmBBAm6/ULa0NNhp1copaqA04FnlVIvuttHKaWeA9BaB4BbgBeBLcBjWutNPau2IAi9yebToc7TW6gJRorBjkUwIs9ZnpAb3r6nqU+qJ/QSEqksCEJCLvgA/uExG2UCrUsjyxirnB5CvgH1H+q7ugkdkUhlQRB6jVeifAhtOL2FX3qynH5jlPPd0N18GcKAQARBEIS4+FaFfQN3T4TPecZ5fHlPOLX2j6f2abWEXkIEQRCEmGSsCifI+/ZouLUMfuumvwg1HKHU2l/zZDu9uLwvaymkEhEEQRA6kPsahDJTXFsCK6dE7g8uhR94woR+fjA8W9uLdX1RQ6E36DQwTRCEk4viN6DZtRNdVAgPz+5YZuxbUBU1T8JjM+CSTjOUCQMZEQRBEE4w5q1wXqLFufCCJ61FQwMUvx/uOYQY5Yfqs9wVEYS0RkxGgiAAMOMdqHZb+6lZ8O6pzvKaesdPUBAlBhe5KS9OiIGH54/0cmWFXkF6CIIgcMYa2NrqLI/0w7Yl8IPd8P9Vdiz7iwlwy7jY5/EDAeCGbXBQMs+kHSIIgnCSc+la+JebsrrIB9OyI6OSQ2w9BablJz7X6fnwRkM4SZ6QXoggCMJJzOc3wbO14fXaIKxqCK/HikpOxAszIPe9FFVO6HNEEAThJOW7O+HBOLb+eKmvOyMnp2d1EvoXcSoLwknIzyrgv6s6bv+UO79yd8Qgmrdren4OoW+RHoIgnITcVhG5/vBUuHZUas5t4EQ4X70NKktSc06hb5AegiCchHykCDKA6lOcHkGqxABgvpsOOzpwTRj4SA9BEE5Cnp7fe+d+cRYMe08mzElHpIcgCEJKKRHHctoigiAIQq+xVWZQSytEEARBSDlu4lMuWd+v1RC6iAiCIAgpZ1a2872nrX/rIXQNEQRBEFLOs3Ocb3EspxciCIIgpJwycSynJSIIgiD0KpXN/V0DIVlEEARB6BVCjuWLxLGcNoggCILQK0zKdL63tfZvPYTkEUEQBKFX+Ptc51scy+mDCIIgCL3C9Nz+roHQVUQQBEHodWrEsZwWiCAIgtBrhBzLF2zs12oISdIjQVBKfVwptUkpZSulFiUoV6GU2qCUWquUWtOTawqCkD6MzXC+10sPIS3oaQ9hI/Ax4PUkyp6rtZ6vtY4rHIIgDC6enOF82/1bDSFJejQfgtZ6C4BSqrOigiCkKQfaoMQEsxuvjwuHpL4+Qu/RVxPkaOAlpZQGfq21vj9eQaXUcmC5u9qmlBqs1scSYDDPOiv3l96k/P4G2GvjYH5+07p7YKeCoJT6BzAixq7vaK3/luR1ztRa71dKlQIvK6W2aq1jmplcsbjfvfaawWpiGsz3BnJ/6Y7cX/rSEz9tp4KgtT6/uyf3nGO/+31YKfUksJjk/A6CIAhCH9Hrw06VUrlKqfzQMnAhjjNaEARBGED0dNjpFUqpKuB04Fml1Ivu9lFKqefcYsOBN5VS64D3gGe11i8keYm4voZBwGC+N5D7S3fk/tKXbt+b0loyjQiCIAgSqSwIgiC4iCAIgiAIwAAShMGeBqML93exUmqbUmqnUmpFX9axJyilhiqlXlZK7XC/Y4Ykpdvz6+x5KIe73f3rlVIL+6Oe3SGJe1uqlKpzn9VapdR3+6Oe3UUp9aBS6nC8WKY0f3ad3Vv3np3WekB8gBk4ARWrgEUJylUAJf1d3964P8AH7AImAhnAOmBmf9c9yfu7A1jhLq8AfpTuzy+Z5wEsA57HibtaArzb3/VO4b0tBZ7p77r24B4/BCwENsbZn5bPLsl769azGzA9BK31Fq31tv6uR2+R5P0tBnZqrXdrrduBR4HLer92KeEy4A/u8h+Ay/uvKikjmedxGfCQdngHKFJKjezrinaDdP6/lhTaCX49lqBIuj67ZO6tWwwYQegCoTQY77tpLgYTo4F9nvUqd1s6MFxrfQDA/S6NUy6dnl8yzyNdn1my9T5dKbVOKfW8UmpW31Stz0jXZ5csXX52fZXLCOj7NBh9TQruL1a6lwEzLjjR/XXhNAP2+cUgmecxoJ9ZApKp9wfAOK11o1JqGfAUMKW3K9aHpOuzS4ZuPbs+FQQ9yNNgpOD+qoCxnvUxwP4enjNlJLo/pdQhpdRIrfUBt9t9OM45Buzzi0Eyz2NAP7MEdFpvrXW9Z/k5pdQvlVIlWuvBkhQuXZ9dp3T32aWVyegkSIOxGpiilJqglMoArgae7uc6JcvTwGfd5c8CHXpEafj8knkeTwPXuSNWlgB1IdPZAKfTe1NKjVDKyW2vlFqM014c7fOa9h7p+uw6pdvPrr+95R6v+BU4it0GHAJedLePAp5zlyfijIZYB2zCMcX0e91TdX/u+jJgO84IkHS6v2LgFWCH+z10MDy/WM8D+BLwJXdZAfe6+zeQYITcQPskcW+3uM9pHfAOcEZ/17mL9/cIcACw3N/e5wfRs+vs3rr17CR1hSAIggCkmclIEARB6D1EEARBEARABEEQBEFwEUEQBEEQABEEQRAEwUUEQRAEQQBEEARBEASX/x+PL48bI2cWIQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, clear_output\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "filenames = []\n",
    "\n",
    "for i in range(N):\n",
    "    print(i)\n",
    "    for j in range(numofmol):\n",
    "\n",
    "        new_y = origtraj[:i,1,j]\n",
    "        new_x = origtraj[:i,0,j]\n",
    "\n",
    "\n",
    "        #ax.plot(new_x,new_y,'b')\n",
    "        if where[j] == 1:\n",
    "            ax.plot(new_x,new_y,'deepskyblue')\n",
    "        else:\n",
    "            ax.plot(new_x,new_y,'lightskyblue')\n",
    "\n",
    "\n",
    "        plt.xlim((-1.5, 1.5))   # set the xlim to left, right\n",
    "        plt.ylim((-1.5, 1.5)) \n",
    "\n",
    "\n",
    "    plt.savefig(f'test_{i}.png')\n",
    "    filenames.append(f'test_{i}.png')\n",
    "    #ax.clear()\n",
    "    #display(fig)\n",
    "    clear_output(wait=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "import imageio.v2 as imageio\n",
    "with imageio.get_writer(f'{numofmol}mol_{theta}degstick_origin.gif', mode='I') as writer:\n",
    "    for filename in filenames:\n",
    "        image = imageio.imread(filename)\n",
    "        writer.append_data(image)\n",
    "        os.remove(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trying 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "brownian() implements one dimensional Brownian motion (i.e. the Wiener process).\n",
    "\"\"\"\n",
    "\n",
    "from math import sqrt\n",
    "from scipy.stats import norm\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def brownian_singlestep(x0, n, dt, delta, radius, theta, out=None):\n",
    "    \"\"\"\n",
    "    Generate an instance of Brownian motion (i.e. the Wiener process):\n",
    "\n",
    "        X(t) = X(0) + N(0, delta**2 * t; 0, t)\n",
    "\n",
    "    where N(a,b; t0, t1) is a normally distributed random variable with mean a and\n",
    "    variance b.  The parameters t0 and t1 make explicit the statistical\n",
    "    independence of N on different time intervals; that is, if [t0, t1) and\n",
    "    [t2, t3) are disjoint intervals, then N(a, b; t0, t1) and N(a, b; t2, t3)\n",
    "    are independent.\n",
    "    \n",
    "    Written as an iteration scheme,\n",
    "\n",
    "        X(t + dt) = X(t) + N(0, delta**2 * dt; t, t+dt)\n",
    "\n",
    "\n",
    "    If `x0` is an array (or array-like), each value in `x0` is treated as\n",
    "    an initial condition, and the value returned is a numpy array with one\n",
    "    more dimension than `x0`.\n",
    "\n",
    "    Arguments\n",
    "    ---------\n",
    "    x0 : float or numpy array (or something that can be converted to a numpy array\n",
    "         using numpy.asarray(x0)).\n",
    "        The initial condition(s) (i.e. position(s)) of the Brownian motion.\n",
    "    n : int\n",
    "        The number of steps to take.\n",
    "    dt : float\n",
    "        The time step.\n",
    "    delta : float\n",
    "        delta determines the \"speed\" of the Brownian motion.  The random variable\n",
    "        of the position at time t, X(t), has a normal distribution whose mean is\n",
    "        the position at time t=0 and whose variance is delta**2*t.\n",
    "    out : numpy array or None\n",
    "        If `out` is not None, it specifies the array in which to put the\n",
    "        result.  If `out` is None, a new numpy array is created and returned.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    A numpy array of floats with shape `x0.shape + (n,)`.\n",
    "    \n",
    "    Note that the initial value `x0` is not included in the returned array.\n",
    "    \"\"\"\n",
    "\n",
    "    x0 = np.asarray(x0)\n",
    "\n",
    "    # For each element of x0, generate a sample of n numbers from a\n",
    "    # normal distribution.\n",
    "    coordhold = np.copy(x0)\n",
    "    coordfinal = np.zeros((3,n))\n",
    "    #coordfinal = []\n",
    "\n",
    "    #c, s = np.cos(theta), np.sin(theta)\n",
    "    #rot = np.array(((c, -s), (s, c)))\n",
    "\n",
    "\n",
    "    if x0[0] < -radius:\n",
    "        print('less than -')\n",
    "        for i in range(n):\n",
    "            #print(i)\n",
    "            r = norm.rvs(size=(3,),scale=delta*sqrt(dt))\n",
    "            if (r[0] + coordhold[0]) < -radius:\n",
    "                coordhold += r\n",
    "                coordfinal[:,i] = coordhold\n",
    "            else:\n",
    "                while (r[0] + coordhold[0]) > -radius:\n",
    "                    r = norm.rvs(size=(3,),scale=delta*sqrt(dt))\n",
    "                coordhold += r\n",
    "                coordfinal[:,i] = coordhold\n",
    "    elif x0[0] > radius:\n",
    "        print('greater than +')\n",
    "        for i in range(n):\n",
    "            #print(i)\n",
    "            r = norm.rvs(size=(3,),scale=delta*sqrt(dt))\n",
    "            if (r[0] + coordhold[0]) > radius:\n",
    "                coordhold += r\n",
    "                coordfinal[:,i] = coordhold\n",
    "            else:\n",
    "                while (r[0] + coordhold[0]) < radius:\n",
    "                    r = norm.rvs(size=(3,),scale=delta*sqrt(dt))\n",
    "                coordhold += r\n",
    "                coordfinal[:,i] = coordhold\n",
    "    else:\n",
    "        print('inside')\n",
    "        for i in range(n):\n",
    "            #print(i)\n",
    "            r = norm.rvs(size=(3,),scale=delta*sqrt(dt))\n",
    "            if (r[0] + coordhold[0]) > -radius and (r[0] + coordhold[0]) < radius and (r[1] + coordhold[1]) > -radius and (r[1] + coordhold[1]) < radius:\n",
    "                coordhold += r\n",
    "                coordfinal[:,i] = coordhold\n",
    "            else:\n",
    "                while (r[0] + coordhold[0]) < -radius or (r[0] + coordhold[0]) > radius:\n",
    "                    r = norm.rvs(size=(3,),scale=delta*sqrt(dt))\n",
    "                coordhold += r\n",
    "                coordfinal[:,i] = coordhold\n",
    "\n",
    "    coordfinal = np.squeeze(np.array(coordfinal))\n",
    "\n",
    "    # Add the initial condition.\n",
    "    #coordfinal += np.expand_dims(x0, axis=-1)\n",
    "\n",
    "    #coordfinal = rot @ coordfinal\n",
    "\n",
    "    return coordfinal"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "e554f0a28730f968725323226da918ce8c300a334dc9c93b2dac8e99d3c39b00"
  },
  "kernelspec": {
   "display_name": "Python 3.9.5 ('venv')",
   "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.9.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}